0 引言

非均匀正交介质中地震波速度与波传播的方向相关，即地震各向异性现象(Thomsen，2002).多位学者的研究表明，具有垂直对称轴的横向各向同性背景介质(VTI，Transverse Isotropic with a Vertical symmetry axis)中发育一组平行排列的裂缝，或者纯各向同性背景介质中与VTI背景介质中发育两组垂直正交的裂缝，或者两组共轭裂缝均可等效为正交各向异性介质(Schoenberg and Helbig, 1997；Tsvankin，1997；Rüger，1996；Pšenčik and Martins, 2001；Bakulin et al., 2000, 2002；Shaw and Sen, 2004；Far et al., 2013；Liu and Martinez, 2012；Bachrach，2015).

由于非均匀正交介质中地震波传播的复杂性，很难推导其线性的纵波反射系数解析表达式.Rüger(1996)和Tsvankin(1997)引入了描述正交各向异性介质特征的各向异性参数表征记号.在此基础上，Pšenčik和Martins(2001)推导了正交各向异性介质的纵波反射系数方程.然而，其推导的方程受限于水平层状假设，而当地下存在裂缝等非均质目标时，地震波产生散射(印兴耀等，2013b)，此时水平层状假设不再满足非均匀各向异性介质参数的反演需求.弹性逆散射理论是非均匀介质参数反演的有效途径(Wu and Aki, 1985；Weglein，1993)，其常采用Born或Rytov近似(Stolt and Weglein, 1985；Mora，1987).基于弹性逆散射理论，Shaw和Sen(2004)借助稳相法推导了非均匀介质的纵波散射系数方程.裂缝弱度参数是储层裂缝描述的特征参数(Hsu and Schoenberg, 1993；陈怀震等，2014a)，有助于指导地下的裂缝识别，因此探索裂缝弱度参数与地震反射之间的关系尤为重要.Schoenberg和Helbig(1997)将OA介质的刚度矩阵建模为VTI背景介质矩阵与裂缝弱度参数之和.为此，我们使用弹性逆散射原理推导了以Thomsen各向异性参数(Thomsen，1986)与裂缝弱度参数表征的正交各向异性介质纵波散射系数方程.

由于分方位地震数据的低信噪比，正交各向异性介质参数的反演仍然是一个难题.由于角度叠加数据较高信噪比的优势，弹性阻抗(EI，Elastic Impedance)反演已经广泛应用于各向同性介质中(Connolly，1999；Whitcombe，2002；马劲风，2003；印兴耀等，2013b；张峰和李向阳，2015)，目前也扩展到各向异性介质(Martins，2006；陈天胜等，2006；Cui et al., 2010；陈怀震等，2014b).然而，其推导的EI方程的各向异性扰动部分存在指数校正项，与各向同性背景部分的表达式存在明显差异，可能增加各向异性参数反演的不稳定性，导致各向异性参数反演的精度相比于各向同性参数较低.为此，我们提出了一种新颖的正交各向异性参数化方法，并推导了相应的正交各向异性方位弹性阻抗方程.为了提高反演的稳定性与横向连续性，我们同时发展了贝叶斯框架下的正交各向异性方位弹性阻抗反演方法，同时考虑了柯西稀疏约束正则化和平滑模型约束正则化(Alemie and Sacchi, 2011；宗兆云等, 2012a, 2012b)，并最终使用非线性的迭代重加权最小二乘策略(Scales and Smith, 1994)实现了各向异性参数的稳定估算.

1 方法原理 1.1 非均匀正交介质有效弹性刚度矩阵

介质参数随空间位置变化，具有明显的非均匀特征.非均匀介质的模型参数化可使用一阶扰动理论，将非均匀介质的有效弹性刚度矩阵看为背景介质刚度矩阵C_b与扰动介质刚度矩阵ΔC的叠加，即

(1)

背景介质可为各向同性，也可为其他类型介质如各向异性、黏弹性等，同样可得到相应的非均匀介质刚度矩阵表达式.进一步，在不同的等效近似情况下，扰动介质刚度矩阵又可表示为

(2)

