基于MFCC樣本熵和灰狼算法優(yōu)化支持向量機的天然地震與人工爆破自動識別

龐 聰,江 勇,廖成旺,吳 濤,丁 煒

(1.中國地震局地震研究所,湖北 武漢 430071;2.地震預(yù)警湖北省重點實驗室,湖北 武漢 430071;3.湖北省地震局,湖北 武漢 430071)

0 引言

天然地震與振動幅值較大的人工地震動事件具有許多相似的特性,區(qū)分研究它們的事件對地震目錄生成、地震預(yù)警算法改進、強震動儀器抗干擾等至關(guān)重要。支持向量機方法被廣泛應(yīng)用于天然地震與人工爆破、塌陷等非天然地震動事件的分類識別研究中,取得了良好的研究成果與行業(yè)應(yīng)用。如陳潤航等[1]利用SVM識別186個首都圈天然地震事件和174個人工爆破事件;范曉易等[2]以MATLAB Lib-SVM工具箱為基礎(chǔ),進行天然地震、爆破及塌陷事件的多分類研究;黃漢明等[3]從天然地震和人工爆破事件信號中提取小波能量熵,應(yīng)用支持向量機識別上述事件。支持向量機[1-3]在處理非線性問題上具有天然的優(yōu)勢,其分類機理在于選取一個核函數(shù)(線性核函數(shù)Linear、多項式核函數(shù)Poly、徑向基核函數(shù)RBF以及神經(jīng)元非線性作用核函數(shù)Sigmoid等)將低維樣本數(shù)據(jù)映射到高維空間,構(gòu)建并尋求一個最優(yōu)分類超平面,但是其分類性能受核函數(shù)參數(shù)影響較大,采取算法優(yōu)化的方式尋求一個合適的核函數(shù)參數(shù)值是當下SVM算法改進的主要方向之一。

對于事件性質(zhì)辨識而言,特征判據(jù)與分類模型皆是影響識別效果的關(guān)鍵因素,尋找良好分類性能的判據(jù)也是眾多研究人員不斷嘗試突破的工作方向之一。梅爾頻率倒譜系數(shù)是一種應(yīng)用于語音識別領(lǐng)域的數(shù)據(jù)處理方法,在近些年已被國內(nèi)外相關(guān)人員應(yīng)用到地震數(shù)據(jù)處理[4-5]中,更在2018年被陳潤航等[6]應(yīng)用至地震事件識別中,取得了良好分類效果,但是其直接利用梅爾頻率倒譜系數(shù)(Mel Frequency Cepstral Coefficent，MFCC)各系數(shù)值組成的40維系數(shù)向量作為輸入樣本,訓練過程較為復(fù)雜,存在較多的冗余或多重共線性特征參數(shù),直接影響學習器的分類效果,限制了該特征判據(jù)的推廣與應(yīng)用。

本文基于灰狼算法優(yōu)化SVM中的徑向基核函數(shù)參數(shù)值,并從樣本數(shù)據(jù)中提取梅爾頻率倒譜系數(shù)、一階差分系數(shù)及二階差分系數(shù)的樣本熵作為機器學習樣本數(shù)據(jù),并利用2013年四川蘆山7.0級地震和人工爆破事件波形記錄驗證GWO-SVM分類器和MFCC樣本熵特征的辨識效果。

1 基于灰狼算法改進的支持向量機

1.1 灰狼算法

灰狼算法[7](Grey Wolf Optimization Algorithm,GWO)是一種模擬狼群狩獵活動的新型群體仿生優(yōu)化算法,在2014年由Mirjalili首次提出,它模仿狼群的四層等級制度,由頭狼α、次級狼β、三級狼δ及普通狼ω等構(gòu)成,狼群狩獵活動步驟主要分為追蹤目標、圍捕獵物及捕捉等,狩獵過程中的指揮命令優(yōu)先級為:α>β>δ>ω。GWO算法具體步驟為:

首先狼群發(fā)現(xiàn)并確定目標位置,對目標進行包圍,二者的距離表示為:

