基于STA/LTA改進(jìn)的CEEMD-SVD微震信號(hào)降噪算法

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

(1.河北工程大學(xué) 機(jī)械與裝備工程學(xué)院，河北 邯鄲 056038；2.河北省煤炭生態(tài)保護(hù)開采產(chǎn)業(yè)技術(shù)研究院，河北 邯鄲 056038；3.邯鄲市智能車輛重點(diǎn)實(shí)驗(yàn)室，河北 邯鄲 056038)

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

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

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

1 基本理論

1.1 CEEMD

CEEMD算法是對(duì)EMD算法的進(jìn)一步優(yōu)化，在原始信號(hào)中加入n對(duì)不同的正負(fù)白噪聲，用于抵消信號(hào)中加入的白噪聲。對(duì)兩個(gè)信號(hào)分別進(jìn)行EMD模態(tài)分解。求得各IMF分量的均值，再對(duì)n組IMF集總平均得到重構(gòu)信號(hào)。具體算法流程如圖1所示。

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

算法的實(shí)現(xiàn)步驟如下：

在原微震信號(hào)s(t)中分別加入和減去n個(gè)均值為0、幅值和標(biāo)準(zhǔn)差為常數(shù)的高斯白噪聲n(t)，得到新的重構(gòu)信號(hào)

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

(1)

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

(2)

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

1.2 SVD

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

對(duì)于信號(hào)x={x1,x2,…,xi}，i=1,2,…,M?？梢詷?gòu)建如下Hankel矩陣

(3)

式中，m+n-1=M,m≥n。根據(jù)文獻(xiàn)[20]可知，當(dāng)m=M/2時(shí)，使用Hankel矩陣進(jìn)行SVD降噪效果最佳。

對(duì)D進(jìn)行奇異值分解，得到

(4)

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

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

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

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

降噪實(shí)現(xiàn)步驟如下：

(1) 對(duì)原含噪聲信號(hào)進(jìn)行CEEMD分解，得到最終分解的一系列IMF分量。

(2) 設(shè)定STA/LTA算法的閾值ε，ε的取值根據(jù)實(shí)際信號(hào)以及時(shí)窗長(zhǎng)度進(jìn)行選取[21]。

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

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

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

3 仿真實(shí)驗(yàn)分析

為了驗(yàn)證本文提出方法的有效性和降噪效果，采用實(shí)測(cè)微震信號(hào)通過(guò)LabVIEW平臺(tái)進(jìn)行試驗(yàn)，并與EMD、EEMD和SVD降噪效果進(jìn)行對(duì)比、分析。

3.1 降噪效果評(píng)價(jià)

為了更好的說(shuō)明算法的有效性，本文采用信噪比(SNR)、能量保持百分比(E)以及信號(hào)標(biāo)準(zhǔn)差(RMSE)作為評(píng)價(jià)指標(biāo)。其定義

(5)

E=Ex/Ey

(6)

(7)

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

3.2 仿真分析

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

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

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

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

(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算法產(chǎn)生的偽分量現(xiàn)象，通過(guò)結(jié)合STA/LTA算法，可以去除包含大量噪聲信號(hào)的IMF分量。求各分量的STA/LTA值如圖5所示，通過(guò)分析長(zhǎng)短時(shí)窗的長(zhǎng)度，本文設(shè)定閾值為2，重構(gòu)信號(hào)預(yù)選取IMF2～I(xiàn)MF7及IMF9等7個(gè)分量為重構(gòu)信號(hào)的子頻段信號(hào)。

(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

計(jì)算各IMF分量與原微震信號(hào)的相關(guān)性如圖6所示，可知MF3為分界IMF，進(jìn)而確定信號(hào)主導(dǎo)的IMF分量為IMF3～I(xiàn)MF11，與STA/LTA算法相結(jié)合選取IMF3～I(xiàn)MF7重構(gòu)信號(hào)。

圖6 相關(guān)系數(shù)Fig.6 Correlation coefficient

對(duì)預(yù)選的IMF分量中的噪聲主導(dǎo)的IMF1，IMF2進(jìn)行SVD降噪。最后將噪聲主導(dǎo)的IMF降噪后與IMF3～I(xiàn)MF7及殘余分量進(jìn)行信號(hào)重構(gòu)，得到降噪信號(hào)，如圖7所示。

(a)

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

通過(guò)本文算法得到的重構(gòu)信號(hào)見(jiàn)圖7，重構(gòu)信號(hào)的頻率范圍與原信號(hào)頻率接近。從表1可以看出重構(gòu)后的信號(hào)信噪比、能量百分比及信號(hào)標(biāo)準(zhǔn)差較加噪信號(hào)都用明顯提升。說(shuō)明本文方法可以有效的降除噪聲。

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

3.3 實(shí)驗(yàn)分析

以某礦的微震監(jiān)測(cè)系統(tǒng)采集到的微震信號(hào)為分析對(duì)象，對(duì)本文方法的降噪效果進(jìn)行驗(yàn)證，實(shí)測(cè)含噪微震信號(hào)的波形圖及幅頻特性如圖8所示。

(a)

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

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

(a) 本文

(b) EMD

(c) EEMD

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

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

通過(guò)對(duì)各方法降噪前后的信噪比、能量百分比及信號(hào)標(biāo)準(zhǔn)差進(jìn)行計(jì)算，結(jié)果如表2所示。從表2可以看出，微震信號(hào)經(jīng)過(guò)4種方法降噪后，本文方法較其他方法的信噪比均高出5 dB以上，且信號(hào)的能量百分比高達(dá)95.586 7%，信號(hào)標(biāo)準(zhǔn)差達(dá)到了0.014 8，優(yōu)于其他降噪方法。因此，應(yīng)用本文方法對(duì)微震信號(hào)進(jìn)行降噪不僅是有效的，而且較其他方法優(yōu)勢(shì)明顯，具有更好的應(yīng)用效果。

表2 微震信號(hào)降噪對(duì)比Tab.2 Noise reduction comparison of microseismic signals

4 結(jié) 論

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

(2) 通過(guò)模擬噪聲信號(hào)進(jìn)行仿真實(shí)驗(yàn)，驗(yàn)證了頻率約束的可行性，經(jīng)過(guò)對(duì)降噪后信號(hào)的信噪比、能量百分比及信號(hào)標(biāo)準(zhǔn)差計(jì)算，分別從5.013 32 dB，75.402%及0.034 5達(dá)到了15.670 1 dB，95.227%及0.021 9，信號(hào)降噪效果明顯，能有效的去除信號(hào)中存在的隨機(jī)噪聲，證明了本文方法對(duì)微震信號(hào)進(jìn)行降噪是有效的。

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