1 引言

目前，复杂地下介质、复杂岩性油气藏成为世界石油天然气的勘探目标，为了更准确地圈定有利目标区块和对地下介质进行清晰照明，叠前去噪预处理，特别是多次波等相关噪声的压制，成为至关重要的研究课题.其原因在于:速度分析、CMP 叠加以及偏移等技术都是基于地震数据只含有反射波的假设条件.多次波的存在会影响速度模型的建立、地震成像效果、反演和解释过程中的介质属性提取等^[1].

为了有效压制多次波而保留有效波，自20世纪50年代以来，多次波处理技术经历了从一维到二维再到三维，从水平介质到复杂介质，从基于几何射线到基于波动方程，从需要先验信息到不需要先验信息，从消除到利用的过程.其压制方法可以归结为两大类:基于信号分析的直接滤波方法^[2-5]和基于波动方程的预测相减方法^[6-9].

基于波动方程的预测相减法^[10-18]始于20 世纪70年代.Morley^[5]、Wiggins^[6]、Berryhill研究和实现了利用波场延拓外推方法来压制多次波；Berkhout和Verschuur提出反馈迭代方法；Cervalho、Coates、Weglein^[7]、金德刚等^[8]提出或完善了基于点散射模型的逆散射级数法.这些方法能适用于复杂的地下结构，需要较少或不需要关于地下结构的假定，可以预测出所有类型的多次波.但是，由于子波特性的差异，预测出的多次波记录和原始记录中的多次波存在振幅、相位、到时上的差异，需要对两者进行自适应匹配滤波才能达到真正消除多次波的目的.

针对自适应匹配相减，各种方法被提出:时间域维纳滤波方法、频率域最小能量优化方法、基于独立变量分析方法、预测误差滤波及基于模式匹配滤波等.时间域维纳滤波具有操作简单、效率高等优点，但是需要满足多次波和一次波的正交性条件，基于2范数匹配还需要满足大值条件；频率域最小能量优化方法和时间域维纳滤波方法的物理机理是一样的，都是应用误差能量最小，但其优化搜索速度较慢，且涉及到非线性化求解问题，效率较低；基于独立分量分析是一种盲信号分离技术^[19-20]，刘喜武等^[21]、陆文凯等^[22]采用了输出信号非高斯性最大准则，并利用高阶统计量表征非高斯性，改善了正交性约束条件，但是该方法建立在两者没有时移的假设基础上；基于模式匹配^[23]的方法不需要进行近偏移距道外推，但是对于多次波和一次波具有相同的时移及横向振幅变化时，是失效的.

由于基于时间域维纳滤波方法在多次波自适应匹配滤波中是应用最广泛的一种方法，因此，对该方法进行了改进.时间域维纳滤波方法中，基于2范数最小准则是最常见的模式，人们采取各种手段试图对它的大值条件和正交性条件进行改善.李鹏等^[24]、刘华峰^[25]在常规2范数匹配滤波中引入均衡拟多道的概念，但是，由于基于2 范数最小准则，没有克服2范数自身的大值条件约束，当多次波和一次波能量相差较大时，存在求解误差.为了避免大值条件，Guitton^[26]提出单道1范数最小自适应匹配滤波方法，虽然1范数对大的异常保持稳健，但是单道1范数最小自适应匹配滤波方法仍然对一次波和多次波正交的假设存在一定依赖.本文在前人基础上，基于1范数最小准则进行计算，加入均衡的统计效应及拟多道对相位的修正作用，对常规1 范数的正交性条件进行改善.该方法不需要满足大值条件，只需要满足加权后的一次波和加权后的多次波道及其变换道正交即可，最大程度上提高匹配效果.

2 方法原理

基于波动方程多次波压制方法中，将原始地震记录和预测出的多次波记录进行相减，就可以得到多次波压制后的一次波记录:

(1)

其中，d(t)表示原始地震记录，d₀(t)表示多次波压制后的一次波记录，m(t)表示预测出的多次波记录.

但是，由于地震子波的不一致性，基于上述方法预测得到的多次波记录和原始记录中的多次波是不匹配的，在到达时间、振幅、相位上存在差异，因此，需要对两者进行自适应匹配滤波.设s(t)表示自适应匹配滤波器，则自适应匹配的一次波记录在时间域表示为

