CN103310072A - 基于力反馈的股骨生物力学有限元分析系统 - Google Patents

Abstract

基于力反馈的股骨生物力学有限元分析系统，涉及股骨建模软件开发领域。本发明为了解决现有的股骨建模分析软件因操作极其复杂，用其指导外部手术的成功率较低等问题。股骨CT图像灰度值提取模块用于导入多种格式的CT图像；股骨弹性模量赋值模块用于获取股骨体网格单元内心坐标；股骨生物力学有限元分析模块用于通过每个股骨体网格单元的弹性模量计算出每个股骨体网格单元刚度矩阵；力反馈模块用于根据上述应变或应力数据模拟出测试点的力学特性：力的大小和力的方向，力反馈模块还用于将力的大小和力的方向传递给外部的力感知设备的操纵手柄上，经数据转换后的力反馈器可以输出股骨实际变形时所产生的力，能真实地感知股骨模型的生物力学特性。

Description

基于力反馈的股骨生物力学有限元分析系统

技术领域

本发明涉及一种基于力反馈的股骨生物力学有限元分析系统，涉及股骨建模软件开发领域。

背景技术

近年来，国内外的虚拟现实技术在医学领域中得到了充分的发展。例如日本东京工业大学开发了一套新型微创虚拟手术系统，医生可以基于自主研发的该软件，结合虚拟现实眼镜、虚拟力反馈系统等外部硬件设备进行模拟的手术训练。CHAI3D开源库由美国开发，是针对虚拟手术中的力反馈关键技术所研发的一套触觉库，该触觉库已经被世界100多家顶级医院所应用。在澳大利亚也开发了一套被称为“触觉手术台”的虚拟手术系统，为实习医生开发模拟训练的环境。在国内也有香港中文大学、清华大学、国防科技大学、浙江大学等许多高校单位逐步进行研究并取得了一些成果，但与国外相比还是相对较为落后。

当前，在股骨建模的过程中遇到的主要问题是股骨生物力学特性的分析上，现有的股骨建模分析软件因操作极其复杂，要将研究成果应用到临床成本很高，用其指导外部手术的成功率较低，不能普及应用。

发明内容

本发明的目的是提供一种基于力反馈的股骨生物力学有限元分析系统，以解决现有的股骨建模分析软件因操作极其复杂，用其指导外部手术的成功率较低，不能普及应用等问题。

本发明为解决上述技术问题采取的技术方案是：

一种基于力反馈的股骨生物力学有限元分析系统，所述系统包括：

股骨CT图像灰度值提取模块股骨弹性模量赋值模块、股骨生物力学有限元分析模块和力反馈模块；

股骨CT图像灰度值提取模块，用于导入多种格式的CT图像，并把该图像存成以灰度值为元素的空间三维数据矩阵；股骨CT图像灰度值提取模块还用于将获得的CT灰度值信息进行预处理：利用灰度门限法筛选出股骨，而且能手动擦除噪点和边缘上的黏连组织的数据，存留以股骨作为研究对象的三维数据矩阵；

股骨弹性模量赋值模块用于获取股骨体网格单元内心坐标；

股骨弹性模量赋值模块还用于根据股骨CT图像灰度值提取模块得到的股骨三维数据矩阵和股骨体网格单元内心坐标计算出每个股骨体网格单元的内心灰度值；

股骨弹性模量赋值模块还用于根据股骨体网格单元的内心灰度值计算得到股骨体网格单元的表观密度，进一步计算出每个股骨体网格单元的弹性模量，建立股骨有限元模型；

股骨生物力学有限元分析模块用于通过每个股骨体网格单元的弹性模量计算出每个股骨体网格单元刚度矩阵，进一步地，通过给出的股骨外加载荷和固定约束得到股骨模型的整体刚度矩阵，进而求出股骨模型的每个股骨体网格单元的位移矩阵；最后通过股骨模型的每个股骨体网格单元的位移矩阵计算出股骨模型上所要测试点的应变或应力数据；

力反馈模块用于根据上述应变或应力数据模拟出测试点的力学特性：力的大小和力的方向，力反馈模块还用于将力的大小和力的方向传递给外部的力感知设备的操纵手柄上，使医生能真实地感知股骨模型的生物力学特性。

本发明的有益效果是：

在MATLAB环境下开发有限元分析程序对股骨生物力学建模的研究与今后的有限元软件的开发提出新的思路和方法，也可以在虚拟手术系统建模中的生物力学建模提供理论和方法上的基础。

做出来的软件与力反馈设备结合，将股骨生物力学特性反馈到医生的手中，增加手术过程中沉浸感、身临其境的感觉。本发明可以应用到骨科虚拟手术系统中，指导医生做手术，或者可以让医生进行罕见类型的手术规划，从而提高实际手术的成功率，同时可有效地解决在传统手术训练或规划时，因为尸体的数量少、成本高而无法让医生进行反复的手术训练的问题。

