小波分解-STFT方法在地形變觀測數據中的應用*

呂品姬 趙 斌 陳志遙 張 燕 李正媛

1)中國地震局地震研究所，武漢 430071 2)地殼運動與地球觀測實驗室，武漢 430071 3)中國地震臺網中心，北京100086

小波分解-STFT方法在地形變觀測數據中的應用*

呂品姬1，2)趙 斌1，2)陳志遙1，2)張 燕1，2)李正媛3)

1)中國地震局地震研究所，武漢 430071 2)地殼運動與地球觀測實驗室，武漢 430071 3)中國地震臺網中心，北京100086

把小波分解與短時傅里葉變換(STFT)相結合，實現了對信號高頻部分和低頻部分的精細分解，同時給出信號高頻部分的時頻譜，結果直觀明確，計算過程簡單。這種方法不僅可以作為連續形變觀測數據的常規分析方法，也可用于其他連續觀測數據的分析。

小波分解;短時傅里葉變換;連續形變;高頻信號;時頻譜

1 引言

以潮汐形變為主體的地殼形變連續觀測的優勢頻段是周期為數秒至數十天的中頻段地殼形變信息。該頻段信息中包含了能從理論上精確計算的傾斜、應變、重力固體潮汐，以及地震前驅波等暫態事件的豐富信息。這些信息傳遞到地表的量級大多小于10-8，惟有潮汐形變連續觀測才能夠捕捉和分辨［1］。

分鐘采樣的數字化觀測數據增加了觀測數據中的較高頻率成分，保留了其中原來的低頻成分。數據分析方法除了基于固體潮理論的分析方法以外，還增加了基于數字信號處理的分析方法［2-4］。特別是近兩年新型寬頻帶垂直擺傾斜儀的使用，產出了秒采樣的形變觀測數據，在觀測結果中表現出了更加豐富的信息［5，6］。為了對其信息本質的認識，本文將運用小波分解與短時傅里葉變換相結合的方法，對連續形變觀測數據進行時頻分析。

2 小波變換與短時傅里葉變換

小波變換在實際計算中是通過Mallat算法實現的，它由小波濾波器H、G對信號進行分解，算法如下。

式中，j為層數，j=1，2，…，J，J=log2N(N為信號長度)。Aj為信號f(t)在第j層的近似部分(即低頻部分)，Dj為信號f(t)在第j層的細節部分。短時傅里葉變換(STFT)公式為：

其中：g(t)是給定的具有緊支集的窗函數，起著時限的作用;e-iwt起頻限的作用;S(ω，τ)則大致反映f(t)在τ時刻，頻率為ω的“信號成分”的相對含量。

以黃梅臺寬頻帶垂直擺東西分量2008年5月11日的秒采樣觀測數據為例予以說明。圖1給出了小波變換與短時傅里葉變換結合求高頻信號時頻譜的過程。首先通過小波變換將信號分離成趨勢和細節兩部分，再對細節進行短時傅里葉變換，得到高頻信號的時頻譜圖。本文在進行STFT計算時采用Hamming窗。

圖2給出了當窗長分別為64 s、128 s和256 s的時頻圖。從圖2可以看到，窗長越長，對信號頻率的分辨力越高，考慮到計算機的運算速度，一般采用128 s或256 s即可。

從圖2細節部分可以看到，當天12點以前，在0.2 Hz和0.1 Hz兩個頻段上出現了兩個能量較強的信號，而12點之后，這兩個信號逐漸靠近，并于18時之后重合到0.15 Hz的頻段上。

圖1 小波變換與短時傅里葉變換結合求高頻信號時頻譜的過程Fig.1 Process of high-frequency signal spectrum calculation by use of the wavelet transform and short time Fourier transform spectrum

圖2 不同窗長計算出的時頻譜Fig.2 Time-frequency spectrum calculated by using different window length

3 連續形變觀測數據處理

連續形變觀測數據中既包含有高于小時為周期的、主要由風、氣壓、小振幅震顫等引起的微震動信號，也包含有周期為小時到一天之間固體潮半日波和全日波，同時還包含一些一天以上低頻波動的響應。根據研究對象不同，用小波分解進行濾波時所取的參數也不同。

