基于Vold-Kalman廣義解調(diào)的變轉(zhuǎn)速軸承和齒輪復(fù)合故障診斷

趙德尊, 李建勇, 程衛(wèi)東, 溫偉剛

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

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

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

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

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

Vold-Kalman階比跟蹤是Vold等[11]基于Kalman濾波器提出的可以用于估計振動信號中某一諧波成分的算法。與傳統(tǒng)的帶通濾波或經(jīng)驗?zāi)Ｊ椒纸馑惴ㄏ啾?，該算法可在時域中直接提取某一特定的近似、交叉的諧波分量，避免了由時域至頻域變換帶來的相位偏差和能量泄露[11-12]，因此受到了諸多學(xué)者關(guān)注。其基本原理從第一代發(fā)展到第二代，具備了成熟的理論體系，并廣泛的應(yīng)用于商業(yè)軟件中。然而現(xiàn)有的Vold-Kalman濾波算法均是與階比分析相結(jié)合處理變轉(zhuǎn)速旋轉(zhuǎn)機械振動信號[13-14]。廣義解調(diào)算法是Olhede等[15]提出的一種可以將時頻分布是傾斜、非線性的非平穩(wěn)信號轉(zhuǎn)換成時頻分布是線性的且平行于時間軸的平穩(wěn)信號的分析方法?；趶V義解調(diào)算法對于調(diào)幅-調(diào)頻信號的處理優(yōu)勢，很多學(xué)者將其應(yīng)用于旋轉(zhuǎn)機械設(shè)備的故障診斷[16]。

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

1 理論基礎(chǔ)

1.1 Vold-Kalman濾波

旋轉(zhuǎn)機械中產(chǎn)生的調(diào)制信號可表示為幅值和載波信號的乘積

(1)

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

(2)

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

(3)

設(shè)定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)，系統(tǒng)的觀測方程可表示為

(8)

其矩陣形式為

Y-CA=δ

(9)

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

(10)

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

進一步求得目標(biāo)分量的幅值包絡(luò)矩陣A

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

(11)

通過載波矩陣C和幅值包絡(luò)A即可重構(gòu)目標(biāo)分量

X=AC

(12)

1.2 廣義解調(diào)變換

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

(13)

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

(14)

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

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

(15)

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

(16)

式中：f(t)為特定時變非平穩(wěn)成分的頻率方程；f0為頻率方程對應(yīng)的頻率點；s0(t)為頻率方程對應(yīng)的相位函數(shù)；ds0(t)/dt為對s0(t)求導(dǎo)。

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

2 本文算法

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

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

步驟1同步測取振動信號和轉(zhuǎn)速脈沖信號，并根據(jù)轉(zhuǎn)速信號計算轉(zhuǎn)頻方程。

步驟2基于轉(zhuǎn)頻方程、機械結(jié)構(gòu)傳動比、目標(biāo)軸承故障特征系數(shù)，計算軸承和齒輪故障特征頻率方程，并根據(jù)式(16)計算頻率方程對應(yīng)的相位函數(shù)和頻率點。

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

步驟4基于各個濾波分量對應(yīng)的相位函數(shù)，分別對濾波分量進行廣義解調(diào)，得到解調(diào)信號。

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

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

3 仿真驗證

為驗證本文算法的有效性，構(gòu)造了變轉(zhuǎn)速條件下滾動軸承和齒輪復(fù)合故障振動信號。

變轉(zhuǎn)速故障軸承信號的仿真模型為

(17)

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

(18)

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

變轉(zhuǎn)速下故障齒輪仿真模型如下

(19)

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

(20)

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

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

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

(21)

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

表1 變轉(zhuǎn)速下軸承和齒輪復(fù)合故障仿真模型參數(shù)

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

圖2 仿真信號Fig.2 Simulated signal

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

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

表2 仿真復(fù)合故障特征頻率方程及其對應(yīng)的相位函數(shù)和頻率點

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

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

4 實驗驗證

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

圖4 軸承頻譜Fig.4 Bearing spectrum

圖5 齒輪頻譜Fig.5 Gear spectrum

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

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

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

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

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

圖8 實測振動信號對應(yīng)的轉(zhuǎn)頻曲線Fig.8 Rotational frequency trend

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

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

表4 實測復(fù)合故障特征頻率方程及其對應(yīng)的相位函數(shù)和頻率點

圖9 濾波信號頻率譜Fig.9 Spectra of filter signal

圖10 軸承頻率譜Fig.10 Bearing spectrum

圖11 齒輪頻率譜Fig.11 Gear spectrum

圖12 基于Si(t)的軸承廣義解調(diào)頻率譜Fig.12 Spectrum of generalized demodulated signal by Si(t)

5 結(jié) 論

本文提出了基于Vold-Kalman廣義解調(diào)的變轉(zhuǎn)速軸承和齒輪復(fù)合故障特征提取方法。研究得出以下結(jié)論：

(1)通過轉(zhuǎn)頻方程預(yù)設(shè)故障特征頻率方程，直接利用Vold-Kalman濾波算法實現(xiàn)特定的時變非平穩(wěn)故障特征成分的提取。無需計算最優(yōu)濾波參數(shù)以提取原始信號中由故障激起的高頻共振帶。

(2)計算故障特征頻率方程對應(yīng)的相位函數(shù)和頻率點，利用廣義解調(diào)變換將提取的時變非平穩(wěn)故障特征成分進行平穩(wěn)化重置。無需角域重采樣即可消除轉(zhuǎn)速變化的影響。

(3)基于Vold-Kalman 廣義解調(diào)的復(fù)合故障特征提取方法直接在時域提取、平穩(wěn)化重置特定的故障特征成分，避免了能量的泄露。同時利用快速傅里葉變換實現(xiàn)故障特征的定量表達，使得頻率譜更簡潔，增加了易讀性。