面波是陆上地震勘探中的主要干扰波之一, 尤其在我国西部山地地区, 复杂的地质构造使得面波干扰现象尤为突出。面波是一种具有频散特性的规则干扰波, 在炮集上呈线性分布, 其特点是低频率、强振幅, 而且衰减较慢, 严重影响中、深层的有效反射, 降低了地震资料的信噪比。在叠前地震数据处理过程中, 必须有效压制面波干扰, 压制效果直接影响到地震数据的后续处理, 乃至偏移成像效果, 最终影响地震资料的地质解释结果。

为了有效压制面波干扰, 研究人员提出了各种解决方案。李文杰等^[1]提出了基于频率衰减特性的面波压制方法。先在地震记录上确定面波区域, 然后针对不同的炮检距选用不同的低截频率, 在选定的面波区域里进行高通滤波, 有效地衰减了面波区域里的面波能量。但是低截滤波或高通滤波会严重损失深层的低频有效反射波信息。岳龙等^[2]研究了连续小波变换在面波压制中的应用, 根据面波和有效波在小波域能量和分布区域的不同, 提出了利用连续小波变换自动压制面波的方法; 陈文超等^[3]提出了基于连续小波变换的自适应面波压制方法, 根据面波干扰的影响随炮检距变化等特点, 采用连续小波变换实现自适应面波衰减。基于小波变换实现面波压制的有效性, 很大程度上取决于小波函数, 选取不同的小波函数会取得不同的面波压制效果。李杨等^[4]提出了基于广义S变换(GST)的面波压制技术。采用广义S变换将时间域地震信号变换到时频域, 在时频域内确定面波干扰区域, 利用时频滤波方法压制面波, 达到分频分时衰减面波的目的。由于面波的实际频率难以精确标定, 为了彻底去除面波, 在时频域内确定面波干扰区域时往往存在扩大化的现象, 导致在衰减面波的同时有效反射波也会受到损害。传统压制面波的方法还有F-K滤波切除、τ-p域变换^[5-6]等, 这些方法压制面波都有一定的效果, 但也都有其局限性^[7]。譬如F-K滤波要求有规则的空间采样间隔, 因而复杂条件下面波去除效果不佳, 而且F-K滤波是将若干道记录作为一个数据整体进行处理, 在切除面波的同时, 也会将包含在面波记录中的有效低频信息一并去除, 造成对低频有效反射波的损害。τ-p域变换是依据有效波和干扰波视速度符号和大小的不同来压制干扰, 由于面波的速度和频率从浅到深都会发生变化, 经过τ-p域变换后, 面波数据不会是一个点, 因此很难将其完全去除。近年来出现了一些新的采用多种变换联合压制面波的方法。马见青等^[8]提出了采用S变换和TT变换联合压制面波的方法。包乾宗等^[9]提出了采用连续小波变换和脊波变换联合压制面波的方法。孔庆丰^[10]提出了基于Hilbert-Huang变换压制面波的方法。毕云云等^[11]提出了基于离散曲波变换字典和二维局部离散余弦变换字典组合的面波压制方法。这些利用联合变换压制面波的方法取得了较好效果。

广义S变换具有可变形态的窗函数, 可以获得比S变换更好的局部频谱以及更好的时频聚集性。利用这一优势, 对时域内含有面波干扰的地震记录逐道进行广义S变换, 将单道地震记录的一维时间域信号映射成二维时频图。具有低频特征的面波信号, 其能量在时频图上沿时间轴聚集, 呈长条带状分布, 反射波能量则沿频率轴聚集, 并朝频率轴正方向收敛。根据面波与反射波各自能量聚集特征, 对时频图进行分析, 确定出面波信号的时频区域, 针对该区域设计时频滤波器, 将位于该区域内的信号清零, 得到去除面波之后的地震记录, 再与原始地震记录相减, 可将面波从整个地震记录中分离出来^[12-13]。由于面波与有效反射波在低频范围内可能存在重叠, 在用广义S变换分离面波的过程中, 会损害重叠区内低频有效反射波。为减少这种损害, 选择有利于面波压制的小波函数, 将分离得到的面波记录再次进行二维离散小波分解, 得到4个小波系数分量, 对面波主要集中的低频高波数分量进行高通滤波, 滤除面波成分, 然后再进行小波系数分量重构, 重构结果中将不再包含面波成分, 但有效反射波信号得以保留。最后将重构结果叠加到已实现面波分离的地震记录中, 实现在压制面波干扰的同时, 最大限度地减少对低频有效反射波的损害。本文在前人研究的基础上提出了一种新的广义S变换与二维离散小波变换联合压制面波的方法, 并通过理论模型测试和实际地震资料处理验证了本文方法的有效性。

