第46卷 第1期2007年 1月中山大学学报(自然科学版)ACT A S C I E NTI A RUM NAT URAL I U M UN I V ERSI T ATI S S UNY ATSE N I Vol 146 No 11Jan 1 2007经验模式分解算法的探讨和改进3郑天翔,杨力华(中山大学科学计算与计算机应用系,广东广州510275)摘 要:对经验模式分解算法中的滤波停止条件和端点延拓问题进行了研究。
在改进的E MD 算法基础上,通过对本征模函数使用“新的滤波停止条件”,获得了更好的实验分解结果,同时,由于改进的E MD 算法假定信号是无限长的,回避了B 样条插值中节点延拓的固有问题,研究了有限长度信号的端点延拓问题,给出了端点延拓算法,从而弥补了已有方法的不足,使之更具实用性。
实验表明,文中提出的算法是有效的。
关键词:经验模式分解;端点延拓;本征模函数;滤波停止条件;B 样条插值中图分类号:TP274 文献标识码:A 文章编号:052926579(2007)0120001206 H ilbert 2Huang Transfor m (简称HHT )是近年来发展起来的一种新的时间序列信号分析方法[1](以下简称H98)。
其核心是经验模式分解(E mp ir 2icalMode Decompositi on,E MD ),它把复杂的信号分解成若干个本征模式函数(I ntrinsic Mode Func 2ti on,I M F )之和。
由于E MD 是自适应的,故其分解非常有效,尤其适用于非线性和非平稳过程分析。
HHT 自1998年由N 1Huang 及其合作者提出以来,一直受到国内外学者的关注,并取得了一系列的研究成果。
Huang 所提出的E MD 是算法型的,虽然该算法在实际信号分解中十分有效,但迄今为止并没有关于该算法的收敛性结果。
实际上,人们在利用E MD 进行信号分解时,有两个方面是采取了主观的规则:其一是根据人们对零均值条件的主观理解,使用了特定的门限作为I M F 滤波停止条件;其二是利用三次样条计算信号的上、下包络时,根据人们对信号两端走势的主观经验,使用了特定的端点延拓方法。
当使用E MD 时,在上述两点上使用不同的规则将导致不同的分解结果。
Huang 等[1]在提出E MD 算法时给出了较好的IM F 滤波停止条件,然而该算法依然存在某些方面的不足,为了使用尽量合理的I M F 滤波停止条件,2003年R illing 等[2]对文[1]中的E MD 算法进行了改进,提出一种“新的I M F 滤波停止条件”。
实验结果表明,该改进算法可以获得更好的分解结果。
为了从理论上有效解决E MD 算法的边界效应,许多学者对端点延拓问题作了研究,这些工作包括2001年邓拥军等[3]提出的神经网络方法、2003年黄大吉等[4]提出的镜像闭合法和极值点延拓法及2004年刘慧婷等[5]提出的多项式拟合算法等。
另外,为了得到E MD 算法的解析表示,2004年Chen 等[6]提出了“直接采用基于极值点滑动平均的B 样条函数的线性组合作为均值”(滑动平均)的方法代替传统的“用极值点插值的三次样条函数分别得到信号的上下包络从而求得均值”(包络平均)的方法。
该方法获得了较好的实验结果,尤其重要的是,借助B 样条函数已有的良好性质,可以为E MD 算法中信号的低频走势(其定义参见本文§111)给出明确的解析表达式,从而为建立E MD 方法的理论基础进行了有益的探索。
但文[6]并没有讨论I M F 的滤波停止条件问题,也没有考虑E MD 算法的端点延拓问题,在那里,信号被假定是无限长的,这对实际的信号分析和处理带来不便。
本文将在这两方面对文[6]中的算法进行研究。
1 E MD 方法简介111 原始E MD 算法的基本思想E MD 算法本质上是一个有限次的滤波过程(sifting p r ocess ),使得信号具有如下两个特性:①极值点(极大值和极小值)数目与跨零点数目相等或最多相差一个(以下简称过零点条件);②由局部极大值构成的上包络和由局部极小值构成的下包络的平均值为零(以下简称均值条件)。
满足上述特征的信号就称为一个I M F 。
E MD 方法的滤波过程[1]可写成如下的算法:3收稿日期:2006202222基金项目:国家自然科学基金资助项目(60475042,10631080)作者简介:郑天翔(1979年生),男,博士生;通讯联系人:杨力华;E 2mail :mcsylh@mail 1sysu 1edu 1cn中山大学学报(自然科学版)第46卷 算法11令r (t )=x (t );2如果余量r (t )为单调函数或者其幅度小于预先给定的阈值,则算法停止,否则:3通过以下滤波过程得到各I M F,记为c (t ); 311令h (t )=r (t ); 312判断h (t )是否一个I M F?这包括过零点条件和均值条件的检验。
如果两个条件都满足,则转第4步,否则: 31211求出信号h (t )的低频走势m (t ); 31212h (t )=h (t )-m (t ),转312;4 c (t )=h (t );5 r (t )=r (t )-c (t ),转第2步;这里需要指出的是,“低频走势”并非是一个专用术语。
根据E MD 的思想,该分解是从高至低,逐层抽取出信号在各个局部的频率分量,即首先抽取信号的高频信息,剩下的就是低频信息了。
而“走势”一词反映的是信号在不同局部的差异。
112 现有几个算法变种的简单介绍在Huang 等[1]的原始算法(H98)中,m (t )是简单地取为上下包络的均值e (t )=(u (t )+v (t ))/2,其中u (t )和v (t )分别是信号的上包络和下包络。
而在判断信号h (t )是否一个I M F 时,其均值条件是加在以下这个物理量上的:连续两次滤波结果的标准偏差SD=∑Tt =0e 2(t )r 2(t )=∑Tt =0|h k (t )-h k -1(t )|2h 2k -1(t ),其中h k -1(t )、h k (t )分别表示第k 次滤波前、后的信号。
相应的滤波停止条件为SD <α,其中α∈[012,013]。
2003年R illing 等[2]对[1]中的滤波停止条件给出了改进,在R illing 的算法中,m (t )仍然是简单地取为上下包络的均值e (t ),所不同的是,其IM F 均值条件的判断是加在以下这个物理量上:σ(t )=e (t )a (t ),其中a (t )=12(u (t )-v (t ))。
相应的滤波停止条件有两个:其一是满足σ(t )<θ1的时刻个数与全部持续时间之比以至少1-α成立,即#{t ∈D |σ(t )<θ1}#{t ∈D }≥1-α,其中D 是信号的持续范围,#A 表示集合A 中元素的个数,θ1=0105,α=0105;其二是对每个时刻t ,有σ(t )<θ2,其中θ2=10θ1。
与Huang 的方法相比,σ(t )更能反映I M F 的均值特性,且两个条件相互补充,使得信号只能在某些局部出现较大的波动,从而保证了整体均值为零。
对此,在§211中将以例子详细说明。
在Chen 等[6]提出的改进算法中,对I M F 均值条件的判断仍与Huang 的方法相同,也就是使用标准偏差SD ,但m (t )取为基于极值点滑动平均的B 样条函数的线性组合,即m (t )=∑j ∈Z12k -2∑k -1l =1k -2l -1x (τj+l )B j ,k (t ),其中{τj :j ∈Z }是信号x (t )的所有极值点,B j ,k (t )是第j 个k 阶B 样条函数,其节点为{τj :j ∈Z }。
Chen 等[6]指出,除了第一个I M F 外,其余所有I M F 均是一个三次样条函数,而B 样条函数可以组成样条空间的基底,所以使用B j ,k (t )来表示m (t )是合理的。
此外,系数αj =12k -2∑k -1l =1k -2l -1x (τj+l )可近似看作是对极值点序列进行低通滤波的结果[7],代表了原信号的低频信息。
这样避免了求包络,把“极值点→包络→均值”变成更直接的“极值点→均值”。
2 对Chen 算法的探讨和改进211 R illing 改进算法中滤波停止条件的分析首先用一个例子说明R illing 改进算法的优点。
原始信号如图1(a )所示。
对原信号x (t )进行第一次滤波,得到“低频走势”m (t ),如图1(b )所示。
然后使用Huang 和R illing 方法分别产生两条均值曲线φ(t )=e (t )2r (t )2和σ(t ),并显示在图1(c )和图1(d )中。
经过多次滤波以后,最终两种方法都分别得到第一个I M F,此时的均值曲线分别如图1(e )和图1(f )所示。
从图1各图可以看出,R illing 的均值曲线(图1(d ))比Huang 的(图1(c ))更接近原信号的低频走势(图1(b )),且函数值更小,因此滤波次数大大减少,以此信号为例,使用Huang 的方法需要进行22次循环滤波才得到第一个I M F,即从图1(c )所示的结果到图1(e )所示的结果需要22次循环,而使用R illing 的方法只需要进行2次循环,也就是说,图1(f )所示的正是图1(d )的下一次滤波结果。
212 应用R illing 算法中的滤波停止准则于Chen 算法中由于R illing 的新的停止准则更优于Huang 原始的滤波停止条件,本节将它应用于Chen 提出的E MD 算法中,即把I M F 均值条件的判断加在σ(t )上,而m (t )就取为基于极值点滑动平均的三次B 样条函数的线性组合(下简写作b (t )),即2 第1期郑天翔等:经验模式分解算法的探讨和改进图1 R illing 改进算法的优点示意图Fig 11 Illustrati on of the merits for the modified algorith m devel oped by R illingm (t )=b (t )=∑j ∈Z14[x (τj+1)+2x (τj+2)+x (τj+3)]B j ,4(t )值得注意的是,由于Chen 的算法并不需要求上下包络,而σ(t )是依赖于包络的,因此,仿照R illing 的思想,我们直接使用b (t )作为σ(t ),也就是说,我们将使用如下的均值条件:#{t ∈D |b (t )<θ1}#{t ∈D }≥1-α和b (t )<θ2,其中α=0105,θ1=0105,θ2=10θ1。
