單通道腦電信號中眼電干擾的自動分離方法

吳明權(quán) 李海峰 馬 琳

1 引言

近年來，非侵入式腦機接口(Brain-Computer Interface, BCI)技術(shù)深入發(fā)展，為人與機器、人與人提供了新的交流方式。目前已經(jīng)可以利用腦電信號(ElectroEncephaloGram, EEG)來實現(xiàn) 3維空間的直升飛機飛行控制[1]。在美國軍方“Silent Talk”項目中，通過檢測分析特殊詞的神經(jīng)信號實現(xiàn)了用戶之間的腦-腦直接通信[2]。干電極腦電采集技術(shù)的日益成熟，加快了便攜式BCI產(chǎn)品的商品化，使其深入到人們學(xué)習(xí)、娛樂、生活的方方面面。然而，便攜式BCI產(chǎn)品通常是少通道甚至是單通道的，所獲取的 EEG信號極易受到肌電、心電尤其是眼電(ElectroOculo Gram, EOG)的干擾，導(dǎo)致真實腦電特征難以提取。對各種偽跡特別是眼電的剔除方法研究，始終是腦電信號處理領(lǐng)域的難點問題。

目前主要的眼電偽跡干擾去除技術(shù)有兩類：偽跡排除和偽跡校正。前者是將包含眼電干擾的腦電時段簡單地排除；后者則采用各種方法來估計眼電信號，并在原始EEG中減去估計的偽跡成分，以獲得純凈的EEG。當(dāng)前偽跡校正的主流方法有平均偽跡回歸分析方法[3]和基于獨立成分分析(Independent Component Analysis, ICA)的方法[4,5]。平均偽跡回歸分析方法假設(shè)眼電電極和各頭皮電極之間的傳導(dǎo)系數(shù)不變，利用眼電通道與其他多個通道的相關(guān)性來估計傳導(dǎo)系數(shù)，并從各個通道中按傳導(dǎo)系數(shù)減去眼電信號而獲得正常的 EEG信號。ICA則恰恰相反，它假設(shè)多個信號源相互獨立，并通過機器學(xué)習(xí)的方法來最大化源信號間的獨立性，分離出眼電成分并予以剔除，再利用其余的成分重構(gòu)EEG信號。平均偽跡回歸分析方法需要獨立的眼電通道，而ICA方法則要求不同通道在空間具有一定的分布，且通道數(shù)要大于信號源的數(shù)目。可見，以上兩種主流方法在便攜式BCI中均難以適用，如何自動高效地去除眼電干擾已成為少通道/單通道腦電信號處理的難點問題。

小波分析可以對信號進行多尺度細化分析，具有良好的時域、頻域分辨率和適應(yīng)性。近年來小波變換在肌電、神經(jīng)動作電位、腦電等生理信號的壓縮[6,7]、去噪[8,9]、特征提取[10]，分類[11]和預(yù)測[12]等方面廣泛應(yīng)用。由于正常腦電頻率較高、幅值較小、能量譜相對分散，而眼電信號卻具有低頻、高能量、時間有限和能量相對集中等特點，分解后的小波系數(shù)稀疏，合理選擇重構(gòu)閾值，就能實現(xiàn)單通道腦電信號的眼電分離。因此，本文提出了一種在眼電初步定位后，利用小波變換實現(xiàn)眼電分離的方法。

2 基于小波變換分離眼電方法的原理

眼電污染的單通道腦電信號可簡單表示為

X[n]是一段受眼電干擾的單通道腦電信號，S[n]是正常腦電信號，A[n]是眼電信號。由于眼電信號是低頻信號(<5 Hz)[13]，具有比正常腦電信號更高的幅值和更光滑的波形，因此，眼電信號提取問題可以認為是尋找一個與 X[n]低頻特性相似的光滑函數(shù)A'[n]作為A[n]估計，并使 S '[n ] = X [ n ] - A '[n]的平均方差盡可能小[14]。由于正交多分辨率分析中的小波基是 Banach空間中的無條件正交基，所以小波系數(shù)的衰減會使重構(gòu)函數(shù)的模增大為原來的常數(shù)倍，其中細節(jié)小波系數(shù)的衰減會使重構(gòu)函數(shù)比原函數(shù)更為平滑。鑒于眼電信號小波能量譜相對集中于低頻區(qū)域，因此合理地衰減眼電部分細節(jié)小波系數(shù)就能估計出眼電 A '[n]。

