基于改进的希尔伯特-黄变换的地震强震动特征提取

郭亚茹,丰继华,杨黎薇,钟玉盛

(1.云南省地震局,昆明 650224;2.云南民族大学,昆明 650031)

0 引言

近半个世纪以来,伴随新出现时频分析方法被不断成功地应用于地球物理、环境科学和地质学等领域,现代时频分析理论得到了长足发展。在信号频谱研究和分析过程中,使用最多的是快速傅里叶变换(Fast Fourier Transform,FFT)。根据信号处理过程,可将以傅里叶理论为基础的时频分析方法归纳为三类[1]：第一类首先对信号加窗,其次进行傅里叶变换,前提条件为假设信号局部平稳,典型代表为STFT和Gabor变换等[2-3];第二类是直接更换傅里叶变换的基函数,典型代表为小波变换等[4-5];第三类是从信号中获得双线性函数,然后执行傅里叶变换,但存在交叉项的问题,如Winger-Ville分布等。

傅里叶分析的理论是将信号分解为正弦信号的加权和,当信号仅由几个有限的正弦信号组成时,使用傅里叶分析是可行的,但如果信号没有规律,使用傅里叶分解时需要大量的虚假的正弦信号去拼凑。此外,对于非线性分布的信号,始终存在来自交叉项的干扰。再者,受到不确定性原理的限制,传统以傅里叶理论为基础的时频分析方法在时间和频率分辨率并不理想,所以需要一种能自适应准确反映频率随时间变化规律的新方法。1998年,美国华裔科学家HUANG博士将经验模态分解(Empirical Mode Decomposition,EMD)算法引入到希尔伯特变换中,提出了希尔伯特-黄变换(Hilbert-Huang Transform,简称“HHT”)。该方法被认为是近二个世纪应用数学领域中最重大的发现,一经提出,便受到了研究者们的普遍关注。

HHT旨在服务于非平稳和非线性信号分析,如在地震波分析、脑电波分析研究以及机械故障检测等。在地震研究领域内,希尔伯特-黄变换的独特性在于分析信号过程中不需要基函数[6-7],其基本原理是将原始信号自适应分解为多个本征模态函数(IMF),再将IMF进行Hilbert变换得出具有物理意义的物理参量。2001年HUANG[8]利用HHT对1999年集集地震中Tcu129地震台站的加速度记录进行分析,与短期傅里叶变换和Morlet小波分析相比,得到的希尔伯特谱能准确反映地震动能量的时频分布;2003年ZHANG等[9]使用 EMD对地震波传播进行研究;2009年刘强等[10]通过HHT变换计算结构地震反应的能量;2012年WANG等[11]运用EMD方法对地震信号进行时频分析,证实了该算法的有效性;2016年梁岳[12]将改进HHT运用于储层预测性能中。

本文选取了部分云南漾濞Ms6.4地震中不同震中距的强震动台站所获取的三分量加速度记录数据,基于改进后的希尔伯特黄算法对其进行时频分析,以不同角度客观分析强震动的尺度和频域特征。通过对实际地震动记录进行特征提取后,获取的希尔伯特时频谱信号可清晰反映出强震能量随时频变化分布特征的能量谱,同时实现了时频局部化的目的。与梁宏等[13]利用模态混叠问题及集合经验模式分解方法(EEMD)分析九寨沟M7.0地震加速度记录相比,本文改进后的模态混叠问题及集合经验模式分解方法,弥补了其幅值系数不确定性的不足,其对地震强震动特征的提取过程可为类似破坏性浅源地震的强震动特征分析提供更有效的信息参数,为现代信号应用分析于强震动特征提供另一种新思路。

1 希尔伯特-黄变换的改进

1.1 引用数据来源

根据中国地震台网正式测定,北京时间2021年5月21日21时48分,云南省大理州漾濞县发生Ms6.4级地震,震中位于东经99.87°,北纬25.67°漾濞县苍山西镇境内,震源深度为8 km[14]。地震共计造成了3人遇难,34人受伤,部分房屋破裂甚至倒塌,同时云南多地州有强烈震感。此次地震是云南继Ms6.5鲁甸地震和Ms6.6景谷地震后近10 a以来发生的又一次破坏性浅源地震。震中位于滇西北地区维西-乔后断裂带附近,震源机制结果显示：此次地震是一次走滑型破裂事件,与区域构造特征一致。从地震序列的类型上可以看出：此次漾濞地震属于“前震-主震-余震型”地震[15-18]。

