郭遠哲，王正勇，胡瀅濱，梁 鑫，何小海*

（1.四川大學電子信息學院，成都 610065；2.成都錦城學院電子信息學院，成都 611731）

0 引言

腦電（electroencephalogram，EEG）作為一種非線性、非平穩(wěn)的生理電信號包含豐富的人體信息，但由于幅值微弱、信噪比（signal to noise ratio，SNR）低等特性易受到多種偽跡的干擾，若不有效去除其中的干擾偽跡，將對后續(xù)的EEG分析工作帶來很大影響。EEG信號采集過程中出現(xiàn)的干擾偽跡主要分為兩大類：第一類是來自外部的非生理性偽跡，如采集過程中的50 Hz工頻干擾、電極接觸不良干擾、電磁干擾等[1]，這些物理性干擾會使得EEG信號形成多個尖銳的波形，在測試中可以遠離供電線，減少電極線擺動來規(guī)避以上干擾；第二類是由受試者測試過程中生理行為產生的生理性偽跡，如眼電（electrooculogram，EOG）、肌電、心電等偽跡。其中，EOG偽跡振幅往往是EEG信號振幅的上百倍，頻率與EEG信號部分頻段吻合，且由于受試者眨眼的隨機性，EOG偽跡會隨機出現(xiàn)在EEG信號中從而造成嚴重干擾[2]，因此對于EOG偽跡的排查與剔除一直是EEG信號研究中的重要內容。

當前對EEG信號中EOG偽跡的去除方式有以下幾種[3]：一是基于時域回歸的方式，在采集EEG信號的同時采集EOG信號并進行偽跡去除。但在實際采集過程中可能面臨EOG信號與EEG信號雙向污染的問題，不能采集到完全純凈的EOG信號。因此偽跡去除過程必然會造成EOG的過估計，從而去除掉有用的EEG信號。二是采用經驗模態(tài)分解（empirical mode decomposition，EMD）的方式。EMD是一種非線性信號分解技術，作為數(shù)據(jù)驅動型及自適應型算法，無需任何先驗知識便可對信號進行分解[4-5]，但該算法對噪聲敏感度極高，會造成模態(tài)混疊的問題。Wu等[6]通過添加白噪聲來解決模態(tài)混疊的問題，由此提出了集成EMD。但該類方法計算量較大，算法實時性不足。三是通過盲源分離（blind source separation，BSS）算法去除EOG偽跡。BSS算法能在源信號和傳輸系統(tǒng)特性均未知的情況下，對混合信號進行成分分離，剔除偽跡成分后再重建信號。目前有許多基于BSS算法設計的標準化EEG信號預處理工具箱[7-9]。常見的BSS算法有主成分分析（principal component analysis，PCA）、獨立成分分析（independent component analysis，ICA）等[10-11]，但該類算法需要手動挑選EOG成分[12]，且未對識別出的EOG分量進行分析就直接置零或者剔除，從而導致混疊于EOG分量中的少量有效EEG信號丟失[13]。針對上述問題，本研究提出一種新的EOG偽跡自動去除方法。該方法在BSS算法的基礎上引入自適應樣本熵，通過EEG信號與EOG信號時間序列復雜度上的差異，對不同被試進行獨立成分的篩選，在自動化篩選的基礎上有效避免因不同個體差異性導致的篩選偏差。同時結合經驗模態(tài)變換防止混疊于EOG分量中的有效EEG信號丟失，利用自相關系數(shù)來區(qū)分EEG信號與EOG信號在時域上的振幅差異及產生的隨機性差異，最后進行EOG偽跡的篩選與剔除。

1 基于快速獨立成分分析（fast independent component analysis，F(xiàn)astICA）識別包含EOG偽跡的獨立成分

ICA可以看成是PCA和因子分析（factor analysis，F(xiàn)A）的一種延伸擴展，該算法的作用是從多維統(tǒng)計信號中分解出內在因素和成分，僅需觀測信號便可分解出各路原信號。假定觀測信號X=[x1，x2，…，xN]T，未知源信號S=[s1，s2，…，sM]T，ICA的數(shù)學模型如下：

式中，A為N×M階滿秩常數(shù)的混合矩陣；ε（t）為高斯白噪聲。為簡化模型，分析過程可忽略由高斯白噪聲帶來的干擾。在進行ICA前，必須滿足源信號是平穩(wěn)隨機信號，相互間獨立且服從非高斯分布。

