龐 聰 廖成旺 江 勇 程 誠(chéng) 吳 濤 舒 鵬 丁 煒

1 中國(guó)地震局地震研究所，武漢市洪山側(cè)路40號(hào)，430071 2 地震預(yù)警湖北省重點(diǎn)實(shí)驗(yàn)室，武漢市洪山側(cè)路48號(hào)，430071 3 湖北省地震局，武漢市洪山側(cè)路48號(hào)，430071 4 運(yùn)城學(xué)院數(shù)學(xué)與信息技術(shù)學(xué)院，山西省運(yùn)城市復(fù)旦西街1155號(hào)，044031 5 南華大學(xué)計(jì)算機(jī)學(xué)院，湖南省衡陽市常勝西路28號(hào)，421001

隨著人類生產(chǎn)活動(dòng)的急劇增加，距離地震臺(tái)站較近或震感強(qiáng)烈的人工爆破事件有一定概率被地震觀測(cè)儀器視作天然地震事件所記錄，產(chǎn)生的誤記錄操作將會(huì)污染已有的地震目錄甚至產(chǎn)生地震災(zāi)害應(yīng)急誤判。天然地震事件與人工爆破事件性質(zhì)辨識(shí)是當(dāng)下地震學(xué)研究及地震觀測(cè)技術(shù)應(yīng)用的重要內(nèi)容，結(jié)合機(jī)器學(xué)習(xí)或深度學(xué)習(xí)的分類決策方法在其中占據(jù)重要地位[1]。

地震事件性質(zhì)辨識(shí)的主要手段是尋找合適的地震波形特征及特征識(shí)別方法。在地震特征提取領(lǐng)域，已有學(xué)者基于時(shí)域或頻域特征深入研究，如香農(nóng)熵小波系數(shù)[2]、波形對(duì)稱度[3]、HHT提取IMF特征[4]等；在地震特征識(shí)別領(lǐng)域，常見的識(shí)別方法有決策樹[5]、支持向量機(jī)[6]、Boosting[7]、線性判別分析[8]等，但這些方法存在易過擬合、分類錯(cuò)誤率較高或?qū)Χ嘀毓簿€性特征敏感等弊端。

本文首先對(duì)原始記錄進(jìn)行歸一化處理，提取排列熵、近似熵及香農(nóng)熵等特征值形成實(shí)驗(yàn)數(shù)據(jù)集，利用基于SVM改進(jìn)的LSSVM二分類方法對(duì)天然地震、非天然地震事件進(jìn)行性質(zhì)分類，采用機(jī)器學(xué)習(xí)混淆矩陣中的準(zhǔn)確率、召回率及F-measure指標(biāo)對(duì)分類結(jié)果進(jìn)行有效性評(píng)價(jià)，結(jié)合其他機(jī)器學(xué)習(xí)方法的辨識(shí)效果對(duì)比分析LSSVM算法的優(yōu)劣性能。

1 事件性質(zhì)分類方法

最小二乘支持向量機(jī)(least square SVM，LSSVM)擅長(zhǎng)處理大規(guī)模非線性問題，數(shù)據(jù)訓(xùn)練速度較快，屬于支持向量機(jī)的改進(jìn)型方法[9]，其變化有3點(diǎn)：1)優(yōu)化目標(biāo)函數(shù)采用2NF(第2范式)；2)將不等式約束條件轉(zhuǎn)化為等式約束條件，進(jìn)而求解線性方程組，加快了求解速度，更有利于處理較大規(guī)模問題；3)對(duì)排序后的Lagrange乘子進(jìn)行裁剪處理，增強(qiáng)解的稀疏性和決策函數(shù)的適用性。

LVSVM的分類過程為：

1)對(duì)于LSSVM分類問題，約束條件是不等式約束：

yk[ωTφ(xk)+b]=1-ek,k=1,2,…,N

(1)

式中，γ為權(quán)重參數(shù)；ω為超平面法向量；ek為偏離超平面的松弛變量；φ(xk)為核函數(shù)定義的非線性映射；b為超平面的截距；yk表示線性分類一般模型{yk|yk=ω·xk+b}的分類類別，xk為訓(xùn)練樣本數(shù)據(jù)。

2)采用Lagrange乘數(shù)法轉(zhuǎn)換為求解α的極大值問題：

L(ω,b,e;α)=J(ω,e)-

(2)

式中，αk為L(zhǎng)agrange乘子。分別對(duì)ω、αk、b、ek求導(dǎo)，使其等于0：

(3)

根據(jù)式(3)列出關(guān)于α和b的線性方程組：

(4)

定義核矩陣Ωkl：

Ωkl=ykylφ(xk)Tφ(xl)=ykylK(xk,xl)，

k,l=1,2,…,N

(5)

