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

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

（1. 中國礦業大學（北京）力學與建筑工程學院，北京 100083；2. 河北省土木工程診斷、改造與抗災重點實驗室，河北 張家口 075000；3. 北旺建設集團有限公司，河北 承德 067400）

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

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

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

EMD 分解[12]無需像小波分析一樣設定小波基函數及閾值函數，就能將信號分解為一系列本征模態分量（IMF）。EEMD 分解針對EMD 分解中存在的模態混疊現象[13]，對原始信號添加白噪聲信號進行分析。許多學者通過經驗模態分解去除高頻模態分量的方式實現去噪，然而這種“一刀切”式的降噪方式會在降噪的同時消除信號中的主要成分，造成信號失真。

小波分析與經驗模態類分解的聯合去噪方法可以綜合兩種方法的優點，在消噪的同時保持原始信號的完整性和特殊性；EMD-小波閾值降噪、EEMD-小波閾值降噪等方法已經應用于爆破振動信號的處理中。然而，EMD-小波閾值降噪、EEMD-小波閾值等降噪方法仍存在不少缺陷，例如：雖然EEMD 分解通過在原始信號中添加高斯白噪聲解決了EMD 分解存在的模態混疊的問題，但是EEMD 并未對添加的白噪聲信號進行隔離，造成殘留噪聲由高頻過渡到低頻，影響降噪效果。

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

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

1.1 CEEMDAN 分解原理闡述

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

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

（2）對產生的N 個模態分量進行加總平均得到：

（3）計算去除第1 個模態分量后的殘差：

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

（5）計算去除第二個模態分量的殘差：

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

1.2 小波包降噪原理

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

可以得到小波包分解結果：

1.3 CEEMDAN-小波包分析聯合降噪

圖1 為CEEMDAN-小波包分析聯合降噪流程圖。如圖所示，首先，利用CEEMDAN 算法將爆破振動信號分解，得到多個本征模態分量；然后，通過比較各個模態分量與原始信號的相關系數篩選出包含噪聲的本征模態分量，并通過各本征模態分量的頻譜分布特征及方差貢獻率進行校核；最后，利用小波包分析對篩選出含有噪聲的模態分量進行閾值去噪，將去噪完成的模態分量與未經處理的模態分量重構得到純凈信號。

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

2 仿真信號分析

為驗證本文引入的CEEMDAN-小波包閾值降噪方法的可行性，且不失一般性，首先對仿真信號進行降噪分析。近年來，有學者研究發現可以近似地利用正余弦函數模擬爆破振動信號進行分析[19]。因此，本文亦利用正余弦函數疊加得到的仿真信號進行研究。仿真信號表達式如下：

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

圖 2 仿真信號及模態分量波形圖Fig. 2 Simulation signal and modal component waveform

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

圖 4 草帽山隧道進口工區[20]Fig. 4 Caomaoshan tunnel entrance area[20]

3 工程實例分析

3.1 工程背景

本文以新建京張高鐵草帽山隧道爆破施工為工程背景進行研究。草帽山隧道地形起伏較大，穿越草帽山主脈，穿越第四系上更新統新黃土，粉質黏土地等地層，穿越全風化凝灰巖，具有中等膨脹性，對隧道爆破施工造成較大的安全隱患。本文使用的原始數據來自于此隧道進口段鉆爆法施工產生的爆破振動速度信號。圖4 為隧道進口工區場地布置圖。

3.2 爆破振動速度信號采集

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

圖 5 測點布置[21]Fig. 5 Layout of measuring points

圖6 為某次爆破施工采集到的Z 方向爆破振動速度原始信號，從圖6 可知，受爆破施工現場復雜環境的影響，爆破振速時程曲線中夾雜著各種毛刺噪聲信號。

圖 6 爆破振動速度原始信號Fig. 6 Original signal of blasting vibration speed

3.3 CEEMDAN-小波包閾值降噪實現

經CEEMDAN 分解，原始爆破信號被分解為14 個本征模態分量。依據EEMD 分解的取值范圍，利用“試錯法”得到本次試驗中CEEMDAN分解的特征參數為：正負高斯白噪聲標準差為0.2，加入噪聲的次數為100，允許的最大篩選迭代次數為3 000。

利用式（12）計算得到每個模態分量與原始速度信號的相關系數（ri），計算結果見表1。

式中：xi表示各本征模態分量，y 表示原始爆破信號。

表 1 本征模態分量（IMF）的相關系數Table 1 Correlation coefficients of modal components (IMF)

從表1 中可以看出，各個本征模態分量與原始信號的相關系數差異較大。其中C6、C7、C8、C9的相關系數均大于0.2，可以看作是優勢模態分量；而C1～C5及C10～C14的相關系數均小于0.14，初步認定C1～C5及C10～C14這兩組模態分量中含有較多的噪聲信號。

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

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

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

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

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

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

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

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

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

3.4 降噪效果評價

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

為了評價本文方法的降噪效果，利用小波包降噪、EMD-小波包聯合降噪及EEMD-小波包聯合降噪三種方法分別對原始爆破信號進行處理。四種方法的降噪指標見表3，降噪后的時程曲線如圖10 所示。

對比四種方法的客觀降噪指標發現，經驗模態分解與小波包分析聯合降噪的總體效果均優于小波包降噪。同時，與EMD-小波包聯合降噪及EEMD-小波包聯合降噪相比，CEEMDAN-小波包降噪方法得到的信噪比（94.080 2）最大，得到的均方根差（2.400 1×10-5）最小。對比圖9 及圖10 可以看出，與另外三種方法相比，經CEEMDAN-小波包閾值降噪后的純凈信號在保留原始信號完整信息的同時，可以更有效地去除信號中的噪聲分量，有助于提高后續對隧道爆破信號分析的準確性。

表 3 去噪效果對比Table 3 Comparison of noise reduction effects

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

4 小波包能量譜分析

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

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

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

5 結 論

本文以新建京張高鐵草帽山隧道爆破施工為工程背景，引入一種基于CEEMDAN 分解聯合小波包分析的降噪方法對實測爆破振動信號進行處理，得到如下結論：

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

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

（3）通過小波包能量譜分析得到，經過CEEMDAN-小波包降噪的信號可以在消除噪聲分量的同時有效保留原始信號中的能量特征信息，為后續對爆破信號進行準確分析奠定了基礎。