基于MOMEDA與Teager能量算子的滾動軸承故障診斷

祝小彥， 王永杰

(華北電力大學 能源動力與機械工程學院，河北 保定 071003)

旋轉機械設備中滾動軸承是重要的組成零部件之一，由于軸承故障造成機毀人亡的事故時有發生，因此對滾動軸承的狀態檢測和故障診斷有著重要意義。但因工況的復雜性所測振動信號具有非平穩性，且早期故障軸承振動信號中反應故障特征的沖擊成分較微弱，極易被噪聲覆蓋，較難直接由時域或頻域判斷故障類型[1]。因此，滾動軸承的早期故障診斷一直以來都是研究的難點和熱點。

滾動軸承故障源信號的傳遞過程可以看作是源信號與信道的一個線性卷積混合過程，提取故障原始沖擊信號則可以看作是一個解卷積的過程。從這個角度出發，Wiggins[2]提出了最小熵解卷積(Minimum Entropy Deconvolution, MED)并成功應用到了地震波的處理中。近年來，不少研究人員將MED解卷積算法引入到滾動軸承故障診斷中。柳玉昕等[3]通過將最小熵與Teager能量算子相結合并在包絡譜中找到了多個倍頻，但干擾頻率幅值較大，解卷積效果欠佳。冷軍發等[4]將MED應用于滾動軸承早期故障診斷只得到了基頻與2倍頻，且包絡譜中噪聲影響嚴重。王宏超等[5]提出了基于最小熵解卷積與稀疏分解的滾動軸承微弱故障特征提取方法成功提取到了3倍頻，但從包絡譜中可以看出噪聲的影響仍然較大。研究發現最小熵解卷積算法在處理滾動軸承故障信號時并不完全適用，主要包括以下三點：①MED定義故障信號第一個采樣點之前的幅值為零，造成了MED解卷積信號中偽沖擊成分的出現；②MED算法采用迭代的方式只能找到局部最優濾波器，并不一定是全局最優濾波器；③對于多點連續性沖擊的滾動軸承故障信號，MED解卷積信號中往往只有一個或者幾個沖擊成分，顯然這樣的解卷積結果并不能反映軸承故障時的真實情況。

針對以上問題，Mcdonald等[6]提出了一種新的解卷積方法——多點最優調整的最小熵解卷積(Multipoint Optimal Minimum Entropy Deconvolution Adjusted，MOMEDA)。MOMEDA針對旋轉機械故障信號的特點對解卷積的定義作了改進，并引入了目標向量和多點D-范數，不僅解決了最優濾波器的設計問題，而且實現了對軸承故障信號中的連續多點沖擊成分準確提取。

本文首次將多點最優調整的最小熵解卷積算法應用到滾動軸承早期故障診斷中，提出了基于MOMEDA和Teager能量算子的滾動軸承早期故障診斷方法。仿真和實測信號的分析表明該方法在滾動軸承早期故障的診斷中具有一定的有效性和實用性。

1 MOMEDA算法

假設y為故障軸承的一個沖擊信號，h為系統頻響函數，x為傳感器采集到的振動信號，e為隨機噪聲。則沖擊信號由信源到傳感器的傳輸過程可以近似地表達為

x=h*y+e

(1)

MOMEDA算法的核心部分是要通過非迭代的方式找到一個最優濾波器f，實現對原沖擊信號y的重構，并盡量削減噪聲對提取沖擊信號的影響。解卷積過程為

(2)

式中：k=1,2,…,N-L。MOMEDA算法針對旋轉機械故障信號中存在周期性沖擊的特點，在D-范數的基礎上提出了多點D-范數。即

多點D-范數

(3)

(4)

式中：t為目標向量，定義了解卷積目標沖擊成分的位置和權重。當目標向量t與原沖擊信號y完全契合時，解卷積效果達到最佳。此時多點D-范數取到最大值，與之對應的濾波器就是一組最優濾波器f。

式(4)的求解問題等價于求解方程

(5)

式中：f=f1,f2,f3,…,fL，t=t1,t2,t3,…,tN-L。

