目前, 中国西部地区已经成为油气勘探的重点地区之一, 但是在该地区油气勘探深度内普遍存在低速夹层, 例如四川盆地二叠系和奥陶系的低速泥页岩层、塔里木盆地侏罗系煤层等。一系列高低速地层组合产生大的波阻抗差异, 地震波信号在这些波阻抗差异大的界面之间来回反射, 使得层间多次波发育。由于上覆地层生成的层间多次波会掩盖下伏地层的一次反射波信号, 因而不利于地震数据准确成像, 会产生虚假构造, 影响解释的精度^[1]。层间多次波压制是陆上地震勘探资料处理最困难的环节之一。目前对层间多次波的生成机制和响应特征还不是很清楚, 因此利用数值模拟方法研究层间多次波非常必要。

地震波场数值模拟方法主要分为两大类, 一类以射线理论为基础, 另一类以波动方程为基础。其中以波动方程为基础的地震波正演模拟方法应用最为广泛, 包括有限差分法^[2-3]、有限元法^[4-5]、谱元法^[6]和间断有限元方法^[7-8]等。尽管这些数值模拟方法具有很大的优势, 如有限差分法计算效率高, 占用内存小; 有限元类方法对地形描述更准确, 计算精度高, 但它们模拟的是总地震波场信号, 无法将一次波和层间多次波信号单独分离, 缺乏灵活性, 对研究层间多次波的响应机理极其不方便。在实际应用中, 研究人员利用预测方式压制多次波, 如扩展的自由表面多次波压制(SRME)方法^[9-10]、拉东变换方法^[11]、虚同相轴方法^[12-13]等, 而不是直接研究多次波的生成机理特征。当前大多数学者对多次波的研究也主要集中在表面多次波, 如LIU等^[14]提出最小二乘表面多次波偏移方法, 叶月明等^[15]提出分阶表面多次波偏移成像, 李强等^[16]提出分阶表面多次波数值模拟方法, 石颖等^[17]基于波动理论预测三维表面多次波等。对层间多次波的讨论深度和广度却明显不如表面多次波。主要原因在于层间多次波的产生机理更加复杂, 其研究难度也比表面多次波更大。

层间多次波模拟方法可以将一次波和层间多次波直接分离, 且模拟波场的相位和振幅与地质属性保持一致。KENNETT^[18]最早提出基于反射率法的层间多次波模拟方法, 但该方法目前无法很好地适应介质的横向变化; COVEY等^[19]提出一种基于射线追踪的层间多次波模拟方法, 该方法运算效率较高, 能处理简单介质中的层间多次波模拟问题, 但无法应对复杂变化的模型; BERKHOUT^[20]提出一种基于波场延拓的全波场模拟方法, 弥补了有限差分等只能一次性模拟总波场的缺陷, 可以分阶次模拟波场, 且能够处理相对复杂的地质模型; 匡伟康等^[21]提出自适应可控层的层间多次波模拟与识别技术, 进一步提高了效率, 同时可忽略其它层位的影响, 可以针对特定层的层间多次波进行模拟; 戴晓峰等^[22]提出叠前反射率正演模拟层间多次波方法, 对川中地区的层间多次波进行了有效识别。

在上述研究的基础上, 我们进一步提出了一种层间多次波模拟方法, 即分别在不同维度下分阶模拟层间多次波, 并采用傅里叶有限差分方法作为波场延拓算子, 提高地震波的横向传播精度。多次波分阶的过程完全可控, 一次循环得到一次反射波, 二次循环得到一次反射波和一阶层间多次波。通过设置对应的迭代次数, 可得到对应阶次的层间多次波。循环次数越高, 波场信息越丰富。本文研究了层间多次波在一维、二维和三维模型中的传播规律, 并对比分析了叠前和叠后层间多次波地震数据与全波场数据的关系。

1 基本原理

层间多次波数值模拟方法可单独将不同阶次的层间多次波分离出来, 主要包含两个过程: 一个是层间多次波分阶过程, 满足反射和透射条件; 另一个是地震波传播过程, 基于波场延拓算子传播地震波。

1.1 层间多次波分阶模拟

全波场模拟(FWMod)方法最早由BERKHOUT^[20]提出, 该方法将地震数据看成是一种全波场数据, 由一次反射波和各阶层间多次波叠加而成(本文主要研究层间多次波, 假设已经消除了表面多次波的影响)。下面具体描述波场分阶模拟的基本原理。

