地震前兆數據時頻分析

戴 勇 高立新 陳立峰 楊彥明 格根

（中國呼和浩特010010內蒙古自治區地震局）

地震前兆數據時頻分析

戴 勇 高立新 陳立峰 楊彥明 格根

（中國呼和浩特010010內蒙古自治區地震局）

采用平滑偽魏格納分布方法，相繼對烏海地震臺洞體應變、哈圖烏素地震臺體應變、三號地井水位、大甸子井水位、翁牛特地震臺地電阻率和寶昌地震臺地電阻率等地震前兆數據進行時頻分析，結果顯示，時頻分析可確定地震前兆數據包含的主要諧波成分及頻段；時頻分布能清晰顯示地震前兆數據干擾短期變化時段和頻段；時頻分析方法的選取和對地震前兆數據的預處理工作，將影響時頻結果信度。

時頻分析；平滑偽魏格納分布；前兆數據；能量密度

0 引言

隨機信號可分為平穩信號和非平穩信號，其中平穩信號滿足線性、高斯性和平穩性特征，可以運用傅里葉變換等從時域或頻域對其進行分析（張賢達等，1998）。但在自然界和工程領域，包括地震波、地震前兆數據等在內的信號，由于受到各種隨機因素的影響，通常是非平穩的，其統計量是時變函數（丁風和等，2007；戴勇等，2012；丁風和等，2015）。此時僅了解信號在時域或頻域的全局特性是遠遠不夠的，還需知道信號頻譜隨時間變化的情況。為了分析和處理非平穩信號，研究者對傅里葉分析進行推廣甚至根本性變革，提出并發展了一系列新的信號處理理論（戴勇等，2012），時頻分析即為其中一種重要手段。時頻分析的基本思想是，設計時間和頻率的聯合函數，描述信號在不同時間和頻率的能量密度和強度。時頻分析手段已經廣泛應用于通信、自動化、雷達、聲納、生物、天文、醫學、地球物理和故障診斷等領域，常用時頻分析方法包括短時傅里葉變換、Wigner-Ville分布、小波變換等（葛學哲等，2006）。

時頻分析在地學中的早期應用主要是對天然地震余震時頻特征的研究（Utsu T，1962）、古地磁變化的時頻研究（Chant I J，1992）以及人工地震信號的除噪及能量補償（高軍等，1996）等。隨著時頻理論的發展及計算機運算能力的增強，時頻分析在地學中應用的深度和廣度也在不斷拓寬（王培茂，2008），在地震勘探、地震前兆數據處理及異常提取等方面發揮著重要作用。地震前兆數據主要指由形變、流體、電磁3大傳統觀測體系及由其衍生、擴展的GNSS、熱紅外等產出的觀測數據，采用時頻分析進行處理，可以充分挖掘有用信息。例如，范瑩瑩等（2010）應用最大熵譜估計等方法處理2008年汶川8.0級地震震中周圍電磁臺觀測數據，研究地電、地磁場變化，發現青藏高原東北緣地震臺記錄的地電場等在汶川8.0級地震前具有功率譜值增大現象；戴勇等（2015）采用基于自適應最優核的時頻方法，處理2013年莫力達瓦、嫩江交界5.0級地震震中區格點OLR渦度值，顯示地震前存在能量密度增強現象。本文嘗試對形變、流體、電磁等學科中典型數據進行時頻分析，研究其頻譜結構及隨時間的變化特征，以促進時頻分析在地震前兆數據中的應用。

1 方法原理

采用希爾伯特—黃變換方法對地震前兆預處理數據進行分解和重構，在此基礎之上，利用平滑偽魏格納分布方法進行時頻分析。

1.1 希爾伯特—黃變換

希爾伯特—黃變換（Hilbert-Huang Transform，簡稱HHT）方法是美國工程院院士Norden Huang等提出的一種全新信號分析方法。希爾伯特—黃變換由經驗模態分解（Empirical Mode Decomposition，簡稱EMD）和希爾伯特譜分析（Hilbert Spectrum Analysis，簡稱HSA）兩部分組成。經驗模態分解對信號進行非線性自適應分解，得到不同固有模態函數(Intrinsic Mode Function，簡稱IMF)，是希爾伯特—黃變換的核心部分（王黎黎，2009；賈春花，2013）。

為保證IMF是單分量函數，必須滿足下列條件：①極值點和零交叉點的數目相同或至多相差一個；②固有模態函數在任意點，由局部極值定義的包絡均值為零。第2個條件很難滿足，一般取近似。在此采用Huang（2000）提出的近似條件，設SD是連續兩個分解結果的標準偏差，于是有

其中，hk-1(t)、hk(t)分別表示第k次分解前后的信號，constant∈[0.2，0.3]。

