適于泥石流除噪的EMD聯合小波閾值除噪方法

朱鳳杰，焦瑞莉，滕鵬曉

?

適于泥石流除噪的EMD聯合小波閾值除噪方法

朱鳳杰1，焦瑞莉1，滕鵬曉2

(1. 北京信息科技大學信息與通信工程學院，北京 100085；2. 中國科學院聲學研究所，北京 100190)

次聲傳感器采集到的泥石流次聲信號中包含有大量的無關干擾信號，嚴重影響信號的分析與評估。針對含噪泥石流信號中無法準確確定噪聲頻段的特點，以及傳統經驗模態分解(Empirical Mode Decomposition, EMD)聯合小波閾值去噪方法無法智能分辨噪聲所在頻段的缺點，提出了信號經EMD分解后，基于相關性選擇噪聲頻段的方法。首先利用EMD分解獲取信號的固有模態函數(Intrinsic Mode Function, IMF)分量，然后計算各個IMF分量與原始信號的相關性，根據相關性大小確定IMF噪聲頻段，然后采用小波閾值去噪方法對噪聲頻段進行處理，最后對處理后的信號進行重構得到去噪泥石流信號。通過模擬實驗分析，證明該方法具有智能選擇噪聲頻段的能力，是一種更適于泥石流信號的去噪方法。

泥石流次聲信號；經驗模態分解；小波閾值去噪；相關性

0 引言

在泥石流發生過程中，泥土土質的破壞以及巖石的破碎會產生次聲波信號。次聲波擁有衰減慢，穿透力強的聲學信號特點，可以利用該特點實現次聲發生源的準確定位[1]。近幾年來，基于次聲波的泥石流檢測技術已經逐漸成為了泥石流檢測研究方向的一大熱點，隨著次聲探測傳感器性能的提高以及信號處理技術的提升，國內外也開始不斷給出了泥石流次聲研究的新成果，已經可以更加準確的定位泥石流次聲信號的能量分布范圍以及主頻頻段等信號特征[2]。

在泥石流次聲傳感器采集到的信號中，噪聲可能存在于任何一個頻段中[3]。傳統的經驗模態分解(Empirical Mode Decomposition, EMD)聯合小波閾值去噪方法中只能對某一個頻段的噪聲進行硬性去噪[4]，缺乏靈活性并且容易對真實信號造成損傷，由此提出了適合于泥石流信號處理的EMD分解聯合小波閾值去噪方法。

1 方法及步驟

1.1 EMD分解

EMD分解是一種不需要先驗基底的自適應分解方法，信號可以分解為從高頻到低頻排列的若干個固有模態函數(IMF)。

EMD分解算法步驟[5]如下：

1.2 小波閾值去噪

傳統的濾波去噪方法在信號處理過程中都存在著比較嚴重的缺陷，比如，理想低通濾波器會出現嚴重的抖動現象；切比雪夫濾波器在通帶和阻帶范圍內的幅頻特性有波動；巴特沃斯濾波器的幅頻特性是單調下降的[6]。這些經典的信號去噪方法在實現噪聲平滑的過程中必然會引起信號的模糊而降低其清晰度，使其濾波后的熵值增高，從而無法很好地刻畫信號的非平穩性。而小波變換可以對信號進行多尺度細化，適合時變和非平穩信號的時頻分析需求[5]，已成為信號去噪的主要方法之一。

小波閾值去噪的主要步驟如下[7]：

(1) 選取合適的小波以及確定合適的小波分解層數，得到相應的小波分解系數。

(2) 選取合適的閾值，并對上一步得到的小波系數進行閾值處理。

(3) 小波重構，通過得到的新的閾值，對信號進行重構得到去噪后的信號。

2 改進的EMD聯合小波閾值去噪

2.1 方法介紹

傳統的EMD聯合小波閾值去噪方法，只能憑借經驗對一些固定的含噪IMF頻段進行除噪處理，而在泥石流次聲信號中，無法準確確定某一個頻段為噪聲頻段。因此需要通過一種更加靈活的方式去確定含噪IMF頻段，然后再對其進行除噪處理。

針對傳統的EMD聯合小波閾值去噪方法在泥石流信號應用上的不足，對其進行了改進。首先對次聲傳感器測得的含噪泥石流信號進行EMD分解，由于泥石流信號會在含噪泥石流信號中占據較大的比重，因此可以通過計算每一個IMF分量與含噪信號的相關性來確定每一個IMF分量是否包含重要的泥石流信息。每一個IMF分量與原始信號的相關系數可以表示為

