吳紅安 呂 勇 易燦燦 袁 銳

1.武漢科技大學(xué)冶金裝備及其控制教育部重點實驗室，武漢，4300812.武漢科技大學(xué)機械傳動與制造工程湖北省重點實驗室，武漢，430081

0 引言

非平穩(wěn)信號廣泛存在于自然界和人類社會生活中，如機械振動信號、通信信號、生物信號等都屬于典型的非平穩(wěn)信號。傳統(tǒng)的時頻域分析方法難以兼顧信號在時頻域的局部化特征，不能有效分析非平穩(wěn)信號[1]。所以既能反映信號頻率隨時間的變化規(guī)律，又可以描述頻域能量信息的時頻分析方法[2-3]成為分析非平穩(wěn)信號的有效工具，諸如短時傅里葉變換(short-time Fourier transform，STFT)、連續(xù)小波變換(continuous wavelet transform，CWT)、魏格納-威利分布(Vigner-Ville distribution，WVD)等時頻分析方法已被廣泛應(yīng)用[4]。STFT與CWT受海森伯格不確定性原則的限制，無法在頻率與時間方向同時獲得最佳分辨率，因此無法產(chǎn)生能量集中的時頻表示(time-frequency representation，TFR)；WVD利用信號的自相關(guān)函數(shù)進行時頻變換，在很大程度上提高了時頻表示的頻率分辨率，但這只針對于處理單分量信號，當(dāng)處理多分量信號時，會不可避免地產(chǎn)生交叉項干擾。時頻重排(reassignment method，RM)[5-6]通過調(diào)整時頻平面的能量定位來提高時頻表示的聚集性和可讀性，但卻無法進行信號重構(gòu)。

DAUBECHIES等[7]提出了一種類似于時頻重排的同步壓縮小波變換(synchrosqueezing transform，SST)方法，在CWT的基礎(chǔ)上，對變換后的小波系數(shù)進行重排，在提高信號時頻分辨率的同時很好地保持了信號重構(gòu)的特性。類似地，基于STFT的同步壓縮(SST-based STFT，F(xiàn)SST)[8]也可以提高時頻聚集性并支持可逆重構(gòu)，但是同步壓縮本質(zhì)上是基于信號由純諧波疊加的假設(shè)，因此只適用于處理平穩(wěn)或緩變信號，對于頻率快速變化的信號，其時頻表示的能量仍會在一定范圍內(nèi)發(fā)散,導(dǎo)致在時間方向上的能量模糊問題依然嚴(yán)重。

基于SST、FSST等的時頻分析方法存在抗噪性差、僅對諧波信號有較好的時頻聚集性等問題[9-10]，目前關(guān)于時頻同步壓縮的研究主要有時頻聚集性的提高，瞬頻估計精度的提高，針對多分量、變工況、強調(diào)頻信號分析的時頻分辨率的提高等[11-13]。陳向民等[14]基于自適應(yīng)時頻濾波提取變轉(zhuǎn)速齒輪故障特征來提高瞬時故障特征估計效果；YU等[15-18]先后提出了同步提取變換、多重同步壓縮變換等方法，不斷提高時頻聚集性。盡管這些方法均在一定程度上提高了時頻分辨率，但也都延續(xù)了STFT固有的帶寬大、能量發(fā)散等缺陷[17]。OBERLIN等[19]、PHAM等[20]進一步提出了高階同步壓縮變換方法，在SST的前提下進一步將頻率軸細(xì)分，從而提取更加銳化的時頻脊線，得到了更精確的瞬頻估計模型。

針對現(xiàn)有的STFT、SST和FSST等時頻方法窗口參數(shù)固定、抗噪性能差、分析非平穩(wěn)調(diào)頻信號時無法為軸承故障診斷提供足夠準(zhǔn)確的信息等問題，本文依據(jù)信號的可分性，針對頻率快速變化的非平穩(wěn)信號，提出一種窗口伸縮優(yōu)化同步壓縮短時傅里葉變換(SST-based STFT for window extension and compression optimization，WECOFSST)算法，并推廣至二階和高階的時頻變換。

1 基于短時傅里葉的時頻同步壓縮變換

