脊柱腰段解剖结构复杂，负重较大且运动灵活，是临床病变的好发部位。其中腰椎间盘承担椎体间力的传导，并完成脊柱各向弯曲及扭转运动，易发生退变及损伤。通过分析脊柱腰段的生物力学特性，可以为临床分析损伤原因及疗效评估提供理论指导[1-2]。现有的脊柱腰段有限元模型研究重点集中于模型的结构和材料属性的仿真度，而对运动载荷多通过简单的力矩或力偶模拟，对腰椎稳定性骨骼肌和运动性骨骼肌的实际作用考虑较少[3]。本研究利用CT及MRI断层数据建立正常脊柱腰段三维有限元模型，施加载荷模拟腰椎屈伸、侧弯和扭转运动，旨在为脊柱腰段损伤及病变的发病机制、手术方案及疗效评价等建立数字化研究模型。
1 资料与方法
1.1 数据采集 选取1位男性健康青年志愿者，28岁，身高178 cm，体质量68 kg，既往身体健康，无脊柱相关性疾病及外伤史，行腰椎X射线正侧位片和CT预扫描排除脊柱腰段的器质性病变。实验前将实验的相关内容告知志愿者，并征得其同意。
志愿者取仰卧位体位，腰部放松。首先行螺旋CT扫描(GE公司，Light speed 16排螺旋CT机)，最终得到了247幅断层扫描二维图像。以相同体位行MRI扫描(Philips公司，Achieva 3.0T MRI机)，最终获得到57幅T2WI序列矢状面二维扫描断层图像。所有断层图像数据均以DICOM格式存储。
1.2 脊柱腰段三维有限元模型构建 ①椎骨及髓核模型：用Mimics 15.0软件打开CT数据文件，提取椎骨轮廓，自动生成蒙板；对蒙板进行手动修整，填补空洞、去除无用部分，计算生成各椎骨3D模型，导入3-Matic软件中，进行光滑处理，最终得到结构完整且表面光滑的椎骨模型；以相同的方法基于MRI数据生成椎间盘髓核的3D模型。②椎间盘及韧带模型：将椎骨及髓核3D模型以STL格式输入到3-Matic中利用坐标系进行对位组装；用选择工具确定相邻椎骨的下表面和上表面，翻转法线，使用fix工具连接2个表面创建椎间盘整体外形；根据实际解剖结构在椎体上下表面利用Move surface生成终板，最后通过椎间盘整体、终板和髓核之间进行布尔运算生成纤维环；根据脊柱腰段各韧带的解剖学参数，确定前、后纵韧带在脊柱腰段表面的范围，通过Move surface生成前韧带、后纵韧带模型[4]；以相同方法生成黄韧带、棘间韧带、棘上韧带、横突间韧带、关节囊韧带等结构；在3-Matic软件中对各3D模型组件进行光滑处理，然后进行面网格和体网格划分，以*.cdb格式输出保存。
将cdb文件导入Ansys 15.0，在Finite Element Modeler分割表面后，利用Static Structural模块中的Engineering Date定义材料属性，通过弹性模量和泊松比设置骨密质、骨松质、韧带等各部分参数，各部分单元数目及材料参数见表1。
表1 有限元模型单元数量和材料属性
Table 1 Element number and material parameters of
finite element model
材料 单元数量杨氏模量(E/MPa)泊松比(μ)骨密质 62 74512 0000.3骨松质 103 7451000.2终板 10 0445000.4髓核 14 35810.5纤维环 53 462 4.20.5前纵韧带 17 721 7.80.3后纵韧带 11 661 100.3黄韧带 6 772 150.3棘间韧带 5 581 100.3棘上韧带 2 292 100.3横突间韧带4 592 100.3关节囊韧带46 191 100.3
设定模型各组件之间的接触关系，关节突关节面之间设定为No Separation，其他位置均设定为Bond。在骶骨下面添加Fixed Support，限制各向自由度。至此，人体脊柱腰段的三维有限元模型建立完成。
1.3 负载与骨骼肌附着点 脊柱的负重载荷加载部位确定为椎体上表面和上关节突关节面，以垂直载荷300 N模拟身体上部的重力；在椎骨表面根据解剖结构确定脊柱腰段主要稳定性骨骼肌的附着点。左右稳定肌包括横突间肌、回旋肌，前后稳定肌包括棘间肌、多裂肌。
1.4 模型有效性验证 在椎间盘损伤最为好发的L4～L5节段上进行模型有效性验证。约束L5椎体下面和下关节突关节面，限制其所有的自由度，在L4椎体上表面及关节突关节面上分别给予500 N、1 000 N、1 500 N、2 000 N的轴向压缩载荷[5]，测量轴向位移。在L4椎体上表面中部指定一节点，并在节点上施加300 N的垂直载荷和10 N·m的力矩[6]，分别模拟腰椎屈伸、侧弯和扭转，并测量不同运动工况下腰椎的角位移范围，观察椎间盘应力的分布趋势及特点，将结果与其他研究进行比较，从而验证模型的有效性。
2 结 果
2.1 脊柱腰段有限元模型的结构 利用CT及MRI数据成功建立了外形结构真实准确、生物力学仿真度高的正常脊柱腰段的三维有限元模型，包括L1～L5椎骨及骶骨上部的皮质和髓质、纤维环、髓核、终板、前纵韧带、后纵韧带、黄韧带、棘间韧带、棘上韧带、关节突关节囊等重要结构，总节点数为104 190个、总单元数为339 165个。
2.2 模型有效性 后处理计算结果显示，有限元模型L4～L5节段在500 N、1 000 N、1 500 N、2 000 N的轴向压力下，有限元模型L4轴向位移分别是0.25 mm、0.54 mm、0.81 mm、1.04 mm，描绘出轴向压力-位移曲线图并与文献数据比较，该结果与在相同条件下离体实验和有限元分析的结果相近。模型在正常受力状态下表现出弹性特性。测得不同运动工况下腰椎的角位移平均值数据分别为前屈4.12 °、后伸2.83 °、侧弯3.71 °、扭转1.64 °，与文献中标本和有限元实验测量值相近。
2.3 脊柱腰段L4～L5腰椎椎间盘的应力分析 为了进一步验证模型的有效性，对脊柱腰段有限元模型L5椎体和下关节突关节面限制其所有的自由度，在L4椎体上表面中部指定一节点，并在节点上施加300 N的垂直载荷和10 N·m的力矩，模拟了屈伸、侧弯、扭转等运动工况，观察椎间盘在不同工况下的受力情况。在运动工况下，应力主要集中于腰椎间盘边缘，前屈、后伸时，应力分别集中于椎间盘前、后两侧；侧弯时，椎间盘受压一侧存在着明显的应力集中，且向椎间盘中心有逐渐减小的趋势；扭转时，纤维环受到扭力后发生倾斜至牵张，且应力集中于轴向扭转方向的侧后方。这与人体脊柱腰段的生理特性相符。
3 讨 论
随着计算机技术和有限元软件的发展，利用三维有限元方法研究人体脊柱腰段生物力学特点的研究不断深入[7-8]。三维有限元法具有简便、快速、经济的特点，可以模拟并运算复杂条件下各种材料的力学特点，且实验具有易调整和可重复等特点，可模拟脊柱等复杂结构，使其在腰椎的力学研究中得到迅速的推广。构建胸腰椎三维有限元模型，是脊柱生物力学研究的有效手段；利用CT扫描图像建立模型，从多角度、不同方法可验证三维有限元模型的准确和实用性。吴小辉等[9]研究发现，三维有限元模型能很好地评估人体胸腰椎的受力状况；同时为脊柱内固定系统的稳定性提供理论依据。姜伟等[10]通过扫描健康成人腰椎体建立了L3～L5关节突关节未融合和融合的有限元模型，未融合模型L3/4、L4/5节段活动度与既往文献中腰椎活动度趋势一致。众多学者在构建有限元模型时主要利用CT扫描数据建立椎骨模型，再通过软件按一般解剖结构特点手动建立椎间盘、韧带等软组织[11-13]。这样的模型中软组织较真实情况有一定差距。而MRI对软组织有较好的分辨率，有学者利用MRI数据建立单独的椎间盘结构，但MRI数据对椎骨建模有欠缺。
本研究建立完整而准确的脊柱腰段三维有限元模型，该模型包含椎骨、椎间盘和韧带等结构。本研究总结了以往建模的不足之处，对建模进行了优化，椎骨部分利用CT数据进行构建，这是由于CT对复杂形态和各种密度的组织均有较高的分辨率，尤其对高密度的骨组织成像清晰准确，确保构建出的椎骨模型精准。椎间盘部分利用MRI数据进行建模，由于MRI对软组织显示良好，提高了椎间盘模型的精确性。同时为了保证通过CT和MRI 2组数据建立模型的吻合度，在数据采集和处理上采取了以下措施：①2组数据均采用DICOM格式读取，保证了不同设备之间的兼容性，单位及坐标系参数统一，距离数据等价；②2组数据均采集自同一志愿者，采集时采用相同体位；③在Mimics软件中利用MRI T2序列数据建立髓核模型，这是由于该序列对含水量多的组织显示清晰，髓核边界明显，伪影相对较少，利用2组数据建立的椎骨上下表面重合对2种模型进行吻合，确定髓核的空间位置；④将建立的所有模型导入Mimics中，平滑、组装、检查模型质量后，利用软件中自带的3-Matic组件对模型进行面网格和体网格的划分。这一方法比直接在Ansys软件前处理单元划分网格的效率明显提高，同时保证了有限元模型的网格质量。
脊柱的载荷来源主要有2个途径，分别是负重和骨骼肌的牵拉。现有的模型在生物力学模拟时，载荷形式单一，不能很好地反映脊柱受力的真实情况。有学者尝试探讨不同的肌力方向对人体有限元模型预测结果的影响[4]。但绝大部分的标本和有限元生物力学实验均没有反映脊柱腰段骨骼肌的影响。直立时，躯干重力线经过L4椎体中心腹侧，故脊柱常处于一种持续向前弯曲的运动状态。通过背侧肌的力量和韧带的牵拉来对抗维持脊柱平衡。当人体处于不同体位时，椎间盘内压力的变化除与负重有关外，同椎旁肌的牵拉关系密切。放松直立位时L3～L4椎间盘上的负荷约为测量平面以上体重的2倍。无负荷状态下椎间盘内存在大约10 N/cm2的内压力，如果单纯给予上半身体重载荷，会引起腰椎模型发生前倾，而引入稳定性骨骼肌载荷后，很好地维持了腰椎模型的原有姿态，更接近人体腰椎的正常状态。腰椎的姿态维持和运动主要依靠周围肌的作用力实现，但这些力的大小和方向是在动态调节的，单纯力的加载很难模拟真实状态。只能在脊柱静态稳定状态，在载荷中引入稳定肌力因素。
椎间盘是脊柱功能重要的载荷中心和缓冲结构，在脊柱运动、承载、传递各种载荷中具有关键作用[14]，也是腰椎疾病的好发部位。腰椎间盘应力云图结果显示，在前屈、后伸、侧弯和扭转4种工况下，应力主要集中在椎间盘边缘，且受压侧应力较集中，并向受拉侧扩散消释，这与腰椎的生物力学特性相符。本研究结果发现腰椎椎间盘在轴向、前屈、后伸及侧弯4种工况下，纤维环形变较大，且出现明显的应力集中。当扭转时，纤维环应力主要集中于侧后方，这也间接证实了椎间盘纤维环受力过大导致的纤维环破裂、髓核脱出即腰椎间盘突出，常发生在后外侧，是人体腰椎疾患的重要病因之一。
本研究成功构建了脊柱腰段的三维有限元模型，且有效模拟了腰椎的生物力学特性，但也存在一些不足，需要进一步完善和改进。首先，有限元分析受到了各种客观因素的限制，如骨骼肌、肌腱、韧带的非线性特点难以实现准确的模拟。本研究中虽然引入了肌力因素，但由于脊柱活动时，肌力的大小和方向是在动态变化的，在进行有限元分析时，对骨骼肌等结构的生物力学特性进行了线性化简化处理，故结果和真实情况还存在一定误差，希望在后续研究中通过编写载荷参数曲线来进行模拟。其次，三维有限元模型各项参数的处理，受人为因素的影响，最终会出现运算数据值和量的差异。当有限元运算数据量增加时可能引起误差的增大，从而影响有限元方法分析结果的可靠性。最后，由于数据获取方便、无创性和可重复等特点，目前基于三维有限元模型进行生物力学研究较多。基于尸体标本测试获得的生物力学数据更接近真实情况，有较高的临床可靠性，但标本不易获取且不能重复实验。因此，将2种方法结合使用、相互验证可以提高其结果的可信度，能为各类腰椎疾病的手术设计和评估提供有效参考。
[1] 刘治华,许伟超,徐新伟,等.腰椎牵引角度有限元分析及优化[J].郑州大学学报：医学版,2015,50(4):507-511.
[2] Schollum M,Wade K,Robertson P,et al. A microstructural investigation of disc disruption induced by low frequency cyclic loading[J]. Spine(Phila Pa1976),2018,43(3):E132-142.
[3] 冯其金,赵玲娟,谷福顺,等.三维有限元动物椎体模型的建立及应力分析[J].中国中西医结合外科杂志,2018,24(1):63-68.
[4] Zhu R,Niu WX. The effect of muscle direction on the predictions of finite element model of human lumbar spine[J]. Biomed Res Int,2018,2018(58):1-6.
[5] 苏春涛,眭承志.全腰椎三维有限元模型的构建及仰卧位屈曲模式下椎间盘蠕变的生物力学分析[J].中国康复医学杂志,2016,31(10):1136-1138.
[6] 王宏卫,刘新宇,万熠.人体腰椎L4～L5段有限元模型建立及力学有效性验证[J].医学与哲学,2017,38(10):50-53.
[7] Wang L,Zhang B,Chen S,et al. A Validated finite element analysis of facet joint stress in degenerative lumbar scoliosis[J]. World Neurosurg,2016,95(11):126-133.
[8] Newcomb AG,Baek S,Kelly BP,et al. Effect of screw position on load transfer in lumbar pedicle screws:a non-idealized finite element analysis[J]. Comput Methods Biomech Biomed Engin,2017,20(2):182-192.
[9] 吴小辉,刘小聪,陶平,等.人体胸腰段三维有限元力学的仿真研究[J].微创医学,2017,12(3):315-319.
[10] 姜伟,李威,袁峰,等.L4/5关节突关节融合后椎间盘应力变化的三维有限元分析[J].中国脊柱脊髓杂志,2017,27(5):441-448.
[11] Du CF,Yang N, Guo JC,et al. Biomechanical response of lumbarfacet joints under follower preload:a fnite element study[J]. BMC Musculoskelet Disord,2016,17(1):126.
[12] Claeson AA,Barocas VH. Computer simulation of lumbar flexion shows shear of the facet capsular ligament[J]. Spine J,2017,17(1):109-119.
[13] 刘治华,许伟超,张新民,等.颈椎C2～7三维有限元模型的建立与最优角度牵引仿真研究[J].郑州大学学报：医学版,2016,51(3):359-363.
[14] Zhu Q,Gao X,Levene HB,et al. Influences of nutrition supply and pathways on the degenerative patterns in human intervertebral disc[J]. Spine(Phila Pa1976),2016,41(7):568-576.