长距离跑后“足外翻”姿态增加膝关节内侧接触力：基于OpenSim肌骨建模及机器学习预测的研究

2019-11-16梅齐昌相亮亮李建设JustinFernandez顾耀东

体育科学 2019年9期

梅齐昌，相亮亮，孙冬，李建设，Justin Fernandez，顾耀东*

（1.宁波大学体育学院，浙江宁波315211；2.宁波大学大健康研究院，浙江宁波315211；3.奥克兰大学生物医学工程研究所，奥克兰1010）

0 前言

近年以来，跑步运动因其普适性及对疾病的预防功效广受欢迎且参与人数快速增长（Hulme et al.，2017；van Gent et al.，2007）。随着跑步相关损伤的概率不断增大，据报道，每1 000 h跑步中约有33例损伤（Videbæk et al.，2015），其中，膝损伤比例高达 79.3%（van Gent et al.，2007）。与外界接触的足部在跑步时常会出现形态和姿态及应力的变化（顾耀东等，2019；李蜀东等，2018；梅齐昌等，2015；Mei et al.，2015）；研究发现长跑后足纵弓出现塌陷（Mei et al.，2018），从足解剖结构和人体动力链角度分析（张燊等，2018；赵晓光，2018），足弓下塌常伴随足外翻，是导致跑步相关损伤的一个关键内部因素（Barnes et al.，2008；Mei et al.，2018；Nigg 2011；Nigg et al.，2015）。除通过三维扫描评估足形态外（Mei et al.，2015），还可结合临床及运动中已广泛应用的足部姿态指数（Foot Posture Index，FPI）（Redmond et al.，2006），通过观察及触摸快速地评估足姿态，避免繁琐昂贵的影像学测试及分析。

超过90%的业余马拉松跑者是后跟落地的跑姿（Larson et al.，2011），实时监测的研究发现，业余跑者力竭性5 km跑后足部的局部动态稳定性下降显著且低于专业跑者（Hoenig et al.，2019），表明业余跑者对远端环节的控制能力减弱，可能与控制足部的肌群功能较弱有关（张燊等，2018；赵晓光，2018），表现出足弓下塌及足外翻；另一个表现是主要关节输出功率从远端的踝关节向近端的髋关节转移（Sanno et al.，2018）。关节冲击及足外翻是可预测跑步相关损伤的因子，相关跑步训练及实时监测的研究通过反馈调节冲击相关的参数，如胫骨最大冲击及加载率（Bowser et al.，2018；Yong et al.，2018），来探究跑步相关损伤风险的变化。虽然跑步鞋对足外翻的功效性尚未达成一致定论（郝琦等，2012），但业余跑者穿着动作控制鞋能够缓解足外翻而降低潜在损伤风险的观点已被证实（Malisoux et al.，2016）。个体化肌骨建模技术在临床生物力学中有广泛应用（郝卫亚，2011），除计算关节角度和力矩外，可通过优化算法预估肌肉激活及肌肉力用于模拟计算关节接触力和韧带组织的受力（黄尚军等，2018；李琳杰等，2013；戎科等2015）。

监测跑步时关节负荷并动态调节跑步量可预防跑步相关损伤（Bertelsen et al.，2017）。本研究运用个体化肌骨建模技术揭示跑者具体的关节负荷，如膝关节接触力的分布及幅值变化，并与足姿态参数通过偏最小二乘回归（partial least squares regression，PLSR）算法作线性关联。前期长距离跑的研究发现，相关变化均发生在5 km左右（Hanley et al.，2018；Hoenig et al.，2019；Sanno et al.，2018）。因此，本研究选用5 km跑为业余跑者的长距离跑测试，旨在揭示长跑后足姿态及下肢关节角度、力矩和膝关节接触力的变化，并构建多元回归PLSR机器学习算法模型，训练、验证并测试足姿态参数与膝关节总接触力和膝关节内、外侧接触力的线性关系，为通过足姿态参数快速预测膝关节接触力负荷的分布及幅值变化提供科学参考依据。

1 研究对象与方法

1.1 研究对象