FastICA是一種在獨立成分求解過程中引入固定點迭代遞推的分析方式[14]，與傳統(tǒng)的ICA相比，它具有收斂速度快、效率高等優(yōu)勢，因而被廣泛運用于日常使用場景中。在進行ICA前，需對信號進行中心化處理和白化處理以減小估計參數(shù)量，提升分解速率，降低計算復雜度。中心化處理即對觀測信號去除均值，其表達式如下：

式中，E{·}為均值運算。白化處理目的是去除相關性，白化后矩陣Z與白化矩陣V的關系為Z=V。為保障白化后觀測信號的協(xié)方差矩陣為單位矩陣，需滿足下式：

式中，D為協(xié)方差矩陣特征值組成的對角矩陣；P為協(xié)方差矩陣特征向量組成的正交矩陣?？汕蟮冒谆仃嚢谆幚砗蟮玫降木仃嘮便具有單位方差、均值為零、各分量獨立的特性，然后根據(jù)牛頓-拉夫遜方法[15]采用如下公式迭代后得到分離矩陣W：

式中，g（X）為任意的非線性函數(shù)，本文為雙曲正切函數(shù)。對W（m+1）矩陣進行歸一化并不斷進行重復迭代和矩陣單位化至相鄰2次W收斂在預設閾值內，或達到最大迭代次數(shù)，即完成獨立成分的分解。在對原始EEG信號完成ICA后，需要從源信號集合中識別出包含EOG偽跡的獨立成分（electrooculogram independent component，EOG-IC）。傳統(tǒng)回歸算法和頻域濾波算法在EOG偽跡去除過程中存在識別難度大的問題，而基于單一閾值的篩選由于個體差異的原因容易對于不同被試出現(xiàn)較大的偏差，造成有用的EEG信號丟失。

本研究采取自適應樣本熵閾值篩選的方式來判斷該分量中是否含有EOG信號。樣本熵是一種通過度量信號中產生新模式的概率來衡量時間序列復雜度的方法。通過與原始信號的部分頻段功率失真率的對比，實現(xiàn)自適應樣本熵閾值的選取，具體計算方法如下：

（1）以m為窗體大小、N為長度的時間序列X可表示為

（2）將X（i）與X（j）間對應元素最大的差值定義為兩者的距離，記為d[X（i），X（j）]。給定閾值r，統(tǒng)計d[X（i），X（j）]＜r的數(shù)量與總矢量數(shù)N-m的比值，記為。對所有求得的計算平均值，得到該窗體下的Bm（r）。

（3）將窗體大小m加1，重復上述步驟。本文中窗體大小m為2，r為0.25std（X），其中std（X）表示標準差。當前m和r下的樣本熵計算公式如下：

（4）將獨立成分的樣本熵由小到大進行排序，依次遍歷每一個獨立分量的樣本熵，將當前樣本熵作為樣本熵閾值。由于EOG分量處于低頻部分，因此基于當前樣本熵閾值下去除EOG偽跡后的EEG信號的α波段與β波段應與原始EEG信號中的α波段與β波段功率相匹配。

（5）選取8~30 Hz部分作為對比頻段，將重構信號設為Q、原始信號設為C，計算信號重構前后EEG信號該頻段的功率損失ΔP并進行對比。ΔP的計算公式如下：

式中，PQ為重構信號當前頻段的功率；PC為原始信號當前頻段的功率。

（6）找到丟失信號最小情況及求得功率損失ΔP最小時的樣本熵作為閾值，將所有樣本熵小于該閾值的獨立成分分量認定為帶有EOG偽跡的分量。

2 基于經驗小波變換（empirical wavelet transform，EWT）算法對EOG偽跡進行剔除

因基于FastICA找到的EOG-IC中仍然包含部分EEG信號，為精確消除EOG偽跡，需要進一步對EOG-IC進行分解操作，找到并剔除其中的EOG偽跡成分。本研究運用EWT算法，結合自相關系數(shù)對EOG-IC進行自適應分解與剔除。EWT是由Gilles[16]在小波變換思路的基礎上結合EMD提出的一種自適應信號處理方法，其不僅具有EWT的科學性還具有EMD的自適應性。由于EWT基于小波框架建立，因此比傳統(tǒng)EMD計算量更小、分解更精確，解決了實時性不足、本征模函數(shù)模態(tài)混疊等問題。

