CN102262699A - 基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法 - Google Patents

Abstract

本发明涉及一种物体形变实时模拟图形处理技术，特别是涉及一种基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法。在前处理过程中，为软组织建立线性粘弹性生物力学模型；在计算形变过程中，根据软组织所承载的载荷动态划分无网格区域和质点弹簧区域，并在无网格区域与质点弹簧区域之间的连接区域建立过渡单元，构造过渡单元近似位移函数，实现无网格伽辽金方法与质点弹簧方法的自适应耦合；在后处理过程中，根据形变计算结果，将形变过程每个时间步长的质点或节点的状态输出到屏幕上，并进行光照渲染，最终在屏幕上显示软组织器官在受力情况下的实时形变过程，实现动态形变可视化效果。利用质点弹簧法效率高和无网格伽辽金法精度高、无需重构网格的优势，弥补伽辽金法不适宜求解大规模问题的缺陷，从而有效降低软组织形变仿真中的计算复杂度，提高运算效率。

Description

基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法

技术领域

本发明涉及一种物体形变实时模拟图形处理技术，特别是涉及一种基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法。

背景技术

人体软组织形变计算模型的研究可以追溯到上世纪80年代，早期的形变模型来自于计算机辅助几何设计(CAGD)领域，是一些运用纯几何技术的非物理模型，这种模型不考虑真实形变的物理规律，一年后，被Terzopoulos提出的基于物理特征的形变模型所替代。目前，基于物理特性的形变计算模型主要分为质点弹簧模型(Mass-Spring)、有限元模型(Finite-Element Model，简写为FEM)和边界元模型(Boundary Element Model，简写为BEM)三大类。质点弹簧模型将物体离散成若干个点，使用有质量的点来表达物体质量，用无质量的弹簧和阻尼器来表达点之间的相互影响。该模型原理简单、计算量小，实时性好，被学者广泛用于软组织的变形、切割、缝合等虚拟手术仿真。但是该模型本身也具有其自身无法克服的缺点，由于假定软组织之间用弹簧连接，而弹簧弹性系数的设置没有理论依据，只能来源于实验员的经验和不断地调试，并且在软组织大变形时会出现失真问题，因此，能够克服这一缺点的有限元模型逐渐成为研究的热点。

有限元法是求解弹性力学问题的经典方法，它的基本思想是将连续的求解区域离散为一组有限个、且按一定方式相互联结在一起的单元的组合体。用每个单元内所假设的近似函数来分片地表示求解域内待定的应变量，利用变分原理建立求解应变量的代数方程组，从而计算形变。它的优点是具有较高的精度，适应性较强，能够适用于各种几何结构与材料特性的仿真。它的缺点是：对于三维复杂结构的前处理即三维网格构建困难，耗时较长；其计算结果通常位移连续，应力、应变不连续，需要后处理修匀；当结构材料不可压缩时出现“体积自锁”现象；在大变形时会出现网格畸变问题。由于软组织器官具有不可压缩性，所以在使用有限元模型时会出现“体积自锁”现象；在虚拟手术中，经常发生大变形情况或切割、缝合操作，这时器官的拓扑结构会发生变化，利用有限元法需要不断进行网格重构，不仅耗时耗力，而且容易出现网格畸变问题；由于有限元计算的应力、应变结果不连续，需要后处理修匀，会影响虚拟手术中力反馈的精确性和实时性。由此可见，有限元法并不是软组织形变计算的最佳方案，需要探索其它更适合的计算方法。

边界元法是指以控制方程的基本解为基础，将区域的边界问题化为边界面上的方程，然后在边界面上划分单元，再用Galerkin法或其它数值计算方法求解。它能够把三维问题变为二维问题，二维问题变为一维问题，可以使求解的自由度下降，但由于边界元法也是基于网格的数值方法，故也可能出现网格畸变和扭曲。

常用的有网格物理模型需要耗费很大的精力来构造网格模型，后续的计算过程大都紧密依赖于这种网络结构；在软组织大变形时有可能会发生网格畸变；在发生切割、缝合等拓扑改变时，需要重新构造网格；因此有网格方法难以准确描述软组织大变形和拓扑改变情况。为解决这些问题，自然产生了在数值处理过程中摆脱网格的想法，无网格法应运而生。无网格法将有限元法中的网格结构去除，完全代之以一系列的节点排列，不需利用预定义的网格对场变量进行插值和近似，而是基于点的近似，采用权函数来表征节点信息，并且某个域上的节点可以影响研究对象上任何一点的力学特性。这样，摆脱了不连续性对问题的束缚(如网格重构等)，保证了求解的精度，特别适合虚拟手术中切割、缝合操作引起的拓扑改变。在无网格方法中，伽辽金法只需要节点信息和边界描述就可以将其弱式形式转换成一组代数系统方程，并采用一种光滑的能与各取样点的值达到最佳近似的移动最小二乘近似函数作为位移近似函数，以此来表征节点信息，因此计算精度远高于无网格强式法，非常稳定。

