基于自適應(yīng)VMD-MPE算法的礦山爆破地震波信號降噪方法研究

彭亞雄， 劉廣進(jìn)， 蘇 瑩， 陳春暉， 劉運思

(1. 湖南科技大學(xué) 巖土工程穩(wěn)定控制與健康監(jiān)測湖南省重點實驗室， 湖南 湘潭 411201；2. 中國地質(zhì)大學(xué)(武漢) 工程學(xué)院， 武漢 430074)

爆破工程中采用實測信號對爆破有害效應(yīng)進(jìn)行詳細(xì)分析與評價。然而受礦山環(huán)境復(fù)雜性、監(jiān)測設(shè)備電磁干擾、巖體介質(zhì)的反射及折射等原因的影響，大量高頻噪聲導(dǎo)致信號波形畸變，掩蓋了真實信號成分[1-2]。為了掌握爆破地震波波形、強度和頻率特征，必須通過信號處理技術(shù)對實測信號進(jìn)行降噪處理。

爆破地震波信號屬于典型的非平穩(wěn)隨機信號，目前針對這類信號的降噪算法主要包括小波類算法[3-4]和經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)算法[5-6]。在爆破信號降噪研究領(lǐng)域，熊正明等[7]通過多次使用平移不變小波對爆破信號進(jìn)行平移并進(jìn)行閾值降噪，實現(xiàn)了擺脫小波基的信號降噪處理。路亮等[8]提出了利用提升小波包變換的最優(yōu)基搜索改進(jìn)算法，有效去除爆破振動信號噪聲。黃智剛等[9]提出了基于CEEMDAN-MPE算法的隧道爆破地震波信號降噪方法，不僅能夠去除高頻噪聲，對信號所含主要信息影響較小。Peng等[10]針對水下鉆孔爆破振動信號所含噪聲特點，提出了CEEMDAN光滑降噪模型，有效降低了實測信號噪聲含量。這些降噪算法的提出與工程應(yīng)用，充分說明了信號處理技術(shù)對實測爆破信號降噪處理的有效性。然而，小波算法在降噪過程中小波基函數(shù)和分解層次難以確定，使得這類方法的自適應(yīng)性不強，降噪效果難以保證；EMD及其改算法具有較好的自適性，但其端點效應(yīng)和模態(tài)混疊問題無法有效解決，使得這兩類算法在處理爆破地震波信號的魯棒性較低。

變分模態(tài)分解(variational mode decomposition, VMD)是一種非遞歸的信號分解方法，迭代搜尋變分的最優(yōu)解，以確定信號各部分頻率中心及帶寬，完成頻域劃分和分量分離。該方法避免了端點效應(yīng)和模態(tài)混疊問題[11]，對各類信號具有普遍適應(yīng)性和魯棒性，但是存在模態(tài)數(shù)K難以確定的問題。因此，本文提出了基于能量差參數(shù)ξ的自適應(yīng)VMD分解方法，將分解得到的本征模態(tài)函數(shù)(intrinsic mode function, IMF)進(jìn)行多尺度排列熵(multi-scale permutation entropy, MPE)隨機性檢測，去除噪聲IMF達(dá)到降噪目的，并應(yīng)用于實測礦山爆破地震波信號降噪處理。

1 VMD分解原理

VMD算法是在信號頻域內(nèi)利用迭代計算獲得變分模態(tài)最優(yōu)解，變換模態(tài)函數(shù)和中心頻率，得到不同帶寬的本征模態(tài)函數(shù)。將信號模態(tài)估計變?yōu)樽兎謫栴}進(jìn)行求解，基本步驟如下[12]。

將原信號x(t)分解為K個中心頻率為ωk的模態(tài)函數(shù)uk，其中K為模態(tài)數(shù)，需人為設(shè)定[13]。模態(tài)函數(shù)表示為

uk(t)=Ak(t)cos(φk(t))

(1)

式中:uk(t)為第k個分量；Ak(t)和φk(t)分別為瞬時幅值和相位。

對uk(t)進(jìn)行Hibert變換，得到解析信號和單邊頻譜ψ；解析信號加入估計中心頻率，將模態(tài)頻譜轉(zhuǎn)換至基頻帶；計算解析信號梯度的平方L2范數(shù)，估計模態(tài)信號帶寬，成為約束性變分問題，可以求解出x(t)的IMF

(2)

式中:{uk}為信號分解得到的K個IMF分量；{ωk}為各分量對應(yīng)的中心頻率；δ(t)為脈沖函數(shù)。

