地球物理学报  2019, Vol. 62 Issue (1): 276-288   PDF    
基于改进反射系数近似方程的纵横波阻抗同步反演
付欣1,2, 张峰1,2, 李向阳1,2, 钱忠平3, 陈海峰3, 梅璐璐3     
1. 中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249;
2. 中国石油大学(北京)CNPC物探重点实验室, 北京 102249;
3. 中国石油天然气集团公司东方地球物理勘探有限责任公司, 河北涿州 072751
摘要:基于地震波反射系数近似公式的叠前反演是油气勘探的重要工具.本文在已有研究的基础上,推导了一个改进的射线参数域地震纵波反射系数近似方程.该方程建立了地震纵波反射系数与纵波阻抗和横波阻抗的非线性关系,在中、小角度的范围内较现有的反射系数线性近似公式精度更高.另外,由于该方程仅包含纵波和横波阻抗反射系数项,因此基于新方程的反演能够有效地降低同步反演纵波速度、横波速度、密度三个参数的不适定性.在此基础上,结合广义线性反演法(GLI)理论和贝叶斯理论,相应地发展了一种叠前地震同步反演方法.模型测试和实际资料的应用表明,基于新方程的反演方法能够利用有限角度(偏移距)的数据稳定地反演纵波和横波阻抗,由于在反演过程中,不需要假设纵横波速度为常数,因此该方法还能有效地提高反演结果的精度.
关键词: 反射系数      叠前反演      纵波阻抗      横波阻抗      贝叶斯广义线性反演     
Simultaneous inversion of P-and S-wave impedances based on an improved approximation equation of reflection coefficients
FU Xin1,2, ZHANG Feng1,2, LI XiangYang1,2, QIAN ZhongPing3, CHEN HaiFeng3, MEI LuLu3     
1. State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China;
2. CNPC Key Laboratory of Geophysical Prospecting, China University of Petroleum, Beijing 102249, China;
3. BGP INC., China National Petroleum Corporation, Heibei Zhuozhou 072751, China
Abstract: Prestack inversion utilizing approximate formulas of reflection coefficients is an important tool for exploration of oil and gas. In this paper, we derivate an improved approximation of reflection coefficients in the ray parameter domain. This equation provides a nonlinear relationship between seismic reflection coefficients and P-and S-wave impedance reflectivity, and is more accurately than the conventional linearized equations in the range of small and middle incident angles. Besides, because the new equation contains only the terms of P-and S-wave impedances, the inversion utilizing this equation can efficiently reduce the ill-posedness of the simultaneous inversion of P-wave velocity, S-wave velocity and density. We then develop a prestack seismic inversion method by combining the generalized linear inversion (GLI) method and the Bayes theory. Both of synthetic and real data tests show the feasibility of this inversion scheme, and for the reason of the assumption that P-and S-wave velocity ratio is constant does not use in the inversion process, the accuracy of inverted results is improved effectively.
Keywords: Reflectivity    Prestack inversion    P-wave impedance    S-wave impedance    Bayes generalized linear inversion    
0 引言

叠前同步反演作为计算岩石物性参数的重要技术, 已广泛运用于油气的勘探与开发中.相比于其他反演方法, 如全波形反演等, 基于地震振幅随偏移距变化(AVO)关系的叠前反演, 能够综合测井、岩石物理等多尺度资料, 在计算速度方面也具有明显的优势.AVO分析和反演为油气的识别和储层的评价提供了很大的帮助, 所以在生产应用中备受青睐.

叠前AVO反演技术经历了长时间的发展, 已经奠定了较为坚实的基础.精确的反射系数公式最初由Zoeppritz(1919)提出, 之后便作为AVO反演的基础理论.

许多学者通过对Zoeppritz方程进行简化, 得到不同形式的反射系数近似公式, 并相应地发展了AVO反演方法(Bortfeld, 1961Richards and Frasier, 1976Aki and Richards, 1980Mallick, 1993Wang, 1999Yin and Zong, 2013).Aki和Richards(1980)提出了包含纵波速度反射率, 横波速度反射率和密度反射率的三项线性近似公式, 在弱阻抗差的情况下具有较高的精度, 但是基于该公式的AVO反演具有较大程度的不适定性.很多学者在此基础上进行了改进研究, 并拓展了待反演参数的种类.Wiggins等(1983)对Aki-Richards的三项线性近似公式进行了重新排列, 将其写成关于截距(A), 梯度(B)和曲率(C)的三项线性方程.Shuey(1985)将横波反射率项用泊松比来替代, 也得到了与Wiggins类似的近似方程.Smith和Gidlow(1987)以及Fatti等(1994)都得到关于纵、横波阻抗的近似方程, 减少了待求参数的个数, 降低了反演的不适定性.