云南地震局强震动台网获取强震动数据25组75条强震动数据(三分向记录)。获取强震动数据的25个强震动台站,土层台有23个,基岩台有2个。按照国家强震动观测常规数据的处理流程对25个台站的记录进行数据处理[19-21]。为了阐明提取记录加速度特性效果,本文根据峰值加速度以及震中等因素,分别选取最具代表性的4个台站：震中距为8.0 km的漾濞台(YBX)、震中距为31.1 km的月溪井台(DLY)、震中距为66.8 km的太和台(BTH)、震中距为129.6 km的施甸台(SDX)。其中：太和台(BTH)为基岩台。图1为该地震中获取加速度记录的强震动台站位置分布,具体记录情况详见表1。

图1 云南强震动台站位置分布图Fig. 1 Location distribution map of strong motion stations in Yunnan

表1 云南漾濞Ms6.4强震动记录Table 1 Ms6.4 strong ground motion records of Yangbi,Yunnan

1.2 希尔伯特-黄变换简要概述

希尔伯特-黄变换(HHT)主要包含两个部分：经验模态分解(EMD)与希尔伯特变换(Hilbert transform)[22]。其原理为：先利用EMD方法分解原始信号得出本征模态函数,再将本征模态函数进行希尔伯特变换得出具有实际物理意义的参量。

1) 经验模态分解(EMD)基本原理

EMD是希尔伯特-黄变换的核心步骤,可将原始混合信号分解为一系列本征模态函数(IMF)组件。EMD不同于快速傅里叶变换(FFT),其以原有信号为基础直接进行分解,不需要基函数,适用于任意形式的信号,能够应用于多种领域。EMD的分解过程一般基于以下条件进行：

(a)用于分解的数据信号是由几种不同的本征模态函数(IMF)构成的,至少含有两个以上的极值点。

(b)IMF同时具有线性与非线性的特点。

(c)如果数据信号没有极值点只有拐点,可以利用二次或者多次求导取得极值,再进行积分。

EMD分解法是先通过寻找原始信号的极大极小值,分别取得极大值包络mmax与极小值包络mmin,再求其包络均值m(t)[22]。设原始信号为x(t),h(t)为准模态函数,在第一次分解中的第一次筛分为：

m1(t)=(mmax(t)+mmax(t))/2

(1)

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

(2)

如果第一次筛分结果h1(t)不满足IMF定义,便开始第二次筛分,此时h1(t)作为被筛分信号重复以上过程,直到满足黄锷博士等提出的停止准则,即可跳出筛分循环,完成第一次分解,得到第一阶本征模态函数c1(t),采用的停止准则为：

(3)

当Sd满足0.2≤Sd≤0.3,即可跳出第一次分解,获得第一阶本征模态函数c1(t)：

c1(t)=h1k(t)

(4)

此时令：r1(t)=x(t)-c1(t)作为第二次分解的原始信号进行筛分循环：

r2(t)=r1(t)-c2(t),……,rn(t)=rn-1(t)-cn(t)

当rn为一个常数或者一个单调递增(递减)的函数再或者只含有一个极点的函数时结束分解,此时rn为剩余量,

(5)

在IMF生成过程中,EMD在每个时间段提取频率最高的分量。因此,每个组成部分都不能间歇,否则,将发生模态混叠现象。

2) 希尔伯特变换(Hilbert transform)基本概念

假设本征模态函数x(t)的希尔伯特变换为ĉi(t),那么：

(6)

式中：P表示柯西(Cauchy)主值。根据这个定义,ci(t)和ĉi(t)形成一个复共轭对,可以得到一个解析信号zi(t)：

zi(t)=ci(t)+jĉi(t)=ai(t)·ejθi(t)

(7)

式中：

(8)

式中：ai(t)为反映ci(t)能量随时间变化的瞬时幅度,θi(t)为ci(t)的瞬时相位。而相位求导即可获得瞬时频率,如公式(9)所示：

(9)

原始数据x(t)(不包含余量rn(t))可表示为：

(10)

希尔伯特振幅谱H(ω,t)可表示为：

(11)

由希尔伯特振幅谱H(ω,t)进而获得：边际谱h(ω)和瞬时能量密度级IE：

(12)

(13)

3) 模态混叠问题及集合经验模式分解(EEMD)原理

HHT最明显的问题是由过包络、欠包络、模式混合以及干扰信号频域特征提取的噪声导致的模态混叠。集合经验模式分解(Ensemble Empirical Mode Decomposition,EEMD)由HUANG等针对模态混叠现象提出的,其原理为：利用相位-幅度耦合方法来量化高频带和低频带之间的相互作用,分析独立成份消除模态混合,改善每个本征模态分量的正交性[23]。

设x(t)为分析信号。EEMD的原理过程如下：

(a)信号xi(t)为：

xi(t)=x(t)+βω(i)(t)

(14)

式中：β为增加的白噪声的方差,ω(i)(i=1,2,……)i是EMD传导的数量,表示均值为0,方差为1的白噪声。

(c)计算每个最终的IMF：

