基于术后髋骨CT数据重建个性化骨肌多体动力学模型及步态仿真

doi:10.3969/j.issn.2095-4344.2017.35.018

摘要/Abstract

摘要：

文章快速阅读：

文题释义：

骨肌多体动力学模型：利用AnyBody软件建立人类的肌肉骨骼系统，肌肉骨骼系统通常被认为是一个刚体系统，允许应用多体动力学的方法来进行研究。

髋骨应力分布：髋骨中的应力集中主要发生在髋臼上部，包括髋臼上边缘。髋骨后上方应力值均大于后下与前上应力值，由于应力的主要分布在髋骨上方的区域，力的传递工作由皮质骨与松质骨的共同参与，皮质骨的弹性模量较大，所以皮质骨区域比松质骨区域传递更多的力。在正常行走时，力的传递主要由髋臼周围区域的皮质骨承担。由于髋骨后上方区域是主要的力的传递区域，所以髋骨后上方的松质骨与皮质骨均承担力的传递。

摘要

背景：髋骨体内受力复杂，需要更高效、准确的方法进行有限元分析。

目的：构建与患者术后髋骨相对应的个性化骨肌多体动力学模型并进行步态仿真，探究患者术后髋骨的生物力学行为。

方法：根据患者术后髋骨CT数据，在MIMICS中进行三维重建，在Geomagic Studio完成实体化和位置匹配并在Hypermesh中网格化。利用AnyBody软件建立患者个性化骨肌多体动力学模型，模拟正常人平地行走时的动作，导出髋骨在整个运动过程中受到的肌肉力、关节力和关节力矩，将导出的数据作为有限元分析的边界条件，在软件Abauqus上计算步态过程中髋骨受到的应力、应变的大小及其集中部位。

结果与结论：①关节力在2.1 s时，髋关节3个方向关节力的最大值分别为600，2 000，100 N；②髂骨关节几乎只有在Y方向受力，最大值为1 000 N；③髋关节所受最大应力在30%步态周期处约为12 MPa，最大相对位移为0.5 mm；④结果表明，患者术后髋骨所受应力较大的部位主要集中在骶髂关节附近、髋臼后上方、髂骨中央、耻骨上方，术后髋骨所受应力较正常髋骨偏大，患者在康复过程中应注意活动强度。

中国组织工程研究杂志出版内容重点：人工关节；骨植入物；脊柱；骨折；内固定；数字化骨科；组织工程
ORCID: 0000-0002-3776-8000(李兴)

关键词: 骨科植入物, 数字化骨科, 髋骨, 有限元, 骨肌多体动力学, 步态周期, 国家自然科学基金

Abstract:

BACKGROUND: The biomechanical properties of the hip bone are complex in vivo that requires a more efficient and accurate method for the finite element analysis.

OBJECTIVE: To develop a musculoskeletal multi-body dynamic model of the patient-specific hip joint after surgery and simulate gait, and to explore the biomechanics of the hip joint.

METHODS: CT data of the patient’s hip joint postoperatively were collected, three-dimension reconstruction underwent in MIMICS, and materialization and position matched on Geomagic Studio and meshed on Hypermesh. A musculoskeletal multi-body dynamic model of the patient-specific hip joint was established using AnyBody software and simulated the normal movements during walking, then export the muscle force, oint forces, joint moment as the Boundary conditions in finite element analysis. Then, the muscle forces, joint forces and torque were obtained and imported to be the boundary conditions for finite element analysis. The stress and strain values and concentration parts were measured on Abauqus during walking.

RESULTS AND CONCLUSION: (1) The maximum value of the joint force in the three directions of the hip joint was 600, 2 000, 100 N at 2.1 seconds. (2) The ilium force was almost only in the Y direction, and the maximum value was 1 000 N. (3) The maximum stress on the hip joint was about 12 MPa at 30% of the gait cycle and the maximum relative displacement was 0.5 mm. (4) To conclude, the stress of the hip after surgery is mainly concentrated on the vicinity of the sacroiliac joint, postersuperior of the acetabulum, middle region of iliac bone, and upper zone of pubic bone; the stress of the hip is higher than that of the normal hip, and the activity intensity should be paid attention during rehabilitation.

