基于CEEMDAN-小波包分析的隧道爆破信號(hào)去噪方法*

王海龍，趙 巖，王海軍，彭嬋媛，仝 瀟

（1. 中國礦業(yè)大學(xué)（北京）力學(xué)與建筑工程學(xué)院，北京 100083；2. 河北省土木工程診斷、改造與抗災(zāi)重點(diǎn)實(shí)驗(yàn)室，河北 張家口 075000；3. 北旺建設(shè)集團(tuán)有限公司，河北 承德 067400）

作為鐵路、公路隧道掘進(jìn)最常用的破巖方式，爆破施工帶來的結(jié)構(gòu)振動(dòng)極有可能對(duì)在建隧道及周圍構(gòu)筑物產(chǎn)生損害。目前，對(duì)爆破振動(dòng)的分析研究大多是通過對(duì)爆破振動(dòng)信號(hào)的處理展開的。隧道施工現(xiàn)場環(huán)境復(fù)雜，從現(xiàn)場獲得的原始數(shù)據(jù)往往夾雜著各種各樣的無用信號(hào)（噪聲），尤其隨著隧道掘進(jìn)深度的增加，施工現(xiàn)場各個(gè)工序交叉進(jìn)行，噪聲與有效信號(hào)相互交織，對(duì)爆破信號(hào)處理分析產(chǎn)生極大的不利影響[1]。因此，對(duì)采集的原始隧道爆破振動(dòng)數(shù)據(jù)進(jìn)行降噪預(yù)處理顯得尤為重要。

目前，應(yīng)用于爆破振動(dòng)信號(hào)消除干擾分量的方法主要包括小波閾值方法[2]、小波包閾值方法[3]、最小二乘擬合方法[4]、經(jīng)驗(yàn)?zāi)B(tài)分解（empirical mode decomposition, EMD）方法[5]、集合經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition, EEMD)方法[6]、變分模態(tài)分解（variational mode decomposition,VMD）方法[7]、EMD-小波閾值聯(lián)合方法及EEMD-小波閾值聯(lián)合方法[8-9]等。其中針對(duì)爆破信號(hào)中的高頻噪聲干擾分量，傅里葉分析、小波分析及小波包分析均可應(yīng)用于去除噪聲的領(lǐng)域；同時(shí)，EMD、EEMD 以及其與小波分析、小波包分析的聯(lián)合方法也可以取得更為良好的效果。針對(duì)信號(hào)中存在的低頻趨勢(shì)項(xiàng)，最小二乘法、EMD 分解、EEMD 分解及VMD 分解等方法都能在保持信號(hào)特征的同時(shí)較為準(zhǔn)確地去除趨勢(shì)項(xiàng)。

作為爆破振動(dòng)信號(hào)降噪的常用方法，小波降噪首先通過選擇合適的小波基函數(shù)及分解層數(shù)對(duì)原始信號(hào)進(jìn)行分解，然后選擇合適的閾值函數(shù)對(duì)小波分解系數(shù)進(jìn)行量化，對(duì)于超過閾值的小波分解系數(shù)予以保留，對(duì)于小于閾值的小波分解系數(shù)通過置零操作予以消除；然后，將小波系數(shù)進(jìn)行重構(gòu)即完成了信號(hào)消噪。小波包分析作為對(duì)小波算法的改進(jìn)算法，通過對(duì)小波變換過程中忽略的高頻部分進(jìn)行精細(xì)分解，提高了小波算法的高頻分辨率[10-11]。但是，小波及小波包算法均受制于小波基函數(shù)、閾值類型及閾值函數(shù)的選擇，給復(fù)雜的隧道爆破振動(dòng)信號(hào)的降噪處理帶來困難。