Smith和Gidlow(1987)通过引入Gardner方程(Gardner et al., 1974)来消除密度项, 之后Simmons和Backus(1996)同时应用Castagna方程(Castagna et al., 1985)和Gardner方程(Gardner et al., 1974)来改善反演问题的稳定性, 进行线性AVO反演.Goodway等(1997)Zong等(2012a)推导了关于弹性模量的线性AVO近似.Gray等(1999)提出了两个线性方程, 建立了反射系数与体积模量, 剪切模量, 以及反射系数与拉梅参数之间的关系.Russell等(2011)利用孔隙介质理论, 讨论了流体因子与反射系数之间的关系.Zong等(宗兆云等,2012; Zong et al., 2012b)以及Zong和Yin(2016), 直接反演了杨氏模量和泊松比, 用于对页岩储层的识别和评价, 以及流体的判定.越来越多的岩石物性参数与反射系数之间的关系被推导出(Zong et al., 2013Yin and Zhang, 2014; Yin et al., 2014Li et al., 2014Zong et al., 2015).

但AVO反演目前仍存在一些问题, 主要包括以下三个方面:(1)反射系数近似公式在大入射角、强阻抗差等情况下, 与精确公式之间还存在一定误差;(2)在反演过程中, 需要预先估算纵、横波速度比;(3)在同步反演三个参数(例如纵波速度、横波速度、密度)时, 反演问题病态程度较高, 需要合理的正则化.为解决上述问题, 笔者在现有研究的基础上, 提出了一个改进的地震纵波反射系数近似方程, 与现有的反射系数线性近似方程相比, 该方程在中、小角度范围内具有更高的精度, 能够很好地表示四种不同储层类型的AVO响应.在进行AVO反演时, 不需要假设纵横波速度比是常数, 避免了纵横波速度比误差对反演结果的影响;改进的AVO方程还降低了反演时灵敏度矩阵条件数的大小, 减小了反演的病态程度.基于改进的近似方程, 还发展了相应的反演方法, 模型测试和实际数据应用都体现了本文方法的稳定性和可靠性.

1 改进的反射系数方程

在相邻的地下两层各向同性介质界面上, 纵波反射系数近似公式以射线参数的形式可表示为包含Rf(p)和Rg(p)的两项关系式(Mallick, 1993Wang, 1999):

(1)

式中, 射线参数p=sinθi/αi=sinφi/βiRf(p(θi))=也被称流体项, 表示第i个反射界面上与孔隙流体相关的反射系数, 其中θi, θi+1分别表示纵波入射角和波透射角, αi, αi+1, βi, βi+1, ρi, ρi+1分别为反射面上层和下层的纵波速度、横波波速度、密度;也被称刚体项, 其中Δμβ2Δρ+2ρβΔβ是反射介面上下层介质的剪切模量之差, β=(βi+1+βi)/2, ρ=(ρi+1+ρi)/2, Δβ=βi+1-βi, Δρ=ρi+1-ρiφ=(φi+1+φi)/2为P-SV波反射角φi与P-SV波透射角φi+1的平均.

根据Potter等(1998), 可将密度与横波速度的反射系数之比表示为r=(Δρ/ρ)/(Δβ/β), 同时结合近似公式(Aki and Richards, 1980Wang, 1999), 则Rg(p)可表示为

(2)

根据Ma(2003)Wang(2003)可得以下近似:

(3)

式中Δcosφ=cosφi+1-cosφ, cosφ=(cosφi+1+cosφi)/2.将式(3)代入方程(2)中, 可得对数形式刚体项:

(4)

对于某个反射界面上下层介质而言cos2φ为定值, 因此对数项内的cos2φ为一常数.为进一步化简该公式, 将对数项内分子中的cos2φ近似为cos2φi+1, 分母中的cos2φ近似为cos2φi:

(5)

根据斯奈尔定理:

(6a)

(6b)

再将式(6a)和(6b)代入方程(5)中可得:

(7)

