近年来, 非常规油气在能源领域逐渐受到人们的重视。然而, 不同于以寻找圈闭为主的常规油气勘探, 非常规油气属于“连续聚集型”油气藏, 需要找到富集且易于开发的区域, 即“甜点”区域^[1]。在“甜点”预测中, 杨氏模量和泊松比的预测具有重要意义, 低泊松比、高杨氏模量往往意味着页岩脆性好。还可以利用泊松比和杨氏模量计算得到岩石的脆性指数^[2-3], 从而评价“甜点”, 为后续的水平井部署、水力压裂提供指导依据。

叠前地震反演技术利用叠前道集的AVO信息和测井数据提供的先验信息直接提取多种弹性参数数据体(纵波阻抗、横波阻抗及密度), 进而结合弹性参数之间的转换关系, 间接获取杨氏模量和泊松比。然而由于弹性参数转换引入的累计误差, 使得常规叠前地震反演技术预测的杨氏模量和泊松比精度不能满足解释要求。为此, 宗兆云等^[4]对Aki-Richards近似公式进行了重新推导, 建立了关于角度反射系数与杨氏模量(Young's modulus)反射率、泊松比(Poisson's ratio)反射率及密度(Density)反射率的YPD近似公式, 在此基础上通过整合贝叶斯线性反演框架, 直接从叠前道集中提取杨氏模量和泊松比。

然而, Zoeppritz方程的线性近似公式的使用需要满足中小入射角、常背景纵横波速度比以及弱弹性参数变化率等诸多假设, 限制了其在不同地质条件下的应用, 为此许多学者针对精确Zoeppritz方程的叠前反演开展了研究。李录明等^[5]针对精确Zoeppritz方程开展了广义线性多波反演研究; 魏修成等^[6]将Zoeppritz方程表达为关于射线参数和弹性参数的函数, 并提出了基于射线参数道集的叠前地震反演技术; 张丰麒等^[7]和FANG等^[8]为了进一步减少反演参数, 将精确Zoeppritz方程改写为关于上、下地层的纵波速度比、横波速度比、密度比以及上层纵横波速度比的函数, 然而该方法仍需要假设常背景纵横波速度比; 张广智等^[9-10]基于精确Zoeppritz方程, 结合贝叶斯理论和马尔科夫链蒙特卡罗(MCMC)算法从后验概率密度函数中抽取反演的解, 该方法属于随机地震反演, 需要大量反复的随机模拟才能实现。ZHOU等^[11]基于各向同性假设, 将精确Zoeppritz方程转换成含有杨氏模量(Y)、泊松比(P)和密度(ρ)的形式, 并进行直接反演。宋建国等^[12]推导了比值均方根形式的杨氏模量、泊松比和密度的Zoeppritz方程, 并基于广义线性反演构建了杨氏模量和泊松比直接反演方法。

本文以贝叶斯广义线性反演为框架, 将精确Zoeppritz方程重组为关于杨氏模量、泊松比和密度自然对数的函数, 即YPD-Zoeppritz方程, 结合一阶差分矩阵和三变量柯西分布, 构建了更为合理的反射系数稀疏约束项, 在此基础上, 进一步引入低频软约束项。最后利用模型试算和实际数据试算, 对方法进行了验证。

1 YPD-Zoeppritz方程的推导

为了能够直接反演出杨氏模量、泊松比, 对Zoeppritz方程进行重新推导。考虑到弹性参数之间存在如下转换关系:

$ P = \left( {v_{\rm{P}}^2 - 2v_{\rm{S}}^2} \right)/\left( {2v_{\rm{P}}^2 - 2v_{\rm{S}}^2} \right) $

(1a)

$ Y = 2v_{\rm{S}}^2\rho \left[ {\left( {3v_{\rm{P}}^2 - 4v_{\rm{S}}^2} \right)/\left( {2v_{\rm{P}}^2 - 2v_{\rm{S}}^2} \right)} \right] $

(1b)

式中:P表示杨氏模量; Y表示泊松比; ρ表示密度。

公式(1b)可以转化为:

$ Y = 2v_{\rm{S}}^2\rho \left( {\frac{{v_{\rm{P}}^2 - 2v_{\rm{S}}^2}}{{2v_{\rm{P}}^2 - 2v_{\rm{S}}^2}} + 1} \right) = 2v_{\rm{S}}^2\rho \left( {P + 1} \right) $

(2)

由公式(2)得到横波速度v_S关于杨氏模量、泊松比和密度的表达式:

$ {v_{\rm{S}}} = \sqrt {\frac{Y}{{2\rho \left( {P + 1} \right)}}} $

(3)

将(3)式代入公式(1a), 得到纵波速度v_P关于杨氏模量、泊松比和密度的关系式:

$ {v_{\rm{P}}} = \sqrt {\frac{Y}{{\rho \left( {P + 1} \right)}}\frac{{\left( {P - 1} \right)}}{{\left( {2P - 1} \right)}}} $

(4)

