電機軸承故障的多尺度排列熵特征提取與GK識別*

周永強，卜文紹

(1.洛陽職業技術學院機電工程學院，河南 洛陽 471000；2.河南科技大學電氣工程學院，河南 洛陽 471023)

0 引言

軸承是旋轉機械設備中的關鍵部件，軸承異常輕則影響產品質量、導致生產過程停滯，重則引發機毀人亡等災難性事故[1]。研究軸承故障診斷技術可以準確把握故障的演變規律，對防止軸承故障引發的系列問題具有重要意義。

從軸承故障診斷機理的角度分類，故障診斷方法包括聲發射檢測法[2]、溫度檢測法[3]、噪聲分析法[4]、振動分析法[5]等，其中振動信號分析法應用最廣、精度最高，且可以實現實時監測。軸承振動信號中包含大量的噪聲，傳統的傅里葉變換沒有局部分析能力，難以有效提取故障信息。時頻分析法同時具有時域和頻域的局部分析能力，在故障診斷中得到廣泛應用，包括短時傅里葉變換、小波變換、EMD分解和局部均值分解等方法。文獻[6]使用短時傅里葉變換提取振動信號的時頻譜作為特征信號，并使用神經網絡進行故障模式識別，經驗證此方法具有較高的識別精度。文獻[7]研究了小波包與小波變換在軸承故障診斷中的區別和性能差異，結果表明小波包的降噪和識別性能優于小波變換。文獻[8]將振動信號使用EMD分解為若干模態分量，與原始信號組成二維特征圖，使用卷積神經網絡進行辨識，實現了多故障不同損傷程度的高準確率識別。文獻[9]提出了快速自適應局部分解方法的故障診斷方法，成功提取了軸承的故障特征。短時傅里葉變換由于窗口長度固定，使得時間和頻率分辨率無法兼顧；小波變換在時頻域上都具有多分辨率分析的能力，但是基小波的選擇需要經驗和先驗知識；EMD分解可以將信號自適應分解為多個本征模態，但是存在模態混疊和端點問題；局部均值分解也存在模態混疊問題。因此有必要研究具有自適應能力、能夠抑制模態混疊等問題的時頻分析法。

本文研究了軸承故障診斷問題，從故障特征提取和模式識別兩個方面入手。在故障特征提取方面，使用自適應局部迭代濾波對振動信號進行分解，選取與原始信號相關性較大的分量，計算分量信號的復合多尺度散布熵偏均值組成特征向量，經GK聚類，實現了故障類型的準確識別。

1 故障診斷方案

本文將故障診斷分為兩個方面，一是故障特征提取，二是故障模式識別，如圖1所示。整個故障診斷流程可以描述為4個步驟，其中前3步為故障特征提取過程，最后1步為故障模式識別過程。即：

步驟1：使用自適應局部迭代濾波將原始振動信號分解為多個I分量；

步驟2：計算各I分量與原始信號的相關性系數，由于主要的故障特征信息隱藏在相關性較大的I分量中，因此選擇相關性靠前的3個I分量做進一步分析；

步驟3：計算3個I分量的復合多尺度散布熵，組成3為故障特征向量，用于描述故障特征；

步驟4：使用GK聚類算法將特征向量進行分類，實現故障類型識別。

圖1 故障診斷方案

2 自適應局部迭代濾波

自適應局部迭代濾波(Adaptive Local Iteration Filter, ALIF)是在迭代濾波(Iteration Filter, IF)基礎上提出的，其自適應性體現在使用Fokker-Planck微分方程的基礎解系自適應構造了濾波函數。

2.1 迭代濾波

迭代濾波可以分為I分量提取和算法終止兩個方面進行說明，首先介紹I分量提取方法。原始振動信號記為x(t)，低通濾波函數記為ω(t)，兩者卷積得到滑動算子Γ(x)為：

(1)

振動信號x(t)減去滑動算子Γ(x)，得到波動算子k(x)為：

k(x)=x(t)-Γ(x)

(2)

一次I分量提取的終止條件定義為：

(3)

式中，θ為一次提取終止參數，ki,n為第i個I分量提取時第n次迭代的波動算子，‖ ‖2為2范數。

通過式(1)和式(2)的反復迭代，提取出I分量。算法迭代過程為：

7步驟1：首先計算濾波區間長度l，然后由式(1)計算滑動算子；

步驟2：對于第i個I分量的第n次迭代，波動算子為kn(x)=xn(t)-Γn(x)，此時進行I分量迭代終止判斷，若θ小于設定閾值則提取結束，kn(x)為第i個I分量(即Ii(t)=kn(x))；若θ大于設定閾值則令xn+1(t)=kn(x)，返回步驟1繼續迭代；

步驟3：計算殘余信號r(t)為r(t)=x(t)-∑Ii(t)，若r(t)表現出明顯的趨勢特征則算法結束；否則令x(t)=r(t)進行下一I分量的提取。

2.2 自適應局部迭代濾波