3)解以上方程組可以得到α和b，最終得到LSSVM的分類表達(dá)式為：

(6)

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

特征提取是展開地震事件屬性識(shí)別的重要基礎(chǔ)工作，提取的特征參數(shù)直接關(guān)系到事件性質(zhì)判斷的準(zhǔn)確性和辨識(shí)模型的合理性。本文利用信息論中的熵來表征地震事件記錄的特征。熵(entropy)是一種描述系統(tǒng)內(nèi)部混亂狀態(tài)的度量值，由克勞修斯于1865年在解決熱力學(xué)問題時(shí)提出。常見的熵值類型包括香農(nóng)熵、交叉熵、條件熵、相對(duì)熵等，已在應(yīng)用物理、生命科學(xué)、應(yīng)用數(shù)學(xué)及人工智能等學(xué)科領(lǐng)域得到應(yīng)用。地震儀器記錄到的震級(jí)較大的地震記錄具有明顯的波形統(tǒng)計(jì)特征和不規(guī)則性(局部混亂)，采用熵概念中的排列熵[10]、近似熵[11]、香農(nóng)熵[12]等表征地震事件性質(zhì)是一種新的嘗試和探討。

2.1 排列熵

排列熵(permutation entropy，Hp)是2002年由Bandt等[10]提出的，其值的大小表示時(shí)間序列的復(fù)雜程度或隨機(jī)性，值越小，信號(hào)的復(fù)雜性越低；值越大，信號(hào)的不規(guī)則性越嚴(yán)重。

設(shè)信號(hào)時(shí)間序列為{X(i),i=1,2,…,n}，進(jìn)行相空間重構(gòu)后得到如下矩陣：

(7)

其中，m為嵌入維數(shù)，τ為延時(shí)因子，k=n-(m-1)τ,j=1,2,…,k。矩陣共有k個(gè)重構(gòu)分量，每個(gè)重構(gòu)分量有m維嵌入元素。

將式(7)中的第j個(gè)分量按數(shù)值大小升序排列，得到：

x(i+(j1-1)τ)≤x(i+(j2-1)τ)≤…

≤x(i+(jm-1)τ)

(8)

每個(gè)重構(gòu)分量都可以得到一個(gè)重構(gòu)符號(hào)序列：

S(l)=(j1,j2,…,jm)

(9)

其中，l=1,2,…,k，滿足k≤m!。每個(gè)重構(gòu)分量是m維空間映射到m維的符號(hào)序列，共有m!種排列方式。

計(jì)算每種m維符號(hào)序列的概率P1,P2,…,Pk，根據(jù)香濃熵的定義，信號(hào)時(shí)間序列X(i)的排列熵定義為：

(10)

當(dāng)Pj=1/m!時(shí)，排列熵達(dá)到最大值ln(m!)。

2.2 近似熵

近似熵(approximate entropy，ApEn)是1991年由Pincus[11]提出的，描述m維向量擴(kuò)展至m+1維向量時(shí)的相似性或周期性度量值。作為一個(gè)非負(fù)值，ApEn的值越大，表示時(shí)序數(shù)據(jù)具有越強(qiáng)的隨機(jī)性或非周期性，可用來描述信號(hào)的復(fù)雜程度。近似熵的計(jì)算步驟為：

1)按照式(11)構(gòu)建m維向量，用y(i)表示，即{y(i),i=1,2,…,M,M=N-m+1}，其中y(i)={x(i),x(i+1),…,x(i+m-1)}，m為嵌入維數(shù)，是重構(gòu)序列的長(zhǎng)度。計(jì)算y(i)與y(j)任意分量之間的歐氏距離d{y(i),y(j)}，并將各分量之間最大的距離定義為最大貢獻(xiàn)成分距離D{y(i),y(j)}，得到：

D{y(i),y(j)}=max{|y(i+k)-y(j+k)|}

(11)

其中，j,i=1,2,…,N-m+1；k=0,1,…,m-1。

(12)

(13)

ApEn(m,r)=φm(i)-φm+1(i)

(14)

其中，當(dāng)嵌入維數(shù)m選取2或3時(shí)，有效閾值可設(shè)定為0.15倍或0.2倍的時(shí)序標(biāo)準(zhǔn)差值。

2.3 香農(nóng)熵

香農(nóng)熵H是包含平均信息量多少的度量指標(biāo)，熵值越大，說明數(shù)據(jù)中的信息量越多[12]。其定義為：

(15)

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

本文實(shí)驗(yàn)數(shù)據(jù)來源于中國(guó)地震局工程力學(xué)研究所記錄的2021年青?，敹郙S7.4地震事件與云南漾濞MS6.4、MS5.6地震事件及人工爆破事件，其中源于地震事件的三分量波形數(shù)目為158，源自人工爆破干擾事件的波形數(shù)目為342，共500份波形數(shù)據(jù)。事件及數(shù)據(jù)描述如下：

