复合下肢想象动作脑电的能量特征提取方法

摘要

本发明属于生物医学工程及计算机领域，涉及一种复合下肢想象动作脑电的能量特征提取方法，包括①脑电信号预处理。②通过经验模态分解方法，将复合下肢想象动作脑电分解为频率从高到低的各个固有振荡模式分量。③通过功率谱密度分析方法分析各固有振荡模式的频谱分布特征，确定反映脑电特征节律的特征振荡模式。④通过对特征振荡模式进行希尔伯特变换，得到反映其能量变化的去同步化特征。⑤模式识别：将特征振荡模式的能量变化去同步化特征作为分类器输入，实现复合下肢想象动作的模式识别。本发明充分考虑到脑电信号的非平稳性，最高识别率为87.8％，相对于传统方法的82.3％有显著提高。

权利要求

1.一种复合下肢想象动作脑电的能量特征提取方法，包括下列步骤：

①利用脑电导联电极采集复合下肢想象动作脑电信号；

②对所采集到的位于初级运动区及辅助功能区的脑电信号进行空间滤波，提高脑电信号的信噪比；

③针对经过步骤②处理后的脑电信号，进行经验模态分解，将复合下肢想象动作脑电信号分解为频率从高到低的各个固有振荡模式分量；

④对步骤③获得的各导联脑电的固有振荡模式分量计算功率谱密度，根据其功率谱密度的频率分布范围，确定脑电特征节律-alhpa节律对应的固有振荡模式分量，即所要获取的特征振荡模式；

⑤ 定 义 特 征 振 荡 模 式 的 事 件 相 关 去 同 步 化 系 数，

其中REPt1为复合下肢动作想象开始后t1时间内

的脑电信号平均幅度，REFt2为复合下肢动作想象开始前t2时间内的脑电信号平均幅度，对脑电信号特征振荡模式计算事件相关去同步化系数；

⑥选择基于径向基核函数的支持向量基分类器，将步骤⑤获得的各个脑电导联事件相关去同步化系数作为分类器输入样本，对脑电信号进行模式识别。

2.根据权利要求1所述的复合下肢想象动作脑电的能量特征提取方法，其特征在于，步骤②处理的脑电信号包括F3P、FzP、F4P、C5A、C3A、C1A、C2A、C4A、C6A、C5、C3、C1、Cz、C2、C4、C6、C5P、C3P、C1P、C2P、C4P、C6P、P3A、PzA、P4A导联信号。

3.根据权利要求1或2任意一项所述的复合下肢想象动作脑电的能量特征提取方法，其特征在于，步骤⑤对初级运动区及辅助功能区的脑电信号特征振荡模式计算事件相关去同步化系数，按照导联位置顺序构成事件相关去同步化系数向量f＝[λF3P，λFzP，…，λP4A]，这里λF3P，λFzP，，…，λP4A依次为导联F3P、FzP、F4P、C5A、C3A、C1A、C2A、C4A、C6A、C5、C3、C1、Cz、C2、C4、C6、C5P、C3P、C1P、C2P、C4P、C6P、P3A、PzA、P4A位置的事件相关去同步化系数。

4.根据权利要求1或2任意一项所述的复合下肢想象动作脑电的能量特征提取方法，其特征在于，步骤⑤取t1为动作想象期的4-7秒，t2为动作想象前的0-2秒。

5.根据权利要求1或2任意一项所述的复合下肢想象动作脑电的能量特征提取方法，其特征在于，将步骤⑤获得的事件相关去同步化系数向量f＝[λF3P，λFzP，…，λP4A]作为分类器输入样本，分类器公式为y＝sgn(W·f+b)，W是权函数，b是阈值，阈值b通过训练样本得到，y是两分类的结果标记±1，+1表示识别出左手左腿同时动作想象，-1表示识别出右手右腿同时动作想象。

说明书

复合下肢想象动作脑电的能量特征提取方法

技术领域

[0001] 本发明属于生物医学工程及计算机领域，涉及一种复合下肢想象动作脑电的能量特征提取方法。

背景技术

[0002] 脑-机接口(Brain-Computer Interface，BCI)是在人脑和计算机或其他电子设备之间建立一种不依赖于常规大脑输出通路(外周神经和肌肉组织)的直接信息交流和控制通道，是一种全新的人-机交互系统。最早被应用于脑-机接口系统的脑电信号主要是自发脑电信号，比如脑电中的alpha(α)波。但是这类脑电信号模式单一，无法真正做到“意识控制行动”，严重制约了脑-机接口系统的发展。近年来，各国学者对不同思维意识下脑电信号的研究逐步开展起来，这为脑-机接口的发展带来了新的曙光。

