王砚麟， 王克义， 章梓星， 马春龙,2， 莫宗骏

(1.哈尔滨工程大学 机电工程学院，黑龙江 哈尔滨 150001； 2.哈尔滨职业技术学院 汽车学院，黑龙江 哈尔滨 150001)

目前，我国有近3 000万肢体存在运动障碍的患者，且只有30%左右的人能接受到不同程度的康复训练，但康复效果各不相同。

对人体运动生物力学的研究是康复医疗技术和康复设备研发的基础，目前，下肢存在运动障碍的人口数占比较大，且对下肢运动生物力学的相关研究报告也比较多：Deschamps等[1]通过动力学链的方法对下肢包括髋关节、膝关节和踝关节的下肢运动和关节力矩的变化进行研究，主要以将健康的运动员运动参数为例进行研究，但并没有考虑下肢肌肉力的变化规律。Hume等[2]建立了人体下肢缩比模型，研究了结构负载和关节力学的预测方法，对受伤患者、术后等下肢肌肉力和力矩负载的预测具有重要实际意义。在步行过程中，速度、步长、频率和不同坡度对下肢的运动学和动力学及肌肉力均有不同程度的影响[3-10]，例如：大腿屈曲/伸展-大腿屈曲/伸展的vCRP(连续相对相位(continuous relative phase, CRP))值对人体运动速度很敏感，而且跑步的速度增加，下肢的协调性下降，相当于对下肢添加了一定的约束。另外，CRP值的减小是由于人体质心在垂直方向的运动所导致的，并非膝关节的屈曲角度[3]。Resende等[11]分析了单侧脚内旋可能造成的下肢伤害，对下肢康复和下肢矫正提供了的依据，但没有涉及下肢肌肉力的变化，无法从机理上研究下肢康复策略。Cheng等[12]分别针对健康青年、步态稳定老年人和步态不稳定的老年人在完成“坐到站”动作时肌肉功率的变化进行了分析，发现健康青年消耗的功率比老年人消耗的功率要大，且完成动作时间短。Koyama等[13]研究了在光脚和穿鞋的情况下进行向下跳跃运动的下肢力学特征和肌肉激活度，研究结果表明：在光脚和穿鞋2种情况下，下肢运动的关节角度差异较小，但在人体着地时，下肢肌肉激活度存在较大的差异。另外，在穿鞋的情况下向下跳跃是比较安全的。Carolina等[14]设计了一种混合神经的假肢用来实现下肢髋关节的康复训练，其中下肢运动关节力矩和运动数据在OpenSim软件环境中获取，基于上述数据研究了假肢的控制方法。

目前报道的文献未见单曲膝行走时下肢生物力学特性的研究，当人体下肢受伤后，所表现出来的保护性跛行的步态称之为单曲膝行走步态。因此，本文将针对青年正常行走步态和单曲膝病态行走下的运动学、动力学和下肢肌肉力及肌肉激活度进行分析，研究结果将为下肢康复医疗提供数据，以填补单曲膝下肢康复数据的空白。

1 下肢运动学、动力学模型

图1 人体下肢骨骼结构及简化模型Fig.1 Skeleton structure and simplified model of human lower limb

1.1 人体下肢骨骼运动学分析

(1)

(2)

根据角速度叠加原理可知，下肢骨骼运动绝对角速度可表示为：

(3)

式(3)可进一步表达为：

(4)

对式(4)两边微分得到人体下肢骨骼运动绝对角加速度：

(5)

1.2 人体下肢骨骼动力学分析

图2 下肢小腿骨骼受力Fig.2 Schematic diagram of the force on the shank

对图2所示的小腿骨骼受力进行分析，下肢其他骨骼受力和小腿受力一致，因此根据牛顿欧拉方程可得到人体下肢骨骼运动的动力学模型：

(6)

且关节力约束力和力矩之间存在关系：

(7)

下肢骨骼运动的动力学模型(6)可进一步表示为：

(8)

式中：τ为关节驱动力矩；X为下肢骨骼广义质量矩阵；B为Coriolis力和离心力；f、L分别表示肌肉力映射的广义力矩阵和对应的力臂矩阵；G为重力产生的广义力矩阵。

1.3 肌肉力力学分析

图3 肌肉力驱动模型Fig.3 Muscle-driven simulations model

1.3.1 肌肉激励动力学

肌肉的激活和神经信号的刺激之间存在一定滞后，因此，采用简化的一阶微分方程表示激活度和神经刺激之间的关系[15]：

(9)

(10)

(11)

式中：Tact、Tdeact分别表示肌肉激活程度的上升和下降时间常数，Tact减小，由于低效钙的释放和扩散会使激活度增加，类似地，Tdeact减小，由于供肌浆网吸收的钙离子不足会使激活度减小[16]；amin表示肌肉激活程度的最小值。为了避免肌肉腱数值计算的奇异性，需要通过修正传统的肌肉激活度状态方程[16]使激活度可以光滑变化。

1.3.2 肌肉腱力平衡方程

