层间多次波傅里叶有限差分偏移成像

黄建东， 胡天跃， 王尚旭

1 北京大学地球与空间科学学院， 北京 100871 2 中国石油大学(北京)油气资源与探测国家重点实验室， 北京 102249 3 清华大学数学科学系， 北京 100084

0 引言

要实现深层复杂构造的高精度成像，多次波是一个不可忽视的问题.在勘探地震学领域，多次波指的是在地层中发生过一次以上向下反射且被地表检波器接收的地震波.按照反射层位不同，多次波可分为表面多次波和层间多次波.在传统成像方法中，多次波通常被认为是一种噪声，会使反射波发生畸变.在深层地区，层间多次波能量高于反射波能量导致拾取不准确的速度谱，严重影响地震偏移成像和解释精度.因此，地震数据中的多次波必须被压制，这个过程也称为去噪.多种方法被用于压制层间多次波，有逆散射级数法(Weglein et al., 2003；Araújo et al., 1994；Yang and Weglein, 2017)，拉东变换法(Sacchi and Ulrych, 1995；Sacchi and Porsani, 1999；Abbad et al., 2011)，预测反褶积方法(Peacock and Treitel, 1969；Taner, 1980)，虚同相轴方法(Ikelle, 2006；吴静等，2013；Liu et al., 2018)等.

研究发现多次波携带的构造信息比反射波更加丰富，包括小角度信息，且有更高的横向照明孔径和覆盖范围(如图1所示).因此，许多学者致力于研究多次波成像.多次波成像研究大致可分为以下几类：一种是将表面多次波转化为反射波，再利用传统的逆时偏移方法进行成像(Berkhout and Verschuur, 2003; Verschuur and Berkhout, 2005).或者是将炮集数据作为正传波场，基于预测方法获得表面多次波数据作为反传波场，再利用互相关方法进行偏移成像(Liu et al., 2011；Liu et al., 2015).还有一些学者基于反演思想，利用最小二乘方法对表面多次波进行成像，提高成像质量(Liu et al., 2016, 2020).叶月明等(2019)提出只利用一阶和二阶表面多次波进行成像，进一步提高多次波的成像效果.Wapenaar等(2014)详细推导了Marchenko方程，创新性地提出利用Marchenko方法对层间多次波进行聚焦成像，但是目前只是停留在理论阶段.Berkhout (2014a,b)提出利用伴随方法对层间多次波进行成像.Li和Wang(2019)提出基于稀疏约束反演的表面多次波成像方法.刘伊克等(2018)详细推导并分析了多次波成像的基本原理，指出多次波成像的巨大挑战在于不同阶次多次波对应的地震波场之间互相关将产生大量的串扰噪声.目前有效降低串扰噪声方法有最小二乘逆时偏移成像(刘学建和刘伊克，2016)，或是将一次波和多次波进行联合反演成像(Lu et al., 2011).上述可知，多次波成像研究在表面多次波方面比较完善，在于表面多次波传播机理比较简单，容易进行预测和识别，与一次波具有一定的周期性关系.而层间多次波形成机理复杂，研究难度大，对层间多次波研究比较少.但是层间多次波的存在严重影响深层地震勘探研究，尤其对于陆地勘探，层间多次波极为发育，严重影响地震数据的处理和解释.在中国西北和西南地区，碳酸盐岩储层的层间多次波发育，严重降低该地区储层成像和解释的准确度.利用层间多次波成像可以增加复杂构造区域的有用信息，辅助勘探区块解释.对于油气勘探来说，层间多次波成像是非常重要的.但是，由于其成像的复杂性，科学界对利用层间多次波进行成像的研究也相对较少.

图1 地震波传播示意图Fig.1 Diagram of seismic wave propagation

本文通过修改传统傅里叶有限差分偏移方法，使得该方法适用于层间多次波成像.将层间多次波作为反传数据，利用修改的相关成像条件对层间多次波进行偏移成像.针对不同阶层间多次波之间产生的串扰问题，基于波场延拓和全波场方法获得对应阶的上下行波场，对分离波场分别成像.在数值部分，利用不同的模型验证层间多次波傅里叶有限差分偏移方法的有效性和适用性.

1 层间多次波成像原理

1.1 方法和理论

传统地震波成像方法是在每一个深度点利用一次反射波的正传波场和反传波场在t=0时刻进行互相关成像(如图2a所示).但是接收到的实际地震资料不仅包含一次反射波，还有不同阶次的多次波.多次波成像条件与一次反射波成像条件不同，利用传统互相关方法对含有多次波的地震数据进行偏移成像，会形成许多虚假的同相轴，影响地震剖面的构造解释.表面多次波逆时偏移成像基本原理是利用一次反射波作为震源激发正传波场与表面多次波作为震源激发反传波场进行互相关成像(如图2b所示).

