页岩气是蕴藏于页岩地层中可供开采的天然气资源, 中国的页岩气可采储量较大, 具有良好的勘探、开发潜力。随着油气勘探开发逐渐由常规油气向非常规油气尤其是页岩气发展, 以及逐步走向精细化, 储层精细描述至关重要。地震反演技术可揭示地下油气藏的精细分布特征, 因此受到广泛关注^[1]。页岩气储层的水平薄互层构造和水平层理具有明显的VTI各向异性特征^[2], 大量研究表明, 页岩气储层的水平层理缝引起的各向异性影响地震反射特征^[3-7], 裂缝弱度参数作为非常规页岩气储层水平层理缝表征的关键参数^[8], 有助于页岩气储层各向异性特征的识别与描述, 因此十分有必要开展基于VTI等效介质的叠前多参数反演方法研究。

建立VTI等效介质弹性参数及各向异性参数随地震振幅的变化关系是地震反演的基础。针对VTI等效介质, 众多学者推导了反射系数和透射系数的精确表达式^[9-11], 由于方程形式复杂, 限制了VTI等效介质反射系数精确方程的实际应用。为此, 众多学者针对不同类型的储层特征参数推导出相应的VTI等效介质反射系数近方程^[12-15]。THOMSEN^[12]基于弱各向异性假设, 推导了VTI等效介质的P波反射系数线性化近似方程。URSIN等^[13]基于弱阻抗差和弱各向异性假设, 推导了VTI等效介质的反射系数和透射系数线性化近似方程。RÜGER^[14]进一步完善了THOMSEN^[12]提出的P波反射系数近似方程, 以适应大入射角度的情况, 该公式即为经典的VTI等效介质反射系数方程。对于具有任意对称性的弱各向异性介质, SHAW等^[16]基于稳相法推导出了PP波和PS波反射系数方程。本文基于波恩积分法, 推导了VTI等效介质纵波反射系数近似方程, 构建了储层弹性模量、裂缝弱度参数与地震反射系数之间的量化关系, 为后续叠前反演奠定基础。

稳定、可靠的地震反演是页岩气储层精细表征的关键。贝叶斯反演理论常被用于求解多元目标函数^[17-27]。从目前的研究趋势看, 地震反演方法研究主要包括叠前AVO反演和全波形反演。基于精确Zoeppritz方程的地震反演方法能够利用丰富的叠前地震信息预测模型参数, 但其非线性强、稳定性差, 全波形反演方法则计算量大、耗时长, 反演尺度和计算效率不能满足储层精细表征要求。针对VTI等效介质理论模型, 传统的近似反射系数方程为包含弹性参数和各向异性参数的五项模型参数表达式, 由于反演参数较多, 造成的反演病态问题限制了方法的应用^[28-33], 且未考虑裂缝诱导的各向异性效应。

针对传统VTI等效介质的叠前多参数反演存在的问题, 本文开展了面向页岩气储层的叠前多参数反演方法研究, 提出了一种稳健的基于VTI等效介质AVA反演的储层弹性模量与裂缝弱度参数评价方法。首先, 基于散射理论以及弱各向异性近似假设条件, 推导出含纵、横波模量、密度和裂缝弱度参数的VTI等效介质纵波反射系数近似方程, 经过整合与简化, 获得含4个组合属性参数(A、B、C、D属性)的线性化纵波反射系数近似方程。进一步, 基于贝叶斯反演理论提出了VTI等效介质四项组合属性参数反演方法, 实现了组合属性参数的稳定、可靠预测。然后, 通过建立组合属性参数与储层参数之间的内在关系, 获得弹性模量和裂缝弱度等储层参数, 实现页岩气储层弹性特征与裂缝发育特征的地震准确描述。最后, 利用VTI等效介质反射系数方程合成AVA叠前道集和实际数据, 验证所提出方法的正确性和有效性。

1 方法理论

针对页岩气储层精细表征与描述问题, 开展基于VTI等效介质AVA反演的储层弹性模量与裂缝弱度参数反演方法研究。首先基于推导的五项式VTI介质纵波反射系数近似方程, 经过四项组合属性参数表征的新VTI介质纵波反射系数方程, 以减少反演的模型参数, 提高反演的稳定性; 进一步结合贝叶斯反演理论以及属性之间的数学关系, 实现储层弹性模量和裂缝参数的稳定、可靠反演, 有效刻画页岩气储层的弹性特征和裂缝发育特征。