短時傅里葉變換是用一個時頻局部化的固定窗函數(shù)來截取信號序列，然后進行一系列加窗數(shù)據(jù)幀的快速傅里葉變換來處理信號的時頻分析方法，通過窗口隨時間“滑動”表示信號任意時刻的頻譜特征。為便于分析非平穩(wěn)信號s(t)，可將其等效為有限個具有時變幅值和時變或固定頻率的分量信號的線性組合，設(shè)信號模型如下：

(1)

m式中，Ak(t)為信號各分量的瞬時幅值，且Ak(t)>0；φk(t)表示瞬時頻率的相位函數(shù)，且φk(t)>0；K為分量個數(shù)且為有限值。

信號s(t)的STFT為

(2)

式中，τ為信號的時移;g(τ-t)為時移窗函數(shù)。

(3)

式中，ωx(t,η)為相位時移函數(shù)，瞬時頻率估計Fx(t,η)=Re(ωx(t,η))。

FSST就是在時頻域?qū)⒆儞Q系數(shù)X(t,η)按照(t,η)→(t,ωx(t,η))的映射關(guān)系對頻率變量進行壓縮重排，表達(dá)式如下：

(4)

式中，δ(ωx(t,η)-α)為頻率重心壓縮重排函數(shù)，δ為狄拉克函數(shù)。

圖1給出了時刻t下頻率點η到α的映射過程，即對頻率點壓縮重排，遍歷所有時刻，所得結(jié)果就是FSST的結(jié)果。

圖1 FSST中時刻t下的頻率分配過程示意圖Fig.1 Schematic diagram of the frequency allocationprocess at time t in FSST

2 窗口伸縮優(yōu)化的時頻同步壓縮算法

2.1 窗口伸縮優(yōu)化的STFT和FSST

傳統(tǒng)的STFT時頻窗寬是固定不變的，使得其獲得的時頻圖能量聚集性差、脊線模糊、時頻分辨率低，導(dǎo)致時頻分析效果較差。但是常見工程信號的瞬時頻率并非平穩(wěn)變化，尤其是當(dāng)非平穩(wěn)多分量信號的各個分量比較接近時，很難確定一個固定的窗口參數(shù)使得全局的FSST時頻效果很好。

由不確定性原理可知，對于給定信號，如果能量有限，則它對應(yīng)的時寬與帶寬乘積等于常數(shù)，滿足以下關(guān)系：

(5)

式中，Δt、Δf分別為時間和頻率分辨率，若信號為Gauss信號，則等式成立。

不確定性原理表明，信號的時寬與帶寬不能同時趨于無限小，時頻分辨率相互制約，所以時頻分析應(yīng)依據(jù)不同信號選定不同的時窗函數(shù)來分析信號的局部特性[11]。本文選用Gauss窗函數(shù)來進行實驗和分析。

考慮單一頻率信號y(t)=Aei2πmt，其中，m為常數(shù)。對y(t)進行窗口伸縮優(yōu)化的STFT如下：

(6)

式中，σ(t)為t的函數(shù)，用YT-STFT(t,η,σ(t))表示時變伸縮窗口σ參數(shù)下信號y(t)的STFT。

(7)

綜上所述，對于一般信號s(t)，根據(jù)式(7)估計信號瞬時頻率函數(shù)ωs(t,η)，即

(8)

ST-STFT(t,η,σ(t))≠0

所以，信號s(t)的窗口伸縮優(yōu)化同步壓縮短時傅里葉變換(SST-based STFT for window extension and compression optimization，WECOFSST)定義為

δ(ωs(t,η,σ(t))-ξ)dη

(9)

如果g(·)為Gauss函數(shù)，則式(8)中右側(cè)第一項為零，信號瞬時頻率估計函數(shù)為

(10)

ST-STFT(t,η,σ(t))≠0

2.2 二階和高階的窗口伸縮優(yōu)化的STFT和FSST

FSST僅適用于分析由緩變純諧波組成的信號，在分析頻率快變信號時，分辨率會大幅降低，時頻脊線變得模糊。為彌補這個缺陷，研究更精確的瞬頻估計的二階同步壓縮變換，通過線性調(diào)頻模型來進行信號瞬頻的精確估計，類似于重新定義了同步壓縮中的相變換算子。

