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

吳明權 李海峰 馬 琳

1 引言

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

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

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

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

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

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

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

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

3.1 長時差分信號的計算

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

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

3.2 長時差分信號振幅包絡的提取

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

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

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

圖2 振幅包絡提取流程示意圖

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

σ是正常腦電振幅的標準差， λh, λl為經驗性常數。由于正常腦電長時差分信號近似服從N ( 0,σ)正態分布，因此，一般可取λh≥ 3以越過高幅值的正常腦電包絡，而取λl≥ 2確定眼電起止點。本文選取 λh=5, λl=3來確定眼電區間。

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

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

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

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

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

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

4.1 腦電信號的Mallat小波分解

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

4.2 基于Birgé.-.Massart策略的自適應小波系數閾值確定方法

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

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

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

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

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

4.3 眼電估計的小波重構和與EEG分離

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

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

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

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

5 實驗結果及分析評價

5.1 評測本文方法所使用的數據

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

5.2 量化評價指標

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

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

5.3 實驗結果及分析

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

表1 數據集1的眼電信號分離方法結果統計

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

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

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

表2 數據集2的眼電分離方法結果統計

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

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

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

表3 數據集3上的眼電分離方法結果統計

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

6 結論

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

[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]. 電子與信息學報, 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插值的小波模極大值重構濾波的肌電信號消噪方法[J]. 電子與信息學報, 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.