HHT變換在地磁數據處理的應用

向 陽孫吉澤，賴愛京史勇軍高小其張建國

1）中國烏魯木齊830011新疆維吾爾自治區地震局

2）中國北京100086中國地震局地球物理研究所

3）中國新疆維吾爾自治區843000阿克蘇中心地震臺

4）中國河北056000邯鄲中心地震臺

HHT變換在地磁數據處理的應用

向 陽1）孫吉澤1），2）賴愛京3）史勇軍1）高小其1）張建國4）

1）中國烏魯木齊830011新疆維吾爾自治區地震局

2）中國北京100086中國地震局地球物理研究所

3）中國新疆維吾爾自治區843000阿克蘇中心地震臺

4）中國河北056000邯鄲中心地震臺

采用改進的EMD閾值進行分解，將地磁數據進行閾值濾波后重整；做頻率域分析，認為EMD閾值濾波可以有效去除高頻干擾，盡可能保留低頻信息。對磁靜日和磁擾日數據進行EMD消噪，做時頻分析，在HHT譜和邊際譜清晰呈現地磁信號的動態變化特征及對應的能量變化。

地磁；HHT變換；EMD分解；時頻分析；閾值濾波

0 引言

大量觀測表明，地震的孕育和發生伴隨著電磁場的變化（韓鵬，2009；徐文耀，2009；邢西淳等，2005）。采用不同信號處理方法分析地磁信號，提取時頻信息，發現數據中蘊含的有用信息是地震地磁學研究的重要內容之一。地震的孕育和發生是一個復雜的物理過程，相關的地震地磁物理機制目前還在探索；以往針對地磁信號的研究（張亮娥等，2007；吳利輝等，2009；謝凡，2012）證明，地磁信號中地震信息微弱，容易淹沒在地磁背景場噪聲中。因此，研究地磁數據中地震信息的提取方法具有重要意義。

Hilbert—Huang變換（HHT）是近年出現的一種自適應的非平穩、非線性信號處理方法，通過非平穩信號的線性化，建立以瞬時頻率為表征信號變化的基本量，創造新的時頻分析方法體系（李彥軍等，2010），對于在地磁資料中提取有用的地震前兆信息具有重要意義。

1 HHT變換

HHT變換（NE Huang，1998），其核心是利用經驗模態分解（EMD），依據自身的時間尺度定出源信號x（t）上的所有極值點，再利用擬合函數將極值點分別進行曲線擬合，得到上下包絡線（圖1），將兩條包絡曲線的平均值記為m，定義h=x（t）-m，將h視為新的源信號，重復操作直到h滿足給定的終止條件時篩選終止，記c1= h，r =x（t）-c1，由此得到x（t）的分解式

式中r為殘余函數，代表信號的平均趨勢。

圖1 EMD分解示意Fig.1 The schematic diagram of Empirical Mode Decomposition

對給定的源信號x（t）的n階IMF分量分別進行Hilbert變換，求得瞬時頻譜，綜合得到H（ω，t），對時間進行積分，獲得信號的邊際譜H（ω），即

2 改進的EMD閾值消噪

傳統的EMD去噪是將含噪的IMF分量在重構時排除，以達到去噪的目的，容易造成某些有價值的頻率成分被舍棄。根據小波變化中的閾值去噪，對EMD分解后的IMF進行閾值去噪，去除噪聲成分，可以最大程度保留原信號中的頻率成分（Donoho D L et al，1995）。

經過EMD分解得到一個或幾個IMF分量，按照一定條件濾除噪聲，并保證有用信號的完整性。給定信號x（t）經EMD分解后得到n個IMF，對每一層IMF（C1-Cn）選擇一個合適的閾值，對每個Ci進行截斷，即，然后進行信號重構。

根據Donoho等（1995）的理論，消噪閾值為

式中，σi為第i層IMF的噪聲水平，MADi為第i層IMF的絕對中值偏差，定義為

軟閾值方法為

