基于STA/LTA改進的CEEMD-SVD微震信號降噪算法

史艷楠， 齊朋磊， 王 裕， 王毅穎,2， 張沖沖

(1.河北工程大學 機械與裝備工程學院，河北 邯鄲 056038；2.河北省煤炭生態保護開采產業技術研究院，河北 邯鄲 056038；3.邯鄲市智能車輛重點實驗室，河北 邯鄲 056038)

近年來，隨著我國大型水利水電工程不斷建設及礦產資源開采進入深部，深部工程造成的巖石爆裂等地質災害問題日益顯著，嚴重影響了高效施工及安全生產。巖爆是深部工程開挖卸荷過程中，巖石中儲存的彈性能突然釋放導致巖體結構發生變化的現象[1-2]，因此獲取巖爆發生的信息并掌握其規律是十分必要的。微震監測通過對巖爆過程中產生的信號進行檢測和分析，可以實現對微震事件的超前預測與預報，被廣泛應用于深埋隧道、大型水電站、礦井深部開采等深部工程[3-4]。因此，迅速、精確的識別巖石爆破所形成的微震信號尤為重要。但深部工程環境條件惡劣，檢波器采集到的微震信號中往往會含有大量噪聲信號，這些噪聲與有效微震信號混疊在一起，會極大的干擾到微震信號的識別、初至拾取和震源定位。微震信號降噪方法研究對提高微震事件監測預警的準確性具有十分重要的現實意義。

目前，信號降噪主要以采用基于傅里葉變換的方法為主，但這種方法主要針對平穩信號具有良好的降噪效果。由于采集到的微震信號為非線性時變信號，采用傳統方法時降噪效果較差，針對此問題文獻[5]提出將卡爾曼濾波的方法應用于微震信號的降噪中，解決了傳統降噪方法處理非線性時變信號效果差的問題。同時，隨著小波技術的發展，利用其多分辨率特性，自適應的對低頻部分和高頻部分的信號進行辨別，成為有效處理非平穩信號局部特征的工具[6-7]。文獻[8-9]將小波技術應用于微震信號的預處理，對采集到的信號進行了很好的濾波和降噪，證明小波技術可以實現微震信號的降噪；文獻[10]將小波閾值算法應用于巖體結構監測的降噪，不僅提取出了低信噪比的有效微震信號，而且還實現了對頻率范圍廣的巖體結構微震信號中有用信號和噪聲信號的分離。但小波算法的基函數及閾值等參數需要根據工作環境確定，增加了算法調試時間，適應性不好。因此，國內外許多學者將Hilber-Huang(希爾伯特黃)的EMD應用到微震信號的降噪中，解決了小波算法預設參數的問題，該方法具有良好的時頻特性，能夠處理非線性、非平穩信號，通過將微震信號分解為不同的IMF，更能反映信號本身固有的屬性，自適應好[11-12]。文獻[13-14]利用EMD算法的自適應性，將振動信號分解成IMF，根據IMF固有屬性，重組振動信號，提高了信號的信噪比，獲得了較為明顯的降噪效果。文獻[15]采用EMD分解與獨立成分分析結合的方法，優化了對IMF分量的選取，并對分界分量進一步去噪，充分保留有效微震信號的特征。在使用EMD降噪時，一個IMF分量中存在不同時間尺度的特征，出現模態混疊的現象，降低去噪效果。文獻[16]提出了基于EEMD的算法，該算法加入均勻分布的噪聲信號，使原信號頻率分布趨于均勻，有效解決了模態混疊的現象，卻無法在多次分解中得到的相同的固有模態函數的數量，需要更多次數的計算，效率下降。文獻[17]通過計算EEMD各IMF分量的STA/LTA值來選取適合的IMF分量重構信號。分量的STA/LTA值可以反應信號的部分特征，可以有效的去除噪聲分量。文獻[18]提出EEMD-Hankel-SVD的方法對微震信號進行降噪，通過對高頻分量的去除和低頻分量的降噪，有效保留了信號原始特征，提高了信噪比和初至拾取效果。文獻[19]針對傳統消噪方法的不足，提出改進型CEEMD-SVD消噪方法對高邊坡變形數據進行處理，消噪效果良好，該方法在一定程度上可以反映原始信號中的細節信息，但會造成有用信息損失。

