第60 卷 第4 期 2024 年7 月

Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 60, No. 4 (July 2024)

doi: 10.13209/j.0479-8023.2024.050

国家自然科学基金(62073038)和北京理工大学高层次人才研究基金资助

收稿日期: 2023–06–19;

修回日期: 2023–09–22

基于肌肉协同和肌肉力相对关系的人体行走下肢肌肉力计算方法

张岭1 黄岩1,2,3,† 王启宁4

1.北京理工大学机电学院, 北京 100081; 2.北京理工大学智能机器人与系统高精尖创新中心, 北京 100081;3.仿生机器人与系统教育部重点实验室, 北京 100081; 4.北京大学工学院先进制造与机器人系, 北京100871;† 通信作者, E-mail: yanhuang@bit.edu.cn

摘要 针对现有计算肌肉力的方法中缺乏对运动过程中肌肉协同机制和肌肉力相对关系考量的问题, 提出一种基于行走过程中肌肉力相对关系与肌肉协同的肌肉力计算方法。该方法在二次优化的基础上, 将人体行走过程中一段时间的肌肉力以及不同肌肉收缩力之间的关系纳入逆动力学计算中。在所建立的肌肉骨骼模型的基础上, 利用该方法得到行走过程中下肢主要肌肉的肌肉力, 计算结果可以有效地反映人体运动规律以及行走过程中的肌肉协同机制。该计算方法可为人体运动形成机制、运动康复与训练以及仿人机器控制等提供理论依据。

关键词 人体肌肉骨骼模型; 逆动力学计算; 肌肉协同; 肌肉力相对关系

人体运动过程中肌肉力的计算对人体运动规律研究、双足运动机理分析、人体运动康复与训练指导、穿戴式机器人的设计以及仿人机器人的运动控制等具有重要意义[1]。但是, 目前关于肌肉力的计算依旧是一个难题[2]

首先, 难以对肌肉进行全面的动力学建模和测量。因为肌肉组织是柔性体, 且位于人体上皮组织下方, 附着在骨骼之上, 形态和结构比较复杂, 很难得到精确的动力学模型。在现有的技术条件下,难以在人体运动过程中进行侵入式测量和分析。因此, 难以利用力传感器直接测量肌肉力[3]

第二, 在人体运动过程中, 神经系统将电生理信号传给肌肉细胞, 每个肌肉细胞在电信号的刺激下收缩, 使肌肉整体产生收缩力[4]。因此, 肌肉力的产生涉及复杂的神经生理学过程, 目前对相关机制的研究还有很多待解决的问题。因此, 难以通过神经生理学分析等方法计算肌肉力。

第三, 关节力矩由肌肉力的合力形成。一般来说, 肌肉的数量远大于关节数量, 因此, 从关节力矩求解肌肉力的过程是一个冗余问题[5], 有无数组肌肉力的组合可以形成特定的关节力矩。因此, 即使得到关节力矩, 也难以通过解动力学方程的方式计算肌肉力。

本文提出一种基于肌肉协同机制和肌肉力相对关系的人体肌肉力计算方法, 可以由关节力矩和运动学数据得到符合人体运动规律的肌肉力。将该方法应用在人体行走运动的分析中, 计算得到人体行走过程中下肢主要肌肉的肌肉力, 并对计算结果进行验证和讨论。

1 人体肌肉力计算方法

目前, 肌肉力计算方法大致分为基于表面肌电信号的计算方法、正动力学法以及逆动力学法[6]

基于表面肌电信号的方法是较常见的肌肉力计算方法[7], 该方法通过测量人体运动过程中的表面肌电信号(sEMG), 并建立从肌电信号到肌肉力的数学模型来进行肌肉力的计算。Bogey 等[8]提出的方法是利用 sEMG 信号计算肌肉力的代表性研究成果, 该方法首先通过肌电信号采集装置, 得到行走过程中踝关节肌肉的表面肌电信号, 并将表面肌电信号转化成肌肉活动性, 然后利用希尔(Hill)肌肉模型[9]建立肌肉活动性到肌肉力的数学模型, 最终得到行走过程中踝关节的肌肉力。该方法是从神经到肌肉的角度得到肌肉力的估计值, 有较强的神经生理学依据。但是, 肌电信号只能反映行走过程中肌肉的收缩程度, 与肌肉力的关系还没有权威的理论支撑。此外, 肌电信号十分不稳定, 易受外界条件影响, 在测量过程中也会有噪声对肌电信号产生干扰, 在处理过程中会丢失很多信息[10]

