吳義文，程鐵棟，易其文，趙 奎

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

通過微震監(jiān)測系統(tǒng)監(jiān)測、分析巖體損傷破裂過程中產(chǎn)生的微震信號，對微震事件進(jìn)行定位，進(jìn)而可判斷開挖過程中的巖體狀態(tài)和巖體的力學(xué)行為，并估測巖體的穩(wěn)定性[1]。然而礦山現(xiàn)場環(huán)境復(fù)雜，震源種類繁多，導(dǎo)致有效的微震事件中夾雜著許多干擾信號，其中爆破震動信號占比較大。信號在傳播時(shí)，由于信號會在不同介質(zhì)體內(nèi)衰減并同時(shí)受外界因素干擾等，使得爆破震動信號與微震信號在時(shí)域波形上較為相似，而當(dāng)前的微震監(jiān)測系統(tǒng)自動識別微震信號的效果較差。若主要依靠人工方法對微震信號進(jìn)行識別并處理，則經(jīng)常會出現(xiàn)漏處理、誤處理的情況。因此，對微震信號和爆破震動信號的自動識別研究具有重要的現(xiàn)實(shí)意義[2-3]。

目前，使用較多的微震信號特征信息提取方法有兩種，即多參數(shù)分析法和時(shí)頻分析法。時(shí)頻分析法主要包括小波變換(wavelet transform，簡稱WT)[4]、傅里葉變換(fourier transform，簡稱FT)[5]、小波包變換(wavelet packet decomposition，簡稱WPD)[6-7]、經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,簡稱EMD)[8]等方法。然而有些方法本身存在一定缺陷，會帶來不準(zhǔn)確的分析結(jié)果，如：傅里葉變換法不能很好地刻畫信號的局部特征，只能做初步的識別；小波變換和小波包變換法會受到小波基和分解層次的限制。多參數(shù)分析法對微震信號的識別效果較好[9-10]，但其對采集信號的準(zhǔn)確度有較高的要求，由于特征參數(shù)復(fù)雜，實(shí)現(xiàn)過程相對困難。

經(jīng)驗(yàn)小波變換(empirical wavelet transform,簡稱EWT)在頻譜上進(jìn)行適當(dāng)?shù)姆指睿⒃诟鱾€(gè)頻帶上建立小波濾波器組，從而將信號分解成多個(gè)頻率特征信息相異的分量[11]。與EMD算法相比較，EWT算法具有良好的理論基礎(chǔ)，并解決了EMD等算法出現(xiàn)的模態(tài)混疊問題，且計(jì)算量要遠(yuǎn)小于EMD算法。目前，EWT算法已經(jīng)被廣泛應(yīng)用于各個(gè)領(lǐng)域，如醫(yī)學(xué)領(lǐng)域[12]、故障診斷[13]等，但基于EWT算法的微震信號識別研究還鮮見報(bào)道。

基于此，嘗試將EWT算法引入到微震信號分析領(lǐng)域。奇異值分解(singular value decomposition，簡稱SVD)能有效壓縮矩陣數(shù)據(jù)，并精細(xì)描述信號序列本征屬性，將其與時(shí)頻分析算法進(jìn)行組合能夠得到不同的特征提取模型，如WPD-SVD[7]、EMD-SVD[8]等。EMD-SVD模型將分解出的模態(tài)分量作為矩陣的行向量，然后按頻率高低依次排列構(gòu)造復(fù)合矩陣，對矩陣進(jìn)行SVD計(jì)算，算出的奇異值是從大到小排列的，但各個(gè)分量所對應(yīng)的奇異值與分量的頻率并不是一一對應(yīng)的，容易使各個(gè)分量的特征信息丟失[14]。鑒于此，利用EWT算法分解得到多個(gè)分量，再利用互信息對分量進(jìn)行篩選，進(jìn)而分別對保留的分量構(gòu)建Hankel矩陣，利用SVD計(jì)算Hankel，提取各個(gè)Hankel矩陣的奇異值平均值、方均根值、標(biāo)準(zhǔn)差作為模型識別的特征量，最后通過構(gòu)建的支持向量機(jī)(support vector machine, 簡稱SVM)自動識別網(wǎng)絡(luò)對微震和爆破信號進(jìn)行準(zhǔn)確識別，以期得到一種有效的微震信號識別方法。

1 經(jīng)驗(yàn)小波變換

2013年，法國學(xué)者Gilles提出了經(jīng)驗(yàn)小波變換，將模態(tài)看作調(diào)幅-調(diào)頻(AM-FM)信號。針對AM-FM窄頻帶、高幅值的特征，EWT算法可實(shí)現(xiàn)顯著模態(tài)的搜索及提取。其主要流程如下：