中国组织工程研究杂志出版内容重点：人工关节；骨植入物；脊柱；骨折；内固定；数字化骨科；组织工程

Key words: Hip Fractures, Finite Element Analysis, Kinetics, Gait, Tissue Engineering

中图分类号:

R318

李兴，王洪福，罗伟，成博，王俊元，刘峰. 基于术后髋骨CT数据重建个性化骨肌多体动力学模型及步态仿真[J]. 中国组织工程研究, 2017, 21(35): 5685-5690.

Li Xing, Wang Hong-fu, Luo Wei, Cheng Bo, Wang Jun-yuan, Liu Feng. Musculoskeletal multi-body dynamic simulation on patient-specific hip bone after surgery and gait simulation [J]. Chinese Journal of Tissue Engineering Research, 2017, 21(35): 5685-5690.

图/表（结果） 1

参考文献

[1] 李春蓉,杨莲欢.优质护理对髋骨骨折患者术后并发症发生率及满意度的影响[J]. 中国当代医药,2016,23(11):193-194.

[2] 曾丽. 优质护理对髋骨骨折患者术后并发症发生率及满意度的影响[J]. 中国医药导刊, 2015,17(9):955-956.

[3] 何学安, 罗兵. 数字骨科在骨科临床的研究进展[J]. 中国社区医师, 2016,32(6):19.

[4] 章莹,尹庆水,万磊,等. 数字技术在创伤骨科的应用临床数字骨科(一)[J]. 中国骨科临床与基础研究杂志, 2011, 3(2):113-119.

[5] 李国庆.髋臼骨折内固定数字化模型的构建及有限元分析[D]. 华中科技大学,2013.

[6] Phillips AT,Pankaj P,Howie CR,et al. Finite element modelling of the pelvis: inclusion of muscular and ligamentous boundary conditions. Med Eng Phys. 2007;29(7):739-748.

[7] 黄进成,刘曦明,蔡贤华,等.复杂髋臼骨折内固定术后有限元分析[J]. 中华实验外科杂志,2014,31(7):1454-1456.

[8] 张景僚. 骨盆三维有限元模型的建立及其分析[D]. 第一军医大学, 2007.

[9] 魏帅帅. 髋臼骨折术后个体化有限元分析[D]. 华中科技大学, 2012.

[10] 邱长茂. 髋臼骨折术后股骨头内移对臼顶负重区应力影响的有限元分析[D]. 天津医科大学, 2015.

[11] Chu A, Hughes RE. A method to determine whether a musculoskeletal model can resist arbitrary external loadings within a prescribed range. Comput Methods Biomech Biomed Engin. 2010;13(6):795-802.

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

[13] 杨挺,郑建河,姚子龙,等.行走中股骨生物力学特性的有限元分析[J].广东医学,2016,37(4):512-515.

[14] 吕国敏.基于AnyBody的汽车驾驶员坐姿力学特性建模及坐姿支撑设计[D].山东大学,2016.

[15] 华猛.基于AnyBody生物力学仿真的驾驶姿势舒适机理研究[D].吉林大学,2015.

[16] 刘述芝,胡志刚,张健.冲击载荷作用下运动员下肢动态响应的逆向动力学仿真[J].医用生物力学,2015,30(1):30-37.

[17] 陈瑱贤,王玲,李涤尘,等.全膝关节置换个体化患者右转步态的骨肌多体动力学仿真[J]. 医用生物力学,2015,30(5):397-403.

[18] 单丽君,胡忠安.基于AnyBody的髋关节康复训练肌肉力的分析[J].大连交通大学学报,2014,35(1):50-52.

[19] Damsgaard M, Rasmussen J, Christensen ST, et al. Analysis of musculoskeletal systems in the AnyBody Modeling System. Simulat Modell Pract Theor. 2006;14(8):1100-1111.

[20] 丁光兴.基于个性化的长管骨有限元参数化建模及实验验证[D].南京理工大学,2012.

[21] 魏峰.基于CT图像的人体腰椎有限元模型构建与力学分析[D].南京理工大学,2015.

[22] Janko DJ, Miomir LJ. Finite element modeling of the vertebra with geometry and material properties retrieved form CT-scan data. Mech Eng. 2010;8(1):19-26.

