爆破振動信號分析中模態(tài)混疊和虛假分量消除的改進(jìn)方法

李 清， 徐文龍， 張 迪， 李 娜， 馮丹丹

(1.中國礦業(yè)大學(xué)(北京) 力學(xué)與建筑工程學(xué)院，北京 100083; 2.廣州大學(xué) 土木工程學(xué)院，廣州 510006)

爆破振動信號是爆破設(shè)計參數(shù)的外在表現(xiàn)。通過分析信號的峰值振速、主頻、能量在不同時間和不同頻帶下的分布特征能有效的控制爆破振動強度，減小爆破危害。以傅氏變換為基礎(chǔ)的信號變換方法主要針對平穩(wěn)信號分析，顯然對于爆破振動信號分析存在局限性。Huang等[1]提出了HHT變換，并廣泛應(yīng)用于爆破振動等非平穩(wěn)信號的分析。HHT變換分成兩步，首先對信號進(jìn)行EMD分解，得到一組頻率由高到低的模態(tài)函數(shù)IMF，然后對IMF進(jìn)行Hilbert變換，得到信號的時頻圖譜、能量譜等。但是EMD中存在一些需要改進(jìn)的缺陷，如模態(tài)混疊和虛假分量。即任一IMF里含有明顯不相同的時間尺度或同一范圍的時間尺度出現(xiàn)在不同的IMF中。

針對模態(tài)混疊問題，Yang等[2]分析得到，在間斷信號處包絡(luò)線產(chǎn)生的過沖是產(chǎn)生模態(tài)混疊的重要原因。胡愛軍等[3]采用非線性形態(tài)濾波器，通過布爾運算和相應(yīng)的加減計算來分離信號內(nèi)的脈沖分量，保留原始信號的局部特性，減少信號中的脈沖噪聲來抑制模態(tài)混疊現(xiàn)象。高云超等[4]采用對混合分量一階差分的方法對分量進(jìn)行再分解，對分解得到的新的待驗證IMF分量進(jìn)行累計求和，消除模態(tài)混疊現(xiàn)象。周穎濤等[5]通過改進(jìn)EMD篩分停止閥值，增加單次篩分，使包絡(luò)平均值小于原始幅值抑制信號分解的的模態(tài)混疊現(xiàn)象。鄭近德等[6]檢測各IMF分量排列熵，剔除前幾階異常分量，將剩余分量再次進(jìn)行EMD分解抑制模態(tài)混疊現(xiàn)象。Ryan等[7]提出了對原始信號添加合適的掩膜信號來抑制模態(tài)混疊。趙玲等[8]、凌立鵬等[9]分別對Ryan的方法進(jìn)行了改進(jìn)，認(rèn)為由于EMD能量泄露的影響,僅通過第一個IMF分量的包絡(luò)幅值和瞬時頻率來定義s(t)的幅值和頻率對模態(tài)混疊的抑制效果較差，改進(jìn)的方法中將表征能量泄露分量的IMF的幅值與頻率和IMF1相加獲得新的s(t)。當(dāng)信號內(nèi)噪聲成分較多時，掩膜信號的確定十分困難，不便于操作。秦品樂等[10]對信號進(jìn)行先期小波閥值消噪，消除信號內(nèi)的突變成分，降低包絡(luò)線的波動。

關(guān)于虛假分量的判斷，黃迪山[11]根據(jù)EMD分解的完備性，利用能量守恒準(zhǔn)則，認(rèn)為當(dāng)EMD分解產(chǎn)生虛假分量時，各分量的均方值和要大于原始信號的均方值。將第一階本征模態(tài)分量與虛假分量相加，消除虛假分量的同時更新第一階模態(tài)分量，得到最終的各階分量。胡愛軍[12]基于互信息理論，將各IMF分量與原始信號進(jìn)行互信息計算判斷虛假分量。李紀(jì)永等[13]通過改進(jìn)的SVD對原始信號降噪，對降噪信號進(jìn)行EMD分解，利用頻譜比值法判斷虛假分量。此外，判斷虛假分量方法還有K-L散度法[14]，相關(guān)系數(shù)法[15-16]等。但對于虛假分量如何處理，目前關(guān)于這方面的研究并不多見。

