唐 刚,程荣山,胡 雄,张博一,吴 钢

(上海海事大学 物流工程学院,上海 201306)

基于Matlab的人体膝关节动力学仿真平台设计

唐 刚,程荣山,胡 雄,张博一,吴 钢

(上海海事大学 物流工程学院,上海 201306)

在对小腿运动学和膝关节动力学分析的基础上,开发了一款基于Matlab的各种步态下人体膝关节动力学仿真平台. 基于反向动力学原理,编写了Matlab计算程序,设计了仿真平台的图形用户界面(GUI),通过读取运动捕捉系统采集的下肢运动信息及三维测力台的支反力,实现了各种步态下膝关节关节力、关节力矩的快速计算. 以正常步态为例,通过对一名健康人体右侧膝关节分析,并与Visual 3D软件结果进行对比,验证了该平台的可靠性,从而为人工膝关节置换和膝关节康复运动评定提供了一种新的分析工具.

膝关节动力学;Matlab;批量读取;快速计算;步态

随着计算机技术的迭代更新,应用计算机对人体下肢膝关节进行建模和动态仿真变得越来越重要[1-7],如下肢膝关节置换前的力学分析、患者下肢膝关节的康复运动评定、膝关节风湿病所引起的非正常行走的分析等. 膝关节是下肢运动的最主要关节之一,对各种步态运动的研究起着非常显著的作用. 同时,膝关节又承载着人体的绝大部分重量,由于它的受力活动频繁,膝关节更易受到伤害. 因此,研究膝关节的动力学参数对人体生物力学、医疗等领域也具有重要意义[8-10]. 目前,国外成熟的运动分析软件主要有Visual 3D步态分析软件等. 虽然这些软件精度高,但是需要人工手动建模,建模速度慢,且需要专业的建模知识,耗费的时间较长,大大减慢了数据处理的速度,因此,设计一个可以快速自动建模,且精度在许可范围内相对较高的仿真平台是本文要解决的问题.

本文基于反向动力学原理,利用Matlab编写了人体膝关节的动力学程序,设计了仿真平台的图形用户界面(GUI),通过读取运动捕捉系统采集的下肢运动信息及三维测力台的支反力数据,实现了膝关节动力学参数的输出.

1 试验数据采集

1.1 试验对象与试验设备

试验对象为一名健康男性,身高1.80 m,体重67.0 kg,无关节或肌肉损伤病史.

试验采用OPTOTRAK CERTUS运动捕捉系统、两台AMTI测力平台同步测量. 运动捕捉系统主要由位置传感器、控制中心和Marker点组成,采样频率为128 Hz. AMTI测力平台主要由两块测力板组成,采样频率为512 Hz. 试验前首先将Marker点粘贴在加热弯曲的有机玻璃板上,然后将带有Marker点的有机玻璃板贴附在试验对象大腿、小腿和足部的侧前部,并用弹性绷带固定. 试验过程中按要求完成5次符合试验要求的步态(如图1所示).

图1 步态测量Fig.1 Gait measurement

1.2 试验测量

由于运动捕捉系统捕捉的Marker点要用全局坐标系表示,需要建立统一的全局坐标系. 本文则以测力板的某直角交点为圆心,两直角边分别为x轴方向和y轴方向,竖直向上为z轴方向. 全局坐标系确定后,对试验对象进行正常步态测量. 试验要求: 试验前对试验对象进行适应性训练直至按节拍行走,试验中则要求试验对象的行走节奏与节拍器相同及步态的步长保持稳定. 为保证测力板上的支反力准确,试验过程中仅单脚接触测力板. 试验测量中,运动捕捉系统每次采集640帧数据,测力板每次采集2 560帧数据,数据采集结束后,运动捕捉系统采集的下肢Marker点的空间坐标原始数据和测力板采集的支反力保存为相应名称的xls格式文件.

2 基于Matlab GUI人体膝关节动力学仿真平台的实现