根据地震波场的传播原理, 当地震波传播至不均匀界面时, 会生成反射波和透射波继续传播。利用单程波方程对地震波进行数值模拟, 地震波场按照深度传播方向可分为下行波场和上行波场。不同界面之间上行波和下行波通过波场延拓算子进行传递。在不同维度的模型中进行层间多次波模拟, 需要将模型按照深度方向分成一系列离散的规则子区域(图 1), 对于一维模型, 其分成一系列的点; 对于二维模型, 其分成一系列的直线; 对于三维模型, 其分成一系列的薄平板。上行波场或下行波场在横向上的传播是在波数域中完成。本文以三维模型为例论述地震波场在传播过程中的关系(图 2), 其中上标n表示第n个地层界面, 下标u表示上行波, 下标d表示下行波, 则A_dⁿ和A_uⁿ分别表示在地层界面z_n上的下行波场和上行波场, B_dⁿ和B_uⁿ分别表示在界面z_n下的下行波场和上行波场。在界面z_n上, 满足关系:

$ \boldsymbol{B}_{d}^{n} =\boldsymbol{T}_{d}^{n} \boldsymbol{A}_{d}^{n}+\boldsymbol{R}_{u}^{n} \boldsymbol{B}_{u}^{n} $

(1)

$ \boldsymbol{A}_{u}^{n} =\boldsymbol{T}_{u}^{n} \boldsymbol{B}_{u}^{n}+\boldsymbol{R}_{d}^{n} \boldsymbol{A}_{d}^{n} $

(2)

式中: T表示透射系数矩阵; R表示反射系数矩阵。引入W表示波场延拓算子, 在下行延拓的过程中, 界面z_n－1下方的下行波场B_d^n－1经过下行延拓到界面z_n上方, 得到入射波场A_dⁿ:

$ \boldsymbol{A}_{d}^{n}=\boldsymbol{W}_{d}\left(z_{n}, z_{n-1}\right) \boldsymbol{B}_{d}^{n-1} $

(3)

同时, 在上行延拓过程中, 界面z_n+1上方的上行波场A_uⁿ⁺¹经过上行延拓到界面z_n下方, 得到入射波场B_uⁿ:

$ \boldsymbol{B}_{u}^{n}=\boldsymbol{W}_{u}\left(z_{n}, z_{n+1}\right) \boldsymbol{A}_{u}^{n+1} $

(4)

单程波数值模拟以递归的方式完成波场传播。当波场传播至某一层界面时, 均要满足方程(1)至方程(4)的关系式, 如地震波从界面z_m传播至界面z_n时, 可以得到关于界面z_n的下行波和上行波表达式分别为:

$ \boldsymbol{A}_{d}^{n} =\boldsymbol{W}\left(z_{n}, z_{m}\right)\left(\boldsymbol{S}_{d}^{m}+\boldsymbol{R}_{d}^{m} \boldsymbol{A}_{d}^{m}+\boldsymbol{R}_{u}^{m} \boldsymbol{B}_{u}^{m}\right) $

(5)

$ \boldsymbol{B}_{u}^{n} =\boldsymbol{W}\left(z_{n}, z_{m}\right)\left(\boldsymbol{S}_{u}^{m}+\boldsymbol{R}_{d}^{m} \boldsymbol{A}_{d}^{m}+\boldsymbol{R}_{u}^{m} \boldsymbol{B}_{u}^{m}\right) $

(6)

式中: S表示一次源; R_d^mA_d^m+R_u^mB_u^m可以理解为在界面z_m激发的二次震源。在推导公式(5)和公式(6)过程中隐含了如下的条件:

$ \boldsymbol{T}_{d}^{m}=\boldsymbol{R}_{d}^{m}+\boldsymbol{I} $

(7)

$ \boldsymbol{T}_{u}^{m}=\boldsymbol{R}_{u}^{m}+\boldsymbol{I} $

(8)

公式(7)和公式(8)给出了透射矩阵与反射矩阵之间的关系, I表示单位对角矩阵。根据公式(5)和公式(6), 波场在不同的界面之间不断地传播和散射, 在地表被观测系统接收, 形成地震记录。

公式(5)和公式(6)仅仅表示在一次循环过程中的上行波场和下行波场计算。地震波场的分阶模拟是一个不断循环迭代的过程, 可通过控制其循环次数来获得对应阶次的波场。在本文中, 一次反射波被认为是0阶层间多次波。层间多次波模拟循环过程如图 3所示, 通过第一次循环, 模拟得到0阶层间多次波场(图 3a), 图中的圆点表示延拓过程中激发的二次源。每增加一次循环可获得更高一阶的层间多次波。一阶层间多次波的生成机制如图 3b中的绿色路径所示, 二阶层间多次波的生成机制如图 3c所示(蓝色路线)。若要模拟更高阶次的层间多次波, 则需要更多的循环迭代过程。

