基于EEMD閾值處理的腦電信號降噪方法

郭曉梅，朱曉軍

(太原理工大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，山西 晉中 030600)

0 引 言

BCI運(yùn)動想象(motor-imagery，MI)的一系列研究中，如何能夠獲得精準(zhǔn)信號對專家們探究大腦分析和加工信息的運(yùn)行機(jī)制具有至關(guān)重要的意義[1]。近幾年，越來越多的研究者們將自己的精力致力于尋找提高對類似腦電等非平穩(wěn)信號信噪比的更為有效的方法中。小波閾值法[2]、經(jīng)驗(yàn)?zāi)B(tài)分解法[3](empirical mode decomposition，EMD)、集合經(jīng)驗(yàn)?zāi)B(tài)分解法[4](ensemble empirical mode decomposition，EEMD)等方法的相繼出現(xiàn)，均在去噪領(lǐng)域取得了一定的效果。其中，小波閾值法在進(jìn)行小波分解前需選擇小波基函數(shù)以及設(shè)定分解層數(shù)，盲目性大。針對這一局限性，Wu和Huang提出了總體經(jīng)驗(yàn)?zāi)B(tài)分解，其可以根據(jù)信號本身特性來確定分解基函數(shù)與分解層數(shù)，自適應(yīng)性良好。在此基礎(chǔ)上發(fā)展起來的EEMD消噪法有時(shí)空濾波法[5]、基于EEMD分解的小波閾值法[6]等，雖然去噪性能有所改善，但仍留有不足。EMD閾值法[7]借鑒小波閾值法的原理，在信號處理中取得了一定的進(jìn)展，現(xiàn)已被廣泛應(yīng)用在對地震信號、振動信號等帶噪信號的除噪以及故障診斷[8-10]中。而在腦電研究中，為了避免EMD分解出現(xiàn)的端點(diǎn)效應(yīng)現(xiàn)象，有必要將閾值化思想引入到EEMD中作進(jìn)一步探討[11]。而平移不變小波去噪[12]利用平移不變算法的思想，通過平移噪聲信號使斷點(diǎn)移位，對閾值處理時(shí)引起的偽吉布斯現(xiàn)象有極大改善。

本文受到平移不變小波閾值去噪的啟發(fā)，將EEMD閾值去噪與平移不變算法相結(jié)合，成功應(yīng)用在對MI EEG的消噪中。通過對仿真信號和真實(shí)信號進(jìn)行實(shí)驗(yàn)驗(yàn)證，結(jié)果均表明：自適應(yīng)性良好的EEMD閾值化方法與平移不變算法的有效結(jié)合，不僅能夠抑制偽Gibbs現(xiàn)象，而且可以進(jìn)一步消除EEMD分解產(chǎn)生的模態(tài)混疊，使降噪性能達(dá)到更優(yōu)。

1 算法介紹

1.1 平移不變小波去噪算法

平移不變小波去噪法的主要思想是利用對帶噪信號移位的方式使得信號中所含斷點(diǎn)產(chǎn)生位移。該方法對于由于振蕩所產(chǎn)生的偽吉布斯現(xiàn)象能夠起到良好的緩解與抑制作用。

平移不變小波去噪法基本分為以下幾個(gè)步驟：

(1)首先設(shè)置好平移的范圍，按照既定位數(shù)對含噪信號進(jìn)行循環(huán)式平移；

(2)將經(jīng)過(1)操作后的含噪信號采取小波閾值處理；

(3)將經(jīng)過(2)操作后的信號再依照步驟(1)設(shè)定的平移位數(shù)作逆循環(huán)平移，移動到與原含噪信號相統(tǒng)一的相位；

(4)步驟(1)～(3)重復(fù)n次，將n次操作后的結(jié)果求平均值，即可得到降噪后的信號。

1.2 EEMD分解法原理

EEMD分解法同經(jīng)典的EMD算法一樣不存在選擇基函數(shù)的困擾，不同之處在于前者需要分多次向原始信號加入不一的白噪聲，從而形成新的含噪信號，然后分別對其進(jìn)行EMD分解。EEMD算法流程如圖1所示。其具體步驟為：