为了便于后续表达, 将公式(4)和公式(3)分别简化为v_P=f(Y, P, ρ)和v_S=g(Y, P, ρ), 并将其代入精确Zoeppritz方程得到其关于杨氏模量、泊松比和密度的表达式:

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {\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{{f\left( {{Y_1},{P_1},{\rho _1}} \right)}}{{g\left( {{Y_1},{P_1},{\rho _1}} \right)}}\cos 2{\beta _1}}&{\frac{{f\left( {{Y_1},{P_1},{\rho _1}} \right)}}{{g\left( {{Y_2},{P_2},{\rho _2}} \right)}}\frac{{{g^2}\left( {{Y_2},{P_2},{\rho _2}} \right)}}{{{g^2}\left( {{Y_1},{P_1},{\rho _1}} \right)}}\frac{{{\rho _2}}}{{{\rho _1}}}\sin 2{\alpha _2}}&{\frac{{{\rho _2}}}{{{\rho _1}}}\frac{{f\left( {{Y_1},{P_1},{\rho _1}} \right)g\left( {{Y_2},{P_2},{\rho _2}} \right)}}{{{g^2}\left( {{Y_1},{P_1},{\rho _1}} \right)}}\cos 2{\beta _2}}\\ {\cos 2{\beta _1}}&{ - \frac{{f\left( {{Y_1},{P_1},{\rho _1}} \right)}}{{g\left( {{Y_1},{P_1},{\rho _1}} \right)}}\sin 2{\beta _1}}&{ - \frac{{{\rho _2}}}{{{\rho _1}}}\frac{{f\left( {{Y_2},{P_2},{\rho _2}} \right)}}{{f\left( {{Y_1},{P_1},{\rho _1}} \right)}}\cos 2{\beta _2}}&{ - \frac{{{\rho _2}}}{{{\rho _1}}}\frac{{g\left( {{Y_2},{P_2},{\rho _2}} \right)}}{{f\left( {{Y_1},{P_1},{\rho _1}} \right)}}\sin 2{\beta _2}} \end{array}} \right] \cdot \\ \left[ {\begin{array}{*{20}{c}} {{R_{{\rm{PP}}}}}\\ {{R_{{\rm{PS}}}}}\\ {{T_{{\rm{PP}}}}}\\ {{T_{{\rm{PS}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \sin {\alpha _1}}\\ {\cos {\alpha _1}}\\ {\sin 2{\alpha _1}}\\ { - \cos 2{\beta _1}} \end{array}} \right] \end{array} $

(5)

另外, 反演过程中需要引入正则化项来提高解的稳定性。通常的做法是结合贝叶斯原理假设反演参数服从某一先验分布。先验分布分为高斯分布和长尾巴分布:高斯分布对应一致性加权矩阵, 其本质是阻尼最小二乘解对应的阻尼系数项, 在参数提取过程中会产生一个相对平滑的结果, 而高分辨率反演需要假设弹性参数反射率服从长尾巴分布, 从而构建稀疏约束项。

然而由公式(5)可以看出, 反演参数为弹性参数而非弹性参数反射率, 因此直接假设其服从长尾巴分布无法达到提高分辨率的目的, 这是因为非一致性加权矩阵会对偏离背景值较大的弹性参数产生较小的加权系数, 对偏离背景值较小的弹性参数产生较大的加权系数, 这样就会凸显偏离背景值较大的弹性参数, 压制偏离背景值较小的弹性参数, 而并不是产生“块化”的反演效果, 因此在理论上也无法达到提高分辨率的目的。

弹性参数(E)与其反射率(ΔE/E)的自然对数存在如下关系:

$ \frac{{\Delta E}}{E} \approx \Delta \ln E = \ln {E_{i + 1}} - \ln {E_i} $

(6)

利用公式(6)建立弹性参数反射率与弹性参数自然对数的关系, 再借助一阶差分矩阵, 间接使得弹性参数反射率服从长尾巴分布, 进而构建更为合理的约束项。

基于以上分析可知, 需要建立以杨氏模量自然对数Y、泊松比自然对数P和密度自然对数ρ为未知参数的Zoeppritz方程。对公式(3)和公式(4)中的杨氏模量、泊松比和密度取对数后再指数化, 可以得到:

$ {v_{\rm{S}}} = \sqrt {\frac{{\exp \bar Y}}{{2\exp \left( {\bar \rho + \bar P} \right) + 2\exp \left( {\bar \rho } \right)}}} $

(7)

$ {v_{\rm{P}}} = \sqrt {\frac{{\exp \left( {\bar Y - \bar \rho } \right)}}{{\left[ {\exp \left( {\bar P} \right) + 1} \right]}}\frac{{\left[ {\exp \left( {\bar P} \right) - 1} \right]}}{{\left[ {2\exp \left( {\bar P} \right) - 1} \right]}}} $

(8)

将公式(8)和公式(7)分别简化为v_P=γ(Y, P, ρ)和v_S=ζ(Y, P, ρ), 并将其代入精确Zoeppritz方程得到:

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {\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{{\gamma \left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}{{\zeta \left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}\cos 2{\beta _1}}&{\frac{{\gamma \left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}{{\gamma \left( {{{\bar Y}_2},{{\bar P}_2},{{\bar \rho }_2}} \right)}}\frac{{{\zeta ^2}\left( {{{\bar Y}_2},{{\bar P}_2},{{\bar \rho }_2}} \right)}}{{{\zeta ^2}\left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}\frac{{{\rho _2}}}{{{\rho _1}}}\sin 2{\alpha _2}}&{\frac{{{\rho _2}}}{{{\rho _1}}}\frac{{\gamma \left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)\zeta \left( {{{\bar Y}_2},{{\bar P}_2},{{\bar \rho }_2}} \right)}}{{{\zeta ^2}\left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}\cos 2{\beta _2}}\\ {\cos 2{\beta _1}}&{ - \frac{{\gamma \left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}{{\zeta \left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}\sin 2{\beta _1}}&{ - \frac{{{\rho _2}}}{{{\rho _1}}}\frac{{\gamma \left( {{{\bar Y}_2},{{\bar P}_2},{{\bar \rho }_2}} \right)}}{{\gamma \left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}\cos 2{\beta _2}}&{ - \frac{{{\rho _2}}}{{{\rho _1}}}\frac{{\zeta \left( {{{\bar Y}_2},{{\bar P}_2},{{\bar \rho }_2}} \right)}}{{\gamma \left( {{{\bar Y}_1},{{\bar P}_1},{{\bar \rho }_1}} \right)}}\sin 2{\beta _2}} \end{array}} \right] \cdot \\ \left[ {\begin{array}{*{20}{c}} {{R_{{\rm{PP}}}}}\\ {{R_{{\rm{PS}}}}}\\ {{T_{{\rm{PP}}}}}\\ {{T_{{\rm{PS}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \sin {\alpha _1}}\\ {\cos {\alpha _1}}\\ {\sin 2{\alpha _1}}\\ { - \cos 2{\beta _1}} \end{array}} \right] \end{array} $

(9)

公式(9)建立了角度反射系数与杨氏模量自然对数、泊松比自然对数及密度自然对数之间的关系, 称之为YPD-Zoeppritz方程。该公式的推导并未引入任何假设和近似, 因此其精度没有任何损失。为了更为直观地展现该公式的精确性, 我们利用一个双层介质模型(表 1)进行了试算。

首先利用精确Zoeppritz方程计算0~90°入射角对应的PP波反射系数与PS波反射系数, 如图 1和图 2中红色实线所示, 然后将表 1中的弹性参数转换为杨氏模量自然对数、泊松比自然对数以及密度自然对数, 再利用YPD-Zoeppritz方程计算0~90°入射角对应的PP波反射系数与PS波反射系数, 如图 1和图 2中蓝圈所示。可以看出, 在0~90°范围内, YPD-Zoeppritz方程与精确Zoeppritz方程所计算的反射系数完全重合。

2 贝叶斯广义线性反演

公式(9)建立了从杨氏模量、泊松比和密度到角度反射系数的正演过程, 由于该方程是关于这3个弹性参数的非线性方程, 通常利用模拟退火、MCMC等全局寻优算法进行反演, 但是由于该类寻优算法需要反复模拟, 计算量较大, 并且反演结果的高频成分具有一定的随机性, 造成单一逐道反演模式下其横向连续性较差(随机反演的横向连续性需要进一步结合空间约束或水平变差等地质统计学加以改善, 常规的逐道反演并不能保证其横向连续性)等原因, 并未得到大规模应用。针对非线性反演问题, 还可以借助于广义线性反演的思想, 利用迭代的思路逐步获得合理的反演结果。

将某一入射角θ的角度部分叠加数据d(θ)在初始解m₀=[Y₀ P₀ ρ₀]^T附近进行泰勒展开, 并忽略二阶以上的高阶项, 得:

$ \begin{array}{l} \mathit{\boldsymbol{d}}\left( \theta \right) \approx \mathit{\boldsymbol{d}}\left( \theta \right)\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right. + \frac{{\partial \mathit{\boldsymbol{d}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar Y}} + \frac{{\partial \mathit{\boldsymbol{d}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar P}} + \\ \;\;\;\;\;\;\;\;\;\;\frac{{\partial \mathit{\boldsymbol{d}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar \rho }} \end{array} $

(10)

式中:d(θ)|_m0表示由初始解m₀得到的合成地震数据; ΔY, ΔP和Δρ分别表示杨氏模量自然对数、泊松比自然对数以及密度自然对数的更新量。利用褶积模型, 可将d(θ)表示为:

$ \mathit{\boldsymbol{d}}\left( \theta \right) = \mathit{\boldsymbol{W}}\left( \theta \right){\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( \theta \right) $

(11)

式中:W(θ)表示角度子波褶积矩阵; R_PP(θ)表示PP波角度反射系数, 可由YPD-Zoeppritz方程正演得到。将公式(11)代入公式(10)中, 得:

$ \begin{array}{l} \mathit{\boldsymbol{d}}\left( \theta \right) \approx \mathit{\boldsymbol{d}}\left( \theta \right)\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right. + \frac{{\partial \mathit{\boldsymbol{W}}\left( \theta \right){\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar Y}} + \\ \frac{{\partial \mathit{\boldsymbol{W}}\left( \theta \right){\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar P}} + \frac{{\partial \mathit{\boldsymbol{W}}\left( \theta \right){\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar \rho }} \end{array} $

(12)

由于角度子波为常向量, 并不随杨氏模量、泊松比和密度的变化而变化, 因此可以将子波褶积矩阵移到偏导项外面:

$ \begin{array}{l} \Delta \mathit{\boldsymbol{d}}\left( \theta \right) = \mathit{\boldsymbol{W}}\left( \theta \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar Y}} + \mathit{\boldsymbol{W}}\left( \theta \right) \cdot \\ \;\;\;\;\;\;\;\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar P}} + \mathit{\boldsymbol{W}}\left( \theta \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( \theta \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar \rho }} \end{array} $

(13)

式中:Δd(θ)表示地震数据的残差向量。如要得到反演参数m的更新量, 至少需要联立3个角度进行求解:

$ \left\{ \begin{array}{l} \Delta \mathit{\boldsymbol{d}}\left( {{\theta _1}} \right) = \mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _1}} \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar Y}} + \mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _1}} \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar P}} + \mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _1}} \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar \rho }}\\ \Delta \mathit{\boldsymbol{d}}\left( {{\theta _2}} \right) = W\left( {{\theta _2}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _2}} \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar Y}} + \mathit{\boldsymbol{W}}\left( {{\theta _2}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _2}} \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar P}} + \mathit{\boldsymbol{W}}\left( {{\theta _2}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _2}} \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar \rho }}\\ \Delta \mathit{\boldsymbol{d}}\left( {{\theta _3}} \right) = \mathit{\boldsymbol{W}}\left( {{\theta _3}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _3}} \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar Y}} + \mathit{\boldsymbol{W}}\left( {{\theta _3}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _3}} \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar P}} + \mathit{\boldsymbol{W}}\left( {{\theta _3}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _3}} \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{\bar \rho }} \end{array} \right. $

