改進EWT_MPE模型在礦山微震信號特征提取中的應(yīng)用

程鐵棟，易其文，吳義文，戴聰聰，蔡改貧，楊麗榮，尹寶勇

(江西理工大學(xué)，江西 贛州 341000)

隨著深部礦產(chǎn)資源的進一步開采及地下空間的開發(fā)利用，沖擊地壓、頂板大面積垮落、礦震等地質(zhì)災(zāi)害越發(fā)頻繁,礦巖條件與開采環(huán)境日益惡化，礦井安全問題日益突出[1-3]。能夠提前監(jiān)測礦山地壓災(zāi)害可能發(fā)生的位置，是對深部巖土工程巖爆進行有效預(yù)測的重要途徑之一。微震監(jiān)測技術(shù)是目前較為先進和有效的地壓監(jiān)測手段，在國內(nèi)外金屬、煤巖礦山都有諸多應(yīng)用[4-8]。微震監(jiān)測將采集到的信號進行波形分析，從而獲取微震事件的時間、位置及震源參數(shù)信息[9]，然而傳感器所采集到的信號混雜了爆破、機械噪聲等干擾信號，目前國外具備自動識別功能的微震監(jiān)測系統(tǒng)，對巖體微震與爆破振動信號進行自動識別效果較差，一般還需要依靠人工識別，受個體因素影響可能存在誤判，且效率不高[10]。如何從諸多波形信號中提取出有效的微震信號用于巖爆預(yù)測顯得尤為重要。

時頻分析是微震與爆破信號波形識別的主要方法之一，能夠?qū)ξ⒄鸷捅菩盘柕念l譜特征進行有效提取和分析。曹業(yè)安等[11-14]采用傅里葉變換對礦山微震和爆破信號的頻譜特征進行了研究，為沖擊地壓的預(yù)測預(yù)報工作提供了一種思路。吳響等[15]在經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)的基礎(chǔ)上分別結(jié)合奇異值分解和Volterra級數(shù)對微震與爆破信號特征提取。朱權(quán)潔等[16]從能量分布的角度采用小波包分解對微震事件進行分析。唐守鋒等[17]提出了以一種基于小波理論的小波能譜系數(shù)分析方法，對煤巖破裂微震與噪聲信號進行分析，為后期模式識別打下了重要基礎(chǔ)。朱權(quán)潔等[18]結(jié)合小波分析與分形理論對微震信號進行分析，以分形盒維數(shù)作為特征向量進行特征提取，并通過支持向量機(support vector machines,SVM)進行分類識別，識別準確率高達94%。李偉提出基于局部均值分解(local mean decomposition LMD)和SVM對微震與爆破信號進行特征提取和分類，準確率達93%。以上方法雖然在微震監(jiān)測領(lǐng)域廣泛應(yīng)用，并取得了不錯的效果，但是也存在一些不足：傅里葉變換僅能通過頻譜分析對微震與爆破信號進行初步判定；傳統(tǒng)EMD缺乏有力的理論基礎(chǔ)，且在信號分解過程中易出現(xiàn)模態(tài)混疊和端點效應(yīng)，依靠經(jīng)驗選取有效IMF分量的過程中存在著誤判現(xiàn)象；小波分析和小波包分解對信號的自適應(yīng)分解性能較差。

經(jīng)驗小波變換(empirical wavelet transform,EWT)是由Gilles[19]提出的一種新的自適應(yīng)信號分解方法，該方法在小波分析理論框架下建立，理論基礎(chǔ)完備，適用于非線性、非平穩(wěn)信號分析。EWT不但可以避免EMD存在的模態(tài)混疊和端點效應(yīng)的發(fā)生，而且提高了信噪比、縮短了計算時間，目前已廣泛應(yīng)用于故障檢測[20-22]。然而，EWT在以往處理復(fù)雜信號頻譜時易出現(xiàn)過切分問題，導(dǎo)致分解效果不好，還需進一步優(yōu)化。多尺度排列熵(multi-scale permutation entropy,MPE)是基于排列熵基礎(chǔ)上改進的一種非線性分析算法，用于衡量時間序列在不同尺度下的復(fù)雜性和隨機性,具有穩(wěn)定性高、運算量小、抗噪能力強等優(yōu)點，目前已被廣泛應(yīng)用于故障診斷領(lǐng)域[23-25]，但MPE在微震領(lǐng)域未見應(yīng)用。鑒于此，本文嘗試利用EWT對巖體微震與爆破振動信號進行分解，并在傳統(tǒng)EWT的基礎(chǔ)上提出了新的改進方法，利用相關(guān)性分析篩選出最優(yōu)IMF分量進行重構(gòu)，再結(jié)合MPE提取特征向量，最后采用GK(Gustafson-Kessel)模糊聚類算法進行分類識別。