最后將去噪后的IMF分量與直接保留的相關性較高的IMF分量進行重構得到去噪后的信號：

2.2 小波基的選擇與閾值的選取

小波閾值去噪過程中最重要的就是小波基的選擇，常用的小波基函數的特征如表1所示。表中CWT為連續傅里葉變換(Continuous Wavelet Transform)，DWT為離散傅里葉變換(Discrete Wavelet Transform)，為消失矩階數。由于泥石流信號是一種存在大量突變的非線性、非平穩信號，因此一般選擇具有正則性、緊支性和消失矩的小波作為小波基。由于Symlet小波函數在正則性、緊支性等方面都要優于其它小波，并且Symlet小波是對Daubechies (dbN)小波的一種改進，不僅具備小波光滑、誤差不容易被察覺、信號重構過程比較光滑的特點，同時具備了較好的正則性。Symlet小波在連續性、支集長度、濾波器長度等方面與Daubechies小波一致[8]，并且Symlet小波具有更好的對稱性，即一定程度上能夠減少對信號進行分析和重構時的相位失真[9]。因此選用Symlet小波作為小波閾值去噪的小波基。

表1 常用小波函數的主要特征

小波閾值去噪過程中閾值的選取分為硬閾值和軟閾值兩種方式[10]，硬閾值的表達式為

軟閾值的表達式為

由于硬閾值的去噪方法在連續性和去噪效果上都沒有軟閾值去噪方法好[11]，因此選擇軟閾值的處理方式。

3 模擬泥石流信號仿真分析

為了驗證EMD分解與小波閾值聯合去噪在泥石流次聲信號去噪中的性能，根據泥石流次聲信號的非線性、非平穩特性，通過仿真實驗獲得的模擬泥石流信號如圖1(a)所示，為該信號添加的真實采集的泥石流原始背景噪聲如圖1(b)所示，加噪后的信號如圖1(c)所示。

圖1 原始信號與加噪信號對比圖

采用小波閾值去噪、傳統EMD分解加小波閾值去噪與本文方法去噪三種方式，分別對加噪后的信號進行去噪處理，得到去噪信號的時域對比圖見圖2，信號的信噪比及與原始信號的相關性見表2。

圖2 不同方法去噪后的時域信號對比圖

表2 三種去噪算法信噪比與相關性表

從表2中數據分析可知，小波閾值去噪和傳統EMD聯合小波閾值去噪得到的去噪信號，在信噪比上要略低于本文方法去噪后得到的去噪信號的信噪比。去噪后的信號與原始信號的相關性也略低于本文方法的結果，即用改進的EMD分解與小波閾值聯合去噪方法得到的信號，在信噪比和原始信號的相似程度上都要好于前兩種去噪方法。

對去噪后的信號進行進一步的頻譜分析，并求得每一種去噪方法去噪后的信號與原始信號的頻譜差，如圖3所示。

從圖3可以明顯看出，傳統的EMD聯合小波閾值去噪的方式得到的頻譜差，雖然在高頻部分擁有很好的去噪效果，但是在低頻部分卻失去了更好的去噪性能，相對而言，改進的方法能智能選擇噪聲頻段，在高頻和低頻段均擁有很好的去噪性能，從而更適于泥石流次聲信號的處理。

4 真實泥石流信號分析處理

4.1 信號數據來源

實驗中用到的數據來源于2014年8月4日云南東川蔣家溝地區的泥石流事件。圖4為原始信號的時域分布圖。

4.2 信號的EMD分解與分量選擇

通過對原始信號進行EMD自適應分解可以得到如圖5所示的20個IMF分量，頻率由IMF1到IMF20依次遞減。

圖4 原始信號時域分布

計算每一個IMF分量與原始信號的相關性，得到相關性分布折線圖如圖6所示。

以模擬泥石流信號作為分析依據，以不同的、值計算去噪重構信號與原始信號的相關性，通過定量分析，發現當、值將數據頻段切割成20%左右的噪聲頻段和10%左右的保留信號頻段時，去噪重構信號與原始信號的相關性最高。

圖6 IMF分量的相關性分布圖