[23] GB 10000-1988 中国成年人人体尺寸[S].

[24] 张迪.基于运动捕捉数据的三维人体运动合成[D].长安大学, 2015.

[25] 赵明光.人工全髋关节置换三维有限元建模及其在体生物力学研究[D].河南科技大学,2013.

[26] 聂涌,马俊,康鹏德,等.正常步态周期中髋臼周围区域的应力分布及其在THA髋臼重建中的指导[J].医用生物力学, 2014,29(1): 31-37.

[27] 荀福兴. 不同材料属性分配方法对椎体有限元模型力学性能的影响[D]. 南方医科大学, 2013.

[28] 张国栋,廖维靖,陶圣祥,等. 股骨有限元分析赋材料属性的方法[J]. 中国组织工程研究,2009,13(43):8436-8441.

[29] 张国栋,廖维靖,陶圣祥,等. 股骨颈有限元分析的赋材料属性方法探讨及有效性验证[J]. 中国组织工程研究, 2009, 13(52): 10263-10268.

[30] Hao Z, Wan C, Gao X, et al. The effect of boundary condition on the biomechanics of a human pelvic joint under an axial compressive load: a three-dimensional finite element model. J Biomech Eng. 2011;133(10):101006.

[31] 刘欣伟,闫寒,刘中洋,等. 包含肌肉、韧带组织的骨盆、髋臼3D有限元模型的构建[J]. 临床军医杂志,2014,42(4):331-335.

[32] Kaku N, Tsumura H, Taira H, et al. Biomechanical study of load transfer of the pubic ramus due to pelvic inclination after hip joint surgery using a three-dimensional finite element model. J Orthop Sci. 2004;9(3):264-269.

[33] Li J, Stewart TD, Jin Z, et al. The influence of size, clearance, cartilage properties, thickness and hemiarthroplasty on the contact mechanics of the hip joint with biphasic layers. J Biomech. 2013;46(10):1641-1647.

[34] Coultrup OJ, Hunt C, Wroblewski BM, et al. Computational assessment of the effect of polyethylene wear rate, mantle thickness, and porosity on the mechanical failure of the acetabular cement mantle. J Orthop Res. 2010;28(5): 565-570.

[35] Shi D , Wang F , Wang D , et al. 3-D ?nite element analysis of the in-?uence of synovial condition in sacroiliac joint on the load transmission in human pelvic system. Med Eng Phys. 2014;36:745-753 .

[36] Cilingir AC, Ucar V, Kazan R. Three-dimensional anatomic finite element modelling of hemi-arthroplasty of human hip joint three-dimensional anatomic finite element modelling of hemi-arthroplasty of human hip joint. Trends Biomat Artific Organs. 2007;21(1):79550B-79550B-8.

[37] Kainz H, Modenese L, Lloyd DG, et al. Joint kinematic calculation based on clinical direct kinematic versus inverse kinematic gait models. J Biomech. 2016;49(9):1658-1669.

引言

髋骨骨折是骨科的常见疾病，在骨折手术治疗中，通常是对骨折部位进行加固或者替换，但由于术后恢复缓慢，严重影响患者的生活质量和健康，进而对病情康复和生活质量造成严重影响^[1-2]。而对患者术后髋骨进行不同场景下的生物力学分析，有利于对患者康复锻炼提供准确建议，规避不必要的风险，并为手术方案的改进提供基础。而有限元法是生物力学分析中常用的技术手段。有限元法通过利用计算机把不规则的、复杂力学分析对象离散化成有限个几何单元体进行分解计算，能够反映机体内部的应力变化情况^[3-4]。由于骨骼、肌肉难以分离且结构复杂，所以为了便于分析求解往往会简化分析模型。髋骨区域结构和力学环境复杂，所以在髋骨往常的有限元研究中边界条件设置得到简化^[5-10]。骨科有限元分析的效果(即分析结果与真实的相似程度)很大程度上依赖于边界条件的设置，所以需要一个更加准确、简便的方法。Anybody软件是人体计算分析的革命性的技术，它可以通过计算机语言的编辑，来模拟人体不同的动作，如行走、奔跑、上下阶梯、抬举重物等，可以分析完整的肌肉与骨骼，所以能真实地模拟人体的生物力学环境^[11-12]，而且有很多学者利用AnyBody软件来进行生物力学分析，并证实了AnyBody软件在生物力学研究中的准确性和科学性^[13-19]。

