一种多类运动想象脑电信号的识别方法与流程

本发明涉及脑电信号识别领域，特别涉及一种多类运动想象脑电信号的识别方法。

背景技术

随着计算机技术和脑科学的发展，人们开始尝试构建一种大脑与外界的交流通路，此通路不依赖外周神经和肌肉组织的参与，能将大脑信号解读成相应的命令来实现与外部世界的交流和控制，此通路被命名为脑-机接口(braincomputerinterface,bci)。bci这种人与外界交流通路，在康复工程与智能辅助机器人上的应用也引起了人们的重视，但是脑-机接口是一项多学科交叉技术，而且脑电信号具有非线性、非平稳性以及随机性强等特点，并且在采集信号的过程中极易受到各种噪声的干扰，那么如何有效、准确地对脑电信号进行特征提取和分类便是脑-机接口系统的关键所在。目前，eeg特征提取主流算法包括时域特征提取(自适应自回归模型)、频域特征提取(小波包分解)和空域下特征提取(公共空间模式)等，分类方法包括线性分类器、支持向量机和神经网络等。

脑电信号的识别过程一般是先对原始信号进行预处理，然后采用上述算法进行特征提取，再使用相关分类算法进行模式识别。bci系统利用此方法，取得了良好的分类效果，但对于运动想象脑电信号的识别，上述方法仍然存在多个不足，当使用自适应自回归模型、小波包分解或者公共空间模式进行运动想象脑电信号特征提取时，结果会使得脑电信号特征单一且所包含的信息描述不足；当使用线型分类器对多维脑电数据进行分类时，会使分类器适应度低，鲁棒性很差，极易受到噪声的干扰；当使用单核支持向量机进行分类时，会导致单核不能够很好的映射出脑电信号的分界面，造成识别精度低；当使用神经网络时，由于脑电信号维度的数量级很高，造成神经网络结构复杂，需要训练的参数较多，并且需要大量的训练样本，而现实中脑电信号的样本较少，结果得到的神经网络模型性能很差。

因此，开发一种不仅能够克服脑电信号特征单一、信息描述不足、对多维脑电数据适应度低及识别精度低等缺陷，而且具有更好鲁棒性和分类性能的脑电信号识别方法极具现实意义。

技术实现要素：

本发明的目的是为了解决上述方法存在的问题，提供一种多类运动想象脑电信号的识别方法。本发明有效地克服了传统单域特征提取算法信息描述不足以及单核分类器识别率较低的缺陷，为多类运动想象脑电信号识别提供了一种新思路。

为了达到上述目的，本发明采用如下的技术方案予以解决：

一种多类运动想象脑电信号的识别方法，将待识别的多类运动想象脑电信号经预处理和提取融合后得到的一维特征向量输入到多核学习支持向量机中，由其输出分类结果实现识别；

所述预处理是指去除噪声和眼电伪迹；所述提取融合是指分别利用离散小波变换和一对多公共空间模式提取时频域特征和空域特征后将各个域的特征首尾相连形成一维特征向量；

对待识别的多类运动想象脑电信号进行识别时需要的参数通过免疫遗传算法优化得到，所述参数包括使用一对多公共空间模式时保留的特征向量的个数m以及多核学习支持向量机基核函数的核参数γ和惩罚系数c，具体步骤如下：

(1)将类型已知的多类运动想象脑电信号进行如上所述的预处理后划分为训练集样本和测试集样本，同时以m、γ和c为抗体，初始化抗体种群chrom、记忆库bestindividuals、最优抗体bestfitness、当前代数gen＝0和停止代数maxgen；

(2)判断gen<maxgen是否成立，如果是，则进入下一步，反之，则输出最优抗体bestfitness即为优化得到的参数；

(3)对所有的样本进行如上所述的提取融合得到一维特征向量；

(4)通过simplemkl算法以训练集样本的一维特征向量和训练集样本的类型为输入和输出训练γ；

(5)将测试集样本的一维特征向量输入γ训练后的多核学习支持向量机中进行分类识别，获得识别率作为当前种群下抗体的适应度值fitness；

(6)计算抗体间的相互亲和力cv即抗体浓度concentration以及每个抗体的繁殖率p即繁殖系数excellence；

