何正祥，彭平安，廖智勤

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

0 引言

微震監(jiān)測技術(shù)是礦山生產(chǎn)活動過程中地壓安全監(jiān)測分析重要的技術(shù)環(huán)節(jié)[1-2]。微震監(jiān)測系統(tǒng)通過將其監(jiān)測區(qū)域內(nèi)的各類微震信號進行實時拾取，并分析巖體破裂活動來判斷巖體狀態(tài)。實際采礦生產(chǎn)過程中，微震信號存在著多種類型[3]，如巖體破裂、爆破振動、電磁噪聲和機械鉆鑿等。傳統(tǒng)的人工分類處理效率較低且經(jīng)常由于各種干擾導致識別錯誤。因此，如何快速準確地識別不同的信號類型是實現(xiàn)微震監(jiān)測系統(tǒng)實時化、智能化的重要環(huán)節(jié)。

在目前用于微震信號識別分類的方法中，主要應用選取波形的相關(guān)計算的指標作為數(shù)據(jù)參數(shù)構(gòu)建識別模型，如Vallejos和董隴軍等[4-5]通過分析波形和事件選取得到的若干個參數(shù)作為識別指標構(gòu)建了識別模型。這類方法具有一定的識別效果，但是參數(shù)的選取人為因素較高，屬于非穩(wěn)定參數(shù)。對微震信號進行時頻分析技術(shù)總結(jié)其特征，如曹安業(yè)、陸菜平等[6-7]利用時頻分析從信號的頻譜特性進行總結(jié)，得出巖體破裂信號有較快的衰減速度、較發(fā)育的尾波、震幅大和較窄的頻帶分布等重要特征，但是該方法僅僅為人工識別提供了有效的輔助手段，無法實現(xiàn)信號的自動識別分類。另外，小波分析對波形更加細致的處理為微震信號的識別分類提供了新的思路，如郭濤、唐禮忠、趙國彥、朱權(quán)潔等[8-11]利用小波變換對微震信號進行了更深層次的分解，用以對微震信號進行識別，但是小波分析的方法仍然存在著小波基的選取、分解后的特征參數(shù)選擇等方面的問題需要解決。

上述研究做了大量的分析探討，并取得了較多成果，但是仍然在微震高效識別分類上存在一定缺陷。本文利用梅爾倒譜系數(shù)法提取微震信號的特征參量，結(jié)合混合高斯隱馬爾可夫模型作為識別分類模型，實現(xiàn)了對4類常見的礦山復雜微震信號的初步識別。該方法具有較高的識別準確率和自動化程度，為實現(xiàn)微震信號自動識別提供了新的方法。

1 基于梅爾倒譜系數(shù)的特征提取

1.1 梅爾倒譜系數(shù)

Davis和 Mermelstein首先提出了梅爾倒譜系數(shù)，是對信號波形的短時能量譜的一種表示，是將微震信號的對數(shù)功率譜通過線性余弦變換運算投影至非線性梅爾標度中所得[12]。梅爾標度與頻率的轉(zhuǎn)換關(guān)系如式(1)所示：

(1)

式中：f為信號的頻率，Hz。

梅爾倒譜系數(shù)對信號的失真有很好的補償能力[13]。

1.2 梅爾倒譜系數(shù)提取過程

梅爾倒譜系數(shù)參數(shù)的提取方法[14-15]如下：

1)高通濾波：將微震信號波形通過高通濾波器可以有效地增強高頻部分，使信號波形的頻譜起伏減小，使得高低頻的任意的頻帶都可以基于近似的信噪比求頻譜。

2)波形分幀：每隔N個采樣點將微震信號波形分段，形成新的波形單位，稱為幀。根據(jù)特征提取的要求和信號的長度一般在256或512中選取N值，同時，使每幀的前后幀都與其有一小部分交叉，目的是為了避免連續(xù)的幀之間的變化過大。

3)加漢明窗：為了增加微震信號波形分幀后每一幀與相鄰幀的連續(xù)性，讓波形的每一幀與漢明窗相乘。假設微震信號為S(n),n=0,1,…,N-1，乘以漢明窗后則有：

(2)

式中：N為微震信號分幀后的幀數(shù)，0≤n≤N-1。

4)快速傅里葉變換：能量分布的不同，可以表示不同信號的特性，因此將微震信號轉(zhuǎn)換成頻域上的能量譜進行對比。在將微震信號進行連續(xù)重疊分幀后，對分解后的每一信號幀作快速傅里葉變換，借以計算頻率域上的能量譜。

