眼电伪迹自动识别与去除的新方法

李明爱，郭硕达，田晓霞，杨金福，郝冬梅

(1.北京工业大学电子信息与控制工程学院，北京 100124； 2.计算智能与智能系统北京市重点实验室，北京 100124； 3.北京工业大学生命科学与生物工程学院，北京 100124)



眼电伪迹自动识别与去除的新方法

李明爱1,2，郭硕达1，田晓霞1，杨金福1,2，郝冬梅3

(1.北京工业大学电子信息与控制工程学院，北京 100124； 2.计算智能与智能系统北京市重点实验室，北京 100124； 3.北京工业大学生命科学与生物工程学院，北京 100124)

为了改善脑电中的眼电伪迹过估计问题及环境干扰耦合引起的非线性混合对眼电去除效果的影响，提出一种基于快速核独立成分分析(Fast Kernel Independent Component Analysis,FastKICA)与离散小波变换(Discrete Wavelet Transform,DWT)的眼电自动去除方法，即(Fast Kernel Independent Wavelet Transform ,FKIWT)方法.首先，利用FastKICA方法对脑电信号进行分离得到独立成分，并以相关系数为依据识别出眼电伪迹；进而，基于DWT对眼电伪迹进行多分辨率分析，将逼近分量置零，而细节分量保持不变，使得重构所得眼电伪迹成分保留更多有用脑电信号；最后，利用FastKICA逆变换重建眼电去除后的脑电信号.实验结果表明：FKIWT不仅有效改善了眼电过估计问题，增强了抗干扰能力和鲁棒性，而且在线性混合和非线性混合情况下，均得到较好的伪迹去除效果，特别是在非线性混合时优势更为明显，适合于实际在线应用.

非线性混合模型；快速核独立成分分析；离散小波变换；眼电过估计；鲁棒性

1　引言

脑电信号(Electroencephalogram,EEG)是通过电极在头皮表面采集到的反映大脑内部状态的生物电信号，在神经科学、脑科学、临床医学与康复工程等领域具有十分重要的作用.脑电信号采集过程中极易受到眼电、肌电、心电等多种噪声干扰[1]，而眼电信号由于幅值较大，严重影响脑电信号的分析和应用，如何有效去除脑电信号中的眼电伪迹尤为重要[2].

盲源分离(Blind Source Separation,BSS)方法是在源信号信息及混合过程都未知的情况下，仅需对观测信号进行处理就可实现对源信号和系统的辨识，近年来在信息处理领域备受青睐而发挥着越来越重要的作用[3].目前，用于去除脑电信号中眼电伪迹的常用BSS方法有：(1)基于矩阵联合对角化的预白化算法(Joint Approximate Diagonalization of Eigenmatrixes,JADE).JADE算法引入多变量数据的四维累积量矩阵，通过特征分解使其得以简化，并提高了结果的稳健性[4,5]；(2)PICA算法(Pearson Independent Component Analysis,PICA).该算法将ICA中的固定非线性对比函数方法与最大似然估计方法相结合，提高了对不同分布的源信号的分离能力，但计算量较大[6]；(3)快速独立向量分析(Fast Independent Component Analysis,FastICA)算法.该算法是在传统ICA算法的基础上发展起来的一种快速寻优迭代算法，因采用了定点迭代优化算法，具有较快的收敛速度[7,8]；(4)离散小波变换(DWT)和FastICA相结合的一种算法，记为DWICA算法[9].该算法的本质是在小波域上利用FastICA完成眼电去除，不仅收敛速度更快，而且增加了算法的稳健性.上述四种方法具有两个共同特点：(1)将分离或判断得到的眼电伪迹直接去除.因算法本身不能保证伪迹信号完全分离，而且噪声对分离效果也有一定影响，因而，分离出的眼电伪迹中依然包含有用的EEG成分，直接或全部去除会造成眼电伪迹的过估计[10]；(2)实际采集的EEG信号符合线性混合模型的假设[11,12].由于EEG信号作为一种生物电信号十分微弱，且易受外界环境多种干扰耦合的影响，实际采集到的EEG不可避免的会发生非线性畸变，因此，脑电和眼电信号的非线性混合更加符合实际，这使得上述BSS算法在眼迹去除问题中具有一定局限性，导致分离能力和稳定性有所下降.

