地震波形的希爾伯特黃變換特征提取與隱馬爾科夫模型識(shí)別研究

周海軍，李 磊

滁州學(xué)院電子與電氣學(xué)院,滁州,239000

隨著科學(xué)技術(shù)的發(fā)展，地震震動(dòng)記錄儀能檢測(cè)到頻率更高或更低、能量更微弱的地下震動(dòng)信號(hào)，實(shí)際記錄到的信號(hào)往往同時(shí)包含多種不同來(lái)源的同時(shí)發(fā)生或發(fā)生時(shí)刻接近的多個(gè)天然地震，或其他如礦塌、化學(xué)爆炸及核爆等人工震源。檢測(cè)到的震動(dòng)信號(hào)本身是非平穩(wěn)、非線性的時(shí)變信號(hào)，分析這種信號(hào)需要了解某一頻率的時(shí)間信息或者某一時(shí)間的頻率信息。傳統(tǒng)的傅里葉變換、瓦格納-Ville分布和小波變換在處理這種震動(dòng)信號(hào)時(shí)存在局限性。美國(guó)科學(xué)家提出希爾伯特黃變換(HHT)分析方法，其擴(kuò)張的基礎(chǔ)是自適應(yīng)性，能通過(guò)譜分析得到瞬時(shí)頻率，從而得到信號(hào)的頻率時(shí)間分布信息。

1 地震和爆破波形的差異性

天然地震和人工爆破都發(fā)生在地下，以震動(dòng)波的形式迅速釋放能量，從時(shí)域和頻域上分析波形，發(fā)現(xiàn)地震和爆破存在諸多差異[1]。

(1)爆破源持續(xù)的時(shí)間一般在3秒以內(nèi)，是極短時(shí)間內(nèi)能量膨脹釋放，P波初動(dòng)方向向上。天然地震源持續(xù)的時(shí)間在10秒以上，是板塊位置移動(dòng)，P波方向可能向上，也可能向下。

(2)天然地震發(fā)生在地底幾公里，甚至幾十公里深處，介質(zhì)密度高，波形衰減慢，頻率成分較多。人工爆破多發(fā)生在地表淺處，土質(zhì)松軟，波形衰減快，頻率成分較少。

(3)相同介質(zhì)中，低頻信號(hào)比高頻信號(hào)衰減慢，因此，天然地震波形的頂峰頻率高于人工爆破的頂峰頻率。

(4)天然地震P波和S波的振幅比值要小于人工爆破。

基于這些差異，國(guó)內(nèi)外在天然地震和人工爆破的識(shí)別上進(jìn)行了深入研究，提出震動(dòng)波形的多種特征[2]，諸如：P波和S波的振幅比、AR模型系數(shù)、震源深度、倒譜、S波的對(duì)數(shù)譜以及自相關(guān)系數(shù)等，也提出多種方法判斷震源類型，諸如：人工神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)和隱馬爾科夫模型(HMM)等。

2 HHT基本原理

HHT方法是一種非線性非平穩(wěn)數(shù)據(jù)的分析方法，特別是對(duì)時(shí)頻能量的分析，由經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和Hilbert譜分析(HSA)兩部分組成[3]。

2.1 震動(dòng)信號(hào)的經(jīng)驗(yàn)?zāi)B(tài)分解

大多數(shù)的震動(dòng)數(shù)據(jù)都不是本征模函數(shù)，包含多個(gè)波動(dòng)模式，因此，直接對(duì)原始信號(hào)Hilbert變換不能完全表征地震信號(hào)的頻率特性，這時(shí)需要對(duì)信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解獲得本征模函數(shù)(IMF)。

原始震動(dòng)信號(hào)為x(t)，找到信號(hào)的所有局部極值點(diǎn)，利用三次樣條插值函數(shù)和極大值點(diǎn)，擬合形成地震信號(hào)的上包絡(luò)線ymax(t)，利用三次樣條插值函數(shù)和極小值點(diǎn)，擬合形成地震信號(hào)的下包絡(luò)線ymin(t)，計(jì)算上下包絡(luò)線的均值得到m1(t)：

(1)

用x(t)減去m1(t)得到h1(t)：

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

(2)

