基于迭代反演的层间多次波压制方法

包培楠， 王孝， 谢俊法， 王维红， 石颖， 徐嘉亮， 杨育臣

1 东北石油大学地球科学学院， 大庆 163318 2 中国石油天然气股份有限公司勘探开发研究院西北分院， 兰州 730070 3 东北石油大学环渤海能源研究院， 秦皇岛 066004 4 东北石油大学非常规油气研究院， 大庆 163318 5 “陆相页岩油气成藏及高效开发”教育部重点实验室， 大庆 163318 6 密苏里科技大学， 罗拉 65401， 美国

0 引言

层间多次波是地震资料中的规则干扰波，其动校正时差和有效波差别较小，难以识别和有效压制，对有效波的成像和识别带来很大的困难(Berkhout, 2006)，同时层间多次波的存在和其他噪声一样，对后续地震资料的解释、反演和应用都会带来不利影响，导致构造和油气识别的精度下降，钻探成功率降低，对油气的勘探开发带来极大的风险(孙小东等，2020).因此，有效压制层间多次波，进而提高地震资料成像的精度，是油气地震勘探领域的一个重要研究方向.

多次波压制方法主要分为两类，一类是滤波方法，另一类是预测相减方法.滤波方法主要利用一次波与多次波的时差关系和周期特征进行多次波的识别和压制(Taner, 1980; Lokshtanov, 1999)，如抛物Radon滤波(Hampson, 1986)和双曲Radon滤波(Foster and Mosher, 1992; 石颖和王维红, 2012; Shi et al., 2017)等，该类方法计算效率高，算法容易实现，当有效波和多次波之间的动校正时差较大时，可取得满意的多次波压制效果.但对于复杂介质，如速度梯度较小(或速度反转)的情况，或构造变化剧烈的介质，应用滤波法难以有效识别有效波和多次波，往往得不到理想的多次波压制结果(Berkhout and Verschuur, 1997).预测相减方法基于波动理论，能更好的适应于复杂介质的情况，诸多地球物理工作者都对该类方法进行了系统深入的研究.该类多次波压制方法包括逆散射级数法(Araujo et al., 1994; Weglein et al., 1997, 2003; Otnes et al., 2004; Kristopher, 2017)、共聚焦点技术(Berkhout and Verschuur, 2006)以及近年来发展的基于虚同相轴的层间多次波压制方法(吴静等, 2013)和基于Marchenko的层间多次波压制方法(Behura et al., 2014; Verschuur and Berkhout, 2015；匡伟康等，2018).上述方法中，基于散射理论的逆散射级数法无需地下的先验信息，预测一次可以得到与所有界面相关的同阶层间多次波，在没有有效手段区分有效波和多次波时，是层间多次波压制较为有效的方法，但是该方法存在计算量大，预测远偏移距多次波效果差的缺点.共聚焦点方法适合复杂介质条件的多次波压制，但该方法需要全波场数据，当地震资料的近偏移距缺失时需要相应的波场重建方法进行数据重建(王锦妍等，2020)，且一次只能预测得到与某一界面相关的层间多次波，在一定程度上依赖初始速度模型以求取准确的聚焦算子.虚同相轴方法能够较准确地预测层间多次波，但对观测系统要求较高，尚未摆脱对人工操作的依赖.Marchenko层间多次波压制方法在模型数据上取得了很好的效果，也可应用于实际数据(Zhang and Slob, 2020)，但是对于层间多次波发育且信噪比较低的陆上地震资料多次波压制还存在计算不稳定的情况.

Jakubowicz(1998)提出了一种数据驱动构造层间多次波的新方法，刘战等 (2019) 对其进行了详细推导，巧妙地将地下散射点移到表面，但是对于实际数据无法进行精确的层间多次波压制.Ikelle(2006)、Ikelle等(2009)等通过引入虚源点的概念较为有效地解决了该方法应用于实际地震资料层间多次波压制的问题.Ikelle认为通过一步预测相减就可有效压制某一界面产生的所有层间多次波，但是该方法对多次波振幅预测精度不够.提高多次波压制效果的方法主要有两种方法，一是增大自适应相减的滤波算子长度以校正多次波的振幅和相位，另一种是提高多次波的预测精度.第一种方法对有效波的振幅产生较大伤害，针对提高多次波预测精度，在数据驱动和CFP层间多次波压制方法的基础上，提出一种基于迭代反演的层间多次波压制方法(MSI)，该方法直接利用地震数据本身进行层间多次波预测，将多次波压制看做优化问题，避免构建CFP道集所需的聚焦运算和地下速度的估计，该方法的优点有两个：一是计算效率高，二是多次波预测精度高.理论模型数据测试表明该方法具有高精度的特点.

