CN102835955B - 一种无需设定阈值的脑电信号中眼电伪迹自动去除的方法 - Google Patents

Abstract

本发明提出一种无需人工设定阈值的脑电信号中眼电伪迹自动去除的方法，属于生物信息技术领域，主要应用于脑电信号的预处理过程中。具体包括：对采集的含有眼电伪迹的脑电信号进行独立分量分解；并且求取每个独立分量的峭度、序列renyi熵和样本熵作为特征向量，进而使用k均值聚类分析的方法自动识别出含有眼电伪迹的独立分量，并将其置零，其余分量不变，对信号进行重构，得到纯净的脑电信号。本发明解决了传统的眼电伪迹去除过程中需要人工对伪迹进行识别、费时费力、工作量大的问题，并且本方法无需人工设定阈值就可以实现自动识别并去除眼电伪迹的目的，弥补了以往方法中设定阈值时需要研究人员具备一定的先验知识、主观性强的不足。

Description

一种无需设定阈值的脑电信号中眼电伪迹自动去除的方法

技术领域：

本发明涉及生物信息技术领域，特别涉及脑电信号（Electroenc-

ephalogram，EEG）的预处理技术。具体是涉及一种无需人工设定阈值的脑电信号中眼电伪迹（Ocular Artifact，OA）的自动去除技术。

技术背景：

脑电信号是一种微弱的电生理信号，其中包含了大量的生理与疾病信息，在临床医学、认知心理学等研究领域中发挥着重要作用。头皮表面记录脑电的方式是目前脑电信号研究中所采用的最主要的信号采集方式，该方式记录技术简单、易操作，但信号易受噪声的干扰。脑电干扰给脑电信号的解释与分析造成了很大的困难。因此，有效识别并去除脑电信号中的伪迹，得到纯净的脑电信号，对进一步研究大脑的真实活动具有重要意义。

影响脑电信号的伪迹主要包括眼电伪迹、心电伪迹、肌电伪迹、50Hz工频干扰等。相对于其它伪迹成分而言，眼电伪迹是使EEG信号失真的一种最主要的干扰成分。眼电伪迹的幅值较大，最大可达到100mV，并且眼电信号源距离测量电极较近，对脑电信号造成了很大的干扰。眼电的产生是由于跨膜间存在静息电位，角膜侧为正，巩膜侧为负，这可以看作是有一个电偶极子存在。眼球转动过程中，该偶极子改变了眼睛周围的电场，从而影响头皮电场，眼球运动方向一侧的电极记录到眼电伪迹为正，眼球运动相反方向一侧的电极记录得到的眼电伪迹为负，由于大脑是电的良导体，眼电信号从前额处向后传播，遍历整个头部，放置在头部用于采集脑电信号的电极可以采集到这些信号，造成对脑电信号的干扰。而且，眼电伪迹的频带与脑电信号的频带相混叠，很难使用传统的滤波方法进行去除，因此，脑电信号预处理阶段最重要的工作就是如何有效地去除眼电伪迹。

从脑电信号中有效地去除眼电伪迹一直是生物医学信号处理领域的研究热点，其中一些主要的研究成果有：

1实验控制法：这种方法要求研究人员给予受试者一定的实验指导，例如要求受试者在实验过程中始终保持闭眼状态或尽量避免眨眼和眼球转动以减少眼电伪迹的干扰，这样做虽然可以在一定程度上抑制干扰，但是前者会引入α波的干扰，后者在长时间的实验过程中极易引起疲劳，并且不易做到。

2伪迹剔除法：它的主要思想是将受到眼电伪迹影响脑电信号剔除以获得较纯净的脑电信号，但是这样做的同时剔除了大量的脑电信号，造成了信息的遗漏，并且人工剔除干扰信号的劳动量过大，主观性强，而且不适应在线脑机接口系统的需求。

3回归方法：研究表明眼电伪迹在测量电极处的衰减程度与电极到眼部的距离有关，此方法利用这一原理通过计算得到眼电信号在每个头皮电极处的衰减因子，然后将采集到的每导脑电信号减去乘以衰减因子的眼电信号，从而得到伪迹去除后的数据。它的模型为

$\mathrm{EEG}\left(\mathrm{\ϵ}\right)={\mathrm{EEG}}_s\left(\mathrm{\ϵ}\right)-\mathrm{\ηEOG}---\left(1\right)$

式中，EOG表示测量得到的眼电信号；EEG_s(ε)表示原始脑电信号；EEG(ε)表示伪迹去除后的脑电信号；ε表示相应的测量通道；η表示比例系数，即眼电信号在第ε导处的衰减因子。η值常用回归分析的方法进行估计。这种伪迹去除方法的优点是直观易懂，物理意义明确，缺点是没有考虑到脑电信号的双向传播性，即实际测得的眼电信号中也包含一些脑电成分。因此，回归方法去除伪迹的同时不可避免的去除某些脑电成分。

4小波变换法：小波变换(Wavelet Transform,  WT)是数学领域的一个重要分支，近年来越来越受到各个学科领域的极大重视。由于其在时频两域都具有表征信号局部特征的能力和多分辨率分析的特点，因此被誉为“数学显微镜”。小波变换的基本思想是将原始信号通过伸缩和平移后，分解为一系列具有不同时间分辨率、不同频率特性的子带信号，这些子带信号具有良好的时域、频域等局部特征，进而可以实现对信号在时域和频域的局部化分析，从而克服了傅里叶分析在处理非平稳信号和复杂图像时所存在的局限性。

