龔思琪，李舜酩，庾天翼

（南京航空航天大學(xué) 能源與動力學(xué)院，南京 210016）

變轉(zhuǎn)速是旋轉(zhuǎn)機械設(shè)備的常見復(fù)雜工況，其核心旋轉(zhuǎn)零部件在變轉(zhuǎn)速工況下的故障振動信號具有明顯的非平穩(wěn)特點［1］?；谡駝有盘柼幚淼墓收显\斷方法對設(shè)備狀態(tài)監(jiān)測和故障診斷都至關(guān)重要［2］。

傳統(tǒng)的基于振動信號處理的故障診斷方法大多以不同狀態(tài)下的時域、頻域及時頻域信息作為參考，通過對它們的統(tǒng)計特征指標(biāo)對比分析設(shè)備當(dāng)下的健康狀態(tài)及潛在的故障特征［3－4］。在變轉(zhuǎn)速下的故障振動信號的信號表征與故障模式之間的對應(yīng)關(guān)系更復(fù)雜，導(dǎo)致傳統(tǒng)的傅里葉變換并不能完全揭示變轉(zhuǎn)速信號的局部細(xì)節(jié)特征。在近20年來，隨著非平穩(wěn)信號處理理論的發(fā)展和電子硬件計算技術(shù)的快速成長［5］，對變轉(zhuǎn)速下機械設(shè)備動態(tài)信號處理及故障診斷方法開展了廣泛的研究。研究者對小波變換［6］和高精度的頻譜與時頻分析方法［7－8］等故障診斷技術(shù)進(jìn)行了廣泛研究，并綜合運用統(tǒng)計物理學(xué)方面的信號特征提取方法［9］、盲源分離方法［10］、變分模式分解［11］等方法提高診斷準(zhǔn)確率。瞬時頻率變化脊線是振動信號的重要特征之一。通常利用峰值檢測法提取時頻脊線［12］，文獻(xiàn)［13］提出了非線性擠壓時頻變換方法提取旋轉(zhuǎn)機械故障振動信號中的多分量非穩(wěn)態(tài)信號。但在受到強噪聲干擾的變轉(zhuǎn)速工況下的振動信號時頻脊線提取和瞬時轉(zhuǎn)速提取時，準(zhǔn)確率較低，辨識效果較差［13］，且對于含有多個頻率分量的信號，脊線難以辨別。文獻(xiàn)［14］提出了一種雙向搜索時頻脊融合方法，提取時頻變工況行星齒輪箱高速軸的瞬時轉(zhuǎn)速，可靠性高，但程序編寫困難，工作量大，而且需要豐富的理論知識。

針對噪聲干擾下的變轉(zhuǎn)速旋轉(zhuǎn)機械振動信號時頻脊線提取和沖擊特征識別問題，本文在利用短時傅里葉變換得到轉(zhuǎn)速波動信號的時頻分布圖；再將其沿時間軸分割，利用提取中心頻率的方法逐段提取時頻脊線；然后，利用卡爾曼濾波技術(shù)提取整段信號的最優(yōu)脊線。通過分析脊線特點，可以準(zhǔn)確有效地判斷轉(zhuǎn)子故障發(fā)生的時間。進(jìn)一步，利用本文所提方法提取含有多個頻率分量信號的時頻脊線特征；最后，利用變轉(zhuǎn)速故障軸承的振動信號驗證本文所提方法的有效性和準(zhǔn)確性。結(jié)果表明：所提方法的時頻脊線提取準(zhǔn)確率高，操作簡單，易于實現(xiàn)，對轉(zhuǎn)速波動信號的特征提取有重要的參考價值。

1 轉(zhuǎn)軸故障振動信號的頻譜分析

1．1 實驗信號采集

本研究所使用的試驗臺如圖1所示，由電機驅(qū)動，在600～1 500 r／min的轉(zhuǎn)速下工作，采用油脂潤滑。驅(qū)動電機的額定功率為0．75 kW，額定電壓為380 V，額定轉(zhuǎn)矩是5．0 kN·m。

圖1 傳動軸變轉(zhuǎn)速故障振動試驗臺及傳感器放置位置

為了測量轉(zhuǎn)軸在變轉(zhuǎn)速下的振動信號，在轉(zhuǎn)軸驅(qū)動端和輸出端各安裝一個振動三向加速度傳感器采集傳動變轉(zhuǎn)速振動信號，傳感器安裝位置如圖1（a）、（b）所示。本次變轉(zhuǎn)速以勻加速方式實現(xiàn)，使電機以600 r／min的轉(zhuǎn)速開始工作，勻加速至1 500 r／min。測試所用的振動傳感器型號和基本參數(shù)如表1所示。

