荆钰霏，李川涛，王 伟，于旭东*

（1.海军军医大学海军特色医学中心特种作战医学研究室，上海 200433；2.海军军医大学海军特色医学中心航空医学研究室，上海 200433）

0 引言

脑电信号由于其非侵入性、灵敏性和准确性应用于脑科学研究具有一定优势[1]。其中前额无发区脑电由于方便采集，在未来可穿戴脑电采集如疲劳监测等领域中有很大的应用前景[2-3]。然而前额无发区距离人眼较近，在采集脑电过程中不可避免地引入较强的眼电伪迹，包含眼动和眨眼产生的电信号，大大地影响了脑电信号的分析处理。一种好的眼电伪迹去除方法有助于提高信号的信噪比和后续疲劳特征识别的准确性[4]。

人工剔除法[5]是在临床脑电图检查中较为常用的方法，医师们往往选择舍弃含有大量眼电信号的数据，保留较为干净的脑电信号进行分析。这种方法需要人为剔除眼电信号，工作量大，并且损失了大量有用信号。盲源分离法[6]由于能将混合信号有效分离成包含噪声和不包含噪声的独立成分，成为近年来热门的伪迹去除方法。但该方法存在分离信号与源信号关系不确定的问题，往往需要手动挑选剔除眼电成分[7]或是引入其他的算法来进行眼电分量的识别。基于频域滤波的伪迹去除方法因其简单有效一直被广泛使用[8]，比如利用无限冲激响应滤波器（infinite impulse responsefilter，IIR）和有限冲激响应滤波器（finite impulse responsefilter，FIR）去除伪迹，然而使用IIR 滤除眼电，眼电尖峰会造成较长时间的振荡而引入新的噪声，使用FIR 滤除眼电，需要较长的阶数以获得较高的噪声抑制能力。基于时频域的小波变换[9]因其良好的时频特性成为最常被使用的伪迹去除方法之一，需要注意的是，眼电伪迹和脑电信号有频带重叠，采用频域滤波的方式在滤除眼电信号的同时也损失了很大一部分脑电信号，因此研究者们在传统小波变换的基础上又联合使用了其他判别方法滤除眼电，比如李洋[9]采用小波变换结合卡尔曼滤波器的方法去除眼电伪迹，周卫东等[10]采用小波软门限法去除脑电中的噪声。

与小波分解相比，经验模态分解（empirical mode decomposition，EMD）吸取了小波变换多分辨力的优势，却没有基函数的选择问题。EMD 在1998 年由Huang 等[11]提出，因其良好的自适应性和短数据序列处理能力被广泛应用于非线性、非平稳数据的分析。传统EMD 算法在信号不属于白噪声时，在分解中会因为时间尺度丢失而造成模态混叠[11]。Torres 等[12]提出的基于自适应噪声完备集成经验模态分解（complete ensemble empirical mode decomposition with adaptive noise，CEEMDAN）方法吸收了集成经验模态分解（ensemble empirical mode decomposition，EEMD）的优点，能有效解决模态混叠问题，并通过逐层添加自适应噪声的分解方法减少了整体的噪声引入，被认为是一种性能良好的信号分解方法。恒虚警率（constant false alarm rate，CFAR）算法是一种雷达动目标检测方法，能在不同的环境噪声背景下保持灵敏的动目标检测能力。

基于此，本研究提出了一种在CEEMDAN 的基础上使用CFAR 算法找到眼电尖峰并去除眼电伪迹的预处理方法，即CEEMDAN-CFAR。首先采用CEEMDAN 的方法分解原始信号，分离出低频噪声信号和含有明显眼电脉冲信号的本征模态函数（intrinsic mode function，IMF）分量。然后再采用CFAR检测算法识别和去除眼电伪迹，重构得到干净的脑电信号。

1 方法

1.1 CEEMDAN

EMD 是将信号分解成一组IMF 的一种信号分解方法。IMF 为满足以下2 个条件的函数：

（1）整个时间历程内，穿越零点的次数与极值点数相等或至多相差1。

（2）信号上任意一点，由局部极大值定义的上包络线和由局部极小值定义的下包络线的均值为0，即信号关于时间轴局部对称。