本发明依靠开放性软件MATLAB自主研发，降低了软件的成本，设计出来的程序操作起来较简单，运行所需占的内存较小，而且MATLAB软件平台是完全免费开放。结合上述的软件与力反馈器操纵杆结合起来，把有限元程序分析出来的应变数据传送感到力反馈器的接口上，经数据转换后的力反馈器可以输出股骨实际变形时所产生的力，这样就可以亲手感觉到股骨的生物力学特性。

本发明可通过MATLAB平台来实现编程，设计出了可以对股骨生物力学有限元分的程序。因为MATLAB平台是完全免费开放的。与ANSYS、ABAQUS等目前市面上比较流行的有限元分析软件相比较，具有软件的设计成本比较低，程序操作起来较简单，运行所需占的内存较小等特点。

附图说明

图1是本发明的流程框图(其中虚线框1是股骨CT图像灰度值提取模块，虚线框2是股骨弹性模量赋值模块，虚线框3是股骨生物力学有限元分析模块，虚线框4是力反馈模块)，图2是股骨三维重建模型图(三维几何模型，图中的箭头表示同一股骨模型的不同方位立体图)，图3是股骨体网格有限元模型(图中的箭头表示同一股骨模型的不同方位立体图)，图4是CT图像的阈值筛选前后对比图(其中图4a是原始CT图像，图4b是经过阈值筛选后的图像)，图5是对选取一侧股骨和手动擦除后的图像(其中图5a是选取一侧股骨后的图像，图5b是经过手动擦除后的图像)图5是对股骨事假载荷及约束的示意图。图6是股骨施加载荷及约束示意图。图7为两种股骨弹性模量与表观密度关系模型的比较图。

具体实施方式

具体实施方式一：如图1所示，本实施方式所述的一种基于力反馈的股骨生物力学有限元分析系统包括：

股骨CT图像灰度值提取模块1、股骨弹性模量赋值模块2、股骨生物力学有限元分析模块3、力反馈模块4；

股骨CT图像灰度值提取模块1，用于导入多种格式的CT图像，并把该图像存成以灰度值为元素的空间三维数据矩阵；股骨CT图像灰度值提取模块1还用于将获得的CT灰度值信息进行预处理：利用灰度门限法筛选出股骨，而且能手动擦除噪点和边缘上的黏连组织的数据，存留以股骨作为研究对象的三维数据矩阵；

股骨弹性模量赋值模块2用于获取股骨体网格单元内心坐标；

股骨弹性模量赋值模块2还用于根据股骨CT图像灰度值提取模块1得到的股骨三维数据矩阵和股骨体网格单元(四面体)内心坐标计算出每个股骨体网格单元的内心灰度值；

股骨弹性模量赋值模块2还用于根据股骨体网格单元的内心灰度值计算得到股骨体网格单元的表观密度，进一步计算出每个股骨体网格单元的弹性模量，也就是建立股骨有限元模型；

股骨生物力学有限元分析模块3用于通过每个股骨体网格单元的弹性模量计算出每个股骨体网格单元刚度矩阵，进一步地，通过给出的股骨外加载荷和固定约束得到股骨模型的整体刚度矩阵，进而求出股骨模型的每个股骨体网格单元的位移矩阵；最后通过股骨模型的每个股骨体网格单元的位移矩阵计算出股骨模型上所要测试点(可能会包括很多个股骨体网格单元)的应变或应力数据；

力反馈模块4用于根据上述应变或应力数据模拟出测试点的力学特性：力的大小和力的方向，力反馈模块4还用于将力的大小和力的方向传递给外部的力感知设备的操纵手柄上，使医生能真实地感知股骨模型的生物力学特性。

具体实施方式二：如图1～7所示，本实施方式中，所述股骨弹性模量赋值模块2实现其功能的具体过程为：

1)计算四面体内心

由有限元网格模型可以得到节点坐标，通过四面体的四个节点坐标计算内心坐标；设四面体A₁A₂A₃A₄的顶点A_i所对的侧面面积为S_i(i＝1，2，3，4)，顶点A_i的坐标为(x_i，y_i，z_i)(i＝1，2，3，4)，四面体内心的坐标为(x₁，y₁，z₁)，则

$\left.\begin{array}{c}{x}_1=\frac{s_1{x}_1+s_2{x}_2+s_3{x}_3+s_4{x}_4}{s_1+s_2+s_3+s_4}\\ y_1=\frac{s_1{y}_1+s_2{y}_2+s_3{y}_3+s_4{y}_4}{s_1+s_2+s_3+s_4}\\ z_1=\frac{s_1{z}_1+s_2{z}_2+s_3{z}_3+s_4{z}_4}{s_1+s_2+s_3+s_4}\end{array}\right\}---\left(1\right)$

其中节点坐标已知，侧面面积应用海伦定理求得；

2)计算内心灰度值

利用式(2)将股骨有限元网格模型与空间三维数据矩阵进行坐标转换(空间灰度值数据矩阵是以像素和层数为坐标轴存储的，而有限元网格数据是空间的直角坐标系，二者虽然是一一对应的，但计算时需要进行坐标转换)：

$\left.\begin{array}{c}{x}^{\′}=x+a/p\\ y^{\′}=y+b/p\\ z^{\′}=z+c\end{array}\right\}---\left(2\right)$

