解調同步壓縮變換及其應用于變轉速滾動軸承故障診斷

趙結昂，包俊義，金嵩，方琦，商高亮，劉俊

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

0 前言

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

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

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

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

1 理論描述

1.1 基于STFT的同步壓縮變換

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

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

(1)

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

(2)

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

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

(3)

(4)

(5)

1.2 解調同步壓縮變換

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

(6)

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

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

(7)

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

(8)

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

(9)

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

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

(2)對信號進行檢索得到對應時刻t的峰值數據P(t)，再對峰值數據P(t)曲線進行多項式擬合，得到信號的IF估計曲線I(t)，然后把瞬時頻率估計曲線I(t)的一階導數賦值給下一次迭代時的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。通過上面幾個步驟就能消除調制成分的影響，較為精確地得到IF估計值。

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

(11)

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

(12)

(13)

1.3 基于DSST的變轉速滾動軸承故障診斷算法

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

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

2 數值模擬仿真

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

(14)

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

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

圖3 理想IF曲線

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

圖4 STFT與DSTFT方法提供的結果比較

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

圖5 不同TFA方法提供的結果比較

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

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

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

3 試驗數據分析

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

表2 試驗軸承參數

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

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

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

圖8 STFT與DSTFT方法提供的結果比較

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

圖9 不同TFA方法提供的結果比較

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

圖10 DSST計算結果時頻脊線提取

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

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

4 結論

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