通过将信号不断分解成IMF 分量的方式，EMD将原始信号的各尺度分量不断从高频到低频进行提取。大量实验证明[11]，EMD 方法类似于一个自适应二进滤波器，能将白噪声分解为一系列中心频率呈二进制关系的IMF[13]。但是当数据不是纯的白噪声时，分解中的时间尺度丢失会导致模态混叠。EEMD 就是为了解决这种模态混叠而产生的[14]。EEMD 在原始信号中加入白噪声以补充原始信号分解中时间尺度不足的问题，并利用白噪声统计特性，通过多次加入噪声求得特征模态函数均值来抵消噪声。EEMD固然在一定程度上解决了模态混叠的问题，然而在重构信号时也引入了新的噪声。同时，不同的噪声添加方式会带来不同的EEMD 分解层数。CEEMDAN通过EEMD 求得高频分量，将剩余的所有分量看成新的信号加入新的自适应噪声再求得高频分量，不断分解到最后一层[12]。这种信号分解方式在降低计算量的同时提供了更精确的重构信号，本文使用CEEMDAN 的方式分解脑电信号，可以有效获得眼动分量和低频噪声分量。假设有原始脑电信号X（t），X（t）的CEEMDAN 具体实现如下：

（1）输入原始脑电信号X（t），令j=1（j 为循环计数变量）。

（2）判断X 的极值个数是否大于2，若大于2 则继续执行（3），否则转到（7）。

（3）令i=1（i 为循环计数变量）。

（4）判断i 是否小于N+1（N 为自定义的重复添加噪声次数），若小于N+1 则执行（5），否则转到（6）。

（5）在原始脑电信号X 中加入白噪声，进行EMD得到第一层IMF1，令i=i+1，转到（4）。

（6）退出循环，计算N 个IMF1 的均值得到IMF（j），得到新的待处理脑电信号X=X-IMF（j），令j=j+1，转到（2）。

（7）退出循环得到全部IMF 分量。

算法流程图如图1 所示。

图1 CEEMDAN 算法流程图

1.2 基于CFAR 算法确定眼动位置

原始脑电信号经由CEEMDAN 方法分解后得到多层IMF 分量，眼电信号在IMF 中表现为振幅突然增大的、类似脉冲的信号。CFAR 算法过去一般被用于雷达动目标检测[15]，本文参考Lazaro 等[16]提到的关于使用方差作为识别单元检测突变信号的方法，将原始信号序列X 以一定窗宽和步长求得各单元方差生成新的序列，如图2 所示。然后通过检测异常方差值来确定眼电信号的位置。将生成的方差序列看成新的信号序列使用CFAR 算法检测异常值。在该算法中保护窗和训练窗长度取值皆为1，CFAR 检测原理如图3 所示。

图2 CFAR 算法识别单元序列生成方法

图3 CFAR 检测原理图

根据经验将CFAR 的阈值T 设为10，识别结果H 为

当H 大于0 时为H1，表示有眼电；当H 小于0 时为H0，表示没有眼电。

2 实验过程及结果

2.1 实验过程

2.1.1 实验条件

本研究以疲劳实验过程中采集到的睁眼脑电信号作为原始信号进行预处理实验。疲劳实验采集被试前额FP1 处脑电信号，以人体左右乳突分别作为参考和地端。脑电采集使用课题组自己研制的双通道前额柔性脑电采集设备，系统采样频率250 Hz，模数转换数据精度为24 bit，短路噪声小于3 μV，10 μV～100 mV 信号测量误差不大于10%，输入阻抗680 MΩ。使用该设备采集一名34 岁健康男性无任务状态下的睁眼脑电数据3 min。

2.1.2 CEEMDAN-CFAR 算法信号处理过程

为具体介绍CEEMDAN-CFAR 方法的信号处理过程，在健康成年人的3 min 睁眼脑电数据中挑选一段10 s 的脑电信号并以此为例进行预处理。10 s含有眼电伪迹的脑电信号如图4 所示，由图可知该信号包含4 个明显的眨眼眼电和低频眼动眼电，是具有代表性的复杂脑电信号。首先，采用CEEMDAN 算法将原始脑电信号分解成多层IMF，CEEMDAN 算法中添加白噪声权重为0.2，重复添加噪声次数为20，将EMD 内部最大包络次数设定为100。CEEMDAN算法分解示意图如图5 所示。

