中国国家数字地震台网目前在广东广州、甘肃高台、广西灵山等17个台站安装了JCZ-1超宽频带地震计，该型号地震计频带为DC-20Hz，范围覆盖了从短周期地震波至固体潮汐^[1]，能记录到非常丰富的震动信息^[2-3]。在以往数字地震观测中，学者们往往倾向于研究震时数据，而忽略了震前数据，海量的数据虽然被仪器记录下来，却没有被使用。

在地震信号处理领域，傅里叶变换是最常用也是最直接的方法，但是地震数据属于非平稳信号，傅里叶变换只能对整个频率或时间跨度内的波形和频谱特性进行分析，当想要分析信号的某一具体时刻或时间段内的频率成分时，不能完整准确地表征信号，影响结果分析。Morlet ^[4]首次采用小波变换对地震波的局部性进行分析。在确定特定频率的振幅异常方面，谱分解是强大有效的方法，通过Morlet小波变换检测联合时频谱相位分量中的不连续性可以改进地震谱分解^[5]。高静怀等^[6]、敬少群等^[7]采用Morlet小波变换方法进行地震波谱分解，取得了较好效果。本文首先详细阐述Morlet小波变换的基本原理，然后以2008-05-12汶川地震为例，对2008-05-01~05-12 JCZ-1超宽频带地震计所记录到的连续波形数据进行时频分析，给出相应的处理结果以及讨论。

1 Morlet小波

Morlet小波变换是连续小波变换的一种，将连续的时间信号s(t)与Morlet小波φ(t, f)进行卷积，从而获得随时间变化的时频能量分布。连续小波变换对时域和频域都进行加窗处理，窗口面积固定但其形态可以改变，通过改变窗长进而改变分辨率的大小。

1.1 Morlet小波基本定义

Morlet小波φ(t, f)是一种复值调制的高斯函数，虚部是实部的Hilbert变换，两者之间相位差为90°。φ(t, f)在时域和频域都是高斯分布，标准差分别为σ_t和σ_f。某频率f₀的一般数学表达式为：

$ \varphi (t, {f_0}) = A{{\rm{e}}^{{\rm{j2}}\pi {f_0}t}}{{\rm{e}}^{ - {t^2}/2\sigma _t^2}} $

(1)

式中，A=(σ_tπ^1/2)^-1/2，σ_t=1/2πσ_f，A为归一化因子，保证小波基本身的能量为1。将式(1)母小波伸缩和平移后，即得到一个小波序列：

$ {\varphi _{a, b}}(t, {f_0}) = \frac{1}{{\sqrt a }}\varphi (\frac{{t{\rm{ - }}b}}{a}, {f_0}){\rm{, }}a, b \in R, a \ne 0 $

(2)

式中，a为伸缩因子，b为平移因子。

φ(t, f₀)的傅里叶变换为：

$ \mathit{\Phi }\left( \omega \right) = A\sqrt {2 \pi} {{\rm{e}}^{ - \frac{1}{2}{{\left( {\omega {\rm{ - }}{\omega _0}} \right)}^2}}} $

(3)

式中，ω₀=2πf₀为小波中心频率。

Morlet小波家族中，f/σ_f比值恒定，是时间分辨率和频率分辨率的权衡，不同的频率f所对应的σ_t和σ_f不同，也就是说在整个时频平面上，分辨率能随频率大小而变化，在实际应用中一般取值大于5。

1.2 Morlet小波变换的时频分析特点

首先，Morlet小波变换可以在时间-频率平面表征信号能量分布，观测具有非平稳特点的地震信号的瞬时频谱特征。能量分布取的是模或者模的平方，没有负值，叠加后不存在正负相抵的问题。不同于傅里叶变换的整体频谱，小波变换后的时频谱不仅涵盖了时间分析的内容，还能展现频率随时间变化的信息，有助于分析地震信号随时间的变化趋势。