式中，a、b、c为有限元网格模型的坐标轴与空间三维数据矩阵的之间的x、y、z轴上的差值，p表示像素尺寸，p由CT图像信息可知；

根据式(2)，构建起节点直角坐标和空间像素点阵之间的关系；由此可以确定每个单元内心在像素矩阵中的位置，在一个四面体网格单元当中，单元同时在空间像素矩阵中每一单元内心均由8个像素点包围，而且8个像素点位置规则；通过空间8节点正六面体的内插插值公式求得单元内心处的灰度值；

$G_i(x_i,y_i,z_i)=\underset{k=1}{\overset{8}{\mathrm{\Σ}}}{N}_j(x,y,z)g(i,j,k)---\left(3\right)$

式中，G_i为所求的内心灰度值，N_j为8节点正六面体单元的形函数，g(i，j，k)为8个像素点处的灰度值；

3)计算表观密度

表观密度：单位体积的质量，指不包括髓脂时骨的质量除以试件体积得到的参变量，用ρ表示；

表观密度与CT值存在线性关系，用灰度值代替CT值计算表观密度；关系式如下：

CT＝G-1024          (4)

$\mathrm{\ρ}=\frac{{\mathrm{\ρ}}_1-{\mathrm{\ρ}}_0}{G_1-G_0}\·G---\left(5\right)$

ρ₁＝2，G₁＝2663；ρ₀＝0，G₀＝0          (6)

式中ρ₁，ρ₀为股骨最大最小表观密度，G₁，G₀为股骨最大最小灰度值；类似于温度的计量单位，为了使用的方便，人们将水的灰度值值定义为0，按照这个标度，骨松质的G为100-300，骨皮质的G约为2000；

计算完内心灰度值后，通过式(5)和(6)计算单元表观密度；

4)计算网格单元的弹性模量

股骨在建立各向异性力学模型时被看做是弹性体；针对股骨有限元单元，应用弹性力学基本方程建立单元应力-应变关系式：

$\left[\begin{array}{c}{\mathrm{\σ}}_X\\ {\mathrm{\σ}}_Y\\ {\mathrm{\σ}}_Z\\ {\mathrm{\τ}}_{\mathrm{XY}}\\ {\mathrm{\τ}}_{\mathrm{YZ}}\\ {\mathrm{\τ}}_{\mathrm{XZ}}\end{array}\right]=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& C_{14}& C_{15}& C_{16}\\ C_{21}& C_{22}& C_{23}& C_{24}& C_{25}& C_{26}\\ C_{31}& C_{32}& C_{33}& C_{34}& C_{35}& C_{36}\\ C_{41}& C_{42}& C_{43}& C_{44}& C_{45}& C_{46}\\ C_{51}& C_{52}& C_{53}& C_{54}& C_{55}& C_{56}\\ C_{61}& C_{62}& C_{63}& C_{64}& C_{65}& C_{66}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ϵ}}_X\\ {\mathrm{\ϵ}}_Y\\ {\mathrm{\ϵ}}_Z\\ {\mathrm{\γ}}_{\mathrm{XY}}\\ {\mathrm{\γ}}_{\mathrm{YZ}}\\ {\mathrm{\γ}}_{\mathrm{XZ}}\end{array}\right],$ 即{σ}=[C]·{ε}   (7)

式中：[C]为三维各向异性刚度矩阵；C_ij为独立的刚性系数。

根据格林公式及广义胡克定律，有下式成立：

$\frac{{\∂U}_0}{\∂{\mathrm{\ϵ}}_X}={\mathrm{\σ}}_X=C_{11}{\mathrm{\ϵ}}_X+C_{12}{\mathrm{\ϵ}}_Y+C_{13}{\mathrm{\ϵ}}_Z+C_{14}{\mathrm{\γ}}_{\mathrm{XY}}+C_{15}{\mathrm{\γ}}_{\mathrm{YZ}}+C_{11}{\mathrm{\γ}}_{\mathrm{XZ}}---\left(8\right)$

$\frac{{\∂U}_0}{\∂{\mathrm{\γ}}_{\mathrm{XY}}}={\mathrm{\τ}}_{\mathrm{XY}}=C_{41}{\mathrm{\ϵ}}_X+C_{42}+{\mathrm{\ϵ}}_Y+C_{43}{\mathrm{\ϵ}}_Z+C_{44}{\mathrm{\γ}}_{\mathrm{XY}}+C_{45}{\mathrm{\γ}}_{\mathrm{YZ}}+C_{46}{\mathrm{\γ}}_{\mathrm{XZ}}---\left(9\right)$

分别对(8)式中的切应变γ_XY和(9)式中的正应变ε_X求偏导，则有等式：

$\frac{{\∂}^2{U}_0}{\∂{\mathrm{\ϵ}}_X\∂{\mathrm{\γ}}_{\mathrm{XY}}}=C_{14},\frac{{\∂}^2{U}_0}{\∂{\mathrm{\γ}}_{\mathrm{XY}}\∂{\mathrm{\ϵ}}_X}=C_{41}---\left(10\right)$

