基于改進(jìn)的HHT-SVM微震信號特征提取及分類識別研究*

熊璐偉,鐘曉陽,李庶林,葉龍珍,方 鑫,鄭宗檳

(1.廈門大學(xué) 建筑與土木工程學(xué)院,福建 廈門 361005;2.江西鐵山垅鎢業(yè)有限公司,江西 贛州 341000;3.自然資源部丘陵山地地質(zhì)災(zāi)害防治重點(diǎn)實驗室,福建 福州 350002)

0 引言

隨著淺層礦產(chǎn)資源的日益減少,我國礦山開采逐漸向深部發(fā)展,工程上面臨的沖擊地壓等巖石破裂災(zāi)害日益頻繁,嚴(yán)重威脅著礦山的生產(chǎn)安全。相較于傳統(tǒng)的監(jiān)測手段,微震監(jiān)測技術(shù)可監(jiān)測到巖體的微破裂活動并能動態(tài)地揭示巖體的損傷程度[1],且已在礦山、邊坡等巖體破裂災(zāi)害預(yù)警中取得不錯的成果[2-4],保障大量人員的生命財產(chǎn)安全,為各相關(guān)行業(yè)的生產(chǎn)安全提供有力支持。由于礦山工作環(huán)境復(fù)雜,爆破、電磁、機(jī)械設(shè)備、車輛等信號不時產(chǎn)生,很大程度上影響對微震信號的識別。其中電磁、機(jī)械設(shè)備、車輛信號具有較為固定的頻率、幅值和波形特征,容易剔除,而爆破致使巖石破裂產(chǎn)生的爆破信號在波形上與微震信號有著較強(qiáng)的相似性,難以準(zhǔn)確辨識。對有效信號進(jìn)行提取是微震監(jiān)測預(yù)警的關(guān)鍵一步,因此如何提高對微震信號、爆破信號的識別成功率具有重要意義。

近年來,在微震信號識別與分析領(lǐng)域,已有大量學(xué)者進(jìn)行研究。裴琳[5]對隧洞開挖時搜集到的信號進(jìn)行時頻和頻域分析,將信號分為巖石破裂信號、爆破信號及施工干擾信號3類。朱權(quán)潔等[6]對礦山爆破振動、巖石破裂信號進(jìn)行小波包分解,求取指定特征頻帶上重構(gòu)信號的分形盒維數(shù),建立特征向量,使用SVM實現(xiàn)對信號的識別。陸菜平等[7]使用快速傅里葉變換對微震信號的幅頻特性與功率譜密度進(jìn)行研究。何正祥等[8]運(yùn)用梅爾倒譜系數(shù)法實現(xiàn)對微震事件信號的準(zhǔn)確識別。Li等[9]使用多分形方法對微震與爆破信號進(jìn)行處理,指出2個信號在簡單分形維數(shù),奇異指數(shù)范圍與多分形頻譜寬度存在一定差異,可作為信號識別特征。

以上許多學(xué)者從頻率角度入手對微震信號進(jìn)行分類識別。本文對微震信號、爆破信號采用經(jīng)ICEEMDAN改進(jìn)的HHT算法進(jìn)行處理,分析2種信號的瞬時能量分布規(guī)律,尋找2種信號的差異點(diǎn),并對差異點(diǎn)進(jìn)行量化處理,提取相應(yīng)特征,構(gòu)建多維特征向量,建立經(jīng)功率譜熵特征加權(quán)與GSWOA算法改進(jìn)后的SVM支持向量機(jī)微震信號識別網(wǎng)絡(luò),研究結(jié)果可為辨識微震信號與爆破信號提供1種新的方法與思路。

1 基本理論

1.1 支持向量機(jī)

支持向量機(jī)(support vector machine,SVM)算法在解決小樣本、非線性和高維的二分類問題方面有較大優(yōu)勢,已被廣泛應(yīng)用于信號識別、文本分類等領(lǐng)域。其原理是通過使用核函數(shù)將低維特征空間轉(zhuǎn)化為高維度特征空間,尋求能夠?qū)颖具M(jìn)行分類的最優(yōu)超平面。假設(shè)訓(xùn)練樣本集為(xi,yi),i=1,2,…,l,xi∈Rn,yi∈R,SVM通過映射核函數(shù)φ(x)將xi映射到高維特征空間,超平面表達(dá)式如式(1)所示:

ω(x)=wφ(x)+b

(1)

