Vold-Kalman濾波和高階能量分離在時變工況行星齒輪箱故障診斷中的應(yīng)用研究

秦嗣峰，馮志鵬，LIANG Ming

（1.北京科技大學(xué)機(jī)械工程學(xué)院，北京100083；2.渥太華大學(xué)機(jī)械工程系，安大略渥太華K1N 6N5，加拿大）

Vold-Kalman濾波和高階能量分離在時變工況行星齒輪箱故障診斷中的應(yīng)用研究

秦嗣峰1，馮志鵬1，LIANG Ming2

（1.北京科技大學(xué)機(jī)械工程學(xué)院，北京100083；2.渥太華大學(xué)機(jī)械工程系，安大略渥太華K1N 6N5，加拿大）

時變工況下的行星齒輪箱振動信號具有明顯的時變調(diào)制特點(diǎn)，常規(guī)的頻譜分析方法難以識別齒輪故障特征頻率。提出了基于Vold-Kalman濾波和能量分離的時頻分析方法，識別行星齒輪箱的時變特征頻率，診斷齒輪故障。與傳統(tǒng)的時頻分析方法相比，基于Vold-Kalman濾波和能量分離的時頻表示具有良好的時頻分辨率，而且沒有交叉項干擾，能夠有效提取非平穩(wěn)信號中的時變頻率成分。分析了行星齒輪箱時變工況下的實(shí)驗信號，準(zhǔn)確地診斷了齒輪故障，驗證了該方法的有效性。

故障診斷；行星齒輪箱；高階能量分離；Vold-Kalman濾波；時頻分析

引 言

行星齒輪箱在風(fēng)力發(fā)電機(jī)、直升機(jī)和車輛等裝備中應(yīng)用廣泛，實(shí)際運(yùn)行中，其工況經(jīng)常發(fā)生變化，產(chǎn)生非平穩(wěn)的振動響應(yīng)。目前的研究主要集中在恒定工況方面［1-2］，而針對時變工況方面的研究并不多見，其中，Williams和Zalubas［3］考慮直升機(jī)傳動系統(tǒng)振動信號的非平穩(wěn)特點(diǎn)，應(yīng)用Wigner-Ville分布檢測行星齒輪箱故障。Bartelmus和Zimroz［4-5］發(fā)現(xiàn)行星齒輪箱故障對載荷變化敏感，應(yīng)用齒輪嚙合振動幅值與工況之間的線性依賴關(guān)系監(jiān)測時變工況下的行星齒輪箱運(yùn)行狀態(tài)。

行星齒輪箱故障診斷的關(guān)鍵是根據(jù)故障特征頻率及其幅值的變化識別故障原因。行星齒輪箱嚙合頻率和齒輪故障特征頻率與轉(zhuǎn)速直接相關(guān)。工況（轉(zhuǎn)速和負(fù)荷）的變化將引起轉(zhuǎn)速波動，嚙合頻率和齒輪故障特征頻率也會隨之發(fā)生變化。因此，如何從非平穩(wěn)信號中有效地識別齒輪故障特征頻率及其時變特征，是時變工況下行星齒輪箱故障診斷的關(guān)鍵問題。

對于時變工況下的行星齒輪箱的非平穩(wěn)振動信號，時頻分析是一種有效的分析手段。但是，常見的時頻分析方法各有優(yōu)缺點(diǎn)［6-7］。例如：線性時頻表示（如短時Fourier變換和小波變換）雖然不存在交叉項干擾問題，但是其時頻分辨率受Heisenberg不確定性原理的限制，不能同時達(dá)到最佳；Wigner-Ville分布雖然時頻分辨率高，卻存在固有的交叉項干擾問題，不適合分析復(fù)雜多分量信號；以Wigner-Ville分布為基礎(chǔ)的雙線性時頻分布（包括Cohen類分布和仿射類分布）通過各種核函數(shù)進(jìn)行平滑處理，雖然抑制了交叉項干擾，但是會造成信號自項畸變，降低時頻分辨率；Hilbert-Huang變換雖然對信號的時變具有自適應(yīng)性，而且具有良好的局部時頻聚集能力，但是對信號中的奇異點(diǎn)敏感，容易產(chǎn)生模式混淆，得到虛假的本質(zhì)模式函數(shù)，影響分析結(jié)果。

