基于FFT－MA谱模拟的快速随机反演方法研究

王保丽，印兴耀，丁龙翔，2，张广智，孙瑞莹

1 中国石油大学（华东）地球科学与技术学院，青岛 266580

2 广州海洋地质调查局，广州 510075

1 引言

地震反演的实质是将地震观测数据转化为储层参数，用反演的储层参数信息进行地下目的储层的预测，从而为油田的勘探开发提供技术支持（Bosch，2010；Doyen，2007）.其中，确定性反演和随机反演是地震反演的两大类别（杨晓春等，2001；Sams和Saussus，2010）.与确定性反演方法相比，以地质统计学为基础、以测井资料为条件数据的随机反演方法能够模拟地震频带以外的信息，具有比常规确定性反演更高的分辨率（Francis，2005；Moyen和Doyen，2009；Helgesen，2000；Sancevero，2005；Sams和Saussus，2008），因此，受到了国内外地球物理学家的广泛关注与研究.

在国外，Bortoli等（1992）首先提出了基于序贯高斯模拟和单道匹配最优原则的随机反演方法，该方法能够在横、纵向变差以及井位置处测井资料的约束下，有效地融入地震频带以外的测井信息，提高反演解的分辨率.Haas和Dubrule（1994）发展了单道模拟的随机反演方法，捕获了超出地震带宽的频率信息.单道模拟方法虽然可以获得高分辨率的反演结果，但是计算效率低，限制了其在实际资料中的应用.Debeye（1996）、Sams（1999）和 Contreras（2005）等通过引入模拟退火和马尔科夫链蒙特卡洛（MCMC）技术，用点扰动的方式取代了单道反复模拟，扩展了随机反演算法.Srivastava和Sen（2010）提出了基于分形初始模型的随机反演方法，该方法能够反演出纵横波阻抗等多个参数数据体.以上方法仍然没有解决计算耗时的问题.在国内，近年来，地球物理学家也在随机反演方面做了很多研究.黄捍东等（2009）研究了以测井、地质信息作为约束条件的非线性随机反演方法.张世鑫等（2010）针对目的储层薄互层发育的特点，运用随机反演方法成功实现了储层的精细预测.Huang等（2012）研究指出随机反演在薄储层地震反演中具有更大的优势，并对随机反演中关键参数的优选及其效果进行了分析.杨锴等（2012）将测井、井间地震与地面地震数据同时用于介质参数的随机反演，大大降低了模拟出的介质参数的不确定性，提高了反演精度.

尽管高分辨率的随机反演工作在实际资料应用中已取得了一定的进展，但仍普遍存在计算耗时的问题.鉴于此，本文针对常规随机反演计算效率低，占用内存大的问题，用傅里叶滑动平均（Fast Fourier Transform－Moving Average，FFT－MA）模拟方法进行频率域的快速模拟（Pardo－Igúzquiza，1993；Ravalec和Noetinger，2000；Oliver，1995），并利用逐步变形算法（Gradual Deformation Method，GDM）（Hu，2000；Hu和Ravalec，2004；Ravalec和 Noetinger，2002）确保反演解的快速收敛，构建了基于FFT－MA谱模拟的快速随机反演方法.模型试算与实例分析表明，与常规随机反演相比，该方法能够在保证高分辨反演结果的情况下，确保反演解的快速收敛，大大提高随机反演的计算效率.

2 方法原理

2.1 FFT－MA非条件谱模拟方法

FFT－MA方法可进行快速的地质统计模拟，计算效率高，且不受网格大小限制，是一种非条件模拟方法.其基本公式为

式中，y表示模拟值，m表示模拟的均值，g表示协方差矩阵C的共轭根，是g的共轭，＊为卷积符号，z表示满足模拟网格大小的随机高斯白噪.

该方法实现的关键之一是协方差矩阵C的构建，其中一维协方差矩阵C（h）和变差函数之间满足C（h）＝σ2－γ（h），σ2为模拟的方差，γ（h）为变差函数，h为空间距离.

2.2 条件化处理

由于FFT－MA谱模拟是一种非条件模拟，不能满足硬数据（已知井数据）.为解决该问题，本文采用克里金条件化方法（Francis，2006a，2006b）对其进行了条件化处理：