(3)

2 自適應(yīng)VMD-MPE降噪方法

2.1 自適應(yīng)VMD算法

VMD算法中必須人為設(shè)定模態(tài)數(shù)K(即IMF分量個數(shù))和懲罰因子α。當(dāng)模態(tài)數(shù)K過少，將導(dǎo)致原信號缺失；過多則會產(chǎn)生頻率混疊、過分解等問題。當(dāng)懲罰因子α較小時，IMF帶寬較大易產(chǎn)生混疊現(xiàn)象，當(dāng)α較大則影響收斂速度[14]。

通過觀察不同模態(tài)數(shù)K分解得到的各IMF中心頻率，可以選出最合適的模態(tài)數(shù)K。然而該方法在判斷中心頻率變化規(guī)律時，受到人為因素影響較大[15]。為此，本文提出了基于信號能量的自適應(yīng)VMD算法。爆破地震波信號能量定義[16]為

(6)

式中:E為地震波信號能量;Es為原信號能量；第k個IMF能量表示為Ek；v(t)為信號序列；t為信號采集長度。

不同模態(tài)數(shù)K的IMF能量和與原信號能量不同，為了確定最適合模態(tài)數(shù)K，定義了重構(gòu)信號與原信號能量差參數(shù)ξ的計算式

(7)

能量差參數(shù)ξ隨著模態(tài)數(shù)K增加而增加；信號過分解前，ξ的增長幅度很小，而當(dāng)信號出現(xiàn)分解時，ξ將大幅增長。因此，由K=2,3,…進(jìn)行信號分解，計算能量差參數(shù)ξK；根據(jù)文獻(xiàn)[17-18]，當(dāng)|ξK+1-ξK|>0.1時，認(rèn)為模態(tài)數(shù)K+1分解的信號處于過分解狀態(tài)，則確定最合適模態(tài)數(shù)為K。

在模態(tài)數(shù)K值確定后，分別計算不同懲罰因子α條件下信號VMD分解的信噪比SNR

(8)

2.2 MPE算法

MPE算法是基于排列熵(permutation entropy, PE)的改進(jìn)算法，可用于檢測信號的隨機和動力突變特性。將時間序列進(jìn)行多尺度粗粒化，再計算PE。具體步驟如下[19]

① 時間序列x={x1,x2,…,xL}的多尺度粗?；?/p>

(9)

(10)

(13)

2.3 自適應(yīng)VMD-MPE降噪

自適應(yīng)VMD-MPE降噪方法是對原信號進(jìn)行自適應(yīng)VMD分解，計算最適合的懲罰因子α和模態(tài)數(shù)K，進(jìn)而獲得K個IMF分量，其中IMF1～I(xiàn)MFK的主頻逐漸增加。然后對各IMF分量進(jìn)行MPE隨機性檢測，利用IMF分量的MPE均值判斷是否為噪聲成分。對信號主成分IMF分量進(jìn)行重構(gòu)，得到降噪信號。對于爆破地震波信號，真實信號成分的MPE閾值為0.6[20]。自適應(yīng)VMD-MPE降噪流程如圖1所示。

圖1 自適應(yīng)VMD-MPE降噪流程Fig.1 Denoising process with adaptive VMD-MPE

3 工程應(yīng)用

3.1 信號降噪處理

邦山脈鐵礦位于利比里亞中部，邦州西南地區(qū)，年開采量約為1 000萬噸。礦區(qū)內(nèi)巖層整體向北傾斜，傾角較陡65°～85°，主要巖性為含鐵豐富片麻巖、片巖和石英巖，巖層完整，弱風(fēng)化和中風(fēng)化，裂隙不發(fā)育。礦山開采采用中深孔臺階爆破方式，爆破參數(shù)如表1所示。

表1 臺階爆破參數(shù)Tab.1 Bench blasting parameters

礦區(qū)范圍內(nèi)存在多個村莊，為評價礦山爆破對周圍環(huán)境的影響，采用M20型測振儀進(jìn)行爆破地震波監(jiān)測。選取一條典型爆破地震波信號為研究對象(圖2)。該實測信號的采樣頻率為4 000 sps，根據(jù)Nyquist采樣定理，實測信號頻率為2 000 Hz；采用時間為2 s，共采集8 000個采樣點。

圖2 爆破地震波信號(s1)Fig.2 Blasting seismic wave signal (s1)

