弹性阻抗反演作为叠前弹性参数反演的一个重要分支, 具有计算效率高、实现方便并且抗噪性较好的特点, 由反演获取的弹性阻抗体提取的弹性参数, 可以定量表征储层岩性或孔隙内流体的变化。

自CONNOLLY^[1]提出弹性阻抗的概念以来, 人们对弹性阻抗反演进行了广泛和深入的研究。VERWEST^[2]对弹性阻抗进行标准化处理消除了不同角度弹性阻抗之间的量纲差异; 马劲风^[3]对Zoeppritz方程进行简化, 推导出广义弹性阻抗的表达式; 王保丽等^[4]借助Gray近似式推导出基于拉梅参数的弹性阻抗表达式, 并直接提取了拉梅系数, 提高了拉梅系数的反演精度; 印兴耀等^[5]提出了流体弹性阻抗, 通过该弹性阻抗表达式可以直接提取对流体变化更为敏感的Gassmann流体项; 宗兆云等^[6]重新推导Aki-Richards近似式, 给出了基于纵横波模量的AVO近似式, 并推导了基于纵横波模量的弹性阻抗表达式; 刘晓晶等^[7]通过基追踪反演获取弹性阻抗, 随后利用f-μ弹性阻抗表达式提取Gassmann流体项和剪切模量。

前人的研究大都集中在重新推导弹性阻抗表达式, 然而弹性阻抗分解也同样至关重要。本文针对基于最小二乘原理的传统弹性阻抗分解方法在地震数据含有噪声时无法提取稳定的弹性参数这一问题, 借助贝叶斯理论引入模型参数的先验分布构建正则化项, 提高弹性阻抗分解的稳定性。通过联立三变量柯西分布和一阶差分矩阵来描述模型参数的先验分布, 构建更为合理的稀疏约束项, 并引入低频软约束项来进一步稳定弹性参数的低频成分, 消除反演剖面上的“门帘”效应。模型试算和实际数据测试结果表明, 基于双项约束的弹性阻抗分解方法提取的弹性参数在垂向上具有较高的分辨率, 横向上具有较强的连续性。

以CONNOLLY^[1]提出的标准弹性阻抗表达式为基础, 考虑到不同入射角的弹性阻抗量纲不一致, 因此标准化后的弹性阻抗表达式为:

(1)

其中,

式中:mean(*)表示对“*”求平均; θ表示入射角。

利用对数将(1) 式由乘积形式转换为加法形式:

(2)

考虑到N个时间采样点, 需要M≥3个角度的标准化弹性阻抗才能采用如下方程组提取弹性参数v_P, v_S, ρ:

(3)

其中,

为了后续表达方便, (3) 式可以简化为:

(4)

其中,

在不考虑随机噪声的情况下, 直接利用最小二乘算法便可以提取弹性参数:

(5)

然而将(5) 式直接应用于含噪声模型或实际数据时, 提取的弹性参数不稳定。主要原因是G^TG的特征值矩阵中包含非常小的特征值元素, 随机噪声的存在会扩大弹性参数提取的误差。

针对上述问题, 可以利用阻尼最小二乘算法, 即在G^TG矩阵的对角线加上一个阻尼项, 可以提高弹性阻抗分解的稳定性; 另外还可以借助贝叶斯方法, 通过描述模型参数的先验分布, 引入更为合理的先验约束项, 进而将提取弹性参数的过程转化为求取弹性参数的最大后验概率解。先验约束的引入, 可以有效提高弹性参数提取的稳定性和精度。

模型参数的先验分布分为高斯分布和长尾巴分布两种。高斯分布对应一致性加权矩阵, 其本质是阻尼最小二乘解对应的阻尼系数项在参数提取过程中会产生一个相对平滑的结果。这样就会造成弹性参数的垂向分辨率略低于弹性阻抗的垂向分辨率。