(7)将当代gen下最佳抗体bestchrom，判断|bestchrom-bestfitness|≥1e-7是否成立，如果是，则令gen＝0后进入下一步，反之，则令gen＝gen+1后进入下一步；

(8)更新最优抗体bestfitness、记忆库bestindividuals和父代；

(9)对父代进行选择、交叉、变异和重组，进而获得子代，之后转去执行步骤(2)。

作为优选的技术方案：

如上所述的多类运动想象脑电信号的识别方法，利用离散小波变换提取时频域特征的方法如下：首先，设多类运动想象脑电信号第p类x^p∈(x¹,x²,…,x^b)中的第e样本为信号的采样频率为fs，其中，b表示运动想象脑电信号的类别总数，n表示的电极导联数，xi表示第i电极导联的脑电信号，[x1,x2,…,xi,…,xn]t表示[x1,x2,…,xi,…,xn]矩阵转置；对进行l层离散小波变换，获得[8hz,30hz]的频率带的h个细节系数cdj∈(cd1,cd2,…,cdl)，j表示细节系数的所属层数，j∈(1,l)；利用离散小波变换分析原始eeg信号时，选择合适的分解层次以及小波函数尤为重要，分解层次是依据信号的采样率以及所需eeg信号的频率来确定。α波(8-13hz)和β波(14-30hz)被公认为是反映多类运动想象eeg信号特征的频段，所以只需提取8-30hz的信号特征。

然后，为进一步降低通过dwt提取的特征维数，根据h个细节系数在c3、c4及cz通道出现的事件相关同步(ers)和事件相关去同步(erd)现象，对上述三个通道的h个细节系数进行dwt特征提取；

最后，在此基础上进行统计学处理，令为编号为cl∈(c3、c4、cz)的通道中提取的cdj的细节系数，r表示实数域；分别以的均值、能量均值以及均方差作为的时频域特征，最终得到一个大小为9h的时频域特征向量ff；其中：

的均值

能量均值

均方差

其中表示的第k个值，n是指细节系数cdj的大小。

如上所述的多类运动想象脑电信号的识别方法，利用一对多公共空间模式提取空域特征的方法如下：

将多类运动想象脑电信号中的第1类所有样本x¹看作一类，剩余的其它类的所有样本(x²,……,x^b)则归为一类，计算两类的滤波器，依次分别将x²,……,x^b看作单独的一类，分别计算得到两类的滤波器，即获得b个滤波器投影矩阵，z^qn×m为第q个滤波器投影矩阵，q∈[1,b]，m为在使用一对多公共空间模式时保留的空域特征向量的个数；

再针对使用滤波器投影矩阵进行投影后，计算投影所得矩阵每行数据的方差值，并对所有行的方差值进行归一化处理，即得一个大小为4m的空域特征向量pf。

如上所述的多类运动想象脑电信号的识别方法，所述将各个域的特征首尾相连形成一维特征向量的方法如下：

使用4个径向基核分别对时频域特征向量ff、空域特征向量pf和时频空域特征向量[ff,pf](即特征ff和特征pf两者直接相连的组合)进行d1,g、d2,g和d3,g加权映射得到向量的1～9h维时频域特征、(9h+1)～(9h+4m)维空域特征和(9h+4m+1)～(18h+8m)维融合特征，然后将其首尾相连形成一维特征向量，其中g为径向基核的序号，g＝1,2,3,4。

如上所述的多类运动想象脑电信号的识别方法，所述训练集样本与测试集样本的数量比为4:1，训练集样本与测试集样本的划分随机进行。

如上所述的多类运动想象脑电信号的识别方法，一次迭代过程中所述通过simplemkl算法以训练集样本的一维特征向量和训练集样本的类型为输入和输出是指以训练集样本中所有样本的一维特征向量和类型为输入和输出；所述将测试集样本的一维特征向量输入γ训练后的多核学习支持向量机是指将测试集样本中所有样本的一维特征向量输入γ训练后的多核学习支持向量机。

如上所述的多类运动想象脑电信号的识别方法，所述计算抗体间的相互亲和力cv即抗体浓度concentration以及每个抗体的繁殖率p即繁殖系数excellence的具体公式如下：

