匹配壓縮脊線提取在齒輪箱故障診斷中的應用*

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

(上海交通大學機械系統與振動國家重點實驗室 上海，200240)

引 言

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

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

1 DSRE基本原理

1.1 脊線提取技術與包絡線濾波技術

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

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

(1)

通常選用標準差為σ的高斯窗函數

(2)

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

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

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

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

(3)

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

2) 定義路徑優化函數

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

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

(4)

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

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為衡量脊線連續性的因子。

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

(6)

3) 對上述路徑優化函數進行求解，得到的設計變量即為每一時刻點對應的瞬時頻率。

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

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

1.2 旋轉算子濾波-匹配解調技術

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

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

其中：φ′(t)為信號的瞬時頻率，由脊線經多項式擬合得到；u為t時刻鄰域對應時間點。

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

圖1 匹配解調過程旋轉變換及平移變換Fig.1 Rotating step and shifting step in the moduling process

經過旋轉平移變換之后，信號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)

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

匹配解調方法的表達式為

(11)

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

圖2 匹配解調-旋轉算子濾波Fig.2 The filtering step in the moduling process

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

1.3 2階匹配同步壓縮技術

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

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

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

(12)

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

(13)

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

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

(14)

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

再假設信號頻率1階導數遠大于2階導數的情況下，將式(14)經過對時間的求導，可得

R(?tSTFT(t,w))=

(15)

傅里葉變換的導數可以表示為

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

icTSTFT(t,ω)

對比上面兩式發現

(16)

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

(17)

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

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

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

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

2) 用多項式擬合所得脊線，得到多項式各參數，然后采用匹配解調技術，經過旋轉算子變換后濾波，并進行平移變換得到新的時頻圖；

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

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

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

2 仿真分析

2.1 仿真信號構造

為了檢驗該信號處理方法，采用含有高斯白噪聲的非平穩仿真信號，仿真信號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 能量最大信號分量提取

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

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

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

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

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

圖5 旋轉變換信號濾波前后傅里葉變換圖譜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

筆者采用范數來定量表現所提脊線精確程度

(19)

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

經計算可得本研究算法所得誤差為0.001 0，而文獻[8]所述算法誤差為0.003 0，故本算法所提脊線精度比文獻[8]脊線準確。

2.3 能量次大信號分量提取

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

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

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

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

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

3 行星齒輪箱故障診斷

本試驗的振動數據采集自行星齒輪箱故障模擬試驗臺。如圖9所示，試驗臺包括驅動電機、交流電機、行星齒輪箱、固定軸齒輪箱和制動器。采用轉速傳感器和加速度傳感器分別采集電機軸轉速和行星齒輪箱振動信號，電機軸轉速可由驅動電機自由調節。本試驗中電機為變轉速，包括一段升速(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為行星齒輪箱的齒輪參數，表2為行星齒輪箱內各零件的特征階次。

表1 行星齒輪箱齒輪參數

表2 行星齒輪箱故障階次

Tab.2 Fault orders of the planetary gearbox

10.1672.500.83

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

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

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

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

將所提故障分量信號與基頻信號分量作商，即可得到故障信號分量，如圖13(a)所示。由圖可知，故障頻率為2.5階，為太陽輪故障。將所得故障信號分量與基頻信號分量于一張時頻圖上呈現，亦可直觀地表現故障信號的變化過程及其強弱，如圖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 結 論

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

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

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