本文以CEEMD為基礎(chǔ)，分解含有間斷信號和脈沖信號的仿真信號，分析并討論CEEMD中存在的模態(tài)混疊和虛假分量問題，提出了一種改進(jìn)的解相關(guān)CEEMD方法，有效地消除了仿真信號和爆破振動信號的模態(tài)混疊和虛假分量現(xiàn)象。

1 信號分解的CEEMD方法

1.1 CEEMD基本原理

CEEMD是針對EMD和EEMD (Ensemble Empirical Model Decomposition)在抑制信號分解中模態(tài)混疊問題的不足提出的改進(jìn)方法。由于EMD僅通過信號局部的極值來不斷的構(gòu)造上下包絡(luò)線，因此即使非常小的信號波動(如幅值很小的間歇成分)，也會造成EMD分解產(chǎn)生新的極值點，甚至新的IMF。這樣產(chǎn)生的IMF即不滿足唯一性。采用EMD分解無法解決模式混合問題。Wu等[17]對不同長度的白噪聲做了大量的EMD分解，并提出了噪聲輔助EMD分解的方式，即EEMD。EEMD可以在一定程度上抑制模態(tài)混疊的現(xiàn)象。EEMD是解決模式混合的一大創(chuàng)新，然而越來越多的研究指出EEMD存在著不足，比如殘留噪聲的疊加會導(dǎo)致嚴(yán)重的重構(gòu)誤差。為此，Yeh等[18]提出了噪聲互補的改進(jìn)算法CEEMD，主要步驟為：

對x(t)添加白噪聲時，是成對的添加組，即

(1)

(2)

最終得到的每個模態(tài)分量為

(3)

IMFj是對x(t)做CEEMD分解所得到的固態(tài)模式函數(shù)，它一般應(yīng)該有M=log2N-1個分量，即c1(t)、c2(t)、…、cM(t)，N是x(t)離散化以后的長度。

1.2 CEEMD的模態(tài)混疊和虛假分量現(xiàn)象

雖然CEEMD是一種較好的改進(jìn)算法，可以很好的抑制信號分解過程中存在的模態(tài)混疊現(xiàn)象。但根據(jù)以往的文獻(xiàn)和作者構(gòu)造的仿真信號分解結(jié)果表明，當(dāng)信號內(nèi)噪聲種類較多或信號突變點較多時，CEEMD的分解效果并不理想，同時在分解的過程中可能產(chǎn)生多個虛假分量。構(gòu)造仿真信號D(t)說明CEEMD的分解效果。

D(t)=cos(2π·70t)+0.3cos(2π·30t)+

v1(t)+u(t)

(4)

式中：v1(t)是200 Hz高頻間歇振蕩信號，振蕩區(qū)間為t=0.1～0.15 s、0.3～0.35 s，其余時間范圍幅值為0，兩個余弦信號頻率分別為70 Hz和30 Hz。u(t)為一脈沖信號，表達(dá)式如式(5)

(5)

CEEMD分解中IMF=ceemd (Y,Nstd,NE,TNM)，其中Nstd是噪聲標(biāo)準(zhǔn)偏差，取0.18，NE是集總次數(shù)，研究確定最佳集總平均次數(shù)為300次。TNM是分解獲得的分量個數(shù)，取12。得到的IMF分量見圖1。共分解出12個分量，僅分析前10個主分量，最后一個是余項rt(n)。

為抑制端點效應(yīng)，CEEMD分解中采用了張迪[19]提出的適用于爆破振動信號的端點處理方法。從圖1看到，各IMF分量在端點均未出現(xiàn)大的擾動。CEEMD將信號按頻率從高到低依次分解出來，但分解過程中出現(xiàn)了模態(tài)混疊和過分解現(xiàn)象。IMF3和IMF4均為間歇振蕩信號；IMF5為70 Hz的余弦分量，分解效果較好，信號內(nèi)未出現(xiàn)模態(tài)混疊現(xiàn)象；IMF6為30 Hz的余弦分量，其中出現(xiàn)了模態(tài)混疊現(xiàn)象，對IMF6求頻譜發(fā)現(xiàn)其混入了脈沖信號成分，且能量較小。由于分解過程中出現(xiàn)能量泄露，IMF5和IMF6的幅值比原始信號要小。IMF7和IMF8為脈沖分量，IMF9是否為脈沖分量需要做進(jìn)一步的判斷。

