小波分解-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