基于反馈迭代模型的多次波压制方法综述

贺紫林,李振春,李志娜,徐夷鹏

(中国石油大学(华东) 地球科学与技术学院，山东 青岛 266580)

0 引言

多次波压制问题一直是石油勘探领域中非常重要的一个课题，也是地震数据处理过程中的重点和难点之一。随着勘探工作的深入和研究目标的复杂化，地下构造的成像精度要求也越来越高。多次波的存在会降低地震资料的信噪比，在地震成像过程中产生假象，影响成像结果的可靠性和真实性。而成像结果的质量直接关系着地震资料解释结果的准确性，进一步影响着井位部署和油气勘探的成功率[1]。

随着油气勘探的进步、计算能力的不断提升，多次波压制方法也在不断发展。目前，多次波的压制方法大致可分为两类。一类是基于一次波和多次波之间周期性差异和动校正时差差异的滤波方法，简称为滤波法，如f-k滤波、Radon变换、预测反褶积等。此类方法简单，但是在地下复杂介质中适应性差，对多次波压制效果不理想[2]。一类是基于波动方程的预测减去方法，这种方法通过波动方程模拟实际波场或反演地震数据来预测多次波，然后将其从原始地震数据中减去，简称预测减去法，其中包括反馈迭代法、恒定内插法、波场外推法、逆散射级数法等。这类方法克服了滤波法的缺点，对复杂地下构造具有很好的适应性。根据产生多次波最浅的下行反射边界，多次波被分为两类：自由表面多次波(surface-related multiple)和层间多次波(internal multiple)，随之出现了自由表面多次波衰减和层间多次波衰减的概念。自由表面多次波发育广泛，在海洋地震勘探资料中尤为明显，它的存在影响着地震资料一次波识别的准确性，从而使处理的地震剖面存在假象，得到错误的解释结果。因此，在地震数据处理过程中自由表面多次波压制工作是必不可少的。鉴于波动方程理论的高精度特性，目前，采用预测减去方法中基于反馈迭代模型理论的一类方法实现对自由表面多次波的压制是一个很好的选择。

1 反馈迭代模型

20世纪60年代，Anstey和Newman提出有效波可以通过地震道信号的褶积转化为多次波，这为之后的反馈迭代法奠定了理论基础。Riley等[3]于1976年提出了二维垂直入射时的自由表面多次波模拟算法，但未找到合适的反演方案。Kennett[4]于1979年实现了一维自由表面多次波的模拟和反演，但由于该算法对地下介质和资料采集方式过多的简化和假设，在实际应用中效果并不是很成功。

20世纪80年代，Berkhout[5]提出了自由表面多次波的反馈理论框架，将震源和滤波器的特性考虑其中，突破了对介质和波入射方向的要求，并对自由表面多次波达到了较好的反演效果，奠定了基于反馈迭代模型的多次波压制方法的数学物理基础。

反馈迭代模型是在单频分量矩阵之间的计算上建立的，假设将观测到的地震数据记为p(xr,xs,t)，将地震数据转换到频率域变为p(xr,xs,ω)，提取单频分量的地震数据矩阵记为P(如图1)[6]。一次波(即无自由表层多次波的期望数据)单频分量矩阵记为P0，自由表面多次波单频分量矩阵记为M。

图1 提取单频分量过程Fig.1 Extraction process of single frequency component matrix

将地震波在地层中的传播看作是大地滤波过程，若无自由表面，一次波可表示为：P0=X0S(如图2a)[7-8]。而自由表面的存在使得这个滤波器变成一个带反馈的滤波器(如图2b)[7-8]，其反馈模型表达式总结为：

P=X0(S+RP)=X0S+X0RP=P0+M,

(1)

a—不含自由表面多次波;b—含自由表面多次波a—forward model of seismic data without surface-related multiple;b—forward model of seismic data with surface-related multiple图2 反馈模型示意Fig.2 Feedback model

式(1)反馈模型的提出，在没有过多简化的情况下提供了一种新的自由表面多次波压制思路，人们在此理论上，首先发展出了SRME(surface-related multiple elimination)自由表面多次波压制方法，并在其基础上，扬长避短，在发现问题、解决问题的过程中不断提出更加优化的自由表面多次波压制方法。下面将对基于反馈迭代模型的四种常用方法进行简要的介绍。

