基于非線性短時傅里葉變換階次跟蹤的變速行星齒輪箱故障診斷

王友仁 王 俊 黃海安

南京航空航天大學(xué)自動化學(xué)院，南京，211106

0 引言

行星齒輪箱廣泛用于直升機(jī)主減速器、風(fēng)力發(fā)電機(jī)組等。行星齒輪箱常工作在復(fù)雜多變的環(huán)境下，其太陽輪、行星輪、齒圈等部件的故障發(fā)生概率高、易損壞［1］。行星齒輪箱發(fā)生故障時，振動信號的故障特征微弱、非平穩(wěn)、有噪聲與干擾，不僅受故障、多個傳遞路徑引起的調(diào)頻、調(diào)幅和調(diào)相作用，還受到轉(zhuǎn)速變化引起的調(diào)制、多個激勵振動源之間相互耦合作用，使得傳統(tǒng)的信號頻譜分析技術(shù)難以提取有效故障特征［2］。

變轉(zhuǎn)速下旋轉(zhuǎn)機(jī)械振動信號分析主要采用無鍵相階次跟蹤方法（non?bonding phase order tracking method，NPOTT），通過角度域等間隔采樣技術(shù)將時間域的非平穩(wěn)信號轉(zhuǎn)化為角度域的平穩(wěn)或循環(huán)平穩(wěn)信號，消除轉(zhuǎn)速波動帶來的頻率模糊現(xiàn)象。李蓉等［3］提出基于線調(diào)頻小波路徑追蹤（chirplet path pursuit，CPP）算法［4］與總體平均經(jīng)驗?zāi)B(tài)分解（ensemble empirical mode decomposi?tion，EEMD）的齒輪箱復(fù)合故障診斷方法，用CPP算法得到振動信號的轉(zhuǎn)頻曲線，對轉(zhuǎn)頻信息重采樣并進(jìn)行EEMD分解，獲得了變轉(zhuǎn)速齒輪箱復(fù)合故障特征，但CPP算法復(fù)雜、效率低。RODOPOULOS等［5］提出基于諧波信號分解的瞬時轉(zhuǎn)速估計方法，將瞬時轉(zhuǎn)速估計問題轉(zhuǎn)化為信號模型特征值求取問題，但該方法對噪聲敏感，不適用于強(qiáng)噪聲情況下的瞬時轉(zhuǎn)速估計。郭瑜等［6］提出一種基于瞬時頻率（instantaneous frequency，IF）估計的旋轉(zhuǎn)機(jī)械階次跟蹤方法，利用短時傅里葉變換（short?time Fourier transform，STFT）譜峰值進(jìn)行振動信號瞬時頻率估計，該方法在瞬時頻率變化較小情況下結(jié)果較好，但不適用于瞬時頻率變化劇烈的情況。非線性短時傅里葉變換（non?liner short?time Fouri?er transform，NLSTFT）方法［7］能準(zhǔn)確估計振動信號瞬時頻率，對瞬時頻率曲線積分獲得瞬時相位曲線，通過重采樣完成階次跟蹤過程。NLSTFT方法在大轉(zhuǎn)速變化、強(qiáng)噪聲情況下效果較好，值得進(jìn)一步深入研究。

變分模態(tài)分解（variational mode decomposi?tion，VMD）［8］是一種新的非平穩(wěn)信號自適應(yīng)分解方法，相比于 EEMD、LMD［9］和 EWT［10］等傳統(tǒng)信號分解方法，VMD不存在模態(tài)混疊且能分解出頻率相近成分。然而，VMD算法有兩個不足之處：①模態(tài)分解結(jié)果受到二次懲罰參數(shù)α和模態(tài)個數(shù)K影響，難以選擇合適的數(shù)值；②當(dāng)噪聲方差超過一定閾值時，VMD算法可能會失效。為此，本文設(shè)計一種改進(jìn)型VMD信號分解方法，采用粒子群優(yōu)化（particle swarm optimization，PSO）［11］算法進(jìn)行全局優(yōu)化,確定參數(shù)α和K，并在VMD信號分解前基于非凸重疊組收縮（non?convex overlapping group shrinkage，NCOGS）［12］算法進(jìn)行信號降噪處理，提高信噪比，提取有效的故障特征信息。