但是，无网格伽辽金的缺点是计算量较大，大规模的使用伽辽金方法会大大增加计算时间，因此设计一个性能更加良好的模型来实现软组织的形变仿真，并尽可能的满足关于实时性、健壮性、精确性的要求，已成为软组织形变仿真面临的首要问题。实际上，在虚拟手术中只有在手术器械与器官的接触区域才会产生大的变形及切割缝合等拓扑改变，在其它区域形变较小，因此没有必要在整个区域使用无网格伽辽金方法。因此，我们采用无网格伽辽金与质点弹簧耦合的方法，在手术器械与器官的接触区域采用无网格伽辽金方法，在其它区域采用质点弹簧方法。

发明内容

本发明针对现有技术不足，提出一种基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，将无网格伽辽金与质点弹簧耦合来解决虚拟手术中软组织形变仿真问题，目的是利用质点弹簧法效率高和无网格伽辽金法精度高、无需重构网格的优势，弥补伽辽金法不适宜求解大规模问题的缺陷，从而有效降低软组织形变仿真中的计算复杂度，提高运算效率，较好解决软组织形变仿真中精确性与实时性的矛盾。

本发明所采用的技术方案：

一种基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，包括前处理过程、计算形变过程、后处理过程三个步骤，在前处理过程中，为软组织建立线性粘弹性生物力学模型；在计算形变过程中，根据软组织所承载的载荷动态划分无网格区域和质点弹簧区域，并在无网格区域与质点弹簧区域之间的连接区域建立过渡单元，构造过渡单元近似位移函数，实现无网格伽辽金方法与质点弹簧方法的自适应耦合；在后处理过程中，根据形变计算结果，将形变过程每个时间步长的质点或节点的状态输出到屏幕上，并进行光照渲染，最终在屏幕上显示软组织器官在受力情况下的实时形变过程，实现动态形变可视化效果；所述计算形变的过程，具体包括如下四个步骤：

1)设计载荷与距离之间的函数，作为划分无网格区域的依据，并使无网格区域足够大以保证不连续边界都在无网格区域内，其它区域为质点弹簧区域；

2)对于无网格区域和质点弹簧区域，分别建立有效数据结构，分类管理数据；

3)在无网格区域和质点弹簧区域之间的连接区域建立过渡单元：

(1)将已划分的无网格区域与质点弹簧区域作为两个实体，两实体相接触部分作为过渡边界，在背景网格的基础上，细分过渡边界处的背景网格并作为子单元，使每个子单元中最多有一个节点或质点存在，以不包含节点或质点的单元作为空子单元；

(2)以空子单元作为搜索对象，分别向上下左右前后六个方向搜索子单元，经过迭代，逐步将无网格区域内部符合转换条件的空子单元转化为无网格单元，将质点弹簧区域内部符合转换条件的空子单元转化为质点弹簧单元，直至空子单元不再变化，并将剩余的空子单元作为过渡单元；

4)在过渡单元内建立过渡节点，确定过渡节点的近似位移函数，实现无网格区域与质点弹簧区域之间的光滑过渡；过渡节点的建立满足两个条件：

(1)过渡节点分别受到无网格区域和质点弹簧区域的作用力，其合力为零；

(2)过渡节点分别受到无网格区域和质点弹簧区域的作用力，其位移相等；

根据上述条件，分别构建基于无网格伽辽金和质点弹簧的线性粘弹性动态运动方程，求解微分方程，得出同一过渡节点的两个不同近似位移，遍历所有过渡节点，得到两种模型构建出的近似位移函数，建立二者之间的函数关系，作为过渡节点的近似位移函数。

所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，无网格区域内部空子单元转化为无网格单元的转换条件为：

1)上下、左右、前后三组方向中存在两组及两组以上方向的子单元为无网格单元；

2)上下、左右、前后六个方向上的六个子单元都为空子单元；

3)对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，每组方向取一个子单元，共三个方向上的三个子单元为无网格单元；

4)对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中一组方向上的子单元都取到，另两组方向中每组取其中一个方向上的子单元，共计四个方向上的四个子单元为无网格单元；

5)对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中两组方向上的子单元都取到，另一组方向上取其中一个方向上的子单元，共计五个方向上的五个子单元为无网格单元。

所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，将质点弹簧区域内部符合转换条件的空子单元转化为质点弹簧单元，其转换条件为：

1)上下、左右、前后三组方向中存在两组及两组以上方向的子单元为质点弹簧单元；

2)上下、左右、前后六个方向上的六个子单元都为空子单元；

3)对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，每组方向取一个子单元，共三个方向上的三个子单元为质点弹簧单元；

4)对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中一组方向上的子单元都取到，另两组方向中每组取其中一个方向上的子单元，共计四个方向上的四个子单元为质点弹簧单元；

5)对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中两组方向上的子单元都取到，另一组方向上取其中一个方向上的子单元，共计五个方向上的五个子单元为质点弹簧单元。

所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，剩余的空子单元作为过渡单元，是指：

1)不符合无网格单元转换条件，也不符合质点弹簧单元转换条件的空子单元；

