基于模糊遞歸和最優(yōu)硬閾值的局部投影降噪算法

王 東, 崔忠馬, 陳文東, 舒 勤,*

(1. 四川大學電氣工程學院, 四川 成都 610065; 2. 北京遙感設備研究所, 北京 100084)

0 引 言

實際測得的信號易受噪聲干擾,從而破壞信號的結(jié)構(gòu)特征,因此降噪是后續(xù)各種信號處理的基礎[1-2]。傳統(tǒng)信號降噪方法基于線性系統(tǒng)的假設,認為信號和噪聲的頻譜不完全重合。但是，由混沌系統(tǒng)產(chǎn)生的信號往往具有寬頻譜、內(nèi)在偽隨機等混沌特性,信號和噪聲頻譜重合,使得傳統(tǒng)線性濾波方法失效[3]。

Cawley等[4]和Sauer等[5]提出的基于混沌動力學理論的局部投影降噪算法為具有寬頻譜特性的混沌信號降噪提供了一種新方法。該方法首先根據(jù)Takens定理重構(gòu)信號,得到一個與原始動力系統(tǒng)微分同胚的相空間[6]。期望信號的吸引子被限制在一個低維流形上,而噪聲的吸引子分散在流形周圍。局部投影降噪算法根據(jù)信號與噪聲在相空間中局部動力學特性與局部幾何特性的不同，區(qū)分信號子空間和噪聲子空間,利用幾何投影去除噪聲分量,再將相空間反重構(gòu)為時間信號,從而達到降噪的目的。在實際情況中，許多信號都有混沌特性,目前該算法成功應用的場景包括人類語音信號處理[7]、振動信號分析[8]，還包括生物信號處理,如腦電信號處理[9]、心電信號處理[10]、呼吸聲音信號處理[11],以及激光數(shù)據(jù)處理[12]、故障檢測[13]等。

鄰域選取和子空間劃分是局部投影算法的兩個研究重點。鄰域選取的方法很多,但都存在不足。Cawley等[4]和Sauer等[5]直接指定鄰點數(shù)目,這種方式易受人為因素影響。Kantz等[12]利用遞歸圖估計鄰域半徑,效果優(yōu)于原始方法,但計算復雜。馮飛龍等[14]利用小波分解估計初始鄰域半徑,再進行鄰點搜索,徐禮勝等[15]采用經(jīng)驗模態(tài)分解(empirical mode decomposition, EMD)估計鄰域半徑[15],但這些方法都需要預先估計噪聲水平。Przybya等[16]和Kotas等[17]利用K-means聚類算法確定鄰域,但聚類數(shù)難以確定。對于子空間劃分,現(xiàn)有方法直接指定子空間的維數(shù),或根據(jù)特征值的大小對子空間進行劃分。Chelidze等[18]通過重構(gòu)信號的短時軌跡,利用平滑正交分解識別子空間。但考慮到混沌吸引子本身具有分數(shù)維的性質(zhì)[19],系統(tǒng)局部鄰域的動力學特性以及每個鄰域內(nèi)的噪聲分量占據(jù)的空間也不盡相同,且在實際情況中可能并不知道原始動力系統(tǒng)的相空間維數(shù),所以每個局部鄰域應該進行不同的子空間劃分。

針對以上問題,本文提出基于模糊遞歸圖與最優(yōu)硬閾值準則的局部投影降噪算法。首先，根據(jù)模糊遞歸圖對鄰域進行選擇。為避免計算鄰域協(xié)方差矩陣,直接將鄰域矩陣進行奇異值(singular value decomposition, SVD)分解；然后，根據(jù)最優(yōu)硬閾值對局部鄰域的信號子空間和噪聲子空間進行劃分,避免人為因素的影響；最后，針對高斯白噪聲,采用本文所提方法分別對仿真Lorenz信號與實測含噪心電圖(electrocardiogram, ECG)信號進行仿真研究,并與其他局部投影算法以及其他ECG信號降噪方法進行對比,仿真結(jié)果驗證了本文方法的有效性。

1 局部投影算法原理

1.1 相空間重構(gòu)

