自动化学报  2017, Vol. 43 Issue (10): 1726-1735   PDF    
单通道脑电信号眼电伪迹去除算法研究
刘志勇1,2, 孙金玮1, 卜宪庚3     
1. 哈尔滨工业大学电气工程及自动化学院 哈尔滨 150001;
2. 61135 部队 北京 102211;
3. 哈尔滨医科大学基础医学院 哈尔滨 150081
摘要: 由眨眼和眼动产生的眼电伪迹(Electrooculography,EOG)信号是脑电信号(Electroencephalography,EEG)中的主要噪声信号之一.目前,多通道脑电信号中眼电伪迹的去除算法已经较为成熟.而在单通道脑电信号的眼电伪迹去除中,由于采集通道数量较少且缺乏参考眼电信号,目前尚无十分有效的去除方法.本文提出一种基于小波变换(Wavelet transform,WT)、集合经验模态分解(Ensemble empirical mode decomposition,EEMD)和独立成分分析(Independent component analysis,ICA)的WT-EEMD-ICA单通道脑电信号眼电伪迹去除算法.实验表明:WT-EEMD-ICA算法有效地解决了单通道WT-ICA算法中的超完备问题,能够有效去除单通道脑电信号中的眼电伪迹,并且分离出的眼电伪迹成分与参考通道采集的眼电信号相关性较强.
关键词: 脑电信号     眼电伪迹     小波变换     集合经验模态分解     独立成分分析    
EOG Artifact Removing Method for Single-channel EEG Signal
LIU Zhi-Yong1,2, SUN Jin-Wei1, BU Xian-Geng3     
1. School of Electrical Engineering and Automation, Harbin Institute of Technology, Harbin 150001;
2. 61135 PLA Troop, Beijing 102211;
3. School of Basic Medical, Harbin Medical University, Harbin 150081
Manuscript received : March 9, 2016, accepted: July 18, 2016.
Foundation Item: Supported by National Natural Science Foundation of China (61301012), Sci-tech Innovation Foundation of Harbin (2015RA XXJ038), and Fundamental Research Funds for the Central Universities (2013004, 2013005)
Author brief: LIU Zhi-Yong  Ph. D. candidate at the School of Electrical Engineering and Automation, Harbin Institute of Technology. He received his bachelor degree from Harbin Institute of Technology in 2008. His main research interest is biomedical signal processing;
BU Xian-Geng  Professor at the School of Basic Medical, Harbin Medical University. His research interest covers internet of things and its application
Corresponding author. SUN Jin-Wei  Professor at the School of Electrical Engineering and Automation, Harbin Institute of Technology. His research interest covers biomedical sensors and active noise control theory. Corresponding author of this paper, E-mail:jwsun@hit.edu.cn
Recommended by Associate Editor TIAN Jie
Abstract: Electrooculography (EOG) artifacts generated by eye movements and blinks are the major artifacts in electroencephalography (EEG). There are many common effective methods for removing the multi-channel EEG artifacts. However, due to the limitation of input channel number and the lack of reference EOG signal, there is no very effective artifact removing method for single-channel EEG signal. In the present study, a novel EOG artifact removing algorithm WT-EEMD-ICA for single-channel EEG signal is proposed based on wavelet transform (WT), ensemble empirical mode decomposition (EEMD) and independent component analysis (ICA) technologies. The result shows that the WT-EEMDICA method, which successfully solves the overcomplete problem of WT-ICA in single channel artifact removal, can separate the EOG and EEG successfully only from one single-channel EEG, and that the useful information involved in original EEG signal can be greatly saved.
Key words: Electroencephalography(EEG)     electrooculography(EOG)     wavelet transform(WT)     ensemble empirical mode decomposition(EEMD)     independent component analysis(ICA)    

脑电信号(Electroencephalography, EEG)中包含了大量的生理和病理信息.从大脑头皮采集到的脑电信号包含了大量的生理伪迹.常见的伪迹种类包括心电伪迹(Electrocardiography, ECG)[1]、肌电伪迹(Electromyography, EMG)[2]和眼电伪迹(Electrooculography, EOG)[3]等.由于这些伪迹信号的存在, 使得采集到的脑电信号受到了严重的噪声污染.如果伪迹成分得不到有效去除, 会导致在后续的脑电信号特征参数提取与分类过程中出现严重的错误.特别是在脑机接口(Brain computer interface, BCI)[4]和麻醉深度监测等领域, 过多的伪迹干扰会造成分类结果的不准确, 并带来严重的后果[5].眼电伪迹的频带范围分布较广[6], 且眼电伪迹与脑电信号频带严重重叠[7], 这使得传统线性滤波的方法难以去除脑电信号中的眼电伪迹[8].

近年来, 许多眼电伪迹去除算法相继问世.根据采集的通道数量, 可以将伪迹去除算法分为多通道眼电伪迹去除算法和单通道眼电伪迹去除算法.多通道伪迹去除算法主要包括伪迹减法[9]、主成分分析法(Principal components analysis, PCA)[10]和独立成分分析法(Independent component analysis, ICA)[11]等. Gomez等基于盲源分离算法提出了多通道眼电伪迹的自动去除算法, 在不采用眼电参考通道的情况下成功地对脑电信号中的眼电伪迹进行了去除[12], 优点是能够处理高阶统计量[13], 缺点是需要伪迹的先验知识[14].