2)上下、左右、前后六个方向上的六个子单元同时包含无网格单元和质点弹簧单元，或者同时包含无网格单元、质点弹簧单元和空子单元。

所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，在前处理过程中，假定软组织为均匀、各向同性材料，具有准不可压缩性、线性粘弹性特征，并确定初始化数据，定义初始条件、位移边界条件和应力边界条件；然后用弹簧模拟粘弹性体的弹性，阻尼器模拟其粘性，并假设弹簧的变形与负载成比例，而阻尼器的变形与变形的速度成比例；其次各向同性材料的变形可以分离成体积改变和等体积的形状畸变两部分，故把应力和应变张量分解成它的球形张量和偏斜张量，进而构造Kelvin模型的本构方程，模拟软组织应力松弛和蠕变特性；最后建立几何方程、平衡方程，并施加边界条件和初始条件，为软组织建立线性粘弹性生物力学模型。

所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，在计算形变过程中，首先生成初始的背景网格，根据外部载荷划分无网格区域和质点弹簧区域，建立从无网格到质点弹簧数据之间的数据结构，以降低自适应过程中数据重用现象发生的频率，其次在连接区域建立过渡单元，在过渡单元内建立过渡节点及过渡节点的近似位移函数，实现无网格区域与质点弹簧区域之间的光滑过渡，然后施加边界条件，结合Kelvin标准线性粘弹性力学模型，采用增量法求解位移，根据位移求应力、应变。

本发明的有益积极效果：

1、在大变形和拓扑改变区域用无网格伽辽金方法，不需要节点的连接信息，有效避免目前常用的有限元和边界元方法中的网格畸变与重构，非常适合手术仿真中经常发生大变形和拓扑改变的情况；无网格伽辽金方法具有良好连续性和形式灵活的场函数，具有较高的稳定性和计算精度；无网格伽辽金方法的应力应变结果具有较好的光滑性，能够很好支持力反馈设备。

2、采用基于无网格伽辽金与质点弹簧耦合方法，在大变形和拓扑改变区域采用高精度的无网格伽辽金方法，在其它形变较小区域采用高效率的质点弹簧方法，不仅可发挥上述无网格伽辽金方法的诸多优点，而且提高了形变仿真效率。

3、使用基于Kelvin标准线性粘弹性模型可有效模拟软组织在形变过程中所表现出的蠕变和松弛特性。这种基于粘弹性的耦合计算模型较目前普遍使用的弹性力学模型更加符合人体软组织的生物力学特性。

附图说明

图1为本发明基于无网格伽辽金与质点弹簧耦合的软组织形变仿真执行方法流程图；

图2为本发明动态划分无网格区域与质点弹簧区域的剖面图，图中，1为手术器械影响区域，该区域使用无网格伽辽金模型计算软组织形变，2为过渡区域，3为质点弹簧区域，该区域使用质点弹簧模型计算软组织形变；

图3为本发明动态划分无网格区域与质点弹簧区域的立体图；

图4为本发明的无网格伽辽金与质点弹簧耦合模型的剖面图；

图5为肝脏器官分别在受到拉力和压力情况下的形变效果图。

具体实施方式

实施例一：参见图1，本发明基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，包括前处理过程、计算形变过程、后处理过程三个部分：

在前处理过程中，为软组织建立线性粘弹性生物力学模型。假定软组织为均匀、各向同性材料，具有准不可压缩性、线性粘弹性特征，构造本构方程、几何方程和平衡方程，施加边界条件和初始条件；

在计算形变过程中，根据软组织所承载的载荷动态划分无网格区域和质点弹簧区域，并在无网格区域与质点弹簧区域之间的连接区域建立过渡单元，构造过渡单元近似位移函数，实现无网格伽辽金方法与质点弹簧方法的自适应耦合；

在后处理过程中，根据形变计算结果，将形变过程每个时间步长的质点或节点的状态输出到屏幕上，并进行光照渲染，最终在屏幕上显示软组织器官在受力情况下的实时形变过程，实现动态形变可视化效果。

实施例二：本实施例基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，与实施例一不同的是：计算形变的过程具体包含如下四个步骤：

1)设计载荷与距离之间的函数，作为划分无网格区域的依据，并使无网格区域足够大以保证不连续边界都在无网格区域内，其它区域为质点弹簧区域；

2)对于无网格区域和质点弹簧区域，分别建立有效数据结构，分类管理数据；

3)在无网格区域和质点弹簧区域之间的连接区域建立过渡单元；

4)在过渡单元内建立过渡节点及过渡节点的近似位移函数，实现无网格区域与质点弹簧区域之间的光滑过渡。

实施例三：本实施例基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，与实施例二不同的是：在无网格区域和质点弹簧区域之间的连接区域建立过渡单元，包含如下两个步骤：1)将已划分的无网格区域与质点弹簧区域作为两个实体，两实体相接触部分作为过渡边界，在背景网格的基础上，细分过渡边界处的背景网格并作为子单元，使每个子单元中最多有一个节点或质点存在，以不包含节点或质点的单元作为空子单元；

