基于VMD-SVD能量標(biāo)準(zhǔn)譜-Teager能量算子的軸承故障診斷方法

周 洋，向 陽，黃陳哲，黃進(jìn)安

（1.武漢理工大學(xué) 能源與動(dòng)力工程學(xué)院，武漢 430063；2.武漢理工大學(xué) 船舶動(dòng)力系統(tǒng)運(yùn)用技術(shù)交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，武漢 430063）

滾動(dòng)軸承作為減小摩擦和磨損的元件被廣泛應(yīng)用于各種旋轉(zhuǎn)機(jī)械設(shè)備。當(dāng)滾動(dòng)軸承表面出現(xiàn)損傷時(shí)，會(huì)和配合元件之間發(fā)生碰撞，從而激起旋轉(zhuǎn)系統(tǒng)的各階固有振動(dòng)，使得產(chǎn)生的振動(dòng)信號(hào)頻譜中包含多個(gè)共振峰值，使得機(jī)械設(shè)備產(chǎn)生異常的振動(dòng)與噪聲[1]。然而輕微故障引起的沖擊常被強(qiáng)背景噪聲淹沒，導(dǎo)致難以提取故障特征頻率，因此對(duì)軸承進(jìn)行故障診斷，尤其是早期微弱故障診斷一直是學(xué)者研究的熱點(diǎn)。軸承故障產(chǎn)生的振動(dòng)信號(hào)具有典型的周期性、非平穩(wěn)特性，包含了豐富的工作狀態(tài)信息[2]，同時(shí)故障信號(hào)中常摻雜大量背景噪聲，常用的頻譜分析方法處理此類信號(hào)效果差，難以提取故障特征頻率。因此，Huang等[3]提出了經(jīng)驗(yàn)?zāi)B(tài)分解（EMD），此方法具有自適應(yīng)分解特性，處理非線性和非平穩(wěn)信號(hào)具有較高效率。董文智等[4]為解決提取軸承故障特征頻率問題，提出了基于集總經(jīng)驗(yàn)?zāi)B(tài)分解（EEMD）和奇異值差分譜理論的故障診斷方法。但經(jīng)驗(yàn)?zāi)B(tài)分解缺少嚴(yán)格的數(shù)學(xué)推導(dǎo)，同時(shí)存在模態(tài)混疊、端點(diǎn)效應(yīng)和受采樣頻率影響較大等缺點(diǎn)[5]?；诖?，Dragomiretskiy等[6]提出了變分模態(tài)分解(VMD)，該方法將信號(hào)分解轉(zhuǎn)化為非遞歸、變分形式的分解模態(tài)，具有可靠的理論基礎(chǔ)，表現(xiàn)出更好的噪聲魯棒性。Mohanty等[7]利用VMD提取軸承故障特征頻率，驗(yàn)證了該方法具有比EMD更好診斷效果。Wan等[8]利用改進(jìn)的快速譜峭度和VMD相結(jié)合的方法，對(duì)軸承出現(xiàn)的復(fù)合故障進(jìn)行診斷。劉尚坤等[9]以互信息準(zhǔn)則確定VMD分量個(gè)數(shù)，并與Teager能量算子結(jié)合對(duì)轉(zhuǎn)子發(fā)生局部碰擦和油膜失穩(wěn)進(jìn)行時(shí)頻故障診斷，但在設(shè)定互信息的閾值時(shí)缺乏理論依據(jù)。

針對(duì)滾動(dòng)軸承出現(xiàn)微弱故障時(shí)，振動(dòng)信號(hào)表現(xiàn)為調(diào)幅調(diào)頻特性，本文提出了基于VMD-SVD能量標(biāo)準(zhǔn)譜-Teager能量算子聯(lián)合診斷方法，證明了該方法能有效地提取滾動(dòng)軸承微弱故障的特征頻率。

1 基礎(chǔ)理論

1.1 變分模態(tài)分解（VMD）

VMD是由Dragomiretskiy等提出信號(hào)自適應(yīng)分解估計(jì)方法，其實(shí)質(zhì)是用多個(gè)維納濾波組對(duì)信號(hào)進(jìn)行濾波。VMD算法利用了EMD中本征模態(tài)函數(shù)(IMF)概念，但I(xiàn)MF被重新定義為AM-FM信號(hào)

