地球物理学报  2011, Vol. 54 Issue (11): 2926-2932   PDF    
引用本文
张广智, 王丹阳, 印兴耀, 等. 基于MCMC的叠前地震反演方法研究. 地球物理学报, 2011, 54(11): 2926-2932, doi: 10.3969/j.issn.0001-5733.2011.11.022.  
ZHANG G Z, WANG D Y, YIN X Y, et al. Study on prestack seismic inversion using Markov Chain Monte Carlo. Chinese J. Geophys. (in Chinese), 2011, 54(11): 2926-2932, doi: 10.3969/j.issn.0001-5733.2011.11.022.  
基于MCMC的叠前地震反演方法研究
张广智, 王丹阳, 印兴耀, 李宁     
中国石油大学(华东) 地球科学与技术学院,青岛 257061
收稿日期 2010-04-06, 2010-06-28 收修定稿
基金项目: 国家油气重大专项(2011ZX05014-001-010HZ)、中国石油科技创新项目(2011D-5006-0301)、中国石油大学(华东)自主创新科研计划项目(11CX05006A)联合资助
作者简介: 张广智,男,1971年生,博士,中国石油大学(华东)副教授,主要研究方向为地震储层预测与描述方法研究.E-mail:zhanggz@upc.edu.cn
摘要: 马尔科夫链蒙特卡洛方法(MCMC)是一种启发式的全局寻优算法[1].它在贝叶斯框架下,利用已有资料进行约束,既可使最优解满足参数的统计特性,又通过融入的先验信息,提高解的精度;寻优过程可跳出局部最优,得到全局最优解.利用MCMC方法,可以得到大量来自于后验概率分布的样本,不仅可以得到每个未知参数的估计值,而且可以得到与之相关的各种不确定性信息.此外,由于算法并不是利用有单一最优解的目标函数,所以结果对初始值的依赖不强.通过对简单一维层状介质模型的处理,和实际资料的应用,说明利用基于Metropolis-Hastings算法的MCMC方法进行地震反演,通过对解空间的随机搜索能够得到较好的效果.
关键词: 非线性反演      MCMC      叠前地震反演      Metropolis-Hastings算法     
Study on prestack seismic inversion using Markov Chain Monte Carlo
ZHANG Guang-Zhi, WANG Dan-Yang, YIN Xing-Yao, LI Ning     
Shool of Geosciences, China University of Petroleum, Qingdao 257061,China
Abstract: Markov Chain Monte Carlo method is a heuristic global optimization method. In the Bayesian framework, the optimal solution meets the statistical properties of the parameters with the constraints of data(eg. seismic data, well log); the accuracy of solution is improved with the prior information joined in; optimization process can jump out of local optimum, and ultimately obtain the global optimal solution. Using MCMC methods, we draw a large number of samples from the posterior distribution function. With these samples, we obtain not only the estimates of each unknown variable, but also various types of uncertainty information associated with the estimation. In addition, because of not making use of a objective functions with single optimal solution, the results obtained with MCMC method are independent of the choice of initial values. By testing a 1D layered model and the application of real seismic data in South China, shows the MCMC method based on Metropolis-Hastings algorithm is available and can obtain good results by random searching solution space.
Key words: non-linear inversion      Markov Chain Monte Carlo      prestack seismic inversion      Metropolis-Hastings algorithm     
1 引言

叠后地震反演使用全角度多次叠加的地震数据，只能提供种类很少的纵波波阻抗等参数，不能给出反映物性、流体特征的参数，在研究储层物性、流体方面受到了限制，因此，常规的叠后反演已经不能满足日益增长的精细储层描述的要求[2].与叠后反演相比，叠前反演[3~6]具有良好的保真性，可以充分利用叠前地震数据丰富的振幅和旅行时信息，提供研究岩性和储层变化规律的更多、更有效的数据体，适合进行含油气性反演[7~12].笔者根据Fatti近似公式，提出一种基于MCMC 的叠前反演方法.

1 MCMC 方法简介

MCMC方法属于蒙特卡洛方法.Hammersley and Handscomb(1964)将蒙特卡洛方法定义为“有关随机数的实验的实验数学的分支"，其基础是对概率分布的直接模拟[13].

