储层流体识别是油气勘探的关键环节, 识别精度关系到油气藏开发的成败。地震数据中包含的诸多信息都被用于指示储层流体性质, 弹性参数因其明确的物理意义和稳定的估计方法成为较常用的一类参数。地下岩石的速度、密度、阻抗等参数由于较容易获得, 因而最早用于判别储层性质。20世纪七八十年代, 泊松比因对含气砂岩储层的良好指示而被广泛应用, 至今仍然在含气储层预测中发挥重要作用^[1]。GOODWAY^[2]提出拉梅参数和密度的乘积是储层流体性质的良好指示参数。GRAY^[3]认为拉梅参数本身就可以区分流体性质, 推导了可以用于拉梅参数直接反演的方程。QUAKENBUSH等^[4]指出在复杂储层条件下可以通过组合纵、横波阻抗得到泊松阻抗来区分储层的岩性和流体性质。RUSSELL等指出流体项ρf和Gassmann流体项f对于流体性质的变化都比较敏感, 而且可以利用地震反演技术对其进行估计^[5-6]。ZONG等^[7]建立了利用纵横波模量表示的孔隙弹性理论和反射系数近似方程, 并发展了基于纵横波模量反演的流体识别方法。郝前勇等^[8]提出基于体积模量和剪切模量的组合弹性模量流体因子可以降低流体识别的多解性。刘立峰等^[9]利用多种地震属性的不同角度叠加数据体组合流体识别因子进行碳酸盐岩储层流体识别。李超等^[10]提出利用孔隙流体模量作为流体因子可以降低固体骨架的影响, 提高流体识别精度。

精确估计流体因子是流体识别过程必不可少的步骤^[11]。叠前地震反演技术极大地促进了地震流体识别技术的发展。BULAND等^[12]和DOWNTON^[13]分别提出了不同的AVO反演方法, 用于直接反演纵横波波阻抗; YOUNG等^[14]采用LMR反演方法将拉梅参数用于储层流体识别; ZONG等^[15]提出了一种叠前同步反演方法用于Gassmann流体项的估计; ZONG等^[16]提出了泊松比和杨氏模量的直接反演方法, 用于页岩储层的脆性评价; LI等^[17]基于入射角近似方程提出了Russell流体项的叠前反演方法; 桂金咏等^[18]提出了致密储层流体因子的叠前同步反演方法。这些常用的叠前反演方法大都利用三项反射系数近似方程, 因此需要从小角度到大角度的叠前道集数据才能获得稳定的反演结果。

Zoeppritz方程的线性近似方程是目前叠前反演中较为常用的反射系数方程。这些近似方程以三项方程为主, 且第3项通常是密度项。要进行密度的稳定反演, 则需要大角度的叠前地震数据^[19-20]。深部地层中蕴含有大量的石油天然气, 但是由于埋深较大, 深层油气藏勘探采集到的地震数据质量远不及中浅层, 最突出的表现就是大角度(大偏移距)的叠前地震道集品质低下, 从而限制了常规的基于叠前地震反演的油气预测技术的应用。

本文提出了一种可以在缺少大角度地震道集数据的情况下精确估计流体因子的叠前地震反演方法。基于Zoeppritz方程的线性近似方程和岩石物理理论推导出一个两项反射系数近似方程, 其包含Gassmann流体项和剪切模量。该近似方程仅有两项, 而且不包含密度项, 因而利用其进行叠前反演不需要大角度叠前地震数据。基于贝叶斯反演理论构建了叠前地震反演的目标方程和相应的反演算法来直接反演流体项, 并利用得到的流体项实现深部储层的流体识别。最后用模型数据和实际地震数据验证了方法的正确性和实用性。

1 流体因子反演方法 1.1 反射系数方程

Zoeppritz反射系数方程可以精确描述反射、透射波的振幅, 但是由于本身的强非线性性, 用其进行反演稳定性较差。因此在叠前地震反演中, 常采用反射系数方程的近似方程。Aki和Richards给出的Zoeppritz方程近似公式如下:

$ \begin{array}{l} R\left( \theta \right) = \frac{1}{2}\left( {1 + {{\tan }^2}\theta } \right)\frac{{\Delta {v_{\rm{P}}}}}{{{v_{\rm{P}}}}} - 4{\gamma ^2}{\sin ^2}\theta \frac{{\Delta {v_{\rm{S}}}}}{{{v_{\rm{S}}}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\left( {1 - 4{\gamma ^2}{{\sin }^2}\theta } \right)\frac{{\Delta \rho }}{\rho } \end{array} $

(1)

其中, Δv_P=v_P2-v_P1, Δv_S=v_S2-v_S1, Δρ=ρ₂-ρ₁, v_P=(v_P1+v_P2)/2, v_S=(v_S1+v_S2)/2, ρ=(ρ₁+ρ₂)/2, γ=v_S/v_P。v_P1, v_S1, v_P2, v_S2, ρ₁, ρ₂分别为界面两侧地层的纵波速度、横波速度和密度, θ为入射角和透射角的平均值。

GARDNER等^[21]研究了地下岩石的纵波速度与密度之间的关系, 将其提出的方程改写为差分形式:

$ \frac{{\Delta \rho }}{\rho } = \frac{1}{c}\frac{{\Delta {v_{\rm{P}}}}}{{{v_{\rm{P}}}}} $

(2)

其中, c为常系数。

根据孔隙弹性介质理论, Gassmann流体项对孔隙流体性质的变化非常敏感, 可以作为储层孔隙流体的指示因子^{[6, 11, 15]}。Gassmann流体项(f)、剪切模量(μ)与纵、横波速度、密度之间存在如下关系:f=(v_P²-γ_dry²v_S²)ρ, μ=v_S²ρ。表示成差分形式为:

$ \frac{{\Delta f}}{f} = \left( {\frac{2}{{1 - \gamma _{{\rm{dry}}}^2{\gamma ^2}}} + \frac{1}{c}} \right)\frac{{\Delta {v_{\rm{P}}}}}{{{v_{\rm{P}}}}} - \frac{{2\gamma _{{\rm{dry}}}^2{\gamma ^2}}}{{1 - \gamma _{{\rm{dry}}}^2{\gamma ^2}}}\frac{{\Delta {v_{\rm{S}}}}}{{{v_{\rm{S}}}}} $

(3)

$ \frac{{\Delta \mu }}{\mu } = 2\frac{{\Delta {v_{\rm{S}}}}}{{{v_{\rm{S}}}}} + \frac{{\Delta \rho }}{\rho } $

(4)

其中, γ_dry为干燥岩石纵横波速度比。

综合方程(2), (3), (4)可以得到:

$ \frac{{\Delta {v_{\rm{P}}}}}{{{v_{\rm{P}}}}} = \frac{{c\gamma _{{\rm{dry}}}^2{\gamma ^2}}}{{2c + 1}}\frac{{\Delta \mu }}{\mu } + \frac{{c\left( {1 - \gamma _{{\rm{dry}}}^2{\gamma ^2}} \right)}}{{2c + 1}}\frac{{\Delta f}}{f} $

(5)

$ \begin{array}{l} \frac{{\Delta {v_{\rm{S}}}}}{{{v_{\rm{S}}}}} = \frac{1}{{2\left( {2c + 1} \right)}}\left[ {\left( {2c + 1 - \gamma _{{\rm{dry}}}^2{\gamma ^2}} \right)\frac{{\Delta \mu }}{\mu } - } \right.\\ \;\;\;\;\;\;\;\;\;\;\left. {\left( {1 - \gamma _{{\rm{dry}}}^2{\gamma ^2}} \right)\frac{{\Delta f}}{f}} \right] \end{array} $

(6)

将方程(2), (5), (6)代入方程(1), 可以得到:

$ \begin{array}{*{20}{c}} {R\left( \theta \right) = \frac{{\left( {1 - \gamma _{{\rm{dry}}}^2{\gamma ^2}} \right)}}{{2\left( {2c + 1} \right)}}\left( {c{{\sec }^2}\theta + 1} \right)\frac{{\Delta f}}{f} + \frac{{{\gamma ^2}}}{{2\left( {2c + 1} \right)}} \cdot }\\ {\left( {c\gamma _{{\rm{dry}}}^2{{\sec }^2}\theta + \gamma _{{\rm{dry}}}^2 - 8c{{\sin }^2}\theta - 4{{\sin }^2}\theta } \right)\frac{{\Delta \mu }}{\mu }} \end{array} $

(7)

方程(7)是一个两项近似方程, 基于此方程建立叠前地震反演方法, 可以在不需要大角度地震数据的情况下反演Gassmann流体项和剪切模量。

需要注意的是, 方程(7)中的参数都需要根据实际工区的具体地质情况确定, 这些参数的准确与否对于方法的精度有很大影响。其中, γ_dry²和γ²需要根据岩石骨架模量、流体模量进行计算, 具体计算公式为:

$ \gamma _{{\rm{dry}}}^2 = \frac{{{K_{{\rm{dry}}}}}}{\mu } + \frac{4}{3} $

(8)

$ {\gamma ^2} = \frac{f}{\mu } + \gamma _{{\rm{dry}}}^2 $

(9)

其中, K_dry为干燥岩石的体积模量。

1.2 方程精度验证

为了检验方程(7)的精度, 选用如表 1所示的两个经典模型(Ostrander模型和Goodway模型), 计算模型中地层界面处的纵波反射系数, 并与Zoeppritz精确方程和Russell近似方程(包含Gassmann流体项的三项近似方程^{[6, 22]})的计算结果进行对比。

图 1a和图 1b分别给出了Ostrander模型和Goodway模型利用不同方程的计算结果。从图 1可以看出, 蓝色线和粉色线均与红色线接近, 而且蓝色线与红色线的距离更近, 说明方程(7)的精度略高于Russell近似方程, 因此其精度满足叠前反演的要求。

1.3 叠前地震反演算法

叠前地震反演是综合利用叠前地震数据、测井数据和地质先验信息的地震反演方法。基于贝叶斯反演理论框架, 建立一种可以稳定估计流体因子的叠前反演方法^[23]。首先考虑1个界面、m个角度的情况, 将方程(7)写成矩阵形式:

$ \left[ {\begin{array}{*{20}{c}} {a\left( {{\theta _1}} \right)}&{b\left( {{\theta _1}} \right)}\\ {a\left( {{\theta _2}} \right)}&{b\left( {{\theta _2}} \right)}\\ \vdots&\vdots \\ {a\left( {{\theta _m}} \right)}&{b\left( {{\theta _m}} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\frac{{\Delta f}}{f}}\\ {\frac{{\Delta \mu }}{\mu }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {R\left( {{\theta _1}} \right)}\\ {R\left( {{\theta _2}} \right)}\\ \vdots \\ {R\left( {{\theta _m}} \right)} \end{array}} \right] $

(10)

其中, a(θ)=[(1-γ_dry²γ²)]/[2(2c+1)](csec²θ+1), b(θ)=γ²/[2(2c+1)](cγ_dry²sec²θ+γ_dry²-8csin²θ-4sin²θ)。

在n个界面、m个角度的情况下, 基于褶积理论添加子波矩阵W, 将矩阵变成块化矩阵方程:

$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{A}}_1}}&{\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{B}}_1}}\\ {\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{A}}_2}}&{\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{B}}_2}}\\ \vdots&\vdots \\ {\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{A}}_m}}&{\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{B}}_m}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{m}}_f}}\\ {{\mathit{\boldsymbol{m}}_\mu }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{d}}_1}}\\ {{\mathit{\boldsymbol{d}}_2}}\\ \vdots \\ {{\mathit{\boldsymbol{d}}_m}} \end{array}} \right] $

