文章快速检索    
  地震地磁观测与研究  2019, Vol. 40 Issue (5): 70-78  DOI: 10.3969/j.issn.1003-3246.2019.05.010
0

引用本文  

朱杰. 元谋地电场观测数据小波阈值去噪[J]. 地震地磁观测与研究, 2019, 40(5): 70-78. DOI: 10.3969/j.issn.1003-3246.2019.05.010.
Zhu Jie. Wavelet threshold de-noising of Yuanmou geoelectric field observation data[J]. Seismological and Geomagnetic Observation and Research, 2019, 40(5): 70-78. DOI: 10.3969/j.issn.1003-3246.2019.05.010.

作者简介

朱杰(1986—), 男, 助理工程师, 本科, 主要从事地震监测、系统运维工作。E-mail:583805982@qq.com

文章历史

本文收到日期:2018-05-30
元谋地电场观测数据小波阈值去噪
朱杰     
中国昆明 650224 云南省地震局
摘要:为了解决元谋地电场观测数据的噪声干扰问题,采用小波阈值去噪的方法,根据元谋地电场观测数据特征,选取适当的阈值及不同的小波函数、分解层数进行研究。综合各小波函数去噪效果,通过对相对误差、均方根误差、信噪比等指标进行评价,确定适合元谋地电场观测数据去噪的阈值为软阈值函数、分层阈值;最优小波函数分别为Sym4、Sym5、db6和db8;分解层数分别为5、7、9。最后应用研究结果处理实际地电场观测数据,取得了较好的去噪效果。
关键词地电场    小波函数    分解层数    阈值去噪    误差分析    
Wavelet threshold de-noising of Yuanmou geoelectric field observation data
Zhu Jie     
Yunnan Earthquake Agency, Kunming 650224, China
Abstract: In order to solve the problem of noise interference in Yuanmou geoelectric field observation data, this paper adopts the method of wavelet threshold denoising to study by choozing appropriate threshold, different wavelet functions and decomposition layers according to the characteristics of Yuanmou geoelectric field data. Based on the evaluation of relative error, root mean square error and signal-to-noise ratio, the threshold of Yuanmou geoelectric field denoising is determined to be soft threshold function and hierarchical threshold. The optimal wavelet functions are Sym4, Sym5, db6 and db8, and the decomposition layers are 5, 7 and 9 respectively. Finally, the actual geoelectric field data are processed with the results of the study and good denoising effect is achieved.
Key words: geoelectric field    wavelet function    decomposition layers    threshold de-noising    error analysis    
0 引言

地电场是一种重要的地球物理场,主要包括两部分:大地电场和自然电场。前者源于外太空电离层的电流体系及其变化,分布在整个地表的较大区域;后者源于地下场源,具有相对的局部稳定性(傅承义等,1985)。我国将地电场用于地震地球物理观测始于1966年邢台地震之后,自此积累了一定的观测资料。但是由于地电场观测的场地范围较大,容易受到场地环境因素以及其他不明因素的干扰(席继楼等,2002),而这些干扰的频段与地电场信号的频段有大量的重叠和交复,因此通过物理手段滤波难以有效消除。邱颖等(2009)通过小波方法研究了地电场观测中去除地电阻率观测干扰的方法。李伟等(2013)利用小波变换有效地滤除了上海地电场观测中的地铁噪声干扰。梁跃等(2015)利用小波包去噪对密山地电场观测数据去噪,有效还原了地震发生后地电场的变化信息。孙雷等(2013)通过FFT、小波分析等方法研究了新沂台地电场的频谱特征,并得出其优势周期。

通过小波阈值去噪方法处理地电场观测数据,可以有效地滤除地电场信号中的高频噪声,同时信号中不同频率的异常干扰信号不会遭到破坏(李伟等,2013)。然而,目前对于研究地电场所用小波函数的选取和分解层数没有系统的理论方法,需要根据所选研究数据的特点,依据经验不断比较测试,才能找到适合的小波函数和分解层数(李月琴等,2008)。而小波函数的选取也是小波分析运用中的难点,不同数据需要选择不同的小波函数,而相同数据选取不同小波函数的分析结果也不尽相同(梁跃等,2015)。本文利用元谋地电站ZD9A-2B型地电仪所记录的分钟采样数据,对常见的几个小波函数族进行比较分析,并利用相对误差、均方根误差、信噪比等手段比较不同小波函数、分解层数和阈值函数在信号去噪时的优劣,最终确定了适合元谋台地电场观测数据小波阈值去噪的阈值、最优小波函数和分解层数。

