基于脑源域空间的MI-EEG识别方法与流程

本发明属于基于脑源空间的运动想象脑电信号(mi-eeg)识别与处理技术领域，具体涉及：采用标准化低分辨率脑电磁断层扫描成像(sloreta)方法将头皮层脑电信号逆变换到脑皮层，基于数据驱动(data-driven)及连续小波变换(cwt)进行偶极子优选及最优时间段的确定，并采用一对一共空间模式算法(ovo-csp)以及支持向量机(svm)实现mi-eeg在脑源空间的特征提取和分类。

背景技术：

脑机接口bci(brain-computerinterface)是一种不依赖于外周神经和肌肉组织等常规大脑信息输出通路，而运用工程技术手段在大脑和计算机或者其他机电设备之间建立连接，实现“让思想变成行动”的对外信息交流和控制的新途径。因此，bci技术在多种医疗与康复领域中获得广泛应用。bci主要分为侵入式和非侵入式两类，特别是非侵入式bci，使用头皮上的电极测量脑电信号(eeg)以反映大脑皮层活动，且该信号因具有采集便捷、价格低廉、时间分辨率高等特点，在bci系统中广泛采用。

运动想象脑电信号(mi-eeg)具有空间分布的特点，即不同运动想象任务，其在脑皮层对应的激活区域不同，因此由头皮层测得的脑电信号在空间分布上也会产生差异。而共空间模式算法(commonspatialpatterns,csp)能够利用mi-eeg在空间上具有分布差异性的特点，通过设计一组空间滤波器，使不同想像任务下的脑电信号滤波后的方差差异最大化，从而达到提取mi-eeg的特征、提高其可分离性的目的。然而，mi-eeg信号容易受到来自传感器以及头骨传导结构等噪声的严重污染，导致信噪比(snr)降低；同时，因体积传导效应及头皮电极数量的限制，使得头皮层测得的eeg信号的空间分辨率较低。这些因素均会造成头皮层mi-eeg在传感器域难以获得高质量的特征信息。

脑电源成像(esi)指的是利用头皮层记录的脑电数据对其在脑皮层源信号进行估计的一种常用技术，即将头皮层传感器空间转换为脑皮层源空间，获得脑皮层偶极子的源分布估计。近年来有研究表明，在脑皮层空间上应用csp算法对偶极子进行特征提取，对降低脑电信号在头骨传导过程中噪声的干扰、弥补头皮层脑电信号空间分辨率不足有一定作用，取得了比头皮层mi-eeg更好的特征提取效果。

然而，采用esi与csp相结合对运动想象脑电信号进行特征提取仍存在以下问题：

(1)求解脑电逆问题得到的不同类别运动想象任务偶极子时间序列杂乱无章，很少有规律可循，尤其是复杂的同侧肢体运动想象时，相应的皮层激活区域可能存在重叠，或者彼此非常接近，仅从时域上分析偶极子偶极矩时间序列，没有突显与运动想象有关的节律特点，因此很难提取到可分性较好的特征；

(2)脑源成像后得到的偶极子数量庞大，对所有偶极子进行特征提取将会造成特征信息冗余，大量与想象任务不相关的特征混杂其中，影响分类准确率；而且通过脑电正问题计算零导域矩阵时，由于实验设备的限制，无法对每个受试者进行fmri扫描，通常基于通用模板采用数值计算方式构建头模型。因数值方法存在近似性，这使求解正问题时对于不同受试者无法自适应地获得头皮电极和皮质神经元之间的精确对应关系。因此，通过脑电逆问题求得偶极子分布估计同样存在误差，从而影响基于神经生理学分区方法(brodmann分区)选取roi区域的准确性；

(3)对整个运动想象期间或者某个特定时段的偶极子进行分析，没有考虑在运动想象期间偶极子的时变特性，及不同受试者和不同运动想象任务造成的偶极子变化的差异性对分类识别率的影响。

技术实现要素：

针对上述不足之处，本发明对现有的esi与csp相结合对mi-eeg进行特征提取的算法进行了改进，提出了一种基于偶极子和连续小波变换的脑源域mi-eeg识别方法。具体涉及：

(1)从脑源域对偶极子时间序列进行时频分析，并采用ovo-csp算法对偶极子进行特征提取，实现脑皮层偶极子时域、频域、空域的特征有机融合，改善体积传导效应以及头皮mi-eeg空间分辨率低带来的不利影响。

(2)采用data-driven方法选择roi区域，初步减小偶极子量级；进而，对偶极子进行时频分析，精选偶极子，再次减小偶极子数量，精确地获得个性化的偶极子数量与偶极子信息，避免特征信息冗余的同时，提高计算效率。

