地磁HHT頻譜分析在地震預測中的應用

徐強 張建國 王同利

1）中國石家莊 050021 河北省地震局

2）中國河北 056001 邯鄲地震監測中心站

3）中國北京 100080 北京市地震局

0 引言

地磁觀測用于地震預測研究已有多年，且取得了較好的進展。但地磁相對觀測中哪個分量與地震孕育發生間關系較密切還不太明確。鑒于此，我們利用希爾伯特-黃變換方法（以下簡稱HHT 方法）（Huang et al，1998，2003）的頻譜方法對地磁相對觀測數據進行處理和分析。

HHT 方法是Hilbert Huang 等提出的一種新的信號處理方法。該方法首先對信號進行經驗模態分解（empirical mode decomposition，簡稱EMD）（羅奇峰等，2003），有效地把各種頻率成分以本征模態函數（intrinsic mode function，簡稱IMF）的形式從原始信號中分離出來，然后對上述分離出來的IMF 利用Hilbert 變換進行計算，并對計算結果在整個時間區間內進行積分，最終得到希爾伯特邊際譜（公茂盛等，2003）。類似于通常頻譜分析方法，利用該方法也可以得到研究信號的優勢頻率。安張輝等（2011）利用該方法對地電場觀測數據中受城市軌道交通的干擾進行了數據剝離和分析，將數據受城市鐵路干擾的程度降到了最低。

1 HHT 方法計算原理

利用Hilbert-Huang（HHT）方法分析數據，即通過經驗模態分解，將信號分解成一系列本征模態函數，再通過Hilbert 變換得到三維時頻空間的Hilbert 譜（羅奇峰等，2003；石春香等，2003；段生全等，2005；葸曉宇等，2006，2007；張立等，2007；周摯等，2008，安張輝等，2010）。

對于任意1 個時間序列X(t)，其Hilbert 變換定義為

通過式（1）變換，X(t)、Y(t)可以組合成1 個復數信號Z(t)

為了使瞬時頻率具有意義，進行Hilbert 變換的信號序列必須是單組分的，即時間與頻率間必須是一一對應的關系。經過經驗模態分解所得到的本征模態函數可滿足該要求，對X(t)的n階本征模態函數分別進行Hilbert 變換后，原始時間序列可以表示為

將式（4）稱為X(t)的Hilbert 譜，記為H(ω,t)。式（4）中省略了殘余項Rn(t)。aj(t)和ωj(t)均為時間的函數，把振幅顯示在頻率時間平面上便可得到Hilbert 譜H(ω,t)。如果對H(ω,t)進行時間積分，則可以得到Hilbert 邊際譜

邊際譜通常代表某固定頻率在整個時間區域上所有振幅的總和。

2 震例分析

2.1 汶川MS 8.0 地震

成都地磁臺始建于1971 年，1997 年地磁相對觀測儀器正式觀測，2000 年該臺進行了數字化改造，M15 觀測儀器于2007 年開始運行，觀測數據準確可靠，臺站周邊環境良好，干擾較少。成都地磁臺距汶川地震震中約100 km。

重慶地磁臺始建于1982 年2 月，1992 年成為國家Ⅱ類地磁臺。2000 年增上DI 儀觀測；2007 年6 月、7 月先后安裝2 套GM4 磁通門磁力儀，完成數字化改造。后因臺站環境遭到破壞，2011 年重慶地磁臺停測。重慶地磁臺距汶川地震震中約100 km。

由前述可見，利用HHT 數據處理方法進行分析，不受數據采集頻率和振幅的影響，可隨意定制處理，一般干擾在IMF1—IMF4 就已經分解完畢。因此，在數據干擾無法判別時，可以用IMF5—IMF8 的數據正常分析（限于篇幅，后面圖件略），且不受其他因素的影響。圖1、2 分別為成都地磁臺2007 年7 月至2008 年10 月和重慶地磁臺2007 年9 月至2008 年6 月地磁觀測HHT 頻譜。由圖1、2 可見，頻率基本穩定不變，發生頻率擾動的時刻反映了HHT 所具有的瞬時特性，頻譜圖中較好地顯示出分量優勢頻率為0.01 Hz，且各測向優勢頻率基本一致。由圖1、2 還可見汶川MS8.0 地震前兆異常的同震變化情況，成都地磁臺、重慶地磁臺地磁相對觀測分量不同，HHT 譜能量積累過程也不同，磁偏角D分量震前能量積累過程較集中、明顯。

