黃維新， 劉敦文

(中南大學(xué) 資源與安全工程學(xué)院，長沙 410083)

微震監(jiān)測作為一種先進(jìn)的地壓監(jiān)測手段，在礦山領(lǐng)域得到了越來越廣泛的應(yīng)用[1-4]。微震監(jiān)測借助埋設(shè)的傳感器接收巖體破裂等微震信號，再通過處理微震信號展開后續(xù)工作，因而微震信號的質(zhì)量直接關(guān)系到后續(xù)工作的有效性。然而，礦井工作環(huán)境復(fù)雜，微震信號常與汽車運(yùn)輸、風(fēng)機(jī)振動、鑿巖沖擊、放礦等噪音混合，這些噪音將對P波、S波初至拾取、震源定位和震源機(jī)制反演等造成影響。因此，開展微震信號去噪的研究具有重要意義。

微震信號具有隨機(jī)性、非平穩(wěn)的特點(diǎn)，基于傅里葉變換的降噪對周期信號效果較好，而對含尖峰和突變的微震信號適應(yīng)性較差。近年來，小波和經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)在微震信號降噪中的應(yīng)用越來越廣泛。例如：徐宏斌等[5-7]采用小波閥值對微震信號進(jìn)行了降噪，而曹思遠(yuǎn)等[8]和Mousavi等[9]分別借助二代小波變換+閥值、同步壓縮小波變換+閥值對微震信號去噪；梁喆等[10]和秦晅等[11]、Shang等[12]分別采用EMD與互信息熵、EMD與閥值對微震信號去噪。然而，基于小波變換的降噪，需選取合適的基函數(shù)和閾值才能達(dá)到較好的降噪效果；EMD是一種基于經(jīng)驗性的分解，其存在模態(tài)混疊、虛假分量等缺陷。梁喆等[13]采用基于LMS(Least Mean Square)自適應(yīng)濾波器的LEMD(Empirical Mode Decomposition Based on LMS)方法對微震信號去噪，降低了EMD的端點(diǎn)效應(yīng)。Han等[14]采用集合經(jīng)驗?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition, EEMD)和閥值對微震信號去噪，降低了EMD模態(tài)混疊的影響。

VMD(Variational Mode Decomposition)作為一種新的信號分解方法，其具有嚴(yán)格的數(shù)學(xué)推導(dǎo)，能降低EMD模態(tài)混疊、虛假分量的影響[15-16]。然而，由于微震信號的復(fù)雜性，VMD得到的部分模態(tài)分量仍為噪音與有用信號的混合，直接舍棄或不做處理均不能取得較好的降噪效果。獨(dú)立成分分析(Independent Component Analysis, ICA)能根據(jù)多個觀測信號分離出相互獨(dú)立的原始信號，然而其觀測信號數(shù)須大于等于分離出的獨(dú)立信號數(shù)目。為此，可將噪音與有用信號混合的多個模態(tài)分量作為ICA的觀測信號，進(jìn)而分離出有用信號，再與低頻模態(tài)分量相加得到VMD_ICA降噪信號。此外，礦山微震信號常含有工頻噪音，本文采用了一種正弦函數(shù)擬合的方法去除工頻噪音。

1 降噪方法

1.1 變分模態(tài)分解

VMD理論假定待分解信號x(t)由有限個帶寬有限、中心頻率不同的模態(tài)分量組成。在各模態(tài)分量之和等于待分解信號的約束下，以求取各模態(tài)分量的估計帶寬之和最小而構(gòu)建的變分模型。與EMD得到的模態(tài)分量不同，VMD的模態(tài)分量均為調(diào)幅-調(diào)頻信號，即

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

(1)

式中：uk(t)為第k個模態(tài)分量；Ak(t)和φk(t)分別為uk(t)的瞬時幅值和相位。

VMD變分模型的建立過程如圖1所示，其核心思想是構(gòu)建式(2)所示的變分模型，并采用式(3)所示的Lagrange方程求解該變分模型。

圖1 VMD工作原理流程圖Fig.1 Principle flow diagram of the VMD method

(2)

(3)

式中：wk為uk的頻率中心；δ(t)為狄拉克函數(shù)；α為懲罰因子；λ為Lagrange乘法算子。

1.2 獨(dú)立成分分析

1.2.1 ICA

獨(dú)立成分分析從觀測信號出發(fā)，借助信號的高階統(tǒng)計特性，對混合信號進(jìn)行估計分析，進(jìn)而得到相互獨(dú)立的近似源信號。假設(shè)M個觀測信號X=[x1,x2, …,xM]由N個獨(dú)立信號S=[s1,s2, …,sN]線性組合而成，則有

X=AS

(4)

式中：A為M×N階矩陣，M≥N。

1.2.2 FastICA

本文采用Hyvarinen[17]提出的FastICA分離混合信號，其主要步驟如下：

步驟1 對觀測信號中心化，使其均值為0；進(jìn)而白化，得到標(biāo)準(zhǔn)化后的數(shù)據(jù)z；

步驟2 設(shè)定獨(dú)立成分的個數(shù)N和收斂閥值ε；