(11)

其中, A, B, m, d分别为方程(10)中对应项的矩阵形式。

矩阵方程(11)可以简写成以下形式:

$ \mathit{\boldsymbol{Gm}} = \mathit{\boldsymbol{d}} $

(12)

根据贝叶斯理论, 待反演数据的后验概率密度分布可以表示为:

$ P\left( {\mathit{\boldsymbol{m}}\left| \mathit{\boldsymbol{d}} \right.} \right) = \frac{{P\left( {\mathit{\boldsymbol{d}}\left| \mathit{\boldsymbol{m}} \right.} \right)P\left( \mathit{\boldsymbol{m}} \right)}}{{P\left( \mathit{\boldsymbol{d}} \right)}} $

(13)

其中, P(m|d)是后验概率密度, P(d|m)是似然函数, P(m)是先验概率密度, P(d)是边缘概率密度。

由于观测到的地震数据中包含的噪声通常服从高斯分布, 因此假设似然函数服从高斯分布。为了更好地保护弱小反射, 得到高分辨率的反演结果, 假设先验分布服从两参数柯西分布, 则后验概率密度函数可以表示为:

$ \begin{array}{l} P\left( {\mathit{\boldsymbol{m}}\left| \mathit{\boldsymbol{d}} \right.} \right) = a\exp \left[ {\left( { - \frac{3}{2}} \right)\sum\limits_{i = 1}^n {\ln \left( {1 + {\mathit{\boldsymbol{m}}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^i}\mathit{\boldsymbol{m}}} \right)} } \right] \cdot \\ \;\;\;\;\;\;\;\exp \left[ {\frac{{ - {{\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Gm}}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Gm}}} \right)}}{{2\sigma _{\rm{m}}^2}}} \right] \end{array} $