2)以空子单元作为搜索对象，分别向上下左右前后六个方向搜索子单元，经过迭代，逐步将无网格区域内部符合转换条件的空子单元转化为无网格单元，将质点弹簧区域内部符合转换条件的空子单元转化为质点弹簧单元，直至空子单元不再变化，并将剩余的空子单元作为过渡单元。

本发明基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，在计算形变过程中，首先确定无网格区域和质点弹簧区域，将两个区域中的节点或质点信息分别存储；在无网格区域和质点弹簧区域之间的连接区域细分背景网格，保证每个子单元内最多只包含一个节点或质点；然后，以空子单元作为搜索对象，应用搜索算法将部分空子单元转换为无网格单元或质点弹簧单元，并将剩余的空子单元作为过渡单元。

实施例四：本实施例基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，与实施例三不同的是：在过渡单元内建立过渡节点，按照如下方法确定过渡节点的近似位移函数。过渡节点满足两个条件：

1)过渡节点分别受到无网格区域和质点弹簧区域的作用力，其合力为零；

2)过渡节点分别受到无网格区域和质点弹簧区域的作用力，其位移相等；

根据上述条件，分别构建基于无网格伽辽金和质点弹簧的线性粘弹性动态运动方程，求解微分方程，得出同一过渡节点的两个不同近似位移，遍历所有过渡节点，得到两种模型构建出的近似位移函数，建立二者之间的函数关系，并作为过渡节点的近似位移函数。

实施例五：参见图1，本实施例基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，前处理过程、计算形变过程、后处理过程分别如下：

前处理过程：首先为软组织建立线性粘弹性生物力学模型。假定软组织为均匀、各向同性材料，具有准不可压缩性、线性粘弹性特征，用弹簧模拟粘弹性体的弹性，阻尼器模拟其粘性，并假设弹簧的变形与负载成比例，而阻尼器的变形与变形的速度成比例；其次各向同性材料的变形可以分离成体积改变和等体积的形状畸变两部分，故把应力和应变张量分解成它的球形张量和偏斜张量，进而构造Kelvin模型的本构方程，较好的模拟软组织应力松弛和蠕变特性；最后建立几何方程、平衡方程，并施加边界条件和初始条件。

计算形变过程：手术器械接触软组织时，不管是按压还是提拉或者切割，均根据手术器械接触的器官局部表面所受的载荷，设计载荷与距离之间的函数关系，从而动态划分无网格区域和质点弹簧区域，使无网格区域和质点弹簧区域满足预定的范围，并且不连续边界最终仍在无网格区域内，连续边界都在质点弹簧区域内；其次分别在无网格区域与质点弹簧区域建立有效的数据结构，并分类管理数据；然后将已划分的无网格区域与质点弹簧区域作为两个实体，两实体相接触部分作为过渡边界，细分过渡边界处的背景网格，并将细分后的背景网格作为子单元，子单元满足至多存在一个节点或质点的条件，这样得到两个实体的边界；由于细分后的背景网格存在不包含节点或质点的单元，故将不包含节点或质点的单元定义为空子单元；以空子单元作为搜索对象，选用先上后下，从左到右，从前到后的顺序逐步搜索子单元，通过迭代，逐步将符合转换条件的空子单元转化为无网格单元或质点弹簧单元，直至空子单元不再变化，并将剩余的空子单元作为无网格和质点弹簧之间的过渡单元；最后在过渡单元内建立过渡节点，根据过渡节点满足的位移和力平衡条件，建立过渡节点的近似位移函数。

所述的无网格伽辽金模型，用来计算手术器械与软组织接触的形变较大区域，采用Pascal图建立三维基函数p(x)，确定基函数的个数m；构造全局弱式法的背景网格，设定积分点x_Q为影响域中的计算点x，计算影响域的尺寸d_s，确定参与构造计算点的MFree形函数的场节点个数n，保证m＝n；利用影响域所选择的场节点，建立权函数W及其导数；采用移动最小二乘法构造影响域中所有场节点的节点形函数φ(x)及其偏导数

由粘弹性力学控制方程对应的能量泛函得到伽辽金弱变分形式，并进行总体离散，取驻值，得EFG法的总体离散控制方程，施加本质边界条件，利用二维等带宽存储的高斯循环消去法获得离散的节点位移。

所述的质点弹簧模型，用来计算形变较小区域，按照胡克定律，建立单个质点运动的拉格朗日运动方程：

$m\frac{{\∂}^2{u}^{\mathrm{MS}}}{\∂t^2}+\mathrm{\γ}\frac{\∂u^{\mathrm{MS}}}{\∂t}+f_{\mathrm{int}}^i=f_{\mathrm{ext}}^i$

式中，u^MS表示质点i的位置矢量，m表示质点i的质量，γ表示相邻质点间的黏性密度，

表示弹簧对质点i所施加的内力，

表示质点i所受的外力，根据质点所受外力，解微分方程得质点位移。