步驟3 令分離矩陣W=[w1,w2, …,wM]T，初始化wi使W得的模為1；

步驟5 對矩陣W正交化,W←(WWT)-1/2W；

步驟6 判斷W是否收斂，若1-min{abs[diag(W(k+1)T)W(k)]}小于收斂閥值ε，則W為所求分離矩陣；反之，則繼續(xù)步驟3～步驟6；

步驟7 由公式Y(jié)=WX求取獨(dú)立分量。

1.3 VMD_ICA_Sine去噪

VMD能將信號分解至不同中心頻帶，最后幾個模態(tài)分量通常為高頻噪音，可采用相關(guān)系數(shù)法直接去除；中間的模態(tài)分量為有用信號與噪音的混合，采用ICA分離噪音和有用信號。此外，VMD_ICA不能去除工頻噪音，提出采用正弦函數(shù)擬合去噪。相關(guān)系數(shù)和Sine函數(shù)定義分別為

(5)

x′(t)=A0sin(w0t+φ0)+c

(6)

本文將上述VMD_ICA和Sine函數(shù)聯(lián)合降噪簡稱為VMD_ICA_Sine去噪，其流程如圖2所示，具體步驟如下：

步驟1 導(dǎo)入待去噪信號x(t)，t=1, 2, …,n, 其中，n為信號點(diǎn)的個數(shù)；

步驟2 采用VMD將x(t)分解為不同個數(shù)模態(tài)分量，再用主頻差別法確定最佳模態(tài)分量個數(shù)；

步驟3 采用式(5)計算各模態(tài)分量與x(t)的相關(guān)系數(shù)，去掉相關(guān)系數(shù)較小的模態(tài)分量；

步驟4 采用相關(guān)系數(shù)適中的模態(tài)分量構(gòu)建矩陣，再用FastICA分離噪音和有用信號；

步驟5 步驟4中提取的有用信號與剩余相關(guān)系數(shù)較大的模態(tài)分量相加，得到VMD_ICA去噪信號；

步驟6 采用式(6)擬合步驟5去噪信號的工頻噪音，再用步驟5去噪信號減去正弦擬合信號得到VMD_ICA_Sine去噪信號。

圖2 VMD_ICA_Sine去噪流程圖Fig. 2 Flow diagram of the VMD_ICA_Sine denoising method

2 微震信號測試

2.1 VMD和ICA應(yīng)用

以某礦微震信號為例(圖3(a)中最上側(cè)信號)，信號采樣頻率為6 000 Hz。采用VMD將該信號分解為2～8個模態(tài)分量，每個模態(tài)分量的主頻如表1所示。由表1知：主頻隨模態(tài)分量序號增加而增大，模態(tài)數(shù)為7和8時，u2的主頻不再下降。為此，本文確定將該微震信號分解為6個模態(tài)分量。

表1 不同模態(tài)數(shù)對應(yīng)的主頻

VMD分解得到的6個模態(tài)分量及其FFT幅值譜，如圖3所示。采用式(5)計算得到u1～u6與原信號的相關(guān)系數(shù)分別為0.83,0.77,0.73,0.71,0.58和0.32，可見u6與原信號的相關(guān)系數(shù)明顯小于其他模態(tài)分量的相關(guān)系數(shù)。此外，由圖3易得u6幾乎為隨機(jī)噪音，而u2～u5仍含有部分有用信號(圖中圈出部分及不太明顯的尾部信號)。由此，直接舍掉u6，對u2～u5進(jìn)一步提取特征。

將u2～u5作為觀測信號，進(jìn)而構(gòu)建FastICA特征矩陣X，F(xiàn)astICA分析得到的獨(dú)立分量IC1～I(xiàn)C4及其對應(yīng)幅值譜，如圖4所示。由圖4可知，F(xiàn)astICA有效的將噪音與有用信號分離開來，提取到了P波初至信息特征。再將IC4與矩陣A對應(yīng)的分量相乘得到具有實際幅值的有用信號，并與u1相加得到VMD_ICA降噪信號。

圖3 微震信號VMD分解結(jié)果及其幅值譜Fig.3 VMD decomposition results of the microseismic signal and their corresponding amplitude spectra

圖4 ICA獨(dú)立分量及其幅值譜Fig.4 Independent modes obtained from the ICA and their corresponding amplitude spectra

2.2 去噪效果

為測試VMD_ICA_Sine降噪方法的優(yōu)越性，本文將舎棄噪音模態(tài)分量u6和噪音與有用信號混合模態(tài)分量u2～u5、舍棄噪音模態(tài)分量u6、VMD_ICA法和VMD_ICA_Sine法去噪信號分別繪于圖5中。為定量的評價降噪效果，引入信噪比(Signal to Noise Ratio, SNR)，其定義為[18]

SNR=10 lg(σsignal/σnoise)

(7)

圖5 不同方法降噪效果及其幅值譜Fig.5 Denoising performance of different methods and their corresponding amplitude spectra