其次，小波变换在时频分析过程中的分辨率可变，当分析低频信号时，其时间窗口增大；当分析高频信号时，其时间窗口减小。而短时傅里叶变换的时间窗口都是固定不变的，舍弃了时间信息，缺乏局域信息，不适合非平稳信号的处理。

2 基于Morlet小波变换的地震数据处理算法 2.1 数据来源

选取2008-05有连续波形记录的10个台站(广东广州、甘肃高台、广西灵山等)进行分析，格式为拷贝自中国地震台网中心的Miniseed的连续波形数据。

2.2 数据处理算法

对基于Morlet小波变换的地震数据处理的算法由Python编程实现。

1) 下采样。从2005年开始，中国地震台网中心所安装的JCZ-1地震计的采样频率都改为100 Hz，而地震信号是长周期的，每秒100个点的采样率过高，故先将原始数据进行下采样，把采样频率降为2 Hz。

2) Morlet小波参数设置。取σ_t=1，ω₀=8，则σ_f=1/2π，A=π^-1/4，式(2)变为：

$ \varphi \left( t \right) = {\pi ^{ - 1/4}}{{\rm{e}}^{{\rm{j8}}t}}{{\rm{e}}^{ - {t^2}/2}} $

(4)

当取σ_t=1时，f/σ_f=ω₀，通过改变ω₀的值即可改变时间分辨率和频率分辨率的权衡，由上文可知f/σ_f一般取值大于5，故而取ω₀=8。

3) Morlet小波变换。根据帕塞瓦尔定理和卷积定理，设f(t)的傅里叶变换为F(ω)，g(t)的傅里叶变换为G(ω)，IF表示傅里叶逆变换，则时域中2个信号的卷积等价于2个信号的傅里叶变换的乘积进行傅里叶逆变换的结果，表示为：

$ f\left( t \right)*g\left( t \right) = {\rm{IF}}\left[ {F\left( \omega \right)\cdot G\left( \omega \right)} \right] $

(5)

因此，在程序编写时，绕开时域的卷积，把时域计算转换到纯频域中进行，可以充分利用快速傅里叶变换运算速度快的优点。将输入的地震信号和Morlet母小波分别进行快速傅里叶变换，再利用式(5)，通过傅里叶反变换取其模，即可求得输入信号的时频分布。由于计算是在整个时间段上进行，平移因子b自动映射到整个时间段上，不用专门设置，计算结果只是伸缩因子a的函数，反映频率信息，只需对a离散化，即可对应一个频段的Morlet小波变换。

3 数据处理及结果 3.1 单天数据

首先，对10个台站2008-05-01~05-12的数据进行单独分析处理，每个台站都有东西、南北、垂直3个分向，本文选择垂直向的数据。将全天数据作为输入，经Morlet小波变换得到当天各频率分量能量随时间的变化图，选取有代表性的4个台站的时频图进行说明。

图 1的4个时频图中，上半部分是时频分布图，下半部分是连续波形图。可以看出，4个时域图像在03：00~04:00、11:00~12:00和13:00~14:00之间都出现了不同程度的抖动。在中国地震台网中心(CENC)以及美国地质勘探局(USGS)官网上查询可知，这3个时间段内都有地震发生。相应地，在时频分布图上，这3个时间段里出现了不同的幅值极点，13:00~14:00之间的幅值最大，集中于0.02~0.07 Hz频段内。此外，在0.2~0.4Hz频段内，全天都有持续的抖动。

图 2的4个时频图中，由于当天发生M_S8.0汶川地震，震时的频率能量幅值过大，其余时间的能量幅值相对较小，在图像中被掩盖。地震发生时，主要的能量集中在0.08~0.3 Hz和0.01~0.05 Hz这2个频段内。

对于同一天，不同的台站记录到的波形数据都显示出相同的时频特征，表明了基于Morlet小波变换的时频分析方法的准确性，说明其可以应用于震前波形处理，表征震前各频率分量的特点。

3.2 汶川地震震前8 d数据