表1 振動傳感器的基本參數(shù)

在試驗過程中，采用LMS振動信號采集系統(tǒng)采集被測轉(zhuǎn)軸的故障振動信號。根據(jù)轉(zhuǎn)速和采樣頻率可知信號一個周期，即試驗轉(zhuǎn)軸旋轉(zhuǎn)一周所采得的數(shù)據(jù)點個數(shù)為2 048～5 120。

在測試過程中，選擇5 s內(nèi)的振動信號作為分析目標(biāo)，測試得到的變轉(zhuǎn)速轉(zhuǎn)軸的不平衡故障振動信號波形如圖2所示。

圖2 轉(zhuǎn)速波動信號時域波形

1．2 信號的頻譜特點

根據(jù)式（1）所示的傅里葉變換（fourier trans form，F(xiàn)T）公式，通過改進(jìn)可將信號通過快速傅里葉變換，得到圖3所示的轉(zhuǎn)速波動信號的頻譜圖。

從圖3轉(zhuǎn)速波動信號的頻譜圖可以看到：信號主要集中在小于1 000 Hz的范圍內(nèi)，并伴有噪聲干擾，可以看到頻率在0～400 Hz出現(xiàn)了最大幅值，且有多個頻率分量。

圖3 轉(zhuǎn)速波動信號的頻譜圖

2 信號的時頻分析

根據(jù)式（2）所示的短時傅里葉變換（short time fourier transform，STFT）的表達(dá)式，可得到時頻分布圖。

式中：STFT（t，w）表示STFT分析結(jié)果；f（t）為待分析信號；g（t）表示STFT窗函數(shù)。

窗口函數(shù)g（t）的選擇決定了分辨率的高低。根據(jù)不同窗函數(shù)的特點可知：具有帶通特性的漢寧窗或漢明窗雖然對信號濾波有更好的效果［15］，但是會使能量幾乎集中在主瓣。而高斯窗函數(shù)可以取得時間和頻率分辨率的最佳折衷［16］。因此，選取高斯窗函數(shù)對信號進(jìn)行STFT。

對變轉(zhuǎn)速故障振動信號進(jìn)行STFT分析，得到如圖4所示的二維時頻圖。可以看出：頻率在0～1 000 Hz和4 000 Hz附近的能量高，頻率隨時間基本保持不變，它們是由轉(zhuǎn)子的固有振動模態(tài)引起的。但在6～15 s的時間范圍內(nèi)，如圖4中的頻率線2和3以及它們附近的多個頻率出現(xiàn)了隨時間變化而不斷上升的趨勢。這表示隨著采樣時間的增加，轉(zhuǎn)軸的轉(zhuǎn)動頻率也逐漸增加，這是由于轉(zhuǎn)子在6～8 s產(chǎn)生波動導(dǎo)致的頻率沿采樣時間的增加而升高。

圖4 轉(zhuǎn)速波動信號的二維時頻圖

由于轉(zhuǎn)子轉(zhuǎn)速在600～1 500 r／min波動，轉(zhuǎn)子轉(zhuǎn)頻最高在25 Hz附近，通過圖3的轉(zhuǎn)速波動信號的頻譜圖可以判斷轉(zhuǎn)速波動信號的頻率主要集中在0～400 Hz。從圖5可以看到：頻率不隨時間變化的線4，以及在時間大于8 s時，頻率隨時間變化的線5、線6。其中，線4的頻率是56 Hz，大約是轉(zhuǎn)子頻率的3～4倍。

圖5 時頻圖4的局部放大圖

3 信號的脊線提取

對信號進(jìn)行STFT后，一些與機械系統(tǒng)零部件的動力學(xué)特性相關(guān)的曲線會清晰或隱約地出現(xiàn)，這些曲線在三維空間內(nèi)的分布形態(tài)與“山脊”很像，是脊線標(biāo)志或沿脊頂延伸的線。因此，在時頻圖中把三維空間內(nèi)沿功率最高點延伸的線稱之為時頻脊線［17］，表現(xiàn)在時頻面上，是明亮程度和趨勢各異的線條，它能夠表征信號頻率變化的各種模式［18］。

信號在時域和頻域的能量可表示為：

能量分布還滿足邊緣性質(zhì)：

從時頻表示上直接提取的時頻脊線受噪聲影響大，脊線特征不明顯，從峰峰值上得到的脊點，噪聲難以消除，導(dǎo)致所需提取的脊線淹沒在噪聲的峰值點中，如圖6所示。

