基于擴展卡爾曼濾波和希爾伯特-黃變換瞬時頻率的齒輪嚙合剛度辨識算法

闕紅波， 高 揚， 吳成攀， 栗 華

(中車戚墅堰機車車輛工藝研究所有限公司， 常州 213000)

齒輪傳動是機械傳動中應(yīng)用最廣的一種傳動形式，該傳動形式具有傳動比精確、效率高、結(jié)構(gòu)緊湊等特點。由于齒輪在嚙合時參與嚙合的齒數(shù)周期變化，因此齒輪副的嚙合剛度具有時變特性。時變嚙合剛度會對齒輪的支承結(jié)構(gòu)施加動剛度激勵，產(chǎn)生振動，振動的大小和形式與嚙合狀態(tài)有關(guān)[1]，因此通過振動響應(yīng)可以辨識齒輪的嚙合剛度。由于齒輪的傳遞誤差會影響嚙合剛度，因此通過監(jiān)測齒輪的嚙合剛度能夠控制傳遞誤差，提高傳動效率[2]。目前通過系統(tǒng)輸出響應(yīng)辨識系統(tǒng)參數(shù)的方法有遞推最小二乘算法[3-4]、擴展卡爾曼濾波算法[5-6]、遺傳算法[7]、神經(jīng)網(wǎng)絡(luò)算法[8]等。

楊恒等[9]通過擴展卡爾曼濾波(extended Kalman filter，EKF)參數(shù)識別算法對高速列車車下懸掛橡膠彈簧的阻尼和剛度進行了辨識，針對卡爾曼濾波算法跟蹤速度不足的問題，提出了自適應(yīng)強跟蹤算法，結(jié)果表明提出的算法能夠快速對參數(shù)進行跟蹤。李亞偉等[10]應(yīng)用擴展卡爾曼濾波器，結(jié)合機匣上的振動信號，采用衰減記憶濾波的方法實現(xiàn)了轉(zhuǎn)子-支承-機匣模型的參數(shù)識別，準確識別了系統(tǒng)的不對中和不平衡參數(shù)。鄒玥等[11]根據(jù)擴展卡爾曼濾波算法提出了一種磁懸浮軸承剛度和阻尼的參數(shù)識別方法，辨識結(jié)果表明激勵頻率在200 Hz以內(nèi)時辨識準確率達到90%以上，但誤差率隨轉(zhuǎn)速上升而上升。

由于轉(zhuǎn)軸的瞬時旋轉(zhuǎn)速度在很大程度上反映了傳動系統(tǒng)的工作狀態(tài)，因此一些研究人員開展了使用振動信號估計轉(zhuǎn)軸瞬時旋轉(zhuǎn)頻率的研究?？档碌萚12]利用變分模態(tài)分解(variational mode decomposition，VMD)變換對振動信號進行降噪重構(gòu)，然后對重構(gòu)信號頻譜進行Viterbi瞬時頻率估計，該方法提高了瞬頻的估計精度。陳建新等[13]結(jié)合希爾伯特變換和傅里葉變換微分性質(zhì)提出了一種針對聲信號的瞬時頻率特征提取方法。程衛(wèi)東等[14]采用經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)降噪方法優(yōu)化了基于瞬時故障特征頻率的滾動軸承瞬時轉(zhuǎn)頻估計方法[15]。希爾伯特-黃變換(Hilbert Huang transform，HHT)結(jié)合經(jīng)驗?zāi)B(tài)分解算法和希爾伯特(Hilbert)變換將信號分解為多個本征模態(tài)函數(shù)(intrinsic mode function，IMF)，每個IMF滿足非負頻率計算條件，通過Hilbert變換計算各IMF的瞬時頻率[16]。梁明等[17]對含噪聲的瞬變電磁信號進行了HHT分析，獲取了信號的瞬時幅度和瞬時頻率，結(jié)果表明HHT提取信號瞬時頻率分辨率較高，優(yōu)于傳統(tǒng)分析方法。胡杰等[18]對內(nèi)燃機瞬時轉(zhuǎn)速信號進行了EMD分解，將分解后的各個IMF分量作了HHT變換，得到了各IMF分量的瞬時頻率，研究結(jié)果表明HHT方法能夠有效分離瞬時轉(zhuǎn)速信號中的各個頻段的信息，并且各IMF分量都有其物理意義。

