Document
一种用于OpenSim的非线性缩放方法
伍子豪 , 王正陶 , 马晓宇 , 王冬梅

《医用生物力学》 2023年 39卷 第3期 021
中图分类号:R 318.01
全文 图表 参考文献 作者 出版信息
摘要
关键词
1 材料和方法
1.1 位姿调整
1.2 单调插值函数计算
1.3 生成缩放模型
1.4 方法验证
2 结果与分析
2.1 非线性缩放计算结果
2.2 逆运动学计算结果
2.3 膝关节接触力计算结果
2.4 蒙特卡洛分析结果
3 讨论
3.1 插值曲线意义及缩放误差来源
3.2 方法优点及缺点
4 结论

摘要

目的 提高OpenSim生物力学模型静态与动态配准精度,进而提高运动学和动力学参数计算的可靠性。方法 通过位姿调整、构建插值函数,实现基于运动学实验数据的模型非线性缩放。使用两组公开实验数据(GC3和GC5)进行建模分析,计算模型肢段长度和关节作用力,并与解剖标记点缩放方法 (anatomical landmark scale, ALS)和线性缩放方法的计算结果及实际测量值比对,验证方法的有效性。结果 缩放所得模型的肢段长度与实际测量值的最大偏差14.74 mm,在文献给定范围(4.0±13.8) mm内;缩放和逆运动学计算的标记点误差满足OpenSim给定要求;计算所得关节接触力均方根误差(GC3:0.40 BW,GC5:0.34 BW)较ALS方法(GC3:0.64 BW)和OpenSim线性缩放方法(GC5:0.40 BW)小;蒙特卡洛分析结果同时表明,对于模型标记点初始位置的变化,该方法计算关节接触力误差范围更小,肢段长度波动小于5%。结论 该非线性缩放方法可行,且在目前验证条件下可以提高OpenSim运动学和动力学建模效率及仿真分析结果的精度。

关键词: 非线性缩放 个性化建模 模型配准 运动学 动力学

骨肌生物力学分析是研究人体运动、辅助疾病诊断和治疗方案制定的方法,在进行计算分析前需要先建立与受试者相匹配的模型,缺乏医学影像数据时可以通过对标准模型进行缩放获得 。生物力学分析软件OpenSim提供了一种依据光学标记点数据的线性缩放方法:通过指定的标记点对在模型上的距离和实际距离间的比值确定缩放因子,再对模型上各点的三维坐标值分别乘以相应缩放因子,获得缩放后模型。该方法较为直观,但所得到的模型以及后续运动学计算结果因操作者而异,不同操作者所得结果一致性较差
关于线性缩放方法,有研究基于Andersen等 提出的针对运动学正定或超定问题的计算方法,以及Rasmussen等 提出的基于质量的经验公式计算缩放参数,该方法可以保证结果一致性,且所得模型的标记点与实际标记点误差较小。但Lund等 认为,通过该方法获得的模型一般仅适用于缩放时所使用的动态数据。
与线性缩放相比,非线性缩放方法具有更高的准确性。Nolte等 研究认为,在计算大腿和小腿长度时,非线性缩放方法误差显著小于线性缩放。统计学形状模型(statistical shape model,SSM)是非线性缩放的一种方法,Zhang等 基于此建立了骨肌模型图集项目 (musculoskeletal atlas project, MAP);同时为适应实际采集标记点,Zhang等 提出一种基于稀疏标记点求取主成分系数的优化方法;SSM可以实现高匹配度,但其生成模型前需要先采集既有模型数据建立训练集。有研究团队完成了部分肢段的建立 [8-10] ,但未遍及全身,因而基于SSM的缩放方法前期准备要求较高,更适用于局部模型建立。此外,Lund等 提出通过径向基函数插值对模板模型进行调整,该方法仅需基础模型和标记点数据,但所需标记点数量多于常用标记点集,需要添加标记点或结合医学影像数据。并且,上述方法均未能应用于OpenSim缩放建模。
本文提出一种可用于OpenSim的非线性缩放方法,可以依据常用骨性标记点缩放获得与受试者匹配的计算模型,且免于建立模型训练集,缩放结果可以满足OpenSim后续计算要求 :均方根误差(root mean square error, RMSE)小于10 mm,标记点最大误差小于20 mm,动力学计算结果与实际测量值具有较好一致性。