由式(2)、式(4)、式(5)可以求得

(6)

式中：k=1,2,…,N-L。

令X0=[M1，M2， …，Mk]，則式(6)簡記為

‖y‖-1X0t-‖y‖-3tTyX0y=0

(7)

整理得

(8)

(9)

取其特解作為一組最優濾波器，記為

(10)

2 Teager 能量算子

Teager能量算子是一種非線性差分算子，相對于傳統的信號能量定義，它增加了和頻率平方根的乘積，更能夠突出沖擊的瞬時特征[7]。

對于調幅調頻信號x(t)=a(t)cos[Φ(t)]，Teager能量算子Ψ的定義為

(11)

式(11)展開整理得

Ψ[x(t)]=[a(t)φ′(t)]2+a2(t)φ″(t)×
sin[2φ(t)]/2+cos2[φ(t)]Ψ[a(t)]

(12)

由于通常調制信號比載波信號變化慢得多，所以其幅頻值可近似視為常數處理[8]。令Ψ[a(t)]=0,φ″(t)=0，可得

Ψ[x′(t)]≈[a(t)φ′(t)]2=a2(t)ω2(t)

(13)

同理可得

Ψ[x′(t)]≈a2(t)ω4(t)

(14)

由式(13)和式(14)分別計算x(t)的瞬時幅值和瞬時相位

(15)

(16)

Teager能量算子對信號能量的定義中包含了信號的動能和勢能。與傳統能量的定義方式相比，Teager能量算子能夠較大程度上增強信號的幅值。然而研究發現，在實際應用中Teager能量算子對沖擊信號與噪聲信號并沒有太大的分辨能力，算法在增強沖擊信號幅值的同時，信號中噪聲的幅值也會相應的增加，使得在包絡譜中難以出現故障特征頻率。因此有必要對原始故障信號進行濾波處理，提高信噪比。

3 故障診斷流程

滾動軸承早期故障信號中，沖擊成分往往比較微弱，在強背景噪聲的影響下故障檢測更加困難。MOMEDA算法不僅能夠找到最優濾波器而且能夠準確提取出故障信號中周期性的沖擊信號，尤其適合于處理滾動軸承早期故障信號。然而筆者研究發現，直接對MOMEDA解卷積信號包絡分析其故障特征頻率并不十分突出，為此本文進一步提出利用Teager能量算子增強故障信號，從而使故障特征頻率在包絡譜中更加明顯。本文中嘗試將MOMEDA算法應用到滾動軸承故障診斷領域，并針對軸承早期故障信號的特點，提出了MOMEDA算法與Teager能量算子相結合的診斷方法。具體實現步驟：①利用MOMEDA算法提取滾動軸承早期故障信號中連續性周期沖擊成分；②借助于Teager能量算子對解卷積輸出信號中的沖擊成分進一步突出；③通過比較包絡譜中主導頻率成分與滾動軸承各元件故障特征頻率判斷故障軸承的故障類型。主要流程如圖1所示。

圖1 故障診斷流程圖Fig.1 The flowchart of fault diagnosis

4 仿真實驗

利用滾動軸承故障模型[9-11]對內圈故障時產生的沖擊信號進行模擬，并添加強烈白噪聲模擬軸承內圈早期故障信號。仿真信號為

(17)

式中：s(t)為周期性沖擊成分；n(t)為高斯白噪聲； 幅值為A0為0.5；τi為第i次沖擊相對于周期T的微小波動； 衰減系數C為800； 共振頻率fn為4 000 Hz； 轉頻fr為20 Hz，內圈故障特征頻率fi=1/T=110 Hz，隨機波動服從零均值正態分布，標準差為轉頻的0.5%，仿真信號信噪比為-13 dB，采樣頻率fs為12 000 Hz，采樣點數N為8 192。仿真信號時域圖、幅值譜、Teager能量算子包絡譜如圖2所示。

