劉永強， 郝高巖， 廖英英， 楊紹普

(1.石家莊鐵道大學機械工程學院 石家莊，050043)(2.石家莊鐵道大學土木工程學院 石家莊，050043)

10.16450/j.cnki.issn.1004-6801.2017.05.025

一種基于Viterbi法的改進瞬時轉速估計算法

劉永強1， 郝高巖1， 廖英英3， 楊紹普1

(1.石家莊鐵道大學機械工程學院 石家莊，050043)(2.石家莊鐵道大學土木工程學院 石家莊，050043)

針對無轉速計的瞬時轉速估計問題及現(xiàn)有方法在抗噪、抗臨近階次和實時性方面的不足，基于Viterbi算法提出了一種改進型瞬時頻率估計(instantaneous frequency estimation,簡稱IFE)方法，并將其應用于變轉速工況下滾動軸承的瞬時轉速估計。將隱馬爾科夫模型中的Viterbi算法引入轉頻估計，分析了某一時頻平面中代價函數(shù)的計算次數(shù)；根據(jù)歐幾里得距離函數(shù)的性質，提出了代價函數(shù)迭代循環(huán)停止的新型判定準則。該準則的優(yōu)點在于可以快速搜索時頻平面，尋找到最優(yōu)的局部路徑，提高了IFE精度和計算速度。通過仿真信號和實驗數(shù)據(jù)對該算法進行驗證。結果表明：改進后算法效率明顯提高，和基于峰值搜索IFE方法相比，改進的Viterbi-IFE方法具有較高的精度和穩(wěn)定性。

滾動軸承； 振動信號； 瞬時頻率估計； Viterbi算法； 代價函數(shù)

引 言

階次分析為變轉速工況下滾動軸承故障診斷的重要方法之一[1]，該方法中的關鍵問題就是準確追蹤參考軸轉速。一般情況下可采用轉速計等硬件裝置測量轉速，但在某些復雜工況下轉速計等鍵相裝置不易安裝，難以跟蹤參考軸轉速。針對該問題本研究改進了瞬時頻率估計方法，從測量的振動信號中提取出轉頻信息，繼而實現(xiàn)無需轉速計的階次分析。

陳平等[2]研究了基于MUSIC的瞬時頻率估計方法，但存在相位重疊且只適用于單分量信號；為實現(xiàn)無轉速計的階次跟蹤，郭瑜等[3-5]提出對采集的振動信號進行短時傅里葉變換得到時-頻圖，再通過局部峰值搜索獲得瞬時頻率(instantaneous frequency,簡稱IF)，但該方法的抗噪性及抗臨近階次的能力不足；趙曉平等[6-7]提出基于Viterbi算法的瞬時頻率估計方法，但由于該算法較復雜，運算時間較長，無法滿足實時分析的要求。Leclere等[8]提出了一種瞬時轉速估計算法，不需要進行先驗假設，而是由信號瞬時譜產生的概率密度函數(shù)來確定轉速信息?；谒矔r時間尺度因子估計，Combet等[9]提出了一種瞬時轉速相對波動估計方法，該方法采用峭度最大值法獲得分析信號的最佳長度，基于短時尺度變換估計瞬時轉速信息。Urbanek等[10]基于相空間重構和時頻分析，提出一種兩步法用于估計瞬時頻率，可處理劇烈波動的轉速或復雜譜特征信號，但這些算法大多對抗臨近階次和實時性方面沒有過多要求。

針對目前算法的不足，筆者研究了Viterbi算法的路徑尋優(yōu)過程及其迭代次數(shù)，改進了該算法的尋優(yōu)方法，并提出了迭代循環(huán)停止的新型判定準則。改進后的方法提高了瞬時頻率的估算精度和計算速度，為實現(xiàn)無轉速計的實時階次跟蹤打下基礎。

1 基于Viterbi算法IFE的理論基礎

Viterbi算法是針對隱馬爾科夫模型提出的[11]：給出一個觀測序列O1，O2，O3，…，希望找到觀測序列背后的隱藏狀態(tài)序列P1，P2，P3，…，可采用動態(tài)規(guī)劃的方法來尋找出現(xiàn)概率最大的隱藏狀態(tài)序列。

