秦 毅， 張清亮， 趙 月

(1. 重慶大學(xué) 機(jī)械傳動(dòng)國家重點(diǎn)實(shí)驗(yàn)室， 重慶 400044； 2. 重慶大學(xué) 機(jī)械工程學(xué)院， 重慶 400044)

行星齒輪箱具有結(jié)構(gòu)緊湊、體積小、質(zhì)量輕、承載能力強(qiáng)、傳動(dòng)比大、傳動(dòng)效率高等優(yōu)點(diǎn)[1]，已被廣泛應(yīng)用于諸多機(jī)械傳動(dòng)系統(tǒng)。行星齒輪箱大部分使用在低速重載的惡劣工況，太陽輪、行星輪、內(nèi)齒圈和行星架等部件很容易出現(xiàn)故障。對(duì)行星齒輪箱進(jìn)行故障診斷可避免因事故造成的巨大經(jīng)濟(jì)損失，具有重大的意義。

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

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

1 奇異值分解理論

1.1 奇異值分解

設(shè)離散數(shù)字信號(hào)X=[x(1),x(2),…,x(N)]，則利用此信號(hào)可以構(gòu)造Hankel矩陣如下

(1)

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

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

A=U∑VT

(2)

重構(gòu)信號(hào)就是要選擇Hankel矩陣A的奇異值之中的前面k個(gè)有效奇異值并將之后r-k的個(gè)奇異值置零，并進(jìn)行奇異值分解逆變換，即

A*=U∑*VT

(3)

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

x*(i)=

(4)

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

含噪聲的離散數(shù)字信號(hào)x(i)可表示為

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

(5)

式中:s(i)為理想信號(hào);ε(i)為噪聲信號(hào);N為信號(hào)長(zhǎng)度。

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

A=As+Aε

(6)

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

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

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

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

λh(B)] 1≤t≤q

(7)

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

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

(8)

顯然，對(duì)于含噪聲的混合信號(hào)構(gòu)造的Hankel矩陣，其后面的q-k個(gè)奇異值明顯小于前k個(gè)，奇異值在第k點(diǎn)發(fā)生突變，而前k個(gè)奇異值代表了須提取的理想信號(hào)。因此只要選擇前面k個(gè)奇異值進(jìn)行重構(gòu)，就可以獲得降低了噪聲的信號(hào)。

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

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

(9)

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

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

2 自適應(yīng)奇異值分解故障診斷方法

MDSVD雖然能夠自動(dòng)選擇有效奇異值數(shù)目，但是如果僅選取奇異值差分譜最大峰值點(diǎn)對(duì)應(yīng)奇異值數(shù)目來重構(gòu)的信號(hào)會(huì)丟失一些弱沖擊信號(hào)特征。而根據(jù)奇異值差分譜峰值群法又不能自動(dòng)地選擇有效奇異值數(shù)目。因此，為更好地自動(dòng)保留信號(hào)的沖擊特征成分，最大限度地抑制噪聲成分，本文在傳統(tǒng)奇異值分解方法的基礎(chǔ)上，提出了一種能夠自動(dòng)選擇有效奇異值數(shù)目的ASVD方法。

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

(10)

圖1 偏態(tài)絕對(duì)值對(duì)比圖Fig.1 The comparison chart of the absolute values of skewness

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

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

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

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

步驟4分別根據(jù)這p個(gè)有效奇異值數(shù)目得到p個(gè)重構(gòu)信號(hào)；

步驟5分別計(jì)算出這p個(gè)重構(gòu)信號(hào)的偏態(tài)絕對(duì)值，將偏態(tài)絕對(duì)值最大的重構(gòu)信號(hào)作為ASVD的結(jié)果；

步驟6將分解結(jié)果進(jìn)行Hilbert變換，再將得到的包絡(luò)幅值進(jìn)行Fourier變換，從而得到其包絡(luò)譜；

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

綜上所述，本文處理行星齒輪箱振動(dòng)信號(hào)的具體流程如圖2所示。

3 仿真信號(hào)分析

根據(jù)行星齒輪箱太陽輪故障振動(dòng)模型[19]，其仿真信號(hào)可以表示為

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

(11)

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

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

(a) 時(shí)域波形

(b) 包絡(luò)譜圖3 仿真信號(hào)的時(shí)域波形及包絡(luò)譜

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

(a) 時(shí)域波形

(b) 包絡(luò)譜圖4 MDSVD之后的仿真信號(hào)的時(shí)域 波形及包絡(luò)譜

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

(a) 時(shí)域波形

(b) 包絡(luò)譜圖5 ASVD之后的仿真信號(hào)時(shí)域 波形及包絡(luò)譜

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

4 實(shí)測(cè)行星齒輪箱振動(dòng)信號(hào)分析

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

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

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

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

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

然后，將振動(dòng)信號(hào)進(jìn)行ASVD，所得結(jié)果如圖8所示。顯然，在圖8(b)的包絡(luò)譜中，第二級(jí)太陽輪故障特征頻率和絕對(duì)旋轉(zhuǎn)頻率的2倍頻處的幅值十分突出。因此，可以診斷出第二級(jí)太陽輪存在局部故障。比較圖7和圖8可知，本文提出的ASVD能準(zhǔn)確識(shí)別出相應(yīng)故障，而MDSVD不能。

(a) 時(shí)域波形

(b) 包絡(luò)譜圖6 振動(dòng)信號(hào)時(shí)域波形及包絡(luò)譜

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

(a) 時(shí)域波形

(b) 包絡(luò)譜圖7 MDSVD之后的振動(dòng)信號(hào)時(shí)域 波形及包絡(luò)譜

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

(a) 時(shí)域波形

(b) 包絡(luò)譜圖8 ASVD之后的振動(dòng)信號(hào)時(shí)域 波形及包絡(luò)譜

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

5 結(jié) 論

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