基于EMD的IMF時(shí)域統(tǒng)計(jì)特征提取及其應(yīng)用于震動(dòng)事件源類型識(shí)別研究

薛思敏, 黃漢明, 施佳鵬, 袁雪梅, 黎炳君(廣西師范大學(xué)計(jì)算機(jī)科學(xué)與信息工程學(xué)院, 廣西 桂林 54100)

0 引言

人工爆炸和天然地震振動(dòng)產(chǎn)生的波形非常相似。各地震觀測(cè)臺(tái)站記錄的事件波形數(shù)據(jù)含有人工爆炸、泥石流、山體隧道等其他事件的波形信息,這些混有非天然地震事件的波形記錄會(huì)混淆地震目錄,若對(duì)這些波形處理不當(dāng)會(huì)對(duì)地震學(xué)的研究造成影響。如何快速、準(zhǔn)確地從震源事件震動(dòng)波形中識(shí)別出天然地震事件與人工爆炸事件是地震學(xué)研究及相關(guān)應(yīng)用領(lǐng)域亟待解決的問(wèn)題。

天然地震與人工爆炸信號(hào)都屬于非線性、非平穩(wěn)信號(hào)。自二十世紀(jì)五十年代開(kāi)始,國(guó)內(nèi)外學(xué)者進(jìn)行了許多提取蘊(yùn)藏在噪聲時(shí)間序列的該信號(hào)的特征來(lái)進(jìn)行天然地震和非天然地震分類識(shí)別方法研究,并提出了多種識(shí)別依據(jù),如希爾伯特-黃變換譜在地震信號(hào)分析中的應(yīng)用[1],將HHT與Fourier分析進(jìn)行比較,顯示了HHT方法的優(yōu)越性;基于經(jīng)驗(yàn)?zāi)B(tài)分解的地震P波初至自動(dòng)識(shí)別研究[2],將前三個(gè)IMF分量和趨勢(shì)項(xiàng)合成新的地震信號(hào),采用遞歸STA/LTA方法識(shí)別地震P波初至點(diǎn),其精確度較高;天然地震與人工爆炸波形信號(hào)HHT特征提取和SVM識(shí)別研究[3],提取分解后前3個(gè)IMF分量的最大幅值對(duì)應(yīng)周期和倒譜平均值,采用SVM分類識(shí)別,識(shí)別率較好;基于聚類經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)的汶川MS8.0強(qiáng)震動(dòng)記錄時(shí)頻特性分析[4],將數(shù)據(jù)進(jìn)行EEMD和Hilbert變換及譜分析,發(fā)現(xiàn)采用EEMD可以抑制EMD中存在的模態(tài)混疊問(wèn)題;地震波形的HHT特征提取和GMM識(shí)別研究[5],使用GMM對(duì)提取前三個(gè)IMF分量的MFCC特征分別進(jìn)行分類識(shí)別,得出第一個(gè)分量提取的MFCC特征識(shí)別效果最好,可以作為分類識(shí)別有效判據(jù)。

傳統(tǒng)的信號(hào)處理方法如小波變換、傅里葉變換在使用時(shí),都要預(yù)先選取合適的基函數(shù)。經(jīng)驗(yàn)?zāi)B(tài)分解方法依據(jù)數(shù)據(jù)自身的時(shí)間尺度特征來(lái)進(jìn)行信號(hào)分解,無(wú)須預(yù)先設(shè)定任何基函數(shù)。正是由于這個(gè)特點(diǎn),EMD方法在理論上可以應(yīng)用于任何類型的信號(hào)的分解,因而在處理非平穩(wěn)及非線性數(shù)據(jù)上,具有非常明顯的優(yōu)勢(shì)。

1 經(jīng)驗(yàn)?zāi)B(tài)分解和內(nèi)模函數(shù)

