一种实用化在线脑电伪迹剔除方法与流程

本发明涉及一种实用化在线脑电伪迹剔除方法，属于生物医学信息处理技术领域。

背景技术：

脑电(electroencephalography,eeg)蕴含了大量的心理、生理和病理信息，目前被广泛运用于脑部疾病诊断和脑-机接口(brain-computerinterface,bci)等诸多研究领域中。正常脑电信号主要分布在δ频带(0-4hz)，θ频带(4-8hz)，α频带(8-13hz)，β频带(13-20hz)，γ频带(30-80hz)五个频带范围内，而eeg的幅度十分微弱，易受眼电(electrooculogram,eog)、肌电(electromyogram,emg)、50hz工频以及线性漂移等各种伪迹的干扰，这些伪迹存在严重影响了bci的性能。

以往的一些伪迹剔除方法中，伪迹拒绝无可避免使脑电受到eog、emg等伪迹的污染；线性滤波无法适用于伪迹与正常脑电信号的频带发生重叠的情况；线性回归的方法依赖于一个好的伪迹参考信号，通常我们很难找到一个适合于emg伪迹和非生理伪迹的参考信号。一些研究者直接使用带通滤波器或设定阈值的方法来剔除原始脑电中高频肌电、噪声和眼电信号，这种方法虽然简单，但是往往在剔除伪迹的同时也剔除了部分有用的脑电信号，导致预处理效果不理想。盲源分离(blindsignalseparation,bbs)技术是一种最有前景的脑电伪迹剔除方法，独立成分分析(independentcomponentanalysis,ica)就是bbs方法中的一种。

ica是信号处理研究中常用的一种统计方法，它旨在寻找能最大化数据中各个分量独立性的线性投影。它假设从电极阵列采集到的n维脑电信号x＝[x1,x2,…,xn]^t，该信号是由m个未知且统计独立的信号源s＝[s1,s2,…,sm]^t的线性组合。ica算法的目标是寻找分离矩阵w，使得：s＝wx。一些改进的ica算法已经被实现和免费使用，其中，fastica相比于其它改进的ica算法而言主要有以下优点：a、fastica算法收敛速度快，能够提高算法的实时性；b、该算法使用方便，数据选取时不需要设置步长参数；c、它具有类似于神经网络算法并行、分布的优点，计算复杂度低，占用内存空间小，适合于在线脑电信号的分析。通常，ica要求观测信号的数量要不少于源信号的数量，当观测信号的数量少于源信号数量时则会出现超完备ica现象。

目前，多通道脑电信号中伪迹的剔除算法已经较为成熟。而在单/少通道脑电信号的伪迹剔除中，由于脑电采集通道数量较少，所包含的有效信息较少，且缺乏伪迹参考信号，目前尚无十分有效的伪迹剔除方法。此外，多导联脑电信号设备由于价格昂贵，体型巨大，操作繁琐，往往不能满足便捷采集、处理脑电的要求。近年来，为了解决ica算法用于单/少通道脑电信号分析中存在的超完备问题，通道变换技术随着便携式脑电设备的研究而兴起，bogdanmijovic针对单通道脑电信号提出经验模态分解(empiricalmodedecomposition,emd)和ica相结合的脑电伪迹剔除方法，即emd-ica方法，而李明爱针对少通道脑电信号提出将离散小波变换(discretewavelettransform,dwt)和独立成分分析相结合的方法，即dwt-ica方法，该方法主要针对脑电信号中的眼电伪迹，而且没有在线运用，本发明主要是在该方法的基础上进行了改进，以提高伪迹成分自动识别和剔除能力，并进行在线运用。

小波变换(wavelettransform,wt)通过将信号分解成多尺度的小波函数来获取信号在时域和频率的局部特征，而dwt的计算速度快，广泛的应用于脑电这种非线性、非平稳信号的分析。传统的ica算法直接对含高斯噪声的观测信号进行分解时容易产生较大的误差，还增大了计算量，大幅降低了分离效果，为此，本发明在ica分离之前引入了dwt，dwt在将少通道脑电转化为适合ica分析的多通道脑电的同时增强了ica的分离效果。虽然ica能够成功的将多通道原始脑电信号中的真实脑电和伪迹成分分离，但是在这种方法中伪迹分量的识别往往依赖于专业人员的视觉检查，费时且带有很大的主观性，单纯的ica算法还远远达不到在线自动高效识别伪迹成分的要求。