行星齒輪箱振動信號具有明顯的調(diào)制特點(diǎn)，在時變工況下，將表現(xiàn)出時變調(diào)制特征。非線性信號處理中的能量算子方法為分析這種時變調(diào)制信號提供了有效的途徑。能量算子通過信號的瞬時幅值及其微分的非線性組合估計信號的“能量”，在其基礎(chǔ)上提出的能量分離算法可以計算任意調(diào)制信號的幅值包絡(luò)和瞬時頻率［8-12］。能 量算 子 和 能 量分離在信號及其導(dǎo)數(shù)的基礎(chǔ)上進(jìn)行計算，算法簡單、時間分辨率高，對信號的瞬態(tài)變化具有良好的適應(yīng)性。但是，該方法僅適用于單分量信號。對于成分復(fù)雜的行星齒輪箱振動信號，需要將其分解為單分量再進(jìn)行分析。

近年來提出的Vold-Kalman濾波方法［13-15］在最小化結(jié)構(gòu)和數(shù)據(jù)方程誤差的基礎(chǔ)上，可以有效分解復(fù)雜多分量信號的頻率成分，為解決頻率交錯重疊的行星齒輪箱振動信號單分量分解問題提供了一種有效手段。

針對時變工況下行星齒輪箱故障診斷面臨的復(fù)雜非平穩(wěn)振動信號分析問題，考慮Vold-Kalman濾波方法在單分量分解以及能量分離方法在瞬時頻率計算方面的獨(dú)特優(yōu)勢，提出了基于二者的時頻分析方法。

1 Vold-Kalman濾波

Vold-Kalman濾波可以將復(fù)雜多分量信號分解為單分量。與傳統(tǒng)的濾波方法相比，該方法直接從時域中提取感興趣的信號分量，避免了由時域至頻域變換帶來的相位偏差；與經(jīng)驗?zāi)Ｊ椒纸夥椒ㄏ啾龋琕old-Kalman濾波的中心頻率可以根據(jù)瞬時頻率自適應(yīng)調(diào)節(jié)，能夠有效分離在時頻域內(nèi)鄰近甚至交叉的信號分量［14-15］，避免 了 由經(jīng)驗 模 式分 解 方法帶 來的模式混淆現(xiàn)象。

任意調(diào)制信號可以表示為

式中k為階次，Ak（t）為第k個分量的幅值包絡(luò)，Θk（t）為載波信號。

式中ω（τ）為瞬時角頻率（τ）dτ為瞬時相位。

與載波信號相比，幅值包絡(luò)以較低頻率緩慢變化，可以用低階多項式表示。對于離散信號

式中 ▽為差分算子，s為差分階次，εk（n）為非齊次項。設(shè)多項式階次為2，根據(jù)式（3）得

當(dāng)n=1時，Ak（0）-2Ak（1）+Ak（2）=εk（1）。由于實(shí)際信號中，Ak（0）不存在，因此將Ak（0）視為0，可得-2Ak（1）+Ak（2）=εk（1）。同樣，當(dāng)n=N時，Ak（N+1）不存在，將Ak（N+1）視為0，可得Ak（N-1）-2Ak（N）=εk（N），式中N為采樣點(diǎn)數(shù)。

將式（4）展開成矩陣形式

不同階次的多項式具有相同的矩陣形式

式中M為N×N矩陣。

實(shí)際測試信號y（n）由各階分量與噪聲或誤差δ（n）組成，則

式中y（n）是實(shí)測振動信號，δ（n）為噪聲或誤差項。

為了獲得某一感興趣的信號分量xk（n），由式（7）得

其矩陣形式為

根據(jù)常規(guī)的時頻分析（如短時Fourier變換）的時頻脊線，可以估計感興趣信號分量的瞬時頻率ωk（n），從而獲得其載波矩陣C。

根據(jù)式（6）和（9），應(yīng)用最小二乘法使非齊次項和噪聲或誤差項的平方和最小［16］，即

