EMD局部积分均值增密插值改进算法及其在转子故障诊断中的应用

蔡艳平, 徐光华, 李艾华, 王 涛

(1.西安交通大学 机械工程学院，西安 710049；2.西安交通大学 机械制造系统工程国家重点实验室，西安 710049；3.火箭军工程大学，西安 710025)



EMD局部积分均值增密插值改进算法及其在转子故障诊断中的应用

蔡艳平1，3, 徐光华1，2, 李艾华3, 王 涛3

(1.西安交通大学 机械工程学院，西安 710049；2.西安交通大学 机械制造系统工程国家重点实验室，西安 710049；3.火箭军工程大学，西安 710025)

包络拟合是EMD(Empirical Mode Decomposition)算法的一个关键环节，针对EMD包络拟合问题，提出一种基于局部积分均值增密约束三次样条插值的EMD改进算法。该方法利用定比分点法来增密EMD样条插值型值点，利用分段约束三次样条插值来抑制传统EMD包络过冲与欠冲问题，利用内禀模态函数直接筛选方法来减少EMD计算局部均值时的样条插值次数。仿真分析验证了改进算法的有效性，将改进算法应用在转子故障振动信号的实例分析中，结果表明，改进算法提高了EMD的分解精度，更加准确地提取了油膜涡动信号的故障特征。

旋转机械; 故障诊断; 经验模态分解; 样条插值; 包络拟合

EMD[1]是由美国航空航天局 HUANG于1998年提出，被认为是近年来应用最广泛的时间序列分析方法。由于EMD是完全自适应的，尤其适用于非线性和非平稳过程分析，因此被广泛应用在机械的故障诊断中[2-5]。EMD包络拟合是EMD算法中的一个关键问题，EMD筛分出的IMF几乎都与包络拟合方法有关。因此，计算各层信号极值点包络均值曲线是EMD的核心所在，它的好坏直接影响到EMD的分解结果。

在标准EMD中，EMD包络曲线拟合是采用三次样条插值求取的，存在过冲、欠冲问题和筛分速度慢等问题，为解决EMD包络拟合问题，许多学者提出了改进方法，大致分为以下二类：

(1) 包络平均法。用极大、极小极值点分别插值求包络，再采用包络平均作为极值点局部均值的方法，如三次样条插值法[1]、阿克玛(Akima)插值法[6]、分段幂函数插值法[6-7]、高次样条插值法[8]等。这些方法均各有所长，也各有所短，如阿克玛插值法是一种自适应性比三次样条插值好的非均匀点曲线插值算法，计算量小，应用灵活，但是其光滑性较差，拟合的曲线有时显得“过于柔软”，甚至出现明显的“折点”。分段幂函数法与三次样条插值和阿克玛插值法相比，分段幂函数法光滑性较好，同时具有较好的抗过冲和欠冲性能，但其光滑性较三次样条插值差。高次样条插值法虽然可以提高曲线拟合和EMD分解精度，但计算复杂性大，耗时多，而且算法精度的提高并不是与插值次数成正比，而是到一定的阶数后，精度不再提高。另外，该方法并未解决三次样条插值带来的过冲和欠冲问题，只是提高了算法的精度而已。

(2) 局部均值直接求取法。通过对极值点直接求取局部均值时间序列，再通过对局部均值时间序列进行一次插值来求取局部均值的方法，如B样条插值法[9]、自适应时变滤波法[10]、局域中值定理分解方法[11]、稳定点均值筛法[12]。上述这些方法虽减少了运算量和1次三次样条的过冲及欠冲问题，但是通过极值点构建的局部均值时间序列中，局部均值位置tξ的粗略估计虽然简化了算法，却影响了加权系数的准确性，导致局部均值的估计仍存在误差。另外，对于局部均值的包络曲线拟合，文献[10-11]仍采用了传统三次样条，而三次样条对非均匀插值点，易产生过冲和欠冲现象，进而破坏信号的模态结构，造成分解误差，而且还会造成EMD的过分解，从而产生虚假的频率成分，影响瞬时频率的计算，因此EMD中的包络拟合仍需要进一步完善。