基于Viterbi算法的IFE方法的基本思想來源于數(shù)字圖像處理中的邊線追蹤(edge-following)算法[12]，其原理是連接地圖(時頻平面)上的點，使這些點形成的路徑的長度值和高度變化值之和盡可能小，這和隱馬爾科夫模型中尋找觀測序列背后最可能的隱藏狀態(tài)序列相同[13]。

在一個由N組時間點和M組頻率點組成的時頻網(wǎng)格中，求出從第一組時間點到最后一組時間點的最優(yōu)路徑。設相鄰時間點n1，n2之間的所有路徑集合為K。n1，n2之間的最佳路徑表達式為

(1)

其中：p(k(n);n1,n2)為沿路徑k(n)的代價函數(shù)，其為g(x,y)和f(x)之和，時間從n1到n2；函數(shù)g(x,y)=g(|x-y|)為關于|x-y|的非增函數(shù)，其中x和y分別為路徑k(n)和k(n-1)所處的網(wǎng)格位置頻率編號，則g(x,y)表示x、y兩個連續(xù)時間點之間的距離；f(x)為關于x=Sf(n,k(n))的非減函數(shù)，其中函數(shù)Sf(x,y)為短時傅里葉變換(shorttimeFouriertransformation,簡稱STFT)后所得時頻平面的幅值。

通過這樣的定義，時頻平面中的頻率點所對應的STFT值越大，其成為瞬時頻率估計點的概率就越大。函數(shù)f(x)定義如下：考慮時刻n，把該時刻的STFT值Sf(n,w)按從大到小的順序排列

Sf(n,ω1)≥Sf(n,ω2)≥…≥Sf(n,ωj)≥

…≥Sf(n,ωM)

(2)

其中：j∈[1,M]，M為信號頻率點數(shù)。

f(x)定義為f(Sf(n,ωj))=j-1，該定義使IFE值第j個最大值的概率隨著j的增大而線性降低，可保證Sf(n,w)的較大者作為估計結果[14]。

下面對g(x,y)的形式進行討論。假設信號在持續(xù)時間范圍內IF值沒有發(fā)生突變，可將g(x,y)的形式表達為

(3)

若整個信號持續(xù)時間范圍內IF有突變點，g(x,y)的形式為

(4)

分析式(3)，(4)，當兩連續(xù)時間點的IF變化小于Δ1，此時距離函數(shù)g(x,y)為零，沒有距離函數(shù)引起的代價。當Δ1→∞時，此時g(x,y)恒為零，估計結果為STFT的最大值。當兩連續(xù)點的瞬時頻率有突變時，即|x-y|>Δ2，g(x,y)=c(Δ2-Δ1)，此時g(x,y)不為|x-y|的函數(shù)，而為定值c(Δ2-Δ1)，將距離函數(shù)限制在一定范圍內，準確追蹤到突變點。另外，c和Δ的選擇需要根據(jù)預提取信號的波動情況而定，波動越大Δ越大。

2 基于Viterbi算法的改進

2.1 改進的迭代算法

在一個由M組頻率點和Q組時間點組成的時頻網(wǎng)格中，即T={(ni,ωj)|i∈[1,Q],j=[1,M]}，整個時頻平面需要計算的路徑總數(shù)為MQ，但在點數(shù)比較多的情況下，計算量非常大，實現(xiàn)起來難度很大。文中采用迭代方法簡化計算，并提出代價函數(shù)循環(huán)停止的判定準則，改進的迭代過程如下。

1) 假設從時刻n1到時刻(ni,ωj)的最優(yōu)路徑Li(n,ωj)已經(jīng)確定，最優(yōu)路徑表達式為

(5)

其中：n∈[n1,ni]；j∈[1,M]；Kij為從n1到(ni,ωj)的路徑的集合；p(k(n);n1,(ni,ωj))為沿路徑k(n)的代價函數(shù)，最優(yōu)路徑為代價函數(shù)最小的那條。