2.1 平台运算的主要参数

平台运算的主要参数包括小腿运动学参数(小腿转动角、小腿角速度、小腿角加速度)和膝关节动力学参数(膝关节力、膝关节力矩). 各参数的数学符号和在本试验中的物理意义如表1所示.

表1 小腿运动学参数和膝关节动力学参数Table 1 Kinematics parameters for shank and dynamics parameters for knee joint

2.2 平台的算法设计

本平台的算法设计基于反向动力学原理,引入空间向量坐标编写Matlab程序求解小腿的运动学参数,利用小腿运动学参数进而求得膝关节动力学参数[11]. 仿真平台算法设计主要围绕膝关节动力学算法展开.

为了获得膝关节动力学参数首先对小腿进行运动学研究,主要包括对小腿角速度和角加速度的分析. 小腿段对应踝关节和膝关节,小腿段转动角即为局部坐标系下小腿围绕x轴、y轴、z轴的旋转角度. 每一帧对应一个小腿转动角度值,相邻两帧的角度变化量dθ与两帧时间间隔dt的商即为该时间内的小腿平均角速度ω,角加速度αi的计算与角速度基本相同. 小腿角速度ωi、角加速度αi涉及的算法以及对应的Matlab的程序(图2)如下所述.

图2 小腿角速度和角加速度的Matlab程序图Fig.2 Matlab program diagram of shank angular velocity and angular acceleration

踝关节中心的向量表示为

R1=(RLA+RMA)/2

(1)

其中:RLA为踝外侧Marker点坐标;RMA为踝内侧Marker点坐标.

膝关节中心的向量表示为

R2=(RLK+RMK)/2

(2)

其中:RLK为膝外侧Marker点坐标;RMK为膝内侧Marker点坐标.

在小腿的局部坐标系中,za轴由踝关节中心指向膝关节中心:

za=R1-R2

(3)

ya轴垂直于由踝关节中心、RLK、RMK三点确定的平面:

ya=za×(RMK-RLK)

(4)

根据右手螺旋定则进一步确定xa轴:

xa=ya×za

(5)

由于每一帧的坐标已知,小腿绕x轴旋转α角,旋转后x轴方向的坐标值不变,y轴和z轴方向的坐标值变化相当于在yOz平面内作正α角旋转:

(6)

小腿绕y轴旋转β角,旋转后y轴方向的坐标值不变,z轴和x轴方向的坐标值变化相当于在zOx平面内作正β角旋转:

(7)

小腿绕z轴旋转γ角,旋转后z轴方向的坐标值不变,x轴和y轴方向的坐标值变化相当于在xOy平面内作正γ角旋转:

(8)

综上,小腿的旋转角为

(9)

小腿角速度ωi、角加速度αi的公式为

(10)

利用上文获得的小腿运动学分析参数对膝关节进行动力学研究,主要包括对膝关节关节力和关节力矩的分析. 通常在人体步态运动研究中会采用多刚体动力学方法将人体简化为待研究的关节链,即将关节链的各段近似为刚体,从而将人体简化成为具有有限个自由度的多刚体模型[12-15],因此,可以通过多刚体模型计算出与人体实际运动等效的膝关节的关节力和关节力矩. 基于牛顿动力学和欧拉动力学方程,在全局坐标系下膝关节的关节力(F)和关节力矩(M)涉及的算法以及对应的Matlab的程序(图3)如下所述.

(11)

其中:mF和mS分别为足部、小腿的质量;aF和aS分别为足部、小腿的全局坐标系下的质心加速度;FD为测力台测得的足底反力;MG为足部和小腿的重力产生的力矩;M惯性力为小腿和足部的惯性力产生的惯性力矩;MJ为小腿转动角加速度和小腿质心相对膝关节中心坐标产生的力矩;MD为足底反力产生的力矩.

图3 膝关节力和关节力矩的Matlab程序图Fig.3 Matlab program diagram of knee joint force and knee joint torque

