許偉
(上海電氣集團股份有限公司 中央研究院,上海 200070)
滾動軸承是旋轉機械中極為關鍵的承載部件,廣泛應用于航天、冶金、交通、新能源等領域。如果不能及時發現并消除意外故障,可能會導致機械設備進一步劣化,造成異常停機等嚴重后果,因此滾動軸承早期故障診斷具有重要研究意義。在滾動軸承失效進程的早期,故障沖擊振動信號往往十分微弱,容易被背景噪聲和其他干擾所淹沒,失效特征難以準確提取,需進行早期故障特征增強方法的研究。
基于濾波的方法可有效提高故障特征的提取效果,其關鍵在于設計濾波器或通過解卷積優化濾波器。解卷積方法具有靈活性和自適應性,最小熵解卷積(Minimum Entropy Deconvolution,MED)是一種經典的解卷積方法[1],其將峭度作為信號沖擊性的度量,構建逆濾波器并提取故障脈沖。盡管MED在多種應用中都有不錯的效果,但易受到背景噪聲和沖擊干擾的影響[2-7]。為進一步改善MED的應用效果并增強周期性故障特征提取能力,文獻[8]提出了一種基于最大相關峭度解卷積(Maximum Correlated Kurtosis Deconvolution,MCKD)方法,其引入了相關峭度(CK)作為判據,通過考慮故障脈沖的周期性有效抑制了不同故障脈沖之間的相關性:文獻[9]在MCKD中加入故障特征周期的自適應估計,提高了該算法的效率和魯棒性;文獻[10]將MCKD作為后處理方法,結合自適應局部迭代濾波提取故障脈沖特征,并用粒子群算法對解卷積參數進行優化。盡管MED和MCKD有多種改進方法,但本質上都基于峭度最大化。由于峭度是高階統計量指標,利用目標函數法只能通過迭代進行求解,上述方法效率有限且求解結果可能僅是局部最優解。
在實際應用中,解卷積方法的性能主要取決于判據的選擇。文獻[11]將最小廣義Lp/Lq解卷積與譜峭度結合用于滾動軸承故障特征提取,文獻[12]提出了基于基尼指數的盲解卷積方法框架,但上述方法仍然被迭代問題所困擾。文獻[13]提出了一種稱為D范數的解卷積判據,并提出了最優最小熵解卷積(Optimal MED,OMED)方法,其與MED在幾何上等價,但卻無需迭代過程。文獻[14]對OMED方法進行了改進,通過調整Toeplitz矩陣消除了由數據不連續帶來的虛假脈沖并將其命名為最優最小熵反褶積調整(OMED Adjusted,OMEDA)。與MCKD類似,OMEDA中進一步針對故障特征的周期性提出了多點最優最小熵解卷積調整(Multipoint OMEDA,MOMEDA)方法,雖然相對于MCKD并不能有效增強提取效果,但效率有明顯提高。在實際應用中可通過優化濾波器長度、實際故障周期等參數以獲得更好的提取效果。
上述方法多應用于工況穩定的情況,而大多數機械裝備都很難保證平穩的服役工況,也導致基于周期性的解卷積方法在很多情況下特征提取能力有限。為解決這個問題,文獻[15]基于故障特征的二階循環平穩性最大化準則開發了最大二階循環平穩盲解卷積(Maximum Second-order Cyclostationarity Blind Deconvolution,CYCBD)及其改進方法CYCBDang,此類方法可以使用測得的瞬時轉速在轉速非恒定時精確提取故障脈沖,與其他解卷積方法相比有很大的優勢。文獻[16]提出了最大平均峭度解卷積(Maximum Average Kurtosis Deconvolution,MAKD),其根據循環頻率對濾波后的信號進行分割并將信號片段的平均峭度作為判據,不僅適用于變轉速工況,而且能很好地抑制隨機干擾的影響。與經驗模態分解、階次跟蹤等[17-18]經典變轉速故障診斷方法相比,上述解卷積方法可以在不進行重采樣或不受模態耦合影響的情況下突出故障分量,但仍需通過迭代求解。此外,階次跟蹤、廣義解調等非平穩信號處理方法在恢復信號平穩性的同時會導致信號中時不變共振特征的缺失或損壞,給后續的特征提取帶來很大困擾。因此,如何在變轉速工況下快速有效地識別滾動軸承早期故障仍然是一項具有挑戰性的任務。
為滿足上述要求,本文提出了一種改進多點最優最小熵解卷積調整(Improved Multipoint Optimal Minimum Entropy Deconvolution Adjusted,IMOMEDA)方法,在不進行重采樣,不破壞時不變共振特征的前提下,首先,根據脈沖在角度域的等角度分布特性計算故障脈沖的時域間隔并構成目標向量,使循環頻率隨轉速動態變化,實現故障特征的表征;然后,在將其擴展到變轉速工況的同時,通過分段計算多D范數保證原有的高計算效率;最后,引入包絡諧波強度進一步自動定位濾波信號的循環調制階次,為故障診斷提供依據,以實現滾動軸承早期故障特征提取。
在恒定轉速條件下,局部故障引起的脈沖力周期性地激勵在旋轉部件上使機器產生共振,這種共振可以被振動傳感器獲取,但易被各種干擾所埋沒。為增強故障信息,MOMEDA采用新的解卷積方式設計濾波器,使濾波信號中的脈沖分量最大化。
假設給定的振動信號為x=[x1,x2,…,xN]T,MOMEDA的目標在于找到最優逆濾波器,使濾波信號y的多D范數最大化。多D范數的表達式為
(1)
基于最大化目標函數,假定逆濾波器f可表示為長度L的FIR濾波器,則MOMEDA可變為以下優化問題,即
(2)
式中:t為確定脈沖位置和權重的目標向量。通過設計合適的目標向量t,可以對任意位置和任意權重的目標脈沖進行解卷積。通過求解多D范數的最大值獲得最優逆濾波器系數,進而獲得最終的濾波信號。
對逆濾波器系數求導可得
(3)
其中
(4)
Mk=[xk+L-1,xk+L-2,…,xk],
(5)
(6)
則(3)式可寫為
(7)
t1M1+t2M2+…+tN-LMN-L=X0t。
(8)
通過導數等于0求解極值可得
(9)
即
(10)

