雙自尋優(yōu)小波去噪方法及其在滾動(dòng)軸承故障診斷中的應(yīng)用

湯 偉, 肖方煜, 陳 曦

(陜西科技大學(xué) 自動(dòng)化應(yīng)用研究所， 陜西 西安 710021)

0 引言

滾動(dòng)軸承在各類機(jī)械中應(yīng)用很廣泛，且其運(yùn)行狀態(tài)直接影響整臺(tái)機(jī)器的性能、壽命、功能、功率及安全.在旋轉(zhuǎn)機(jī)械中滾動(dòng)軸承發(fā)生故障的概率很高，據(jù)統(tǒng)計(jì)大約有30% 的機(jī)械故障都是由滾動(dòng)軸承引起的[1].但在實(shí)際的故障診斷過(guò)程中，所獲取的滾動(dòng)軸承故障信號(hào)往往夾雜著各種各樣的噪聲，給故障的判斷造成了很大的影響.所以在診斷故障前，必須去除噪聲.小波閾值去噪方法由于原理簡(jiǎn)練、可操作性強(qiáng)與去噪效果顯著的特點(diǎn)，已在信號(hào)去噪領(lǐng)域中獲得了成功應(yīng)用.

小波閾值去噪方法中，決定去噪效果最為關(guān)鍵的兩個(gè)因素是閾值和閾值函數(shù)的選取[2].但是如何選取閾值和閾值函數(shù)，并且它們之間如何搭配，使其去噪的效果最優(yōu)成為小波閾值去噪方法最關(guān)鍵的問(wèn)題[3].本文正是基于這個(gè)問(wèn)題，提出一種閾值和閾值函數(shù)各自尋優(yōu)又相互結(jié)合的選取策略予以解決.通過(guò)模擬信號(hào)消噪仿真與滾動(dòng)軸承故障信號(hào)消噪仿真實(shí)例表明，該方法無(wú)論在信噪比和去噪圖像的直觀效果上都優(yōu)于傳統(tǒng)的去噪方法.

1 現(xiàn)有閾值和閾值函數(shù)選取方法簡(jiǎn)介

關(guān)于閾值和閾值函數(shù)的選取思想分為3類：第一，閾值和閾值函數(shù)的選取相互獨(dú)立.第二，閾值的選取基于閾值函數(shù).第三，閾值函數(shù)的選取基于閾值.第一類中閾值和閾值函數(shù)相互獨(dú)立選取策略忽略了閾值和閾值函數(shù)之間耦合的關(guān)系.因?yàn)楦鶕?jù)風(fēng)險(xiǎn)函數(shù)的定義[4]:

(1)

去噪的目標(biāo)是使估計(jì)信號(hào)與真實(shí)信號(hào)的風(fēng)險(xiǎn)最小化[5].而估計(jì)信號(hào)是含噪信號(hào)通過(guò)由閾值和閾值函數(shù)映射處理后得到的，因此可以把閾值和閾值函數(shù)當(dāng)做一個(gè)耦合的系統(tǒng).第一類思想的選取方法忽略了兩者之間的耦合關(guān)系，因此很難使去噪效果達(dá)到最優(yōu)化.第二類思想有基于某種閾值函數(shù)的閾值選取策略，代表有基于史坦的無(wú)偏似然估計(jì)方法、文獻(xiàn)[6]提出的基于硬閾值函數(shù)的閾值選取策略等.這種策略認(rèn)同了閾值和閾值函數(shù)的耦合關(guān)系，會(huì)提高去噪的效果.但是，每種閾值函數(shù)都有它們的缺點(diǎn)，文獻(xiàn)[7]中證明硬閾值函數(shù)存在較大的方差，軟閾值函數(shù)存在較大的偏差，基于第二類思想選取策略勢(shì)必也存在著缺點(diǎn)，會(huì)影響去噪的效果.基于閾值函數(shù)存在的這些缺點(diǎn)，很多學(xué)者創(chuàng)造出了含有變量參數(shù)的閾值函數(shù)，典型的有半軟閾值函數(shù)、指數(shù)型閾值函數(shù)等.第三類方法是在這些含有參數(shù)的閾值函數(shù)中建立起來(lái)的，通過(guò)某一確定閾值選取策略，對(duì)閾值函數(shù)按照試湊或自尋優(yōu)的方法得到能使去噪效果達(dá)到最優(yōu)的參數(shù)[8].