EMD的基本過程為：對于任意給定信號x(t)，將極大值點和極小值點分別進行曲線擬合，使兩條曲線間包含所有信號數據，從而得到x(t)的上、下兩條包絡線。m(t)記作包絡線的平均值，令h(t)=x(t) -m(t)，則h(t)為近似的IMF。將h(t)作為新的x(t)，重復以上操作，直到h(t)滿足IMF條件，得到第1階IMF分量，記作c1(t)，即：c1(t)=h(t)，令：r(t)=x(t) -c1(t)，將r(t)作為新的x(t)，重復操作，可依次得到第2階IMF分量c2(t)，第3階IMF分量c3(t)，…，最終得到分解式

式中r(t)稱為殘余函數。

1.2 平滑偽魏格納時頻分布

解析信號（復信號）z(t)定義為

其中，H(s(t))是s(t)的Hilbert變換。

平滑偽魏格納時頻分布（SPWVD）為

其中，g(u)、h(τ)是兩個實的偶窗函數，且g(0)=h(0)=1（葛學哲等，2006）。

2 時頻分析

對烏海地震臺洞體應變、哈圖烏素地震臺體應變、三號地井水位、大甸子井水位、翁牛特地震臺地電阻率和寶昌地震臺地電阻率等觀測質量較好的地震前兆數據，進行去突跳、去階變等預處理，通過希爾伯特—黃變換方法進行分解和重構，獲取去除趨勢變化和高頻噪聲的數據，選取聚集性高且交叉項小的平滑偽魏格納分布方法，對上述數據進行時頻分析。

2.1 烏海地震臺洞體應變

選取2015年7月1日至9月30日烏海地震臺洞體應變整點值觀測數據，采用希爾伯特—黃變換方法，對NS、EW測項進行經驗模態分解，對獲取的6個模態和一個剩余分量中的模態1、2、3、4進行重構，采用SPWVD方法，對去除趨勢后的體應變數據進行時頻處理，獲取時頻聯合域中能量密度分布，見圖1。由圖1清晰可見：①兩測項時頻結果均顯示，在歸一化頻率為0.042和0.083附近存在高能量密度分布，說明烏海臺洞體應變存在日潮和半日潮；②在2015年9月10日前后，日潮能量密度NS向顯著減弱，EW向有所增強，說明兩測向日潮形態均出現畸變。EW向同期在歸一化頻率為0—0.02低頻區間內出現能量密度增強現象，經調查由降水等因素引起。

圖1 烏海臺洞體應變整點值曲線（重構）及時頻分析結果(a) NS測項； (b) EW測項Fig.1 Hourly value reconstruction curve of cave strain at Wuhai Seismic Station and its time-frequency result

2.2 哈圖烏素地震臺體應變

選取哈圖烏素地震臺2015年7月1日至9月30日體應變整點值數據，通過采用希爾伯特—黃變換方法，進行經驗模態分解，并對分解后獲得的4個模態和一個剩余分量中的模態1、2、3進行重構，在此基礎之上，采用SPWVD方法對去除趨勢后的體應變數據進行時頻處理，獲取能量密度在時—頻聯合域中的分布，見圖2。由圖2清晰可見：①在歸一化頻率為0.042和0.083附近存在高能量密度分布，說明體應變存在周期為24 h和12 h的成分，應為日潮和半日潮；②半日潮能量密度未隨時間發生顯著變化，但日潮能量密度隨著時間有所減弱。

2.3 井水位

（1）三號地井：采用希爾伯特—黃變換方法，對三號地井2015年7月1日至9月30日水位整點值數據進行經驗模態分解，并對分解后獲得的7個模態和一個剩余分量中的模態1、2、3、4進行重構，在此基礎之上，采用SPWVD方法，對去除趨勢后的水位數據進行時頻處理，見圖3(a)。由圖3(a)可見：①在歸一化頻率為0.042和0.083附近存在高能量密度分布，說明三號地井水位存在周期為24時的日潮和12時的半日潮；②2015年9月1日—15日，能量密度在歸一化頻率為0.042附近有所減弱，而在歸一化頻率0—0.02區間有所增強，說明低頻因素影響水位日潮形態，進而發生畸變。

（2）大甸子井：采用希爾伯特—黃變換方法，對大甸子井水位2001年5月1日至2003年11月30日日值數據進行經驗模態分解，并對分解后獲得的7個模態和一個剩余分量中的模態1至模態7進行重構，在此基礎之上，采用SPWVD方法對，去除趨勢后的水位數據進行時頻處理。由圖3(b)可見，2002年12月27日—2003年3月17日大甸子井水位能量密度在歸一化頻率0—0.005區間顯著增強，該異常出現后5個月內，距大甸子井194 km處發生2003年8月16日內蒙古巴林左旗、阿魯科爾沁旗間5.9級地震。

