潘高元， 李舜酩， 杜華蓉， 朱彥祺

(南京航空航天大學 能源與動力學院，南京 210016)

汽車的齒輪箱是傳動系統(tǒng)中的重要組成部分，而齒輪箱中有60%的失效形式與齒輪有關[1]。齒輪常見的故障有齒面磨損、齒根裂紋、斷齒等，其中斷齒是最嚴重的損傷形式，直接影響機構的正常運轉。齒輪元件在工作中發(fā)生局部故障，相應的特征將會表現在振動信號中。因此，通過對齒輪箱振動信號的采集和分析，找出故障源及其可能產生的原因，對確保齒輪箱及時、準確的維修和調試，具有重要的意義。在理想情況下，齒輪箱斷齒故障信號中的沖擊特征出現頻率即為斷齒旋轉的頻率。但是齒輪的故障信號往往是周期嚙合振動、調幅、調頻及附加脈沖的綜合表現[2]，是典型的非平穩(wěn)信號。這導致了故障源振動信號的沖擊特征通常淹沒在背景信號與噪聲中，比較難以直接識別。因此，如何從帶有噪聲背景的信號中準確提取故障特征，進而進行齒輪箱的故障診斷是目前研究的熱點問題。

S變換是最近發(fā)展起來的一種時頻分析方法，在地球物理學領域已經證明了S變換的有效性，如分析內在的大氣波包、大氣研究、地震信號的鑒定、全球海洋表面溫度分析、電氣工程、機械工程、數字信號處理等。它和連續(xù)小波變換相同的是漸進分辨率，但和小波變換不同的是，S變換保留了絕對參考相位信息，且S變換不但可以估計局部的功率譜，還可以估計局部的相位譜，這也適用一般的復值時間序列。并且S變換克服了短時Fourier變換窗函數以及連續(xù)小波變換基函數固定不變的缺點，是一種高效的自適應信號時頻分析方法，適合處理與分析非平穩(wěn)信號，尤其是包含沖擊特征的信號，其逆變換完全無損[3]。

奇異值分解(Singular Value Decomposition,SVD)降噪方法是一種非線性濾波方法，可以有效抑制信號中的寬帶隨機噪聲。利用選取奇異值閾值的方式，可以減少信號中的噪聲，提高信號識別的準確性[4]。而利用SVD對一維時域信號進行處理與分析，關鍵問題之一是構造合適的數據矩陣[5-6]。顯然，需要得到能夠明顯識別沖擊信號的矩陣，而S變換作為對沖擊信號敏感的時頻變換方法，在提取齒輪箱斷齒特征方面非常合適于SVD的構造矩陣。

因此，基于S變換時頻譜SVD降噪的沖擊特征提取方法，能夠很好的提取出信號中的沖擊特征。而奇異值閾值的選擇，則是此方法的關鍵之一。朱怡等[7-8]利用奇異值差分譜，選取奇異值差分譜最前面峰值群最后一個峰值的方式確定閾值位置，有效地提取軸承故障振動信號中的沖擊特征頻率。但是利用差分譜選取閾值的方法在一定條件下具有試探性，有時不能直觀的挑選出閾值位置。此外，利用時域特征的進行沖擊特征的識別，對軸承的某些故障特征頻率識別已經足夠，而行星齒輪箱行星輪斷齒沖擊特征比軸承的沖擊特征更復雜，單獨分析時域信息會導致特征識別不全。因此本文提出使用奇異值比值譜的方法，結合奇異值差分譜與奇異值比值譜選取閾值，能夠在識別出沖擊特征的情況下，直觀的選擇出奇異值閾值。并且在識別齒輪斷齒沖擊特征時，同時考慮時域特征和S變換的時頻譜圖，從而能夠正確地識別出齒輪箱斷齒故障。

1 S變換時頻譜SVD降噪理論

1.1 齒輪斷齒時域頻域特征

齒輪斷齒時的時域表現為幅值很大的沖擊型振動，沖擊頻率理論上等于存在斷齒軸的轉頻。在頻域上表現為在嚙合頻率及其高次諧波附近出現間隔為斷齒軸轉頻的邊頻帶；邊頻帶一般數量多、幅值較大、分布較寬[9]。

