少通道下运动想象脑电特征的提取方法

摘要

本发明提供了一种少通道下运动想象脑电特征的提取方法，该方法利用每个通道时间序列的多个坐标延迟来同时优化空间滤波器与一个高阶多参数FIR滤波器，并用非参数化方法估计特征量与类别间信息熵的方法来选取最优个数的时间坐标延迟，基于多个时间坐标延迟的高阶FIR滤波器比一个参数的FIR滤波器有更好的幅频特性，该方法提供了在少通道情况下提取想象运动类别特征的有效途径。

说明书

技术领域

本发明涉及生物信息技术领域，尤其涉及一种少通道下运动想象脑电特征 的提取方法。

背景技术

目前，由于病变或意外造成的大脑与神经肌肉通路损伤的疾病在现有的医 疗条件下还无法治愈。如在中国大约有六百万脑瘫患者，十二岁以下的儿童就 占有一百八十万左右，这些患者失去肢体控制、自理等能力。而脑-机接口 （Brain-Computer Interface，BCI）提供了制造治疗和改善这些患者疾病的医疗 器械设备的技术基础，其主要机理是：训练患者进行相关的想象运动，同时采 集其脑电信号，对信号进行相应的特征提取与识别，然后将识别结果作为控制 量输入到相应的控制器以控制相应器械帮助患者实现日常的一些基本活动，从 而提高他们的生活质量。

由于大脑皮层对躯体运动具有交叉支配的特点，即一侧的皮质运动区支配 对侧的躯体骨骼肌运动。人脑在有运动的意愿后，无论是实际输出运动还是仅 仅想象运动（例如左手做抓、握运动或是右手做抓、握运动），都会导致对侧 大脑运动感觉区mu和beta节律的脑电信号失去同步的节律性活动而导致能量 减弱，而其同侧运动感觉区的mu和beta节律的脑电信号节律性活动加强。这 种现象称为事件相关去同步（Event-Related Desynchronization,ERD）和事 件相关同步（Event-Related Synchronization,ERS）的生理现象。这种现象 构成区分运动想象脑电信号的生理基础。

在研究大脑生理现象时，现有的仪器如Neuroscan等商用脑电采集仪具有 32导、64导和128导等多通道采集信号的能力。较多的通道数在一定程度上 可以提高脑电信号的空间分辨率，但存在成本高、不易便携等不足。而BCI在 实际应用中，易用性、便携性、成本低、抗干扰能力强等是实用化主要关注的 问题。本领域的技术人员考虑到如果仅采集几个通道，甚至是2-3个通道（除 参考极、参考地外）的脑电信号，采集设备的通道数就大大降低，从而使得设 备简单、经济，同时操作方便、鲁棒性也增强。

共空域模式（Common Spatial Pattern，CSP）方法在通道数比较多的情 况下是基于运动想象脑-机接口中非常有效的一种空间滤波方法，而在只有几 个通道EEG（脑电图描记器）信号的情况下，空间域信息量本身有限，直接应 用CSP的效果甚至不如直接用带通能量法（Band Power，BP）。现有技术中存 在一种想象单侧肢体运动的脑电特征的提取方法，该方法基于共空域模式CSP 方法的脑电特征提取方法。这种方法在少通道脑电信号的情况下就存在不足， 无法提供较高的分类正确率。

最近的一些以减少脑电通道数的研究，如Lal、Arvaneh等学者通过研究 通道优化选择问题，优化结果表明对于运动想象的脑-机接口，并非通道数目 越多越好，在他们各自提出的优化目标下不同的受试者有不同数目和位置的最 优通道布局。利用最优通道布局可以在用较少通道的脑电信号下有较高的分类 正确率。然而，针对少通道数目脑电信号情况下，如何在保持分类正确率最高 的情况下提取运动想象相关的特征这一关键问题并没有得到关注和解决。

因此，本领域的技术人员致力于开发一种脑-机接口（Brain-Computer  Interface，BCI）系统采集少量通道数目脑电信号（Electroencephalogram， EEG）情况下，运动想象脑电特征的提取方法。

发明内容

