基于MLS-MWS混合模型的软组织建模

陈卫东 刘 波 朱奇光*

1(燕山大学信息科学与工程学院，河北 秦皇岛 066004)2(河北省特种光纤与光纤传感重点实验室，河北 秦皇岛 066004)

基于MLS-MWS混合模型的软组织建模

陈卫东1,2刘 波1朱奇光1,2*

1(燕山大学信息科学与工程学院，河北 秦皇岛 066004)2(河北省特种光纤与光纤传感重点实验室，河北 秦皇岛 066004)

针对无网格模型计算量大和边界处理繁琐的问题，提出一种基于改进最小二乘近似函数的无网格强弱势混合模型(MLS-MWS)。首先，基于场理论将模型划分为内部域和外部域，在内部域中，不与自然边界相连接的节点通过无网格强势法来构建平衡方程，该方法无数值积分，从而达到计算量的减少；在外部域中，通过无网格局部弱势法来弱化边界问题。其次，在强势方法中，需利用最小二乘近似式表示位移的二阶导数，针对其计算量大和效率低的问题，提出改进移动最小二乘近似函数和导数，通过将节点值进行两次MLS导数插值，求出该节点的二阶导数，进一步简化计算；最后采用PHANTOM触觉交互设备，对肝脏模型进行拉伸与挤压的形变模拟实验。 实验结果表明，MLS-MWS混合模型相比无网格伽辽金模型，单步执行时间从17.3 ms减少到13.8 ms，频率从56.80帧/s提高到72.46帧/s。该混合模型结合无网格强弱势法和双重无网格配点法，达到优势互补，有效提高模型的实时性。

软组织建模；无网格法；双重配点无网格法；MLS-MWS模型

引言

虚拟手术仿真系统在医学教学与实践中占据重要地位，通过该仿真系统，医生可以进行大量实践来丰富自己的临床经验[1]。软组织形变模型的构建是虚拟手术仿真系统的核心步骤，软组织模型运用计算机图形学与处理、可视化技术和生物医学等相关学科知识，获得可视化图像和力学反馈数据。通过对所获得数据的处理，可以模拟出软组织模型的组织形态特性与力学形态特性，因此软组织建模研究具有重要的理论依据与实际意义。

软组织模型通常分为弹簧质点模型和无网格模型。Vollinger等提出了改进的质点弹簧模型[2]，模型参数的选择是该模型的难题，在某些特殊的仿真环境中，系统会因参数选择不当而引起震荡的不稳定现象。闫雒恒等提出了一种改进四面体网格的弹簧质点模型[3]，通过从新定义弹性力学方程来描述形变过程中应力的变化，提高了仿真的逼真度和实时性，但该模型还存在稳定性差的问题。陈卫东等提出了一种简单的弹簧质点混合模型[4]，通过选择有效的节点，计算各层形变区域，相比传统的弹簧质点模型提高了逼真度，但稳定性相对较差。宋爱国等提出了一种高效的层状菱形链接连接模型[5]，通过搭建双通道力/触觉交互的虚拟肺仿真系统，提高了逼真度和计算速度，但是该模型的稳定性相对较差。Francois等提出了运用超弹性质点连接来快速计算软组织的形变与实时仿真[6]，其仿真的逼真度比较高，但模型的稳定性相对较差，容易变形。Zhang等采用全拉格朗日近似的无网格伽辽金方法建立软组织模型[7]，通过使用大量节点来建立形函数，在大形变中可以显著延迟网格畸变，提高了仿真的逼真度，但该方法计算过程复杂。Steinemann等提出了一种基于无网格和有限元方法的混合模型[8]，但该方法计算复杂，计算效率低。刘雪梅等提出了无网格弱势法与质点弹簧模型[9]，将形变区域分为手术区域和其他区域，手术区域采用无网格弱势法，具有较好的精度与稳定性，但该模型的缺点是计算量大。许少剑等提出了一种最小二乘法的无网格方法对模型进行研究[10]，减少了畸形网格生成对虚拟手术仿真系统所产生的不稳定现象，无需网格，减少了计算量，但是如何处理边界问题没有得到解决。

