改进EEMD方法及混沌降噪应用研究

位秀雷， 林瑞霖， 刘树勇， 杨庆超

(1. 中国人民解放军91404部队，河北 秦皇岛 066000; 2. 海军工程大学 动力工程学院，武汉 430033)

改进EEMD方法及混沌降噪应用研究

位秀雷1， 林瑞霖2， 刘树勇2， 杨庆超2

(1. 中国人民解放军91404部队，河北 秦皇岛 066000; 2. 海军工程大学 动力工程学院，武汉 430033)

在总体平均经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)降噪过程中，对本征模态分量(Intrinsic Mode Function，IMF)的有效处理一直是影响降噪效果的关键。为此，提出一种基于改进EEMD的去噪方法。基于“3σ”法则和奇异值分解(Singular Value Decomposition，SVD)提取第一个IMF分量中有用信号细节。利用连续均方误差准则对剩余IMF分量进行高低频区分，分别使用SVD和S-G算法提取高低频分量的有用信号，可以有效避免了高频部分有用信号的流失，同时剔除低频分量中的部分噪声，克服了EEMD去噪时IMFs难以有效处理的不足。为了验证该方法的有效性，进行了数字仿真与双势阱混沌振动试验，结果表明，该方法的降噪效果优于小波加权和EEMD去噪方法。

总体平均经验模态分解；混沌信号；奇异值分解；降噪；S-G滤波

混沌振动信号在采集过程中不可避免受到复杂噪声干扰，极大干扰了信号真实信息的解读，从而影响了混沌振动信号控制、识别、预测等进一步分析，并且混沌复杂的动力学行为使其具有功率谱宽带性和似噪声性，其频带与叠加噪声的频带往往全部或部分重叠，传统的线性去噪方法很难实现有效滤波[1]。

Huang等[2]提出了一种新的信号处理方法-EMD(Empirical Mode Decompomposition)方法。与小波变换方法相比，EMD无需信号的先验知识，其分解完全依赖信号本身，数据分解真实可靠，因此被广泛应用于机械振动信号分析[3-6]。但是，在脉冲强干扰的影响下，EMD分解出来的本征模态分量(Intrinsic Mode Function，IMF)会发生畸变，导致信号失真[7]，且EMD本身存在一些不足，如模式混叠、端点效应、停止条件等[8]。为了抑制模式混叠，Wu等[9]提出了总体平均经验模态分解方法，有效地克服了这一缺陷。在IMF处理过程中，文献[10]将EEMD 方法用于对疲劳应变信号降噪，其方法是对信号做EEMD 分解后，选取IMF 分量来重构信号以对疲劳信号降噪，但文中也未说明IMF分量的选取方法；陈仁祥等基于能量密度与其平均周期的乘积为一常量这一特点，设计了自动选择IMF 分量重构信号的算法。但是，EEMD的时空域算法因简单地去掉1个或多个IMF分量实现去噪，会导致相应分量上的有效信号一起被剔除，进而导致信号失真，为了在分离噪声的过程中尽可能地保留有用信号，本文提出了改进EEMD降噪方法。首先对信号进行EEMD分解，分析各IMF分量特性，然后分别使用“3σ”法则、奇异值分解[11]和S-G滤波方法对IMF分量进行处理，克服了EEMD去噪时IMFS选择的难题。为了验证本文方法的有效性，进行了数字仿真与双势阱混沌振动试验，结果表明，本文方法优于小波加权和传统的EEMD去噪方法。

1 EEMD原理

利用EMD进行信号处理时，由于异常事件干扰导致极值点分布不均匀产生了模式混叠现象。为此，Wu等将白噪声加入待分解信号来抑制异常事件，利用白噪声频谱的均匀分布来使不同尺度的信号自动分布到合适的参考尺度上。同时，利用白噪声的零均值特性，经过多次平均使噪声相互抵消，从而抑制甚至完全消除噪声的影响。EEMD的本质就是叠加高斯白噪声的多次经验模式分解，其步骤如下：

