经验模态分解算法中端点问题的处理(1)

x=[0 30 60 90 120 150 180 210 240 270 300 330 360];y=[-0.0167 -1.0927 -1.8725 -2.3586 -2.3061 -1.9576 -0.9574 -0.0080 0.8896 1.3877 1.1139 0.8517 -0.0167];fun=@(a,t) a(1)+a(2)*sind(t+a(3)) %matlab7.0以上版本，否则用inline%fun=inline('a(1)+a(2)*sind(t+a(3))','a','t')a0=[-0.5 -1.9 -0.079];a=nlinfit(x,y,fun,a0)t=0:5:360;yf=fun(a,t);plot(x,y,'o',t,yf)结果：fun =@(a,t) a(1)+a(2)*sind(t+a(3))a =-0.5239 -1.8995 -14.2382经验模态分解算法中端点问题的处理摘要：经验模态分解(EMD)方法就是对非线性、非平稳信号运用时间区域序列的上下包络线的均值得到瞬时平衡位置，将被分析信号分解成一组相互独立的稳态和线性的固有模态函数(IMF)数集。

经验模态分解(EMD)方法是基于原始信号本事出发，经过筛选先把频率高的IMF 分量分离出来，然后在分离频率较低的IMF分量。

其实质就是利用时间特征尺度来获取原始信号数据中的振荡模态，本文对经验模态分解算法中端点问题的处理进行研究。

关键词：经验模态分解算法端点函数经验模态分解(EMD)方法被提出后在各个领域普遍的应用，其具有直观、简单、自适应、完备性和正交性以及调制特性等一系列良好的特点。

(1)自适应性经验模态分解(EMD)方法的自适应性表现为自适应生成基函数。

在整个筛选分解过程中是根据原始信号自己的时间特征尺度实现的，不需要事先设立任何基函数。

这与傅立叶变换和小波变换有着根本性的不同。

傅立叶变换和小波变换需要事先设定谐波基函数和小波基函数，他们是先验性的。

可以说在理论上，经验模态分解(EMD)方法适用任何信号的分解，其在对非线性、非平稳信号的处理上的优越性是其他时频分析方法无法比拟的。

经验模态分解(EMD)方法的自适应性还表现为自适应滤波特性。

经验模态分解方法是基于原始信号本身出发，经过筛选先把频率高的IMF分量分离出来，然后在分离频率较低的工MF分量。

这些不同频率成分以及带宽都是随原始信号的变化而改变的。

因此，EMD方法可被视为是一组具有自适应性能的带通滤波器，它的截止频率和带宽均随着原始信号的变化而自动改变的，随着信号分析的目的改变而自动变化的。

这些尺度范围和频率成分均随着原始信号的变化而自动改变的，这样原始信号的特征可在不同分辨率下被表示，实现自适应多分辨率。

(2)正交性与完备性所谓信号分解方法的完备性，是指可以从被分解后的信号的各个分量还原出原始信号的性质。

经验模态分解(EMD)方法的本身就已经证明了其是完备的。

可以得到证明EMD方法的完备性。

且从EMD整个分解过程和结果都说明EMD方法的完备性。

所谓信号的正交性指的是被分解后的信号的各个分量之间相互正交的性质。

例如频率不同的两个正弦信号它们是相互正交的。

在EMD方法中，根据IMF的概念，每一个IMF分量应该是在局部应该是相互正交的。

此处的正交性是在局部意义上而言的。

对于有一些特殊数据，有可能会出现两个相邻的分量在不同的时间段内含有相同的频率成分，因此在全局意义上不正交。

由于一般在实际进行经验模态分解时采用的都是截取数据的有限长度，这样即使对于不同频率的纯正弦波形叠加信号进行分解也会有严重泄漏。

泄漏的程度一般与数据的长度以及分解的结果是密切相关的。

黄通过大量的实验数据证明EMD的泄漏一般小于1%;对于极短的数据为5%,与正弦型傅立叶分解在同一数量级上。

据此可以认为EMD分解得到的各个IMF分量近似正交的。

(3) IMF分量的调制特性由固有模态函数的概念可知，对于任意信号被分解为有限工MF分量，这些分量可以是幅值和频率调制的。

