程鐵棟,王運來,張志釗,易其文,尹寶勇,袁海平

(1.江西理工大學 電氣工程與自動化學院，江西 贛州 341000；2.江西理工大學 理學院，江西 贛州 341000；3.合肥工業(yè)大學 土木與水利工程學院，安徽 合肥 230009)

巖石聲發(fā)射(acoustic emission，AE)是指巖石在受到外部載荷或內力作用下，其內部含有的大量微孔洞、裂隙發(fā)生了聚合、擴展及相互貫穿，儲存的能量以彈性波形式釋放的過程[1]。聲發(fā)射技術作為一種動態(tài)檢測方法，被廣泛應用于巖石力學工程領域的災害監(jiān)測、預警方面[2]。通過對聲發(fā)射傳感器采集的巖石聲發(fā)射信號進行分析，可以獲取巖石破壞的時間、地點、機制等信息。由于背景噪聲、機械噪聲等因素的干擾，傳感器所采集的聲發(fā)射信號波形極為復雜，導致聲發(fā)射信號的識別異常困難。因此，巖石聲發(fā)射信號的準確識別對巖體失穩(wěn)災變的預警研究具有重要意義。

隨著聲發(fā)射技術在各個領域的廣泛應用，對信號去噪能力的要求逐漸提高。劉東瀛等[3]提出一種基于相關系數(shù)原理的經(jīng)驗模態(tài)分解(empirical mode decomposition，EMD)去噪方法；戴聰聰?shù)萚4]提出一種將總體經(jīng)驗模態(tài)分解法(ensemble empirical mode decomposition，EEMD)與云相似度理論相結合的去噪方法；姜長泓等[5]通過對小波閾值法的改進，提出基于平移不變小波的聲發(fā)射信號去噪方法。上述去噪方法雖然取得了不錯的效果，但也存在不足之處。例如，EMD在分解信號的過程中易出現(xiàn)模態(tài)混疊及端點效應等問題[6]，EEMD在分解信號的過程中會存在殘余的白噪聲[7]，小波閾值去噪在分解層數(shù)的確定和閾值的選擇等方面存在困難[8]。

特征提取是聲發(fā)射信號識別的關鍵，常用的聲發(fā)射信號特征提取方法主要包括基于參數(shù)的分析方法和基于波形的分析方法[9]。參數(shù)分析法通過對信號的統(tǒng)計特征參數(shù)進行分析，是目前聲發(fā)射信號領域應用最廣泛的方法。姚旭龍等[10]利用能量貢獻率構建了巖石破裂關鍵事件聲發(fā)射信號的優(yōu)選方法；司莉等[11]通過對聲發(fā)射信號的遞歸定量分析，提出基于遞歸定量分析與支持向量機相結合的聲發(fā)射信號識別方法。這些方法為聲發(fā)射信號的識別提供了多個途徑并取得不錯的效果，但也存在一些不足：能量貢獻率作為單一特征，提供的聲發(fā)射源信息有限；遞歸定量分析存在計算量大的問題。

變分模態(tài)分解(variational mode decomposition，VMD)是Dragomiretskiy等[12]提出的一種新的非遞歸式模態(tài)分解算法，該方法可有效避免EEMD在分解信號過程中產(chǎn)生殘余的白噪聲。然而，VMD對信號進行分解時，需要提前設置本征模態(tài)分量(intrinsic mode function，IMF)的個數(shù)K，且K的取值影響算法的分解效果[13]。目前時頻域特征參數(shù)已被廣泛應用于故障診斷領域，但多個特征融合參數(shù)在聲發(fā)射領域應用較少。

鑒于此，本研究利用改進的VMD算法對巖石聲發(fā)射信號進行分解，結合排列熵篩選出最優(yōu)IMF分量并進行重構，然后提取重構信號的多個特征參數(shù)，將多個特征參數(shù)融合并輸入到多元宇宙優(yōu)化算法(multi-verse optimizer，MVO)優(yōu)化的支持向量機(support vector machine，SVM)，對巖石聲發(fā)射信號進行識別。

