面向隨機(jī)振動(dòng)功率譜估計(jì)的小波變換去噪算法理論分析*

羅中良,汪華斌,陳治明,楊發(fā)權(quán)

(1.惠州學(xué)院電子科學(xué)系,廣東 惠州 5160072.佛山科學(xué)技術(shù)學(xué)院電子與信息工程學(xué)院，廣東 佛山 528000)

振動(dòng)試驗(yàn)的目的在于確定所設(shè)計(jì)的設(shè)備在運(yùn)輸、工作過(guò)程中能承受外來(lái)或者自身產(chǎn)生的振動(dòng)而不被破壞，能正常發(fā)揮性能，并達(dá)到預(yù)定的使用壽命。振動(dòng)大約占整個(gè)環(huán)境因素的27％,可見(jiàn)振動(dòng)試驗(yàn)在環(huán)境試驗(yàn)中的重要地位[1-3]。在可靠性振動(dòng)試驗(yàn)中，隨機(jī)振動(dòng)譜均衡控制中首先是通過(guò)對(duì)功率譜的估計(jì)，再利用均衡控制算法與參考譜進(jìn)行比較迭代，最終使控制譜收斂到參考譜。實(shí)驗(yàn)證明譜估計(jì)的誤差會(huì)引起自調(diào)整算法誤差顯著增大甚至不能收斂[4]。傳統(tǒng)基于FFT 的周期圖法對(duì)功率譜作估計(jì)由于缺少對(duì)估計(jì)值的統(tǒng)計(jì)平均，因此產(chǎn)生方差大[4-10]，韋斯(Welch)法，可以減小譜估計(jì)的誤差，卻降低了譜估計(jì)的頻率分辨率[11]。增大處理數(shù)據(jù)量提高頻率分辨率的同時(shí)造成系統(tǒng)回路時(shí)間變大，影響到系統(tǒng)的均衡性能。小波變換是一種時(shí)間-尺度分析方法，具有多分辨分析的特點(diǎn)，且具有時(shí)頻兩域表征信號(hào)局部特征的能力[6.9]，很適合于振動(dòng)頻譜均衡控制，如B.Williams[12]、LIU的研究[4]。

在通過(guò)對(duì)小波系數(shù)進(jìn)行閾值處理去除低幅度噪聲時(shí)，由于振動(dòng)譜對(duì)數(shù)域中噪聲統(tǒng)計(jì)模型在各頻率點(diǎn)是相關(guān)的，這為小波去噪方法直接應(yīng)用振動(dòng)功率譜估計(jì)提出了難題。另外，振動(dòng)功率譜通常存突變和間斷現(xiàn)象，在采用軟閾值方法時(shí)得到的估計(jì)信號(hào)在間斷處會(huì)顯得比較模糊,且整體誤差較大。而硬閾值方法在間斷點(diǎn)附近會(huì)產(chǎn)生偽Gibbs現(xiàn)象。實(shí)際上,軟、硬閾值方法在去噪和保護(hù)間斷方面處于兩個(gè)極端,若適當(dāng)將這兩種方法折衷,則能夠取得更好的去噪效果，從而提高譜均衡的整體性能。為此，本文通過(guò)對(duì)隨機(jī)振動(dòng)譜的統(tǒng)計(jì)模型進(jìn)行分析，推導(dǎo)出非線性閾值與小波變換尺度之間的關(guān)系，為濾波去噪提供依據(jù)，也為自適應(yīng)閾值算法提供理論支撐。

1 面向振動(dòng)譜估計(jì)的小波變換去噪閾值理論分析

1.1 譜估計(jì)統(tǒng)計(jì)特征及其噪聲模型

對(duì)于零均值的平穩(wěn)高斯過(guò)程{x(t),t:0,1,...,2N-1}，加w(t)窗的傅里葉變換為：

(1)

(2)

(3)

(4)

,

(5)

(6)

若信號(hào)序列{x(t),t:0,1,...,2N-1}是零均值的高斯白噪聲，方差為σ2，則式(5)可寫為：

,

(7)

對(duì)于高斯白噪聲或非高斯的平穩(wěn)隨機(jī)過(guò)程，MTSA估計(jì)器的統(tǒng)計(jì)分布為：

(8)

(9)

1.2 對(duì)數(shù)功率譜模型

1.2.1 周期圖譜估計(jì)器的對(duì)數(shù)功率譜模型 對(duì)周期圖估計(jì)器進(jìn)行對(duì)數(shù)變換，將功率譜估計(jì)值變換到對(duì)數(shù)坐標(biāo)，即：