針對噪聲情況下大轉(zhuǎn)速變化行星齒輪箱故障診斷問題，本文提出一種基于NLSTFT階次跟蹤和帶前置降噪自適應(yīng)變分模態(tài)分解的行星齒輪箱故障診斷方法。通過NLSTFT無鍵相階次跟蹤消除振動信號的變轉(zhuǎn)速頻率模糊現(xiàn)象，利用NCOGS算法降低信號噪聲，進(jìn)行VMD自適應(yīng)模態(tài)分解，提取各模態(tài)分量信號包絡(luò)譜特征頻率進(jìn)行故障診斷。

1 振動信號瞬時頻率估計與階次跟蹤

NLSTFT是一種新型的非平穩(wěn)信號時頻分析方法，它在STFT基礎(chǔ)上額外引入一個時變解調(diào)算子，用于解決調(diào)制成分帶來的時頻聚集性不高、時頻能量發(fā)散問題，可用于信號瞬時頻率估計。

根據(jù)Ville對瞬時頻率的定義：瞬時頻率等于信號瞬時相位的導(dǎo)數(shù)?？紤]一輸入信號x（t），基于Hilbert變換構(gòu)造的解析信號可以表示為

式中，s（t）為x（t）的解析信號；A（t）為幅值；ω（t）為解析信號的瞬時頻率。

由于信號中有較強(qiáng)的調(diào)制成分e-iωu，瞬時頻率的STFT幅值會小于諧波幅值，導(dǎo)致時頻譜模糊，且瞬時頻率在時頻譜上會出現(xiàn)時頻擴(kuò)散現(xiàn)象，時頻分辨率不高。為了消除調(diào)制成分帶來的負(fù)面影響，NLSTFT算法在STFT基礎(chǔ)上，乘以一個額外的時變解調(diào)算子e-ic(t)(u-t)2/2，具體表達(dá)式為

標(biāo)準(zhǔn)STFT的表達(dá)式為

其中，c（t）是瞬時頻率的一階導(dǎo)數(shù)。

如果解調(diào)算子和調(diào)制成分一致，則可以徹底消除調(diào)制影響，瞬時頻率幅值能夠達(dá)到最大化，使得時頻譜上的時頻聚集度及時頻分辨率高。

實際情況中，一般不知道瞬時頻率變化規(guī)律，c（t）的值也未知。本算法采用遞歸策略，在第一次計算 NLSTFT 時 c（t）=0，即 NLSTFT 退化為STFT算法。進(jìn)行信號時頻譜分析后，利用頻域峰值搜索算法得到對應(yīng)時刻t的峰值數(shù)據(jù)P(t)，并通過4階多項式擬合P(t)曲線得到瞬時頻率估計曲線I（t），然后把I（t）的一階導(dǎo)數(shù)作為第二次迭代時c（t）的值。由此，依次進(jìn)行迭代運(yùn)算，直到瞬時頻率曲線沒有明顯改變?yōu)橹埂Ｒ缘趎次與第n-1次迭代獲得的瞬時頻率曲線I（t）之間的平均誤差ξ作為迭代終止條件：

其中，In（t）、In-1（t）分別是第n次和第n-1次迭代計算得到的瞬時頻率估計值，一般選取δ=0.05。Lt為In（t）曲線時間長度，例如Lt=1 s。

利用NLSTFT得到估計瞬時頻率曲線后，對瞬時頻率求積分可得瞬時相位曲線，然后利用瞬時相位信息進(jìn)行振動信號重采樣，將與轉(zhuǎn)速有關(guān)的非平穩(wěn)時域信號轉(zhuǎn)化為與轉(zhuǎn)速無關(guān)的平穩(wěn)角域信號，完成無鍵相階次跟蹤過程。階次跟蹤算法流程如圖1所示。

NLSTFT算法具有復(fù)雜度低、運(yùn)算速度快、精度高等優(yōu)點，NLSTFT瞬時頻率估計精度高和抗噪能力好，能有效提高階次跟蹤精度和效率。

2 變分模態(tài)信號分解與故障診斷

VMD將信號分解過程轉(zhuǎn)化為通過迭代搜尋變分模型最優(yōu)解，使得每個模態(tài)估計帶寬之和最小。VMD算法對本征模態(tài)函數(shù)（intrinsic mode function，IMF）進(jìn)行了重新定義，它是一個調(diào)幅-調(diào)頻信號：

