宋 婷，舒智林，孫玉波，韓建達(dá)，于寧波?

(1.南開大學(xué)人工智能學(xué)院，天津 300350；2.南開大學(xué)天津市智能機器人技術(shù)重點實驗室，天津 300350)

腦電圖(Electroencephalogram，EEG)是大腦神經(jīng)細(xì)胞群活動在大腦皮層的綜合反映，是頭皮表面記錄到的大腦神經(jīng)元產(chǎn)生的電活動，用于量化大腦電活動狀態(tài)，廣泛應(yīng)用于腦機接口和疾病診斷等[1－4]。隨著技術(shù)的發(fā)展，腦電信號采集設(shè)備更加便攜，通道數(shù)目也在不斷減少，甚至為單通道，研究單通道腦電信號偽跡去除算法對于準(zhǔn)確有效地分析大腦活動是至關(guān)重要的[5－6]。

在腦電信號采集過程中，位于額葉區(qū)域的腦電通道極易受到眼電偽跡(Electrooculogram，EOG)的干擾[7－8]。眼電信號的幅值遠(yuǎn)高于腦電信號，且頻譜與腦電信號的δ頻段重疊[9]。針對上述問題，研究者們已經(jīng)提出了從腦電信號中去除眼電偽跡的各種方法，包括回歸方法、盲源分離(Blind Source Separation，BSS)、小波變換(Wavelet Transform，WT)和經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)等[10－11]。回歸方法是單通道腦電信號中去除眼電偽跡最直接的方法，需要同時采集腦電和眼電信號，增加了硬件復(fù)雜度，且存在腦電和眼電雙向污染的問題[12]。盲源分離是將不同來源的多通道信號分解為多個獨立分量，去除偽跡相關(guān)的獨立分量。該算法對于多通道腦電信號具有較好的效果，但由于其要求通道的數(shù)量必須大于等于源的數(shù)量，所以對于單通道腦電信號并不適用[13]。小波變換是腦電信號中偽跡分離的重要方法，但是存在母波選擇、分解層數(shù)設(shè)置以及在頻譜重疊情況下去噪效果不佳的問題[14－15]。經(jīng)驗?zāi)B(tài)分解是將腦電信號分解為有限個不同時間尺度的本征模態(tài)函數(shù)(Instrinsic Mode Functions，IMFs)，去除以眼電偽跡成分為主的IMFs。該算法可以在無任何先驗知識的情況下對信號進(jìn)行自適應(yīng)分解，對于腦電等非平穩(wěn)信號的處理具有突出優(yōu)勢，但是對噪聲比較敏感，容易造成模態(tài)混疊的問題[16]。

針對單一方法的優(yōu)缺點，研究人員開發(fā)了結(jié)合兩種或多種方法的算法。Maddirala 等[17]將奇異譜分析(Singular Spectrum Analysis，SSA)得到的參考信號用于眼電偽跡的自適應(yīng)噪聲去除(Adaptive Noise Cancellation，ANC)，用于解決單通道腦電信號眼電偽跡去除問題。該算法存在模態(tài)混疊及嵌入維度等參數(shù)選擇的問題。羅志增等[18]提出了一種將自適應(yīng)完備經(jīng)驗?zāi)B(tài)分解(Complete Ensemble Empirical Mode Decomposition with Adaptive Noise，CEEMDAN) 和獨立成分分析 (Independent Component Analysis，ICA)相結(jié)合的算法，用于單通道腦電信號中眼電偽跡的去除。雖然CEEMDAN 將單通道信號分解為多維數(shù)據(jù)，滿足了BSS 的先驗條件，但是腦電和眼電信號存在頻譜混疊，會將腦電和眼電信號分解到同一個IMF，直接去除IMF 會導(dǎo)致腦電信號失真。張銳等[19]結(jié)合小波變換和集合經(jīng)驗?zāi)B(tài)分解(Ensemble EMD，EEMD)兩種方法自動去除單通道腦電信號中的眼電偽跡。雖然算法復(fù)雜度低，但是沒有給出含眼電偽跡的小波成分選擇標(biāo)準(zhǔn)，且EEMD 算法會存在噪聲殘留的問題。