希爾伯特-黃變換(Hilbert-Huang Transform,簡(jiǎn)稱HHT)是美籍華人N.E Huang等[6-10]提出的一種信號(hào)分析方法,其核心思想是將時(shí)間序列信號(hào)通過(guò)經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,簡(jiǎn)稱EMD),得到有限個(gè)內(nèi)模函數(shù)(intrinsic mode function,簡(jiǎn)稱IMF),再對(duì)內(nèi)模函數(shù)進(jìn)行希爾伯特變換,從而得到該信號(hào)的希爾伯特譜、時(shí)頻能量譜等,以便對(duì)信號(hào)進(jìn)行分析。該方法對(duì)非線性、非平穩(wěn)信號(hào)有較好的分析和處理效果。

經(jīng)驗(yàn)?zāi)B(tài)分解是假設(shè)任何復(fù)雜的信號(hào)都是由簡(jiǎn)單固有模態(tài)函數(shù)組成,是以局部時(shí)間尺度為基礎(chǔ)的。經(jīng)驗(yàn)?zāi)B(tài)分解是將數(shù)據(jù)集分解為個(gè)數(shù)有限的內(nèi)模函數(shù)IMF的線性疊加。

每一個(gè)IMF不論是線性或者是非線性、平穩(wěn)或是非平穩(wěn)的,都具備以下兩個(gè)條件:在整個(gè)信號(hào)波形中,極點(diǎn)個(gè)數(shù)和過(guò)零率的個(gè)數(shù)必須相等,或最多相差不能超過(guò)一個(gè);在任意時(shí)刻,由局部極大值點(diǎn)形成的上包絡(luò)線和由局部極小值點(diǎn)形成的下包絡(luò)線的平均值為零,即上、下包絡(luò)線相對(duì)于時(shí)間軸局部對(duì)稱。我們可以在上述的基礎(chǔ)上,采用模態(tài)分解方法通過(guò)下面的步驟對(duì)任何信號(hào)x(t)進(jìn)行分解[11-16]:找到信號(hào)x(t)所有的極大值點(diǎn)并利用三次樣條函數(shù)插值擬合成原數(shù)據(jù)序列的上包絡(luò)線v1(t);找出信號(hào)x(t)所有的極小值點(diǎn)并利用三次樣條函數(shù)插值擬合成原數(shù)據(jù)序列的下包絡(luò)線v2(t),并求出其上包絡(luò)線和下包絡(luò)線的平均值m1(t)。其公式為:

(1)

將原數(shù)據(jù)序列減去包絡(luò)平均值后即可得到一個(gè)去掉低頻的新數(shù)據(jù)序列。其公式為:

h1(t)=x(t)-m1(t)

(2)

判斷h1(t)是否為IMF函數(shù)。若不滿足IMF條件,將h1(t)看作新的信號(hào)序列,重復(fù)上述處理過(guò)程,直到h1(t)滿足條件時(shí),記c1(t)為IMF1。其公式為:

c1(t)=h1(t)

(3)

令:

r(t)=x(t)-c1(t)

(4)

看作新的信號(hào)序列,重復(fù)h1(t)的判斷步驟,即可得到不同結(jié)果的IMF分量,直到滿足給定的終止條件時(shí)篩選結(jié)束。原始的數(shù)據(jù)序列即可由這些IMF分量以及一個(gè)剩余項(xiàng)r(t)表示,即:

(5)

式(5)說(shuō)明,可以將信號(hào)x(t)分解成頻率從大到小也即變化由快到慢排列的n項(xiàng)IMF分量與一個(gè)剩余項(xiàng)r(t)之和。由于每一個(gè)IMF分量代表一個(gè)特征尺度的數(shù)據(jù)序列,因此篩選的過(guò)程實(shí)際上是將原來(lái)數(shù)據(jù)序列分解為各種不同特征波動(dòng)序列的疊加。