任何频率随时间的变化都可定义为频率调制。

频率调制有两种概念:一是波间调制;二是波内调制。

EMD分解得到的各个IMF分量不仅含有原始信号的非线性以及非平稳特性，而且工MF分量有波内调制特性，能用一个IMT表示由不同傅立叶频率描述的同一分量。

1 经验模态分解和主成分分析1.1 经验模态分解经验模态分解算法的主要目的是将待分析信号分解为一系列表征时间尺度的IMF 分量,要求IMF 分量必须满足两个条件：IMF 的极值点个数与过零点个数不超过1；由极大值点和极小值点确定的包络线均值为零.对信号()x t 进行EMD 分解的步骤如下[6]：Step1 确定()x t 的所有极大值和极小值，分别对极大值点和极小值点进行三次样条插值,构造()x t 的上下包络线()x up t 和()x low t ，计算上下包络线的均值1()(()())/2up low t t t =+m x x ；Step2 计算()x t 和()m t 之间的差值, 11()()()t t t =-g x m ； Step3 判断1()t g 是否是一个IMF,3.1) 如果1()t g 符合IMF 的定义条件，是一个IMF ，则抽取1()t g 作为第一个IMF 分量,令11()()imf t g t =，并求原信号与1imf 之间的差值1()r t ， 11()()()rt x t imf t =-。

3.2) 如果1()t g 不是一个IMF,则将1()t g 视为一个新的信号序列，重复步骤1和步骤2，求其包络均值11()t m 及1()g t 与11()m t 间的差值11()g t 。

对11()g t 重复上述过程n 次，直到1()n g t 符合IMF 的定义条件，则令11(1)1()()n n n g t g t m -=-为()x t 的第一个IMF 分量，并求原信号与1imf 之间的差值1()r t ， 11()()()rt x t imf t =-； Step4 将1()r t 作为一个新的“原始”信号，重复步骤（step1-step3），抽取第2个内蕴模态函数分量2imf ，令212()()()r t r t im f t =-；将2()r t 作为一个新的“原始”信号，抽取第3个内蕴模态函数分量3i mf ， ；以此类推，直到第K 次的余项(1)K K k r r imf -=-满足终止条件，则停止迭代，()x t 的EMD 分解完成。

