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

(重慶大學(xué) 機(jī)械傳動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，重慶 400030)

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

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

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

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

1 Vold-Kalman階比跟蹤原理

階比定義為參考軸每轉(zhuǎn)內(nèi)發(fā)生的旋轉(zhuǎn)振動(dòng)次數(shù)，可以設(shè)定為一個(gè)幅值和頻率隨時(shí)間變化的正弦函數(shù)，階比的頻率和參考軸的旋轉(zhuǎn)頻率相關(guān)聯(lián)。階比信號(hào)可以表示為幅值和載波的乘積［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的主要目標(biāo)是獲得復(fù)包絡(luò)ak(t)或者離散形式的ak(n)

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

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

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

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

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

其矩陣形式為:

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

其矩陣形式為:

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

2 基于能量重心校正法的瞬時(shí)頻率估計(jì)

2.1 瞬時(shí)頻率原理

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

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

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

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

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

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

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

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

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

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

比較式(13)和式(14)發(fā)現(xiàn):式(14)是式(13)的近似離散表達(dá)形式，也就是說(shuō)通過(guò)能量重心法校正獲得的頻率近似等于瞬時(shí)頻率?；谀芰恐匦姆↖FE的基本思路是根據(jù)假定的瞬時(shí)頻率初始值進(jìn)行整周期采樣，同時(shí)做加窗的DFT;然后利用能量重心法對(duì)該段信號(hào)做離散頻譜校正，得到新的瞬時(shí)頻率;依次迭代重復(fù)，最后得到準(zhǔn)確的瞬時(shí)頻率值。該方法的誤差主要來(lái)自信號(hào)的截?cái)嗾`差，所以采用整周期采樣來(lái)減小誤差，同時(shí)對(duì)信號(hào)加窗，以便提高精度和減小各相鄰譜線之間的干擾。

下面給出基于能量重心法瞬時(shí)頻率估計(jì)的計(jì)算步驟。

(1)對(duì)振動(dòng)信號(hào)進(jìn)行分段，數(shù)據(jù)長(zhǎng)度為L(zhǎng)，每M個(gè)點(diǎn)為一段，重疊長(zhǎng)度為 l，共分為 m段，m∈［0，Num－1］，Num=int［L/(M －1)］=int(L/dN)，式中 int表示對(duì)計(jì)算結(jié)果取整;

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

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

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

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

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

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

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

3 基于能量重心法瞬時(shí)頻率估計(jì)的VKF-OT

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

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

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

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

4.1 仿真分析

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

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

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

圖2 仿真信號(hào)時(shí)域波形圖Fig.2 Time waveform of synthetic signal

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

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

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

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

圖4 能量重心法IFE和STFT瞬時(shí)頻率估計(jì)對(duì)比Fig.4 Comparison of IFE by energy centrobaric method and STFT

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

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

圖5 基于能量重心法瞬時(shí)頻率估計(jì)的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階比分量與理想信號(hào)的前1000點(diǎn)對(duì)比結(jié)果Fig.7 Comparison between order－1 of axis－I and designed signal for n=1000

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

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

圖8 電機(jī)升速振動(dòng)試驗(yàn)裝置Fig.8 Test setup for vibration measurement during motor's run-up

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

5 結(jié)論

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

圖10 偏心電機(jī)升速階段振動(dòng)信號(hào)第4階比分量提取Fig.10 Extracted order-4 waveform of motor's vibration signal

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

［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］孔慶鵬，宋開(kāi)臣，陳 鷹.最小二乘自適應(yīng)濾波旋轉(zhuǎn)機(jī)械階比跟蹤研究［J］.浙江大學(xué)學(xué)報(bào)(工學(xué)版)，2003，40(9):1648－1651.

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

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

［7］趙曉平，張令彌，郭勤濤.基于瞬時(shí)頻率估計(jì)的自適應(yīng)Vold-Kalman階比跟蹤研究［J］.振動(dòng)與沖擊，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］.振動(dòng)工程學(xué)報(bào)，2001，14(3):354－358.