由上式可知C₁₄＝C₄₁；即通过对能量守恒定律与应变位能的考察，可以证明C_ij＝C_ji。由此可以得到[C]为对称矩阵的结论，矩阵中的36个独立系数减少到21个。又因为骨骼是正交各向异性材料，由此刚度系数矩阵中的力学参数被简化为9个。

股骨在建立各向异性力学模型时被看做是弹性体；针对股骨有限元单元，应用弹性力学基本方程，并进行一系列简化后，建立单元应力-应变关系式：

$\left[C\right]=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& & & \\ C_{12}& C_{22}& C_{23}& & 0& \\ C_{13}& C_{23}& C_{33}& & & \\ & & & C_{44}& & \\ & 0& & & C_{55}& \\ & & & & & C_{66}\end{array}\right]---\left(11\right)$

在本发明中采用的求解方向是已知载荷情况求解位移情况，因此用应力分量表达应变分量的形式为：{ε}＝[S]·{σ}；其中[S]为柔度矩阵，由刚度矩阵[C]求逆变换得到；

$\left[S\right]=\left[\begin{array}{cccccc}{S}_{11}& S_{12}& S_{13}& & & \\ S_{12}& S_{22}& S_{23}& & 0& \\ S_{13}& S_{23}& S_{33}& & & \\ & & & S_{44}& & \\ & 0& & & S_{55}& \\ & & & & & S_{66}\end{array}\right]---\left(12\right)$

长期的临床实践证明，股骨沿骨干轴向的材料特性与骨干径向和环向的材料特性相差较大，而径向和环向则很接近，因此本发明将股骨材料看做是横观各向同性体，即柔性矩阵[S]中的力学参数只有5个，如下式：

$\left[S\right]=\left[\begin{array}{cccccc}\frac{1}{E_{\mathrm{XY}}}& -\frac{{\mathrm{\μ}}_{\mathrm{XY}}}{E_{\mathrm{XY}}}& -\frac{{\mathrm{\μ}}_{\mathrm{XY}}}{E_Z}& & & \\ -\frac{{\mathrm{\μ}}_{\mathrm{XY}}}{E_{\mathrm{XY}}}& \frac{1}{E_{\mathrm{XY}}}& -\frac{{\mathrm{\μ}}_{\mathrm{XY}}}{E_Z}& & 0& \\ -\frac{{\mathrm{\μ}}_{\mathrm{XZ}}}{E_{\mathrm{XY}}}& -\frac{{\mathrm{\μ}}_{\mathrm{XZ}}}{E_{\mathrm{XY}}}& \frac{1}{E_Z}& & & \\ & & & \frac{1}{G_{\mathrm{XY}}}& & \\ & 0& & & \frac{1}{G_{\mathrm{XZ}}}& \\ & & & & & \frac{1}{G_{\mathrm{XZ}}}\end{array}\right]---\left(13\right)$

其中，5个独立的力学参数是横向和轴向的弹性模量E，剪切模量G、横向和轴向的泊松比μ；这5个参数之间的转化关系如下：

$G_{\mathrm{XY}}=\frac{E_{\mathrm{XY}}}{2(1+{\mathrm{\μ}}_{\mathrm{XY}})}---\left(14\right)$

$\frac{{\mathrm{\μ}}_{\mathrm{ZX}}}{E_Z}=\frac{{\mathrm{\μ}}_{\mathrm{XZ}}}{E_{\mathrm{XY}}}---\left(15\right)$

通过上式的换算，5个独立的力学参数减少到3个；在实际计算中，剪切模量和泊松比常取定值，而弹性模量的计算模型有很多种。

5)股骨有限元模型

本发明中建立了两种有限元模型，分别是股骨弹性模量与表观密度关系的胞元理论模型(以下简称为理论模型)和股骨弹性模量与表观密度关系的力学实验模型(以下简称为实验模型)；两种模型在一定范围内(表观密度在0.3g/cm³～0.7g/cm³之间)弹性模量是略显差异，但不是很明显。因此可以认为在相同条件下，理论模型和实验模型之间是可以相互代替，如图。

a.股骨弹性模量与表观密度关系的胞元理论模型：

基于吉林大学力学系的朱兴华等提出的股骨弹性模量与表观密度关系的胞元理论模型，得到弹性模量与表观密度的分段函数式：

