一种基于微观结构参数个体化关节软骨仿真方法

摘要

一种基于微观结构参数个体化关节软骨仿真方法，本发明涉及基于微观结构参数个体化关节软骨仿真方法。本发明是为了解决现有获取专门患者关节软骨参数的手段为生理切片等破坏性手段的问题。本发明步骤为：步骤一：关节软骨微观结构参数赋值；步骤二：根据步骤一建立个体化线弹性关节软骨微观模型；步骤三：将粘弹性材料参数整合到步骤二建立的个体化线弹性关节软骨微观模型中，建立粘弹性关节软骨微观模型；步骤四：对步骤三建立的粘弹性关节软骨微观模型进行有限元仿真。本发明应用于生物医学工程领域。

说明书

技术领域

本发明涉及基于微观结构参数个体化关节软骨仿真方法。

背景技术

关节软骨是一种特殊的结缔组织，位于关节两端的关节面位置，并在其所处的生理部位为关节的活动提供低磨损和摩擦的光滑界面，起着缓冲震荡，传递载荷等不可替代的作用。鉴于以上功能，为适应多变的力学环境，关节软骨具有复杂的结构。经过多年国内外广大学者的不懈努力，人类对于关节软骨这一特殊组织的了解越来越深入，研究者通过各种实验手段去探究关节软骨的力学特性并尝试用各种仿真模拟方法把已发现的关节软骨特性表现出来。

但是目前的研究都是基于实验-宏观力学特性分析-宏观力学特性模拟的研究路线进行的，由于关节软骨力学性能测试实验必须利用关节软骨的生理切片在体外进行，通过这样的研究方法想要获得专门患者的关节软骨力学特性是不可能的，如何利用非破坏手段(无需生理切片)获取专门患者关节软骨参数，进而建立起专门患者关节软骨模型，这一问题的解决对于关节软骨的特性研究和临床诊断具有重要价值。

发明内容

本发明是为了解决现有获取专门患者关节软骨参数的手段为生理切片等破坏性手段和宏观模型不能表达专门患者个体化关节软骨特性的问题，而提出的一种基于微观结构参数个体化关节软骨仿真方法。

一种基于微观结构参数个体化关节软骨仿真方法按以下步骤实现：

步骤一：关节软骨微观结构参数赋值；

步骤二：根据步骤一建立个体化线弹性关节软骨微观模型；

步骤三：将粘弹性材料参数整合到步骤二建立的个体化线弹性关节软骨微观模型中，建立粘弹性关节软骨微观模型；

步骤四：对步骤三建立的粘弹性关节软骨微观模型进行有限元仿真。

发明效果：

(1)传统关节软骨力学特性模型建立主要是通过获取关节软骨的生理切片进行体外实验，依据实验数据反推出关节软骨应具有的力学特征，进而寻找能表征这样宏观特性的数学模型，依据这样的途径所建立的关节软骨模型为宏观模型，宏观模型显然不能实现应用个体化参数进行关节软骨建模。本发明依据复合材料微观力学思想，建立个体化关节软骨微观模型。

(2)传统的复合材料微观模型中很少考虑粘弹性的影响，微观力学中更没有能够兼顾关节软骨材料两种力学特性即与时间无关的弹性材料特性和与时间有关的粘弹性特性方面的模型。本发明通过结合宏观模型与微观模型的研究方法，将粘弹性整合进关节软骨力学微观模型中。

(3)传统的关节软骨模型仿真都是基于宏观模型设计的，没有考虑个体关节软骨模型仿真计算问题，本发明不仅提出了建立个体化弹性关节软骨模型，还在此模型基础上建立个体化粘弹性关节软骨模型，不仅引进个体微观结构参数同时兼顾全面的表现出关节软骨复杂力学性能，本发明有效解决了新模型的计算问题。

本发明通过建立个体化关节软骨模型，该模型可以用于医疗诊断和软骨退行性病变病情预测。关节软骨一旦发生损伤或病变，很难自我痊愈，1990年美国CDC疾病预防控制中心数据显示，关节软骨造成的经济负担仅次于高血压、心脏病和精神病，而我国每年用于关节软骨治疗的费用高达1500亿人民币，个体化关节软骨模型能帮助医生掌握软骨特性损伤程度，对于是否手术等进一步的治疗方案做出正确判断。