其中nc为抗体总数，lc为抗体长度，kv,s表示抗体v与抗体s中相同的位数，threshold为抗体浓度阈值，ps为种群多样性参数，fitnessv为抗体v的适应度值。

如上所述的多类运动想象脑电信号的识别方法，所述最佳抗体bestchrom是指当代gen下适应度值fitness最高的抗体。

如上所述的多类运动想象脑电信号的识别方法，所述交叉为低位交叉，所述变异为高位变异。

如上所述的多类运动想象脑电信号的识别方法，所述m、γ及c采用151位0/1进行表示，1～5位用于表示m，6～19位表示c，20～151位表示γ；m∈(1,30)，精度为1；c∈(2^-10,2¹⁰)，精度为0.01；γ∈(2^-10,2¹⁰)，精度为0.1。

有益效果：

(1)本发明的一种多类运动想象脑电信号的识别方法，将多域的特征进行加权融合，综合考虑脑电信号在时频空域下的特征，克服了传统单域特征提取算法信息描述不足的特点；

(2)本发明的一种多类运动想象脑电信号的识别方法，该方法使用多个核对应融合特征，而且多个不同核对应不同域的特征，使得分类器具有更好的鲁棒性，能够自适应的改变每个域特征所占的比重，获得最适合受试者的多域特征组合，使得算法具有很好的分类性能和较高的识别精度；

(3)本发明的一种多类运动想象脑电信号的识别方法，应用广泛，主要应用于大脑工作记忆、脑电信号特征提取、脑-机接口系统和康复工程与智能辅助机器人等领域。

附图说明

图1为待识别的多类运动想象脑电信号进行识别时所需参数的免疫优化算法流程图；

图2为不同小波函数下的单核svm分类器分类正确率图；

图3为不同基核个数的训练耗时对比图；

图4为采集的840组样本分别在四种算法下的进化代数和识别率关系图；

图5为在本发明算法下每位受试者进化代数和识别率关系图；

图6为在csp+dwt-svm算法下每位受试者进化代数和识别率关系图；

图7为在csp-svm算法下每位受试者进化代数和识别率关系图；

图8为在dwt-svm算法下每位受试者进化代数和识别率关系图；

图9是在四种算法下k3b受试者进化代数和识别率关系图。

具体实施方式

下面结合具体实施方式，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

一种多类运动想象脑电信号的识别方法，具体步骤如下：

(一)对待识别的多类运动想象脑电信号进行识别时需要的参数通过免疫遗传算法进行优化，流程图如图1所示，具体为：

(1)将类型已知的多类运动想象脑电信号预处理进行去除噪声和眼电伪迹后划分为训练集样本和测试集样本，同时以m、γ和c为抗体，初始化抗体种群chrom、记忆库bestindividuals、最优抗体bestfitness、当前代数gen＝0和停止代数maxgen，其中训练集样本与测试集样本的数量比为4:1，训练集样本与测试集样本的划分随机进行，m表示使用一对多公共空间模式时保留的特征向量的个数，γ表示多核学习支持向量机基核函数的核参数，c表示惩罚系数，m、γ及c采用151位0/1进行表示，1～5位用于表示m，6～19位表示c，20～151位表示γ；m∈(1,30)，精度为1；c∈(2^-10,2¹⁰)，精度为0.01；γ∈(2^-10,2¹⁰)，精度为0.1；

(2)判断gen<maxgen是否成立，如果是，则进入步骤(3)，反之，则输出最优抗体bestfitness即为优化得到的参数；

(3)对所有的样本分别利用离散小波变换和一对多公共空间模式提取时频域特征和空域特征后将各个域的特征首尾相连形成一维特征向量提取融合得到一维特征向量；具体如下：

(3.1)利用离散小波变换提取时频域特征的方法如下：(3.1.1)设多类运动想象脑电信号第p类x^p∈(x¹,x²,…,x^b)中的第e样本为信号的采样频率为fs，其中，b表示运动想象脑电信号的类别总数，n表示的电极导联数，xi表示第i电极导联的脑电信号，[x1,x2,…,xi,…,xn]^t表示[x1,x2,…,xi,…,xn]矩阵转置；(3.1.2)对进行l层离散小波变换，获得[8hz,30hz]的频率带的h个细节系数cdj∈(cd1,cd2,…,cdl)，j表示细节系数的所属层数，j∈(1,l)；(3.1.3)根据h个细节系数在c3、c4及cz通道出现的事件相关同步和事件相关去同步现象，对上述三个通道的h个细节系数进行dwt特征提取；(3.1.4)令为编号为cl∈(c3、c4、cz)的通道中提取的cdj的细节系数，r表示实数域；分别以的均值、能量均值以及均方差作为的时频域特征，最终得到一个大小为9h的时频域特征向量ff；其中：