目前MCMC 方法已经广泛应用于理论物理、社会学、生物、医学、人工智能、金融统计、信号通信等领域.在过去的十五年间，地球物理学家逐渐开始使用MCMC 方法进行反演[14~22]，直接模拟后验概率分布函数(PDF)，也就是抽取符合后验PDF 的随机样本，并用这些样本计算解的贝叶斯估计.该方法也可适用于解决具有复杂的多峰多维PDF的问题.

Alberto Malinverno(2002)利用直流电阻率测深数据反演一维地球模型，这是首次将MCMC 方法应用于完全非线性反演问题[23].

1.1 基本原理

马尔科夫链蒙特卡洛(MCMC)方法用于对贝叶斯推理中的后验概率分布进行抽样，它通过抽取收敛于贝叶斯后验分布的随机样本，再对这些样本进行统计，来间接得到后验分布的一些性质[24].

马尔科夫链是数学中具有马尔科夫性质的离散时间过程.若{Xt:t≥0}为一随机序列.随机序列所有可能取到的值组成的集合记为S，称为状态空间.如果对于$\forall $t≥0及任意状态si，sj，si0，…，sit-1都有

(1)

则称{Xt:t≥0}为一马尔科夫链[25].直观上看，对于马尔科夫过程，要预测将来的唯一有用的信息就是过程当前的状态，而与以前的状态无关.

用πj(t)=P(Xt=sj)表示马氏链在t时刻处于状态sj的概率.则对于在时刻t+1马尔科夫链有状态si的概率πi(t+1)可以由Chapman-Kolomogorov方程得到:

(2)

它可写为更紧凑的矩阵形式:π(t+1)=π(t)P.连续应用可得到:π(t)=π(0)Pt.Chapman-Kolomogorov方程的这种连续迭代描述了马尔科夫链的发展[26].

如果所构造的转移核P(s，s′)满足π(s′)P(s′，s)=π(s)P(s，s′)，那么由其构造的马尔科夫链将具有唯一的不变分布.这个条件称为细致平衡条件.对于一条存在唯一不变分布的马尔科夫链，当它经过充分长的迭代后，马尔科夫链将收敛于此，π(s)就可称为平稳分布.对于本文的反演问题，我们希望得到所估参数的后验概率分布，所以各个参数的马尔科夫链应收敛至所估计参数的后验概率分布.

1.2 Metropolis-Hastings算法

常用的构造转移核的方法有Metropolis-Hastings算法和Gibbs采样算法.本文选用Metropolis-Hastings算法，其构造马尔科夫链的方法如下[27]:

若要使π(x)为平稳分布.首先由建议分布q(·|x)产生一个潜在的转移x￫x*，然后根据概率α(x，x*)来决定是否转移，也就是说在潜在转移点x* 找到后，以概率α(x，x*)接受x* 作为链在下一时刻的状态值.于是，在有了x* 后，可从[0， 1]的均匀分布上，抽取一个随机数u，则

常用的选择是:

易证接受概率的此种形式，可使构造出的转移核满足细致平衡条件.

2 MCMC 叠前反演方法原理

在本文提出的方法中，正演部分使用的是Zoeppritz方程的Fatti波阻抗近似式[28， 29].

(3)

其中珔γ是横纵波速度的比值的均值，rp、rs、rd 分别是纵、横波阻抗反射系数，密度反射系数，并且有

(4)

IP =α·ρ 是纵波阻抗，Is =β·ρ 是横波阻抗，珋θ是界面两侧入射角的均值.与其它近似公式相比，该公式适用于高角度的入射角(小于50°).

本文中的反问题可以表述如下:

(5)

其中d为地震记录.未知参数rp, rs, rd 的后验概率密度为:

(6)

其中第一项为先验概率密度，假设rp, rs, rd 相互独立，且服从高斯分布，其均值分别为μrp, μrs, μrd, 标准差分别为σrp, σrs, σrd, 则有

(7)

第二项为似然函数，假设噪音e为均值为0，方差为σe的高斯随机噪音，可得到观测地震记录的似然函数为:

(8)

为了得到所估参数的后验概率分布p(rp, rs, rd|d)，利用Metropolis-Hasting 算法来生成马尔科夫链.其建议分布函数取为均匀分布:

该建议分布函数满足对称的随机游走，其中r代表rp, rs, rd.

所以转移概率为:

(9)

为了使马尔科夫链最终收敛于未知参数r的后验概率分布，平稳分布应为p(rp, rs, rd)p(d|rp, rs, rd).令

