闕紅波， 高 揚(yáng)， 吳成攀， 栗 華

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

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

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

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

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

(6)

式(6)中：

(7)

(8)

(9)

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

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

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

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

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

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

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

2 嚙合剛度識(shí)別算法

2.1 EKF算法

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

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

(10)

Yk=h(Xk,nk)+nk

(11)

式中：V、n分別為系統(tǒng)輸入噪聲和觀測(cè)噪聲；u為系統(tǒng)的輸入量；X為系統(tǒng)狀態(tài)參數(shù)；Y為觀測(cè)值；f(·)和h(·)分別為系統(tǒng)的狀態(tài)轉(zhuǎn)移矩陣和觀測(cè)矩陣。首先考慮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)，為過(guò)程噪聲的協(xié)方差。卡爾曼濾波器使用預(yù)測(cè)誤差來(lái)修正原來(lái)的狀態(tài)估計(jì)結(jié)果，因此更新后的狀態(tài)向量為

(19)

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

(20)

(21)

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

2.2 瞬時(shí)頻率計(jì)算

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

(22)

式中：s(t)為原始信號(hào)；rn(t)為EMD分解殘差；Ci為分解得到的各階IMF信號(hào)。對(duì)各IMF進(jìn)行希爾伯特(Hilbert)變換，可以求取信號(hào)的瞬時(shí)相位角，希爾伯特變換的頻域?yàn)V波器為

(23)

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

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

(24)

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

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

(25)

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

(26)

3 齒輪嚙合剛度識(shí)別

3.1 仿真分析

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

表1 齒輪副參數(shù)

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

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

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

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

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

3.2 實(shí)測(cè)信號(hào)分析

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

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

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

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

圖9 齒輪箱測(cè)試現(xiàn)場(chǎng)Fig.9 Gearbox test site

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

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

KC=0.8×103b/q

(27)

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

(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。當(dāng)式(29)中x等于1時(shí)，表示下一對(duì)輪齒進(jìn)入嚙合。根據(jù)齒輪箱內(nèi)齒輪的重合度，可以計(jì)算得到主動(dòng)齒輪在360°旋轉(zhuǎn)范圍內(nèi)，齒輪副的理論嚙合剛度，360°范圍內(nèi)的齒輪副嚙合剛度曲線如圖10所示。

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

齒輪箱測(cè)試振動(dòng)加速度信號(hào)波形如圖11所示。

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

使用EMD算法對(duì)測(cè)得的加速度信號(hào)進(jìn)行分解，分解得到的第1到第5個(gè)IMF如圖12所示。

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

使用Hilbert變換計(jì)算各IMF的瞬時(shí)頻率，瞬時(shí)頻率曲線如圖13所示，各個(gè)IMF的平均頻率如表3所示。

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

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

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

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

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

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

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

5 結(jié)論

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

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

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