单通道脑电信号中眼电干扰的自动分离方法

吴明权 李海峰 马 琳

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.