一种基于多通道频率特征定制的左右手运动想象脑电信号分类方法与流程

本发明属于运动想象脑电信号识别领域，尤其涉及一种基于多通道频率特征定制的左右手运动想象脑电信号分类系统。

背景技术：

脑卒中，俗称“中风”，是最典型的脑血管病之一。调查显示，世界受中风影响的绝对人数正在上升。而我国脑卒中的患病率一直高于全球，是国民第一死亡原因。这些患者通过正常的大脑输出途径来完成与外界环境的交流十分困难。因此，为了提高脑卒中患者的生活质量，减轻疾病对家庭与社会的负担，我们期望可以通过人工智能技术辅助并提高脑卒中等手部偏瘫患者的康复治疗。

脑机接口(brain-computer-interface,bci)是目前国际上广泛研究的一种脑卒中康复方式，使用bci技术的机器人康复设备，利用脑信号采集装置和智能算法自动检测人脑主动活动的电信号，并通过外接机电设备对人体相应的部位产生刺激或引导活动，人工建立起神经-肌肉的运动通路，刺激原神经回路的重建。这种方式通过提高对剩余脑回路的利用辅助运动功能的恢复与改善。采集头皮表面脑电信号的eeg(electroencephalography)技术是一种比较方便、安全和低廉的bci技术。基于运动想象(motorimagery,mi)的eeg-bci系统由于其依赖于使用者想象肢体的运动动作，涉及了实际运动的控制机制和神经基质，可以促进大脑运动感觉皮层(sensorimotorcortex,sm1)的激活，尤其适用于脑卒中患者的手部康复治疗。实现一个运动想象分类系统包括信号采集设计、脑电信号数据预处理、脑电信号特征提取以及分类算法设计等关键技术。但由于eeg信号的非稳态特性且存在个性化差异的问题，目前大多数方法忽略了用户数据的多样性，往往在实验中达到的高效率无法推广到实际应用中。

技术实现要素：

针对现有技术存在的以上问题，本发明旨在设计一种基于多通道频率特征定制的左右手运动想象脑电信号分类方法，通过将相关通道去噪信号的功率谱密度作为特征在时频二维平面上分区，并分别进行svm分类，根据分类的交叉验证得分来评估特征区的分类贡献度，并选择贡献度大于阈值0.75的通道频率作为定制频率，相应的通道频率特征作为定制特征用于svm分类器的训练和在线识别。由此改善脑电信号的个性化差异对运动想象分类在实际应用中的影响。

为实现本发明的目的，本发明所采用的技术方案如下：一种基于多通道频率特征定制的左右手运动想象脑电信号分类系统，具体包括：通过fir滤波器和独立成分分析去除采集信号的噪声，对其进行epochs数据段截取后，提取10-20电极分布系统中fc(fc1、fc2和cz)、c(c3、c4和cz)和cp(cp1、cp2和cz)三组电极通道的信号数据，并使用multitaper与morlet小波结合的方法计算其在分析频率段foi的功率谱密度，作为信号特征在foi-epochs时频平面上进行分区；分别做svm分类得到每个特征区的分类得分用来衡量其分类贡献度，标记贡献度大于0.75的通道频率作为定制通道频率，并将定制通道频率特征输入svm分类器进行调参和模型训练。在线识别直接提取定制通道频率特征输入训练好的svm分类器，得到识别结果。

为了实现上述目的，本发明的技术方案如下：一种基于多通道频率特征定制的左右手运动想象脑电信号分类方法，该方法包括以下几个步骤：

步骤1：采集用户的脑电信号；

步骤1.1：开启脑电数据采集设备，将电脑usb接口连接数据采集设备的digital口，打开电脑实验范式程序(即实验引导gui)；程序运行过程中，电脑通过usb接口发送两种串口信号给数据采集设备，标记两种想象状态的触发时间点；

步骤1.2：用户正确佩戴上符合10-20电极分布系统的脑电帽，并在电极腔中注射导电膏，使得电极导通(电阻≤5kω)；在用户手腕、眼眶上下分别贴上电极片，作为参考电极记录伪迹信号；

步骤1.3：采集实验中，用户须端正地坐在扶手靠椅上，身体放松，双手自然搭在扶手上方，掌心向上；实验过程中须注视距离其头部约40-60cm远的显示屏，按照屏幕指令操作；当用户准备好时，即可开始实验；一轮(run)实验包含40次想象尝试(trials)，每次尝试包含1个手部运动想象任务，左/右手握拳运动想象任务各20个，出现次序随机；执行运动想象任务之前，屏幕上呈现3s“十”字提示用户集中注意力；之后屏幕上显示1s“准备”并伴有“滴”提示音；提示音停止随即开始4s的运动想象任务，此时屏幕的提示为“十”字上添加向左或向右的箭头，表示受试者进行左/右手运动想象任务；运动想象任务结束之后是1.5s的休息时间；在进行下一次尝试之前，还有小于1s的随机时间间隙；离线训练需要每位用户完成3轮数据采集，共120个trials；