5)三角帶通濾波器：利用一系列的三角帶通濾波器在梅爾標度上轉(zhuǎn)換快速傅里葉計算得到的能量譜，得到1組系數(shù)。該系列濾波器為一系列的三角窗，均勻重疊地排列在梅爾頻率軸上。

6)計算對數(shù)能量譜：對每個濾波器組進行計算，并對結(jié)果取對數(shù)，得到的值為其對數(shù)能量，即可得到相應頻帶的對數(shù)功率譜，計算公式如下：

(3)

式中：s(m)為對數(shù)能量；X(k)為微震信號能量譜；Hm(k)為濾波器組，其中m=1,2,…,M，M為濾波器的個數(shù)。

7)離散余弦變換：利用離散余弦變換把頻譜變換到時域上，所得結(jié)果就是標準的梅爾頻率倒譜系數(shù)。倒譜系數(shù)的數(shù)學計算公式為：

(4)

式中：n為計算的幀數(shù)，0≤n≤N；m為梅爾倒譜系數(shù)個數(shù)，0≤m≤M,M取12。

在實際用于識別分類時，為了改善識別效果，顯示梅爾倒譜系數(shù)在時間軸上的變化，將12個特征參數(shù)加上差量倒譜參數(shù)產(chǎn)生24維的特征參數(shù)向量。經(jīng)過差分運算得到的差量倒譜參數(shù)顯示了原有的12個梅爾倒譜系數(shù)在時間域上的變化。差量倒譜參數(shù)Dt(n)的計算方法如式(5)所示：

(5)

式中：C(n)為信號第n幀計算出的梅爾倒譜系數(shù)；Dt(n)為第t個一階差量倒譜參數(shù)；為式中一階導數(shù)的時間差，一般取=2, 1≤θ≤Θ 。

根據(jù)上述流程，可以從原始的復雜微震信號中成功提取出24維特征參數(shù)向量。

2 微震信號的自動識別分類

2.1 自動識別方法

實現(xiàn)信號的自動識別分類，基本方法是通過算法提取出波形的特征，利用這些特征結(jié)合機器學習實現(xiàn)自動識別分類。本文使用梅爾倒譜系數(shù)法從微震波形中提取出24維的特征向量，并且使用混合高斯隱馬爾可夫模型實現(xiàn)對微震信號的自動識別分類。

隱馬爾可夫模型(Hidden Markov Model, HMM)是關(guān)于時序的概率模型[16]。1個HMM通?？梢杂?個參數(shù)表示：λ=(N,M,π,A,B)，其中，N為HMM中馬爾可夫狀態(tài)鏈的大小，實際使用中該值為固定值。設N個狀態(tài)為θ1,θ2,…,θN，則n時刻的狀態(tài)為qn，qn∈(θ1,θ2,…,θN)。M為馬爾可夫狀態(tài)鏈中與每個狀態(tài)可能相對應的觀察值數(shù)。設M個觀察值為V1,V2,…,VM，則n時刻的觀察值為on，on∈(V1,V2,…,VM)。π為初始狀態(tài)概率分布，π∈(π1,π2,…,πN)，其中πi=P(q1=θi), 1 ≤i≤N。參數(shù)A為狀態(tài)轉(zhuǎn)移概率矩陣:A= [aij]N×N，其中aij=P(qn+k=θj|qn=θi)，1 ≤i,j≤N，表示在任意時刻n，若狀態(tài)為θi，則在下一時刻狀態(tài)為θj的概率。B為觀察值概率矩陣:B=[bij]N×M，其中bij=P(on=Vk|qn=θi)，1≤I≤N, 1≤j≤M，表示在任何時刻n，狀態(tài)為θi時觀測值Vk被獲取的概率。

而為了更好地對復雜微震信號進行識別，混合高斯隱馬爾可夫模型在原有的隱馬爾可夫模型的技術(shù)上對觀察值概率密度函數(shù)利用混合高斯模型進行建模，將bjk修改為當前狀態(tài)和觀察值之間的高斯分布概率密度函數(shù)，即：

(6)

式中：μjm為均值；Ujm為方差；cjm為高斯分布權(quán)重。構(gòu)成了混合高斯隱馬爾可夫模型。

2.2 微震信號自動識別分類流程

基于梅爾倒譜系數(shù)的微震信號自動識別分類主要分為3個步驟：

1)微震信號的特征提取。本文根據(jù)1.2節(jié)中的提取過程，將微震信號進行分幀，并計算每幀的梅爾倒譜系數(shù)，得到12個梅爾倒譜系數(shù)值，再通過對此12個梅爾倒譜系數(shù)值做差分計算。故1個微震信號可以提取出1個24維的特征向量，該特征向量是自動識別的數(shù)據(jù)基礎。

