多通道脑电信号中眼电伪迹的自动识别与去除方法研究

郭远哲，王正勇，胡滢滨，梁 鑫，何小海*

（1.四川大学电子信息学院，成都 610065；2.成都锦城学院电子信息学院，成都 611731）

0 引言

脑电（electroencephalogram，EEG）作为一种非线性、非平稳的生理电信号包含丰富的人体信息，但由于幅值微弱、信噪比（signal to noise ratio，SNR）低等特性易受到多种伪迹的干扰，若不有效去除其中的干扰伪迹，将对后续的EEG分析工作带来很大影响。EEG信号采集过程中出现的干扰伪迹主要分为两大类：第一类是来自外部的非生理性伪迹，如采集过程中的50 Hz工频干扰、电极接触不良干扰、电磁干扰等[1]，这些物理性干扰会使得EEG信号形成多个尖锐的波形，在测试中可以远离供电线，减少电极线摆动来规避以上干扰；第二类是由受试者测试过程中生理行为产生的生理性伪迹，如眼电（electrooculogram，EOG）、肌电、心电等伪迹。其中，EOG伪迹振幅往往是EEG信号振幅的上百倍，频率与EEG信号部分频段吻合，且由于受试者眨眼的随机性，EOG伪迹会随机出现在EEG信号中从而造成严重干扰[2]，因此对于EOG伪迹的排查与剔除一直是EEG信号研究中的重要内容。

当前对EEG信号中EOG伪迹的去除方式有以下几种[3]：一是基于时域回归的方式，在采集EEG信号的同时采集EOG信号并进行伪迹去除。但在实际采集过程中可能面临EOG信号与EEG信号双向污染的问题，不能采集到完全纯净的EOG信号。因此伪迹去除过程必然会造成EOG的过估计，从而去除掉有用的EEG信号。二是采用经验模态分解（empirical mode decomposition，EMD）的方式。EMD是一种非线性信号分解技术，作为数据驱动型及自适应型算法，无需任何先验知识便可对信号进行分解[4-5]，但该算法对噪声敏感度极高，会造成模态混叠的问题。Wu等[6]通过添加白噪声来解决模态混叠的问题，由此提出了集成EMD。但该类方法计算量较大，算法实时性不足。三是通过盲源分离（blind source separation，BSS）算法去除EOG伪迹。BSS算法能在源信号和传输系统特性均未知的情况下，对混合信号进行成分分离，剔除伪迹成分后再重建信号。目前有许多基于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信号时间序列复杂度上的差异，对不同被试进行独立成分的筛选，在自动化筛选的基础上有效避免因不同个体差异性导致的筛选偏差。同时结合经验模态变换防止混叠于EOG分量中的有效EEG信号丢失，利用自相关系数来区分EEG信号与EOG信号在时域上的振幅差异及产生的随机性差异，最后进行EOG伪迹的筛选与剔除。

1 基于快速独立成分分析（fast independent component analysis，FastICA）识别包含EOG伪迹的独立成分

ICA可以看成是PCA和因子分析（factor analysis，FA）的一种延伸扩展，该算法的作用是从多维统计信号中分解出内在因素和成分，仅需观测信号便可分解出各路原信号。假定观测信号X=[x1，x2，…，xN]T，未知源信号S=[s1，s2，…，sM]T，ICA的数学模型如下：

式中，A为N×M阶满秩常数的混合矩阵；ε（t）为高斯白噪声。为简化模型，分析过程可忽略由高斯白噪声带来的干扰。在进行ICA前，必须满足源信号是平稳随机信号，相互间独立且服从非高斯分布。

FastICA是一种在独立成分求解过程中引入固定点迭代递推的分析方式[14]，与传统的ICA相比，它具有收敛速度快、效率高等优势，因而被广泛运用于日常使用场景中。在进行ICA前，需对信号进行中心化处理和白化处理以减小估计参数量，提升分解速率，降低计算复杂度。中心化处理即对观测信号去除均值，其表达式如下：

式中，E{·}为均值运算。白化处理目的是去除相关性，白化后矩阵Z与白化矩阵V的关系为Z=V。为保障白化后观测信号的协方差矩阵为单位矩阵，需满足下式：