步骤2：对采集数据进行预处理，即fir滤波和独立成分分析(ica)；

步骤2.1：使用窗函数为hamming窗的有限冲击响应(fir)滤波器分别对采集到的原始数据进行1hz高通滤波和45hz的低通滤波，其中低通过渡带为min(max(l_freq*0.25,2),l_freq)，高通过渡带为这里的l_freq、h_freq和sfreq分别表示下通带频率、上通带频率和采样频率；滤波器长度设定为最短过渡带倒数的6.6倍；这里可以去除电源噪声(50hz或60hz)和高频率的肌电伪迹(约60hz)；

步骤2.2：滤波后，使用独立成分分析(ica)的fastica算法对32个脑电(eeg)通道的信号提取10个独立成分，将10个独立成分与眼部伪迹电极信号中的眼动成分进行相关性分析，设定阈值2.5，自适应z得分高于阈值的独立成分被认为眼动伪迹并剔除，将剩下的独立成分逆映射回电极空间；至此，预处理完成，将得到相对洁净的数据；

步骤3：对步骤2中得到的去噪数据根据同步信号的标记点分割出运动想象开始前2s至想象任务结束后2s时间内的数据段epochs，并提取出位于大脑运动感觉区的fc、c和cp三组电极通道的epoch数据，每组数据都包括10-20电极分布系统中的左右两个电极(fc1和fc2、c3和c4、cp1和cp2)和中心电极cz；对三组数据计算psd，这里使用的multitaper算法其实是morlet小波与倍数dpss(discreteprolatespheroidal)multitaper窗相结合的；令采样周期为δ，则算法的具体步骤如下：

(1)选择f＝[5hz,37hz]作为分析频率段(frequenciesofinterest,foi)，每个频率段morlet小波的长度ncycles设定为7.0，就是说对于频率f∈f，小波窗(和taper向量)的长度为

(2)对每个频率f∈f求dpss序列：对φr进行主成分分析(pca)，φr为nf×nf的矩阵，φr的每个元素

其中，r是频率分辨率，且α是时间带宽积，这里取α＝2；因为taper的数量l需要满足

所以l取3，则得到的3个主成分即为3个相互正交的tapers向量，它们组成的序列就是dpss，或称slepian序列；

(3)将每个频率f∈f且长度为nf的morlet小波函数分别与dpss序列中每个taper向量做点乘，得到关于每个频率f的morlet-dpss序列，序列维度为(3×33×nf)；

(4)对每个通道的每段epoch数据的快速傅里叶变换(fft)与morlet-dpss序列中每个taper向量的快速傅里叶变换相乘并做快速傅里叶逆变换(ifft)，得到每个频率f的multitaper系数(c表示通道序号，k表示epoch序号)；

(5)与其共轭相乘得到频率f的功率谱密度

至此，便得到三组数据在foi-epoch时间频率二维平面上的功率谱密度psd分布；

步骤4：对步骤3中得到的三张foi-epoch时频二维特征平面，将每组通道数据的特征分为面积为1hz×0.125s的2112个bins，共6336个bins；每一个bin的特征分别输入svm分类器中进行训练，并使用5×5foldcross-validation验证每个bin特征的分类auc得分；这里，svm核函数为径向基函数rbf，n_features表示特征数量，当采样频率sfreq＝1000hz时，n_features＝125，惩罚因子c设定为15；通过上述过程，可以得到三组通道数据在foi-epoch特征平面的svm分类得分分布featurescores，其维度为(3×33×64)；

步骤5：对步骤4中得到的三组通道数据foi-epoch特征平面的svm分类得分featurescores设定阈值thresholdscores＝0.75，筛选出每组通道数据大于阈值的通道频率特征，组成定制通道频率特征向量z0,并标记定制的通道特征频率fo；可以注意到，由于在线识别无法对时间特征进行选择，所以这里仅仅需要提取符合要求的频率特征；令定制通道频率特征z为输入，以rbf为核函数的svm的参数(c,gamma)作为调参对象，选择在5×5foldcross-validation中auc得分最高的一组参数(co,gammao)，并将此时训练的svm模型作为在线识别的分类器；

步骤6：离线训练模型用于在线识别，即实时采集用户脑电数据，依次进行预处理、提取定制通道频率fo的数据并计算通道频率特征z、输入svm模型进行识别，最终得到识别结果反馈给用户；

