基于參數(shù)優(yōu)化的MCKD的滾動軸承早期故障診斷

潘昕怡，岳建海

(北京交通大學 機械與電子控制工程學院，北京100044)

滾動軸承是旋轉(zhuǎn)機械的關鍵零部件，以振動信號為基礎的診斷方法，常用于滾動軸承的故障的檢測[1]。當滾動軸承出現(xiàn)早期故障時，由于受到工作環(huán)境噪聲干擾、振動傳輸路徑復雜多變以及信號衰減等因素的影響，故障沖擊信號特征一般十分微弱，因此對于滾動軸承早期故障特征的識別和提取十分困難。為了準確、有效地提取故障特征，需要降低故障信號噪聲，突出信號中的沖擊成分。

近年來，不少學者針對此問題進行了研究。蘇文勝等[2]針對譜峭度法對滾動軸承早期故障診斷效果不明顯的問題，提出了基于經(jīng)驗模態(tài)分解(Empirical Mode Decomposition，EMD)降噪與譜峭度法結(jié)合的方法，實現(xiàn)了滾動軸承早期故障診斷。馬倫等[3]利用最小Shannon熵優(yōu)化的Morlet小波，結(jié)合Morlet小波的降噪原理為滾動軸承早期故障診斷提供了新方法。馬增強等[4]將變分模態(tài)分解(Variational Mode Decomposition，VMD)與Teager 能量算子相結(jié)合，實現(xiàn)了軸承故障的精確診斷。

Donald 等[5]在最小解卷積的基礎上提出了最大相關峭度解卷積(Maximum Correlated Kurtosis Deconvolution，MCKD)，其以相關峭度最大化作為評價指標，通過迭代突出信號中的脈沖沖擊成分。文獻[6]將經(jīng)最大相關峭度解卷積后的信號與1.5維譜相結(jié)合判斷故障類型，有效地提取了早期故障。文獻[7]利用VMD對經(jīng)MCKD增強后的信號進行重構，重構后的信號經(jīng)包絡解調(diào)后可有效地提取故障特征頻率。

但MCKD 算法的移位數(shù)及濾波器長度均依靠人為選取，其準確性及適用性無法得到保證。文獻[8]利用單步長網(wǎng)格搜索法對濾波器長度和解卷積周期進行了優(yōu)化。文獻[9]利用粒子群優(yōu)化算法對MCKD算法中3個參數(shù)進行了優(yōu)化。因此本文提出了一種基于參數(shù)優(yōu)化的最大相關峭度解卷積算法。首先針對滾動軸承的實際工況，利用仿真信號對移位數(shù)M的選取進行了討論，選出了最優(yōu)移位數(shù)M。以形態(tài)能量熵作為評價函數(shù)通過分階段網(wǎng)格搜索法實現(xiàn)了濾波器長度L的自適選取，可在獲得最優(yōu)參數(shù)的同時提高運算效率，更有利于工程實際的應用。通過實測信號驗證，該方法可有效地增強被噪聲淹沒的早期故障信號，實現(xiàn)滾動軸承的故障診斷。

1 最大相關峭度解卷積

MCKD 算法的本質(zhì)是尋找一個FIR 濾波器，以相關峭度作為評價指標，使濾波后的信號的峭度值最大，突出信號中被噪聲淹沒的周期性沖擊。

信號y的相關峭度表達式為：

式中：M為移位數(shù)；T為故障信號周期；N為輸入信號樣本數(shù)。

逆濾波器的表達式為：

式 中：=[f1f2…fL]T為濾波器系數(shù)；L為濾波器長度。

MCKD的目標函數(shù)為：

為尋找最優(yōu)濾波器以使得CKM(T)最大，對式(3)求導得：

將濾波器結(jié)果以矩陣的形式表示為

式中：

將濾波器系數(shù)f代入式(2)，可得實測信號x的解卷積信號y。

最大相關峭度解卷積算法的迭代過程如下：

(1) 確定解卷積周期T、濾波器長度L、移位數(shù)M；

