基于Hilbert-Huang变换在九寨沟、宁强地震的研究

张 健， 乔 波， 何 杨， 王卫平

(陕西省地震局，陕西 西安 710000)

随着地震监测台网的数字化、网络化大力发展，选用新方法实现地震数据精度和采样率双高的分析是目前急需的技术。地震波形是一种随时间变化的典型的地球物理数据信号。对于这类信号的研究，与传统的频谱分析方法相比，Hilbert-Huang变换(简称HHT)不受Heisenberg测不准原理的制约，它完全摆脱了线性和平稳性的制约，非常适合用于分析非线性非平稳信号。

1 HHT的实现过程

从物理角度来看，并非所有的信号都能用瞬时频率来分析，只有当信号只包含一种振动模式，而没有叠加波的情况时才能利用瞬时频率分析。但是根据一般的瞬时频率求解方法，求得的瞬时频率既有正的也有负的。而负的频率是没有任何物理意义。为了获得有物理意义的瞬时频率，N.E.Huang等人研究并认为任何一种信号都可以由基础信号即本征模态函数IMF(Intrinsic Mode Function)组成，该本征模态函数是一种局部限制而非全局限制的函数形式，实现这一目的要用一种特别的分解方法，将信号分解为瞬时频率能够合理定义的分量形式；因此定义了一种基于信号的局部特性的函数—IMF。根据以上定义，一个本征模态函数需满足以下两个条件：

(1)信号中的极值点的数目与过零点的数目相等或者至多相差一个；

(2)信号的极大值点构成的上包络与极小值点构成的下包络关于时间轴对称。

基于这样的想法，他们提出了经验模式分解(EMD)：可以用经验模态分解方法将信号的固有模态筛选出来的这种方法具体就是个筛选过程，实现振动模式的提取。

EMD分解过程可概括如下：

①给定信号x(t)，求出所有信号中的每个区域内的极大值和极小值，为更好保留原始信号的性质，极大值被认为是时间序列中的某个时刻的值，它只需要满足大于前、后时刻的值就可以了。局部极小值的获取同理。然后，使用三次样条插值函数进行拟合，获得上包络线xmax(t)和下包络线xmin(t)；

②计算分解原始信号所需的均值m(t)=[xmax(t)+xmin(t)]/2；③用原始信号x(t)减去均值m(t)，得到第一组h(t)=x(t)-m(t)；但是，h(t)不一定就是一个IMF，如果h(t)不满足IMF两个条件，就把h(t)替换为原始信号，重复(1)～(3)的步骤，直到满足SD条件为止。SD的表达式为：

(1)

通常情况下SD<0.3。这时满足固有模态函数条件的h(t)作为一个IMF，令c1(t)=h(t)，至此第一个IMF已经成功的提取了。剩余的r(t)=x(t)-c1(t)仍然包含有用信息，因此可以把它看成新原始信号；重复上述过程，依次得到第二个c2(t)，第三个c3(t)，…，当r(t)满足函数单调性或为常数条件时，终止筛选，这样就完成了提取IMF。由此可得x(t)的表达式为：

(2)

即原始序列是由n个IMF与一个趋势项r(t)组成。

完成了经验模式分解过程就得到了所有可提取的IMF，在Hilbert变换的基础上计算瞬时频率就可以了。

用下面方式表示信号的前提是对固有模态函数进行Hilbert变换，即：

(3)

因为希尔伯特变换的振幅和频率都是和时间有关，用三维立体图形表达幅值、频率和时间的关系，就可以得到Hilbert谱H(w,t)。

如果把H(w,t)对时间积分，就得到希尔伯特边际谱h(w)：

(4)

可以用式(4)定义的希尔伯特瞬间的能量：

(5)

事实上，如果振幅的平方对时间积分，可以得到希尔伯特能量谱：

(6)

希尔伯特能量谱表达了每个频率在整个时间长度内所累积的能量。

2 HHT在数据分析中的初步运用

HHT不仅在滤波方面有出色的表现，还在数字信号处理方面有独特的方法，根据公式(4)与公式(5)可以得出，HHT还可以研究在像地震这种非平稳信号在发展过程中频率在时间上的分布以及能量随时间的具体变化情况等一系列隐藏在地震波背后的信息。

以2017年8月8日四川省北部阿坝州九寨沟县发生的Ms7.0级地震(以下简称九寨沟地震)和2018年9月12日陕西汉中市宁强县Ms5.3级地震(以下简称宁强地震)为例，研究HHT方法在地震波形中的更深一步的运用。通过收集陕西省南部地区15个地震台站(见表1)所记录到的2次地震的地震波形数据进行分析。

表1 各台震中距统计表

其中第二列为各台到九寨沟地震震中距，第三列为各台到宁强地震震中距。

我们首先以青木川台所记录到的垂直向九寨沟地震波数据为例，如图1所示的波形进行分解分析。

图1 青木川台记录的垂直向九寨沟地震波

使用EMD方法对上图的波形数据进行分解得到图2所示。

将地震波信号分解为共12道IMF分量，对以上各分量进行希尔伯特变换后可以得出每道分量的瞬时频率，如图3所示。

求出各分量的瞬时频率，根据公式(3)计算得出希尔伯特谱，如图4所示。

通过图2以及图4综合分析，我们认为地震发展过程可分为以下几个过程：