在伪迹的在线自动识别方面，聚类技术是一种高效的方法，它通常基于每个独立分量的一些特定特征，用于自动从各个独立成分中分离出伪迹相关的成分。一些研究者基于两个独立成分之间的相似度运用了k均值算法和模糊c均值聚类算法进行聚类，自动识别并剔除伪迹。而这两种聚类算法都是一种迭代的方法，都需要事先给定聚类的目标数目来确定迭代的次数，实际聚类过程中聚类的目标数目往往是未知的。层次聚类(hierarchicalclustering,hc)算法相比于上述两种聚类算法而言主要有两大优点：首先，层次聚类过程中的树状图不仅包含各个类的组成信息，而且在每个类中还以节点高度的形成提供了每个类中元素的亲密信息；其次，整个聚类过程不需要事先设定聚类的数目。

基于运动想象(motorimagery,mi)的bci相对于基于稳态视觉诱发电位(stady-statevisualevokedpotential,ssvep)和p300电位的bci而言，它不依赖于屏幕的同步视觉刺激，可以实现异步bci，而且被试想象左右手运动意图可以通过分析被试初级运动区c3和c4导联的脑电信号加以识别，因而基于该类范式的bci更容易应用于便携式bci设备中。近年，新兴的便携式bci、便携式麻醉深度监测系统、便携式睡眠监测、便携式情绪状态识别等应用领域通常针对在线单/少通道的脑电信号进行分析。因此，研究实用的在线少通道运动想象脑电伪迹剔除算法对于便携式脑电设备的发展具有重要的意义。

技术实现要素：

本发明提供了一种结合dwt、ica和hc的实用化在线脑电伪迹剔除方法，即dwicac方法，以用于解决目前运动想象脑电预处理算法中，多种常规脑电伪迹无法在单通道的情况下，在线自动识别并剔除的问题。

本发明的技术方案是：一种实用化在线脑电伪迹剔除方法，所述方法的具体步骤如下：

a、对各通道的原始脑电信号进行降采样、工频陷波和线性漂移校正；

b、对校正之后各通道的脑电信号进行多尺度小波分解和小波系数单支重构；

c、分别对各通道的小波系数单支重构信号进行ica分解，计算各个分量的时域特征、频谱特征以及序间相似度特征；

d、对各通道的独立分量分别进行层次聚类，自动识别并剔除伪迹分量，并对剩余的独立分量进行fastica逆变换重构和小波系数重构，得到干净的脑电信号。

所述步骤a中的降采样、工频陷波和线性漂移校正具体方法如下：

a1、利用matlab信号处理工具箱中的downsample函数分别对c3、c4通道实时的原始脑电信号进行250hz降采样；

a2、利用matlab信号处理工具箱中的dlsim函数对c3、c4通道降采样之后的脑电数据进行50hz工频陷波处理；

a3、利用matlab信号处理工具箱中的detrend函数对陷波处理之后的信号进行线性漂移校正；

经过上述处理得到校正之后的脑电信号x(t)＝[x1(t),x2(t)，…，xm(t)]^t，其中：x(t)是一组m维随机向量，m为脑电采集通道数，t为转置运算，t＝1,2,…,n，n为脑电信号x(t)的采样点个数。

所述步骤b中，多尺度小波分解和小波系数单支重构的具体步骤如下：

b1、选择“db4”小波基，利用matlab软件离散小波变换工具箱中的wavedec函数对某一通道的信号x1(t)进行7层离散小波分解，得到信号x1(t)的逼近系数分量ca7和细节系数分量cd7、cd6、cd5、cd4、cd3、cd2、cd1；

b2、采用离散小波变换工具箱中的wrcoef函数分别对步骤b1中得到的逼近系数分量和细节系数分量进行单支重构，得到对应于b1中ca7,cd7，cd6,…,cd1的小波重构信号r(t)＝[r1(t),r2(t),…,r8(t)]^t。

所述步骤c中，ica分解和各个分量的时域特征、频谱特征和序间相似度计算方法如下：

