小波包分解结合异常值检测自动去除眼电中眨眼干扰的方法

王 田， 李玉榕， 陈建国， 陈东毅

(福州大学电气工程与自动化学院， 福建省医疗器械和医药技术重点实验室， 福建 福州 350116)

0 引言

眼睛是人体重要的视觉器官， 由眼球和眼的附属器官组成. 眼电图是将电极片放置在眼睛周围采集视网膜和角膜之间的电势差， 从而测量水平和垂直眼电信号的方法. 眼电图可以反映眼睛的运动情况. 如今眼电图已应用于脑机接口、 人机交互、 疲劳评估、 睡眠监测、 行为识别等研究领域. 例如， 文献提出用一种基于差分进化的混合特征选择技术， 利用EOG进行活动识别， 该特征选择技术可以识别眼动特征中最具眼动特征的子集， 它可以识别读书、 浏览网页、 写字、 看视频和复制等活动 . 文献用模糊学习的方法对眼电信号进行分类， 将分类的眼电信号作为轮椅的控制信号， 使轮椅能够做直线运动， 根据眼球运动强度的不同， 轮椅运动的距离也会有差距. 文献应用眼电图分析疲劳、 焦虑和任务负荷之间的关系， 从而判断受试者的状态.

无意识眨眼作为一种正常的生理反应， 其幅值较大且无法避免， 会使基于眼电的控制系统产生误操作， 例如： 文献等提出通过EOG信号识别阅读状态和非阅读状态， 但研究发现眨眼信号会影响阅读行为的识别率. 为了确保基于眼电信号的分析和控制系统的应用效果， 准确去除无意识眨眼是十分必要的. 因此， 文献提出利用ICA方法从原始多通道EOG信号中分离出眨眼信号并将此通道置零， 再将信号重构从而去除眨眼信号， 提高了阅读行为的识别率. 文献等提出离散小波变换与盲源分离结合的方法， 将分解后与眨眼信号相关系数最高的通道置零， 最后将信号进行重构得到去眨眼信号. 上述方法是将原始信号分解， 找到与眨眼信号相关性最大的通道直接置零， 这样会造成原始信号中部分有意义的信号损失. 为减少信号的损失， 文献[10-11]将原始信号在滑动窗口内求导从而确定无意识眨眼的位置， 但导数数值作为眨眼的判别比较粗糙， 并且滑窗的大小也会影响结果. 文献[12]提出用小波包将原始信号在频域上进行分解， 结合统计学理论确定阈值， 对小波包分解的低频分量采用阈值判定准则和重构策略去除眨眼信号， 但是由于样本数目有限， 基于统计信息确定的阈值有效性差， 无法准确确定眨眼区域.

在确保眨眼信号准确率的情况下， 为了减少去除无意识眨眼对有用信号造成的损失， 提出一种去除无意识眨眼的新方法. 该方法采用小波包分析结合异常值检测方法， 即首先将预处理后的信号进行小波包分解， 去除高频噪声， 然后用异常值检测方法准确确定无意识眨眼的起点和终点位置， 将该区间信号置零， 从而保证原始信号的损失尽可能小.

1 EOG信号的产生与获取

1.1 EOG信号的产生

医学研究表明， 眼球是一个双极性的球体， 其中角膜呈现正电位， 视网膜呈现负电位. 当眼球转动时角膜和视网膜之间会存在电势差. 这个电势差称为眼电信号[13]. 眼电信号随眼球转动角度而变化. 将眼电信号的幅值反映在时间轴上构成一条连续的曲线， 称为眼电图.

眼球运动大致分为三类： 扫视， 静止和眨眼. 扫视运动是指目光迅速向四周移动掠过. 眨眼信号是在不自觉的情况下， 眼睑自主的张开和闭合的过程， 眼睑的张开和闭合会带动眼球运动， 从而形成眨眼信号的EOG波形曲线. 无意识眨眼信号是常见的人体生理活动， 人每分钟无意识眨眼次数一般为12～19次.

1.2 数据采集

本次实验共有7名(4男3女)受试者， 年龄均在24～26岁之间， 且眼部健康未患有任何眼部疾病. 信号采集设备选用美国Texas Instruments的ADS1299脑电采集模块. 通常眼电信号的采集方式有两种， 分别为单级导联方式和双极导联方式. 其中单极导联方式采集到的是混合信号， 所以需要高性能的采集记录系统， 而双极导联方式可以克服单极导联方式的缺点， 两电极接触的信号一样多， 不会因为异常信号而影响差值， 所以本文采用双极导联方式[14]. 采用5个电极片， 电极片的位置如图1(a)所示. 其中ch1， ch2测量水平眼电信号， ch3， ch4测量垂直眼电信号， ch5为参考电极放于乳突处.