图2 互相关成像示意图(a) 一次反射波成像； (b) 表面多次波成像.Fig.2 Diagram of cross-correlation imaging(a) Migration imaging for primaries; (b) Migration imaging for surface multiples.

本文提出利用傅里叶有限差分偏移方法进行层间多次波成像，基本原理如图3所示，在每一个深度点对正传多次波波场和反传波场进行互相关成像.对于M1成像路径来说，正向传播波场传播至x1点与M1检波点的反向传播波场相遇成像；对于M2成像路径，正向传播波场传播至x2点与M2检波点的反向传播波场相遇成像.对于其他类型层间多次波成像路径，可以此类推.图中的红色圆点表示层间多次波成像的位置.

根据Liu等(2011)提出的多次波成像理论，本文进一步发展完善，提出一种只针对层间多次波的成像公式，其表达式为

(1)

其中，Im表示成像结果，MF表示正向传播波场，MB表示反向传播波场，tmax表示最大的接收时间.符号M表示所有阶次的层间多次波，进一步将其展开可得：

图3 层间多次波成像示意图Fig.3 Diagram ofmigration imaging for internal multiples

M(x,z;t)=M1(x,z;t)+M2(x,z;t)+M3(x,z;t)+…+Mn(x,z;t)，

(2)

其中，M1表示一阶层间多次波，M2表示二阶层间多次波，Mn表示n阶层间多次波.把公式(2)代入公式(1)中展开可得：

+…，

(3)

式中，第一项表示不同阶次层间多次波的真实成像结果，第二项和第三项表示假成像结果.根据公式(3)可知，假成像结果的组合数明显多于真成像结果.因此利用层间多次波进行成像，由于不同阶数多次波之间的成像条件不同，对多次波使用单一成像条件时，会导致不同阶数多次波之间形成串扰噪声，降低成像质量.刘伊克等(2018)、刘学建和刘伊克(2015)文章中也指出过表面多次波成像的串扰噪声问题.但是对于层间多次波来说，其成像条件更加复杂，不仅不同阶数的多次波之间会形成串扰，即使是同阶层间多次波之间也会形成串扰.这是层间多次波区别于表面多次波的特征之一，也是目前无法对实际地震资料中层间多次波进行有效成像的根本原因.但是进一步的试验发现，利用层间多次波成像形成的真实构造同相轴能量比虚假构造的同相轴能量更强，且具有一定周期性，可利用一些叠后去噪方法消除虚假构造，如叠后虚同相轴方法，滤波方法等.本文不考虑去噪问题，对所有成像结果没有进行压制处理.

1.2 傅里叶有限差分偏移

横向非均匀介质中的二维频率域单程波方程表达式为

(4)

(5)

(6)

和

(7)

其中，方程(5)是相移算子(Gazdag, 1978)；方程(6)是慢度修正算子(Stoffa et al.,1990),Δl=1/v-1/v0表示慢度；方程(7)是横向上强速度差异的修正算子(Huang and Fehler, 2000)，a和b为优化系数.

基于傅里叶有限差分延拓算子和全波场(FWM)方法(本文不再赘述，具体原理可参考Berkhout，2014a，b；匡伟康等，2020；Huang et al., 2021)，可以得到不同阶次层间多次波对应的上行波场和下行波场.再由上一部分提供的层间多次波成像理论依据，推导得到多炮频率域层间多次波互相关成像条件为

R(zm,xn)=

(8)

其中，R表示反射系数，即成像结果，N表示炮数，σ表示稳定常数，[·]*表示[·]的共轭.基于层间多次波的互相关成像条件，可在每个深度点上对层间多次波进行偏移成像.而上下行波场分别成像可在一定程度上压制层间多次波成像中的串扰噪声.

2 数值实验

本节中，若干个模型被用于验证本文所提出层间多次波成像方法的有效性和适用性.第一个模型用于测试方法的正确性；第二个模型用于测试方法对横向不均匀介质的适应性.最后，基于SEG/EAGE盐丘模型和Sigsbee模型产生的层间多次波进行偏移成像，验证了利用层间多次波偏移成像对一些特殊构造如断层，孔洞，裂缝等拥有更高的成像精度.

2.1 层状介质模型