1 改進的VMD算法

1.1 VMD算法

VMD算法的實質是通過構建并求解約束變分模型，得到多個具有特定稀疏性的IMF分量。VMD算法構建的約束變分模型為：

(1)

式中：uk為VMD分解得到的各個IMF分量；ωk為各個IMF分量對應的中心頻率；δ(t)為單位脈沖函數(shù)；j為虛數(shù)單位；f(t)為原始信號。

為了求解該約束變分模型，通過引入二次懲罰因子和拉格朗日乘子項將有約束變分問題轉換為如下無約束變分問題：

(2)

1.2 VMD算法的改進

VMD算法需要提前設置IMF分量個數(shù)，且通過主觀經(jīng)驗難以合理設置該參數(shù)。本研究利用奇異值最佳有效秩階次自動獲得VMD分解的最優(yōu)K值。排列熵是一種檢測動力學突變和時間序列隨機性的方法[14]，其大小反映時間序列的隨機程度，可以根據(jù)排列熵值評定IMF分量的含噪程度[15]。改進VMD去噪算法的主要步驟如下。

1) 信號的奇異值分解

假設信號序列為X={x1,x2,…,xN}，將X構造為p×q階Hankel矩陣：

(3)

對Hankel矩陣H進行奇異值分解，得到：

(4)

式中：U∈Rp×p稱為左奇異矩陣，V∈Rq×q稱為右奇異矩陣，且均為正交矩陣；Σ=diag(σ1,σ2,…,σr),為對角矩陣。對角矩陣Σ中的元素σi(i=1,2,…,r)稱作H的奇異值，且奇異值為降序排列，即σ1≥σ2≥…≥σr≥0，r=min(p，q)為矩陣的H秩。

2) 求相鄰兩個奇異值之間的差值

若k階奇異值與k-1階奇異值之差最大，且與k+1階奇異值之差最小，即max(σk-σk-1), min(σk-σk+1)時，則將此時的k值確定為奇異值最佳有效秩階次。

3) 信號的VMD分解

將k設定為VMD分解的最優(yōu)K值并對信號進行分解，得到K個IMF分量。

4) 計算K個IMF分量的排列熵值，對一個長度為N的IMF分量進行相空間重構，得到矩陣

(5)

式中：m為嵌入維數(shù)，t為延遲時間，M=N-(m-1)t。

將矩陣Y中的每行按升序重新排列，得到向量中各元素位置的列索引構成的一組符號序列:

S(l)={j1,j2,…,jm}，l=1,2,…,d。

(6)

式中:m維相空間映射不同的符號序列總共有m!種，d≤m!。

每一種符號序列出現(xiàn)的次數(shù)t與m!種符號序列的總次數(shù)M的比值作為該符號序列出現(xiàn)的概率，即：

(7)

進而得到排列熵的計算式為：

(8)

將排列熵值進行歸一化處理，即：

(9)

對每一個IMF分量都進行排列熵的計算，可以得到K個IMF分量的排列熵值。

5) 信號重構

(10)

2 多特征融合提取方法

2.1 梅爾頻率倒譜系數(shù)

梅爾頻率倒譜系數(shù)(Mel frequency cepstral coefficients，MFCC)是一種廣泛應用于音頻場景和語音識別的特征參數(shù)，其計算簡單、區(qū)分能力較為突出[16]。MFCC提取信號特征的過程如下。

1) 語音信號一般需要經(jīng)過預加重、分幀以及加窗處理，變成單幀的短時信號，而本研究的聲發(fā)射信號在進行特征提取之前已經(jīng)是單幀的信號，可將每一幀信號視為一組信號進行處理。

2) 對聲發(fā)射信號作快速傅里葉變換，得到頻域數(shù)據(jù)

(11)

3) 對頻域數(shù)據(jù)取模后再平方，得到信號的功率譜

E(h)=|X(h)|2。

(12)