c1、将步骤b2中的r(t)作为fastica的输入，利用matlab软件下fastica工具箱中的fastica函数对r(t)进行ica分解，得到分离矩阵w和8个独立分量(s1,s2,…,s8)；

c2、对步骤c1得到的8个独立分量计算峰度值作为各个独立分量的时域特征，式中，si表示第i个独立分量，i∈[1,8]，e()为统计期望函数，式中，l为独立分量si的数据量，为独立分量si中第l个数据点的m次方，和为独立分量i的4阶矩和2阶矩；

c3、计算步骤c1得到的8个独立分量在5个频带范围内的平均频带功率：fspectral＝[f(δ)f(θ)f(α)f(β)f(γ)]，式中，f()是平均频带功率；

c4、对步骤c1中8个独立分量s1,s2,…,s8设置滑动窗时间长度为l的通道独立分量数据作为一个子序列，总共得到nep个子序列，一个时序中通道的8独立分量最后就分成了nep个子序列的8×nep个独立分量，然后计算这8×nep个子序列分量的序间相似度，式中，nep是一个时序中的子序列个数，是分别从子序列j和k中计算出来的两个独立分量和之间的相关系数,式中，和分别为两个独立分量和中n个数据点的均值。

所述步骤d中，层次聚类、伪迹的自动识别和剔除、fastica逆变换和小波系数重构的具体方法如下：

d1、利用matlab软件信号处理工具箱中的zscore函数对特征向量f进行数据标准化；利用linkage函数计算特征向量中任意向量之间的相关距离，以距离最小准则为标准，利用linkage函数根据内平方距离法计算聚类树；

d2、根据步骤d1的聚类结果自动给8个独立分量s1,s2,…,s8贴上相应的类别标签，对伪迹分量进行置零处理，其余分量保持不变，得到剔除伪迹分量的独立分量矩阵z(t)；

d3、利用fastica逆变换将z(t)进行投影映射，得到小波变换系数u(t)＝w^-1z(t)，利用离散小波变换工具箱中的waverec函数对u(t)进行重构，得到干净的脑电信号。

本发明的有益效果是：

1、解决了现有运动想象脑电预处理方法需要训练和伪迹参考电极，以及无法在线实现的问题。

2、通过离散小波变换可以将少通道(c3和c4通道)的运动想象脑电信号转换为多通通信号用于后续fastica的分析。

3、小波变换系数比原始信号的超高斯性更强，峰度更大，从小波域进行ica在迭代算法的收敛速度、抗噪能力等方面具有显著优势。

4、利用层次聚类算法基于独立分量的时域、频谱和序间相关性特征，可高效的自动识别和剔除多种常规伪迹。

附图说明

图1为本发明的在线脑电伪迹剔除方法流程图；

图2为脑电采集电极的放置示意图；

图3为被试进行左右手运动想象任务的实验时序图；

图4为c3和c4通道原始脑电信号的时域图和频谱图，其中，(a)为c3和c4通道原始脑电信号的时域图，(b)为c3通道原始脑电信号的频谱图；

图5为校正之后c3和c4通道脑电信号的时域图和频谱图，其中，(a)为校正之后c3和c4通道脑电信号的时域图，(b)为校正之后c3通道脑电信号的频谱图；

图6为c3通道小波系数单支重构信号；

图7为c3通道ica分解后各个独立分量的波形图；

图8为c3通道ica分解后眼电分量的时域特征和频谱特征，其中，(a)为眼电分量的时域特征，(b)为眼电分量的频谱特征；

图9为c3通道ica分解后肌电分量的时域波形图和频谱特征，其中，(a)为肌电分量的时域波形图，(b)为肌电分量的频谱特征；

图10为c3通道各个分量进行层次聚类的树状图；

图11为预处理之后c3和c4通道的脑电信号波形图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

实施例1：如图1-11所示，一种实用化在线脑电伪迹剔除方法，实验设计和数据说明如下：

本文提出的预处理方法一共在10名健康被试的左右手运动想象脑电中进行了测试，其中被试为6男4女；均为右利手，年龄范围：19岁～25岁；本科及硕士学历，均无影响脑功能的病史和脑电采集经历，所有被试在实验前签署了实验研究知情同意书。

