郭遠晶，魏燕定，周曉軍，唐 昉

(浙江大學 現(xiàn)代制造工程研究所，杭州 310027)

旋轉機械設備出現(xiàn)某些故障時，如齒輪斷齒、齒面磨損、軸承滾道或者滾動體表面金屬剝落等，很大程度上會出現(xiàn)具有沖擊特征的振動狀態(tài)。但由于旋轉機械設備傳動系統(tǒng)的復雜性以及工作環(huán)境的多樣性，各種激勵源產(chǎn)生的振動信號相互耦合，導致故障源振動信號通常淹沒在強背景信號與噪聲中，比較難以識別。因此，若能成功提取所采集振動信號中的沖擊特征，即可方便有效地對旋轉機械設備相關故障進行診斷。

沖擊特征信號是典型的非線性、非高斯信號。對于沖擊特征的提取，國內(nèi)外學者已經(jīng)進行了一些相關的研究。LEE等［1－2］利用高階統(tǒng)計量，如峭度、偏斜度等，對信號沖擊成分具有高度敏感性的特點，提出以高階統(tǒng)計量為優(yōu)化目標函數(shù)的盲解卷積算法，從而提取沖擊特征。李繼猛等［3］提出了基于粒子群算法的多參數(shù)同步優(yōu)化自適應單穩(wěn)態(tài)隨機共振方法，最優(yōu)地檢測出了原始信號中的沖擊成分。章立軍等［4］提出一種基于多尺度Top-Hat變換的形態(tài)非抽樣小波分解方法，成功提取出了振動信號中的沖擊特征。劉小峰等［5］采用Chirplet作為基函數(shù)對非平穩(wěn)的瞬態(tài)沖擊信號進行自適應分解，然后分別提出了時頻濾波與模型逼近、Chirplet函數(shù)直接重構的沖擊特征提純方法。針對非高斯信號的處理，Simon等［6］首先提出基于高階累積量和緊縮原理的時域迭代盲解卷積算法，王宇等［7］則在該算法的基礎上，提出了盲解卷積和分層聚類的改進算法，并在機械聲信號弱沖擊成分提取中獲得了較好的效果。楊富春等［8］利用信號的零時滯四階累積量即峰態(tài)對非高斯信號非常敏感的特點，提出了滑動峰態(tài)算法的信號沖擊特征提取方法。于明月等［9］利用所提出的雙自適應小波局部極大模方法，有效地提取出了滾動軸承振動信號中的微弱沖擊特征。

從噪聲混合的信號中提純沖擊特征，國內(nèi)外學者針對各自研究的對象信號提出了各種行之有效的方法，但最直接的方法依然是對信號進行降噪處理，以凸顯沖擊特征。因此，本文提出一種基于S變換譜閾值去噪的沖擊特征提取方法。該方法以包含沖擊特征的噪聲混合信號為研究對象，先利用S變換獲得其時頻譜，然后根據(jù)譜系數(shù)的模值大小進行閾值去噪。去噪過程中，閾值函數(shù)分別采用了硬閾值函數(shù)和軟閾值函數(shù)，最優(yōu)閾值則利用步長迭代算法獲取。最后利用S逆變換重構得到時域沖擊特征。

1 S變換

1.1 連續(xù)S變換

S變換由 Stockwell在分析地球物理數(shù)據(jù)時提出［10－11］。對于給定的連續(xù)時間信號x(t)，其S變換定義為

式中，S(τ，f)表示 x(t)的連續(xù) S變換，t表示時間，f表示頻率，參數(shù)τ控制著分析窗ω在t軸上位置。分析窗ω在時域定義為高斯窗

同時將控制窗寬的參數(shù)σ定義為信號頻率的倒數(shù)，即

由此可見，S變換所用高斯窗的時寬與信號頻率成反比關系，這使得S變換在信號低頻段頻率分辨率高，高頻段時間分辨率高，也即S變換具有多分辨率特性，克服了短時Fourier分辨率固定的缺點，且對于頻率突變且又暫態(tài)的信號沖擊特征具有較高的敏感性。此外，S變換還滿足線性疊加原理，不存在交叉項的干擾。假設噪聲混合信號x(t)=s(t)+n(t)，n(t)為噪聲信號，則有