圖1 D(t)經(jīng)CEEMD分解得到的前10個IMF分量

對于易造成模態(tài)混疊的間歇信號和脈沖信號，容易出現(xiàn)兩個甚至多個虛假分量，需要對各IMF分量進(jìn)行虛假分量的判斷。采用求取分量與原始分量的互相關(guān)系數(shù)和頻譜判斷圖1中的虛假分量。IMF3、IMF4與原始分量的互相關(guān)系數(shù)分別為0.62和0.753，IMF7、IMF8和IMF9與原始分量的互相關(guān)系數(shù)分別為0.881、0.546和0.158。求取5個IMF分量的頻譜，IMF3和IMF4的主頻率均為200 Hz，并混有少量的70 Hz的成分，且IMF3的頻譜幅值較?。籌MF7和IMF8的主頻率均在10～13 Hz，IMF8的頻譜幅值較低，而IMF9的主頻率為4 Hz。判斷可知，IMF4和IMF7和原始單分量最為接近，為主要IMF分量，IMF3和IMF8與原始分量存在一定程度相關(guān)性，但能量相對較小，判斷這兩個分量為“偽”虛假分量。IMF9從互相關(guān)系數(shù)和主頻率上和原始分量都不同，為真實虛假分量。

1.3 CEEMD中虛假分量的處理

目前沒有虛假分量的統(tǒng)一處理方法，文中對分解得到的間歇信號和脈沖信號的虛假分量分別進(jìn)行合成與直接剔除處理，并與原始單分量對比分析，合成結(jié)果及其頻譜如圖2和圖3所示。對比結(jié)果如表1。

Fig.2 The original intermittency signal, mixture intermittency signal, intermittency signal by decomposed and

Fig.3 The original pulse signal, mixture pulse signal, pulse signal by decomposed and their frequency spectrum

從圖2和圖3和表1看出，間歇振蕩信號的頻率較高，最先被分離，因此合成后的IMF分量的幅值和原始單分量相差不大，能量損失較小。而由于端點效應(yīng)和間歇信號內(nèi)的奇異點的干擾，低頻IMF的分解誤差越來越大導(dǎo)致IMF7和IMF8的合成分量的能量損失較多，最大幅值約為原始分量的2/3。將“偽”虛假分量與主分量合并更有利于消除模態(tài)混疊，而且減輕了能量損失。

基于以上分析，由于爆破振動信號本身也是一種脈沖信號，而且含有多種高頻噪聲，CEEMD的分解效果不會特別理想。

2 改進(jìn)的解相關(guān)CEEMD方法

上述分析可知，當(dāng)信號內(nèi)噪聲種類和信號突變點較多時，僅用CEEMD分解還是會產(chǎn)生模態(tài)混疊和過分解現(xiàn)象，本質(zhì)原因是分解得到的IMF不滿足EMD的原始定義，即各個IMF沒有完全正交，存在信息耦合問題。從以下兩方面解決IMF正交性問題。

(1) 減少原始信號的奇異點、噪聲的干擾，即對信號進(jìn)行去噪處理。

(2) 將分解得到的IMF分量進(jìn)行后處理，使得到的最終IMF分量滿足正交性定義。

肖瑛等[20]提出了對EMD分解后進(jìn)行解相關(guān)操作抑制各分量間的模態(tài)混疊問題。研究發(fā)現(xiàn)，對頻率接近的平穩(wěn)信號，采用解相關(guān)的EMD算法可以較好抑制模態(tài)混疊現(xiàn)象。但在分解含脈沖噪聲和間歇信號噪聲時會出現(xiàn)嚴(yán)重的模態(tài)混疊問題。對解相關(guān)的EMD方法進(jìn)行改進(jìn)，使其能更有效的分離出爆破振動信號中的間歇振蕩信號和脈沖信號成分。