EMD 分解[12]無需像小波分析一樣設(shè)定小波基函數(shù)及閾值函數(shù)，就能將信號(hào)分解為一系列本征模態(tài)分量（IMF）。EEMD 分解針對(duì)EMD 分解中存在的模態(tài)混疊現(xiàn)象[13]，對(duì)原始信號(hào)添加白噪聲信號(hào)進(jìn)行分析。許多學(xué)者通過經(jīng)驗(yàn)?zāi)B(tài)分解去除高頻模態(tài)分量的方式實(shí)現(xiàn)去噪，然而這種“一刀切”式的降噪方式會(huì)在降噪的同時(shí)消除信號(hào)中的主要成分，造成信號(hào)失真。

小波分析與經(jīng)驗(yàn)?zāi)B(tài)類分解的聯(lián)合去噪方法可以綜合兩種方法的優(yōu)點(diǎn)，在消噪的同時(shí)保持原始信號(hào)的完整性和特殊性；EMD-小波閾值降噪、EEMD-小波閾值降噪等方法已經(jīng)應(yīng)用于爆破振動(dòng)信號(hào)的處理中。然而，EMD-小波閾值降噪、EEMD-小波閾值等降噪方法仍存在不少缺陷，例如：雖然EEMD 分解通過在原始信號(hào)中添加高斯白噪聲解決了EMD 分解存在的模態(tài)混疊的問題，但是EEMD 并未對(duì)添加的白噪聲信號(hào)進(jìn)行隔離，造成殘留噪聲由高頻過渡到低頻，影響降噪效果。

隧道爆破振動(dòng)信號(hào)中的噪聲分量主要來自于鉆孔裝藥、初期支護(hù)及運(yùn)輸渣土等施工工序中產(chǎn)生的高頻干擾信號(hào)；同時(shí)，信號(hào)采集過程中，儀器溫度變化造成的零點(diǎn)漂移及測試儀器的安裝松動(dòng)也會(huì)造成采樣信號(hào)產(chǎn)生基線偏離及形狀不規(guī)則。針對(duì)上述隧道爆破信號(hào)處理中存在的噪聲干擾問題，本文引入一種基于CEEMDAN 分解及小波包分析的聯(lián)合去噪方法：首先，對(duì)采集得到的原始信號(hào)進(jìn)行CEEMDAN分解，獲取多個(gè)本征模態(tài)分量（IMF）；然后，利用相關(guān)系數(shù)篩選出包含噪聲的模態(tài)分量[14]，并對(duì)含噪分量進(jìn)行小波包閾值降噪處理；最后，將經(jīng)過處理的模態(tài)分量與未處理的分量重構(gòu)得到純凈的爆破振動(dòng)信號(hào)。通過CEEMDAN 分解來解決添加白噪聲由高頻傳遞到低頻的問題[15-17]，提高高頻成分的頻率分辨率，以期更為精細(xì)地去除噪聲成分。

1 CEEMDAN-小波包分析降噪原理

1.1 CEEMDAN 分解原理闡述

EEMD 與CEEMD 分解算法均是在待分解信號(hào)中加入高斯白噪聲或者成對(duì)正負(fù)高斯白噪聲來消除EMD 分解的模態(tài)混疊問題。然而，這兩種算法均未對(duì)殘留噪聲進(jìn)行隔離，造成添加的白噪聲信號(hào)總是可以從高頻轉(zhuǎn)移到低頻。因此，分解得到的模態(tài)本征分量中總會(huì)殘留一定的白噪聲信號(hào)，影響后續(xù)信號(hào)分析處理。

CEEMDAN 分解從兩個(gè)方面解決了上述問題：(1) 加入經(jīng)EMD 分解后含輔助噪聲的IMF 分量，而不是將高斯白噪聲信號(hào)直接添加在原始信號(hào)中；(2) EEMD 分解和CEEMD 分解是將經(jīng)驗(yàn)?zāi)B(tài)分解后得到的模態(tài)分量進(jìn)行總體平均，CEEMDAN 分解與之不同，則是在得到的第一階IMF 后就進(jìn)行總體平均計(jì)算，得到最終的第一階IMF，然后對(duì)殘余部分重復(fù)進(jìn)行如上操作，這樣便有效地解決了噪聲從高頻到低頻的轉(zhuǎn)移傳遞問題。接下來闡述CEEMDAN 分解原理。