本研究共招募了20名业余男性跑者（年龄：25.8±1.6岁；体重：67.8±5.3 kg；身高：1.73±0.05 m），跑姿均为后跟落地，均有跑台、跑道或公路跑步史，每周跑量约30 km。均无任何下肢畸形或损伤及心血管系统疾病。在招募阶段，受试者的足部均按照Foot Posture Index（足部姿态指数）进行筛查，足部畸形（如拇外翻、足过度内翻或外翻、扁平足或高弓足等）均被排除。测试前，受试者均签署了相关知情说明书，并了解测试的目的、方法及步骤。该测试获得了所在单位的科学研究伦理委员会批准。

1.2 测试方案

测试内容包括裸足静态站立时的足部姿态评估，该测试参照国际通用的Foot Posture Index操作指南（Redmond et al.，2006），包括距骨触诊、外侧踝尖、跟骨内外翻、距舟关节突起、内侧足纵弓及前足展收的观察6项指标评分（-2-112），依据总评分划分足姿态为过度内翻（-12）、内翻（-5）、正常（0）、外翻（5）及过度外翻（12）；8个摄像头的Vicon三维运动捕捉系统（Vicon Metrics Ltd.，Oxford，UK）及1块位于20 m跑道中间的AMTI三维测力台（AMTI，Watertown，Massachusetts，USA）被用于记录37个反光Mark点的运动轨迹（Hamner et al.，2013；Rajagopal et al.，2016）和跑者1个支撑期的地面反作用力，测试频率分别是200 Hz和1 000 Hz；16通道的Delsys表面肌电测试系统（delsys，boston，Massachusetts，US）被用于收集股直肌（rectus femoris）、股外侧肌（vastus lateralis）、股内侧肌（vastus medialis）、股二头肌（biceps femoris）、半腱肌（semitendinosus）、胫骨前肌（tibialis anterior）、内侧和外侧腓肠肌（medial gastrocnemius and lateral gastrocnemius）8块肌肉的活动情况，频率设定在1 000 Hz，以供模型验证。

5 km跑台跑前，受试者均先进行热身活动及实验室环境熟悉及步态适应，后开展上述测试，分别以肩宽站立时评估足部姿态，并记录以自选舒适速度跑步时Mark点运动轨迹和地面反作用力及肌肉的表面肌电活动情况，每位受试者收集5次有效数据；结束5 km跑后的5 min内，按照上述步骤重复5 km跑后的跑步测试。

1.3 肌骨模型构建、数据处理及模型验证

本研究升级了标准OpenSim 2392肌骨模型（23个活动度及92个肌肉-肌腱单位），添加无摩擦力的髌骨以引导股四头肌的肌力矩，股四头肌末端包绕该髌骨并止于胫骨粗隆（图1）（DeMers et al.，2014）。鉴于该研究主要目的在于揭示膝关节接触力的分布及幅值变化，模型编辑过程中仅保留了膝关节的屈伸自由度，避免展收及回旋活动度对接触力的可能影响。该模型采用Besier等（2003）的测试方案确定了髋关节（hip joint centre）和膝关节（knee joint centre）中心点，结合Saxby等（2016）已发表的方案，可计算出股骨的长度（femur length）（公式1）、膝关节在股骨上的中心点（knee joint centre femur）（公式2）；结合膝关节内外上髁（medial condyle、lateral condyle）的距离计算出膝关节展收的力臂（knee add/abd moment arm）（公式3），并运用膝关节展收力臂计算股胫关节的股骨外侧髁和内侧髁在胫骨平台上的外侧（knee lateral compartment）（公式 4）及内侧（knee medial compartment）（公式 5）接触点，从而用于计算外侧（Flc）和内侧（Fmc）点的接触力（图1）。

图1 肌骨模型构建及膝关节内外侧接触力示意图Figure 1.Illustration of Subject Specific Modelling and Medial and Lateral Knee Contact Force

