基于CEEMD-MPE-IMPROVEDWTD的爆炸振動信號噪聲抑制方法

童曉玲, 曾斐宇, 江 軻, 梁 磊

(武漢理工大學 機電工程學院, 武漢 430070)

爆破振動是爆炸后炸藥的部分能量以彈性波的形式向外傳播,使巖體發生振動的現象。爆破振動分析是研究爆破振動損傷控制技術的重要基礎之一[1]。在防護工程試驗中,為了準確定位爆炸源,需要分析爆炸產生的振動信號。但由于爆破現場環境的復雜性和現場設備的誤差,爆破測量信號中含有大量的噪聲成分,給信號分析帶來了困難。為了確定采集到信號的特征,獲得真實的信號,需要對被測信號進行降噪處理。

爆破振動信號是典型的非平穩、隨機信號,可用經驗模態分解類算法進行去噪處理。khaldi等[2]采用EMD(empirical mode decomposition)算法對語音信號進行去噪。Kopsinis等[3]設計了一種基于經驗模態分解并結合頻率分析的新去噪方案,該方案在原始信號能量水平較低時,去噪效果較好。但EMD分解會出現明顯的端點效應和模態混疊,需要對算法進行改進。Chen等[4]利用WT(wavelet transform)和EEMD(ensemble empirical mode decomposition)對電磁聲信號去噪,發現EEMD的適應性更強,優于WT。雖然EEMD分解可以改善模態混疊現象,但去噪效果仍不理想。CEEMD(complementary ensemble empirical mode decomposition)算法是EMD的一種改進算法,既保留了EMD處理非平穩信號的優點,又能有效克服EMD的模態混疊問題。但CEEMD去噪方法會直接丟棄高頻IMF分量,直接重構剩余的低頻IMF(intrinsic mode function)分量,從而導致高頻[5]分量高頻噪聲抑制不完全或有效信號丟失。

小波算法同樣適用于處理非平穩信號[6-11],Jiang等[12]提出了小波模量最大算法,通過試驗證明小波變換是一種有效的信號去噪算法。Zhang等[13]將小波閾值去噪方法應用于爆破振動信號的預處理,并對四種去噪閾值進行比較。小波閾值函數分為硬閾值函數和軟閾值函數。傳統的硬閾值和軟閾值去噪效果都存在不足。硬閾值會在某些點產生不連續,導致較大的均方誤差和波形振蕩,而軟閾值容易造成高頻信息的部分丟失。改進的小波閾值函數可以在一定程度上消除前兩個閾值函數的缺陷。面對單獨使用模態分解方法和小波方法所遇到的問題,兩種方法可以結合起來抑制信號[14-17]的噪聲。與單一CEEMD方法和分離小波變換方法相比,該組合方法對微震信號去噪效果更佳。

本文提出了基于MPE的CEEMD和小波閾值相結合的噪聲抑制方法,對爆破振動信號進行去噪處理。首先對被測信號進行CEEMD分解,得到各階IMF分量;其次,得到各階IMF的MPE值。根據MPE閾值確定高頻分量和低頻分量。最后,利用改進的WTD對高頻分量進行處理,同未去噪的IMF分量進行重構,完成信號去噪,并應用于震源定位的預處理。

1 CEEMD-MPE-WTD算法

1.1 CEEMD原理

CEEMD算法是EMD的一種改進算法,它既保留了EMD在處理非平穩信號方面的優點,又能有效克服EMD的模態混疊問題。CEEMD的主要步驟如下：

(1) 向原始信號中加入n組正負成對的白噪聲,生成2n組信號。

(2) 對2n組信號分別進行EMD分解,得到2n組IMF分量cij。其中i表示添加輔助噪聲的第i個信號,j表示對該信號進行EMD分解得到的第j個IMF分量。

(3) 對2n組IMF分量對應尺度取平均,得到CEEMD的對應分解結果

(1)

1.2 多尺度排列熵原理

多尺度排列熵(multiscale permutation entropy, MPE)是計算時間信號經過不同尺度粗粒度運算后的排列熵,可以分析信號的隨機性和復雜性。步驟如下：

(1) 對時間序列X={x(i),i=1,2,…,N}進行粗粒化處理得到粗粒化后的序列

(2)

(3)

式中：m為嵌入維數;τ為延遲因子。

(4)

(5)

(4) 對尺度因子為s的序列的排列熵進行計算并歸一化處理

(6)

通過MPE方法分析CEEMD分解得到的各階IMF的復雜性,判斷MPE值大于0.6的IMF為高頻分量,MPE值小于0.6的為低頻分量[18]。

1.3 小波閾值原理

