朱鳳杰，焦瑞莉，滕鵬曉

?

適于泥石流除噪的EMD聯(lián)合小波閾值除噪方法

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

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

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

泥石流次聲信號；經(jīng)驗?zāi)B(tài)分解；小波閾值去噪；相關(guān)性

0 引言

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

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

1 方法及步驟

1.1 EMD分解

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

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

1.2 小波閾值去噪

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

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

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

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

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

2 改進(jìn)的EMD聯(lián)合小波閾值去噪

2.1 方法介紹

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

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

最后將去噪后的IMF分量與直接保留的相關(guān)性較高的IMF分量進(jìn)行重構(gòu)得到去噪后的信號：

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

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

表1 常用小波函數(shù)的主要特征

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

軟閾值的表達(dá)式為

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

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

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

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

采用小波閾值去噪、傳統(tǒng)EMD分解加小波閾值去噪與本文方法去噪三種方式，分別對加噪后的信號進(jìn)行去噪處理，得到去噪信號的時域?qū)Ρ葓D見圖2，信號的信噪比及與原始信號的相關(guān)性見表2。

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

表2 三種去噪算法信噪比與相關(guān)性表

從表2中數(shù)據(jù)分析可知，小波閾值去噪和傳統(tǒng)EMD聯(lián)合小波閾值去噪得到的去噪信號，在信噪比上要略低于本文方法去噪后得到的去噪信號的信噪比。去噪后的信號與原始信號的相關(guān)性也略低于本文方法的結(jié)果，即用改進(jìn)的EMD分解與小波閾值聯(lián)合去噪方法得到的信號，在信噪比和原始信號的相似程度上都要好于前兩種去噪方法。

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

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

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

4.1 信號數(shù)據(jù)來源

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

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

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

圖4 原始信號時域分布

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

以模擬泥石流信號作為分析依據(jù)，以不同的、值計算去噪重構(gòu)信號與原始信號的相關(guān)性，通過定量分析，發(fā)現(xiàn)當(dāng)、值將數(shù)據(jù)頻段切割成20%左右的噪聲頻段和10%左右的保留信號頻段時，去噪重構(gòu)信號與原始信號的相關(guān)性最高。

圖6 IMF分量的相關(guān)性分布圖

因此，根據(jù)當(dāng)?shù)丨h(huán)境的噪聲特點，將噪聲頻段設(shè)置為原始數(shù)據(jù)的20%，即將舍去IMF分量的閾值下限設(shè)置為0.05；將真實泥石流信號頻段設(shè)置為原始信號的10%，即將保留IMF分量的閾值上限設(shè)置為0.37。

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

4.3 小波閾值除噪與信號重構(gòu)

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

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

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

將小波閾值去噪后的IMF分量與相關(guān)性最高的兩個IMF分量進(jìn)行重構(gòu)，可以得到如圖8所示的重構(gòu)信號。

圖8 原始信號和除噪后重構(gòu)信號對比圖

4.4 重構(gòu)信號時頻分析

經(jīng)過EMD分解和小波閾值聯(lián)合去噪之后，得到了去除噪聲之后的泥石流信號。為進(jìn)一步對泥石流信號的頻段信息進(jìn)行研究，在此對重構(gòu)信號進(jìn)行小波時頻分析。

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

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

圖9 原始信號和重構(gòu)信號的小波時頻分析圖

5 結(jié)論

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

[1] 謝濤, 徐小林, 陳洪凱. 泥石流攔擋壩研究現(xiàn)狀及發(fā)展趨勢[J]. 中國地質(zhì)災(zāi)害與防治學(xué)報, 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] 李朝安, 胡卸文, 王良瑋. 山區(qū)鐵路沿線泥石流次聲監(jiān)測預(yù)警方法[J]. 聲學(xué)技術(shù), 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] 鄭菲. 次聲波源產(chǎn)生的機(jī)理及有限元模擬[D]. 成都: 成都理工大學(xué), 2015.

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

[4] 趙宇翔, 劉文會, 蘇亮亮, 等. 小波分析在振動信號去噪中的應(yīng)用研究[J]. 吉林建筑大學(xué)學(xué)報, 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]. 計算機(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]. 智能計算機(jī)與應(yīng)用, 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信號處理中小波基函數(shù)的選擇[J]. 中國錳業(yè), 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] 王劍平, 張捷. 小波變換在數(shù)字圖像處理中的應(yīng)用[J]. 現(xiàn)代電子技術(shù), 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]. 北京: 中國地質(zhì)大學(xué)(北京), 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] 曾憲偉, 趙衛(wèi)明, 師海闊, 等. 利用小波包變換對地震信號進(jìn)行時頻分析時小波基函數(shù)的選取木[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

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

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

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