現(xiàn)通過齒輪傳動系統(tǒng)振動響應(yīng)信號的HHT瞬時頻率積分結(jié)果，應(yīng)用EKF，對齒輪嚙合動剛度進行跟蹤辨識。

1 齒輪傳動動力學(xué)離散模型

假設(shè)主動齒輪和從動齒輪僅存在旋轉(zhuǎn)自由度，典型的齒輪傳動動力學(xué)模型如圖1所示。

φ為主動齒輪的旋轉(zhuǎn)角度，θ為從動齒輪的旋轉(zhuǎn)角度；r1為主動齒輪節(jié)圓半徑；r2為從動齒輪節(jié)圓半徑；K(t)為齒輪副之間的嚙合剛度圖1 齒輪嚙合動力學(xué)模型Fig.1 Dynamic model of gear meshing

由于齒輪副在不同的嚙合角度下，參與嚙合的輪齒對數(shù)不同，因此嚙合剛度為一時變參數(shù)。由于齒輪副的嚙合阻尼比較小，阻尼比的取值范圍一般在0.005～0.075，因此模型中忽略了嚙合阻尼比。圖1所示的動力學(xué)系統(tǒng)微分方程為

(1)

式(1)中：J1為主動齒輪的轉(zhuǎn)動慣量；J2為從動齒輪的轉(zhuǎn)動慣量。假設(shè)初始狀態(tài)下齒輪靜止，通過拉氏變換可以推導(dǎo)系統(tǒng)的傳遞函數(shù)為

(2)

應(yīng)用沖擊響應(yīng)不變法，可以將系統(tǒng)傳遞函數(shù)由拉式變換域轉(zhuǎn)換到z變換域，系統(tǒng)的兩個極點為

(3)

式(3)中：j表示虛數(shù)。應(yīng)用部分分式展開法可以將傳遞函數(shù)化為

(4)

轉(zhuǎn)換后，傳遞函數(shù)的系統(tǒng)傳遞函數(shù)為

(5)

利用z變換的性質(zhì)，可以將連續(xù)域的傳遞函數(shù)轉(zhuǎn)換為離散差分方程形式，以便于在計算機中通過離散點的數(shù)值進行參數(shù)辨識，系統(tǒng)的離散模型為

θ(n)=Aθ(n-1)-θ(n-2)+Bφ(n-1)

(6)

式(6)中：

(7)

(8)

(9)

在齒輪故障的初期，嚙合動剛度的變化量相對較小，但是隨著齒輪故障的加劇，齒輪嚙合剛度的變化將會越發(fā)明顯。不同類型故障對嚙合剛度的影響不同，所有故障類型中，對嚙合剛度影響最大的是齒根裂紋，而齒面點蝕、齒形誤差等其他類型故障對剛度的影響在2%以內(nèi)[19]。根據(jù)齒輪時變嚙合剛度的計算理論[20]，對于一般的齒輪，當齒根裂紋長度分別為10、20、30 mm時，裂紋深度為3 mm時，含裂紋的輪齒通過嚙合區(qū)間時齒輪副嚙合剛度如圖2所示。

當齒根裂紋長度為30 mm，裂紋深度分別為1、2、3 mm時，含裂紋的輪齒通過嚙合區(qū)間時齒輪副嚙合剛度如圖3所示。

由圖2和圖3可以看出，齒根裂紋對嚙合剛度的影響較大，不同裂紋長度和深度造成剛度減小的程度不同。不同裂紋長度和深度下齒輪嚙合剛度最大減小量如圖4所示。

圖2 不同裂紋長度下嚙合剛度Fig.2 Meshing stiffness under crack of different length

圖3 不同裂紋深度下嚙合剛度Fig.3 Meshing stiffness under crack of different depth

圖4 不同裂紋長度和深度下剛度最大減小百分比Fig.4 Maximum percentage reduction of stiffness under different crack length and depth

