CEEMDAN-小波閾值法在爆破振動(dòng)信號(hào)處理中的應(yīng)用

費(fèi)鴻祿，山 杰

(遼寧工程技術(shù)大學(xué) 爆破技術(shù)研究院，阜新 123000)

爆破振動(dòng)信號(hào)分析是礦業(yè)、巖土等工程界研究爆破振動(dòng)效應(yīng)的的主要方法，但在施工現(xiàn)場(chǎng)采集的的樣本信息中，受施工現(xiàn)場(chǎng)復(fù)雜環(huán)境的影響，采樣信息中含有大量高頻噪聲，同時(shí)還受到信號(hào)采樣儀器松動(dòng)或者溫度變化產(chǎn)生零點(diǎn)漂移的情況，采樣信息可能會(huì)出現(xiàn)形狀不規(guī)則和基線偏移情況，即采樣信息中產(chǎn)生了趨勢(shì)項(xiàng)[1]。兩者共同作用導(dǎo)致采樣信號(hào)失真，使爆破信號(hào)時(shí)頻能量的進(jìn)一步分析產(chǎn)生較大誤差。因此，采樣信號(hào)的去噪及趨勢(shì)項(xiàng)處理就成了不可或缺的前處理工作。

目前，國(guó)內(nèi)外常用的爆破振動(dòng)信號(hào)的去噪方法主要有小波閾值去噪[2]、小波包閾值去噪[3]、經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，EMD)方法[4]、集合經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition，EEMD)方法[5]、互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解(complementary ensemble empirical mode decomposition，CEEMD)方法[6]，以及兩者共同使用的EMD-小波閾值法[7]、EEMD-小波閾值法[8]、CEEMD-小波閾值法等[9]。小波類(lèi)方法對(duì)不同的頻率成分提供不同的分析分辨率，同時(shí)在時(shí)頻域的局部化性質(zhì)中有明顯優(yōu)勢(shì)，在非平穩(wěn)信號(hào)處理中備受青睞，但存在小波基難以選取，分解層數(shù)難以確定的缺點(diǎn)。目前對(duì)信號(hào)的趨勢(shì)項(xiàng)消除方法主要有EMD法[10]、EEMD法[11]、變分模態(tài)分解(variational mode decomposition，VMD)法等[12]，而同時(shí)對(duì)兩者進(jìn)行處理的方法中[13,14]，EMD方法主要問(wèn)題為端點(diǎn)效應(yīng)和模態(tài)混疊，EEMD方法雖然避免了EMD方法端點(diǎn)效應(yīng)和模態(tài)混疊現(xiàn)象，但同時(shí)引入了均勻分布的白噪聲。CEEMD方法通過(guò)加入一對(duì)互為相反數(shù)的輔助白噪聲解決EEMD方法存在的白噪聲問(wèn)題，但處理本質(zhì)沒(méi)變，分解結(jié)果難免受到殘余白噪聲影響。自適應(yīng)噪聲完備集合經(jīng)驗(yàn)?zāi)B(tài)分解(complete ensemble empirical mode decomposition with adaptive noise，CEEMDAN)方法是從EMD的基礎(chǔ)上加以改進(jìn)[15]，同時(shí)借用了EEMD方法中加入高斯白噪聲和通過(guò)多次疊加并集成平均以抵消噪聲的思想，能減輕端點(diǎn)效應(yīng)和模態(tài)混疊現(xiàn)象，避免殘余噪聲影響，但處理效果還是不夠理想，容易丟失部分真實(shí)信息。

為解決CEEMDAN分解丟失信息以及小波基難以選取、分解層數(shù)難以確定的問(wèn)題，將自適應(yīng)噪聲完備集合經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMDAN)和小波閾值法的優(yōu)勢(shì)結(jié)合起來(lái)，輔以頻譜篩選并去除趨勢(shì)項(xiàng)、互相關(guān)與自相關(guān)函數(shù)分析與校核噪聲分量，提出了一種處理爆破振動(dòng)信號(hào)趨勢(shì)項(xiàng)與高頻噪聲的方法，對(duì)爆破振動(dòng)響應(yīng)研究具有重要意義。

1 爆破振動(dòng)信號(hào)處理理論

1.1 CEEMDAN分解原理

EEMD、CEEMD分解方法是將添加白噪聲后的信號(hào)直接做EMD分解，然后相對(duì)應(yīng)的IMF分量直接求均值。CEEMDAN算法是每求完一階IMF分量，又重新給殘余分量加入正態(tài)分布的高斯白噪聲并求此時(shí)的IMF分量均值并逐次迭代。其詳細(xì)實(shí)現(xiàn)過(guò)程如下[16]：