图4 含有眼电伪迹的脑电信号

图5 CEEMDAN 算法分解示意图

由1.1 章节中CEEMDAN 算法的原理介绍可知，该算法具有类似小波的二进滤波器特性，第n 层的中心频率可由公式（2）计算得出：

式中，f（n）为第n 层的中心频率，n 为IMF 所在层数；fs为系统采样频率。为了直观地看到各层IMF 分量的频率分布，本研究采用功率谱估计算法Lomb[17]对各层IMF 分量进行频域转换，结果如图6 所示。

由图6 可以明显看到，各层中心频率逐层降低。正常人的眨眼时间为0.2～0.4 s[18]，眼电信号在频带0～10 Hz 之间，根据公式（2）计算可得IMF4、IMF5、IMF6、IMF7 的中心频率分别为15.62、7.81、3.91、1.95 Hz，本文依据多段信号CEEMDAN 分解后观察到的眼电分布特征以及中心频率同眼电频率接近的IMF5、IMF6、IMF7 作为需要去除眼电的通道。

图6 各IMF 分量Lomb 转换图

依据本文算法流程，使用CFAR 检测法找到眼电尖峰所在的位置。将原始信号序列X 以fs/2 为窗宽、以fs/4 为步长求得各单元方差后生成的新序列Y（n），即图7 中红色曲线（采样频率fs为250 Hz）。CFAR的阈值取值为10，使用公式（1）进行方差检测，方差检测结果异常标记为1（即眼电所在位置），其余标记为0，如图7 所示。图7（a）、（b）、（c）分别为IMF5、IMF6、IMF7 分量的CFAR 检测结果。图7 表明，CFAR算法能准确地识别出眼电尖峰所在的位置。

图7 CFAR 检测结果图

确定眼电尖峰位置后，首先将该层分量中以尖峰位置为中心的小段信号置零，置零信号的长度计算方法为

式中，L 为清零信号段的长度。假设落在第n 层IMF 上的眼电信号频率正好为中心频率，则眼电信号长度为中心频率信号周期长度，即fs/f（n），又由于在EMD 的过程中使用三次样条插值的方法来构造信号的上、下包络线，眼电脉冲信号在分解的过程中会引起周围信号变形，所以计算置零信号长度为2 倍的中心频率信号周期长度。

最后，舍弃低频分量，即把IMF8 及其以后的分量置零，将其余分量重构成不含伪迹的纯净脑电信号。整体的信号处理流程如图8 所示。

图8 基于CEEMDAN-CFAR 算法的眼电伪迹去除流程图

2.2 实验结果

以随机挑选的健康成年人10 s 放松状态下的睁眼脑电数据为例进行预处理，含有眼电伪迹的脑电信号（如图4 所示）经CEEMDAN-CFAR 处理后的信号如图9 所示。由图9 可以直观地看到眼电伪迹全部被去除。将原始信号和处理后的信号局部放大后进行比较，如图10 所示。由图10 可以看出，预处理前后信号一致性较强，证实了使用该预处理方法能在去除伪迹的同时很好地保留局部细节和有用信息。

图9 CEEMDAN-CFAR 算法处理后的信号

图10 原始信号和处理后的信号局部对比图

为进一步验证预处理算法的有效性，使用该方法处理完整的3 min 睁眼脑电数据。预处理前后对比如图11 所示，具体细节如图12 所示。由图11、12可以看出，本文提出的预处理方法针对真实的脑电数据具有良好的眼电噪声滤除效果。

图11 3 min 脑电数据预处理前后效果图

图12 3 min 脑电数据预处理细节图

2.3 算法性能