(2)

2.1 2范数匹配滤波法

按照消除多次波后的地震记录2范数最小能量准则，自适应滤波器s(t)的求取是通过最小化如下的目标函数得到的:

(3)

设矩阵

则上式向量褶积可写为矩阵乘向量形式:

(4)

上式对s求偏微商，并令其等于零，可得:

(5)

(6)

通过(6)式就可以计算得到自适应滤波器.而(5)式相当于:

(7)

即要求多次波和有效波满足正交性；

同时可以看到，2 范数最小要求去除多次波后的地震记录满足能量最小准则，即大值条件^[4-5].在地震记录中的多次波比其附近的有效波能量弱很多的情况下，消除多次波后的地震记录不满足能量最小，从而自适应滤波器s的求取不准确，多次波的压制就会带来很大的误差.

2.2 1范数匹配滤波法

基于对2范数大值条件的改善，本文采用1范数最小准则进行匹配滤波:

(8)

对于上式的求解，通常的优化方法如牛顿迭代及共轭梯度法等，都是基于目标函数一阶连续可导.本文采用Huber于1973提出来的方法近似求解1范数最小值.该方法将求解1 范数转换为迭代加权的2范数最小.Bube^[27]、李鹏^[24]，李学聪^[28]等就该方法进行过讨论.

(9)

其中，W为加权矩阵:

(10)

ε为先验值，r_i(i=1，…，N)为第i采样点压制多次波后的剩余值.

利用能量最小化原则，对于任意给定的r_i，可以得到:

(11)

其中，N是每个地震道的采样点数.从(11)式可以看出:当r_i比较小时，求解的是2范数最小，而当r_i比较大时求解的是1范数解，它们的过渡点为ε.

对(9)式两边求s的偏导数并令其零，则转化为求解下面的方程组:

(12)

变形可得:

(13)

上式相当于:

(14)

其中，$\tilde{d}$₀(t)是压制多次波后的一次波d₀(t)的加权，·表示内积，(14)式表示多次波和加权后的一次波的内积为零，即单道1 范数最小自适应匹配滤波要求多次波和加权后的一次波是正交的.

上面的1范数匹配方法，是计算加权后的2范数，当多次波和一次波能量差异较大时，求解2范数最小；否则，求解的是1 范数最小.因此该滤波方法不存在滤波后的大值约束条件；另外，1范数准则只要求加权后的多次波和加权后的一次波正交.因此，1范数滤波方法大大改善了2 范数的约束条件，使得求解更加稳健准确.

2.3 改进的1范数匹配滤波法

尽管1范数匹配滤波方法改善了2范数准则下的正交性，但当这种加权后的正交性不满足时，求解依然出现误差.本文提出改进的1 范数匹配滤波方法是对参与匹配的地震道进行空间上的均衡，使得地震记录的非高斯性增强.另外，对输入的多次波模型m(t)进行Hilbert变换及求导运算，进行多个角度的相位旋转.输入多道原始地震记录及相应多次波记录的变换道数据进行匹配，具体计算如下:设参与匹配的地震道数为k，则

(15)

其中，m_i为第i个参与匹配的多次波模型记录；m_i为第i个参与匹配的多次波模型的求导道；m_i^H为第i个参与匹配的多次波模型的Hilbert变换道；m_i^H为第i个参与匹配的多次波模型的Hilbert变换导数道；s₁，s₂，s₃，s₄ 分别为1范数匹配自适应滤波器；d_0i表示第i个压制多次波后得到的记录.

为了方便运算，将(15)式写为矩阵相乘形式:

(16)

根据1范数最小准则，求解自适应滤波器s₁，s₂，s₃，s₄ 的目标函数为

(17)

其中，，j=1，2，…，N，N为每道参与匹配运算的采样点数；r_i为匹配后的第i道一次波记录.

通过对自适应滤波器s₁，s₂，s₃，s₄ 分别求偏微商为零，利用莱文森方法就可进行求解.

(18)

其中，

(19)

(20)

n=1表示多次波记录的加权计算值；n=2 表示多次波导数道的加权计算值；n=3 表示多次波Hilbert变换道的加权计算值；n=4 表示多次波Hilbert变换导数道的加权计算值.

