基于時頻分析法的地磁數據研究

王 靜，賀 巍，汪偉明

（陜西省地震局 榆林綜合地震臺，陜西 西安 710068）

0 引 言

目前，我國的地震前兆觀測主要有地下流體觀測、形變觀測和地磁觀測3種手段[1]。其中，地磁觀測是基于地球主磁場及其長期變化、巖石磁學和古地磁學、地磁場起源的研究。

榆林綜合地震臺的地磁相對觀測系統采用GM4型磁通門磁力儀，主要是對磁偏角D、水平強度H和垂直強度Z三要素的測量。磁通門磁力儀記錄的是秒數據，地磁觀測數據易受地球磁場、觀測場地環境和觀測設備故障等因素的影響，導致觀測數據出現臺階、突跳、毛刺和缺數現象，而這些干擾和噪音都屬于隨機信號范疇。時頻分析法是針對時域和頻域聯合的分布信息，主要分析信號頻率隨時間變化的關系。文中通過使用短時傅里葉變換和小波變換兩種時頻分析方法對地磁數據進行對比分析研究。

1 短時傅里葉變換變換

短時傅里葉變換（Short Time Fourier Transform，STFT）通過選擇一個時頻局部化的窗函數，通過移動窗函數計算各個不同時刻的功率譜。窗的長度決定頻譜圖的時間分辨率和頻率分辨率。窗長越長，截取的信號越長，傅里葉變換后頻率分辨率越高，時間分辨率越差。但是，為了獲得信號的平穩性，需要一個較短的窗函數。在短時間內信號是平穩的，時間分辨率越高，頻率分辨率越低。使用短時傅里葉變換不能兼顧時間分辨率和頻率分辨率，應根據具體需求進行選擇。

離散時間信號x(n)的短時傅里葉變換定義為[2]：

其中，ω(n)是窗函數，窗的長度為Nω。選擇中心對稱的滑動窗，能量集中在m或者原點（m=0）附近。

2 小波變換

受塔式算法的影響，1989年Mallat提出了著名的Mallat算法，可以快速實現信號的塔式多分辨率分析分解與重構。小波變換是原始信號與經過伸縮后的小波函數族的相關運算，通過尺度伸縮和平移等運算實現對信號進行多尺度細化分析。原始信號進行濾波后，高通濾波器輸出信號的高頻部分，而高頻成分具有信號的細節或差別。低通濾波器輸出信號的低頻部分，低頻成分蘊含著信號的近似特征。小波變換能較好地解決時間和頻率分辨率的矛盾[3]：小波變換的窗是可調時頻窗，在高頻時使用短窗口，在低頻時則用寬窗口，以不同的尺度觀察信號，以不同的分辨力分析信號。

離散小波變換（Discrete Wavelet Transform，DWT）是對連續小波變換的尺度、位移按照2的冪次進行離散化得到的，也稱之為二進制小波變換。取小波變換的尺度伸縮參數a=a0m，小波變換的平移參數b=nb0a0m，一維離散小波函數可表示為[4]：

3 計算結果與分析

分析榆林臺GM4型磁通門磁力儀GM4-1的地磁秒數據，選擇2019年01月02日榆林臺GM4-1受高壓直流干擾數據作為實驗數據，上臨線日00:31至11:09時Z分量累計變幅16.2 nT，9:35至9:39變幅6.1 nT。

3.1 地磁秒數據的短時傅里葉變換結果

對地磁秒數據進行傅里葉變換，結果如圖1所示。圖1（a）為地磁秒數據的時間域曲線圖，圖1（b）為地磁數據H分量在頻率域的Fourier變換，地磁數據集中在低頻段。對地磁秒數據進行短時傅里葉變換，選擇長度為41的hamming窗，得到圖1（c）的時頻平面分布的能量譜圖，圖1（d）則是地磁秒數據短時傅里葉變換的時頻分布的立體表示?？梢姡虝r傅里葉變換不僅能看出頻率的分布區間，還能看出頻率隨時間的變化。

圖1 原始數據及其對應的短時Fourier變換譜圖

3.2 地磁秒數據的小波變換結果

對榆林臺地磁秒數據進行小波變換，結果如圖2所示。其中，圖2（a）為地磁數據的時間域曲線圖，圖2（b）為地磁數據的一層小波分解的低頻信息，與原始信號的特征相似。圖2（c）為地磁數據的一層小波分解的高頻信息，具有原始信號的細節或差別，可看出信號整體受噪聲干擾，在高直干擾時段臺階波形明顯。圖2（d）是地磁數據原始經過一層小波分解后的重構信號，濾除了部分噪聲。

圖2 原始秒數據及其對應的小波變換

4 結 論

本文采用短時傅里葉變換和小波變換對地磁數據進行分析研究。使用短時傅里葉變換不僅可以知道頻率分布，還可以知道在哪個時間點上出現了什么頻率但是無法使時間分辨率和頻率分辨率最優。小波變換可以在高頻分量采用高的時間分辨率和低的頻率分辨率，在低頻分量采用高的頻率分辨率和低的時間分辨率，這種多分辨分析方法對于地磁信號有很好的分析效果。