上述的3種方法在實(shí)際應(yīng)用中都展露出了較為出色的效果，但這3種方法都沒(méi)能把閾值和閾值函數(shù)的選取完全聯(lián)合在一起.不僅從閾值和閾值函數(shù)之間耦合的角度來(lái)看存在著片面性，而且從數(shù)與數(shù)的組合(閾值和閾值函數(shù)的參數(shù)都是具體的數(shù)值)角度來(lái)看也存在片面性，因?yàn)樯鲜?類方法中皆都運(yùn)用數(shù)理統(tǒng)計(jì)決策論的方法估算閾值和閾值函數(shù)，這樣從數(shù)與數(shù)的組合角度來(lái)看估算到最優(yōu)的搭配的概率極低，會(huì)對(duì)產(chǎn)生最優(yōu)去噪效果埋下伏筆.

2 雙自尋優(yōu)方法

雙自尋優(yōu)方法的目標(biāo)是所有閾值和所有閾值函數(shù)參數(shù)都能互相結(jié)合，從各種結(jié)合中尋找出最優(yōu)的搭配.這樣的好處不僅解決了耦合問(wèn)題，而且也涵蓋了所有數(shù)與數(shù)之間組合性的所有可能性.

2.1 閾值自尋優(yōu)區(qū)間及高頻死去值現(xiàn)象

閾值是一個(gè)不小于零的數(shù)[9]，它的自尋優(yōu)區(qū)間是從0到一個(gè)未知數(shù)，數(shù)學(xué)表示為[0,x].在確定自尋優(yōu)區(qū)間這個(gè)未知數(shù)之前，先介紹在試驗(yàn)中得到的一個(gè)現(xiàn)象--高頻死去值現(xiàn)象[10].

高頻死去值現(xiàn)象的表述：小波閾值去噪方法中，含噪信號(hào)經(jīng)過(guò)小波分解后，對(duì)任意一層閾值使其從零按某一較小步長(zhǎng)遞增，當(dāng)閾值增大到某一數(shù)值之后，一切衡量去噪效果的數(shù)值(信噪比等)不會(huì)再隨閾值的增大而改變，定義此閾值為這層的高頻死去值.

證明：設(shè)vj表示小波分解后信號(hào)的低頻部分，Wj表示對(duì)應(yīng)的高頻部分，{vj}j∈z是L2(R)(L2(R)是全體能量有限信號(hào)的集合)的一個(gè)OMRA，φ(t)是相應(yīng)的生成元，相應(yīng)的雙尺度方程為：

(2)

(3)

φjk(t)=2j/2φ(2jt-k)

(4)

ψjk(t)=2j/2ψ(2jt-k)

(5)

而Vj與Wj分別是由{φjk(t),k∈z}和{ψjk(t),k∈z}生成，即

(6)

(7)

設(shè)fj(t)和wj(t)分別是f(t)在vj和wj中的逼近與細(xì)節(jié)，即

(8)

(9)

因?yàn)閂j+1=Vj⊕Wj,Vj⊥Wj,所以

cj,k=〈fj,φjk〉=〈fj+wj,φjk〉=〈fj+1,φjk〉

(10)

而

(11)

代入(10)式，得

(12)

同理，得

(13)

經(jīng)過(guò)分解后信號(hào)有如下表達(dá)式：

(14)

假定對(duì)fJ(t)在第j空間上的細(xì)節(jié)小波分量wj(t) 進(jìn)行映射處理，也就是對(duì)dj,k進(jìn)行閾值處理.因?yàn)閐j,k是一個(gè)具體的值，所以在第j空間上必定可以找到一個(gè)閾值λ，使λ≥dj,k.因此，定義λ=dj,k的λ值為高頻死去值.

根據(jù){φj+1,k(t),k∈z}的正交性，得

cj+1,k=〈fj+1,φj+1,k〉=〈fj,φj+1,k〉+〈wj,φj+1,k〉

(15)

因?yàn)檫x擇高頻死去值λ=dj,k作為閾值，使得dj,k被置零，即〈wj,φj+1,k〉 項(xiàng)為零.因此

cj+1,k=〈fj+1,φj+1,k〉

(16)

證畢.

2.2 高頻死去值的獲取及閾值自尋優(yōu)區(qū)間的確定

