許偉

(上海電氣集團(tuán)股份有限公司 中央研究院,上海 200070)

1 概述

滾動軸承是旋轉(zhuǎn)機械中極為關(guān)鍵的承載部件,廣泛應(yīng)用于航天、冶金、交通、新能源等領(lǐng)域。如果不能及時發(fā)現(xiàn)并消除意外故障,可能會導(dǎo)致機械設(shè)備進(jìn)一步劣化,造成異常停機等嚴(yán)重后果,因此滾動軸承早期故障診斷具有重要研究意義。在滾動軸承失效進(jìn)程的早期,故障沖擊振動信號往往十分微弱,容易被背景噪聲和其他干擾所淹沒,失效特征難以準(zhǔn)確提取,需進(jìn)行早期故障特征增強方法的研究。

基于濾波的方法可有效提高故障特征的提取效果,其關(guān)鍵在于設(shè)計濾波器或通過解卷積優(yōu)化濾波器。解卷積方法具有靈活性和自適應(yīng)性,最小熵解卷積(Minimum Entropy Deconvolution,MED)是一種經(jīng)典的解卷積方法[1],其將峭度作為信號沖擊性的度量，構(gòu)建逆濾波器并提取故障脈沖。盡管MED在多種應(yīng)用中都有不錯的效果,但易受到背景噪聲和沖擊干擾的影響[2-7]。為進(jìn)一步改善MED的應(yīng)用效果并增強周期性故障特征提取能力,文獻(xiàn)[8]提出了一種基于最大相關(guān)峭度解卷積(Maximum Correlated Kurtosis Deconvolution,MCKD)方法,其引入了相關(guān)峭度(CK)作為判據(jù),通過考慮故障脈沖的周期性有效抑制了不同故障脈沖之間的相關(guān)性:文獻(xiàn)[9]在MCKD中加入故障特征周期的自適應(yīng)估計,提高了該算法的效率和魯棒性;文獻(xiàn)[10]將MCKD作為后處理方法,結(jié)合自適應(yīng)局部迭代濾波提取故障脈沖特征,并用粒子群算法對解卷積參數(shù)進(jìn)行優(yōu)化。盡管MED和MCKD有多種改進(jìn)方法,但本質(zhì)上都基于峭度最大化。由于峭度是高階統(tǒng)計量指標(biāo),利用目標(biāo)函數(shù)法只能通過迭代進(jìn)行求解,上述方法效率有限且求解結(jié)果可能僅是局部最優(yōu)解。

在實際應(yīng)用中,解卷積方法的性能主要取決于判據(jù)的選擇。文獻(xiàn)[11]將最小廣義Lp/Lq解卷積與譜峭度結(jié)合用于滾動軸承故障特征提取,文獻(xiàn)[12]提出了基于基尼指數(shù)的盲解卷積方法框架,但上述方法仍然被迭代問題所困擾。文獻(xiàn)[13]提出了一種稱為D范數(shù)的解卷積判據(jù),并提出了最優(yōu)最小熵解卷積(Optimal MED,OMED)方法,其與MED在幾何上等價,但卻無需迭代過程。文獻(xiàn)[14]對OMED方法進(jìn)行了改進(jìn),通過調(diào)整Toeplitz矩陣消除了由數(shù)據(jù)不連續(xù)帶來的虛假脈沖并將其命名為最優(yōu)最小熵反褶積調(diào)整(OMED Adjusted,OMEDA)。與MCKD類似,OMEDA中進(jìn)一步針對故障特征的周期性提出了多點最優(yōu)最小熵解卷積調(diào)整(Multipoint OMEDA,MOMEDA)方法,雖然相對于MCKD并不能有效增強提取效果,但效率有明顯提高。在實際應(yīng)用中可通過優(yōu)化濾波器長度、實際故障周期等參數(shù)以獲得更好的提取效果。