针对EMD中三次样条插值的过冲、欠冲问题和筛分速度等问题，本文拟从三个不同角度出发，提出一种基于局部积分均值增密约束三次样条插值的EMD改进算法，并成功应用于转子油膜涡动的故障监测诊断中。

1 EMD包络拟合问题分析

在标准EMD中，EMD包络曲线拟合是采用三次样条插值求取的，即先找出信号的所有极大值和极小值点，然后分别用三次样条插值法来对极大值点、极小值点序列插值，从而获得上、下包络，再由上、下包络线求信号的局部均值，以筛出信号的IMF。这里可通过一个仿真信号的EMD分解来分析其“包络拟合“的问题，仿真信号的表达式为：

x(t)=cos(2π×3t)+0.4cos(2π×10t)+

0.9sin(2π×25t)

(1)

仿真信号采样频率为1 024，数据长度为100。仿真信号的三次样条包络曲线插值见图1，其中实线是原始信号，虚线是基于三次样条插值求取的包络均值。

图1 EMD中三次样条插值的过冲与欠冲问题Fig.1 The problem of the cubic spline interpolation overshoot and undershoot in EMD

从图1可以看出，三次样条插值产生了过冲和欠冲现象，破坏了EMD的模态结构，造成了较大误差，且其上、下包络曲线还产生了端点效应问题。传统EMD上下包络均值求解存在的两大不足：①三次样条插值的过冲欠冲问题；②一次包络均值的计算需2次三次样条拟合，计算量大。

2 EMD包络拟合的局部积分均值增密插值改进算法

在EMD包络拟合问题分析基础上，基于插值拟合算法改进的思想，从消除三次样条插值的过冲、欠冲影响出发，本文在EMD包络曲线拟合中首先引入一种基于约束条件的分段约束三次样条插值，约束三次样条插值算法详见文献[13]。为了介绍分段插值算法，这里给出一个定理。

定理1[14]已知P0(x0,y0),P1(x1,y1),… ,Pn(xn,yn),n≥3，且x0≤x1≤…≤xn，S(x;xi-1,xi,xi+1)为经过插值点Pi-1,Pi,Pi+1，1≤i≤n-1的插值函数，且S(x;xi-1,xi,xi+1)一阶连续可导，令：

x∈[xk,xk+1],1≤k≤n-2

(2)

则S(x)是关于点P0(x0,y0),P1(x1,y1),… ,Pn(xn,yn),n≥3的插值函数，且S(x)∈C1(x0,xn)。

定理1的证明见文献[14]，定理1所描述的实际上是一种分段插值法，每个分段为一个插值函数的加权平均，本文将约束三次样条插值应用在分段插值法中，提出一种分段约束三次样条插值算法，其算法如下：

步骤1 设插值点P0(x0,y0),P1(x1,y1),…,Pn(xn,yn),n≥3，且x0≤x1≤…≤xn，采用约束三次样条插值分别求出关于插值点Pi-1,Pi,Pi+1，1≤i≤n-1的约束三次样条插值函数；

步骤2 用式(2)求取P0(x0,y0),P1(x1,y1),…,Pn(xn,yn),n≥3的插值函数S(x)；

步骤3 输出分段约束三次样条插值算法。

图2给出了采用传统三次样条插值和分段约束三次样条插值进行EMD包络曲线拟合的对比，其中采用分段约束三次样条求得的包络曲线为虚线表示，实线包络为传统三次样条插值的结果。从图中可以看出，本文提出的分段约束三次样条插值算法对EMD包络的过冲及欠冲问题得到了抑制。

图2 不同三次样条插值EMD包络曲线对比Fig.2 Comparing the EMD envelope curve by different cubic spline interpolation

在插值算法改进基础上，本文同时将从另外二个角度继续对EMD包络拟合问题进行改进，以期进一步提高EMD包络拟合的精度和速度：① 基于局部均值点直接求取的思想，通过信号的极值点先计算出信号的局部均值点，再由局部均值点直接通过样条插值拟合得到局部均值，以此减少传统EMD计算局部均值时的样条插值次数；② 由于样条插值的精度取决于其插值型值点密集程度，因此通过型值点增密算法提高样条插值精度，以此减少三次样条插值包络拟合中的过冲和欠冲。

