基于CT影像的骨质疏松性胸腰椎压缩骨折PKP术后三维有限元力学模型的建立与验证

权祯，张晓刚，，秦大平，，宋敏，，张华，张宏伟，赵希云，王志鹏，马涛，陈钵

1.甘肃中医药大学中医临床学院，甘肃兰州730000；2.甘肃中医药大学附属医院脊柱外科，甘肃兰州730020

前言

骨质疏松性胸腰椎压缩骨折（Osteoporotic Thoracolumbar Compression Fracture,OVCF）是在骨质疏松的前提下，暴力释放后胸腰椎局部应力集中，导致椎体发生骨质疏松性骨折，轻者椎体塌陷，高度丢失，脊柱后凸畸形，重者椎体爆裂骨折，压迫脊髓神经，甚至截瘫［1］。胸腰段在整个脊柱解剖结构中特殊，处于活动度较小的胸椎和活动范围较大的腰椎的应力移行区域，不同模载下的生物力学环境复杂，在脊柱损伤中发生率高。胸腰椎骨折发生生物力学机制核心为“力学失衡”，应用有限元生物力学途径研究压缩骨折发生的力学机制及特点，完善骨质疏松性椎体压缩骨折的防治策略至关重要，有限元分析作为生物力学研究的一种新兴方法，在脊柱骨折生物力学研究中逐渐被重视，本文就围绕OVCF经过经皮椎体后凸成形术（PKP）治疗后三维有限元模型的建立展开，为OVCF内置物椎体强化治疗的生物力学研究提供模型支撑。

有限元分析在20 世纪70年代初逐渐被应用于骨科生物力学研究［2］。1974年，Belytschko 等［3］首次应用有限元方法进行脊柱生物力学分析，使得有限元分析在脊柱生物力学研究中得到应用与拓展［4］。有限元分析是一种离散化的数理求解方法，建立高仿真的生物力学模型，对模型有限单元、自由度、节点添加设定限制，应用计算机软件辅助模拟人体真实受力模式［5］。其核心是将无限个质点构成的自由连续体分割成有限的单元及节点，明确模型不同结构的材料属性和受力特性，计算出每个单元的刚度矩阵，每个小单元的刚度矩阵最后集合成总刚度矩阵，运用数字表达出来［6-7］。人体的生理环境和承载特征极其复杂，传统尸体标本建模很难适应脊柱生物力学研究，对椎体的应力测试仅局限于椎体表面应力，无法测量椎体内部的生物力学变化，而有限元分析弥补了这一缺陷，对椎体内部的力学分布特性进行测试。近年来随着计算机软件的开发及生物力学数理研究的不断深入，有限元分析作为生物力学研究的新兴模式而被广泛应用，其具备模型逼真、精准测量、重复性强等优势，弥补了传统骨生物力学研究的诸多干扰因素和不足，在骨质疏松性椎体压缩骨折、内固定置入治疗及PKP 椎体强化治疗的生物力学研究领域逐渐被深入应用［8］。

1 研究对象

选取甘肃中医药大学附属医院脊柱外科收治的T12 骨质疏松性椎体压缩骨折患者1 例，女性，72 岁，身高165 cm，体质量52 kg，骨密度测量T＜-2.5，且并发T12 椎体压缩骨折，确诊为骨质疏松性T12 椎体压缩骨折。排除T12 椎体压缩骨折以外的脊柱其他病变。本研究通过医院医学伦理委员会审核，患者表示知情，并自愿参加实验。

2 实验设备及软件

本实验选用西门子64 排高分辨率螺旋CT（甘肃中医药大学附属医院放射科提供）；高配置台式计算机辅助建模平台（电脑配置：处理器：锐龙AMD R5 2600X；硬盘：金士顿480 G；8 GB 运行内存；8 GB 显卡）。Mimics 19.0 建模软件（Materialise 公司，鲁汶，比利时）进行原始模型的提取与转化；Geomagic Warp 2017软件（Raindrop公司，马布希尔希，美国）进行特征去除、光顺、曲面实体拟合；Solidworks 2017软件（达索公司，马萨诸塞州，美国）进行零件结构的组装与附属结构生成；Ansys Workbench 17.0 软件（Ansys 公司，宾州匹兹堡，美国）添加材料属性、边界条件、坐标及载荷设定及进行生物力学分析。