因?yàn)椴煌撝岛瘮?shù)的映射性質(zhì)不同，相同的含噪信號(hào)經(jīng)過(guò)不同閾值函數(shù)產(chǎn)生的信噪比數(shù)值也就不同.根據(jù)高頻死去值性質(zhì)中大于高頻死去值后的小波系數(shù)不含有信號(hào)成分這一特性，可以使用兩個(gè)不同閾值函數(shù)進(jìn)行自尋優(yōu)搜索，使它們自尋優(yōu)得到的信噪比數(shù)值相減.當(dāng)相減的結(jié)果在某一值之后為零，則這個(gè)值就為高頻死去值.圖1為基于同一含噪信號(hào)兩種不同閾值函數(shù)消噪的信噪比-閾值示意圖.圖2為經(jīng)過(guò)matlab的仿真數(shù)據(jù)圖，此圖是基于Blocks[11]信號(hào),小波基為sym8，信噪比為4 db，在第一層時(shí)模平方閾值函數(shù)和軟閾值函數(shù)產(chǎn)生的信噪比數(shù)值矩陣相減結(jié)果圖.此圖表明當(dāng)閾值小于54時(shí)，其差值曲線一直變化，當(dāng)大于54之后，其差值曲線為零.又因?yàn)樗x的步長(zhǎng)為十分位，所以這個(gè)信號(hào)在第一層的高頻死去值為5.4.

圖1 基于同一含噪信號(hào)兩種不同閾值函數(shù)消噪的信噪比-閾值圖

圖2 基于不同閾值函數(shù)的信噪比相減結(jié)果圖

2.3 基于軟硬閾值折衷法閾值函數(shù)的自尋優(yōu)區(qū)間的確定

針對(duì)含有參數(shù)的閾值函數(shù)的尋優(yōu)過(guò)程方法比較簡(jiǎn)單，只需對(duì)其參數(shù)在某一區(qū)間進(jìn)行自尋優(yōu)的過(guò)程.例如針對(duì)軟硬閾值折衷法閾值函數(shù)[12].

其定義為：

(17)

參數(shù)α的自尋優(yōu)區(qū)間是[0,1]，假設(shè)以0.1為其步長(zhǎng)進(jìn)行尋優(yōu)過(guò)程，因此經(jīng)過(guò)11次去噪算法后將得到11個(gè)衡量去噪效果的數(shù)值(例如信噪比)，在這十一個(gè)數(shù)值中經(jīng)過(guò)大小的比較，可以得到一個(gè)最優(yōu)的數(shù)值.

因?yàn)檐涢撝岛瘮?shù)代表收縮力最強(qiáng)的閾值函數(shù)，而硬閾值函數(shù)代表保留性最強(qiáng)的閾值函數(shù)[13].因此，閾值函數(shù)中參數(shù)選取的區(qū)間范圍，是使閾值函數(shù)能從軟閾值函數(shù)變化到硬閾值函數(shù)為標(biāo)準(zhǔn)的.

2.4 雙自尋優(yōu)方法程序的實(shí)現(xiàn)

整個(gè)雙自尋優(yōu)的過(guò)程方法如下.

首先，根據(jù)確定的分解層數(shù)，運(yùn)用2.2節(jié)的方法確定最高分解層的高頻死去值.由此確定閾值的自尋優(yōu)范圍是0到最高分解層高頻死去值.(因?yàn)榉纸鈱訑?shù)越高，則對(duì)應(yīng)的高頻部分能量越大[14]，則高頻死去值也越大，所以這里是最高層的高頻死去值)

其次，確定含有參數(shù)的軟硬閾值折衷法閾值函數(shù)的自尋優(yōu)區(qū)間.這個(gè)范圍是[0,1].

最后，編寫(xiě)整個(gè)自尋優(yōu)去噪程序.

程序如下(本文程序?yàn)閙atlab語(yǔ)言程序,且以分解層數(shù)為一層,含有參數(shù)的閾值函數(shù)是軟硬閾值折衷法來(lái)編寫(xiě))：

for i=0∶1:最高層高頻死去值(閾值自尋優(yōu)過(guò)程)

for j=0∶1：已知值(閾值函數(shù)自尋優(yōu)過(guò)程)

t=i/步長(zhǎng)1；

u=j/步長(zhǎng)2；

e=含噪信號(hào);e=e';

[c,l]=wavedec(e,1,'sym8');

a1=appcoef(c,l,'sym8',1);

d1=detcoef(c,l,1);

softd1=wthresh(d1,'p',t);

c2=[a1 softd1];

s3=waverec(c2,l,'sym8');

p1=1/length(e)*norm(e)^2;

p2h=1/length(e)*norm(e-s3)^2;

snrh=10*log(p1/p2h);

rmseh=sqrt(p2h);

plot(s3),axis([0 1200 -3 3]);

