趙結(jié)昂，包俊義，金嵩，方琦，商高亮，劉俊

(1.金華市特種設(shè)備檢測中心，浙江金華 321000；2.武漢科技大學機械自動化學院，湖北武漢 430081)

0 前言

滾動軸承在機械設(shè)備傳動系統(tǒng)中具有廣泛應用，由于其工作環(huán)境特點以高負荷、高速、變工況、強噪聲干擾為主，故其故障率較高。因此，對滾動軸承的潛在故障以及故障類型進行挖掘和判斷，是指導維修工作和提高維修效率的重要前提，是保障機械設(shè)備安全性和可靠性的必要途徑[1]。

變轉(zhuǎn)速工況下的非平穩(wěn)信號是工程實際中最常見的信號，對它準確分析是實現(xiàn)機械設(shè)備故障診斷和健康管理的前提。但非平穩(wěn)信號與平穩(wěn)信號不同，采用傳統(tǒng)的頻域分析方法并不能直接找到故障特征頻率，這是因為信號的頻率具有時變性，使得特征頻率的能量分散到了頻率軸上，難以像平穩(wěn)信號的頻譜一樣突出某個特征頻率的能量[2-3]。時頻分析(Time-Frequency Analysis，TFA)是進行變轉(zhuǎn)速滾動軸承故障診斷的有效方法，使用時頻域聯(lián)合分布描述信號的瞬態(tài)特征，并通過瞬時頻率估計來表征信號特征頻率隨時間變化的趨勢[4]。同步壓縮變換(Synchrosqueezing Transform，SST)是由DAUBECHIES等[5]提出的一種經(jīng)典的時頻分析方法，在連續(xù)小波變換(Continuous Wavelet Transform，CWT)的基礎(chǔ)上估計瞬時頻率(Instantaneous Frequency，IF)，然后在時頻面上將能量沿頻率方向進行重排，其基本思想是時頻能量在時間方向上同步，在頻率方向上壓縮[6-7]。

SST能夠在一定程度上揭示非平穩(wěn)信號復雜的時頻特性，因此眾多學者對它進行研究。最初的SST是基于CWT進行推導的，但是由于短時傅里葉變換(Short Time Fourier Transform，STFT)具有突出的時頻物理含義以及其簡單直觀的公式，之后有學者提出了基于STFT的FSST(STFT based Synchrosqueezing Transform，F(xiàn)SST)，并將它應用于變轉(zhuǎn)速滾動軸承故障診斷[8]。此后，BEHERA等[9]在FSST中引入群延遲的思想，并對STFT的計算結(jié)果展開二階瞬時頻率估計計算，由此提出了二階同步壓縮變換(Second-order Synchrosqueezing Transform，SST2)，在調(diào)制多分量信號的分析中驗證了其有效性。FSST計算的核心在于實現(xiàn)對IF的準確估計，然后利用同步壓縮算子將時頻能量往IF曲線上進行壓縮。以上方法均基于STFT進行后處理，然而STFT在進行瞬時頻率估計時，由于調(diào)制項的存在使得估計的瞬時頻率產(chǎn)生偏差，且瞬時頻率在時頻譜上出現(xiàn)時頻擴散現(xiàn)象[10-11]，這不利于后續(xù)的同步壓縮操作。

本文作者基于經(jīng)典的STFT理論進行改進研究，針對STFT公式中存在的調(diào)制現(xiàn)象，提出一種解調(diào)短時傅里葉變換(Demodulation Short Time Fourier Transform，DSTFT)方法，在STFT的基礎(chǔ)上引入時變解調(diào)算子解決由調(diào)制成分帶來的時頻能量發(fā)散、時頻聚集性不高的問題，并通過迭代算法獲得接近理想的IF估計。在DSTFT的基礎(chǔ)上，提出解調(diào)同步壓縮變換(Demodulation Synchrosqueezing Transform，DSST)，該方法作為一種時頻分析后處理的手段用以改善傳統(tǒng)SST瞬時頻率估計不夠準確，從而導致壓縮后的時頻面出現(xiàn)能量發(fā)散的問題。因此，利用該方法進行時頻分析后處理，能夠獲得能量集中的時頻表達。分別通過仿真信號分析、實測試驗臺滾動軸承故障信號分析對該方法進行驗證，試驗結(jié)果進一步驗證了DSST的有效性。

1 理論描述

1.1 基于STFT的同步壓縮變換

定義一個調(diào)幅調(diào)頻(AM-FM)信號為

s(t)=A(t)eiΦ(t)

(1)

其中：A(t)和Φ(t)分別對應調(diào)幅調(diào)頻信號的瞬時幅值和瞬時相位，φ(t)=Φ′(t)可以視為信號的理想瞬時頻率。對于信號s(t)，其STFT寫為

(2)

式中：g(t)為窗函數(shù)，g*(t)為g(t)的復共軛。

基于STFT的FSST算法可以用下面的公式表示：

(3)

(4)

(5)

1.2 解調(diào)同步壓縮變換

根據(jù)Ville理論[12]，公式(1)所示的調(diào)幅-調(diào)頻信號同樣可以寫成：