經(jīng)驗(yàn)?zāi)B(tài)分解的過(guò)程其實(shí)就是一個(gè)篩選的過(guò)程,在篩選的過(guò)程中不僅消除了模態(tài)波形的疊加,而且使波形輪廓更加對(duì)稱。在天然地震和人工爆炸信號(hào)分析中,利用經(jīng)驗(yàn)?zāi)Ｊ椒纸獠恍枰潭ǖ幕瘮?shù),可將原始信號(hào)分解為若干個(gè)頻率自高到低排列的內(nèi)模函數(shù),其分解過(guò)程具有高效性、自適應(yīng)性。

2 KL距離

KL距離能夠有效衡量?jī)蓚€(gè)概率分布之間的距離。若這兩個(gè)概率分別為離散隨機(jī)變量P(xi)和Q(xi),xi是隨機(jī)變量,則KL距離可定義為:

(6)

式中:KL(P(xi)‖Q(xi))表示兩個(gè)概率密度函數(shù)P(xi)與Q(xi)似然比的對(duì)數(shù)期望,也稱相對(duì)熵(relative entropy)。當(dāng)且僅當(dāng)兩概率分布完全相同時(shí),相對(duì)熵KL(P(xi)‖Q(xi))最小且為零。由于相對(duì)熵不具有對(duì)稱性,概率分布P(xi)與Q(xi)的相對(duì)熵距離通常并不等于從Q(xi)到P(xi)的距離。

定義一種對(duì)稱的KL距離來(lái)衡量這兩個(gè)概率分布的相似度[17-19]:

d(P(xi),Q(xi))=KL(P(xi)‖Q(xi))+

KL(Q(xi)‖P(xi))

(7)

式中:P(xi)表示提取天然地震或人工爆炸數(shù)據(jù)特征值按概率分布形式來(lái)表達(dá)的矩陣;Q(xi)表示提取天然地震或人工爆炸數(shù)據(jù)特征均值按概率分布形式來(lái)表達(dá)的矩陣;d(P(xi),Q(xi))表示P(xi)與Q(xi)概率分布之間的相似度。

本文利用KL距離對(duì)微小異常的敏感性,采用對(duì)稱的KL相似度,考慮天然地震數(shù)據(jù)特征值按概率分布形式表達(dá)的矩陣與天然地震數(shù)據(jù)特征均值按概率分布形式表達(dá)的矩陣之間的距離(d1),以及天然地震數(shù)據(jù)特征值按概率分布形式表達(dá)的矩陣與人工爆炸數(shù)據(jù)特征均值按概率分布形式表達(dá)的矩陣之間的距離(d2),將兩個(gè)距離進(jìn)行比較,若d1

3 實(shí)驗(yàn)與結(jié)論

3.1 數(shù)據(jù)集的選取

實(shí)驗(yàn)中所用的天然地震記錄樣本和人工爆炸記錄樣本均選取于三河市2010—2016年發(fā)生的地震記錄。其中有16個(gè)天然地震事件,共131條數(shù)據(jù)記錄;有17個(gè)人工爆炸事件,共186條記錄。在實(shí)際的天然記錄樣本和人工爆炸記錄樣本中有一些記錄采樣點(diǎn)的數(shù)據(jù)幅值為0的,本文作者認(rèn)為是無(wú)效數(shù)據(jù)。最終本文共得到307條天然地震和人工爆炸波形記錄的數(shù)據(jù)作為實(shí)驗(yàn)原始數(shù)據(jù)。

3.2 特征提取

本文對(duì)真實(shí)的天然地震和人工爆炸事件的三分量波形提取特征。我們對(duì)波形數(shù)據(jù)進(jìn)行歸一化預(yù)處理,采用EMD方法將原波形信號(hào)分解為7個(gè)內(nèi)模函數(shù)和1個(gè)殘差函數(shù),對(duì)原波形、每個(gè)內(nèi)模函數(shù)和殘差函數(shù)分別提取26個(gè)時(shí)域統(tǒng)計(jì)特征,組成9個(gè)特征組(Q0、Q1、…、Q8);再計(jì)算7個(gè)內(nèi)模函數(shù)的幅度能量比得到7個(gè)能量比特征,選擇前4個(gè)內(nèi)模函數(shù)的26個(gè)時(shí)域統(tǒng)計(jì)特征中的8個(gè)特征共32個(gè)特征組成一個(gè)有39個(gè)特征的特征組。如圖1為序號(hào)37的天然地震事件HE.XIL臺(tái)站記錄的垂直波形數(shù)據(jù)分解得到的內(nèi)模函數(shù)(IMF)和殘差函數(shù),圖2為序號(hào)147的人工爆炸事件BU.LBG臺(tái)站記錄的垂直波形數(shù)據(jù)分解得到的內(nèi)模函數(shù)(IMF)和殘差函數(shù)。