(1)將白噪聲序列n(t)添加到原始信號x(t)上，得到帶噪信號s(t)，即

s(t)=x(t)+n(t)

(1)

已知：n(t)為高斯白噪聲，且有n(t)～N(0，σ2)。

(2)對步驟(1)得到的帶噪信號s(t)采用EMD算法進(jìn)行分解，可以得到一系列IMF固有模態(tài)分量和剩余殘量rc(t)，即

(2)

式中：c為EMD分解后得到的IMF分量的數(shù)目。

(3)重復(fù)進(jìn)行步驟(1)～(2)m次，且步驟(1)中每次在原信號x(t)上加入的白噪聲序列n(t)的幅值是不一致的，即

(3)

(4)將m次步驟(2)中經(jīng)EMD分解生成的IMF求平均，即為最終IMF

(4)

在EEMD算法中所添加的高斯白噪聲需要滿足以下公式

(5)

其中：白噪聲的幅值由ε表示；所疊加的高斯白噪聲n(t)的次數(shù)由N表示；εn代表誤差大小，其是指EMD分解得到的各階IMF分量相加后與原始信號兩者之間的誤差。文獻(xiàn)[13]指出，當(dāng)N的取值落在100～300之間，且ε取信號標(biāo)準(zhǔn)偏差的0.01～0.5倍時(shí)，可以發(fā)現(xiàn)，此時(shí)的白噪聲誤差已然降到很低。于是，考慮到本文研究對象的特殊性，規(guī)定ε取值0.2，N為100。

圖1 EEMD分解法流程

2 EEMD分量閾值去噪法

根據(jù)Flandrin所提出的，含噪信號在經(jīng)過自適應(yīng)性良好的EEMD算法分解后，會生成多個(gè)IMF固有模態(tài)分量，并且每一個(gè)IMF分量當(dāng)中除了包含用于實(shí)驗(yàn)與分析的有用信號分量，還不可避免地?cái)y帶影響實(shí)驗(yàn)效果的噪聲分量。一般而言，在EEMD分解后，低頻段IMF分量中有效信號成分高占主導(dǎo)，噪聲信號大量存于高頻IMF分量中。依照這一理論，傳統(tǒng)的方法有EEMD時(shí)空濾波法，其認(rèn)為高頻IMF分量噪聲成分大，予以剔除。然而通常而言，存有大量噪聲的高頻IMF分量的判定還沒有統(tǒng)一標(biāo)準(zhǔn)，仍然需要進(jìn)一步考究，且部分高頻IMF分量中還留有少量對實(shí)驗(yàn)起重要作用的有用信號，所以使用EEMD時(shí)空濾波法對含噪信號進(jìn)行處理會導(dǎo)致信號失真，增大實(shí)驗(yàn)誤差，影響分析效果。于是學(xué)者們將研究重點(diǎn)轉(zhuǎn)移到EEMD分解產(chǎn)生的IMF分量上，繼續(xù)提出基于EEMD分解的小波閾值算法，EEMD分解后再采用小波閾值法對各階IMF分量做進(jìn)一步處理。該方法在提高信噪比方面確實(shí)有所突破，但是小波分解在基函數(shù)及分解層數(shù)上存在的不足依然無法避免。已有研究表明，借鑒小波閾值算法的閾值化處理思想，將其應(yīng)用在EEMD分解后的IMF分量間，這樣產(chǎn)生的EEMD分量閾值去噪法不僅保留了與EMD一致的自適應(yīng)性，而且有效克服了小波去噪依賴性大的缺陷。

2.1 互相關(guān)系數(shù)法確定信噪模態(tài)臨界點(diǎn)