式中，D为协方差矩阵特征值组成的对角矩阵；P为协方差矩阵特征向量组成的正交矩阵。可求得白化矩阵白化处理后得到的矩阵Z便具有单位方差、均值为零、各分量独立的特性，然后根据牛顿-拉夫逊方法[15]采用如下公式迭代后得到分离矩阵W：

式中，g（X）为任意的非线性函数，本文为双曲正切函数。对W（m+1）矩阵进行归一化并不断进行重复迭代和矩阵单位化至相邻2次W收敛在预设阈值内，或达到最大迭代次数，即完成独立成分的分解。在对原始EEG信号完成ICA后，需要从源信号集合中识别出包含EOG伪迹的独立成分（electrooculogram independent component，EOG-IC）。传统回归算法和频域滤波算法在EOG伪迹去除过程中存在识别难度大的问题，而基于单一阈值的筛选由于个体差异的原因容易对于不同被试出现较大的偏差，造成有用的EEG信号丢失。

本研究采取自适应样本熵阈值筛选的方式来判断该分量中是否含有EOG信号。样本熵是一种通过度量信号中产生新模式的概率来衡量时间序列复杂度的方法。通过与原始信号的部分频段功率失真率的对比，实现自适应样本熵阈值的选取，具体计算方法如下：

（1）以m为窗体大小、N为长度的时间序列X可表示为

（2）将X（i）与X（j）间对应元素最大的差值定义为两者的距离，记为d[X（i），X（j）]。给定阈值r，统计d[X（i），X（j）]＜r的数量与总矢量数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算法，结合自相关系数对EOG-IC进行自适应分解与剔除。EWT是由Gilles[16]在小波变换思路的基础上结合EMD提出的一种自适应信号处理方法，其不仅具有EWT的科学性还具有EMD的自适应性。由于EWT基于小波框架建立，因此比传统EMD计算量更小、分解更精确，解决了实时性不足、本征模函数模态混叠等问题。

EWT算法包含信号频谱的分割、构造经验小波和对所有独立段应用经验小波处理3个重要方面。EWT的分割方式是将目标信号的傅里叶频谱规范化到[0，π]后，分割为N段连续且独立的部分。分割后频带集合为{Λ1，Λ2，…，ΛN}，所有的独立段用Λn=[ωn-1，ωn]表示，其中n=1，2，…，N。EWT频谱分割示意图如图1所示。

图1 EWT频谱分割示意图

在完成分割后，Gilles[16]基于Meyer小波的构建思路，将经验小波函数Ψn（ω）和尺度函数Φn（ω）定义如下：

式中，v为辅助函数，定义为

将 过 渡 区 的 宽 度 设 置 为τn=γωn，0＜γ＜minn在确定好经验小波函数和经验尺度函数后，经验小波变换的近似系数与细节系数可分别用经验尺度函数Φn（ω）和经验小波函数Ψn（ω）同信号f（t）的内积表示[17]。信号f（t）通过EWT可得到经验模态分量：

由EOG-IC通过EWT算法分解而来的经验模态分量依据EEG与EOG的特性，运用信号的自相关系数阈值法来进行两者的辨识，自相关系数（AutoCC）作为一种信号与时延信号相关程度的评价指标，其计算公式如下：

式中，x（t）为分解的经验模态分量信号；x（t-τ）为时延τ的经验模态分量信号；Cov（x）为协方差；Var（x）为方差。EEG信号振幅较低、复杂度较高、随机性更强，因此所求得的自相关系数值会在较低水平。而EEG信号振幅大、频率低、瞬时性等特点使自相关系数值会更高。依据自相关系数阈值对经验模态分量进行划分，将EOG模态分量去除后进行EWT逆变换，合成当前成分分量，并与其他独立成分分量进行重构，合成不含EOG的纯净EEG信号。EOG伪迹识别与去除整体流程图如图2所示。

图2 EOG伪迹识别与去除整体流程图

3 试验结果与分析

3.1 数据采集

为保障数据来源的可靠性，本研究选取9名健康被试者（均为男性，右利手，年龄为22~24岁）在睁眼静息状态下的EEG测试数据（以下EEG波形均来自被试者1）。采集设备为64通道的脑电仪设备，在确认仪器工作正常与被试者精神状态正常后，进行5 min的EEG信号采集，采样频率设置为512 Hz。采集过程中要求被试者尽量减少身体活动，避免EEG信号被除EOG伪迹外的其他生理性伪迹影响。

