基于Vold-Kalman廣義解調的變轉速軸承和齒輪復合故障診斷

2019-04-03 01:17:42趙德尊李建勇程衛東溫偉剛

振動與沖擊 2019年6期

趙德尊, 李建勇, 程衛東, 溫偉剛

(1.清華大學機械工程系，北京 100084； 2.北京交通大學機械與電子控制工程學院,北京 100044; 3.載運工具先進制造與測控技術教育部重點實驗室,北京 100044)

軸承和齒輪是機械結構中用于支撐和傳遞力矩的重要部件[1]。由于復雜的工作環境和運行狀態，軸承和齒輪也是最易受損的零部件。實際工況中，軸承和齒輪往往同時發生故障。受轉速波動的影響，多個故障特征相互影響、彼此干擾的同時其自身重復頻率也表現出時變非平穩性，給故障診斷帶來困難。因此，探究變轉速條件下軸承和齒輪復合故障特征的提取具有重要意義。

對于故障軸承，隨著故障點與其表面間的接觸，測量的振動信號中將產生一系列短時、快速衰減的沖擊脈沖。沖擊脈沖的重復頻率即故障特征頻率是確定軸承故障點位置的重要指標，另外，由軸承故障引起的高頻共振將調制故障特征頻率[2]。對于故障齒輪，有文獻指出，其振動信號中除了含有具有幅值優勢的嚙合頻率外還包含由齒輪故障激起的共振帶，該共振帶同樣會調制齒輪轉頻[3]。轉頻可作為有效指標診斷齒輪故障[4]。據此，計算最優的濾波參數，區分提取軸承和齒輪分別激起的共振帶，是其復合故障特征提取的關鍵步驟，也引起了諸多學者的關注。

基于譜峭度對振動信號中沖擊脈沖的敏感性，Antoni[5]提出了快速 Kurtogram 算法，該算法可以有效確定最優帶通濾波參數，通過帶通濾波提取由機械部件局部故障引起的高頻共振。為進一步提高快速Kurtogram算法抗干擾能力，Wang等[6]提出了基于二進制小波包分解的快速 Kurtogram增強算法，并成功的應用于滾動軸承故障診斷。Combet等[7]提出了利用譜峭度值確定最優濾波參數的齒輪早期故障特征提取方法。基于譜峭度的最優濾波參數計算方法得到諸多學者認可和廣泛應用的同時，其對于復合共振帶提取的無能無力也引起了學者的重視。Wang等[8]改進了譜峭度濾波算法，即通過設定嚙合指標確定只含有齒輪共振帶的嚙合共振譜峭度圖，繼而與包含軸承和齒輪共振帶信息的Kurtogram圖對比分析完成軸承和齒輪復合共振帶的區分。然而，該方法只適用于恒轉速條件，對于變轉速復合故障特征的提取無能為力。

變轉速模式普遍存在于旋轉機械的實際運行中[9]。變工況下，表征軸承故障的故障特征頻率和表征齒輪故障的轉頻將表現出時變非平穩的特點。現有的以FFT(Fast Fourier Transform)為核心的故障診斷方法將因頻譜模糊現象而失去作用。針對傳統的計算階比分析方法在計算效率和精度方面的不足[10]，探索無需角域重采樣的復合故障特征提取方法對提高故障診斷技術的有效性和實用性具有重要意義。

Vold-Kalman階比跟蹤是Vold等[11]基于Kalman濾波器提出的可以用于估計振動信號中某一諧波成分的算法。與傳統的帶通濾波或經驗模式分解算法相比，該算法可在時域中直接提取某一特定的近似、交叉的諧波分量，避免了由時域至頻域變換帶來的相位偏差和能量泄露[11-12]，因此受到了諸多學者關注。其基本原理從第一代發展到第二代，具備了成熟的理論體系，并廣泛的應用于商業軟件中。然而現有的Vold-Kalman濾波算法均是與階比分析相結合處理變轉速旋轉機械振動信號[13-14]。廣義解調算法是Olhede等[15]提出的一種可以將時頻分布是傾斜、非線性的非平穩信號轉換成時頻分布是線性的且平行于時間軸的平穩信號的分析方法。基于廣義解調算法對于調幅-調頻信號的處理優勢，很多學者將其應用于旋轉機械設備的故障診斷[16]。