式中:w為權(quán)重向量;φ(x)為SVM的核函數(shù);b代表原點(diǎn)與超平面間的距離。

上述求解超平面問題可轉(zhuǎn)化為求解式(1)的最小值問題,如式(2)所示:

(2)

式中:C為懲罰因子;l為樣本集數(shù)量;Lε為損失函數(shù)。

為了考慮擬合誤差,引入松弛變量,使得優(yōu)化目標(biāo)變?yōu)槿缡?3)所示:

(3)

式中:ξi為松弛變量。

為了求解式(3)的最優(yōu)解,引入拉格朗日函數(shù)將上述問題轉(zhuǎn)換成對偶問題,則可得到最終的超平面表達(dá)式,如式(4)所示:

(4)

本文選用徑向基核函數(shù)(RBF),SVM的核心是求解1個二次規(guī)劃問題,使得目標(biāo)函數(shù)f(xi)最小,如式(5)所示:

(5)

式中:f(xi)為目標(biāo)函數(shù);標(biāo)簽yi=±1。

1.2 功率譜熵特征加權(quán)算法

支持向量機(jī)在對樣本進(jìn)行分類時,將樣本的每一維特征都賦予相同的權(quán)重,但事實上,樣本的各維度特征對分類效果的貢獻(xiàn)大小不一定一致[10]。若對分類效果貢獻(xiàn)較低的特征參數(shù)賦予與對分類效果貢獻(xiàn)較高的特征參數(shù)相同的權(quán)重,則必然會影響最終的分類效果。因此本文使用信息學(xué)中熵的概念對特征向量賦予權(quán)重。在信息學(xué)中,熵能有效度量隨機(jī)變量的不確定性程度,能反映隨機(jī)變量包含的信息量。本文采用功率譜熵對樣本特征向量賦予權(quán)重,若樣本特征的功率譜熵越大,則說明這一維度特征的不確定性越高,對分類效果的貢獻(xiàn)較小,應(yīng)分配較小的權(quán)重;反之,則應(yīng)分配較大的權(quán)重。

若訓(xùn)練樣本集數(shù)量為N,特征維度為d,則特征向量矩陣為Φ=[Φ1,Φ2,…,Φd],其具體表達(dá)式如式(6)所示:

(6)

樣本集的第k維度特征可表達(dá)為Φk,對樣本的第k維度特征的功率譜熵進(jìn)行計算,計算公式如式(7)所示:

(7)

式中:H(k)為第k維度特征的功率譜熵;pl為離散序列信號的第l頻段的能量占總能量的比重大小。

H=[H(1),H(2),…,H(d)]表示d維特征的功率譜熵,第k個維度的特征權(quán)重計算公式如式(8)所示:

(8)

式中:ωq(k)為第k個維度的特征權(quán)重。

對樣本特征向量加權(quán)后,SVM的優(yōu)化問題轉(zhuǎn)化為如式(9)所示:

(9)

式中:ωq為特征權(quán)重向量。

1.3 全局搜索策略的鯨魚算法

根據(jù)鯨魚的捕食習(xí)慣,文獻(xiàn)[11]中提出WOA算法,主要包括包圍獵物、旋轉(zhuǎn)搜尋及隨機(jī)搜尋3部分。但由于此算法存在收斂速度慢且易陷入局部最優(yōu)解的問題,本文采用劉磊等[12]提出的全局搜索策略的鯨魚算法(global search whale optimization algorthm,GSWOA),其原理如下:

1)包圍獵物

鯨魚將目前搜捕獵物的信息進(jìn)行共享,然后鯨魚向離獵物最近的鯨魚靠近,逐漸縮小包圍圈,實現(xiàn)對獵物的包圍,此時鯨魚位置更新表達(dá)式如式(10)所示:

(10)

式中:t為迭代搜尋次數(shù);X為鯨魚位置;X*是全局最優(yōu)位置;A和B為系數(shù)矩陣;D為權(quán)重系數(shù)。

2)自定義權(quán)重與變螺旋搜尋

在對鯨魚位置進(jìn)行更新時,增加自適應(yīng)權(quán)重來減少處于最佳位置鯨魚對其鯨魚位置更新的影響,減少算法陷入局部最優(yōu)的風(fēng)險;同時為了使鯨魚能以多種方式搜尋獵物,引入變螺旋搜尋方法,自適應(yīng)權(quán)重及位置更新公式如式(11)～(12)所示:

(11)

(12)