圖4表明，裂紋深度的增加更能夠?qū)е聡Ш蟿偠鹊南陆?，裂紋深度在1 mm以內(nèi)且裂紋長度在10 mm以下時，最大嚙合剛度的下降率在1%以下，當裂紋深度達到3 mm且裂紋長度達到30 mm時，最大嚙合剛度的下降率約為10%。由于輪齒在含有不同深度和不同長度的裂紋時最大剛度的下降量存在差異，因此通過比較齒輪系統(tǒng)正常和當前狀態(tài)下的最大嚙合剛度差異就能夠估計輪齒是否存在裂紋，以及齒根裂紋的深度和長度。

2 嚙合剛度識別算法

2.1 EKF算法

考慮如下離散系統(tǒng)，方程(10)表示系統(tǒng)的狀態(tài)方程，方程(11)表示系統(tǒng)的觀測方程。

Xk=f(Xk-1,Vk-1)+uk+Vk-1

(10)

Yk=h(Xk,nk)+nk

(11)

式中：V、n分別為系統(tǒng)輸入噪聲和觀測噪聲；u為系統(tǒng)的輸入量；X為系統(tǒng)狀態(tài)參數(shù)；Y為觀測值；f(·)和h(·)分別為系統(tǒng)的狀態(tài)轉(zhuǎn)移矩陣和觀測矩陣。首先考慮f(·)為線性方程，在將狀態(tài)向量增廣后，狀態(tài)一步更新方程為

(12)

(13)

(14)

誤差矩陣的一步更新算法為

Pk,k-1=FPFT+Q

(15)

式(15)中：

(16)

(17)

(18)

式(18)中：Qw=E(Wk·Wj)，Qn=E(nk·nj)，為過程噪聲的協(xié)方差?？柭鼮V波器使用預(yù)測誤差來修正原來的狀態(tài)估計結(jié)果，因此更新后的狀態(tài)向量為

(19)

誤差協(xié)方差矩陣更新算法為

(20)

(21)

若系統(tǒng)的狀態(tài)轉(zhuǎn)移矩陣中含有非線性項，將f(x)用泰勒級數(shù)展開并忽略2階及以上項可以將f(x)化為線性函數(shù)，此時得到的濾波值為近似值。利用式(12)、式(15)、式(19)、式(20)以及式(21)，可以通過系統(tǒng)的輸入和輸出信號對系統(tǒng)參數(shù)進行跟蹤辨識。

2.2 瞬時頻率計算

由于動力學(xué)方程(6)中的輸入、輸出參數(shù)皆為轉(zhuǎn)軸的轉(zhuǎn)動角度，而一般測試系統(tǒng)中對轉(zhuǎn)軸旋轉(zhuǎn)角度的測量難度較大，由于旋轉(zhuǎn)角度是轉(zhuǎn)軸轉(zhuǎn)速的積分，因此可以利用對信號瞬時頻率的積分實現(xiàn)旋轉(zhuǎn)角度的計算。Hilbert變換可以求解實信號的共軛部分從而構(gòu)成解析信號，通過對解析信號的相位角求導(dǎo)就可以得到信號的瞬時頻率，然而一般的信號不滿足帶限要求，因此需要使用HHT算法將系統(tǒng)響應(yīng)信號進行處理。HHT算法使用EMD算法將信號分解成若干個IMF，各個IMF滿足帶限要求，對IMF進行Hilbert變換得到的瞬時頻率不會出現(xiàn)負數(shù)，因此通過此方法可以獲得信號的瞬時頻率。EMD分解過程為

(22)

式中：s(t)為原始信號；rn(t)為EMD分解殘差；Ci為分解得到的各階IMF信號。對各IMF進行希爾伯特(Hilbert)變換，可以求取信號的瞬時相位角，希爾伯特變換的頻域濾波器為

(23)

信號的希爾伯特(Hilbert)變換通過Hilbert濾波器在頻域與信號的頻域表達相乘實現(xiàn)，即

CH(t)=F-1[H(f)C(f)]

(24)