圖1 天然地震的EMD分解Fig.1 EMD of natural earthquake

本文前述的Q0、Q1、…、Q8特征組中的26個(gè)時(shí)域特征都是某種統(tǒng)計(jì)特征,它們分別是:描述樣本集中程度的量值:算數(shù)均值、中位數(shù)、模態(tài)數(shù)、截尾均值、調(diào)和均值;描述樣本散布程度的量值:四分位差、標(biāo)準(zhǔn)方差、平均絕對(duì)值、各階中心距(3～9階中心距,共7個(gè)中心距);描述樣本分布偏離正態(tài)性的量值:偏斜度、峰度、分位數(shù)(共9個(gè)分位數(shù))。其公式如下:

均值

(8)

調(diào)和均值

(9)

四分位差

(10)

平均絕對(duì)值

(11)

圖2 人工爆炸的EMD分解Fig.2 EMD of artificial explosion

標(biāo)準(zhǔn)方

(12)

中心距值

mi=E(x-E(x))i,i=3,4,…,9

(13)

偏斜度:

(14)

峰度:

(15)

式中:x表示信號(hào)數(shù)據(jù);i表示第i個(gè)數(shù)值;n表示一段地震波段的采樣點(diǎn)數(shù);E(x)表示期望值。

為了避免事件強(qiáng)度差異的影響(事件強(qiáng)度差異會(huì)造成波形幅值范圍顯著不同,從而使事件類型識(shí)別結(jié)果很可能只依賴于事件的強(qiáng)度),計(jì)算這些時(shí)域統(tǒng)計(jì)特征前,原觀測(cè)波形需先進(jìn)行歸一化,即每一份波形取值范圍規(guī)范化為界于[-1,1]。

本文闡述的Q9特征組的39個(gè)特征,分別是:由7個(gè)內(nèi)模函數(shù)的幅度能量比得到7個(gè)能量比特征,如式(16);前4個(gè)內(nèi)模函數(shù)的26個(gè)時(shí)域統(tǒng)計(jì)特征中的8個(gè)特征,分別為算數(shù)均值、中位數(shù)、四分位差、標(biāo)準(zhǔn)方差、偏斜度、峰度、3階中心距及4階中心距。

(16)

3.3 用對(duì)稱KL距離對(duì)天然地震和人工爆炸信號(hào)進(jìn)行識(shí)別

本文采用對(duì)稱KL距離分類識(shí)別方法,首先由于提取的各個(gè)特征矩陣可能存在權(quán)重不等的問(wèn)題,也為了剔除不同記錄波形的樣本點(diǎn)幅值有顯著差異的影響,應(yīng)將特征數(shù)據(jù)集中的所有特征進(jìn)行預(yù)處理,所以以列為單位對(duì)天然地震和人工爆炸分別歸一化處理,采用的公式為:

(17)

式中:X表示m行n列的特征值矩陣;Xmax、Xmin也是m行n列的矩陣,這2個(gè)矩陣任一個(gè)的每一行與同矩陣的其他行都是完全相同的,Xmin/Xmax矩陣的每一行的各元素分別是X矩陣相應(yīng)列的最小值/最大值。

