基于參數優化的MCKD的滾動軸承早期故障診斷

潘昕怡，岳建海

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

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

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

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

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

1 最大相關峭度解卷積

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

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

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

逆濾波器的表達式為：

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

MCKD的目標函數為：

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

將濾波器結果以矩陣的形式表示為

式中：

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

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

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

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

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

(4)根據y計算；

(5)更新濾波器系數

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

2 最大相關峭度解卷積參數的確定

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

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

2.1 移位數M的討論

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

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

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

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

圖1 不同移位數的MCKD結果

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

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

網格搜索法[10]是將目標參數歷經在一定取值范圍內所劃分的所有網格，并計算各網格的評價函數值。

對于非線性規劃問題網格搜索法可表示為：

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

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

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

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

在數學形態學的基礎上，經多尺度形態學分解后的各尺度分量如式(12)和式(13)所示：

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

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

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

3 診斷流程

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

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

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

(3) 根據尋優結果設置MCKD 算法的參數，并用MCKD算法對故障信號進行預處理。

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

4 實測信號分析

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

表1 滾動軸承結構參數

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

圖2 內圈故障時域波形及包絡譜

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

4.1 內圈故障分析

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

圖4 內圈故障MEE隨L變化曲線

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

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

4.2 外圈故障分析

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

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

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

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

4.3 濾波器長度優化結果及運算效率對比

分別利用單步長網格搜索法和粒子群優化算法對上節中內圈、外圈故障降噪時的濾波器長度進行優化。操作系統處理器為Inter(R) Core(TM) i5-5257U CPU@2.70 GHz，安裝內存8.00 GB。利用軟件MATLAB2018b，分別設置單步長網格搜索法的濾波器長度搜索范圍[20，500]，搜索步長Lstep=1；粒子群優化算法濾波器長度搜索范圍[20，500]，種群規模N=20，學習因子c1=c2=1.5，最大迭代次數20。與本文所提出方法的濾波器長度優化結果及運算時間進行對比，如表2 所示?？梢钥闯鰧τ诓煌墓收项愋? 種方法所獲得優化結果相同，且在獲得效果相同的去噪結果時，本文提出的方法的運算時間小于單步長網格搜索法和粒子群優化算法所用時間。

表2 3種優化算法的L優化結果及運算時間

5 結語

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

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

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