運用解相關(guān)方法處理CEEMD分解得到的各IMF分量的算法如下。

信號經(jīng)EMD(CEEMD)分解后，得到n個模態(tài)分量dn(t)(n=1,2,…)。任意兩相鄰分量的相關(guān)系數(shù)rc為(以d1(t)和d2(t)為例)

(6)

式中:上標(biāo)T代表信號的轉(zhuǎn)置。rc·d2(t)為d1(t)和d2(t)中的共同部分，將這部分從d1(t)中剔除，即得到與d2(t)中不相關(guān)的部分dv1(t)

dv1(t)=d1(t)-rcd2(t)

(7)

由HHT分解特性可知，IMF分量的頻率和幅值都隨著階數(shù)的增大逐漸減小。因此通過式(6)得到的rc不保證恒小于1，必須對兩個相鄰分量做歸一化處理，即每個分量內(nèi)所有數(shù)據(jù)均除以數(shù)列中的最大值，得到歸一化分量。歸一化后的相關(guān)系數(shù)rn與rc的關(guān)系為

(8)

n1和n2分別為d1(t)和d2(t)的歸一化系數(shù)，歸一化后得到的|rn|≤1，rn值越小，說明兩個分量獨立性越強；反之則表明分量間的共同部分越多。在歸一化條件下，式(7)改寫為

dv1(t)=d1(t)-rnd2(t)

(9)

若不對信號進(jìn)行端點處理，由于端點處的飛翼現(xiàn)象，信號最大值可能出現(xiàn)在端點附近，此時歸一化處理得到的歸一化信號與實際信號的誤差相差較大，因此必須對信號進(jìn)行端點處理后再進(jìn)行歸一化計算。采用張迪提出的適合于爆破振動信號的端點處理方法。

基于上述算法，提出了改進(jìn)的解相關(guān)CEEMD方法，步驟如下：

(1) 對原始信號端點處理，并進(jìn)行CEEMD分解。

(2) 利用互相關(guān)系數(shù)和頻譜判斷虛假分量，將虛假分量和主分量相加，得到新的合成分量。

(3) 將噪聲分量和含突變成分的前幾階分量去除，通過互相關(guān)系數(shù)法判斷第一個主分量，從這一個分量開始依次設(shè)為d1(t)、d2(t)、d3(t)…。

(4) 求取d1(t)和d2(t)的最大幅值并做歸一化處理，保證相關(guān)系數(shù)|rc|≤1。根據(jù)式(6)對d1(t)和d2(t)做解相關(guān)，求取兩個分量的相關(guān)系數(shù)，并根據(jù)式(9)得到d1(t)與d2(t)的不相關(guān)分量dv1(t)。

(5) 得到的解相關(guān)后的信號余量dr(t)為

dr(t)=d1(t)+d2(t)-dv(t)

(10)

對dr(t)再次進(jìn)行CEEMD分解，根據(jù)步驟(2)和(3)得到新的模態(tài)分量d11(t)、d22(t)、d33(t)…。將dv1(t)與d22(t)繼續(xù)做解相關(guān)處理，得到不相關(guān)分量dv2(t)。

(6) 重復(fù)(2)～(5)，直到混疊部分消除為止。其余的IMF分量參照前兩個分量的計算方法逐一分離。

由于CEEMD存在能量泄露，得到的IMF分量小于原始分量，因此分解次數(shù)并不是越多越好。針對不同信號類型和模態(tài)混疊程度，選擇合理的計算次數(shù)。由于CEEMD的分解效果比EMD有著較大的提升，不需要對每個IMF分量都進(jìn)行解相關(guān)計算。對得到的各個IMF分量求取頻譜，僅對包含多個頻帶成分的IMF分量進(jìn)行解相關(guān)操作，有利于對工程信號處理的實時性。

3 仿真分析