二階FSST首先定義一個二階局部調(diào)制系數(shù)qt,x(t,f)：

(11)

即調(diào)制系數(shù)qt,x(t,f)是FSST瞬頻估計關(guān)于時間t的一階導(dǎo)數(shù),所以二階的瞬時頻率估計為

(12)

則信號s(t)的二階同步壓縮變換如下：

(13)

高階的同步壓縮變換方法(FSST-N)基于信號的幅值與相位的泰勒展開，定義新的瞬時頻率估計為

(14)

類似于FSST-2，F(xiàn)SST-N的表達(dá)式為

XFSST-N(t,ω)=

(15)

(16)

SWECOFSST-2(t,ω,σ(t))=

(17)

類似于SWECOFSST-2(t,ω,σ(t))，N階窗口伸縮優(yōu)化的FSST的表達(dá)式為

SWECOFSST-N(t,ω,σ(t))=

(18)

2.3 窗口伸縮優(yōu)化確定最優(yōu)窗口參數(shù)

對于任意給定的一個信號，如果能確定一個時變窗口參數(shù)σ(t)，它可以比較清晰地反映出信號局部信息，那么就可通過最優(yōu)窗寬的FSST獲得較好的時頻分布，以確定更精確的瞬時頻率，但是實際信號并沒有任何已知的先驗信息，所以需要對適合信號時變特性的窗口參數(shù)σ(t)進行伸縮優(yōu)化并確定。

信息熵是用來度量信息不確定性的一個函數(shù)[11,13]。本文通過使用信息熵值大小來評判局部時刻的最優(yōu)窗寬。當(dāng)信號時變平穩(wěn)，則伸長對應(yīng)窗寬，以加大信號數(shù)據(jù)量，提高頻率分辨率；當(dāng)信號時變非平穩(wěn)，則壓縮對應(yīng)窗寬，提高時間分辨率，從而有效反映信號的突變情況。所以，按照信號的局部特性進行合適截取，以最小信息熵準(zhǔn)則計算截取的各段信號的最優(yōu)窗寬，能夠有效提高信號的時頻分辨率。

根據(jù)信號時頻分布和概率的相似性，各窗內(nèi)信號的概率分布如下：

(19)

則對應(yīng)的信息熵為

(20)

依據(jù)最小信息熵準(zhǔn)則，該時刻對應(yīng)的最佳窗長可計算如下：

σbest=argmin(Hs(σ))

(21)

運用局部最小信息熵，估計各個時刻的最優(yōu)窗口參數(shù)σ的流程如圖2所示。

圖2 基于最小信息熵準(zhǔn)則進行窗口伸縮優(yōu)化的時變參數(shù)估計流程示意圖Fig.2 Process diagram of time-varying parameterestimation process based on minimum informationentropy criterion for window extension andcompression optimization

3 數(shù)值仿真

3.1 合成信號仿真

為了說明提出的方法在處理非平穩(wěn)信號上的優(yōu)勢，建立了一個具有不同的調(diào)幅-調(diào)頻規(guī)律的兩個分量信號組成的多分量信號:

s(t)=s1(t)+s2(t)+n(t)

(22)

(23)

其中，n(t)是信噪比為3 dB的Gauss白噪聲。取采樣頻率100 Hz，采樣時間4 s，則多分量信號的瞬時頻率為

(24)

由圖3可以看出，信號分量s2比s1具有更強的調(diào)頻特性。

(a)信號s