通常情況下，實(shí)驗(yàn)采集到的EEG信號中包含有頻率較高的干擾噪聲。而經(jīng)EEMD分解后產(chǎn)生的眾多IMF分量會按照從高頻至低頻的排列順序聚集在某一固定頻段上。其中，高頻IMF分量中背景噪聲強(qiáng)度較大，而隨著分解層數(shù)的逐層增加，噪聲能量會越來越低，直到低頻IMF分量中有用信號占主導(dǎo)。若能夠確定前k個(gè)IMF分量主要包含噪聲成分，則僅對這前k個(gè)分量采取IMF閾值化法做進(jìn)一步處理。所以，如何準(zhǔn)確判定出噪聲為主的IMF與信號為主的IMF分量是EEMD分量閾值去噪法的關(guān)鍵，即需要找出信號與噪聲模態(tài)的臨界值k。通常認(rèn)為，EEMD分解后低頻段以信號為主要成分的IMF分量同原始信號的相似度更大，反之，噪聲能量大的高頻IMF分量與原始信號存在明顯差異，相似度較小。因此，本文選用互相關(guān)系數(shù)法，根據(jù)式(6)計(jì)算出各階IMF分量與原始信號的互相關(guān)程度大小

(6)

且有

(7)

可以根據(jù)計(jì)算出的互相關(guān)程度值繪制出每階IMF分量與原始信號的互相關(guān)系數(shù)曲線，并按照曲率計(jì)算公式得出該曲線的曲率值

(8)

R″=Ri+1-2Ri+Ri-1

(9)

R′=Ri-Ri-1

(10)

從曲線上不難看出，互相關(guān)系數(shù)首次發(fā)生轉(zhuǎn)折時(shí)橫軸對應(yīng)的位置就是信噪模態(tài)臨界點(diǎn)k的取值，或者可以理解為曲線第一次從緩慢下降變換到迅速增加的位置。在經(jīng)過曲率值的計(jì)算后，可以認(rèn)為IMF分量閾值去噪的臨界點(diǎn)即為第一個(gè)曲率極小值對應(yīng)的點(diǎn)，表示從第k+1個(gè)IMF分量開始均視為有用信號。

2.2 EEMD可導(dǎo)閾值函數(shù)

2.2.1 可導(dǎo)閾值函數(shù)的構(gòu)造

方便起見，令cj(k)=imfj(k)， 1≤k≤N，N代表信號長度。IMF分量閾值算法是受小波閾值的啟發(fā)而產(chǎn)生的。針對軟硬閾值各自的缺陷，采用一種連續(xù)且高階可導(dǎo)的新的EEMD閾值函數(shù)，閾值截?cái)嗪蟮腎MF如公式所示

(11)

2.2.2 自適應(yīng)閾值的選取

經(jīng)過EEMD分解后的閾值計(jì)算公式為

(12)

(13)

其中，cj(k)為第j層IMF系數(shù)值，length(cj(k))表示各層IMF的長度。該閾值大小可隨階數(shù)j的變化自適應(yīng)地進(jìn)行改變。

3 TI集成經(jīng)驗(yàn)?zāi)B(tài)分解自適應(yīng)閾值去噪法

為了預(yù)防奇異點(diǎn)在閾值處理時(shí)引發(fā)由于非自然原因?qū)е碌娜藶檎鹗幀F(xiàn)象，更深入地降低奇異點(diǎn)會對EEMD分解造成的影響，利用平移不變算法的思想，通過平移噪聲信號使斷點(diǎn)移位，對閾值處理時(shí)引起的偽吉布斯現(xiàn)象有極大改善，避免異常事件對EEMD造成困擾。

若將噪聲序列n(i)疊加在長度為N的原始信號s(i)上，形成的含噪信號f(i)為

f(i)=s(i)+n(i)，i=0,1,…,N-1

(14)

式中：n(i)為服從N(0,σ2)的高斯白噪聲。

最終設(shè)計(jì)的TI集成經(jīng)驗(yàn)?zāi)B(tài)分解自適應(yīng)閾值去噪法可以分為以下幾個(gè)步驟：

(1)首先設(shè)置好平移的范圍h，按照既定位數(shù)對含噪信號f(i)進(jìn)行循環(huán)式(左/右)平移；

(2)將經(jīng)過步驟(1)循環(huán)平移后的帶噪信號采用EEMD分解處理，得到按照從高頻到低頻的順序排列的IMF分量；

(3)用2.1節(jié)提出的互相關(guān)系數(shù)法確定IMF信噪模態(tài)臨界點(diǎn)，從而挑選出以噪聲為主的IMF分量，再應(yīng)用2.2節(jié)提出的EEMD可導(dǎo)閾值函數(shù)對其做量化處理；