1)根據(jù)中國(guó)地震臺(tái)網(wǎng)中心測(cè)定，2021-05-21 21:21云南省大理州漾濞縣發(fā)生MS5.6地震，震源深度10 km，共有28個(gè)臺(tái)站(表1)記錄到此次強(qiáng)震動(dòng)事件，其中漾濞臺(tái)震中距最小，為5.4 km，東西、南北及垂直向加速度峰值分別為-339.2 cm/s2、267.2 cm/s2、-220.1 cm/s2，計(jì)算儀器地震烈度為7.7度。

表1 云南漾濞地震事件臺(tái)站信息

2)2021-05-21 21:48云南省大理州漾濞縣(25.67°N，99.87°E)發(fā)生MS6.4地震，震源深度8 km，其中漾濞臺(tái)震中距最小，為7.9 km，東西、南北、垂直向加速度峰值分別為-379.9 cm/s2、720.3 cm/s2、-448.4 cm/s2，速度峰值分別為30.4 cm/s、-29.8 cm/s、-7.2 cm/s，計(jì)算儀器地震烈度為8.3度。

3)2021-05-22 02:04青海省果洛州瑪多縣(34.59°N，98.34°E)發(fā)生MS7.4地震，震源深度17 km，共有16個(gè)臺(tái)站(表2)記錄到此次強(qiáng)震動(dòng)事件，其中大武臺(tái)震中距175.6 km，東西、南北、垂直向加速度峰值分別為46.0 cm/s2、40.6 cm/s2、-19.1 cm/s2，計(jì)算儀器地震烈度為6.0度。

表2 青?，敹嗟卣鹗录_(tái)站信息

4)人工爆破事件數(shù)據(jù)來自中國(guó)水利水電科學(xué)研究院巖土工程研究所，使用PCB 350B01型加速度計(jì)和PCB 350D02型加速度計(jì)進(jìn)行采集，頻率響應(yīng)為1 Hz～10 kHz，得到共計(jì)342條加速度波形記錄，記錄長(zhǎng)度皆為40 000 m/s2。

由于每條記錄的長(zhǎng)度(或采樣時(shí)間)及采樣頻率并不一致且差異較大，為消除幅值大小對(duì)特征工程的影響，本文統(tǒng)一取全波段數(shù)據(jù)進(jìn)行歸一化處理，并按照排列熵、近似熵、香農(nóng)熵的公式計(jì)算得到地震辨識(shí)實(shí)驗(yàn)的機(jī)器學(xué)習(xí)特征集，其中熵值計(jì)算的嵌入維數(shù)設(shè)定為2，時(shí)間延遲為1，容忍度(標(biāo)準(zhǔn)差std)為0.15，計(jì)算結(jié)果如圖1所示，統(tǒng)計(jì)結(jié)果見表3。設(shè)計(jì)多個(gè)辨識(shí)子實(shí)驗(yàn)，其中訓(xùn)練集比例分別為10%、20%、40%、60%、80%，皆為完全隨機(jī)抽取。

圖1 熵值計(jì)算結(jié)果Fig.1 Calculation results of entropy

表3 樣本特征集統(tǒng)計(jì)結(jié)果

結(jié)合圖1和表3可知，排列熵、近似熵、香農(nóng)熵等指標(biāo)均具有明顯的特征區(qū)分效果：在排列熵指標(biāo)上，非天然地震信號(hào)的熵值分布具有明顯的平穩(wěn)特性，標(biāo)準(zhǔn)差為0.004 9；而天然地震信號(hào)熵值則有較大的起伏變化，較不規(guī)律，標(biāo)準(zhǔn)差為0.424 6；在近似熵與香農(nóng)熵指標(biāo)上，非天然地震信號(hào)與天然地震的熵值都具有振蕩現(xiàn)象，標(biāo)準(zhǔn)差均超過0.2，但非天然地震信號(hào)的熵值期望皆比天然地震信號(hào)的大，證明該熵可用于區(qū)分二者波形。

采用機(jī)器學(xué)習(xí)混淆矩陣中的準(zhǔn)確率(accuracy)、召回率(recall)、特效度(specificity)、精確度(precision)及F-measure等評(píng)價(jià)指標(biāo)，衡量LS-SVM模型分類的有效性，核函數(shù)采用RBF函數(shù)，分類子實(shí)驗(yàn)相關(guān)結(jié)果如表4所示。

表4 LSSVM辨識(shí)實(shí)驗(yàn)結(jié)果