考虑到弹性参数提取过程的模型参数是三参数自然对数而非三参数反射率, 因此直接假设其服从长尾巴分布也不能达到提高分辨率的目的, 这是因为非一致性加权矩阵会对偏离背景值较大的弹性参数产生较小的加权系数, 对偏离背景值较小的弹性参数产生较大的加权系数, 这样就会凸显偏离背景值较大的弹性参数, 并压制偏离背景值较小的弹性参数, 而并不是产生“块化”效果, 因此理论上也无法达到提高分辨率的目的。

但是考虑到模型参数和三参数反射率存在如下近似:

(6)

式中:i表示第i个采样点。

同理可推导出v_S, ρ。(6) 式说明可以利用模型参数的一阶差分近似表征三参数反射率。在此基础上, 假设其服从长尾巴分布, 便可以达到提高分辨率的目的。同时在考虑地质背景的情况下, 同一时间采样点的三参数反射率并非独立, 因此选择三变量柯西分布更为合适。

考虑N个时间采样点的情况下, 三参数反射率的联合先验分布为:

(7)

式中:ψ表示尺度矩阵, 可以利用最大期望(EM)算法从测井曲线中获取; φ_i的具体表达式见文献[8]; , 为三参数反射率。利用一阶差分矩阵D建立模型参数m和三参数反射率r的关系:

(8)

其中,

(9)

因此(7) 式可以改写为:

(10)

考虑到地震数据中的随机噪声会传递到弹性阻抗反演结果中, 而且从弹性参数到弹性阻抗的正演过程不精确, 其似然函数可以用多变量高斯分布表示:

(11)

式中:C_d表示噪声的协方差矩阵。

结合贝叶斯定理, 在考虑弹性阻抗的先验分布P(d)为常数的情况下:

(12)

因此:

(13)

对公式(13) 两边取对数, 将求模型参数后验概率极大值转换为求目标函数极小值的过程, 得到:

(14)

(14) 式右边第一项控制弹性参数提取的精度, 第二项为弹性参数提取的正则化项, μ表示超参数, 控制弹性参数提取的“块化”效果, μ越大, 弹性参数“块化”效果越明显, 相反提取的弹性参数将出现明显的带限效果。

然而直接利用(14) 式提取弹性参数时, 会在一定程度上破坏其低频成分, 表现在弹性参数剖面上就会出现明显的“门帘”效应(由低频成分不准确引起的反演剖面横向不连续、数条类似门帘的现象)。

在(14) 式右边引入低频软约束项, 使提取的弹性参数的低频成分逼近测井曲线的低频成分, 确保弹性参数提取过程中低频成分准确, 因此:

(15)

其中,

(16)

式中:λ_P, λ_S和λ_ρ分别表示纵波速度、横波速度和密度的低频软约束权重系数; L_P, L_S, L_ρ分别代表对模型参数中的ln(v_P/v_P0)^T, ln(v_S/v_S0)^T, ln(ρ/ρ₀)^T起作用的低通滤波矩阵, 其具体表达式见附录A。m_{P_Trend}, m_{S_Trend}, m_{ρ_Trend}分别代表ln(v_P/v_P0)^T, ln(v_S/v_S0)^T, ln(ρ/ρ₀)^T的低频成分。对(15) 式两边求导并令其为0:

(17)

其中,

(18)

式中:[*]_xy表示取矩阵“*”第x行、第y列元素。由于非一致性加权矩阵Q中包含了待求解参数m, 因此, 需采用迭代重加权最小二乘算法迭代求解m。

为了检验本文算法的有效性和精度, 进行如下模型试算。在模型试算中, 利用如下反演目标函数来分别获取各个角度的弹性阻抗:

(19)

式中:r(θ)为角度反射系数; S(θ)表示入射角为θ的角度叠加数据; W(θ)表示角度子波; σ_r(θ)表示角度反射系数的均方差。

