利用肌肉协同研究上肢简单和复杂运动的关联性

黄永欣,文斌,毛菁菁,姚瑞,尹孟奇,徐进

(1. 西安交通大学生命科学与技术学院健康与康复科学研究所,710049,西安;2. 西安交通大学生物医学信息工程教育部重点实验室,710049,西安;3. 西安交通大学第一附属医院康复科,710061,西安)

肢体运动是人体重要的生理功能之一,中枢神经系统(central nervous system,CNS)通过控制相关肌肉运动的协同,以完成日常生活中的各种运动,由于运动模式多样,使得运动控制问题变得极为复杂[1-2]。

在人体运动控制研究中,肌肉协同学说被广泛用于研究神经肌肉控制问题,即CNS将相应部位的多块肌肉以不同强度的协同作用结合在一起,完成运动控制的任务[3-4]。CNS通过肌肉协同机制实现对运动控制的简化,协同被认为是神经运动模块化控制的基本单元[5-6]。因此,研究运动过程中相关肌肉的协同作用对于神经肌肉控制模式研究具有十分重要的意义。

表面肌电信号(surface electromyography, sEMG)是经皮肤表面电极采集的生物电信号,其蕴含丰富的神经肌肉活动和运动意图等信息,同时具有提取方便、准确和无创的优点,已经成为分析肢体运动过程中肌肉协同的有利手段[3,7]。

从多通道sEMG信号中获取肌肉协同信息属于盲源分离的过程,常见的矩阵分解算法有主成分分析(principal components analysis, PCA)、独立成分分析(independent component analysis, ICA)和非负矩阵分解(non-negative matrix factorization, NMF)[6,8]等。因为NMF的结果具有非负约束和鲁棒性强等优点,使得基于NMF的肌肉协同分析更具优势[8]。

肌肉协同分析作为一种探索神经肌肉控制模式的有效方法,已经被广泛用于临床研究和机器人控制等领域。该方法可用于探索神经功能障碍患者的肢体运动控制模式,进而为康复治疗提供指导。Steele等对比脑瘫患者和健康人在步态活动中提取的肌肉协同模式,发现脑瘫患者在步态活动中的肌肉变化用较少的肌肉协同即可表示[9]。Allen等通过比较帕金森患者在短期康复治疗前后的步行与平衡的肌肉协同模式,发现了步行与平衡存在相似的协同模式与特有的协同模式,且康复治疗后相似的协同模式数量增加[10]。

肌肉协同分析可以提取运动意图,从而实现运动控制和识别。桂奇政等利用肌肉协同理论和支持向量回归构建模型,实现了从sEMG信号到关节运动的连续估计[11]。郑楠等利用肌肉协同和最小二乘法构建手势识别算法,从神经协同的角度提取运动意图,减少训练数据的同时提高用户无关手势识别的正确率[12]。

研究特定肢体运动下的肌肉协同模式和不同肢体运动间肌肉协同模式的相似性以及差异,将能够更好地了解、掌握肢体运动中的神经肌肉控制模式。Saito等研究健康人在不同速度和坡度的跑步机上跑步的肌肉协同差异,发现不同条件下的跑步均可由4个肌肉协同来解释并存在相似的肌肉协同模式[13]。

与下肢运动相比,上肢运动模式多样且复杂,运动控制问题更为复杂。上肢运动的肌肉协同分析开展较少,且已有研究大多是针对简单运动。谢平等在研究健康人的腕屈伸运动的肌肉协同中发现了与腕部运动相关肌群的协同关系、共享的基本肌肉协同模式以及同一肌肉协同模式中协同性较高的肌肉间耦合关系较强[14]。Tang等通过对健康人的肘关节和肩关节的3个相似运动的肌肉协同分析,发现了3个相似运动中存在共享的基本肌肉协同模式[15]。王洪安等在伸手运动的肌间神经耦合关系研究中发现,卒中患者的肌肉协同数下降、β频段内肌间耦合强度具有显著差异以及γ频段内大部分协同肌肉的耦合关系较强[16]。

