朱鳳杰，焦瑞莉，滕鵬曉

?

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

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

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

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

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

0 引言

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

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

1 方法及步驟

1.1 EMD分解

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

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

1.2 小波閾值去噪

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

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

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

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

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

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

2.1 方法介紹

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

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

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

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

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

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

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

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

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

3 模擬泥石流信號(hào)仿真分析

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

圖1 原始信號(hào)與加噪信號(hào)對(duì)比圖

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

圖2 不同方法去噪后的時(shí)域信號(hào)對(duì)比圖

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

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

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

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

4 真實(shí)泥石流信號(hào)分析處理

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

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

4.2 信號(hào)的EMD分解與分量選擇

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

圖4 原始信號(hào)時(shí)域分布

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

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

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

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

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

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

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

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

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

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

圖8 原始信號(hào)和除噪后重構(gòu)信號(hào)對(duì)比圖

4.4 重構(gòu)信號(hào)時(shí)頻分析

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

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

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

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

5 結(jié)論

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

[1] 謝濤, 徐小林, 陳洪凱. 泥石流攔擋壩研究現(xiàn)狀及發(fā)展趨勢(shì)[J]. 中國地質(zhì)災(zāi)害與防治學(xué)報(bào), 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)測(cè)預(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] 趙宇翔, 劉文會(huì), 蘇亮亮, 等. 小波分析在振動(dòng)信號(hào)去噪中的應(yīng)用研究[J]. 吉林建筑大學(xué)學(xué)報(bào), 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] 徐潔, 王阿明, 鄭小鋒. 基于小波閾值去噪的心電信號(hào)分析[J]. 計(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ì)算機(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信號(hào)處理中小波基函數(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] 劉雪勇. 次聲信號(hào)特征提取與分類識(shí)別研究[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)明, 師海闊, 等. 利用小波包變換對(duì)地震信號(hào)進(jìn)行時(shí)頻分析時(shí)小波基函數(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－), 男, 山東日照人, 碩士研究生, 研究方向?yàn)榇温曅盘?hào)處理。

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