肌肉腱驱动包含主收缩元、被动弹性元件和弹性肌腱，如图4所示，其中被动弹性元件的力和归一化肌纤维长度之间的关系可以表示为：

(12)

图4 肌肉腱力学简化模型Fig.4 Simplified model of musculo-tendon contraction mechanics

主收缩力和归一化的肌纤维长度之间的关系表示为高斯函数的形式：

(13)

式中ν表示力-速度曲线变异系数。

肌肉力fM可以表示为[17]：

(14)

(15)

方程式(14)、(15)表示了肌肉力和肌腱力与肌肉纤维长度、速度和力之间的关系。

(16)

肌肉长度和速度与肌肉腱长度之间存在运动学关系：

lMT=lT+lMcosα

(17)

对式(17)两边分别对时间求导数，可以得到肌肉、肌腱和肌肉腱驱动速度之间的关系：

(18)

1.3.3 肌肉力的优化控制

上面建立了下肢的肌肉骨骼模型和肌肉力学模型，但还要通过约束优化解决冗余肌肉力的问题，对肌肉力的优化应用最广泛的方法为静态优化。选取式(19)的静态优化目标函数:

(19)

式中：N为肌肉数量；ai为第i块肌肉的激活度；n的选取要根据不同的要求进行选择，本文选取n=2，选取肌肉激活度平方和最小函数J作为优化目标，可使人体运动感到很好的舒适性。

图5 CMC算法框图Fig.5 Schematic of CMC algorithm applied to gait

(20)

且线性反馈控制器所计算得到的肌肉激励为：

u=a*+Ku(a*-a)

(21)

通过合理设置PD控制参数，可以实现下肢生物力学参数的最优化和下肢运动障碍患者的运动学、生物力学等参数的仿真，为下肢运动障碍患者的康复训练模式提供依据。

2 人体下肢运动数据采集

为了分析和比较单屈膝行走步态与正常行走步态的下肢运动生物特性差异及对主要肌肉的疲劳损伤的影响，需要对2种行走步态的运动学进行测量(在无任何外界干扰的情况下)，将测量髋关节、膝关节和踝关节的运动数据进行函数拟合，并将拟合函数作为在OpenSim建立的3-D人体肌骨模型的运动驱动，最后通过仿真2种行走步态，对下肢肌肉收缩量等问题进行分析。

本文以8名男性青年 (身高165±10 cm，体重60±5 kg，年龄24±1岁)作为受试对象，通过FAB采集所需的人体运动数据，FAB系统的惯性传感器结合了加速度传感器、方位传感器和足底力传感器，将各部位的惯性传感器安装在受试者身体相应的位置上，受试者分别按2种步态行走，FAB会通过实时无线步态检测系统的无线步态数据接收端检测和接收人体运动相关的数据，步态数据接收端会通过USB接口将数据传输到计算机分析软件中，软件就可以实时记录人体运动相关的所有运动学和动力学数据，且采样频率为100 Hz。传感器在人体上的安装和FAB采集系统如图6所示。

图6 传感器安装和FAB采集系统Fig.6 Installation positions of sensors and FAB collection system

受试者行走100 s，在单曲膝行走步态行走过程中，右腿受伤屈膝行走，左腿正常，因此，只研究右腿的运动生物力学特性，忽略膝关节和踝关节的内翻/外翻和内旋/外旋的现象，对记录的下肢运动参数以步态周期为步长取运动参数的平均值，并利用傅里叶级数对关节运动数据进行拟合，傅里叶级数的拟合函数表达式为：

(22)

式中：θ(t)为关节运动角度；t为时间；a0、ai、bi和ω均为傅里叶级数方程的系数，n∈Z+的大小取决于拟合度的大小，本研究中拟合度要求不低于99.4%。

针对2种行走步态在步态周期内的关节角度由FAB系统采集记录，通过傅里叶级数进行拟合，其傅里叶级数拟合函数的参数如表1所示。

表1 傅里叶级数拟合函数的参数Table 1 The parameters of Fourier series fitted function

由表1中参数构成函数的拟合度均大于99.4%，FAB采集到的下肢关节运动角度和傅里叶级数拟合角度的变化曲线如图7所示，结果表明拟合函数具有很高的跟踪性能，可以直接作为3-D人体肌骨模型关节运动的驱动函数。其中NWG为正常行走步态数据，SCWG为单曲膝行走步态数据。

图7 下肢关节角度变化曲线Fig.7 Joint angle change curves of lower extremity

由图7可知：人体正常步态行走时，髋关节的屈伸角变化范围为[-21.046°,21.534°]，内收外展角变化范围为[-6.9.5°,6.991°]，内旋外旋角变化范围为[-10.572°,6.325°]，膝关节变化角度变化范围为[-69.27°,0.647°] ，踝关节的屈伸角度变化范围为[-6.85°,16.05°]，角度变化速度较大，且髋关节的屈伸角和膝关节角度的变化范围很大；而在单屈膝步态行走时，髋关节的屈伸角变化范围[9.891°,45.812°]，内收外展角变化范围为[11.935°,14.956°]，内旋外旋角变化范围为[16.335°,20.089°]，膝关节变化角度变化范围为[-47.223°,-28.911°]，踝关节的屈伸角度变化范围为[-35.595°,-3.391°]，角度变化速度较小，因此，在下肢单曲膝步态行走时，关节角度变化范围较小，变化速度较小。从变化规律来看，髋关节内旋/外旋角的变化规律和正常步态由较大差别，其他关节角度变化规律和正常步态基本一致。

