童曉玲, 曾斐宇, 江 軻, 梁 磊
(武漢理工大學(xué) 機(jī)電工程學(xué)院, 武漢 430070)
爆破振動(dòng)是爆炸后炸藥的部分能量以彈性波的形式向外傳播,使巖體發(fā)生振動(dòng)的現(xiàn)象。爆破振動(dòng)分析是研究爆破振動(dòng)損傷控制技術(shù)的重要基礎(chǔ)之一[1]。在防護(hù)工程試驗(yàn)中,為了準(zhǔn)確定位爆炸源,需要分析爆炸產(chǎn)生的振動(dòng)信號(hào)。但由于爆破現(xiàn)場(chǎng)環(huán)境的復(fù)雜性和現(xiàn)場(chǎng)設(shè)備的誤差,爆破測(cè)量信號(hào)中含有大量的噪聲成分,給信號(hào)分析帶來(lái)了困難。為了確定采集到信號(hào)的特征,獲得真實(shí)的信號(hào),需要對(duì)被測(cè)信號(hào)進(jìn)行降噪處理。
爆破振動(dòng)信號(hào)是典型的非平穩(wěn)、隨機(jī)信號(hào),可用經(jīng)驗(yàn)?zāi)B(tài)分解類(lèi)算法進(jìn)行去噪處理。khaldi等[2]采用EMD(empirical mode decomposition)算法對(duì)語(yǔ)音信號(hào)進(jìn)行去噪。Kopsinis等[3]設(shè)計(jì)了一種基于經(jīng)驗(yàn)?zāi)B(tài)分解并結(jié)合頻率分析的新去噪方案,該方案在原始信號(hào)能量水平較低時(shí),去噪效果較好。但EMD分解會(huì)出現(xiàn)明顯的端點(diǎn)效應(yīng)和模態(tài)混疊,需要對(duì)算法進(jìn)行改進(jìn)。Chen等[4]利用WT(wavelet transform)和EEMD(ensemble empirical mode decomposition)對(duì)電磁聲信號(hào)去噪,發(fā)現(xiàn)EEMD的適應(yīng)性更強(qiáng),優(yōu)于WT。雖然EEMD分解可以改善模態(tài)混疊現(xiàn)象,但去噪效果仍不理想。CEEMD(complementary ensemble empirical mode decomposition)算法是EMD的一種改進(jìn)算法,既保留了EMD處理非平穩(wěn)信號(hào)的優(yōu)點(diǎn),又能有效克服EMD的模態(tài)混疊問(wèn)題。但CEEMD去噪方法會(huì)直接丟棄高頻IMF分量,直接重構(gòu)剩余的低頻IMF(intrinsic mode function)分量,從而導(dǎo)致高頻[5]分量高頻噪聲抑制不完全或有效信號(hào)丟失。
小波算法同樣適用于處理非平穩(wěn)信號(hào)[6-11],Jiang等[12]提出了小波模量最大算法,通過(guò)試驗(yàn)證明小波變換是一種有效的信號(hào)去噪算法。Zhang等[13]將小波閾值去噪方法應(yīng)用于爆破振動(dòng)信號(hào)的預(yù)處理,并對(duì)四種去噪閾值進(jìn)行比較。小波閾值函數(shù)分為硬閾值函數(shù)和軟閾值函數(shù)。傳統(tǒng)的硬閾值和軟閾值去噪效果都存在不足。硬閾值會(huì)在某些點(diǎn)產(chǎn)生不連續(xù),導(dǎo)致較大的均方誤差和波形振蕩,而軟閾值容易造成高頻信息的部分丟失。改進(jìn)的小波閾值函數(shù)可以在一定程度上消除前兩個(gè)閾值函數(shù)的缺陷。面對(duì)單獨(dú)使用模態(tài)分解方法和小波方法所遇到的問(wèn)題,兩種方法可以結(jié)合起來(lái)抑制信號(hào)[14-17]的噪聲。與單一CEEMD方法和分離小波變換方法相比,該組合方法對(duì)微震信號(hào)去噪效果更佳。
本文提出了基于MPE的CEEMD和小波閾值相結(jié)合的噪聲抑制方法,對(duì)爆破振動(dòng)信號(hào)進(jìn)行去噪處理。首先對(duì)被測(cè)信號(hào)進(jìn)行CEEMD分解,得到各階IMF分量;其次,得到各階IMF的MPE值。根據(jù)MPE閾值確定高頻分量和低頻分量。最后,利用改進(jìn)的WTD對(duì)高頻分量進(jìn)行處理,同未去噪的IMF分量進(jìn)行重構(gòu),完成信號(hào)去噪,并應(yīng)用于震源定位的預(yù)處理。
CEEMD算法是EMD的一種改進(jìn)算法,它既保留了EMD在處理非平穩(wěn)信號(hào)方面的優(yōu)點(diǎn),又能有效克服EMD的模態(tài)混疊問(wèn)題。CEEMD的主要步驟如下:
(1) 向原始信號(hào)中加入n組正負(fù)成對(duì)的白噪聲,生成2n組信號(hào)。
(2) 對(duì)2n組信號(hào)分別進(jìn)行EMD分解,得到2n組IMF分量cij。其中i表示添加輔助噪聲的第i個(gè)信號(hào),j表示對(duì)該信號(hào)進(jìn)行EMD分解得到的第j個(gè)IMF分量。
(3) 對(duì)2n組IMF分量對(duì)應(yīng)尺度取平均,得到CEEMD的對(duì)應(yīng)分解結(jié)果
(1)
多尺度排列熵(multiscale permutation entropy, MPE)是計(jì)算時(shí)間信號(hào)經(jīng)過(guò)不同尺度粗粒度運(yùn)算后的排列熵,可以分析信號(hào)的隨機(jī)性和復(fù)雜性。步驟如下:
(1) 對(duì)時(shí)間序列X={x(i),i=1,2,…,N}進(jìn)行粗?;幚淼玫酱至;蟮男蛄?/p>
(2)

