鄧 蕾，傅煒娜，董紹江，湯寶平

(重慶大學(xué) 機械傳動國家重點實驗室，重慶 400030)

階比跟蹤(OT)技術(shù)是旋轉(zhuǎn)機械升降速過程非平穩(wěn)振動信號分析的有效方法，在旋轉(zhuǎn)機械的非平穩(wěn)振動信號分析和故障診斷中占有重要地位。

1993年Vold和Leuridan［1］首次提出了基于卡爾曼濾波的階比跟蹤的方法，開創(chuàng)了通過時域濾波方法進行階比跟蹤的新思維。目前，丹麥的B＆K公司將VKF-OT方法作為其階比分析商業(yè)化軟件的主要技術(shù)，并得到了廣泛的應(yīng)用。文獻［2］對VKF-OT理論作了進一步的詳細研究，提出了角速度和角位移兩種VKFOT方法，對兩者之間的區(qū)別進行了詳細地描述。在此基礎(chǔ)上文獻［3］提出和實現(xiàn)了自適應(yīng)VKF-OT方法。另外，國內(nèi)對此也作了一些研究［4，5］。

然而，VKF-OT方法仍需要通過安裝轉(zhuǎn)速計來獲取轉(zhuǎn)速信號，在不方便安裝鑒相裝置的場合限制了應(yīng)用。針對簡化階比跟蹤的實現(xiàn)途徑和對硬件安裝要求的問題，文獻［6］提出了基于短時傅里葉變換的峰值估計瞬時頻率的階比跟蹤方法，結(jié)合Gabor時頻濾波實現(xiàn)了階比分量的提取。文獻［7］提出了基于STFT-VA瞬時頻率估計的自適應(yīng)VKF-OT方法。但這些方法主要是在旋轉(zhuǎn)機械振動信號的時頻譜圖中通過峰值搜索獲得某階比分量的瞬時頻率，進而得到所需要的參考軸轉(zhuǎn)速，所以分析精度受到頻率分辨率的限制，只有在特定條件下才能應(yīng)用于多分量信號。

本文從瞬時頻率估計理論出發(fā)，結(jié)合離散頻譜校正技術(shù)中的能量重心法，通過研究兩者之間的聯(lián)系，改進了瞬時頻率估計方法，再將其應(yīng)用于VKF-OT方法中，提出了一種無轉(zhuǎn)速計的基于能量重心法瞬時頻率估計的VKF-OT方法。

1 Vold-Kalman階比跟蹤原理

階比定義為參考軸每轉(zhuǎn)內(nèi)發(fā)生的旋轉(zhuǎn)振動次數(shù)，可以設(shè)定為一個幅值和頻率隨時間變化的正弦函數(shù)，階比的頻率和參考軸的旋轉(zhuǎn)頻率相關(guān)聯(lián)。階比信號可以表示為幅值和載波的乘積［4］:

式中，k為基本頻率或參考軸轉(zhuǎn)速的倍數(shù)，表示為被跟蹤的階比;ak(t)表示復(fù)包絡(luò)，代表了階比幅值的變化，a－k(t)是 ak(t)的復(fù)共軛 ak(t)=a－k(t)*。Θk(t)為載波，表示為:

式中，ω(τ)為參考軸的角頻率，∫t0ω(τ)dτ 為角位移。式(2)的離散形式表示為:

VKF-OT的主要目標是獲得復(fù)包絡(luò)ak(t)或者離散形式的ak(n)

由式(1)知，復(fù)包絡(luò)ak(t)是載波Θk(t)的低頻幅值調(diào)制，由于系統(tǒng)的固有慣量，使得階比幅值平滑變化，可以用低階多項式表示，則系統(tǒng)的狀態(tài)方程為:

式中，▽表示不同的算子，s為給定階數(shù)，ε(n)為非一致項。

實際測得的振動信號y(n)由各階比分量xk(n)的總和加上測量誤差和噪聲ξ(n)組成，則系統(tǒng)的觀測方程為:

為了跟蹤某一階比分量xk(n)，用VKF-OT方法來計算復(fù)包絡(luò)ak(n)，則狀態(tài)方程和觀測方程給出以下形式:

根據(jù)式(6)，展開一階狀態(tài)方程:

其矩陣形式為:

其中，M為N×N－1的矩陣，N為采樣點數(shù)。同樣，根據(jù)式(7)，系統(tǒng)的觀測方程為:

其矩陣形式為:

根據(jù)系統(tǒng)的狀態(tài)方程(8)和觀測方程(10)，通過最小二乘法，使得非一致相ε(n)以及測量誤差和噪聲ξ(n)的平方和最?。?］，求得未知量 ak(n)。

2 基于能量重心校正法的瞬時頻率估計

2.1 瞬時頻率原理

瞬時頻率(IF)是非平穩(wěn)信號分析中的一個重要概念，Ville于1948年提出的瞬時頻率定義式是目前在學(xué)術(shù)界最為常用且得到了普遍認可［6，9］，即若實信號為s(t)=a(t)cosφ(t)，則 IF 定義為:

另一種等價形式為時頻分布(TFD)的一階矩，即:

2.2 能量重心離散頻譜校正法

能量重心校正法是由三點卷積幅值校正法發(fā)展起來的。該方法可以同時對離散頻譜的頻率、幅值和相位進行校正［10］。能量重心法是通過利用各種窗函數(shù)離散頻譜的能量重心無窮逼近坐標原點的特點進行離散頻譜校正的。這種方法在分析平穩(wěn)或準平穩(wěn)信號時，具有較高的分析精度。

設(shè)Gi為加Hanning窗的諧波信號的功率譜第i條譜線值，Gk為主瓣內(nèi)譜線最大值，k為幅值最大點對應(yīng)的譜線號，采樣頻率為fs，采樣點數(shù)為N，f0為主瓣中心，得到能量校正法校正頻率的通用公式:

根據(jù)帕賽瓦定理，設(shè)能量恢復(fù)系數(shù)為Kt，則校正后的幅值為:

離散頻率校正方法中的能量法對所有對稱窗函數(shù)都適合，校正精度與窗函數(shù)有關(guān)，加Hanning窗時具有較高的校正精度。負頻率成分和間隔較近的多頻率成分產(chǎn)生的干涉現(xiàn)象所帶來的誤差對精度的影響小。不考慮信號中噪聲的影響，是一種精度較高的近似校正方法。

2.3 基于能量重心法的瞬時頻率估計(IFE)的原理

在STFT時頻面上進行峰值搜索估計瞬時頻率是目前獲取旋轉(zhuǎn)機械升降速階段振動信號的瞬時頻率的主要方法［6］。但是這種方法的估計精度依賴于頻率分辨率，鄰近頻率成分的干擾也會對估計精度產(chǎn)生影響，存在時頻集聚性不是很理想的問題。另外，在求取多分量信號的IF時，一般通過時頻濾波將其它分量遮掩，提取某一分量的IF，依次遞推，增加了算法的計算量。

離散頻譜校正方法可以改善由于時域截斷產(chǎn)生的能量泄漏的問題，提高分析精度。而旋轉(zhuǎn)機械升降速過程的振動信號屬于非平穩(wěn)信號，不能直接使用離散頻譜校正方法。因此，根據(jù)STFT算法原理，把短時間間隔的信號看成是準平穩(wěn)信號，就能夠利用離散頻譜校正方法來校正信號的頻譜。這樣可以在頻域中直接通過振動信號的頻譜進行瞬時頻率估計，并通過離散頻譜校正方法提高分析精度，進而實現(xiàn)階比跟蹤。

比較式(13)和式(14)發(fā)現(xiàn):式(14)是式(13)的近似離散表達形式，也就是說通過能量重心法校正獲得的頻率近似等于瞬時頻率。基于能量重心法IFE的基本思路是根據(jù)假定的瞬時頻率初始值進行整周期采樣，同時做加窗的DFT;然后利用能量重心法對該段信號做離散頻譜校正，得到新的瞬時頻率;依次迭代重復(fù)，最后得到準確的瞬時頻率值。該方法的誤差主要來自信號的截斷誤差，所以采用整周期采樣來減小誤差，同時對信號加窗，以便提高精度和減小各相鄰譜線之間的干擾。

