廣義瞬時速度同步化分步解調(diào)變換及其對旋轉(zhuǎn)機械振動信號分析

石娟娟，花澤暉，沈長青,，江星星，馮毅雄，朱忠奎，孔 林

(1.蘇州大學 軌道交通學院，江蘇 蘇州 215131；2.浙江大學 流體動力與機電系統(tǒng)國家重點實驗室，杭州 310027；3.長光衛(wèi)星技術(shù)有限公司，長春 130102)

旋轉(zhuǎn)機械如汽輪機、風機、發(fā)動機被廣泛用于電力和航空航天等領域。齒輪和軸承作為最常見的旋轉(zhuǎn)機械零部件，在機械系統(tǒng)中發(fā)揮著連接和支撐機械結(jié)構(gòu)的重要作用。隨著實際需求的發(fā)展，旋轉(zhuǎn)機械系統(tǒng)常被要求服役于高速重載等惡劣的復雜工況下，而這些關(guān)鍵零部件一旦發(fā)生意外故障，會直接影響到整個系統(tǒng)的平穩(wěn)運行，嚴重時甚至還會造成經(jīng)濟損失和安全事故。同時，機械系統(tǒng)設計結(jié)構(gòu)和工藝的復雜性以及可能存在的安裝誤差也會引起機械系統(tǒng)的不良振動并導致故障。因此，實時監(jiān)測機械系統(tǒng)中關(guān)鍵零部件的健康狀態(tài)可更加快速和精確地定位故障并做出相應的診斷和維護措施。通常旋轉(zhuǎn)機械的故障診斷方法可通過振動、噪聲、溫度、聲發(fā)射和油液等方法來進行，其中，振動檢測涉及信號處理方面，而旋轉(zhuǎn)機械振動信號常包含揭示系統(tǒng)健康狀況的有用信息，是最常用的診斷方法。當旋轉(zhuǎn)機械發(fā)生故障時，常表現(xiàn)為振動頻率的變化，可通過檢測振動頻率、加速度、相位等參數(shù)進行分析。當旋轉(zhuǎn)機械運行于變速工況下，信號的故障特征頻率也是隨時間而變化的，因此需要采用時頻分析方法來揭示信號中頻率隨時間的變化。加上信號采集過程中不可避免地存在噪聲干擾，還有加速度計有時無法安裝來獲取瞬時頻率的信息等問題都給變轉(zhuǎn)速工況下振動信號處理帶來了新的問題和挑戰(zhàn)。因此，需要提出更加先進和有效的時頻分析方法來提高旋轉(zhuǎn)機械振動信號中特征提取的準確性以及時獲知系統(tǒng)關(guān)鍵零部件的健康狀態(tài)，從而保障整個機械系統(tǒng)的安全運行和維護。

時頻分析方法因能同時揭示信號的時間和頻率信息，廣泛用于非平穩(wěn)信號的處理[1]。在時頻分析方法中，最常見的是線性分析方法，如短時傅里葉變換(short-time Fourier transform, STFT)和小波變換(wavelet transform, WT)等[2]。之所以稱為線性變換方法，是因為這類方法都利用了加窗信號和基函數(shù)的內(nèi)積。這類方法受有限時頻分辨率影響存在嚴重的時頻模糊問題。類似的，還有一類雙線性變換方法，如維格維納分布、Cohen經(jīng)典分布等。這類方法雖能提供較線性變換更高的時頻分辨率，但存在嚴重的交叉項干擾問題。對線性變換得到的時頻表示進行時頻重映射這一類方法稱為時頻后處理，可進一步解決時頻分析中因固定的時頻分辨率導致的時頻模糊問題。如Auger等[3]提出沿著時間和頻率方向重分配時頻分布的RM(reassignment method)方法來提升時頻表示的可讀性；Daubechies等[4]提出了同步壓縮變換(synchrosqueezing transform, SST)結(jié)合小波變換實現(xiàn)信號分解。與RM在時間和頻率同時壓縮不同，SST僅沿頻率方向重新分配時頻表示，保留了信號的重構(gòu)能力。雖然SST將時頻分布集中在目標時頻脊線周圍，但是當信號的頻率變化較快時，由時頻分辨率本身存在的模糊問題仍未解決。并且當待分析信號包含較大的噪聲時，噪聲引起的時頻模糊也被擠壓到了時頻脊線上，導致時頻表示中目標頻率分量處的能量聚集性降低，從而使時頻脊線不清晰。類似的，Yu等[5]提出了同步提取變換(synchroextracting transform, SET)，通過將得到的時頻表示與一個脈沖函數(shù)相乘，達到只保留目標頻率處時頻分布的效果從而提升時頻表示可讀性。這些時頻重分配方法都是在不改變時頻分辨率的情況下通過后處理從而進一步提升時頻表示。除此之外，還有一類信號分解方法，如經(jīng)驗模態(tài)分解(empirical mode decomposition, EMD)、局部均值分解(local mean decomposition, LMD)和改進的變分模態(tài)分解(variational mode decomposition, VMD)等，這些方法旨在從原始的多分量信號中分離出目標頻率分量(如旋轉(zhuǎn)機械振動信號中的與故障特征相關(guān)的頻率分量)，然后使用時頻分析方法再對單個分量進行分析。這一類方法更適用于定轉(zhuǎn)速工況下的軸承故障診斷，因為定轉(zhuǎn)速工況下故障特征頻率是一個定值，不隨時間而變化，因此可根據(jù)頻帶劃分而分離出來特征頻率分量。因此這類方法還常與濾波和譜峭度結(jié)合，譜峭度指導選取含有最多相關(guān)特征所在的頻帶范圍，然后經(jīng)過帶通濾波后再次分析可得到更加清晰的結(jié)果[6-8]。