在便携式脑机接口(BCI)、便携式麻醉深度监测系统、便携式睡眠监测和便携式情绪状态识别等应用领域通常针对单通道的脑电信号进行分析.因此, 研究可靠的单通道脑电信号伪迹去除算法对于便携式脑电设备的发展具有重要意义.在单通道脑电信号的伪迹去除中, 由于通道数量较少, 包含的有效信息有限, 且缺少参考通道, 使得二者的分离较为困难.最初, 人们采用阈值的方法限制眼电伪迹成分的幅度, 但会损失有用的脑电数据成分[15], 特别是当眼电伪迹幅度与脑电信号幅度接近时[16].线性滤波方法能够去除部分眼电伪迹成分[17], 而对眼电伪迹与脑电信号频带重叠的部分无能为力[18].吴明权等提出了一种基于长时差分振幅包络与小波变换(Wavelet transform, WT)的眼电干扰自动分离方法, 通过长时差分结合双门限法确定眼电伪迹起始位置, 随后采用sym 5小波对脑电信号进行分解, 同时引入Birge-Massart策略确定小波重构系数阈值对眼电信号进行估计, 进而实现单通道脑电信号的眼电伪迹分离[19].该方法的缺点是当眼电伪迹与脑电信号幅度较为接近时, 对眼电信号进行准确估计十分困难.近年来, 不少学者尝试将单通道脑电信号转换为多通道信号, 然后采用PCA或ICA等较为成熟的多通道伪迹去除算法进行处理. James等采用动态嵌入的方法将单通道信号转换为多通道信号, 结合ICA算法实现了单通道脑电信号的眼电伪迹去除算法[20-21]. Mammone等采用ICA对小波分解的结果进行独立成分分析, 提出了一种基于WT-ICA方法的眼电伪迹去除算法[22], 在多通道脑电信号中能够取得较好的效果, 但在单通道伪迹去除中存在超完备ICA问题.为解决该问题, 本文在其基础上引入集合经验模态分解(Ensemble empirical mode decomposition, EEMD)算法, 解决了超完备ICA问题, 实现了单通道脑电信号的眼电伪迹去除.

1 WT-EEMD-ICA算法

图 1为WT-EEMD-ICA算法的框图.其中, EEG$_{\rm ac}$为采集到的包含眼电伪迹的脑电信号, EEG$_{\rm nc}$和EOG$_{\rm nc}$分别为从EEG$_{\rm ac}$中分离出的纯净的脑电信号和眼电信号.该算法首先采用小波变换(WT)对EEG$_{\rm ac}$进行分解.经过WT后, 单通道脑电信号被分解为多个小波系数.由于脑电信号与眼电伪迹存在严重的频带重叠, 因此, 脑电信号和眼电伪迹会同时出现在多个小波系数中, 而且分解后的各小波系数之间是相互独立的, 各小波系数包含的脑电分量和眼电分量也是相互独立的.各小波系数中的脑电分量和眼电分量可以认为是原始脑电信号分量和眼电分量的子集.随后, 分别采用集合经验模态分解(EEMD)依次将各个小波系数进一步分解为若干个内蕴模态函数(Intrinsic mode functions, IMFs), 每个小波系数的所有IMFs函数全部输入独立成分分析(ICA)系统进行脑电分量和眼电分量的分离, 进而采用近似熵阈值的方法对二者进行识别.最后, 通过ICA逆变换将脑电分量和眼电分量转化为EEG$_{\rm nc}$和EOG$_{\rm nc}$进行输出.

图 1 WT-EEMD-ICA算法框图 Figure 1 The block of WT-EEMD-ICA algorithm
1.1 小波变换

小波变换是继傅里叶分析后的又一重要的信号处理技术, 具有里程碑式的重要意义.该技术按照不同的时间尺度对时域信号进行多分辨分析, 可以同时在时域和频域获得较好的局部特性.信号的小波分解需要设置相应的小波基函数$\Psi$和分解层数M.通过对小波基函数进行伸缩和平移来生成标准正交基函数.函数$f(x)$的小波变换定义为

$ \begin{align} &{W_f}(a, b) = \langle f(x), {\Psi _{a, b}}\rangle =\notag\\ &\qquad \frac{1}{{\sqrt a }}\int_R {f(x)\Psi \left(\frac{{x - b}}{a}\right){\rm d}x} \end{align} $ (1)

${W_f}(a, b)$的运算过程相当于采用一系列带通滤波器对信号进行多个频率范围内的滤波, 从而得到不同频带内的信息, 即实现按照不同的频带范围对信号进行分解.

1.2 经验模态分解

经验模态分解是由Huang等提出的一种信号分解技术, 能够将信号分解为若干内蕴模态函数IMFs[23].所有的内蕴模态函数都满足以下两个条件: 1) 极值点和零点数相等或最多相差1; 2) 由极大值点定义的包络线和由极小值点定义的包络线均值为0.信号$x(t)$的经验模态分解如式(2) 所示.相对于小波变换, 经验模态分解不需要设置基函数.但有时候会造成模态混叠, 因此, 为了避免模态混叠, Wu和Huang又提出了集合经验模态分解(EEMD).

$ \begin{align} x(t) = \sum\limits_{n = 1}^N {{c_n}(t)} + {r_n}(t) \end{align} $ (2)

其中, ${c_n}(t)$是第n阶IMFs函数, ${r_n}(t)$为余项, N为IMFs函数个数.

1.3 独立成分分析

假设$X=[x_1, x_2, \cdots, x_N]$N维观测信号, S = $[s_1, s_2, \cdots, s_M]$M维源信号.独立成分分析算法ICA的模型可以表示为

