地震波形的希尔伯特黄变换特征提取与隐马尔科夫模型识别研究

周海军，李 磊

滁州学院电子与电气学院,滁州,239000

随着科学技术的发展，地震震动记录仪能检测到频率更高或更低、能量更微弱的地下震动信号，实际记录到的信号往往同时包含多种不同来源的同时发生或发生时刻接近的多个天然地震，或其他如矿塌、化学爆炸及核爆等人工震源。检测到的震动信号本身是非平稳、非线性的时变信号，分析这种信号需要了解某一频率的时间信息或者某一时间的频率信息。传统的傅里叶变换、瓦格纳-Ville分布和小波变换在处理这种震动信号时存在局限性。美国科学家提出希尔伯特黄变换(HHT)分析方法，其扩张的基础是自适应性，能通过谱分析得到瞬时频率，从而得到信号的频率时间分布信息。

1 地震和爆破波形的差异性

天然地震和人工爆破都发生在地下，以震动波的形式迅速释放能量，从时域和频域上分析波形，发现地震和爆破存在诸多差异[1]。

(1)爆破源持续的时间一般在3秒以内，是极短时间内能量膨胀释放，P波初动方向向上。天然地震源持续的时间在10秒以上，是板块位置移动，P波方向可能向上，也可能向下。

(2)天然地震发生在地底几公里，甚至几十公里深处，介质密度高，波形衰减慢，频率成分较多。人工爆破多发生在地表浅处，土质松软，波形衰减快，频率成分较少。

(3)相同介质中，低频信号比高频信号衰减慢，因此，天然地震波形的顶峰频率高于人工爆破的顶峰频率。

(4)天然地震P波和S波的振幅比值要小于人工爆破。

基于这些差异，国内外在天然地震和人工爆破的识别上进行了深入研究，提出震动波形的多种特征[2]，诸如：P波和S波的振幅比、AR模型系数、震源深度、倒谱、S波的对数谱以及自相关系数等，也提出多种方法判断震源类型，诸如：人工神经网络、支持向量机和隐马尔科夫模型(HMM)等。

2 HHT基本原理

HHT方法是一种非线性非平稳数据的分析方法，特别是对时频能量的分析，由经验模态分解(EMD)和Hilbert谱分析(HSA)两部分组成[3]。

2.1 震动信号的经验模态分解

大多数的震动数据都不是本征模函数，包含多个波动模式，因此，直接对原始信号Hilbert变换不能完全表征地震信号的频率特性，这时需要对信号进行经验模态分解获得本征模函数(IMF)。

原始震动信号为x(t)，找到信号的所有局部极值点，利用三次样条插值函数和极大值点，拟合形成地震信号的上包络线ymax(t)，利用三次样条插值函数和极小值点，拟合形成地震信号的下包络线ymin(t)，计算上下包络线的均值得到m1(t)：

(1)

用x(t)减去m1(t)得到h1(t)：

h1(t)=x(t)-m1(t)

(2)

若h1(t)还存在负的局部极大值或正的局部极小值，则需要继续按照上面过程进行筛选。这时将h1(t)当作原始地震信号，参照公式(1)和(2)计算，直到满足本征模函数条件，得到的h11(t)就是第一个IMF分量。

r1(t)=x(t)-h11(t)

(3)

将r1(t)当作新的信号，重复上面三个公式，得到多个IMF分量[4-5]，直到函数单调或者小于预设值，分解结束。图1是某地震台站记录的原始天然

图1 天然地震的原始波形

图2 天然地震的EMD分解波形

地震波形，长度是10 000个点，前期和后期有部分无效信息，以P波初至点为起始点，截至2 000个点长度的波形作为EMD分解波形。图2是2 000点长度的天然地震波形经过EMD分解之后的9个IMF分量波形。

图3是某地震台站记录的原始人工爆破波形，长度是10 000个点，前期和后期有些无效信息，以P波初至点为起始点，截至2 000个点长度的波形作为EMD分解波形。图4是2 000点长度的人工爆破波形经过EMD分解之后的9个IMF分量波形。

图3 人工爆破的原始波形

图4 人工爆破的EMD分解波形

2.2 Hilbert谱分析

(4)

其中，P为柯西主值[4]。

瞬时振幅：

(5)