圖1 NLSTFT無鍵相階次跟蹤算法流程圖Fig.1 Flow chart of NLSTFT non-bonding phase order tracking algorithm

其中，Ak(t)為瞬時幅值；φk(t)為非單調(diào)遞減的瞬時相位；k表示第k個模態(tài)，k=1,2,…,K。

通過 φk（t）對 t求導(dǎo)，可得瞬時頻率 ωk(t)=φ˙k(t)。Ak(t)和ωk(t)相對于φk(t)來說變化非常緩慢，即在一段較長的時間間隔[t-δ,t+δ ]內(nèi),δ ≈ 2π φ˙k(t)，uk(t)可看作是幅值 Ak(t)、相位φk(t)的諧波信號。VMD自適應(yīng)模態(tài)分解過程是變分問題求解過程，包括變分問題構(gòu)造和變分問題求解兩部分。

2.1 變分問題構(gòu)造

變分問題可以描述為尋求K個模態(tài)函數(shù)uk(t)，使得每個模態(tài)的估計帶寬之和最小，約束條件為各模態(tài)之和等于輸入信號x（t），其中，帶寬的定義為

其中，Δ f為瞬時頻率的最大偏差；fFM為瞬時頻率的偏移率；fAM為Ak(t)的最高頻率。具體構(gòu)造過程如下：

（1）計算模態(tài)分量uk(t)的解析信號，通過Hil?bert變換得到單邊頻譜；

（2）通過指數(shù)修正，將一個預(yù)估中心頻率混合到解析信號中，使得每個模態(tài)分量的單邊頻率調(diào)制到各自估算的中心頻率處；

（3）估計各模態(tài)信號帶寬，使其帶寬之和最小:

其中，{uk}:={u1,u2,…,uk}為各模態(tài)函數(shù)集合，{ωk}:={ω1,ω2,…,ωk}為各中心頻率集合，?t表示函數(shù)對t求偏導(dǎo)。

2.2 變分問題求解

式（7）所示的優(yōu)化問題是帶約束的變分問題，引入二次懲罰因子α和拉格朗日乘積算子λ，將約束變分問題轉(zhuǎn)化為無約束變分問題，修改后的表達(dá)式為

其中，f(t)為原始信號的頻域表示。

利用交替方向乘積算子（alternate direction method of multipliers，ADMM）可以求解上述變分問題。VMD算法實現(xiàn)步驟如下：

（2）更新所有模態(tài)信號，即

（3）更新所有模態(tài)信號的中心頻率，即

（4）更新λ

其中，τ為噪聲容限參數(shù)，可以設(shè)定為0。

（5）重復(fù)步驟（2）～（4），直到滿足收斂條件為止。收斂條件：

2.3 VMD算法改進(jìn)

VMD算法分解結(jié)果受到二次懲罰參數(shù)α和模態(tài)個數(shù)K影響。α越小，計算得到的各個IMF分量的帶寬越大，則各個IMF分量中噪聲就越多，且各IMF成分間可能出現(xiàn)交叉。α越大，各個IMF分量的帶寬越小，可能會出現(xiàn)中心頻率偏移，導(dǎo)致模態(tài)分解錯誤。模態(tài)個數(shù)K決定中心頻率個數(shù)，K選取過小，模態(tài)分量間會出現(xiàn)混疊；K選取過大，則會導(dǎo)致某個模態(tài)分解分布在多個模態(tài)中，增加了計算時間，且會產(chǎn)生虛假分量。

VMD算法在低噪聲情況下信號分解魯棒性較好［13］，例如，當(dāng)噪聲方差為0.005時，可以達(dá)到100%的分解成功率，而當(dāng)噪聲方差大于0.1時，分解成功率降為0，算法就不收斂，無法正確分解信號。

針對VMD算法的上述不足，本文采用PSO算法全局尋優(yōu)確定K和α，且在VMD分解前，采用NCOGS算法進(jìn)行信號降噪，減小噪聲方差，解決VMD失效問題。

2.3.1 VMD參數(shù)優(yōu)選