本发明通过建立个体化关节软骨模型，设计师可以根据特定病人数据而不是标准的解剖学几何数据来设计并制作种植体，这样就极大地减少了种植体设计的出错空间，同时，也能帮助预测假体骨与关节软骨性能的匹配程度，为假体骨的设计提供有价值的参考。

本发明通过建立个体化关节软骨模型，能够替代组织机械力学测试实验，可以实现通过非破坏性的或者微创的手段获取组织性能。

本发明通过建立个体化关节软骨模型，可以实现针对专门患者进行的外科手术仿真与规划。可以利用非破坏性手段实现对移植修复和组织工程化组织的评价。

附图说明

图1为特征体元图；

图2为全局坐标系与局部坐标系示意图；

图3为特征向量E₃在球面坐标中的取向示意图；

图4为基于微观结构参数个体化关节软骨建模与仿真分析实施方法流程图，图中的CTO为一致切线算子。

具体实施方式

具体实施方式一：一种基于微观结构参数个体化关节软骨仿真方法包括以下步骤：

步骤一：关节软骨微观结构参数赋值；

步骤二：根据步骤一建立个体化线弹性关节软骨微观模型；

步骤三：将粘弹性材料参数整合到步骤二建立的个体化线弹性关节软骨微观模型中，建立粘弹性关节软骨微观模型；

步骤四：对步骤三建立的粘弹性关节软骨微观模型进行有限元仿真。

随着医学图像获取方法的发展，结合连续介质宏观力学理论、复合材料微观力学理论、有限元仿真计算方法，本发明按照材料-参数-性能的全新研究路线建立了专门患者关节软骨的模型与仿真分析系统。

基于微观结构参数个体化关节软骨建模与仿真分析实施方法流程图如图4所示。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中关节软骨微观结构参数赋值具体为：

基于微观结构参数个体化关节软骨模型参数分为两类：一类是组分材料参数，关节软骨的组分分为纤维和基质两部分，与粘弹性纤维相关联的材料参数是λ_F,μ_F，与弹性基质相关联的材料参数为λ_M,μ_M；λ_F为纤维拉梅常数，μ_F为纤维剪切模量，λ_M为基质拉梅常数，μ_M为基质剪切模量；第二类参数是专门患者结构参数，包括特征体元的纤维体积分数V_F和特征体元在整个关节软骨中纤维取向分布函数Φ(θ,φ)；

第一类组分材料参数通过力学试验获取；第二类结构参数获取包括以下步骤：

步骤一一：关节软骨微观结构参数的选取；

选择纤维体积分数和纤维取向分布函数作为建立个体化关节软骨模型依据的微观结构参数，所述纤维体积分数为纤维所占面积与观察面积之比，纤维取向分布函数是纤维三维空间分布的一种方位表达形式；

步骤一二：关节软骨微观结构参数的计算；基于纤维占光学图像的比例确定纤维体积分数，分为两个步骤；

步骤一二一：应用全局算法和局部算法，选择阈值使光学成像分为纤维和背景两个区域；

步骤一二二：纤维占据图像区域的计算；

$>V_F=\frac{\underset{k=1}{\overset{k=N}{\Sigma }}{A}_{F_K}}{{\mathrm{NA}}_P}$>

式中：为第k层图像中纤维所占的面积；N为断层图像的总数；A_P为图像的总面积；

步骤一三：采用基于正交滤波器的方法，通过滤波器定义获取滤波输出及基于取向张量获得纤维取向分布，确定纤维取向分布函数；具体确定方法如下：方向函数为：

式中：ω为频率向量，符号“^”是常规化的相关矢量，是常规化的频率向量，为定义滤波器k方向的单位向量；

取向张量的建立需要滤波器的输出q_k，即图像与每个滤波器卷积后的结果；滤波器F_k定义在频域上，通过卷积定理可知滤波器的输出q_k为：

式中：为傅里叶逆变换；I′(ω)为与图像相关的像素强度函数；

R(ω)为滤波器的径向函数；

得到滤波器输出q_k后，每个像素处得取向张量T可由下式得到：

