基于QSTFT的同步壓縮變換及在軸承故障診斷中的應用

葉杰凱，湯小明，林青云，梁登，吳華盛

(麗水市特種設備檢測院，浙江 麗水 323000)

0 引言

信號處理是一種利用數學算子從信號中提取特征信息的反向處理方法。傳統的傅里葉方法只能處理平穩(wěn)信號，并且是基于全局的頻譜分析，對于頻率隨時間變化的非平穩(wěn)信號往往無能為力。而時頻分析方法能夠有效地揭示時變性特征，已被廣泛應用于從非平穩(wěn)信號中提取特征信息[1-2]。

傳統的時頻分析方法有連續(xù)小波變換(CWT)和短時傅里葉變換(STFT)，常被用來描述信號在時頻平面上的特性[3]。然而，它們都有同樣的局限性，就是所謂的“不確定性原理”，即不能在時間和頻率上任意精確地定位一個信號，也就是常說的時頻模糊問題。這種局限性導致后續(xù)的瞬時頻率提取和信號重構等方面會存在不足，難以較為清晰地提取出故障信號及其分量。因此，針對時頻平面去模糊化的研究一直都是一個熱門方向，時頻重排和同步壓縮變換(synchrosqueezing transform, SST)就是其中一個較好的思路。

自從DAUBECGES I等[4]提出了同步壓縮變換來揭示非平穩(wěn)信號復雜的時頻特性，在醫(yī)學、地球物理學和音頻等多個領域對其都有所應用。SST是一種特殊的重分配算法，它對應于一種提高時頻分辨率的非線性算子。在SST中，時頻系數僅在頻率軸上重新分配，使其既能夠簡化分配過程又能與原始參量保持聯系[5-6]。然而，SST只能使瞬時頻率“慢變”信號的時頻表示銳化，并且當處理強時變信號如非線性調頻信號時，SST不能產生集中的時頻表達。這是因為SST中的頻率重分配算子不能為強時變情況提供準確的無偏估計。為此，PHAM D H等[7]提出了高階同步壓縮變換方法(high-order synchrosgueezing transform, HSST)。在高階振幅和相位近似的基礎上定義了新的同步壓縮算子，并對快速變化的瞬時頻率(IF)信號產生高度集中的時頻表達。但是，HSST是以傳統短時傅里葉變換STFT為基礎，在時頻變換窗寬上無法做到自適應選擇。同時，BERRIAN A等[8]提出了基于縫式短時傅里葉變換(quilted STFT, QSTFT)的同步壓縮變換方法，即SST-QSTFT。該方法對復雜多組分信號的瞬時頻率無法做到準確估計，再加上對噪聲極為敏感，因此在強背景噪聲下，很難獲得清晰的時頻表示[9-10]。但是，QSTFT重新定義了一組隨時頻變化的自適應窗口函數集合，具有明顯的自適應特性，其可以使用最合適的窗函數來適應不同的時變信號，自動匹配信號的局部變化，以此來提高STFT的時頻分辨率。

綜上，本文引入基于縫式短時傅里葉變換QSTFT的自適應窗口計算方法，并將該方法應用到HSST中，提出一種基于縫式短時傅里葉變換的自適應高階同步壓縮變換(adaptive high-order synchrosqueezing transform based on a quilted short-time fourier transform, AHSST-QSTFT)。

1 理論描述

1.1 基于QSTFT的同步壓縮變換

定義一個調幅調頻(AM-FM)信號，令其表達式為

s(t)=A(t)eiφ(t)

(1)

式中A(t)和φ(t)分別是其瞬時幅值和瞬時相位。

繼而信號的STFT變換可以寫成下面的形式：

(2)

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

那么，SST-STFT可以用下面的公式表示：

(3)

(4)

將時頻面任意一個局部區(qū)域定義為Ω?R2,且該區(qū)域對應一個縫窗口函數集合hΩ。令ht,ω=hΩ，對其中每一個(t,ω)∈Ω。則函數集合h的表達式為

h(τ,t,ω)=ht,ω(τ)

(5)

式中t,τ∈R且ω∈R+。

對信號s(t)，本文研究的QSTFT的表達式為

(6)

因此，由式(3)和式(6)可得SST-QSTFT表達式為

(7)

1.2 高階同步壓縮變換

SST采用零階估計來對信號的頻率進行估計，僅在處理慢時變的中頻信號時效果較好。而高階SST方法則是基于信號振幅和相位的高階Taylor展開，對信號的IF進行更加精確的估計，以此來提高在特征頻率上面的能量，達到去模糊化的目的。

AM-FM信號的Taylor展開式如下：

(8)

式中φ(k)(t)表示φ(t)的第k階導數。

那么，式(2)可以改寫成如下形式：

(9)

(10)

(11)

此時，HSST的表達式可以寫成

(12)

1.3 基于QSTFT的自適應高階同步壓縮變換

(13)

最后，根據式(7)、式(11)和式(13)可知，AHSST-QSTFT的表達式為

(14)

因此，本文提出的基于QSTFT的自適應高階同步壓縮變換算法AHSST-QSTFT的具體實現步驟如下：

1)定義縫窗口函數，對STFT的窗函數進行改進，得到自適應窗函數集合；