所述的无网格伽辽金与质点弹簧耦合模型，用来实现无网格区域与质点弹簧区域之间的自然过渡。在过渡单元内建立的过渡节点满足两个条件：1)过渡节点分别受到无网格区域和质点弹簧区域的作用力，其合力为零；2)过渡节点分别受到无网格区域和质点弹簧区域的作用力，其位移相等。现根据上述条件，分别构建基于无网格和质点弹簧的线性粘弹性动态仿真运动方程，求解微分方程得同一过渡节点的两个不同近似位移，遍历所有过渡节点，得到使用两种模型构建出的近似位移函数，建立二者之间的函数关系，并作为过渡节点的近似位移函数。

假定t＝0时刻的位移、速度和加速度已知，将时间求解域0～T等分为n个时间间隔Δt(＝T/n)，利用中心差分法或Newmark法计算每一时间步长的有效载荷和位移，根据位移求应变，并由应变计算新的体积应变和形状畸变，根据粘弹性本构方程，计算体积应力、偏应力和质点应力；

后处理过程：根据形变计算的结果进行光照的渲染，最终在屏幕上实时绘制出每个时间步长的质点或节点的状态，这样，就可以实现动态形变可视化效果。

实施例六：本实施例结合图1～图5，对本发明技术方案作详细说明：

如图1所示，当模拟软组织形变时，主要分为三大模块进行实施：

1、前处理过程

首先，为软组织建立线性粘弹性生物力学模型。假定软组织为均匀、各向同性材料，具有准不可压缩性、线性粘弹性特征，该模型由Maxwell模型与线性弹簧并联组合而成，并且假定线性弹簧即刻产生与载荷成正比的变形，阻尼器可于任一瞬间产生与载荷成正比的速度，因此本构方程为：

$\mathrm{\σ}+\frac{\mathrm{\η}}{E_2}\stackrel{\·}{\mathrm{\σ}}=E_1\mathrm{\ϵ}+\mathrm{\η}(1+\frac{E_1}{E_2})\stackrel{\·}{\mathrm{\ϵ}}$

式中σ表示应力，ε表示应变，E表示线性弹簧的弹性系数，E₁表示与Maxwell模型并联的线性弹簧的弹性系数，E₂表示Maxwell模型中线性弹簧的弹性系数，η表示阻尼器的粘性系数，和

分别表示应力与应变关于时间的导数。

Kelvin标准线粘弹性模型的材料参数获取较为容易，且可以实现活体测量，因此能够通过调整参数来描述软组织发生形变时所具有的松弛和蠕变特性，其中蠕变函数为：

$C\left(t\right)=\frac{1}{E_1}[1-(1-\frac{{\mathrm{\τ}}_{\mathrm{\ϵ}}}{{\mathrm{\τ}}_{\mathrm{\σ}}}{e}^{-t/{\mathrm{\τ}}_{\mathrm{\σ}}})]I\left(t\right)$

松弛函数为：

$K\left(t\right)=E_1[1-(1-\frac{{\mathrm{\τ}}_{\mathrm{\σ}}}{{\mathrm{\τ}}_{\mathrm{\ϵ}}}{e}^{-t/{\mathrm{\τ}}_{\mathrm{\ϵ}}})]I\left(t\right)$

式中 ${\mathrm{\τ}}_{\mathrm{\ϵ}}=\frac{\mathrm{\η}}{E_2},$ ${\mathrm{\τ}}_{\mathrm{\σ}}=\frac{\mathrm{\η}}{E_1}(1+\frac{E_1}{E_2}),$

一般情况下，各向同性材料的变形可以分离成体积改变和等体积的形状畸变两部分，即体积应力只改变体积，应力偏量导致等体积的形状畸变，故把应力和应变张量分解成它的球形张量和偏斜张量部分，即

${\mathrm{\σ}}_{\mathrm{mn}}=s_{\mathrm{mn}}+\frac{1}{3}{\mathrm{\δ}}_{\mathrm{mn}}{\mathrm{\σ}}_{\mathrm{ii}},$ ${\mathrm{\ϵ}}_{\mathrm{mn}}=e_{\mathrm{mn}}+\frac{1}{3}{\mathrm{\δ}}_{\mathrm{mn}}{\mathrm{\ϵ}}_{\mathrm{ii}},$ m，n＝1，2，3

其中，δ_mn为Kronecker符号；ε_ii和σ_ii分别为体积应变和体积应力；e_mn和s_mn分别为应变偏量和应力偏量的分量；e_mm＝0和s_mn＝0。

由于两种情形下的粘弹特性与效应可以分别考虑，故根据Kelvin标准线性粘弹性力学模型，偏应力张量和偏应变张量之间、体积应力和体积应变之间的三维粘弹本构方程为：

$s_{\mathrm{mn}}+\frac{\mathrm{\η}}{E_2}{\stackrel{\·}{s}}_{\mathrm{mn}}=E_1{e}_{\mathrm{mn}}+\mathrm{\η}(1+\frac{E_1}{E_2}){\stackrel{\·}{e}}_{\mathrm{mn}},$ m，n＝1，2，3