2)自動識別分類模型的訓練。選取一定時間段已標記的微震信號通過梅爾倒譜法提取特征向量，并按照類別將該特征向量作為對應事件類型訓練集數(shù)據(jù)，依次得到4個訓練集。作為混合高斯隱馬爾可夫模型的輸入數(shù)據(jù)，通過迭代求出相應的4組不同的模型參數(shù)N,M,π,A,B。每組模型參數(shù)構(gòu)建的混合高斯隱馬爾可夫模型即為對應微震事件類型的識別模型。微震信號特征提取及識別模型訓練流程如圖1所示。

圖1 微震信號特征提取及識別模型訓練流程Fig.1 Process of feature extraction and recognition model training of microseismic signal

3)基于識別模型的微震事件識別分類應用。根據(jù)訓練得到的4種微震事件(巖體破裂、爆破振動、鉆機鑿巖和電磁干擾)識別模型，對該礦山生產(chǎn)活動中產(chǎn)生的微震事件做梅爾倒譜系數(shù)計算，將得到的特征向量分別代入4種微震識別模型中計算相應的混合高斯隱馬爾可夫模型概率值，取概率值最大的模型對應的事件類型，作為待判斷微震事件的類型，以此實現(xiàn)生產(chǎn)活動中微震信號的自動識別分類。其自動識別分類流程如圖2所示。

3 微震信號自動識別分類應用

3.1 冬瓜山銅礦微震監(jiān)測系統(tǒng)

冬瓜山銅礦位于銅陵市以東處, 礦床處于獅子山銅礦區(qū)的深部,是我國首座開采深度達到千米并具有明顯巖爆傾向性的硬巖金屬礦山。冬瓜山銅礦具有礦床儲量較大、埋藏深度較高的賦存特點，設計開采模式為高產(chǎn)強化開采，該方法的特點為開采盤區(qū)和采場較多，這導致了采場的分布范圍廣泛，并且多采場并行開采的模式也導致了推進速度過快，巖爆事件過多，且大范圍分布。因此，引入微震監(jiān)測為該礦山開采過程中的巖爆活動進行預測預警，2005年 8月28日冬瓜山銅礦正式運行微震監(jiān)測系統(tǒng)，該系統(tǒng)對冬瓜山生產(chǎn)活動進行實時監(jiān)測，并在監(jiān)測數(shù)據(jù)的基礎上對微震活動做了大量的研究[17]。

3.2 基于冬瓜山銅礦微震數(shù)據(jù)的識別模型訓練

選取冬瓜山銅礦2017年6月1日-2017年12月13日之間的部分微震數(shù)據(jù)作為訓練集，選取已人工標記的巖體破裂事件、爆破振動事件、鉆機鑿巖事件和電磁干擾事件各350件，共計1 400件微震事件，用以訓練識別分類模型。

圖2 微震信號自動識別分類流程Fig.2 Automatic identification and classification of microseismic signals

根據(jù)圖1所示流程對選取的數(shù)據(jù)進行模型訓練。首先，對每個事件波形進行梅爾倒譜系數(shù)提取，圖3～6為從每種事件類型中隨機選取1個事件的波形及其梅爾倒譜系數(shù)值；然后，將手動標記的事件類型與梅爾倒譜系數(shù)值組成的24維特征向量作為輸入數(shù)據(jù)，確定混合高斯隱馬爾科夫模型參數(shù)N,M,π,A,B；最后，通過迭代算法得出模型的最優(yōu)參數(shù)，即為訓練完畢的微震信號自動識別分類模型。

圖3 某巖體破裂事件的波形及其梅爾倒譜系數(shù)值Fig.3 The waveform of a rock burst event and Mel-frequency cepstral coefficients

圖4 某爆破振動事件的波形及其梅爾倒譜系數(shù)值Fig.4 The waveform of a blast vibration event and Mel-frequency cepstral coefficients

圖5 某鉆機鑿巖事件的波形及其梅爾倒譜系數(shù)值Fig.5 The waveform of a rig drilling event and Mel-frequency cepstral coefficients

圖6 某電磁干擾事件的波形及其梅爾倒譜系數(shù)值Fig.6 The waveformof a electromagnetic interference event and Mel-frequency cepstral coefficients

3.3 自動識別分類結(jié)果與分析

利用冬瓜山銅礦微震監(jiān)測系統(tǒng)于2017年12月13日-2018年1月17日(其中2017年12月23日與2017年12月24日由于系統(tǒng)檢修，監(jiān)測系統(tǒng)沒有數(shù)據(jù))監(jiān)測到的微震事件，對建立的混合高斯隱馬爾科夫識別分類模型進行測試，所有事件均已經(jīng)人工反復確認并進行類型標記，可作為計算模型識別準確率的依據(jù)。冬瓜山銅礦這一期間每日監(jiān)測到的各類型微震事件的數(shù)量如圖7所示。