2 基于反馈迭代模型的多次波压制方法

2.1 SRME方法

SRME是最经典的反馈迭代方法，Verschuur[9]于1991年在反馈迭代模型理论上通过引入一次波能量最小假设实现了多次波的自适应相减，很大程度推动了SRME的发展。Berkhout和Verschuur[7-8]在1997年提出了迭代SRME方法，将自由表面算子的估计从一个非线性问题转化为线性问题，并通过实际应用说明收敛速度比较快，几次迭代就可以较好地实现自由表面多次波压制。

该方法是由式(1)入手，推得M=P0S-1RP=P0AP，于是可概括为：

(2)

式中：m(xr,xs,t)为时间域的自由表面多次波；a(t)为时间域的自由表面算子。在实际生产中，为使预测的多次波更接近实际多次波，可以将单炮或者单接收点道集应用最小二乘滤波对a更准确地求取，并且预测因子a的长度在一定范围内，多次波的预测会随着因子长度的增大而更为准确，但当超过此范围时不仅对预测结果影响不大，反而会增加计算量。因此，为了在最少的计算量下得到令人满意的自由表面多次波压制效果，多次波预测过程中需要选取合适的预测因子长度。一次波能量最小假设的引入虽然实现了多次波的自适应相减，但在一次波和多次波能量存在干涉时，此假设并不成立，从而会造成一次波能量的损伤。

SRME方法包括两部分，一部分是自由表面多次波的预测，另一部分是把预测出的自由表面多次波有效地从原始数据中减去。虽然预测的多次波在传播路径上做到了与实际的匹配，但在波的相位、时间、振幅等方面缺少充分的考虑，预测得到的多次波与实际的多次波相比具有一定差异，从地震资料中直接减去预测的多次波可能会得到错误的一次波信息，所以需要对预测得到的多次波与原始数据先进行匹配再进行相减。

在长期的研究进程中，国内外学者们提出了不同的匹配相减方法，Verschuur等[6]于1992年提出频率域最小能量滤波方法来消除自由表面多次波，并在1997年[7-8]提出了时间域单道维纳滤波方法。为达到更好的匹配相减结果，Monk[10]在1993年提出了伪多道维纳滤波。1999年Spitz等[11]提出了模式匹配滤波法。2003年，Guitton[12]提出了预测误差滤波，同年，Wang[13]对常规匹配方法进行了总结，并在此基础上提出了扩展伪多道维纳滤波，将匹配滤波器的约束由二维扩展到了三维。2004年陆文凯等[14]提出了基于独立变量分析方法，使用高阶能量统计函数来进行多次波的匹配相减，得到较好的效果。

之后许多学者进一步发展了SRME，Dedem和Verschuur[15-18]将SRME发展到了三维数据处理中，并在2005年提出将稀疏反演与三维SRME相结合。Anatoly等[19]于2006年通过共偏移距道集的DMO和反DMO操作进行数据重建，实现了三维SRME多次波压制。石颖等[20]在2013年提出全三维多次波预测方法，通过GPU并行的方式加速，在精度和速度上有了明显提升。井洪亮等[21]于2015年在GPU并行下实现基于L1、L2范数联合的SRME自适应相减，克服了L2范数的条件限制，也降低了L1范数的计算成本。

图3展示了SRME方法对三层水平层状模型正演数据的自由表面多次波压制结果[22]。由图3c可见，SRME在一定程度上压制了自由表面多次波能量，但仍存在明显的能量残余(如黑色箭头所示)，并在一次波与多次波能量干涉处，存在一次波能量损失的情况(如橙色箭头所示)。

a—原始单炮数据;b—真实一次波;c—多次波压制结果a—original shot record;b—real primaries;c—result of surface-related multiple suppression图3 某炮数据SRME方法对自由表面多次波的压制结果Fig.3 Suppression results of surface-related multiples of a shot record by SRME method