基于以上分析，本文提出了一種基于STA/LTA改進的CEEMD-SVD微震信號降噪方法，本文算法采用將正負白噪聲加入原始信號后，進行EMD分解，計算各IMF分量的STA/LTA值，通過閾值來判定是否對其進行SVD分解。集總平均IMFi，得到降噪后信號。CEEMD是對EEMD進一步改進，集總平均次數會減少，從EEMD的百量級減小到CEEMD的幾十的量級，重建后的信號噪聲明顯減少。本文方法提高了信號預處理的效率，減小了重構誤差，有效的保留了原始信號的特征，很好的壓制了噪聲信號，為信號后續研究提供了高信噪比的數據，為微震活動特征及規律性研究分析奠定基礎。

1 基本理論

1.1 CEEMD

CEEMD算法是對EMD算法的進一步優化，在原始信號中加入n對不同的正負白噪聲，用于抵消信號中加入的白噪聲。對兩個信號分別進行EMD模態分解。求得各IMF分量的均值，再對n組IMF集總平均得到重構信號。具體算法流程如圖1所示。

圖1 CEEMD流程圖Fig.1 The flow chart of CEEMD

算法的實現步驟如下：

在原微震信號s(t)中分別加入和減去n個均值為0、幅值和標準差為常數的高斯白噪聲n(t)，得到新的重構信號

xia(t)=si(t)+ni(t)

(1)

xib(t)=si(t)-ni(t)

(2)

式中:xia(t)為加入第i個高斯白噪聲信號后的復合信號;xib(t)為減去第i個高斯白噪聲信號后的復合信號;ni(t)為加入的第i個高斯白噪聲信號。

1.2 SVD

SVD算法是通過把時間序列信號矩陣化分解的方法，將一維信號轉換為矩陣進行運算分解，而構建矩陣方算法眾多，主要有Toeplitz矩陣和Cycle矩陣及Hankel矩陣，使用不同的矩陣構建方法，SVD分解處理信號的效果是不一樣的。本文采用構建Hankel矩陣，然后進行SVD分解時，得到一系列的奇異值，每個奇異值對應不同的頻率幅值，想要提取的特征信號，只需提找出對應的奇異值，并把剩余奇異值歸0即可。根據這一特性，通過加入與原信號相關的頻率約束信號，放大該頻率的k個奇異值，將剩余奇異值歸0，再進行奇異值分解逆變換，并將加入的約束信號去除，即可提取出信號中有用信息。構建Hankel矩陣進行SVD分解的過程如下：

對于信號x={x1,x2,…,xi}，i=1,2,…,M?？梢詷嫿ㄈ缦翲ankel矩陣

(3)

式中，m+n-1=M,m≥n。根據文獻[20]可知，當m=M/2時，使用Hankel矩陣進行SVD降噪效果最佳。

對D進行奇異值分解，得到

(4)

式中:A、B為正交矩陣;AT為A的轉置;Σ為對角方陣，其對角元素由奇異值組成，且奇異值存在如下關系：β1>β2>…>βk>0，式中βk為奇異值分解得到的第k個奇異值，xi、yi是矩陣A、B的第i個列向量，r為矩陣D的秩。

2 基于STA/LTA的CEEMD-SVD算法

CEEMD算法與EMD算法有一些相同的特征，在其IMF分量中，在排列靠前的幾個分量是高頻分量，包含大量的隨機噪聲，通過STA/LTA閾值可以去除大量的含燥分量。而舍棄的IMF分量中仍包含有用的信號，直接舍去會影響到重構之后的信號波形，無法完全體現原有的特征信息，相應評價指標無法達到最佳。若對含噪聲信號直接采用傳統SVD降噪，則對奇異值的選取困難信號頻率范圍達不到理想狀態。因此采用基于STA/LTA的CEEMD-SVD方法進行降噪，選擇不符合閾值要求的IMF分量進行SVD分解降噪，以充分保留高頻IMF成分中的有用信息，降低了有用信息的損失，使重構后波形與各項指標達到最優化。算法流程圖如圖2所示。

