基于分段冪函數(shù)的LMD方法及其在軸承故障診斷中的應(yīng)用

肖 露，熊邦書(shū)，莫 燕，趙平均

(1.無(wú)損檢測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(南昌航空大學(xué))，南昌 330063;2.中航工業(yè)洪都航空工業(yè)集團(tuán)有限責(zé)任公司飛機(jī)設(shè)計(jì)研究所，南昌 330024）

0 引言

滾動(dòng)軸承是旋轉(zhuǎn)機(jī)械中最重要的零部件之一，在運(yùn)行過(guò)程中，若出現(xiàn)故障，則會(huì)造成嚴(yán)重后果;因此，滾動(dòng)軸承的故障診斷對(duì)保障旋轉(zhuǎn)機(jī)械安全工作具有重要意義。

由于滾動(dòng)軸承在工作過(guò)程中受載荷、摩擦以及噪聲等因素影響，實(shí)際采集到的振動(dòng)信號(hào)是一種非平穩(wěn)的多分量調(diào)幅—調(diào)頻信號(hào)，因此時(shí)頻分析方法是目前滾動(dòng)軸承故障診斷的主流方法?，F(xiàn)有時(shí)頻分析方法可分為3類，第1類為線性時(shí)頻變換，如Potter等首次提出了通過(guò)滑動(dòng)時(shí)窗來(lái)計(jì)算頻譜的短時(shí)傅里葉變換［1］，Morlet提出了基于時(shí)間－尺度的小波分析［2］等，這類方法受到測(cè)不準(zhǔn)原理的約束，并缺少自適應(yīng)性［3］。第2類為非線性時(shí)頻變換，起源于量子力學(xué)的 Wigner分布［4］，因其良好的時(shí)頻聚集性無(wú)需加窗，得到了廣泛應(yīng)用［5］。第3類為自適應(yīng)時(shí)頻分析方法，其中最經(jīng)典的方法有Huang提出的經(jīng)驗(yàn)?zāi)Ｊ椒纸?Empirical Mode Decomposition，EMD)［6］，該方法先將信號(hào)分解為各個(gè)本征模態(tài)函數(shù)，再進(jìn)行Hilbert分解，可得到清晰的時(shí)頻分布;2005年Smith提出局域均值分解(Local Mean Decomposition，LMD)［7］，LMD將非平穩(wěn)多分量信號(hào)分解為一系列瞬時(shí)頻率有物理意義的PF(Production Function)和一個(gè)殘余分量之和，每個(gè)PF為單分量信號(hào)，結(jié)合所有PF可得到整個(gè)信號(hào)的時(shí)頻分布［8］，該方法改善了EMD的端點(diǎn)效應(yīng)，避免了EMD中的無(wú)意義負(fù)頻率，是時(shí)頻分析領(lǐng)域的重要突破。

傳統(tǒng)LMD采用滑動(dòng)平均計(jì)算局域均值函數(shù)與局域包絡(luò)函數(shù)，若滑動(dòng)跨度選取不當(dāng)會(huì)引起函數(shù)不收斂［9］，且導(dǎo)致過(guò)平滑，影響算法精確度;因此，本研究提出一種改進(jìn)算法，采用分段冪函數(shù)法分別求取信號(hào)極大值點(diǎn)集和極小值點(diǎn)集的包絡(luò)線，上下包絡(luò)線之和的均值為局域均值函數(shù)，上下包絡(luò)線之差絕對(duì)值的均值為局域包絡(luò)函數(shù)，避免過(guò)平滑現(xiàn)象，得到的局域包絡(luò)與局域均值函數(shù)更精確，并通過(guò)仿真信號(hào)分析及其在滾動(dòng)軸承故障診斷的應(yīng)用，以驗(yàn)證所設(shè)計(jì)方法的有效性。

1 LMD算法原理

LMD將多分量的非平穩(wěn)信號(hào)分解為若干個(gè)單分量的PF和一殘余分量。每個(gè)PF由一個(gè)包絡(luò)信號(hào)和一個(gè)純調(diào)頻信號(hào)相乘組成。PF的瞬時(shí)幅值由包絡(luò)信號(hào)得到，瞬時(shí)頻率由純調(diào)頻信號(hào)得到。具體方法如下:

Step 1:根據(jù)找出的信號(hào)x(t)的極值點(diǎn)集n(i)，求出第i段的局域均值點(diǎn)m(i)和包絡(luò)估計(jì)點(diǎn)a(i)，如式(1)和式(2)所示:

Step 2:分別對(duì)局域均值點(diǎn)和包絡(luò)估計(jì)點(diǎn)滑動(dòng)平均處理得到局部均值函數(shù)m11(t)與局域包絡(luò)函數(shù)a11(t)，滑動(dòng)平均計(jì)算公式為