3 建模方法

3.1 基于Mimics软件进行原始三维模型的提取与转化

3.1.1 骨质疏松性胸腰段CT影像数据的重构与导入对胸腰段T11椎体上缘至L2椎体下缘进行多序列扫描，层厚0.625 mm，电压140 kV，电流200 mA，获取512×512像素的高清连续性图层数据，提取骨组织窗层，并以DICOM格式保存［9］。将DICOM格式的骨质疏松性胸腰椎（T11-L2）的二维CT数据导入Mimics 19.0软件，利用导入（Import Image）命令，对图像的前后、左右、上下不同方位进行调整与重置。图层位置调整完毕后，转化完毕的影像资料会自动重组显示为断面影像，影像断面序列在Mimics操作界面中分别以横断面、矢状面及冠状面的形式显示。

3.1.2 灰度阈值设定及“骨骼－软组织”边界分割二维层面数据导入成功后，进行阈值的选择与界定，点击灰度阈值（Thresholding）命令，选择并调整CT的灰度阈值，设定合理的窗宽，选择阈值区间为239～1 809 HU。只有选择恰当的阈值区间，才能将骨骼和软组织进行精确分离，最大限度保留模型的原始形态（图1）。实际操作中，原始模型提取过程中由于分割对象主要为骨骼组织，则常规选择系统默认的骨骼（Bone）的灰度阈值。再利用编辑蒙版（Edit Masks）命令对图层进行逐层分离，分别对T11、T12、L1、L2、T12 骨水泥与毗邻软组织进行准确分割。选择恰当的分离工具圆形（Circle）、拉索（Lasso）及方形（Square）等，提取过程中不同的解剖部位选用不同的分离工具，应用上述分离工具准确擦除椎体与相邻组织的衔接点［10］。要求精准细致，尽量保留所需结构，删除冗杂结构，若擦除与保留操作衔接不准，则易导致下一步模型重建的失真或表面缺失。

图1 T12椎体的阈值选择与蒙版分离界面Fig.1 Threshold selection of T12 vertebral body and mask separation interface

3.1.3 渐进化区域增长及三维原始单结构模型提取转化单结构模型提取与转化的前提是渐进化的区域增长（Region Growing），具体而言就是将椎体与毗邻组织结构精确分割，最大限度地保持单结构模型的原始形态（图2）。Mimics软件的3D 视图模型预览（Toggle Mask 3D Preview）命令可对模型椎的分割情况进行即时追踪与预览（图3），分析模型椎是否被准确分离提取。若模型椎与毗邻组织有限接触均可导致区域增长失败，蒙版分离不足则结构粘连，分割过度则模型结构丢失，内部缺损需用腔隙填充（Cavity Fill）命令进行编辑，边缘缺损需用边缘绘制（Draw）命令进行修补［11］。待模型椎区域增长成功后，利用3D 模型计算（Caulatate 3D）命令进行原始模型提取与转化，分别构建T11、T12、L1、L2 及骨水泥的三维原始模型，最后将提取完毕的模型椎以STL 格式储存（图4）。

图2 模型椎的区域增长与转化Fig.2 Regional growth and transformation of model vertebrae

图3 模型椎二维CT平面中三维视图Fig.3 Three-dimensional view of model vertebrae in two-dimensional CT plane

图4 提取完毕后的T12高质量原始模型Fig.4 High-quality original model of T12 after extraction

3.2 基于Geomagic软件的T11-L2节段三维曲面模型优化与构建

从Mimics 软件提取的三维模型表面粗糙，内部冗杂，无法进行曲面片划分、拟合实体等操作，不满足脊柱生物力学分析条件，因此，需对原始三维模型的结构缺陷进行优化处理。