式中C*為C共軛轉(zhuǎn)置；A*為A共軛轉(zhuǎn)置；r為加權(quán)因子，決定了濾波器的頻率分辨力。r越大，濾波器的分辨力越高，但幅值包絡(luò)的收斂速度越慢；相反，r越小，濾波器的分辨力越低，但幅值包絡(luò)的收斂速度越快。為了保證幅值包絡(luò)的收斂性及計算的可行性，r的取值不能過大，其最大值如表1所示［17］。

表1 加權(quán)因子最大值Tab.1 The maximum value of the weighting factor

由式（10）求得各信號分量的幅值包絡(luò)矩陣A為

由載波矩陣C和幅值包絡(luò)矩陣A可以重構(gòu)感興趣的各信號分量

2 高階能量分離算法

能量算子和能量分離方法完全由信號數(shù)據(jù)驅(qū)動，無需構(gòu)造任何基函數(shù)，能夠適應(yīng)信號的瞬時變化，時間分辨率高，適合檢測信號幅值和頻率的非平穩(wěn)變化特征，結(jié)合時頻分析方法，能夠分析處理任意時變調(diào)制信號。高階能量算子不僅具有普通二階能量算子的各種優(yōu)點(diǎn)，而且對噪聲干擾的魯棒性更好［11，12］，因此 本 文應(yīng)用高 階 能量 算 子進(jìn)行 能 量分離，計算信號的幅值包絡(luò)和瞬時頻率。

對于任意連續(xù)信號x（t），高階能量算子定義 為［10］

式中k為高階能量算子的階次，為任意整數(shù)；x（k）為x（t）的第k階導(dǎo)數(shù)

當(dāng)k為負(fù)數(shù)時，x（k）（t）由微分變?yōu)榉e分，將削弱高階能量算子檢測瞬態(tài)成分的能力，因此本文只關(guān)注k為正數(shù)時的高階能量算子。通常，高階能量算子可由低階能量算子通過遞歸推導(dǎo)得到

由式（13），（14）可知，當(dāng)k=2時，γ2（x）=（˙x）2-x¨x，即為二階Teager能量算子［8］；零階能量算子為γ0（x）=x-x2；一階微分能量算子對于任何信號都為零。在高階能量算子中，三階和四階能量算子比較實(shí)用，分別定義為

其中

對于簡諧信號x（t）=Acos（ωt+θ），可以得到高階能量算子輸出的遞歸方程

初始條件為E0=-A2，E1=0。由式（18）遞推，可得

對于簡諧信號x（t）=Acos（ωt+θ），由四階和二階能量算子可以估計信號的頻率和幅值

對于調(diào)幅-調(diào)頻信號

若幅值包絡(luò)a（t）和瞬時頻率ω（t）變化的速率和幅度相對于載波頻率來說不大，則依然可以由四階和二階能量算子估計信號的時變頻率和包絡(luò)

3 分析過程

（1）應(yīng)用短時Fourier變換等常規(guī)時頻分析方法對信號進(jìn)行初步時頻分析，判斷感興趣的信號成分。根據(jù)時頻脊線估計感興趣信號成分的瞬時頻率。

（2）根據(jù)已獲得的瞬時頻率曲線，構(gòu)造其載波信號，從而獲得感興趣信號成分的載波矩陣。

（3）根據(jù)感興趣信號成分的載波矩陣，應(yīng)用Vold-Kalman濾波將感興趣的信號成分分解為單分量，滿足能量分離方法的要求。

（4）應(yīng)用高階能量分離方法計算各單分量信號成分xi（t）的四階“能量”Ei（t）和瞬時頻率fi（t）。

（5）根據(jù)各單分量信號成分的四階“能量”和瞬時頻率，構(gòu)造信號的高階時頻“能量”分布

式中Ei（t）為第i個單分量在t時刻的四階“能量”，fi（t）為第i個單分量在t時刻的瞬時頻率，δ（·）為Dirac函數(shù)。

4 實(shí)驗信號分析

4.1 實(shí)驗說明