图4 膝关节动力学仿真平台的GUIFig.4 GUI of knee dynamics simulation platform

2.3 平台的GUI设计

人体膝关节动力学仿真平台界面设计以Matlab为编程语言,平台界面主要由数据输入(Import Data)、仿真运行(Simulation)、小腿段运动学参数输出(Shank)、膝关节动力学参数输出(Knee joint)、工具箱(Pipeline)、下肢棍棒行走演示(Begin、Pause、Continue)共6个功能模块构成,如图4所示.

数据输入模块包括手动输入试验对象的身高和体重,试验采样频率以及读取的xls数据单元格位置等. 仿真运行模块主要用于对平台算法的Matlab程序的运行. 小腿运动学参数输出模块包括本平台算法输出的小腿角速度、角加速度曲线图，以及两者与Visual 3D输出的小腿角速度、角加速度对比曲线图. 膝关节动力学参数输出模块包括本平台算法输出的膝关节关节力、关节力矩曲线图，以及两者与Visual 3D输出的膝关节关节力、关节力矩对比曲线图. 工具箱模块包括优化工具箱(Optimtool)和曲线拟合工具箱(Sftool),分别用于对平台算法的优化以及平台输出曲线的拟合. 下肢棍棒行走演示模块主要通过Matlab对下肢在一个完整步态周期内的行走过程的再现.

3 试验结果分析

3.1 平台结果输出及分析

将试验对象的身高与体重、试验采样频率、保存的试验采集的xls数据单元格位置等参数，通过数据输入模块输入到平台内(随机选择多次采集中的5组比较稳定的、完整的右侧下肢步态数据作为本次试验数据),之后运行仿真模块,结束后即可获得试验对象右小腿的角速度、角加速度和右膝关节的关节力、关节力矩的曲线图和数值矩阵. 曲线图则可以直观地反映膝关节正常步态下的关节力和关节力矩的变化,而数值矩阵则可以直接保存使用或用作下一步计算的原始数据.

图5和6分别为本平台和Visual 3D软件在一个运动周期内沿x轴、y轴、z轴3个方向输出的膝关节的关节力、关节力矩的对比曲线图(实线表示本平台的输出曲线,点划线表示Visual 3D软件的输出曲线， 5次试验输出曲线按照线宽由细到粗加以区分).

(a) x轴方向

(b) y轴方向

(c) z轴方向图5 膝关节关节力对比曲线图Fig.5 Contrast curves of knee joint force

(a) x轴方向

(b) y轴方向

(c) z轴方向图6 膝关节力矩对比曲线图Fig.6 Contrast curves of knee joint torque

由图5可以看出,本平台和Visual 3D 5次试验的关节力曲线基本吻合.x轴方向上,运动周期在10%左右时,出现了明显的波谷(图5(a));y轴方向上,一个运动周期内出现了明显的波峰和波谷(图5(b));z轴方向上,一个运动周期内出现了两个波谷(图5(c)). 由图6可以看出,本平台和Visual 3D 5次试验的关节力矩曲线基本吻合.x轴方向上,运动周期在5%左右时,出现了小的波谷,运动周期在20%左右时,出现了明显的波峰(图6(a));y轴方向上,一个运动周期内出现了一次波峰和两次波谷(图6(b));z轴方向上,一个运动周期内出现了两次波峰和一次波谷(图6(c)).

综上所述,膝关节在一个运动周期中,运动周期从 0到 40%为试验对象右腿接触到测力板的过程,因右腿屈伸过程中急剧减速使得膝关节沿3个方向所受关节力和关节力矩出现了明显的波峰和波谷;运动周期从40%到60%为试验对象右腿在测力板上保持平衡的过程,膝关节沿3个方向所受关节力、关节力矩无明显波动;运动周期从60%到100%为试验对象逐渐离开测力板过程,膝关节沿3个方向所受关节力和关节力矩再一次出现了明显的波峰和波谷并直至趋于0.