但是使用小波变换法去除伪迹要求信号和噪声的频带不能混叠。而脑电信号与眼电伪迹的频带相混叠，于是研究者开始将传统的小波变换法与其它去除伪迹的方法相结合，取长补短，期望获得更好的效果。

5主成分分析法（Principal Components Analysis, PCA）：主成分分析是一种多通道的信号处理方法，其基本思想是利用正交原理将一组相关变量转化成为另一组相互独立的变量(即主成分)，分解出来的主成分按大小排序，通常前面几个主成分就包括了信号80%的信息，因此PCA也被用来数据降维，去除冗余。PCA可以把多导联的脑电信号分解为相互正交的成分，去掉噪声和伪迹成分，再对其余成分进行重构，就可以得到较为纯净的脑电信号。

PCA的主要缺点在于：1信号和噪音波形相似时，难以从EEG中去除伪迹；2它只对正交的信号源比较有效，在去除伪迹时只能去除与脑电信号相正交的干扰成分；3它只适合于分解低阶相关信号，难以分解高阶的相关信号。

6独立分量分析(Independent Component Analysis, ICA)： 独立分量分析理论的发展可追溯到上个世纪80年代初期，法国学者J.Hernah和C.Juttne等人首次提出了ICA分析的基本概念，最初被用来解决“鸡尾酒会”问题，随着近年来在ICA方面研究兴趣的增加，该理论已被广泛应用在图像处理、天线阵列、控制科学、生物医学信号处理等领域，其研究热潮方兴未艾。

ICA是实现盲源分离 (Blind Source Separation, BSS)最主要的方法之一，是一种基于高阶统计量的分析方法，处理对象为非高斯信号。理论上可以认为脑电信号、眼电伪迹、心电伪迹、肌电伪迹以及其他干扰源所产生的干扰信号都是由相互独立的源产生的，而且可以将观测信号视为是伪迹与脑电信号的线性组合，并且伪迹源的数目通常比头皮上测量的脑电通道数要少，因而独立成分分析方法适用于分离脑电信号中的各种伪迹。在实际应用中，ICA方法可以很好的分离频带相重叠的信号，弥补了小波分析去除伪迹的不足。但是ICA方法也有缺陷：在分离出信号和伪迹后，需要人工识别并筛选出伪迹成分再进行去除，费时费力，且具有很大的主观性。

发明内容

本发明针对传统脑电信号中眼电伪迹去除的过程中需要人工对眼电伪迹进行识别、费时费力、工作量大的缺陷，提出了一种眼电伪迹自动去除的方法。

其主要特征包括：首先利用独立分量分析对含有眼电伪迹的脑电信号进行分解，得到若干个独立分量，然后求取每个独立分量的峭度、序列renyi熵和样本熵作为区分脑电信号和眼电伪迹成分的特征，进而采用k均值聚类分析——一种无监督的分类算法进行特征分类，通过将含有眼电伪迹的独立分量置零重构出纯净脑电信号。

具体步骤如下：

1独立分量分解

本发明采用1999年芬兰学者Aapo Hyvarinen等人提出的基于负熵的固定点ICA算法，简称fastica算法，进行独立分量分解。设n通道观测信号x(t)=[x₁(t),x₂(t),…,x_n(t)]，经过ICA分解得到m个独立分量y(t)=[y₁(t),y₂(t),…,y_m(t)]。

2特征提取

经过ICA分解得到m个独立分量y(t)=[y₁(t),y_i(t),…,y_m(t)]，分别求取它们的峭度、序列renyi熵和样本熵作为分类特征。

2.1 峭度：峭度是时间序列的四阶累计量，可以表征信号的波动情况。眼电伪迹的峭度值要显著大于脑电信号，因此，利用该参数作为识别眼电伪迹的特征。

2.2 renyi熵：renyi熵是一种反映信号随机性的参数，眼电伪迹的renyi熵值明显小于脑电成分,因而使用其作为区分眼电伪迹和脑电成分的特征。

2.3 样本熵：样本熵是一种测量时间序列复杂度的方法，而且眼电伪迹的样本熵值比脑电成分低，因此使用样本熵作为特征区分眼电伪迹和脑电成分。

3特征分类：本发明使用k均值聚类分析的方法进行分类，这是一种动态的、无监督的分类方法。区分眼电伪迹和脑电成分显然是一个两分类问题，因此先随机取两个样本作为聚类中心，计算每个样本到这两个中心的距离，把样本归入最近的聚类中心，然后用本类所有样本的均值代替原有的聚类中心的值，这样就得到了新的聚类中心，再重新计算各样本到新的聚类中心的距离，重新分类、修改中心点，直到新的聚类中心与上一次的中心点一致时，算法停止，然后再按照最小距离原则将独立分量分为两类，算法流程如图2所示。

4信号的重构

将含有眼电伪迹的独立分量置零，其它独立分量不变，对信号进行重构，得到纯净的脑电成分。

具体如下：