在汶川地震时期，西太平洋上空出现威马逊强台风。根据中央气象台的报道，威马逊台风形成于2008-05-07，05-08增强为强热带风暴，05-09增强为台风，05-11风力达到最大，此后风力逐渐减弱，05-13转化成温带气旋。由中国气象台网站查询可知，威马逊的风暴中心与中国海岸的距离很远，最近时也超过1 300 km，对中国沿海地区的气象没有造成影响。故将05-05~05-12的波形数据连接起来，作为一个连续时间信号进行处理。由上文可知，0.5~1 Hz频段内没有出现明显的幅值变化，此节中把研究重点放在0.01~0.5 Hz频段。

为了更好地观察不同频段信号分量的能量幅值变化情况，对信号进行带通滤波处理后，分3个频段进行Morlet小波变换，即0.01~0.1 Hz、0.1~0.25 Hz和0.25~0.5 Hz。此外，汶川地震震时的能量幅值远大于于其他时刻，反映到图中会导致其他时刻的幅值不明显甚至不可见，故将汶川地震震时的能量幅值重置为0，同一台站不同频段的处理标准相同。部分结果见图 3~5。

在中国地震台网中心和美国地质勘探局网站查询后可知，图 3~5中出现的3条亮线代表 05-08和05-10在日本附近发生的M_S7.1、M_S4.5地震以及05-11在台湾以东海域发生的M_S5.6地震。除去这3个时刻的能量幅值极值，在图 3和图 5中时频图里没有发现明显的幅值变化趋势，带通滤波后的时域图像也没有明显异常。由此可得，0.0~0.1 Hz和0.25~0.5 Hz这2个频段没有出现震前的趋势特征。

图 4中，4个台站都在震前3 d出现能量幅值增大趋势。有研究指出，西太平洋上的台风气象对中国大陆地脉动的扰动集中在0.15~0.25 Hz^[8]。胡小刚等^[9]在此基础上将这一频带细分成0.1~0.18 Hz和0.2~0.25 Hz进行研究。结果发现，0.2~0.25 Hz的扰动与台风相关，0.1~0.18 Hz是非台风扰动。夏英杰等^[10]同样发现0.2 Hz左右频段出现异常，采用单台法进行信号定向后认为，地脉动信号是从我国东部海域传来，方向的变化与同期的威马逊台风路径吻合良好。在图 4中，0.1~0.18 Hz频段和0.2~0.25 Hz频段的能量幅值都有明显的增大趋势，由Morlet小波变换所得到的时频分布图良好地吻合了胡小刚等的研究结果，并且将0.1~0.18 Hz的非台风扰动频段变化特征清晰地表现出来。

从时域图像来看，无论是原始信号还是带通滤波后各个频段的信号分量，都不能直观地得到能量幅值趋势变化的特征，从侧面反映出通过Morlet小波变换分析时频特征比直接由时域对比更直观方便。

4 结语

JCZ-1超宽频带地震计的频带范围宽、噪声低、灵敏度高，记录了大量的连续波形数据，为地震研究提供了丰富的数据资料。Morlet小波变换结合了傅里叶变换能把信号表示成各频率成分的叠加和的形式的优点，又与快速傅里叶变换的固定时频窗形状不同，在分析过程中可以根据信号的频率改变尺度因子的大小，即改变时频窗的大小。运用小波变换处理低频分量多且非平稳的地震信号时，低频时的频率分辨率高，有利于分离各频率分量。利用时频分析图来观察信号，有利于分析各频率分量的能量随时间变化的趋势，对识别震前异常非常有益。另外，通过对真实的宽频带地震仪数据的处理，可以得到地震时的优势频段位于0.01~0.07 Hz，在0.08~0.3 Hz范围内也有能量极值出现，但弱于前者。其次，在汶川地震发生前的数十h中，0.1~0.18 Hz频段的能量幅值出现明显的递增趋势，验证了非台风扰动的存在。

Morlet小波成功分离了宽频带地震仪所记录的震前信号的频率分量，与已有结论符合良好，是一种有效的震前数据处理和特征分析方法。

致谢 感谢国家地震科学数据共享中心(http://data.earthquake.cn/)提供的数据。