Step 3:從信號(hào)x(t)中分離出局域均值函數(shù)m11(t)，得

Step 4:已剔除均值的信號(hào)h11(t)除以包絡(luò)信號(hào)，得

Step 5:若新信號(hào)s11(t)對(duì)應(yīng)的局域包絡(luò)函數(shù)a12(t)=1，則其為符合要求的純調(diào)頻信號(hào);若不為純調(diào)頻信號(hào)，則重復(fù)3，4步的迭代，直至其為純調(diào)頻信號(hào)s1n(t)，這里

Step 6:由迭代結(jié)果得到PF1分量

其中:s1n(t)為純調(diào)頻信號(hào)，可表示PF分量的瞬時(shí)頻率;a1(t)可表示為PF分量的瞬時(shí)幅值，表達(dá)式為

Step 7:分離出原始信號(hào)x(t)中的PF1(t)，得到u1(t)，將u1(t)作為原始信號(hào)重復(fù)前六步直至uk(t)的極點(diǎn)數(shù)不大于1，分離出所有PF分量。x(t)的最終表達(dá)式為

LMD分離出若干個(gè)PF，能夠完全表示整個(gè)信號(hào)的時(shí)頻特性，并用于多分量信號(hào)的分析和特征提取。

2 基于分段冪函數(shù)的LMD算法

LMD算法通過(guò)找到所有極值點(diǎn)集，應(yīng)用滑動(dòng)平均法求取局域均值函數(shù)與局域包絡(luò)函數(shù)?；瑒?dòng)平均法會(huì)引起過(guò)平滑，且普遍選擇最大局域均值的1/3為滑動(dòng)平均跨度，易造成沖擊類信號(hào)處理失真。針對(duì)上述問(wèn)題提出一種分別找出信號(hào)的極大極小值點(diǎn)集，用分段冪函數(shù)求取包絡(luò)線以得到局域均值函數(shù)和局域包絡(luò)函數(shù)的算法。

2.1 分段冪函數(shù)方法原理

分段冪函數(shù)插值方法原理為:設(shè)有插值點(diǎn)P1，P2，P3，首先連接 P1，P3，過(guò) P2作平行于 P1P3的直線，該直線與x=x1和x=x3分別相交于A，B點(diǎn)(圖1)。以為x軸，點(diǎn)P2為原點(diǎn)建立新的坐標(biāo)系，則點(diǎn)A在新坐標(biāo)中映射為x軸上坐標(biāo)為(x1－x2，0)的點(diǎn)A'，點(diǎn)B在新坐標(biāo)中映射為 x軸上坐標(biāo)為(x3－x2，0)的點(diǎn)B'。取與x軸垂直的向量得到點(diǎn)P1在新坐標(biāo)中的映射點(diǎn)，取得到點(diǎn)P2在新坐標(biāo)中的映射點(diǎn)(圖2)。在圖2中作冪函數(shù)曲線y=xβ(β＞1)，同分別位于x軸左右兩側(cè)的與相交于C，D 點(diǎn)。

圖1 分段冪函數(shù)算法中A，B的含義Fig.1 The meaning of A and B in piecewise power function algorithm

圖2 分段冪函數(shù)算法中C，D的含義Fig.2 The meaning of C and D in piecewise power function algorithm

則由此構(gòu)建函數(shù)［9］:

其中，

設(shè)過(guò)點(diǎn)P2與P1P3平行的直線為g(x)，插值公式為

則過(guò)點(diǎn) P1，P2，P3的曲線為

同理，P2，P3，P4的插值結(jié)果為 s2(x)，則 P2，P3之間的曲線為

因此可得，對(duì)于點(diǎn)集 Pi(i=1，2，…，n)，將Pi－1，Pi，Pi+1分段冪函數(shù)插值得到 si(x)，Pi，Pi+1，Pi+2插值得到 si+1(x)，則依照式(15)得 Pi與Pi+1之間最終的曲線表達(dá)形式s(x)。函數(shù)s(x)在定義域范圍內(nèi)一階連續(xù)可導(dǎo)，即曲線s(x)為單值光滑曲線。文獻(xiàn)［10］證明了分段冪函數(shù)誤差與β成反比，誤差范圍比3次樣條插值等要小，說(shuō)明該方法適合于非平穩(wěn)非線性振動(dòng)信號(hào)的包絡(luò)生成。

2.2 基于分段冪函數(shù)的LMD算法過(guò)程

基于分段冪函數(shù)的LMD算法與傳統(tǒng)LMD算法在求取局域均值函數(shù)和局域包絡(luò)函數(shù)上不同，其他步驟相同。計(jì)算局域均值和局域包絡(luò)函數(shù)的主要步驟如下:

Step 1:分別計(jì)算振動(dòng)信號(hào)的極大值點(diǎn)集maxi(i=1，2，3，…，M)和極小值點(diǎn)集 mini(i=0，1，2，3，…，N)。

Step 2:采用分段冪函數(shù)法，即利用式(12)和式(15)對(duì)極大值點(diǎn)集擬合得到包絡(luò)線Emax(t)，對(duì)極小值點(diǎn)集擬合得到包絡(luò)線Emin(t)。

Step 3:由包絡(luò)線Emax(t)和Emin(t)得到局域均值函數(shù)m(t)和局域包絡(luò)函數(shù)a(t)，其計(jì)算公式為

與傳統(tǒng)LMD相比，本算法不需重復(fù)計(jì)算，且其誤差界較小，故在計(jì)算效率和精確度上有較大的優(yōu)勢(shì)。

3 仿真信號(hào)分析

為了驗(yàn)證本改進(jìn)LMD算法的性能，采用一個(gè)非線性信號(hào)進(jìn)行仿真實(shí)驗(yàn)，其表達(dá)式為

仿真信號(hào)由3個(gè)分量相加而成，其中:x1為調(diào)頻調(diào)幅信號(hào)，基頻為100 Hz，調(diào)制頻率為5 Hz;x2為純調(diào)頻信號(hào)，基頻為30 Hz，調(diào)制頻率為7.5 Hz;x3為頻率為7 Hz的正弦信號(hào)。信號(hào)x及3個(gè)分量x1，x2，x3的波形圖如圖3所示。

圖3 各組成分量的仿真信號(hào)Fig.3 The simulation signals of all components

圖4、圖5分別給出了采用傳統(tǒng)LMD算法和改進(jìn)算法對(duì)仿真信號(hào)的分解結(jié)果?？梢钥闯?圖4中傳統(tǒng)LMD分解的殘余分量Res起伏較大，而圖5中殘余分量Res比較平緩，無(wú)明顯波動(dòng)極值點(diǎn)，殘余分量值較小。為了量化評(píng)估2種算法的性能，分解信號(hào)的相關(guān)系數(shù)、均方誤差(Mean Squared Error，MSE)如表1所示，從表1可以看出:1)改進(jìn)算法分解的3個(gè)PF分量的均方誤差明顯小于傳統(tǒng)LMD方法，其分解精度更高，分解更接近于理想結(jié)果;2)改進(jìn)算法所得Res的相關(guān)系數(shù)小于0.1，為偽分量，可不予考慮，且Res的相關(guān)系數(shù)的值小于傳統(tǒng)LMD，說(shuō)明改進(jìn)算法分解得更徹底。計(jì)算2種算法運(yùn)行時(shí)間，傳統(tǒng)LMD算法為3.36 s，改進(jìn)算法為3.08 s，節(jié)省了運(yùn)算時(shí)間。

圖4 傳統(tǒng)LMD算法分解結(jié)果Fig.4 The decomposition result of the traditional LMD algorithm

圖5 改進(jìn)算法的分解結(jié)果Fig.5 The decomposition result of the improved algorithm

表1 兩種算法性能對(duì)比Table 1 Performance comparison of two kinds of algorithm

4 滾動(dòng)軸承故障診斷實(shí)驗(yàn)

本研究采用美國(guó)凱斯西大學(xué)電氣工程與計(jì)算機(jī)科學(xué)系軸承實(shí)驗(yàn)數(shù)據(jù)［11］進(jìn)行實(shí)驗(yàn)，測(cè)試軸承為6502—2RS SKF深溝球，其鋼球直徑d=8 mm，球組節(jié)圓直徑D=40 mm，鋼球數(shù)N=9，接觸角α=0°。數(shù)據(jù)的采樣頻率為12 kHz，采樣點(diǎn)N=1024，發(fā)動(dòng)機(jī)轉(zhuǎn)速n=1797 rpm/min，旋轉(zhuǎn)頻率fr=n/60。

圖6、圖7分別給出了信號(hào)經(jīng)傳統(tǒng)LMD算法和改進(jìn)算法分解的結(jié)果，可見(jiàn):傳統(tǒng)LMD經(jīng)過(guò)4次分解才能得到平緩的殘余分量Res，PF4及殘余分量Res與原信號(hào)相關(guān)系數(shù)分別為0.028 0和0.0241，遠(yuǎn)小于0.1，為不含有效信息的偽分量，可剔除，但根據(jù)Res單調(diào)的終止條件仍需繼續(xù)分解，分解效率低;改進(jìn)算法只需進(jìn)行3次分解就可得到有效分量以及符合要求的殘余分量Res，分解效率高。