此外,协同模式的关联性研究较少。Cheung等在卒中患者的上肢健侧运动与患侧运动的肌肉协同模式关联性研究中,通过构建协同模式的线性模型并利用非负最小二乘法求解模型系数,发现了肌肉协同模式的保存、合并和分离[17]。Hashiguchi等利用相同方法研究亚急性期卒中患者在康复过程中的步态活动,发现肌肉协同模式的合并可以作为亚急性期卒中患者运动协调异常的重要标志[18]。

迄今,复杂运动的肌肉协同模式研究较少,上肢简单运动和复杂运动的肌肉协同模式之间的关联性研究缺乏,给上肢运动的神经肌肉控制机制的全面掌握带来困难。

因此,本文分别设计了上肢简单和复杂运动范式,采集不同运动执行过程中的sEMG信号,通过肌肉协同分析探究简单运动和复杂运动执行过程中肌肉协同模式之间的关联性,从而为肢体运动的神经肌肉控制模式和肌-机接口技术研究提供参考。

1 实验设计和信号采集

1.1 实验设计

选择日常生活中常见的上肢动作设计简单和复杂运动范式。简单运动范式包括:手伸展(HE)和手屈曲(HF)、腕屈曲(WF)和伸展(WE)、肘屈曲(EF)和伸展(EE)、肩前屈(SF)和外展(SA);复杂运动包括两个,复杂运动Ⅰ(ComplexF)主要由肩前屈、肘屈伸和手屈曲组成,复杂运动Ⅱ(ComplexA)主要由肩外展、肘屈伸和手屈曲组成,如图1、图2所示。

(a)手伸展

(a)复杂运动Ⅰ

实验利用Matlab软件和Psychtoolbox工具箱呈现运动指示图片以规范数据采集流程,如图3所示。

首先在屏幕上出现一个“+”注视点,持续2 s;当运动指示图片出现后,根据图片提示做相应的动作并保持到图片结束,简单运动图片持续10 s,复杂运动图片持续7 s;间隔休息30 s。每个动作完成之后,必须恢复初始位置状态再进行下一个动作。每种动作会重复3次,然后换下一种动作,每种动作结束后,间隔30 s用于休息和切换采集姿态,再继续进行实验。

另外,采集尺侧腕屈肌(flexor carpi ulnaris,FCU)、指浅屈肌(superficial flexor digitorum,SFD)、桡侧腕伸肌(extensor carpi radialis,ECR)、指伸肌(extensor digitorum,ED)、肱二头肌(biceps brachii,BB)、肱三头肌(triceps brachii,TB)、三角肌前束(anterior bundle of deltoid,ABD)和三角肌中束(middle bundle of deltoid,MBD)8个与运动相关肌肉的最大自主收缩(maximal voluntary contraction, MVC),用于后续肌电信号的MVC归一化,MVC动作模式如图4所示。

(a)尺侧腕屈肌

1.2 表面肌电信号采集

使用美国Delsys公司的无线肌电信号采集系统进行8通道肌电信号采集,采样频率为1 259 Hz。传感器顺着肌肉纤维的方向,贴附于肌腹隆起处。同步采集实验过程中右侧上肢的8通道sEMG信号,肌电电极位置及相关肌肉如图5所示。

图5 上肢运动相关肌肉与肌电采集位置Fig.5 Related muscles of the movements of upper limbs and electrode position of sEMG signal

来自西安交通大学的20名健康青年参加了实验,实验参与者均无上肢肌肉和中枢神经等疾病,受试前24 h未进行剧烈运动,均为右利手,男女比例为11∶9,年龄为25.1±3.093(均值±标准差)。

2 基于非负矩阵分解的肌肉协同分析

2.1 信号预处理

为了更好地分析数据,首先,提取sEMG原始信号进行预处理并提取信号包络[19-21],具体步骤如下:①通过截止频率为20 Hz的4阶巴特沃斯高通滤波器,用于去除信号的运动伪迹和基线漂移;②进行全波整流,使得信号满足非负要求,符合生理上的肌肉激活度;③信号通过截止频率为5 Hz的4阶巴特沃斯低通滤波器去除高频噪声,完成信号包络的提取。