核独立成分分析法(Kernel Independent Component Analysis,KICA)是近几年为解决BSS问题而提出的一种新方法，它利用核函数将信号从低维空间映射到高维空间，从而将非线性问题转换到高维空间的线性问题，对非线性混合问题有较好的分离能力和更好鲁棒性[13,14].FastKICA算法是在KICA的基础上发展而来的，具有更快的运算速度[15].

本文提出一种基于FastKICA和DWT算法的眼电伪迹自动识别与去除方法，即FKIWT方法.采用FastKICA算法对脑电混合信号进行盲源分离，并以相关系数识别出眼迹成分，进而对分离出的眼电伪迹基于DWT进一步分析和处理，再利用FastKICA逆算法重建眼迹去除后的脑电信号.实验研究表明，FKIWT方法在有效改善眼电过估计问题、增强方法自身抗干扰能力和鲁棒性等方面均得到较好效果.

2　基本原理

2.1核独立成分分析

2.1.1KICA原理

在源信号分量独立且最多只有一个分量服从高斯分布的前提下，KICA模型如式(1)所示.

X=HS

(1)

其中，X=[x1,x2,…,xn]T∈Rn×M表示观测信号，S=[s1,s2,…,sn]T∈Rn×M表示源信号，n为源信号和观测信号分量个数，M为每个分量的样本点数，H∈Rn×n表示未知混合矩阵.

KICA方法首先利用核函数K(xi,xj)(i,j∈{1,2,…,n}且i≠j)来代替向量xi和xj间的内积，以实现将观测信号从低维空间映射到满足可再生核希尔伯特空间(ReproducingKernelHilbertSpace，RKHS)特性的高维特征空间F；然后，以RKHS内的非线性函数作为对比函数，并运用对Gram矩阵的低阶近似等方法在RKHS内搜索对比函数的最小值[13]，以求解分离矩阵W=H-1；进而，求得对源信号的估计Y=WX.

对比函数最小值的求解等价于式(2)所示最小广义特征值的求解问题：

Kα=λFDα

(2)

(3)

显然，通过最小化C(W)，即可获得分离矩阵W.

2.1.2FastKICA工作原理

FastKICA算法是在KICA的基础上发展而来的，该算法利用Hilbert-Schmidt独立性判决准则(Hilbert-Schmidt Independence Criterion,HSIC)作为衡量变量统计独立性的对比函数，用牛顿类法对对比函数进行优化，通过极小化对比函数，获取分离矩阵，并采用不完全Cholesky分解方法来提高计算性能[15].

FastKICA的基本工作过程如下：

(1)确定观测信号X=[x1,x2,…,xn]T和核函数K(·,·).

(2)对观测信号进行去中心化和白化处理.

(3)基于HSIC准则定义对比函数C(W).

(4)对对比函数C(W)的Hessian矩阵进行不完全的Cholesky估计.

(5)用牛顿类法优化求解分离矩阵W.

2.2离散小波变换

小波变换是一种有效的时频分析方法，具有多分辨分析和对信号的自适应特点，广泛应用于非平稳信号分析中.

对任意离散函数f(t)∈L2(R)，其离散小波变换定义为：

(4)

式中，j,k分别代表频率分辨率和时间平移量，j,k∈Z,Ψj,k(t)为离散小波函数，满足

(5)

相应的DWT逆变换定义为

(6)

Mallat将计算机的多分辨率分析思想引入到小波分析中，统一了正交小波基的构造，提出了离散小波变换的快速分解和重构算法，即Mallat塔式分解算法，显著减少了DWT的运算数据量[16].

图1给出了Mallat算法的信号分解过程示意图.

(7)

其中，Aj,k、Dj,k分别表示j尺度空间的逼近分量和细节分量，且A0=f(t).

Mallat算法信号重构满足下式：

(8)

随着空间尺度逐级减小至1，完成信号的重构过程.

3　基于FKIWT的眼电伪迹自动去除方法

基于FKIWT进行眼电伪迹自动识别与去除，主要步骤如下：

(1)基于FastKICA算法的眼电伪迹分离

