基于VMD_MPE和FCM聚類的變轉速工況下轉子不平衡故障診斷方法

鐘志賢， 馬李奕， 蔡忠侯， 段一戩， 陳金華

(桂林理工大學 機械與控制工程學院，廣西 桂林 541006)

轉子作為旋轉設備的核心部件，其正常運轉是工業生產的重要保證，一旦旋轉設備的轉子失衡，其零部件會受到額外的應力和變形，引起旋轉設備的振動，影響運行安全，還會造成產品質量下降，浪費能量，惡化工作環境等問題；因此，為了保證轉子在服役全壽命過程中平穩安全地運行，進行轉子不平衡故障的診斷尤為重要[1]。

目前，許多學者在轉子不平衡故障診斷領域進行了卓有成效的研究，文獻[2-7]采用了如嵌入式譜隨機有限元法、多尺度拉普拉斯特征映射法、等效載荷最小化法、多源異構信息融合深度信念網絡法、K-最近鄰分類器法以及決策樹等方法對轉子不平衡故障進行了理論分析與試驗驗證，為轉子不平衡故障的診斷提供了極具參考意義的指導。但文獻[2-7]的研究對象均為恒定轉速工況下的轉子不平衡故障，而在實際工業生產中，旋轉機械的轉速往往不是恒定的，需根據具體情況調整、改變轉子轉速以滿足相應工況，而變轉速問題也一直是故障診斷領域的難點問題[8]。

時至今日，許多學者在旋轉機械變轉速故障診斷領域取得了豐富的研究成果。Hu等[9]提出一種自適應無轉速階次分析方法，采用了一種新型基于動態路徑優化的山脊提取算法來估計瞬時頻率，仿真和試驗結果表明：該方法對噪聲具有較強的魯棒性，對變轉速工況下軸承故障的檢測是有效的。唐貴基等[10]提出一種基于奇異譜分解-希爾伯特變換(singular spectrum decomposition & Hilbert transform, SSD-HT)時頻階比跟蹤的轉子故障診斷方法來解決變轉速工況下轉子故障特征難提取的問題，運用奇異譜分解轉子故障振動信號，得到奇異譜分量，再用希爾伯特變換計算各奇異譜分量的瞬時頻率，以得到時頻分布，進而采用時頻分布中的轉頻信息對振動信號進行階比跟蹤分析，以提取直觀階次特征。最后，采用仿真與油膜和碰膜故障的試驗分析，驗證了該方法的有效性。郭俊[11]提出一種通過單一振動信號分析實現變轉速工況下永磁同步發電機軸承故障診斷的方法。該方法首先從振動信號中估計發電機的轉速，再用估計轉速對原始振動信號進行重采樣實現階次分析和故障診斷，降低了傳感器的需求，提高了故障診斷的效率。任勇[12]在分析轉速變化對峭度分類指標影響的基礎上，結合信號階譜特性提出階譜相關峭度(order spectrum corrlated kurtosis, OSCK)指標，并采用OSCK替換快速譜峭度圖中的譜峭度指標，提出了基于快速階譜相關峭圖(fast order spectrum kutogram, FOSCK)的信號故障敏感頻帶定位方法，實現了不同轉速運行工況下轉子、軸承和齒輪不同故障對應敏感頻帶的準確定位，有效消除了轉速變化對信號敏感頻帶評判指標帶來的影響。

文獻[9-12]采用了如自適應無轉速階次分析法，SSD-HT時頻階比跟蹤法，基于FOSCK的信號故障敏感頻帶定位法，基于瞬時轉速分析等方法，為旋轉機械變轉速故障診斷的進一步研究提供了極具價值的理論與試驗基礎，但依然存在如下問題：

(1)進行變轉速轉子斷條故障研究時，只側重于轉子斷條故障的定性判斷，未進行定量分析研究，奇異值分解(singular value decomposition, SVD)方法分解出的矩陣解釋性不強，且無理論基礎，不可研究信號在多重尺度上的特征信息；所提方法的實施效果與電機的先驗知識和故障診斷人員的經驗密切相關，故不具推廣性且只研究了轉子斷條一種故障類型，未涉及轉子不平衡故障；