(14)

将其表示为矩阵形式:

$ \left[ {\begin{array}{*{20}{c}} {\Delta \mathit{\boldsymbol{d}}\left( {{\theta _1}} \right)}\\ {\Delta \mathit{\boldsymbol{d}}\left( {{\theta _2}} \right)}\\ {\Delta \mathit{\boldsymbol{d}}\left( {{\theta _3}} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _1}} \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.}&{\mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _1}} \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.}&{\mathit{\boldsymbol{W}}\left( {{\theta _1}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _1}} \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.}\\ {\mathit{\boldsymbol{W}}\left( {{\theta _2}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _2}} \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.}&{\mathit{\boldsymbol{W}}\left( {{\theta _2}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _2}} \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.}&{\mathit{\boldsymbol{W}}\left( {{\theta _2}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _2}} \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.}\\ {\mathit{\boldsymbol{W}}\left( {{\theta _3}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _3}} \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.}&{\mathit{\boldsymbol{W}}\left( {{\theta _3}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _3}} \right)}}{{\partial \mathit{\boldsymbol{\bar P}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.}&{\mathit{\boldsymbol{W}}\left( {{\theta _3}} \right)\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _3}} \right)}}{{\partial \mathit{\boldsymbol{\bar \rho }}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\Delta \mathit{\boldsymbol{\bar Y}}}\\ {\Delta \mathit{\boldsymbol{\bar P}}}\\ {\Delta \mathit{\boldsymbol{\bar \rho }}} \end{array}} \right] $