本文對(duì)這10組特征樣本集以事件為識(shí)別單元進(jìn)行了單組、多組的特征組合事件類型識(shí)別實(shí)驗(yàn),采用隨機(jī)選擇部分(30%,50%,70%,90%)事件即作為訓(xùn)練集,同時(shí)也作為測(cè)試集,對(duì)每次選定的百分比都反復(fù)進(jìn)行了1 000次的實(shí)驗(yàn)(每次實(shí)驗(yàn),16個(gè)地震和17個(gè)爆炸事件必須各被選中一個(gè),還有其他215=32 768和216=65 536個(gè)可能被選中的事件組合)。由于每次都是隨機(jī)劃分得到訓(xùn)練樣本和測(cè)試集樣本,為了更好地識(shí)別結(jié)果以及檢驗(yàn)特征的有效性,本文分別求正確識(shí)別率的1 000次隨機(jī)抽樣的平均值、方差、最大值和最小值四個(gè)統(tǒng)計(jì)參數(shù),其中各統(tǒng)計(jì)量的計(jì)算公式為:

平均值

(18)

方差

σr=std(r1,…r1 000)

(19)

最大值

r(max)=max(r1,…r1 000)

(20)

最小值

r(min)=min(r1,…r1 000)

(21)

其中:r1…r1 000表示每次的特征識(shí)別結(jié)果。其識(shí)別效果統(tǒng)計(jì)量如表1所列,與其他實(shí)驗(yàn)對(duì)比結(jié)果如表2所列。

其中,Q2(90%)和Q2、Q9(90%)表示隨機(jī)選取事件的90%進(jìn)行實(shí)驗(yàn)。Q2(90%)特征組的識(shí)別率高于2016年王永康[20]和2012年畢明霞[3],穩(wěn)定性也高于王永康的實(shí)驗(yàn)結(jié)果;Q2、Q9(90%)特征組的識(shí)別率高于2012年畢明霞的實(shí)驗(yàn)結(jié)果,穩(wěn)定性高于王永康和畢明霞的實(shí)驗(yàn)結(jié)果。

表1 識(shí)別效果的統(tǒng)計(jì)量

表2 實(shí)驗(yàn)對(duì)比結(jié)果

4 結(jié)論與討論

本文實(shí)驗(yàn)所用的數(shù)據(jù)集包含有16個(gè)地震事件的121條波形信號(hào),17個(gè)爆炸事件的186 條波形信號(hào)。采用隨機(jī)選擇部分(30%,50%,70%,90%)事件即作為訓(xùn)練集同時(shí)也作為測(cè)試集,反復(fù)進(jìn)行了1 000次的實(shí)驗(yàn),結(jié)果表明第2個(gè)內(nèi)模函數(shù)提取的時(shí)域統(tǒng)計(jì)特征在選擇90%事件作為訓(xùn)練集和測(cè)試集的識(shí)別效果最好。從表1容易看出,無(wú)論哪個(gè)特征組,訓(xùn)練集合所選的事件數(shù)目占全部事件的比例越高識(shí)別效果越好;Q2、Q9特征組的組合明顯比其中任意一個(gè)特征組的識(shí)別效果更好。本文的實(shí)驗(yàn)所用數(shù)據(jù)集的事件數(shù)目不多,只在原信號(hào)和IMF的時(shí)間域上進(jìn)行了特征提取,但采用了隨機(jī)有回放的采樣實(shí)驗(yàn),結(jié)果比較理想,且具有一定的推廣性,值得在取得更多事件的數(shù)據(jù)集后進(jìn)一步深入探討。另外,IMF也可認(rèn)為是某種頻率分解成分,所以本文的時(shí)域特征里也已蘊(yùn)涵有頻率信息,且分解得到的IMF具有自適應(yīng)濾波性,這很可能就是 2ndIMF具有明顯較好識(shí)別性能的原因。把IMF作為普通信號(hào),再?gòu)钠涮崛☆l域或時(shí)頻特征,這也應(yīng)該值得探討。

致謝:中國(guó)地震局地球物理研究所“國(guó)家數(shù)字測(cè)震臺(tái)網(wǎng)數(shù)據(jù)備份中心”為本研究提供地震波形數(shù)據(jù)[21]。