基于CEEMDAN-CFAR 的單通道腦電信號眼電偽跡去除方法研究

荊鈺霏，李川濤，王 偉，于旭東*

（1.海軍軍醫(yī)大學(xué)海軍特色醫(yī)學(xué)中心特種作戰(zhàn)醫(yī)學(xué)研究室，上海 200433；2.海軍軍醫(yī)大學(xué)海軍特色醫(yī)學(xué)中心航空醫(yī)學(xué)研究室，上海 200433）

0 引言

腦電信號由于其非侵入性、靈敏性和準(zhǔn)確性應(yīng)用于腦科學(xué)研究具有一定優(yōu)勢[1]。其中前額無發(fā)區(qū)腦電由于方便采集，在未來可穿戴腦電采集如疲勞監(jiān)測等領(lǐng)域中有很大的應(yīng)用前景[2-3]。然而前額無發(fā)區(qū)距離人眼較近，在采集腦電過程中不可避免地引入較強的眼電偽跡，包含眼動和眨眼產(chǎn)生的電信號，大大地影響了腦電信號的分析處理。一種好的眼電偽跡去除方法有助于提高信號的信噪比和后續(xù)疲勞特征識別的準(zhǔn)確性[4]。

人工剔除法[5]是在臨床腦電圖檢查中較為常用的方法，醫(yī)師們往往選擇舍棄含有大量眼電信號的數(shù)據(jù)，保留較為干凈的腦電信號進行分析。這種方法需要人為剔除眼電信號，工作量大，并且損失了大量有用信號。盲源分離法[6]由于能將混合信號有效分離成包含噪聲和不包含噪聲的獨立成分，成為近年來熱門的偽跡去除方法。但該方法存在分離信號與源信號關(guān)系不確定的問題，往往需要手動挑選剔除眼電成分[7]或是引入其他的算法來進行眼電分量的識別?；陬l域濾波的偽跡去除方法因其簡單有效一直被廣泛使用[8]，比如利用無限沖激響應(yīng)濾波器（infinite impulse responsefilter，IIR）和有限沖激響應(yīng)濾波器（finite impulse responsefilter，F(xiàn)IR）去除偽跡，然而使用IIR 濾除眼電，眼電尖峰會造成較長時間的振蕩而引入新的噪聲，使用FIR 濾除眼電，需要較長的階數(shù)以獲得較高的噪聲抑制能力?；跁r頻域的小波變換[9]因其良好的時頻特性成為最常被使用的偽跡去除方法之一，需要注意的是，眼電偽跡和腦電信號有頻帶重疊，采用頻域濾波的方式在濾除眼電信號的同時也損失了很大一部分腦電信號，因此研究者們在傳統(tǒng)小波變換的基礎(chǔ)上又聯(lián)合使用了其他判別方法濾除眼電，比如李洋[9]采用小波變換結(jié)合卡爾曼濾波器的方法去除眼電偽跡，周衛(wèi)東等[10]采用小波軟門限法去除腦電中的噪聲。

與小波分解相比，經(jīng)驗?zāi)B(tài)分解（empirical mode decomposition，EMD）吸取了小波變換多分辨力的優(yōu)勢，卻沒有基函數(shù)的選擇問題。EMD 在1998 年由Huang 等[11]提出，因其良好的自適應(yīng)性和短數(shù)據(jù)序列處理能力被廣泛應(yīng)用于非線性、非平穩(wěn)數(shù)據(jù)的分析。傳統(tǒng)EMD 算法在信號不屬于白噪聲時，在分解中會因為時間尺度丟失而造成模態(tài)混疊[11]。Torres 等[12]提出的基于自適應(yīng)噪聲完備集成經(jīng)驗?zāi)B(tài)分解（complete ensemble empirical mode decomposition with adaptive noise，CEEMDAN）方法吸收了集成經(jīng)驗?zāi)B(tài)分解（ensemble empirical mode decomposition，EEMD）的優(yōu)點，能有效解決模態(tài)混疊問題，并通過逐層添加自適應(yīng)噪聲的分解方法減少了整體的噪聲引入，被認(rèn)為是一種性能良好的信號分解方法。恒虛警率（constant false alarm rate，CFAR）算法是一種雷達動目標(biāo)檢測方法，能在不同的環(huán)境噪聲背景下保持靈敏的動目標(biāo)檢測能力。

