受制于常规地震数据处理理论技术的发展,一般海上地震勘探常用的采集方式为单炮激发,且为防止单炮之间的地震信号相互干扰,不同单炮的激发需要有一定的时间间隔,这直接导致了海上地震数据采集的低效率性,尤其是在高密度、宽方位角等震源较密集的观测系统中。为提高地震数据采集效率,多震源激发采集技术逐渐应用到地震数据的采集中。早在1987年,Lynn等^[1]在海洋地震采集中将来自其他震源的干扰定义为交叉噪声,并通过叠加来压制交叉噪声;随后Beasley等^[2]提出在海洋采集中采用同步源激发可提高采集效率,采集过程中将气枪震源对称放置于海底电缆两端同时激发震源,在共中心点道集中应用倾角相关滤波分离混合波场;Bagaini^[3]证实多震源同步采集技术能更快更好地采集数据,该技术不仅能显著提高采集效率,还能通过更高的覆盖次数获取更高质量的地震数据。Ikelle^[4]在海洋采集中引入了震源编码和解码,并对海洋采集中使用瞬时震源进行编码和非编码的优劣进行讨论;Berkhout^[5]随后通过矩阵运算的方式来表示地震数据的混合过程和分离处理过程,并提出设计更好的混合采集方式以优化源分离并减少交叉噪声的干扰。Vaage、Hampson以及 Dragosset 等^[6,7,8]在第二震源的激发时间中引入随机延迟方式,将此激发方式采集得到的混合数据按照各个单震源的激发时间进行时间上的延迟校正,再从共炮点道集分选到共检波点域、共偏移距域和共中心点域等道集时,第二震源(次震源)的地震信号表现为随机干扰,只要去除这些噪声,就可实现多震源地震数据的分离。

影响多震源激发采集地震数据成像的关键因素与常规单炮激发技术有所不同,其能否将混合波场分离直接影响这一技术的应用效果。目前多震源分离的主要方法分为两种方法,一种是反演法,一种是去噪法。基于反演理论方法将多震源混合数据的分离转化为一个带约束的最优化问题,通过将地震数据转化到曲波域、小波域等稀疏域中来进行约束,不断迭代完成混合数据的分离。Akerberg等^[9]通过将混合数据分选到共中心点道集,将高分辨 Radon 变换作为稀疏性提升方法引进反演中进行分离多震源混合波场;Moore^[10]将时间域线性 Radon变换引进混合波场反演分离中; Abma 等^[11]提出一种 Fourier 变换稀疏约束反演混合数据分离方法,最终得到了高质量的分离数据体;Chen 等^[12]提出了一种迭代正则化约束反演分离混合波场的框架,通过对比 Seislet 域阈值、频率—波数域阈值、f-x 预测滤波3种分离算法在框架中的应用,得出 Seislet 域阈值算法的分离效果最优。去噪法主要通过基于迭代预测相减法来实现。迭代预测相减法首先预测的是干扰波场,即预测第二震源的主波场,然后从总的波场中减去便得到了所需的波场;Moore等^[13]使用最小二乘法得到伪分离道集,然后将伪分离道集分选到共炮检道集进行去噪;Huo等^[14]将伪分离炮集分选到共中心点道集,利用全方位中值滤波器压制混叠噪声;Mahdad等^[15]进行了迭代分离算法,使用阈值分离法提取混叠噪声,每一次迭代都从分离后的炮集记录中减去混叠噪声,最终得到分离结果。韩立国等^[16]经研究发现,Mahdad使用的分离法适用于混合度较低的混合采集数据,对高混合度混叠数据分离效果欠佳,并提出了在迭代中预测有效信号的方法,结合中值滤波和曲波变换实现了混合数据的分离;周丽等^[17]提出利用自适应迭代多级中值滤波预测有效信号,对海上3种不同观测系统下的多震源混合波场进行分离。

文中采用一种改进的迭代预测去噪法,同时引入多方向矢量中值滤波(MD-VMF)来去除混叠噪声,通过理论模型的合成混合记录和海上实际数据的分离来验证本文方法的有效性和正确性。

多震源混合数据可用矩阵可表示为:

其中,P_blend为多震源混合数据,P为单震源激发所得数据,Γ为震源随机激发时间算子,可采用随机时间激发的方式。假如一个混合炮中震源个数为N_s个,则式(1)可写为:

其中,T表示转置矩阵。多震源混合数据的分离问题就是由混合数据P_blend求解单震源数据P的过程。由于混合炮中的震源个数至少为两个,即单炮记录个数多于混合数据个数,使方程(2)求解成为一个严重的欠定问题。基于最小二乘原理^[13],震源激发算子Γ的最小二乘逆可表示为:

其中,H表示复共轭变换。当混合算子采用随机时间延迟激发时,Γ^HΓ=I,I为单位矩阵,则式(3)可以写为:

P'即为分离后的炮集记录,该分离数据并不是准确的单震源数据,因此可将其称之为伪分离记录^[18]。从式(4)可以看出,混合数据的分离是将混合炮数据复制N_s次,通过对每份混合数据按照各个单震源的激发时间进行时间上的延迟校正,再从共炮点道集分选到共检波点域、共偏移距域和共中心点域等道集时,第二震源(次震源)的地震信号表现为随机干扰,只要去除这些噪声,就可实现多震源地震数据的分离。图1为常规单震源激发采集与多震源混合激发采集得到的某一共偏移距道集的对比。

迭代预测相减法在每一次迭代中预测交叉噪声,不断的从混合数据中减去,最终交叉噪声逐渐减少,信噪比提高,完成混合数据的分离。该迭代过程可以用表示为^[15]:

其中,Pⁱ为第i次分离的结果,F为去噪算子。式(5)右边第二项(Γ^HΓ-I)Pⁱ刚好为次震源信号引起的交叉噪声,将此从混合数据中减去,然后再进行滤波处理。笔者对式(5)右边进行改进,迭代过程可表示为:

其中,N_s为混合度,即一个混合炮里面的震源个数。式(6)与式(5)相比较,每次迭代时滤波器作用的地震信号分成两部分,一部分是上一次分离结果Pⁱ,一部分是减去随机交叉噪声的结果,两个部分的加权系数由混合度决定。

中值滤波被广泛应用于地震随机噪声的去除中。Liu等^[19]提出将矢量中值滤波(VMF)应用于向量值函数常遇到的地球物理问题,并证明矢量中值滤波(VMF)比传统标量中值滤波(SMF)更适合处理地球物理问题。文中使用的是由VMF扩展的多方向矢量中值滤波(MD-VMF)来实现混合地震数据分离。标量中值滤波(scalar median fliter,SMF)

可定义为一般情形下对标量序列进行排序处理,也可按最小距离和定义,定义如下:

式(7)为对一标量数列x_i∈{x_i|i=1,2,…,N}进行标量中值滤波(SMF)距离定义,D(x_s)即为滤波结果与每一个元素x_i的距离和,当其值最小时x_e= argminxsD(x_s),x_e为滤波结果。

在标量中值滤波的基础上,发展成矢量中值滤波(vector median filter,VMF),可将其定义为应求矢量为与其他矢量距离和最小的矢量,其距离定义如下:

其滤波结果为:X_e= argminXsD(X_s);式(8)较式(7)相比将标量换成矢量,可以说标量中值滤波为矢量中值滤波的特例。在矢量中值滤波中,当多个矢量满足这一条件时,需定义其他输出条件。

多方向矢量中值滤波基于矢量中值滤波,由Huo等^[14]于2009年提出,并将多方向矢量中值滤波用于压制混叠噪声以实现混合波场分离。该方法中,用滑动时窗截取部分地震数据获得一组矢量值,所截数据可被定义为矢量元素。对一组矢量值进行滤波,输出结果为该组矢量值滤波后的输出,将此方法称为矢量中值滤波。上述定义下,若想达到最佳滤波效果,地震数据中的同相轴须为水平,但在实际资料处理中难以实现,即便对共中心点(CMP)道集使用正常时差(NMO)校正后的地震数据进行滤波也难以达到理想的滤波效果。基于上述情况,多方向矢量中值滤波(multidirectional vector median filter,MD-VMF)得以发展,可沿不同地层倾角方向进行矢量中值滤波。MD-VMF可用下式表示:

式中,p代指一方向,X_j(p)为沿着p方向截得的一矢量。由上式,可用空间道数j与方向p的函数表示矢量间距离,所求矢量中值即为最小距离。在多方向矢量中值滤波中,多个地震道段被当作矢量元素——被沿着一时窗截取的一组矢量。所得矢量中值随方向p不同而不同,每一个方向p都能求得一矢量中值,在所有方向p中,最终的多方向矢量中值滤波结果取最小距离向量。通常,矢量长度需大于半个子波长度,具体选择可通过测试试验而定。多方向标量中值滤波是在向量长度为1的情况下基于多方向矢量中值滤波定义而来。多方向矢量中值滤波因其特性在较多领域均可被开发。多方向矢量中值滤波中矢量中值求取结果源于基础数据,生成的数据存在有效,在将多方向矢量中值滤波应用于去噪及多震源混合波场分离上也均有较好的效果。

基于上述基本原理,本文采用的多震源混合数据分离流程如下:

1)计算伪分离炮记录P'=Γ^HP_blend,结合实际资料,设置相应的多方向矢量中值滤波参数,如时窗长度、方向范围、方向变化率、地震道数,设置迭代次数N_iter,令i=0,初始分离结果P⁰=0。

2)i=i+1。将属于次级震源的上次分离结果Pⁱ按各自随机时间混合,再进行伪分离,并与P'进行相减得到去噪交叉噪声后的地震信号Rⁱ=Γ^-1P_blend-(Γ^HΓ-I)Pⁱ。

3)将上述地震信号与上次分离结果按照混合度进行加权叠加 Raddi= 1NsRaddi+ 1-1NsPⁱ。

4)将 Raddi由共炮集记录分选到其他非共炮集中,进行多方向矢量中值滤波,得到压制交叉噪声后的有效信号,即第i+1次炮分离结果Pⁱ⁺¹。

