李康強，馮志鵬，孫宏圖，神克常

(1.魯東大學 蔚山船舶與海洋學院，山東 煙臺 264025；2.北京科技大學 機械工程學院，北京 100083)

行星齒輪箱由于其獨特的結(jié)構和運行方式，相比于其他的機械旋轉(zhuǎn)機構，其振動信號具有明顯的調(diào)幅調(diào)頻特性。此外，實際工程測試中，行星齒輪箱自身的多自由度和環(huán)境中的多振源也造成了信號的復雜非平穩(wěn)多分量特性，常規(guī)方法對于提取行星齒輪箱的故障特征存在困難[1-2]。

目前，基于振動信號分析的行星齒輪箱的故障診斷，大多數(shù)研究集中在振動信號的時域統(tǒng)計指標分析和頻域故障特征提取。Skrimpas等[3]利用不同階嚙合頻率之間的殘差信號的統(tǒng)計特征值檢測行星齒輪與相關軸承的故障。Yoon等[4]利用單壓電應變傳感器的應變信號分析行星齒輪箱并提取故障特征。Wang等[5]利用K均值聚類方法對于不同的故障類型進行分類并判斷設備健康狀態(tài)。Feng等[6]利用本征時間尺度分解結(jié)合幅頻解調(diào)提取行星齒輪箱的故障特征。針對于現(xiàn)有研究，常規(guī)方法在故障診斷分析中對創(chuàng)新性提出了更高的挑戰(zhàn)。

通常行星齒輪箱在運行中，太陽輪以軸向自轉(zhuǎn)作為動力輸入傳遞扭矩，行星輪與太陽輪和固定在箱體上的齒圈同時嚙合，設備運行時多個行星輪既有公轉(zhuǎn)又有自轉(zhuǎn)，種種因素造成行星齒輪箱振動信號的復雜性，同時給故障診斷帶來了困難。在行星齒輪箱的穩(wěn)定運行過程中，隨著齒輪的轉(zhuǎn)動，齒輪對的嚙合位置不斷變化，齒輪的變形量的不斷變化引起的剛度隨之變化，換言之，即系統(tǒng)的瞬時阻尼比是不斷變化的。假想狀態(tài)下的健康齒輪箱運行時，前一齒輪對的嚙合狀態(tài)與后一齒輪對是完全相同的，因此整個持續(xù)變化狀態(tài)具有周期性。假如嚙合齒面發(fā)生損傷時，系統(tǒng)運行會產(chǎn)生額外的能量沖擊，因此，故障狀態(tài)引起的能量變化區(qū)別于正常狀態(tài)的能量變化，而能量變化會反映在瞬時阻尼比的變化上?；诖?，我們提出齒輪齒面故障會引起阻尼比的瞬態(tài)變化，換言之，阻尼比瞬態(tài)變化與故障特征頻率相關。那么如何從振動信號中準確獲取瞬時阻尼比成為行星齒輪箱故障診斷的關鍵問題。

對于阻尼比的研究，一般針對于模態(tài)阻尼比識別方法，如對數(shù)衰減法、ITD(ibrahim time domain)法和STD(spare time domain)法等[7]，對于瞬時阻尼比的研究非常有限。Feldman[8]提出利用Hilbert變換估計瞬時阻尼比的方法。但是，該方法必須綜合計算幅值、頻率和瞬時無阻尼固有頻率，多種類型的誤差因素勢必影響到瞬時阻尼比的準確度。Loh[9]利用經(jīng)驗模式分解(empirical mode decomposition,EMD)與Hilbert譜相結(jié)合來估算瞬時阻尼比，為識別非線性非平穩(wěn)信號提供了一種潛在方法。但是上述研究只是提出估算瞬時阻尼比的方法，并沒有提出具體應用。在現(xiàn)今研究中，尚未發(fā)現(xiàn)利用瞬時阻尼比進行故障診斷的案例。

Zayezdny等[10]提出了基于生成微分方程(generating differential equation,GDE)的信號表示方法，利用信號、信號微分以及信號各種變換組合之間的關系描述信號的結(jié)構和性質(zhì)。李康強等[11-12]通過深入研究發(fā)現(xiàn)，生成微分方程在解調(diào)信號方面有其獨特的優(yōu)勢，并可以應用于旋轉(zhuǎn)機械故障診斷。通過繼續(xù)研究發(fā)現(xiàn)，利用信號的各階微分函數(shù)的非線性組合運算，可以方便高效的計算瞬時阻尼比等動力學參數(shù)信息。但是生成微分方程的非線性組合運算過程需要信號滿足單分量要求，而實際工程測試的振動信號由于振源多雜造成其成分復雜。為了準確提取故障特征，本文發(fā)揮EMD在復雜信號分解方面的自適應性和高計算性，將行星齒輪箱振動信號分解為多個本質(zhì)模式函數(shù)并從中挑選出滿足生成微分方程需要的分量作為敏感單分量，繼而利用生成微分方程準確估計敏感分量的瞬時阻尼比，最后通過識別譜圖中的故障特征頻率來診斷故障。利用不同的仿真信號與試驗實測信號驗證方法的有效性與實用性。

