轉(zhuǎn)軸轉(zhuǎn)速波動(dòng)信號(hào)的時(shí)頻脊線提取方法

龔思琪，李舜酩，庾天翼

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

變轉(zhuǎn)速是旋轉(zhuǎn)機(jī)械設(shè)備的常見(jiàn)復(fù)雜工況，其核心旋轉(zhuǎn)零部件在變轉(zhuǎn)速工況下的故障振動(dòng)信號(hào)具有明顯的非平穩(wěn)特點(diǎn)［1］?；谡駝?dòng)信號(hào)處理的故障診斷方法對(duì)設(shè)備狀態(tài)監(jiān)測(cè)和故障診斷都至關(guān)重要［2］。

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

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

1 轉(zhuǎn)軸故障振動(dòng)信號(hào)的頻譜分析

1．1 實(shí)驗(yàn)信號(hào)采集

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

圖1 傳動(dòng)軸變轉(zhuǎn)速故障振動(dòng)試驗(yàn)臺(tái)及傳感器放置位置

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

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

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

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

圖2 轉(zhuǎn)速波動(dòng)信號(hào)時(shí)域波形

1．2 信號(hào)的頻譜特點(diǎn)

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

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

圖3 轉(zhuǎn)速波動(dòng)信號(hào)的頻譜圖

2 信號(hào)的時(shí)頻分析

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

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

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

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

圖4 轉(zhuǎn)速波動(dòng)信號(hào)的二維時(shí)頻圖

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

圖5 時(shí)頻圖4的局部放大圖

3 信號(hào)的脊線提取

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

信號(hào)在時(shí)域和頻域的能量可表示為：

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

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

圖6 轉(zhuǎn)速波動(dòng)信號(hào)提取時(shí)頻脊線

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

圖7 轉(zhuǎn)速波動(dòng)信號(hào)脊線圖

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

根據(jù)脊線過(guò)渡平滑的特點(diǎn)，以及圖7直接提取的中心頻率的劇烈變化，將提取區(qū)域按時(shí)間劃分為3個(gè)部分：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)脊線，通過(guò)對(duì)時(shí)間軸分段提取，將提取的頻率閾值設(shè)置在接近一階中心頻率附近，再根據(jù)能量最高點(diǎn)提取頻率范圍內(nèi)的曲線。

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

由圖7的信號(hào)脊線圖可以看出：56 Hz的頻率大概是轉(zhuǎn)速信號(hào)的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也是一個(gè)轉(zhuǎn)速上升的過(guò)程，但相比第2段，其過(guò)渡較為平穩(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)脊線在這一時(shí)間段的頻率隨時(shí)間波動(dòng)比較劇烈，結(jié)合圖7中該段時(shí)間內(nèi)的能量集中在56 Hz附近，因此在圖8中該段時(shí)間內(nèi)在該處能量比較分散會(huì)引起較大的波動(dòng)。目標(biāo)脊線在5～10 s之間的過(guò)渡十分復(fù)雜，時(shí)間在10～20 s之間的目標(biāo)脊線波動(dòng)較小，但頻率隨時(shí)間增加的趨勢(shì)也更為明顯。

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

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

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

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

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

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

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

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

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

式中，Kg為卡爾曼增益（Kalman Gain），可通過(guò)式（13）估計(jì)：

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

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

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

圖9 轉(zhuǎn)速波動(dòng)信號(hào)最優(yōu)脊線

5 結(jié)論

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

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

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