${\mathrm{\σ}}_{\mathrm{ii}}+\frac{\mathrm{\η}}{E_2}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{ii}}=E_1{\mathrm{\ϵ}}_{\mathrm{ii}}+\mathrm{\η}(1+\frac{E_1}{E_2}){\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ii}}$

其中

分别为偏应变率和体积应变率。

2、计算形变过程

1)建立无网格伽辽金模型

所述的无网格伽辽金模型，用来计算手术器械与软组织接触的形变较大区域，采用Pascal图建立三维基函数p(x)，并确定基函数的个数m，

p^T(x)＝[1 x y z]^T m＝4，p＝1，为线性基函数；

p^T(x)＝[1 x y z x² xy y² yz z² zx]^T m＝10，p＝2，为二次基函数；

p^T(x)＝[1 x y z x² xy y² yz z² zx x³ x²y xy² y³ y²z yz² z³ z²x zx² xyz]^T m＝20，p＝3，为三次基函数；

构造全局弱式法的背景网格，设定积分点x_Q为影响域中的计算点x，通过设定影响域的尺寸d_s

$d_s=a_s{d}_c,d_c=\frac{\sqrt[3]{V_s}}{\sqrt[3]{n_{v_s}}-1}$

确定参与构造计算点的MFree形函数的场节点n，保证m＝n。式中的a_s为该影响域的无量纲尺寸，d_c为位于点x_Q附近的节点间距，V_s为预估计影响域的体积；

为包含在体积V_s中的预估计影响域中节点数。

利用影响域所选择的场节点，建立权函数W及其导数，现以三次样条函数作为权函数为例：

$W_i\left(x\right)=\left\{\begin{array}{cc}2/3-4r_i^2+4r_i^3& r_i\&le;0.5\\ 4/3-4r_i+4r_i^2-4/3r_i^3& 0.5<r_i\&le;1\\ 0& r_i>1\end{array}\right.$

式中d_i＝|x-x_i|为节点x_i与计算点x之间的距离，r_w为权函数支持域尺寸。

对权函数在各方向求导，得权函数的一阶导数为：

$W_i\left(x\right)=\left\{\begin{array}{cc}-8r_i+12r_i^2& r_i\&le;0.5\\ -4+8r_i-4r_i^2& 0.5<r_i\&le;1\\ 0& r_i>1\end{array}\right.$

采用移动最小二乘法构造影响域中所有场节点的节点形函数φ(x)及其偏导数

由粘弹性力学控制方程对应的能量泛函得到伽辽金弱变分形式，

$\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{\&delta;}{\left(\mathrm{Lu}\right)}^{T}D\left(\mathrm{Lu}\right)\mathrm{d\&Omega;}-\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{\&delta;}{u}^T\mathrm{bd\&Omega;}-\underset{{\mathrm{\&Gamma;}}_t}{\&Integral;}\mathrm{\&delta;}{u}^T\stackrel{\&OverBar;}{t}\mathrm{d\&Gamma;}-\mathrm{\&delta;}\underset{{\mathrm{\&Gamma;}}_u}{\&Integral;}\frac{1}{2}{(u-\stackrel{\&OverBar;}{u})}^T\mathrm{\&alpha;}(u-\stackrel{\&OverBar;}{u})\mathrm{d\&Gamma;}=0$

将上述方程进行总体离散，取驻值，得EFG法的总体离散控制方程，

[K+K^α]U＝F+F^α

式中，U为整个问题域中所有节点的节点位移参数，K为由节点刚度矩阵组装而成的总体刚度矩阵，F为由节点力向量组装而成的总体外力向量，K^α为由总体惩罚刚度矩阵引起的附加矩阵，F^α为由本质边界条件所引起的附加力向量，K、K^α、F和F^α通过背景网格中的高斯积分点做三维高斯积分或Irons积分得到；施加本质边界条件，利用二维等带宽存储的高斯循环消去法获得离散的节点位移。

2)建立质点弹簧模型

所述的质点弹簧模型，用来计算形变较小区域，按照胡克定律，建立单个质量点运动的拉格朗日运动方程

$m\frac{{\&PartialD;}^2{u}^{\mathrm{MS}}}{\&PartialD;t^2}+\mathrm{\&gamma;}\frac{\&PartialD;u^{\mathrm{MS}}}{\&PartialD;t}+f_{\mathrm{int}}^i=f_{\mathrm{ext}}^i$

其中：u^MS表示质点i的位置矢量；m表示质点i的质量，γ表示相邻质点间的黏性密度，

表示弹簧对质点i所施加的内力，

表示质点i所受的外力，根据质点所受外力，解偏微分方程得质点位移。

3)建立无网格伽辽金与质点弹簧耦合模型