图1是FFT－MA谱模拟条件化处理的示意图，首先求取硬数据点处数值（蓝色圆圈）的克里金得到待条件化点处的克里金值ydk（x）（红色圆圈），然后再求取硬数据点处的非条件模拟值（橙色三角）的克里金得到待条件化点处的克里金值yusk（x）（绿色三角），yusk（x）与ydk（x）的差值再加上待条件化点处的非条件模拟值yus（x），就可以得到待条件化点处的条件化模拟值ycs（x）.

2.3 GDM算法寻优

从式（1）中可以看出，FFT－MA谱模拟的关键问题是高斯白噪z的更新，传统的随机反演通常采用反复模拟的方式来寻找最优结果，不但计算效率

其中，zpropose（t）为更新后的白噪，zcurrent为当前待更新的白噪，znew为加入的新的白噪.当t＝0时，zpropose＝zcurrent，当t从0变化为π／2时，zpropose＝znew，即逐步替换为新的白噪，设置合适的变形参数t，便能够搜索zcurrent和znew之间的参数空间，每次加入一个znew便构成一个新的搜索链，通过多个链的搜索最终可以找到目标函数的最小值.

需要指出的是，高斯白噪z的存在使FFT－MA算法具有随机性，但由于协方差结构C不改变，所以这种扰动不会影响数据的空间变差结构.低，而且不能保证每次的模拟结果都比上一次更接近实际地震记录.因此，需要寻找一种合适的优化算法来改善反演解的收敛问题.

本文在FFT－MA算法扰动模拟的基础上，通过引入GDM不断更新模拟结果，使每次的模拟结果都比上一次结果更逼近实际地震数据，从而加快反演的求解速度.

GDM算法的核心思想是将一个高斯白噪的局部或整体加入另一个高斯白噪的对应部分，其结果仍是一个高斯白噪.具体可以表示为：

3 算法测试与分析

3.1 FFT－MA谱模拟与序贯高斯模拟的对比

为验证FFT－MA谱模拟的计算效率，本文基于相同的硬数据和变差，分别进行了条件化FFT－MA谱模拟和序贯高斯模拟，其结果如图2所示.条件化FFT－MA谱模拟比序贯高斯模拟的分辨率高.图3是二者的计算效率对比，可以看出，FFT－MA谱模拟的时间从154.84s降至11.20s，速度提高了14倍.因此，与序贯高斯模拟相比，频率域的FFT－MA谱模拟方法不但分辨率高，而且能够有效解决计算耗时、耗内存的问题.

3.2 FFT－MA扰动模拟特性分析

由方程（1）可知，FFT－MA谱模拟的最终结果受协方差C和高斯白噪z的影响，其中模拟的随机性来自于加入的随机高斯白噪z.图4是一个50点的模拟，采用球状变差模型，变程值为10，在更新第25点的z值时，模拟结果在相关长度10以内均发生了扰动.图5和图6分别是网格大小均为80×60时的点和线扰动的情况，图5是对位置（30，5）的点扰动，图6是对第11列的线扰动，从图5和图6最右边的两者残差图可以看出，在相关半径内，模拟结果均发生了扰动.

模拟分析结果表明，局部（一个点或者一条线）高斯白噪z值的扰动，会引起FFT－MA谱模拟结果的局部扰动.正由于FFT－MA谱模拟的这种特性，在对实际地震数据进行模拟时，通过对z值在某道位置的扰动，可以产生相应道位置处的模拟波阻抗的扰动，因此，可以对实际地震数据按道进行模拟，实现模拟模型的连续修改，最终达到与实际地震记录的匹配.

图1 条件化处理示意图Fig.1 Sketch map of conditional process

图2 条件化FFT－MA谱模拟与序贯高斯模拟的结果对比图Fig.2 Simulation results comparison of conditional FFT－MA spectrum simulation with sequential Gauss simulation

图3 条件化FFT－MA谱模拟与序贯高斯模拟的计算效率对比图Fig.3 Computational efficiency comparison of conditional FFTMA spectrum simulation with sequential Gauss simulation

图4 点（25）扰动前后模拟结果对比图Fig.4 Comparison of simulation results before and after point（25）disturbance

图5 点（30，5）扰动前后模拟结果对比图Fig.5 Comparison of simulation results before and after point（30，5）disturbance

图6 线（第11列）扰动前后模拟结果对比图Fig.6 Comparison of simulation results before and after line（Line 11）disturbance

4 基于FFT－MA谱模拟的快速随机反演