(19) 式右边由3部分组成, 分别是地震数据拟合差、基于单变量柯西分布先验约束的反射系数稀疏约束项以及低频软约束项。公式(15) 中也有低频软约束项, 但是意义和目的却不相同:(19) 式中的低频软约束项是为了在弹性阻抗反演过程中融合测井的低频成分, 因为地震资料中缺失0~10 Hzx的低频分量, 因此意义在于补充低频成分; 而在弹性阻抗分解提取弹性参数阶段, 弹性阻抗的低频成分已经融合了测井低频成分, 所以低频成分在理论上应该正确, 但是三变量柯西约束项会在一定程度上破坏该低频成分, 因此其意义在于修正和稳定低频成分。

由于弹性阻抗反演的模型参数是r(θ), 因此表达式也略有不同, 需要引入积分矩阵, 将角度反射系数通过积分转化为弹性阻抗。

(20)

式中:λ表示弹性阻抗低频软约束权重系数; C表示积分矩阵。

式中:E_Ii(θ), i=1, 2, …, N, 表示第i个采样点的弹性阻抗; smooth(*)表示对“*”进行低通滤波, 取其低频成分。

采用图 1所示的弹性参数模型结合Zoeppritz方程计算角度反射系数并和角度子波褶积得到合成角道集(图 2a), 其中角度子波均为30Hz的雷克子波, 其反射角分别是10°, 20°和30°。为了检验本文算法的抗噪性, 在合成角道集中加入了信噪比为2的随机噪声(图 2b)。

图 1(Fig. 1)

图 1 弹性参数模型

图 2(Fig. 2)

图 2 不含噪(a)和信噪比为2(b)的合成角道集

对图 2a所示角道集采用公式(19) 进行弹性阻抗反演, 结果如图 3所示。可以看出, 弹性阻抗的反演结果非常精确。在此基础上, 利用基于双项约束的弹性阻抗分解方法提取弹性参数。图 4给出了提取的弹性参数结果, 其中黑实线代表模型数据, 红实线代表反演结果, 从两者的对比可以看出, 提取的弹性参数与模型的吻合度相对较好。为了检验提取的弹性参数低频成分的准确性, 图 4中用黑虚线和蓝虚线分别表示模型的低频成分和提取弹性参数的低频成分, 可以看出两者基本重合, 说明获取的低频成分完全正确。图 5给出了仅基于三变量柯西约束的弹性阻抗分解方法提取的弹性参数, 可以看出黑虚线和蓝虚线已不再重合, 说明仅利用三变量柯西约束项时会在一定程度上破坏低频成分。图 6和图 7分别给出了迭代1次和迭代5次提取的三参数反射率与模型的对比。可以看出, 随着迭代次数的增加, 三变量柯西约束在迭代过程中会使得三参数反射率逐渐“尖脉冲化”, 相应的弹性参数也就逐渐“块化”。

图 3(Fig. 3)

图 3 弹性阻抗反演结果 a 小角度; b 中角度; c 大角度

图 4(Fig. 4)

图 4 基于双项约束的弹性参数提取结果

图 5(Fig. 5)

图 5 仅基于三变量柯西约束的弹性参数提取结果

图 6(Fig. 6)

图 6 迭代1次提取的弹性参数反射率与模型数据

图 7(Fig. 7)

图 7 迭代5次提取的弹性参数反射率与模型数据

图 8给出了含噪声合成角道集的弹性阻抗反演结果。在此基础上, 图 9给出了弹性参数的提取结果。对比模型与反演结果可以看出, 弹性阻抗的反演结果和后续提取的弹性参数与井曲线吻合度均较好, 说明了本文的弹性阻抗反演算法的抗噪性较好。

图 8(Fig. 8)

图 8 含噪声合成角道集的弹性阻抗反演结果 a 小角度; b 中角度; c 大角度

图 9(Fig. 9)

图 9 含噪声合成角道集的弹性参数提取结果

实际资料选自海上某工区, 目的层为薄互层砂泥岩储层。经过前期真振幅处理后, 对共反射角道集进行角度部分叠加, 分别形成小角度、中角度和大角度部分叠加剖面, 如图 10到图 12所示。根据目标函数公式(19), 角度反射系数的先验分布——单变量柯西分布对应的尺度参数σ利用EM算法从井曲线中计算获取, 低频软约束项中的弹性阻抗低频趋势通过对井曲线计算的弹性阻抗进行10Hz的低通滤波获取。