D=|C×Xp(t)-X(t)|

(1)

第t+1次迭代后灰狼的位置為:

X(t+1)=Xp(t)-A×D

(2)

式中:X(t+1)為第t+1次迭代后灰狼的最新位置坐標;A、C皆為系數(shù)調(diào)節(jié)因子;Xp(t)為獵物位置向量。

各個等級的灰狼實施圍捕計劃時,狼會帶領(lǐng)狼對獵物實施追捕,追捕動作直接影響狼和獵物的位置變化。該過程的位置更新策略為:

(3)

(4)

Xp(t+1)=(X1+X2+X3)/3

(5)

式中:Dα、Dβ、Dδ分別表示α、β、δ狼與狼群中其他狼的距離;Xp(t+1)表示最接近獵物的灰狼位置,即最優(yōu)解。

追捕到何時階段,狼群開始攻擊、捕捉獵物(最優(yōu)解)。攻擊速度的快慢通過調(diào)節(jié)參數(shù)a來完成,且A的取值范圍為[-a,a],滿足迭代停止條件后即得到最優(yōu)解。

1.2 GWO優(yōu)化的SVM識別算法

在處理非線性問題時,基于徑向基核函數(shù)的支持向量機(SVM)分類效果很大程度上取決于參數(shù)懲罰系數(shù)c和核函數(shù)半徑σ的賦值,通過引入灰狼算法優(yōu)化上述兩個參數(shù),得到最優(yōu)解,從而改進SVM分類機(優(yōu)化流程如圖1所示)。

圖1 GWO-SVM辨識流程Fig.1 Identification process of GWO-SVM

(1) 通過某特征提取方法提取出若干個特征向量,建立特征數(shù)據(jù)集,同時按照一定比例劃分訓練集和測試集;

(2) 初始化狼群數(shù)目N、最大迭代次數(shù)t、目標獵物位置自變量數(shù)等參數(shù);

(3) 遍歷灰狼種群,計算每個個體的適應(yīng)度fi:

fi=m/n

(6)

式中:m為被SVM分類器準確識別的樣本數(shù)量;n為參與SVM分類的總樣本量。將適應(yīng)度值最大的三個灰狼位置分別記為Xα、Xβ、Xδ;

(4) 計算α、β、δ狼與狼群中其他狼的距離Dα、Dβ、Dδ,同時更新上述狼的位置,以及刷新參數(shù)a、A、C值;

(5) 循環(huán)計算(3)～(4),直至迭代次數(shù)達到迭代次數(shù)臨界值t,并得到最佳目標位置,即最優(yōu)解;

(6) 設(shè)置SVM參數(shù),將得到的最優(yōu)懲罰系數(shù)和最優(yōu)核函數(shù)半徑最優(yōu)值替換RBF核函數(shù)中默認的初始參數(shù)值,以構(gòu)建新的GWO-SVM分類器;

(7) 應(yīng)用測試集與GWO-SVM模型進行地震屬性辨識,得到天然地震與人工爆破的識別結(jié)果。

2 數(shù)據(jù)集與特征提取

2.1 數(shù)據(jù)集

選取2013年蘆山7.0級地震事件的360條強震動記錄和39條人工爆破記錄,其中蘆山地震震中位于30.3°N、103.0°E,震源深度13 km;爆破事件位于中國水利水電科學研究院巖土所試驗場,信號采樣率為1 000 Hz,采用1.58倍TNT當量的RDX球型炸藥。原始記錄統(tǒng)一截取長度為12 000的數(shù)據(jù),并進行數(shù)據(jù)歸一化、NaN去除等數(shù)據(jù)預(yù)處理操作,消除臺站所處環(huán)境、儀器精度、波形幅值等因素可能給數(shù)據(jù)分析帶來的不利影響。圖2(a)為天然地震事件震動波形圖;圖2(b)為不同人工爆破事件的固定分量信號。

圖2 天然地震與人工爆破信號Fig.2 Signals of natural earthquakes and artificial blasting

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