图1 电极片位置图以及眼电信号采集模块图Fig.1 Electrode plate position and EOG aquistion module

为避免受试者在实验过程中受到外界环境干扰， 本实验设置在一个安静的房间进行， 同时为了排除眼睛疲劳的影响， 实验不宜过长， 每组实验时间设置为20 s， 每人采集8组实验， 共采集56组数据， 每次实验不规定受试者眼睛运动的方向， 受试者听到提示音后保持头部及身体不动， 眼睛做垂直扫视或者水平扫视运动， 在实验期间受试者可以自主进行无意识眨眼， 不规定受试者眨眼时刻以及眨眼次数. 实验实物图如图1(b)～1(d)所示.

2 小波包结合异常值检测方法

采用小波包结合异常值检测进行无意识眨眼的检测. 首先将采集的原始信号进行预处理. 将预处理后的垂直眼电信号进行小波包分解滤除高频噪声， 然后基于异常值检测的方法准确确定无意识眨眼的起点和终点， 最后将此区域置零. 算法流程图如图2所示.

图2 WPT-OD去除眨眼算法流程图Fig.2 Flow chart of blink signal removing algorithm based on WPT-OD

2.1 信号预处理

将采集到的原始眼电信号进行信号预处理， 首先将垂直两通道的信号和水平两通道的信号分别进行差分计算得到垂直和水平眼电信号， 同时去除基线漂移. 眼电图信号是微伏级别的低频复杂信号， 硬件采集过程中易受工频50 Hz的干扰， 因此采用巴特沃斯滤波去除工频干扰.

2.2 小波包分解

小波包分析是在小波分析的基础上提出的， 它弥补了小波分解在高频段的频率分辨率较差， 而在低频段的时间分辨率较差的缺点. 小波包分析不仅能够多层次地划分信号频带， 而且可依据信号的特点自适应地选择与之相对应的频带， 使信号的频带与频谱相匹配， 从而提高信号在时域与频域的分辨率[15-16].

在小波多分辨率分析的基础上， 把尺度函数由φ(t)改为u0(t)， 把小波函数由Ф(t)改为u1(t)， 于是就有如下的两个尺度方程：

(1)

由式(1)定义的函数集合un(t)，n∈Z， 为由u0(t)=Ф(t)所确定的小波包， 其中g(k)=(-1)kh(1-k)，h(k)，g(k)相互正交， 分别是多分辨率分析中的高通滤波系数和低通滤波系数.

设x(t)为L2(R)空间函数， 对其离散采样序列x(p)，p=1, 2, …,N， 小波包分解算法为：

(2)

式(2)体现了小波包分解的过程， 其本质上是利用一组低通、 高通滤波器组合成共轭正交滤波器h，g在不同频带上对信号进行分解.

由于EOG信号频率较低， 所以本文采用小波包分解滤波得到低频分量. 实验采样频率为1 000 Hz， 将原始信号进行7层小波包分解， 根据Nyquist原理知第7层的前三个分量的频率主要处于11.7 Hz以下， 由于眼电信号的主要频率处于10 Hz以下， 7层小波分解后的前三个分量刚好包含所有有用的眼电信号. 所以将小波包分解的前三个分量进行重构作为异常值检测的原始信号.

2.3 异常值检测确定眨眼信号的起点和终点的方法

修正箱线图法是一种基于中位数的运算方法. 首先计算信号的中值， 并依据从小到大的顺序将信号分为四个部分， 每个部分称为四分位数.A1为最小值与信号中值之间的中位数，A2为信号的中值，A3为信号中值与最大值之间的中位数，A4为A3与A1的差值. MC是信号的偏度估计， 通过A1，A2，A3，A4和MC可以确定信号的阈值, 上限(upper limit)用UL表示， 下限(lower limit)用LL表示.

如果信号中的值高于上限或低于下限则标记为异常值. 而本文中由于将信号进行小波包重构后的信号幅值都大于零， 所以只需确定上限即可.

为了确定连续眼电信号中所有无意识眨眼信号的起点和终点， 本文用以下步骤标记眨眼起点和终点的位置， 流程图如图3所示.