1 方法原理

基于迭代反演的层间多次波压制方法，可以看作是自由表面多次波压制方法(SRME，Surface-Related Multiple Elimination)的扩展.在CFP层间多次波压制方法的基础上，首先借助于数据矩阵表示法和反馈迭代模型(WRW)，建立反射波的数学模型.然后利用“反馈迭代模型”表示出一次波与多次波之间存在的关系，预测出层间多次波.最终利用自适应匹配滤波法将预测的多次波从地震数据中减去.

1.1 CFP层间多次波产生及压制

假定在地表z0处激发震源产生地震波，如图1a所示，地震波在地下第zm层处发生反射，置于地表的检波器接收一次反射波.依据图1b所示的WRW模型，检波器接收到的一次反射波地震记录可表示为

ΔP(z0,z0)=D(z0)ΔX(z0,z0)S(z0)

=D(z0)ΔP-(z0,z0)，

(1)

其中:

ΔP-(z0,z0)=ΔX(z0,z0)S(z0)，

(2)

(3)

式中，ΔP(z0,z0)为只包含一次反射波的地震数据矩阵，D(z0)为检波点矩阵，ΔX(z0,z0)为一次反射波的传递矩阵，S(z0)为震源矩阵，ΔP-(z0,z0)为反射波场，W-(z0,zm)为上行传播矩阵，W+(zm,z0)为下行传播矩阵，R(zm,zm)为地下第zm层反射系数矩阵，m为地下反射层序号，(z0,z0)表示地表激发地表接收.

若地震波在地表z0处发生至少一次下行反射后，被置于地表的检波器所接收，其生成表面多次波的射线路径如图2a所示，依据图2b所示包含表面多次波的WRW模型，则地震记录可表示为

P(z0,z0)=D(z0)ΔX(z0,z0)[S(z0)+R(z0,z0)

×P-(z0,z0)]，

(4)

式中，P-(z0,z0)表示在不考虑检波点对地震波场影响情况下的含表面多次波的反射波场，R(z0,z0)表示向下延拓算子.则有:

P-(z0,z0)=ΔP-(z0,z0)+ΔX(z0,z0)[R(z0,z0)

×P-(z0,z0)]，

(5)

P(z0,z0)=D(z0)P-(z0,z0)，

(6)

由于总波场P(z0,z0)含一次波和表面多次波两部分，根据公式(4)，表面多次波{δM(z0,z0)}0可表示为

{δM(z0,z0)}0=D(z0)ΔX(z0,z0)[R(z0,z0)

×P-(z0,z0)].

(7)

引入界面算子A(z0,z0)，假定其与震源矩阵、检波点矩阵以及向下延拓算子有关，且存在关系式：

A(z0,z0)=S-1(z0)R(z0,z0)D-1(z0)，

(8)

假定在自由表面处，向下延拓算子R(z0,z0)可近似为-I(I为单位矩阵)，则有：

A(z0,z0)≈-[D(z0)S(z0)]-1，

(9)

若等式(7)右边乘以一个单位矩阵I=-A(z0,z0)[D(z0)S(z0)]，结合公式(1)和公式(5)，表面多次波又可表示为

{δM(z0,z0)}0=ΔP(z0,z0)A(z0,z0)P(z0,z0)，

(10)

从总波场中去除表面多次波，可得到压制多次波后的地震波场，即：

图1 一次波传播路径及其WRW模型(a) 一次波传播路径； (b) 产生一次波的WRW模型.Fig.1 Propagation path of primary waves and WRW model(a) Propagation path of primary waves; (b) WRW model of primary wave generation.

ΔP(z0,z0)=P(z0,z0)-{δM(z0,z0)}0.

(11)

若令X(z0,z0)为含一次波和多次波的传递矩阵，根据公式(1)，包含一次波和多次波的地震记录P(z0,z0)可表示为

P(z0,z0)=D(z0)X(z0,z0)S(z0).

(12)

