爆破振動信號分析中模態混疊和虛假分量消除的改進方法

2019-09-17 06:46:32徐文龍馮丹丹

振動與沖擊 2019年17期

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

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

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

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

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

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

1 信號分解的CEEMD方法

1.1 CEEMD基本原理

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

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

(1)

(2)

最終得到的每個模態分量為

(3)

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

1.2 CEEMD的模態混疊和虛假分量現象

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

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

v1(t)+u(t)

(4)

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

(5)

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

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

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

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

1.3 CEEMD中虛假分量的處理

目前沒有虛假分量的統一處理方法，文中對分解得到的間歇信號和脈沖信號的虛假分量分別進行合成與直接剔除處理，并與原始單分量對比分析，合成結果及其頻譜如圖2和圖3所示。對比結果如表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分量的幅值和原始單分量相差不大，能量損失較小。而由于端點效應和間歇信號內的奇異點的干擾，低頻IMF的分解誤差越來越大導致IMF7和IMF8的合成分量的能量損失較多，最大幅值約為原始分量的2/3。將“偽”虛假分量與主分量合并更有利于消除模態混疊，而且減輕了能量損失。

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

2 改進的解相關CEEMD方法

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

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

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

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

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

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

(6)

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

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

(7)

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

(8)

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

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

(9)

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

基于上述算法，提出了改進的解相關CEEMD方法，步驟如下：

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

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

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

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

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

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

(10)

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

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

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

3 仿真分析

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

(a) IMF分量

(b) 頻譜

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

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

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

(11)

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

(12)

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

(13)

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

表2 IMF分量和原始分量的相關系數和能量均方根差值

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

評價方法模態分量EMDCEEMDD-CEEMDIMF10.0890.760.929相關IMF20.9240.9990.998系數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沒有規律性。從CEEMD和改進的解相關CEEMD方法得到的ρi和γm看到，改進方法得到的4個分量的精度明顯高于只用CEEMD分解的計算精度，相關系數均在0.92以上，且利用能量均方根差值得到的結果均比CEEMD的小。從能量差值上也可以看到，兩種方法得到的IMF的γm值逐漸增大，這是由于分解次數的增加，能量泄露也逐漸增大導致。

4 爆破振動信號分析

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

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

圖5 典型的爆破振動信號

如圖6，前兩個IMF分量頻率最高、但所占能量最小，為現場環境中高頻噪聲，IMF3～IMF7分量頻率逐漸降低，但相比IMF1，幅值很高，包含信號大部分能量，為原始信號的主頻段。頻譜圖中看出IMF3～IMF7之間模態混疊現象最嚴重，求取各分量與原始信號相關系數，IMF3～IMF10分別為0.76、0.79、0.39、0.24、0.15、0. 01、0.001 5、0.003 6。對比相關系數和頻譜圖可見IMF3～IMF7為主要分量，IMF8為IMF7的“偽”虛假分量，IMF9為真實虛假分量。使用改進的解相關CEEMD方法，將虛假分量與主分量相加并對各分量進行處理后，剔除噪聲分量和趨勢項，共得到5個IMF分量。如圖7，可以清楚看出，信號內按頻率從高到低主要包含100 Hz、77 Hz、62 Hz、43 Hz、25 Hz、11 Hz和6 Hz共7種頻率成分。除了100 Hz和62 Hz的成分未分量出來外，其它5個IMF中的主頻率能量均較高，非頻率成分得到較好抑制。改進方法得到的各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)均有一定程度提高，對比頻譜圖中主頻對應的峰值亦是如此，可見改進的CEEMD方法在處理模態混疊的同時可以減輕一部分能量損失。