3 基于長時差分振幅包絡(luò)的雙門限眼電定位方法

相對于正常腦電，眼電污染的單通道腦電幅值較高，具有近似正弦波形(如圖 1(a)所示)，難以使用簡單閾值方法直接確定眼電區(qū)間。雖然眼電偽跡波動較大，但由于高頻腦電的影響，使用一階差分確定眨眼區(qū)間也難以適用。本文綜合閾值方法和差分方法兩者的優(yōu)勢，采用腦電信號的長時差分作為眼電波動程度的度量，然后利用低通濾波獲取信號的振幅包絡(luò)，最后使用雙門限端點檢測方法確定眼電干擾的區(qū)間。

3.1 長時差分信號的計算

式(1)所表示單通道腦電信號的長時差分 Dk[n]可表示為

k是時間延遲，{S[ n ] - S [ n - k ]}為腦電成分的長時差分，{A[n ] - A [n - k ]}為眼電成分的長時差分。由于正常腦電具有短時相關(guān)性，且呈現(xiàn)近似 N ( 0,σ)正態(tài)分布，所以只要k大于普通腦電的相關(guān)時間，則可認為 S [ n]與 S [ n - k ]不相關(guān)，{S[ n ] - S [ n - k ]}服從N ( 0,2σ)分布。由于眼電具有緩變信號特點，當(dāng) k較小時，眼電的差分{A[n ] - A [n - k ]}仍具有時變性。所以合理選擇k可以做到既能獲得較大的眼電長時差分，又能有效地去除正常腦電的相關(guān)性。通常自然眨眼活動時間在300~500 ms之間，而普通腦電在50 ms以后相關(guān)性已經(jīng)很不明顯。本文在512 Hz采樣率下，選取k=80，相當(dāng)于半個眨眼周期(約160 ms)，以獲得較大的眼電差異。提取的單通道長時差分信號如圖1(b)中所示。另外，長時差分還有去除由放大器引入的直流成分和緩慢電位漂移等作用，更有利于眼電區(qū)間的檢測。

3.2 長時差分信號振幅包絡(luò)的提取

振幅包絡(luò)提取過程圖2所示。首先，運用平方運算計算長時差分能量信號，并進行2倍幅值提升，以補償平方運算造成的低頻能量損失；然后，通過M階降采樣，以減小后續(xù)低通濾波器的階數(shù)，提升運算速度。其次，能量信號通過N階具有線性相位延遲的有限沖擊響應(yīng)低通濾波器 h[n]，獲得能量信號包絡(luò)；再實施平移(N-1)/2個單位以補償?shù)屯V波時間延遲，做到與原信號的時間對齊。最后運用開方運算，將能量包絡(luò)轉(zhuǎn)換為振幅包絡(luò) E[n]。提取的長時差分信號包絡(luò)如圖1(b)中所示。

3.3 雙門限眼電干擾區(qū)間端點檢測方法

圖1 單通道腦電眼電分離效果圖

圖2 振幅包絡(luò)提取流程示意圖

由于眼電部分的振幅包絡(luò)具有較高幅值，所以可以首先設(shè)置較高的閾值Th，排除高幅值的正常腦電，準(zhǔn)確定位到眼電的位置，然后設(shè)置較低的閾值Tl，通過前后向搜索小于閾值Tl的振幅位置，以精確確定眼電干擾區(qū)間。由于正常腦電長時差分信號的振幅包絡(luò)幅值正比于正常腦電的幅值，因此，可通過統(tǒng)計腦電的標(biāo)準(zhǔn)差獲得合適的閾值為