式中：相位φk(t)非遞減，φk’(t)≥0；瞬時(shí)幅值A(chǔ)k(t)≥0；并且瞬時(shí)幅值A(chǔ)k(t)和瞬時(shí)頻率wk(t)=φ’k(t)相對(duì)于φk(t)變化緩慢，即在[t-δ,t+δ] (δ=2π∕φk′(t))范圍內(nèi)，uk(t)能視作為幅值為Ak(t)、頻率為wk(t)的諧波信號(hào)。

VMD是一個(gè)變分問題的求解過程，假設(shè)每個(gè)IMF為具有中心頻率的有限帶寬，其中心頻率和帶寬在迭代求解過程中不斷更新。變分問題可表示為求K個(gè)本征模態(tài)函數(shù)uk(t)，使得所有uk(t)的估計(jì)帶寬之和最小，并且各模態(tài)函數(shù)之和等于分解信號(hào)的約束條件。具體分解過程如下：

（1）將各模態(tài)函數(shù)uk(t)進(jìn)行希爾伯特(Hilbert)變換，求得解析函數(shù)

（2）各模態(tài)解析函數(shù)分別與預(yù)估中心頻率e-jωkt相乘，將各模態(tài)頻譜移至基帶

（3）利用H1高斯平滑估計(jì)解調(diào)信號(hào)的帶寬，即梯度L2范數(shù)平方根，構(gòu)造受約束變分問題

式中：uk={u1,…,uK}和ωk={ω1,…,ωK}分別為分解得到的各模態(tài)分量和中心頻率為所有IMF之和；K為迭代次數(shù)；δ(t)為沖激函數(shù)。

（4）為求解約束變分問題，須引入二次懲罰因子α和Lagrangian乘法算子λ(t（)當(dāng)信號(hào)中混有高斯白噪聲時(shí)，二次懲罰因子α可確保重構(gòu)信號(hào)的精度；而Lagrangian乘法算子確保約束條件的嚴(yán)格性）將其轉(zhuǎn)變?yōu)闊o約束變分問題，增廣Lagrangian表達(dá)式為

利用交替方向乘子算法(ADMM)對(duì)增廣Lagrangian表達(dá)式中{uk}、{wk}和λ交替迭代，尋找“鞍點(diǎn)”，獲得最優(yōu)解，即為各模態(tài)分量{uk}和中心頻率{wk}。

1.2 奇異值能量標(biāo)準(zhǔn)譜

對(duì)于一維離散時(shí)間序列X={x(1),…,x(N)}，構(gòu)造m×n階Hankel矩陣如下

矩 陣 中 ：1＜n＜N，m=N-n+1，延 時(shí) 值 為 1，則H∈Rm×m。

存在矩陣U∈Rm×m，V∈Rm×m使得

奇異值標(biāo)準(zhǔn)能量譜定義為[10]

信號(hào)中的各成分在奇異值能量標(biāo)準(zhǔn)譜中表現(xiàn)為：有用信號(hào)譜線幅值大而陡峭，噪聲譜線幅值小而平緩。因此，奇異值能量標(biāo)準(zhǔn)譜線出現(xiàn)的拐點(diǎn)，即為有用信號(hào)和噪聲的分界點(diǎn)。保留分界點(diǎn)及前k個(gè)奇異值，將其后奇異值置為零，得到降噪后矩陣H′

式中：{u1…,uk}為U的前k個(gè)列向量；{v1…,vk}為V的前k個(gè)列向量。

將矩陣H′中相應(yīng)位置的元素相加求平均值[11]，即為降噪后信號(hào)。

1.3 Teager能量算子

Teager能量算子(Teager Energy Operator，TEO)是Teager在研究非線性語音建模時(shí)提出的信號(hào)分析算法[12]，記作ψ。設(shè)有連續(xù)時(shí)間信號(hào)x(t)，則Teager能量算子定義為

式中：x˙(t)和x¨(t)分別為連續(xù)信號(hào)的一階和二階微分，ψ[x(t)] 能表達(dá)振動(dòng)信號(hào)的瞬時(shí)能量值。

對(duì)于離散信號(hào)x(n)，利用差分代替微分形式，則式(10)可以近似表示為