3.2 平台可靠性验证

目前处理人体运动测量数据较为成熟的C-Motion公司Visual 3D步态/体态分析软件,可以利用运动捕捉系统采集人体运动信息作为原始数据用于人体运动学、动力学参数的计算. 利用SPSS 20.0软件对平台输出的关节力、关节力矩与Visual 3D软件的处理结果进行相关性分析，结果如表2和3所示. 二者的相关程度越高,说明计算结果越趋于一致. 二者的相关程度由相关系数决定,当相关系数取值范围为0.7～1.0,则二者呈高度相关;当相关系数取值范围为0.40～0.69,则二者呈中度相关;当相关系数取值为0～0.39,则二者呈低度相关. 由表2可知,二者输出的关节力在x轴、y轴、z轴3个方向均呈现出高度相关. 由表3可知,二者输出的关节力矩在x轴和y轴方向呈高度相关,而在z轴方向呈现中度相关,这主要因为二者模型具有一定的差异性. 综上所述,就5次试验的相关性和总体趋势而言,二者的结果基本吻合,从而验证了本仿真平台的可靠性.

表2 平台与Visual 3D的关节力在3个坐标轴方向的相关性分析
Table 2 Correlation analysis of joint force on the platform and Visual 3D in three axes direction

变量第1次第2次第3次第4次第5次xyzxyzxyzxyzxyzPearson相关性0.720**0.975**0.999**0.946**0.997**1.00**0.914**0.994**1.00**0.973**0.974**1.00**0.898**0.998**1.00**显著性(双侧)000000000000000N100100100100100100100100100100100100100100100

注: **. 在 0.01 水平(双侧)上显著相关.

表3 平台与Visual 3D的关节力矩在3个坐标轴方向的相关性分析
Table 3 Correlation analysis of joint torque on the platform and Visual 3D in three axes direction

变量第1次第2次第3次第4次第5次xyzxyzxyzxyzxyzPearson相关性0.979**0.689**0.030**0.988**0.982**0.496**0.992**0.964**0.437**0.993**0.984**0.647**0.992**0.971**0.502**显著性(双侧)000.767000000000000N100100100100100100100100100100100100100100100

注: **. 在 0.01 水平(双侧)上显著相关.

3.3 平台结果的实际应用

本平台设计主要应用于医疗领域,如人工膝关节的假体设计、膝关节的康复训练等,其中，针对人工膝关节的假体设计,如通过本平台输出膝关节假体的关节力和关节力矩,则可以为膝关节假体的设计提供依据. 而对于那些已经安装膝关节假体的患者,可以通过采集患者的运动参数输入本平台对关节力和关节力矩的计算,判断其膝关节是否过载. 如果过载较轻,可以通过康复训练对其调节;如果过载较为严重,可以考虑为其更换假体.

4 结 语

本文阐述了基于Matlab所设计的人体膝关节动力学仿真平台的开发过程,并通过对行走步态下膝关节动力学参数即关节力和关节力矩的分析,验证了本平台的可靠性,为人工膝关节置换和膝关节康复运动评定提供了一种新的分析工具. 与目前常用且以C++或Java开发的运动分析软件相比,本平台利用Matlab编程,程序代码比较简洁,没有复杂的建模环节,通过快速自动建模有效地解决了目前常用运动分析软件的人工手动建模、建模速度慢、较高的专业知识要求等带来的不便,同时对不同人体运动数据批量处理和分析的效率提高起着显而易见的作用. 但本平台与目前常用运动分析软件相比,还存在许多不足,如精度相对较低、缺少对人体其他体段运动学和关节动力学参数的分析等,在以后的工作中还需对本平台的功能不断完善以及计算结果做更深一步的分析.

[1] 孙瑜. 人下肢运动链的运动建模仿真分析[J]. 计算机仿真,2015,32(10): 365-368,388.