因此，根據當地環境的噪聲特點，將噪聲頻段設置為原始數據的20%，即將舍去IMF分量的閾值下限設置為0.05；將真實泥石流信號頻段設置為原始信號的10%，即將保留IMF分量的閾值上限設置為0.37。

從圖6可以看出，固有模態函數IMF15、IMF17～IMF20與原始信號的相關性低于0.05，完全可以作為無關干擾信號進行舍棄。IMF1和IMF7分量與原始信號的相關性較高，均達到0.37以上，故包含有重要的泥石流信號信息，為防止重要信息丟失，對其進行完全保留，不再進行小波閾值去噪處理。

4.3 小波閾值除噪與信號重構

通過對EMD分解得到的IMF分量進行選擇，舍棄了低于閾值下限的4個IMF分量，完全保留了高于閾值上限的2個IMF分量，對于剩余的13個IMF分量，由于其既包含有相當量的泥石流信號信息，又包含有大量的干擾噪聲，因此需要對其進行去噪處理。

對每一個IMF分量利用Symlet5小波作為小波基進行去噪后得到的信號如圖7所示。

與圖5的各含噪IMF分量相比，經小波閾值去噪后，信號在不同頻段的噪聲含量均有所降低，在保留泥石流重要信息的同時，可以智能地分辨噪聲所在頻段，不再拘泥于消除某一固定頻段的噪聲，因此是針對泥石流次聲信號處理的一種更靈活、可靠的去噪方式。

將小波閾值去噪后的IMF分量與相關性最高的兩個IMF分量進行重構，可以得到如圖8所示的重構信號。

圖8 原始信號和除噪后重構信號對比圖

4.4 重構信號時頻分析

經過EMD分解和小波閾值聯合去噪之后，得到了去除噪聲之后的泥石流信號。為進一步對泥石流信號的頻段信息進行研究，在此對重構信號進行小波時頻分析。

這里選擇cmor小波作為小波時頻分析的小波基，小波帶寬參數為3，中心頻率為3 Hz，得到原始信號的小波時頻圖與重構信號的小波時頻圖如圖9所示。

通過圖9中的時頻分析對比可以發現，經EMD分解與小波閾值聯合去噪后的信號濾除了大量的高頻噪聲，同時去除了一些無關的低頻細節分量，從而可以獲得更加可靠的泥石流信號。對圖9進一步分析可以看出，在泥石流發生時(1 000～1 500 s)，次聲信號能量開始出現明顯的聚集，且能量主要集中在2～6 Hz頻段，在整個泥石流事件中，信號在10 Hz頻段以下都有明顯的能量變化。

圖9 原始信號和重構信號的小波時頻分析圖

5 結論

實驗證明本文提出的方法更適合于泥石流信號的去噪處理，其在濾除了高頻段信號噪聲的同時還能對低頻段的信號進行相應的選擇處理，由于對與原始信號具有高相關性的IMF分量不進行去噪處理，因此又能更好地保護信號細節信息。通過對去噪后的泥石流信號的進一步時頻分析，得到泥石流信號主要頻率范圍為2～6 Hz，且在泥石流發生過程中，整個10 Hz以下頻段都有明顯的能量變化。

[1] 謝濤, 徐小林, 陳洪凱. 泥石流攔擋壩研究現狀及發展趨勢[J]. 中國地質災害與防治學報, 2017, 28(2): 137-145.

XIE Tao, XU Xiaolin, CHEN Hongkai. Research status and development trend of debris flow dam[J]. The Chinese Journal of Geological Hazard and Control, 2017, 28(2): 137-145.

[2] 李朝安, 胡卸文, 王良瑋. 山區鐵路沿線泥石流次聲監測預警方法[J]. 聲學技術, 2012, 31(4): 351-356.

LI Chao’an , HU Xiewen, WANG Liangwei. Infrasound monitoring and early warning of debris flow along montanic railway line[J]. Acoustic Technique, 2012, 31(4): 351-356.

[3] 鄭菲. 次聲波源產生的機理及有限元模擬[D]. 成都: 成都理工大學, 2015.

ZHENG Fei. Infrasound source generation mechanism and finite element simulation[D]. Chengdu: Chengdu University of Technology, 2015.

[4] 趙宇翔, 劉文會, 蘇亮亮, 等. 小波分析在振動信號去噪中的應用研究[J]. 吉林建筑大學學報, 2016, 33(5): 27-30.