(1) 在原始信号x(t)中叠加均值为0，幅值和标准差为常数的高斯白噪声ni(t)，i=1～M，叠加次数为M(M>1)，即：

xi(t)=x(t)+ni(t)

(1)

(2) 对xi(t)进行EMD分解，得到N个IMF记为aij(t)，j=1～N，余项表示为ri(t)。其中aij(t)表示第i次叠加高斯白噪声后，分解得到的第j个IMF分量。

(3) 由于不相关随机序列的统计均值为0，所以将以上步骤所得的IMF进行平均运算，即可消除多次叠加高斯白噪声对真实IMF的影响，平均后得到的IMF为

(2)

式中，aj(t)表示对原始信号进行EEMD分解后所得的第j个IMF分量。

2 阈值决策方法

2.1利用“3σ”法则和SVD提取IMF1中的有用信号

信号经EEMD去噪处理，通常认为第一层IMF分量全部由噪声构成，但随着学者深入研究发现，IMF1中仍含有一定量的信号成分[12]。提取IMF1中的信号细节成分，会提高去噪效果，减少信号失真。但由于缺乏先验知识，因此，对IMF1处理仍是个难题。在IMF1中，噪声占绝大部分，仅有少量的信号细节成分，且所含噪声仍近似服从零均值正态分布，文献[12]利用“3σ”法则进行细节信息提取。

(3)

由式(3)可以看出，“3σ”法则提取的细节分量是一种窄时间段的高幅值冲击信号，图1是利用“3σ”法则对信噪比为5 dB的Lorenz时间序列经EEMD分解后IMF1提取的细节信息，显然，直接把其作为IMF1提取的有用信号是不合理的，并且经大量仿真实验证明：信号含噪声强度越大，“3σ”法则提取的细节分量就越多；反之，越少，甚至当噪声强度小至系数幅值都小于3σ，就提取不到任何细节分量，图2是利用“3σ”法则对信噪比为15 dB的Lorenz时间序列经EEMD分解后IMF1提取的细节信息。因此需分两种情况对进行处理，第一是“3σ”法则提取的细节分量不为零，即对象是噪声强度大的混沌信号；第二是是“3σ”法则提取的细节分量为零，即对象噪声强度小的混沌信号。

图1 SNR为5 dB的IMF1细节

图2 SNR为15 dB的IMF1细节

对于第二种情况，噪声强度较小，利用传统的非线性去噪方法即可取得较好的降噪效果，本文不再详述。第一种情况，即低信噪比、高强度的混沌含噪信号，先利用“3σ”法则提取IMF1的细节分量，再采用奇异值分解(SVD)方法对有用信号进行提取，SVD降噪的关键在于选择合适的奇异值进行重构，由于“3σ”法则提取IMF1的细节分量可以作为IMF1所含有用信号估计的先验知识，因此可以利用有用信号能量相近原则选取SVD重构的奇异值个数。具体的步骤如下：

1) 利用式(3)对进行IMF1细节提取；

2) 利用下式估计IMF1所含有用信号的能量；

(4)

式中：M为满足“3σ”法则条件的本征模态分量个数。

3) 对IMF1进行奇异值分解，假设分解后奇异值为λ=[λ1,λ2,…,λn]，若选择前k个特征值重构得到降噪后的信号S=[S1,S2,…,Sk]，则根据下式便可求出k的值：

(5)

2.2SVD提取高频分量中的有用信号

由于含噪信号经EEMD分解，IMF1相当于一个高通滤波器，其频谱占据了原始信号所含频谱的一半，白噪声的各阶IMF的能量密度按照2倍关系逐渐递减，即随着IMF阶数的增加白噪声的比重大幅度减少，而有用信号的比重大幅度增加，因此对于高频IMF，IMF1经奇异值分解，各个特征量比较平均，而其他IMF经奇异值分解后，信号的特征值会明显大于噪声的特征值，依据这一特性，便可利用差分谱的方法选择合理的奇异值个数进行重构，达到去噪目的。