σ是正常腦電振幅的標(biāo)準(zhǔn)差， λh, λl為經(jīng)驗性常數(shù)。由于正常腦電長時差分信號近似服從N ( 0,σ)正態(tài)分布，因此，一般可取λh≥ 3以越過高幅值的正常腦電包絡(luò)，而取λl≥ 2確定眼電起止點。本文選取 λh=5, λl=3來確定眼電區(qū)間。

本文眼電區(qū)間端點檢測方法的時延T是影響整個眼電分離方法實時性的關(guān)鍵因素，主要來自低通濾波的時延，可以由式(5)計算得到

M為降采樣因子，N為濾波器階數(shù)，本文取M=8,N=21，即最遲約160 ms(80個采樣點，采樣率512 Hz)后就能檢測到眼電信號，具有較強的實時性。

4 基于小波變換的眼電估計與分離方法

眼電估計和分離過程如圖3所示。首先，對包含眼電的原始腦電信號進行Mallat小波分解。之后，根據(jù)各層系數(shù)分布自適應(yīng)調(diào)整小波系數(shù)。使用新系數(shù)重構(gòu)眼電估計信號后，從原始腦電信號中減去眼電估計信號，實現(xiàn)眼電的分離。最后使用中值濾波進行端點修正，分離出純凈連續(xù)的腦電信號。

圖3 基于小波變換的腦電眼電分離過程示意圖

應(yīng)用小波分解的關(guān)鍵在于選擇合適的小波基函數(shù)和確定小波分解的層數(shù)。由于人眼球可被視為一個角膜端為正極、視網(wǎng)膜端為負極的雙極性球體，眼球轉(zhuǎn)動和眨眼時，會改變眼球附近的電位分布形成復(fù)雜的眼電信號。典型眨眼眼電波形表現(xiàn)為持續(xù)大約 200~400 ms(<5 Hz)的單相偏移(記錄方式不同也可表現(xiàn)出雙相偏移)[15]。為了獲得眼電信號的光滑近似，本文選用與其具有較高相似性，且具有良好的對稱性和光滑性的 sym5小波基函數(shù)，以使分解稀疏，并降低相位損失。小波分解層數(shù)過多，經(jīng)閾值處理后會損失較多的局部眼電信息；反之分解層數(shù)過少，重構(gòu)眼電中會混入過多的腦電信號。本文采取6層小波分解，使頂層的細節(jié)小波系數(shù)對應(yīng)于0~16 Hz，約為典型眼電波形頻率范圍的3倍，有效地平衡對小波分解層數(shù)與信號細節(jié)的不同要求。此外，為減小小波變換端點效應(yīng)的影響，本文選取[N1-80, N2+80]范圍內(nèi)的信號段進行小波分析。

4.1 腦電信號的Mallat小波分解

Mallat方法是正交多分辨率小波變換的高效計算方法。首先，將包含眼電的原始腦電信號作為最底層，即c[0, k]=X[k]。然后，自底向上逐層分解本層近似系數(shù)c[j,k](j為層數(shù)06j≤≤; k為該層第k個小波系數(shù))，得到上一層近似系數(shù) c[j+1,k]和細節(jié)系數(shù)d[j+1,k]。其分解式(6)可表示為c[j, k]分別通過h[n], g[n]的互為正交共軛鏡像的低、高通濾波器，并施行因子為2的降采樣處理。

4.2 基于Birgé.-.Massart策略的自適應(yīng)小波系數(shù)閾值確定方法

準(zhǔn)確提取小波系數(shù)是影響眼電重構(gòu)質(zhì)量的關(guān)鍵。通常采取在全部保留最上層的近似系數(shù)基礎(chǔ)上，每層保留模大于規(guī)定閾值的細節(jié)系數(shù)，而將模小于規(guī)定閾值的細節(jié)系數(shù)置零的方法來選擇重構(gòu)小波系數(shù)。本文使用Birgé.-.Massart策略自適應(yīng)地確定小波系數(shù)閾值：

(1)保留最高層J=6的全部近似系數(shù)；

(2)確定第j層(1)jJ≤≤保留系數(shù)的個數(shù)為

L通常為 LJ≤ L ≤ 2 LJ，其中LJ為最高層近似系數(shù)的個數(shù)，2≤α≤3，本文取α=2；