对于一个软组织来讲，当某一位置受力，则与受力点距离较近的点发生的形变较剧烈，而远离受力点的组织发生的形变比较微小。为了保证计算精度，提高计算效率，形变较大区域采用无网格伽辽金模型，形变较小区域采用质点弹簧模型。现根据手术器械所接触的器官局部表面所承受的载荷，设计载荷与距离之间的函数，并将其作为划分无网格区域与质点弹簧区域的依据，其中无网格区域和质点弹簧区域要满足预定的范围，以保证不连续边界在无网格区域内，连续边界在质点弹簧区域内；其次在无网格区域和质点弹簧区域，分别建立有效数据结构，并分类管理数据；然后将已划分的无网格区域与质点弹簧区域作为两个实体，两实体相接触部分作为过渡边界，在背景网格的基础上，细分过渡边界处的背景网格，并将细分后的背景网格作为子单元，子单元满足的条件为单元中至多一个节点或质点存在；定义不包含节点或质点的单元为空子单元，现以空子单元作为搜索对象，按照上下左右前后的顺序搜索子单元，经过迭代，逐步将无网格区域内部符合转换条件的空子单元转化为无网格单元，将质点弹簧区域内部符合转换条件的空子单元转化为质点弹簧单元，并将剩余的空子单元作为过渡单元，其中

(1)无网格区域内部的空子单元转化为无网格单元的转换条件为：

①上下、左右、前后三组方向中存在两组及两组以上方向的子单元为无网格单元；

②上下、左右、前后六个方向上的六个子单元都为空子单元；

③对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，每组方向取一个子单元，共三个方向上的三个子单元为无网格单元；

④对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中一组方向上的子单元都取到，另两组方向中每组取其中一个方向上的子单元，共计四个方向上的四个子单元为无网格单元；

⑤对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中两组方向上的子单元都取到，另一组方向上取其中一个方向上的子单元，共计五个方向上的五个子单元为无网格单元；

(2)质点弹簧区域内部的空子单元转化为质点弹簧单元的转换条件为：

①上下、左右、前后三组方向中存在两组及两组以上方向的子单元为质点弹簧单元；

②上下、左右、前后六个方向上的六个子单元都为空子单元；

③对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，每组方向取一个子单元，共三个方向上的三个子单元为质点弹簧单元；

④对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中一组方向上的子单元都取到，另两组方向中每组取其中一个方向上的子单元，共计四个方向上的四个子单元为质点弹簧单元；

⑤对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中两组方向上的子单元都取到，另一组方向上取其中一个方向上的子单元，共计五个方向上的五个子单元为质点弹簧单元。

(3)剩余的空子单元的特征为：

①不符合无网格单元转换条件，也不符合质点弹簧单元转换条件的空子单元；

②上下、左右、前后六个方向上的六个子单元同时包含无网格单元和质点弹簧单元，或者同时包含无网格单元、质点弹簧单元和空子单元。

确定过渡单元后，在过渡单元内建立过渡节点，并且过渡节点满足两个平衡条件：

(1)过渡节点分别受到无网格区域和质点弹簧区域的作用力，其合力为零；

(2)过渡节点分别受到无网格区域和质点弹簧区域的作用力，其位移相等；

根据上述条件，分别构建基于无网格伽辽金和质点弹簧的线性粘弹性动态仿真运动方程，求解微分方程，得出同一过渡节点的两个不同近似位移，遍历所有过渡节点，得到两种模型构建出的近似位移函数，建立二者之间的函数关系，并作为过渡节点的近似位移函数，实现无网格区域与质点弹簧区域之间的光滑过渡。

最后假定t＝0时刻的位移、速度和加速度已知，将时间求解域0～T等分为n个时间间隔Δt(＝T/n)，利用中心差分法或Newmark法计算每一时间步长的有效载荷和位移，根据位移求应变，并由应变计算新的体积应变和形状畸变，根据Kelvin标准线粘弹性生物力学模型的本构方程，计算体积应力、偏应力和质点应力。

3、后处理过程

根据形变计算的结果进行光照渲染，最终在屏幕上实时绘制出每个时间步长的质点或节点的状态，实现动态形变可视化效果。

如图2、3、4所示，首先根据手术器械所接触的器官局部表面所承受的载荷，设计载荷与距离之间的函数，动态划分无网格区域与质点弹簧区域，然后在无网格区域与质点弹簧区域之间建立一系列的过渡单元，最后构造出相应的位移近似函数，把两种方法耦合到一起。当对肝脏施加向上或向下的载荷时，通过基于无网格伽辽金与质点弹簧耦合的软组织形变计算模型，逼真的模拟了手术区域软组织形变效果，如图5所示。相对于目前普遍采用的有限元模型，本模型降低了整体复杂度，加快了计算速度。此外，本模型便于实现切割、缝合等手术操作仿真，能够有力支持力反馈设备，可以作为虚拟手术中软组织物理模型。

Claims (6)

1.一种基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，包括前处理过程、计算形变过程、后处理过程三个步骤，在前处理过程中，为软组织建立线性粘弹性生物力学模型；在计算形变过程中，根据软组织所承载的载荷动态划分无网格区域和质点弹簧区域，并在无网格区域与质点弹簧区域之间的连接区域建立过渡单元，构造过渡单元近似位移函数，实现无网格伽辽金方法与质点弹簧方法的自适应耦合；在后处理过程中，根据形变计算结果，将形变过程每个时间步长的质点或节点的状态输出到屏幕上，并进行光照渲染，最终在屏幕上显示软组织器官在受力情况下的实时形变过程，实现动态形变可视化效果；其特征是：所述计算形变的过程，具体包括如下四个步骤：