(15)

考虑到第i层和第i+1层分界面对应的角度反射系数是第i层和第i+1层弹性参数的函数, 因此对于N层介质来说, 角度反射系数向量关于某一弹性参数向量的一阶偏导数是一个N-1行、N列的双对角矩阵, 如$\frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _1}} \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.$为:

$ \frac{{\partial {\mathit{\boldsymbol{R}}_{{\rm{PP}}}}\left( {{\theta _1}} \right)}}{{\partial \mathit{\boldsymbol{\bar Y}}}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right. = {\left[ {\begin{array}{*{20}{c}} {\frac{{\partial R_{{\rm{PP}}}^1\left( {{\theta _1}} \right)}}{{\partial {{\bar Y}_1}}}}&{\frac{{\partial R_{{\rm{PP}}}^1\left( {{\theta _1}} \right)}}{{\partial {{\bar Y}_2}}}}&{}&{}&{}\\ {}&{\frac{{\partial R_{{\rm{PP}}}^2\left( {{\theta _1}} \right)}}{{\partial {{\bar Y}_2}}}}&{\frac{{\partial R_{{\rm{PP}}}^2\left( {{\theta _1}} \right)}}{{\partial {{\bar Y}_3}}}}&{}&{}\\ {}&{}& \cdots &{}&{}\\ {}&{}&{}&{\frac{{\partial R_{{\rm{PP}}}^{N - 1}\left( {{\theta _1}} \right)}}{{\partial {{\bar Y}_{N - 1}}}}}&{\frac{{\partial R_{{\rm{PP}}}^{N - 1}\left( {{\theta _1}} \right)}}{{\partial {{\bar Y}_N}}}} \end{array}} \right]_{\left( {N - 1} \right) \times N}} $

(16)

式中:$\left[ {\partial R_{{\rm{PP}}}^i\left( {{\theta _1}} \right)} \right]/\left( {\partial {{\bar Y}_i}} \right)$表示第i层介质和第i+1层介质分界面对应的角度反射系数关于第i层介质的杨氏模量自然对数的一阶偏导数, 一阶偏导数可以用解析方式求取^[4-5], 也可以利用数值方式求取^[2]。为了后续表达方便, 公式(15)可以简化为:

$ \Delta \mathit{\boldsymbol{d}} = \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\Delta \mathit{\boldsymbol{m}} $

(17)

通过给定初始解m₀, 利用(17)式反复修正反演参数, 直至获取合理准确的解。然而由于子波带限、角度反射系数向弹性参数分解过程的病态性等原因, 造成Jacobian矩阵G|_m0的奇异性, 进而导致反演解的不稳定和非唯一。这一问题的解决, 可以借助贝叶斯原理, 通过假设反演参数服从某一先验分布, 进而引入反演的正则化项, 有效降低广义线性反演的不适定性。通常反演参数m的先验分布利用测井数据获取, 但是反演参数的更新量Δm的先验分布难以获取, 因此需要将公式(17)转化为以m为待求参数的形式, 并通过在公式(17)左右两边加上G|_m0m₀的方式加以解决:

$ \mathit{\boldsymbol{H}} = \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}} $