EMD 分解结束后,原始信号()x t 可被表示为各IMFs 和一个余项之和：1()()()x imf r Kk K k t t t ==+∑其中imf k 表示第k 个IMF 分量.如果()x t 被零均值高斯白端点污染,则imf k 中所含端点仍近似服从零均值正态分布[14],即可设imf ky n k k=+(1)其中y k 表示没被污染的原始信号,n k 表示所含端点,且n k ～2(0,)k N σ.对于含随机端点的信号,先分解出的IMF 分量通常对应于信号的高频端点,若去除几个先分解出的IMF,把剩余IMF 进行重构,则可减弱信号的端点.1.2 主成分分析主分量分解(PCA)是统计分析中常用的多通道数据分析方法[15],被广泛应用于数据的降维和端点处理中.设M 道原始数据构成一个N M ⨯的数据阵()12,,,T T T T M =X x x x,令()x x x i i i E =-,其中()x i E 为x i 的期望,记()12,,,TT T TM =X x x x ,则X 的协方差矩阵为())T M M M N N M E ⨯⨯⨯=C (X)(X ,通过奇异值分解(SVD),协方差矩阵C M M ⨯可被写为)T M M M M M M M M ⨯⨯⨯⨯=C U Λ(U这里()12,,,ΛM diag λλλ= 为对角矩阵,M λλλ≥≥≥ 21为C 的特征值,()12,,,M =U u u u 为特征值所对应的特征向量组成的正交矩阵.令(T M N M M M N⨯⨯⨯=P U )X(2)称()12,,,TT T TM=P p p p 的各行为X 的主分量,它们在P 中依贡献率大小排序.p i 对应的特征值与特征值总和的比∑=Mi i i1λλ称为该主分量的贡献率,表征该主分量代表原始信号能量的百分比.如果仅取前H 个主分量()12,,,,,,T T T T T T H =P p p p 00 重构原始数据,则重构的数据为X U P M N M M M N⨯⨯⨯= 端点信号在经PCA 处理时,由于1,,H M λλ+ 对应的主分量1,,p p H M + 包含了信号中的大部分端点,在重构时直接丢掉；而代表了信号主要特征的主分量1,,p p H 被保留,因而PCA 可以有效的端点处理.2 利用3σ准则提取1imf 中的细节信息在最初的EMD 端点算法中,通常认为第一层端点IMF 全部由端点构成,并由此出发推导其他层IMF 中所含端点的能量.但随着研究的深入,逐渐发现1imf 中仍含有一定量的信号细节信息[12,13].对1imf 进行适当的处理,提取其所含的信号细节信息并加以保留,会提高端点效果,利用处理后的1imf 估计其余IMF 中所含端点的能量也更准确.但由于先验知识很少,所以对1imf 进行处理是个难题.由EMD 的研究可知,在1imf 中端点占绝大部分,而仅含有少量的信号细节信息,而且所含端点仍近似服从零均值正态分布,所以非常适合采用“3σ法则”进行细节信息提取[16].由公式(1)可知,1imf 满足加性端点模型1imf 11y n =+ 且1n ～21(0,)N σ.根据“3σ”法则,端点1n 的分布满足{}11|[]|399.73%n P i σ≤=即端点1n 落在11[3,3]σσ-之间的概率为0.9973,而落在31σ之外的概率仅约为0.003.因此如果1[]imf i 的值没有落在11[3,3]σσ-之内,则可认为1[]imf i 中必然含显著误差,也即有必然含有信号信息1y ,需要予以保留.利用“3σ法则”对1imf 进行细节信息提取可表示为：{11111[],([])3[]0([])3d 1imf imf imf ,imf i if abs i i if abs i σσ≥=<其中d1imf 表示从1imf 提取出的信号细节；端点方差21σ采用文献[17]中提出的方法进行估计,即21()0.6745HH Median σ=,这里HH 表示1imf 的高频子带小波系数.3 各层IMF 中所含端点能量的估计利用“3σ法则”对1imf 进行细节信息提取后,可求出1imf 中所含端点的能量]1[W [18]：21111111[1]()()([][])imf imf imf imf imf imf Mdd T d i W i i ==--=-∑(3)假设imf k (2≥k )中所含端点的能量为][k W ,则][k W 是未知的,由于信号和端点混杂在一起,因此一般情况下并不能求出][k W .但通过含噪信号经EMD 分解后的端点能量模型,可对][k W 进行近似计算.被白端点污染的信号经EMD 分解后,如果第一层内蕴模态函数1imf 中所含端点的能量为]1[W ,则imf k 中所含的端点的能量][k W 可由下式求出[18],][k W kW ρβ]1[=,2≥k(4)其中719.0≈β,01.2≈ρ.因此，求出1imf 中所含端点的能量]1[W 后, 即可通过公式(4)估计imf k 中所含端点的能量][k W .4 根据端点能量利用PCA 去除imf k (2≥k )中的端点PCA 是一种自适应的分解方法,信号经PCA 分解后各主分量间互不相关,而且按照贡献率选择合适主分量重构后,能有效端点处理并保留信号绝大部分的主特征信息. imf k(2≥k )经PCA 分解后,信号和端点能够被有效分离,如果想较好去除imf k 中端点,必须要选择合适个数的主分量进行重构,通常根据前H 个主分量的累计贡献率⎪⎭⎫ ⎝⎛=∑∑==Ni i Hi ir 11λλ来确定保留的主成分分量的个数.但累积贡献率r 的选择并不是一个简单的事情：1)累计贡献率r 取得太大,会残留较多的端点,导致端点不能完整的去除；累计贡献率r 取得太小,又会损失较多的信号细节信息；2) 每一层端点项IMF 中所含端点的强度并不相同,因此在对不同层的IMF 进行PCA 端点时,累计贡献率r 不能取固定的值.为了更好地端点处理,本文根据k imf 中所含端点能量的比例,提出了一种自适应确定累计贡献率r 的方法.由公式(1)可知 imf y n k k k =+,()0n k E =,(2≥k ).为了表示方便,设X =imf T k k ,()()X =imf imf X X T T k k k k k E E -=-,()X C X X k Tk k E =为X 的协方差矩阵,其中M 表示k imf 的长度,显然X k 所含端点与imf Tk 所含端点相同.假设X k 经PCA 分解后的主分量为()12,,,TT T TM=P p p p ,如果选择前H个主分量()12,,,,,,T T T T T T H =P p p p 00 进行重构以去除X k中的端点,则可得到端点后的信号 X UP k= ()()1211,,,,,,,,HTTT TT M Hi ii ===∑u u u p p 00u p(5)此时从X k 中删除的端点为1ΔX X X u pMk k kiii H =+=-=∑(6 )设X k 和ΔX k 的能量分别为()X ε、()ΔX ε.在PCA 端点中,通常认为前几项主分量包含了信号的主要特征信息,而比较靠后的主分量主要由端点构成,而且按照主分量对总能量的贡献率选择应保留的主分量个数.因此在利用PCA 对X k 端点时,如果选择合适的H ,使删除的端点ΔX k 的能量与X k 本身所含的端点能量相同,也即使得()ΔX k ε[]W k =则可认为X k 中的端点基本被全部去除,达到了较好的端点效果.上式等价于()[]()()ΔX X X k k k W k εεε=(7)由公式(2)和U 的正交性可知：X k UP =1u pMi ii ==∑,p u X Ti i k =,所以信号X k 的能量为：()X X X TTk kε==11()()u p u p NNTi i j j i j ==∑∑11p u u p N N T T iij j i j ===∑∑21p Nii ==∑1p p NTi i i ==∑1u X X u NT T ik kii ==∑1(1)X u C u NT i ii N ==-∑1(1)Nii N λ==-∑(8)而所删除的端点ΔX k 的能量可表示为()()ΔX ΔX ΔX Tk k k ε==11p u u pN NT T iijji H j H =+=+∑∑1p pNT iii H =+=∑1uX X u NT T i k k ii H =+=∑1(1)Xu Cu NTiii H N =+=-∑1(1)Nii H N λ=+=-∑(9)由公式（8）和（9）可知,被删除端点ΔX k 的能量与X k 的能量之比为11()()ΔX X NNk i i i H i k ελλε=+=⎛⎫= ⎪⎝⎭∑∑.所以,为了使X k 中的端点被完整去除,应选择合适的H 使11NNii i H i λλ=+=⎛⎫ ⎪⎝⎭∑∑[]()X kW k ε=成立 .但在选择H 时,很难保证使得11NNii i H i λλ=+=⎛⎫ ⎪⎝⎭∑∑[]()X kW k ε=恰好成立,本文中对H 按照以下方法进行取值：如果存在β使得式(10)成立, 则令β=H ,11NNi i i i βλλ=+=⎛⎫ ⎪⎝⎭∑∑[]()X k W k ε≤≤1N Ni i i i βλλ==⎛⎫ ⎪⎝⎭∑∑(10)应保留的主分量个数H 确定后,根据公式(5)可求出X k 端点后的信号X k,因为()X =i m f i m f TTk kk E -,所以k imf 端点后的值为：()()imf X imf Td T k d kE =+ .本文所提出的基于PCA 的EMD 端点算法具体步骤为：Step1 对信号()x t 进行EMD 分解,设分解后的IMF 为1imf ,…, k imf ,余项为r K ； Step2 对1imf 采用“3σ法则”提取信号细节信息,设提取的细节信息为1imf d； Step3 根据公式(3)求1imf 所含端点的能量,并利用公式(4)估计k imf (2≥k )中所含端点的能量；Step4 对k imf (2≥k )进行PCA 分解,根据公式(10)选择合适个数的主分量进行重构端点,设k imf 端点后的值为imf dk ；Step5 累加全部imf dk (K k ≤≤1)和余项r K ,得到端点后信号()xt . 5 实验结果与分析当一维信号长度M 较大时,其协方差矩阵M M C ⨯的规模较大,直接进行PCA 变换运算量很大.为了降低运算复杂度,本文按照文献[19]中所提出的嵌入方法首先将一维信号转换为多维信号后,再进行PCA 分解,此时分解结果的特征和性质均保持不变.为了分析所提出算法的端点性能,分别对模拟信号和真实信号进行端点实验.模拟信号利用Matlab 中的wnoise 函数生成,分别生成“Blocks ”,“Bumps ”,“Heavy sine ”和“Doppler ”等四类具有典型特征的测试样本；真实信号选择来自Bell 实验室的一段电力系统信号. 为了对比端点效果,对端点信号分别采用基于Bayesian 阈值的小波端点法(Bayesian-Wavelet)、基于模态单元的EMD 阈值法(Mode-EMD)和本文提出的基于PCA 的EMD 端点法 (PCA-EMD) 进行端点. 在利用Bayesian-Wavelet 法进行端点时，小波基选用“db8”小波,分解层数取为10；在Bayesian-Wavelet 和Mode-EMD 的端点中,均采用硬阈值法. 本文采用均方误差MSE 和信噪比SNR 来评估算法的性能：信噪比越大,均方误差越小,表明端点效果越好.5.1 对模拟信号的端点在对模拟信号的实验中,首先利用wnoise生成信噪比分别为SNR=0,5dB,10,15,20dB的测试信号,信号长度L取为4096. 图1是SNR=5dB时,测试样本“Doppler”及不同方法端点后的实验结果,在图1(c)-(f)中,虚线为原始信号,实线为端点后信号.对3种方法端点后的结果分别计算MSE和SNR(见表1).可以看出,在当SNR=5dB时,“Doppler”信号经PCA-EMD 方法端点后效果相对较好,与Bayesian-Wavlet算法相比,SNR提高了2.509,MSE约减小了0.053；与Mode-EMD算法先比,SNR约提高了1.529,MSE约降低了0.017.不同信噪比的模拟信号经三种方法端点后的MSE和SNR如表1所示,通过比较可知,本文提出的PCA-EMD方法总体端点效果要优于Bayesian-Wavelet方法和Mode-EMD方法,端点后的信号更接近原始信号；但当信噪比增大时,PCA-EMD与Mode-EMD之间的差距在逐渐减小,当信噪比增加到20dB时,本文方法与Mode-EMD算法的端点结果已非常接近.(a)(b)(c)(d)(e)Fig.1 Doppler signal de-noising results comparison (a)Original Doppler signal(b)noisy Doppler signal (d) Bayesian Wavelet de-noising result (d)Mode-EMDde-noising result (e) PCA-EMD de-noising result 图1 Doppler信号的端点结果比较 (a)原始Doppler信号 (b)含噪Doppler信号 (c) Bayesian Wavelet端点 (d) Mode-EMD端点 (e)PCA-EMD端点Table 1 Results of experiments using noisy simulated signals5.2 对实际电力信号的端点图2(a)是在有端点干扰的环境下采集到的一段电力系统信号,对该信号分别利用三种方法进行端点,端点后的结果如图2(b)-(d)所示.从图2可以看出,PCA-EMD端点后信号的细节和突变部分保持较好,与Bayesian-Wavelet和Mode-EMD端点结果相比,端点明显有所减少.分别计算各端点信号的MSE和SNR,结果如下：Bayesian-wavelet端点后,MSE=81.8073,SNR=32.023；Mode-EMD端点后,MSE=73.0587, SNR=32.5142,；PCA-EMD端点后,MSE=70.9869,SNR=32.6391.可以看出,PCA-EMD方法端点后的MSE最小,而SNR最大.这表明电力信号经PCA-EMD方法端点后,信号细节的保留程度和信号的复原程度都更好一些.可见,对实际信号进行端点时,PCA-EMD方法与Bayesian-Wavelet和Mode-EMD方法相比,端点能力也有一定程度的提高,是一种比较有效的端点方法.实验在matlab7.8.0环境下进行, 计算机内存为2G，cpu主频为3.06GHz， EMD采用flandrin 提供的pack_EMD (http://perso.ens-lyon.fr/patrick.flandrin/emd.html)。