上面对三维情况下层间多次波的生成机理和特征作了详细介绍, 据此可实现层间多次波的分阶数值计算。对于一维和二维模型中的层间多次波模拟也可以按照类似的思路进行数值计算。进行层间多次波分阶模拟时, 波场延拓算子要保持好的稳定性, 其在深度域满足递推关系(图 4), 波场延拓算子表达式为:

$ W_{d}\left(z_{n}, z_{m}\right)=\prod\limits_{i=m}^{n-1} W\left(z_{i+1}, z_{i}\right) \ \ m<n $

(9)

$ W_{u}\left(z_{n}, z_{m}\right)=\prod\limits_{i=m}^{n+1} W\left(z_{i-1}, z_{i}\right) \ \ m>n $

(10)

公式(9)和公式(10)表明波场延拓算子的累乘传递性质。

1.2 傅里叶有限差分波场延拓算子

波场延拓算子的精度会影响数值模拟的精度, 因而其备受关注。本文选取傅里叶有限差分方法作为波场延拓算子, 该方法可以利用有限差分校正项提高地震波在横向上的传播精度。横向非均匀介质中的三维频率域单程波方程表达式为^[23]:

$ \frac{\partial \hat{p}(x, y, z ; \omega)}{\partial z}=\mathrm{i} k_{z} \hat{p}(x, y, z ; \omega) $

(11)

式中: $ {\hat p}$(x, y, z; ω)表示频率域的波场; ω是角频率; x, y和z分别代表空间中的三个坐标轴方向; i表示虚数单位; k_z表示均方根算子, 且k_z=$\sqrt {{\omega ^2}/{v^2} + {\partial ^2}/\partial {x^2} + {\partial ^2}/\partial {y^2}} $(v表示介质速度)。由于无法直接利用均方根算子对波场${\hat p} $(x, y, z; ω)进行延拓, 需要将均方根算子k_z设置在参考速度v₀中进行计算。令k_z=k_z0+(k_z－k_z0)(k_z0为参考速度中的均方根算子), 并将k_z和k_z0进行泰勒级数展开和使用连分式, 方程(11)可以进一步分解为以下三个串联方程:

$ \frac{\partial \hat{p}_{1}(x, y, z ; \omega)}{\partial z}=\mathrm{i} k_{z_{0}} \hat{p}(x, y, z ; \omega) $

(12)

$ \frac{\partial \hat{p}_{2}(x, y, z ; \omega)}{\partial z}=\mathrm{i} \omega \Delta l \hat{p}_{1}(x, y, z ; \omega) $

(13)

$ \frac{\partial \hat{p}(x, y, z ; \omega)}{\partial z}=\frac{b\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right)}{1+a\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right)} \hat{p}_{2}(x, y, z ; \omega) $

(14)

方程(12)表示相移算子; 方程(13)表示慢度修正项, 其中Δl表示慢度, 其值为Δl=1/v－1/v₀; 方程(14)是横向上强速度差异的修正项, a和b为常系数。实际上, 方程(12)、方程(13)和方程(14)可以等效为:

$ k_{z} \approx k_{z_{0}}+\omega \Delta l+\frac{b\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right)}{1+a\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right)} $

(15)

方程(15)是傅里叶有限差分延拓算子。利用波场延拓算子和迭代循环, 可以实现高精度层间多次波的分阶数值模拟。

2 数值实例

利用层间多次波模拟技术在不同维度的模型中进行数值模拟, 从而识别单炮记录中的层间多次波和一次反射波, 并对叠前和叠后的地震数据进行分析。

2.1 一维模型

首先, 构建一个一维模型(图 5a), 该模型含有4个强波阻抗差异反射界面, 2个速度反转区。模型密度为常数2.0g/cm³。由于存在强速度差异, 导致层间多次波发育, 且多次波振幅能量较强。使用主频30Hz的雷克子波作为激发震源, 时间采样间隔为5ms, 记录时间为3.5s。图 5b展示了总的波场信号, 图 5c至图 5f是0~3阶层间多次波信号。从模拟的结果中可以得到4个一次反射波以及各阶层间多次波, 一阶层间多次波的幅值量级与一次反射波大致相当, 而二阶和三阶层间多次波幅值与一次反射波相差大约10²倍。四阶以上的多次波相对于一次反射波幅值能量很弱。从走时来看, 一次反射波与多次波相互混杂, 低阶多次波与高阶多次波也相互混杂, 它们之间都没有明显的时间分界线。

2.2 二维模型

