虚拟手术中软组织形变建模及力反馈算法研究

陈卫东 赵成龙，3 朱奇光* 关永贞1(燕山大学信息科学与工程学院，秦皇岛 066004)

2(河北省特种光纤与光纤传感重点实验室，秦皇岛 066004)

3(海军91821 部队，潮州 515729)

引言

虚拟手术仿真系统，是虚拟现实(virtual reality，VR)技术在医学领域的一个典型应用。在虚拟手术训练中，引入力触觉反馈，可以使医生在训练时不仅能够看到而且还能感觉到手术器官;通过与虚拟模型和环境进行交互，形成对虚拟模型的一个完整的认识，如同操作真实物体一样，这无疑可使训练更真实、准确、可靠［1］。

人体软组织形变技术是虚拟手术仿真系统的关键技术，在虚拟环境力触觉交互过程中，建立基于物理意义的力触觉模型是非常重要的。大多数软组织的物理形变模型都是基于弹性理论研究的。Majumder 等基于有限元素法(FEM)提出了具有物理意义的软组织形变模型，虽然精确，但要求网格的节点数多、计算成本高［2］。Peterlik 等提出非线性有限元模型，提高了实时性，但仍无法满足力触觉反馈上千赫兹的频率要求［3］。Wang 等用边界元法建立模型，虽对模型的边界进行了离散，简化了计算，但在稳定性方面却存在一定的困难［4］。Zhong等提出具有物理化学特性的反应扩散表面模型，难于进行复杂的手术操作［5］。Chen 在2006 年提出了基于弹簧-质点的三角形表面网格拓扑结构模型［6］，模型原理简单，但稳定性有待提高。徐少平等基于无网格的软组织切割模型的研究现状进行了分析［7］，指出其未来发展方向，但稳定性较差。吴涓等提出了一种基于沿径向方向分割为呈同心圆分布的弹簧-质点模型及实时力觉响应算法［8］，但该模型各点计算不一致，不易扩展。弹簧-质点模型原理简单，易建模，能够满足实时性交互操作要求，但精度有待提高。因此，在保证变形模型准确性的同时，提高计算的实时性，是当前虚拟力触觉交互系统亟待解决的首要问题［9］。

传统的基于表面网络拓扑结构的模型，在模拟软组织形变过程中稳定性较差。本研究对虚拟人体软组织表面模型在形变过程中设置了虚拟弹簧，使得表面模型形变过程中有了体的信息，以实现更加真实的形变效果，形变过程中算法简单，能够很好地满足交互系统的真实性、稳定性和实时性的要求。

1 软组织形变模型

1.1 模型原理及结构

从生物力学的角度来看，皮肤、脂肪等软组织可近似认为是各向同性且均质的［10］。但由于软组织在受力状态下，受力面的大小和方向都在实时变化，这使得基于运动学模型的解必然会产生震荡。传统的基于表面网络拓扑结构的模型，在模拟软组织形变过程中稳定性较差，不能很好地表现软组织的生物力学特性。本研究中的模型是基于表面正六边形模型的改进，为了用表面正六边形模型近似地仿真它们所表示的软组织结构，得到一个比较逼真的含有体行为的仿真效果，在模型中为每个质点设计了一个虚拟体弹簧［11］。虚拟体弹簧与受力方向一致，在仿真开始之前它是不存在的，仿真开始之后它才动态地产生，主要用于力反馈的计算和增加系统的稳定性。

虚拟体弹簧的变形基本思想是:施加在软组织表面上的力，通过与用户直接交互的那个表面接触点，传递到周围的与之相连的其他各个质点，这些质点产生的应力作用在其他相邻质点，这样把力向周围传递，并同时发生变形。在系统中，各质点表面弹簧形变量的叠加可等效为物体表面形变，各层质点虚拟体弹簧弹性力的合力可等效为物体表面的接触力。

