Hilbert_Huang变换与地震信号的时频分析_图文

1、文章编号100124683(2005022207209收稿日期2004206208;修定日期2005203230。项目类别北京市自然科学基金项目(8041001、地震科学联合基金项目(604022、中国地震局三结合项目。第一作者简介武安绪,男,生于1967年,副研究员,研究方向为地震预报、地震波形处理与应用。Hilbert 2Huang 变换与地震信号的时频分析武安绪1,2吴培稚1兰从欣1徐平1,3林向东11北京市地震局,北京市苏州街28号1000802中国地震局地球物理研究所,北京1000813中国科学院地质与地球物理研究所,北京100029摘要本文介绍了HHT 时频分析方法及瞬时频率的概念

2、,给出了已知信号的经验模态分解及其时频分布,并对实际地震波形信号进行了HHT 时频处理与剖析。结果表明,HHT 方法能准确描述地震波形信号的非线性时变特征,是地震信号时频分析的有效工具。关键词:H ilbert 2H u ang 变换瞬时频率地震波形时频分析中图分类号P315文献标识码A0引言随着数字化地震台网的不断建设,采用新方法对高精度、高采样率地震数据的分析研究具有重要的现实意义(吴书贵等,2003。地震波形是具有时变特性(或称非稳态性质的典型信号(沈萍等,1999;刘希强等,2000,对于这类信号,不仅需要从总体上了解它的频率成分,而且还需要了解每一时刻信号中所包含的频率成分。目前对地

3、震信号进行分析的主要工具是傅里叶变换(胡广书,1997;郑治真,1998、现代谱估计(皇甫堪等,2003、G abor 变换(郑治真等,1996;科恩,1988、Wigner 2Ville 分布(沈萍等,1999;科恩,1988;郑治真等,1993、小波变换(刘希强等,2000;章珂等,1996;李宪优等,1999;Mallat ,1989;Daubechies ,1988;C oifman ,1990等。傅里叶变换和现代谱估计是建立在稳态信号处理基础之上,它仅能给出信号总体所包含的各种频率成分,不能解决何时出现何频率的问题;G abor 变换是在傅里叶变换基础上发展起来的一种时频分析方法,它

4、在时域和频域采用固定的分辨率,或对不同的频率成分,在时域上的取样步长不变,有时会出现G ibbs 现象;Wigner 2Ville 分布也是一种局部化的时频分析方法,它因交叉项的引入会导致时频平面上的伪影现象,使分辨率受到一定影响;小波理论与以往其他研究方法相比,在频谱分辨率和时间定位精度方面都得到了显著提高,但小波变换的有效性依赖于小波函数的选取,有时还会存在随着尺度增大相应正交基函数的频谱局部性变差的缺陷,使其对信号的更精细分解受到了一定的限制(秦前清等,1995;李世雄等,1995。第21卷第2期(2072152005年6月中国地震E ARTH QUAKE RESE ARCH I N C

5、HI NA V ol.21N o.2Jun.2005近年发展起来的希尔伯特-黄变换(Hilbert 2Huang T rans form ,HHT 是一种新的非线性信号处理方法(Huang et al.,1998;Y ue et al.,2001;Deng Y ongjun et al.,2001。HHT 方法将信号作经验模态分解(Em pirical M ode Decom position ,E MD ,能有效地将信号的各种频率成分以固有模态函数(Intrinsic M ode Function ,I MF 形式从时间曲线中分离出来。不同的I MF 分量是平稳信号或简单的非线性信号,具有简

6、单的非线性特征。其缓变波包特征意味着不同特征尺度波动的波幅随时间变化,因而也具有时间上的局域化特征,而I MF 则属于窄带信号,正好满足Hilbert 变换(胡广书,1997的要求。对I MF 序列进行Hilbert 变换,得到包含时间、频率、振幅的三维离散骨架谱,可提供更加清晰的局部细节时频特征。HHT 方法具有较好的客观性、内在性与自适应性,对信号的非线性反映能力较好,适合于对具有非线性和非平稳动态变化特征的地震信号的描述与刻划。本文试图在地震波形分析与研究中引入HHT 时频分析方法。首先通过对已知信号的仿真试验,探讨该方法的有效性;然后对地震波形信号进行多尺度经验模态分解和希尔伯特变换,

7、剖析地震波形在不同I MF 尺度上的变化特征,研究地震波形在时频谱上的演化特性;最后讨论HHT 方法的适用性及其对非平稳动态变化的地震波形信号的时频刻划能力。1基于HHT 方法的时频分析HHT 方法由Huang 变换和Hilbert 变换两部分组成。111H uang 变换Huang 变换,即E MD ,该方法的本质是通过特征时间尺度获得本征震荡模式,然后由本征震荡模式来分解时间序列资料。设地震波形信号时间序列为X (t ,则:(1找出X (t 的所有极大值点和极小值点,将其用三次样条函数分别拟合为原数据序列的上、下包络线。上、下包络线的均值为平均包络线m 10,将原序列减去m 10便可得到一

8、个去掉低频的新序列h 10,即h 10=X (t -m 10(1一般h 10不一定是一个平稳序列,为此需要对它重复上述过程。如果h 10的平均包络线为m 11,则去除该包络线所代表的低频成份后的序列为h 11,即h 11=h 10-m 11(2重复上述过程,经k 次循环后,使得到的平均包络m 1k 趋向于零,此时的h 1k 为第一阶I MF 序列,定义为分量c 1,它表示信号数据序列X (t 中最高频成份。(2用X (t 减去c 1,得到一个去掉高频成份的新序列r 1。再对r 1进行(1中的分解,便得到第二阶I MF 分量c 2。如此重复直到最后一个序列r n 不可再被分解时为止。这时的r n

9、 代表序列X (t 的残余项,即通常为X (t 的趋势项或均值。上述过程可以表述为r 1=x (t -c 1,r 2=r 1-c 2,r j =r j -1-c j ,r n =r n -1-c n(3由式(3可得X (t =n j =1c j (t +r n (t (4802中国地震21卷式(4表明原始数据序列X (t 可表示为I MF 分量和一个残余项之和。从上述的E MD 分解可看出,越是早分解出来的I MF 频率越高,第一个分解出来的I MF 序列代表原信号的最高频率成份。E MD 分解出的I MF 序列是多通带滤波的结果,每个I MF 序列都应是稳态的,可对其进行Hilbert 变换

10、或其它稳态方法的进一步分析处理。E MD 方法不但是Hilbert 变换的前提,也是其它方法进一步分析与处理的基础,这一点与多尺度小波变换是不同的,尤其在刻划非线性和动态资料方面。112H ilbert 变换与瞬时频率对于任意一个时间序列X (t ,都能得到它的希尔伯特变换结果Y (t ,即Y (t =1P +-X (t -d (5其中,P 是柯西主分量。通过这个变换,X (t 和Y (t 可以组成一个复数信号Z (t ,即Z (t =X (t +iY (t =a (t ei (t (6其中,a (t =X 2(t +Y 2(t ,(t =tan -1Y (t X (t 定义瞬时频率(t 为(

11、t =d (td t (7由(7式可看出,(t 是时间t 的单值函数,即某一时间对应某一频率。为了使瞬时频率具有意义,作Hilbert 变换的序列必须是单组分的,即某一时间只有一个频率,而经验模态分解后的固有模态函数序列恰好满足该要求,每个I MF 序列在每一点的频率唯一。把(5(7式所表示的变换用于每个固有模态函数序列,可表示为下式X (t =realnj =1a j (t e i j (t d t (8式(8忽略了残余项r n ,因为它不过是单调函数或常数值。同样的数据其傅立叶变换表示为X (t =realj =1a j e i j t (9由(8式和(9式的比较可看出,(8式是(9式的广

12、义傅立叶表达。(9式的频率j 和幅值a j 是常量,j 、a j 可构成二维傅立叶幅值频谱图;(8式的频率j (t 和幅值a j (t 是时间的变量,可构成时间、频率、幅值的三维时频谱图H (,t 。由三维时频谱图H (,t 可以定义边缘谱h (Huang et al.,1998为h (=T 0H (,t d t (10其中T 为序列的时间长度。把三维时频谱经过对时间的积分,便形成了只有频率和幅值的二维谱图。边缘谱从统计意义上表征了整组数据每个频率点的积累幅值分布,而傅里叶频谱的某一点频率上的幅值表示在整个信号里有一个含有此频率的三角函数组分。所以,边缘谱和傅氏谱的意义是不同的,这从理论上也证

13、明了二者的差别。9022期武安绪等:Hilbert 2Huang 变换与地震信号的时频分析2已知信号的经验模态分解与希尔伯特变换时频分析211已知信号用已知信号的E MD 分解及其由此产生的一系列固有模态函数I MF 重构和Hilbert 变换进行时频分析,是检验HHT 方法实际应用效果的一种最简单、直观、有效的方法。为此我们设计了含有两种成分的已知信号,其解析表达式为x (t =1+0.2sin (27.5t cos230t +0.5sin (215t +sin (2120t (11信号由一基频30H z 、调制频率15H z 的调频调幅成分,以及120H z 的正弦波叠加而成。按E MD

14、原理,上述信号的一阶I MF 应为120H z 的信号,二阶I MF 应为30H z 的调频调幅信号,更高阶的I MF 将是E MD 分解的残差,其值越小,表明E MD 分解效果越好。212经验模态分解进行Hilbert 变换的前提是E MD 分解,即Huang 变换,E MD 分解质量的好坏是时频谱估计的关键。图1由(a (f 等6 条曲线组成。其中,图1(a 的曲线代表(11式的已知信号。图1(b 代表一阶I MF 序列(c 1,是最先被分解出来的高频信号,对应120H z 的正弦波。图1(c 代表二阶I MF 序列(c 2,是第二次被分解出来的信号,对应30H z 的调频调幅信号。图1(

15、d 和图1(e 分别代表三、四阶I MF 序列(c 3、c 4,从图中可以看出,其变化幅度已是原始幅值的1197,接近一条直线,说明了E MD 分解与预计的结果十分吻合。图1(f 是c 1、c 2阶I MF 重构图形,它与图1(a 相比几乎相同,两者相关系数达019996,相对残差为01019%,引起残差的原因可能是三次样条函数等分析中所产生的误差。图1已知信号的时序、分解与重构图Fig 11The curves of I MFs reconstruction about origin simulation signal对于以上分解与重构,总体上可清楚地看到2种成分的分离情况和I MF 变化特

16、征,特别是调频调幅信号的频带与振幅变化特性在高频部分更加清晰和准确。012中国地震21卷213希尔伯特变换HHT 时频估计对已知数据进行E MD 分解,得到多阶I MF 序列(图1后,进行Hilbert 变换,并经频率整理形成HHT 时频谱阵图(图2 。图2已知信号的HHT 时频谱阵图Fig 12The curves of time 2frequency distribution about origin simulation signal从时频谱阵图(图2中可明显看出2个已知成份(分别为带状和齿状的能量分布随时间和频率的动态变化特征,谱值能很好地区分出能量随时间和频率的细微变化。能量的变化不

17、是连续的,而是离散的。带状和齿状图形表示能量的分布相对集中,集中程度和形成的谱值空间图像正好反应了调频调幅成分与正弦信号成分的实际变化特征。在没有能量的区域,谱值一般为0,而不会像其他方法那样,因为加窗造成能量泄露形成连续的空间时频谱值分布,并因此产生虚假信息,给解释工作带来困难。已知数据的HHT 时谱值分布图像特征说明该方法对于时间域和频率域分辨率都很高,不受测不准原理限制,能量突变的定位与检测能力较强,具有很强的非稳态动态变换的时频刻划能力。3地震波形分析应用实例311实际地震波形信号图3(a 是北京一个地震台记录到的一个原始波形记录,该地震持续了大约30秒,我们以此为例,检验HHT 方法

18、处理波形资料的有效性和刻划动态非线性的能力。312地震波形信号的经验模态分解图3中(b (m 曲线是图3(a 的E MD 分解图,按频率由高到低排列,分别给出了各频段I MF 序列( c 1c 12变化曲线。在高频部分,震相初动表现更清晰,特别是P 波初动,它为初动检测、震相识别、走时测量等提供更好的基础数据。从图中的超低频部分还可以看到线性趋势、非线性趋势变化,它在一定程度上可替代波形分析中的多种基线校正方法,如去均值、线性和非线性拟合及滤波等,具有消除波形背景趋势干扰的能力, 可通过选择一定的I MF 序列进行Hilbert 变换,获得特定的时频谱估计结果。E MD 分解实际上反应了该方法

19、的1122期武安绪等:Hilbert 2Huang 变换与地震信号的时频分析212 中 国 地 震 21 卷 自适应滤波能力 ,获得的 IMF 时序曲线大体上为单一谐波 , 是 Hilbert 变换进行时频估计的 基础 。 图3 波形信号的原始曲线 、 曲线图 IMF Fig13 The curves of IMFs reconstruction about actual digital seismic signal 分别表示 5 个变化阶段 313 基于经验模态分解的地震波形信号希尔伯特变换 HHT 时频分析 图 4 是图 3 ( b ( m 经过 Hilbert 变换得到的时频谱值图 。为

20、了能清楚地看到台站接收 到的地震波随时频演化的过程 ,采用二维平面等值线图表示能量变换特征的时频演化过程 。 图中横坐标为时间 ,纵坐标为频率 ,图中各点表示能量 。图 3 可分为 5 个变化比较明显的阶 段: 地震发生前的地脉动记录阶段 。此时地震还没有发生 , 在时频图上能量谱值接近于 0 。 P 波阶段 ,能量从无到有 ,特别对 P 波初动的突变性 ,被精确地刻划出来 。P 波阶段能 量值不是很大 。 S 波阶段 。这个阶段是一个渐变过程 ,没有 P 波初动突变现象 ,但能量点 密集 ,是地震能量释放的主要阶段 。 尾波产生和持续阶段 ( 张天中等 ,1990 ,含有一定能 量 ,相对较

21、弱 ,变化不稳 ,可能是介质散射造成的 。 地震波结束阶段 。此阶段地震能量基 本消失 ,标志着地震信号的结束 ,又回到正常地脉动记录阶段 。在这 5 个阶段中 , 均可明显 看到长周期背景趋势变化 ,与原始波形曲线中看到的长周期变化特征一致 。 从上述地震波的时频谱演化特征可以看出 ,地震波形信号的非线性动态变化特征可被 HHT 方法较好地刻划出来 ,演化过程可以揭示出以下信息 : 台站记到的地震波动过程是 非线性的 、 动态的 ,可能存在一定的混沌性 。 HHT 方法对地震波形的动态变化过程刻划 得比较清楚 ,反应了地震波运动的阶段性 ,每一阶段都有各自的频率特性 、 能量差异 ,其它方

22、法难于揭示出这些细微性变化 。这种变换特征为震源和介质的研究提供一种新的思路与方 法 。 时频谱阵图上的能量主要集中在 S 波阶段进行释放 ,这与实际地震发生情况相符 。 对突变点的检测能力较高 ,验证了 HHT 方法具有完全局部时频特性 。 在 IMF 序列曲 线上和时频谱值图像中 ,地震数据的突变点 、 持续时段和频带能量分布均能够清晰地显现 。 2期 利用此特性可以为诸如初动检测 、 震相识别 、 频散曲线计算 、 值分频 、 Q 多尺度反演等等提 供条件 。上述分析显然是初步的 ,大量的实际地震波形记录利用 HHT 方法进行深入的分析 与研究后 ,相信会得到更进一步的认识 。 4 讨论

23、 通过理论分析及已知数据与地震波形信号的 EMD 分解和 Hilbert 变换 ,可看到 HHT 方 法的高精度分解优势与非线性动态数据时频刻划能力 。但要真正有效地用好 HHT 方法 ,还 应注意以下几点 : ( 1 HHT 方法算得的结果是离散骨架谱 。在 EMD 分解基础上 , 进行 Hilbert 变换得到的 时频谱阵是离散谱 ,它与其它方法算得的时频谱不一样 ,由此会引起频谱值的解释方式可能 有所不同 。更好地 、 合理地解释 HHT 时频谱值 ,需要进一步深入研究与分析 。 ( 2 波动现象 。HHT 方法多次用到三次样条函数 , 这可能会引起过冲或欠冲现象 。如 高阶 IMF 序

24、列 ( 低频部分 可能会出现不应有的波动现象 。选择一定的 IMF 序列进行 Hilbert 变换 ,可在一定程度上避免样条函数造成的波动影响 。对于高频部分 ,结果应该是比较可靠 的。 (3 端点效应 。端点效应的解决程度直接影响 HHT 方法的应用效果 。解决得好 , 时频 5 结论 分析效果就好 ,否则结果可能不理想甚至是完全错误的 。消除端点效应的影响是 HHT 方法 一个重要研究课题 。本文利用简单的小波串方法 ( Huang et al1 ,1998 可较好地避免序列两 此小波串替代原数据端部相同数据长度的波 。 前文介绍了 HHT 方法 、 仿真试验和地震波形的时频估计与分析结果

25、 ,可得如下结论 : 端 “飞点” 现象 : 确定序列的两个极值点 , 确定极值间的点数 ; 用最后两个连续极值间 ( 极大值与极小值之间或极小值与极大值之间 的波的 115 个周期成份组成一个小波串 , 用 Fig14 The curves of time2frequency distribution about actual digital seismic signal 武安绪等 : Hilbert2Huang 变换与地震信号的时频分析 图4 波形信号的 HHT 时频谱值图 213 214 中 国 地 震 21 卷 (1 经 EMD 分解变换得到的 IMF 序列是直接从原始时序数据中分离出

26、来的 , 事先无需 确定分解阶次 ,不受人为因素影响 , 不存在机械分解 。因此 IMF 序列能更好反映原始数据 固有的物理特性 ,其分解是客观的 , 内在的和自适应的 。每阶 IMF 序列都代表了某种特定 意义的频带信息 , 给实际应用与解释工作带来了方便 , 是取得高精度时频估计的前提和基 础。 ( 2 原始信号经过 EMD 分解获得的 IMF 序列具有稳态性 , 不但可用 Hilbert 变换计算三 维时频谱值 ,而且也可用于多尺度建模 、 谱分析等工作 ,这是其他方法所不具备的优势 。 (3 实际工作中 ,可根据需要选择不同阶次的 IMF 序列 ( 不一定相邻 进行 Hilbert 变

27、换 时频估计 ,具有很大的灵活性 。 (4 将 HHT 方法用于对非平稳动态变化的地震波形信号的分析 , 能够敏感捕捉到地震 波形信号随时间和频率动态变化的不同阶段的主要特征 ,其图像清晰 ; 能够反应能量突变点 信息 ,其时频局部定位能力较强 。HHT 时频谱的大部分能量主要集中在一定的时间和频率 范围内 ,而不是整个时频空间 ,这可能真实地反应了非线性序列的本质特征 。 HHT 方法克服了其它一些方法的缺陷 , 完全取消了窗函数的作用 , 其结果不受核函数 影响与时频测不准原理限制 ,具有完全的局部时频特性 ,可准确描述地震信号的时变特征 。 HHT 方法可在初动检测 、 震相识别 、 信

28、号分频 、 频散计算 、 序列类型识别等方面具有潜在的 应用价值 。 致谢 : 本文在修改过程中得到了张天中研究员和魏富胜副研究员的帮助 。 参考文献 胡广书 ,1997 ,数字信号处理理论 、 算法与实现 ,北京 : 清华大学出版社 。 皇甫堪 、 陈建文 、 楼生强编著 ,2003 ,现代数字信号处理 ,北京 : 电子工业出版社 。 科恩 L. ( 美 著 ,白居宪译 ,1988 ,时 - 频分析 : 理论与应用 ,西安 : 西安交通大学出版社 。 李宪优 、 蒲明辉 ,1999 ,小波变换实现地震信号的滤波去噪与压缩 ,现代电子技术 , (3 :810 。 李世雄 、 刘家琦 ,1995

29、 ,小波变换和反演数学基础 ,北京 : 地质出版社 。 刘希强 、 周蕙兰 、 ,2000 ,基于小波包变换的地震数据时频分析方法 ,西北地震学报 ,22 (2 :143147 。 李 红 沈 、 萍 郑治真 ,1999 ,瞬态谱在地震与核爆识别中的应用 ,地球物理学报 ,42 (2 :233240 。 秦前清 、 杨宗凯 ,1995 ,实用小波分析 ,西安 : 西安电子科技大学出版社 。 吴书贵 、 蒋秀琴主编 ,2003 ,数字遥测地震台网建设与运行 ,北京 : 地震出版社 。 章 、 珂 刘贵忠 ,1996 ,二进小波变换方法的地震信号分时分频去噪处理 ,地球物理学报 ,39 (2 :2

30、65271 。 张天中 、 高龙生 、 张卫平 ,1990 ,滇西试验场区的 Q 值及其随时间窗的变化 ,地震学报 , (12 :1221 。 郑治真 ,1998 ,数字信号处理基础 ,北京 : 地震出版社 。 郑治真 、 张少芬 ,1993 ,瞬态谱估计理论及其应用 ,北京 : 地震出版社 。 郑治真 、 、 ,1996 ,从 Cabor 变换到小波分析 ,中国地震 ,12 (3 :237242 。 沈 萍 谢 永 Coifman R. R. , 1990 , Adapted Multiresolution Analysis , Computation , Signal Processing

31、 and Operator Theory , Proceeding of the International Congress of Mathematicians , Kyoto ,Math Soc Japan :887897. Daubechies I. ,1988 ,Orthogonal Bases of Compactly Supported Wavelets , Communication Pure And Application Mathematics ,41 :909 996. Deng Y ongjun ,Wang Wei ,Qian Chengchun et al ,2001

32、,Oundary2Processing2Technique in EMD Method and Hilbert Transform , Chinese Science Bulletin ,46 (11 :954960. 2期 Mallat S. G. ,1989 ,Multifrequency Channel Decompositions of Images and Wavelet Models , IEEE Trans ASS P ,37 (12 :20912110. Geoscience and Remote Sensing Symposium ,I2 G ARSS , IEEE 2001

33、 International , (5 :20612063. 01 1 ,2 1 Earthquake Administration of Beijing Municipality ,Beijing 100080 ,China 2 Institute of Geophysics ,China Earthquake Administration ,Beijing 100081 ,China 3 Institute of Geology and Geophysics ,Chinese Science Academy ,Beijing 100029 ,China frequency are introduced