基于此，本研究提出了一種在CEEMDAN 的基礎(chǔ)上使用CFAR 算法找到眼電尖峰并去除眼電偽跡的預(yù)處理方法，即CEEMDAN-CFAR。首先采用CEEMDAN 的方法分解原始信號，分離出低頻噪聲信號和含有明顯眼電脈沖信號的本征模態(tài)函數(shù)（intrinsic mode function，IMF）分量。然后再采用CFAR檢測算法識別和去除眼電偽跡，重構(gòu)得到干凈的腦電信號。

1 方法

1.1 CEEMDAN

EMD 是將信號分解成一組IMF 的一種信號分解方法。IMF 為滿足以下2 個條件的函數(shù)：

（1）整個時間歷程內(nèi)，穿越零點的次數(shù)與極值點數(shù)相等或至多相差1。

（2）信號上任意一點，由局部極大值定義的上包絡(luò)線和由局部極小值定義的下包絡(luò)線的均值為0，即信號關(guān)于時間軸局部對稱。

通過將信號不斷分解成IMF 分量的方式，EMD將原始信號的各尺度分量不斷從高頻到低頻進行提取。大量實驗證明[11]，EMD 方法類似于一個自適應(yīng)二進濾波器，能將白噪聲分解為一系列中心頻率呈二進制關(guān)系的IMF[13]。但是當(dāng)數(shù)據(jù)不是純的白噪聲時，分解中的時間尺度丟失會導(dǎo)致模態(tài)混疊。EEMD 就是為了解決這種模態(tài)混疊而產(chǎn)生的[14]。EEMD 在原始信號中加入白噪聲以補充原始信號分解中時間尺度不足的問題，并利用白噪聲統(tǒng)計特性，通過多次加入噪聲求得特征模態(tài)函數(shù)均值來抵消噪聲。EEMD固然在一定程度上解決了模態(tài)混疊的問題，然而在重構(gòu)信號時也引入了新的噪聲。同時，不同的噪聲添加方式會帶來不同的EEMD 分解層數(shù)。CEEMDAN通過EEMD 求得高頻分量，將剩余的所有分量看成新的信號加入新的自適應(yīng)噪聲再求得高頻分量，不斷分解到最后一層[12]。這種信號分解方式在降低計算量的同時提供了更精確的重構(gòu)信號，本文使用CEEMDAN 的方式分解腦電信號，可以有效獲得眼動分量和低頻噪聲分量。假設(shè)有原始腦電信號X（t），X（t）的CEEMDAN 具體實現(xiàn)如下：

（1）輸入原始腦電信號X（t），令j=1（j 為循環(huán)計數(shù)變量）。

（2）判斷X 的極值個數(shù)是否大于2，若大于2 則繼續(xù)執(zhí)行（3），否則轉(zhuǎn)到（7）。

（3）令i=1（i 為循環(huán)計數(shù)變量）。

（4）判斷i 是否小于N+1（N 為自定義的重復(fù)添加噪聲次數(shù)），若小于N+1 則執(zhí)行（5），否則轉(zhuǎn)到（6）。

（5）在原始腦電信號X 中加入白噪聲，進行EMD得到第一層IMF1，令i=i+1，轉(zhuǎn)到（4）。

（6）退出循環(huán)，計算N 個IMF1 的均值得到IMF（j），得到新的待處理腦電信號X=X-IMF（j），令j=j+1，轉(zhuǎn)到（2）。

（7）退出循環(huán)得到全部IMF 分量。

算法流程圖如圖1 所示。

圖1 CEEMDAN 算法流程圖

1.2 基于CFAR 算法確定眼動位置

