基于樣本熵的改進小波包閾值去噪算法*

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

(西南科技大學制造過程測試技術教育部重點實驗室 綿陽，621000)

引 言

長期以來，小波包濾波是對信號進行分析預處理的主要工具，該方法對信號進行小波包變換，然后利用小波包閾值函數對小波包系數進行閾值收縮處理以達到去噪目的。然而由于含噪信號中的噪聲分布往往不均勻，傳統的軟、硬閾值方法對信號小波包系數進行固定格式的閾值處理，無法很好地滿足信號去噪要求[1]。另外，在閾值的選擇方面，常用的Heursure閾值、Donoho閾值等能夠對信噪比較高的信號實現噪聲與信號的最優分離，而對于強噪聲信號，去噪效果并不理想[2]。近年來，針對這兩個問題，許多學者進行了研究。

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

基于以上分析，筆者提出將樣本熵與小波包閾值去噪算法相結合，且將其應用于高速深溝球軸承振動信號去噪分析，通過分析信號的小波包系數噪聲分布情況及其對應的樣本熵值，使閾值去噪算子具有噪聲分布自適應性以達到最優去噪效果；同時以噪聲估計信號樣本熵值為基準，提出了一種最優閾值估計方法。仿真分析以及實驗結果皆對所提方法進行了驗證。

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

1.1 樣本熵的計算

設有長度為N的時間序列Xi={x1,x2,…,xN}，其樣本熵的計算方法如下：

1) 確定嵌入維數為m，對Xi的元素按順序進行排列，即可得到一組維數為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)]為兩向量之間對應元素求差的絕對值的最大值，即

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

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

(2)

3) 對于固定的Xm(i)，統計Xm(i)與Xm(j)之間距離小于等于參數r的j(1≤j≤N-m,j≠i)的個數，且記為Bi，則當1≤i≤N-m時定義

(3)

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

(4)

(5)

6) 定義Am(r)為

(6)

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

(7)

實際信號中N無法趨近于無窮，因此可將樣本熵設為

(8)

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

1.2 樣本熵與噪聲

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

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

(9)

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

圖1 樣本熵與幾種常見噪聲大小關系Fig.1 Correlation between sample entropy and several common noise

由圖1(a)看出，信號的樣本熵值與噪聲標準差成正相關，且在采樣點相差巨大的情況下其樣本熵值變化依然相近。圖1(b)，(c)較圖1(a)而言，樣本熵隨噪聲標準差變化有少許波動，熵值大小區間有變化，但總體趨勢仍呈正相關，且數據長度對其趨勢影響不大。這說明，雖然對仿真信號施加不同種類的噪聲，其樣本熵變化曲線有所不同，但由于對信號增加噪聲的結果導致了信號的隨機度與復雜度增加，其結果仍然是樣本熵增加。圖1所示說明了樣本熵算法可用來表征時間序列的含(多種)噪聲情況且受數據長度影響較小。

2 改進小波包閾值去噪算法

2.1 傳統小波包去噪及其不足

傳統的閾值函數主要有軟閾值函數和硬閾值函數兩種。軟閾值函數

(10)

硬閾值函數

(11)

工程實踐中常用的閾值函數還有形如文獻[8]中提出的一種介于軟、硬閾值函數之間的改進閾值函數

(12)

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

2.2 自適應改進閾值函數

根據以上分析，筆者對傳統閾值函數作改進，使其不受限于固定去噪形式，得到的改進閾值函數為

(13)

其中：0

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

由1.2節分析可知，樣本熵能較好地反應時間序列的噪聲變化情況，且適用于分析短時間序列，因此設定調節參數s的確定方法如下：

1) 對信號的小波包系數序列(w1,w2,…,wn)按順序分割為n-l個相互之間最大重疊且長度均為l的子序列(k1,k2,…,kn-l)，即ki向后移動一位數據得到ki+1；

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

2.3 閾值估計

(14)

圖2 閾值與的關系Fig.2 Correlation between threshold and

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

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

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

表1 不同噪聲尺度下的去噪結果

3 實 驗

3.1 振動實驗裝置與信號獲取

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

(15)

采集到的軸承振動時域波形如圖5所示。由圖5可知，從該時域信號無法直觀地得出任何軸承信號特征。

圖4 軸承振動實驗平臺Fig.4 Bearing vibration experiment platform

圖5 軸承振動時域信號Fig.5 Bearing vibration time domain signal

3.2 去噪分析與對比

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

圖6 閾值函數調節參數Fig.6 Adjust parameters of threshold function

圖7 不同去噪算法結果對比Fig.7 Denoising results comparison of different algorithm

根據以上分析利用本方法對實驗信號進行去噪，為進行直觀對比，將去噪前后信號進行功率譜分析且與其他閾值去噪方法相對比得到如圖7所示的結果。由于軸承外圈故障而產生周期性脈沖信號，其表現形式為調制信號，即以軸承的轉動頻率為調制頻率，外圈故障頻率為載波頻率，形成邊頻帶。通常由于故障特征微弱且存在噪聲干擾，軸承故障特征及其調制特征無法清晰顯露。通過信號去噪處理后，利用功率譜分析得到軸承振動信號頻譜，且由于不同的去噪方法得到的功率譜分析效果不同，因此為了直觀地分析軸承振動信號去噪效果，以去噪后信號特征頻率成分的還原情況作為去噪效果評判標準，具體分析如下。

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

事實上，通過對比信號去噪前后0.01 s的細致波形(如圖8所示)也可以直觀地發現，信號去噪前其波形毛刺較多，即存在較多高頻噪聲干擾，而去噪后信號細致波形毛刺較少，信號波形較為規律。

通過實驗發現，本算法雖然去噪效果優于傳統去噪算法，但用時較長，計算速度較傳統算法更慢，不適用于在線實時去噪分析與故障診斷。因此，在實際應用中，可利用傳統的軟、硬閾值去噪算法作為在線實時初步分析工具，而本算法用于離線的進一步精細分析。

圖8 去噪前后信號細致波形對比Fig.8 Detail waveform comparison of noisy and denoised signal

4 結 論

1) 樣本熵可以表征信號中不同種類的噪聲含量的大小，樣本熵越大，則信號含噪越多。將樣本熵應用于小波包系數序列中，能得到其噪聲分布情況，據此對閾值函數進行自適應調整，使得其能夠對含噪較多的小波包系數進行大尺度壓縮，而使含信號較多的小波包系數得到盡可能的保護。

3) 利用本方法對滾動軸承振動實驗信號進行分析，能夠獲得較其他方法更好的去噪效果，且有效地還原了信號的轉動特征頻率與故障特征頻率，提高了軸承狀態監測的準確度。

4) 本方法去噪效果較傳統算法更好，但計算速度較慢，實際應用中可將傳統算法作為在線初步去噪算法，本方法作為離線精細算法進行配合分析。