向北平， 周 建， 倪 磊， 艾攀華

(西南科技大學(xué)制造過(guò)程測(cè)試技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室 綿陽(yáng)，621000)

引 言

長(zhǎng)期以來(lái)，小波包濾波是對(duì)信號(hào)進(jìn)行分析預(yù)處理的主要工具，該方法對(duì)信號(hào)進(jìn)行小波包變換，然后利用小波包閾值函數(shù)對(duì)小波包系數(shù)進(jìn)行閾值收縮處理以達(dá)到去噪目的。然而由于含噪信號(hào)中的噪聲分布往往不均勻，傳統(tǒng)的軟、硬閾值方法對(duì)信號(hào)小波包系數(shù)進(jìn)行固定格式的閾值處理，無(wú)法很好地滿足信號(hào)去噪要求[1]。另外，在閾值的選擇方面，常用的Heursure閾值、Donoho閾值等能夠?qū)π旁氡容^高的信號(hào)實(shí)現(xiàn)噪聲與信號(hào)的最優(yōu)分離，而對(duì)于強(qiáng)噪聲信號(hào)，去噪效果并不理想[2]。近年來(lái)，針對(duì)這兩個(gè)問(wèn)題，許多學(xué)者進(jìn)行了研究。

文獻(xiàn)[3]提出了一種自適應(yīng)對(duì)數(shù)小波閾值函數(shù)去噪算法，結(jié)合自適應(yīng)對(duì)數(shù)閾值函數(shù)對(duì)每一層小波系數(shù)設(shè)置最優(yōu)閾值，且應(yīng)用于動(dòng)壓信號(hào)，增加了信號(hào)信噪比且減少了計(jì)算時(shí)間，但其去噪算子形式依然固定不變。文獻(xiàn)[4]分析了傳統(tǒng)軟、硬閾值方法的局限，對(duì)閾值函數(shù)進(jìn)行改進(jìn)，使其具有能量分布自適應(yīng)性，但信號(hào)小波包系數(shù)的能量分布并不能準(zhǔn)確表達(dá)小波包系數(shù)的噪聲分布，因此該方法仍然存在一些不足。文獻(xiàn)[5]對(duì)Donoho閾值方法中的噪聲標(biāo)準(zhǔn)差估計(jì)方法進(jìn)行改進(jìn)，進(jìn)行仿真且取得了良好的去噪效果。Richman等[6]提出一種新的時(shí)間序列復(fù)雜度表征參數(shù)即樣本熵(sample entropy,簡(jiǎn)稱SE)算法，被廣大學(xué)者所關(guān)注且近年來(lái)被常用于機(jī)械信號(hào)分析與故障診斷領(lǐng)域[7]。

基于以上分析，筆者提出將樣本熵與小波包閾值去噪算法相結(jié)合，且將其應(yīng)用于高速深溝球軸承振動(dòng)信號(hào)去噪分析，通過(guò)分析信號(hào)的小波包系數(shù)噪聲分布情況及其對(duì)應(yīng)的樣本熵值，使閾值去噪算子具有噪聲分布自適應(yīng)性以達(dá)到最優(yōu)去噪效果；同時(shí)以噪聲估計(jì)信號(hào)樣本熵值為基準(zhǔn)，提出了一種最優(yōu)閾值估計(jì)方法。仿真分析以及實(shí)驗(yàn)結(jié)果皆對(duì)所提方法進(jìn)行了驗(yàn)證。

1 樣本熵算法及其噪聲表征

1.1 樣本熵的計(jì)算

設(shè)有長(zhǎng)度為N的時(shí)間序列Xi={x1,x2,…,xN}，其樣本熵的計(jì)算方法如下：

1) 確定嵌入維數(shù)為m，對(duì)Xi的元素按順序進(jìn)行排列，即可得到一組維數(shù)為m的向量{xm(1),…,xm(Ν-m+1)}，且

Xm(i)={x(i),x(i+1), …,x(i+m-1)}

(1≤i≤N-m+1)

(1)

2) 定義向量Xm(i)與Xm(j)之間的間隔d[Xm(i),Xm(j)]為兩向量之間對(duì)應(yīng)元素求差的絕對(duì)值的最大值，即

d[Xm(i),Xm(j)]=

maxk=0,…,m-1(|x(i+k)-x(j+k)|)

(2)

3) 對(duì)于固定的Xm(i)，統(tǒng)計(jì)Xm(i)與Xm(j)之間距離小于等于參數(shù)r的j(1≤j≤N-m,j≠i)的個(gè)數(shù)，且記為Bi，則當(dāng)1≤i≤N-m時(shí)定義