正动力学方法是将神经信号输入肌肉骨骼模型, 结合肌肉动力学和刚体动力学, 使肌肉骨骼模型生成运动, 并结合动态优化和最优控制等方法,使肌肉骨骼模型的能量耗散最小, 以便得到能使模型稳定行走的肌肉力结果。Davy 等[11]提出的计算肌肉力的方法是正动力学法的代表性研究成果, 该方法建立了摆动腿的肌肉骨骼模型, 通过动态优化的方法, 以行走过程中肌肉力之和最小以及能量耗散最小为目标, 生成稳定行走时摆动腿的运动学数据, 并将此时的肌肉力作为人体行走时摆动腿肌肉力的计算结果。但是, 该方法的计算量太大, 计算速度十分缓慢, 并且在某些情况下得不到可行解。

逆动力学方法一般先基于人体行走的运动学数据和约束力计算出关节力矩, 再用优化算法进行肌肉力的求解。de Groote 等[12]提出的方法是逆动力学法的代表性研究成果, 该方法利用运动采集设备得到人体的的运动学数据和行走时的足底压力, 将这些数据导入肌肉骨骼模型中, 求出行走时的关节力矩, 再利用静态优化法, 通过关节力矩及运动学数据求出肌肉力。该方法的优点是计算速度快, 得到的肌肉力能满足设定的优化目标和约束条件。但是, 这种方法没有考虑肌肉协同机制和行走过程中肌肉力的一些约束条件。村井昭彦[13]建立了一个全身肌肉骨骼模型, 并提出一种将电生理信号和逆动力学优化结合的肌肉力计算方法。该方法按照关节和自由度将肌肉分组, 并在每组划分主要肌肉和次要肌肉, 主要肌肉的肌肉力由 sEMG 信号计算,次要肌肉的肌肉力用逆动力学优化的方法计算, 得到人体全身肌肉的收缩力。但是, 该方法依然需要测量表面肌电信号, 因此主要肌肉的肌肉力会受到肌电信号不稳定的影响。此外, 该方法的肌肉协同关系仅仅通过肌肉分组来实现, 对肌肉协同关系的描述不够充分, 且没有考虑运动过程中不同肌肉之间肌肉力的相对大小, 也没有对下肢节律性运动的肌肉力计算结果进行详细的对比。

针对现有肌肉力计算方法中存在的问题, 本文提出一种基于肌肉协同机制和肌肉力相对关系的肌肉力计算方法。该方法在计算时不需要采集人体表面肌电信号, 可以针对人体全身复杂运动进行求解,能够得到唯一解, 得到的肌肉力符合人体运动规律。将该方法应用到人体行走运动分析中, 根据关节力矩和运动学数据, 得到行走过程中下肢主要肌肉的肌肉力, 并与文献报道的相关结果进行比较,对本文方法进行验证。

2 肌肉骨骼模型的建立与步态划分

人体的肌肉骨骼系统十分复杂, 有些肌肉在行走过程中起到非常重要的作用, 有些肌肉只起到次要的作用, 甚至已经在人类进化的过程中逐渐退化[14]。因此, 需要对人体的结构和肌肉的数量进行合理的简化。此外, 人体行走过程不同阶段的肌肉协同和肌肉力相对关系也各不相同, 不能针对整个步态周期都用同一种肌肉特征来描述, 必须进行相应的步态划分。

2.1 肌肉骨骼模型的建立

本文的目标是计算人体行走时下肢的肌肉力,而行走主要靠下肢两条腿的节律性运动来完成, 因此我们保留人体腿部的自由度而忽略上身以及手臂的自由度, 将上身用一根杆来代替。脚趾在人的行走过程中不起主要作用, 所以忽略脚趾的自由度。此外, 虽然侧向自由度和侧面的肌肉对行走过程中的稳定性也有重要的作用, 但平地行走主要是前后运动的过程, 因此我们忽略髋关节和踝关节的侧向自由度, 重点研究人体在矢状面的自由度。

简化的人体骨骼模型如图1 所示。本研究将人体骨骼系统简化成一个七连杆模型, 用来描述人体行走在矢状面的全部自由度, 连杆分别代表上身、大腿、小腿和足部。

图1 人体骨骼模型示意图
Fig.1 Simplified model of human skeleton system

