解調(diào)同步壓縮變換及其應(yīng)用于變轉(zhuǎn)速滾動(dòng)軸承故障診斷

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

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

0 前言

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

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

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

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

1 理論描述

1.1 基于STFT的同步壓縮變換

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

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

(1)

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

(2)

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

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

(3)

(4)

(5)

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

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

(6)

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

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

(7)

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

(8)

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

(9)

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

(1)初始計(jì)算時(shí)設(shè)c(t)=0，進(jìn)行DSTFT,此時(shí)的DSTFT相當(dāng)于STFT。

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

(3)重復(fù)步驟(2)依次進(jìn)行迭代運(yùn)算。

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

(10)

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

基于DSTFT的IF估計(jì)可以表示為

(11)

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

(12)

(13)

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

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

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

2 數(shù)值模擬仿真

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

(14)

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

圖2 模擬信號(hào)時(shí)域圖與頻譜

圖3 理想IF曲線

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

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

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

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

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

表1 幾種時(shí)頻分析算法的Renyi熵值

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

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

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

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

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

圖7 試驗(yàn)信號(hào)時(shí)域圖與頻譜

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

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

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

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

進(jìn)一步地，對圖9(c)中的DSST分析結(jié)果進(jìn)行脊線提取，結(jié)果如圖10所示?？芍撼晒μ崛〕鲂盘?hào)的故障特征頻率fc及其倍頻成分2fc、3fc。

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

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

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

4 結(jié)論

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