根據(jù)Takens嵌入定理,對于無限長、無噪的d維混沌吸引子的標量時間信號{x(t)},在拓撲不變的意義下可以找到一個m維的嵌入相空間。其中,m≥2d+1。而對于有限長、含噪的信號，可采用坐標延遲重構(gòu)[20]。

設x1,x2,…,xN為一長度為N的含高斯白噪聲的單變量時間信號,選定一個時間延遲τ和嵌入維數(shù)m,構(gòu)造如下的相空間矢量:

Xi=[xi,xi+τ,…,xi+(m-1)τ]T,i=1,2,…,M

(1)

式中:M=N-(m-1)τ。τ由平均互信息量法[21]確定,即

(2)

式中：yi=xi+τ。取I(τ)第一個極小值點對應的時間作為時間延遲。m由Cao[22]所提的方法確定,即

(3)

1.2 局部投影算法

基于相空間重構(gòu)局部投影降噪算法,利用信號與噪聲在相空間中軌線的動力學特性與幾何特性的不同,保留信號分量,抑制噪聲分量,最大程度地恢復原始信號的吸引子流形。

對于相點Xn,將其在鄰域內(nèi)線性化展開:

(4)

(5)

(6)

第n個相點Xn的鄰域加權(quán)矩陣為

(7)

對式(7)進行SVD分解,得到:

(8)

式中:Vn=[a1,a2,…,aK]由信號子空間與噪聲子空間構(gòu)成。將SVD分解得到的奇異值按從小到大進行排列，有σ1≥σ2≥…≥σK-1≥σK,其中大的奇異值對應信號子空間,小的奇異值對應噪聲子空間。根據(jù)Vn得到由噪聲引起的分量,減去第n個相點在噪聲子空間中的投影,即得到修正后的相點:

(9)

為了進一步提高降噪效果,避免局部線性化產(chǎn)生較大誤差,本文采用Moore等[23]提出的質(zhì)心修正方法。修正后的質(zhì)心為

(10)

2 改進的局部投影算法

2.1 鄰域選擇

在局部投影降噪算法中,鄰域的選擇十分重要。鄰域選擇得過小,會損失有效信息,受噪聲干擾嚴重;鄰域選擇得過大,會使得線性逼近的效果不好。

遞歸圖[24]可以用來提取時間信號中的相關信息。但是，遞歸圖的缺點是,動力系統(tǒng)的遞歸模式可視化對相似閾值的選取十分敏感[25]。Pham等[26]在遞歸圖基礎上提出了模糊遞歸圖。相比于傳統(tǒng)遞歸圖,模糊遞歸圖不需要選取相似閾值,且模糊遞歸圖以灰度圖的形式展示,能夠為模式分析提供更豐富的信息。本文采用模糊遞歸圖來對鄰域進行選擇。

(11)

式中：μ(Xi,Xj)∈[0,1]表示Xi和Xj之間的一種模糊相似性的度量,其具有如下性質(zhì)。

(1) 自反性:

μ(Xi,Xi)=1

(12)

(2) 對稱性:

μ(Xi,Vj)=μ(Vj,Xi)

(13)

(3) 傳遞性:

μ(Xi,Xj)=max[min{μ(Xi,Vk),μ(Xj,Vk)}]

(14)

式中：i=1,2,…,M;k=1,2,…，c；c為聚類數(shù),1

μ(Xi,Xj)的值根據(jù)模糊C均值(fuzzy C-Mean, FCM)聚類算法[27]計算得到。FCM算法通過最小化如下模糊目標函數(shù)實現(xiàn):

(15)

式中：ω∈[1,+∞)為模糊度參數(shù);U=[μij](i=1,2,…,M;j=1,2,…,c)是劃分矩陣;V=(V1,V2,…,Vc)是聚類中心矩陣;Vj表示第j個聚類中心;d(Xi,Vj)表示某一范數(shù),本文采用歐式范數(shù)。上述模糊目標函數(shù)滿足

(16)

為了得到最優(yōu)的U和V,通過迭代過程數(shù)值求解目標函數(shù)的最小值,迭代過程如下:

(17)

(18)

當滿足‖Ut-Ut+1‖<ε時,停止迭代。其中,t為迭代次數(shù),ε為給定的精度水平。

模糊遞歸圖具有對稱性,可以看作是相空間狀態(tài)矢量之間的一種模糊關系。模糊遞歸圖用灰度圖表示,灰度值代表了狀態(tài)矢量對之間的模糊關系,這與遞歸圖是相互兼容的?；叶戎翟叫?則表示這兩個狀態(tài)矢量越相似。一個灰度值為0的像素點代表了兩個狀態(tài)矢量完全相似,即代表動力系統(tǒng)中一個100%的遞歸事件[28]。

2.2 信號與噪聲子空間劃分

在進行投影修正之前,需要根據(jù)SVD分解得到的奇異值對信號子空間與噪聲子空間進行劃分。本文利用Gavish等[29]提出的最優(yōu)硬閾值準則對信號子空間與噪聲子空間進行劃分,不需要預先估計噪聲水平。最優(yōu)硬閾值γ為

γ=ω(β)σmed

(19)

式中:σmed為奇異值的中位數(shù);ω(β)=λ(β)/μβ,λ(β)為

(20)

對于前述的鄰域加權(quán)矩陣Zn∈Rm×n;當m=n時,β=1;當m>n時,β=n/m;當m

(21)

可以通過數(shù)值方法求解得到式(21)積分方程中μβ的值,從而求得閾值γ。

根據(jù)最優(yōu)硬閾值準則,將大于等于閾值γ的奇異值所對應的奇異向量所形成的子空間作為信號子空間,小于閾值γ的奇異值所對應的奇異向量所張成的子空間作為噪聲子空間,由此實現(xiàn)了子空間的劃分。

2.3 改進算法的基本步驟

基于以上分析,本文所提的局部投影降噪算法步驟如算法1所示。

算法1 改進的局部投影降噪算法步驟輸入:含噪信號x(t)輸出:降噪信號x^(t)開始 相空間重構(gòu):利用式(2)計算平均互信息量并取第一個局部極小值對應的時間作為時間延遲τ;利用式(3)計算E1(m),取停止變化時對應的維數(shù)作為嵌入維數(shù)m;對x(t)進行相空間重構(gòu),得到式(1)表示的相空間X。 循環(huán)(i=1∶M) 1. 選定參考相點Xi。 2. 鄰域選擇:根據(jù)式(11)計算重構(gòu)相空間的模糊遞歸圖,得到第i個參考相點的鄰域。 3. 計算鄰域質(zhì)心和鄰域矩陣:由式(10)和式(7)分別計算第i個參考相點的鄰域質(zhì)心和鄰域矩陣。 4. SVD分解:根據(jù)式(8)對鄰域矩陣進行SVD分解,得到奇異值與右奇異向量。 5. 子空間劃分:根據(jù)式(19)的最優(yōu)硬閾值準則對信號子空間與噪聲子空間進行劃分。 6. 投影修正:根據(jù)式(9)進行投影修正。 結(jié)束 反重構(gòu):將相空間恢復為時間信號x^(t),恢復方式采用式(22)[30],以減小誤差。結(jié)束

為減小反重構(gòu)所產(chǎn)生的誤差，進行如下操作：

(22)

為達到較好的效果，需將以上步驟重復迭代幾次。

3 仿真實例與工程應用

3.1 Lorenz混沌信號仿真

本文首先采用Lorenz混沌系統(tǒng)信號進行仿真。Lorenz混沌系統(tǒng)由如下偏微分方程組描述：

(23)

當參數(shù)取σ=10, r=28, b=8/3時,系統(tǒng)呈現(xiàn)出混沌特性。采用四階Runge-Kutta方法數(shù)值求解上述偏微分方程組,初始值取x0=10, y0=10, z0=10,積分步長為0.01。取X分量中的2 000個數(shù)據(jù)點進行仿真。對信號添加高斯白噪聲,添加噪聲后的信噪比為8 dB,然后使用本文所提的局部投影方法進行降噪處理。降噪效果如圖1所示。