该改进后的方法只需要满足:

(21)

(21) 式相当于:

(22)

其中，m_0i为第i个匹配后的多次波模型记录；m_0i为第i个匹配后的多次波模型的求导道；m_i^H0为第i个匹配后的多次波模型的Hilbert变换道；m_i^H0为第i个匹配后的多次波模型的Hilbert变换导数道.

(22) 式表示同一时窗内匹配加权后的一次波和匹配加权后的多次波、匹配加权后的多次波Hilbert变换道、匹配加权后的多次波求导道、匹配加权后的多次波Hilbert变换求导道内积之和为零.由此可以看出，该方法进一步弱化了其对多次波和一次波正交性的要求.

3 实例应用分析

为了验证本文改进的1范数匹配滤波方法的应用效果，对Pluto模型数据和实际南海深水盆地海洋数据进行了计算，并与2 范数等滤波结果进行了对比.

3.1 模型数据

选取Pluto模型进行试算.该模型是SMAART组织设计的一个具有复杂地下构造的模型，包含起伏较大的海底、较大反射系数的盐丘构造界面，具有非常发育的表面多次波，是SEG 指定用来验证多次波压制算法的测试模型.图 1 显示了经过不同的匹配滤波方法压制多次波之后的对比效果.

总体而言，经过各种方法匹配滤波后，多次波得到很大程度的压制，被掩盖的一次波得到显露.对比而言，图 1c中的箭头所示位置的多次波同相轴由于大值条件的限制没有得到彻底的压制，仍然有较为明显的残留痕迹；相比于常规2 范数，图 1d中经过1范数匹配滤波后的结果有所改善，尤其是箭头所示位置处多次波同相轴能量得到较大消减.图 1e中，无论是箭头所示还是白色线所示位置处的多次波，都得到了很好的压制，而且，和多次波相交的一次波能量也得到很好的保持.

对比图 1c、1d与1e, 常规匹配滤波由于正交性条件的约束，复杂地区形成的多次波不能得到很好的压制；另外，由于改进后的1范数方法克服了大值条件的约束并对正交性条件进行了弱化处理，能够对整个记录中出现的复杂多次波都进行合理有效压制的同时，合理有效地保护一次波.因此，本文的方法在复杂构造条件下对比于其他方法具有一定的优越性.

3.2 实际数据

该数据为南海深水盆地某二维海洋地震记录，其海水深度变化范围较大，从300~1800m起伏变动，地下构造复杂，表面相关多次波发育强烈.图 2a中出现明显的一阶和二阶表面多次波能量.由于排列长度较大，远偏移距处多次波周期特性降低，预测反褶积方法无法达到较好的多次波压制效果.图 2b为利用改进后的1范数匹配滤波后的记录，直剪头所指位置处的多次波能量得到彻底压制，而弯箭头所指位置处的一次反射波同相轴得以显现，这一点在图 3中的局部对比中表现得更为明显，一阶、二阶表面多次波能量几乎看不到.从图 4 中的速度分析看出，多次波速度较低(椭圆框标定)；滤波后，椭圆框中的能量团几乎不存在，证明多次波压制成功.图 5a和5b为多次波压制前后叠加剖面对比，可以看出，该地区地下结构复杂，存在较多的小型起伏构造，利用本文方法进行多次波匹配滤波压制后，多次波的能量得到最大程度衰减，被多次波掩盖的一次波能量得到显现.

4 结论

本文从理论上讨论了改进后的1范数匹配滤波法对约束条件的改善，通过模型及实际资料的处理，在同其他滤波方法对比的基础上，验证了该方法的有效性和优越性.

从实际复杂地下构造的资料处理经验看出，该方法比常规2范数、常规1范数滤波法更有效.但同时对于一些质量较差的资料，该方法无法实现多次波的完全压制，会在记录中存在明显的痕迹.

将该方法和常规方法如预测反褶积、拉东变换等进行多方法的综合处理，取长补短，将会取得更好的压制效果.另外，在复杂地质条件下，多次波可以提供比一次反射波更好的对地下信息的照明，伴随相应理论的提出和计算机技术的发展，尝试将多次波进行有利信息的提取和利用，或者直接对多次波进行成像，都是有益的探索，是未来的研究方向.