若h1(t)還存在負(fù)的局部極大值或正的局部極小值，則需要繼續(xù)按照上面過(guò)程進(jìn)行篩選。這時(shí)將h1(t)當(dāng)作原始地震信號(hào)，參照公式(1)和(2)計(jì)算，直到滿足本征模函數(shù)條件，得到的h11(t)就是第一個(gè)IMF分量。

r1(t)=x(t)-h11(t)

(3)

將r1(t)當(dāng)作新的信號(hào)，重復(fù)上面三個(gè)公式，得到多個(gè)IMF分量[4-5]，直到函數(shù)單調(diào)或者小于預(yù)設(shè)值，分解結(jié)束。圖1是某地震臺(tái)站記錄的原始天然

圖1 天然地震的原始波形

圖2 天然地震的EMD分解波形

地震波形，長(zhǎng)度是10 000個(gè)點(diǎn)，前期和后期有部分無(wú)效信息，以P波初至點(diǎn)為起始點(diǎn)，截至2 000個(gè)點(diǎn)長(zhǎng)度的波形作為EMD分解波形。圖2是2 000點(diǎn)長(zhǎng)度的天然地震波形經(jīng)過(guò)EMD分解之后的9個(gè)IMF分量波形。

圖3是某地震臺(tái)站記錄的原始人工爆破波形，長(zhǎng)度是10 000個(gè)點(diǎn)，前期和后期有些無(wú)效信息，以P波初至點(diǎn)為起始點(diǎn)，截至2 000個(gè)點(diǎn)長(zhǎng)度的波形作為EMD分解波形。圖4是2 000點(diǎn)長(zhǎng)度的人工爆破波形經(jīng)過(guò)EMD分解之后的9個(gè)IMF分量波形。

圖3 人工爆破的原始波形

圖4 人工爆破的EMD分解波形

2.2 Hilbert譜分析

(4)

其中，P為柯西主值[4]。

瞬時(shí)振幅：

(5)

瞬時(shí)相位：

(6)

由θi(t)求導(dǎo)可以得到瞬時(shí)頻率：

(7)

圖5是圖2的IMF分量的Hilbert變換波形，

圖5 天然地震的Hilbert變換

圖6 人工爆破的Hilbert變換

圖6是圖4的IMF分量的Hilbert變換波形。這樣地震信號(hào)x(t)表示為：

(8)

從公式(8)中得知，地震信號(hào)可利用“時(shí)間-頻率-幅值”的三維圖描述，這種地震信號(hào)幅值的時(shí)頻特性就是Hilbert譜。

3 數(shù)據(jù)集選擇和特征提取

3.1 數(shù)據(jù)集選擇

地震記錄儀采集到的地震數(shù)據(jù)包括垂直UD、東西EW、南北NS三個(gè)方向，其中垂直方向UD的信號(hào)更能反映地震信號(hào)的波形特征，實(shí)驗(yàn)選取了首都圈地震帶記錄的天然地震和人工爆破事件，事件波形數(shù)據(jù)格式：EVT，記錄事件2002年5月到2007年12月，震級(jí)小于2.7，天然事件35個(gè)，人工爆破事件27個(gè)。110個(gè)臺(tái)站記錄垂直方向天然地震數(shù)據(jù)3 805份，人工爆破數(shù)據(jù)2 849份，部分臺(tái)站未記錄數(shù)據(jù)。再由人工選擇波形較好的120個(gè)人工爆破和120個(gè)天然地震數(shù)據(jù)。為了消除震級(jí)對(duì)地震信號(hào)識(shí)別的影響，通常對(duì)數(shù)據(jù)進(jìn)行歸一化，將震動(dòng)信號(hào)變?yōu)?到1之間的數(shù)據(jù)。

3.2 分幀與加窗

地震信號(hào)在短時(shí)間內(nèi)可認(rèn)為是平穩(wěn)的信號(hào)，將每個(gè)短時(shí)間內(nèi)采集的信號(hào)看成一幀。這樣地震信號(hào)可被分為m段，處理數(shù)據(jù)時(shí)只需要對(duì)分段的信號(hào)進(jìn)行處理。通常利用窗函數(shù)來(lái)截取地震信號(hào)，在長(zhǎng)度n點(diǎn)的窗函數(shù)中，漢明窗的平滑低通特性較好[6]。漢明窗的定義如下式。

(9)

對(duì)地震信號(hào)加漢明窗分幀，則每幀信號(hào)表示為：

x′(n)=x(n)×w(n)

(10)