SRME方法理论简单，只需要通过地震记录的褶积即可预测自由表面多次波，不需要额外的先验信息，可在较少的迭代次数下得到较好的压制结果。但在实际资料处理中，预测的多次波信息并非十分准确，在匹配相减后仍会有部分多次波能量残余存在。且在一次波和多次波干涉时，不再成立的一次波能量最小的前提假设会造成有效信号的损伤。虽然人们提出不少的匹配方法，但这些方法只降低了对一次波能量的损伤程度，未从根本上解决此问题。

2.2 反数据域方法

SRME方法中一次波能量最小的前提假设从方法上导致了它的局限性，而反数据域方法的提出，提供了一种避免该假设的多次波压制思路，以相同的理论基础、不同的推导方式解决了SRME存在的问题。Berkhout[23]在2006年将SRME理论发展到了反数据域，从而避免自适应相减的这一过程。首先通过式(1)得到地震数据P显式表达式：

P=(I-X0R)-1(X0S) ，

(4)

对式(4)两边求逆，转换到反数据域：

P-1=(X0S)-1(I-X0R)=(X0S)-1-S-1R。

(5)

由式(1)、式(2)可知P0=X0S，A=S-1R，则原始地震记录P的逆为：

(6)

之后，许多学者对反数据域方法进行了扩展，Berkhout[24]于2006年在反数据域多次波压制过程中加入了反褶积处理方法。基于波动方程的数据驱动算法虽然在海洋地震多次波压制中相当有效，但在陆地数据中应用时需要特殊的预处理，Kelaims[25]在2006年实现了反数据域中陆地表面多次波的压制。Luo等[26]在2007年将反数据域方法发展到了层间多次波的压制上，给出了一维实例。Berkhout等[27]在2007年提出了一种利用反数据空间进行时移地震处理的新方法。同年，Berkhout等[28]实现了反数据域二维实际数据的层间多次波压制方法，并给出实例。马继涛等[29]在2009年介绍了两种矩阵求逆的方法，并给出了在二维情况下反数据域处理的实例。同年，马继涛等[30-31]结合平面波域多次波产生机理，实现了平面波域反数据处理压制自由表面多次波方法，在保证一次波能量不受损害的同时，更有效地达到了压制自由表面多次波的目的。

图4展示了反数据域法对比较复杂的水平层状模型正演数据的自由表面多次波压制结果[24]。由图4b可以看到在反数据域中自由表面算子聚焦在零时间处(如红框所示)，通过线性Radon变换，将此聚焦点转换为τ=0这条直线，然后通过简单切除直线τ=0上的数据达到去除自由表面算子信息。再通过转换到正数据域后(如图4c)可以非常直观地看到自由表面多次波得到有效压制。

a—原始单炮数据;b—反数据域的地震炮记录;c—反数据域法多次波压制结果a—original shot record;b—seismic shot record in inverse data domain;c—result of multiple suppression in inverse data domain图4 某炮数据反数据域法对自由表面多次波的压制结果Fig.4 Suppression result of surface-related multiple of a shot record in inverse data domain

与正数据域中的表达公式相比，反数据域公式更加简单明了。反数据域方式算法简单，不需要大量的矩阵相乘运算，具有较大的计算速度优势，该方法利用了一次波和多次波在反数据域的可分离性，通过简单的切除操作就达到了自由表面多次波压制的目的，避免了自适应相减过程，解决了SRME在预测相减时对一次波能量的损伤问题。但由于需要做矩阵求逆运算，求逆过程容易引入噪声影响结果，因此需要选取合适的矩阵求逆方法[30]。

2.3 EPSI方法

为了更准确地得到一次波信息，寻找一种既能避免一次波能量最小假设， 又无求逆过程的自由表面多次波压制方法是很有必要的，于是EPSI方法顺势而生。2009年，Van Groenestijn等[32]提出了稀疏反演一次波估计方法(EPSI)，将SRME中的一次波能量最小假设改变为使估计的一次波和多次波数据之和不断逼近地震原始数据的目的，即将地震数据和估计的一次波及多次波差值的能量作为目标函数J，使其最小化。

(7)

(8)

