基于STFT時頻譜系數(shù)收縮的信號降噪方法*

郭遠晶，魏燕定，周曉軍

(浙江大學浙江省先進制造技術重點研究實驗室 杭州，310027)

?

基于STFT時頻譜系數(shù)收縮的信號降噪方法*

郭遠晶，魏燕定，周曉軍

(浙江大學浙江省先進制造技術重點研究實驗室 杭州，310027)

針對旋轉機械故障振動信號的降噪問題，提出一種基于短時Fourier變換(short time Fourier transform,簡稱STFT)時頻譜系數(shù)收縮的信號降噪方法。先將信號進行STFT，得到其時頻譜。由于譜系數(shù)為復數(shù)，故根據(jù)模值大小進行譜系數(shù)收縮，并利用步長迭代算法在0到譜系數(shù)最大模值的區(qū)間內(nèi)估計最優(yōu)閾值。迭代運算過程中，首先，分別采用基本的硬閾值函數(shù)和軟閾值函數(shù)進行系數(shù)收縮；然后，以改進風險函數(shù)為閾值評價標準，估計最優(yōu)閾值；最后，利用最優(yōu)閾值重新進行譜系數(shù)收縮，對得到的新譜進行STFT逆變換，重構降噪后的時域信號。仿真信號與試驗數(shù)據(jù)的處理結果表明，利用所估計的最優(yōu)閾值，STFT時頻譜系數(shù)硬、軟閾值函數(shù)收縮方法均能夠實現(xiàn)噪聲混合信號的降噪。

故障診斷； 信號降噪； 短時Fourier變換； 步長迭代算法； 改進風險函數(shù)； 最優(yōu)閾值估計； 譜系數(shù)收縮

引 言

對于旋轉機械故障振動信號的降噪，較為實用有效的方法是閾值去噪，其中小波閾值去噪的研究與應用最為廣泛。Donoho等[1-2]提出了基于小波變換系數(shù)的“VisuShrink”閾值去噪方法，引入了硬閾值函數(shù)和軟閾值函數(shù)這兩種基本的系數(shù)收縮函數(shù)(系數(shù)收縮函數(shù)又可稱為閾值函數(shù))，并建立了通用的閾值估計公式。隨后，又提出了小波閾值去噪的“SUREShrink”方法[3]以及Minimax閾值估計方法[4]。國內(nèi)外學者提出了各種形式的閾值函數(shù)以及閾值估計方法[5-10]，改進了小波閾值去噪的性能，擴展了其應用領域，并研究了信號閾值去噪另一種實現(xiàn)方式的聯(lián)合時頻域去噪。文獻[11-14]借鑒了小波閾值去噪的思想，研究了基于Gabor變換時頻譜系數(shù)的閾值去噪方法。Parolai[15]針對地震圖的消噪處理，提出了基于S變換譜系數(shù)的自適應閾值函數(shù)的去噪方法，但其主要是現(xiàn)有的小波閾值函數(shù)和通用閾值估計公式的應用。Hon等[16]針對超聲彈性成像，提出了基于短時Fourier變換譜系數(shù)閾值濾波的圖像去噪算法，該算法本質(zhì)上是時頻濾波。聯(lián)合時頻域閾值去噪研究的核心和難點依然是設計合適的閾值函數(shù)以及估計最優(yōu)閾值，以期獲得最大限度接近理想無噪信號的估計信號。

筆者針對一維時域信號的降噪，提出基于STFT時頻譜系數(shù)收縮的信號降噪方法。該方法先利用STFT獲取一維時域信號的時頻譜系數(shù)，然后將譜系數(shù)按模值大小、根據(jù)閾值函數(shù)進行收縮，最后利用STFT逆變換，重構得到降噪后的時域信號。系數(shù)收縮過程中，閾值函數(shù)可以為硬閾值函數(shù)、軟閾值函數(shù)或者其他的閾值函數(shù)。對于最優(yōu)閾值估計，以所提出的改進風險函數(shù)為評價指標，采用步長迭代算法獲取。

1 STFT

給定一個時間寬度很短的窗函數(shù)γ(t)，令窗滑動，則連續(xù)時間信號x(t)的STFT定義[17]為

(1)

其中：“*”表示復數(shù)共軛。

STFT既是時間的函數(shù)，又是頻率的函數(shù)，它可以看作是信號x(t)在分析時刻τ附近的局部頻譜。

STFT是一種可逆變換。若窗函數(shù)γ(t)滿足完全重構條件

(2)

即能量歸一化條件，則STFT逆變換定義為

(3)

對于任何實際應用而言，都需將STFT(τ,f)離散化?？紤]STFT在等間隔時頻網(wǎng)格點(mT,nF)處采樣，其中T>0和F>0分別是時間變量和頻率變量的采樣周期，m和n為整數(shù)。為了簡便，記STFT(mT,nF) = STFT (m,n)，則很容易得到離散時間信號x(k)的STFT離散化形式及STFT逆變換的離散化形式