根据以上所述的方法原理，本文建立了基于FFT－MA谱模拟的快速随机反演方法，其基本步骤为：在进行地质统计模拟时，首先通过地质统计方法计算测井曲线的变差函数和协方差，以确定测井数据的空间结构特征，然后利用方程（1）和（2）重构出满足指定空间结构的条件化的谱模拟结果，最后采用逐步变形算法不断更新模拟结果，确保模拟结果与实际地震数据的匹配.图7表示其具体反演流程.

4.1 目标函数的建立

为使随机模拟的地震波阻抗与实际地震数据达到最佳匹配，本文将目标函数定义为：

其中，下标i为各个地震、模型参数的采样点，w为加权系数，g为正演算子，msim为模拟波阻抗，dobs为实际地震记录，Nobs为采样点数，m0是先验估计的阻抗，v是加权系数，α是正则化系数.

图7 基于FFT－MA谱模拟的快速随机反演流程图Fig.7 Flow chart of the Fast Stochastic Inversion based on FFT－MA spectrum simulation

一般地震波阻抗反演的目标函数是N维的，本文通过引入逐步变形算法，将反问题的求解转变为对最佳参数t的搜索，从而变成了一维问题，简化了问题的求解，提高了运算效率.因此，即使采用迭代的思想也可达到目标函数的快速求解.

4.2 目标函数的求解

假设一次新的znew的产生即一个新链，每个链上有m个t值，如果在状态tn下得到的目标函数J的数值相对状态tn－1下降了，即J（tn）＜J（tn－1），则跳出最优t的搜索，然后将下一个链的zcurrent更新为当前状态tn下所得的zpropose（公式3），如果目标函数不降低，则保持该链的zcurrent不变，带入下一个链的计算，如此反复迭代直至目标函数满足收敛条件.

具体步骤可描述为：

①：生成新的znew，与上一步产生的zcurrent组合产生一条新的链；

②：搜索链上的t值，如果存在tn，停止搜索t，将下个链的zcurrent更新为该链上tn时的zpropose，转到①；

③：搜索该链上的t值，如果目标函数一直不下降，zcurrent保持不变，转到①；

④：重复①—③直到目标函数满足收敛条件.

5 模型试算及结果分析

5.1 模型试算

为验证算法的有效性，本文首先采用8000个链（即8000次迭代）对单道数据进行了数值试算，t的取值范围为0到π／2，间隔为π／20，所得单道反演结果如图8所示，从图中可看出，模型波阻抗与反演波阻抗吻合比较好.图9是其合成地震记录和实际地震记录的对比图，可以看出，模型波阻抗与反演波阻抗吻合非常好.图10是归一化目标函数的单道收敛效果，可以看出，达到迭代次数时，归一化后的目标函数值为0.0379，即下降了96.3%，在实际计算时，迭代次数为1000多次时即能够使目标函数达到很好的收敛，从而使运算时间大大缩短.

本文还对图11所示的二维模型进行了测试分析，该模型在第250个采样点位置处有两个薄层.具体反演时，从该模型中抽取5道作为伪井，将模型波阻抗的合成地震记录作为实际地震记录.图12为反演结果，其中图12a为稀疏脉冲反演方法所得剖面，图12b为本文方法所得剖面.与确定性稀疏脉冲方法相比，本文方法能够得到与初始波阻抗（图11）更加相似的结果，薄层的可辨识性得到提高，尤其在黑色椭圆标出区域，基于FFT－MA谱模拟的快速随机反演可以明显识别出两个薄层，而稀疏脉冲反演只能识别出一个层.因此，利用本文提出的基于FFTMA谱模拟的快速随机反演方法，能够有效融合测井资料的高频信息，提高地震反演分辨率，识别薄层，为复杂薄互储层的识别提供技术支持.从模型合成地震记录与由反演结果所得合成记录的相关系数对比图（图13）可以看出，大部分相关系数都大于0.95，体现出二者良好的相关性.

5.2 实例分析

实际资料来自国内某老油田，该工区目的层段为河流相沉积，河道发育，储集层砂岩纵横向变化大、连通性差、油水关系复杂.工区目的层段砂体发育，纵向上砂岩与泥岩呈互层分布.

本次试算使用的地震资料共有157道，如图14所示，纵向上采样率2ms，从1.1s到1.3s，提取的子波主频为30Hz，采用了三口井，分别位于10、88、142地震道位置处.