1)由于在截取記錄時(shí)，樣本長度選擇不合適或者某些外界原因(如傳感器的基礎(chǔ)運(yùn)動、零點(diǎn)漂移等造成的波形偏移)會導(dǎo)致趨勢項(xiàng)的產(chǎn)生，因此通過最小二乘法對原始信號進(jìn)行去趨勢項(xiàng)處理[14];

2)通過傅里葉變換法計(jì)算輸入信號的頻譜，并歸一化到[0,π]，設(shè)置模態(tài)個(gè)數(shù)M；

3)計(jì)算2個(gè)連續(xù)局部極大值之間的中間頻率點(diǎn)，作為劃分頻譜的分界ωn(n=1,2，…，M-1)；

4)依據(jù)分割的頻譜構(gòu)建小波函數(shù)ψn(ω)和尺度函數(shù)φn(ω)[11]；

5)應(yīng)用傅里葉反變換計(jì)算F(ω)×ψn(ω)和F(ω)×φn(ω)，從而可得到各個(gè)分量的時(shí)域表達(dá)式。

2 奇異值分解

原始信號經(jīng)EWT變換后，將得到的各個(gè)分量(f1，f2，…，fJ)依尺度按行排列，構(gòu)成矩陣M，然后利用SVD對矩陣M進(jìn)行分解，將算出的奇異值作為識別的特征參數(shù)，這是最自然、最直接、應(yīng)用最廣泛的時(shí)頻分析法與SVD組合的特征提取模型。但相關(guān)研究表明，除非已從分量信號中提取出了差異明顯的特征信息，否則對矩陣M直接進(jìn)行奇異值分解并不能對特征信息的提取效果起到較大的改進(jìn)作用[15]。為了能更有效地利用SVD對微震信號進(jìn)行特征提取，需要先對得到的各個(gè)分量分別進(jìn)行處理。對每個(gè)分量分別構(gòu)建Hankel矩陣Ha，可得下式：

(1)

式中1≤a≤J，且1

特征矩陣Ha行列維數(shù)的選擇應(yīng)滿足[16]：當(dāng)N為偶數(shù)時(shí)，行數(shù)m=N/2+1；當(dāng)N為奇數(shù)時(shí)，行數(shù)和列數(shù)均為m=n(N+1)/2。

對特征矩陣H進(jìn)行奇異值分解，可得下式：

(2)

3 仿真信號研究

為了驗(yàn)證EWT有效分解信號的能力，構(gòu)建仿真信號x(t)進(jìn)行考查。仿真信號x(t)由調(diào)頻信號、調(diào)幅信號和正弦信號組成，t∈[0,1]，并將分解結(jié)果與EMD法進(jìn)行比較，仿真信號如下：

EWT的輸出是由N-1個(gè)小波函數(shù)和1個(gè)尺度函數(shù)分別進(jìn)行濾波的結(jié)果。EWT的頻譜分割圖及其小波濾波器組圖如圖1所示，仿真信號EMD、EWT的分解結(jié)果如圖2、圖3所示。

(a)頻譜分割圖 (b)小波濾波器組圖

圖2 仿真信號EMD分解結(jié)果

圖3 仿真信號EWT分解結(jié)果

由圖1可見，仿真信號EWT在頻譜中計(jì)算出全部的極大值，將所有極大值按降序排列，根據(jù)前期假設(shè)的模態(tài)個(gè)數(shù)確定保留的極大值，以2個(gè)相鄰極大值之間的中間頻率點(diǎn)作為劃分頻譜的界限，再通過尺度函數(shù)和小波函數(shù)對仿真信號進(jìn)行濾波。

由圖2可以看出，imf1、imf2、imf3出現(xiàn)了模態(tài)混疊情況，另外得到了多個(gè)低頻分量，這些分量原本屬于同一主分量，但由于EMD算法采用不合理的終止條件，出現(xiàn)過分解，產(chǎn)生了虛假模態(tài)，這樣既消耗計(jì)算時(shí)間又影響算法性能，而且虛假模態(tài)的產(chǎn)生會影響后期的特征提取。

由圖3可知，通過EWT算法得到的分量f1～f3分別與原信號中50、25、16 Hz的分量相對應(yīng)，且各個(gè)分量的吻合度和原始信號相比非常接近，不存在虛假分量。由此可見，EWT算法在提取微震信號特征時(shí)，具有明顯優(yōu)勢。

4 工程實(shí)例研究

4.1 EWT分解

