地磁HHT頻譜分析在地震預(yù)測(cè)中的應(yīng)用

徐強(qiáng) 張建國(guó) 王同利

1）中國(guó)石家莊 050021 河北省地震局

2）中國(guó)河北 056001 邯鄲地震監(jiān)測(cè)中心站

3）中國(guó)北京 100080 北京市地震局

0 引言

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

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

1 HHT 方法計(jì)算原理

利用Hilbert-Huang（HHT）方法分析數(shù)據(jù)，即通過經(jīng)驗(yàn)?zāi)B(tài)分解，將信號(hào)分解成一系列本征模態(tài)函數(shù)，再通過Hilbert 變換得到三維時(shí)頻空間的Hilbert 譜（羅奇峰等，2003；石春香等，2003；段生全等，2005；葸曉宇等，2006，2007；張立等，2007；周摯等，2008，安張輝等，2010）。

對(duì)于任意1 個(gè)時(shí)間序列X(t)，其Hilbert 變換定義為

通過式（1）變換，X(t)、Y(t)可以組合成1 個(gè)復(fù)數(shù)信號(hào)Z(t)

為了使瞬時(shí)頻率具有意義，進(jìn)行Hilbert 變換的信號(hào)序列必須是單組分的，即時(shí)間與頻率間必須是一一對(duì)應(yīng)的關(guān)系。經(jīng)過經(jīng)驗(yàn)?zāi)B(tài)分解所得到的本征模態(tài)函數(shù)可滿足該要求，對(duì)X(t)的n階本征模態(tài)函數(shù)分別進(jìn)行Hilbert 變換后，原始時(shí)間序列可以表示為

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

邊際譜通常代表某固定頻率在整個(gè)時(shí)間區(qū)域上所有振幅的總和。

2 震例分析

2.1 汶川MS 8.0 地震

成都地磁臺(tái)始建于1971 年，1997 年地磁相對(duì)觀測(cè)儀器正式觀測(cè)，2000 年該臺(tái)進(jìn)行了數(shù)字化改造，M15 觀測(cè)儀器于2007 年開始運(yùn)行，觀測(cè)數(shù)據(jù)準(zhǔn)確可靠，臺(tái)站周邊環(huán)境良好，干擾較少。成都地磁臺(tái)距汶川地震震中約100 km。

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

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

圖1 成都地磁臺(tái)地磁觀測(cè)HHT 頻譜（a）D分量；（b）經(jīng)驗(yàn)?zāi)B(tài)分解圖；（c）干擾剔除示意圖；（d）H分量；（e）Z分量Fig.1 HHT spectrum of the geomagnetic station of Chengdu

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

2.2 汶川MS 8.0 地震余震

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

2007 年10 月23 日拉薩地磁臺(tái)安裝了地磁GM4 數(shù)字磁通門磁力儀和FHDZ-M15 地磁組合觀測(cè)系統(tǒng)，實(shí)現(xiàn)了數(shù)字化觀測(cè)，同時(shí)撤掉了57 型記錄儀、CB3 型記錄儀、Schmidt偏角磁力儀。臺(tái)站使用M15 組合觀測(cè)系統(tǒng)、GM4 磁通門磁力儀2 套儀器對(duì)地磁場(chǎng)D、H、Z、F分量進(jìn)行數(shù)字連續(xù)記錄，其中，CTM-DI 地磁經(jīng)緯儀觀測(cè)地磁D、I分量，G856 儀、FHDZ-M15 地磁觀測(cè)系統(tǒng)對(duì)F分量進(jìn)行絕對(duì)觀測(cè)。

圖3 是拉薩地磁臺(tái)2009 年3 月至2010 年12 月地磁觀測(cè)HHT 頻譜圖。由圖3 可見，頻率基本穩(wěn)定不變，發(fā)生頻率擾動(dòng)的時(shí)刻反映了HHT 所具有的瞬時(shí)特性，大部分地磁分量的優(yōu)勢(shì)頻率為0.01 Hz；地震前地磁總場(chǎng)F分量能量積累過程較顯著，其他分量則較弱，能量變化較分散。

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

3 討論與結(jié)論

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

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

（2）理論的完備性。在處理非平穩(wěn)信號(hào)時(shí)，HHT 方法表現(xiàn)出了明顯的優(yōu)越性，但在基礎(chǔ)理論方面，相關(guān)研究滯后，且缺少對(duì)其應(yīng)用的研究，這包括對(duì)EMD 和Hilbert 變換的數(shù)學(xué)解釋和數(shù)學(xué)證明。

通過對(duì)不同臺(tái)站、地磁不同觀測(cè)分量震例HHT 頻譜進(jìn)行分析可見，震前有一定前兆反應(yīng)，尤其是4 級(jí)以上的中強(qiáng)度地震，大部分地震異常均發(fā)生在震前3 個(gè)月內(nèi)，且各分量的優(yōu)勢(shì)頻率均為0.01 Hz。但因觀測(cè)數(shù)據(jù)有限，并非震中周邊300 km 內(nèi)的所有臺(tái)站HHT頻譜分析結(jié)果都反映出前兆異常，這可能還與臺(tái)站所處的構(gòu)造位置有關(guān)，即沿同一斷裂帶排列的臺(tái)站才能監(jiān)測(cè)到地震前兆現(xiàn)象。該問題有待進(jìn)一步研究。