从产生原因角度, 肌肉力分为主动部分和被动部分, 主动部分由神经信号刺激肌细胞产生, 被动部分由关节角度变化牵拉肌纤维产生。不论是主动部分还是被动部分, 肌肉都只在收缩方向产生拉力,本文只需要研究某块肌肉产生的合力而忽略从肌细胞到肌肉力的过程。所以, 假设肌肉力为一个向量,向量的起点位于肌肉与骨骼的连接点, 方向与肌肉收缩的方向相同。

在肌肉的形态方面, 我们假设肌肉有与骨骼连接的起点和终点, 跨过一个关节的肌肉称为单关节肌肉, 跨过多个关节的肌肉称为多关节肌肉。在人体中, 肌肉一般附着在骨骼上, 如果只设置起点和终点, 不仅会使有些肌肉的力臂过大, 还会使有些肌肉在运动过程中出现不符合生理学的状态。所以, 我们在简化过程中设置中间点, 中间点的位置与跨过的关节的距离保持一个定值, 使肌肉的形态更加符合人体生理学原理。

人体有约 300 多块下肢肌肉, 但是对行走起主要作用的肌肉只有一部分。我们选取对行走有重要作用的肌肉, 将其他肌肉忽略或与重要肌肉合并,以肌群的方式对肌肉进行建模。因此, 在本文模型中, 一块肌肉代表包含其本身的一组肌群, 计算得出的肌肉力代表整个肌群的作用。

按照功能, 人体骨骼肌可以分为主动肌和拮抗肌。在人体运动过程中, 主动肌收缩完成动作, 拮抗肌舒张保持关节的稳定性。一般来说, 每个关节都要通过一对主动肌和拮抗肌来完成运动, 本文按照此思路建立下肢运动的主要肌肉。

在踝关节前侧的肌肉主要是胫骨前肌(TA), 其起点位于小腿, 终点位于脚掌前端, 在行走时主要控制摆动腿踝关节的背屈和支撑腿向前的转动。与胫骨前肌对应的拮抗肌有腓肠肌(GAS)和比目鱼肌(SOL), 这两块肌肉的主要部分都位于小腿。比目鱼肌在小腿内侧, 只跨过踝关节, 是一块单关节肌肉, 一般控制踝关节的跖屈; 腓肠肌的主体部分在比目鱼肌外侧, 起点位于膝关节上方, 除控制踝关节的跖屈外, 还参与控制膝关节的屈曲。

人体前侧的大肌肉肌群称为股肌群, 主要由 4 块肌肉组成: 股内侧肌、股外侧肌、股中间肌和股直肌。股内侧肌和股外侧肌分别位于大腿的内侧和外侧, 主要控制髋关节的偏摆, 与本研究关注的二维行走模式关系不大。股直肌和股中间肌都位于大腿前侧, 股直肌位于内部, 被股中间肌完全覆盖。这两块肌肉的作用十分相似, 且位置几乎相同,因此我们将其合并, 用股直肌(RF)命名。股直肌的起点位于膝关节上方, 与髋关节十分接近, 主要作用是控制膝关节的伸展。

与股肌群相对, 大腿的后方也有一块肌肉群,统称为腘绳肌, 主要由半膜肌、半腱肌和股二头肌长头组成。这 3 块肌肉的位置和作用十分相似, 因此我们在大腿后侧建立腘绳肌(HAM), 其起点位于上身, 靠近髋关节, 终点位于膝关节下方, 是一块双关节肌肉, 主要作用是控制膝关节的屈曲。

人体臀部附近跨过髋关节的肌群叫做臀部肌群, 主要包括臀大肌、臀中肌和臀小肌等。臀大肌包裹臀中肌和臀小肌, 在髋关节的伸展中起主要作用, 因此在躯干附近建立臀大肌(GLU)。臀大肌的起点位于上身, 终点位于髋关节下方, 主要作用是维持上身稳定, 并控制髋关节的伸展。与臀大肌对应的位于人体髋关节前侧的肌群称为髂腰肌, 主要包括髂肌和腰大肌, 这两块肌肉位置和作用十分相似, 因此在模型中添加髂腰肌(ILPSO), 其起点位于上身, 终点位于大腿前侧, 主要作用是维持上身稳定并控制髋关节屈曲。

本文在骨骼模型的基础上, 每条腿添加上述 7块具有代表性的肌肉(图2), 从左到右分别为添加小腿附近肌肉、大腿附近肌肉和躯干部分肌肉, 并标注单关节肌肉、双关节肌肉以及肌肉中间点。这 7块肌肉与骨骼的组合可以产生二维行走矢状面上的所有运动。