为了解决上述技术问题，本发明提供了一种少通道下运动想象脑电特征的 提取方法，该方法针对采集少通道脑电信号的情况，提出一种基于信息熵优化 选取多个坐标延迟的共空域频谱模式的特征提取方法，以提高基于运动想象的 脑电特征分类正确率。

本发明的基本思想是：由于大脑是公认的复杂的非线性动力学系统，将时 间坐标延迟利用低维时间序列重构多维相空间是重构复杂动力学系统的等价 状态空间的常用方法之一。最初，共空域频谱模型的提出，是基于每个通道时 间序列的一个坐标延迟提出的。该方法等价于同时优化空间滤波器与一个只有 单个可变参数的有限冲击响应（Finite Impulse Response，FIR）滤波器的原 理。由于一个参数的FIR滤波器的性能有限，本发明利用每个通道时间序列的 多个坐标延迟来同时优化空间滤波器与一个高阶多参数FIR滤波器，并用非参 数化方法估计特征量与类别间信息熵的方法来选取最优个数的时间坐标延迟。 基于多个时间坐标延迟的高阶FIR滤波器比一个参数的FIR滤波器有更好的幅 频特性，该方法提供了在少通道情况下提取想象运动类别特征的有效途径。

为实现上述目的，本发明提供的少通道下运动想象脑电特征的提取方法包 括以下步骤：

步骤（1）脑电信号采集：根据不同运动想象实验范式，采集运动想象区 的少数目通道脑电信号。

优选地，步骤（1）中的所述不同运动想象实验范式为左右手运动想象实 验范式。

优选地，步骤（1）中的少数目通道脑电信号为2～10路通道。

优选地，步骤（1）中所述采集的采样频率为250Hz。

步骤（2）数据预处理：对采集到的脑电信号用滤波器做带通滤波，在运 动想象实验范式中“提示”出现后的信号上做滑动窗记信号为矩阵X_t，时间间 隔为一个采样点。

优选地，步骤（2）中在运动想象实验范式中提示出现后的信号上做2秒 钟的滑动窗记信号为矩阵X_t。

优选地，步骤（2）中所述滤波器为四阶巴特沃斯滤波器。

优选地，步骤（2）中所述带通滤波为通频带8～30Hz带通滤波。

优选地，步骤（2）中所述带通滤波为通频带7～32Hz带通滤波

步骤（3）信号延时：δ^τ表示信号的延时操作，即δ^τ(X_t)=X_t-τ；将采集到 的每个通道的脑电信号做τ、2τ、…nτ的延时操作增加到原始输入信号中，这 样少通道的原始脑电信号扩展为多通道的较高维数的脑电信号。也即原信号延 时nτ后并入原信号成为一个新的高维信号：

$X_i=\left(\begin{array}{c}{X}_t\\ {\delta }^{\tau }\left(X_t\right)\\ {\delta }^{2\tau }\left(X_t\right)\\ ...\\ {\delta }^{\mathrm{n\tau }}\left(X_t\right)\end{array}\right).$

步骤（4）对新的高维信号进行CSP特征提取，分别由延时τ（也即n=1）、 延时2τ（也即n=2）至延时Nτ构成的N个高维信号提取到N个特征集{f₁(i)}_i、 {f₂(i)}_i…{f_N(i)}_i，其中f_n(i)∈R^2*p×1,i=1,...,L，L为训练集中训练样本的总数目， 2*p为利用CSP算法降维后特征向量的维数。

步骤（5）利用非参数化方法估计特征集{f_n(i)}_i与类别间的互信息熵

I(f_n;C)=H(C)-H(C|f_n)

其中I(f_n;C)是特征集{f_n(i)}_i与类别间的互信息熵，H(C)是类别标签集合C 的熵函数值，H(C|f_n)是在给定特征集f_n的情况下类别标签集合C的条件熵函 数值。

步骤（6）选取为最优特征集，n^*为最优个数延迟 因子。

步骤（7）对测试集脑电信号做τ、2τ...n^*τ的延时操作后，扩展成新的较 高维数的脑电信号，用CSP算法提取特征。

步骤（8）输出分类结果。

优选地，利用线性分类器LDA或是支持向量机做分类器输出分类结果。