采用改進(jìn)的解相關(guān)CEEMD方法對D(t)求解，觀察其對模態(tài)混疊和虛假分量問題的處理效果，得到的4個IMF分量，采用傅里葉變換求取頻譜見圖4。

(a) IMF分量

(b) 頻譜

間歇振蕩信號的振蕩部分是由高頻余弦分量構(gòu)成，本質(zhì)上也是平穩(wěn)信號，因此使用解相關(guān)計算后的分解效果較好，從頻譜圖中看到，IMF1和IMF2內(nèi)已經(jīng)沒有模態(tài)混疊現(xiàn)象，由于使用的解相關(guān)次數(shù)較少，兩個信號的能量泄露也較小。對于IMF3和IMF4，脈沖信號屬于非平穩(wěn)信號，分解較困難。兩個分量經(jīng)多次解相關(guān)和CEEMD計算后，能量泄露較IMF1和IMF2多，而且IMF3和IMF4內(nèi)還存在少量模態(tài)混疊現(xiàn)象，但模態(tài)混疊部分能量較小，在頻譜圖中幾乎無法觀察到。對比圖4和圖1，經(jīng)解相關(guān)計算后得到的30 Hz余弦分量和脈沖分量與僅通過CEEMD分解得到的兩個分量相比，模態(tài)混疊現(xiàn)象得到極大改善，能量泄露也較小，因此信號經(jīng)CEEMD分解后再進(jìn)行解相關(guān)操作，分解的效果最好。

為了更直觀的評價改進(jìn)的解相關(guān)計算對消除模態(tài)混疊和虛假分量問題的效果，對D(t)分別進(jìn)行EMD、CEEMD和解相關(guān)CEEMD方法的分解，根據(jù)EMD、EEMD和CEEMD定義，得到的各IMF分量不可避免的存在一定程度的能量泄露，計算各IMF分量與原始單分量的能量差值也可以判斷抑制效果。對分解結(jié)果和原始分量計算相關(guān)系數(shù)ρi和能量均方根差值γm。

原始信號X(t)的均方根值SO

(11)

式中，N為X(t)的采樣點數(shù)。若X(t)經(jīng)CEEMD分解得到k個IMF，趨勢項rn(t)為第k+1個IMF分量，分解后各分量的均方根值SR為

(12)

ci(t)為分解得到的各階IMF。定義γm為能量均方根差值，即分解前后的均方根差值比值

(13)

γm值越小，說明得到的IMF分量越接近原始信號，分解效果越好；反之效果越差。計算結(jié)果見表2，其中D-CEEMD指端點處理、解相關(guān)計算和CEEMD聯(lián)合處理方法。

表2 IMF分量和原始分量的相關(guān)系數(shù)和能量均方根差值

Tab.2 The correlation coefficient and the difference value of energy root mean square between IMFs components and the original components

評價方法模態(tài)分量EMDCEEMDD-CEEMDIMF10.0890.760.929相關(guān)IMF20.9240.9990.998系數(shù)IMF30.6170.8890.985IMF40.7380.8730.978IMF10.860.0940.033能量IMF20.030.1260.116差值法IMF30.720.0980.122IMF40.160.5030.103

EMD分解中，4個信號主要分布在后3個IMF分量中，計算的ρi和γm沒有規(guī)律性。從CEEMD和改進(jìn)的解相關(guān)CEEMD方法得到的ρi和γm看到，改進(jìn)方法得到的4個分量的精度明顯高于只用CEEMD分解的計算精度，相關(guān)系數(shù)均在0.92以上，且利用能量均方根差值得到的結(jié)果均比CEEMD的小。從能量差值上也可以看到，兩種方法得到的IMF的γm值逐漸增大，這是由于分解次數(shù)的增加，能量泄露也逐漸增大導(dǎo)致。

4 爆破振動信號分析

圖5是深圳某地鐵隧道爆破振動信號，采用普通毫秒延期雷管起爆。圖中看出，爆破前期，由于圍巖夾制作用和較差自由面條件，質(zhì)點振動速度往往較大，前期段位微差間隔較小，相鄰段位之間波形會發(fā)生混疊，隨后期段位間隔增大，段與段之間基本可看成獨立爆破。