选取中国西北某地区的测井曲线进行二维建模, 模拟分析该地区的层间多次波特征。图 6为测井曲线和据此设计的速度模型。其中图 6a中蓝色曲线为实际的声波测井曲线, 黑色曲线为根据实际测井曲线构建的近似层状速度模型, 图 6b为根据图 6a中黑色曲线构造的二维层状模型, 模型大小为7.5km×12.5km。使用雷克子波作为激发震源。模拟数据道间距为0.02km, 最大偏移距为6.25km, 采用中间放炮两边接收的观测系统。本文方法模拟的炮集数据如图 7所示, 其中图 7a是模拟的总波场, 图 7b是分离的一次反射波, 图 7c是分离的一阶层间多次波。从图 7可以看到实现了一次反射波和一阶层间多次波的有效分离。同时为了研究层间多次波对地震数据处理的影响, 对模拟的地震资料利用Geoeast软件进行速度分析, 拾取的速度谱如图 8所示, 其中图 8a是一次反射波速度谱, 图 8b是总波场速度谱。对比图 8a和图 8b可知, 增加层间多次波后, 速度谱中产生了很多靠近一次反射波能量团的层间多次波能量团, 这些层间多次波能量团会影响一次波能量团的拾取, 从而影响最终的速度谱拾取精度。并且层间多次波的叠加速度小于一次反射波。这种现象在实际地震资料处理中经常出现。

为了进一步研究叠后地震数据之间的关系, 构建了一个弱横向非均匀介质的地堑模型, 模型的水平长度为9.6km, 深度为4.0km。在深度为1.8km处有一个凹陷构造, 具体的速度结构如图 9所示。假设介质密度为常数2.0g/cm³。记录时间长度为3.5s, 时间采样间隔为5ms, 模拟数据的炮间距为0.02km, 最大偏移距4.3km, 总共模拟481炮。本文对不同的数据体进行了叠加处理, 获得相应的零偏移距剖面, 如图 10所示, 其中图 10a是基于有限差分法得到的叠加剖面, 图 10b至图 10d依次是0阶到二阶层间多次波叠加剖面。黑色箭头指的是与有限差分数据对应的一次反射波叠加结构, 橙色箭头表示与有限差分数据对应的一阶层间多次波叠加结构, 红色箭头表示对应的二阶层间多次波叠加结构。经过对比可发现不同阶次层间多次波的叠加剖面构造与有限差分计算的剖面构造相对应, 证明了单程波数值模拟层间多次波的正确性, 同时说明观测的地震记录是不同阶次的层间多次波和一次反射波的组合。

两个数值模拟结果表明不管在叠前还是叠后, 地震数据可以理解为不同阶次层间多次波(一次反射波被认为是0阶层间多次波)的叠加。利用本文层间多次波数值模拟方法, 联合测井资料, 可辅助实际地震数据层间多次波的识别, 指导其速度谱分析, 帮助拾取一次反射波能量团。利用本文方法还可以确定层间多次波的聚集区域, 便于指导后续的层间多次波的压制。

2.3 三维逆掩冲断层模型

利用本文方法对三维复杂模型进行层间多次波数值模拟。逆掩冲断层模型是公认的复杂模型, 其含有断层、背斜、逆掩冲断层等复杂构造, 横向非均匀性强。图 11是逆掩冲断层模型的反射率模型, 模型的大小为101×401×401个网格点。在模型的上表面(2km, 2km)位置设置一个20Hz的雷克子波震源。模型的网格尺度为10m, 采样间隔为2ms, 记录时长为2s。图 12是过震源点沿Crossline方向合成的单炮记录, 其中图 12a是一次反射波单炮记录, 图 12b和图 12c分别是含有一阶和二阶层间多次波的单炮记录。从模拟结果中可以看出本文方法有效分离了三维复杂构造模型中产生的层间多次波。对比有限差分方法和本文层间多次波模拟方法对地震记录的近偏移距道集与远偏移距道集的计算结果(图 13), 可以看出两种方法对近偏移距道集的地震记录计算结果较好, 但对远偏移距道集的计算结果有一定的差别。这是由于不同方法具有不同的属性导致。

3 讨论与结论

本研究提出了一种层间多次波模拟方法, 可在不同维度下对一次反射波和层间多次波进行有效分离, 并使用傅里叶有限差分方法作为波场延拓算子, 提高了在二维和三维情况下地震波的横向传播精度。波场分阶数值模拟分为散射和传播过程, 两个过程相分离使得波场分阶模拟具有很大的灵活性, 且本文方法具有迭代性, 循环次数越多, 波场信息越丰富。数值实例验证了本文方法可以有效地分离不同阶的层间多次波, 并且叠前和叠后合成数据表明地震资料是各阶层间多次波的叠加。数值模拟结果还表明, 本文层间多次波分阶模拟方法可以辅助识别实际地震资料中的层间多次波, 指导速度分析, 便于后续层间多次波的压制处理, 从而提高深部储层解释的准确度。