(3)根据各个受试者偶极子的平均小波系数功率曲线图，选择不同想象任务之间差别较大的时间段作为最优时间段，充分体现不同受试者在不同想象任务类别下偶极子的时变特性及激活时间存在的个体差异性，进一步提高分类精度。

因此，本发明采用的技术方案为：基于脑源域空间的mi-eeg识别方法，首先对头皮层采集的脑电信号进行预处理，利用共平均参考降低原始信号基线漂移现象，同时对脑电信号进行8-32hz带通滤波，得到与运动节律相关的脑电信号；然后，基于sloreta的脑源成像逆算法将预处理之后的脑电信号转换到脑皮层空间，得到偶极子偶极矩时间序列，并根据脑皮层偶极子偶极矩幅值分布图对偶极子进行初选；进而，利用cwt对初选的偶极子进行时频分析，根据α波段小波系数功率序列精选偶极子并确定最优时间段；最后，在最优时间段内采用ovo-csp算法对精选的偶极子的α波段小波系数功率序列进行特征提取，并利用svm进行特征分类，采用十次十折交叉验证得到平均分类识别率，验证本方法的正确性和有效性。

基于上述分析，本发明的具体实现步骤如下：

s1头皮层eeg信号预处理；

s1.1假设为第i次实验采集的运动想像期间的头皮脑电信号，其中，m∈{1，2，3，4}表示运动想象任务的类别标签，i＝{1，2，3，…，nm}，nm表示采集实验次数；nc代表导联数；ns代表采样点数。

s1.2针对每个导联，选取运动想象前2s，即静息态时段内的脑电信号作为基准参考线，对xm，i逐行进行共平均参考处理，获得共平均参考处理信号x′m，i。

s1.3根据神经生理学理论，构造8-32hz带通滤波器对脑电信号x′m，i进行滤波，获得与μ、β节律(μ节律为8-12hz，β节律为13-30hz)相关的运动想象脑电信号，记为(r代表实数集)。

s2基于sloreta算法的脑电逆变换；

s2.1脑源成像正问题模型建立；

脑源成像正问题，指在人受到外部环境的刺激或进行特定的思维活动时，大脑皮层中的神经元细胞产生突触后电位，经过脑容积效应和颅骨的传导，产生了头皮层所测量得到的脑电信号，该脑电信号的数学模型建立如下：

m＝gd+n(1)

式中，d表示脑皮层电流偶极子的偶极矩，模拟了脑皮层中的神经元细胞的分布和其突触后电位的变化情况；g表示零导域矩阵，代表了脑皮层源空间与头皮层脑电信号的一种映射关系，通过bem边界元法在多层头模型上计算得到；m表示头皮层脑电信号；n为传导噪声。

s2.2根据正问题模型，假设偶极子源模型为分布式偶极子模型，方向垂直于脑皮层，数量为nd，位置已知且固定；令m＝mm，i，采用sloreta算法计算公式(1)的逆解，求得nd个偶极子源估计序列其中，第k个偶极子在t时刻的偶极矩幅值依下式计算：

式中，为由最小范数估计算法mne估计的偶极子源序列中第k个偶极子在t时刻的偶极矩幅值，由下式确定：

式中，α为正则化参数，为nc×nc的单位矩阵，为零导域矩阵。为的简写；

为偶极子电流密度方差矩阵的第k个对角元素，的计算公式如下：

s3基于连续小波变换cwt的偶极子选择与时频分析；

s3.1对于某类想象任务m，绘制偶极子时间序列曲线及脑皮层偶极矩幅值分布图，找到运动想象期间偶极子的偶极矩幅值最大值发生的时刻tmax，将此最大值记为dmax；在t＝tmax时选择偶极矩幅值超过的偶极子，其数量记为则该个偶极子覆盖的区域构成第m类想象任务的roi区域，记为qm；进而，求得四类想象任务的roi区域q＝q1∪q2∪q3∪q4，q内偶极子即为初选的偶极子，其数量记为

s3.2对个偶极子源估计时间序列依次进行cwt，得到频率为f时的偶极子小波系数序列，记为：其中，且f∈z，fs为采样频率。则第q个偶极子在t时刻的小波系数记为cm，i，f(q，t)。

s3.3依式(5)计算个偶极子在覆盖α波段的5个频率，即f＝8、9、10、11、12hz时对应小波系数序列平均值：

其中，cm，i(q，t)代表第q个偶极子在t时刻的小波系数在α波段的平均值，为α波段小波系数平均值矩阵；

s3.4对于cm，i(q，t)，利用公式(6)取平方运算求得其功率值：

pm，i(q，t)＝|cm，i，(q，t)|²(6)

其中，pm，i(q，t)代表α波段下第q个偶极子在t时刻的小波系数功率值，为α波段小波系数功率矩阵。

s3.5偶极子的精选及最优时间段的确定；