（2）對(duì)產(chǎn)生的N 個(gè)模態(tài)分量進(jìn)行加總平均得到：

（3）計(jì)算去除第1 個(gè)模態(tài)分量后的殘差：

（4）在r1(t)中加入正負(fù)成對(duì)高斯白噪聲得到新信號(hào)，以新信號(hào)為載體進(jìn)行EMD 分解，得到第1 階模態(tài)分量，由此可以得到第2 階模態(tài)分量：

（5）計(jì)算去除第二個(gè)模態(tài)分量的殘差：

（6）重復(fù)上述步驟，直到獲得的殘差信號(hào)為單調(diào)函數(shù)，不能繼續(xù)分解，算法結(jié)束。此時(shí)得到的本征模態(tài)分量數(shù)量為K，則原始信號(hào)z(t)被分解為：

1.2 小波包降噪原理

小波降噪是通過小波變換對(duì)分解得到的小波系數(shù)進(jìn)行分類處理。經(jīng)小波分解后，有用信號(hào)的小波分解系數(shù)較大，噪聲的小波分解系數(shù)較小。因此，通過設(shè)置閾值函數(shù)進(jìn)行去噪處理，保留分解系數(shù)大于閾值的信號(hào)，而通過置零操作消除分解系數(shù)小于閾值的信號(hào)。小波包分析以小波變換為基礎(chǔ)，對(duì)小波變換中未處理的高頻分量再次細(xì)化分解，相比于小波降噪，小波包分析具有更高的頻率分辨率，可以進(jìn)一步消除高頻部分存在的噪聲余量，提高去噪精度。

可以得到小波包分解結(jié)果：

1.3 CEEMDAN-小波包分析聯(lián)合降噪

圖1 為CEEMDAN-小波包分析聯(lián)合降噪流程圖。如圖所示，首先，利用CEEMDAN 算法將爆破振動(dòng)信號(hào)分解，得到多個(gè)本征模態(tài)分量；然后，通過比較各個(gè)模態(tài)分量與原始信號(hào)的相關(guān)系數(shù)篩選出包含噪聲的本征模態(tài)分量，并通過各本征模態(tài)分量的頻譜分布特征及方差貢獻(xiàn)率進(jìn)行校核；最后，利用小波包分析對(duì)篩選出含有噪聲的模態(tài)分量進(jìn)行閾值去噪，將去噪完成的模態(tài)分量與未經(jīng)處理的模態(tài)分量重構(gòu)得到純凈信號(hào)。

圖 1 去噪流程Fig. 1 Flow of de-noising

2 仿真信號(hào)分析

為驗(yàn)證本文引入的CEEMDAN-小波包閾值降噪方法的可行性，且不失一般性，首先對(duì)仿真信號(hào)進(jìn)行降噪分析。近年來，有學(xué)者研究發(fā)現(xiàn)可以近似地利用正余弦函數(shù)模擬爆破振動(dòng)信號(hào)進(jìn)行分析[19]。因此，本文亦利用正余弦函數(shù)疊加得到的仿真信號(hào)進(jìn)行研究。仿真信號(hào)表達(dá)式如下：

針對(duì)上述模擬信號(hào)，利用CEEMDAN 分解將其分解為一系列本征模態(tài)分量，并利用MATLAB 中互相關(guān)函數(shù)的概念計(jì)算得到各個(gè)模態(tài)分量與原始信號(hào)的相關(guān)性，并依據(jù)相關(guān)性絕對(duì)值的大小篩選得到需要進(jìn)行處理的含噪分量。仿真信號(hào)z(t)及CEEMDAN 分解得到的本征模態(tài)分量（C1～C12）見圖2。