這里對于ni時刻的不同的頻率點，均存在一條最優(yōu)路徑，共M條最優(yōu)路徑。則在時間段[n1,ni]中，估計的瞬時頻率曲線就取M條最優(yōu)路徑中代價函數(shù)最小的那條。表達式為

(6)

2) 迭代到下一點ni+1的局部最優(yōu)路徑可以表示為

Li+1(n，ωj)=[Li(n，ωl),(ni+1,ωj)]

(j∈[1,M]；l∈[1,M])

(7)

新的代價函數(shù)為

p(li+1(n,ωl);n1,(ni+1,ωj))=

p(Li(n,ωl);n1,(ni,ωj))+g(ωl,ωj)+

f(Sf(ni+1,ωj))

(8)

最優(yōu)路徑取代價函數(shù)值最小的路徑。這里Li+1(n;ωj)對于不同頻率點也存在M個局部最優(yōu)路徑。任一時刻均有M個頻率點，計算每個頻率點的局部最優(yōu)路徑需要計算M次，因此計算出該時刻所有點的局部最優(yōu)路徑共計算M2次，對于整個時頻平面Q組時間點，共計算(Q-1)×M2次，遠遠小于之前MQ的計算量。

為了進一步減小計算量，觀察代價函數(shù)表達式(8)，p(Li(n;ωl);n1,(ni,ωj))為前ni組時間點局部最優(yōu)路徑的代價函數(shù)，g(ωl,ωj)為歐幾里得距離|ωl-ωj|的單調增函數(shù)。計算局部最優(yōu)路徑時，如果每組時間點中的每個頻率點ωj(j∈[1,M])不變，那么所對應的函數(shù)值f(Sf(ni+1,ωj))不變。因此采用以下步驟進一步減小計算量。

①令ρ=Δ；

3) 按步驟2)計算每個頻率點的局部最優(yōu)路徑。將該時間點的最優(yōu)路徑代價函數(shù)作為結果，帶入下一時間點代價函數(shù)的計算中，依次迭代，求出整個時頻平面的最優(yōu)路徑。

2.2 迭代算法仿真算例

假設由M=8和Q=3組成的時頻網(wǎng)格的一個元素為fij=f(Sf(ni,wj))，如圖1所示。

圖1 M=8，Q=3迭代過程仿真算例Fig.1 Iterative process when M=8，Q=3

由式(3)計算g(ωi,ωj)=gij，其中c=2.5，Δ=1。連接某時刻的某一頻率點與其前一時刻距離最近的3個頻率點，即|i-j|≤1時，gij=0。當|i-j|>1時，gij=2.5×|i-j|。

3 算法驗證

3.1 線性掃頻信號驗證

在軸承故障診斷中，軸承振動信號為含噪聲的多分量信號，經(jīng)過預處理后其噪聲仍然很大，且可能出現(xiàn)臨近階次。為檢驗基于Viterbi算法的IFE方法的抗噪和抗臨近階次能力，采用仿真信號進行研究，其瞬時頻率如式(9)所示。

ω(t)=2.5πt+30π

(9)

通過該頻率調制規(guī)律來仿真滾動軸承升速工況振動信號，其多分量信號模型[15]為

(10)

其中：η(t)為高斯白噪聲；信噪比為-9 dB。

時域波形如圖2所示。轉頻從15 Hz線性變化到40 Hz，持續(xù)時間20 s。仿真模型中包含轉頻的一倍頻，幅值為1；0.66倍頻，幅值為0.8；0.5倍頻，幅值為0.7。STFT時頻圖如圖3所示，可以觀察到由于噪聲影響，各倍頻分量并不突出。

圖2 線性掃頻時域信號Fig.2 Linear sweep frequency signal in time domain

分別采用局部峰值搜索算法和基于Viterbi算法的IFE方法提取轉頻，結果如圖4所示。圖中星點線為基于Viterbi算法的IFE方法得到的轉頻曲線，短畫線為局部峰值搜索法得到的轉頻曲線，實線為理論轉頻曲線。

圖4 線性掃頻信號IFE結果對比圖Fig.4 IFE results comparison for linear sweep frequency signal

