CN104462723A - 一种基于拓扑优化及骨重建仿真的个性化椎间融合器设计方法 - Google Patents

Abstract

本发明涉及一种基于拓扑优化及骨重建仿真的椎间融合器个性化设计方法。所述方法将椎体进行CT/MRI扫描，然后将椎体CT/MRI连续断层图像进行分割及建立三维模型；通过网格划分及平滑、材料赋值、边界条件设定及力学加载等步骤完成有限元模型的建立；采用变密度法来对椎间融合器进行拓扑优化设计；最后通过基于应变能密度的适应性骨重建理论对椎体及融合器填充骨区域进行骨重建数值仿真模拟，来评价设计是否合理。

Description

一种基于拓扑优化及骨重建仿真的个性化椎间融合器设计方法

技术领域

本发明涉及一种个性化椎间融合器设计方法，更具体地说是涉及一种基于拓扑优化及骨重建仿真评价对椎间融合器进行个性化设计的方法。

背景技术

在临床椎间融合术中，椎间融合器的植入方式、大小、空间形状与神经损伤、植入稳定性、椎间植骨量密切相关。但是，因为患者生理参数和体内环境个性化特征的差异较大，标准化、批量化、序列化生产的椎间融合器与人体的相互作用难以最优化，导致一些椎间融合器进入体内后，发生过度磨损、疲劳失效等。因此，椎间融合器的设计必须是个性化的。

目前，市面上的椎间融合器大多为螺纹状融合器与箱形融合器，其外部均为实体，中间为上下贯通的孔，用以填充移植骨。在椎间融合术中，椎体间融合器一方面要求要有足够的力学强度来维持椎间隙高度，也就是实体部分体积尽可能大，以提供足够的力学支撑；另一方面要求能尽可能多地植入骨便于融合，也就是实体部分体积尽可能小，用于植入骨的孔足够大。为了解决上述矛盾，有必要引入拓扑优化方法对融合器进行优化设计。

对于椎间融合器的有限元仿真评价研究已有很多，现有研究基本均关注于融合器植入后的即时受力状态。然而，即刻状态下力学性能较好的融合器或放置方式在术后长期的恢复过程中所产生的效果不一定较好。实际上融合器植入后会改变原有的生物力学环境，从而诱发骨重建过程，引起植入物周围骨质的变化。很多情况下，力学信号在适当的情况可以刺激骨的形成，骨组织可以展现出愈合能力来恢复功能。因此，为了更好的评价融合器的长期效果，将骨重建理论加入椎间融合器有限元仿真评价中是很有必要的。

发明内容

为了满足目前大多数的椎间融合器的个性化设计需求，本发明公开了一种椎间融合器设计方法，所述方法是通过利用CT/MRI图像三维重建技术，在数字化虚拟环境下还原脊椎的解剖学结构特征，在此基础上运用拓扑优化方法及骨重建仿真进行椎间融合器的个性化设计及评价。

本发明根据实际需要，在所开发的软件AdvancedCageDesigner中导入CT/MRI等图像数据，自动地进行脊柱椎体的分割，建立脊柱的几何模型，自动地进行模型的有限元网格划分，然后采用拓扑优化方法进行椎间融合器的个性化设计，最后采用骨重建仿真方法对所设计的融合器进行分析评价及改进。

本发明公开了一种个性化椎间融合器设计方法，具体包括以下步骤:

1)图像采集及预处理

将椎骨进行CT/MRI扫描，然后将CT/MRI扫描获取的椎骨的原始数据导入软件中并按照图像的编号对图像进行自动排序和判断。然后对导入的CT图像进行中值滤波、高斯滤波及二值化处理。

2)椎体图像分割及三维模型的建立

首先在软件水平面视图区手动挑选一个初始片层来启动分割过程，手工绘 制初始片层的椎体大致轮廓，然后采用水平集方法经过迭代获取初始片层的椎体精确轮廓。此后的分割过程将会以参考图层为起点从上、下两个方向对所有片层进行分割操作。由于椎体CT图像相邻片层目标的变化幅度较小，故可得用已处理片层的分割结果所得轮廓作为相邻片层的初始轮廓，然后运用水平集方法获取所有椎体的精确轮廓。最后，可采用手工方法对所获取的轮廓进行修正。对于分割的结果，采用体绘制方法进行绘制，并在窗口中进行显示。

3)有限元模型的建立

采用基于体素的网格划分方法对椎体进行网格划分。待优化设计的融合器模型可以是市面上已有的融合器模型，也可以是***自动生成的初始融合器模型。***自动生成的初始融合器模型根据椎间隙内纤维环包围的区域形成。对于融合器的网格划分，***自动生成的初始融合器模型采用基于体素的网格划分方法进行网格划分，市面上已有的融合器模型直接导入inp格式的融合器网格模型。

