1、收稿日期: 修订日期: (日期待编辑部填写) 基金项目: 国家优秀青年科学基金(51222503)、教育部新世纪优秀人才支持计划(NCET-11-0421)和陕西省自然科学基础研究计划项目资助。一种改进的 EEMD 方法及其应用研究孔德同 1，范 炜 1，雷亚国 2，丁小川 1，王 志 1（1. 华电电力科学研究院，浙江 杭州 310030；2. 西安交通大学机械制造系统工程国家重点实验室，陕西 西安 710049）摘要：针对集成经验模式分解(Ensemble empirical mode decomposition, EEMD)中协助噪声幅值大小需要人为经验确定的不足，基于经验模式分解(Em
2、pirical mode decomposition, EMD)二进滤波器特性，讨论了 EMD 出现模式混淆的原因，研究了 EEMD 中协助噪声幅值大小的确定原则，提出基于极值点分布特性的改进 EEMD 方法，通过遍态历经，以极值点分布特性为评价参数，自适应确定EEMD 方法中高斯白噪声优化幅值。通过数据仿真，验证了其有效性。最后，应用于转子早期故障诊断中，结果显示可以自适应确定噪声幅值，避免参数人为选择导致分解结果的盲目性，有效抑制了传统 EMD 方法的模式混淆现象，可有效识别转子早期碰磨引起的故障特征。关键词：改进 EEMD；极值点分布特性；故障诊断中图分类号：TH17 文献标识码： A引
3、言EMD是美国工程院N.HUANG博士于1998年提出的一种自适应数据处理方法，在非线性、非平稳信号分析中具有显著优势。与FFT、小波分解等不同，EMD方法无需选择基函数，其分解完全基于信号本身极值点分布，通过多次筛选，将信号分解为多个表征信号中某种单一模态的本征模式分量(Instrinsic mode function, IMF)与一个趋势项，得到国内外广泛关注 2。然而，当信号的极值点分布不均时，EMD 分解结果会出现“过冲 ”、 “欠冲”现象，导致模式混淆 3。针对以上问题，Zhao等 42009年提出 EEMD方法，通过对信号加入高斯白噪声，改善信号极值点分布，减小模式混淆。EEMD方
4、法已被成功应用于转子 5、轴承 6、电机 7等机械设备的故障诊断中。2009年陈略等 8指出EEMD 加入高斯白噪声的幅值不能改变原始信号高频成分极值点分布，然而文章提出的高频成分依赖于EMD分解结果，当 EMD分解得到的第一个IMF存在模式混淆时，导致EEMD 加入噪声不准确。2010年Jian Zhang等 9研究了EEMD方法中加入的高斯白噪声幅值及总体平均次数两个参数，从加入的高斯白噪声与原始信号能量比的角度规定了加入噪声的原则，但是该方法仅仅考虑了由两种成分组成的信号，对于存在多种模式分量的信号未展开研究。2010年，雷亚国通过改进Hilbert-Huang变换，提出敏感IMF的选择
5、方法，成功诊断出转子早期碰摩故障 10，2013年雷亚国等 11以有色噪声代替高斯白噪声，有效改善了信号的极值点分布，但该研究工作并未建立加入噪声大小的准则。针对EEMD方法中加入的高斯白噪声幅值大小问题，本文提出基于极值点分布特性的改进EEMD方法。研究了EEMD 加入噪声的准则，利用全局寻优，以极值点分布特性为评价函数，建立了EEMD方法中信号极值点分布特性与加入噪声大小的对应关系，评估不同噪声大小对原始信号极值点的改善程度，可实现噪声幅值的自适应优化选取，消除参数人为选择导致分解结果的盲目性与主观性。1 基本原理1.1 EEMD 方法基本理论EMD/EEMD方法本质上是基于极值点的筛选过
6、程。以仿真信号为例，EMD分解时，如图 1所示，利用信号的极大值与极小值，通过三次样条拟合出极大值包络与极小值包络，得到局部均值。将信号减去局部均值，重复以上过程，直到筛选出的成分满足IMF条件。最后，信号依次被分解为若干个频率由高到低的IMF与一个趋势项，如图2所示。当极值点分布不均时，通过插值得到的局部均值发生扭曲，导致筛选结果出现模式混淆现象。因此，极值点分布特性是EMD 分解结果的决定性因素。针对极值点分布问题，EEMD对信号加入高斯白噪声，利用白噪声的频率均布特性，改善信号极值点分布，通过多次分解取平均，减小加入噪声对分解结果的影响，得到EEMD分解的最终结果。其算法流程如图3所示。
7、时间 t/s幅值 A/g0 0.2 0.4 0.6 0.8 1-20-1001020304050上上上上上上上上上上上上图 1 EMD筛选示意图Fig.1 EMD filter schematic-20020-202-101-100100 0.2 0.4 0.6 0.8 101020时间 t/s幅值 A/gIMF1IMF2IMF3IMF4R图 2 EMD筛选结果Fig.2 EMD sifting results开始输入信号 x(t)初始化噪声大小 e、总体平均次数 M第 1 次分解 第 2 次分解 第 M 次分解EMD 分解 EMD 分解 EMD 分解加入高斯白噪声 x2(t)= x(t)+n
8、2(t)加入高斯白噪声xM(t)= x(t)+nM(t)加入高斯白噪声x1(t)= x(t)+n1(t)结束得到一系列IMF： c1,i 得到一系列IMF： c2,i得到一系列IMF：c M,iMmii1,计算总体平均图 3 EEMD算法流程图Fig.3 Flowchart of EMD algorithm 1.2 EEMD 方法加入噪声大小准则当信号不连续时，极值点分布不均，EMD 分解结果会出现模式混淆现象，具体表现为：不同尺度的振动信号被分解在同一个IMF中；同一尺度的信号被分解在不同IMF 中。EEMD对信号加入高斯白噪声，改善信号的极值点分布，通过借助高斯白噪声的频率均布特性，在频域
9、为信号构建频率均布的尺度，信号中的相应成分自动映射到高斯白噪声建立的不同尺度中，克服模式混淆。然而，加入的噪声幅值大小需要人为经验决定。噪声幅值优化问题一直是影响EEMD的瓶颈。当噪声幅值较小时，由于噪声振动微弱，难以改善信号高频成分的极值点分布，导致EEMD分解结果与EMD 类似，无法克服模式混淆不足。另一方面，为了避免同一尺度的振动信号被分解在不同IMF中，应避免加入的高斯白噪声幅值较大，噪声幅值较大时，噪声振动剧烈，剧烈振动的噪声破坏了低频成分的极值点分布，噪声协助信号x m(t)的极值点均由高斯白噪声诱发，导致EEMD分解结果会产生多余成分，出现模式混淆。尽管Zhao等人 4指出按照下
10、式确定加入的高斯白噪声大小e。或者 式1N0ln2Ne其中为加入加入幅值为e 的高斯白噪声平均 N次后分解误差。通过大量仿真验证，一般建议e取0.2然而，该方法并未考虑加入高斯白噪声后信号的极值点分布特性，在实际应用中往往需要根据分解结果多次调整加入的高斯白噪声大小再次分解。因此，加入的高斯白噪声大小应在未破坏原始信号低频成分极值点分布特性基础上尽可能改善高频成分极值点分布。2 基于极值点分布特性的改进 EEMD 方法通过以上分析，不难看出，EEMD方法加入噪声后应保证极值点在各尺度分布最具离散性，进而增强各尺度的可区分性，减小各尺度间的模式混淆。因此，若以极值点分布特性为评价参数，通过测试加
11、入不同大小高斯白噪声后的极值点分布情况，遍态历经，可以自适应确定高斯白噪声优化幅值，保证在未破坏原始信号低频成分极值点分布特性基础上，使得信号高频成分的极值点分布最为均匀。根据以上思路，本文提出基于极值点分布特性的改进EEMD方法，其实现过程如下：1）初始确定加入高斯白噪声的最大幅值e及网格数I，确定网格i对应的噪声大小e i=i/I*e，其中i=1,2,I。如图4所示。针对最大幅值e及网格数I ，e过小，可能导致高斯白噪声优化幅值没有落在寻优范围内；e过大，加入的高斯白噪声会破坏高频成分的极值点分布，导致分解结果失真；I过小，相邻网格对应的噪声幅值相差较大，寻优结果与高斯白噪声优化幅值存在较
12、大差异；I过大，相邻网格对应的噪声幅值相差很小，对应的EEMD分解结果差异性很小。通过上述分析及大量数据仿真结果表明，一般e取0.5，I 取1000即可满足要求。网格编号 i0 10000.5噪声幅值eI 0 1000噪声幅值 eiFi幅值e优化幅值 ey图 4 改进EEMD算法的噪声寻优方法Fig.4 Optimization method for the improved EEMD algorithm2）对信号x(t) 依次加入幅值大小为e i的高斯白噪声ni(t)，得到 xi(t)=x(t)+ni(t)，遍态历经，得到噪声协助信号x i(t)的极大值与极小值对应的位置Max(j i)、M
13、in(j i)，以信号x i(t)极值点分布特性F i为评价参数： 22(1()(1iiiii knDjajD(ji=1,2,Ji，k i=1,2,Ki)其中，D代表标准差，J i为信号x i(t)的极大值点数，Ki为信号 xi(t)的极小值点数；一般来说，随着噪声e i幅值的不断增大，原始信号的高频成分极值点分布逐步被改善，F i逐渐增大（若原始信号的极值点分布较为均匀，则F i基本保持不变） ；当噪声e i幅值增大到一定程度时，加入的高斯白噪声破坏了高频成分的极值点分布，导致信号极值点均由高斯白噪声诱发，F i迅速减小，直至趋于稳定。因此， Fi迅速减小时所对应的高斯白噪声e i即为优化幅
14、值e y。3） 、按照加入噪声大小准则，自适应确定高斯白噪声优化幅值e y，EEMD分解。图5为改进EEMD方法的流程图。初始化噪声最大值 e， 网格数 I对信号 x(t)依次加入幅值为 ei 的噪声 ni(t)得到xi(t)= x(t)+ ei(t)确定网格 i 对应的噪声大小 ei=i/I*e输入信号 x(t)确定 xi(t)所有极大值点与极小值所在位置 Max(ji)、Min(ji)计算 Fi得到信号的极值点分布与噪声大小的曲线自适应确定高斯白噪声优化幅值 ey ，EEMD 分解结束图 5 改进EEMD方法流程图Fig.5 Flowchart of improved EEMD metho
15、d 3 仿真与工程应用3.1 仿真验证旋转机械正常工作时一般呈现周期性的运动，用谐波分量代表设备的旋转振动信号。由于旋转机械的损伤故障通常表现为高频冲击响应特性，且在时域上不连续，呈现出冲击与调制。因此，仿真一组包括冲击c 1、调制c 2、谐波分量c 3等 3个成分的信号S表示故障信号，示于图6中。时间 t/s幅值A/g-101-303-5050 0.1 0.2 0.3 0.4 0.5-808c1c2c3S图 6 仿真信号Fig.6 Simulated signal利用本文提出的基于极值点分布特性的改进EEMD方法对信号S进行分解。初始化噪声最大值e=0.5，网格数1000，得到信号的极值点分
16、布与噪声大小的曲线如图7。0 0.1 0.2 0.3 0.4 0.505101520噪声大小 A幅值 A/g高斯白噪声优化幅值 0.008图 7 仿真信号极值点与噪声大小的对应关系Fig.7 Relationship between the extreme points of the simulated signal and noise level从图7中可以看出，当加入的噪声较小时(e i0.008)，极值点分布的不均匀性随噪声的增加而减小，最后趋于稳定，主要是由于随着噪声较大，振动剧烈，剧烈振动的噪声破坏了高频成分的振动特性，噪声协助信号的极值点全部是由加入噪声诱发。根据极值点分布特性与噪
17、声大小的对应曲线可以确定，最优噪声大小为0.008，将确定大小的高斯白噪声加入仿真信号s进行EEMD分解，分解得到的IMF 如图8。-101-303-5050 0.1 0.2 0.3 0.4 0.5-505时间 t/s幅值 A/gIMF1IMF 2IMF 3R图 8 仿真信号的改进EEMD分解结果Fig.8 Improved EEMD decomposition results of the simulated signal可以看出，分解得到的前三个 IMF 分别是冲击(IMF1)、调制(IMF 2)、谐波分量(IMF 3)，依次对应原始信号中的各组成成分 c1、c 2、 c3，每个 IMF
18、的幅值与构成原始信号的各成分几乎相等，各 IMF 之间不存在模式混淆现象。为了进行对比，分别利用EMD 方法与传统EEMD方法对该仿真信号也进行分解，分解结果如图9、10所示。-404-20020-100100 0.1 0.2 0.3 0.4 0.5-606时间 t/s幅值 A/gIMF1IMF2IMF3R图 9 仿真信号的EMD分解结果Fig.9 EMD decomposition results of the simulated signal时间 t/s-0.500.5-0.500.5-3030 0.1 0.2 0.3 0.4 0.5-606幅值A/gRIMF3IMF2IMF1图 10 仿
19、真信号的传统EEMD分解结果Fig.10 Traditional EEMD decomposition results of the simulated signal从图9中不难看出，EMD分解得到的 IMF1既包括了冲击成分，也含有部分调制成分；IMF 3相比原始信号的谐波成分c 3，存在严重失真。从图10中不难看出，传统EEMD 分解得到的IMF 1为噪声成分，IMF 2相比原始信号的冲击成分c 1，幅值衰减严重。相比改进EEMD分解结果，EMD及传统EEMD分解结果均出现了不同程度的模式混淆，得到的各IMF丧失了其固有意义。3.2 实际信号验证某炼油厂重催机组由电机、齿轮箱、压缩机、烟机
20、等组成。机组工作时，电机驱动，通过齿轮箱减速，带动压缩机及烟机工作。在线监测系统对烟机、风机和齿轮轴进行实时监测，烟机工作转速为5859r/min，即97.65Hz，齿轮箱低速轴转频为25.35Hz，采样频率2000Hz。机器工作时，监测系统显示振动过大报警，传感器安装在振动最大测点处，拾取的振动信号如图11所示。时间 t/s0 0.1 0.2 0.3 0.4 0.5-50050幅值 A/g图 11 某转子时域信号Fig.11 Waveform of the rotor vibration0 200 400 600 800 100001020频率 f/Hzf2=97.65Hzf3=193.36
21、Hzf1=25.35Hz幅值 A/g图 12 某转子信号频谱Fig.12 Spectrum of the rotor vibration从时域波形中可以看出信号中的主要成分是周期大约为0.01s的谐波分量，约为烟机的转动周期，从时域波形中无法判别引起振动过大的原因。其频谱如图12所示，频谱中的三个主要频率为f1=25.35Hz，f 2=97.65Hz与f 3=193.36Hz，分别对应齿轮箱低速轴的转频、机组转频和二次谐波。因此从信号的时域波形及频谱中均无法判别振动过大的原因。利用本文提出的改进EEMD方法对该信号进行分析。初始化噪声最大值e=0.5，网格数1000，得到信号的极值点分布与噪声
22、大小的曲线如图13。从图中可以看出，当噪声幅值较小时(e i0.02)，随着加入噪声的不断增加，极值点的分布不均匀性逐渐减少。在噪声大小为0.02时，信号极值点分布不均匀特性趋于平稳，之后，随着噪声的不断增大，信号极值点分布不均匀性急剧下降，最后趋于平稳，这是由于信号的极值点特性已被加入的噪声破坏。因此，通过利用本文提出的方法确定最佳噪声大小为0.02.0 0.1 0.2 0.3 0.4 0.5246810噪声大小 A幅值A/g高斯白噪声优化幅值 0.02图 13 实际信号极值点与噪声大小的对应关系Fig.13 Relationship between the extreme points o
23、f the actual signal and noise level将确定大小的高斯白噪声加入实际信号进行EEMD分解，分解得到的前三个IMF及趋势项分别为IMF1、IMF 2、IMF 3及R，如图14所示。-10010-40040-5050 0.1 0.2 0.3 0.4 0.5-6010时间 t/s幅值 A/gIMF1IMF2IMF3RT=0.031s图 14 实际信号的改进EEMD分解结果Fig.14 Improved EEMD decomposition results of the actual signal为了更清晰描述IMF 1周期性成分，对其包络解调分析，包络解调谱如图15所
24、示。包络谱中存在明显的频率为33Hz的谱线及其倍频，对应机组1/3的转频（97.65Hz/3=32.55Hz） 。通过以上分析表明，本文提出的改进EEMD方法分解结果IMF 1蕴含丰富故障信息，冲击信号频率是1/3 倍的机组工频，表明轴与轴瓦之间产生早期碰摩故障 12。幅值A/g0 200 400 600 800 100000.20.40.60.811.21.4频率 f/Hzf1=33Hzf2=65Hz图 15 实际信号的改进EEMD分解结果IMF 1解调谱Fig.15 IMF1 demodulation spectrum of improved EEMD 为了进行对比，分别利用EMD 方法、
25、传统EEMD方法对该信号进行分解，分解结果如图16、17所示。相比本文提出的改进EEMD分解结果， EMD分解得到的IMF 1中周期性冲击成分不明显，同时IMF 2、IMF 3存在模式混淆。传统EEMD方法分解得到的 IMF1中周期性冲击成分被背景噪声淹没，同一尺度信号被分解在了IMF 2、IMF 3两个基本模式分量中，分解结果出现严重的模式混淆。EMD方法、传统 EEMD方法分解结果无法诊断故障原因。-20020-40040-100100 0.1 0.2 0.3 0.4 0.5-804时间 t/s幅值 A/gIMF1IMF2IMF3R图 16 实际信号的EMD分解结果Fig.16 EMD d
26、ecomposition results of the actual signal-505-30030-200200 0.1 0.2 0.3 0.4 0.5-10-505时间 t/sRIMF3IMF2IMF1幅值A/g图 17 实际信号的传统EEMD分解结果Fig.17 Traditional EEMD decomposition results of the actual signal4 结论EEMD 是针对 EMD 存在模式混淆而提出的一种非线性、非平稳信号处理方法。但是，其分解结果较大程度的依赖于人为选择的噪声大小。为减小参数人为选择带来的盲目性，提出改进 EEMD 方法。通过研究协助噪
27、声大小对极值点的影响，建立了噪声协助大小的选择原则，以极值点分布为评价参数，通过全局寻优，自适应确定 EEMD 方法中高斯白噪声大小。通过仿真验证了其有效性。最后将该方法应用于转子振动信号分析，结果显示该方法可以有效诊断转子早期碰摩故障。参考文献：1HUANG Norden E, SHEN Zheng, LONG Steven R., et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysisJ. Proceeding of the Royal Society: A Mathematical Physical Distribution of signal extreme points; Fault diagnosis作者简介: 孔德同(1988 ),男,工程师。电话:(0571)85246726; E-mail: dt_