一种无需人工设定阈值的脑电信号中眼电伪迹自动去除的方法，其特征在于：首先，对采集的脑电信号进行独立分量分解，得到m个独立分量；其次，对于每一个独立分量，求取它们的峭度、序列renyi熵和样本熵作为特征向量；然后，使用k均值聚类分析——一种无监督的分类算法自动识别出含有眼电伪迹的独立分量；最后，将含有眼电伪迹的独立分量置零，其余独立分量不变，对信号进行重构，得到纯净的脑电信号；包括以下步骤：

1）独立分量分解

使用fastICA算法对信号进行分解，它的基本模型表达式为：

$x\left(t\right)=\mathrm{As}\left(t\right)---\left(1\right)$

式中，x(t)=[x₁(t),x₂(t),…,x_n(t)]^T∈R^n×N是n通道的观测向量，为测量得到，1≤t≤N，N表示采样点总数，s(t)∈R^m×N为源信号，混合矩阵R^n×m，独立分量分析的目标是求解分离矩阵w，通过w从观测信号x(t)中恢复出未知的源信号s(t)，源信号s(t)使用下式估计得到:

$y\left(t\right)=\mathrm{wx}\left(t\right)---\left(2\right)$

式中，y(t)是源信号s(t)的估计，y(t)=[y₁(t),y_i(t),…,y_m(t)]^T∈R^m×N是求得的m个独立分量，其中y_i(t)是第i个独立分量，1≤i≤m，i为整数，t表示第t个采样点，分离矩阵w∈R^m×n；

2）特征提取

2.1 峭度：对于每一个独立分量y_i(t)，求取它的峭度；

$k=m_4-3m_2^2---\left(3\right)$

式中，k为独立分量y_i(t)的峭度值，1≤i≤m，i为整数；m_β为y_i(t)的β阶中心矩，m_β=E{(y_i(t)-m₁)^β}，β={2,4}，E(·)表示对括号里的表达式求均值，m₁为y_i(t)的均值；

2.2 renyi熵：对于每一个独立分量y_i(t)，其计算公式为

$H\left(y_i\left(t\right)\right)=\frac{1}{1-\mathrm{\α}}\log \frac{1}{N^2}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{h=1}{\overset{N}{\mathrm{\Σ}}}G(y_i\left(j\right)-y_i\left(h\right),2{\mathrm{\σ}}^2)---\left(4\right)$

式中，

，θ=y_i(j)-y_i(h)；方差σ=0.2～0.3，α=1.8，N为独立分量y_i(t)的采样点个数；y_i(j)与y_i(h)是独立分量y_i(t)在采样点j和h的采样值；

2.3 样本熵：对于每一个独立分量y_i(t)，其样本熵的计算过程如下：

第一步：设独立分量y_i(t)有N个采样点，t表示第t个采样点，1≤t≤N，t为整数；

第二步：从独立分量y_i(t)中提取b维矢量，b=2，即：

$Y_b\left(t\right)=[y_i\left(t\right),y_i(t+1)]---\left(5\right)$

式中，y_i(t)和y_i(t+1)是独立分量y_i(t)在采样点t和t+1的采样值，1≤t≤N-1，t为整数；

第三步：使用第二步的方法得到两个b维矢量为Y_b(t)和Y_b(t₁)，定义它们之间的距离为：

$d[Y_b\left(t\right),Y_b\left(t_1\right)]=\underset{c=\left\{\mathrm{0,1}\right\}}{\max }\left[\right|y_i(t+c)-y_i(t_1+c)\left|\right]---\left(6\right)$

式中，1≤t,t₁≤N-1,t≠t₁，t、t₁为整数，Y_b(t)和Y_b(t₁)为公式（5）中Y_b(t)在采样点t和t₁的值；y_i(t+c)和y_i(t₁+c)为独立分量y_i(t)在采样点t+c和t₁+c的采样值，c∈{0,1}，求取差值的最大值，即为两个b维矢量的距离；

第四步：设定一个阈值r，r=(0.1~0.25)std(y_i(t))，std(y_i(t))为y_i(t)的标准差；对于每一个t值，统计d[Y_b(t),Y_b(t₁)]小于r的数目，记为T，1≤t,t₁≤N-1,t≠t₁，t、t₁为整数，t₁可以取N-2个值，因此，总的距离数为N-2，计算T与总的距离数N-2的比值，记为

，即

$B_t^b\left(r\right)=\frac{T}{N-2}---\left(7\right)$

第五步：计算

的均值，即

$B^b\left(r\right)=\frac{1}{N-1}\underset{t=1}{\overset{N-1}{\mathrm{\Σ}}}{B}_t^b\left(r\right)---\left(8\right)$

第六步：再从同一个独立分量y_i(t)中提取b₁维矢量，b₁=3，即：

$Y_{b_1}\left(t\right)=[y_i\left(t\right),y_i(t+1),y_i(t+2)]---\left(9\right)$

式中，y_i(t)、y_i(t+1)和y_i(t+2)是独立分量y_i(t)在采样点t、t+1和t+2的采样值，1≤t≤N-2，t为整数；

第七步：使用第六步的方法得到两个b₁维矢量为

和

，定义它们之间的距离为：

$d[Y_{b_1}\left(t\right),Y_{b_1}\left(t_1\right)]=\underset{e=\left\{\mathrm{0,1,2}\right\}}{\max }\left[\right|y_i(t+e)-y_i(t_1+e)\left|\right]---\left(10\right)$