(6)

式(6)中信號的瞬時頻率可以表示為fIF=φ(t)。在一個很短的時間間隔ξ內(nèi)，信號s(t)的瞬時頻率可以近似地視為線性方程，將φ(t)進行一階泰勒展開可以得到：

φ(ξ)=φ(t)+φ(t)(ξ-t)

(7)

則信號s(t)的STFT可以展開為

(8)

通過上式可以看出，STFT由于eiφ(t)(ξ-t)2/2項的存在，信號發(fā)生了調(diào)制，使得STFT的能量發(fā)散。為此，引入時變解調(diào)算子e-ic(t)(ξ-t)2/2消去eiφ(t)(ξ-t)2/2的影響。引入時變解調(diào)算子后的DSTFT可以寫成：

(9)

其中：c(t)是瞬時頻率的一階導數(shù)，由于c(t)的值未知，文中采用遞歸的策略，通過算法多步逼近，從而在精度允許的范圍內(nèi)得到其近似值。c(t)的計算流程如下：

(1)初始計算時設(shè)c(t)=0，進行DSTFT,此時的DSTFT相當于STFT。

(2)對信號進行檢索得到對應時刻t的峰值數(shù)據(jù)P(t)，再對峰值數(shù)據(jù)P(t)曲線進行多項式擬合，得到信號的IF估計曲線I(t)，然后把瞬時頻率估計曲線I(t)的一階導數(shù)賦值給下一次迭代時的c(t)，即c(t)=I′(t)。

(3)重復步驟(2)依次進行迭代運算。

(4)迭代的終止條件為第n次與第n-1次迭代時的IF曲線c(t)之間的平均誤差ε小于一個預定的閾值μ，即：

(10)

其中：In(t)和In-1(t)分別為進行第n次與第n-1次計算所得到的瞬時頻率估計值；Lt為In(t)曲線的時間長度；μ一般選擇為0.05。通過上面幾個步驟就能消除調(diào)制成分的影響，較為精確地得到IF估計值。

基于DSTFT的IF估計可以表示為

(11)

因此，本文作者提出的解調(diào)同步壓縮變換(DSST)的計算公式可以描述為

(12)

(13)

1.3 基于DSST的變轉(zhuǎn)速滾動軸承故障診斷算法

滾動軸承可能存在著外圈、內(nèi)圈等多種類型的故障，每一種故障在頻域中所展現(xiàn)出來的特征頻率是不同的，其大小與滾動軸承的結(jié)構(gòu)參數(shù)和轉(zhuǎn)頻相關(guān)。在變轉(zhuǎn)速工況下，由于轉(zhuǎn)頻的變化，故障特征頻率在時頻平面中呈現(xiàn)出的是一條曲線。實現(xiàn)變轉(zhuǎn)速工況下滾動軸承的故障診斷，需要對故障特征頻率曲線進行準確提取，通過對此提取的脊線與理論計算的軸承各類型故障特征曲線，實現(xiàn)故障類型判斷。文中提出的基于DSST的診斷算法流程如圖1所示。

圖1 文中提出的DSST算法流程

2 數(shù)值模擬仿真

為說明DSST在分析時變信號時的優(yōu)越性，首先進行數(shù)值仿真試驗，定義模擬信號為

(14)

其中：Φ(t)表示信號相位，是一個隨時間t變化的參數(shù)；φ(t)為Φ(t)的導數(shù)，即信號的理想IF。對信號添加SNR為5 dB的噪聲n(t)，圖2(a)所示為仿真信號的時域圖，圖2(b)所示為信號頻譜，仿真信號的理想IF曲線如圖3所示。

圖2 模擬信號時域圖與頻譜

圖3 理想IF曲線

由于傳統(tǒng)傅里葉變換理論上的不足，圖2(b)所示的頻譜不能夠有效展示信號包含的時變信息，因此利用STFT和DSTFT分別對信號進行分析，結(jié)果如圖4所示。對比局部放大圖可知：STFT結(jié)果時頻脊線附近發(fā)散嚴重，時頻能量不集中，而DSTFT經(jīng)過迭代得到的時頻脊線與理論值較為接近，時頻能量集中性相比于STFT更好。

圖4 STFT與DSTFT方法提供的結(jié)果比較

從脊線提取的角度看，希望在STFT和DSTFT的基礎(chǔ)上，獲得能量更加集中的時頻脊線，同步壓縮變換方法則在STFT時頻面的基礎(chǔ)上進行了能量重排，進一步銳化了時頻脊線，更加利于脊線提取。本文作者利用FSST、SST2與提出的DSST分別對仿真信號進行分析，得到的計算結(jié)果及各結(jié)果的局部放大如圖5所示。

圖5 不同TFA方法提供的結(jié)果比較

由圖5可以看出：DSST方法獲得的時頻分析結(jié)果是最佳的，它與理想時頻曲線最接近，時頻能量發(fā)散現(xiàn)象較少，可以得到清晰集中的時頻脊線。對以上幾種時頻分析方法的處理效果進行定量評判，Renyi熵是一種有效的評價指標[13]，表1所示為幾種時頻分析方法的Renyi熵值對比。