构建一个三层层状模型，如图4所示.模型的大小为2000 m×2000 m，网格大小为10 m，纵向和横向上分别有201个网格点.炮间距为100 m，总共21炮.模型速度从上到下依次为2300，1500，3000 m·s-1，产生两个强反射界面.介质的密度为2000 kg·m-3，震源为30 Hz的雷克子波，时间采样间隔为 1 ms.利用Berkhout(2014a,b)提出的全波场方法，可以得到全波场地震记录，包括反射波，各阶表面多次波和各阶层间多次波.

图4 层状介质模型Fig.4 Layered media model

本文主要是研究层间多次波，假设模型表面为吸收边界条件.图5是被地表检波器记录的反射波，一阶和二阶层间多次波炮集数据.基于本文修改的互相关成像条件对层间多次波进行偏移成像.对于该层状模型来说，不同阶次的层间多次波成像路径如图6所示，红色圆点表示成像点.利用层间多次波的正传波场和反传波场在每个层位点进行互相关成像.图7是反射波和不同阶次层间多次波偏移成像结果对比，其中(a)是一次反射波偏移结果，(b)和(c)分别是一阶和二阶层间多次波偏移结果.从偏移结果可知基于图6成像路径，应用层间多次波层上成像，由于层间多次波在第一个反射界面上没有成像点，无法在该层位进行成像(图7b和c).和原始速度模型对比可以看出第二个反射界面构造通过层间多次波偏移成像准确恢复.该实验结果验证了本文修改的层间多次波成像方法的准确性和可用性.

图5 合成地震记录(a) 全波场记录； (b) 一次波； (c) 一阶层间多次波； (d) 二阶层间多次波.Fig.5 Synthetic seismic record(a) Full wavefield record; (b) Primaries; (b) First-order internal multiples; (b) Second-order internal multiples.

2.2 倾斜界面模型

简单的层状介质模型验证了本文方法的正确性.为了进一步验证其对复杂模型的适用性，一个倾斜界面模型被构建，速度结构如图8所示，该模型含有一个较陡的倾斜界面，两个平层界面.模型大小为2000 m×2000 m，速度从上到下依次为3000、4000、4500和3500 m·s-1.利用上面所提出的成像条件对该模型的层间多次波进行偏移成像.图9是层间多次波偏移成像结果，从图中可以观察到除了真实的地层层位，还有一些虚假的层位结构产生，并且第一个反射界面缺失.橙色箭头所指的是真实层位的同相轴，而黑色箭头所指的是虚假层位的同相轴.通过前面理论部分对层间多次波成像的原理分析，可以知道虚假构造是由于同阶层间多次波之间相互串扰导致的.但是进一步观察发现，真实层位对应的同相轴比虚假层位的同相轴能量更强，且它们之间存在一定的周期性关系.和原始的速度模型对比，可知恢复的真实层位同相轴很好地反映了模型的地层结构.

图6 层间多次波偏移成像示意图(a) 一阶层间多次波相关成像； (b) 二阶层间多次波相关成像.Fig.6 Diagram ofcross-correlation imaging for internal multiples(a) First-order internal multiples cross-correlation condition; (b) Second-order internal multiples cross-correlation condition.

图7 平层模型地震波偏移成像结果(a) 一次波偏移结果； (b) 一阶层间多次波偏移结果； (c) 二阶层间多次波偏移结果.Fig.7 Seismic wave migration of horizontal layer model(a) Migration imaging for primaries; (b) Migration imaging for first-order internal multiples; (c) Migration imaging for second-order internal multiples.

图8 倾斜界面模型Fig.8 Tilt interface model

图9 层间多次波偏移成像结果 橙色箭头表示真实层位构造(能量强)，黑色箭头 表示虚假层位构造(能量弱)Fig.9 Migration imaging for internal multiples, where the orange arrow indicates the real structures (strong energy), and the black arrow indicates the false structures (weak energy)

2.3 SEG/EAGE盐丘模型

由于地下盐丘高速体的存在，对地震波传播有一定屏蔽作用.常规的反射波成像无法恢复盐下构造，导致盐下构造成像不清楚，而层间多次波成像对盐下构造恢复拥有一定的优越性.本文利用SEG/EAGE盐丘模型进行层间多次波成像.图10是盐丘速度模型，该模型含有多种复杂的构造—高速盐丘体、复杂的断层和裂缝.使用主频30 Hz的雷克子波作为震源，最长的记录时间为6 s，1ms采样间隔.该模型横向上有1341个网格点，纵向上有301个网格点.

图10 SEG/EAGE声波速度模型Fig.10 SEG/EAGE acoustic velocity model