(3)

4) 定義B(m)(r)為

(4)

(5)

6) 定義Am(r)為

(6)

據(jù)上述分析可知，B(m)(r)是兩個(gè)序列在相似容限r(nóng)下匹配m個(gè)點(diǎn)的概率，而Am(r)是兩個(gè)序列匹配m+1個(gè)點(diǎn)的概率。則該時(shí)間序列樣本熵定義為

(7)

實(shí)際信號(hào)中N無(wú)法趨近于無(wú)窮，因此可將樣本熵設(shè)為

(8)

上述算法中的嵌入維數(shù)m和相似容限r(nóng)通常取為m=1或2，r=0.1，Std～0.25Std(Std為原始數(shù)據(jù)Xi={x1,x2,…,xN}的標(biāo)準(zhǔn)差)。文中取m=2，r=0.2Std。

1.2 樣本熵與噪聲

當(dāng)信號(hào)受到噪聲干擾后，其狀態(tài)取值不確定性增加，即信號(hào)無(wú)序程度與復(fù)雜程度增加。而樣本熵作為信號(hào)復(fù)雜度表征參數(shù)，當(dāng)信號(hào)中噪聲增加，樣本熵值也應(yīng)增加。為驗(yàn)證以上分析，設(shè)定1 kHz采樣率與1 000總采樣數(shù)的仿真信號(hào)f(t)=y(t)+n(t)，式中n(t)為標(biāo)準(zhǔn)差σ∈[0,3]的帶限500 Hz高斯白噪聲信號(hào)。

y(t)=5sin(10πt)sin(2πt) (0≤t≤1)

(9)

分析信號(hào)f(t)的樣本熵值在不同采樣數(shù)下與噪聲標(biāo)準(zhǔn)差的關(guān)系，其結(jié)果如圖1(a)所示(圖中100采樣點(diǎn)為信號(hào)前100采樣點(diǎn)，下同)。由于實(shí)際工程信號(hào)中的噪聲很少為單純的白噪聲，因此，利用反冪律濾波器對(duì)功率譜密度分布均勻的白噪聲n(t)上色，選定反頻譜指數(shù)為1，此時(shí)理想數(shù)字濾波器的幅度平方響應(yīng)為1/f1。由于該濾波器處理后導(dǎo)致噪聲幅值衰減，因此，對(duì)加色后的噪聲信號(hào)進(jìn)行10倍處理，得到有色噪聲n1(t)及含噪仿真信號(hào)f1(t)=y(t)+n1(t)。設(shè)n2(t)為1 kHz采樣率，1 000采樣數(shù)的幅值區(qū)間為[0,3]的帶限500Hz周期性隨機(jī)噪聲，得到含噪仿真信號(hào)f2(t)=y(t)+n2(t)。同樣分析信號(hào)f1(t)與f2(t)的樣本熵值隨噪聲大小的變化關(guān)系，其結(jié)果分別如圖1(b)，(c)所示。

圖1 樣本熵與幾種常見(jiàn)噪聲大小關(guān)系Fig.1 Correlation between sample entropy and several common noise

由圖1(a)看出，信號(hào)的樣本熵值與噪聲標(biāo)準(zhǔn)差成正相關(guān)，且在采樣點(diǎn)相差巨大的情況下其樣本熵值變化依然相近。圖1(b)，(c)較圖1(a)而言，樣本熵隨噪聲標(biāo)準(zhǔn)差變化有少許波動(dòng)，熵值大小區(qū)間有變化，但總體趨勢(shì)仍呈正相關(guān)，且數(shù)據(jù)長(zhǎng)度對(duì)其趨勢(shì)影響不大。這說(shuō)明，雖然對(duì)仿真信號(hào)施加不同種類的噪聲，其樣本熵變化曲線有所不同，但由于對(duì)信號(hào)增加噪聲的結(jié)果導(dǎo)致了信號(hào)的隨機(jī)度與復(fù)雜度增加，其結(jié)果仍然是樣本熵增加。圖1所示說(shuō)明了樣本熵算法可用來(lái)表征時(shí)間序列的含(多種)噪聲情況且受數(shù)據(jù)長(zhǎng)度影響較小。

2 改進(jìn)小波包閾值去噪算法

2.1 傳統(tǒng)小波包去噪及其不足

傳統(tǒng)的閾值函數(shù)主要有軟閾值函數(shù)和硬閾值函數(shù)兩種。軟閾值函數(shù)

(10)