3.2.1 原始三维模型的导入及“多边形”命令下模型的优化处理Mimics 软件中提取的三维原始模型存在结构缺陷，不符合生物力学分析特性，需将原始模型导入Geomagic Warp 2017 软件进行前处理。开始/导入命令将STL 格式的原始模型导入优化程序，开始命令仅可导入一个文件，进行单椎体优化操作时，可以使用此命令，导入命令则可同时导入多个模型椎，便于多模型同时优化［12］。模型导入后，软件会自动推荐网格医生命令，由于模型椎的内部冗杂及边缘缺失，暂不推荐使用网格医生命令进行模型检测（图5）。

图5 原始三维模型椎导入与“多边形”修饰界面Fig.5 Import of original 3D model vertebrae and"polygon"modification interface

模型优化处理中椎体内部冗杂删减、边缘孔弥合是关键，由于在CT解剖断面提取层面的边缘空洞未被完全弥合，部分边缘空洞所形成的表面边界向椎体内部发生“逃逸”，造成椎体内部多余的冗杂结构（图6），限制了曲面实体的拟合。多边形/填充单个孔命令将对模型边缘孔和内部孔进行修补，并通过曲面/斜率追踪命令，进行边缘孔“搭桥”处理（图7），降低模型弥合后的变形系数［13］。椎体内部冗杂及边缘孔修复完毕后，模型椎表面结构仍凹凸不平，呈现为钉状物、凹陷及粗糙特性，可利用多边形/去除特征/删除钉状物/光顺命令进行模型椎表面优化，待模型表面光滑、轮廓清晰、结构完整后进入精确曲面程序拟合实体。

图6 模型椎内部冗杂剔除Fig.6 Elimination of miscellaneous inside model vertebrae

图7 “搭桥”辅助修复边缘孔Fig.7 "Bridging"assisted repair of the edge hole

3.2.2 对优化模型进行精确曲面、构造格栅及拟合实体建模曲面实体拟合的前提是绘制合理的轮廓线，轮廓线引导曲面片的绘制至关重要，规整的轮廓线有利于有限元分析前处理阶段接触关系的设定和求解（图8）。待模型外观修正与真实状态完全一致后进行曲面片分割，使用精确曲面命令中的探测轮廓线工具，此时软件本身会自动根据模型的边界轮廓绘制轮廓线，但是软件本身绘制的轮廓线比较粗糙扭曲，不规整的轮廓线对曲面片的生成造成影响，此处则需要手动编辑轮廓线（图9）。将曲率敏感度降低至20，软件自动探测的轮廓线将消失，需要手动根据椎体表面边界绘制轮廓线，为曲面片的生成提供线性引导（图10）。至生成规整、正确的曲面片后，重画网格（图11）。构造格栅可以检查到曲面分割的缺陷，通过编辑曲面命令可以进行修改，格栅越均匀则曲面表面越光滑，格栅构造的均匀程度将会影响到后期的曲面提取和网格划分等诸多操作（图12）［14］。生成的皮质骨的曲面模型分别保存为STP格式（图13）。松质骨位于皮质骨的内部，通过多边形/偏移命令进行建模，将皮质骨向内部偏移1.00 mm，到达松质骨层面，此时松质骨模型的外观形态及边缘较为锐利，仍需用多边形命令中的特征去除、磨砂、光顺等工具对模型进行处理，待模型前处理接近于椎体松质骨真实形态时，进行网格重画、绘制轮廓线、构造曲面片、编造格栅及拟合实体等顺序性建模（图14～图16），最终构建T11-L2节段的椎体及骨水泥（图17～图18）的拟合实体模型。

图8 探测轮廓线Fig.8 Detecting profile

图9 编辑轮廓线Fig.9 Editing profile

图10 构造曲面片Fig.10 Constructing curved surface

图11 重画网格Fig.11 Redrawing grid

图12 构造格栅Fig.12 Constructing grid

图13 拟合曲面实体Fig.13 Fitting curved surface entity

图14 松质骨轮廓线Fig.14 Contour of cancellous bone

3.3 基于Solidworks 软件的T11-L2 节段模型组装与结构重建

图15 松质骨曲面片构造Fig.15 Construction of cancellous bone curved plate

图16 松质骨构造格栅Fig.16 Structural grid of cancellous bone