1 EWT算法改進

EWT的核心是通過對信號的傅里葉頻譜進行自適應(yīng)劃分，構(gòu)造合適的正交小波濾波器組，將單一信號分解并提取分解信號中不同頻率特征信息的模態(tài)分量。然而所采用的局部極大值分割法，對簡單信號頻譜進行劃分時效果較好，處理復(fù)雜信號時易出現(xiàn)過切分或欠切分問題，使得結(jié)果產(chǎn)生誤差。如何對傅里葉頻譜進行有效劃分極其重要，這一步?jīng)Q定了EWT對分解信號的自適應(yīng)效果。為此，本文采用最大熵譜估計法[26-27]對EWT進行改進，其主要步驟流程如圖1所示，其中虛線框為改進部分。

圖1 改進EWT算法流程圖Fig.1 Improved EWT algorithm flowchart

2 多尺度排列熵(MPE)

多尺度排列熵是Aizz等[28]在排列熵基礎(chǔ)上改進的一種非線性分析算法，其基本思想是將時間序列首先進行多尺度粗粒化，再計算各尺度子序列的排列熵，具體步驟如下：

(1)

式中，s為尺度因子，s=1,2,...；[N/s]表示對N/s取整。

(2)

其中l(wèi)為第l個重構(gòu)分量l=1,2,…,N-(m-1)τ；τ為延遲時間；m為嵌入維數(shù)。

將重構(gòu)序列按升序排列，可得到符號序列S(r)=(l1,l2,…,lm)。其中r=1,2,…,R且R≤m!。計算每一種符號序列出現(xiàn)的概率Pr。

根據(jù)等式(3)計算每個粗?；蛄械呐帕徐兀瑲w一化后得到時間序列在多尺度下的排列熵值。

(3)

3 基于改進EWT_MPE的信號特征提取及GK分類方法

3.1 基于改進EWT_MPE的信號特征提取

由于改進EWT分解得到的IMF分量中，包含了如信號的局部特征、信號不同時間尺度特征等諸多信息，一些分量對微震信號的特征提取貢獻不大，為此，采用分量與原始信號的相關(guān)系數(shù)對IMF分量進行篩選。相關(guān)系數(shù)越大，其對應(yīng)的分量就越重要[29]。其定義如下：

(i=1,2,…,n)

(4)

式中:co(i)為各分量與原始信號x(i)的相關(guān)系數(shù);N為微震事件的總采樣點數(shù)。

3.2 GK分類方法

GK模糊聚類算法[30]是利用協(xié)方差矩陣的自適應(yīng)距離進行度量，通過最小化目標函數(shù)來獲得隸屬度矩陣U=[μij]c×n和聚類中心V=[v1,v2,…,vn]Τ。其中vi(i=1,2,…,c)為第i個聚類中心，n為樣本個數(shù)，μij表示j個元素屬于第i類隸屬度，且滿足：

(5)

GK模糊聚類的目標函數(shù)為

(6)

(7)

(8)

(9)

GK模糊聚類的聚類效果可用分類系數(shù)和平均模糊熵進行檢驗。

分類系數(shù)

(10)

平均模糊熵

(11)

其中，C越接近1，分類效果越清晰，反之分類越模糊；E的值越接近0，聚類效果越好。

4 仿真信號分析

為驗證本文所提算法的可行性和準確性，先采用仿真信號進行對比分析試驗。由于原信號經(jīng)過EWT分解可得到N+1個IMF分量，那么原信號也可看成由多個分辨率不同的信號分量疊加而成[31]。為此，仿真信號fsig(t)取單調(diào)趨勢信號f1(t)，調(diào)頻信號f2(t)和兩個不同時間段頻率單一的信號f3(t)組成。另外，考慮到礦山微震信號采集過程中存在噪聲干擾，仿真信號中添加了白噪聲n(t)，信噪比為3，t∈[0,1]。fsig(t)時域波形如圖2所示。

圖2 仿真信號fsig(t)的時域波形Fig.2 Time domain waveform of the simulated signal fsig(t)