1.1 含裂缝弱度的VTI等效介质P波反射系数线性化推导

VTI对称介质存在5个独立的弹性刚度分量, 其6×6的弹性刚度矩阵C_VTI可以表示为:

$ \begin{aligned} & \boldsymbol{C}_{\mathrm{VTI}}= \\ & {\left[\begin{array}{cccccc} C_{11} & C_{11}-2 C_{66} & C_{13} & 0 & 0 & 0 \\ C_{11}-2 C_{66} & C_{11} & C_{13} & 0 & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{array}\right]} \end{aligned} $

(1)

式中: C₁₁, C₁₃, C₃₃, C₄₄, C₆₆分别为弹性刚度分量。

在VTI等效介质的弱各向异性假设下, Thomsen简化了上述复杂关系式, 定义了5个弹性常数:

$ \left\{\begin{array}{l} v_{\mathrm{P}}=\sqrt{\frac{C_{33}}{\rho}} \quad v_{\mathrm{S}}=\sqrt{\frac{C_{44}}{\rho}} \\ \varepsilon=\frac{C_{11}-C_{33}}{2 C_{33}} \quad \gamma=\frac{C_{66}-C_{44}}{2 C_{44}} \\ \delta=\frac{\left(C_{13}+C_{44}\right)^2-\left(C_{33}-C_{44}\right)^2}{2 C_{33}\left(C_{33}-C_{44}\right)} \end{array}\right. $

(2)

式中: v_P, v_S和ρ分别为各向同性背景部分的纵、横波速度和密度; ε, δ, γ为弱各向异性参数, 也称为Thomsen各向异性参数。根据弱各向异性近似假设条件, RÜGER^[14]推导了VTI等效介质纵波反射系数近似方程, 具体形式如下:

$ \begin{gathered} R_{\mathrm{PP}}^{\mathrm{VTI}}(\theta)=\frac{1}{2}\left(\frac{\Delta \rho}{\bar{\rho}}+\frac{\Delta v_{\mathrm{P}}}{\bar{v}_{\mathrm{P}}}\right)+\frac{1}{2}\left\{\frac{\Delta v_{\mathrm{P}}}{\bar{v}_{\mathrm{P}}}-\left(2 \frac{\bar{v}_{\mathrm{S}}}{\bar{v}_{\mathrm{P}}}\right)^2 \cdot\right. \\ \left.\left(\frac{\Delta \rho}{\bar{\rho}}+\frac{2 \Delta v_{\mathrm{S}}}{\bar{v}_{\mathrm{S}}}\right)\right\} \sin ^2 \theta+\frac{1}{2} \frac{\Delta v_{\mathrm{P}}}{\bar{v}_{\mathrm{P}}} \sin ^2 \theta \tan ^2 \theta+ \\ \frac{1}{2} \sin ^2 \theta \Delta \delta+\frac{1}{2} \sin ^2 \theta \tan ^2 \theta \Delta \varepsilon \end{gathered} $

(3)

式中: θ为入射角; Δε, Δδ分别为上下两层介质的Thomsen弱各向异性参数差值。

根据线性滑移理论, 裂缝诱导的VTI等效介质刚度系数可以表示为^[8]:

$ \boldsymbol{C}_{\mathrm{VTI}}^{\prime}=\left[\begin{array}{cccccc} M\left(1-\xi^2 \delta_{\mathrm{N}}\right) & \lambda\left(1-\xi \delta_{\mathrm{N}}\right) & \lambda\left(1-\delta_{\mathrm{N}}\right) & 0 & 0 & 0 \\ \lambda\left(1-\xi \delta_{\mathrm{N}}\right) & M\left(1-\xi^2 \delta_{\mathrm{N}}\right) & \lambda\left(1-\delta_{\mathrm{N}}\right) & 0 & 0 & 0 \\ \lambda\left(1-\delta_{\mathrm{N}}\right) & \lambda\left(1-\delta_{\mathrm{N}}\right) & M\left(1-\delta_{\mathrm{N}}\right) & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu\left(1-\delta_{\mathrm{T}}\right) & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu\left(1-\delta_{\mathrm{T}}\right) & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{array}\right] $

(4)

式中: M和μ分别表示纵波、横波弹性模量, M=λ+2μ, ξ=λ/M; δ_N和δ_T分别表示法向和切向裂缝弱度参数。

基于散射理论和玻恩近似, 在弱各向异性假设条件下, 推导出基于VTI等效介质的纵波反射系数五项式近似方程(具体推导过程见附录A):

$ \begin{gathered} R_{\mathrm{PP}}^{\mathrm{VTI}}(\theta)=\frac{1}{4} \sec ^2 \theta \frac{\Delta M}{M}-2 g \sin ^2 \theta \frac{\Delta \mu}{\mu}+\frac{\cos 2 \theta}{4 \cos ^2 \theta} \cdot \\ \frac{\Delta \rho}{\rho}+\frac{1}{4 \cos ^2 \theta} \Delta \delta_{\mathrm{N}}\left(-1-4 g^2 \sin ^4 \theta+4 g \sin ^2 \theta\right)+ \\ g \sin ^2 \theta \Delta \delta_{\mathrm{T}} \end{gathered} $

(5)

式中: ΔM/M和Δμ/μ是纵、横波弹性模量的反射率; Δρ/ρ是密度反射率; g是纵横波速度比。

为了减少待反演模型参数的数量, 将公式(5)进一步改写为如下形式(具体简化过程见附录B):

$ \begin{gathered} R_{\mathrm{PP}}^{\mathrm{VTI}}(\theta)=\frac{1}{4}\left(\frac{\Delta M}{M}+\frac{\Delta \rho}{\rho}-\Delta \delta_{\mathrm{N}}\right)-2 g \sin ^2 \theta\left(\frac{\Delta \mu}{\mu}-\right. \\ \left.\frac{1}{2} g \Delta \delta_{\mathrm{N}}-\frac{1}{2} \Delta \delta_{\mathrm{T}}\right)+\frac{1}{4} \tan ^2 \theta\left(\frac{\Delta M}{M}-\frac{\Delta \rho}{\rho}-\Delta \delta_{\mathrm{N}}\right)- \\ g(g-1) \sin ^2 \theta \tan ^2 \theta \Delta \delta_{\mathrm{N}} \end{gathered} $

(6)

在界面两侧弹性特征差异较小的情况下, 对公式(6)进行BORTFELD^[34]和OLDENBURG等^[35]提出的近似, 也即:

$ \left\{\begin{array}{l} \frac{\Delta M_d}{M_d} \approx \ln \left[\frac{(M)_{i+1}}{(M)_i}\right] \\ \frac{\Delta \mu}{\mu} \approx \ln \left(\frac{\mu_{i+1}}{\mu_i}\right) \\ \frac{\Delta \rho}{\rho} \approx \ln \left(\frac{\rho_{i+1}}{\rho_i}\right) \end{array}\right. $

(7)

式中i+1和i分别表示第i+1层和第i层。将公式(7)代入公式(6), 得到如下表达式:

$ R_{\mathrm{PP}}^{\mathrm{VTI}}(\theta)=\frac{1}{4} \ln \left[\frac{\left(M \rho \mathrm{e}^{-\delta_{\mathrm{N}}}\right)_{i+1}}{\left(M \rho \mathrm{e}^{-\delta_{\mathrm{N}}}\right)_i}\right]-2 g \sin ^2 \theta \cdot\\ \ln \left[\frac{\left(\mu \mathrm{e}^{-\frac{1}{2} g \Delta \delta_{\mathrm{N}}-\frac{1}{2} \Delta \delta_{\mathrm{T}}}\right)_{i+1}}{\left(\mu \mathrm{e}^{-\frac{1}{2} \mathrm{g} \Delta \delta_{\mathrm{N}}-\frac{1}{2} \Delta \delta_{\mathrm{T}}}\right)_i}\right]+\frac{1}{4} \tan ^2 \theta \ln \left[\frac{\left(\frac{M \mathrm{e}^{-\delta_{\mathrm{N}}}}{\rho}\right)_{i+1}}{\left(\frac{M \mathrm{e}^{-\delta_{\mathrm{N}}}}{\rho}\right)_i}\right]-\\ \;\;\left.g(g-1) \sin ^2 \theta \tan ^2 \theta \ln \left[\frac{\left(\mathrm{e}^{\delta_{\mathrm{N}}}\right)_{i+1}}{\left(\mathrm{e}^\delta \mathrm{N}\right.}\right)_i\right] $

(8)

设:

$ \left\{\begin{array}{l} A=M \rho \mathrm{e}^{-\delta_{\mathrm{N}}} \\ B=\mu \mathrm{e}^{-\frac{1}{2} g \Delta \delta_{\mathrm{N}}-\frac{1}{2} \Delta \delta_{\mathrm{T}}} \\ C=\frac{M \mathrm{e}^{-\delta_{\mathrm{N}}}}{\rho} \\ D=\mathrm{e}^{\delta_{\mathrm{N}}} \end{array}\right. $

(9)

则有以下表达式:

$ \begin{gathered} R_{\mathrm{PP}}^{\mathrm{VTI}}(\theta)=\frac{1}{4} \ln \left[\frac{A_{i+1}}{A_i}\right]-2 g \sin ^2 \theta \ln \left[\frac{B_{i+1}}{B_i}\right]+\frac{1}{4} \tan ^2 \cdot \\ \theta \ln \left[\frac{C_{i+1}}{C_i}\right]-g(g-1) \sin ^2 \theta \tan ^2 \theta \ln \left[\frac{D_{i+1}}{D_i}\right] \end{gathered} $

(10)

公式(10)被称为“含四项组合属性参数的新方程”, 新属性A, B, C和D是用于VTI等效介质的弹性模量和裂缝弱度预测的四项模型组合参数, A=Mρe^－δ_N, B=μe^{－1/2gΔδ_N－1/2Δδ_T}, C=(Me^－δ_N)/ρ, D=e^δ_N, A为纵波模量、密度与裂缝弱度参数相组合的属性, B为剪切模量与裂缝弱度参数相组合的属性, C为剪切模量与裂缝弱度参数、密度相组合的属性, D为指数的裂缝弱度参数。4个属性的物理意义可以解释为: A和C属性代表弹性性质和裂缝弱度的项, 两者相除与密度相关, 两者的乘积与纵波模量和裂缝弱度有关; 在裂缝弱度较小的情况下, 属性B近似等于剪切模量; D属性与裂缝弱度直接相关。

设计两层反射系数模型(表 1), 对比近似公式与精确公式的纵波反射系数。反射系数误差计算公式为: error=|R_精确－R_近似|/R_精确, 由图 1可算出, 在入射角达到30°时, 反射系数计算误差为3.26%, 可见近似公式的精度完全满足反演的要求。

1.2 基于贝叶斯反演理论的VTI等效介质新属性参数估计

根据褶积理论, 地震振幅响应可表示为如下矩阵形式:

$ \boldsymbol{d}=\boldsymbol{G} \boldsymbol{m} $

(11)

式中: d是观测数据; 模型向量m =[A, B, C, D]^T, T表示转置; 正演算子G = WP, 其中, W是子波矩阵; P是敏感度矩阵, 具体表现形式如下:

$ \begin{aligned} & \boldsymbol{P}= \\ & {\left[\begin{array}{cccc} \frac{1}{4} & -2 g \sin ^2 \theta_1 & \frac{1}{4} \tan ^2 \theta_1 & -g(g-1) \tan ^2 \theta_1 \\ \vdots & \vdots & \vdots & \vdots \\ \frac{1}{4} & -2 g \sin ^2 \theta_m & \frac{1}{4} \tan ^2 \theta_m & -g(g-1) \tan ^2 \theta_m \end{array}\right]} \end{aligned} $

(12)

式中: θ₁, θ₂, ⋯, θ_m表示m个入射角。基于贝叶斯定理, 密度p(m|d)与先验概率密度p(m)和p(d|m)的乘积成正比^[17]:

$ p(\boldsymbol{m} \mid \boldsymbol{d}) \propto p(\boldsymbol{m}) p(\boldsymbol{d} \mid \boldsymbol{m}) $

(13)

当观测地震资料包含高斯随机噪声时, 似然函数p(d|m)可表示为:

$ \begin{aligned} & p(\boldsymbol{d} \mid \boldsymbol{m})=\frac{1}{\left(2 \pi \sigma_e^2\right)^{\frac{N}{2}}} \\ & \exp \left\{-\sum \frac{[\boldsymbol{d}-\boldsymbol{G}(\boldsymbol{m})]^{\mathrm{T}}[\boldsymbol{d}-\boldsymbol{G}(\boldsymbol{m})]}{2 \sigma_e^2}\right\} \end{aligned} $

(14)

式中: σ_e²是噪声的方差; N是输入数据采样点数; G(m)为基于模型参数的正演合成记录。假设模型参数之间相互独立, 先验概率p(m)的具体表现形式为:

$ p(\boldsymbol{m})=\frac{1}{(2 \pi)^{\frac{N}{2}}\left|C_\boldsymbol{m}\right|^{1 / 2}} \exp \left(-1 / 2 \boldsymbol{m}^{\mathrm{T}} C_\boldsymbol{m}^{-1} \boldsymbol{m}\right) $

(15)

式中: C_m为模型参数的方差。根据贝叶斯理论, 求解反演问题可以转化为求解观测数据与拟合数据之间的残差最小值, 即:

$ \begin{gathered} J(\boldsymbol{m})=\frac{1}{2} \exp \left\{-\frac{1}{2 \delta_d^2}[\boldsymbol{d}-\boldsymbol{G}(\boldsymbol{m})]^{\mathrm{T}}[\boldsymbol{d}-\right. \\ \boldsymbol{G}(\boldsymbol{m})]+\boldsymbol{R}(\boldsymbol{m})\} \end{gathered} $

(16)

$ \boldsymbol{R}(\boldsymbol{m})=\frac{1}{2} \boldsymbol{m}^{\mathrm{T}} C_ \boldsymbol{m}^{-1} \boldsymbol{m} $

(17)

公式(16)和公式(17)中, J(m)是待求解的目标泛函, R(m)是模型参数正则化项, σ_d²是数据方差。根据最优化原理, 求解能量残差极小值等价于能量残差对未知模型参数的导数为零, 即:

$ \frac{\partial J(\boldsymbol{m})}{\partial \boldsymbol{m}}=\frac{1}{\sigma_d^2} \boldsymbol{L}^{\mathrm{T}}(\boldsymbol{L m}-\boldsymbol{d})+\frac{\partial R(\boldsymbol{m})}{\partial \boldsymbol{m}}=0 $

(18)

因此, m的解可表示为^{[23, 36-37]}:

$ \boldsymbol{m}=\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{G}+\frac{1}{2} \mu_h \boldsymbol{C}_\boldsymbol{m}^{-1}\right)^{-1}\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{d}+\frac{1}{2} \mu_h C_\boldsymbol{m}^{-1}\right) $

(19)

式中: μ_h是模型参数扰动的权重。

2 实际应用 2.1 合成数据算例

对实际测井曲线进行预处理后, 首先根据构建的页岩储层岩石物理模型, 以矿物组分、孔隙度、含水饱和度等储层参数为输入, 通过岩石物理正演, 获得弹性参数和裂缝弱度参数(图 2), 通过计算进一步获得待估新属性参数(图 3)。

利用雷克子波(主频为30 Hz)与由公式(10)计算的反射系数进行褶积得到合成叠前角度道集(图 4)。为了验证反演方法在含噪情况下的可靠性和稳定性, 在无噪合成地震数据上添加信噪比分别为5和2的高斯随机噪声, 采用频谱估算法计算随机噪声, 如图 4a和图 4b所示, 叠前道集的入射角范围为3°~30°, 每隔3°为一道, 共有10道。反演结果如图 5所示, 从图 5中可看出, 在信噪比分别为5和2的情况下, 反演值(黑色)与真实值(红色)之间均具有良好的一致性。因此, 验证了本文所提出的反演方法的正确性和有效性。

2.2 实际数据反演

过井地震剖面如图 6所示, X井位于CDP为2 122处, 在采样时刻2 200~2 400 ms位置钻遇含气储层。基于实际工区构建的页岩岩石物理模型, 通过岩石物理正演获得表征储层弹性特征和裂缝特征的参数, 测井曲线如图 7所示, 在采样时刻2 235 ms处, 速度和密度的值相对较低, 裂缝弱度的值相对较高, 结合测井解释成果, 明确该处为含气储层。利用上述测井曲线进一步计算待反演的组合属性参数, 为后续反演奠定模型基础。

图 8为基于贝叶斯反演理论获得的四项组合属性参数的反演剖面。叠合测井曲线显示, 可看出反演结果与测井曲线吻合度高、一致性好。为了进一步对比反演方法的效果, 抽取过井处的单道反演结果, 开展井震对比与分析(图 9), 可看出在目的层(采样时刻为2 235 ms处)反演的组合属性参数(红线)与测井真实值(黑线)吻合度高。为了进一步验证反演的合理性和稳定性, 对比CDP为2 122处的原始道集、经预处理后的叠前道集和利用估计的新属性参数合成道集(图 10), 可看出三者具有很好的一致性。

在获得四项组合属性参数反演剖面的基础上, 首先预测纵波模量和剪切模量, 通过对比反演属性B(红线)、实际剪切模量(黑线)和实际测井资料计算获得的属性B(蓝线)(图 11a), 可看出属性B(蓝线)和真实剪切模量(黑线)之间差异较小, 因此, 反演属性B可近似表征各向同性背景介质的剪切模量。进一步利用反演属性A, C和D预测纵波模量, 计算A与C的乘积为A·C=M²e^－2δ_N, 代入由属性D计算获得的δ_N, 则可计算得到纵波模量, 对比反演(红线)和实际(黑线)的纵波模量(图 11b)可见, 两者吻合度高。进一步根据反演属性A和属性C, 计算过渡属性(A/C=ρ²), 开方即可获得密度预测结果, 对比真实值(黑线)与反演值(红线)(图 11c)可见, 两者趋势一致。利用反演获得的D属性预测裂缝弱度, 对比反演(红线)和实际(黑线)裂缝弱度(图 11d)也发现, 两者一致性较好。

将上述计算单井的方法应用于地震数据, 利用反演获得的属性A和属性C(图 8a和图 8c)计算纵波模量和密度(图 12a和图 12b); 利用反演获得的属性D(图 8d)计算裂缝弱度(图 12c)。从反演剖面可看出, 目的层高产储层的纵波模量和密度表现为低值异常, 裂缝弱度表现为高值异常, 与测井解释结论一致, 同时X井处的预测结果与测井曲线具有良好的一致性, 充分验证了该方法的合理性和有效性。

3 结论

针对页岩气储层精细表征与描述问题, 提出了一种稳定的基于VTI等效介质AVA反演的储层弹性模量与裂缝弱度参数评价方法, 并利用合成数据和实际资料进行了测试与应用, 得出如下结论和认识。

1) 基于五项式VTI等效介质纵波反射系数近似方程, 构建了四项组合属性参数表征的新VTI等效介质纵波反射系数近似方程, 有效减少了反演的模型参数, 为VTI等效介质叠前AVA反演奠定基础。

