一种基于多特征融合的多类运动想象脑电信号识别方法

本发明涉及运动想象脑电信号分类识别，具体涉及一种基于多特征融合的多类运动想象脑电信号识别方法。

背景技术：

目前，脑机接口（braincomputerinterface，bci）是一种不依赖于正常的由外围神经和肌肉组成的输出通路的通讯控制系统，其中最重要的是基于运动想象（mi）的bci系统，它可以将一个人对运动的想象转化为控制信号来控制外部设备。

脑电信号具有高度的受试者特异性和非平稳性，不同运动想象脑电信号识别的准确性仍然是亟待解决的问题，特征提取在运动想象脑电信号的分类中起到了至关重要的作用，因为信号具有维数高的特点，许多学者提出了多种特征提取方法，包括功率谱密度，基于小波的时频分析、希尔伯特－黄变换的时频能量和其他复杂特征来提取具有良好表示能力的特征，然而上述大多数研究都是基于单个领域进行特征提取，只能描述有限的信息，结果会使得脑电信号特征单一而且包含的信息描述不足，导致分类性能不理想，并且可能对不同的运动想象模式具有不同的敏感性，因此有必要对它们进行融合来增强对运动想象脑电信号的分类能力。

同时，为了充分利用不同领域的更多特征来提高分类性能，需要提取多域特征，然而多域特征通常具有较高维数，如果这些具有高维数的特征直接输入到分类器，就需要更复杂的模型来训练，甚至有可能降低分类性能，因此，开发一种不仅能够克服脑电信号提取特征单一、信息描述不足、识别精度低等缺陷，而且具有更好的鲁棒性和分类性能的脑电信号识别方法极具现实意义。

技术实现要素：

本发明的主要目的在于提供一种基于多特征融合的多类运动想象脑电信号识别方法，该方法将基于时－频－空的多域特征提取技术引入到多类运动想象脑电信号分类识别研究领域中，结合“信号预处理”、“特征提取及特征融合”、“模式识别”三个流程，实现多类运动想象脑电信号的识别。有效的提升了针对多类运动想象脑电信号的识别准确率，克服了传统单域特征提取算法信息描述不足的问题。

本发明采用的技术方案是：一种基于多特征融合的多类运动想象脑电信号识别方法，包括以下步骤：

第一步，进行运动想象脑电信号数据采集，对采集的信号进行数据预处理；

第二步，对预处理后的脑电信号采用ar模型法分析，提取脑电信号的时域特征；

第三步，使用变分模态分解－双谱分析法，对运动想象脑电信号进行频域特征提取；

第四步，使用共空间模式分析法提取运动想象脑电信号在空间域内的特征；

第五步，将二、三、四步提取的时域、频域、空域特征联合，构造一类运动想象的联合特征向量；

第六步，重复进行二到五步，提取各类运动想象所有样本信号的联合特征向量，构造联合特征矩阵；

第七步，对联合特征矩阵进行核主元分析，按照累计贡献率大于85%的原则选取主元个数，将样本数据向新的主元空间投影得到降维后的融合特征矩阵；

第八步，将融合特征作为训练数据输入到孪生支持向量机，并利用蝙蝠算法优化的孪生支持向量机进行参数寻优，完成分类器的训练；

第九步，针对测试数据，按二、三、四、五的处理过程构造联合特征向量，并向新的主元空间投影，得到融合特征向量，将融合特征向量输入到训练好的分类器实现运动想象脑电信号的识别。

进一步地，所述第一步中对运动想象脑电信号的数据预处理过程，是在类运动想象实验中，选取第类单次运动想象的个通道采集到的脑电信号，进行分析，的表达式如下：

利用共平均参考和butterworth带通滤波对采集到的多类运动想象脑电信号进行数据预处理，共平均参考计算方法是从所选通道中减去所有通道的平均值，计算公式为：

式中：为经过共平均参考处理后的单通道脑电信号；经过预处理后的脑电信号为：

。

更进一步地，所述第二步时域特征提取方法的具体步骤为：

步骤1：采用阶ar模型拟合，公式为：

其中为ar模型的系数，为均值为0，方差为的白噪声残差，提取模型系数作为该通道特征；

步骤2：对于选取的m个通道，求取每一通道的ar模型，并将模型系数作为该通道信号的特征，构造对应运动想象的时域特征向量，表达式如下：

。