的均值

能量均值

均方差

其中，表示的第k个值，n是指细节系数cdj的大小；

(3.2)利用一对多公共空间模式提取空域特征的方法如下：(a)将多类运动想象脑电信号中的第1类所有样本x¹看作一类，剩余的其它类的所有样本(x²,……,x^b)则归为一类，计算两类的滤波器，依次分别将x²,……,x^b看作单独的一类，分别计算得到两类的滤波器，即获得b个滤波器投影矩阵，z^qn×m(q∈[1,b])为第q个滤波器投影矩阵，m为在使用一对多公共空间模式时保留的空域特征向量的个数；(b)针对使用滤波器投影矩阵进行投影后，计算投影所得矩阵每行数据的方差值，并对所有行的方差值进行归一化处理，即得一个大小为4m的空域特征向量pf；

(3.3)将各个域的特征首尾相连形成一维特征向量的方法如下：使用4个径向基核分别对时频域特征向量ff、空域特征向量pf和时频空域特征向量[ff,pf]进行d1,g(g＝1,2,3,4)、d2,g(g＝1,2,3,4)和d3,g(g＝1,2,3,4)加权映射得到向量的1～9h维时频域特征、(9h+1)～(9h+4m)维空域特征和(9h+4m+1)～(18h+8m)维融合特征，然后将其首尾相连形成一维特征向量，其中g为径向基核的序号；

(4)通过simplemkl算法以训练集样本中所有样本的一维特征向量和训练集样本的类型为输入和输出训练γ；

每种特征对应的径向基核数量较少时，simplemkl算法的收敛速度会较慢，每种特征对应的径向基核数量与simplemkl算法训练收敛速度如图3所示，当选择每种特征对应的基核数量为4时达到最佳。故可确定rbf的个数为12个，每种特征对应4个rbf径向基核，而每个rbf具有一个核参数γ，且γ的取值直接影响到最终的分类性能，因此也会有12个γ需要优化。

(5)将测试集样本中所有样本的一维特征向量输入γ训练后的多核学习支持向量机中进行分类识别，获得识别率作为当前种群下抗体的适应度值fitness；

(6)计算抗体间的相互亲和力cv即抗体浓度concentration以及每个抗体的繁殖率p即繁殖系数excellence，具体公式如下：

其中nc为抗体总数，lc为抗体长度，kv,s表示抗体v与抗体s中相同的位数，threshold为抗体浓度阈值，ps为种群多样性参数，fitnessv为抗体v的适应度值；

(7)将当代gen下适应度值fitness最高的最佳抗体bestchrom，判断|bestchrom-bestfitness|≥1e-7是否成立，如果是，则令gen＝0后进入下一步，反之，则令gen＝gen+1后进入下一步；

(8)更新最优抗体bestfitness、记忆库bestindividuals和父代；

(9)对父代进行选择、低位交叉、高位变异和重组，进而获得子代，之后转去执行步骤(2)；

(二)将待识别的多类运动想象脑电信号经上述预处理和提取融合后得到的一维特征向量输入到参数优化后的多核学习支持向量机中，由其输出分类结果实现识别。

本发明利用分类正确率作为分类结果的评定标准，为验证本发明的有效性，选取的运动想象eeg数据来源于第三届国际bgi竞赛数据集datasetⅲ，该数据记录了三位受试者的基于视觉提示的四类(分别是左手、右手、脚和舌头)运动想象思维活动。采样频率是250hz，信号来自60个导联通道，保存数据之前对数据进行了1-50hz的带通滤波和50hz的陷波滤波。其中受试者k3b四种想象任务分别进行90次，共进行了360次实验，其他两名受试者k6b与l1b四种想象任务分别进行60次，各进行了240次实验，三位受试者共提供840组样本。本实施例中各参数具体设置如表1所示：

表1