其次,为了避免肌电幅值个体差异对肌肉协同分析的影响,对sEMG信号进行MVC归一化处理[6,19]。MVC归一化反映了该位置肌肉的激活程度,且保留了各通道之间的比例关系[22-23]。

MVC归一化的具体步骤为:对被试者MVC动作时的sEMG信号做上述的提取包络处理,然后计算滑动窗内最大平均值,窗长设置为250 ms。

最后,利用Teager-Kaiser能量算子[24]自动检测运动起始点和终止点,提取sEMG信号片段,降采样到200 Hz,并将处理后的片段进行串联。

2.2 肌肉协同提取

1999年,Lee和Seung提出了NMF算法,该算法在非负的约束下可以实现非线性的降维[25]。因为NMF的分解结果具有非负约束性,其分解的结果更具实际意义[26]。因此,本文选择NMF算法进行肌肉协同分析。

对于预处理后的sEMG信号X,肌肉协同提取算法需要找到r个肌肉协同及其对应的线性加权系数,使之能较好地重构信号X,即

Xm×n=Wm×rHr×n+Em×n

(1)

式中:W是肌肉协同矩阵;H是肌肉协同的加权系数矩阵;E是误差矩阵;r为肌肉协同数;m是肌电导联数;n为肌电信号长度。

NMF是一个最值优化问题,在肌肉协同研究中的求解算法主要有基于欧式距离的乘性迭代算法、基于交替非负最小二乘法的算法和基于梯度下降的加性算法[27]等。基于欧氏距离的乘性迭代算法具有较优的分解效果,但收敛性不稳定,可能陷入局部最优;基于梯度下降的加性算法虽然简单但收敛速度慢,迭代过程易受迭代参数影响;基于交替非负最小二乘法的算法即固定一个矩阵对另一个矩阵进行优化,不断交替进行,可收敛至稳定点[28]。

因此,为了避免陷入局部最优,本文使用了由两步骤组成的NMF算法,具体如下。

(1)随机选择W和H,并重复使用20次基于欧式距离的乘性迭代算法运行NMF(假设噪声符合高斯分布),选择初步收敛最好的结果作为矩阵分解结果。

(2)使用初步收敛结果最好的W和H作为第二轮NMF分解的初始值,使用基于交替非负最小二乘法的算法运行NMF。

2.3 肌肉协同数确定

矩阵X经NMF后得到W和H矩阵,但NMF算法需要确定肌肉协同数r。通常计算原始矩阵X和重构矩阵X′=WH的变异性占比(RVAF)[9,16,29-30]来确定协同数r,计算公式为

(2)

RVAF越大,表示重构矩阵X′与原始矩阵X越接近。计算不同协同数r下的所有被试者平均的RVAF,具体结果如图6所示,图6中误差线为标准差,Th为阈值。

图6 不同动作下r与RVAF之间的关系Fig.6 The relationship between r and RVAF under different movements

以往的研究确定肌肉协同数时,RVAF阈值常设在0.80～0.95之间[21]。因此,为了更好地重构信号并保留原始信息,本文设定当满足RVAF的均值大于等于95%,且随着协同数r的增加RVAF的增量小于等于1%时[26],所分解的肌肉协同数r作为最佳分解数。满足上述条件,得到简单、复杂运动的最佳分解数和模式,见表1、表2,以手伸展分解出的协同模式HE-1为例进行解释,其中HE表示手伸展,1为编号,其他类同。

表1 简单运动的最佳分解数和协同模式

表2 复杂运动的最佳分解数和协同模式

2.4 协同模式排序

对同一动作下不同被试者提取到了类似的协同矩阵W,但提取的协同模式Wi,j以随机顺序出现,因此需要对协同模式根据相似度进行排序,具体步骤如下。

(1)选取被试者k的协同模式作为初始模板,其顺序为Wk,1,Wk,2,…,Wk,r。

(2)计算每个被试者提取的协同模式与初始模板的余弦相似度[30],计算公式为

(3)

式中:Wi,j为第i个被试者提取的第j个协同模式;Wk,r为初始模板的第r个协同模式;cosθi,k,j,r为Wi,j与Wk,r的余弦相似度。