EWT算法包含信號頻譜的分割、構造經驗小波和對所有獨立段應用經驗小波處理3個重要方面。EWT的分割方式是將目標信號的傅里葉頻譜規(guī)范化到[0，π]后，分割為N段連續(xù)且獨立的部分。分割后頻帶集合為{Λ1，Λ2，…，ΛN}，所有的獨立段用Λn=[ωn-1，ωn]表示，其中n=1，2，…，N。EWT頻譜分割示意圖如圖1所示。

圖1 EWT頻譜分割示意圖

在完成分割后，Gilles[16]基于Meyer小波的構建思路，將經驗小波函數(shù)Ψn（ω）和尺度函數(shù)Φn（ω）定義如下：

式中，v為輔助函數(shù)，定義為

將 過 渡 區(qū) 的 寬 度 設 置 為τn=γωn，0＜γ＜minn在確定好經驗小波函數(shù)和經驗尺度函數(shù)后，經驗小波變換的近似系數(shù)與細節(jié)系數(shù)可分別用經驗尺度函數(shù)Φn（ω）和經驗小波函數(shù)Ψn（ω）同信號f（t）的內積表示[17]。信號f（t）通過EWT可得到經驗模態(tài)分量：

由EOG-IC通過EWT算法分解而來的經驗模態(tài)分量依據(jù)EEG與EOG的特性，運用信號的自相關系數(shù)閾值法來進行兩者的辨識，自相關系數(shù)（AutoCC）作為一種信號與時延信號相關程度的評價指標，其計算公式如下：

式中，x（t）為分解的經驗模態(tài)分量信號；x（t-τ）為時延τ的經驗模態(tài)分量信號；Cov（x）為協(xié)方差；Var（x）為方差。EEG信號振幅較低、復雜度較高、隨機性更強，因此所求得的自相關系數(shù)值會在較低水平。而EEG信號振幅大、頻率低、瞬時性等特點使自相關系數(shù)值會更高。依據(jù)自相關系數(shù)閾值對經驗模態(tài)分量進行劃分，將EOG模態(tài)分量去除后進行EWT逆變換，合成當前成分分量，并與其他獨立成分分量進行重構，合成不含EOG的純凈EEG信號。EOG偽跡識別與去除整體流程圖如圖2所示。

圖2 EOG偽跡識別與去除整體流程圖

3 試驗結果與分析

3.1 數(shù)據(jù)采集

為保障數(shù)據(jù)來源的可靠性，本研究選取9名健康被試者（均為男性，右利手，年齡為22~24歲）在睜眼靜息狀態(tài)下的EEG測試數(shù)據(jù)（以下EEG波形均來自被試者1）。采集設備為64通道的腦電儀設備，在確認儀器工作正常與被試者精神狀態(tài)正常后，進行5 min的EEG信號采集，采樣頻率設置為512 Hz。采集過程中要求被試者盡量減少身體活動，避免EEG信號被除EOG偽跡外的其他生理性偽跡影響。

在去除偽跡前需對信號進行濾波。本研究在采集階段通過采集軟件去除50 Hz工頻干擾，并用巴特沃斯濾波器進行低通濾波，頻率設置范圍為0.1~60.0 Hz。根據(jù)EOG偽跡形成的條件與其衰減的規(guī)律，偽跡干擾較為嚴重的區(qū)域是位于被試者大腦頭皮前部區(qū)域。在FP1、FPz、FP2、AF7、AF3、AF4、AF8這7處通道，EOG偽跡對EEG信號干擾較大，其他通道由于與眼睛相隔較遠，EOG在傳播過程中逐步衰減，對EEG信號的干擾較小。被試者1的7通道EEG信號和獨立成分分量波形圖如圖3所示。其中，圖3（a）經過FastICA可得到獨立成分分量波形，如圖3（b）所示。

圖3 被試者1的7通道EEG信號和獨立成分分量波形圖

計算各獨立成分的樣本熵，由IC1至IC7獨立成分分量對應的樣本熵依次為1.352、1.228、1.355、0.351、1.279、1.046、0.963。

3.2 數(shù)據(jù)處理