文章首先利用AnyBody建立与患者术后髋骨相对应的骨肌多体动力学模型，模拟患者术后在平地行走的情况，然后输出AnyBody中的肌肉附着点、约束、肌肉力、骶骨/髋关节力和力矩作为有限元分析的边界条件，利用有限元软件Abaqus对患者术后髋骨进行生物力学分析。

中国组织工程研究杂志出版内容重点：人工关节；骨植入物；脊柱；骨折；内固定；数字化骨科；组织工程

材料方法

1.1 设计 有限元分析实验。

1.2 时间及地点 于2016年9至12月在中北大学完成。

1.3 对象 实验数据来自于成年男性髋骨骨折患者，58岁，身高175 cm，体质量70 kg。在征得患者及家属同意后，按照医疗机构管理条例的相关要求与患者签署知情同意书，报备伦理委员会批准同意

1.4 方法

1.4.1 所用仪器及软件 16层螺旋CT(德国Siemens公司)，Mimics 17.0医学建模软件(比利时Materialise公司)，Geomagic Studio 2013(美国Geomagic)，Hypermesh 14.0软件(美国Altair公司)，Abaqus 6.13软件(美国达索公司)。

1.4.2 髋关节模型的前处理

髋关节模型的建立：将术后骨盆CT的 DICOM格式平扫图像导入进Mimics 17.0，经过软件自动处理后，获得髋骨矢状面及冠状面图片。在Mimics 17.0中建立右髋骨3D模型，并以STL格式导出。将髋骨的三维模型图像以STL格式导入Geomagic 2013中，通过松弛曲面、删除钉状物、填补漏洞、删除岛等操作，进行曲面的光滑处理，完成三角面片的处理，并在Geomagic 2013中构建NURBS曲面片，完成髋骨模型的实体化，并进行与AnyBody中的髋骨模型进行位置匹配，输出IGES文件。

髋关节模型的网格划分：在HyperMesh 14.0选择C3D10M四面体单元对髋骨进行网格划分，由于髋骨上带有螺钉孔，螺钉孔直径为5 mm，所以选择网格的尺寸为 0.5 mm。由于髋臼内侧含有一层皮质骨，在对髋骨进行网格划分时，也需要对皮质骨部分进行划分。在Geomagic 2013中，选择髋臼内表面进行1.5 mm的加厚偏移，以模拟皮质骨并对其划分网格。网格划分完成后，松质骨单元数目为190 700，皮质骨单元数目为32 500。

材料赋值：人体骨骼的各部分材料呈现出异质性，这一特性就需要利用尽可能多的材料数目来准确的对骨骼进行材料赋值，从而来满足骨骼材料的多样性。原则上来说，灰度值区间分的越细以及材料数目越多，建立的材料模型与真实骨骼模型越接近，但是过多的材料数目会加大分析的难度并且计算时间过长。参考丁光兴等^[20-21]关于材料数目对有限元分析的影响，证明了小数目的松质骨以及密质骨材料能满足有限元仿真分析的需要。故本文将松质骨设定为一种材料，皮质骨设定为一种材料。根据Janko等^[22]的文献，髋骨的泊松比选择0.3。皮质骨的弹性模量选择为18 GPa，松质骨的弹性模量为1.37 GPa。在hypermesh 14.0进行完网格划分和材料赋值后，输出名为hip bone的inp文件。

1.4.3 髋骨骨肌多体动力学模型的构建运用骨肌多体动力学软件Any Body，通过修改软件提供的通用模型，构建与患者髋骨相对应的个性化骨肌多体动力学模型。

髋骨位置匹配：因为髋骨模型的参考坐标系要与AnyBody软件MoCap Model中人体模型参考坐标系一致，所以需将AnyBody导出的模型与患者髋骨的模型进行位置匹配，固定AnyBody模型与浮动髋骨模型，进行N点注册，若模型不匹配，则需要在AnyBody对髋骨模型的大小进行修改，直至两个模型位置高度吻合，并建立相一致的坐标系，如图1所示。