式中：ΔX0和ΔS表示为一次波响应X0和震源子波函数S的梯度值；H表示复共轭转置运算。根据X0和S的梯度值得到X0和S的更新值：

(9)

其中，α和β为通过ΔX0和ΔS对应搜索步长，求取公式如下:

其中，

(10)

为约束X0的反演过程，在时间域将稀疏性引入X0的反演过程中，从而加速目标函数的收敛。通过上式(8)～(10)可计算得到一次波响应X0和震源子波函数S的估计值，然后通过两者频率域相乘(或时间域褶积)得到一次波数据，若不满足要求，则继续迭代反演。可以看出对X0和S的估计值都对一次波估计值的准确度具有直接影响，反演过程中的多参数是影响EPSI方法速度与精度的重要因素。

随后，经过一系列改进，Van Groenestijn和Verschuur[33]在2010年将其发展到了被动源地震中，并实现了海底OBC数据EPSI一次波估计[34-39]。刘国昌等[40]在2013年实现了海底OBS数据EPSI一次波估计。Lin等[41]在2013年提出了L1范数下的EPSI一次波反演方法，提高了多次波压制结果的精度。为了提高多次波去除的准确性和效率，Lopez等[42]在2012年结合SRME和EPSI方法各自的优点，提出了一种SRME和EPSI联合的多次波压制方法。2013年，Ypma等[43]在此基础上扩展到了广义的EPSI方法，对海洋实际数据层间多次波进行了有效压制；Angarita等[44]在聚焦变换域实现了EPSI方法，重建了大范围缺失的地震数据；Feng等[45]提出L1范数约束下三维曲波变换稀疏反演一次波估计方法。程浩等[46]在2015年结合范数约束的准确性和优越性，提出了L1约束被动源数据稀疏反演一次波估计方法。

图5展示了EPSI方法对双层水平层状模型正演数据的一次波估计结果，并与SRME多次波压制结果进行对比[47]。EPSI对一次波的估计是通过X0和S进行褶积得到的，因此需要先对上述两者进行反演,再通过褶积得到估计的一次波信息(图5c)，通过图5b～c 对比，可以看到，SRME方法在一次波和多次波具有能量干涉处(如橘色箭头所示)会产生一次波能量的损失，并且在局部地方存在明显的多次波能量残余(如黑色箭头所示)。而EPSI方法通过反演的方法成功避免了上述问题，比较准确地反演出了一次波能量。

a—原始单炮数据;b—SRME多次波压制结果;c—EPSI压制结果a—original shot record;b—result of surface-related multiple suppression by SRME;c—result of surface-related multiple suppression by EPSI图5 某炮数据EPSI方法对自由表面多次波的压制结果Fig.5 Suppression results of surface-related multiples of a shot record by EPSI method

EPSI通过一个波形反演替换了SRME的预测和相减，一次波被直接反演得到，克服了SRME一次波能量损伤问题和反数据域的求逆过程。但是由于EPSI一次波脉冲响应的稀疏性假设，在存在波散等情况的复杂地下构造下无法选取准确的地下一次波响应初值，导致了不能得到令人满意的一次波估计效果[48]。并且EPSI需要额外的调优因子和较高的迭代次数，大大地增加了方法的计算量。

2.4 CL-SRME方法

SRME和EPSI各有优缺点，于是一种结合两种方法优势，避免双方劣势的自由表面多次波压制新方法被提出。Lopez等在2014年[22]提出了CL-SRME(Closed-loop SRME)方法，它和EPSI相似，也是通过波形反演的方式替代了SRME的预测相减，将传统的多次波预测方法变成了闭合循环的一次波反演方法。

为使估计的一次波及多次波之和尽可能接近地震记录，通过式(2)建立目标函数：

(12)

通过梯度值ΔP0对一次波P0进行第i次迭代数据更新：

(13)

其中α为搜索步长，公式如下：

其中，

(14)

式中：Re表示取实部运算；Tr表示求矩阵的迹，即求矩阵主对角线数据之和。表面算子A与SRME中的假设和计算相同。可以看到，CL-SRME是直接求取一次波的梯度和步长，对一次波信息不断更新，直至符合要求。也正由于该方法是直接对一次波数据进行操作的，所以与EPSI方法相比CL-SRME方法对一次波反演更准确且更稳定。