式中，ΔC_iso，ΔC_ani与ΔC_Q分别表示介质的各向同性、各向异性及黏滞性扰动刚度矩阵.

如图 1所示，水平细分层地层中发育一组垂直排列的裂缝时可等效为有效的正交各向异性介质，其某个裂缝对称面与裂缝面一致.水平层状地层可等效为VTI背景介质，其刚度矩阵C_VTI表示为

(3)

式中，C_ijb表示VTI背景介质的弹性刚度张量元素，且满足C_12b=C_11b－2C_66b.

在弱各向异性假设条件下，Thomsen(1986)引入了一组描述VTI介质各向异性特征的参数：

(4)

Hood(1991)研究了任意各向异性背景介质中发育任意方向的平行排列裂缝的柔度矩阵表征形式，若背景介质的固有各向异性及由裂缝引起的各向异性是较弱的，则等效介质的柔度矩阵认为是裂缝柔度矩阵与背景各向异性介质柔度矩阵的叠加.基于线性滑移理论(Schoenberg, 1980, 1983)，一组以x₁轴为法向的平行排列垂直裂缝的柔度矩阵s_f表征为(Schoenberg and Douma, 1988)

(5)

式中，三个正定的裂缝参数Z_N，Z_V与Z_H 分别是法向裂缝柔度、垂直切向裂缝柔度及水平切向裂缝柔度.因此，由VTI背景介质中发育一组平行垂直裂缝等效的正交各向异性介质的有效弹性刚度矩阵C_OA可表示为

(6)

裂缝弱度参数是储层裂缝描述的特征参数，有助于指导地下的裂缝识别.Hsu和Schoenberg(1993)引入了一组无量纲的描述裂缝的特征参数：

(7)

式中，M_b和μ_b是背景介质的纵波模量与横波模量；Δ_N，Δ_V与Δ_H分别表征法向裂缝弱度、垂直切向裂缝弱度及水平切向裂缝弱度.

因此，等效正交各向异性介质的刚度矩阵可用裂缝弱度参数表示为

(8)

基于一阶扰动理论，正交各向异性介质刚度矩阵C_OA可认为是VTI背景介质刚度矩阵C_VTIb与裂缝扰动介质刚度矩阵ΔC_HTI的叠加，即

(9)

式中，VTI背景介质的刚度矩阵又可认为是均匀各向同性背景介质刚度矩阵与各向同性扰动介质及VTI扰动介质刚度矩阵的叠加，即

(10)

因此，等效的正交各向介质刚度矩阵可表征为

(11)

也即，

(12)

式中，ΔM与Δμ分别表示界面两侧纵波模量与横波模量的差值；表示背景横、纵波速度的平方比；Δε，Δδ与Δγ分别表示界面两侧Thomsen各向异性参数的差值；δ_ΔN，δ_ΔV与δ_ΔH分别表示界面两侧裂缝弱度参数的差值.

1.2 非均匀正交介质纵波散射系数

在弱散射假设条件下，非均匀介质频率空间域弹性波散射波场表示为(Červený and Pšenčík，2001)

(13)

式中，ω表示圆频率；u_i⁰(·)与G_ni⁰(·)分别表示激发源和格林函数引起的弹性波场；Δρ=ρ－ρ_b表示非均匀介质与背景介质密度之差.弱反射的非均匀正交各向异性介质在均匀各向同性上的扰动刚度矩阵可表示为

(14)

考虑纵波传播条件，激发点源的格林函数可简单表示为

(15)

式中，r表示激发点到接收点的距离；N_n与N_i分别表示激发点和接收点的方向.

基于Born近似和稳相法，非均匀正交介质纵波散射系数可表示为(Shaw and Sen, 2004)：

(16)

式中，θ是入射角；r₀是水平界面一点，而S(r₀)表征散射函数(Eaton and Stewart, 1994；Burridge et al., 1998)，其表示为

(17)

其中，

(18)

p和t分别表示慢度矢量和极化矢量.

在此基础上，我们推导了正交各向异性介质的纵波散射系数方程(详见附录A)：

(19)

其中，

(20)