基于上述三方面的改进思路，本文提出一种基于增密约束三次样条插值的EMD改进算法(APCEMD)，其基本思想是利用定比分点法来增密样条插值型值点，利用基于约束条件的分段约束三次样条插值来抑制过冲与欠冲问题，利用IMF直接筛选方法来避免计算一次局部均值时需要使用两次三次样条插值的不足。通过上述三方面的同时改进，以此来提高标准EMD算法的包络曲线拟合精度，其算法流程见图3。

图3 EMD包络拟合改进算法流程Fig.3 Theimproved EMD envelope fitting algorithm process

APCEMD具体算法流程如下：

步骤1 找出待分解信号x(t)的所有极值点，并用(tk,xk)表示x(t)的第k个极值点，其中极大值点序列为xmax，极小值点序列为xmin，极大值点对应的时间序列为tmax，极小值点对应的时间序列为tmin；

(3)

步骤6：采用分段约束三次样条插值算法对排列后的局部均值点进行插值，将插值得到的包络曲线作为EMD中信号的局部均值曲线；

步骤7：根据EMD双阈值终止准则循环运算，得到IMF。

3 仿真分析

为了比较EMD与本文提出的分段约束三次样条插值EMD算法(CEMD)和APCEMD算法的不同，建立仿真信号如下：

x(t)=cos(2π×400t+0.5cos(2π×25t))

(4)

仿真信号的采样频率为2 kHz，信号长度为1 024。对x(t)进行APCEMD分解和EMD分解，分解结果分别如图4～图6所示。从图4中可以看出，APCEMD分解结果为3个IMF分量和1个残余项Res.，而在图5、图6中，EMD将仿真信号分解为5个IMF分量和1个残余项Res.。不难看出，在图5、图6的EMD、CEMD分解中，IMF4和IMF5事实上是没有任何物理意义的伪分量，而在图4的APCEMD分解中，其IMF分量准确分离了信号中内在的振荡模式，也就是说，其每一个IMF分量都具有物理意义。

图4 仿真信号的APCEMD分解Fig.4 The simulation signal APCEMD decomposition

图5 仿真信号的CEMD分解Fig.5 The simulation signal CEMD decomposition

图6 仿真信号的EMD分解Fig.6 The simulation signal EMD decomposition

与传统EMD算法相比，APCEMD虽然在筛选IMF过程中，增加了局部积分均值的求解和局部均值的增密计算，表面上看，APCEMD的计算比传统EMD算法复杂，但由于APCEMD通过增密局部均值点，增加了分段约束三次样条插值型值点，提高了样条插值精度。同时采用分段约束三次样条插值抑制了传统三次样条插值的过冲与欠冲问题。另外，由于APCEMD基于局部均值点直接筛选IMF方法，避免了传统EMD中计算一次局部均值时需要使用2次三次样条插值的不足，在某种程度上也减小1次三次样条插值的过冲与欠冲问题。因此，APCEMD每次求得的局部均值曲线都更接近于理想曲线，这样APCEMD只需要较少次“筛”的过程就可以形成IMF，其分解效率高。相反，CEMD虽然改进了包络插值算法，但由于分段后插值型值点密度不够等原因，仍出现了分解误差。相比之下，传统EMD在“筛”的过程中存在过冲与欠冲问题，需要多次“筛”的过程来形成IMF，甚至还会出现多个伪分量“筛”的过程，而在EMD分解过程中，最消耗时间的就是每“筛”出一个IMF都要进行若干次三次样条插值，因而其分解效率低。计算表明：仿真信号的EMD分解耗时1.89 s，CEMD分解耗时1.76 s，而APCEMD分解耗时1.16 s，EMD运算速度提高了近1倍。因此APCEMD的分解误差更小一些，其IMF物理意义更明确一些，运算速度也更快一些。

为了进一步分析APCEMD、CEMD、EMD的不同，对仿真信号经APCEMD、CEMD与EMD分解得到的IMF进行HHT时频分析。图7给出了仿真信号APCEMD分解的HHT时频图，从图7中可以清楚地看出400 Hz载波频率成分和25 Hz的调制频率成分，且载波频率成分和调制频率成分表现出的波动性也被反映出来。