在去除伪迹前需对信号进行滤波。本研究在采集阶段通过采集软件去除50 Hz工频干扰，并用巴特沃斯滤波器进行低通滤波，频率设置范围为0.1~60.0 Hz。根据EOG伪迹形成的条件与其衰减的规律，伪迹干扰较为严重的区域是位于被试者大脑头皮前部区域。在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 数据处理

根据自适应样本熵阈值的选择原理，各阈值下8~30 Hz频段的功率失真率如图4所示。结果表明，以IC7的样本熵作为阈值时，该频段功率失真率最小及此阈值下的原始EEG有用信号得到保留情况较好，因此可认定当前EEG信号的独立成分样本熵的阈值应选为0.963。在分离出被试者1的7个独立分量中，由于自适应样本独立成分4下的样本熵（SEIC4）小于此时的样本熵阈值，因此可以认为IC4为EOG-IC。根据EWT算法并结合频率特性，对挑选出的EOG-IC进行频谱分割，分解为4个经验模态分量，如图5~6所示。

图5 EOG-IC频谱分割图

图6 各经验模态分量图

本研究在延时τ为1的情况下计算自相关系数，当自相关系数阈值设置为0.99时，可将EEG成分与EOG伪迹区分开[18]。求得EOG-IC的4个经验模态分量自相关系数依次为0.999、0.994、0.975、0.698。可以看到在S1分量与S2分量处，自相关系数值超过阈值，因此判定为EOG分量。剔除EOG伪迹过后，进行逆EWT运算得到去除EOG伪迹后的IC4。将新的IC4与其他独立成分分量进行重构，得到伪迹去除后的EEG信号。运用本研究方法，被试者1的EOG伪迹去除前后的EEG信号对比图如图7所示。

3.3 数据对比

本试验在自采集数据集与Klados等[19]公开的半模拟EOG/EEG数据集（选取27位被试者）的基础上，选择对应通道段的EEG信号进行算法验证，利用均方根误差（root mean square error，RMSE）及SNR作为评价指标对所提出的EOG伪迹消除方式的性能进行比较。将选取的7个通道原始信号记为XGround、去除伪迹后的信号记为XClean，RMSE计算公式如下：

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

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

被试者1不同通道、自采集人员和数据集平均评估数据详见表1~3。由表中可以看出，本文算法的RMSE及SNR均优于近2 a运用较广泛、典型性较强的ICA算法与wICA（wavelet-enhanced ICA）算法。在不同被试者以及不同通道下本文算法去除EOG伪迹效果更好，伪迹去除后的EEG信号与纯净EEG更为贴合，有效EEG信号丢失更少，且对不同被试者的稳定性也优于另外2种算法。由于加入了自适应样本熵，EOG伪迹的筛选也具有更强的鲁棒性，从而证明了本研究提出的多通道EOG伪迹去除算法的可靠性和有效性。

表1 3种方法下被试者1不同通道的评估数据

表2 3种方法下不同被试者的评估数据

表3 3种方法下不同数据集的平均评估数据

4 结语

本研究提出了一种FastICA与EWT相结合的方法，对EEG信号中EOG伪迹进行自动识别与去除，改善了传统算法在EOG伪迹去除中的有用信息丢失、伪迹辨别难度大、个体特异性强、去除不彻底等情况。通过对比发现，在伪迹识别方面通过引入自适应样本熵阈值，有效提高了自动识别EOG伪迹的效率，节省了人工判定伪迹成分的时间，提升了EOG信号自动识别的准确性和鲁棒性。同时在伪迹去除方面引入EWT而非直接丢弃EOG-IC分量，进一步保留了原始EEG信号中的有用信息。从试验结果与对比结果可以看出，对于不同被试者及不同通道下的伪迹去除效果，本研究提出方法均更优，从而证明了其有效性。本研究的不足之处在于该方法在自适应样本熵求解过程中时间复杂度较高，下一步将在EOG伪迹识别流程中采用速度更快的自适应特征进行伪迹识别，以更高效地获得纯净的EEG信号。