图2 肌肉骨骼模型
Fig.2 Musculoskeletal model

2.2 行走步态划分

步态周期是行走过程中两个相同事件发生的一段时间[15], 一般将一条腿的脚跟开始触地作为一个步态周期的开始, 将该脚跟再次触地作为步态周期的结束。一个步态周期可分为单脚支撑阶段(single support phase)和双脚支撑阶段(double support phase,DS)[16]。其中, 单脚支撑阶段只有一只脚接触地面,另一条腿处于摆动状态, 双脚支撑阶段两只脚都接触地面, 双腿都作为支撑腿。

本文从行走过程中的运动学数据和肌肉力特点出发, 进一步划分步态(图3), 以便提取不同阶段的肌肉协同和肌肉力相对关系, 有助于肌肉力的计算。在双脚支撑阶段, 主要进行支撑腿和摆动腿的相互转换以及重心的向前移动, 由于整个阶段的生物力学特点类似, 因此本文不对双脚支撑阶段进行进一步的划分。

图3 步态阶段的划分
Fig.3 Division of the human gait

在单脚支撑阶段, 摆动腿先抬起并向前摆动,髋关节屈曲, 关节角逐渐减小, 膝关节先屈曲后伸展, 关节角先减小后变大, 此时腿部肌肉的主要作用是将摆动腿抬起并向前摆动。摆动腿向前摆动一段时间之后, 髋关节会达到行走时的最大角度并开始伸展, 膝关节也会伸展到最大角度, 并维持在这个角度继续运动。此阶段, 肌肉的主要作用是将重心前移。因此, 我们将单脚支撑阶段进一步划分为单脚支撑–膝关节打开阶段(single support-knee open,SSKO)和单脚支撑–膝关节固定阶段(single supportknee close, SSKC)。如图3 所示, 以腿 1 脚跟着地作为一个步态周期的开始和结束, 以前腿脚跟的状态作为步态阶段的分界线, 将一个步态周期划分为6 个阶段。

3 肌肉力计算方法

t 时刻关节的力矩和肌肉力的关系用下式表示:

其中, τ (t)代表 t 时刻的关节力矩, f (t)代表 t 时刻的肌肉力, J(t)代表 t 时刻从肌肉力到关节力矩之间的雅克比矩阵。行走过程中关节角度随时间变化, 因此每一时刻的肌肉力与关节力矩都需满足式(1)。

虽然关节力矩与肌肉力之间有十分明确的关系, 但由于肌肉个数大于自由度的个数, 因此从关节力矩求解肌肉力的过程是一个冗余问题, 式(1)不存在唯一解, 不能通过一般的解动力学方程的方法(如牛顿欧拉法和拉格朗日法等)来求解肌肉力。

对于冗余问题, 可以将其转化成一个动态优化问题, 即设定一个随时间变化的目标函数, 找到每一时刻都能使目标函数最优的解。对于优化目标,首先考虑得到的肌肉力需满足当前的关节力矩, 即找到合适的 f (t), 使得|JT(t)f(t)-(t)|最小。此外,人体往往倾向于消耗最少的能量来完成运动, 除满足关节力矩误差最小的条件外, 还应该使肌肉力之和尽量小, 以便符合人体行走时的生理特征。因此,目标函数的第二部分为其中 n 为参与运动的肌肉的数量, fi(t)为 t 时刻第 i 个肌肉的肌肉力。

除了设定合理的目标函数, 还应该找到优化过程中需要满足的约束条件, 对优化目标进行限制,使解出的肌肉力满足运动规律。在人体运动的过程中, 肌肉通常存在一个最大等距收缩力, 即肌肉在等距收缩时能够产生的最大力, 一般用来代表行走过程中肌肉力的最大值。因此, 肌肉力在运动过程中必须时刻满足0≤f(t)≤Fmax 。表1 列出本文使用的肌肉最大等距收缩力[17]

表1 肌肉最大等距收缩力
Table 1 Maximum isometric contract force of muscle

肌肉名称 最大等距收缩力/N胫骨前肌(TA) 750腓肠肌(GAS) 2900比目鱼肌(SOL) 4000股直肌(RF) 2000腘绳肌(HAM) 2000臀大肌(GLU) 1800髂腰肌(ILPSO) 3000

对于优化方法的选择, 综合考虑计算速度、准确度以及算法的简洁程度, 我们选择顺序二次规划法(SQP)作为计算肌肉力使用的优化算法。

综上所述, 本文计算肌肉力的方法可以总结为找到合适的肌肉力, 使以下目标函数最优:

其中, f(t)为优化变量, 即当前时刻的肌肉力; τj(t)为当前时刻的关节力矩; kf 为目标函数的权重系数, 改变 kf 可以调整目标函数两部分的权重。在优化过程中需满足以下约束条件:

虽然利用优化算法可以得出某一时刻满足优化目标和约束条件的肌肉力, 但从生物力学的角度来看, 在不同步态阶段, 支撑腿和摆动腿有各自的肌肉协同关系, 肌肉力之间也存在约束。虽然目前还没有较为准确的肌肉力计算数据, 但我们综合文献[18–23]的结果, 可以初步将肌肉协同和不同肌肉收缩力之间的关系纳入动态优化过程, 使计算出的肌肉力更加符合人体运动规律。

3.1 优化算法中的肌肉力相对关系

在优化过程中, 设当前时刻 t 的优化变量为 f(t),取 t 之前 n 个时刻((tn)~(t–1))的优化结果进行采样, 得到两组满足特定关系的肌肉力序列 f1(t–1),f1(t–2), …, f1(tn)和 f2(t–1), f2(t–2), …, f2(tn)。n 为采样数, 即取前 n 个时刻的优化结果来描述肌肉力之间的相对关系。本文取最小采样数为 3, 最大采样数为 10。在优化过程中, 设优化变量满足

式(4)描述当前时刻两个特定肌肉力的优化结果与采样结果的平均值之间满足一定的关系, 取平均值的意义是使优化过程更加连续, 避免某一时刻特殊值的影响。kg 为描述两组肌肉的肌肉力相对关系的系数, 本文通过人体行走过程中各肌肉表面肌电信号的数据, 并综合文献[18–19]的肌肉力计算结果来确定 kg 的取值。表2 展示行走一个步态周期 3 个阶段的肌肉力相对关系。

表2 肌肉力相对关系
Table 2 Relationship between specific muscle force

腿 步态阶段DS阶段 SSKO阶段 SSKC阶段ff 2 SOLGAS ff SOL2GAS 腿1 3 2 ff HAMGLU ff SOL2GAS f≤f HAMGLU 3 2 ff 3 2 HAMSOA ff SOL2GAS 腿2 ff —SOL2GAS 1 5 f≤f TAGLU

3.2 优化算法中的肌肉协同机制

在行走过程中, 除不同肌肉的肌肉力满足一定的关系外, 肌肉之间还存在一定的协同关系[24], 而现有的逆动力学优化算法往往欠缺此方面的考量。本文在逆动力学优化的基础上, 将行走过程中的肌肉协同关系纳入求解肌肉力的优化过程, 使肌肉力的计算更加符合人体运动时的生理特性。

肌肉协同指在运动过程中某些肌肉力呈现相同的变化趋势。例如, 在摆动腿抬起时, 股直肌和胫骨前肌的肌肉力都呈现逐渐增大的过程。本文用两组肌肉序列之间相关系数代表肌肉在一段采样时间内的协同关系, 衡量运动过程中两组肌肉力变化趋势的相似程度, 反映两组肌肉力协同性的强弱。

设当前时刻优化结果与之前时刻优化结果采样的序列为

其中, 两块肌肉的采样序列为

将已知变量和优化变量分离:

可以得到 f1f2 的协方差和方差分别为

通过式(13)~(15), 计算得到 f1f2 的相关系数:

相关系数越接近 1, 表明在该采样阶段两组肌肉力的变化趋势越接近; 相关系数越接近 0, 表明两组肌肉力的变化越不相关。

我们将相关系数约束纳入肌肉力的计算中, 即在优化过程中, 具有协同关系的两组肌肉力序列f1(tt-n)和f2(tt-n)满足

其中, kρ 为相关系数的约束系数, 本文认为当两组肌肉力的相关系数大于 0.7 时, 两组肌肉具有较高的协同性。表3 列出本文在优化过程中使用的相关系数约束。

表3 优化过程中的相关系数
Table 3 Coefficient of correlation in optimization

步态阶段腿DS 阶段 SSKO 阶段 SSKC 阶段腿1 (TA,)ffk(,)(,) GASSOL HAM GASSOL ffkffk(,)ffkILPSOGAS(,)ffk腿2 GASSOL GASSOL (,)ffk RFHAM (,)ffk(,)ffkILPSOGAS ILPSOTA (,)ffk TAGLU (,)ffk

3.3 肌肉力计算方法总结