在迭代濾波中，低通濾波函數ω(t)根據個人經驗設定，不具備對振動信號x(t)自適應性設置。為了解決這個問題，文獻[10]使用Fokker-Planck微分方程構造了自適應濾波函數。

如果在區間[a,b]上g(x)、h(x)有：(1)g(a)=g(b)=0，且?x∈(a,b)，有g(x)>0；(2)h(a)<0

(4)

式中，α,β∈(0,1)為系數，p(x)為微分方程Pt(x)=0的解方程，且其形式滿足：

p(x)即為自適應低通濾波函數ω(t)，隨著窗口[a,b]在振動信號x(t)上的滑動，低通濾波函數ω(t)實現了自適應變化。

2.3 ALIF分解與EMD分解對比

為了對比ALIF與EMD分解的模態混疊情況，仿真出由調幅調頻信號、弦信號、隨機噪聲組成的混合信號，為：

(5)

式中，y(t)為混合信號，y1(t)為調幅調頻信號，y2(t)為弦信號，y3(t)為弦信號，n為采樣點數量。采樣時間為1 s，采樣頻率為1000 Hz。仿真信號如圖2所示。

圖2 仿真信號

使用ALIF濾波算法進行分解，得到I1分量、I2分量、趨勢項r(t)如圖3a所示。使用EMD對仿真信號進行分解，得到本征模態及趨勢項如圖3b所示。

(a) 自適應局部迭代濾波分解結果

(b) EMD分解結果圖3 兩種方法分解結果

圖3中自適應局部迭代濾波的I1分量對應原仿真信號中的調幅調頻信號y1(t)，I2分量對應原仿真信號中的弦信號y2(t)，原仿真信號中的弦信號y3(t)在分解過程中被濾除，由分解結果可知，自適應局部迭代濾波不存在模態混疊問題。EMD分解中IMF1分量與調幅調頻信號y1(t)頻率和幅值相近，IMF2與弦信號y2(t)相似，但是兩個分量均存在嚴重的模態混疊問題，已用紅色圓圈圈出。對比兩種方法的分解結果可知，ALIF濾波有效解決了模態混疊問題。

3 特征參數提取方法

3.1 散布熵

散布熵是用于描述時間序列復雜程度或分布混亂程度的指標，對于長度為N的序列x={xi}，其散布熵計算過程為：

步驟1：參數映射。將時間序列x映射為整數序列zi∈[1,2,…c]，映射方法為：

(6)

式中，yi為中間轉換變量，σ為正態分布的標準差，μ為正態分布均值，t為積分變量，R為四舍五入取整函數。

步驟2：抽取m維嵌入向量，并統計各散布模式數量。從{zi,i=1,2,…,N}中等間隔抽取m維嵌入向量，得：

zi,m={zi,zi+d,…,zi+(m-1)d},i=1,2,…,N-(m-1)d

式中，zi,m為m維嵌入向量，m為嵌入維度，d為延遲時間，i為嵌入向量編號。經過抽樣過程，得到了N-(m-1)d個嵌入向量，每個嵌入向量的散布模式記為πv0v1…vm-1，其中v0=zi，v1=zi+d，…，vm-1=zi+(m-1)d。統計每種散布模式的數量，記為Num(πv0v1…vm-1)。

步驟3：計算每種散布模式的概率，為：

式中，p(πv0v1…vm-1)為散布模式πv0v1…vm-1的概率。

步驟4：計算時間序列x的散布熵，為：

(7)

式中，DE(x)表示散布熵，cm表示所有可能存在散布模式的數量。

3.2 復合多尺度散布熵

散布熵是在單一尺度上計算序列的DE值，但是軸承振動信號的故障特征也可能隱藏在其他尺度上，因此有必要進行多尺度分析。序列x的復合多尺度散布熵計算方法為[11]：

步驟1：對于給定的時間序列x={xi}，經過映射后為{zi},i=1,2,…,N，其對尺度因子為τ的第k個粗?；蛄械牡趈個元素為：

(8)

式中，?」為向下取整函數。

步驟2：對于具有相同尺度因子τ的τ個粗粒序列，計算其散布熵的平均值即為復合多尺度散布熵CMDE(x)，為：

(9)

分析式(8)和式(9)可知，復合多尺度散布熵求取了時間序列在多個時間尺度上的DE值，能夠提取更多隱藏在振動信號中的特征信號。

3.3 復合多尺度散布熵偏均值

步驟1：記時間序列x在不同尺度下的CMDE值序列為CMDE(τ)={CMDE(1),CMDE(2),…,CMDE(τmax)}，則CMDE序列的偏斜度為：

(10)

式中，Ske(CMDE)為CMDE序列的偏斜度，VCMDE為CMDE序列的均值，MCMDE為CMDE序列的中值，SCMDE為CMDE序列的標準差。偏斜度的取值范圍為[-3,3]，當Ske=0時表示數據分布對稱，當Ske>0時表示數據分布右偏，當Ske<0時表示數據分布左偏，偏斜度絕對值越大則偏斜度越大。