$\left\{\begin{array}{c}1007{\mathrm{\&rho;}}^2,\mathrm{\&rho;}\&le;0.25\\ 255\mathrm{\&rho;},0.25<\mathrm{\&rho;}\&le;0.4\\ 2972{\mathrm{\&rho;}}^2-933\mathrm{\&rho;},0.4<\mathrm{\&rho;}\&le;1.2\\ 1736{\mathrm{\&rho;}}^{3.25},\mathrm{\&rho;}>1.2\end{array}\right.---\left(16\right)$

式中，弹性模量的单位是MPa，表观密度的单位是g/cm3。

b.股骨弹性模量与表观密度关系的力学实验模型：

基于美国斯坦福大学的Beaupre等和德国于利希研究中心的Fischer等研究，将骨骼表观密度与弹性模量的关系取为如下分段函数：

$E=\left\{\begin{array}{c}2014{\mathrm{\&rho;}}^{2.5},\mathrm{\&rho;}\&le;1.2\\ 1736{\mathrm{\&rho;}}^{3.5},\mathrm{\&rho;}>1.2\end{array}\right.---\left(17\right)$

其它组成及连接关系与具体实施方式一相同。

具体实施方式三：如图1～7所示，本实施方式中，所述股骨生物力学有限元分析模块3实现其功能的具体过程为：

有限元方法是近似求解一般连续域问题的数值方法，利用有限元方法解决问题得到的是数值解。其求解的过程包括离散、集合、建立微分方程、线性代数方程及求解，列方程时常选择位移为未知量，而求解时常应用最小势能原理和虚功原理。具体步骤如下：

1)结构离散化：

有限元分析的第一步就是结构离散化。在处理不规则表面的组织结构时，将其离散成有限个单元体，利用有限个单元体代替原来的组织结构进行计算。通过研究单元体与单元体之间的结构关系来研究整体的结构关系，包括节点位移、单元应变等，这在一定程度上也算是简化模型。

2)选择形函数：

有限元方法就是要建方程解方程，方程中的未知量就是位移和力，而目前大多数研究位移关系。选择合适的位移关系可以轻松的解开方程，本文选择线性的位移函数：

{f}＝[N]{δ}^ε          (18)

式中：{f}表示单元内任意一点的位移列阵；{δ}^ε表示节点位移列阵；[N]表示形状函数矩阵。

3)单元力学特性：

离散后的组织结构通过单元相连接，单元的力学特性能够体现整体的力学特性。单元刚度矩阵就是表现力学特性的一个矩阵：矩阵中的每个元素都是一个刚度系数，表示位移的变化引起的力的变化。其计算方法如下：

首先应用几何方程建立位移和应变的关系：

{ε}＝[B]{δ}^ε          (19)

式中，{ε}是单元内任意一点的应变列阵，[B]是系数矩阵。

然后将上式代入广义胡克定律关系式(物理方程)得到：

{σ}＝[D][B]{δ}^ε          (20)

式中，{σ}是单元内任意一点的应力列阵，[D]是与材料有关的弹性矩阵。

再由平衡微分方程得到力和位移的关系：

{R}^ε＝[K]^e{δ}^ε          (21)

式中，[K]^e就是单元刚度矩阵，它的计算公式如下：

[K]^e＝∫∫∫[B]^T[D][B]dxdydz          (22)

4)整体平衡方程：

研究完单元的力学特性，我们要通过单元的特性研究整体的，这就需要将离散的单元整合到一起，组装整体刚度矩阵。单元内的位移和力都是通过节点相互传递的，如果外载荷或约束的作用点不在节点上，就需要通过等效的方法将作用点转移到节点上来，单元内部的传递有时也需要等效转移，这在整合单元是很常用的方法。在整体上建立组织结构的平衡方程为：

[K]{δ}＝{R}          (23)

式中，[K]为整体刚度矩阵，{δ}为整体节点位移，{R}为整体节点载荷。

5)解方程：解方程后得出应力与应变的数据；

通过上面几步建立了有限元微分方程，其中整体刚度矩阵是奇异对称矩阵，这在一定程度上减少了计算量。同时利用混合边界条件将外载列阵和约束矩阵引入辅助求解，使得微分方程转化成线性代数方程，降低了方程求解的难度。

其它组成及连接关系与具体实施方式一或二相同。

针对本发明再进行如下阐述：

作为有限元分析软件，实现股骨CT图像灰度值提取、股骨弹性模量赋值、股骨生物力学模型有限元分析等功能。

1CT图像灰度值提取以及图像的预处理

1)MATLAB中可以导入多种格式的CT图像，并把该图像存成以灰度值为元素的空间三维点阵。

2)在上步中获得的CT灰度值信息中进行预处理：利用灰度门限法筛选出股骨，而且可以手动擦除噪点和边缘上的黏连组织。

2网格单元弹性模量的赋值

MATLAB可以通过查找节点坐标精确的对每个网格单元进行弹性模量赋值，

其中，运用的算法有：计算四面体内心算法；提取灰度值算法；内心灰度值计算算法；弹性模量计算算法。最后把该模型以数据矩阵的形式保存。

3股骨生物力学有限元分析模块

1)分析过程中所需要的节点坐标、单元编号、单元密度值等保存成数据矩阵，分别命名为：fna_nodes.mat、fna_elements.mat、bgmd.mat，需要时函数调用。与上步的CT图像的灰度值信息结合，就可以实现灰度值与表观密度之间的转换。

2)在MATLAB有限元分析程序中包含三个计算公式，分别在编译调试成功后保存成M文件：T3D4_Stiffness.m(单元刚度矩阵)、T3D4_Assembly.m(组装整体刚度矩阵)、T3D4_Strain.m(单元应变)，需要时直接调用运行；

3)手动输入载荷和约束条件，构建外载列阵；

4)矩阵运算采用高斯消元法，这样有利于微分方程求解和兼顾运算速度；