图17 松质骨拟合实体Fig.17 Cancellous bone fitting entity

图18 弥散型骨水泥Fig.18 Dispersive bone cement

3.3.1 单椎体零件结构的生成将STP 格式优化修饰后的曲面实体导入Solidworks软件中，进行模型组装及单元结构的添加。利用菜单工具栏中的“文件打开”命令，对已经生成的皮质骨拟合曲面和松质骨拟合曲面及骨水泥逐一打开，分别导入后将生成单一椎体皮质骨、松质骨、骨水泥的零件文件，将生成文件以零件（SLDPRT）格式进行保存，如此反复，将所有整体三维模型的各个部分均以零件（SLDPRT）格式保存，生成T11-L2 节段的皮质骨、松质骨、骨水泥的单椎体零件。

3.3.2 基于“原点重合”理论的装配体零件合成及干涉检查生成的模型结构零件是松散的，具有若干单独“坐标”，需按照CT 三维坐标“原点重合”理论将单零件进行装配，然后将其组装为装配体零件，方便下一步整体操作。点击“新建装配体”命令，以“装配体原点重合原理”对离散的T11、T12、L1、L2 及T12 椎体内部的骨水泥进行结构整体化装配（图19），建立T11-L1 节段包含皮质骨、松质骨、骨水泥的整体模型，并以SLDPRT 格式保存。待整体模型组装完毕后，需要对模型局部接触关系进行干涉检查，尤其是关节突关节的局部接触构建，避免模型曲面生成过程中“边界逃逸”后使得模型局部变形，如果整体模型干涉检查阳性，则会导致模型失真和力学测试失败，需要进一步返回Geomagic 软件中对失真模型局部进行重建，以确保模型的高仿真效应和精准性（图20）。

图19 “离散”的T11-L2单结构零件Fig.19 "Discrete"single-structure part of T11-L2

图20 组装后T11-L2节段整体模型Fig.20 Overall model of T11-L2 after assembly

3.3.3 椎间盘、终板、关节突关节软骨的立体模型构建椎间盘、终板、关节突关节软骨由于其结构特性的限制，柔性结构在CT 骨窗层面的分辨率较低，提取较为困难，要依靠高质量的MRI 数据进行提取。本实验是基于胸腰椎CT 数据的仿真模拟，无法直接提取，需要在Solidworks 实体建模软件中进行重构［15］。

椎间盘的生成以椎体表面网格结构基准面为基础，“三点一面”理论结合草图框制，在椎体表面的网络结构上选择同一平面的不同三点，构建平行于椎体上下表面的基准面，借助于测量命令，绘制椎间盘的平面轮廓，利用“拉伸-凸台”命令，初步构建椎间盘的立体形态，此时椎间盘与上椎体的下表面和下椎体的上表面立体结构重合，不符合脊柱的生理特性，需要对模型局部结构进行微调。“等距曲面”命令可以删除重合的冗杂结构，并且等距曲面刀具距离须设置为零，至此T11-12、T12-L1、L1-2 胸腰段椎间盘的立体结构构建完毕（图21）。

图21 T11-L2节段椎间盘组织Fig.21 Intervertebral disc tissue at T11-L2

关节突关节的构建类似于椎间盘建模，同样需要借助于“基准面-草图-等距曲面”命令，在上位椎体的下关节突和下位椎体的上关节突结构重合部位选择基准面并绘制关节软骨平面结构，再利用“拉伸-凸台”命令，构建关节软骨的立体结构，将等距曲面厚度设置为零，对冗杂结构进行删减，最终建立T11-12、T12-L1、L1-2 胸腰段左右关节突关节的立体结构（图22）。

图22 T11-L2节段关节突关节Fig.22 Facet joint at T11-L2