ZHAO Yuxiang, LIU Wenhui, SU Liangliang, et al. Application of wavelet analysis in vibration signal denoising[J]. Journal of Jilin Jianzhu University, 2016, 33(5): 27-30.

[5] 徐潔, 王阿明, 鄭小鋒. 基于小波閾值去噪的心電信號分析[J]. 計算機仿真, 2011, 28(12): 260-263.

XU Jie, WANG Aming, ZHENG Xiaofeng. ECG signal analysis based on wavelet threshold denoising, 2011, 28(12): 260-263.

[6] 曾敬楓. 基于MATLAB不同小波基的小波閾值圖像去噪算法[J]. 智能計算機與應用, 2016, 6(4): 75-77.

ZENG Jingfan. Wavelet threshold image denoising algorithm based on different wavelet bases of MATLAB[J]. Intelligent Computer and Applications, 2016, 6(4): 75-77.

[7] 陶偉, 李文堯, 張登, 等. 基于小波變換的TEM信號處理中小波基函數的選擇[J]. 中國錳業, 2016, 34(6): 175-176.

TAO Wei, LI Wenyao, ZHANG Deng, et al. Selection of wavelet basis functions in TEM signal processing based on wavelet transform[J]. China’s Manganese Industry, 2016, 34(6): 175-176.

[8] 王劍平, 張捷. 小波變換在數字圖像處理中的應用[J]. 現代電子技術, 2011, 34(1): 91-94.

WANG Jianping, ZHANG Jie. Application of wavelet transform in digital image processing[J]. Modern Electronics Technique, 2011, 34(1): 91-94.

[9] 劉雪勇. 次聲信號特征提取與分類識別研究[D]. 北京: 中國地質大學(北京), 2015.

LIU Xueyong. Research on infrasound signal feature extraction and classification recognition[D]. Beijing: China University of Geosciences (Beijing), 2015.

[10] LIU D L. Monitoring and recognition of debris flow infrasonic signals[J]. Journal of Mountain Science, 2015, 12(4): 797-815.

[11] 曾憲偉, 趙衛明, 師海闊, 等. 利用小波包變換對地震信號進行時頻分析時小波基函數的選取木[J]. 地震研究, 2010, 33(4): 323-328.

ZENG Xianwei, ZHAO Weiming, SHI Haikuo, et al. Selection of wavelet basis function in time-frequency analysis of seismic signals by wavelet packet transform[J]. Journal of Seismological Research, 2010, 33(4): 323-328.

EMD decomposition and wavelet threshold denoising method for removing noise from debris flow signals

ZHU Feng-jie1, JIAO Rui-li1, TENG Peng-xiao2

(1. School of information and Communication Engineering, Beijing Information Science and Technology University, Beijing 100085, China; 2.The Institute of Acoustics of the Chinese Academy of Sciences,Beijing 100190, China)

Infrasound signals collected by infrasound sensor contain a large number of irrelevant interference signals, which seriously affect the analysis and evaluation of the signals. In view of the characteristics that the noise frequency band can not be accurately identified in the noisy debris flow signal and the shortcoming that the traditional method of empirical mode decomposition (EMD) combined with wavelet threshold denoising can not intelligently distinguish the frequency band where the noise is located, a correlation based method of selecting the noise frequency band is proposed after decomposing the EMD signal. Firstly, the EMD decomposition is used to obtain the intrinsic mode function (IMF) components of the signal, and then the correlation between each IMF component and the original signal is calculated. The frequency band of IMF noise components is selected according to the level of correlation, and then processed by the wavelet threshold de-noising method. Finally, the processed IMF components are reconstructed to get the denoised infrasound signal of debris flow. The simulation results show that this method has the ability to select the noise frequency band intelligently, and is a more suitable denoising method for debris flow signals.

debris flow infrasound signal; empirical mode decomposition (EMD); wavelet threshold denoising; correlation

O425+.3

A

1000-3630(2019)-01-0083-08

10.16300/j.cnki.1000-3630.2019.01.014

2018-01-14;

2018-02-17

北京信息科技大學橫向課題(9161624104)。

朱鳳杰(1991－), 男, 山東日照人, 碩士研究生, 研究方向為次聲信號處理。

朱鳳杰, E-mail: fengjie_zhu @qq.com