基于貝葉斯估計(jì)的小波腦電信號(hào)去噪算法研究

王宏旭，張晨潔，劉勇，郭濱

（長春理工大學(xué) 電子信息工程學(xué)院，長春 130022）

大腦是由億萬個(gè)神經(jīng)元組成的復(fù)雜系統(tǒng)，負(fù)責(zé)身體的各個(gè)功能的協(xié)調(diào)運(yùn)作，通過大腦皮層上電極記錄下的大腦細(xì)胞群的電位活動(dòng)稱為腦電信號(hào)[1]。通過對(duì)腦電信號(hào)進(jìn)行研究可以獲得豐富的心理及生理的疾病信息，所以腦電信號(hào)的分析及去噪算法的研究無論是在臨床的診斷和急病治療上都是十分重要的。一般來說，腦電信號(hào)具有背景噪聲強(qiáng)、信號(hào)微弱等特點(diǎn)，所以如何消除腦電數(shù)據(jù)的噪聲，更好地獲取對(duì)大腦有用的信息是當(dāng)今研究的熱門話題。

近年來小波變換廣泛應(yīng)用于信號(hào)處理和圖像融合等領(lǐng)域。小波去噪的方法主要有模極大值法、空間濾波、閾值法，其中閾值法國內(nèi)外學(xué)者廣泛應(yīng)用的方法[2]。閾值的選取是根據(jù)信號(hào)的長度設(shè)定，其設(shè)定的方法具有單一性，無法一次性對(duì)信號(hào)參數(shù)進(jìn)行選定，由于選取過于零散，無法滿足信號(hào)復(fù)雜等特點(diǎn)。所以針對(duì)這一問題，很多學(xué)者提出閾值選取的方法，文獻(xiàn)[3]提出一種全局閾值選取方法，根據(jù)分解層數(shù)的不同對(duì)腦電信號(hào)進(jìn)行去噪，雖然在信號(hào)中混入基線漂移，但是還原了腦電特征；文獻(xiàn)[4]利用小波變換對(duì)EEG信號(hào)進(jìn)行高頻和低頻系數(shù)分解，并且利用改進(jìn)的閾值選取方法，根據(jù)分解層數(shù)，對(duì)小波變換后的小波系數(shù)進(jìn)行自適應(yīng)閾值處理，去噪效果優(yōu)于硬閾值、軟閾值、Garrote閾值；文獻(xiàn)[5]設(shè)計(jì)了閾值選取器，結(jié)合最小均方誤差和小波去噪方法，提出一種LMS算法來自適應(yīng)控制閾值參數(shù)，可較大程度減少噪聲，提高信噪比，且該方法比傳統(tǒng)去噪方法更準(zhǔn)確。

實(shí)驗(yàn)研究表明腦電去噪的效果好壞對(duì)大腦功能的分析和疾病診斷具有重要的意義，所以本文利用美國波恩醫(yī)學(xué)中心的癲癇患者實(shí)測(cè)的腦電時(shí)間序列進(jìn)行去噪處理，通過觀察噪聲特征模型，采用拉普拉斯最大后驗(yàn)概率進(jìn)行估計(jì)，考慮到小波系數(shù)的長拖尾性和噪聲干擾性，提出了一種貝葉斯估計(jì)的自適應(yīng)算法，相比于小波閾值去噪和其他的去噪方法，本文的去噪算法有效的避免了常用的小波去噪造成的細(xì)節(jié)模糊和信息丟失，實(shí)驗(yàn)研究顯示，本文的算法可較大程度保留原始腦電信號(hào)特征同時(shí)具有良好的去噪性能，并且可以有效地提高峰值信噪比。

1 小波基本理論

1.1 小波變換

小波變換是對(duì)基本小波進(jìn)行伸縮和平移得到的，小波系數(shù)是原信號(hào)與小波基函數(shù)的相似系 數(shù)[6]。 設(shè)函數(shù)Φ（t）為可積函數(shù)且滿足t=L2（R），則傅里葉變換滿足如下條件：

則Φ（t）為小波母函數(shù)，將小波母函數(shù)進(jìn)行伸縮和平移得到：

式中，a為尺度系數(shù)；b為平移系數(shù)；Φa，b（t）是經(jīng)過伸縮平移后得到的連續(xù)小波基函數(shù)。信號(hào)f（t）在式（2）中進(jìn)行展開，這樣的變換為連續(xù)小波變換，簡稱CWT，表達(dá)式為：

式中，a大于0，在CWT中尺度參數(shù)和平移參數(shù)進(jìn)行離散采樣后，得到離散小波變換，其表達(dá)式為：

1.2 小波基的選擇