图 10(Fig. 10)

图 10 小角度部分叠加剖面

图 11(Fig. 11)

图 11 中角度部分叠加剖面

图 12(Fig. 12)

图 12 大角度部分叠加剖面

经过井震标定, 并提取角度子波后, 分别反演获得3个角度的弹性阻抗, 其结果如图 13到图 15所示。对图 13到图 15的结果进行仅基于三变量柯西分布约束的弹性阻抗分解, 提取的弹性参数如图 16到图 18所示。可以看出, 单项约束无法保证提取的弹性参数剖面横向连续, P波速度、S波速度以及密度剖面均出现了较为明显的竖条纹——“门帘”现象。图 19, 图 20, 图 21给出了基于双项约束提取的弹性参数, 由于低频软约束项的引入, 消除了图 16到图 18中的竖条纹, 剖面的横向连续性得到了有效增强。

图 13(Fig. 13)

图 13 小角度弹性阻抗

图 14(Fig. 14)

图 14 中角度弹性阻抗

图 15(Fig. 15)

图 15 大角度弹性阻抗

图 16(Fig. 16)

图 16 仅基于柯西分布约束的P波速度提取结果

图 17(Fig. 17)

图 17 仅基于柯西分布约束的S波速度提取结果

图 18(Fig. 18)

图 18 仅基于柯西分布约束的密度提取结果

图 19(Fig. 19)

图 19 基于双项约束的P波速度提取结果

图 20(Fig. 20)

图 20 基于双项约束的S波速度提取结果

图 21(Fig. 21)

图 21 基于双项约束的密度提取结果

利用纵横波速度比和P波阻抗交会图可以识别储层岩性或含油气性, 因此解释人员更为关心纵横波速度比的反演结果。该参数反演的精度, 可以作为叠前反演算法稳健性的一个重要指标。图 22给出了纵横波速度比的反演结果, 从剖面上看, 横向变化连续且丰富, 可有效揭示储层的岩性和含油气性变化。

图 22(Fig. 22)

图 22 基于双项约束的纵横波速度比提取结果

为了进一步对比井曲线和反演结果, 图 23和图 24分别给出了未参与反演的A井井旁道弹性阻抗与弹性参数反演结果与井曲线的对比, 其中黑线为经过250Hz Backus滤波后, 再重采样到时间域的井曲线, 红线则表示井旁道反演结果。可以看出, 井旁道弹性阻抗和弹性参数的反演结果与井曲线吻合度均较高, 可以较为有效地揭示地下储层弹性参数的垂向变化, 说明了该弹性阻抗反演算法精度较高。

图 23(Fig. 23)

图 23 井旁道弹性阻抗反演结果(红线)与测井曲线(黑线)对比 a 小角度; b 中角度; c 大角度

图 24(Fig. 24)

图 24 井旁道弹性参数提取结果(红线)与测井曲线(黑线)对比

针对常规弹性阻抗分解算法在存在噪声时不稳定这一问题, 本文介绍了利用双项约束来稳定弹性阻抗分解, 提高了弹性参数提取过程的鲁棒性。

三变量柯西约束项为稀疏约束项, 通过结合一阶差分矩阵将模型参数转换至三参数反射率, 利用迭代算法可以有效“尖脉冲”化三参数反射率, 提高岩性分界面的识别能力。

低频软约束项通过使得反演低频逼近测井低频来恢复和稳定三参数的低频趋势, 可以有效消除弹性参数剖面的“门帘”现象, 使之横向更加连续, 因此具有明确的地质意义。

模型测试和实际数据应用表明, 基于双项约束的弹性阻抗分解算法具有较高的稳定性和精度。但是低频软约束项是借助离散傅里叶变换实现低通滤波的矩阵表达, 因此计算效率不高, 需进一步优化算法。

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

(A1)

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

(A2)

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

(A3)

(A4)

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