本研究主要分析右侧下肢在支撑期内的关节运动学、动力学和接触力的数据，将实验收集的Mark点运动轨迹及地面反作用力数据采用自编的MATLAB程序进行处理及格式转换（C3D文件转换成trc及mot文件），并进行6 Hz（Mark点轨迹）和30 Hz（地面反作用力）的Butterworth零滞后四阶低通滤波。采用上述公式1～5的算法获得模型膝关节的内外侧接触点，依照Delp等（2007）发表的分析流程：1）对模型进行缩放标定（scaling）以获得符合各受试者人体测量学的模型，肌肉起止点及肌力臂均与受试者肢体长度相一致；2）运用逆向运动学（inverse kinematics）算法（能够缩小模型虚拟点与实验记录Mark点间的误差）计算关节角度；3）运用逆向动力学（inverse dynamics）算法计算关节的力矩；4）运用静态优化（static optimization）算法计算肌肉的激活度（activations）（用于对比实验记录的表面肌电活动验证模型）及肌肉力（forces）；5）运用分析（analysis）工具关节反作用（joint reaction）算法计算膝关节的总接触力（total knee contact force）、膝关节内侧（Fmc）及外侧（Flc）接触力。

除记录支撑期内的肌肉活动外，测试前均记录了上述肌肉最大肌力收缩（maximum voluntary contraction，MVC）时的肌肉活动水平；在Delsys EMGworks分析软件中对肌电数据进行均方根（root mean square，RMS）算法处理，用支撑期内的肌肉活动水平除以最大肌力收缩时肌肉活动水平获得标准化后的肌肉激活水平（0表示无激活，1表示完全激活），将实验测得的肌肉激活水平与静态优化算法获得的肌肉激活度进行对比从而完成模型验证。

1.4 统计分析及偏最小二乘回归算法模型

统计学分析5 km跑前和跑后的跑步测试中右侧下肢支撑期内的相关数据，包括足部姿态指数（FPI）、跑速（m/s）、支撑期时长（s）、下肢关节角度（°）、下肢关节力矩（Nm/kg）及膝关节的接触力（body weight，BW）。统计学测试前对数据进行正态分布检验，后采用配对样本t检验对跑前和跑后的相关参数进行统计学检验。支撑期被划分为支撑初期（0%～50%）、支撑中期（～50%～）及蹬离期（50%～100%）（Novacheck et al.，1998）；鉴于关节角度、关节力矩和接触力随时间变化的一维特征，Statistical Parametric Mapping（SPM）1d（一维统计参数映射）程序是现通用的连续性数据分析方法，可检验包含极值以外其他数据的统计学差异（Pataky 2010）；本文采用SPM1d v0.4对上述数据进行配对样本t检验，所有分析均在MATLAB软件（MATLAB R2018a，The MathWorks，MA，USA）中完成，显著性水平均设定在0.05。

偏最小二乘回归（partial least squares regression，PLSR）是一种机器学习算法，通过多因变量（反应集Y）对多自变量（预测集X）的线性回归建模，可运用至新的自变量数集预测新的因变量数集，在骨骼三维形态、有限元模拟应力分布及肌肉力学等研究中已有大量应用（Fernandez et al.，2016，2019；Zhang et al.，2016）。本研究自编了MATLAB偏最小二乘回归（PLSR）算法程序，包含预测（predictors）矩阵集（X）和反应（responses）矩阵集（Y），结合“留一交叉检验法”（leave-one-out cross validation）训练并测试预测集（足姿态指数）与反应集（包括膝关节总接触力、膝关节内侧接触力及膝关节外侧接触力等）的关联。

预测矩阵集（X）（公式6）和反应矩阵集（Y）（公式7），如下公式：

其中，N是数据总量，M是预测因子总数，P是反应因子总数，L是主成分数量。T和U分别是预测集和反应集的投影矩阵，P和Q分别是预测集和反应集的转置正交载荷矩阵，E和F分别是预测集和反应集的缀余量。投影矩阵的矢量关系用如下线性方程（公式8）表示，其中，H是赘余矩阵的矢量。

该偏最小二乘回归（PLSR）算法模型采用数据集中80%样本量用作训练集，20%样本量用作测试集。首先，通过“留一法”交叉检验该模型的训练集；其次，运用测试集的新数据集检验该模型。预测因子（足姿态参数）的均值是4.65（最大可达过度外翻值：12），为验证足姿态参数对膝关节总接触力及膝关节内、外侧接触力的影响，通过灵敏度分析（sensitivity analysis）分别按照 5%、10%、15% 至250%的递增摄动（perturbation）预测因子，评估反应集（膝关节总接触力及膝关节内、外侧接触力）的变化。