硬閾值函數(shù)

(11)

工程實(shí)踐中常用的閾值函數(shù)還有形如文獻(xiàn)[8]中提出的一種介于軟、硬閾值函數(shù)之間的改進(jìn)閾值函數(shù)

(12)

由以上函數(shù)式可知：硬閾值由于其函數(shù)不連續(xù)性可能產(chǎn)生偽Gibbs現(xiàn)象且導(dǎo)致有用信號(hào)缺失；軟閾值函數(shù)雖然連續(xù)性好，但存在較大偏差[9-10]。文獻(xiàn)[8]閾值函數(shù)雖進(jìn)行了改進(jìn)，但只是軟硬閾值函數(shù)的折中算法，通過(guò)下文分析可知，該閾值函數(shù)僅為本研究提出的自適應(yīng)閾值函數(shù)的一種取值情況(即當(dāng)調(diào)節(jié)參數(shù)s=0.5時(shí))。在實(shí)際應(yīng)用中，閾值函數(shù)的選擇并無(wú)固定標(biāo)準(zhǔn)，對(duì)于不同的信號(hào)選用不同的去噪算子達(dá)到的去噪效果也不同[11]，因此，研究一種能依據(jù)信號(hào)小波包系數(shù)的含噪情況自適應(yīng)調(diào)整的新閾值函數(shù)具有實(shí)際的工程意義。

2.2 自適應(yīng)改進(jìn)閾值函數(shù)

根據(jù)以上分析，筆者對(duì)傳統(tǒng)閾值函數(shù)作改進(jìn)，使其不受限于固定去噪形式，得到的改進(jìn)閾值函數(shù)為

(13)

其中：0

由于閾值函數(shù)偏硬時(shí)對(duì)大于閾值的小波包系數(shù)保留較好，而偏軟時(shí)對(duì)小波包系數(shù)壓縮更大。因此，對(duì)于含噪較多的小波包系數(shù)，閾值處理方式應(yīng)偏軟即s較大；而受噪聲影響小的小波包系數(shù)閾值函數(shù)應(yīng)偏硬，即s應(yīng)較小。

由1.2節(jié)分析可知，樣本熵能較好地反應(yīng)時(shí)間序列的噪聲變化情況，且適用于分析短時(shí)間序列，因此設(shè)定調(diào)節(jié)參數(shù)s的確定方法如下：

1) 對(duì)信號(hào)的小波包系數(shù)序列(w1,w2,…,wn)按順序分割為n-l個(gè)相互之間最大重疊且長(zhǎng)度均為l的子序列(k1,k2,…,kn-l)，即ki向后移動(dòng)一位數(shù)據(jù)得到ki+1；

3) 將該序列進(jìn)行極值歸一化后帶入改進(jìn)閾值函數(shù)中，可使其具有小波包系數(shù)噪聲分布自適應(yīng)性。

2.3 閾值估計(jì)

(14)

圖2 閾值與的關(guān)系Fig.2 Correlation between threshold and

如圖2所示，按本方法計(jì)算得到的閾值為λ=0.6，為進(jìn)行對(duì)比，設(shè)定其他閾值分別為λ=0.5,0.7，利用文中自適應(yīng)改進(jìn)閾值函數(shù)(sym8小波，分解層數(shù)為3層)分別選定上述3個(gè)閾值進(jìn)行去噪分析，得到結(jié)果如圖3所示。由圖3可知，當(dāng)閾值λ=0.5時(shí)，閾值過(guò)小，信號(hào)毛刺較多，噪聲去除不完全；當(dāng)閾值取0.7時(shí)，原始信號(hào)被過(guò)度壓縮；而通過(guò)本方法確定的閾值取得了更好的去噪效果，證明了該閾值估計(jì)方法的有效性。

圖3 不同閾值去噪效果對(duì)比Fig.3 Denoising effect comparison of different threshold

為了驗(yàn)證該去噪方法受噪聲標(biāo)準(zhǔn)差的影響，同樣分析上述信號(hào)，分別加入不同標(biāo)準(zhǔn)差的混合噪聲，對(duì)其進(jìn)行去噪，計(jì)算去噪前后的信噪比(signal to noise ratio,簡(jiǎn)稱SNR)與均方根誤差(root mean square error,簡(jiǎn)稱RMSE)，且將結(jié)果列于表1。由表1可知，對(duì)不同噪聲標(biāo)準(zhǔn)差情況下的含噪信號(hào)進(jìn)行去噪，去噪后信號(hào)的SNR與RMSE變化幅度不大，證明該方法受噪聲標(biāo)準(zhǔn)差影響較小，適用于對(duì)不同噪聲尺度情況下的信號(hào)進(jìn)行去噪。