4) 將信號的功率譜通過Q個Mel濾波器組進行計算，得到：

(13)

式中，Hi(h)為第i個濾波器的傳遞函數(shù)，其計算式為：

(14)

式中，f(i)為中心頻率。

5) 將S(i)取對數(shù)并進行離散余弦變換，可以得到：

(15)

式中，a為Mel頻率倒譜系數(shù)的個數(shù)，即特征向量的維數(shù)。a的取值范圍通常為12～16，增大該值會使得計算量相應增加[17]，因而本研究中MFCC作為特征向量時的維數(shù)確定為12維。

2.2 短時能量

短時能量(short-time energy，SE)是對聲音信號在時域上的分析。由于聲發(fā)射信號是一種非線性、非平穩(wěn)信號，能量會隨時間變化，因此引入短時能量作為聲發(fā)射信號在時域方面的特征。短時能量的計算式為：

(16)

式中：xn(m)為輸入的第n幀信號；N為一幀信號的長度；SEn為第n幀信號的能量，即短時能量。

2.3 譜質心

譜質心(spectral centroid，SC)是描述音色屬性的重要物理參數(shù)之一，是頻率成分的重心[18]。對信號在一定頻率范圍內的能量進行加權平均，便可得到該信號的譜質心，計算式為：

(17)

式中：fi為信號的頻率，E(fi)為信號的能量譜，SC為信號做離散傅里葉變換的譜質心。

2.4 特征融合

將上述3個特征向量以列向量的形式進行組合，最終得到一個融合的特征向量

F=[F1,F2, …,F13,F14]。

(18)

式中：F1,F2, …,F12為MFCC，表征信號的頻域特征；F13為SE，表征信號的時域特征，與MFCC在時頻域起互補作用；F14為SC，反映信號的頻率分布和能量分布情況。

3 建立MVO-SVM模型

SVM是Vapnik[19]在統(tǒng)計學習理論和結構風險最小化的基礎上提出的一種通用學習方法，該方法適合解決小樣本和非線性問題。MVO是Mirjalili等[20]提出的一種基于物理機理的啟發(fā)式優(yōu)化算法，具有待調節(jié)的參數(shù)少、搜索效率高及優(yōu)化能力強等優(yōu)點，已經(jīng)被成功應用于許多工程問題中。SVM在應用過程中，需要提前設置合適的懲罰參數(shù)γ。本研究利用MVO對SVM的參數(shù)進行尋優(yōu)，可以避免因人工選取參數(shù)不合適導致模型分類精度降低問題。利用訓練集樣本作為SVM的輸入數(shù)據(jù)，并引入MVO算法優(yōu)化SVM，從而建立MVO-SVM模型。如圖1所示，MVO-SVM模型的建立過程包括5個步驟。

圖1 MVO-SVM模型的建立過程

2) 定義SVM對訓練數(shù)據(jù)的分類誤差作為MVO優(yōu)化SVM的適應度函數(shù)，當適應度函數(shù)達到最小值時，對應的宇宙分量作為當前最優(yōu)分量。

3) 依據(jù)式(20)更新宇宙。

(20)

4) 當算法運行達到最大迭代次數(shù)時，停止更新宇宙，此時輸出的全局最優(yōu)宇宙即為懲罰參數(shù)C和核函數(shù)參數(shù)γ的最優(yōu)值。

5) 以懲罰參數(shù)C和核函數(shù)參數(shù)γ的最優(yōu)值重構SVM模型。

4 實驗分析

4.1 實驗方案

實驗室模擬巖石聲發(fā)射信號，聲發(fā)射信號采集設備的工作原理如圖2所示。

圖2 聲發(fā)射信號采集設備的工作原理

聲發(fā)射信號的采集使用江西理工大學自主研發(fā)的PIE-3-L型傳感器如圖3所示，傳感器的采樣頻率為145 kHz，聲發(fā)射信號的采樣點數(shù)為4 000。實驗時，在傳感器與巖石的接觸面之間涂抹少量黃油，然后將傳感器放置于巖石的固定位置。