為了解決時頻分析中時間和頻率上存在的時頻模糊問題，廣義解調(diào)(generalized demodulation, GD)因其能將變化的時頻曲線解調(diào)到固定頻率從而提升能量聚集性這一特性可有效提升時頻表示的可讀性。為了解決SST在分析頻率快速變化信號時存在的頻率模糊問題，提出了GD和SST結(jié)合的新方法[9-10]。該組合的優(yōu)勢在于廣義解調(diào)可增強時頻表示的能量聚集性，SST用于進一步銳化時頻脊線，從而可進一步提升時頻圖的可讀性。Wang等[11]提出迭代使用匹配解調(diào)變換方法(matching demodulation transform, MDT)來改善提取的瞬時頻率的準確性；Shi等[12-13]提出了基于分形維數(shù)和廣義解調(diào)的變速軸承故障診斷方法，避免了處理變速信號時的重采樣過程。上述方法都已成功用于仿真分析和試驗驗證。但是，前述方法依賴瞬時頻率相關(guān)先驗知識，需要有迭代過程來改善瞬時頻率估計的準確性和處理多分量信號從而提升時頻表示[14]。

為了不依賴瞬時頻率等先驗知識并成功地運用廣義解調(diào)以增強時頻表示的能量聚集性，本文提出了廣義瞬時頻率同步化分步解調(diào)變換。該變換不需對瞬時頻率進行預估計，并且避免了在處理多分量信號時所需的迭代操作。在足夠小的時間窗內(nèi)信號的瞬時頻率脊線可近似為直線，該直線的斜率以及與橫坐標軸的夾角也是唯一確定的。然后根據(jù)使用真實解調(diào)因子進行解調(diào)變換可將信號解調(diào)至水平直線并得到最高的能量聚集性這一特點，提出了最大峭度指導的自適應角度選取方法。為有效處理多分量信號，對原變換基函數(shù)進行擴展，得到與瞬時頻率同步的線性變換基函數(shù)，由此可同步增強多個頻率分量的時頻表示，適合處理旋轉(zhuǎn)機械所產(chǎn)生的多分量振動信號并且避免了算法中的迭代過程。此外，考慮到窗長的變化也會對時頻圖可讀性產(chǎn)生影響，因此同時變化角度和窗長，利用峭度最大準則，提出在每個時刻選取合適參數(shù)組合的方法，使所提方法可在一定的噪聲下自適應選擇合適的參數(shù)并避免對先驗知識的依賴。根據(jù)提升的時頻表示，可更加準確地從非平穩(wěn)振動信號中提取出相關(guān)時頻特征并用于變轉(zhuǎn)速旋轉(zhuǎn)機械的健康狀態(tài)監(jiān)測。

1 廣義瞬時速度同步化分步解調(diào)變換

1.1 廣義分步解調(diào)變換

廣義解調(diào)可處理非平穩(wěn)信號并增強時頻聚集性，在于該方法可將時變的頻率曲線解調(diào)到固定頻率上。通過這一步解調(diào)操作，還可有效地把多分量信號分解成多個代表不同分量的單分量信號。給定信號x(t)，廣義解調(diào)變換可表示成

(1)

式中：d(t)為解調(diào)因子；e-j2πd(t)為變換基函數(shù)?？梢园l(fā)現(xiàn),x(t)的廣義解調(diào)變換等同于對x(t)e-j2πd(t)的標準傅里葉變換。當d(t)=0，式(1)就等價于標準的傅里葉變換。類似地，根據(jù)傅里葉變換的性質(zhì)，可定義逆廣義傅里葉變換為

(2)

觀察式(2)，可以發(fā)現(xiàn)如果XG(f)≡δ(f-f0)，x(t)可以被寫成

x(t)=ej2π[f0t+d(t)]=ejφ(t)

(3)

式中,φ(t)為信號的瞬時相位。由此可發(fā)現(xiàn)，信號原先的瞬時頻率f(t)=f0+d′(t)=dφ(t)/2πdt被投影到一固定頻率f=f0，其中d′(t)表示d(t)的一階導。根據(jù)式(3)可知，如需將時變頻率曲線映射到一固定頻率，只需計算相關(guān)的解調(diào)因子并對信號進行解調(diào)。根據(jù)上述分析，可總結(jié)出廣義解調(diào)的兩個主要特性：①時變的瞬時頻率曲線可被解調(diào)到某一恒定頻率；②解調(diào)后信號能量集中在頻率f=f0，較解調(diào)之前的信號具有更高的能量集中性。

