0 引言

多震源同时采集方法由Silverman^[1]首先提出，因效率高而广泛应用于实际生产。该方法因采集过程中有多个震源激发，采集的地震数据存在严重的混叠噪声，严重影响地震资料的后续处理。多震源地震数据的分离，目前主要有反演^[2-5]和滤波^[6]两种方法。

滤波分离方法是根据地震信号和混叠噪声在不同域中的分布特性差异分离混叠噪声。如Huo等^[6]在非炮域使用矢量中值滤波法实现混叠噪声分离；Mahdad等^[7]根据反射信号和混叠噪声在不同时域相干性差异，基于干扰噪声估计和自适应去噪的迭代法分离混叠噪声；Doulgeris等^[8]通过调整阈值在共检波点域使用FK滤波分离混叠噪声；韩立国等^[9]联合Curvelet变换和多级中值滤波分离混叠噪声；Gan等^[10]通过估算局部时窗内同相轴斜率及拉平同相轴进行中值滤波实现混叠噪声分离。上述各种滤波分离方法，主要是在共炮检距道集或共检波点道集内分离反射信号与混叠噪声。中值滤波类的分离方法，根据道集内各地震道间反射信号的相干性与混叠噪声分布的随机性，对地震信号和混叠噪声进行分离；FK等频谱滤波类的分离方法依据不同域内反射信号与混叠噪声的分布区域差异，通过切除混叠噪声分布区域的频谱实现信号与混叠噪声的分离。中值滤波类分离方法，存在破坏信号波组特征、分离不干净等问题；频谱滤波类分离方法会存在混叠噪声频谱与信号频谱在变换域内分布范围重叠的问题，造成混叠噪声无法完全去除或反射信号被切除，还因方法的局限性造成混波、假频等现象。

Stockwell等^[11]提出的S变换方法是一种介于短时傅里叶变换和小波变换之间的非平稳信号分析方法。相比小波变换、短时傅里叶变换等时频分析方法，S变换时频谱具有多尺度聚焦性，且与傅里叶谱直接相关。NMO-CMP道集内因各道相同时刻地震信号代表同一反射点信息，故它们的振幅、相位基本一致，具有较高的相干性^[12-13]；混叠噪声在NMO-CMP道集中随机噪声分布。根据叠加统计原理，叠加道是由各道有效信号叠加生成，各道中噪声的叠加结果趋于零，故叠加道的振幅、相位与各道间也具有很好的相干性。因此，对NMO-CMP道集内各地震道和叠加道做S变换，将其变换到时频域，以叠加道S谱为参考，对道集内各道S谱进行自适应滤波，可有效去除数据中的混叠噪声。

本文以NMO-CMP道集的叠加道S谱为参考设计自适应滤波器，对道集内地震道进行多级滤波，并以多级滤波的结果重新建立叠加道，对地震数据进行迭代滤波，以提取数据中的有效反射信号、分离混叠噪声。本文方法是依据叠加道与地震道S谱间的相干性提取数据中的反射信号，实现信号与混叠噪声的分离，可有效避免假频、混波等现象。理论数据和实际数据模拟的多震源混叠数据试算，验证了本文方法的有效性。

1 方法原理 1.1 S变换

连续信号x(t)的S变换S(τ，f)定义^[11]为

$ \begin{array}{*{20}{l}} {S(\tau ,f) = \int_{ - \infty }^\infty x (t)\frac{{|f|}}{{\sqrt {2\pi } }} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{exp}}\left[ { - \frac{{{{(t - \tau )}^2}{f^2}}}{2}} \right]{\rm{exp}}( - {\rm{i}}2\pi ft){\rm{d}}t} \end{array} $

(1)

式中τ为窗函数的中心点，控制高斯窗在时间轴上的位置。由上式可知，S变换的高斯窗函数和基本小波分别为

$ g(t) = \frac{{|f|}}{{\sqrt {2\pi } }}{\rm{exp}}\left( { - \frac{{{t^2}{f^2}}}{2}t} \right) $

(2)

$ \begin{array}{*{20}{l}} {w(t) = \frac{{|f|}}{{\sqrt {2\pi } }}{\rm{exp}}\left( { - \frac{{{t^2}{f^2}}}{2} - {\rm{i}}2\pi ft} \right)}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = g(t){\rm{exp}}({\rm{i}}2\pi ft)} \end{array} $

(3)

反S变换表达式为

$ x(t) = \int_{ - \infty }^{ + \infty } {\left[ {\int_{ - \infty }^{ + \infty } S (\tau ,f){\rm{d}}\tau } \right]} \exp ({\rm{i}}2\pi ft){\rm{d}}f $