本实施例使用的eeg信号数据的采样频率为250hz，则dwt分解后不同层次所对应的频率信息如表2所示：

表2

可知，α波和β波主要集中在d3、d4子频带中，因此选择进行4层小波分解，每个eeg样本所具有的时频特征维度是18。小波基函数常用的有四种：daubechies(db4)、coiflet(coif3)、symlet(sym2)和biorthog-onal(bior3.1)。分别用这4种小波函数进行18维时频特征提取，之后使用单核svm做分类，其结果如图2所示，结果表明使用db4小波基函数的识别率优于其它三种，因而本发明最终选择db4作为小波基函数。

首先对三人的840组数据集(左手、右手、脚和舌头各210次)采用hold-out交叉验证法来计算识别率，即每次实验随机将每人每类数据的80％作为训练集，剩余20％作为测试集，进行分类识别。其次，再测试模型对每个人的识别率。模型参数寻优时，取每组参数对应3次识别率的均值作为最终识别率进行评判。

在实验过程中，将本发明方法与加入免疫遗传参数优化的csp-svm、dwt-svm、csp+dwt-svm三种算法进行脑电信号分类结果对比，其中对比csp-svm、dwt-svm算法是用于说明单域下特征的信息描述与多域下特征的信息描述对分类结果的影响，而对比csp+dwt-svm算法是为了反映单核svm分类器和多核svm分类器对解决多域特征信息融合分类问题的差异。实验首先获得模型对k3b、k6b、l1b三位受试者共840组样本的识别率，四种算法的进化过程中进化代数与识别率关系如图4所示；其次获得模型对每个受试者的识别率，对每位受试者的四种算法的进化过程如图5～8所示。多次实验的最优识别率如表3所示。对三位受试者共840组样本，本文算法的识别率为94.21％，利用csp-svm算法获得的分类正确率为83.21％，dwt-svm为77.11％，csp+dwt-svm为82.69％，本发明方法相比上述三种方法分别提高11.0％、17.1％、11.52％。由于每位受试者提供的样本数不同以及信号质量存在差异，使得三位受试者在相同算法下出现识别率不相同，但是本文所提出的算法都在92％以上，其中k3b受试者达到98％，而其它三种算法都未达到90％。

表3

从表3四种算法识别率对比结果可以看出，不同受试者对不同域下的特征适应度不同，即有些受试者的eeg信号更适合在空域特征下分类，例如受试者k6b在空域特征下的识别率为82.55％，相比其它域特征提高了8％；有些受试者可能在时频空域特征下更容易区分，例如l1b在时频空域融合特征下取得了83.13％的识别率，相比其它域提高5％。由于受试者k3b的样本比另外两位受试者的样本多120，因此在各域特征下的识别率相对较高，本发明又针对k3b受试者样本在不同域下做了多次实验，试验中的样本数与其他两位受试者的样本数相同，并且每次实验是从360组样本中随机选择出240组样本。其中一组实验的进化图如图9，多次实验的结果如表4：

表4

表中，1表示k3b受试者，2表示k6b受试者，3表示l1b受试者。从表4实验结果可以看出k3b在不同域下的识别率虽然也有差异，但是差别较小，说明k3b受试者的eeg信号对不同域下的特征适应度相似，也进一步说明不同受试者获得最优识别率对应不同的域特征，但是3位受试者在时频域、空域、时频空域下的分类结果还是较低，反映出单域特征或者直接将多域特征简单融合都具有信号特征单一、信息描述不足的特点。

综上各实验结果表明，本发明算法识别率有较大提升，原因在于，一方面本发明算法将多个域特征进行融合，使得描述eeg信号的特征不再单一，获得较充足的描述信息；另一方面本发明算法不是将多个域的特征直接简单融合，由于不同受试者对不同域下的特征适应度不同，如果只是简单融合多域特征可能会降低识别率，例如受试者k6b，为了使算法具有更好的鲁棒性，即无论哪个域的特征更利于受试者进行eeg信号分类，算法都能获得较高的识别率，本发明通过对每个域特征对应的4个rbf核加权映射，在训练样本对权系数d1,g、d2,g或者d3,g自适应训练后，会获得最适合受试者的多域特征组合，使得算法具有很好的鲁棒性和较高的识别率。