(18)

其中, H=Δd+G|_m0m₀。

贝叶斯反演通过引入反演参数的先验约束, 可以有效减少解空间的搜索范围, 提高解的稳定性和合理性。利用极大化反演参数后验概率密度函数得到相应的反演目标函数, 反演参数的先验分布等价反演目标函数的正则化项, 因此能够降低反演问题的不适定性。将广义线性反演的迭代思想与贝叶斯反演相结合, 引入先验分布, 可以有效提高广义线性反演在迭代过程中的稳定性。

通常假设弹性参数反射率服从长尾巴分布, 可以产生稀疏脉冲反演的效果。在有地质背景的情况下, 同一时间采样点的杨氏模量反射率、泊松比反射率以及密度反射率之间存在一定的相关性, 因此三变量柯西分布更适合作为其先验分布。N层介质对应N-1个分界面, 其弹性参数反射率的联合概率密度函数为:

$ p\left( \mathit{\boldsymbol{r}} \right) = \prod\limits_{i = 1}^{N - 1} {\frac{{2{{\left| \mathit{\boldsymbol{\psi }} \right|}^{ - \frac{1}{2}}}}}{{{{\rm{ \mathsf{ π} }}^2}{{\left( {1 + {\mathit{\boldsymbol{r}}^{\rm{T}}}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{r}}} \right)}^2}}}} $

(19)

式中:r表示弹性参数反射率, $\mathit{\boldsymbol{r}} = {\left[ {\begin{array}{*{20}{c}} {\frac{{\Delta \mathit{\boldsymbol{Y}}}}{\mathit{\boldsymbol{Y}}}}&{\frac{{\Delta \mathit{\boldsymbol{P}}}}{\mathit{\boldsymbol{P}}}}&{\frac{{\Delta \mathit{\boldsymbol{\rho }}}}{\mathit{\boldsymbol{\rho }}}} \end{array}} \right]^{\rm{T}}}$。根据公式(6), 弹性参数反射率向量可以表示为一阶差分矩阵和弹性参数自然对数向量相乘:

$ \mathit{\boldsymbol{r}} = \mathit{\boldsymbol{Dm}} $

(20)

其中, D为1阶差分矩阵。

$ \mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} { - 1}&1&{}&{}&{}\\ {}&{ - 1}&1&{}&{}\\ {}&{}& \cdots &{}&{}\\ {}&{}&{}&{ - 1}&1 \end{array}} \right] $

(21)

因此(19)式可以变换为以弹性参数自然对数向量m为自变量的函数:

$ p\left( \mathit{\boldsymbol{m}} \right) = \prod\limits_{i = 1}^N {\frac{{2{{\left| \mathit{\boldsymbol{\psi }} \right|}^{ - \frac{1}{2}}}}}{{{{\rm{ \mathsf{ π} }}^2}{{\left( {1 + {\mathit{\boldsymbol{m}}^{\rm{T}}}{\mathit{\boldsymbol{D}}^{\rm{T}}}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{Dm}}} \right)}^2}}}} $

(22)

式中:ψ表示描述弹性参数反射率之间相关性的尺度矩阵, 可以利用最大似然估计从测井曲线中获取; ALEMIE等^[13]给出了φ_i的具体表达式。

通过假设噪声服从多变量高斯分布来表征和描述合成地震数据与实测地震数据之间的相似性, 称之为似然函数:

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{H}}\left| \mathit{\boldsymbol{m}} \right.} \right) = \frac{1}{{{{\left( {2{\rm{ \mathsf{ π} }}} \right)}^{\left( {3 \times N - 3} \right)/2}}{{\left| {{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{d}}}} \right|}^{1/2}}}} \cdot \\ \;\;\;\;\;\exp \left[ { - \frac{{{{\left( {\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}}} \right)}}{{2{\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{d}}}}}} \right] \end{array} $

(23)

式中:C_d表示噪声的协方差矩阵。由于不同角度、不同时间采样点的噪声独立不相关, 因此C_d=σ_d²I, I表示单位矩阵, σ_d²表示噪声的方差。

根据贝叶斯公式, 反演参数的后验概率可以表示为:

$ p\left( {\mathit{\boldsymbol{m}}\left| \mathit{\boldsymbol{H}} \right.} \right) = \frac{{p\left( \mathit{\boldsymbol{m}} \right)p\left( {\mathit{\boldsymbol{H}}\left| \mathit{\boldsymbol{m}} \right.} \right)}}{{p\left( \mathit{\boldsymbol{H}} \right)}} $