(4)从相似度矩阵Ri,k找到最大值maxj,r,即第i个被试者提取的协同模式Wi,j的顺序位置为初始模板中的协同模式Wk,r的位置,并将相似度矩阵对应的第j行和第r列置为0,重复步骤(4),依次确定其他协同模式顺序。

(5)计算所有被试者排序后的协同模式与该初始模板的平均相似度,选择其他被试者作为初始模板并重复步骤(2)～(5)。

(6)找出平均相似度最高的初始模板作为最终排序模板。

2.5 关联性分析

为了研究简单、复杂运动肌肉协同模式的关联性,首先构建线性组合模型如下式

(4)

其次,利用非负最小二乘法[17,31]求解系数mi,k,并对求解后的系数mi,k进行标准化,并选取0.2作为贡献度阈值[17,31]。

最后,选取大于0.2阈值的非负线性系数及其对应的简单运动协同模式构建重构运动模式,并通过余弦相似度量化重构运动模式与复杂运动协同模式之间的相似性,以探讨简单运动和复杂运动的关联性。

3 结果及讨论

3.1 简单运动与复杂运动的协同模式

基于排序后的协同模式,将所有被试者在简单运动和复杂运动协同模式下相关肌肉的权值系数进行平均,并对简单运动和复杂运动的协同模式进行分析。

3.1.1 简单运动分析

(1)手部动作。手部动作协同模式下相关肌肉的权值系数(均值±标准差),如图7所示。手部动作分解的协同主要由FCU、SFD、ECR、ED 4块肌肉构成。手伸展动作主要分解为一个协同,模式HE-1以ED为主,ED为前臂后群肌,其可以使得后四指掌指关节伸展,同时HE-1中显示了FCU、SFD、ECR和ED 4块肌肉的协同配合。手屈曲动作分解为两个协同,HF-1以SFD为主,SFD为前臂前群肌,负责第2～5指的掌指关节和近侧指骨间关节的屈曲;HF-2中FCU激活度最高,表明FCU对手屈曲也有一定作用。简单手部动作的协同模式中均显示了4块肌肉的协同关系,表明主要屈肌和主要伸肌会相互配合。

(a)手伸展

(2)腕部动作。腕部动作协同模式下相关肌肉的权值系数(均值±标准差),如图8所示。腕部动作分解的协同也主要涉及FCU、SFD、ECR、ED 4块肌肉。腕伸展动作主要分解为两个协同,模式WE-1以ED为主,体现该模式与手伸展运动有关,表明ED除了参与四指掌指关节的伸展,还参与腕关节的伸展运动;模式WE-2以ECR和ED为主,ECR为前臂后群肌,其主要功能就是伸展腕关节。腕屈曲动作分解为两个协同,模式WF-1主要由SFD、FCU和ED贡献,尤其是SFD,体现该模式与手屈曲运动有关,表明SFD虽然主要负责指关节屈曲,但是也会支持腕部屈曲;模式WF-2以FCU为主,FCU为前臂前群肌,是腕关节屈曲的主要肌肉之一。

(a)腕伸展

(3)肘部动作。肘部动作协同模式下相关肌肉的权值系数(均值±标准差),如图9所示。

(a)肘伸展

肘伸展动作分解为两个协同,模式EE-1以TB为主导,TB属于上臂后侧肌群,主要负责伸肘;模式EE-2中ED存在较大的肌肉激活度,表明ED对肘伸展也有一定作用。在肘屈曲动作中,模式EF-1中以ECR、FCU、ED和BB为主,表明4块肌肉协同配合完成曲肘;模式EF-2以BB为主,BB位于上臂前侧,当BB收缩时会使肘关节屈曲。

(4)肩部动作。肩部动作协同模式下相关肌肉的权值系数(均值±标准差),如图10所示。在肩前屈动作的协同模式SF-2中,体现了ABD在肩前屈运动中起到的重要作用。与此类似,在模式SA-2中,表明MBD在肩外展动作起到了主要作用,同时TB和ABD协同配合进行肩部外展。此外,在简单肩部动作的协同模式SF-1和SA-1中表现出MBD、ABD和TB等肌肉之间的协同关系。