(15)

4) 改进后的模态混叠问题及集合经验模式分解(EEMD)原理及方法

模态混叠问题及集合经验模式分解方法(EEMD)具有幅值系数不确定性,EEMD中EMD执行的总次数也是添加噪声的次数为N,添加白噪声的幅值系数为ε。WU和HUANG曾在论文中讨论：幅度的选择会直接影响分解效果,如果增加的噪声幅度太小,其可能不会引入EMD所依赖的极端变化,如果增加的噪声幅度太大,则分解结果可能不是原始数据,并且可能出现冗余IMF。根据算法提出的研究人员的建议与实验经验,建议在实验过程中选用N=100,ε取值范围为0.01～0.5较为合适。但是ε取值的具体情况应根据实验中的具体情况而定,并没有遵循某种特定的规则。

由此尝试在两个强震记录的加速度信号中加入不同幅值系数(信噪比)的白噪声,通过更改所加白噪声的标准差观察相关系数的变化,进而得出最佳标准差下最能消除模态混叠现象的白噪声的分布。

在本文中,经过对比分析发现对于连续信号,正态分布信号,具有线性关系的信号,皮尔逊相关系数作为评判指标最为合适,当然用斯皮尔曼相关系数也可以,但是就效率而言没有皮尔逊相关系数高。所以最终我们采用了皮尔逊相关系数,即ρX,Y定义为：

(16)

式中：cov(X,Y)是X和Y的协方差,σX为X的标准差,σY为Y的标准差,μX、μY分别表示X和Y的数学期望,所以相关系数可以表示为：

(17)

由图2加入噪声后的相关系数曲线可以直观得到对于漾濞台(YBX)、月溪井台(DLY)、太和台(BTH)和施甸台(SDX)各个分量的最佳系数,选取结果即表2所列。为便于叙述将添加特定幅值系数白噪声的EMD定义为AEEMD。为更为直观比较EMD、EEMD与AEEMD分解效果,以太和台(BTH)为例,对太和台(BTH)EW向量进行分解,分别计算每种方法分解得到的各个IMF的方差σ与贡献率mσ,以此作为模态混叠现象消除程度的评判标准。

图2 加入噪声后的相关系数曲线Fig.2 Correlationcoefficientcurveafteraddingnoise表2 噪声系数选取Table2 Noisefigureselection台站代码方向最佳系数YBXEW0.30NS0.50UD0.50DLYEW0.43NS0.45UD0.49BTHEW0.46NS0.46UD0.38SDXEW0.38NS0.50UD0.11

第k层本征模态函数ck(t)的方差σk为：

σk=E(ck(t)2)-(E(ck(t)))2

(18)

则第k层IMF的贡献率为：

(19)

式中：σk与mσk越大说明所得IMF分量与原始信号最为相近,分解效果越好,反之效果则越差。分解效果如图3、表3所示。

表3 三种分解方法中IMF的方差(BTHEW)Table 3 IMFs’ variance of three decomposition methods (BTHEW)

图3 分解谱(BTHEW)Fig. 3 Decomposition spectrum(BTHEW)

图3(a)为EMD分解结果,信号分为10层,图中Data为原始信号,C1-C10为10层IMF,R10为分解10层之后留下来的残余信息。从图中可以看到IMF1到IMF10都包含原始信号,我们无法分辨哪个层具有最原始信号的信息。

图3(b)为EEMD分解结果,添加根据算法提出的研究人员的建议与实验经验：100次幅值系数为0.2的高斯白噪声。

图3(c)为AEEMD即调整幅值系数为0.46后分解结果。

表3-表4所列为三种分解方法中IMF的方差与贡献率。对于EMD方法：表3-表4均显示IMF2的方差最大与贡献率最高,仅从理论可知IMF1应与原始信号最为相近,因此最为合理的解释为EMD方法出现了模态混合现象。对于EEMD方法：IMF1方差最大和贡献率最高,应该具有原始信号信息最多,图3(b)中IMF1与原始信号最为相似恰巧印证这一点。AEEMD方法同样显示IMF1方差最大和贡献率最高,并且AEEMD中IMF1的贡献率高于EEMD分解方法,所以AEEMD相较于EEMD更为突出的改善了模态混叠现象。

表4 三种分解方法中IMF的贡献率(BTHEW)Table 4 IMFs’ variance of three decomposition methods (BTHEW)

上述实验表明：添加的白噪声幅值系数与所分解数据有关,只有添加适当的白噪声,才可以最大程度改善模态混叠现象,提高分解的准确性。

2 特征提取

Hilbert-Huang变换(HHT)类似于傅里叶变换,两者均是将时域中的信号进行解析,以便于频域分析。在此,我们用希尔伯特-黄变换对强震动加速度进行分析：一是弥补傅里叶变换在信号时频特性提取中的不足;二是通过HHT量化提取信号能量的时频分布特征。