圖2 算法流程圖Fig.2 Algorithm flow chart

降噪實現步驟如下：

(1) 對原含噪聲信號進行CEEMD分解，得到最終分解的一系列IMF分量。

(2) 設定STA/LTA算法的閾值ε，ε的取值根據實際信號以及時窗長度進行選取[21]。

(3) 計算各IMF分量的STA/LTA值；舍去不符合要求的分量。

(4) 將保留的IMF分量的和作為頻率約束條件，加入到舍棄的IMF分量中，分別進行SVD分解，提取有用信號。

(5) 將提取到的有用信號與最終殘余分量求和，即得到降噪后的信號。

3 仿真實驗分析

為了驗證本文提出方法的有效性和降噪效果，采用實測微震信號通過LabVIEW平臺進行試驗，并與EMD、EEMD和SVD降噪效果進行對比、分析。

3.1 降噪效果評價

為了更好的說明算法的有效性，本文采用信噪比(SNR)、能量保持百分比(E)以及信號標準差(RMSE)作為評價指標。其定義

(5)

E=Ex/Ey

(6)

(7)

式中:y(n)為降噪后信號;x(n)為有用信號;N為采樣點數;Ex為有用信號能量;Ey為含噪聲信號能量。

3.2 仿真分析

為驗證本文算法降噪效果，在原有用信號中加入噪聲進行仿真模擬。本文采用某礦采樣時間段長度為4 000個采樣點的微震信號進行模擬仿真分析[22]。圖3(a)為原微震信號的幅頻特性圖3(b)含5 dB噪聲微震信號的幅頻特性。

(a) 原信號的波形及幅頻特性

(b) 含噪微震信號的波形及幅頻特性圖3 原信號及含噪微震信號的波形及幅頻特性Fig.3 Waveform and amplitude-frequency characteristics of the original signal and the noisy microseismic signa

對圖3(b)含燥微震信號進行CEEMD分解，得到11個IMF分量，如圖4所示。對圖3(a)原信號的幅頻特性分析可知不含噪聲信號的頻率主要咋0～500 Hz之間，從圖4(a)和圖4(b)中可以看出IMF1、IMF2的波形噪聲污染較為嚴重，頻率分布寬泛與原信號相差較大。IFM3以后的分量波形較為平滑，振動趨勢與原信號相似，且頻率范圍主要在0～500 Hz之間，集中在原信號頻率范圍附近。

(a) IMF1的波形及幅頻特性

(b) IMF2的波形及幅頻特性

(c) IMF3的波形及幅頻特性

(d) IMF4的波形及幅頻特性

(e) IMF5的波形及幅頻特性

(f) IMF6的波形及幅頻特性

(g) IMF7的波形及幅頻特性

(h) IMF8的波形及幅頻特性

(i) IMF9的波形及幅頻特性

(j) IMF10的波形及幅頻特性

(k) IMF11的波形及幅頻特性圖4 各IMF分量的波形及幅頻特性Fig.4 Waveform and amplitude-frequency characteristics of each IMF

為了解決CEEMD算法產生的偽分量現象，通過結合STA/LTA算法，可以去除包含大量噪聲信號的IMF分量。求各分量的STA/LTA值如圖5所示，通過分析長短時窗的長度，本文設定閾值為2，重構信號預選取IMF2～IMF7及IMF9等7個分量為重構信號的子頻段信號。

(a) IMF2

(b) IMF3

(c) IMF4

(d) IMF5

(e) IMF6

(f) IMF7

(g) IMF9圖5 各IMF分量的STA/LTA值Fig.5 STA/LTA values of each IMF

計算各IMF分量與原微震信號的相關性如圖6所示，可知MF3為分界IMF，進而確定信號主導的IMF分量為IMF3～IMF11，與STA/LTA算法相結合選取IMF3～IMF7重構信號。

圖6 相關系數Fig.6 Correlation coefficient

對預選的IMF分量中的噪聲主導的IMF1，IMF2進行SVD降噪。最后將噪聲主導的IMF降噪后與IMF3～IMF7及殘余分量進行信號重構，得到降噪信號，如圖7所示。

(a)

