基于贝叶斯理论的叠前多波联合反演弹性模量方法

侯栋甲，刘洋＊，胡国庆，魏修成，陈天胜

1中国石油大学（北京）油气资源与探测国家重点实验室，北京 102249

2中国石油大学（北京）CNPC物探重点实验室，北京 102249

3新疆油田勘探开发研究院，克拉玛依 834000

4中国石化勘探开发研究院，北京 100083

1 引言

随着油气资源勘探程度的增加，勘探目标已由构造型油气藏转向岩性油气藏，岩性识别及油气检测显得尤为重要.相对于叠后地震数据反演而言，叠前反演保留了地震反射振幅随偏移距变化特征，能够得到密度、纵波速度、横波速度、泊松比等多种弹性参数，因而能更可靠地揭示地下储层的物性及含油气性（Downton，2005）.由于叠前反演问题都是高维的和非适定的，因此获得可靠稳定的解对叠前反演至关重要.多波地震勘探技术的出现为提高反演结果的精度、提高岩性识别和油气预测的可靠性提供了支持.随着多波地震勘探技术的发展，充分发挥多波资料含有丰富岩性信息的优势，利用多波资料反演地层岩性参数已成为一项重要的研究课题（赵邦六，2007）.

国内外许多学者对多波资料叠前联合反演做了大量的研究.在国外，Stewart（1990）首次给出了实际的纵、横波速度联合反演方法，他将Smith和Gidlow（1987）提出的用于纵波的加权叠加方法推广到联合反演，获得到了纵波速度比和横波速度比等弹性参数.Larsen等（1999）在Stewart（1990）研究的基础之上做了更进一步的研究，提出了同时反演纵波波阻抗和横波波阻抗的方法.Jin（1999）分析了纵、横波联合反演的稳定性，并指出单一纵波反演存在不稳定.Margrave等（2001）利用加权叠加的方法从PP波和PS波振幅资料反演出波阻抗比和泊松比.Zhang等（2003）对Pikes Peak油田的多波地震数据体进行了联合反演，得到了良好的应用效果.Downton和Lines（2004）考虑到地震子波的影响引入了地震褶积模型，利用PP波资料做了基于贝叶斯原理的三参数反演方法研究.Veire和Landrø（2006）利用奇异值分解技术进行了带限的纵、横波联合三参数反演，并且给出了实际数据的应用效果.在国内，陈天胜等（2006）提出了一种基于方向加速度最优方法的速度比值扫描的纵、横波联合反演方法.印兴耀等（2007，2008a，2008b）提出了叠前三参数反演的新方法，该方法在反演过程中引入贝叶斯理论，并假设待求参数服从改进的Cauchy分布以改善反演效果.张春涛等（2010）通过PP波PS波AVO联合二参数反演进行了波阻抗反演研究.李录明等（2010）进行了各向异性介质的三维纵波、横波叠前联合反演方法研究及应用.胡国庆等（2011）利用纵波和转换波进行联合反演来估算密度比和速度比，在反演过程中利用贝叶斯理论加入先验信息，改善反演的不适定性.魏超等（2011）提出一种反向加权系数的方法来提高非线性AVO反演的精度.邴萍萍等（2011）提出利用支持向量机来求解AVO非线性反演问题，获得高精度的岩性参数.宗兆云等（2012a，2012b）提出一种反演杨氏模量和泊松比的方法来识别页岩气“甜点”.

Yin等（2012）提出了一种利用纵波资料、基于贝叶斯理论来反演密度比和模量比的方法.鉴于模量比与流体因子有关，可以更好地指示油气；基于速度比和模量比之间的关系，将Zoeppritz方程的近似形式Aki-Richards公式中的速度比替换为模量比，可以得到与模量比有关的反射系数公式，从而来反演密度比等参数.本文在Yin等（2012）研究的基础上，提出了一种联合纵波和转换波反演密度比和模量比的方法，利用多波资料联合反演密度比、剪切模量比、体积模量比三参数.在反演过程中引入贝叶斯理论，假定先验信息服从高斯分布，待求参数服从改进的Cauchy分布，并去除待求参数之间的相关性.理论模型数据和实际数据的测试结果表明该方法可以有效提高反演结果的稳定性、精度和抗噪能力.

2 基本原理

2.1 AVO近似方程