图7 仿真信号APCEMD分解的HHT时频图Fig.7 HHT time-frequency diagram of the simulation signal APCEMD decomposition

图8为仿真信号CEMD分解的HHT时频图，图9为仿真信号EMD分解的HHT时频图，从图8、图9中可以清楚地看出25 Hz的调制频率成分，但400 Hz载波频率成分却表现出来为多种信息的叠加结果，真正的400 Hz载波频率成分被“伪频率分量”所掩盖，其中图9的“伪频率分量”显得幅度更大些。以上仿真对比分析说明，APCEMD方法相比CEMD、EMD方法来说，其分解出的分量更准确，误差更小。

图8 仿真信号CEMD分解的HHT时频图Fig.8 HHT time-frequency diagram of the simulation signal CEMD decomposition

图9 仿真信号EMD分解的HHT时频图Fig.9 HHT time-frequency diagram of the simulation signal EMD decomposition

4 转子油膜振荡故障诊断实例

图10为转子试验台测试系统图。在转子实验台上模拟油膜振荡故障，用电涡流位移传感器测量转子的径向振动量，进行整周期等相位采样，采样频率为2 kHz。图11为转子油膜振荡故障振动信号的时域波形，采样点数为1 024。这里用EMD方法与CEMD方法、APCEMD方法分别对转子油膜振荡故障振动信号进行分析，图12为应用EMD方法的分解结果，图13为应用CEMD方法的分解结果，图14为应用本文提出的APCEMD分解结果。

图10 转子试验台测试系统图Fig.10 Test system of rotor experiment platform

图11 转子油膜振荡故障振动信号时域波形Fig.11 Rotor oil film oscillation fault vibration signal time domain waveform

图12 转子故障振动信号标准EMD分解结果Fig.12 Rotor fault vibration signal standard EMD decomposition results

图13 转子故障振动信号CEMD分解结果Fig.13 Rotor fault vibration signal CEMD decomposition results

图14 转子油膜振荡故障振动信号APCEMD分解结果Fig.14 Rotor fault vibration signal standard APCEMD decomposition results

在图12中，EMD由上至下依次为分解得到的8个IMF分量和1个残余量，在图13、图14中，经过CEMD方法、APCEMD方法的信号分解结果中，IMF的个数明显少于标准EMD。为更准确地反映出图12～图14中的IMF特征信息，对上述油膜振荡故障振动信号经APCEMD、CEMD与EMD分解得到的IMF进行HHT时频分析。图15为用EMD方法得到的HHT时频分布图，其HHT时频图虽然基本能看出1倍频和0.5倍频振动模态，但1倍频的频率调制现象没有反映出来。另外，在HHT时频图中还出现了严重的伪时频谱现象，这是由于EMD三次样条包络拟合导致的分解偏差，致使EMD分解的结果中出现了许多伪分量造成的，这从图14中1倍频附近出现的伪频率分量可以看出。

图15 基于EMD方法的HHT时频图Fig.15 The HHT time-frequency diagrams based on the EMD method

图16为用CEMD方法得到的HHT时频分布图，其HHT时频图比图13中的HHT时频分布图更能清楚地看出1倍频和0.5倍频振动模态，但1倍频的频率调制被伪时频谱所干扰，调制现象不理想，且0.5倍频振动模态也太不清晰，这说明CEMD虽能改进EMD分解结果，但包络拟合的插值精度仍不够导致的分解偏差，致使EMD分解的结果中出现了许多伪分量造成的。

图16 基于EMD方法的HHT时频图Fig.16 The HHT time-frequency diagrams based on the EMD method

图17为APCEMD方法得到的HHT时频分布图。从图17中可以看出：HHT时频图上主要频率分布为1倍频及0.5倍频成分，且0.5倍频振动模态振幅较大，1倍频显示出明显的频率调制现象，这说明本文所提方法得到的HHT时频分析结果能较理想地提取转子油膜振荡的故障特征。

图17 基于APCEMD方法的HHT时频图Fig.17 The HHT time-frequency diagrams based on the APCEMD method