[0003] 已有研究表明：人在想象某个肢体动作时，与该动作相关的大脑运动皮层区域会发生与该动作实施时相似的电生理响应，如诱发电位(evoked potential，EP)或事件相关电位(event related potential，ERP)，特称之为想象动作电位(motor imaginary potentials)。临床上已通过功能性磁共振成像(functional magnetic resonance image，fMRI)观察脑局部血液图的方法确认：想象和实施动作时所激发的大脑运动皮层区域相同。

[0004] 最先发现想象动作电位的Jasper等人在研究脑电(Electroencephalograph，EEG)信号的过程中，注意到肢体运动准备或规划时能够引起皮层运动中枢大量神经细胞的活动状态变化，导致EEG中某些频率成分的同步增强或者同步减弱，即
所谓事件相关同步(event-related synchronization，ERS)和事件相关去同步
(event-relateddesynchronization，ERD)现象。Pfurtscheller和Aranibar通过实验证实了上述现象，指出ERD/ERS现象主要集中在EEG中alpha节律和beta节律段，并且提出了针对ERS/ERD现象的量化理论。近些年来，各国学者对各种想象动作电位进行了大量的思维模式(多围绕ERD/ERS现象)提取的研究，其中以舌部想象动作和左右手想象动作电位的模式提取最为常见。特别需要指出的是，先期开展这方面工作的奥地利Graz大学，已经成功地将提取到的左右手想象动作模式作为控制命令，操纵脑-机接口系统并试用于肢瘫患者的上肢康复训练中，取得了令人振奋的阶段性成果。

[0005] 但到目前为止，下肢想象动作电位的模式提取进展缓慢，判断准确度难以提高。其主要原因是：下肢运动所映射的脑皮层功能区为头顶部沟回内较狭小的区域，其空间结构的区分度已是十分有限，加之头皮电极提取的脑电信号存在极大弥散性和混叠性，非常不利于源信号获取和识别。这个关键因素导致了应用于上肢想象动作脑电的特征提取算法在下肢想象动作脑电的特征提取方面适用性有限。而在现实生活中，下肢想象动作脑电特征的有效提取，是实现真正完全由四肢瘫痪病人或截瘫患者的意志控制的下肢康复助行系统的核心技术：通过提取患者想象下肢动作的脑电信号，转换为相应的下肢康复助行系统的外部控制命令，以此帮助患者行走和肌肉刺激恢复。这种完全靠患者自主意志控制的系统，不仅可以实现这部分残疾病人正常生活的回归，而且可以提升患者的自信心，因此具有广泛的应用前景和重要的社会意义。由于脑-机接口实现康复助行系统的急迫需求以及单纯下肢想象动作脑电模式的困境，希望通过复合下肢想象动作脑电实现脑-机接口系统的研究脱颖而出，逐步成为脑-机接口实现康复助行系统的研究重点之一。经过近几年的发展，该研究方向已被认为是脑-机接口进一步发展的必由之路，相应地，复合下肢想象动作脑电的能量特征提取也成为实现下肢康复助行系统任务中极其关键的研究内容之一。

发明内容

[0006] 本发明的主旨是提出一种复合下肢想象动作脑电的能量特征提取方法，以此解决基于脑-机接口的下肢康复助行系统中的基本问题：实现截瘫患者自主控制下肢动作实现正常行走的关键动作，左右迈步动作的自主控制。本发明解决了复合下肢想象动作脑电的能量特征准确提取问题，从而为正确识别复合下肢动作模式，有效转换为应用于下肢康复助行系统的控制命令，实现截瘫患者的自主行走提供了有力支持。

[0007] 为此，本发明采用如下的技术方案：

[0008] 一种复合下肢想象动作脑电的能量特征提取方法，包括下列步骤：

[0009] ①利用脑电导联电极采集复合下肢想象动作脑电信号；

[0010] ②对所采集到的位于初级运动区及辅助功能区的脑电信号进行空间滤波，提高提高脑电信号的信噪比；

[0011] ③针对经过步骤②处理后的脑电信号，进行经验模态分解，将复合下肢想象动作脑电信号分解为频率从高到低的各个固有振荡模式分量。

[0012] ④对步骤③获得的各导联脑电的固有振荡模式分量计算功率谱密度，根据其功率谱密度的频率分布范围，确定脑电特征节律-alhpa节律对应的固有振荡模式分量，即本发明所要获取的特征振荡模式。

[0013] ⑤ 定 义 特 征 振 荡 模 式 的 事 件 相 关 去 同 步 化 系 数，其中REP为复合下肢动作想象开始后t1时间内的
脑电信号平均幅度，REF为复合下肢动作想象开始前t2时间内的脑电信号平均幅度，对脑电信号特征振荡模式计算事件相关去同步化系数；

[0014] ⑥选择基于径向基核函数的支持向量基分类器，将步骤⑤获得的各个脑电导联事件相关去同步化系数作为分类器输入样本，对脑电信号进行模式识别。