更进一步地，所述第三步频域特征提取方法的具体步骤为：

步骤1：对进行变分模态分解，分解为个imf分量；变分模态分解是通过迭代搜索获取约束变分模型的最优解，自动获取固有模态函数的中心频率及带宽，的约束变分模型具体如下：

式中：和分别为imf分

量集和imf分量集对应的中心频率集；用交替方向乘数法迭代更新，，求解增广拉格朗日函数的鞍点，不断更新，，满足给定的判定表达式则迭代终止；

步骤2：对imf分量个数基于中心频率进行筛选，根据不发生过分解现象的原则，确定imf分量个数，具体为当分解为个分量时，第个分量的中心频率与第个分量的中心频率的差值，和第个分量的中心频率与第个分量的中心频率的差值在同一个数量级时，则分解为个分量出现了过分解，选择分解为个imf分量；

步骤3：对变分模态分解出的个imf分量进行双谱分析，分析过程如下：

设为变分模态分解出的imf分量，根据信号的离散傅里叶变换，双谱定义为：

式中为频率变量；为中每个imf分量的离散傅里叶变换；

；为数学期望，采用间接法进行双谱估计，首先估计所得的imf分量的三阶累积量，然后对该累积量进行二维的离散傅里叶变换，这样就得到了imf分量的双谱估计；

步骤4：计算双谱平面主域中所有的绝对对数之和作为频域特征，双谱特征的数学表达公式如下：

步骤5：利用变分模态分解－双谱分析法提取所选m个通道的频域特征，并且对提取到的频域特征进行信号平滑以及归一化处理，归一化处理公式如下所示：

其中是经过信号平滑及归一化后的双谱特征，是所选m个通道提取出来的双谱特征向量，表示的数据长度，min和max是计算最大值和最小值的运算符；构造的频域特征向量形式如下所示：

。

更进一步地，所述第四步空域特征提取方法的具体步骤为：

步骤1：求取的规范化协方差矩阵为：

其中表示矩阵的对角线元素的和；

步骤2：对每一类所有的样本数据的归一化协方差矩阵求平均值作为该类数据的平均归一化空间协方差矩阵，则任意两类样本数据的混合空间协方差矩阵为：

其中是规范化特征值的对角矩阵，并且对角线上的特征值按照递减的顺序排列，是每列对角矩阵元素对应的特征向量矩阵；构造白化矩阵并转化为协方差矩阵，首先构造白化矩阵和空间系数矩阵：

其中和有共同的特征向量；最后对和进行特征值分解：

其中是特征向量矩阵，它的每列是相应的特征向量，中的特征值按升序排列时，的特征值就是按降序排列的，两个特征值对角矩阵之和为单位矩阵，即，选取特征向量矩阵的前行和后行构成矩阵，其中；

步骤3：构造空间滤波器，由于两类任务的特征值和恒等于1，矩阵特

征值最大时，的特征值最小，所以：

对类问题进行分类时，对任意两类样本数据都计算一个投影矩阵

，对于，将通过构造的空间滤波器矩阵与相乘，就得到变换到共空间中的特征：

表示单个样本的模式特征矩阵，其中一行表示一个通道上的特征分布情况，提取每个通道的脑电信号的特征，再对特征值进行对数运算，空域特征向量的计算公式如下：

式中表示方差；将每类动作构建的空间滤波器得到的空域特征向量横向拼接，得到空域特征向量：

。

更进一步地，所述第五步具体步骤为：

步骤1：构造联合特征向量，将时、频、空域特征向量依次首位相连，

构造如下所示的联合特征向量：

联合特征向量为同一个样本数据，采用不同信号处理方法提取的各域内特征的联合，为一维特征向量。

更进一步地，所述第六步具体步骤为：

构造联合特征矩阵，将选取的个样本数据按步骤一的方法构造联合特征矩阵：

矩阵的每行代表一个样本的联合特征向量。

更进一步地，所述第七步核主元分析的具体步骤为：

步骤1：通过核函数将联合特征矩阵投影到线性可分的高维空间中进行核主元分析，获得主成分系数矩阵：

在主成分分析的同时可以得到样本协方差矩阵，它由特征值构成列向量，并且特征值按降序排列，可得特征值和相应的特征向量；

步骤3：比值反映了生成特征空间中第个主元的贡献率，根据特征值方差累计贡献率公式确定核主元个数，计算公式如下：