瞬时相位：

(6)

由θi(t)求导可以得到瞬时频率：

(7)

图5是图2的IMF分量的Hilbert变换波形，

图5 天然地震的Hilbert变换

图6 人工爆破的Hilbert变换

图6是图4的IMF分量的Hilbert变换波形。这样地震信号x(t)表示为：

(8)

从公式(8)中得知，地震信号可利用“时间-频率-幅值”的三维图描述，这种地震信号幅值的时频特性就是Hilbert谱。

3 数据集选择和特征提取

3.1 数据集选择

地震记录仪采集到的地震数据包括垂直UD、东西EW、南北NS三个方向，其中垂直方向UD的信号更能反映地震信号的波形特征，实验选取了首都圈地震带记录的天然地震和人工爆破事件，事件波形数据格式：EVT，记录事件2002年5月到2007年12月，震级小于2.7，天然事件35个，人工爆破事件27个。110个台站记录垂直方向天然地震数据3 805份，人工爆破数据2 849份，部分台站未记录数据。再由人工选择波形较好的120个人工爆破和120个天然地震数据。为了消除震级对地震信号识别的影响，通常对数据进行归一化，将震动信号变为0到1之间的数据。

3.2 分帧与加窗

地震信号在短时间内可认为是平稳的信号，将每个短时间内采集的信号看成一帧。这样地震信号可被分为m段，处理数据时只需要对分段的信号进行处理。通常利用窗函数来截取地震信号，在长度n点的窗函数中，汉明窗的平滑低通特性较好[6]。汉明窗的定义如下式。

(9)

对地震信号加汉明窗分帧，则每帧信号表示为：

x′(n)=x(n)×w(n)

(10)

3.3 线性预测倒谱系数(LPCC)特征参数

地震信号的LPCC参数是处理后的2 000点长的地震信号经过分帧和加汉明窗，然后对信号进行Z变换，取变换后的对数，再进行逆Z变换，从得到线性预测系数变换而成。若线性预测系数的参数是ak，LPCC的参数是ck，则：

(11)

其中，k为LPCC的阶数，p为线性预测系数的阶数。

图7 特征参数提取过程

4 HMM模型

HMM是用来描述一个隐含未知参数的马尔可夫过程。利用可观察的参数确定该过程的隐含参数，然后利用这些参数作进一步分析。天然地震和人工爆破模型建立和识别采用的是连续高斯混合密度隐马尔可夫模型(CGHMM)。

4.1 震源信号中CGHMM的初始化

4.2 模型建立

4.3 震源类型识别

随机选取数据集中的20个天然地震数据和20个人工爆破数据作为训练数据，建立天然地震和人工爆破HMM模型。再从余下的天然地震数据集和人工爆破数据集中分别选取20个数据作为识别数据进行识别；HMM模型的状态个数为4，每个状态有3个高斯函数；每个帧的长度为64，帧移为50；将正确识别的个数除以总的识别数据的个数作为识别单次概率。按照上述操作，随机10次求得平均识别率见表1。

表1 LPCC参数的识别率

5 总结与分析

虽然地震信号是非平稳、非线性的随机信号，但在短时间内可看成是平稳的、线性的信号。因此，本文采用了分帧和加窗的方式来得到短时信号特性，能较好地提高天然地震和人工爆破的识别率。

国内外在利用特征识别震源信号时，一般采用多特征识别，识别率在90%以上。根据地震信号的波形特征，实验利用语言信号的识别方法识别地震信号。由于地震信号存在P波、S波和面波，有两个尖峰信号，故通过HHT处理，选择适当的维度，通过单一特征识别震源类别，识别率在85%，利用同一波形的多个分量，能将识别率提升到93.1%。

从表1可知：HMM模型识别率随着LPCC特征的维数增大而增大，达到一定程度时，识别率反而会减小。第一个IMF分量波形信息更丰富，LPCC特征参数识别率在三个分量中最高。利用三个IMF分量各自的识别结果综合判断，有两个或两个以上的认为是某种震源，系统就确定是这种震源类型。震源类型的识别实质是对震动性质和传播路径的分析判断。文中震源选择的区域是首都圈地震带，选取的区域范围较小，今后应进一步对我国及周边大型区域的较大震级的震动信号进行识别研究。