(3)
式中:m為嵌入維數(shù);τ為延遲因子。
(4)
(5)
(4) 對(duì)尺度因子為s的序列的排列熵進(jìn)行計(jì)算并歸一化處理
(6)
通過(guò)MPE方法分析CEEMD分解得到的各階IMF的復(fù)雜性,判斷MPE值大于0.6的IMF為高頻分量,MPE值小于0.6的為低頻分量[18]。
小波閾值算法廣泛應(yīng)用于信號(hào)去噪和圖像去噪中。其原理是對(duì)信號(hào)進(jìn)行小波分解,得到一系列小波系數(shù),判斷大于閾值的小波系數(shù)為有用信號(hào),而確定小于閾值的小波系數(shù)為噪聲,可以去除[19],從而達(dá)到抑制噪聲的目的。小波閾值算法的具體步驟如下:
(1) 對(duì)信號(hào)進(jìn)行小波分解。在分解時(shí),需要選擇合適的小波基和分解層數(shù)。常用的小波基有dB小波和sym小波。經(jīng)過(guò)多次試驗(yàn)選擇了sym8小波。選擇3層作為分解層的數(shù)量。
(2) 閾值處理,選擇固定閾值。
(3) 閾值函數(shù)分為硬閾值函數(shù)和軟閾值函數(shù)。但軟硬閾值法本身也存在一些缺陷。硬閾值函數(shù)在閾值處不連續(xù),對(duì)于每一層細(xì)節(jié)分量,經(jīng)過(guò)閾值處理后得到的小波系數(shù)都會(huì)產(chǎn)生額外的振蕩。軟閾值函數(shù)是連續(xù)的,但經(jīng)過(guò)閾值處理后的小波系數(shù)與小波系數(shù)之間存在恒定的偏差。為了提高小波閾值去噪(wavelet threshold denoising, WTD)效果,采用軟硬閾值折衷函數(shù)作為閾值函數(shù)。
軟硬閾值折衷函數(shù)的表達(dá)式是
(7)
式中:ω為閾值處理前的小波系數(shù);ωthr為閾值處理后的小波系數(shù);thr為小波閾值;α為調(diào)整參數(shù),0≤α≤1。
(4) 信號(hào)重構(gòu)。重構(gòu)閾值處理后的近似分量和細(xì)節(jié)分量,得到小波閾值去噪后的信號(hào)。
改進(jìn)的小波閾值函數(shù),如圖1所示。