1 材料和方法

本文所述非线性缩放方法分为3步:位姿调整、单调插值函数计算、生成缩放模型(见 图1 )。建立模型后进行逆运动学计算、静态优化肌肉力计算、关节接触力计算,并与实际测量数据比对验证模型有效性。
图1 非线性缩放方法流程示意图

1.1 位姿调整

调整缩放肢段的初始位姿以调整标记点位置,确保插值计算时模型标记点与实际标记点位置的对应性,缩放后的模型可以与实际相匹配。位姿调整方法参考最近点迭代法(iterative closest point,ICP) ,以肢段旋转和位移量为变量,模型和实际数据中对应标记点的距离平方和为目标函数,考虑关节约束,进行最优化计算。
为完成位姿调整,肢段上标记点数量需满足:
(1)
式中: l 为调整肢段数; n [ { "name": "text", "data": "i" } ] 为第 i 肢段上标记点数; C 表示组合运算;DOF为肢段总自由度。

1.2 单调插值函数计算

在缩放模型时,为避免模型上各点相对位置发生改变而导致原有骨骼模型点云连接关系不合理,加入肢段上模型标记点间和实际标记点间的相对位置相同的假设,并对模型坐标值为自变量、实际坐标值为因变量建立的点对作单调插值,确保模型上各点位置相对关系不变。
计算插值函数分轴向进行,先提取位姿调整后的肢段模型和实际标记点某一轴向坐标,然后分别对其作升序排列,再一一对应建立点对,以模型坐标值为自变量,实际坐标值为因变量,作单调插值计算,求取插值函数。此处采用Fritsch等 提出的单调分段3次插值方法,根据模型和实际坐标点对确定函数系数,以肢段 X 轴向坐标为例,其插值函数形如:
p ( x )= Ax 3 + Bx 2 + Cx + D
(2)
式中: p ( x )为实际标记点坐标; x 为模型标记点坐标。
超出标记点对范围的部分采用通过边界点斜率为1的线性函数:
p ( x )= x +( p 0 - x 0 )
(3)
式中: p 0 x 0 分别为左侧边界实际和模型标记点坐标。
插值曲线上各点切线斜率表示模型在该坐标相关点处的缩放因子,详细讨论见3.1节。

1.3 生成缩放模型

将几何模型点云点、关节中心点、肢段质心分别代入3个轴向的插值函数,计算得到缩放后对应点的坐标,并依据位姿调整计算结果调整回模型初始位置状态,生成几何模型,再依据关节中心点坐标进行组合。
由于计算单调插值函数时假定模型标记点与实际标记点的位置相对关系相同,且在计算前就依据该设定对坐标值进行重新排列和对应,故需要适当调整缩放后模型标记点的位置。调整方法为将实际标记点依据位姿调整计算结果映射至缩放模型位置下,以此作为缩放后的模型标记点。
缩放模型中的肌肉参数(最优肌纤维长度)参考OpenSim的缩放方法 进行计算。

1.4 方法验证