原始腦電信號經(jīng)由CEEMDAN 方法分解后得到多層IMF 分量，眼電信號在IMF 中表現(xiàn)為振幅突然增大的、類似脈沖的信號。CFAR 算法過去一般被用于雷達動目標(biāo)檢測[15]，本文參考Lazaro 等[16]提到的關(guān)于使用方差作為識別單元檢測突變信號的方法，將原始信號序列X 以一定窗寬和步長求得各單元方差生成新的序列，如圖2 所示。然后通過檢測異常方差值來確定眼電信號的位置。將生成的方差序列看成新的信號序列使用CFAR 算法檢測異常值。在該算法中保護窗和訓(xùn)練窗長度取值皆為1，CFAR 檢測原理如圖3 所示。

圖2 CFAR 算法識別單元序列生成方法

圖3 CFAR 檢測原理圖

根據(jù)經(jīng)驗將CFAR 的閾值T 設(shè)為10，識別結(jié)果H 為

當(dāng)H 大于0 時為H1，表示有眼電；當(dāng)H 小于0 時為H0，表示沒有眼電。

2 實驗過程及結(jié)果

2.1 實驗過程

2.1.1 實驗條件

本研究以疲勞實驗過程中采集到的睜眼腦電信號作為原始信號進行預(yù)處理實驗。疲勞實驗采集被試前額FP1 處腦電信號，以人體左右乳突分別作為參考和地端。腦電采集使用課題組自己研制的雙通道前額柔性腦電采集設(shè)備，系統(tǒng)采樣頻率250 Hz，模數(shù)轉(zhuǎn)換數(shù)據(jù)精度為24 bit，短路噪聲小于3 μV，10 μV～100 mV 信號測量誤差不大于10%，輸入阻抗680 MΩ。使用該設(shè)備采集一名34 歲健康男性無任務(wù)狀態(tài)下的睜眼腦電數(shù)據(jù)3 min。

2.1.2 CEEMDAN-CFAR 算法信號處理過程

為具體介紹CEEMDAN-CFAR 方法的信號處理過程，在健康成年人的3 min 睜眼腦電數(shù)據(jù)中挑選一段10 s 的腦電信號并以此為例進行預(yù)處理。10 s含有眼電偽跡的腦電信號如圖4 所示，由圖可知該信號包含4 個明顯的眨眼眼電和低頻眼動眼電，是具有代表性的復(fù)雜腦電信號。首先，采用CEEMDAN 算法將原始腦電信號分解成多層IMF，CEEMDAN 算法中添加白噪聲權(quán)重為0.2，重復(fù)添加噪聲次數(shù)為20，將EMD 內(nèi)部最大包絡(luò)次數(shù)設(shè)定為100。CEEMDAN算法分解示意圖如圖5 所示。

圖4 含有眼電偽跡的腦電信號

圖5 CEEMDAN 算法分解示意圖

由1.1 章節(jié)中CEEMDAN 算法的原理介紹可知，該算法具有類似小波的二進濾波器特性，第n 層的中心頻率可由公式（2）計算得出：

式中，f（n）為第n 層的中心頻率，n 為IMF 所在層數(shù)；fs為系統(tǒng)采樣頻率。為了直觀地看到各層IMF 分量的頻率分布，本研究采用功率譜估計算法Lomb[17]對各層IMF 分量進行頻域轉(zhuǎn)換，結(jié)果如圖6 所示。

由圖6 可以明顯看到，各層中心頻率逐層降低。正常人的眨眼時間為0.2～0.4 s[18]，眼電信號在頻帶0～10 Hz 之間，根據(jù)公式（2）計算可得IMF4、IMF5、IMF6、IMF7 的中心頻率分別為15.62、7.81、3.91、1.95 Hz，本文依據(jù)多段信號CEEMDAN 分解后觀察到的眼電分布特征以及中心頻率同眼電頻率接近的IMF5、IMF6、IMF7 作為需要去除眼電的通道。

圖6 各IMF 分量Lomb 轉(zhuǎn)換圖