圖1 改進(jìn)的小波閾值函數(shù)
本文提出的去噪算法是將CEEMD、MPE和改進(jìn)的WTD算法相結(jié)合。算法流程圖如圖2所示。噪聲抑制方法的實(shí)現(xiàn)步驟如下:
(1) 對(duì)被測(cè)信號(hào)進(jìn)行CEEMD分解,得到各階IMF。
(2) 計(jì)算IMF的MPE值。當(dāng)IMF的MPE大于0.6時(shí),可以認(rèn)為是一個(gè)有噪聲的高頻分量,需要進(jìn)行去噪處理。當(dāng)IMF的MPE小于0.6時(shí),可以認(rèn)為是低頻分量。
(3) 采用改進(jìn)的WTD算法對(duì)高頻分量進(jìn)行去噪,并與低頻分量重構(gòu)。
由于由正弦和余弦信號(hào)組成的疊加信號(hào)與爆破振動(dòng)信號(hào)[20]具有相似的特性,我們利用疊加信號(hào)來(lái)驗(yàn)證所選方法的性能。仿真信號(hào)可以表示為

圖2 基于CEEMD、MPE和改進(jìn)WTD的噪聲抑制流程圖
(8)
式中,y4為高斯白噪聲。
分別使用EEMD-MPE方法、CEEMD-MPE方法、CEEMD-MPE-WTD方法和CEEMD-MPE-IMPROVED WTD對(duì)模擬信號(hào)y(t)進(jìn)行去噪處理,采用信噪比(signal-to-noise ratio, SNR)、均方根誤差(root mean square error, RMSE)和皮爾遜相關(guān)系數(shù)(Pearson correlation coefficient, PCC)作為抑制噪聲的指標(biāo),處理結(jié)果如表1所示。

表1 四種不同方法的噪聲抑制效果比較
由表1可知,CEEMD-MPE-IMPROVED WTD處理得到的信號(hào)信噪比最高,均方根誤差最小。通過(guò)計(jì)算四種去噪方法下純凈信號(hào)和去噪信號(hào)的PCC值,發(fā)現(xiàn)該方法去噪方法的PCC值最接近1,表明該方法中的去噪信號(hào)與純信號(hào)最相關(guān),能夠最有效地保留原始信號(hào)的有效信息。
結(jié)果表明,CEEMD-MPE-IMPROVED WTD算法效果最佳,為后續(xù)對(duì)實(shí)測(cè)隧道爆破振動(dòng)信號(hào)的降噪處理提供了理論依據(jù)。
CEEMD-MPE-IMPROVED WTD算法研究是基于江蘇省南京市的隧道爆炸試驗(yàn)。土地為相對(duì)平坦的黃土層,含一定程度的砂土,相對(duì)干燥堅(jiān)硬,表面有少量礫石。
在實(shí)際爆炸試驗(yàn)中,以隧道為爆炸主體,在隧道口放置炸藥進(jìn)行爆炸,在隧道外鋪設(shè)光纜,采用分布式光纖聲傳感系統(tǒng)(distributed acoustic sensing, DAS)結(jié)合光纜采集爆炸振動(dòng)信號(hào)。DAS是一種最先進(jìn)的光纖聲學(xué)檢測(cè)技術(shù),它使用光纖作為傳感器來(lái)檢測(cè)振動(dòng)信號(hào)。通過(guò)檢測(cè)光纖中瑞利后向散射的變化,幾乎可以在光纖上的任何點(diǎn)檢測(cè)振動(dòng)信息[21]?,F(xiàn)場(chǎng)平面測(cè)試圖如圖3所示?,F(xiàn)場(chǎng)試驗(yàn)的場(chǎng)地存在約為1 m的垂直高度差,測(cè)試場(chǎng)地長(zhǎng)度約為100 m,寬度約為50 m。本次爆炸試驗(yàn)中使用的DAS系統(tǒng)的空間分辨率為1 m,采樣頻率為20 000 Hz,最大測(cè)量長(zhǎng)度為20 km。