圖2 哈圖烏素臺體應變整點值曲線（重構）及時頻結果Fig.2 Hourly value reconstruction curve of body strain at Hatuwusu Seismic Station and its time-frequency result

圖3 井水位重構曲線及時頻結果(a)三號地井水位整點值； (b)大甸子井水位日值Fig.3 Reconstruction curve of well water level and its time-frequency result

2.4 臺站地電阻率

（1）翁牛特地震臺地電阻率NS測道：采用希爾伯特—黃變換方法，對翁牛特地震臺地電阻率NS測道2008年以來日值數據進行經驗模態分解，并對分解后獲得的8個模態和一個剩余分量中的模態5、6、7、8進行重構，在此基礎之上，采用SPWVD方法，對去除高頻噪聲和低頻趨勢變化后的地電阻率數據進行時頻處理，見圖4(a)。結果顯示：①在歸一化頻率0.002 8附近存在高能量密度分布，說明翁牛特地電阻率存在周期為1年的變化；②2015年能量密度在歸一化頻率0—0.005區間顯著增強，說明該臺地電阻率2015年出現幅度增大的年變畸變。

圖4 地電阻率重構曲線及時頻結果(a) 翁牛特臺地電阻率日值； (b) 寶昌臺地電阻率整點值Fig.4 Time-frequency results of earth resistivity reconstruction curves

（2）寶昌地震臺地電阻率：采用希爾伯特—黃變換方法，對寶昌地震臺地電阻率NS、EW兩測道2015年7月1日至9月30日整點值數據進行經驗模態分解，并對分解后獲得的模態分量（NS測道分解為10個模態，EW測道分解為9個模態）和1個剩余分量中的模態2、3、4、5進行重構，在此基礎之上，采用SPWVD方法，對去除趨勢和高頻噪聲后的地電阻率數據進行時頻處理，見圖4(b)。結果顯示：①在歸一化頻率0.042附近存在高能量密度分布，說明寶昌臺地電阻率存在周期為24小時的日變化；②在2015年8—9月，兩測道能量密度在歸一化頻率0—0.03區間顯著增強，在歸一化頻率0.083附近顯著變化，說明地電阻率日變形態由于低頻因素影響發生畸變。

3 結論

（1）對烏海地震臺洞體應變、哈圖烏素地震臺體應變、三號地井水位、大甸子井水位、翁牛特地震臺地電阻率和寶昌地震臺地電阻率等地震前兆數據進行時頻分析，結果清晰顯示出數據的主要諧波成分及其所分布的頻段。烏海臺洞體應變、哈圖烏素臺體應變和三號地井水位整點值時頻分布中，能量密度較高區域主要位于周期為24小時和12小時附近，且隨時間變化能量密度增強或減小，其物理意義為上述前兆數據存在日潮汐和半日潮汐，且隨著時間變化潮汐現象增強或弱化；翁牛特臺地電阻率日值時頻分布中，能量密度較高區域主要位于周期為365天附近，說明存在年變；寶昌臺地電阻率整點值時頻分布中，能量密度較高區域主要位于周期為24小時附近，說明存在日變。

（2）時頻結果能清晰顯示地震前兆數據的異常變化時段和頻段。例如，翁牛特臺地電阻率日值時頻分布中，2015年能量密度在0—0.005頻段內異常增強，這說明其在2015年存在幅度增大的年變畸變；大甸子井水位能量密度在2002年12月27日至2003年3月17日0—0.005頻段存在增強現象，異常出現后5個月內距大甸子井194 km處發生內蒙古巴林左旗—阿魯科爾沁旗5.9級地震。

（3）在進行時頻分析時，由于計算機性能限制，高采樣率數據所能分析的時間段短，低采樣率數據所能分析的時間段長。研究人員在對地震前兆數據進行時頻分析時，需依據研究目的，從采樣率、時間段長度等角度選取合適數據進行處理。在此舉例說明，寶昌臺地電阻率既存在年變，也存在日變，若要研究年變等時頻特征，則適宜選取日值數據作為處理對象，若要研究周期更短的日變等時頻特征，則適宜選取整點值數據作為處理對象。

（4）在對地震前兆數據進行時頻分析時，預處理工作是不可或缺的。若預處理不充分，則會出現以下情況：①突跳、階變等可引起時頻分布中出現上述干擾對應的高能量密度分布；②形變等前兆數據中通常存在由觀測系統“零漂”等引起的趨勢變化，若在時頻分析前未去除，則時頻分布中趨勢變化的能量密度大，淹沒其他有用信息；③對于存在噪聲的地震前兆數據，若時頻分析前未濾除噪聲，則時頻結果中噪聲能量密度將影響有用信息能量密度分布。但若過度預處理，則可能濾除有效信息。因此，在對數據進行時頻分析前，應結合研究需求，對數據進行去突跳和階變等處理，選取小波變換、希爾伯特—黃變換、數字濾波、線性或非線性擬合等方法，將噪聲、趨勢等原因明晰且對提取有效信息有影響的變化予以剔除。