圖3 聲發(fā)射信號采集裝置

在距離聲發(fā)射傳感器20 cm處進行4種不同信號的采集測試：①爆破信號。將鞭炮置于巖石底下，用鞭炮爆炸產(chǎn)生的信號模擬礦山現(xiàn)場巖石爆破產(chǎn)生的信號；②斷鉛信號。按國際無損檢測界規(guī)定，選擇0.5 mmHB鉛筆芯進行斷鉛實驗，以模擬巖石破斷聲發(fā)射信號；③電鉆信號。利用電鉆對巖石進行鉆孔，模擬礦山現(xiàn)場的工況聲源；④敲擊信號。用金屬棒保持一定的力度敲擊巖石，模擬礦山現(xiàn)場的另一工況聲源。實驗分別采集不同類型的信號各105組，總計采集到數(shù)據(jù)420組。

4.2 改進的VMD去噪

從采集的聲發(fā)射信號中取一幀長度為4 000點的信號，通過計算信號的奇異值最佳有效秩階次，得到VMD分解預設的分量個數(shù)K=3，再用VMD算法對信號進行分解得到IMF1、IMF2、IMF3，進而計算3個IMF分量的排列熵值分別為0.442 1、0.504 1、0.845 5，最后確定IMF1和IMF2為最優(yōu)分量并重構，重構后的信號如圖4所示。由圖4可見，改進的VMD去噪后聲發(fā)射信號波形較為光滑，在剔除大部分噪聲的同時保留了聲發(fā)射信號的主要特征。對采集的多組聲發(fā)射信號分別運用該方法去噪，均能達到較好的效果。

圖4 改進的VMD去噪結果

原始信號分別經(jīng)EEMD、小波閾值和改進的VWD去噪處理后的頻譜如圖5所示。圖5(a)為原始聲發(fā)射信號的頻譜，信號包含了部分高頻噪聲，頻譜分布范圍寬廣。圖5(b)為EEMD去噪后的信號頻譜，可以看出EEMD去噪效果一般，無法抑制大部分噪聲。圖5(c)為小波閾值去噪后的信號頻譜圖，表明小波閾值去噪能抑制高頻噪聲，對低頻噪聲也能起到較好的抑制作用，但該方法需要人為設定閾值，在信號重構時會丟失部分有用信息。圖5(d)為改進的VMD去噪頻譜圖，可以看出，該方法對低頻和高頻噪聲均能達到良好的抑制作用，且去噪后的頻譜集中分布于某一頻率，與EEMD去噪和小波閾值去噪算法相比，改進的VMD去噪算法效果更好。

圖5 原始聲發(fā)射信號及各方法去噪后的頻譜

為定量評價改進的VMD去噪算法的性能，用信噪比(signal-to-noise ratio，SNR)對3種算法的去噪效果進行比較，結果見表1。

表1 3種去噪算法參數(shù)對比

由表1可見，相比于EEMD、小波閾值去噪，改進的VMD去噪算法的SNR值分別提高了5.23、0.84 dB，說明改進的VMD降噪方法相對較好。

4.3 融合特征向量提取

利用融合的特征向量對采集的爆破、斷鉛、電鉆和敲擊信號進行特征提取，再對所有信號樣本的特征值進行歸一化處理。為展示不同類別信號樣本的差異，表2列出了每一類信號中4個樣本的14維特征值。其中，F(xiàn)1～F12為MFCC的12維特征參數(shù)，F(xiàn)13為SE的一維特征參數(shù)，F(xiàn)14為SC的一維特征參數(shù)。由表2可見，不同類別信號的特征參數(shù)存在較大差異。

表2 4類信號的特征參數(shù)對比

4.4 聲發(fā)射信號識別

1) 信號預處理。利用改進的VMD方法對采集的4類信號共計420組數(shù)據(jù)進行去噪處理。

