多次波是地震数据处理,特别是海洋地震数据处理中最突出的问题之一.多次波与有效波相互混叠会干扰对有效波的正确识别,降低地震资料的信噪比和分辨率,影响地震偏移的精度.多次波压制方法主要有滤波方法和基于波动方程的预测相减法[1].滤波方法主要利用多次波的周期性和一次波与多次波之间的可分离性[2-5];基于波动方程的预测相减法首先通过模拟实际波场或者反演地震数据估计出多次波,然后从原始地震信号中将预测多次波自适应减去,从而得到一次波[6-9].因此,基于波动方程的预测相减法分为多次波预测和多次波自适应相减两个步骤.其中,多次波自适应相减的主要方法有基于匹配滤波器设计的方法[10]、基于模式的方法[11-16]和基于盲信号分离的独立成分分析的方法[17-21].
基于匹配滤波器设计的多次波自适应相减方法,通过优化与一次波有关的目标函数,得到有限长的匹配滤波器,从而利用预测多次波逼近真实多次波,达到压制多次波的目的.常用的准则有基于一次波能量最小化的方法[22],此方法隐含了一次波和多次波的正交性假设.因此,在一次波和多次波相互正交时,能够很好地压制地震数据中的多次波;而在非正交时,则常常会损伤有效信号.因此,为了克服这一问题,基于增广匹配滤波器的一次波L2 范数最小化方法[23-24]、基于一次波L1 范数最小化的方法[25-26]和非平稳正则化L2范数方法[27]等被相继提出.
基于模式的多次波自适应相减方法主要利用地震信号的一次波和多次波各自具有不同的空间模式,比如一次波和多次波的预测误差滤波器系数不同或者一次波和多次波的曲线拟合系数不同等等,以此作为求解一次波和多次波的约束条件,并构造与预测误差或拟合误差相关的优化目标,从而得到多次波压制结果[11-16].
基于盲信号分离的ICA 方法[17-21],将多次波自适应相减问题表示为多个信号瞬时混合的盲分离问题,利用几何ICA 的方法估计得到混合矩阵,并根据地震信号的超高斯性构造优化目标,实现一次波和多次波的分离.基于盲信号分离的ICA 方法利用了地震信号的超高斯性,避免了基于一次波能量最小化方法需要采用的多次波和一次波正交的假设.需要指出的是,上述方法[17-21]的实现过程包括子波差异消除和ICA 分离两步,且所采用的混合ICA 模型只考虑预测和实际多次波在时间方向的差异.
考虑到预测和实际多次波不只在时间上,而且在空间上也存在差异,本文提出了一种基于多道卷积信号盲分离的多次波自适应相减方法.在所提出的方法中,利用2D 卷积模型,表示预测多次波和实际多次波之间的差异,并采用分离出的一次波信号的非高斯性最大化作为优化目标,估计2D 卷积核系数.和传统的基于L2 范数的多道匹配滤波器相比,所提出的方法采用了非高斯性最大化准则,避免了多次波和一次波正交的假设.和基于瞬时混合ICA 的多次波自适应相减方法不同,基于卷积信号盲分离的多次波自适应相减方法在实现一次波和多次波分离的过程中,同步消除预测和实际多次波在时间-空间上的差异.
下面首先介绍基于卷积信号盲分离的多次波自适应相减方法的基本原理,在给出优化目标的基础上推导迭代求解算法,最后给出模型数据和实际数据的处理结果.
2 基本原理基于滤波器设计的多次波自适应相减方法通常要对地震数据进行分块处理,并假定同一分块数据中的预测多次波与真实多次波之间,具有相同的尺度差异、子波差异、时延差异等,因此同一分块中的数据只需要一个滤波器就能实现预测多次波和真实多次波之间的匹配.此时多次波自适应相减模型可以表示为[10]
(1) |
其中,s为原始地震信号,p0 为真实一次波,m是真实多次波,$\tilde m$为预测多次波,* 表示卷积算子,h表示预测多次波与真实多次波之间差异的滤波器.这里,h可以是单道算子,也可以是多道算子.根据公式(1)表示的模型,多次波自适应相减问题就是如何设计匹配滤波器f来消除预测多次波与真实多次波之间的差异:
(2) |
得到滤波器f后,一次波的估计p可以由原始地震信号、预测多次波和匹配滤波器获得,即
(3) |
其中,$\tilde M$是由预测多次波构成的数据矩阵.
滤波器f设计的目标就是使预测多次波经过f滤波后与真实多次波之间的差异达到最小,从而使相减得到的一次波p的某个统计量最小.常用的基于能量最小化的多次波相减法就是假定一次波能量最小,即最小能量准则:
(4) |
由式(4)得到f的最小二乘解
(5) |
地震数据是超高斯分布的信号[28],因此利用卷积信号盲分离中常用的非高斯性最大化准则来求取滤波器f和地震数据的统计特性更吻合.原始地震信号是由一次波和多次波相加得到的混合信号,而且在同一分块中的一次波和多次波接近统计独立且均符合超高斯分布[17].根据概率论的中心极限定理,两个满足超高斯分布的相互独立随机变量的和,其分布比其中任何一个随机变量的分布都要更接近高斯分布[29].因此分块数据中的原始地震信号减去多次波后,所得到的一次波具有最大的超高斯性.负熵是盲信号处理中一种常用的衡量非高斯性的准则,具有较高的鲁棒性.根据负熵的定义,我们可建立如下的优化目标[29]:
(6) |
其中$\bar p$为对一次波估计p进行零均值化和方差归一化后得到的数据,ν 是一个零均值,单位方差的高斯分布的随机信号,E表示求均值,G为非线性函数,可以采用函数G的一种形式如下:
(7) |
可利用多种方法对式(6)中的优化目标求解.我们借鉴文献[4, 30]中的求解思路,得到类似迭代加权最小二乘(IRLS)的求解算法.
附录A 中给出了单道滤波器求解的详细推导过程,附录B 推导了多道滤波器的求解过程.针对多道匹配滤波器首先根据式(3)将滤波方程写为[10]
(8) |
其中,pi,j为滤波后得到的一次波估计,i= 1,2,…,T,T为时间方向分块的大小,j= 1,2,…,S,S为分块的最大道数,si,j为原始地震数据,$\tilde m$i,j为相应的预测多次波数据,fk,l为多道滤波器系数,2K+1为滤波器在时间方向的长度,2L+1为滤波器在空间方向的道数.
将式(8)代入优化目标式(6)中并求导,最终可以得到如下方程(附录B):
(9) |
其中,矩阵A中的第(k+K)×(2L+1)+(l+L+1)行,第(k1 +K)× (2L+1)+ (l1 +L+1)列元素ak,k1,l,l1为
列向量b中第(k+K)× (2L+1)+ (l+L+1)行元素bk,l为
列向量f中第(k+K
在公式(9)中,矩阵A和列向量b都需要知道一次波估计p.而一次波估计p需要通过滤波器f得到.这里我们采用迭代最小二乘方法得到滤波器f的估计:
(1) 设置最大迭代次数,并给定滤波器f的初始值;
(2) 由式(8)得到当前滤波器下的一次波估计p,若一次波估计p的能量为零,则停止迭代;
(3) 由p,f和s构造式(9);
(4) 求解式(9)得到更新后的滤波器f,然后转入步骤(2)进行下一轮的迭代,直到滤波器f的系数基本不变或者迭代次数达到最大迭代次数.
在第(1)步中,通常滤波器的初始值可以设为随机数,或者采用基于能量最小化的方法得到初始滤波器系数.第(2)步中为了减少迭代求解的次数,可以给定一个阈值,使得当滤波后的一次波能量的变化小于给定的阈值时,就停止迭代.方程(9)可以通过直接矩阵求逆求解f=A-1b,也可采用Leivensen递归求解方程(ATA)f=ATb,或者采用共轭梯度法进行求解,以减少由于矩阵求逆运算而带来的计算量问题.在下面的实验中,卷积信号盲分离方法的最大迭代次数设为5次,就可以得到满意的结果,所以,和采用同样滤波器参数的能量最小方法相比,所提出方法的计算量大概是前者的5倍.
4 应用实例为了验证所提方法的有效性,我们处理了人工合成数据和实际数据,并和基于能量最小化的多道匹配滤波器方法进行了比较.
首先处理了一个简单的人工合成模型数据.模型数据是由一个一次波和一个多次波组成的共中心点道集,用来合成地震数据的子波在空间域上是变化的.时间采样率是4ms.预测多次波由下式得到:
(10) |
其中,$\tilde m$为预测多次波,m为真实多次波.这里,预测和实际多次波之间存在时间和空间上的差异.
图 1a显示了原始地震数据,图 1b为预测多次波,图 1c和图 1d分别是单道和多道滤波器的能量最小方法得到的一次波估计.图 1e和图 1f分别是单道和多道滤波器的卷积信号盲分离方法得到的一次波估计.在这个实验中,单道滤波器参数为K=2,多道滤波器参数为K=2,L=1.数据分块参数为100采样点和40 道.可以看出,当预测和实际多次波之间存在空间上的差异时,单道的方法得不到满意的结果.需要指出的是,单道滤波器可以利用更长的滤波器(即参数K取更大的值)来减少多次波残余,但会造成一次波的畸变.多道滤波器的采用,可以在时间方向使用较少的滤波样点,从而能更好地保护一次波.
为了定量评价基于卷积信号盲分离的多次波自适应相减方法的性能,采用如下定义的信噪比:
(11) |
其中,p0 是真实一次波,p是估计出的一次波.表 1给出了简单合成模型数据的信噪比结果,从图 1 和表 1可以看出,采用单道滤波器和多道滤波器时,本文所提出的方法的性能都优于基于能量最小的方法.
将卷积信号盲分离方法应用于Pluto数据,同时考虑到预测多次波和真实多次波之间可能具有空间上的差异,因此采用多道滤波器的设计方法.对整个数据采用分块处理,以保证其满足匹配滤波器设计的一致性假设条件.作为对比,多道能量最小化方法也用来处理同样的数据.在本实验中,两种方法采用同样的处理参数.每块数据中包含30 道,每道包含40个采样点.多道滤波器参数为K=7,L=1.图 2是对Pluto数据分别用基于能量最小化的方法和本文方法得到的结果.图 2a是原始地震数据,图 2b是预测出的多次波.图 2c是基于能量最小化方法得到的一次波估计,图 2d是本文方法得到的一次波估计.图 3给出了图 2 所给出结果的局部放大图,图 3a和图 3b分别是图 2c和图 2d左侧方框内数据的放大图,图 3c和图 3d是图 2c和图 2d右侧方框内数据的放大图.从图 3 中箭头所指区域可以看出,本文方法的处理结果较基于能量最小化方法,能够更好地压制多次波,且估计出的一次波具有更高的连续性和信噪比.
我们将本文算法应用于一个海上数据集.在本实验中,每块数据中包含30道,每道包含40个采样点.多道滤波器参数为K=3,L=1.图 4a是原始的共偏移距道集(偏移距为208m),图 4b是预测出的多次波,图 4c和图 4d分别是用某商业软件的多道基于能量最小化方法和本文方法进行自适应相减所得的结果,黑色箭头所指区域表明本文方法相对于能量最小化方法可更好地压制多次波,白色箭头所指区域表明本文方法可更好地保护一次波,压制结果中一次波具有更好的连续性.
在基于匹配滤波器设计的多次波自适应相减问题中,原始地震信号可以表示为一次波和预测多次波的卷积混合,因此多次波自适应相减问题可以看作是从两个观测信号中分离出一个源信号的卷积盲分离问题.非高斯性最大化是盲信号分离中一个重要准则,认为独立非高斯性信号在混合后的非高斯性要低于混合前信号的非高斯性.本文中,我们将盲信号分离理论中的非高斯性最大化技术和基于滤波器设计的多次波自适应相减技术相结合,提出一种基于多道卷积信号盲分离的多次波自适应相减算法.由于避免了采用一次波与多次波的正交性假设,因此与基于能量最小化的方法相比,所提出的方法能够更有效地压制多次波和保留一次波的能量,使单道的一次波估计结果具有更高的信噪比和连续性.对人工合成数据和实际数据的处理结果表明了算法的有效性.
一般预测多次波和真实多次波之间可能会存在空间上的差异,仅利用单道滤波器的设计方式很难校正空间差异,因此我们推荐使用多道滤波器进行多次波的自适应相减.多道滤波器参数K,L的取值会影响处理效果,对于某一待处理数据,我们首先利用几个道集进行不同参数的实验,然后选取效果好的参数来处理整个数据体.非线性函数G的选取对衡量非高斯有重要的作用,本文中选取的是负指数函数.需要指出的是,非线性函数G也可以选取双曲函数、对数双曲函数等多种类型的函数.
附录A对于单道滤波器,滤波方程为
(A1) |
其中,pi为一次波估计,si为原始地震信号,$\tilde m$i为预测多次波.fk是单道滤波器,2K+1为滤波器在时间方向的长度.
由式(6)中的优化目标
(A2) |
对f求导,得
(A3) |
由于p为超高斯信号,v为高斯信号,因此
(A4) |
令
(A5) |
将上式展开,并求导:
(A6) |
由非线性函数式(7)得
(A7) |
另一方面,由滤波方程式(A1)得d
(A8) |
由于地震数据的均值近似为零,因此我们通常只对一次波估计p进行方差归一化后得到珚p,σ表示一次波估计p的方差:
(A9) |
于是有
(A10) |
将式(A7)和式(A10)代入式(A6),于是链式求导为零的结果变为
(A11) |
再将式(A1)代入上式,得到
(A12) |
将上式改写为
(A13) |
整理上式就能得到如下方程组
(A14) |
其中
附录A 中给出了单道滤波器的求解算法,多道滤波器的求解算法与此类似,只是迭代求解的方程组中的元素有所不同而已,附录A 中1D 信号变为2D 信号,如pi,$\tilde m$i分别变为pi,j和$\tilde m$i,j.
对于多道滤波器,滤波方程为
(B1) |
其中,pi,j为滤波后得到的一次波估计,si,j为原始地震数据,$\tilde m$i,j为预测多次波,i=1,2,…,T,T为时间方向分块的大小,j=1,2,…,S,S为分块的最大道数.fk,l为多道滤波器系数,2K+1为滤波器在时间方向的长度,2L+1为滤波器在空间方向的道数.
由式(6)中的优化目标
(B2) |
对f求导,得
(B3) |
由于珚p为超高斯信号,v为高斯信号,因此
(B4) |
令
(B5) |
将上式展开,并求导:
(B6) |
由非线性函数式(7)得d
(B7) |
另一方面,由滤波方程式(B1)得
(B8) |
由于地震数据的均值近似为零,因此我们通常只对一次波估计p进行方差归一化后得到p,即
(B9) |
于是有
(B10) |
将式(B7)和式(B10)代入式(B6),于是链式求导为零的结果变为
(B11) |
再将式(B1)代入上式,得到整理上式可得如下方程组:
(B12) |
将上式改写为
(B13) |
整理上式可得如下方程组:
(B14) |
其中,矩阵A中的第(k+K)×(2L+1)+(l+L+1)行,第(k1 +K)× (2L+1)+ (l1 +L+1)列元素ak,k1,l,l1为
列向量b中第(k+K)× (2L+1)+ (l+L+1)行元素bk,l分别为
列向量f中第(k+K)× (2L+1)+ (l+L+1)行元素为fk,l.
[1] | Weglein A B. Multiple attenuation: an overview of recent advance and the road ahead. The Leading Edge , 1999, 18(1): 40-44. DOI:10.1190/1.1438150 |
[2] | Foster D J, Mosher C C. Suppression of multiple reflections using the Radon transform. Geophysics , 1992, 57(3): 386-395. DOI:10.1190/1.1443253 |
[3] | Lu W K, Zhang X G, Li Y D. Multiple removal based on detection and estimation of localized coherent signal. Geophysics , 2003, 68(2): 745-750. DOI:10.1190/1.1567244 |
[4] | Liu J, Lu W K. An improved predictive deconvolution based on maximization of non-Gaussianity. Applied Geophysics , 2008, 5(3): 189-196. DOI:10.1007/s11770-008-0027-1 |
[5] | 薛亚茹, 陈小宏, 陆文凯. 压制多次波的正交多项式谱减法. 地球物理学报 , 2009, 52(3): 817–823. Xue Y R, Chen X H, Lu W K. Orthogonal polynomial spectrum subtraction for multiple attenuation. Chinese J. Geophys. (in Chinese) , 2009, 52(3): 817-823. |
[6] | Wiggins J W. Attenuation of complex water-bottom multiples by wave-equation-based prediction and subtraction. Geophysics , 1988, 53(12): 1527-1539. DOI:10.1190/1.1442434 |
[7] | Weglein A B, Gasparotto F A, Carvalho P M, et al. An inverse-scattering series method for attenuating multiples in seismic reflection data. Geophysics , 1997, 62(6): 1975-1989. DOI:10.1190/1.1444298 |
[8] | Verschuur D J, Berkhout A J, Wapenaar C P A. Adaptive surface-related multiple elimination. Geophysics , 1992, 57(9): 1166-1177. DOI:10.1190/1.1443330 |
[9] | Ma J T, Sen M K, Chen X H. Free-surface multiple attenuation using inverse data processing in the coupled plane-wave domain. Geophysics , 2009, 74(4): V75-V81. DOI:10.1190/1.3137059 |
[10] | 庞廷华. 基于匹配滤波器的多次波自适应相减. 北京: 清华大学自动化系, 2009 . Pang T H. Adaptive multiple subtraction based on matching filters (in Chinese). Beijing: Automation Department of Tsinghua University, 2009 . |
[11] | Manin M, Spitz S. 3D attenuation of targeted multiples with a pattern recognition technique. //57th Annual International Meeting, European Association of Geoscientists and Engineers, Extended Abstracts. 1995, B046. |
[12] | Guitton A, Cambois G. Prestack elimination of complex multiples: A Gulf of Mexico subsalt example. 68th Annual International Meeting, SEG, Expanded Abstracts , 1998, 17: 1329-1332. |
[13] | Brown M, Clapp R G. T-x domain, pattern-based ground roll removal. //70th Annual International Meeting, SEG, Expanded Abstracts. 2000: 2103-2106. |
[14] | Guitton A. Coherent noise attenuation using inverse problems and prediction-error filters. First Break , 2002, 20(3): 161-167. DOI:10.1046/j.1365-2397.2002.00246.x |
[15] | Fomel S. Applications of plane-wave destruction filters. Geophysics , 2002, 67(6): 1946-1960. DOI:10.1190/1.1527095 |
[16] | Guitton A. Multiple attenuation in complex geology with a pattern-based approach. Geophysics , 2005, 70(4): V97-V107. DOI:10.1190/1.1997369 |
[17] | 陆文凯, 骆毅, 赵波, 等. 基于独立分量分析的多次波自适应相减技术. 地球物理学报 , 2004, 47(5): 886–891. Lu W K, Luo Y, Zhao B, et al. Adaptive multiple wave subtraction using independent component analysis. Chinese J. Geophys. (in Chinese) , 2004, 47(5): 886-891. |
[18] | Lu W K, Mao F. Adaptive multiple subtraction using independent component analysis. The Leading Edge , 2005, 24(3): 282-284. DOI:10.1190/1.1895313 |
[19] | Lu W K. Adaptive multiple subtraction using independent component analysis. Geophysics , 2006, 71(5): S179-S184. DOI:10.1190/1.2243682 |
[20] | 李艳东. 约束盲信号分离算法及应用研究. 北京: 清华大学自动化系, 2006 . Li Y D. Studies on algorithm and application of constrained blind source separation (in Chinese). Beijing: Automation Department of Tsinghua University, 2006 . |
[21] | Lu W K, Liu L. Adaptive multiple subtraction based on constrained independent component analysis. Geophysics , 2009, 74(1): V1-V7. DOI:10.1190/1.3000600 |
[22] | Luo Y, Kelamis P G, Wang Y. Simultaneous inversion of multiples and primaries: inversion versus subtraction. The Leading Edge , 2003, 22(9): 814-818. DOI:10.1190/1.1614151 |
[23] | Wang Y H. Multiple subtraction using an expanded multichannel matching filter. Geophysics , 2003, 68(1): 346-354. DOI:10.1190/1.1543220 |
[24] | 李鹏, 刘伊克, 常旭, 等. 均衡拟多道匹配滤波法在波动方程法压制多次波中的应用. 地球物理学报 , 2007, 50(6): 1844–1853. Li P, Liu Y K, Chang X, et al. Application of the equipoise pseudomultichannel matching filter in multiple elimination using wave-equation method. Chinese J. Geophys. (in Chinese) , 2007, 50(6): 1844-1853. |
[25] | Guitton A, Verschuur D J. Adaptive subtraction of multiples using the L1-norm. Geophysical Prospecting , 2004, 52(1): 27-38. DOI:10.1046/j.1365-2478.2004.00401.x |
[26] | Pang T H, Lu W K, Ma Y J. Adaptive multiple subtraction using a constrained L1-norm method with lateral continuity. Applied Geophysics , 2009, 6(3): 241-247. DOI:10.1007/s11770-009-0028-8 |
[27] | Fomel S. Adaptive multiple subtraction using regularized nonstationary regression. Geophysics , 2009, 74(1): 25-33. |
[28] | Walden A T. Non-Gaussian reflectivity, entropy, and deconvolution. Geophysics , 1985, 50(12): 2862-2888. DOI:10.1190/1.1441905 |
[29] | Hyvarinen A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks , 1999, 10(3): 626-634. DOI:10.1109/72.761722 |
[30] | Liu J. Echo cancellation using predictive deconvolution based on high order statistics. //International Congress on Image and Signal Processing. 2008: 324-327. |