據(jù)此，Shi等提出廣義分步解調(diào)來研究振動信號時頻特性，針對每一段截取信號將所提方法分為兩步：①前向解調(diào)增強時頻聚集性；②后向解調(diào)將時頻能量聚集性得到增強的瞬時頻率脊線返回至真實頻率。具體介紹如下：短時窗內(nèi)的信號頻率可寫成

f(t)=f0+d′(t)

(4)

前向解調(diào)因子可計算為

xd f(t)=x(t)·df(t)=ej2π[f0t+d(t)]e-j2πd(t)=ej2πf0t

(5)

式中：df(t)為前向解調(diào)因子；xd f為前向解調(diào)后的信號。式(5)的物理含義即為信號在乘上前向解調(diào)因子后，原先變化的頻率f(t)被變換到了f0處。

為更清晰說明上述過程，廣義分步解調(diào)變換的示意圖，如圖1所示。從圖1中可以看出，在時刻t=τ原先的頻率曲線f0+d′(t)被前向解調(diào)到了頻率f0處，比較圖1(b)中的幅值和頻率的關(guān)系可發(fā)現(xiàn)解調(diào)后的時頻脊線能量更加集中，即時頻分布的能量相比解調(diào)之前更加聚集在頻率f0周圍，但在時間段[τ-Δt,τ+Δt]內(nèi)，信號前向解調(diào)后的頻率并非信號原先的真實頻率，因此還需要后續(xù)操作(如后向解調(diào))來恢復分析信號的真實頻率，如圖1(c)所示。在后向解調(diào)中，時頻聚集性與前向解調(diào)之后一致，僅發(fā)生頻移，由此后向解調(diào)可定義為

FD為前向解調(diào)；BD為后向解調(diào)。

xd(t)=xd f(t)db(t)=ej2πf0tej2πd′(τ)t=ej2π[f0t+d′(τ)t]

(6)

式中：db(t)為后向解調(diào)因子；d′(τ)>0為脊線上移；d′(τ)<0為下移。通過上述分步解調(diào)過程，可得到能量分布更加集中的時頻表示。但是，上述的解調(diào)過程依賴瞬時頻率的先驗知識，即需首先對瞬時頻率進行預估計，據(jù)此構(gòu)建正確的解調(diào)因子來對信號進行解調(diào)變換。在不知道瞬時頻率的前提下，目前已有文獻大多是從得到的時頻表示中先粗略地估計瞬時頻率來對信號進行解調(diào)，然后利用迭代操作來提升瞬時頻率估計的準確性。一方面，解調(diào)因子的準確性極大地依賴瞬時頻率的預估計；另一方面，涉及的迭代操作使算法本身更復雜，并且在處理多分量諧波信號時，還需要再次涉及迭代操作來提升多分量的時頻表示。

1.2 廣義瞬時速度同步化分步解調(diào)變換

為避免對瞬時頻率預估計的依賴并簡化算法對多分量信號所需的迭代過程，提出了本文的研究方法。

與廣義分步解調(diào)變換相同，在短時窗截取的信號，該段時間窗內(nèi)的瞬時頻率可近似看成線性變化，根據(jù)泰勒展開多項式，該瞬時頻率可寫成

f(t)≈f′(τ)(t-τ)+υ

(7)

式中：f(t)為信號的瞬時頻率；f′(τ)為瞬時頻率在t=τ處的斜率；υ為該時刻t=τ的頻率。該瞬時頻率還可寫成

(8)

引入?yún)?shù)α，定義為

(9)

將式(9)代入式(8)，瞬時頻率可表示成

f(t)=υ[tanα(τ)(t-τ)+1]

(10)

式(10)表明，引入隨時間變化的參數(shù)α來對斜率進行描述，由此可用角度有限的變化范圍(-π/2, π/2)來映射斜率的變化范圍。根據(jù)式(4)，解調(diào)因子滿足

(11)

將式(10)代入式(11)，即可計算得到前向解調(diào)因子的表達式，寫成

(12)

式中：f0為解調(diào)后的頻率；υ為時刻t=τ的頻率?？煽闯觯?12)中前向解調(diào)因子是與tanα有關(guān)的函數(shù)。只要確定tanα(τ)，便可對信號進行前向解調(diào)，前向解調(diào)后的信號可表示成

(13)

式(13)表明，只要角度α(τ)選取正確，能正確匹配待分析信號的真實頻率變化，該加窗信號即可被解調(diào)到固定頻率f=f0處。

然而，前向解調(diào)后信號所處頻率并不是待分析信號在該時刻的真實頻率，因此還需對前向解調(diào)后的信號進行后向解調(diào)以恢復其真實頻率。根據(jù)傅里葉變換性質(zhì)，頻域內(nèi)發(fā)生頻移，等同于在時域內(nèi)乘上復指數(shù)函數(shù)，可表示為

ej2πυ(t-τ)

(14)