(2)計算輸入信號x的XT，X0X0T，(X0X0T)-1；

(3)計算經(jīng)濾波后的輸出信號y；

(4)根據(jù)y計算；

(5)更新濾波器系數(shù)

(6)如果濾波前后的迭代誤差ΔCKM(T)＜ε則停止迭代，否則回到步驟(3)繼續(xù)迭代。

2 最大相關峭度解卷積參數(shù)的確定

使用最大相關峭度解卷積算法前需確定解卷積周期T、濾波器長度L、移位數(shù)M3個參數(shù)。其中，周期T由采樣頻率和軸承故障特征頻率確定，表達式為：

式中：fs為采樣頻率；fi為軸承故障特征頻率，由軸承尺寸及實際轉(zhuǎn)速決定。

2.1 移位數(shù)M的討論

移位數(shù)M影響經(jīng)解卷積處理后信號的脈沖個數(shù)，一般選取1～7，當M＞7 時由于迭代方法超出浮點指數(shù)范圍導致精度降低。

文獻[5]中利用周期性尖脈沖信號進行仿真，隨移位數(shù)M的增加解卷積信號脈沖個數(shù)增加，從而提高算法故障檢測能力，因此選用M=7。但周期性尖脈沖信號各脈沖間隔相同，而實際軸承故障信號沖擊間隔受運行工況的影響存在微小波動，因此其仿真結(jié)果不適用于實際軸承故障診斷。本文將對仿真信號引入一個微小隨機波動τ，采用如下信號對移位數(shù)M的選取進行討論：

式中：h(t)為周期性沖擊；Ai為沖擊信號幅值；n(t)為高斯白噪聲；τi為第i次沖擊的微小隨機波動；轉(zhuǎn)頻fr=20 Hz；系統(tǒng)固有頻率fn=4 000 Hz；衰減指數(shù)C=800；幅值A0=0.5。

設置采樣頻率fs=12 000 Hz，故障特征頻率fi=160 Hz，采樣點數(shù)8 000。取濾波器長度L=200，解卷積周期T=fs/fi=75，改變移位數(shù)M，仿真結(jié)果如圖1所示。當M取l時，周期性沖擊特征明顯，解卷積效果最佳；隨著M的增加解卷積效果不佳，基本失去意義。

圖1 不同移位數(shù)的MCKD結(jié)果

仿真結(jié)果表明，MCKD 算法應用于滾動軸承故障診斷時移位數(shù)取1 效果最佳，當移位數(shù)取值增加后不但解卷積結(jié)果不佳，同時增加了計算復雜程度，因此本文最終取M=1。

2.2 濾波器長度L的自適選取

網(wǎng)格搜索法[10]是將目標參數(shù)歷經(jīng)在一定取值范圍內(nèi)所劃分的所有網(wǎng)格，并計算各網(wǎng)格的評價函數(shù)值。

對于非線性規(guī)劃問題網(wǎng)格搜索法可表示為：

式中：f(x)為目標函數(shù)；X=(X1，X2，…，Xn)為n個自變量；Xs、Xf分別為的下限和上限。

將 區(qū) 間[Xs，Xf] 劃分為N段，搜索步長t=(Xf-Xs)/N。從Xs開始，每次增加t至下一點Xs+t，直至歷經(jīng)整個區(qū)間形成全部網(wǎng)格，計算所有網(wǎng)格點的評價函數(shù)值，確定最優(yōu)評價函數(shù)值對應的點為最優(yōu)點。

形態(tài)能量熵[11](Morphological Energy Entropy，MEE)是指在對信號進行多尺度形態(tài)學分解后，不同尺度的信號分量的能量分布，可以反映信號分布的復雜程度。多尺度形態(tài)學的理論基礎是數(shù)學形態(tài)學，其中多尺度形態(tài)開、閉運算是兩種基礎運算，其計算公式如式(9)和式(10)所示：

