在石油地球物理勘探中, 深层偏移成像效果一直受限, 其中主要原因就在于地层对地震波的吸收衰减效应。由于地下介质本身具有十分强烈的粘性特征, 而传统偏移成像方法以胡克定律作为地下的岩石物理假设, 与实际情况存在出入, 因而导致成像结果出现偏差。

逆时偏移成像^[1-3]以双程波动方程为理论基础, 大幅度提高了叠前深度偏移结果的精确性和可信度, 但是, 由于早期所采用的波动方程均为没有考虑地下吸收衰减特性的声波或者弹性波方程, 在实际应用中, 该方法成像结果的分辨率与精确度不是特别理想^[4]。

补偿因地层吸收而衰减的能量主要有两种方法。第一种是反Q滤波^[5-6], 在偏移成像前直接对地震道进行反Q粘性补偿, 该方法无论是在时域还是在频域都可以在一定程度上直接提高地震资料的分辨率。然而, 振幅的衰减与相位的频散是随地震波的传播逐渐变化的, 直接采用反Q滤波对地震道本身做补偿, 忽视了地震波传播的几何特征^[7]。第二种是在偏移成像过程中, 对延拓的波场进行粘性补偿, 采用具有粘性特征的波动方程进行延拓^[8-10], 其中, 粘性介质逆时偏移^[11-13]由于其较高的精度, 深受学术界与工业界的关注。

有两种方程被广泛应用于粘性介质逆时偏移。其一是基于标准线性体模型的带有记忆变量的波动方程^[14], 其二是基于常Q模型的含有分数阶拉普拉斯算子的波动方程^[14-16]。前者具有计算效率高、编程简便、易于实现等优势, 只要将记忆变量的符号变负, 便可实现对成像结果的补偿^[17-18]。但是其振幅衰减与相位频散耦合, 改变记忆变量的符号将会导致震源波场与检波点波场传播速度出现差异, 使得成像结果同相轴发生错位^[19-20]。而后者虽然需要利用傅里叶变换等方法求解分数阶拉普拉斯算子, 增加了计算量, 但是由于其振幅衰减项与相位频散项是分开的, 仅改变振幅衰减项的符号, 因而不会影响到震源波场与检波点波场的相位特征, 传播速度一致, 保证了成像结果中同相轴位置的准确性^{[12, 21-22]}。

然而, 基于分数阶拉普拉斯算子的波动方程在进行振幅补偿的过程中, 由于其补偿时间算子的补偿比例随着波数的增加呈指数上升, 相当于一种高通滤波, 会导致在延拓过程中高频噪声比例持续上升, 最终淹没有效信号^{[12, 19, 23]}。

对于高频噪声的压制, 最初采用的方法是低通滤波, 包括巴特沃斯低通滤波器^[13]、Turkey窗低通滤波器^[19]等, 但是, 由于这些滤波器滤波参数固定, 在较长延拓时间后仍然会有部分高频噪声残留。许多学者对波场稳定性补偿传播进行了研究, XIE等^[24]提出了稳定性因子法, 分别进行未衰减延拓以及衰减延拓, 然后计算稳定性补偿因子, 从而进行波场的精确补偿。田坤等^[25]基于非相位频散粘声方程, 提出了一种添加正则项的稳定性传播方程, 避免了补偿不稳定问题。SUN等^[26]将稳定性补偿因子法拓展到成像条件, 分别提出了两种稳定性补偿成像条件。GUO等^[27]将最小二乘理论中的Born近似理论用于计算粘性反传伴随方程, 由于其具有衰减特性, 不存在不稳定问题, 然后根据反偏移计算等步骤获得经过补偿的成像结果。WANG等^[23]借鉴反Q滤波的思想, 提出了一种自适应稳定传播算子。

本文提出了一种稳定性补偿策略, 与田坤等^[25]所提出的方法相似, 均属于正则化方法, 不同的是, 本文所采用的粘声方程同时考虑了振幅衰减与相位频散。首先介绍了考虑吸收衰减效应的逆时偏移补偿原理, 然后推导出正则化的稳定传播的粘声波动方程, 最后进行两种模型的逆时偏移测试, 对相位与振幅解耦的成像结果进行了比较, 进一步验证了本文方法的实用性。

粘声介质逆时偏移的实现过程与声波逆时偏移类似:①通过子波激发延拓震源波场; ②加载观测地震记录延拓检波点波场; ③根据时间一致性原理, 将震源波场与检波点波场进行相关计算获得偏移成像结果。与声波逆时偏移不同的地方在于, 粘声介质逆时偏移要求在延拓的过程中进行波场补偿。偏移原理如图 1所示。

图 1(Fig. 1)

