0 引言

对于“亮点”型气层，可以利用地震数据的高振幅异常进行识别^[1]；对于振幅差异不明显的储层，则需要利用地震数据的其他特征，如弹性参数进行识别。泊松比在AVO理论出现的伊始就被用于指示含气储层，并广泛应用于碎屑岩储层的预测^[2]。Connolly^[3]提出弹性阻抗的概念，并证明了利用不同角度弹性阻抗的差异可以指示含气层的存在；Gray^[4]利用拉梅参数探测油气储层，并建立了反射系数与拉梅参数之间的关系；Russell等^[5]基于孔隙弹性介质理论，从Gassmann方程中推导出流体项；Sena等^[6]在页岩储层的研究中发现杨氏模量和泊松比可以作为有利储层的指示参数。

叠前地震反演是基于地震资料获取此类参数的有效途径^[7]，但是地震数据的反演通常是一个病态问题，为了得到合理稳定的解，需要对问题进行正则化，稀疏正则化是目前解决地震反演问题常用方法之一^[8]。L₀范数正则化虽然有较好的稀疏性，但求解L₀范数稀疏最优化是公认的难度较大的问题，需要采用贪婪迭代的策略进行求解^[9]；匹配追踪算法是一种经典的贪婪算法，Nguyen^[10]利用匹配追踪算法分解地震信号、反演反射系数；张繁昌等^[11]基于匹配追踪算法提取地震数据的时频谱特征识别三角洲砂岩层的尖灭线；刘晓晶等^[12]利用基于正交匹配追踪算法的叠前AVO反演方法实现了块化反演。虽然匹配追踪算法在地震反演领域应用广泛，但缺点是空间分辨率低及横向连续性差。为此，Zhang等^[13-14]提出了利用基追踪算法开展地震反演，该方法不仅有较好的横向连续性，而且计算效率高。基追踪可以得到块化的反演结果，Theune等^[15]认为块化反演结果的分辨率更高，对地层的解释能力更强。为了进一步提高基追踪反演的横向连续性，Zhang等^[16]提出了基于空间正则化的多道基追踪反演方法。

致密气储层是一种非常规储层，通常具有低孔、低渗的特征，常用的储层预测方法难以对其进行有效识别^[17-18]。本文借鉴页岩气勘探的成功经验，利用杨氏模量和泊松比表征致密砂岩储层含油气特征，并从Zoeppritz方程的近似公式出发推导出以杨氏模量和泊松比表示的地震反射系数公式。基于该方程开展叠前稀疏层基追踪反演，为了提高反演结果的稳定性和横向连续性，在基追踪的目标函数中加入模型约束项。利用谱投影梯度法算法进行求解，得到稳定的反演结果。模型实验证明了方法的有效性和稳定性。实际工区数据应用不仅得到了稳定连续的反演结果，同时块化反演结果为储层边界的确定提供了数据基础。

1 方法原理

叠前地震反演是进行弹性参数预测常用的方法之一，其基础是地震反射系数方程。由于Zoeppritz方程形式复杂，很少直接用于地震反演。Gray^[4]提出了利用体积模量和剪切模量表示的反射系数近似公式

$ \begin{array}{*{20}{c}} {{R_{{\rm{PP}}}}\left( \theta \right) = \left( {\frac{1}{4} - \frac{1}{3}{\gamma ^2}} \right){{\sec }^2}\theta \frac{{\Delta K}}{K} + {\gamma ^2}\left( {\frac{1}{3}{{\sec }^2}\theta - } \right.}\\ {\left. {2{{\sin }^2}\theta } \right)\frac{{\Delta \mu }}{\mu } + \left( {\frac{1}{2} - \frac{1}{4}{{\sec }^2}\theta } \right)\frac{{\Delta \rho }}{\rho }} \end{array} $

(1)

式中：R_PP(θ)表示入射角θ对应的反射系数；K、ΔK、μ、Δμ、ρ、Δρ分别表示界面两侧介质的体积模量、剪切模量及密度的均值和差值；γ为横、纵波速度比。

根据岩石物理理论，体积模量K、剪切模量μ与杨氏模量E、泊松比σ之间的关系可表示为

$ K = \frac{E}{{3\left( {1 - 2\sigma } \right)}} $

(2)

$ \mu = \frac{E}{{2\left( {1 + \sigma } \right)}} $

(3)

对式(2)和式(3)两边进行差分运算，可得

$ \frac{{\Delta K}}{K} = \frac{{\Delta E}}{E} + \frac{{2\sigma }}{{1 - 2\sigma }}\frac{{\Delta \sigma }}{\sigma } $

(4)

$ \frac{{\Delta \mu }}{\mu } = \frac{{\Delta E}}{E} - \frac{\sigma }{{1 + \sigma }}\frac{{\Delta \sigma }}{\sigma } $

(5)

将式(4)和式(5)代入式(1)，整理得到用杨氏模量、泊松比和密度的相对变化率表示的地震反射系数近似方程