s3.5.1对每类想象任务m的所有样本，i＝{1，2，3，…，nm}，依式(7)求取偶极子q的平均小波系数功率矩阵

并绘制各偶极子的平均小波系数功率曲线获得想象任务m下的平均小波系数功率曲线图，这里，综合考虑四类想象任务，即m＝{1，2，3，4}，选取平均小波系数功率曲线差异较明显的最优时间段t，其包含的采样点数为nt，即

s3.5.2在时间段t内，基于同类想象任务m的平均小波系数功率矩阵依照采样点求第q个偶极子的小波系数功率均值并将降序排列，对每类任务m，选取前个偶极子；再从各类任务的前个偶极子中精选出共同的偶极子，将其数量记为

s3.6基于将精选的个偶极子按最优时间段t的采样点范围进行截取，获得最优小波系数功率矩阵，记为

s4基于一对一共空间模式ovo-csp的特征提取；

s4.1构造ovo-csp空间滤波器：对于四类想象任务中的任意两类运动想像脑电信号，根据经典csp算法构建空间滤波器a＝{1，2，3，4，5，6}，进而获得ovo-csp算法空间滤波器：

s4.2特征提取：将s3.6得到的作为ovo-csp空间滤波器sf的输入，并对滤波器输出信号的方差取对数运算和标准化处理，得到特征向量构成m类想像任务的特征集

s5基于支持向量机svm的分类；

利用svm对四类运动想象脑电信号特征集进行特征分类，采用十次十折交叉验证得出平均分类准确率。

与现有技术相比，本发明具有以下优点：

(1)针对运动想象脑电信号在头皮层具有低空间分辨率、低信噪比的缺点，本发明基于神经生理学理论基础，利用sloreta算法将头皮层脑电信号映射到脑皮层高维空间，并采用cwt对偶极子进行时频分析，选用ovo-csp算法提取偶极子小波系数功率特征，有利于降低脑电信号传导过程中体积传导效应及采集过程中的噪声干扰，综合利用偶极子的时-频-空特点，在脑源域实现了mi-eeg的多域特征融合。

(2)本发明克服了基于神经生理学分区选择roi区域不准确的缺点，采用基于受试者的data-driven的方法初选偶极子，再根据对脑皮层偶极子时频分析得到的小波系数功率序列精选偶极子，大量排除了与运动想象任务无关的偶极子冗余信息，减小计算量的同时，增强了对受试者的自适应性。

(3)本发明根据各个受试者偶极子的平均小波系数功率曲线图选择不同想象任务之间差别较大的时间段作为最优时间段，体现了不同受试者和不同运动想像任务下偶极子激活时间的差异性与偶极子的时变特性，有利于分类正确率的提高。

本发明适合多导联、复杂运动想象任务的bci系统，将为bci技术提供更加广阔的应用前景。

附图说明

图1为本发明采用的算法流程图；

图2.1为10-20导联系统电极分布图；

图2.2为mi-eeg采集实验时序图；

图3.1为受试者s1左手想象任务下tmax时刻的脑皮层偶极子偶极矩幅值分布图；

图3.2为受试者s1右手想象任务下tmax时刻的脑皮层偶极子偶极矩幅值分布图；

图3.3为受试者s1脚想象任务下tmax时刻的脑皮层偶极子偶极矩幅值分布图；

图3.4为受试者s1舌头想象任务下tmax时刻的脑皮层偶极子偶极矩幅值分布图；

图4.1为受试者s1左手想象任务下运动想象提示期与想象期的平均小波系数功率曲线图；

图4.2为受试者s1右手想象任务下运动想象提示期与想象期的平均小波系数功率曲线图；

图4.3为受试者s1脚想象任务下运动想象提示期与想象期的平均小波系数功率曲线图；

图4.4为受试者s1舌头想象任务下运动想象提示期与想象期内的平均小波系数功率曲线图；

具体实施方式

本发明具体实验是在windows8(64位)操作系统下的matlabr2017a仿真环境中进行的。

本发明采用的mi-eeg数据集来源于bcicompetitioniv的datasets2a公开数据库，由开发者使用国际标准的10-20导联系统下均匀分布的22个导联采集完成，采样频率为250hz，经过0.5-100hz的带通滤波。头皮层电极位置分布如图2.1所示。

采集试验时序图如图2.2所示，每次实验持续7.5s。0～2s为静息状态期，屏幕出现一个十字光标，同时t＝0s时发出短促的报警声；2s～3.5s为运动想象任务提示期，屏幕出现箭头，指向左、右、上、下，分别代表左手、右手、脚、舌头四种运动想象任务；3s～6s为运动想象期，受试者根据屏幕上的提示箭头进行运动想象；6s～7.5s为休息期，屏幕处于黑屏状态，受试者休息；随后进行下一次实验。datasets2a数据集共采集了九个受试者的运动想象脑电信号，每个受试者共576次实验(四种运动想象任务各144次)，一次实验共1875个采样点。