采用第3、5届“Grand Challenge Competition to Predict in vivo Knee Loads”提供的公开数据(后简称为GC3和GC5),其测试对象信息分别为:SC(身高167 cm,体质量78.4 kg,女性),PS(身高180 cm,体质量75 kg,男性),均曾进行左脚膝关节置换术,膝关节假体中置入6轴力传感器用于测量运动过程中的关节作用力,采样频率为120 Hz。步态运动及力数据采集使用8摄像头Vicon动作捕捉系统(Vicon Motion Systems公司, 英国)以及3块Bertec(Bertec 公司,美国)/AMTI(Advanced Mechanical Technology公司,美国)测力板,采样率分别为0.12、1.2 kHz。标记点集为修正的Cleveland Clinic markerset。受试者静态标定时采用双手下垂姿势。
方法验证使用Lai等 全身模型进行计算,将左侧膝关节屈伸相关肌肉的最大等长力缩减35% 。使用GC3和GC5中脚尖向前静态站立数据以及正常步态数据,由静态站立数据通过前述非线性缩放方法建立模型,对比所得模型与数据中提供的三维模型股骨、胫骨长度,后者通过 3-matic 13.0软件(Materialise公司, 比利时)测量获得,再通过标记点误差判断是否满足OpenSim计算要求;接着进行OpenSim逆运动学计算,肌肉力计算采用静态优化方法,最后计算关节接触力,对比关节接触力结果与力传感器数据,通过决定系数 R 2 衡量其相关性、RMSE评价其误差大小。步态数据使用GC3中的ngait5~9,GC5中的ngait(1、3、8、9、11),最终所得 R 2 及RMSE为两组数组中各运动数据所得结果的平均值。
通过蒙特卡洛(Monte Carlo)分析评估缩放前模型标记点位置对最终所得关节接触力的影响,对下肢标记点产生随机方向10 mm偏移 ,重新缩放并计算肢段长度和关节接触力,统计1 000次计算后其数值波动范围。

2 结果与分析

2.1 非线性缩放计算结果

模型初始标记点位置参考文献[16]确定,并根据静态标记点数据和公开数据中的实验照片进行调整,模型标记点与实际标记点一一对应。使用Lai等 全身模型的计算结果表明,误差值满足OpenSim缩放结果要求 ,说明该方法适用于建立OpenSim计算模型(见 表1 )。肢段长度绝对误差值在文献[17]给出的误差范围[(4.0±13.8) mm]内,说明该非线性缩放方法可以获得与文献[17]中线性缩放方法一致的准确性(见 表2 )。
表1 缩放后模型标记点误差
表2 股骨、胫骨长度对比

2.2 逆运动学计算结果

OpenSim提出逆运动学计算结果可靠条件 :RMSE<20 mm,标记点最大误差小于40 mm。结果表明,非线性缩放后得到模型的逆运动学计算结果可以满足要求(见 图2 )。
图2 不同数据逆运动学计算标记点误差

2.3 膝关节接触力计算结果

对比模型计算值与测量值,计算 R 2 与RMSE,并与文献结果比较。为保证计算所用数据中受试者所受支持力可完全由测力板反映,GC3截取的步态周期为右脚摆动相开始至下一次摆动相开始,GC5截取的步态周期为左脚支撑相开始至下一次支撑相开始。选用周期不同,故关节力曲线变化存在差异(见 图3 )。本文认为,关节接触力计算值与测量值变化趋势一致,具有一定可靠性。
图3 不同数据关节接触力计算值与测量值对比
表3 对比结果可见,相较于使用径向基函数插值的解剖标记点缩放方法(anatomical landmark scale,ALS)方法 ,本研究提出的非线性缩放方法获得的模型的关节接触力计算结果具有更小的RMSE,而决定系数相当;相较于OpenSim线性缩放方法,则具有更小的RMSE与更大的决定系数。以上结果说明,该非线性缩放建模方法可以获得比上述方法更优的动力学计算结果。
表3 关节接触力计算结果与对比

2.4 蒙特卡洛分析结果

使用GC3的ngait5数据进行蒙特卡洛分析,其所得关节接触力波动范围,以及关节接触力RMSE和 R 2 、股骨长度和胫骨长度分布情况如 图4 所示。
图4 蒙特卡洛分析结果
经过1 000次随机更改模型标记点坐标并进行非线性缩放、逆运动学、静态优化、关节力计算后,所得关节接触力RMSE变化范围为0.19 BW,该变化范围以及RMSE的最大值和最小值较ALS方法所得结果 小,说明本研究方法受缩放前模型标记点位置变化的影响相对较小,且关节接触力曲线变化趋势与实际值保持一致,肢段长度变化范围相对平均值的比例小于5%(股骨:4.92%,胫骨:3.67%),因而该非线性缩放方法具有一定的稳定性(见 表4 )。
表4 蒙特卡洛分析结果统计

3 讨论

3.1 插值曲线意义及缩放误差来源