$ {R_{{\rm{pp}}}}\left( \theta \right) = A\left( \theta \right)\frac{{\Delta E}}{E} + B\left( \theta \right)\frac{{\Delta \sigma }}{\sigma } + C\left( \theta \right)\frac{{\Delta \rho }}{\rho } $

(6)

式中：E、σ和ΔE、Δσ分别表示界面两侧杨氏模量和泊松比的平均值和差值; $A(\theta ) = \frac{1}{4}{\sec ^2}\theta - 2{\gamma ^2}{\sin ^2}\theta ;\;\;$ $B(\theta ) = \frac{1}{4}{\sec ^2}\theta \frac{{\left( {2{\gamma ^2} - 3} \right){{\left( {2{\gamma ^2} - 1} \right)}^2}}}{{{\gamma ^2}\left( {4{\gamma ^2} - 3} \right)}} + 2{\gamma ^2}{\sin ^2}\theta \frac{{{{\left( {1 - 2{\gamma ^2}} \right)}^2}}}{{3 - 4{\gamma ^2}}};$ $C(\theta ) = \frac{1}{2} - \frac{1}{4}{\sec ^2}\theta $。

利用式(6)进行叠前数据反演，可以在地震尺度上估计杨氏模量和泊松比，据此进行有利储层的预测和评价。

2 反演方法

根据褶积理论，地震信号可以表示为子波与反射系数的褶积形式。在有噪条件下，三者的关系可以表示为

$ d = WR + e $

(7)

式中：d表示地震数据；W表示子波；e表示噪声。

以地层为研究对象，进行反射系数奇偶分解有利于得到更加清晰的地层边界。地层顶、底界面的反射系数可以表示为奇、偶反射系数对的线性组合^[19]

$ R = \sum\limits_{i = 1}^I {\sum\limits_{j = 1}^J {\left[ {{a_{i,j}}{R_{\rm{e}}}\left( {i,j,\Delta t} \right) + {b_{i,j}}{R_{\rm{o}}}\left( {i,j,\Delta t} \right)} \right]} } $

(8)

式中R_o(i, j, Δt)和R_e(i, j, Δt)分别为奇、偶系数对，表达式如下

$ {R_{\rm{e}}}\left( {i,j,\Delta t} \right) = \delta \left( {t - i\Delta t} \right) + \delta \left( {t - i\Delta t + j\Delta t} \right) $

(9)

$ {R_{\rm{o}}}\left( {i,j,\Delta t} \right) = \delta \left( {t - i\Delta t} \right) + \delta \left( {t - i\Delta t + j\Delta t} \right) $

(10)

式中：Δt表示地震数据采样间隔；i、j表示反射点位置。

式(8)可以表示为矩阵形式