图15显示了实际资料的反演结果，其中图15a为稀疏脉冲反演方法所得结果，图15b为本文方法所得结果，通过对比可以看出，虽然图15a所示反演结果基本满足反演要求，但纵向分辨率低，对薄互储层展示效果不好，不能真实反映地下地质情况，因此难以满足油田生产需求.图15b所示反演结果的纵向分辨率明显提高，砂泥岩互层得到很好的区分，砂体横向连续性也比较合理.从图中黑色椭圆标出的位置可以看出，在1.25s目的层附近的两个薄层得到了很好的区分，1.2s附近的薄单砂体也能够清晰地分辨，验证了本文提出的基于谱模拟的快速随机反演方法在薄储层识别中的优势以及本文研究方法的有效性.

图8 单道波阻抗反演结果与模型对比图Fig.8 Comparison of impedance inversion result with model data

图9 反演结果合成记录与实际地震记录的对比图Fig.9 Comparison of Synthetic seismic data of inverted impedance and model data

图10 归一化目标函数的收敛效果图Fig.10 Convergence effect of normalized objective function

6 结论

图11 波阻抗模型Fig.11 Impedance model

地震随机反演是一种可有效提高反演分辨率的方法，但存在计算效率低的问题，本文构建的基于FFT－MA谱模拟快速随机反演方法，不仅能够有效融合测井资料的高频信息，获得高分辨率反演结果，而且能够借助FFT－MA在频率域的快速模拟和GDM算法的快速优化更新，一方面加快随机反演的速度，有效提高反演的运算效率，另一方面使反演解得到快速收敛.模型试算和实例分析表明，快速随机反演方法不但可以提高计算效率，而且能够得到与模型吻合的反演结果，比稀疏脉冲反演方法更有效的识别出薄层.另外，由于新方法不受序贯因素的影响，因此，可以在今后的研究中对该方法进行并行运算，将储层分为多个部分同时进行反演，进一步提高计算效率.

图12 波阻抗反演结果 （a）稀疏脉冲反演，（b）基于FFT－MA谱模拟的快速随机反演方法Fig.12 Impedance inversion results of （a）Sparse Impulse Inversion and （b）Fast Stochastic Inversion based on FFT－MA spectrum simulation

图13 模型合成记录与反演合成记录的相关系数对比图Fig.13 Correlation coefficient comparison of model synthetic with inversion synthetic

图14 实际地震资料Fig.14 Real Seismic data

图15 波阻抗反演结果 （a）稀疏脉冲反演；（b）基于FFT－MA谱模拟的快速随机反演方法Fig.15 Impedance inversion results of（a）Sparse Impulse Inversion and（b）Fast Stochastic Inversion based on FFT－MA spectrum simulation

Bortoli L J，Alabert F，Haas A，et al.1992.Constraining stochastic images to seismic data.／／Soares A ed.Geostatistics Tróia′92.Netherlands：Springer，325－337.

Bosch M，Mukerji T，Gonzalez F E.2010.Seismic inversion for reservoir properties combining statistical rock physics and geostatistics：A review.Geophysics，75（5）：75A165－75A176.Contreras A，Torres－Verdin C，Kvien K，et al.2005.AVA stochastic inversion of pre－stack seismic data and well logs for 3D reservoir modeling.／／EAGE 67th Conference and Exhibition.

Debeye H W，Sabbah E，Made P M.1996.Stochastic inversion.／／66Ann.Internat Mtg.，Soc.Expi.Geophys.Expanded Abstracts，1212－1215.

Doyen P M.2007.Seismic Reservoir Characterization：An Earth Modeling Perspective.Houten，The Netherlands：European Association of Geoscientists and Engineers，Publications.

Francis A.2005.Limitations of deterministic and advantages of stochastic seismic inversion.CanadianSoc.ExpiGeophys.Recorder，30（2）：5－11.

Francis A M.2006a.Understanding stochastic inversion：Part 1.FirstBreak，24（11）：69－77.

Francis A M.2006b.Understanding stochastic inversion：Part 2.FirstBreak，24（12）：79－84.

Haas A，Dubrule O.1994.Geostatistical inversion—a sequential method of stochastic reservoir modeling constrained by seismic data.FirstBreak，12（11）：561－569.