(24)

由于p(H)只与实测地震数据d和初始模型m₀有关, 并且这2种数据均为已知, 因此p(H)是一常数, 不影响后验概率密度函数的形状。因此:

$ p\left( {\mathit{\boldsymbol{m}}\left| \mathit{\boldsymbol{H}} \right.} \right) \propto p\left( \mathit{\boldsymbol{m}} \right)p\left( {\mathit{\boldsymbol{H}}\left| \mathit{\boldsymbol{m}} \right.} \right) $

(25)

将公式(22)和公式(23)代入公式(25), 并忽略其中的常系数, 可得:

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{m}}\left| \mathit{\boldsymbol{H}} \right.} \right) \propto \exp \left[ { - \frac{{{{\left( {\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}}} \right)}}{{2\sigma _d^2\mathit{\boldsymbol{I}}}}} \right] \cdot \\ \;\;\;\;\;\;\prod\limits_{i = 1}^N {\frac{1}{{{{\left( {1 + {\mathit{\boldsymbol{m}}^{\rm{T}}}{\mathit{\boldsymbol{D}}^{\rm{T}}}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{Dm}}} \right)}^2}}}} \end{array} $

(26)

贝叶斯反演通过求取后验概率密度的极大值来获取反演的解, 这一过程可以通过对公式(26)左右两边取对数后再取负号, 转为极小化反演目标函数, 即$\max \left[ {\mathit{\boldsymbol{p}}\left( {\mathit{\boldsymbol{m}}|\mathit{\boldsymbol{H}}} \right)} \right] \Leftrightarrow \min \left[ {{\rm{obj}}\left( \mathit{\boldsymbol{m}} \right)} \right]$, 亦即:

$ \begin{array}{l} {\rm{obj}}\left( \mathit{\boldsymbol{m}} \right) = \frac{1}{2}{\left( {\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}}} \right) + \\ \sigma _d^2\sum\limits_{i = 1}^N {\ln \left( {1 + {\mathit{\boldsymbol{m}}^{\rm{T}}}{\mathit{\boldsymbol{D}}^{\rm{T}}}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{Dm}}} \right)} \end{array} $

(27)

公式(27)右边第一项表示地震数据的拟合差, 第二项表示反射系数的稀疏约束项, 由于地震数据中缺失0~10Hz的低频分量, 直接利用公式(27)无法得到全频带的反演结果, 因此需要进一步引入低频软约束项, 使反演结果的低频成分在L2范数准则下逼近测井数据的低频趋势, 从而使反演结果拥有正确的低频信息。引入低频软约束项后的反演目标函数为:

$ \begin{array}{l} {\rm{obj}}\left( \mathit{\boldsymbol{m}} \right) = \frac{1}{2}{\left( {\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right.\mathit{\boldsymbol{m}}} \right) + \\ \sigma _d^2\sum\limits_{i = 1}^N {\ln \left( {1 + {\mathit{\boldsymbol{m}}^{\rm{T}}}{\mathit{\boldsymbol{D}}^{\rm{T}}}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{Dm}}} \right)} + {C_{{\rm{LFC}}}}\left( \mathit{\boldsymbol{m}} \right) \end{array} $

(28)

其中, C_LFC(m)表示低频软约束项, 表达式为:

$ \begin{array}{l} {C_{{\rm{LFC}}}}\left( \mathit{\boldsymbol{m}} \right) = \frac{{{\lambda _Y}}}{2}{\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{Y}}}\mathit{\boldsymbol{m}} - {{\mathit{\boldsymbol{\bar Y}}}_{{\rm{Trend}}}}} \right)^{\rm{T}}}\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{Y}}}\mathit{\boldsymbol{m}} - {{\mathit{\boldsymbol{\bar Y}}}_{{\rm{Trend}}}}} \right) + \\ \frac{{{\lambda _P}}}{2}{\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{P}}}\mathit{\boldsymbol{m}} - {{\mathit{\boldsymbol{\bar P}}}_{{\rm{Trend}}}}} \right)^{\rm{T}}}\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{P}}}\mathit{\boldsymbol{m}} - {{\mathit{\boldsymbol{\bar P}}}_{{\rm{Trend}}}}} \right) + \frac{{{\lambda _\rho }}}{2}\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{\rho }}}\mathit{\boldsymbol{m}} - } \right.\\ {\left. {{{\mathit{\boldsymbol{\bar \rho }}}_{{\rm{Trend}}}}} \right)^{\rm{T}}}\left( {{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{\rho }}}\mathit{\boldsymbol{m}} - {{\mathit{\boldsymbol{\bar \rho }}}_{{\rm{Trend}}}}} \right) \end{array} $

(29)

式中:λ_Y, λ_P, λ_ρ分别为杨氏模量、泊松比和密度的低频软约束权重; Y_Trend, P_Trend, ρ_Trend分别表示杨氏模量自然对数、泊松比自然对数以及密度自然对数的低频趋势, 可以通过在层位和沉积模式双重控制下的测井曲线横向插值, 并结合低通滤波获取; L表示低通滤波算子。低通滤波算子L由张丰麒等^[14]提出, 具体表达式见附录A。

