一种基于个体自适应的运动想象脑电特征表征方法与流程

1.本发明属于脑电信号处理领域，具体涉及一种基于个体自适应的运动想象脑电特征表征方法。


背景技术：

2.随着疾病、自然灾害和安全事故等原因造成残疾人员数量的增加，以及人口老龄化的增速，各类失能人群逐年增长。运动失能是一种常见的失能症状，这类人群通常由于各种疾病或事故造成肌肉组织或者神经传输通道严重受损并导致局部甚至全身瘫痪，因此如何拓展运动失能患者与外界的联系，实现生活自理成为全社会重点关注的问题之一。
3.脑-机接口(brain-computer interface,bci)技术让人类与计算机之间出现了一道通信的桥梁，失能患者可以通过这种交互方式与外界设备(轮椅等)建立连接，逐渐实现行动自主化。由于脑电信号(electroencephalography,eeg)较复杂，解码理论和方法还处于探索阶段。当前bci系统中广泛使用的运动脑电模式有诱发式和自发式，其中诱发式有稳态视觉诱发电位(steady-state visually evoked potential,ssvep)和p300，自发式的有运动想象(motor imagery,mi)。相比较诱发式eeg模式，自发式的mi无需外界刺激信号激发，更贴近生活，应用范围更广泛，受到了较高关注。科学研究表明，当人们在想象左手或者右手运动时，大脑皮层相关的运动感知区域会发生显著的变化，具体表现形式为对侧的运动感知区域在(8～12hz)节律和(16～28hz)节律上会出现信号振幅降低，即事件相关去同步化(event related desynchronization,erd)现象；而同侧的运动感知区域会出现相对明显的静息状态，即事件相关同步化(event related synchronization,ers)现象；因此利用这一现象，可以根据想象不同身体部位的运动进行eeg信号分类识别，从而产生不同的控制命令。
4.近年来，研究者们提出了大量关于运动脑电辨识的相关方法，包括基于时域、频域、时频域和空域等方面的特征提取算法以及用于分类的各种分类器。其中一种针对二分类模型高效的空间滤波算法，共空间模式(common space pattern,csp)得到了较广的应用，它通过寻找出一组最优的空间线性变换，从而获得两个类别之间差异最大化的空间特征，进而达到最优的区分效果。然而，传统的csp只考虑了脑电信号的空域信息，未能考虑到脑电信号在时域和频域内的可分性，且csp需要多通道分析，易受噪声信号干扰；因此许多研究者提出了改进的方法，比如正则化csp(regularized csp,rcsp)、稀疏csp(space csp,scsp)等方法。运动脑电信号辨识研究，其核心工作之一是相关通道和频段的选择。das等人通过有效值和z值的选择方法，挑选出区分度较高的通道，同时去除了较为嘈杂的部分。齐飞飞等人根据不同受试者的特点提出基于时空滤波的通道选择，利用脑电数据的时空信息自动识别出所需的通道，挑选出不同个体的最佳通道集，一定程度上减小了个体之间的差异性。par等人在频段划分的基础上，引入时域特征和相关系数来进行滤波器组共空间模式的选择，弥补了csp在单一频段上的缺陷。目前在运动脑电信号辨识方面，研究者们在信号通道选择和频段选择上未能充分考虑个体的差异性。不同个体在运动想象时，erd/ers现象
较为明显的通道往往不尽相同，即使同一个体最佳通道的最佳频段也受自身内外在因素的影响，因此最终特征差异性较大，影响分类辨识准确性。


技术实现要素：

