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

(1.中國(guó)地震局地震研究所,湖北 武漢 430071;2.地震預(yù)警湖北省重點(diǎn)實(shí)驗(yàn)室,湖北 武漢 430071;3.湖北省地震局,湖北 武漢 430071)

0 引言

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

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

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

1 基于灰狼算法改進(jìn)的支持向量機(jī)

1.1 灰狼算法

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

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

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

(1)

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

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

(2)

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

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

(3)

(4)

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

(5)

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

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

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

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

圖1 GWO-SVM辨識(shí)流程Fig.1 Identification process of GWO-SVM

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

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

(3) 遍歷灰狼種群,計(jì)算每個(gè)個(gè)體的適應(yīng)度f(wàn)i:

fi=m/n

(6)

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

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

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

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

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

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

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

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

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

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

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

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

H(z)=1-μz-1

(7)

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

0≤n≤N-1

(8)

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

(9)

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

(10)

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

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

(11)

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

(12)

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

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

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

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

通過(guò)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ù)向量的長(zhǎng)度。

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

(13)

(14)

(15)

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

(16)

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

(17)

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

3 實(shí)驗(yàn)與分析

本實(shí)驗(yàn)硬件條件為:Intel(R) Core(TM) i5-8400 @ 2.80GHz,內(nèi)存為8 G的Windows10 64位操作系統(tǒng),實(shí)驗(yàn)平臺(tái)為MATLAB 2019a。采用國(guó)家地震科學(xué)數(shù)據(jù)共享中心(data.earthquake.cn)提供的2013年蘆山7.0級(jí)地震事件的360條強(qiáng)震動(dòng)三分量記錄(EW、NS、UD),及中國(guó)水利水電科學(xué)研究院巖土工程研究所(www.geoeng.iwhr.com)與陳祖煜院士團(tuán)隊(duì)提供的39條人工爆破事件數(shù)據(jù),共399條加速度幅值記錄,并根據(jù)本文特征提取方法從長(zhǎng)度一致的歸一化三分量數(shù)據(jù)中得到一個(gè)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,人工爆破信號(hào)(樣本號(hào):361～399)與天然地震信號(hào)(樣本號(hào):1～360)存在較大的區(qū)分度,爆破信號(hào)波形特征較為一致,樣本熵變化幅度較小,天然地震信號(hào)成分較復(fù)雜,不同信號(hào)波形特征差異較大,其MFCC系數(shù)樣本熵值變化明顯。MFCC靜態(tài)系數(shù)首維向量樣本熵(記作Mfcc0_SE)、MFCC一階差分系數(shù)首維向量樣本熵(記作Mfcc1_SE)以及MFCC二階差分系數(shù)首維向量樣本熵(記作Mfcc2_SE)的t檢驗(yàn)結(jié)果如表1所列,顯著性水平設(shè)置為0.05,3個(gè)特征參數(shù)的假設(shè)檢驗(yàn)結(jié)果均為1,符合特征參數(shù)有明顯地震事件區(qū)分能力的預(yù)期目標(biāo)。

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

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

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

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

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

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

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

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

圖6 4個(gè)辨識(shí)模型的識(shí)別率圖Fig.6 Identification rates of four models

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

4 結(jié)論

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

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

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

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

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