齒輪斷齒的主要特征為：①以齒輪嚙合頻率及其高次諧波為載波頻率，齒輪所在軸轉頻及其倍頻為調制頻率的嚙合頻率調制，調制邊頻帶寬而高；②以齒輪各階固有頻率為載波頻率，齒輪所在軸轉頻及其倍頻為調制頻率的齒輪共振頻率調制，調制邊頻帶寬而高；③振動能量(包括有效值、峰值等)有較大幅度的增加；④時域圖中沖擊頻率等于有斷齒軸的轉頻。

1.2 S變換

S變換是地球物理學家Stockwell提出的一種將信號從一維時域信號變換到二維時頻的一種信號處理方法，其具有可逆性。其基本思想是對短時F變換和以Morlet小波為基本小波的連續(xù)小波變換的組合和延伸?？梢詫變換認為是一種加了高斯窗函數且窗寬度與信號頻率成反比的特殊的短時傅里葉變換STFT(Short-Time Fourier Transform)，因此在低頻段具有較高的頻率分辨率，在高頻段具有較高的時間分辨率，并且對非平穩(wěn)信號例如沖擊信號十分敏感，滿足線性疊加原理，不存在交叉項的干擾，是一種無損的且可逆的時頻變換。

信號x(t)的S變換為

(1)

式中：f為頻率；t為時間；b為時間軸上的位移參數。

1.3 SVD降噪理論

設矩陣A為m×n的矩陣，矩陣的秩為r，則必存在m×n的正交矩陣U和正交矩陣V，使得

A=USVT

(2)

(3)

即除了前r階對角線元素外，矩陣S其他元素都為0，并且稱σi(i=1,2,…，r)為矩陣A的奇異值。奇異值跟特征值類似，在矩陣中也是從大到小排列，奇異值大小其主要反映各元素的能量集中情況。而且奇異值減少特別的快，在很多情況下，前10%甚至1%的奇異值的和就占了全部的奇異值之和的99%以上了。一般使用Hankel矩陣作為SVD的構造矩陣。由于沖擊特征的幅值通常都呈指數衰減形式，此情況下Hankel矩陣不適合作為構造矩陣。而沖擊信號在S變換的時頻譜中比較集中，噪聲信號則比較分散，因此S變換更適合作為SVD的構造矩陣。這樣，奇異值中較大奇異值就主要反映信號的沖擊成分，而較小的奇異值則反映信號的噪聲與小的沖擊部分。因此選擇一定的奇異值設置閾值，將值較小的奇異值置零就可以消除部分噪聲，使得沖擊成分凸顯出來。

2 奇異值閾值降噪方法與步驟

2.1 奇異值比值譜閾值降噪方法

奇異值閾值的選取決定了SVD降噪的效果，閾值數值過大會丟失重要的沖擊信號，閾值數值過小則會使得降噪效果不明顯。因此選擇合適的閾值，在噪聲降低的情況下不丟失沖擊特征，是準確提取出沖擊特征的關鍵之一。一種方法是利用奇異值差分譜的奇異值閾值確定方法來尋找合適的奇異值閾值[10]。定義奇異值差分譜為

P=(σ1-σ2,σ2-σ3,σ3-σ4,…,σr-1-σi)

(4)

由于奇異值按照降序排列，因此差分譜全部為正數。若僅選擇奇異值差分譜最大峰值點對應的奇異值進行信號重構，會丟失一些弱沖擊信號特征。為更好地保留信號的沖擊特征成分，最大限度地抑制噪聲成分，將奇異值差分譜中峰值群中最后一個峰值點序號對應的奇異值作為奇異值閾值[11]，能夠很好地確定閾值的位置。所謂的峰值群，即為奇異值差分譜最前面部分較為集中的一組峰值點，并且其幅值顯著大于后續(xù)峰值點的幅值。

但是這種方法確定閾值時在某些情況下具有試探性[12]，且有一定的人為選擇的差異性，即在某些情況下峰值群的幅值無法顯著大于后續(xù)峰值點的峰值，不能夠直觀的選取合適的位置。

