基于自適應奇異值分解的行星齒輪箱故障診斷方法

秦 毅， 張清亮， 趙 月

(1. 重慶大學 機械傳動國家重點實驗室， 重慶 400044； 2. 重慶大學 機械工程學院， 重慶 400044)

行星齒輪箱具有結構緊湊、體積小、質量輕、承載能力強、傳動比大、傳動效率高等優點[1]，已被廣泛應用于諸多機械傳動系統。行星齒輪箱大部分使用在低速重載的惡劣工況，太陽輪、行星輪、內齒圈和行星架等部件很容易出現故障。對行星齒輪箱進行故障診斷可避免因事故造成的巨大經濟損失，具有重大的意義。

目前，利用振動信號進行診斷是故障診斷技術中最有效的方法之一。設備的故障振動信號包含尖峰、突變的非平穩信號，還可能包含非平穩的白噪聲，這些是信號的高頻部分[2]。受到設備運行環境、電磁干擾和采樣儀器的限制，采集到的數據中不可避免地含有噪聲。由于振動信號中包含的故障信號一般較弱，常淹沒在噪聲信號中，所以必須對信號進行消噪處理，需要在消除噪聲所表現的高頻量的同時，保留反映信號突變部分的高頻量。為了去除振動信號中的噪聲，目前已經提出了很多方法，如自適應濾波[3-4]、時域同步平均[5-6]、EEMD分解[7-8]、小波閥值降噪[9-10]等。自適應濾波雖然只需要含噪信號就能進行，但是要根據具體信號來確定延遲時間、濾波器階數和收斂因子。時域同步平均需要將信號分成N段，然后將這N段信號相加再除以N。雖然其降噪效果隨著N的增大而增加，但是也會減小信號的長度從而減小分辨率。EEMD分解是EMD分解的改進，簡單高效、高分辨率、能夠抑制模式混疊，但是分解時加入的高斯白噪聲標準差只有一個范圍而沒有一定的計算公式。小波閥值降噪雖然實現簡單、計算量小，但是小波基、分解層數、閥值等參數選取上存在很大不確定性。

奇異值分解(Singular Value Decomposition, SVD)是一種非線性濾波。其消噪具有小相移[11]，無時延等特點，被廣泛應用于振動[12]信號消噪。由于實際振動信號中信噪比難以估計，因此峰值信噪比[13]很難使用。而奇異值差分譜[14]能自適應地選擇用于重構的奇異值個數。但是僅選取奇異值差分譜最大峰值點對應奇異值數目重構會丟失一些弱沖擊信號特征。峰值群[15]的提出是最大奇異值分解的改進使其克服了會丟失一些弱沖擊信號特征的缺點，但是該方法要求根據觀察奇異值差分譜的分布來確定奇異值數目，不能做到自動選擇。本文在現有奇異值分解的基礎上，提出了一種能夠自動選擇有效奇異值數目的自適應奇異值分解(Adaptive Singular Value Decomposition，ASVD)方法，以便更加準確地提取行星齒輪箱的微弱故障特征。

1 奇異值分解理論

1.1 奇異值分解

設離散數字信號X=[x(1),x(2),…,x(N)]，則利用此信號可以構造Hankel矩陣如下

(1)

式中:N=m+n-1,m≤n,A∈Rm×n且矩陣A的反對角線上的元素相同。由文獻[14]知，Hankel矩陣A的行數和列數最好相等或者相差不超過1。因此，當N為偶數時，應取m=N/2、n=N/2+1；當N為奇數時，應取m=(N+1)/2、n=(N+1)/2。

Hankel矩陣A的奇異值分解可定義為

A=U∑VT

(2)

重構信號就是要選擇Hankel矩陣A的奇異值之中的前面k個有效奇異值并將之后r-k的個奇異值置零，并進行奇異值分解逆變換，即

A*=U∑*VT

(3)

由于矩陣A*的反對角上的元素不同，所以需要先進行反對角線平均[16]，從而轉化成Hankel矩陣，之后取反對角線平均之后矩陣的第一行和消去第一行元素的最后一列作為重構信號，其計算公式如下

x*(i)=

(4)

1.2 奇異值差分譜與最大差分奇異值分解

含噪聲的離散數字信號x(i)可表示為

x(i)=s(i)+ε(i)i=1,2,…,N

(5)

式中:s(i)為理想信號;ε(i)為噪聲信號;N為信號長度。

則x(i)構造Hankel矩陣A可表示為

A=As+Aε

(6)

式中:As和Aε分別為s(i)和ε(i)所構造的Hankel矩陣，As,Aε∈Rm×n。