(2)在進行瞬時轉速分析診斷法進行研究時，只進行了仿真分析，設定參數過于理想，未進行實際工況的試驗驗證。

(3)在進行變轉速工況下永磁同步發電機軸承故障診斷的研究時，只研究了發電機軸承的故障，并未對除軸承外的其他旋轉部件，例如轉子進行研究。

(4)在采用文獻中所提到的各方法進行振動信號處理時，都不可避免地遇到不自適應、模態混疊、端點效應等問題，從而影響最終的故障診斷效果。

為解決上述問題，本文提出一種基于變分模態分解(variational mode decomposition,VMD)與多尺度排列熵(multiscale permutation entropy,MPE)和模糊C均值(fuzzy C means,FCM)聚類的變轉速工況下轉子不平衡故障診斷的方法。首先，采用VMD非線性自適應信號處理法對轉子振動位移信號進行分解，得到K個本征模態分量(intrinsic mode function,IMF)，根據轉子不平衡故障時1×處振幅劇烈增加的現象[13]，從VMD所得的IMF頻譜圖中篩選出最能表征不平衡故障特征(1×處振幅最大)的IMF，然后采用MPE法對所篩選的IMF進行量化，將作為特征向量的多尺度排列熵輸入FCM，得到不同轉速工況的標準聚類中心，通過試驗數據從四種常用的模糊貼近算法中選取出最適合變轉速工況下轉子不平衡故障識別的算法，再用篩選出的模糊貼近算法計算待測試數據與各標準聚類中心的貼近度，進而實現變轉速工況下轉子不平衡故障的識別。最后，采用VMD_MPE-FCM法對轉子試驗臺變轉速工況下轉子不平衡故障位移信號數據進行了分析試驗，結果表明：VMD_MPE方法能夠精準地提取出變轉速工況下轉子不平衡位移信號的故障特征，FCM聚類方法能有效的識別出轉子不平衡故障。因此VMD_MPE-FCM法應用于轉子變轉速工況下不平衡故障的診斷是有效的，具有一定的理論意義。

1 VMD_MPE-FCM原理

VMD_MPE-FCM聚類結合方法的原理流程如圖1所示，對轉子振動位移進行處理的流程分為特征提取(VMD_MPE)和模式識別(FCM聚類)兩個部分。

圖1 VMD_MPE-FCM原理流程圖Fig.1 Flow chart of VMD_MPE-FCM

1.1 VMD原理

VMD算法中，重新定義IMF為一個調幅，調頻信號[14]

μk(t)=Ak(t)cos[φk(t)]

(1)

VMD的實質是為了尋找模態μk的集合，這些μk可以通過最小二乘法重構給定的輸入信號f(t)。同時，每個μk都限制在一個在線估計的中心頻率ωk內，即VMD方法建立的約束變分模型為

(2)

式中：{μk}={μ1,μ2，…,μk}為分解得到的K個模態分量；{ωk}={ω1,ω2,…,ωk}為各分量的頻率中心；δ(t)為脈沖函數；*為卷積符號。

建立該模型時，首先通過Hilbert變換得到μk(t)的解析信號，從而得到μk(t)的單邊頻譜，然后進行頻率混合，對各模態解析信號混合一個預先估計的中心頻率e-jωkt，將各模態頻譜變換到基頻帶上，最后計算解調信號梯度的平方L2的范數，估計出模態信號帶寬[15]。

為將約束問題轉化為非約束問題，引入代表每個模態初始中心約束強度的二次乘法因子α和拉格朗日乘法算子λ(t)，可得擴展的拉格朗日函數表達式為

(3)

(4)

(5)

VMD將原始信號分解為K個模態分量的具體步驟總結如下：

(6)

步驟4重復步驟2和步驟3，直到滿足迭代終止條件式(7)

(7)