本研究针对无网格模型计算量大和边界处理繁琐的问题，提出了一种基于改进最小二乘近似函数的无网格强弱势混合模型。首先，基于场理论将模型划分为内部域和外部域，在内部域中所有节点通过无网格强势方法建立平衡方程，减小计算量，在外部域中利用无网格弱势方法弱化边界问题；其次，在强势方法中，利用双重无网格配点法改进近似函数二阶导数的求导过程，进一步简化计算。该混合模型继承了无网格强弱势法和双重无网格配点法的优点，从而有效提高了模型的实时性。

1 模型构建方法

在建模过程中，无网格算法只需对场内节点建立离散方程，从而构建出模型的形状，但是传统的无网格算法存在计算量大和边界问题，因此本研究提出了一种改进最小二乘近似函数的无网格强弱势法，包括内部域与外部域的划分和利用改进的最小二乘函数建立无网格强弱势离散方程。MLS-MWS方法的核心思想是：基于场理论知识，将场内任意节点分别通过无网格强势和无网格局部弱势来建立离散方程。

改进最小二乘近似函数的无网格强弱势方法具有相同的离散形式，因此在整体离散域Ω上，MLS-MWS的离散方程为

(1)

式中，n为求解域的节点数，u为形函数的函数值，K为系统的刚度矩阵，F载荷向量。

如图1所示，求解域Ω和边界Γ由一些任意分布的节点组成，Ωq是积分域的任意节点。如果Ωq不与自然边界相连接，则采用无网格强势方法；反之，采用无网格局部弱势方法[11]。将求解域划分为内部求解域Ωin和边界域ΩB，有

(2)

式中,CI表示以场点为场域圆心、半径为κdIS的圆，dIS为节点支持域的半径，κ为改变圆形场域的控制参数。

利用调整控制参数κ，通过改变内部域Ωin与边界域ΩB的节点数，可以有效改变内部域与外部域的大小。

图1 MLS-MWS方法示意图Fig.1 A schematic diagram of MLS-MWS method

在软组织建模时，时间的快慢将严重影响虚拟仿真系统整体性能。Γu、Γt分别表示位移边界和面力边界，考虑强势边界真值问题时，其控制方程表示如下：

在Ω中，平衡方程为

(3)

在Γu和Γt上，边界条件分别为

(4)

(5)

1.1 无网格内部域的算法

在内部域Ωin中，数值积分和边界问题不做考虑，只需要满足平衡方程。平衡方程之间采用强势方法离散，使得内部域Ωin中的节点精确满足平衡方程，位移μ的二阶导数利用最小二乘法近似式表示，即可得到Ωin中节点的离散方程，有

(6)

在z=0时，KIJ,in、uJ和FI,in分别为

(8)

式中，KIJ,in、FI，in分别为系统刚度矩阵[12]与载荷向量[12]，E、μ分别表示弹性模量和泊松比，φJ,ij是通过式(33)得出的。

1.2 无网格边界域的算法

尽管无网格局部弱势法已经向理想无网格法的发展迈出了重要一步，尤其在复杂形状的边界问题上，数值积分仍然是一项繁琐的任务。因此，应该尽量减少使用数值积分。利用边界区域场点少，没有切割和缝合等复杂问题。根据无网格局部弱势法，可以弱化边界问题。在边界域中，根据弹性力学中的本构关系、几何方程以及应力-面力关系，可以求得应力δij和面力ti的移动最小二乘近似式，可以得到ΩB中节点的离散方程，有

(9)

在z=0时，KIJ,B和FI,B[12]分别为

(10)

(11)

其他相关矩阵表达式如下：

(12)

(13)

(14)

(15)

在边界域ΩB中，子域通过加权残值法，将本质边界条件通过罚函数引入到残值方程，得到