S變換的以上特點使得其非常適合于處理與分析非平穩(wěn)信號，尤其是包含沖擊特征的信號。S變換譜表征的是信號局部頻譜，因此，對信號x(t)的S變換譜沿時間軸積分，可以得到其Fourier譜，即

從而由S(τ，f)完全重構時域信號x(t)，即

1.2 離散S變換

對于時域離散信號 x(k)(k=0，1，2，…，N －1)，其離散S變換的計算公式如下所示

其中，j，n=0，1，…，N － 1，X(n)為 X(k)的傅里葉變換。

離散S逆變換為

2 S變換譜閾值去噪

為了實現(xiàn)對沖擊特征信號x(t)的S變換時頻譜閾值去噪，以提純其沖擊特征，首先需要將x(t)進行S變換，得到其時頻譜 S(j，n)(j，n=0，1，2，…，N －1)。由于S(j，n)為一個N×N的復數(shù)矩陣，對譜S(j，n)進行閾值去噪的方式便有兩種:① 對譜系數(shù)的實部和虛部分別進行去噪;② 根據(jù)譜系數(shù)的模值大小進行去噪?？紤]到前種方式會改變系數(shù)的相位，從而使S逆變換重構得到的時域信號產(chǎn)生附加噪聲，所以本文采用后種方式。