小波基的選擇對(duì)腦電去噪效果至關(guān)重要，一般對(duì)稱性好的小波基不會(huì)產(chǎn)生相位畸變[8]。正則性好的小波基在進(jìn)行腦電重構(gòu)時(shí)腦電信號(hào)更加平滑，緊支撐的小波在腦電去噪時(shí)處理的速度更快。

本文選擇的小波基為dbNv小波，可以滿足以上的條件，db1（簡稱Haar）小波是小波基當(dāng)中最簡單的小波形式，局域性較差，因此本文選擇db6小波，db6的特點(diǎn)是處理速度更快，MATLAB軟件應(yīng)用方便，其波形和腦電波形相似，對(duì)腦電去噪效果最好。db6小波曲線如圖1所示。

圖1 db6小波函數(shù)

1.3 小波去噪原理

小波去噪原理是將小波進(jìn)行分解，再對(duì)小波進(jìn)行多尺度變換，盡可能提取有用的腦電信號(hào)[7]。然后再根據(jù)波恩腦電信號(hào)特征和噪聲特點(diǎn)，選擇合適的去噪模型，利用貝葉斯估計(jì)后的系數(shù)進(jìn)行小波重構(gòu)，從而得到去噪后的腦電信號(hào)。其基本原理如圖2所示，在進(jìn)行小波變換前，腦電信號(hào)的突變處峰值很高，突變處集中了噪聲，很多平穩(wěn)有用的腦電信號(hào)幅值很小，從而大大增加了去噪難度，所以本文以拉普拉斯為去噪模型，結(jié)合貝葉斯進(jìn)行小波估計(jì)，已達(dá)到去噪的目的。

圖2 腦電去噪原理

2 拉普拉斯最大后驗(yàn)概率模型

由于拉普拉斯分布模型具有尖峰脈沖特性以及嚴(yán)重的拖尾特性、概率密度分布與癲癇患者腦電信號(hào)的概率密度分布相似，所以本文采用拉普拉斯為去噪模型，在利用貝葉斯估計(jì)方法進(jìn)行腦電去噪，信號(hào)模型的先驗(yàn)準(zhǔn)確性直接影響去噪效果，所以模型的選取至關(guān)重要，拉普拉斯分布去噪模型為：

式中，g為含噪的腦電信號(hào)；s為真實(shí)的腦電信號(hào)；ε為噪聲信號(hào)。其中 ε服從N（0，?2）高斯噪聲分布，對(duì)于含噪的腦電信號(hào)g進(jìn)行相應(yīng)的小波變換，得到：

式中，y為含噪腦電信號(hào)；w為不含噪聲的理想腦電信號(hào)；n為噪聲的小波系數(shù)。對(duì)于小波變換n服從N(0，σ2)分布，拉普拉斯概率分布與噪聲的高斯分布的概率密度函數(shù)非常相似，所以對(duì)于理想腦電信號(hào)w在0點(diǎn)附近有很長的峰值，并且有拖尾現(xiàn)象發(fā)生，這和高斯分布相符合，其拉普拉斯概率密度公式為：

式中，w為小波系數(shù)的先驗(yàn)?zāi)Ｐ?；σ為尺度參?shù)。該公式可以有效的估計(jì)噪聲模型，對(duì)去噪效果至關(guān)重要。

3 貝葉斯估計(jì)基礎(chǔ)理論

3.1 貝葉斯原理

在腦電信號(hào)中通常認(rèn)為噪聲都是相互獨(dú)立的，經(jīng)過小波變換后的噪聲n的概率密度函數(shù)為：

最大后驗(yàn)估計(jì)（MAP）在貝葉斯理論中是一種常用的方法，在給定的腦電數(shù)據(jù)中，p（y|x）為似然密度函數(shù)，p（x）被稱為先驗(yàn)密度函數(shù)，所以為了求p（y|x）的最大值，根據(jù)貝葉斯估計(jì)方法，可以將貝葉斯公式轉(zhuǎn)化為：

在給定腦電信號(hào)y的條件下，使概率密度函數(shù)pw|y(w|y)最大的w，其表示含噪腦電信號(hào)為：

根據(jù)公式（6）的貝葉斯準(zhǔn)則變換得到：

對(duì)公式（11）兩邊取對(duì)數(shù)，則得到如下公式：

將式（7）和式（8）代入上式，得到：

腦電信號(hào)的高斯白噪聲符合均值μ=0，標(biāo)準(zhǔn)差為σ2=1，對(duì)上式w進(jìn)行求導(dǎo)，得到公式：

另|w|=sign(w)，并設(shè)w=y則公式可轉(zhuǎn)化為：

3.2 貝葉斯自適應(yīng)去噪算法