(16)

式中：α为罚因子；Ωq为局部积分域；Γqu为局部积分域边界∂Ωq与位移边界Γu相交的部分；γi为权函数，任意形状的局部积分域均可使用，本研究的局部积分域采用圆形域。

在局部积分域中，边界Γq可以包括三部分：内部边界Γqi、局部积分域边界∂Ωq与位移边界Γu相交的部分Γqu、自然边界Γqt。文献[13]对式(16)利用散度定理进行处理，可以得到边界域ΩB中控制方程的局部弱形式，有

(17)

1.3 移动最小二乘近似函数和导数

假设待求函数u(x)在求解域Ω中N个节点XI(I=1,2,3，…，n)的函数值已知，即uI=u(xI)，则u(x)局部近似函数为

(18)

二维空间单项式基函数中的线性基和二次基分别为

(19)

(20)

三维空间单项式基函数中的线性基和二次基分别为

(21)

(22)

(23)

令J取最小值，即

(24)

通过式(24)得

(25)

(26)

式中

(27)

由式(26),可得待定系数向量α(x),即

(28)

将式(28)代入式(18)，得

(29)

其中，形函数

(30)

在配点无网格法中，利用近似函数的二阶导数来表示位移的二阶导数，由于其计算两次形函数导数，使得计算量大、效率低。本研究利用双重无网格配点法[14]来进行计算。首先引入两组点，即节点xI和计算点xke。

(31)

(32)

式中

(33)

利用改进的最小二乘近似函数求二阶导数，降低了计算的复杂度，避免出现病态方程，使方程的稳定性增强，且通过对节点进行两次MLS导数插值，使结果具有光滑的连续性，在软组织形变过程中更能逼真地进行形变。

2 模型实验

以肝脏模型为例，其泊松比为0.35。采用SensAble科技公司的PHANTOM触觉交互设备，基于酷睿i5处理器，4GB内存，NVIDIAGeForceGTX960M显卡台式电脑。在Windows7操作系统下，采用VisualC++ 6.0 开发环境，并结合OpenGL三维图形标准，搭建了虚拟软组织实时力反馈交互系统。

在实验中，分别用无网格伽辽金方法和MLS-MWS方法构建肝脏模型，并在相同场点施加相同载荷的拉力和压力，然后计算出各自的位移大小，进行对比实验。同时，通过设置参数κ(无明确的物理意义)的大小来控制内部域与边界域大小。由于MLS-MWS法的计算效率随着边界域节点数占总节点数的比例值减小而提高，参数κ通过几组数据测试来进行设置，当0.5≤κ≤2.0时，使得MLS-MWS法减少了计算复杂度，提高了实时性。

在实验中，无网格伽辽金模型和MLS-MWS模型在质点个数为293时，施加力F=0.3N，分别对无网格伽辽金模型和MLS-MWS模型进行平行对比实验。

3 形变仿真结果

在虚拟手术仿真系统中，逼真度与实时性是决定仿真成功的指标。视觉逼真度是描述模型仿真的一个能力，描述能力越强，模型的逼真度越高。图2为肝脏模型的拉伸仿真，图3是肝脏模型压迫仿真图。

图2 拉力作用下的受力Fig.2 The liver rendering under action of pull force

图3 压力作用下的受力Fig.3 The liver rendering under action of pressure force

在实验中，分别运用无网格伽辽金法[9]和所提出的MLS-MWS方法构造肝脏模型，并在相同场点施加相同载荷的拉力和压力，然后计算出各自的位移大小，进行对比实验。图4、5分别是拉力与压力作用下的位移曲线。通过仿真实验结果可知，在施加相同拉力或压力下，MLS-MWS混合模型与无网格模型的位移误差为±0.12。从数值上可以看出，该混合模型的精度与无网格模型的精度基本相同，可以说明MLS-MWS混合模型有效。