對于秒采樣或分采樣數據進行分析時，可將小波分解層數設為6層，將原始數據中周期為2～128 s)或2～128分鐘之間的信號提取出來，作為準平穩信號進行短時傅里葉變換分析。對于小時采樣數據進行日固體潮分析時，可將小波分解層數設為4層，將原始數據中周期為2～32小時之間的信號提取出來，作為準平穩信號進行STFT分析，這樣做的目的是讓研究對象中包含需要研究的頻率。

3.1 高頻信號分析

臺風或熱帶氣旋易在1 Hz以及更高頻率采樣的儀器中引起響應，本文運用小波變換-STFT法對黃梅臺秒采樣的寬頻帶垂直擺數據進行計算。圖3給出了用小波變換-STFT法對湖北黃梅臺寬頻帶垂直擺數據2008年第1號臺風(浣熊)影響期間的時頻分析結果。從圖中可以看出，臺風浣熊對黃梅臺寬頻帶垂直擺觀測數據在0.15～0.35 Hz頻段內的影響最大。

圖4為黃梅臺寬頻帶垂直擺傾斜儀2008年5月7—12日的時頻譜圖，從圖中可以看到：黃梅臺垂直擺在2008年5月9日受到臺風威馬遜的影響，在0.2～0.33 Hz出現了能量增強的擾動信號，5月10日15時左右垂直擺兩分量在0.1 Hz的刻度線上似乎同時出現了一個新的信號，到汶川地震發震時該信號和臺風威馬遜產生的信號重合。這個頻率為0.1 Hz的信號是不是文獻［6］所指的地震前由臺風觸發的、與地震發生有關的“非臺風信號”，還需要更多的研究結果來證明。時頻分析圖4將得觀測數據中所有微小的信號都一覽無疑，得到的結果更加豐富直觀。

大風產生的地脈動信號，經常包含在分鐘采樣的連續形變觀測數據中。圖5是2009年10月15—22日北京地區大風對延慶臺、北京臺和香山臺分鐘采樣傾斜觀測資料的影響表現。從其時頻譜圖可以看出，局部大風引起的地脈動信號頻域較寬，最短周期可達2分鐘，優勢頻率不明顯。

圖3 臺風浣熊在秒采樣寬頻帶垂直擺觀測數據中的時頻特征Fig.3 Time-frequency characteristics of the sampled by second broadband vertical pendulum observations caused by Typhoon Neoguri

圖4 汶川地震前及震時黃梅臺寬帶垂直擺傾斜儀高頻段的頻譜圖Fig.4 High frequency spectrum of vertical pendulum tilt observation before and during Wenchuan earthquake

3.2 固體潮分析

除了可以對秒采樣和分鐘采樣的形變觀測數據進行時頻分析外，對整時值采樣的形變觀測數據也可以進行時頻分析。圖6為觀測環境條件好、觀測精度高的甘肅白銀洞體應變北南分量2010年全年的小時采樣觀測值用小波分解-STFT做出的結果。將觀測數據的細節部分通過短時傅里葉變換得到的時頻譜，并在頻率軸中標出了固體潮能量較強的1/3日波、半日波和全日波的中心頻率。從時頻圖中可以看到，在觀測環境條件好的地區可以獲得信噪比強的固體潮信號，基本上沒有非固體潮頻段的雜波信號出現。

用相同的方法對觀測精度較差的昆明洞體應變東西分量2010年全年的小時采樣觀測數據做出的結果見圖7，從時頻圖可以看出，該臺的觀測資料在固體潮頻段的信號頻率分布較為凌亂，存在較大的噪聲，固體潮主要波群頻率的集中性差。通過圖6和圖7的比較，很容易區分兩個臺站觀測資料的優劣。

圖5 分鐘采樣連續觀測數據的小波分析-STFT計算結果Fig.5 Wavelet and STFT analysis of continuous observations sampled minute

圖6 白銀臺2010年洞體應變北南分量觀測數據的時頻分析結果Fig.6 Time-frequency analysis result of north-south component of cave extensor meter at Baiyin station in 2010

圖7 昆明臺2010年洞體應變東西分量觀測數據的時頻分析結果Fig.7 Time-frequency analysis results of east-west component of cave extensor meter at Kunming station in 2010

4 結果與討論