圖6 傳統(tǒng)LMD算法實(shí)例信號(hào)分解結(jié)果Fig.6 The instance signal decomposition result of the traditional LMD algorithm

圖7 改進(jìn)算法實(shí)例信號(hào)分解結(jié)果Fig.7 The instance signal decomposition result of the improved algorithm

將本研究設(shè)計(jì)的改進(jìn)算法應(yīng)用于軸承故障診斷。由于損傷軸承零部件會(huì)周期碰撞損傷位置，形成具有不同通過(guò)頻率的減幅振蕩，可以此作為判別不同故障的依據(jù)。采用滾珠故障數(shù)據(jù)，依據(jù)式(19)結(jié)合上述的實(shí)驗(yàn)滾動(dòng)軸承參數(shù)可得滾珠故障頻率 fd=141.17 Hz。

分別對(duì)改進(jìn)算法的各PF分量進(jìn)行分析。根據(jù)計(jì)算PF3與殘余分量的相關(guān)系數(shù)均小于0.1，為偽分量，可不予考慮。PF1和PF2頻譜圖分別如圖8a、8b所示:

圖8 故障診斷頻譜圖Fig.8 The spectrum diagram of fault diagnosis

圖8中PF2分量在140.8 Hz處存在明顯的譜線，與式(19)的理論故障頻率結(jié)果相近，可確定為故障頻率;在 PF1 分量的 58.59、117.2、210.9 Hz及PF2分量94.19 Hz處存在明顯譜線，分別與軸承基頻的2倍頻59.9 Hz、4倍頻119.8 Hz、7 倍頻209.65 Hz及3 倍頻89.85 Hz接近，與信號(hào)存在基頻倍頻這一實(shí)際相符，由此可以判定軸承出現(xiàn)了以140.8 Hz為特征頻率的滾珠故障。上述實(shí)驗(yàn)驗(yàn)證了本研究所設(shè)計(jì)的算法可應(yīng)用于滾動(dòng)軸承故障診斷。

5 結(jié)論

1)本研究提出了一種基于分段冪函數(shù)法的LMD改進(jìn)算法;

2)將改進(jìn)算法分別采用信號(hào)上下極值點(diǎn)包絡(luò)線計(jì)算LMD中局域均值函數(shù)與局域包絡(luò)函數(shù)，避免了傳統(tǒng)LMD算法的過(guò)平滑，分解結(jié)果更精確，運(yùn)算時(shí)間縮短;

3)通過(guò)仿真信號(hào)分析和軸承故障診斷實(shí)驗(yàn)，驗(yàn)證了方法的優(yōu)越性和有效性。

［1］Auslander L，Buffalano C，Orr R S，et al.Comparison of the Gabor and short-time Fourier transforms for signal detection and feature extraction in noisy environment［J］.Proceedings of the International Society for Optical Engineering，1990，1348:230－247.

［2］羅光坤.Morlet小波變換理論與應(yīng)用研究及軟件實(shí)現(xiàn)［D］.南京:南京航空航天大學(xué)，2007:13－23.

［3］孫云嶺，樸甲哲，張永祥.Wigner-Ville時(shí)頻分布在內(nèi)燃機(jī)故障診斷中的應(yīng)用［J］.振動(dòng)與沖擊，2004，15(6):505 －507.

［4］馬輝，趙鑫，趙群超，等.時(shí)頻分析在旋轉(zhuǎn)機(jī)械故障診斷中的應(yīng)用［J］.振動(dòng)與沖擊，2007，26(3):61－64.

［5］向玲，唐基貴，胡愛(ài)軍.旋轉(zhuǎn)機(jī)械非平穩(wěn)振動(dòng)信號(hào)的時(shí)頻分析比較［J］.振動(dòng)與沖擊，2010，29(2):42－46.

［6］Huang N E，Zheng S，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 of London，1998，454(1971):903 －995.

［7］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.

［8］任千達(dá).基于局域均值分解的旋轉(zhuǎn)機(jī)械故障特征提取方法及系統(tǒng)研究［D］.浙江:浙江大學(xué)，2008:60－65.

［9］鐘佑明，金濤，秦樹(shù)人.希爾伯特－黃變換中的一種新包絡(luò)線算法［J］.數(shù)據(jù)采集與處理，2005，20(1):13－17.

［10］劉霖雯，劉超，江成順.EMD新算法及其應(yīng)用［J］.系統(tǒng)仿真學(xué)報(bào)，2007，19(2):446 －448.

［11］Case Western Reserve University Bearing Data Center Website.Bearing Date Center Seeded Fault Test Date［DB/OL］.［2011 －05 －03］.http://www.eecs.case.edu/laboratory/bearing/.