Teager能量算子對(duì)離散信號(hào)相鄰的采樣點(diǎn)進(jìn)行計(jì)算時(shí)對(duì)調(diào)制信號(hào)具有良好的時(shí)頻分辨率，能檢測出信號(hào)中的瞬態(tài)沖擊成分。

2 算法步驟及流程

針對(duì)滾動(dòng)軸承早期微弱故障信號(hào)常受噪聲干擾、不易提取特征頻率的問題。本文提出了基于VMD-SVD能量標(biāo)準(zhǔn)譜-Teager能量算子的軸承故障診斷方法，其具體步驟和流程如圖1所示。

(1)對(duì)原始信號(hào)采取進(jìn)行去除趨勢項(xiàng)處理，消除基線漂移的干擾；

(2)對(duì)預(yù)處理后的信號(hào)進(jìn)行VMD分解，設(shè)置初始模態(tài)數(shù)K=2，懲罰因子為2000以及帶寬為0。模態(tài)數(shù)K依次增加1，觀察不同模態(tài)數(shù)K下，各IMF分量的頻譜和中心頻率，確定過分解時(shí)的模態(tài)數(shù)K，將其減1即為最優(yōu)模態(tài)數(shù)K；

(3)計(jì)算最優(yōu)模態(tài)數(shù)K下各IMF分量的峭度和相關(guān)系數(shù)，選取包含故障狀態(tài)信息的敏感IMF分量；

(4)構(gòu)造敏感IMF分量的Hankel矩陣，并進(jìn)行SVD分解。利用奇異值能量標(biāo)準(zhǔn)譜，確定拐點(diǎn)，選擇有效奇異值個(gè)數(shù)，對(duì)信號(hào)進(jìn)行降噪重構(gòu)，提高信噪比；

(5)求解重構(gòu)信號(hào)的Teager能量算子并進(jìn)行FFT變換即得到包絡(luò)譜；

(6)提取微弱故障特征頻率，確定故障原因。

圖1 具體步驟和流程

3 實(shí)驗(yàn)驗(yàn)證及分析

3.1 實(shí)驗(yàn)裝置

為了驗(yàn)證本文提出方法的有效性和準(zhǔn)確性，將該聯(lián)合方法應(yīng)用于實(shí)際的滾動(dòng)軸承故障診斷。本文所使用的數(shù)據(jù)來源于美國Case Western Reserve大學(xué)的軸承中心數(shù)據(jù)庫[13]。該試驗(yàn)平臺(tái)主要由1臺(tái)1.5 kW驅(qū)動(dòng)電機(jī)、1個(gè)扭矩傳感器、1臺(tái)測功儀和電機(jī)控制單元組成，試驗(yàn)平臺(tái)如圖2所示。

圖2 試驗(yàn)測試平臺(tái)

驅(qū)動(dòng)端軸承是深溝球軸承，型號(hào)為6205-2-RS JEM SKF，具體規(guī)格參數(shù)如表1所示。

采用電火花技術(shù)對(duì)電機(jī)軸承的內(nèi)圈、外圈和滾動(dòng)體進(jìn)行點(diǎn)蝕來模擬常見故障，采集電機(jī)驅(qū)動(dòng)端、風(fēng)扇端和基座在不同工況下的振動(dòng)加速度信號(hào)。

3.2 滾動(dòng)軸承內(nèi)圈單點(diǎn)故障分析

試驗(yàn)中設(shè)計(jì)4種故障尺寸，分別為0.007、0.014、0.021和0.028 inch(其中1 inch=25.4 mm)；設(shè)定4種載荷，分別為0、1、2和3 hp(1 hp=0.75 kW)。選取驅(qū)動(dòng)端振動(dòng)信號(hào)，其中驅(qū)動(dòng)軸轉(zhuǎn)速為1772 r/min，負(fù)載為1 hp，故障為微弱尺寸0.007 inch，采樣頻率fs=12kHz，信號(hào)長度為6000個(gè)采樣點(diǎn)。滾動(dòng)軸承內(nèi)圈故障理論公式為

表1 滾動(dòng)軸承結(jié)構(gòu)參數(shù)

式中：fi為內(nèi)圈故障頻率，fn為轉(zhuǎn)頻，n為滾動(dòng)體數(shù)量，d為滾動(dòng)體直徑，D為節(jié)徑，α滾動(dòng)體接觸角。