综上所述, 本文提出的肌肉力计算流程如图4所示。计算流程可以分为 3 个部分, 分别为优化参数设定、逆动力学求解和肌肉力计算。优化参数设定部分主要通过人体运动采集装置采集人体行走过程中表面肌电信号等数据, 提取人体行走过程中肌肉力的相对关系和肌肉协同, 确定优化过程中相应约束条件的参数。逆动力学求解部分基于动作捕捉设备采集的运动学数据和肌肉骨骼模型, 用来计算关节力矩。肌肉力计算部分基于关节力矩及得到的约束条件, 优化得到人体行走过程中下肢主要肌肉的肌肉力。

图4 肌肉力计算流程
Fig.4 Process of muscle force calculation

4 仿真验证

为了消除不同被试行走差异的影响, 本文选取人体行走时的关节角度和力矩作为参考数据[22], 对提出的肌肉力计算方法的计算效果进行仿真验证。其中, 力矩为通过人体全身质量归一化后的数值,可以忽略个体之间的差异。

将数据导入肌肉骨骼模型, 利用本文提出的方法进行计算, 得到人行走时一个步态周期一条腿 7块主要肌肉的肌肉力计算结果(图5~11)。肌肉力的单位为%BW, 以人体全身重量为基准进行归一化,并用百分比表示。该单位表示肌肉力相对于个体重量的比例, 用来抵消不同个体肌肉力结果之间的差异, 便于进行比较。

图5 胫骨前肌肌肉力计算结果
Fig.5 Muscle force calculation result of tibialis anterior

(a) 仅考虑肌肉协同的计算结果, kρ=0.3 (协同关系比较弱); (b) 仅考虑肌肉协同的肌肉力优化结果, kρ=0.7 (协同关系比较强);(c) 仅考虑肌肉相对关系; (d) 既考虑肌肉协同又考虑肌肉相对关系; (e) 文献报道的肌肉力计算结果[19–20]; (f) 相关肌肉的肌电信号测量结果[23]。虚线为一个步态周期中各个步态阶段的分界, 步态阶段的顺序与图3 相同。下同

图6 腓肠肌肌肉力计算结果
Fig.6 Muscle force calculation result of gastrocnemius

图7 比目鱼肌肌肉力计算结果
Fig.7 Muscle force calculation result of soleus

图8 股直肌肌肉力计算结果
Fig.8 Muscle force calculation result of rectus foremis

图9 腘绳肌肌肉力计算结果
Fig.9 Muscle force calculation result of hamstrings

图10 臀大肌肌肉力计算结果
Fig.10 Muscle force calculation result of gluteus maximus

图11 髂腰肌肌肉力计算结果
Fig.11 Muscle force calculation result of iliopsoas

首先分析既考虑肌肉协同, 同时考虑肌肉之间相对关系的完整方法计算结果(图5~11 中(d)图), 可以发现, 得到的肌肉力基本上符合人体运动规律。在腿1 脚跟着地到脚趾着地这段时间, 踝关节的角度几乎不变, 胫骨前肌处于等距收缩阶段, 胫骨前肌的肌肉力逐渐增大。之后, 踝关节的角度逐渐增大, 胫骨前肌开始舒张, 肌肉力逐渐减小, 直至减小到 0。

腓肠肌和比目鱼肌都位于小腿后侧, 主要作用是控制踝关节的跖屈。在步态周期开始阶段, 两组肌肉处于放松状态, 肌肉力几乎为 0; 当步态进入支撑中期, 即腿 2 摆动到最高点时, 比目鱼肌和腓肠肌开始进行离心收缩, 肌肉力逐渐增大, 在腿 2脚跟着地之前达到最大值, 之后开始向心收缩, 肌肉力逐渐减小。直到腿 2 脚跟着地, 这两块肌肉的肌肉力才减小到 0。

股直肌在步态周期的开始阶段先进行离心收缩, 肌肉力逐渐增大, 当腿 1 脚趾着地时, 肌肉力达到最大, 之后逐渐减小到 0。当腿 1 脚后跟再次抬起时, 股直肌开始进行向心收缩, 肌肉力也出现一个小的峰值。

在步态周期开始阶段, 髋关节处于伸展状态,腘绳肌离心收缩, 肌肉力逐渐增大。当腿 1 脚尖着地时, 肌肉力达到最大, 之后一直减小, 直至为 0。

臀大肌在步态周期开始阶段先进行等距收缩,肌肉力从零增大到一个较小的峰值。当腿 1 脚趾着地时, 髋关节开始伸展, 臀大肌开始进行向心收缩,肌肉力逐渐减小, 直至为 0。