理想的時頻表示應(yīng)該高度集中在瞬時頻率軌跡上，即在每個時間點上只出現(xiàn)一個頻點在時頻平面上，以描述這種單組分模式。圖4a、圖4b中STFT時頻帶寬較大，時頻譜圖的能量發(fā)散在瞬時頻率軌跡周圍，并在瞬時頻率中心位置達(dá)到最大值，而這正是由STFT的固定窗寬導(dǎo)致的，STFT具有較差的時頻分辨率。RM和FSST方法產(chǎn)生的時頻結(jié)果如圖4c～圖4h所示，RM、FSST和二階FSST(FSST-2)三種方法均顯示出比STFT頻譜圖能量更加集中的結(jié)果。從放大視圖中可以看出，對于分量s1的時頻描述，RM、FSST和FSST-2反映了相似的集中時頻表示。但對于強調(diào)頻信號分量s2，F(xiàn)SST方法能量發(fā)散嚴(yán)重，RM和二階FSST能生成較為緊湊的時頻描述，但能量仍會散布在一定范圍內(nèi)，且RM在頻率峰值變化處能量發(fā)散現(xiàn)象較二階FSST更嚴(yán)重。

(a)STFT時頻圖(b)STFT時頻圖局部放大

該仿真信號通過本文提出的窗口伸縮優(yōu)化方法(WECOFSST)進行處理，結(jié)果如圖5a～圖5d所示。對比圖4局部放大的時頻圖可以清晰地觀察到：RM、FSST和二階FSST在描述強調(diào)頻的非平穩(wěn)信號時頻分布時均會存在不同程度的能量發(fā)散現(xiàn)象，F(xiàn)SST在信號頻率快變區(qū)間，能量聚集性較差，RM和FSST-2雖能更準(zhǔn)確地描述信號的調(diào)頻時變規(guī)律，但仍存在一定的能量發(fā)散現(xiàn)象，時頻脊線模糊，無法精確估計瞬時頻率。鑒于固定窗口分辨率單一、難以反映信號局部特性的缺陷，本文提出的基于窗口伸縮優(yōu)化的FSST(WECOFSST)和二階WECOFSST方法WECOFSST-2通過伸縮優(yōu)化時變窗口來適應(yīng)信號的局部變化，使得信號在頻率變化慢的時刻具有寬的時窗，在頻率變化快的時候具有窄的時窗，從而提高信號的時頻分辨率，得到更好的時頻分布。WECOFSST和二階WECOFSST結(jié)果顯示出比RM、FSST、二階FSST結(jié)果更加能量集中的時頻表示，如圖5所示。對比觀察圖5b、圖5d可以發(fā)現(xiàn)，WECOFSST-2在描述強調(diào)頻信號時頻分布時的能量脊線比WECOFSST更加集中、更加清晰。由圖5e、圖5f可以發(fā)現(xiàn)，二階WECOFSST提取出的瞬頻與實際瞬頻曲線高度重合，瞬時頻率估計結(jié)果精確度較高。

(a) WECOFSST時頻圖(b) WECOFSST時頻圖局部放大

圖6 不同時頻分析方法下的Rényi熵Fig.6 Rényi entropy under different time-frequencyanalysis methods

為了客觀地比較每個結(jié)果，使用Rényi熵來定量評估不同方法的性能，Rényi熵值越小表示時頻表示越集中。如圖6所示，Rényi熵曲線表明，WECOFSST和二階WECOFSST結(jié)果比RM、FSST及二階FSST結(jié)果顯示出更集中的時頻表示。當(dāng)信噪比為10時，相應(yīng)的Rényi熵分別為15.11(STFT)、11.61(FSST)、11.34(RM)、11.33(FSST-2)、9.899(WECOFSST)和9.824(WECOFSST-2)。顯然，在上述時頻分析方法中，窗口伸縮優(yōu)化的二階FSST(WECOFSST-2)針對瞬時頻率快速變化的信號時頻描述效果最好，能量集中度最高。

3.2 時變轉(zhuǎn)速工況下行星齒輪箱太陽輪故障信號仿真

在本小節(jié)中，首先對變轉(zhuǎn)速工況下行星齒輪箱故障信號進行仿真分析，驗證所提方法提取敏感共振頻帶內(nèi)故障特征頻率的有效性。

為模擬太陽輪故障，建立信號模型[21]如下：

(25)

(a)信號的時域波形