式中,db(t)為后向解調(diào)因子。通過觀察式(14)可發(fā)現(xiàn)，兩次解調(diào)后的真實頻率得以保持，而解調(diào)后信號頻譜在當前時刻的能量聚集性得到了提高，即原先線性的頻率直線先被解調(diào)至固定頻率以聚集能量，然后再通過第二次解調(diào)恢復信號的真實頻率。整合前、后向解調(diào)因子，可得完整的解調(diào)因子為

d(t)=df(t)·db(t)=

(15)

根據(jù)式(15)定義的解調(diào)因子，可完成對信號的前、后向解調(diào)。原始信號的加窗時頻表示為

XSTFT(τ,υ)=〈x(t)g(t-τ),ej2πυ(t-τ)〉=

(16)

式中：g(t-τ)為窗函數(shù)；υ為時刻t=τ的頻率?？紤]式(15)定義的解調(diào)因子，則解調(diào)后信號時頻表示可寫成

XG[τ,υ,α(τ)]=

(17)

式(17)即為引入了傾斜角α的廣義分布解調(diào)變換。

所提變換還可滿足對信號中目標頻率分量的時域重構(gòu)，具體證明為

2πg(shù)(0)x(τ)

(18)

式中，δ為脈沖函數(shù)。

根據(jù)式(15)～式(17)可發(fā)現(xiàn)，所提變換關(guān)鍵在于不依賴先驗知識而找到正確的角度參數(shù)α以成功解調(diào)待分析信號并提升其時頻聚集性。同時，為對解調(diào)之后信號的時頻聚集性進行量化描述，引入峭度指標。通過離散化角度α，分別對所截取的信號進行分析，并計算出每個角度所對應的峭度值，其中，正確的角度因能匹配真實頻率(與真實頻率最接近)，可將該分段信號解調(diào)至一固定頻率并具有最高的時頻聚集性，即對應峭度的最大值。因此，可根據(jù)峭度最大值來選取準確的角度α對原始信號進行解調(diào)，從而獲得更佳的時頻聚集性。相反，若所選取角度α不能反映該截取信號真實的瞬時頻率，則解調(diào)之后信號的頻率不為定值(解調(diào)后的時頻脊線仍無法平行于時間軸)。根據(jù)廣義解調(diào)的性質(zhì)，此時該時頻表示中信號能量相對分散，對應的峭度值也更小。

為對上述分析進行說明，定義仿真信號x(t)，寫成

(19)

式中,f(t)為瞬時頻率，寫成

f(t)=40sin(0.36t+0.6)+10sin(2t-2.5)+

5.4sin(3.6t-7.6)

(20)

所定義信號原先的時頻表示及廣義解調(diào)后信號的時頻表示如圖2(a)所示，時變的瞬時頻率首先被變換到了固定頻率f0。圖2(b)中，解調(diào)因子由真實角度計算得到，而在圖2(c)中，解調(diào)因子是由有別于真實值的角度計算所得，對比兩個時頻表示的譜圖，如圖2(d)所示，用真實解調(diào)因子解調(diào)所得信號的時頻能量更集中(實線代表采用真實角度計算得到的解調(diào)因子對信號進行解調(diào)所得時頻幅值；虛線為與真實值不同的角度得到的結(jié)果)，同時將兩種解調(diào)結(jié)果與真實值對比(原先真實的時頻能量如圖2(d)中點劃線所示)。頻譜分析結(jié)果也說明用真實的解調(diào)因子可將頻率曲線解調(diào)至水平線并獲得較高的能量，即解調(diào)后能量的聚集性可用于指導角度的選取，從而實現(xiàn)在無需預先知曉瞬時頻率(瞬時速度)的前提下對信號進行廣義解調(diào)變換以增強時頻聚集性。

圖2 采用不同角度對信號進行解調(diào)分析的結(jié)果

對于含N個點的頻譜S，峭度K可被計算成

(21)

為了正確選取式(9)定義的角度，對角度進行離散化，寫為

(22)

式中，I決定了角度的個數(shù)。將式(22)代入式(12)，前向解調(diào)因子可寫成

(23)

式(23)表明，每個角度α都可計算得到一個特定的解調(diào)因子。角度離散化的個數(shù)增多會提高解調(diào)因子估計的準確性，但同時也會導致計算量的增大。角度和最大峭度的對應關(guān)系可表示為

(24)

式中，α*(τ)為時刻τ處由頻譜計算得到的最大峭度值所對應的角度。