步驟2：復合多尺度散布熵偏均值定義為：

(11)

式中，PM(CMDE)為CMDE序列的偏均值。

4 GK聚類原理及評價指標

4.1 GK聚類原理

對于任意一個待聚類序列S=(s1,s2,…,sn)，其聚類的目標函數定義為[12]：

(12)

式中，J(S,V,U)為聚類目標函數；m為模糊指數，一般取為1或2；Dij為樣本j與聚類中心i的馬氏距離，即：

(13)

式中，Zi是正定對稱矩陣，其定義為：

(14)

式中，Fi為模糊協方差矩陣。

對式(12)優化求解可得隸屬度矩陣U、聚類中心向量V分別為：

(15)

給出GK聚類的迭代過程為：

步驟1：設置參數初值，包括：模糊指數m，類數d，隸屬度矩陣初值U，迭代次數計數量l=0；

步驟2：根據式(15)更新聚類中心向量V=(v1,v2,vd)T，根據式(13)和式(14)更新距離Dij；

4.2 聚類評價指標

聚類效果一般使用分類系數和平均模糊熵進行評價，兩個評價參數的定義如下：

(16)

式中，PC為分類系數，其值越接近于1說明分類效果越好；CE為平均模糊熵，其值越接近于0說明分類效果越好。

5 實驗驗證

使用美國某大學軸承數據中心的數據對故障診斷方法進行驗證。軸承型號選擇SKF6205深溝球軸承，電機轉速選為1797 r/min，軸承狀態分為正常狀態、內圈故障、外圈故障和滾動體故障等4種，故障直徑為0.177 8 mm。數據采樣頻率為12 kHz，每個數據樣本長度為2048。每種狀態下選擇40組樣本，共4×40=160組樣本數據。以一組樣本數據為例，介紹故障診斷過程。原始數據如圖4所示。

(a) 正常狀態 (b) 滾動體故障

(c) 外環故障 (d) 內環故障圖4 原始數據

為了形成比較，同時使用EMD分解和ALIF濾波在時頻域分解原始數據。由于篇幅有限，在此僅給出ALIF濾波分解的前5個I分量和EMD分解的前5個IMF分量，如圖5所示。

(a) ALIF分解結果

(b) EMD分解結果圖5 分解結果

ALIF濾波分解與EMD分解的5個分量與原始數據的互相關系數如表1所示。

表1 各分量互相關系數

由表中各分量與原始數據的互相關系數可知，相關系數由分量1到分量5呈單調遞減，且前3個分量的相關系數較大，從第4個分量開始互相關系數開始小于0.1。以上數據說明原始數據的故障特征信息主要包含在前3個分量中，因此選擇前3個分量構造特征向量。

根據以上分析，一組原始數據經過分解，選擇前3個信號分量計算其復合多尺度散布熵偏均值作為特征向量，那么160組原始數據就相應得到160組3維的特征向量，由于數據量過大，在此不進行展示。將此160組3維特征向量輸入到GK模糊聚類中，聚類中心設置位d=4，模糊指數m=2，迭代終止閾值為ε=10-4。得到三維空間聚類結果如圖6所示。

(a) ALIF提取特征的聚類結果

(b) EMD提取特征的聚類結果圖6 聚類效果

對比圖6中兩種方法的聚類效果可知，ALIF提取的特征向量在聚類時類與類之間區分明顯，樣本分布極為緊湊，完全不存在混疊現象；且類與類中心之間的距離較大，說明類與類之間區分明顯，特征參數能夠極好地描述故障特征。EMD提取的特征向量聚類后，類與類之間存在邊界不清，或者交叉混疊現象，尤其是滾動體故障與內圈故障兩種情況；且樣本圍繞類中心分布較為分散，說明特征參數難以完全將軸承狀態區分開來。為了進一步比較聚類效果，計算兩個聚類結果的分類系數和平均模糊熵如表2所示。

表2 聚類參數

由表2中數據可知，使用ALIF提取特征的聚類分類系數大于EMD提取特征的聚類分類系數，使用ALIF提取特征的聚類平均模糊熵遠小于EMD提取特征的聚類平均模糊熵。以上數據也說明了ALIF提取特征的聚類效果優于EMD提取特征的聚類效果。

綜上所述，使用ALIF濾波進行分解，計算前3個分量的復合多尺度散布熵偏均值，并使用GK聚類進行模式識別的診斷方法更加有效，這是因為ALIF分解能夠自適應給出低通濾波函數，且不存在模態混疊問題，能夠較好的從原始數據中分解出特征分量。

6 結論

本文針對軸承故障診斷問題，從故障特征提取和模式識別兩個角度進行了研究。經過驗證得出以下結論：

(1)ALIF分解法能夠自適應給出低通濾波函數，且不存在模態混疊問題，能夠更加有效地提取特征信息。

(2)使用ALIF分解，計算分量復合多尺度散布熵偏均值，并使用GK聚類的方法，能夠完全識別出故障類型，類與類之間區分明顯，不存在混疊問題。