骨骼模型建立：人体部件模型主要指肌肉和骨骼模型的建立。骨肌多体动力学中将人体的每块骨骼都she为刚体，不考虑任何的弹塑性变形，每块骨骼都有质心和转动惯量，以模拟骨骼的惯性特征。肌肉模型采用起止点肌肉模型，即肌肉束直线穿过肌肉的起点和终点。骨骼形状及肌肉在骨骼上的附着点根据Horseman发布的解剖数据建立。根据患者的CT数据，并参考根据中国成年人人体尺寸标准^[23]，如表1所示。

关节模型的建立：关节模型主要包括关节之间的连接和其驱动，模型中包含了每块骨骼之间的相对运动关系，骨骼之间的连接用刚性的机械联接代替。研究中不考虑髋关节、膝关节以及踝关节的复杂运动，忽略了关节面之间的分离，只考虑其主要的运动。髋关节的运动分别是屈曲/伸展、内收/外展、内/外旋，用球型铰链模拟其运动；模型中只考虑髋关节的屈曲/伸展，用铰链模拟其运动；同样的，踝关节也只考虑了屈曲/伸展；由于上肢(头部、手臂、躯干)的运动对下肢运动与力学的影响很小，不予考虑。

模型的运动驱动：Anybody模型的驱动可以通过运动方程的定义，也可以运用带有标记点的C3D文件格式进行驱动。C3D驱动文件的获取，首先要对人体进行标记，在Any Body软件的数据库中可以找到标记点离骨头的距离。运动生成的数据会生成一个标准的C3D格式数据文件^[24]。C3D格式是一个公认的二进制格式文件，这些格式文件是步态分析领域记录同步的3D和模拟数据^[25]。AnyBody可以直接从C3D文件中读取数据，运用这些数据继续驱动。模型中标记点以及C3D数据文件已经确定了模型运行的空间位置和时间参数，如图2所示。

1.4.4 有限元分析将前处理完成后名为hip bone的inp文件导入到AnyBody中。患者个性化骨肌多体动力学模型已构建完成，模拟患者在平地行走时的情况，完成一个完整步态周期内的仿真。AnyBody运行完成后，生成包含有髋关节关节力和关节力矩和肌肉力等信息且格格式为.Xml的文件，AnyBody中的AnyFE2Abq.exe插件则可以将有限元模型的网格信息和xml代码进行整合，输出包括肌肉附着点、约束、肌肉力、骶骨/髋关节力和力矩以及有限元模型网格信息的inp文件，将其命名为output。肌肉附着点、约束、肌肉力、骶骨/髋关节力和力矩则在有限元分析中作为髋骨模型的边界条件。将hip bone和output的inp文件同时导入到Abaqus中，Abaqus自动识别inp文件中的信息，建立任务，提交计算，完成髋骨模型的准静态分析。Abaqus读取完inp文件，生成的模型的边界条件如图3所示。

1.5 主要观察指标 基于CT数据建立患者髋骨术后个性化骨肌多体动力学模型，分析患者术后髋骨的生物力学特征。

中国组织工程研究杂志出版内容重点：人工关节；骨植入物；脊柱；骨折；内固定；数字化骨科；组织工程

讨论

文章基于患者术后髋骨的CT数据，成功在MIMICS软件中建立了患者术后髋骨的三维模型，并利用Geomagic Studio和Hypermesh对模型进行处理，然后在AnyBody中建立了患者个性化髋骨骨肌动力学模型，最后通过有限元求解出术后髋骨受到的应力、应变及分布区域。髋骨中的应力集中主要发生在髋臼上部，包括髋臼上边缘。髋骨后上方应力值均大于后下与前上应力值，由于应力的主要分布在髋骨上方的区域，力的传递工作由皮质骨与松质骨的共同参与，皮质骨的弹性模量较大，所以皮质骨区域比松质骨区域传递更多的力，髋臼后上骨表面的应力值大于后下与前上表面应力值。在正常行走时，力的传递主要由髋臼周围区域的皮质骨承担。由于髋骨后上方区域是主要的力的传递区域，所以髋骨后上方的松质骨与皮质骨均承担力的传递。聂涌等^[26]研究的正常步态周期中髋臼周围区域应力最大值接近7 MPa，而本次研究中髋骨区域最大应力值为12 MPa。由此可以看出，本文中得到的髋骨所受应力最大值的结果比聂涌等的研究中髋骨所受应力最大值的结果偏大。