太陽輪故障仿真信號的分析結(jié)果如圖8a～圖8d所示?？梢郧逦赜^察到，STFT時頻圖能量集中效果較差，僅能大致描述出信號的瞬頻變化趨勢，時頻分辨率低，時變故障特征難以提??；RM和FSST時頻圖雖獲得了較為清晰的瞬時頻率變化趨勢，但是由于邊頻帶成分密集，在瞬頻快變區(qū)間內(nèi)，沿時間方向能量依然比較模糊，如圖8b、圖8c、圖8d中標(biāo)記的紅色區(qū)域內(nèi)存在部分混疊現(xiàn)象，影響故障特征識別，所以上述方法難以精準(zhǔn)識別密集分布的時變邊頻帶。

(a)STFT時頻分布圖及局部放大圖

圖9a、圖9b所示時頻圖為本文所提方法WECOFSST和WECOFSST-2方法分析的變轉(zhuǎn)速工況下行星齒輪箱模擬太陽輪故障的仿真信號時頻表示。與FSST相比，WECOFSST和WECOFSST-2使用了更加精確的瞬時頻率估計公式，進一步細(xì)分了頻率軸，與圖8對比可以發(fā)現(xiàn)，WECOFSST和WECOFSST-2能更有效地提取快變信號的瞬時頻率特性，顯著提高時頻分辨率，使得時頻面脊線能量更加銳化集中，各個瞬時頻率變化曲線更加清晰，充分展示了算法的優(yōu)越性。

(a)WECOFSST時頻分布

4 實驗仿真分析

4.1 實測蝙蝠信號仿真

針對非平穩(wěn)多分量信號，通過引用美國伊利諾伊大學(xué)Beckman實測的大棕色蝙蝠發(fā)出的數(shù)字化回聲定位脈沖信號[22]來驗證所提方法在蝙蝠信號中的分析效果。蝙蝠通過產(chǎn)生調(diào)頻和向下掃頻信號，并借助聲波的反彈來實現(xiàn)回聲定位識別。該信號包括400個采樣點，采樣頻率為143 kHz，采樣時間間隔為7 μs。為了驗證本文所提方法的有效性，對比STFT、FSST和FSST-2分析結(jié)果，如圖10所示。蝙蝠回波信號實際上是一個多分量信號，它包含了四個能量與持續(xù)時間均不同的非線性調(diào)頻信號分量。

(a)時域波形(b)STFT時頻分布

(a)WECOFSST時頻分布

圖11所示為由WECOFSST和WECOFSST-2得到的時頻分析結(jié)果。對比圖10可以發(fā)現(xiàn)，STFT有著最差的頻率分辨率，F(xiàn)SST和FSST-2方法得到的時頻分布雖然具有較高的能量集中度，但僅限于前三個分量，第四個分量(高頻分量)幾乎無法被描述(對比圖10、圖11紅色標(biāo)記框)。相較而言，本文所提方法根據(jù)信號時變特性伸縮優(yōu)化窗寬得到了頻率帶寬更小、能量集中度更高的時頻表示，且WECOFSST-2時頻圖能量聚集性最高。WECOFSST和WECOFSST-2在FSST和FSST-2的基礎(chǔ)上不僅有能量聚集性更高的時頻表示，同時還可以清晰地表示出蝙蝠信號的高頻分量。

4.2 實驗信號分析

本文以電機啟停運轉(zhuǎn)過程來模擬實際變轉(zhuǎn)速工況，實驗裝置如圖12所示，主要由控制臺、伺服電機、齒輪變速箱(減速比Z1∶Z2=16∶24)、磁粉制動器和外置軸承座組成，加速度傳感器安裝在外置軸承座上。外置軸承座內(nèi)安裝深溝球軸承作為振動信號測試對象，軸承型號為6202，其具體參數(shù)如表1所示。

(a)實驗臺實物圖

表1 故障軸承的幾何參數(shù)

在實驗中設(shè)置伺服電機運行工況為“加速啟動—平穩(wěn)運行—減速停機”，電機啟動后轉(zhuǎn)速呈線性增大至1500 r/min，然后開始線性降低至停止運行。設(shè)置信號采樣頻率為512 Hz，采樣時間為11 s。振動加速度傳感器測得的時域信號和頻譜如圖13所示。

圖13 實驗采集的信號的時域波形和頻譜Fig.13 Signals collected in the experiment