式中：g代表單位結(jié)構元素；n為分析尺度；°代表開算子，·代表閉算子；Θ、⊕分別代表腐蝕算子和膨脹算子，離散函數(shù)f(n)(n=0，1，…，N-1)在結(jié)構元素g(m)(m=0，1，…，M-1)下的腐蝕、膨脹運算定義為：

在數(shù)學形態(tài)學的基礎上，經(jīng)多尺度形態(tài)學分解后的各尺度分量如式(12)和式(13)所示：

式中：di為在第i尺度下的信號分量；hi為不同尺度下的組合形態(tài)濾波器；g代表單位結(jié)構元素；·代表閉算子，°代表開算子。

得到不同尺度信號分量di(i=1，2，…，n)后，計算各尺度的信號能量Ei，建立能量譜E=[E1，E2，…，En]后計算信號的形態(tài)能量熵，如式(14)所示：

形態(tài)能量熵反映了信號能量分布的復雜程度，當信號能量分布越簡單，形態(tài)能量熵越??；當信號內(nèi)包含尺度較多，信號越復雜，形態(tài)能量熵越大[12]。因此，本文將形態(tài)能量熵作為網(wǎng)格搜索法的評價函數(shù)，以經(jīng)解卷積后信號的最小MEE值作為尋優(yōu)目標，進行濾波器長度的自適選取。

3 診斷流程

本文根據(jù)滾動軸承故障仿真信號討論的最佳移位數(shù)M值，及利用以經(jīng)解卷積后信號的形態(tài)能量熵作為評價函數(shù)的網(wǎng)格搜索法搜尋MCKD 算法的濾波器長度，提出了一種基于參數(shù)優(yōu)化的最大相關峭度解卷積的滾動軸承早期故障診斷方法，具體診斷流程如下：

(1)根據(jù)滾動軸承結(jié)構尺寸計算故障特征頻率，根據(jù)故障特征頻率及采樣頻率設置MCKD 算法解卷積周期T及其他相關參數(shù)。

(2)利用以解卷積后形態(tài)能量熵最小作為評價函數(shù)的網(wǎng)格搜索法，對MCKD算法的濾波器長度進行搜索。首先根據(jù)最佳濾波器的一般選取范圍設置第1 階段搜索范圍[20，500]及步長Lstep=10。得到最小MEE值所對應的濾波器長度bestL后設置新的搜索范圍[bestL-Lstep，bestL+Lstep]，更新搜索步長Lstep=1，得到最終的最優(yōu)濾波器長度。

(3) 根據(jù)尋優(yōu)結(jié)果設置MCKD 算法的參數(shù)，并用MCKD算法對故障信號進行預處理。

(4)對經(jīng)解卷積后的故障信號做Hilbert 包絡解調(diào)，獲得包絡譜，將故障特征頻率理論值與包絡譜中的突出頻率進行對比，判斷故障類型。

4 實測信號分析

為驗證本文提出方法的有效性，利用美國凱斯西儲大學軸承數(shù)據(jù)中心的滾動軸承故障數(shù)據(jù)進行分析。實驗采用軸承型號為SKF6205-2RS，其軸承結(jié)構參數(shù)如表1所示。利用電火花切割機在滾動軸承內(nèi)圈、外圈分別加工損傷直徑為0.177 8 mm的缺陷，以模擬內(nèi)圈、外圈故障。為更有效地模擬滾動軸承早期故障故障特征微弱的特點，分別對內(nèi)圈、外圈故障信號加入信噪比為-3 dB 和-7 dB 的高斯白噪聲。實驗中電機轉(zhuǎn)速為1 772 r/min，采樣頻率為12 kHz。根據(jù)滾動軸承故障特征頻率計算公式得內(nèi)圈故障特征頻率fi=159.93 Hz，外圈故障特征頻率fo=105.87 Hz。

表1 滾動軸承結(jié)構參數(shù)