小波閾值算法廣泛應用于信號去噪和圖像去噪中。其原理是對信號進行小波分解,得到一系列小波系數,判斷大于閾值的小波系數為有用信號,而確定小于閾值的小波系數為噪聲,可以去除[19],從而達到抑制噪聲的目的。小波閾值算法的具體步驟如下：

(1) 對信號進行小波分解。在分解時,需要選擇合適的小波基和分解層數。常用的小波基有dB小波和sym小波。經過多次試驗選擇了sym8小波。選擇3層作為分解層的數量。

(2) 閾值處理,選擇固定閾值。

(3) 閾值函數分為硬閾值函數和軟閾值函數。但軟硬閾值法本身也存在一些缺陷。硬閾值函數在閾值處不連續,對于每一層細節分量,經過閾值處理后得到的小波系數都會產生額外的振蕩。軟閾值函數是連續的,但經過閾值處理后的小波系數與小波系數之間存在恒定的偏差。為了提高小波閾值去噪(wavelet threshold denoising, WTD)效果,采用軟硬閾值折衷函數作為閾值函數。

軟硬閾值折衷函數的表達式是

(7)

式中：ω為閾值處理前的小波系數;ωthr為閾值處理后的小波系數;thr為小波閾值;α為調整參數,0≤α≤1。

(4) 信號重構。重構閾值處理后的近似分量和細節分量,得到小波閾值去噪后的信號。

改進的小波閾值函數,如圖1所示。

圖1 改進的小波閾值函數

1.4 信號去噪算法

本文提出的去噪算法是將CEEMD、MPE和改進的WTD算法相結合。算法流程圖如圖2所示。噪聲抑制方法的實現步驟如下：

(1) 對被測信號進行CEEMD分解,得到各階IMF。

(2) 計算IMF的MPE值。當IMF的MPE大于0.6時,可以認為是一個有噪聲的高頻分量,需要進行去噪處理。當IMF的MPE小于0.6時,可以認為是低頻分量。

(3) 采用改進的WTD算法對高頻分量進行去噪,并與低頻分量重構。

2 仿真信號分析

由于由正弦和余弦信號組成的疊加信號與爆破振動信號[20]具有相似的特性,我們利用疊加信號來驗證所選方法的性能。仿真信號可以表示為

圖2 基于CEEMD、MPE和改進WTD的噪聲抑制流程圖

(8)

式中,y4為高斯白噪聲。

分別使用EEMD-MPE方法、CEEMD-MPE方法、CEEMD-MPE-WTD方法和CEEMD-MPE-IMPROVED WTD對模擬信號y(t)進行去噪處理,采用信噪比(signal-to-noise ratio, SNR)、均方根誤差(root mean square error, RMSE)和皮爾遜相關系數(Pearson correlation coefficient, PCC)作為抑制噪聲的指標,處理結果如表1所示。

表1 四種不同方法的噪聲抑制效果比較

由表1可知,CEEMD-MPE-IMPROVED WTD處理得到的信號信噪比最高,均方根誤差最小。通過計算四種去噪方法下純凈信號和去噪信號的PCC值,發現該方法去噪方法的PCC值最接近1,表明該方法中的去噪信號與純信號最相關,能夠最有效地保留原始信號的有效信息。

結果表明,CEEMD-MPE-IMPROVED WTD算法效果最佳,為后續對實測隧道爆破振動信號的降噪處理提供了理論依據。

3 工程實例

CEEMD-MPE-IMPROVED WTD算法研究是基于江蘇省南京市的隧道爆炸試驗。土地為相對平坦的黃土層,含一定程度的砂土,相對干燥堅硬,表面有少量礫石。

在實際爆炸試驗中,以隧道為爆炸主體,在隧道口放置炸藥進行爆炸,在隧道外鋪設光纜,采用分布式光纖聲傳感系統(distributed acoustic sensing, DAS)結合光纜采集爆炸振動信號。DAS是一種最先進的光纖聲學檢測技術,它使用光纖作為傳感器來檢測振動信號。通過檢測光纖中瑞利后向散射的變化,幾乎可以在光纖上的任何點檢測振動信息[21]。現場平面測試圖如圖3所示。現場試驗的場地存在約為1 m的垂直高度差,測試場地長度約為100 m,寬度約為50 m。本次爆炸試驗中使用的DAS系統的空間分辨率為1 m,采樣頻率為20 000 Hz,最大測量長度為20 km。

圖3 測試系統整體示意圖

DAS系統采集的震源時空數據可以組成矩陣。假設有L個傳感節點,傳感節點個數由DAS系統的空間分辨率和光纖長度決定,DAS系統采集時間序列長度為N=fs×t;其中t為時間,t=1 s;fs為采樣率,fs=20 kHz。由此可得N行L列的DAS信號矩陣