表1 不同噪聲尺度下的去噪結(jié)果

3 實(shí) 驗(yàn)

3.1 振動(dòng)實(shí)驗(yàn)裝置與信號(hào)獲取

為了對(duì)本算法進(jìn)一步驗(yàn)證，搭建了軸承振動(dòng)實(shí)驗(yàn)平臺(tái)如圖4所示，該實(shí)驗(yàn)平臺(tái)主要包括中間的小型高速直流電機(jī)，兩端測(cè)試軸承為分子泵中使用的微小型高速深溝球軸承QC0011286(軸承參數(shù)為：內(nèi)徑為4 mm，外徑為13 mm，軸承節(jié)圓直徑D為8.5 mm，滾珠直徑d為2.3 mm，滾珠個(gè)數(shù)N為7個(gè)，接觸角β為0°)以及固定于軸承外圈的加速度傳感器。測(cè)試時(shí)電機(jī)轉(zhuǎn)速設(shè)定為60 kr/min，采樣率為20 kHz，采樣時(shí)間為0.1 s，采樣時(shí)前置加上10 kHz的低通抗混濾波。為模擬軸承外圈故障，在軸承外圈內(nèi)壁加工寬為0.15 mm、深為0.2 mm的橫向溝槽，此時(shí)軸承基頻為fr=1 kHz，軸承外圈故障頻率根據(jù)公式(15)可得

(15)

采集到的軸承振動(dòng)時(shí)域波形如圖5所示。由圖5可知，從該時(shí)域信號(hào)無(wú)法直觀地得出任何軸承信號(hào)特征。

圖4 軸承振動(dòng)實(shí)驗(yàn)平臺(tái)Fig.4 Bearing vibration experiment platform

圖5 軸承振動(dòng)時(shí)域信號(hào)Fig.5 Bearing vibration time domain signal

3.2 去噪分析與對(duì)比

對(duì)圖5所示的信號(hào)選用sym8母小波進(jìn)行四層小波包分解，獲得最大尺度上的小波包系數(shù)序列w(i)，為進(jìn)行樣本熵計(jì)算，對(duì)其分割為若干個(gè)子序列，為在不失統(tǒng)計(jì)分布性的前提下盡可能準(zhǔn)確地表達(dá)小波包系數(shù)噪聲分布情況，進(jìn)行多次實(shí)驗(yàn)，最終取子序列長(zhǎng)度為l=127，求得小波包系數(shù)序列所對(duì)應(yīng)的樣本熵序列，且對(duì)其進(jìn)行歸一化即可得到閾值函數(shù)調(diào)節(jié)參數(shù)序列s如圖6所示。由圖6可知，該小波包系數(shù)序列中的噪聲分布不均勻，因此在對(duì)低頻與中頻小波包系數(shù)去噪時(shí)，調(diào)節(jié)參數(shù)較小，去噪形式偏硬，而對(duì)其他含噪較多(調(diào)節(jié)參數(shù)較大)的系數(shù)則去噪形式偏軟。

圖6 閾值函數(shù)調(diào)節(jié)參數(shù)Fig.6 Adjust parameters of threshold function

圖7 不同去噪算法結(jié)果對(duì)比Fig.7 Denoising results comparison of different algorithm

根據(jù)以上分析利用本方法對(duì)實(shí)驗(yàn)信號(hào)進(jìn)行去噪，為進(jìn)行直觀對(duì)比，將去噪前后信號(hào)進(jìn)行功率譜分析且與其他閾值去噪方法相對(duì)比得到如圖7所示的結(jié)果。由于軸承外圈故障而產(chǎn)生周期性脈沖信號(hào)，其表現(xiàn)形式為調(diào)制信號(hào)，即以軸承的轉(zhuǎn)動(dòng)頻率為調(diào)制頻率，外圈故障頻率為載波頻率，形成邊頻帶。通常由于故障特征微弱且存在噪聲干擾，軸承故障特征及其調(diào)制特征無(wú)法清晰顯露。通過(guò)信號(hào)去噪處理后，利用功率譜分析得到軸承振動(dòng)信號(hào)頻譜，且由于不同的去噪方法得到的功率譜分析效果不同，因此為了直觀地分析軸承振動(dòng)信號(hào)去噪效果，以去噪后信號(hào)特征頻率成分的還原情況作為去噪效果評(píng)判標(biāo)準(zhǔn)，具體分析如下。