圖 2 仿真信號(hào)及模態(tài)分量波形圖Fig. 2 Simulation signal and modal component waveform

圖 3 仿真信號(hào)及降噪處理后的純凈信號(hào)Fig. 3 Simulated signal and pure signal after de-noising

圖 4 草帽山隧道進(jìn)口工區(qū)[20]Fig. 4 Caomaoshan tunnel entrance area[20]

3 工程實(shí)例分析

3.1 工程背景

本文以新建京張高鐵草帽山隧道爆破施工為工程背景進(jìn)行研究。草帽山隧道地形起伏較大，穿越草帽山主脈，穿越第四系上更新統(tǒng)新黃土，粉質(zhì)黏土地等地層，穿越全風(fēng)化凝灰?guī)r，具有中等膨脹性，對(duì)隧道爆破施工造成較大的安全隱患。本文使用的原始數(shù)據(jù)來自于此隧道進(jìn)口段鉆爆法施工產(chǎn)生的爆破振動(dòng)速度信號(hào)。圖4 為隧道進(jìn)口工區(qū)場地布置圖。

3.2 爆破振動(dòng)速度信號(hào)采集

此次爆破信號(hào)采集使用中科測控公司研發(fā)的TC-4850 爆破測振儀。該爆破測振儀可以同時(shí)采集X、Y 及Z 方向的爆破振動(dòng)速度，同時(shí)利用Blasting vibration analysis 的分析軟件可以對(duì)采集到的振速數(shù)據(jù)進(jìn)行初步分析。采集過程中，利用配置好的不銹鋼夾片將速度傳感器固定于隧道初期支護(hù)上。其中X 方向朝向掌子面的掘進(jìn)方向，Y 方向朝向隧道徑向，Z 方向垂直于XY 平面垂直向上。測點(diǎn)布置見圖5。

圖 5 測點(diǎn)布置[21]Fig. 5 Layout of measuring points

圖6 為某次爆破施工采集到的Z 方向爆破振動(dòng)速度原始信號(hào)，從圖6 可知，受爆破施工現(xiàn)場復(fù)雜環(huán)境的影響，爆破振速時(shí)程曲線中夾雜著各種毛刺噪聲信號(hào)。

圖 6 爆破振動(dòng)速度原始信號(hào)Fig. 6 Original signal of blasting vibration speed

3.3 CEEMDAN-小波包閾值降噪實(shí)現(xiàn)

經(jīng)CEEMDAN 分解，原始爆破信號(hào)被分解為14 個(gè)本征模態(tài)分量。依據(jù)EEMD 分解的取值范圍，利用“試錯(cuò)法”得到本次試驗(yàn)中CEEMDAN分解的特征參數(shù)為：正負(fù)高斯白噪聲標(biāo)準(zhǔn)差為0.2，加入噪聲的次數(shù)為100，允許的最大篩選迭代次數(shù)為3 000。

利用式（12）計(jì)算得到每個(gè)模態(tài)分量與原始速度信號(hào)的相關(guān)系數(shù)（ri），計(jì)算結(jié)果見表1。

式中：xi表示各本征模態(tài)分量，y 表示原始爆破信號(hào)。

表 1 本征模態(tài)分量（IMF）的相關(guān)系數(shù)Table 1 Correlation coefficients of modal components (IMF)

從表1 中可以看出，各個(gè)本征模態(tài)分量與原始信號(hào)的相關(guān)系數(shù)差異較大。其中C6、C7、C8、C9的相關(guān)系數(shù)均大于0.2，可以看作是優(yōu)勢(shì)模態(tài)分量；而C1～C5及C10～C14的相關(guān)系數(shù)均小于0.14，初步認(rèn)定C1～C5及C10～C14這兩組模態(tài)分量中含有較多的噪聲信號(hào)。