髓核结构外于椎间盘内部，利用测量工具测量髓核体积，髓核大小占椎间盘体积的50%，重建草图，将草图大小牵拉为接近于髓核体积，利用等距曲面进行切割，生成纤维环和髓核［16］图23。上下软骨终板从立体椎间盘凸台结构中进行分离，软骨终板的厚度约为1.5 mm，利用“等距曲面”命令，调整等距曲面厚度为1.5 mm，以接近于真实解剖结构中软骨终板厚度（图24）。皮质骨外壳与松质骨内核的分离借助于“移动复制实体/组合”命令，对T11、T12、L1、L2松质骨实体及骨水泥团块进行复制，“组合”命令实质是将皮质骨与内部的松质骨、骨水泥重合部分进行删减，得到单纯的皮质骨外壳和内部的松质骨、骨水泥团块，以提高模型仿真度。

图23 T11-L2节段髓核、纤维环Fig.23 Nucleus pulposus and annulus fibrosus at T11-L2

图24 T11-L2节段软骨终板、骨水泥结构Fig.24 Cartilage endplate and bone cement structure at T11-L2

3.4 基于Ansys软件构建有限元分析前处理模型

3.4.1 弹性模量和泊松比的设定弹性模量和泊松比是表示在纵向压缩和横向拉伸下材料结构强度和刚度的物理量（表1），就椎体而言，弹性模量和泊松比越大，则材料强度及刚度越强，相反则越低（图25）。骨质疏松模型的材料属性均低于正常骨结构，因此将骨质疏松模型松质骨的所有骨结构降低67%，将皮质骨、终板和后部的弹性模量降低34%［17］。

表1 三维有限元模型材料属性Tab.1 Material properties of three-dimensional finite element model

图25 模型结构的材料参数添加视图Fig.25 Structural material properties addition view

3.4.2 模型坐标系的重置与接触关系设置坐标系为力的加载与释放提供了约束，通过坐标系（construc-tion geometry）命令，可以对模型的坐标进行修正与重置，以适应人体真实静力测试的要求，Z 轴设定为垂直应力加载，模拟垂直压缩受力，同时模型围绕垂直力模拟左右旋转受力，Y轴设定为水平向后的应力加载，模拟前屈、后伸受力，后伸时为正向力，前屈时为反向力，X 轴设定为水平向左力的加载，模拟左右侧弯时受力载荷，左侧弯曲时为正向力，右侧弯曲时为反向力［18］。

合理的结构间接触是应力传递的必要前提，其接触关系主要有固定接触和有限接触，“面-面接触”、“点-面接触”是模型接触的主要形式，可以通过接触（connection）命令对接触关系进行设定，将椎体-终板、皮质骨-松质骨、终板-椎间盘设定为绑定（bonded）模式，关节软骨具有一定的活动度但又不会脱离活动轨迹，其接触为一定无摩擦位移滑动范围内的有限接触，故将其接触模式设定为不绑定的（no sperarete）模式。

3.4.3 模型网格化处理模型网格化将结构精确化，划分为有限个利于求解和控制的单元和节点（表2），不同结构的差异性与复杂性导致模型网格化也存在差异，适当的网格大小影响着模型求解精细度，利用Generate Mesh（网格生成）命令进行网格重置，皮质骨、松质骨结构网格大小设定为2.0 mm，椎间盘、关节软骨、终板、骨水泥等附属结构网格大小设定为1.5 mm。

表2 结构的节点与单元Tab.2 Nodes and units of different structures

3.5 弹簧韧带的添加及刚度系数设定

腰椎韧带结构复杂，在CT 骨窗层面难以提取，需要在Ansys17.0 软件中用弹簧韧带对生理韧带结构进行替代，利用弹簧韧带（spring elements）命令进行添加，并将其生理特性设置为仅拉伸（tensiononly）模式，依次添加前纵韧带、后纵韧带、黄韧带、棘间韧带、棘上韧带、横突间韧带等，并将其模拟特性设置为线性弹性材料。

3.6 力学加载及边界条件设定