從圖2(a)中可以看出，由于在沖擊信號中添加了較重的噪聲信號，沖擊成分已經完全淹沒在噪聲信號中，時域波形中難以找到明顯的周期性沖擊特征。圖2(b)仿真信號頻譜中由于沖擊信號特征頻率及其倍頻與轉頻及其倍頻之間調制現象的發生，以及信號中強烈的噪聲信號的干擾，雖然能夠在頻譜中找到較突出的頻率成分但發現都是干擾頻率，并不能反映滾動軸承內圈的故障信息。圖2(c)是仿真信號直接經過Teager能量算子解調之后得到的Teager能量算子包絡譜，可以看出由于信號中強烈噪聲的影響包絡譜中仍然難以找到明顯的故障特征頻率成分。通過以上分析可知：軸承早期故障信號中微弱的沖擊信號常常會被背景噪聲所淹沒，單純地利用Teager能量算子包絡譜難以實現對滾動軸承內圈的準確識別，因此有必要在Teager能量算子包絡分析之前對故障信號進行濾波處理，提高信噪比。

圖2 仿真信號時域波形、頻譜及其Teager能量算子包絡譜Fig.2 Time domain waveform, spectrum and Teager energy operator envelope spectrum of simulation signal

圖3(a)是MOMEDA對仿真信號解卷積(濾波器階數為1 000，窗函數[1 000,1]，周期為109)后的時域波形，可以明顯看到有周期性的沖擊成分存在。圖3(b)中利用Teager能量算子對MOMEDA解卷積信號中的沖擊成分進一步增強并作包絡分析，可以清晰的看到軸承內圈故障特征頻率及其倍頻fi～9fi。

圖3 MOMEDA解卷積信號時域波形及其Teager能量算子包絡譜Fig.3 Time domain waveform and Teager energy operator envelope spectrum of MOMEDA demodulation

滾動軸承內圈早期故障信號模擬實驗表明：軸承故障發生早期故障信號中的沖擊成分被強烈的噪聲信號所淹沒，僅僅利用時域波形、幅值譜難以發現故障相關的頻率特征。而Teage能量算子在增強故障信號中的沖擊信號的同時噪聲信號也隨之增強，因此在Teage能量算子包絡譜中仍然難以發現明顯的故障特征頻率。通過故障仿真發現，MOMEDA算法能夠對故障信號中沖擊成分進行有效提取，通過Teager能量算子對故障信號中沖擊成分的進一步增強作包絡分析后可以清楚地看到軸承內圈故障特征頻率及其多個倍頻成分，證明了本文所提方法在提取故障信號中沖擊特征的有效性。

5 實測故障信號分析

5.1 人工植入軸承故障

實驗數據采用美國凱斯西儲大學電氣工程實驗室滾動軸承滾動體故障信號。實驗中所用驅動端滾動軸承型號為SKF 6205，具體參數如下表所示。采用電火花技術在軸承上加工單點凹痕模擬故障早期故障，人為加工的軸承損傷直徑分為0.177 8 mm、0.355 6 mm和0.533 4 mm，為體現本文所提方法的有效性，選用故障程度最輕的0.177 8 mm。傳感器采樣頻率為12 kHz，電動機轉速為1 797 r/min。滾動體擾動頻率為4.713 5 Hz，計算可知滾動體故障特征頻率fb為141.2 Hz。電動機轉頻fr為29.95 Hz。

表1 SKF 6205軸承結構參數

驅動端滾動軸承故障時域波形如圖4所示，由于原始故障信號中的沖擊成分被強噪聲淹沒，因此在時域波形中難以觀察到明顯而有規律性的沖擊特征。圖5原始信號頻譜圖中雖然能夠找到多條明顯譜線，但都不能反映滾動軸承的故障特征。對原始故障信號作Teage能量算子包絡分析，結果如圖6所示，包絡圖中只能觀察到電動機轉頻fr及2fr，并沒有滾動軸承故障特征頻率及其倍頻出現，所以難以判斷軸承故障類型。

圖4 原始故障信號時域波形Fig.4 Time domain waveform of fault signal

圖5 原始信號頻譜圖Fig.5 Spectrum of original fault signal