無噪聲理想信號所構造的Hankel矩陣As下一行矢量比上一行矢量僅僅滯后1個數據點，因此其相鄰兩行將高度相關，是一種病態的矩陣。這種病態矩陣前面k個奇異值比較大，而之后的奇異值則非常接近于零，奇異值曲線在第k點有明顯的突變，而k實際上就是此矩陣的秩。因此無噪聲理想信號構造的矩陣As的奇異值矢量可近似表示為：λ(As)=(λs1，λs2,…,λsk,0,…,0),其長度為q=min(m,n)，而信號則由前面k個奇異值所代表。

噪聲信號構造的Hankel矩陣Aε盡管相鄰兩個行矢量同樣只滯后一位，但是它們卻并不相關，矩陣是一個良態的滿秩矩陣。而它的所有的奇異值大小均勻，都不為零。因此理想噪聲構造的矩陣Aε的奇異值矢量可表示為：λ(Aε)=(λε，λε,…,λε),其長度為q=min(m,n)。

對于含噪聲的混合信號構造的Hankel矩陣A，其奇異值矢量并不是各分量的奇異值矢量的簡單相加。關于矩陣和的奇異值，理論上只有如下的不等式結論[17]。設矩陣A，B∈Rm×n,而且q=min(m,n)，對于矩陣之和的奇異值，有如下關系

λh(B)] 1≤t≤q

(7)

式中:λh(A)表示矩陣A的奇異值矢量中的第h個奇異值，其余符號類似，而t表示奇異值累加數。其中，當q大到一定程度，矩陣和的奇異值矢量將會有如下近似關系[18]

λ(A)≈(λs1+λε,λs2+λε,…,λsk+λε,λε,…,λε)

(8)

顯然，對于含噪聲的混合信號構造的Hankel矩陣，其后面的q-k個奇異值明顯小于前k個，奇異值在第k點發生突變，而前k個奇異值代表了須提取的理想信號。因此只要選擇前面k個奇異值進行重構，就可以獲得降低了噪聲的信號。

在文獻[14]中提出了一種奇異值差分譜理論，奇異值的差分譜定義為

bk=λk-λk+1k=1,2,…,q-1

(9)

則所有bk形成的序列B=(b1,b2,bq-1)稱為奇異值的差分譜，它描述了兩兩相鄰奇異值的變化情況，當兩相鄰奇異值差別較大時，在差分譜中必將產生一個峰值。

最大差分奇異值分解(Maximized Differential Singular Value Decomposition，MDSVD)就是一種利用差分譜峰值的坐標位置來選擇奇異值的方法：當差分譜最大峰值的坐標位置為1時，說明原始信號中存在一個較大的直流信號，此時以第二最大峰值的坐標位置來確定有效奇異值的數目，否則就以最大峰值的坐標位置來確定有效奇異值數目。

2 自適應奇異值分解故障診斷方法

MDSVD雖然能夠自動選擇有效奇異值數目，但是如果僅選取奇異值差分譜最大峰值點對應奇異值數目來重構的信號會丟失一些弱沖擊信號特征。而根據奇異值差分譜峰值群法又不能自動地選擇有效奇異值數目。因此，為更好地自動保留信號的沖擊特征成分，最大限度地抑制噪聲成分，本文在傳統奇異值分解方法的基礎上，提出了一種能夠自動選擇有效奇異值數目的ASVD方法。

由前述可知，有效奇異值數目在奇異值差分譜上前后差分相差較大的位置。因此，可以設定一個合適閾值d并將滿足式bk/bk+1>d的所有的有效奇異值數目k值求出，得到有效奇異值數目序列K=(k1,k2,…,kp)。這個閥值表示前后差分相差的程度，若閾值太大，則可能丟失一些弱沖擊信號特征；弱閾值太小，則達不到將噪效果。因此，閾值需要針對信號來選取。本文通過不同的d值對比實驗發現，當d=100時特征提取效果最好。然后，根據有效奇異值數目序列K就可以重構出p個重構信號。其中，偏態絕對值計算公式如下所示

(10)

圖1 偏態絕對值對比圖Fig.1 The comparison chart of the absolute values of skewness

于是，本文提出的基于ASVD的行星齒輪箱故障診斷方法步驟如下：

步驟1對含噪信號x(i)構造的Hankel矩陣A奇異值分解，得到其全部奇異值λi；

步驟2將這些相鄰奇異值兩兩相減，得到其差分譜序列B=(b1,b2,…,bq-1)；

步驟3得到有效奇異值數目序列K=(k1,k2,…,kp)；