假定地震波在地面以下的zn界面发生至少一次下行反射后，被置于地表的检波器所接收，产生层间多次波的射线路径如图3a所示，WRW模型如图3b所示，则只包含一次反射波和层间多次反射波的地震记录可表示为

{P(z0,z0)}0=D(z0){X(z0,z0)}0S(z0)，

(13)

式中，{}0表示与界面z0相关的多次波都已被压制，{P(z0,z0)}0表示与表面相关的多次波压制后的地震记录，{X(z0,z0)}0表示与表面相关的多次波压制后的传递矩阵.对于地表以下的界面，做类似定义：

{P(z0,z0)}n=D(z0){X(z0,z0)}nS(z0)，

(14)

式中，n=0,1,2,…,∞，当n=0时，表示与自由表面相关的情况.{}n表示关于界面z≤zn的多次波都已被压制，只有关于界面z>zn的层间多次波存在.

CFP方法是将表面多次波压制方法引入到层间多次波压制中，该方法假设与层z≤zn-1相关的所有多次波在之前都已被压制，得到向下外推炮记录{P(zn,z0)}n-1作为压制层间多次波程序的输入数据，表达式为

{P(zn,z0)}n-1=Γ(zn,z0){P(z0,z0)}n-1，

(15)

式中，{P(zn,z0)}n-1表示震源在地表z0处，检波器在地下zn处的地震记录.Γ(zn,z0)表示检波器一侧的向下延拓算子.根据公式(3)，关于界面zn的传递矩阵为

(16)

根据公式(7)，与界面zn相关的层间多次波为

{δM(z0,z0)}n=

(17)

图2 表面多次波传播路径及其WRW模型(a) 表面多次波传播路径； (b) 产生表面多次波的WRW模型.Fig.2 Propagation path of surface multiples and WRW model(a) Propagation path of surface multiples; (b) WRW model of generation of surface multiples.

图3 层间多次波传播路径及其WRW模型(a) 边界zn处产生层间多次波的波传播路径； (b) 一次波(m>n)和层间多次波正演的WRW反馈模型.Fig.3 Propagation path of internal multiples and WRW model(a) Propagation path of internal multiples at boundary zn; (b) WRW feedback model of forward modeling for primary waves (m>n) and internal multiples.

{δM(z0,z0)}n=

(18)

A(zn,zn)=S-1(zn)R(zn,zn)D-1(zn)

≈-[D(zn)S(zn)]-1，

(19)

则压制与该界面相关的层间多次波公式为

{P(z0,z0)}n={P(z0,z0)}n-1-{δM(z0,z0)}n

×{P(zn,z0)}n-1.

(20)

通过以上分析可知，层间多次波压制算法与表面多次波压制算法相似，只是对于层间多次波而言，需要由检波器在深度zn时的炮记录来代替检波器位于表面时的炮记录.

1.2 MSI层间多次波压制算法

CFP层间多次波压制方法需要求取地震数据的延拓算子，在很大程度上增加了层间多次波预测过程的计算量.为了降低计算成本，由前面阐述的CFP方法入手，假设经过k次迭代压制层间多次波后的地震记录{P(z0,z0)}′k可以由上一次迭代结果{P(z0,z0)}′k-1得到，引入一个类似于界面算子A(zn,zn)的卷积因子T(k-1)，则经过k次迭代的层间多次波{M(z0,z0)}′k表示为

{M(z0,z0)}′k=T(k-1)[P(z0,z0)-{P(z0,z0)}′k-1]，

(21)

由公式(20)可知，界面算子A(zn,zn)为

(22)

将公式(22)代入公式(20)，得到：

×({P(z0,z0)}n-1-{P(z0,z0)}n).

(23)

在实际应用中，矩阵的逆可以根据最小平方的形式进行求取(Berkhout，2006)，则公式(23)又可表示为