Helgesen J，Magnus I，Prosser S，et al.2000.Comparison of constrained sparse spike and stochastic inversion for porosity prediction at Kristin field.TheLeadingEdge，19（4）：400－407.

Hu L Y.2000.Gradual deformation and iterative calibration of Gaussianrelated stochastic models.Math.Geology，32（1）：87－108.

Hu L Y，Ravalec－Dupin M L.2004. An improved gradual deformation method for reconciling random and gradient searches in stochastic optimizations.Math.Geol.，36（6）：703－719.

Huang H D，Zhang R W，Wei S P.2009.Research on application of seismic nonlinear random inversion to reservoir prediction in the thin sandstone of continental deposits.ActaPetroleiSinica（in Chinese），30（3）：386－390.

Huang Z Y，Gan L D，Dai X F，et al.2012.Key parameter optimization and analysis of stochastic seismic inversion.AppliedGeophysics，9（1）：49－56.

Moyen R，Doyen P M.2009.Reservoir connectivity uncertainty from stochastic seismic inversion.／／79th Ann.Internat Mtg.，Soc.Expi.Geophys..Expanded Abstracts，2378－2381.

Oliver D S.1995.Moving averages for Gaussian simulation in two and three dimensions.Math.Geol.，27（8）：939－960.

Pardo－Igúzquiza E，Chica－Olmo M.1993.The Fourier integral method：An efficient spectral method for simulation of random fields.Math.Geol.，25（2）：177－217.

Ravalec M L，Noetinger B，Hu L Y.2000.The FFT moving average（FFT－MA）generator：An efficient numerical method for generating and conditioning Gaussian simulations.Math.Geol.，32（6）：701－723.

Ravalec－Dupin M L，Noetinger B.2002.Optimization with the gradual deformation method.Math.Geol.，34（2）：125－142.

Sams M，Saussus D.2008.Comparison of uncertainty estimates from deterministic and geostatistical inversion.／／78th Ann.Internat Mtg.，Soc.Expi.Geophys..Expanded Abstracts，1486－1490.

Sams M S，Atkins D，Siad N，et al.1999.Stochastic inversion for high resolution reservoir characterisation in the central Sumatra basin.／／Society of Petroleum Engineers Asia Pacific Improved Oil Recovery Conference，257－260.

Sams M S，Saussus D.2010.Comparison of lithology and net pay uncertainty between deterministic and geostatistical inversion workflows.FirstBreak，28（2）：35－44.

Sancevero S S，Remacre A Z，Portugal R S，et al.2005.Comparing deterministic and stochastic seismic inversion for thin－bed reservoir characterization in a turbidite synthetic reference model of Campos basin，Brazil.TheLeadingEdge，24（11）：1168－1172.

Srivastava R P，Sen M K.2010.Stochastic inversion of prestack seismic data using fractal－based initial models.Geophysics，75（3）：R47－R59.

Yang K，Ai D F，Geng J H.2012.A new geostatistical inversion and reservoir modeling technique constrained by well－log，crosshole and surface seismic data.ChineseJ.Geophys.（in Chinese），55（8）：2695－2704，doi：10.6038／j.issn.0001－5733.2012.08.022.

Yang X C，Li X F，Zhang M G.2001.Process and Mathematical basis of investigations of seismic inversion.Process in Geophysics（in Chinese），16（4）：96－109，doi：10.3969／j.issn.1004－2903.2001.04.013.

Zhang S X，Zhang F C，Li S P，et al.2010.Detailed prediction of thin reservoir in Andian area of Biyang Depression.Geophysical Prospecting for Petroleum （in Chinese），49（1）：34－41.

附中文参考文献

黄捍东，张如伟，魏世平.2009.地震非线性随机反演方法在陆相薄砂岩储层预测中的应用.石油学报，30（3）：386－390.

杨锴，艾迪飞，耿建华.2012.测井、井间地震与地面地震数据联合约束下的地质统计学随机建模方法研究.地球物理学报，55（8）：2695－2704，doi：10.6038／j.issn.0001－5733.2012.08.022.

杨晓春，李小凡，张美根.2001.地震波反演方法研究的某些进展及其数学基础.地球物理学进展，16（4）：96－109，doi：10.3969／j.issn.1004－2903.2001.04.013.

张世鑫，张繁昌，李少鹏等.2010.泌阳凹陷安店区块薄储层精细预测.石油物探，49（1）：34－41.