选择前个主元的累计贡献率超过85%，因此主成分系数矩阵变为：

将每个样本构造的联合特征向量向新的主元投影，就得到降维后的融合特征向量，投影方法如下：

。

更进一步地，所述第八步具体步骤为：

将所有样本降维后的特征向量作为训练集输入到孪生支持向量机中进行分类器的训练，并使用蝙蝠算法对孪生支持向量机的参数惩罚因子、惩罚因子、和核参数进行寻优，并以优化后的、、建立分类模型；设置相关参数，包括蝙蝠种群数、最大迭代次数、觅食空间维数；随机产生蝙蝠的位置和速度，蝙蝠的位置代表twsvm的三个参数惩罚因子、惩罚因子、和核参数，评估蝙蝠的个体适应度，并且找出当前时刻的最优解，适应度最大的蝙蝠个体即为全局最优解，本发明将惩罚因子和核参数交叉验证所得的识别准确率作为评价适应度的标准，定义为：

将蝙蝠种群位置代入评价适应度函数中，不断更新最优蝙蝠个体对应的位置，直到满足迭代终止条件，记录下最优的参数值，将全局最优参数带代入孪生支持向量机，完成分类器模型的训练。

更进一步地，所述第九步测试数据识别具体步骤为：

步骤1：将采集的数据分为训练与测试两类，用训练数据训练孪生支持向量机，测试数据用于测试分类效果；对于测试数据，按照第一、二、三、四步的方法对测试数据进行预处理和多域特征提取；

步骤2：构造出联合特征向量，然后向第五步得到的主成分系数矩阵投影，得到融合特征向量，具体公式如下：

步骤3：将得到的融合特征向量输入到训练好的分类器中进行多类运动想象脑电信号的识别。

本发明的优点：

本发明的方法将预处理后的脑电信号分别采用ar模型法、变分模态分解－双谱分析法以及共空间模式分析法提取脑电信号在时域、频域以及空间域内的特征，构造联合多特征向量，然后利用核主元分析法对联合多特征向量进行融合降维处理，通过提取贡献率大于85%的主元，生成融合特征，将融合特征作为训练数据输入到孪生支持向量机中训练分类器，并利用蝙蝠算法优化的孪生支持向量机进行参数寻优，然后通过参数优化后的孪生支持向量机实现测试数据的实时识别。本发明方法克服了传统单域特征提取算法信息描述不足的问题，有效地提升了针对多类运动想象脑电信号的识别准确率。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

图1是本发明的基于多特征融合的多类运动想象脑电信号识别方法流程示意图；

图2是本发明的多类运动想象脑电信号采集实验过程示意图；

图3是本发明实施例中选取a受试者动作1运动想象f3通道变分模态分解图。

图4是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第一双谱分析图；

图5是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第二双谱分析图；

图6是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第三双谱分析图；

图7是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第四双谱分析图；

图8是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第五双谱分析图；

图9是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第六双谱分析图；

图10是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第七双谱分析图；

图11是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第八双谱分析图；

图12是本发明实施例中对a受试者动作1运动想象f3通道imf分量的第九双谱分析图；

图13是本发明实施例中不同主元的第一三维散度图；

图14是本发明实施例中不同主元的第二三维散度图；

图15是本发明实施例中不同主元的第三三维散度图；

图16是本发明实施例中不同主元的第四三维散度图；

图17是本发明实施例中单域特征提取方法与融合特征识别率比较图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

参照图1，一种基于多特征融合的多类运动想象脑电信号识别方法，包括以下步骤：

第一步，实验采集三类运动想象eeg信号，使用emotivepoc+进行运动想象脑电信号数据采集，该设备包含14个电极通道（af3，f7，f3，fc5，t7，p7，o1，o2，p8，t8，fc6，f4，f8，af4）和2个参考电极（cms，drl），信号采样频率为128hz，实验采集肩部前屈、后伸以及外展三类运动想象eeg信号。本次实验研究共有七名被试者参与，所有被试者均身体健康，精神状态良好，实验时被试者在外界提示下进行运动想象，单次实验数据采集时间持续5秒，每个动作实验重复进行20次，实验流程示意图如图2所示。首先对运动想象脑电信号的数据进行预处理，选取fc5、f3、f4、fc6四个电极通道采集的运动想象过程中1～4s脑电信号进行分析，利用共平均参考和butterworth带通滤波（8~30hz）对采集到的多类运动想象脑电信号进行数据预处理。共平均参考计算方法是从所选通道中减去所有通道的平均值，计算公式为：