$ \begin{align} X = A \times S \end{align} $ (3)

其中, A$M \times N$维矩阵.

ICA的主要原理就是在AS均为未知的情况下, 通过一定的判定准则对源信号S进行估计.通常, ICA要求观测信号的数量不少于源信号的数量.当观测信号的数量少于源信号数量时会出现超完备ICA现象.本文在ICA运算时采用FastICA算法. FastICA也称固定点算法[24], 是一个在批处理模式下的快速优化迭代算法, 可基于四阶统计量、最大似然估计或最大负熵等参数进行计算.相对于传统ICA算法, FastICA收敛速度更快且不需要设置步长参数, 这使得该算法相比于基于梯度的算法更容易.但FastICA结果中各独立成分的顺序不确定, 不能够确定独立成分对应的信号源.并且结果幅度并不是各独立成分的真实幅度.

1.4 基于近似熵的脑电分量识别

ICA结果中的各分量可以看作各信号源, 重构脑电信号的过程是将其他非脑电信号分量置零, 然后通过式(3) 得到重构的EEG$_{\rm nc}$信号.在信息论中, 熵指数能够表征系统的复杂程度.脑电信号是大脑细胞内外生物电活动的外在表现, 能够表征大脑状态.而眼电信号源自眨眼活动或眼球运动.相对于眼电信号, 脑电信号则更为复杂.本文采用近似熵[25]作为脑电分量与眼电分量的判据.

近似熵的具体算法为:

步骤 1. 对系统进行等间隔采样, 得到时间序列$u(1) $, $u(2) , \cdots, u(N)$;

步骤 2. 计算时间序列$u(i)$的标准差$SD$;

步骤 3. 设定比较序列长度参数m以及阈值参数r的值, 取$m=2$, $r=0.1\times SD$;

步骤 4. 构造新的m维矢量$X(1) , X(2) , \cdots$, $X(N- m +1) $, 其中$X(i)=[u(i), u(i+1) , \cdots, u(i$ $+$ $m-1) ]$;

步骤 5. 构造距离矩阵$d[X(i), X(j)]$, $1 \leq i$, j $\leq$ $N-m+1$,

$ \begin{align} d[X(i), X(j)] = \mathop {\max }\limits_{k = 0, \cdots, m - 1} [\left| {u(i + k)-u(j + k)} \right|] \end{align} $ (4)

步骤 6. 对于新的时间序列$X(i)=[u(i), u(i+1) , \cdots, u(i+m-1) ]$, $1 \leq i \leq N-m+1$, 记为

$ \begin{align} C_i^m(r) = \frac{1}{{N - m + 1}}\left\{ d[x(i), x(j)] \le r ~{\rm numbers} \right\} \end{align} $ (5)

步骤 7. 定义$\Phi^{m}(r)$

$ \begin{align} {\Phi ^m}(r) = \frac{1}{{N - m + 1}}\sum\limits_{i = 1}^{N - m + 1} {\ln C_i^m(r)} \end{align} $ (6)

步骤 8. 将比较序列长度参数m加1, 即$m=$ $m+ 1$, 重复步骤3 ~ 7, 得到$\Phi^{m+1}(r)$;

步骤 9. 计算近似熵$ApEn$ \begin{align} ApEn = {\Phi ^m}(r) -{\Phi ^{m + 1}}(r) \end{align}

2 实验 2.1 实验设计

图 2为本文实验框图.实验对象为一名健康成年男子, 年龄为24岁.实验设备采用基于ADS1299的多通道生物信号采集开发板. EEG$_{\rm ac}$和EOG$_{\rm ac}$分别为采集到的脑电信号和水平眼电信号. EEG$_{\rm ac}$为标准10 ~ 20导联系统中的$Fp1$$A1$两个位置的差分输出. EOG$_{\rm ac}$为左右眼水平电极的差分输出.信号采样率为250 Hz, 采集时间约为30分钟.实验开始前, 首先对面部皮肤进行处理, 去除油脂和角质降低表面肌电伪迹成分.实验过程中, 实验对象要保持平静状态并在自然状态下眨眼和水平眼动若干次.因此, 采集到EEG$_{\rm ac}$包含的其他伪迹成分相对较少, 可以近似认为只包含脑电信号和眼电伪迹成分.而采集到的EOG$_{\rm ac}$可近似认为是纯净的眼电信号.

图 2 实验设计框图 Figure 2 The block of experiment design
2.2 效果评价

EEG$_{\rm ac}$中包含明显的眼电伪迹.同时, 水平眼电通道采集到的EOG$_{\rm ac}$也包含少量的脑电信号成分.这就使得对新算法的伪迹去除效果评价较为困难.通常, 对伪迹去除算法的评价包含如下三个方面: 1) 脑电信号中的眼电伪迹成分明显去除; 2) 在眼电伪迹不明显的部分, 脑电信号中的有价值信息得到最大程度的保留; 3) EEG$_{\rm ac}$中分离出的眼电伪迹EOG$_{\rm nc}$与对比通道采集到的眼电信号EOG$_{\rm ac}$高度相关.眼电伪迹去除效果可以从主观进行判断.脑电信号中有价值信息的保留程度以及分离出的EOG$_{\rm nc}$和EOG$_{\rm ac}$的相关性可以通过相关系数进行评价.计算公式如下:

$ R{\rm{ = }}\frac{{\sum\limits_{i = 1}^n {({x_i} - \bar x)({y_i} - \bar y)} }}{{\sqrt {\sum\limits_{i = 1}^n {{{({x_i} - \bar x)}^2}} } \sqrt {\sum\limits_{i = 1}^n {{{({y_i} - \bar y)}^2}} } }} $ (8)

为了进一步评价WT-EEMD-ICA算法, 对比了小波算法以及文献[21]中提出的WT-ICA算法.同时, 为了突出小波变换在整个算法中的作用, 对比分析了EEMD-ICA算法.此外, 为了更合理验证本文算法的有效性, 对比分析了两种其他眼电伪迹去除算法.所有的对比实验包括:

1) WT算法.采用小波变换对原始脑电信号进行分析, 去除含有眼电伪迹的小波系数, 用所有不含眼电伪迹的小波系数重构EEG$_{\rm nc}$信号.

2) WT-ICA算法.采用ICA算法对所有小波系数进行独立成分的计算, 通过脑电分量的ICA逆变换重构EEG$_{\rm nc}$信号.

3) EEMD-ICA算法.首先采用EEMD算法将原始脑电信号分解为若干IMFs, 然后采用ICA算法对所有IMFs进行计算, 最后, 通过脑电分量的ICA逆变换重构EEG$_{\rm nc}$信号.

4) DE-ICA算法.文献[20]提出的动态嵌入(Dynamical embedding, DE)结合ICA算法.动态嵌入的计算如式(9) 所示, 目的是为了将单通道的脑电信号转化为多个子信号.

$ \begin{array}{l} X(t) = \\ \quad \{ x(t),x(t{\rm{ + }}\tau ),x(t{\rm{ + }}2\tau ), \cdots ,x(t{\rm{ + }}(m - 1)\tau )\} \end{array} $ (9)

5) 希尔伯特黄变换(Hilbert-Huang transform, HHT)的眼电伪迹去除算法.详细算法见文献[26].

3 结果与分析

图 3分别为从采集对象头部获得的长度为8秒的EEG$_{\rm ac}$信号和水平EOG$_{\rm ac}$信号, 分别记为$P1$$P2$.从图 3(a)中可以看出明显的眼电伪迹成分, 从图 3(b)中可以看出存在两处较大的眼电伪迹信号.

图 3 EEG$_{\rm ac}$和EOG$_{\rm ac}$波形 Figure 3 Wave of EEG$_{\rm ac}$和EOG$_{\rm ac}$

图 4$P1$时间段的脑电信号和眼电信号波形图.可以发现眼电伪迹成分主要集中于61 ~ 62秒内.相对来讲, 59 ~ 61秒($S1$)和62 ~ 63秒($S2$)两个时间段的脑电信号中不包含明显的眼电伪迹成分.因此, 采用$S1$段和$S2$段脑电信号经过WT-EEMD-ICA算法处理前后的相关系数来评价该算法对脑电信号中所包含有价值信息的保留程度.

图 4 $P1$段EEG$_{\rm ac}$和EOG$_{\rm ac}$波形 Figure 4 Wave of EEG$_{\rm ac}$和EOG$_{\rm ac}$ in $P1$
3.1 预处理

EEG$_{\rm ac}$和EOG$_{\rm ac}$中包含一定的直流成分, 高频噪声以及工频干扰.在进行下一步处理前, 需要进行预处理.脑电信号和眼电信号的频率成分主要集中在64 Hz以下, 因此首先采用通带频率为0.5 ~ 64 Hz的带通滤波器对EEG$_{\rm ac}$和EOG$_{\rm ac}$进行滤波处理. 图 5$P1$段EEG信号和EOG信号经带通滤波处理前后的时域和频域波形对比.图 6$P1$段EEG信号和EOG信号经过50 Hz带阻滤波处理前后的时域和频域波形对比.

图 5 带通滤波前后EEG与EOG信号对比 Figure 5 EEG and EOG processed by band pass filter
图 6 带阻滤波前后EEG与EOG信号对比 Figure 6 EEG and EOG processed by band stop filter
3.2 WT-EEMD-ICA结果

图 7为经过7层小波分解后的各小波系数图.对经过预处理的EEG$_{\rm ac}$进行小波变换, 首先要选择小波基函数.为了使得分解后的小波系数能够更明显地反映眼电伪迹成分, 最佳的小波基函数要与眼电伪迹波形较为接近, 并且具有正交性、双正交性以及紧支撑等特点.经过对比, 小波基函数选择sym7小波, 分解层数为7层.图 7中, 多个小波系数都出现了较为明显的EOG伪迹成分, 特别是在$d6$, $d5$$d4$中.

图 7 EEG信号小波分解图 Figure 7 Wavelet coefficients of EEG

随后, 采用EEMD算法分别对每个小波系数进行分解.以$d4$为例, 图 8$d4$小波系数经过EEMD分解后的IMFs结果.最后将所有IMFs函数输入到ICA系统中计算脑电和眼电的分量. 图 9$d4$小波系数的ICA结果.

图 8 $d4$小波系数EEMD结果图 Figure 8 The EEMD results of $d4$
图 9 $d4$小波系数ICA结果图 Figure 9 The ICA results of $d4$