下面給出基于能量重心法瞬時頻率估計的計算步驟。

(1)對振動信號進行分段，數(shù)據(jù)長度為L，每M個點為一段，重疊長度為 l，共分為 m段，m∈［0，Num－1］，Num=int［L/(M －1)］=int(L/dN)，式中 int表示對計算結(jié)果取整;

(2)旋轉(zhuǎn)機械升降速過程中總有一穩(wěn)定轉(zhuǎn)速(相對恒定)階段，比如，最高轉(zhuǎn)速為9000 r/min的降速階段，對應(yīng)的一階最大轉(zhuǎn)頻為150 Hz。所以選擇最高速階段作為初始分析段，根據(jù)最大轉(zhuǎn)頻和該段的頻譜成分設(shè)定初始轉(zhuǎn)頻f0，局部極大值搜索范圍q和整周期采樣系數(shù)K;

(3)根據(jù)當前m，計算當前分段數(shù)據(jù)的起始點位置:p=m×dN。在采樣信號中截取M點數(shù)據(jù){s(i)}，i∈［p，M －1+p］};

(4)根據(jù)f0和系數(shù)R計算出整周期采樣點數(shù)N0，以分析段中點為中心重新選取N0點作為新的分析段，并對其做加窗的DFT，求得搜索范圍q內(nèi)最大頻率譜線位置以及其左右鄰近的n條譜線的功率值。其中，n的取值與加窗的類型有關(guān)，這里加Hanning窗，n=2;

(5)利用能量重心法對搜索到的峰值頻率譜線進行頻譜校正，得到校正頻率f，令f為新的轉(zhuǎn)頻，重新計算出整周期采樣所需的采樣點數(shù)N;

(6)如果N=N0或者達到最大迭代次數(shù)j時，f為分析段中點所對應(yīng)的頻率估計值，否則，f0=f，轉(zhuǎn)步驟(4);

(7)令m=m+1，返回(3)，直到m=Num－1。

最終根據(jù)迭代求得的參考軸轉(zhuǎn)頻f(m)，插值擬合后得到瞬時頻率擬合值。

3 基于能量重心法瞬時頻率估計的VKF-OT

從前面的討論可以看出，通過能量重心法得到的參考軸轉(zhuǎn)頻近似等于瞬時頻率，由此可以獲得參考軸轉(zhuǎn)速，進而實現(xiàn)階比跟蹤。本文方法的實現(xiàn)過程如圖1所示。

獲取旋轉(zhuǎn)機械升降速階段的振動信號，對該振動信號作基于能量重心法的瞬時頻率估計，獲得參考軸瞬時頻率估計值;插值擬合后得到連續(xù)的參考軸瞬時頻率，并對各個振動采樣點進行階比相位估計;然后，設(shè)置加權(quán)因子r，選取階次p等參數(shù)，進行Vold-Kalman自適應(yīng)濾波;最后得到所提取的第p階信號分量。

圖1 基于能量重心法的瞬時頻率VKF-OT流程圖Fig.1 Flow diagram of VKF-OT based on the improved IFE with energy centrobaric method

4 仿真與應(yīng)用實例分析

4.1 仿真分析

為了驗證算法，設(shè)計一個包含兩個獨立參考軸轉(zhuǎn)速信號的仿真信號來模擬旋轉(zhuǎn)機械升速階段的振動信號。采樣點數(shù)N=10 k，采樣頻率fs=1 kHz。仿真信號的時域波形如圖2所示。

式中，η(n)為高斯白噪聲［η(n)～N(0，1)］。此仿真信號由兩組信號分量組成，第Ⅰ組包含第1階、第3階分量，基準頻率從15 Hz到150 Hz線性變化，·Δt表示第n個采樣時刻第Ⅰ轉(zhuǎn)軸轉(zhuǎn)過的相位角。第1、3階比分量的幅值分別呈0到10和5到15線性變化。第Ⅱ組為添加了高斯白噪聲的正弦信號，旋轉(zhuǎn)頻率f2(n)固定為 200 Hz，幅值固定為 10，·Δt表示第n個采樣時刻第Ⅱ轉(zhuǎn)軸轉(zhuǎn)過的相位角。