綜上所述，針對傳統的以譜峭度為核心的帶通濾波和計算階比分析算法的不足，結合Vold-Kalman濾波和廣義解調變換的優勢，提出了滾動軸承和齒輪時變非平穩復合故障特征提取方法。該方法舍棄了傳統的高頻共振帶提取和角域重采樣的故障診斷策略，直接利用Vold-Kalman濾波提取時變非平穩的故障特征成分，其次通過廣義解調實現提取的時變非平穩故障特征成分的平穩化重置，最后基于頻譜分析完成重置成分的定量表達。仿真和實驗分析結果證實了本文算法提取、重置和定量表達故障特征成分的有效性；同時避免了傳統帶通濾波算法帶來的相位偏差和能量泄露，頻譜中故障特征成分具有更明顯的幅值優勢。

1 理論基礎

1.1 Vold-Kalman濾波

旋轉機械中產生的調制信號可表示為幅值和載波信號的乘積

(1)

式中：Ak為第k個諧波成分的幅值；φk(t)為載波成分，其表達式如式(2)所示

(2)

由于機械系統的固有慣量，調制信號中的幅值包絡有別于載波信號以較低的頻率緩慢變化，因此可以用低階多項式表示。對于離散信號其狀態方程可表示為

(3)

設定s=2，式(3)可改寫為

Ak(n-1)-2Ak(n)+Ak(n+1)=εk(n)

(4)

將式(4)展開為矩陣

(5)

因此，不同階次的多項式具有相同的矩陣形式

MA=ε

(6)

式中：M為N×N矩陣。

實測信號包含各諧波成分和噪聲或誤差成分，其觀測方程可表示為

(7)

式中：y(n)為實測振動信號；δ(n)為噪聲或誤差項。

為提取特定的諧波分量xk(t)，系統的觀測方程可表示為

(8)

其矩陣形式為

Y-CA=δ

(9)

根據系統狀態方程式(6)和觀測式(9)以及預先設定的目標分量的載波矩陣C，應用最小二乘法，使得非齊次項ε(n)和噪聲或誤差δ(n)的平方和最小[17]，即

(10)

式中：C*為C的共軛轉置；A*為A的共軛轉置；r為加權因子，具體取值可參考文獻[18]。

進一步求得目標分量的幅值包絡矩陣A

A=(r2MTM+E)-1C*Y

(11)

通過載波矩陣C和幅值包絡A即可重構目標分量

X=AC

(12)

1.2 廣義解調變換

廣義解調變換的本質為廣義傅里葉變換。對于任意的解析信號x(t)，其廣義傅里葉變換的定義為

(13)

式中：s0(t)為隨時間t變化的實值函數，即相位函數；廣義傅里葉變換實際上是對x(t)e-2πjs0(t)做標準的傅里葉變換。同樣對XG(f)進行逆傅里葉變換可得

(14)

根據式(14)，如果XG(f)=δ(f-f0)，x(t)可改寫為

x(t)=e-2πj[f0t+s0(t)]

(15)

通過式(15)可知，解析信號中某一特定時變非平穩成分的相位函數s0(t)只要滿足式(16)，其即可被轉換成線性的、平行時間軸的平穩成分。

(16)

式中：f(t)為特定時變非平穩成分的頻率方程；f0為頻率方程對應的頻率點；s0(t)為頻率方程對應的相位函數；ds0(t)/dt為對s0(t)求導。

基于上述廣義解調變換的理論介紹，給出對原始信號x(t)中感興趣的某一時變非平穩成分的解調過程：①對原始信號x(t)進行希爾伯特變換，獲得解析信號y(t)=x(t)+jH[x(t)]，其中H[x(t)]是x(t)的希爾伯特變換結果；②根據感興趣的時變非平穩成分的頻率方程和式(16)計算其對應的相位函數s0(t)，基于相位函數得到變換結果d(t)=y(t)e-2πjs0(t)；③對變換結果d(t)進行希爾伯特變換得到新的解析信號z(t)=d(t)+jH[d(t)]，即原始信號x(t)的廣義解調信號。

2 本文算法

故障特征頻率和轉頻分別是判定軸承和齒輪健康狀況的關鍵性指標。本文利用Vold-Kalman濾波算法和廣義解調變換對時變非平穩的故障特征成分進行提取和平穩化重置，進而通過頻率譜對重置的故障特征成分進行定量表達，最終基于理論頻率點和頻譜中的峰值判斷軸承和齒輪的健康狀況。計算故障特征成分方程及其對應的相位函數和頻率點是所提算法的關鍵。本文利用轉頻方程、機械系統中各轉頻間的比例參數以及目標軸承故障特征系數對上述參量進行計算。本文方法的流程圖如圖1所示。具體的步驟總結如下：

圖1 算法流程圖Fig.1 Flowchart of proposed diagnosis strategy