[0015] 作为优选实施方式，所述的复合下肢想象动作脑电的能量特征提取方法，步骤②处理的脑电信号包括F3P、FzP、F4P、C5A、C3A、C1A、C2A、C4A、C6A、C5、C3、C1、Cz、C2、C4、C6、C5P、C3P、C1P、C2P、C4P、C6P、P3A、PzA、P4A导联信号；步骤⑤对初级运动区及辅助功能区的脑电信号特征振荡模式计算事件相关去同步化系数，按照导联位置顺序构成事件相关去同步化系数向量f＝[λF3P，λFzP，…，λP4A]，这里λF3P，λFzP，…，λP4A依次为导联F3P、FzP、F4P、C5A、C3A、C1A、C2A、C4A、C6A、C5、C3、C1、Cz、C2、C4、C6、C5P、C3P、C1P、C2P、C4P、C6P、P3A、PzA、P4A位置的事件相关去同步化系数；取t1为动作想象期的4-7秒，t2为动作想象前的0-2秒；将步骤⑤获得的事件相关去同步化系数向量f＝[λF3P，λFzP，…，λP4A]作为分类器输入样本，分类器公式为y＝sgn(W·f+b)，W是权函数，b是阈值，通过训练样本得到，y是两分类的结果标记(±1)，+1表示识别出左手左腿同时动作想象，-1表示识别出右手右腿同时动作想象。

[0016] 本发明通过经验模态分解方法结合功率谱密度分析及希尔伯特变换方法获取alpha节律能量变化特征，该方法的最高识别率为87.8％，相对于传统方法的82.3％有显著提高，其优势在于充分考虑到脑电信号的非平稳性，并且基于动作想象脑电的事件相关同步/去同步特征是由群体神经元同步振荡产生的结论，因此在基于复合肢体动作的脑-机接口系统中将具有更为广阔的应用前景，这将为基于脑-机接口的下肢康复助行系统的发展提供有力支持。

[0017] 本发明致力于解决基于脑-机接口的下肢康复助行系统中的关键问题，即复合下肢想象动作脑电的能量变化特征提取问题。下肢想象动作脑电的能量特征的有效提取，是实现真正完全由四肢瘫痪病人或截瘫患者的意志控制的下肢康复助行系统的核心技术：通过提取患者想象下肢动作的脑电信号，转换为相应的下肢康复助行系统的外部控制命令，以此帮助患者行走和肌肉刺激恢复。这种完全靠患者自主意志控制的系统，不仅可以实现这部分残疾病人正常生活的回归，而且可以提升患者的自信心，因此具有广泛的应用前景和重要的社会意义。

附图说明

[0018] 图1特征振荡模式的能量变化特征提取过程；

[0019] 图2复合肢体想象动作实验时段分配图；

[0020] 图3脑电采集所采用的41导联分布示意图；

[0021] 图4左手左腿协同动作想象时的脑电信号；

[0022] 图5左手左腿协同动作想象脑电的经验模态分解结果(a)C3导联脑电的固有振荡模式(b)C4导联脑电的固有振荡模式；

[0023] 图6左手左腿协同动作想象时，C3导联脑电信号的固有振荡模式的功率谱密度分布(a)第一固有振荡模式的功率谱密度分布(主要频带22-24Hz)(b)第二固有振荡模式的功率谱密度分布(主要频带8-13Hz，属于alpha节律)(c)第三固有振荡模式的功率谱密度分布(主要频带4-8Hz)(d)第四固有振荡模式的功率谱密度分布(主要频带2-4Hz)

[0024] 图7左手左腿协同动作想象时，C4导联脑电信号的固有振荡模式的功率谱密度分布；(a)第一固有振荡模式的功率谱密度分布(主要频带22-25Hz)(b)第二固有振荡模式的功率谱密度分布(主要频带10-13Hz，属于alpha节律)(c)第三固有振荡模式的功率谱密度分布(主要频带2-8Hz)(d)第四固有振荡模式的功率谱密度分布(主要频带2-5Hz)[0025] 图8右手右腿协同动作想象时，导联C3，C4位置特征振荡模式的功率谱密度分布对比；

[0026] 图9导联C3，C4位置特征振荡模式的希尔伯特包络。(a)左手左腿同时动作想象；(b)右手右腿同时动作想象；

[0027] 图10基于经验模态分解方法和基于传统滤波方法的模式识别准确率对比。

具体实施方式

[0028] 针对复合下肢想象动作脑电的征提取难点：下肢运动所映射的脑皮层功能区为头顶部沟回内较狭小的区域，其空间结构的区分度十分有限，加之复合下肢想象动作涉及多个脑功能区，模式复杂，进一步突出了脑电信号的非平稳特性。这就导致了传统的基于脑电信号短时平稳假设的信号处理方法在下肢想象动作脑电的特征提取方面适用性有限。为此，本发明针对脑电信号的非平稳特性，并基于动作想象脑电的事件相关同步/去同步特征是由群体神经元同步振荡产生的结论，提出应用经验模态分解方法提取复合下肢想象动作脑电的特征振荡模式，进而获取特征振荡模式的能量变化特征，以此作为脑电信号模式识别的输入参数，实现复合下肢想象动作的有效识别。

