基于Hilbert-Huang變換在九寨溝、寧強地震的研究

張 健， 喬 波， 何 楊， 王衛平

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

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

1 HHT的實現過程

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

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

(2)信號的極大值點構成的上包絡與極小值點構成的下包絡關于時間軸對稱。

基于這樣的想法，他們提出了經驗模式分解(EMD)：可以用經驗模態分解方法將信號的固有模態篩選出來的這種方法具體就是個篩選過程，實現振動模式的提取。

EMD分解過程可概括如下：

①給定信號x(t)，求出所有信號中的每個區域內的極大值和極小值，為更好保留原始信號的性質，極大值被認為是時間序列中的某個時刻的值，它只需要滿足大于前、后時刻的值就可以了。局部極小值的獲取同理。然后，使用三次樣條插值函數進行擬合，獲得上包絡線xmax(t)和下包絡線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)替換為原始信號，重復(1)～(3)的步驟，直到滿足SD條件為止。SD的表達式為：

(1)

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

(2)

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

完成了經驗模式分解過程就得到了所有可提取的IMF，在Hilbert變換的基礎上計算瞬時頻率就可以了。

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

(3)

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

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

(4)

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

(5)

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

(6)

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

2 HHT在數據分析中的初步運用

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

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

表1 各臺震中距統計表

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

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

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

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

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

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

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

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

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

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

第四階段：第107～401 s之間，在此階段，地震波主要持續階段已經結束。從分解到的波形來看，這個階段含有一定振幅，但幅值很小，頻率時大時小，變化不穩定。根據該臺的震中距以及所處的地質構造判斷，產生這樣的波形是由于地震波在當地的地質體中發生散射，或者是地震波的尾波造成。

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

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

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

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

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

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

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

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

表2 各臺九寨溝地震反應相關參數匯總表

續表2

臺站名稱主要振動持續時間/s最大振幅值/×106瞬時頻率分布范圍/Hz西鄉臺46111.50.0313～2.5鎮巴臺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譜

具體各臺對該次地震的反應情況見下表3。可以看出在該次地震過程中地震主要持續時間較短，頻率成份分布廣。

表3 各臺寧強地震反應相關參數匯總表

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

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

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

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

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

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

5 結論

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