图 1 粘声介质逆时偏移原理 a理想情况下地震数据的接收; b实际情况下地震数据的接收; c理想情况下的逆时偏移成像; d实际情况下不考虑地震数据衰减的逆时偏移成像; e实际情况下考虑地震数据衰减的逆时偏移成像

图 1a是在无吸收时所接收的理想地震记录, 其中e⁰代表无吸收衰减传播; 图 1b是实际情况接收地震记录, 其中e^－Lt为衰减时间算子, 代表波场在从激发到接收的过程中, 共经历了下行与反射后上行两个阶段的吸收衰减; 图 1c是对无吸收理想地震记录进行声波逆时偏移的原理, 其得到的成像结果是完全理想的结果; 图 1d是实际地震记录的声波逆时偏移原理, 其得到的成像结果会受到地震记录本身所携带衰减特征的不利影响, 从而导致最终成像结果存在衰减现象; 图 1e是实际地震数据的粘声补偿逆时偏移原理, 由于考虑了地震记录在接收过程中的衰减特性, 采用与接收时一致的地下参数所计算出的补偿时间算子e^+Lt进行波场延拓补偿, 从而可以在理论上得到与理想状态下高度接近的逆时偏移结果^[19]。

基于以上原理, 本文采用震源波场与检波点波场共同补偿的互相关成像条件^[28]:

式中:I(x, z)为最终成像结果; T为地震记录的最大时间; S^comp(x, z, t), R^comp(x, z, t)分别是经过补偿的震源波场与检波点波场。为了能够提高深层构造的成像质量, 可以在逆时偏移后进行照明补偿, 通过计算声波情况下震源波场的自相关得到照明结果, 将互相关成像结果与照明结果相除, 便可以得到深浅能量更加均衡的偏移剖面。

基于常Q模型的分数阶拉普拉斯算子粘声波动方程^[17]表示为如下形式:

式中:u为声波波场; γ=arctan(1/πQ); c=c₀cos(πγ/2);η=－c₀^2γω₀^－2γcos(πγ); τ=－c₀^2γ－1ω₀^－2γsin(πγ); Q为品质因子; c₀为定义在参考频率ω₀处的速度。

方程(2)等号右侧第一项为相位延迟项, 若将其更改为拉普拉斯算子, 即▽²u, 则变成了不考虑相位延迟的方程; 第二项为振幅衰减项, 若将其去掉, 则变成了不考虑振幅衰减的方程。为了能够在延拓过程中补偿振幅而使相位不受影响, 则只需要改变第二项的符号, 便可以得到振幅补偿的粘声波动方程:

采用方程(3)进行震源波场以及检波点波场的延拓, 便可得到吸收补偿的逆时偏移结果。

对于方程(2)和方程(3)中的分数阶拉普拉斯算子, 可以采用波数域求取的方式^[29]进行求解:

式中:F和F^－1分别代表傅里叶正变换和逆变换; k为波数。

粘声吸收作用于波场的比例是伴随着波数增加而呈自然指数递减, 而方程(3)中的补偿作用比例伴随着波数增加呈自然指数递增, 这将会使波场中的高频成分被严重放大, 从而导致数值计算不稳定, 直至淹没有效信号的低频成分^[23]。所以, 选择一种使补偿传播变得稳定的方法至关重要。

基于在求解对流问题偏微分方程中常用的预测-校正算法(predictor-corrector algorithm, PCA)策略^[30-31], 本文提出如下稳定性补偿粘声传播方程:

式中:w为无物理意义的中间变量; Δt为时间采样间隔; α为正则化参数, 取值范围在0~1, 如果α=0, (5)式退化为正常形式的粘声补偿波动方程。对于正则化参数的选择, 田坤等^[25]以及QU等^[32]认为可以通过模型的截止频率和最大速度计算得到, 参照这两篇文献, 本文将正则化参数设定为0.02。在方程(5)中, η(－▽²)^γ+1控制着相位延迟特征, 若修改为拉普拉斯算子▽², 则不考虑相位延迟; $\tau\left(-\nabla^{2}\right)^{\gamma+\frac{1}{2}}$控制着振幅衰减特征, 若将其去掉, 则不考虑振幅衰减。

PCA策略最初被运用于求解对流问题, 其后WANG等^[33]将其引入地震波传播问题, 用于解决位于高波数的地震波数值频散问题; ZHOU等^[34]将其引入常Q模型粘声地震波重建中, 但是采用的是非相位频散粘声方程, 并将声波情况下的正则化项直接引入, 没有进行重新推导。方程(5)基于相位频散振幅衰减粘声方程进行推导, 带有正则化参数α的项便为正则化项, 同时考虑了相位与振幅两方面的粘声特征。

