随机噪声具有强烈的随机性, 在地震资料中普遍存在^[1], 因此无论是天然地震还是人工地震, 随机噪声的处理都是提高地震资料信噪比的关键之一。常用的地震信号随机噪声压制方法有拉东变换、经验模态分解、曲波变换与小波变换等^[2], 其中小波变换具有多分辨率分析的特点, 能够同时表征信号在时域和频域的局部特性^[3]。经小波分解得到小波系数后, 根据有效信号与随机噪声小波系数的不同而设定阈值, 对随机噪声的小波系数进行处理后再对小波系数进行重构, 可达到信号去噪的目的。

提升小波算法被称为“第二代小波”, 是一种基于时域运算的信号分析方法, 可克服经典小波降噪方法局部信息丢失的缺点^[4]。对于小波阈值去噪, 阈值函数与阈值的选取对地震资料去噪的结果具有很大影响。经典的阈值函数包含硬阈值和软阈值, 其中硬阈值函数是当小波系数的绝对值大于设定的阈值时系数保持不变、小于阈值时系数为零, 这会导致其在阈值处不连续, 从而使重构信号产生一定的震荡; 而软阈值函数可克服硬阈值函数在阈值处不连续的缺点, 当小波系数的绝对值大于阈值时减去阈值, 小于阈值时则进行置零处理, 但经过软阈值处理的小波系数与实际小波系数存在恒定偏差, 从而影响重构信号的精度。

本文在经典软、硬阈值函数的基础上提出一种改进的阈值函数, 从理论上证明该阈值函数可克服经典阈值函数的不连续与存在恒定偏差的缺点, 并通过信号仿真证实经过该阈值函数处理后可很好地压制高斯白噪声, 提高含噪信号的信噪比。最后将改进阈值函数的提升小波算法运用于气枪震源, 分析其在人工震源去噪方面的可行性。

1 改进阈值函数的提升小波变换 1.1 提升小波

Mallat算法利用多分辨率分析分解信号, 从而去除信号间的相关性, 而提升小波算法则是将信号分解为偶数采样点和奇数采样点^[5]。由于相邻数据的相关性较强, 可通过奇数采样点估计偶数采样点, 再去除包含在偶数采样点中关于奇数采样点的信息, 从而消除信号间的相关性。提升小波算法的实现包含分裂(spilit)、预测(predict)和更新(update)。其中, 分裂是将原始信号S_i按照序列的奇偶性分裂为2个互不相关的子集, 即偶数序列S_e与奇数序列S_o。预测是利用原始信号相邻数据间的相关性, 用偶数序列S_e与1个预测算子P来预测奇数序列S_o, 得到P(S_o), 并将奇数序列S_o与奇数序列预测值P(S_o)的差值d_i-1作为小波系数, 其中d_i-1=S_o-P(S_o)。更新则是为使分裂后得到的偶数序列S_e保持原始信号的某些特性, 构造一个更新算子U(d_i-1)对其进行更新, 得到尺度系数a_i-1^[6]：

$ a_{i-1}=S_{e}+U\left(d_{i-1}\right) $

(1)

1.2 阈值去噪

对含噪地震信号进行一维离散小波分解, 得到一系列小波系数。由于有效信号在时间域或空间域具有一定的连续性, 而噪声信号则表现为随机性, 从而导致小波域中有效信号对应的小波系数值较大, 而噪声信号对应的小波系数值较小。根据该特点设置一个阈值, 将小于该阈值的小波系数去除, 即将噪声主导对应的小波系数去除, 再对剩余的小波系数进行重构, 达到去除噪声的目的, 因此阈值与阈值函数的选择具有重要意义：

$ \lambda=\sigma \sqrt{2 \log N} $

(2)

$ \sigma=\frac{\operatorname{median}\left(c d_{n}\right)}{0.674\;5} $

(3)

研究表明^{[3, 6-7]}, 随着分解层数的增加, 各层小波细节系数中噪声含量会逐渐减小, 因此阈值应当有所变化。经过一系列仿真实验, 本文使用一种改进的阈值, 具体公式为：