式中，1≤t,t₁≤N-2,t≠t₁，t、t₁为整数，

和

为公式（9）中

在采样点t、t₁的采样值；y_i(t+e)和y_i(t₁+e)为独立分量y_i(t)在采样点t+e和t₁+e的采样值，e∈{0,1,2}，求取差值的最大值，即为两个b₁维矢量的距离；

第八步：在第四步中已计算得到阈值r，对于每一个t值，统计

小于r的数目，记为T₁，1≤t,t₁≤N-2,t≠t₁，t、t₁为整数，t₁可以取N-3个值，因此总的距离数为N-3，计算T₁与总的距离数N-3的比值，记为

，即

$B_t^{b_1}\left(r\right)=\frac{T_1}{N-3}---\left(11\right)$

第九步：计算

的均值

$B^{b_1}\left(r\right)=\frac{1}{N-2}\underset{t=1}{\overset{N-2}{\mathrm{\Σ}}}{B}_t^{b_1}\left(r\right)---\left(12\right)$

第十步：

样本熵的计算公式为

$\mathrm{sampEn}=-\ln [B^{b_1}\left(r\right)/B^b\left(r\right)]---\left(13\right)$

3）特征分类

经过上面的计算，对于每个独立分量y_i(t) {1≤i≤m}，可计算得三个特征量，将其依序构成如下特征向量，即：

$f_i={\left[\begin{array}{ccc}{f}_{1i}& f_{2i}& f_{3i}\end{array}\right]}^T---\left(14\right)$

式中，f_1i、f_2i、f_3i分别表示y_i(t)的峭度、renyi熵和样本熵；分别针对每一个特征进行一次k均值聚类分析，共进行三次分类，以f_1i为例，给出计算步骤:

3.1给定样本种类，选定初始聚类中心；显然，这是一个两分类问题，初始聚类中心选为f₁₁,f₁₂两个分量，即center₁(0)=f₁₁,center₂(0)=f₁₂； center_q(0)表示第q类初始聚类中心，q=1或2，第p次迭代计算的聚类中心用center_q(p)表示；

3.2将样本按最小距离原则分为两个类，即在第p次迭代时，若||f_1i-center₁(p)||＜||f_1i-center₂(p)||，则f_1i∈c₁(p)，否则f_1i∈c₂(p)；||·||表示对其中的两个向量求取它们之间的欧式距离；c_q(p)表示在第p次迭代计算中第q类的聚类域，f_1i∈c_q(p)表示将样本f_1i划分到第q类，q=1或2；

3.3重新计算聚类中心；以本类所有样本的均值代替原有的聚类中心值；即

，n₁为第q类样本总数；

3.4若两次迭代计算的聚类中心相同，即center_q(p+1)=center_q(p)，则计算完毕，否则令p=p+1，返回3.2；

4信号的重构

眼电伪迹的峭度值明显高于脑电信号，而renyi熵和样本熵的数值低于脑电信号，因此，综合三次分类的结果，将所有独立分量分为c₁和c₂两类，设c₁类表示含有眼电伪迹的独立分量，c₂类表示脑电独立分量；将c₁类中所有成分置零，c₂不变，得到去除眼电伪迹后的独立分量

，利用式（15）进行重构，即可得到纯净的脑电成分

；

$\tilde{x}=w^{-1}\tilde{y}---\left(15\right)$

式中，w^-1为第1部分求得的分离矩阵w的伪逆矩阵。

本方法无需人工设定阈值就可以达到自动识别并去除眼电伪迹的目的，改善了以往方法中设定阈值时需要研究人员具备一定的先验知识、主观性强的不足附图说明

图1 脑电信号中眼电伪迹自动去除的算法流程图

图2 K均值聚类分析算法流程图

图3 国际10-20系统标准电极放置图

图4 C3导联眼电伪迹去除效果图

图5 Cz导联眼电伪迹去除效果图

图6 C4导联眼电伪迹去除效果图

具体实施方式

本发明提出一种无需人工设定阈值的脑电信号中眼电伪迹的自动去除方法，完整的算法流程包括下面四（1.2.3.4）个部分。其中，第1部分是现有方法，本发明申请的特征包括三(2.3.4)个部分：

完整的算法流程如图1所示。

1独立分量分解

本发明采用fastICA算法，该算法有如下优点：①收敛速度比批处理和自适应处理都快；②负熵作为高斯性度量的效果优于累计量；③采用牛顿迭代法，收敛有保证；

它的基本模型表达式为：

$x\left(t\right)=\mathrm{As}\left(t\right)---\left(2\right)$

式中，x(t)=[x₁(t),x₂(t),…,x_n(t)]^T∈R^n×N是n通道的观测向量，为测量得到，1≤t≤N，N表示采样点总数，s(t)∈R^m×N为源信号，混合矩阵A∈R^n×m，独立分量分析的目标是求解分离矩阵w，通过w从观测信号x(t)中恢复出未知的源信号s(t)，源信号s(t)使用下式估计得到:

$y\left(t\right)=\mathrm{wx}\left(t\right)---\left(3\right)$

式中，y(t)是源信号s(t)的估计，y(t)=[y₁(t),y_i(t),…,y_m(t)]^T∈R^m×N是求得的m个独立分量，其中y_i(t)是第i个独立分量，1≤i≤m，i为整数，t表示第t个采样点，分离矩阵w∈R^m×n。

