基于OpenSim人体步行腓肠肌静态生物力学分析

赵 闯， 蔡玉强

(华北理工大学机械工程学院，唐山 063000)

近年来，在伤病预防以及人体下肢康复机器人的设计中需要考虑人体下肢的生物力学状态[1-2]。在研究人体下肢康复机器人的模型构建与动力学分析中，下肢肌肉的运动状态是进行仿真的第一环节。然而受到人体生物系统的复杂性、外界环境的干扰以及伦理等的限制，现有的数据采集方式难以良好地获取腿部肌肉受力的情况[3]。

OpenSim是Stanford University开发的一个开源软件平台[4]，用于创建、更改和分析骨骼肌肉系统的计算机模型以及运动的动态仿真。文献[3]基于OpenSim软件搭建了人体肌肉骨骼仿真平台，验证了人体在直立状态下会根据外界环境的扰动选择性的使人体保持平衡。文献[5]基于OpenSim建立肌肉骨骼模型，提出了计算肌肉拉力和膝关节刚度的方法。文献[6]构建简化人体下肢力学模型，建立相应动力学方程，经过实验数据分析得到人体下肢膝关节与踝关节的受力情况，但在模型受力曲线上会出现数据不连续的情况。文献[7]基于ADAMS 动力学仿真软件中建立人-机耦合系统模型并实现该系统的运动仿真，得到特定人体按照不同足部轨迹运动时下肢肌肉伸缩量的变化情况。文献[8]研究人体在运动过程中，肌肉长度变化、肌肉长度变化速率、肌肉力臂等肌肉功能参数以及肌肉功能的动态变化。

总之，中外有不少学者对人体下肢关节运动的形态特征进行研究和分析，但对于下肢运动状态下的静态力学分析等方面研究不足。本文首先分析人体下肢受力机理，简述仿真所需骨骼肌肉模型，利用人工合成的运动学方法依赖于生物力学知识，并结合正向和逆向运动学知识[9]，使用OpenSim的Static Optimization Tool模块，简化腰部肌肉模型为扭矩执行器，改善肌肉激活度和静态优化输出的力，利用逆向运动学(IK)、残差缩减(RRA)和肌肉计算控制(CMC)对人体下肢肌肉进行静态优化分析，为人体下肢康复机器人的研究与研发提供基础。

1 人体肌肉骨骼模型

1.1 人体下肢受力

人体的膝关节和踝关节与小腿部肌肉有着复杂的肌骨结构，涉及肌肉、韧带、骨骼、神经和血管等，是一个复杂的神经运动系统。人体通过前庭、视觉、感受和触觉的综合分析处理，控制着肌肉完成身体运动的平衡[5]。因此，在研究人体行走过程中下肢的受力状态，不仅要考虑足底受力和其他质量，还要考虑关节肌肉的受力[10-11]。

图1 人体下肢骨骼-肌肉模型Fig.1 Human lower limb bone-muscle model

人体站立、行走和慢跑等保持身体平衡的运动小腿部分主要是通过腓肠肌、比目鱼肌和胫骨前肌之间的相互调节作用来实现的，如图1人体下肢骨骼-肌肉模型。腓肠肌的作用是曲踝关节和曲膝关节，于皮肤表面，其深方是比目鱼肌。比目鱼肌起于胫骨和腓骨头部，尾下为跟腱，止于跟底，受胫神经支配。在运动过程中，腓肠肌和比目鱼肌作为一体，维持身体运动过程中小腿部肌肉的肌力，在站姿或膝关节伸直、蹬足时发力的作用。胫骨前肌主要参与膝关节和踝关节的背曲运动，与腓肠肌和比目鱼肌协调作用，维持人体平衡。

1.2 肌肉反射模型

将人体的下肢运动受力模型简化基于图2骨骼-肌肉模型简图，该模型适用于腿部骨骼-肌肉近似的平面运动，上肢同时保持近似垂直的状态。其中人体上半身简化为固定在臀部的质量块；大腿、小腿和脚简化为3个连杆，髋关节、膝关节和踝关节简化为铰接点；右腿肌肉简化为2对作用在连杆两端的相对作用力，分布于大腿和小腿两侧。