為驗證最大峭度指導的角度是否正確，對式(20)提出的仿真信號進行詳細分析，得到的結(jié)果如圖3所示。該單分量信號的短時傅里葉分析結(jié)果如圖3(a)所示，矩形代表分步解調(diào)的時間區(qū)間，虛線代表t=1.155 s時刻的時頻表示切片。根據(jù)式(21)和式(23)，圖3(a)中虛線處的時頻分布在經(jīng)不同角度構(gòu)建的解調(diào)因子作用后所得時頻表示的峭度變化規(guī)律如圖3(b)所示，在圖3(b)中三角標記出了一系列角度對應的峭度的最大值。通過滑動時間窗，同步估計每個時刻點最大峭度值對應的角度，結(jié)果如圖3(c)所示。另外，圖3(c)還給出了由式(20)直接計算得到的真實角度來驗證峭度指導角度選取的準確性，經(jīng)對比可以發(fā)現(xiàn)，由峭度指導選取的角度變化規(guī)律與真實情況基本吻合，說明峭度可用于判斷所選的角度是否合適。在數(shù)值上，分析了所估計的角度與真實值的平均相對誤差，經(jīng)計算誤差為0.98%，說明了所提峭度指導的角度估計策略的有效性。最終，所提方法得到的時頻表示如圖3(d)所示。通過與STFT結(jié)果(見圖3(a))對比可以發(fā)現(xiàn)所提方法可得到時頻能量更加集中的時頻脊線和更高的時頻聚集性，尤其是在待分析信號的頻率變化較快時，如2～4 s。

圖3 單分量信號的角度估計

除了常見的峭度指標，平滑因子(smoothness index, SI)也常被用于衡量信號的沖擊程度。為進行比較，將平滑因子也用于指導角度的選取。選用式(21)中的頻譜S，平滑因子可定義為

(25)

式(25)表明，平滑因子SSI可看成是幾何均值與算術(shù)均值的比值，范圍在0～1。當頻譜中存在明顯的脈沖聚集時，SSI更接近于0?？紤]到時頻表示中信號更聚集的時頻能量在頻譜中表現(xiàn)為幅值更高的脈沖成分，如圖2(b)～圖2(d)所示。因此，所估計角度和平滑因子的對應關(guān)系為

(26)

式中，α*(τ)為時刻τ處由頻譜計算得到的最小平滑因子所對應的角度。

再次對圖3中所用信號進行分析，峭度和平滑因子指導的角度結(jié)果，如圖4所示。計算平滑因子和峭度所估計角度的平均相對誤差，分別為4.86%和0.99%。經(jīng)對比可發(fā)現(xiàn)平滑因子雖能揭示角度大致變化趨勢，但在信號首尾兩端估計精度沒有峭度高且同時存在振蕩現(xiàn)象。因此可得出結(jié)論：峭度指標可更準確地指導角度的選取。

圖4 平滑因子和峭度所估計角度

該信號的重構(gòu)結(jié)果如圖5所示，經(jīng)計算，重構(gòu)信號與原始信號的平均相對誤差為1.07%，說明了重構(gòu)的準確性。

圖5 單分量信號時域重構(gòu)結(jié)果

所提方法的簡要流程如圖6所示。在每個時刻依據(jù)最大峭度準則估計得到合適的角度，通過構(gòu)建合適的解調(diào)因子匹配待分析信號的頻率變化趨勢，從而提升時頻表示。

圖6 最大峭度指導下的角度估計流程圖

1.3 改進的線性變換基函數(shù)

上述分析是以假設信號中僅有單個頻率分量為前提進行的，但是在實際應用中，信號常常包含多個分量，如：旋轉(zhuǎn)機械的振動信號中就包含多個分量，有軸轉(zhuǎn)速及其諧波和故障特征頻率及其諧波等多個分量。因此在分析多分量信號時，傳統(tǒng)的分析方法(如可將STFT看成用角度α=0來解調(diào)信號)還可能會產(chǎn)生額外的時頻模糊問題。為了能同步增強多個頻率分量的時頻表示，研究擴展后的線性變換基函數(shù)。先定義一個含有多個頻率分量的信號，寫成

(27)

式中,fk(t)=Ckf(t)，Ck=2.0,Ck=3.2,Ck=3.5,Ck=3.8,Ck=5.0，f(t)為瞬時頻率并在式(20)定義。式(17)中定義的線性變換基函數(shù)，還滿足如下關(guān)系

(28)

對該多分量信號進行時頻分析，STFT結(jié)果如圖7(a)所示，圖7(a)中矩形區(qū)域經(jīng)放大后的結(jié)果如圖7(b)所示。結(jié)合這兩個圖可以看出，STFT的時頻表示受較嚴重的交叉項干擾，圖中的三個相對較近的頻率分量難以被區(qū)分開，并且得到的時頻聚集性欠缺。采用所提方法對該信號進行分析，得到的分析結(jié)果如圖7(c)所示。圖7(c)中矩形區(qū)域?qū)臅r頻局部放大結(jié)果如圖7(d)所示。與圖7(a)和圖7(b)中STFT的結(jié)果相比，可以發(fā)現(xiàn)，所提方法可得到能量更加集中的時頻表示和更容易分辨的時頻脊線，并且在一定程度上解決了STFT中因基函數(shù)頻率與待分析信號時頻特征不匹配而引起的交叉項干擾。