式中：为经过共平均参考处理后的单通道脑电信号。

经过预处理后的脑电信号为：

第二步，对预处理后的脑电信号采用ar模型法分析，提取脑电信号的时域特征，具体步骤如下：

步骤1：利用fpe准则对ar模型定阶，得到9阶为最佳的ar模型阶数，采用阶ar模型拟合，，公式为：

其中为ar模型的系数，为均值为0，方差为的白噪声残差。提取模型系数作为该通道特征。

步骤2：对于选取的4个通道，求取每一通道的ar模型，并将模型系数作为该通道信号的特征，构造对应运动想象的时域特征向量，为一36维特征向量，表达式如下：

第三步，使用变分模态分解－双谱分析法，对运动想象脑电信号进行频域特征提取，具体步骤如下：

步骤1：对进行变分模态分解，分解为3个imf分量。变分模态分解是通过迭代搜索获取约束变分模型的最优解，自动获取固有模态函数的中心频率及带宽，的约束变分模型具体如下：

式中：和分别为imf分量集和imf分量集对应的中心频率集。

用交替方向乘数法迭代更新，，求解增广拉格朗日函数的鞍点，不断更新，，满足给定的判定表达式则迭代终止。

步骤2：对imf分量个数基于中心频率进行筛选，根据不发生过分解现象的原则，确定imf分量个数3，具体为当分解为4个分量时，开始出现相邻模态的中心频率之差小于0。判定为出现过分解现象，因此分解层数选为3，如图3所示是a受试者动作1运动想象f3通道变分模态分解图。

步骤3：对变分模态分解出的3个imf分量进行双谱分析，分析过程如下：

设为变分模态分解出的imf分量，根据信号的离散傅里叶变换，双谱定义为：

式中为频率变量；为中每个imf分量的离散傅里叶变换；；为数学期望，采用间接法进行双谱估计，首先估计所得的imf分量的三阶累积量，然后对该累积量进行二维的离散傅里叶变换，这样就得到了imf分量的双谱估计，如图4至图12所示，是a受试者动作1运动想象f3通道imf分量的双谱分析图。

步骤4：计算双谱平面主域中所有的绝对对数之和作为频域特征，双谱特征的数学表达公式如下：

步骤5：利用变分模态分解－双谱分析法提取所选4个通道的频域特征，并且对提取到的频域特征进行信号平滑以及归一化处理，归一化处理公式如下所示：

其中是经过信号平滑及归一化后的双谱特征，是所选4个通道提取出来的双谱特征向量，表示的数据长度，min和max是计算最大值和最小值的运算符。

构造的频域特征向量为一12维的特征向量，形式如下所示：

第四步，使用共空间模式分析法提取运动想象脑电信号在空间域内的特征，具体步骤如下:

步骤1：求取的规范化协方差矩阵为：

其中表示矩阵的对角线元素的和。

步骤2：对每一类所有的样本数据的归一化协方差矩阵求平均值作为该类数据的平均归一化空间协方差矩阵，则任意两类样本数据的混合空间协方差矩阵为：

其中是规范化特征值的对角矩阵，并且对角线上的特征值按照递减的顺序排列。是每列对角矩阵元素对应的特征向量矩阵。

构造白化矩阵并转化为协方差矩阵，首先构造白化矩阵和空间系数矩阵：

其中和有共同的特征向量。

最后对和进行特征值分解：

其中是特征向量矩阵，它的每列是相应的特征向量，中的特征值按升序排列时，的特征值就是按降序排列的，两个特征值对角矩阵之和为单位矩阵，即，选取特征向量矩阵的前2行和后2行构成矩阵。

步骤3：构造空间滤波器，由于两类任务的特征值和恒等于1，矩阵特征值最大时，的特征值最小，所以：

对类问题进行分类时，对任意两类样本数据都计算一个投影矩阵，对于，将通过构造的空间滤波器矩阵与相乘，就得到变换到共空间中的特征：

表示单个样本的模式特征矩阵，其中一行表示一个通道上的特征分布情况，提取每个通道的脑电信号的特征，再对特征值进行对数运算，空域特征向量的计算公式如下：

式中表示方差，为一4维的向量。

将每类动作构建的空间滤波器得到的空域特征向量横向拼接，得到8维的空域特征向量：

