費鴻祿,山 杰
(遼寧工程技術(shù)大學 爆破技術(shù)研究院,阜新 123000)
爆破振動信號分析是礦業(yè)、巖土等工程界研究爆破振動效應(yīng)的的主要方法,但在施工現(xiàn)場采集的的樣本信息中,受施工現(xiàn)場復雜環(huán)境的影響,采樣信息中含有大量高頻噪聲,同時還受到信號采樣儀器松動或者溫度變化產(chǎn)生零點漂移的情況,采樣信息可能會出現(xiàn)形狀不規(guī)則和基線偏移情況,即采樣信息中產(chǎn)生了趨勢項[1]。兩者共同作用導致采樣信號失真,使爆破信號時頻能量的進一步分析產(chǎn)生較大誤差。因此,采樣信號的去噪及趨勢項處理就成了不可或缺的前處理工作。
目前,國內(nèi)外常用的爆破振動信號的去噪方法主要有小波閾值去噪[2]、小波包閾值去噪[3]、經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)方法[4]、集合經(jīng)驗模態(tài)分解(ensemble empirical mode decomposition,EEMD)方法[5]、互補集合經(jīng)驗模態(tài)分解(complementary ensemble empirical mode decomposition,CEEMD)方法[6],以及兩者共同使用的EMD-小波閾值法[7]、EEMD-小波閾值法[8]、CEEMD-小波閾值法等[9]。小波類方法對不同的頻率成分提供不同的分析分辨率,同時在時頻域的局部化性質(zhì)中有明顯優(yōu)勢,在非平穩(wěn)信號處理中備受青睞,但存在小波基難以選取,分解層數(shù)難以確定的缺點。目前對信號的趨勢項消除方法主要有EMD法[10]、EEMD法[11]、變分模態(tài)分解(variational mode decomposition,VMD)法等[12],而同時對兩者進行處理的方法中[13,14],EMD方法主要問題為端點效應(yīng)和模態(tài)混疊,EEMD方法雖然避免了EMD方法端點效應(yīng)和模態(tài)混疊現(xiàn)象,但同時引入了均勻分布的白噪聲。CEEMD方法通過加入一對互為相反數(shù)的輔助白噪聲解決EEMD方法存在的白噪聲問題,但處理本質(zhì)沒變,分解結(jié)果難免受到殘余白噪聲影響。自適應(yīng)噪聲完備集合經(jīng)驗模態(tài)分解(complete ensemble empirical mode decomposition with adaptive noise,CEEMDAN)方法是從EMD的基礎(chǔ)上加以改進[15],同時借用了EEMD方法中加入高斯白噪聲和通過多次疊加并集成平均以抵消噪聲的思想,能減輕端點效應(yīng)和模態(tài)混疊現(xiàn)象,避免殘余噪聲影響,但處理效果還是不夠理想,容易丟失部分真實信息。
為解決CEEMDAN分解丟失信息以及小波基難以選取、分解層數(shù)難以確定的問題,將自適應(yīng)噪聲完備集合經(jīng)驗模態(tài)分解(CEEMDAN)和小波閾值法的優(yōu)勢結(jié)合起來,輔以頻譜篩選并去除趨勢項、互相關(guān)與自相關(guān)函數(shù)分析與校核噪聲分量,提出了一種處理爆破振動信號趨勢項與高頻噪聲的方法,對爆破振動響應(yīng)研究具有重要意義。
EEMD、CEEMD分解方法是將添加白噪聲后的信號直接做EMD分解,然后相對應(yīng)的IMF分量直接求均值。CEEMDAN算法是每求完一階IMF分量,又重新給殘余分量加入正態(tài)分布的高斯白噪聲并求此時的IMF分量均值并逐次迭代。其詳細實現(xiàn)過程如下[16]:
(1)將標準正態(tài)分布的高斯白噪聲wi(t)加入到爆破振動信號x(t)中,則待處理信號變?yōu)閤(t)+ζ0wi(t)。使用EMD方法對待處理信號進行I次重復分解后對其進行算數(shù)平均可得到IMF1分量為
(1)
去除第一階分量IMF1后的殘余分量
r1(t)=x(t)-IMF1(t)
(2)
(2)定義從待處理信號的EMD分解中獲得的第j個模態(tài)函數(shù)分量的算子為Ej(·)。再次將標準正態(tài)分布的高斯白噪聲wi(t)加入到分量r1(t)中,則對信號r1(t)+ζ1E1[wi(t)]再次通過EMD方法進行I次重復分解,可得到IMF2分量為
(3)
去除第二階分量IMF2后的殘余分量
r2(t)=r1(t)-IMF2(t)
(4)
(3)對于k=2,3,…,K,去除IMFk分量后的殘余分量為
rk(t)=rk-1(t)-IMFk(t)
(5)
(4)提取信號rk(t)+ζkEk[wi(t)]的第一個IMF分量,定義k+1個模態(tài)函數(shù)分量為
(6)
(5)重復執(zhí)行步驟(3)和(4),直到殘余分量信號不可再繼續(xù)被分解為止,最終可以得到K個模態(tài)函數(shù)分量。
原始爆破振動信號可以表示為
(7)
式中,R(t)為最終殘余分量。
從上述CEEMDAN的基本理論可以看出,此方法在較小的平均次數(shù)下就可以有很好的完備性與模態(tài)分解結(jié)果,計算速度上也有明顯優(yōu)勢,且可以改變每一階段的噪聲系數(shù)ζi,自適應(yīng)的選擇不同信噪比的高斯白噪聲。
小波閾值法去噪的本質(zhì)是對信號的濾波。小波閾值法是先通過小波變換把信號轉(zhuǎn)換為不同尺度的尺度空間,然后通過閾值處理篩選出噪聲并過濾掉噪聲,最后再重構(gòu)信號。小波閾值法在爆破振動信號濾波去噪中,閾值處理至關(guān)重要。閾值處理包括閾值選擇與閾值函數(shù)選擇。目前閾值選擇方法有4種,RigSure閾值、Heursure閾值、Sqtwolog閾值、Minmax閾值。常用的閾值量化規(guī)則一般分為硬閾值與軟閾值,硬閾值函數(shù)處理信號后會在閾值λ處產(chǎn)生局部突變。軟閾值函數(shù)處理得到的信號雖然丟失了一小部分高頻信號,但信號更平滑,消除了硬閾值函數(shù)處理引起的局部突變。
小波去噪硬閾值處理函數(shù)
(8)
小波去噪軟閾值處理函數(shù)
(9)
式中:w為小波系數(shù);λ為閾值。
CEEMDAN-小波閾值方法先將爆破振動信號進行CEEMDAN分解處理,在判斷出趨勢項分量后,對剩余IMF分量進行互相關(guān)系數(shù)和自相關(guān)函數(shù)分析確定噪聲IMF分量,對其進行小波閾值處理過濾掉噪聲信息,提取CEEMDAN分解丟失的部分有效信息。最后,將小波閾值濾波處理后的IMF分量與純凈IMF分量重構(gòu)得到處理后信號。完整流程如圖1所示。
圖 1 信號處理流程圖Fig. 1 Signal processing flowchart
選取武家塔露天礦實測的某次爆破振動徑向信號,如圖2所示??梢钥闯鲈撔盘栍捎诂F(xiàn)場環(huán)境復雜受到較多噪聲污染且有很明顯的趨勢項情況,信號波形偏離基線中央,出現(xiàn)失真情況。為了不影響信號時頻能量的進一步分析,需要對其進行處理。
通過反復實驗,加入白噪聲的噪聲標準差與振動徑向信號標準差的比值(幅值)為0.2,信號的平均次數(shù)為100,最大迭代次數(shù)為2000時,CEEMDAN分解處理效果最好。11個IMF分量(IMF1~IMF11)和一個殘余分量Res由原始爆破振動徑向信號通過CEEMDAN分解得到,結(jié)果如圖3所示。
圖 2 振動徑向信號Fig. 2 Vibration radial signal
2.2.1 趨勢項處理
對各個IMF分量進行FFT(Fast Fourier Transfer)變換,得到各個IMF分量的主頻,結(jié)果如表1所示。由于趨勢項為低頻干擾,對疑似趨勢項IMF9、IMF10、IMF11和Res分量進行頻譜主頻帶分析,結(jié)果如圖4所示。從圖4可以看出,IMF11、Res分量的頻帶主要集中分布在0~5 Hz的低頻范圍內(nèi),且主頻低于爆破測振儀監(jiān)測的有效范圍(2~200 Hz),可能為零漂或信號變化趨勢[17]。相比之下,IMF9、IMF10頻帶更寬,主頻在有效監(jiān)測區(qū)間之內(nèi),因此去掉IMF11、Res趨勢項分量。
圖 3 CEEMDAN分解結(jié)果Fig. 3 CEEMDAN decomposition results
表 1 各主頻對應(yīng)表
圖 4 疑似趨勢項頻譜Fig. 4 Spectrum of suspected trend items
2.2.2 噪聲處理
將剩余IMF分量與振動徑向信號通過MATLAB中的corrcoef函數(shù)求出互相關(guān)系數(shù)R,結(jié)果如表2所示。從表2可以看出IMF1、IMF2、IMF3、IMF4和IMF5分量的互相關(guān)系數(shù)均小于0.1,初步認定為含噪IMF分量[18]。針對IMF1、IMF2、IMF3、IMF4和IMF5疑似高頻噪聲分量,通過MATLAB中的xcorr函數(shù)對其進行自相關(guān)函數(shù)分析,結(jié)果如圖5。
表 2 剩余分量互相關(guān)系數(shù)對應(yīng)表
圖 5 疑似噪聲分量自相關(guān)函數(shù)Fig. 5 Autocorrelation function of suspected noise component
從圖5(a)~(d)分析可知,分量IMF1~IMF4在零點處的自相關(guān)函數(shù)達到最大值,其他時間的自相關(guān)函數(shù)趨近于0,因此定義分量IMF1~IMF4為高頻噪聲分量[19]。由圖5(e)可知,IMF5分量雖然在零點處自相關(guān)函數(shù)為最大值,但其余時間處并不完全符合噪聲特性,由表1也可知IMF5分量主頻為180.7 Hz,主頻在有效監(jiān)測區(qū)內(nèi),因而IMF5分量同時也包含了一些信號的有效信息,需要對其進行小波閾值濾波去噪處理。由于本文中爆破測振儀器設(shè)置的最小采樣頻率為2 Hz,根據(jù)采樣定理[20],設(shè)置的采樣頻率為4000 Hz,奈奎斯特頻率為2000 Hz,選取小波基函數(shù)“db8”對IMF5分量進行8層分解,其最低頻帶為0~7.8125 Hz。選取基于Stein的無偏移風險估計原理獲得小波閾值和具有光滑性的軟閾值函數(shù)進行濾波去噪處理。
將純凈IMF分量與濾波去噪處理后分量疊加得到處理后信號。處理(去除趨勢項及噪聲)前后信號波形及頻譜如圖6和7所示。分析圖6可知:處理后的信號波形回到基線中央,曲線更加光滑,消除了噪聲污染影響,保留了信號真實信息。此外,由圖7(a)可知,處理前的信號頻譜低頻段存在突高的相對幅值,高頻段存在明顯的噪聲相對幅值,嚴重影響了信號頻譜分析的準確性。由圖7(b)可知,處理后的信號在0~2 Hz低頻段較大相對幅值已經(jīng)消除,表明去趨勢項效果顯著。2~200 Hz頻段振動信號無明顯變化,表明處理方法消除噪聲干擾外,并不會影響信號的真實信息,而200 Hz以上頻段相對幅值基本為零,說明處理方法去除了高頻噪聲污染。
圖 6 處理前后信號波形Fig. 6 Signal waveform before and after processing
圖 7 處理前后信號頻譜Fig. 7 Signal spectrum before and after processing
為了反應(yīng)處理方法的效果,將信噪比(SNR)、均方根誤差(RMSE)作為處理效果的衡量指標,一般認為SNR越大,RMSE越小,其處理效果越好。其計算公式如下
(10)
(11)
式中:Xi為處理前信號;xi為處理后信號;N為信號長度。
將本文處理方法與單純的CEEMDAN方法、小波閾值方法以及EMD-小波閾值法、EEMD-小波閾值法和CEEMD-小波閾值法在衡量指標上進行對比,對比結(jié)果如表3。
從表3可以看出,在信噪比方面,CEEMDAN方法、小波閾值法與本文方法相比,本文處理方法最高,CEEMDAN方法次之,小波閾值方法最小,說明CEEMDAN方法處理不夠理想,容易丟失部分真實信息,而小波閾值處理效果不如其他兩種方法。本文較EMD-小波閾值法、EEMD-小波閾值法和CEEMD-小波閾值法相比信噪比也有所提高,說明本文處理效果更好;在均方根誤差方面,CEEMDAN-小波閾值法較其他幾種方法相比也顯示出更好的可靠性與穩(wěn)定性。綜上所述,經(jīng)過與目前廣泛采用的多種信號處理方法比較,本文應(yīng)用的CEEMDAN-小波閾值處理方法最優(yōu)。
表 3 衡量指標對比
為了進一步反映CEEMDAN-小波閾值方法的處理效果,使用Hilbert變換獲得時間、頻率和能量的三維圖來反應(yīng)處理前后的變化情況。處理前后的三維圖如圖8所示。
圖 8 處理前后三維圖Fig. 8 Three-dimensional image before and after processing
從圖8(a)分析可知,由于趨勢項的干擾,源信號在0~0.6 s的高頻段(2000Hz左右)以及0.6~1 s的低頻段(0~2 Hz)出現(xiàn)較大能量的分量,而由于噪聲干擾,在200~2000 Hz頻段上也出現(xiàn)了較小的噪聲能量,這就導致了信號的失真,從而增大爆破振動信號分析的誤差甚至致使其結(jié)果嚴重不符合實際。從圖8(b)可以看出,上述由于趨勢項和噪聲引起的能量分量通過CEEMDAN-小波閾值方法處理后被成功剔除,而0~200 Hz頻段能量沒有明顯變化,說明CEEMDAN-小波閾值方法不僅能消除趨勢項和噪聲干擾,而且對2~200 Hz頻段能量信息保留完整。
針對爆破施工中采集的振動信號的高頻噪聲和趨勢項問題,使用CEEMDAN分解和小波閾值法共同處理信號。通過實例分析,得出以下結(jié)論:
(1)結(jié)合CEEMDAN分解與小波閾值法的優(yōu)勢處理爆破振動信號的高頻噪聲與趨勢項問題,為趨勢項的準確篩選與去除做了新的研討,相比于目前廣泛采用的方法信噪比更高,均方根誤差更小,提高了爆破振動信號處理的精度。
(2)對比CEEMDAN-小波閾值方法處理前后的信號波形、頻譜及三維時頻能量分析,證明本文方法能有效處理高頻噪聲及趨勢項的干擾,而且不會影響信號的真實頻率、能量信息,更利于爆破振動響應(yīng)研究中的應(yīng)用。