（5）由采樣定理（又稱奈奎斯特定理）可知，若采樣頻率為f，則信號包含的有效信息頻段最大頻率可達f/2，說明時頻結果歸一化頻率范圍一般為0—0.5（萬永革，2012），考慮到本文所得時頻結果中潮汐、年變和其他典型變化等有效信息對應的高能量密度并未分布在此歸一化頻率區間內，而是集中于某一小頻段內，因此在繪制時頻結果時，根據實際情況對歸一化頻率范圍進行合理選取，以達到時頻分布圖既包含有效信息對應的能量密度分布，同時突出主要變化對應的能量密度分布細節的目的。

（6）時頻分析能夠清晰、準確地反映信號頻譜結構隨時間的變化，在地震前兆數據處理方面具有廣闊的應用前景。但時頻分析也具有時間分辨率與頻率分辨率相矛盾、高聚集性和低交叉項相矛盾等缺點，在眾多線性、非線性時頻分析方法中，如何從數據自身特征、分辨率、聚集性和交叉性等角度選取適合地震前兆數據分析的方法，需進一步研究。

戴勇，高立新，尹戰軍，等.莫力達瓦、嫩江交界 5.0 級地震長波輻射異常研究[J].華南地震，2015，35（1）：103-106.

戴勇，高立新，尹戰軍，等.基于快速傅立葉方法的地震前兆振幅譜分析[J].地震地磁觀測與研究，2012，33（2）：51-58.

丁風和，韓曉雷，哈媛媛，等.承壓井含水層孔隙度與固體骨架和水的體積壓縮系數之間的關系[J].地球科學（中國地質大學學報），2015，40（7）：1248-1253.

丁風和，趙鐵鎖，尹占軍，等.大甸子井水位的氣壓系數及其震前異常[J].西北地震學報，2007，29（2）：174-176.

范瑩瑩，杜學彬，Jacques Zlotnicki，等.汶川MS8.0大震前的電磁現象[J].地球物理學報，2010，53（12）：2 887-2 898.

高軍，凌云，周興元，等.時頻域球面發散和吸收補償[J].石油地球物理勘探，1996，31（6）：856-866.

葛學哲，陳仲生.MATLAB時頻分析技術及其技術應用[M].北京：人民郵電出版社，2006.

賈春花.希爾伯特—黃變換及其在信號處理中的應用研究[J].電力學報，2013，28（2）：148-151.

萬永革.數字信號處理的MATLAB實現[M].北京：科學出版社，2012.

王黎黎.基于希爾伯特-黃變換的時頻分析算法研究[D].西安電子科技大學，2009.

王培茂.地震信號的時頻特征表示方法及應用[D].長春：吉林大學，2008.

張賢達，保錚.非平穩信號分析與處理[M].北京國防工業出版社，1998.

Utsu T.On the nature of three Alaskan aftershock sequences of 1957 and 1958[J].Bulletin of the Seismological Society of America, 1962, 52（2）: 279-297.

Chant I J, Hastie L M.Time-frequency analysis of magnet-telluric data[J].Geophysical Journal International, 1992, 111（2）: 399-413.

Time-frequency analysis application in earthquake precursor data processing

Dai Yong，Gao Lixin，Chen Lifeng，Yang Yanming and Gegen
(Earthquake Administration of Inner Mongolia Autonomous Region,Hohhot010010,China)

In this paper, time-frequency analysis of the precursor data of cave strain at Wuhai Seismic Station, body strain at Hatuwusu Seismic Station, water level at Sanhaodi well, water level at Dadianzi well, earth resistivity at Wengniute Seismic Station and earth resistivity at Baochang Seismic Station is studied by using smoothed pseudo Wigner-Ville distribution method.Timefrequency result can clearly shows time interval and frequency range of main components and short-term changes of earthquake precursor data.How to select time-frequency method and how to preprocess precursor data will directly affect the time-frequency results.

time-frequency analysis，smoothed pseudo Wigner-Ville distribution，precursor data，energy density

10.3969/j.issn.1003-3246.2016.06.017

戴勇（1981—），男，安徽巢湖人，工程師，主要從事數據處理及地震預測研究工作。

E-mail： daiyong06@mails.ucas.ac.cn

陳立峰（1984—），男，內蒙古赤峰人，工程師，主要從事形變資料分析工作。E-mail：lfchen@zju.edu.cn

中國地震局監測、預測、科研三結合課題（項目編號：160504）；震情跟蹤定向工作任務（項目編號：2016010407）

本文收到日期：2016-01-11