由于故障所在軸轉(zhuǎn)速隨時間變化，傳統(tǒng)的包絡(luò)分析無法為準(zhǔn)確的軸承故障診斷提供足夠的信息，因此分別使用STFT、RM、FSST、FSST以及WECOFSST和WECOFSST-2這六種方法對該信號的幅值包絡(luò)進行時頻分析，解調(diào)出瞬時故障特征頻率，以驗證所提出方法的有效性及相對其他同類方法的優(yōu)越性。

圖14所示為STFT、RM、FSST、FSST-2四種方法時頻表示及其提取到的瞬時故障特征頻率fIF及其倍頻曲線。從圖14中可以觀察到：在時頻表示上，使用STFT和RM僅能粗略地識別出瞬時故障特征頻率，尤其是STFT在能量集中方面效果很差，難以準(zhǔn)確識別瞬時頻率；FSST和FSST-2相較于STFT，無論是在能量集中方面還是在平面清晰度等方面均優(yōu)于STFT方法，結(jié)合峰值搜尋算法可以較為準(zhǔn)確地提取出瞬時故障特征頻率，但是在識別瞬時故障特征頻率的二倍頻、三倍頻時誤差較大，精度較低。

(a)STFT時頻圖(b)STFT時頻脊線提取

圖15為本文提出的新的窗口伸縮優(yōu)化方法WECOFSST和WECOFSST-2時頻表示及其提取到的瞬時故障頻率及其倍頻曲線。從圖15中可以清晰地發(fā)現(xiàn)，窗口伸縮優(yōu)化后的方法結(jié)合峰值搜尋算法能夠較為準(zhǔn)確地提取出瞬時故障特征頻率及其二倍頻、三倍頻曲線，且精度較高，符合理論公式估算結(jié)果。

(a)WECOFSST時頻圖(b)WECOFSST時頻脊線提取

分別計算STFT、RM、FSST、FSST-2以及WECOFSST和WECOFSST-2六種方法在輸入信噪比R=20時的Rényi熵值，如表2所示。可以發(fā)現(xiàn)WECOFSST和WECOFSST-2方法的Rényi熵值較其他方法的Rényi熵值更小，其中WECOFSST-2方法Rényi熵值最小，表明本文所提新方法的瞬時頻率估計更加精確。綜合對比圖14和圖15可以發(fā)現(xiàn)，窗口伸縮優(yōu)化方法在頻率緩變時刻伸長時窗來截取信號，在頻率快變時刻壓縮時窗來截取信號，可以提高信號的時頻分辨率，得到更好的時頻分布，提取的時頻脊線較FSST方法時頻脊線更加清晰準(zhǔn)確，該方法無論是在能量集中還是在時頻平面清晰度方面均優(yōu)于STFT、RM、FSST和FSST-2方法，充分驗證了算法的可靠性和有效性。

表2 不同時頻分析方法在信噪比R=20時Rényi熵值

5 結(jié)論

(1)本文在同步壓縮變換的基礎(chǔ)上，考慮時變非平穩(wěn)信號的局部變化特性，將提出的窗口伸縮優(yōu)化短時傅里葉變換應(yīng)用于同步壓縮理論中，得到新的窗口伸縮優(yōu)化的同步壓縮短時傅里葉變換方法，并應(yīng)用于變轉(zhuǎn)速工況瞬時故障特征頻率的識別提取，提出的方法可以優(yōu)化窗口參數(shù)，較好地反映信號的局部特性，增強信號在局部和全局的能量聚集效果。

(2)算法通過窗口伸縮優(yōu)化來高效匹配信號的時變特征，利用最小信息熵值確定最優(yōu)時變窗口參數(shù)，使信號在任意時刻都具有最優(yōu)的時頻分辨率和較高的魯棒性，在時頻表示的清晰度上與現(xiàn)有方法相比具有較大的優(yōu)勢。

(3)該方法應(yīng)用在旋轉(zhuǎn)機械變轉(zhuǎn)速工況瞬時頻率估計上，適用于處理非平穩(wěn)的多分量信號，相較于其他方法，能獲得較高精度的估計結(jié)果，且具有較高的可靠性。