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

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

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

1 HHT的實現(xiàn)過程

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

(1)信號中的極值點的數(shù)目與過零點的數(shù)目相等或者至多相差一個；

(2)信號的極大值點構(gòu)成的上包絡(luò)與極小值點構(gòu)成的下包絡(luò)關(guān)于時間軸對稱。

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

EMD分解過程可概括如下：

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

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

(1)

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

(2)

即原始序列是由n個IMF與一個趨勢項r(t)組成。

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

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

(3)

因為希爾伯特變換的振幅和頻率都是和時間有關(guān)，用三維立體圖形表達幅值、頻率和時間的關(guān)系，就可以得到Hilbert譜H(w,t)。

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

(4)

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

(5)

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

(6)

希爾伯特能量譜表達了每個頻率在整個時間長度內(nèi)所累積的能量。

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

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

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

表1 各臺震中距統(tǒng)計表

其中第二列為各臺到九寨溝地震震中距，第三列為各臺到寧強地震震中距。

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

圖1 青木川臺記錄的垂直向九寨溝地震波

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

將地震波信號分解為共12道IMF分量，對以上各分量進行希爾伯特變換后可以得出每道分量的瞬時頻率，如圖3所示。

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

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

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

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

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

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

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

圖2 青木川臺記錄的地震信號分解

圖3 每道IMF分量的瞬時頻率

圖4 九寨溝地震青木川臺Hilbert譜

3 地震波頻率在時間域中的分布

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

圖5 九寨溝地震青木川臺三維Hilbert譜

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

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

續(xù)表2

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

同理，我們以寧強臺為例，繪制寧強地震的三維Hilbert譜，見圖6。

圖6 寧強地震寧強臺三維Hilbert譜

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

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

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

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

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

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

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

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

5 結(jié)論

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