為此，本文針對腦電信號和眼電偽跡存在頻譜混疊的問題，創(chuàng)新性地提出一種單通道腦電信號中眼電偽跡自動去除算法。采用經(jīng)驗小波變換(Empirical Wavelet Transform，EWT)[20]將單通道腦電信號分解為δ頻段和高頻段信號，采用改進(jìn)的自適應(yīng)噪聲完備經(jīng)驗?zāi)B(tài)分解(Improved CEEMDAN，ICEEMDAN)[21]將得到的δ頻段腦電信號分解為多維本征模態(tài)函數(shù)IMFs，設(shè)置樣本熵閾值自動去除以眼電偽跡成分為主的IMFs，最后重構(gòu)得到濾波后的腦電信號，使用半模擬腦電數(shù)據(jù)和真實的腦電數(shù)據(jù)對算法進(jìn)行驗證，并與其他方法對比。

1 方法

本文提出的單通道腦電信號眼電偽跡去除算法EWT-ICEEMDAN 的流程如圖1 所示。該算法首先采用EWT 算法將含眼電偽跡的單通道腦電信號分解為δ頻段和高頻段信號，然后采用ICEEMDAN 算法將δ頻段信號分解為有限個IMFs，進(jìn)一步設(shè)置樣本熵閾值去除含眼電偽跡的IMFs，最后重構(gòu)得到去除眼電偽跡的單通道腦電信號。

圖1 單通道腦電信號眼電偽跡去除算法流程圖

1.1 基于經(jīng)驗小波變換的腦電信號分解

EWT 通過生成自適應(yīng)濾波器組來提取輸入非平穩(wěn)信號的模態(tài)分量，不僅能克服EMD 噪聲殘留的問題，還能根據(jù)信號頻譜的不同調(diào)整濾波器組的帶寬，克服離散小波變換(Discrete Wavelet Transform，DWT)分解信號濾波器組恒定的問題[20]。研究顯示眼電偽跡主要影響腦電信號的δ頻段，當(dāng)腦電信號中存在眼電偽跡時，δ頻段的功率譜密度較高[17]。因此，本文采用EWT 提取腦電信號的δ頻段和高頻段信號。EWT 分解過程如下:

首先，對于給定信號f(t)通過傅里葉變換得到傅里葉頻譜，假設(shè)將傅里葉支撐區(qū)間[0，π]分割成N個連續(xù)的區(qū)間，各個區(qū)間可表示為Λn＝[ωn－1，ωn]，其中，n＝1，2，…，N，并且滿足＝[0，π]。以ωn為中心，定義過渡相Tn＝2τn，其中，τn＝γωn(0<γ<1)，γ滿足式(1):

式中:β(x)的表達(dá)式為[22]:

最后，得到近似小波系數(shù)Wf(0，t)和細(xì)節(jié)小波系數(shù)Wf(n，t)，公式如下:

原始信號f(t)可被重構(gòu)為:

EWT 通過對信號頻譜進(jìn)行顯著的分割來創(chuàng)建自適應(yīng)濾波器組，然后將單通道腦電信號的傅里葉頻譜分割為δ頻段和高頻段，對應(yīng)的頻率范圍分別為[0，4]Hz 和[4，40]Hz，最后得到δ頻段信號fδ(t)和高頻段信號fhigh(t)。

1.2 基于ICEEMDAN 的δ 頻段腦電信號分解

ICEEMDAN 不僅能解決模態(tài)混疊的問題，還能解決CEEMDAN 中的殘留噪聲和偽模態(tài)問題[21]。本文采用ICEEMDAN 對δ頻段的腦電信號fδ(t)進(jìn)行分解。設(shè)Ek(?)為EMD 分解得到第k個IMF 分量的算子，w(t)(i)為服從N(0，1)的白噪聲。ICEEMAN 算法的流程如圖2 所示，對于δ頻段的腦電信號分解過程如下:

圖2 基于EMD 的ICEEMDAN 分解IMF迭代流程示意圖

①在信號fδ(t)中添加噪聲構(gòu)成信號fδ(t)(i)＝fδ(t)＋ε0E1(w(t)(i))，i＝1，…，N。對每個fδ(t)(i)進(jìn)行EMD 分解得到第一個余量信號:

②計算第一個IMF:

③在第一個余量信號中添加噪聲構(gòu)成信號r1(t)(i)＝r1(t)＋ε1E2[w(t)(i)]，估計第二個余量信號，得到第二個IMF:

④對于k＝3，…，K，在第k－1 個余量信號中添加噪聲構(gòu)成信號rk－1(t)(i)＝rk－1(t)＋εk－1Ek[w(t)(i)]，計算第k個余量信號:

⑤計算第k個IMF:

⑥重復(fù)步驟④～⑤，直到所獲取的余量信號是一個單調(diào)函數(shù)或常量時，分解過程停止。此時，δ頻段的腦電信號fδ(t)最終被分解為:

式中:εk－1為第k次分解中添加噪聲的系數(shù)，通常取值范圍為[0.1，0.3]。

1.3 基于樣本熵的眼電分量識別

樣本熵是一種衡量時間序列復(fù)雜度的方法[23－25]。腦電信號中包含大量的信息，復(fù)雜度高，樣本熵值較高；而眼電信號頻率低，復(fù)雜度低，樣本熵值較低。本文通過設(shè)置樣本熵閾值可實現(xiàn)腦電分量與眼電分量的判別。樣本熵通常表示為SampEn(m，r，N)，其中m表示嵌入維數(shù)，r表示相似容限，N表示數(shù)據(jù)長度，其可以定義為:

式中:Am(r)和Bm(r)分別表示預(yù)處理后的信號與根據(jù)嵌入維度n＋1 和n構(gòu)造成的新序列在相似容限r(nóng)下匹配的概率。本文選取r＝0.2std[x(i)]，(std[x(i)]表示時間序列x(i)的標(biāo)準(zhǔn)差)，m＝2。

計算得到各個IMF 分量的樣本熵值之后，需要設(shè)置閾值才能區(qū)分腦電和眼電偽跡成分。圖3 所示為半模擬數(shù)據(jù)集中54 組FP1 通道純凈腦電信號及眼電信號的樣本熵分布情況。從圖中可以發(fā)現(xiàn)，純凈腦電信號的樣本熵均在0.4 以上，因此將樣本熵閾值設(shè)置為0.4，可以區(qū)分腦電分量和眼電分量。

圖3 純凈腦電信號與眼電信號的樣本熵分布曲線

2 試驗

2.1 數(shù)據(jù)獲取

為了驗證本文提出的EWT-ICEEMDAN 算法的有效性，分別使用半模擬的腦電數(shù)據(jù)和真實腦電數(shù)據(jù)進(jìn)行了評估。

半模擬數(shù)據(jù)集使用的是Klados 和Bamidis 于2016 年公開發(fā)表的數(shù)據(jù)集[26]。該數(shù)據(jù)集記錄了27名健康受試者在閉眼期間的腦電數(shù)據(jù)，包括14 名男性(平均年齡:28.2±7.5 歲)和13 名女性(平均年齡:27.1±5.2 歲)。腦電電極按照國際標(biāo)準(zhǔn)10－20系統(tǒng)在頭皮上放置，信號的采樣頻率為200 Hz，經(jīng)過0.5 Hz～40 Hz 的帶通濾波和50 Hz 的陷波濾波。本文做單通道腦電信號的眼電偽跡去除，選取距離眼睛最近受眼電影響較大的電極Fp1 進(jìn)行眼電偽跡去除。從這27 例受試者中，總共采集54 組腦電數(shù)據(jù)，采集到的腦電數(shù)據(jù)不含眼電偽跡，可以作為純凈的腦電信號。此外，記錄受試者在睜眼期間的眼電信號，分別為垂直眼電(Vertical EOG，VEOG)和水平眼電(Horizontal EOG，HEOG)信號，采樣頻率同樣為200 Hz，經(jīng)過0.5 Hz～5 Hz 的帶通濾波。根據(jù)式(21)生成含眼電偽跡的半模擬腦電信號[27]。

式中:Scon表示含眼電偽跡的腦電信號，SEEG表示純凈的腦電信號，SVEOG和SHEOG分別表示垂直和水平眼電信號，a和b分別表示垂直和水平眼電的污染系數(shù)。

真實數(shù)據(jù)集采集12 名受試者(10 名男性，2 名女性，平均年齡:24.1±1.1 歲)在睜眼情況下FP1 電極的腦電信號，采樣頻率為2 000 Hz，降采樣到200 Hz，再經(jīng)過0.5 Hz～40 Hz 的帶通濾波和50 Hz的陷波濾波，從每名受試者的數(shù)據(jù)中選取10 段含眼電偽跡的腦電信號，每段數(shù)據(jù)段長度為10 s，共得到120 段數(shù)據(jù)。

2.2 效果評價

