地球物理反演的本质是依据地震数据获取我们所需地层的物性、弹性参数, 揭示地下地层空间展布和含油气性等特征。OSTRANDER^[1]初次提出AVO技术, 揭示地震反射波振幅反映储层特征的奥秘, 由于Zoeppritz方程^[2]精确求解的复杂性, 许多学者开展了基于Zoeppritz方程的AVO方法研究。AKI等^[3]在《定量地震学》中, 对Zoeppritz方程进行了简化近似。拉梅参数(λ, μ)是储层预测和流体识别中最常用的弹性参数, 分别表征岩石的压缩性和刚性, 相较纵、横波速度和密度参数具有更好的预测及识别能力^[4]。GOODWAY等^[5]将Fatti近似式改写为模量组合和密度表示的形式, 通过提取纵、横波阻抗相对变化量求取纵、横波阻抗信息, 进而依据公式计算λρ和μρ, 但其本质仍是一种间接反演方法, 而间接计算产生的累积误差会降低预测精度^[6-7]。XU等^[8]提出利用体积模量K和拉梅常数λ, μ表征的完全隐含速度信息的近似方程。在密度相对变化较小时, 能够利用两项参数反映入射角不超过45°的反射特性。GRAY等^[9]基于Aki-Richards近似方程提出了由体积模量/压缩模量、剪切模量和密度的相对变化率表示的反射系数近似公式, 实现利用叠前地震数据直接反演出λ和μ, 可以更加准确地识别岩性和流体。CONNOLLY^[10]率先提出弹性阻抗的概念, 能够获得纵、横波速度和密度信息, 进而计算其它弹性参数。王保丽等^[11]基于Gray近似方程, 推导了由λ、μ和ρ三参数表达式表示的弹性阻抗方程, 并进行标准化处理, 可以直接从弹性阻抗数据中提取拉梅参数。宗兆云等^[12]基于平面波入射假设条件, 推导了基于杨氏模量Y、泊松比σ和密度D的YPD反射系数近似方程, 实现目标参数的直接反演, 能够有效识别页岩气“甜点”。

URSIN等^[13]认为利用直接从常规纵波地震数据反演纵、横波速度和密度3个参数难度较大。张春涛等^[14]认为利用纵波地震数据一般很难反演得到密度比。叠前AVO反演能够利用更多地下真实数据, 其本质仍是不适定问题, 一般通过删除第三项提高反演问题的稳定性^[15]。SMITH等^[16]利用Gardner经验关系式^[17]改写Aki-Richards近似方程, 消除密度项, 虽然稳定性有所提高, 但经验关系式与实际地层的吻合情况严重影响反演的精确度。

在Zoeppritz方程的基础上, 提出一种关于剪切模量和密度的AVO直接反演方法, 该方法直接提取剪切模量和密度, 避免间接计算引起的误差。本方法仅涉及两参数反演, 且未引入经验关系式, 具有比常规三参数反演更好的稳定性。通过模型数据试算与方法对比, 验证了此方法正、反演问题的有效性与稳定性。最后, 将其应用于实际叠前地震数据, 以验证方法在岩性识别与储层预测中的正确性。

1 方法原理 1.1 AVO正演公式

ZOEPPRITZ^[2]基于平面波理论, 根据斯奈尔(Snell)定理、界面两侧位移和应力的连续性, 建立了用位移振幅比表示的反射透射系数, 对Zoeppritz方程组进行简化推导:

$ \left[\begin{array}{cccc} \sin \alpha_1 & \cos \beta_1 & -\sin \alpha_2 & \cos \beta_2 \\ \cos \alpha_1 & -\sin \beta_1 & \cos \alpha_2 & \sin \beta_2 \\ \sin 2 \alpha_1 & \frac{v_{\mathrm{P} 1}}{v_{\mathrm{S} 1}} \cos 2 \beta_1 & \frac{v_{\mathrm{P} 1} v_{\mathrm{S} 2}^2 \rho_2}{v_{\mathrm{P} 2} v_{\mathrm{S} 1}^2 \rho_1} \sin 2 \alpha_2 & -\frac{v_{\mathrm{P} 1} v_{\mathrm{S} 2} \rho_2}{v_{\mathrm{S} 1}^2 \rho_1} \cos 2 \beta_2 \\ \cos 2 \beta_1 & -\frac{v_{\mathrm{S} 1}}{v_{\mathrm{P} 1}} \sin 2 \beta_1 & -\frac{v_{\mathrm{P} 2} \rho_2}{v_{\mathrm{P} 1} \rho_1} \cos 2 \beta_2 & -\frac{v_{\mathrm{S} 2} \rho_2}{v_{\mathrm{P} 1} \rho_1} \sin 2 \beta_2 \end{array}\right]\left[\begin{array}{c} R_{\mathrm{PP}} \\ R_{\mathrm{PS}} \\ T_{\mathrm{PP}} \\ T_{\mathrm{PS}} \end{array}\right]=\left[\begin{array}{c} -\sin \alpha_1 \\ \cos \alpha_1 \\ \sin 2 \alpha_1 \\ -\cos 2 \beta_1 \end{array}\right] $

(1)

式中: R_PP、T_PP、R_PS、T_PS分别为纵波反射系数、纵波透射系数、转换横波反射系数和转换横波透射系数; v_P1、v_P2、v_S1、v_S2分别为界面两侧的纵、横波速度; α₁、β₁、α₂、β₂分别为界面两侧的纵、横波的反射和透射角; ρ₁、ρ₂分别为界面上、下地层的密度。(1)式左右两侧矩阵中的速度、角度项均满足Snell定律。则有:

$ \frac{\sin \alpha_1}{v_{\mathrm{P} 1}}=\frac{\sin \alpha_2}{v_{\mathrm{P} 2}}=\frac{\sin \beta_1}{v_{\mathrm{S} 1}}=\frac{\sin \beta_2}{v_{\mathrm{S} 2}}=P $

(2)

式中: P为射线参数。记x=sinα₁, a=v_P2/v_P1, b=v_S2/v_S1, c=ρ₂/ρ₁, γ₁=v_S1/v_P1, 代入(1)式整理可得:

$ \begin{gathered} \left[\begin{array}{cccc} x & \sqrt{1-\gamma_1^2 x^2} & -a x & \sqrt{1-b^2 \gamma_1^2 x^2} \\ \sqrt{1-x^2} & -\gamma_1 x & \sqrt{1-a^2 x^2} & b \gamma_1 x \\ 2 \gamma_1 x \sqrt{1-x^2} & 1-2 \gamma_1^2 x^2 & 2 b^2 c \gamma_1 x \sqrt{1-a^2 x^2} & -b c\left(1-2 b^2 \gamma_1^2 x^2\right) \\ 1-2 \gamma_1^2 x^2 & -2 \gamma_1^2 x \sqrt{1-\gamma_1^2 x^2} & -a c\left(1-2 b^2 \gamma_1^2 x^2\right) & -2 b^2 c \gamma^2 x \sqrt{1-b^2 \gamma_1^2 x^2} \end{array}\right]\left[\begin{array}{l} R_{\mathrm{PP}} \\ R_{\mathrm{PS}} \\ T_{\mathrm{PP}} \\ T_{\mathrm{PS}} \end{array}\right]\\ =\left[\begin{array}{c} -x \\ \sqrt{1-x^2} \\ 2 \gamma_1 x \sqrt{1-x^2} \\ 2 \gamma_1^2 x^2-1 \end{array}\right] \end{gathered} $

(3)

简化矩阵形式为:

$ \boldsymbol{A}\left[\begin{array}{c} R_{\mathrm{PP}} \\ R_{\mathrm{PS}} \\ T_{\mathrm{PP}} \\ T_{\mathrm{PS}} \end{array}\right]=\boldsymbol{B} $

(4)

主要对纵波反射系数进行重新推导, 根据克莱姆法可知:

$ R_{\mathrm{PP}}=\frac{\operatorname{det} \boldsymbol{A}_{11}}{\operatorname{det} \boldsymbol{A}} $

(5)

det表示求取矩阵的行列式, 其中:

$ \boldsymbol{A}_{11}=\left[\begin{array}{cccc} -x & \sqrt{1-\gamma_1^2 x^2} & -a x & \sqrt{1-b^2 \gamma_1^2 x^2} \\ \sqrt{1-x^2} & -\gamma_1 x & \sqrt{1-a^2 x^2} & b \gamma_1 x \\ 2 \gamma_1 x \sqrt{1-x^2} & 1-2 \gamma_1^2 x^2 & 2 b^2 c \gamma_1 x \sqrt{1-a^2 x^2} & -b c\left(1-2 b^2 \gamma_1^2 x^2\right) \\ 2 \gamma_1^2 x^2-1 & -2 \gamma_1^2 x \sqrt{1-\gamma_1^2 x^2} & -a c\left(1-2 b^2 \gamma_1^2 x^2\right) & -2 b^2 c \gamma^2 x \sqrt{1-b^2 \gamma_1^2 x^2} \end{array}\right] $

(6)

式中: A₁₁是将A中的第一列用B替代所得到的矩阵。借鉴AVO公式常用各类参数相对变化率表示反射系数的方式, 引入横波、密度的相对变化率, 记R_S=(v_S2－v_S1)/(v_S2+v_S1), R_ρ=(ρ₂－ρ₁)/(ρ₂+ρ₁), 可对(3)式中参数进行转换:

$ \begin{aligned} &a=\frac{v_{\mathrm{P} 2}}{v_{\mathrm{P} 1}}=\frac{v_{\mathrm{S} 2} / \gamma_2}{v_{\mathrm{S} 1} / \gamma_1}=\frac{\gamma_1}{\gamma_2} b \\ &b=\frac{v_{\mathrm{S} 2}}{v_{\mathrm{S} 1}}=\frac{1+\frac{v_{\mathrm{S} 2}-v_{\mathrm{S} 1}}{v_{\mathrm{S} 2}+v_{\mathrm{S} 1}}}{1-\frac{v_{\mathrm{S} 2}-v_{\mathrm{S} 1}}{v_{\mathrm{S} 2}+v_{\mathrm{S} 1}}}=\frac{1+R_{\mathrm{S}}}{1-R_{\mathrm{S}}} \\ &c=\frac{\rho_2}{\rho_1}=\frac{1+\frac{\rho_2-\rho_1}{\rho_2+\rho_1}}{1-\frac{\rho_2-\rho_1}{\rho_2+\rho_1}}=\frac{1+R_\rho}{1-R_\rho} \end{aligned} $

(7)

对(7)式泰勒展开可得:

$ \begin{array}{l} b=\frac{1+R_S}{1-R_S}=\frac{2}{1-R_S}-1=\frac{2}{1-\frac{R_\mu-R_\rho}{2}}-1= \\ 1+\left(R_\mu-R_\rho\right)+\frac{\left(R_\mu-R_\rho\right)^2}{2}+\frac{\left(R_\mu-R_\rho\right)^3}{4}+\cdots \\ c=\frac{1+R_\rho}{1-R_\rho}=\frac{2}{1-R_\rho}-1=1+2 R_\rho+2 R_\rho^2+2 R_\rho^3+\cdots \\ \;\;\;\;\;\;a=T b=T\left[1+\left(R_\mu-R_\rho\right)+\frac{\left(R_\mu-R_\rho\right)^2}{2}+\right. \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left.\frac{\left(R_\mu-R_\rho\right)^3}{4}+\cdots\right] \end{array} $

(8)

式中: R_μ=(μ₂－μ₁)/(μ₂+μ₁), 表示界面上、下剪切模量的相对变化率; T=γ₁/γ₂, 表示界面上、下地层横、纵波速度比的比值, 可通过测井曲线获得。

同样展开角度项, 则有:

$ \begin{array}{l} \sqrt{1-x^2}= 1-\frac{1}{2} x^2-\frac{1}{2 \times 4} x^4-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \frac{1 \times 3}{2 \times 4 \times 6} x^6-\cdots \\ \sqrt{1-\gamma^2 x^2}= 1-\frac{1}{2} \gamma^2 x^2-\frac{1}{2 \times 4} \gamma^4 x^4-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \frac{1 \times 3}{2 \times 4 \times 6} \gamma^6 x^6-\cdots \\ \sqrt{1-a^2 x^2}= 1-\frac{1}{2} a^2 x^2-\frac{1}{2 \times 4} a^4 x^4-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \frac{1 \times 3}{2 \times 4 \times 6} a^6 x^6-\cdots \\ \sqrt{1-b^2 \gamma^2 x^2}= 1-\frac{1}{2} b^2 \gamma^2 x^2-\frac{1}{2 \times 4} b^4 \gamma^4 x^4-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \frac{1 \times 3}{2 \times 4 \times 6} b^6 \gamma^6 x^6-\cdots \end{array} $

(9)

将(8)式、(9)式代入(5)式, 为保持精度, 角度项保留至4次方项, 整理可得:

$ \begin{gathered} R_{\mathrm{PP}}(\theta)=\frac{1}{2(T+1)}\left[\left(A_1+A_2 \sin ^2 \theta_1+\right.\right. \\ \left.A_3 \sin ^4 \theta_1\right) R_\mu+\left(B_1+B_2 \sin ^2 \theta_1+\right. \\ \left.\left.B_3 \sin ^4 \theta_1\right) R_\rho+C_1\right] \end{gathered} $

(10)

其中,

$ \begin{gathered} A_1=B_1=T+1 \\ A_2=\left(-2 T+4 T^2\right)+(-7-9 T) \gamma_1^2 \\ A_3=\left(-T^2+2 T^3\right)+\left(T-T^2\right) \gamma_1^2+\frac{(1-T) \gamma_1^4}{2} \\ B_2=(-2 T)+(-1+T) \gamma_1^2 \\ B_3=\left(-T^2\right)+\left(T^2-T\right) \gamma_1^2+\frac{(T-1) \gamma_1^4}{2} \\ C_1=2(T-1)\left(1+T \sin ^2 \theta_1+\frac{T^2}{2} \sin ^4 \theta_1\right) \end{gathered} $

(10) 式为由剪切模量、密度的相对变化率表示的反射系数方程, 记为μ-ρ方程, 对其进行化简得:

$ R_{\mathrm{PP}}(\theta)=A R_\mu+B R_\rho+C $

(11)

其中,

$ \begin{gathered} A=\frac{A_1+A_2 \sin ^2 \theta_1+A_3 \sin ^4 \theta_1}{2(T+1)} \\ B=\frac{B_1+B_2 \sin ^2 \theta_1+B_3 \sin ^4 \theta_1}{2(T+1)} \\ C=\frac{C_1}{2(T+1)} \end{gathered} $

假定叠前地震数据有m个入射角, 每一道地震数据有n个时间采样点, 结合褶积模型, 反射系数与子波褶积得到合成地震记录, 即正演方程为:

$ \left[ \begin{matrix} {{d}^{1}}\left( {{\theta }_{1}} \right) \\ \vdots \\ {{d}^{n}}\left( {{\theta }_{1}} \right) \\ \vdots \\ {{d}^{1}}\left( {{\theta }_{m}} \right) \\ \vdots \\ {{d}^{n}}\left( {{\theta }_{m}} \right) \\ \end{matrix} \right]=\left[ \begin{matrix} W\left( {{\theta }_{1}} \right) & {} & {} \\ {} & \ddots & {} \\ {} & {} & W\left( {{\theta }_{1}} \right) \\ {} & \vdots & {} \\ W\left( {{\theta }_{m}} \right) & {} & {} \\ {} & \ddots & {} \\ {} & {} & W\left( {{\theta }_{m}} \right) \\ \end{matrix} \right]\left[ \begin{matrix} R_{\text{PP}}^{1}\left( {{\theta }_{1}} \right) \\ \vdots \\ R_{\text{PP}}^{n}\left( {{\theta }_{1}} \right) \\ \vdots \\ R_{\text{PP}}^{1}\left( {{\theta }_{m}} \right) \\ \vdots \\ R_{\text{PP}}^{n}\left( {{\theta }_{m}} \right) \\ \end{matrix} \right] $