采用基于能量重心法的瞬時頻率估計方法對仿真的振動信號進行瞬時頻率估計。分析中，加Hanning窗，n=2，相鄰兩段時移點數(shù)為500，結(jié)果如圖3所示。

圖2 仿真信號時域波形圖Fig.2 Time waveform of synthetic signal

圖3 基于能量重心法的瞬時頻率估計Fig.3 IFE based on energy centrobaric method

圖3中圓點表示采用基于能量重心法的瞬時頻率估計得到的估計值，實線表示設(shè)計的參考軸Ⅰ的瞬時頻率曲線，虛線為對估計值采用插值擬合后得到的擬合曲線。從圖中可以看出，采用基于能量重心法得到的轉(zhuǎn)頻理論上近似于瞬時頻率，與仿真信號的理想值非常吻合。

分別對基于STFT峰值搜索法和本文方法獲得的瞬時頻率估計值求相對于理想瞬時頻率值的百分比誤差。

式中，f(i)為瞬時頻率的理想值，～f(i)為f(i)的擬合值，根據(jù)式(17)求得，基于STFT峰值搜索得到的瞬時頻率估計值相對于理想瞬時頻率值的百分比誤差ζ1=0.37%，基于能量重心法的瞬時頻率估計得到的擬合值相對于理想瞬時頻率值的百分比誤差ζ2=0.03%，對比結(jié)果表明，由本文提出的基于能量重心法IFE誤差小，估計的瞬時頻率逼近真實的瞬時頻率值。取移動步長為1，分別用這兩種方法計算得到瞬時頻率估計值，取某一段的數(shù)據(jù)如圖4所示。圖中結(jié)果再次說明，本文方法連續(xù)性好，分析精度受頻率分辨率的影響小。

圖4 能量重心法IFE和STFT瞬時頻率估計對比Fig.4 Comparison of IFE by energy centrobaric method and STFT

然后，采用基于能量重心法IFE獲得的瞬時頻率擬合值作為參考軸對應(yīng)的瞬時頻率，并求得相應(yīng)的瞬時角位移信號，對振動信號進行VKF-OT分析。本方法中使用一階自適應(yīng)Vold-Kalman濾波器分離各個階比分量。在圖5中，圖(a)、圖(b)分別為第Ⅰ參考軸的第1、3階比分量提取結(jié)果，圖(c)為第Ⅱ參考軸的200 Hz常頻分量提取結(jié)果。圖6為兩個軸各自階比分量的幅值，可以看出由VKF-OT分析得到的幅值與設(shè)計值相符合。

為了能夠清楚地看到提取結(jié)果與理想信號之間的差異，取第Ⅰ參考軸的第1階分量的前1000個點與設(shè)計信號進行對比。如圖7所示，實線表示用本文方法所提取的第1階比分量，虛線表示理想信號，從圖中可以看到，除了起始點處有一些偏差外，提取的階比幅值都和設(shè)計值吻合。

圖5 基于能量重心法瞬時頻率估計的VKF-OT各階比分量提取結(jié)果Fig.5 Extracted orders components using VKF-OT based on the improved IFE with energy centrobaric method

圖6 VKF-OT分析得到的Ⅰ、Ⅱ軸各階比分量幅值Fig.6 Amplitudes of orders about axes－ ⅠandⅡtracked by using VKF-OT

圖7 第I參考軸的第1階比分量與理想信號的前1000點對比結(jié)果Fig.7 Comparison between order－1 of axis－I and designed signal for n=1000

4.2 應(yīng)用實例分析

下面對某一偏心直流電機在升速階段引起的懸臂梁結(jié)構(gòu)的振動信號進行測試。試驗裝置如圖8所示，直流電機安裝在簡支梁振動臺上，振動信號用YD－37型加速計拾取。采樣頻率為10 kHz，采樣長度為75 k點。首先，通過基于能量重心離散頻譜校正法的瞬時頻率估計法提取振動信號的瞬時頻率估計值，并通過三次樣條擬合獲得瞬時頻率曲線。