根據(jù)式(11)求兩種算法的瞬時頻率估計值相對于理論值的百分比誤差。

(11)

局部峰值搜索算法的誤差ξ=43.4%，基于Viterbi算法的IFE誤差ξ=3.11%?？梢钥闯鲈诟咴肼暫统霈F(xiàn)臨近階次的情況下，局部峰值搜索算法已經(jīng)失效，抗噪和抗臨近階次的能力有限。而基于Viterbi算法的IFE結果具有很高的精度，在較高噪聲干擾下仍可精確進行估計。

3.2 正弦掃頻信號

仿真強噪聲環(huán)境下的正弦掃頻信號，并計算其相對于理論值的百分比誤差。瞬時頻率由式(12)給出。

ω(t)=2π(7.5-2.5cos(t))t=[0,2π]

(12)

正弦掃頻信號模型如式(13)所示，其中包含1，2，4階三個分量[15]。

(13)

其中：η(t)為高斯白噪聲；信噪比為-3 dB；采樣率為100 Hz。

轉頻呈正弦規(guī)律變化，仿真模型中包含轉頻的一倍頻，即一階分量，幅值為1；二階分量，幅值為0.6；四階分量，幅值為0.4。STFT時頻譜如圖5所示。

圖5 正弦掃頻信號STFT譜圖Fig.5 STFT for sinusoidal sweep frequency signal

現(xiàn)采用基于Viterbi算法的IFE方法提取第一個分量，結果如圖6所示。點線為采用Viterbi算法估計的瞬時頻率連線；粗實線為將瞬時頻率點采用5次多項式擬合得到的頻率曲線；短劃線為理論頻率曲線?？梢钥吹近c線相對于理論值在局部存在一定誤差，出現(xiàn)突變點，但轉速通常為平穩(wěn)連續(xù)變化，因此這里采用低階多項式進行擬合，得到比較平滑的估計曲線。擬合曲線相對于理論曲線的百分比誤差為ξ=3.45%，可以得出，基于Viterbi算法的IFE方法可以精確估計瞬時轉頻，具有較好的抗噪性和抗臨近階次的能力。

圖6 正弦掃頻信號IFE結果Fig.6 IFE results for sinusoidal sweep frequency signal

3.3 實驗信號

采用QPZZ-Ⅱ型旋轉機械振動及故障模擬實驗平臺進行變轉速振動測試試驗，實驗臺如圖7所示。

圖7 旋轉機械振動及故障模擬實驗臺Fig.7 Test rig of rotating machinery vibration and fault simulation

實驗采用NI PXIe 4496數(shù)據(jù)采集設備，采用PCB 355B03型壓電加速度傳感器測量振動信號，采樣頻率為25.6 kHz；采用PCB LaserTach ICP型激光轉速傳感器測量參考軸轉速，采樣頻率為1 kHz。根據(jù)轉速脈沖信號采用5點公式法[16]計算參考軸的瞬時轉頻，并將該方法測得的轉頻作為參考值。

實驗中軸承節(jié)圓直徑D=38.5 mm，滾子直徑d=7.2 mm，滾子個數(shù)Z=13，接觸角α=0。根據(jù)外圈故障特征頻率計算公式，得外圈故障特征階次為

Obpo=Z[1-d/Dcos(α)]/2≈5.284

(14)

模擬轉速劇烈波動工況，原始振動信號如圖8所示。根據(jù)軸承故障診斷實驗平臺振動信號的特點，采用共振解調方法進行故障特征提取。采集信號通過包絡解調就可以得出隨轉速變化的低頻幅值包絡信號。

圖8 原始時域振動信號Fig.8 Original vibration signal in time domain

為提取低頻區(qū)段內的轉速信息，對信號進行低通濾波和降采樣處理，低通濾波截止頻率設定為1 500 Hz，降采樣倍數(shù)為5，降采樣后的采樣率為5 120 Hz。然后對降采樣后的數(shù)據(jù)進行Hilbert包絡并去趨勢項得預處理振動信號。預處理振動信號STFT譜圖如圖9所示。

圖9 振動信號STFT時頻譜圖Fig.9 STFT for the rolling bearing vibration signal