设有虚拟体弹簧的正六边形结构模型结构如图1 所示。图中以接触点为中心，用辐射状的表面同心圆分层。把整个形变物体分割成一系列均匀的同心环层，半径依次为r，2r，3r，…，虽然每个同心环层中质点数和半径不同，但由于表面正六边形结构的对称性，为了计算方便，可近似地认为其半径相同，且同一环层内质点刚性连接。相邻层间弹簧弹性系数相同。这样第一层有6 个质点，6 个弹簧。每增加一层质点数增加6 个，弹簧数增加12个。第m 层有质点6m 个，弹簧6(2m -1)个。

图1 设有虚拟体弹簧的正六面体结构模型Fig． 1 Regular hexagon structure model with a virtual spring

1.2 模型形变的力反馈算法

当力F 作用于质点n0时，n0在力的作用下离开初始位置发生形变，此时n0处的虚拟体弹簧发生作用。同时向下移动的质点n0通过表面弹簧带动相邻层质点发生位移，发生位移的各质点处虚拟体弹簧也发生作用。当达到平衡时处质点n0运动到n'0。这里，任取受力法平面一个方向进行分析，如图2 所示。

图2 形变法平面某一方向截面Fig．2 The schematic diagram in normal plane

在受力作用点处，根据力平衡关系有

对于第i 层质点，有

式中，fTi表示切向力，fNi表示法向力，KNi表示ni质点处的虚拟弹簧弹性系数。ΔZi表示ni质点的位移量，即虚拟弹簧的形变长度。θi为形变后ni质点与ni-1质点连线与水平方向夹角。

当i = m 时，由式(1)和式(3)得

即

设相邻两层质点的切向和法向形变，如图3所示。

图3 相邻层质点形变位置Fig．3 The schematic diagram of surface spring and virtual spring between adjacent layers

li为第i 层质点法向虚拟弹簧的形变量，r 为表面弹簧的原始长度，Δri表示切向相邻质点表面弹簧的形变量。

根据图2 ～图3 的关系，可知各层质点的法向形变量，即发生形变的穿刺深度ΔZ，有

第i 层网格切向表面弹簧的弹性系数KTi分布是影响接触力与变形计算的关键，这里假设切向表面弹簧的弹性系数KTi均相等。由式(2)得

令εi=为第i 层弹簧切向表面弹簧的形变率，αi= 2i - 1，则有

法向上第i 层质点法向虚拟弹簧的形变量li为

在计算软组织受力变形的模型中，计算量的大小在很大程度上取决于传播深度m 的大小。这里m 是一个边界参数，它表示在虚拟物体上以受力点为中心划分的同心圆层数。根据文献［12］的推导，有

式中:lm/r 代表计算区域边界的形变量，也是边界参数;ε1则代表了第一层切向表面弹簧的形变率。

计算中选定m 后，可认为sinθm+1≈0 。此时，式(5)变为

以上的公式推导是基于交互点n0得出，根据模型的对称性，可以看出模型中各点的计算都是一致的。式(6)和式(13)中穿刺深度与传播深度的大小决定了整个模型形变区域大小和反馈力的大小。

2 虚拟手术力反馈系统仿真

2.1 虚拟软组织实时力反馈交互系统

虚拟软组织实时力反馈交互系统的建立，首先要建立软组织形变的几何模型。由医学影像数据出发，基于OpenGL 建立和渲染软组织的几何模型，然后在三维重构的几何模型中加入了弹簧-质点模型，使软组织能够真实地反映器官的物理形变过程。在虚拟场景中，建立虚拟手术器械模型和软组织模型，导入Matlab 的Simulink 中建立仿真模型。在碰撞检测［13］模块中，虚拟手术器械的交互定位用力反馈设备由触觉笔尖控制。当发生碰撞时，程序根据穿刺深度ΔZ 计算出反馈力的大小，操纵者能够实时体验到作用的力反馈。同时，模型形变模块依据穿刺深度更新接触点的位置，依据传播深度m计算出形变区域，在虚拟场景中更新模型渲染图，提供视觉反馈。虚拟手术仿真实现的具体流程如图4 所示。