(3)根據(jù)第j層的保留個數(shù)nj，設(shè)定第j層的閾值為第nj大的系數(shù)模|d[j,k]|。

4.3 眼電估計的小波重構(gòu)和與EEG分離

使用Mallat重構(gòu)算法重構(gòu)眼電估計信號，用式(8)自頂向下逐層得到下一層的近似系數(shù)，直到第 0層即為重構(gòu)眼電信號，A'[k]=c[0,k]。

4.4 基于中值濾波的端點校正方法

原始腦電去除眼電后，可能會在端點處形成間斷點，故本文采用5點的中值濾波來補償端點效應(yīng)，具體做法如式(9)所示。med表示取中值操作。

應(yīng)用小波變換眼電分離后的眼電信號和腦電信號如圖1(d)和圖1(e)所示。

5 實驗結(jié)果及分析評價

5.1 評測本文方法所使用的數(shù)據(jù)

為了評測所提出眼電分離方法的效果，本文使用了大量公開腦電數(shù)據(jù)和本單位認知實驗的腦電數(shù)據(jù)。按照測試目標(biāo)，可以歸納3類：數(shù)據(jù)集1為單通道 BCI實驗數(shù)據(jù)集，共計 3 h，使用 NeuroSky的一款發(fā)帶式單通道腦電設(shè)備MindBand采集。參考電極和接地電極都連接在耳掛上，采集腦電的干電極傳感器置于前額相當(dāng)于FP1位置，AD轉(zhuǎn)換分辨率為12 bit，采樣率為512 Hz。數(shù)據(jù)集2為本單位心理實驗數(shù)據(jù)集，共計3 h，使用Neuroscan公司SynAmps2系統(tǒng)采集，包括雙極記錄的 VEOG和HEOG的66通道腦電數(shù)據(jù)，AD轉(zhuǎn)換分辨率為24 bit，采樣率為500 Hz。數(shù)據(jù)集3是公開的聽、視覺注意轉(zhuǎn)換實驗數(shù)據(jù)集，包括單極記錄的左、右EOG的36通道腦電數(shù)據(jù)[16]，共計28 h，采樣率為250 Hz。其中數(shù)據(jù)集1為本文方法應(yīng)用的主要目標(biāo)場合，數(shù)據(jù)集2, 3為與平均偽跡回歸分析方法和ICA方法的比較數(shù)據(jù)，其長時差分時間延遲k根據(jù)實際采樣率來相應(yīng)設(shè)置。

5.2 量化評價指標(biāo)

眼電分離效果具有鮮明的直觀性，穩(wěn)健的眼電分離方法不僅能夠去除高幅值的眨眼偽跡，還能保證信號在細節(jié)上與原信號高度相似。本文選用相關(guān)系數(shù)(R)作為提取的眼電估計信號的量化評價指標(biāo)，選用信噪比(SNR)作為分離眼電后的腦電信號的量化評價指標(biāo)，其定義為

式(10)中的X1, X2為兩個待比較信號，在本文中對應(yīng)包含眼電的原始腦電時段X[n]和分離后的眼電時段A'[n]。式(11)中的S1, S2分別代表信號和包含噪聲的信號，在本文中分別對應(yīng)原始腦電信號X[n]和分離后的純凈腦電信號S'[n]。

5.3 實驗結(jié)果及分析

5.3.1 數(shù)據(jù)集1上眼電分離效果 圖1是數(shù)據(jù)集1上的眼電分離效果圖，可以看出提取的眼電與原信號趨勢基本一致，分離后的腦電與原信號細節(jié)非常相似，基本達到了眼電分離的要求。表1給出了數(shù)據(jù)集1上本文的小波變換眼電去除方法(*-WT)與檢測眼電區(qū)間后使用偽跡排除方法(*-OR)分離的腦電信噪比和眼電的檢出率。可以得出以下結(jié)論：本文方法具有相當(dāng)高的眼電檢出率，分離的眼電信號與包含眼電的原始腦電段高度相關(guān)，具有很好的相似性，分離的腦電信號較OR方法具有更高的信噪比，有效減小了信號失真。