上述方法多應(yīng)用于工況穩(wěn)定的情況,而大多數(shù)機械裝備都很難保證平穩(wěn)的服役工況,也導(dǎo)致基于周期性的解卷積方法在很多情況下特征提取能力有限。為解決這個問題,文獻(xiàn)[15]基于故障特征的二階循環(huán)平穩(wěn)性最大化準(zhǔn)則開發(fā)了最大二階循環(huán)平穩(wěn)盲解卷積(Maximum Second-order Cyclostationarity Blind Deconvolution,CYCBD)及其改進(jìn)方法CYCBDang,此類方法可以使用測得的瞬時轉(zhuǎn)速在轉(zhuǎn)速非恒定時精確提取故障脈沖,與其他解卷積方法相比有很大的優(yōu)勢。文獻(xiàn)[16]提出了最大平均峭度解卷積(Maximum Average Kurtosis Deconvolution,MAKD),其根據(jù)循環(huán)頻率對濾波后的信號進(jìn)行分割并將信號片段的平均峭度作為判據(jù),不僅適用于變轉(zhuǎn)速工況,而且能很好地抑制隨機干擾的影響。與經(jīng)驗?zāi)B(tài)分解、階次跟蹤等[17-18]經(jīng)典變轉(zhuǎn)速故障診斷方法相比,上述解卷積方法可以在不進(jìn)行重采樣或不受模態(tài)耦合影響的情況下突出故障分量,但仍需通過迭代求解。此外,階次跟蹤、廣義解調(diào)等非平穩(wěn)信號處理方法在恢復(fù)信號平穩(wěn)性的同時會導(dǎo)致信號中時不變共振特征的缺失或損壞,給后續(xù)的特征提取帶來很大困擾。因此,如何在變轉(zhuǎn)速工況下快速有效地識別滾動軸承早期故障仍然是一項具有挑戰(zhàn)性的任務(wù)。

為滿足上述要求,本文提出了一種改進(jìn)多點最優(yōu)最小熵解卷積調(diào)整(Improved Multipoint Optimal Minimum Entropy Deconvolution Adjusted,IMOMEDA)方法,在不進(jìn)行重采樣,不破壞時不變共振特征的前提下,首先,根據(jù)脈沖在角度域的等角度分布特性計算故障脈沖的時域間隔并構(gòu)成目標(biāo)向量,使循環(huán)頻率隨轉(zhuǎn)速動態(tài)變化,實現(xiàn)故障特征的表征;然后,在將其擴(kuò)展到變轉(zhuǎn)速工況的同時,通過分段計算多D范數(shù)保證原有的高計算效率;最后,引入包絡(luò)諧波強度進(jìn)一步自動定位濾波信號的循環(huán)調(diào)制階次,為故障診斷提供依據(jù),以實現(xiàn)滾動軸承早期故障特征提取。

2 MOMEDA解卷積方法原理

在恒定轉(zhuǎn)速條件下,局部故障引起的脈沖力周期性地激勵在旋轉(zhuǎn)部件上使機器產(chǎn)生共振,這種共振可以被振動傳感器獲取,但易被各種干擾所埋沒。為增強故障信息,MOMEDA采用新的解卷積方式設(shè)計濾波器,使濾波信號中的脈沖分量最大化。

假設(shè)給定的振動信號為x=[x1,x2,…,xN]T,MOMEDA的目標(biāo)在于找到最優(yōu)逆濾波器,使濾波信號y的多D范數(shù)最大化。多D范數(shù)的表達(dá)式為

(1)

基于最大化目標(biāo)函數(shù),假定逆濾波器f可表示為長度L的FIR濾波器,則MOMEDA可變?yōu)橐韵聝?yōu)化問題,即

(2)

式中:t為確定脈沖位置和權(quán)重的目標(biāo)向量。通過設(shè)計合適的目標(biāo)向量t,可以對任意位置和任意權(quán)重的目標(biāo)脈沖進(jìn)行解卷積。通過求解多D范數(shù)的最大值獲得最優(yōu)逆濾波器系數(shù),進(jìn)而獲得最終的濾波信號。

對逆濾波器系數(shù)求導(dǎo)可得

(3)

其中

(4)

Mk=[xk+L-1,xk+L-2,…,xk],

(5)

(6)

則(3)式可寫為

(7)