1 小波阈值去噪原理及评价方法 1.1 小波变换及其重构

f (t)∈L2(R),即f (t)为平方可积函数,其连续小波变换为

$ W{T_f}(a, b) = \frac{1}{{\sqrt {|a|} }}\int_{ - \infty }^{ + \infty } f (t)\mathit{\Psi }*\left({\frac{{t - b}}{a}} \right){\rm{d}}t\quad (a \ne 0) $ (1)

式中,$\frac{1}{{\sqrt {|a|} }}\mathit{\Psi }*\left({\frac{{t - b}}{a}} \right) = {\mathit{\Psi }_{a, b}}(t)$,被称为由基本小波Ψ(t)生成的尺度伸缩和平移。其中,a为伸缩因子;b为平移因子。

f (t)∈L2(R),其傅里叶变换为${\mathit{\hat \Psi }(\omega)}$,当${\mathit{\hat \Psi }(\omega)}$满足容许条件${C_\mathit{\Psi }}(\omega) = \int_0^\infty {\frac{{|\mathit{\hat \Psi }(\omega){|^2}}}{{|\omega |}}} {\rm{d}}\omega < \infty $时,则其重构公式为

$ f(t) = \frac{1}{{{C_\mathit{\Psi }}}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\frac{1}{{{a^2}}}} } {W_f}(a, b)\mathit{\Psi }\left({\frac{{t - b}}{a}} \right){\rm{d}}a{\rm{d}}b $ (2)
1.2 小波阈值去噪方法

现实中,受外界环境的干扰,信号的采集和传输过程中,无法避免地会产生噪声。为了提高信号的可信度,在使用信号前必须进行降噪处理。地电场观测信号是非平稳信号,其中含有的低频信号为地电场信号,高频信号为噪声干扰信号。小波阈值去噪过程就是利用小波分解将地电场信号(含噪声)分解到各尺度中,然后设定一个阈值,对各尺度小波系数进行阈值处理,把不同尺度中属于噪声的模值较小的系数去除,保留并增强属于地电场信号的模值较大的系数,再利用小波逆变换将处理后的系数重构得到去噪后的地电场信号。小波阈值去噪流程见图 1

图 1 小波阈值去噪流程 Fig.1 Wavelet threshold de-noising flow chart
1.3 小波阈值去噪质量评价

本文使用相对误差(error)、均方根误差(RMSE)和信噪比(SNR)等3个指标(陶珂等,2012)对去噪的质量进行评价,找出适合元谋台地电场观测数据处理的最优结果。

相对误差(error)指重构后的信号与原始信号之间的差异

$ {\rm{ error }} = \frac{{|\hat f(i) - f(i)|}}{{f(i)}} $ (3)

式中,f(i)为原始信号;${\hat f(i)}$为重构信号。相对误差越小,则表示重构信号与原始信号间的相似度越高。

均方根误差(RMSE)指重构后的信号与原始信号的均方根误差

$ {\rm{RMSE}} = \sqrt {\frac{1}{n}\sum\nolimits_{i = 1}^n {{{[f(i) - \hat f(i)]}^2}} } $ (4)

RMSE越小,则表示去噪效果越好。

信噪比(SNR)指信号能量与噪声能量的比值,公式如下

$ {\rm{SNR}} = 101{\rm{g}}\frac{{{P_{\rm{S}}}}}{{{P_{\rm{n}}}}} $ (5)

式中,${P_{\rm{S}}} = \frac{1}{n}\sum\nolimits_n {{f^2}} (n)$,为信号的功率;${P_{\rm{n}}} = \frac{1}{n}\sum\nolimits_n {{{[f(n) - \hat f(n)]}^2}} $为噪声的功率。SNR越大,去噪效果越好。

2 小波函数、分解层数及阈值 2.1 小波函数基本性质

分析不同信号需要使用不同的小波函数。小波函数一般具有以下几种基本性质。

(1)支撑长度。是指当时间或频率趋于无穷大时,小波函数、尺度函数从一个有限值收敛到0的长度。支撑长度过长或过短均不利于计算和信号能量集中,需要对其适当选择(朱吉昌,2016)。

(2)消失矩。其定义为

$ \int {{t^p}} \mathit{\Psi }(t){\rm{d}}t = 0 $ (6)

