基于Hilbert-Huang變換在九寨溝、寧強(qiáng)地震的研究

張 健， 喬 波， 何 楊， 王衛(wèi)平

(陜西省地震局，陜西 西安 710000)

隨著地震監(jiān)測(cè)臺(tái)網(wǎng)的數(shù)字化、網(wǎng)絡(luò)化大力發(fā)展，選用新方法實(shí)現(xiàn)地震數(shù)據(jù)精度和采樣率雙高的分析是目前急需的技術(shù)。地震波形是一種隨時(shí)間變化的典型的地球物理數(shù)據(jù)信號(hào)。對(duì)于這類(lèi)信號(hào)的研究，與傳統(tǒng)的頻譜分析方法相比，Hilbert-Huang變換(簡(jiǎn)稱(chēng)HHT)不受Heisenberg測(cè)不準(zhǔn)原理的制約，它完全擺脫了線(xiàn)性和平穩(wěn)性的制約，非常適合用于分析非線(xiàn)性非平穩(wěn)信號(hào)。

1 HHT的實(shí)現(xiàn)過(guò)程

從物理角度來(lái)看，并非所有的信號(hào)都能用瞬時(shí)頻率來(lái)分析，只有當(dāng)信號(hào)只包含一種振動(dòng)模式，而沒(méi)有疊加波的情況時(shí)才能利用瞬時(shí)頻率分析。但是根據(jù)一般的瞬時(shí)頻率求解方法，求得的瞬時(shí)頻率既有正的也有負(fù)的。而負(fù)的頻率是沒(méi)有任何物理意義。為了獲得有物理意義的瞬時(shí)頻率，N.E.Huang等人研究并認(rèn)為任何一種信號(hào)都可以由基礎(chǔ)信號(hào)即本征模態(tài)函數(shù)IMF(Intrinsic Mode Function)組成，該本征模態(tài)函數(shù)是一種局部限制而非全局限制的函數(shù)形式，實(shí)現(xiàn)這一目的要用一種特別的分解方法，將信號(hào)分解為瞬時(shí)頻率能夠合理定義的分量形式；因此定義了一種基于信號(hào)的局部特性的函數(shù)—IMF。根據(jù)以上定義，一個(gè)本征模態(tài)函數(shù)需滿(mǎn)足以下兩個(gè)條件：

(1)信號(hào)中的極值點(diǎn)的數(shù)目與過(guò)零點(diǎn)的數(shù)目相等或者至多相差一個(gè)；

(2)信號(hào)的極大值點(diǎn)構(gòu)成的上包絡(luò)與極小值點(diǎn)構(gòu)成的下包絡(luò)關(guān)于時(shí)間軸對(duì)稱(chēng)。

基于這樣的想法，他們提出了經(jīng)驗(yàn)?zāi)Ｊ椒纸?EMD)：可以用經(jīng)驗(yàn)?zāi)B(tài)分解方法將信號(hào)的固有模態(tài)篩選出來(lái)的這種方法具體就是個(gè)篩選過(guò)程，實(shí)現(xiàn)振動(dòng)模式的提取。

EMD分解過(guò)程可概括如下：

①給定信號(hào)x(t)，求出所有信號(hào)中的每個(gè)區(qū)域內(nèi)的極大值和極小值，為更好保留原始信號(hào)的性質(zhì)，極大值被認(rèn)為是時(shí)間序列中的某個(gè)時(shí)刻的值，它只需要滿(mǎn)足大于前、后時(shí)刻的值就可以了。局部極小值的獲取同理。然后，使用三次樣條插值函數(shù)進(jìn)行擬合，獲得上包絡(luò)線(xiàn)xmax(t)和下包絡(luò)線(xiàn)xmin(t)；

②計(jì)算分解原始信號(hào)所需的均值m(t)=[xmax(t)+xmin(t)]/2；③用原始信號(hào)x(t)減去均值m(t)，得到第一組h(t)=x(t)-m(t)；但是，h(t)不一定就是一個(gè)IMF，如果h(t)不滿(mǎn)足IMF兩個(gè)條件，就把h(t)替換為原始信號(hào)，重復(fù)(1)～(3)的步驟，直到滿(mǎn)足SD條件為止。SD的表達(dá)式為：