對不同采樣固體潮觀測數據的處理運用小波濾波和STFT時頻分析方法，可以克服小波變換“看不清”信號的時頻譜和短時傅里葉變換“看不見”非平穩信號頻譜結構的缺點。將小波變換使用在信號濾波的環節，實現對信號高頻部分和低頻部分的精細分解，再將短時傅里葉變換使用在濾波后高頻、準平穩信號的整體時頻分析環節，計算過程簡單、物理意義明確。

用該方法對秒采樣的連續形變觀測數據進行分析，可以得到以地震波、臺風激發的震顫波為主的高頻事件在頻域和時域上的表現;對分采樣連續形變觀測數據進行分析時，在時域表現明顯的風擾事件，在頻域上表現并不集中;小時采樣的連續形變觀測數據在1/3日波、半日波和全日波段都能表現出明顯的優勢頻率，在環境干擾小、臺基巖性好的臺站，固體潮優勢頻段的信號突出，在環境干擾大、臺基巖性較差的臺站固體潮頻段的信號不突出。

致謝 感謝周碩愚、杜學彬教授的悉心指導及寬頻帶地震儀器研制組提供觀測數據!

1 陳德福，李正媛，陳鵬.定點潮汐形變觀測與GPS大地測量［J］.大地測量與地球動力學，2003，(1)：107-113.(Chen Defu，Li Zhenyuan and Chen Peng.Observation of fixed point tidal deformation and GPS geodesy［J］.Journal of Geodesy and Geodynamics，2003，(1)：34-39)

2 張燕，等，潮汐形變資料中地震前兆信息的識別與提取.大地測量與地球動力學，2003，(4)：34-39(Zhang Yan，et al.Identification and extraction of earthquake precursor from tidal deformation data［J］.Journal of Geodesy and Geodynamics，2003，(4)：107-113)

3 陳濤，等.Hilbert-Huang變換在固體潮分析中的應用［J］.大地測量與地球動力學，2009，(4)：131-134.(Chen Tao，et al.Application of Hilbert-Huang transform in analysis of earth tide deformation［J］.Journal of Geodesy and Geodynamics，2009，(4)：131-134)

4 呂品姬，等.一種新的固體潮觀測數據特征量提取方法［J］.大地測量與地球動力學，2011，(2)：76-79.(Lü Pinji，et al.A new method to extract feature signal from earth tide observations［J］.Journal of Geodesy and Geodynamics，2011，(2)：76-79)

5 馬武剛，等.新型寬頻帶垂直擺傾斜儀的設計及應用［J］.測繪信息與工程，2010，35(5)：28-30.(Ma Wugang，et al.Design and application of new wide frequency band vertical pendulum tiltmeter［J］.Journal of Geomatics，2010，35 (5)：28-30)

6 胡小剛，郝曉光，薛秀秀.汶川大地震前非臺風擾動現象的研究［J］.地球物理學報，2010，53(12)：2 875-2 886 (Hu Xiaogang，Hao Xiaoguang and Xue Xiuxiu.The analysis of the non-typhoon-induced microseisms before the 2008 Wenchuan earthquake［J］.Chinese J Geophys，2010，53 (12)：2 875-2 886)

APPLICATION OF WAVELET-DECOMPOSITION AND STFT METHOD IN CONTINUOUS DEFORMATION OBSERVATION ANALYSIS

Lü Pinji1，2)，Zhao Bin1，2)，Chen Zhiyao1，2)，Zhang Yan1，2)and Li Zhengyuan3)

(1)Institute of Seismology，CEA，Wuhan 430071 2)Crustal Movement Laboratory，Wuhan 430071 3)China Earthquake Network Center，Beijing100045)

By combining wavelet transform with STFT，finely decomposing the part of signal of high frequency from the part of low frequency has been realized and the time-frequency spectrum of the high frequency of the signal is given in the same time.Its results are direct and clear and the computational course is simple.This method can be used not only for analyzing continuous deformation observation data as a conventional method，but also can be used for other continuous observations.

wavelet decomposition;Short Time Fourier Transform(STFT);continuous deformation observation; high frequency signal;time-frequency spetrum

1671-5942(2011)05-0136-05

2011-04-13

中國地震局科研運行專項(201101014);中國地震局地震研究所所長基金(IS201026010)，地震行業專項(200808033)

呂品姬，女，1979年生，碩士，助研，主要從事地殼形變觀測與管理、地殼形變資料分析研究.E-mail：pinjilv@163.com

P207;P315.72+5

A