5)若i<N_iter,则回到第2步,继续循环,否则退出迭代,分离结束。

为了验证本文方法的有效性与正确性,首先进行理论数据的多震源分离试验。按照海洋地震采集模式模拟出一条二维测线,其混合激发方式和随机激发时间如图2所示,即等间距相同的两个单炮按照一定的激发延迟时间混合激发。该激发方式可看作两条震源船在一条测线上同时采集。常规单震源激发一共进行120炮,每一炮均有120个检波点接收,炮点与检波点重合。按照震源激发延迟时间,将第1炮与第61炮混合成第1混合炮,第2炮与第62炮混合成第2混合炮,以此类推,一共模拟出60个混合炮数据。图3为单震源激发时第50炮、第110炮以及这两炮混合激发的第50个混合炮。

分别采用迭代预测相减法与本文方法对模拟出来的数据进行去噪分离。多方向矢量中值滤波参数选择为:时窗长度选择5个点,地震道选为7道,矢量方向范围设置为(-8,+8)。图4给出了两种方法去噪后的某一共检波道集,其中第0~120道为混合采集时的共检波道集,第120~240道为去噪后的有效数据,240~360道为两者之差,即次震源导致的随机噪声。从图4红圈区域对比来看,本文去噪效果更好。图5给出了两种方法的信噪比变化曲线,可以看出随着迭代次数的增加,本文方法分离出的共检波点道集信噪比更高,信噪比随着迭代次数的增加不断提升并最终稳定,但预测相减法在后期迭代过程中出现明显的震荡现象。

为了进一步验证上述结论,改变多方向矢量中值滤波参数,对两种方法进行另外3组重复试验。表1给出了3组不同实验参数的信噪比,图6为对应的迭代曲线。可以看到当参数不同时,采用本文方法去噪后的信噪比均要高于传统的迭代预测相减法,且信噪比曲线在迭代后期不会出现震荡现象,本文方法更为强健。

图7给出了两种方法得到的第50混合炮的分离结果,该混合炮是由第50单炮与第110炮单炮数据混合而来,将分离结果与图5对比,其共炮道集与未混合前的单炮共炮道集较一致,且本文方法的分离结果比常规迭代预测相减法的分离结果与理论数据更加吻合。

实际数据使用的是2017年石油地球观测系统(PGS)为东京的IEEE信号处理竞赛提供的海洋多震源采集数据^[20]。海上实际多震源数据采集使用二维拖缆采集,共激发256炮,256道接收,数据采样率为4 ms。多震源激发过程中,相邻的左右源混合激发,震源船和拖缆同时沿着测线移动,进行下一个混合震源的激发。由于该实际数据中,单震源激发时间相差较大,难以看出分离效果。本文通过将主震源数据复制后设计一随机激发延迟时间当作次震源激发,将数据混合后作为实际多震源采集数据,如图8所示,取前100道,用本文方法对实际多震源采集数据进行先混合,然后分离。

对其中100个混合炮进行混合炮集的分离计算,这里仅显示每炮近道100道、时窗范围为0~3 000 ms。由于使用的海上采集方式为检波器与空气枪距离固定,并且随着船的移动两者皆向前移动,故将共炮集记录分选到共偏移距道集上来进行去噪。图9a、c为第32左右源对应的共偏移距道集,道集上有明显的交叉噪声,图9b、d为对应的去除交叉噪声后的共偏移距道集。将偏移距道集上交叉噪声压制后,再将其重新分选至共炮点道集,即得到了分离了后的地震记录,实现了多震源分离,图10为第32混合炮分离结果。将第32混合炮分离后的主震源对比海上单震源激发实际数据图8a可以发现,实际数据浅层反射信息能量强,导致深层能量显示较弱,分离出的数据在采集时间1 s处有明显反射波波组。将原始数据与分离后数据绘制出相应局部,采集时间在0.75~1.5 s处,如图11,对应位置有一组明显反射波组,反应出对应地层反射信息,即分离后的数据有效信号明显,可用于做滤波、静校正计算及速度分析等标准叠前处理,分离出的各单炮记录特征明显,原始单炮记录中的有用信号保存较完整,交叉噪声被明显消除。

多震源激发采集作为现阶段及未来地震勘探中采集方法研究的重要发展方向,能极大提高野外采集效率,能否将混合波场分离直接影响这一技术的应用效果。将多震源混合数据经最小二乘求逆得到的伪分离记录分选到共炮道集外的其他道集,通过去噪将次震源产生的交叉噪声去除以实现混合炮的分离。文中采用基于多方向矢量中值滤波的改进迭代去噪法对交叉噪声进行压制,在传统迭代预测去噪法的基础上,通过加权叠加有效信号来提高分离数据信噪比,最终实现多震源混合炮集的分离。从理论数据及海上实际多震源数据处理结果可以看出,本文方法能更好地分离混合数据,不仅在信噪比方面效果更好,在迭代上有更加优越的强健性,不会因参数设置问题出现结果的波动。