王簫劍， 涂曉彤， 李鴻光， 李富才， 包文杰

(上海交通大學(xué)機(jī)械系統(tǒng)與振動國家重點實驗室 上海，200240)

引 言

行星齒輪箱是由太陽輪、行星輪、大齒圈以及行星架等元件組成的機(jī)構(gòu)，具有降低轉(zhuǎn)速比、放大電機(jī)扭轉(zhuǎn)力等功能。因其傳遞平穩(wěn)、承載力大、體積小巧等優(yōu)點，廣泛應(yīng)用于船舶、航天、汽車等領(lǐng)域中。然而此類行業(yè)對行星齒輪箱的安全平穩(wěn)要求極高，所以對行星齒輪箱的工作情況進(jìn)行實時監(jiān)測是十分必要的。目前，越來越多的行星齒輪箱在變轉(zhuǎn)速的情況下工作，比如汽車的升降速、飛機(jī)的起落過程等等。傳統(tǒng)的頻域信號處理方法僅能用于平穩(wěn)信號的分析，所以一種既能定量、又能直觀地描述信號故障特征的非平穩(wěn)信號處理方法是極其必要的。時頻分析對于提取非平穩(wěn)信號的時變特性很有優(yōu)勢。通常使用一些傳統(tǒng)的方法，例如連續(xù)小波變換或者短時傅里葉變換，在時頻平面上呈現(xiàn)非平穩(wěn)信號特征[1-2]。但是，由于不確定性定理，它們并不能很好地保證時頻平面的分辨率。為了處理這一問題，匹配解調(diào)技術(shù)[3-5]和匹配同步壓縮技術(shù)[6-7]應(yīng)運而生，然而，每一時刻瞬時頻率的預(yù)提取仍然是目前面臨的問題。對于這些瞬時頻率(時頻圖譜脊線)提取問題，國內(nèi)外研究很多[8-10]，但在強(qiáng)噪聲、多分量情況下，因噪聲、無關(guān)分量對目標(biāo)分量的干擾，對于弱分量提取較為困難。

筆者提出了匹配壓縮脊線提取(demodulated synchrosqueezing ridge extraction，簡稱DSRE)，該方法適用于處理強(qiáng)噪聲、多分量信號，并可應(yīng)用于變轉(zhuǎn)速旋轉(zhuǎn)機(jī)械故障診斷與在線監(jiān)測。DSRE可分為兩個步驟：a.根據(jù)時頻圖進(jìn)行脊線提?。籦.根據(jù)所提脊線，進(jìn)行匹配解調(diào)變換和二階同步壓縮得到新的時頻圖和新的脊線。為了處理多分量信號，提出了時變?yōu)V波方法，用于分量剔除與提取，得到每一階分量的瞬時頻率。將該方法應(yīng)用于行星齒輪箱的故障診斷中，能夠直觀地提取基頻及故障頻率的特征，便于行星齒輪箱實時監(jiān)測與故障診斷。

1 DSRE基本原理

1.1 脊線提取技術(shù)與包絡(luò)線濾波技術(shù)

對于非平穩(wěn)信號，由于每一時刻的瞬時頻率不同，需要對其進(jìn)行時頻變換，例如小波變換、短時傅里葉變換等。

對于任意一個信號x(t)=A(t)ei2π(f0t+φ(t))，f0為其載波頻率，φ(t) 為其瞬時相位角，對應(yīng)φ′(t) 為其瞬時頻率。其短時傅里葉變換可以表示為

(1)

通常選用標(biāo)準(zhǔn)差為σ的高斯窗函數(shù)

(2)

由于短時傅里葉變換的時域窗長度保持恒定，因而無法有效處理時變信號。

在時頻圖上每一時刻具有最大能量的點的序列，被定義為脊線。脊線的縱坐標(biāo)，對應(yīng)著信號的瞬時頻率。

對于脊線提取的方式有很多種[8-10]，筆者用于脊線提取的過程[8]如下。

1) 對信號經(jīng)過時頻變換得到的圖譜，計算每一時間點幅值旳極值Qm(t)及其所對應(yīng)的瞬時頻率νm(t)

(3)

其中：TFD(t,ω)為t時間點ω頻率處對應(yīng)的時頻變換幅值；Np(t)對應(yīng)t時刻的最大極值點個數(shù)；m為t時刻對應(yīng)的第m個最大極值點。

2) 定義路徑優(yōu)化函數(shù)

{mc(τ1),mc(τ2),…,mc(τN)}=

{νm1(τ1),…,νmN(τN)]

(4)

F由當(dāng)前時刻的時頻譜幅值Qm(t)(當(dāng)前時刻的能量)，及當(dāng)前時刻的頻率與上一時刻的頻率之差νm(tn)-νk(tn-1) (時頻圖譜脊線的連續(xù)性)共同決定