(1)

通常情況下SD<0.3。這時(shí)滿(mǎn)足固有模態(tài)函數(shù)條件的h(t)作為一個(gè)IMF，令c1(t)=h(t)，至此第一個(gè)IMF已經(jīng)成功的提取了。剩余的r(t)=x(t)-c1(t)仍然包含有用信息，因此可以把它看成新原始信號(hào)；重復(fù)上述過(guò)程，依次得到第二個(gè)c2(t)，第三個(gè)c3(t)，…，當(dāng)r(t)滿(mǎn)足函數(shù)單調(diào)性或?yàn)槌?shù)條件時(shí)，終止篩選，這樣就完成了提取IMF。由此可得x(t)的表達(dá)式為：

(2)

即原始序列是由n個(gè)IMF與一個(gè)趨勢(shì)項(xiàng)r(t)組成。

完成了經(jīng)驗(yàn)?zāi)Ｊ椒纸膺^(guò)程就得到了所有可提取的IMF，在Hilbert變換的基礎(chǔ)上計(jì)算瞬時(shí)頻率就可以了。

用下面方式表示信號(hào)的前提是對(duì)固有模態(tài)函數(shù)進(jìn)行Hilbert變換，即：

(3)

因?yàn)橄柌刈儞Q的振幅和頻率都是和時(shí)間有關(guān)，用三維立體圖形表達(dá)幅值、頻率和時(shí)間的關(guān)系，就可以得到Hilbert譜H(w,t)。

如果把H(w,t)對(duì)時(shí)間積分，就得到希爾伯特邊際譜h(w)：

(4)

可以用式(4)定義的希爾伯特瞬間的能量：

(5)

事實(shí)上，如果振幅的平方對(duì)時(shí)間積分，可以得到希爾伯特能量譜：

(6)

希爾伯特能量譜表達(dá)了每個(gè)頻率在整個(gè)時(shí)間長(zhǎng)度內(nèi)所累積的能量。

2 HHT在數(shù)據(jù)分析中的初步運(yùn)用

HHT不僅在濾波方面有出色的表現(xiàn)，還在數(shù)字信號(hào)處理方面有獨(dú)特的方法，根據(jù)公式(4)與公式(5)可以得出，HHT還可以研究在像地震這種非平穩(wěn)信號(hào)在發(fā)展過(guò)程中頻率在時(shí)間上的分布以及能量隨時(shí)間的具體變化情況等一系列隱藏在地震波背后的信息。

以2017年8月8日四川省北部阿壩州九寨溝縣發(fā)生的Ms7.0級(jí)地震(以下簡(jiǎn)稱(chēng)九寨溝地震)和2018年9月12日陜西漢中市寧強(qiáng)縣Ms5.3級(jí)地震(以下簡(jiǎn)稱(chēng)寧強(qiáng)地震)為例，研究HHT方法在地震波形中的更深一步的運(yùn)用。通過(guò)收集陜西省南部地區(qū)15個(gè)地震臺(tái)站(見(jiàn)表1)所記錄到的2次地震的地震波形數(shù)據(jù)進(jìn)行分析。

表1 各臺(tái)震中距統(tǒng)計(jì)表

其中第二列為各臺(tái)到九寨溝地震震中距，第三列為各臺(tái)到寧強(qiáng)地震震中距。

我們首先以青木川臺(tái)所記錄到的垂直向九寨溝地震波數(shù)據(jù)為例，如圖1所示的波形進(jìn)行分解分析。

圖1 青木川臺(tái)記錄的垂直向九寨溝地震波

使用EMD方法對(duì)上圖的波形數(shù)據(jù)進(jìn)行分解得到圖2所示。

將地震波信號(hào)分解為共12道IMF分量，對(duì)以上各分量進(jìn)行希爾伯特變換后可以得出每道分量的瞬時(shí)頻率，如圖3所示。

求出各分量的瞬時(shí)頻率，根據(jù)公式(3)計(jì)算得出希爾伯特譜，如圖4所示。

通過(guò)圖2以及圖4綜合分析，我們認(rèn)為地震發(fā)展過(guò)程可分為以下幾個(gè)過(guò)程：