髂腰肌在步态开始阶段一直处于舒张状态, 直至腿 2 向前摆动到最远, 髂腰肌开始进行离心收缩,肌肉力逐渐增大, 直到一个较大的峰值; 当腿 1 脚趾着地时, 髋关节伸展到极限, 髂腰肌开始进行向心收缩, 肌肉力逐渐减小, 直至为 0。

除在运动过程中符合人体运动规律外, 肌肉力计算结果还反映行走过程中的肌肉协同关系和不同肌肉之间肌肉力的相对关系。例如, 通过图5~11的(a)与(b)图对比可以发现, 当相对系数约束系数增大, 即肌肉协同约束设定更加严格时, 肌肉力会受行走过程中肌肉协同的影响, 从而避免出现一些不符合运动规律的状态。例如, 胫骨前肌肌肉力在有些步态阶段降至 0, 髂腰肌在肌肉力下降阶段的斜率也发生改变, 这些都与人体运动规律相符。

通过对比(c)图和其他结果可以发现, 在添加肌肉力相对关系约束之后, 得到的肌肉力计算结果更符合(e)和(f)图中的肌肉力相对关系。例如, 对比(b)和(c)图, 在添加肌肉力相对关系之前, 虽然腓肠肌和比目鱼腓肠肌的协同关系比较强, 但不符合文献中肌肉力计算结果和 sEMG 信号结果中体现的相对关系。添加相对关系约束后, 在整个步态周期,腓肠肌的肌肉力大致是比目鱼肌肌肉力的一半, 与现有研究结果[20]相符。此外, 在添加肌肉力相对关系之后, 臀大肌的肌肉力明显增大, 髂腰肌肌肉力由一个波峰变为两个波峰, 更加符合人体运动规律以及文献报道的结果[20,24]

综上所述, 本文肌肉力计算方法得到的行走过程中下肢肌肉力结果符合人体运动规律, 能反映人体行走过程中的肌肉协同和肌肉力相对关系, 并且可以与文献报道的结果相互印证。

5 结论

针对现有计算肌肉力的方法中缺乏对运动过程中肌肉协同机制和肌肉力相对关系考量的问题, 本文提出一种基于人体行走过程中肌肉协同和肌肉力相对关系的肌肉力计算方法, 可以通过关节力矩和运动学数据计算出行走过程的肌肉力。针对人体肌肉骨骼系统的结构和功能以及行走过程的步态特征, 本文建立了人体肌肉骨骼模型, 并划分行走步态。利用行走过程中的运动学和动力学数据, 在优化过程中考虑了肌肉力之间的关系以及肌肉协同关系, 得到符合人体运动学规律的肌肉力计算结果。

该方法有助于研究人体运动过程中的生物力学机制, 对探究人体神经运动控制原理、开发仿人机器人等也具有十分重要的作用。

参考文献

[1] Tsirakos D, Baltzopoulos V, Bartlett R.Inverse optimization: functional and physiological considerations related to the force-sharing problem.Critical Reviews™ in Biomedical Engineering, 1997, 25(4/5): 34–43

[2] Erdemir A, McLean S, Herzog W, et al.Model-based estimation of muscle forces exerted during movements.Clinical Biomechanics, 2007, 22 (2): 131–154

[3] Finni T, Komi P V, Lukkariniemi J.Achilles tendon loading during walking: application of a novel optic fiber technique.European Journal of Applied Physiology and Occupational Physiology, 1998, 77(3): 289–291

[4] Dietz V.Spinal cord pattern generators for locomotion.Clinical Neurophysiology, 2003, 114(8): 1379–1389

[5] Crowninshield R D.Use of optimization techniques to predict muscle forces.Journal of Biomechanical Engineering, 1978, 100(2): 88–92

[6] Trinler U, Hollands K, Jones R, et al.A systematic review of approaches to modelling lower limb muscle forces during gait: applicability to clinical gait analyses.Gait & Posture, 2018, 61(1): 353–361

[7] Sutherland D H, Cooper L, Daniel D.The role of the ankle plantar flexors in normal walking.The Journal of Bone and Joint Surgery: American Volume, 1980, 62(3): 354–363

[8] Bogey R A, Perry J, Gitter A J.An EMG-to-force processing approach for determining ankle muscle forces during normal human gait.IEEE Transactions on Neural Systems and Rehabilitation Engineering,2005, 13(3): 302–310

