改進樣本熵最優小波包閾值選擇算法在信號降噪中的應用

何 勇，張祥金，姚宗辰

(南京理工大學 機械工程學院， 南京 210094)

在信號降噪處理中小波包閾值的選擇顯得尤為重要，傳統的閾值估計方法如軟閾值、硬閾值等在信噪比較高時能得到信號與噪聲的最優分離，而對于受噪聲污染嚴重的信號，其去噪效果往往不甚理想[1]。羅元[2]提出基于Teager能量算子的改進閾值函數去噪算法；康玉梅[3]提出了小波包分數冪閾值的新方法；劉沖沖[4]提出了一種基于自適應閾值和新閾值函數的小波包語音增強算法；周建[5]利用樣本熵作為信息價值函數以確定最優小波包，且以樣本熵為判據，對不同的分解層數設置不同的閾值，選取使得去噪后得到的噪聲信號樣本熵值最大的閾值作為最優閾值，并且與傳統閾值選擇進行了比較分析。

文獻[5]給出了小波包樣本熵最優閾值選擇原理，但是沒有考慮閾值步長大小對樣本熵曲線的影響，也沒有詳細說明怎樣設置閾值步長。本文在此基礎上，分析了初始閾值以及閾值增加步長的大小對噪聲序列樣本熵產生的影響。為了改進算法，提出將閾值作為變量，閾值步長作為新的環境特征，構造目標函數，結合粒子群動態優化算法估計小波包最優閾值。

1 相關理論

1.1 小波包理論

小波包分解是在全頻帶內對信號進行多分辨率分解系數的分析方法。設采樣到的離散數據為{X(k)}k=1,2,…,N，對其進行小波包的分解，分解算法[6]為

(1)

(2)

1.2 樣本熵理論

樣本熵反映時間序列的復雜度與無規律程度，樣本熵值越小，則時間序列的自相似度越高，序列越規則; 熵值越大，表示序列越為復雜，取值也更加隨機[7]。

設有時間序列為Xi={x1,x2,…,xN},長度為N,則樣本熵計算如下。

① 設定嵌入維數為m和相似容限r,考慮維數為m的向量組{xm(1),…,xm(N-m-1)},其中：

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

1≤i≤N-m+1

(3)

② 設兩個向量Xm(i)和Xm(j)之間距離d[Xm(i),Xm(j)]為向量之間對應元素之差絕對值的最大值。則有：

(4)

③ 設定固定的Xm(i),對Xm(i)與Xm(j)之間的距離小于等于相似容限r的j(1≤i≤N-m，j≠i)的數目進行統計，并記作Bi。則當1≤i≤N-m時，定義：

④ 設B(m)(r)為：

(5)

⑤ 將嵌入維數增加到m+1，按照上述計算Xm+1(i)與之間距離小于等于r的個數，并記作Ai。則有：

(6)

⑥ 同樣地

(7)

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

(8)

實際信號中N通常為有限值，樣本熵可以估計為:

(9)

一般r=0.1Std～0.25Std(Std為時間序列Xi={x1，x2，…，xN} 的標準差),取m=1或2。在本文的分析中，取m=2，r=0.2Std。

1.3 動態粒子群優化算法

粒子群算法，也稱粒子群優化算法或鳥群覓食算法(PSO)。PSO算法從隨機解出發，通過追隨當前搜索到的最優值來尋找全局最優，以適應度來評價解的品質。

系統初始化為一組隨機解，通過迭代搜尋最優值，粒子在解空間追隨最優的粒子進行搜索。粒子速度公式為[8]：

vid(t+1)=ω*vid(t)+c1r1(pid-xid(t))+

c2r2(pgd-xid(t))

(10)

位置更新公式：

xid(t+1)=xid(t)+vid(t+1)

(11)

式(10)～式(11)中：vid(t)，xid(t)分別表示t時刻粒子i的飛行速度和位置；pid表示粒子i的個體歷史最佳位置；pgd表示群體中最佳位置；ω表示慣性權重；c1和c2表示加速因子；r1和r2表示在[0,1]范圍內的隨機數。