分析产生上述结果的原因主要如下：标准EMD采用三次样条插值进行包络曲线拟合时，易产生过冲和欠冲现象，破坏信号的模态结构，造成较大误差。CEMD虽然采用的是分段约束三次样条插值算法，改进了1次三次样条的过冲及欠冲问题，但由于插值过程中的插值点不够会造成部分插值误差，从而引起EMD的过分解，导致产生虚假的频率成分，影响瞬时频率的计算，因此CEMD中的包络拟合仍不完善。本文提出的APCEMD方法，较现有包络拟合方法的优势在于，改进了包络平均法中的插值函数，提出了一种基于约束条件的分段约束三次样条插值来来抑制三次样条插值法、Akima插值法、分段幂函数插值法、高次样条插值法存在的问题，同时利用定比分点法来增密样条插值型值点，改进了局部均值直接求取法中的插值型点问题，减小了局部均值直接求取法的插值误差。同时还采用IMF直接筛选方法来避免计算一次局部均值时需要使用两次三次样条插值的不足，吸取了局部均值直接求取法的优势，因此，通过上述三个方面的同时改进，APCEMD具有较大的优势。

但从图14和图17中可以看出，APCEMD相对于标准EMD来说，尽管通过解决EMD包络拟合问题最终能诊断出转子油膜振荡故障，但其APCEMD分解图和HHT时频分布图仍能看出改进算法中存在部分端点效应问题和模态混叠问题，但相对于标准EMD、CEMD来说上述问题已有改善。因此，这说明APCEMD除了提高包络拟合精度和运算速度外，还能改善标准EMD中的部分端点效应问题和模态混叠问题，这说明EMD包络拟合算法的好坏，与EMD最终分解结果好坏有直接关系。当然，不考虑算法实时性的情况下，为了避免因端点效应、模态混叠等因素对EMD结果的影响，以得到更高精度的EMD分解效果，在APCEMD的应用中应尽量对信号的边界进行延拓，并对EMD的模态混叠问题进行改进。

5 结 论

(1) 分析了包络曲线拟合中三次样条插值的不足，从抑制三次样条插值的过冲、欠冲影响出发，在EMD包络曲线拟合中提出一种基于约束条件的分段约束三次样条插值，仿真表明，EMD包络的过冲及欠冲问题得到了改善。

(2) 提出一种基于局部积分均值增密分段约束三次样条插值的EMD改进算法，利用定比分点法来增密样条插值型值点，利用分段约束三次样条插值来抑制过冲与欠冲问题，利用IMF直接筛选方法来避免计算一次局部均值时需要使用两次三次样条插值的不足。通过仿真对比表明，在分解精度上均要好于EMD、CEMD分解方法。

(3) 将本文改进方法应用于转子油膜涡动的故障监测诊断中，成功地提取出转子油膜涡动的故障特征，并与EMD、CEMD算法进行了对比，结果表明，改进EMD能够更准确地提取转子油膜振荡故障特征。

[1] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis[J]. Proceedings of the Royal Society, 1998, 454(1971): 903-995.

[2] LIU Hui, TIAN Hongqi, LIANG Xifeng, et al. New wind speed forecasting approaches using fast ensemble empirical model decomposition, genetic algorithm, Mind Evolutionary Algorithm and Artificial Neural Networks[J]. Renewable Energy, 2015, 83:1066-1075.

[3] 向丹, 岑健. 基于EMD熵特征融合的滚动轴承故障诊断方法[J]. 航空动力学报, 2015, 30(5): 1149-1155.

XIANG Dan, CEN Jian. Method of roller bearing fault diagnosis based on feature fusion of EMD entropy[J]．Journal of Aerospace Power, 2015, 30(5): 1149-1155.

[4] LEI Yaguo,LIN Jing,HE Zhengjia．A review on empirical mode decomposition in fault diagnosis of rotating machinery [J]．Mechanical Systems and Signal Processing, 2013, 35(1/2):108-126．

[5] WANG Yanxue, MARKERT R, XIANG Jiawei, et al. Research on variational mode decomposition and its application in detecting rub-impact fault of the rotor system[J]. Mechanical Systems and Signal Processing, 2015,60/61:243-251.

[6] QIN S R, ZHONG Y M.A new envelope algorithm of Hilbert-Huang transform[J]. Mechanical Systems and Signal Processing, 2006, 20: 1941-1952.