經(jīng)驗小波算法分解所輸出的結(jié)果，實質(zhì)是由一個尺度函數(shù)和N個小波濾波函數(shù)濾波所得。試驗取N=4對信號頻譜進行劃分。圖3(a)和(b)分別為改進EWT和傳統(tǒng)EWT頻譜劃分結(jié)果及其對應(yīng)的小波濾波器組；圖4(a)和(b)是fsig(t)信號的改進EWT和傳統(tǒng)EWT分解結(jié)果。

(a) 改進EWT

(a) 改進EWT

由圖4分解結(jié)果可知,傳統(tǒng)EWT和改進EWT均能將噪聲及4個分量分解出來，其中改進EWT能較好地將各模態(tài)分量分解出來，且與原始信號非常吻合。相比之下，圖4(b)所示傳統(tǒng)EWT所分解出的各模態(tài)分量吻合度并不高，IMF2～IMF4分量更是出現(xiàn)畸形。這是由于傳統(tǒng)EWT對傅里葉頻譜劃分時，受噪聲和非平穩(wěn)信號分量影響產(chǎn)生偽峰和毛刺[32](圖3(b)所示)，引起傅里葉譜波動而不平滑，從而導(dǎo)致所檢測到的部分邊界結(jié)果出現(xiàn)錯誤，由此產(chǎn)生的無用極值被保留在峰值序列中，而一些有用極值未被保留在峰值序列中，造成錯誤的頻譜劃分結(jié)果。由于改進EWT是在功率譜中檢測邊界，信號的功率譜更為平滑，幾乎沒有毛刺，波峰波谷也很清晰(如圖3(a)所示)，因而可成功檢測到邊界，對頻譜進行有效劃分。

與EMD相比，EWT允許預(yù)先設(shè)定模態(tài)數(shù)量，因此信號分解結(jié)果更加穩(wěn)定[33]。如圖5所示為EMD分解結(jié)果，雖然也能將噪聲區(qū)分開來，但EMD由于終止條件不合理而出現(xiàn)過分解現(xiàn)象，即產(chǎn)生虛假分量，這些虛假分量并沒完備的理論可以解釋，多次迭代使得分解時間過長。EWT可將屬于同一部分的兩個模態(tài)分解出來，但是分解出的兩個模態(tài)具有更實質(zhì)的原信號特征，可看做兩個獨立模態(tài)。相比傳統(tǒng)EWT和EMD，改進EWT能對信號進行更加有效地分解，避免了模態(tài)混疊現(xiàn)象，算法簡單、計算速度快，分解結(jié)果更加穩(wěn)定。

圖5 EMD分解結(jié)果Fig.5 EMD decomposition result

5 工程實例分析

5.1 微震信號波形分析

為進一步驗證本文所提算法的有效性，采用千秋煤礦采集到的微震事件數(shù)據(jù)進行試驗分析,該礦安裝了波蘭ARAMIS微震監(jiān)測系統(tǒng)，傳感器采樣頻率設(shè)置為500 Hz?，F(xiàn)場監(jiān)測到的信號主要有兩類：有效的煤巖破裂微震事件；放炮、卸壓炮、松動炮等人工爆破事件。

圖6所示為系統(tǒng)采集到的兩個典型巖體微震信號和爆破振動信號，通過對比可發(fā)現(xiàn)兩種波形存差異：① 微震是礦巖變形、裂紋擴展過程中的伴生現(xiàn)象，其呈現(xiàn)的波形較為單一、衰減緩慢、尾波發(fā)育；② 由于礦山采用微差爆破，而導(dǎo)致爆破信號出現(xiàn)重復(fù)波形，其信號波形的衰減速度較快、尾波不發(fā)育；③ 通常情況下，爆破信號的峰值振幅比微震信號大。雖然通過人工可以辨別這些差異，但是很難通過上述差異對信號進行自動識別。因此，必須通過特征提取對兩種信號的差異性進行量化，進而達到有效自動識別的目的。

(a) 巖體微震信號

5.2 EWT信號分解

將典型巖體微震信號與爆破振動信號進行分解，圖7是經(jīng)過改進EWT分解后得到的結(jié)果，分別計算出各IMF分量的最大幅值。微震信號各IMF分量的最大幅值分別為0.047 8,0.129 8,0.066 5,0.040 8,0.022 1,0.001 6 mm/s；爆破振動信號各IMF分量的最大振幅為0.097 0,0.198 3,0.185 2,0.107 5,0.033 8,0.012 2 mm/s，可見巖體微震與爆破振動信號的主要振動強度集中于分量IMF1～IMF5。

(a) 巖體微震信號