表1 數(shù)據(jù)集1的眼電信號分離方法結(jié)果統(tǒng)計

5.3.2 數(shù)據(jù)2上與主流眼電分離方法效果比較 在數(shù)據(jù)集2上選擇受眼電影響強、中和弱的3個電極FP1,CZ和OZ，用于本文方法與普遍采用的平均偽跡回歸分析方法(*-RA)和 ICA方法(*-ICA)進行眼電分離效果比較。眼電分離效果如圖4所示，虛線左側(cè)為眼電偽跡去除對比，右側(cè)為其他較大的偽跡去除對比。其中，F(xiàn)P1, CZ和OZ為原始信號，VEOG為垂直眼電信號?？梢钥闯?，對于眼電偽跡，平均偽跡回歸分析方法和本文方法都能有效地去除，而ICA方法則引入了大幅值的高頻噪聲，在沒有單電極記錄的EOG數(shù)據(jù)上ICA方法不能有效地分離和去除眼電干擾。另外，在FP1位置上，平均偽跡回歸分析方法不僅去除了眼電，而且去除了低頻的線性腦電成分。對比3種方法的眼電分離效果，本文方法無論是對受眼電影響較大或較小的電極，都能很好地去除眼電。對于其他較大的偽跡，與另外兩種方法相比，本文方法不僅能檢測出偽跡，而且還能有效地去除。

表 2給出了本文方法與平均偽跡回歸分析和ICA方法眼電分離效果的量化評價，可以得出如下結(jié)論：本文方法與平均偽跡回歸分析方法相比，在FP1電極上眼電去除效果明顯較好，而在CZ和OZ電極上眼電去除效果相當(dāng)；與 ICA方法比較則在CZ, OZ效果明顯較好，在FP1上效果相當(dāng)。因此，本文方法在總體上均好于平均偽跡回歸分析方法和ICA方法。

圖4 數(shù)據(jù)集2與主流眼電分離方法的效果比較圖

表2 數(shù)據(jù)集2的眼電分離方法結(jié)果統(tǒng)計

5.3.3 公開數(shù)據(jù)集3上與ICA眼電分離方法效果比較

由于公開數(shù)據(jù)集3沒有雙極記錄的單獨垂直眼電信號，不能使用平均偽跡回歸分析方法，因此本文僅與 ICA眼電分離方法進行效果比較如圖 5所示?？梢钥闯鯥CA方法除了眼電去除不徹底外，還存在引入了新的偽跡成分問題。更難以接受的是ICA眼電去除方法造成沒有眼電污染時段較嚴(yán)重的形變，使某些電極波形細節(jié)面目全非。另外，ICA眼電去除方法是半人工參與的方法，需要具有豐富經(jīng)驗的人員參與，通常不同人或不同試次分解的眼電成分?jǐn)?shù)量和順序并不一致，眼電成分選擇不充分則眼電去除不干凈，選擇過多則去除了正常腦電成分，這都增加了人工去除的難度。而且，當(dāng)多個眼電成分在時序上不完全一致時，使用ICA眼電去除方法會造成嚴(yán)重的損失，導(dǎo)致方法崩潰。

表3比較了本文方法和ICA方法去除眼電后的信噪比。可以看出本文方法在3個電極上都比ICA方法有更高的信噪比，眼電去除后的腦電保留了更多的細節(jié)成分。

表3 數(shù)據(jù)集3上的眼電分離方法結(jié)果統(tǒng)計

圖5 公開數(shù)據(jù)集3上與ICA方法眼電分離效果比較圖

6 結(jié)論