图11是一阶层间多次波下行波场偏移成像结果.从图中可观察到盐丘模型主要层位都得到了准确的恢复，但是由于不同层位产生的同阶层间多次波存在成像串扰，所以会出现一些虚假同相轴.受高速盐丘体的影响(速度差异大)，在盐丘体附近的层间多次波能量非常强，多次波偏移成像结果更加清楚.而远离盐丘体地区，由于地层之间波阻抗差异小，导致多次波不发育，多次波成像效果不明显.图12是反射波偏移成像结果，可以观察到浅层以及盐丘体上面的构造恢复非常好.但是盐丘体下面构造的同相轴能量非常弱，甚至缺失，成像非常不清楚，这是由于高速盐丘体的屏蔽作用造成的.图13是一阶层间多次波上行波场偏移结果，可以看到盐丘体附近的构造恢复非常好，尤其是盐下地层和裂缝构造非常清楚.这是由于盐丘体附近波阻抗差异大，产生的多次波能量强，成像效果也更好.但是不可避免地产生了串扰噪声.通过对比图12和13，可以观察到层间多次波对于盐丘体以及断层构造的成像具有良好的优越性.在反射波偏移结果中，由于高速盐丘体的存在，掩盖了盐下构造，层位恢复不明显，甚至缺失.而对于层间多次波，盐下构造成像非常明显，断层和裂缝同相轴信息得到准确恢复.红色箭头表示地层的同相轴对比，黄色框内的同相轴表示裂缝构造.结果说明层间多次波成像可以弥补反射波成像的不足，增加偏移成像的构造信息，对于复杂的构造，也比反射波拥有更高的成像优越性.多次波成像也为地下成像提供了额外的照明度.图14是反射波偏移和层间多次波偏移的照明度对比，在右边使用相同炮数的叠加结果.可以看出层间多次波在横向上拥有更宽的照明度.

图11 层间多次波下行波场偏移成像Fig.11 Migration imaging for first-order internal multiples based on down-going wavefield

图13 层间多次波上行波场偏移成像Fig.13 Migration imaging for internal multiples based on up-going wavefield

图14 地震波照明度对比(a) 反射波； (b) 层间多次波.Fig.14 Comparison of seismic wave illumination(a) Primaries; (b) Internal multiples.

2.4 Sigsbee 2B模型

Sigsbee 2B模型是非常复杂的储层模型，包含有不同尺度和不同组合的孔洞，高速盐丘体，断层等构造，速度结构如图15所示.利用本文的层间多次波傅里叶有限差分偏移方法对其进行成像.层间多次波成像结果如图16所示，由于地层结构非常复杂，利用层间多次波进行成像，串扰噪声非常严重.噪声严重掩盖了真实的地层层位，但是高速的盐丘体形态大小非常清楚，以及盐丘体附近的能量非常强，孔洞组合形态也很清晰.

图15 Sigsbee 2B 速度模型Fig.15 Sigsbee 2B velocity model

图16 层间多次波偏移结果Fig.16 Migration imaging for internal multiples

图17 局部放大图(a) 区域1； (b) 区域2.Fig.17 Zoomed view(a) Area 1 ofFigure 16; (b) Area 2 ofFigure 16.

图17是图16局部区域1和2的放大图，和原始的速度模型对比，发现利用层间多次波成像方法得到的构造与真实层位构造都可以准确对上.即使多次波的串扰噪声严重，地层的结构也是非常清楚.不同地层之间会产生成像串扰噪声，但是真实层位恢复的能量比虚假层位能量高.数值实验的结果也验证了本文方法的可靠性.

3 讨论和结论

层间多次波傅里叶有限差分偏移是单程波偏移方法进一步的发展，通过循环方式可获得不同阶次的层间多次波场，同时将层间多次波进行反向传播，与正传多次波场相关成像.和反射波偏移方法相比，本文通过修改传统成像条件获得适用于层间多次波偏移的成像条件.

层间多次波傅里叶有限差分偏移方法利用层间多次波进行偏移成像，可以准确恢复地层的结构，同时也会产生许多的串扰噪声.数值实例证明层间多次波成像方法对盐下构造恢复、裂缝和断层成像效果更加明显.而一次反射波对简单构造成像精度高，对盐下构造成像精度不足，利用层间多次波成像正好可以弥补这一缺点.层间多次波在波阻抗差异大的地区更加发育，成像也更加清楚.层间多次波偏移成像可以帮助恢复盐丘体下方的地层结构，为下地表成像提供额外的照明.同时层间多次波成像也为地震波成像提供了一个可用的思路，将地震波中的一次反射波和层间多次波进行分离，分别进行成像.但是层间多次波的形成机理非常复杂，不同地层的层间多次波波场会相互之间形成串扰，导致其成像困难，这也是我们下一步解决的问题之一.同时该方法也可进一步发展到三维层间多次波成像.