為了使腦電信號(hào)去噪算法具有自適應(yīng)性，可以將小波系數(shù)看作拉普拉斯的分布模型，根據(jù)公式（14）可估計(jì)真實(shí)腦電信號(hào)的小波系數(shù)w，在知道標(biāo)準(zhǔn)差σ和噪聲方差σ2的情況下，噪聲方差可描述噪聲的統(tǒng)計(jì)特性，為了對(duì)噪聲方差做出更好的估計(jì)，1995年Donoho和Johnstone等人[9]提出小波域魯棒中值法對(duì)噪聲方差進(jìn)行估計(jì)。表達(dá)式是：

因信號(hào)與噪聲相互獨(dú)立，所以公式中y(i)是第一級(jí)小波系數(shù)分解后對(duì)子帶的分解系數(shù)。要得到真實(shí)的腦電信號(hào)，首先對(duì)每個(gè)含噪信號(hào)進(jìn)行方差估計(jì)，其公式為：

對(duì)含噪腦電信號(hào)的標(biāo)準(zhǔn)差估計(jì)為：

去噪后的小波系數(shù)具有不同的拉普拉斯分布，信號(hào)的邊緣標(biāo)準(zhǔn)差具有相關(guān)的隨機(jī)性，所以領(lǐng)域的估計(jì)系數(shù)標(biāo)準(zhǔn)差為：

根據(jù)公式推導(dǎo)基于貝葉斯自適應(yīng)去噪算法流程為：

（1）對(duì)含噪的腦電信號(hào)進(jìn)行小波變換；

（2）利用公式（15）進(jìn)行噪聲方差估計(jì)；

（3）按照公式（16）和公式（17）得到拉普拉斯分布模型的標(biāo)準(zhǔn)差估計(jì)，再利用公式（14）進(jìn)行處理，最后得到真實(shí)小波系數(shù)的貝葉斯估計(jì)；

（4）進(jìn)行小波反變換進(jìn)行重構(gòu)，得到去噪后的腦電信號(hào)。

4 腦電去噪性能的評(píng)價(jià)方法

4.1 信噪比和均方根誤差

腦電去噪效果在主觀和客觀評(píng)價(jià)主要依靠觀察者的主觀感覺和客觀評(píng)價(jià)方法。為了改進(jìn)算法的優(yōu)越性，本文采用信噪比（SNR）和均方根誤差（RMSE）兩個(gè)指標(biāo)來檢驗(yàn)貝葉斯估計(jì)去噪算法，其公式為：

式中，n為信號(hào)的長度；s(i)為原始信號(hào)；x(i)為去噪后的腦電信號(hào)。通常情況下信噪比越高去噪效果越好，而均方誤差越低表明去噪后的腦電數(shù)據(jù)精確度越高，去噪效果越好。

4.2 功率譜估計(jì)

經(jīng)典的功率譜估計(jì)需要掌握概率論的先驗(yàn)方法，而這些先驗(yàn)知識(shí)和方法往往是由人的主觀思維決定的，所以往往會(huì)導(dǎo)致頻率低等缺點(diǎn)[10]?，F(xiàn)在的功率譜估計(jì)方法主要以隨機(jī)過程的參數(shù)模型為基礎(chǔ)，對(duì)加窗函數(shù)進(jìn)行假設(shè)，腦電功率譜估計(jì)方法大致分為以下幾個(gè)步驟：

（1）被估計(jì)的隨機(jī)過程要選擇一個(gè)模型，本文選自回歸模型，即AR模型；

（2）根據(jù)觀測(cè)數(shù)據(jù)模型進(jìn)行白化處理；

（3）根據(jù)估計(jì)后的參數(shù)模型進(jìn)行功率譜估計(jì)。

5 實(shí)驗(yàn)結(jié)果及分析

本文的腦電去噪實(shí)驗(yàn)在MATLAB2017a中進(jìn)行，建立算法的M文件，通過軟件對(duì)波恩醫(yī)學(xué)中心腦電數(shù)據(jù)進(jìn)行去噪實(shí)驗(yàn)，腦電去噪用的樣本采樣頻率大致為173.61 Hz，采樣時(shí)間為23.6 s，采樣數(shù)據(jù)點(diǎn)個(gè)數(shù)為4 096個(gè)，為了計(jì)算方便，本文截取1 024個(gè)點(diǎn)進(jìn)行分析。如圖2所示為癲癇患者腦電數(shù)據(jù)圖。從圖中可以看出，癲癇患者的腦電波形持續(xù)時(shí)間小于70 ms則被稱為棘波，70～200 ms之間常被稱為尖波。

5.1 實(shí)驗(yàn)分析