[0029] 下面分几个方面对本发明的脑电信号提取方法做详细说明。

[0030] 1经验模态分解(Emprical Mode Decomposition，EMD)理论的提出及其适用性[0031] 在传统的Fourier分析中，频率被定义为整个分析数据长度中具有一定幅度的正、余弦函数。受这种固有观念的影响，人们在认识和接受瞬时频率的意义和概念时，总是从正、余弦函数的有关角度来分析。这样当人们定义局部频率值时就需要多于一个周期的正、余弦波动，基于这个逻辑，少于一个周期长度的信号将无法给出其频率的定义。而对于非线性和非平稳信号来说，其主要特征频率是时变的，即仅仅是在某一局部时间内存在或曾在某一时刻出现过，在描述频率随时间的变化关系上，Fourier变换显然已经无能为力。为了弥补Fourier变换对时变信号分析的不足，人们对原始信号加窗，认为在某个“窄带”内的信号是平稳的或近似平稳的，然后再对窗内的信号进行分析，如短时傅立叶变换、小波分析等，这些方法不同程度上对非线性和非平稳信号的时变性进行了描述，大大改进了Fourier变换的不足。但是由于受Heisenberg不确定原理的制约，在时间和频率上的分辨率不能同时达到最小，因此，所得的结果是窗内信号的平均结果，同样也摆脱不了Fourier变换的局限性。

[0032] 为了把信号的频谱分析精确到每一个时间点、频率点上，美国国家航空航天局，美国工程院院士Huang N E等人提出了经验模态分解(Emprical Mode Decomposition，EMD)法。1996年，Huang在一次国际学术会议上首次提出这种适于非平稳信号分析的新方法的设想-基于经验的模式分解方法。Huang认为，对于瞬态与非平稳现象，频率与能量一般都是时间的函数，因此需要给出瞬时频率和瞬时能量的定义。目前信号瞬时能量的概念已被广泛接受，但是瞬时频率的概念和意义却一直存有争议。当可以使数据解析化的Hilbert变换出现之后，人们根据Hilbert变换提供的能够完全表达原始数据幅度和相位的函数，给出了瞬时频率的统一定义，从定义可以看出瞬时频率是时间的单值函数，即在任意时刻只存在一个振荡模式。所以在使用瞬时频率这一概念时，对应的数据受到了一定的限制。这主要因为任何一个时刻，数据中可能包含多个振荡模式，此时Hilbert变换不能给出该信号完全的频率内容，所得到的结果只是多个振荡模式的平均效果，从而瞬时频率的意义变得模糊。为了从复杂信号中得到有意义的瞬时频率，Huang根据瞬时频率物理意义上的必要条件，提出把含有多个振荡模式的数据分解为满足一定条件的多个单一振荡模式分量的线性叠加，每个单一振荡模式分量又叫做一个基本模式分量，并提出了一种基于经验的模式分解方法。每一个单一模式分量都满足Hilbert变换的必要条件，使得用Hilbert变换求解信号的瞬时频率成为可能。

[0033] 经验模态分解法的意义在于：在信号分析中信号频率的定义是基于波形的局部特征和瞬时特征，它能在信号数据的每个时间点上，从点与点之间的变化特征来给出瞬时频率值，而不是需要多个振荡周期的波形才能给出一个频率值。经验模态分解法分析中若存在一个频率，仅仅表示该频率对应的信息在某一局部时间内存在或者曾在某一时刻出现过。这样无论从概念上还是从信号分析本质上来看，这种分析方法打破了传统频率思想，给出了一个全新的频率概念。经验模态分解法对信号分析具有重大的意义，同时也为非平稳信号处理提供了新思路，开辟了新途径。

[0034] 2经验模态分解算法

[0035] 为了能把一般数据分解成固有特征振荡模式，Huang N E提出了经验模态分解的方法。和以前几乎所有的分解方法不同，该新方法是直观的、直接的、后验的和自适应的，其分解的基函数立足于数据并且来源于数据本身。

[0036] 经验模态分解方法是建立在以下的假设之上的：(1)信号至少有两个极值点，一个最大值和一个最小值；(2)特征时间尺度是通过两个极值点之间的时间间隔来定义；(3)若数据缺乏极值点但有变形点，则可通过数据微分一次或几次获得极值点，然后再通过积分来获得分解结果。