假设X(t)=[x1(t),x2(t),…,xn(t)]T∈Rn×M为n导观测信号，其中，M表示每导信号的样本点数，xi(t)(i=1,2,…,n-1)为n-1导脑电信号，xn(t)为眼电参考信号.利用FastKICA算法对X(t)进行核独立成分分析，得到分离矩阵W；进而，根据Y(t)=WX，(t)得到n个独立成分：

Y(t)=[y1(t),y2(t),…,yn(t)]T∈Rn×M

(9)

(2)利用相关系数识别眼电伪迹

相关系数用以描述两变量之间的相关性，其绝对值越大表明两个变量相似度越高.依式(10)计算每个独立成分yi(t)(i=1,2,…,n)与眼电参考信号xn(t)的相关系数：

(10)

Y(t)=[y1(t),…,yeog(t),…,yn(t)]T

(11)

(3)基于DWT的眼电伪迹分析与重构

(4)利用FastKICA逆变换重建眼电去除后的脑电信号

(12)

(13)

4　实验研究

4.1实验数据来源

本实验纯净脑电信号采用BCI Competition II的Data set Ⅲ数据库，该数据库记录了一位25岁的女性受试者左右手运动想象的EEG信号.该实验包含140组训练数据，其中左右手运动想象各70组.每组实验持续9s，实验开始后0～2s休息，2～3s出现“+”准备提示符，3～9s出现箭头(向左或向右)运动想象提示，实验时序如图3所示.

该实验利用G.tec脑电采集仪和Ag/AgCl电极采集了C3、Cz和C4三导的EEG数据，采样电极位置如图4所示，采集频率为128Hz，并对数据经过0.5～30Hz滤波.

眼电数据来源于BCI Competition IV的Data sets 2b中采集的垂直眼电参考信号.由于采集的眼电信号中不可避免的会混有脑电信号，为了还原纯净的眼电信号，本文根据眼电信号的频率范围，对眼电数据采用FIR低通滤器，进行0～10Hz的低通滤波，去掉其中的高频信号后，作为纯净眼电信号.

4.2基于线性混合模型的实验研究

本部分将基于线性混合模型，利用FKIWT方法从眼电过估计、伪迹去除及鲁棒性等多方面开展眼迹去除实验研究，以展示本文提出方法的有效性.实验环境为matlab2014a.

4.2.1线性混合模型

根据眼电和脑电间的双向传递特性，构建线性混合模型如下：

(14)

4.2.2性能评价指标

均方误差(Mean Squared Error,MSE)是普遍采用的一种伪迹去除评价指标，其计算式如下：

(15)

4.2.3眼电伪迹过估计实验

本节将从FastKICA算法和FKIWT方法的工作过程入手，通过计算每种方法去除掉的眼迹成分与纯净眼电信号的相关系数，检验本文方法改善眼电伪迹过估计能力和抗干扰能力.FastKICA算法中，核函数选用高斯径向基函数，核函数宽度取值为1，迭代精度设为0.0001，最大迭代次数设为10000.FKIWT方法中，FastKICA算法的参数选定同前；DWT选用coif小波基函数，小波分解层数为3.

图5给出任取一组脑电数据的眼电过估计实验结果.图中，EOGcle表示纯净眼电信号，yeog(t)为FastKICA算法分离和去除掉的眼电伪迹，ydel(t)则为基于FKIWT方法去除掉的眼电伪迹.

由图可见，yeog(t)依然包含较多的脑电信号，而ydel(t)与EOGcle的波形更为相似，说明FKIWT方法相对FastKICA算法而言，去除掉的眼电伪迹包含更少的脑电成分，有效改善了眼电过估计问题.

下面对140组脑电数据进行实验，并在线性混合模型中加入不同强度的白噪声，用来模拟脑电采集过程中受到的心电、肌电、出汗等其它干扰的影响，计算yeog(t)和ydel(t)分别与EOGcle间的相关系数reog和rdel.图6展现了从无噪声到有噪声，且噪声强度从-5dBw逐渐增强至10dBw时，基于140组数据所得平均相关系数的变化曲线.

从图6可知，随着噪声的增强，相关系数rdel变化很小，reog变化较大且减小趋势明显，说明本文FKIWT方法相对FastKICA算法具有更好的眼迹去除效果，不仅有效保留了有用的脑电信号，有利于减弱眼电过估计影响，而且具有更强的抗干扰能力.