完成有限元网格划分后，设置相应的材料属性及边界条件，施加相应的力学加载，完成有限元模型的建立。

4)拓扑优化设计 

椎间融合器拓扑优化设计的目标是在给定设计体积分数约束下，使应变能最小(即刚度最大)。采用变密度法来进行拓扑优化设计。变密度法的基本思想是人为引进了一种假想的密度可变的材料，其相对密度和弹性模量之间的关系也是假定的，每个单元的相对密度为设计变量，将结构拓扑优化问题转化为材料最优分布设计问题。因此，在设计体积分数为m_f的约束条件下，椎间融合器的拓扑优化设计可如下表示：

目标函数： 

$\min C\left(x\right)=\frac{1}{2}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{d}_i^T{\left(x_i\right)}^3{k}_0{d}_i,$

约束条件： 

$\left\{\begin{array}{c}V\left(x\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{v}_i{x}_i\≤m_f{V}_0\\ {10}^{-3}=x_i^{\min }\≤x_i\≤1\end{array}\right.$

其中x_i为单元相对密度，是自变量，C(x)为应变能，d_i为单元位移，k₀为初始单元刚度矩阵，v_i为单元体积，V₀为初始体积，V(x)为设计后体积，N为单元个数，m_f为体积分数。x_i的值在(0,1)之间，x_i接近0代表该单元需移除，x_i接近1代表需要保留。

建立上述的优化模型后，采用优化准则法求解优化模型。优化准则法是一种间接的优化方法，它不直接优化目标函数，而是基于Kuhn-Tucker条件，通过构造Lagrange函数来形成设计变量的更新方案，通过迭代完成优化模型的求解。优化模型的求解完成后，将x_i接近0的单元移除，实现融合器的优化设计。

5)骨重建仿真评价

骨所处的力学环境发生变化时，骨组织会产生适应性改变，采用基于应变 能密度的适应性骨重建理论对椎体及融合器填充骨区域进行骨重建数值仿真模拟。

在模拟研究中，采用分段函数形式的骨重建控制方程来分析预测骨重建结果，将骨密度的变化定义为骨重建刺激的函数，公式为：

其中Δρ为骨密度增量；B为骨重建速率常数，取经验值1(g/cm3)/MPa；U为骨重建刺激值，即单位质量内的应变能密度，表示为：

$U=\frac{0.5\·\mathrm{\σ}\·\mathrm{\ϵ}}{\mathrm{\ρ}},$

其中ε为von Mises应变张量，σ为von Mises应力张量，ρ为骨密度，s是惰性区域的参考宽度(在惰性区域内，骨生成和骨吸收呈动态平衡，骨密度不发生变化)，取经验值0.1；k为骨重建参考激励，取经验值0.004；Δt为时间增量。

骨重建仿真的流程为：首先从有限元模型中计算出刺激量U即应变能密度，进而判断U与骨生长及骨吸收阈值的关系。判断完成后根据相应重建公式进行骨密度的数值计算。若经过多次迭代最终骨密度值达到收敛状态则终止迭代，不收敛则更新骨的材料属性继续计算直至骨密度达到收敛。骨组织材料属性的更新中，骨密度与弹性模量之间的公式为：

E＝3790·ρ³

其中E单位为MPa，对应ρ的单位为g/cm3。

根据骨重建仿真的结果来确定设计是否合理，对于不合理的设计，返回到步骤4对设计进行修改。最终设计的个性化椎间融合器以STL格式输出。

本发明个性化椎间融合器设计方法优点包括：

1、在设计中全面考虑到融合器在体积尽可能小的情况下刚度最大的设计要求，将拓扑优化技术引入到椎间融合器的个性化设计中。

2、将骨重建理论加入椎间融合器有限元仿真评价中，能更好的评价融合器在椎间融合术中的长期效果；

3、设计方法简便实用，操作时间短，易于临床推广。

附图说明

图1是根据本发明的一个实施例的个性化椎间融合器设计具体实施流程图；

图2是根据本发明的一个实施例的拓扑优化具体实施流程图；

图3是根据本发明的一个实施例的骨重建仿真具体实施流程图；

图4是本发明的一个实施例的初始椎间融合器模型；

图5是本发明的一个实施例的植有椎间融合器的椎体有限元模型；

图6是本发明的一个实施例的拓扑优化后的椎间融合器模型相对密度分布图；

图7是本发明的一个实施例的优化后的椎间融合器模型。

具体实施方式

如图1所示，是根据本发明的一个实施例的个性化椎间融合器设计具体实施流程图，包括以下步骤：

步骤101：图像采集及预处理。

对椎骨进行螺旋CT扫描，获取椎骨CT图像。将获取的椎骨的原始图像数据导入软件中并按照图像的编号对图像进行自动排序和判断。然后对导入的CT图像进行中值滤波、高斯滤波及二值化处理。

步骤102：椎体图像分割及三维模型的建立。

首先在软件水平面视图区手动挑选一个初始片层来启动分割过程，手工绘制初始片层的椎体大致轮廓，然后采用水平集方法经过迭代获取初始片层的椎体精确轮廓。此后的分割过程将会以参考图层为起点从上、下两个方向对所有片层进行分割操作。由于椎体CT图像相邻片层目标的变化幅度较小，故可得用已处理片层的分割结果所得轮廓作为相邻片层的初始轮廓，然后运用水平集方法获取所有椎体的精确轮廓。最后，可采用手工方法对所获取的轮廓进行修正。对于分割的结果，采用体绘制方法进行绘制，并在窗口中进行显示。

步骤103：有限元模型的建立。

采用基于体素的网格划分方法对椎体进行网格划分。首先将各部分的像素用不同的编码表示；然后将上下层相邻的8个像素作为节点连接成1个体素，去除冗余的体素，将剩余的体素作为一个八节点的六面体单元。在模型表面和不同材料界面处的单元的八个节点的像素编码是不一样的，此时采用镜像剖分的方法把每个边界六面体单元分解为5个四面体单元，然后通过移动这些四面体单元的节点来进行界面的平滑；为了减少计算量，对于材料一致的一些单元将根据需要进行合并。

待优化设计的融合器模型可以是市面上已有的融合器模型，也可以是***自动生成的初始融合器模型。***自动生成的初始融合器模型根据椎间隙内纤维环包围的区域形成。如图4所示为本发明的一个实施例的***自动生成的初始椎间融合器模型。对于融合器的网格划分，***自动生成的初始融合器模型采用基于体素的网格划分方法进行网格划分，市面上已有的融合器模型直接导入inp格式的融合器网格模型。

利用单元节点灰度计算出单元的弹性模量，设置模型皮质骨厚度和小关节摩擦系数，假设松质骨材料为各向同性，椎间融合器与椎体之间采用面-面接触。下椎体的下终板所有自由度全部固定，选取上椎体的上终板中心点R，通过在R上施加向下压力模拟椎体承重，施加转矩模拟前屈、后伸、侧弯、扭转。通过上述步骤，完成植有椎间融合器的椎体有限元模型的建立(图5)。

步骤104：拓扑优化设计。

椎间融合器拓扑优化设计的目标是在给定设计体积分数约束下，使应变能最小(即刚度最大)。采用变密度法来建立拓扑优化模型，采用优化准则法求解优化模型。

如图2所示，是根据本发明的一个实施例的拓扑优化具体实施流程图,步骤如下：

a)、首先是调入步骤103所建立的有限元模型，并输入设计的体积分数；