[0037] 这种方法的本质是通过数据的经验特征时间尺度来获得其内在的振荡模式，然后据此分解数据。根据Drazin的经验，数据分析的第一步是人工观察数据，有两种能直接区分不同尺度振荡模式的方法：观察依次交替出现的极大、极小值点间的时间间隔；和观察依次出现的过零点的时间间隔。交替的局部极值点和过零点形成了复杂的数据：一个波动骑在另一个波动上，同时它们又骑在其他的波动上，依此类推，每个波动都定义了数据的一个特征尺度，这个特征尺度是内在的。采取依次出现的极值点的时间间隔作为振荡模式的时间尺度，因为这个方法对振荡模式不但有更高的分辨率，而且能应用于非零均值的数据，例如没有过零点的，全部数据点是正的或负的数据。为了把各种振荡模式依次从数据中提取出来，使用一种系统的方法，即经验模态分解方法，或形象地称之为“筛选”的过程，对实信号s(t)进行经验模态分解的步骤为：

[0038] 1)确定s(t)的所有极大值和极小值；

[0039] 2)根据极大值和极小值作三次样条差值来构造s(t)的上下包络线；

[0040] 3)根据上下包络线，计算出s(t)的局部均值(上下包络的平均值)m1(t)，以及s(t)和m1(t)的差值h1(t)＝s(t)-m1(t)；

[0041] 4)以h1(t)代替原始信号s(t)，重复以上三步骤k次，直到所得的平均包络趋于零为止(h1，(k-1)(t)与h1，k(t)之间的方差小于设定值)，即认为h1，k(t)是一个IMF分量，记c1(t)＝h1，k(t)，r1(t)＝s(t)-c1(t)，s(t)＝r1(t)；第一个IMF分量代表原始数据中最高频率的分量。将原始数据序列s(t)减去第一个IMF分量c1(t)，可以得到一个去掉高频分量的差值数据序列r1(t)。

[0042] 5)对r1(t)重复以上四步平稳化处理过程，可以得到第二个IMF分量c2(t)，如此重复下去直到最后一个差值序列rn(t)不可分解为止(rn(t)小于一设定值，或者变成一个单调函数时)，原始信号的经验模态分解结束，得到s(t)的分解式如下：

[0043]

[0044] 其中rn(t)为残余函数，代表原始数据的趋势或均值。

[0045] 由于每一个IMF分量都是代表一组特征尺度的数据序列，因此这个平稳化的处理过程实际上是将原始数据序列分解为不同特征波动的叠加。需要说明的是，每一个IMF分量既可以是线性的也可以是非线性的。

[0046] 经验模态算法实际是一个筛选的过程，首先将信号中频率最高的成份筛选出来，而后从原信号中将该成份去除，再从新的信号中选出频率最高的成份，依此类推，直到信号不可分解为止。这种过程可以看作为一系列的滤波器族，从算法的以上描述可知，信号可以通过有限的筛选步骤完成分解。每次筛选时，新信号的最大、最小值的个数都在减少。为了减少提取IMF的筛选步骤，定了了参数SD：

[0047]

[0048] 当SD小于某一常数时停止筛选，一般SD的值在0.2至0.3之间。另外在筛选过程中，由于该算法采用的是三次样条插值，所以当信号的极大值或者极小值的个数小于2时，停止筛选。这些有限的IMF经过希尔伯特变换后产生了以时间为参数的具有实际意义的瞬时能量和瞬时频率。这样可以在时域和频域同时对信号进行有效的分析，这种特性是以往以傅立叶变换为基础的信号分析方法所不具有的。

[0049] 由于经验模态方法是依据数据本身的时域信息进行的时域分解，得到的IMF通常个数是有限的和平稳的，而且是具有实际意义的窄带信号，基于这些IMF分量进行的Hilbert变换其结果反映了真实的物理信息，而且基于经验模态的Hilbert变换得到的每个IMF的振幅和频率是随时间变化的，消除了经典频谱分析法中为反映非线性、非平稳过程而引入的多余无物理意义的简谐波。因此，其Hilbert谱也能够准确反映出信号能量、频率在空间或时间尺度上的分布。经验模态方法是基于信号的局部特征时间尺度实现分解的。与小波分析法相比，经验模态-Hilbert方法具有小波分析的全部优点，并克服了小波变换的非自适应性，因此这种基于经验模态的Hilbert频谱分析方法，在非线性和非平稳过程的分析中具有很高的应用价值。

[0050] 从本质上说，经验模态是对数据在进行Hilbert变换之前作的一个预处理。通过经验模态，数据被分解为若干固有模态函数(IMF)的集合，每一个IMF刻画了信号的一个简单振荡模式。从表达形式上来看，一个IMF类似于Fourier分解表达式中的一个简谐振动，但是，它比简谐振动更加一般化。虽然经验模态分解方法的提出，为非线性非平稳算法的分析提供了有力的工具，但这种新方法还处在发展阶段，在实际应用中遇到很多问题，目前已有大量学者致力于该方法的研究和应用。

[0051] 3本发明的基于希尔伯特变换的时域能量特征分析方法