最后将式(7)中, 横波与纵波速度之比βi+1/αi+1βi/αi分别写成横波阻抗SIi+1(βi+1ρi+1)与纵波阻抗AIi+1(αi+1ρi+1)之比SIi+1/AIi+1, 以及SIi(βiρi)与AIi(αiρi)之比SIi/AIi的形式, 则第i个反射界面的反射系数方程可写作纵波阻抗与横波阻抗的函数:

(8)

由于(8)式建立了纵波反射系数与纵波阻抗和横波阻抗之间的关系, 我们称之为ASI方程, 纵波反射系数记作RASI.

2 反射系数近似方程对比分析

为了说明上文推导的反射系数方程的精度, 下面将新方程与已有的经典公式进行对比分析.Fatti等(1994)对Aki-Richards方程(Aki and Richards, 1980)进行整理, 推导出与纵波阻抗、横波阻抗、密度有关的线性反射系数近似方程:

(9)

式中k=β/α, 在反演过程中设定为常数, AI=(AIi+1+AIi)/2, ΔAI=AIi+1-AIi, SI=(SIi+1+SIi)/2, ΔSI=SIi+1-SIi.由于密度项在假设k接近0.5, 入射角度不大的情况下, 可近似抵消, 因此在反演过程中, 也常用两项Fatti近似公式(Downton et al., 2001).其表达式为

(10)

为对新方程进行分析, 建立了如表 1所示的四类AVO模型, 将ASI方程与Zoeppritz精确方程, 公式(1), Fatti近似方程进行比较.笔者从多方面对ASI方程进行了分析, 并与Fatti方程进行比较, 以说明本文推导的ASI方程与传统纵横波阻抗近似方程相比的优劣性.

表 1 四类AVO模型参数 Table 1 Elastic parameters of four classes AVO models
2.1 反射系数近似公式精度对比

图 1(a—d)分别显示了四类AVO模型反射系数随纵波入射角变化的曲线.从图 1a1b中可知, 对于第一类和第二类AVO模型, 方程(1)与Fatti公式非常相似, 均随着入射角的增大而远离精确方程;而ASI方程的误差是先增大后减小的, 在大约40°时, 与真实值相等, 然后继续远离真实值;可以看出在入射角小于40°时, ASI具有较好的精确度.对于第三类AVO模型, 四个近似公式都具有良好的精度, ASI方程与方程(1)的反射系数曲线与精确AVO曲线几乎重合.对于第四类AVO模型,除两项Fatti方程误差较大外,其他方程均具有良好的反射系数表达精度.总体而言, ASI方程与其他几个近似公式相比, 在0~40°的范围内具有更高的精度.

图 1 四种AVO分类模型的反射系数曲线 (a)第一类; (b)第二类; (c)第三类; (d)第四类. Fig. 1 Reflection coefficients of four classes AVO models (a) Class Ⅰ; (b) Class Ⅱ; (c) Class Ⅲ; (d) Class Ⅳ.
2.2 固有参数敏感性的分析

在实际反演过程中, Fatti方程中包含的固有参数k(β/α)需要预先估算并假设为常数, ASI方程中包含的固有参数r((Δρ/ρ)/(Δβ/β))也需要类似的估算.这一部分,笔者讨论了当上述两个固有参数不准确时,分别对Fatti方程和ASI方程反射系数表达精度的影响(使用第三类AVO模型).从图 2中可知, Fatti方程对于k值的变化是非常敏感的, 在入射角为50°时, 28%的k值估算误差引了起近10%的反射系数误差, 而ASI方程, 在r值误差为130%时, 50°处只引起了约1%的误差, 仍能非常准确地表达反射系数.

图 2 第三类AVO模型 利用不同近似方程, 设定不同r值和k值所得的反射系数曲线.真实kr值分别为0.3914和0.1722.(a)Fatti方程, (b)ASI方程. Fig. 2 AVO curves of the Class Ⅲ AVO model for different r and k values The exact values of k and r are 0.3914 and 0.1722 respectively. (a) Fatti equation; (b) ASI equation.

对于一个反演问题:

(11)

其中R表示数据项(反射系数), m为模型参数项, 而F是数据与模型直接的映射关系, 通过分析其映射关系矩阵F可讨论反演问题的不适定性.F也被称为灵敏度矩阵, 其条件数的大小反映了反演方程病态问题的严重程度.上式两边同时对m求导可得灵敏度矩阵的表达式:

(12)

对应三项和两项Fatti方程的灵敏度矩阵分别为

而对于ASI方程而言, 笔者采用了广义线性反演方法(GLI)(在后文中会详细讲解), 因此灵敏度矩阵写作.

图 3a3b中, 笔者用表 1中的第一个模型的参数, 计算了ASI方程和Fatti方程灵敏度矩阵的特征值和条件数(最大特征值与最小特征值之比).在图 3a中, 蓝色曲线表示的是ASI方程的第一、二特征值, 绿色曲线表示的是两项Fatti方程的第一、二特征值, 红色曲线表示的是三项Fatti方程的第一、二、三特征值;图 3b则是各近似方程对应的条件数随最大入射角变化的曲线.由图 3a可知, 除ASI方程的第一特征值不变和第二特征值先增大后减小外,其他各特征值均随最大入射角增大而增大;在图 3b中, 条件数总体趋势均随最大入射角的增大而减小, 意味着, 反演问题的病态性随着最大入射角的增大而减小.值得注意的是,在20°最大入射角附近,ASI方程的第二特征值等于第一特征值,说明此时纵、横波阻抗的反演难度相当,且此时的条件数达到极小值1,反演问题的稳定性最好.从图 3b中可看出, ASI方程灵敏度矩阵的条件数远远低于三项Fatti方程, 比两项Fatti方程也有明显减小.相比于两项Fatti方程,ASI方程对应的条件数在10几度最大入射角之时已下降到较低水平,此后便在低值区域小范围波动,而即使是两项Fatti方程,其对应条件数要达到这样的低值水平,所需要的最大入射角也应超过40°.灵敏度矩阵条件数的降低, 说明ASI方程一定程度上改善了反演问题的病态情况.

图 3 灵敏度矩阵的特征值(a)和条件数(b) Fig. 3 Eigenvalues (a) and condition numbers (b) of sensitivity matrices
3 基于贝叶斯理论的同步反演

由于ASI方程是一个非线性的方程, 在此笔者利用Cook和Schneider(1983)提出的广义线性反演方法(GLI), 对ASI方程进行多元泰勒展开, 并保留一阶项:

(13)

式中R0, AI0SI0分别为初始模型的反射系数, 纵波阻抗和横波阻抗.设ΔR=R-R0, ΔAI=AI-AI0, ΔSI=SI-SI0,则在多入射角资料下可得如下矩阵方程:

(14)

式中ΔR表示的是反射系数的残差项, G是反射系数的雅可比矩阵, Δm是参数的扰动量.

本文采用贝叶斯反演理论对反演方程进行稀疏约束, 假设残差项符合零均值的高斯分布, 扰动项符合改进的Cauchy分布, 以突出弱小反射, 提高反演结果的分辨率(杨培杰和印兴耀, 2008Alemie and Sacchi;2011), 则可构造如下目标函数:

(15)

式中, σn为误差项εR-G·Δm的标准差;M为模型参数扰动量Δm的总长度;σm为扰动量的标准差;Δmi对应第i个待反演参数扰动量.再通过对(15)式求导最小化目标函数, 整理后, 可得最终反演方程:

(16)

式中λ=2σn2/σm2, 称为改进的柯西分布因子;λQ为改进的柯西约束, 用来约束反演参数的稀疏程度, 是在柯西分布的基础上对分母项取平方而得到的, 改进后的对角加权矩阵.最后, 待反演参数可通过下式求得:

(17)

由于ASI方程中, 同时包含了纵波入射角θi和纵波透射角θi+1, 在共射线参数道集中, 上一反射界面对应的透射角与下一反射界面的入射角是相同的.但是一般所采用的叠前反演道集是共反射角域的, 而在共入射角道集中, 上一反射面对应的纵波透射角并不等于下一反射面对应的入射角.所以在获得了每一采样点对应的纵波入射角的同时, 还需要知道其对应的透射角才能进行准确的反演.在实际资料中, 一般可通过下式, 将偏移距域资料转换到入射角域:

(18)

其中θi为第i层对应的纵波入射角, αi为第i层的纵波层速度, αrmsi为第i层的均方根速度, t0为双程走时, x为炮检距.同样的, 也可以结合斯奈尔定理, 从偏移距域资料中得到透射角信息:

(19)