實際問題中，有些目標函數和約束條件也會隨著變量的變化而變化，最終導致其問題的最優解改變，是一個動態優化問題，采用改進的粒子群優化算法。種群不隨算法迭代而變化，而目標函數和約束函數會隨變量改變而改變，那么求解的問題也改變了。因此引入環境檢測算子，檢測環境變化計算公式為[9]：

(12)

式(12)中，n為重新評價個體的數目，一般為種群大小的10%；fj(x,t-1)為環境變化前的適應值；fj(x,t) 為環境變化后的適應值。對計算結果ε(t)進行如下分類

(13)

式(13)中,ε1為劇烈強度變化的分界點；ε2為中等程度分界點，通常10-3<ε2≤10-2。若ε(t)>ε1，為劇烈的環境變化，此時環境變化強度較大，最優值及其位置將會偏離原始值；當ε2<ε(t)≤ε1，認為是中等強度的環境變化，最優解及其位置應該在原始值附近波動。

2 粒子群動態優化算法的樣本熵最優小波包閾值選擇原理

目前，已經獲取了某炮彈脈沖激光引信發射過載信號。如圖1所示為回收的原始信號。

圖1 加速度原始信號時域頻譜圖

由圖1可知，盡管回收系統硬件濾波電路濾除大于5 kHz的信號，但是無法完全濾除內部電子噪聲。以該回收的原始過載信號為研究對象，利用改進的小波包閾值選擇算法分析降噪過程。

2.1 樣本熵最優閾值選擇原理

設含噪的原始信號為yn(t),小波包設定閾值重構信號為y(t),消除的噪聲信號為n(t)，即有n(t)=yn(t)-y(t)。若小波包閾值選擇較小時，重構的信號y(t)中必然含有噪聲，而噪聲n(t)此時序列較為簡單，熵值較低；隨著閾值增高，噪聲序列變得復雜，熵值增高，當與之增大到n(t)恰好包含所有噪聲時，熵值最大；再隨著閾值增高，噪聲信號n(t)中必包含有效信號(有效信號序列較規則)，熵值會開始降低；當閾值增加到一定值后，此時噪聲序列同時包括了噪聲和有用信號，此后無論閾值如何增加，樣本熵曲線都不會再改變，呈水平線。因此，噪聲序列n(t)熵值大體趨勢是先增加后減少，一定閾值后，呈水平線。取熵值最大的對應的閾值作為小波包重構最優的閾值。

文獻[5]僅從原理上分析樣本熵曲線而未考慮閾值步長設置的影響以及噪聲序列樣本熵曲線在局部的隨機性。若按照文獻[5]進行閾值選擇，初始閾值設定為0，閾值區間為[0,4]。設置不同的閾值步長step，計算不同閾值所對應的噪聲序列樣本熵。如圖2所示。

整體來看，無論閾值步長如何變化，噪聲序列的樣本熵曲線總體趨勢都是先增加后減小。而當閾值增加的步長變化時，最大熵值對應的閾值也在不斷變化；閾值步長大，則噪聲序列樣本熵對應的最優閾值誤差就大。圖2中閾值步長設置不同時，最優閾值也不同(0.4、0.1或者0.15)，因此無法快速獲得最優閾值。

圖2 不同步長對應噪聲序列樣本熵曲線

局部來看，由于噪聲的隨機性，噪聲序列局部樣本熵大小排列也變得具有隨機性，如圖3所示，閾值范圍為[0,0.5]的噪聲序列樣本熵曲線。若閾值步長較小(如0.0005)，雖然樣本熵曲線能夠在整體上呈現先增后減的趨勢，但是局部走勢會趨于平緩，有可能會出現多個峰值,甚至由于噪聲的隨機性，會過早的陷入局部最優解,而且運算代價高。

圖3 不同步長噪聲序列局部樣本熵曲線

由圖2、圖3可以看出，閾值步長的調整是局部與整體的組合關系，本質上是噪聲序列的隨機性與樣本熵規律性的反應。無論步長設置多大，都無法準確快速地找出最大熵值所對應的閾值，且不同的閾值步長，最優小波包閾值也會有所不同。

