ROBINSON等^[1-2]提出了地震勘探中重要的褶积以及反褶积模型, 即自激自收的地震波场近似相当于子波与地层反射系数的褶积。考虑到地震信号的带限特征, 传统反褶积只能在有限的频带内拓宽频带, 因此, 假设地下模型参数稀疏特征的稀疏地震反演方法应运而生。针对不同变换域中的稀疏地震反演分类, 地震反演可以分为时间域地震反演和频率域地震反演(即谱反演)。

CASTAGNA等^[3]和PORTNIAGUINE等^[4]在谱分解的基础上, 在频率域损失函数的约束下, 发展了一种仅仅采用部分频谱分量去反演稀疏反射系数和层厚度的谱反演方法。CHOPRA等^[5]应用谱分解得到的局部频谱资料反演薄层反射系数, 并且认为在不需要任何先验假设的情况下可分辨调谐厚度以内的薄层。YUAN等^[6]基于贝叶斯学习理论发展了频率域稀疏反射系数反演方法, 并通过模型测试证明其可以分辨小于调谐尺度的薄层。CHAI等^[7]基于非平稳褶积模型发展了基追踪频率域地震高分辨率反射系数反演方法。目前, 基于褶积模型的频率域地震反演仍然聚焦于高分辨率反射系数反演方面, 由于频率域反演缺少背景模型约束信息, 仍存在稳定性差且对噪声敏感等问题。

相对于频率域地震反演来说, 时间域地震反演发展更为迅速, 利用模型参数先验信息和正则化约束可以有效缩小模型参数求解空间。近年来, 匹配追踪和基追踪算法成为高分辨率稀疏反演的研究热点^[8-11], 块化地层假设有助于恢复宽频带的反射系数序列。杨培杰^[12]根据贝叶斯参数估计理论, 并假设模型参数服从柯西分布, 发展了基于高阶统计量的叠后地震盲反演方法。张繁昌等^[13]从贝叶斯稀疏反演理论出发, 并结合波阻抗先验约束信息, 发展了修正柯西约束地震盲反演方法。BULAND等^[14]研究了贝叶斯线性AVO/AVA反演方法, 直接反演得到纵、横波速度及密度信息。ZONG等^[15]和YIN等^[16-17]发展了基于平滑模型约束的贝叶斯叠前AVA反演方法, 提高了反演算法的稳定性和预测结果的可靠性。在储层物性参数和流体因子时间域反演方面, 胡华锋等^[18]和张世鑫等^[19]分别研究了贝叶斯分类储层物性联合反演方法和基于孔隙弹性理论的Russell流体因子贝叶斯直接反演方法。

由于时间域地震反演相比频率域反演的地层分辨率有所降低, 考虑到时间域反演和频率域反演在分辨率和抗噪性方面独有的特征, 联合贝叶斯稀疏反演框架、时间域反演理论及频率域稀疏反演理论, 提出了时频联合域贝叶斯地震反演方法, 并引入平滑模型约束准则来降低模型空间的多解性。模型测试表明, 时频联合域稀疏反演在地震资料高信噪比的情况下保留了频率域反演分辨率高的优良特性, 同时在低信噪比的情况下延续了时间域反演抗噪能力强的优势, 在提高地震资料分辨率和压制弱小随机噪声之间达到平衡。实际资料处理验证了该方法的稳定性和可靠性。

Robinson褶积模型^[1-2]可以表示为地层反射系数序列与子波核矩阵的乘积, 即:

(1)

式中:d为实际观测地震数据; G为不同时移地震子波构建的核矩阵; m为地层模型参数组成的列向量; n为非有效反射信息组成的向量。地震反演的目的就是在包含噪声的带限地震记录中恢复地层稀疏反射系数, 因此, 求解公式(1) 的最直观求解方法可以表示为:

(2)

由于地震波场仅仅为有限频带支撑, 带限性质是造成反演多解性的主要原因。方程(2) 最小二乘估计(LSQ)的求解将表现出不稳定、可靠性低及抗噪能力差等缺点。为了解决LSQ地震反演存在的众多问题, 正则化约束准则被率先应用于不适定问题的求解过程^[9-15], 归纳其主要思想是通过将模型参数的先验信息加入反演目标函数, 进一步缩小模型参数的求解空间, 从而得到符合地下介质模型的可靠解。