F[Qm(tn),νm(tn),νk(tn-1)]=

logQm(tn)+w1(νm(tn)-νk(tn-1))

(m∈{1,2,…,NP(tn)}；k∈{1,2,…，NP(tn-1)})

(5)

其中：w為衡量脊線連續(xù)性的因子。

w1(Δξ)=-σ2fs|Δξ|

(6)

3) 對上述路徑優(yōu)化函數(shù)進(jìn)行求解，得到的設(shè)計變量即為每一時刻點對應(yīng)的瞬時頻率。

此路徑優(yōu)化函數(shù)的優(yōu)點是既考慮了脊線當(dāng)前點的能量，又考慮了脊線的連續(xù)性。在處理含噪聲的信號時，從脊線連續(xù)性角度，可以有效地消除噪點對所提信號信息的干擾。

對于多分量信號，可經(jīng)過包絡(luò)線濾波技術(shù)去除所提脊線分量信號。在時頻圖上，找到對應(yīng)所提脊線的上下包絡(luò)線，將上下包絡(luò)線中間部分幅值均重置為零。

1.2 旋轉(zhuǎn)算子濾波-匹配解調(diào)技術(shù)

多項式匹配解調(diào)算法的思路是在處理信號每個時間點鄰域的時候，消除信號脊線在時間點鄰域的高次項，從而在每一個時間點的鄰域，頻率都不隨時間改變，即在所有高斯窗里，信號頻率均時不變。

將上文所提取的脊線，作為瞬時頻率參考，據(jù)此提出了旋轉(zhuǎn)算子dR(t)和平移算子dS(t,u)

其中：φ′(t)為信號的瞬時頻率，由脊線經(jīng)多項式擬合得到；u為t時刻鄰域?qū)?yīng)時間點。

旋轉(zhuǎn)算子的作用是去除整個時域內(nèi)瞬時頻率脊線的1階及其以上導(dǎo)數(shù)，而平移算子的作用是將每個時間點鄰域賦予恒等同于該時間點的瞬時頻率，如圖1所示。

圖1 匹配解調(diào)過程旋轉(zhuǎn)變換及平移變換Fig.1 Rotating step and shifting step in the moduling process

經(jīng)過旋轉(zhuǎn)平移變換之后，信號Sd(t,u)為

SdR(t)=x(t)dR(t)=A(t)ei2πf0t

(9)

Sd(t,u)=sdR(t)dS(t)=A(t)ei2π[f0t+φ′(u)t]

(10)

此信號經(jīng)過短時傅里葉變換之后，即得到較為精確的時頻圖。

匹配解調(diào)方法的表達(dá)式為

(11)

對于多分量信號，可通過旋轉(zhuǎn)算子濾波得到目標(biāo)分量。該分量脊線提取旋轉(zhuǎn)算子后，對信號旋轉(zhuǎn)變換，信號將近似于時不變信號，此時濾波，即可在時域信號上去除其他分量和噪聲，如圖2所示。

圖2 匹配解調(diào)-旋轉(zhuǎn)算子濾波Fig.2 The filtering step in the moduling process

對濾波之后的信號進(jìn)行平移變換和短時傅里葉變換，即可得到所需信號的時頻圖。

1.3 2階匹配同步壓縮技術(shù)

匹配同步壓縮技術(shù)就是把信號的各時頻點在信號的瞬時頻率處進(jìn)行重排，從而使能量匯集在信號的瞬時頻率附近。本研究匹配同步壓縮的目的是使匹配解調(diào)所得的時頻變換圖譜上面的能量更集中，從而有利于接下來脊線的提取。

時頻信號的壓縮方法如下。

以一個信號的短時傅里葉變換為例

(12)

首先，將此短時傅里葉變換經(jīng)過從u-t到v的換元可得

(13)

對于一個恒定幅值的信號，將信號瞬時相位泰勒級數(shù)展開，并保留二次項系數(shù)

x(v+t)=Ae(a+b(t+v)+0.5c(t+v)2)i

(14)

其中：a，b，c分別為泰勒展開常數(shù)項、一次項、二次項系數(shù)。

再假設(shè)信號頻率1階導(dǎo)數(shù)遠(yuǎn)大于2階導(dǎo)數(shù)的情況下，將式(14)經(jīng)過對時間的求導(dǎo)，可得

R(?tSTFT(t,w))=

(15)

傅里葉變換的導(dǎo)數(shù)可以表示為

R(?tSTFT(t,w))=i(b+ct)STFT(t,ω)+

icTSTFT(t,ω)

對比上面兩式發(fā)現(xiàn)

(16)