(4)

(5)

2 STFT時頻譜系數(shù)收縮的信號降噪方法

在STFT時頻譜系數(shù)收縮的信號降噪方法實現(xiàn)過程中，首先需要將采樣得到的長度為N的離散時間信號x(k)進行STFT變換，得到其時頻譜STFT(m,n)?？紤]到STFT(m,n)為復數(shù)矩陣，如果對譜系數(shù)的實部和虛部分別進行收縮，則必然改變系數(shù)的相位，從而使STFT逆變換重構得到的時域信號中產(chǎn)生附加噪聲，因此，筆者根據(jù)模值大小進行系數(shù)收縮。

1) 將目標模值區(qū)間或者模值系數(shù)目標區(qū)間分成M個步長，則每個系數(shù)步長Δd=(β-α)/M，對應的閾值步長ΔD=Δdmax{|STFT(m,n)|}。閾值變量δth初始化為αmax{|STFT(m,n)|} 。

2) 使閾值變量δth按閾值步長ΔD增大，每增大一個步長，獲得一個新閾值δnew。

3) 以δnew為閾值(初始值為αmax{|STFT(m,n)|})，根據(jù)模值大小對譜系數(shù)進行收縮。其中，閾值函數(shù)可以為硬閾值函數(shù)、軟閾值函數(shù)或者其他的閾值函數(shù)。最后對系數(shù)收縮后得到的新譜進行STFT逆變換，重構時域信號。

假設信號x(k)中原始無噪的信號為s(k)，并混有加法性噪聲u(k)，則x(k)表示為

x(k)=s(k)+u(k)

(6)

令STFT時頻譜經(jīng)過系數(shù)收縮后得到的新譜為STFT′[m,n]，那么新譜經(jīng)STFT逆變換重構得到的時域信號為

(7)

(8)

但是，在處理實際信號時， 原始無噪信號s(k)是不可知的，L2-RISK無法計算，因此，筆者提出一種改進風險函數(shù)RM(r)，其定義為

(9)

這里，xr(k)為閾值變量δth增加了r個閾值步長ΔD，即δth=αmax{|STFT(m,n)|}+rΔD時，經(jīng)過STFT時頻譜系數(shù)收縮后重構得到的時域信號。這樣，RM(r)就定義為第r步迭代計算得到的時域信號xr(k)與第r-1步迭代計算得到的時域信號xr-1(k)的均方誤差。因此，重構得到時域信號后，進一步計算改進風險函數(shù)RM(r)。

4) 重復步驟2和3，直到閾值變量δth增加到βmax{|STFT(m,n)|}，此時可獲取改進風險函數(shù)RM(r)關于步數(shù)M或閾值變量δth的關系曲線。

3 仿真研究

為驗證STFT時頻譜系數(shù)收縮方法對于信號降噪的有效性，構造仿真信號x(t)進行試驗。仿真信號x(t)由頻率為余弦信號分量s(t)疊加噪聲分量u(t)得到，其信噪比為0，即

(10)

仿真振動信號x(t)與其中的余弦信號分量s(t)分別如圖1、圖2所示。

圖1 仿真振動信號x(t)

圖2 余弦信號分量s(t)

下面對信號x(t)進行基于STFT時頻系數(shù)收縮的降噪，以獲取余弦信號分量s(t)的最佳估計信號s*(t)。

(11)

軟閾值函數(shù)的表達式為

(12)

對于硬閾值函數(shù)的STFT時頻譜系數(shù)收縮，RM(r)關于步長數(shù)M的關系曲線如圖3所示。

圖3 仿真振動信號硬閾值函數(shù)系數(shù)收縮得到的改進風險函數(shù)RM(r)曲線

對于軟閾值函數(shù)的STFT時頻譜系數(shù)收縮，RM(r)關于步長數(shù)M的關系曲線如圖5所示。

仿真試驗中，對比圖2、圖4和圖7可知，基于STFT時頻譜系數(shù)硬、軟閾值函數(shù)收縮的信號降噪方法均能夠從噪聲混合信號x(t)中獲取余弦信號分量s(t)的最佳估計信號s*(t)。雖然這兩種方法得到的s*(t)幅值存在一定程度上的失真，且后者幅值失真更為嚴重，但這主要是由閾值函數(shù)所決定

圖4 硬閾值函數(shù)系數(shù)收縮得到余弦信號分量s(t)的最佳估計信號s*(t)

Fig.4 The optimal estimated signals*(t)of cosine signal components(t) obtained by hard threshold coefficients shrinkage