圖3 測(cè)試系統(tǒng)整體示意圖
DAS系統(tǒng)采集的震源時(shí)空數(shù)據(jù)可以組成矩陣。假設(shè)有L個(gè)傳感節(jié)點(diǎn),傳感節(jié)點(diǎn)個(gè)數(shù)由DAS系統(tǒng)的空間分辨率和光纖長(zhǎng)度決定,DAS系統(tǒng)采集時(shí)間序列長(zhǎng)度為N=fs×t;其中t為時(shí)間,t=1 s;fs為采樣率,fs=20 kHz。由此可得N行L列的DAS信號(hào)矩陣
(9)
光纜響應(yīng)圖即為原始數(shù)據(jù)X所繪制瀑布圖,其中橫軸代表光纜的位置信息,單位為m??v軸代表時(shí)間信息,單位為s。光纜響應(yīng)圖的顏色代表信號(hào)的強(qiáng)度,顏色越深,代表該該位置該時(shí)刻的光相位變化越明顯。其中某一位置的信號(hào)波形如圖4所示。光纜響應(yīng)圖如圖5所示。

圖4 爆破振動(dòng)信號(hào)

圖5 光纜響應(yīng)圖
在對(duì)采集到的爆炸信號(hào)進(jìn)行去噪時(shí),首先對(duì)一個(gè)信號(hào)進(jìn)行CEEMD分解,得到各階IMF,并繪制其頻譜圖,結(jié)果如圖6所示,我們可以看到通過(guò)測(cè)量信號(hào)的CEEMD分解獲得的各階IMF分量及其頻譜??梢钥闯?一階和二階IMF包含頻率約3 000 Hz的設(shè)備共振噪聲、環(huán)境噪聲和有用信息。如果直接丟棄一階和二階IMF分量,信號(hào)中的有用信息將丟失,信號(hào)將失真。因此,為了獲得真實(shí)信號(hào),有必要在盡可能抑制噪聲分量的同時(shí),將有用信息保留在高階IMF分量中。
在信號(hào)通過(guò)CEEMD方法分解得到各階IMF后,我們使用MPE方法來(lái)判斷各階IMF的類(lèi)型。IMF1~I(xiàn)MF2的MPE值大于0.6,這被視為噪聲分量,需要使用改進(jìn)的WTD進(jìn)行去噪;IMF3~I(xiàn)MF11的MPE值小于0.6,將其視為信號(hào)分量,并在去噪后用噪聲分量重構(gòu)。
為了說(shuō)明改進(jìn)WTD對(duì)高頻分量中噪聲分量的抑制效果,將CEEMD分解得到的IMF1分量與經(jīng)過(guò)改進(jìn)小波閾值處理后的IMF1分量進(jìn)行了比較,繪制了它們的波形和頻譜圖,如圖7所示。
如圖7所示,經(jīng)過(guò)改進(jìn)的小波閾值去噪后,時(shí)域波形的噪聲分量顯著減少,頻譜IMF1的噪聲分量的幅度從0.248降低到0.147,表明IMF1經(jīng)過(guò)此方法處理后確實(shí)可以抑制噪聲分量。然而,信號(hào)的有用成分也存在一定程度的降低。


圖6 CEEMD分解結(jié)果及其對(duì)應(yīng)的譜


(a) CEEMD分解得到IMF1


(b) 改進(jìn)小波閾值處理得到的IMF1
為了分析該去噪算法中有效成分的損失程度。我們將信噪比為25 dB的高斯白噪聲添加到采集的信號(hào)中,并分別通過(guò)EEMD-MPE、CEEMD-MPE、CEEMD-MPE-WTD和CEEMD-MPE-IMPROVED WTD算法進(jìn)行去噪。每個(gè)信號(hào)的波形圖和AOK時(shí)頻分析圖如圖8所示。X是主頻能量,Y是主頻。



(a) 原始信號(hào)



(b) EEMD-MPE 方法

(c) CEEMD-MPE 方法



(d) CEEMD-MPE-WTD 方法