步驟1同步測取振動信號和轉速脈沖信號，并根據轉速信號計算轉頻方程。

步驟2基于轉頻方程、機械結構傳動比、目標軸承故障特征系數，計算軸承和齒輪故障特征頻率方程，并根據式(16)計算頻率方程對應的相位函數和頻率點。

步驟3通過故障特征頻率方程計算該故障特征對應的載波矩陣C，并根據載波矩陣分別對原始信號的包絡進行Vold-Kalman濾波，即利用式(11)求取載波矩陣C對應的幅值包絡矩陣A，通過式(12)實現時變非平穩故障特征成分的重構。

步驟4基于各個濾波分量對應的相位函數，分別對濾波分量進行廣義解調，得到解調信號。

步驟5對各個解調成分進行頻譜分析。通過頻譜峰值和理論頻率點的對比，初步確定故障點位置。

步驟6基于初步診斷結果，提取其二倍和三倍諧波成分。通過頻率譜的定量表達進一步驗證步驟5診斷結果的正確性。

3 仿真驗證

為驗證本文算法的有效性，構造了變轉速條件下滾動軸承和齒輪復合故障振動信號。

變轉速故障軸承信號的仿真模型為

(17)

式中：Ai,b為第i個故障沖擊幅值；βb為結構衰減系數；ωb為由軸承引起的共振帶中心頻率；τj為由滾動體滑移帶來的故障沖擊間隔誤差，其取值一般為0.01～0.02；μ(t)為單位階躍函數；ti,b為第i個沖擊發生時刻，可由式(18)計算

(18)

式中：fcc為軸承故障特征系數；fb(t)為軸承轉頻。

變轉速下故障齒輪仿真模型如下

(19)

式中：第一行為齒輪嚙合振動；第二行為由齒輪故障引起的高頻共振；n為齒輪嚙合頻率序號；本文取值1，2和3；An,g為第n個嚙合頻率幅值；z為齒輪齒數；Fg(t)為齒輪轉頻的積分函數；Am,g,βg,ωg和tm,g的含義等同于式(17)中的Ai,b,βb,ωb和ti,b。其中tm,g的計算公式為

(20)

式中：fg(t)為齒輪轉頻。

因此，變轉速條件下軸承和齒輪復合故障振動信號的仿真模型可表示為

x(t)=xb(t)+xg(t)+n(t)

(21)

式中：n(t)為高斯白噪聲。仿真模型的參數見表1。

表1 變轉速下軸承和齒輪復合故障仿真模型參數

根據仿真模型得到的信號的時域表達，如圖2所示。

圖2 仿真信號Fig.2 Simulated signal

根據式(16)計算仿真模型中預設的轉頻方程fb(t)對應的相位函數sb(t)=1.875t2和頻率點fpb=39 Hz。基于轉頻方程fb(t)及其對應的相位函數sp(t)和頻率點fpb、傳動比和軸承故障特征系數分別計算軸承和齒輪故障特征頻率方程及其對應的相位函數和頻率點，詳見表2。此外，為說明本文算法不會因不存在的特征成分，在分析結果中出現相應峰值以造成誤判，隨機選取了故障特征系數fcc1=2.8作為對比分析。

對原始信號進行希爾伯特變換得到包絡信號。首先，分別利用表2中的故障特征頻率方程對包絡信號進行Vold-Kalman濾波，得到四個濾波信號。其次，根據表中的相位函數分別對相應的濾波信號進行廣義解調，得到四個解調信號。再次，對解調信號分別做FFT得到頻率譜依次如圖3中四個小圖所示。最后，通過頻率譜中峰值與表2中頻率點的對比確定故障點位置。圖3中只有兩張圖中存在明顯峰，即只有分別利用表2中第一行和第三行參數獲得的頻率譜中存在明顯峰值，且兩圖中的峰值144.3和25.98分別與表2中頻率點fp0和fpd1近似或相同。因此可以初步判定該仿真模型中軸承存在與故障特征系數3.7關的故障，以及主動輪存在故障。而分別利用表2中第二行和第四行參數獲得的頻率譜中在相應的頻率點fpc和fpd2處并沒有出現峰。說明上述信號中不存在與故障特征系數fcc1=2.8有關的故障特征，而被動輪也不存在故障。

表2 仿真復合故障特征頻率方程及其對應的相位函數和頻率點

圖3 Vold-Kalman廣義解調成分頻率譜Fig.3 Spectra of Vold-Kalman generalized demodulation components