则转移概率α(ri，r*)= exp(α′(ri，r*))，其中α′(ri，r*)= min 0，f(r*)-f(ri{)}.

3 数值分析 3.1 模型试算

随机生成51个横波速度β 的值，利用经验公式得到纵波速度α 和密度ρ 的值[30， 31].

(10)

(11)

入射角为0~30°，间隔为3°，计算得到随角度变化的PP波反射系数R(珋θ)，与采样间隔为4 ms, 50Hz的雷克子波褶积，合成精确的地震记录.加上一定量的噪音，得到用于反演的观测地震记录，其信噪比为6.48(代表精确的地震记录，dob-代表用于反演的观测地震记录).

250000次迭代后的估计值，如图 1所示，其中3幅图从左至右分别是零偏移距纵波阻抗、横波阻抗和密度曲线rp, rs, rd, 红色虚线为反演值，黑色实线线为实际值.从图中可以看出，反演得到的曲线与实际曲线重合较好，反演效果较好.

图 1 用合成数据进行反演的参数(红色虚线）与实际参数(黑色实线）的对比. (a)纵波速度反射系数，（b)横波速度反射系数，（c)密度反射系数. Fig. 1 The comparison between estimations(red dashed line)，and true values (black solid line). (a) P-wave reflection coefficient； (b)S-wave reflection coefficient; (c) the density reflection coefficient
图 2 地震记录数据的对比 (a)用于反演的地震记录数据，(b)用反演结果合成的地震记录. Fig. 2 The comparison of seismic data (a) The seismic data used for nversion;(b) The synthetical seismic data using the estimation
3.2 实际叠前地震资料反演

实际资料共有200道，井位于第100道，初始模型为rp, rs, rd 的低频趋势，用上述方法进行反演.反演测线上的前10个CMP道集如图 3所示.

图 3 实际资料的前10个叠前CMP道集 Fig. 3 Real prestack dataCFirst ten)

在反演过程中，统计井的信息，作为参数的先验信息.图 3由左至右分别是由井曲线统计出来的rp, rs, rd 的直方图.图 4 是第100 道的反演值的直方图，也就是后验概率分布.图 3中拟合出的高斯分布的标准差分别为0.0162，0.0167，0.0067.图 4 为0.0051，0.0064，0.0022.显然后验分布比先验分布的范围窄，这是由于似然函数的参与，所以增加了约束改变了先验分布.

图 4 井曲线的直方图 (a) rp , (b) rs, (c) rd. Fig. 4 The histogram of rp ,rs ,rd according to the logging trace statistics

图 6(a)是第100道，2190 ms处纵波阻抗反射系数的马尔科夫链，(b)是2176ms处横波阻抗反射系数的马尔科夫链，(c)是2176ms处密度反射系数的马尔科夫链，可以看到随着迭代的进行，马氏链逐渐收敛，直至平稳.对这些马尔科夫链进行统计，就可以得到图 5的后验概率分布.

图 5 反演值的直方图 (a) rp , (b) rs, (c) rd. Fig. 5 The histogram of estimations in CDP 100
图 6 第100道某点的马尔科夫链 (a) rp的马氏链，（b)rs的马氏链，（c）rd的马氏链 Fig. 6 The Markov Chains in CDP 100

反演出的反射系数剖面如图 7，依次为纵波阻抗反射系数剖面，横波阻抗反射系数剖面，密度反射系数剖面.在CDP100处插入的是相应的测井曲线，反演结果与井曲线相符合.图 7的红色圆圈，表明了气层所在，在纵波阻抗反射系数剖面上含气储层的位置比较清楚.图 8 为实际的叠后地震剖面与合成地震剖面，可见二者几乎一致.

图 7 实际资料叠前反演结果 (a)纵波阻抗反射系数剖面；（b)横波阻抗反射系数剖面；（c)密度反射系数剖面. Fig. 7 Real data prestack inversion results
图 8 地震剖面的对比 (a)实际的地震剖面；（b)由反演结果合成的地震剖面 Fig. 8 The comparison of seismic profiles
4 结论

本文在AVO 的理论基础上，使用了Zoeppritz方程的近似式.反演时所采用的是叠前动校后的道集，它可以反映地震波振幅随角度的变化情况，根据振幅的变化可以反演出纵、横波速度和密度.纵波速度和密度是随横波速度变化的.为了使得到的参数(如，泊松比等)在合适的范围内变化，所以在建立模型时采用了经验公式.