5)为了清理内存减少数据垃圾，计算过程中的过程量均保存成数据矩阵。

4力反馈模块

把上述股骨有限元分析后的具体数据输入到力感知设备--虚拟力反馈器，应变数据包括应变的大小和方向。这里需要有碰撞检测单元，该单元用于判断操纵杆末端与股骨之间是否发生碰撞，而且可以在碰撞时对股骨模型的位置信息进行定位。位置信息定位单元用于计算质点位移的变化，并把每个质点的位移变化信息传送到反馈力计算单元当中，经力反馈计算获得力补偿信息，最后输出到操纵手柄上，这样就实现了软件与外部力反馈设备的结合。

实施例：

为了说明本系统的使用方法，下面具体举一个例子来阐述本发明的操作过程。

1.股骨有限元模型的建立

1)建立几何模型

把CT图像数据导入到几何三维重建软件中进行几何建模。本发明可采用mimics、amirad等三位建模软件，以mimics为例，可以导入CT图像的DICOM、BMP、TIFF、JPEG等多种格式的数据。对图像的蒙版进行擦除、修补、图像预处理后，计算生成出三维几何模型，如图2所示。

2)网格划分

上述的几何模型可以继续在mimics中面网格划分，单元可以是三角形，并保存成体网格划分软件(如：Ansys、Abaqus、Tetra、Hexa、GIobal等)能够识别的格式的文件。以Ansys为例，把mimics中生成的面网格文件保存成Ansys能够识别的lis格式文件。Lis格式文件导入到ANSYS后，可将面网格等二维网格通过拉伸、旋转等变换延伸成三维网格即体网格。如果面网格单元式三角形，那么体单元网格可以是四面体，如图3所示。

3)股骨弹性模量赋值

股骨的生物力学模型就是添加材料属性的有限元模型，本发明中导入的模型有两个，分别是：股骨生物力学胞元理论模型和股骨生物力学实验模型。股骨有限元网格单元的材料属性的过程由市面上常见的计算机编程语言来实现(如VC、VB、MATALB等)，因MATLAB语言的编起来比较简单，下面就以MATLAB为例来说明其具体实现的步骤：

①提取CT图像灰度值：利用命令函数dicomread编写一个循环程序读入CT图像后，用灰度门限法设置阈值，并进行预处理擦除噪点和边缘上的粘连组织，最后保存三维数据矩阵，命名为HDZ.mat。

②计算网格单元弹性模量：股骨有限元网格输出数据包括网格元编号和节点坐标，MATLAB可以通过查找节点坐标精确的对每个网格单元进行弹性模量的赋值，模型以数据矩阵的形式保存，理论模型命名为LLFEA.mat，实验模型命名为SYFEA.mat。

2.在MATLAB环境下股骨有限元分析的实现

在MATLAB环境下，首先需要设置文件存储路径和文件的读取路径。运行MATLAB软件读取和保存文件，需要将带读取和保存的文件放入指定位置。利用MATLAB程序进行有限元分析的步骤如下：

(1)将已知数据包括节点坐标、单元编号、单元密度值等保存成数据矩阵，分别命名为为fna_nodes.mat、fna_elements.mat、bgmd.mat，需要时函数调用；

(2)有限元分析中几大步骤包括单元刚度矩阵计算、组装整体刚度矩阵、应力应变计算等，将其分别编译调试成功后保存成M文件，文件名为T3D4_Stiffness.m(单元刚度矩阵)、T3D4_Assembly.m(组装整体刚度矩阵)、T3D4_Strain.m(单元应变)，需要时直接调用运行；

(3)手动输入载荷和约束条件，构建外载列阵；

(4)矩阵运算采用高斯消元法，这样有利于微分方程求解和兼顾运算速度；

(5)为了清理内存减少数据垃圾，计算过程中的过程量均保存成数据矩阵。

3.外载荷和约束的固定

本研究采用一种较精确的简化受力模型，在股骨头顶端一节点处作用垂直向下的集中载荷，模拟股骨头处的关节力，而臀肌肌群肌力、髂胫束肌力等附加肌力均不考虑。预先测量股骨模型总长，设定股骨干下端固定约束长度，即将股骨下端约束的那段所有节点的自由度均赋值为0，包括3个方向的位移和3个方向的旋转。载荷及约束情况如图6所示。

4.有限元分析结果

在股骨生物力学建模中，股骨力学特性的各向异性被简化为横观各向同性，所以只需测量横向和轴向的应变数据，最后保存载荷与应变数据。

5.与硬件的配合

本程序与虚拟力反馈器相结合。力反馈器与软件的接口相接后，将有限元分析后的股骨应变数据输入到力反馈器，虚拟力反馈器把接受到的股骨应变数据以真实力的形式输出到操纵杆上，这样操纵者就可以真实的感觉出股骨生物力学的特性。

Claims (3)

1.一种基于力反馈的股骨生物力学有限元分析系统，其特征在于：所述系统包括：

股骨CT图像灰度值提取模块(1)、股骨弹性模量赋值模块(2)、股骨生物力学有限元分析模块(3)和力反馈模块(4)；