(14)

式中:σ_m²是噪声的方差; Φⁱ=(Dⁱ)^TΨ^-1Dⁱ, [D_qlⁱ]= $\left\{ \begin{array}{l} 1\;\;q = 1\;{\rm{且}}l = i\\ 1\;\;q = 2\;且l = i + n\\ 0\;\;其它 \end{array} \right.$, Ψ为一个2×2的系数矩阵, 体现各反演参数之间的相关性; a=1/[(2π)^mn/2σ_m^mn]·1/(π²ⁿ|Ψ|^n/2)。

对方程(14)待反演参数求导并取极值, 可以得到如方程(15)的初始目标函数:

$ J\left( \mathit{\boldsymbol{m}} \right) = {\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Gm}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Gm}}} \right) + \frac{3}{2}\mathit{\boldsymbol{Qm}} $

(15)

式中: $Q = \sum\limits_{i = 1}^n {\frac{{2\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{kl}^i}}{{1 + {\mathit{\boldsymbol{m}}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^i}\mathit{\boldsymbol{m}}}}} $。

为了提高反演结果的横向连续性并补偿低频信息, 在目标函数中加入模型约束项, 得到:

$ \begin{array}{l} J\left( \mathit{\boldsymbol{m}} \right) = {\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Gm}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Gm}}} \right) + \frac{3}{2}\mathit{\boldsymbol{Qm}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\alpha {\left\| {\mathit{\boldsymbol{Cm}} - \xi } \right\|^2} \end{array} $

(16)

式中:α是权重因子; C和ξ由低频模型计算得到。

方程(16)的最大后验概率解为:

$ \left( {{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{G}} + \lambda \mathit{\boldsymbol{Q}} + \alpha {\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{C}}} \right)\mathit{\boldsymbol{m}} = {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{d}} + \alpha {\mathit{\boldsymbol{C}}^{\rm{T}}}\xi $

(17)

利用迭代重加权最小二乘方法可以求解方程(17), 然后利用方程(18)对结果积分即可得到最终反演结果。

$ \left\{ \begin{array}{l} f\left( t \right) = f\left( {{t_0}} \right)\exp \left[ {2\int_0^t {{m_f}\left( \tau \right){\rm{d}}\tau } } \right]\\ \mu \left( t \right) = \mu \left( {{t_0}} \right)\exp \left[ {2\int_0^t {{m_\mu }\left( \tau \right){\rm{d}}\tau } } \right] \end{array} \right. $

(18)

式中:t₀为第1个采样点的时间。

两项近似方程的优势不仅在于其不需要大角度地震道集数据, 而且由于方程的项数减少, 其系数矩阵的稳定性也大大提高。计算方程(7)与三项方程(Russell方程)在不同角度情况下系数矩阵的条件数, 来比较两个方程的稳定性(图 2)。图 2中实线和虚线分别是方程(7)和三项方程的条件数, 可以看出前者的条件数更小, 即稳定性更好, 抗噪能力更强。

2 模型试算

为了检验本文反演方法的有效性, 建立典型模型, 利用合成地震记录进行试算。图 3给出了利用实际测井数据建立的模型数据。利用褶积理论和精确Zoeppritz方程可以合成不同角度的叠前地震数据。图 4是无噪声的合成地震记录。对叠前数据进行分角度叠加得到小角度叠加数据(1°~10°)和中角度叠加数据(11°~20°), 利用本文方法进行反演, 结果如图 5所示。图 5中黑色曲线为原始模型, 蓝色曲线为反演结果。从图 5可以看出, 反演得到的Gassmann流体项和剪切模量与模型数据十分接近。

为了进一步说明本文方法的有效性, 利用含噪声的数据测试无大角度数据的叠前地震反演。首先利用前述正演方法得到1°~30°的叠前地震数据, 然后在1°~20°的地震数据中加入信噪比为2:1的高斯噪声, 用于模拟正常包含噪声的中小角度数据。实际生产中采集到的大角度地震数据质量往往较差, 因而在21°~30°的地震数据中加入信噪比为1:1的高斯噪声。加噪声后的地震记录如图 6所示。将1°~20°的地震数据采用本文方法进行反演, 为了进行对比说明, 将1°~30°的地震数据用基于Russell近似方程的贝叶斯反演方法进行了反演。反演结果如图 7a和图 7b所示, 图中黑色曲线是原始模型, 蓝色曲线是本文方法反演结果, 红色曲线是基于Russell近似方程的贝叶斯反演方法的反演结果。从图 7可以看出, 红色曲线与真实模型曲线相差比较大, 说明大角度数据的低信噪比对反演结果的影响较大, 这种情况下依然利用大角度数据进行反演得到的Gassmann流体项不能用于流体识别; 而蓝色曲线误差较小, 说明虽然受噪声的影响, 本文方法反演结果的精度略有下降, 但是仍然在可以接受的范围内, 仅利用中小角度叠前地震数据进行反演也可以得到准确的Gassmann流体项, 从而进行流体识别, 同时也证明本文方法具有较好的稳定性。

3 实际应用

利用东部某油田实际工区地震资料进行应用研究。研究的目标储层位于东营组二段, 埋深约为3400m。由于储层埋深较大, 且受采集条件限制得到的大角度道集信噪比较低, 不适合用于地震反演。为了开展储层特性的研究, 将本文提出的流体识别方法应用于该工区储层流体识别研究。

开展叠前地震反演前各种数据都要进行必要的预处理, 地震数据主要的预处理是保幅处理, 测井数据的主要预处理是曲线校正。之后, 将地震道集转换到角度域, 数据角度范围为3°~23°, 将其分别叠加为中心角为8°的小角度部分角度叠加数据(3°~13°)和中心角为18°的中角度部分角度叠加数据(13°~23°), 如图 8所示。利用工区的测井、地质、岩心资料开展岩石物理研究, 计算出本文方法中所需的γ_dry², γ², c参数值。

在图 8所示地震剖面的CDP136处有一口井(A井), 井位置处有一条测井曲线, 该柱状曲线表示测井解释结果, 红色、绿色和蓝色部分分别表示含油砂岩、含水砂岩和泥岩。利用工区内井数据进行Gassmann流体项和剪切模量的交会分析, 结果如图 9所示。从图 9可以看出, Gassmann流体项可以较好地识别含油砂岩和含水砂岩, 而剪切模量识别效果差于前者。反演前, 综合工区的地震数据、测井数据和地质资料进行了井震标定, 图 10是A井的井震标定结果。

对该工区资料利用本文方法进行叠前地震反演, 得到的Gassmann流体项和剪切模量分别如图 11a和图 11b所示。图 11a中在2.65s处有低值异常, 与井上的含油砂岩对应, 而2.49s处有另一低值异常, 与井上的含水砂岩对应。而且前者Gassmann流体项值明显低于后者, 因此Gassmann流体项可以清晰地区分出含油、含水储层, 与测井解释结果和钻井资料一致。图 11b中在2.49s和2.65s处也有低值异常, 但是二者差别不大, 这是因为孔隙流体的性质对于剪切模量的影响不大。另外, 使用未参与反演过程的B井作为验证井, 过B井的Gassmann流体项和剪切模量剖面分别如图 12a和图 12b所示。图 12a中2.46s和2.56s处的含水砂岩和含油砂岩所对应的低值异常有明显差别, 而且与测井、钻井资料相符合, 说明反演得到的Gassmann流体项可以识别含油砂岩, 与图 11中A井过井剖面结论一致。

4 结论

为了在缺乏大角度叠前地震数据的情况下实现深部储层的流体识别, 本文推导了包含Gassmann流体项的两项地震反射系数近似方程, 并基于此方程构建了贝叶斯叠前地震反演方法, 用于直接、稳定地反演Gassmann流体项。利用该方法可以使用不含大角度的叠前道集数据进行流体项的直接反演, 并开展深部储层的流体识别。模型实验和实际工区应用结果均证明了该方法的有效性、稳定性和实用性。

本文提出的流体识别策略可以扩展到其它类型储层的预测中, 只需要选择合适的弹性参数作为流体因子, 推导包含该弹性参数的两项近似方程, 建立相应的反演方法进行反演和流体识别。另外需要注意的是, 本文方法综合多种理论提出, 涉及的参数较多, 在实际应用中需要根据工区实际情况计算出准确的参数, 如γ_dry², c等, 才可以使方法得到有效应用。