图4 拉力作用下的位移曲线Fig.4 The displacement curve under action of pull force

图5 压力作用下的位移曲线Fig.5 The displacement curve under action of pressure force

为了验证MLS-MWS法的计算效率，在实验中，当κ=0.8、1.0、1.5、2.0时，分别对无网格伽辽金方法和MLS-MWS方法进行了平行对比实验，只有当κ=1.0时效果最好。在实验中，当κ=1.0时，内部域节点数为232，边界域节点数为61，分别施加不同载荷。由表1的实验数据可以得出，MLS-MWS混合模型的单步执行时间比无网格模型的单步执行时间少了3.8ms，频率提高了15.64帧。结果说明，MLS-MWS混合模型的计算速度快，减少了计算时间，实时性优越。同时，在不影响精度的前提下，MLS-MWS混合模型单步执行时间比无网格模型的单步执行时间约提高了20%，频率提高近15.64帧，进一步说明该混合模型实时性和逼真度优于无网格模型。

表1 无网格伽辽金法与MLS-MWS法计算效率对比

4 讨论和结论

虚拟手术是虚拟现实技术在医学方面上的一个重要应用，它利用医学图像数据和计算机图形学构建出软组织模型，模拟出虚拟的手术环境，并利用触觉交互设备与之进行交互。软组织形变模型的构建是虚拟手术仿真系统的核心步骤，由于软组织模型在生物黏弹性、各向异性等方面与软组织器官有一定差距，因此构建一个实时性与逼真度高的软组织模型至关重要。随着虚拟手术系统被引入临床教学中，医生通过大量的训练并结合临床经验，减少了手术的风险系数。因此，通过虚拟手术仿真系统可以减少手术的风险，降低培养一个优秀外科医生的费用。同时，通过典型病例并利用该仿真系统对医生进行培训，对提高医生的操作能力、质量和医学水平有着深远的意义。

本研究提出MLS-MWS混合模型，通过调节参数κ的大小来控制内部域与外部域的范围。在内部域中，采用无数值积分的无网格强势方法，减少了计算的复杂程度；在外部域中，采用无网格局部弱势方法，使形状函数具有Kroneckerdelta函数性质，可以有效处理边界问题，使计算结果准确稳定。同时，在强势方法中，利用双重无网格配点法改进近似函数二阶导数的求导过程，进一步简化了计算。实验结果表明，MLS-MWS混合模型单步执行时间比无网格模型的单步执行时间约提高了20%，频率提高了近15.642帧，极大地提升了医疗手术过程中软组织的实时形态及其逼真程度，降低了手术风险，提高了手术成功几率。

本研究提出了一种MLS-MWS混合模型，减少了建模过程中的计算量，弱化了边界问题，大幅地提高了模型的实时性和逼真程度。但是，在无网格软组织建模中，边界域节点数占总节点数的比重值越小，软组织模型的计算效率越高。本研究提出MLS-MWS混合模型，通过控制参数κ的大小来控制边界域的节点个数，但是κ值的选取需进行大量实验，从而导致该方法选取的κ值存在一定误差，因此如何有效控制内部域与外部域的范围是下一步研究的重点工作。

[1] 吴峥, 谢叻, 马浩博. 虚拟手术实时物体碰撞检测和软组织变形研究[J]. 计算机仿真, 2010, 27(2):255-259.

[2] Völlinger U, Stier H, Priesnitz J, et al. Evolutionary optimization of mass-spring models[J]. Cirp Journal of Manufacturing Science & Technology, 2009, 1(3):137-141.

[3] 闫雒恒. 基于改进弹簧振子模型的软组织形变仿真[J]. 计算机仿真, 2012(7):330-333.

[4] 陈卫东, 陈攀攀, 朱奇光. 基于改进弹簧-质点模型的形变建模及力反馈算法研究[J]. 生物医学工程学杂志, 2015，32(5)：989-996.

