基于HHT和AR模型的步行想象和空闲状态识别研究

程昭立，王晓琳，李红权，李鸿云，伏云发，3

1(昆明理工大学 信息工程与自动化学院，云南 昆明 650500)

2(昆明理工大学 脑认知与脑机智能融合创新团队，云南 昆明 650500)

3(云南省计算机技术应用重点实验室，云南 昆明 650500)

1 引 言

脑机接口(Brain Computer Interface，BCI)是一种变革性的人机交互，属于人机交互的子领域，旨在大脑与外设之间建立一条通路，而不受大脑外周神经和肌肉的控制[1，2].该技术可为运动障碍患者在特定情况下提供可选的、新的通信或控制方式，目前，基于运动想象的BCI是一种新兴的可用于运动障碍康复的技术[3].

脑卒中是一种脑组织损伤的疾病，造成该疾病的原因是脑部血管突然破裂或因血管阻塞导致血液不能流入大脑.近些年来我国每年新发脑卒中病例250万，其中50%～80%的脑卒中患者会遗留有不同程度的肢体运动功能障碍[4]，从而严重影响患者的生活自理能力和社会参与能力，给社会和家庭带来沉重的负担.已有研究表明，运动想象和实际肢体运动能激活相同的大脑皮层运动区[5]，由于人的大脑具有极强的可塑性，运动想象将有助于脑卒中偏瘫患者大脑受损神经的康复[6].

综上所述，本研究针对运动协调想象，设计了一种较新的实验范式，包括步行想象和空闲状态，前者为被试想象正常的步行动作，后者为被试保持清醒放松的状态.HHT是一种适用于非线性非平稳信号的有效时频分析方法[7-9]，考虑到脑电信号的非线性非平稳性，已有研究把HHT应用到MI-BCI中，且具有较好的效果，但是该方法尚没有在复合肢体参与的MI中验证其有效性.AR模型已用于MI-BCI中，该模型是一种时间序列分析方法，模型参数凝聚了重要信息，同时研究表明，AR模型的自回归参数对状态变化反映很敏感[10-12].因此，AR模型的自回归参数作为特征向量来表征被试状态相关的脑电变化具有较好的效果[13，14].本研究拟采用HHT提取平均瞬时能量特征，并与AR模型系数结合构造特征向量，最后采用支持向量机(Support Vector Machine，SVM)对步行想象和空闲状态两类任务进行分类，以验证其有效性.

2 材料及方法

2.1 被试

本研究招募了18名被试参与本实验，均为男性，本科以上学历，右利手，所有被试身体健康且无任何心理疾病，年龄在22～28岁之间，此前均无参与过相关实验.由于不同被试进行运动想象时激活脑区的强度有较大差异，故在实验前，依据运动觉-视觉想象问卷(The Kinesthetic and Visual Imagery Questionnaire，KVIQ)[15]对每名被试的运动想象能力进行评估.本实验获得昆明理工大学医学院医学伦理委员会的批准并遵照《赫尔辛基宣言》，每名被试在实验前都自愿签署了知情同意书.

2.2 实验范式

运动协调想象至少包含两种肢体参与协作，是复合肢体运动想象，其运动感觉的心理排演和传统单侧肢体简单运动想象的心理排演是不一样的.传统单侧肢体简单运动想象可能激活有限数目的脑区，而运动协调想象可能需要更多的脑区参与，协同相关的感知觉资源.其中，上肢和手的协调运动可以训练动作的精确性、反应速度以及动作的节奏性，而下肢的协调运动可以练习正确的步态和上下肢动作的配合、协调等.为探索脑卒中偏瘫患者运动障碍康复的方法，针对运动协调想象，本实验设计了一种新的步行想象范式，被试想象正常的步行动作，步行动作涉及人体四肢的协调配合，同时又设计了任务以识别被试的空闲状态.

对于步行想象，一次完整试验的时序如图1所示，在进行试验前，先告知被试试验任务为进行步行想象以及保持空闲状态，步行想象为在大脑中想象正常的行走，空闲状态为全身心清醒放松的状态.为了提高被试对两种任务的适应度，在第一个trial开始前，被试先观看正常步行动作视频，并认真体会肢体感觉，以强化步行动作在大脑中的预演，以及练习保持空闲状态，每名被试仅在第一个trial开始前进行两种任务的预演.在t=0s时，屏幕上呈现“固定十字”，提示被试准备开始实验，持续时间为3s；在t=3s时，屏幕上的固定十字消失，且呈现试验任务的提示，提示被试准备进行对应试验任务，持续时间为1s；在t=4s时，提示消失，屏幕为空屏，要求被试执行试验任务，持续时间为4s；t=8s时，任务结束，屏幕呈现星形提示，提示被试进行休息，持续时间为2s.休息结束，开始下一个trial.每名被试对两种任务分别进行80个trials，两种任务则为160个trials.