图 9可知, 眼电伪迹成分主要出现在$IC1$$IC5$中. $IC1$ ~ $IC5$分量的近似熵分别为: 0.60, 1.65, 1.13, 1.04, 0.18.可以发现脑电分量的近似熵值相对眼电分量明显较高. $P1$段脑电信号所有小波系数的脑电分量和眼电分量的近似熵均值与标准差如表 1所示.通过设置适当的近似熵判定阈值, 可以准确区分脑电分量和眼电分量.因此, 采用$IC2$ ~ $IC4$等成分的ICA逆变换重构EEG$_{\rm nc}$, 采用$IC1$$IC5$重构EOG$_{\rm nc}$.需要指出的是, $IC1$ ~ $IC5$所示的幅度并不是各独立分量的真实幅度, 经过ICA逆变换后才是真实的信号幅度.

表 1 脑电分量与眼电分量近似熵对比 Table 1 Comparison of $ApEn$ between EEG and EOG

图 10为从$d4$小波系数中分离出的脑电信号和眼电信号与$d4$小波系数和对比通道EOG$_{\rm ac}$信号的波形对比图. 图 10(a)10 (b)中的虚线曲线分别为从$d4$小波系数中重构的EEG信号和EOG眼动信号, 图 10(a)中的实线曲线为$d4$小波系数波形, 图 10(b)中的实线曲线为对比通道采集的EOG信号波形.从图 10(a)可知, $d4$小波系数中EOG伪迹成分的去除效果十分明显, 并且在前2秒和最后1秒内的EEG信号处理前后一致性较好.图 10(b)从小波系数$d4$中重构的EOG伪迹信号波形与对比通道的EOG波形变化趋势较为接近, 且能够准确反映出在第2 ~ 3秒内的眼动伪迹.从图 10(b)还可以发现, 仅从$d4$小波系数重构出的EOG眼动伪迹波形与对比通道的EOG波形还存在一定差距, 因此需要对其他小波系数进行同样处理, 并将最后结果相融合得到最终的伪迹去除结果.

图 10 $d4$小波系数WT-EEMD-ICA结果 Figure 10 The WT-EEMD-ICA results of $d4$

图 11$d3$小波系数的处理结果.与$d4$小波系数相似, $d3$小波系数的分析结果能较准确地去除EEG信号中的EOG伪迹, 去除EOG伪迹后, EEG中的有效信息保留较好.值得注意的是, 从图 7的小波系数波形图中很难发现$d3$小波系数中包含的EOG伪迹成分.因此, 图 11的结果证明仅依靠直观判断信号中是否包含EOG伪迹成分并不是十分可靠的, 需要对每个小波系数进行单独处理才能分辨是否包含EOG伪迹成分.

图 11 $d3$小波系数WT-EEMD-ICA结果 Figure 11 The WT-EEMD-ICA results of $d3$

图 12为采用WT-EEMD-ICA算法对$P1$段脑电信号进行处理后的结果. $S1$段和$S2$段脑电信号经过WT-EEMD-ICA算法处理前后的相关系数分别为0.82和0.74, 证明该算法能够有效去除EOG伪迹, 并保留非伪迹部分的EEG信号包含的有价值信息.图 12(b)中, 从单通道EEG信号中分离出的EOG伪迹与对比通道EOG波形较为接近, 能够准确识别眼电伪迹的位置, 特别是在第2 ~ 3秒时间段内, 大部分的EOG眼动伪迹峰值都较为吻合.

图 12 $P1$段信号WT-EEMD-ICA结果 Figure 12 The WT-EEMD-ICA results of signals in $P1$
3.3 WT-EEMD-ICA算法稳定性

前文采用WT-EEMD-ICA算法对$P1$段脑电信号和眼电信号进行了分析.为了验证WT-EEMD-ICA算法的稳定性, 对$P2$段脑电信号和眼电信号进行了同样处理, 并对$P1 + P2$段8秒长度的脑电信号和眼电信号进行WT-EEMD-ICA处理, 如图 13图 14所示.

图 13 $P2$段信号WT-EEMD-ICA结果 Figure 13 The WT-EEMD-ICA results of signals in $P2$
图 14 $P1+P2$段信号WT-EEMD-ICA结果 Figure 14 The WT-EEMD-ICA results of signals in $P1+P2$

图 13(a)图 14(a)可知, 重构的$P2$段EEG信号以及$P1 + P2$段EEG信号中EOG伪迹已经被明显去除.在图 13(b)图 14(b)中, 从单通道EEG信号中重构的EOG信号相对于参考通道的EOG信号定位较为准确, 且二者幅度较为接近, 证明了采用WT-EEMD-ICA算法从单通道EEG信号中重构的EEG和EOG信号的有效性, 以及WT-EEMD-ICA算法的有效性.

图 15中, 方法1为分别对$P1$$P2$段信号进行处理, 方法2为直接对$P1 + P2$段信号进行处理.图 15(a)为从两种方法分离出的脑电信号波形对比图, 二者相关系数为0.86.图 15(b)为采用同样方法对眼电信号的处理结果, 二者相关系数为0.84.

图 15 WT-EEMD-ICA算法稳定性实验结果 Figure 15 Stability results of WT-EEMD-ICA
3.4 对比实验结果与分析

图 16 ~ 20分别为采用WT算法、WT-ICA算法、EEMD-ICA算法、DE-ICA算法和HHT算法对$P1$段脑电信号和眼电信号的处理结果.其中, 由于WT-ICA算法存在超完备问题, 导致伪迹去除结果不稳定, 多次运行WT-ICA的结果重复性较差.图 17(a)为其中一次运行的结果.从图 16(a)17(a)可以看出, WT和WT-ICA算法明显地去除了$P1$段脑电信号眼电伪迹, 而EEMD-ICA算法没有实现眼电伪迹的有效去除.