然后，將短時傅里葉圖譜對于式(16)解得的瞬時頻率再賦值，將每一時間點全部能量聚集到瞬時頻率周圍，可以得到壓縮之后的圖譜

(17)

其中：Ssd(t,f)為1.2節(jié)變換得到的時頻分布；g*(0)為信號能量加權(quán)因子；γ為能量閾值。

1.4 多分量信號提取故障分量方法實施

對于多分量信號，首先，提取能量最大信號分量的脊線，時頻變換得到所需信息后，去除該分量，并提取下一階信號分量；其次，依次提取余下信號分量。具體步驟如下：

1) 進(jìn)行短時傅里葉變換，用脊線提取法估計出能量最大分量的脊線，同時采用包絡(luò)法估計脊線的上下邊界；

2) 用多項式擬合所得脊線，得到多項式各參數(shù)，然后采用匹配解調(diào)技術(shù)，經(jīng)過旋轉(zhuǎn)算子變換后濾波，并進(jìn)行平移變換得到新的時頻圖；

3) 根據(jù)得到的瞬時頻率進(jìn)行2階同步壓縮，得到該分量的時頻圖，并進(jìn)行第2次脊線提取提出修正之后的脊線及其上下邊界；

4) 在時頻圖上根據(jù)脊線及其上下邊界去除能量最大的分量，重復(fù)步驟1～3，直到提取出故障信號。

對每一階信號分量的提取，該法的計算代價為O(3NfN)+O(M2N)，其中：N為信號時域采樣點；Nf為頻域點數(shù)目；M為1.1節(jié)中Np(t)最大值。

2 仿真分析

2.1 仿真信號構(gòu)造

為了檢驗該信號處理方法，采用含有高斯白噪聲的非平穩(wěn)仿真信號，仿真信號f(t)如下

(18)

該仿真信號的采樣頻率為1 024 Hz，由3種瞬時頻率隨時間改變的信號分量疊加而成，同時加入高斯白噪聲η(t)，該白噪聲平均值為0，信噪比為-1 dB。時域波形如圖3(a)所示，頻譜圖如圖3(b)所示。

圖3 仿真信號特征Fig.3 Characteristics of the simulation signal

2.2 能量最大信號分量提取

對信號進(jìn)行時頻變換并提取脊線，所得結(jié)果如圖4(a)所示。圖4(a)為采用短時傅里葉變換得到的時頻圖，紅線代表采用文獻(xiàn)[8]提取的脊線。如圖所示，所提脊線雖然連續(xù)，并能反應(yīng)信號特征，但由于噪聲的存在，脊線存在一定波動，和真實情況相比具有誤差。

對信號進(jìn)行匹配解調(diào)旋轉(zhuǎn)變換，得到時頻圖見圖4(b)。如圖所示，信號經(jīng)過匹配解調(diào)的旋轉(zhuǎn)變換后，其第1階分量的瞬時頻率基本保持不變。

對所得信號進(jìn)行傅里葉變換及帶通濾波，濾去噪聲信號和其他階分量信號，如圖5所示。

為了驗證所提脊線的準(zhǔn)確性，比較本研究方法所提脊線的精確程度與文獻(xiàn)[8]中所述方法的精確程度，如圖6所示。其中：點劃線為本研究方法所提脊線；虛線為文獻(xiàn)[8]所提脊線；實線為真實脊線。由圖可以看出，本研究算法所提脊線與真實脊線基本重合，而文獻(xiàn)[8]中所提脊線在真實脊線附近存在波動。此波動誤差會隨著所提分量的階數(shù)的增高而逐漸積累，從而影響弱能量分量提取的精度。

圖4 旋轉(zhuǎn)變換前后時頻圖譜Fig.4 Time-frequency transform before and after the rotating process

圖5 旋轉(zhuǎn)變換信號濾波前后傅里葉變換圖譜Fig.5 The Fourier transform of the rotated signal before and after the filtering process

圖6 第1階分量本研究所提脊線與真實脊線對比Fig.6 Comparation between the extracted 1 st ridge curves and true 1 st ridge curve

筆者采用范數(shù)來定量表現(xiàn)所提脊線精確程度

(19)

其中：EIF為估計瞬時頻率；TIF為真實瞬時頻率。

經(jīng)計算可得本研究算法所得誤差為0.001 0，而文獻(xiàn)[8]所述算法誤差為0.003 0，故本算法所提脊線精度比文獻(xiàn)[8]脊線準(zhǔn)確。

2.3 能量次大信號分量提取

在提取出第1階信號瞬時頻率之后，在時頻圖上將第1階信號頻率成分濾去，并進(jìn)行第2次脊線提取，如圖7所示。

圖7 去除1階信號分量后時頻變換圖譜及脊線預(yù)提取Fig.7 The time-frequency transform and previous ridge extraction after removing the first signal