因此，在樣本熵最優小波包閾值估計的理論基礎上，提出結合粒子群動態優化算法來確定最優小波包閾值。

2.2 基于粒子群動態算法閾值優化模型

假設噪聲序列閾值區間為[ai,bi]，設定閾值步長為step，則閾值序列thr=[ai,…,ai+k*step,…,bi],即有：

thr=g(k;step),k=1,2,3，…,n

(14)

設閾值序列對應重構信號序列y(t)=[yai(t),…,yai+k*step(t),…,ybi(t)],即有：

y(t)=h(thr,t)

(15)

若yn(t)原始信號，則thr對應的噪聲序列即有：

n(t)=yn(t)-y(t)

(16)

噪聲序列對應的樣本熵曲線為

ysamp=Q(n(t))

(17)

而由前面分析，噪聲序列樣本熵曲線是隨著閾值步長的變化而變化的，即目標函數ysamp=Q(n(t))是隨著step的變化而變化，目標函數和約束條件也會隨著步長的變化而變化，最終導致其問題的最優解改變，當step→0時(step≠0，可以趨于0，否則閾值就不會改變，就是固定閾值)，樣本熵序列就會變為樣本熵曲線。因此，這是一個動態優化問題，需要增加環境變化因子，該模型中step作為環境變化的因子。一般step初始設置較大，不斷減小，以便有較快的收斂速度。

設step=[1,…,0.4,…,0.04,…,0.004,…]是隨時間t變化，則有：

step=f(s;t),s=1,2,3,…,n

(18)

結合式(9)、式(14)～式(18),此時，求解最優閾值問題轉化為：

max:ysamp=Q(yn(t)-h(g(k;f(s;t)),t))

s.t.ai≤ai+k*step≤bi;

k=1,2,3,…,n；

step=f(t;s)且step≠0；

s=1,2,3,…,n；

其中，ysamp為目標函數,aibi則分別是閾值thr的上下界。

采用改進粒子群算法，在算法迭代過程中根據不同的環境變化采取不同的種群多樣性方法。如圖4所示，其算法如下：

圖4 粒子群動態優化算法

① 建模；將問題空間映射到算法空間，設定相關參數從c1、c2、ω以及粒子數量。

② 初始化種群和最優解存儲空間，隨機生成群體中粒子候選解的位置和速度。

③ 計算每個粒子對應的適應度值，存儲當前粒子最優解，檢測粒子的適應度值。

④ 根據式(12)計算step改變前后粒子適應度值的變化ε(t)。比較ε(t)、ε1、ε2；根據實際情況，本文取ε1=0.1、ε2=0.001。

⑤ 若ε(t)>ε1，返回步驟②；若ε2<ε(t)≤ε1，則當前最優解的粒子位置隨機生成部分種群，其余粒子重新初始化，返回到③；若ε(t)≤ε2，則更新所有尋優粒子的位置和速度。

⑥ 判斷迭代是否完成；若沒有完成返回②；若完成則計算結束。

2.3 改進樣本熵最優小波包閾值去噪算法原理

由以上分析，對改進樣本熵最優小波包閾值的去噪算法如圖5所示。

圖5 樣本熵去噪算法原理

原理如下：

① 對原始信號選擇合適的小波包分解層數和最優小波基。

② 對小波包分解的最低層開始，給一個較大閾值步長，計算不同閾值對應的噪聲序列樣本熵，確定最底層小波包閾值的大致范圍thr∈[ai,bi]。

③ 將閾值作為變量，步長作為自變量，構造閾值函數thr=g(k;step)，由于步長是環境變化量，構造步長step隨時間變化的函數step=f(s;t)(s=1,2,…)。

④ 建立最優小波包閾值數學模型，將問題空間映射到算法空間，利用粒子群動態優化算法不斷尋找全局最佳閾值(噪聲序列樣本熵最大對應的閾值)。