[0052] 首先考虑传统的基于频带功率谱密度定义的事件相关同步/去同步特征方法，该方法的基础是对脑电信号进行短时傅立叶变换，以获取指定频带范围下功率谱变化规律来分析事件相关同步/去同步特征。假定待分析脑电信号为f(t)，由于采集的脑电信号是离散时间序列，并不具备对称性，因此经过短时傅立叶变换后，可以获得其频域表示形式：

[0053] F(ω，t)＝a(ω，t)+jb(ω，t)＝A(ω，t)ejφ(ω，t) (3)

[0054] 式中F(ω，t)为f(t)的短时傅立叶变换结果。

[0055] 那么其对应的功率谱密度为：

[0056]

[0057] 对于本发明的方法获得的特征振荡模式，其本身处于一定的频带范围下，比如alpha节律下的特征振荡模式必然局限于8-14Hz，而后在该频带范围内表现出特定的功率谱分布。也希望获得类似于(3)的表达式，以获取sα(t)的能量变化特征。

[0058] 对于一个单组份实信号，Hilbert变换提供了一种可以获得类似于式(3)的解析信号表达式的方法。假设单组份信号为s(t)，s(t)的解析信号定义为：

[0059]

[0060] 解析信号的虚部为信号s(t)的Hilbert变换，计算如下：

[0061]

[0062] Hilbert变换器的频率响应为：

[0063] H(jω)＝-jsgn(ω) (7)

[0064] 可以获取(7)式在频域的表达式：

[0065]

[0066] Hilbert变换获取的是信号的瞬时值，B(t)为信号的能量变化趋势，ψ(t)为信号的瞬时相位，对于多组份复杂信号，瞬时相位没有意义，因此Hilbert变换仅适用于单组份信号。对于本发明获取的alpha节律下的特征振荡模式，根据EMD方法的定义，完全满足单组份信号的要求，因此可以采用(5)式获取特征振荡成份的解析信号。

[0067] 由此，对于特征振荡模式sα(t)，其能量变化特征可以用Bα(t)来反映，Bα(t)又被称为sα(t)的包络，即：

[0068]

[0069] 式中 为sα(t)的解析信号。Bα(t)即反映了alpha节律下特征振荡模式的能量变化特征。

[0070] 4能量特征提取算法

[0071] 本发明的复合下肢想象动作脑电的能量特征提取算法流程描述如下：

[0072] ①利用脑电导联电极采集复合下肢想象动作脑电；

[0073] ②对所采集的脑电信号进行滤波消噪和归一化处理；

[0074] ③复合下肢想象动作脑电的振荡模式分解：通过经验模态分解方法，将复合下肢想象动作脑电分解为频率从高到低的各个固有振荡模式分量。

[0075] ④特征振荡模式识别：通过功率谱密度分析方法分析各固有振荡模式的频谱分布特征，确定反映脑电特征节律的特征振荡模式。

[0076] ⑤特征振荡模式的去同步化特征获取：通过对特征振荡模式进行希尔伯特变换，获取信号包络，从而得到反映其能量变化的去同步化特征。

[0077] ⑥模式识别：将特征振荡模式的去同步化特征作为分类器输入，实现复合下肢想象动作的模式识别。

[0078] 特征振荡模式的能量变化特征的获取过程可以用图1表示。

[0079] 5实施例

[0080] 本发明采用奥地利EMSPHOENIX公司生产的128导数字脑电记录仪来采集脑电数据的设备。受试者在一个电磁屏蔽良好、隔音良好的房间内进行实验，房间内的背景噪声约2
为31dB，背景光照为2cd/m。受试者以感觉舒适但不影响数据采集的姿势坐在扶手椅中。
距离实验对象1米左右的正前方是一台19英寸显示器，用于显示实验对象进行复合下肢动作想象的提示符。每个复合下肢想象动作子实验持续8秒。第一时段为放松阶段，屏幕为无显示的黑屏状态，持续时间2秒。第二时段为准备期，此时段屏幕正中央显示一个十字提示符，提示受试者做好准备，该时段持续时间1秒。第三阶段为想象动作期，该时段持续时间
4秒，计算机屏幕上随机显示向左或者向右的箭头方位提示符，要求受试者进行相应的左腿和左手协同动作想象或者右手和右腿协同动作想象。第四阶段为恢复期，时段内显示器保持为无显示的黑屏状态。实验过程中要求受试者保持放松状态，不允许有任何的实际动作，并且为避免受试者由于视觉刺激引起的脑电波动，显示器以黑屏灰字的方式显示提示符。
整个实验方案要求每个受试者完成3组实验，每一组实验(run)包含2类复合下肢想象动作各30次子实验(trail)。在每两组实验之间，留有足够长的休息时间用于受试者进行疲劳恢复。如图2所示。