根据地震波传播理论，当纵波非垂直入射到一个波阻抗界面时，会产生反射纵波、反射转换横波、透射纵波、透射转换横波四种波.反射系数和透射系数可以用Zoeppritz方程来精确表述（Zoeppritz，1919），在上下层介质的密度和速度差异不太大的情况下，可以将Zoeppritz方程近似简化为Aki-Richards公式（Aki和Richards，1980），纵波反射系数和转换波反射系数可分别表述为以下形式：

其中，rpp表示纵波反射系数，rps表示转换波反射系数，

根据文献Yin等（2012）的分析，在小角度入射的情况下，公式（7）和（8）具有较高的精度，其与精确解的误差可以忽略不计.

2.2 多波联合反演方程

为了从地震记录中反演出密度比和模量比，应该先建立多波联合反演方程.首先利用最小二乘原理构建多波联合反演目标函数（Veire和Landrø，2006）：

PP波记录经过反褶积、叠前偏移等之后的结果，rpp、rps分别表示PP波和PS波反射系数序列，i、j分别表示PS波和PP波资料道序号，M表示参加反演的转换波道数，N表示参加反演的纵波道数，η（0≤η≤1）是转换波加权因子，通过分析PS波和PP波资料的品质可以确定加权因子大小，品质好的资料加权系数相应较大.为使目标函数取得最小值，将目标函数对待求参数求导，并令其为零，可得

其中G是一个3×3的矩阵，D是一个3维列向量，其元素分别为

X是一个3维列向量，

由于待求参数之间是统计相关的，因此需要应用三者之间的协方差矩阵来去除参数之间的相关性.首先将待反演参数的协方差矩阵表示为

协方差矩阵的主对角线上的元素表示每一个随机变量的方差，非对角线上的元素是两个变量的协方差.对协方差矩阵（21）进行奇异值分解可得

利用分解的结果，将方程（10）进行如下变换：

变换后待求参数之间是相互独立的，则可得

一般来说，在给定初值的情况下，采用合适的反演算法循环迭代求解方程（24），即可得到＾X，通过转换就能得到待求参数X.但为了使结果更可靠，反演更稳定，利用贝叶斯方法将估计的先验信息加入反演公式中，以克服反演问题的不适定性.

2.3 基于贝叶斯理论的反演

基于贝叶斯理论，假设不同测量条件得到的测量结果的高斯分布是独立的，利用乘积准则得到似然函数的解为

其中，δer为误差的标准差，M+N为地震记录的总道数.

对于三参数反演的情况，改进的Cauchy分布描述待求参数具有更高的分辨率（Alemie和Sacchi，2012），故采用（26）式来描述待求参数的分布：

根据贝叶斯理论，待求参数的后验概率分布为

对式（28）右边取对数可得目标函数为

其中，λ为常数加权因子，Q是对角加权矩阵：

应用预条件共轭梯度法迭代求解非线性方程组（30），即可以得到要反演的参数，但要得到最终的结果还需做如下转化：

对实际地震资料进行叠前反演之前，需要对叠前地震数据进行保幅处理，包括多次波衰减、几何扩散校正、振幅补偿、地表一致性反褶积、保幅叠前偏移等，并假设处理后层间多次波、各向异性的影响可以忽略不计.偏移距道集要通过共角度偏移或射线追踪方法转换为角道集.Buland和More（2003）用模型证明，偏移距道集到角道集之间的转换误差可以忽略不计.

3 模型数据测试及实际资料反演

3.1 模型参数及模型数据测试

将某口井的实测数据从深度域转换到时间域，将转换后的密度和速度作为地层模型参数，对模型参数进行平滑处理即可得到初始模型.地层模型参数如图1（a，b，c）所示，分别为地层模型的密度、横波速度、纵波速度，其中蓝线表示地层模型参数的真实值，绿线表示经过平滑处理的结果，将其作为初始模型.

采用图1所示的地层模型参数真实值，依据Snell定律，将炮检距通过射线追踪转换为入射角.利用精确的Zoeppritz方程计算出纵波反射系数和转换波反射系数.制作合成记录的主要参数为：偏移距为60m，道间距为60m，一共40道；子波为主频为40Hz的Ricker子波；时间采样率为1ms.通过褶积模型得到经过动校正的多波叠前道集数据如图2（a，b）所示，分别为转换波剖面和纵波剖面.截取反射时间为500～800ms的反射层作为反演目的层段.对于模型数据而言，考虑到两种资料所受的影响因素是相同的，品质是等同的，故在进行联合反演时，将纵、横波加权系数均取为0.5.