為進一步驗證上述診斷結果，分別對圖3(a)和圖3(c)中峰的二倍及三倍頻進行提取、重置以及定量表達。其中Vold-Kalman濾波用到的頻率方程分別為2f0(t)，3f0(t)，2fd1(t)和3fd1(t)。廣義解調用到的相位函數分別為2s0(t)，3s0(t)，3sd1(t)和3sd1(t)。結合初步驗證結果得到完整頻率譜分別如圖4和圖5所示。圖4中第二、第三個峰值分別為第一個峰值的2倍、3倍。同樣，圖5中的峰值也存在上述比例關系。因此得出結論：仿真信號中軸承存在與故障特征系數3.7相關的故障以及主動齒輪存在故障。診斷結果與仿真模型相一致。

4 實驗驗證

本節在機械故障仿真試驗臺上測取了減速條件下軸承和齒輪復合故障特征振動信號，進一步驗證所提算法的有效性。試驗臺如圖6所示，其中目標軸承存在內圈故障，其故障為人為破壞形成的裂紋故障，故障程度為輕度(寬：0.2 mm 深：0.4 mm)。齒輪箱中被動輪存在斷齒故障，即其中一個齒被人為破壞，出現斷裂。試驗臺機械結構參數見表3。本試驗中采樣頻率設置為48 000 Hz，截取2 s時長的信號(見圖7)進行分析。上述信號對應的轉頻曲線如圖8所示。該轉頻曲線由轉速計測取的軸的轉速脈沖信號計算。對轉頻曲線進行擬合得到擬合方程：fa(t)=-0.042t3+0.029t2-2.767t+30.86。基于式(16)，對上述方程進行改寫，計算其對應的相位函數sa(t)=-0.011t4+0.01t3-1.38t2和頻率點fpa=30.86 Hz。

圖4 軸承頻譜Fig.4 Bearing spectrum

圖5 齒輪頻譜Fig.5 Gear spectrum

根據轉頻地fa(t)及其對應的相位函數sa(t)和頻率點fpa，以及表3中的參數計算所有故障特征成分方程及其相位函數和頻率點如表4所示。

圖6 機械故障仿真試驗臺Fig.6 Machinery fault simulator

參數值軸承類型ER-10k滾珠直徑d/mm7.94節圓直徑D/mm33.5滾珠數量nb8接觸角α0fcco3.052fcci4.95fccb1.99主動輪齒數zg18帶傳動比ib2.56齒輪傳動比ig1.5

首先根據該工況下外圈故障時對應的故障特征頻率方程和相位函數對原始信號的包絡進行處理。即，基于fo(t)對原始信號的包絡進行Vold-Kalman濾波獲取表征軸承外圈故障的時變非平穩成分，基于相位函數so(t)對濾波成分進行廣義解調，最后利用FFT計算解調信號的頻率譜。同樣分別利用表4中內圈、滾珠、主動輪和被動輪故障時對應參數獲得相應的頻率譜。所有頻譜中，只有利用軸承內圈和主動輪故障時對應的參數獲得的頻譜中存在明顯峰，且峰值153.1和8.057分別與頻率點fpi和fpdrive近似，如圖9所示。因此可以初步判定該軸承內圈存在故障，齒輪箱中被動輪存在故障。

圖7 實測振動信號Fig.7 Measured signal

圖8 實測振動信號對應的轉頻曲線Fig.8 Rotational frequency trend

分別利用方程2fi(t)，3fi(t)，2fdriven(t)和相位函數2si(t)，3si(t)，2sdriven(t)對原始信號的包絡進行Vold-Kalman廣義解調。整合上一步驟中獲得的頻率譜，得到軸承和被動輪完整頻率譜分別如圖10和11所示。圖中峰值近似等于頻率點fpi，2fpi，3fpi，fpdriven和2fpdriven，進一步驗證了診斷結果的正確性。

為驗證本文算法在避免能量泄露和消除無關項干擾方面的優越性，利用帶通濾波與廣義解調相結合的算法對本節所測量的原始振動信號進行分析。處理過程如下：首先基于原始信號的時頻表達確定帶通濾波參數提取軸承故障激起的高頻共振帶；其次利用相位函數Si(t)對原始信號的包絡進行廣義解調；最后通過FFT獲得頻率譜如圖12所示。本文方法獲得相對應的頻譜圖如圖9(a)所示。對比分析圖9(a)和圖12可知本文算法獲得頻譜圖更加簡潔、無其他干擾成分。另外，本文算法獲得的頻譜峰值(7.6×106)相比圖12中的(1.4×105)更有優勢。