2 结果

2.1 模型验证

将实验记录的表面肌电（EMG）信号通过RMS处理转化为激活度（0表示完全无激活，1表示完全激活），与OpenSim优化算法处理获得的激活度（activation）相对比，包括股直肌（rectus femoris）、股外侧肌（vastus lateralis）、股内侧肌（vastus medialis）、股二头肌长头（biceps femoris lh）、半腱肌（semitendinosus）、胫骨前肌（tibialis anterior）、内侧腓肠肌（medial gastrocnemius）和外侧腓肠肌（lateral gastrocnemius）。如图2所示，在支撑期内5 km跑前及跑后计算获得的肌肉激活度与实验记录的表面EMG信号活动基本相似，表明本研究中的OpenSim模型数据较为可靠（Hamner et al.，2013；Rajagopal et al.，2016）；通过Open-Sim逆向运动学及逆向动力学算法计算出的关节角度及关节力矩也与其他相关研究做对比，趋势及幅值均相近（Bonacci et al.，2013；Hamner et al.，2013）。

图2 EMG信号与5 km跑前（Pre 5k）及跑后（Post 5k）OpenSim优化算法获得肌肉激活水平对比Figure 2.Comparison of Muscle EMG Signals and Activations from OpenSim Optimization at Pre and Post 5 km Running

2.2 足姿态及步态参数

5公里跑前与跑后的足姿态指数表现出显著差异，呈现出外翻姿态；测试时的跑速细微增快且表现出显著性；支撑期时长未表现显著变化（表1）。

表1 5 km跑前（Pre 5k）及跑后（Post 5k）的足姿态参数（FPI）、跑速及支撑时间Table 1 FPI Scores，Speed and Contact Time at Pre and Post 5km Running（M±SD，95%CI）

2.3 下肢关节角度及力矩

5 km跑后髋关节在蹬离期82%～100%伸的角度显著增大（P=0.031），在支撑初期0%～8%外旋角度显著增大（P=0.036）；踝关节在蹬离期76%～92%跖屈角度显著增大（P=0.006）；其他运动面的关节角度未呈现出显著性差异（图3）。

5 km跑后，髋关节在支撑期20%～76%伸的力矩显著增大（P＜0.001），在支撑初期 12%～18%（P=0.003）、28%～30%（P=0.03）、34%（P=0.05）、38%～48%（P＜0.001）时段外展力矩显著增大，在支撑期10%～20%时段髋关节外旋力矩显著增大（P＜0.001）；膝关节的伸力矩在支撑初期8%（P=0.05）、22%（P=0.047）及支撑中期到蹬离期的 44%～98%（P＜0.001）时段均表现出显著性；踝关节跖屈力矩在整个支撑期6%～98%时间均表现出差异性（P＜0.001）；其他运动面的关节力矩未呈现出显著性差异（图4）。

2.4 膝关节接触力

总接触力在蹬离期70%～76%表现出显著性差异（P=0.008）；内侧接触力在支撑初期4%（P=0.05）、12%～20%（P=0.001）、26%～30%（P=0.015）时段显著增大；外侧接触力在支撑期 8%～16%（P=0.002）、22%（P=0.046）及 40%～44%（P=0.016）及蹬离期 68%～76%（P=0.003）表现出显著性差异（图5）。

图3 5 km跑前（Pre 5k）及跑后（Post 5k）下肢关节角度对比及统计学分析Figure 3.Statistical Comparions of Joint Angles in the Lower Extremity at Pre and Post 5 km Running

图4 5 km跑前（Pre 5k）及跑后（Post 5k）下肢关节力矩对比及统计学分析Figure 4.Statistical Comparions of Joint Moments in the Lower Extremity at Pre and Post 5 km Running

2.5 偏最小二乘回归（PLSR）模型预测及灵敏性分析