式中，ϕ=ϕ_obs－ϕ_sym是方位角，表示观测方位ϕ_obs与裂缝对称轴方位ϕ_sym的差值.为了获得线性的正交各向异性介质纵波散射系数方程，需首先计算裂缝对称轴方位ϕ_sym.在本文中，我们假设测线与裂缝对称轴方位平行，即ϕ_sym=0.R_M=ΔM/2M_b，R_μ=Δμ/2μ_b与R_ρ=Δρ/2ρ_b分别表示纵波模量、横波模量与密度的散射系数；R_ε=Δε/2与R_δ=Δδ/2分别表示Thomsen各向异性参数差值；R_ΔN=δ_ΔN/2，R_ΔV=δ_ΔV/2与R_ΔH=δ_ΔH/2分别表示裂缝弱度参数差值.

1.3 正交各向异性方位弹性阻抗参数化

Connolly(1999)定义弹性阻抗为

(21)

式中，n-1与n分别代表下层或上层介质；EI表示介质的弹性阻抗属性.对于方位各向异性介质而言，R_PP≡R_PP(θ, ϕ)且EI≡EI(θ, ϕ).

单界面的纵波散射系数方程(19)可扩展为时间连续的散射系数函数(Stolt and Weglein, 1985；Buland and Omre, 2003)：

(22)

对方程(22)进行积分计算，可得

(23)

对方程(23)取指数，则非均匀正交各向异性介质的弹性阻抗可表示为

(24)

为了消除弹性阻抗方程(24)中入射角变化对量纲尺度的影响，对其进行归一化(Whitcombe，2002)，可得

(25)

式中，EI₀=(M₀ρ₀)^1/2，M₀，μ₀与ρ₀分别表示介质的固有波阻抗及弹性参数属性；ε₀与δ₀表示介质的VTI特征属性；Δ_N0，Δ_V0与Δ_H0分别表示介质不同方向的裂缝弱度特征属性，其均可由测井信息或岩石物理实验信息获取.

观察方程(25)可以发现，推导的方位各向异性EI方程的各向异性扰动指数校正项与各向同性背景项的表达式存在明显差异，增加了各向异性参数反演的不稳定性，导致各向异性参数反演的精度相比于各向同性参数较低.为了消除各向异性扰动项与各向同性背景项的表达差异，我们引入了一组无量纲的正交各向异性特征参数(详见附录B)：

(26)

因此，非均匀正交各向异性的纵波散射系数方程(19)可表示为

(27)

式中，与表示Thomsen各向异性特征参数的散射系数；，与表示为裂缝弱度特征参数的散射系数.

同样地，按照上述弹性阻抗方程的推导步骤，可得到新的非均匀正交各向异性方位弹性阻抗表达式为

(28)

1.4 正交各向异性方位弹性阻抗反演

非均匀正交各向异性介质的方位弹性阻抗反演主要包括方位各向异性弹性阻抗数据体的反演及介质弹性参数与各向异性特征参数的提取.已知不同方位角度的部分角度叠加地震数据体，从中反演得到不同方位角的方位各向异性弹性阻抗数据体，以测井数据和岩石物理预测结果为约束，最终反演获得正交各向异性参数.然而，正交各向异性参数的反演到目前为止仍然是一个难题.

为了提高参数反演的稳定性与抗噪性，我们发展了一种贝叶斯框架下的正交各向异性方位弹性阻抗稳定反演方法，同时考虑柯西稀疏约束正则化和平滑模型约束正则化，最终使用非线性的迭代重加权最小二乘策略实现了各向异性参数的稳定估算.

对归一化的正交各向异性方位弹性阻抗方程(28)两边取对数，可得到用于正交各向异性介质弹性参数与各向异性参数估算的线性化方位各向异性弹性阻抗表达式：

(29)

式(29)写成矩阵形式为

(30)

其中，

(31)

M与N分别为入射角与方位角的个数；L_EI(t, θ₁, ϕ₁)=ln[EI(t, θ, ϕ)/EI₀]表示为对数域方位各向异性弹性阻抗散射系数，类似地，m表示为对数域弹性参数与各向异性特征参数散射系数.