地磁信號去噪時選擇軟閾值進行處理，得到的波形更接近實際，而硬閾值處理信號比較粗糙。以烏什地震臺地磁記錄的單日分鐘值數據［圖2（a）］為例，通過EMD利用信號自身極值進行分解，將一個復雜的信號分解為n個有規律的本征模函數（IMF）和一個殘余信號，反映源信號自身頻率特征；對每個IMF分量進行軟閾值去噪，進行重構，得到去噪后的信號［圖2（b）］。為驗證EMD閾值濾波的有效性，通過傅里葉變換，對比消噪前后傅里葉譜（圖3），結果發現，前后頻率無太大差別，優勢頻率集中在低頻段，說明消噪去除高頻影響，而對低頻段未造成損失，可以作為對干擾抑制效果的一種評判。

圖2 烏什臺地磁分鐘值EMD閾值消噪（a）烏什臺地磁分鐘值數據；（b）EMD閾值消噪后重構數據Fig.2 Reconstructed data by EMD threshold fltering at Wushi Seismic Station

圖3 消噪前后傅里葉譜比較Fig.3 The compare of FFT spectra with original data and de-nosing data

3 HHT變換分析

地磁信號各頻率成分是分段甚至毫無規律的，建立在穩態信號處理基礎上的傅里葉變換必須采用額外的諧波分量來彌補或均衡，采用該方法估計頻譜勢必帶來偏差，造成信號能量泄漏，無法反映各個時段上的頻譜變化特征。HHT方法能很好地刻畫地磁信號的動態變化特征，信號的每個突跳、持續時段和頻帶分布能夠清晰呈現出來。從時頻譜可以看出頻率成分及能量變化。選取磁靜日和磁擾日連續3天數據，對比分析HHT時頻譜。

對磁靜日信號消噪后，做HHT時頻譜和邊際譜分析，見圖4。從時頻譜上可以看出，通過EMD閾值消噪只是將源信號中蘊含的噪聲信號去除，剩余的日變長周期信息保留下來。時頻譜中清楚顯現信號各時段的不同頻率，從邊際譜也可以看出信號中存在的優勢頻率范圍（圖4），可用于日常地磁數據篩選。

對磁擾日信號，進行消噪處理和時頻分析，見圖5，可以看出，時頻譜能較好反映原始信號中高頻干擾時間段。

圖4 磁靜日地磁信號HHT時頻譜和邊際譜（a）消噪前后信號；（b）消噪后HHT時頻譜和邊際譜Fig.4 HHT spectrum and marginal spectrum of data of magnetically quiet day

圖5 磁擾日地磁信號的HHT時頻譜和邊際譜（a）磁擾日消噪前后信號；（b）消噪后磁擾日地磁信號的HHT時頻譜和邊際譜Fig.5 HHT spectrum and marginal spectrum of data of magnetically disturbed day

3 結論

隨著斷層蠕動和巖體變形，地下壓電巖體由于壓電效應產生一定電場或磁場，當能量足夠大時，被地磁儀記錄下來而成為地磁信號中的地震信息，這些信號往往比較微弱且頻率較低，容易淹沒在地磁背景場噪聲中。地磁數據分析是地磁學研究基礎，選擇合適的分析方法提取信號主要特征，對識別信號中的地震前兆信息意義重大。利用HHT方法，對新疆烏什地震臺地磁數據進行消噪處理和頻譜分析，得到以下結論。

（1）采用EMD分解，分解地磁數據信號，通過閾值濾波后進行重整，結果與原信號在頻率域上比較接近，說明EMD閾值濾波，消噪效果較好，可盡量保留原信號的優勢頻率成分。

（2）通過對磁靜日和磁擾日數據進行消噪處理和時頻分析，在HHT頻譜中，清楚可見地磁信號的動態變化特征，頻譜中信號的每個突跳、持續時段、能量變化和頻帶分布能夠清晰呈現出來。

隨著經濟發展，地磁臺站越來越受到場源的擾動、環境和噪聲的干擾，地磁數據中伴隨大量噪聲和干擾信號，呈現非穩定、非線性特征，利用HHT方法，在地磁數據處理中引入EMD濾波方法，可以有效提高信號利用率，增強數據處理質量。從時頻譜上對信號進行分析，尋找大震前信號頻譜的相關變化規律和異常特征，以期能將地磁信號更好地用于地震預報。

安張輝，元麗華，等.HHT方法在地電場數據處理中的應用［J］.地球物理學進展，2010，25（2）：525-532.