式中θi+1为第i层对应的纵波透射角, αi+1是第i+1层的纵波层速度, 其他参数则与公式(18)是一致的.在将共偏移距道集资料转换成角道集之后, 还需记录每一个采样点对应的透射角信息, 之后便可利用ASI方程进一步开展纵横波阻抗同步反演.

本文在进行反演时采用叠前地震反演方法:(1)先进行部分角道集叠加, 包括共入射角道集叠加和透射角信息的叠加;(2)再对每个部分叠加剖面进行子波估算;(3)然后是模型约束的反射系数反演;(4)最后采用本文方法从反射系数中, 同步反演出纵横波阻抗(详细推导过程见附录A).

4 模型验证

为验证用本文方法进行纵横波阻抗同步反演的稳定性和可行性, 笔者利用北海地区的实际测井数据(Avseth et al., 2005)进行合成数据试算(图 4a, b).图 5a5b显示了在输入数据无噪的情况下, 分别用ASI方程和Fatti方程反演的纵横波阻抗.分别比较了当k=β/αr=(Δρ/ρ)/(Δβ/β)在不同取值的情况下(真实值分别为0.45和0.16)的反演结果.反演中使用的资料为合成的共入射角道集, 所使用的所有初始模型均为用滑动平均滤波法(选取窗口大小为100)由真实模型滤波后得而到的.对于k值和r值的变化, 纵波阻抗反演结果都很稳定, 基本不发生较大改变.但对于横波反演结果, 图 5a中利用Fatti公式在不同的k值情况下反演的横波阻抗差异较大, 说明Fatti方程反演的横波阻抗对k值变化较敏感, 尤其是在阻抗差较大的地方, 如时间大约为2090 ms, 2270 ms和2380 ms处, 其反演结果易出现较大误差;而在图 5b中, 对于不同的r值, ASI方程反演的结果基本不变, 都能很好地与真实值相匹配.在图 6a6b6c中, 使用本文所述方法, 分别用加入不同比例高斯噪声的合成资料进行反演, 使用的合成地震资料信噪比分别为4: 1, 2: 1, 1: 1和0.5: 1.图 6d中显示了, 不同信噪比资料反演结果对应的模型参数的相对误差(反演结果和真实模型的绝对误差与真实模型之比).在无噪声时, 纵横波阻抗反演结果的相对误差较小.而随着噪声的出现, 两者的相对误差水平有所提高, 且随着噪声水平的增加而增加.纵波阻抗反演结果的相对误差增加量相对较小, 在信噪比为0.5时, 总体相对误差仍小于5%.而有噪声时, 横波阻抗反演结果的相对误差相对较大, 在信噪比为0.5时, 相对误差有时接近20%.但总体而言, 在不同信噪比情况下, 都能得到很好的纵横波阻抗反演结果, 说明本文方法具有很好的稳定性和可靠性.

图 4 北海地区实际测井曲线 (a)纵波速度; (b)横波速度; (c)密度; (d)无噪声合成地震记录. Fig. 4 Real model extracted from well data of North Sea and synthetic seismic
图 5 无噪声时, (a)利用Fatti方程, k取不同值的反演结果;(b)利用ASI方程, r取不同值的反演结果 Fig. 5 Free of noise, P- and S-wave impedance inversions (a) utilizing Fatti equation with different k values; (b) utilizing ASI equation with different r values.And the black curves are corresponding real values, the dash blue curves are initial models, and the other curves are inversions.
图 6 不同信噪比(SNR)时, 用本文方法反演的结果 (a) SNR=4; (b) SNR=2; (c) SNR=1; (d) SNR=0.5;其中黑色曲线为真实值, 蓝色虚线为初始模型, 红色曲线为反演结果; (e)不同信噪比时, 反演结果的模型参数的相对误差. Fig. 6 Utilizing the method presented in this paper, P- and S-wave impedance inversions in the case of different signal to noise ratios (SNR) (a) SNR=4; (b) SNR=2; (c) SNR=1; (d) SNR=0.5. The black curves are corresponding real values, the dash blue curves are initial models and the red curves are inversions. And (e) model relative error for different noise levels.