根據(jù)表1中參數(shù)和式(12)可求得：fn=29.5 Hz，fi=159.7 Hz

圖3(a)、圖3(b)和圖3(c)分別為內(nèi)圈單點(diǎn)微弱故障振動(dòng)加速度信號(hào)、頻譜和局部頻譜圖，從圖3(a)中可以看出，當(dāng)軸承出現(xiàn)故障時(shí)有明顯的沖擊成分和強(qiáng)背景噪聲干擾；而圖3(c)中存在較多峰值點(diǎn)，無法有效確定故障類型。

利用本文提出的方法分析故障信號(hào)，首先對(duì)原始數(shù)據(jù)采取去趨勢項(xiàng)預(yù)處理，消除基線漂移。對(duì)預(yù)處理后的信號(hào)進(jìn)行VMD分解時(shí)應(yīng)先確定模態(tài)數(shù)K，通過觀察不同模態(tài)數(shù)K下各IMF的中心頻率fc及頻譜是否出現(xiàn)重疊來確定最優(yōu)的模態(tài)數(shù)K。表2給出了不同模態(tài)數(shù)K時(shí)各IMF的中心頻率fc。

從表2可以看出：從K=5開始，求得的IMF之間出現(xiàn)了中心頻率相近情況。例如，當(dāng)K=5時(shí)，IMF4和IMF5對(duì)應(yīng)的中心頻率僅相差262 Hz，并通過計(jì)算這兩個(gè)IMF分量的頻譜，可以判斷取K=5時(shí)出現(xiàn)了過分解。因此，最優(yōu)的模態(tài)數(shù)K=4。對(duì)內(nèi)圈單點(diǎn)故障信號(hào)取4個(gè)模態(tài)分量進(jìn)行VMD分解，得到4個(gè)IMF及其相對(duì)應(yīng)的頻譜，如圖4(a)和圖4(b)所示。

表2 不同K值時(shí)各IMF的中心頻率

圖3 內(nèi)圈微弱故障振動(dòng)信號(hào)與頻譜

圖4 各IMF分量的時(shí)域波形與頻譜

觀察圖4(a)，可以發(fā)現(xiàn)各IMF分量時(shí)域波形均存在明顯沖擊成分，無法直接選取敏感IMF分量；觀察圖4(b)，可以發(fā)現(xiàn)各IMF圍繞其中心頻率分布，未出現(xiàn)模態(tài)重疊現(xiàn)象。

計(jì)算各IMF峭度和與預(yù)處理后信號(hào)之間的相關(guān)系數(shù)，結(jié)果如表3所示。

表3 各IMF峭度與相關(guān)系數(shù)（K=4）

從表3可以看出IMF3和IMF4的峭度值和相關(guān)系數(shù)均較大，包含的故障信息也更豐富。因此選取IMF3和IMF4對(duì)應(yīng)位置的元素進(jìn)行相加，可得到重構(gòu)信號(hào)u(t)，其時(shí)域波形和頻譜如圖5(a)和圖5(b)所示。

圖5 重構(gòu)振動(dòng)加速度信號(hào)與頻譜

構(gòu)造信號(hào)u(t)的Hankel矩陣并進(jìn)行SVD分解，求得的奇異值和奇異值能量標(biāo)準(zhǔn)譜分別如圖6(a)和圖6(b)所示。

從圖6(b)中可以看到，振動(dòng)加速度信號(hào)沖擊成分主要集中在前面的奇異值能量譜線，當(dāng)r＞18時(shí)，Er＜0.1。為了在保證故障信息不丟失同時(shí)最大程度降噪，取前18個(gè)奇異值對(duì)信號(hào)進(jìn)行降噪重構(gòu)，得到信號(hào)v(t)，圖7(a)和圖(b)為其時(shí)域和頻譜。

將圖7(a)和圖7(b)與圖3(a)和圖3(b)相比，可以看出降噪后時(shí)域波形毛刺明顯減少，頻譜上的峰值易于區(qū)分，便于后續(xù)處理。

對(duì)降噪信號(hào)v(t)進(jìn)行Teager能量算子解調(diào)，圖8(a)和圖8(b)分別為Teager能量算子時(shí)域和頻譜。

