希尔伯特_黄变换谱及其在地震信号分析中的应用

第34卷第2期福州大学学报(自然科学版)Vol.34No.2 2006年4月Journal of Fuzhou University(Natural Science)Apr.2006文章编号:1000-2243(2006)02-0260-05希尔伯特-黄变换谱及其在地震信号分析中的应用陈子雄,吴琛,周瑞忠(福州大学土木建筑工程学院,福建福州350002)摘要:介绍了希尔伯特-黄变换(HHT)这一非线性、非平稳信号处理方法,并利用HHT处理了地震工程中常用的El Centro地震波,得到了该信号的Hilbert谱、边际谱和能量谱,提取了该信号的主要动力特性,并与该信号的Fourier分析结果进行了对比,显示出HHT这一方法的优越性.关键词:希尔伯特-黄变换;经验模态分解;固有模态函数;地震信号中图分类号:TU311.3文献标识码:AHilbert-Huang transform spectru m and its application in seismic signal analysisCHEN Zi-xiong,W U Chen,ZHOU Rui-zhong(College of Civil Engineering and Architecture,Fuzhou University,Fuzhou,Fujian350002,China) Abstract:HHT is a ne w method to deal with non-linear and non-stationary data.El Centro earth-quake wave is analyzed by HHT.Through the way,Hilbert spectrum,marginal spectrum and energyspec trum are got and dynamic property is extrac ted.The comparison between HHT spectrum and Fourierspec trum is made and the superiority of HHT is demonstrated.Keyw ords:Hilbert-Huang transform;empirical mode decomposition;intrinsic mode function;seismic signal地震信号具有短时、突变等特点,是一种典型的非平稳随机信号,必须对其进行分析与处理,才可以提取信号的主要特征.传统的Fourie r变换能够表述信号的频率特性,但不提供任何时域信息[1],而小波分析虽然在时域和频域都具有很好的局部化性质,但本质上仍是一种窗口可调的Fourier变换,在小波窗内的信号必须是平稳的,因而没有根本摆脱Fourier分析的局限[2].小波基的选择也是信号分析中的一个重要问题,另外,小波基的有限长会造成信号能量的泄漏,使信号的能量-频率-时间分布很难定量表述.Hilbert-Huang变换(HH T)的信号处理方法被认为是近年来对以Fourier变换为基础对线性和稳态谱分析的一个重大突破[2].它由经验模态分解(E mpirical Mode Decomposition,E MD)方法和Hilbert变换(H T)两部分组成,其核心是E MD分解.该方法采用了固有模态函数(Intrinsic Mode Function,I MF)概念以及将任意信号分解为I MF组成的思想,即E MD法,使得瞬时频率具有实际的物理意义[3].它不受Fourier分析的局限,可依据数据本身的时间尺度特征进行模态分解,分解过程中保留了数据本身的特性,再对各I MF分量进行Hilbert变换,得到信号能量在时间尺度上的分布规律,实现地震动力特性的提取.1Hilbert-Huang变换1.1经验模态分解和固有模态函数经验模态分解(EMD)的目的是通过对非线性非平稳信号的分解获得一系列表征信号特征时间尺度的固有模态函数(I MF),使得各个I MF是窄带信号,可以进行Hilbert分析.首先设定两个条件:¹整个时间序列的极大极小值数目与过零点数目相等或最多相差一个;º时间序列的任意点上,由极大值确收稿日期:2005-07-27作者简介:陈子雄(1981-),男,硕士研究生;通讯联系人:周瑞忠,教授.基金项目:教育部博士点专项科研基金资助项目(20040386004)定的包络与由极小值确定的包络的均值始终为零.能满足以上两个条件的信号称为I MF 信号,用Hilbert 变换求I MF 信号的瞬时参数,其结果是准确的.求取I MF 的具体步骤如下:1)找出原始序列x (t)的各个局部极大值.这里,为更好地保留原序列的特性,局部极大值定义为时间序列中的某个时刻的值,其前一时刻的值不比它大,后一时刻的值也不比它大.然后用三阶样条函数进行插值,得到原序列的上包络序列值x max (t).同理,可以得到下包络序列值x min (t ).2)对每个时刻的x ma x (t)和x min (t)取平均,得到瞬时平均值m(t):m (t)=[x max (t)+x min (t)]/2(1) 3)用原序列x (t)减去瞬时平均值m(t),得到一个去掉低频的新序列h(t ):h(t)=x (t)-m(t)(2)此时,如果h(t)满足I MF 的两个条件,那它就是固有模态函数.否则,把h (t)当作原序列,重复以上步骤,直至满足I MF 的定义,得到第一个固有模态函数,并以c 1(t)记之.一般来说,c 1(t)代表了原始序列中的高频部分.然后,用原序列减去c 1(t),得到剩余值序列r 1(t):r 1(t)=x (t)-c 1(t)(3)至此,提取第1个固有模态函数的过程全部完成.然后,把r 1(t)作为一个新的原序列,按照以上步骤,依次提取第2,第3,,,,直至第n 个固有模态函数c n (t).之后,r n (t)变成一个单调序列,再也没有内在模函数能被提取出来.如果把分解后的各分量合并起来,就得到原序列x (t):x (t )=E nj =1c j (t)+r n (t )(4)1.2 Hilbert 变换和Hilbert 谱通过E MD 分解得到的I MF 分量非常适合作Hilbert 变换,进而求出瞬时频率,得到HHT 谱.Hilbert 变换是一种线性变换,如果输入信号是平稳的,那么输出信号也应该是平稳的;Hilbert 变换强调局部属性,这避免了Fourier 变换时为拟合原序列而产生的许多多余的、事实上并不存在的高、低频成分.对于任一固有模态函数c(t),其Hilbert 变换^c (t)定义为[4]:^c (t )=1P P V Qc(S )t -Sd S (5)式中:P V 代表柯西主值,则对于c(t)的解析信号z (t)为:z (t)=c(t)+i^c (t )=a(t )ei H (t )(6)式中:a(t)和H (t)分别为信号x (t)的瞬时振幅和瞬时相位,按下式计算:a (t)=c 2(t)+^c 2(t)(7)H (t)=arctan (^c (t)/c(t))(8)由瞬时相位可得到信号的瞬时频率:X (t)=d H (t)/d t (9)可见,由Hilbert 变换得到的振幅和频率都是时间的函数,如果把振幅显示在频率-时间平面上,就可以得到Hilbert 幅值谱H (X ,t),简称Hilbert 谱,记作:H (X ,t)=Re E nj=1a j (t)e i Q X j(t)d t(10)H (X ,t)精确地描述了信号的幅值随时间和频率的变化规律.将H (X ,t)对时间积分,就得到Hilbert边际谱h (X ):h(X )=Q TH (X ,t)d t (11)边际谱提供了对每个频率的振幅量测,表达了在整个时间长度内振幅的累积.将振幅的平方对频率积#261#第2期陈子雄,等:希尔伯特-黄变换谱及其在地震信号分析中的应用分,可定义Hilbert 瞬时能量I E (t),如下式所示:I E (t)=QXH 2(X ,t)d X (12)Hilbert 瞬时能量提供了信号能量随时间的变化情况.将振幅的平方对时间积分,可得到Hilbert 能量谱E S (X ),如下式所示:E S (X )=Q TH 2(X ,t)d t (13)Hilbert 能量谱提供了对于每个频率的能量量测,表达了每个频率在整个时间长度内所累积的能量.2 EMD 前的准备工作2.1 端点效应的抑制在运用EMD 方法对非线性资料进行分解时,必须进行端点抑制.否则,要么会因在端点处弃值而严重影响资料的完整性;要么会因在端点处的发散而使运算溢出;如果抑制的不好,又会因/污染0度过大而使分解严重失真.本文采用了文献[5]介绍的在端点处依据端点附近数据变换插入样条插值起点的方法.在每一次I MF 分解前,先判断待分解系列起始段的变化趋势,若为上升,则将第一个局部极大值作为上包络样条插值的起点,而直接把端点作为下包络样条插值的起点.若为下降,直接把端点作为上包络样条插值的起点,将第一个局部极小值作为下包络样条插值的起点,信号另一端按同样方法处理.该法虽具有一定强制性,也缺乏明确的理论依据,但确能有效抑制端点效应.2.2 确定固有模态函数的判据要分解出固有模态函数,必须确定I MF 的判据,否则,无限重复的结果只可能得到定常振幅的调频波,这就失去了应有的物理意义.习惯上,用前后h(t)的标准差S 的大小作为I MF 的判据:S =ETt=0h k-1(t)-h k (t)2ETt=0h k-1(t)2(14)一般说来,S 的值越小,所得的固有模态函数的线性和稳定性就越好,能够分解出的固有模态函数的个数就越多.实践表明,当S 的值介于0.2~0.3时,既能保证固有模态函数的线性和稳定性,又能使所得的固有模态函数具有相应的物理意义[6].3 HHT 在地震信号分析中的应用作为Hilbert-Huang 变换在地震信号分析中的一个简单实例,对地震工程中最常用的一条地震记录El Centro 地震波(图1)的南北向分量进行分析.El Centro 波经E MD 分解,形成8个I MF 分量和一个残量,如图2所示.将分解的I FM 分量实行Hilbert 变换,得到地震波的时间-频率-振幅分布图(图3),#262#福州大学学报(自然科学版)第34卷利用公式(11)、(12)、(13)分别得到边际谱(图4)、瞬时能量谱(图5)、Hilbert 能量谱(图6),为比较还作了地震信号的Fourier 频谱(图7)、功率谱密度(图8).表1、表2分别列出了信号能量集中区的一些能量值.表1 信号能量集中区各时段所对应的能量Tab.1 Energy of the every period of centralized district of signal energyt /s 0~11~22~33~44~55~66~77~88~99~10能量412802990900325290074387012072006382807491294053482610584550t /s 10~1111~1212~1313~1414~1515~1616~1717~1818~1919~20能量22562013448001919202142802321901368201524001444108677299118t /s 20~2121~2222~2323~2424~2525~2626~2727~2828~2929~30能量174330954289803460215257140460500294810551894921530926#263#第2期陈子雄,等:希尔伯特-黄变换谱及其在地震信号分析中的应用表2 信号能量集中区各频段所对应的能量Tab.2 Energy of the every frequency band of centralized district of signal energyf /Hz 0~11~22~33~44~55~66~77~88~99~10能量78036004567400121920040473027728025035011326015446023282368814 分析与讨论1)由图2可知,E MD 首先将地震波中的高频分量分离出来,然后,所分解的IMF 分量的频率依次降低,这说明EMD 分解的效率是非常高的.同时,由于它不依赖如Fourier 、小波那样所需的基函数,自适应性好,因此更为有效地反映信号的内部特征,避免了信号能量的扩散与泄漏.2)由图2可见,E MD 方法清晰地说明了地震波中无周期性的正弦分量,如果有的话,必然满足I MF 条件分离出来,事实上并没有.这一事实证明:Fourier 分析对地震波这种非平稳信号的分析是不适用的,对非平稳信号的解释是不合理的.3)从所分解出来的各I MF 分量可以看到,前几个IMF 分量体现了原始信号的主要特征,该频率是地震信号的主要贡献频率.4)由图3可见,Hilbert 谱是随着频率分布而瞬时变化的,图中显示了该信号的大部分能量集中在10Hz 以下,30s 以内,这在边际谱、瞬时能量谱和Hilbert 能量谱中也得到体现.这在以往的信号处理方法中是不能得到的.5)由图7的Fouier 频谱和图4的边际谱对比可见,Fourier 谱给出的结果往往会低估低频部分的幅值,而高频部分又往往高于实际的幅值,这应该引起地震工程界的重视.6)图5所示的瞬时能量谱表达了信号的能量随时间的变化过程,其峰值高达155860,而在大部分能量集中的前30s 内的波能平均值只有9676,比值超过16倍,如此巨大的瞬时峰值能量势必对结构安全产生不利影响.由图6也可见,信号的能量基本上集中在10Hz 以下,这与信号的功率谱的表现情况基本一致.7)对图5所示的瞬时能量谱每时刻能量求和,可得其结果为1.4787@107;同理,对图6所示的Hilbert 能量谱的求和,其结果仍然相同,相对误差不到十万分之一,其中的误差是由于计算截断引起的,这也证实了Paraseval 公式的正确性.参考文献:[1]大崎顺彦(日).地震动的谱分析入门[M].吕敏申,谢礼立,译.北京:地震出版社,1980.[2]Huang N E,Shen Z,Long S R,et al.The empirical mode decomposi tion and Hilbert spectrum for nonlinear and non-stationaryti me series analysis[J].Proc R Soc Lond,1998,A454:903-995.[3]钟佑明,秦树人,汤宝平.一种振动信号新变换法的研究[J].振动工程学报,2002,15(2):233-238.[4]郑君里,应启珩,杨为理.信号与系统[M].北京:高等教育出版社,2000.[5]杜寿昌.基于Hilbert-Huang 变换的心音信号时频分析研究[D].昆明:云南大学,2003.[6]熊学军,郭炳火,胡筱敏,等.E MD 方法和Hilbert 谱分析法的应用与探讨[J].黄渤海海洋,2002,20(2):12-21.#264#福州大学学报(自然科学版)第34卷。