2特征提取

2.1 峭度：对于每一个独立分量y_i(t)，求取它的峭度。

$k=m_4-3m_2^2---\left(4\right)$

式中，k为独立分量y_i(t)的峭度值，1≤i≤m，i为整数；m_β为y_i(t)的β阶中心矩，m_β=E{(y_i(t)-m₁)^β}，β={2,4}，E(·)表示对括号里的表达式求均值，m₁为y_i(t)的均值。

2.2 renyi熵：对于每一个独立分量y_i(t)，其计算公式为

$H\left(y_i\left(t\right)\right)=\frac{1}{1-\mathrm{\α}}\log \frac{1}{N^2}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{h=1}{\overset{N}{\mathrm{\Σ}}}G(y_i\left(j\right)-y_i\left(h\right),2{\mathrm{\σ}}^2)---\left(5\right)$

式中，

，θ=y_i(j)-y_i(h)；方差σ=0.2～0.3，α=1.8，N为独立分量y_i(t)的采样点个数；y_i(j)与y_i(h)是独立分量y_i(t)在采样点j和h的采样值。

2.3 样本熵：对于每一个独立分量y_i(t)，其样本熵的计算过程如下：

第一步：设独立分量y_i(t)有N个采样点，t表示第t个采样点，1≤t≤N, t为整数。

第二步：从独立分量y_i(t)中的提取b维矢量，b=2，即：

$Y_b\left(t\right)=[y_i\left(t\right),y_i(t+1)]---\left(6\right)$

式中，y_i(t)和y_i(t+1)是独立分量y_i(t)在采样点t和t+1的采样值，1≤t≤N-1，t为整数。

第三步：使用第二步的方法得到两个b维矢量为Y_b(t)和Y_b(t₁)，定义它们之间的距离为：

$d[Y_b\left(t\right),Y_b\left(t_1\right)]=\underset{c=\left\{\mathrm{0,1}\right\}}{\max }\left[\right|y_i(t+c)-y_i(t_1+c)\left|\right]---\left(7\right)$

式中，1≤t,t₁≤N-1,t≠t₁，t、t₁为整数，Y_b(t)和Y_b(t₁)为公式（6）中Y_b(t)在采样点t和t₁的值。y_i(t+c)和y_i(t₁+c)为独立分量y_i(t)在采样点t+c和t₁+c的采样值，c∈{0,1}，求取差值的最大值，即为两个b维矢量的距离。

第四步：设定一个阈值r，r=(0.1~0.25)std(y_i(t))，std(y_i(t))为y_i(t)的标准差。对于每一个t值，统计d[Y_b(t),Y_b(t₁)]小于r的数目，记为T，1≤t,t₁≤N-1,t≠t₁，t、t₁为整数，t₁可以取N-2个值，因此，总的距离数为N-2，计算T与总的距离数N-2的比值，记为，即

$B_t^b\left(r\right)=\frac{T}{N-2}---\left(8\right)$

第五步：计算

的均值，即

$B^b\left(r\right)=\frac{1}{N-1}\underset{t=1}{\overset{N-1}{\mathrm{\Σ}}}{B}_t^b\left(r\right)---\left(9\right)$

第六步：再从同一个独立分量y_i(t)中的提取b₁维矢量，b₁=3，即：

$Y_{b_1}\left(t\right)=[y_i\left(t\right),y_i(t+1),y_i(t+2)]---\left(10\right)$

式中，y_i(t)、y_i(t+1)和y_i(t+2)是独立分量y_i(t)在采样点t、t+1和t+2的采样值，1≤t≤N-2，t为整数。

第七步：使用第六步的方法得到两个b₁维矢量为

和

，定义它们之间的距离为：

$d[Y_{b_1}\left(t\right),Y_{b_1}\left(t_1\right)]=\underset{e=\left\{\mathrm{0,1,2}\right\}}{\max }\left[\right|y_i(t+e)-y_i(t_1+e)\left|\right]---\left(10\right)$

式中，1≤t,t₁≤N-2,t≠t₁，t、t₁为整数，

和

为公式（10）中在采样点t、t₁的采样值。y_i(t+e)和y_i(t₁+e)为独立分量y_i(t)在采样点t+e和t₁+e的采样值，e∈{0,1,2}，求取差值的最大值，即为两个b₁维矢量的距离。

第八步：在第四步中已计算得到阈值r，对于每一个t值，统计小于r的数目，记为T₁，1≤t,t₁≤N-2,t≠t₁，t、t₁为整数，t₁可以取N-3个值，因此总的距离数为N-3，计算T₁与总的距离数N-3的比值，记为

，即

$B_t^{b_1}\left(r\right)=\frac{T_1}{N-3}---\left(11\right)$

第九步：计算

的均值

$B^{b_1}\left(r\right)=\frac{1}{N-2}\underset{t=1}{\overset{N-2}{\mathrm{\Σ}}}{B}_t^{b_1}\left(r\right)---\left(12\right)$

第十步：

样本熵的计算公式为

$f_i={\left[\begin{array}{ccc}{f}_{1i}& f_{2i}& f_{3i}\end{array}\right]}^T---\left(14\right)$

3特征分类

经过上面的计算，对于每个独立分量y_i(t) {1≤i≤m}，可计算得到三个特征量，将其依序构成如下特征向量，即：

