王 靜，賀 巍，汪偉明

（陜西省地震局 榆林綜合地震臺(tái)，陜西 西安 710068）

0 引 言

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

榆林綜合地震臺(tái)的地磁相對(duì)觀測(cè)系統(tǒng)采用GM4型磁通門(mén)磁力儀，主要是對(duì)磁偏角D、水平強(qiáng)度H和垂直強(qiáng)度Z三要素的測(cè)量。磁通門(mén)磁力儀記錄的是秒數(shù)據(jù)，地磁觀測(cè)數(shù)據(jù)易受地球磁場(chǎng)、觀測(cè)場(chǎng)地環(huán)境和觀測(cè)設(shè)備故障等因素的影響，導(dǎo)致觀測(cè)數(shù)據(jù)出現(xiàn)臺(tái)階、突跳、毛刺和缺數(shù)現(xiàn)象，而這些干擾和噪音都屬于隨機(jī)信號(hào)范疇。時(shí)頻分析法是針對(duì)時(shí)域和頻域聯(lián)合的分布信息，主要分析信號(hào)頻率隨時(shí)間變化的關(guān)系。文中通過(guò)使用短時(shí)傅里葉變換和小波變換兩種時(shí)頻分析方法對(duì)地磁數(shù)據(jù)進(jìn)行對(duì)比分析研究。

1 短時(shí)傅里葉變換變換

短時(shí)傅里葉變換（Short Time Fourier Transform，STFT）通過(guò)選擇一個(gè)時(shí)頻局部化的窗函數(shù)，通過(guò)移動(dòng)窗函數(shù)計(jì)算各個(gè)不同時(shí)刻的功率譜。窗的長(zhǎng)度決定頻譜圖的時(shí)間分辨率和頻率分辨率。窗長(zhǎng)越長(zhǎng)，截取的信號(hào)越長(zhǎng)，傅里葉變換后頻率分辨率越高，時(shí)間分辨率越差。但是，為了獲得信號(hào)的平穩(wěn)性，需要一個(gè)較短的窗函數(shù)。在短時(shí)間內(nèi)信號(hào)是平穩(wěn)的，時(shí)間分辨率越高，頻率分辨率越低。使用短時(shí)傅里葉變換不能兼顧時(shí)間分辨率和頻率分辨率，應(yīng)根據(jù)具體需求進(jìn)行選擇。

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

其中，ω(n)是窗函數(shù)，窗的長(zhǎng)度為Nω。選擇中心對(duì)稱(chēng)的滑動(dòng)窗，能量集中在m或者原點(diǎn)（m=0）附近。

2 小波變換

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

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

3 計(jì)算結(jié)果與分析

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

3.1 地磁秒數(shù)據(jù)的短時(shí)傅里葉變換結(jié)果

對(duì)地磁秒數(shù)據(jù)進(jìn)行傅里葉變換，結(jié)果如圖1所示。圖1（a）為地磁秒數(shù)據(jù)的時(shí)間域曲線圖，圖1（b）為地磁數(shù)據(jù)H分量在頻率域的Fourier變換，地磁數(shù)據(jù)集中在低頻段。對(duì)地磁秒數(shù)據(jù)進(jìn)行短時(shí)傅里葉變換，選擇長(zhǎng)度為41的hamming窗，得到圖1（c）的時(shí)頻平面分布的能量譜圖，圖1（d）則是地磁秒數(shù)據(jù)短時(shí)傅里葉變換的時(shí)頻分布的立體表示?？梢?jiàn)，短時(shí)傅里葉變換不僅能看出頻率的分布區(qū)間，還能看出頻率隨時(shí)間的變化。

圖1 原始數(shù)據(jù)及其對(duì)應(yīng)的短時(shí)Fourier變換譜圖

3.2 地磁秒數(shù)據(jù)的小波變換結(jié)果

對(duì)榆林臺(tái)地磁秒數(shù)據(jù)進(jìn)行小波變換，結(jié)果如圖2所示。其中，圖2（a）為地磁數(shù)據(jù)的時(shí)間域曲線圖，圖2（b）為地磁數(shù)據(jù)的一層小波分解的低頻信息，與原始信號(hào)的特征相似。圖2（c）為地磁數(shù)據(jù)的一層小波分解的高頻信息，具有原始信號(hào)的細(xì)節(jié)或差別，可看出信號(hào)整體受噪聲干擾，在高直干擾時(shí)段臺(tái)階波形明顯。圖2（d）是地磁數(shù)據(jù)原始經(jīng)過(guò)一層小波分解后的重構(gòu)信號(hào)，濾除了部分噪聲。

圖2 原始秒數(shù)據(jù)及其對(duì)應(yīng)的小波變換

4 結(jié) 論

本文采用短時(shí)傅里葉變換和小波變換對(duì)地磁數(shù)據(jù)進(jìn)行分析研究。使用短時(shí)傅里葉變換不僅可以知道頻率分布，還可以知道在哪個(gè)時(shí)間點(diǎn)上出現(xiàn)了什么頻率但是無(wú)法使時(shí)間分辨率和頻率分辨率最優(yōu)。小波變換可以在高頻分量采用高的時(shí)間分辨率和低的頻率分辨率，在低頻分量采用高的頻率分辨率和低的時(shí)間分辨率，這種多分辨分析方法對(duì)于地磁信號(hào)有很好的分析效果。