经过研究对比发现两者载荷、有限元材料属性设置以及边界条件均有差异。本文中加载作用在髋骨垂直方向上的关节力达到了2 000 N，而聂涌等加载的关节接触力最大为1 400 N，而且只考虑了髋臼与股骨头之间的接触力，并没有考虑肌肉力和关节力矩的作用，而本文加载的载荷包括了肌肉力，关节力和关节力矩。本文中设置的皮质骨和松质骨的材料属性为固定值术，而聂涌等^[26]是根据髋骨的CT灰度值针对不同单元赋予独立的值。本文中设置的边界条件来自于AnyBody软件，包括肌肉、韧带对髋骨的约束，边界条件设置更为复杂。术后髋骨在步态过程中所受应力比正常时期偏大，所以患者骨折手术后应注意活动强度，尽量避免剧烈活动，在行走过程中可借助拐杖等，以减轻髋骨所受应力。

有限元分析作为分析生物力学的主要辅助技术手段，提高有限元分析的准确性和效率十分有必要。材料属性的设置和边界条件的设置对有限元精度的求解有着重要的影响。通常对材料属性赋值有均匀赋值法、查表赋值法、蒙板赋值3种方法，而人体骨骼是由羟基磷灰石和胶原纤维组成的复合材料，骨的力学性质(弹性模量、剪切模量、黏弹性、破坏时的极限应力和应变等)不仅与复合材料本身性质有关，还与骨的结构组成、形态和胶原纤维蛋白如何组合等有关^[3]。张国栋等^[27-29]认为利用MIMICS软件，根据CT灰度值对骨骼进行赋值从而贴近人骨复杂材料属性的状况，CT灰度幅值比本文所用的手工均匀赋值得出的有限元结果精度要高。

前人对髋骨有限元分析中如何设置准确边界条件的研究一直没有停止^[30-36]。Hao等^[30]构建了三个髋骨有限元模型，分析边界条件对完整髋关节、无骶髂韧带的髋关节和无股骨近端骨盆关节的影响，发现髋骨生物力学有限元分析受到边界条件设置的影响极大。刘欣伟、黄进成和Phillips等^[6-7^，^21]在有限元软件中重建髋骨周围的肌肉和韧带等作为髋骨有限元分析的边界条件，但是在建模过程中重建肌肉和韧带耗时耗力，而且重建点和真实附着点之间有一定的位置出入。本文的边界条件来源于AnyBody软件的输出数据，包括肌肉附着点、约束、肌肉力、骶骨/髋关节力和力矩，可以更好地模拟髋骨在人体所处的真实的生物力学环境。在Kainz、陈瑱贤和刘述芝等^[16-17^，^37]的研究中，经过与临床试验的对比，AnyBody软件仿真结果的准确性和实用性得到肯定。而且通过Anybody 仿真的有限元分析，可以精确、快捷地获得某个模型的边界条件^[13]，利用AnyBody软件和有限元软件联合进行生物力学分析，可以使分析更高效、准确。

本文只研究了患者在日常行走中术后髋骨所受的最大应力和最大相对位移，由于时间关系仿真实验并没有涵盖太多场景，如日常坐姿起立等一些生活场景，并且没有和临床实验结合实为遗憾。但由于计算机技术的发展，仿真方法不啻为一个快捷、有效、经济的选择，为患者康复或治疗过程提供数据支撑，使治疗过程更高效。

本文通过骨肌动力学软件AnyBody与有限元软件Abaqus的联用，得到了髋骨在患者体内所受的应力大小及应力区域，为患者的康复提供积极的指导。而且此方法不需要在有限元软件中重建模型的肌肉和骨骼附着点作为有限元分析的边界条件，节省了仿真分析的时间；模拟环境更贴合人体真实生物力学环境，提高了有限元分析的精度，为生物力学分析提供了新的借鉴。

中国组织工程研究杂志出版内容重点：人工关节；骨植入物；脊柱；骨折；内固定；数字化骨科；组织工程