圖5 仿真振動信號軟閾值函數(shù)系數(shù)收縮得到的改進風險函數(shù)RM(r)曲線

Fig.5 Emulational vibration signal modified risk functionRM(r) curve obtained by soft threshold function coefficients shrinkage

圖6 由仿真振動信號改進風險函數(shù)RM(r)獲得的差分譜dRM(r)曲線的，而作為最重要的信號頻率信息以及信號形態(tài)，兩種方法均能理想地提取。因此，仿真試驗結果驗證了基于STFT時頻譜系數(shù)收縮方法對于信號降噪的有效性。

Fig.6 Difference spectrum dRM(r) curve obtained from emulational vibration signal modified risk function RM(r)

圖7 軟閾值函數(shù)系數(shù)收縮得到的余弦信號分量s(t)最佳估計信號s*(t)

4 試驗數(shù)據(jù)降噪處理

為進一步驗證基于STFT時頻譜系數(shù)收縮方法對于實際故障振動信號降噪的有效性與實用性，采用美國PHM Society 2009年故障預測與健康管理挑戰(zhàn)賽的齒輪箱相關故障數(shù)據(jù)進行試驗[19]。以斜齒輪為試驗齒輪，選擇齒輪箱的兩種故障狀況：輸入軸彎曲和24T齒輪的一個齒上有破損。試驗時的齒輪箱輸入轉速為3 kr/min，重載。振動信號原始采樣頻率為66.666 7 kHz，再以6.666 7 kHz的頻率對采集到的信號數(shù)據(jù)進行重采樣。

對于第1種故障狀況，即輸入軸彎曲，所采集到的待分析振動信號x1(t)如圖8所示。

圖8 輸入軸彎曲故障振動信號x1(t)

信號x1(t)的STFT時頻譜系數(shù)最大模值max{|STFT(m,n)|}=0.060 3。在運用步長迭代算法的過程中，取模值系數(shù)區(qū)間[α,β]=[0.2,0.7]，并將該區(qū)間分成M=100個步長。去噪過程中同樣分別采用硬閾值函數(shù)和軟閾值函數(shù)。

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

圖9 故障振動信號x1(t)硬閾值函數(shù)系數(shù)收縮得到的改進風險函數(shù)RM(r)曲線

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

對于第2種故障狀況，即24T齒輪的一個齒上有破損，所采集到的待分析振動信號x2(t)如圖14所示。信號x2(t)的STFT時頻譜系數(shù)最大模值max{|STFT(m,n)|}=0.028 4。在運用步長迭代算法的過程中，取模值系數(shù)區(qū)間[α,β]=[0.3,0.8]，并將該區(qū)間分成M=100個步長。去噪過程中同樣分別采用硬閾值函數(shù)和軟閾值函數(shù)。

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

圖11 故障振動信號x1(t)軟閾值函數(shù)系數(shù)收縮得到的改進風險函數(shù)RM(r)曲線

Fig.11 Fault vibration signalx1(t) modified risk functionRM(r)curve obtained by soft threshold function coefficients shrinkage

圖12 由故障振動信號x1(t)改進風險函數(shù)RM(r)獲得的差分譜dRM(r)曲線

Fig.12 Difference spectrum dRM(r) curve obtained from fault vibration signalx1(t) modified risk functionRM(r)

圖13 故障振動信號x1(t)經(jīng)由軟閾值函數(shù)系數(shù)收縮得到的時域降噪信號

圖14 輪齒破損故障振動信號x2(t)

圖15 故障振動信號x2(t)硬閾值函數(shù)系數(shù)收縮得到的改進風險函數(shù)RM(r)曲線

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

圖17 故障振動信號x2(t)軟閾值系數(shù)收縮得到的改進風險函數(shù)RM(r)曲線

Fig.17 Fault vibration signalx2(t) modified risk functionRM(r) curve obtained by soft threshold function coefficients shrinkage

圖18 故障振動信號x2(t)改進風險函數(shù)RM(r)獲得的差分譜dRM(r)曲線

Fig.18 Difference spectrum dRM(r) curve obtained from bearing fault vibration signalx2(t) modified risk functionRM(r)

圖19 故障振動信號x2(t)經(jīng)由硬閾值函數(shù)系數(shù)收縮得到的時域降噪信號

5 結 論

2) STFT時頻譜系數(shù)收縮方法能夠從噪聲混合信號中恢復時域降噪信號。將此方法應用于齒輪箱相關故障振動信號的處理，可以有效提取出與故障特征信號一致的時域降噪信號，從而指導相關故障的診斷。

[1] Donoho D L， Johnstone I M． Ideal spatial adaptation by wavelet shrinkage[J]．Biometrika，1994，81(3)：425-455．

[2] Donoho D L． De-noising by soft-thresholding [J]．IEEE Transactions on Information Theory，1995，41(3)：613-627．