自適應(yīng)VMD算法的主要設(shè)定參數(shù)包括：模態(tài)數(shù)K、懲罰因子α、保真度τ和收斂條件，針對不同類型信號需要設(shè)置合適參數(shù)值。對于隨機非平穩(wěn)信號，保真度τ和收斂條件ε采用推薦值[21]；計算不同模態(tài)數(shù)K分解的能量差參數(shù)ξ和|ξK+1-ξK|如表2所示。根據(jù)自適應(yīng)VMD算法的判斷依據(jù)，該爆破地震波信號的最適合模態(tài)數(shù)K為6，其IMF分量如圖3所示。然后通過計算不同懲罰因子α條件下信號VMD分解的信噪比SNR，得到適合本工程α=2 000，減少了分解過程的細(xì)節(jié)特征損失。

表2 不同模態(tài)數(shù)K分解的能量差參數(shù)ξTab.2 Energy difference parameter ξ with different modes number K

由表2和圖3可知，采用VMD算法分解得到實測信號共6個IMF分量，IMF1～I(xiàn)MF6的中心頻率逐漸增加。由IMF分量的波形和中心頻率變化，可以推斷IMF4～I(xiàn)MF6的幅值較小且頻率較高可能是爆破地震波信號中的高頻噪聲分量，而IMF1～I(xiàn)MF3則為真實信號成分。

為準(zhǔn)確區(qū)分真實信號成分和噪聲，計算各IMF分量的MPE值。經(jīng)過多次試算確定合適的嵌入維數(shù)m=6，時間延遲τ=1，尺度因子s=5。得到各IMF分量的MPE平均值如表3所示。

表3 IMF分量的MPE均值Tab.3 Mean MPE of IMF

由表3可知，IMF1～I(xiàn)MF6的MPE均值逐漸增加，IMF分量中所含噪聲成分逐漸減少，說明自適應(yīng)VMD算法能夠很好地將真實信號與噪聲信號成分分離。對于爆破地震波信號，真實信號成分的MPE閾值為0.6，IMF4～I(xiàn)MF6的MPE均值大于閾值為噪聲信號，與上述波形分析結(jié)果一致。因此，將IMF4～I(xiàn)MF6從原信號中剔除，保留爆破地震波真實信號成分，達(dá)到降噪的目的。同時，自適應(yīng)算法避免了VMD分解過程中的人工干預(yù)，降低了人為誤差。

降噪后的爆破地震波信號如圖4所示。與原信號波形(圖2)相比，二者的波形一致，噪聲得到明顯消除，且對爆破地震波峰值振速影響較小。

圖4 降噪后的爆破地震波信號Fig.4 Blasting seismic wave signal after denoising

3.2 時頻能量分析

自適應(yīng)最優(yōu)核(adaptive optimal kernel, AOK)時頻技術(shù)采用了短時模糊函數(shù)和時變自適應(yīng)最優(yōu)核函數(shù)，能夠獲得更多的細(xì)節(jié)信息，較好反映了爆破地震波信號的時頻特性。為了進(jìn)一步驗證自適應(yīng)VMD-MPE算法的降噪效果，對原信號和降噪信號進(jìn)行AOK時頻譜對比，如圖5所示。圖中X為峰值能量，Y為信號主頻。

(a) 原信號

(b) 降噪信號圖5 AOK時頻譜Fig.5 AOK time-frequency spectrum

由圖5可知，實測礦山爆破地震波信號的時頻能量主要分布在頻率0～250 Hz和時間0.2～0.9 s范圍內(nèi)，信號主頻為10.71 Hz。原信號在頻率500 Hz以上含有明顯的噪聲成分，是導(dǎo)致信號失真的主要原因，如圖5(a)中虛線框中所示。利用自適應(yīng)VMD-MPE算法獲得的降噪信號，較好地剔除了原有噪聲成分；降噪信號與原信號的主頻一致，峰值能量和主頻帶能量(0～250 Hz)均沒有明顯降低。說明自適應(yīng)VMD-MPE算法不僅能有效剔除實測礦山爆破地震波信號中的高頻噪聲，且對低頻信號能量影響較小。

4 降噪效果對比分析

4.1 評價指標(biāo)

實測信號在VMD分解過程中不可避免導(dǎo)致原有信號成分的丟失，為反映原信號與重構(gòu)信號的差異性，定義了重構(gòu)標(biāo)準(zhǔn)差ESD作為評價信號丟失和失真情況的評價指標(biāo)。采用降噪信號與原信號的信噪比SNR、均方根誤差ε作為降噪效果評價指標(biāo)，分析實測礦山爆破地震波信號的降噪效果。重構(gòu)標(biāo)準(zhǔn)差ESD和均方根誤差ε定義如下