1）设计载荷与距离之间的函数，作为划分无网格区域的依据，并使无网格区域足够大以保证不连续边界都在无网格区域内，其它区域为质点弹簧区域；

2）对于无网格区域和质点弹簧区域，分别建立有效数据结构，分类管理数据；

3）在无网格区域和质点弹簧区域之间的连接区域建立过渡单元：

（1）将已划分的无网格区域与质点弹簧区域作为两个实体，两实体相接触部分作为过渡边界，在背景网格的基础上，细分过渡边界处的背景网格并作为子单元，使每个子单元中最多有一个节点或质点存在，以不包含节点或质点的单元作为空子单元；

（2）以空子单元作为搜索对象，分别向上下左右前后六个方向搜索子单元，经过迭代，逐步将无网格区域内部符合转换条件的空子单元转化为无网格单元，将质点弹簧区域内部符合转换条件的空子单元转化为质点弹簧单元，直至空子单元不再变化，并将剩余的空子单元作为过渡单元；

4）在过渡单元内建立过渡节点，确定过渡节点的近似位移函数，实现无网格区域与质点弹簧区域之间的光滑过渡；过渡节点的建立满足两个条件：

（1）过渡节点分别受到无网格区域和质点弹簧区域的作用力，其合力为零；

（2）过渡节点分别受到无网格区域和质点弹簧区域的作用力，其位移相等； 

根据上述条件，分别构建基于无网格伽辽金和质点弹簧的线性粘弹性动态运动方程，求解微分方程，得出同一过渡节点的两个不同近似位移，遍历所有过渡节点，得到两种模型构建出的近似位移函数，建立二者之间的函数关系，作为过渡节点的近似位移函数。

2.根据权利要求1所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，其特征是：无网格区域内部空子单元转化为无网格单元的转换条件为： 

1）上下、左右、前后三组方向中存在两组及两组以上方向的子单元为无网格单元；

2）上下、左右、前后六个方向上的六个子单元都为空子单元；

3）对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，每组方向取一个子单元，共三个方向上的三个子单元为无网格单元；

4）对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中一组方向上的子单元都取到，另两组方向中每组取其中一个方向上的子单元，共计四个方向上的四个子单元为无网格单元；

5）对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中两组方向上的子单元都取到，另一组方向上取其中一个方向上的子单元，共计五个方向上的五个子单元为无网格单元。

3.根据权利要求1所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，其特征是：将质点弹簧区域内部符合转换条件的空子单元转化为质点弹簧单元，其转换条件为：

1）上下、左右、前后三组方向中存在两组及两组以上方向的子单元为质点弹簧单元；

2）上下、左右、前后六个方向上的六个子单元都为空子单元；

3）对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，每组方向取一个子单元，共三个方向上的三个子单元为质点弹簧单元；

4）对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中一组方向上的子单元都取到，另两组方向中每组取其中一个方向上的子单元，共计四个方向上的四个子单元为质点弹簧单元；

5）对于边界上的单元，上下、左右、前后三组方向上的每组子单元不是全部存在的情形，其中两组方向上的子单元都取到，另一组方向上取其中一个方向上的子单元，共计五个方向上的五个子单元为质点弹簧单元。

4.根据权利要求1所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，其特征是：剩余的空子单元作为过渡单元，是指： 

1）不符合无网格单元转换条件，也不符合质点弹簧单元转换条件的空子单元；

2）上下、左右、前后六个方向上的六个子单元同时包含无网格单元和质点弹簧单元，或者同时包含无网格单元、质点弹簧单元和空子单元。

5.根据权利要求1～4任一项所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，其特征是：在前处理过程中，假定软组织为均匀、各向同性材料，具有准不可压缩性、线性粘弹性特征，并确定初始化数据，定义初始条件、位移边界条件和应力边界条件；然后用弹簧模拟粘弹性体的弹性，阻尼器模拟其粘性，并假设弹簧的变形与负载成比例，而阻尼器的变形与变形的速度成比例；其次各向同性材料的变形可以分离成体积改变和等体积的形状畸变两部分，故把应力和应变张量分解成它的球形张量和偏斜张量，进而构造Kelvin模型的本构方程，模拟软组织应力松弛和蠕变特性；最后建立几何方程、平衡方程，并施加边界条件和初始条件，为软组织建立线性粘弹性生物力学模型。

6.根据权利要求5所述的基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法，其特征是：在计算形变过程中，首先生成初始的背景网格，根据外部载荷划分无网格区域和质点弹簧区域，建立从无网格到质点弹簧数据之间的数据结构，以降低自适应过程中数据重用现象发生的频率，其次在连接区域建立过渡单元，在过渡单元内建立过渡节点及过渡节点的近似位移函数，实现无网格区域与质点弹簧区域之间的光滑过渡，然后施加边界条件，结合Kelvin标准线性粘弹性力学模型，采用增量法求解位移，根据位移求应力、应变。

Images