基于GLCT和FPO的時變工況行星齒輪箱故障特征提取*

孫斌 李宏坤 劉艾強 張孔亮

（大連理工大學機械工程學院）

0 引言

行星齒輪箱在工業應用廣泛，如風力發電，航天航空等[1]。行星齒輪箱作為旋轉機械的主要傳動部件，長期工作在沖擊、交變載荷和工況大幅度變化的情況下，其關鍵零部件如齒輪，極易發生故障，從而引起巨大的經濟損失和惡性事故[2-3]。因此，對行星齒輪箱的齒輪故障診斷具有重大的經濟與安全意義。

行星齒輪箱在變轉速工況下，所測得的振動信號具有典型的非平穩和非線性特點。如果利用傳統的信號處理方法，如快速傅里葉變化，得到的頻譜模糊，無法獲得故障特征頻率，難以識別變轉速下的齒輪故障[4]。

時變轉速下對行星齒輪箱的故障診斷，階次跟蹤被認為是最有效的工具之一[5]。在階次跟蹤中，原始信號以恒定的角度間隔重新被采樣，原始的時域非平穩信號轉化為角域的平穩信號。階次跟蹤的方法可以消除由于轉速的變化帶來的頻譜模糊的現象，在傳統的階次跟蹤中，需要安裝硬件設備獲得轉速信息用于信號的重采樣。但是，由于空間和成本有限，安裝硬件設備在很多場合都難以實現[6]。

為了克服上述問題，從振動信號的時頻表示中提取到瞬時速度信息是比較常用和節約成本的方法。但是，時頻聚集性對所要提取的瞬時速度精度影響較大。許多學者將短時傅里葉變換（STFT）和連續小波變換（CWT）的傳統時頻分析方法，應用于機械故障診斷領域。雖然這些傳統的時頻方法能夠處理非平穩信號，但是受海森堡（Heisenberg）不確定性原理的限制，時間和頻率分辨率不能同時達到最佳，同時也對噪聲比較敏感[7]。此外，由于所測得行星齒輪箱的信號是典型的調制信號，使用傳統的時頻分析會出現能量混疊和時頻脊線模糊，導致難以準確提取到轉頻信息[8]。為了提高時頻聚集性，I Daubechies 等提出了同步壓縮變換（SST）,利用時頻重排的方法，獲得了聚集性良好的時頻表示[9]。Gang Yu等人提出了同步提取變換(SET)，與SST的壓縮方式不同，只保留與時變特征最相關的時頻信息，大大提高了時頻的能量聚集性[10]。但是上述時頻分析方法存在對噪聲敏感，不能很好的表征非線性特征的多分量信號的問題。本文運用一種新型的時頻分析方法，廣義線性調頻小波變換（GLCT）,可以克服現有時頻方法的一些缺陷，并且對噪聲不敏感[11]。

從時頻表示中，提取到轉頻曲線是無轉速階次分析中不可或缺的一部分。Zhao等使用峰值搜索算法從時頻圖中提取瞬時頻率，實現了變轉速下軸承的故障診斷[12]。陳劍等人提出一種改進的Crazy Cliber算法，從STFT 時頻脊線提取脊線，并成功實現了變轉速下的軸承故障特征提取[13]。但是，峰值搜索算法和Crazy Cliber算法容易陷入局部最優，考慮不到整個信號的輪廓，如果選擇不合適的起點，提取到的瞬時頻率將會是錯誤的，而且這種方法在低信噪比中是無效的[14]。

針對上述存在時頻聚集性不高和提取到的瞬時頻率準確性低的問題，本文提出了一種基于GLCT和FPO相結合的方法，用于診斷在時變轉速下行星齒輪箱齒輪的故障。采用GLCT對振動信號進行時頻表示，該時頻方法比傳統的時頻分析方法如STFT，時頻聚集性高，并且對噪聲不敏感；再采用FPO 從時頻矩陣中提取瞬時速度信息，該提取算法克服了峰值搜索算法的噪聲干擾和局部最優[15]，然后利用所提取到的瞬時速度信息進行信號重采樣，得到角域平穩信號，再進行包絡階次分析，得到故障特征階次，達到故障診斷的目的。仿真和實驗的振動信號證明了該方法在變轉速運行下行星齒輪箱齒輪故障診斷的有效性。