圖6 轉(zhuǎn)速波動信號提取時頻脊線

從時頻能量最高點上提取信號的中心主頻率，遍歷時間軸，會得到圖7中的線1、2、3。在0～6 s的時間段內(nèi)，轉(zhuǎn)速波動信號的頻率主要在56 Hz附近有小范圍的波動；在6～9 s時間段內(nèi)，中心頻率與時間的關(guān)系變得很復(fù)雜；在9～12 s時間段內(nèi)，中心頻率集中在二階頻率分量處有上升趨勢；在12～20 s的時間段內(nèi)，中心主頻率在25 Hz附近隨著時間的增加而升高；在12 s處頻率有突變，如圖7中區(qū)域3所示。導(dǎo)致這種突變產(chǎn)生的主要原因是在6～8 s時間段內(nèi)產(chǎn)生故障后，信號的頻率分量與能量的對應(yīng)關(guān)系模糊，使得部分頻率分量的能量突然間提高，而其他頻率分量的能量降低。

圖7 轉(zhuǎn)速波動信號脊線圖

時頻面內(nèi)脊線附近局部區(qū)域內(nèi)的時頻點能量突出，具有相對較高的信噪比，機械系統(tǒng)運行過程中狀態(tài)一般不會出現(xiàn)突變，出現(xiàn)的時頻脊線具有連續(xù)平滑的特點。因此，直接通過能量最高點來提取時頻脊線的方法失效。但該方法表示出了該轉(zhuǎn)速波動信號的能量變化，為后面沿時間域分割信號的分割點提供了依據(jù)。

根據(jù)脊線過渡平滑的特點，以及圖7直接提取的中心頻率的劇烈變化，將提取區(qū)域按時間劃分為3個部分：0～6 s，6～9 s，9～20 s。目標(biāo)頻率范圍設(shè)置為（f－Δf，f＋Δf），Δf＝5 Hz由于轉(zhuǎn)子的轉(zhuǎn)速為600～1 500 r／min，將圖7中線2的12～20 s所示的脊線即表示一階頻率作為即將要提取的目標(biāo)脊線，通過對時間軸分段提取，將提取的頻率閾值設(shè)置在接近一階中心頻率附近，再根據(jù)能量最高點提取頻率范圍內(nèi)的曲線。

圖8 轉(zhuǎn)速波動信號目標(biāo)脊線

由圖7的信號脊線圖可以看出：56 Hz的頻率大概是轉(zhuǎn)速信號的3～4倍頻。由此可推算出：一階頻率在14 Hz附近，在0～6 s將目標(biāo)頻率f設(shè)置為14 Hz，6～9 s頻率變化范圍大，且初步確定是轉(zhuǎn)速增加的，故將目標(biāo)頻率設(shè)置為f＝18 Hz。9～20 s也是一個轉(zhuǎn)速上升的過程，但相比第2段，其過渡較為平穩(wěn)，最高轉(zhuǎn)速1 500 r／min，頻率為25 Hz，將其目標(biāo)頻率設(shè)置為f＝25 Hz，得到如圖8線1所示的脊線。

從圖8中的線1可以看到：在0～5 s目標(biāo)脊線在這一時間段的頻率隨時間波動比較劇烈，結(jié)合圖7中該段時間內(nèi)的能量集中在56 Hz附近，因此在圖8中該段時間內(nèi)在該處能量比較分散會引起較大的波動。目標(biāo)脊線在5～10 s之間的過渡十分復(fù)雜，時間在10～20 s之間的目標(biāo)脊線波動較小，但頻率隨時間增加的趨勢也更為明顯。

4 最優(yōu)脊線提取

卡爾曼濾波器是一種最優(yōu)化自回歸數(shù)據(jù)處理算法，對于解決很大部分問題都會是最優(yōu)、效率最高甚至是最有用的［20］方法。因此，本文采用卡爾曼濾波在時頻面內(nèi)提取目標(biāo)脊線的最優(yōu)脊線［21］。

Kalman濾波器的系統(tǒng)測量值、狀態(tài)運動方程和觀測方程分別為［22］：

式中：X（k）是k時刻的系統(tǒng)狀態(tài)；U（k）是k時刻對系統(tǒng)的控制量；A是狀態(tài)轉(zhuǎn)移矩陣；Z（k）是k時刻的觀測值；H是觀測矩陣；W（k）和V（k）表示上述兩個方程中的噪聲。