圖1 成都地磁臺地磁觀測HHT 頻譜（a）D分量；（b）經驗模態分解圖；（c）干擾剔除示意圖；（d）H分量；（e）Z分量Fig.1 HHT spectrum of the geomagnetic station of Chengdu

圖2 重慶地磁臺地磁觀測HHT 頻譜（a）D分量；（b）F分量；（c）H分量；（d）Z分量Fig.2 HHT spectrum of the geomagnetic station of Chongqing

2.2 汶川MS 8.0 地震余震

拉薩地磁臺位于拉薩市西郊“七一”農場內，其海拔高度約3655 m，是中國地震局地磁基本臺網Ⅰ類基準臺。臺站建設于1957 年4 月完工，臺基為礫石沖積層，厚約70 m，周圍無明顯干擾源，觀測區磁場梯度為1 nT/m，分布較均勻。地磁絕對觀測室（2006 年“十五”數字地震觀測網絡項目新建）和相對記錄室（2008 年“子午”工程新建）為無磁或弱磁性材料的地面建筑。

2007 年10 月23 日拉薩地磁臺安裝了地磁GM4 數字磁通門磁力儀和FHDZ-M15 地磁組合觀測系統，實現了數字化觀測，同時撤掉了57 型記錄儀、CB3 型記錄儀、Schmidt偏角磁力儀。臺站使用M15 組合觀測系統、GM4 磁通門磁力儀2 套儀器對地磁場D、H、Z、F分量進行數字連續記錄，其中，CTM-DI 地磁經緯儀觀測地磁D、I分量，G856 儀、FHDZ-M15 地磁觀測系統對F分量進行絕對觀測。

圖3 是拉薩地磁臺2009 年3 月至2010 年12 月地磁觀測HHT 頻譜圖。由圖3 可見，頻率基本穩定不變，發生頻率擾動的時刻反映了HHT 所具有的瞬時特性，大部分地磁分量的優勢頻率為0.01 Hz；地震前地磁總場F分量能量積累過程較顯著，其他分量則較弱，能量變化較分散。

圖3 拉薩地磁臺地磁觀測HHT 頻譜（a）D分量；（b）F分量；（c）H分量；（d）Z分量Fig.3 HHT spectrum of the geomagnetic station of Lasa

3 討論與結論

如何消除實際信號中的噪聲，從混有噪聲的信號中提取與地震有關的異常信息一直是地震學研究的熱點之一。傳統的數據處理方法，如傅立葉變換只能處理線性非平穩的信號，小波變換可處理非線性非平穩信號。而HHT 與傳統的信號或數據處理方法相比，具有能分析非線性非平穩信號、完全自適應性、不受Heisenberg 測不準原理制約——適合突變信號、采用求導得到瞬時頻率等特點。因此，HHT 方法具有更廣闊的應用前景。盡管該方法得到了越來越多的關注，但仍存在以下未能很好解決的問題。

（1）邊界處理問題。在HHT 方法中，邊界處理問題是EMD 過程的核心問題，在邊界常常出現數據端部的“飛冀”，這使得邊界一定區域的數據分解后精度很差。

（2）理論的完備性。在處理非平穩信號時，HHT 方法表現出了明顯的優越性，但在基礎理論方面，相關研究滯后，且缺少對其應用的研究，這包括對EMD 和Hilbert 變換的數學解釋和數學證明。

通過對不同臺站、地磁不同觀測分量震例HHT 頻譜進行分析可見，震前有一定前兆反應，尤其是4 級以上的中強度地震，大部分地震異常均發生在震前3 個月內，且各分量的優勢頻率均為0.01 Hz。但因觀測數據有限，并非震中周邊300 km 內的所有臺站HHT頻譜分析結果都反映出前兆異常，這可能還與臺站所處的構造位置有關，即沿同一斷裂帶排列的臺站才能監測到地震前兆現象。該問題有待進一步研究。