之后，Lopez等[49-50]在2015年将CL-SRME方法发展到了聚焦域，完成了部分数据缺失的地震记录的重构和一次波估计。Lopez等[51]在2015年将聚焦域CL-SRME发展到三维。Vrolijk等[52]在2016年将CL-SRME用于鬼波的压制。Wang等[53]于2017年提出了三维L1范数下稀疏约束反演CL-SRME方法。Vrolijk等[54]在同年实现了CL-SRME方法对鬼波和多次波的一体化压制。王铁兴[48]在2020年提出在被动源数据中的三维稀疏约束CL-SRME一次波反演方法和混采数据CL-SRME一次波反演方法。

图6展示了CL-SRME方法对盐丘模型正演数据的一次波估计结果，并与SRME、EPSI两种方法对多次波压制结果进行对比[22]。可以看出，3种方法在不同程度上对自由表面多次波具有压制效果，通过图6b～d的比较可见，CL-SRME方法估计的一次波信息明显更为准确，更加全面(如橙色箭头所示)，多次波能量的泄露也是最小的(如黑色箭头所示)。对于复杂构造，CL-SRME方法的适应性相比较EPSI来说更强。因此，CL-SRME方法是目前看来比较具有优势的自由表面多次波压制方法。

a—原始单炮数据;b—SRME多次波压制结果;c—EPSI多次波压制结果;d—CL-SRME多次波压制结果a—original shot record;b—result of surface-related multiple suppression by SRME;c—result of surface-related multiple suppression by EPSI;d—result of surface-related multiple suppression by CL-SRME图6 某炮数据CL-SRME方法对自由表面多次波压制的结果Fig.6 Suppression results of surface-related multiples of a shot record by CL-SRME method

CL-SRME一次波估计方法避免了SRME的一次波能量损伤和EPSI稀疏性问题，也不需要逆计算，并具有较少的调优因子，只需不断更新自由表面算子A就可以达到一次波估计的效果。且迭代次数较少，在稳定性和准确性上比EPSI都有所提高。

2.5 SRME、反数据域法、EPSI、CL-SRME方法对比

本文介绍了反馈迭代模型和在此基础上发展的4种常用自由表面多次波压制方法的原理及研究进展，通过模型示例(表1)展示了4种方法压制自由表面多次波的效果，也直观地看出了每种方法的优缺点。4种方法以地震数据为预算因子实现了对自由表面多次波的压制，不需要额外的先验条件，但都需要全波场信息，对观测系统和观测数据要求比较严格；SRME和反数据域方法理论简单，但SRME在一次波与多次波具有能量干涉时，匹配相减会对一次波能量产生损伤，因此，对于如何做到更准确的匹配相减是SRME的发展方向；反数据域法的求逆运算容易引入噪声，寻找合适的求逆方法则是反数据域法的重点研究内容；EPSI避免了上述两种方法存在的问题，但由于稀疏性的引入，对复杂构造的适应性不强，并且需要较多的调优参数和较高的迭代次数；CL-SRME方法则避免了前3种方法的问题，具有较好的稳定性和准确性。但由于EPSI和CL-SRME是基于波形反演的思路，其计算成本要高于SRME和反数据域法。因此，如何提高计算精度与计算效率是EPSI和CL-SRME的研究重点。

表1 基于反馈模型多次波压制方法优缺点Table 1 The advantages and disadvantages of multiple suppression methods based on feedback model

3 总结

本文从方法理论、发展现状、实例比较、优缺点分析等方面深入研究了基于反馈迭代模型的SRME、反数据域法、EPSI、CL-SRME四种自由表面多次波压制方法。可以看到四种方法各有利弊，在实际生产中可结合不同方法的优缺点选择合适的自由表面多次波压制方法。对于成本问题，SRME和反数据域法是比较理想的可考虑方案。而在多次波压制精度上，EPSI和CL-SRME方法更具优势。当然，在已有方法的基础上寻求可以兼顾成本与精度的新方法一直是多次波压制方法研究的追求目标。