hh中的包络线和分段幂函数法

hh中的包络线和分段幂函数法

本征模函数方法信号分析方法广泛应用于振动工程中。1998年美国科学家Norden.E.Huang等人提出了一种主要用于分析非平稳信号的新信号分析方法——希尔伯特-黄变换(Hilbert-Huangtransform,HHT),其主要创新是本征模函数(Intrinsicmodefunction,IMF)概念的提出和经验模态分解(Empiricalmodefunction,EMD)方法的引入:通过EMD,将信号分解成IMF(一般为有限数目)的和,对每个IMF分别进行Hilbert变换谱分析就可以获得有意义的瞬时频率,从而给出非平稳信号中频率随时间变化的精确表达。目前这一方法已得到许多应用,正吸引越来越多的人对其进行研究。HHT的关键问题是通过信号极值点拟合信号的包络线,但目前采用的包络线(拟合)算法还存在一些不足,需要进一步研究与完善。1信号极值的估计HHT中的IMF是满足以下两个条件的信号:(1)整个信号中,零点数与极点数相等或至多相差1;(2)信号上任意一点,由局部极大值点确定的包络线和由局部极小值点确定的包络线的均值为零。用EMD将信号分解成IMF的和的步骤是:先识别出信号(设为s(t))的所有极值点,然后用所有极大值点和所有极小值点分别拟合出s(t)的上包络线esup(t)和下包络线elow(t),满足elow(t)≤s(t)≤esup(t)(1)elow(t)≤s(t)≤esup(t)(1)记上、下包络线的平均值为m(t)m(t)=(esup(t)+elow(t))/2(2)m(t)=(esup(t)+elow(t))/2(2)将s(t)减去m(t)得c(t)=s(t)-m(t)(3)c(t)=s(t)−m(t)(3)将c(t)视为新的s(t),反复重复以上过程就可以筛选出原信号中的第一阶IMFc1(t)(由判止准则判断)。m1(t)=s(t)-c1(t)(4)m1(t)=s(t)−c1(t)(4)将m1(t)视为新的s(t),用同样的方法可以依次筛选出原信号中的其他IMF。s(t)最终可表示为s(t)=n∑k=1ck(t)+mn(t)(5)s(t)=∑k=1nck(t)+mn(t)(5)其中mn(t)称为s(t)的余项,代表信号的平均趋势。对式(5)中的每个IMF应用式(6)fik(t)=12πddt[argzk(t)](6)fik(t)=12πddt[argzk(t)](6)就可以分别获得它们的有物理意义的瞬时频率,其中上标“i”表示瞬时,zk(t)为分量ck(t)的解析信号,它通过Hilbert变换计算,即zk(t)=ck(t)+jΗ[ck(t)](7)zk(t)=ck(t)+jH[ck(t)](7)其中H[ck(t)]为ck(t)的Hilbert变换,即Η[ck(t)]=ck(t)*1πt=Ρ.V.∫∞-∞ck(t-τ)πτdτ(8)这里P.V.表示柯西主值积分。可以看出,通过信号极值点拟合信号包络线的算法影响着EMD的全过程,最终也影响着瞬时频率的计算,因而是HHT中的关键问题。但如何从理论上严格确定HHT中的包络线算法是迄今还未解决的问题,目前采用的是经验算法,其中最常用的经验算法是三次样条插值法。实践证明,三次样条插值法的确是一种相对较理想的包络线算法,在许多情况下能获得比较满意的分析效果。但三次样条插值法容易出现过冲问题,如图1(a)中A,B和C处所示。“过冲”不仅容易使上、下包络线的平均值曲线发生偏移,甚至使EMD过程所要求的条件式(1)也得不到满足,如图1(a)中的B处所示,因此需对HHT中的包络线算法进行深入探索和改进。作者在研究中发现,三次样条插值法容易造成过冲的根本原因是其“柔性”不够。事实上,对于一个有经验的人(如木工师傅)来说,对一条给定的曲线,可以凭经验画出其上、下两条“很好”的包络线。例如,对于图1(a)中A,B和C处所出现的过冲,木工师傅们在该处会设法使包络线更具“柔性”一些。但在算法中,“柔性”和“光滑性”是一对矛盾,因此算法中使包络线“更具柔性”意味着要牺牲“光滑性”,而三次样条插值曲线是二阶光滑的(即二阶连续可导),因此作者想到用只有一阶光滑性(即插值曲线一阶连续可导)的典型插值法——Akima插值法进行实验。Akima插值法依据如下原理:对于一个有经验的人来说,当它在考虑由一系列数据点所确定的曲线时,往往只要注意所需的一段曲线附近的几个数据点就行,因为这几个数据点的位置完全可以描绘出这段曲线的形状。所以用Akima插值法来作为HHT中的包络线算法具有一定的可行性。但作者经过实验发现,虽然用Akima插值法是可行的,但有时所拟合的曲线又显得“过于柔软”,甚至出现明显的“折点”,如图1(b)中的A,B和C处所示。2分段物理法的原理在HHT中的包络线算法问题上,三次样条插值法柔性不够,而Akima插值法“过于柔软”(即光滑性不够)。下面先介绍曲线拟合中的抛物线参数样条插值法。对插值点P1,P2,…,Pn进行抛物线参数样条插值的步骤是:将所有插值点中任意3个相邻点Pi-1,Pi,Pi+1用抛物线参数样条插值,将Pi,Pi+1,Pi+2也用抛物线参数样条插值。设所得抛物线分别为Si(t),Si+1(t),则Pi与Pi+1之间的曲线用Si(t),Si+1(t)重新插值即得抛物线参数样条曲线。Pi与Pi+1之间的最终插值曲线的参数形式为S(t)=(1-2t)Si(t+0.5)+2tSi+1(t)(9)或S(t)=(-t+4t2-4t3)Ρi-1+(t-10t2+12t3)Ρi+(t+8t2-12t3)Ρi+1+(-2t2+4t3)Ρi+2(10)其中0≤t≤0.5。已经证明这样的曲线是一阶光滑曲线。抛物线参数样条插值法是一种常用的曲线拟合算法,如MicrosoftWord中的曲线绘图就是使用的该算法。用抛物线参数样条法所拟合出的曲线既具有很好的柔性,又具有很好的光滑性,因此具有HHT中包络线算法所希望的性质,但抛物线参数样条插值法不能保证所插值曲线是单值的,如图2中A,B,C,D为插值点,插值曲线在x0处显然是多值的。作者利用抛物线参数样条插值法的原理提出一种新的单值插值法——分段幂函数法。分段光滑定理:设P1(x1,y1),P2(x2,y2),P3(x3,y3),P4(x4,y4)为平面内4点,且x1≤x2≤x3≤x4,S2(x),S3(x)分别为经P1,P2,P3和P2,P3,P4所插值出的单值光滑曲线(即其一阶连续可导),则曲线S(x)=x3-xx3-x2S2(x)+x-x2x3-x2S3(x)(11)在区间(x1,x4)内一阶光滑。证明:显然只需证明S(x)在点P2,P3一阶连续可导即可。由式(11)得S′(x)=1x3-x2(-S2(x)+(x3-x)S′2(x)+S3(x)+(x-x2)S′3(x))所以S′(x2)=1x3-x2(-S2(x2)+(x3-x2)S′2(x2)+S3(x2))=1x3-x2(-y2+(x3-x2)S′2(x2)+y2)=S′2(x2)故S(x)在点P2一阶连续可导,同理可证S(x)在点P3一阶连续可导。分段幂函数法的步骤与抛物线参数样条插值法类似,只是分段幂函数法采用一种幂函数的方法对插值点P1,P2,…,Pn中的相邻点Pi-1,Pi,Pi+1和Pi,Pi+1,Pi+2分别插值,再依据式(11)对两段曲线进行“拼接”。在幂函数法中对任意3个插值点P1,P2,P3插值的方法如下:如图3(a)所示,设过点P2且与P1P3平行的直线分别交直线x=x1和x=x3于点A和B,先构造幂函数y=xβ(β>1且β的取值使y=xβ的定义域为(-∞,+∞),值域为[0,+∞)),再构造函数y={→A′Ρ′1→A′Cxβx≤0→B′Ρ′3→B′Dxβx≥0={→AΡ1(x1-x2)βxβx≤0→BΡ3→(x3-x2)βxβx≥0={(xx1-x2)β→AΡ1x≤0(xx3-x2)β→AΡ1x≥0(12)其中A′P′1=AP1,B′P′3=BP3,C,D分别为曲线y=xβ与直线x=x1-x2和x=x3-x2的交点,如图3(b)所示。设过点P2且与P1P3平行的直线的方程为y=g(x),式(12)最右边的函数为f(x),则对插值点P1,P2,P3采用式(13)进行插值y=f(x-x2)+g(x)(13)即S2(x)={(x-x2x1-x2)β→AΡ1+y3-y1x3-x1(x-x2)+y2x≤x2(x-x2x3-x2)β→AΡ1+y3-y1x3-x1(x-x2)+y2x≥x2(14)其中→AΡ1=(x3-x2)y1-(x2-x1)y3x3-x1。容易证明,由式(14)确定的S2(x)为单值光滑曲线。分段幂函数法插值中参数β的取值很关键,它可以是整数和小数,可以有无穷多种取值,但作者经过大量实验发现当β取2.5时效果较为理想,这大概是因为β略低于三次样条插值中各分段函数的次数(3次),因而更具柔性。3包络线算法hhd分析图4给出了用分段幂函数法对图1(a,b)中同一信号s所拟合出的包络线,从图中A,B,C,D处可以看出,新包络线合算法使三次样条插值法的过冲问题和Akima插值法的“折点”问题均得到了较明显的改善。图5(a～c)分别给出了用三次样条插值法、分段幂函数法和Akima插值法作为包络线算法时对图1中信号s进行EMD分解的结果,图6则给出了它们各自前两阶IMF分量的瞬时频率(由于第三阶IMF分量不是主要分量,故这里不对其瞬时频率进行分析)。从图5,6可以看出,用3种包络线算法所分析出的结果在总体上很一致,表明3种算法总体上都比较合理;但在细节部分则存在明显差异。从图5中还难以说明这种差异代表了孰优孰坏,但从图6(a,c)可以看出,采用三次样条插值法时所得到的第一阶IMF分量的瞬时频率在图中A和B处出现了明显跳跃,采用Akima插值法时所得到的第二阶IMF分量的瞬时频率在图中C和D处出现了较大波动,而从图1中却看不出原信号具有这种急剧变化的频率(实际上,图1中的信号是用三次样条曲线所拟合出的曲线,理论上也不存在急剧变化的频率),这说明,用三次样条插值法或Akima插值法作为包络线算法进行HHT分析时会出现明显的虚假频率,其原因正是由于三次样条插值法存在过冲,而Akima插值法存在“折点