$ R = \left[ {\begin{array}{*{20}{l}} {{R_{\rm{e}}}}&{{R_{\rm{o}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} a\\ b \end{array}} \right] $

(11)

根据式(11)可将式(8)表示为

$ \left[ {\begin{array}{*{20}{c}} {{r_E}}\\ {{r_\sigma }}\\ {{r_\rho }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{R_{\rm{e}}}}&{{R_{\rm{o}}}}&0&0&0&0\\ 0&0&{{R_{\rm{e}}}}&{{R_{\rm{o}}}}&0&0\\ 0&0&0&0&{{R_{\rm{e}}}}&{{R_{\rm{o}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{a_{{\rm{e}}E}}}\\ {{b_{{\rm{o}}E}}}\\ {{a_{{\rm{e}}\sigma }}}\\ {{b_{{\rm{o}}\sigma }}}\\ {{a_{{\rm{e}}\rho }}}\\ {{b_{{\rm{o}}\rho }}} \end{array}} \right] = \mathit{\boldsymbol{Dm}} $

(12)

式中r_E、r_σ、r_ρ分别表示杨氏模量、泊松比和密度的相对变化率。

将式(12)代入式(7)，可以得到

$ \mathit{\boldsymbol{d}} = \mathit{\boldsymbol{WDm}} + \mathit{\boldsymbol{e}} = \mathit{\boldsymbol{Gm}} + \mathit{\boldsymbol{e}} $

(13)

式中W表示具体的子波矩阵。

为了得到稀疏性较强的反演结果，基于基追踪理论建立反演算法^[20-23]。常规的基追踪反演目标函数由L₂模的误差项和L₁模的约束项组成

$ \left\{ {\begin{array}{*{20}{l}} {F\left( \mathit{\boldsymbol{m}} \right) = \left\| {\mathit{\boldsymbol{Gm}} - \mathit{\boldsymbol{d}}} \right\|_2^2 + \lambda {{\left\| \mathit{\boldsymbol{m}} \right\|}_1}}\\ {\mathit{\boldsymbol{G}} = \mathit{\boldsymbol{WD}}} \end{array}} \right. $

(14)

式中λ表示权重系数。

L₁模约束无法提供合理的低频信息^[24]，且逐道的基追踪反演结果缺乏横向连续性，需要结合测井资料建立低频模型对式(14)进行约束，目标函数变为

$ \left\{ \begin{array}{l} F\left( \mathit{\boldsymbol{m}} \right) = \left\| {\mathit{\boldsymbol{Gm}} - \mathit{\boldsymbol{d}}} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{m}} \right\|_1} + \mathit{\boldsymbol{Y}}\\ \mathit{\boldsymbol{Y}} = {\lambda _E}{\left( {\mathit{\boldsymbol{C}}{r_E} - {\mathit{\boldsymbol{\zeta }}_E}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{C}}{r_E} - {\mathit{\boldsymbol{\zeta }}_E}} \right) + {\lambda _\sigma }{\left( {\mathit{\boldsymbol{C}}{r_\sigma} - {\mathit{\boldsymbol{\zeta }}_\sigma }} \right)^{\rm{T}}} \times \\ \;\;\;\left( {\mathit{\boldsymbol{C}}{r_\sigma } - {\mathit{\boldsymbol{\zeta }}_\sigma }} \right) + {\lambda _\rho }{\left( {\mathit{\boldsymbol{C}}{r_\rho } - {\mathit{\boldsymbol{\zeta }}_\rho }} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{C}}{r_\rho } - {\mathit{\boldsymbol{\zeta }}_\rho }} \right) \end{array} \right. $

(15)

式中：积分矩阵${\mathit{\boldsymbol{C}}} = \int_{{t_0}}^t {\rm{d}} \tau $，可以表示为元素为1的下三角矩阵；${\zeta _E} = \frac{1}{2}\ln \frac{{{\mathit{\boldsymbol{E}}}(t)}}{{E\left( {{t_0}} \right)}};{\zeta _\sigma } = \frac{1}{2}\ln |\frac{{{\mathit{\boldsymbol{ \pmb{\mathit{ σ}} }}}(t)}}{{\sigma \left( {{t_0}} \right)}};{\zeta _\rho } = \frac{1}{2}\ln \frac{{{\mathit{\boldsymbol{ \pmb{\mathit{ ρ}} }}}(t)}}{{\rho \left( {{t_0}} \right)}}$。利用谱投影梯度法^[25]即可求解式(15)。

由于低频信息来自实际工区的测井资料，因此最终的反演结果更多包含了工区内的真实背景信息，反演的结果则更加合理，并且模型约束能够有效地改善反演结果横向的连续性。

3 模型试算

为了检验本文反演方法的正确性，根据实际工区测井数据建立了一维模型进行试验。图 1为根据工区测井数据建立的杨氏模量、泊松比和密度模型曲线。将模型曲线代入Zoeppritz方程，利用30Hz雷克子波进行褶积得到合成地震道集，作为输入的地震数据(图 2)，进行叠前地震反演。得到的反演结果如图 3所示，可以看出，模型数据与反演结果基本吻合。为了进一步检验方法的抗噪性，在合成地震道集中加入信噪比为2的高斯白噪声(图 4)，反演结果如图 5所示。可以看出，反演结果的精度较图 3有所下降，但是仍然比较准确，证明反演方法具有一定的抗噪性。

4 应用实例

为了检验本文提出方法的实际应用效果，在中国M油田开展方法的实际应用。工区发育石盒子组含气致密砂岩储层，目的层埋深约1000~1500m。首先利用工区内A井的测井数据计算得到杨氏模量和泊松比曲线，进行岩石物理交会分析(图 6)。可以看出，含气储层具有较高的杨氏模量值和较低的泊松比，对这两种弹性参数进行反演可以为含气储层的预测提供数据支持。

将处理后的地震道集转换到角度域，并进行分角度叠加，得到小角度(1°~13°)叠加道集、中角度道集(14°~26°)和大角度道集(27°~39°)(图 7)。分别利用常规叠前反演方法和本文方法进行叠前地震反演，得到的杨氏模量和泊松比如图 8所示。由图可见，在井旁道的0.89、0.91、0.93和0.96s处有四套储层，显示出杨氏模量高值异常和泊松比低异常。对比可以发现，图 8中两种弹性参数的异常位置与测井曲线(含水饱和度)吻合度更高，显示的储层边界更加清楚，横向展布更加合理。

5 结论

杨氏模量和泊松比可以表征岩石的脆性，对于致密砂岩油气藏的评价和描述非常有利。本文提出了一种基于模型约束的叠前基追踪反演方法，可以从地震数据直接反演得到杨氏模量和泊松比，用来指示致密砂岩油气藏。模型数据测试表明，该方法在低信噪比的情况下仍可以得到高质量的反演结果；实际数据计算结果证明了利用本文方法可以在生产中得到准确、稳定的杨氏模量和泊松比反演结果，为致密砂岩气藏的评价提供可靠的数据支持。