(1)將標(biāo)準(zhǔn)正態(tài)分布的高斯白噪聲wi(t)加入到爆破振動(dòng)信號(hào)x(t)中，則待處理信號(hào)變?yōu)閤(t)+ζ0wi(t)。使用EMD方法對(duì)待處理信號(hào)進(jìn)行I次重復(fù)分解后對(duì)其進(jìn)行算數(shù)平均可得到IMF1分量為

(1)

去除第一階分量IMF1后的殘余分量

r1(t)=x(t)-IMF1(t)

(2)

(2)定義從待處理信號(hào)的EMD分解中獲得的第j個(gè)模態(tài)函數(shù)分量的算子為Ej(·)。再次將標(biāo)準(zhǔn)正態(tài)分布的高斯白噪聲wi(t)加入到分量r1(t)中，則對(duì)信號(hào)r1(t)+ζ1E1[wi(t)]再次通過(guò)EMD方法進(jìn)行I次重復(fù)分解，可得到IMF2分量為

(3)

去除第二階分量IMF2后的殘余分量

r2(t)=r1(t)-IMF2(t)

(4)

(3)對(duì)于k=2,3,…,K，去除IMFk分量后的殘余分量為

rk(t)=rk-1(t)-IMFk(t)

(5)

(4)提取信號(hào)rk(t)+ζkEk[wi(t)]的第一個(gè)IMF分量，定義k+1個(gè)模態(tài)函數(shù)分量為

(6)

(5)重復(fù)執(zhí)行步驟(3)和(4)，直到殘余分量信號(hào)不可再繼續(xù)被分解為止，最終可以得到K個(gè)模態(tài)函數(shù)分量。

原始爆破振動(dòng)信號(hào)可以表示為

(7)

式中，R(t)為最終殘余分量。

從上述CEEMDAN的基本理論可以看出，此方法在較小的平均次數(shù)下就可以有很好的完備性與模態(tài)分解結(jié)果，計(jì)算速度上也有明顯優(yōu)勢(shì)，且可以改變每一階段的噪聲系數(shù)ζi，自適應(yīng)的選擇不同信噪比的高斯白噪聲。

1.2 小波閾值方法

小波閾值法去噪的本質(zhì)是對(duì)信號(hào)的濾波。小波閾值法是先通過(guò)小波變換把信號(hào)轉(zhuǎn)換為不同尺度的尺度空間，然后通過(guò)閾值處理篩選出噪聲并過(guò)濾掉噪聲，最后再重構(gòu)信號(hào)。小波閾值法在爆破振動(dòng)信號(hào)濾波去噪中，閾值處理至關(guān)重要。閾值處理包括閾值選擇與閾值函數(shù)選擇。目前閾值選擇方法有4種，RigSure閾值、Heursure閾值、Sqtwolog閾值、Minmax閾值。常用的閾值量化規(guī)則一般分為硬閾值與軟閾值，硬閾值函數(shù)處理信號(hào)后會(huì)在閾值λ處產(chǎn)生局部突變。軟閾值函數(shù)處理得到的信號(hào)雖然丟失了一小部分高頻信號(hào)，但信號(hào)更平滑，消除了硬閾值函數(shù)處理引起的局部突變。

小波去噪硬閾值處理函數(shù)

(8)

小波去噪軟閾值處理函數(shù)

(9)

式中：w為小波系數(shù)；λ為閾值。

1.3 CEEMDAN-小波閾值處理方法

CEEMDAN-小波閾值方法先將爆破振動(dòng)信號(hào)進(jìn)行CEEMDAN分解處理，在判斷出趨勢(shì)項(xiàng)分量后，對(duì)剩余IMF分量進(jìn)行互相關(guān)系數(shù)和自相關(guān)函數(shù)分析確定噪聲IMF分量，對(duì)其進(jìn)行小波閾值處理過(guò)濾掉噪聲信息，提取CEEMDAN分解丟失的部分有效信息。最后，將小波閾值濾波處理后的IMF分量與純凈IMF分量重構(gòu)得到處理后信號(hào)。完整流程如圖1所示。

圖 1 信號(hào)處理流程圖Fig. 1 Signal processing flowchart

2 工程實(shí)例分析及效果評(píng)價(jià)

選取武家塔露天礦實(shí)測(cè)的某次爆破振動(dòng)徑向信號(hào)，如圖2所示?？梢钥闯鲈撔盘?hào)由于現(xiàn)場(chǎng)環(huán)境復(fù)雜受到較多噪聲污染且有很明顯的趨勢(shì)項(xiàng)情況，信號(hào)波形偏離基線中央，出現(xiàn)失真情況。為了不影響信號(hào)時(shí)頻能量的進(jìn)一步分析，需要對(duì)其進(jìn)行處理。

2.1 CEEMDAN分解