以訓(xùn)練集與測(cè)試集數(shù)據(jù)分別占比總數(shù)據(jù)量80%和20%為例，地震事件LSSVM模型二分類結(jié)果如圖2所示，計(jì)算得到的核函數(shù)參數(shù)為2.071 1，懲罰因子為2.307，并得到一條劃分事件波形類別較明顯的非線性超平面。LSSVM模型的輸入往往只需要二維特征矩陣即可，若為多維輸入矩陣，可由模型選擇2條特征向量繪制二分類平面結(jié)果，并自動(dòng)添加與核函數(shù)類型相匹配的超平面(如采用線性核函數(shù)得到的超平面多為直線，采用多項(xiàng)式核函數(shù)得到的超平面多為光滑曲線，基于RBF核函數(shù)得到的超平面多為封閉曲線)。如圖2所示，橫坐標(biāo)和縱坐標(biāo)分別為測(cè)試集的排列熵與香農(nóng)熵，由超平面分割開的淺灰色地帶代表LSSVM超平面劃分的非天然地震分類區(qū)，深灰色地帶為L(zhǎng)SSVM超平面劃分的天然地震分類區(qū)。超平面近似為凹曲線，這與非天然地震事件在分類平面上的特征映射散點(diǎn)在凹曲線上側(cè)且較為聚集的現(xiàn)象保持一致，說明模型對(duì)二者的辨識(shí)較為合理。

圖2 LSSVM模型二分類結(jié)果Fig.2 Two-class results of LSSVM model

為獲悉LSSVM模型在同類機(jī)器學(xué)習(xí)模型中的優(yōu)勢(shì)，選取決策樹、樸素貝葉斯、二次判別分析(QDA)、線性判別分析(LDA)及集成學(xué)習(xí)中的LogitBoost與RobustBoost等作為辨識(shí)實(shí)驗(yàn)對(duì)照方法，其中訓(xùn)練集與測(cè)試集的數(shù)據(jù)量比例均設(shè)定為4∶1，具體結(jié)果見表5。

表5 機(jī)器學(xué)習(xí)模型辨識(shí)效果對(duì)比

由表5可知，在相同的訓(xùn)練集與測(cè)試集比例設(shè)置情況下，QDA、LDA、樸素貝葉斯、決策樹、LogitBoost、RobustBoost等6種機(jī)器方法的準(zhǔn)確率、召回率、特效度、精確度及F-measure值均小于或等于本文方法，證明3種熵結(jié)合最小二乘向量機(jī)的地震辨識(shí)方法具有較好的分類效果，對(duì)天然地震事件與非天然地震事件波形的區(qū)分皆達(dá)到預(yù)期，在同類方法中也有一定的先進(jìn)性。

4 結(jié) 語

本文通過對(duì)預(yù)處理后的全波段波形信號(hào)提取排列熵、近似熵、香農(nóng)熵等熵值，組建三維樣本特征矩陣，并采用最小二乘支持向量分類機(jī)方法進(jìn)行地震辨識(shí)，從而實(shí)現(xiàn)對(duì)地震事件的有效識(shí)別，得到以下結(jié)論：

1)將排列熵、近似熵、香農(nóng)熵引入地震識(shí)別中，具有一定新穎性，其特征提取結(jié)果也具有明顯的事件區(qū)分效果。

2)LSSVM模型收斂速度較快，在訓(xùn)練量較少的情況下依然可以保持較理想的識(shí)別效果，且整體識(shí)別效果要優(yōu)于QDA、LDA、樸素貝葉斯、決策樹、LogitBoost及RobustBoost等經(jīng)典機(jī)器學(xué)習(xí)方法，基于準(zhǔn)確率、召回率、特效度、精確度、F-measure等指標(biāo)的多重評(píng)價(jià)結(jié)果也證明LSSVM模型在地震識(shí)別領(lǐng)域具有一定的優(yōu)越性。

3)基于RBF核函數(shù)的LSSVM模型性能會(huì)受懲罰因子、RBF等參數(shù)的影響，在樣本學(xué)習(xí)和目標(biāo)泛化方面仍有一定的改進(jìn)空間，下一步將結(jié)合智能優(yōu)化算法對(duì)其正則化超參數(shù)進(jìn)行尋優(yōu)，增強(qiáng)LSSVM模型的穩(wěn)健性，并提高其對(duì)含噪地震樣本的非線性求解能力。另外，本文研究?jī)H限于識(shí)別地震波形，若用于地震預(yù)警等，需要考慮辨識(shí)模型的運(yùn)算速度與預(yù)警時(shí)效，將在后續(xù)研究中對(duì)此進(jìn)行補(bǔ)充。

致謝：中國(guó)地震局工程力學(xué)研究所為本文提供數(shù)據(jù)支持，中國(guó)水利水電科學(xué)研究院巖土工程研究所提供爆破數(shù)據(jù)，在此一并表示感謝!