(10)

對(duì)應(yīng)的的期望與方差為：

(11)

(12)

其中，ψ0(·)是Digamma函數(shù)；ψ1(·)是Trigamma函數(shù)；γ是Euler-Mascheroni常數(shù)，其值約為0.577 216。式(12)表明經(jīng)對(duì)數(shù)變換的PSD，其方差不再依賴PSD的真值，而是一常量。

,≤f≤1/2

(13)

其中，{u(f),0

(14)

γ=lnP(f)+ε(f), 0≤f≤1/2

(15)

(16)

由此可見(jiàn)，對(duì)數(shù)周期圖譜估計(jì)(與常數(shù)γ之和)等于真實(shí)的對(duì)數(shù)PSD與零均值、方差為π2/6的非高斯隨機(jī)變量之和，因?qū)?duì)數(shù)功率譜密度中的噪聲項(xiàng)ε(f)濾除可以提高譜估計(jì)的精度，減小譜估計(jì)的方差。

1.2.2 MTSA譜估計(jì)器對(duì)數(shù)功率譜模型 對(duì)式(8)進(jìn)行對(duì)數(shù)變換得：

(17)

對(duì)應(yīng)的期望與方差為：

(18)

(19)

類似于1.2.1節(jié)分析式(15)結(jié)果，可得：

ψ0(K)+lnK=lnP(f)+η(f),

0≤f≤1/2

(20)

(21)

當(dāng)對(duì)式(20)以傅里葉頻率進(jìn)行離散頻率采樣時(shí)，可知η(f)的傅里葉r頻率點(diǎn)的值是相關(guān)的。為了后面描述方便，現(xiàn)定義：

sη(ω)≡cov{η(f),η(f+ω)}=

cov{ln[W(f)],ln[W(f+ω)]}

(22)

sη(ω)=

(23)

1.3 噪聲小波變換系數(shù)的方差與小波變換尺度的關(guān)系

由于實(shí)際中MTSA譜估計(jì)器較周期圖譜估計(jì)器應(yīng)該廣泛的多，因此，僅討論MTSA譜估計(jì)器噪聲小波變換系數(shù)的方差與小波變換尺度的關(guān)系，從而為非線性閾值處理提供理論依據(jù)。為了表述上方便，記(20)式記為：

≡ψ0(K)+lnK=

(24)

對(duì)上中的f以傅里葉頻率進(jìn)行離散采樣，得到

(25)

其中離散傅里葉頻率fk=k/2N。

對(duì)式(25)進(jìn)行尺度為J的小波變換

(26)

以下分析如何采用閾值法將噪聲序列的小波變換系數(shù)濾除，從而得到對(duì)數(shù)功率譜lnP(fk)的估計(jì)。

由式(23)可以得到η(fk)各頻率間的協(xié)方差陣(自相關(guān)矩陣)為

(27)

(28)

在式(26)中，設(shè)gT為WJ中產(chǎn)生Εj,t的行向量，即Εj,t=bTη，由于WJ是實(shí)矩陣，有g(shù)H=gT，因此Εj,t的方差可表示為：

var{Εj,t}=gTRηg=[FTg]HDFTg

(29)

(30)

由于在同一尺度下，gT是由序列的移位而產(chǎn)生的，其傅里葉變換的模相等，因此：

(31)

上式說(shuō)明噪聲η的小波系數(shù)的方差僅與小波變換的尺度有關(guān)，而與時(shí)移參數(shù)無(wú)關(guān)，可得到與小波變換尺度相關(guān)的非線性閾值處理，其尺度為j時(shí)的閾值為

(32)

2 面向振動(dòng)譜的小波變換自適應(yīng)噪聲消除快速算法及仿真

根據(jù)以上分析，基于小波變換的MTSA譜估計(jì)非線性(自適應(yīng))閾值算法具體步驟為：

1){x(t),t:0,1,...,2N-1}為一平穩(wěn)隨機(jī)過(guò)程的2N點(diǎn)采樣序列，通過(guò)預(yù)處理使其均值為零；

2)選擇合適的K值，計(jì)算多窗口序列；

5)依據(jù)式(23)、(27)、(31)和(32)分別計(jì)算尺度j上的噪聲η的小波系數(shù)的方差及對(duì)應(yīng)的閾值；

6)選擇合適的小波系數(shù)非線性處理方法，如硬閾值處理法、軟閾值處理法或中值法對(duì)小波系數(shù)進(jìn)行非線性處理；

