小波分解-STFT方法在地形變觀測(cè)數(shù)據(jù)中的應(yīng)用*

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

1)中國(guó)地震局地震研究所，武漢 430071 2)地殼運(yùn)動(dòng)與地球觀測(cè)實(shí)驗(yàn)室，武漢 430071 3)中國(guó)地震臺(tái)網(wǎng)中心，北京100086

小波分解-STFT方法在地形變觀測(cè)數(shù)據(jù)中的應(yīng)用*

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

1)中國(guó)地震局地震研究所，武漢 430071 2)地殼運(yùn)動(dòng)與地球觀測(cè)實(shí)驗(yàn)室，武漢 430071 3)中國(guó)地震臺(tái)網(wǎng)中心，北京100086

把小波分解與短時(shí)傅里葉變換(STFT)相結(jié)合，實(shí)現(xiàn)了對(duì)信號(hào)高頻部分和低頻部分的精細(xì)分解，同時(shí)給出信號(hào)高頻部分的時(shí)頻譜，結(jié)果直觀明確，計(jì)算過(guò)程簡(jiǎn)單。這種方法不僅可以作為連續(xù)形變觀測(cè)數(shù)據(jù)的常規(guī)分析方法，也可用于其他連續(xù)觀測(cè)數(shù)據(jù)的分析。

小波分解;短時(shí)傅里葉變換;連續(xù)形變;高頻信號(hào);時(shí)頻譜

1 引言

以潮汐形變?yōu)橹黧w的地殼形變連續(xù)觀測(cè)的優(yōu)勢(shì)頻段是周期為數(shù)秒至數(shù)十天的中頻段地殼形變信息。該頻段信息中包含了能從理論上精確計(jì)算的傾斜、應(yīng)變、重力固體潮汐，以及地震前驅(qū)波等暫態(tài)事件的豐富信息。這些信息傳遞到地表的量級(jí)大多小于10-8，惟有潮汐形變連續(xù)觀測(cè)才能夠捕捉和分辨［1］。

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

2 小波變換與短時(shí)傅里葉變換

小波變換在實(shí)際計(jì)算中是通過(guò)Mallat算法實(shí)現(xiàn)的，它由小波濾波器H、G對(duì)信號(hào)進(jìn)行分解，算法如下。

式中，j為層數(shù)，j=1，2，…，J，J=log2N(N為信號(hào)長(zhǎng)度)。Aj為信號(hào)f(t)在第j層的近似部分(即低頻部分)，Dj為信號(hào)f(t)在第j層的細(xì)節(jié)部分。短時(shí)傅里葉變換(STFT)公式為：

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

以黃梅臺(tái)寬頻帶垂直擺東西分量2008年5月11日的秒采樣觀測(cè)數(shù)據(jù)為例予以說(shuō)明。圖1給出了小波變換與短時(shí)傅里葉變換結(jié)合求高頻信號(hào)時(shí)頻譜的過(guò)程。首先通過(guò)小波變換將信號(hào)分離成趨勢(shì)和細(xì)節(jié)兩部分，再對(duì)細(xì)節(jié)進(jìn)行短時(shí)傅里葉變換，得到高頻信號(hào)的時(shí)頻譜圖。本文在進(jìn)行STFT計(jì)算時(shí)采用Hamming窗。

圖2給出了當(dāng)窗長(zhǎng)分別為64 s、128 s和256 s的時(shí)頻圖。從圖2可以看到，窗長(zhǎng)越長(zhǎng)，對(duì)信號(hào)頻率的分辨力越高，考慮到計(jì)算機(jī)的運(yùn)算速度，一般采用128 s或256 s即可。

從圖2細(xì)節(jié)部分可以看到，當(dāng)天12點(diǎn)以前，在0.2 Hz和0.1 Hz兩個(gè)頻段上出現(xiàn)了兩個(gè)能量較強(qiáng)的信號(hào)，而12點(diǎn)之后，這兩個(gè)信號(hào)逐漸靠近，并于18時(shí)之后重合到0.15 Hz的頻段上。

圖1 小波變換與短時(shí)傅里葉變換結(jié)合求高頻信號(hào)時(shí)頻譜的過(guò)程Fig.1 Process of high-frequency signal spectrum calculation by use of the wavelet transform and short time Fourier transform spectrum