式(24)中：C(f)為IMF信號的頻譜；F-1表示傅里葉逆變換；CH(t)為信號C(t)的希爾伯特變換信號。信號s(t)和sH(t)共同構(gòu)成信號s(t)的解析信號z(t)。將解析信號z(t)用歐拉方程寫作指數(shù)形式為

z(t)=a(t)ejφ(t)

(25)

式(25)中：a(t)為解析信號的幅度；φ(t)為解析信號的瞬時相位角。信號s(t)的瞬時頻率定義為

(26)

3 齒輪嚙合剛度識別

3.1 仿真分析

假設(shè)主動齒輪以20 r/s的速度旋轉(zhuǎn)，齒輪副的物理參數(shù)如表1所示。通過式(5)計算得到的主、從動齒輪的轉(zhuǎn)動角度曲線，如圖5所示。

表1 齒輪副參數(shù)

圖5 主、從動齒輪旋轉(zhuǎn)角度曲線Fig.5 Angle curve of driving and driven gears

由圖5可知，從動齒輪的旋轉(zhuǎn)角度隨著主動齒輪旋轉(zhuǎn)角度增加而上升，由于主、從動齒輪的節(jié)圓半徑比例為1∶2，因此從動齒輪的旋轉(zhuǎn)角度為主動齒輪的一半。將主、從動齒輪旋轉(zhuǎn)角度分別作為輸入信號和觀測信號，使用EKF算法進行辨識，辨識結(jié)果如圖6所示。

利用式(9)可以求得辨識結(jié)果對應(yīng)的齒輪嚙合剛度，在0.1 s內(nèi)，動剛度的辨識結(jié)果如圖7所示。

由圖7可知，辨識剛度結(jié)果與仿真剛度并不完全重合，但在大部分時間點兩者幅值接近。對仿真剛度曲線和辨識剛度曲線進行頻譜分析，計算結(jié)果如圖8所示。

結(jié)果表明辨識剛度頻譜峰值所在頻率與仿真剛度的預(yù)設(shè)變化頻率相同，為300 Hz。辨識剛度曲線頻譜在300 Hz處幅值為1.96×105N/m，仿真剛度曲線頻譜在300 Hz處幅值為2×105N/m，辨識準確率為98%。因此在仿真條件下，擴展卡爾曼濾波算法能夠很好地跟蹤嚙合剛度的變化，辨識準確率較高。

3.2 實測信號分析

通過在某型高速列車中使用的齒輪箱上安裝振動加速度傳感器，測試了轉(zhuǎn)速在5 400 r/min下齒輪旋轉(zhuǎn)產(chǎn)生的振動加速度，現(xiàn)場測試情況如圖9所示。被測齒輪箱中齒輪的齒數(shù)比為73∶29，傳動比為2.517，測試系統(tǒng)采樣頻率為20 000 Hz，齒輪副的基本參數(shù)符號如表2所示。

圖6 參數(shù)辨識結(jié)果Fig.6 Parameters identification results

圖7 動剛度辨識結(jié)果Fig.7 Identification results of dynamic stiffness

圖8 模擬剛度曲線與辨識剛度曲線頻譜Fig.8 Spectra of simulation stiffness curve and identified stiffness curve

圖9 齒輪箱測試現(xiàn)場Fig.9 Gearbox test site

表2 測試齒輪參數(shù)

基圓半徑等于分度圓半徑乘以壓力角的余弦值。當變位系數(shù)Xn1大等于Xn2，且0.5≤Xn1+Xn2≤2時，單齒的節(jié)點嚙合剛度[21]為

KC=0.8×103b/q

(27)

式(27)中：b為齒輪的齒寬;q的計算式為

(28)

假設(shè)齒輪的重合度為ε，單齒在0～ε范圍內(nèi)的嚙合剛度為

K(x)=Ax2+Bx+C,x∈[0,ε]

(29)

式(29)中：A=-2K′/ε2；B=2K′/ε2；C=K′；K′=KC/1.5。當式(29)中x等于1時，表示下一對輪齒進入嚙合。根據(jù)齒輪箱內(nèi)齒輪的重合度，可以計算得到主動齒輪在360°旋轉(zhuǎn)范圍內(nèi)，齒輪副的理論嚙合剛度，360°范圍內(nèi)的齒輪副嚙合剛度曲線如圖10所示。