7)通過(guò)小波逆變換，得到對(duì)數(shù)功率譜lnP(fk)的估計(jì)。

為了檢驗(yàn)上述算法的性能，選用周期圖功率譜估計(jì)和N=1 024、K=4的多正弦窗口序列MTSA法功率譜估計(jì)以及本文算法進(jìn)行比較。在本文算法仿真時(shí)使用Coiflet小波，并采用文中推導(dǎo)的閾值在對(duì)數(shù)功率譜估計(jì)時(shí)進(jìn)行非線性去噪。三種算法對(duì)隨機(jī)振動(dòng)功率譜估計(jì)性能的比較如圖1所示，圖1(a)為經(jīng)濾波、抽樣后的2 048點(diǎn)隨機(jī)振動(dòng)時(shí)域信號(hào)，抽樣后的采樣頻率為8.192 kHz。從圖1(b)、(c)、(d)看出，周期圖譜估的功率譜波動(dòng)范圍為21.9 dB，MTSA法得到的功率譜計(jì)法得波動(dòng)范圍為11.3 dB，而本文算法得到的功率譜波動(dòng)范圍降低為6.9 dB，并保持了振動(dòng)譜的邊緣陡峭度，表明了本文方法的有效性。

3 結(jié) 論

1)構(gòu)建了小波域振動(dòng)譜估計(jì)噪聲統(tǒng)計(jì)模型及對(duì)應(yīng)的對(duì)數(shù)域振動(dòng)譜噪聲的統(tǒng)計(jì)模型；

2)在噪聲模型基礎(chǔ)上，理論推導(dǎo)了噪聲小波系數(shù)方差與尺度關(guān)系，并得到閾值與尺度的數(shù)值關(guān)系；

3)構(gòu)建了通用的小波自適應(yīng)去噪譜估計(jì)方法，在抑制譜振蕩和保護(hù)譜突變之間取得平衡；

4)通過(guò)仿真對(duì)研究成果進(jìn)行了驗(yàn)證，證明算法的有效性。

參考文獻(xiàn)：

[1]侯瑞,陳國(guó)平.振動(dòng)臺(tái)虛擬試驗(yàn)的建模和仿真研究[J].力學(xué)季刊,2008,29(2):254-258.

[2]關(guān)廣豐,叢大成,韓俊偉.隨機(jī)振動(dòng)功率譜復(fù)現(xiàn)迭代算法的研究[J].地震工程與工程振動(dòng),2006,26(6):71-76.

[3]宦海祥.電動(dòng)振動(dòng)設(shè)備的發(fā)展及展望[J].環(huán)境試驗(yàn),2006,28(4):28-31,2006.

[4]LIU Xiaoyong,YAO Xianghua,SHI Ren.Research on improving the performance of random vibration spectrum equalization using wavelet analysis[J].Mini-Micro Systemss,2004,25(12):2273-2276.

[5]LIU Xiaoyong,LUO Zhongliang,YI Mingzhu.New frequency-domain method for shock vibration control[J].中山大學(xué)學(xué)報(bào)：自然科學(xué)版,2008,47(5):49-53.

[6]周有平,羅中良.基于振動(dòng)工程相位檢測(cè)的故障診斷方法研究[J].中山大學(xué)學(xué)報(bào)：自然科學(xué)版,2009,48(2):45-50.

[7]CRISTAN A C,WALDEN A T.Multitaper power spectrum estimation and thresholding: wavelet packets versus wavelets[J].IEEE Transactions on Signal Processing,2005,50:2976-2986.

[8]羅中良,李昕.基于小波變換的紅外弱小目標(biāo)檢測(cè)新方法[J].紅外技術(shù),2006,28(8):456-559.

[9]張峰，石現(xiàn)峰.多正弦窗譜估計(jì)應(yīng)用于振動(dòng)信號(hào)頻譜分析[J].西安工業(yè)大學(xué)學(xué)報(bào),2010,30(4)：387-391.

[10]余訓(xùn)鋒,馬大瑋,魏琳.改進(jìn)周期圖法功率譜估計(jì)中的窗函數(shù)仿真分析[J].計(jì)算機(jī)仿真,2008,3,25(3):111-114.

[11]吳紅衛(wèi),吳鎮(zhèn)揚(yáng),趙鶴鳴.多正弦窗譜估計(jì)的性能分析[J].信號(hào)處理,2007,23(6):932-936.

[12]KARSHENAS A,DUNNIGAN M,WILLIAMS B.Wavelet power spectrum smoothing for random vibration control[J].IEEE Transactions on Industrial Electronics,1999,46：466-467.