實(shí)驗(yàn)1：通過第2層、第3層、第4層和第5層進(jìn)行近似分解，如圖4所示為全局閾值去噪波形圖，小波分解過程中分解層數(shù)決定去噪效果，通過圖3和圖4可以看出，全局閾值去噪效果明顯，分解層數(shù)越多，去噪效果越好，但是計(jì)算量會(huì)相應(yīng)的增大，癲癇患者腦電信號(hào)去噪邊緣相對(duì)平滑，邊緣丟失了腦電數(shù)據(jù)信息，且不能很好的保留原始腦電信號(hào)特征。

圖3 波恩醫(yī)學(xué)中心腦電數(shù)據(jù)圖

圖4 全局閾值去噪

實(shí)驗(yàn)2：小波變換對(duì)腦電信號(hào)進(jìn)行高頻和低頻系數(shù)分解，并且利用改進(jìn)的閾值選取方法進(jìn)行癲癇患者腦電去噪。如圖5所示為改進(jìn)閾值方法進(jìn)行腦電信號(hào)降噪，該方法不僅抑制了高斯噪聲，相比于全局閾值去噪，改進(jìn)的閾值方法保留了大部分有用的細(xì)節(jié)信息，信噪比較高，但是計(jì)算量大，不能保留全部的腦電信號(hào)。

圖5 改進(jìn)閾值法

實(shí)驗(yàn)3：采用了自適應(yīng)閾值選取器，根據(jù)腦電信號(hào)噪聲強(qiáng)度的不同，在去噪過程中自適應(yīng)的選取閾值，已達(dá)到最優(yōu)閾值的目的。如圖6所示為自適應(yīng)閾值去噪方法，該方法有計(jì)算量低、運(yùn)算速度快等特點(diǎn)，且實(shí)驗(yàn)研究表明，此方法優(yōu)于傳統(tǒng)方法，能較好地提取有用的癲癇患者腦電信號(hào)，更好地保留原始信號(hào)的細(xì)節(jié)信息。

圖6 自適應(yīng)閾值法

針對(duì)以上實(shí)驗(yàn)去噪方法存在的不足，本文設(shè)計(jì)了貝葉斯估計(jì)自適應(yīng)去噪方法，如圖7所示，從圖中可以看出該方法幾乎保留了癲癇患者腦電的高頻部分和低頻部分的所有細(xì)節(jié)信息，從肉眼可以看出該波形與原始腦電信號(hào)的波形相似度較高，光滑性較好，計(jì)算量低。

圖7 貝葉斯自適應(yīng)法

5.2 仿真結(jié)果分析

通過實(shí)驗(yàn)的對(duì)比可知貝葉斯算法去噪效果明顯優(yōu)于其他算法，在實(shí)驗(yàn)中也可對(duì)SNR和RMSE的值進(jìn)行測(cè)試。從表1中對(duì)比SNR和RMSE的值，其中貝葉斯算法的值最高，同時(shí)RMSE的值最低，所以從數(shù)值和評(píng)價(jià)指標(biāo)上表明了貝葉斯算法對(duì)腦電去噪效果優(yōu)于其他算法。

表1 不同的算法去噪效果的SNR和RMSE結(jié)果對(duì)比

在對(duì)腦電進(jìn)行去噪過程中，可以通過能量比對(duì)去噪效果進(jìn)行評(píng)估。從圖8可以得知貝葉斯算法的能量比高于其他算法，表明貝葉斯算法去噪效果最好。

圖8 能量比指標(biāo)對(duì)比

本文分別對(duì)不同算法的功率譜進(jìn)行估計(jì)，結(jié)果如圖9所示，從圖可以看出全局閾值、改進(jìn)閾值、自適應(yīng)閾值的功率譜形狀呈鋸齒形，譜峰點(diǎn)位置不確定，而貝葉斯方法的功率譜更加平坦，去噪性能更好。

圖9 不同算法的功率譜估計(jì)

6 結(jié)論

本文針對(duì)波恩癲癇腦電數(shù)據(jù)及噪聲模型提出了一種以拉普拉斯為最大后驗(yàn)概率估計(jì)，再對(duì)腦電信號(hào)進(jìn)行小波變換，利用貝葉斯算法估計(jì)小波變換后的系數(shù)，從而達(dá)到腦電去噪的目的。從貝葉斯估計(jì)算法、全局閾值算法、改進(jìn)閾值算法和自適應(yīng)閾值算法對(duì)比來看，貝葉斯估計(jì)算法取得的信噪比最高，均方根誤差最低，從能量比和功率譜估計(jì)來看，貝葉斯算法不僅結(jié)構(gòu)簡單，計(jì)算速度快，而且?guī)缀醣Ａ羲械母哳l部分和低頻部分的細(xì)節(jié)信息，沒有對(duì)有用的腦電信號(hào)頻譜造成傷害，并取得了較大的信噪比增益。有利于未來推動(dòng)臨床醫(yī)學(xué)和疾病診斷的研究和預(yù)防。