图 16 $P1$段信号WT算法结果 Figure 16 The WT algorithm results of signals in $P1$
图 17 $P1$段信号WT-ICA算法结果 Figure 17 The WT-ICA algorithm results of signals in $P1$
图 18 $P1$段信号EEMD-ICA算法结果 Figure 18 The EEMD-ICA algorithm results of signals in $P1$
图 19 $P1$段信号DE-ICA算法结果 Figure 19 The DE-ICA algorithm results of signals in $P1$
图 20 $P1$段信号HHT算法结果 Figure 20 The HHT algorithm results of signals in $P1$

DE-ICA算法进行单通道信号伪迹去除时, 首先需要确定嵌入维数与延迟常数.采用G-P算法计算时间序列相空间重构的最佳嵌入维数$m=12$和延迟常数$\tau=5$.图 19为采用动态嵌入的方法将原始EEG信号转化为多个子信号并进行ICA处理的结果.从图 19可知, DE-ICA的结果中, 每个独立成分中均含有明显眼电伪迹, 可以证明该方法没能有效区分脑电信号分量和眼电伪迹分量, 无法有效去除单通道脑电信号中的眼电伪迹.

根据文献[26]的建议, HHT伪迹去除算法中的$fl$$fh$分别设置为1 Hz和30 Hz. 图 20为采用HHT算法对$P1$段信号进行伪迹去除的结果.图 20(b)中重构的EOG相对其他几种方法高频分量明显偏高.

表 2$S1$$S2$段脑电信号经过四种对比算法去除眼电伪迹前后的相关系数对比表, 还包括$P1$段脑电信号中分离出的眼电信号EOG$_{\rm nc}$$P1$段对比通道眼电信号EOG$_{\rm ac}$的相关系数.

表 2 重构信号与原始信号相关系数 Table 2 The correlation coefficients between reconstructed and original signal
4 讨论

从头皮采集的脑电信号和眼电信号不可避免地会引入一些高频噪声、直流成分、工频干扰和一些生理伪迹.本文对脑电和眼电信号的采集方式进行合理设计, 从源头降低了信号中的干扰成分.经过带通滤波和带阻滤波后的脑电信号可以近似认为仅包含纯净的脑电成分和少量的眼电成分.由图 5(b)图 5(d)可知, 经过带通滤波后, 脑电信号和眼电信号高于64 Hz成分的功率谱成分均小于0 dB.在图 6(b)6 (d)中, 50 Hz工频干扰成分经过带阻滤波后明显降低.值得注意的是, 脑电信号中50 Hz工频干扰的功率谱比眼电信号中的工频干扰功率谱幅度相对较高.造成这种现象的原因是脑电信号的采集方式是单极采集, 而眼电信号的采集方式是双极采集.

图 7中, 经过小波变换, 眼电伪迹成分出现在部分小波系数中, 特别是$d6$, d5和$d4$小波系数中的眼电伪迹成分较为明显.但并不是所有的小波系数都能够直接观察到其中的眼电伪迹成分.例如, 在对$d3$小波系数进行EEMD和ICA处理后, ICA的结果中含有较为明显的眼电伪迹成分.如表 1所示, 通过近似熵值可以准确区分脑电分量和眼电分量.由图 10(a)可知, 眼电伪迹成分去除效果非常明显并且对于$S1$$S2$段眼电伪迹成分不明显的部分, 脑电信号保持较好.图 10(b)中, 仅从$d4$小波系数中分离出眼电信号大致上与对比通道的EOG$_{\rm ac}$信号趋势保持一致.由于缺乏其他小波系数中的眼电伪迹成分, 结果与对比通道的EOG$_{\rm ac}$信号还存在一定差距.

图 12为采用WT-EEMD-ICA算法对$P1$段脑电信号和眼电信号的处理结果. $S1$$S2$段脑电信号经过WT-EEMD-ICA算法处理前后的相关系数分别为0.82和0.74.证明该算法对EEG信号中的有价值信息保留较好.图 13(b)中, 从$P1$段脑电信号中分离出的EOG$_{\rm nc}$与对比通道EOG$_{\rm ac}$波形较为相似, 能够准确识别眼电伪迹的位置, EOG$_{\rm nc}$与EOG$_{\rm ac}$峰值较为接近, 由于EOG$_{\rm ac}$中含有一些高频的脑电信号, 二者相关系数仅为0.54.值得注意的是, EOG$_{\rm ac}$中的一些其他伪迹成分没有分离出来.例如在1.6秒左右的眼电信号.造成这种现象的原因是该时刻的伪迹成分并不是真正的眼电信号, 而是受到了其他干扰信号的影响, 而这种影响并没有体现在从$Fp1-A1$采集的EEG$_{\rm ac}$信号中.

为了验证WT-EEMD-ICA算法在不同数据长度下的稳定性, 本文分别对$P1$段(长度4秒)、$P2$段(长度4秒)和$P1 + P2$段(长度8秒)的信号进行WT-EEMD-ICA算法处理, 结果分别如图 12~ 14所示.将分别对$P1$段和$P2$段脑电信号去除眼电伪迹后的结果和对$P1 + P2$段脑电信号去除眼电伪迹后的结果进行对比, 结果如图 15所示.二者的EEG$_{\rm nc}$相关系数为0.86, EOG$_{\rm nc}$相关系数为0.84.可以证明在不同序列长度下, WT-EEMD-ICA算法具有较好的稳定性.