(12)

式中: d(θ)为合成地震数据; W(θ)为子波矩阵; R_PP(θ)为PP波反射系数; 上标均代表对应时间采样点; 下标均代表对应入射角。

1.2 与Aki-Richards近似方程的比较

AKI等^[3]在《定量地震学》中给出了Zoeppritz方程的经典近似式为:

$ R(\bar{\theta}) \approx \sec ^2 \bar{\theta} R_{\mathrm{P}}-8 \bar{\gamma}^2 \sin ^2 \bar{\theta} R_{\mathrm{S}}+\left(1-4 \bar{\gamma}^2 \sin ^2 \bar{\theta}\right) R_\rho $

(13)

式中: θ和γ分别表示波在界面处的入射角和透射角的平均角度和横、纵波速度比的均值。

由于v_P=v_S/γ, 结合纵波速度与剪切模量的关系, 得:

$ R_{\mathrm{P}}=R_{\mathrm{S}}-R_\gamma, R_{\mathrm{S}}=\frac{R_\mu-R_\rho}{2} $

(14)

将(14)式代入(13)式, Aki-Richards近似方程可改写为:

$ \begin{gathered} R(\bar{\theta}) \approx \frac{\sec ^2 \bar{\theta}-8 \bar{\gamma}^2 \sin ^2 \bar{\theta}}{2} R_\mu+\frac{2-\sec ^2 \bar{\theta}}{2} R_\rho- \\ \sec ^2 \bar{\theta} R_\gamma \end{gathered} $

(15)

其中,

$ \bar{\theta}=\frac{\theta_1+\theta_2}{2}=\frac{\theta_1+\arcsin \left(a \sin \theta_1\right)}{2} $

(16)

对(16)式两端取正弦得:

$ \sin \bar{\theta}=\sqrt{\frac{1-\cos (2 \bar{\theta})}{2}}=\sqrt{\frac{1-\cos \theta_1 \cos \theta_2+\sin \theta_1 \sin \theta_2}{2}} $

(17)

对(17)式右端余弦项泰勒展开

$ \begin{gathered} \cos \theta_1=\sqrt{1-\sin ^2 \theta_1}=1-\frac{1}{2} \sin ^2 \theta_1-\frac{1}{2 \times 4} \sin ^4 \theta_1- \\ \frac{1 \times 3}{2 \times 4 \times 6} \sin ^6 \theta_1-\cdots \end{gathered} $

(18)

$ \begin{aligned} \cos \theta_2=& \sqrt{1-a^2 \sin ^2 \theta_1}=1-\frac{1}{2} a^2 \sin ^2 \theta_1-\frac{1}{2 \times 4} \cdot \\ & a^4 \sin ^4 \theta_1-\frac{1 \times 3}{2 \times 4 \times 6} a^6 \sin ^6 \theta_1-\cdots \end{aligned} $

(19)

将(18)式和(19)式代入(17)式且忽略高阶项, 则有:

$ \sin \bar{\theta}=\sin \theta_1\left(\frac{1}{2}+\frac{1}{2} a\right) $

(20)

将(8)式代入(20)式

$ \sin \bar{\theta}=\sin \theta_1\left[\frac{1+T}{2}+\frac{T}{2}\left(R_\mu-R_\rho\right)\right] $

(21)

$ \sec ^2 \bar{\theta}=\frac{1}{1-\sin ^2 \bar{\theta}}=1+\sin ^2 \bar{\theta}+\sin ^4 \bar{\theta}+\sin ^6 \bar{\theta}+\cdots $

(22)

将(21)式和(22)式代入(15)式, 保留相对变化率的一阶近似, 有:

$ \begin{aligned} &R(\bar{\theta}) \approx \\ &\frac{1+\left(1-8 \gamma^2\right)\left(\frac{1+T}{2}\right)^2 \sin ^2 \theta_1+\left(\frac{1+T}{2}\right)^4 \sin ^4 \theta_1}{2} R_\mu+ \\ &\frac{1-\left(\frac{1+T}{2}\right)^2 \sin ^2 \theta_1-\left(\frac{1+T}{2}\right)^4 \sin ^4 \theta_1}{2} R_\rho- \\ &{\left[1+\left(\frac{1+T}{2}\right)^2 \sin ^2 \theta_1+\left(\frac{1+T^4}{2} \sin ^4 \theta_1\right] R_\gamma\right.} \end{aligned} $

(23)

其中,

$ R_\gamma=\frac{\gamma_2-\gamma_1}{\gamma_2+\gamma_1}=\frac{1-T}{1+T} $

考虑Aki-Richards近似方程成立的假设条件为相邻两层介质的弹性参数变化较小, 则

$ \begin{gathered} \bar{\gamma}=\frac{v_{\mathrm{S} 1}+v_{\mathrm{S} 2}}{v_{\mathrm{P} 1}+v_{\mathrm{P} 2}}=\frac{2 v_{\mathrm{S} 1}+\Delta v_{\mathrm{S}}}{2 v_{\mathrm{P} 1}+\Delta v_{\mathrm{P}}}=\gamma_1+ \\ \frac{\Delta v_{\mathrm{S}}-\gamma_1 \Delta v_{\mathrm{P}}}{2 v_{\mathrm{P} 1}+\Delta v_{\mathrm{P}}} \approx \gamma_1 \end{gathered} $

(24)

(23) 式可表示为:

$ \begin{gathered} R(\bar{\theta})=\frac{1+\left(1-8 \bar{\gamma}^2\right) \sin ^2 \theta_1+\sin ^4 \theta_1}{2} R_\mu+ \\ \frac{1-\sin ^2 \theta_1-\sin ^4 \theta_1}{2} R_\rho \end{gathered} $

(25)

比较(10)式和(25)式可以发现, 假设相邻两层介质的弹性参数变化较小, Aki-Richards近似方程在外推过程中, γ的平均值用界面上层的γ替代, 从而得到Aki-Richards近似方程等价于μ-ρ方程, 但μ-ρ方程在推导过程中没有进行平均值的近似, 避免了这部分近似的误差。

2 正演模拟

为测试μ-ρ方程的可行性, 基于CASTAGNA等^[18]提出的AVO分类方法建立了4类AVO模型(表 1至表 4), 且都利用Zoeppritz方程和μ-ρ方程(图例以“μ-ρ”表示)、Aki-Richards近似方程和Gray近似方程计算4种模型界面的纵波反射系数并计算各类方法与精确值之间的误差(图 1至图 4)。图 1a、图 2a、图 3a和图 4a分别是第Ⅰ、Ⅱ、Ⅲ、Ⅳ类AVO模型的正演结果; 图 1b、图 2b、图 3b和图 4b分别是三类近似方法求得的第Ⅰ、Ⅱ、Ⅲ、Ⅳ类AVO模型的反射系数与Zoeppritz方程求得的准确值之间的误差(已校正法向入射时初始误差)。对比3种方法的正演结果可以看出, 在0入射时, μ-ρ方程与Gray近似方程的正演结果均与准确值存在不同程度的误差, 说明在法向入射时两种方法均存在直接引入弹性参数的转化误差。但就整体趋势而言, μ-ρ方程正演结果在4类模型中均与实际模型趋势吻合最好, 误差图更加直观地表明, 误差基本来自法向入射时计算误差, 在校正初始误差后, 本文提出的方法明显比其它两种方法更为精确。此外, Aki-Richards近似方程和Gray近似方程在4类模型中的正演结果具有相同的趋势, 且随着入射角度的增加两者逐步偏离精确Zoeppritz方程计算结果, 其根本原因是Gray近似是基于Aki-Richards方程参数转换所得, 本质上基于一种近似的不同表示, 两种方法涉及的角度项仅保留至二阶, 在大角度入射时易产生较大误差, 新方法的角度项保留至四阶, 能够较好地避免此类误差。

为验证μ-ρ方程的稳定性, 对建立的4类AVO模型(表 1至表 4)进行条件数分析。图 5为4类AVO模型采用3种方法正演时, 条件数随最大入射角度的变化情况, 其中最小入射角度保持为0, 最大入射角度范围为5°~40°。由图 5可以看出, 当最大入射角度比较小时, 3种方法的条件数都较大, 表明过小角度入射时反演问题非常不稳定, 因此需要获取较大的偏移量信息以提高反演的稳定性。同时可以看出, μ-ρ方程所得条件数比以往3项参数的近似方程条件数成指数级降低, 表明μ-ρ方程的稳定性远高于Aki-Richards近似方程和Gray近似方程。

3 反演分析 3.1 一维层状模型模拟

为验证μ-ρ方程反演的可行性, 构建一维五层模型进行反演模拟。反演算法基于最小二乘法, 模型数据如表 5所示。利用Zoeppritz方程获得模型精确数据, 分别利用μ-ρ方程、Aki-Richards近似方程和Gray近似方程对精确数据进行μ、ρ反演并求取各方法计算结果与精确值之间的误差, 结果如图 6和图 7所示。对比3种方法的反演结果, 两种三参数近似方程在881~1040m层段与真实数据拟合相对较好。当深度较大时, 反演误差增大, 认为是反演迭代中误差累积所致。对比3种方法的结果认为, μ-ρ方程直接反演得到的μ、ρ与真实模型的拟合误差最小, 原因在于此方法直接提取弹性参数, 减少了累积误差且两参数反演具有更高的稳定性。Aki-Richards近似方程和Gray近似方程在层状模型反演时仍具有相同的趋势, 再次验证两种近似式本质上是基于一种近似的不同表示。

假设井数据丰富, 可以通过井位速度信息求取泊松比, 结合μ-ρ方程直接反演得到μ, 可以通过拉梅参数与泊松比的关系获得λ。在本次模型模拟中代入每层真实泊松比, 按上述方法计算λ, 并利用Gray近似方程对精确数据进行直接反演得到λ, 结果如图 8所示。由图 8可以看出, 此方法求解出的λ与真实数据吻合程度更高, 表明在井数据丰富条件下, μ-ρ方程能准确获得λ。

3.2 反演流程

在实际地震反演中, 通常采用离散方程进行计算, 正演公式(12)可直接简写为矩阵形式:

$ \boldsymbol{d}=\boldsymbol{W r} $

(26)

式中: d为叠前地震数据, W为子波矩阵, r为反射系数矩阵, 即(11)式。

基于最小二乘原则约束反演误差, 考虑L₀范数约束下的反演问题直接求解是NP-hard^[19], 采用L₁范数构建稀疏脉冲反演的目标函数为:

$ J(r)=\|\boldsymbol{d}-\boldsymbol{W} \boldsymbol{r}\|_2^2+l_1\|\boldsymbol{r}\|_1 $

(27)

式中: ‖·‖₂表示测量反演误差的L₂范数, ‖·‖₁表示反射系数的L₁范数, 用于测量反射系数序列的稀疏性; l₁是L₁范数的正则化因子。

公式(27)中, L₁范数测量了反射系数序列的垂向稀疏性, 并没有考虑多道反射系数序列的横向变化和连续性, 因此引入核模约束^[20], 以约束r为低秩, 保证r相邻道具有相关性, 其定义为矩阵r所有奇异值的和, 即

$ \|\boldsymbol{r}\|_*=\|\varLambda(\boldsymbol{r})\|_1=\sum\limits_i \sigma(i) $

(28)

式中: σ(i)为矩阵r的第i个奇异值, 则(27)式可以改写为:

$ J(\boldsymbol{r})=\|\boldsymbol{d}-\boldsymbol{W} \boldsymbol{r}\|_2^2+l_1\|\boldsymbol{r}\|_1+l_2\|\boldsymbol{r}\|_* $

(29)

式中: l₂为核模约束的正则化参数。在实际数据处理中, 仅通过L₁范数的稀疏约束和核模约束并不能完全准确地恢复低频分量。因此, 将模型参数的先验信息r_priori作为正则化约束项加入到目标函数中, 目标函数(29)可写为

$ \begin{gathered} J(\boldsymbol{r})=\|\boldsymbol{d}-\boldsymbol{W} \boldsymbol{r}\|_2^2+l_1\|\boldsymbol{r}\|_1+l_2\|\boldsymbol{r}\|_*+ \\ l_3\left\|\boldsymbol{r}-\boldsymbol{r}_{\text {priori }}\right\|_2^2 \end{gathered} $

(30)

式中: l₃为先验模型的正则化约束系数。(30)式即为最终的叠前反演目标函数。可以通过迭代重加权最小二乘算法^[21]进行求解。

3.3 应用实例

为验证μ-ρ方程在实际资料反演中的可行性, 选取某工区的实际地震数据进行测试。该工区的叠前地震数据包括3个不同角度范围的部分角度叠加数据体, 叠加剖面如图 9所示, 其中, 图 9a为0°~10°叠加剖面, 图 9b为8°~17°叠加剖面, 图 9c为15°~25°叠加剖面。首先, 对工区内的测井资料进行预处理。然后, 利用求解基于μ-ρ方程的反演目标函数实现剪切模量和密度参数的反演, 反演结果如图 10和图 11所示。图 10为剪切模量μ的反演剖面与井位反演曲线。图 10a中红色曲线为井的位数据经滤波后计算的μ曲线; 左侧为井的岩性图。由图 10可以看出, 反演剖面的高μ区与井的砂岩段吻合较好, 说明通过反演μ可以较好地识别岩性。图 10b为井位位置的反演结果。其中红色曲线为μ-ρ方程反演结果, 证实基于μ-ρ方程的反演结果与测井曲线吻合较好。图 11为密度ρ的反演剖面与井位反演曲线, 井位反演结果与测井曲线基本一致。图 10和图 11表明, 基于μ-ρ方程得到的反演剖面与测井解释结果基本吻合, 能得到较为精确的μ和ρ值, 将其与其它参数联合分析、构建流体因子, 可更精确地识别岩性和流体。

4 结论

拉梅参数是储层预测和流体识别中最常用的弹性参数, 目前常用的提取方法是通过反演纵、横波速度间接计算获得, 这会导致误差累积。同时, 利用叠前地震数据实现三参数的精确反演仍是一个难题。我们提出了一种基于剪切模量和密度的两参数纵波反射系数方程, 并将该方程应用于实际资料处理, 能够较好地直接反演出剪切模量和密度。正演模拟和反演分析结果表明:

1) 从Zoeppritz方程出发, 直接推导由剪切模量和密度相对变化率表征的反射系数方程, 未引入经验关系式, 且角度项保留至四阶, 在较大角度入射时同样具有较好的正演效果。避免了反演出现累积误差。

2) μ-ρ方程涉及两参量, 条件数远低于Aki-Richards方程与Gray近似方程, 具有更高的稳定性。

3) 在井资料丰富的情况下, 基于μ-ρ方程精确反演得到μ, 进而结合井位信息以及拉梅参数与泊松比的关系式获得λ, 实现三参数的获取。

4) 引入核模约束项以及先验信息约束项构建反演目标函数, 能够增强反演的横向连续性并还原原始数据的频率信息, 实际资料处理证实本方法能较好地拟合井数据, 反演得到较为精确的μ, ρ值, 并用于岩性识别。

5) 正演模拟与反演分析结果均表明, μ-ρ方程具有较高的稳定性和准确性, 基于μ-ρ方程反演得到的弹性参数可与其它参数联合分析、构建流体因子, 可更精确地识别岩性和流体。