t1M1+t2M2+…+tN-LMN-L=X0t。

(8)

通過導(dǎo)數(shù)等于0求解極值可得

(9)

即

(10)

(11)

由于任意倍數(shù)均為(11)式的解,所以MOMEDA的逆濾波器和輸出解可表示為

(12)

(13)

3 IMOMEDA方法

使用MOMEDA提取周期性脈沖時,目標(biāo)向量t中的非零元素是周期性的,從而可以捕獲故障脈沖的位置;但當(dāng)機械設(shè)備轉(zhuǎn)速具有明顯時變時,故障部件的時域周期性就不再成立,MOMEDA的提取效果不佳。盡管可以在逐步遍歷整個信號時每次只對幾個連續(xù)脈沖進(jìn)行解卷積以減小轉(zhuǎn)速變化的影響,但該過程并沒有充分利用信號的全局信息來提高結(jié)果的有效性。為充分利用信號的全局信息,本文利用先驗轉(zhuǎn)速信息提出用于處理變轉(zhuǎn)速工況下解卷積問題的IMOMEDA方法。

IMOMEDA的關(guān)鍵在于重新構(gòu)建目標(biāo)向量t,使其可以隨轉(zhuǎn)速變化自適應(yīng)地調(diào)整非零元素的位置。首先,定義v=[v1,v2,…,vN]T為與輸入振動信號x相對應(yīng)的轉(zhuǎn)速信號,單位為r/min,則累計轉(zhuǎn)角可近似計算為

(14)

根據(jù)脈沖的故障特征階次OFCO,角度域脈沖間隔為

Δθ=360/OFCO,

(15)

則角度域中的目標(biāo)向量可定義為

(16)

式中:fs為轉(zhuǎn)頻;K為脈沖總數(shù);θ0為第1個脈沖所在角位移;δ(·)為階躍函數(shù)。

角度域目標(biāo)向量tn中的非零元素具有周期性,可以通過調(diào)整θ0構(gòu)造一個目標(biāo)向量來計算最優(yōu)逆濾波器;但由于離散采樣的原因,在tn非零元素點可能沒有對應(yīng)的采樣點,因此采用最近的采樣點作為近似,將(16)式松弛為

(17)

通過(17)式構(gòu)造的角度域目標(biāo)向量tn使多D范數(shù)具有較好的突出脈沖的能力。但振動信號是固定時間間隔采樣,在角度域上并非等間距,這意味著相鄰脈沖之間的采樣點數(shù)不是常數(shù),因此時域上目標(biāo)函數(shù)Nt可定義為

(18)

在實際應(yīng)用中,不僅需要提取故障脈沖,更需要對提取的故障脈沖進(jìn)行度量以作為預(yù)測性維護(hù)或診斷的依據(jù)。在MOMEDA中,將多點峭度作為故障脈沖的度量指標(biāo),其定義為

(19)

基于多點峭度的定義,采用故障周期、整數(shù)倍周期和分?jǐn)?shù)倍周期構(gòu)造目標(biāo)向量同樣可以成功提取故障脈沖,多點峭度最大時對應(yīng)的采樣點數(shù)可能為故障周期或其整數(shù)倍或分?jǐn)?shù)倍,這使得故障脈沖的確切特征階次難以直接識別。為解決這個問題,本文基于階次平方包絡(luò)譜定義包絡(luò)諧波強度并形成包絡(luò)諧波強度譜圖,以突出信號在角度域的周期性調(diào)制,定位主要分量的故障特征階次。包絡(luò)諧波強度的表達(dá)式為

(20)

式中:AOSEA(O)為階次平方包絡(luò)譜在第O階的幅值;REHIn(O)為前n個幅值的幾何平均值。

在包絡(luò)諧波強度譜圖中可以識別到的潛在故障特征階次為

Od=argmaxREHIn(O),

(21)

當(dāng)Od接近某個故障特征階次時,可認(rèn)為發(fā)生了相應(yīng)故障,從而實現(xiàn)故障特征階次的自動定位。IMOMEDA方法的具體流程如圖1所示。