(1) 重構(gòu)標(biāo)準(zhǔn)差ESD

(14)

(2) 均方根誤差ε

(15)

此外，對比降噪前后信號波形特征，能夠定性判斷降噪效果，以確保特征波形的一致性和明顯噪點去除干凈。

4.2 降噪效果對比

為了進(jìn)一步分析自適應(yīng)VMD-MPE算法的降噪效果，對3組礦山爆破地震波信號(如圖6)進(jìn)行降噪處理，并與EEMD-MPE、CEEMDAN-MPE降噪算法進(jìn)行對比。三種算法獲得的降噪信號如圖6所示，并分別計算重構(gòu)標(biāo)準(zhǔn)差ESD、信噪比SNR、均方根誤差ε，如表4所示。

圖6 原信號與降噪信號對比Fig.6 Comparison between original signals and denoised signals

圖6 原信號與降噪信號對比Fig.6 Comparison between original signals and denoised signals

表4 爆破地震波重構(gòu)和降噪效果指標(biāo)Tab.4 Reconstruct and denoised effect index of blasting seismic wave signals

對比3組實測礦山爆破地震波信號的降噪結(jié)果，由波形曲線(圖6)可知，EEMD-MPE、CEEMDAN-MPE和自適應(yīng)VMD-MPE三種算法均取得了一定程度降噪效果，去除了部分高頻噪聲，其中EEMD-MPE算法的降噪信號仍然具有較為明顯的噪聲，降噪效果相對較差。由表3可知，在信號分解與重構(gòu)過程中，原信號與重構(gòu)信號的標(biāo)準(zhǔn)差ESD在0.008 3～0.020 3，均屬于較低范圍，說明三種算法的保真性均較好。然而由于EEMD分解過程中，以增加白噪聲方式消除模態(tài)混疊問題，導(dǎo)致了其重構(gòu)標(biāo)準(zhǔn)差ESD均最大。對比兩個降噪效果指標(biāo)，自適應(yīng)VMD-MPE算法的信噪比SNR均最大，CEEMDAN-MPE算法略小于EMD自適應(yīng)VMD-MPE算法；自適應(yīng)VMD-MPE算法的均方根誤差ε最低。

因此，自適應(yīng)VMD-MPE算法的降噪波形和降噪效果指標(biāo)均優(yōu)于其他兩種算法，能夠更好地去除高頻噪聲信號，同時能夠更好地保留信號中的有效信息，在礦山爆破地震波信號降噪處理上具有明顯的優(yōu)越性。

5 結(jié) 論

受到礦山復(fù)雜環(huán)境的影響，礦山爆破地震波信號所含高頻噪聲較多，導(dǎo)致真實信號內(nèi)容被噪聲掩蓋。針對這一問題提出了基于自適應(yīng)VMD-MPE算法降噪方法，用于礦山爆破地震波信號降噪處理。

(1) 對爆破地震波信號進(jìn)行VMD分解，得到不同頻率的IMF，有效避免了模態(tài)混疊與端點效應(yīng)問題；提出了基于能量差參數(shù)ξ確定模態(tài)數(shù)K的自適應(yīng)方法，實現(xiàn)了最優(yōu)分解，其結(jié)果比主觀決策更為可靠。

(2) 對實測礦山爆破地震波信號進(jìn)行處理表明，自適應(yīng)VMD分解可以獲得精細(xì)的IMF，MPE熵值識別出噪聲IMF，達(dá)到了較好去除降噪的目的。對比分析降噪前后振動信號的AOK時頻能量譜，自適應(yīng)VMD-MPE降噪方法不僅能成功地去除500 Hz以上的高頻噪聲能量，而且對低頻信號能量影響較小。

(3) 對比EEMD-MPE、CEEMDAN-MPE和自適應(yīng)VMD-MPE三種降噪算法，三種方法都具有較好的保真性和一定的降噪效果。自適應(yīng)VMD-MPE算法的降噪信號波形和降噪效果指標(biāo)均優(yōu)于EEMD-MPE、CEEMDAN-MPE算法，驗證了基于自適應(yīng)VMD-MPE算法的有效性，對礦山爆破地震波信號降噪和爆破有效效應(yīng)分析具有指導(dǎo)意義。