首先，需要对IMF分量进行高低频区分，噪声主要分布在高频IMF分量上，而信号主要分布在低频IMF分量上，因此可以利用连续均方误差准则[13](Consecutive Mean Square Error，CMSE)对信号分量起主导作用模态与噪声起主导作用模态进行区分，即找到一个索引值js，使得从该索引开始往后的IMF分量对信号进行重构的误差最小。

(6)

式中：N为信号长度；n为IMF分量的个数；IMFk(ti)表示第k个IMF的第ti个分量的重构误差，基于该准则，索引值js可由下式给出：

(7)

式中，arg min表示重构误差取最小的函数。

2.3S-G滤波方法平滑低频分量

由于低频IMF分量中大部分都是有用信号，噪声成分比较少，所以奇异值分解后特征值分布比较平均，不适合利用SVD提取有用信号。赵志宏等指出Savitaky-Golay滤波算法对低频IMF分量有很好的降噪效果，因此，本文低频IMF分量按式(8)做SG平滑处理。设ci是其中的一个IMF分量系数，在ci附近以nl+nr+1个点在最小二乘意义下拟合一个M次多项式pi(c)，多项式pi(c)在ci的值，即光滑数值gi表示为

(8)

式中：nl为ci左边点的个数；nr为ci右边点的个数；bp为多项式的系数。设实测数据为yi，为了使pi(c)拟合测试数据，必须定义系数bp，使得式(9)达到最优。

(9)

3 阈值决策降噪仿真分析

实验信号为Lorenz方程产生的混沌信号，如图3所示。

(10)

其中参数分别为：a=10，r=34，c=8/3。

(a) 原始时间序列

(b) 信噪比SNR=5 dB的含噪序列

图4是Lorenz时间序列的EEMD结果，包括7个IMF分量(IMF1～7)和一个余项(res)，根据式(6)，(7)计算索引值js=4，即IMF1～3为高频噪声主导，IMF4～7为低频信号主导，因此，EEMD阈值决策降噪的步骤具体为：

1) 对IMF1依据“3σ”法则和SVD提取有用信号细节d1。图5为IMF1经SVD后的特征值，可以看出，特征值分布比较均匀，利用差分谱选取特征值个数无法取得理想效果，而先利用“3σ”法则估计IMF1中信号分量的能量，为SVD特征值的选取提供先验知识，则可达到较好的降噪效果。去噪后的信号d1如图6所示，相比于图1，提取的信号比较均匀，更为合理，由于提取的信号仅占第一个IMF分量的0.003，因此，信号幅度很小。

图4 Lorenz信号EEMD结果

图5 IMF1SVD的特征值

图6 去噪后的IMF1

2) 利用SVD提取高频分量的有用细节。由于分解得到的IMF分量频率按2的指数降幂排列，IMF分量频率骤降，利用连续均方误差准则确定高低频分界点js=4，即噪声起主导作用的模态为IMF2、IMF3，图7为IMF2、IMF3经SVD后的特征值，可以看出，前面几个特征值比较大，后面的比较均匀，利用差分谱选取特征值个数即可取得理想效果，重构选定的奇异值，得到d2、d3，如图8所示。

(a) IMF2 SVD的特征值

(b) IMF3SVD的特征值

(a) 去噪后的IMF2

(b) 去噪后的IMF3

3) 对剩余IMF分量利用Savitaky-Golay滤波算法进行平滑处理。SG算法的拟合阶数取3，数据窗口为15，得到d4～d8。

4) 重构信号x_denoised(t)=d1+d2+…+d8。