結合Kanjilal等[13-14]僅使用前兩個奇異值σ1/σ2的奇異值比譜(Singular Value Ratio,SVR)的方法，為了考慮全部奇異值，本文提出奇異值比值譜的方法對奇異值閾值進行選擇。定義奇異值比值譜為

Q=(σ1/σ2,σ2/σ3,σ3/σ4,…,σi-1/σ1)

(5)

由于奇異值按照降序排列，因此比值譜全部大于一。同樣的，選取奇異值比值譜中最前面峰值群的最后一個峰值點序號作為奇異值閾值。由于奇異值代表的是信號的能量集中情況，相鄰的奇異值之間比值能夠反映出變化趨勢，在識別沖擊特征上，和差分譜方法相比有異曲同工之妙。但是由于沖擊特征的幅值通常都呈指數衰減形式，而且奇異值的前幾個數都特別大，而后面的數值比較小，而S變換時頻譜沖擊特征更為集中，因此對于衰減較快的沖擊特征，奇異值比值變化相比較差值變化更為明顯，更容易選取出奇異值閾值。

此外，由于使用SVD降噪的方法對信號進行重構，使用信號的時域信息進行沖擊信號分析，對沖擊周期的倒數作為沖擊信號的頻率特征，因此實驗數據的時間分辨率需要有足夠的精度才能夠準確獲取沖擊特征。如果信號的時間分辨率較差，當沖擊信號時間間隔比較小時，一個時間分辨率誤差就能引起較大的頻率誤差。因此若分析完信號發(fā)現特征的頻率比較高，則需要考慮是否對信號重采樣，重新計算以獲取更準確的頻率。若頻率比較低，則對故障特征識別已經有足夠的精度。

2.2 信號沖擊特征提取的步驟

按照上述思路，為了從包含噪聲的沖擊信號x(t)中提取出沖擊特征，對其進行S變換時頻譜SVD降噪處理，其主要步驟如下：

步驟1對信號x(t)進行S變換，得到S變換時頻譜系數矩陣A。

步驟2對S變換得到的時頻譜矩陣A進行奇異值分解，得到矩陣A的奇異值，并將奇異值按照遞減的順序排列，即σ1≥σ2≥…≥σr。

步驟3求解其奇異值差分譜和比值譜并比較，選取奇異值比值譜中最前面峰值群的最后一個峰值點序號作為奇異值閾值，并設置為閾值σh，將奇異值序列中σh后的奇異值置為零，即對信號x(t)的S變換時頻譜進行SVD降噪處理，然后重建S變換時頻譜系數矩陣。

步驟4對重建后的時頻譜矩陣進行S逆變換，得到x(t)降噪后的時域沖擊特征，根據沖擊特征的時間差值得到沖擊的周期，從而得到沖擊信號的特征頻率。

步驟5結合判斷S變換后的二維時頻矩陣，判斷x(t)時域特征的信息是否正確完整，并根據采樣的時間分辨率以及得到沖擊特征頻率值大小判斷是否需要重新采樣。

3 仿真信號分析

3.1 較高信噪比信號仿真分析

仿真設計的沖擊信號x(t)的沖擊成分由以下的指數衰減形成

P(t)=e-ξωnπt×cos(2πft)

(6)

圖1 沖擊信號p(t)Fig.1 Impact signal p(t)

首先，對信號x(t)進行S變換，得到其時頻矩陣，其云圖如圖2所示。從圖2中可以看出，在相對頻率為8 Hz左右出現十分明顯的沖擊信號，說明S變換識別沖擊特征的效果還是十分明顯的。

圖2 信號x(t)的S變換時頻譜圖Fig.2 S-transform time-frequency spectrum of the signal x(t)

其次，對信號x(t)的S變換時頻矩陣進行SVD奇異值分解，求得其奇異值，并將奇異值按照降序排列。由于奇異值主要信息存儲在前面較大的奇異值之中，因此只顯示其前一百個奇異值，如圖3所示。從圖3可以看出在第75個奇異值之后趨勢變化已經非常緩慢，基本接近于零。