圖1 IMOMEDA故障診斷流程

4 仿真結(jié)果分析

為評價IMOMEDA方法的性能,通過仿真信號進(jìn)行研究。仿真信號x(t)由故障沖擊信號x1(t)、離散諧波信號x2(t)以及高斯白噪聲n(t)組成,即

x(t)=x1(t)+x2(t)+n(t),

(22)

(23)

(24)

(25)

仿真信號中:x1(t)表示滾動軸承故障引起的周期脈沖信號,采樣頻率為20 kHz,采樣時間為1 s;故障脈沖參數(shù)A1,A2分別為1,0.25;s(t)為單個脈沖;a為諧振阻尼系數(shù),取500;fn為共振頻率,取2 000;Ti為故障脈沖位置;f0(t)為瞬時轉(zhuǎn)速。x2(t)表示來自軸等其他部件的離散諧波信號,該部分假定出現(xiàn)了2個諧波,參數(shù)分別為C1=0.2,R1=1以及C2=0.1,R2=20;fk(t)為對應(yīng)的瞬時頻率;高斯白噪聲n(t)的標(biāo)準(zhǔn)差為0.4。

為充分說明IMOMEDA方法在時變轉(zhuǎn)速下提取故障特征的能力,采用4種不同的速度曲線(拋物線、線性曲線、指數(shù)曲線和余弦曲線)產(chǎn)生信號并進(jìn)行分析,轉(zhuǎn)速變化曲線如圖2所示。將CYCBDang作為對比方法,同時在Intel(R)Core i7-7700K CPU和32GB內(nèi)存的計算機上運行,數(shù)據(jù)處理軟件為MATLAB R2020a。

圖2 仿真工況轉(zhuǎn)速變化曲線

轉(zhuǎn)速曲線線性變化時模擬信號所含分量的波形如圖3所示,循環(huán)脈沖相對于轉(zhuǎn)速的階數(shù)設(shè)為5?；旌虾笮盘栔袩o法觀察到循環(huán)脈沖,幾乎完全被干擾所淹沒。

(a) 故障沖擊信號

IMOMEDA與CYCBDang在這些信號上的處理時間如圖4所示,可以明顯看出IMOMEDA的處理時間遠(yuǎn)小于CYCBDang。CYCBDang需要迭代計算直到得到最優(yōu)解或局部最優(yōu)解,而IMOMEDA則無需迭代過程,因此可以在很短的時間內(nèi)完成。

圖4 IMOMEDA和CYCBDang的計算時間

4種速度條件下,IMOMEDA和CYCBDang的提取結(jié)果如圖5所示,2種方法都得到了相似的結(jié)果,說明IMOMEDA在提取效率提高的同時不影響其提取效果的有效性。

圖5 不同方法獲得的濾波信號

故障類型的識別是機械設(shè)備預(yù)測性維護(hù)的另一個關(guān)鍵問題。對于MOMEDA,CYCBDang,CYCBD和MCKD等解卷積方法,即使輸入的周期/階次是準(zhǔn)確故障周期/階次的倍數(shù)或分?jǐn)?shù)時,仍然可以提取故障特征,但此時輸入的周期/階次也可能對應(yīng)其他部件的故障,會產(chǎn)生誤導(dǎo)。因此,準(zhǔn)確識別故障特征的循環(huán)頻率或階次至關(guān)重要。

MOMEDA采用多點峭度進(jìn)行故障脈沖的度量,但并不適用于故障脈沖特征階次的自動定位。不同故障特征階次信號的多點峭度譜圖如圖6所示:對于故障特征階次為10,20的信號,多點峭度難以直接定位故障特征階次的準(zhǔn)確值,在10階和20階處均存在明顯峰值,甚至在分?jǐn)?shù)倍階次的5階處峰值更大,易造成故障特征階次的誤判。

(a) OFCO=10

IMOMEDA處理不同故障特征階次信號時的包絡(luò)諧波強度譜圖如圖7所示:當(dāng)故障特征階次為10,20階時,EHI譜圖中最大值對應(yīng)的階次同樣為10,20階,并且對應(yīng)幅值相對于其他階次更大,不易造成誤判。因此,可以直接在包絡(luò)諧波強度譜圖中進(jìn)行故障特征階次的定位,且方法簡單。

