第15卷第2期 智能系统学报 Vol.15 No.2 2020年3月 CAAI Transactions on Intelligent Systems Mar.2020 D0:10.11992/tis.202001013 网络出版地址：http:/kns.cnki.net/kcms/detail/23.1538.TP.20200319.1521.002.html 利用肌电信号求解关节力矩的研究及应用综述 姜峰，尹逊锋2，衣淳植，杨炽夫2 (1.哈尔滨工业大学计算机科学与技术学院，黑龙江哈尔滨150001；2.哈尔滨工业大学机电工程学院，黑龙 江哈尔滨150001) 摘要：表面肌电信号(surface electromyography,sEMG)是人体的易于检测的神经信号，其富含大量人体运动信 息。利用肌电信号作为输入信号，结合相关生物学模型分析肌电信号同肌肉力和对应关节力矩之间的关系，对 于深入理解分析人体动力学具有重要意义。本文详细归纳总结了利用肌电信号求解人体关节力矩方法的研究 成果，同时介绍神经肌肉骨骼模型的计算及优化过程，给出部分模型生理参数为之后的研究提供参照，并给出 现阶段该方法在人体关节力矩求解中的应用。通过分析该求解过程中所面临的一些问题，总结出该方法的发 展展望，为之后的研究提供参考。 关键词：表面肌电信号：神经激活：信号处理：肌肉激活：肌肉骨骼模型：关节力矩：关节动力学：参数辨识 中图分类号：TP391.4文献标志码：A文章编号：1673-4785(2020)02-0193-11 中文引用格式：姜峰，尹逊锋，衣淳植，等.利用肌电信号求解关节力矩的研究及应用综述J儿.智能系统学报，2020,15(2)： 193-203. 英文引用格式：JIANG Feng.,YIN Xunfeng,YI Chunzhi,,etal.A review of the research and application of calculating joint torque by electromyography signalsJ.CAAI transactions on intelligent systems,2020,15(2):193-203. A review of the research and application of calculating joint torque by electromyography signals JIANG Feng,YIN Xunfeng',YI Chunzhi,YANG Chifu2 (1.College of Computer Science and Technology,Harbin Institute of Technology,Harbin 150001,China;2.School of Mechatronics Engineering.Harbin Institute of Technology,Harbin 150001,China) Abstract:Surface electromyography is a neuronal signal of the human body that is easily detected.It extensively provides information about human motion.For a deeper understanding of human body dynamics,the use of elec- tromyographic signals and biological models to examine the relationship between myoelectric signals and muscle forces or corresponding joint torques,is of great importance.This paper summarizes:the research results of solving human joint torque using electromyography signals;introduces the process of measuring and optimizing the neuromusculoskeletal model;gives some model physiological parameters to provide reference for future research;and provides application of the method in solving human joint torque in the current stage.Some of the problems encountered in the solution process are then evaluated,and the development prospect of the method is summarized,providing a reference for future research. Keywords:surface electromyography;neural activation;signal processing;muscle activation;musculoskeletal model; joint moment;joint dynamics,parameter identification 由于人体运动学信息在教学、医疗)、运动 两种：l)采取人体表面肌电信号(surface ectrom- 康复B、人体外骨骼5刀等领域的广泛应用，人体 yography)结和人体肌肉骨骼模型求解相应的关 关节力矩的求解方法一直受到学者的广泛关注。 节力矩⑧；2)运用惯性测量元件测量人体运动学 目前对于关节力矩的求解方法根据原理主要分为 信息，通过构建人体逆动力学模型来解算关节力 收稿日期：2020-01-08.网络出版日期：2020-03-20 矩。在人体逆动力学模型方面，哈尔滨工业大学 基金项目：国家重点研发计划项目(2018YFC0806800,2018YFC 0832105). 的郭伟等将单腿支撑相的人体建立为刚体模 通信作者：尹逊锋.E-mail:mr yinxf@hit.edu.cn 型，通过光学动捕系统测量人体运动学信息，对
DOI: 10.11992/tis.202001013 网络出版地址: http://kns.cnki.net/kcms/detail/23.1538.TP.20200319.1521.002.html 利用肌电信号求解关节力矩的研究及应用综述 姜峰1 ，尹逊锋2 ，衣淳植2 ，杨炽夫2 （1. 哈尔滨工业大学 计算机科学与技术学院，黑龙江 哈尔滨 150001; 2. 哈尔滨工业大学 机电工程学院，黑龙 江 哈尔滨 150001） 摘 要：表面肌电信号 (surface electromyography, sEMG) 是人体的易于检测的神经信号，其富含大量人体运动信 息。利用肌电信号作为输入信号，结合相关生物学模型分析肌电信号同肌肉力和对应关节力矩之间的关系，对 于深入理解分析人体动力学具有重要意义。本文详细归纳总结了利用肌电信号求解人体关节力矩方法的研究 成果，同时介绍神经肌肉骨骼模型的计算及优化过程，给出部分模型生理参数为之后的研究提供参照，并给出 现阶段该方法在人体关节力矩求解中的应用。通过分析该求解过程中所面临的一些问题，总结出该方法的发 展展望，为之后的研究提供参考。 关键词：表面肌电信号；神经激活；信号处理；肌肉激活；肌肉骨骼模型；关节力矩；关节动力学；参数辨识 中图分类号：TP391.4 文献标志码：A 文章编号：1673−4785(2020)02−0193−11 中文引用格式：姜峰, 尹逊锋, 衣淳植, 等. 利用肌电信号求解关节力矩的研究及应用综述 [J]. 智能系统学报, 2020, 15(2): 193–203. 英文引用格式：JIANG Feng, YIN Xunfeng, YI Chunzhi, et al. A review of the research and application of calculating joint torque by electromyography signals[J]. CAAI transactions on intelligent systems, 2020, 15(2): 193–203. A review of the research and application of calculating joint torque by electromyography signals JIANG Feng1 ，YIN Xunfeng2 ，YI Chunzhi2 ，YANG Chifu2 (1. College of Computer Science and Technology, Harbin Institute of Technology, Harbin 150001, China; 2. School of Mechatronics Engineering, Harbin Institute of Technology, Harbin 150001, China) Abstract: Surface electromyography is a neuronal signal of the human body that is easily detected. It extensively provides information about human motion. For a deeper understanding of human body dynamics, the use of electromyographic signals and biological models to examine the relationship between myoelectric signals and muscle forces or corresponding joint torques, is of great importance. This paper summarizes:the research results of solving human joint torque using electromyography signals; introduces the process of measuring and optimizing the neuromusculoskeletal model;gives some model physiological parameters to provide reference for future research; and provides application of the method in solving human joint torque in the current stage. Some of the problems encountered in the solution process are then evaluated, and the development prospect of the method is summarized, providing a reference for future research. Keywords: surface electromyography; neural activation; signal processing; muscle activation; musculoskeletal model; joint moment; joint dynamics; parameter identification 由于人体运动学信息在教学[1] 、医疗[2] 、运动 康复[3-4] 、人体外骨骼[5-7] 等领域的广泛应用，人体 关节力矩的求解方法一直受到学者的广泛关注。 目前对于关节力矩的求解方法根据原理主要分为 两种：1) 采取人体表面肌电信号 (surface ectromyography) 结和人体肌肉骨骼模型求解相应的关 节力矩[8-9] ；2) 运用惯性测量元件测量人体运动学 信息，通过构建人体逆动力学模型来解算关节力 矩。在人体逆动力学模型方面，哈尔滨工业大学 的郭伟等[10] 将单腿支撑相的人体建立为刚体模 型，通过光学动捕系统测量人体运动学信息，对 收稿日期：2020−01−08. 网络出版日期：2020−03−20. 基金项目：国家重点研发计划项目 (2018YFC0806800, 2018YFC 0832105). 通信作者：尹逊锋. E-mail：mr_yinxf@hit.edu.cn. 第 15 卷第 2 期 智 能 系 统 学 报 Vol.15 No.2 2020 年 3 月 CAAI Transactions on Intelligent Systems Mar. 2020
·194· 智能系统学报 第15卷 人体下肢关节力矩进行解算；韩国科学技术院结 2神经激活求解 合人体上肢质心的运动规律，建立了弹簧单摆模 型进而解算了人体踝关节的关节力矩。由于对 从对于肌肉骨骼模型的叙述可知，肌电信号 动力学模型精确度的依赖较高，上述方法难以提 是由神经信号刺激肌肉激活所产生的，因此肌电 供准确度较高的关节力矩求解结果。与之相比， 信号可以用来衡量肌肉激活及神经激活。然而其 以表面肌电信号作为输入的人体肌肉骨骼模型， 间的映射关系是高度非线性的，同时肌电信号的 由于肌肉相对于关节的高度冗余，具有能够预测 幅度也受诸多因素影响)，（例如放大器的增益、 肌肉出力、体现拮抗肌肉驱动、反算肌肉激活等 所选用的电极类型、电极相对于肌肉运动点的位 方面的优势，能够体现出更丰富的人体运动及肌 置、肌肉与电极之间的组织关系等)，所以很难直 肉激活的信息。因此近些年来，很多学者提出了 接利用肌电信号精准的计算和预测肌肉激活以及 使用肌电信号求解关节力矩的方法。该方法主要 分为4个过程：1)从肌电信号(sEMG)求解肌肉 肌肉力。 激活(Muscle activation):2)根据神经肌肉骨骼模 为了从肌电信号中获取有效信息以求解关节 型求解肌肉力；3)求解关节力臂；4)参数辨识。 力矩，必须滤除噪声，量化神经及肌肉的激活程 现阶段所使用的模型例如Zajac在1989年提出的 度。为此，在求解肌肉激活和肌肉力之前，先要 模型a或Huxley的更复杂的生物物理模型3 对肌电信号进行预处理，再将预处理后的肌电信 又如Zahalak的模型s-16(Zahalak,1986,2000年) 号转化为神经激活()。预处理过程可以描述 是基于H的经典著作7。本文主要总结了近年 为：1)5~30Hz高通滤波：由于放大器的精度不足 来使用肌电信号以及相关生理学模型对关节力矩 或者电极相对肌肉发生串动会导致采集到的肌电 的求解方法，并阐述了利用神经肌肉骨骼模型对 信号的平均值产生时变的漂移，而这部分信号并 关节力矩求解的基本原理及过程，给出了模型中 不属于肌电信号，因此要消除这些由于采集误差 涉及的一些参数的生理学统计数据。又给出了以 导致的直流偏移或低频噪声。2)全波整流和归 求解关节划分的该研究的主要应用，最后总结了 化（除以最大自主收缩期间获得的峰值整流EMG 这一求解方法目前所存在的问题和之后发展的方向。 值)处理刘。3)低通滤波：肌肉的电信号的频率 基本原理 成分超过了100Hz,但是肌肉力的频率要比这低 得多，因此为了使肌电信号和肌肉力量相关，滤除 人体运动过程中，骨骼肌的收缩是由我们的 高频成分是非常重要的2,182川。4)求解神经激 神经系统进行控制的，当我们的大脑产生的运动 指令经过脊柱中枢神经传递到骨骼肌时，骨骼肌 活：求取神经激活u的方法如下： 开始收缩，与此同时，骨骼肌会产生一定的动作 0=a×t-T d -B1×u(t-1)-B2×ut-2) (1) 电位，该电位在皮肤表面表现即为肌电信号。 我们将通过表面电极获取的肌电信号进行一 式中：u(t)为第1个采样点的神经激活；e为肌电 定的处理。之后利用神经激活模型求解神经激 信号；a、B,、B2分别为神经激活系数；d为时间延 活()。将()代人到肌肉激活模型来求解肌肉 迟；TE为采样时间间隔。其中上述神经激活系数 受到神经信号后的肌肉激活程度α(t),将其代入 满足： 神经肌肉骨骼模型进而来计算肌肉力和关节力矩 B=1+y2yl<1;y2l<1) (2) 臂，进一步得到关节力矩M,之后根据实际测量 B2=YY(lyil 1:lyal <1) (3) 关节力矩M。,对整体肌肉骨骼模型的参数进行优 a-B1-p2=1 (4) 化，最终得到准确的关节力矩。整个流程如图1所示。 J 肌肉激活求解 EMG 神经激活 u(1) 肌肉激活求解 肌电信号 求解 从肌电信号得到神经激活后，便可以对肌肉 神经激活 激活进行计算和估计。其具体过程为 a(t) a(t)=dln(cu(t)+1),0≤(t)<0.3 M R,Ft) a(t)=mu(t)+b, 0.3≤(t)<1 (5) 参数辨识 力矩臂 神经肌肉骨酪模型 肌肉力 式中：(1)是神经激活；a(t)是肌肉激活。肌肉激 运动学参数 活系数c、d、m和b可以根据曲线过渡点及导数 图1关节力矩求解流程 关系同时求解。一些学者又将进行了一定的化 Fig.1 Joint moment solution 简822
人体下肢关节力矩进行解算；韩国科学技术院结 合人体上肢质心的运动规律，建立了弹簧单摆模 型进而解算了人体踝关节的关节力矩[11]。由于对 动力学模型精确度的依赖较高，上述方法难以提 供准确度较高的关节力矩求解结果。与之相比， 以表面肌电信号作为输入的人体肌肉骨骼模型， 由于肌肉相对于关节的高度冗余，具有能够预测 肌肉出力、体现拮抗肌肉驱动、反算肌肉激活等 方面的优势，能够体现出更丰富的人体运动及肌 肉激活的信息。因此近些年来，很多学者提出了 使用肌电信号求解关节力矩的方法。该方法主要 分为 4 个过程：1) 从肌电信号 (sEMG) 求解肌肉 激活 (Muscle activation)；2) 根据神经肌肉骨骼模 型求解肌肉力；3) 求解关节力臂；4) 参数辨识。 现阶段所使用的模型例如 Zajac 在 1989 年提出的 模型[12] 或 Huxley 的更复杂的生物物理模型[13-14] ， 又如 Zahalak 的模型 [15-16] (Zahalak，1986，2000 年) 是基于 Hill 的经典著作[17]。本文主要总结了近年 来使用肌电信号以及相关生理学模型对关节力矩 的求解方法，并阐述了利用神经肌肉骨骼模型对 关节力矩求解的基本原理及过程，给出了模型中 涉及的一些参数的生理学统计数据。又给出了以 求解关节划分的该研究的主要应用，最后总结了 这一求解方法目前所存在的问题和之后发展的方向。 1 基本原理 人体运动过程中，骨骼肌的收缩是由我们的 神经系统进行控制的，当我们的大脑产生的运动 指令经过脊柱中枢神经传递到骨骼肌时，骨骼肌 开始收缩，与此同时，骨骼肌会产生一定的动作 电位，该电位在皮肤表面表现即为肌电信号。 我们将通过表面电极获取的肌电信号进行一 定的处理。之后利用神经激活模型求解神经激 活 u(t)。将 u(t) 代入到肌肉激活模型来求解肌肉 受到神经信号后的肌肉激活程度 a(t)，将其代入 神经肌肉骨骼模型进而来计算肌肉力和关节力矩 臂，进一步得到关节力矩 M，之后根据实际测量 关节力矩 Mc，对整体肌肉骨骼模型的参数进行优 化，最终得到准确的关节力矩。整个流程如图1所示。 神经激活 求解 EMG 肌肉激活求解 参数辨识 神经肌肉骨骼模型 Mc M 运动学参数 肌电信号 神经激活 a(t) u(t) R, F(t) 力矩臂， 肌肉力 图 1 关节力矩求解流程 Fig. 1 Joint moment solution 2 神经激活求解 从对于肌肉骨骼模型的叙述可知，肌电信号 是由神经信号刺激肌肉激活所产生的，因此肌电 信号可以用来衡量肌肉激活及神经激活。然而其 间的映射关系是高度非线性的，同时肌电信号的 幅度也受诸多因素影响[18] ，(例如放大器的增益、 所选用的电极类型、电极相对于肌肉运动点的位 置、肌肉与电极之间的组织关系等)，所以很难直 接利用肌电信号精准的计算和预测肌肉激活以及 肌肉力。 为了从肌电信号中获取有效信息以求解关节 力矩，必须滤除噪声，量化神经及肌肉的激活程 度。为此，在求解肌肉激活和肌肉力之前，先要 对肌电信号进行预处理，再将预处理后的肌电信 号转化为神经激活 u(t)。预处理过程可以描述 为：1) 5~30 Hz 高通滤波：由于放大器的精度不足 或者电极相对肌肉发生串动会导致采集到的肌电 信号的平均值产生时变的漂移，而这部分信号并 不属于肌电信号，因此要消除这些由于采集误差 导致的直流偏移或低频噪声。2) 全波整流和归一 化 (除以最大自主收缩期间获得的峰值整流 EMG 值) 处理[19-20]。3) 低通滤波：肌肉的电信号的频率 成分超过了 100 Hz，但是肌肉力的频率要比这低 得多，因此为了使肌电信号和肌肉力量相关，滤除 高频成分是非常重要的[12, 18, 21]。4) 求解神经激 活：求取神经激活 u 的方法如下： u(t) = α×e ( t− d TE ) −β1 ×u(t−1)−β2 ×u(t−2) (1) α、β1、β2 TE 式中：u(t) 为第 t 个采样点的神经激活；e 为肌电 信号； 分别为神经激活系数；d 为时间延 迟； 为采样时间间隔。其中上述神经激活系数 满足： β1 = γ1 +γ2(|γ1| < 1;|γ2| < 1) (2) β2 = γ1 · γ2 (|γ1| < 1;|γ2| < 1) (3) α−β1 −β2 = 1 (4) 3 肌肉激活求解 从肌电信号得到神经激活后，便可以对肌肉 激活进行计算和估计。其具体过程为 { a(t) = d ln(cu(t)+1), 0 ⩽ u(t) < 0.3 a(t) = mu(t)+b, 0.3 ⩽ u(t) < 1 (5) 式中：u(t) 是神经激活；a(t) 是肌肉激活。肌肉激 活系数 c、d、m 和 b 可以根据曲线过渡点及导数 关系同时求解。一些学者又将进行了一定的化 简 [8, 22-23] ： ·194· 智 能 系 统 学 报 第 15 卷
第2期 姜峰，等：利用肌电信号求解关节力矩的研究及应用综述 ·195· eAur)-1 a()= (6) 单元CE和一个被动收缩单元PE来近似骨骼肌 e4-1 的工作情况。 其中，A(-3,0)是非线性因子。对于肌肉激活的 整个模型对于骨骼肌力量的求解可以由式(T) 求解还有一些方法主要如表1所示。 进行表示： 表1肌肉激活方法列举 Fm(0,)=F=[F项+FP]cosφ= Table 1 Muscle activation method [fa(D)f(v)a(t)F+fp(D)F]cos (7) 序号相关文献 公式 式中：F为肌腱力；F为最优肌肉纤维下肌肉 Manal3 e4w-1 力；f(④为肌肉主动力-长度关系：f(④为肌肉被 1 a(t)= eA-1 动力-肌肉长度关系B;w为肌肉主动力-速度 A0-1 关系23划 2 Cavallaros a(t)= A1-1 w=03085-Am =03085+Asm周 Manal25] a6-1 m= 4o-1 c=1-m eA:-1 A2 In(Bu(t)+1). 0≤(t)<o B= mu(t)+c, lo≤(t)<1 Chadwick27 u(t) 1-u(t) a(t)= 十 (u(t)-a(t) AL Az 2 a(t)u(t) 十 u(t)≥a(t) RengifoP图 a(t)= A 5 a(t)u(t) (t)<a(t) A2A2 图2肌肉骨骼模型 6 Chadwick a(t)=(A1u(t)+A2)(u(t)-a(t) Fig.2 Musculoskeletal model u(t)-a(t) a(t)= 4.1肌肉力与肌肉纤维长度关系 Ta Thelen A1(0.5+1.5a(t). u(t)>a(t) 根据肌肉骨骼模型来看，肌肉力的大小和肌 A2 肉纤维长度是有一定关系的，这里将这种关系分 0.5+1.5a(0) u(t)≤a(t) 别描述为f)和f()。对于这两种肌肉力同肌肉 表1中列举的公式中Manal在文献[23]中提 纤维长度关系是非线性的，相关学者给出了这两 出的方法相对简便，求解的肌肉激活准确性较 种关系的描述形式。根据相关文献[18,20,32-34] 高，受到广泛的使用。Cavallaro在文献[25]中提 主动收缩力同肌肉纤维长度的关系可以描述为 出的方法是对Manal方法的简化，其中参数A,是 fa(I)=sin(b P+b2l+b3) Manal方法中指数常数项的估计，这里会进一步 (b1=-1.317,b2=-0.403,b3=2.454) (8) 引入误差，因此准确性较Manal方法低，至于 另外一种求解方式为 Manal在文献[26]中提出的方法是针对于特定肌 f(0= ∫6o+61l+62P0.5≤1≤1.5 0 其他 (9) 肉激活求解所提出的方法，该方法进一步提高了 计算精度，但是相对比较繁琐，适用广泛性不强。 被动收缩力同肌肉纤维长度之间的关系可以 至于Chadwich、Rengifo、Chadwick、Telen提出的 表述为位3 方法为常微分方程，求解比其他方法繁琐，使用 fi0=Ap-e学-1] (10) 相对较少。 Ap=0.129,Ke=4.525 可以简化为 4神经肌肉骨骼模型 e10-1) ()e-el-is (11) 在获得肌肉激活及相关人体运动学信息后， 这里所提到的I即肌肉纤维长度是经过归一 便可以利用肌肉骨骼模型对肌肉力以及关节力矩 化之后的长度，即1="/，式中的严是指肌肉纤 进行求解。肌肉骨骼模型如图2所示，其描述了 维长度，：是指最佳肌肉纤维长度，最佳纤维长度 肌肉肌腱的力的作用关系，其中用一个主动收缩 同肌肉激活度具有一定的关系，该关系可以描述为
a(t) = e Au(t)−1 e A −1 (6) 其中，A(−3，0) 是非线性因子。对于肌肉激活的 求解还有一些方法主要如表 1 所示[24]。 表 1 肌肉激活方法列举 Table 1 Muscle activation method 序号 相关文献 公式 1 Manal[23] a(t) = e A1 u(t) −1 e A1 −1 2 Cavallaro[25] a(t) = A u(t) 1 −1 A1 −1 3 Manal[26] u0 = 0.308 5− A1 cos( π 4 ) a0 = 0.308 5+ A1 sin( π 4 ) m = a0 −1 u0 −1 c = 1−m β = e a0 /A2 −1 u0 { A2 ln(βu(t)+1), 0 ⩽ u(t) < u0 mu(t)+c, u0 ⩽ u(t) < 1 4 Chadwick[27] a˙(t) = ( u(t) A1 + 1−u(t) A2 ) (u(t)−a(t)) 5 Rengifo[28] a˙(t) = − a(t) A1 + u(t) A1 u(t) ⩾ a(t) − a(t) A2 + u(t) A2 u(t) < a(t) 6 Chadwick[29] a˙(t) = (A1u(t)+ A2)(u(t)−a(t)) 7 Thelen[30] a˙(t) = u(t)−a(t) Ta Ta = A1(0.5+1.5a(t)), u(t) > a(t) A2 0.5+1.5a(t) , u(t) ⩽ a(t) 表 1 中列举的公式中 Manal 在文献 [23] 中提 出的方法相对简便，求解的肌肉激活准确性较 高，受到广泛的使用。Cavallaro 在文献 [25] 中提 出的方法是对 Manal 方法的简化，其中参数 A1 是 Manal 方法中指数常数项的估计，这里会进一步 引入误差，因此准确性较 Manal 方法低，至于 Manal 在文献 [26] 中提出的方法是针对于特定肌 肉激活求解所提出的方法，该方法进一步提高了 计算精度，但是相对比较繁琐，适用广泛性不强。 至于 Chadwich、Rengifo、Chadwick、Telen 提出的 方法为常微分方程，求解比其他方法繁琐，使用 相对较少。 4 神经肌肉骨骼模型 在获得肌肉激活及相关人体运动学信息后， 便可以利用肌肉骨骼模型对肌肉力以及关节力矩 进行求解。肌肉骨骼模型如图 2 所示，其描述了 肌肉肌腱的力的作用关系，其中用一个主动收缩 单元 CE 和一个被动收缩单元 PE 来近似骨骼肌 的工作情况。 整个模型对于骨骼肌力量的求解可以由式 (7) 进行表示： F mt(θ,t) = F t = [F m A + F m P ] cos ϕ = [fA(l)f(v)a(t)F m o + fP(l)F m o ] cos ϕ (7) F t F m o fA(l) fP(l) f(v) 式中： 为肌腱力； 为最优肌肉纤维下肌肉 力； 为肌肉主动力−长度关系； 为肌肉被 动力−肌肉长度关系[31] ； 为肌肉主动力−速度 关系[12-13, 32]。 CE PE F mt l mt F mt l m l t /2 l t /2 FP m FA m F m ϕ 图 2 肌肉骨骼模型 Fig. 2 Musculoskeletal model 4.1 肌肉力与肌肉纤维长度关系 fA(l) fP(l) 根据肌肉骨骼模型来看，肌肉力的大小和肌 肉纤维长度是有一定关系的，这里将这种关系分 别描述为 和 。对于这两种肌肉力同肌肉 纤维长度关系是非线性的，相关学者给出了这两 种关系的描述形式。根据相关文献 [18, 20, 32-34] 主动收缩力同肌肉纤维长度的关系可以描述为 fA(l) = sin(b1l 2 +b2l+b3) (b1 = −1.317,b2 = −0.403,b3 = 2.454) (8) 另外一种求解方式为 fA(l) = { δ0 +δ1 l+δ2 l 2 0.5 ⩽ l ⩽ 1.5 0 其他 (9) 被动收缩力同肌肉纤维长度之间的关系可以 表述为[32, 35] fP(l) = AP ·[eKpe l m−l m o lm o −1] AP = 0.129,Kpe = 4.525 (10) 可以简化为 fP(l) = e 10(l−1) e 5 = e 10l−15 (11) l = l m /l m o l m l m o 这里所提到的 l 即肌肉纤维长度是经过归一 化之后的长度，即 ，式中的 是指肌肉纤 维长度， 是指最佳肌肉纤维长度，最佳纤维长度 同肌肉激活度具有一定的关系，该关系可以描述为 第 2 期 姜峰，等：利用肌电信号求解关节力矩的研究及应用综述 ·195·
·196· 智能系统学报 第15卷 ()=l(d(1-a(t)+1) (12) 4.3肌腱力求解 式中：为相关系数；。为在最大肌肉激活程度时 求得肌肉激活之后，便可以利用肌肉激活对 的最佳肌肉纤维长度。 肌肉力F(代表肌肉个数)进行求解。肌肉产生 4.2肌肉力与肌肉纤维收缩速度关系 力量的同时，和它相连的肌腱也开始产生承载负 肌肉力量的大小除了和肌肉纤维长度有关系 荷，并将力从肌肉传递到骨骼，这个力称为肌腱 外，它还受到肌肉纤维收缩速度的影响2,3536。 力F。因为肌腱与肌肉是串联的，所以任何通过 Hl在1938年建立肌肉骨骼模型最初是为了研 肌肉的力也必须通过肌腱，反之亦然。肌腱是具 究肌肉收缩和热量之间的关系切，H通过实验 有弹性的被动元件。在肌腱松弛长度飞以下，肌 发现，当肌肉缩短一个距离x时，伴随肌肉释放出 腱不承载任何负荷。然而，当肌腱束长度超过松 一个收缩热量（川，因此有 弛长度时，它产生的力与拉伸距离成比例。当肌 H=ax (13) 肉产生最大等距力F时，肌腱的应变为33%，当 式中：α是与肌肉横截面相关的热常数。同时 力为3.5F时，肌腱的应变为10%。肌腱束应变 Hill还推测a/F是一个近似于0.25的常值，F 可定义为 是指在最佳纤维长度下的最大肌肉力。从能量守 E=-4 (20) 恒的角度考虑，总的能量消耗E等于肌肉耗散的 式中：ε为肌腱应变；为肌腱长度；为肌腱松弛长度 热量和肌肉所做功的总和，即 这里值得注意的是，肌腱应变不会为负值，由 E=F"x+H=(Fm+a)x (14) 等式右边对时间求导得 式(20)可以看到，只有当肌腱束长度大于肌腱束 (F+ar」 松弛长度时，力才会随应变变化，否则肌腱束力 di=(F+a)v (15) 为零。肌腱应变与肌腱力的变化关系由式(21)及图3 式中：m为肌肉纤维收缩速度，之后Hl提出能 所示。 量的消耗应当与力的变化成正比，于是有 0, E≤0 1480.3s2, 0<ε<0.0127 (Fm+a)vm=b(F-Fm) (21) 37.5s-0.2375,g≥0.0127 整理得 Fm=Fab-avm 2 (16) b+ym 到这里肌肉力和肌肉纤维收缩速度的关系是 z1.5 「非线性区 线性区 建立在最佳纤维长度的情况下的，没有考虑到纤 F 维长度变化所带来的影响。下面给出变纤维长度 的肌肉力和肌肉收缩速度之间的关系： 哭0.5 Fm= Fb-avm" b+四(0 (17) 1.27 3.33 2 4 56 对于肌肉的舒张时的情况为 肌腱应变% =FF-(Fi-1)btav (18) 图3肌腱力与应变关系 Fig.3 Relationship between tendon stress and strain 式中：d和b是变心系数；Fc是肌肉力的乘数， 4.4羽状角求解 Epstein和Herzog在1998年提出其值在1.1~l.8之 间”。对上述肌肉力同速度的关系进行简化 如图2所示，由于肌腱与肌肉纤维串联，肌腱 力F和肌肉力F"存在下述关系： 可得： Fr=Fmcos中 (22) 其中，中是羽状角，即肌腱和肌肉纤维之间的角 f(v)= -<0 度，对于大多数的肌肉，羽状角对于肌肉力的影 ⑧*03 Vak 响是可以忽略不计的，但是对于一些转动角度大 (19) 的关节骨骼肌，羽状角对其力的大小还是有很大 2.34 -+0.039 f= 影响的。 13 -≥0 羽状角的大小随关节的活动及肌肉的收缩而 -+0.039 Vmax 变化。KAWAKAMIY等在1998年利用超声波发现， 其中x为最大肌肉收缩力。 内侧腓肠肌羽状角可根据关节角度和肌肉激活量
l m o (t) = lo(λ(1−a(t))+1) (12) 式中： λ 为相关系数； lo 为在最大肌肉激活程度时 的最佳肌肉纤维长度。 4.2 肌肉力与肌肉纤维收缩速度关系 肌肉力量的大小除了和肌肉纤维长度有关系 外，它还受到肌肉纤维收缩速度的影响[12, 35-36]。 Hill 在 1938 年建立肌肉骨骼模型最初是为了研 究肌肉收缩和热量之间的关系[17] ，Hill 通过实验 发现，当肌肉缩短一个距离 x 时，伴随肌肉释放出 一个收缩热量 (H)，因此有 H = ax (13) a/F m o F m o 式中： a 是与肌肉横截面相关的热常数。同时 Hill 还推测 是一个近似于 0.25 的常值， 是指在最佳纤维长度下的最大肌肉力。从能量守 恒的角度考虑，总的能量消耗 E 等于肌肉耗散的 热量和肌肉所做功的总和，即 E = F m x+ H = (F m +a)x (14) 等式右边对时间求导得 (F m +a) dx dt = (F m +a)v m (15) v 式中： m 为肌肉纤维收缩速度，之后 Hill 提出能 量的消耗应当与力的变化成正比，于是有 (F m +a)v m = b(F m o − F m ) 整理得 F m = F m o b−avm b+v m (16) 到这里肌肉力和肌肉纤维收缩速度的关系是 建立在最佳纤维长度的情况下的，没有考虑到纤 维长度变化所带来的影响。下面给出变纤维长度 的肌肉力和肌肉收缩速度之间的关系： F m = ( F m o b−avm b+v m ) fA(l) (17) 对于肌肉的舒张时的情况为 F m = ( F m EccF m o −(F m Ecc −1) F m o b ′ +a ′ v m b ′ −v m ) fA(l) (18) a ′ b ′ F m 式中： 和 是变心系数； Ecc 是肌肉力的乘数， Epstein 和 Herzog 在 1998 年提出其值在 1.1~1.8 之 间 [ 3 7 ]。对上述肌肉力同速度的关系进行简化 可得： f(v) = 0.3 [ v m v m max +1 ] − v m v m max +0.3 , v m v m max < 0 f(v) = 2.34 v m v m max +0.039 1.3 v m v m max +0.039 , v m v m max ⩾ 0 (19) v m 其中 max 为最大肌肉收缩力。 4.3 肌腱力求解 F m i F t l t s F m o F m o 求得肌肉激活之后，便可以利用肌肉激活对 肌肉力 (i 代表肌肉个数) 进行求解。肌肉产生 力量的同时，和它相连的肌腱也开始产生承载负 荷，并将力从肌肉传递到骨骼，这个力称为肌腱 力 。因为肌腱与肌肉是串联的，所以任何通过 肌肉的力也必须通过肌腱，反之亦然。肌腱是具 有弹性的被动元件。在肌腱松弛长度 以下，肌 腱不承载任何负荷。然而，当肌腱束长度超过松 弛长度时，它产生的力与拉伸距离成比例。当肌 肉产生最大等距力 时，肌腱的应变为 3.3%，当 力为 3.5 时，肌腱的应变为 10%。肌腱束应变 可定义为 ε t = l t −l t s l t s (20) ε t l t l t 式中： 为肌腱应变； 为肌腱长度； s为肌腱松弛长度 这里值得注意的是，肌腱应变不会为负值，由 式 (20) 可以看到，只有当肌腱束长度大于肌腱束 松弛长度时，力才会随应变变化，否则肌腱束力 为零。肌腱应变与肌腱力的变化关系由式 (21) 及图 3 所示。 F t = 0, ε ⩽ 0 1 480.3ε 2 , 0 < ε < 0.012 7 37.5ε−0.237 5, ε ⩾ 0.012 7 (21) 标准化肌腱力 F m o /N 肌腱应变/% 2 0 1 2 3 4 5 6 1.5 1 0.5 0 非线性区 线性区 Fo m 1.27 3.33 图 3 肌腱力与应变关系 Fig. 3 Relationship between tendon stress and strain 4.4 羽状角求解 如图 2 所示，由于肌腱与肌肉纤维串联，肌腱 力 F t 和肌肉力 F m 存在下述关系： F t = F m cos ϕ (22) 其中， ϕ 是羽状角，即肌腱和肌肉纤维之间的角 度，对于大多数的肌肉，羽状角对于肌肉力的影 响是可以忽略不计的，但是对于一些转动角度大 的关节骨骼肌，羽状角对其力的大小还是有很大 影响的。 羽状角的大小随关节的活动及肌肉的收缩而 变化。KAWAKAMI Y 等 [38] 在 1998 年利用超声波发现， 内侧腓肠肌羽状角可根据关节角度和肌肉激活量 ·196· 智 能 系 统 学 报 第 15 卷
第2期 姜峰，等：利用肌电信号求解关节力矩的研究及应用综述 ·197· 从22°变为67°。虽然有一些简单的模型描述了羽 F-fF cos f(v)= (28) 状角随肌肉激活的变化，但很少有人通过体内成 fA(0a(t)Fm cos中 像研究（如超声）来证实这些变化。通过对动物 在设定一个初始的肌肉纤维长度后，便可根 肌肉的研究，一些研究人员已经提出了相对比较 据式(8)~(11)求出肌肉力和肌肉纤维长度的关 复杂的模型（如Woittiez等）和一些简化的模型o 系；根据式(12)和(23)可以对羽状角进行相应的 (Scotte&Winter)。为便于计算，利用简化模型来预 估计，又可以根据式(28)求出肌肉力-肌肉收缩 测收缩肌肉的摆动角度，记采样点为1时的羽状 速度关系，之后根据式(19)得到肌肉纤维速度， 角为()，其求解方程为 这时候利用数值积分的方法可以推出下一采样 sin o 点的肌肉纤维长度。将该时刻的肌肉纤维长度 (0）=sin (23) (t) 以及关节角度代入到式(24)和(26)后，可以求得 式中中。为肌肉最佳纤维长度时的羽状角。 此时的肌腱长度，再利用式(20)和(21)未求出该 4.5肌肉肌腱长度 时刻的F,再代入到式（⑦），进行下一步的计算， 从上文对于肌肉肌腱力的论述来看，对于肌 不断重复此过程，最终得到测试时间段的肌肉力。 肉力的求解来说，肌腱长度的大小的求解是至关 4.7关节力矩求解 重要的，其值可以根据式(24)得到： M'(0,t)= m=r+严cos中 (24) c..Fo.m (29) i=1 即 式中：M(0,0)代表关节j在关节角度为0，采样时 =m-"cos中 (25) 刻为t时的关节力矩，：(0为第i块骨骼肌对关节 从式(21)和式(25)可以看出，肌肉肌腱长度 j在关节角度为0时的关节力矩臂，F"(0,)为第 对于求解肌腱长度和肌腱力是必不可少的，但是 i块骨骼肌对关节j在关节角度为0时的肌肉力， 肌肉肌腱长度属于人体生物参数，求解较为困 其中对于力矩臂的求解可见式(30)： 难。Menegaldo等a在2004年提出了一种用于计 算下肢肌肉肌腱长度和关键力矩臂的方法。其研 0=O (30) 80 究参照了Dep等24提出的生理模型所求出的 生理学参数。其具体的计算方法如下： 5参数辨识 L(Q1,Q2,Q3,Q4)=a1+a2fi(Q1,Q2,Q3,Q4)+ 如上所述，d、1、2、A是求解肌肉激活的4个 a3f5(Q1,Q2,Q3,Q4)+…+anf-1(21,Q2,Q3,Q) (26) 主要参数，假设所有参与的肌肉都具有相同的肌 式中：L表示肌肉肌腱长度：a,表示为待定系数： 肉激活求解过程。因此，只需要调整4个激活参 ∫是指泛型非线性系数，根据拟合曲面的形状进 数。在肌肉肌腱力求解过程中，、：、、mx、中 行选择。而Q1、Q2、Q、Q4等则是根据Delp模 是5个反映肌肉力大小的生理学参数，而在这些 型所假设的关节角度数。之后利用测定的关节角 参数中F。、%、?对肌肉力的影响更为显著444， 度信息，进行拟合求解的方式即可求出相应的肌 由于对整体肌肉力的影响较小，可以近似表 肉肌腱长度。另一种较为简便的拟合求解方式为 示为=10g。由于不同肌肉的这3个参数差别 Im =bo+b10 (27) 较大，因此每增加一块肌肉就需要求解3个参数， 式中：b和b1是待定参数；0为关节角度，可以使 诸多的参数使得整个过程变得更加繁琐，因此对 用惯性传感器或光学动捕等方式进行测得。 参数进行优化和化简的工作就显得尤为重要。 4.6肌肉力整合求解 Ai等进行了两组优化，第1组根据受试者 通过对整个过程的论述，在这里重新讨论最 特定的肌肉骨骼几何模型进行估计：第2组采用 初给定的式（⑦），这时式中的大部分未知量已经在 通用的肌肉骨骼模型。之后将两组原始肌电信号 上述内容给出了求解方法，然而现在依然不能依 和运动学数据都输入到优化模型（非线性最小二 靠式(7)直接给出对于肌肉力量的一个估计，因 乘和levenberg-marquardt算法)中，使得估计力矩 为肌肉纤维长度的大小并不知道。 和参考力矩差值最小。Menegaldo等4I采用式 Fm (0,t)=F=[F+F]cos= [fa(⑩f(v)a(t)F+fp(OF]cos中 (31)对参数进行优化评价指标： 对于肌肉纤维长度的求解需要用到数值积分 (TM()-TS(0)2 的方法，因此这里首先根据上述方程（⑦）进行对 RMSE= 肌肉力-肌肉收缩速度关系进行求解： TMMAX ×100%(31)
ϕ(t) 从 22°变为 67°。虽然有一些简单的模型描述了羽 状角随肌肉激活的变化，但很少有人通过体内成 像研究 (如超声) 来证实这些变化。通过对动物 肌肉的研究，一些研究人员已经提出了相对比较 复杂的模型[39] (如 Woittiez 等) 和一些简化的模型[40] (Scott&Winter)。为便于计算，利用简化模型来预 测收缩肌肉的摆动角度，记采样点为 t 时的羽状 角为 ，其求解方程为 ϕ(t) = sin−1 ( l m o sinϕ0 l m(t) ) (23) 式中 ϕo 为肌肉最佳纤维长度时的羽状角。 4.5 肌肉肌腱长度 从上文对于肌肉肌腱力的论述来看，对于肌 肉力的求解来说，肌腱长度的大小的求解是至关 重要的，其值可以根据式 (24) 得到： l mt = l t +l m cos ϕ (24) 即 l t = l mt −l m cos ϕ (25) 从式 (21) 和式 (25) 可以看出，肌肉肌腱长度 对于求解肌腱长度和肌腱力是必不可少的，但是 肌肉肌腱长度属于人体生物参数，求解较为困 难。Menegaldo 等 [41] 在 2004 年提出了一种用于计 算下肢肌肉肌腱长度和关键力矩臂的方法。其研 究参照了 Delp 等 [42-43] 提出的生理模型所求出的 生理学参数。其具体的计算方法如下： L(Q1,Q2,Q3,Q4) = a1 +a2 f1(Q1,Q2,Q3,Q4)+ a3 f2(Q1,Q2,Q3,Q4)+···+an fn−1(Q1,Q2,Q3,Q4) (26) Q1、Q2、Q3、Q4 式中：L 表示肌肉肌腱长度；ai 表示为待定系数； fi 是指泛型非线性系数，根据拟合曲面的形状进 行选择。而 等则是根据 Delp 模 型所假设的关节角度数。之后利用测定的关节角 度信息，进行拟合求解的方式即可求出相应的肌 肉肌腱长度。另一种较为简便的拟合求解方式为 l mt = b0 +b1θ (27) 式中：b0 和 b1 是待定参数；θ 为关节角度，可以使 用惯性传感器或光学动捕等方式进行测得。 4.6 肌肉力整合求解 通过对整个过程的论述，在这里重新讨论最 初给定的式 (7)，这时式中的大部分未知量已经在 上述内容给出了求解方法，然而现在依然不能依 靠式 (7) 直接给出对于肌肉力量的一个估计，因 为肌肉纤维长度的大小并不知道。 F mt(θ,t) = F t = [F m A + F m P ] cosϕ = [fA(l)f(v)a(t)F m o + fP(l)F m o ] cos ϕ 对于肌肉纤维长度的求解需要用到数值积分 的方法，因此这里首先根据上述方程 (7) 进行对 肌肉力−肌肉收缩速度关系进行求解： f(v) = F t − f l pF m o cos ϕ fA(l)a(t)Fm o cos ϕ (28) 在设定一个初始的肌肉纤维长度后，便可根 据式 (8)~(11) 求出肌肉力和肌肉纤维长度的关 系；根据式 (12) 和 (23) 可以对羽状角进行相应的 估计，又可以根据式 (28) 求出肌肉力−肌肉收缩 速度关系，之后根据式 (19) 得到肌肉纤维速度， 这时候利用数值积分的方法可以推出下一采样 点的肌肉纤维长度。将该时刻的肌肉纤维长度 以及关节角度代入到式 (24) 和 (26) 后，可以求得 此时的肌腱长度，再利用式 (20) 和 (21) 未求出该 时刻的 F t ，再代入到式 (7)，进行下一步的计算， 不断重复此过程，最终得到测试时间段的肌肉力。 4.7 关节力矩求解 Mj (θ,t) = ∑m i=1 (ri(θ)· F mt i (θ,t)) (29) Mj (θ,t) ri(θ) F mt i (θ,t) 式中： 代表关节 j 在关节角度为 θ，采样时 刻为 t 时的关节力矩， 为第 i 块骨骼肌对关节 j 在关节角度为 θ 时的关节力矩臂， 为第 i 块骨骼肌对关节 j 在关节角度为 θ 时的肌肉力， 其中对于力矩臂的求解可见式 (30)： r(θ) = ∂l mt(θ) ∂θ (30) 5 参数辨识 F m o l m o l t s v m max ϕ F m o l m o l t s v m max v m max = 10l m o 如上所述，d、γ1、γ2、A 是求解肌肉激活的 4 个 主要参数，假设所有参与的肌肉都具有相同的肌 肉激活求解过程。因此，只需要调整 4 个激活参 数。在肌肉肌腱力求解过程中， 、 、 、 、 是 5 个反映肌肉力大小的生理学参数，而在这些 参数中 、 、 对肌肉力的影响更为显著[44-45] ， 由于 对整体肌肉力的影响较小，可以近似表 示为 。由于不同肌肉的这 3 个参数差别 较大，因此每增加一块肌肉就需要求解 3 个参数， 诸多的参数使得整个过程变得更加繁琐，因此对 参数进行优化和化简的工作就显得尤为重要。 Ai 等 [46] 进行了两组优化，第 1 组根据受试者 特定的肌肉骨骼几何模型进行估计；第 2 组采用 通用的肌肉骨骼模型。之后将两组原始肌电信号 和运动学数据都输入到优化模型 (非线性最小二 乘和 levenberg-marquardt 算法) 中，使得估计力矩 和参考力矩差值最小。Menegaldo 等 [47] 采用式 (31) 对参数进行优化评价指标： RMSE = 1 TMMAX vuuuuut∑N t=1 (TM(t)−TS(t))2 N ×100% (31) 第 2 期 姜峰，等：利用肌电信号求解关节力矩的研究及应用综述 ·197·