在本文中分別從時域和頻域兩方面對EWTICEEMDAN 算法的性能進(jìn)行評估。在時域方面，計算純凈的腦電信號SEEG和去除眼電偽跡后的腦電信號Sclean之間的相關(guān)系數(shù)(Correlation Coefficient，CC)和相對均方根誤差(Relative Root-Mean-Squared Error，RRMSE)[28]。CC 的定義為:

RRMSE 的定義為:

式中:Cov(?)表示協(xié)方差，Var(?)表示方差，RMS(?)表示均方根。

相關(guān)系數(shù)評估兩個信號之間的相關(guān)性，CC 值越大，表示二者的相關(guān)性越大，即去除眼電偽跡后的腦電信號信息保留越完整，失真越小。相對均方根誤差評估信號的幅度失真情況，RRMSE 值越小，表示眼電偽跡去除越完全，與純凈的腦電信號越接近。

在頻域方面，使用兩個腦電信號質(zhì)量指標(biāo)來評估所提出算法的性能，分別為含眼電偽跡的腦電信號Scon和去除眼電偽跡后的腦電信號Sclean在δ頻段能量比(Energy Ratio，ER)的變化以及在δ、θ、α、β頻段的功率譜(Power Spectral Density，PSD)的平均絕對誤差(Mean Absolute Error，MAE)[29]。Scon在δ頻段的ER 定義為:

同樣，Sclean在δ頻段的ER 定義為:

進(jìn)一步，計算含眼電偽跡的腦電信號Scon和去除眼電偽跡后的腦電信號Sclean在各個頻段的功率譜的平均絕對誤差，以α頻段為例，的定義為:

對于真實的腦電信號，無法提取純凈的腦電信號，因此不能使用時域指標(biāo)CC 和RRMSE 對眼電偽跡去除效果進(jìn)行評價。在頻域方面，計算眼電偽跡去除前后腦電信號在δ頻段能量比的變化ΔERδ以及在各個頻段的功率譜的平均絕對誤差

為了進(jìn)一步驗證本文算法的有效性，本文對比了DWT、文獻(xiàn)[19]提出的WT-EEMD 以及文獻(xiàn)[29]提出的FBSE-EWT-LPATV 算法。所有的對比實驗包括:

①DWT 算法:首先使用基函數(shù)將腦電信號分解成L 級系數(shù)，然后將高于閾值的每一級的系數(shù)設(shè)為零，最后用逆DWT 重構(gòu)腦電信號。在文獻(xiàn)[14]中研究了四個基函數(shù):haar、coif3、sym3 和bior4.4 以及通用閾值和統(tǒng)計閾值。作者認(rèn)為，基于統(tǒng)計閾值的bior4.4 基函數(shù)是去除眼電偽跡的最佳選擇。

②WT-EEMD 算法:采用EEMD 算法將含眼電偽跡的小波成分分解為若干IMFs，通過設(shè)置自相關(guān)系數(shù)閾值去除含眼電偽跡的IMFs，最后重構(gòu)得到去除眼電偽跡后的腦電信號[19]。

③FBSE-EWT-LPATV 算法:在文獻(xiàn)[29]中首先基于傅里葉－貝塞爾級數(shù)展開的經(jīng)驗小波變換(Fourier-Bessel Series Expansion based Empirical Wavelet Transform，F(xiàn)BSE-EWT)將腦電信號分為δ頻段，再基于局部多項式逼近的全變分算法(Local Polynomial Approximation based Total Variation，LPATV)去除δ頻段的眼電偽跡，最后與其他頻段信號線性重構(gòu)得到去除眼電偽跡的腦電信號[29]。

3 結(jié)果和討論

3.1 半模擬數(shù)據(jù)測試結(jié)果

在本文中，對半模擬數(shù)據(jù)使用本文算法進(jìn)行眼電偽跡去除處理，并與DWT、WT-EEMD 和FBSEEWT-LPATV 三種算法進(jìn)行對比。

圖4 分別顯示了含有眼電偽跡的腦電信號以及使用EWT 提取得到的δ頻段和高頻段信號。在δ頻段信號中，眼電偽跡清晰可見。隨后，采用ICEEMDAN 算法對δ頻段腦電信號進(jìn)行分解。圖5所示為δ頻段腦電信號經(jīng)過ICEEMDAN 分解后的IMFs 結(jié)果及其對應(yīng)的樣本熵值?？梢园l(fā)現(xiàn)IMF4～I(xiàn)MF10 的樣本熵值小于0.4，即為眼電偽跡，去除相應(yīng)的IMF。最后重構(gòu)得到去除眼電偽跡后的腦電信號，如圖6 所示。