(4)將未加工過的低頻IMF分量與經(jīng)過步驟(3)閾值量化后的IMF疊加以完成重構(gòu)；

(5)將經(jīng)過步驟(2)～步驟(4)操作后的信號再依照步驟(1)設(shè)定的平移位數(shù)h做逆循環(huán)(右/左)平移，移動到與原含噪信號相統(tǒng)一的相位；

(6)步驟(1)～步驟(5)重復(fù)n次，將n次操作后的結(jié)果求平均值，即可得到降噪后的信號。

假設(shè)原始信號為f(i)(0≤i≤N)，設(shè)置平移位數(shù)h，按照該既定位數(shù)循環(huán)平移得到Sh，規(guī)定快速TI算法中h等于1。整個(gè)過程可以表示為

(15)

圖2 本文去噪方法流程

4 實(shí)驗(yàn)與分析

4.1 去噪效果評定指標(biāo)

本小節(jié)旨在尋找能夠更加有效地刻畫去噪性能的判定標(biāo)準(zhǔn)，從較為全面更為直觀的角度來分析每一種方法的降噪效果。不僅僅局限于信噪比(SNR)和均方根誤差(RMSE)這兩種常規(guī)的評定指標(biāo)，為了防止單一化，進(jìn)一步引入了皮爾遜相關(guān)系數(shù)(ρ)和最大峰值誤差(MPE)兩種非常規(guī)性能指標(biāo)。

可以得出：SNR和ρ與去噪性能呈正相關(guān)，即值越大，去噪效果越好；RMSE和MPE與去噪性能呈負(fù)相關(guān)，即值越小，去噪效果則越好。SNR、RMSE、ρ、MPE的計(jì)算方法分別如公式(16)～式(19)所示

(16)

(17)

(18)

(19)

4.2 實(shí)驗(yàn)仿真

該小節(jié)旨在驗(yàn)證第2節(jié)提出的EEMD可導(dǎo)閾值函數(shù)的可行性以及第3節(jié)提出的TI集成經(jīng)驗(yàn)?zāi)B(tài)分解自適應(yīng)閾值去噪法的可靠性，現(xiàn)選擇在Matlab平臺上，按照真實(shí)EEG的波形特征，構(gòu)造采樣率Fs=250 Hz，采樣時(shí)間t=0：1/Fs：4，頻率范圍在2 Hz～30 Hz間的標(biāo)準(zhǔn)信號

(20)

對應(yīng)(δ、θ、α、β)4個(gè)腦電圖關(guān)鍵節(jié)律信號。疊加噪聲后的信號為

(21)

圖3 原始信號和加噪信號

圖6為使用不同算法對噪聲信號做消噪處理后的去噪圖。按照標(biāo)準(zhǔn)信號、加噪信號，依次由傳統(tǒng)硬、軟閾值法、EEMD時(shí)空濾波法、改進(jìn)小波閾值算法、基于EEMD的改進(jìn)小波閾值算法、EEMD可導(dǎo)閾值算法、TI集成經(jīng)驗(yàn)?zāi)B(tài)分解自適應(yīng)閾值算法等7種方法除噪后的信號圖的順序自上而下進(jìn)行排列。從圖6來看，整體而言，7種方法的消噪性能不論是依光滑度，還是從與標(biāo)準(zhǔn)信號特征逼近度出發(fā)，均得到了不同程度的改進(jìn)。對比來看，傳統(tǒng)硬閾值法是7種方法中效果最差的，而TI集成經(jīng)驗(yàn)?zāi)B(tài)分解自適應(yīng)閾值法可以更好地逼近未疊加噪聲前的標(biāo)準(zhǔn)信號。

圖4 EEMD分解后的IMF分量

圖5 各IMF分量與原始信號的互相關(guān)系數(shù)及對應(yīng)的曲率曲線

圖6 不同方法下帶噪信號去噪

為了更加直觀精確地評估本文提出的方法的去噪性能，并與傳統(tǒng)硬、軟閾值法、改進(jìn)小波閾值算法、基于EEMD的改進(jìn)小波閾值算法等以上5種算法的去噪結(jié)果進(jìn)行比較分析，對于閾值處理時(shí)需要關(guān)心的閾值函數(shù)及閾值的選取原則均按照本文提及到的方法來操作。