图2 骨骼-肌肉右腿模型简图Fig.2 Skeleton-muscle right leg model diagram

由于肌肉运动的非线性，希尔模型[12-13]包含主动收缩单元(CE)和被动弹簧单元(SE)表示的肌肉模型，用非线性弹簧表示肌腱，如图3所示。lMT为肌肉-肌腱长度，lT为肌腱长度，lM为标准化的肌纤维长度，α为羽状角，它们之间存在lMT=lT+lMcosα的关系。肌纤维收缩产生的拉力通过肌腱传递给所附着的骨骼，引起关节运动，力量的传递与α有直接关系，α越小传递越明显，反之亦然；但是，一些主要肌肉为羽状肌，主因在于相同体积羽状肌的生理横截面积更大[14]。

图3 希尔模型Fig.3 Hill model

研究单步步态试验中执行静态优化分析，并在简化的示例模型中评估腓肠肌激活度和力量产生。静态优化结果的质量在很大程度上取决于输入：模型、运动和力等参数。静态优化输入为肌纤维长度的变化值v(t)，肌肉反射控制模型的表达式为

R(t)=Kl[l(t)-l0]+Kvv(t)

(1)

v(t)=[l(tn+1)-l(tn)]/(tn+1-tn)

(2)

式中：R(t) 为t时刻肌肉激活度控制值；Kl为肌肉纤维长度变化系数；l0为静态下肌腱长度；Kv为肌肉纤维速度变化系数；v(t) 为肌肉纤维长度变化速率；l(t)为t时刻肌肉纤维长度变化量；tn为步态的时间。

为了建立小腿部腓肠肌的控制方程，对脚部进行受力分析，几何标注如图4所示。

图4 脚部受力示意图Fig.4 Schematic diagram of the force on the foot

对于脚部质心A分别建立x、y方向上的移动分量[xA(t)、yA(t)]和转动分量[αy(t)]的3个运动方程。

(3)

(4)

(5)

式中：mf为脚部质量；Fax、Fay为踝关节x、y方向的分力；F1x、F1y为腓肠肌x、y方向的分力；F2x、F2y为胫骨前肌x、y方向的分力；Fgx、Fgy为脚部与地面之间摩擦力分力；Gf为脚部自身的重力；Gt为单脚承受人体50%的重力；Jf为脚部转动惯量；Max、May为踝关节x、y方向的力矩；M1、M2为腓肠肌和胫骨前肌力矩；Mgx、Mgy为x、y方向摩擦力矩。

其中脚部质量mf和转动惯量Jf等生物力学参数可由经验方程[15]得到，脚部压力可通过仪器得到[16]。式(3)～式(5)可以计算踝关节在各矢量面上Fa的分力，并且可以推导出腓肠肌受力F1的表达式：

(6)

将得到的数据代入式(6)，即得到腓肠肌受力与时间之间的对应关系。

1.3 生物力学参数确定

由于难以获取脚部的生物力学数据，在模型中将质心作为脚部几何中心，并将看作刚体，脚部质心位置lf由人体脚部大小确定，转动惯量Jf、脚部重量Gt和人体重量Gf由人体解剖学参数确定。

2 模型动态仿真试验

OpenSim作为建立和分析骨骼-肌肉的开源软件，可通过软件建立人体骨骼-肌肉模型，设置人体参数来计算和利用OpenSim的plot功能绘制肌肉力、膝、踝关节的运动位移曲线和力矩等数据，进行人体下肢步行状态下肌肉静态力学分析。

2.1 人体肌肉骨骼搭建

研究采用的骨骼-肌肉模型如图5所示，具有10个自由度和18个肌肉的平面缩小步态模型，模拟正常人体步行运动状态。由于研究不涉及上肢运动，经简化保留躯干部分，并在分析过程中简化为质量块处理。模型下肢经简化后保留支撑人体站立的胫骨前肌、腓肠肌和比目鱼肌[16]。

图5 人体骨骼-肌肉模型Fig.5 Human bone-muscle model