为了进一步验证WT-EEMD-ICA算法的伪迹去除效果, 本文对比了WT、WT-ICA、EEMD-ICA、DE-ICA和HHT等五种算法. WT算法将原始信号进行小波分解, 去除含有眼电伪迹的小波系数, 采用所有不含眼电伪迹成分的小波系数重构EEG$_{\rm nc}$信号.由于含有眼电伪迹的小波系数中含有较多的脑电信号, 直接去除含有眼电伪迹的小波系数必定会造成有用脑电信号的损失.由图 16(a)可知, 采用WT算法重构的EEG$_{\rm nc}$信号, 与原始EEG$_{\rm ac}$波形存在较大差距, 意味着WT伪迹去除方法虽然去除了眼电伪迹成分, 但对原始信号中的脑电信号有价值成分也损失较多.

文献[21]提出的WT-ICA算法在应用于单通道脑电信号伪迹去除时, 经过小波变换, 单通道脑电信号转换为若干个小波系数, 似乎满足了传统ICA的输入输出条件, 即观测信号的数目不小于源信号数量.但是, 由于每个小波系数仍然包含多个独立成分.因此, 该条件下的ICA仍然是一个超完备的输入输出问题.这也是导致该算法结果不稳定的主要原因.实验结果表明, 该算法每次运行的结果均不一样.图 17为其中一次的运行结果图.与WT算法类似, 其结果中EEG$_{\rm nc}$与EEG$_{\rm ac}$信号在眼电伪迹不明显的部分却发生了较大的变化, 损失了较多有价值的脑电信息.图 18中, EEMD-ICA算法的眼电伪迹去除效果更差, 重构的EEG$_{\rm nc}$中包含非常明显的眼电伪迹成分, 未达到伪迹去除的效果.经过对比, WT-EEMD-ICA算法克服了WT-ICA和EEMD-ICA算法中伪迹成分去除不彻底, 有价值脑电信号损失严重的问题, 实现了较好的单通道眼电伪迹去除.并且, 在处理不同长度的脑电信号时得到的结果较为稳定.

文献[20]提出的DE-ICA算法在多通道脑电信号的伪迹去除中效果较好, 但在去除单通道脑电信号的伪迹时效果较差.图 19的处理结果显示, DE-ICA的处理结果中, 眼电伪迹出现在所有独立分量中, 说明该方法不能够有效区分脑电信号和眼电伪迹成分. 图 20中基于HHT的眼电伪迹去除算法有效去除了$P1$段脑电信号的眼电伪迹, 但重构的$S1$$S2$段脑电信号存在一定程度的信息损失, 导致该段信号处理前后的相关系数较低, 分别为0.64和0.57.从表 2的对比中可知, 本文提出的WT-EEMD-ICA算法相对于其他四种算法具有明显的优势.

5 结论

由于受采集通道数量少, 缺乏参考眼电通道等因素的制约, 单通道脑电信号的眼电伪迹去除问题一直没有得到较好解决.本文提出的WT-EEMD-ICA单通道眼电伪迹去除算法融合了小波变换、集合经验模态分解和独立成分分析等算法, 解决了ICA算法在去除单通道眼电伪迹过程中的超完备问题, 实现了单通道脑电信号中眼电伪迹去除.相对于WT、WT-ICA和EEMD-ICA三种算法, WT-EEMD-ICA在眼电伪迹去除方面的效果更优, 并且对脑电信号中的有价值信息保留程度较好.经过对比验证, WT-EEMD-ICA算法在处理不同长度脑电信号时, 具有较好的稳定性.