(9)

光纜響應圖即為原始數據X所繪制瀑布圖,其中橫軸代表光纜的位置信息,單位為m。縱軸代表時間信息,單位為s。光纜響應圖的顏色代表信號的強度,顏色越深,代表該該位置該時刻的光相位變化越明顯。其中某一位置的信號波形如圖4所示。光纜響應圖如圖5所示。

圖4 爆破振動信號

圖5 光纜響應圖

在對采集到的爆炸信號進行去噪時,首先對一個信號進行CEEMD分解,得到各階IMF,并繪制其頻譜圖,結果如圖6所示,我們可以看到通過測量信號的CEEMD分解獲得的各階IMF分量及其頻譜。可以看出,一階和二階IMF包含頻率約3 000 Hz的設備共振噪聲、環境噪聲和有用信息。如果直接丟棄一階和二階IMF分量,信號中的有用信息將丟失,信號將失真。因此,為了獲得真實信號,有必要在盡可能抑制噪聲分量的同時,將有用信息保留在高階IMF分量中。

在信號通過CEEMD方法分解得到各階IMF后,我們使用MPE方法來判斷各階IMF的類型。IMF1～IMF2的MPE值大于0.6,這被視為噪聲分量,需要使用改進的WTD進行去噪;IMF3～IMF11的MPE值小于0.6,將其視為信號分量,并在去噪后用噪聲分量重構。

為了說明改進WTD對高頻分量中噪聲分量的抑制效果,將CEEMD分解得到的IMF1分量與經過改進小波閾值處理后的IMF1分量進行了比較,繪制了它們的波形和頻譜圖,如圖7所示。

如圖7所示,經過改進的小波閾值去噪后,時域波形的噪聲分量顯著減少,頻譜IMF1的噪聲分量的幅度從0.248降低到0.147,表明IMF1經過此方法處理后確實可以抑制噪聲分量。然而,信號的有用成分也存在一定程度的降低。

圖6 CEEMD分解結果及其對應的譜

(a) CEEMD分解得到IMF1

(b) 改進小波閾值處理得到的IMF1

為了分析該去噪算法中有效成分的損失程度。我們將信噪比為25 dB的高斯白噪聲添加到采集的信號中,并分別通過EEMD-MPE、CEEMD-MPE、CEEMD-MPE-WTD和CEEMD-MPE-IMPROVED WTD算法進行去噪。每個信號的波形圖和AOK時頻分析圖如圖8所示。X是主頻能量,Y是主頻。

(a) 原始信號

(b) EEMD-MPE 方法

(c) CEEMD-MPE 方法

(d) CEEMD-MPE-WTD 方法

(e) CEEMD-MPE-IMPROVED WTD方法

圖8顯示了原始信號和通過各種方法去噪的信號的波形圖和頻譜圖。信號的采樣時間為1 s。從圖8(a)中的頻譜圖可以看出,原始信號的主頻率能量為64.760 4,并且在3 053 Hz的頻率處存在峰值能量為3.84的噪聲分量。

從信號波形可以看出,CEEMD-MPE方法由于直接丟棄高頻IMF分量,3 053 Hz頻率處的噪聲分量消失,信號波形最平滑。似乎CEEMD-MPE方法后的信號去噪效果最好。但將原始信號的主頻能量與通過各種方法去噪的信號的主頻能量進行比較時,如表2所示。本文提出方法的主頻能量僅減少了0.000 2,其他方法獲得的主頻能量損失ΔX=0.000 5～0.053 8,與本文提出方法相比能量損耗相差很多倍。且使用本文提出的去噪方法,3 053 Hz處信號噪聲分量的能量從3.84降低到1.53。由此可以判斷出,CEEMD-MPE-IMPROVED WTD算法相比于其他算法,不僅能很好地抑制信號中的噪聲,而且在很大程度上保留了信號的有效成分。

表2 各方法的信號主頻能量X和主頻Y的結果比較

為了進一步證明CEEMD-MPE-IMPROVED WTD算法的優越性,用上述各種方法對信號進行去噪處理,并計算信號的SNR和RMSE。

表3為EEMD-MPE、CEEMD-MPE、CEEMD-MPE-WTD和CEEMD-MPE-IMPROVED WTD算法的噪聲抑制結果指標。CEEMD-MPE-IMPROVED WTD算法相比于其他方法信噪比明顯提高,均方根誤差顯著降低。證明了噪聲抑制方法的合理性和有效性。

表3 四種方法的噪聲抑制結果比較