步骤6.1：用户佩戴好采集装置后，跟随显示屏指示进行左右手握拳的运动想象；电脑实时接收采集设备输送的脑电信号，并在接收到触发信号后，对每个采样时刻的数据进行fir高通(1hz)和低通(45hz)滤波后，进行独立成分分析剔除伪迹噪声；

步骤6.2：对去噪信号进行定制特征提取和分类识别操作：

(1)根据定制通道特征频率fo，对特定电极通道的信号做带通滤波，中心频率为f∈fo，带宽为1hz，提取全部通道特征频率的数据；

(2)对(1)中提取的数据使用multitaper算法计算psd特征，即定制通道频率特征z；

(3)将z输入svm分类器识别想象内容；

(4)将识别结果输出给反馈装置，形成闭环。

与现有技术相比，本发明的技术方案具有以下的有益技术效果：

1、本发明提出了一种在运动想象最相关的三组通道(10-20电极分布系统中的fc(fc1、fc2和cz)、c(c3、c4和cz)和cp(cp1、cp2和cz))上做psd特征提取，能够有效消除信号的冗余，降低数据计算量，提高在线识别的速度；

2、针对运动想象脑电信号的个性化差异问题，对三组通道的信号特征在时频平面上作分区，对特征区的分类贡献打分，筛选得分大于0.75的特征频率作为定制通道频率，仅计算该频率范围内的特征，使得特征的维数降低，但质量提高，在提高在线识别的速度同时，改善在线识别的效果；

3、不同用户可以得到不同的定制通道频率，但由于阈值的限定，定制通道频率特征的质量不会受到用户差异的影响。因此，在兼顾脑电信号个性化差异的同时，也能保证系统分类的效果。

4、定制特征从通道和频率两方面出发，与离线训练的数据无关，避免了很多现有技术处理脑电信号个性化问题时由于过分关注训练数据而导致的过拟合问题，本方法具有较好且稳定的泛化能力。

附图说明

图1是本发明一种基于多通道频率特征定制的左右手运动想象脑电信号分类系统的框架示意图；

图2是本发明采用通道频率特征定制的方法流程图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

实施例1：如图1所示，本发明提出了一种基于多通道频率特征定制的左右手运动想象脑电信号分类方法，该方法的详细步骤为：

步骤1：采集用户的脑电信号；

步骤1.1：开启脑电数据采集设备，将电脑usb接口连接数据采集设备的digital口，打开电脑实验范式程序(即实验引导gui)。程序运行过程中，电脑通过usb接口发送两种串口信号给数据采集设备，标记两种想象状态的触发时间点。

步骤1.2：用户正确佩戴上符合10-20电极分布系统的脑电帽，并在电极腔中注射导电膏，使得电极导通(电阻≤5kω)。在用户手腕、眼眶上下分别贴上电极片，作为参考电极记录伪迹信号。

步骤1.3：采集实验中，用户须端正地坐在扶手靠椅上，身体放松，双手自然搭在扶手上方，掌心向上。实验过程中须注视距离其头部约40-60cm远的显示屏，按照屏幕指令操作。当用户准备好时，即可开始实验。一轮(run)实验包含40次想象尝试(trials)，每次尝试包含1个手部运动想象任务，左/右手握拳运动想象任务各20个，出现次序随机；执行运动想象任务之前，屏幕上呈现3s“十”字提示用户集中注意力；之后屏幕上显示1s“准备”并伴有“滴”提示音；提示音停止随即开始4s的运动想象任务，此时屏幕的提示为“十”字上添加向左或向右的箭头，表示受试者进行左/右手运动想象任务；运动想象任务结束之后是1.5s的休息时间；在进行下一次尝试之前，还有小于1s的随机时间间隙。离线训练需要每位用户完成3轮数据采集，共120个trials。

步骤2：对采集数据进行预处理，即fir滤波和独立成分分析(ica)；

步骤2.1：使用窗函数为hamming窗的有限冲击响应(fir)滤波器分别对采集到的原始数据进行1hz高通滤波和45hz的低通滤波，其中低通过渡带为min(max(l_freq*0.25，2)，l_freq)，高通过渡带为这里的l_freq、h_freq和sfreq分别表示下通带频率、上通带频率和采样频率；滤波器长度设定为最短过渡带倒数的6.6倍。这里可以去除电源噪声(50hz或60hz)和高频率的肌电伪迹(约60hz)。

步骤2.2：滤波后，使用独立成分分析(ica)的fastica算法对32个脑电(eeg)通道的信号提取10个独立成分，将10个独立成分与眼部伪迹电极信号中的眼动成分进行相关性分析，设定阈值2.5，自适应z得分高于阈值的独立成分被认为眼动伪迹并剔除，将剩下的独立成分逆映射回电极空间。至此，预处理完成，将得到相对洁净的数据。