$f_i={\left[\begin{array}{ccc}{f}_{1i}& f_{2i}& f_{3i}\end{array}\right]}^T---\left(15\right)$

式中，f_1i、f_2i、f_3i分别表示y_i(t)的峭度、renyi熵和样本熵。分别针对每一个特征进行k均值聚类分析，共进行三次分类，以f_1i为例，给出计算步骤:

3.1给定样本种类，选定初始聚类中心；显然，这是一个两分类问题，初始聚类中心选为f₁₁,f₁₂两个分量。即center₁(0)=f₁₁,center₂(0)=f₁₂；center_q(0)； center_q(0)表示第q类初始聚类中心，q=1或2，第p次迭代计算的聚类中心用center_q(p)表示。

3.2将样本按最小距离原则分为两个类，即在第p次迭代时，若||f_1i-center₁(p)||＜||f_1i-center₂(p)||，则f_1i∈c₁(p)，否则f_1i∈c₂(p)。||·||表示对其中的两个向量求取它们之间的欧式距离。c_q(p)表示在第p次迭代计算中第q类的聚类域，f_1i∈c_q(p)表示将样本f_1i划分到第q类，q=1或2。

3.3重新计算聚类中心；以本类所有样本的均值代替原有的聚类中心值。即，n₁为第q类样本总数。

3.4若两次迭代计算的聚类中心相同，即center_q(p+1)=center_q(p)，则计算完毕，否则令p=p+1，返回3.2。

4信号的重构

眼电伪迹的峭度值要明显高于脑电信号，而renyi熵和样本熵的数值要低于脑电信号，因此，综合三次分类的结果，将所有独立分量分为c₁和c₂两类，设c₁类表示含有眼电伪迹的独立分量，c₂类表示脑电独立分量。将c₁类中所有成分置零，c₂不变，得到去除眼电伪迹后的独立分量

，使用下式进行信号重构

$\tilde{x}=w^{-1}\tilde{y}---\left(16\right)$

式中，

为重构后纯净的脑电信号，w^-1为第1部分求得的分离矩阵w的伪逆矩阵。

为了给出详细具体的实施方式和步骤，将使用本方法对国际BCI（Brain Computer Interface, BCI）竞赛数据进行处理，给出实验结果并验证本方法去眼电伪迹的有效性。

数据介绍

脑电数据来源于第四届国际BCI竞赛数据集中的Data sets 2b数据包。采集的脑电数据经过0.5Hz~100Hz的带通滤波，50Hz的工频陷波，采样频率250Hz，幅值范围为±100μV。电极按照国际标准导联10－20系统安放在C3、Cz、C4三个导联，如图3所示。另外，还同步采集了3导眼电信号。

实验步骤：

1.对数据进行ICA分解：输入数据如表1所示。

表1 原始数据表

得到的6个独立分量，数据如表2所示。

表2 独立分量数据表

2、以第一个独立分量y₁(t)为例，分别计算它的峭度、renyi熵和样本熵。

根据公式（4），m₁=E(y₁(t))=-0.0176,m₂=E{(y₁(t)-m₁)²}=1.000,

m₄=E{(y₁(t)-m₁)⁴}=9.0068,则y₁(t)峭度值为

。

将已知参数带入公式（5），得到y₁(t)的renyi熵的计算公式为

$H\left(y_1\left(t\right)\right)=\frac{1}{1-1.8}\log \frac{1}{{3001}^2}\underset{j=1}{\overset{3001}{\mathrm{\Σ}}}\underset{h=1}{\overset{3001}{\mathrm{\Σ}}}\frac{1}{0.25\sqrt{2\mathrm{\π}}}{e}^{-\frac{{(y_1\left(j\right)-y_1\left(h\right))}^2}{2{\left(0.25\right)}^2}}=1.4117$

式中，y₁(j)与y₁(h)是独立分量y₁(t)在采样时刻j和h的值。

根据第2部分中样本熵的计算步骤可以得到独立分量y₁(t)的样本熵为

$\mathrm{sampEn}=1.7003$

分别求取每个独立分量的峭度、renyi熵和样本熵作为特征向量，数据如表3所示。

表3 各独立分量的特征提取结果

独立分量	峭度	Renyi熵	样本熵
				1	6.00678780713594	1.41174521697880	1.70030048041082
2	-0.588623024270334	1.44455029109777	0.52937952495178
				3	0.31548101801454	1.57832435990357	1.54817199374019
4	0.44292110594564	1.57436808397370	1.75928182468544
				5	0.39335875263220	1.56190706120367	0.96911146049469
6	0.24823401414372	1.58037659833919	1.50219022964535

3、眼电伪迹的识别：利用k均值聚类分析进行特征分类。结果如表4所示。

表4 各独立分量的分类结果

4、眼电伪迹的峭度值要明显高于脑电信号，而renyi熵和样本熵的数值要低于脑电信号。因此将1、2、5三个独立分量判别为眼电伪迹分量并将其置零，利用式（16）进行重构。眼电伪迹去除结果如图4-6所示。