[3] Donoho D L，Johnstone I M．Adapting to unknown smoothness via wavelet shrinkage[J]． Journal of the American Statistical Association，1995，90(432)：1200-1224．

[4] Donoho D L，Johnstone I M．Minimax estimation via wavelet shrinkage[J]． The Annals of Statistics，1998，26(3)：879-921．

[5] 曲天書，戴逸松，王樹勛．基于SURE無偏估計的自適應小波閾值去噪[J]．電子學報，2002，30(2)：266-268．

Qu Tianshu，Dai Yisong，Wang Shuxun．Adaptive wavelet thresholding denoising method based on SURE estimation[J]．Acta Electronica Sinica，2002，30(2)：266-268．(in Chinese)

[6] Nezamoddin N K，Paul F，Edward J．BayesShrink ridgelets for image denoising [C]∥Proceedings of International Conference on Image Analysis and Recognition (ICIAR) . Verlag Berlin Heidelberg, Germany:Springer, 2004：163-170．

[7] 李方，李友榮，王志剛．基于Morlet小波與最大似然估計方法的降噪技術[J]．振動、測試與診斷，2005，25(1)：40-42．

Li Fang， Li Yourong，Wang Zhigang．De-noising technology based on Morlet wavelet transform and maximum likelihood estimation[J]．Journal of Vibration，Measurement & Diagnosis，2005，25(1)：40-42．(in Chinese)

[8] 侯新國，劉開培，魏建華．最佳小波包基改進軟閾值的消噪方法及應用[J]．振動、測試與診斷，2008，28(4)：366-368.

Hou Xinguo，Liu Kaipei，Wei Jianhua．Application of improved soft threshold noise eliminating method based on optimal wavelet packet[J]．Journal of Vibration，Measurement & Diagnosis，2008，28(4)：366-368．(in Chinere)

[9] 張弦，王宏力．進化小波消噪方法及其在滾動軸承故障診斷中的應用[J]．機械工程學報，2010，46(15)：76-81．

Zhang Xian，Wang Hongli．Evolutionary wavelet denoising and its application to ball bearing fault diagnosis[J]．Journal of Mechanical Engineering，2010，46(15)：76-81．(in Chinese)

[10]劉文藝，湯寶平，蔣永華．一種自適應小波消噪方法[J]．振動、測試與診斷，2011，31 (1)：74-77．

Liu Wenyi，Tang Baoping，Jiang Yonghua．Research on an adaptive wavelet denoising method[J]．Journal of

Vibration，Measurement & Diagnosis，2011，31(1)：74-77．(in Chinese)

[11]Xia Xianggen．System identification using chirp signals and time-variant filters in the joint time-frequency domain [J]． IEEE Transactions on Signal Processing，1997，45(8)：2072-2084．

[12]Christiansen O．Time-frequency analysis and its applications in denoising [D]．Norway：University of Bergen，2002．

[13]Nezamoddin N K，Paul F．A gabor based technique for image denoising[C]∥Canadian Conference on Electrical and Computer Engineering．Waterloo，Canada：IEEE，2005：980-983．

[14]Walker J S，Chen Y J．Denoising gabor transforms [EB/OL]．[2013-12-05] http:∥www. uwec. edu./walkerjs．

[15]Parolai S．De-noising of seismograms using the S-transform [J]．Bulletin of the Seismological Society of America，2009，99(1)：226-234．

[16]Hon T K，Subramaniam S R，Georgakis A，et al． STFT-based Denoising of Elastogram[C]∥Proceedings of IEEE International Conference on Acoustics，Speech and Signal Processing(ICASSP)． Prague，Czech Republic：IEEE，2011：677-680．

[17]張賢達，保錚．非平穩(wěn)信號分析與處理[M]．北京：國防工業(yè)出版社，1998：20-25．

[18]郭遠晶，魏燕定，周曉軍，等.基于S變換譜閾值去噪的沖擊特征提取方法[J].振動與沖擊，2014，33(21)：44-50.

Guo Yuanjing, Wei Yanding, Zhou Xiaojun, et al. An impact feature extracting method based on S-transformation spectrum threshold denoising[S].Journal of Vibration and Shock,2014,33(21):44-50.(in Chinese)

[19]IEEE PHM Society.2009 PHM challenge competition data set[EB/OL]．http:∥[2013-12-05].www.phmsociety.org/references/datasets．

10.16450/j.cnki.issn.1004-6801.2015.06.014

*國家自然科學基金資助項目(51275453)

2013-10-23；

2013-12-20

TN911.7； TH133； TH165.3

郭遠晶，男，1987年10月生，博士研究生。主要研究方向為機械動力學、機械信號處理與故障診斷。E-mail:gyjyn@126.com