snr= rmse=[rmseh]

end

end

其中wthresh調(diào)用的程序是

function y = wthresh(x,sorh,t)

switch sorh

case 'p'

tmp=(abs(x)-u*t);

tmp=(tmp+abs(tmp))/2;

y=sign(x)*tmp;

最后，由max(max(snr)就可得到衡量去噪效果最優(yōu)的信噪比值，并可以對(duì)應(yīng)的找出具體的u和t是多少.從而得到最優(yōu)的去噪效果.

3 模擬信號(hào)去噪

通過(guò)構(gòu)造一個(gè)數(shù)據(jù)長(zhǎng)度為1 024且受到強(qiáng)噪聲污染的沖擊信號(hào)，以模擬滾動(dòng)軸承故障信號(hào)，利用Symmlet 8小波對(duì)含噪沖擊信號(hào)進(jìn)行3尺度靜態(tài)小波分解.本文的自尋優(yōu)算法將分別基于軟閾值函數(shù)、硬閾值函數(shù)與經(jīng)典的Heursure、Sureshrink、Visushrink、Mnimaxi方法進(jìn)行比較，比較結(jié)果在表1中，表1所得數(shù)據(jù)為模擬信號(hào)與所選閾值與閾值函數(shù)共同決定的信噪比,單位為db；產(chǎn)生的圖形如圖3所示.

表1 所得信噪比試驗(yàn)結(jié)果

圖3 模擬信號(hào)去噪結(jié)果圖形的比較

4 滾動(dòng)軸承故障信號(hào)去噪

試驗(yàn)數(shù)據(jù)由美國(guó)凱斯西儲(chǔ)大學(xué)電氣工程實(shí)驗(yàn)室網(wǎng)站下載[15].試驗(yàn)裝置由一個(gè)1.5 kW 電動(dòng)機(jī)、一個(gè)扭矩傳感器/譯碼器、一個(gè)功率測(cè)試計(jì)與相應(yīng)電器控制裝置組成.實(shí)驗(yàn)中使用加速傳感器采集振動(dòng)信號(hào)，加速度傳感器安裝在風(fēng)扇端6點(diǎn)鐘方向.振動(dòng)信號(hào)通過(guò)16通道的DAT記錄器采集的MAT格式，數(shù)字信號(hào)的采集信號(hào)為120 000 S/s.采用6205-2RS JEM SKF 型軸承，軸承的內(nèi)徑直徑為25.001 2 mm ，外徑直徑為51.998 9 mm ，厚度為0.590 6 mm，節(jié)徑為39.039 8 mm，滾動(dòng)體直徑為15.001 2 mm，其外圈損傷直徑為0.177 8 mm，故障特征頻率為104 Hz.選擇時(shí)間長(zhǎng)度1 s(對(duì)應(yīng)的是1 200個(gè)數(shù)據(jù)點(diǎn))的故障數(shù)據(jù)，如圖4所示.

圖4 原始信號(hào)

圖4所顯示的信號(hào)難以判斷其表示故障特征的沖擊信號(hào).圖5為本文方法去噪后的故障信號(hào)，其去噪的時(shí)域圖形中故障沖擊信號(hào)特征明顯.可以清晰的分辨出由故障引起的沖擊信號(hào)大約為115個(gè)點(diǎn)位周期出現(xiàn)，這與故障特征頻率為104 Hz(104 Hz對(duì)應(yīng)的點(diǎn)為115個(gè)點(diǎn))比較一致.特別是在400～600處的點(diǎn)和1 100～1 200處的點(diǎn)信號(hào)恢復(fù)的情況尤為明顯，可以明顯辨認(rèn)出故障沖擊信號(hào).圖6為經(jīng)minimaxi與模平方閾值函數(shù)去噪后的故障信號(hào).500～600點(diǎn)處的信號(hào)恢復(fù)較好但600～700和1 100～1 200點(diǎn)處的信號(hào)仍然難以辨認(rèn).圖7為經(jīng)rigrsure和軟閾值函數(shù)去噪后的故障信號(hào)，600～700和1 100～1 200點(diǎn)處的信號(hào)難以辨認(rèn).圖8為經(jīng)rigrsure和軟硬閾值折衷法函數(shù)去噪后的故障信號(hào).其較好的得到了故障特征信號(hào)，但在1 100～1 200點(diǎn)處的信號(hào)仍然難以辨認(rèn).

圖5 本文方法

圖6 minimaxi與模平方閾值函數(shù)方法

圖7 rigrsure和軟閾值函數(shù)方法

圖8 rigrsure和軟硬閾值折衷法函數(shù)

上述三種方法中，第一種minimaxi與模平方閾值函數(shù)的結(jié)合是本文第一章中提到的第一類方法，閾值與閾值函數(shù)相互獨(dú)立的選取方法.因?yàn)殚撝蹬c閾值函數(shù)相互獨(dú)立的選取方法忽略了兩者之間的耦合關(guān)系因此去噪效果不佳.第二種rigrsure和軟閾值函數(shù)的結(jié)合是第二類方法，因?yàn)檐涢撝岛瘮?shù)自身存在的偏差大的特點(diǎn)使得去噪效果也不是很好.第三種rigrsure和軟硬閾值折衷法函數(shù)結(jié)合的方法屬于第三類方法，因?yàn)檫@種方法只是一個(gè)因素在變化，所以去噪效果仍然不是很好.