由圖7(a)可知，原始含噪軸承振動(dòng)信號(hào)功率譜中僅能分辨出軸承基頻1 kHz及其倍頻2 kHz的頻率成分，其他有用的頻率成分被噪聲所淹沒(méi)；通過(guò)自相關(guān)去噪法去噪后得到的結(jié)果如圖7(b)所示，由圖可知，自相關(guān)去噪法去除了部分噪聲，還原了信號(hào)的3 kHz頻率成分，但故障頻率依然難以分辨，這表明該信號(hào)所含噪聲并非單純的高斯白噪聲，想要得到更好的去噪結(jié)果，有必要考慮其他方法；通過(guò)軟、硬閾值函數(shù)且選用基于Stein無(wú)偏似然估計(jì)的閾值確定規(guī)則(文獻(xiàn)[8]同樣采取該規(guī)則)，得到去噪結(jié)果如圖7(c)與(d)所示，硬閾值函數(shù)去噪后，能夠識(shí)別軸承故障特征頻率2 552.9 Hz，但在故障頻率周圍存在著無(wú)效頻率成分，無(wú)法判斷該故障特征的有效性，而軟閾值函數(shù)雖然濾除了大部分噪聲成分，但原始頻率成分如基頻1 kHz，故障特征頻率2 552.9 Hz也被扼殺嚴(yán)重；文獻(xiàn)[8]閾值函數(shù)去噪結(jié)果如圖7(e)所示，該結(jié)果較軟硬閾值函數(shù)更好，信號(hào)的基頻及其倍頻1,2,3 kHz與軸承外圈故障頻率foc=2 552.9 Hz及其調(diào)制成分3 552.9 Hz被很好的還原，然而在3 500 Hz與高頻區(qū)域仍然存在一些無(wú)關(guān)頻率成分，去噪效果有待提升；圖7(f)顯示自適應(yīng)閾值函數(shù)與閾值估計(jì)方法去除噪聲明顯且有效地還原了原始信號(hào)的頻率特征(基頻及其諧波1,2，3 kHz；故障頻率及其倍頻與調(diào)制成分2 552.9,3 552.9,5 150.8 Hz)，信號(hào)失真較少，去噪效果較好，提升了軸承故障診斷的準(zhǔn)確性。

事實(shí)上，通過(guò)對(duì)比信號(hào)去噪前后0.01 s的細(xì)致波形(如圖8所示)也可以直觀地發(fā)現(xiàn)，信號(hào)去噪前其波形毛刺較多，即存在較多高頻噪聲干擾，而去噪后信號(hào)細(xì)致波形毛刺較少，信號(hào)波形較為規(guī)律。

通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)，本算法雖然去噪效果優(yōu)于傳統(tǒng)去噪算法，但用時(shí)較長(zhǎng)，計(jì)算速度較傳統(tǒng)算法更慢，不適用于在線實(shí)時(shí)去噪分析與故障診斷。因此，在實(shí)際應(yīng)用中，可利用傳統(tǒng)的軟、硬閾值去噪算法作為在線實(shí)時(shí)初步分析工具，而本算法用于離線的進(jìn)一步精細(xì)分析。

圖8 去噪前后信號(hào)細(xì)致波形對(duì)比Fig.8 Detail waveform comparison of noisy and denoised signal

4 結(jié) 論

1) 樣本熵可以表征信號(hào)中不同種類的噪聲含量的大小，樣本熵越大，則信號(hào)含噪越多。將樣本熵應(yīng)用于小波包系數(shù)序列中，能得到其噪聲分布情況，據(jù)此對(duì)閾值函數(shù)進(jìn)行自適應(yīng)調(diào)整，使得其能夠?qū)胼^多的小波包系數(shù)進(jìn)行大尺度壓縮，而使含信號(hào)較多的小波包系數(shù)得到盡可能的保護(hù)。

3) 利用本方法對(duì)滾動(dòng)軸承振動(dòng)實(shí)驗(yàn)信號(hào)進(jìn)行分析，能夠獲得較其他方法更好的去噪效果，且有效地還原了信號(hào)的轉(zhuǎn)動(dòng)特征頻率與故障特征頻率，提高了軸承狀態(tài)監(jiān)測(cè)的準(zhǔn)確度。

4) 本方法去噪效果較傳統(tǒng)算法更好，但計(jì)算速度較慢，實(shí)際應(yīng)用中可將傳統(tǒng)算法作為在線初步去噪算法，本方法作為離線精細(xì)算法進(jìn)行配合分析。