式中,Ψ(t)为基本小波;当0≤pN时,则称小波具有N阶消失矩。消失距和支撑长度一样,需要根据信号特征折中选取(梁跃等,2015)。

(3)对称性。若小波函数关于某点对称,则称其具有对称性。对称性可以有效地避免信号分解、重构中的相位畸变。

(4)正则性。小波函数的正则性代表小波函数连续可微。正则性好,则在信号的重构中可以获得较好的平滑效果。其与消失距有一定关系,同样需要折中选取。

(5)相似性。所选小波函数的波形与分析信号的波形相似度越高,信号去噪的适应性越强。

常用的几种小波函数参数特性见表 1

表 1 常用小波基参数特性比较 Table 1 Parameter comparison of different wavelet basis
2.2.1 适合元谋地电场观测信号去噪的小波函数选取

元谋台地电场观测场地位于农田内,场地内机井、金属管道、铁丝网等干扰源较多,农业生产活动频繁,信号中含有较多的高频噪声。所以,适合信号去噪的小波函数需要具有较高的消失矩、适中的支撑长度、正则性较好、对称性好等特性。通过分析,可以在常用的小波函数中选取Haar、dbN、CoifN和SymN几组小波函数族来作进一步分析:Haar函数即为N=1的db小波,它是支撑域在t∈[0, 1]范围内的单个矩形波,在时域上是不连续的,不适用于地电场信号的分析。Daubechies(dbN)函数,随着db小波序列N的增大,平滑性变好,在频域上局部细节分辨率增强,划分频带的效果好。Coiflets(Coif N)函数,只有Coif NN =1,2,3,4,5)这一系列具有比db小波更好的对称性。Symlets(SymN)函数,是在db小波基础上改进的一种小波函数。

综合各小波的性质以及表 1中支撑长度和消失距的分析,可以初步确定,dbN小波、SymN小波使用N=3—9阶及CoifN小波使用N=1—3阶进行地电场观测信号的分解重构。

2.2.2 元谋地电场观测信号去噪小波分解层数的确定

小波分解过程是迭代的,理论上可以无限分解下去(邱颖,2008刘明才,2013)。随着分解层数的增加,信号与噪声的细节特征差异会越显著,更便于从信号中分离噪声。但是,如果分解层数取值过大,重构信号中的部分有效信号也会受到影响而失真,从而影响去噪的效果。所以,分解层数的选择需要根据信号的频率特点综合考虑。小波分解层数的频段范围与采样频率有关。若分解为N层,则各层频段大小为Nyquist/2N,其中,Nyquist为奈奎斯特频率(万永革,2012),其大小为采样频率的一半。

选取元谋地电站2016年NS向长极距测道分钟采样数据作快速傅里叶变换(FFT),结果见图 2,可见信号能量主要集中在1.157×10-5 Hz、2.314×10-5 Hz、3.473×10-5 Hz、4.63×10-5 Hz处,这几个峰值点对应的周期分别为24 h、12 h、8 h、6 h,优势周期较明显。由采样定理可知,该信号的奈奎斯特频率为1/120 Hz,所以元谋地电场观测信号的能量为1.000×10-5—1.000×10-3 Hz。由公式Nyquist/2N可以算出各分解层数所逼近的细节频率(表 2)。由表 2可以初步筛选得到分解层数为4—9时较适合元谋地电场观测信号的处理。

图 2 元谋地电场NS向长极距频谱 Fig.2 The long polar distance spectrum diagram of NS direction of the Yuanmou geoelectric field
表 2 各分解层数逼近频率 Table 2 The approximate frequency of each decomposition layers
2.2.3 元谋地电场信号去噪最优小波函数和分解层数确定

综上分析,本文使用相对误差(error)指标来确定元谋地电场信号去噪的最优小波函数和分解层数。选取2016年NS向长极距数据,对dbN小波族、SymN小波族和Coif N小波族分别进行4—9层分解,通过相对误差评估最优结果,分析结果见表 3表 5

表 3 db小波族相对误差分析 Table 3 Analysis of relative error of db wavelet family
表 4 Sym小波族相对误差分析 Table 4 Analysis of relative error of Sym wavelet family
表 5 Coif小波族相对误差分析 Table 5 Analysis of relative error of Coif wavelet family