图1 地层模型参数.（a）密度；（b）横波速度；（c）纵波速度Fig.1 Stratum model parameters.（a）Density；（b）S-wave velocity；（c）P-wave velocity

对合成地震记录，分别采用纵波反演、纵波和转换波联合反演方法来估算三参数.纵波反演的流程和方法与联合反演是一致的，将转换波加权因子η取为0即可.初始模型采用图1绿线所示的结果，通过射线追踪计算角度.反演的结果如图3、4、5所示，分别为密度比、剪切模量比、体积模量比，其中第一列的蓝线表示纵波和转换波联合反演的结果，第二列的黑线表示单独利用纵波反演的结果，第三列的红线表示待反演参数的真实值，最后一列表示反演结果的绝对误差，其中蓝线表示联合反演结果的误差，黑线表示纵波反演结果的误差.从反演结果误差统计来看，联合反演结果的误差要明显小于单一纵波反演，精度明显提高，且联合反演的条件数减小，因而反演过程更加稳定，降低了反演的不适定性.

图2 合成地震记录.（a）转换波剖面；（b）纵波剖面Fig.2 Synthetic data.（a）S-wave section；（b）P-wave section

图3 密度比.（a）联合反演；（b）纵波反演；（c）由模型参数计算的真实值；（d）误差比较Fig.3 Density contrast.（a）Joint inversion result；（b）PP inversion result；（c）Real value；（d）Error comparison

我们利用公式（1）和（2），采用联合反演的方法计算出密度比、横波速度比、纵波速度比，然后分别利用公式（5）和（6）换算出剪切模量比、体积模量比，其中剪切模量比通过密度比和横波速度比得到，体积模量通过密度比和纵波速度比得到.通过该方法得到的密度比精度明显低于本文联合反演的结果，且在换算模量比时会将密度比和速度比的误差累加，导致换算出的结果精度降低.结果如图6、7、8所示，分别为直接反演出的密度比、换算出的剪切模量比和体积模量比，其中各图中的（a）表示换算出的参数结果，（b）是参数真实值，（c）是两者之间的误差.由图可见，相对联合反演的结果而言，通过替换的方法得到的结果误差总体较大.为了定量比较二者之间的误差，我们计算了误差能量.联合反演的密度比误差能量为0.1093，直接反演的密度比误差能量为0.1990，联合反演的剪切模量比误差能量为0.0696，转换得到的剪切模量比误差能量为0.3479，联合反演的体积模量比误差能量为0.0474，转换得到的体积模量比误差能量为0.1065.可见联合反演的结果其误差能量要比转换的结果小.

图4 剪切模量比.（a）联合反演结果；（b）纵波反演结果；（c）由模型参数计算的真实值；（d）误差比较Fig.4 S-wave molulus contrast.（a）Joint inversion result；（b）PP inversion result；（c）Real value；（d）Error comparison

图5 体积模量比.（a）联合反演；（b）纵波反演；（c）由模型参数计算的真实值；（d）误差比较Fig.5 P-wave modulus contrast.（a）Joint inversion result；（b）PP inversion result；（c）Real value；（d）Error comparison

3.2 含有噪声的合成数据测试

对合成的地震数据加入随机噪声，将信噪比设为2，加入噪声的地震记录如图9（a，b）所示，分别为转换波剖面和纵波剖面.分别采用纵波反演方法、纵波与转换波联合反演方法来估算三参数.反演的结果如图10、11、12所示，分别为密度比、剪切模量比、体积模量比，各条线的意义与图3、4、5中相同.由图可见，在信噪比较低的情况下，联合反演依然能较好地估算出地层模型参数.从反演结果误差统计来看，联合反演的误差要明显小于单一纵波反演，抗噪性能和精度明显优于常规纵波反演.

3.3 实际资料反演

实际数据为某工区多波地震资料.图13（a，b）为抽取的经过动校正之后的叠前转换波数据道集和纵波数据道集.在反演之前先进行纵波和转换波在时间域内的匹配：利用速度分析获得的纵波和转换波速度资料，根据转换波反射时间与纵波反射时间之间的关系，将转换波资料从转换波时间域转换到纵波时间域.将纵波反射时间为300～1100ms的层位作为反演的目标层位.根据实际资料的品质，在联合反演时，纵波加权系数为0.7，转换波加权系数为0.3.在此剖面上有一口井，该井测井数据作为待求参数的先验信息，并用于验证反演效果.