图4 虚拟手术仿真系统流程Fig．4 The schematic diagram of force feedback interaction system

2.2 模型参数的选择

软组织实时力反馈交互系统的变形仿真一般要求比较高的实时性，这就对该系统的计算速度和刷新频率有很高的要求。要使视觉、力觉再现具有连续性，视觉显示的软组织形变刷新率应该不小于30 Hz，力觉再现的刷新频率应该在1 kHz 以上［14］。即要求力反馈的计算具有快速性的特点。

触觉交互设备是用 SensAble 科技公司的PHANTOM 设备，它能够准确跟踪人手的三维运动，并将由虚拟模型计算的虚拟力实时地反馈给操作者，在力觉上提供逼真的沉浸感。PHANTOM 设备的最大输出反馈力为3.3 N。

传播深度m 是控制形变区域的一个参数，它的大小不仅影响着计算量的大小，而且影响着计算精度。当形变穿刺深度一定时，传播深度m 越大，精度越高，反馈力越接近真实力。考虑到计算量的问题，实验中可根据具体情况确定传播深度m 的大小。根据式(14)，取虚拟体弹簧弹性系数KN0= 50 N/m，r = 0.002 m 进行模拟仿真，得到反馈力与传播深度m 的关系曲线，如图5 所示。

图5 反馈力F 与传播深度m 的关系Fig．5 The relation between the force feedback and the propagation depth m

由式(13)可知，m 的大小由lm/r 和ε1共同决定，选取m =25。表面弹簧原始长度r 为固定长度。根据式(6)和式(12)，受力时软组织虚拟体弹簧形变量与传播深度m 的关系如图6 所示。

图6 虚拟体弹簧形变量与传播深度m 的关系Fig． 6 The deformation of virtual spring with propagation depth m

当接触点形变率ε1取不同值时，得到不同的形变曲线。图6 中表明:在相同力的作用下，不同组织的形变率ε1不同，软组织的形变也不同。ε1越大，表明组织越柔软，受力变形越大。当然，软组织的柔软度还与虚拟体弹簧弹性系数KN的大小［15］相关，由二者共同决定。

虚拟体弹簧弹性系数KN的选取不仅要考虑模型形变的大小，还要考虑力反馈设备反馈力上限的制约因素。根据式(14)，令输出力F 小于设备最大输出值。

2.3 系统仿真结果

使用SensAble 科技公司的PHANTOM 触觉交互设备，基于AMD Athlon(tm)Ⅱ2.80 GHz 双核处理器，1.75GB 内存，NVIDIA GeForce7025 显卡的台式电脑，在 WindowS XP 操作系统下，采用Visua1C++6.0 开发环境并结合OpenGL 三维图形标准，搭建了虚拟软组织实时力反馈交互系统。

基于该实验系统，实现了虚拟软组织网格(皮肤、脂肪等)在虚拟手术器械(虚拟注射器、虚拟手术刀等)单点交互下的按压变形仿真。当操作者控制触觉交互设备运动时，虚拟手术器械也随之运动。当碰撞检测算法检测到虚拟手术器械虚拟软组织发生碰撞时，在给定接触点处运用本研究建立的形变模型及力反馈算法，提供给操作者以视觉与力觉信息，产生真实的交互。实验结果如图7 所示，图中显示的是软组织受外力在单点按压作用下的形变效果。

图7 实验结果。(a)交互图;(b)形变前;(c)形变后Fig．7 The simulation results． (a)Interaction diagram;(b)Before deformation;(c)After deformation

在实时力反馈模块中，提取出实时交互过程中按压操作的反馈力，如图8 所示。

图8 交互过程中按压操作的实时反馈力Fig．8 The real-time feedback force under the pressure