3.3 線性預(yù)測(cè)倒譜系數(shù)(LPCC)特征參數(shù)

地震信號(hào)的LPCC參數(shù)是處理后的2 000點(diǎn)長(zhǎng)的地震信號(hào)經(jīng)過(guò)分幀和加漢明窗，然后對(duì)信號(hào)進(jìn)行Z變換，取變換后的對(duì)數(shù)，再進(jìn)行逆Z變換，從得到線性預(yù)測(cè)系數(shù)變換而成。若線性預(yù)測(cè)系數(shù)的參數(shù)是ak，LPCC的參數(shù)是ck，則：

(11)

其中，k為L(zhǎng)PCC的階數(shù)，p為線性預(yù)測(cè)系數(shù)的階數(shù)。

圖7 特征參數(shù)提取過(guò)程

4 HMM模型

HMM是用來(lái)描述一個(gè)隱含未知參數(shù)的馬爾可夫過(guò)程。利用可觀察的參數(shù)確定該過(guò)程的隱含參數(shù)，然后利用這些參數(shù)作進(jìn)一步分析。天然地震和人工爆破模型建立和識(shí)別采用的是連續(xù)高斯混合密度隱馬爾可夫模型(CGHMM)。

4.1 震源信號(hào)中CGHMM的初始化

4.2 模型建立

4.3 震源類型識(shí)別

隨機(jī)選取數(shù)據(jù)集中的20個(gè)天然地震數(shù)據(jù)和20個(gè)人工爆破數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，建立天然地震和人工爆破HMM模型。再?gòu)挠嘞碌奶烊坏卣饠?shù)據(jù)集和人工爆破數(shù)據(jù)集中分別選取20個(gè)數(shù)據(jù)作為識(shí)別數(shù)據(jù)進(jìn)行識(shí)別；HMM模型的狀態(tài)個(gè)數(shù)為4，每個(gè)狀態(tài)有3個(gè)高斯函數(shù)；每個(gè)幀的長(zhǎng)度為64，幀移為50；將正確識(shí)別的個(gè)數(shù)除以總的識(shí)別數(shù)據(jù)的個(gè)數(shù)作為識(shí)別單次概率。按照上述操作，隨機(jī)10次求得平均識(shí)別率見(jiàn)表1。

表1 LPCC參數(shù)的識(shí)別率

5 總結(jié)與分析

雖然地震信號(hào)是非平穩(wěn)、非線性的隨機(jī)信號(hào)，但在短時(shí)間內(nèi)可看成是平穩(wěn)的、線性的信號(hào)。因此，本文采用了分幀和加窗的方式來(lái)得到短時(shí)信號(hào)特性，能較好地提高天然地震和人工爆破的識(shí)別率。

國(guó)內(nèi)外在利用特征識(shí)別震源信號(hào)時(shí)，一般采用多特征識(shí)別，識(shí)別率在90%以上。根據(jù)地震信號(hào)的波形特征，實(shí)驗(yàn)利用語(yǔ)言信號(hào)的識(shí)別方法識(shí)別地震信號(hào)。由于地震信號(hào)存在P波、S波和面波，有兩個(gè)尖峰信號(hào)，故通過(guò)HHT處理，選擇適當(dāng)?shù)木S度，通過(guò)單一特征識(shí)別震源類別，識(shí)別率在85%，利用同一波形的多個(gè)分量，能將識(shí)別率提升到93.1%。

從表1可知：HMM模型識(shí)別率隨著LPCC特征的維數(shù)增大而增大，達(dá)到一定程度時(shí)，識(shí)別率反而會(huì)減小。第一個(gè)IMF分量波形信息更豐富，LPCC特征參數(shù)識(shí)別率在三個(gè)分量中最高。利用三個(gè)IMF分量各自的識(shí)別結(jié)果綜合判斷，有兩個(gè)或兩個(gè)以上的認(rèn)為是某種震源，系統(tǒng)就確定是這種震源類型。震源類型的識(shí)別實(shí)質(zhì)是對(duì)震動(dòng)性質(zhì)和傳播路徑的分析判斷。文中震源選擇的區(qū)域是首都圈地震帶，選取的區(qū)域范圍較小，今后應(yīng)進(jìn)一步對(duì)我國(guó)及周邊大型區(qū)域的較大震級(jí)的震動(dòng)信號(hào)進(jìn)行識(shí)別研究。