(4)

1.2 NMO-CMP道集中信号和噪声的分布特征

地震信号x(t)由信号s(t)和随机噪声n(t)两部分组成

$ x(t) = s(t) + n(t) $

(5)

多震源数据的地震信号x(t)由三部分组成，包括其本身的信号s(t)、随机噪声n(t)和混叠噪声x_b(t)，即

$ x(t) = s(t) + n(t) + {x_{\rm{b}}}(t) $

(6)

x_b(t)在CMP道集中随机分布，可看作随机噪声(图 1)。根据叠加统计原理，噪声n(t)和混叠噪声x_b(t)叠加后趋于0，则NMO-CMP道集的叠加道为

$ \begin{array}{*{20}{l}} {{x_{\rm{s}}}(t) = \frac{1}{M}\sum\limits_{i = 1}^M {{s_i}} (t) + \frac{1}{M}\sum\limits_{i = 1}^M {{n_i}} (t) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{M}\sum\limits_{i = 1}^M {{x_{{\rm{b}}i}}} (t) \approx \frac{1}{M}\sum\limits_{i = 1}^M {{s_i}} (t)} \end{array} $

(7)

式中M为CMP道数。

考虑到异常振幅噪声和混叠噪声对叠加道的影响，叠加时应根据数据信噪比情况，剔除同时刻各采样点的振幅较大和较小部分，本文取中间三分之一部分生成叠加道，即

$ {x_{\rm{s}}}(t) = \frac{3}{M}\sum\limits_{i = M/3}^{2M/3} {s_i^\prime } (t) $

(8)

式中s′_i(t)为重新排序后的各道振幅。

图 1为合成的多震源理论数据NMO-CMP道集，在1000、1700、2400、2900ms处有四个反射层，混叠噪声在数据中呈随机分布，数据中呈现双曲线形态分布的为转换波。图 2a为图 1数据的叠加道，图 2b~图 2d分别为叠加道、图 1数据第34道、第41道的S谱。图 2b的S谱能量分布与道集和叠加道中的反射层分布一致，图 2c、图 2d的S变换谱中存在严重的混叠噪声干扰。三道的S变换谱对比证明，用式(8)建立叠加参考道可以有效抑制异常振幅噪声和混叠噪声对参考道的影响。

1.3 滤波器设计

在不考虑各道信号的各向异性属性情况下，NMO-CMP道集中各地震道与叠加道S谱的振幅和相位差异可代表其受噪声污染程度，在地震道中消除这种振幅差即可去除地震道中的噪声；在考虑各向异性属性情况下，地震道之间以及地震道与叠加道之间的S谱存在振幅和相位差异，滤波结果要保留这种各向异性属性^[14-15]。因此在设计滤波器时，既要考虑地震道S谱与叠加道S谱之间的一致性，消除数据中的噪声；也需考虑地震道S谱与叠加道S谱之间的非一致性，保留地震数据中的各向异性信息^[16]。所以滤波器设计需要满足以下几点：

(1) 当地震道S谱与叠加道S谱振幅差异较小时，滤波器对地震道的S谱的振幅调整幅度要小，以保证地震道的各项异性差异；反之则滤波器对地震道的S谱的振幅调整幅度要大，以消除存在于S谱中的噪声^[17-19]。

(2) 对地震数据滤波时，滤波器以叠加道S谱为参考对地震道的振幅和相位进行自适应滤波，会对地震道的振幅和相位进行归一化调整，在一定程度上会破坏地震信号的各向异性属性，使部分各向异性信息被当作噪声去除，可通过多级滤波的方法从去除的噪声中提取这些各向异性信息，尽可能地保留地震信号中的各向异性信息。

(3) 地震数据中的各种噪声会影响叠加道的准确性，以受污染的叠加道为参考道进行滤波也会影响滤波结果的准确性，因此以多级滤波结果的叠加道作为参考道，对地震数据重新进行多级滤波(参考道的更新迭代)，可提高滤波结果的准确性。

本文根据滤波器设计需求，参考非局部均值滤波方法^[20]，设计了一个自适应滤波器。首先计算各地震道与参考道S谱之间的相对高斯距离

$ {d_{i,j,k}} = \frac{{|{S_{i,j,k}} - {S_{{\rm{c}},j,k}}|}}{{{S_{{\rm{c}},j,k}}}} $

(9)