根據(jù)圖2的微震信號自動識別分類流程，對圖7中所示的事件進行類型識別分類，并與人工標記的事件類型相對比，測試本文識別模型的準確率。

對于構(gòu)建的自動識別分類模型，需要根據(jù)礦山實際情況確定相關(guān)模型參數(shù)：混合高斯方法中的高斯方程數(shù)目Q；馬爾科夫鏈狀態(tài)數(shù)N；馬爾科夫鏈中觀測值數(shù)M。其中，觀測值M就是本文的梅爾倒譜系數(shù)的維數(shù)，即M=24。對于高斯方程數(shù)目Q，通過測試不同數(shù)目下的識別準確率進行選擇。圖8為識別模型中高斯方程數(shù)對識別準確率的影響分析，如圖8所示，當Q≥3時，即可達到較高的識別準確率，但較多的高斯方程會加大模型的計算量，因此本文選用Q=4。同樣，對于狀態(tài)數(shù)N，通過測試不同的狀態(tài)數(shù)N值進行對比。圖9為識別模型中馬爾科夫鏈狀態(tài)數(shù)對識別準確率的影響分析，如圖9所示，識別準確率隨著N的增大而增大，但后續(xù)增加幅度較小，考慮到計算效率，本文實例中選用N=6。

根據(jù)上述選取的模型參數(shù)，對實際微震信號進行識別。每種事件識別的情況如圖10所示。

由圖7可知，2017年12月13日-2018年1月17日之間的微震事件總數(shù)為981個，其中：巖體破裂事件467件，爆破振動事件138件，鉆機鑿巖事件108件，電磁干擾事件268件。圖10為識別模型識別后的每日各類型微震事件數(shù)量，由圖10可知，巖體破裂事件識別錯誤數(shù)為26個，爆破振動事件識別錯誤數(shù)為32個，鉆機鑿巖事件識別錯誤數(shù)為2個，電磁干擾事件識別錯誤數(shù)為14個。因此，基于梅爾倒譜系數(shù)法的微震事件自動識別分類模型在該測試期間的識別準確率為92.46%，識別準確率較高。

圖7 冬瓜山銅礦2017年12月12-2018年1月17日每日微震事件數(shù)量Fig.7 Total daily microseismic events of Dongguashan copper mine from December 12, 2017 to January 17, 2018

圖8 識別模型中高斯方程數(shù)對識別準確率的影響分析Fig.8 Influence analysis of gaussian equation number on recognition accuracy in recognition model

圖9 識別模型中馬爾科夫鏈狀態(tài)數(shù)對識別準確率的影響分析Fig.9 Influence of the number of markov chain states on the recognition accuracy in the recognition model

從識別效果來看，基于梅爾倒譜系數(shù)與混合高斯隱馬爾科夫識別模型的微震信號識別分類方法具有較高的準確率，可以用于常規(guī)微震監(jiān)測系統(tǒng)的信號自動分類，以期提高微震信號分類的時效性和準確性，為礦山開采活動中的微震監(jiān)測系統(tǒng)真正實現(xiàn)實時監(jiān)測、實時分析提供了技術(shù)支撐。

圖10 識別模型識別后的每日各類型微震事件數(shù)量Fig.10 The number of each type of microseismic events per day after recognition of the model

4 結(jié)論

1)利用梅爾倒譜系數(shù)對波形進行特征提取，其原理是將頻率域上的信號轉(zhuǎn)化為梅爾頻率域上的非線性頻譜，并對其做離散余弦變換得到梅爾倒譜系數(shù)，結(jié)合梅爾倒譜系數(shù)差分運算，可以從波形中提取出24組特征向量，為信號的識別分類提供特征數(shù)據(jù)。

2)基于梅爾倒譜系數(shù)提取出的信號特征，訓練構(gòu)建4種微震事件對應的混合高斯隱馬爾科夫識別模型，利用冬瓜山銅礦實際生產(chǎn)活動中監(jiān)測的微震數(shù)據(jù)進行測試，識別準確率達到92.46%，具有較高的識別準確率。

3)基于梅爾倒譜系數(shù)的復雜微震信號自動識別分類方法適用于等長采樣的微震監(jiān)測數(shù)據(jù)和一定的數(shù)據(jù)基礎。同時，該自動識別分類方法實現(xiàn)了礦山微震監(jiān)測系統(tǒng)的高時效性，為實現(xiàn)微震監(jiān)測系統(tǒng)的實時性分析提供了技術(shù)支持。