本发明结合峭度、renyi熵和样本熵三个分别表征序列峰值、随机性和复杂度的参数作为特征向量，为说明这样提取特征向量的优势，下面分别使用单一特征和三种特征对眼电伪迹进行去除，得到4段去除眼电伪迹后的脑电信号，然后分别计算这4段信号与一段该受试者在实验过程中无眼电伪迹干扰的纯净脑电信号之间的均方误差，结果如表5所示。从表5可以得出，三种特征相结合去除眼电伪迹后得到的脑电信号与纯净脑电信号之间的均方误差最小，波形最为接近，这说明三种特征相结合去除眼电伪迹最为有效，并且对脑电信号破坏较小。

表5 不同特征去除伪迹后的均方误差对比表

	峭度	Renyi熵	样本熵	三种特征
					C3	60.0824114780632	15.9543462159107	16.4989572483297	13.6996491011524
Cz	156.878955654143	84.3315031246680	53.1042782036112	50.5866649065173
					C4	60.1610176410670	22.2264340506635	19.6937866400084	18.0135610303600

将本发明与希尔伯特—黄变换(Hilbert—Huang Transform，HHT)去除眼电伪迹的方法做一个对比，进一步说明本发明的有效性和优势。分别使用本方法和HHT法去除同一段脑电信号中的眼电伪迹，得到两段处理后的脑电信号，然后再分别计算这两段信号与一段该受试者在实验过程中无眼电伪迹干扰的纯净脑电信号之间的均方误差，结果如表6所示。从表6中可以看出使用本方法得到的均方误差要明显低于HHT法，证明使用本方法去除眼电伪迹后的脑电信号波形与纯净的脑电信号波形比较接近，即本方法可以更有效地去除眼电伪迹成分，并且对脑电成分破坏较小。

表6 两种方法均方误差对比

	C3	Cz	C4
				原始信号	62.8809	159.3931	61.8404
HHT法	19.6196	64.3851	18.7603
				本方法	13.6996	50.5867	18.0136

Claims (1)

1.一种无需人工设定阈值的脑电信号中眼电伪迹自动去除的方法，其特征在于：首先，对采集的脑电信号进行独立分量分解，得到m个独立分量；其次，对于每一个独立分量，求取它们的峭度、序列renyi熵和样本熵作为特征向量；然后，使用k均值聚类分析——一种无监督的分类算法自动识别出含有眼电伪迹的独立分量；最后，将含有眼电伪迹的独立分量置零，其余独立分量不变，对信号进行重构，得到纯净的脑电信号；包括以下步骤：

1）独立分量分解

使用fastICA算法对信号进行分解，它的基本模型表达式为：

x(t)=As(t)     （1）

式中，x(t)=[x₁(t),x₂(t),…,x_n(t)]^T∈R^n×N是n通道的观测向量，为测量得到，1≤t≤N，N表示采样点总数，s(t)∈R^m×N为源信号，混合矩阵A∈R^n×m，独立分量分析的目标是求解分离矩阵w，通过w从观测信号x(t)中恢复出未知的源信号s(t)，源信号s(t)使用下式估计得到:

y(t)=wx(t)     （2）

式中，y(t)是源信号s(t)的估计，y(t)=[y₁(t),…,y_i(t),…,y_m(t)]^T∈R^m×N是求得的m个独立分量，其中y_i(t)是第i个独立分量，1≤i≤m，i为整数，t表示第t个采样点，分离矩阵w∈R^m×n；

2）特征提取

2.1峭度：对于每一个独立分量y_i(t)，求取它的峭度；

$k=m_4-3m_2^2---\left(3\right)$

式中，k为独立分量y_i(t)的峭度值，1≤i≤m，i为整数；m_β为y_i(t)的β阶中心矩，m_β=E{(y_i(t)-m₁)^β}，β={2,4}，E(·)表示对括号里的表达式求均值，m₁为y_i(t)的均值；

2.2renyi熵：对于每一个独立分量y_i(t)，其计算公式为

$H\left(y_i\left(t\right)\right)=\frac{1}{1-\mathrm{\α}}\log \frac{1}{N^2}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{h=1}{\overset{N}{\mathrm{\Σ}}}G(y_i\left(j\right)-y_i\left(h\right),2{\mathrm{\σ}}^2)---\left(4\right)$

式中，

θ=y_i(j)-y_i(h)；方差σ=0.2～0.3，α=1.8，N为独立分量y_i(t)的采样点个数；y_i(j)与y_i(h)是独立分量y_i(t)在采样点j和h的采样值；

2.3样本熵：对于每一个独立分量y_i(t)，其样本熵的计算过程如下：

第一步：设独立分量y_i(t)有N个采样点，t表示第t个采样点，1≤t≤N，t为整数；

第二步：从独立分量y_i(t)中提取b维矢量，b=2，即：

Y_b(t)=[y_i(t),y_i(t+1)]      （5）

式中，y_i(t)和y_i(t+1)是独立分量y_i(t)在采样点t和t+1的采样值，1≤t≤N-1，t为整数；

第三步：使用第二步的方法得到两个b维矢量为Y_b(t)和Y_b(t₁)，定义它们之间的距离为：

$d[Y_b\left(t\right),Y_b\left(t_1\right)]=\underset{c=\left\{\mathrm{0,1}\right\}}{\max }\left[\right|y_i(t+c)-y_i(t_1+c)\left|\right]---\left(6\right)$