對(duì)C1～C5及C10～C14兩組分量進(jìn)行快速傅里葉分析得到頻譜見圖7。從圖7 可以看出，C1～C3分量的主頻主要集中在800～1 000 Hz，C4分量的主頻集中在400 ～600 Hz，C5分量的主頻主要分布在300 Hz左右；而C10～C14分量的主頻率均小于5 Hz。由于本文振動(dòng)測試儀的采集頻率主要范圍為5～250 Hz，因此從頻譜分布特征來看，C1～C5及C10～C14分量確實(shí)含有某些干擾信號(hào)。

定義模態(tài)分量平方的算術(shù)平均值與其均值的平方之差為模態(tài)分量的方差（D(j)），通過各個(gè)模態(tài)分量的方差貢獻(xiàn)率（e(j)）[22-23]校核上述選擇的可行性。模態(tài)分量的方差及方差貢獻(xiàn)率計(jì)算如下：

計(jì)算得到各個(gè)模態(tài)分量的方差貢獻(xiàn)率見表2。由表2 可知，C6～C9模態(tài)分量的方差貢獻(xiàn)率較大，均大于8，而剩余模態(tài)分量的方差貢獻(xiàn)率較小，其中貢獻(xiàn)率最大的C10也只有1.59。因此可以看出，由方差貢獻(xiàn)率得到的結(jié)論與相關(guān)系數(shù)篩選及頻譜分析結(jié)果保持一致。綜上，可以認(rèn)定C1～C5和C10～C14兩組分量包含噪聲信號(hào)。

將C1～C5及C10～C14重構(gòu)得到新的分量C15，以C15為新信號(hào)進(jìn)行小波包分析。由于本文中爆破測振儀TC-4 850 設(shè)置的最小工作采樣頻率為5 Hz，根據(jù)采樣定理[24]，信號(hào)采集頻率為5 000 Hz，Nyquist 頻率為2 500 Hz，因此根據(jù)小波包分解算法，利用具有良好的緊支撐性、光滑性及近似對(duì)稱性的“db8”小波基函數(shù)對(duì)新信號(hào)進(jìn)行8 層分解，對(duì)應(yīng)的最低頻帶為0～9.765 625 Hz；然后，利用默認(rèn)的軟閾值函數(shù)進(jìn)行去噪處理。圖8 為小波包降噪處理前后C15信號(hào)時(shí)程，從圖8 可以看出，經(jīng)過小波包降噪重構(gòu)后的C15分量的噪聲信號(hào)基本被消除。

圖 7 C1～C5 及C10～C14 分量頻譜Fig. 7 Spectra of C1-C5 and C10-C14 components

表 2 模態(tài)分量（IMF）的方差貢獻(xiàn)率Table 2 Variance contribution rate of modal component (IMF)

將經(jīng)過降噪處理的C15與優(yōu)勢(shì)分量C6、C7、C8、C9重構(gòu)得到最終的純凈信號(hào)如圖9 所示。由圖9 可知，與原始信號(hào)相比，經(jīng)CEEMDAN-小波包降噪后的純凈信號(hào)在保證局部波形特征及峰值不變的基礎(chǔ)上，基本消除了原始信號(hào)中存在的噪聲分量。

圖 8 小波包降噪處理前后的C15 信號(hào)Fig. 8 signal (C15) before and after wavelet packet noise reduction

圖 9 CEEMDAN-小波包閾值降噪后的信號(hào)Fig. 9 CEEMDAN-wavelet packet threshold signal after noise reduction

3.4 降噪效果評(píng)價(jià)

爆破振動(dòng)信號(hào)降噪效果常利用信噪比（η）、均方根差（σ）進(jìn)行衡量，信噪比越大，均方根差越小，去噪效果越好。信噪比、均方根差的計(jì)算公式[25-26]如下：