通过训练集数据对PLSR模型的“留一法”交叉检验，该模型可预测87.15%的膝关节接触力数据，膝关节总接触力、内侧接触力及外侧接触力的均方根误差（root mean square error，RMSE）分别为1.07、0.73及0.51；结合全新的数据集对该PLSR模型进行测试，发现该模型可预测新数据集中94.54%的膝关节接触力数据，预测准确率提高；其中，总接触力、内侧接触力及外侧接触力的均方根误差分别是0.67、0.38及0.41（图6）。基于摄动足姿态参数的PLSR模型的灵敏性分析发现（图7），随着足姿态的外翻程度加大，膝关节内侧接触力呈快速增大的变化，膝关节总接触力增大趋势略缓；相反，膝外侧接触力表现出缓慢减小的趋势。

图5 5 km跑前（Pre 5k）及跑后（Post 5k）膝关节接触力对比及统计学分析Figure 5.Statistical Comparions of Knee Contact Force at Pre and Post 5 km Running

图6 PLSR模型训练数集与测试数集预测准确性Figure 6.The Prediction Accuracy of Training and Testing Dataset in the PLSR Model

图7 PLSR模型灵敏分析Figure 7.Sensitivity Analysis of the PLSR Model

3 讨论与分析

本研究采用个体化OpenSim肌骨建模揭示了5 km跑后业余跑者下肢关节的关节角度、力矩及膝关节接触力的变化，并结合足姿态指数量化了足姿态的外翻变化；构建且验证了偏最小二乘回归（PLSR）算法运用足姿态参数预测膝关节总接触力及内、外侧最大接触力的机器学习模型。研究发现，髋关节蹬离期伸角度显著增大，支撑初期外旋角度显著增大，踝关节在蹬离期跖屈角度显著增大；髋关节的伸力矩在支撑期及外展力矩在支撑初期均显著增大；膝关节屈力矩在支撑初期及蹬离期均显著降低；踝关节跖屈在支撑期显著增大；膝关节总接触力在蹬离期显著减小，内侧接触力在支撑初期显著增大而外侧接触力在支撑初期及蹬离期均显著减小；机器学习算法PLSR模型发现，随着足外翻程度增大，膝关节总接触力及内侧接触力快速增大，外侧接触力缓慢减小。

5 km跑后，业余跑者的支撑时间未表现出任何显著性，与近期的一项探究下肢环节做功的研究结果相一致（Sanno et al.，2018）。本文中跑者在5 km跑前及跑后的跑速尽管表现出统计学差异，仍在3.1 m/s左右（Pre 5k：3.068 m/s；Post 5k：3.137 m/s），与其他研究中的业余跑者长跑时5 km的跑速3.2 m/s（Chan-Roper et al.，2012）及3.3 m/s（Hoenig et al.，2019）均相近，但低于竞速跑者的配速（4.7 m/s）（Hanley et al.，2018）。足姿态在跑后向外翻姿态的改变与前期相关研究也相一致，表明业余跑者的足部在长距离跑步重复冲击后会出现足弓塌陷及足外翻的变化（Fukano et al.，2018；Mei et al.，2018）。

本文结合当前国际生物力学研究最常用的一维统计参数映射（one-dimensional statistical parametric mapping）分析方法（Pataky，2010），对下肢关节角度、力矩及膝关节接触力进行统计分析；不同于传统的统计学方法主要分析特征值（极值）的差异性，该方法基于随机场理论（Radom Field Theory，RFT）对多维度平滑连续性变化的数据进行拓扑分析（Pataky，2012），可揭示出包含特征值及以外的其他连续性数据的差异性。运动科学（特别是运动生物力学）中的数据常随时间变化（连续性），传统的数据检验仅分析极值或某特征值的显著性，运用该方法可分析整个时间周期内的数据差异性，能够发现特征值以外的差异信息。本研究发现，髋关节在蹬离期伸的角度显著增大，在支撑早期外旋的角度显著增大，与业余跑者的10 km跑台跑步测试时5 km所表现出的变化相近（Sanno et al.，2018）；而髋关节的伸力矩、外展力矩及外旋力矩均在跑后显著增大，与业余跑者中已发现的髋关节屈、外展及外旋肌力变化相关（Kollock et al.，2016），导致关节力矩的变化，且与业余跑者长跑后下肢环节做功向近端环节转移的结论相一致（Sanno et al.，2018）。对长跑运动员特定的肌力训练已证实可改变跑者的跑步表现且提高经济性（薛锋等，2014）。