经验和实际观测表明，介质的弹性参数、各向异性参数间存在很强的相关性，大大影响参数反演的合理性(Downton，2005).参数的协方差矩阵可用于待反演参数的去相关，其可统计计算或由岩石物理的关系获取.正交各向异性介质的协方差矩阵C_m可表示为

(32)

式中，协方差矩阵的对角线元素表示对数域弹性参数与各向异性特征参数的方差，而非对角线元素表示参数间的协方差.对协方差矩阵进行奇异值分解(SVD，Singular Value Decomposition)可得：

(33)

式中，V表示介质参数协方差矩阵的特征矢量，Σ表示协方差矩阵的特征值.

对式(30)中的参数进行去相关，得

(34)

式中，G′=GV与m′=V^－1m分别表示去相关后的系数矩阵及对数域参数的散射系数.

假定地震噪声服从高斯分布，而待反演的参数服从柯西分布以突出地层的弱散射，则后验概率密度函数可表示为

(35)

式中，σ_n²与σ_m²分别表示噪声方差及待反演参数的方差；K表示参数的个数.对方程(35)取对数，则目标函数可写为

(36)

式中，λ_Cauchy表示柯西稀疏约束正则化因子.

由于地震数据的低频分量相对较弱，易受噪声影响.平滑模型约束正则化能够补充反演参数的低频信息，提高其稳定性及横向连续性.对方程(36)添加平滑模型约束正则化项，得

(37)

式中，

(38)

λ_i(i=M, μ, ρ, ε^r, δ^r, Δ_N^r, Δ_V^r, Δ_H^r)为弹性参数与各向异性特征参数平滑模型约束的正则化参数；P表示积分算子矩阵.

柯西分布的引入使得方程(37)具有弱非线性，采用重加权最小二乘反演方法(IRLS，Iteratively Reweighted Least Squares)能够稳定的获得参数的反演结果.

2 模型试算与实际资料处理 2.1 模型试算

为了使用常规测井资料计算正交各向异性介质的弹性参数及各向异性特征参数，从而建立参数的背景平滑模型，本文根据Xu和Payne(2009)及张广智(2013, 2015)等人的研究，初步建立了适用的正交各向异性岩石物理等效模型，主要步骤如下：

(1) VTI背景介质等效弹性模量的获取：使用Voigt-Reuss-Hill模型(Mavko et al., 2009)计算混合矿物的弹性模量，然后使用Backus平均(Backus，1962)混合VIT背景介质，计算最终的VTI背景介质等效弹性模量；

(2) 考虑孔隙的“干”骨架等效弹性模量的求取：利用各向异性SCA-DEM模型(Hornby et al., 1994)向VTI背景介质中加入孔隙系统，得到考虑孔隙的“干”骨架弹性模量；

(3) 考虑裂缝的“干”骨架等效弹性模量的求取：在弱各向异性假设的条件下(Thomsen，1986)，利用Hudson模型(Hudson, 1980, 1981)与Schoenberg线性滑移模型(Schoenberg, 1980, 1983)向VTI背景介质中加入裂缝系统，构建正交各向异性等效介质模型，得到考虑裂缝的“干”骨架弹性模量；

(4) 饱和岩石等效弹性模型的计算：使用Brown-Korringa方程(Brown and Korringa, 1975)将流体加入构建的考虑孔隙与裂缝的“干”骨架中，计算饱和岩石弹性模量，进而计算正交各向异性等效介质的弹性参数及各向异性特征参数.

基于建立的正交各向异性岩石物理等效模型，本文首先使用测井解释数据估测了正交各向异性介质的弹性参数及各向异性特征参数，如图 2所示.为了验证非均匀正交介质各向异性参数地震散射反演方法的可行性与稳定性，本文采用实际的单井模型进行验证.模型测试使用主频35 Hz的Ricker子波.图 3是使用精确Zoeppritz方程进行正演合成的信噪比(SNR，Signal-to-Noise-Ratio)分别为4和2的角度域叠前方位道集.图 4与图 5分别是SNR为4和2的方位叠前道集反演获得的参数结果.反演结果表明，在适当噪声的情况下，弹性参数与各向异性特征参数的反演结果与真实值仍有较高的吻合度.从反演精度而言，由于密度项对纵波散射系数的贡献较小，密度参数相较纵波模量与横波模量反演精度较低.同时，由于方位叠前信息的使用，裂缝弱度特征参数相较Thomsen各向异性特征参数反演精度高.