廣義瞬時頻率同步化分步解調(diào)變換在同步增強多分量時頻表示時，峭度指導下的角度選取方法與單分量一致，估計的角度如圖7(e)和圖7(f)所示。在時刻t=1.155 s處峭度隨角度的變化如圖7(e)所示，圖7(e)中，三角表示最大峭度，表明該最大峭度值對應的角度為該時刻最終選取的角度。同步估計所有時刻的角度，結(jié)果如圖7(f)所示，其中，實線代表估計所得角度，虛線代表真實角度，可見所估計的角度與真實值基本吻合。經(jīng)計算，角度估計的平均相對誤差為0.94%，說明在最大峭度指導下所提方法能同時提升多個頻率分量的時頻聚集性。對所得時頻表示進行瞬時頻率估計，然后根據(jù)估計的目標頻率分量進行時域信號重構(gòu)，結(jié)果如圖8所示。經(jīng)計算，重構(gòu)的信號與原始信號的平均相對誤差為1.28%，說明所提方法擁有信號的重構(gòu)性能[16]。

圖7 多分量信號角度估計

圖8 多分量信號的重構(gòu)

為比較所提方法對噪聲的抗干擾能力，對圖7中所分析信號添加適量噪聲以進一步分析。含噪信號的信噪比從-5 dB變化至20 dB，估計角度的平均相對誤差如圖9所示。從圖9中可以看出，隨著信噪比降低，誤差也逐漸增大，說明了噪聲對角度估計的準確性有一定的影響。從圖9可以看出，在分析信噪比低為-5 dB的信號時，所提方法仍能比較準確地估計角度，計算的誤差均小于5%。

圖9 所估計角度的平均誤差隨輸入信噪比變化

1.4 窗長的確定

上述分析中，單分量和多分量信號分析結(jié)果均是采用固定長度的窗函數(shù)分析所得到的，考慮到窗長的變化對時頻表示的分辨率也會有較大影響，故本節(jié)進一步分析窗長變化對最終時頻表示的影響。同樣地，根據(jù)峭度最大準則來選取每個時刻合適的分析窗長。同時考慮窗長和角度變化的時頻表示記為XG[τ,ω,g(τ),α(τ)]。根據(jù)式(24)，在同時考慮窗長與角度變化情況下最大峭度對應的參數(shù)選擇可表示為

[g*(τ),α*(τ)]=

(29)

式中，g*(τ)和α*(τ)分別為τ時刻峭度指導選取的窗長和角度。再次分析1.3節(jié)中的多分量信號，結(jié)果如圖10所示。圖10(a)給出了所提方法在t=1.155 s處角度、窗長與峭度的對應關(guān)系，并用三角標出了計算得到的峭度最大值，通過該點的坐標可知與最大峭度所對應的角度與分析窗長。依據(jù)式(29)，所分析信號在整個時段選擇的角度和窗長如圖10(b)和圖10(c)所示。從圖中可以看出所選取角度基本與真實值吻合，經(jīng)計算平均相對誤差為1.24%。在圖10(c)中，分析窗長在前2 s內(nèi)經(jīng)歷了較多的波動，但仍可以看出，在頻率變化相對緩慢時，可選用相對較長的窗來分析信號，比如0～2 s左右的窗相對3 s更長(3 s左右信號頻率變化相對0～2 s更快)。最終結(jié)果如圖10(d)所示，與圖7(c)獲得的時頻表示(固定窗長下的多分量分析結(jié)果)相比，可以發(fā)現(xiàn)同時考慮角度和窗長可得到更加集中的時頻脊線，也說明了選取合適的窗長可進一步增強時頻表示。

圖10 所提方法中角度和窗長的選擇

2 軸承仿真信號分析

實際采集的信號不可避免會包含噪聲，并且噪聲也會對峭度指導的參數(shù)選擇產(chǎn)生影響。因此本章中在添加噪聲后進一步驗證所提方法的有效性并將結(jié)果與常見時頻分析方法所得結(jié)果進行對比。

考慮旋轉(zhuǎn)機械振動信號常包含多個分量，定義的含局部故障的軸承仿真多分量信號為

x(t)=

(30)

式中：fISRF為仿真信號的軸轉(zhuǎn)頻；n(t)為噪聲。

該仿真信號的采樣頻率設為8 000 Hz，共振頻率設為3 500 Hz，信號持續(xù)1 s。設定的故障特征系數(shù)為2.7,即fIFCF=2.7fISRF。

添加適量噪聲，將信噪比設定為3 dB進行分析。含噪信號的時域波形如圖11(a)所示，用所提方法對信號進行分析，在t=0.5 s角度、窗長與峭度值的對應關(guān)系，如圖11(b)所示，圖11(b)中三角代表該時刻計算得到的最大峭度值及其對應的角度和窗長。同樣，對整個信號分析角度和窗長變化的影響，得到的結(jié)果分別如圖11(c)和圖11(d)所示。圖11(c)中，所估計角度的平均相對誤差為1.62%。在圖11(d)中可看出，前0.5 s選取的窗長比在1.5～2.0 s左右更長，這與時頻分析方法對窗長的要求基本吻合。

圖11 含噪多分量信號分析