本次试验的数据采集设备为16导联eeg放大器(mipower-uc,eegcollectionv2，清华大学神经工程实验室)，0hz～250hz信号频带，采样频率为1000hz，24位a/d转换器，无工频陷波；采用国际标准导联10-20系统定制的16导联脑电帽(ag-agcl粉末电极，武汉格林泰克科技有限公司)。脑电采集电极的放置示意图如图2所示，实验主要采集运动功能区c3和c4两个电极的脑电信号，分别为x1(t)和x2(t)；左侧乳突m1作为记录参考电极，fpz为接地电极；实验之前要确保电极与被试头皮之间的阻抗低于5kω。

实验在安静的环境中进行，实验时序图如图3所示，实验过程中，被试坐在舒适的座椅上，眼睛距离电脑屏幕约1米，每次实验开始时刻会有蜂鸣提示，同时屏幕正中央出现“+”，提示被试本次实验已经开始，该状态持续2s；开始提示结束后“+”消失，同时屏幕正中央会随机出现任务图片提示，提示被试接下来进行哪种运动想象任务，该状态持续1.5s；提示图片消失之后屏幕正中央会出现“*”，提示被试开始进行之前图片所提示的左手/右手运动想象任务，这个任务时间持续6s，任务结束之后会进入休息状态，屏幕白屏，该状态持续时间为3s。该刺激范式由广泛应用于生物医学工程领域的刺激软件e-prime(版本version1.1)来完成。被试进行每次实验之间的休息时间为5分钟，每个被试进行5次实验，每次实验包括6个时序，10个被试的所有实验均在一天内完成。

本发明选取了被试5在第3次实验中某一次想象右手运动时的脑电数据用于后续分析。

对于上述脑电数据结果进行在线脑电伪迹剔除的方法包括如下具体步骤：

a、对c3、c4通道的原始脑电信号进行降采样、工频陷波和线性漂移校正；

b、对校正之后两个通道的脑电信号进行多尺度小波分解和小波系数单支重构；

c、分别对两个通道的小波系数单支重构信号进行ica分解，计算各个分量的时域特征、频谱特征以及序间相关性；

d、对两个通道的独立分量进行层次聚类，自动识别并剔除伪迹分量，并对剩余的独立分量进行fastica逆变换重构和小波系数重构，得到干净的脑电信号。

所述步骤a中的降采样、工频陷波和线性漂移校正具体方法如下：

a1、利用matlab信号处理工具箱中的downsample函数分别对c3、c4通道的原始脑电信号进行250hz降采样，以降低数据总量，提高数据处理的实时性；

降采样之后两个通道的时域波形图和频谱图如图4所示，可以看出两个通道的数据具有明显的数据漂移和50hz工频干扰现象；

a2、利用matlab信号处理工具箱中的dlsim函数对c3、c4通道降采样之后的脑电数据进行50hz工频陷波处理，剔除50hz工频干扰；

a3、利用matlab信号处理工具箱中的detrend函数对陷波处理之后的信号进行线性校正，消除数据线性漂移带来的伪迹；

经过a2和a3两步校正之后，两个通道的时域波形图和频谱图如图5所示，可以看出两个通道的数据已经去除了50hz工频干扰和数据漂移现象，但是数据中还存在严重的eog和emg等伪迹干扰；

经过上述处理得到校正之后的脑电信号x(t)＝[x1(t),x2(t)]^t∈r^m×n，其中：x(t)是一组m维随机向量，m为脑电采集通道数，这里m为2；t为转置运算，x1(t)和x2(t)分别代表c3和c4通道的脑电信号；t＝1,2,…,n，为脑电信号x(t)的n个采样点，这里n为1500；经过步骤a处理之后得到两个通道分别6s的脑电数据。以下分析均以c3通道为例，c4通道的处理方法和步骤同理。

所述步骤b中，多尺度小波分解和小波系数单支重构的具体步骤如下：

b1、选择“db4”小波基，利用matlab软件离散小波变换工具箱中的wavedec函数对c3通道的信号x1(t)进行7层离散小波分解，得到信号x1(t)的逼近系数分量ca7和细节系数分量cd1、cd2、cd3、cd4、cd5、cd6、cd7，此时，单通道的信号就变成了多通道的信号；