步驟4分別根據這p個有效奇異值數目得到p個重構信號；

步驟5分別計算出這p個重構信號的偏態絕對值，將偏態絕對值最大的重構信號作為ASVD的結果；

步驟6將分解結果進行Hilbert變換，再將得到的包絡幅值進行Fourier變換，從而得到其包絡譜；

步驟7通過包絡譜故障特征頻率譜線，識別行星齒輪箱的故障類型。

綜上所述，本文處理行星齒輪箱振動信號的具體流程如圖2所示。

3 仿真信號分析

根據行星齒輪箱太陽輪故障振動模型[19]，其仿真信號可以表示為

cos[2πfmt+Bsin(2πfst+φ)+θ]+n(t)

(11)

圖2 故障診斷流程圖Fig.2 The flow chart of fault diagnosis

φ=φ=θ=0，A=B=0，信號x(t)的信噪比為-20 dB，采樣頻率為640 Hz，采樣時間為5 s。

(a) 時域波形

(b) 包絡譜圖3 仿真信號的時域波形及包絡譜

Fig.3 Simulated signal in time domain waveform and envelope spectrum

(a) 時域波形

(b) 包絡譜圖4 MDSVD之后的仿真信號的時域 波形及包絡譜

Fig.4 Simulated signal in time domain waveform and envelope spectrum after MDSVD

(a) 時域波形

(b) 包絡譜圖5 ASVD之后的仿真信號時域 波形及包絡譜

Fig.5 Simulated signal in time domain waveform and envelope spectrum after ASVD

4 實測行星齒輪箱振動信號分析

本文通過傳動系統診斷綜合實驗臺(DDS)模擬了齒輪缺齒故障。實驗臺包含一個2級行星齒輪箱和一個2級平行軸齒輪箱，行星齒輪箱齒輪的分布情況和特征頻率分別如表1和表2所示。其中，加速度傳感器位于行星齒輪箱頂部，輸入轉速為40 r/s，負載為25.2 Nm，采樣頻率為3 200 Hz，故障齒輪為第二級傳動的太陽輪。

表1 行星齒輪箱齒輪的分布情況Tab.1 Distributed situation of gears in planetary gear box

實際行星齒輪箱太陽輪故障振動信號的時域波形及包絡譜如圖6所示。雖然圖6(b)的包絡譜中存在太陽輪故障特征頻率和絕對旋轉頻率的組合頻率，但是由于噪聲的干擾，這些特征譜線并不太突出。因此，需要對該振動信號進行故障特征提取。

表2 行星齒輪箱齒輪的特征頻率Tab.2 Characteristic frequencies of gears in planetary gear box

為了進行對比，首先將振動信號進行MDSVD，所得結果如圖7所示。在圖7(b)的包絡譜中，可以明顯看表2中所計算的所有特征頻率處的幅值幾乎都為零。顯然該方法把故障特征相關成分同噪聲一起消除了。

然后，將振動信號進行ASVD，所得結果如圖8所示。顯然，在圖8(b)的包絡譜中，第二級太陽輪故障特征頻率和絕對旋轉頻率的2倍頻處的幅值十分突出。因此，可以診斷出第二級太陽輪存在局部故障。比較圖7和圖8可知，本文提出的ASVD能準確識別出相應故障，而MDSVD不能。

(a) 時域波形

(b) 包絡譜圖6 振動信號時域波形及包絡譜

Fig.6 Vibration signal in time domain waveform and envelope spectrum

(a) 時域波形

(b) 包絡譜圖7 MDSVD之后的振動信號時域 波形及包絡譜

Fig.7 Vibration signal in time domain waveform and envelope spectrum after MDSVD

(a) 時域波形

(b) 包絡譜圖8 ASVD之后的振動信號時域 波形及包絡譜

Fig.8 Vibration signal in time domain waveform and envelope spectrum after ASVD

5 結 論

在MDSVD的基礎上，本文提出的ASVD根據前后差分比選擇出有效奇異值數目序列K=(k1,k2,…,kp)并根據偏態絕對值最大原則在這p個重構信號中選擇選擇出最佳的重構信號，使得奇異值分解更加科學合理。仿真信號和實際行星齒輪箱振動信號分析結果表明：ASVD能夠很好地消除噪聲，進而有效地提取出行星齒輪箱振動信號中的微弱故障特征，而傳統的MDSVD方法不能正確提取到故障特征。因此，本文提出基于ASVD的故障特征提取方法能夠有效地用于行星齒輪箱中微弱故障的檢測。此外，該方法還可推廣應用于其他旋轉機械振動信號(如滾動軸承)的故障診斷。