在實際問題中，通常假設(shè)A（k）＝A，H（k）＝H。設(shè)W和V為噪聲（white gaussian noise），協(xié)方差分別為Q，R（這里假設(shè)它們不隨系統(tǒng)狀態(tài)變化而變化）。首先，利用狀態(tài)方程預(yù)測下一狀態(tài)的系統(tǒng)。假設(shè)當(dāng)前時刻為k，根據(jù)系統(tǒng)的模型，可以基于系統(tǒng)的上一狀態(tài)而預(yù)測出現(xiàn)在狀態(tài)：

式中：X（k｜k－1）為利用上一狀態(tài)的預(yù)測結(jié)果；X（k－1｜k－1）為上一狀態(tài)的最優(yōu)結(jié)果；U（k）為現(xiàn)在狀態(tài)的控制量，若沒有控制量，它可以為0。

對應(yīng)于X（k｜k－1）的協(xié)方差，即上面預(yù)測值相對于真實狀態(tài)的誤差方差陣更新公式為：

式中：P（k｜K－1）是X（k｜k－1）對應(yīng)的協(xié)方差；P（k－1｜k－1）是X（k－1｜k－1）對應(yīng)的協(xié)方差；A′表示A的轉(zhuǎn)置矩陣；Q是系統(tǒng)過程的協(xié)方差。

當(dāng)?shù)趉時刻的觀測值Z（k）到達(dá)以后，利用預(yù)測值中沒有的信息去修正第k個時刻的狀態(tài)預(yù)測值。結(jié)合預(yù)測值和測量值，可以得到狀態(tài)（k）時最優(yōu)化估算值X（k｜k）：

式中，Kg為卡爾曼增益（Kalman Gain），可通過式（13）估計：

為了使Kalman濾波器運行直到系統(tǒng)過程結(jié)束，需要更新k狀態(tài)下X（k｜k）的協(xié)方差：

式中，I為單位矩陣，對于單模型單測量，I＝1。當(dāng)系統(tǒng)進(jìn)入k＋1狀態(tài)時，P（k｜k）就是式（11）的P（k－1｜k－1）。

通過以上卡爾曼濾波器的5個核心公式不斷進(jìn)行迭代優(yōu)化運算，可以完整實現(xiàn)Kalman濾波。將原始脊線的初始坐標(biāo)作為卡爾曼濾波的初始狀態(tài)輸入，將P的初值設(shè)為1，保證卡爾曼濾波能夠穩(wěn)定運行。通過多次試驗可得協(xié)方差Q與R分別為10－8，10－5時，Kalman濾波的濾波效果最佳。根據(jù)脊線平滑且表示信號的頻率隨時間變化趨勢的特點，得到圖9所示的轉(zhuǎn)速波動信號的最優(yōu)目標(biāo)脊線。圖9中藍(lán)色曲線表示經(jīng)過沿時間軸分割后提取的原始脊線，紅色曲線表示經(jīng)過卡爾曼濾波技術(shù)后得到的最優(yōu)脊線，橙色曲線表示兩曲線之間的誤差。為了使脊線信息更加明顯，將縱坐標(biāo)范圍縮小至0～45 Hz。相比原始脊線，最優(yōu)脊線更加平滑，且其表達(dá)的頻率隨時間變化的趨勢沒有改變。在0～4 s時間段內(nèi)，可以看到經(jīng)卡爾曼濾波器提取的脊線頻率更加集中，這也更加真實地反映了在這段時間內(nèi)轉(zhuǎn)子正常運轉(zhuǎn)。而誤差曲線主要集中在0 Hz附近，說明最優(yōu)脊線與實際信號誤差小，最優(yōu)脊線能反映實際信號波動情況的真實信息。

圖9 轉(zhuǎn)速波動信號最優(yōu)脊線

5 結(jié)論

1）在STFT時頻分析的基礎(chǔ)上，結(jié)合能量最高點提取中心主頻率的方法，提取轉(zhuǎn)速波動信號的時頻脊線，可以觀察到信號頻率在不同時刻發(fā)生的變化，并可以判斷轉(zhuǎn)子系統(tǒng)是否發(fā)生故障以及故障發(fā)生的時間。相比直接提取時頻脊，其所受噪聲影響更小。

2）通過對不平穩(wěn)信號的時頻脊線進(jìn)行卡爾曼濾波，可以有效地找到目標(biāo)脊線的最優(yōu)解，得到表達(dá)信號頻率隨時間變化的趨勢的平滑的時頻脊線。

3）通過對轉(zhuǎn)軸試驗臺實測信號變換結(jié)果的觀察分析可知：從故障振動信號的時頻譜圖提取時頻脊線，并利用卡爾曼濾波估計最優(yōu)脊線，可以有效地濾除背景噪聲，剔除干擾分量，得到平滑且清晰的時頻脊線。