HHT变换,固有模态函数IMF,经验模式分解EMD,端点延拓,EMD结束准则,Hilbert谱

HHT变换1.1简介传统的信号处理方法，如傅立叶分析是一种纯频域的分析方法。它用频率不同的各复正弦分量的叠加来拟合原函数，也即用F来分辨f。而F在有限频域上的信息不足以确定在任意小范围内的函数f，特别是非平稳信号在时间轴上的任何突变，其频谱将散布在整个频率轴上。而且，非平稳动态信号的统计特性与时间有关，对非平稳信号的处理需要进行时频分析，希望得到时域和频域中非平稳信号的全貌和局域化结果。在傅立叶变换中，人们若想得到信号的时域信息，就得不到频域信息。反之亦然。后来出现的小波（Wavelet）变换通过一种可伸缩和平移的小波对信号变换，从而达到时频局域化分析的目的。但这种变换实际上没有完全摆脱傅立叶变换的局限，它是一种窗口可调的傅立叶变换，其窗内的信号必须是平稳的。另外，小波变换是非适应性的，小波基一旦选定，在整个信号分析过程中就只能使用这一个小波基了。HHT（Hilbert-HuangTransform）技术是(1998年由NASA的NordenEHuang等提出的新的信号处理方法。该方法适用于非线性非平稳的信号分析，被认为是近年来对以傅立叶变换为基础的线性和稳态谱分析的一个重大突破。目前HHT技术已用于地球物理学和生物医学等领域的研究，并取得了较好的结果。存在的问题尽管HHT技术在处理非线性、非稳态信号方面有很大的优势，但是这个方法本身还是有许多的问题有待进一步研究。正如Huang在文章中指出的那样，对于这种新的信号处理方法，其基的完备性还需要严密的证明。另外，在做Hilbert变换时出现的边界效应也需要更好的方法来解决。但是，HHT技术中最严重，也是现今研究的最多的是EMD分解中的包络过程。从对EMD分解方法的介绍可以看出，包络线的构造影响着整个分解的结果，也决定了后面的Hilbert变换。Huang采用的三次样条插值来拟和包络线。在实际应用中，发现这样做会产生严重的边界效应，污染了原始数据。特别是对短数据而言，这种影响使分析所得的结果失去了原有的意义。对此，Huang等采用的是根据信号端点处的振幅和频率，分别增加两组特征波的方式进行数据延拓。Huang的这种延拓方法已经向NASA申请了专利。除此之外，还有人提出了其它方法进行端点延拓。比如国家海洋局的黄大吉等提出的镜像闭合延拓法和极值延拓法。其中镜像闭合延拓法是根据信号的分布特性，把镜子放在具有对称性的极值位置，通过镜像法把镜内信号映射成一个周期性的环形信号，不存在端点，从根本上避免了端点效应。而极值延拓法具有镜像闭合延拓法相当的效果，而不增加信号序列本身的长度，计算较快，尤其在处理非对称波形信号时比较优越。还有青岛海洋大学的邓拥军等，利用神经网络分析方法来延拓数据端。这些方法对边界效应的抑制都有一定的作用，但是也都需要更深一步的研究。1.2瞬时频率频率是个极其重要的物理量，定义为信号周期倒数，其物理含义显而易见。对于正弦信号，它的频率为恒值。但是对于大部分信号，它的频率是随时间变化的函数，故提出瞬时频率概念。瞬时频率即表征信号在局部时间点上瞬态频率特性，整个持续期上的瞬时频率反映了信号频率的时变规律。对于随机时间序列X(t)，对其进行Hilbert变换[39]，可以得到Y(t)如下：1()XYtPVdt(1-1)其中，PV为柯西主值（Cauchyprincipalvalue)。该式表示Y(t)是X(t)与1的卷积。通过这个定义，X(t)和Y(t)组成了一个共轭复数对，于是可以得到一个解析信号Z(t)如下itZtXtiYtate(1-2)其中1222arctanatXtYtYttXt(1-3)从理论上讲，虚部的定义方法有很多种。但是Hilbert变换为其提供了一个唯一的虚部值，这就使得其结果成为一个解析函数。得到了相位，就可以得到瞬时频率，因为瞬时频率就是相位导数。dtdt(1-4)从本质上说，公式(1-1)将Hilbert变换定义为X(t)与1的卷积，因此它强调了X(t)的局部特性。公式(1-2)的极坐标表达式更进一步澄清了这个表达式的局部特性：它的幅值和相位的最佳局部匹配随三角函数变化。1.3固有模态函数IMF（IntrinsicModeFunction）由瞬时频率的物理意义可知，并不是任意的信号都能用瞬时频率来讨论。只有当信号满足只包括一种振动模式，而没有复杂叠加波的情况时才行。实际上，定义一个有意义的瞬时频率的必要条件就是要求函数关于局部零平均值对称，并且零交叉点和极值点数量相同。基于此种原因，提出了固有模态函数的概念。固有模态函数满足以下两个条件：(1)整个数据范围内，极值点和过零点的数量相等或者相差一个；(2)在任意点处，所有极大值点形成的包络线和所有极小值点形成的包络线的平均值为零。第一个条件是显而易见的，它类似于平稳过程中传统的稳定且满足高斯分布的窄带信号条件。第二个条件把传统的全局条件调整到局部情况。只有满足了这个条件，得到的瞬时频率才不会因为不对称波形的存在而引起不规则波动。所以这一点是得到正确瞬时频率的必要条件。而这一点是必须的，因为这样瞬时频率就可以不包含由于不对称波形造成的波动。为了使用瞬时频率定义，必须要把随机数据归结为IMF组件，这样才可以为每个IMF组件定义瞬时频率。为了将数据归结为所需的IMF组件，接下来引入经验模式分解方法。1.4经验模式分解EMD（EmpiricalModeDecomposition）经验模态分解方法的大体思路是利用时间序列上下包络的平均值确定“瞬时平衡位置”，进而提取固有模态函数[39]。这种方法基于如下假设：（1）信号至少有两个极点—－个极大值和一个极小值；（2）信号特征时间尺度是由极值间的时间间隔来确定的；（3）如果数据没有极值而仅有拐点，可以通过微分、分解、再积分的方法获得IMF。在此假设基础上，Huang等人进一步指出：可以用经验模态分解方法将信号的固有模态筛选出来。经验模式分解过程就是个筛选过程，实现振动模式的提取。该方法的基本思路是用波动上、下包络的平均值去确定“瞬时平衡位置”，进而提取出固有模态函数。上、下包络线是由三次样条函数对极大值点和极小值点进行拟合得到的。经验模式分解过程的基本过程可概括如下：(1)寻找信号x(t)所有局部极大值和局部极小值，为更好保留原序列的特性，局部极大值定义为时间序列中的某个时刻的值，它只要满足既大于前一时刻的值也大于后一时刻的值即可。局部极小值的提取同理，即该时刻的值满足既小于前一时刻的值也小于后一时刻的值。使用三次样条函数进行拟合，获得上包络线xmax(t)和下包络线xmin(t)；(2)计算上、下包络线的均值m(t)=[xmax(t)+xmin(t)]/2；(3)用原信号x(t)减去均值m(t)，得到第一个组件h(t)=x(t)-m(t)；由于原始序列的差异，组件h(t)不一定就是一个IMF，如果h(t)不满足固有模态函数两个条件，就把h(t)当成原始信号，重复(1)-(3)，直到满足条件为止，这时满足固有模态函数条件的h(t)作为一个IMF，令I1(t)=h(t)，至此第一个IMF已经成功的提取了。由于剩余的r(t)=x(t)-I1(t)仍然包含具有更长周期组件的信息，因此可以把它看成新的信号，重复上述过程，依次得到第二个I2(t)，第三个I3(t)，…，当r(t)满足单调序列或常值序列条件时，终止筛选过程，可以认为完成了提取固有模态函数的任务，最后的r(t)称为余项，它是原始信号的趋势项。由此可得x(t)的表达式1niixtItrt，即原始序列是由n个IMF与一个趋势项组成。如上所述，整个过程就像筛选过程，根据时间特性把固有模态函数从信号中提取出来。1.5端点延拓HHT方法的分析质量很大程度上取决于EMD分解的质量，而在应用EMD方法时的一个非常棘手的问题是，由于信号两端不可能同时处于极大值和极小值，这样，“筛”过程中构成上下包络的三次样条函数在数据序列的两端就会出现发散现象（如图1.1(a)所示，原始信号是中心频率为1KHz、5个波峰、以100个数据描绘一周期的窄带信号），并且这种发散的结果会随着“筛”过程的不断进行逐渐向内“污染”整个数据序列而使所得结果严重失真（如图1.1(b)所示）。另外，在进行Hilbert变换时，信号的两端也会出现严重的端点效应。对于一个较长的数据序列来讲，可以根据极值点的情况不断抛弃两端的数据来保证所得到的包络的失真度达到最小。但对于一个数据点数少的序列来讲，这样的操作就变得完全不可行。因此，必须对信号或其极值向外进行延拓，以确保包络线抵达端点。端点延拓的目的是确保上、下包络都与端点相交，以便有与每一个信号点相对应的局部平均值。而上、下包络是由极大值和极小值连结而成的，因此只要对极大值和极小值进行延拓，而不必对信号本身进行延拓。极大值和极小值是相间分布的，同时考虑到样条插值的要求，所以只要在信号左、右两端分别延拓两个极大值和两个极小值即可。由于端点以外没有信号，任何延拓都是人为的，通过对不同形式的信号的多次试验，发现一种端点延拓的方法以端点的一个特征波为依据进行延拓具有较好的效果，具体做法如下[40]：设离散信号：12[(1),(2),,()][,,,]nttttnttt(1-5)1212()[(),(),,()][,,,]nnXtxtxtxtxxx(1-6)其采样步长为t,()Xt有M个极大值和N个极小值，对应的序列下标（mI,nI）、时间（mT,nT）和函数值（U，V）记为：[(1),(2),,()]mmmmIIIIM(1-7)[(1),(2),,()]nnnnIIIIN(1-8)II(),(),1,,mmmTitUixiM(1-9)II(),(),1,,nnnTitVixiN(1-10)具体作法如下：1.左端信号左端第一个特征波包含的信号点数为1k（b）EMD分解(a)拟合结果原始信号上下包络线图3.1未作端点延拓的分析效果(a)拟合结果（b）EMD分解原始信号上下包络线图3.2以端点的一个特征波为依据进行延拓的分析效果1(2)(1)(1)(1)(2)(1)(1)(1)2|(1)(1)|mmmnnnmnmnIIIIkIIIIIIMN，当，当，当==1(1-11)向外延拓的两个极值的位置（mT,nT）和数值（U，V）为：1(0)(1),(0)(1)mmTTktUU(1-12)1(1)(1)2,(1)(1)mmTTktUU(1-13)1(0)(1)(0)(1)nnTTktVV(1-14)1(1)(1)2(1)(1)nnTTktVV(1-15)2.右端信号右端第一个特征波包含的信号点数为k22()(1)()()()(1)()()2|()()|mmmnnnmnmnIMIMIMINkININIMINIMINMN，当，当，当==1(1-16)向外延拓的两个极值的位置（mT,nT）和数值（U，V）为：2(1)(),(1)()mmTMTMktUMUM(1-17)2(2)()2,(2)()mmTMTMktUMUM(1-18)2(1)(),(1)()nnTNTNktVNVN(1-19)2(2)()2,(2)()nnTNTNktVNVN(1-20)3.当端点的数值比近端点的第一个极大值大或极小值小时，要进行特殊的处理，以避免信号落到包络线之外，111(0),(0),(1)mTtUxxU当(1-21)111(0),(0),(1)nTtVxxV当(1-22)(1),(1),()mnnnTMtUMxxUM当(1-23)(1),(1),()nnnnTNtVNxxVN当(1-24)由图1.2可见，以端点的一个特征波为依据进行延拓，分别在信号两端增加