2019,
Vol. 43
2. 中国石油集团东方地球物理勘探有限责任公司, 河北涿州 072750
2. Bureau of Geophysical Prospecting INC, China National Petroleum Corporation, Zhuozhou 072750, China
依据地震数据的流体识别, 即利用地下介质单一弹性参数、组合形式或孔隙流体模量等一系列敏感流体指示参数对孔隙含流体特征进行描述和定量表征, 而叠前地震AVO/AVA(amplitude variation with offset/angle)和EI(elastic impedance)反演是弹性参数预测的有效途径[1]。复杂双相介质情况下, 岩石孔隙流体的存在会引起地震波的衰减和频散, 尤其在含烃储层中传播时, 这种衰减现象更为明显[2]。为充分利用地震资料蕴含的振幅和频率异常信息, 低频阴影[3-4]、衰减特征[5]、频散属性[6]等逐渐被用来进行储层含流体性质检测。国内外学者提出多种孔裂隙含流体岩石物理模型, 试图从机制上明确储层中流体的存在对地震波衰减及频散现象的影响, 同时为开展频变流体识别奠定理论基础[7-9]。Wilson等[10-11]依据Chapman多尺度裂缝频散岩石物理模型, 利用拓展的两项频变Smith-Gidlow近似式, 率先提出了频变AVA反演的概念(F-AVA); 张世鑫等[12-13]将叠后反射系数和Shuey近似公式拓展为频率依赖, 提出纵波速度频散梯度的概念, 研究了基于地震数据的频散属性孔隙流体识别方法。程冰洁等[14]研究了频变AVA属性之间的关系, 将频变AVA流体识别技术应用到低孔低渗碎屑岩储层的含气性识别中。张震等[15]基于Russell近似公式构建了频变流体项f, 研究了基于连续小波多尺度分解的频变流体因子AVA反演方法。李坤等[16]提出了基于匹配追踪高分辨率谱分解算法的时频域F-AVA (Frequency dependent AVA)反演方法, 旨在改善“子波叠印”去除误差和频变属性反演的分辨率。然而, 当前的频散属性反演和流体识别应用仍集中于完全弹性介质假设情况下的叠前AVA反射方面, 大多数是常规AVA非频变反射公式的一种广义频变拓展, 并未真正依据衰减介质的频变AVA反射特征方程展开。笔者利用Futterman常Q衰减模型来考虑地下介质的黏弹性[17-18], 且构建频变黏弹性流体因子(FDV-FI), 推导含频变黏弹性流体因子的F-AVA反射特征方程; 依据贝叶斯理论和先验正则化推导黏弹性F-AVA反演的映射方程和目标泛函; 利用叠前部分角度叠加道集实现基于F-AVA反演的砂岩储层含油气性识别。
1 频变黏弹性流体因子叠前地震F-AVA反演原理 1.1 黏弹性介质频变反射系数考虑介质黏滞性来模拟地震波的衰减和频散效应, 运用Futterman近似常Q模型描述黏弹性介质中纵横波复速度(α, β)随频率的变化特征[17-18], 其中利用相速度和品质因子定量表征的复速度如下:
| $ {\alpha ^{ - 1}} = v_{\rm{p}}^{ - 1}\left( {1 - \frac{1}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}} + \frac{{\rm i}}{{{\rm{2}}{Q_{\rm{p}}}}}} \right), $ | (1) |
| $ {\beta ^{ - 1}} = v_{\rm{s}}^{ - 1}\left( {1 - \frac{1}{{{\rm{ \mathsf{ π} }}{Q_{\rm{s}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}} + \frac{{\rm{i}}}{{{\rm{2}}{Q_{\rm{s}}}}}} \right). $ | (2) |
式中, Qp和Qs分别为纵波品质因子和横波品质因子; vp和vs分别为参考频率ωr对应的纵横波相速度。黏弹性介质的平面波散射方程与完全弹性介质的反射形式一致, 差异在于5个波(入射P波、反射P波和SV波、透射P波和SV波)的传播速度、波数、振幅均为复数, 使非弹性介质的平面波传播特征异于弹性介质[18], 非弹性介质Aki-Richard纵波反射系数为
| $ \begin{array}{l} {m_{{\rm{pp}}}}\left( {\theta ,\omega } \right) = a\left( \theta \right)\left[ {\Delta \alpha /\alpha \left( \omega \right)} \right] + b\left( \theta \right)\left[ {\Delta \beta /\beta \left( \omega \right)} \right] + \\ c\left( \theta \right)\left[ {\Delta \rho /\rho } \right]. \end{array} $ | (3) |
其中
| $ a\left( \theta \right) = \left( {1 + {{\tan }^2}\theta } \right), $ |
| $ b\left( \theta \right) = - 8\gamma _{{\rm{sat}}}^{ - 2}{\sin ^2}\theta , $ |
| $ c\left( \theta \right) = 1 - 4\gamma _{{\rm{sat}}}^{ - 2}{\sin ^2}\theta , $ |
| $ \gamma _{{\rm{sat}}}^{ - 2} = {\beta ^2}{\alpha ^{ - 2}}. $ |
式中, γsat为饱和岩石横纵波速度平方比。令黏弹性流体因子fvis=ρ[α2-γdry2β2], μvis=ρβ2, ρvis=ρ, 假设干岩石骨架纵横波速度比γdry不受衰减影响且Qs≫Qp的情况:
| $ \begin{array}{l} {m_{{\rm{pp}}}}\left( {\theta ,\omega } \right) = \frac{1}{4}\left( {1 - \frac{{\gamma _{{\rm{dry}}}^2}}{{\gamma _{{\rm{sat}}}^2}}} \right){\rm{se}}{{\rm{c}}^2}\theta \frac{{\Delta {f_{{\rm{vis}}}}}}{{{f_{{\rm{vis}}}}}} + \\ \left( {\frac{{\gamma _{{\rm{dry}}}^2}}{{4\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta - \frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta } \right)\frac{{\Delta {\mu _{{\rm{vis}}}}}}{{{\mu _{{\rm{vis}}}}}} + \left( {\frac{1}{2} - \frac{{{{\sec }^2}\theta }}{4}} \right)\frac{{\Delta \rho }}{\rho }. \end{array} $ | (4) |
其中
| $ \begin{array}{l} \gamma _{{\rm{sat}}}^{ - 2} = {\beta ^2}{\alpha ^{ - 2}} = \frac{{v_{\rm{s}}^2}}{{v_{\rm{p}}^2}}{\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}} + \frac{{\rm{i}}}{{2{Q_{\rm{p}}}}}} \right)^2} \approx \\ \gamma {'}_{{\rm{sat}}}^{ - 2}\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}} + \frac{{\rm{i}}}{{{Q_{\rm{p}}}}}} \right), \end{array} $ | (5) |
| $ \begin{array}{l} {f_{{\rm{vis}}}} = \rho \left( {{\alpha ^2} - \gamma _{{\rm{dry}}}^2{\beta ^2}} \right) \approx \rho v_{\rm{p}}^2 - \gamma _{{\rm{dry}}}^2\rho v_{\rm{s}}^2 + \\ \rho v_{\rm{p}}^2\left( {\frac{2}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}} - \frac{{\rm{i}}}{{{Q_{\rm{p}}}}}} \right). \end{array} $ | (6) |
将式(5)、(6)带入式(4)中, 可得
| $ \begin{array}{l} {m_{{\rm{pp}}}}\left( {\theta ,\omega } \right) = \frac{1}{4}\left[ {1 - \frac{{\gamma _{{\rm{dry}}}^2}}{{\gamma {'}_{{\rm{sat}}}^2}}\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}} + } \right.} \right.\\ \left. {\left. {\frac{{\rm{i}}}{{{Q_{\rm{p}}}}}} \right)} \right]{\sec ^2}\theta \frac{{\Delta {f_{{\rm{vis}}}}}}{{{f_{{\rm{vis}}}}}} + \left[ {\frac{{\gamma _{{\rm{dry}}}^2}}{{4\gamma {'}_{{\rm{sat}}}^2}}\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}} + \frac{{\rm{i}}}{{{Q_{\rm{p}}}}}} \right){{\sec }^2}\theta - } \right.\\ \left. {\frac{2}{{\gamma {'}_{{\rm{sat}}}^2}}\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}} + \frac{{\rm{i}}}{{{Q_{\rm{p}}}}}} \right){{\sin }^2}\theta } \right]\frac{{\Delta {\mu _{{\rm{vis}}}}}}{{{\mu _{{\rm{vis}}}}}} + \\ \left( {\frac{1}{2} - \frac{{{{\sec }^2}\theta }}{4}} \right)\frac{{\Delta \rho }}{\rho }. \end{array} $ | (7) |
由于黏弹性纵波反射mpp(θ, ω)含虚部项, 假设弱非弹性情况, 则虚部项为相对小量, 将其忽略得实反射系数为
| $ \begin{array}{l} {m_{{\rm{pp}}}}\left( {\theta ,\omega } \right) = \\ \frac{1}{4}\left[ {1 - \frac{{\gamma _{{\rm{dry}}}^2}}{{\gamma {'}_{{\rm{sat}}}^2}}\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}}} \right)} \right]{\sec ^2}\theta \frac{{\Delta {f_{{\rm{vis}}}}}}{{{f_{{\rm{vis}}}}}} + \\ \frac{1}{{\gamma {'}_{{\rm{sat}}}^2}}\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{\rm{p}}}}}\log \frac{\omega }{{{\omega _{\rm{r}}}}}} \right)\left( {\frac{{\gamma _{{\rm{dry}}}^2}}{4}{{\sec }^2}\theta - 2{{\sin }^2}\theta } \right)\frac{{\Delta {\mu _{{\rm{vis}}}}}}{{{\mu _{{\rm{vis}}}}}} + \\ \left( {\frac{1}{2} - \frac{{{{\sec }^2}\theta }}{4}} \right)\frac{{\Delta \rho }}{\rho }. \end{array} $ | (8) |
式(8)即基于Futterman黏弹模型构建的黏弹性AVAF线性近似方程。为验证该黏弹性反射特征方程的精确度, 本文中采用双层模型计算单界面黏弹性AVAF反射系数, 模型参数见表 1。基于Futterman黏弹模型的精确Zoeppritz方程、常规Russell近似方程和本文推导的黏弹AVAF近似方程计算得到的反射系数如图 1和2所示。
|
表 1 黏弹性介质频变双层模型 Table 1 Frequency-dependent two-layer model for viscoelastic media |
|
图 1 参考频率fr=20 Hz处对应的AVA反射系数精度分析 Fig.1 Precision analysis of AVA reflection coefficient at reference frequency fr=20 Hz |
|
图 2 黏弹性AVA近似方程的频变反射特征与精度分析 Fig.2 Frequency-dependentcharacteristics and accuracy analysis of viscoelastic AVA approximation equation |
图 1(a)为参考频率fr=
黏弹性介质中叠前AVA反射系数随着频率的变化发生变化, 如图 2(a)、(b)所示(其中, mZoepp(θ, ω)表示黏弹性介质精确Zoeppritz方程), 当频率ω≠ωr时, 黏弹性反射的实部和虚部都是频率依赖的, 然而虚部的反射振幅相比实部AVA反射振幅较小, 因此本文中忽略弱小虚部反射。图 2(c)、(d)所示为黏弹性AVO近似方程的频变反射特征与误差分析。由图 2可见, 本文中构建的黏弹性频变AVA反射系数(式(8))与黏弹Zoeppritz方程的实部保持较好的一致性, 进而说明该近似方程在黏弹性频变流体因子提取方面的可行性。
1.2 黏弹性频变流体因子反演目前, 利用叠前地震振幅和频率信息的频变AVA反演均是建立在完全弹性介质假设情况下进行的[10-16], 即平面谐纵波入射到两半无限弹性介质的Zoeppritz近似公式的的频变AVA拓展应用。本文中在黏弹性AVA近似方程(式(8))推导的基础上, 提出了基于F-AVA反演的黏弹性频变流体因子提取方法。针对式(8)求取频率ω的偏导数得
| $ \begin{array}{l} \frac{{\partial {m_{{\rm{pp}}}}\left( {\theta ,\omega } \right)}}{{\partial \omega }} = A\left( {\theta ,\omega } \right)\frac{\partial }{{\partial \omega }}\left( {\frac{{\Delta {f_{{\rm{vis}}}}}}{{{f_{{\rm{vis}}}}}}} \right) + B\left( {\theta ,\omega } \right)\frac{\partial }{{\partial \omega }}\left( {\frac{{\Delta {\mu _{{\rm{vis}}}}}}{{{\mu _{{\rm{vis}}}}}}} \right) + \\ \frac{{\partial A\left( {\theta ,\omega } \right)}}{{\partial \omega }}\frac{{\Delta {f_{{\rm{vis}}}}}}{{{f_{{\rm{vis}}}}}} + \frac{{\partial B\left( {\theta ,\omega } \right)}}{{\partial \omega }}\frac{{\Delta {\mu _{{\rm{vis}}}}}}{{{\mu _{{\rm{vis}}}}}}. \end{array} $ | (9) |
式中, A(θ, ω)和B(θ, ω)为黏弹AVA近似方程的频变系数。令黏弹性反射系数的频率偏导数Rωθ=
| $ \begin{array}{l} R_\omega ^\theta \left| {_{{\omega ^ * }}} \right. = \left[ {A\left( {\theta ,\omega } \right){F_\omega } + B\left( {\theta ,\omega } \right){U_\omega } + A_\omega ^\theta R{f^\omega } + } \right.\\ {\left. {B_\omega ^\theta R{\mu ^\omega }} \right]_{{\omega ^ * }}}. \end{array} $ | (10) |
式中, Rεθ|ω*为频率ω*位置处的反射系数频率偏导数。为了建立频变黏弹性反射系数与叠前地震数据间的直接关系, 令式(10)与主频为ω*的地震子波褶积, 可得
| $ \begin{array}{l} W\left( {\theta ,{\omega ^ * }} \right) * R_\omega ^\theta \left| {_{{\omega ^ * }}} \right. = W\left( {\theta ,{\omega ^ * }} \right) * \left[ {A\left( {\theta ,\omega } \right){F_\omega } + } \right.\\ {\left. {B\left( {\theta ,\omega } \right){U_\omega } + A_\omega ^\theta R{f^\omega } + B_\omega ^\theta R{\mu ^\omega }} \right]_{{\omega ^ * }}}. \end{array} $ | (11) |
令偏导数Rωθ在频率ω*时的偏导数近似表示为
| $ R_\omega ^\theta \left| {_{{\omega ^ * }}} \right. = \frac{{{m_{{\rm{pp}}}}\left( {\theta ,{\omega _{i + 1}}} \right) - {m_{{\rm{pp}}}}\left( {\theta ,{\omega _i}} \right)}}{{{\omega _{i + 1}} - {\omega _i}}}. $ | (12) |
当Δω=ωi+1-ωi在较小范围波动, 假设Sθ, ωi+1≈W(θ, ω*)*mpp(θ, ωi+1), Sθ, ωi≈W(θ, ω*)*mpp(θ, ωi), 则式(11)可以重新组织为
| $ \begin{array}{l} {S^{\theta ,{\omega _{i + 1}}}} - {S^{\theta ,{\omega _i}}} \approx \Delta {\omega _i}W\left( {\theta ,{\omega ^ * }} \right) * \left[ {A\left( {\theta ,{\omega ^ * }} \right){F_{{\omega ^ * }}} + } \right.\\ \left. {B\left( {\theta ,{\omega ^ * }} \right){U_{{\omega ^ * }}} + A_{{\omega ^ * }}^\theta R{f^{{\omega ^ * }}} + B_{{\omega ^ * }}^\theta R{\mu ^{{\omega ^ * }}}} \right]. \end{array} $ | (13) |
式中, Fω*为在频率ω*处的频变黏弹性流体因子; Uω*为在频率ω*处的频变黏弹性剪切模量; Rfω*为频率ω*处的黏弹流体因子反射率; Rμω*为频率ω*处的黏弹剪切模量反射率; A(θ, ω*)和B(θ, ω*)为频率ω*处的常系数; Aω*θ和Bω*θ为频率ω*处的频变系数。令L为入射角个数, M为选择频率对[ωi+1, ωi]的个数。考虑式(13)中系数Aω*θ和Bω*θ所引起的频变AVA效应相比黏弹性低Q频散效应要小得多, 故式(13)在含噪音N干扰情况下的矩阵形式可简化为
| $ \begin{array}{l} \underbrace {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{S}}^{{\theta _1},{\omega _2}}} - {\mathit{\boldsymbol{S}}^{{\theta _1},{\omega _1}}}}\\ {{\mathit{\boldsymbol{S}}^{{\theta _1},{\omega _4}}} - {\mathit{\boldsymbol{S}}^{{\theta _1},{\omega _3}}}}\\ \vdots \\ {{\mathit{\boldsymbol{S}}^{{\theta _1},{\omega _{2M}}}} - {\mathit{\boldsymbol{S}}^{{\theta _1},{\omega _{2M - 1}}}}}\\ \vdots \\ {{\mathit{\boldsymbol{S}}^{{\theta _L},{\omega _{2M}}}} - {\mathit{\boldsymbol{S}}^{{\theta _L},{\omega _{2M - 1}}}}} \end{array}} \right]}_{\Delta \mathit{\boldsymbol{S}}} = \\ \underbrace {\left[ {\begin{array}{*{20}{c}} {\Delta {\omega _1}\mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\mathit{\boldsymbol{A}}\left( {{\theta _1},{\omega ^ * }} \right)}&{\Delta {\omega _1}\mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\mathit{\boldsymbol{B}}\left( {{\theta _1},{\omega ^ * }} \right)}\\ {\Delta {\omega _2}\mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\mathit{\boldsymbol{A}}\left( {{\theta _1},{\omega ^ * }} \right)}&{\Delta {\omega _2}\mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\mathit{\boldsymbol{B}}\left( {{\theta _1},{\omega ^ * }} \right)}\\ \vdots&\vdots \\ {\Delta {\omega _M}\mathit{\boldsymbol{W}}\left( {{\theta _L}} \right)\mathit{\boldsymbol{A}}\left( {{\theta _L},{\omega ^ * }} \right)}&{\Delta {\omega _M}\mathit{\boldsymbol{W}}\left( {{\theta _L}} \right)\mathit{\boldsymbol{B}}\left( {{\theta _L},{\omega ^ * }} \right)} \end{array}} \right]}_{{{\mathit{\boldsymbol{\tilde G}}}_F}} \cdot \\ \underbrace {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{F}}_{{\omega ^ * }}}}\\ {{\mathit{\boldsymbol{U}}_{{\omega ^ * }}}} \end{array}} \right]}_\mathit{\boldsymbol{R}} + \mathit{\boldsymbol{N}}. \end{array} $ | (14) |
其中
| $ \mathit{\boldsymbol{A}}\left( {{\theta _i},{\omega ^ * }} \right) = {\rm{diag}}\left[ {\mathit{\boldsymbol{A}}\left( {{\theta _i},{\omega ^ * }} \right)} \right], $ |
| $ \mathit{\boldsymbol{B}}\left( {{\theta _i},{\omega ^ * }} \right) = {\rm{diag}}\left[ {\mathit{\boldsymbol{B}}\left( {{\theta _i},{\omega ^ * }} \right)} \right]. $ |
式中, Δ
| $ \begin{array}{l} p\left( {\mathit{\boldsymbol{R}}\left| {\Delta \tilde S} \right.} \right) = \frac{{{p_{{\rm{Cauchy}}}}\left( \mathit{\boldsymbol{R}} \right){p_{{\rm{Gauss}}}}\left( {\Delta \mathit{\boldsymbol{\tilde S}}\left| \mathit{\boldsymbol{R}} \right.} \right)}}{{\int {{p_{{\rm{Cauchy}}}}\left( \mathit{\boldsymbol{R}} \right){p_{{\rm{Gauss}}}}\left( {\Delta \mathit{\boldsymbol{\tilde S}}\left| \mathit{\boldsymbol{R}} \right.} \right)d\left( \mathit{\boldsymbol{R}} \right)} }} \propto \\ {p_{{\rm{Cauchy}}}}\left( \mathit{\boldsymbol{R}} \right){p_{{\rm{Gauss}}}}\left( {\Delta \mathit{\boldsymbol{\tilde S}}\left| \mathit{\boldsymbol{R}} \right.} \right). \end{array} $ | (15) |
为了提高频变黏弹性流体因子反演算法的稳定性, 本文中利用常规L2范数正则化方法来引入低频模型约束||ηj-PjRj||22来缓解该不适定性问题, 后验概率分布的最大值与如下目标泛函F(R)等价[19-20]:
| $ \begin{array}{l} F\left( \mathit{\boldsymbol{R}} \right) = L_2^2\left[ \mathit{\boldsymbol{N}} \right] + {F_{{\rm{pri}}}}\left[ \mathit{\boldsymbol{R}} \right] = \sigma _n^{ - 2}\left\| {\Delta \mathit{\boldsymbol{\tilde S}} - {{\mathit{\boldsymbol{\tilde G}}}_{\rm{F}}} \cdot \mathit{\boldsymbol{R}}} \right\|_2^2 + \\ 2\sum\limits_{i = 1}^N {\ln \left[ {1 + r_i^2/\sigma _{\rm{r}}^2} \right]} + \sum\limits_{j = 1}^2 {{\lambda _j}\left\| {{\mathit{\boldsymbol{\eta }}_j} - {\mathit{\boldsymbol{P}}_j}{\mathit{\boldsymbol{R}}^j}} \right\|_2^2} . \end{array} $ | (16) |
式中, L22[N]为不同角度道集间频率差异信息的拟合度; Fpri[R]为稀疏正则项和低频模型约束项; σn2为叠前角道集间差异数据的方差; η1=[Fω*]pri为频变黏弹性流体因子Fω*的背景先验; η2=[Uω*]pri为频变黏弹性剪切模量Uω*的背景先验; Pj为正则化对角阵。通过最小化目标泛函F(R), 即得最大化后验概率密度解R*(MAP解):
| $ \begin{array}{l} {\mathit{\boldsymbol{R}}^ * } = {\left[ {\sigma _n^{ - 2}\mathit{\boldsymbol{\tilde G}}_{\rm{F}}^{\rm{T}}{{\mathit{\boldsymbol{\tilde G}}}_{\rm{F}}} + \sum\limits_{j = 1}^2 {{\lambda _j}\mathit{\boldsymbol{P}}_j^{\rm{T}}{\mathit{\boldsymbol{P}}_j}} + {\varepsilon ^2}\mathit{\boldsymbol{Q}}} \right]^{ - 1}} \times \\ \left[ {\sigma _n^{ - 2}\mathit{\boldsymbol{\tilde G}}_{\rm{F}}^{\rm{T}} \cdot \mathit{\boldsymbol{ \boldsymbol{\varDelta} \tilde S}} + \sum\limits_{i = 1}^2 {{\lambda _i}\mathit{\boldsymbol{P}}_i^{\rm{T}}{\mathit{\boldsymbol{\eta }}_i}} } \right]. \end{array} $ | (17) |
其中
| $ \mathit{\boldsymbol{Q}} = {\rm{diag}}\left\{ {{{\left[ {1 + r_i^2/\sigma _{\rm{r}}^2} \right]}^{ - 1}}} \right\}. $ |
式中, Q为Cauchy概率密度分布引入的稀疏正则化矩阵; σr2为待反演参数R的方差; λi为不同先验模型在黏弹性F-AVA叠前反演中所占权重; ε2为稀疏约束系数。当叠前地震数据信噪比(RSN)较低时, 应适当提高ε2和λi正则化系数有助于提高反演结果的稳定性。针对式(17)的弱非线性问题, 本文中采用反复重加权最小二乘算法(IRLS)计算最优解, 其算法如下:
(1) 设置初始解R0, 迭代阈值κ2, 迭代次数K, 数据协方差σn2, 参数协方差σr2, 模型系数λj, 稀疏约束系数ε2。
(2) 输入R0、σr2, 计算初始稀疏正则化对角阵Q0=diag{[1+ri(0)2/σr2]-1}。
(3) 采用共轭梯度算法求解方程(16)的最优解Rk-1。
(4) 利用Rk-1、σr2, 更新稀疏正则化对角阵Qk-1=diag{[1+ri(k-1)2/σr2]-1}。
(5) 判断迭代终止条件, 迭代残差
(6) 输出频变黏弹性属性的最优模型参数Fω*、Uω*。
2 理论数据分析为了验证频变黏弹性流体因子叠前地震F-AVA反演方法的可行性和抗噪声能力, 本文中利用重采样后的实测弹性参数、30 Hz零相位Ricker子波及黏弹性AVA近似方程式(8)合成叠前地震数据, 图 3所示为理论流体因子f、剪切模量μ, 密度ρ和纵波品质因子Qp(注:横波品质因子Qs＞100)。图 4所示为含不同噪声情况下的合成地震数据和黏弹性频散属性反演结果(无噪声、RSN=5、RSN=2及RSN=1), 其中黑色为依据Futterman近似常Q模型计算的频变黏弹性流体因子和频变剪切模量, 绿色虚线为初始模型, 红色虚线为黏弹性频散属性Fω*和Uω*的反演结果。
|
图 3 理论流体因子f、剪切模量μ, 密度ρ和纵波品质因子Qp模型 Fig.3 Theoretical fluid factor f, shear modulus μ, density ρ and P-wave quality factor Qp |
|
图 4 不同信噪比RSN情况下的合成地震数据和黏弹性频散属性反演测试 Fig.4 Synthetic seismic data and estimated results of viscoelastic dispersion attributes with different signal to noise ratio RSN |
从图 4(a)、(b)中可见, 理论叠前地震反射振幅不仅随着入射角度的变化发生变化, 且在黏弹性介质的假设下, 其同样随着频率的变化发生明显变化, 故此为黏弹性频散属性的F-AVA反演预测提供了有效的数据基础。同时, 依据本文中提出的黏弹性频散属性F-AVA反演算法预测得到的频变黏弹性流体因子Fω*与理论数值保持较高的一致性, 即使在RSN=2和RSN=1噪声情况下也可以实现理论频散模型的高吻合度评价。然而, 对比本算法中频变黏弹性参数Fω*和Uω*的反演精度可以看出, 频变黏弹性剪切模量Uω*的反演精度比频变黏弹性流体因子Fω*的预测精度相对较低, 尤其当地震数据在受到强噪声(RSN=2和RSN=1)干扰时, 频变黏弹性剪切模量Uω*的预测精度将大大降低, 究其原因是由于横波品质因子Qs＞100, 其横波的衰减和频散程度相对较小。本文中主要关注频变黏弹性流体因子Fω*在含油气储层流体识别中的应用, 可靠的频变黏弹性流体因子Fω*预测结果验证了该黏弹性F-AVA反演算法的有效性和抗噪能力。
3 实际资料处理为了验证该频变黏弹性流体因子F-AVA反演方法的实用性, 在理论模型验证该方法可行性的基础上, 针对中国东部K勘探工区的实际地震资料进行了处理。如图 5所示为过LM井的的部分角度叠加实际地震剖面, 其中入射角度依次为12°、21°和30°, 采样间隔为2 ms, 剖面中黑色曲线所示为电阻率Rs测井曲线(浅双侧向电阻率测井结果在油藏发育位置为异常高值, 如图中黑色箭头指示位置), H1和H2为研究区目的层。
|
图 5 中国东部K研究区内的部分角度叠加地震道集 Fig.5 Partial-angle stack seismic gathersfrom K explorationarea in east of China |
首先通过对LM井旁道数据频谱分析可知, 高信噪比有效频带范围为10.2~63.5Hz。本文中选取f0=31 Hz为最优参考频率, 利用MPQ-WVD(快速匹配追踪谱分解方法)对中角度地震数据进行谱分解, 提取f=15、25、40 Hz的单频瞬时谱分量, 如图 6(a)~(c)所示。由图 6中可见, 15Hz单频谱中在油藏下方(白色椭圆区域)存在明显的“低频阴影”异常现象, 随着频率分量的不断升高, 油藏下方低频强能量团逐渐消失, “低频阴影”现象的产生为后续频变黏弹性流体因子的反演提供佐证, 但单频剖面中存在很多“虚假亮点”, 故仅依靠“低频阴影”对孔隙流体进行判识将会产生干扰。
|
图 6 中角度地震数据的匹配追踪单频谱剖面 Fig.6 Iso-frequency instantaneous spectrums of middle-angle seismic gathers based on matching pursuitspectral decomposition |
频变黏弹性流体因子反演综合利用了叠前反射数据的频率和振幅异常响应, 图 7(a)~(c)依次为频变黏弹性流体因子、频变黏弹性剪切模量、常规频变流体因子[16]及归一化高频瞬时谱衰减梯度的反演结果。其中, 高频瞬时谱衰减梯度是利用中角度地震数据(图 5(b))的单频瞬时谱f=56、40、31 Hz计算得到, 即求取单频振幅谱的高频能量随频率的衰减速度。由图 7中可见, 归一化的频变黏弹性流体因子(图 7(a))、常规频变流体因子(图 7(c))和高频瞬时谱衰减梯度属性(图 7(d))在油藏区域(如箭头所示)均呈现异常高值, 与实测电阻率解释结果吻合较好, 且反演剖面的频变异常不明显, 验证了该方法的有效性。然而, 高频瞬时谱衰减梯度和常规频变流体因子的反演剖面(图 7(c))中却存在较多“虚假亮点”干扰, 尤其是高频瞬时谱衰减梯度(图 7(d))相比频变黏弹性流体因子FAVA地震反演的纵向分辨率较低, 在精确识别油气藏发育位置方面存在挑战, 在保证相同叠前反演算法的情况下, 黏弹性流体因子fvis的频散程度Fω*的油层指示敏感性将大幅度提高(如图 7(a)和(c)中黑色椭圆和箭头所示), 进而说明了频变黏弹属性相比常规频散属性能够更好地表征储层孔隙内的含烃性质。
|
图 7 频变黏弹性流体因子和常规频散衰减属性的实际应用效果分析 Fig.7 Analysis of practical application effect of frequency dependent viscoelastic fluid factor and conventional dispersion-attenuation attributes |
通过在频变AVO/AVA叠前地震反演中引入Futterman近似常Q模型, 提出了基于贝叶斯F-AVA反演算法的频变黏弹性流体因子预测方法。频变黏弹性流体因子叠前反演是常规频变AVA反演的一种拓展, 其综合考虑了地下介质黏滞性对地震波AVA反射特征的影响, 且黏弹性频率依赖的AVA近似为该F-AVA反演的关键理论基础。通过理论模型将本文中推导的黏弹性近似方程与基于Futterman常Q模型的精确Zoeppritz方程进行了对比, 验证了该F-AVA近似公式具有较高的精度。关于反演算法方面, 联合贝叶斯公式和先验模型正则化项构建了反演映射矩阵和目标泛函, 模型测试说明了该方法在提取黏弹性频散属性方面的可行性和抗噪性。最后本文中将该F-AVA反演方法应用于实际地震数据的处理和储层流体识别测试, 通过与匹配追踪单频地震剖面、常规频变流体因子及高频瞬时谱衰减梯度属性的应用效果对比发现, 该频变黏弹性流体因子的反演结果中具有更少的“虚假亮点”干扰, 在储层孔隙烃类检测方面将具有更高的预测精度。
| [1] |
印兴耀, 曹丹平, 王保丽, 等. 基于叠前地震反演的流体识别方法研究进展[J]. 石油地球物理勘探, 2014, 49(1): 22-34. YIN Xingyao, CAO Danping, WANG Baoli, et al. Research progress of fluid discrimination with prestack seismic inversion[J]. Oil Geophysical Prospecting, 2014, 49(1): 22-34. |
| [2] |
YIN Xingyao, ZONG Zhaoyun, WU Guochen. Research on seismic fluid identification driven by rock physics[J]. Science China: Earth Sciences, 2015, 58(2): 159-171. DOI:10.1007/s11430-014-4992-3 |
| [3] |
TANER M T, KOEHLER F, SHERIFF R E. Complex seismic trace analysis[J]. Geophysics, 1979, 44(6): 1041-1063. DOI:10.1190/1.1440994 |
| [4] |
CASTAGNA J P, SUN S, SIEGFRIED R W. Instantaneous spectral analysis: detection of low frequency shadows associated with hydrocarbons[J]. The Leading Edge, 2003, 22(2): 120-127. DOI:10.1190/1.1559038 |
| [5] |
LI F, VERMA S, ZHOU H, et al. Seismic attenuation attributes with applications on conventional and unconventional reservoirs[J]. Interpretation, 2016, 4(1): SB63-SB77. DOI:10.1190/INT-2015-0105.1 |
| [6] |
郝前勇, 张世鑫, 张峰, 等. 基于频变AVO反演的频散属性估算方法及其应用[J]. 石油地球物理勘探, 2013, 48(2): 255-261. HAO Qianyong, ZHANG Shixin, ZHANG Feng, et al. Dispersion attributes estimation based on frequency-dependent AVO inversion and its application in hydrocarbon detection[J]. Oil Geophysical Prospecting, 2013, 48(2): 255-261. |
| [7] |
印兴耀, 刘倩. 一致密储层各向异性地震岩石物理建模及应用[J]. 中国石油大学学报(自然科学版), 2016, 40(2): 52-58. YIN Xingyao, LIU Qian. Anisotropic rock physics modeling of tight sandstone and applications[J]. Journal of China University of Petroleun (Edition of Natural Science), 2016, 40(2): 52-58. DOI:10.3969/j.issn.1673-5005.2016.02.006 |
| [8] |
CHAPMAN M, ZATSEPIN S V, CRAMPIN S. Derivation of a microstructural poroelastic model[J]. Geophysical Journal International, 2002, 151(2): 427-451. DOI:10.1046/j.1365-246X.2002.01769.x |
| [9] |
CHAPMAN M. Frequency-dependent anisotropy due to meso-scale fractures in the presence of equant porosity[J]. Geophysical Prospecting, 2003, 51(5): 369-379. DOI:10.1046/j.1365-2478.2003.00384.x |
| [10] |
WILSON A, CHAPMAN M, LI Xiangyang. Frequency-dependent AVO inversion[C/OL]. 2009 SEG Annual Meeting, Houston City, October 341-345, 2009: SEG Technical Program Expanded Abstracts[2016-11-16]. https://library.seg.org/doi/abs/10.1190/1.3255572.
|
| [11] |
WU X, CHAPMAN M, WILSON A, et al. Estimating seismic dispersion from pre-stack data using frequency-dependent AVO inversion[C/OL]. 2010 SEG Annual Meeting, Denver City, October 425-429, 2010: SEG Technical Program Expanded Abstracts[2016-12-08]. https://library.seg.org/doi/abs/10.1190/1.3513759.
|
| [12] |
张世鑫, 印兴耀, 张广智, 等. 纵波速度频散属性反演方法研究[J]. 石油物探, 2011, 50(3): 219-224. ZHANG Shixin, YIN Xingyao, ZHANG Guangzhi, et al. Inversion method for the velocity dispersion-dependent attribute of P-wave[J]. Geophysical Prospectingfor Petroleum, 2011, 50(3): 219-224. DOI:10.3969/j.issn.1000-1441.2011.03.002 |
| [13] |
ZHANG Shixin, YIN Xingyao, ZHANG Guangzhi. Dispersion-dependent attribute and application in hydrocarbon detection[J]. Journal of Geophysics and Engineering, 2011, 8(4): 498-507. DOI:10.1088/1742-2132/8/4/002 |
| [14] |
程冰洁, 徐天吉, 李曙光. 频变AVO含气性识别技术研究与应用[J]. 地球物理学报, 2012, 55(2): 608-613. CHENG Bingjie, XU Tianji, LI Shuguang. Research and application of frequency dependent AVO analysis for gas recognition[J]. Chinese Journal of Geophysics, 2012, 55(2): 608-613. |
| [15] |
张震, 印兴耀, 郝前勇. 基于AVO反演的频变流体识别方法研究[J]. 地球物理学报, 2014, 57(12): 4171-4184. ZHANG Zhen, YIN Xingyao, HAO Qianyong. Frequency-dependent fluid identification method based on AVO inversion[J]. Chinese Journal of Geophysics, 2014, 57(12): 4171-4184. DOI:10.6038/cjg20141228 |
| [16] |
李坤, 印兴耀, 宗兆云. 基于匹配追踪谱分解的时频域FAVO流体识别方法[J]. 石油学报, 2016, 37(6): 777-786. LI Kun, YIN Xingyao, ZONG Zhaoyun. Time frequency domain FAVO fluid discrimination method based on matching pursuit spectrum decomposition[J]. Acta Petrolei Sinica, 2016, 37(6): 777-786. |
| [17] |
FUTTERMAN W I. Dispersive body waves[J]. Journal of Geophysical Research, 1962, 67(13): 5279-5291. DOI:10.1029/JZ067i013p05279 |
| [18] |
郭俊超.叠前地震"3+Q"反演研究[D].青岛: 中国石油大学(华东), 2014. GUO Junchao. Research on pre-stack seismic inversion for "3D"parameters[D]. Qingdao: China University of Petroleum(East China), 2014. http://cdmd.cnki.com.cn/Article/CDMD-10425-1016711699.htm |
| [19] |
李坤, 印兴耀, 宗兆云. 利用平滑模型约束的频率域多尺度地震反演[J]. 石油地球物理勘探, 2016, 51(4): 760-768. LI Kun, YIN Xingyao, ZONG Zhaoyun. Seismic multi-scale inversion in the frequency domain based on smooth model constraint[J]. Oil Geophysical Prospecting, 2016, 51(4): 760-768. |
| [20] |
李坤, 印兴耀, 宗兆云, 等. 基于快速匹配追踪的混合域地震稀疏反演方法[J]. 中国石油大学学报(自然科学版), 2018, 42(1): 50-59. LI Kun, YIN Xingyao, ZONG Zhaoyun, et al. Seismic sparse inversion in mixed-domain utilizing fast matching pursuit algorithm[J]. Journal of China University of Petroleum (Edition of Natural Science), 2018, 42(1): 50-59. DOI:10.3969/j.issn.1673-5005.2018.01.006 |
| [21] |
LI Kun, YIN Xingyao, ZONG Zhaoyun. Pre-stack Bayesian cascade AVA inversion in complex-Laplace domain and its application to the broadband data acquired at East China[J]. Journal of Petroleum Science and Engineering, 2017, 158(4): 751-765. |