2.1 时域分解特征

图4表示经过改进后希尔伯特黄变换之后的漾濞台(YBX)、月溪井台(DLY)、太和台(BTH)和施甸台(SDX)各个分量的三维采样点-IMF层数-幅值谱,图形将4个台站每分量的采样点-IMF层数-幅值进行了立体展现,从图中可以看出：在任意一层本征模态函数均存在不同的频率成分,各频率成分的幅值(即其的相对含量)也不尽相同,与每个IMF分量对应的瞬时频率之间往往没有明显的界限;区别于小波分析中频带必须严格划分,并且随着震中距的逐渐增大,每个本征模态函数的振幅波峰差距逐渐减小,而基岩台加速度记录相对更为平稳。获取的希尔伯特时域谱信号可清晰反映出强震能量随时频变化分布特征的能量谱,同时实现了时频局部化的目的,为类似破坏性浅源地震的强震动特征分析提供更有效的信息参数。

图4 三维采样点-IMF层数-幅值谱Fig. 4 3D sample point-IMF layer number-amplitude spectrums

2.2 Hilbert边际谱分析

图5是经过改进希尔伯特黄变换后得到的漾濞台(YBX)、月溪井台(DLY)、太和台(BTH)和施甸台(SDX)各个分量的边际谱与FFT谱归一化后的对比图。以YBXEW为例,图5 YBXEW为漾濞台东西分量处理结果,震中距为8.0 km。对比FFT 谱和HHT边际谱可以直观发现：FFT谱在低频时会低估信号幅值,而在高频时放大信号幅值。漾濞台(YBX)其他两个分量(北南向NS、垂直向UD)、月溪井台(DLY)三分量(东西EW、南向NS、垂直向UD)、施甸台(SDX)三分量(东西EW、南向NS和垂直向UD)同样印实了这一点。

图5 Hilbert边际谱与傅里叶谱对比Fig. 5 Comparison of Hilbert Marginal Spectrum and Fourier Spectrum

对比4个站点Hilbert边际谱与傅里叶谱可以发现：随着震中距的增加,最大地震动的峰值加速度出现明显的衰减,而加速度记录的主频分量也出现相应增加。与梁宏等[13]利用EEMD分析九寨沟M7.0地震加速度记录相比,抑制模态混叠现象效果更为突出,经本文改进过的希尔伯特黄变换所得到的Hilbert边际谱与傅里叶谱对比图,更为清晰刻画出在震中距影响下FFT谱对信号低频低估高频放大现象的变化：随着震中距不断增加FFT谱与边际谱差异性逐渐变小。

2.3 Hilbert变换与谱分析

图6为经过改进希尔伯特黄变换后得到的漾濞台(YBX)、月溪井台(DLY)、太和台(BTH)和施甸台(SDX)三分量的谱分析结果二维显示,其中太和台(BTH)为基岩台,另外三个为土层台,时频谱横坐标为频率,纵坐标为HHT的瞬时能量,图像显示了能量分布特征：当震中距为50 km以内时能量集中于20 Hz左右,随着震中距的增大,最大加速度峰值出现急速衰减,能量的分布规律也显现出特定形态。联系时域分解特征与边际谱特征分析,该方法在强震记录数据处理时频分析方面展现了特有的实效性。

图6 HHT时频谱二维显示Fig. 6 2D HHT spectrums

3 结语

本文选取了云南漾濞Ms6.4地震中不同震中距的强震动台站所获取的三分量加速度记录数据,采用了基于改进后的AEEMD算法对其进行时频分析。通过分析可知：

1) 该方法可在保证信号分解稳定性的同时,抑制模态混叠现象。通过本文对强震动数据的分析拟合,不但成功将加速度数据分解为本征模态函数,还获取了希尔伯特时频谱信号。

2) 本文获得的边际谱可更加准确的反映信号的实际频率成分,更适于处理非平稳信号。

3) 本文获得的能量谱可清晰反映出强震能量随时频变化分布特征,实现了时频同时局部化的目的。

经分析对比发现：传统的FFT谱普遍存在着低频时低估信号幅值并且高频时又过度放大信号幅值的问题,本文所改进的方法能在一定程度上修正这种误差。实验结果表明：基于改进的希尔伯特-黄变换在地震强震动特征提取过程中具有显著优点,与前人利用EEMD分析地震加速度记录相比,抑制模态混叠现象效果更为突出,其得到的边际谱和时频谱可在最大程度上保留原始数据中的主要特性,可为类似破坏性浅源地震的强震动特征分析提供更有效的信息参数,更适于处理非平稳信号,为现代信号应用分析于强震动特征提供另一种新思路。