实验结果显示，在施加外力作用之后，软组织形变模型就像一个实体一样很有体积感，在形变过程中不会发生形变失真现象，撤销外力之后表面变形几乎完全恢复。反馈力的计算简单，一次刷新小于1 ms 也满足刷新频率的要求。在交互过程中，力觉感受平稳，模拟效果逼真，富有弹性。

3 结论

针对实时软组织力触觉交互系统的需要，提出了一种具有体信息的表面正六边形结构模型。与文献［12］相比，该模型中每层的各个质点受力与变形情况近似相同，计算量小，保证了力触觉反馈计算的高刷新频率的要求，模拟效果逼真，各个质点的相对位移和弹性力的叠加等效为软组织的变形和接触力。在PHANTOM 触觉交互设备上，通过虚拟手术工具，对虚拟软组织进行单点的按压变形仿真。结果表明，该模型能够很好地满足交互系统的真实性、稳定性和实时性的要求。

［1］ 谢叻，张艳，张天宇，等. 虚拟手术中的力学变形和力觉感知［J］. 医用生物力学，2006，21(3):241 -245.

［2］ Santanu M，Amit R，Subrata P. Simulation of hip fracture in sideways fall using a 3D finite element model of pelvis-femur-soft tissue complex with simplified representation of whole body［J］.Medical Engineering ＆ Physics，2007，29:1167 -1178.

［3］ Igor P，Mert S，Cagatay B，et al. Real-time visio-haptic interaction with static soft tissue models having geometric and material nonlinearity［J］. Computers ＆ Graphics，2010，34:43-54.

［4］ Wang P，Becker AA，Jones IA，et al. Virtual reality simulation of surgery with haptic feedback based on the boundary element method［J］. Computers and Structures，2007，85:331 -339.

［5］ Zhong Yongmin，Bijan S，Julian S，et al. Soft tissue deformation with reaction-diffusion process for surgery simulation［J］. Journal of Visual Languages ＆ Computing，2012，23:1 -12.

［6］ Chen P，Barner KE，Steiner KV. A displacement driven realtime deformable model for haptic surgery simulation ［C］ //Proceedings of 14th Symposium on Haptics Interfaces for Virtual Environment and Teleoperator Systems. Alexandria: IEEE，2006:499 -505.

［7］ 徐少平，李春泉，江顺亮，等. 基于无网格的软组织切割模型的研究进展［J］. 中国生物医学工程学报，2012，31(1):123-128.

［8］ 吴涓，宋爱国，李建清. 用于虚拟力觉再现的弹簧- 质点力/变形模型研究［J］. 系统仿真学报，2006，18(11):3152 -3156.

［9］ 张小瑞，宋爱国，孙伟，等. 用于软组织变形仿真的层状菱形链连接模型［J］. 计算机辅助设计与图形学学报，2010，22(3):473 -479.

［10］ 王彦臻，熊岳山，徐凯，等. 改进的采用表面网格的弹簧振子模型［J］.计算机辅助设计与图形学学报，2007，19(2):168-171

［11］ 鲍春波，王博亮，刘卓，等. 一种用于软组织变形仿真的动态质点弹簧模型［J］. 系统仿真学报，2006，18(4):847 -851.

［12］ 吴涓，宋爱国，李建清. 一种基于物理意义的快速力反馈形变模型及实时力觉响应算法［J］. 传感技术学报，2005，18(1):90 -94.

［13］ Asma AE，Elsayed EH，Magda BF. Rapid collision detection for deformable objects using inclusion-fields applied to cloth simulation［J］. Journal of Advanced Research，2011，16:25.

［14］ Chang Yeonghwa， Chen Yungte， Chang Chiawen， et al.Development scheme of haptic-based system for interactive deformable simulation［J］. Computer-Aided Design，2010，42:414 -424.

［15］ Suriya N，Cenk C. Determination of elasticity parameters in lumped element (mass-spring)models of deformable objects［J］.Graphical Models，2010，72:61 -73.