为了对比降噪效果，对Lorenz含噪序列分别采用小波加权阈值去噪方法[14]、传统EEMD和所提方法进行降噪处理，小波加权阈值去噪方法基小波选取“sym8”小波，分解层数为5；传统EEMD降噪中，采用硬阈值方法，3种方法对5 dB的Lorenz含噪序列的去噪效果如图9所示。对于模型已知的Lorenz时间序列，采用信噪比SNR和均方误差MSE来评价降噪效果，SNR越大，MSE越小，降噪效果越好。为了全面检验3中方法对不同噪声强度的含噪信号的去噪效果，分别计算对SNR=0 dB、5 dB、10 dB、15 dB的Lorenz含噪序列降噪后的SNR和MSE，详见表1。

由图9可以看出，小波加权阈值法和EEMD方法去噪后的信号含有很多毛刺，整体不够平滑，后者尤为突出，而本文方法去噪后的信号较为光滑平整更接近原始信号。由表1的计算结果看出，本文所提方法整体上要优于小波加权阈值和EEMD方法，并且对较低信噪比的Lorenz含噪序列尤为突出，随信噪比逐渐增大，3种方法的差距减小，这和2.1节中的分析一致，也是由于低强度噪声对混沌时间序列的轨迹造成干扰较小，传统的非线性降噪方法都比较有效。

(a) 小波加权阈值法

(b) EEMD方法

(c) 本文所提方法

降噪方法SNR/MSE0dB5dB10dB15dB小波阈值加权法14.273/0.225716.536/0.167820.305/0.07624.578/0.032EEMD14.095/0.237516.475/0.176220.257/0.08324.557/0.037本文方法15.237/0.209717.335/0.109420.828/0.05524.732/0.025

4 机械式混沌振动信号处理

混沌振动是混沌科学研究的重要课题，在混沌隔振、混沌振动压路机等方面有着广泛的应用，国内外学者围绕混沌振动的应用、控制和识别等问题开展了深入的研究。但是，在混沌振动信号采集过程中，难免受到周围环境噪声的干扰，严重影响了关联维数、Lyapunov指数等混沌特征量的计算，为混沌参数区域确定和混沌控制的研究带来了困难。

为了进一步验证本文所提方法对于模型未知的混沌信号的去噪效果，本文基于双势阱理论的单端磁吸式混沌振动试验装置产生的振动信号为对象，双势阱单端磁吸式混沌振动装置如图10所示。

实验本质为正弦信号的慢速频率扫描实验，扫描频率范围为5～25 Hz，采样频率为2 kHz，数据采集时长为5 s。调节功率放大器增益为1、激励频率为13 Hz时，重构相空间参数中嵌入维数为3，延迟时间为10，得到的重构相图如图11(a)所示。分别使用小波阈值加权法、EEMD和本文方法对采集的振动信号进行降噪处理。

图10 双势阱单端磁吸式混沌振动装置

以及本文方法对原始信号进行降噪处理，去噪后的二维相图如图11(b)、(c)、(d)所示，通过比较可以看出，本文方法降噪后的二维相图曲线更加光滑，更能清晰地展现原信号吸引子的几何结构。由于原始信号为机械式振动信号，其信号和噪声未知，不能利用SNR和MSE进行定量比较，但是对于混沌信号，其自相关函数值比较大，且远远大于噪声的自相关函数值，因此，可以利用3种方法去噪后的自相关函数值以评价去噪效果。

(a)原始信号(b)小波加权阈值方法(c)EEMD方法(d)本文方法

图11 振动信号去噪前后二维相图

Fig.11 The phase diagrams of vibration signals

3种方法去噪前后的部分自相关函数值如表3所示，可以看出，原始信号受到噪声的干扰，其自相关函数值相比去噪后要小很多，而本文所提方法降噪后序列的自相关函数值最大，进一步展现了其降噪的优越性。

通过扫频试验，不断调节激振器激励频率(以步长0.5 Hz)和功率放大器的增益(以步长0.1 V)，利用本文方法对采集到的信号进行降噪处理，然后分别计算降噪前后的信号的Lyapunov指数(LE)和关联维数。当激励增益为1时，部分去噪前后信号特征量如表3所示，混沌运动的LE为正值且关联维数为分数维，去噪后信号的两个特征量都证实系统运动为混沌运动，而f=14 Hz、15 Hz时，原始信号的关联维数接近整数维，说明了噪声干扰了系统特征量的计算，也证实了本文方法为实现混沌振动信号进一步分析及其应用提供了重要的前提条件。通过对本文方法降噪后信号特征量计算判断振动的运动形式，并以激励频率为横坐标，激励增益为纵坐标作出系统不同响应的参数区域图，如图12所示，从图中可以明显获得混沌运动的参数区域。