2.2 实际资料处理

实际方位地震道集来自于中国路上某勘探工区，进行叠前反演前，需对地震数据进行保幅等预处理.笔者利用本文提出的方法对不同角度叠加剖面进行了实际资料的试处理.本文反演方法输入的是不同方位的部分角度叠加道集，图 6(a)、(b)及(c)分别是四个方位的小、中、大角度部分叠加剖面.图 7(a)、(b)及(c)分别是基于CSSI反演的不同方位的弹性阻抗剖面，图 8、图 9及图 10分别是反演获得的纵波模量、横波模量及密度弹性参数、两个Thomsen各向异性特征参数及三个裂缝弱度特征参数剖面，其反演结果与实际钻遇结果一致(红色虚线椭圆为裂缝型气藏发育区)，验证了方法的有效性.

3 结论

正交各向异性特征参数的可靠估算有助于非均匀正交各向异性介质的特征描述，利用方位各向异性弹性阻抗反演获取裂缝弱度特征参数成为利用方位叠前资料进行裂缝发育区识别的重要手段.地震散射理论是非均匀介质参数反演的有效途径，本文从弹性逆散射原理出发，首先推导了非均匀正交介质中纵波散射系数方程，并通过正交各向异性特征参数的引入，提出了一种新颖的正交各向异性方位弹性阻抗参数化方法，奠定了非均匀正交各向异性介质方位弹性阻抗反演的理论基础.在此基础上，发展了贝叶斯框架下柯西稀疏正则化及平滑模型正则化约束的正交各向异性方位弹性阻抗反演方法，实现了正交各向异性特征参数的稳定估算.模型和实际资料处理表明，该方法可从分方位叠前地震资料中获取正交各向异性特征参数，为非均匀正交介质的各向异性预测提供了一种高可靠性的地震反演方法.

然而，单纯使用方位纵波数据预测正交各向异性特征具有一定的不确定性.在实际应用中，需充分利用地质、地震、测井及岩心资料进行综合考虑，在方位振幅信息充分利用的前提下，综合方位速度及旅行时等信息，提高正交各向异性预测的精度.

附录A 非均匀正交各向异性介质纵波散射系数方程推导

方程(18)中慢度矢量p和极化矢量t的下标m，n，i，j，k和l可表示为

(A1)

式中，δ_ij与δ_kl均表示Kronecker符号.

纵波慢度矢量p表示为

(A2)

纵波极化矢量t表示为

(A3)

因此，

(A4)

(A5)

将式(A4)与(A5)代入式(17)中，则纵波散射函数可表示为

(A6)

故非均匀正交各向异性介质的纵波散射系数方程可最终表示为

(A7)

当介质为均匀各向同性介质时，方程(A7)退化为各向同性介质的纵波散射系数方程(宗兆云等，2012 ).同样地，当介质分别为VTI介质或HTI介质时，方程(A7)退化为VTI介质或HTI介质的纵波散射系数方程(Rüger，1996 ；Shaw and Sen, 2006).

附录B 无量纲正交各向异性特征参数的引入

假定R_ε^r=R_ε，即，也即

(B1)

对方程(B1)左右两侧进行积分，得

(B2)

也即，

(B3)

式中，C是一常数.各向异性特征参数ε^r可写为ε的一阶扰动形式，即

(B4)

式中，m与n均为两个未知量.将方程(B4)代入方程(B3)可得：

(B5)

同样，将方程(B5)代入方程(B4)可得：

(B6)

取C=1，则

(B7)

同样地，对Thomsen各向异性参数δ及裂缝弱度参数Δ_N，Δ_V与Δ_H作类似的变形处理，则

(B8)