表1 幾種時頻分析算法的Renyi熵值

Renyi熵值越低，說明時頻能量發(fā)散部分越少，脊線能量越集中。由表1可以發(fā)現(xiàn)：DSST方法計算結(jié)果的Renyi熵值最低，表明其獲得的時頻譜中能量最集中，獲得的時頻脊線與理想IF曲線最吻合，驗證了DSST的優(yōu)越性。

3 試驗數(shù)據(jù)分析

為進一步驗證文中提出的DSST在實測信號分析中的有效性，利用Ottawa大學公開數(shù)據(jù)集的故障軸承數(shù)據(jù)進行試驗驗證[14-15]。該試驗在SpectraQuest機械故障模擬器(MFS-PK5M)上進行，試驗裝置如圖6所示。轉(zhuǎn)軸由電機驅(qū)動，轉(zhuǎn)速由變頻控制器控制，試驗臺左側(cè)軸承為健康軸承，試驗軸承在最右側(cè)，用不同健康狀況的軸承替換，軸承型號均為ER16K，表2為其結(jié)構(gòu)參數(shù)。在試驗軸承的上方放置加速度傳感器(ICP加速度計，型號623C01)以采集軸承振動數(shù)據(jù)。此外，還安裝了轉(zhuǎn)速計(EPC，型號775)測量軸轉(zhuǎn)速。

表2 試驗軸承參數(shù)

文中分析選用數(shù)據(jù)集中內(nèi)圈故障軸承數(shù)據(jù)，文件名為I-A-2.mat。試驗數(shù)據(jù)采樣時長為10 s、采樣頻率為200 kHz，由于數(shù)據(jù)點數(shù)過大，為減少MATLAB的計算壓力，對數(shù)據(jù)進行10倍降采樣后再截取時間長度為2 s的數(shù)據(jù)進行分析。計算的試驗數(shù)據(jù)時域圖及頻譜如圖7所示。

圖7 試驗信號時域圖與頻譜

圖7的時域圖和頻譜難以表征信號的時變特征，因此利用STFT及DSTFT對信號進行分析，結(jié)果如圖8所示?？梢钥闯觯篋STFT提取的時頻脊線能量明顯比STFT集中，分析的效果明顯優(yōu)于STFT。

圖8 STFT與DSTFT方法提供的結(jié)果比較

進一步地，利用FSST、SST2以及DSST對信號進行分析，結(jié)果如圖9所示。可以發(fā)現(xiàn)：基于STFT的FSST及SST2的時頻面能量發(fā)散嚴重，僅能模糊看到故障特征頻率的變化趨勢，而不能準確有效地提取出故障特征頻率曲線；DSST能夠準確估計出瞬時頻率，并將時頻能量往估計的瞬時頻率曲線上壓縮，計算結(jié)果時頻面清晰，故障特征頻率更加突出。

圖9 不同TFA方法提供的結(jié)果比較

進一步地，對圖9(c)中的DSST分析結(jié)果進行脊線提取，結(jié)果如圖10所示?？芍撼晒μ崛〕鲂盘柕墓收咸卣黝l率fc及其倍頻成分2fc、3fc。

圖10 DSST計算結(jié)果時頻脊線提取

圖11所示為試驗時轉(zhuǎn)速計測得的轉(zhuǎn)速n計算得到的轉(zhuǎn)頻曲線,fr=n/60。根據(jù)轉(zhuǎn)頻fr及內(nèi)圈故障特征系數(shù)δ計算軸承理論內(nèi)圈故障特征頻率fi=δ·fr。將從DSST計算結(jié)果提取出的故障特征頻率fc與理論內(nèi)圈故障特征頻率fi曲線繪制于同一時頻圖，如圖12所示?？梢园l(fā)現(xiàn)：fc與fi高度吻合，判斷該數(shù)據(jù)軸承的故障為內(nèi)圈故障，這與試驗選擇的故障軸承類型一致。通過對數(shù)據(jù)集軸承內(nèi)圈故障數(shù)據(jù)的分析，驗證了DSST在實際振動信號分析時的有效性。

圖11 轉(zhuǎn)速計測得轉(zhuǎn)頻曲線 圖12 DSST提取故障特征頻率與理論值比較

4 結(jié)論

本文作者提出一種變轉(zhuǎn)速工況下滾動軸承的故障診斷策略。針對STFT存在的調(diào)制現(xiàn)象，引入時變解調(diào)算子消去調(diào)制項的影響，通過迭代獲得接近理想的IF估計，提出了解調(diào)短時傅里葉變換(DSTFT)的計算公式。在DSTFT的基礎(chǔ)上，提出解調(diào)同步壓縮變換(DSST)，用以改善傳統(tǒng)FSST瞬時頻率估計不夠準確導致壓縮后時頻面能量發(fā)散的問題，獲得能量更加集中的時頻表達。通過仿真，驗證了該方法的有效性，并成功實現(xiàn)了對實測試驗臺滾動軸承內(nèi)圈故障的分析診斷。