首先使用前文提到的异常值检测方法确定阈值. 找到重构信号中的所有局部极大值及局部极大值的位置. 取出所有信号的局部极大值Mn(局部极大值的总数为N)， 用相邻局部极大值进行相减(ΔM=Mn+1-Mn)并取绝对值. 将相邻的差值绝对值与所求阈值相比较. 如果ΔM为正且大于阈值， 则标记为眨眼信号的起点. 如果ΔM为负且绝对值大于阈值， 则标记为眨眼信号的终点. 这些步骤应用于剩余的ΔM值， 并且将满足条件的区域标记起点或终点. 确定出起点和终点区域之后， 将该区域置零， 从而去除眨眼信号.

图3 起点与终点确定过程流程图Fig.3 The flow chart of starting and ending point determine process

3 实验结果与分析

3.1 实验结果

将实验采集的一例数据进行预处理后结果如图4所示， 图中0～1.6 s眼睛处于静止状态， 1.6 s左右是一次无意识眨眼， 2～6 s眼睛处于静止状态， 6～7 s水平扫视， 9 s左右是无意识眨眼， 12～14 s是垂直扫视， 15和18 s左右是无意识眨眼， 18.5～20 s是水平扫视. 从图4中可以看出垂直眼电信号中的眨眼信号十分明显， 图5展示了7层小波包分解后的前三个分量的波形以及前三个分量重构后的波形， 利用重构后的信号进行异常值检测. 分别采用CC， PC和修正箱线图法确定眨眼信号的阈值， 利用图3流程图确定所有眨眼信号的起点和终点. 取重构的前10 s眼电信号分别用CC、 PC和修正箱线图法进行检测， 并对比三种方法的检测结果， 结果如图6所示. 图6中黑色实线代表EOG信号， 红色虚线、 蓝色短线、 绿色实线分别代表CC， PC和调整框图法的检测结果. 结果表明， 在9 s处幅值较大的眨眼信号， 三种WPT-OD方法都能检测眨眼信号的起点和终点. 但在1.5 s左右处幅值相对较小的眨眼信号， 只有修正箱线图法可以检测出来, 可见三种方法中修正箱线图法的效果相对较好. 用修正箱线图法检测20 s的实验数据结果如图7所示， 该方法能够准确检测所有眨眼信号.

图4 原始水平和垂直眼电信号图Fig.4 Horizontal and vertical EOG

图5 7层小波包分解前三个分量以及重构信号结果图Fig.5 The first three components of the 7 layers WPT and reconstruction of signal result

图6 CC， PC和调整框图法结果对比图Fig.6 CC, PC and adjusted box plot method comparison chart

图7 调整框图法无意识眨眼信号检测结果图Fig.7 Adjusted box plot method unconscious blink signal detection result

3.2 性能分析

对WPT-OD方法中的三种异常值检测准则进行对比分析， 为验证不同准则对眨眼信号的检测效果， 将实验得到的56组实验数据进行处理， 每组实验数据都为20 s. 采用三种指标分析算法的性能： WPT-OD方法中三种标准对眨眼检测的准确率(accuracy， AC)； 相关系数(correlation coefficient， CC)表示去除眨眼后的信号与原始信号之间的相关性； 去眨眼信号与原始信号之间的损失率(loss).

式中：i代表检测出的眨眼信号；I代表眨眼信号的总数.

WPT-OD方法三种准则的性能与文献[11]中提出的小波包结合统计参数算法的性能对比见表1, 包含各种方法的准确率及其标准差(standard deviations ， SD)， 相关性及其标准差和损失率及其标准差的结果. 如表1所示， 三种WPT-OD的准确率、 相关性和损失率均优于小波包和统计参数算法， 且基于调整框图法的三种指标均达到最优. 表2对比滑窗为1 s和2 s时， WPT-OD方法检测眨眼信号的准确率.

表1 WPT-OD， WPT统计法及MSDW方法性能对比

表2 WPT-OD方法在不同窗口下的准确率对比

4 结语

为去除眼电信号中的无意识眨眼信号， 并且使原始信号损失尽可能少， 提高基于眼电信号控制的人机交互系统的正确率， 提出一种去除眼电信号中无意识眨眼干扰的方法. 实验结果验证了算法的有效性， 此方法去除眼电信号中的无意识眨眼信号， 可为后续眼电信号的研究做基础. 但此算法确定出无意识眨眼信号区间后， 将该区间置零， 后续可以研究合适的插值方法将该区间重构去除无意识眨眼， 从而进一步减少损失率， 提高相关性. 另外， 目前脑机接口去除眼电伪迹是常见的研究内容， 由于眼电伪迹中存在大量的无意识眨眼， 因此后续也可以将此方法推广应用到脑机接口的研究中.