式中，ε為判別精度，且ε>0。

步驟5輸出結果，得到n個模態分量。

1.2 MPE原理

(8)

(9)

式中：l為第l個重構分量，l=1,2,…,N-(m-1)τ；m為嵌入維數；τ為延遲時間。

(10)

S(r)=(l1,l2,···,lm)

(11)

式中,r=1,2,…，R且R≤m！，嵌入維數為m的時間重構序列共有m!種排列，符號序列S(r)是其中的一種排列。計算出每一種符號序列出現的概率Pr(r=1,2,…,R)后，用信息熵的形式定義不同符號序列的排列熵Hp(m)為

(12)

對Hp(m)進行歸一化處理，得

(13)

式中,Hp為歸一化處理后的排列熵值。其中0≤Hp≤1。其值越小，說明時間序列越規則平滑；反之，說明時間序列越復雜越有隨機性。采用多尺度排列熵的這一特性，可以有效的反映出轉子振動位移信號在多個尺度下的隨機性。

1.3 FCM原理

FCM法是一種基于目標函數的模糊聚類算法，其以極小化所有數據點與聚類中心的模糊隸屬度的加權和為目標，不斷修正聚類中心和分類矩陣到符合終止準則，將類似特征的數據樣本聚為一類[17]。

假定樣本數據集為D={d1,d2,…,dn}，把這些數據劃分為c類，即對應c個標準聚類中心為C={c1,c2,…,cn}，每個樣本j屬于某一類i的隸屬度為uij，由此可定義FCM目標函數及其約束條件為

(14)

(15)

式中,m為隸屬度因子。

采用拉格朗日乘數法將約束條件式(15)并入目標函數式(14)中，并把約束條件式(15)的所有j展開，則被處理過的目標函數可化為

(16)

式中,λ=[λ1,…，λn]為拉格朗日乘法算子，它將帶約束條件的極值問題轉化為無約束問題。

為了求解目標函數的極值，分別對式(16)中的uij和ci求導，使目標函數式為最小值，可得最終隸屬度uij的迭代公式為

(17)

最終聚類中心ci的迭代公式為

(18)

根據以上FCM法的基本原理，FCM聚類算法步驟如下。

步驟1確定分類數c，隸屬度因子m，以及迭代次數l=0；

步驟2初始化隸屬度，其中要求隸屬度uij的和為1；

步驟3根據式(18)計算樣本所有聚類中心C；

步驟4根據式(17)、式(18)更新目標函數J；

步驟5給定判別精度ξ>0，直到滿足條件，‖Jl+1-Jl‖<ξ，否則置l=l+1，根據聚類中心C返回計算隸屬度uij，再返回步驟3，直至滿足條件。

2 模糊模式識別

2.1 模糊模式識別直接方法

假定待識別對象的數據集為U={u1,u2,…,un}。當采用某一種模糊識別算法作用于診斷對象集U時，會識別出n個不同的故障類別，記識別出的故障類別數據集為Q={q1,q2,…,qn}；同時，也會產生一組對象uj隸屬于類別qi程度的參數，記為隸屬度μqi(uj)，如果對?uj∈U，有i∈{1,2,…,n}，使式(19)成立。

μqi(uj)=max{μq1(uj),μq2(uj),···,μqn(uj)}

(19)

即符合最大隸屬度原則，認為uj優先隸屬于qi。

上述方法被稱為故障診斷模糊模式識別的直接方法。雖可以直觀地實現模糊識別。但是，本文中轉子故障的特征是多元素的，而最大隸屬度原則無法進行多元素識別，所以模糊模式識別直接方法在此失效。

2.2 模糊模式識別間接方法

針對模糊識別直接方法不能識別多元素故障特征的特性，在此引入模糊識別間接方法—采用擇近原則和貼近度實現識別及分類[18]。

假定A={a1,a2,…,an}，B為A內的一個模糊子集，σ∈[0,1]為某種貼近度，式(20)成立

σ(B,ai)=max[σ(B,a1),···,σ(B,an)]