[2] 于佳彬,郝卫亚,周兴龙. 纵跳落地动作地面反作用力计算机仿真方法的研究[J]. 天津体育学院学报,2013,28(6): 497-501.

[3] 刘书朋,司文,严壮志,等. 基于AnyBodyTM 技术的人体运动建模方法[J]. 生物医学工程学进展,2010,31(3): 131-134.

[4] 郝卫亚. 人体运动的生物力学建模与计算机仿真进展[J]. 医用生物力学,2011,26(2): 97-104.

[5] 唐刚,魏高峰,聂文忠,等. 人体下肢关节坐标系的一种简单定义方法[J]. 北京生物医学工程,2009,28(6): 606-609.

[6] ARNOLD E,WARD S,LIEBER R,et al. A model of the lower limb for analysis of human movement [J]. Annals of Biomedical Engineering,2010,38 (2): 269-279.

[7] HOLZBAUR K R,MURRAY W M,DELP S L. A model of the upper extremity for simulating musculoskeletal surgery and analyzing neuromuscular control [J]. Ann Biomed Eng,2005,33 (6): 829-840.

[8] ALKNER B A ,TESCH P A. Knee extensor and plantar flexor muscle size and function following 90 days of bed rest with or without resistance exercise [J]. Eur J Appl Physiol,2004,93(3): 294-305.

[9] 唐刚,白雪岭,王洪生,等. 基于关节坐标系的人体运动学仿真[J]. 计算机仿真,2011,28(8): 94-97,241.

[10] FUKASHIRO S,HAY D C,NAGANO A. Biomechanical behavior of muscle-tendon complex during dynamic human movements[J]. Journal of Applied Biomechanics,2006,22 (2): 131-147.

[11] TANG G,WANG C T. A muscle-path-plane method for representing muscle contraction during joint movement [J]. Computer Methods in Biomechanics and Biomedical Engineering,2010,14(1): 59-69.

[12] HATZE H. A comprehensive model for human monition simulation and its application to the take-off phase of the long jump [J]. J Biomechanics,1981,14(3): 135-142.

[13] PANDY M G,ZAJAC F E. Optimal muscular coordination strategies for jumping [J]. J Biomechanics,1991,24(1): 1-10.

[14] MILLER D I. A computer simulation model of the airborne phase of diving [D]. Pennsylvania: Pennsylvania State University,1970.

[15] PASSERELLO C E,HUSTON R L. Human attitude control [J]. J Biomechanics,1971,4(2): 95-102.

Simulation Platform Design for Human Knee Joint Dynamics Based on Matlab

TANGGang,CHENGRong-shan,HUXiong,ZHANGBo-yi,WUGang

(School of Logistics Engineering,Shanghai Maritime University,Shanghai 201306,China)

On the basis of the analysis of shank kinematics and knee joint dynamics,a Matlab-based under various kinds of gait of simulation platform for human knee joint dynamics was designed. Basing on inverse dynamics principle to write the Matlab program,graph user interface (GUI) of simulation platform was designed,through reading lower limb movement information from motion capture system and reaction force from 3D load stage,knee joint force and joint torque for various gaits were calculated rapidly. With a normal gait as an example,the right knee joint of a health subject was analyzed,the reliability of the platform was verified by comparing the results with those generated by Visual 3D software,and the platform also provided new analysis tools for the artificial knee joint replacement and knee rehabilitation exercise evaluation.

knee joint dynamics;Matlab;batch read;fast calculation;gait

1671-0444 (2016)05-0725-07

2016-01-19

国家高技术研究发展计划(863)资助项目(2013A2041106);国家自然科学基金资助项目(31300783);中国博士后科学基金资助项目(2014M561458);教育部博士点基金联合资助项目(20123121120004);上海海事大学科研基金资助项目(20130474);上海高校一流学科-管理科学与工程资助项目

唐 刚 (1982—),男,重庆人,副教授,博士,研究方向为人机工程. E-mail: gangtang@shmtu.edu.cn

TP 391

A