人体脊柱的整体运动状态及机制是复杂的，我们无法在软件中进行真实还原，高仿真下模拟不同载荷求解工况，在围绕Z 轴在T11 椎体表面施加垂直向下500 N 的压力，以模拟压缩载荷（图26），-7.5 N·m的反向扭力，模拟前屈载荷（图27），围绕Y轴施加水平向后的7.5 N·m的正向扭力，模拟后伸载荷（图28）；X 轴设定为水平向左的7.5 N·m 扭力和向右的-7.5 N·m扭力，模拟左右侧弯时受力载荷（图29～图30）；围绕Z轴施加7.5 N·m左旋力和-7.5 N·m的反向右旋力，模拟左右侧弯载荷（图31～图32）。为了使其仿真度达到最优，需要对T11-L2节段的运动轨迹进行约束，实验中固定了L2椎体下表面及下关节区域（图33），限制其自由度［19］。

图26 轴向受力加载Fig.26 Axial force-loading

图27 前屈受力加载Fig.27 Anterior flexion force-loading

图28 后伸受力加载Fig.28 Posterior extension force-loading

图29 左侧弯受力加载Fig.29 Left flexion force-loading

图30 右侧弯受力加载Fig.30 Right flexion force-loading

图31 右旋转受力加载Fig.31 Right rotation force-loading

图32 左旋转受力加载Fig.32 Left rotation force-loading

图33 L2椎体约束自由度Fig.33 Constraint freedom of L2 vertebral body

4 OVCF有限元模型的有效性验证

模型建立成功后，利用有限元软件进行计算，在进行生物力学实验之前，对三维有限元模型进行生物力学验证是必要前提，以验证模型的精确性、可靠性。脊柱有限元模型有效性验证的方法主要有直接验证和间接验证，直接验证需要与离体模型进行验证，可靠性及难度较大。间接验证是将有限元模型的角位移与既往生物力学实验进行对比验证，本模型采用了精确化建模路径，精细划分网格及单元，模型仿真度高，准确性强，有效性验证合格，可进行骨质疏松性胸腰椎压缩骨折PKP术后生物力学测试［20］。记录模型在各个运动方向上的最大角位移即活动范围（Range of Motion,ROM）。之后参照标准脊柱有限元实验和生物力学实验的结果，对比所构建模型的屈曲度变化判断所构建模型的有效性，评判其是否可以适用于进一步的模型处理和有限元分析（表3）。

表3 T11-L2节段模型ROM与既往数据T11-L2胸腰段ROM活动度均值（°）Tab.3 Average range of motion(°)of T11-L2 between model and previous data

5 结语

随着计算机应用技术和有限元分析软件的开发与升级，有限元作为一种生物力学分析的有效方法被深入应用于骨科生物力学研究当中，通过仿真模拟途径揭示骨科诸多现象的生物力学机制，就骨质疏松性椎体压缩骨折而言，其在椎体压缩骨折及慢性疼痛的生物力学机制应用广泛［23］。以往脊柱骨折的生物力学研究以尸体建模为主，通过机械压力机测试材料生物力学性能，模型单一、可重复性差，与人体真实复杂的生物力学受力情况差异性大，而有限元分析作为一种数据仿真模拟，具有可重复性、仿真度高、精细性强等优势。

骨质疏松性椎体压缩骨折的生物力学模型仍需要优化升级，目前建模数据是基于DICOM格式的CT数据的建模，对CT 数据质量要求比较高，高精度的三维有限元模型必须以多序列、薄层、高分辨率的CT数据为基础，但只能对骨窗组织进行提取与转化，柔性结构很难被提取，椎间盘、软骨终板、关节软骨及韧带等组织在原始模型中依然缺失，需要在Solidworks软件中进行附属结构的重建及Ansys软件中进行韧带重构。随着生物力学测试的精准度不断提高，基于3.0T 高精度MRI 数据的模型提取与建立将是有限元建模的热点与难点，由于核磁对软组织的识别度比较高，不同组织在MRI 序列中均有不同信号的显影成像，依托于腰椎的MRI 数据有望直接提取椎间盘、韧带、关节软骨、软骨终板等柔性结构，部分学者对基于MRI 数据建模进行了探索［24］，但腰椎韧带、肌肉结构混杂，单纯利用腰椎MRI 数据提取腰椎包括韧带在内的整体模型案例较少，后期随着有限元技术的开发与应用，为提高模型仿真度，利用核磁构建腰椎三维整体模型已成为必然趋势。