局部均值分解方法中乘積函數(shù)判據(jù)問題研究

張 亢，程軍圣，楊 宇

(湖南大學(xué) 汽車車身先進設(shè)計制造國家重點實驗室，長沙 410082)

局部均值分解(Local mean decomposition，LMD)是一種新的自適應(yīng)時頻分析方法［1］。該方法能自適應(yīng)地將一個復(fù)雜的多分量信號分解為若干個瞬時頻率具有物理意義的乘積函數(shù)(Product function，PF)分量之和，其中每一個PF分量是一個包絡(luò)信號與一個純調(diào)頻信號的乘積，包絡(luò)信號就是該PF分量的瞬時幅值，PF分量的瞬時頻率則可由純調(diào)頻信號直接求出，進一步將所有PF分量的瞬時頻率和瞬時幅值組合，便可以得到原始信號完整的時頻分布。LMD方法避免了EMD方法中的過包絡(luò)、欠包絡(luò)、Hilbert變換帶來的邊緣效應(yīng)以及無法解釋的負頻率等問題［2］，已在腦電信號分析和機械故障診斷領(lǐng)域得到了應(yīng)用［3］。

在LMD方法中，PF分量定義為一個純調(diào)頻信號和一個包絡(luò)信號的乘積，而純調(diào)頻信號需要滿足相應(yīng)的要求，即本身的絕對值不能大于1，而且其包絡(luò)信號恒等于1，這在實際應(yīng)用中很難滿足。在實際的算法中，由于PF分量的純調(diào)頻信號是需要經(jīng)過多次的迭代才能獲得的，因此需要找到一個合適的PF分量判據(jù)，它既是可以實現(xiàn)的，同時也能盡量滿足純調(diào)頻信號的理論條件，從而使得PF分量具有物理意義。由于PF分量與經(jīng)驗?zāi)B(tài)分解［4］(Empirical mode decomposition，EMD)中的內(nèi)稟模態(tài)函數(shù)(Intrinsic mode function，IMF)分量具有相似性，而IMF分量之間已從數(shù)值上驗證了具有正交的特點［5］，因此本文擬從PF分量的正交性入手，將 IMF分量的正交性判據(jù)［6］(Orthogonality criterion，OC)引入到LMD方法中，作為PF分量迭代終止的判據(jù)。通過對仿真信號以及實際信號的分析，結(jié)果表明由正交性判據(jù)所確定的PF分量是完備的且具有物理意義，能夠?qū)崿F(xiàn)信號的正確分解。

1 LMD方法

LMD方法的本質(zhì)是通過迭代從原始信號中分離出純調(diào)頻信號和包絡(luò)信號，然后將純調(diào)頻信號和包絡(luò)信號相乘便可以得到一個瞬時頻率具有物理意義的PF分量，循環(huán)處理至所有的PF分量分離出來，便可以得到原始信號的時頻分布。對于任意信號x(t)，其分解過程如下［1］。

(1)確定原始信號所有的局部極值點ni，計算相鄰兩個極值點ni和ni+1的平均值mi:

將所有平均值點mi在其對應(yīng)的極值點ni和ni+1間延伸成一條線段，然后用滑動平均法進行平滑處理，得到局部均值函數(shù)m11(t)。

(2)采用局部極值點ni計算局部幅值ai:

將所有局部幅值點ai在其對應(yīng)的極值點ni和ni+1間延伸成一條線段，然后采用滑動平均法進行平滑處理，得到包絡(luò)估計函數(shù)a11(t)。

(3)將局部均值函數(shù)m11(t)從原始信號x(t)中分離出來，即去掉了一個低頻成分，得到:

(4)用h11(t)除以包絡(luò)估計函數(shù)a11(t)以對h11(t)進行解調(diào)，得到:

對s11(t)重復(fù)上述步驟便能得到s11(t)的包絡(luò)估計函數(shù)a12(t)，假如a12(t)不等于1，說明s11(t)不是一個純調(diào)頻信號，需要重復(fù)上述迭代過程n次，直至s1n(t)為一個純調(diào)頻信號，也即s1n(t)的包絡(luò)估計函數(shù)a1(n+1)(t)=1，所以，有:

理論迭代終止的條件為:

(5)把迭代過程中產(chǎn)生的所有包絡(luò)估計函數(shù)相乘便可以得到包絡(luò)信號(瞬時幅值函數(shù)):

(6)將包絡(luò)信號a1(t)和純調(diào)頻信號s1n(t)相乘便可以得到原始信號的第一個PF分量:

它包含了原始信號中最高的頻率成分，是一個單分量的調(diào)幅-調(diào)頻信號，其瞬時幅值就是包絡(luò)信號a1(t)，其瞬時頻率f1(t)則可由純調(diào)頻信號s1n(t)求出，即:

(7)將第一個PF分量PF1(t)從原始信號x(t)中分離出來，得到一個新的信號u1(t)，將u1(t)作為原始數(shù)據(jù)重復(fù)以上步驟，循環(huán)k次，直到uk為一個單調(diào)函數(shù)為止。

原始信號x(t)能夠被所有的PF分量和uk重構(gòu)，即:

將所有PF分量的瞬時幅值和瞬時頻率組合便可以得到原始信號x(t)完整的時頻分布。

2 PF分量迭代終止條件

由上述LMD分解算法知，每個PF分量都是通過多次迭代過程得到的，式(7)表示的是理論的迭代終止條件，由它所確定的PF分量的瞬時頻率和瞬時幅值具有最準確的物理意義，并且各個PF分量之間具有正交性，最終分解結(jié)果也是完備的。但在LMD方法的實際運用中該迭代終止條件是不可能實現(xiàn)的［1］，因此需要另外確定一種合理的、可以實現(xiàn)的判據(jù)替代理論的迭代終止條件。EMD分解中“篩分”IMF分量的過程與PF分量的迭代過程類似，而在IMF分量的判據(jù)問題上，已有一些學(xué)者做了研究，如Huang等［4］提出了標準差判據(jù)(Standard deviation，SD)，即通過計算連續(xù)兩次迭代結(jié)果的差值，并與預(yù)先確定的閾值進行比較，以此決定迭代的終止點;Rilling等［7］對SD判據(jù)進行了改進，提出了3-閾值判據(jù)(Three-threshold criterion);Xie等［8］從信號的頻率和相位信息入手，提出了帶寬判據(jù)(Bandwidth criterion，BC);Lin 等［6］根據(jù) IMF 分量正交的特點，提出了正交性判據(jù)(Orthogonality criterion，OC)等，這些判據(jù)在EMD方法的實際運用中都取得了一定的效果。借鑒上述確定IMF分量的判據(jù)的思想，并根據(jù)PF分量的正交性質(zhì)，將正交性判據(jù)OC引入了LMD方法中。

理論上每一個PF分量是一個單分量的調(diào)幅-調(diào)頻信號，而LMD方法的實質(zhì)就是將一個包含多種時間尺度特征的信號分解成有限個單分量的調(diào)幅-調(diào)頻信號之和，其中每個單分量的調(diào)幅-調(diào)頻信號代表了原信號中的一種時間尺度特征，因此不同PF分量之間是不相關(guān)的，也即正交的，即有下式成立:

但式(13)并不是嚴格準確的，這是由于LMD方法中局部均值函數(shù)和包絡(luò)估計函數(shù)是通過平滑得到的，并不是信號真正的均值和包絡(luò)值，再加上端點效應(yīng)［9］的影響，不可避免地會發(fā)生能量泄漏現(xiàn)象，因此不可能具備嚴格的正交性質(zhì)，也即式(13)的右邊不可能等于0。但可以得出如果該值越接近于0，則說明PF分量的正交性越好，分解的效果也越好的結(jié)論。

根據(jù)以上PF分量的正交性質(zhì)，將正交性判據(jù)OC引入LMD方法中作為PF分量的迭代終止判據(jù)，其定義如下:

3 應(yīng)用

考察式(15)所示的仿真信號x(t)，它由兩個調(diào)幅-調(diào)頻分量組成，其時域波形如圖1所示。

圖1 仿真信號x(t)Fig.1 Simulated signal

圖2 仿真信號x(t)的LMD分解結(jié)果Fig.2 The PFs of the simulated signal from LMD

引入正交索引參數(shù)(Index of orthogonality，IO)［10］評價LMD分解的整體的正交性，其表達式如下:

式中i，j=k+1時，表示分解得到的余量。理論上PF分量間如果完全正交，應(yīng)該有IO=0，但由第2節(jié)分析知，實際中不可能完全正交。經(jīng)計算上述仿真信號x(t)的分解結(jié)果的IO=0.002 1，很接近于0，幾乎滿足完全正交。另外x(t)分解前后的能量變化情況見表1，可以看出，總能量以及各個分量的能量分解前后幾乎一致，體現(xiàn)了分解的完備性與正確性，說明了正交性判據(jù)是合理的。

表1 仿真信號x(t)分解前后能量的變化情況Tab.1 The energy change situation of simulated signal before and after decomposition

對于實際信號，選取了某內(nèi)圈發(fā)生了局部損傷的滾動軸承故障振動加速度信號x(t)進行分析。軸承型號為6311，轉(zhuǎn)頻為20 Hz，采樣頻率為4 096 Hz，計算的理論故障特征頻率為f0=99.2 Hz。x(t)的時域波形如圖3所示。

圖3 滾動軸承內(nèi)圈故障振動信號x(t)Fig.3 The vibration signal of roller bearing with inner race defect

采用LMD方法對x(t)進行分解，并以正交性判據(jù)準則作為PF分量的迭代終止條件(采用自適應(yīng)波形匹配延拓法［9］解決端點效應(yīng))。每迭代一次，計算OC的值，最終完全分解后得到7個PF分量和1個余量，如圖4所示。每個PF分量的迭代次數(shù)以及OC值見表2。

圖4 滾動軸承內(nèi)圈故障振動信號x(t)的LMD分解結(jié)果Fig.4 The PFs of the vibration signal of roller bearing with inner race defect from LMD

表2 各個PF分量的迭代次數(shù)以及OC值Tab.2 Iterative times and OC value of each PF component

同樣為了評價分解的效果，按式(16)計算本次分解的正交索引參數(shù)IO的值。經(jīng)計算有IO=0.047 5，幾乎滿足完全正交。x(t)分解前后的能量對比情況見表3，其中Ex為信號x(t)的能量，Ei(i=1，2，…，7)為各個PF分量的能量，Eu7為余量的能量，E'x為7個PF分量與余量的能量之和，可以看出分解前后的能量差很小，體現(xiàn)了分解的完備性。進一步對前兩個PF分量PF1(t)和PF2(t)的瞬時幅值a1(t)和a2(t)作頻譜分析，結(jié)果如圖5和圖6所示。圖5中，在內(nèi)圈故障特征頻率99.2 Hz處存在著明顯的譜線，同時存在著轉(zhuǎn)頻的2倍頻、3倍頻成分;圖6中也可以找到99.2 Hz的譜線，但主要成分是轉(zhuǎn)頻的2倍頻。因此通過幅值譜可以判定軸承存在內(nèi)圈局部故障，與實際情況相符。以上分析說明本次LMD分解結(jié)果是合理的，正交性判據(jù)是合適的。

表3 軸承故障振動信號x(t)分解前后能量變化的情況Tab.3 The energy change situation of the vibration signal with inner race defect before and after decomposition

標準差判據(jù)(Standard deviation，SD)是一種在EMD分解中得到廣泛應(yīng)用的判據(jù)準則，其定義如下:

其中hik(t)為求第i個IMF分量時的第k次“篩分”結(jié)果。對比EMD和LMD方法的理論可知，如果將式(17)中的hik(t)用hij(t)=Si(j-1)(t)-mij(t)替換，其中Si(j-1)(t)為分解第i個PF分量時的第j-1次迭代結(jié)果，mij(t)為Si(j-1)(t)的局部均值函數(shù)，便可以得到LMD方法中的SD判據(jù)。為了與OC判據(jù)進行比較，將SD判據(jù)同樣應(yīng)用于圖3所示的滾動軸承故障振動信號，最終分解得到了7個PF分量和1個余量，如圖7所示。對比圖4分解效果沒有很大差別。圖8為第一個PF分量PF1(t)的幅值譜，同樣能找到內(nèi)圈故障特征頻率。但是進一步計算了本次分解的正交索引值為IO=0.085 3，以及能量值見表4，可以看出IO值和分解前后的能量差值都要比采用OC判據(jù)分解時的大，說明從總體上分析采用SD判據(jù)分解時的正交性和完備性都不如OC判據(jù)好。另外各個PF分量的迭代次數(shù)見表5，其值明顯要大于采用OC判據(jù)分解時的情況，說明計算效率OC判據(jù)較SD判據(jù)高。

表4 軸承故障振動信號x(t)分解前后能量變化的情況Tab.4 The energy change situation of the vibration signal with inner race defect before and after decomposition

表5 各個PF分量的迭代次數(shù)Tab.5 Iterative times of each PF component

4 結(jié)論

利用PF分量的正交性質(zhì)，將正交性判據(jù)OC引入到LMD方法中，解決了LMD方法中PF分量的判據(jù)問題。通過對仿真信號和實際信號的分析，結(jié)果表明該判據(jù)不但充分保證了LMD分解的正交性質(zhì)，而且穩(wěn)定性好、收斂速度快，得到的分解結(jié)果具有真實的物理意義并且是完備的，這對LMD方法的工程實際應(yīng)用是非常有意義的。但值得提出的是，由于平滑步長會對局部均值函數(shù)和包絡(luò)估計函數(shù)的形狀造成影響，因此勢必會影響到該判據(jù)的值，而影響的規(guī)律如何，這需要做進一步研究。

［1］Smith J S.The local mean decomposition and its application to EEG perception data［J］.Journal of the Royal Society Interface，2005，2(5):443-454.

［2］Boashash B.Estimating and interpreting the instantaneous frequency of a signal-part 1:Fundamentals［J］.Proc.IEEE，1992，80(4):520-538.

［3］程軍圣，楊 宇，于德介.局部均值分解方法及其在齒輪故障診斷中的應(yīng)用［J］.振動工程學(xué)報，2009，22(1):76-84.

［4］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］.Proc.R.Soc.Lond.A，1998，454:903-995.

［5］ Wu Z，Huang N E.Ensemble empirical mode decomposition:a noise-assisted data analysis method［J］.Advances in Adaptive Data Analysis，2009，1(1):1-41.

［6］Li L，Ji H B.Signal feature extraction based on an improved EMD method［J］.Measurement，2009，42:796-803.

［7］Rilling G，F(xiàn)landrin P，Goncalves P.On empirical mode decomposition and its algorithms［C］.Proceedings of IEEE EURASIP Workshop on NonlinearSignaland Image Processing.Grado，Italy，2003.

［8］Xie Q W，Xuan B，Li J P，et al.EMD algorithm based on bandwidth and the application on one economic data analysis［C］. European Signal Processing Conference，Poznan，Poland，2007:2419-2423.

［9］張 亢，程軍圣，楊 宇.基于自適應(yīng)波形匹配延拓的局部均值分解端點效應(yīng)處理方法［J］.中國機械工程，2010，21(4):457-462.

［10］邱綿浩，劉 菁，叢 華.基于B樣條插值曲線直接篩選的EMD及其在機械振動信號處理中的應(yīng)用［J］.裝甲兵工程學(xué)院學(xué)報，2007，21(3):29-33.