股骨CT图像灰度值提取模块(1)，用于导入多种格式的CT图像，并把该图像存成以灰度值为元素的空间三维数据矩阵；股骨CT图像灰度值提取模块(1)还用于将获得的CT灰度值信息进行预处理：利用灰度门限法筛选出股骨，而且能手动擦除噪点和边缘上的黏连组织的数据，存留以股骨作为研究对象的三维数据矩阵；

股骨弹性模量赋值模块(2)用于获取股骨体网格单元内心坐标；

股骨弹性模量赋值模块(2)还用于根据股骨CT图像灰度值提取模块(1)得到的股骨三维数据矩阵和股骨体网格单元内心坐标计算出每个股骨体网格单元的内心灰度值；

股骨弹性模量赋值模块(2)还用于根据股骨体网格单元的内心灰度值计算得到股骨体网格单元的表观密度，进一步计算出每个股骨体网格单元的弹性模量，建立股骨有限元模型；

股骨生物力学有限元分析模块(3)用于通过每个股骨体网格单元的弹性模量计算出每个股骨体网格单元刚度矩阵，进一步地，通过给出的股骨外加载荷和固定约束得到股骨模型的整体刚度矩阵，进而求出股骨模型的每个股骨体网格单元的位移矩阵；最后通过股骨模型的每个股骨体网格单元的位移矩阵计算出股骨模型上所要测试点的应变或应力数据；

力反馈模块(4)用于根据上述应变或应力数据模拟出测试点的力学特性：力的大小和力的方向，力反馈模块(4)还用于将力的大小和力的方向传递给外部的力感知设备的操纵手柄上，使医生能真实地感知股骨模型的生物力学特性。

2.根据权利要求1所述的一种基于力反馈的股骨生物力学有限元分析系统，其特征在于：所述股骨弹性模量赋值模块(2)实现其功能的具体过程为：

1)计算四面体内心

由有限元网格模型可以得到节点坐标，通过四面体的四个节点坐标计算内心坐标；设四面体A₁A₂A₃A₄的顶点A_i所对的侧面面积为S_i(i＝1，2，3，4)，顶点A_i的坐标为(x_i，y_i，z_i)(i＝1，2，3，4)，四面体内心的坐标为(x₁，y₁，z₁)，则

$\left.\begin{array}{c}{x}_1=\frac{s_1{x}_1+s_2{x}_2+s_3{x}_3+s_4{x}_4}{s_1+s_2+s_3+s_4}\\ y_1=\frac{s_1{y}_1+s_2{y}_2+s_3{y}_3+s_4{y}_4}{s_1+s_2+s_3+s_4}\\ z_1=\frac{s_1{z}_1+s_2{z}_2+s_3{z}_3+s_4{z}_4}{s_1+s_2+s_3+s_4}\end{array}\right\}---\left(1\right)$

其中节点坐标已知，侧面面积应用海伦定理求得；

2)计算内心灰度值

利用式(2)将限元网格模型与空间三维数据矩阵进行坐标转换：

$\left.\begin{array}{c}{x}^{\&prime;}=x+a/p\\ y^{\&prime;}=y+b/p\\ z^{\&prime;}=z+c\end{array}\right\}---\left(2\right)$

式中，a、b、c为有限元网格模型的坐标轴与空间三维数据矩阵的之间的x、y、z轴上的差值，p表示像素尺寸，p由CT图像信息可知；

根据式(2)，构建起节点直角坐标和空间像素点阵之间的关系；由此可以确定每个单元内心在像素矩阵中的位置，在一个四面体网格单元当中，单元同时在空间像素矩阵中每一单元内心均由8个像素点包围，而且8个像素点位置规则；通过空间8节点正六面体的内插插值公式求得单元内心处的灰度值；

$G_i(x_i,y_i,z_i)=\underset{k=1}{\overset{8}{\mathrm{\&Sigma;}}}{N}_j(x,y,z)g(i,j,k)---\left(3\right)$

式中，G_i为所求的内心灰度值，N_j为8节点正六面体单元的形函数，g(i，j，k)为8个像素点处的灰度值；

3)计算表观密度

表观密度：单位体积的质量，指不包括髓脂时骨的质量除以试件体积得到的参变量，用ρ表示；

表观密度与CT值存在线性关系，但是由于灰度值更容易获取，所以用灰度值代替CT值计算表观密度；关系式如下：

$\mathrm{\&rho;}=\frac{{\mathrm{\&rho;}}_1-{\mathrm{\&rho;}}_0}{G_1-G_0}\&CenterDot;G---\left(5\right)$

ρ₁＝2，G₁＝2663；ρ₀＝0，G₀＝0          (6)

式中ρ₁，ρ₀为股骨最大最小表观密度，G₁，G₀为股骨最大最小灰度值；骨松质的G为100-300，骨皮质的G为2000；

计算完内心灰度值后，通过式(5)和(6)计算单元表观密度；

4)计算网格单元的弹性模量

股骨在建立各向异性力学模型时被看做是弹性体；针对股骨有限元单元，应用弹性力学基本方程，并进行一系列简化后，用应力分量表达应变分量的形式为：{ε}＝[S]·{σ}；其中[S]为柔度矩阵；