(20)

就認為B歸屬于集合A，此為擇近原則基本思想，而擇近原則最基礎的步驟為計算兩個模糊子集之間的貼近度。貼近度可用來表示模糊子集間彼此相近的程度，兩個模糊子集之間的貼近度越大，則相近程度越高，反之越小。

為了更好地表征待檢測數據與聚類中心的貼近度，引入了四種最為常見的模糊貼近算法，分別為：格貼近算法、最小最大貼近算法、海明貼近算法及歐幾里得貼近算法，原理分述如下：

(1)格貼近算法

(21)

式中：B·A為B與A的內積；B?A為B與A的外積。

(2)最小最大貼近算法

(22)

(3)海明貼近算法

(23)

(4)歐幾里得貼近算法

(24)

對給定的同一域中的兩個模糊子集，采用不同的模糊貼近算法所求出的貼近度不相同，可根據經驗和實際需要選擇合適的計算公式。

3 試驗分析與驗證

3.1 變轉速工況下轉子不平衡故障特征提取

為驗證本文提出的VMD_MPE-FCM方法的有效性，采用DHRMT轉子試驗臺進行變轉速工況下轉子不平衡故障試驗，試驗采用一個光電傳感器和兩個電渦流傳感器分別采集試驗時的轉速和振動位移，圖2為轉子試驗臺結構圖，圖3為轉子試驗臺，圖4、圖5分別為光電轉速傳感器和電渦流位移傳感器的安裝位置。

圖2 轉子試驗臺結構圖Fig.2 Schematic diagram of test system

圖3 轉子試驗臺Fig.3 The rotor experimental platform

圖4 光電轉速傳感器安裝位置Fig.4 Position of photoelectric speed sensor

圖5 電渦流位移傳感器安裝位置Fig.5 Position of eddy current displacement sensor

為模擬轉子不平衡故障，在試驗臺轉子圓盤45°位置的螺孔內安裝M5×16，質量為2 g的螺釘，如圖6所示。經現場動平衡測試系統測得：轉子動不平衡量為61.46 g·mm，位置為36°。通過圖7下方的轉子控制器，調整轉子轉速分別為600 r/min、1 200 r/min、1 800 r/min、2 400 r/min、3 000 r/min，經光電轉速傳感器測得的各轉速為：592 r/min、1 217 r/min、1 800 r/min、2 386 r/min、2 984 r/min。用圖7上方的信號采集儀，采樣頻率為2 000 Hz，采樣時間為1 s，采集每種轉速工況下的轉子不平衡故障數據各100組，共計500組。因篇幅限制，圖8～圖12為任意選取的各轉速中某一組數據的時域波形。

圖6 配重螺釘安裝位置Fig.6 Position of counterweight screw

圖7 轉子控制器與信號采集儀Fig.7 Rotor controller & signal acquisition instrument

圖8 600 r/min工況下不平衡故障信號時域波形Fig.8 Time-domain waveform of unbalance fault signals under 600 r/min

圖9 1 200 r/min工況下不平衡故障信號時域波形Fig.9 Time-domain waveform of unbalance fault signals under 1 200 r/min

圖10 1 800 r/min工況下不平衡故障信號時域波形Fig.10 Time-domain waveform of unbalance fault signals under 1 800 r/min

圖11 2 400 r/min工況下不平衡故障信號時域波形Fig.11 Time-domain waveform of unbalance fault signals under 2 400 r/min

圖12 3 000 r/min工況下不平衡故障信號時域波形Fig.12 Time-domain waveform of unbalance fault signals under 3 000 r/min

VMD算法中，K值的選取會影響分解的結果，如K值過小，則波形信息會顯示不全；如K值過大，則會出現過分解現象。本文通過觀察各模態中心頻率的方法確定K值，若出現中心頻率相近的模態分量，則認為出現VMD過分解現象。規定兩個模態中心頻率差值的絕對值小于等于中心頻率較大值的10%時為過于接近，即認定為過分解。本文以600 r/min工況下不平衡故障信號為例，經VMD后，不同K值各模態分量的中心頻率見表1。