圖2、圖3 分別為滾動軸承內(nèi)圈、外圈故障時域波形及包絡譜。故障信號的時域波形中沖擊特征被噪聲淹沒；在包絡譜中由于故障特征頻率的多倍頻不突出，難以作出準確判斷。

圖2 內(nèi)圈故障時域波形及包絡譜

圖3 外圈故障時域波形及包絡譜

4.1 內(nèi)圈故障分析

利用本文的方法對內(nèi)圈故障信號進行分析，圖4為兩次尋優(yōu)后得到的最佳濾波器長度L=157。根據(jù)第二節(jié)設置移位數(shù)M=1，解卷積周期T=75.03。

圖4 內(nèi)圈故障MEE隨L變化曲線

經(jīng)參數(shù)優(yōu)化后的MCKD結(jié)果如圖5所示。由圖5(a)可以看出，經(jīng)解卷積后信號的時域波形具有明顯的周期性脈沖成分；圖5(b)為經(jīng)解卷積后內(nèi)圈故障信號的包絡譜，可以清晰地看出內(nèi)圈故障特征頻率及多倍頻成分，且無其它突出故障特征成分，根據(jù)包絡譜基本可以判斷該軸承為內(nèi)圈故障。

圖5 MCKD降噪后的內(nèi)圈故障時域波形及包絡譜

4.2 外圈故障分析

利用本文的方法對外圈故障信號進行分析，圖6為兩次尋優(yōu)后得到的最佳濾波器長度L=56。根據(jù)第二節(jié)設置移位數(shù)M=1，解卷積周期T=113.35。

圖6 外圈故障MEE隨L變化曲線

經(jīng)參數(shù)優(yōu)化后的MCKD結(jié)果如圖7所示。由圖7(a)可以看出，經(jīng)解卷積后信號的時域波形具有明顯的周期性脈沖成分；由圖7(b)為經(jīng)解卷積后外圈故障信號的包絡譜，可以清晰地看出外圈故障特征頻率及多倍頻成分，且無其它突出故障特征成分，根據(jù)包絡譜基本可以判斷該軸承為外圈故障。

圖7 MCKD降噪后的外圈故障時域波形及包絡譜

4.3 濾波器長度優(yōu)化結(jié)果及運算效率對比

分別利用單步長網(wǎng)格搜索法和粒子群優(yōu)化算法對上節(jié)中內(nèi)圈、外圈故障降噪時的濾波器長度進行優(yōu)化。操作系統(tǒng)處理器為Inter(R) Core(TM) i5-5257U CPU@2.70 GHz，安裝內(nèi)存8.00 GB。利用軟件MATLAB2018b，分別設置單步長網(wǎng)格搜索法的濾波器長度搜索范圍[20，500]，搜索步長Lstep=1；粒子群優(yōu)化算法濾波器長度搜索范圍[20，500]，種群規(guī)模N=20，學習因子c1=c2=1.5，最大迭代次數(shù)20。與本文所提出方法的濾波器長度優(yōu)化結(jié)果及運算時間進行對比，如表2 所示。可以看出對于不同的故障類型3 種方法所獲得優(yōu)化結(jié)果相同，且在獲得效果相同的去噪結(jié)果時，本文提出的方法的運算時間小于單步長網(wǎng)格搜索法和粒子群優(yōu)化算法所用時間。

表2 3種優(yōu)化算法的L優(yōu)化結(jié)果及運算時間

5 結(jié)語

(1)本文利用接近滾動軸承實際故障工況的仿真信號，確定了MCKD算法中更適用于滾動軸承故障診斷的移位數(shù)。

(2)利用以形態(tài)能量熵為評價函數(shù)的網(wǎng)格搜索法，實現(xiàn)了MCKD 算法中濾波器長度的自適選取，且具有較高的運算效率，提高了算法中參數(shù)的準確性及適用性。

(3)通過軸承故障實測數(shù)據(jù)分析表明，基于參數(shù)優(yōu)化的最大相關峭度解卷積算法，可以實現(xiàn)滾動軸承早期故障的診斷，具有一定工業(yè)應用價值。