本发明选取了daubechies小波c3通道脑电信号的分解与重构。它通过最大频响平坦频率来刻画，不同阶的daubechies小波基被用在目前的一些研究中；然而，研究表明，本发明最终选取的4阶daubechies小波(db4)能够更好的正确表达eeg信号和尖峰。

b2、采用离散小波变换工具箱中的wrcoef函数对上一步分解得到的高低频小波系数进行单支重构，得到对应于b1中小波系数ca7,cd7，cd6,…,cd1的重构信号r(t)＝[r1(t),r2(t),…,r8(t)]^t，其中，8个小波重构信号与x1(t)的数据点数相同，其波形如图6所示。

所述步骤c中，ica分解和各个分量的时域特征、频谱特征以及序间相关性计算方法如下：

c1、将步骤b2中的小波重构信号r(t)作为fastica的输入，利用matlab软件下fastica工具箱中的fastica函数，通过不断迭代得到分离矩阵w，此处设定迭代精度为0.0001。当权值向量达到此精度时，就结束迭代过程，此时的权值即为最终的权值向量。经过迭代，此处求出的分离矩阵w为：

r(t)经过fastica分解后得到的8个独立分量(s1,s2,…,s8)，其波形如图7所示；

c2、分别对fastica分解后得到的8个独立分量计算峰度值作为各个独立分量的时域特征，式中，si表示第i个独立分量，e()统计期望函数,式中，l为独立分量si的数据量，这里l＝1500，为独立分量si中第l个数据点的m次方，当m为2和4时就得到和如果kurt(si)的值越大，则表示信号的动态分布越陡峭。各个独立分量的峰度值如图8(a)所示，从图中可以看出，s1和s2的峰度值排在前两位(kurt(s1)＝11.08,kurt(s2)＝8.878)，显著大于其它分量的峰度值，s3的峰度值为-0.7788，是8个独立分量中峰度值唯一为负值的分量，而以往的大量研究表明，脑电信号的伪迹成分中眼电伪迹往往具有最高的峰度值，噪声伪迹作为亚高斯信号，峰度值往往为负值。

c3、分别计算fastica分解后8个独立分量在5个频带范围内的平均频带功率：fspectral＝[f(δ)f(θ)f(α)f(β)f(γ)]，式中，f()是平均频带功率，它通过matlab软件中的pwelch函数来计算，将平均频带功率作为各个独立分量的频谱特征；其中，s1的频谱分布如图8(b)所示，从图中可以看出：s1的频率较低，主要分布在0-12hz范围内，且具有较高的幅值，s2也表现出和s1同样的频谱特征；而脑电信号的伪迹成分中只有眼电伪迹往往频率最低、幅值较大；s6的时域波形图和频谱分布分别如图8(a)和8(b)所示，从图8(b)中可以看出：s6的频谱主要分布在30-80hz的频带范围内，且在这一范围内具有相对较高的幅值，而脑电信号的伪迹成分中，由于头部转动、手部晃动等身体运动所导致的肌电伪迹往往频率最高(通常≥30hz)、幅值相对较大；

c4、在8个独立分量s1,s2,…,s8中设置滑动窗时间长度为1s(即250个采样点)的c3通道独立分量数据作为一个子序列，一个时序中c3通道的8独立分量最后就分成了6个子序列的48(6×8)个独立分量，然后计算这48个子序列分量的序间相似度，式中，nep是一个时序中的子序列个数，这里nep为6，是分别从子序列j和k中计算出来的两个独立分量(和)之间的相关系数,式中，和分别为两个独立分量和中n个数据点的均值，这里n为250；在这48个子序列分量中，当某些子序列分量和大多数子序列分量(这里设定为60％的子序列分量)都表现出极低相关系数时(其中：|·|为取绝对值运算)，取这些子序列分量中的最小相关系数作为原8个独立分量s1,s2,…,s8中对应独立分量的序间相似度；当某些子序列分量和大多数子序列分量(这里设定为60％的子序列分量)都表现出很高的相关系数时取这些子序列分量中的最大相关系数作为原8个独立分量s1,s2,…,s8中对应独立分量的序间相似度，从而得到一个时序中c3通道8个独立分量s1,s2,…,s8的序间相似度特征向量。

