HHT变换讲义EMD剖析.docx
《HHT变换讲义EMD剖析.docx》由会员分享，可在线阅读，更多相关《HHT变换讲义EMD剖析.docx（14页珍藏版）》请在冰点文库上搜索。
HHT变换讲义EMD剖析
HHT变换讲义
1.1简介
传统的信号处理方法，如傅立叶分析是一种纯频域的分析方法。
它用频率不同的各复正弦分量的叠加来拟合原函数，也即用
来分辨
。
而
在有限频域上的信息不足以确定在任意小范围内的函数
，特别是非平稳信号在时间轴上的任何突变，其频谱将散布在整个频率轴上。
而且，非平稳动态信号的统计特性与时间有关，对非平稳信号的处理需要进行时频分析，希望得到时域和频域中非平稳信号的全貌和局域化结果。
在傅立叶变换中，人们若想得到信号的时域信息，就得不到频域信息。
反之亦然。
后来出现的小波（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-1）
其中，PV为柯西主值（Cauchyprincipalvalue）。
该式表示Y（t）是X（t）与
的卷积。
通过这个定义，X（t）和Y（t）组成了一个共轭复数对，于是可以得到一个解析信号Z（t）如下
（1-2）
其中
（1-3）
从理论上讲，虚部的定义方法有很多种。
但是Hilbert变换为其提供了一个唯一的虚部值，这就使得其结果成为一个解析函数。
得到了相位，就可以得到瞬时频率，因为瞬时频率就是相位导数。
（1-4）
从本质上说，公式（1-1）将Hilbert变换定义为X（t）与
的卷积，因此它强调了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）的表达式
，即原始序列是由n个IMF与一个趋势项组成。
如上所述，整个过程就像筛选过程，根据时间特性把固有模态函数从信号中提取出来。
1.5端点延拓
HHT方法的分析质量很大程度上取决于EMD分解的质量，而在应用EMD方法时的一个非常棘手的问题是，由于信号两端不可能同时处于极大值和极小值，这样，“筛”过程中构成上下包络的三次样条函数在数据序列的两端就会出现发散现象（如图1.1（a）所示，原始信号是中心频率为1KHz、5个波峰、以100个数据描绘一周期的窄带信号），并且这种发散的结果会随着“筛”过程的不断进行逐渐向内“污染”整个数据序列而使所得结果严重失真（如图1.1（b）所示）。
另外，在进行Hilbert变换时，信号的两端也会出现严重的端点效应。
对于一个较长的数据序列来讲，可以根据极值点的情况不断抛弃两端的数据来保证所得到的包络的失真度达到最小。
但对于一个数据点数少的序列来讲，这样的操作就变得完全不可行。
因此，必须对信号或其极值向外进行延拓，以确保包络线抵达端点。
端点延拓的目的是确保上、下包络都与端点相交，以便有与每一个信号点相对应的局部平均值。
而上、下包络是由极大值和极小值连结而成的，因此只要对极大值和极小值进行延拓，而不必对信号本身进行延拓。
极大值和极小值是相间分布的，同时考虑到样条插值的要求，所以只要在信号左、右两端分别延拓两个极大值和两个极小值即可。
由于端点以外没有信号，任何延拓都是人为的，通过对不同形式的信号的多次试验，发现一种端点延拓的方法以端点的一个特征波为依据进行延拓具有较好的效果，具体做法如下[40]：
设离散信号：
（1-5）
（1-6）
其采样步长为
有
个极大值和
个极小值，对应的序列下标（
）、时间（
）和函数值（
，
）记为：
（1-7）
（1-8）
（1-9）
（1-10）
具体作法如下：
1.左端
信号左端第一个特征波包含的信号点数为
（1-11）
向外延拓的两个极值的位置（
）和数值（
，
）为：
（1-12）
（1-13）
（1-14）
（1-15）
2.右端
信号右端第一个特征波包含的信号点数为k2
（1-16）
向外延拓的两个极值的位置（
）和数值（
，
）为：
（1-17）
（1-18）
（1-19）
（1-20）
3.当端点的数值比近端点的第一个极大值大或极小值小时，要进行特殊的处理，以避免信号落到包络线之外，
（1-21）
（1-22）
（1-23）
（1-24）
由图1.2可见，以端点的一个特征波为依据进行延拓，分别在信号两端增加两个极大值和两个极小值，从而使原始信号被延长的包络线所限制，有效地抑制了端点效应，使得EMD分解得到合理的各个IMF模态。
另外，国家海洋局的黄大吉等还提出了镜像闭合延拓法，即根据信号的分布特性，把镜子放在具有对称性的极值位置，通过镜像法把境内信号映射成一个周期性的环形信号。
这样做就不存在端点，从根本上避免了端点效应[40]。
青岛海洋大学的邓拥军等人利用神经网络分析的方法来进行端点延拓，在数据的两端各得到两个附加的极大值点和两个附加的极小值点，由此有效的抑制了端点效应[41]。
1.6EMD结束准则
根据理论分析IMF必须满足两个条件，但是在实际筛选过程中，严格满足这两个条件的信号几乎不存在，因此如果以这两个条件为依据判定IMF，可能根本就得不到结果或者要以冗长的程序执行时间为代价。
因此筛选过程要小心进行，为了保证IMF组件的调幅和调频都具有物理意义，同时考虑到程序可行性，必须制定一个结束筛选的准则。
习惯上，可以通过标准偏差
来完成，它可以由连续两个筛选结果得到[18]：
（1-25）
h1（k-1）（t）和h1k（t）分别表示两个连续的筛选结果，一般来说，SD值越小，所得的IMF分量的线性和稳定性就越好，典型的SD值为0.2-0.3，这样对于不同筛选过程就有严格的SD限制了。
但是在筛选过程中发现，SD的这个范围过于死板，对很多信号而言，以0.2-0.3为限制的
值并不科学，因此在实际筛选过程中，可以以0.2-0.3为参考值，根据实际情况进行适当调整。
另外，如果不是进行信号精确处理，可以考虑控制筛选次数来决定筛选结果，这样不仅保障了程序的有效性，也控制了程序的执行时间，在具体问题中取得了很好的效果。
1.7Hilbert谱
完成了经验模式分解过程就得到了所有可提取的固有模态函数，在Hilbert变换的基础上，只要根据公式（1-1）-（1-4）计算瞬时频率就可以了。
对固有模态函数进行Hilbert变换后，可以用下面方式表示信号
（1-26）
在这儿没有考虑余项rn，因为它只是单调函数或常量。
尽管Hilbert变换可以把单调函数看成是振动的一部分，但是其余项中的能量很小，一般情况可以不用考虑。
由希尔伯特变换得出的振幅和频率都是时间的函数，如果用三维图形表达幅值、频率和时间之间的关系，或者把振幅用灰度的形式显示在频率－时间平面上，就可以得到Hilbert谱
。
如果把
对时间积分，就得到希尔伯特边际谱
：
（1-27）
边际谱提供了对每个频率的总振幅的量测，表达了整个时间长度内累积的振幅。
另外，作为希尔伯特边际谱的附加结果，可以得到如式（1-28）定义的希尔伯特瞬时能量：
（1-28）
瞬时能量提供了信号能量随时间的变换情况。
事实上，如果振幅的平方对时间积分，可以得到希尔伯特能量谱：
（1-29）
希尔伯特能量谱提供了对于每个频率的能量的量测，表达了每个频率在整个时间长度内所累积的能量。
1.8HHT方法的软件实现
HHT变换的整个过程包括EMD与Hilbert变换，整个过程是通过matlab软件编写的M文件实现的，其程序流程图如图1.3所示。
1.9HHT方法的优越性
1、HHT方法是一种全新的信号分析方法，能够描绘出信号的时频谱图、边际谱，能量谱等，是一种更具有适应性的时频局域化分析方法[28]。
而传统的信号处理方法，如傅立叶变换是一种纯时域的分析方法，它用频率从零到无穷大的各复正弦分量的叠加来拟合原函数f（t）在每个时刻的值，也即用F（ω）在有限频域上的信息不足以确定在很小范围内的函数f（t），特别是非平稳信号在时间轴上的任何突变，其频谱将散布在整个频率轴上[39]。
所以，这种分析方法适用于确定性的平稳信号，而在非线性、非平稳过程的处理上，傅立叶变换将只作为一种数学变化手段。
而且，非平稳信号的统计特性与时间有关，对非平稳信号的处理需要进行时频分析，希望得到时域和频域中非平稳信号的全貌和局域化结果。
但在傅立叶变换中，若想得到信号的时域信息，就得不到频域信息，反之亦然。
后来出现的小波变换通过一种可伸缩和平移的基小波对信号变换，从而达到时域局域化分析的目的，但这种变化实际上并没有完全摆脱傅立叶变换的局限，它是一种窗口可调的傅立叶变换，其窗内的信号必须是平稳的[12]。
2、HHT局部性能良好而且是自适应的，对平稳信号和非平稳信号都能进行分析。
它没有固定的先验基底，分解完全基于数据本身进行。
固有模态函数是基于序列数据的时间特征尺度得出的不同的数据有不同的固有模态函数，每个固有模态函数可以认为是信号中固有的一个模态，所以通过Hilbert变换得到的瞬时频率具有清晰的物理意义，能够表达信号的局部特征。
而傅立叶变换是以余弦函数为基底进行信号分解的，这样难免产生许多实际上并不存在的虚假分量。
小波变换是选定适当的小波基，再对信号进行分解的，且小波基一但选定，在整个信号分析过程中就只能用这个基。
另外，小波变化的解释也不直观，用小波变化时，为了确定某一信号的局部变化，即使这一局部变化只是发生在低频范围内，也必须从高频范围内开始去寻找这一结果，频率越高，小波变换得局域化性能就越好[28]。
3、式（1-26）中信号幅值和频率都是时间函数，但同样的信号如果用傅立叶变换表示的话，可得
，（1-30）
其中，
和
是常量。
公式（1-26）和（1-30）的差别在于，IMF表示了广义的傅立叶扩展，可变的幅值和瞬时频率不仅强化了信号信息还使之适用于非平稳信号。
通过IMF可以清楚的区分调幅和调频，这样就打破了傅立叶变换中固有幅值和固有频率的限制，允许分解出的IMF的幅值随时间变化，使信号分析更加灵活方便。
同时，由于Hilbert变化通过微分法来定义瞬时频率，因此不需要大量的信号点来定义振动。
4、Hilbert谱把各IMF分量的幅值以灰度的形式表示在频率—时间图上，其中幅值以点的灰度表示：
点越亮，幅值越大；点越暗，幅值越小。
Hilbert谱三维灰度图的形式对各分量的瞬时频率和瞬时振幅都进行了比较确切的刻画[42]。
比传统的频谱图更直观清晰。
5、很多信号不是围绕水平位置振动的，这就需要从测量数据中去掉这个背景场，使之成为水平直线附近的振动[43]。
EMD方法就可以有效地提取或取出数据序列的均值，消除序列的趋势项，把复杂的数据分解成若干线性、平稳的模态，且不改变原数据的物理特性。
即式（1-31）。
而小波变换和傅立叶变换，它们的基底都是以水平轴为平衡位置的波动现象。
（1-31）
升级会员 