图1 实验范式

2.3 实验设置

本研究使用的实验设备是北京中科新拓仪器有限责任公司生产的NT9200系列EEG采集设备，电极帽为32通道，采样率设为1000Hz，带通滤波设为0.1-100Hz之间，同时采用50Hz陷波滤波器避免电源线干扰.电极帽的电极位置按照国际标准10-20导联系统设置[16]，参考电极为A1和A2，所有电极的阻抗都调至10kΩ以下.在本实验中，采集了FC3、FCz、FC4、C3、Cz、C4、CP3、CPz、CP4、O1、Oz、O2共12个通道的脑电数据进行研究.电极布局示意图如图2所示.

图2 电极布局示意图

2.4 数据分析

2.4.1 数据预处理

现有研究表明，运动想象能够激活相关大脑感觉皮层运动区的脑电活动，并且主要在μ(8-13Hz)和β(14-30Hz)节律段表现最为显著[17，18].脑电信号是一种很微弱的信号，幅值大概在10-100uV之间，很容易受到其他电极或者环境造成的干扰，因此需要对脑电信号进行预处理，从而提高信噪比[19，20].本文首先采用椭圆滤波器对原始脑电数据进行8-30Hz的数字带通滤波，通带衰减为0.5dB，阻带衰减为50dB.然后对滤波后的脑电信号进行基线漂移校正，以消除脑电信号相对于基线的偏离.最后采用独立成分分析(Independent Component Analysis，ICA)去除眼电、肌电等伪迹.

2.4.2 希尔伯特-黄变换

特征提取是脑机接口系统中的一个关键部分，通过特征提取可以将表征被试进行运动想象任务的大脑活动特征量提取出来.对于大脑活动中差异显著的特征量，可以降低模式识别分类器的错误分类率，从而有效提高特征分类正确率.

希尔伯特-黄变换(Hilbert-Huang Transform，HHT)由美国工程院院士黄锷于1998年提出[21]，HHT是一种适用于处理非线性、非平稳信号的自适应信号处理方法，它能够很好地描述信号中存在的局部振荡频率分量[22].HHT主要由两个部分组成：1)经验模态分解(Empirical Mode Decomposition，EMD)；2)Hilbert谱分析(Hilbert Spectral Analysis，HSA).

由于任意信号都可以看作是由一系列的本征振荡模态组成，则经验模态分解能够根据信号特征自适应地将任意复杂信号分解成一系列本征模函数(Intrinsic Mode Function，IMF)，且EMD具有自适应性、正交性、完备性以及对IMF分量的调制特性[23].经验模态分解需要满足两个条件[24]：1)初始信号中极值点的数量和过零点数量之差不超过1；2)初始信号上下包络关于时间轴局部对称.

EMD分解过程如下[25，26]：

1)求出原始数据x(t)中所有的极大值点和极小值点，并用Spline插值拟合形成原数据的上包络线，以及将所有的极小值点通过Spline插值拟合形成数据的下包络线，上包络线和下包络线的均值记为h1，将原始数据x(t)减去该平均包络h1，得到一个新的数据序列k1：

x(t)-h1=k1

(1)

如果k1满足IMF的两个条件，则k1为x(t)的第1个分量.

2)若k1不满足IMF的两个条件，则将k1作为原始数据，重复步骤1)，直到循环第i次，得到k1(i-1)-h1i=k1i，使得k1i满足IMF条件.记d1=k1i，则d1为原始数据x(t)的第1个IMF分量.

3)从x(t)中分离出d1，得到：

r1=x(t)-d1

(2)

将r1作为原数据重复步骤1)-步骤2)，求得x(t)的第2个IMF分量d2，重复循环n次，得到n个IMF分量.

(3)

若rn成为单调函数且不能再从中提取IMF分量时，则循环结束.由式(2)和式(3)得到：

(4)

因此，对于任何信号x(t)，都可以将其分解为n个IMF分量和一个残量rn之和.其中，n个IMF分量d1，d2，…，dn分别包含了数据信号从高到低不同频段的成分，且其都是平稳的.

经EMD分解后，再进行Hilbert谱分析，下式对每个IMF分量进行Hilbert谱变换：

(5)

解析信号为：

zj(t)=dj(t)+iyj(t)=aj(t)eiθj(t)