4.2.4多种方法的眼电伪迹去除对比实验

本节将基于线性混合模型，以相关系数和均方误差为性能评价指标，将本文方法与其它常用方法进行比对实验研究，检验FKIWT方法的眼迹去除能力.

图7呈现了某组纯净的脑电信号和眼电信号，图8则进一步给出了依据式(14)进行线性混合并加入-5dBw的高斯白噪声的脑电信号和眼电信号.这里，影响因子k1、k2和k3设定为0.2到0.4之间的随机数，k4、k5和k6设定为0.02到0.3之间的随机数，高斯白噪声用以模拟采集过程中眼迹之外的其它干扰.由图清晰可见，在采样点300和800附近，眼电信号对脑电信号产生非常强烈的扰动.

利用本文FKIWT方法对图8所示混有眼迹和干扰的脑电信号进行处理，实验结果如图9.和图7对比可知，相应导联脑电信号波形非常接近，去噪效果良好.

进而，随机产生20个线性混合矩阵H，采用多种方法基于140组脑电数据各进行2800次眼迹去除实验.表1显示了基于JADE、PICA、FastICA和DWICA四种常用方法及FKIWT法去除眼电伪迹后的脑电信号与纯净脑电信号的平均相关系数(r)和均方误差(MSE).这里，FastICA及DWICA方法中的ICA均选用基于负熵判据的FastICA算法，其迭代精度和最大迭代次数与PICA及FKIWT中的FastKICA设置相同，DWICA中DWT与FKIWT中DWT参数选择相同.由表1可知，FKIWT方法得到三导脑电信号的平均MSE较JADE、PICA、FastICA和DWICA分别减少了46.1%、28.8%、13.4%和8.3%；而平均相关系数则分别提高了0.0029、0.0013、0.0005和0.0003.显然，FKIWT法相对其它方法在C3、Cz和C4三导脑电信号上均取得了最大相关系数和最小均方误差，效果优势明显.

表1　基于多种方法的眼迹去除实验结果

JADEPICAFastICADWICAFKIWTrMSErMSErMSErMSErMSEC30.99320.93150.99500.67280.99560.57060.99600.52220.99620.4805Cz0.99310.94650.99460.72030.99550.60130.99550.58650.99580.5389C40.99320.89440.99470.70440.99560.55180.99580.52000.99620.4740

4.2.5FKIWT方法鲁棒性实验

方法的鲁棒性对于其能否获得稳定的实验结果及在线应用十分重要.为此，在4.2.4节实验的基础上，进一步计算眼电伪迹去除后的脑电信号与纯净脑电信号的相关系数的方差(var)，以比较不同方法眼电分离效果波动的大小，从而评价各方法的鲁棒性.图10给出了五种方法在三导脑电信号上获得的平均实验结果.

从图10可知，FKIWT方法在三导脑电上均得到了最小相关系数方差，特别是相对于JADE和PICA方法优势尤为明显，说明FKIWT方法具有较好的鲁棒性，更适合应用于实时脑机接口系统中.

4.3基于非线性混合模型的实验研究

4.3.1非线性混合模型

脑电信号采集过程中易受外界环境多种干扰耦合的影响，导致实际采集到的EEG会不可避免地发生非线性畸变.为此，本节选用后置非线性混合模型作为脑电与眼电信号的混合模型进行实验研究，如图11所示.

4.3.2非线性畸变函数

考虑生物电信号的自身特性及干扰耦合的特点，本节选择以下三种非线性函数模拟非线性畸变：

(1)多项式函数f1(·)

(16)

(2)三角函数f2(·)

(17)

(3)双曲函数f3(·)

(18)

4.3.3基于多种方法的实验结果

下面将基于多种非线性畸变函数的非线性混合模型，针对FKIWT和其它常用方法的眼迹去除问题展开实验研究.为便于对比，实验环境及各方法中参数取值同4.2节.

为了简化非线性混合模型并排除随机线性混合矩阵H对分离结果的影响，实验中将H设定为固定值，影响因子k1到k6分别取值为：0.2、0.2、0.2、0.05、0.1和0.15；根据生物信号非线性畸变的特点，有关非线性函数系数分别取值为a=0.01，b=0.001，c=30，d=π/120，g=π/120，p=30.将140组脑电数据进行后置非线性混合，获得含有眼迹的脑电数据，再利用FKIWT和其他四种方法完成眼迹去除，以均方误差和相关系数为评价指标，平均实验结果如图12和图13所示.