表1 不同K值對應的中心頻率Tab.1 Center frequencies corresponding to different K

由表1知，當K為7時，第三模態的中心頻率為35 Hz，第四模態的中心頻率為38 Hz，兩個模態中心頻率差的絕對值為3 Hz，小于中心頻率較大值的10%，為3.8 Hz，故認定為過分解，即兩個模態之間的中心頻率開始出現中心頻率相近的模態分量。所以，600 r/min工況下的K值為6。用同樣的方法對其他轉速的信號進行分析，可以得到1 200 r/min、1 800 r/min、2 400 r/min、3 000 r/min工況下的不平衡故障信號的K值分別為6、6、7、8。懲罰因子α取默認值2 000。當K取對應值時，用前文從每種轉速信號中選取的數據繪制VMD時域圖，分別如圖13～圖17所示，各圖中，f為原始信號的時域圖，IMFi(i=1,2,…,n)為分解得到的各IMF時域圖。

圖13 600 r/min工況下不平衡故障信號分解Fig.13 Decomposed unbalance signals under 600 r/min

圖14 1 200 r/min工況下不平衡故障信號分解Fig.14 Decomposed unbalance signals under 1 200 r/min

圖15 1 800 r/min工況下不平衡故障信號分解Fig.15 Decomposed unbalance signals under 1 800 r/min

圖16 2 400 r/min工況下不平衡故障信號分解Fig.16 Decomposed unbalance signals under 2 400 r/min

圖17 3 000 r/min工況下不平衡故障信號分解Fig.17 Decomposed unbalance signals under 3 000 r/min

為讓故障診斷效果更準確，計算效率更高，根據轉子不平衡故障時1×處振幅劇烈增加的故障現象，對各轉速工況求取并繪制IMF頻譜圖，并篩選出最能表征不平衡故障特征的IMF，圖18～圖22分別為上文中選取的600 r/min、1 200 r/min、1 800 r/min、2 400 r/min、3 000 r/min工況下的轉子不平衡故障信號與其IMF頻譜圖。其中，X為原始信號頻譜圖，IMFi(i=1,2,…,n)為經VMD分解得到的IMF頻譜圖，圖中數據游標處即為1×劇烈增長的地方，從頻譜圖中篩選出的最能表征不平衡故障的IMF見表2，其他組數據也采用相同方法進行IMF篩選。

圖18 600 r/min工況下不平衡故障信號頻譜圖Fig.18 Spectrum of unbalance signals under 600 r/min

圖19 1 200 r/min工況下不平衡故障信號頻譜圖Fig.19 Spectrum of unbalance signals under 1 200 r/min

圖20 1 800 r/min工況下不平衡故障信號頻譜圖Fig.20 Spectrum of unbalance signals under 1 800 r/min

圖21 2 400 r/min工況下不平衡故障信號頻譜圖Fig.21 Spectrum of unbalance signals under 2 400 r/min

圖22 3 000 r/min工況下不平衡故障信號頻譜圖Fig.22 Spectrum of unbalance signals under 3 000 r/min

表2 各轉速頻譜圖中篩選出的IMFTab.2 Selected IMF from spectrums under each speed

在進行MPE量化篩選出的IMF時，需要設定三個參數，分別為：尺度因子s、嵌入維數m和時延λ；一般來說，嵌入維數m取3～7；因為，如m過小，則重構序列中包含的狀態過少，不能反映時間序列的動力學變化；如m過大，不僅消耗時間，而且無法反映時間序列的微小變化。時延λ對時間序列的影響可忽略不計，通常來說，尺度因子s的取值大于10即可。本文參考了文獻[19]中的參數取值，s取12，m取6，λ取1。分別對五種轉速，共計500組數據求取MPE，并計算每種轉速工況下各尺度因子排列熵的均值，結果如圖23所示，五種轉速的轉子振動信號在每一尺因子上的MPE均值區分度很明顯。這是因為，轉子不平衡故障處于不同轉速時，信號的隨機性會發生改變，從而導致其對應的多尺度排列熵發生了變化。MPE可以高效、準確地衡量多尺度情況下微小的信號變化。至此，MPE法量化IMF成功。