式中:w(t)為自適應(yīng)權(quán)重;Xrand(t)為隨機(jī)鯨魚位置;p為[0,1]均勻分布隨機(jī)數(shù)。

3)最優(yōu)鄰域擾動

為了提高收斂速度以及避免陷入局部最優(yōu)解,引入了最優(yōu)鄰域擾動策略,在當(dāng)前最優(yōu)值附近進(jìn)行隨機(jī)搜尋,尋找更佳的全局值,其定義如式(13)所示:

(13)

按照式(10)～(13)完成種群位置更新,在預(yù)定迭代次數(shù)下求解最佳適應(yīng)度對應(yīng)的參數(shù)值。

1.4 改進(jìn)的希爾伯特-黃變換

希爾伯特-黃變換(hilbert-huang transform,HHT)是目前分析爆破、振動信號的主要方法。其主要由經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition,EMD)與Hilbert變換組成。其中EMD算法分解的本征模態(tài)分量(intrinsic mode function,IMF)存在模態(tài)混疊現(xiàn)象,導(dǎo)致分解效果不理想。改進(jìn)自適應(yīng)完備經(jīng)驗?zāi)Ｊ椒纸?improved complete ensemble empirical mode decomposition with adaptive noise,ICEEMDAN)算法通過改進(jìn)自適應(yīng)分解過程,大幅減少了信號分解后產(chǎn)生的殘余分量,獲得更加良好的重構(gòu)信號,一定程度上改善了模態(tài)混疊問題。因此本文采用ICEEMDAN算法[13]對HHT算法進(jìn)行改進(jìn)。ICEEMDAN算法的計算流程可分為以下5個步驟:

1)基于原始信號S構(gòu)造N個含可控噪聲信號,如式(14)所示:

S(i)=S+β0E1(ω(i))i=1,2,…,N

(14)

式中:S(i)表示第i個構(gòu)造信號;β0為第1次分解時信號的噪聲標(biāo)準(zhǔn)偏差;β0的取值范圍一般為[0.1,0.3],本文取0.2;ω(i)表示第i個零均值單位方差白噪聲;E1(·)表示第1個IMF算子。

2)計算S(i)的局部均值并求平均,得到第1個殘余分量,計算公式如式(15)所示:

(15)

式中:M(·)為局部均值函數(shù)。

3)將原信號S減去第1個殘余分量r1得到第1個模態(tài)(k=1),如式(16)所示:

(16)

4)計算第k個模態(tài)(k≥2),計算公式如式(17)～(18)所示:

(17)

(18)

5)返回4),計算第k+1個模態(tài),直到殘余分量不可再分或迭代次數(shù)滿足終止條件,本文將迭代次數(shù)設(shè)置為100次。

使用互信息法對分解后的模態(tài)進(jìn)行篩選,計算各IMF分量與原信號的互信息量MIi,并定義平均互信息量MIm,計算公式如式(19)所示:

(19)

式中:K為ICEEMDAN算法的分解個數(shù)。

若MIi>MIm,則認(rèn)為該分量為有效分量。

希爾伯特變換主要原理為先將信號使用ICEEMDAN算法分解成多個IMF分量,使用互信息法進(jìn)行篩選,對有效分量進(jìn)行Hilbert變換得到希爾伯特譜,其表達(dá)式如式(20)所示:

(20)

式中:Re表示取實部進(jìn)行運(yùn)算;F為分解信號頻率;T為各信號采樣點(diǎn)對應(yīng)時間;n為分解信號個數(shù);ai(T)表示瞬時幅值;φi(T)表示瞬時頻率。

將Hilbert譜的平方對頻率進(jìn)行積分,可得到相應(yīng)的Hilbert瞬時能量譜,如式(21)所示:

(21)

式中:IE(T)為希爾伯特譜瞬時能量。

2 基于改進(jìn)HHT微震信號瞬時能量特征分析

2.1 工程背景簡介