第五步，构造联合多特征向量，然后利用核主元分析法对联合多特征向量进行融合降维处理，通过提取贡献率大于85%的主元，生成融合特征，具体步骤如下：

步骤1：构造联合特征向量，将时、频、空域特征向量依次首位相连，构造如下所示的联合特征向量：

联合特征向量为同一个样本数据，采用不同信号处理方法提取的各域内特征的联合，为一维特征向量。

第六步，重复进行二到五步，提取各类运动想象所有样本信号的联合特征向量，构造联合特征矩阵，具体步骤如下：

构造联合特征矩阵，将选取的个样本数据按步骤一的方法构造联合特征矩阵：

矩阵的每行代表一个样本的联合特征向量。

第七步，对联合特征矩阵进行核主元分析，按照累计贡献率大于85%的原则选取主元个数，将样本数据向新的主元空间投影得到降维后的融合特征矩阵。

步骤1：通过核函数将联合特征矩阵投影到线性可分的高维空间中进行核主元分析，获得主成分系数矩阵：

在主成分分析的同时可以得到样本协方差矩阵，它由特征值构成列向量，并且特征值按降序排列，可得特征值和相应的特征向量。

步骤3：比值反映了生成特征空间中第个主元的贡献率，根据特征值方差累计贡献率公式确定核主元个数，计算公式如下：

前个主元的累计贡献率超过85%，因此主成分系数矩阵变为：

将每个样本构造的联合特征向量向新的主元投影，就得到降维后的融合特征向量，投影方法如下：

如图13至图16所示是不同主元的分布散度图，由图可以看出，经过kpca降维后的融合特征具有较好的聚类能力。

第八步，将融合特征作为训练数据输入到孪生支持向量机，并利用蝙蝠算法优化的孪生支持向量机进行参数寻优，完成分类器的训练，具体步骤为：

将所有样本降维后的特征向量作为训练集输入到孪生支持向量机中进行分类器的训练，并使用蝙蝠算法对孪生支持向量机的参数惩罚因子、惩罚因子、和核参数进行寻优，并以优化后的、、建立分类模型。设置相关参数，包括蝙蝠种群数、最大迭代次数、觅食空间维数等。随机产生蝙蝠的位置和速度，蝙蝠的位置代表twsvm的三个参数惩罚因子、惩罚因子、和核参数，评估蝙蝠的个体适应度，并且找出当前时刻的最优解，适应度最大的蝙蝠个体即为全局最优解，本发明将惩罚因子和核参数交叉验证所得的识别准确率作为评价适应度的标准，定义为：

将蝙蝠种群位置代入评价适应度函数中，不断更新最优蝙蝠个体对应的位置，直到满足迭代终止条件，记录下最优的参数值，将全局最优参数带代入孪生支持向量机，完成分类器模型的训练。

第九步，针对测试数据，按二、三、四、五的处理过程构造联合特征向量，并向新的主元空间投影，得到融合特征向量，将融合特征向量输入到训练好的分类器实现运动想象脑电信号的识别，具体步骤如下：

步骤1：将采集的数据分为“训练”与“测试”两类，利用训练数据训练孪生支持向量机，“测试”数据用于测试分类效果。对于测试一个训练数据，按照第一、二、三、四步的方法对测试数据进行预处理和多域特征提取。

步骤2：构造出联合特征向量，然后向第五步得到的主成分系数矩阵投影，得到融合特征向量，具体公式如下：

步骤3：将得到的融合特征向量输入到训练好的分类器中进行多类运动想象脑电信号的识别。识别结果图如图17所示。可以观察到，多特征融合的特征提取方法相较于单域特征提取方法具有更高的识别率。

本发明的方法将预处理后的脑电信号分别采用ar模型法、变分模态分解－双谱分析法以及共空间模式分析法提取脑电信号在时域、频域以及空间域内的特征，构造联合多特征向量，然后利用核主元分析法对联合多特征向量进行融合降维处理，通过提取贡献率大于85%的主元，生成融合特征，将融合特征作为训练数据输入到孪生支持向量机中训练分类器，并利用蝙蝠算法优化的孪生支持向量机进行参数寻优，然后通过参数优化后的孪生支持向量机实现测试数据的实时识别。本发明方法克服了传统单域特征提取算法信息描述不足的问题，有效地提升了针对多类运动想象脑电信号的识别准确率。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。