根據(jù)軸承故障特征頻率計算公式可看出，故障特征階次一倍頻的峰值在時頻譜圖中最突出，因此將其作為搜索目標。

圖10 Viterbi法、峰值搜索法、轉速計測得結果對比圖Fig.10 Comparison of the results from Viterbi-IFE、Peaksearch-IFE and the tachometer

對圖9進行基于Viterbi算法的轉頻估計，得到估計結果，如圖10所示。圖中實線為采用激光轉速計測量轉速脈沖計算得到的參考軸瞬時頻率曲線，星點線為基于Viterbi算法所估計出的轉頻曲線，點虛線為采用峰值搜索算法估計的轉速曲線。可以看出，Viterbi-IFE和轉速計測得轉速曲線基本一致，而峰值搜索在局部出現(xiàn)較大偏差。求估計瞬時頻率值相對于參考值的百分比誤差，根據(jù)公式(9)計算，峰值搜索誤差為8.5%，Viterbi-IFE估計誤差僅為0.71%。

分別采用改進前與改進后的Viterbi算法對實驗信號STFT譜圖進行瞬時頻率估計，其代價函數(shù)迭代過程所用時間如表1所示。計算時統(tǒng)一采用MATLAB編程，硬件為DELL筆記本電腦，4核Core i5處理器，主頻1.80 GHz。在表1中，將式(10)所示的線性掃頻信號和式(13)所示的正弦掃頻信號同時考慮在內，分別采用改進前與改進后的Viterbi算法對該兩組信號STFT譜圖進行瞬時頻率估計。從表1中可知，在仿真信號和實驗信號驗證中，使用相同計算資源的前提下，改進后的Viterbi算法計算效率得到明顯提升。尤其是實驗信號，無判定準則的Viterbi-IFE用時1.45 s，有判定準則的改進Viterbi-IFE用時降至0.52 s。整個轉速變化過程持續(xù)時間為38 s，而進行轉頻估計的迭代計算過程用時0.52 s，從響應的及時性分析，該算法可以滿足實時分析的要求。

表1 改進前后迭代計算所用時間

Tab.1 The time used for the iterative calculation before and after the improvement s

4 結束語

基于振動信號的瞬時轉速估計使無轉速計的軸承狀態(tài)監(jiān)測和故障診斷成為可能，其中基于振動信號STFT譜圖的局部峰值搜索IFE方法是比較早的方法，但其抗噪性和抗臨近階次的能力有限。筆者研究了基于Viterbi算法的IFE方法，并分析了代價函數(shù)的迭代過程和次數(shù)，在此基礎上，提出了代價函數(shù)迭代循環(huán)停止的新型判定準則，該準則可以實現(xiàn)時頻平面最優(yōu)路徑的快速搜索，提高了瞬時頻率的估算精度和計算速度。通過仿真信號和實驗數(shù)據(jù)進行驗證，具有判定準則的改進后算法較改進前算法的效率明顯提高，為無轉速計的軸承狀態(tài)監(jiān)測和故障診斷奠定了基礎。

[1] 趙曉平，張令彌，郭勤濤. 旋轉機械階比跟蹤技術研究進展綜述[J]. 地震工程與工程振動, 2008, 28(6): 213-219.

Zhao Xiaoping, Zhang Lingmi, Guo Qintao. Advances and trends in rotational machine order tracking methodology[J]. Journal of Earthquake Engineering and Engineering Vibration, 2008, 28(6): 213-219. (in Chinese)

[2] 陳平. 信號瞬時頻率的估計方法及其應用[D].濟南：山東大學, 2007.

[3] 郭瑜，秦樹人，湯寶平，等. 基于瞬時頻率估計的旋轉機械階比跟蹤[J]. 機械工程學報，2003, 39(3): 32-36.

Guo Yu, Qin Shuren, Tang Baoping, et al. Order tracking of rotating machinery based on instantaneous frequency estimation[J]. Journal of Mechanical Engineering, 2003, 39(3): 32-36. (in Chinese)