黃沙鎢礦位于贛州市于都縣,地處南嶺東西構(gòu)造帶大余——會昌隆起帶與新華復(fù)系于山隆起帶的反接復(fù)合區(qū),礦區(qū)構(gòu)造復(fù)雜,裂隙發(fā)育。由于黃沙礦區(qū)開采時間長,礦脈密集,采深大,坑場體積大,特別是礦脈密集的采空區(qū)或大采區(qū),有利于形成直通地表的大型倒楔形不穩(wěn)固結(jié)構(gòu)體,使通過采空區(qū)的斷層不穩(wěn)定性加大,導(dǎo)致地壓的增強(qiáng),致使局部區(qū)域出現(xiàn)冒頂、片幫等問題。針對日益嚴(yán)重的地壓問題,于2023年引入具有大范圍、立體監(jiān)測特點(diǎn)的實時在線多通道微震監(jiān)測系統(tǒng)。微震系統(tǒng)采用加速度型傳感器,其響應(yīng)頻率范圍為50～5 000 Hz,數(shù)據(jù)采集儀的最高采集頻率為30 kHz。從微震系統(tǒng)數(shù)據(jù)庫內(nèi),選取典型的爆破信號與微震信號進(jìn)行分析。

2.2 信號瞬時能量變化規(guī)律分析

通過改進(jìn)的HHT變換對拾取波形進(jìn)行處理分析,揭示信號瞬時能量的變化規(guī)律特征,典型微震、爆破信號的時域及瞬時能量圖如圖1～2所示。

圖1 典型微震信號時域及瞬時能量圖Fig.1 Time domain and transient energy diagram of typical microseismic signal

通過分析峰前能量的變化規(guī)律,爆破信號的瞬時能量在爆破事件開始后迅速達(dá)到能量峰值,而由于微震信號來源于巖體內(nèi)部微破裂的產(chǎn)生、發(fā)展、貫通過程中產(chǎn)生的應(yīng)力波[14],導(dǎo)致微震信號從事件開始發(fā)生至能量達(dá)到峰值的過程相對緩慢,即事件孕育過程速度較慢。分析峰后瞬時能量的變化規(guī)律,從峰值能量點(diǎn)后瞬時能量線的下降斜率可以明顯看出爆破信號在達(dá)到能量峰值后衰減的速度超過微震信號的峰后衰減速度。從圖2可以看出爆破信號的后期部分衰減波動較大,這是由于爆破信號往往會伴隨著2～3個地震子波,且爆破信號的地震子波能量較高,對爆破信號后期的幅值衰減有著較大的影響,導(dǎo)致爆破信號具有持續(xù)時間相對較長的特點(diǎn)[15]。

圖2 典型爆破信號時域及瞬時能量圖Fig.2 Time domain and transient energy diagram of typical blasting signal

2.3 信號瞬時能量變化規(guī)律研究及特征提取

當(dāng)信號瞬時能量初次達(dá)到1%峰值能量時標(biāo)志著事件開始,當(dāng)信號能量最后一次衰減至1%峰值能量時標(biāo)志著事件的終止,如圖1(b)所示,A點(diǎn)為事件的起始點(diǎn),C點(diǎn)為事件的終止點(diǎn)。通過峰值能量點(diǎn)B將整個事件(AC段)劃分為事件孕育階段(AO段)與事件衰減階段(OC段)。

通過將信號事件按瞬時能量劃分為3個階段進(jìn)行分析,發(fā)現(xiàn)微震信號、爆破信號在孕育速度、持續(xù)時間、峰后衰減速度具有一定差異。對差異點(diǎn)進(jìn)行定量分析,將事件孕育階段時間與整個事件持續(xù)時間之比記作量化結(jié)果1,對事件的孕育速度進(jìn)行描述;將整個事件持續(xù)時間記作量化結(jié)果2;將事件衰減階段按時間分為20個區(qū)域,并統(tǒng)計每個區(qū)域內(nèi)的累計能量值即瞬時能量對該區(qū)域時間進(jìn)行積分,將峰值能量后的6個區(qū)域的累計能量值進(jìn)行記錄,并將6個區(qū)域間的累計能量的衰減率記作量化結(jié)果3～7,對信號峰后衰減速度進(jìn)行描述,衰減率計算公式如式(22)所示:

(22)

式中:Ei表示第i個部分的區(qū)間累計能量。

從微震監(jiān)測系統(tǒng)中隨機(jī)抽取微震信號與爆破信號各200組,對這400組信號采用上述量化方法進(jìn)行處理后,得到各量化結(jié)果概率密度分布如圖3所示。

3 改進(jìn)的SVM信號識別網(wǎng)絡(luò)

3.1 改進(jìn)的SVM信號識別網(wǎng)絡(luò)建立

基于功率譜熵特征加權(quán)與GSWOA算法改進(jìn)的SVM信號識別網(wǎng)絡(luò)的具體實現(xiàn)步驟如下:

1)構(gòu)造特征向量,建立訓(xùn)練、測試數(shù)據(jù)集,設(shè)置信號類型標(biāo)簽。選取量化結(jié)果1～7構(gòu)建7維特征向量。從信號數(shù)據(jù)庫中隨機(jī)抽取200組微震信號樣本、200組爆破信號樣本,將60%樣本作為訓(xùn)練集,40%樣本作為測試數(shù)據(jù)集。將微震信號類型標(biāo)簽設(shè)置為1,爆破信號類型標(biāo)簽設(shè)置為2。

2)數(shù)據(jù)預(yù)處理。對訓(xùn)練、測試集數(shù)據(jù)使用[0,1]區(qū)間歸一化方法進(jìn)行預(yù)處理,計算訓(xùn)練樣本集特征向量的功率譜熵及權(quán)重,計算結(jié)果如表1所示。

表1 特征參數(shù)功率譜熵值及權(quán)重Table 1 Entropy values and weights of power spectrum of feature parameters

3)SVM網(wǎng)絡(luò)參數(shù)選取。SVM的分類器選用徑向基核函數(shù)RBF,經(jīng)GSWOA算法對RBF中懲罰系數(shù)C,核函數(shù)參數(shù)g進(jìn)行優(yōu)化。設(shè)置鯨魚數(shù)量為30,迭代次數(shù)100,懲罰系數(shù)C與核函數(shù)參數(shù)g的臨界區(qū)域分別設(shè)置為[0,100]、[0,2.00]。通過100次迭代計算后,得到優(yōu)化參數(shù)組合[54,0.21]。

3.2 分類識別結(jié)果

將經(jīng)改進(jìn)HHT提取的訓(xùn)練集特征向量及其權(quán)重與經(jīng)GSWOA算法尋優(yōu)后得到的C與g的優(yōu)化參數(shù)組合[54,0.21]輸入SVM識別網(wǎng)絡(luò)中,隨后輸入測試集數(shù)據(jù)并進(jìn)行分類識別,最終的綜合識別成功率為96.250%,微震、爆破信號具體識別情況如圖4所示。

圖4 信號識別效果圖Fig.4 Signal identification effect

為了驗證本文提出的基于功率譜熵特征加權(quán)算法與GSWOA算法改進(jìn)后的SVM信號識別網(wǎng)絡(luò)的可行性與優(yōu)越性,將其與傳統(tǒng)SVM、GSWOA-SVM、基于功率譜熵特征加權(quán)-SVM方法對比,未進(jìn)行參數(shù)優(yōu)化的SVM中C與g參數(shù)組合取[10,0.10],分類識別結(jié)果如表2所示。結(jié)果表明改進(jìn)后的HHT的識別結(jié)果整體優(yōu)于傳統(tǒng)HHT算法,說明改進(jìn)HHT算法能更加充分地提取微震、爆破信號特征。對于同一特征向量,經(jīng)功率譜熵特征加權(quán)算法及GSWOA算法改進(jìn)的SVM與改進(jìn)HHT相結(jié)合的方法對的識別成功率達(dá)到96.250%,驗證了本文建立的信號識別網(wǎng)絡(luò)的優(yōu)越性與可行性,表明本文所提方法能有效識別微震與爆破信號。

表2 分類識別結(jié)果Table 2 Classification and identification results

4 結(jié)論

1)使用ICEEMDAN改進(jìn)HHT算法對微震信號與爆破信號進(jìn)行處理,分析信號的瞬時能量變化規(guī)律,并對400組信號進(jìn)行定量研究,提取信號事件孕育速度、事件持續(xù)時間及峰后衰減速度的特征并構(gòu)建多維特征向量。

2)通過定量研究分析,構(gòu)建多維特征向量的各量化結(jié)果分類識別效果具有差異,使用功率譜熵特征加權(quán)算法對多維特征向量賦予權(quán)重,能有效提高對微震信號的識別成功率。

3)改進(jìn)的HHT分類效果優(yōu)于傳統(tǒng)HHT算法;使用改進(jìn)HHT對信號進(jìn)行特征提取并構(gòu)建特征向量,輸入至經(jīng)功率譜熵特征加權(quán)與GSWOA算法改進(jìn)的SVM信號識別網(wǎng)絡(luò)中,最終測試集的綜合識別成功率為96.250%,表明本文提出的改進(jìn)算法是有效可行的,能夠?qū)ξ⒄?、爆破信號進(jìn)行有效分類識別,為后續(xù)利用微震信號進(jìn)行監(jiān)測預(yù)警奠定基礎(chǔ)。