[9] Zajac F E.Muscle and tendon: properties, models,scaling, and application to biomechanics and motor control.Critical Reviews in Biomedical Engineering,1989, 17(4): 359–411

[10] Anderson K G, Behm D G.Maintenance of EMG activity and loss of force output with instability.The Journal of Strength & Conditioning Research, 2004,18(3): 637–640

[11] Davy D T, Audu M L.A dynamic optimization technique for predicting muscle forces in the swing phase of gait.Journal of Biomechanics, 1987, 20(2): 187–201

[12] De Groote F, Pipeleers G, Jonkers I, et al.A physiology based inverse dynamic analysis of human gait: potential and perspectives.Computer Methods in Biomechanics and Biomedical Engineering, 2009, 12(5):563–574

[13] 村井昭彦.Modeling and computation of human neuromusculoskeletal system and their application with visualization [D].東京: 東京大学, 2009

[14] Webb R C.Smooth muscle contraction and relaxation.Advances in Physiology Education, 2003, 27(4): 201–206

[15] Racic V, Pavic A, Brownjohn J M W.Experimental identification and analytical modelling of human walking forces: literature review.Journal of Sound and Vibration, 2009, 326(1): 1–49

[16] Mummolo C, Mangialardi L, Kim J H.Quantifying dynamic characteristics of human walking for comprehensive gait cycle.Journal of Biomechanical Engineering, 2013, 135(9): 091006

[17] Hoy M G, Zajac F E, Gordon M E.A musculoskeletal model of the human lower extremity: the effect of muscle, tendon, and moment arm on the moment-angle relationship of musculotendon actuators at the hip,knee, and ankle.Journal of Biomechanics, 1990, 23(2):157–169

[18] Trinler U, Schwameder H, Baker R, et al.Muscle force estimation in clinical gait analysis using AnyBody and OpenSim.Journal of Biomechanics, 2019, 86(6): 55–63

[19] Gomes A A, Ackermann M, Ferreira J P, et al.Muscle force distribution of the lower limbs during walking in diabetic individuals with and without polyneuropathy.Journal of NeuroEngineering and Rehabilitation, 2017,14(1): no.111

[20] Lin Y C, Dorn T W, Schache A G, et al.Comparison of different methods for estimating muscle forces in human movement.Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, 2012, 226(2): 103–112

[21] Fraysse F, Dumas R, Cheze L, et al.Comparison of global and joint-to-joint methods for estimating the hip joint load and the muscle forces during walking.Journal of Biomechanics, 2009, 42(14): 2357–2362

[22] Winter D A.The biomechanics and motor control of human gait: normal, elderly and pathological.2nd ed.Waterloo: University of Waterloo Press, 1991: 112–123

[23] Cappellini G, Ivanenko Y P, Poppele R E, et al.Motor patterns in human walking and running.Journal of Neurophysiology, 2006, 95(6): 3426–3437

[24] Barroso F O, Torricelli D, Moreno J C, et al.Shared muscle synergies in human walking and cycling.Journal of Neurophysiology, 2014, 112(8): 1984–1998

A Muscle Force Calculation Method on Lower Limb during Human Walking Based on Muscle Synergy and Muscle Force Relationships

ZHANG Ling1, HUANG Yan1,2,3,†, WANG Qining4

1.School of Mechatronical Engineering, Beijing Institute of Technology, Beijing 100081; 2.Beijing Advanced Innovation Center for Intelligent Robotics and Systems, Beijing Institute of Technology, Beijing 100081; 3.Key Laboratory of Biomimetic Robots and Systems, Ministry of Education, Beijing 100081; 4.Department of Advanced Manufacturing and Robotics, College of Engineering, Peking University, Beijing 100871; † Corresponding author, E-mail: yanhuang@bit.edu.cn

Abstract In response to the issue of existing muscle force calculation methods lacking consideration for the muscle synergy and muscle force relationships during motions, a biomechanical calculation method based on the muscle force relationships and muscle synergy during walking is proposed.This method incorporates the muscle synergy and the relative force of muscles into inverse dynamics calculation based on the quadratic optimization.The proposed method is applied for the calculation of forces of major lower-extremity muscles of human walking.Muscle forces are obtained based on a musculoskeletal model, which effectively reflects the human motion mechanisms and muscle synergy during walking.This method provides significance for exploring the laws of human movement,biomechanical principles, and humanoid robot motion control.

Key words human musculoskeletal model; inverse dynamics calculation; muscle synergy; muscle force relationship