通過(guò)反復(fù)實(shí)驗(yàn)，加入白噪聲的噪聲標(biāo)準(zhǔn)差與振動(dòng)徑向信號(hào)標(biāo)準(zhǔn)差的比值(幅值)為0.2，信號(hào)的平均次數(shù)為100，最大迭代次數(shù)為2000時(shí)，CEEMDAN分解處理效果最好。11個(gè)IMF分量(IMF1～I(xiàn)MF11)和一個(gè)殘余分量Res由原始爆破振動(dòng)徑向信號(hào)通過(guò)CEEMDAN分解得到，結(jié)果如圖3所示。

圖 2 振動(dòng)徑向信號(hào)Fig. 2 Vibration radial signal

2.2 振動(dòng)信號(hào)處理

2.2.1 趨勢(shì)項(xiàng)處理

對(duì)各個(gè)IMF分量進(jìn)行FFT(Fast Fourier Transfer)變換，得到各個(gè)IMF分量的主頻，結(jié)果如表1所示。由于趨勢(shì)項(xiàng)為低頻干擾，對(duì)疑似趨勢(shì)項(xiàng)IMF9、IMF10、IMF11和Res分量進(jìn)行頻譜主頻帶分析，結(jié)果如圖4所示。從圖4可以看出，IMF11、Res分量的頻帶主要集中分布在0～5 Hz的低頻范圍內(nèi)，且主頻低于爆破測(cè)振儀監(jiān)測(cè)的有效范圍(2～200 Hz)，可能為零漂或信號(hào)變化趨勢(shì)[17]。相比之下，IMF9、IMF10頻帶更寬，主頻在有效監(jiān)測(cè)區(qū)間之內(nèi)，因此去掉IMF11、Res趨勢(shì)項(xiàng)分量。

圖 3 CEEMDAN分解結(jié)果Fig. 3 CEEMDAN decomposition results

表 1 各主頻對(duì)應(yīng)表

圖 4 疑似趨勢(shì)項(xiàng)頻譜Fig. 4 Spectrum of suspected trend items

2.2.2 噪聲處理

將剩余IMF分量與振動(dòng)徑向信號(hào)通過(guò)MATLAB中的corrcoef函數(shù)求出互相關(guān)系數(shù)R，結(jié)果如表2所示。從表2可以看出IMF1、IMF2、IMF3、IMF4和IMF5分量的互相關(guān)系數(shù)均小于0.1，初步認(rèn)定為含噪IMF分量[18]。針對(duì)IMF1、IMF2、IMF3、IMF4和IMF5疑似高頻噪聲分量，通過(guò)MATLAB中的xcorr函數(shù)對(duì)其進(jìn)行自相關(guān)函數(shù)分析，結(jié)果如圖5。

表 2 剩余分量互相關(guān)系數(shù)對(duì)應(yīng)表

圖 5 疑似噪聲分量自相關(guān)函數(shù)Fig. 5 Autocorrelation function of suspected noise component

從圖5(a)～(d)分析可知，分量IMF1～I(xiàn)MF4在零點(diǎn)處的自相關(guān)函數(shù)達(dá)到最大值，其他時(shí)間的自相關(guān)函數(shù)趨近于0，因此定義分量IMF1～I(xiàn)MF4為高頻噪聲分量[19]。由圖5(e)可知，IMF5分量雖然在零點(diǎn)處自相關(guān)函數(shù)為最大值，但其余時(shí)間處并不完全符合噪聲特性，由表1也可知IMF5分量主頻為180.7 Hz，主頻在有效監(jiān)測(cè)區(qū)內(nèi)，因而IMF5分量同時(shí)也包含了一些信號(hào)的有效信息，需要對(duì)其進(jìn)行小波閾值濾波去噪處理。由于本文中爆破測(cè)振儀器設(shè)置的最小采樣頻率為2 Hz，根據(jù)采樣定理[20]，設(shè)置的采樣頻率為4000 Hz，奈奎斯特頻率為2000 Hz，選取小波基函數(shù)“db8”對(duì)IMF5分量進(jìn)行8層分解，其最低頻帶為0～7.8125 Hz。選取基于Stein的無(wú)偏移風(fēng)險(xiǎn)估計(jì)原理獲得小波閾值和具有光滑性的軟閾值函數(shù)進(jìn)行濾波去噪處理。

2.3 波形重構(gòu)

將純凈IMF分量與濾波去噪處理后分量疊加得到處理后信號(hào)。處理(去除趨勢(shì)項(xiàng)及噪聲)前后信號(hào)波形及頻譜如圖6和7所示。分析圖6可知：處理后的信號(hào)波形回到基線中央，曲線更加光滑，消除了噪聲污染影響，保留了信號(hào)真實(shí)信息。此外，由圖7(a)可知，處理前的信號(hào)頻譜低頻段存在突高的相對(duì)幅值，高頻段存在明顯的噪聲相對(duì)幅值，嚴(yán)重影響了信號(hào)頻譜分析的準(zhǔn)確性。由圖7(b)可知，處理后的信號(hào)在0～2 Hz低頻段較大相對(duì)幅值已經(jīng)消除，表明去趨勢(shì)項(xiàng)效果顯著。2～200 Hz頻段振動(dòng)信號(hào)無(wú)明顯變化，表明處理方法消除噪聲干擾外，并不會(huì)影響信號(hào)的真實(shí)信息，而200 Hz以上頻段相對(duì)幅值基本為零，說(shuō)明處理方法去除了高頻噪聲污染。