表2 3种方法去噪后部分自相关函数值

表3 部分实测信号的特征指数

5 结 论

基于总体平均经验模态分解，提出一种阈值决策的去噪方法。信号经总体模态分解得到固有模态函数分量，首先根据“3σ”法则和SVD提取第一个IMF分量中有用信号细节。然后，利用连续均方误差准则对剩余IMF分量进行高低频区分，分别使用SVD和S-G算法提取高低频分量的有用信号，可以有效避免了高频部分有用信号的流失，同时剔除低频分量中的部分噪声，克服了EEMD去噪时IMFs选择的难题。最后通过对多普勒仿真信号和实测振动信号进行去噪分析，结果表明，本文方法是非常有效的。

图12 系统不同响应的参数区示意图

致谢：感谢杨庆超讲师和Southampton University Solent Institution of Acoustic: Jian Jiang James, Chris, Lee, Lawrance的讨论。

[1] 赵志宏，杨绍普，申永军. 一种改进的EMD降噪方法[J]. 振动与冲击，2009，28(12)：35-37.

ZHAO Zhihong，YANG Shaopu，SHEN Yongjun. Improved EMD based de-noising method J]. Journal of Vibration and Shock，2009，28(12)：35-37.

[2] HUANG N E，SHEN Z，LONG S R，et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society of London Series A-Mathematical Physical and Engineering Science，1998，454：903-995.

[3] MANUEL B V，WENG B W，BARNER K E. ECG signal denoising and baseline wander correction based on the empirical mode decomposition[J]. Computers in Biology and Medicine，2008，38(1)：1-13.

[4] 王红军，万鹏.基于EEMD和小波包变换的早期故障敏感特征获取[J]. 北京理工大学学报，2013，33(9)：945-950.

WANG Hongjun，WAN Peng. Sensitive features extraction of early fault based on EEMD and WPT[J].Transactions of Beijing Institute of Technology，2013，33(9)：945-950.

[5] 甘雨，隋立芬，王冰. 经验模态分解阈值消噪方法及其在惯性导航系统数据处理中的应用[J].测绘学报，2012，41(4)：504-509.

GAN Yu，SUI Lifen，WANG Bing. EMD threshold de-noising and its applications in INS data processing[J].Acta Geodaetica et Cartographica Sinica，2012，41(4)：504-509.

[6] 陈忠，符和超. 基于EEMD方法的永磁同步电机电磁噪声诊断[J]. 振动与冲击，2013，32(20)：124-128.

CHEN Zhong，FU Hechao. Electromagnetic noise diagnosis for a permanent magnet synchronous motor based on EEMD[J]. Journal of Vibration and Shock，2013，32(20)：124-128.

[7] 柏林，刘小峰，秦树人. 小波-形态-EMD综合分析法及其应用[J]. 振动与冲击，2008，27(5)：1-4.

BO Lin，LIU Xiaofeng，QIN Shuren. Hybrid wavelet-morphology-EMD analysis and its application[J]. Journal of Vibration and Shock，2008，27(5)：1-4.

[8] 陈仁祥，汤宝平，马靖华. 基于EEMD的振动信号自适应降噪方法[J]. 振动与冲击，2012，31(15)：82-86.

CHEN Renxiang，TANG Baoping，MA Jinghua. Adaptive de-noising method based on ensemble empirical mode decomposition for vibration signal[J]. Journal of Vibration and Shock，2012，31(15)：82-86.

[9] WU Z，HUANG N E. A study of the characteristics of white noise using the empirical mode decomposition method[J]. Proc. R. London A，2004，460：1597-1611.

