现代信号处理第八章 基于EMD的时频分析方法及其应用

2020/2/8机械工程及自动化研究所现代信号处理技术及应用第八章基于EMD的时频分析方法及其应用西安交通大学机械工程学院研究生学位课程第八章基于EMD的时频分析方法及其应用8.1EMD的基本理论和算法8.2EMD实用化技术研究8.3基于EMD的Laplace小波结构模态参数识别方法研究8.4EMD方法在机械设备故障诊断中的应用8.1EMD的基本理论和算法8.1.1EMD方法的基本概念8.1.2EMD方法的基本原理8.1.3EMD方法的完备性和正交性8.1.4基于EMD的Hilbert变换(HHT)的基本原理和算法瞬时频率的概念1()()xytdt()()()()()itztxtiytate22()()()atxtyt()()arctan()yttxt()()dttdt时间序列的Hilbert变换为:构造解析函数其中幅值函数相位函数相位函数的导数即为瞬时频率()xt(8.1.2)(8.1.1)(8.1.3)(8.1.5)(8.1.4)1()()2dtftdt(8.1.6)瞬时频率的概念(8.1.7)(8.1.8)(8.1.10)(8.1.9)(8.1.11))(2121)()()()(21tjtjtjetAeAeAtxtxtx)()()(2211AAXtAtAtAtAt22112211coscossinsinarctan)(tAAAAtA)cos(2)(122122212)()(21)(21)()(221221212'tAAAtt然而按上述定义求解的瞬时频率在某些情况下是有问题的,考虑如下信号这个信号是解析的，按式(8.1.3)和(8.1.4)可以求解其相位和幅值，得到假设信号幅值是恒定的,频率是正的，信号的频谱瞬时频率的概念当两个正弦频率取，两个频率时，幅值的取值不同，其瞬时频率亦有很大的不同。如图8.1.1(a)所示，，时，其瞬时频率是连续的。而在图8.1.1(b)中，，，虽然信号是解析的，瞬时频率却出现了负值。1012022.01A12A2.11A12A而我们已知信号的频率是离散的和正的。可见，对任一信号做简单的Hilbert变换可能会出现无法解释的频率成分。图8.1.1两个正弦波叠加的瞬时频率（a）（b）基本模式分量(IMF)的概念NordenE.Huang等人对瞬时频率进行深入研究后发现，只有满足一定条件的信号才能求得具有物理意义的瞬时频率，并将此类信号称之为基本模式分量(IMF)。基本模式分量需要满足的两个条件为：在整个数据序列中，极值点的数量(包括极大值点和极小值点)与过零点的数量必须相等，或最多相差不多于一个。在任一时间点上，信号局部极大值确定的上包络线和局部极小值确定的下包络线的均值为零。同时还提出了将任意信号分解为基本模式分量组成的经验模式分解方法(EmpiricalMODEDecomposition,EMD)基本模式分量(IMF)的概念图8.1.2一个典型的基本模式分量图8.1.2所示，是一个纯调频调幅正弦波，它满足上述两个条件，是一个典型的基本模式分量。EMD方法的基本原理和算法图中曲线：黑色—原始信号，蓝色—上包络线红色—下包络线，粉色—包络线均值()xt()xt()mt第一步确定时间序列的所有局部极值点，然后将所有极大值点和所有极小值点分别用样条曲线连接起来，得到的上、下包络线。记上、下包络线的均值为EMD方法的基本原理和算法第二步：用原始时间序列减去包络线的均值，得到,检测是否满足基本模式分量的两个条件。如果不满足，使作为待处理数据，重复第一步，直至是一个基本模式分量，记第一个基本模式分量()xt()mt1()ht1()ht1()ht1()()()htxtmt11()()ctht1()ctEMD方法的基本原理和算法第三步用原始时间序列分解出第一个基本模式分量之后，用减去，得到剩余值序列。把当作一个新的“原始序列”，重复上述步骤，依次提取出第2、第3、直至第n个基本模式分量。最后剩下原始信号的余项剩余值序列由此，时间序列可表示成n个基本模式分量和一个余项的和，即：(8.1.17)()xt1()ct1()ct()xt1()xt11()()()xtxtct()nrt11()()()xtxtct()xt()ict1()()()ninixtctrtEMD分解过程停止准则理论准则①当最后一个基本模式分量或剩余分量，变得比预期值小时便停止；②当剩余分量变成单调函数，从而从中不能再筛选出基本模式分量为止实际准则筛选过程的停止准则可以通过限制两个连续的处理结果之间的标准差的大小来实现，通常取0.2~0.32(1)20(()())()TkkdtkhthtShtEMD方法的完备性和正交性信号分解方法的完备性就是指把分解后的各个分量相加就能获得原信号的性质。通过经验模式分解方法的过程，方法的完备性已经给出，如式(8.1.17)所示。到目前为止，经验模式分解的正交性在理论上还难以严格地进行证明[17]，只能在分解后在数值上进行检验。文献[2]和[11]分别用某一齿轮箱的振动信号和某一风波信号模式分解的正交性进行了检验，结果证明EMD方法基本上是正交的，或者称是近似正交的。EMD方法的完备性表征整体正交性的指标IO(IndexofOrthogonal)定义为112011IO(()()/())Tnniktikctctxtki,220()()IO()()TikiktikctctctctOREMD方法的完备性和正交性图8.1.5小波变换与EMD方法划分信号频带（a）小波变换二进划分信号频带（b）EMD方法自适应划分信号频带常用的二进小波在对信号进行分解时，每次分解都会平分被分解信号的频带。而EMD方法则是根据信号本身具有的特性对其频带进行自适应划分，每个基本模式分量所占据的频带带宽不是人为决定的，而是取决于每个基本模式分量所固有的频率范围。xD3xD2xD1xA1xA2xA3)(tx)(3tr)(3tc)(1tc)(tx)(1tr)(2tc)(2trEMD特点EMD方法得到了一个自适应的广义基，基函数不是通用的，没有统一的表达式，而是依赖于信号本身，是自适应的，不同的信号分解后得到不同的基函数，与传统的分析工具有着本质的区别。因此可以说，经验模式分解方法是基函数理论上的一种创新。HHT方法的基本原理()()11()Re()Re()innjtdtjtiiiixtateate以上基于EMD的希尔伯特变换分析方法也称为Hilbert-Huang变换(Hilbert-HuangTransformation,HHT)。式(8.1.25)称为信号的Hilbert幅值谱，简称Hilbert谱，记做(8.1.24)(8.1.25)对式(8.1.17)中的每个IMF进行Hilbert变换可以得到其中Re表示取实部，在推导中省去了，因为它是一个单调函数或是一个常量。nr()1Re()()(,)0injtdtiiiatetHt其他8.2EMD实用化技术研究EMD分解过程的一个重要步骤就是求解信号的局部均值，这表明该方法是基于信号的局部特征的，在信号分解方法的体系中是一种概念性的创新。同时，也为我们指出了两个值得研究的方向：一是如何进一步提高局部均值的求解精度，二是如何有效地消除因边界不连续而产生的边界效应。局部均值的求解EMD方法以信号的局部极大值和局部极小值定义的包络线的均值作为信号的局部均值，只利用了信号中极值点的信息，局部均值的精度较低，且包络的求取需要两次三次样条插值，计算速度较慢。我们可以采用其它的方法来求解局部均值以提高计算的精度和速度，不同的方法对应着不同的分解过程，我们将之通称为信号模式分解技术。EMD方法中以局部极大值与极小值的包络线的均值代替信号局部均值并不是唯一的求解方法，其他求解方法有：自适应时变滤波法(ATVFD)极值域均值模式分解法(EMMD)改进的极值域均值模式分解法(IEMMD)改进的极值域均值模式分解法(IEMMD)改进的极值域均值模式分解方法(ImprovedExtremumfieldMeanModeDecomposition,IEMMD)，取消了极值域均值模式分解方法中“两极值点间的数据是均匀变化的”这一假设。首先，求出原始数据中所有局部极值点组成极值点序列)(tx)(ite再按式(8.2.1)计算出两相邻极值点间的局部均值序列im1)(1)(1iitttiiitxttttmki,,2,1其中(8.2.1)改进的极值域均值模式分解法(IEMMD)且图8.2.1信号、极值点与局部均值的关系设在原始数据中介于和之间，此处im)(jtx)(1jtx11kjjjjjjijixxttxmtt11)(],[1jjittt此时可以按式（8.2.2）求得对应的时间im（8.2.2）改进的极值域均值模式分解法(IEMMD)然后就可以用两个相邻的局部均值和加权平均求处极值点的局部均值，即)(1tmi)(21tmi1it)(1itm(8.2.3)111)()()(iiiiimtkmtktm式中和是通过相似梯形得到的加权系数，即)(itk)(1itk1212)(tttttkii12111)(tttttkii(8.2.4)求得极值点处的局部均值之后，就可以用这些点来拟合数据的局部均值曲线，进而分解出IMF。改进的极值域均值模式分解法(IEMMD)端点效应处理方法经验模式分解方法虽然能够有效的分析和处理非平稳信号，但在实际应用中存在一个比较重要的问题，就是在应用EMD方法对非平稳信号进行分解时，在数据的两端会产生发散现象，并且这种发散的结果会逐渐向内“污染”整个数据序列而使所得分解结果严重失真，这就是所谓EMD分解过程中产生的端点效应问题[2,14]。边界效应严重影响着模式分解的效果，为了解决这个问题，Huang在提出EMD方法的同时，还提出了根据特征波对原始数据进行延拓以抑制边界效应的方法，并在美国申请了专利。该特征波是由信号两端两个连续的极值点及其频率与幅值决定的，但在相关文献中并没有给出确定特征波的具体方法。端点效应处理方法目前，人们已经提出了一些抑制端点效应的方法，包括直接对原始数据进行简单延托的方法、采用神经网络对数据延托法、在端点出按照端点数据变化的“平衡位置”附加两条平行线段的方法、边界波形匹配预测法、极值点延托法、基于AR模型的时间序列线性预测方法、神经网络等，这些方法对抑制端点效应都有一定的效果。作为一种新的非线性时间序列预报方法，支持向量机(SupportVectorMachine，SVM)具有更高的预测精度[16]，可以利用该方法对时间序列进行双边延拓，在数据两端各得到若干个附加的局部极大值点和极小值点，再对模式分解后得到的各基本模式分量进行截取，从而将边界效应释放到原始数据的支撑区域外端，不影响原始数据的分析和处理。端点效应处理方法8.3基于EMD的Laplace小波结构模态参数识别方法研究8.3.1基于EMD的Laplace小波模态参数识别方法8.3.2应用实例直接采用Laplace小波相关滤波法的不足构造式(8.3.1)所示的仿真信号，来模拟单自由度结构前三阶模态的响应信号：)(txNtxtxtxtx01.0)(2.0)()(2)(321(8.3.1)其中表示第个脉冲响应信号：)(txii)12sin()(22iitfitfetxii3,2,1i(8.3.2)它们得频率分别为Hz，Hz，Hz；阻尼比分别为，，。冲击发生的时刻为0.05s，N表示幅值为1的白噪声。10001f4002f603f005.01010.02020.03直接采用Laplace小波相关滤波法的不足最