2) 基于贝叶斯反演理论及属性之间的数学关系, 实现储层弹性模量和裂缝弱度参数的稳定、可靠反演, 有效刻画页岩气储层的弹性特征和裂缝发育特征。

3) 合成数据和实际资料的测试与应用显示了良好的效果, 因此本文提出的方法在页岩气的勘探开发中具有广泛的潜在应用前景。

需要指出的是, 本文的方法基于VTI等效介质理论假设, 实际储层中常发育高角度裂缝, 因此仍需进一步研究基于正交各向异性介质的叠前多参数反演方法。

附录A 含裂缝弱度参数的线性化纵波反射系数近似方程的推导

根据裂缝诱导的VTI介质线性滑移理论, 其刚度系数可以表示为:

$ \left\{\begin{array}{l} C_{11}=M-\xi^2 M \delta_{\mathrm{N}} \\ C_{12}=\lambda-\xi \lambda \delta_{\mathrm{N}} \\ C_{13}=\lambda-\lambda \delta_{\mathrm{N}} \\ C_{33}=M-M \delta_{\mathrm{N}} \\ C_{44}=\mu-\mu \delta_{\mathrm{T}} \\ C_{66}=\mu \end{array}\right. $

(A1)

在此基础上, 进一步求取刚度参数的扰动矩阵:

$ \left\{\begin{array}{l} \Delta C_{11}=\Delta M-\xi^2 M \Delta \delta_{\mathrm{N}}-\xi^2 \Delta M \delta_{\mathrm{N}} \\ \Delta C_{12}=\Delta \lambda-\xi \lambda \Delta \delta_{\mathrm{N}}-\xi \Delta \lambda \delta_{\mathrm{N}} \\ \Delta C_{13}=\Delta \lambda-\lambda \Delta \delta_{\mathrm{N}}-\Delta \lambda \delta_{\mathrm{N}} \\ \Delta C_{33}=\Delta M-M \Delta \delta_{\mathrm{N}}-\Delta M \delta_{\mathrm{N}} \\ \Delta C_{44}=\Delta \mu-\mu \Delta \delta_{\mathrm{T}}-\Delta \mu \delta_{\mathrm{T}} \\ \Delta C_{66}=\Delta \mu \end{array}\right. $

(A2)

式中: ΔM, Δμ, Δδ_N和Δδ_T分别表示纵波模量、横波模量、法向裂缝弱度参数和切向裂缝弱度参数在两个VTI介质分界面上的变化量。

假设裂缝弱度参数较小并且界面两侧纵波和横波模量差较小的情况下, 我们忽略了与ΔMδ_N, Δλδ_N和Δμδ_T成比例的项, 则刚度参数的扰动简化为:

$ \left\{\begin{array}{l} \Delta C_{11} \approx \Delta M-\xi^2 M \Delta \delta_{\mathrm{N}} \\ \Delta C_{12} \approx \Delta \lambda-\xi \lambda \Delta \delta_{\mathrm{N}} \\ \Delta C_{13} \approx \Delta \lambda-\lambda \Delta \delta_{\mathrm{N}} \\ \Delta C_{33} \approx \Delta M-M \Delta \delta_{\mathrm{N}} \\ \Delta C_{44} \approx \Delta \mu-\mu \Delta \delta_{\mathrm{T}} \\ \Delta C_{66} \approx \Delta \mu \end{array}\right. $

(A3)

利用散射理论和玻恩近似, 任意对称各向异性介质的P波反射系数可以写成散射函数S的表达式^[16]:

$ R_{\mathrm{P}}=\frac{1}{4 \rho \cos ^2 \theta} S $

(A4)

式中: ρ是各向同性背景的密度; θ是入射角。散射函数S表示为：

$ S=\Delta \rho \cos 2 \theta+\Delta C \xi $

(A5)

式中: Δρ表示穿过界面的密度扰动; ΔC是刚度元素的扰动; ξ是关于慢度和偏振的函数^[16]。

$ \begin{aligned} & \xi_{11}=\frac{\rho \sin ^4 \theta \cos ^4 \varphi}{M} \quad \xi_{12}=\frac{\rho \sin ^4 \theta \sin ^2 \varphi \cos ^2 \varphi}{M} \quad \xi_{13}=\frac{\rho \sin ^2 \theta \cos ^2 \theta \cos ^2 \varphi}{M} \\ & \xi_{21}=\xi_{12} \quad \xi_{22}=\frac{\rho \sin ^4 \theta \sin ^4 \varphi}{M} \quad \xi_{23}=\frac{\rho \sin ^2 \theta \cos ^2 \theta \sin ^2 \varphi}{M} \quad \xi_{31}=\xi_{13} \\ & \xi_{32}=\xi_{23}, \xi_{33}=\frac{\rho \cos ^4 \varphi}{M} \quad \xi_{44}=\frac{-4 \rho \sin ^2 \theta \cos ^2 \theta \sin ^2 \varphi}{M} \\ & \xi_{55}=\frac{-4 \rho \sin ^2 \theta \cos ^2 \theta \cos ^2 \varphi}{M} \quad \xi_{66}=\frac{4 \rho \sin ^4 \theta \sin ^2 \varphi \cos ^2 \varphi}{M} \end{aligned} $

(A6)

式中: φ是方位角。结合方程(A3)至方程(A6), 可进一步推导与弹性参数、裂缝弱度有关的线性化P波反射系数方程如下:

$ R_{\mathrm{pp}}^{\mathrm{VTI}}(\theta)=\frac{1}{4} \sec ^2 \theta \frac{\Delta M}{M}-2 g \sin ^2 \theta \frac{\Delta \mu}{\mu}+\frac{\cos 2 \theta}{4 \cos ^2 \theta} \frac{\Delta \rho}{\rho}+\frac{1}{4 \cos ^2 \theta} \Delta \delta_{\mathrm{N}}\left(-1-4 g^2 \sin ^4 \theta+4 g \sin ^2 \theta\right)+g \sin ^2 \theta \Delta \delta_{\mathrm{T}} $

(A7)

附录B 含裂缝弱度参数的四项式简化公式推导

对方程(A7)进行改写, 根据三角函数关系sec²θ=1+tan²θ以及cos2θ=2cos²θ－1, 分别将系数项tan²θ, sin²θ, sin²θtan²θ和1进行合并, 则可得:

$ \begin{gathered} R_{\mathrm{PP}}^{\mathrm{VTI}}(\theta)=\frac{1}{4}\left(1+\tan ^2 \theta\right) \frac{\Delta M}{M}-2 g \sin ^2 \theta \frac{\Delta \mu}{\mu}+\left[\frac{1}{2}-\frac{1}{4}\left(1+\tan ^2 \theta\right)\right] \frac{\Delta \rho}{\rho}+\frac{1}{4}\left(1+\tan ^2 \theta\right) \Delta \delta_{\mathrm{N}} \cdot \\ \left(-1-4 g^2 \sin ^4 \theta+4 g \sin ^2 \theta\right)+g \sin ^2 \theta \Delta \delta_{\mathrm{T}} \\ =\frac{1}{4} \frac{\Delta M}{M}+\frac{1}{4} \tan ^2 \theta \frac{\Delta M}{M}-2 g \sin ^2 \theta \frac{\Delta \mu}{\mu}+\frac{1}{4} \frac{\Delta \rho}{\rho}-\frac{1}{4} \tan ^2 \theta \frac{\Delta \rho}{\rho}+ \\ \left(\frac{1}{4}+\frac{1}{4} \tan ^2 \theta\right) \Delta \delta_{\mathrm{N}}\left(-1-4 g^2 \sin ^4 \theta+4 g \sin ^2 \theta\right)+g \sin ^2 \theta \Delta \delta_{\mathrm{T}} \\ =\frac{1}{4} \frac{\Delta M}{M}+\frac{1}{4} \frac{\Delta \rho}{\rho}-\frac{1}{4} \Delta \delta_{\mathrm{N}}-2 g \sin ^2 \theta \frac{\Delta \mu}{\mu}+ \\ g \sin ^2 \theta \Delta \delta_{\mathrm{N}}+g \sin ^2 \theta \Delta \delta_{\mathrm{T}}+\frac{1}{4} \tan ^2 \theta \frac{\Delta M}{M}-\frac{1}{4} \tan ^2 \theta \frac{\Delta \rho}{\rho}- \\ \frac{1}{4} \tan ^2 \theta \Delta \delta_{\mathrm{N}}-g(g-1) \sin ^2 \theta \tan ^2 \theta \Delta \delta_{\mathrm{N}} \\ =\frac{1}{4}\left(\frac{\Delta M}{M}+\frac{\Delta \rho}{\rho}-\Delta \delta_{\mathrm{N}}\right)-2 g \sin ^2 \theta\left(\frac{\Delta \mu}{\mu}-\frac{1}{2} \Delta \delta_{\mathrm{N}}-\frac{1}{2} \Delta \delta_{\mathrm{T}}\right)+ \\ \frac{1}{4} \tan ^2 \theta\left(\frac{\Delta M}{M}-\frac{\Delta \rho}{\rho}-\Delta \delta_{\mathrm{N}}\right)-g(g-1) \sin ^2 \theta \tan ^2 \theta \Delta \delta_{\mathrm{N}} \end{gathered} $

(B1)

整理方程(B1)可得方程(B2), 方程可应用于具有弱各向异性特征地质条件的叠前多参数反演。

$ \begin{gathered} R_{\mathrm{PP}}^{\mathrm{VTI}}(\theta)=\frac{1}{4}\left(\frac{\Delta M}{M}+\frac{\Delta \rho}{\rho}-\Delta \delta_{\mathrm{N}}\right)-\\ 2 g \sin ^2 \theta\left(\frac{\Delta \mu}{\mu}-\frac{1}{2} \Delta \delta_{\mathrm{N}}-\frac{1}{2} \Delta \delta_{\mathrm{T}}\right)+\frac{1}{4} \tan ^2 \theta\left(\frac{\Delta M}{M}-\frac{\Delta \rho}{\rho}-\Delta \delta_{\mathrm{N}}\right)- \\ g(g-1) \sin ^2 \theta \tan ^2 \theta \Delta \delta_{\mathrm{N}} \end{gathered} $

(B2)

致谢: 感谢中南大学潘新朋教授在论文撰写过程中的帮助。