由圖1可知,本文方法具有較好的降噪效果,在噪聲較強的情況下也能夠恢復信號。由圖1可以看出，降噪后，信號比較光滑,看不到噪聲的影響,同時基本保持了Lorenz混沌吸引子的幾何形狀。

3.2 參數(shù)對算法的影響

本小節(jié)討論鄰域選擇和閾值這兩個參數(shù)對算法的影響。首先是鄰域選擇對算法的影響。Lorenz信號添加噪聲后的信噪比為8 dB。選取不同鄰域?qū)π盘栠M行降噪,每個鄰域仿真200次，然后取輸出信噪比的平均值作為最終輸出信噪比。由圖2可知,當鄰域選取合適時,降噪后信噪比達到最大。鄰域選擇過小或過大都會使得信噪比下降。這是因為鄰域選擇過小,信號受噪聲影響明顯;而鄰域選擇過大,對于相空間軌線的逐段線性逼近效果差。

圖1 Lorenz混沌信號降噪前后時域波形圖和相空間圖Fig.1 Time-domain waveform and phase space diagram of Lorenz chaotic signal before and after denoising

圖2 鄰域選擇對算法的影響Fig.2 Effect on the algorithm of neighborhood selection

仍然選取上述信號,在計算得到的閾值上疊加一個隨機誤差后再進行降噪處理。由圖3可以看出,在有誤差和無誤差時，其降噪效果基本相同,這表明算法具有較好的魯棒性。

圖3 閾值誤差對算法的影響Fig.3 Effect of threshold error on the algorithm

3.3 不同局部投影方法對比

為了進一步研究本文所提方法的降噪效果,將其他局部投影降噪算法與本文所提方法進行對比。選取降噪后的信噪比、均方誤差、復雜度和耗時作為衡量降噪效果的評價指標。信噪比的計算公式為

(24)

(25)

首先,在信噪比為8 dB的情況下,計算混沌信號重要的特征量:李雅普諾夫指數(shù)、關聯(lián)維數(shù),然后再統(tǒng)計每個方法計算的耗時。每種方法設定相同參數(shù),迭代8次,仿真200次,然后計算平均耗時。由表1可以看出,由本文方法降噪后的信號計算得到的李雅普諾夫指數(shù)和關聯(lián)維數(shù)最接近原始Lorenz信號的對應值,這表明Lorenz混沌吸引子中的確定性結(jié)構(gòu)得到了較好的保留。另外,由表1可以看出，本文方法較為耗時，這是因為本文方法在計算過程中涉及了模糊遞歸圖的求解以及數(shù)值求解積分方程。但是本文方法在犧牲計算時間的情況下,具有更好的降噪效果。

表1 不同方法降噪后信號的李雅普諾夫指數(shù)、關聯(lián)維數(shù)以及計算耗時(SNR=8 dB)

同樣選定Lorenz混沌信號,添加0～20 dB的高斯白噪聲,然后分別用這些方法進行降噪處理。為減小偶然誤差,每種方法仿真200次,最后取200次的平均值作為最終的結(jié)果。

不同降噪方法的降噪后輸出信噪比和均方誤差的結(jié)果對比分別如圖4和圖5所示?；谛〔ǖ木植客队芭c基于EMD的局部投影的降噪效果接近,因為這兩種方法均是采用預先估計噪聲水平方式確定鄰域的。平滑子空間局部投影通過平滑正交分解識別子空間,效果優(yōu)于標準局部投影方法。遞歸局部投影采用遞歸圖計算最優(yōu)鄰域半徑,K均值局部投影利用K均值聚類方法選取鄰域,能夠?qū)崿F(xiàn)更好的去噪效果。而本文方法同時考慮了鄰域選取、子空間劃分以及質(zhì)心修正問題,相比上述方法,本文所提方法具有最高的輸出信噪比和最低的均方誤差,這說明方法較好地保留了信號時域波形的形狀,具有較好的降噪性能。

圖4 不同局部投影降噪方法的輸出信噪比Fig.4 Output signal to noise ratio of different local projective denoising methods