然后，對奇異值進行差分譜和比值譜的求解，分別按照式(4)、式(5)所定義的公式進行計算，得到的奇異值差分譜和比值譜的曲線如圖4和圖5所示。

圖3 S變換時頻譜矩陣的奇異值譜Fig.3 Singular value spectrum of S-transform spectrum time-frequency matrix

圖4 S變換時頻矩陣的奇異值譜的差分譜Fig.4 Singular value difference spectrum of S-transform spectrum time-frequency matrix

圖5 S變換時頻矩陣的奇異值譜的比值譜Fig.5 Singular value ratio spectrum of S-transform spectrum time-frequency matrix

最后，對信號進行降噪，重構時域信號。由于信號信噪比高，因此奇異值應該保留的多而置零的少。可以看出，奇異值最前面的峰值群上的最后一個峰值點，根據差分譜的結果，由于峰值一直衰減，在選取閾值時選第26個還是第31個還是第46個，存在人為選擇的差異，具有不確定性，很難直接按照峰值群進行判斷。而奇異值比值譜則明顯能看出，在第46個奇異值位置，比值譜的峰值就已經都能夠明顯的與后面的峰值分辨開。

因此，根據奇異值差分譜的結果，本次選則將第31個奇異值之后的奇異值置零，重建S變換時頻矩陣，然后進行S逆變換得到降噪后的信號時域圖如圖6(a)所示。根據比值譜的結果，將第46個奇異值之后的奇異值置零，重建S變換時頻矩陣，然后進行S逆變換得到降噪后的信號時域圖如圖6(b)所示。可以看出根據奇異值比值譜的方法選取閾值。盡管保留的奇異值比差分譜方法多，但是已經能夠提取出信號的沖擊特征，噪聲在時域特征上已經很少，相比較差分譜選擇閾值的方式，比值譜方法在提取沖擊特征的同時保留了更多的沖擊特征，并且比值譜方法更容易判斷出閾值的位置。

圖6 x(t)選取不同閾值降噪后的沖擊信號Fig.6 The denoised impact signal using different thresholds

3.2 低信噪比信號仿真分析

仿真沖擊信號x(t)的沖擊成分p(t)與不含噪聲的沖擊信號相同，在信號中添加的幅值序列均值為1，標準差為0.3的高斯隨機序列噪聲。含噪聲的沖擊信號x(t)的時域圖如圖7所示。

圖7 含高斯噪聲的x(t)沖擊信號Fig.7 Impact signal x(t) with Gaussian noise

同樣按照前面所述的信號分析步驟，首先對信號x(t)其進行S變換，得到的S變換時頻圖如圖8所示，可以看出在8 Hz左右仍然能看到沖擊特征，但是沖擊特征很不明顯，說明沖擊信號淹沒在噪聲中。其次對S變換的時頻矩陣進行SVD奇異值分解，并將奇異值按照降序排列，選取其中前100個奇異值如圖9所示。然后同樣的計算出奇異值差分譜和比值譜曲線，如圖10和圖11所示。

圖8 信號x(t)的S變換時頻譜圖Fig.8 S-transform time-frequency spectrum of the signal x(t)

圖9 S變換時頻譜矩陣的奇異值譜Fig.9 Singular value spectrum of S-transform spectrum time-frequency matrix

圖10 S變換時頻矩陣的奇異值譜的差分譜Fig.10 Singular value difference spectrum of S-transform spectrum time-frequency matrix

圖11 S變換時頻矩陣的奇異值譜的比值譜Fig.11 Singular value ratio spectrum of S-transform spectrum time-frequency matrix

從圖10和圖11可以看出，差分譜和比值譜中最前面部分出現峰值群，都以峰值群中最后一個峰值點的序號16作為奇異值閾值的位置坐標，即選擇閾值σh=16，對后面的奇異值進行置零。但是比值譜曲線的最后一個峰值變化顯然比差分譜曲線更明顯。然后重建S變換降噪后的時頻矩陣，對其進行S逆變換，從而得到SVD降噪處理后的時域信號，如圖12所示。

圖12 x(t)降噪后沖擊信號Fig.12 The denoised impact signal of x(t)