5.发明目的：针对以上问题，本发明提出一种基于个体自适应的运动想象脑电特征表征方法，通过引入基于relieff的自适应通道选择方法和子频带优劣等级化共空间模式(sub-band hierarchy csp,sbh-csp)，完成各频带最优通道集挑选及最佳特征数量的获取，进而对运动脑电达到最有效的表征。
6.技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种基于个体自适应的运动想象脑电特征表征方法，具体包括如下步骤：
7.步骤1，获取受试者多类别多通道运动想象脑电信号并进行预处理，利用巴特沃兹滤波器对预处理后的脑电信号带通滤波后进行频带划分，得到多个子频带下的脑电信号样本数据；
8.步骤2，对于步骤1得到的各子频带下的脑电信号样本数据，通过relieff的自适应通道选择方法，完成各频带下最优通道集的选择；
9.步骤3，对于步骤2所述各频带下的最优通道集，利用sbh-csp方法完成各频带的有效特征提取，将提取到的有效特征进行特征融合，得到各频带下的融合特征；
10.步骤4，将各频带下的融合特征按照一定比例划分为训练特征集和测试特征集；利用svm算法构建网络模型，并利用训练特征集对所述网络模型进行训练，利用测试特征集对网络模型进行测试，得到各频带下的融合特征的分类识别结果。
11.进一步地，所述步骤1的具体方法如下：
12.步骤1.1，获取受试者多类别多通道运动想象脑电信号并进行预处理，包括小波去噪、共平均参考和基线校正；
13.步骤1.2，利用6阶巴特沃兹滤波器对预处理后的脑电信号进行8～32hz的带通滤波，并对滤波后的脑电信号进行频带平均划分，得到6个子频带脑电信号；
14.步骤1.3，对每一个子频带脑电信号均进行样本划分，得到每一个子频带下的高维脑电信号样本数据；其中，第i个子频带的高维脑电信号样本数据表示为：
[0015][0016]
式中，为第i个子频带脑电信号的第j个通道，i＝1,2,3,4,5,6，j＝1,2,
…
,n，t＝1,2,
…
,n；n为脑电信号采样点总数，n为脑电信号通道总数，t为第t个采样点。
[0017]
进一步地，步骤2所述通过relieff的自适应通道选择方法，完成各频带下最优通道集的选择，具体方法如下：
[0018]
步骤2.1，对于高维脑电信号样本数据，提取每个通道的5个时频域观察特征，即频带能量、功率、功率谱密度、瞬时能量和小波包能量；
[0019]
步骤2.2，利用relieff的自适应通道选择方法计算每个观察特征的权重，并根据权重由大到小进行排列；
[0020]
步骤2.3，根据观察特征权重大小的通道分布情况，挑选出每个子频带符合要求的通道，得到各频带下最优通道集。
[0021]
进一步地，所述步骤2.1中，频带能量的计算公式为：
[0022][0023]
式中，b为频带能量，n为脑电信号采样点总数，t为第t个采样点，为第i个子频带脑电信号的第j个通道，i＝1,2,3,4,5,6，j＝1,2,
…
,n，t＝1,2,
…
,n；n为脑电信号通道总数，j表示虚数单位，ω表示角频率；
[0024]
功率的计算公式为：
[0025][0026]
式中，t为功率；var()为协方差函数；
[0027]
功率谱密度的计算公式为：
[0028][0029]
式中，p'为功率谱密度；
[0030]
瞬时能量的计算公式为：
[0031][0032]
式中，i为瞬时能量，*表示卷积运算；
[0033]
小波包能量的计算公式为：
[0034][0035]
式中，w为小波包变换能量，表示经过小波包分解重构后的信号。
[0036]
进一步地，步骤2.2所述利用relieff的自适应通道选择方法计算每个观察特征的权重，具体方法如下：
[0037]
对步骤1所述脑电信号样本数据中的每一个样本标注类别标签；
[0038]
对于任一观察特征，首先，从步骤1所述脑电信号样本数据中随机选择一个样本rq；
[0039]
然后，在与样本rq相同类别标签的样本数据中选择k个最近邻样本，并在与样本rq不同类别标签的样本数据中选择k个最近邻样本；
[0040]
最后，根据权重单次迭代公式不断更新所述观察特征对应的权重，循环计算m次，直至所述样本数据中全部样本均完成计算，得到该观察特征的最终权重；所述权重单次迭代公式为：
[0041][0042]
式中，w
q 1
(a
l
)为第q 1个样本的第l个观察特征a的权重；wq(al)为第q个样本的第l个观察特征a的权重；a
l
为第l个观察特征a；hr为与样本rq相同类别标签的k个最近邻样本中
的第r个样本；p(c)为在样本中属于类别c的样本所占比值；p(class(rq))为与样本rq相同类别标签的样本占样本数据总样本的比值，其中class(rq)为样本rq的标签；mr(c)为与样本rq不同类别标签k个最近邻样本中的第r个样本，r＝1,2,3
…
k；函数diff(a,r1,r2)表示样本r1和样本r2在观察特征a上的差，计算方法如下：
[0043][0044]
r1[a]表示在样本r1上观察特征a的值；r2[a]表示在样本r2观察特征a上的值；max(a)表示在样本r1和r2上观察特征a的最大的值；min(a)表示在样本r1和r2上观察特征a的最小的值。
[0045]
进一步地，步骤3所述利用sbh-csp方法完成各频带的有效特征提取，具体包括如下步骤：
[0046]
步骤3.1，针对各频带构建各自对应的csp模型，所述csp模型具体为：
[0047]
记x1和x2为步骤1中任意两种类别运动想象脑电信号矩阵；
[0048]
求上述两种类别运动想象脑电信号矩阵归一化后的协方差矩阵ri和r
ii
分别为：
[0049]049][0050]
式中，上标t表示矩阵的转置，trace( )表示矩阵的迹；
[0051]
求解混合空间协方差矩阵r：
[0052][0053]
式中，分别表示两种类别运动想象脑电信号的平均协方差矩阵；
[0054]
对混合空间协方差矩阵r进行特征值分解：
[0055]
r＝uλu
t
[0056]
式中，u是矩阵r的特征向量矩阵，λ是混合空间协方差矩阵r对应的特征值构成的对角矩阵；其中且λ1,λ2…
λ
l
分别为l个特征值，矩阵u中每个列向量u
α
是矩阵r对应于特征值λ
α
的特征向量，α＝1,2,3
…
l；
[0057]
将求解得到的对角矩阵λ按照其特征值进行降序排列后，求白化矩阵p：
[0058][0059]
式中，λ-1
表示矩阵λ的逆矩阵；
[0060]
分别对两种类别运动想象脑电信号矩阵归一化后的协方差矩阵ri和r
ii
进行白化处理并进行主成分分解：
[0061][0061][0062]
式中，s1表示矩阵ri白化处理后的矩阵；s2表示矩阵r
ii
白化处理后的矩阵；b1表示矩阵s1的特征向量矩阵，b2表示矩阵s2的特征向量矩阵；λ1、λ2分别对应矩阵s1和s2特征值所构成的对角矩阵；
[0063]
矩阵ri和r
ii
白化处理后，满足b1＝b2，λ1 λ2＝e，则得到s1的最大特征值所对应的特征向量使s2有最小的特征值；其中e表示单位矩阵；
[0064]
则求得空间滤波器所对应的投影矩阵w为：
[0065]
w＝b
t
p
[0066]
式中，b＝b1＝b2；
[0067]
将步骤1所述对应频带下的样本脑电信号矩阵通过空间滤波器w，得到特征向量矩阵z：
[0068]
z＝wx
[0069]
式中，x为不同类别运动想象脑电信号矩阵，x＝x1，x2。
[0070]
接着在特征向量矩阵z中，提取前m个和后m个特征值所对应的特征向量，并构成空间滤波信号z
p
，p＝1,2,
…
,2m，并归一化求得特征向量：
[0071][0072]
式中，y为对应频带下的脑电信号样本归一化后的特征向量，z
po
表示空间滤波信号z
p
的第po行信号；var()为协方差函数；
[0073]
步骤3.2，计算各频带下特征向量矩阵z中对应的前四个特征值之和，并选择特征值之和最大的两个频带，设置这两个频带的特征参数为4；计算剩余四个频带下的特征向量矩阵z对应的前三个特征值之和，选择特征值之和最大的两个频带，设置这两个频带的特征参数为3；将最终剩余的两个频带设置其特征参数为2；
[0074]
步骤3.3，利用csp对各频带下脑电信号进行特征提取，并根据步骤3.2中各频带特征参数获得有效特征数量集。
[0075]
有益效果：与现有技术相比，本发明技术方案具有以下有益技术效果：
[0076]
本发明基于个体差异和频带信号特点，提出基于relieff的自适应通道选择方法。通过提取各频带5个时频域观察特征来表征通道的erd/ers情况，运用relieff算法对各频带信号通道进行有效评估，避免了全通道信息噪声多、普适性差等缺点；在使用csp特征提取的同时根据频带优劣性确定各自频带的最佳特征参数m，在一定程度上优化了特征筛选，更能实现运动想象脑电状态的有效表征，从而提高运动脑电的辨识率。
附图说明
[0077]
图1是一种实施例下本发明所述一种基于个体自适应的运动想象脑电特征表征方法流程图；
[0078]
图2是一种实施例下单次运动想象任务时序实验图；
[0079]
图3是一种实施例下基于relieff的自适应通道选择方法流程图；
[0080]
图4是一种实施例下sbh-csp特征提取流程图。
具体实施方式
[0081]
下面结合附图和实施例对本发明的技术方案作进一步的说明。
[0082]
本发明所述的一种基于个体自适应的运动想象脑电特征表征方法，在一定程度上可以减小因个体差异导致脑电特征识别率低的问题，如图1所示，具体包括如下步骤：
[0083]
本发明采用公用数据库中的eeg movement/imagery dataset公共数据集进行说明。该数据库由美国纽约卫生署wadsworth中心组织开发，使用bci2000采集系统，共采集109人的64导联通道数据，信号采样频率为160hz。实验执行左拳、右拳、双拳和双脚这四种类别任务，每个受试者在上述四类运动想象次数共90次。每次实验持续约8s，其中0～2s时屏幕黑屏，受试者处于放松准备状态，2～6s时在电脑屏幕上下左右区域会随机出现目标作为提示(左：左拳运动；右：右拳运动；上：双拳运动；下：双脚运动)，受试者根据屏幕提示完成运动想象实验，然后进入2s的休息时间，随后进入下一个实验，具体单次实验步骤时序示意图如图2所示。
[0084]
单次实验数据运动想象时间为4s，运动想象间有640个采样点，然后对每次运动想象数据进行预处理及样本划分。首先通过小波去噪、共平均参考、基线校正等进行前期预处理；再使用6阶巴特沃斯滤波器对预处理后的eeg信号进行8～32hz的带通滤波，并将其划分为6个子频带：a1频带(8～12hz)、a2频带(12～16hz)、a3频带(16～20hz)、a4频带(20～24hz)、a5频带(24～28hz)和a6频带(28～32hz)；最后对于每一个子频带的样本，采用160点窗口重叠50％进行样本滑窗划分。因此，对于每一位受试者的每个频带得到共630个样本数据，其数据存储形式为160x64x630(数据点x通道x样本数)。
[0085]
将所得的样本按照7:3划分成训练集和测试集。
[0086]
步骤2、基于对步骤2所得的信号样本，使用基于relieff的自适应通道选择方法，完成各频带最优通道集选择。
[0087]
基于relieff的自适应通道选择方法如图3所示，主要包括3个步骤：
[0088]
a1、针对每个频带样本，提取64个通道的5个时频域观察特征，即频带能量(band energy，be)，功率(power)、功率谱密度(powerspectral density,psd)、瞬时能量(instantaneous energy，ie)和小波包能量(wavelet packet transform,wpt)。
[0089]
其中定义5种时频域观察特征公式频带能量b，功率t，功率谱p'，瞬时能量为i，小波包变换能量w分别为：
[0090][0091][0092]
[0093][0094]
其中*表示卷积运算。
[0095][0096]
其中，表示信号经过小波包分解重构后的信号。
[0097]
初步得到每个频带样本320维特征。
[0098]
a2、使用relieff算法计算各观察特征的权重，根据权重大小进行由大到小的排。relieff算法具体计算过程如下：
[0099]
对任意特征，首先从训练集中随机选择一个样本rq，然后从和rq同类(类别标签相同)的样本中选择k个最近临样本(与rq距离最近)，从与rq不同类的样本中选择出k个最近邻样本；最后根据权重单次迭代公式不断更新该特征对应的权重，循环计算m次直至所有样本依次计算完毕，得到单个特征的最终权重。relieff算法单次特征权重迭代公式为：
[0100][0101]
其中wq(a
l
)为第q个样本中第l个特征a的权重；a
l
为第l个观察特征a；hr(r＝1,2,3
…
k)为与rq同类的k个最近邻样本中的第r个样本；p(c)为在训练样本中属于类别c的样本所占比值；p(class(ri))为与rq同类的样本占总样本的比值，其中cl
ass
(rq)为样本rq的标签；mr(c)(r＝1,2,3
…
k)为与rq不同类的k个最近邻样本中的第r个样本。函数diff(a,r1,r2)表示样本r1和样本r2在特征a上的差，其计算方法如下：
[0102][0103]
r1[a]表示在样本r1上观察特征a的值；r2[a]表示在样本r2上观察特征a的值；max(a)表示在样本r1和r2上观察特征a的最大的值；min(a)表示在样本r1和r2上观察特征a的最小的值。
[0104]
a3、挑选出各子频带权重前70％的观察特征，统计分析其所处通道分布情况，当通道分布特征数不少于3个时，该通道即为此子频带信息活跃通道，该通道保留，从而完成各子频带通道的挑选。
[0105]
步骤3、通过使用sbh-csp完成各子频带最有效特征提取，各频带特征融合后组成训练特征集和测试特征。其中参考图4，sbh-csp特征提取具体流程为：
[0106]
b1、使用训练集构建csp模型，其csp模式具体为：
[0107]
记x1和x2为任意两类运动想象下信号矩阵，首先求两类信号归一化后的协方差矩
阵ri和r
ii
分别为：
[0108][0109]
(1)式中：x
t
表示x矩阵的转置，trace()表示矩阵的迹。
[0110]
然后求混合空间协方差矩阵r：
[0111][0112]
式中：分别为两类信号的平均协方差矩阵。
[0113]
对混合空间协方差矩阵r进行特征值分解
[0114]
r＝uλu
t
ꢀꢀꢀ
(10)
[0115]
式中：u是矩阵r的特征向量矩阵，λ是对应的特征值构成的对角矩阵。其中且λ1,λ2…
λ
l
分别为l个特征值，矩阵u中每个列向量u
α
(α＝1,2,3
…
l)是矩阵r对应于特征值λ
α
的特征向量。
[0116]
将求解得到的对角矩阵λ按照其特征值进行降序排列，求得白化矩阵p为：
[0117][0118]
对ri和r
ii
进行白化处理：
[0119]
s1＝prip
t
,s2＝pr
ii
p
t
ꢀꢀꢀ
(12)
[0120]
然后对s1和s2做主成分分解，得到：
[0121][0122]
通过上面的式子可以证明矩阵s1的特征向量和矩阵s2的特征向量矩阵是相等的，即：
[0123]
b1＝b2＝b
ꢀꢀꢀ
(14)
[0124]
与此同时，两个特征值的对角矩阵λ1和λ2之和为单位矩阵，即：
[0125]
λ1 λ2＝e
ꢀꢀꢀ
(15)
[0126]
同时两类矩阵的特征向量对应的特征值相加总是为1，因此s1的最大特征值所对应的特征向量使s2有最小的特征值，反之亦然。白化eeg到与λ1和λ2中的最大特征值对应的特征向量的变换对于分离两个信号矩阵的方差是最佳的。因而求得空间滤波器所对应的投影矩阵w为：
[0127]
w＝b
t
p
ꢀꢀꢀ
(16)
[0128]
接着将原始eeg信号矩阵通过空间滤波器w，得到对应的特征向量矩阵z：
[0129]
z＝wx
ꢀꢀꢀ
(17)
[0130]
最后提取前m个和后m个特征值所对应的特征向量，构成空间滤波后的信号z
p
，p＝1,2,
…
,2m，并归一化求得特征向量：
[0131][0132]
其中y为对应频带下脑电信号样本归一化后的特征向量，z
po
为空间滤波信号z
p
的第po行信号，log(i)为对数运算符。
[0133]
b2、确定不同频带的优劣等级参数；根据csp原理，csp模型中λ越大与之对应特征向量的变换对于分离两类信号矩阵的方差是最佳的。因此。计算各频带下特征向量矩阵z中对应的前四个特征值之和，并选择特征值之和最大的两个频带，设置这两个频带的特征参数为4；计算剩余四个频带下特征向量矩阵z中对应的前三个特征值之和，选择特征值之和最大的两个频带，设置这两个频带的特征参数为3；将最终剩余的两个频带设置其特征参数为2。
[0134]
b3、使用csp对各频带脑电信号进行特征提取，并根据步骤b2中选择各频带特征参数，从而获得有效的特征集，然后将各频带特征进行融合。
[0135]
步骤4、采用ovo-csp的方法，并使用svm构建训练模型并完成分类识别以及分类识别性能评估。
[0136]
为了更好的说明本发明对于运动脑电特征表征的有效性，采用数据库中的前十个人(s001,s002,s003,s004,s005,s006,s007,s008,s009,s010)进行具体实验说明，参考表1所示，相关每类实验重复5次交叉验证以减小偶然性造成的影响。
[0137]
表1个体各频带通道选择结果
[0138][0139]
表1给出了基于relieff的自适应通道选择方法前10位受试者不同频带挑选出来的通道个数。从表中可以看出，采用本文所提方法每位受试者各频带所挑选出来的通道数大多为40左右，所得数量较为适中且所选通道会因个体特点和频带信号特点存在差异。
[0140]
表2各个体受试者不同方法辨识精度
[0141][0142]
表2给出了10位受试者在三种方法下的具体分类识别精度，从个体来看，采用relieff-sbhcsp方法后10位个体的辨识率皆有提升，相比较传统的csp前后精度提升最高和最低量分别为17.7％和2.8％，单个个体最高辨识准确率达到了89.7％。整体看来，通过综合采用本文所提的relieff-sbhcsp方法平均准确率达到了82.5％，提高了6.9％，由此可见本发明在运动脑电特征表征上具有明显的效果。
[0143]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。