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

彭亞雄， 劉廣進， 蘇 瑩， 陳春暉， 劉運思

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

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

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

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

1 VMD分解原理

VMD算法是在信號頻域內利用迭代計算獲得變分模態最優解，變換模態函數和中心頻率，得到不同帶寬的本征模態函數。將信號模態估計變為變分問題進行求解，基本步驟如下[12]。

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

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

(1)

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

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

(2)

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

(3)

2 自適應VMD-MPE降噪方法

2.1 自適應VMD算法

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

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

(6)

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

不同模態數K的IMF能量和與原信號能量不同，為了確定最適合模態數K，定義了重構信號與原信號能量差參數ξ的計算式

(7)

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

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

(8)

2.2 MPE算法

MPE算法是基于排列熵(permutation entropy, PE)的改進算法，可用于檢測信號的隨機和動力突變特性。將時間序列進行多尺度粗?；儆嬎鉖E。具體步驟如下[19]

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

(9)

(10)

(13)

2.3 自適應VMD-MPE降噪

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

圖1 自適應VMD-MPE降噪流程Fig.1 Denoising process with adaptive VMD-MPE

3 工程應用

3.1 信號降噪處理

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

表1 臺階爆破參數Tab.1 Bench blasting parameters

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

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

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

表2 不同模態數K分解的能量差參數ξTab.2 Energy difference parameter ξ with different modes number K

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

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

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

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

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

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

3.2 時頻能量分析

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

(a) 原信號

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

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

4 降噪效果對比分析

4.1 評價指標

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

(1) 重構標準差ESD

(14)

(2) 均方根誤差ε

(15)

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

4.2 降噪效果對比

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

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

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

表4 爆破地震波重構和降噪效果指標Tab.4 Reconstruct and denoised effect index of blasting seismic wave signals

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

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

5 結 論

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

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

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

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