蔡劍華，湯井田，等.基于Hilbert—Huang變換的大地電磁信號譜估計方法［J］.石油地球物理學報，2010，5（5）：762-767.

蔡劍華.基于Hilbert—Huang變換的大地電磁信號處理方法與應用研究［D］.中南大學，2010.

曹志磊，周瓊，等.蒙城地震臺地磁Z分量日變化特征分析［J］.地震地磁觀測與研究，2008，29（5）：48-51.

韓鵬.地震地磁數據處理方法研究［D］.中國地震局地球物理研究所，2009.

李彥軍，胡祥云，等.HHT變換在地球物理中的應用現狀及前景［J］.工程地球物理學報，2010，7（5）：537-544.

湯井田，化希瑞，等. Hilbert—Huang變換與大地電磁噪聲壓制［J］.地球物理學報，2008，51（2）：603-610.

湯井田，蔡劍華，等.Hilbert—Huang變換與大地電磁信號的時頻分析［J］.中南大學學報（自然科學版），2009，10（5）：1 399-1 405.

吳利輝，滕云田，等.南京地磁臺地鐵干擾特征分析與抑制處理［J］.地震地磁觀測與研究，2009，30（6）：32-39.

謝凡.地磁觀測中干擾抑制方法的發展及展望［J］.地球物理學進展，2012，27（3）：967-976.

解用明，武連祥，等.地磁Z分量日變化特征［J］.地震地磁觀測與研究，2004，25（1）：52-56.

熊學軍，郭炳火，胡筱敏，等.EMD方法和Hilbert譜分析方法的應用與討論［J］.黃渤海海洋，2002，20（2）：12-21.

邢西淳，邵輝成，等.小波變換在地磁數據分析中的應用［J］.地震地磁觀測與研究，2005，26（2）：38-42.

徐文耀.地球電磁現象物理學［M］.合肥：中國科學技術大學出版社，2009

楊培杰，印興耀，張廣智，等.希爾伯特—黃變換地震信號時頻分析與屬性提取［J］.地球物理學進展，2007，22（5）：1 585-1 590.

張亮娥，殷志剛，等.太原郝莊5.0級地震前后地磁異常研究［J］.地震地磁觀測與研究，2007，28（6）：24-31.

D L Donoho.De-noising by soft-thresholding［J］.IEEE transaction on information theory，1995，41（3）：613-627.

D L Donoho，I MJohnstone et al.Wavelet Shrinkage: Asymptopia？［J］.Journal of the Royal Statistical Society，1995，57（2）：301-369.

NE Huang，Z Shen，SR Long et al.The empirical mode decomposition and the Hilbert Spectrum for nonlinear and

non-stationary time series analysis［J］.Royal Society of London Proceedings，1998，454（1 971）：903-995.

Application of HHT method in geomagnetic data handling

Xiang Yang1），Sun Jize1），2），Lai Aijing3），Shi Yongjun3），Gao Xiaoqi1）and Zhang Jianguo4）
1） Earthquake Administration of Xinjiang Uygur Autonomous Region，Urumqi 830011，China
2） Institute of Geophysics，China Earthquake Administration，Beijing 100081，China
3） Akesu Central Seismic Station，Xinjiang Uygur Autonomous Region 843000，China
4） Handan Central Seismic Station，Hebei Province 056000，China

Using the modified EMD threshold method，geomagnetic data is de-noised by threshold method and then reconstructed.By analysis in frequency domain，EMD threshold method is more effective in de-noise high frequency interference from the low frequency information.After denoising by EMD method，we make time frequency analysis to the data of magnetically quiet day and disturbed day.And the result of HHT spectrum and marginal spectrum displayed the change of dynamic frequency characteristics of geomagnetic data and the change of energy clearly.

geomagnetism，HHT，EMD，time frequency analysis，threshold flter

10.3969/j.issn.1003-3246.2015.05.009

向陽（1989—），女，助理工程師，畢業于喀什師范學院物理系，主要從事地下流體監測、 預報和科研工作

國家自然科學基金（41274079）、財政部行業專項（201408008）和新疆地震科學基金（201309） 聯合資助

本文收到日期：2014-07-29