由圖12結果可以看出，利用奇異值差分譜和比值譜，都可以保留源信號的主要沖擊成分，將其從噪聲中提取出來，而且兩者選擇閾值是一致的。而奇異值比值譜的方式，第16個位置處的幅值比后面差別更大，選取閾值時候更直觀。

綜上，當信噪比較低時，比值譜和差分譜都能夠篩選出接近一致的閾值，提取出信號中的沖擊信號；而當信噪比較高時，比值譜能夠保留更多的沖擊信息。相比較之下，使用奇異值比值譜選取閾值的方式，其峰值群的變化更大，更容易直接選取合適的閾值。

4 齒輪箱斷齒振動信號分析

4.1 實驗設備與基本參數

為驗證奇異值比值譜選取閾值的方法在齒輪箱斷齒故障信號中提取沖擊特征的有效性與實用性，下面采用實驗進行方法的驗證。

實驗采用電機-齒輪箱-負載轉盤的形式進行連接，斷開轉盤與發(fā)動機的連接，部件之間用聯軸器連接，設備的連接圖如圖13所示。電機型號為WF2-180L-4，運行工況電機恒定轉速。齒輪箱采用降速的形式進行連接，即輸入軸太陽輪轉速等于電機轉速n=1 500 r/min。傳感器采用振動加速度傳感器，連接如圖14所示。齒輪箱為一級行星齒輪箱，為3個行星輪，1個太陽輪，1個固定內齒圈，行星架為輸出，主要參數參見表1。齒輪故障類型為其中一個行星輪斷一個齒，如圖15所示。

圖13 齒輪箱設備連接圖Fig.13 Gearbox equipment connection

圖14 齒輪箱加速度傳感器位置Fig.14 Position of gearbox acceleration sensor

圖15 齒輪箱行星輪故障齒輪Fig.15 Fault gear of gearbox planetary gear

表1 測試齒輪箱參數Tab.1 Parameters of the test gearbox

為了獲得完整的沖擊信號頻率特征，需要獲取至少兩個周期的時間序列；為了提高計算效率，則時間序列長度不能過大。實驗設置采樣頻率為fs=12.8 kHz，時間分辨率為t=7.812 5×10-5s，根據電機的轉速可知，截取采樣點數N=4 096就已經達到了分析的要求。

4.2 實驗信號分析

選取輸入端太陽輪處的加速度傳感器的進行分析，加速度信號沿齒輪軸的徑向方向。實驗所測得的故障齒輪箱的振動加速度時域信號y(t),如圖16所示。

圖16 故障齒輪的原振動加速度信號y(t)Fig.16 The original vibration acceleration signal y(t) of the fault gear

對時域信號進行自相關處理，然后進行FFT(Fast Fourier Transformation)計算得到頻譜圖，可以得到主要信號的基頻為fn=25 Hz。由實驗參數可知太陽輪基頻為1 500/60=25 Hz，與實驗參數電機轉速相符。根據齒輪的傳動比以及齒輪齒數可知，其太陽輪基頻為25 Hz，輸出端行星架基頻為8.33 Hz。

首先對所得信號進行FFT，得到信號的自功率譜圖，如圖17所示。從自功率譜圖可以看出，齒輪的整個信號主要集中在基頻25 Hz以及其倍頻段和500～850 Hz的頻段。因此取0～1 000 Hz的自功率譜圖的局部放大圖,如圖18所示。

圖17 信號y(t)的自功率圖Fig.17 Autopower of the signal y(t)

圖18 信號y(t)自功率圖的局部放大圖Fig.18 Partial enlargement of the autopower diagram of the signal y(t)

在500 Hz以下的頻率，除了太陽輪自轉的基頻25 Hz外，主要由50 Hz，75 Hz，200 Hz，300 Hz，400 Hz等基頻的倍頻組成，加速度振動幅值較小。而在600.13 Hz等處能量比較集中，在583.50 Hz，591.75 Hz，608.38 Hz，616.75 Hz，625.13 Hz，633.50 Hz，641.75 Hz，650.13 Hz，691.75 Hz，716.75 Hz等頻率，相鄰的幅值的差值基本上在8.25～8.38 Hz，在誤差范圍內基本為輸出軸的轉動的基頻。使用FFT主要獲取以上頻率信息。