利用Notepad++文本编辑器软件修改原有模型的内部参数，使该新模型subject01.osim具有代表实验参与者的质量、人体身高尺寸和强度；自由度和肌肉几何形状适合于本次研究的肌肉群——腓肠肌；在腰部添加扭矩执行器，保证人体上肢有足够的驱动控制；在腓肠肌添加备用执行器，肌肉不能产生足够的加速度时可以为步态周期的部分期间增加额外的动力；添加残差执行器，测量分析程中运动和力之间的关系，确保整个分析满足牛顿第二定律。模拟运动包含平滑、真实的加速度，并能精确地测量运动期间的肌肉的激活度与力矩。图6为腓肠肌受力曲线。

图6 腓肠肌受力大小Fig.6 The size of the gastrocnemius muscle

2.2 肌肉静态激活度分析

人体在运动状态下，调动休眠的肌肉参与活动，这些被调动的肌肉称为激活。参与活动的肌肉群越多，其力量和感觉等方面比较高，激活度就越高，反之则激活度越低。静态优化的产生需要输入运动学中加速度的肌肉力，如果这些加速度不稳定，则肌肉激活和力量会产生波动。研究两种改进运动数据的方法：运动学反解(IK)和使用残差缩减算法(RRA)来平滑运动曲线。

过滤运动学将对运动数据进行样条拟合并过滤坐标位置，从而使加速度趋于稳定。如果使用IK运算结果作为静态优化的输入，则应始终在外部使用MATLAB或Python或使用OpenSim过滤器进行过滤，本次分析使用OpenSim进行过滤。

腓肠肌过滤加速度使用运动学反解和残差缩减算法，如图7所示。在人体步行时，IK从0.85 s开始激活度从0.01增大到0.37，总时长用了0.35 s；RRA从0.9 s开始，激活度从0.01增大到0.34，总时长相对IK时间缩短0.05 s。即利用RRA方法更能模拟人体肌肉激活度，又由于RRA使用逆动力学模拟，因此输出运动学将是一致的加速度，非常适合在静态优化中使用。

图7 腓肠肌不同优化方法激活度Fig.7 Different optimization methods for the gastrocnemius

图8 腓肠肌Fx，Fy方向上残差Fig.8 The residual of gastrocnemius in Fx, Fy direction

从图8中可以看出，利用IK法在x和y分量上的残差区间分别为[-44.9,22.0]、[-55.6,53.5]；利用RRA法在x和y分量上的残差区间分别为[-0.6,8.1]、[-14.0,21.5]；即在x分量上RRA比IK的波动范围缩小了近9倍，在y分量缩短了近3倍，且曲线更加平滑。通过对比，可以明显地看出利用RRA法比IK法提高了静态分析准确度。

2.3 调整计算控制分析

优化器将优先使用肌肉来产生关节扭矩并且仅在肌肉不能产生足够扭矩时才使用备用执行器。降低肌肉控制计算(CMC)的强度，并观察备用执行器扭矩和腓肠肌激活度和力的最终变化。对于2.2节分析过程中使用备用执行器为100 N的最佳力，并且对于每个扭矩执行器具有100 N·m的最佳扭矩。利用Notepad++将编辑执行器文件的每个执行器的最佳力从100 N更改为1 N。

从图9中可以看出腓肠肌在最佳执行力为100 N时，激活都峰值为0.34，此时模拟腓肠肌未被完全激活进行分析，备用执行器被调用，导致分析结果可能出现误差；将最佳执行力降到1 N时，此时腓肠肌激活度峰值为0.82，但相对于肌肉而言，成本较高。对比图10和图8(b)，可以看出最佳执行力的改变对于x和y分量上的残差影响并不明显。

图9 腓肠肌不同优化强度激活度Fig.9 Different optimal intensity activation of gastrocnemius

图10 CMC最小执行力Fig.10 CMC minimum execution

3 结论

利用OpenSim分析人体在行走过程中腿部肌肉静态力学分析，以腓肠肌为例。通过对比IK和RRA分析方法得出如下结论。

(1)RRA法提高了静态分析准确度，适用于之后人体下肢康复机器人的设计分析。同时分析出激活度a的变化曲线，防止康复机器人设计激发激活度a变化过快，导致人体肌肉痉挛。

(2)降低外设备的辅助制动力，最大化提升自身肌肉的激活度，可以为进一步揭示膝关节和踝关节的生物力学特性和人体下肢康复机器人的设计奠定基础。