梅爾頻率倒譜系數(shù)[8](MFCC)是語音識別領(lǐng)域中被廣泛應(yīng)用的特征提取算法。步驟為:

(1) 預(yù)加重:采用數(shù)字濾波器補償初始信號中受抑制的高頻信號,該濾波器的傳遞函數(shù)為：

H(z)=1-μz-1

(7)

(2) 加窗:窗函數(shù)一般采用漢明窗,使得相鄰兩窗口之間有重疊區(qū)域,即

0≤n≤N-1

(8)

(3) 時頻域轉(zhuǎn)換:對各窗口信號進行離散傅里葉變換(DFT),將信號從時域變換到頻域,其計算公式為:

(9)

(4) 計算Mel濾波器組對數(shù)能量:將每個窗口的頻譜取模的平方得到功率譜,并除以該窗口采樣長度,然后經(jīng)Mel濾波器組濾波,得到第m個梅爾濾波器組輸出能量的對數(shù),即得到對數(shù)能量s(m)。

(10)

式中:Hm(k)為第m個三角濾波系數(shù)。

(5) 得到MFCC各階系數(shù):對數(shù)能量s(m)經(jīng)離散余弦變換(DCT)后,即可得到第n階梅爾頻率倒譜系數(shù)

(11)

經(jīng)過上式計算,只能得到MFCC靜態(tài)L階系數(shù),包含原始信號的靜態(tài)特性,而研究信號的動態(tài)特性就需要單獨提取動態(tài)系數(shù)(一階差分系數(shù)和二階差分系數(shù)等),即

(12)

式中:dt為第t個一階差分;Ct即第t個倒譜系數(shù);Q為倒譜系數(shù)Ct的階數(shù);K是一階導數(shù)時差值。

以天然地震和爆破事件中隨機各取出一條記錄為例,得到的MFCC靜態(tài)系數(shù)、MFCC一階差分系數(shù)、MFCC二階差分系數(shù),其中Mel濾波器的階數(shù)設(shè)置為24,離散余弦變換系數(shù)維度為12,信號長度為8 000,預(yù)加重濾波器參數(shù)取為0.937 5,窗函數(shù)采用漢明窗。由于MFCC一階差分系數(shù)及其二階差分系數(shù)的首尾兩幀皆為0,這里作刪除處理。同時,由于MFCC系數(shù)的幅值與維數(shù)一般呈現(xiàn)顯著的負相關(guān)性,維數(shù)越低,系數(shù)的幅值越大,為降低地震事件類型辨識所需的特征向量維數(shù)與特征矩陣復(fù)雜度,提高辨識效率,可只提取MFCC靜態(tài)系數(shù)矩陣、MFCC一階差分系數(shù)矩陣及MFCC二階差分系數(shù)矩陣的第一維系數(shù),具體如圖3所示。

圖3 天然地震與爆破信號的MFCC系數(shù)首維特征提取結(jié)果Fig.3 The first-dimensional extraction results of MFCC from signals of natural earthquakes and blasting

2.3 提取MFCC各階系數(shù)樣本熵

通過MFCC分析法,將得到的MFCC靜態(tài)系數(shù)向量組、MFCC一階差分系數(shù)向量組及MFCC二階差分系數(shù)向量組的首維系數(shù)向量分別記為{M1(i),i=1,2,…,N}、{M′1(i),i=1,2,…,N}及{M″1(i),i=1,2,…,N},其中N為MFCC每一維系數(shù)向量的長度。

樣本熵[9](Sample Entropy,簡稱SampEn)是一種可表征一維離散樣本內(nèi)部混亂狀態(tài)的特征值,由學者Richman等于2000年提出。由于天然地震事件與人工爆破事件的震源機制、傳播波形、瞬時最大能量、能量衰減規(guī)律等皆有明顯差異,將樣本熵應(yīng)用至MFCC各系數(shù)向量的狀態(tài)特征描述中是一種新的嘗試。樣本熵計算步驟為:

(13)

(14)

(15)