×({P(z0,z0)n-1-{P(z0,z0)}n),

(24)

公式(24)中{δM(z0,z0)}n表示关于界面zn的层间多次波，若考虑所有的层间多次波，公式(24)可写为

{M(z0,z0)}′k={P(z0,zn)}′k({P(z0,z0)}′k-1)H

×[{P(z0,z0)}′k-1({P(z0,z0)}′k-1)H]-1

×(P(z0,z0)-{P(z0,z0)}′k),

(25)

式中，H表示共轭转置，则卷积T(k-1)因子为

T(k-1)={P(z0,z0)}′k-1({P(z0,z0)}′k-2)H

×[{P(z0,z0)}′k-2({P(z0,z0)}′k-2)H]-1，

(26)

最终压制的多次波结果为

{P(z0,z0)}′k=P(z0,z0)-{M(z0,z0)}′k.

(27)

由以上可知，在MSI算法中，不需要显式的表面算子和显式的震源矩阵，表面算子被原始数据P(z0,z0)及先前两步的多次波压制结果{P(z0,z0)}′k-1和{P(z0,z0)}′k-2所替代.MSI算法在多次波预测中隐含的考虑了表面算子的空间变化，并在多次波预测不断迭代反演更新中实现层间算子的空间变化，因此在迭代反演过程中，未被预测或者预测不完整的层间多次波将随迭代次数的增加逐渐完善，同时在一定程度上解决了后续多次波自适应相减方法中的非线性问题.

在实际应用中，将去除直达波的地震数据输入到程序中进行迭代以求取卷积因子并对多次波模型进行更新，根据数据的具体情况及计算效率的综合考虑选择合适的迭代次数(一般迭代2～3次)得到最终的多次波模型，然后利用自适应匹配滤波进行相减.其实现流程如图4所示.

图4 MSI层间多次波压制算法实现流程图Fig.4 Implementation flowchart of MSI internal multiple suppression algorithm

2 模型试算

为了证明上文所述方法的实用性和有效性，采用两个不同的地质模型进行测试.模型一为含楔状构造的地质模型，模型二是针对我国西北某地区的地质特征，通过测井资料插值的方式生成速度模型正演的数模数据.

2.1 模型一

首先对图5所示模型进行测试，其中图5a为速度模型，包含三个强阻抗界面，第三层为高速层，地震波遇到层界面将产生较强的层间多次波.图5b为密度模型，第三层为高密度层.模型大小为4500 m×1400 m，采用的观测系统为中间放炮两边接收，炮点和检波点均匀分布于地表，炮间距和道间距均为20 m，震源激发201炮，每炮201个检波器接收.激发震源采用主频为25 Hz的雷克子波，采样间隔为4 ms，地震反射记录为2.5 s.为了避免表面多次波的影响，地震波场正演计算时模型四周都采用了吸收边界.

图6为第101炮地震波场正演模拟记录，黑色箭头指示的同相轴分别是三个反射界面的一次波，其余均为层间多次波.图7为预测的层间多次波，相比迭代1次，经过迭代3次预测的层间多次波信息更完整，其相位和能量都更与实际的层间多次波相符，尤其在1.3～1.5 s之间，经过3次迭代后，这个时间段的层间多次波基本上都被预测出.注意在层间多次波预测过程中，为了避免产生噪声而伤害有效波，会设置时窗来控制多次波预测范围.图8为经过3次迭代压制层间多次波后地震记录，其中图8a为迭代1次压制层间多次波后的地震记录，图8b为迭代3次压制层间多次波后的地震记录，观察可知，迭代1次后，层间多次波得到部分压制，一次反射波的原始波场特征完好；迭代3次即获得较为理想的层间多次波压制效果，多次波残余能量很少，同时一次波的振幅和相位信息得到较好的保护.图9为3次迭代层间多次波压制后的零偏移距剖面，与原始模型地质特征相符，层位特征明显，无虚假构造出现.

图5 (a) 速度模型； (b) 密度模型Fig.5 (a) Velocity model； (b) Density model

图6 含层间多次波的地震记录Fig.6 Seismic record with internal multiples

图10为单道数据进行层间多次波压制前后的对比图，黑色虚曲线为原始单道数据，绿色实曲线为压制层间多次波后的，红色实曲线为预测的层间多次波.从图中可以看到预测的层间多次波与实际的层间多次波在到时、相位上都有很好的一致性，只在振幅上有微弱差异，此差异可利用后续的自适应匹配相减来消除.黑色与红色实曲线的对比也可以看出自适应匹配相减后层间多次波得到了有效的压制，且不伤害有效波，说明迭代反演的层间多次波压制方法对模型数据的压制效果较好.

图7 预测的层间多次波(a) 迭代1次； (b) 迭代3次.Fig.7 Predicted internal multiples(a) By one iteration; (b) By three iterations.

2.2 模型二

图11为通过层位数据和测井资料控制的方式得到西北地区研究区的速度模型，速度模型大小为9300 m×3200 m，共有17个波阻抗分界面，图中白色箭头所指向的四个层为低速薄煤层，煤层速度在2800 m·s-1左右，厚度约10 m，与上下高速地层(4000 m·s-1左右)形成强反射界面，层间多次波较发育，因此层间多次波与一次波的差异很小且能量较强，压制时很容易伤害有效波，压制难度较大.

依据速度模型，采用有限差分正演方法模拟地震记录，图12为炮数据层间多次波压制前后对比图，每炮400道，道间距为10 m，炮检距为10 m，采样间隔为0.2 ms，模型的网格间距在水平方向和垂直方向上均为2 m，由于MSI方法是一个迭代反演的过程，直达波的存在会影响层间多次波的压制效果，在压制层间多次波之前首先需去除直达波.图12a、b为压制层间多次波前后的单炮记录对比，观察图12a发现，层间多次波主要分布在2.7～4.0 s之间.从图12b易知，在反射时间3.0～4.0 s之间层间多次波得到明显压制，而有效波得到较好地保护.对层间多次波压制前后的数据进行频谱分析，如图12c、d所示，主频范围稍有拓宽，拓宽5 Hz左右，能量整体有所抬升，这在一定程度上提高了分辨率，说明该方法未伤及有效波，只有层间多次波被有效压制(图12中白色箭头所示都为层间多次波).

图13为图12a、b实际数模单炮记录的局部放大图，对比发现，2.7～4.0 s之间的层间多次波得到有效压制.

图14为层间多次波压制前后叠加剖面对比图，从图中可以看出层间多次波压制效果较为明显，为了能够更清楚的看出层间多次波压制效果，将整个地震剖面(图14a、b)上白色框所示部分放大，分别得到图14c、d.从整个叠加剖面及局部放大图中可以进一步看出，3.0～4.0 s之间的层间多次波基本被压制，压制效果较好.

3 结论

在共聚焦点层间多次波压制方法基础上，引入迭代反演策略，将层间多次波压制的计算转化为反演过程，用卷积因子代替CFP算法中的聚焦算子，形成基于迭代反演的层间多次波压制方法，进而应用该方法实现数据驱动的高精度层间多次波预测.MSI方法通过将地下散射点移至地表，使得该方法的实用性大大增强，而且是完全数据驱动，不需要地下介质的任何先验信息，层间多次波预测精度较高，MSI方法的试算和应用表明，一般通过三次迭代即可得到高精度的多次波预测结果，具有很高的计算效率.采用自适应匹配滤波方法将预测得到的层间多次波进行压制，可得到多次波压制结果.理论模拟算例表明MSI方法数据适应性强、计算精度高、计算成本低，具有良好的实际地震资料层间多次波压制的应用前景.

图8 层间多次波压制(a) 迭代1次； (b) 迭代3次.Fig.8 Internal multiples suppression(a) By one iteration; (b) By three iterations.

图9 经3次迭代压制层间多次波的零偏移距剖面(a) 原始零偏移距剖面； (b) 迭代3次.Fig.9 Zero-offset profiles after internal multiples suppression by three iterations(a) Original zero-offset profile; (b) After three iterations.

图10 单道数据层间多次波压制前后对比图Fig.10 Comparison of single-trace data before and after internal multiples suppression

图11 速度模型Fig.11 Velocity model

图12 层间多次波压制前后对比图(a) 层间多次波压制前的单炮记录； (b) 层间多次波压制后的单炮记录； (c) 层间多次波压制前频谱； (d) 层间多次波压制后频谱.Fig.12 Comparison before and after internal multiples suppression(a) Shot gather before internal multiples suppression; (b) Shot gather after the internal multiples suppression; (c) Spectrum before internal multiple suppression; (d) Spectrum after internal multiple suppression.

图13 单炮部分数据层间多次波压制前后对比图(a) 层间多次波压制前； (b) 层间多次波压制后.Fig.13 Comparison of part shot gather before and after internal multiples suppression(a) Before multiples suppression; (b) After internal multiple suppression.

图14 叠加数据层间多次波压制前后对比图(a) 层间多次波压制前； (b) 层间多次波压制后； (c) 部分层间多次波压制前； (d) 部分层间多次波压制后.Fig.14 Comparison of stacked data before and after multiples suppression(a) Stacked data before internal multiples suppression; (b) Stacked data after multiples suppression; (c) Part stacked data before multiples suppression; (d) Part stacked data after internal multiple suppression.