(e) CEEMD-MPE-IMPROVED WTD方法
圖8顯示了原始信號(hào)和通過(guò)各種方法去噪的信號(hào)的波形圖和頻譜圖。信號(hào)的采樣時(shí)間為1 s。從圖8(a)中的頻譜圖可以看出,原始信號(hào)的主頻率能量為64.760 4,并且在3 053 Hz的頻率處存在峰值能量為3.84的噪聲分量。
從信號(hào)波形可以看出,CEEMD-MPE方法由于直接丟棄高頻IMF分量,3 053 Hz頻率處的噪聲分量消失,信號(hào)波形最平滑。似乎CEEMD-MPE方法后的信號(hào)去噪效果最好。但將原始信號(hào)的主頻能量與通過(guò)各種方法去噪的信號(hào)的主頻能量進(jìn)行比較時(shí),如表2所示。本文提出方法的主頻能量?jī)H減少了0.000 2,其他方法獲得的主頻能量損失ΔX=0.000 5~0.053 8,與本文提出方法相比能量損耗相差很多倍。且使用本文提出的去噪方法,3 053 Hz處信號(hào)噪聲分量的能量從3.84降低到1.53。由此可以判斷出,CEEMD-MPE-IMPROVED WTD算法相比于其他算法,不僅能很好地抑制信號(hào)中的噪聲,而且在很大程度上保留了信號(hào)的有效成分。

表2 各方法的信號(hào)主頻能量X和主頻Y的結(jié)果比較
為了進(jìn)一步證明CEEMD-MPE-IMPROVED WTD算法的優(yōu)越性,用上述各種方法對(duì)信號(hào)進(jìn)行去噪處理,并計(jì)算信號(hào)的SNR和RMSE。
表3為EEMD-MPE、CEEMD-MPE、CEEMD-MPE-WTD和CEEMD-MPE-IMPROVED WTD算法的噪聲抑制結(jié)果指標(biāo)。CEEMD-MPE-IMPROVED WTD算法相比于其他方法信噪比明顯提高,均方根誤差顯著降低。證明了噪聲抑制方法的合理性和有效性。

表3 四種方法的噪聲抑制結(jié)果比較
為了進(jìn)一步驗(yàn)證該方法的適用性,我們隨機(jī)選擇了采集的40個(gè)爆炸信號(hào),添加25 dB的高斯白噪聲,采用該方法應(yīng)用于這些數(shù)據(jù)的去噪,求取SNR和RMSE并作圖。
圖9顯示了采用四種不同方法對(duì)隨機(jī)選擇的40個(gè)振動(dòng)信號(hào)噪聲抑制后的SNR值。從圖9可以看出,CEEMD-MPE-IMPROVED WTD算法獲得的SNR值明顯高于其他算法。圖10顯示了采用四種不同方法對(duì)隨機(jī)選擇的40個(gè)振動(dòng)信號(hào)噪聲抑制后的RMSE值,CEEMD-MPE-IMPROVED WTD算法得到的RMSE值最低,證明了本文提出算法的優(yōu)勢(shì)。

圖9 四種方法的信噪比值的折線(xiàn)圖

圖10 四種方法的均方根誤差值的折線(xiàn)圖
為了驗(yàn)證本文提出的去噪方法在震源定位中的作用,分別利用去噪前后的數(shù)據(jù)進(jìn)行震源定位前的信號(hào)到時(shí)提取。已知DAS系統(tǒng)和光纜起始通道在(0,0)坐標(biāo)處,真實(shí)震源位置為(60,25)。
采用Geiger定位方法,首選確定一個(gè)目標(biāo)函數(shù),利用接收點(diǎn)到估計(jì)震源和到真實(shí)震源的走時(shí)差最小確定,并采用優(yōu)化算法使目標(biāo)函數(shù)值最小[22]。
目標(biāo)函數(shù)如下所示

(10)
式中:Ti為光纜第i個(gè)傳感點(diǎn)接收到振動(dòng)信號(hào)的時(shí)刻,即信號(hào)到達(dá)時(shí)間,利用赤池信息準(zhǔn)則(Akaike information criterion, AIC)算法求得;t0為炸藥爆炸時(shí)刻;v為振動(dòng)波在地下的傳播速度;(xi,yi,zi)為光纜第i個(gè)傳感點(diǎn)的坐標(biāo);(x,y,z)為假定震源位置;n為光纜上傳感點(diǎn)的個(gè)數(shù),且n=300。
由于現(xiàn)場(chǎng)試驗(yàn)的場(chǎng)地并非完全水平,存在約為1 m的垂直高度差。但是現(xiàn)場(chǎng)測(cè)試時(shí),測(cè)試場(chǎng)地長(zhǎng)度約為100 m,寬度約為50 m,由于高度引起的誤差可忽略。所以本稿件中目標(biāo)函數(shù)可簡(jiǎn)化為