根據(jù)自適應樣本熵閾值的選擇原理，各閾值下8~30 Hz頻段的功率失真率如圖4所示。結果表明，以IC7的樣本熵作為閾值時，該頻段功率失真率最小及此閾值下的原始EEG有用信號得到保留情況較好，因此可認定當前EEG信號的獨立成分樣本熵的閾值應選為0.963。在分離出被試者1的7個獨立分量中，由于自適應樣本獨立成分4下的樣本熵（SEIC4）小于此時的樣本熵閾值，因此可以認為IC4為EOG-IC。根據(jù)EWT算法并結合頻率特性，對挑選出的EOG-IC進行頻譜分割，分解為4個經驗模態(tài)分量，如圖5~6所示。

圖5 EOG-IC頻譜分割圖

圖6 各經驗模態(tài)分量圖

本研究在延時τ為1的情況下計算自相關系數(shù)，當自相關系數(shù)閾值設置為0.99時，可將EEG成分與EOG偽跡區(qū)分開[18]。求得EOG-IC的4個經驗模態(tài)分量自相關系數(shù)依次為0.999、0.994、0.975、0.698?？梢钥吹皆赟1分量與S2分量處，自相關系數(shù)值超過閾值，因此判定為EOG分量。剔除EOG偽跡過后，進行逆EWT運算得到去除EOG偽跡后的IC4。將新的IC4與其他獨立成分分量進行重構，得到偽跡去除后的EEG信號。運用本研究方法，被試者1的EOG偽跡去除前后的EEG信號對比圖如圖7所示。

3.3 數(shù)據(jù)對比

本試驗在自采集數(shù)據(jù)集與Klados等[19]公開的半模擬EOG/EEG數(shù)據(jù)集（選取27位被試者）的基礎上，選擇對應通道段的EEG信號進行算法驗證，利用均方根誤差（root mean square error，RMSE）及SNR作為評價指標對所提出的EOG偽跡消除方式的性能進行比較。將選取的7個通道原始信號記為XGround、去除偽跡后的信號記為XClean，RMSE計算公式如下：

RMSE值越小表明XClean與XGround越接近，消除性能越好。SNR計算公式如下：

式中，RMS（X）表示均方根算子。SNR值越大表明XClean與XGround越接近，消除性能越好。

被試者1不同通道、自采集人員和數(shù)據(jù)集平均評估數(shù)據(jù)詳見表1~3。由表中可以看出，本文算法的RMSE及SNR均優(yōu)于近2 a運用較廣泛、典型性較強的ICA算法與wICA（wavelet-enhanced ICA）算法。在不同被試者以及不同通道下本文算法去除EOG偽跡效果更好，偽跡去除后的EEG信號與純凈EEG更為貼合，有效EEG信號丟失更少，且對不同被試者的穩(wěn)定性也優(yōu)于另外2種算法。由于加入了自適應樣本熵，EOG偽跡的篩選也具有更強的魯棒性，從而證明了本研究提出的多通道EOG偽跡去除算法的可靠性和有效性。

表1 3種方法下被試者1不同通道的評估數(shù)據(jù)

表2 3種方法下不同被試者的評估數(shù)據(jù)

表3 3種方法下不同數(shù)據(jù)集的平均評估數(shù)據(jù)

4 結語

本研究提出了一種FastICA與EWT相結合的方法，對EEG信號中EOG偽跡進行自動識別與去除，改善了傳統(tǒng)算法在EOG偽跡去除中的有用信息丟失、偽跡辨別難度大、個體特異性強、去除不徹底等情況。通過對比發(fā)現(xiàn)，在偽跡識別方面通過引入自適應樣本熵閾值，有效提高了自動識別EOG偽跡的效率，節(jié)省了人工判定偽跡成分的時間，提升了EOG信號自動識別的準確性和魯棒性。同時在偽跡去除方面引入EWT而非直接丟棄EOG-IC分量，進一步保留了原始EEG信號中的有用信息。從試驗結果與對比結果可以看出，對于不同被試者及不同通道下的偽跡去除效果，本研究提出方法均更優(yōu)，從而證明了其有效性。本研究的不足之處在于該方法在自適應樣本熵求解過程中時間復雜度較高，下一步將在EOG偽跡識別流程中采用速度更快的自適應特征進行偽跡識別，以更高效地獲得純凈的EEG信號。