(a)肩前屈

3.1.2 复杂运动分析

复杂运动协同模式下相关肌肉的权值系数(均值±标准差),如图11所示。复杂运动Ⅰ分解出3个协同模式,模式ComplexF-1中FCU、SFD、ECR和ED有较大贡献,体现了4块肌肉之间的协同配合,反映了该范式主要与手部动作或腕部动作有关;模式ComplexF-2体现了FCU、ECR、ED和ABD的主要作用,表明与手部动作、腕部动作和肩部动作有关;模式ComplexF-3中ABD和MBD贡献较大,尤其是ABD,说明该模式与简单运动的肩前屈有关。

复杂运动Ⅱ分解出3个协同模式,模式ComplexA-1中FCU、SFD、ECR和ED有较大贡献,尤其是FCU,表明与手屈曲、腕屈曲动作有关;模式ComplexA-2和模式ComplexA-3体现了MBD、ED、FCU和SFD等肌肉的主要作用,尤其是MBD,说明这两个模式与手部动作、腕部动作和肩外展动作有关。

3.2 关联性研究

通过线性拟合,并对系数进行标准化,线性组合的关系如图12所示。以复杂运动模式ComplexF-1为例,其可以由非零系数(0.19、0.06、0.51、0.05、0.04、0.04和0.11)及对应的简单运动模式(HF-2、WE-2、WF-1、EE-1、EF-2、SF-2和SA-2)构成。同理,可以得到其他复杂运动协同模式的构成。之后,选取大于0.2的系数及其对应的简单运动协同模式,用于后续构建重构运动模式。

红色方框中数据为大于0.2阈值的非负线性系数图12 肌肉协同模式的线性组合关系Fig.12 Linear combination of muscle synergy patterns

在以往协同模式的相似性研究中,余弦相似度大于0.75定性为相似,大于0.9定性为非常相似[17,30]。

通过复杂运动协同模式与重构运动模式之间的相关性分析,发现了协同模式的保留和合并,如图13所示。

(a)复杂运动Ⅰ(Fi,i=1,…,5,Fi为大于0.2阈值的非负线性系数)

在复杂运动Ⅰ中,模式ComplexF-1与对应的重构运动模式的余弦相似度为0.899,具有较高的相似度,同时,该重构运动模式仅由WF-1与非线性系数构成。因此,模式ComplexF-1可以由WF-1构成,体现了协同模式的保留;模式ComplexF-2与对应的重构运动模式非常相似,两者的余弦相似度为0.984,因此模式ComplexF-2可以通过组合HF-2、WE-1和SF-2来重构,体现了复杂运动协同模式为简单运动协同模式的合并;同样,模式ComplexF-3与对应的重构运动模式的相似度为0.957,说明模式ComplexF-3可以由SF-2来重构,体现了复杂运动协同模式为简单运动协同模式的保留。

重构运动模式与复杂运动Ⅱ的肌肉协同模式也具有很高的相似度(在0.909～0.995之间)。在复杂运动Ⅱ中,模式 ComplexA-1主要由HF-2和WF-1组成,表明与手部屈曲和腕部屈曲运动有关;模式ComplexA-2由HF-2、WE-1和SA-2组成,体现了协同模式的合并;模式 ComplexA-3主要由HF-2、WF-1和SA-2组成,表明与手部屈曲、腕部屈曲和肩部外展运动的组合有关,体现了复杂运动协同模式为简单运动协同模式的线性组合。

4 结 论

本文基于表面肌电信号,通过采用非负张量分解和余弦相似度分析等方法,分别对简单运动和复杂运动执行过程中的肌肉协同状况进行比较分析。研究发现复杂运动中存在明显的简单运动肌肉协同模式的保留和合并过程。从肌肉协同的角度证明了人体复杂运动模式是从简单运动提取的基本肌肉协同模式的基础上线性组合而来,复杂运动与简单运动之间还存在共享的基本肌肉协同模式。本文研究结果有助于进一步增强对神经肌肉运动控制机制的认识,为肌-机控制技术相关领域研究提供理论基础。