2) 構建融合特征向量，將多個特征參數(shù)融合作為分類識別的特征向量。劃分訓練數(shù)據(jù)和測試數(shù)據(jù)，從已有的4類信號共420組數(shù)據(jù)中，按照訓練數(shù)據(jù)與測試數(shù)據(jù)為4∶1的劃分方式對信號進行隨機分層采樣，以保證訓練數(shù)據(jù)和測試數(shù)據(jù)中的類別數(shù)目均衡。設定爆破信號的標識類別為1，斷鉛信號的標識類別為2，電鉆信號的標識類別為3，敲擊信號的標識類別為4。

3) 將訓練集和測試集分別作歸一化處理。目的是將所有的特征值大小限定在[0，1]，防止某個特征值過大或者過小。

4) 建立MVO-SVM模型。把MVO算法尋優(yōu)得到的模型參數(shù)作為懲罰參數(shù)、核函數(shù)參數(shù)。采用RBF核函數(shù)訓練SVM模型，利用訓練好的SVM模型對測試集分類，結果如圖6所示。

圖6 MVO-SVM分類識別結果

圖6中，測試集樣本數(shù)量為84個，類別標簽為1、2、3、4，樣本實際標簽與樣本預測標簽重合表示該樣本被正確識別。由圖6可見，類別1樣本中有1個樣本被誤識別為類別3、類別2樣本中有1個樣本被誤識別為類別4，類別3樣本中有1個樣本被誤識別為類別2，類別4樣本中有1個樣本被誤識別為類別1。綜上，MVO-SVM對爆破、斷鉛、電鉆及敲擊信號的識別準確率為95.24%。

為了測試去噪方法的有效性與融合特征向量的優(yōu)越性，采用改進的VMD去噪、EEMD去噪、小波閾值去噪等方法，對比不同特征提取方法在MVO-SVM模型上的識別準確率，如表3所示。

表3 不同去噪方法下各種特征提取方法的識別準確率

由表3可知：①在不同的去噪方法下，采用融合特征向量提取聲發(fā)射信號的特征，可以獲得較高的識別準確率，而單一的特征向量則無法獲得較好的效果，且僅SE和SC作為特征向量時，聲發(fā)射信號的識別準確率均較低；②在不同的特征提取方法下，利用改進的VMD去噪后識別效果整體上優(yōu)于EEMD和小波閾值去噪方法；③在EEMD和小波閾值去噪方法上，單一特征向量的識別準確率相同，而融合特征向量的識別準確率顯著提高，且在兩種去噪方法上的識別效果存在差異。因此，基于改進的VMD與融合特征向量構建的聲發(fā)射信號識別方法準確率高，對4種不同類別信號均能進行有效識別。

5 結論

1) 提出一種改進的VMD算法，利用奇異值最佳有效秩階次確定VMD分解的IMF分量個數(shù)K，并利用排列熵篩選出最優(yōu)IMF分量并重構。將改進的VMD算法應用于聲發(fā)射信號去噪，提高了信號的信噪比。

2) 利用MFCC、SE和SC等3種特征構建融合特征向量F，以此計算信號的特征值，使各類信號樣本的特征值差異更加明顯，為MVO-SVM模型識別聲發(fā)射信號提供了新的方法。

3) 構建融合特征向量的特征提取方法，在EEMD去噪和小波閾值去噪的基礎上均可以取得較好的識別效果，而改進VMD去噪的融合特征向量的識別準確率最高，可以達到95.24%。表明基于改進的VMD去噪與融合特征向量構建的巖石聲發(fā)射信號識別方法是可行的，能有效識別出聲發(fā)射信號。

本研究提出的基于改進VMD去噪與多特征融合的聲發(fā)射信號識別方法對室內模擬巖石聲發(fā)射信號可以達到較高的識別準確率，但對于實際礦山巖石聲發(fā)射信號的識別還需進一步驗證，后續(xù)將針對現(xiàn)場巖石聲發(fā)射信號進行研究，建立泛化性能更好的分類模型，以期得到較好的識別效果，為實際工程應用提供參考。