圖2 不同窗長(zhǎng)計(jì)算出的時(shí)頻譜Fig.2 Time-frequency spectrum calculated by using different window length

3 連續(xù)形變觀測(cè)數(shù)據(jù)處理

連續(xù)形變觀測(cè)數(shù)據(jù)中既包含有高于小時(shí)為周期的、主要由風(fēng)、氣壓、小振幅震顫等引起的微震動(dòng)信號(hào)，也包含有周期為小時(shí)到一天之間固體潮半日波和全日波，同時(shí)還包含一些一天以上低頻波動(dòng)的響應(yīng)。根據(jù)研究對(duì)象不同，用小波分解進(jìn)行濾波時(shí)所取的參數(shù)也不同。

對(duì)于秒采樣或分采樣數(shù)據(jù)進(jìn)行分析時(shí)，可將小波分解層數(shù)設(shè)為6層，將原始數(shù)據(jù)中周期為2～128 s)或2～128分鐘之間的信號(hào)提取出來(lái)，作為準(zhǔn)平穩(wěn)信號(hào)進(jìn)行短時(shí)傅里葉變換分析。對(duì)于小時(shí)采樣數(shù)據(jù)進(jìn)行日固體潮分析時(shí)，可將小波分解層數(shù)設(shè)為4層，將原始數(shù)據(jù)中周期為2～32小時(shí)之間的信號(hào)提取出來(lái)，作為準(zhǔn)平穩(wěn)信號(hào)進(jìn)行STFT分析，這樣做的目的是讓研究對(duì)象中包含需要研究的頻率。

3.1 高頻信號(hào)分析

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

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

大風(fēng)產(chǎn)生的地脈動(dòng)信號(hào)，經(jīng)常包含在分鐘采樣的連續(xù)形變觀測(cè)數(shù)據(jù)中。圖5是2009年10月15—22日北京地區(qū)大風(fēng)對(duì)延慶臺(tái)、北京臺(tái)和香山臺(tái)分鐘采樣傾斜觀測(cè)資料的影響表現(xiàn)。從其時(shí)頻譜圖可以看出，局部大風(fēng)引起的地脈動(dòng)信號(hào)頻域較寬，最短周期可達(dá)2分鐘，優(yōu)勢(shì)頻率不明顯。

圖3 臺(tái)風(fēng)浣熊在秒采樣寬頻帶垂直擺觀測(cè)數(shù)據(jù)中的時(shí)頻特征Fig.3 Time-frequency characteristics of the sampled by second broadband vertical pendulum observations caused by Typhoon Neoguri

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

3.2 固體潮分析

除了可以對(duì)秒采樣和分鐘采樣的形變觀測(cè)數(shù)據(jù)進(jìn)行時(shí)頻分析外，對(duì)整時(shí)值采樣的形變觀測(cè)數(shù)據(jù)也可以進(jìn)行時(shí)頻分析。圖6為觀測(cè)環(huán)境條件好、觀測(cè)精度高的甘肅白銀洞體應(yīng)變北南分量2010年全年的小時(shí)采樣觀測(cè)值用小波分解-STFT做出的結(jié)果。將觀測(cè)數(shù)據(jù)的細(xì)節(jié)部分通過(guò)短時(shí)傅里葉變換得到的時(shí)頻譜，并在頻率軸中標(biāo)出了固體潮能量較強(qiáng)的1/3日波、半日波和全日波的中心頻率。從時(shí)頻圖中可以看到，在觀測(cè)環(huán)境條件好的地區(qū)可以獲得信噪比強(qiáng)的固體潮信號(hào)，基本上沒(méi)有非固體潮頻段的雜波信號(hào)出現(xiàn)。

用相同的方法對(duì)觀測(cè)精度較差的昆明洞體應(yīng)變東西分量2010年全年的小時(shí)采樣觀測(cè)數(shù)據(jù)做出的結(jié)果見(jiàn)圖7，從時(shí)頻圖可以看出，該臺(tái)的觀測(cè)資料在固體潮頻段的信號(hào)頻率分布較為凌亂，存在較大的噪聲，固體潮主要波群頻率的集中性差。通過(guò)圖6和圖7的比較，很容易區(qū)分兩個(gè)臺(tái)站觀測(cè)資料的優(yōu)劣。