基于上述的运动想象脑电数据集，采用图1所示算法流程，以受试者s1为例，本发明具体的实施步骤如下：

s1四类mi-eeg信号预处理

s1.1根据每类任务类别标签(左手m＝1、右手m＝2、脚m＝3、舌头m＝4)提取受试者s1每类运动想象单次实验xm，i∈r^22×1875，其中i＝{1，2，3，…，144}，共获得576组mi-eeg。

s1.2截取每次实验0～2s静息状态期的脑电信号对mi-eeg进行基线校正和8-32hz带通滤波，得到预处理后的运动想象脑电信号，记为mm，i∈r^22×1875。

s2基于sloreta算法的脑电逆变换

s2.1选择ibm52模板头模型，采用边界元法(bem)计算正问题零导域矩阵g∈r^22×15002。

s2.2采用sloreta算法对步骤1.2得到的mi-eeg信号mm，i∈r^22×1875进行脑电逆变换，由式(2)求得15002个偶极子偶极矩时间序列估计

s3基于连续小波变换(cwt)的偶极子选择与时频分析

s3.1对于某类想象任务m，绘制偶极子时间序列曲线及偶极矩幅值幅值分布图，找到运动想象期，即第500到第1500个采样点之间偶极子的偶极矩幅值最大值发生的时刻tmax，将此最大值记为dmax，图3.1至图3.4分别展示了受试者s1在想象左手、右手、脚、舌头运动时，在tmax时脑皮层偶极子偶极矩幅值分布图。在t＝tmax时选择偶极矩幅值超过的偶极子，其数量记为则该个偶极子构成第m类想象任务的roi区域，记为qm；进而，可求得四类想象任务的roi区域q＝q1∪q2∪q3∪q4，q内偶极子即为初选的个偶极子，对于受试者s1，

s3.2对受试者s1的1285个初选偶极子源估计时间序列依次进行cwt，得到频率为f时的偶极子小波系数序列，记为：cm，i，f∈r^1285×1875。

s3.3依式(5)计算1285个初选偶极子在覆盖α波段的5个频率，即f＝8、9、10、11、12hz时对应小波系数平均值序列cm，i∈r^1285×1875。

s3.4对于s3.3得到的小波系数平均值序列cm，i∈r^1285×1875，利用公式(6)计算小波系数功率矩阵pm，i∈r^1285×1875。

s3.5偶极子的精选及最优时间段的确定

s3.5.1对每类想象任务m的144组实验，依式(7)求取平均小波系数功率矩阵然后，绘制各偶极子的平均小波系数功率曲线，获得想象任务m下的平均小波系数功率曲线图。图4.1至图4.4分别展示了受试者s1左手、右手、脚、舌头四类想象任务下的平均小波系数功率曲线图。从四类想象任务的平均小波系数功率曲线图可以看出，四类波形差异较明显的时间段为2.5～4s，因此受试者s1的最优时间段t＝2.5～4s，其包含的采样点数为nt＝375。

s3.5.2在最优时间段t＝2.5～4s内，基于同类想象任务m的平均小波系数功率矩阵依照采样点求第q个偶极子的小波系数功率均值并将降序排列，对每类任务m，选取前个偶极子；再从各类任务的前个偶极子中精选出共同的偶极子，将其数量记为对于受试者s1，

s3.6基于s3.4得到的576组小波系数功率矩阵pm，i∈r^1285×1875，m∈{1，2，3，4}，i＝{1，2，3，…，144}，将精选的32个偶极子按最优时间段t＝2～3.5s的采样点范围进行截取，获得最优小波系数功率矩阵，记为rm，i∈r^32×375。

s4基于一对一共空间模式(ovo-csp)的特征提取

采用ovo-csp算法对s3.6得到的576组四类运动想象任务最优小波系数功率矩阵rm，i进行特征提取，得到特征向量fm，i∈r^12×1，构成m类想像任务的特征集fm＝{fm，1，fm，2，…，fm，288}∈r^12×144，m∈{1，2，3，4}。

s5基于支持向量机(svm)的分类

利用svm对四类运动想象脑电信号特征集fm∈r^12×144(m∈{1，2，3，4})进行特征分类。

针对其他受试者，即s2、s3、…、s9的mi-eeg，采用以上类似的过程进行处理，获得十次十折交叉验证平均识别率如表1所示。

表1平均分类识别率

从表1可以看出，9个受试者采用本发明中的脑源空间mi-eeg识别方法获得的平均分类识别率达到了82.0％，与基于全部偶极子或采用brodmann方法选择roi，且无时间段优选和偶极子时频分析的识别结果相比，分别提高了7.4％和22.3％，显示了本发明的显著优势。