現(xiàn)將信噪比分別為1 dB、5 dB、10 dB、15 dB、20 dB這5個(gè)級別的噪聲輪流疊加到原始信號上。因?yàn)榧词乖诏B加同一信噪比的情況下，實(shí)驗(yàn)多次得到的去噪結(jié)果都不會相同，故考慮在同一信噪比下重復(fù)做20組實(shí)驗(yàn)，將不同方法下消噪后的信噪比、均方根誤差、皮爾遜相關(guān)系數(shù)、最大峰值誤差4個(gè)指標(biāo)計(jì)算并統(tǒng)計(jì)出來，然后對每種方法下的20組實(shí)驗(yàn)結(jié)果取平均值。各種算法在不同信噪比下的消噪效果見表1。方法(一)～方法(七)對應(yīng)圖6中的7種算法。

分析表1的數(shù)據(jù)可知，方法(一)～方法(三)有效地驗(yàn)證了改進(jìn)后的可導(dǎo)閾值函數(shù)較傳統(tǒng)閾值函數(shù)去噪性能有所提高。疊加信噪比為10 dB時(shí)，不同方法之間性能指標(biāo)有明顯差距：單從去噪后信噪比的角度來看，EEMD可導(dǎo)閾值法比基于EEMD的改進(jìn)小波閾值法提高約0.3409 dB，而本文提出的方法在此基礎(chǔ)上SNR又提高了約0.136 dB；在皮爾遜相關(guān)系數(shù)指標(biāo)上，本文較基于EEMD改進(jìn)小波閾值法提高約0.0007，提升幅度大于EEMD可導(dǎo)閾值法；同樣，在均方根誤差RMSE和最大峰值誤差MPE兩大標(biāo)準(zhǔn)上，本文去噪方法較其他幾種算法降低幅度最大。綜合而言，同一信噪比下，采用本文提出的方法對帶噪信號進(jìn)行處理后，SNR和ρ值最大，RMSE和MPE值最小，消噪效果最為理想，可以在有效剔除干擾的同時(shí)最大化保留信號的原始特征。

4.3 實(shí)際腦電實(shí)驗(yàn)

經(jīng)過4.2節(jié)對仿真信號所進(jìn)行的實(shí)驗(yàn)，驗(yàn)證了本文提出的降噪方法已具備一定的可行性。為了再進(jìn)一步驗(yàn)證其在真正經(jīng)過實(shí)驗(yàn)采集到的腦電信號上的有效性，現(xiàn)選用奧地利GRAZ大學(xué)公開共享的左右手運(yùn)動想象腦電(MI EEG)數(shù)據(jù)集III作為本節(jié)重點(diǎn)研究對象[14]。該數(shù)據(jù)集分140個(gè)有標(biāo)識的訓(xùn)練集與140個(gè)未被標(biāo)記的測試集兩部分。受試者是一名25歲女性，頭部佩戴電極帽，要求其在規(guī)定的時(shí)間內(nèi)完成對左手和右手兩種運(yùn)動任務(wù)的想象。該實(shí)驗(yàn)采樣率Fs=128 Hz，按照每組40次實(shí)驗(yàn)共執(zhí)行7組，共280次實(shí)驗(yàn)。每一次實(shí)驗(yàn)對應(yīng)有C3、C4和Cz 這3個(gè)通道采集到的MI腦電數(shù)據(jù)。

本文采用不同降噪方法分別對帶標(biāo)識訓(xùn)練集中3個(gè)通道上的MI EEG信號做預(yù)處理。然后將每一通道下的140組數(shù)據(jù)在經(jīng)過各種算法去噪后的SNR、RMSE、ρ、MPE

表1 不同算法去噪效果定量評價(jià)指標(biāo)結(jié)果對比

這4個(gè)指標(biāo)值統(tǒng)計(jì)出來，最后再計(jì)算出每個(gè)評價(jià)指標(biāo)所對應(yīng)的均值和標(biāo)準(zhǔn)差大小。