(11)
利用遺傳算法優(yōu)化得到使誤差值ε最小的假定震源坐標(biāo)(x0,y0),認(rèn)為遺傳算法優(yōu)化得到的震源坐標(biāo)(x0,y0)即為真實(shí)震源位置。分別利用去噪前后信號(hào)進(jìn)行10次定位,得到光纜的實(shí)際響應(yīng)圖與定位結(jié)果圖。
如圖11所示為原始信號(hào)的到時(shí)提取結(jié)果以及震源定位結(jié)果圖。 圖11(a)中每個(gè)通道上的點(diǎn)代表算法捕捉的每個(gè)信號(hào)的起始時(shí)間。已知震源真實(shí)坐標(biāo)為(60,25),由于噪聲成分的干擾使得地震波初至?xí)r間判斷存在較大誤差,據(jù)此計(jì)算得到的震源位置與真實(shí)震源位置差距較大,定位結(jié)果如圖11(b)所示。

(a) 原始信號(hào)的光纜響應(yīng)圖

(b) 原始信號(hào)的定位結(jié)果圖
利用去噪后的信號(hào)重復(fù)定位過(guò)程,如圖12所示,到時(shí)提取的時(shí)間能較為準(zhǔn)確的反映地震波到達(dá)的起始時(shí)間,通過(guò)去噪后得到的到時(shí)計(jì)算的震源位置為(64,27),與震源真實(shí)位置(60,25)誤差為5 m左右。

(a) 去噪后信號(hào)的光纜響應(yīng)圖

(b) 去噪后信號(hào)的定位結(jié)果圖
為了進(jìn)一步確認(rèn)論文提出方法優(yōu)越性,分別采用原始數(shù)據(jù)、EEMD-MPE、CEEMD-MPE、CEEMD-MPE-WTD以及本文提出算法去噪后,利用相同的定位算法求出定位結(jié)果以及定位誤差,并繪制表格如表4所示。

表4 不同算法的定位結(jié)果對(duì)比
根據(jù)定位結(jié)果可以得出結(jié)論,該去噪算法對(duì)爆炸振動(dòng)信號(hào)有較好的去噪效果,可用于震源定位前對(duì)信號(hào)的預(yù)處理。
本文提出了一種基于互補(bǔ)集成經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)、多尺度置換熵(MPE)和改進(jìn)小波閾值去噪(WTD)的融合降噪方法對(duì)爆炸振動(dòng)信號(hào)進(jìn)行降噪處理。結(jié)論如下:
(1) CEEMD,MPE,以及改進(jìn)WTD的融合算法在處理隧道的爆破振動(dòng)信號(hào)時(shí)效果較好,可有效地去除隧道爆破振動(dòng)信號(hào)的高頻噪聲成分。
(2) 基于AOK時(shí)頻分析技術(shù),分析CEEMD-MPE-IMPROVED WTD算法得到的信號(hào)波形及主頻能量可知,該算法在有效去除噪聲成分的同時(shí),由于保留了高頻IMF分量,信號(hào)的主成分幾乎沒(méi)有損失。為后續(xù)分析爆源位置及炸藥威力奠定了基礎(chǔ)。
(3) 采用EEMD-MPE、CEEMD-MPE、CEEMD-MPE-WTD和CEEMD-MPE-IMPROVED WTD算法分別處理爆炸振動(dòng)信號(hào)時(shí),CEEMD-MPE-IMPROVED WTD算法處理后的信號(hào)的去噪效果最優(yōu)。對(duì)隨機(jī)40個(gè)信號(hào)進(jìn)行重復(fù)試驗(yàn),得出相同結(jié)論。驗(yàn)證了該算法的有效性和適用性,適合用于隧道爆炸振動(dòng)信號(hào)的去噪處理。
(4) 將CEEMD-MPE-IMPROVED WTD算法處理后的信號(hào)用于震源定位應(yīng)用中,定位誤差不超過(guò)5 m,遠(yuǎn)優(yōu)于原始信號(hào)定位結(jié)果,說(shuō)明該方法適合用于爆炸信號(hào)定位前的信號(hào)預(yù)處理過(guò)程。