圖6 原始故障信號Teager能量算子包絡譜Fig.6 Teager energy operator Envelope spectrum of original fault signal

如圖7所示，MOMEDA解卷積信號(濾波器階數為1 500，窗函數[1 500,1]，周期為85.1)時域波形及其包絡譜。在圖7(a)MOMEDA解卷積信號中可以明顯看到信號中的周期性沖擊成分。圖7(b)中對MOMEDA解卷積信號直接作包絡，結果顯示，包絡譜中可以找到多個相關的滾動體故障特征頻率譜線，但并不十分突出。

圖7 MOMEDA解卷積信號時域波形及其包絡譜Fig.7 Time domain waveform and envelope spectrum of MOMEDA demodulation

圖8中利用Teager能量算子對MOMEDA解卷積信號中的沖擊成分作進一步增強，然后對增強解卷積信號作包絡分析，可以看到包絡譜中故障特征頻率更加明顯、突出，與圖7(b)相比故障特征頻率更明顯。由包絡譜中故障特征頻率可以斷定滾動軸承的滾動體發生了故障，分析結果與實際情況一致，證明了本文所提方法的有效性。

圖8 MOMEDA解卷積信號Teager能量算子包絡譜Fig.8 Teager energy operator envelope spectrum of MOMEDA demodulation

為了證明本文所提方法的優越性，利用MED算法(濾波器階數為1 500，迭代次數為100，迭代終止條件為0.001)提取原始故障信號中的沖擊成分，并利用Teager能量算子包絡譜對MED解卷積信號中的滾動體故障特征頻率進行包絡分析，結果如圖9所示。

如圖9(a)所示，MED解卷積信號中只有少量沖擊成分出現，解卷積效果并不理想，難以準確提取出旋轉機械故障信號中包含的連續性周期沖擊成分。圖9(b)中，Teager能量算子包絡譜中無法看到與滾動體故障相關的頻率成分。對比發現，本文所提方法在提取故障信號連續性周期沖擊特征方面有著明顯的優勢。

圖9 MED解卷積信號時域波形及其Teager能量算子包絡譜Fig.9 Time domain waveform and Teager energy operator envelope spectrum of MED demodulation

5.2 全壽命周期加速試驗信號

試驗分析數據來自NSFI/UCR智能維護系統中心的滾動軸承全壽命周期加速度試驗[12-14]，試驗臺布局如圖10所示。試驗臺轉軸上安裝有4個型號為ZA2115滾動軸承，其結構參數如表2所示。并利用彈性系統在軸承和轉軸上加載約2 671 N的徑向載荷，轉軸轉速為2 000 r/min。滾動軸承軸向和徑向分別安裝有353B33型高靈敏度ICP加速度傳感器。試驗過程中共進行3組試驗，利用NI DAQCard-6062采集卡采集試驗振動信號，采樣頻率為20 kHz。其中第二組試驗持續時間為164 h，共采集數據文件984個，采樣間隔為10 min，采樣點數為20 480。本文中數據選用第二組試驗數據對一號滾動軸承外圈故障進行分析。經計算，軸承外圈故障特征頻率fo為236.4 Hz。

圖10 試驗臺示意圖Fig.10 Schematic diagram of experiment platform

軸承型號軸承節徑滾動體直徑滾動體數接觸角ZA211571．5mm8．4mm1615．17°

如圖11所示，試驗軸承在0～9 790 min全壽命周期內，表征軸承故障程度的均方根值發生了顯著的變化。約5 100 min之后軸承故障開始有所增加，但波動幅度并不大，這一階段一般稱之為軸承故障的早期階段；而在7 020 min軸承振動信號的均方根值發生了突變，軸承故障進一步加劇，直到試驗最后均方根值達到最大軸承失效。試驗結束后在軸承外圈上發現了明顯的剝蝕現象。

為了驗證本方法對滾動軸承早期故障的有效性，本文選擇5 310 min時故障剛開始發生時的試驗數據