從圖8(b)中可以看出在30 Hz、58 Hz和88 Hz處存在明顯峰值，這些數(shù)值的最大公約數(shù)近似為30，因此可以判斷它們分別對(duì)應(yīng)驅(qū)動(dòng)軸的1倍、2倍和3倍轉(zhuǎn)頻。由于滾動(dòng)軸承在運(yùn)轉(zhuǎn)過程中滾動(dòng)體和滾道之間會(huì)發(fā)生相對(duì)滑動(dòng)以及頻率分辨率的影響，使得提取的轉(zhuǎn)頻和理論轉(zhuǎn)頻（29.5 Hz）有一定的差別。同時(shí)可以較容易看到160 Hz、318 Hz、478 Hz、638 Hz、798 Hz和956 Hz存在峰值，且這些數(shù)值的最大公約數(shù)近似為160 Hz（內(nèi)圈故障特征頻率為159.7 Hz），因此可以判斷它們分別對(duì)應(yīng)1到6階內(nèi)圈故障特征諧波頻率。在圖8(b)中的218 Hz、260 Hz等處出現(xiàn)峰值，其對(duì)應(yīng)各階故障特征頻率的調(diào)制邊頻帶。

圖6 奇異值大小與能量標(biāo)準(zhǔn)譜

圖7 降噪信號(hào)波形與頻譜

圖8 Teager能量算子波形與包絡(luò)譜

圖9 根據(jù)EEMD分解IMF分量的時(shí)域波形與包絡(luò)譜

3.3 與基于EEMD方法對(duì)比

采用EEMD算法對(duì)與上節(jié)相同的故障振動(dòng)信號(hào)進(jìn)行分解，得到11個(gè)IMF分量和1個(gè)趨勢項(xiàng)。根據(jù)EEMD算法的特性[5]，取前3個(gè)包含故障信息較多的IMF分量并進(jìn)行Hilbert變換解調(diào)，各IMF分量時(shí)域波形和包絡(luò)譜如圖9(a)和圖9(b)所示。

從圖9(b)能提取1倍轉(zhuǎn)頻30 Hz、2倍轉(zhuǎn)頻58 Hz和故障特征頻率160 Hz，但各IMF分量的包絡(luò)譜差別不太明顯，即驗(yàn)證了EEMD分解存在嚴(yán)重的模態(tài)混疊現(xiàn)象，同時(shí)在0 Hz附近出現(xiàn)了端點(diǎn)效應(yīng)，不利于故障特征頻率的提取。

根據(jù)本文提出的方法，在診斷過程中利用VMD-SVD能量標(biāo)準(zhǔn)譜對(duì)預(yù)處理后信號(hào)進(jìn)行最大限度的降噪重構(gòu)，同時(shí)采取Teager能量算子增強(qiáng)故障沖擊特征，最后對(duì)得到的離散時(shí)間周期信號(hào)進(jìn)行傅里葉變換，能提取故障特征頻率及其諧波頻率。與EEMD方法相比，本文提出的方法可以提取故障特征頻率及其高階次諧波頻率，利用最大公約數(shù)能準(zhǔn)確提取故障特征頻率，有效地避免了模態(tài)混疊和端點(diǎn)效應(yīng)所導(dǎo)致的不利于提取故障特征頻率的問題。

4 結(jié)語

(1)將文中所提方法應(yīng)用于滾動(dòng)軸承內(nèi)圈微弱故障信號(hào)診斷，根據(jù)VMD分解得到各IMF分量中心頻率及頻譜，能有效確定最優(yōu)模態(tài)數(shù)K；利用峭度和相關(guān)系數(shù)指標(biāo)能準(zhǔn)確選取敏感IMF分量。

(2)利用SVD能量標(biāo)準(zhǔn)譜能夠最大限度抑制噪聲，對(duì)降噪后的信號(hào)采取Teager能量算子解調(diào)能增

?強(qiáng)故障沖擊成分，便于提取轉(zhuǎn)頻及高階故障特征頻率。

(3)與EEMD分解方法相比，本文所提方法避免了模態(tài)混疊和端點(diǎn)效應(yīng)所導(dǎo)致的不利于提取故障特征頻率的問題。