从图可见，基于非线性混合模型，采用三种非线性畸变函数时，本文FKIWT法相对其它四种方法均取得了最佳实验效果.针对每种非线性函数情况，MSE均有所减小，相关系数都有所增大.其中，基于三种非线性函数在三导脑电信号上的平均MSE较JADE、PICA、FastICA和DWICA分别减少了36.1%、24.1%、20.4%和18.3%；而平均相关系数则分别提高了0.0051、0.0020、0.0008和0.0005.说明本文方法在非线性混合时的分离效果上较其他四种方法有较大的改善，且优于线性混合模型情况下的实验结果.

5　结论

本文提出一种FastKICA和DWT相结合的眼迹自动识别与去除方法，即FKIWT法.该方法的特色主要体现在从脑电与眼电信号时频特性的差别入手，利用DWT的多分辨特性仅对眼电伪迹成分进行处理与分析，通过置近似小波系数为零、细节小波系数不变，使得去除的眼电伪迹成分更加逼近纯净的眼电信号.基于国际标准数据库的大量实验研究表明，FKIWT法不仅有效减少了眼电过估计问题，具有较强的抗干扰能力和鲁棒性，而且在线性混合模型和非线性混合模型下，相对其它常用方法而言，FKIWT法均取得了最佳的眼迹去除效果，特别在非线性混合模型下优势更为明显，这为FKIWT法真正应用于脑电信号的实际在线分析与处理奠定了基础，另外，本研究对于更好地提示脑电混合方式和脑电建模具有重要理论意义.

[1]Akhtar M T,Mitsuhashi W,James C J.Employing spatially constrained ICA and wavelet denoising,for automatic removal of artifacts from multichannel EEG data[J].Signal Processing,2012,92(2):401-416.

[2]Vijila C,et al.Artifacts removal in EEG signal using adaptive neuro fuzzy inference system[A].2007 International Conference of Signal Processing,Communications and Networking[C].Chennai:IEEE,2006.589-591.

[3]Jung T P,Makeig S,Humphries C,et al.Removing electroencephalographic artifacts by blind source separation[J].Psychophysiology,2000,37(2):163-178.

[4]Cardoso J F,Souloumiac A.Blind beamforming for non-Gaussian signals[J].IEE Proceedings-F Padar and Signal Processing,1993,140(6):362-370.

[5]Pham D T,Cardoso J F.Blind separation of instantaneous mixtures of nonstationary sources[J].IEEE Transactions on Signal Processing,2001,49(9):1837-1848.

[6]Karvanen J,Koivunen V.Blind separation methods based on Pearson system and its extensions[J].Signal Processing,2002,82(4):663-673.

[7]Shen H,Kleinsteuber M,Huper K.Local convergence analysis of FastICA and related algorithms[J].IEEE Transactions on Neural Networks,2008,19(6):1022-1032.

[8]Dermoune A,Wei T W.FastICA algorithm:Five criteria for the optimal choice of the nonlinearity function[J].IEEE Transactions on Signal Processing,2013,61(8):2078-2087.

[9]李明爱,崔燕,杨金福.脑电信号中眼电伪迹自动去除方法的研究[J].电子学报,2013,41(06):1207-1213.

LI Ming-ai,CUI Yan,YANG Jin-fu.Research on removing ocular artifact automatically from EEG signals[J].Acta Electronica Sinica,2013,41(06):1207-1213.(in Chinese)

[10]Klados M A,Papadelis C L,Bamidis P D.REG-ICA:A new hybrid method for EOG artifact rejection[A].2009 9thInternational Conference on Information Technology and Applications in Biomedicine[C].Larnaka:IEEE,2009.278-281.

[11]李晓欧.基于独立分量分析和共同空间模式的脑电特征提取方法[J].生物医学工程学杂志,2010,27(06):1370-1374.

LI Xiao-ou.EEG feature extraction basedon ICA and CSP algorithms[J].Journal of Biomedical Engineering,2010,27(06):1370-1374.(in Chinese)