针对贝叶斯参数估计的时间域反演理论, 即考虑地下介质反射系数的稀疏性度量准则^[12-15], 首先将某种确定性概率密度分布赋给似然函数和先验分布, 然后通过求取最大后验概率分布(MAP)来实现最优解估计。贝叶斯地震反演中的“长尾巴状”先验概率分布主要用来提高反演过程的抗噪性和稀疏性, 如柯西分布、Huber分布、t-分布、指数分布等^{[14, 20-27]}。从本质上分析, 先验分布约束从数学意义上可以理解为一系列模型参数的重加权正则化方法, 此时小的模型参数对应大的正则化系数, 大的模型参数对应小的正则化系数^{[12, 20, 28-30]}, 进而压制弱小反射突出强反射系数。其中, 贝叶斯框架下反问题具体表示为:

(3)

式中:J(m)为贝叶斯框架下的目标泛函; J_fitting(m)=‖d-Gm‖₂²为实际地震数据和模型正演数据的L₂范数残差度量(似然函数引入); J_prior(m)为待反演参数的先验约束项; λ₁, λ₂分别为似然函数和模型参数先验分布在贝叶斯估计中所占权重。研究结果表明, 基于贝叶斯反演框架的稀疏脉冲反演方法具有较高的稳定性和可靠性^{[15, 22]}。

图 1为时间域Robinson平稳褶积模型。时间域反演即在满足模型参数先验信息的前提下搜索与时间域带限数据(图 1a)相匹配的最优反射系数模型(图 1c)。

图 1(Fig. 1)

图 1 时间域平稳地震Robinson褶积模型 a地震数据; b子波核矩阵; c反射系数

将Robinson褶积模型变换到频率域, 得到对应的表达式为:

(4)

式中:ω为角频率; W(ω)为地震子波ω(t)的频谱; R(ω)为反射系数时间序列r(t)的频谱; N(ω)为背景噪声n(t)的频谱。为了更直观表达且适合讨论反演问题, 展开(4) 式中反射系数频谱:

(5)

式中:r(τ)为介质反射系数; exp(·)为自然指数运算。将(5) 式写成矩阵形式:

(6)

式中:r_k(k=1, 2, …, N)为N个反射系数; ω_i(i=1, 2, …, m)为参与反演的m个有效频率成分。将公式(6) 进一步改写成矩阵的形式:

(7)

式中:S=[S(ω₁) S(ω₂) … S(ω_m)]^T; m为待反演模型参数; G为频率域复数正演算子。由于在频率域正演算子构建过程中, 应考虑各频率成分信噪比(S/N)的影响, 选择具有较高信噪比的频率段(ω_low, ω_high)参与反演, ω_low和ω_high分别代表最低、最高截止频率。

由于复数域讨论反问题不够简洁, 将公式(7) 中的复数矩阵展开为实部和虚部表达式:

(8)

令S′=[Re(S) Im(S)]^T, G′=[Re(G) Im(G)]^T, 则公式(8) 可以改写为:

(9)

由于频率域地震反演是利用带限的地震频率域响应恢复地层真实反射系数序列, 通常在没有任何先验信息的前提下求解公式(9) 是不适定的。同样, 根据贝叶斯反演框架可以将模型参数的稀疏性先验概率分布作为假设条件, 得到贝叶斯反演的统一表达式:

(10)

式中:λ₁和λ₂为权衡地震频谱和先验信息反演权重的系数。图 2为平稳褶积的混合域示意图。混合域模型(图 2b)即在正则化过程中搜索拟合带限地震频谱(图 2a)的最优时间域反射系数(图 2c), 其在时频联合域反演过程中起到了重要的桥梁作用, 将时间域数据、频率域响应和时间域反射系数连接起来。

图 2(Fig. 2)

图 2 平稳褶积模型的混合域表达 a数据频谱; b混合域核矩阵; c反射系数