圖 6 處理前后信號(hào)波形Fig. 6 Signal waveform before and after processing

圖 7 處理前后信號(hào)頻譜Fig. 7 Signal spectrum before and after processing

2.4 對(duì)比分析

為了反應(yīng)處理方法的效果，將信噪比(SNR)、均方根誤差(RMSE)作為處理效果的衡量指標(biāo)，一般認(rèn)為SNR越大，RMSE越小，其處理效果越好。其計(jì)算公式如下

(10)

(11)

式中：Xi為處理前信號(hào)；xi為處理后信號(hào)；N為信號(hào)長(zhǎng)度。

將本文處理方法與單純的CEEMDAN方法、小波閾值方法以及EMD-小波閾值法、EEMD-小波閾值法和CEEMD-小波閾值法在衡量指標(biāo)上進(jìn)行對(duì)比，對(duì)比結(jié)果如表3。

從表3可以看出，在信噪比方面，CEEMDAN方法、小波閾值法與本文方法相比，本文處理方法最高，CEEMDAN方法次之，小波閾值方法最小，說(shuō)明CEEMDAN方法處理不夠理想，容易丟失部分真實(shí)信息，而小波閾值處理效果不如其他兩種方法。本文較EMD-小波閾值法、EEMD-小波閾值法和CEEMD-小波閾值法相比信噪比也有所提高，說(shuō)明本文處理效果更好；在均方根誤差方面，CEEMDAN-小波閾值法較其他幾種方法相比也顯示出更好的可靠性與穩(wěn)定性。綜上所述，經(jīng)過(guò)與目前廣泛采用的多種信號(hào)處理方法比較，本文應(yīng)用的CEEMDAN-小波閾值處理方法最優(yōu)。

表 3 衡量指標(biāo)對(duì)比

為了進(jìn)一步反映CEEMDAN-小波閾值方法的處理效果，使用Hilbert變換獲得時(shí)間、頻率和能量的三維圖來(lái)反應(yīng)處理前后的變化情況。處理前后的三維圖如圖8所示。

圖 8 處理前后三維圖Fig. 8 Three-dimensional image before and after processing

從圖8(a)分析可知，由于趨勢(shì)項(xiàng)的干擾，源信號(hào)在0～0.6 s的高頻段(2000Hz左右)以及0.6～1 s的低頻段(0～2 Hz)出現(xiàn)較大能量的分量，而由于噪聲干擾，在200～2000 Hz頻段上也出現(xiàn)了較小的噪聲能量，這就導(dǎo)致了信號(hào)的失真，從而增大爆破振動(dòng)信號(hào)分析的誤差甚至致使其結(jié)果嚴(yán)重不符合實(shí)際。從圖8(b)可以看出，上述由于趨勢(shì)項(xiàng)和噪聲引起的能量分量通過(guò)CEEMDAN-小波閾值方法處理后被成功剔除，而0～200 Hz頻段能量沒(méi)有明顯變化，說(shuō)明CEEMDAN-小波閾值方法不僅能消除趨勢(shì)項(xiàng)和噪聲干擾，而且對(duì)2～200 Hz頻段能量信息保留完整。

3 結(jié)論

針對(duì)爆破施工中采集的振動(dòng)信號(hào)的高頻噪聲和趨勢(shì)項(xiàng)問(wèn)題，使用CEEMDAN分解和小波閾值法共同處理信號(hào)。通過(guò)實(shí)例分析，得出以下結(jié)論：

(1)結(jié)合CEEMDAN分解與小波閾值法的優(yōu)勢(shì)處理爆破振動(dòng)信號(hào)的高頻噪聲與趨勢(shì)項(xiàng)問(wèn)題，為趨勢(shì)項(xiàng)的準(zhǔn)確篩選與去除做了新的研討，相比于目前廣泛采用的方法信噪比更高，均方根誤差更小，提高了爆破振動(dòng)信號(hào)處理的精度。

(2)對(duì)比CEEMDAN-小波閾值方法處理前后的信號(hào)波形、頻譜及三維時(shí)頻能量分析，證明本文方法能有效處理高頻噪聲及趨勢(shì)項(xiàng)的干擾，而且不會(huì)影響信號(hào)的真實(shí)頻率、能量信息，更利于爆破振動(dòng)響應(yīng)研究中的應(yīng)用。