圖5 不同局部投影降噪方法的均方誤差Fig.5 Mean square error of different local projective denoising methods

綜上所述,本文所提的基于模糊遞歸圖和最優(yōu)硬閾值準則的局部投影方法能夠有效濾除混沌信號中的噪聲。

3.4 實測ECG信號仿真

ECG可以記錄心臟的生理活動,在醫(yī)療臨床領域應用廣泛。但是，ECG信號在儀器測量過程中容易被噪聲污染,從而影響各種生理特征的檢測、提取和識別,容易造成誤診,因而ECG信號降噪對實際臨床診斷具有重要的意義和價值。目前，常見的ECG信號去噪方法有奇異譜分析法、EMD方法、小波閾值方法等[31]。

ECG信號屬于非平穩(wěn)非線性信號,研究表明ECG信號具有復雜的動力學特性,表現(xiàn)出一定的混沌特性,因此非線性動力系統(tǒng)的方法作為一種研究心臟病的有效工具,被廣泛應用于ECG信號的研究當中。

本文選取MIT-BIH(Massachusettes Institute of Technology-Boston’s Beth Israel Hospital)噪聲壓力測試數(shù)據(jù)庫中編號為118e12的實測數(shù)據(jù),其采樣頻率為f=360 Hz,采集到的信號含有基線漂移、噪聲等影響。

首先，利用最小二乘擬合去掉信號內(nèi)的基線漂移；然后，再使用本文方法、奇異譜分析方法、EMD方法和小波閾值方法分別對ECG信號進行處理及對比。處理結(jié)果如圖6所示。

由圖6對比去噪前后的相圖,可以發(fā)現(xiàn)噪聲大大減少，本文方法去噪后的相圖表現(xiàn)出了較為清晰的吸引子結(jié)構(gòu),說明原來含有隨機噪聲的ECG信號在經(jīng)過降噪處理后，其確定性成分得以加強,展現(xiàn)出了ECG信號原來的特征。對比本文方法,奇異譜分析(singular spectrum analysis, SSA)方法、EMD方法和小波閾值(wavelet threshold, WT)方法雖然也有較好的去噪效果,但時域波形還是不夠光滑,而且相空間軌線細節(jié)部分呈現(xiàn)出雜亂堆疊的現(xiàn)象。

圖6 ECG信號降噪前后波形圖和相空間圖Fig.6 Waveform and phase space diagram of ECG signal before and after denoising

為進一步說明降噪效果,計算如圖7所示的原始信號和降噪信號在嵌入維數(shù)不同時的偽最近鄰點比例。當偽最近鄰點比例等于0或者小于某個值而保持不變,可認為不含偽最近鄰點;當時間信號受到噪聲影響時,混沌吸引子軌線受到影響,偽最近鄰點的比例也會受到影響。由圖7可知,降噪前偽最近鄰點的比例在嵌入維數(shù)為4時降到較低,但沒有降至0。降噪后,當嵌入維數(shù)為7時，偽最近鄰點的比例就降至0。這進一步說明了應用本文所提方法對實測ECG信號可進行有效的降噪處理。

圖7 ECG降噪前后的虛假最近鄰點比例Fig.7 Proportion of false nearest neighbor points before and after ECG noise reduction

4 結(jié) 論

本文提出的基于模糊遞歸圖和最優(yōu)硬閾值的局部投影算法,主要針對原始局部投影算法中的鄰域選擇問題和子空間劃分問題。在相空間重構(gòu)的基礎上,首先使用模糊遞歸圖對鄰域進行選擇,同時為避免求取原始方法中的協(xié)方差矩陣,本文直接使用SVD分解,然后采用最優(yōu)硬閾值準則對信號子空間和噪聲子空間進行劃分。通過對Lorenz混沌信號進行仿真,并與其他局部投影算法相比較,表明該方法具有一定的優(yōu)越性。最后，將其應用于實測含噪的ECG信號,有效地降低了原始信號中的噪聲,使信號的特征更加明顯。同時，將本文所提算法與其他類型的ECG信號降噪方法進行了對比,證明了本文方法的有效性,為后續(xù)的信號處理奠定了基礎。