本實(shí)驗在Ottawa大學(xué)機(jī)器監(jiān)控實(shí)驗室完成。圖1為該實(shí)驗室的風(fēng)電機(jī)組行星齒輪傳動實(shí)驗臺，其中行星齒輪箱和定軸齒輪箱均為兩級結(jié)構(gòu)，齒輪參數(shù)分別見表2和3所示。該實(shí)驗臺的具體傳動路線為：電機(jī)→定軸齒輪箱→第1級行星齒輪箱→第2級行星齒輪箱。為了模擬行星齒輪箱的齒輪故障，在第1級太陽輪的某個輪齒上加工了剝落損傷，如圖2所示。在行星齒輪箱箱體頂部安裝了加速度傳感器，采集行星齒輪箱的振動信號。在自然停機(jī)過程中采集了振動信號，采樣頻率為20 k Hz，行星齒輪箱第2級輸出軸承受的負(fù)荷為16.284 N·m。為了計算的可行性，僅取60至40 Hz降速過程內(nèi)的振動信號進(jìn)行分析。

圖1 風(fēng)電機(jī)組行星齒輪傳動實(shí)驗系統(tǒng)Fig.1 Wind turbine planetary gearbox experimental system

表2 行星齒輪箱參數(shù)Tab.2 Planetarygearboxconfigurationparameter

表3 定軸齒輪箱參數(shù)Tab.3 Fixed-shaftgearboxconfigurationparameter

圖2 太陽輪損傷Fig.2 Sun gear damage

根據(jù)行星齒輪箱和定軸齒輪箱的參數(shù)及電機(jī)轉(zhuǎn)速，可計算出第1級行星齒輪箱齒輪特征頻率與電機(jī)轉(zhuǎn)速之間的關(guān)系：

嚙合頻率

太陽輪旋轉(zhuǎn)頻率

行星架旋轉(zhuǎn)頻率

太陽輪局部故障特征頻率

行星輪局部故障特征頻率

齒圈局部故障特征頻率

4.2 實(shí)驗信號分析

4.2.1 正常信號

圖3為正常行星齒輪箱降速過程中驅(qū)動電機(jī)的轉(zhuǎn)速變化曲線。在降速過程中，驅(qū)動電機(jī)的最高轉(zhuǎn)速為60 Hz，相應(yīng)地，第1級行星齒輪箱的最高嚙合頻率為222.222 Hz。由于高階倍頻處特征頻率的變化規(guī)律與基頻處相同，不失一般性，以下分析的頻率范圍為0～400 Hz，覆蓋第1級行星齒輪箱最高嚙合頻率的3/2倍，振動信號中包含該齒輪箱的健康狀態(tài)信息。

圖3 正常狀態(tài)電機(jī)轉(zhuǎn)速Fig.3 Motor speed of normal gearbox

圖4（a）為正常信號的短時Fourier變換譜圖。從圖中可以看出信號成分主要包括嚙合頻率fm、嚙合頻率與太陽輪故障特征頻率之差fm-fs、以及嚙合頻率、太陽輪故障特征頻率的2倍頻和太陽輪旋轉(zhuǎn)頻率的組合fm±2fs-fsr，但是時頻分辨率較低。

圖4 正常信號Fig.4 Normal gearbox signal

應(yīng)用Vold-Kalman濾波分離上述4個主要信號分量，然后應(yīng)用高階能量分離算法計算各分量的4階能量和瞬時頻率，最后應(yīng)用各分量的瞬時頻率和4階能量構(gòu)造高階時頻能量分布，結(jié)果如圖4（b）所示。該方法得到的時頻分布的分辨率高，而且沒有交叉項干擾。可以清晰地識別出嚙合頻率fm、嚙合頻率與太陽輪故障特征頻率之差fm-fs、以及嚙合頻率、太陽輪故障特征頻率的2倍頻和太陽輪旋轉(zhuǎn)頻率的組合fm±2fs-fsr，但它們的幅值均較小。齒輪制造和安裝誤差以及微小的損傷在所難免，導(dǎo)致振動信號中出現(xiàn)上述頻率成分，這些現(xiàn)象符合正常信號的理論預(yù)期結(jié)果。

4.2.2 太陽輪故障信號

圖5為太陽輪故障狀態(tài)下驅(qū)動電機(jī)的轉(zhuǎn)速變化曲線。在降速過程中，驅(qū)動電機(jī)的最高轉(zhuǎn)速仍為60 Hz，因此，以下分析頻率范圍仍為0～400 Hz。