提取10-20电极分布系统中的三组通道数据，即电极fc1、fc2和cz、电极c3、c4和cz，以及电极cp1、cp2和cz，使用multitaper方法计算信号功率谱密度(psd)；

步骤3：对步骤(2)中得到的去噪数据根据同步信号的标记点分割出运动想象开始前2s至想象任务结束后2s时间内的数据段epochs，并提取出位于大脑运动感觉区的fc、c和cp三组通道的epoch数据，每组数据都包括左右两个电极(fc1和fc2、c3和c4、cp1和cp2)和中心电极cz。对三组数据计算psd，这里使用的multitaper算法其实是morlet小波与倍数dpss(discreteprolatespheroidal)multitaper窗相结合的。令采样周期为δ，则算法的具体步骤如下：

1)选择f＝[5hz,37hz]作为分析频率段(frequenciesofinterest,foi)，每个频率段morlet小波的长度ncycles设定为7.0，就是说对于频率f∈f，小波窗(和taper向量)的长度为

2)对每个频率f∈f求dpss序列：对φr进行主成分分析(pca)，φr为nf×nf的矩阵，φr的每个元素

其中，r是频率分辨率，且α是时间带宽积，这里取α＝2。因为taper的数量l需要满足

所以l取3，则得到的3个主成分即为3个相互正交的tapers向量，它们组成的序列就是dpss，或称slepian序列；

3)将每个频率f∈f且长度为nf的morlet小波函数分别与dpss序列中每个taper向量做点乘，得到关于每个频率f的morlet-dpss序列，序列维度为(3×33×nf)；

4)对每个通道的每段epoch数据的快速傅里叶变换(fft)与morlet-dpss序列中每个taper向量的快速傅里叶变换相乘并做快速傅里叶逆变换(ifft)，得到每个频率f的multitaper系数(c表示通道序号，k表示epoch序号)；

5)与其共轭相乘得到频率f的功率谱密度

至此，便得到三组数据在foi-epoch时间频率二维平面上的功率谱密度psd分布。

将三组通道信号的psd在时-频二维平面上进行分区，并逐一输入svm分类器,得到psd特征在时-频平面上的分类得分分布；

步骤4：对步骤(3)中得到的三张foi-epoch时频二维特征平面，将每组通道数据的特征分为面积为1hz×0.125s的2112个bins，共6336个bins。每一个bin的特征分别输入svm分类器中进行训练，并使用5×5foldcross-validation验证每个bin特征的分类auc得分。这里，svm核函数为径向基函数rbf，n_features表示特征数量，当采样频率sfreq＝1000hz时，n_features＝125，惩罚因子c设定为15。

通过上述过程，可以得到三组通道数据在foi-epoch特征平面的svm分类得分分布featurescores，其维度为(3×33×64)。

设定得分阈值thresholdscores＝0.75，获得分类贡献较高的一批通道频率特征z0，并标记定制的通道特征频率fo，将z0输入svm分类器进行调参，训练得到最终的svm分类器模型；

步骤5：对步骤(4)中得到的三组通道数据foi-epoch特征平面的svm分类得分featurescores设定阈值thresholdscores＝0.75，筛选出每组通道数据大于阈值的通道频率特征，组成定制通道频率特征向量z0,并标记定制的通道特征频率fo。可以注意到，由于在线识别无法对时间特征进行选择，所以这里仅仅需要提取符合要求的频率特征。具体操作过程由图2所示。

令定制通道频率特征z为输入，以rbf为核函数的svm的参数(c,gamma)作为调参对象，选择在5×5foldcross-validation中auc得分最高的一组参数(co,gammao)，并将此时训练的svm模型作为在线识别的分类器。

步骤6：离线训练模型用于在线识别，即实时采集用户脑电数据，依次进行预处理、提取定制通道频率fo的数据并计算通道频率特征z、输入svm模型进行识别，最终得到识别结果反馈给用户。

步骤6.1：用户佩戴好采集装置后，跟随显示屏指示进行左右手握拳的运动想象。电脑实时接收采集设备输送的脑电信号，并在接收到触发信号后，对每个采样时刻的数据进行fir高通(1hz)和低通(45hz)滤波后，进行独立成分分析剔除伪迹噪声。

步骤6.2：对去噪信号进行定制特征提取和分类识别操作：

1)根据定制通道特征频率fo，对特定电极通道的信号做带通滤波，中心频率为f∈fo，带宽为1hz，提取全部通道特征频率的数据；

2)对(1)中提取的数据使用multitaper算法计算psd特征，即定制通道频率特征z；

3)将z输入svm分类器识别想象内容；

4)将识别结果输出给反馈装置，形成闭环。

需要说明的是上述实施例仅仅是本发明的较佳实施例，并没有用来限定本发明的保护范围，在上述技术方案的基础上做出的等同替换或者替代均属于本发明的保护范围。