首先进行简单的水平层状模型测试, 同时对于方程的相位与振幅解耦特征进行成像说明。模型网格大小为300 m×300 m, 水平和垂直方向的网格间距均为10 m, 介质参数如图 2所示。从0开始在地表均匀布置震源, 炮间距为150 m, 共20炮, 采用检波点随炮移动接收的方式, 检波点间距为10 m, 共100个, 炮两侧各放置50个进行接收, 时间采样间隔为1 ms, 记录长度为2.5 s。正演地震记录时同时考虑地下吸收衰减导致的相位与振幅变化。成像结果如图 3所示。

图 2(Fig. 2)

图 2 水平层状模型

图 3(Fig. 3)

图 3 水平层状模型成像结果 a声波数据声波逆时偏移结果; b粘声数据声波逆时偏移结果; c粘声数据粘声逆时偏移结果

由图 3可以看出, 使用本文方法进行粘声数据的补偿成像, 可以达到利用理想状态下声波数据获得的成像结果同样的效果。图 4为抽取水平方向1.5 km处的振幅曲线, 可以明显看出, 相对于非补偿成像, 粘声补偿成像的结果不仅相位校正到正确的位置, 同时振幅能量也得到了十分良好的补偿。将抽取的振幅曲线进行傅里叶变换, 得到如图 5所示的振幅谱, 可以明显看出, 补偿成像结果的频带相对于非补偿结果有了明显的补偿与展宽, 可以达到理想状态下声波数据声波成像的频带水平。

图 4(Fig. 4)

图 4 水平方向1.5 km处振幅曲线

图 5(Fig. 5)

图 5 水平方向1.5 km处曲线振幅谱

同时, 本文还进行了相位与振幅解耦成像测试, 比较结果见图 6。可以看出, 如果仅在方程中考虑相位问题, 成像结果的同相轴将会被校正到正确的位置, 但振幅却无法得到补偿; 而仅考虑振幅问题, 则与文献[25]中的情况一样, 虽然振幅得到了很好的补偿, 但是同相轴位置却出现了错位。可见, 只有同时考虑两种因素, 成像结果才能够在振幅水平与同相轴位置方面同时达到精确校正与补偿。

图 6(Fig. 6)

图 6 相位与振幅解耦成像结果比较

为了进一步验证本文方法在复杂情形中的适用性, 从Marmousi模型中截取一部分进行成像测试, 模型参数如图 7所示。模型网格大小为400 m×300 m, 横向以及纵向网格间距均为10 m, 从0开始在地表均匀布置震源, 炮间距为200 m, 共20炮。采用检波点随炮移动接收的方式, 检波点间距为10 m, 共100个, 炮两侧各放置50个进行接收。时间采样间隔为1 ms, 记录长度为3 s。成像结果如图 8所示。

图 7(Fig. 7)

图 7 局部Marmousi模型 a速度模型; b品质因子模型

图 8(Fig. 8)

图 8 局部Marmousi模型成像结果 a声波数据声波逆时偏移结果; b粘声数据声波逆时偏移结果; c粘声数据粘声逆时偏移结果

由图 8可以看出, 在复杂模型的情况下, 本文方法依旧可以得到理想的成像结果, 相对于粘声数据声波逆时偏移结果, 本文方法在保证了浅层成像效果的同时, 深层成像更加清晰, 由水平方向2 km处抽取出的振幅曲线(图 9)也可以进一步看出本文方法对于深层成像的优势。求取所抽取曲线的振幅谱, 结果如图 10所示, 可以看出, 本文方法对于复杂模型的频谱也有极佳的恢复与展宽效果。

图 9(Fig. 9)

图 9 水平方向2 km处振幅曲线

图 10(Fig. 10)

图 10 水平方向2 km处振幅曲线的振幅谱

本文提出了一种正则化的稳定传播的粘声补偿方程, 通过理论推导以及模型测试证明了本文方法的精确性和有效性, 得出如下结论和认识。

1) 通过对方程添加正则化项或者进行正则化处理, 可以有效避免补偿过程中的数值不稳定问题, 保证了补偿逆时偏移的有效性;

2) 地层吸收衰减会同时影响地震波的振幅与相位特征, 只有在进行补偿成像的过程中同时考虑这两种因素, 才能保证最终成像结果的精确性与可信度;

3) 采用互相关成像条件进行粘声补偿成像, 需要对震源波场以及检波点波场都进行补偿处理, 得到的结果无论是振幅还是频带都能够得到良好的恢复。

对于正则化参数的选取, 本文暂不能提出一种量化的方式, 这也是在今后需要继续探索的方向。