式中，1≤t≤N-1，1≤t₁≤N-1，t≠t₁，t、t₁为整数，Y_b(t)和Y_b(t₁)为公式（5）中Y_b(t)在采样点t和t₁的值；y_i(t+c)和y_i(t₁+c)为独立分量y_i(t)在采样点t+c和t₁+c的采样值，c∈{0,1}，求取差值的最大值，即为两个b维矢量的距离；

第四步：设定一个阈值r，r=(0.1～0.25)std(y_i(t))，std(y_i(t))为y_i(t)的标准差；对于每一个t值，统计d[Y_b(t),Y_b(t₁)]小于r的数目，记为T，1≤t≤N-1，1≤t₁≤N-1，t≠t₁，t、t₁为整数，t₁可以取N-2个值，因此，总的距离数为N-2，计算T与总的距离数N-2的比值，记为

即

$B_t^b\left(r\right)=\frac{T}{N-2}---\left(7\right)$

第五步：计算

的均值，即

$B^b\left(r\right)=\frac{1}{N-1}\underset{t=1}{\overset{N-1}{\mathrm{\Σ}}}{B}_t^b\left(r\right)---\left(8\right)$

第六步：再从同一个独立分量y_i(t)中提取b₁维矢量，b₁=3，即：

$Y_{b_1}\left(t\right)=[y_i\left(t\right),y_i(t+1),y_i(t+2)]---\left(9\right)$

式中，y_i(t)、y_i(t+1)和y_i(t+2)是独立分量y_i(t)在采样点t、t+1和t+2的采样值，1≤t≤N-2，t为整数；

第七步：使用第六步的方法得到两个b₁维矢量为

和

定义它们之间的距离为：

$d[Y_{b_1}\left(t\right),Y_{b_1}\left(t_1\right)]=\underset{e=\left\{\mathrm{0,1,2}\right\}}{\max }\left[\right|y_i(t+e)-y_i(t_1+e)\left|\right]---\left(10\right)$

式中，1≤t≤N-2，1≤t₁≤N-2，t≠t₁，t、t₁为整数，

和

为公式（9）中在采样点t、t₁的采样值；y_i(t+e)和y_i(t₁+e)为独立分量y_i(t)在采样点t+e和t₁+e的采样值，e∈{0,1,2}，求取差值的最大值，即为两个b₁维矢量的距离；

第八步：在第四步中已计算得到阈值r，对于每一个t值，统计

小于r的数目，记为T₁，1≤t≤N-2，1≤t₁≤N-2，t≠t₁，t、t₁为整数，t₁可以取N-3个值，因此总的距离数为N-3，计算T₁与总的距离数N-3的比值，记为即

$B_t^{b_1}\left(r\right)=\frac{T_1}{N-3}---\left(11\right)$

第九步：计算

的均值

$B^{b_1}\left(r\right)=\frac{1}{N-2}\underset{t=1}{\overset{N-2}{\mathrm{\Σ}}}{B}_t^{b_1}\left(r\right)---\left(12\right)$

第十步：

样本熵的计算公式为

$\mathrm{sampEn}=-\ln [B^{b_1}\left(r\right)/B^b\left(r\right)]---\left(13\right)$

3）特征分类

经过上面的计算，对于每个独立分量y_i(t){1≤i≤m}，可计算得三个特征量，将其依序构成如下特征向量，即：

f_i=[f_1i f_2i f_3i]^T     （14）

式中，f_1i、f_2i、f_3i分别表示y_i(t)的峭度、renyi熵和样本熵；分别针对每一个特征进行一次k均值聚类分析，共进行三次分类，以f_1i为例，给出计算步骤:

3.1给定样本种类，选定初始聚类中心；显然，这是一个两分类问题，初始聚类中心选为f₁₁,f₁₂两个分量，即center₁(0)=f₁₁,center₂(0)=f₁₂；center_q(0)表示第q类初始聚类中心，q=1或2，第p次迭代计算的聚类中心用center_q(p)表示；

3.2将样本按最小距离原则分为两个类，即在第p次迭代时，若||f_1i-center₁(p)||<||f_1i-center₂(p)||，则f_1i∈c₁(p)，否则f_1i∈c₂(p)；||·||表示对其中的两个向量求取它们之间的欧式距离；c_q(p)表示在第p次迭代计算中第q类的聚类域，f_1i∈c_q(p)表示将样本f_1i划分到第q类，q=1或2；3.3重新计算聚类中心；以本类所有样本的均值代替原有的聚类中心值；即 ${\mathrm{center}}_q(p+1)=\frac{1}{n_1}\underset{f_{1i}\&Element;c_q\left(p\right)}{\mathrm{\&Sigma;}}{f}_{1i},$ n₁为第q类样本总数；

3.4若两次迭代计算的聚类中心相同，即center_q(p+1)=center_q(p)，则计算完毕，否则令p=p+1，返回3.2；

4）信号重构

眼电伪迹的峭度值明显高于脑电信号，而renyi熵和样本熵的数值低于脑电信号，因此，综合三次分类的结果，将所有独立分量分为c₁和c₂两类，设c₁类表示含有眼电伪迹的独立分量，c₂类表示脑电独立分量；将c₁类中所有成分置零，c₂不变，得到去除眼电伪迹后的独立分量利用式（15）进行重构，即可得到纯净的脑电成分

$\tilde{x}=w^{-1}\tilde{y}---\left(15\right)$

式中，w^-1为第1部分求得的分离矩阵w的伪逆矩阵。

Images