第一階段：0～17.7 s，此時(shí)處于地震發(fā)生前的地脈動(dòng)階段。此時(shí)地震還沒(méi)有發(fā)生，在圖中表現(xiàn)為振動(dòng)頻率小，主要頻率成分為低頻地脈動(dòng)與背景噪聲。

第二階段：17.7～61 s，此階段地震波首波已經(jīng)到達(dá)臺(tái)站，圖4中用顏色表示其振幅大小，顏色由藍(lán)到紅表示振幅由小到大，清晰的描述了各頻率成分在時(shí)間軸的分布特征以及能量分布特征。由圖4分析出，在第20.3 s為P波的初動(dòng)時(shí)刻。此外，地震信號(hào)頻率分布在0.3～8 Hz之間，屬于近震，所以首先可以判斷Pn震相是首先到達(dá)的。

第三階段:第61～107 s，通過(guò)圖2可以看到此階段波形幅值與周期增大，為S波出現(xiàn)并發(fā)展的階段，也是地震波能量集中釋放的階段。而且通過(guò)圖4，可以分析出，在第80～100 s之間是地震波發(fā)展到頂峰階段，振幅達(dá)到最大但其頻率最低，為0.06～0.1 Hz。符合S波頻率小，振幅大的特征。

第四階段：第107～401 s之間，在此階段，地震波主要持續(xù)階段已經(jīng)結(jié)束。從分解到的波形來(lái)看，這個(gè)階段含有一定振幅，但幅值很小，頻率時(shí)大時(shí)小，變化不穩(wěn)定。根據(jù)該臺(tái)的震中距以及所處的地質(zhì)構(gòu)造判斷，產(chǎn)生這樣的波形是由于地震波在當(dāng)?shù)氐牡刭|(zhì)體中發(fā)生散射，或者是地震波的尾波造成。

第五階段：第401 s之后，為大地回歸平穩(wěn)階段。隨著地震波及其產(chǎn)生的尾波逐漸結(jié)束， 又回到正常地脈動(dòng)記錄階段。

圖2 青木川臺(tái)記錄的地震信號(hào)分解

圖3 每道IMF分量的瞬時(shí)頻率

圖4 九寨溝地震青木川臺(tái)Hilbert譜

3 地震波頻率在時(shí)間域中的分布

在對(duì)15個(gè)臺(tái)站記錄的波形數(shù)據(jù)進(jìn)行分解的基礎(chǔ)之上，對(duì)其每一道IMF分量根據(jù)公式(3)進(jìn)行Hilbert變換，然后根據(jù)公式(4)進(jìn)行計(jì)算，得出2次地震15個(gè)臺(tái)全部的Hilbert譜，我們以青木川臺(tái)為例，用三維圖像呈現(xiàn)出來(lái)用于表達(dá)幅值、頻率和時(shí)間的關(guān)系。如圖5所示。

圖5 九寨溝地震青木川臺(tái)三維Hilbert譜

由圖5可見(jiàn)，三維的Hilbert譜清晰刻畫(huà)出接收到的地震波頻率、振幅隨時(shí)間的變化特征，振幅幅值隨著時(shí)間迅速降低。具體各臺(tái)對(duì)該次地震反應(yīng)情況見(jiàn)下表2。

表2 各臺(tái)九寨溝地震反應(yīng)相關(guān)參數(shù)匯總表

續(xù)表2

臺(tái)站名稱(chēng)主要振動(dòng)持續(xù)時(shí)間/s最大振幅值/×106瞬時(shí)頻率分布范圍/Hz西鄉(xiāng)臺(tái)46111.50.0313～2.5鎮(zhèn)巴臺(tái)4754.10.0313～2.5石泉臺(tái)4028.10.0313～1.9紫陽(yáng)臺(tái)4784.10.0313～2.3漢陰臺(tái)4933.10.0313～1.8略陽(yáng)臺(tái)2628.80.0313～3.1留壩臺(tái)2468.30.0313～2.8佛坪臺(tái)3735.90.0313～2.2鳳縣臺(tái)30110.50.0313～3.5太白臺(tái)4324.30.0313～2.8眉縣臺(tái)3526.10.0313～2.2