依據(jù)本文算法流程，使用CFAR 檢測法找到眼電尖峰所在的位置。將原始信號序列X 以fs/2 為窗寬、以fs/4 為步長求得各單元方差后生成的新序列Y（n），即圖7 中紅色曲線（采樣頻率fs為250 Hz）。CFAR的閾值取值為10，使用公式（1）進行方差檢測，方差檢測結(jié)果異常標(biāo)記為1（即眼電所在位置），其余標(biāo)記為0，如圖7 所示。圖7（a）、（b）、（c）分別為IMF5、IMF6、IMF7 分量的CFAR 檢測結(jié)果。圖7 表明，CFAR算法能準(zhǔn)確地識別出眼電尖峰所在的位置。

圖7 CFAR 檢測結(jié)果圖

確定眼電尖峰位置后，首先將該層分量中以尖峰位置為中心的小段信號置零，置零信號的長度計算方法為

式中，L 為清零信號段的長度。假設(shè)落在第n 層IMF 上的眼電信號頻率正好為中心頻率，則眼電信號長度為中心頻率信號周期長度，即fs/f（n），又由于在EMD 的過程中使用三次樣條插值的方法來構(gòu)造信號的上、下包絡(luò)線，眼電脈沖信號在分解的過程中會引起周圍信號變形，所以計算置零信號長度為2 倍的中心頻率信號周期長度。

最后，舍棄低頻分量，即把IMF8 及其以后的分量置零，將其余分量重構(gòu)成不含偽跡的純凈腦電信號。整體的信號處理流程如圖8 所示。

圖8 基于CEEMDAN-CFAR 算法的眼電偽跡去除流程圖

2.2 實驗結(jié)果

以隨機挑選的健康成年人10 s 放松狀態(tài)下的睜眼腦電數(shù)據(jù)為例進行預(yù)處理，含有眼電偽跡的腦電信號（如圖4 所示）經(jīng)CEEMDAN-CFAR 處理后的信號如圖9 所示。由圖9 可以直觀地看到眼電偽跡全部被去除。將原始信號和處理后的信號局部放大后進行比較，如圖10 所示。由圖10 可以看出，預(yù)處理前后信號一致性較強，證實了使用該預(yù)處理方法能在去除偽跡的同時很好地保留局部細節(jié)和有用信息。

圖9 CEEMDAN-CFAR 算法處理后的信號

圖10 原始信號和處理后的信號局部對比圖

為進一步驗證預(yù)處理算法的有效性，使用該方法處理完整的3 min 睜眼腦電數(shù)據(jù)。預(yù)處理前后對比如圖11 所示，具體細節(jié)如圖12 所示。由圖11、12可以看出，本文提出的預(yù)處理方法針對真實的腦電數(shù)據(jù)具有良好的眼電噪聲濾除效果。

圖11 3 min 腦電數(shù)據(jù)預(yù)處理前后效果圖

圖12 3 min 腦電數(shù)據(jù)預(yù)處理細節(jié)圖

2.3 算法性能

使用Klados 等[19]公開發(fā)表的一個半模擬腦電數(shù)據(jù)集對算法的有效性進行驗證。在純凈的腦電信號中人為地加入眼電偽跡模擬真實的腦電信號，如圖13所示。采集27 名健康受試者閉眼期間的腦電數(shù)據(jù)，信號的采樣頻率為200 Hz，采用0.5～40 Hz 的帶通濾波和50 Hz 的陷波濾波，共獲得54 個數(shù)據(jù)集。每個數(shù)據(jù)集都包含連續(xù)30 s 的腦電信號，對獲得的數(shù)據(jù)集進行了仔細檢查，確保沒有顯著污染，可以作為純凈的腦電信號。在相同情況下采集27 名健康受試者的眼電圖（electro-oculogram，EOG）數(shù)據(jù)，EOG 信號在0.5～5Hz下進行帶通濾波。根據(jù)公式（4）、（5）模擬生成含有眼電偽跡的半模擬腦電信號：

圖13 半模擬數(shù)據(jù)集中的純凈腦電信號和污染腦電信號

式中，Scon表示污染的腦電信號；Spure表示純凈的腦電信號；SEOG表示EOG信號；SVEOG表示垂直EOG信號；SHEOG表示水平EOG信號；α 和β 分別表示垂直EOG 和水平EOG 的污染系數(shù)。