表 3表 5的相对误差数据分析可知:①Sym小波族对信号的去噪效果最优,其中的最优小波函数为Sym4小波,最优分解层数为5层;②同理可得,NS向短极距最优小波函数为db8小波,最优分解层数为9层;EW向长极距最优小波函数为db6小波,最优分解层数为7层;EW向短极距最优小波函数为Sym5小波,最优分解层数为7层,鉴于篇幅,具体数据不一一列出。

2.3 元谋地电场观测信号去噪阈值选取

目前常用的阈值可以分为全局阈值和分层阈值2类。全局阈值是对各层所有的小波系数选用同一个阈值进行处理,而分层阈值是根据各层不同尺度分别选取阈值进行处理。

2.3.1 全局阈值处理中阈值的选取

Matlab中全局阈值处理常用的阈值选择方法有:无偏风险阈值估计(rigrsure)、启发式阈值估计(heursure)、固定阈值估计(sqtwolog)和极值阈值估计(minimaxi)等4种。当噪声在信号的高频段分布较少时,极值估计法和无偏风险估计法去噪效果较好,但针对地电场观测信号的去噪显然不合适。固定阈值估计法和启发式阈值估计法去噪较彻底,二者的阈值取值相同,但固定阈值估计法去噪后平滑性较好,所以在全局阈值处理中选择固定阈值估计法(安蕴,2012)。

2.3.2 全局阈值与分层阈值的对比

选取元谋地电场2016年NS向长极距数据进行分析,使用Sym4小波,分解重构层数为5层,通过均方根误差(RMSE)和信噪比(SNR)2个指标,比较全局阈值方法和分层阈值方法的优劣。其中,全局阈值${\mathop{\rm th}\nolimits} r = \sigma \sqrt {2\lg N} $σ为噪声的标准偏差;N为信号长度。分层阈值thr1由Matlab中wdcbm函数算出,结果见表 6

表 6 不同阈值的比较 Table 6 Comparison of different threshold

表 6可见,分层阈值处理时RMSE较小,SNR较大,相较于全局阈值处理更适合于元谋地电场信号的去噪。

2.3.3 软、硬阈值函数选取

确定了噪声小波系数阈值的门限后,就需要用阈值函数对这个含有噪声的小波系数进行过滤,除去噪声系数。常用的阈值函数有软阈值、硬阈值函数(Donoho,1995)。

软阈值函数的表达式为