$\left[S\right]=\left[\begin{array}{cccccc}{S}_{11}& S_{12}& S_{13}& & & \\ S_{12}& S_{22}& S_{23}& & 0& \\ S_{13}& S_{23}& S_{33}& & & \\ & & & S_{44}& & \\ & 0& & & S_{55}& \\ & & & & & S_{66}\end{array}\right]---\left(11\right)$

将股骨材料看作是横观各向同性体，即柔性矩阵[S]中的力学参数只有5个，如下式：

$\left[S\right]=\left[\begin{array}{cccccc}\frac{1}{E_{\mathrm{XY}}}& -\frac{{\mathrm{\&mu;}}_{\mathrm{XY}}}{E_{\mathrm{XY}}}& -\frac{{\mathrm{\&mu;}}_{\mathrm{XY}}}{E_Z}& & & \\ -\frac{{\mathrm{\&mu;}}_{\mathrm{XY}}}{E_{\mathrm{XY}}}& \frac{1}{E_{\mathrm{XY}}}& -\frac{{\mathrm{\&mu;}}_{\mathrm{XY}}}{E_Z}& & 0& \\ -\frac{{\mathrm{\&mu;}}_{\mathrm{XZ}}}{E_{\mathrm{XY}}}& -\frac{{\mathrm{\&mu;}}_{\mathrm{XZ}}}{E_{\mathrm{XY}}}& \frac{1}{E_Z}& & & \\ & & & \frac{1}{G_{\mathrm{XY}}}& & \\ & 0& & & \frac{1}{G_{\mathrm{XZ}}}& \\ & & & & & \frac{1}{G_{\mathrm{XZ}}}\end{array}\right]---\left(12\right)$

其中，5个独立的力学参数是横向和轴向的弹性模量E，剪切模量G、横向和轴向的泊松比μ；这5个参数之间的转化关系如下：

$G_{\mathrm{XY}}=\frac{E_{\mathrm{XY}}}{2(1+{\mathrm{\&mu;}}_{\mathrm{XY}})}---\left(13\right)$

$\frac{{\mathrm{\&mu;}}_{\mathrm{ZX}}}{E_Z}=\frac{{\mathrm{\&mu;}}_{\mathrm{XZ}}}{E_{\mathrm{XY}}}---\left(14\right)$

通过上式的换算，5个独立的力学参数减少到3个；

5)股骨有限元模型，采用下述方式中的一种；

a、股骨弹性模量与表观密度关系的胞元理论模型为：

$\left\{\begin{array}{c}1007{\mathrm{\&rho;}}^2,\mathrm{\&rho;}\&le;0.25\\ 255\mathrm{\&rho;},0.25<\mathrm{\&rho;}\&le;0.4\\ 2972{\mathrm{\&rho;}}^2-933\mathrm{\&rho;},0.4<\mathrm{\&rho;}\&le;1.2\\ 1736{\mathrm{\&rho;}}^{3.25},\mathrm{\&rho;}>1.2\end{array}\right.---\left(15\right)$

式中，弹性模量的单位是MPa，表观密度的单位是g/cm³；

或b、股骨弹性模量与表观密度关系的力学实验模型：

$E=\left\{\begin{array}{c}2014{\mathrm{\&rho;}}^{2.5},\mathrm{\&rho;}\&le;1.2\\ 1736{\mathrm{\&rho;}}^{3.5},\mathrm{\&rho;}>1.2\end{array}\right.---\left(16\right)$

3.根据权利要求2所述的一种基于力反馈的股骨生物力学有限元分析系统，其特征在于：所述股骨生物力学有限元分析模块(3)实现其功能的具体过程为：

1)结构离散化：

将其离散成有限个单元体，利用有限个单元体代替原来的组织结构进行计算；

2)选择形函数：

选择线性的位移函数作为位移关系：

{f}＝[N]{δ}^ε          (17)

式中：{f}表示单元内任意一点的位移列阵；{δ}^ε表示节点位移列阵；[N]表示形状函数矩阵；

3)单元力学特性：

计算方法如下：

首先应用几何方程建立位移和应变的关系：

{ε}＝[B]{δ}^ε          (18)

式中，{ε}是单元内任意一点的应变列阵，[B]是系数矩阵；

然后将上式代入广义胡克定律关系式(物理方程)得到：

{σ}＝[D][B]{δ}^ε          (19)

式中，{σ}是单元内任意一点的应力列阵，[D]是与材料有关的弹性矩阵；

再由平衡微分方程得到力和位移的关系：

{r}^ε＝[K]^ε{δ}^ε          (20)

式中，[K]^e就是单元刚度矩阵，它的计算公式如下：

[K]^e＝∫∫∫[B]^T[D][B]dxdydz          (21)

4)整体平衡方程：

在整体上建立组织结构的平衡方程为：

[K]{δ}={R}         (22)

式中，[K]为整体刚度矩阵，{δ}为整体节点位移，{R}为整体节点载荷；

5)解方程：

解方程后得出股骨有限元模型上目标测试点的应力与应变的数据。

Images