(b)圖7 重構微震信號波形及幅頻幅頻特性Fig.7 Reconstructed microseismic signal waveform and amplitude-frequency amplitude-frequency characteristics

通過本文算法得到的重構信號見圖7，重構信號的頻率范圍與原信號頻率接近。從表1可以看出重構后的信號信噪比、能量百分比及信號標準差較加噪信號都用明顯提升。說明本文方法可以有效的降除噪聲。

表1 降噪效果分析Tab.1 Analysis of noise reduction effect

3.3 實驗分析

以某礦的微震監測系統采集到的微震信號為分析對象，對本文方法的降噪效果進行驗證，實測含噪微震信號的波形圖及幅頻特性如圖8所示。

(a)

(b)圖8 含噪微震信號的波形及幅頻特性Fig.8 Waveform and amplitude-frequency characteristics of noisy microseismic signals

從圖8可以看出微震信號頻率在0～100 Hz時幅值較大，有用的信號頻率主要分布在該頻率范圍內，但信號中包含大量隨機噪聲，有用微震信號波形受到嚴重干擾。采用本文方法對該信號進行處理，同時對信號使用NAEEMD、EEMD和EMD三種方法進行降噪效果對比分析，采用各方法對微震信號降噪后的波形及幅頻特性分別如圖9所示，其中圖9(a)為本文方法降噪，圖9(b)為EMD分解降噪，圖9(c)為EEMD分解降噪，圖9(d)為NAEEMD分解降噪。

(a) 本文

(b) EMD

(c) EEMD

(d) NAEEMD圖9 微震信號各方法降噪后波形及幅頻特性Fig.9 Waveform and amplitude-frequency characteristics of microseismic signal after denoising by various methods

從圖9可以看出，兩段微震信號波形通過各降噪方法降噪后比降噪前波形更加清晰，且較好的保持了波形的變化趨勢，幅頻特性的優勢頻率分布更加明顯。圖9(a)中降噪后的波形平滑性較好，而圖9(b)中的波形平滑性更好，但其波峰波谷的削弱較大，圖9(c)和圖9(d)的雖然波峰波谷沒有明顯削弱，但波形的平滑度較差；同時圖9中頻率范圍集中在0～100 Hz，與圖9(b)、(c)和(d)對比優勢明顯。

通過對各方法降噪前后的信噪比、能量百分比及信號標準差進行計算，結果如表2所示。從表2可以看出，微震信號經過4種方法降噪后，本文方法較其他方法的信噪比均高出5 dB以上，且信號的能量百分比高達95.586 7%，信號標準差達到了0.014 8，優于其他降噪方法。因此，應用本文方法對微震信號進行降噪不僅是有效的，而且較其他方法優勢明顯，具有更好的應用效果。

表2 微震信號降噪對比Tab.2 Noise reduction comparison of microseismic signals

4 結 論

(1) 對采集到的微震信號經過CEEMD分解，得到固有模態分量，利用STA/LTA確定保留的IMF分量和舍棄的IMF分量，然后對舍棄的IMF分量降噪，有效的獲取IMF分量中信號主導分量，使得重構后的微震信號較好的保留了微震信號的波形特征，同時在CEEMD分解時通過對殘差信號不斷加入白噪聲，并只進行一階模態分解，因此，每次加噪后的誤差不是通過獨立分解得到的，而是由最后一次分解的誤差得到的，很大程度上減少了分解后的重構誤差。

(2) 通過模擬噪聲信號進行仿真實驗，驗證了頻率約束的可行性，經過對降噪后信號的信噪比、能量百分比及信號標準差計算，分別從5.013 32 dB，75.402%及0.034 5達到了15.670 1 dB，95.227%及0.021 9，信號降噪效果明顯，能有效的去除信號中存在的隨機噪聲，證明了本文方法對微震信號進行降噪是有效的。

(3) 采用本文方法對某礦實測微震信號進行降噪分析，并與EMD，EEMD，NAEEMD降噪方法進行對比，結果表明本文方法可以有效的降低煤礦微震信號中的噪聲干擾，且較其他方法優勢明顯，具有更好的應用效果。但仍需在深埋隧道、大型水電站等工程環境背景條件下進行進一步驗證與優化。