b)、采用体素有限元模型求解常用的算法EBE-PCG(element by element precon-ditioned conjugate gradient)算法进行求解，得到位移矩阵d；

c)、位移矩阵d作为输入，求出应变能变化率dc；

d)、dc作为输入参数，利用优化准则法求出更新后的单元相对密度x；

e)、当x达到要求的精度后，停止循环，结束迭代,否则转到步骤b)继续循环； 

f)、输出优化后单元相对密度x及其分布图。拓扑优化后的椎间融合器模型相对密度分布图的一个结果实例如图6所示。

将优化后单元相对密度接近0的单元删除，得到优化后的椎间融合器模型(图7)。

步骤105：骨重建仿真评价。

采用基于应变能密度的适应性骨重建理论对椎体及融合器填充骨区域进行骨重建数值仿真模拟。在模拟研究中，采用分段函数形式的骨重建控制方程来分析预测骨重建结果。

如图3所示，是根据本发明的一个实施例的骨重建仿真具体实施流程图,步骤如下：首先从有限元模型中计算出刺激量U即应变能密度，进而判断U与骨生长及骨吸收阈值的关系，判断完成后根据相应重建公式进行骨密度的数值计算，若经过多次迭代最终骨密度值达到收敛状态则终止迭代，不收敛则更新骨的材料属性继续计算直至骨密度达到收敛。

根据骨重建仿真的结果来确定设计是否合理，对于不合理的设计，返回到步骤104对设计进行修改。最终设计的个性化椎间融合器以STL格式输出。

