• 首页关于本刊投稿须知期刊订阅编委会期刊合作查询检索English
一种小波综合阈值地震信号降噪方法
一种小波综合阈值地震信号降噪方法
杨敬松*, 王晗, 姚振静, 韩智明, 李亚南, 高强
(防灾科技学院,振动精密测量技术与仪器重点实验室,河北三河 065201)
 [收稿日期]: 2017-10-16
摘要

小波阈值方法中硬、软阈值方法是地震信号降噪常用方法,但容易造成信号中高频信息丢失导致地震误判和漏判情况发生。小波综合阈值方法继承和发展了硬、软阈值降噪方法的优点,对信号高频部分用硬阈值方法,以提高高频信号能量,对信号低频部分用软阈值方法,提高信号降噪能力的同时保证信号连续性和光滑性。利用噪声信号小波系数小和地震信号小波系数大的特征,进行雷克子波降噪仿真实验和实际地震信号降噪实验。仿真实验表明,小波综合阈值方法降噪后波形MSE值最小,且降噪后与原信号波形最近似,降噪后波形高频部分能量增强且抑制低频部分能量。最后,对实际采集的地震信号进行降噪处理,处理后信号中能量增强被压制,利用处理后的信号可得到地震的初至时间。



引言

地震信号采集过程中受外界环境干扰,以及传感器自身不确定性影响,信号夹杂次生噪声、环境噪声和仪器噪声,导致地震误判和漏判情况的发生(郑作亚等,2007范涛,2014)。3种噪声中,仪器噪声对地震信号影响较弱,因此主要通过消除次生噪声和环境噪声带来的干扰,降低误判、漏判现象(李英等,2006)。

单纯的傅里叶方法很难从复杂的噪声环境中分离地震信号,而小波阈值方法可以在时频域表征信号变化,与傅里叶变换、窗口傅里叶变换相比,具有细节区分能力(孔祥茜等,2005刘霞等,2010)。典型的小波阈值方法有硬阈值与软阈值降噪方法,硬阈值处理后容易造成信号的不连续,导致有效信号丢失;软阈值处理后信号与原信号相差较大,影响信号重构效果(侯跃伟等,2015)。Gao和Bruce提出半软阈值方法对硬、软阈值方法进行了改进(魏学强等,2016唐守峰等,2011),对阈值函数进行加权平均,将加权因子设为0.5,但仍不能解决信号连续性的问题,且小波系数估计值与真实值偏差较大。新发展起来的基于S变换的软阈值降噪方法(曲中党等,2015)在S变换的基础上结合软阈值方法提高地震信号信噪比,有效提高地震信号降噪水平。小波综合阈值方法继承和发展了硬、软阈值的优点,结合软、硬阈值函数的优势对阈值函数进行改进,对信号的小波系数高频部分用硬阈值方法提高高频信号能量,对小波系数低频部分用软阈值方法保持信号连续性(Xia等,2017曾宪伟等,2010),能够在保留信号连续性的同时提高高频信号能量。

针对小波阈值降噪中存在软、硬阈值函数不能有效消除噪声信号对地震信号的影响等问题,提出小波综合阈值方法对阈值函数进行改进。改进后,小波综合阈值函数的小波系数与真实函数的小波系数无限接近,既保持信号的连续性又能保留高频信号实现降噪。

1 小波综合阈值降噪方法
1.1 小波综合阈值降噪方法原理

阈值函数法也称小波阈值降噪方法。Donoho等人已经证明小波阈值降噪方法优于其它经典降噪方法(Bruni等,2006)。目前,常用的阈值降噪方法包括软、硬阈值降噪方法。硬阈值方法将信号小波系数绝对值与小波系数阈值比较,实现信号高频部分小波系数的保留,但在阈值置零处易出现不连续现象,造成有效信号缺失(耿冠世等,2015Mousavi等,2016),硬阈值函数如公式(1)所示;软阈值方法改善硬阈值方法中出现的信号缺失现象,但损失高频信号能量,软阈值函数如公式(2)所示。