3 单曲膝行走步态对下肢肌肉的影响

根据受试者身高、体重，在OpenSim软件中进行SCALE人体3-D模型设置，以图7所示的关节角度变化数据驱动人体下下肢运动，对下肢运动的动力学模型进行仿真，图8所示为下肢在NWG 和SCWG步态下人体3-D模型的周期运动模型，由图可知：在SCWG行走时，2条腿的运动学特征均有变化，但受伤的右腿变化更为明显。图9所示为下肢关节驱动力矩的变化曲线，其中实线表示正常步态，点划线表示单曲膝行走步态。

图8 NWG 和SCWG的人体3-D运动模型对比Fig.8 Comparison of human 3-D movement models of NWG and SCWG

下肢各关节运动对应驱动力矩如图9所示，和正常行走步态相比，按SCWG行走时，下肢各关节运动驱动力矩波动较大，且在步态周期的48.7%～53.6%，髋关节屈伸运动力矩、髋关节内收外展运动力矩和膝关节运动力矩出现的波动很大，波动范围分别为[-197, 211]N·m、[-149, 161]N·m和[-80, 75] N·m。在步态周期的0～15.9%之间，SCWG的下肢各关节的驱动力矩明显大于正常行走步态的驱动力矩(髋关节屈伸力矩除外)，在步态周期的78.2%～100%，SCWG的下肢膝关节和踝关节的驱动力矩明显小于正常行走步态的驱动力矩。因此，当人体下肢受伤后以SCWG行走时，在步态周期的48.7%～53.6%和0～15.9%会产生明显的不舒适感，对该类患者的康复治疗和训练过程中需要考虑这点。

图9 下肢关节驱动力矩Fig.9 Joint torque change curves of lower extremity

在CMC计算，参数设置为Kp=100，Kv=20和Ku=1。在行走步态周期中，下肢运动主要的12块肌肉力变化如图10所示，其中实线表示正常步态，点划线表示单曲膝行走步态。总的来看，SCWG行走时，半膜肌、半腱肌、股二头肌、长收肌、臀大肌、髂肌、腰肌、股四头肌、股直肌和胫骨前肌的肌肉力明显大于NWG下的肌肉力，而比目鱼肌和胫骨后肌的肌肉力小于NWG下的肌肉力，SCWG行走的12块肌肉的肌肉力的波动幅度均大于NWG的肌肉力的波动幅度，主要是由于步态的不稳定引起的。这表明：SCWG行走时，半膜肌、半腱肌、股二头肌、长收肌、臀大肌、髂肌、腰肌、股四头肌、股直肌和胫骨前肌均起到了主要作用，而比目鱼肌和胫骨后肌起到辅助的作用。

图10 下肢肌肉力变化曲线Fig.10 The change curves of muscle forces of lower extremity

和NWG相比，SCWG在步态周期后期86.1%～100%，半膜肌、半腱肌和股二头肌的肌肉力减小大约一半，在步态后期63.8%～100%，胫骨前肌的肌肉力为NWG的2倍，在步态周期前期0～14%，长收肌、臀大肌、髂肌、腰肌和股四头肌的肌肉力较小。

图11给出了下肢运动过程中肌肉激活度的变化情况，其中实线表示正常步态，点划线表示单曲膝行走步态，由图可知肌肉激活度的变化情况和肌肉力的变化趋势基本一致，激活度的变化情况此处不再赘述。

图11 下肢肌肉激活度Fig.11 The change curves of muscle activation of lower extremity

在SCWG行走时，下肢关节力矩的增加和大幅度波动，会导致关节运动不稳定，增加关节的损伤和肌肉的负担，肌肉活性的增强同样会导致肌肉的疲劳损伤等问题的出现，因此，研究结果为下肢损伤的患者康复和康复设备的研发提供了数据基础。

4 结论

1)在单曲膝行走步态下，下肢各关节角度的变化范围小，特别是髋关节和膝关节的变化范围尤为明显。

2)在步态周期的初期，单曲膝行走步态的下肢关节驱动力矩的波动明显增加，并且下肢关节驱动力矩明显大于正常行走步态下的关节驱动力矩(髋关节屈伸运动的驱动力矩除外)；在步态周期的后期，膝关节和踝关节的驱动扭矩明显小于正常行走步态下的关节驱动力矩。

3)在单曲膝行走步态下，半膜肌、半腱肌、股二头肌、前胫骨肌、臀大肌、股四头肌，股直肌和胫骨的肌肉力和激活度明显增加，其波动值和频率也明显增加。