現(xiàn)從C4導(dǎo)聯(lián)上采集到的原始MI EEG信號中任取一段作為去噪實(shí)驗(yàn)的研究信號。原始腦電信號圖如圖7(a)所示。對該信號依次采用傳統(tǒng)硬、軟閾值、EEMD時(shí)空濾波法、改進(jìn)小波閾值、基于EEMD的改進(jìn)小波閾值、EEMD可導(dǎo)閾值及本文提出的去噪算法進(jìn)行處理，得到消噪后的效果圖如圖7(b)～圖7(h)所示。

圖7 原始MI EEG信號及不同算法消噪效果對比

分析圖7可發(fā)現(xiàn)，對比方法圖7(b)～圖7(d)，方法圖7(h)在對信號的處理中丟失的信息較少，能夠較全面地保留原始信號特征；同方法圖7(e)～圖7(g)相比，方法圖7(h)處理后的信號更加平滑，且從視覺效果上而言，模態(tài)混疊和偽Gibbs現(xiàn)象得到了改善。因此，從整體把握上來看，本文所提算法去噪效果達(dá)到更優(yōu)。同時(shí)，為了使結(jié)果一目了然且更具備說服力，這里統(tǒng)計(jì)出了C4導(dǎo)聯(lián)下的140組數(shù)據(jù)在經(jīng)過各種算法去噪后的SNR、RMSE、ρ、MPE這4個(gè)評價(jià)標(biāo)準(zhǔn)所對應(yīng)的均值和標(biāo)準(zhǔn)差大小，見表2。

分析表2的數(shù)據(jù)可知，改進(jìn)后的可導(dǎo)閾值函數(shù)較傳統(tǒng)閾值函數(shù)去噪性能有所提高：在SNR這一指標(biāo)上，其較硬閾值法提高約5.9807 dB，較軟閾值法提高約5.2435 dB；從皮爾遜相關(guān)系數(shù)的角度來看，改進(jìn)后的閾值函數(shù)較硬閾值法提高約0.0571，較軟閾值法提高約0.0614；同樣，在RMSE和MPE兩大指標(biāo)上，采用新閾值函數(shù)和閾值選取規(guī)則較傳統(tǒng)軟硬閾值函數(shù)有更大幅度降低。基于EEMD的改進(jìn)小波閾值法在采用新的可導(dǎo)閾值函數(shù)的基礎(chǔ)上結(jié)合了自適應(yīng)性較強(qiáng)的EEMD分解法，在去噪性能上較改進(jìn)的小波閾值法更優(yōu)。而EEMD可導(dǎo)閾值法有效避免了基于EEMD的改進(jìn)小波閾值法在小波分解上的困擾，且從數(shù)據(jù)上顯示，其去噪效果良好。本文提出的算法將EEMD可導(dǎo)閾值法與平移不變算法有效結(jié)合，采用本文提出的方法對MI EEG信號進(jìn)行處理后，跟其余方法相比，SNR和ρ值最大，RMSE和MPE值最小，消噪效果最為理想，可以在有效剔除干擾的同時(shí)最大化保留MI EEG信號的原始特征。

表2 采用不同算法對C4通道的140組數(shù)據(jù)進(jìn)行處理的實(shí)驗(yàn)結(jié)果

5 結(jié)束語

為改善對非線性非平穩(wěn)信號的去噪性能，本文受到平移不變方法思想的啟發(fā)，將其與EEMD閾值去噪法相結(jié)合，提出一種TI總體經(jīng)驗(yàn)?zāi)B(tài)分解自適應(yīng)閾值處理的EEG去噪方法，并成功應(yīng)用在對MI EEG的消噪中。通過對仿真信號和真實(shí)信號進(jìn)行實(shí)驗(yàn)驗(yàn)證，結(jié)果均表明：自適應(yīng)性良好的EEMD閾值化方法與平移不變算法的有效結(jié)合，不僅能夠抑制偽Gibbs現(xiàn)象，而且可以進(jìn)一步消除EEMD分解產(chǎn)生的模態(tài)混疊，使降噪性能達(dá)到更優(yōu)。本文方法的提出不僅在腦電信號有效信息的精確保留上占據(jù)優(yōu)勢，為后續(xù)特征提取與模式分類等方面提供關(guān)鍵的數(shù)據(jù)保障，更為去噪領(lǐng)域樹立了新航標(biāo)。