$ \lambda=\frac{\sigma \sqrt{2 \log N}}{j} $

(4)

式(2)~式(4)中, λ为去噪阈值, σ为噪声强度, N为信号长度, cd_n为小波分解后第n层高频小波系数, j为分解层数。

在得到阈值后, 可使用硬阈值函数和软阈值函数对小波系数进行处理。硬阈值函数表达式为：

$ {\hat w_{j, k}} = \left\{ \begin{array}{l} \;\;0, \left| {{w_{j, k}}} \right| < \lambda \\ {w_{j, k}}, \left| {{w_{j, k}}} \right| \ge \lambda \end{array} \right. $

(5)

图 1为硬阈值函数图像, 从图中可以看出, 硬阈值函数在阈值±λ处不连续。

软阈值函数表达式为：

$ \hat{w}_{j, k}=\left\{\begin{array}{l} 0, \left|w_{j, k}\right|<\lambda \\ \operatorname{sign}\left(w_{j, k}\right)\left|w_{j, k}-T\right|, \left|w_{j, k}\right| \geqslant \lambda \end{array}\right. $

(6)

图 2为软阈值函数图像, 从图中可以看出, 软阈值函数可弥补硬阈值函数的缺点, 在阈值±λ处连续, 但存在恒定偏差问题。式(5)~式(6)中, sign为符号函数, w_{j, k}为小波系数, $\hat{w}$_{j, k}为小波阈值处理后的小波系数。

1.3 改进的阈值函数

为弥补传统阈值函数存在的缺点, 本文提出一种改进的阈值函数, 具体表达式为：

$ f\left(w_{j, k}\right)=\hat{w}_{j, k}=\left\{\begin{array}{l} 0, \left|w_{j, k}\right|<\lambda \\ \frac{2}{\pi} \arctan \left(\left|w_{j, k}\right|-\lambda\right) w_{j, k}+\left(1-\frac{2}{\pi} \arctan \left(\left|w_{j, k}\right|-\lambda\right)\right) \operatorname{sign}\left(w_{j, k}\right) \\ \left(\left|w_{j, k}\right|-\frac{\lambda}{\log _{2}\left(\frac{\left(1+\left|w_{j, k}\right|\right)^{2}-\left(2+\left|w_{j, k}\right|\right) \lambda}{1+\left|w_{j, k}\right|-\lambda}\right)+1}\right), \left|w_{j, k}\right| \geqslant \lambda \end{array}\right. $

(7)

本文将从数学方面来分析改进的阈值函数在±λ处的连续性、改进阈值函数的渐近线及克服软阈值函数的恒定偏差。

当w_{j, k}→λ⁺时, $\lim\limits_{w_{j, k} \rightarrow \lambda^{+}} \hat{w}_{j, k}=\frac{2}{\pi} \times 0 \times w_{j, k}$+$1 \times 1 \times\left(\left|w_{j, k}\right|-\frac{\lambda}{\log _{2}\left(\frac{\left(1+\left|w_{j, k}\right|\right)^{2}-\left(2+\left|w_{j, k}\right|\right) \lambda}{1+\left|w_{j, k}\right|-\lambda}\right)+1}\right)$=$\left|w_{j, k}\right|-\lambda=0$; 当w_{j, k}→λ^-时, $\sum\limits_{w_{j, k} \rightarrow \lambda^{-}} \hat{w}_{j, k}$=0, 即$\sum\limits_{w_{j, k} \rightarrow \lambda} \hat{w}_{j, k}=0$。同理, 当w_{j, k}→-λ时, $\sum\limits_{w_{j, k} \rightarrow-\lambda^{+}} \hat{w}_{j, k}=\sum\limits_{w_{j, k} \rightarrow-\lambda^{-}} \hat{w}_{j, k}=0$。综上所述, 改进的阈值函数在±λ处连续, 从而可克服硬阈值函数不连续的缺点, 使去噪后的信号更为平滑, 有利于消除伪吉布斯现象。

当w_{j, k}→+∞时, $\sum\limits_{w_{j, k} \rightarrow+\infty} \frac{\hat{w}_{j, k}}{w_{j, k}}$=$\frac{\frac{2}{\pi} \times \frac{\pi}{2} \times w_{j, k}+0}{w_{j, k}}=1$; 当w_{j, k}→-∞时, $\sum\limits_{w_{j, k}} \rightarrow-\infty \hat{w}_{j, k}=\frac{\frac{2}{\pi} \times \frac{\pi}{2} \times w_{j, k}+0}{w_{j, k}}=1$。综上所述, 当w_{j, k}→∞时, f(w_{j, k})逐渐趋近w_{j, k}, 从而可消除软阈值函数f(w_{j, k})与w_{j, k}之间的恒定偏差。

从图 3可以直观看出, 改进的阈值函数可弥补硬阈值函数在阈值±λ处不连续的缺点, 并且在阈值之后能快速趋近硬阈值函数, 从而弥补了软阈值函数存在恒定偏差的缺点。

2 信号仿真

本文选用MATLAB中4种标准测试信号(Blocks, Bumps, Doppler, Heavy sine)分别加入标准偏差为8的高斯白噪声, 再使用本文改进的阈值函数及传统的软、硬阈值函数进行去噪。

从图 4~7可以看出, 软阈值函数与硬阈值函数去噪均会使原始信号严重失真, 甚至变形。相比于软阈值函数, 硬阈值函数去噪结果存在较多“毛刺”凸起, 这是由硬阈值函数在λ处不连续造成的。相比于硬阈值函数, 虽然软阈值函数可弥补其不连续的缺点, 使得去噪后的信号较为光滑, 但软阈值函数在小波系数大于λ时减去λ, 存在恒定偏差, 从而造成去噪后的结果过于平滑, 影响去噪效果。另外, 改进的阈值函数能在很好地去除高斯白噪声的同时, 改进硬阈值函数的“毛刺”现象及软阈值函数的过度光滑, 使信号不失真。

本文将仿真信号的信噪比(SNR)、均方根误差(RMSE)及能量比例(P_er)作为去噪后的评价指标：

$ \mathrm{SNR}=10 \lg \frac{\sum\limits_{i=1}^{n} f^{2}(t)}{\sum\limits_{i=1}^{n}[f(t)-\hat{f}(t)]^{2}} $

(8)

$ \mathrm{RMSE}=\left[\frac{\sum\limits_{i=1}^{n}[f(t)-\hat{f}(t)]^{2}}{n}\right]^{\frac{1}{2}} $

(9)

$ P_{e r}=\frac{\left[\frac{1}{n} \sum\limits_{i=1}^{n} \hat{f}(t)^{2}\right]^{1 / 2}}{\left[\frac{1}{n} \sum\limits_{i=1}^{n} f(t)^{2}\right]^{1 / 2}} $

(10)

式(8)~式(10)中, f(t)为原始信号, $\hat{f}(t)$为去噪后信号, n为信号长度1 000。信噪比(SNR)越大、均方根误差(RMSE)越小, 表明去噪效果越好, 有效信息保留得越多。

在4种信号中加入标准偏差为8的高斯白噪声, 利用控制变量法, 选择提升小波基为bior6.8, 分解与重建尺度为5, 分别使用改进的阈值函数与传统的阈值函数进行去噪处理, 评价指标结果见表 1。通过SNR与RMSE可以看出, 相比于传统的阈值函数, 改进的阈值函数可大大提高去噪后信号的信噪比, 有效去除高斯白噪声, 同时保证仿真信号不失真, 去噪效果明显优于传统方案。从去噪前后的P_er可以看出, 改进的阈值函数的提升小波在4种仿真信号中均保持较高比例, 说明相对于软阈值函数与硬阈值函数, 改进的阈值函数可保留较多的原始信号能量成分, 而软阈值函数与硬阈值函数在去噪过程中不仅会去除噪声, 同时也会去除较多的有效信号。

3 实际气枪震源去噪 3.1 改进阈值函数与传统阈值函数实际去噪对比

地震波是研究地球结构的有效工具之一^[8], 天然地震学与勘探地震学都面临着提高地震资料信噪比的难题^[9]。气枪震源是应用于地球内部结构探测的一种新型人工震源^[10], 相比于其他人工震源具有可控性好、低频成分丰富、信号重复性高等优点。但气枪震源信号较弱, 受外界噪声干扰较大, 且传播距离有限, 衰减与噪声干扰会使远处台站接收到的气枪信号信噪比极低, 难以分辨。

提高气枪信号信噪比的常用方法是利用气枪震源的重复性进行信号叠加^[11], 并在此基础上进行带通滤波处理。由于气枪信号的有效信号集中在2~5 Hz之间, 通常使用2~6 Hz的带通滤波处理气枪信号^[9], 但该方法存在一定问题, 即不能有效滤除与有效信号具有相同频率的噪声。本文使用云南宾川气枪信号发射台2016年接收到的实验激发数据, 信号接收台站是震中距为151.13 km的53252台站和震中距为115.89 km的53253台站。在叠加气枪信号和进行2~6 Hz带通滤波的基础上对改进阈值函数的提升小波进行去噪处理(图 8~9), 结果表明, 改进的阈值函数不仅可克服传统软、硬阈值函数不连续和存在恒定偏差的缺点, 还可极大地提高气枪震源地震资料的信噪比。

3.2 实际气枪震源信号评价指标分析

为直观对比改进的阈值函数的去噪信号与原始信号, 将2个信号进行叠加, 结果见图 10~11。结合图 8~11可以明显看出, 波软阈值函数与硬阈值函数的提升小波在去噪过程中会滤除部分气枪信号的相位信息, 导致有效信号严重失真; 而改进阈值函数的提升小波阈值去噪则在压制噪声的同时几乎未造成有效信息的丢失。上述分析表明, 改进的阈值函数在实际的人工地震信号去噪中比传统的阈值函数表现更优。

从表 2可以看出, 相比于软、硬阈值函数, 改进阈值函数的提升小波在实际的气枪信号去噪中信噪比有大幅提高, 均方根误差减小且能量比例接近1。传统的阈值函数去噪能量比例较小, 表明在去噪过程中传统阈值函数的提升小波会去除许多有效信号, 这与图 8~9的结果相符。分析结果表明, 在实际的气枪震源信号处理中, 改进阈值函数的提升小波算法可克服软阈值函数存在恒定偏差与硬阈值函数不连续的缺点, 提升远处台站气枪震源信号的信噪比是一种可用于气枪震源信号数据处理的有效方法。

3.3 改进阈值函数去噪与带通滤波对比

图 8~9均为在进行900次气枪信号叠加的基础上通过2~6 Hz带通滤波进行提升小波的阈值去噪处理, 为比较提升小波阈值去噪与带通滤波去噪的效果, 对53253台站接收到的信号进行500次叠加, 并在此基础上分别进行2~6 Hz与2~8 Hz的带通滤波。

从图 12可以看出, 气枪信号在进行叠加后, 改进的阈值函数的提升小波可较好地压制气枪信号的随机噪声, 且未造成气枪信号失真; 而带通滤波虽然可以滤除随机噪声, 但同时也会滤除部分有效信号, 如55~60 s信号, 从而造成一定程度上的信号失真。

4 结语

本文在传统软、硬阈值函数基础上提出一种新的阈值函数, 并从理论上证明其可在一定程度上弥补软阈值函数存在恒定偏差及硬阈值函数不连续的缺点, 同时将其运用到Blocks、Doppler、Heavy sine和Bumps四种加入高斯白噪声的仿真信号中进行去噪实验。从评价指标SNR和RMSE可以看出, 改进的阈值函数去噪效果优于传统阈值函数, 并且仿真信号不失真。将改进的阈值函数运用到气枪人工震源中发现, 相比于传统阈值函数, 改进的阈值函数能在压制随机噪声并提高远处台站信噪比的同时, 保证有效气枪信号的相位和幅值不失真。将气枪信号进行叠加后发现, 改进的阈值函数在保留有效信号方面优于带通滤波。

致谢: 云南省地震局“主动源创新团队”提供气枪震源数据, 在此表示感谢。