为了获取目标函数极小值, 需要对公式(28)左右两边求偏导并令其为0:

$ \begin{array}{l} \left( {\mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}^T} \right.\mathit{\boldsymbol{G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}} \right. + \sigma _d^2\mathit{\boldsymbol{Q}} + {\lambda _Y}\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{Y}}^{\rm{T}}{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{Y}}} + {\lambda _P}\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{P}}^{\rm{T}}{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{P}}} + {\lambda _\rho }\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{\rho }}^{\rm{T}}{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{\rho }}}} \right)\mathit{\boldsymbol{m}}\\ \mathit{\boldsymbol{ = G}}\left| {_{{\mathit{\boldsymbol{m}}_0}}^T} \right.\mathit{\boldsymbol{d}} + {\lambda _Y}\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{Y}}^{\rm{T}}{{\mathit{\boldsymbol{\bar Y}}}_{{\rm{Trend}}}} + {\lambda _P}\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{P}}^{\rm{T}}{{\mathit{\boldsymbol{\bar P}}}_{{\rm{Trend}}}} + {\lambda _\rho }\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{\rho }}^{\rm{T}}{{\mathit{\boldsymbol{\bar \rho }}}_{{\rm{Trend}}}} \end{array} $

(30)

其中, Q表示非一致性加权矩阵, 其表达式为:

$ {\left[ \mathit{\boldsymbol{Q}} \right]_{xy}} = \sum\limits_{i = 1}^N {\frac{{2{{\left[ {{\mathit{\boldsymbol{\varphi }}_i}} \right]}_{xy}}}}{{1 + {\mathit{\boldsymbol{m}}^{\rm{T}}}{\mathit{\boldsymbol{D}}^{\rm{T}}}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{Dm}}}}} $

(31)

式中:[*]_xy表示矩阵*第x行、第y列元素。

针对(30)式的求解, 可以利用迭代法进行:

1) 在层位和沉积模式的双重控制下对测井曲线进行横向插值, 并结合低通滤波, 获取低频趋势Y_Trend, P_Trend, ρ_Trend;

2) 根据附录A, 计算低通滤波算子L_Y, L_P, L_ρ;

3) 根据初始模型(首次迭代可以利用低频趋势作为初始模型, m₀=[Y_Trend P_Trend ρ_Trend]^T)计算Jacobian矩阵G|_m0、H矩阵及非线性一致性加权矩阵Q;

4) 根据公式(30), 利用Cholesky分解和三角矩阵分解计算反演参数m;

5) 令下一次迭代的初始模型为更新后的反演参数m→m₀;

6) 重复步骤3)~步骤5), 直到预定的收敛阈值或者迭代次数;

7) 对最终的反演结果取指数, 获取弹性参数。

3 模型试算

为了检验本方法的有效性和精度, 进行如下模型试算。基于图 3给出的弹性参数模型, 不考虑多次波、透射波以及绕射波的存在, 只考虑一次反射的情况下, 利用精确Zoeppritz方程分别计算入射角分别为10°、20°以及30°时的PP波反射系数, 然后分别与主频为30Hz的雷克子波进行褶积, 得到图 4a所示的合成角道集。图 4b为含随机噪声合成角道集(信噪比为4)。

首先利用本文方法进行反演。将弹性参数转化为杨氏模量和泊松比, 再进行10Hz的低通滤波获取初始模型。并令其反射系数稀疏约束权重为0.05, 杨氏模量、泊松比以及密度的低频软约束项权重均为1, 迭代次数为10次, 其反演结果如图 5所示。可以看出, 反演结果(红线)与实测曲线(黑线)吻合度较好, 特别是杨氏模量和泊松比两项, 密度反演结果精度稍低, 但仍优于常规AVO反演精度。在三变量柯西分布约束下, 反演结果呈现明显的“块化”效果。图 6给出了对应的弹性参数反射率, 可以看出反演的反射率具有明显的尖脉冲特性, 子波的旁瓣效应得到了明显压制, 因此有利于降低后续解释的不确定性。为了进一步检验本文方法的精度, 给出了基于常规AVO三参数反演间接获取的杨氏模量和泊松比反演结果(图 7), 由图 7可以看出, 由于弹性参数转换带来的累计误差以及近似公式引入的正演误差等, 造成反演结果与模型数据的吻合度较差。

为了检验本文方法的抗噪性, 对图 4b所示合成角道集进行反演计算, 在相同的初始模型、子波、稀疏约束权重、低频软约束权重以及迭代次数下, 反演结果如图 8所示。可以看出, 由于随机噪声的影响, 弱反射区域的反演结果与模型数据吻合度较差, 但是整体上吻合度较好, 因此可以证明本文方法具有较强的抗噪性。

4 实际数据应用