圖5 太陽輪故障狀態(tài)電機(jī)轉(zhuǎn)速Fig.5 Motor speed of faulty sun gearbox

圖6（a）為太陽輪故障信號的短時Fourier變換譜圖。從圖中可以看出信號成分主要包括嚙合頻率fm、嚙合頻率與太陽輪故障特征頻率及其倍頻的組合fm±nfs、以及嚙合頻率、太陽輪故障特征頻率及其倍頻和太陽輪旋轉(zhuǎn)頻率及其倍頻的組合fm±nfs-nfsr等7個頻率成分，但是時頻分辨率較低。

應(yīng)用Vold-Kalman濾波和高階能量分離方法構(gòu)造的高階時頻能量分布如圖6（b）所示。嚙合頻率fm、嚙合頻率與太陽輪故障特征頻率之差fmfs、以及嚙合頻率、太陽輪故障特征頻率的2倍頻和太陽輪旋轉(zhuǎn)頻率的組合fm±2fs-fsr依然存在，但在1～2 s階段，嚙合頻率與太陽輪故障特征頻率之差fm-fs的能量明顯增強(qiáng)。此外，還出現(xiàn)了嚙合頻率與太陽輪故障特征頻率及其3倍頻之和fm+fs，fm+3fs以及嚙合頻率、太陽輪故障特征頻率和太陽輪旋轉(zhuǎn)頻率2倍頻的組合fm+fs-2fsr等3個頻率成分。與正常信號相比，該信號的高階時頻能量分布中出現(xiàn)了更多的頻率成分，且這些頻率均和太陽輪故障直接相關(guān)，說明太陽輪出現(xiàn)了故障，符合實(shí)驗中的實(shí)際情況。

圖6 太陽輪故障信號Fig.6 Faulty sun gearbox signal

5 結(jié) 論

風(fēng)力發(fā)電機(jī)組、直升機(jī)和車輛等裝備中的行星齒輪箱的運(yùn)行工況經(jīng)常變化，振動信號具有非平穩(wěn)性。針對時變工況下行星齒輪箱振動信號的時變調(diào)制特點(diǎn)，發(fā)揮Vold-Kalman濾波單分量分解和高階能量分離時間分辨率高的優(yōu)勢，提出了一種具有良好的時頻聚集能力、無交叉項干擾的時頻分析方法。應(yīng)用該方法分析了行星齒輪箱時變工況下的實(shí)驗信號，準(zhǔn)確識別信號中的時變頻率成分，正確診斷了齒輪故障，驗證了該方法的有效性。

［1］ 雷亞國，何正嘉，林京，等.行星齒輪箱故障診斷技術(shù)的研究進(jìn)展［J］.機(jī)械工程學(xué)報，2011，47（19）：59—67. Lei Yaguo，He Zhengjia，Lin Jing，et al.Research advances of fault diagnosis technique for planetary gearboxes［J］.Journal of Mechanical Engineering，2011，47（19）：59—67.

［2］ Lei Yaguo，Lin Jing，Zuo M J，et al.Condition monitoring and fault diagnosis of planetary gearboxes：a review［J］.Measurement，2014，48：292—305.

［3］ Williams W J，Zalubas E J.Helicopter transmission fault detection via time-frequency，scale and spectral methods［J］.Mechanical Systems and Signal Processing，2000，14（4）：545—559.

［4］ Bartelmus W，Zimroz R.Vibration condition monitoring of planetary gearbox under varying external load［J］.Mechanical Systems and Signal Processing，2009，23：246—257.

［5］ Bartelmus W，Zimroz R.A new feature for monitoring the condition of gearboxes in non-stationary operation conditions［J］.Mechanical Systems and Signal Processing，2009，23：1 528—1 534.

［6］ 馮志鵬，范寅夕，Liang Ming，等.行星齒輪箱故障診斷的非平穩(wěn)振動信號分析方法［J］.中國電機(jī)工程學(xué)報，2013，33（17）：105—110.Feng Zhipeng，F(xiàn)an Yinxi，Liang Ming，et al.A nonstationary vibration signal analysis method for fault diagnosis of planetary gearboxes［J］.Proceedings of the CSEE，2013，33（17）：105—110.