$>T=\underset{k}{\Sigma }\left|\right|q_k\left|\right|M_k$>

式中：I为单位矩阵；

从式中可以看出T是2阶3维的，它有3个特征值e₁≥e₂≥e₃和3个相应的特征向量E₁，E₂，E₃；若某点处e₁为最大特征值，则E₁为这一点上最大强度变化的方向，若该点为线性图像特征(如纤维)的一部分，则特征的长轴方向与特征向量E₃(最小特征值e₃对应的特征向量)近似。由此，通过对特征值大小的比较，可得到确定方向的方法：

a)e₁≈e₂>>e₃表示线性特征；

b)e₁>>e₂≈e₃表示平面特征；

c)e₁≈e₂≈e₃表示没有确定取向的各向同性区域。

上述考虑为减少来自于噪音等区域的的伪象影响提供了一条途径，最终由图3和可以得到在球面坐标系中描述特征向量E₃的表达式为：

式中：和分别为E₃在X，Y和Z上的分量；

由此可以从每个像素处估计的取向张量中转换出在球面坐标系中表示的纤维取向的直方图，这些直方图可以常规化为纤维取向分布函数。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二中建立个体化线弹性关节软骨微观模型的具体过程为：

步骤二一：基于微观结构参数线弹性关节软骨微观模型的建立；

关节软骨微观力学模型特征体元为同心圆柱模型，如图1所示，在特征体元上建立直角坐标系(1,2,3)，1沿着纤维轴向方向，称为特征体元的径向，2和3位于垂直于纤维的平面内，称为特征体元的横向；由1、2和3建立的坐标系称为复合材料的局部坐标系，全局坐标系放置于复合材料的中心位置；全局坐标系与局部坐标系如图2所示。

每一个特征体元由一对同轴的圆柱体组成，内部圆柱为内部纤维，外部圆柱为外部基质；

设内部纤维和外部基质是线弹性的，线弹性模型本构方程表述为：

σ＝λeI+2με＝Cε

式中，σ为应力，ε为应变，e为体积应变，λ,μ是弹性模量E和泊松比v的函数；

$>\lambda =\frac{Ev}{(1+v)(1-2v)}$>

$>\mu =\frac{E}{2(1+v)}$>

纤维的线弹性特性由λ_F,μ_F确定，基质的线弹性特性由λ_M,μ_M确定；

步骤二二：计算刚度矩阵；

对于每一个特征体元，横观各向同性刚度矩阵为：

$>C_{ijkl}=\left(\begin{array}{cccccc}{C}_{\left[11\right]}& C_{\left[12\right]}& C_{\left[12\right]}& 0& 0& 0\\ & C_{\left[22\right]}& C_{\left[23\right]}& 0& 0& 0\\ & & C_{\left[22\right]}& 0& 0& 0\\ & & & \frac{C_{\left[22\right]}-C_{\left[23\right]}}{2}& 0& 0\\ & Sym& & & C_{\left[55\right]}& 0\\ & & & & & C_{\left[55\right]}\end{array}\right)$>

式中，Sym代表对称矩阵中的对称项；

将上式中的各系数替换为工程常数的表达式如下：

$>C_{ijkl}=\left(\begin{array}{cccccc}{E}_{11}+4v_{12}^2{K}_{23}& 2K_{23}{v}_{12}& 2K_{23}{v}_{12}& & & \\ & {\mu }_{23}+K_{23}& K_{23}-{\mu }_{23}& & & \\ & & {\mu }_{23}+K_{23}& & & \\ & & & {\mu }_{23}& & \\ & Sym& & & {\mu }_{12}& \\ & & & & & {\mu }_{12}\end{array}\right)$>

式中：E₁₁为纵向的杨氏模量；v₁₂(＝v₁₃)为纵向泊松比；K₂₃为平面应变体积模量；μ₁₂(＝μ₁₃)为面内剪切模量；μ₂₃为特征体元的横观剪切模量；

表征的工程常数与各组分(纤维和基质)材料参数之间的关系式，其表达形式如下：

$>E_{11}=\frac{V_F{\mu }_F(3{\lambda }_F+2{\mu }_F)}{{\xi }_1}+\frac{(1-V_F){\mu }_M(3{\lambda }_M+2{\mu }_M)}{{\xi }_2}+\frac{(1-V_F)V_F{\mu }_M{\xi }_3^2}{{\xi }_4}$>