從54 個數(shù)據(jù)集中隨機選擇12 個數(shù)據(jù)集中的污染腦電信號（如圖14 所示），使用CEEMDAN-CFAR算法對FP1 通道的污染的腦電信號進行處理，處理后得到較為純凈的腦電信號，如圖15 所示。

圖14 原始的半模擬腦電信號

圖15 CEEMDAN-CFAR 算法處理后的半模擬腦電信號

將組成半模擬腦電信號的純凈腦電信號進行CEEMDAN，去除其中低頻成分后再與處理后的半模擬腦電信號進行相關(guān)性分析，得到相關(guān)系數(shù)見表1。由表1 中數(shù)據(jù)可知，經(jīng)過預(yù)處理的信號與純凈腦電信號的相關(guān)性得到了很大提升。與樣本熵-完備經(jīng)驗?zāi)B(tài)分解（sample entropy-complete ensemble empirical mode decomposition with adaptive noise，SE-CEEMDAN）[20]算法的相關(guān)系數(shù)對比：CEEMDAN-CFAR 的相關(guān)系數(shù)（均值）為0.797，方差為0.062，SE-CEEMDAN 的相關(guān)系數(shù)（均值）為0.841，方差為0.005。處理200 Hz采樣頻率下5 s 數(shù)據(jù)所用時間的對比：SE-CEEMDAN運行時間為68 s，CEEMDAN-CFAR 運行時間為12 s。

表1 純凈腦電信號與處理前、后腦電信號的相關(guān)系數(shù)

綜上，CEEMDAN-CFAR 算法在單通道腦電信號預(yù)處理性能上有較高的準(zhǔn)確性，雖然與同類算法SE-CEEMDAN[20]相比保留原始信號的能力較弱，但速度較快，是一種方便快捷的預(yù)處理方法。

3 結(jié)語

腦電信號作為一種微弱的生物電信號，其波幅一般都在μV 的數(shù)量級上，特別容易受到外界干擾。在針對腦電的信號處理中，各種偽跡的存在嚴(yán)重影響了信號特征提取的準(zhǔn)確性。因此，有效的偽跡去除方法是整個腦電信號處理流程中至關(guān)重要的一環(huán)。其中眼電偽跡由于其幅度大、隨機性強、頻率與腦電信號存在頻帶重疊等特點，成為腦電信號偽跡去除中的重點與難點。傳統(tǒng)的眼電偽跡去除方法存在丟失有用信號或引入冗余噪聲等問題，而近年來較為熱門的盲源分離算法大多是針對多通道腦電數(shù)據(jù)而言的。在單通道腦電偽跡去除方法的研究上，雖然也有研究者[21]提出了將單通道腦電信號通過希爾伯特-黃變換等方法轉(zhuǎn)化為多通道數(shù)據(jù)后，再使用獨立成分分析（independent component analysis，ICA）分離眼電偽跡的方法，但從ICA 分離得到的各獨立源信號中判斷眼電成分又需要加入新的算法，這大大增加了信號處理的難度。在此基礎(chǔ)上，本文提出基于CEEMDAN-CFAR算法眼電偽跡去除方法，通過CEEMDAN 算法對原始腦電信號進行分解，得到低頻噪聲成分和含有明顯眼電脈沖信號的IMF 分量，然后采用CFAR 算法準(zhǔn)確定位眼電位置并予以剔除。該方法具有計算簡單、眼電偽跡去除效果好的優(yōu)勢，與同類算法SECEEMDAN 相比，CEEMDAN-CFAR 算法保留原始信號的能力稍弱，但速度快捷。不足之處在于該方法對于低頻信號的保留度較低，不適用于研究深度睡眠等需要關(guān)注低頻腦電信號的領(lǐng)域。在下一步的研究中，將進一步對CEEMDAN 分層后相對低頻的部分進行處理，通過各低頻IMF 分量的頻域特征判別是否需要保留該分量，以此獲得更為純凈的腦電信號。