[发明专利]一种基于组织分化的长骨骨折愈合仿真系统

权利要求书

1.一种基于组织分化的长骨骨折愈合仿真系统，其特征在于，所述系统包括：

骨折区域几何建模模块(1)、骨折区域生物力学有限元分析模块(2)、骨痂单元组织分化模块(3)和程序终止判断模块(4)；

骨折区域几何建模模块(1)用于根据导入的二维断层扫描图像数据，经过图像预处理后进行骨折部位的三维表面几何模型的建立；

所述的骨折区域几何建模模块(1)实现其功能的具体过程为：

采用基于分割的三维医学影像表面重建算法对图像进行三维表面重构，通过阈值筛选、交互式分割和三维重建过程得到三维表面几何模型；

所述的影像由影像设备CT得到，数据存储格式为DICOM；

骨折区域生物力学有限元分析模块(2)用于对建立好的骨折区域模型进行网格划分，施加外部载荷和设置边界条件；

骨折区域生物力学有限元分析模块(2)还用于初始骨折区域环境的设置；

骨折区域生物力学有限元分析模块(2)还用于计算单元力学刺激；

所述的骨折区域生物力学有限元分析模块(2)实现其功能的具体过程为：

1)将骨折区域三维表面几何模型进行网格化分，使连续的几何模型离散化，得到骨折区域有限元模型；

所述的网格划分包括面网格划分和体网格划分两个步骤；面网格划分过程用于将三维表面模型进行优化，包括：表面模型优化，平滑处理，修补漏洞；表面模型的优化通过减小表面模型的三角面片来实现，该过程只需将相邻的两个顶点合并到一个新的顶点上，并延续原有的拓扑关系；平滑处理的过程中，对三维的面网格模型进行去噪；修补漏洞的过程中，通过将模型当中的空洞提取成空间多边形，然后对空洞多边形进行三角化的方法实现；体网格划分的过程是将面网格模型进行拉伸、旋转步骤来实现的；

通过网格划分得到的骨折区域有限元模型包括单元编号和节点坐标两部分；

节点坐标包含三列数据，三列数据分别代表每个节点的空间坐标值；

单元编号包含四列数据，四列数据分别为每个单元的四个节点的节点序号；

2)在有限元模型上施加外加载荷，并设置边界条件；载荷的大小由骨所承受的力的大小决定，实验对象不同，所受的力也不同；

3)对骨折区域有限元模型进行骨折区域初始环境设置；骨折区域由皮质骨和骨痂区域两部分组成；初始骨折区域环境设置包括皮质骨材料属性赋值，皮质骨血供赋值，初始骨痂材料属性赋值，初始骨痂区域血供赋值；

距离骨折断端5mm内的皮质骨血供设置为0％，其余部分皮质骨血供设置为100％；

骨痂外周3mm内血供设置为30％，骨痂内部血供设置为0％；

4)将骨折区域看作双相多孔弹性模型，由多孔弹性理论得到骨痂单元的本构方程，平衡方程和几何方程，并通过有限单元法计算骨痂单元应力刺激S，具体过程为：

a.本构方程

σrrσθθσzzτrθτθzτrz=M11M12M13000M12M11M13000M13M13M33000000M44000000M55000000M55ϵrrϵθθϵzzγrθγθzγrz-ααα′000p---(1)]]>

式中，σ_rr，σ_θθ，σ_zz为正应力，τ_rθ，τ_θz，τ_rz为剪应力；ε_rr，ε_θθ，ε_zz为正应变，γ_rθ，γ_θz，γ_rz为主应变；α，α'分别为各向同性弹性面的Biot系数和轴向Biot系数；p为骨痂单元中的流体压力；M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅分别为脱水的弹性模量矩阵分量；

其中，M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅脱水的弹性模量矩阵分量表达式如下所示：

M11=Er(Er-Ervz2)(1+vr)-1(Ez-Ezvr-2Ervz2)-1---(2)]]>

M12=Er(Ezvr-Ervz2)(1+vr)-1(Ez-Ezvr-2Ervz2)-1---(3)]]>

M13=ErEzvz(Ez-Ezvr-2Ervz2)-1---(4)]]>

M33=Ez2(1-vr)(Ez-Ezvr-2Ervz2)-1---(5)]]>

M₄₄＝E_r/2(1+ν_r) (6)

M₅₅＝G'(7)

式中，E_r，ν_r分别是各同性弹性层的弹性模量和泊松比；E_z，ν_z分别是轴向弹性模量和泊松比；G'为剪切模量；

b.平衡方程

∂σrr∂r+1r∂τrθ∂θ+1r∂τrz∂z+σrr-σθθr=0---(8)]]>

∂τrθ∂r+1r∂σθθ∂θ+1r∂τθz∂z+2τθzr=0---(9)]]>

∂τrz∂r+1r∂τθz∂θ+1r∂σzz∂z+τrzr=0---(10)]]>

式中，σ_rr，σ_θθ，σ_zz为正应力，τ_rθ，τ_θz，τ_rz为剪应力；r为径向半径；