1 基本原理介紹

1.1 廣義線性調頻小波變換(GLCT)

標準的短時傅里葉變換（STFT）如式(1)：

將引進一個解調因子，去消除調制成分的影響?？紤]到調制成分是隨著時間變化的，所以引入的解調因子也必須是時變的。引入時變的解調因子s(t) 的STFT表達式為：

s(t)的STFT幅值用瞬時頻率φ(t')可以表達為：

從上式可以看出，如果解調因子和分析信號的調制成分相一致，則瞬時頻率φ(t')的STFT 幅值達到最大，然后時頻能量在瞬時頻率處具有很高的聚集性。但是在大多數情況下，一個信號的IF 特征不能預先得知，則解調因子e-ic(t')(u-t')2/2難以準確確定，特別是含有復雜的多成分信號中，更加困難。為了解決上述問題，用離散的解調因子無限去逼近最優的解調因子e-ic(t')(u-t')2/2，則用離散化的解調因子表示STFT 的公式如式(4)：

如果離散的解調因子與信號的調制成分相近時，則在時頻表示中，瞬時頻率的時頻點(t',ω)有很高的能量聚集，則|S(t',ω,c) |幅值達到最大值。從而能夠得到最佳的參數c為：

則時頻表達式為：

頻譜的定義式為：

為了得到最優的解調因子e-ic(t')(u-t')2/2，引入參數α，通過旋轉arctan(-c)，對時頻結果產生影響，則時頻平面旋轉為：

則上述的時頻表達為：

因為參數α有N個值，則時頻平面能均分N+1 個部分：

由上述可知，GLCT比STFT多引入了一個參數N,如果N=1時，GLCT即為STFT。

1.2 快速路徑最優算法(FPO)

經過上述GLCT 得到的時頻結果為GS(t',ω)，計算時頻表示中每個時刻幅值的局部極大值，然后得到該時刻所有極大值所對應的頻率值為：

其中，tn(n=1,2,...,N)表示時頻表示的時刻點；Np(tn)表示時刻tn的局部極大值的個數；fm(tn)是時刻tn處第m個極值所對應的頻率值，其幅值的大小記為SPm(tn)。

采用FPO 確定每個時刻點處應該被提取出來的幅值極大值，極大值所對應的頻率連接起來即為目標脊線，具體算法如下：

其中，mc(tn)為時刻tn應該被提取的幅值極大值；[F]表示優化方程，該方程的解就是目標脊線。優化方程為：

其中，Δfd是fd的導數；percp[f(t)]表示f(t) 的第p分位數。則優化方程的解為：

其中，n=1,2,...,N,m=1,2,...,Np(tn)，q(m,tn)確定每個時刻幅值的極大值與前一時刻哪一個極值點構成優化路徑，規劃所有路徑的可能性。U(m,tn)表示優化結果的中間向量，選擇該時刻U(m,tn)最大值的幅值極大值點對應的頻率值作為目標脊線。

2 提取方法流程

為解決時變工況下無轉速信號和傳統時頻方法時頻聚集性不高導致提取的瞬時轉速誤差較大的問題，提出了基于GLCT 和FPO 的時變轉速行星齒輪箱故障診斷方法。所提方法的流程圖如圖1，具體實施步驟如下：

圖1 所提方法流程圖Fig.1 The flowchart of proposed method

步驟1：采集時變轉速下行星齒輪箱的振動信號。

步驟2：對獲取的振動信號進行GLCT 時頻表示GS(t',ω)。

步驟3：設置最大搜索頻率，在GS(t',ω)時頻域中采用FPO算法提取瞬時轉頻趨勢線。

步驟4：采用最小二乘法對提取的趨勢線進行擬合，得到目標脊線；進而利用提取的轉速曲線對原始振動信號進行等角度重采樣，將非平穩信號轉換為角域平穩信號。

步驟5：對角域信號進行包絡階次分析，從包絡階次譜中提取故障特征階次并判斷故障類型。

3 仿真驗證

為了證明上述方法的有效性，仿真一個如下時變轉速下的太陽輪故障信號：

其中，xs(t)為軸的振動信號；xm(t)為太陽輪局部故障振動信號；n(t)為高斯白噪聲，具體數學模型如式(17)和(18)：

式中，(t)為軸旋轉頻率；fs(t)為太陽輪故障特征頻率；fm(t)為齒輪嚙合頻率。

考慮到時變轉速，則上述表達式為：

在這個仿真信號中，采樣頻率Fs=1024Hz,采樣時間為10 秒，n(t)=-5dB。而仿真信號中，太陽輪故障階次Os為：

上述仿真信號的時域波形及其FFT 頻譜和包絡譜如圖2所示。由圖2可知，故障沖擊被完全淹沒在噪聲中，此外受時變轉速影響，故障特征頻率以及調制成分在頻譜圖中出現混疊和模糊現象，從而不能提取到有用的故障特征信息。

圖2 仿真信號時域圖及其相關譜圖Fig.2 Time-domain diagram of simulated signal and its related spectrum diagram

圖3 為本文所提方法的分析結果，圖3(a)為仿真信號GLCT時頻表示，瞬時頻率最低的時頻脊線為理論轉頻曲線。因此，以瞬時頻率最低的時頻脊線為基頻O即為所提取的目標脊線，進行角度重采樣。用FPO從圖3(a)的時頻表示中提取到的目標脊線如圖3(b)的黑色曲線所示，因所提取的目標脊線不夠光滑，本文采用最小二乘法進行曲線擬合，圖3(b)的藍色線為擬合曲線，紅色線為理論曲線。從圖3(b)可以明顯看出擬合后的曲線和理論曲線基本重合。根據上述提取出來的曲線進行角度重采樣包絡階次譜，從圖中可以明顯得到O,Os-O,Os,Os+O,2Os。從上述的分析結果可以得出本文所提出的方法能夠準確提取到變轉速下仿真信號的太陽輪故障特征。

圖3 本文所提方法分析結果Fig.3 The analysis result of the method proposed in this article

為了驗證本文所提方法的優越性，分別采用STFT與FPO 和GLCT 與峰值搜索算法進行對比。仿真信號的STFT時頻表示如圖4(a)所示，用FPO從圖4(a)中提取到的目標脊線如圖4(b)黑色曲線所示，紅色線為最小二乘擬合曲線，藍色線為理論曲線，使用擬合后的曲線得到的包絡階次譜如圖4(c)所示?？梢钥闯?，由于所采用的時頻分析方法時頻聚集性不高，并且對噪聲比較敏感，導致提取出來的目標脊線與理論曲線誤差較大，進而最后的診斷效果不理想。

圖4 STFT與FPO分析結果Fig.4 STFT and FPO analysis results

從圖5可以看出，由于峰值搜索算法考慮不到時頻矩陣的全局信息，容易出現局部最優的現象，導致所提取到的目標曲線與理論曲線誤差較大，最后的診斷效果不理想。

圖5 GLCT與峰值搜索算法結果包絡階次譜Fig.5 GLCT and peak search analysis results envelope order spectrum

綜上，本文所提出的方法，能克服傳統時頻分析方法的時頻聚集性不高，而帶來的目標脊線提取精度不高以及傳統的時頻脊線提取算法由于考慮不到時頻矩陣的全局信息，導致提取到的目標脊線局部最優的問題，最后能夠準確提取時變轉速下太陽輪故障特征，從而驗證了本文所提的方法在時變轉速下對行星齒輪箱齒輪故障特征提取的有效性。

4 實驗驗證

4.1 實驗說明

設計了行星齒輪箱太陽輪和行星輪故障在時變轉速下的模擬實驗，實驗系統如圖6 所示，采集實驗是在圖6所示的實驗臺上進行的，主要由變頻器、驅動電機、行星齒輪箱、磁粉制動器、轉速傳感器、加速度傳感器（型號為美國DYTRAN 公司生產的3035B 型傳感器）和安裝NI9234采集卡的工控機組成。

圖6 時變轉速下行星齒輪箱實驗系統Fig.6 Planetary gearbox experimental system under time-varying speed

在本實驗中，行星齒輪箱的齒圈固定，太陽輪為輸入端，行星架為輸出端，電機的轉速在20s內從0～1800r/min線性增加，為了驗證提取到的目標脊線的準確度，以轉速傳感器測得的瞬時轉速為理論轉速。采樣頻率為12800Hz，為了增加分析的可靠性，本文取5～15秒的數據。磁粉制動器施加的扭矩為0，行星齒輪箱的參數見表1。為了模擬行星齒輪箱齒輪的局部故障，通過線切割技術在太陽輪的某個輪齒上切割一部分和切掉行星輪一個輪齒，這樣形成太陽輪斷齒和行星輪缺齒局部故障，故障齒輪如圖7所示。行星齒輪箱的故障特征階次見表2。

表1 行星齒輪箱的參數Tab.1 Parameters of planetary gearbox

圖7 齒輪故障件Fig.7 Gear failure parts

表2 行星齒輪箱各部分故障特征階次Tab.2 Characteristic order of failure of each part of planetary gearbox

4.2 實驗驗證1：太陽輪故障信號分析

太陽輪缺齒故障振動信號的時域波形及其對應的傅里葉頻譜和包絡譜如圖8所示?？梢钥吹?，隨著轉速的升高，沖擊逐漸變大；同時，受時變轉速影響，譜圖中會出現頻譜模糊現象，因此根據傳統的時頻分析很難提取到有價值的故障特征信息。

圖8 太陽輪斷齒時域圖及其相關譜圖Fig.8 Time domain diagram of sun gear broken tooth and its related spectrum diagram

為提升運算效率，本文對原始信號進行降采樣預處理，降采樣的頻率設置為512Hz。圖9(a)為降采樣信號GLCT 時頻表示，從圖中可以明顯的看出來，瞬時頻率最低的時頻脊線為太陽輪輸入端的瞬時速度曲線，即目標提取脊線。用FPO 從圖9(a)提取到的目標瞬時頻率脊線如圖9(b)的黑色曲線所示，本文采用最小二乘法擬合所提取出來的目標曲線為圖9(b)的紅色曲線，藍色線為理論曲線，從圖中可以看出擬合后的曲線和理論轉頻曲線基本重合。根據上述所擬合后的曲線進行角度重采樣后的包絡階次譜如圖9(d)所示，可以清晰地觀察到O,Os-O,Os,Os+O,2Os,2O的故障特征信息。上述分析結果表明太陽輪出現故障，從而驗證了所提算法的有效性。

圖9 本文所提方法分析結果Fig.9 The analysis result of the method proposed in this article

為了驗證本文所提方法的優越性，分別采用STFT 與FPO 和GLCT 與峰值搜索算法進行對比。圖10(a)為太陽輪斷齒信號的STFT 時頻表示，用FPO 從STFT 時頻表示中提取到的目標脊線如圖10(b)的黑色曲線所示，紅色線為擬合曲線，藍色線為理論曲線，使用擬合曲線得到的包絡階次譜如圖10(c)所示。從最后的診斷結果可以看出，所采用的時頻分析方法時頻聚集性不高，并且對噪聲敏感，導致提取出來的目標曲線與理論的曲線誤差較大，進而最后的診斷效果不理想。

圖10 STFT與FPO分析結果Fig.10 STFT and FPO analysis results

GLCT 和峰值搜索的包絡階次譜如圖11 所示?？梢钥闯?，從最后診斷結果可以看出，因為峰值搜索算法一旦出現選錯頻率起始點，將導致所提取到的目標曲線錯誤，從而使該方法失效。

圖11 GLCT與峰值搜索算法結果包絡階次譜Fig.11 GLCT and peak search analysis results envelope order spectrum

4.3 實驗驗證2：行星輪故障信號分析

行星輪缺齒故障振動信號的時域波形及其對應的傅里葉頻譜和包絡譜如圖12 所示?？梢钥吹?，隨著轉速的升高，沖擊逐漸變大；同時，受時變轉速影響，譜圖中會出現頻譜模糊現象，因此根據傳統的時頻分析很難提取到有價值的故障特征信息。

圖12 行星輪缺齒時域圖及其相關譜圖Fig.12 Time domain diagram of planetary gear missing teeth and its related spectrum diagram

圖13(a)為降采樣信號GLCT時頻表示，從圖13(a)可以明顯的看出來，瞬時頻率最低的時頻脊線為目標提取脊線。用FPO從圖13(a)提取到的目標脊線如圖13(b)黑色曲線所示，紅色線為擬合后的曲線，藍色線為理論曲線，從圖中可以看出擬合后的曲線和理論轉頻曲線基本重合。根據上述所擬合后的曲線進行角度重采樣，最后的包絡階次譜如圖13(d)所示，可以明顯得到O,Op-Oc,Op,Op+Oc,2Op,2O,4Op的故障信息。上述分析結果表明行星齒輪出現故障，從而驗證了所提算法的有效性。

圖13 本文所提方法分析結果Fig.13 The analysis result of the method proposed in this article

為了驗證本文所提方法的優越性，分別采用STFT與FPO 和GLCT 與峰值搜索算法進行對比。行星輪缺齒信號的STFT時頻表示如圖14(a)所示，用FPO從圖14(a)中提取到的目標曲線如圖14(b)黑色曲線所示，紅色線為擬合曲線，藍色線為理論曲線,使用擬合曲線得到最后的包絡階次譜如圖14(c)所示。從最后的診斷結果可以看出，因為所采用的時頻分析方法時頻聚集性不高，并且對噪聲敏感，導致提取出來的時頻曲線與理論的時頻曲線誤差較大，導致最后的診斷效果不理想。

圖14 STFT與FPO分析結果Fig.14 STFT and FPO analysis results

GLCT和峰值搜索的包絡階次譜如圖15所示。從最后診斷結果可以看出，峰值搜索算法一旦選錯頻率起始點，將導致所提取到的目標曲線錯誤，從而使該方法失效。

圖15 GLCT與峰值搜索算法結果包絡階次譜Fig.15 GLCT and peak search analysis results envelope order spectrum

綜上，本文提出的方法，能在強噪聲的背景下和時變轉速工況下，準確提取到行星輪缺齒故障信息。提出的GLCT方法能夠生成時頻聚集性較高的時頻圖，克服了傳統的時頻方法在強噪聲背景下時頻聚集性不高的問題。提出的FPO脊線提取算法能夠考慮時頻矩陣的全局信息，準確提取到目標脊線，克服了傳統脊線提取算法容易陷入局部最優的問題。

5 結論

針對無轉速計時變工況下行星齒輪箱故障特征難以提取的問題，提出了基于GLCT與FPO的故障特征提取算法。通過構造仿真信號和實際應用對其有效性和優越性進行分析，可以得到如下結論：

1）在噪聲干擾比較嚴重的情況下，GLCT的時頻聚集性和噪聲魯棒性優于傳統的STFT 方法，從而可以清晰地展現齒輪箱瞬時頻率曲線，為后續目標脊線提取提供基礎。

2）FPO時頻脊線搜索算法能夠考慮到時頻矩陣的全局信息，提取到的時變轉速曲線與理論曲線能夠完全重合，克服了峰值搜索算法容易陷入局部最優和提取目標脊線不準確的問題。

3）該算法能在不安裝轉速計等硬件的條件下提取轉速曲線，隨后對原始振動信號進行角域重采樣，在其包絡階次譜中能夠準確識別到行星齒輪箱故障特征階次，為變轉速工況下的行星齒輪箱故障診斷提供了一種新思路。