一種基于EWT-ICEEMDAN 的單通道腦電信號眼電偽跡去除算法?

宋 婷，舒智林，孫玉波，韓建達，于寧波?

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

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

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

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

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

1 方法

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

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

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

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

首先，對于給定信號f(t)通過傅里葉變換得到傅里葉頻譜，假設將傅里葉支撐區(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)的表達式為[22]:

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

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

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

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

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

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

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

②計算第一個IMF:

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

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

⑤計算第k個IMF:

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

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

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

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

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

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

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

2 試驗

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

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

半模擬數(shù)據(jù)集使用的是Klados 和Bamidis 于2016 年公開發(fā)表的數(shù)據(jù)集[26]。該數(shù)據(jù)集記錄了27名健康受試者在閉眼期間的腦電數(shù)據(jù)，包括14 名男性(平均年齡:28.2±7.5 歲)和13 名女性(平均年齡:27.1±5.2 歲)。腦電電極按照國際標準10－20系統(tǒng)在頭皮上放置，信號的采樣頻率為200 Hz，經(jīng)過0.5 Hz～40 Hz 的帶通濾波和50 Hz 的陷波濾波。本文做單通道腦電信號的眼電偽跡去除，選取距離眼睛最近受眼電影響較大的電極Fp1 進行眼電偽跡去除。從這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 算法的性能進行評估。在時域方面，計算純凈的腦電信號SEEG和去除眼電偽跡后的腦電信號Sclean之間的相關系數(shù)(Correlation Coefficient，CC)和相對均方根誤差(Relative Root-Mean-Squared Error，RRMSE)[28]。CC 的定義為:

RRMSE 的定義為:

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

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

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

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

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

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

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

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

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

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

3 結果和討論

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

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

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

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

圖5 δ 頻段信號進行ICEEMDAN 得到IMFs 分量

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

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

表1 四種算法的眼電偽跡去除效果評價結果(均值±標準差)

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

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

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

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

表2 三種算法的眼電偽跡去除效果評價結果(均值±標準差)

4 結論與展望

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