1 基本原理 1.1 广义S变换

S变换是由STOCKWELL等^[14]于1996年提出的一种加时窗傅里叶变换, 是介于短时傅里叶变换和连续小波变换之间的一种非平稳信号分析处理方法。传统S变换采用形态固定的高斯窗函数, 给实际应用带来一些不便。许多学者对时窗函数进行改造, 提出了广义S变换^[15-17]。MANSINHA等^[18]提出了如下可变形态窗函数:

$ g\left( {t,f} \right) = \frac{{\left| f \right|}}{{k\sqrt {2{\rm{ \mathsf{ π} }}} }}{{\rm{e}}^{ - \frac{{{f^2}{t^2}}}{{2{k^2}}}}},k \in {\bf{R}},且\;k > 0 $

(1)

式中:t为时间; f为频率; k为时窗尺度。当k取值较大时, 时窗变宽, 可提高信号的频率分辨率, 当k取值较小时, 时窗变窄, 可提高信号的时间分辨率, 从而得到可调节时间、频率分辨率的广义S变换:

$ {S_{{\rm{GST}}}}\left( {\tau ,f} \right) = \int_{ - \infty }^{ + \infty } {h\left( t \right)\frac{{\left| f \right|}}{{k\sqrt {2{\rm{ \mathsf{ π} }}} }}{{\rm{e}}^{ - \frac{{{f^2}{{\left( {t - \tau } \right)}^2}}}{{2{k^2}}}}}{{\rm{e}}^{ - {\rm{i}}2{\rm{ \mathsf{ π} }}ft}}{\rm{d}}t} $

(2)

式中:S_GST(τ, f)为时域信号h(t)的广义S变换结果; f为频率参数; k为时窗尺度参数; τ为时窗的时间中点。

时域信号h(t)可通过GST反变换进行重构, 反变换公式为:

$ \begin{array}{l} h\left( t \right) = \int_{ - \infty }^{ + \infty } {\left[ {\int_{ - \infty }^{ + \infty } {{S_{{\rm{GST}}}}\left( {\tau ,f} \right){\rm{d}}\tau } } \right]{{\rm{e}}^{{\rm{i}}2{\rm{ \mathsf{ π} }}ft}}{\rm{d}}f} \\ \;\;\;\;\;\;\; = \int_{ - \infty }^{ + \infty } {H\left( f \right){{\rm{e}}^{{\rm{i}}2{\rm{ \mathsf{ π} }}f}}{\rm{d}}f} \end{array} $

(3)

式中:H(f)为h(t)的傅里叶变换。

由于GST与其傅里叶谱保持直接联系, 在进行正、反变换的过程中, 可以实现信号的无损转换, 这也为实际应用带来极大方便。GST能够标定信号中不同频率成分在各个时刻的频谱, 精确表达有效信号和干扰信号在时频域内各自的特征, 因而可以利用GST来进行时频滤波。一般时域信号可表示为有效信号s(t)与干扰信号n(t)之和:

$ h\left( t \right) = s\left( t \right) + n\left( t \right) $

(4)

上式两边同时取GST:

$ {S_{{\rm{GST}}}}\left( h \right) = {S_{{\rm{GST}}}}\left( s \right) + {S_{{\rm{GST}}}}\left( n \right) $

(5)

对由(5)式得到的时频图进行分析, 根据干扰信号的时间和频率范围, 确定其时频区域D, 针对该区域设计时频滤波器:

$ M\left( {t,f} \right) = \left\{ \begin{array}{l} 1,\left( {t,f} \right) \not\subset D\\ 0,\left( {t,f} \right) \in D \end{array} \right. $

(6)

用时频滤波器对(5)式进行滤波处理, 使得S_GST(n)项为0, 保留S_GST(s)项, 经GST反变换就可以得到去除干扰之后的有效时域信号。

1.2 二维离散小波变换

对L²(R)空间满足允许条件的小波母函数Ψ(x)进行伸缩或平移, 可得一维连续小波函数:

$ \begin{array}{l} {\mathit{\Psi }_{a,\tau }}\left( x \right) = \frac{1}{{\sqrt a }}\mathit{\Psi }\left( {\frac{{x - \tau }}{a}} \right)\\ a,\tau \in {\bf{R}},且\;a > 0 \end{array} $

(7)

式中:a为与频率相关的尺度因子; τ为与时间相关的平移因子。如果取a=2^j, τ=2^jk, 其中j, k∈Z, 则可得到一维离散小波函数:

$ {\mathit{\Psi }_{j,k}}\left( x \right) = \frac{1}{{\sqrt {{2^j}} }}\mathit{\Psi }\left( {\frac{x}{{{2^j}}} - k} \right) $

(8)

同理, 可定义一维离散尺度函数:

$ {\mathit{\Phi }_{j,k}}\left( x \right) = \frac{1}{{\sqrt {{2^j}} }}\mathit{\Phi }\left( {\frac{x}{{{2^j}}} - k} \right) $

(9)

将一维小波函数推广, 可定义二维离散小波函数及二维离散尺度函数。若h(x)和g(x)分别为与Φ_{j, k}(x)和Ψ_{j, k}(x)对应的镜像滤波器函数, 可得到二维离散小波分解的Mallat算法公式^[19]:

$ \left\{ \begin{array}{l} {A_{j + 1}}\left( {x,y} \right) = \sum\limits_{m,n} {{A_j}\left( {x,y} \right)h\left( {m - 2x} \right)h\left( {n - 2y} \right)} \\ D_{j + 1}^1\left( {x,y} \right) = \sum\limits_{m,n} {{A_j}\left( {x,y} \right)h\left( {m - 2x} \right)g\left( {n - 2y} \right)} \\ D_{j + 1}^2\left( {x,y} \right) = \sum\limits_{m,n} {{A_j}\left( {x,y} \right)g\left( {m - 2x} \right)h\left( {n - 2y} \right)} \\ D_{j + 1}^3\left( {x,y} \right) = \sum\limits_{m,n} {{A_j}\left( {x,y} \right)g\left( {m - 2x} \right)g\left( {n - 2y} \right)} \end{array} \right. $

(10)

实际地震记录是二维数据, 水平方向表示偏移距, 垂直方向表示时间, 利用Mallat公式进行二维离散小波分解, 可得到如图 1所示的4个小波系数分量:低频低波数分量LL、低频高波数分量LH、高频低波数分量HL和高频高波数分量HH, LL分量还可以继续进行下一层分解。二维离散小波变换对信号采取二倍抽取方式, 随着分解层数的增加, 同一类型分量的频率和波数会成倍数递减, 例如LL3的频率和波数是LL2的一半, 而LL2的频率和波数又是LL1的一半。

1.3 联合压制面波的步骤

利用GST和二维离散小波变换联合压制面波的具体步骤如下。

1) 从单炮地震记录中按顺序依次抽取每一道信号进行广义S变换, 得到该道信号对应的二维时频图。

2) 面波和有效反射波在二维时频图上的能量聚集和图像显示明显不同, 面波能量在低频范围内沿时间轴聚集, 其图像呈长条带状分布, 反射波能量则沿频率轴聚集, 其图像向频率轴正方向收敛, 据此可在时频图上找出面波区域。

3) 针对面波区域设计时频滤波器, 将面波对应区域内的数据充0, 并进行GST反变换, 得到去除面波后的单道时域信号。对整个单炮记录重复步骤1)至步骤3), 即可得到整个记录去除面波后的结果。

4) 用原始记录与去除面波后的结果相减, 得到分离出来的面波记录。

5) 对分离出来的面波记录, 进行二维离散小波分解, 得到四个小波系数分量, 对面波成分主要集中的低频高波数分量LH进行高通滤波, 其它分量不变。如果低频低波数分量LL中还含有面波成分, 则对该分量进行下一层次的二维离散小波分解, 同时对其中的低频高波数分量LH进行高通滤波, 重复进行, 直到低频低波数分量LL中不再含有面波成分。