[10] 陈隽，李想. 运用总体经验模态分解的疲劳信号降噪方法[J].振动、测试与诊断，2011，31(1)：15-19.

CHEN Juan, LI Xiang. Application of Ensemble Empirical Mode decomposition to noise reduction of fatigue signal[J].Journal of Vibration, Measurement & Diabnosis，2011，31(1)：15-19.

[11] 雷达，钟诗圣.基于奇异值分解和经验模态分解的航空发动机健康信号降噪[J].吉林大学学报(工学版)，2013，43(3)：764-770.

LEI Da，ZHONG Shisheng. Aircraft engine health signal denoising based on singular value decomposition and empirical mode decomposition methods[J]. Journal of Jilin University (Engineering and Technology Edition)，2013，43(3)：764-770.

[12] 王文波，张晓东，汪祥莉. 基于主成分分析的经验模态分解消噪方法[J]. 电子学报，2013，41(7)：1425-1430.

WANG Wenbo，ZHANG Xiaodong，WANG Xiangli.Empirical mode decomposition de-noising method based on principal component analysis[J]. Acta Electronica Sinica，2013，41(7)：1425-1430.

[13] 张文忠，周蓉，武旭红，等. 利用白噪声分解特征的EEMD阈值降噪方法[J].测绘科学技术学报，2013，30(3)：255-259.

ZHANG Wenzhong，ZHOU Rong，WU Xuhong，et al. A threshold de-noising method based on the characteristics of white noise decomposed by EEMD[J]. Journal of Geomatics Science and Technology，2013，30(3)：255-259.

[14] 位秀雷，林瑞霖，刘树勇，等．基于改进小波变换的混沌信号去噪研究[J]．武汉理工大学学报(交通科学与工程版)，2013，37(5)：1062-1065．

WEI Xiulei，LIN Ruilin，LIU Shuyong，et al. Denoising for chaotic signal based on an improved wavelet transform method[J]．Journal of Wuhan University of Technology(Transportation Science & Engineering)，2013，37(5)：1062-1065．

Ade-noisingmethodforchaoticsignalsbasedonimprovedEEMD

WEI Xiulei1, LIN Ruilin2, LIU Shuyong2, YANG Qingchao2

(1. NO.91404 troops of PLA , Qinhuangdao, 066000, China; 2. College of Power Engineering，Naval University of Engineering，Wuhan 430033, China)

In de-noising process using the ensemble empirical mode decomposition (EEMD), the effective treatment of intrinsic mode functions (IMFs) is a key affecting noise reduction effect. Here, an improved EEMD de-noising method was proposed. Firstly, the useful signal details of the first IMF were extracted based on the “3σ” criterion and the singular value decomposition (SVD). Then the remaining IMFs were divided into higher frequency components and lower ones based on the consecutive mean square error (CMSE). Secondly, useful signals in higher frequency components and lower ones were extracted based on SVD and Savitzky-Golay(S-G) filtering method, respectively. Thus, the loss of useful signals in higher frequency components was avoided while parts of noise in lower frequency components were removed effectively to overcome the shortcoming of IMFs being difficult to treat with EEMD to do de-noising. In order to evaluate the effectiveness of the proposed method, the test rig of leaf spring based on the double-potential well theory was made and the test results showed that the proposed method is better than the wavelet weighted parameters method and the EEMD de-noising one.

ensemble empirical mode decomposition (EEMD); chaotic signal; singular value decomposition (SVD); denoising; Savitzky-Golay (S-G) filtering

国家自然科学基金(51579242；51179197)；国家自然科学基金青年基金(51509253)；海军工程大学科研基金(425517K143)

2016-04-18 修改稿收到日期：2016-07-16

位秀雷 男,博士生，1988年生

刘树勇 男,博士,副教授,硕士生导师，1975年生 E-mail:wxlcln@163.com

TN911

： A

10.13465/j.cnki.jvs.2017.17.006