［7］ Feng Zhipeng，Liang Ming，Chu Fulei.Recent advances in time-frequency analysis methods for machinery fault diagnosis：A review with application examples［J］.Mechanical Systems and Signal Processing，2013，38：165—205.

［8］ Kaiser J F.On a simple algorithm to calculate the ‘energy’of a signal［A］.Proceedings of IEEE International Conference on Acoustics，Speech and Signal Processing［C］.Albuquerque，USA.1990，1：381—384.

［9］ Maragos P，Kaiser J F，Quatieri T F.Energy separation in signal modulations with application to speech analysis［J］.IEEE Transactions on Signal Processing，1993，41（10）：3 024—3 051.

［10］Maragos P，Potamianos A.Higher order differential energy operators［J］.IEEE Signal Processing，1995，2 （8）：152—154.

［11］Diop E H S，Boudraa A O，Salzenstein F.A joint 2D AM-FM estimation based on higher order Teager-Kaiser energy operators［J］.Signal Image and Video Processing，2011，5：61—68.

［12］Salzenstein F，Boudraa A O，Chonavel T.A new class of multi-dimensional Teager-Kaiser and higher order operators based on directional derivatives［J］.Multidimensional Systems and Signal Processing，2013，24 （3）：543—572.

［13］Vold H，Herlufsen H.Order tracking at extreme slew rates，using Kalman tracking filters［J］.SAE Paper Number 931288，1993.

［14］Vold H，Mains M，Blough J.Theoretical foundations for high performance order tracking with the Vold-Kalman tracking filter［J］.SAE Paper Number 972007，1997.

［15］Vold H，Herlufsen H，Mains M.Multi axle order tracking with the Vold-Kalman tracking filter［J］. Journal of Sound and Vibration，1997，13（5）：30—34.

［16］Feldbauer Ch，Holdrich R.Realization of a Vold-Kalman tracking filter—A least squares problem［A］. Proceedings of the COST G-6 Conference on Digital Audio Effects（DAFX-000）［C］.Verona Italy，December 7-9，2000：1—4.

［17］Tuma J.Setting the pass bandwidth in the Vold-Kalman order tracking filter［A］.Twelfth International Congress on Sound and Vibration［C］.Lisbon.2005，Paper 719.

Application of Vold-Kalman filter and higher order energy separation to fault diagnosis of planetary gearbox under time-varying conditions

QIN Si-feng1，F(xiàn)ENG Zhi-peng1，LIANG Ming2
（1.School of Mechanical Engineering，University of Science and Technology Beijing，Beijing 100083，China；2.Department of Mechanical Engineering，University of Ottawa，Ottawa，ON K1N 6N5，Canada）

The vibration signals of planetary gearboxes under nonstationary running conditions have significant time-varying modulation feature.Conventional spectral analysis methods are unable to identify the gear fault characteristic frequencies from such nonstationary signals.A method based on Vold-Kalman filter and higher order energy separation is proposed to analyze the vibration signals of planetary gearboxes under nonstationary conditions in time-frequency domain，thus to identify the timevarying characteristic frequencies and diagnose the gear faults.The time-frequency representation derived from Vold-Kalman filter and higher order energy separation provides nice time-frequency resolution and is free from cross term interference，and thus it can effectively pinpoint the time-varying constituent frequency components of nonstationary signals.The proposed method is validated with analysis of lab experimental signals of a planetary gearbox under time-varying running conditions.

fault diagnosis；planetary gearbox；higher order energy separation；Vold-Kalman filter；time-frequency analysis

TH165+.3；TH132.425

A

1004-4523（2015）05-0839-07

10.16385/j.cnki.issn.1004-4523.2015.05.20

秦嗣峰（1989—），男，碩士研究生。電話：18866226821；E-mail：qsf0312@163.com

馮志鵬（1973—），男，博士，教授，博士生導(dǎo)師。電話：13520978607；E-mail：fengzp@ustb.edu.cn

2014-05-13<； class="emphasis_bold">；修訂日期：2；

2014-11-05

國家自然科學(xué)基金資助項目（11272047，51475038）和教育部新世紀(jì)優(yōu)秀人才支持計劃項目（NCET-12-0775）