由于上文所述反演方法需要预先建立初始模型, 因此笔者通过以下分析说明该反演方法是否对模型依赖.图 7(a—c)还分别显示了由真实模型经过不同滤波窗口(50, 3和1)所得的初始模型进行反演的结果(信噪比为1).在图 7d中对比了使用不同初始模型进行反演的数据残差随迭代次数的下降曲线.可以看出, 使用不同的初始模型进行反演, 残差的收敛情况并不一样, 但是通过多次的迭代(<30次), 最终都能够获得全局最小值.使用不同初始模型的纵横波反演结果非常相近.结果说明, 当初始模型具备较为可靠的低频信息时, 该反演方法并不易陷入局部极小值.因此在实际地震数据反演中, 需要结合测井数据建立可靠的初始模型, 保证反演结果的准确性.

图 7 使用滑动平均滤波在选取不同窗口大小时获得的初始模型后的反演结果 (a)窗口大小为50;(b)窗口大小为3;(c)窗口大小为1; (d)信噪比为1时, 使用不同模型反演结果的数据残差随迭代次数变化曲线.输入数据信噪比为1时. Fig. 7 In the case of SNR=1, Inversions of using different initial models obtained from utilizing filters to the real model with diverse window sizes of (a) 50, (b) 3 and (c) 1 respectively, (d) Reduction curves of data residual using different initial model.
5 实际资料反演

实际数据来自中国西南地区的四川盆地西部, 储层类型为上三叠统的含气砂岩.图 8a8b8c, 分别是小、中、大三个角度范围的部分叠加剖面, 在该工区有一口井, 测井位置对应剖面上的CDP为1396, 即图 8a中黑色测井线所在位置.井位置1550 ms以下的强反射轴指示储层的顶界面反射信息.图 9a9b分别为反演中所使用的纵波阻抗和横波阻抗初始模型, 图 9c9d为用本文方法, 从大、中、小三个部分叠加剖面中反演出的纵波阻抗剖面和横波阻抗剖面.图 10a10b显示的是井旁道地震反演结果与井数据计算结果, 图中黑色曲线为从井资料中计算出的阻抗值, 蓝色曲线为反演中所使用的初始模型, 红色曲线为反演结果.从图 9c图 9d可看出, 井位置1550 ms以下为储层所在位置, 纵波阻抗和横波阻抗都表现为高值异常, 泥岩层夹杂其中, 表现为阻抗低值异常, 从图 10(a—b)可看出, 在目标层位附近, 纵横波阻抗的反演结果与实际测井资料较为吻合, 体现了本文方法在实际数据应用中的稳定性与可靠性.

图 8 部分叠加剖面(a)小角度(1~12°)、(b)中角度(13~24°)和(c)较大角度(24~35°) Fig. 8 Partial stack sections: (a) small incident angle (1~12°), (b) middle incident angle (13~24°) and (c) large incident angle (24~35°)
图 9 (a) 初始纵波阻抗模型;(b)初始纵横阻抗模型;(c)反演的纵波阻抗剖面; (d)反演的横波阻抗剖面 Fig. 9 (a) Initial P-wave impedance model; (b) Initial S-wave impedance model; (c) Inverted P-wave impedance section; (d) Inverted S-wave impedance section.
图 10 地震反演结果(红色曲线)与井资料计算结果(黑色曲线)对比 Fig. 10 The comparison between inverted results from seismic data and these calculated from well data
6 结论

本文提出了一种改进的反射系数近似方程(ASI方程), 在此基础上结合广义线性反演方法和贝叶斯理论, 发展了对应的纵横波阻抗同步反演方法.相比于经典的反射系数近似公式, ASI方程可以较为准确的表示不同类型的反射系数, 并且对方程中所含的固有参数(r)的误差不敏感, 因此能够更好地进行正演.由于ASI方程只包含纵、横波阻抗两个待求参数, 因此基于ASI方程的同步反演很大程度上降低了常规三参数同步反演问题的不适定性.在一系列合成数据试算和分析的基础上, 将该方法应用于实际地震数据, 验证了该方法的稳定性与可靠性.

附录A

利用Cook和Schneider(1983)提出的广义线性反演方法(GLI), 对ASI方程进行多元泰勒展开, 得线性方程:

(A1)

式(A1)中R0, AI0SI0分别为初始模型的反射系数, 纵波阻抗和横波阻抗.设ΔR=R-R0, ΔAI=AI-AI0, ΔSI=SI-SI0, 在多角度地震资料情况下, 上式可写成以下矩阵形式:

(A2)

上式中, DpDs分别是地震道S对纵波阻抗(AI)和横波阻抗(SI)的偏导;上标θk代表不同的纵波入射角.方程(A2)记作:

(A3)

贝叶斯公式可表示为

(A4)

式中, d为观测数据;I为先验信息;x为待求解模型;P(d|I)为常数, P(x|d, I)为后验概率分布函数, P(d|x, I)是似然函数, 表示模型参数确定时地震记录的条件概率分布, P(x|I)为先验概率分布.假设方程(A3)的误差项εR-G·Δm符合零均值的高斯分布, 则似然函数为

(A5)

上式中, σn为误差项ε的标准差.

相比于高斯分布, 待反演参数的柯西分布更有利于突出弱小反射(杨培杰和印兴耀, 2008Alemie and Sacchi, 2011), 因此我们假设待反演参数的扰动项Δm也符合柯西分布:

(A6)

上式中, M为模型参数总数;σm为扰动量的标准差.由公式(A4)可得参数扰动量的后验概率:

(A7)

最大化后验概率后可得目标函数:

(A8)

再通过对扰动量求导最小化目标函数, 可得

(A9)

其中QM×M的对角加权矩阵, 其第i行, 第i列元素为Qi, i=1/[1+(Δmi)2/σm2]2, Δmi对应第i个待反演参数扰动量, Q矩阵对角元素的值是由柯西分布取平方而来, 称作改进的柯西分布因子(Bogdanov et al., 2007).令式(A9)为零, 可得最终反演方程:

(A10)

λ=2σn2/σm2, 称为柯西分布因子, 上式中λQ为改进的柯西约束, 用来约束反演参数的稀疏程度.最后, 待反演参数可通过下式求得.

(A11)