[7] 钟佑明. 希尔伯特-黄变换局瞬信号分析理论的研究[D]. 重庆:重庆大学, 2002.

[8] 杨世锡, 胡劲松, 吴昭同, 等. 基于高次样条插值的经验模态分解方法研究[J]. 浙江大学学报(工学版), 2004, 38(3): 267-270.

YANG Shixi, HU Jinsong, WU Zhaotong, et al. Study of empiricalmode decomposition based on high-order spline interpolation[J]. Journal of Zhejiang University(Engineering Science), 2004, 38(3): 267-270.

[9] CHEN Q H，HUANG N E，Riemenschneider S, et al．A B-spline approach for empirical mode decompositions[J]. Advances in Computational Mathematics, 2006, 24:171-195.

[10] 余泊. 自适应时频分析方法及其在故障诊断中的应用研究[D]. 大连: 大连理工大学, 1998.

[11] 盖强, 马孝江, 张海勇,等. 几种局域波分解方法的比较研究[J]. 系统工程与电子技术, 2002, 24(2): 57-59.

GAI Qiang, MA Xiaojiang, ZHANG Haiyong, et al. Comparing study for decomposing methods using local wave method[J]. Systems Engineering and Electronics,2002,24(2): 57-59.

[12] ZHONG Y M, QIN S R, TANG B P. Research on theoretic evidence and realization of directly-mean EMD method[J].Chinese Journal of Mechanica Engineering, 2004, 17(3):399-404.

[13] 侯云山, 刘宏兵. 一种约束三次样条函数的构造方法及应用[J]. 信阳师范学院学报(自然科学版), 2004, 17(4): 403-404.

HOU Yunshan, LIU Hongbing. The construction and application of a type of constrained cubic spline[J]. Journal of Xinyang Teachers College(Natural Science Edition), 2004, 17(4): 403-404.

[14] 陈一贤. HHT方法分析[D]. 杭州：浙江大学, 2007.

[15] HONG H，WANG X L,TAO Z Y. Local intergal mean-based sifting for empirical mode decomposition[J]. IEEE Signal Processing Letters, 2009, 16(10): 841-844.

A improved EMD algorithm based on local integral average constraint cubic spline interpolation by adding dense points and its application in rotary machine fault diagnosis

CAI Yanping1,3, XU Guanghua1,2, LI Aihua3, WANG Tao3

(1. School of Mechanical Engineering，Xi’an Jiaotong University，Xi’an 710049, China;2. State Key Laboratory for Manufacturing Systems Engineering，Xi’an Jiaotong University，Xi’an 710049, China;3. Department No.5, Rocket Force University of Engineering, Xi’an 710025, China)

Envelope fit is a key link of empirical mode decomposition (EMD). The current methods for solving EMD envelope fit problems have two kinds including the envelope average method and the local mean direct strike method. But, these two methods have cubic spline interpolation overshoot and undershoot, screening speed or local mean estimate errors and other problems. For envelope fit problems of EMD, an improved EMD algorithm was proposed here, based on the constraint cubic spline interpolation with local integral average and adding dense points. The constant proportion division point method was used to add dense points of EMD spline interpolation, the piecewise constrained cubic spline interpolation was employed to suppress EMD envelope overshoot and undershoot problems, and the IMF (intrinsic mode function) direct screening method was adopted to reduce spline interpolation times when calculating the EMD local average. Based on the implement of three above aspects to improve EMD, the proposed method can improved the standard EMD algorithm calculation speed and the envelope curve fitting accuracy. The simulated signals and rotor fault vibration signals analyses showed that the proposed method possesses better advantages both in fitting error and operation speed, the EMD decomposition accuracy and decomposition efficiency increase. Lastly, the proposed method was successfully applied in the fault diagnosis of rotary machines.

rotating machinery; fault diagnosis; empirical mode decomposition (EMD); spline interpolation; envelope fit

国家自然科学基金(51405498);陕西省自然科学基金(2013JQ8023);中国博士后基金(2015M582642)

2015-06-04 修改稿收到日期：2015-08-06

蔡艳平 男，博士后，硕士生导师，1982年生

TN911.7；TH165.3

A

10.13465/j.cnki.jvs.2016.19.014