6) 从最后一层开始逐层进行二维离散小波分量重构, 重构结果中将不再包含面波成分。将重构结果叠加到原始记录去除了面波之后的结果中, 即完成了对面波干扰的联合压制。

上述过程中, 在GST环节要注意面波区域的划分, 它关系到面波分离是否彻底。在二维离散小波变换环节, 要注意选用合适的小波函数, 本文通过分析和反复试验, 结合面波信号特征和小波函数的性质, 从理论分析和实际资料处理两方面考虑, 优选出了Sym5小波函数, 实际应用表明, Sym5小波函数的支撑集、滤波器长度、消失矩等都适合于对面波干扰的压制^[20-21]。

2 理论模型测试

图 2为采用波场模拟方式生成的理论模型合成记录, 单震源中间激发, 二维均匀各向同性弹性波介质, 采用变系数吸收边界条件。理论模型设计了2层反射层, 上层为水平层, 下层为倾斜层, 采样间隔为1ms, 共750个采样点, 70道。面波频率25Hz, 反射波频率30Hz。整个记录不含随机噪声。

图 3为理论模型第6道数据进行GST后得到的时频图, 根据面波与反射波各自在时频图上的能量聚集特点, 很容易区别出来。由于面波频率与反射波频率很接近, 从图 3a可见它们在时频图的低频范围内内存在重叠区域, 因此在进行时频滤波的时候, 对面波区域充0必然会影响到低频反射波, 如图 3b所示。

图 4为对理论模型全部70道数据分别进行GST, 并进行面波区域充0, 将面波信号分离之后得到的结果, 与图 2对比可以发现, 低频反射波的同相轴形状发生了改变, 这正是由于面波区域充0影响了低频反射波的结果, 换句话说, 低频反射波信号在面波分离过程中受到了损害。图 5为理论模型利用GST分离出来的面波记录, 可以看到, 其中除了面波信号被分离了出来之外, 还有一部分低频反射波信号也被分离了出来。

为了减小对低频反射波信号的损害, 将如图 5所示从理论模型中分离出来的面波记录再次进行二维离散小波分解, 第一层分解得到4个小波系数分量如图 6所示, 可见面波成分主要集中在低频低波数分量LL1和低频高波数分量LH1中。对于第一层LL1分量中还含有的面波, 再进行第二层分解, 同理还可进行第三、第四层分解, 每一层分解都会得到LL、LH、HL、HH4个小波系数分量。对各层的LH分量进行高通滤波, 滤除其中的面波成分, 然后将最后一层的4个小波系数分量重构, 作为前一层的LL分量, 依次进行, 直到第一层为止。重构结果如图 7所示, 可见图上已经没有了面波成分, 有效反射波得以保留。将重构结果叠加到理论模型分离面波之后的记录(图 4)中, 得到的最终结果如图 8所示, 可见面波成分得到很好的压制, 同时有效反射波信号的同相轴得到较大改善, 更加接近于原始理论模型合成记录(图 2), 说明本文方法在压制面波干扰的同时起到了保护低频反射波的作用。

3 实际资料处理

图 9为西部某山地原始单炮地震记录, 共184道, 每道采样点数3500, 采样间隔2ms, 图上可见典型的“扫帚状”面波发育, 其下面的有效反射波信号被完全遮盖, 降低了整个记录的信噪比。图 10a为从该地震记录中抽取第61道信号进行GST得到的二维时频图, 图 10b为在时频图上找出面波区域, 并进行数据充0后的结果。对记录中每一道都进行如上相同的处理, 就完成了整个记录中的面波分离。图 11为采用GST从原始记录中分离面波后的结果, 图 12为GST分离出来的面波记录。在利用GST进行面波分离过程中, 为了使面波彻底分离, 往往会扩大数据充0区域, 从而导致低频有效反射波受到损害。为了减少这种损害, 对GST分离出来的面波记录进行二维离散小波分解, 将二维地震数据转换到时间、空间、频率、波数四维域中, 得到4个小波系数分量, 对包含面波成分的低频高波数分量进行高通滤波, 滤除面波, 保留有效反射波, 然后再进行小波系数分量重构, 最后将重构结果叠加到由GST分离面波之后的记录中, 实现在压制面波的同时减少对低频有效反射波的损害。