为充分利用带限地震信号的时间域和频率域信息, 将地震道集经傅里叶分析由时间域变换到频率域, 然后综合时间域和频率域响应特征发展基于贝叶斯框架下的时间-频率域联合反演方法。根据贝叶斯反演理论可知, 似然函数主要是通过背景噪声的概率密度分布来表示, 此外, 假设时间域地震数据和频率域地震数据是互相独立的随机变量, 则独立变量的联合概率分布(似然函数)可表示为:

(11)

式中:d为时间域地震数据; S′为频率域地震数据; m为待反演模型参数; p(d, S′|m)表示反演结果与时频联合域地震数据的匹配程度; p(d|m)代表反演结果与实际时间域地震数据之间的匹配程度; p(S′|m)代表反演结果与频率域地震数据的匹配程度。

假设时间域地震数据d和频率域地震数据S′的似然函数均符合均值为零的高斯概率分布, 则公式(11) 中的各成分可分别表示为:

(12)

(13)

式中:σ_n1²和σ_n2²分别为时间域和频率域背景噪声分布的方差; N, m分别表示时间域地震数据的长度和频率域正演矩阵所包含的频率分量个数。综合(12) 式和(13) 式可知, 时频域联合反演的似然函数p(d，S′|m)可以写成:

(14)

因此, 通过公式(14) 可以将时间域和频率域等不同变换域间的地震响应有效地联合在一起。

模型参数先验分布函数旨在描述待反演模型参数的先验信息, 基于柯西准则的先验概率分布算法稳定, 抗噪性能较强, 本文将先验分布函数表示为柯西准则:

(15)

式中:σ_m²为待反演模型参数m的方差。

贝叶斯理论建立了模型参数后验概率分布与先验分布和似然函数之间的关系, 其具体公式为:

(16)

联合(15) 式和(16) 式, 建立时频联合域反演中模型参数m的后验概率分布函数p(m|d, S′)为:

(17)

通过最大化(17) 式中的联合反演后验概率分布函数p(m|d, S′), 得到等价的目标泛函J(m):

(18)

式中:J_G(m)为时频域联合反演中实际地震响应与正演合成道集差异信息的量度; J_Cauchy(m)为先验概率分布引入的稀疏正则化项。

由于地震数据低频分量相对较弱, 故低频分量更容易受强噪声的干扰。因此实际工区中需要引入先验信息(如井资料)进一步对模型参数施加约束, 进而有效补充地震资料的低频缺失, 提高反演结果的稳定性及横向连续性。针对时频联合域反演方法加入平滑纵波阻抗约束, (18) 式中目标函数J(m)可表示为:

(19)

其中,

(20)

式中:C_smooth代表平滑模型与模型反演结果的匹配程度; λ₀为模型约束系数; 为积分矩阵; η₀为相对平滑纵波阻抗模型。为最小化目标泛函J(m), 求取(19) 式中J(m)对模型参数m的梯度, 并令其为0, 即可得模型参数非线性反演方程。

为验证时频联合域反演方法的可行性和稳定性, 我们根据实际工区中A井设计了一维地震模型进行验证, 采用主频为30Hz的零相位Ricker子波正演得到地震道集, 采样频率为500Hz, 如图 3所示。

图 3(Fig. 3)

图 3 A井合成地震记录 a时间域地震道集; b合成地震频谱分析

对A井合成地震记录添加不同比例的高斯白噪干扰, 其合成地震道集及频谱如图 3所示。图 4, 图 5和图 6分别为时间域、频率域和时频联合域中含不同噪声时的纵波阻抗反演结果, 蓝线为理论模型, 绿线为低频初始平滑模型, 红线为反演结果。分析无噪声干扰情况下的反演结果可知, 时间域反演结果(图 4a)在无噪情况下分辨薄层和高频弱小反射上仍存在误差, 通过控制联合域中频带参数的选择, 频率域(图 5a)和时频联合域(图 6a)反演结果在高频弱小反射位置相比时间域反演结果均有所改善。

图 4(Fig. 4)

图 4 时间域中含不同噪声时的纵波阻抗反演结果 a无噪声; b 10%噪声; c 20%噪声; d 50%噪声

图 5(Fig. 5)