缩放计算过程中由模型标记点和实际标记点坐标获得的不同轴向插值曲线,曲线上各点的切线斜率表示模板模型在该点处的缩放因子[见 图5 (a)],因而不同坐标轴、同一坐标轴方向上的不同位置具有不一样的缩放因子,即实现非线性缩放。
由插值曲线的这一意义,可以将线性缩放整合到该方法中。线性缩放的特点是缩放因子处处相等,在调整模型初始位姿获得调整标记点坐标后,对各轴向模型与实际标记点坐标关系散点图作线性拟合,所得直线斜率即是该轴向的线性缩放因子[见 图5 (b)];这样所得缩放因子会较OpenSim直接指定标记点对计算的缩放因子带来更小的缩放模型标记点误差,因为指定标记点对只能确保该轴向上部分标记点误差最小。
图5 (b)中也可以看出,通过线性拟合获得的缩放因子可以确保误差最小,但标记点对不完全位于拟合直线上,因而与建立插值函数的非线性缩放方法相比,线性缩放存在相对较大的标记点误差。
图5 GC3右侧股骨
本文所提出的非线性缩放方法仍然存在标记点误差,推测原因如下:① 计算过程中为避免出现过大斜率而对部分标记点进行合并处理;② 为不改变点云模型连接关系而对模型与实际中相对位置不一致的标记点进行合并处理。

3.2 方法优点及缺点

相比于OpenSim现有的线性缩放方法,在可以获得同样满足计算要求模型的情况下,该非线性缩放方法由于不用调整用于计算缩放因子的标记点对和模型标记点位置(除与实际标记点位置明显不匹配的情况外),因而缩放操作用时较少;同时不具有操作者依赖性,可以提高计算流程标准性。相比于Zhang等 提出的SSM方法以及Lund等 提出的径向基函数插值方法,该方法可以免去模型训练集准备和优化迭代计算以及可以用于标记点较少情况下的缩放计算,并在目前验证下获得与实际测量值较为一致的动力学计算结果,且优于同类型ALS方法 与OpenSim线性缩放方法 ,同时在针对标记点初始位置变动的蒙特卡洛分析中表现出一定的稳定性。该方法还适用于不同标记点集(如Helen Hayes model、 Cleveland model、 VCM model等),只要其标记点在位姿调整时可以满足数量要求,因而具有一定的灵活性。
该方法也存在一定的不足:① 由于所需输入信息较少,该非线性缩放方法获得模型的细节信息相对不足,且缩放后所得的模型骨块可能出现不合理的变形,可能表现为局部特征被破坏(如接近于球形的股骨头在缩放后失去球形特征),这种情况可以参考Kraevoy等 提出的基于易变性(vulnerability)计算结果对易变局部进行等比例缩放以保留其特征的方法,改进缩放方法;② 单调插值方法会带来一定的局限性,此处所引入的同一肢段上模型标记点间与实际标记点间具有相同的位置相对关系假设可能与实际不相符,这样缩放结果就不能真正反映分析对象的实际情况,此时可能需要借助SSM,对其生成的模板模型进行非等比例缩放,可以获得与实际更为相近的缩放结果;③ 未充分考虑多关节部位缩放,如肩部以及脊柱中的椎关节,如果执行标定姿势时这些部位处关节活动较小,不作考虑也不会引入明显误差,否则应当将其关节运动加入到缩放计算中。
此外,本研究中仅采用两组公开数据进行计算验证及方法对比,说明方法可行且在目前条件下具有有效性,为提高结论的普遍性需对更多数据进行建模分析,作进一步验证;关于不同的肌肉模型和肌肉力计算方法对缩放模型关节接触力计算结果的影响仍需分析讨论 ;该非线性缩放方法也可以参考Kinematically Scale方法 ,加入通过关节运动获得的关节中心和关节转轴信息,作进一步优化。

4 结论

本文提出的非线性缩放方法可基于光学运动捕捉标记点数据建立有效的OpenSim计算分析模型,基于两组公开数据进行计算验证,其中模型的运动学计算结果可以满足OpenSim计算要求,动力学计算结果优于ALS方法和OpenSim线性缩放方法,因此该方法可行,且在目前验证条件下提高了OpenSim建立运动学和动力学分析模型的效率及仿真分析结果的精度。