圖8 電機升速振動試驗裝置Fig.8 Test setup for vibration measurement during motor's run-up

結(jié)果如圖9所示，圓點表示瞬時頻率的估計值，實線表示三次樣條插值擬合后的瞬時頻率曲線。然后，利用參考軸瞬時轉(zhuǎn)速與其瞬時頻率的對應(yīng)關(guān)系，求得參考軸的轉(zhuǎn)速信號。最后，對原始振動信號進行VKFOT分析，提取第4階比分量，結(jié)果如圖10所示。分析結(jié)果表明，通過本文方法在提取振動信號瞬時頻率的同時，能夠在時域中直接提取所感興趣的某階比分量。

5 結(jié)論

圖9 能量重心法瞬時頻率估計的結(jié)果Fig.9 IFE based on energy centrobaric method

圖10 偏心電機升速階段振動信號第4階比分量提取Fig.10 Extracted order-4 waveform of motor's vibration signal

本文分析了瞬時頻率估計理論和離散頻譜校正方法中的能量重心法，根據(jù)兩者之間的理論聯(lián)系，通過用能量重心法對信號的頻譜做校正，得到了信號的瞬時頻率，再采用Vold-Kalman階比跟蹤提取感興趣的某階分量，從而實現(xiàn)無轉(zhuǎn)速計的旋轉(zhuǎn)機械Vold-Kalman階比跟蹤。因為該方法根據(jù)獲得的瞬時頻率對原分段信號進行了整周期采樣處理，使得頻率分辨率能夠自適應(yīng)的變化，同時能量重心校正法的引入可以對誤差進行了修正，提高了分析精度并減小了各相鄰譜線之間的干擾。仿真分析和應(yīng)用實例分析表明該方法的有效性，在精確提取參考軸轉(zhuǎn)速信號的同時，提取了階比分量，是對現(xiàn)有技術(shù)的一個有力補充。

［1］Vold H，Herlufson H，Mains M.Multi axle order tracking with the Vold-Kalman tracking filter［J］.Sound and Vibration Magazine，1997，13(5):30 －34.

［2］Pan M Ch，Lin Y F.Further exploration of Vold-Kalman filtering order tracking with shaft-speed information-I:Theoretical Part，numerical implementation and parameter investigations［J］.Mechanical Systems and Signal Processing，2006，20:1134－1154.

［3］Pan M Ch，Wu C X.Adaptive Vold-Kalman filtering order tracking［J］.Mechanical Systems and Signal Processing，2007，21:2957 －2969.

［4］孔慶鵬，宋開臣，陳 鷹.最小二乘自適應(yīng)濾波旋轉(zhuǎn)機械階比跟蹤研究［J］.浙江大學(xué)學(xué)報(工學(xué)版)，2003，40(9):1648－1651.

［5］孔慶鵬，劉敬彪，章雪挺，等.多軸機械自適應(yīng)濾波交叉階比跟蹤［J］.農(nóng)業(yè)機械學(xué)報，2009，40(1):213－216.

［6］郭 瑜，秦樹人，湯寶平，等.基于瞬時頻率估計的旋轉(zhuǎn)機械階比跟蹤［J］.機械工程學(xué)報，2003，39(3):32－35.

［7］趙曉平，張令彌，郭勤濤.基于瞬時頻率估計的自適應(yīng)Vold-Kalman階比跟蹤研究［J］.振動與沖擊，2008，27(12):112－116.

［8］Feldauer Ch，Holdrich R.Realisation of a Vold-Kalman Tracking Filter-A Least Square Problem［C］.Proceedings of the COST G-6 Conference on Digital Audio Effects(DAFX－410800)，Verona Italy，2000，December，7－9.

［9］Boashash B.Interpreting and estimating the instantaneous frequency of a signal—Part I:Fundamentals［J］.Proceedings of IEEE，1992，80(4):520－538.

［10］丁 康，江利旗.離散頻譜的能量重心校正法［J］.振動工程學(xué)報，2001，14(3):354－358.