⑤ 依據②、③、④確定每層節點系數的最佳閾值。

⑥ 每層系數中小于該層最佳閾值的小波包系數置零，大于該層最佳閾值的小波包系數減去閾值后保留下來。

⑦ 對選擇的系數進行重構得到去噪后的信號。

3 數據分析與誤差評價

為了更好地說明閾值步長設置的影響以及結合粒子群動態優化算法的優越性，同時為了簡化計算，對該信號小波包分解的母樹即原始信號進行動態優化分析。在算法運行時，閾值步長step=f(s;t)隨著一定時間變化，獲得的最佳閾值和最佳樣本熵也在不斷變化。

圖6表示分別進行三次算法尋優，步長變化時，最佳閾值的變化軌跡。圖6中①為步長開始變化時的最佳閾值，隨著步長變化，最佳閾值也開始產生變化，在圖6中②出現之前，環境劇烈程度一直為ε(t)>ε1；圖6中②表示當步長變為某一值時，此時ε2<ε(t)≤ε1，出現劇烈強度變化的分界點，且從②之后，每次計算都會有不一樣的最佳閾值，這是由于隨著步長的減小，噪聲序列樣本熵曲線在局部會顯示多個波峰，因此會產生一定的誤差；圖6中③表示算法運行最終的最優閾值為0.135，在閾值誤差≤0.01時，認為在可接受誤差范圍。沒有出現ε(t)≤ε2的環境變化情況。

圖6 三次算法運行的最優閾值位置

為了進行誤差評價與精度分析，在算法計算過程中選取閾值步長為0.4、0.04、0.004的某一次計算的歷史最優值進行評估,如圖7、8所示。

圖7 歷史最佳適應值(樣本熵)收斂曲線

圖8 歷史最佳位置(閾值)收斂曲線

由圖7和圖8可知，當step=0.4時，開始迭代時就獲得了最佳位置和適應值(0.4，1.99)；當step=0.04時，曲線經過14次迭代到達最優值(0.16，2.3)；當step=0.004時，經過88次迭代到達最優值(0.143，3.37)，與算法運行最終的最優閾值(0.135，2.39)相比較，在閾值誤差允許范圍內(≤0.01)，也可認為此次尋優獲得的最佳閾值為0.143。

在未改進算法前，將閾值步長設置為固定值，可能選取固定的閾值步長為0.4、0.04等，獲得的最優閾值也可能是0.4、0.16等。而改進算法之后，步長作為隨時間的改變量，在誤差范圍內可認為算法最終獲得的最優閾值必然為0.143。

評價算法精度，最終是看小波包閾值降噪的效果。去噪方法的優劣都要用去噪的質量來評價，而信噪比(SNR)和均方根誤差(RMSE)是被廣泛采用的評價標準,信噪比越大且均方根誤差越小，則去噪效果越好。

(19)

(20)

圖9 不同閾值重構信號

列出各個重構信號與原始信號的信噪比和均方根誤差，如表1所示。

表1 去噪效果對比

對比圖1和圖9，結合表1可以看出，在未改進算法前，獲得的最優閾值可能是0.4、0.16等。由圖9可以看出，閾值為0.4和0.16的重構信號雖然去除了噪聲，但是由于閾值設置偏大，同樣也去除了部分有用信號。而改進算法之后，獲得的最優閾值為0.143，對應重構的信號，既去除了噪聲又很好地保留了有用信號，雖然有局部噪聲，但是在誤差允許范圍內可以忽略。且通過表1中的比較，閾值為0.143的重構信號效果好。

4 結論

結合粒子群動態優化算法，將閾值步長設為隨時間變化的因變量，利用樣本熵曲線來確定最優小波包閾值。以激光引信發射過載信號為例，對比分析了人為選擇固定閾值步長而可能獲得的最優閾值和算法選擇的最優閾值的去噪效果。結合粒子群動態優化算法改進了樣本熵尋找最優小波包閾值的缺陷，避免了陷入局部最優解以及人為選擇步長誤差的影響，獲得的最佳閾值去噪效果更好。