[0081] 电极的放置采用国际10/20系统标准，如图3所示，同时记录体四肢动作功能区附近的41个导联的EEG数据。这41导脑电导联包含通常采用的标准19导脑电导联，其余22导联是根据本发明中脑电采集目的，覆盖人体四肢动作的功能映射区，以此获取复合肢体动作想象过程中，更为精细的脑电特征。

[0082] 电极采用Ag/AgCI电极，并以左耳垂(A1)作为参考电平，右耳垂(A2)作为参考地，脑电采样频率为256Hz，滤波通带为0.5～35Hz。电极阻抗小于5000欧姆。

[0083] 数据采集完毕后，为提高后期模式识别的准确性，我们通过多个导联数据的线性组合进行空间滤波，以提高脑电信号的信噪比。首先采用了Hjorth提出的共参考平均方法去除脑电信号的参考电平，即从原始信号中减去所有电极信号的均值，其计算方法如下：

[0084]

[0085] 其中n是电极总数，ViER是电极上采集到的原始信号，通过共参考平均计算，在大部分电极中共有的空间低频成份将被取出，因此共参考平均方法的作用为高通空间滤波，可以突出在空间分布上高度集中的脑电成份。而后将所得信号中仍然存在较多肌电干扰或眼电干扰的子实验过程(trail)去除。

[0086] 图4给出了受试者在进行左手左腿协同动作想象时，以导联C3和C4为例，单次实验记录的脑电信号。

[0087] 对图4所示的信号，进行经验模态分解，结果如图5所示。图5(a)是对图4所示的C3导联采取经验模态分解后获得的7个固有振荡模式(IMF)及一个残留分量(res)。而图5(a)是对图4所示的C4联采取经验模态分解后获得的7个固有振荡模式(IMF)及一个残留分量(res)。

[0088] 为了获取脑电特征节律alpha节律的特征振荡模式，对图5所得的固有振荡模式进行功率谱密度分析。图6和图7分别给出了C3与C4导联脑电信号的前4个固有振荡模式(IMF)的功率谱密度曲线。从图中可以看出，C3，C4导联位置上，alpha节律的特征振荡模式均存在于第二固有振荡模式imf2上，并且C3导联位置上特征振荡模式的功率谱密度峰值要高于C4导联。这说明alpha节律的能量在C4导联受到的削弱要强于C3导联，即当进行复合左腿左手协同动作想象时，C4导联位置的alpha节律去同步化现象比C3导联位置的现象更为显著。

[0089] 采样上述同样的分析步骤，对复合右手右腿协同动作想象时的脑电信号进行经验模态分解。为了获取特征振荡模式，同样采用功率谱密度分析方法。图8给出了右手右腿协同动作想象时在C3，C4导联位置的特征振荡模式的功率谱密度曲线。显然，alpha节律的能量在C3导联受到更强的削弱，即当进行复合右腿右手协同动作想象时，C3导联位置的alpha节律去同步化现象比C4导联位置的现象更为显著。

[0090] 为了获取alpha节律去同步化现象的时变特征，对特征振荡模式进行希尔伯特变换，图9(a)和图9(b)分别给出了2类复合下肢动作想象脑电的alpha节律特征振荡模式的希尔伯特包络。从图9(a)可以看出，在左手左腿协同动作想象时，C3，C4导联都存在着alpha节律的去同步化现象，但是在导联C4位置更为显著，这与图6，图7的结论是一致的。相对的，在右手右腿协同动作想象时，C3导联位置的alpha节律去同步化现象则更显著。

[0091] 本发明根据Pfurtsehelle提出的事件相关同步化/去同步化定义，对上述获得的alpha节律的去同步化现象作进一步量化分析，给出基于特征振荡模式的事件相关去同步化系数的定义式：

[0092]

[0093] 其中REP为复合下肢动作想象开始后t1时间内的信号平均幅度，REF为复合下肢动作想象开始前t2时间内的信号平均幅度。这里，取t1为动作想象期的4-7秒，t2为动作想象前的0-2秒。由于获得的去同步化现象是通过幅度变化来反映能量变化的趋势(如图9)，因此，在采用本发明获得的特征振荡模式的希尔伯特包络进行去同步化量化系数计算时，需要将幅度信号转换为能量信号，即对信号幅度进行平方，而后根据式(10)进行计算。

[0094] 根据去同步化系数的定义，分别对各受试者在两类复合下肢动作想象的去同步化量化系数进行计算，而后进行T统计，统计分析结果如图10所示。结果表明，在进行左手左脚协同动作想象时，C4导联的去同步化现象更为显著，平均去同步化系数为55.6％，而在C3导联，平均去同步化系数仅有10.4％；在进行右手右脚协同动作想象时，C3导联显示出更为显著的去同步化现象，其平均去同步化系数达到了71.6％，而在对侧的C4导联，平均去同步化系数为13.2％。由于参与实验的受试者都是右利手和右利脚，所以在进行右手右脚协同动作想象时，引起的去同步化现象更强(平均71.6％，相对左手左脚同时动作想象时的55.6％)。