[5] 张小瑞, 孙伟, 宋爱国,等. 双通道力/触觉交互的虚拟肺手术仿真系统[J]. 仪器仪表学报, 2012, 33(2):421-428.

[6] Goulette F, Chen Zhuowei. Fast computation of soft tissue deformations in real-time simulation with Hyper-Elastic Mass Links[J]. Computer Methods in Applied Mechanics & Engineering, 2015, 295:18-38.

[7] Zhang GY, Wittek A, Joldes GR, et al. A three-dimensional nonlinear meshfree algorithm for simulating mechanical responses of soft tissue[J]. Engineering Analysis with Boundary Elements, 2014, 42(42):60-66.

[8] Steinemann D, Otaduy MA, Gross M. Splitting meshless deforming objects with explicit surface tracking[J]. Graphical Models, 2009, 71(6):209-220.

[9] 刘雪梅, 毛磊, 李运华,等. 耦合无网格迦辽金与质点弹簧实现软组织形变仿真[J]. 计算机辅助设计与图形学学报, 2013, 25(1):1-6.

[10] 许少剑, 陈国栋. 基于移动最小二乘法的无网格肝脏形变研究[J]. 计算机仿真, 2015, 32(4):408-413.

[11] 肖毅华, 胡德安, 韩旭. 弹性静力问题的无网格弱-强形式结合法[J]. 计算力学学报, 2010, 27(5):764-769.

[12] Gu YTGR. A meshfree weak-strong (MWS) form method for time dependent problems[J]. Computational Mechanics, 2005, 35(2):134-145.

[13] 龙述尧, 姜琛, 郑娟. 三维弹性静力问题的无网格局部Petrov-Galerkin法[J]. 湖南大学学报(自然科学版),2013, 40(12):55-61.

[14] 张雄, 刘岩, 马上. 无网格法的理论及应用[J]. 力学进展, 2009, 39(1):1-36.

Soft Tissue Modeling Based on MLS-MWS Mixing Method

Chen Weidong1,2Liu Bo1Zhu Qiguang1,2*

1(InsituteofInformationScienceandEngineering,YanshanUniversity,Qinhuangdao066004,China)2(TheKeyLaboratoryforSpecialFiberandFiberSensorofHebeiProvince,Qinhuangdao066004,China)

Targeting to problem of the large computation amount of mesh-free model and tedious boundary processing, MLS-MWS based on the improved least square approximation was proposed in this work. At first, the model was divided into the internal domain and external domain based on the field theory. In the internal domain, nodes far from the natural boundary constructed the equation of equilibrium with the mesh-free method. The method has no numerical integration, which may reduce the computation amount. In the external domain, the boundary problem was weakened through the mesh-free partial weak method. Secondly in the strong method, second-order derivative of the displacement should be indicated with the least square approximation equation. Targeting tothe large computation amount of mesh-free model and low efficiency, it was proposed to improve the moving least square approximation function and derivate, and MLS derivative interpolation was carried out for the node value, for working out the second-order derivative and further simplifying the computation. Eventually, deformation simulation experiment of stretching and squeezing the liver model was carried out with the PHANTOM touch interactive equipment. According to the results, the single-step execution time of MLS-MWS combined model was reduced to 13.8 ms from the original 17.3 ms of the mesh-free model, while the frequency was improved from 56.80 frames/s to 72.46 frames/s. The mixed model integrated advantages of the mesh-free weak-strong method and dual mesh-free point-allotting method, which improved the instantaneity of the model effectively.

soft tissue modeling; meshless method; double match point meshless method; MLS-MWS hybrid model

10.3969/j.issn.0258-8021. 2016. 06.009

2016-01-12， 录用日期:2016-05-19

国家自然科学基金(61201112)；河北省自然科学基金(F2016203245)；河北省普通高等学校青年拔尖人才计划(BJ2014056)

R318

A

0258-8021(2016) 06-0699-06

*通信作者(Corresponding author)， E-mail: zhu7880@ysu.edu.cn