实际数据来自国内某非常规油气勘探工区, 工区内岩性主要为致密砂岩, 为了研究本工区脆性指数的空间展布, 需要利用叠前弹性参数反演获取杨氏模量和泊松比数据体。经过前期的线性噪声压制、多道统计反褶积、叠加速度拾取与多次波压制交替进行、真振幅叠前时间偏移、偏移距域向角度域道集的转换以及角度部分叠加等处理流程, 分别形成了0~12°, 10°~22°以及20°~32°叠加的3个角度部分叠加数据体, 图 9给出了对应的全角度叠加剖面。在对工区内测井曲线进行去除野值、环境校正以及一致性校正, 再结合井震标定确定深时关系后, 分别提取了对应的角度子波。在层位和沉积模式的双重控制下, 对时间域测井曲线横向插值并进行低通滤波, 作为反演的初始模型, 滤波参数在频谱分析后确定, 本次滤波参数定为10Hz。设置反演的稀疏约束权重为0.5、低频软约束权重为1、迭代次数为5次, 图 10、图 11以及图 12分别给出了对应的杨氏模量、泊松比、密度反演剖面以及相应的未做低频软约束的反演结果, 可以明显看出低频软约束项提高了剖面的连续性。整体看, 反演剖面具有明显的层次感, 在低频软约束项的约束下, 稳定的低频成分使得反演剖面具有较好的横向连续性, 反演剖面细节丰富, 特别是泊松比反演剖面, 可以展现储层微弱的横向变化。B井作为验证井, 未参与到反演当中, 将其杨氏模量曲线、泊松比曲线以及密度曲线分别进行250Hz Backus滤波后, 再进行深时转换, 分别叠置于对应的反演剖面上(如图 10、图 11、图 12所示), 可以看出反演结果与实测井曲线有较好的吻合度, 验证了本文方法在实际数据应用中的有效性和精度。

5 结论

本文从精确Zoeppritz方程出发, 直接建立角度反射系数与杨氏模量、泊松比和密度的关系, 在推导过程中, 并未引入任何假设, 因此YPD-Zoeppritz方程是全精度的, 适合于大角度下的叠前地震反演。

广义线性反演利用迭代算法可以获取高精度的反演结果, 但通常Jacobian条件数巨大, 造成广义线性反演的不适定性比常规线性反演高, 通过结合贝叶斯理论引入模型参数的先验分布, 可以有效缩小反演寻优的搜索范围, 正则化项的引入可以有效提高广义线性反演的适定性。

将一阶差分矩阵与三变量柯西分布相结合, 在以弹性参数的自然对数为反演参数的情况下, 可以完美构建反射系数稀疏约束项, 进而在迭代过程中, 逐步压缩子波旁瓣, 使得反演结果呈现明显的“块化”效果, 降低子波旁瓣对反演结果的影响。

模型试算和实际数据测试结果表明本文方法具有较高的精度和稳定性, 具有一定的推广应用价值。

6 附录A

低频软约束中的低通滤波矩阵L的具体表达式为:

$ \mathit{\boldsymbol{L}} = {\mathit{\boldsymbol{E}}^{\rm{T}}}{\mathit{\boldsymbol{F}}^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varLambda} FE}} $

(A1)

(A1) 式中E矩阵起到将模型参数进行3倍延长的作用。考虑到反演目的层段时窗过短会引起频域的分辨率不够, 或者由于数据两端的截断效应会引起低通滤波结果不准确, 通过该矩阵可以使得在数据两端各以原数据的倒序进行延长, E矩阵的数学表达式为:

$ \mathit{\boldsymbol{E}} = {\left[ {\begin{array}{*{20}{c}} 0&{}&1\\ {}& {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} &{}\\ 1&{}&0\\ {}& \ddots &{}\\ 0&{}&1\\ {}& {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} &{}\\ 1&{}&0 \end{array}} \right]_{\left( {N - 2} \right) \times N}} $

(A2)

F和F^-1分别为离散傅里叶(DFT)正、反变换矩阵, 起到将数据进行傅里叶正、反变换的作用。

$ \mathit{\boldsymbol{F}} = {\left[ {\begin{array}{*{20}{c}} {{W^0}}&{{W^0}}&{{W^0}}& \cdots &{{W^0}}\\ {{W^0}}&{{W^1}}&{{W^2}}& \cdots &{{W^{N - 1}}}\\ {{W^0}}&{{W^2}}&{{W^4}}& \cdots &{{W^{N - 2}}}\\ \vdots&\vdots&\vdots &{}& \vdots \\ {{W^0}}&{{W^{N - 1}}}&{{W^{2N - 2}}}& \cdots &{{W^{\left( {N - 1} \right)2}}} \end{array}} \right]_{N \times N}} $

(A3)

$ {\mathit{\boldsymbol{F}}^{ - 1}} = \frac{1}{N} \cdot {\left[ {\begin{array}{*{20}{c}} {{W^0}}&{{W^0}}&{{W^0}}& \cdots &{{W^0}}\\ {{W^0}}&{{W^{ - 1}}}&{{W^{ - 2}}}& \cdots &{{W^{1 - N}}}\\ {{W^0}}&{{W^{ - 2}}}&{{W^{ - 4}}}& \cdots &{{W^{2 - 2N}}}\\ \vdots&\vdots&\vdots &{}& \vdots \\ {{W^0}}&{{W^{1 - N}}}&{{W^{2 - 2N}}}& \cdots &{{W^{ - \left( {N - 1} \right)2}}} \end{array}} \right]_{N \times N}} $

(A4)

(A3) 式和(A4)式中的W=cos[－(2π/N)]+jsin[－(2π/N)]。而(A1)式中的Λ矩阵为由汉宁窗函数组成的对角矩阵。