(3) 將維數(shù)增加至m+1,重復(fù)步驟(1)～(2),得到Bm+1(r)

(16)

(4) 樣本熵的計算公式為

(17)

在實際應(yīng)用中,嵌入維數(shù)m常取2,閾值r為原始樣本的標準差×自定義權(quán)值,其中權(quán)值常取0.15或0.2。

3 實驗與分析

本實驗硬件條件為:Intel(R) Core(TM) i5-8400 @ 2.80GHz,內(nèi)存為8 G的Windows10 64位操作系統(tǒng),實驗平臺為MATLAB 2019a。采用國家地震科學數(shù)據(jù)共享中心(data.earthquake.cn)提供的2013年蘆山7.0級地震事件的360條強震動三分量記錄(EW、NS、UD),及中國水利水電科學研究院巖土工程研究所(www.geoeng.iwhr.com)與陳祖煜院士團隊提供的39條人工爆破事件數(shù)據(jù),共399條加速度幅值記錄,并根據(jù)本文特征提取方法從長度一致的歸一化三分量數(shù)據(jù)中得到一個399×3的MFCC系數(shù)樣本熵特征矩陣(圖4),其特征從分別為MFCC靜態(tài)系數(shù)首維向量樣本熵、MFCC一階差分系數(shù)首維向量以及樣本熵MFCC二階差分系數(shù)首維向量樣本熵。

圖4 MFCC系數(shù)樣本熵特征提取結(jié)果Fig.4 Feature extraction results for MFCC coefficient sample entropy

由圖4,人工爆破信號(樣本號:361～399)與天然地震信號(樣本號:1～360)存在較大的區(qū)分度,爆破信號波形特征較為一致,樣本熵變化幅度較小,天然地震信號成分較復(fù)雜,不同信號波形特征差異較大,其MFCC系數(shù)樣本熵值變化明顯。MFCC靜態(tài)系數(shù)首維向量樣本熵(記作Mfcc0_SE)、MFCC一階差分系數(shù)首維向量樣本熵(記作Mfcc1_SE)以及MFCC二階差分系數(shù)首維向量樣本熵(記作Mfcc2_SE)的t檢驗結(jié)果如表1所列,顯著性水平設(shè)置為0.05,3個特征參數(shù)的假設(shè)檢驗結(jié)果均為1,符合特征參數(shù)有明顯地震事件區(qū)分能力的預(yù)期目標。

表1 天然地震與爆破信號的單一特征向量t檢驗結(jié)果Table 1 The t-test results of single eigenvector of signals of natural earthquakes and blasting

實驗方案設(shè)計為:特征矩陣按360∶39比例隨機分成訓練集Tr360×3與測試集Te39×3,對應(yīng)的標簽值向量為訓練集標簽列向量Lr360×1與測試集標簽列向量Le39×1,共含2個辨識子實驗,其循環(huán)次數(shù)分別為100次和1 000次,實驗結(jié)果如圖5和表2。

圖5 GWO-SVM與SVM模型辨識結(jié)果對比Fig.5 Comparison between identification results of GWO-SVM and SVM models

表2 1 000次循環(huán)辨識實驗結(jié)果Table 2 Identification results of experiment under 1 000 cycles

由圖5和表2可看出,GWO-SVM分類機的辨識性能各方面都優(yōu)于SVM分類機,不論是小型實驗和較大型實驗,GWO-SVM算法都保持了較為穩(wěn)定的辨識效果:以1 000次循環(huán)辨識實驗為例,GWO-SVM的識別準確率穩(wěn)健性指標—標準差、極差、變異系數(shù)分別為1.180 2、7.692 3、0.011 9,均遠低于SVM的4.382 5、23.076 9、0.048 6;在100次辨識實驗條件下,GWO-SVM的辨識曲線與SVM曲線邊界區(qū)分明顯,大部分的連續(xù)點位連線體現(xiàn)為水平分布趨勢。這證明:GWO-SVM算法具有更強的非線性目標求解能力與穩(wěn)健性,實現(xiàn)了灰狼算法優(yōu)化支持向量機關(guān)鍵參數(shù)懲罰系數(shù)c與核函數(shù)半徑σ的預(yù)定辨識目標。優(yōu)化參數(shù)列于表3。