3.2 變轉速工況下轉子不平衡故障識別分類

將每種轉速工況下轉子不平衡故障數據前50組作為已知樣本，將已知故障樣本用MPE法構造成特征向量并輸入FCM，以求取每種轉速工況下轉子不平衡故障對應的標準聚類中心，如圖24所示。由圖24可知，五種轉速工況下轉子不平衡故障信號在每一尺度s上都具有明顯差異，可以進行識別與分類。

圖24 五種轉速工況下不平衡故障信號標準聚類中心Fig.24 Standard clustering center of unbalance signals under five speeds

將每種轉速不平衡位移數據的后50組作為待檢測數據，采用格貼近，最小最大貼近，海明貼近，歐幾里得貼近四種模糊算法，求取各轉速工況下待檢測樣本與標準聚類中心的貼近度，再從四種模糊貼近算法中選取出識別變轉速工況下轉子不平衡故障效果最好的算法，如圖25所示。由圖25中的貼近度分布可知，海明貼近度和歐幾里得貼近度的值最接近1，故上述兩種算法作用于變轉速工況轉子不平衡故障的效果最好。

圖25 四種模糊算法的貼近度分布Fig.25 Distribution of closeness degree of four fuzzy closeness algorithms

將各轉速工況下轉子不平衡故障的待檢測數據混合。采用海明與歐幾里得貼近算法求取混合數據與每種轉速標準聚類中心的貼近度。貼近度最大者表明為一類，當兩組測試樣本的貼近度都較大且相差小于0.01時，認為識別失敗[20]。采用此規則對五種不同轉速工況的轉子不平衡故障進行識別分類，其結果如圖26、表3所示。由圖表分析知，600 r/min下轉子不平衡識別正確率平均值為94.8%，1 200 r/min下轉子不平衡識別正確率均值為93.8%，1 800 r/min下轉子不平衡識別正確率均值為80.6%，2 400 r/min下轉子不平衡識別正確率均值為81.1%，3 000 r/min下轉子不平衡識別正確率均值為88.2%。總體而言，采用VMD_MPE-FCM方法進行變轉速工況下轉子不平衡故障的識別效果好。

圖26 混合識別正確率Fig.26 Accuracy of mixed recognition

表3 混合識別正確率值Tab.3 Values of mixed recognition

4 結 論

本文提出一種變轉速工況下轉子不平衡故障的診斷方法。首先，根據不平衡故障時，1×處振幅劇烈增加的現象，從VMD得到的IMF頻譜圖中篩選出最能表征不平衡故障的IMF，再采用MPE法量化篩選出的IMF，將作為特征向量的MPE輸入FCM，得到每種轉速工況下的標準聚類中心。最后，采用擇近原則和模糊貼近算法，計算出待識別數據與標準聚類中心的貼近度，從而實現變轉速工況下轉子不平衡故障的識別與分類。通過在DHRMT轉子試驗臺進行變轉速轉子的不平衡故障試驗分析，得到如下結論：

(1)采用VMD與轉子不平衡故障機理現象結合的方法對變轉速工況下轉子不平衡故障信號進行分析，可以提高故障診斷的正確率及算法的效率。

(2)采用MPE法量化篩選出的IMF，可以有效地檢測出變轉速工況下轉子振動信號的動態變化，并且多尺度可以更好地反映不平衡故障的本質特征。

(3)采用海明貼近和歐幾里得貼近模糊算法，可以更好地表征待檢測數據與聚類中心的貼近度，提高變轉速工況下轉子不平衡故障識別正確率。

(4)VMD_MPE-FCM方法可以有效地提取出轉子不平衡故障的特征，并可實現多轉速工況下轉子不平衡故障的識別與分類。