當(dāng)行星齒輪箱出現(xiàn)故障時，振動信號會呈現(xiàn)沖擊調(diào)制特性。若某個IMF分量中包含故障敏感信息，其波形中會有相應(yīng)的沖擊脈沖，此時IMF分量包絡(luò)熵值較小，則將包絡(luò)譜熵作為PSO算法的適應(yīng)度函數(shù)，實現(xiàn)全局尋優(yōu)參數(shù)K和α。

2.3.2 NCOGS降噪算法

NCOGS是一種利用信號組稀疏特性的平移不變收縮算法，利用信號成組稀疏性將降噪問題轉(zhuǎn)化為凸優(yōu)化問題，算法計算效率高。

假定信號表達(dá)式為

其中，y（i）是含噪檢測信號，i∈ I={0，1，…，N-1}；N為信號點數(shù)；x(i)是不含噪聲的原始信號且已知具有組特性（即大幅值系數(shù)以組形式存在）；η為高斯加性噪聲。

將上述信號降噪問題轉(zhuǎn)化成最優(yōu)化問題：

定義罰函數(shù)為

其中，j={0,1,…,J-1},J為組個數(shù)。

代價函數(shù)轉(zhuǎn)化為

式（16）中，目標(biāo)函數(shù)、罰函數(shù)都是嚴(yán)格凸函數(shù)，能利用受控極小化算法進(jìn)行求解。

2.4 故障診斷算法

設(shè)計一種基于NLSTFT無鍵相階次跟蹤和帶前置降噪VMD分解的行星齒輪箱故障診斷方法，算法實現(xiàn)過程如下：

（1）采用NLSTFT精確估計瞬時頻率，對瞬時頻率曲線進(jìn)行積分，轉(zhuǎn)化為瞬時相位曲線，利用瞬時相位信息進(jìn)行角域重采樣完成階次跟蹤過程，輸出角域信號；

（2）利用NCOGS算法進(jìn)行角域信號降噪；

（3）用PSO算法優(yōu)化確定模態(tài)個數(shù)K和二次懲罰參數(shù)α。對降噪的角域信號進(jìn)行VMD自適應(yīng)模態(tài)分解，獲得K個模態(tài)分量信號；

（4）依據(jù)包絡(luò)譜熵選取對故障敏感的IMF分量，并對其進(jìn)行包絡(luò)譜解調(diào)分析，獲得故障特征頻率。進(jìn)而分析與正常狀態(tài)信號包絡(luò)譜特征頻率的差別，實現(xiàn)故障診斷。

3 仿真結(jié)果分析

3.1 仿真信號

當(dāng)行星齒輪箱發(fā)生局部損傷故障時，振動信號中有以嚙合頻率為載波頻率、齒輪轉(zhuǎn)頻及其倍頻為調(diào)制頻率的穩(wěn)態(tài)調(diào)頻信號，同時還存在周期性沖擊調(diào)幅信號［14］。

根據(jù)振動信號特點，構(gòu)造以下仿真信號：

表1 行星齒輪箱參數(shù)中齒輪的齒數(shù)Tab.1 Number of gear teeth in planetary gearbox

圖2所示分別為仿真信號波形、頻譜與時頻譜，由圖2b可知，轉(zhuǎn)速變化導(dǎo)致頻譜圖中頻率模糊現(xiàn)象，而且嚙合頻率及其邊頻出現(xiàn)混淆，容易誤診。

圖2 仿真信號波形及時頻譜Fig.2 The waveform and time-frequency spectrum of simulated signal

3.2 階次跟蹤

圖3 為NLSTFT對仿真信號第一次迭代運(yùn)算得到的瞬時頻率曲線，因第一次迭代時缺少先驗知識，故設(shè)置解調(diào)算子c（t）=0，該算法退化為ST?FT，此時瞬時頻率估計有較大偏差。

圖3 第一次迭代得到的瞬時頻率曲線Fig.3 IF curve obtained in the first iteration

由圖4可知，由于第二次迭代使用了第一次迭代得到的瞬時頻率先驗信息，故能大幅改善瞬時頻率估計精度。

圖4 第二次迭代得到的瞬時頻率曲線Fig.4 IF curve obtained in the second iteration

迭代終止發(fā)生在第4次迭代，如圖5所示。經(jīng)過4次迭代后得到了精確的瞬時頻率曲線，且時頻分辨率大幅提高。