$>v_{12}=V_M{\mu }_M+V_F{\mu }_F+\frac{V_F{V}_M(v_F-v_M)(\frac{3}{3K_M+{\mu }_M}-\frac{3}{3K_F+{\mu }_F})}{\frac{3V_M}{3K_F+{\mu }_F}+\frac{3V_F}{3K_M+{\mu }_M}+\frac{1}{{\mu }_M}}$>

$>K_{23}={\xi }_2+\frac{V_F}{\frac{1}{{\xi }_1-{\xi }_2}+\frac{1-V_F}{{\lambda }_M+2{\mu }_M}}$>

$>{\mu }_{12}=\frac{{\mu }_M[(1+V_F){\mu }_F+(1-V_F){\mu }_M]}{(1-V_F){\mu }_F+(1+V_F){\mu }_M}$>

$>{\mu }_{23}={\mu }_M+\frac{V_F{\mu }_M}{\frac{{\mu }_M}{{\mu }_F-{\mu }_M}+\frac{(1-V_F)({\lambda }_M+3{\mu }_M)}{2({\lambda }_M+2{\mu }_M)}}$>

式中：

ξ₁＝λ_F+μ_F

ξ₂＝λ_M+μ_M

$>{\xi }_3=\frac{{\lambda }_F}{{\xi }_1}-\frac{{\lambda }_M}{{\xi }_2}$>

$>{\xi }_4=1+\frac{V_F{\mu }_F}{{\xi }_2}+\frac{(1-V_F){\mu }_M}{{\xi }_1}$>

V_α，E_α，v_α，μ_α和K_α分别为α相(α纤维或基质)的体积分数，弹性模量，泊松比，剪切模量和体积模量；

则刚度矩阵中的各刚度系数表达式为：

$>C_{\left[11\right]}=({\lambda }_M+2{\mu }_M)(1-V_F)+V_F[{\lambda }_F+2{\mu }_F-\frac{(1-V_F){({\lambda }_F-{\lambda }_M)}^2}{(1-V_F){\xi }_1+V_F{\xi }_2+{\mu }_M}]$>

$>C_{\left[12\right]}=\frac{(1-V_F){\lambda }_M({\mu }_F+{\mu }_M)+{\lambda }_F({\lambda }_M+2V_F{\mu }_M)}{(1-V_F){\xi }_1+V_F{\xi }_2+{\mu }_M}$>

C_[22]＝λ_M+2μ_M+ξ₅+ξ₆

C_[23]＝λ_M+ξ₅-ξ₆

$>C_{\left[55\right]}={\mu }_{12}=\frac{{\mu }_M[(1+V_F){\mu }_F+(1-V_F){\mu }_M]}{(1-V_F){\mu }_F+(1+V_F){\mu }_M}$>

式中：

$>{\xi }_5=\frac{V_F}{\frac{1}{{\xi }_1-{\xi }_2}+\frac{1-V_F}{{\lambda }_M+2{\mu }_M}}$>

$>{\xi }_6=\frac{V_F{\mu }_M}{\frac{{\mu }_M}{{\mu }_F-{\mu }_M}+\frac{(1-V_F)({\lambda }_M+3{\mu }_M)}{2({\lambda }_M+2{\mu }_M)}}$>

通过上述推导，可以得到局部坐标系(1，2，3)中的刚度矩阵和相应的本构方程；通过适当的变换，本构方程可以在全局坐标系表达为下式的：

σ′_ij＝C′_ijklε'_kl

式中：“′”表示参数在全局坐标系中的表达形式；

C｀_ijkl张量能够在全局坐标系中通过长轴方向的θ和φ角度来描述单个特征体元的刚度；由于假设整个组织是由多个特征体元以随机分布的形式组成的，组织的整体刚度可以通过对θ和φ所有角度范围内的C｀_ijkl积分得到。在真实的组织样本中，纤维的取向在各个方向上是不均匀的，为了考虑这种随机的各向异性就必须引入合适的加权函数。假设样本有足够大的体积，并含有足够多的纤维，则这个加权函数可以表示为组织单位体积内组织纤维取向的统计分布密度函数。使纤维取向分布函数Φ(θ，φ)为这样一个统计分布，组织作为整体的等效刚度矩阵可由下式给出：