实施例 个性化腰椎椎间融合器的设计

将L3-L4腰椎进行螺旋CT扫描，扫描参数:层厚0.63mm，层距0.63mm，球管电压120kV，电流225mAs，分辨率512*512pxl，得到150张腰椎CT图像。将这150张CT数据导入软件中进行中值滤波、高斯滤波及二值化处理。

首先在软件水平面视图区手动挑选一个具有椎体及椎弓的图片作为初始片层来启动分割过程，手工绘制初始片层的椎体大致轮廓，然后采用水平集方法经过迭代获取初始片层的椎体精确轮廓。此后的分割过程将会以参考图层为起点从上、下两个方向对所有片层进行分割操作。对于分割的结果，采用体绘制方法进行绘制，并在窗口中进行显示。

待优化设计的融合器模型采用***根据椎间隙内纤维环包围的区域自动形成(图4)。采用基于体素的网格划分方法对椎体和待优化设计的融合器模型进行 网格划分。设置小关节摩擦系数为0.1，假设松质骨材料为各向同性，椎间融合器与椎体之间采用面-面接触。L4椎体的下终板所有自由度全部固定，选取L3椎体的上终板中心点R，通过在R上施加400N向下压力模拟椎体承重，施加10N-m转矩模拟前屈、后伸、侧弯、扭转。通过上述步骤，完成植有椎间融合器的椎体有限元模型的建立(图5)。

设定体积分数为40％，即优化后的椎间融合器的体积为优化前的40％。采用变密度法来建立拓扑优化模型，采用优化准则法求解优化模型，得到优化后的单元相对密度分布图(图6)。将优化后单元相对密度接近0的单元删除，得到优化后的椎间融合器模型(图7)。采用基于应变能密度的适应性骨重建理论对椎体及优化后的融合器填充骨区域进行骨重建数值仿真模拟，来确定设计是否合理，对于不合理的设计，返回到步骤104对设计进行修改。最终设计的个性化椎间融合器以STL格式输出。

应当理解的是，在以上叙述和说明中对本发明所进行的描述只是说明而非限定性的，且在不脱离如所附权利要求书所限定的本发明的前提下，可以对上述实施例进行各种改变、变形、和/或修正。

附图标号说明

101、图像采集及预处理  102、椎体图像分割及三维模型的建立

103、有限元模型的建立   104、拓扑优化设计

105、骨重建仿真评价。

Claims (5)

1.一种基于拓扑优化及骨重建仿真的个性化椎间融合器设计方法，其特征在于包括以下步骤：

1)图像采集及预处理：将椎骨进行CT/MRI扫描，然后将CT/MRI扫描获取的椎骨原始数据进行图像预处理，从而得到椎体图像；

2)椎体图像分割及三维模型的建立：通过自动或手动分割实现椎体图像的分割，采用体绘制方法完成椎体三维模型的重建及显示；

3)有限元模型的建立：分别对椎体三维模型及椎间融合器模型进行网格划分，设置椎体三维模型各部分的材料属性、边界条件及力学加载，完成植有椎间融合器模型的椎体有限元模型的建立；

4)拓扑优化设计：在给定设计体积分数约束下，使模型应变能最小，完成拓扑优化设计；

5)骨重建仿真评价：采用基于应变能密度的适应性骨重建理论对椎体及融合器填充骨区域进行骨重建数值仿真模拟，并对设计进行评价。

2.根据权利要求1所述的基于拓扑优化及骨重建仿真的个性化椎间融合器设计方法，其特征在于，第1)步所述的图像预处理包括中值滤波、高斯滤波及二值化处理。

3.根据权利要求1所述的基于拓扑优化及骨重建仿真的个性化椎间融合器设计方法，其特征在于，第3)步所述椎体有限元模型的建立中，网格划分的具体步骤为：采用基于体素的网格划分方法对椎体进行网格划分，若选择***自动生成的初始融合器模型，则采用基于体素的网格划分方法进行网格划分；若选择其它已有的融合器模型，则将其直接导入inp格式的网格模型。

4.根据权利要求1所述的基于拓扑优化及骨重建仿真的个性化椎间融合器设计方法，其特征在于，第4)步所述拓扑优化设计具体方式为：采用变密度法建立拓扑优化模型，采用优化准则法求解优化模型。

5.根据权利要求1所述的基于拓扑优化及骨重建仿真的个性化椎间融合器设计方法，其特征在于，第5)步所述骨重建数值仿真模拟采用分段函数形式的骨重建控制方程来预测分析骨重建结果。