再次采用不同時頻分析方法對信號進行分析，包括STFT，廣義線性調(diào)頻變換(general linear chirplet transform, GLCT)[17]，多項式調(diào)頻變換(polynomial chirplet transform, PCT)[18], SST,SET和所提方法，得到的時頻表示如圖12所示。在分析該仿真信號時，所有時頻方法的窗長統(tǒng)一設為0.4 s以方便對比。圖12(a)中所示STFT結(jié)果僅在第1階轉(zhuǎn)頻fISRF以及t=1.7 s左右有較好的時頻聚集性。GLCT和PCT方法的分析結(jié)果如圖12(b)和圖12(c)所示, GLCT算法中所用調(diào)頻率在結(jié)果中導致了較為嚴重的時頻交叉干擾，從而難以獲得準確的瞬時頻率估計且難以觀測得到2倍轉(zhuǎn)頻對應的頻率分量存在；而PCT因采用固定的線性變換基函數(shù)，適合處理具有相同變化趨勢的信號。圖12(c)中，基函數(shù)設置與fISRF一致，因此fIFCF及其諧波分量的能量相對分散。此外,PCT還需依賴轉(zhuǎn)頻相關(guān)的先驗知識。

SST和SET時頻后處理方法的分析結(jié)果如圖12(d)和圖12(e)所示。SST因僅沿頻率方向壓縮，受快變頻率和噪聲的影響，得到的結(jié)果在0.5～0.8 s遭受了較為嚴重的時頻模糊問題，難以觀測得到清晰的瞬時頻率脊線，如2fIFCF；而SET通過類似脈沖函數(shù)的同步提取因子，其結(jié)果為提取因子與STFT時頻表示的乘積，即在STFT所得時頻表示的基礎上提取目標頻率分量(瞬時頻率估計因子)及其周圍的時頻分布，因此不可避免地受到STFT原始時頻圖中時頻模糊問題的影響，只能提供相比圖12(a)更加銳利的時頻脊線，而無法提升能量的聚集性。所提方法的結(jié)果如圖12(f)所示，從圖12(f)中可清晰觀測得到fISRF和fIFCF的倍頻及其諧波存在。雖然得到的結(jié)果中在某些時刻的時頻脊線不如后處理方法(見圖12(d)和圖12(e))集中，但時頻表示的能量聚集性相對更高，且不需要瞬時頻率相關(guān)的先驗知識。因此，所提方法與其他方法相比具有一定的優(yōu)越性，可有效用于提升時頻表示可讀性并具有一定的抗噪性能。

圖12 不同時頻分析方法所得時頻表示(含噪多分量信號)

3 旋轉(zhuǎn)機械故障診斷的應用

變轉(zhuǎn)速工況下旋轉(zhuǎn)機械振動信號中同樣含有多個頻率分量，如旋轉(zhuǎn)頻率、故障特征頻率及噪聲干擾等。相比傳統(tǒng)的廣義解調(diào)方法在分析時需要考慮迭代，所提方法因與轉(zhuǎn)速同步，可同時增強多個頻率分量的能量聚集性從而簡化算法。在獲得時頻表示后，根據(jù)時頻圖中的瞬時頻率脊線可診斷出該旋轉(zhuǎn)機械系統(tǒng)中是否存在故障及故障類型。

3.1 行星齒輪信號分析

將所提方法用于分析試驗信號，該行星齒輪信號的故障試驗臺如圖13所示。具體配置如表1所示，根據(jù)齒輪箱的配置可計算出各自的故障特征頻率系數(shù)。

表1 齒輪箱的參數(shù)配置

1.電機;2.測速計;3.定軸齒輪變速箱;4.行星齒輪箱一級傳動軸;5.行星齒輪箱二級傳動軸;6.加速度計;7.制動器。

該旋轉(zhuǎn)機械系統(tǒng)的轉(zhuǎn)速fd是由交流電機驅(qū)動控制的，轉(zhuǎn)速先增加然后降低。采集信號的時域波形如圖14(a)所示，該信號的采樣頻率設為20 kHz，信號采集時間持續(xù)30 s。轉(zhuǎn)頻fd的變化如圖14(b)所示，先從40 Hz變化到60 Hz，再變化到40 Hz。測速計采集的轉(zhuǎn)速信息主要用來驗證所提角度選擇方法是否準確。對該信號進行頻譜分析，0～400 Hz內(nèi)的頻譜如圖14(c)所示，由于該信號是在變轉(zhuǎn)速工況下采集得到，因此在頻譜中觀測不到定轉(zhuǎn)速工況下特別突出的特征頻率等峰值。通常試驗信號包含故障所引起的共振頻段，為避免失真，采樣頻率設置較大；而當試驗信號采集時間也較長時，分析中若考慮窗長變化會極大地增加計算量。因此，采用固定窗長來比較所獲得時頻表示的能量聚集性。