c.几何方程

ϵrr=∂ur∂r---(11)]]>

ϵθθ=1r∂uθ∂θ+urr---(12)]]>

ϵzz=∂uz∂z---(13)]]>

γrθ=1r∂ur∂θ+∂uθ∂r-uθr---(14)]]>

γθz=1r∂uz∂θ---(15)]]>

γrz=∂uz∂r---(16)]]>

式中，ε_rr，ε_θθ，ε_zz为正应变，γ_rθ，γ_θz，γ_rz为主应变；u_r，u_θ，u_z分别为三个方向上的位移；r为径向半径；

通过上述方程的求解得到骨痂单元的正应变ε_rr，ε_θθ，ε_zz，由正应变可得到骨痂单元受到的畸变应变：

D=12(ϵrr-ϵθθ)2+(ϵrr-ϵzz)2+(ϵθθ-ϵzz)2---(17)]]>

式中，D为骨痂单元受到的畸变应变；ε_rr，ε_θθ，ε_zz分别为各个方向上的正应变；

骨痂单元中液体的流速V为：

V=-ku∂p∂r---(18)]]>

其中，k为骨痂中液体的达西渗透系数；u液体粘度；p为液体压力；

由此可得到骨痂单元所受的力学刺激S为：

S=Da+Vb---(19)]]>

其中，D为骨痂单元受到的畸变应变；V为骨痂单元中液体流速；a，b分别为经验常数；

骨痂单元组织分化模块(3)用于对单元内组织分化进行仿真，使单元内各组织含量得到更新，从而使单元材料属性得到更新，进而得到下一迭代步中所需要的单元力学刺激；

所述的骨痂单元组织分化模块(3)实现其功能的具体过程为：

1)设置输入变量隶属度函数

骨痂单元组织分化过程由六个输入变量协同决定，分别为单元力学刺激，骨痂单元骨含量，骨痂单元软骨含量，骨痂单元血供，周围骨痂单元骨含量，周围骨痂单元血供；输入变量的精确值通过隶属度函数转变为相应输入变量的隶属度；输入变量隶属度函数由梯形函数表示；

2)输入隶属度函数模糊化

将输入变量的隶属度函数用模糊逻辑语言值来表示：

单元力学刺激可分为三个等级，分别为：低，中，高；

单元血供可分为三个等级，分别为：低，中，高；

单元骨含量可分为三个等级，分别为：低，中，高；

单元软骨含量可分为三个等级，分别为：低，中，高；

3)设置输出变量隶属度函数

输出变量共有三个，分别为骨痂单元血供改变量，骨痂单元骨含量改变量和骨痂单元软骨含量改变量；设置输出变量隶属度函数，输出隶属度函数由梯形函数表示，且将输出隶属度函数用如下的模糊逻辑语言值表示：

骨痂单元血供改变量分为三个等级，分别为：低，中，高；

骨痂单元骨含量改变量分为三个等级，分别为：低，中，高；

骨痂单元软骨含量改变量分为三个等级，分别为：低，中，高；

4)模糊规则的编写

骨痂内组织分化主要包括：血管再生，膜内骨化，软骨生成，软骨钙化，软骨骨化五个过程；通过模糊规则分别表示这五个组织分化过程，由此可得到组织分化后，骨痂单元中骨的模糊逻辑语言值，软骨的模糊逻辑语言值以及血供的模糊逻辑语言值；具体过程如下：

过程一、血管再生

如果单元应力刺激为非高且骨痂单元血供为低且周围骨痂单元血供为非低，那么骨痂单元血供增加；

如果单元应力刺激为非高且骨痂单元血供为中且周围骨痂单元血供为高，那么骨痂单元血供增加；

如果单元应力刺激为非高且骨痂单元血供为高，那么骨痂单元血供增加；

过程二、膜内骨化

如果骨痂单元软骨含量为低且单元应力刺激为中且骨痂单元血供为高且周围骨痂单元血供为高，那么骨痂单元骨含量增加；

过程三、软骨生成

如果骨痂单元骨含量为非高且骨痂单元软骨含量为低且单元应力刺激为中，那么骨痂单元软骨含量增加；

如果骨痂单元软骨含量为非低且单元应力刺激为高，那么骨痂单元软骨含量增加；

如果骨痂单元软骨含量为高且单元应力刺激为中，那么骨痂单元软骨含量增加；

过程四、软骨钙化

如果骨痂单元软骨含量为非低且单元应力刺激为中且骨痂单元血供为非低且周围骨痂单元骨含量为非低，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

如果骨痂单元软骨含量为非低且单元应力刺激为高且骨痂单元血供为非低且周围骨痂单元骨含量为非低，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

过程五、软骨骨化

如果骨痂单元骨含量为高且骨痂单元软骨含量为低且单元应力刺激为中且骨痂单元血供为非低且周围骨痂单元骨含量为高，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

如果骨痂单元骨含量为高且骨痂单元软骨含量为低且单元应力刺激为低且骨痂单元血供为非低且周围骨痂单元骨含量为高，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