圖5 第四次迭代Fig.5 The fourth iteration

將本方法與階次跟蹤中應(yīng)用廣泛的瞬時頻率估計算法CPP進(jìn)行對比。CPP參數(shù)設(shè)置為：動態(tài)時間支撐區(qū)設(shè)置為5，在路徑連接后，使用4階多項式對各區(qū)間曲線擬合。由圖6可知，CPP算法可以得到瞬時頻率曲線，但存在一定的瞬時頻率估計誤差，即使增加動態(tài)時間支撐區(qū)長度，效果改善程度也有限。

為驗證本文所提階次跟蹤方法的優(yōu)越性，選擇CPP和STFT算法進(jìn)行對比分析。由圖7可知：NLSTFT算法的相對誤差小于0.5%，明顯優(yōu)于CPP和STFT算法的瞬時頻率估計精度。

圖6 CPP估計瞬時頻率曲線Fig.6 IF curve estimated by CPP

圖7 三種方法對比Fig.7 Comparison of three IF estimation methods in relative error

表2所示為三種算法運(yùn)行時間的結(jié)果，其中仿真運(yùn)行環(huán)境為：MATLAB 2013a，CPU為i7?6700的8核處理器，主頻為3.4 GHZ，內(nèi)存為16 GB?？梢钥闯?，本文算法運(yùn)行時間遠(yuǎn)低于CPP算法。

表2 三種算法運(yùn)行時間Tab.2 Comparison of three methods in running time

利用NLSTFT算法得到瞬時頻率曲線后，通過積分計算得到瞬時相位，再通過重采樣轉(zhuǎn)換為角域信號，完成階次跟蹤過程。圖8為階次跟蹤后的譜圖，可發(fā)現(xiàn)頻率模糊現(xiàn)象和時頻聚集性明顯改善。

圖8 階次跟蹤后信號譜圖Fig.8 Signal spectra after order tracking

進(jìn)一步評估本文算法的噪聲魯棒性，分析不同信噪比條件下的計算測量值和真實值的均方根誤差，如圖9所示，可以看出：在一定的噪聲范圍內(nèi)，本文方法的均方根誤差均小于CPP和STFT算法，抗噪聲能力較強(qiáng)。

圖9 噪聲對瞬時頻率估計的影響Fig.9 The effect of noise on the IF estimation

3.3 信號分解

經(jīng)過階次跟蹤消除轉(zhuǎn)速波動影響后，采用VMD自適應(yīng)模態(tài)分解法分離故障信號，利用PSO算法確定VMD的參數(shù)K和α，為了避免在噪聲情況下VMD失效，采用NCOGS算法進(jìn)行濾波去噪。

圖10所示為對角域加噪信號進(jìn)行NCOGS濾波處理結(jié)果，NCOGS能有效去除加性噪聲，提高信噪比。

圖10 NCOGS降噪前后角域圖Fig.10 Angular domain waveform before and afterNCOGS de-noising method

對NCOGS降噪后的角域信號進(jìn)行VMD分解，由PSO算法得到參數(shù)K=3，α=2 050。VMD將原始信號分解為3個模態(tài)分量，如圖11、圖12所示。模態(tài)分量2包含故障信息，對其進(jìn)行包絡(luò)譜分析，得到故障特征頻率，如圖13所示。其中，太陽輪局部故障特征頻率fs=3.33階幅值遠(yuǎn)大于其他階次幅值，可正確識別太陽輪局部故障，驗證了本文算法有效性。

圖11 VMD分解得到三個模態(tài)分量Fig.11 Three IMFs obtained by VMD

圖12 各模態(tài)分量Fig.12 Each component of IMFs

圖13 模態(tài)分量2包絡(luò)譜Fig.13 Envelop spectrum of the second IMF

4 實驗結(jié)果分析

4.1 實驗平臺與模擬故障

圖14所示為行星齒輪箱故障診斷綜合實驗臺，由驅(qū)動電機(jī)、可編程控制面板、傳動齒輪箱、磁粉制動器、加速度傳感器、轉(zhuǎn)速傳感器、數(shù)據(jù)采集器等組成。行星齒輪箱參數(shù)見表1。實驗中模擬故障為太陽輪切齒故障，故障位置在第二級太陽輪，如圖15所示。第2級行星齒輪箱輸出軸承所受的負(fù)荷為40 N·m，轉(zhuǎn)速變化范圍是2400～3600 r/min。加速度傳感器安裝在行星齒輪箱箱體上，振動信號采樣頻率為10 kHz，采樣時間為3 min。