從漂塘鎢礦微震監(jiān)測系統(tǒng)數(shù)據(jù)中選取典型的爆破震動信號與微震信號進(jìn)行對比分析。該礦安裝的是南非IMS微震監(jiān)測系統(tǒng)，傳感器的采樣頻率為 6 000 Hz。通過對采集到的微震信號和爆破震動信號波形進(jìn)行對比分析可知：微震信號持續(xù)時(shí)間相對較長，衰減較緩慢；爆破震動信號的持續(xù)時(shí)間較短，衰減較快，并且爆破震動信號的幅值要大于微震信號的幅值。這些分析結(jié)果通常是人工識別微震信號的一個(gè)依據(jù)，但地下地質(zhì)條件復(fù)雜，接收到的信號不一定都存在這種差異，因此需要通過特征提取的方式對 2種信號進(jìn)行量化。礦區(qū)監(jiān)測到的巖體微震和爆破震動信號波形如圖4所示。

(a)微震信號波形 (b)爆破震動信號波形

對微震信號和爆破震動信號進(jìn)行EWT分解，微震信號得到7個(gè)分量，爆破震動信號得到8個(gè)分量，如圖5所示。

(a)巖體微震

為了能更有效地提取信號特征，需要利用互信息量(MI，用符號IMI表示)對得到的分量進(jìn)行篩選。IMI描述2個(gè)變量間存在共同信息量的多少，IMI越大則變量間的共同信息越多，相關(guān)性越強(qiáng)[18]。對于信號的第i個(gè)分量ci(t)與原始信號x(t)，二者之間的IMI定義為：

(3)

式中：p(ci)和p(x)分別是ci(t)與x(t)的邊緣概率分布；p(ci，x)是ci(t)與x(t)的聯(lián)合概率分布。

依據(jù)公式(3)計(jì)算各個(gè)分量與原始信號的互信息量。劉長良等[18]認(rèn)為，當(dāng)分量與原始信號的互信息量小于0.02時(shí)，可將該分量作為虛假分量進(jìn)行剔除，所以選取圖5(a)中前7個(gè)分量作為特征分量。各分量與原始信號的互信息量如圖6所示。

圖6 各分量與原始信號的互信息量

4.2 特征提取

利用篩選的主分量f1～f7分別構(gòu)建Hankel矩陣Ha(a=1,2,…,7)，并對各個(gè)Hankel矩陣進(jìn)行奇異值分解(SVD)，從而得到7組不同奇異值σ1,σ2,…,σr，每組奇異值σ1,σ2,…,σr由矩陣Ha唯一確定。由SVD性質(zhì)可知，矩陣Ha這些非零奇異值σ1,σ2,…,σr反映的是矩陣的特征，由這些非零矩陣組成的特征向量σ=(σ1,σ2,…,σr)唯一表征了微震信號的特征。為突顯微震信號和爆破震動信號的差異，對特征向量σ進(jìn)行分析，計(jì)算平均值h1、方均根值h2及標(biāo)準(zhǔn)差h3，并作為微震信號的特征參數(shù)。h1、h2、h3表達(dá)式如下：

(4)

(5)

(6)

根據(jù)微震信號和爆破震動信號各個(gè)分量的特征參數(shù)h1、h2、h3的數(shù)值變化，探討特征參數(shù)h1、h2、h3是否能有效地區(qū)分微震與爆破震動信號。由式(1)構(gòu)建Hankel矩陣，并對Hankel矩陣進(jìn)行奇異值分解，然后通過式(4)～(6)得到各個(gè)分量的奇異值平均值h1、方均根值h2、標(biāo)準(zhǔn)差h3。各分量Hankel矩陣的平均值、方均根值、標(biāo)準(zhǔn)差變化如圖7所示。

(a)平均值

由圖7(a)可看出，爆破震動信號前6個(gè)分量的奇異值平均值均大于微震信號前6個(gè)分量的奇異值平均值，其中前3個(gè)分量的奇異值平均值差異較大；而第6個(gè)分量的奇異值平均值近似相等，爆破震動信號第7個(gè)分量的奇異值平均值要小于微震信號的奇異值平均值。

由圖7(b)、(c)可看出，爆破震動信號各分量Hankel矩陣的方均根值和標(biāo)準(zhǔn)差都要大于微震信號的方均根值和標(biāo)準(zhǔn)差，且前5個(gè)分量的差異性較大；爆破震動信號與微震信號第6個(gè)分量的方均根值近似相等。

由圖7可知，巖體微震信號和爆破震動信號各個(gè)分量的奇異值平均值、方均根值、標(biāo)準(zhǔn)差均存在一定的差異，且前5個(gè)分量差異較大，所以，利用各個(gè)分量的奇異值平均值、方均根值、標(biāo)準(zhǔn)差作為微震信號識別的特征量是可行的。然而，上述過程只是對特征參數(shù)的定性分析，未對微震信號做定量識別。鑒于此，筆者借助機(jī)器學(xué)習(xí)方法進(jìn)一步挖掘數(shù)據(jù)的內(nèi)在信息，利用模式識別方法對微震信號進(jìn)行識別分類。