本文使用MCMC 方法进行叠前地震资料的非线性反演，通过模型试算，理论模型与反演结果基本吻合，以及对实际资料的应用，说明了算法的有效性和可靠性.限于篇幅限制，只从理论上做了简单的分析和应用.

由于MCMC方法是在贝叶斯框架下进行非线性随机反演，它可以结合多种信息，比常规反演方法精度更高，同时还可分析反演结果的不确定性，进而进行风险评估.本文未对此进行论述，还须后续研究.

此外，在利用MCMC 方法进行反演时，不同的建议分布、约束条件及迭代次数等，均会对反演精度及算法效率产生影响.先验分布给出了参数的统计特性；建议分布控制着候选点的产生，也就是随机游走的路径，会对后验概率分布的形状产生影响，也会影响算法的效率；迭代次数与建议分布相配合，若若建议分布范围窄，需要很多次迭代；若建议分布范围宽，则无需大的迭代次数，但同时又无法保证反演精度.所以，在选取这些参数时，需要综合考虑，此方面还需要深入研究和分析.

致谢

感谢中国石油化工股份有限公司石油勘探开发研究院为本项目的研究提供实际地震资料.

参考文献
[1] 王家映. 地球物理资料非线性反演方法讲座(二)蒙特卡洛法. 工程地球物理学报 , 2007, 4(2): 81–85. Wang J Y. Monte Carlo Method. Chinese Journal of Engineering Geophysics (in Chinese) (in Chinese) , 2007, 4(2): 81-85.
[2] 杨培杰, 印兴耀. 基于支持向量机的叠前地震反演方法. 中国石油大学学报 (自然科学版) , 2008, 32(1): 37–41. Yang P J, Yin X Y. Prestack seismic inversion method based on support vector machine. Journal of the University of China University of Petroleum (Edition of Natural Science) (in Chinese) (in Chinese) , 2008, 32(1): 37-41.
[3] 殷八斤, 曾灏, 杨在岩. AVO技术的理论与实践. 北京: 石油工业出版社, 1995 .
[4] 陈天胜, 刘洋, 魏修成. 纵波和转换波联合AVO反演方法研究. 中国石油大学学报 (自然科学版) , 2006, 30(1): 33–37. Chen T S, Liu Y, Wei X C. Joint amplitude versus offset inversion of PP and PSV seismic data. Journal of China University of Petroleum (Edition of Natural Science) (in Chinese) (in Chinese) , 2006, 30(1): 33-37.
[5] 沈凤, 钱绍新, 刘雯林, 等. AVO多参数反演在气层检测和储层定量研究中的应用. 石油学报 , 1994, 15(2): 11–20. Shen F, Qian S X, Liu W L, et al. The application of AVO multiple parameter inversion in gas sand detection and reservoir quantitative examination. Acta Petrolei Sinica (in Chinese) (in Chinese) , 1994, 15(2): 11-20.
[6] [JP2]HampsonD. AVO inversion, theory and practice. Geophysics , 1991, 10(6): 39–42.
[7] 王天禧. 用平面波延拓方程进行地震数据的叠前速度反演. 地球物理学报 , 1993, 36(3): 388–395. Wang T X. Prestack velocity inversion of seismic data using. Chinese J| Geophys| (Acta Geophysica Sinica) (in Chinese) (in Chinese) , 1993, 36(3): 388-395.
[8] 岳碧波, 彭真明, 洪余刚, 等. 基于粒子群优化算法的叠前角道集子波反演. 地球物理学报 , 2009, 52(12): 3116–3123. Yue B B, Peng Z M, Hong Y G, et al. Wavelet inversion of prestack seismic anglegather based on particle swarm optimization algorithm. Chinese J| Geophys| (in Chinese) (in Chinese) , 2009, 52(12): 3116-3123.
[9] 杨培杰, 印兴耀. 非线性二次规划贝叶斯叠前反演. 地球物理学报 , 2008, 51(6): 1876–1882. Yang P J, Yin X Y. Nonlinear quadratic programming. Chinese J| Geophys| (in Chinese) (in Chinese) , 2008, 51(6): 1876-1882.
[10] 路鹏飞, 杨长春, 郭爱华, 等. 改进的模拟退火算法及其在叠前储层参数反演中的应用. 地球物理学进展 , 2008, 23(1): 104–109. Lu P F, Yang C C, Guo A H, et al. Modified simulated annealing algorithm and its application in prestack inversion of reservoir parameters. Progress in Geophysics (in Chinese) , 2008, 23(1): 104-109.
[11] 黄伟传, 杨长春, 王彦飞. 利用叠前地震数据预测裂缝储层的应用研究. 地球物理学进展 , 2007, 22(5): 1602–1606. Huang W C, Yang C C. The application of prestack seismic data in predicting the fractured reservoir. Progress in Geophysics (in Chinese) , 2007, 22(5): 1602-1606.
[12] 尚永生, 杨长春, 王真理, 等. 塔里木盆地卡4区块AVO研究. 地球物理学进展 , 2007, 22(5): 1408–1415. Shang Y S, Yang C C, Wang Z L, et al. AVO research in Ka4 region of Tarim basin. Progress in Geophysics (in Chinese) (in Chinese) , 2007, 22(5): 1408-1415.
[13] Hammersley J M, Handscomb D C. Monte Carlo Methods. New York: Chapman and Hall , 1964.
[14] Grandis H, Menvielle M, Roussignol M. Bayesian inversion with Markov chainsI The magnetotelluric onedimensional case. Geophys| J| Int| , 1999, 138(3): 757-768.
[15] Malinverno A, Leaney S. A Monte Carlo method to quantify uncertainty in the inversion of zerooffset VSP data. In: SEG 70th Annual Meeting, Calgary, Alberta , 2000.
[16] Chen J S. Joint inversion of seismic AVO and EM data for gas saturation estimation using a samplingbased stochastic model. In: SEG 74th Annual Meeting, Denver, Colorado , 2004.
[17] Chen J S, Hoversten G M, Vasco D, et al. A Bayesian model for gas saturation estimation using marine seismic AVA and CSEM data. Geophysics , 2007, 72(2): 85-95.
[18] [CM(26]HongT C, SenM K. A new MCMC algorithm for seismic waveform inversion and corresponding uncertainty analysis. Geophys| J| Int| , 2009, 177(1): 14–32.
[19] Wang D Y, Zhang G Z. AVO simultaneous inversion using Markov Chain Monte Carlo. In: SEG 80th Annual Meeting, Denver, Colorado , 2010.
[20] KeilisBorok V I, Yanovskaya T B. Inverse problems of seismology. Geophys| J| , 1967, 13(1/3): 223-234.
[21] Press F. Earth models obtained by Monte Carlo inversion. J| Geophys| Res| , 1968, 73(16): 5223-5234.
[22] 王辉, 张玉芬. 基于模型的叠前数据多参数非线性反演. 岩性油气藏 , 2008, 20(2): 108–113. Wang H, Zhang Y F. Modelbased multiparameters nonlinear inversion method in prestack domain. Lithologic Reservoirs (in Chinese) (in Chinese) , 2008, 20(2): 108-113.
[23] Malinverno A. Parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem. Geophys| J| Int| , 2002, 151(3): 675-688.
[24] Tarantola A. Inverse Problem Theory. New York: Elsevier Sci , 1987.
[25] 赵 琪. MCMC方法研究. 山东: 山东大学 [博士论文], 2007. 11~18
[26] Ross S M. 统计模拟. 北京: 人民邮电出版社, 2007 : 203 -206.
[27] 朱 嵩. 基于贝叶斯推理的环境水力学反问题研究. 浙江: 浙江大学 [博士论文], 2008. 63~90
[28] Fatti J L, Smith G C, Vail P J, et al. Detection of gas in sandstone reservoirs using AVO analysis: A 3D seismic case history using the Geostack technique. Geophysics , 1994, 59(9): 1362-1376. DOI:10.1190/1.1443695
[29] 郑晓东. Zoeppritz方程的近似及其应用. 石油地球物理勘探 , 1991, 26(2): 129–144.
[30] Castagna J P, Batzle M L, Eastwood R L. Relationships between compressionalwave and shearwave velocities in clastic silicate rocks. Geophysics , 1985, 50(4): 571-581. DOI:10.1190/1.1441933
[31] 马中高, 解吉高. 岩石的纵、横波速度与密度的规律研究. 地球物理学进展 , 2005, 20(4): 905–910. Ma Z G, Xie J G. Relationship among compressional wave, shear wave velocities and density of rocks. Progress in Geophysics (in Chinese) , 2005, 20(4): 905-910.