脑电采集过程中伪迹是随机出现的，难以预料，而且有时候在一些子序列中只会出现眼电伪迹或肌电等随机伪迹。而包含伪迹成分的子序列分量和其它不含伪迹的子序列分量往往没有相似的特征或相似度极低，另一方面，运动想象诱发的脑电子序列分量和其它有用脑电信号的子序列分量表现出极高的相似度。

通过上述计算，得到一个7维的特征向量f＝[f1,f2,f3,f4,f5,f6,f7]^t,其中f1为各个独立分量的时域特征，f2,f3,f4,f5,f6为频谱特征，f7为相似度特征；

所述步骤d中，对c1中得到的8个独立分量s1,s2,…,s8进行层次聚类，自动识别并剔除伪迹，然后依次对剩余的独立分量进行fastica逆变换重构和小波系数重构的具体方法如下：

d1、首先利用matlab软件信号处理工具箱中的zscore函数对步骤c中计算出的特征向量f进行数据标准化，使得特征向量矩阵的平均值为零，标准差为1；然后利用linkage函数计算特征向量中任意向量之间的相关距离，这里实际计算的是两个向量之间的欧氏距离，接下来以距离最小准则利用linkage函数根据内平方距离法计算聚类树，最后利用dendrogram函数根据linkage函数的输出绘制出如图10所示的层次聚类树状图；

d2、如图10所示，眼电分量s1和s2由于区别于其它独立分量的时域、频谱和序间相关性特征而被自动聚合为一类，根据以往的研究可以进一步看出：s1为眨眼伪迹，s2为眼动伪迹；肌电分量s6由于区别于其它独立分量的频谱和序间相似度特征而被单独聚合为一类；正常的脑电信号分量s4、s5、s7、s8首先聚合为一类，虽然噪声分量s3具有区别于独特的时域特征，但是它的总体特征和脑电相似，所以最后和正常脑电聚为了一类；聚类的同时给8个独立分量贴上相应的标签，对伪迹分量s1、s2、s3和s6进行置零处理，其余分量保持不变，即得到无伪迹成分的独立成分矩阵z(t)；

d3、利用fastica逆变换将z(t)进行投影映射，得到小波变换系数u(t)，即：

u(t)＝w^-1z(t)

选择小波基函数“db4”，利用离散小波变换工具箱中的waverec函数对步骤d2中剔除了各种伪迹的小波系数u(t)重构脑电信号，如图11所示，可以看出脑电信号中的eog伪迹、emg伪迹、50hz工频干扰、噪声干扰和数据漂移现象都已经被很好的去除了。

为了进一步验证本发明剔除伪迹的有益效果，分别将本发明中的方法和其它方法用于本实验时的性能进行了对比，10名被试采用李明爱提出的dwt-ica、bogdanmijovie提出的emd-ica和本发明提出的dwicac伪迹剔除方法并计算预处理之后信号的频带能量特征，最后基于支持向量机(supportvectormachine,svm)分类器的分类准确率如表1所示，svm是基于统计学习理论的机器学习方法，适合小样本分类问题，泛化能力强，能很好地适应运动想象脑电信号特征的单次试验分类。由于本发明只使用了c3和cz通道的脑电数据用于分析，分类时只提取了脑电信号的频带能量特征，特征维数相对较低，因此本发明最后选择了基于线性核的线性支持向量机，误差惩罚因子取值为1。

表1采用emd-ica、dwt-ica、dwicac伪迹剔除方法并基于svm分类器的分类准确率

从表1中可以发现，本发明提出的dwicac方法相对于emd-ica方法和dwt-ica伪迹剔除方法而言，几乎能自动识别并剔除脑电中的全部常规伪迹(eog、emg和噪声等)；而在分类准确率方面该方法也显著优于其它两种方法，其中：本发明对于左右手运动想象脑电的平均分类准确率为84.64％，分别高出emd-ica方法和dwt-ica方法8.23％和5.26％，最高分类准确率为89.27％，分别高出emd-ica方法和dwt-ica方法7.64％和3.55％；上述分析进一步体现了dwicac方法在少通道运动想象脑电在线伪迹剔除中的优越性。

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。