(6)

进而，可以分别得到瞬时幅值和瞬时相位为：

(7)

(8)

本研究采用HHT变换提取平均瞬时能量再结合AR模型系数构造特征向量.首先，对每一种任务下的每个脑电信号进行EMD分解，然后，对每一个IMF分量用Hilbert变换进行谱分析，计算出经HHT变换后幅值的平均瞬时能量，再结合AR模型系数构造特征向量，最后采用SVM对测试集分类.

2.4.3 AR模型

对于一般的随机信号，其AR模型可以描述如下：

其中，y(i)是信号的第i个采样值，φi是AR模型的第j个系数，n(i)是白噪声残差，p是AR模型的阶数.

本研究首先通过AIC准则对AR模型进行定阶，得到的AR模型最佳阶次为6.然后采用Burg算法提取6阶的AR模型系数AR1，…，AR6，由C3、Cz和C4这3个通道的AR模型系数组成18维特征向量：{C3AR1，…，C3AR6，CzAR1，…，CzAR6，C4AR1，…，C4AR6}.

最后把C3、Cz和C4这3个通道由AR模型系数构成的18维特征向量与alpha和beta频带能量构成的2维特征组合成20维的特征向量，这20维的特征向量由2.4.4节的支持向量机方法进行分类.整个特征提取过程如图3所示.

图3 特征提取过程

2.4.4 支持向量机

支持向量机(Support Vector Machine，SVM)[27]通过寻求结构风险最小化来提升算法的泛化能力，使得置信范围和经验风险达到最小，从而能够在样本量较少的情况下，也能获得较高的分类准确率.

SVM算法的关键在于[28]：1)当线性可分时，利用所求的支持向量来寻求一个最大间隔的超平面；2)当线性不可分时，通过核函数将在低维空间中线性不可分的特征向量映射到高维空间中使得线性可分.本实验采集的脑电信号正是样本量较少且在低维空间中线性不可分的，因此将SVM算法应用于本实验.

在本研究中，对C3、Cz、C4这3个通道的脑电数据进行分析，实验数据为80个训练样本和80个测试样本，以及训练样本和测试样本的标签.首先用训练样本和标签训练模型，然后用训练好的模型对测试样本进行标签预测，最后将预测标签与真实标签进行对比，从而得到分类正确率.

3 结 果

表1为12个电极通道每个频段的t检验值和步行想象的方差分析p值，方差分析因子有3个：第1个是试验任务，包括步行想象及空闲状态；第2个是脑电信号的频段范围(delta，theta，alpha，beta，gamma)；第3个是脑电电极通道(12个通道).显著性水平设为p<0.05.两种视觉想象任务的t检验结果如表1所示，表中还包括了方差分析得出的每个电极位置的p值和df值.从表中可以看出，最显著的差异出现在C3、Cz和C4的α频带(p<0.001).除此之外，还进行了t检验，以确定在哪个通道有显著性的差异，结果表明，位置C3、Cz、C4的p值小于0.01，这表明步行想象和空闲状态两种任务之间存在99%以上的显著差异.

表1 12个电极通道每个频段的t检验值和方差分析p值

本文计算各个IMF的方差贡献率与相关系数以进行IMF筛选，如表2所示.从表2可以看出，前4阶IMF的方差贡献率之和达到了96.47%，相关系数也很高.因此，本文对前4阶IMF采用Hibert变换并求取平均能量值作为能量特征.

表2 所有IMF分量的方差贡献率与相关系数

表3呈现了采用HHT和AR模型组合方法提取特征(试验任务开始后4s期间)并利用SVM分类器取的平均、最高和最低分类率.

表3 采用HHT和AR模型组合方法提取特征取得的分类正确率

表4呈现了分别采用CSP、AR、HHT和AR模型组合方法提取特征(试验任务开始后4s期间：4～8s)并利用SVM取得的平均、最高和最低分类率.

表4 分别采用CSP、AR、HHT和AR模型组合方法提取特征并利用SVM取得的平均、最高和最低分类率

图4为步行想象脑电信号经过EMD后得到的IMF，由图3可以看出，脑电信号分解成了8个IMF分量.