對于該多分量信號，圖8為第2階分量的脊線對比圖。其中：點劃線為本研究所提脊線；虛線為采用文獻(xiàn)[8]并結(jié)合本研究的脊線流程所提脊線；實線為真實脊線。由圖可以看出，本研究所提脊線擬合程度更高。

對于該脊線，本算法范數(shù)為0.019 1，文獻(xiàn)[8]算法范數(shù)為0.048 4，故本算法精確程度更高。

圖 8 第2階信號分量所提脊線與真實脊線對比Fig.8 Comparation between the extracted 2st ridge curves and true 2st ridge curve

3 行星齒輪箱故障診斷

本試驗的振動數(shù)據(jù)采集自行星齒輪箱故障模擬試驗臺。如圖9所示，試驗臺包括驅(qū)動電機(jī)、交流電機(jī)、行星齒輪箱、固定軸齒輪箱和制動器。采用轉(zhuǎn)速傳感器和加速度傳感器分別采集電機(jī)軸轉(zhuǎn)速和行星齒輪箱振動信號，電機(jī)軸轉(zhuǎn)速可由驅(qū)動電機(jī)自由調(diào)節(jié)。本試驗中電機(jī)為變轉(zhuǎn)速，包括一段升速(1 800～2 700 r/min)和一段減速過程(2 700～1 800 r/min)。信號采樣頻率為12 800 Hz，采集時長為16 s。

圖9 行星齒輪試驗臺Fig.9 The test rig of the planetary gearbox

表1為行星齒輪箱的齒輪參數(shù)，表2為行星齒輪箱內(nèi)各零件的特征階次。

表1 行星齒輪箱齒輪參數(shù)

表2 行星齒輪箱故障階次

Tab.2 Fault orders of the planetary gearbox

10.1672.500.83

對振動信號原始數(shù)據(jù)時頻變換，信號的時域波形如圖10(a)所示，時頻變換圖譜如圖10(b)所示。由于外界噪聲和信號多重分量的疊加，難以從時域信號或時頻變換圖像中發(fā)現(xiàn)故障特征。信號能量最集中的分量是齒輪加速度信號的三倍頻分量，即使采用傳統(tǒng)方法提取能量最集中的信號分量，再進(jìn)行階次分析的方式也難以實現(xiàn)。

圖10 行星齒輪箱振動信號特征Fig.10 Characteristics of the acceleration signal of vibration of a planetary gearbox

首先，采用本方法提取基頻。在這個信號中，基頻并不是能量最大的分量，能量最大的分量是三倍頻，所以先提取信號的三倍頻，如圖11(a)所示。提取三倍頻并濾去三倍頻信號分量之后，即可提取信號一倍頻特征頻率，如圖11(b)所示。

以此類推，將倍頻成分全部去除后，繼續(xù)進(jìn)行脊線提取，這樣即可提取出故障頻率，如圖12所示。實線為采用文獻(xiàn)[8]方法并結(jié)合本研究提脊線流程所提脊線；虛線為本研究提取脊線。由圖可知，本研究所提脊線方法可行，且較為準(zhǔn)確。

將所提故障分量信號與基頻信號分量作商，即可得到故障信號分量，如圖13(a)所示。由圖可知，故障頻率為2.5階，為太陽輪故障。將所得故障信號分量與基頻信號分量于一張時頻圖上呈現(xiàn)，亦可直觀地表現(xiàn)故障信號的變化過程及其強(qiáng)弱，如圖13(b)所示，其中：上方線為故障分量；下方線為基頻分量。若時間足夠長，可通過時頻圖的明暗觀察故障的變化情況。

圖11 行星齒輪箱各階振動信號分量提取Fig.11 Ridge extraction of vibration signals of different orders

圖12 故障頻率脊線提取Fig.12 Ridge extraction of vibration signals of failure

圖13 行星齒輪箱振動信號故障分量提取Fig.13 Characteristic extraction of fault component of the vibration signal

4 結(jié) 論

1) DSRE用于對于含噪聲信號的處理，采用脊線提取與匹配解調(diào)技術(shù)相結(jié)合的方式，比一般直接脊線提取精度更高。

2) 對于多分量信號的處理，提出旋轉(zhuǎn)算子濾波方法，去除了其余分量信號與噪聲信號對所提信號的干擾。

3) 用于齒輪等具有階次頻譜的變轉(zhuǎn)速旋轉(zhuǎn)機(jī)械故障診斷中，可以將故障所在階次與其能量反映出來，監(jiān)測故障階次與其能量幅值的變化過程，以便于對后續(xù)故障機(jī)理的研究，可對故障的實時監(jiān)測提供參考。