由圖5可知，不同降噪方法對SNR均有一定的提升，其中僅去除u6的降噪效果最差，仍含有大量高頻噪音；舎棄u2～u6雖然得到了較好的去噪效果，但其相較于VMD_ICA法去噪，其P波初至等局部特征稍差；而VMD_ICA_Sine法有效的去除了工頻噪音，同時很好的保留了P波初至信息。可見，VMD_ICA_Sine法具有較好的降噪效果。

3 微震應(yīng)用

從某礦山微震監(jiān)測系統(tǒng)中隨機(jī)抽取1 938個微震信號，信號采樣頻率為6 000 Hz。采用PAI-K法拾取P波初至殘差來評價降噪效果，PAI-K法借助P波初至前，微震信號主要為噪音，其高斯性很強(qiáng)，峰度值趨于零；而當(dāng)P波初至?xí)r，非高斯性增強(qiáng)，取峰度值最大點(diǎn)作為P波初至。峰度值計算式[19]

(8)

從上述1 938個微震信號中選取4個典型信號作為分析，其原始波形、VMD_ICA降噪波形和VMD_ICA_Sine降噪波形及其對應(yīng)的峰度值序列，如圖6所示，并將PAI-K拾取結(jié)果與原始波形人工拾取相比較，得到的拾取殘差統(tǒng)計于表2中。

由圖表知：PAI-K法對高信噪比信號(見圖6(d))拾取效果均較好，但其可能受噪音和尾部信號影響(見圖6(a)～圖6(c))；VMD_ICA法能有效去除高頻噪音和毛刺噪音，但其未能去除工頻噪音(見圖6(a))；VMD_ICA_Sine法可有效的去除工頻噪音(見圖6(a))，且其不受低頻信號干擾(見圖6(b)和圖6(c))。

表2 典型微震信號PAI-K法拾取結(jié)果

圖6 典型微震信號PAI-K法拾取效果Fig.6 Picking performance for typical microseismic signals using the PAI-K method

將上述1 938個原微震信號及降噪微震信號的信噪比繪于圖7中，圖中SNR1，SNR2和SNR3分別指原微震信號、VMD_ICA降噪信號和VMD_ICA_Sine降噪信號的信噪比。由圖7可知，VMD_ICA法和VMD_ICA_Sine法有效的提高了原信號的信噪比，且VMD_ICA_Sine法的信噪比優(yōu)于VMD_ICA法。需要指出的是VMD_ICA_Sine法降噪與微震信號中含工頻噪音的比例有關(guān)，筆者統(tǒng)計得到這1 938個微震信號中264個微震信號含有較明顯的工頻噪音。倘若含工頻噪音微震信號的比例增加，VMD_ICA_Sine降噪效果會更好。

上述1 938個微震信號拾取殘差如圖8所示，拾取殘差指PAI-K法拾取減去人工拾取。由圖8可知，VMD_ICA法和VMD_ICA_Sine法降噪信號拾取殘差為0～ 5 ms的微震信號數(shù)目明顯高于原微震信號拾取數(shù)目，而拾取殘差大于30 ms的微震信號數(shù)目小于原微震信號拾取數(shù)目，可見降噪后拾取效果優(yōu)于原信號拾取效果。

進(jìn)一步地，本文采用Li等[20]提出的定量方法評價P波初至拾取效果。該方法認(rèn)為信號拾取誤差越小，則懲罰值越小，并對所有信號的懲罰值求和，總懲罰值越小對應(yīng)的方法越優(yōu)。總懲罰值和懲罰值計算公式分別為

圖7 原信號與降噪信號信噪比比較Fig.7 SNR comparison between original signals and denoised signals

圖8 原信號與降噪信號的PAI-K法拾取殘差統(tǒng)計圖Fig.8 Picking residual statistics for original signals and denoised signals using the PAI-K method

(9)

(10)

式中：Ci為微震信號i拾取的懲罰值，大括號右側(cè)第1列為懲罰值，第2列為拾取誤差范圍。

計算得到原微震信號、VMD_ICA法和VMD_ICA_Sine法降噪信號拾取誤差的總懲罰值TCF分別為611.5,528.9和465.2。由此可知，VMD_ICA_Sine法降噪信號拾取效果優(yōu)于VMD_ICA法降噪拾取，且VMD_ICA降噪拾取優(yōu)于原信號拾取，即VMD_ICA_Sine法能為礦山微震信號降噪提供了一種較好的分析方法。

4 結(jié) 論

為降低礦山微震信號噪音，本文提出了一種基于VMD_ICA_Sine的降噪方法，并將其用于微震信號測試和應(yīng)用中，主要結(jié)論如下：

(1) 相關(guān)系數(shù)法可剔除高頻噪音模態(tài)分量，且ICA能從噪音與有用信號混合的模態(tài)分量中提取有用信號。

(2) VMD_ICA法和VMD_ICA_Sine法能有效保留微震信號的局部特征，且其信噪比大于直接去除部分模態(tài)分量。

(3) VMD_ICA和VMD_ICA_Sine法均有效地提高了PAI-K法P波初至拾取效果，且VMD_ICA_Sine法優(yōu)于VMD_ICA法降噪效果。