图5呈现了分别采用CSP提取特征、AR模型提取特征、HHT和AR模型结合提取特征所得的分类正确率随时间变化的曲线.图4(a)结果表明，测试集在0～1.5s间的分类曲线的趋势是向上的，且包含一个极大值点，在1.5～2.5s之间波动幅度比较大，在2.5～4s趋于平稳，并维持在80%附近.图4(b)结果表明，测试集在整个0～4s期间的分类曲线波动幅度较大，在0～1.5s振荡向上，包含一个极大值点，1.5～3s之间有较大幅度下降，包含一个极小值点，3～4s之间又有较大幅度上升，并在4s达到最大值.图4(c)结果表明，测试集在0～4间的分类曲线整体趋于平稳，0～1.2s之间趋于上升，之后便趋于平稳，并维持在83%附近.

图4 步行想象脑电信号EMD后各阶IMF

图5 分类正确率随时间变化的曲线

4 讨 论

迄今为止，大多数运动想象针对的是单侧肢体的研究，而传统的单侧肢体简单的运动想象不易被患者理解并且康复效果有限.与传统单侧肢体简单运动想象相比，运动协调想象至少包含两种肢体参与协作，是复合肢体运动想象，其运动感觉的心理排演和传统单侧肢体简单运动想象的心理排演是不一样的.因此，有别于传统单侧肢体运动想象研究，本研究针对运动协调想象，设计了一种新的实验范式，包括想象正常的步行以及保持清醒放松的状态，易于被试理解.

本论文在实验范式上选择步行想象和空闲状态作为试验任务对，科学假设是该两种任务诱发的脑电特征具有可分性.由表1的结果可以看出，两种任务诱发的脑电信号之间最显著的差异出现在C3、Cz和C4的α频带；由表2的结果可以看出，前四阶IMF分量的方差贡献率达到了96.47%；由表3的结果可以看出，采用HHT和AR模型相结合的方法提取特征能够较好地区分(平均分类正确率为80.51±4.06%)步行想象和空闲状态两种任务.因此，本研究设计的实验任务具有可分性.可以预见，如果能够找到更好的特征提取和分类方法，还有可能提高分类正确率，这是我们未来要开展的一个工作.

在本文的研究中，共空间模式(CSP)特征提取与自回归(AR)模型特征提取相比，最高分类正确率提高了5.32个百分点，平均分类正确率提高了10.29个百分点，这表明前者优于后者；希尔伯特-黄变换(HHT)与AR模型相结合的特征提取和CSP特征提取相比，最高分类正确率提高了4.12个百分点，平均分类正确率提高了6.93个百分点，这表明该组合方法优于CSP特征提取方法，更优于AR模型.比较图5的结果，HHT和AR模型相结合的方法比CSP具有更佳的分类正确率随时间变化的曲线(其分类曲线整体上更趋于稳定且呈上升状态).

表5呈现了运动想象研究中实验范式、特征提取方法、分类方法及分类正确率.本研究对步行想象和空闲状态的平均分类正确率比Tariq M等人[29]对左右膝关节伸展运动想象的平均分类正确率高10.51个百分点，我们推测可能是由于想象左右膝关节伸展具有一定的难度，导致诱发的脑电信号可分性较差；本研究的平均分类正确率比Hsu W C等人[30]对左右腿踩踏运动想象的平均分类正确率低5.74个百分点，比唐智川等人[31]对左手和脚运动想象的平均分类正确率低7.47个百分点，比任士鑫等人[32]对脚运动想象和空想的平均分类正确率低1.49个百分点，我们推测可能是由于实验范式有较大的差异，且特征提取和分类方法各不相同.

表5 运动想象研究中实验范式、特征提取方法、分类方法及分类正确率

本研究进行的是离线实验任务以及在健康被试中进行的实验，我们未来的研究工作是：1)在线验证和完善所提方法；2)进一步完善实验范式，可以考虑将步行想象分解为数个动作；3)将此实验应用到脑卒中患者，以期步行想象能应用到脑卒中患者运动功能障碍康复中.

5 结 论

有别于传统的单侧肢体简单运动想象，运动协调想象至少包含两种肢体参与协作，是复合肢体运动想象，其激活的脑区数目更多，且更易于被试理解并在大脑中预演.针对运动协调想象，同时为提高试验任务对诱发的脑电特征之间的可分性，本研究提出一种新的且易于被试执行的步行想象和空闲状态实验范式.本研究表明该试验任务对诱发的脑电信号之间最显著的差异出现在C3、Cz和C4这3个通道的α频带.通过不同的特征提取方法，并采用SVM对步行想象和空闲状态两种任务进行分类，发现HHT和AR模型的组合方法优于CSP和AR模型，且平均分类正确率达到80.51±4.06%，具有良好的可分性.本研究可望为基于步行想象的BCI用于脑卒中偏瘫患者运动障碍康复提供思路，同时也可望为构建新型的在线下肢运动协调想象BCI提供策略.