同理，我們以寧強(qiáng)臺(tái)為例，繪制寧強(qiáng)地震的三維Hilbert譜，見(jiàn)圖6。

圖6 寧強(qiáng)地震寧強(qiáng)臺(tái)三維Hilbert譜

具體各臺(tái)對(duì)該次地震的反應(yīng)情況見(jiàn)下表3?？梢钥闯鲈谠摯蔚卣疬^(guò)程中地震主要持續(xù)時(shí)間較短，頻率成份分布廣。

表3 各臺(tái)寧強(qiáng)地震反應(yīng)相關(guān)參數(shù)匯總表

對(duì)比表2與表3，結(jié)合表1可以看出，表2對(duì)應(yīng)的九寨溝地震震級(jí)大，震中距都在190～500 km以?xún)?nèi)，屬于近震，其地震主要振動(dòng)持續(xù)時(shí)間在第190～500s左右的范圍內(nèi)，而且瞬時(shí)頻率范圍在0.0313～4.2 Hz。而對(duì)于表3對(duì)應(yīng)的寧強(qiáng)地震震級(jí)小，震中距在20～300 km以?xún)?nèi)，屬于地方震。其地震主要振動(dòng)持續(xù)時(shí)間大部分在第110～430 s左右的范圍內(nèi)，而且瞬時(shí)頻率范圍在0.0313～12.8 Hz。

通過(guò)對(duì)比表2和表3結(jié)合圖5、圖6可以得出以下結(jié)論：①震中距與瞬時(shí)頻率的范圍呈反比例關(guān)系，隨著震中距增大，頻率范圍逐漸減小。②震級(jí)越大，其地震所對(duì)應(yīng)的主要頻率越集中。

4 地震波能量在頻率域的分布

我們主要以九寨溝地震為例，經(jīng)過(guò)計(jì)算，圖7是15個(gè)臺(tái)的能量譜圖。

圖7 九寨溝地震垂直與水平向能量譜

根據(jù)公式(6)可以得出，該圖譜反映了各個(gè)瞬時(shí)頻率成份在整個(gè)地震波持續(xù)時(shí)間內(nèi)所積累的能量。由圖可以看出地震波能量在隨著瞬時(shí)頻率的增大有一個(gè)現(xiàn)增加到最大值然后迅速衰減直至為零。以圖7為例，在0.0313 Hz處開(kāi)始出現(xiàn)地震波能量然后迅速增大，在0.0938 Hz處達(dá)到最大值，達(dá)到了的能量，在此處也是各個(gè)臺(tái)能量達(dá)到最大。隨后迅速衰減，直到瞬時(shí)頻率大于3 Hz時(shí)，能量衰減至及以下的數(shù)量級(jí)，與最大值相差4個(gè)數(shù)量級(jí)。所以，通過(guò)結(jié)合表3分析我們可以得到，九寨溝地震的優(yōu)勢(shì)瞬時(shí)頻率的范圍在0.0313～2 Hz之間，即該地震主要釋放的能量所對(duì)應(yīng)的頻率范圍。

5 結(jié)論

本文通過(guò)利用HHT方法分別對(duì)九寨溝地震和寧強(qiáng)地震地?cái)?shù)據(jù)進(jìn)行處理分析，利用這種方法可以很好的得到以下結(jié)果：(1)HHT方法很好的解析了接收到的地震波頻率特征、振幅特征、能量特征隨時(shí)間發(fā)展過(guò)程，整個(gè)過(guò)程可以分五個(gè)階段；(2)總結(jié)出兩次地震震中距、幅值和持續(xù)時(shí)間、瞬時(shí)頻率之間的相關(guān)關(guān)系；(3)統(tǒng)計(jì)出地震波優(yōu)勢(shì)瞬時(shí)頻率范圍，并得出各個(gè)頻率所對(duì)應(yīng)的能量大?。壕耪瘻系卣鸬闹饕獌?yōu)勢(shì)頻率為0.0938 Hz，其累計(jì)釋放了約大小的能量，同理可得到寧強(qiáng)地震的主要優(yōu)勢(shì)頻率為0.219 Hz，累計(jì)釋放了大小的能量。