圖14 行星齒輪箱故障診斷實驗平臺Fig.14 Experiment platform for fault diagnosis of planetary gearboxes

根據(jù)行星齒輪箱參數(shù)，計算出相應(yīng)的嚙合頻率和故障特征頻率，如表3所示。其中，驅(qū)動電機(jī)旋轉(zhuǎn)頻率

圖15 太陽輪切齒故障Fig.15 Chipped tooth of sun gear

表3 行星齒輪箱有關(guān)特征頻率Tab.3 Characteristic frequency of planetary gearbox

4.2 振動信號分解與故障診斷結(jié)果分析

振動信號時域圖和頻譜圖見圖16。由圖16b可知，因二級嚙合頻率過低且離傳感器較遠(yuǎn)，導(dǎo)致被一級嚙合頻率完全被覆蓋。同時，轉(zhuǎn)速波動引起Fourier頻譜出現(xiàn)左右漂移現(xiàn)象，難以準(zhǔn)確分析時頻變化規(guī)律，故障容易誤診。

圖16 振動信號時域圖和頻譜圖Fig.16 Vibration signals under the sun gearfault condition

圖17 為NLSTFT階次跟蹤估計得到的轉(zhuǎn)速曲線。圖18所示為本文算法階次跟蹤后的階次譜，可看出，較好地消除了頻率漂移現(xiàn)象，頻率成分以fm1及其諧波、f(r)s1及其倍頻成分為主，而二階嚙合頻率被掩蓋。

圖17 NLSTFT估計轉(zhuǎn)速信號Fig.17 Speed signal estimated by NLSTFT

圖18 NLSTFT獲得的階次譜Fig.18 Order spectrum obtained by NLSTFT algorithm

將NLSTFT階次跟蹤法與STFT和CPP法進(jìn)行比較。由表4可知，NLSTFT算法的瞬時頻率估計相對誤差小于0.6%，優(yōu)于STFT算法，且運(yùn)行時間與STFT算法接近，而CPP算法運(yùn)行時間相對較長。證明了NLSTFT算法的瞬時頻率估計精度高、計算復(fù)雜度較低。

表4 三種瞬時頻率估計算法比較Tab.4 Comparison of three IF estimation methods

因第二級行星齒輪系和特征頻率集中在低頻部分，在故障診斷前，需要去除一級行星齒輪系的影響，故對角域信號進(jìn)行低通濾波，只分析［0，10］階頻段內(nèi)信號。圖19a所示為濾波后角域信號由VMD分解獲得的3個模態(tài)分量，PSO算法得到優(yōu)化參數(shù)K=3，α=2 075。為提取故障特征，對3個模態(tài)分量均進(jìn)行包絡(luò)譜分析，發(fā)現(xiàn)模態(tài)分量2包含了故障信息，如圖19b所示。由表3可知，太陽輪局部故障特征頻率fs2=0.44 Hz。由圖19b可知，故障特征頻率fs2階處幅值明顯大于其他頻率成分，很容易判斷出現(xiàn)了太陽輪局部故障，診斷結(jié)果與實際情況相符合。

5 結(jié)論

（1）給出了一種NLSTFT無鍵相階次跟蹤算法，新算法的時頻聚集性好、計算復(fù)雜度低，適用于轉(zhuǎn)速變化大、噪聲強(qiáng)情況下振動信號瞬時頻率估計和無鍵相階次跟蹤。

圖19 VMD分解模態(tài)分量Fig.19 IMFs obtained by VMD algorithm

（2）設(shè)計了基于PSO的VMD模態(tài)分解算法，并采用非凸重疊組收縮算法進(jìn)行角域信號濾波降噪，避免了噪聲情況下VMD模態(tài)分解失效。

（3）提出了一種基于NLSTFT無鍵相階次跟蹤和粒子群優(yōu)化VMD模態(tài)分解的行星齒輪箱故障診斷方法，提高了變轉(zhuǎn)速與噪聲情況下行星齒輪箱故障診斷性能。