為了評(píng)價(jià)本文方法的降噪效果，利用小波包降噪、EMD-小波包聯(lián)合降噪及EEMD-小波包聯(lián)合降噪三種方法分別對(duì)原始爆破信號(hào)進(jìn)行處理。四種方法的降噪指標(biāo)見表3，降噪后的時(shí)程曲線如圖10 所示。

對(duì)比四種方法的客觀降噪指標(biāo)發(fā)現(xiàn)，經(jīng)驗(yàn)?zāi)B(tài)分解與小波包分析聯(lián)合降噪的總體效果均優(yōu)于小波包降噪。同時(shí)，與EMD-小波包聯(lián)合降噪及EEMD-小波包聯(lián)合降噪相比，CEEMDAN-小波包降噪方法得到的信噪比（94.080 2）最大，得到的均方根差（2.400 1×10-5）最小。對(duì)比圖9 及圖10 可以看出，與另外三種方法相比，經(jīng)CEEMDAN-小波包閾值降噪后的純凈信號(hào)在保留原始信號(hào)完整信息的同時(shí)，可以更有效地去除信號(hào)中的噪聲分量，有助于提高后續(xù)對(duì)隧道爆破信號(hào)分析的準(zhǔn)確性。

表 3 去噪效果對(duì)比Table 3 Comparison of noise reduction effects

圖 10 用于對(duì)比的幾種方法的降噪效果Fig. 10 Noise reduction effect of several methods for comparison

4 小波包能量譜分析

以原始爆破信號(hào)及經(jīng)降噪處理后的純凈信號(hào)為研究對(duì)象進(jìn)行小波包分析，選用“db8”小波基函數(shù)分別進(jìn)行8 層分解，并利用MATLAB對(duì)式(17)～(20)進(jìn)行編程，得到各個(gè)子頻帶的能量及信號(hào)總能量，并求得各頻帶能量百分比如圖11 所示。

由圖11 中可知，與原始信號(hào)相比，CEEMDAN-小波包降噪得到的純凈信號(hào)中低頻帶能量百分比有較大幅度的增加；同時(shí)，降噪后信號(hào)的高頻帶能量百分比呈減小的趨勢(shì)。研究發(fā)現(xiàn)[27]，隧道爆破信號(hào)的特征信息主要集中在中低頻率帶。因此，可以得到：經(jīng)過CEEMDAN-小波包降噪處理的隧道爆破信號(hào)在去除噪聲分量的同時(shí)，信號(hào)能量分布向中低頻帶轉(zhuǎn)移，很好地保留了原始爆破信號(hào)中的特征信息，為后續(xù)對(duì)爆破信號(hào)進(jìn)行深入分析奠定了基礎(chǔ)。

圖 11 小波包能量占有百分比Fig. 11 Signal energy distribution before and after noise

5 結(jié) 論

本文以新建京張高鐵草帽山隧道爆破施工為工程背景，引入一種基于CEEMDAN 分解聯(lián)合小波包分析的降噪方法對(duì)實(shí)測爆破振動(dòng)信號(hào)進(jìn)行處理，得到如下結(jié)論：

（1）CEEMDAN 分解可以有效解決EMD 分解中存在的模態(tài)混疊問題及EEMD 分解中存在的白噪聲傳遞問題；與小波去噪相比，小波包去噪可以精細(xì)地去除信號(hào)中的噪聲分量；

（2）與小波包單獨(dú)降噪、EMD-小波包聯(lián)合降噪及EEMD-小波包聯(lián)合降噪方法相比，CEEMDAN-小波包降噪得到的信噪比最大，均方根差最小，去噪效果較好；

（3）通過小波包能量譜分析得到，經(jīng)過CEEMDAN-小波包降噪的信號(hào)可以在消除噪聲分量的同時(shí)有效保留原始信號(hào)中的能量特征信息，為后續(xù)對(duì)爆破信號(hào)進(jìn)行準(zhǔn)確分析奠定了基礎(chǔ)。