第一阶段：0～17.7 s，此时处于地震发生前的地脉动阶段。此时地震还没有发生，在图中表现为振动频率小，主要频率成分为低频地脉动与背景噪声。

第二阶段：17.7～61 s，此阶段地震波首波已经到达台站，图4中用颜色表示其振幅大小，颜色由蓝到红表示振幅由小到大，清晰的描述了各频率成分在时间轴的分布特征以及能量分布特征。由图4分析出，在第20.3 s为P波的初动时刻。此外，地震信号频率分布在0.3～8 Hz之间，属于近震，所以首先可以判断Pn震相是首先到达的。

第三阶段:第61～107 s，通过图2可以看到此阶段波形幅值与周期增大，为S波出现并发展的阶段，也是地震波能量集中释放的阶段。而且通过图4，可以分析出，在第80～100 s之间是地震波发展到顶峰阶段，振幅达到最大但其频率最低，为0.06～0.1 Hz。符合S波频率小，振幅大的特征。

第四阶段：第107～401 s之间，在此阶段，地震波主要持续阶段已经结束。从分解到的波形来看，这个阶段含有一定振幅，但幅值很小，频率时大时小，变化不稳定。根据该台的震中距以及所处的地质构造判断，产生这样的波形是由于地震波在当地的地质体中发生散射，或者是地震波的尾波造成。

第五阶段：第401 s之后，为大地回归平稳阶段。随着地震波及其产生的尾波逐渐结束， 又回到正常地脉动记录阶段。

图2 青木川台记录的地震信号分解

图3 每道IMF分量的瞬时频率

图4 九寨沟地震青木川台Hilbert谱

3 地震波频率在时间域中的分布

在对15个台站记录的波形数据进行分解的基础之上，对其每一道IMF分量根据公式(3)进行Hilbert变换，然后根据公式(4)进行计算，得出2次地震15个台全部的Hilbert谱，我们以青木川台为例，用三维图像呈现出来用于表达幅值、频率和时间的关系。如图5所示。

图5 九寨沟地震青木川台三维Hilbert谱

由图5可见，三维的Hilbert谱清晰刻画出接收到的地震波频率、振幅随时间的变化特征，振幅幅值随着时间迅速降低。具体各台对该次地震反应情况见下表2。

表2 各台九寨沟地震反应相关参数汇总表

续表2

台站名称主要振动持续时间/s最大振幅值/×106瞬时频率分布范围/Hz西乡台46111.50.0313～2.5镇巴台4754.10.0313～2.5石泉台4028.10.0313～1.9紫阳台4784.10.0313～2.3汉阴台4933.10.0313～1.8略阳台2628.80.0313～3.1留坝台2468.30.0313～2.8佛坪台3735.90.0313～2.2凤县台30110.50.0313～3.5太白台4324.30.0313～2.8眉县台3526.10.0313～2.2

同理，我们以宁强台为例，绘制宁强地震的三维Hilbert谱，见图6。

图6 宁强地震宁强台三维Hilbert谱

具体各台对该次地震的反应情况见下表3。可以看出在该次地震过程中地震主要持续时间较短，频率成份分布广。

表3 各台宁强地震反应相关参数汇总表

对比表2与表3，结合表1可以看出，表2对应的九寨沟地震震级大，震中距都在190～500 km以内，属于近震，其地震主要振动持续时间在第190～500s左右的范围内，而且瞬时频率范围在0.0313～4.2 Hz。而对于表3对应的宁强地震震级小，震中距在20～300 km以内，属于地方震。其地震主要振动持续时间大部分在第110～430 s左右的范围内，而且瞬时频率范围在0.0313～12.8 Hz。

通过对比表2和表3结合图5、图6可以得出以下结论：①震中距与瞬时频率的范围呈反比例关系，随着震中距增大，频率范围逐渐减小。②震级越大，其地震所对应的主要频率越集中。

4 地震波能量在频率域的分布

我们主要以九寨沟地震为例，经过计算，图7是15个台的能量谱图。

图7 九寨沟地震垂直与水平向能量谱

根据公式(6)可以得出，该图谱反映了各个瞬时频率成份在整个地震波持续时间内所积累的能量。由图可以看出地震波能量在随着瞬时频率的增大有一个现增加到最大值然后迅速衰减直至为零。以图7为例，在0.0313 Hz处开始出现地震波能量然后迅速增大，在0.0938 Hz处达到最大值，达到了的能量，在此处也是各个台能量达到最大。随后迅速衰减，直到瞬时频率大于3 Hz时，能量衰减至及以下的数量级，与最大值相差4个数量级。所以，通过结合表3分析我们可以得到，九寨沟地震的优势瞬时频率的范围在0.0313～2 Hz之间，即该地震主要释放的能量所对应的频率范围。

5 结论

本文通过利用HHT方法分别对九寨沟地震和宁强地震地数据进行处理分析，利用这种方法可以很好的得到以下结果：(1)HHT方法很好的解析了接收到的地震波频率特征、振幅特征、能量特征随时间发展过程，整个过程可以分五个阶段；(2)总结出两次地震震中距、幅值和持续时间、瞬时频率之间的相关关系；(3)统计出地震波优势瞬时频率范围，并得出各个频率所对应的能量大小：九寨沟地震的主要优势频率为0.0938 Hz，其累计释放了约大小的能量，同理可得到宁强地震的主要优势频率为0.219 Hz，累计释放了大小的能量。