取時刻t=20 s分析峭度隨所選取的角度的變化規(guī)律，得到的結(jié)果如圖14(d)所示，圖14(d)中用三角標記了最大峭度及其對應角度。同步估計所有時刻的角度，結(jié)果如圖14(e)所示，其中，虛線代表由測速計采集信號直接計算得到的真實角度，實線為峭度指導選取的角度。通過對比可發(fā)現(xiàn)，所選取角度的變化趨勢基本與真實情況吻合。STFT,GLCT和所提方法的時頻分析結(jié)果分別如圖14(f)～圖14(h)所示。對比時頻分析結(jié)果表明：采用相同窗長分析時，所提方法由于找到了合適的角度來同步解調(diào)信號，可得到更高的能量聚集性和更集中的時頻脊線，如時頻表示中轉(zhuǎn)速二倍頻2fd處。GLCT也是通過匹配信號調(diào)頻率來增強時頻表示，但得到的時頻分析結(jié)果同時受窗長和不準確的調(diào)頻率影響，從而導致時頻模糊問題并降低時頻表示可讀性。從圖14(h)中，可清晰地觀測到轉(zhuǎn)頻fd及其前5階諧波分量和第一級齒輪嚙合的頻率分量fmesh1(fmesh1=100/27fd)，而并未發(fā)現(xiàn)與該行星齒輪故障相關(guān)的邊頻帶存在，這些頻率分量可說明該行星齒輪是健康的。

圖14 行星齒輪箱振動信號分析

3.2 軸承內(nèi)圈故障診斷

同樣將所提方法用于實際軸承振動信號分析[19]。該軸承故障模擬試驗臺的示意圖，如圖15所示。該系統(tǒng)中包括兩個軸承，右邊的軸承存在內(nèi)圈故障。信號采樣頻率設為20 kHz，持續(xù)10 s。

圖15 軸承內(nèi)圈故障試驗臺

采集得到的振動信號如圖16(a)所示，該信號的瞬時旋轉(zhuǎn)頻率(電機轉(zhuǎn)頻)fISRF變化如圖16(b)所示，這里采集的fISRF主要用于計算真實角度來驗證所提方法根據(jù)最大峭度選取角度的準確性。該軸承的內(nèi)圈故障特征系數(shù)為5.43，即故障特征頻率fIFCF=5.43fISRF。選取時刻t=4 s分析峭度隨角度的變化規(guī)律，結(jié)果如圖16(c)所示，圖16(c)中用三角標記出了最大峭度值和對應的角度。對整個信號分析的結(jié)果如圖16(d)所示?？梢园l(fā)現(xiàn)，根據(jù)最大峭度選擇的角度基本是與真實值一致的，經(jīng)計算所估計角度的平均相對誤差為1.19%。采用STFT以及SST作對比，分析結(jié)果如圖17(a)和圖17(b)所示。所提方法結(jié)果如圖17(c)所示。經(jīng)對比可以發(fā)現(xiàn)所提方法能夠得到能量更加集中的時頻脊線并表征出信號中微弱的頻率特征成分如fIFCF+fISRF,2fIFCF+fISRF和3fIFCF-2fISRF。因所提方法同時增強多個頻率分量的時頻聚集性，可得到更精確的時頻脊線，如圖17(d)所示。根據(jù)搜尋得到的時頻脊線以及存在的轉(zhuǎn)頻調(diào)制現(xiàn)象均可表明此試驗軸承存在內(nèi)圈故障[20]。因此所提方法可成功用于診斷旋轉(zhuǎn)機械系統(tǒng)潛在的健康狀態(tài)。

圖16 軸承振動信號的角度估計

圖17 不同方法對軸承振動信號的分析結(jié)果

4 結(jié) 論

本文針對變轉(zhuǎn)速工況下旋轉(zhuǎn)機械非平穩(wěn)信號處理存在時頻模糊和脊線提取困難等問題，基于廣義解調(diào)提出了廣義瞬時速度同步化分步解調(diào)變換，用以對非平穩(wěn)信號展開分析，并成功應用于變轉(zhuǎn)速下旋轉(zhuǎn)機械的健康狀態(tài)監(jiān)測。首先，引入角度變量并計算了廣義分步解調(diào)變換中的前向和后向解調(diào)因子；然后，為減少對瞬時頻率預估計的依賴、實現(xiàn)瞬時速度同步化，提出了峭度指導的角度(瞬時頻率)自適應選擇策略，可避免人為干預和對先驗知識的依賴；最后，通過新的解調(diào)因子擴展原先的線性變換基函數(shù)，使所提方法能同時對多個分量進行解調(diào)并同步增強多個頻率分量在時頻表示中的能量聚集性。

針對信號重構(gòu)，展開了證明，并對仿真信號進行了重構(gòu)分析，結(jié)果表明：所提時頻分析方法可根據(jù)提取的目標時頻脊線重構(gòu)出原始時域信號。此外，對窗長變化對時頻分析的影響也展開了探討，說明動態(tài)窗長可進一步提升信號的時頻表示。仿真分析和試驗驗證說明該方法能得到更加集中的時頻脊線和更高的時頻能量聚集性從而提升時頻表示可讀性。對增強后的時頻表示進行分析，可更準確地提取出變轉(zhuǎn)速工況下旋轉(zhuǎn)機械振動信號中與故障相關(guān)的特征，從而實現(xiàn)對變轉(zhuǎn)速旋轉(zhuǎn)機械的故障診斷。