图6 密度比.（a）直接计算的结果；（b）由模型参数计算的真实值；（c）误差Fig.6 Density contrast.（a）Calculated result；（b）Real value；（c）Error

图7 剪切模量比.（a）转换计算的结果；（b）由模型参数计算的真实值；（c）误差Fig.7 S-wave molulus contrast.（a）Converted result；（b）Real value；（c）Error

图8 体积模量比.（a）转换计算的结果；（b）由模型参数计算的真实值；（c）误差Fig.8 P-wave modulus contrast.（a）Converted result；（b）Real value；（c）Error

对实际数据分别采用纵波反演方法、纵波和转换波联合反演方法来估算三参数.反演结果如图14、15所示，分别为联合反演和纵波单独反演的结果.由图可见，联合反演对目标层位刻画更清晰.表1给出了在井点位置处400～800ms内联合反演和单独反演结果与实际井数据的误差均方根能量，可见联合反演精度高于纵波反演，与实际井数据的结果也较吻合.

表1 井位处反演结果误差比较Table 1 Error comparison of inversion results at the well location

图9 加入噪声的合成地震记录 （SNR＝2）.（a）转换波剖面；（b）纵波剖面Fig.9 Synthetic data with noise（SNR＝2）.（a）S-wave section；（b）P-wave section

图10 密度比.（a）联合反演；（b）纵波反演；（c）由模型参数计算的真实值；（d）误差比较Fig.10 Density contrast.（a）Joint inversion result；（b）PP inversion result；（c）Real value；（d）Error comparison

图11 剪切模量比.（a）联合反演；（b）纵波反演；（c）由模型参数计算的真实值；（d）误差比较Fig.11 S-wave modulus contrast.（a）Joint inversion result；（b）PP inversion result；（c）Real value；（d）Error comparison

图12 体积模量比.（a）联合反演；（b）纵波反演；（c）由模型参数计算的真实值；（d）误差比较Fig.12 P-wave modulus contrast.（a）Joint inversion result；（b）PP inversion result；（c）Real value；（d）Error comparison

4 结论

本文发展了一种基于贝叶斯理论的多波联合叠前反演密度比和模量比的方法.该方法具有以下几个方面的特点：

（1）反演的密度比、剪切模量比、体积模量比三参数具有重要的地质意义，可以更好地识别岩性以及指示油气；

（2）相对于常规纵波反演方法，联合反演方法精度更高、稳定性更好、抗噪能力更强；

（3）在反演过程中引入贝叶斯理论和加入井约束先验信息，降低了反演多解性，改善了反演不适定性.

模型数据测试和实际资料应用验证了该方法的可行性和有效性.

图13 经过动校正的实际叠前地震数据.（a）转换波道集；（b）纵波道集Fig.13 Real prestack seismic data after NMO.（a）S-wave data set；（b）P-wave data set

Aki K，Richards P G.1980.Quantitative Seismology.New York：W.H.Freeman ＆Co.

Alemie W，Sacchi M D.2012.High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution.Geophysics，76（3）：R43-R55.

Bing P P，Cao S Y，Lu J T.2012.Non-linear AVO inversion based on support vector machine.Chinese J.Geophys.（in Chinese），55（3）：1025-1032.

Buland A，Omre H.2003.Bayesian linearized AVO inversion.Geophysics，68（1）：185-198.

Chen J J，Yin X Y.2007.Three-parameter AVO waveform inversion based on Bayesian theorem.Chinese J.Geophys.（in Chinese），50（4）：1251-1260.

Cheng T S，Liu Y，Wei X C.2006.Joint amplitude versus offset inversion of P-P and P-SV seismic data.Journal of China University of Petroleum（in Chinese），30（1）：33-37.

Downton J E，Lines L R.2004.Three-term AVO waveform inversion.∥74th Annual International Meeting.SEG，Expanded Abstracts，215-219.

Downton J E.2005.Seismic Parameter Estimation from AVO Inversion.Calgary：University of Calgary.

Hu G Q，Liu Y，Wei X C，etal.2011.Joint PP and PS AVO inversion based on Bayes theorem.Applied Geophysics，8（4）：293-302.

Jin S.1999.Characterizing reservoir by using jointly P-and S-wave AVO analysis.∥Expanded Abstracts of 69th Annual Internat SEG Mtg.Houston：Expanded Abstracts of 69th Annual，687-690.