(11)
由于任意倍數均為(11)式的解,所以MOMEDA的逆濾波器和輸出解可表示為
(12)
(13)
使用MOMEDA提取周期性脈沖時,目標向量t中的非零元素是周期性的,從而可以捕獲故障脈沖的位置;但當機械設備轉速具有明顯時變時,故障部件的時域周期性就不再成立,MOMEDA的提取效果不佳。盡管可以在逐步遍歷整個信號時每次只對幾個連續脈沖進行解卷積以減小轉速變化的影響,但該過程并沒有充分利用信號的全局信息來提高結果的有效性。為充分利用信號的全局信息,本文利用先驗轉速信息提出用于處理變轉速工況下解卷積問題的IMOMEDA方法。
IMOMEDA的關鍵在于重新構建目標向量t,使其可以隨轉速變化自適應地調整非零元素的位置。首先,定義v=[v1,v2,…,vN]T為與輸入振動信號x相對應的轉速信號,單位為r/min,則累計轉角可近似計算為
(14)
根據脈沖的故障特征階次OFCO,角度域脈沖間隔為
Δθ=360/OFCO,
(15)
則角度域中的目標向量可定義為
(16)
式中:fs為轉頻;K為脈沖總數;θ0為第1個脈沖所在角位移;δ(·)為階躍函數。
角度域目標向量tn中的非零元素具有周期性,可以通過調整θ0構造一個目標向量來計算最優逆濾波器;但由于離散采樣的原因,在tn非零元素點可能沒有對應的采樣點,因此采用最近的采樣點作為近似,將(16)式松弛為
(17)
通過(17)式構造的角度域目標向量tn使多D范數具有較好的突出脈沖的能力。但振動信號是固定時間間隔采樣,在角度域上并非等間距,這意味著相鄰脈沖之間的采樣點數不是常數,因此時域上目標函數Nt可定義為
(18)

在實際應用中,不僅需要提取故障脈沖,更需要對提取的故障脈沖進行度量以作為預測性維護或診斷的依據。在MOMEDA中,將多點峭度作為故障脈沖的度量指標,其定義為
(19)
基于多點峭度的定義,采用故障周期、整數倍周期和分數倍周期構造目標向量同樣可以成功提取故障脈沖,多點峭度最大時對應的采樣點數可能為故障周期或其整數倍或分數倍,這使得故障脈沖的確切特征階次難以直接識別。為解決這個問題,本文基于階次平方包絡譜定義包絡諧波強度并形成包絡諧波強度譜圖,以突出信號在角度域的周期性調制,定位主要分量的故障特征階次。包絡諧波強度的表達式為
(20)
式中:AOSEA(O)為階次平方包絡譜在第O階的幅值;REHIn(O)為前n個幅值的幾何平均值。
在包絡諧波強度譜圖中可以識別到的潛在故障特征階次為
Od=argmaxREHIn(O),
(21)
當Od接近某個故障特征階次時,可認為發生了相應故障,從而實現故障特征階次的自動定位。IMOMEDA方法的具體流程如圖1所示。