图 5 频率域中含不同噪声时的纵波阻抗反演结果 a无噪声; b 10%噪声; c 20%噪声; d 50%噪声

图 6(Fig. 6)

图 6 时频联合域含不同噪声时的纵波阻抗反演结果 a无噪声; b 10%噪声; c 20%噪声; d 50%噪声

为了更清楚地明确反演结果的分辨率差异, 本文采用时频拟合准则来综合评价时频联合域反演和时间域反演的分辨率特征, 其中时频谱局部拟合误差F_tf-loc和全局拟合误差F_tf-glob可分别定义为:

(21)

(22)

式中:|W(t, f)|为反演信号时频谱; |W_r(t, f)|为理论信号时频谱。对无噪声情况下的反演结果进行高分辨率匹配追踪(MP)时频谱分析, 结果如图 7所示。由图 4和图 6可知, 反演结果即使在时间域存在微小差别的情况下, 其时频拟合误差的差别仍然很大(图 7b, 图 7c)。时间域反演结果在恢复反射系数高频段能量时存在较大的误差, 而时频联合域反演却可以较好地恢复反射系数高频段能量(图 7b和图 7c中椭圆区域)。

图 7(Fig. 7)

图 7 理论纵波阻抗反射系数时频谱(a)、时间域反演的时频全局拟合误差Ftf-glob(b)和时频联合域反演的时频全局拟合误差Ftf-glob(c)

随着地震记录信噪比的不断降低, 反演分辨高频弱小反射层的能力均有不同程度的降低(图 4b, 图 4c, 图 4d, 图 5b, 图 5c, 图 5d, 图 6b, 图 6c, 图 6d)。时间域反演结果在分辨率上虽有较大幅度的降低, 但是在20%噪声的情况下时间域反演结果仍然保留了相对可靠的纵波阻抗变化趋势; 此外, 从图 5b, 图 5c, 图 5d可以看出, 在含有50%噪声时, 频率域反演的不稳定性就会突显(但仍然保持较高的分辨特性), 反演结果的可靠性大大降低; 建立在此基础上的时频联合域反演方法不但继承了时间域反演抗噪性强的优势, 并且由于在频率域中各频率分量实现自动解耦, 通过在频率域中选择具有较高信噪比的部分有效频率成分作为时频联合域反演的输入, 从而保证了频率域反演的高分辨率特性。

为了讨论本方法中时间域和频率域地震响应的权重系数对反演精度的影响, 对时间域贝叶斯反演、频率域贝叶斯反演及时频联合域反演中的权重系数进行分析, 结果如图 8所示。从图 8可以看出, 随着地震响应权重系数的不断增大, 3种方法的反演结果与理论模型的相关系数呈现先增大后减小的趋势。其中, 随着噪声含量的不断增加, 时间域反演表现出较强的稳定性(反演精度基本一致); 随着噪声含量和频率域响应权重的增加, 频率域反演却表现出反演精度剧烈下降的趋势, 进一步证实了频率域反演对噪声敏感的问题。值得注意的是, 在低噪声情况下, 频率域反演结果比时间域反演结果的精度高。本文方法建立在时间域反演稳定性强的前提下, 改善了频率域反演对噪声敏感的问题, 尤其在强噪声情况下提高了反演结果的精度(图 8c矩形框)。

图 8(Fig. 8)

图 8 不同变换域权重系数对反演结果的影响 a时间域; b频率域; c时频联合域

为了进一步验证本文方法对二维地震数据反演的可行性和横向连续性, 根据某工区地质沉积背景抽象出二维地质沉积模型作为理论纵波阻抗模型(图 9), 其中包括河流相沉积的杂乱无章砂体分布(约3260ms附近)和下覆两套前积薄砂互层(约3290ms附近)。图 9中, 两条黑色曲线分别表示CDP135和CDP265处的纵波阻抗曲线。模型大小为横向500道, 采样时长180ms。采用30Hz零相位雷克子波正演得到合成地震记录, 如图 10a所示, 可见, 由于薄层的调谐作用模糊了地层界面信息, 造成合成道集中难以直观地识别出砂体形状以及砂体尖灭位置。

图 9(Fig. 9)