在上述描述中，训练集是指运动想象实验中按照提示的类别（如左、右手） 进行运动想象过程中采集到的脑电信号集合。测试集是指运动想象实验中想象 类别未知的过程中采集到的脑电信号集合。

大脑是公认的复杂的非线性动力学系统，将时间坐标延迟利用低维时间序 列重构多维相空间是重构复杂动力学系统的等价状态空间的经典方法。最初的 共空域频谱模型是基于每个通道时间序列的一个坐标延迟提出的，其等价于同 时优化空间滤波器与一个只有单个可变参数的有限冲击响应（Finite Impulse  Response，FIR）滤波器。由于一个参数的FIR滤波器的性能有限，本发明利 用每个通道时间序列的多个坐标延迟来同时优化空间滤波器与一个高阶多参 数FIR滤波器，并用非参数化方法估计特征量与类别间信息熵的方法来选取最 优个数的时间坐标延迟。基于多个时间坐标延迟的高阶FIR滤波器比一个参数 的FIR滤波器有更好的幅频特性，该方法提供了在少通道情况下提取想象运动 类别特征的有效途径。

本发明提供的少通道下运动想象脑电特征的提取方法具有如下优点：本方 法利用原始信号的多个延时组成新的高维信号以克服原始信号空间域信息量 有限的不足，同时所求得的空频滤波器等价于空域滤波器中嵌入一个高阶FIR 滤波器，在优化空间域信息的同时优化频域信息；其次，最优延时因子个数由 特征集与类别间的互信息熵来优化，在最优延时因子确定后，特征提取只需要 将原始信号作相应个数延时，然后作线性降维求特征即可，除初始需要滤波， 如8-30Hz带通滤波外，不再需要其他滤波，计算量很小，在线实现方便。

附图说明

图1为本发明的较佳实施例中的少通道运动想象脑电特征提取方法的原理 示意图。

具体实施方式

以下结合附图对本发明较佳实施例的技术方案作进一步描述。

图1为本发明的较佳实施例中的少通道运动想象脑电特征提取方法的原理 示意图。在本实施例中，利用覆盖运动感觉区的少通道电极帽（如覆盖C3、 Cz、C4）采集脑电信号，根据左右手运动想象实验范式进行运动想象输出。

受试者按要求佩戴好脑电帽后坐在椅子上，显示器距受试者双眼约1米左 右。脑-机接口系统中脑电信号的采样频率设为250Hz，设置50Hz陷波器去除 工频干扰。

在实验流程中，每个训练样本开始时屏幕全黑，同时出现白色十字，2s钟 后‘哔’一声提示受试者集中精神准备开始。在第3s钟开始时有黑色条块叠加 在白色十字上提示运动想象方式。当黑色条块在左侧时代表想象左手运动（如 抓或握），当黑色条块在右侧时代表想象右手运动，提示持续1.25s。要求受试 者根据提示做相应的运动想象内容，持续想象3s以上，系统在第7s给出判断 结果。在8.5s钟时十字消失，一个训练样本的实验结束，受试者可以休息，休 息时间为1.5～2.5s，随机产生。左、右手运动想象各150次，顺序随机出现。 取左、右运动想象各75次做训练集，其他训练样本做测试集。

对采集到的脑电信号用四阶巴特沃斯滤波器做8-30Hz带通滤波，在运动 想象实验范式中提示出现后的信号上做2s钟的滑动窗，时间间隔为一个采样 点，为特征提取做准备。

将每个通道原始信号做τ、2τ、…nτ的延时操作增加到原始输入信号中， 将少通道的原始信号扩展为多通道的高维信号：

$X_i=\left(\begin{array}{c}{X}_t\\ {\delta }^{\tau }\left(X_t\right)\\ {\delta }^{2\tau }\left(X_t\right)\\ ···\\ {\delta }^{\mathrm{n\tau }}\left(X_t\right)\end{array}\right),$i=1,...,L，L为训练集中训练样本的总数目。

对新的高维信号进行CSP特征提取，具体按下列步骤完成：