参考文献
1
Marker R J, Maluf K S. Effects of electrocardiography contamination and comparison of ECG removal methods on upper trapezius electromyography recordings. Journal of Electromyography and Kinesiology, 2014, 24(6): 902-909. DOI:10.1016/j.jelekin.2014.08.003
2
Hu Jin, Hou Zeng-Guang, Chen Yi-Xiong, Zhang Feng, Wang Wei-Qun. Lower limb rehabilitation robots and interactive control methods. Acta Automatica Sinica, 2014, 40(11): 2377-2390.
( 胡进, 侯增广, 陈翼雄, 张峰, 王卫群. 下肢康复机器人及其交互控制方法. 自动化学报, 2014, 40(11): 2377-2390.)
3
Zeng H, Song A G, Yan R Q, Qin H Y. EOG artifact correction from EEG recording using stationary subspace analysis and empirical mode decomposition. Sensors, 2013, 13(11): 14839-14859. DOI:10.3390/s131114839
4
Wang Xing-Yu, Jin Jing, Zhang Yu, Wang Bei. Brain control:human-computer integration control based on braincomputer interface. Acta Automatica Sinica, 2013, 39(3): 208-221.
( 王行愚, 金晶, 张宇, 王蓓. 脑控:基于脑-机接口的人机融合控制. 自动化学报, 2013, 39(3): 208-221.)
5
Jie D, Tao W, Zhang A T, Dai H Y. The removal of blink and saccade artifact in EEG recordings by independent component analysis. In:Proceedings of the 3rd International Conference on Biomedical Engineering and Informatics. Yantai, China:IEEE, 2010, 3: 1071-1075.
6
Daly I, Scherer R, Billinger M, Müller-Putz G. FORCe:fully online and automated artifact removal for brain-computer interfacing. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2015, 23(5): 725-736. DOI:10.1109/TNSRE.2014.2346621
7
Mannan M M N, Jeong M Y, Kamran M A. Hybrid ICAregression:automatic identification and removal of ocular artifacts from electroencephalographic signals. Frontiers in Human Neuroscience, 2016, 10: 193-193.
8
Qi W. EOG artifacts removal in EEG measurements for affective interaction with brain computer interface. In:Proceedings of the 8th International Conference on Intelligent Information Hiding and Multimedia Signal Processing. Piraeus-Athens, Greece: IEEE, 2012, 471-475.
9
Chang W D, Cha H S, Kim K, Im C H. Detection of eye blink artifacts from single prefrontal channel electroencephalogram. Computer Methods and Programs in Biomedicine, 2016, 124: 19-30. DOI:10.1016/j.cmpb.2015.10.011
10
Vigon L, Saatchi M R, Mayhew J E W, Fernandes R. Quantitative evaluation of techniques for ocular artefact filtering of EEG waveforms. IEE Proceedings-Science Measurement and Technology, 2000, 147(5): 219-229. DOI:10.1049/ip-smt:20000475
11
Fatourechi M, Bashashatia A, Ward R K, Birch G E. EMG and EOG artifacts in brain computer interface systems:a survey. Clinical Neurophysiology, 2007, 118(3): 480-494. DOI:10.1016/j.clinph.2006.10.019
12
Gomez-Herrero G, Clercq W D, Anwar H, Kara O, Egiazarian K, van Huffel S, van Paesschen W. Automatic removal of ocular artifacts in the EEG without an EOG reference channel. In:Proceedings of the 7th Nordic Signal Processing Symposium. Reykjavik, Iceland: IEEE, 2006, 130-133.
13
Chen Yong-Qiang, Wang Hong-Xia. A method for underdetermined blind source separation based on new mixture model. Acta Automatica Sinica, 2014, 40(7): 1412-1420.
( 陈永强, 王宏霞. 基于新型混合模型的欠定盲分离方法. 自动化学报, 2014, 40(7): 1412-1420.)
14
Jung T P, Makeig S, Westerfield M, Townsend J, Courchesne E, Sejnowski T J. Removal of eye activity artifacts from visual event-related potentials in normal and clinical subjects. Clinical Neurophysiology, 2000, 111(10): 1745-1758. DOI:10.1016/S1388-2457(00)00386-2
15
del R Millan J, Mourino J, Franze M, Cincotti F, Varsta M, Heikkonen J, Babiloni F. A local neural classifier for the recognition of EEG patterns associated to mental tasks. IEEE Transactions on Neural Networks, 2002, 13(3): 678-686. DOI:10.1109/TNN.2002.1000132
16
Croft R J, Barry R J. Removal of ocular artifact from the EEG:a review. Neurophysiologie Clinique, 2000, 30(1): 5-19. DOI:10.1016/S0987-7053(00)00055-1
17
Zhou W D, Gotman J. Removing eye-movement artifacts from the EEG during the intracarotid amobarbital procedure. Epilepsia, 2005, 46(3): 409-414. DOI:10.1111/epi.2005.46.issue-3
18
Sarossy M G, Lee M H A, Bach M. A fast automated method for calculating the EOG arden ratio. Documenta Ophthalmologica, 2014, 128(3): 169-178. DOI:10.1007/s10633-014-9430-5
19
Wu Ming-Quan, Li Hai-Feng, Ma Lin. Automatic electrooculogram separation method for single channel electroencephalogram signals. Journal of Electronics and Information Technology, 2015, 37(2): 367-372.
( 吴明权, 李海峰, 马琳. 单通道脑电信号中眼电干扰的自动分离方法. 电子与信息学报, 2015, 37(2): 367-372. DOI:10.11999/JEIT140602)
20
James C J, Lowe D. Single channel analysis of electromagnetic brain signals through ICA in a dynamical systems framework. In:Proceedings of the 23rd Annual International Conference of the IEEE Engineering in Medicine and Biology Society. Istanbul, Turkey: IEEE, 2001, 1974-1977.
21
Davies M E, James C J. Source separation using single channel ICA. Signal Processing, 2007, 87(8): 1819-1832. DOI:10.1016/j.sigpro.2007.01.011
22
Mammone N, La Foresta F, Morabito F C. Automatic artifact rejection from multichannel scalp EEG by wavelet ICA. IEEE Sensors Journal, 2012, 12(3): 533-542. DOI:10.1109/JSEN.2011.2115236
23
Wu Z H, Huang N E. Ensemble empirical mode decomposition:a noise-assisted data analysis method. Advances in Adaptive Data Analysis, 2009, 1(1): 1-41. DOI:10.1142/S1793536909000047
24
Tang Y, Li J P, Wu H. A simple and accurate ICA algorithm for separating mixtures of up to four independent components. Acta Automatica Sinica, 2011, 37(7): 794-799. DOI:10.1016/S1874-1029(11)60210-3
25
Li H Q, Feng X L, Cao L, Li E B, Liang H A, Chen X L. A new ECG signal classification based on WPD and ApEn feature extraction. Circuits Systems and Processing, 2016, 35(1): 339-352. DOI:10.1007/s00034-015-0068-7
26
Li M A, Yang L B, Yang J F. Separation of EOG artifacts from EEG signals using Hilbert-Huang transform. In:Proceedings of the 2011 International Conference on Electric Information and Control Engineering. Wuhan, China: IEEE, 2011, 4453-4456.