5)输出变量反模糊化

通过模糊规则的作用，得到输出变量的模糊逻辑语言值，利用面积中心法得到输出变量的隶属度；

6)更新输出变量含量

将得到的输出变量隶属度经过转化得到输出变量的改变量，从而得到经过组织分化后，骨痂单元骨组织，软骨组织，结缔组织含量以及骨痂单元血供；具体过程如下：

a.计算经组织分化后骨痂单元骨组织含量

μBone,n(k+1)=μBone,n(k)+ΔtddtμBone,n(k)---(20)]]>

式中，为第n个骨痂单元在第k+1迭代步中骨组织的含量；为第n个骨痂单元在第k迭代步中骨的含量；为单位时间第n个骨痂单元在第k迭代步中骨的生成量；△t为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中骨的改变量由下式得到：

ddtμBone,n(k)=wBone(k)TBone---(21)]]>

式中，为单位时间第n个骨痂单元在第k迭代步中骨的生成量；为骨痂单元中骨含量的隶属度；T_Bone在骨痂单元骨的转化率；

b.计算经组织分化后骨痂单元软骨骨组织含量

μCartilage,n(k+1)=μCartilage,n(k)+ΔtddtμCanilage,n(k)---(22)]]>

式中，为第n个骨痂单元在第k+1迭代步中软骨的含量；为第n个骨痂单元在第k迭代步中软骨的含量；为单位时间第n个骨痂单元在第k迭代步中软骨的生成量；△t为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中软骨的改变量由下式得到：

ddtμCartilage,n(k)=wCartilage(k)TCartilage---(23)]]>

式中，为单位时间第n个骨痂单元在第k迭代步中软骨的生成量；为骨痂单元中软骨含量的隶属度；T_Cartilage在骨痂单元软骨的转化率；

c.计算经组织分化后骨痂单元血供

μPerfusion,n(k+1)=μPerfusion,n(k)+ΔtddtμPerfusion,n(k)---(24)]]>

式中，为第n个骨痂单元在第k+1迭代步中结缔组织的含量；为第n个骨痂单元在第k迭代步中结缔组织的含量；为单位时间第n个骨痂单元在第k迭代步中结缔组织的生成量；△t为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中软骨的改变量由下式得到：

ddtμPerfusion,n(k)=wPerfusion(k)TPerfusion---(25)]]>

式中，为单位时间第n个骨痂单元在第k迭代步中结缔组织的生成量；为骨痂单元中结缔组织含量的隶属度；T_Perfusion在骨痂单元结缔组织的转化率；

d.计算经组织分化后骨痂单元结缔组织含量

骨痂单元由骨组织，软骨组织，结缔组织三部分组成，三者之间的关系如下所示：

μ_Bone+μ_Cartilage+μ_ConnTissue＝1 (26)

式中，μ_Bone为骨痂单元骨组织含量；μ_Cartilage为骨痂单元软骨组织含量；μ_ConnTissue为骨痂单元结缔组织含量；

由此可得组织分化后骨痂单元结缔组织含量；

程序终止判断模块(4)用于判断程序是否终止，若不满足终止条件，程序进行下一迭代步；若满足程序终止条件，则程序结束并输出愈合时间。

2.根据权利要求1所述的一种基于组织分化的长骨骨折愈合仿真系统，其特征在于：所述的程序终止判断模块(4)实现其功能的具体过程为：

1)判断骨痂单元材料属性

判断当前骨痂单元材料属性是否与骨的材料属性相同，若不同，则程序执行下面的2)步骤；若相等，程序执行3)步骤；

2)更新骨痂单元材料属性

若骨痂单元材料属性不等于骨的材料属性，则骨痂单元材料属性进行更新，进入下一个迭代步，骨痂材料属性更新公式如下：

ECallus,n(k+1)=EBoneμBone,n(k+1)3+ECartilageμCartilage,n(k+1)3+EConnTissueμConnTissue,n(k+1)3---(27)]]>

式中，为第n个骨痂单元在第k+1迭代步中的弹性模量；E_Bone为骨的弹性模量；为第n个骨痂单元在第k+1迭代步中骨的含量；E_Cartilage为软骨的弹性模量；为第n个骨痂单元在k+1迭代步中的软骨的含量；E_ConnTissue为结缔组织的弹性模量；为第n个骨痂单元在k+1迭代步中结缔组织的含量；

vCallus,n(k+1)=vBoneμBone,n(k+1)+vCartilageμCartilage,n(k+1)+vConnTissueμConnTissue,n(k+1)---(28)]]>

式中，为第n个骨痂单元在第k+1迭代步中的泊松比；ν_Bone为骨的泊松比；为第n个骨痂单元在第k+1迭代步中骨的含量；ν_Cartilage为软骨的泊松比；为第n个骨痂单元在k+1迭代步中的软骨的含量；ν_ConnTissue为结缔组织的泊松比；为第n个骨痂单元在k+1迭代步中结缔组织的含量；

3)程序结束

若所有骨痂单元材料属性等于骨的材料属性，程序终止并输出愈合时间。