使用Klados 等[19]公开发表的一个半模拟脑电数据集对算法的有效性进行验证。在纯净的脑电信号中人为地加入眼电伪迹模拟真实的脑电信号，如图13所示。采集27 名健康受试者闭眼期间的脑电数据，信号的采样频率为200 Hz，采用0.5～40 Hz 的带通滤波和50 Hz 的陷波滤波，共获得54 个数据集。每个数据集都包含连续30 s 的脑电信号，对获得的数据集进行了仔细检查，确保没有显著污染，可以作为纯净的脑电信号。在相同情况下采集27 名健康受试者的眼电图（electro-oculogram，EOG）数据，EOG 信号在0.5～5Hz下进行带通滤波。根据公式（4）、（5）模拟生成含有眼电伪迹的半模拟脑电信号：

图13 半模拟数据集中的纯净脑电信号和污染脑电信号

式中，Scon表示污染的脑电信号；Spure表示纯净的脑电信号；SEOG表示EOG信号；SVEOG表示垂直EOG信号；SHEOG表示水平EOG信号；α 和β 分别表示垂直EOG 和水平EOG 的污染系数。

从54 个数据集中随机选择12 个数据集中的污染脑电信号（如图14 所示），使用CEEMDAN-CFAR算法对FP1 通道的污染的脑电信号进行处理，处理后得到较为纯净的脑电信号，如图15 所示。

图14 原始的半模拟脑电信号

图15 CEEMDAN-CFAR 算法处理后的半模拟脑电信号

将组成半模拟脑电信号的纯净脑电信号进行CEEMDAN，去除其中低频成分后再与处理后的半模拟脑电信号进行相关性分析，得到相关系数见表1。由表1 中数据可知，经过预处理的信号与纯净脑电信号的相关性得到了很大提升。与样本熵-完备经验模态分解（sample entropy-complete ensemble empirical mode decomposition with adaptive noise，SE-CEEMDAN）[20]算法的相关系数对比：CEEMDAN-CFAR 的相关系数（均值）为0.797，方差为0.062，SE-CEEMDAN 的相关系数（均值）为0.841，方差为0.005。处理200 Hz采样频率下5 s 数据所用时间的对比：SE-CEEMDAN运行时间为68 s，CEEMDAN-CFAR 运行时间为12 s。

表1 纯净脑电信号与处理前、后脑电信号的相关系数

综上，CEEMDAN-CFAR 算法在单通道脑电信号预处理性能上有较高的准确性，虽然与同类算法SE-CEEMDAN[20]相比保留原始信号的能力较弱，但速度较快，是一种方便快捷的预处理方法。

3 结语

脑电信号作为一种微弱的生物电信号，其波幅一般都在μV 的数量级上，特别容易受到外界干扰。在针对脑电的信号处理中，各种伪迹的存在严重影响了信号特征提取的准确性。因此，有效的伪迹去除方法是整个脑电信号处理流程中至关重要的一环。其中眼电伪迹由于其幅度大、随机性强、频率与脑电信号存在频带重叠等特点，成为脑电信号伪迹去除中的重点与难点。传统的眼电伪迹去除方法存在丢失有用信号或引入冗余噪声等问题，而近年来较为热门的盲源分离算法大多是针对多通道脑电数据而言的。在单通道脑电伪迹去除方法的研究上，虽然也有研究者[21]提出了将单通道脑电信号通过希尔伯特-黄变换等方法转化为多通道数据后，再使用独立成分分析（independent component analysis，ICA）分离眼电伪迹的方法，但从ICA 分离得到的各独立源信号中判断眼电成分又需要加入新的算法，这大大增加了信号处理的难度。在此基础上，本文提出基于CEEMDAN-CFAR算法眼电伪迹去除方法，通过CEEMDAN 算法对原始脑电信号进行分解，得到低频噪声成分和含有明显眼电脉冲信号的IMF 分量，然后采用CFAR 算法准确定位眼电位置并予以剔除。该方法具有计算简单、眼电伪迹去除效果好的优势，与同类算法SECEEMDAN 相比，CEEMDAN-CFAR 算法保留原始信号的能力稍弱，但速度快捷。不足之处在于该方法对于低频信号的保留度较低，不适用于研究深度睡眠等需要关注低频脑电信号的领域。在下一步的研究中，将进一步对CEEMDAN 分层后相对低频的部分进行处理，通过各低频IMF 分量的频域特征判别是否需要保留该分量，以此获得更为纯净的脑电信号。