圖1 IMOMEDA故障診斷流程
為評價IMOMEDA方法的性能,通過仿真信號進行研究。仿真信號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;故障脈沖參數A1,A2分別為1,0.25;s(t)為單個脈沖;a為諧振阻尼系數,取500;fn為共振頻率,取2 000;Ti為故障脈沖位置;f0(t)為瞬時轉速。x2(t)表示來自軸等其他部件的離散諧波信號,該部分假定出現了2個諧波,參數分別為C1=0.2,R1=1以及C2=0.1,R2=20;fk(t)為對應的瞬時頻率;高斯白噪聲n(t)的標準差為0.4。
為充分說明IMOMEDA方法在時變轉速下提取故障特征的能力,采用4種不同的速度曲線(拋物線、線性曲線、指數曲線和余弦曲線)產生信號并進行分析,轉速變化曲線如圖2所示。將CYCBDang作為對比方法,同時在Intel(R)Core i7-7700K CPU和32GB內存的計算機上運行,數據處理軟件為MATLAB R2020a。

圖2 仿真工況轉速變化曲線
轉速曲線線性變化時模擬信號所含分量的波形如圖3所示,循環脈沖相對于轉速的階數設為5。混合后信號中無法觀察到循環脈沖,幾乎完全被干擾所淹沒。

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

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

圖5 不同方法獲得的濾波信號
故障類型的識別是機械設備預測性維護的另一個關鍵問題。對于MOMEDA,CYCBDang,CYCBD和MCKD等解卷積方法,即使輸入的周期/階次是準確故障周期/階次的倍數或分數時,仍然可以提取故障特征,但此時輸入的周期/階次也可能對應其他部件的故障,會產生誤導。因此,準確識別故障特征的循環頻率或階次至關重要。
MOMEDA采用多點峭度進行故障脈沖的度量,但并不適用于故障脈沖特征階次的自動定位。不同故障特征階次信號的多點峭度譜圖如圖6所示:對于故障特征階次為10,20的信號,多點峭度難以直接定位故障特征階次的準確值,在10階和20階處均存在明顯峰值,甚至在分數倍階次的5階處峰值更大,易造成故障特征階次的誤判。

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

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

表1 軸承設計參數及各特征階次
軸承外圈故障試驗信號的采樣頻率為20 kHz,共采集4 s,實際工況下的轉速變化如圖8所示,從130 r/min勻減速到73.2 r/min,該過程中采集到的軸承外圈故障信號及其階次包絡譜如圖9所示,從時域波形中可以看到部分瞬時脈沖,但在階次包絡譜中難以觀察到顯著的沖擊成分。
將IMOMEDA和CYCBDang的濾波器長度設定為200,處理所得濾波信號如圖10所示,IMOMEDA和CYCBDang明顯增強了故障脈沖,在階次包絡譜中可以發現顯著的沖擊成分。另外,IMOMEDA的處理時間為1.91 s,遠遠小于CYCBDang的194.39 s,在獲得相似結果的同時大大提高了處理效率。

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

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

圖12 滾動體故障工況的轉速變化曲線

圖13 滾動體故障信號及其包絡譜
將IMOMEDA和CYCBDang的濾波器長度設定為200,處理所得濾波信號如圖14所示,IMOMEDA明顯增強了故障脈沖,在階次包絡譜中可以發現顯著的峰值及諧波成分,效果優于CYCBDang。CYCBDang的處理時間為192.61 s,而IMOMEDA的處理時間僅為2.64 s,同樣大大縮短了處理時間。

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

圖15 濾波信號的包絡諧波強度
1)IMOMEDA方法將MOMEDA中以固定間隔的目標向量替換為隨轉速變化的目標向量,使其獲得時變轉速條件下的振動信號分析能力。
2)相較于CYCBDang等解卷積方法,IMOMEDA無需進行迭代求解,計算效率更高。
3)基于提出的包絡諧波強度可以準確定位故障特征階次,解決了MOMEDA基于多點峭度而無法識別故障脈沖確切特征階次的問題。