進行分析。圖12為早期故障信號的時域波形、頻譜圖及其Teager能量算子包絡譜。原始信號頻譜圖中雖然出現了許多較突出的頻率成分，但觀察發現這些頻率成分都是干擾頻率，并不能正確反映軸承故障。而在Teager能量算子包絡譜中也沒有出現較突出的外圈故障特征頻率。

圖11 軸承故障發展趨勢圖Fig.11 Bearing fault tendency chart

圖12 原始信號的時域波形、頻譜圖及其Teager能量算子包絡譜Fig.12 Time domain waveform, spectrum and Teager energy operator envelope spectrum of original signal

圖13為MED算法的解卷積信號(濾波器階數為1 500，迭代次數為100，迭代終止條件為0.001)及其Teager能量算子包絡譜，可以看出MED解卷積信號中并沒有出現連續性沖擊信號，在其Teager能量算子包絡譜中也沒有發現與軸承外圈故障特征頻率相關的譜線出現。

圖13 MED解卷積信號時域波形及其Teager能量算子包絡譜Fig.13 Time domain waveform and Teager energy operator envelope spectrum of MED demodulation

采用本文所提方法對軸承早期故障信號進行分析，圖14為MOMEDA解卷積信號及其包絡譜。與MED解卷積信號相比，圖14(a)MOMEDA解卷積信號(濾波器階數為1 500，窗函數[1 500,1]，周期為50.8)中出現了連續性周期沖擊成分。直接對MOMEDA解卷積信號進行包絡分析，結果如圖14(b)所示。

圖14 MOMEDA解卷積信號時域波形及其包絡譜Fig.14 Time domain waveform and envelope spectrum of MOMEDA demodulation

包絡譜中可以看到軸承外圈故障特征頻率及其倍頻fo～4fo，但故障頻率譜線并不十分突出。圖15中對MOMEDA解卷積信號作Teager能量算子包絡譜，可以看到包絡譜中出現了更加突出的外圈故障特征頻率及其倍頻，增強效果明顯。對比圖中主導頻率成分與軸承外圈故障特征頻率表明軸承外圈已經發生了故障，與試驗結果故障類型吻合。

圖15 MOMEDA解卷積信號的Teager能量信號包絡譜Fig.15 Teager energy operator envelope spectrum of MOMEDA demodulation

6 結 論

(1) 滾動軸承早期故障信號中信噪比較低，原始故障信號中難以找到規律的周期性沖擊成分，試驗證明MOMEDA算法確實能夠有效地提取出故障信號中的沖擊信號。與MED算法的解卷積效果相比，MOMEDA算法能夠更準確的提取出故障信號中連續性周期沖擊成分。

(2)試驗結果表明直接對MOMEDA解卷積信號作包絡時，包絡譜中也能夠出現故障特征頻率及其倍頻，但主導頻率成分并不突出。相比之下Teager能量算子包絡譜中故障特征頻率及其倍頻譜線則更加清晰明顯，試驗證明Teager能量算子能夠在一定程度上增強故障信號中的沖擊特征。

(3) 通過將MOMEDA算法與Teager能量算子相結合，能夠實現對故障信號中的故障特征的有效提取和分析，與單一方法的分析結果相比效果明顯，證明了本文所提方法的有效性。

[ 1 ] 余發軍,周鳳星,嚴保康.基于字典學習的軸承早期故障稀疏特征提取[J].振動與沖擊,2016, 35(6): 181-186.

YU Fajun, ZHOU Fengxing, YAN Baokang. Bearing initial fault feature extraction via sparse representation based on dictionary learning [J].Journal of Vibration and Shock, 2016, 35(6): 181-186.

[ 2 ] WIGGINS R A. Minimum entropy deconvolution[J]. Geoexploration, 1978, 9(16): 21- 35.

[ 3 ] 柳玉昕,石巖,王美俊,等. 基于最小熵解卷積和能量算子的滾動軸承故障診斷方法[J].組合機床與自動化加工技術, 2016(6): 114-117.

LIU Yuxin, SHI Yan, WANG Meijun, et al. A method of fault diagnosis for rolling bearing based on minimum entropy deconvolution and energy operator[J]. Modular Machine Tool & Automatic Manufacturing Technique, 2016(6): 114-117.