為對S變換譜進行閾值去噪，首先需要獲取最優(yōu)閾值，針對于此，本文提出步長迭代算法。顯然，最優(yōu)閾值在［0，max{|S(j，n)|］］的 S 變換譜系數(shù)模值區(qū)間內(nèi)。為了減少計算量，可以根據(jù)實際情況，定義最優(yōu)閾值估計的一個較小的模值目標區(qū)間［α，β］·max{|S(j，n)|］，其中，［α，β］為模值系數(shù)區(qū)間，且 0≤α ＜ β≤1。最簡單地，設定［α，β］=［0，1］，即 S 變換譜系數(shù)模值的全局區(qū)間;一般地，可以設定如［α，β］=［0.1，0.6］的子區(qū)間。因此，獲取最優(yōu)閾值的步長迭代算法具體步驟如下所述。

步驟1，將模值目標區(qū)間或者系數(shù)區(qū)間分成M個步長，則每個系數(shù)步長Δd=(β－α)/M，模值步長ΔD= Δd·max{|S(j，n)|］;閾值變量 δth初始化為 α·max{|S(j，n)|}。

步驟2，令閾值變量δth按模值步長ΔD增大，每增大一個模值步長，得到一個新閾值δnew。

步驟3，以δnew為閾值(包括初始值α·max{|S(j，n)|})進行S變換譜閾值去噪;去噪過程中所采用的閾值函數(shù)可以為硬閾值函數(shù)、軟閾值函數(shù)或者其它的閾值函數(shù);最后對去噪后的時頻譜進行S逆變換，重構時域信號。

到這里，有必要探討最優(yōu)閾值的評價標準。

假設信號 x(k)(k=0，1，2，…，N －1)中真正的信號為s(k)，并包含有加法性噪聲u(k)，則x(k)表示為x(k)=s(k)+u(k)。再令S變換譜經(jīng)過閾值去噪后的S變換譜為S'(j，n)，那么經(jīng)S逆變換重構得到的時域信號為

最優(yōu)閾值δ*th去噪的目標是使s'(k)與真正信號s(k)的均方誤差(Mean Squared Error)最小，也即使所謂的風險函數(shù)L2－RISK最小化，其定義為

但是，在處理實際信號時，s(k)是不可知的，L2－RISK無法計算，針對于此，本文提出一種改進風險函數(shù)，其定義為

式中，r表示步數(shù)變量，xr(k)為閾值變量δth增加了r個模值步長，即 δth= α·max{|S［j，n］ |］ +r·ΔD 時，經(jīng)過S變換譜閾值去噪后重構得到的時域信號。這樣的話，RM(r)就定義為第r步閾值去噪得到的時域信號xr(k)與前一步閾值去噪得到的時域信號xr－1(k)的均方誤差。

因此，在步驟3中，重構得到時域信號后，隨之計算改進風險函數(shù)RM(r)。

步驟4，重復步驟2、3，直到閾值變量δth增加到β·max{|S［j，n］ |］，此時可以得到改進風險函數(shù) RM(r)關于步數(shù)變量r或者閾值變量δth的關系曲線。

步驟5，對于硬閾值去噪，RM(r)曲線水平趨勢段起始點所對應的閾值變量δth值即可認為是最優(yōu)閾值。然而，此R(r)曲線一般為非平滑不規(guī)則曲線，需M要進一步根據(jù)M個RM(r)離散值進行多項式擬合，以得到RM(r)曲線的水平趨勢。為使最優(yōu)閾值δ*th的獲取更為準確，擬合曲線的水平段范圍應該可能寬，因而擬合次數(shù)也應該盡可能高，但其有限制，次數(shù)過高會使擬合曲線不能明顯呈現(xiàn)甚至沒有水平段。最終以擬合趨勢曲線水平段起始點所對應的閾值變量δth值作為最優(yōu)閾值。

對于軟閾值去噪，RM(r)曲線通常為衰減的平滑曲線，當RM(r)曲線衰減至呈較為平緩趨勢、且其值已經(jīng)足夠小時，其上某一點所對應的閾值變量δth值即可認為是最優(yōu)閾值。如果按此方式不便判斷，則進一步定義RM(r)差分譜

當d RM(r)曲線衰減至呈水平趨勢、且其值已經(jīng)足夠小時，其上某一點所對應的閾值變量δth值即可認為是最優(yōu)閾值δ*th。

需要說明的是，由改進風險函數(shù)RM(r)或者其差分譜d RM(r)確定最優(yōu)閾值δ*th的方法帶有一定的主觀性與試探性，但其有足夠的準確性，并不影響其實用性。

最后，利用步長迭代算法中估計的最優(yōu)閾值δ*th重新對信號x(t)進行S變換譜閾值去噪，再利用S逆變換重構時域信號，得到信號x(t)中的沖擊特征。

3 仿真研究

為驗證S變換譜閾值去噪方法對于噪聲混合信號沖擊特征提取的有效性，構造仿真信號x(t)進行試驗。信號x(t)中的沖擊成分p(t)由下面指數(shù)衰減模型生成，

其中，信號幅值序列 A(i)(i=1，2，…，5)，為均值 μ=2，標準差σ=0.3的高斯隨機序列，沖擊分量p(t)如圖1所示。沖擊分量疊加高斯白噪聲，得到信噪比為0的仿真信號x(t)，如圖2所示。

圖1 沖擊分量p(t)Fig.1 Impact component p(t)

對x(t)進行S變換，得到其時頻譜系數(shù)的最大模值 max{|S(j，n)|］=1.685 3。在運用步長迭代算法的過程中，取模值系數(shù)區(qū)間［α，β］=［0.1，0.6］，將區(qū)間分成M=100個步長。去噪過程中分別采用硬閾值函數(shù)和軟閾值函數(shù)，其中硬閾值函數(shù)的表達式為

圖2 仿真信號x(t)Fig.2 Simulated signal x(t)

對于硬閾值去噪，RM(r)關于步數(shù)變量r的關系曲線及其15次多項式擬合曲線如圖3所示。

圖3 仿真信號x(t)經(jīng)硬閾值去噪得到的改進風險函數(shù)RM(r)曲線及其15次多項式擬合曲線Fig.3 Modified risk function RM(r)curve obtained from the simulated signal x(t)by hard threshold denoising and its 15th degree polynomial curve

可以看出，當閾值變量δth增加的步數(shù)r=75，即δth= α·max{|S(j，n)|］+r·ΔD=0.800 5 時，RM(r)曲線出現(xiàn)水平趨勢段。因此，估計的最優(yōu)閾值=0.800 5，然后以此最優(yōu)閾值對信號x(t)重新進行S變換譜硬閾值去噪，最后利用S逆變換重構得到時域沖擊特征p'(t)，如圖4所示。

對于軟閾值去噪，RM(r)關于步數(shù)變量r的關系曲線如圖5所示。

可以認為，當閾值閾值變量δth增加的模值步數(shù)r=45，即 δth= α·max{|S［j，n］|］+r·ΔD=0.547 7 時，RM(r)曲線趨于水平。進一步地，從圖6所示的差分譜d RM(r)曲線可以看出，當r=45時，差分譜d RM(r)曲線衰減呈水平趨勢，且其值就已經(jīng)足夠小。因此，估計的最優(yōu)閾值=0.547 7，然后以此最優(yōu)閾值對信號x(t)重新進行S變換譜軟閾值去噪，最后利用S逆變換重構得到時域沖擊特征p'(t)，如圖7所示。

圖4 仿真信號x(t)經(jīng)硬閾值去噪得到的時域沖擊特征Fig.4 Time-domain impact feature obtained from the simulated signal x(t)by hard threshold denoising

圖5 仿真信號x(t)經(jīng)軟閾值去噪得到的改進風險函數(shù)RM(r)曲線Fig.5 Modified risk function RM(r)curve obtained from the simulated signal x(t)by soft threshold denoising

圖6 仿真信號x(t)改進風險函數(shù)RM(r)的差分譜d RM(r)曲線Fig.6 Difference spectrum d RM(r)curve obtained from the modified risk function RM(r)of the simulated signal x(t)

圖7 仿真信號x(t)經(jīng)軟閾值去噪得到的時域沖擊特征Fig.7 Time-domain impact feature obtained from the simulated signal x(t)by soft threshold denoising

本次仿真試驗中，S變換譜硬閾值去噪和軟閾值去噪所達到的效果基本一樣，均能夠從噪聲混合信號x(t)中提純時域沖擊特征。對比圖1、圖4和圖7，可以看出，雖然所得到的沖擊特征幅值已經(jīng)失真，這主要是由S變換中的高斯窗函數(shù)以及閾值函數(shù)所決定的，但成功提取出了作為最重要信息的沖擊特征出現(xiàn)周期，進而可以計算沖擊特征頻率，這驗證了所提出的S變換譜閾值去噪方法對于信號沖擊特征提取的有效性。

4 軸承故障信號實例分析

為進一步驗證S變換譜閾值去噪方法在實際故障振動信號中提取沖擊特征的有效性與實用性，下面采用美國Case Western Reserve University電氣工程實驗室中的相關軸承數(shù)據(jù)進行試驗［12］。試驗軸承為試驗臺驅動端的SKF 6205－2RS軸承。在試驗軸承內(nèi)滾道上用電火花加工出直徑為0.177 8 mm，深為0.279 4 mm的凹坑，以模擬軸承內(nèi)圈單點故障。驅動電機功率2.2 kW，試驗轉速1 730 r/min。理論計算該軸承內(nèi)圈故障特征頻率為156 Hz。采樣頻率為12 kHz，采樣數(shù)據(jù)長度為2 048。顯然，該軸承在運行過程中，會產(chǎn)生具有沖擊特征的振動信號，且沖擊特征出現(xiàn)頻率與為軸承內(nèi)圈故障特征頻率一致。待分析的軸承故障振動信號如圖8所示，現(xiàn)采用S變換譜閾值去噪方法將其中的沖擊特征提取出來。

圖8 軸承故障振動信號Fig.8 Bearing fault vibration signal

圖9 軸承故障振動信號經(jīng)硬閾值去噪得到的改進風險函數(shù)RM(r)曲線及其5次多項式擬合曲線Fig.9 Modified risk function RM(r)curve obtained from the bearing fault vibration signal by hard threshold denoising and its 5th degree polynomial curve

圖10 軸承故障振動信號經(jīng)硬閾值去噪得到的時域沖擊特征Fig.10 Time-domain impact feature obtained from the bearing fault vibration signal by hard threshold denoising

圖11 軸承故障振動信號經(jīng)軟閾值去噪得到的改進風險函數(shù)RM(r)曲線Fig.11 Modified risk function RM(r)curve obtained from the bearing fault vibration signal by soft threshold denoising

對該軸承故障振動信號進行S變換，得到其時頻譜系數(shù)的最大模值 max{|S(j，n)|］=0.635 8。在運用步長迭代算法的過程中，取模值系數(shù)區(qū)間［α，β］=［0.1，0.7］，并將區(qū)間分成 M=100 個步長。去噪過程中同樣分別采用硬閾值函數(shù)和軟閾值函數(shù)。

對于硬閾值去噪，RM(r)關于步數(shù)變量r的關系曲線及其5次多項式擬合曲線如圖9所示。

不難看出，當閾值閾值變量δth增加的模值步長數(shù)r=45，即 δth= α·max{|S［j，n］|］+r·ΔD=0.765 9時，RM(r)曲線出現(xiàn)水平趨勢段。因此，估計的最優(yōu)閾值=0.235 3，然后以此最優(yōu)閾值對軸承故障振動信號重新進行S變換譜硬閾值去噪，最后利用S逆變換重構得到軸承故障振動信號中的沖擊特征，如圖10所示。

對于軟閾值去噪，RM(r)關于步長數(shù)M的關系曲線如圖11所示。

可以認為，當閾值閾值變量δth增加的模值步長數(shù)r=26，即 δth= α·max{|S［j，n］|］+M·ΔD=0.1628時，RM(r)曲線衰減趨勢較為平緩。進一步地，從圖12所示的差分譜d RM(r)曲線看出，當M=26時，差分譜d RM(r)曲線趨于水平，且其值已經(jīng)很小。因此，估計的最優(yōu)閾值=0.162 8，然后以此最優(yōu)閾值對軸承故障振動信號重新進行S變換譜軟閾值去噪，最后利用S逆變換重構得到時域沖擊特征，如圖13所示。

圖12 軸承故障振動信號改進風險函數(shù)RM(r)的差分譜d RM(r)曲線Fig.12 Difference spectrum d RM(r)curve obtained from the modified risk function RM(r)of the bearing fault vibration signal

圖13 軸承故障振動信號經(jīng)軟閾值去噪得到的時域沖擊特征Fig.13 Time-domain impact feature obtained from the bearing fault vibration signal by soft threshold denoising

軸承故障振動信號的處理結果顯示，S變換譜硬閾值去噪和軟閾值去噪這兩種方法均能夠成功地從軸承故障振動信號中提取出較為純粹的時域沖擊特征。從圖10和圖13中，可以明顯看出沖擊特征出現(xiàn)的周期為 Δt=0.006 45 s，對應的頻率為 155 Hz，這與該滾動軸承內(nèi)圈的故障特征頻率156 Hz基本一致，從而可以準確作出故障診斷結論。

作為比較，現(xiàn)對軸承故障振動信號分別進行小波硬閾值去噪與軟閾值去噪。

小波硬閾值去噪的閾值采用S變換譜硬閾值去噪所用的最優(yōu)閾值=0.235 3，這樣得到軸承故障振動信號的小波硬閾值去噪結果如圖14所示。

圖14 軸承故障振動信號經(jīng)小波硬閾值去噪得到的時域沖擊特征Fig.14 Time-domain impact feature obtained from the bearing fault vibration signal by wavelet hard threshold denoising

小波軟閾值去噪的閾值采用S變換譜軟閾值去噪所用的最優(yōu)閾值=0.162 8，這樣得到軸承故障振動信號的小波軟閾值去噪結果如圖15所示。

圖15 軸承故障振動信號經(jīng)小波軟閾值去噪得到的時域沖擊特征Fig.15 Time-domain impact feature obtained from the bearing faultvibration signal by wavelet soft threshold denoising

通過比較可以明顯看出，利用步長迭代算法估計的最優(yōu)閾值，S變換譜硬、軟閾值去噪方法所提取出的軸承故障振動信號時域沖擊特征更為純粹，其中的各個沖擊特征之間沒有噪聲成分的干擾。因此，針對含噪信號時域沖擊特征的提取，相比于小波硬、軟閾值去噪，S變換譜閾值去噪方法更具有優(yōu)越性。

5 結論

(1)S變換具有多分辨率特性，滿足線性疊加原理，且對于暫態(tài)的信號沖擊特征具有較高的敏感性，非常適合于處理與分析非平穩(wěn)信號，尤其是包含沖擊特征的信號。

(2)利用所提出的改進風險函數(shù)RM(r)或者其差分譜d RM(r)，步長迭代算法能夠準確有效地從0到S變換譜系數(shù)最大模值的區(qū)間內(nèi)估計S變換譜閾值去噪所需的最優(yōu)閾值，雖然的獲取帶有一定的主觀性，但并不影響其實用性。

(3)S變換譜閾值去噪方法可以成功地從噪聲混合信號中提取出較為純粹的時域沖擊特征，且相比于小波硬、軟閾值去噪，更具有優(yōu)越性。最終根據(jù)所獲取的沖擊特征頻率，可以有效方便地實現(xiàn)滾動軸承相關故障的診斷。

［1］Lee J Y，Nandi A K.Blind deconvolution of impacting signals using higher-order statistics［J］.Mechanical Systems and Signal Processing，1998，12(2):357－371.

［2］ Lee J Y，Nandi A K.Extraction of impacting signals using blind deconvolution ［J］.Journal of Sound and Vibration，2000，232(5):945－962.

［3］李繼猛，陳雪峰，何正嘉.采用粒子群算法的沖擊信號自適應單穩(wěn)態(tài)隨機共振檢測方法［J］.機械工程學報，2011，47(21):58－63.LI Ji-meng， CHEN Xue-feng， HE Zheng-jia. Adaptive monostable stochastic resonance based on PSO with application in impact signal detection ［J］.Journal of Mechanical Engineering，2011，47(21):58 －63.

［4］章立軍，陽建宏，徐金梧，等.形態(tài)非抽樣小波及其在沖擊信號特征提取中的應用［J］.振動與沖擊，2007，26(10):56－59.ZHANG Li-jun， YANG Jian-hong， XU Jin-Wu， et al.Morphological undecimated wavelet and its application to feature extraction of impulsive signal［J］. Journal of Vibration and Shock，2007，26(10):56－59.

［5］劉小峰，柏林，秦樹人.基于自適應時頻分解的瞬態(tài)沖擊信號提純［J］.機械工程學報，2008，44(11):166－170.LIU Xiao-feng，BO Lin，QIN Shu-ren.Transient impulse signal extraction based on self-adaptive time-frequency decomposition［J］. Chinese Journal of Mechanical Engineering，2008，44(11):166 －170.

［6］ Simon C，Loubaton P H，Jutten C，et al.Separation of a class of convolutive mixtures:a contrast function approach［J］.Signal Processing，2001，81(4):883 －887.

［7］王宇，伍星，遲毅林，等.基于盲解卷積和聚類的機械弱沖擊聲信號提?。跩］.振動工程學報，2009，22(6):620－624.WANG Yu，WU Xing，CHI Yi-lin，et al.Weak transient impulse signal extraction based on blind deconvolution and cluster in acoustical machine diagnosis［J］.Journal of Vibration Engineering，2009，22(6):620 －624.

［8］楊富春，周曉軍，張志剛.基于滑動峰態(tài)算法的信號弱沖擊特征提取及應用［J］.振動與沖擊，2009，28(4):103－105.YANG Fu-chun，ZHOU Xiao-jun，ZHANG Zhi-gang.A new method for extracting weak impulse characteristic based on a sliding kurtosis algorithm ［J］.Journal of Vibration and Shock，2009，28(4):103－105.

［9］于明月，陳果.雙自適應小波局部極大模方法及其在信號特征提取中的應用［J］.振動與沖擊，2013，32(1):53－59.YU Ming-yue，CHEN Guo.Double adaptive wavelet partial maximum modulus method and its application in feature extraction of impulse signals［J］.Journal of Vibration and Shock，2013，32(1):53－59.

［10］Stockwell R G，Mansinha L，Lowe R P.Localization of the Complex Spectrum:The STransform［J］.IEEE Transactions on Signal Processing，1996，44(4):998 －1001.

［11］ Stockwell R G.Why use the S-transform? ［J］.Fields Institute Communications，2007，52:279－309.

［12］ Welcome to the Case Western Reserve University Bearing Data Center Website［EB/OL］.http://csegroups.case.edu/bearing dataconter/home