本研究中业余跑者的膝关节屈力矩在支撑早期至蹬离期均显著降低，可归因于业余跑者显著减弱的膝关节屈肌群肌力（Kollock et al.，2016），这也是跑步相关损伤产生的一个预警信号。肌力不足易导致关节代偿性的承受重复冲击负荷，本研究发现，膝关节内侧接触力显著性增大而外侧显著减小，表明负荷向膝内侧转移，与研究发现业余跑者30 min跑后膝关节软骨显著变形的结果相一致（Boocock et al.，2009）。尽管运动学关节角度未表现出差异，与相关研究的趋势及幅值均一致（Bonacci et al.，2013；Hamner et al.，2013），表明常规关节角度和力矩未出现变化，但关节负荷已重新分布（Bertelsen et al.，2017），也是导致膝关节损伤甚至关节炎的早期信号。踝关节在支撑期内的关节角度及关节力矩也均与近期相关研究相一致；5 km跑后，蹬离期表现出显著增大的跖屈角度及整个支撑期内显著增大跖屈力矩可归因于本研究中发现的足姿态外翻及外翻的跟骨，与前期所报道的距下关节（跟骨）外翻能够增大跖屈肌群的活动及肌力从而加大跖屈肌力距的研究发现相一致（Wang et al.，2011）；对比模型验证中肌骨模型计算与实验测得的肌电数据也发现，在跑后内、外侧腓肠肌的激活度比跑前更大。本研究中未出现差异性的距下关节的角度和力矩可能与测试采用的单一置于跟骨后端的Marker点有关，也是本文的一点不足之处，后期研究建议可采用多Marker点集，如Oxford牛津足多环节模型（Sun et al.，2018），以揭示足部各环节间的动态变化。本文的一个创新点是采用了一维统计参数映射（SPM1d）分析法检验了下肢关节的角度、力矩及接触力在支撑期内的连续性数据差异性，并证实了该方法的可行性，为该方法在运动科学研究中的应用铺垫了基础，揭示随时间变化的运动科学数据差异信息可为运动科学工作者提供参考信息。

跑步的支撑期内，足弓是一个动态的能量传递链，能够缓解冲击同时产生势能辅助蹬离（张燊等，2018）。然而，跑步是周期性循环重复的动作，鉴于长跑的重复冲击，已有对足弓塌陷及足外翻姿态的相关报道与本研究中发现跑后足姿态外翻变化的结果相一致。尽管足外翻是否直接导致跑步相关损伤尚无一致定论，然而，高弓足跑者的踝关节损伤概率较高；相反，低足弓及外翻足跑者的膝关节损伤更为常见，尤其是胫骨内侧应力综合征等损伤（Malisoux et al.，2016）。需考虑的一点是，本研究中的受试者均着自选跑鞋，未要求统一控制鞋具，尽管前期研究已发现，不同鞋具或裸足会对跑步的动作控制有一定的影响，本研究不控制穿着鞋具因素的目的旨在揭示大部分业余跑者当前均无具体功能性鞋具需求（仅考虑舒适性）的现状。

本研究的另一个创新点在于构建了机器学习PLSR模型，训练且验证了量化足姿态参数与肌骨建模计算的膝关节负荷之间的线性关联，在长距离跑后外翻足会增加膝关节内侧接触力负荷。前者的足姿态参数可通过实地简易测试获得（避免复杂昂贵的影像学测试），后者的膝关节接触力需繁琐的步态测试及肌骨建模计算。基于该算法模型及灵敏度分析，发现业余跑者随着足姿态参数外翻程度加剧，膝关节总接触力及内侧接触力快速增大，表明外翻足使膝关节负荷向内侧聚集的趋势，与已发现的膝关节内侧应力综合征及内侧膝关节炎相一致（Boocock et al.，2009）。结合该模型，业余跑者可在单次长跑活动后量化记录自身足姿态变化，从而快速预测膝关节负荷幅值及分布的变化，了解并监测自身跑步负荷积累状况，合理安排跑量及训练计划，从而有效避免跑步相关损伤的发生。