[ 4 ] 冷軍發,楊鑫,荊雙喜. 最小熵解卷積在滾動軸承早期故障診斷中的應用[J].機械傳動, 2015, 39(8): 189-192.

LENG Junfa, YANG Xin, JING Shuangxi. Application of minimum entropy deconvolution in the rolling element bearing fault diagnosis[J]. Journal of Mechanical Transmission, 2015, 39(8): 189-192.

[ 5 ] 王宏超, 陳進,董廣明. 基于最小熵解卷積與稀疏分解的滾動軸承微弱故障特征提取[J].機械工程學報, 2013, 49(1): 88-94.

WANG Hongchao, CHEN Jin, DONG Guangming. Fault diagnosis method for rolling bearing’s weak fault based on minimum entropy deconvolution and sparse decomposition[J]. Journal of Mechanical Engineering, 2013, 49(1): 88-94.

[ 6 ] MCDONALD G L. ZHAO Q. Multipoint optimal minimum entropy deconvolution and convolution fix: application to vibration fault detection [J]. Mechanical Systems and Signal Processing, 2017, 82: 461-477.

[ 7 ] 王天金,馮志鵬,郝如江,等.基于Teager能量算子的滾動軸承故障診斷研究[J].振動與沖擊,2012, 31(2): 1-5.

WANG Tainjin, FENG Zhipeng, HAO Rujiang, et al. Fault diagnosis of rolling element bearings based on Teager energy operator[J]. Journal of Vibration and Shock, 2012, 31(2): 1-5.

[ 8 ] 鞠萍華,秦樹人,趙玲,等. 基于LMD的能量算子解調方法及其在故障特征信號提取中的應用[J]. 振動與沖擊, 2011, 30(2): 1-4.

JU Pinghua, QIN Shuren, ZHAO Ling, et al. Energy operator demodulating approach based on LMD and its application in extracting characteristics of a fault signal [J]. Journal of Vibration and Shock, 2011, 30(2): 1-4.

[ 9 ] ANTONI J, BNNARDOT F, RAAD A, et al. Cyclostationary modeling of rotating machine vibration signals [J]. Mechanical Systems and Signal Processing, 2004, 18(6): 1285-1314.

[10] 唐貴基,王曉龍.自適應最大相關峭度解卷積方法及其在軸承早期故障診斷中的應用[J].中國電機工程學報, 2015, 35(6): 1436-1444.

TANG Guiji, WANG Xiaolong. Adaptive maximum correlated kurtosis deconvolution method and its application on incipient fault diagnosis of bearing[J]. Proceedings of The Chinese Society for Electrical Engineering, 2015, 35(6): 1436-1444.

[11] RANDALL R B, ANTONI J, CHOBSAARD S. The relationship between spectral correlation and envelope analysis in the diagnosis of bearing faults and other cyclostationary machine signals [J]. Mechanical Systems and Signal Processing, 2001, 15(5): 945-962.

[12] 馬倫,康建設,孟妍,等. 基于Morlet小波變換的滾動軸承早期故障特征提取研究[J]. 儀器儀表學報, 2013, 34(4): 920-926.

MA Lun, KANG Jianshe, MENG Yan, et al. Research on feature extraction of rolling bearing incipient fault based on morlet wavelet transform[J]. Chinese Journal of Scientific Instrument, 2013, 34(4): 920-926.

[13] QIU H, LEE J, LIN J, et al. Wavelet filter-based weak signature detection method and its application on rolling element bearing prognostics[J]. Journal of Sound and Vibration, 2006, 289(4/5): 1066-1090.

[14] 唐貴基, 王曉龍. 自適應最大相關峭度解卷積方法及其在軸承早期故障診斷中的應用[J].中國電機工程學報,2015,35(6): 1436-1444.

TANG Guiji, WANG Xiaolong. Adaptive maximum correlated kurtosis deconvolution method and its application on incipient fault diagnosis of bearing [J]. Proceedings of the Chinese Society for Electrical Engineering, 2015, 35(6): 1436-1444.