針對單通道腦電信號眼電偽跡去除的難點問題，本文提出了一種基于小波變換的自動分離方法。實驗結(jié)果表明基于腦電長時差分振幅包絡(luò)的雙門限眼電端點檢測方法，能夠高效地檢測出幾乎所有的眼電偽跡和大強度的其他偽跡，且算法時延較小?；谛〔ㄗ儞Q的眼電分離方法所估計的眼電偽跡與原始信號之間高度相似。同主流的偽跡平均回歸分析方法和ICA方法相比，本文方法分離的腦電信號具有更高的信噪比，最終得到的純凈腦電不僅保留了正常腦電的高頻部分，還保留了一部分低頻成分。綜上所述，本文方法具有較強的實時性，較好的精確性和良好的眼電分離效果，適于多種腦機接口應(yīng)用中單通道腦電信號的眼電分離。

[1] Doud A J, Lucas J P, Pisansky M T, et al.. Continuous three-dimensional control of a virtual helicopter using a motor imagery based brain-computer interface[J]. PloS One,2011, 6(10): 1-10.

[2] Pei X, Barbour D L, Leuthardt E C, et al.. Decoding vowels and consonants in spoken and imagined words using electrocorticographic signals in humans[J]. Journal of Neural Engineering, 2011, 8(4): 1-11.

[3] Semlitsch H V, Anderer P, Schuster P, et al.. A solution for reliable and valid reduction of ocular artifacts, applied to the P300 ERP[J]. Psychophysiology, 1986, 23(6): 695-703.

[4] Jung T P, Makeig S, Westerfield M, et al.. Removal of eye activity artifacts from visual event-related potentials in normal and clinical subjects[J]. Clinical Neurophysiology,2000, 111(10): 1745-1758.

[5] Gwin J T, Gramann K, Makeig S, et al.. Removal of movement artifact from high-density EEG recorded during walking and running[J]. Journal of Neurophysiology, 2010,103(6): 3526-3534.

[6] 張冰塵, 戴博偉. 一種基于隨機濾波的神經(jīng)動作電位信號壓縮感知采樣方法[J]. 電子與信息學(xué)報, 2013, 35(9): 2283-2286.Zhang Bing-chen and Dai Bo-wei. Compressed sampling for neural action potentials based on random convolution[J].Journal of Electronics & Information Technology, 2013, 35(9):2283-2286.

[7] Garry H, McGinley B, Jones E, et al.. An evaluation of the effects of wavelet coefficient quantisation in transform based EEG compression[J]. Computers in Biology and Medicine,2013, 43(6): 661-669.

[8] 羅志增, 沈寒霄. 基于Hermite插值的小波模極大值重構(gòu)濾波的肌電信號消噪方法[J]. 電子與信息學(xué)報, 2009, 31(4):857-860.Luo Zhi-zeng and Shen Han-xiao. Hermite interpolationbased wavelet transform modulus maxima reconstruction algorithm’s application to EMG de-noising[J]. Journal of Electronics & Information Technology, 2009, 31(4): 857-860.

[9] An Peng. Research on the EEG signal denoising method based on improved wavelet transform[J]. International Journal of Digital Content Technology and Its Applications,2013, 7(4): 154-163.

[10] Chen G. Automatic EEG seizure detection using dual-tree complex wavelet-Fourier features[J]. Expert Systems with Applications, 2014, 41(5): 2391-2394.

[11] Gandhi T, Panigrahi B K, and Anand S. A comparative study of wavelet families for EEG signal classification[J].Neurocomputing, 2011, 74(17): 3051-3057.

[12] Zoughi T, Boostani R and Deypir M. A wavelet-based estimating depth of anesthesia[J]. Engineering Applications of Artificial Intelligence, 2012, 25(8): 1710-1722.

[13] McFarland D J, McCane L M, David S V, et al.. Spatial filter selection for EEG-based communication[J]. Electroencephalography and Clinical Neurophysiology, 1997, 103(3): 386-394.

[14] Donoho D L and Johnstone J M. Ideal spatial adaptation by wavelet shrinkage[J]. Biometrika, 1994, 81(3): 425-455.

[15] Luck S J. An Introduction to the Event-Related Potential Technique[M]. Second edition, Cambridge, MIT Press, 2014:158-162.

[16] Rapela J, Gramann K, Westerfield M, et al.. Brain oscillations in switching vs. focusing audio-visual attention[C].Proceedings of the 34th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, San Diego, USA, 2012: 352-355.