[0095] 按照上述步骤，对左右初级运动区和辅助功能区位置的脑电导联(包括F3P、FzP、F4P、C5A、C3A、C1A、C2A、C4A、C6A、C5、C3、C1、Cz、C2、C4、C6、C5P、C3P、C1P、C2P、C4P、C6P、P3A、PzA、P4A)分别进行特征振荡模式获取，而后计算特征振荡模式的事件相关去同步化系数，按照导联位置顺序构成去同步化系数向量，经过拼接后构成25维特征向量：

[0096] f＝[λF3P，λFzP，…，λP4A] (12)

[0097] 式中λF3P，λFzP，…，λP4A依次为导联F3P、FzP……P4A位置的事件相关去同步化系数。

[0098] 通过上述步骤，获取特征振荡模式的能量变化特征后，进一步考察该方法的可行性和有效性。将基于经验模态分解方法获取的去同步化特征与传统的通过基于短时傅里叶变换的时频图谱分析方法获取的去同步化特征分别用于分类模式识别。选择初级运动区与辅助运动区各导联脑电信号的去同步化系数作为分类器输入，并采用基于径向基核函数的支持向量基方法作为分类器。

[0099] 复合下肢想象动作电位的两分类分类器可以表示为：

[0100] y＝sgn(W·f+b) (13)

[0101] W和b是支持向量机的参数，通过训练样本得到。y是两分类的结果标记(±1)，+1表示识别出左手左腿同时动作想象，-1表示识别出右手右腿同时动作想象。f为(11)式对应的25维特征向量。

[0102] 图11给出了两种方法对应的分类结果对比，训练数据选择为总样本数目的10％到90％。从图11可以看出，基于经验模态分解方法对应的分类结果要明显优于传统的基于频带滤波的方法。另外一点需要说明，当训练样本只有总样本数的10％时，两种方法对应的分类结果均比较低，当训练样本增加到总样本数的90％以上时，识别准确率趋于稳定。

[0103] 本发明的复合下肢想象动作脑电的能量特征提取方法，具体步骤如下：

[0104] ①利用脑电导联电极采集复合下肢想象动作脑电，主要采集位于初级运动区及辅助功能区的脑电信号(包括F3P、FzP、F4P、C5A、C3A、C1A、C2A、C4A、C6A、C5、C3、C1、Cz、C2、C4、C6、C5P、C3P、C1P、C2P、C4P、C6P、P3A、PzA、P4A)；

[0105] ②对所采集的脑电信号进行滤波消噪和归一化处理；

[0106] ③针对经过步骤②处理后的脑电信号进行经验模态分解，将复合下肢想象动作脑电分解为频率从高到低的各个固有振荡模式分量，使原本非平稳的多组份脑电信号转化为平稳的单组份信号。

[0107] ④对步骤③获得的各导联脑电的固有振荡模式分量计算功率谱密度，根据其功率谱密度的频率分布范围，确定脑电特征节律-alhpa节律对应的固有振荡模式分量，即本发明所要获取的特征振荡模式。

[0108] ⑤ 定 义 特 征 振 荡 模 式 的 事 件 相 关 去 同 步 化 系 数，其中REP为复合下肢动作想象开始后t1时间内的
脑电信号平均幅度，REF为复合下肢动作想象开始前t2时间内的脑电信号平均幅度。这里，取t1为动作想象期的4-7秒，t2为动作想象前的0-2秒；对初级运动区及辅助功能区的脑电信号特征振荡模式计算事件相关去同步化系数，按照导联位置顺序构成事件相关去同步化系数向量f＝[λF3P，λFzP，…，λP4A]，这里λF3P，λFzP，·，λP4A依次为导联F3P、FzP……P4A位置的事件相关去同步化系数。

[0109] ⑥选择基于径向基核函数的支持向量基分类器，将步骤⑤获得的事件相关去同步化系数向量f＝[λF3P，λFzP，…，λP4A]作为分类器输入样本，分类器公式为y＝sgn(W·f+b)，W为权向量，b为阈值，通过训练样本得到通过训练样本得到，y是两分类的结果标记(±1)，+1表示识别出左手左腿同时动作想象，-1表示识别出右手右腿同时动作想象。通过分类器模式识别的结果，用于下一步转换为外部控制命令的输入参数。

[0110] 最终结果表明，通过经验模态分解方法结合功率谱密度分析及希尔伯特变换方法获取alpha节律能量变化特征是完全可行的，该方法的最高识别率为87.8％，相对于传统方法的82.3％有显著提高。该方法的优势在于充分考虑到脑电信号的非平稳性，并且基于动作想象脑电的事件相关同步/去同步特征是由群体神经元同步振荡产生的结论，因此在基于复合肢体动作的脑-机接口系统中将具有更为广阔的应用前景，这将为基于脑-机接口的下肢康复助行系统的发展提供有力支持。