$ f(x) = \left\{ {\begin{array}{*{20}{l}} {{\mathop{\rm sgn}} (x)(|x| - \lambda)}&{|x| \ge \lambda }\\ 0&{|x| < \lambda } \end{array}} \right. $ (7)

其中,sgn(x)为符号函数:${\mathop{\rm sgn}} (x) = \left\{ {\begin{array}{*{20}{l}} 1&{x > 0}\\ { - 1}&{x < 0} \end{array}} \right.$。当小波系数的绝对值小于给定的阈值时,令其为0;大于等于阈值时,令其减去阈值。

硬阈值函数的表达式为

$ f(x) = \left\{ {\begin{array}{*{20}{l}} x&{|x| \ge \lambda }\\ 0&{|x| < \lambda } \end{array}} \right. $ (8)

当小波系数的绝对值小于给定的阈值时,令其为0;大于等于阈值时,令其保持不变。

使用阈值函数的目的是使有效信号的小波系数得以保留,而噪声的系数被剔除,并且在小波域内系数的过渡要平滑。硬阈值函数去噪后得到的信噪比较高,在均方误差层面上优于软阈值函数,但会有局部的抖动现象,会出现附加震荡,产生跳跃点,不具有原始信号的平滑性。软阈值函数去噪后得到的信噪比不如硬阈值函数,但是得到的小波系数整体连续性较好,最后的结果看起来较平滑,重构信号不会出现附加震荡。综合考虑元谋地电场观测信号的特点,为保持其原始信号的平滑性,这里选用软阈值函数进行处理。

通过上述分析可以确定,NS向长极距信号阈值选取为软阈值函数,分层阈值处理;同理,其余3个测道阈值选取也为软阈值函数,分层阈值处理。

3 实际数据处理

选取2016年2—6月元谋地电场NS向长极距数据分析,使用Sym4小波、分解5层、软阈值函数、分层阈值处理,对比去噪前后信号的差异,结果如图 3所示。

图 3 原始信号与去噪信号对比 (a)2016年2—6月原始地电场观测信号;(b)2016年2—6月去噪后地电场观测信号;(c)噪声波形 Fig.3 Comparison between original signal and de-noising signal

图 3可见,去噪重构后信号还原度较高,与原始信号的相似度较好,同时细节处高频噪声得以消除,而信号中较大的突跳等异常干扰信号得以保留。分别对原始地电场观测信号和去噪重构后的信号作小波时频分析,结果如图 4图 5所示。

图 4 地电场原始观测信号小波时频 Fig.4 Wavelet time-frequency graph of original signal
图 5 去噪后信号小波时频 Fig.5 Wavelet time-frequency graph of de-noising signal

图 4可见,2016年5月、6月农忙时节原始信号频域内高频噪声的干扰较为严重,分布区间广,异常干扰频率分辨率较低,低频段信号分辨率较低。由图 5可以看出,去噪后信号在频域内高频噪声得到较好的抑制,同时对各时段内干扰信号的频率成分保留较好,对低频段信号的分辨率有一定程度的提升。

4 结论

(1)小波阈值去噪的最优小波函数、分解层数的选取无固定方式,需要根据具体信号的特点、性质进行选择。

(2)元谋地电场NS向长极距小波阈值去噪最优小波函数为Sym4小波,最优分解层数为5层;NS向短极距最优小波函数为db8小波,最优分解层数为9层;EW向长极距最优小波函数为db6小波,最优分解层数为7层;EW向短极距最优小波函数为Sym5小波,最优分解层数为7层。

(3)元谋地电场观测数据小波阈值去噪选用软阈值函数、分层阈值去噪处理较为合适。

(4)小波阈值去噪后信号中的高频噪声被较好地抑制,同时各时段内干扰信号的频率成分得以较好地保留,对低频段地电场观测信号频率的分辨率有一定程度的提升。

参考文献
安蕴.基于MATLAB的小波分析用于地震信号的去噪研究[D].太原: 中北大学, 2012: 40-43. http://d.wanfangdata.com.cn/Thesis/D315967
傅承义, 陈运泰, 祁贵仲. 1985. 地球物理学基础[M]. 北京: 科学出版社, 447.
李伟, 马钦忠, 宋志平, 等. 2013. 小波变换在地电场数据分析中的应用[J]. 地震学报, 35(1): 26-35. DOI:10.3969/j.issn.0253-3782.2013.01.004
李月琴, 栗苹, 闫晓鹏, 等. 2008. 无线电引信信号去噪的最优小波基选择[J]. 北京理工大学学报, 28(8): 723-726.
梁跃, 董桂菊, 崔天时, 等. 2015. 密山地电场数据小波包去噪研究[J]. 地震地磁观测与研究, 36(1): 84-88. DOI:10.3969/j.issn.1003-3246.2015.01.014
梁跃, 李太岩, 王海龙, 等. 2015. 地电场信号去噪小波函数的选取[J]. 防灾减灾学报, 31(2): 42-46.
刘明才. 2013. 小波分析及其应用[M]. 北京: 清华大学出版社, 30-32.
邱颖.地电场观测中已知源干扰抑制研究[D].北京: 中国地震局地震预测研究所, 2008: 21-22. http://cdmd.cnki.com.cn/Article/CDMD-85405-2008129230.htm
邱颖, 席继楼. 2009. 小波方法在地电场干扰处理中的分析研究[J]. 地震, 29(2): 57-63.
孙雷, 李飞, 杨冯威. 2013. 新沂台地电场频谱特征的分析与研究[J]. 华南地震, 33(2): 93-103. DOI:10.3969/j.issn.1001-8662.2013.02.012
陶珂, 朱建军. 2012. 小波去噪质量评价方法的对比研究[J]. 大地测量与地球动力学, 32(2): 128-133.
万永革. 2012. 数字信号处理的MATLAB实现[M]. 2版. 北京: 科学出版社, 60.
席继楼, 赵家骝, 王燕琼, 等. 2002. 地电场观测技术研究[J]. 地震, 22(2): 67-73.
朱吉昌.测井曲线小波分析方法的改进及其在地层对比中的应用[D].北京: 中国地质大学(北京), 2016: 16-17. http://cdmd.cnki.com.cn/Article/CDMD-11415-1016183925.htm
Donoho D L. 1995. De-noising by soft-thresholding[J]. IEEE Transactions on Information Theory, 41(3): 613-627. DOI:10.1109/18.382009