[4] 胡愛軍，朱瑜. 基于改進峰值搜索法的旋轉機械瞬時頻率估計[J]. 振動與沖擊，2013, 32(7): 113-117.

Hu Aijun, Zhu Yu. Instantaneous frequency estimation of a rotating machinery based on an improved peak search method[J]. Journal of Vibration and Shock, 2013, 32(7): 113-117. (in Chinese)

[5] 楊志堅，丁康，梁茜. 基于頻譜校正理論的階比跟蹤分析[J]. 機械工程學報，2009, 45(12): 41-45.

Yang Zhijian, Ding Kang, Liang Qian. Novel method of order tracking analysis based on spectrum correction[J]. Journal of Mechanical Engineering, 2009, 45(12): 41-45. (in Chinese)

[6] 趙曉平. 旋轉機械階比分析研究與軟件實現(xiàn)[D]. 南京：南京航空航天大學, 2008.

[7] 趙曉平，趙秀莉，侯榮濤，等. 一種新的旋轉機械升降速階段振動信號瞬時頻率的估計算法[J]. 機械工程學報, 2011, 47(7): 104-108.

Zhao Xiaoping, Zhao Xiuli, Hou Rongtao, et al. A new method for instantaneous frequency estimation of run-up or run-down vibration signal for rotating machinery[J]. Journal of Mechanical Engineering, 2011, 47(7): 104-108. (in Chinese)

[8] Leclere Q, André H, Antoni J. A multi-order probabilistic approach for Instantaneous Angular Speed tracking debriefing of the CMMNO′14 diagnosis contest[J]. Mechanical Systems and Signal Processing, 2016(81): 375-386.

[9] Combet F, Zimroz R. A new method for the estimation of the instantaneous speed relative fluctuation in a vibration signal based on the short time scale transform[J]. Mechanical Systems and Signal Processing, 2009, 23(4): 1382-1397.

[10] Urbanek J, Barszcz T, Antoni J. A two-step procedure for estimation of instantaneous rotational speed with large fluctuations[J]. Mechanical Systems and Signal Processing, 2013, 38(1): 96-102.

[11] Viterbi A. Error bounds for convolutional codes and an asymptotically optimum decoding algorithm[J]. IEEE Transactions on Information Theory, 1967,13(2):260-269.

[12] Martell A. Edge detection using heuristic search methods[J]. Comput, Graphic,Image Process, 1982,1(2): 69-182.

[13] 陳源，江修富. Viterbi算法的關鍵問題研究及DSP實現(xiàn)[J]. 裝備指揮技術學院學報, 2005,16(5):77-81.

Chen Yuan, Jiang Xiufu. Study on the main problems of viterbi algor ithm and the realization of DSP[J]. Journal of the Academy of Equipment Command & Technology, 2005,16(5):77-81. (in Chinese)

[14] 胡旭娟，梁紅，蒯繼武. 改進WVD的調頻信號瞬時頻率估計算法[J]. 彈箭與制導學報, 2007, 27(5): 329-334.

Hu Xujuan, Liang Hong, Kuai Jiwu. FM Signal IF estimation algorithm based on improved WVD[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2007, 27(5): 329-334. (in Chinese)

[15] 戴功偉. 基于自適應轉頻跟蹤濾波的旋轉機械階次析研究[D]. 重慶：重慶大學，2013.

[16] 魏玉果. 基于自適應計算階次跟蹤的旋轉機械階次分析系統(tǒng)[D].重慶：重慶大學, 2007.

國家自然科學基金資助項目(11572206，11472179，U1534204，11372199)；河北省自然科學基金資助項目(A2015210005，A2016210099)；河北省人才工程培養(yǎng)經(jīng)費資助科研項目(A2016002036)

2016-09-13；

2017-01-18

TH165.3

劉永強，男，1983年12月生，副教授、博士生導師。主要研究方向為車輛系統(tǒng)動力學、機車狀態(tài)監(jiān)測與故障診斷。曾發(fā)表《一種自適應共振解調方法及其在滾動軸承早期故障診斷中的應用》(《振動工程學報》 2016年第29卷第2期)等論文。

E-mail：liuyq@stdu.edu.cn