圖5 分鐘采樣連續(xù)觀測(cè)數(shù)據(jù)的小波分析-STFT計(jì)算結(jié)果Fig.5 Wavelet and STFT analysis of continuous observations sampled minute

圖6 白銀臺(tái)2010年洞體應(yīng)變北南分量觀測(cè)數(shù)據(jù)的時(shí)頻分析結(jié)果Fig.6 Time-frequency analysis result of north-south component of cave extensor meter at Baiyin station in 2010

圖7 昆明臺(tái)2010年洞體應(yīng)變東西分量觀測(cè)數(shù)據(jù)的時(shí)頻分析結(jié)果Fig.7 Time-frequency analysis results of east-west component of cave extensor meter at Kunming station in 2010

4 結(jié)果與討論

對(duì)不同采樣固體潮觀測(cè)數(shù)據(jù)的處理運(yùn)用小波濾波和STFT時(shí)頻分析方法，可以克服小波變換“看不清”信號(hào)的時(shí)頻譜和短時(shí)傅里葉變換“看不見(jiàn)”非平穩(wěn)信號(hào)頻譜結(jié)構(gòu)的缺點(diǎn)。將小波變換使用在信號(hào)濾波的環(huán)節(jié)，實(shí)現(xiàn)對(duì)信號(hào)高頻部分和低頻部分的精細(xì)分解，再將短時(shí)傅里葉變換使用在濾波后高頻、準(zhǔn)平穩(wěn)信號(hào)的整體時(shí)頻分析環(huán)節(jié)，計(jì)算過(guò)程簡(jiǎn)單、物理意義明確。

用該方法對(duì)秒采樣的連續(xù)形變觀測(cè)數(shù)據(jù)進(jìn)行分析，可以得到以地震波、臺(tái)風(fēng)激發(fā)的震顫波為主的高頻事件在頻域和時(shí)域上的表現(xiàn);對(duì)分采樣連續(xù)形變觀測(cè)數(shù)據(jù)進(jìn)行分析時(shí)，在時(shí)域表現(xiàn)明顯的風(fēng)擾事件，在頻域上表現(xiàn)并不集中;小時(shí)采樣的連續(xù)形變觀測(cè)數(shù)據(jù)在1/3日波、半日波和全日波段都能表現(xiàn)出明顯的優(yōu)勢(shì)頻率，在環(huán)境干擾小、臺(tái)基巖性好的臺(tái)站，固體潮優(yōu)勢(shì)頻段的信號(hào)突出，在環(huán)境干擾大、臺(tái)基巖性較差的臺(tái)站固體潮頻段的信號(hào)不突出。

致謝 感謝周碩愚、杜學(xué)彬教授的悉心指導(dǎo)及寬頻帶地震儀器研制組提供觀測(cè)數(shù)據(jù)!

1 陳德福，李正媛，陳鵬.定點(diǎn)潮汐形變觀測(cè)與GPS大地測(cè)量［J］.大地測(cè)量與地球動(dòng)力學(xué)，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 張燕，等，潮汐形變資料中地震前兆信息的識(shí)別與提取.大地測(cè)量與地球動(dòng)力學(xué)，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變換在固體潮分析中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，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 呂品姬，等.一種新的固體潮觀測(cè)數(shù)據(jù)特征量提取方法［J］.大地測(cè)量與地球動(dòng)力學(xué)，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 馬武剛，等.新型寬頻帶垂直擺傾斜儀的設(shè)計(jì)及應(yīng)用［J］.測(cè)繪信息與工程，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 胡小剛，郝曉光，薛秀秀.汶川大地震前非臺(tái)風(fēng)擾動(dòng)現(xiàn)象的研究［J］.地球物理學(xué)報(bào)，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

中國(guó)地震局科研運(yùn)行專項(xiàng)(201101014);中國(guó)地震局地震研究所所長(zhǎng)基金(IS201026010)，地震行業(yè)專項(xiàng)(200808033)

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

P207;P315.72+5

A