1 故障特征提取方法

1.1 延拓經(jīng)驗模式分解

為滿足生成微分方程單分量的需要和增強故障特征，首先需要將復雜多分量信號進行分解。經(jīng)驗模式分解可以通過“篩選”將包含多個振蕩模式的信號自適應的分解為符合以下條件的本質(zhì)模函數(shù)(intrinsic mode function,IMF)[13]：

(1)整個數(shù)據(jù)長度下，IMF的極值數(shù)量等于過零點的數(shù)量或者至多相差一個；

(2)任意數(shù)據(jù)點，局部極大值的包絡和局部極小值的均值等于0。

經(jīng)驗模式分解的簡化步驟為：

步驟1計算源信號x(t)的所有局部極大值與局部極小值，利用三次樣條插值函數(shù)可以獲得上下包絡線，源信號減去上下包絡線的均值m1可以得到第一個分量h1=x(t)-m1；

步驟2如果上述第一個分量滿足本質(zhì)模函數(shù)的條件，則視為分解出的第一個本質(zhì)模函數(shù)IMF1，如果不滿足，則令源信號x(t)重復步驟1直至計算出的分量滿足條件；

步驟3源信號減去第一個本質(zhì)模函數(shù)得到的余量r1=x(t)-IMF1作為新的源信號重復步驟1和步驟2進行第二次“篩選”過程，重復循環(huán)n次，當殘余量rn成為不能再繼續(xù)分解的單調(diào)函數(shù)時，分解終止，殘余量rn即源信號的趨勢項；

可以看出，經(jīng)驗模式分解滿足自適應性和完備性，而且分解的過程相對簡便高效。但是實際應用中，EMD由于獲取端點極值的擬合包絡時會出現(xiàn)較大跳動，而且隨著分解的進行，這種跳動會逐漸往中間傳遞，造成分解結(jié)果產(chǎn)生端點效應的問題，基于此干擾，我們可以利用信號對稱延拓來解決這個問題。在EMD分解之前，對源信號進行如下對稱延拓

(1)

式中：rnN為源信號x(t)的長度；Ne為延拓值,Nei=1,2,…,Ne。延拓之后數(shù)據(jù)長度增加2Ne個點，所以分解完成后數(shù)據(jù)需要對稱剪切2Ne個點以便恢復原長度。

1.2 瞬時阻尼比計算

生成微分方程可以利用重定義的狀態(tài)函數(shù)對信號進行映射變換，將信號分解為一系列的微分方程的解的合集形式，從而估算信號的動力學參數(shù)。不失一般性，對于一般的調(diào)幅調(diào)頻信號可表示為

x=Acos(ωt+φ)

(2)

式中：A(t)為幅值調(diào)制信號；φ(t)為頻率(相位)調(diào)制信號。

kx-(2δA+δω)δx-[kA-δA(2δA+δω)-ω2]=0

(3)

不失一般性，調(diào)制部分A(t)和φ(t)相對于載波頻率是緩變的，因此可以近似為定值。從而關于其狀態(tài)函數(shù)可以近似δA=0，δω=0，kA=0。繼而推導瞬時頻率為

(4)

系統(tǒng)在振動過程中的能量是周期性持續(xù)衰減的，就振動過程中某一時刻點來說，其能量衰減響應可以表述為形如自由衰減響應形式

x(t)=Ae-ξωntsin(ωdt+φ)

(5)

式中：ξ為阻尼比；ωn為無阻尼固有頻率；ωd為有阻尼固有頻率；φ為初相位。

如上所述，結(jié)合微分與狀態(tài)函數(shù)可以推出式(5)的生成微分方程為

(6)

(7)

對于離散信號而言，實際計算過程可以由對稱差分替代微分，由中心差商替代導數(shù)

(8)

同時離散信號的基本狀態(tài)函數(shù)可以表示為

(9)

代入式(7)可以得到離散信號的瞬時阻尼比序列，其中步長Δt=1。由于計算過程只有數(shù)值的基本運算，因此運算簡便高效。

1.3 特征提取

基于EMD和GDE識別瞬時阻尼比的故障診斷方法的步驟歸納如下，其流程圖如圖1所示。

圖1 流程圖

步驟1對待分析信號進行對稱延拓；

步驟2利用經(jīng)驗模式分解延拓信號，得到數(shù)個本質(zhì)模函數(shù)；

步驟3基于GDE計算各個本質(zhì)模函數(shù)的瞬時頻率，并選擇最先分解出的且瞬時頻率在嚙合頻率或其諧波上下波動的本質(zhì)模函數(shù)作為敏感分量進行下一步分析[14]；

步驟4基于GDE計算敏感分量的瞬時阻尼比，并對稱剪切恢復至原信號長度；

步驟5對瞬時阻尼比進行Fourier變換，得到瞬時阻尼比頻譜；