地震记录中面波信号在时间-空间域的分布范围与有效波不同, 主要集中在低频低波数分量LL和低频高波数分量LH中。对由GST分离出来的面波记录(图 12)进行二维离散小波分解, 得到第一层分量如图 13所示, 可见LL1和LH1分量中存在面波成分。对包含面波成分的LL1分量还可以进行第二层分解。将各层分解得到的LH分量进行高通滤波, 滤除面波成分, 然后从最后一层开始进行小波系数分量重构, 并作为前一层的LL分量, 依次进行, 直到第一层为止。重构结果如图 14所示, 可见明显同相轴, 它们正是在GST分离面波过程中丢失的低频反射波。其中一些同相轴的“蚯蚓化”是由于多层小波变换后的边界效应和累积误差所导致, 选择适当的小波函数及分解层数可以降低累积误差, 重构结果与由GST分离面波之后的结果记录叠加后“蚯蚓化”现象可以得到改善, 但依然存在, 这是本文方法联合压制面波的局限性。

将重构结果叠加到分离面波之后的记录(图 11)中, 实现联合压制面波的最终结果如图 15所示, 与图 11对比可以看到, 图 15中有效反射波的能量得到加强, 说明本文方法在压制面波干扰的同时能够保护低频有效反射波。

作为对比, 采用传统F-K滤波法对图 9所示原始地震记录进行面波压制。图 16为原始单炮地震记录的F-K谱, 可见面波与有效反射波在低频范围内存在重叠区。图中淡蓝色折线以上部分是面波干扰的压制范围, 只要将该范围内的数据充0切除, 就可以实现面波压制, 但是重叠区内面波干扰与有效信号无法完全分开, 在进行F-K滤波切除时必然会损失掉部分有效反射波。图 17为原始地震记录进行F-K滤波后的结果, 由于重叠区的存在, 切除面波的同时也切除了部分有效反射波, 导致原来被面波覆盖的区域内同相轴不够清晰。

为了进一步说明本文方法的优越性, 分别对传统F-K滤波、单纯GST面波分离以及本文提出的GST与二维离散小波变换联合压制面波的结果与原始记录一起作频谱分析, 如图 18所示。

图 18a为采用传统F-K滤波压制面波之后的频谱, 由于面波干扰和有效信号在低频区域内不可能完全分开, 在进行面波压制时, 不可避免会损害重叠部分的有效反射波。

图 18b为采用GST进行面波分离后的频谱, 观察图中蓝色谱线可见, 频率在10Hz以下的面波被压制掉了, 但是, 该范围内有效反射波也一起被压制掉了, 因此单纯采用GST分离面波的方法来压制面波干扰将造成低频区内有效反射波的损失。

图 18c为采用GST+二维离散小波变换联合压制面波后的频谱, 观察图中蓝色谱线可见, 频率在10Hz以下的强振幅面波信号被压制掉了, 但该频率范围内低振幅信号的能量得以保留, 说明本方法在压制面波的同时, 可以减少对低频区内有效反射波的损害。

4 结论

利用广义S变换将地震记录信号从时间域转换到时频域, 采用时频滤波可以方便地实现面波分离, 但分离过程可能损伤低频有效反射波。为了避免或减少这种损伤, 将GST分离出来的面波记录进行二维离散小波分解, 对面波成分主要集中的低频高波数分量进行高通滤波, 滤除面波, 保留有效反射波, 然后再进行小波系数分量重构, 最后将重构结果叠加到已分离面波的地震记录中, 实现在压制面波干扰的同时, 最大限度地减少对低频有效反射波的损害, 整个操作仅对地震记录中的面波区域及面波成分进行分离与压制, 其它域不受影响。通过理论模型测试及实际地震资料处理验证, 本文提出的这种联合压制面波方法切实可行, 面波压制效果明显, 在压制过程中能较好地保护低频有效反射波, 提高了地震记录的信噪比。