3)按照式(14)的方法重置高階瞬時頻率估計，從而獲得改進的自適應高階時頻表達QHTs(t,ω)，并畫出時頻圖。

2 仿真分析

為了驗證該方法的可行性，構造兩個具有AM-FM性質的分量信號：

(15)

其中：信號的采樣頻率設為4 096 Hz；信號長度為4 096；信號的理想瞬時頻率曲線如圖1所示(本刊為黑白印刷，如有疑問請咨詢作者)。為了更加貼近實際情況，使模擬信號更有代表性，將兩個信號進行了疊加，并且添加了SNR=-2的高斯白噪聲(信噪比公式SNR=10lg(Ps/Pn)，得到信號x(t)的時域波形如圖2所示，頻域波形如圖3所示。由FFT計算的結果可以看出，在噪聲干擾較大的情況下，無法從頻譜圖中區(qū)分和辨識時變的信號特征。

圖1 信號x1(t)和x2(t)的頻率曲線

圖2 信號x1(t)、x2(t)和x(t)的時域波形圖

圖3 信號x(t)的頻域波形圖

因此，為了說明本文所提出算法的優(yōu)勢，分別畫出了仿真信號的4階HSST和AHSST-QSTFT的時頻圖，并且采用相同的脊曲線提取算法對分量進行了提取，結果見圖4。其中虛線為真實的瞬時頻率，實線為提取的脊線。

圖4 時頻表示及其脊線提取結果

通過兩組時頻圖的對比，可以得出AHSST-QSTFT方法能夠獲得更好的時頻表達效果，并且提取到的瞬時頻率更加接近于真實情況。

3 實驗分析

故障診斷對于機械設備的預知維修具有重要意義。根據RANDALL R B等提出的滾動軸承模型[11]，為了進一步驗證基于AHSST-QSTFT方法在機械設備故障診斷中的有效性，本文設計了如下實驗。

圖5為試驗裝置的故障模擬器。設備由550 W的交流電機驅動。故障軸承安裝在最右側的皮帶輪上，利用電火花加工在目標軸承加工4個點蝕坑，坑的深度為0.8 mm，并盡可能靠近支撐軸承的機殼上方垂直固定振動加速度傳感器，用以進行數據采集，傳感器位置如圖中箭頭所示。實驗采用B&K3560C數據采集分析儀進行數據采集，采樣頻率為16 384 Hz，轉速1 450 r/min；滾動體個數Z=9，試驗中的滾動軸承型號為6207，壓力角為α=0°，故障類型為軸承外圈故障，可以根據下面的計算得到故障頻率：

圖5 軸承故障模擬實驗臺

(16)

式中：fr為轉頻，Hz；fo為軸承外圈故障特征頻率，Hz。

選擇8 192個采樣數據分析，實際測得的軸承振動信號如圖6所示，在時域波形上能夠明顯看出瞬態(tài)沖擊故障。圖7是該信號的頻譜分析，可以觀察到峰值區(qū)域主要集中在1 000 Hz、4 000 Hz附近。這是因為軸承發(fā)生故障時，當轉子經過缺陷位置時，相當于是在軸承上加上了一個沖擊力所激發(fā)出來的信號，反映的是該系統的固有頻率。

圖6 實測信號的時域波形

圖7 實測信號的頻譜圖

考慮到實際信號的能量都集中在共振區(qū)間，不利于最終故障的識別。因此在本文中首先對實際信號進行希爾伯特變換，得到對應的包絡信號，然后再采用AHSST-QSTFT方法對信號進行時頻分析。

圖8(a)和圖8(b)分別是在相同的參數設置下采用HSST和AHSST-QSTFT方法得到的時頻表示。考慮到軸承外圈的故障頻率為87.01 Hz，一般地，如果能夠找到故障特征頻率對應的多個倍頻就能夠確定故障。圖8(a)雖然能夠較為模糊地看到某些對應的特征頻率，但受到背景噪聲的影響，它們在時頻平面上表示得并不是很突出。圖8(b)為采用提出的AHSST-QSTFT方法得到的時頻表示，相比于圖8(a)，在圖中可以清晰地看到故障特征頻率fo(87.01 Hz附近)及其2倍頻(174 Hz附近)、3倍頻(262 Hz附近)、4倍頻(349 Hz附近)、5倍頻(436 Hz附近)、6倍頻(523 Hz附近)、7倍頻(610 Hz附近)瞬時頻率曲線，并且它們基本上凸顯在時頻平面上，因此，可以準確地判斷故障類型為軸承外圈故障。

圖8 不同方法得到的時頻表示

實驗結果表明，AHSST-QSTFT方法相較HSST具有更好的抗噪聲干擾能力，其最大程度地凸顯了故障特征信息，增強了時頻平面故障特征的辨識度。

4 結語

本文提出了一種基于QSTFT的自適應高階同步壓縮變換算法，充分利用了QSTFT的窗口自適應性和HSST的時頻壓縮特性。該方法提高了經典STFT框架下振動信號的窗口選擇適應性，使重排后頻率估計最大程度地接近信號的真實瞬時頻率，提高了時頻平面的清晰度，進一步增強了HSST方法在處理時變信號時的分辨能力。同時，數值仿真和軸承外圈故障實驗表明，本文提出的方法能夠有效降低噪聲的干擾，并且在提取更為準確的故障特征信息方面更具優(yōu)勢。