步驟6根據(jù)阻尼比頻譜峰值頻率匹配故障特征頻率，定位故障。

2 仿真驗證

2.1 脈沖信號仿真

為直觀顯示瞬時阻尼比變化特征和故障特征，我們首先仿真建模一系列等衰減等間隔脈沖信號。建模生成

u(t-nT)+n(t)

(10)

圖2(a)為仿真信號時域波形，可以看出時域波形為一系列等間隔等衰減的脈沖信號，其中間隔為T=0.1 s，一個間隔周期內(nèi)可以發(fā)現(xiàn)信號從最大幅值衰減為0并且持續(xù)保持一段時間，因此瞬時阻尼比應為一段時間內(nèi)為恒定值然后下降到0值。圖2(b)驗證了瞬時阻尼比的波形結(jié)構為方波，方波頂部為預設定阻尼比ξ=0.06，方波底部為0值。圖2(c)為瞬時阻尼比頻譜，可以明顯看到突出的峰值頻率對應于設定的故障特征頻率及其倍頻kfg，k=1,2,…。仿真結(jié)果證明了通過瞬時阻尼比可以提取故障特征頻率。

2.2 齒輪箱信號仿真

根據(jù)行星齒輪箱的振動機理，其信號可以用調(diào)幅-調(diào)頻信號模型描述[15]。不失一般性，只考慮齒輪嚙合頻率和單一故障齒輪的特征頻率，則振動信號仿真形式設為

x(t)=[1-cos(2πfct)][1+Acos(2πfgt)]×

cos[2πfmt+Bsin(2πfgt)+φ]+n(t)

(11)

式中：fc為行星架旋轉(zhuǎn)頻率；fg為太陽輪故障特征頻率；fm為齒輪嚙合頻率。設置fc=1.5 Hz，fg=27 Hz，fm=200 Hz，調(diào)幅系數(shù)A和調(diào)頻系數(shù)B為1，初相位φ=0，仿真信號采樣頻率為1 500 Hz，加入信噪比為n(t)=6 dB的Gauss白噪聲。

信號分析如圖3所示。同脈沖信號分析流程一致，從圖3(b)細節(jié)放大驗證了瞬時阻尼比的突出幅值間隔為故障周期T=1/fg≈0.037 s。從圖3(c)的瞬時阻尼比頻譜中可以很明顯的得到突出頻率幅值匹配故障特征頻率fg及其倍頻，在噪聲干擾下頻譜底部細節(jié)放大圖可以匹配行星架旋轉(zhuǎn)頻率fc及其倍頻。通過仿真驗證方法可以有效的提取故障，且抗噪性也較為明顯。

圖3 齒輪箱仿真信號分析

3 試驗分析

1.驅(qū)動電機；2.信號調(diào)節(jié)器；3.轉(zhuǎn)矩和速度傳感器；4.行星齒輪箱；5.信號采集器；6.CZ-5型磁粉制動器。

圖5 齒輪模擬故障

3.1 健康信號分析

圖6 正常信號分析結(jié)果

3.2 太陽輪故障信號分析

圖7 太陽輪故障信號分析結(jié)果

3.3 齒圈故障信號分析

由于三個行星輪同時與齒圈嚙合，因此當齒圈齒面出現(xiàn)故障的時候，每個運轉(zhuǎn)周期會出現(xiàn)三次嚙合(三次故障脈沖)，所以我們需要重點關注齒圈故障特征頻率的1/3倍頻的幅值情況，齒圈故障信號的分析圖，如圖8所示。按照分析流程我們依然選擇第一個IMF作為敏感分量。從圖8(d)阻尼比頻譜中可以看出，與正常信號(深色虛線)相比，淺色實線代表的齒圈故障信號顯示出較多的峰值。這些峰值都與齒圈的故障特征頻率fr及其1/3倍頻相關，例如1/3fr，2/3fr，4/3fr，2fr，3fr等較高的突出幅值。通過匹配齒圈故障特征頻率，可以驗證齒圈出現(xiàn)了故障，與此同時頻譜中并沒有出現(xiàn)與太陽輪故障特征頻率和行星輪故障特征頻率的突出幅值，驗證了其他齒輪的健康狀態(tài)。

圖8 齒圈故障信號分析結(jié)果

4 結(jié) 論

瞬時阻尼比能夠反映系統(tǒng)能量的變化狀態(tài)，而且波動特征對應于系統(tǒng)的特征頻率。生成微分方程可以利用各階次微分函數(shù)的非線性組合計算系統(tǒng)的瞬時阻尼比?；诖?，本文提出結(jié)合經(jīng)驗模式分解和生成微分方程的行星齒輪箱故障診斷方法。本方法能夠簡便高效的計算瞬時阻尼比，并通過瞬時阻尼比頻譜可以匹配故障特征頻率，準確定位行星齒輪箱的故障。通過仿真信號驗證和試驗實測信號對比分析，證明了該方法在提取故障特征上的有效性和實用性。