人体运动仿真技术是基于生物力学、计算机科学和机器人学建立的科学研究方法。通过建立人体模型，经过动力学计算或相关优化计算法，得到人体完成动作过程中相关肌肉、关节的受力，以及这些组织与运动过程中运动学指标之间的关系，也可以进行运动中神经肌肉系统对动作协调控制机理等问题的研究。目前，人体运动仿真方法已经被广泛的应用于生物医学、体育科学、航空航天等多个领域。
跟腱是人体最强壮、最厚实的肌腱，它连接比目鱼肌和腓肠肌到跟骨，确保踝关节的跖屈。由于这种生物力学性能，跟腱与人体活动的许多能力有着密切的关联。研究表明，竞技运动中，由于运动项目的性质、强度、场地、环境等因素的影响，跟腱损伤的发病率在不断提升。肌腱的生物力学性能在一定程度上影响着肌肉的收缩力和运动成绩，在竞技体育和体育锻炼中，经常发生的肌腱损伤（如肌腱炎、肌腱断裂等）与肌腱的生物力学性能密切相关。因此，对跟腱力学性能的研究对伤病的预防和治疗有着深远意义。
由于跟腱特殊的生理解剖结构，以往的研究方法存在很多局限性，大部分实验研究的对象也集中在动物上，活体实验较少。研究方法从植入性力学传感器到超声波影像都曾被应用，但都属于探索阶段，且由于实验仪器设备本身的局限性，对于实际体育运动的研究还比较少。
随着生物力学仿真技术的发展，基于骨骼肌肉模型的运动仿真方法为探索跟腱损伤力学机理提供了方法。本研究通过医学图像建立骨骼几何学模型，并采用基于动态优化的方法来计算下肢肌肉的即时受力。结合运动追踪技术和肌肉骨骼模型的生物力学仿真方法来比较跟腱在跨栏和跑步时，肌肉-跟腱单位的拉伸以及跟腱的负荷，分析跨栏跟腱伤病产生的力学机理，探讨该研究方法在分析跟腱力学性能方面的应用，从而为跟腱运动伤病的预防提供科学依据，为进一步跟腱损伤研究奠定基础。
1、研究对象与方法
1.1 研究对象
上海体育学院运动训练专业跨栏专项男性运动员 10 名，其中，2 名国家一级运动员，8 名国家二级运动员。受试者年龄（20.67±1.53）岁，身高（1.83±0.04）m，体重（70.33±2.52）kg，110 m 栏最好成绩 14.5~15 s。所有受试者右腿为惯性起跨腿，且半年内无任何下肢伤病。
1.2 试验数据采集
1.2.1 动作定义 起跨动作是指跨栏跑中起跨腿扒地蹬伸的一步，从起跨腿脚着地准备起跨，到起跨腿经着地缓冲到蹬伸离地结束。本研究将一个完整的起跨动作定义为一个支撑周期（即从受试者起跨脚接触测力台开始到完全离开测力台结束）。
1.2.2 试验流程 每名受试者在跑步机上进行 30 min 的充分热身之后，进行跨栏练习进而熟悉场地。热身完成后，对受试者全身解剖学位置粘贴 Maker 球用来帮助捕捉运动学数据（见图1）。
运动学数据采集之前，每名受试者需要采集一组静态数据，用来确定人体关节解剖学位置和计算关节中心。在静态数据采集结束之后，移除受试者膝关节和踝关节的马克点，然后进行跨栏组和短跑组的数据采集。跨栏组数据，要求受试者在高速助跑后跨越一个栏架；短跑组数据，要求受试者高速完成相同距离的平跑。每名受试者需采集 5 组跨栏数据和 5 组短跑数据，1 组成功数据的定义为受试者右脚完全踩在测力台上，全身的 Marker 清晰可见。
采用 Vicon16 个摄像头红外高速摄影系统采集运动学数据（Vicon Oxford Metrics，Oxford，简称 UK），采集频率为 200 Hz。使用 1 块 Kistler 测力台同步采集地面反作用力数据（KistlerCorporation，Ohio，USA），采集频率为 1 000 Hz。
1.2.3 试验数据处理 使用 Vicon 的 Workstation 和 Bodybuilder软件处理运动学和动力学数据；使用 Butterworth 数字滤波器对Marker 球的轨迹、地面反作用力进行低通滤波，滤波频率分别为 15 Hz 和 55 Hz；使用 OpenSim 仿真软件完成整个仿真处理过程；利用 Excel2007 进行试验结果的图表绘制。
1.3 OpenSim 环境下建模
传统仿真方法研究中，由于所涉及人体运动的复杂性，通常将人体简化为多刚体系统，把人体的肌肉、筋腱等组织处理为各刚体间的作用力及力矩。应用 OpenSim 软件的仿真研究，是通过建立人体骨骼肌肉模型计算肌肉力大小，模型中的人体基本参数根据人体实际测量与仿真计算相结合的方法得到。在进行试验研究时，会根据试验采集的人体静态文件进行相应的计算，即将通用模型调整为符合试验对象人体参数的模型。每个受试者都有专门的模型，从而保证试验结果的精确性。
利用 OpenSim 软件创建三维下肢骨骼肌肉模型（见图 2），用来计算肌肉受力。这个模型由 13 个节段、12 个环节和 23个自由度机械链接组成，包括 54 块肌肉-肌腱单位，可以在矢状面、冠状面和水平面运动。模型的头部、手臂和躯干作为一个刚体结构，相对于骨盆有 3 个旋转的自由度。骨盆在 3 个维度上可以旋转和移动，髋关节是一个球窝接头，膝关节和跖骨关节为铰链结构。踝-距骨关节有 2 个与解剖学轴线平行的关节，每块肌肉的几何学数据由解剖学模型决定。肌肉参数包括羽状角、优化纤维长度、肌腱松弛长度、肌腱应力-长度曲线和最大等速肌力。
跟腱肌肉的生理学参数见表 1。
对于 54 块肌肉-肌腱结构中的每个单位，Hill 肌肉-肌腱模型用来表示动员情况和肌肉-肌腱的收缩力学。Hill 模型包括可收缩的元件，一系列弹性元件和平行弹性元件。收缩元件表示模型中主动产生力，非线性的平行弹性元件和弹性元件属于被动元件，平行弹性元件表示组织支撑和肌肉纤维的连接。平行弹性元件表示肌纤维和骨骼肌腱的连接，在肌肉-肌腱模型中，基于肌肉肌腱的长度和动员程度作为输入参数，然后通过方程式计算的肌纤维长度来描述动员程度和收缩动力学。通过肌肉的动员-受力-长度-速度特征以及肌腱的弹性特征来分配互动的比例。当肌纤维长度确定之后，根据高斯函数计算出表示动员受力-长度关系的肌肉受力。
下肢和背部关节由 54 块肌肉-肌腱的 Hill 模型驱动器驱动，手臂由力矩驱动器驱动。基于每名受试者人体测量学数据，对通用模型按照比例进行缩放，受试者人体测量学数据依据Maker 球的解剖学位置进行计算调整。创建的模型会显示出相应 Marker 球的位置，利用逆动力学算法计算出关节角度数据，从而减小试验中每一帧影像上粘贴的 Marker 球位置与之相应的实际解剖学位置的差异。关节力矩通过 RRA（residualreduction algorithm） 算法计算得到，RRA 允许对关节角度（<1.5°）和躯干质量中点（<5 cm）进行小的改变来减少作用于骨盆的残余力。肌肉的激活情况、动员程度、肌肉力大小决定了这些力矩的结果。Opensim 中，采用 CMC（computed muscle control）算法计算肌肉力。
1.4 生物力学仿真过程
（1）通过静态数据中 Marker 球位置的输入来计算人体测量参数，基于这些人体参数建立下肢生物力学模型。模型包括所有下肢骨的肌肉-骨骼几何学参数，使模型依比例与受试者的人体测量参数匹配，匹配完成之后，得到一个特定受试者骨骼-肌肉模型。该模型包括依比例的躯干、骨盆、胫骨、股骨、腓骨和足的结构，以及依比例的肌肉-肌腱单位的几何学参数（最佳纤维长度和肌腱松弛长度）。利用逆向运动学计算肌肉肌腱长度和关节的运动学数据，如跨栏起跨过程中的关节角度和关节的位移。根据运动学数据、模型的人体测量参数值以及地面反作用力（通过测力台获得），利用这个动态系统中的运动方程计算各个关节的力和力矩。（2）通过动态优化的方法进一步计算关节力矩和支撑期的肌肉力。动力学方程计算分为 2 大部分，基于向前动力学方程和逆动力学方程 2 种，OpenSim 软件的计算基于逆动力学方程，根据试验中采集到的运动学、动力学数据得到关节角度、力矩结果。（3）关节角度和地面反作用力的数据进行 RRA 计算，优化仿真结果，最后通过 CMC 计算，得到向前动力学方程仿真结果（见图 3）。
优化过程中，肌肉力是基于给定的价值函数和生理限制（见表 1）计算得到，这个优化控制使肌肉激发的平方和最小化，可以概括为非线性梯度优化方法；肌肉激发和肌肉肌腱长度接替肌肉力学来计算上述单个肌肉力。在这一程序中，每个身体环节被看成是一个刚体，所有环节以运动链的形式链接起来。基于给定受试者身高、体重、环节长度和惯性参数，通过先前的运动学分析过程得到重心线性加速度、角速度和角加速度，利用牛顿运动学方程对环节系统进行计算。
式中：I 为环节惯性矩阵；a 为重心加速度；ω 为角速度；FD和 TD为远端关节的作用力和力矩；Fp和 Tp是近端关节的关节作用力和力矩；L 是 Fp到重心的力臂；d 是 FD到重心的力臂；W 是重力。方程 3 表明，所有作用于某一系统的外力对系统的作用等于系统动量的变化率，而线性和角速度的变化率可通过运动分析数据得到。作为对逆向动力学的延伸，通过关节力矩数据计算得到每一块肌肉的力量。这是通过一个非线性优化方程所得到的，该优化过程在一个指定的方程和物理限制下进行。
式中：α 为肌肉活动量，计算在 n 个数据桢中所有的肌肉激活的平方和；n 为所有数据点的个数。优化控制器为一个广义的既约控制非线性优化方法，可将 α 最小化。然后，肌肉的激活量、肌肉肌腱长度等数据被输入此力学模型中计算肌肉力量。
1.5 仿真结果准确性验证
整个仿真过程结束后，验证仿真结果的准确性。将 Opensim计算出的运动学（关节角度）和动力学（力矩）结果，同利用三维分析软件（Vicon 系统中的 Bodybuilder 软件） 计算的结果相对比，基本上一致，表明，该模型和仿真结果的可靠性。
该方法计算出的跟腱受力与 LICHTWARK 和 WILSON的研究结果相一致。LICHTWARK 采用运动数据和超声影像数据相结合的方法，使用一种简单的 Hill 模型来计算跟腱的受力。结果显示，跑步中跟腱受力的峰值约为 3 300 N，与本研究计算出的短跑跟腱受力相近。除此之外，本研究仿真结果中，短跑运动腓肠肌与比目鱼肌肌肉动员情况的结果图轮廓上与之前研究的表面肌电信号结果相同。
2、结果与分析
2.1 跨栏起跨动作仿真结果及其准确性验证
运动员在完成跨栏起跨动作时，起跨腿在高速跑动中完成3 个阶段，即着地缓冲、支撑和蹬伸离地（见图 4）。整个仿真结果流畅、完整，同真实运动学影像无明显差异。
仿真结果的准确性验证，是进行仿真研究的基础，也是关键性问题。关于计算机仿真结果准确性验证的方法无统一规定，目前并没有直接的方法检验仿真结果中的肌肉力大小，对于仿真准确性的验证通常采取 2 种间接的方法。（1）将仿真结果与基准数据对比验证。具体是指，相邻肌肉力矩应该等于净力矩，将仿真结果中的关节力矩与试验数据经软件处理后得到的关节力矩相比较。本研究中，Bodybuilder 计算出的关节力矩同 OpenSim得到的力矩无明显差异。（2）基于表面肌电信号技术。在仿真方法中计算肌肉力，取决于肌肉的动员程度，而表面肌电技术可以只获取单个肌肉的动员情况，因而，通过肌肉动员程度验证是多数研究采用的方法，具有较高的准确性。本研究主要涉及的比目鱼肌和腓肠肌的动员情况与之前研究中的数据结果相比对，图形轮廓上高度相似。通过与基准数据和表面肌电信号数据相比较的验证方法可以得出，本研究仿真结果真实可靠，有研究和参考价值。
2.2 肌肉-肌腱单位拉伸的长度变化
跟腱是粘弹性组织，具有粘弹性特质。在人体运动中，跟腱能够承受较强的张力将肌肉收缩产生的力传递至根骨，带动踝关节运动；同时，跟腱也具有组织的柔软特性，能够围绕骨骼的外缘改变肌肉的拉伸方向。正是由于跟腱组织的这些机械特性，在运动中能够承受很大的张力以防止过度拉伸产生损伤。
肌肉-肌腱单位组织在运动支撑期时，承受较大的牵拉力而产生很大形变。（1）腓肠肌肌肉-跟腱单位的形变较大。跨栏运动中，从蹬伸离地时最小的 0.451 mm，到支撑阶段最大拉伸时的 0.482 mm，长度差异为 0.031；短跑中，最大值与最小值分别为 0.468 mm、0.452 mm，长度差异为 0.016；在整个支撑期，腓肠肌肌肉-跟腱单位的拉伸差异主要集中在着地缓冲与支撑阶段，在蹬地阶段开始后，差异性逐渐减少，到趋于离地阶段没有差异。（2）在跨栏和短跑的支撑期，比目鱼肌肌肉-跟腱单位拉伸形变的差异都比较显著。跨栏时的长度峰值分别为 0.325 mm和 0.294 mm，长度差异 0.031；短跑时的长度峰值为 0.318 mm和 0.289 mm，长度差异 0.021。运动支撑期的拉伸程度最大峰值均出现在着地缓冲到支撑阶段（见图 5）。
2.3 跟腱负荷的仿真
跨栏与短跑支撑期，小腿三头肌（比目鱼肌和腓肠肌）牵拉产生较大拉力，拉力传递至跟腱，导致跟腱承受较大的张力。跨栏支撑期比目鱼肌肌肉力峰值（3 130.90 N）出现在支撑后的蹬伸初始阶段，而短跑中的比目鱼肌肌肉力在支撑阶段到达峰值（2 535.48 N）。在着地缓冲阶段，跨栏与短跑中比目鱼肌肌肉力都呈逐步增大趋势，但并未显示出明显的差异，随着身体重心的前移进入支撑阶段，2 种运动形式下比目鱼肌肌肉力出现不同。同为小腿三头肌的腓肠肌肌肉力峰值出现的时间与比目鱼肌同步，峰值分别为 750.91 N 和 759.182 4 N，腓肠肌肌肉力在短跑支撑期的峰值大于跨栏支撑期。跨栏支撑期跟腱所受张力在着地缓冲后不断上升，蹬伸阶段到达峰值（3 850.40 N）；短跑支撑期，跨栏所受张力在支撑阶段到达峰值（3 253.23 N）（见图 6）。
仿真结果显示，在跨栏和短跑运动中，连接跟腱的 2 块肌肉（腓肠肌和比目鱼肌）产生了较大的肌肉力，从而导致跟腱的较大负荷，这是跟腱的解剖学结构决定的。跟腱的主要作用是将肌肉产生的力传递至跟骨从而带动关节运动。肌肉力的数据结果图形显示，在支撑期的后半段，腓肠肌和比目鱼肌出现了峰值。
这是因为，在这个阶段，由于身体重心在垂直和水平方向的移动需要踝关节的跖曲来完成，而踝关节的跖曲肌群（腓肠肌和比目鱼肌）是这种向前运动的主要贡献者。研究显示，踝关节的快速跖屈和膝关节的屈伴随着跟腱的超量负荷，很容易导致跟腱的断裂。因而，与短跑运动相比，跟腱在跨栏起跨阶段更容易出现运动损伤。
研究显示，跟腱损伤的机制是在正常或过度负荷时受到张力的影响，当张力过大导致受伤时，受伤程度便视其张力的速率和力度的大小而定。跟腱是一个节能和储能的弹性结构，由于它自身的机械特性，可以承受巨大的牵拉力，所以很难想象一次简单的拉伸就可以使跟腱断裂。因此，出现了另外一种理论来解释这个现象。GALLOWAY 等认为，跟腱的过度负荷导致过度使用伤病，进而导致跟腱的断裂，而且“过度使用损伤是跟腱断裂的诱因”这一理论得到了许多验证。KANNUS 等的研究发现，97%断裂的跟腱显示出退化疾病的特征，从而进一步支持过度负荷与跟腱断裂的关系。本研究结果支持该结论，短跨运动员跟腱伤病率较高，与跟腱在从事此类运动时高负荷下的过度使用密切相关。试验结果显示，跟腱在起跨支撑期承受着巨大负荷，给运动员长期训练带来损伤隐患。
2.4 仿真方法优劣性分析
跟腱是人体最强壮的肌腱之一，与人体的运动密切相关，由于其特殊的生理解剖构成，对于跟腱负荷的研究受到诸多方面的技术限制。
VAN DEN BOGERT 等的研究是基于一种数字化模型，这种模型的优点在于简单易实现，但是这种方法的局限在于忽略了跟腱本身特殊的结构特征和机械特性，研究结果的准确性和可靠性有待加强。此外，这种建立在数字计算的模型，在人体运动研究应用中对于不同研究对象缺乏针对性，从而造成研究结果的不准确，这种方法目前已经不再应用于跟腱的研究。
LEWIS和 KOMI等的研究是在受试者机体植入力学传感器来计算跟腱负荷大小。利用外科手术的方法，将一种扣式传感器植入人体小腿跟腱处，经过一段恢复期，受试者适应传感器后开始试验，对人体基本动作（包括步态、慢跑、跳跃）的跟腱受力进行研究，应用在尸体试验得到的跟腱伸长负荷曲线来计算跟腱在上述运动中的应力曲线和应变曲线。KOMI 试验结果中，短跑的跟腱受力高达 9 000 N，约 8~12.5 倍的身体重力。对于这种试验方法，有 2 个问题。（1）利用侵入性植入传感器方法会对受试者人体造成伤害，植入性手术存在风险，传感器的敏感性有待验证，漫长的恢复期都会阻碍这种方法的推广。针对本文的研究对象来讲，高水平运动员不可能采用手术方法来进行研究，同时，这种方法也被认为是违背伦理道德的。（2）根据从尸体试验得到的人体跟腱应力时间数据和跟腱与肌肉连接处的横截面积计算，跟腱在运动中的负荷的范围为 1 200~2 000 N，显然力学传感器的试验结果与这一数据相违背，可靠性有待验证。
LICHTWARK采用超声波影像技术研究跟腱在运动中的负荷，跟腱长度定义为跟腱嵌入处到肌肉跟腱连接处。采用数字化影像追踪（Marker）与超声波影像拍摄相结合的技术，得到跟腱在运动中长度拉伸的变化，进而根据跟腱材料特性中长度与应力变化曲线，得出跟腱在运动中的负荷大小。以羊为实验对象，分析走路和爬行动作，得出跟腱伸长与负荷的变化曲线。优点是，对试验对象没有伤害、安全系数较高，而超声波影像技术的应用，准确得到了跟腱在运动中长度的变化，在研究跟腱损伤时该指标价值较高；同时，该方法考虑了跟腱机械特性中的粘弹性特点，增加了试验数据的可靠性。但是，这种方法也存在一定的局限性。（1）基于 Marker 追踪的方法，会因为 Marker 在皮肤上的移动而产生误差，由于跟腱的机械特性，这一误差对试验的最终结果会产生巨大的影响；（2）基于超声波影像技术的研究范围仅应用于人体的简单动作，如步态、原地跳跃、慢跑等不是特别剧烈的动作，但这种方法如果应用于剧烈运动中的跟腱拉伸和负荷，可能会因为试验设备的局限性有所误差，而且，对这类动作的研究目前尚未进行。
3、结 论
（1）与现有跟腱研究方法相比，OpenSim 力学仿真方法在建模和肌肉力计算方面更加先进，仿真结果也比较准确。此外，这种方法对受试者本身没有任何伤害隐患，并且在数据采集和处理方面更加方便、易行，OpenSim 的开源性也为接下来的研究提供了极大的便利。
（2）肌肉-肌腱单位在跨栏运动时比短跑运动产生更大的形变拉伸，且明显差异都出现在着地缓冲阶段，从而导致跟腱在跨栏运动中比短跑运动受到更大的牵拉力，峰值出现的时间有差异，短跑支撑期比跨栏支撑期出现峰值的时间更早，跨栏支撑后期跟腱受力达到峰值。同短跑相比，跨栏跟腱损伤的风险更高，更容易产生伤病。
总而言之，目前的研究提供了新的方法来仿真运动中的跟腱负荷，使用的模型并不局限于计算跟腱受到的牵拉力，同样也可以被用来计算其他肌肉的受力。因而，这个模型可以用来研究其他与运动相关的损伤，进而为伤病的预防提供科学依据。另外，在仿真结果的基础之上，教练员可以以此为科学依据改进训练方法，有针对性地进行相应的肌肉力量和协调性等方面的加强，为运动损伤的预防奠定基础。
参考文献：
[1] 郝卫亚.人体运动的生物力学建模与计算机仿真发展[J].医用生物力学，2011，26（2）：97-103.
[2] HAMILTON B，REMEDIOS D，LOOSEMORE M，et al. Achilles tendonrupture in an elite athlete following multiple injection therapies [J].J SciMed Sport，2008，11（6）：566-573.
[3] AMES P R ，LONGO U G ，DENARO V，et al. Achilles tendon problems：not just an orthopaedic issue [J].Disabil Rehabil，2008，30（20-22）：1646-1650.
[4] DELP S L，ANDERSON F C，ARNOLDA S，et al.OpenSim：Open-sourceSoftware to Create and Analyze Dynamic Simulations of Movement [C].IEEE Trans：Biomed Eng，2007.
[5] DELP S L，LOAN J P.A graphics-based software system to develop andanalyze models of musculoskeletal structures [J].Comput Biol Med，1995，25（1）：21-34.
[6] THELEN D G.Adjustment of muscle mechanics model parameters tosimulate dynamic contractions in older adults [J]. J Biomech Eng，2003，125（1）：70-76.
[7] KUO A D.A least-squares estimation approach to improving the precisionof inverse dynamics computations [J]. Journal of BiomechanicalEngineering，1998，120：148-159.
[8] THELEN D G，ANDERSON F C. Using computed muscle control togenerate forward dynamic simulations of human walking from experimentaldata[J].Journal of Biomechanics，2006，39：1107-1115.
[9] LICHTWARK G A，WILSON A M. A modified Hill muscle model thatpredicts muscle power output and efficiency during sinusoidal lengthchanges[J]. J Exp.Biol，2005，208：2831-2843.
[10] MAGNUSSON S P，KJAER M. Region-specific differences in Achillestendon cross-sectional area in runners and non-runners [J].Eur J ApplPhysiol，2003，90（5-6）：549-563.
[11] HAMNER S R，SETH A，DELP S L.Muscle contributions topropulsion and support during running [J].J biomechanics，2010，43（14）：2709-2725.
[12] GALLOWAY M T，JOKL P，DAYTON O W. Achilles tendon overuseinjuries[J].Clin. Sports Med，1992，11（4）：771-853.
[13] KANNUS P，J魷ZSA L. Histopathological changes preceding spontaneousrupture of a tendon. A controlled study of 891 patients [J]. J Bone JointSurg. Am，1991，73（10）：1507-1531.
[14] VAN DEN BOGERT A J. Computer simulation of locomotion in the horse[D]. Utrecht：University of Utrecht，1989.
[15] LEWIS J L，LEW W D，SCHMIDT J.A note on the application andevaluation of the buckle transducer for the knee ligament forcemeasurement [J].Journal of Biomechanical Engineering，1982，104（2）：125-128.
[16] KOMI P V，SALONEN M，J魧RVINEN M，et al.In vivo registration ofAchilles tendon forces in man I Methodological development [J].International Journal of Sports Medicine，1987（8）1：3-8.
[17] LOUIS-UGBO J，LEESON B，HUTTON W C.Tensile properties of freshhuman calcaneal Achilles tendons[J].Clin Anat，2004，17（1）：30-34.
[18] MAGNUSSON S P，KJAER M.Region-specific differences in Achillestendon cross-sectional area in runners and non-runners [J].Eur J ApplPhysiol，2003，90（5-6）：549-53.
[19] LICHTWARK G A，WILSON A M.Optimal muscle fascicle length andtendon stiffness for maximising gastrocnemius efficiency during humanwalking and running[J].J Theor Biol，2008，252（4）：662-731.
