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

王友仁 王 俊 黃海安

南京航空航天大學自動化學院，南京，211106

0 引言

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

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

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

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

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

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

根據Ville對瞬時頻率的定義：瞬時頻率等于信號瞬時相位的導數。考慮一輸入信號x（t），基于Hilbert變換構造的解析信號可以表示為

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

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

標準STFT的表達式為

其中，c（t）是瞬時頻率的一階導數。

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

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

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

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

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

2 變分模態信號分解與故障診斷

VMD將信號分解過程轉化為通過迭代搜尋變分模型最優解，使得每個模態估計帶寬之和最小。VMD算法對本征模態函數（intrinsic mode function，IMF）進行了重新定義，它是一個調幅-調頻信號：

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

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

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

2.1 變分問題構造

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

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

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

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

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

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

2.2 變分問題求解

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

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

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

（2）更新所有模態信號，即

（3）更新所有模態信號的中心頻率，即

（4）更新λ

其中，τ為噪聲容限參數，可以設定為0。

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

2.3 VMD算法改進

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

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

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

2.3.1 VMD參數優選

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

2.3.2 NCOGS降噪算法

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

假定信號表達式為

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

將上述信號降噪問題轉化成最優化問題：

定義罰函數為

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

代價函數轉化為

式（16）中，目標函數、罰函數都是嚴格凸函數，能利用受控極小化算法進行求解。

2.4 故障診斷算法

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

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

（2）利用NCOGS算法進行角域信號降噪；

（3）用PSO算法優化確定模態個數K和二次懲罰參數α。對降噪的角域信號進行VMD自適應模態分解，獲得K個模態分量信號；

（4）依據包絡譜熵選取對故障敏感的IMF分量，并對其進行包絡譜解調分析，獲得故障特征頻率。進而分析與正常狀態信號包絡譜特征頻率的差別，實現故障診斷。

3 仿真結果分析

3.1 仿真信號

當行星齒輪箱發生局部損傷故障時，振動信號中有以嚙合頻率為載波頻率、齒輪轉頻及其倍頻為調制頻率的穩態調頻信號，同時還存在周期性沖擊調幅信號［14］。

根據振動信號特點，構造以下仿真信號：

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

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

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

3.2 階次跟蹤

圖3 為NLSTFT對仿真信號第一次迭代運算得到的瞬時頻率曲線，因第一次迭代時缺少先驗知識，故設置解調算子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

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

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

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

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

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

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

表2所示為三種算法運行時間的結果，其中仿真運行環境為：MATLAB 2013a，CPU為i7?6700的8核處理器，主頻為3.4 GHZ，內存為16 GB。可以看出，本文算法運行時間遠低于CPP算法。

表2 三種算法運行時間Tab.2 Comparison of three methods in running time

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

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

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

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

3.3 信號分解

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

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

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

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

圖11 VMD分解得到三個模態分量Fig.11 Three IMFs obtained by VMD

圖12 各模態分量Fig.12 Each component of IMFs

圖13 模態分量2包絡譜Fig.13 Envelop spectrum of the second IMF

4 實驗結果分析

4.1 實驗平臺與模擬故障

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

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

根據行星齒輪箱參數，計算出相應的嚙合頻率和故障特征頻率，如表3所示。其中，驅動電機旋轉頻率

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

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

4.2 振動信號分解與故障診斷結果分析

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

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

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

圖17 NLSTFT估計轉速信號Fig.17 Speed signal estimated by NLSTFT

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

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

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

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

5 結論

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

圖19 VMD分解模態分量Fig.19 IMFs obtained by VMD algorithm

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

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