獲得更為有效的信號特征將有助于后期識別準確率的提高，為此需要對分解到的IMF分量進行篩選。隨機選取巖體微震與爆破振動信號各5個進行改進EWT分解，利用式(4)計算出各IMF分量與對應(yīng)原始信號的相關(guān)系數(shù)，其結(jié)果如圖8所示。由誤差棒圖可以看出，微震與爆破信號相關(guān)系數(shù)的變化趨勢大致相同，樣本各IMF分量相關(guān)系數(shù)值誤差上下限均不大，且均在IMF6達到最小。徐峰等[34]經(jīng)多次試驗得出結(jié)論，認為相關(guān)系數(shù)小于0.03的IMF分量應(yīng)該剔除。因此，本文選取IMF1～IMF5作為最優(yōu)模態(tài)分量。

圖8 各分量與對應(yīng)原始信號的相關(guān)系數(shù)Fig.8 Correlation coefficients between each component and the corresponding original signal

5.3 特征提取

利用篩選所得最優(yōu)IMF分量進行信號重構(gòu)，并計算重構(gòu)信號的多尺度排列熵值。多尺度排列熵的計算需設(shè)置嵌入維數(shù)m、尺度因子s和延遲時間τ。其中，嵌入維數(shù)m取值范圍通常為3～7，m取值太小，重構(gòu)序列中包含的狀態(tài)太少，算法失去意義和有效性；m取值太大，計算量越大，導(dǎo)致算法的運行時間變長，并無法反映時間序列的細微變化。鄭近德等經(jīng)過試驗證明延遲時間τ對時間序列的熵值計算結(jié)果影響較小，圖9所示為微震與爆破信號在嵌入維數(shù)不變的情況下信號的排列熵值隨延遲時間τ的變化情況。從圖中可以看出延遲時間τ對信號熵值的影響不大，為此本文取τ=1。

圖9 微震與爆破信號在不同延遲時間下的排列熵Fig.9 The permutation entropy of microseismic and blasting signal at different delay time

多尺度排列熵用于衡量時間序列在不同尺度下的復(fù)雜性和隨機性，可有效檢測振動信號的動態(tài)變化，反映振動信號在不同尺度下的特征。圖10所示為不同嵌入維數(shù)下巖體微震與爆破振動信號的MPE值隨尺度因子s的變化情況。由誤差棒圖可以看出，隨著尺度因子變化，樣本的熵值也在變化，且誤差上下限都不大，因為信號的隨機性發(fā)生了變化，從而導(dǎo)致排列熵值發(fā)生變化；而由于復(fù)雜度不一樣，兩種信號的MPE值也存在明顯差異。由圖10(c)可以看出，m=7時，s=12時信號的MPE值存在重疊部分，此時很難對兩種信號進行區(qū)分；而m=3和m=5時，微震信號的MPE值明顯要小于爆破信號的MPE值，且嵌入維數(shù)m=5時，兩種信號的MPE值差異性最大。因此，本文取嵌入維數(shù)m=5。

尺度因子s的取值一般大于10，如果s取值過小，則只能在較少尺度范圍內(nèi)對時間序列進行分析，不能最大限度地提取信號的特征信息；s取值過大，有可能造成信號之間的復(fù)雜度差異被抹除。由圖10(a)和(b)可以看出，隨著尺度因子s增大，微震信號的MPE值呈遞增趨勢且逐漸接近爆破信號的MPE值。雖然當s=12時微震與爆破信號的MPE值還存在明顯差異，但是當尺度因子s過大時，兩種信號的MPE值有可能會很接近甚至有交叉，這就會影響對微震與爆破信號的自動識別。為了盡可能大的提取信號特征信息，后續(xù)研究中選取尺度因子s=12，嵌入維數(shù)m=5，延遲時間τ=1。

(a) m=3

由上述分析可以發(fā)現(xiàn)，MPE可作為識別微震信號的一個重要特征參數(shù)。然而，上述分析只是對提取的特征做了定性分析，并不能對微震信號進行定量識別。鑒于此，本文采用GK模糊聚類算法對巖體微震與爆破振動信號作進一步分類識別，以期得到較好的分類結(jié)果。

5.4 GK分類識別