為了進一步驗證該方法的適用性,我們隨機選擇了采集的40個爆炸信號,添加25 dB的高斯白噪聲,采用該方法應用于這些數據的去噪,求取SNR和RMSE并作圖。

圖9顯示了采用四種不同方法對隨機選擇的40個振動信號噪聲抑制后的SNR值。從圖9可以看出,CEEMD-MPE-IMPROVED WTD算法獲得的SNR值明顯高于其他算法。圖10顯示了采用四種不同方法對隨機選擇的40個振動信號噪聲抑制后的RMSE值,CEEMD-MPE-IMPROVED WTD算法得到的RMSE值最低,證明了本文提出算法的優勢。

圖9 四種方法的信噪比值的折線圖

圖10 四種方法的均方根誤差值的折線圖

為了驗證本文提出的去噪方法在震源定位中的作用,分別利用去噪前后的數據進行震源定位前的信號到時提取。已知DAS系統和光纜起始通道在(0,0)坐標處,真實震源位置為(60,25)。

采用Geiger定位方法,首選確定一個目標函數,利用接收點到估計震源和到真實震源的走時差最小確定,并采用優化算法使目標函數值最小[22]。

目標函數如下所示

(10)

式中：Ti為光纜第i個傳感點接收到振動信號的時刻,即信號到達時間,利用赤池信息準則(Akaike information criterion, AIC)算法求得;t0為炸藥爆炸時刻;v為振動波在地下的傳播速度;(xi,yi,zi)為光纜第i個傳感點的坐標;(x,y,z)為假定震源位置;n為光纜上傳感點的個數,且n=300。

由于現場試驗的場地并非完全水平,存在約為1 m的垂直高度差。但是現場測試時,測試場地長度約為100 m,寬度約為50 m,由于高度引起的誤差可忽略。所以本稿件中目標函數可簡化為

(11)

利用遺傳算法優化得到使誤差值ε最小的假定震源坐標(x0,y0),認為遺傳算法優化得到的震源坐標(x0,y0)即為真實震源位置。分別利用去噪前后信號進行10次定位,得到光纜的實際響應圖與定位結果圖。

如圖11所示為原始信號的到時提取結果以及震源定位結果圖。 圖11(a)中每個通道上的點代表算法捕捉的每個信號的起始時間。已知震源真實坐標為(60,25),由于噪聲成分的干擾使得地震波初至時間判斷存在較大誤差,據此計算得到的震源位置與真實震源位置差距較大,定位結果如圖11(b)所示。

(a) 原始信號的光纜響應圖

(b) 原始信號的定位結果圖

利用去噪后的信號重復定位過程,如圖12所示,到時提取的時間能較為準確的反映地震波到達的起始時間,通過去噪后得到的到時計算的震源位置為(64,27),與震源真實位置(60,25)誤差為5 m左右。

(a) 去噪后信號的光纜響應圖

(b) 去噪后信號的定位結果圖

為了進一步確認論文提出方法優越性,分別采用原始數據、EEMD-MPE、CEEMD-MPE、CEEMD-MPE-WTD以及本文提出算法去噪后,利用相同的定位算法求出定位結果以及定位誤差,并繪制表格如表4所示。

表4 不同算法的定位結果對比

根據定位結果可以得出結論,該去噪算法對爆炸振動信號有較好的去噪效果,可用于震源定位前對信號的預處理。

4 結 論

本文提出了一種基于互補集成經驗模態分解(CEEMD)、多尺度置換熵(MPE)和改進小波閾值去噪(WTD)的融合降噪方法對爆炸振動信號進行降噪處理。結論如下：

(1) CEEMD,MPE,以及改進WTD的融合算法在處理隧道的爆破振動信號時效果較好,可有效地去除隧道爆破振動信號的高頻噪聲成分。

(2) 基于AOK時頻分析技術,分析CEEMD-MPE-IMPROVED WTD算法得到的信號波形及主頻能量可知,該算法在有效去除噪聲成分的同時,由于保留了高頻IMF分量,信號的主成分幾乎沒有損失。為后續分析爆源位置及炸藥威力奠定了基礎。

(3) 采用EEMD-MPE、CEEMD-MPE、CEEMD-MPE-WTD和CEEMD-MPE-IMPROVED WTD算法分別處理爆炸振動信號時,CEEMD-MPE-IMPROVED WTD算法處理后的信號的去噪效果最優。對隨機40個信號進行重復試驗,得出相同結論。驗證了該算法的有效性和適用性,適合用于隧道爆炸振動信號的去噪處理。

(4) 將CEEMD-MPE-IMPROVED WTD算法處理后的信號用于震源定位應用中,定位誤差不超過5 m,遠優于原始信號定位結果,說明該方法適合用于爆炸信號定位前的信號預處理過程。