$>{\bar{C}}_{ijkl}=\underset{0}{\overset{\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\Phi (\theta ,\phi )C_{ijkl}^{\prime }sin\phi d\phi d\theta$>

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤三中建立粘弹性关节软骨微观模型的具体过程为：

步骤三一：确定基于个体化微观结构参数粘弹性关节软骨粘弹性参数；

将粘弹性材料参数整合进个体化线弹性关节软骨微观模型中，采用Prony级数形式替换原模型(个体化线弹性关节软骨微观模型)中的弹性参数；

在个体化线弹性关节软骨微观模型中，纤维相的参数是λ_F,μ_F，由于胶原纤维和蛋白多糖基质均具有粘弹性特性；胶原纤维粘弹性远高于基质粘弹性；纤维相只有内部出现剪切应力时表现粘弹性。因此，利用粘弹性Prony级数参数替换μ_F；粘弹性Prony级数表示为：

$>\Omega \left(t\right)={\Omega }^{\infty }+\underset{K=1}{\overset{N}{\Sigma }}{\Omega }^K{e}^{-t/{\tau }^K}$>

其中Ω^∞为平衡模量，Ω^K是松弛模量，τ^K是松弛时间常数，t为时间；

步骤三二：确定应力应变关系矩阵；

$>{\sigma }_{ij}\left(t\right)={\bar{C}}_{ijkl}^0{ϵ}_{kl}-\underset{K=1}{\overset{N}{\Sigma }}{D}_{ij}^K$>

式中，将C_ijklε_kl划分为时间相关项和时间无关项两个部分，为时间无关项的系数，为时间相关项；

时间相关相定义：

$>D_{ij}^K=\underset{0}{\overset{t}{\int }}{\bar{C}}_{ijkl}^K(1-e^{(t^{\prime }-t)/x^K})\times \frac{{\mathrm{dϵ}}_{kl}}{{\mathrm{dt}}^{\prime }}{\mathrm{dt}}^{\prime }$>

式中，为全局等效刚度矩阵中与时间相关的项。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤四中对步骤三建立的粘弹性关节软骨微观模型进行有限元仿真的输入变量、输出变量和计算公式为：

输入变量：

材料参数：λ_F、μ_F、λ_M和μ_M，松弛时间常数τ^K；

结构参数：纤维体积分数V_F和纤维取向分布函数Φ(θ,φ)；

输出变量：用一致切线算子表示关节软骨力学特性，一致切线算子为

计算公式：

计算全局等效刚度矩阵：

$>{\bar{C}}_{ijkl}=\underset{0}{\overset{X}{\int }}\underset{0}{\overset{X}{\int }}\Phi (\theta ,\phi )C_{ijkl}^{\prime }\sin \phi d\phi d\theta$>

更新时间相关项：

$>D_{{\mathrm{ij}}_{n+1}}^K=\frac{\frac{D_{{\mathrm{ij}}_n}^K}{\Delta t}+\frac{{\bar{C}}_{ijkl}^K{ϵ}_{{\mathrm{kl}}_{n+1}}}{{\tau }^K}}{\frac{1}{\Delta t}+\frac{1}{{\tau }^K}}$>

更新包含时间相关项和时间无关项的应力应变等式：

$>{\sigma }_{{\mathrm{ij}}_{n+1}}\left(t\right)={\bar{C}}_{ijkl}^0{ϵ}_{{\mathrm{kl}}_{n+1}}-\underset{K=1}{\overset{N}{\Sigma }}{D}_{{\mathrm{ij}}_{n+1}}^K$>

更新一致切线算子：

$>\frac{\partial {\sigma }_{{\mathrm{ij}}_{n+1}}}{\partial {\mathrm{\Delta ϵ}}_{kl}}={\bar{C}}_{ijkl}^0-\underset{K=1}{\overset{N}{\Sigma }}\frac{\Delta t}{\Delta t+{\tau }^K}{\bar{C}}_{ijkl}^K$>

其它步骤及参数与具体实施方式一至四之一相同。