圖4 含眼電偽跡的腦電信號及δ 頻段和高頻段信號

圖5 δ 頻段信號進(jìn)行ICEEMDAN 得到IMFs 分量

圖6 本文算法對半模擬數(shù)據(jù)中眼電偽跡去除

在時域方面，計算純凈的腦電信號和去除眼電偽跡的腦電信號的CC 和RRMSE 的均值及標(biāo)準(zhǔn)差，評估不同算法對于相同腦電信號的眼電去除效果及算法的穩(wěn)定性。表1 所示為本文算法、DWT 和WTEEMD 三種算法的CC 和RRMSE 結(jié)果對比。CC 值越大，表示眼電偽跡被去除的同時，腦電信號失真越小?？梢钥闯觯疚乃惴ǖ腃C 均值高于其他兩種算法。RRMSE 值越小，表示眼電偽跡去除的越完全。本文算法的RRMSE 均值和標(biāo)準(zhǔn)差均低于另外兩種算法。因此，本文算法去除眼電偽跡的同時，信號失真最小。在頻域方面，計算眼電偽跡去除前后的腦電信號的ΔERδ以及在δ，θ，α，β頻段的MAE均值及標(biāo)準(zhǔn)差，并與文獻(xiàn)[29]中所提出的FBSEEWT-LPATV 算法結(jié)果對比，如表1 所示?？梢钥闯觯疚乃惴éRδ的均值在50%以上，說明該算法具有較好的去噪效果。對比四種算法，本文所提出的算法在θ，α，β頻段功率譜的MAE 均值基本小于其他算法。因此，去除眼電偽跡后的腦電信號中保留了θ，α，β頻段的信息。

表1 四種算法的眼電偽跡去除效果評價結(jié)果(均值±標(biāo)準(zhǔn)差)

3.2 真實數(shù)據(jù)測試結(jié)果

圖7 所示為本文算法、DWT 和WT-EEMD 三種算法對同一段真實腦電信號處理前后的時域波形對比?？梢钥闯?，本文算法在去除眼電偽跡后，在原眼電偽跡位置處不存在眼電波動，WT-EEMD 算法處理后存在較少波動，而DWT 算法處理后剩余較多波動。

圖7 三種算法對同一段真實腦電數(shù)據(jù)的眼電偽跡去除效果

在本文中，對每名受試者的10 段含眼電偽跡的腦電數(shù)據(jù)都使用本文算法、WT-EEMD 和DWT 算法進(jìn)行眼電偽跡去除處理，并計算δ，θ，α，β頻段功率譜的平均絕對誤差的均值及標(biāo)準(zhǔn)差，用于評價眼電偽跡去除算法對真實腦電數(shù)據(jù)進(jìn)行處理的腦電數(shù)據(jù)失真情況。表2 為本文算法、WT-EEMD 和DWT 算法去除眼電偽跡后的δ頻段的ΔERδ以及在δ，θ，α，β頻段的功率失真情況。很明顯，所有真實腦電信號的ΔERδ均值為71.20%。同樣，θ，α，β頻段功率譜的MAE 均值分別為0.52、0.07 和0.07。可以看出，本文算法在θ，α，β頻段的失真情況均小于WTEEMD 和DWT 算法，且標(biāo)準(zhǔn)差最小，說明該算法對去除腦電信號中的眼電偽跡是有效的。

表2 三種算法的眼電偽跡去除效果評價結(jié)果(均值±標(biāo)準(zhǔn)差)

4 結(jié)論與展望

針對已有眼電偽跡去除算法存在腦電信息丟失的問題，本文提出一種基于經(jīng)驗小波變換－改進(jìn)的自適應(yīng)噪聲完備經(jīng)驗?zāi)B(tài)分解算法(EWT-ICEEMDAN)，利用樣本熵判別腦電和眼電分量，進(jìn)行單通道腦電信號眼電偽跡去除。通過半模擬的腦電信號和真實腦電數(shù)據(jù)進(jìn)行實驗驗證，相比于DWT、WTEEMD 和FBSE-EWT-LPATV 三種算法，本文提出的單通道腦電信號中眼電偽跡去除算法的效果更優(yōu)，既能去除眼電偽跡，又能最大程度地保留腦電信息。在未來工作中，我們將考慮更多能夠有效區(qū)分腦電分量和眼電分量的判別依據(jù)如模糊熵[18]，并且我們將進(jìn)一步研究能夠滿足實時性的方法。