5 結(jié)束語(yǔ)

閾值和閾值函數(shù)是小波閾值去噪方法中最重要的兩個(gè)因素.為改善小波閾值去噪方法的去噪性能，本文從傳統(tǒng)方法中歸納其優(yōu)缺點(diǎn)，并從閾值和閾值函數(shù)相互耦合的思想出發(fā)，運(yùn)用雙自尋優(yōu)方法對(duì)含噪故障信號(hào)進(jìn)行去噪處理.從模擬的故障沖擊信號(hào)和來(lái)自西儲(chǔ)大學(xué)的滾動(dòng)軸承信號(hào)實(shí)例表明，基于軟硬閾值折衷法的雙尋優(yōu)方法無(wú)論在信噪比還是得到的恢復(fù)圖形中皆優(yōu)于傳統(tǒng)小波閾值去噪方法，具有一定的理論研究與工程實(shí)際應(yīng)用價(jià)值.

[1] 鐘秉林.黃 仁.機(jī)械故障診斷學(xué)[M].北京：機(jī)械工業(yè)出版社，2007.

[2] 張德豐.Matlab小波分析與工程應(yīng)用[M].北京：國(guó)防工業(yè)出版社，2008：42-49.

[3] D.L.Donoho,I.M.Johnstone,Ideal spatial adaptation by wavelet shrinkage[J].Biometrika,1994,81(12):425-455.

[4] Ronggen Yang.Wavelet denoising using principal component analysis[J].Expert Systems with Applications,2011,11(38):1 073-1 076.

[5] 葛哲學(xué).小波分析理論與MATLAB R2007實(shí)現(xiàn)[M].北京：電子工業(yè)出版社，2007：67-69.

[6] 趙治棟，潘 敏.小波收縮中統(tǒng)一閾值函數(shù)及其偏差、方差與風(fēng)險(xiǎn)分析[J].電子與信息學(xué)報(bào)，2005,27(4):536-540.

[7] 張 弦.進(jìn)化小波消噪方法及其在滾動(dòng)軸承故障診斷中的應(yīng)用[J].機(jī)械工程學(xué)報(bào)，2010,46(15):76-81.

[8] Hilton de Oliveira Mota.Partial discharge signal denoising with spatially adaptive wavelet thresholding and support vector machines[J].Electric Power Systems Resarch,2011,23(81):644-659.

[9] D.L.Donoho,I.M.Johnstone,Wavelet Shrinkage:Asymptopia[J].Journal of the Royal Statistical Society,1995,46(57):301-369.

[10] 肖方煜，湯 偉.自尋優(yōu)閾值小波去噪方法[J].信號(hào)處理，2012,28(4)：577-586.

[11] 祝海龍,郭天佑,屈梁生.基于二進(jìn)小波變化和軟閾值改進(jìn)的信號(hào)消噪[J].自動(dòng)化學(xué)報(bào),2004,3(2):199-206.

[12] L.-K.Shark,Gear fault detection using customized multiwavelet lifting schemes[J]. Mechanical Systems and Signal Processing,2010,24(5):1 509-1 528.

[13] Andrew G.Bruce.Understanding waveshrink:variance and bias estimation[J].Biometrika,1996,83(4):727-745.

[14] 孟志強(qiáng),王鼎順,周華安.小波圖像消噪中閾值與信噪比的單峰規(guī)律及閾值試探方法[J].電子測(cè)量與儀器學(xué)報(bào),2006,20(5):63-67.

[15] The Case Western Reserve University Bearing Data Center Website[EB/OL].http：//www.eecs.cwru.edu/laboratory/bearing/.