[12]万柏坤,等.ICA去除EEG中眼动伪差和工频干扰方法研究[J].电子学报,2003,31(10):1571-1574.

WANG Bai-kun,ZHU Xin,YANG Chun-mei,et al.Study on applying independent component analysis to remove blink artifacts and power noise in EEG[J].Acta Electronica Sinica,2003,31(10):1571-1574.(in Chinese)

[13]Marchesi S,Bruzzone L.ICA and kernel ICA for change detection in multispectral remote sensing images[A].2009 IEEE International Geoscience and Remote Sensing Symposium[C].Cape Town:IEEE,2009.1231-1234.

[14]Tateyama T,Nakao Z,Chen Y W.Classification of brain matters in MRI by kernel independent component analysis[A].2008 4thInternational Conference on Intelligent Information Hiding and Multimedia Signal Processing Proceedings[C].Harbin:IEEE Computer Soc,2008.713-716.

[15]Shen H,Jegelka S,Gretton A.Fast kernel-based independent component analysis[J].IEEE Transactions on Signal Processing,2009,57(9):3498-3511.

[16]罗志增,等.基于小波变换和盲信号分离的多通道肌电信号处理方法[J].电子学报,2009,37(04):823-827.

LUO Zhi-zeng,et al.A method of multi-channel EMG disposal based on wavelet transform and blind signal separation[J].Acta Electronica Sinica,2009,37(04):823-827.(in Chinese)

[17]王明祥,宁宇蓉,王晋国.基于Mallat算法的一维离散小波变换的实现[J].西北大学学报(自然科学版),2006,36(03):364-368.

WANG Ming-xiang,NING Yu-rong,WANG Jin-guo.Realization of one dimension discrete wavelet transform based on Mallat algorithm[J].Journal of Northwest University(Natural Science Edition),2006,36(03):364-368.(in Chinese)

李明爱(通信作者)女，2006年于北京工业大学获得博士学位,现为北京工业大学副教授,硕士生导师,主要研究方向为脑机接口技术与智能控制.

E-mail:limingai@bjut.edu.cn

郭硕达男,2012年于武汉纺织大学获得学士学位,现为北京工业大学控制科学与工程专业硕士研究生,主要研究方向为脑机接口技术、信息处理与模式识别.

A Novel Automatic Recognition and Removal Method of Ocular Artifacts

LI Ming-ai1,2,GUO Shuo-da1,TIAN Xiao-xia1,YANG Jin-fu1,2,HAO Dong-mei3

(1.CollegeofElectronicInformation&ControlEngineering,BeijingUniversityofTechnology,Beijing100124,China; 2.BeijingKeyLaboratoryofComputationalIntelligenceandIntelligentSystem,Beijing100124,China; 3.CollegeofLifeScienceandBio-Engineering,Beijing100124,China)

In order to improve the overestimation of ocular artifacts (OA) in electroencephalogram (EEG) and the OA removal effect of nonlinear mixture caused by environmental interference coupling,a novel automatic removal method is proposed based on fast kernel independent component analysis (FastKICA) and discrete wavelet transform (DWT),and it is denoted as FKIWT.The independent components are separated from the mixed EEG by using the FastKICA algorithm,and the correlation coefficient is applied to identify OA component;Then,the Multiresolution analysis of OA is achieved with DWT,the approximation wavelet coefficients are set to zero and the detail wavelet coefficients are not changed.So more useful EEG is remained in the reconstructed OA component;Furthermore,the clean EEG is restored with the inverse algorithm of FastKICA.The experimental results show that FKIWT can effectively improve the overestimation of OA and has perfect anti-interference ability and robustness.Meanwhile,the better effects of OA elimination are also obtained on the condition that the linear or nonlinear mixed model is adopted,and the latter’s advantage is especially obvious.The FKIWT is suitable for on-line application.

nonlinear mixed model;fast kernel independent component analysis;discrete wavelet transform;overestimation of ocular artifacts;robustness

2014-08-14;

2015-05-05 ；责任编辑：马兰英

国家自然科学基金(No.81471770，No.61201362);北京市自然科学基金(No.7132021，No.7132028)

R318

A

0372-2112 (2016)05-1032-08

电子学报URL:http://www.ejournal.org.cn10.3969/j.issn.0372-2112.2016.05.004