Larsen J，Margrave G F，Lu H X.1999.AVO analysis by simultaneous P-P and P-S weighted stacking applied to 3-D seismic data.∥Expanded Abstracts of 69th Annual Internat SEG Mtg.Houston：The Society of Exploration Geophysicists，721-724.

Li L M，Luo S X，Wang M C，etal.2010.The joint inversion method on anisotropy medium and its application.Oil Geophysical Prospecting（in Chinese），45（1）：60-65.

Margrave G F，Stewart R，Larsen J A.2001.Joint PP and PS seismic inversion.The Leading Edge，20（9）：1048-1052.

Smith G C，Gidlow P M.1987.Weighted stacking for rock property estimation and detection of gas.Geophysical Prospecting，35（9）：993-1014.

图15 纵波反演结果（a）密度比；（b）剪切模量比；（c）体积模量比.Fig.15 PP inversion result（a）Density contrast；（b）S-wave modulus contrast；（c）P-wave modulus contrast.

Stewart R.1990.Joint P and P-SV inversion.The CREWES Project Research Report，2：112-115.

Veire H H，LandrøM.2006.Simultaneous inversion of PP and PS seismic data.Geophysics，71（3）：1-10.

Wei C，Zheng X D，Li J S.2011.A study on nonlinear AVO inverse method.Chinese J.Geophys.（in Chinese），54（8）：2110-2116.

Yang P J，Yin X Y.2008.Nonlinear quadratic programming Bayesian prestack inversion.Chinese J.Geophys.（in Chinese），51（6）：1876-1882.

Yin X Y，Yang P J，Zhang G Z.2008.A novel prestack AVO inversion and its application.78th Annual International Meeting，SEG，Expanded Abstracts，2041-2045.

Zhang C T，Wang S X，Li S J，etal.2010.Compressional and shear wave joint inversion technique research and its application.Oil Geophysical Prospecting（in Chinese），45（4）：520-524.

Zhang H，Margrave G F，Brown R J.2003.Joint PP-PS inversion at Pikes Peak oilfield，Saskatchewan.CREWES Research Report，15：1-12.

Zhao B L.2007.The Theory and Practice of Multi-component Seismic Survey Technology（in Chinese）.Beijing：Petroleum Industry Press.

Zoeppritz K.1919.Erdbebenwellen VIII B，Uber Reflexion und Durchgang seismischer Wellen durch Unstetigkeitsflachen.Gottinger.Nachr.，1：66-84.

Zong Z Y，Yin X Y，Wu G C.2012.AVO inversion and poroelasticity with P-and S-wave moduli.Geophysics，77（6）：N17-N24.

Zong Z Y，Yin X Y，Wu G C.2012a.Fluid identification method based on compressional and shear modulus direct inversion.Chinese J.Geophys.（in Chinese），55（1）：284-292.

Zong Z Y，Yin X Y，Zhang F，etal.2012b.Reflection coefficient equation and prestack seismic inversion with Young′s modulus and Poisson ratio.Chinese J.Geophys.（in Chinese），55（11）：3786-3791.

附中文参考文献

邴萍萍，曹思远，路交通.2011.基于支持向量机的非线性AVO反演.地球物理学报，55（3）：1025-1032.

陈建江，印兴耀.2007.基于贝叶斯理论的AVO三参数波形反演.地球物理学报，50（4）：1251-1260.

陈天胜，刘洋，魏修成.2006.纵波和转换波联合AVO反演方法研究.中国石油大学学报，30（1）：33-37.

李录明，罗省贤，王明春等.2010.各向异性介质三维纵横波联合叠前反演方法及应用.石油地球物理勘探，45（1）：60-65.

魏超，郑晓东，李劲松.2011.非线性AVO反演方法研究.地球物理学报，54（8）：2110-2116.

杨培杰，印兴耀.2008.非线性二次规划贝叶斯叠前反演.地球物理学报，51（6）：1876-1882.

张春涛，王尚旭，李生杰等.2010.纵横波联合反演方法研究及应用.石油地球物理勘探，45（4）：520-524.

赵邦六.2007.多分量地震勘探技术理论与实践.北京：石油工业出版社.

宗兆云，印兴耀，吴国忱.2012a.基于叠前地震纵横波模量直接反演的流体检测方法.地球物理学报，2012a，55（1）：284-292.

宗兆云，印兴耀，张峰等.2012b.杨氏模量和泊松比反射系数近似方程及叠前地震反演.地球物理学报，55（11）：3786-3791.