按照S變換SVD降噪的方法，首先對時域數據進行S變換，得到的S變換的時頻圖如圖19所示。可以看出，在信號的相對頻率700 Hz左右出現比較明顯的沖擊特征，在相對頻率1 500 Hz以上基本沒有沖擊信號。

圖19 信號y(t)的S變換時頻譜圖Fig.19 S-transform time-frequency spectrum of the signal y(t)

然后對S變換后的時頻矩陣進行SVD奇異值分解。根據奇異值的特點，選取其前一百個奇異值進行分析，其差分譜和比值譜，分別如圖20和圖21所示。

由奇異值差分譜與比值譜可以判斷，兩者總趨勢是一致的，但是每個圖在不同的局部(例如第10個峰值)卻是不一樣的。綜合考慮兩個曲線，選擇第23個奇異值處作為奇異值置零的閾值，對此奇異值后面的奇異值進行歸零處理，重建S變換時頻矩陣，然后對其進行SVD逆運算，得到降噪后的S變換的時頻矩陣。最后對時頻矩陣進行S變換逆變換，得到降噪后的時域的沖擊信號如圖22所示。

圖20 S變換時頻矩陣的奇異值譜的差分譜Fig.20 Singular value difference spectrum of S-transform spectrum time-frequency matrix

圖21 S變換時頻矩陣的奇異值譜的比值譜Fig.21 Singular value ratio spectrum of S-transform spectrum time-frequency matrix

圖22 y(t)降噪后時域信號Fig.22 The denoised impact signal of y(t)

從降噪后的時域圖可以得到沖擊信號的幾個時刻點，分別在0.064 3 s，0.076 95 s，0.108 7 s，0.184 2 s，0.228 4 s，0.304 3 s處，相鄰的兩個時間差值的倒數從13.17～98.5 Hz不等，與齒輪斷齒的信號頻率完全不符。究其原因，與滾動軸承內圈、外圈故障不同，對于行星齒輪箱而言，行星齒輪和保持架內圈齒以及太陽輪齒輪同時嚙合，因此單個的行星齒輪斷齒會在不同的時間沖擊太陽輪和保持架內齒圈，因此兩個信號的周期是相互重疊的，對太陽輪的兩個沖擊信號時域間會帶有對保持架內齒圈的沖擊信號，對保持架同理。因此，單純的判斷降噪后的時域信息，對行星齒輪而言是不能得到正確的沖擊信號。因此結合S變換的時頻矩陣進行分析，舍去基本無沖擊信號的頻率部分，得到S變換的二維時頻云圖的局部放大圖，如圖23所示。

圖23 信號y(t)的S變換時頻譜圖局部放大圖Fig.23 Partial enlargement of S-transform time-frequency spectrum of the signal y(t)

結合齒輪斷齒的時頻特征，可以判斷出行星齒輪存在斷齒故障。并且使用S變換時頻矩陣SVD降噪，利用奇異值比值譜進行奇異值閾值篩選的方法進行齒輪斷齒特征的識別是可行的。

5 結 論

(1)S變換是一種信號時頻表示方法，對于信號中的沖擊成分具有較高敏感性，可以很好地表征信號的沖擊特征。

(2)選取奇異值比值譜最前面部分峰值群的最后一個峰值點序號作為奇異值序列置零閾值的方法，具有良好的降噪效果，并且比利用差分譜在閾值選取時峰值群特征更明顯一些。

(3)將S變換與SVD降噪相結合的方法用于齒輪箱的振動信號處理中，成功的識別出了頻率在8.33 Hz附近的沖擊振動特征，與行星齒輪的斷齒特征一致，準確的識別出了齒輪的斷齒特征。

(4)使用S變換與SVD降噪相結合的特征提取方法，需要結合判斷S變換時頻圖中的沖擊特征，才能在時域圖中準確的識別沖擊信號的周期，最終得到正確的沖擊特征頻率。