圖10 齒輪理論嚙合剛度Fig.10 Theoretical meshing stiffness of gear

齒輪箱測試振動加速度信號波形如圖11所示。

圖11 齒輪箱測試加速度Fig.11 Measured acceleration of gearbox

使用EMD算法對測得的加速度信號進行分解，分解得到的第1到第5個IMF如圖12所示。

圖12 EMD算法分解結(jié)果Fig.12 EMD algorithm decomposition results

使用Hilbert變換計算各IMF的瞬時頻率，瞬時頻率曲線如圖13所示，各個IMF的平均頻率如表3所示。

根據(jù)齒輪副的齒數(shù)比，第2和第4個IMF的平均頻率分別對應(yīng)主、從動齒輪的旋轉(zhuǎn)頻率，因此第2和第4個IMF對應(yīng)主、從動齒輪旋轉(zhuǎn)引起的振動。使用HHT瞬時頻率計算方法可以計算主、從動齒輪旋轉(zhuǎn)對應(yīng)IMF分量的瞬時頻率，計算得到的瞬時頻率即主、從動齒輪的瞬時轉(zhuǎn)速。使用梯形數(shù)值積分方法可以將轉(zhuǎn)速信息轉(zhuǎn)換為旋轉(zhuǎn)角度信息，以便于系統(tǒng)辨識。使用數(shù)值積分后，主、從動齒輪的旋轉(zhuǎn)角度曲線如圖14所示。

圖14 齒輪旋轉(zhuǎn)角度曲線Fig.14 Rotating angle of gears

使用積分后得到的主、從動齒輪旋轉(zhuǎn)角度信號進行嚙合動剛度辨識，通過式(29)計算得到的理論嚙合剛度曲線以及使用HHT瞬時頻率積分結(jié)果得到的辨識剛度曲線如圖15所示。

圖15 剛度辨識結(jié)果Fig.15 Identification results of stiffness

計算結(jié)果表明基于HHT瞬時頻率的EKF算法能夠跟蹤齒輪剛度的變化，但辨識結(jié)果與理論計算結(jié)果仍有一定差異，最大差異約為8.6×105N/m。對圖15中的辨識結(jié)果進行頻譜分析，分析結(jié)果如圖16所示。

圖16 理論計算嚙合剛度與辨識嚙合剛度頻譜對比Fig.16 Spectrum comparison of theoretical calculating meshing stiffness and identified meshing stiffness

由圖16可知，理論計算嚙合剛度的頻譜峰值所在頻率位置與辨識嚙合剛度頻譜峰值所在頻率位置相同，都為2 610 Hz，該頻率與輸入軸在5 400 r/min轉(zhuǎn)速下的嚙合頻率相同。在嚙合頻率處，辨識嚙合剛度的頻譜幅值為2.81×107N/m，理論計算嚙合剛度的頻譜幅值為2.58×107N/m，辨識準確率約為91%。

5 結(jié)論

應(yīng)用沖擊響應(yīng)不變法對典型齒輪副傳動系統(tǒng)進行了離散差分方程建模，得到了針對離散信號的時間序列差分方程。使用EKF算法和HHT瞬時頻率積分結(jié)果實現(xiàn)了齒輪嚙合剛度的跟蹤辨識，研究得出以下結(jié)論。

(1)使用HHT算法能夠從系統(tǒng)振動響應(yīng)信號中提取主、從動齒輪的旋轉(zhuǎn)瞬時頻率，將瞬時頻率的數(shù)值積分結(jié)果作為辨識系統(tǒng)的輸入?yún)?shù)，使用EKF能夠有效識別齒輪嚙合剛度。

(2)使用齒輪箱振動加速度信號，基于HHT瞬時頻率和EKF算法的辨識算法能夠有效地跟蹤齒輪嚙合剛度變化，辨識得到的嚙合頻率與實際齒輪嚙合頻率相同，但辨識結(jié)果與理論計算結(jié)果存在一定誤差。