4.3 SVM分類識別

SVM是在結(jié)構(gòu)風(fēng)險(xiǎn)最小化的基礎(chǔ)上構(gòu)建的機(jī)器學(xué)習(xí)方法，其克服了小樣本學(xué)習(xí)和局部極大值的問題，并在泛化能力和學(xué)習(xí)精度間進(jìn)行了較好的權(quán)衡。近年來，有學(xué)者將SVM應(yīng)用于地震震相識別[20]、巖爆預(yù)測[21]等方面，并取得了不少研究成果。筆者選擇SVM分類器作為微震信號辨識的分類器，利用微震信號和爆破震動信號各分量的奇異值平均值、方均根值、標(biāo)準(zhǔn)差作為特征量，建立基于SVM的微震信號自動識別體系。

1)從已有數(shù)據(jù)中隨機(jī)選取巖體微震和爆破震動信號各200組，利用上述方法得到微震信號和爆破震動信號前7個(gè)分量的奇異值平均值、方均根值、標(biāo)準(zhǔn)差，并將此作為特征量；將前100組作為訓(xùn)練樣本數(shù)據(jù)，后100組作為測試樣本數(shù)據(jù)。將巖體微震信號的類別設(shè)定為1，爆破震動信號的類別設(shè)定為2。

2)分類器參數(shù)選擇：選用RBF (radial basis function，簡稱RBF)徑向基核函數(shù)ψ(xi,xj)=e-γ‖xi-xj‖2,其中γ>0,懲罰參數(shù)C=2,核函數(shù)參數(shù)γ=1。

依據(jù)構(gòu)建的SVM微震信號自動識別體系，對微震信號和爆破震動信號進(jìn)行識別，為了測試本方法特征提取的有效性，與目前應(yīng)用最廣泛的奇異值分解法進(jìn)行比較。該方法是將EWT分解得到的7個(gè)分量構(gòu)建成復(fù)合矩陣，并直接對復(fù)合矩陣進(jìn)行SVD處理，將得到的奇異值作為特征量進(jìn)行識別，兩種特征提取方法的識別效果如表1所示。

表1 兩種特征提取方法識別效果對比

由表1可知，采用EWT_Hankel和SVD提取微震和爆破震動信號各個(gè)分量的奇異值平均值、方均根值、標(biāo)準(zhǔn)差作為特征量，可以較好地實(shí)現(xiàn)兩種信號的識別，且準(zhǔn)確率達(dá)到92.5%；采用EWT和SVD提取微震和爆破震動信號奇異值作為特征量，對微震信號的識別準(zhǔn)確率較低，僅為83.5%。所以基于EWT_Hankel_SVD和SVM構(gòu)建的微震信號自動識別模型，可以有效地識別爆破震動信號和微震信號。這為獲取礦山信號特征、研究礦山信號分類提供了一種新思路和新方法。

5 結(jié)論

1)經(jīng)驗(yàn)小波變換是一種非線性、非平穩(wěn)信號的分析處理工具，其結(jié)合了小波變換和EMD算法的優(yōu)點(diǎn)，具有充分的理論基礎(chǔ)和自適應(yīng)性。仿真信號結(jié)果表明，EWT算法能夠有效地解決模態(tài)混疊問題，從而更加準(zhǔn)確地分解信號，特別適用于處理微震信號類的非平穩(wěn)信號。

2)時(shí)頻分析法和SVD組合的模式應(yīng)用最廣泛，該模式通常是將時(shí)頻分析得到的各個(gè)主分量構(gòu)建成復(fù)合矩陣，并直接對矩陣進(jìn)行SVD處理；筆者將得到的各個(gè)主分量分別構(gòu)建Hankel矩陣，之后通過SVD計(jì)算每個(gè)Hankel矩陣的奇異值，將各個(gè)Hankel矩陣奇異值平均值、方均根值、標(biāo)準(zhǔn)差這3個(gè)值作為模式識別的特征量，顯著提高了微震信號識別率。

3)通過EWT_Hankel_SVD和SVM網(wǎng)絡(luò)構(gòu)建的微震信號自動識別模型分類識別準(zhǔn)確率達(dá)到92.5%；應(yīng)用最廣泛的奇異值分解特征提取法識別準(zhǔn)確率為83.5%。表明基于EWT_Hankel_SVD的特征提取方法是有效的，且具有較高的準(zhǔn)確率。該方法為研究礦山信號、獲取礦山信號特征提供了一種新思路和新方法。