使用CEEMD方法和改進(jìn)的解相關(guān)CEEMD方法處理上述爆破振動信號，CEEMD分解得到12個IMF分量，列出前10個分量，最后一個為趨勢項。對CEEMD方法和改進(jìn)的CEEMD方法處理得到的波形使用傅里葉變換求頻譜圖，結(jié)果如圖6和圖7所示。

圖5 典型的爆破振動信號

如圖6，前兩個IMF分量頻率最高、但所占能量最小，為現(xiàn)場環(huán)境中高頻噪聲，IMF3～I(xiàn)MF7分量頻率逐漸降低，但相比IMF1，幅值很高，包含信號大部分能量，為原始信號的主頻段。頻譜圖中看出IMF3～I(xiàn)MF7之間模態(tài)混疊現(xiàn)象最嚴(yán)重，求取各分量與原始信號相關(guān)系數(shù)，IMF3～I(xiàn)MF10分別為0.76、0.79、0.39、0.24、0.15、0. 01、0.001 5、0.003 6。對比相關(guān)系數(shù)和頻譜圖可見IMF3～I(xiàn)MF7為主要分量，IMF8為IMF7的“偽”虛假分量，IMF9為真實虛假分量。使用改進(jìn)的解相關(guān)CEEMD方法，將虛假分量與主分量相加并對各分量進(jìn)行處理后，剔除噪聲分量和趨勢項，共得到5個IMF分量。如圖7，可以清楚看出，信號內(nèi)按頻率從高到低主要包含100 Hz、77 Hz、62 Hz、43 Hz、25 Hz、11 Hz和6 Hz共7種頻率成分。除了100 Hz和62 Hz的成分未分量出來外，其它5個IMF中的主頻率能量均較高，非頻率成分得到較好抑制。改進(jìn)方法得到的各IMF振速峰值分別為0.11，0.05，0.021，0.013，0.007。主頻基本分布在100 Hz以下，較之CEEMD得到的分量峰值(0.1，0.04，0.018，0.008 7，0.004 7)均有一定程度提高，對比頻譜圖中主頻對應(yīng)的峰值亦是如此，可見改進(jìn)的CEEMD方法在處理模態(tài)混疊的同時可以減輕一部分能量損失。

由于爆破振動波在傳播過程中的時頻特性與爆源條件、傳播介質(zhì)、地形等因素密切相關(guān)，波形包含重要的爆破參數(shù)信息。使用改進(jìn)的CEEMD方法處理爆破振動信號可以有效提取精確的IMF分量的時頻特征，實現(xiàn)爆破振動信號的精確分析處理，這對分析爆破振動波的作用機理，傳播規(guī)律、振動控制研究有著重要意義。

5 結(jié) 論

(1) 針對CEEMD在分解爆破振動信號存在的模態(tài)混疊和虛假分量問題，提出了一種基于端點處理，CEEMD和解相關(guān)的改進(jìn)方法。

(2) 采用CEEMD對信號進(jìn)行預(yù)處理，可以較好地分離出信號的全部主要分量，通過解相關(guān)計算可以進(jìn)一步解決存在混疊的IMF分量之間的耦合問題，將模態(tài)混疊成分分離。

(3) 采用改進(jìn)的解相關(guān)CEEMD方法，避免了直接剔除IMF虛假分量或者與主分量相加處理造成的能量損失。解決了以往僅能判定處虛假分量而無法對其后處理的問題。

(4) CEEMD仿真結(jié)果表明，改進(jìn)的解相關(guān)CEEMD方法的分解精度比EMD和CEEMD高，處理信號內(nèi)模態(tài)混疊和虛假分量的效果較好。

(5) 將改進(jìn)的解相關(guān)CEEMD方法其應(yīng)用于爆破振動信號，基本保證了各IMF分量中非主頻率信號的能量比例較低，有效抑制了每個IMF內(nèi)的模態(tài)混疊問題。