第一步：将截取后的滑动窗内由原始信号扩展后的高维数字信号组成的数 字矩阵按提示信号的左XL_i、右XR_i类别进行分类，对每个动作类别的信号矩阵 求协方差矩阵，并对所有动作求平均。类别平均协方差矩阵由下式给出：

C_L=Exp(XL_i*XL_i^T)

C_R=Exp(XR_i*XR_i^T)

其中，Exp()为期望，下标L和R分别代表左、右两个类别，C_L和C_R为左、 右手运动想象的协方差矩阵。XL_i表示第i个训练样本是属于左手运动想象类别 的脑电信号，T是转置操作符。

第二步：将求得的两类个别平均协方差矩阵同时对角化，该目标由下述广 义特征值分解实现：

C_Lgw=λC_Rgw

其中w是所求的空域滤波器，λ为对应的广义特征值。

第三步：将广义特征值分解中对应最大的3个特征值以及最小的3个特征 值的特征向量构成CSP降维矩阵（也即空域滤波矩阵）

第四步：利用所求得的降维矩阵构造特征向量

$f\left(i\right)={\left(\begin{array}{ccc}\log \left(w_1^T\mathrm{cov}\left(X_i\right)w_1\right)& \log \left(w_2^T\mathrm{cov}\left(X_i\right)w_2\right)& L\end{array}\right)}^T$

式中w₁、w₂...为对应最大、第二大...等广义特征值的特征向量。

通过上述四个步骤完成CSP降维矩阵求解，以及求得每个训练样本的CSP 特征向量。这样的特征向量放在一起形成特征向量集{f_n(i)}_i。

改变延时操作nτ的数值n=1,2,...N，求得N组CSP特征向量集{f₁(i)}_i、 {f₂(i)}_i、...{f_N(i)}_i，i=1,...,L，其中L为训练集中训练样本的总数目。

利用非参数化方法估计特征量集{f_n(i)}_i与类别间的互信息熵，具体由下式 给出：

I(f_n;C)=H(C)-H(C|f_n)

式中P(ω)为类别的先验概率；给定特征变量f_n下 的条件熵为：

$H\left(C\right|f_n)=-\underset{\omega \in C}{\Sigma }P\left(\omega \right|f_n)\log \left(P\left(\omega \right|f_n)\right)=-\underset{\omega \in C}{\Sigma }\underset{i=1}{\overset{L}{\Sigma }}P\left(\omega \right|f_n\left(i\right))\log \left(P\left(\omega \right|f_n\left(i\right))\right)$

其中f_n(i)为延时因子个数为n的第i个训练样本的特征向量。P(ω|f_n(i))为 给定特征向量f_n(i)时类别ω的条件概率，根据贝叶斯定理可以求得：

$P\left(\omega \right|f_n\left(i\right))=\frac{p\left(f_n\left(i\right)\right|\omega )P\left(\omega \right)}{\underset{\omega }{\Sigma }p\left(f_n\left(i\right)\right|\omega )P\left(\omega \right)}.$

特征f_n(i)的条件概率密度函数p(f_n(i)|ω)用基于Parzen窗的非参数化方法 来估计，具体由下式给出。

$p\left(f_n\left(i\right)\right|\omega )=\frac{1}{n_{\omega }}\underset{y\left(j\right)\in I_{\omega }}{\Sigma }\phi (f_n\left(i\right)-f_n\left(j\right),h)$

其中y(j)∈I_ω表示第j个训练样本的特征向量的 类别标签属于ω的指标集。在高维数据中，Ψ常取对角阵形式，其对角元素为：

Ψ_kk=hσ²，k=1,...,2*p

σ为对应特征分量的方差，h为平滑参数取为用来控制窗函数 的宽度，2*p为特征向量的维数。于是，由上述式可以求得不同延时因子个数 下的I(f_n;C)。

接下来，选取为最优特征集，n^*为最优个数延迟因子。 对测试集中脑电信号做τ、2τ...n^*τ的延时操作后，得到测试集脑电信号的高维 扩展矩阵。然后再用训练集得到的CSP降维矩阵对测试集高维扩展矩阵提取特 征。

最后，利用线性分类器LDA或是支持向量机做分类器判别测试集中特征 向量所属类别，输出分类结果。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术 无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技 术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或 者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围 内。