$ {\rm{hard}}\left({\omega, \lambda } \right) = \left\{ {\begin{array}{*{20}{l}} {\omega, \;\;\left| \omega \right| \ge {\rm{ }}\lambda }\\ {{\rm{0, }}\;\;\;\left| \omega \right| < \lambda {\rm{ }}} \end{array}} \right. $ (1)
$ {\rm{soft}}(\omega, \lambda) = \left\{ {\begin{array}{*{20}{l}} {{\rm{sgn}}(\omega)\left({\left| \omega \right| - \lambda } \right), \left| \omega \right| \ge {\rm{ }}\lambda }\\ {{\rm{0}}, {\rm{ }}\left| \omega \right| < \lambda {\rm{ }}} \end{array}} \right. $ (2)

式中,ω为信号小波系数;λ为阈值(非负值);sgn为符号函数,当含噪信号大于0时,sgn=1;当含噪信号小于0时,sgn=-1。针对软阈值方法在降噪过程中造成高频信号能量损失的不足,小波综合阈值方法结合硬阈值方法提高信号高频部分能量并对信号低频部分用软阈值方法保持信号连续性的优势,提出小波综合阈值函数,如公式(3)所示。

$ {\rm{new}}(\omega, \lambda) = \left\{ {\begin{array}{*{20}{l}} {\omega, {\rm{ }}\left| \omega \right| \ge {\rm{ }}\lambda {\rm{ }}}\\ {{\rm{sgn}}(\omega)(\left| {\omega - \lambda } \right|), {\rm{ }}\left| \omega \right| < \lambda } \end{array}} \right. $ (3)

本文提出的小波综合阈值方法构造的新函数${\rm{new}}(\omega, \lambda) = {\hat \omega _{j, k}} $,需要满足${\hat \omega _{j, k}} = {\omega _{j, k}} $,使改进的小波系数与真实小波系数相近,避免小波重构时损失有效信号。由此,构造出小波综合阈值函数如公式(4)所示:

$ {\hat \omega _{j, k}} = \left\{ {\begin{array}{*{20}{l}} {(1 - b){\rm{sgn(}}{{\hat \omega }_{j, k}})(\left| {{\omega _{j, k}} - \lambda } \right|) + {b_{{\omega _{j, k}}}}, {\rm{ }}\omega < \lambda }\\ {{\omega _{j, k}}, {\rm{ }}\omega \ge {\rm{ }}\lambda } \end{array}} \right. $ (4)

其中,$ {\omega _{j, k}}$为真实小波系数,${\hat \omega _{j, k}} $为改进小波系数,令$b = \frac{{\left| {{\omega _{j, k}}} \right| - \lambda }}{{\left| {{\omega _{j, k}}} \right|}} $,当$ \left| {{\omega _{j, k}}} \right|$= $ \lambda $时,b=0,$ {\hat \omega _{j, k}}$=0;当$ {\omega _{j, k}}$= $\lambda $时,b=0,${\hat \omega _{j, k}} $=0,解决硬阈值方法在阈值置零部分不连续现象。随着$\left| {{\omega _{j, k}}} \right| $的增大,b=1,$ {\hat \omega _{j, k}} \approx {\omega _{j, k}}$,解决软阈值方法高频信号能量损失的问题。小波综合阈值降噪通过改变阈值函数保留高频信号的同时保持信号的连续性。实验结果表明,改进的阈值函数对地震信号降噪有很好的效果。

SNR(信噪比)和MSE(平均方差)是评定降噪方法优劣的一种方式,假定地震信号向量为a=[a0a1a2,…aN-1]T,则有公式(5):

$ {a_i} = {f_i} + {n_i}\;\;\;i = 0, 1, 2, \cdots \cdots, \left({N - 1} \right) $ (5)

其中,fi为函数f的抽样,ni是分布为$N\left({0, \sigma } \right)$的高斯白噪声。降噪的目标是拟合值$ \mathop f\limits^ \wedge $的平均方差MSE最小,MSE的表达式如公式(6)所示:

$ {\rm MSE} = \frac{1}{N}{\left| {\hat f - f} \right|^{\rm{2}}} = \frac{{\rm{1}}}{N}\sum\limits_{i = 0}^{N - {\rm{1}}} {({{\hat f}_i}} - {f_i}{)^2} $ (6)

由式(6)可知,${\hat f_i} - {f_i} $越小,则MSE值越小,降噪后波形与原信号波形越近似。

1.2 小波特征能量谱系数

小波特征能谱系数是降噪方法的表征方式,能直观地观察信号在低频和高频部分的能量分布,便于观察降噪结果快速得出结论。小波特征能谱系数经过i个尺度分解后总能量不变,如公式(7),其中f(n)为地震信号离散采样序列,A为信号中低频部分,D为信号中高频部分,Aif(n)、Dif(n)为尺度变换后各个频率的分量,EiAf(n)、EiDf(n)分别为在分解尺度i上的低频信号分量能量和高频信号分量能量。

$ f(n) = {A_i}f(n) + {D_1}f(n) + \cdots \cdots + {D_N}f(n)\;\;\;i = 1, 2, \cdots \cdots, N $ (7)
$ {E_i}^Af(n) = \sum\limits_{i = 1}^N {({A_i}} f(i){)^2} $ (8)
$ {E_i}^Df(n) = \sum\limits_{i = 1}^N {({D_i}} f(n){)^2}{\rm{ }} $ (9)

定义每个小波分解中分量的能量与总能量之间的比值,即小波特征能谱系数,分别用参数$ h{E_i}^A$$h{E_i}^D $表示,即:

$ h{E_i}^A = \frac{{{E_i}^Af(n)}}{{Ef(n)}}\;\;\;\;i = 1, 2, \cdots \cdots, N $ (10)
$ h{E_i}^D = \frac{{{E_i}^Df(n)}}{{Ef(n)}}\;\;\;i = 1, 2, \cdots \cdots, N $ (11)
2 小波综合阈值降噪实验
2.1 模拟地震信号降噪实验

实验研究处理的信号针对井下近震信号频段,近震信号以coif小波为小波基函数,并计算信号在6次分解后的小波特征能谱系数。在第3次分解尺度上的特征能谱系数中能观察出近震信号能量较强,因此选择3次分解上的小波特征能谱系数。井下近震信号采集过程中包括近震信号和噪声信号,近震信号峰值能量的频率主要集中在3—6Hz。根据随机噪声来源和噪声自身表现规律,将噪声划分为3类(表 1):

表 1 噪声分类 Table 1 Noise classification

为比较小波综合阈值方法与软阈值方法对次生噪声及环境噪声的降噪能力,选取与近震信号具有相似小波系数特征的雷克子波信号进行模拟实验。雷克子波信号添加噪声频率范围为3—30Hz,包括环境噪声和次生噪声。图 1(a)为雷克子波波形及加噪雷克子波信号波形,对加噪雷克子波波形进行软阈值和小波综合阈值降噪实验如图 1(b)所示。


图 1 雷克子波波形及降噪波形 Fig. 1 Ricker wavelet and wavelet after noise reduction

图 1(a)中加噪后的雷克子波信号初至到时为210s,加噪后波形高频信号被压制,无法分辨出地震信号与噪声信号。对比图 1(b)中2种降噪方法的波形,软阈值方法压制高频信号振幅,零频附近噪声与加噪后波形频率相似,降噪作用不明显;小波综合阈值方法提高高频信号振幅,降低噪声在零频时振幅。通过计算SNR和MSE(表 2)可知,小波综合阈值方法在2项指标上有所改进。小波综合阈值方法降噪后MSE值最小,降噪后信号与原信号更近似。

表 2 仿真降噪后所得结果的SNR和MSE值 Table 2 SNR and MSE from the de-noised signal in simulation

雷克子波波形能谱系数如图 2(a)所示,高频信号能量集中在第2次分解,噪声信号能量集中在第5、6、7次分解。加噪后波形能谱系数如图 2(b)所示,波形中高频信号能量被噪声分解,无法分辨高频信号能谱分布。利用软阈值方法对加噪雷克子波信号进行降噪处理(图 2(c)),该方法中高频信号能量集中在第1次分解,与原始信号波形能谱系数分布不符,压制高频信号能量,降噪效果不明显。小波综合阈值方法处理加噪波形结果如图 2(d),该方法中高频信号能谱系数分布与原信号相似,集中在第2次分解。小波综合阈值方法提高原信号中第2次分解的高频信号能量,抑制噪声信号在各次分解中的能量,有效实现降噪。


图 2 雷克子波波形图及降噪方法能谱系数 Fig. 2 Ricker wavelet and energy-spectrum coefficient of noise reduction method
2.2 实际地震信号降噪实验

为验证小波综合阈值方法对实际地震数据处理效果,截取河南省周口市太康县逊母口镇地震波信号进行小波综合阈值滤波实验。太康县逊母口镇地震台站位于河南省周口市2条断裂构造的交会处,台站选择330m井深进行地震监测。

井下高频地震计数据采样频率为1024Hz,采样通道数为6道,记录长度为12s。为方便计算,抽取第2通道0—120s的数据如图 3(a)。地震信号能谱系数如图 3(b),利用软阈值方法对实际地震信号进行降噪处理,信号能谱系数如图 3(c),小波综合阈值和基于S变换的软阈值降噪后波形的能谱系数分别如图 3(d)图 3(e)。实际地震信号、软阈值波形与原信号对比波形如图 3(f)


图 3 逊母口镇软阈值处理图及降噪方法对比图 Fig. 3 Comparison of seismic waveforms with different noise reduction methods in Xunmukou county

图 3(b)中高频信号能量集中在第1次和第2次分解,噪声信号分布在第3次分解后。对比图 3(b)图 3(c),软阈值方法对实际信号降噪后,高频信号能量与原信号高频信号能量分布相似,对噪声信号降噪效果不明显。对比图 3(b)图 3(d),小波综合阈值降噪方法增大了实际信号第1次和第2次分解的高频信号能量,小波综合阈值方法对高频部分作用明显,抑制低频信号能量,实现实际地震信号降噪。对比图 3(d)图 3(e)基于S变换的软阈值方法同样能实现地震信号的降噪。观察图 3(f),由地震波运动学原理可知,初至波由于传播距离较短、能量强、衰减慢,表现为具有高频能量,据此可以判断地震波的初至到时在4.2s左右,地震波中有效信号在4.5—5s之间。6—12s时由于多次阈值分解使信号中夹杂的高频噪声被滤除,从能谱系数分布上可以看到第5次分解后信号能量降低且稳定。通过观察图 3(f)左图得知,软阈值方法丢失信号的初至到时,导致信号失真;观察图 3(f)右图得知,小波综合阈值方法和基于S变换的软阈值方法能保留信号的初至到时且能较完整的重构出原信号。通过实际地震数据试算可知,小波综合阈值降噪能更清晰、直观地反映地震数据的局部信息特征,降低次生噪声和环境噪声对地震信号带来的干扰和误判。

3 结果与讨论

小波综合阈值降噪方法利用硬阈值方法提高高频信号能量,对信号低频部分利用软阈值方法保留信号光滑性,在提高信号降噪能力的同时保证信号连续性。本文分别用软阈值和小波综合阈值方法对加噪后雷克子波信号进行降噪处理,利用软阈值、小波综合阈值方法和基于S变换的软阈值方法对信号进行降噪处理,观察降噪后波形及能谱系数。实验表明,利用小波综合阈值方法降噪后雷克子波波形高频信号得到恢复,噪声能谱系数被压制;小波综合阈值降噪后波形MSE值最小,降噪后波形与原信号波形最近似。此外,通过对实际地震数据进行小波综合阈值降噪分析,能细致判断实际波形的初至时间和有效信号出现的时间范围。地震信号高频部分经过小波综合阈值降噪后能量更集中在第1次和第2次分解。

通过仿真实验和实际波形降噪实验验证,小波综合阈值方法与软阈值降噪方法相比,能有效消除次生噪声和环境噪声对地震信号带来的干扰,降低地震误判和漏判,其降噪效果与基于S变换的软阈值降噪方法同样明显。但是,小波综合阈值方法存在大量数据处理缓慢的问题,应改进算法上存在的冗余问题或进一步提出阈值函数的改进。改进后的小波综合阈值方法应能适应与近震信号有相似特征的信号,增加应用的广泛性与普遍性。

参考文献
耿冠世, 俞言祥, 2015. 中国西部地区震源破裂尺度与震级的经验关系[J]. 震灾防御技术, 10(1): 68-76. DOI:10.11899/zzfy20150107
侯跃伟, 赵兵, 田韬, 2015. 基于Daubechies小波分析的南京数字化钻孔形变震前变化特征研究[J]. 震灾防御技术, 10(2): 388-396. DOI:10.11899/zzfy20150219
孔祥茜, 吴继伟, 岳继光, 2005. 地震信号小波变换的去噪方法[J]. 计算机辅助工程, 14(3): 52-56.
范涛, 赵兆, 吴海, 2014. 矿井瞬变电磁多匝回线电感影响消除及曲线偏移研究[J]. 煤炭学报, 39(05): 932-940.
刘霞, 潘洪屏, 高晓春, 2010. 基于小波阈值的地震信号去噪处理[J]. 科学技术与工程, 10(29): 7251-7254. DOI:10.3969/j.issn.1671-1815.2010.29.030
李英, 张淑贞, 许康生, 2006. 小波降噪方法在地震信号处理中的应用[J]. 西北地震学报, 28(2): 159-162.
曲中党, 吴蔚, 贺日政, 等, 2015. 基于S变换的软阈值滤波在深地震反射数据处理中的应用[J]. 地球物理学报, 58(9): 3157-3168. DOI:10.6038/cjg20150912
唐守锋, 童敏明, 潘玉祥, 等, 2011. 煤岩破裂微震信号的小波特征能谱系数分析法[J]. 仪器仪表学报, 32(7): 1521-1527.
魏学强, 袁洪克, 秦晶晶, 等, 2016. 广义S变换地震信号时频分析[J]. 震灾防御技术, 11(4): 808-813. DOI:10.11899/zzfy20160411
曾宪伟, 赵卫明, 师海阔, 等, 2010. 利用小波包变换对地震信号进行时频分析时小波基函数的选取[J]. 地震研究, 33(4): 323-328.
郑作亚, 卢秀山, 李克行, 2007. 一类小波基函数的构造及其在测量数据处理中的应用探讨[J]. 测绘科学, 32(3): 9-11.
Bruni V., Vitulano D., 2006. Wavelet-based signal de-noising via simple singularities approximation[J]. Signal Processing, 86(4): 859-876. DOI:10.1016/j.sigpro.2005.06.017
Mousavi S. M., Langston C. A., 2016. Hybrid seismic denoising using higher-order statistics and improved wavelet block thresholding[J]. Bulletin of the Seismological Society of America, 106(4): 1380-1393. DOI:10.1785/0120150345
Xia L. M., Zheng W., Liang M. D., 2017. Application of a comprehensive one dimensional wavelet threshold denoising algorithm in blasting signal. In: Balas V. E., Jain L. C., Zhao X. M., Information Technology and Intelligent Transportation Systems. Cham: Springer, 517-527.


A Wavelet Comprehensive Threshold Seismic Signal Denoising Method
Yang Jingsong*, Wang Han, Yao Zhenjing, Han Zhiming, Li Yanan, Gao Qiang
(Vibration Precision Measurement Technology and Instrument Key Laboratory, Institute of Disaster Prevention, Sanhe 065201, Hebei, China)
Abstract

The wavelet threshold method with hard or soft threshold is common for noise reduction of seismic signal, but it is often to cause the loss of high frequency information in the signal to lead in earthquake misjudgment or missing judgment. Taking the advantages of hard threshold and soft threshold denoising method, the wavelet comprehensive threshold improves the high frequency signal energy by the hard threshold treatment, and keeps the signal smoothness and reduce the signal degradation of noise capacity while ensuring signal continuity by lowness of the signal. Based on the small wavelet coefficients of the noise signal and the large wavelet coefficients of the seismic signal, the simulation experiment of the Ricker signal and the noise reduction experiment is carried out. The simulation results show that the MSE of the wavelet comprehensive threshold is the closest to the original signal waveform after noise reduction, and the energy of the high frequency part of the waveform is reduced and the low frequency part is suppressed. Finally, taking the noise reduction of the natural waveform as an example, we found that noise signal energy is suppressed and the first time of the waveform can be obtained.



主办单位:中国地震台网中心
版权所有:《震灾防御技术》编辑部
地址:北京西城区三里河南横街5号,   邮编:100045
邮箱:zzfy2006@126.com   电话:59959251;59959139;59959462
访问人数:1518769
一种小波综合阈值地震信号降噪方法
杨敬松*, 王晗, 姚振静, 韩智明, 李亚南, 高强
《震灾防御技术》, DOI:10.11899/zzfy20180218