從采集到的微震數(shù)據(jù)中隨機抽取200個巖體微震信號和200個爆破振動信號進行EWT分解并重構(gòu)，將重構(gòu)后的兩種信號各構(gòu)成40組樣本。每種信號取前20組數(shù)據(jù)作為已知樣本，后20組數(shù)據(jù)作為待識別樣本。選取參數(shù)優(yōu)化后的MPE值作為特征向量，通過GK模糊聚類算法確定標準聚類中心，再結(jié)合擇近原則利用式(12)計算待識別樣本與標準聚類中心的歐式貼近度，從而實現(xiàn)對兩種信號的分類識別，模式識別流程圖如圖11所示。

圖11 模式識別流程圖Fig.11 Pattern recognition flow chart

(12)

式中：n為樣本個數(shù);μA(ui)、μB(ui)分別為樣本μi對待識別模糊子集A和標準聚類中心模糊子集B的隸屬度。

圖12為巖體微震與爆破振動信號在不同模型下聚類后的二維空間分布以及二維等高線。圖12(a)為改進EWT_MPE_GK模糊聚類結(jié)果，由于同種信號的不同樣本間具有相似性的特點，兩種信號的MPE特征參量按不同信號類別分布在兩個聚類中心的周圍，且不同信號之間區(qū)分較為明顯，這說明基于改進EWT_MPE_GK模糊聚類算法對微震事件的分類識別具有很好的效果。為進一步驗證所提方法的優(yōu)越性，分別采用傳統(tǒng)EWT_MPE_GK模糊聚類和EMD_MPE_GK模糊聚類進行分析，結(jié)果分別如圖12(b)和(c)所示。由圖12可以看出，傳統(tǒng)EWT_MPE_GK和EMD_MPE_GK雖然也能將兩種信號區(qū)分開來，但是聚類效果并不理想。改進EWT_MPE_GK較傳統(tǒng)EWT_MPE_GK和EMD_MPE_GK兩種模型聚類結(jié)果，巖體微震與爆破振動信號交疊區(qū)域數(shù)據(jù)點更少，而且隸屬度等高線區(qū)分微震與爆破信號的誤判數(shù)更少，類內(nèi)MPE特征更緊密，分類效果更理想。

(a) 改進EWT_MPE_GK模糊聚類

表1為各模型對巖體微震與爆破振動信號的分類識別結(jié)果，由表1可知，改進EWT_MPE_GK模糊聚類的分類識別準確率達到93.5%，且改進EWT分解重構(gòu)信號的平均模糊熵(E)更小、分類系數(shù)(C)更接近1，其分類識別結(jié)果具有更高的準確性。對比EMD_MPE_GK模糊聚類和改進EWT_MPE_GK模糊聚類的分類識別結(jié)果，改進EWT的分類系數(shù)、平均模糊熵及準確率均優(yōu)于EMD。因此，基于改進EWT_MPE_GK模糊聚類算法能夠?qū)r體微震信號與爆破振動信號進行有效識別，與其他兩種模型相比較，C值更接近1，E值更小，進一步證明了所提方法對巖體微震信號與爆破振動信號分類識別的優(yōu)越性。

表1 基于MPE_GK模糊聚類分類識別結(jié)果Tab.1 Classification recognition results based on based on MPE_GK fuzzy cluster

6 結(jié) 論

(1) 本文針對EWT在以往處理復(fù)雜信號頻譜出現(xiàn)的過切分問題提出了新的改進方法，試驗結(jié)果表明，與傳統(tǒng)EWT和EMD相比較，改進EWT能夠更加準確有效地將信號分解；相較于EMD，改進EWT避免了端點效應(yīng)和模態(tài)混疊現(xiàn)象，算法簡單、計算速度快，分解結(jié)果更加穩(wěn)定。

(2) 多尺度排列熵用于衡量時間序列在不同尺度下的復(fù)雜性和隨機性，可有效檢測振動信號的動態(tài)變化，反映振動信號在不同尺度下的特征。本文將篩選所得最優(yōu)IMF分量進行重構(gòu)，并計算重構(gòu)信號的MPE值，結(jié)果顯示，巖體微震信號的MPE值要小于爆破振動信號的MPE值，且當嵌入維數(shù)m=5，尺度因子s=12，延遲時間τ=1時，兩種信號的MPE值差異性最大。

(3) 基于改進EWT_MPE_GK模糊聚類算法的分類識別準確率達到93.5%，其平均模糊熵更接近0，分類系數(shù)更接近1。與傳統(tǒng)EWT_MPE_GK模糊聚類和EMD_MPE_GK模糊聚類相比，本文所提方法的聚類效果更優(yōu)、識別準確率更高。因此，改進EWT_MPE_GK模糊聚類模型為礦山微震信號特征提取及分類識別提供了一種新的分析方法。