表3 GWO-SVM辨識結(jié)果中的c與σ(部分)Table 3 c and σ in GWO-SVM identification results

表3數(shù)據(jù)來自1 000次循環(huán)辨識實驗結(jié)果中隨機抽取的12份懲罰系數(shù)與核函數(shù)半徑優(yōu)化結(jié)果,由于每一個辨識過程所需的樣本數(shù)據(jù)都是隨機組合的,這導致灰狼算法優(yōu)化得到的懲罰系數(shù)最優(yōu)值與核函數(shù)半徑最優(yōu)值也相差較大,但這絲毫不影響GWO-SVM分類機的辨識性能。優(yōu)化后的SVM分類器辨識性能得到明顯的提升,識別結(jié)果基本穩(wěn)定在95%以上,證明支持向量機中懲罰系數(shù)與核函數(shù)半徑的最優(yōu)值選取與訓練集數(shù)據(jù)息息相關(guān),也驗證了本文辨識方案的有效性和可行性。

為了驗證GWO-SVM在地震事件類型識別中的優(yōu)越性,在同樣的訓練集與測試集基礎(chǔ)上,本實驗采用RobustBoost集成學習、LDA、PLDA等3種機器學習分類模型與GWO-SVM模型進行辨識效果對比實驗,4個模型的辨識次數(shù)與識別效果如圖6所示。

圖6 4個辨識模型的識別率圖Fig.6 Identification rates of four models

由圖6可知,GWO-SVM模型除了明顯改善SVM分類機識別性能,還在同類型的機器學習方法中有一定優(yōu)勢,具有較突出的曲線平穩(wěn)特性。

4 結(jié)論

本文引用灰狼優(yōu)化算法優(yōu)化傳統(tǒng)支持向量機,并采用MFCC三維系數(shù)樣本熵作為特征樣本,提出一種新型地震事件性質(zhì)辨識方法。核心要點在于使用灰狼算法優(yōu)化支持向量機徑向基核函數(shù)(RBF)中的懲罰系數(shù)和核函數(shù)半徑,使之最大程度地匹配訓練集和支持向量機模型,形成新的GWO-SVM分類器,然后對測試集進行辨識實驗。辨識結(jié)果表明:

(1) GWO-SVM模型性能明顯優(yōu)于SVM分類器,也優(yōu)于RobustBoost、LDA、PLDA等學習機,灰狼算法的參數(shù)優(yōu)化效果高度匹配預(yù)期目標;

(2) MFCC靜態(tài)系數(shù)、MFCC一階差分系數(shù)、MFCC二階差分系數(shù)樣本熵作為地震事件類型辨識判據(jù),積極拓展了天然地震與人工爆破辨識研究在機器學習判據(jù)上的空白。

本文研究可作為一種探索性方法應(yīng)用到中小型地震預(yù)警系統(tǒng)或地震學研究中,以進行天然地震與其他非天然地震事件的分類識別,從而降低地震事件誤判概率。但是在實際監(jiān)測中,天然地震事件性質(zhì)識別除了會受到人工爆破事件的干擾外,也存在監(jiān)測儀器受到電磁干擾后生成異常波形、塌陷地震等各種干擾事件類型。并且,學習樣本數(shù)量和數(shù)據(jù)種類如果較少,只能說明在當前實驗數(shù)據(jù)下是積極可行的,缺乏普適性實驗論證。鑒于以上只是新方法的初步探討,下一步將重點收集各類相關(guān)數(shù)據(jù),研判該方法在不同震級、不同震中區(qū)域、多種地震動事件下的辨識差異性與有效性,挖掘更多有現(xiàn)實意義的地震事件類型辨識方法。

致謝:感謝中國地震局工程力學研究所和中國水利水電科學研究院巖土工程研究所提供數(shù)據(jù)支撐。