References
Aki K I, Richards P G. 1980. Quantitative Seismology:Theory and Methods. New York: W. H. Freeman & Co.
Avseth P, Mukerji T and Mavko G. 2005. Quantitative seismic intrepretation.
Alemie W, Sacchi M D. 2011. High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution. Geophysics, 76(3): R43-R55. DOI:10.1190/1.3554627
Bogdanov A, Knudsen L R, Leander G, et al. 2007. PRESENT: An ultra-lightweight block cipher.//International Workshop on Cryptographic Hardware and Embedded Systems. Berlin Heidelberg: Springer, 450-466.
Bortfeld R. 1961. Approximations to the reflection and transmission coefficients of plane longitudinal and transverse waves. Geophysical Prospecting, 9(4): 485-502. DOI:10.1111/gpr.1961.9.issue-4
Castagna J P, Batzle M L, Eastwood R L. 1985. Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks. Geophysics, 50(4): 571-581. DOI:10.1190/1.1441933
Cooke D A, Schneider W A. 1983. Generalized linear inversion of reflection seismic data. Geophysics, 48(6): 665-676. DOI:10.1190/1.1441497
Downton J E, Pickford S, Lines L R. 2001. Constrained three parameter AVO inversion and uncertainty analysis.//71st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 251-254.
Fatti J L, Smith G C, Vail P J, et al. 1994. Detection of gas in sandstone reservoirs using AVO analysis:A 3-D seismic case history using the Geostack technique. Geophysics, 59(9): 1362-1376. DOI:10.1190/1.1443695
Gardner G H F, Gardner L W, Gregory A R. 1974. Formation velocity and density-The diagnostic basics for stratigraphic traps. Geophysics, 39(6): 770-780. DOI:10.1190/1.1440465
Goodway B, Chen T W, Downton J. 1997. Improved AVO fluid detection and lithology discrimination using Lamé petrophysical parameters; "λρ", "μρ", and "λ/μ fluid stack", from P and S inversions.//67th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 183-186.
Gray D, Goodway B, Chen T. 1999. Bridging the gap-Using AVO to detect changes in fundamental elastic constants.//61st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 852-855.
Li C, Yin X Y, Zhang G Z, et al. 2014. Prestack AVA inversion for pore fluid modulus.//74th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 533-537.
Ma J F. 2003. Forward modeling and inversion methods based on generalized elastic impedance in seismic exploration. Chinese Journal of Geophysics, 46(1): 159-168. DOI:10.1002/cjg2.v46.1
Mallick S. 1993. A simple approximation to the P-wave reflection coefficient and its implication in the inversion of amplitude variation with offset data. Geophysics, 58(4): 544-552. DOI:10.1190/1.1443437
Potter C C, Dey A K, Stewart R R. 1998. Density prediction using P-and S-wave sonic velocities.//68th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 320-320.
Richards P G, Frasier C W. 1976. Scattering of elastic waves from depth-dependent inhomogeneities. Geophysics, 41(3): 441-458. DOI:10.1190/1.1440625
Russell B H, Gray D, Hampson D P. 2011. Linearized AVO and poroelasticity. Geophysics, 76(3).
Shuey R T. 1985. A simplification of the Zoeppritz equations. Geophysics, 50(4): 609-614. DOI:10.1190/1.1441936
Simmons J L Jr, Backus M M. 1996. Waveform-based AVO inversion and AVO prediction-error. Geophysics, 61(6): 1575-1588. DOI:10.1190/1.1444077
Smith G C, Gidlow P M. 1987. Weighted stacking for rock property estimation and detection of gas. Geophysical Prospecting, 35(9): 993-1014. DOI:10.1111/gpr.1987.35.issue-9
Wang Y. 2003. Seismic Amplitude Inversion in Reflection Tomography. Amsterdam: Elsevier.
Wang Y H. 1999. Approximations to the Zoeppritz equations and their use in AVO analysis. Geophysics, 64(6): 1920-1927. DOI:10.1190/1.1444698
Wiggins R, Kenny G S, McClure C D. 1983. A method for determiningand displaying the shear-velocity reflectivities of a geologic formation. European Patent Application 0113944.
Yang P J, Yin X Y. 2008. Non-linear quadratic programming Bayesian prestack inversion. Chinese Journal of Geophysics (in Chinese), 51(6): 1876-1882.
Yin X Y, Zong Z Y, Wu G C. 2013. Improving seismic interpretation:a high-contrast approximation to the reflection coefficient of a plane longitudinal wave. Petroleum Science, 10(4): 466-476. DOI:10.1007/s12182-013-0297-y
Yin X Y, Zong Z Y, Wu G C. 2014. Seismic wave scattering inversion for fluid factor of heterogeneous media. Science China Earth Sciences, 57(3): 542-549. DOI:10.1007/s11430-013-4783-2
Yin X Y, Zhang S X. 2014. Bayesian inversion for effective pore-fluid bulk modulus based on fluid-matrix decoupled amplitude variation with offset approximation. Geophysics, 79(5): R221-R232. DOI:10.1190/geo2013-0372.1
Zoeppritz K. 1919. Erdbebenwellen Ⅷ B, Uber Reflexion und Durchgang seismischer Wellen durch Unstetigkeitsflachen. Gottinger Nachr, 1: 66-84.
Zong Z, Yin X. 2016. Direct inversion of Young's and Poisson impedances for fluid discrimination. Geofluids, 16(5): 1006-1016. DOI:10.1111/gfl.12202
Zong Z Y, Yin X Y, Wu G C. 2012a. AVO inversion and poroelasticity with P-and S-wave moduli. Geophysics, 77(6): N17-N24. DOI:10.1190/geo2011-0214.1
Zong Z Y, Yin X Y, Zhang F, et al. 2012. Reflection coefficient equation and pre-stack seismic inversion with Young's modulus and Poisson ratio. Chinese Journal of Geophysics (in Chinese), 55(11): 3786-3794. DOI:10.6038/j.issn.0001-5733.2012.11.025
Zong Z Y, Yin X Y, Wu G C. 2012b. Elastic impedance parameterization and inversion with Young's modulus and Poisson's ratio. Geophysics, 78(6): N35-N42.
Zong Z Y, Yin X Y, Wu G C. 2013. Direct inversion for a fluid factor and its application in heterogeneous reservoirs. Geophysical Prospecting, 61(5): 998-1005. DOI:10.1111/gpr.2013.61.issue-5
Zong Z Y, Yin X Y, Wu G C. 2015. Geofluid discrimination incorporating poroelasticity and seismic reflection inversion. Surveys in Geophysics, 36(5): 659-681. DOI:10.1007/s10712-015-9330-6
杨培杰, 印兴耀. 2008. 非线性二次规划贝叶斯叠前反演. 地球物理学报, 51(6): 1876-1882. DOI:10.3321/j.issn:0001-5733.2008.06.030
宗兆云, 印兴耀, 张峰, 等. 2012. 杨氏模量和泊松比反射系数近似方程及叠前地震反演. 地球物理学报, 55(11): 3786-3794. DOI:10.6038/j.issn.0001-5733.2012.11.025