图 9 实际工区理论地质沉积纵波阻抗模型

图 10(Fig. 10)

图 10 理论地质沉积模型时频联合域纵波阻抗反演结果 a合成地震记录; b平滑模型; c无噪声时间域反演; d 10%噪声时间域反演; e无噪声联合域反演; f 10%噪声联合域反演

对图 10a所示地震道集开展常规时间域反演和时频联合域反演。为了提高反演的可靠性和稳定性, 在反演过程中添加初始低频平滑模型约束, 如图 10b所示。图 10c和图 10d分别为时间域无噪声和含10%噪声的纵波阻抗反演结果, 可见, 在含10%噪声的情况下, 纵波阻抗反演剖面与理论地层模型仍保持较好的一致性, 超低频模型约束的加入使得反演结果横向连续性较好。而在反演薄砂互层时却存在较大的误差, 如图 10c和图 10d中黑色椭圆区域所示的沉积砂体尖灭位置, 并且在河道砂沉积区也存在边界模糊问题。

时频联合域反演结果如图 10e和图 10f所示, 可见, 即使在含10%噪声的情况下, 该方法同样可以有效分辨出砂体尖灭位置(图中椭圆处), 并且河道砂边界模糊问题有了较大程度的改善, 有助于对优质砂体的展布范围进行精细刻画。

为了验证本文方法的实际应用效果, 将时频域联合地震反演应用于我国东部某探区的纵波阻抗反演, 结果如图 11所示, 图中, 绿线和黑线指示开发井W的实测纵波阻抗曲线, H1和H2为实际工区目的层位。图 11a和图 11b分别为经过精细处理的两条叠后地震过井剖面Ⅰ, Ⅱ。首先, 针对实际地震资料展开频谱分析, 将地震数据和子波通过Fourier变换获取其频率域响应, 进而提取相应频带范围的地震子波; 然后, 将不同变换域中的地震响应作为联合反演的输入信息开展基于贝叶斯理论的时频联合域反演。根据实际测井响应的超低频信息建立初始平滑模型如图 11c和图 11d所示。常规贝叶斯时间域反演得到的纵波阻抗如图 11e和图 11f所示, 利用本文方法得到的纵波阻抗如图 11g和图 11h所示。从图 11中可以看出, 时频联合域反演相比常规时间域反演具有更高的分辨率, 在局部细节部位有较好的辨识能力, 尤其在图 11中的椭圆区域, 与实际测井分层结果保持较高的一致性。此外, 时频联合域反演在保证分辨率的同时保障了反演结果的信噪比, 低频平滑模型的加入使得反演结果横向连续性较好。

图 11(Fig. 11)

图 11 实际地震资料联合域纵波阻抗反演结果 a过井线地震道集(Ⅰ); b过井线地震道集(Ⅱ); c过井线初始平滑模型(Ⅰ); d过井线初始平滑模型(Ⅱ); e过井线常规反演结果(Ⅰ); f过井线常规反演结果(Ⅱ); g过井线联合域反演结果(Ⅰ); h过井线联合域反演结果(Ⅱ)

不同变换域中的地震反演方法具有不同的抗噪性和分辨率特征, 从反演稳定性方面分析, 时间域地震反演要优于频率域地震反演; 从分辨能力来说, 时间域地震反演要逊色于频率域地震反演; 从对噪声的敏感程度来讲, 频率域地震反演要高于时间域反演。时频联合域贝叶斯反演方法充分利用了带限地震信号的时间域和频率域响应。在频率域中不同频率分量会实现自动解耦, 因此可以将不同频率分量的反问题独立约束求解。实际资料测试表明, 时频联合域反演方法可以在提高地震分辨率的同时保证反演结果的信噪比。总体上, 该方法继承了时间域反演抗噪能力强的优势, 且在高信噪比的情况下保留了频率域反演高分辨率的优良特性。时频联合域反演在不同信噪比的情况下参数选择尤为重要, 在地震资料信噪比较高的情况下, 应当首要考虑频率域响应的拟合程度, 而在地震资料信噪比较低的情况下, 应在反演过程中适当增加时间域地震响应的匹配权值。