式中：S_i，j，k为CMP道集中第i道的S变换结果；S_c，j，k为参考道的S变换结果；j为时间样点序号；k为频率样点序号。

利用相对高斯距离设计滤波器的自适应加权函数

$ {R_{i,j,k}} = \frac{{{d_{i,j,k}}}}{{\sqrt {1 + {{\left( {\frac{{{d_{i,j,k}}}}{h}} \right)}^2}} }} $

(10)

式中h为收敛系数(h≥1)，用来控制加权函数的收敛速度。h越小，加权函数的收敛速度越快，反之则越慢。

令

$ {D_{i,j,k}} = {S_{i,j,k}} - {S_{{\rm{c}},j,k}} $

(11)

设计的滤波器为

$ \left\{ {\begin{array}{*{20}{l}} { {\rm{Re}} ({{\hat S}_{i,j,k}}) = {\rm{Re}} ({S_{i,j,k}}) - [ {\rm{Re}} ({D_{i,j,k}}) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {R_{i,j,k}} Re ({D_{i,j,k}})]}\\ { {\rm{Im}} ({{\hat S}_{i,j,k}}) = {\rm{Im}} ({S_{i,j,k}}) - [ {\rm{Im}} ({D_{i,j,k}}) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {R_{i,j,k}} {\rm{Im}} ({D_{i,j,k}})]} \end{array}} \right. $

(12)

式中$\hat{S}_{i, j, k}$为滤波结果。对滤波结果$\hat{S}_{i, j, k}$做反S变换即可得到分离噪声后的地震数据。

图 3为本文所设计滤波器的响应曲线，其中收敛系数h=1。可以看出：在地震道S谱与参考道S谱的高斯距离较小时，滤波结果的振幅调整就较小；当地震道S谱与参考道S谱的高斯距离较大时，滤波结果的振幅调整就较大，并稳定收敛于相对高斯距离1附近；最大收敛高斯距离可根据数据各向异性情况，通过选择h值的大小调节。

表 1对比了使用本文方法、指数函数^[20]、对数函数^[20]设计滤波器的加权函数响应。应用对数函数(d/ln[(abs(d)])设计的滤波器，在相对高斯距离±1和0处，滤波器存在间断点；在相对高斯距离-2处，存在噪声放大(滤波结果的相对高斯距离变大)；当噪声严重时(相对高斯距离大)，滤波结果收敛慢。应用指数函数设计的滤波器(d/exp[abs(d)])，滤波结果收敛速度过快，存在过度收敛问题；当相对高斯距离增大到一定程度，会产生溢出。相比以上两种滤波器，本文设计的滤波器满足了上文所提的要求，且不存在间断点、噪声放大和过度收敛等问题。

1.4 处理流程

本文算法流程(图 4)如下。

(1) 对做过预处理的数据做速度分析和NMO。

(2) 对输入数据做多级滤波。多级滤波时，一级滤波以输入信号的叠加为参考道对多震源地震数据做滤波，其余各级滤波都是将上级滤波输出的噪声叠加，以噪声的叠加道为参考道对上级滤波所得噪声进行滤波提取信号。每一级滤波结束后，将该级滤波所得信号与之前各级滤波所得信号进行叠加，作为该级滤波的信号，对比该级滤波信号与上一级滤波信号，如无明显变化，则进入迭代滤波；若有明显变化，则进入下一级滤波。

(3) 迭代滤波。以上次多级滤波得到的信号做叠加获得参考道，重新对多震源地震数据做多级滤波，直至上下两级滤波结果无变化，本次迭代滤波结束。若迭代滤波结果与迭代前信号相比无变化，则以本次迭代滤波结果为最终滤波信号输出；反之则以本次迭代滤波所得信号做叠加，以该叠加道为参考重新对多震源地震数据做多级滤波。

因本文方法需做叠加获取参考道，故应用前需注意以下事项：

(1) 滤波前需对数据中的有效数据做振幅补偿、去除规则噪声和初至切除，无需针对混叠噪声做去噪处理。这样可提高参考道的信噪比和准确性。

(2) NMO所用速度场需要精准，以保证迭代过程中一直使用。本文方法的理论前提是各道NMO后与叠加道的S谱基本一致，在滤波时对各道的相位和振幅调整，以达到提取有效反射信号、去除噪声的目的，速度是影响本文方法滤波效果的关键因素。

2 理论数据验证

应用理论合成的多震源数据验证本文的多源地震数据分离方法。

图 5a和图 5b是理论单炮记录，图 5c为其合成的多震源数据，其中图 5b为有效数据，图 5a为混叠噪声。图 6为混叠噪声分离前、后炮集记录，可以看出，经两次迭代滤波后，合成数据混叠噪声和有效反射信号得到完全分离，并且原始数据中的三组转换波也得到有效分离。多次波与转换波的分布形态相似，因此本方法亦适用于多次波压制。图 6中第一反射同相轴在混叠噪声分离后变短是由动校拉伸畸变切除引起。本文方法需做NMO和叠加，浅层远炮检距部分动校拉伸畸变严重，会影响参考道的准确性，故噪声分离前需对地震数据做初至切除和动校拉伸切除。

多震源混叠数据NMO后的CMP道集如图 7a所示，其中1000、1700、2400、2900ms处水平同相轴为有效反射波。应用本文设计的滤波器进行滤波，其一级、四级分离后的数据如图 7b、图 7c所示，两次迭代滤波后数据如图 7d所示，图 7e为最终去除的混叠噪声。可见，经两次迭代滤波后，有效分离了混叠噪声和有效反射信号。

图 8为图 7数据混叠噪声分离前、后S谱对比，可见，经过本文方法滤波后，地震道S变换谱中的信号和混叠噪声已得到有效分离。

图 9为混叠噪声分离前、后叠加剖面对比，可见，经本文方法分离混叠噪声后，剖面中的混叠噪声得到了较好分离，与混叠前叠加剖面相比，整体特征一致，各反射层的细节特征都得以保留。

由以上分析可知，本文方法可有效分离理论合成的多震源采集数据中的混叠噪声，恢复有效反射信号。

3 实际地震数据模拟验证

为验证本文方法分离实际多震源数据混叠噪声能力，选取一段传统方法采集的二维地震数据，将其人工合成为多震源地震数据并应用本文方法进行混叠噪声分离。

考虑到一致性、规则噪声和初至等对参考道的影响，在分离前需对多震源数据中的有效数据做振幅补偿、初至切除、动校拉伸切除和适当规则噪声(如面波等)压制，以保证参考道的准确性。

图 10为实际数据合成的多震源采集数据。图 11为混叠噪声分离前、后的CMP道集。图 11a道集为分离前数据，其中类似于面波的噪声为混叠道集中的面波。一级滤波后，原始数据(图 11a)中的混叠噪声得到了相当程度的分离，但肉眼仍可见残余(图 11b)；四级滤波后，一级滤波结果中的残余混叠噪声得到进一步分离，但仍有一定的残余(图 11c)；两次迭代滤波后，数据中的混叠噪声得到有效分离(图 11d)。分离出的噪声中未见有效信号(图 11e)，说明多震源混叠数据分离彻底。

图 12为混叠噪声分离前、后炮集记录。合成数据中混叠噪声明显(图 12a)；一级滤波后，混叠噪声得到了一定程度的分离，但仍残存有较多的混叠噪声(图 12b)；四级滤波后，数据中的混叠噪声得到了进一步的分离，但仍有一定的残余(图 12c)；二次迭代后，混叠噪声得到了有效分离，与合成前有效单炮数据(图 12e)的波组特征一致，并且信噪比得到了一定程度的提高(图 12d)；从分离出的混叠噪声(图 12f)未见有效信号，说明混叠噪声和信号得到了有效分离。

图 13为混叠噪声分离前、后叠加剖面对比。与混叠前叠加剖面(图 13a)相比，混叠后叠加剖面(图 13b)存在严重的混叠噪声；经两次迭代分离混叠噪声后，剖面中的混叠噪声得到了有效分离，信噪比也得到了一定程度的提高(图 13c)；由分离出的混叠噪声叠加剖面(图 13d)可见，数据中的信号和混叠噪声得到了有效分离。

由以上试算、对比可知，本文方法可有效提取多震源数据中的有效反射、分离多震源数据中的混叠噪声。

4 结论

理论和实际地震数据合成的多震源采集数据试验结果表明，本文根据NMO-CMP道集地震数据与叠加参考道S谱之间的差异设计的自适应滤波器，可有效提取多震源数据中的有效反射信号。以各级滤波中的叠加道为参考，在地震数据和去除噪声中提取有效信号，可保留有效反射信号的各向异性特征；通过对多级滤波得到的结果进行迭代，可以有效避免各种噪声在建立参考道时对参考道的污染。

本文方法以叠加道为参考，对有效反射信号进行信噪分离，因此应用前需对多震源数据中的有效数据做初至切除、振幅补偿和规则噪声适当去噪。根据噪声分离过程中转换波的去除效果，本方法还有较好的多次波压制能力。