(a) OFCO=10

5 試驗驗證

在某軌道交通企業(yè)的整車滾動綜合試驗臺進(jìn)行振動測試,試驗對象為某型在役走行部軸箱軸承,機車試驗臺架單次可進(jìn)行單節(jié)車輛的滾動試驗,同時搭載8套CRI-2692試驗軸承,試驗軸承的結(jié)構(gòu)參數(shù)及各特征階次見表1。

表1 軸承設(shè)計參數(shù)及各特征階次

5.1 軸承外圈故障

軸承外圈故障試驗信號的采樣頻率為20 kHz,共采集4 s,實際工況下的轉(zhuǎn)速變化如圖8所示,從130 r/min勻減速到73.2 r/min,該過程中采集到的軸承外圈故障信號及其階次包絡(luò)譜如圖9所示,從時域波形中可以看到部分瞬時脈沖,但在階次包絡(luò)譜中難以觀察到顯著的沖擊成分。

將IMOMEDA和CYCBDang的濾波器長度設(shè)定為200,處理所得濾波信號如圖10所示,IMOMEDA和CYCBDang明顯增強了故障脈沖,在階次包絡(luò)譜中可以發(fā)現(xiàn)顯著的沖擊成分。另外,IMOMEDA的處理時間為1.91 s,遠(yuǎn)遠(yuǎn)小于CYCBDang的194.39 s,在獲得相似結(jié)果的同時大大提高了處理效率。

(a) IMOMEDA

IMOMEDA處理時軸承外圈故障信號的包絡(luò)諧波強度譜圖如圖11所示,從圖中可以明顯看到EHI譜圖中最大值對應(yīng)的階次為外圈故障特征階次,且峰值遠(yuǎn)高于其他階次,可以準(zhǔn)確定位故障特征階次。

圖11 軸承外圈故障濾波信號的包絡(luò)諧波強度

5.2 軸承滾動體故障

軸承滾動體故障試驗信號的采樣頻率為20 kHz,實際工況下的轉(zhuǎn)速變化如圖12所示,從58 r/min勻加速到84 r/min,該過程中采集到的滾動體故障信號及其包絡(luò)譜如圖13所示, 從時域波形中可以看到部分瞬時脈沖,在階次包絡(luò)譜中也可以發(fā)現(xiàn)部分峰值,但諧波成分仍然不顯著。

圖12 滾動體故障工況的轉(zhuǎn)速變化曲線

圖13 滾動體故障信號及其包絡(luò)譜

將IMOMEDA和CYCBDang的濾波器長度設(shè)定為200,處理所得濾波信號如圖14所示,IMOMEDA明顯增強了故障脈沖,在階次包絡(luò)譜中可以發(fā)現(xiàn)顯著的峰值及諧波成分,效果優(yōu)于CYCBDang。CYCBDang的處理時間為192.61 s,而IMOMEDA的處理時間僅為2.64 s,同樣大大縮短了處理時間。

(a) IMOMEDA

IMOMEDA處理時軸承滾動體故障信號的包絡(luò)諧波強度譜如圖15所示，可以明顯看到EHI譜圖中最大值對應(yīng)的階次為滾動體故障特征階次且峰值高于其他階次，可以準(zhǔn)確定位故障特征階次。

圖15 濾波信號的包絡(luò)諧波強度

6 結(jié)論

1)IMOMEDA方法將MOMEDA中以固定間隔的目標(biāo)向量替換為隨轉(zhuǎn)速變化的目標(biāo)向量,使其獲得時變轉(zhuǎn)速條件下的振動信號分析能力。

2)相較于CYCBDang等解卷積方法,IMOMEDA無需進(jìn)行迭代求解,計算效率更高。

3)基于提出的包絡(luò)諧波強度可以準(zhǔn)確定位故障特征階次,解決了MOMEDA基于多點峭度而無法識別故障脈沖確切特征階次的問題。