陈国栋，党琪琪，叶东文

福州大学 物理与信息工程学院，福建福州 350000

基于SPH及形状约束的血流实时仿真

陈国栋，党琪琪，叶东文

福州大学 物理与信息工程学院，福建福州 350000

针对采用光滑粒子流体动力学（SPH）方法仿真血液流动时，由于求解Navier-Stokes方程粘滞力项的运算量大而很难实现实时性仿真的问题，本文提出了一种基于SPH及形状约束的血流实时仿真方法：首先确定血流仿真模型并将整个血流仿真系统离散为一个粒子系统，其次通过SPH方法求解血流控制方程求得血液的流体速度，然后在SPH方法模拟流体的基础上通过形状约束算法控制血粒子的运动规律求得血液的固体速度，最后对SPH方法求得的流体速度和形状约束算法求得的固体速度进行线性插值得到血粒子的实际速度，进而更新血粒子的位置信息并进行可视化处理，最终实现血液流动现象的实时性仿真。实验结果表明，该方法通过形状约束算法可以模拟出血液的粘性效果，避免了传统方法中的复杂计算过程，能够快速地模拟出血流效果，满足实时性要求。

光滑粒子流体动力学；粘滞力；形状约束；血流仿真；实时性

0 前言

随着计算机图形处理技术的进步，虚拟手术仿真系统成为外科医生临床实验和医学院学生接受教育的一个重要平台，医生在进行相应的真人操作之前在虚拟环境中进行反复的训练可提高操作的熟练程度，进而提高手术成功的几率。虚拟手术仿真系统是一个应用前景广阔的平台，它可以在避免病人参与的前提下为学生提供适当的培训环境，提高医学院学生的学习效率。

在手术中，外科医生经常会使用一些工具来执行切口及剥离组织，很可能会导致出血，因此血流的模拟是虚拟手术仿真系统的重要组成部分。血液作为粘弹性流体所具有的特殊物理性质使得血流的实时性模拟成为一大挑战性难题。流体运动规律由一组被称为Navier-Stokes（N-S）方程组的偏微分方程来控制，这些方程的积分求解随着时间的推移在计算上是非常复杂的，实时模拟是很困难的，很多计算机科学家已做出许多尝试来加快其求解速率。

1 相关工作

目前，医学仿真领域的许多国内外学者对于血流的实时仿真已做了大量的工作。Deschamps等[1]使用嵌入式边界法对血管中的血流进行了模拟，取得了一定的仿真效果，但实时性和真实感都有待提高。Cebral等[2]应用计算流体动力学对动脉血管中的血液流动模型进行了研究。Liu等[3]基于粒子系统模拟了动脉血流喷出以及血液在空中自由飞溅的效果，其算法的主要特点是视觉上的逼真性，但在一定程度上忽略了血液精确的物理特性，如果将血液的复杂粘性考虑在内，则很难达到实时性要求。杨礼波等[4]在GPU的CUDA运算平台上采用光滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH)）方法对N-S方程进行求解[5]，实现了血液流动的仿真效果且实时性有了一定的提高，但是N-S方程中粘滞力项的复杂求解过程使得其实时性的提升受到了一定程度的限制。施鹏等[6]通过对传统拉格朗日粒子法模拟血流的模型进行简化实现了动态模拟渗血的实时性要求，但真实感有待进一步提升。

无网格法是一种多功能的仿真方法，既可用于弹性固体的仿真，也可用于粘性流体的仿真。Gerszewski等[7]将无网格法用于仿真弹塑性固体，通过添加弹性力项而实现了弹塑性固体的动画效果。Chang等[8]使用基于粒子的方法对高粘弹性流体进行了实时仿真，并取得了较好的仿真效果。形状约束算法是一种基于无网格的方法，可试图通过恢复物体的原有形状而使物体具有粘性特性。

血液本身具有流体的属性，又因为血液粘性的存在使其也具有了类似弹性固体的属性，因此血液属于一种粘弹性流体。基于以上背景，本文提出了一种简单、快速的仿真血液流动的方法：首先使用SPH方法求解血液流体控制方程获得血粒子的流体速度，其次在SPH方法模拟血液流动的基础上通过形状约束算法来控制血粒子的运动规律从而获得血粒子的固体速度，然后通过插值系数将两种算法的速度进行耦合得到血粒子的真实速度，最后更新血粒子的位置信息并进行可视化处理。本研究首先确定血流仿真模型并进行SPH方法求解，然后将形状约束算法用于血流仿真系统求解血粒子的实际速度，最后进行实验来验证该算法的可行性。可视化部分采用经典的Marching Cubes等值面构造算法。

2 血液流体与固体速度

2.1 血液流体速度

正常生理条件下，血液是不可压缩的，且流动缓慢，因此将其作为不可压牛顿流体来进行模拟。其控制方程组如下所示：

对于压力项，核函数取如下形式：

对于粘滞力项，本研究使用形状约束算法来实现血粒子之间的粘性效果，在此不予考虑。

对于重力项，核函数取如下形式：

其中h为光滑长度。

由牛顿第二定律可知：

由式(5)可以求得血粒子的加速度，进而计算出血粒子的流体速度：

2.2 血液固体速度

2.2.1 形状匹配思想

形状匹配的整体思想[9]相对比较简单，其概述图形见图1。将形状匹配的基本思想转换为数学计算，阐述为两个平移和一个旋转：首先平移原始对象到坐标零点，然后对其进行旋转，最后再平移它到目标点（图2）。

图1 粒子形状匹配思想概述图形

图2 形状匹配旋转、平移示意图

2.2.2 形状匹配算法

依据形状匹配思想，将血流模型仿真为一个简单的粒子系统，这是本研究在2.1节中已经完成的任务，然后以2.1节中所取的光滑半径h为单位长度将整个血粒子系统划分为多个相互重叠的微团（区域），形状约束算法能够维持血液微团的原有形状，以实现血液的粘弹性效果（图3），接下来的任务便是求解最佳的平移矢量和旋转矩阵。

图3 形状约束算法示意图

其中，ωi为加权因子，本研究选择ωi=mi。物体的旋转一般都是围绕它的质心进行的，因此最佳的平移矢量为初始的形状质心到实际形状质心的平移矢量。质心计算公式如下：

其中，t0为原始形状的质心，t为实际形状的质心。定义初始的相对位置为实际的相对位置为，接下来通过找到最优的线性变换A来简化最优旋转矩阵R的求解。公式(7)可以转变为依据最小二乘法则有最优变换矩阵A：

2.2.3 固体速度

在得到了旋转矩阵R以及平移矢量t与t0之后，粒子目标位置gi可通过下式求得：

其中，Nr为粒子i所属区域的个数。如图4所示，粒子b既被划分到区域a又被划分到区域c，所以有Nr=2。

图4 粒子作用区域示意图

基于计算出的目标位置（图5），一个试图恢复物体原有形状的额外速度被引入，如下式所示：

其中，fext为粒子所受到的外力（本文仅包含重力）。Ui为粒子i所属区域的集合。

图5 基于目标位置的形状恢复示意图

3 血流速度求解

由2.1节和2.2节可计算出两种算法下血粒子的速度，然后通过线性插值[11]来求得血粒子的实际速度。所谓线性插值即通过参数ω和(1-ω)赋予两个速度值不同的权重且对两个速度值进行线性叠加，如式(12)所示：当ω=0时，速度部分只有一项，即为形状约束算法求得的速度，该速度试图恢复血液的原始形状而使血液保持固体的属性；当1-ω=0（即ω=1）时，速度部分也只有一项，即为不考虑粘滞力情况下SPH算法求得的速度，该速度因为不考虑粘滞性而具有理想水的属性；当ω的取值为0～1时，线性叠加后的速度具有理想水和弹性固体之间的属性，且ω的值越大越接近理想水的属性，反之越接近弹性固体的属性。血液的属性处于理想流体水与弹性固体之间，是一种粘弹性流体，所以该算法对于血流仿真具有可行性。

在得到血粒子速度信息之后，通过欧拉积分方法更新血粒子的位置信息：

由方程组(1)中的动量守恒方程可以看出，在用SPH算法仿真血流时，粘滞力主要取决于血液粘滞系数μ。在本研究的算法中，则是通过调整插值系数ω进行仿真实验来确定血液粘性所对应的近似ω值。

4 算法思路

本研究算法主要包括以下几个步骤：首先是算法的初始化，将整个血流仿真系统初始化为一个粒子系统，并初始化血粒子的参数信息（比如血粒子的临床参数信息，初始位置、速度以及加速度信息等）。其次是使用SPH算法求解血流控制方程进而求解血流流体速度的部分，在该过程中要循环计算每一血粒子的密度、内部压力以及自身重力来求得粒子所受的合力，根据式(5)计算得出血粒子的加速度，再由式(6)计算得出血粒子的流体速度；然后利用形状约束算法求解血流固体速度的部分，在该过程中要循环计算旋转矩阵R和平移矢量t，根据式(10)求解血粒子的目标位置，进而根据式(11)求得血粒子的固体速度。最后根据式(12)得出血粒子的实际速度并更新血粒子的位置信息，使用Marching Cubes算法可视化处理后即可得出粘弹性流体仿真效果。具体步骤如下：

（1）初始化仿真环境：初始化血液参数信息，初始化血粒子的原始位置信息、原始速度以及原始加速度信息。

（2）确定血流的牛顿流体模型，通过SPH方法将连续的流体动力学方程转变为离散形式，并求解离散方程得到血流的流体速度。

（3）利用形状约束算法控制血粒子的运动规律，求解血流的固体速度。

（4）通过线性插值系数ω得到流体的实际速度ui。

（5）根据血粒子的速度信息更新血粒子的位置信息。

（6）使用Marching Cubes算法实现可视化效果。

（7）调整ω进行重复的仿真实验，直到得到类似血液的仿真结果。

整体算法流程，见图6。

图6 整体算法流程图

5 实验结果与分析

5.1 实验环境

本研究实验的硬件环境：CPU双核1.60 GHz，RAM 4 GB，显卡NVIDIA GeForce GT 620M。

本研究实验的软件环境：操作系统为Windows系统（64位），开发环境为Microsoft Visual Studio 2012，算法实现的语言采用C++，实验结果渲染采用OpenGL来实现。

5.2 实验结果与分析

在仿真实验中，SPH算法中支持域取单倍的光滑长度，形状约束算法区域划分长度也取单倍的光滑长度。为计算方便，血粒子质量mi取值为1.0，粒子信息更新的时间步长Δt取为0.02 s。为了较好地进行实时性的对比，实验中所提到的方法均用来仿真皮肤表面切割后的血液流动效果。实验结果中坐标系O-XYZ为世界坐标系，即虚拟的手术环境，XOZ平面为模拟的皮肤表面。皮肤表面切口效果如图7所示；血流仿真结果如图8～11所示。

图7 皮肤表面切口效果示意图

图8 本研究算法血流仿真效果（粒子总数1000，ω=0．8）

图9 传统SPH算法血流仿真效果（粒子总数1000，血流粘滞系数μ=4 mPa · s）

图10 算法血流仿真效果[2]（粒子总数1000，血流粘滞系数μ=4 mPa · s）

图11 本研究算法血流仿真效果（粒子总数2000，ω=0．8）

由实验结果可知，当ω=0.8时，仿真效果最近似为血流仿真效果。对比图8、图9和图10，可以得到如下结论：在相同的条件下（湿度、温度、血流仿真模型取相同时），相同切口大小（即仿真时所采用的粒子总数相同）时，本研究所采用方法的实时性与传统SPH方法相比有很大的提升，与文献[2]中的方法相比仿真效率也有很大的优势，帧率近似为传统SPH方法的3倍、文献[2]中方法的1.5倍。对比图8和图11，可以得到如下结论：在相同的条件下（湿度、温度、血流仿真模型取相同值），切口大小（即仿真时所采用的粒子总数不相同）不同时，切口越大（即仿真时所采用的粒子总数越多），实时性越低，仿真效率越低，但连续性较强，真实感较好。

本研究算法的实时性主要取决于血粒子目标位置的计算，尤其是当血粒子数目比较多的时候，血粒子位置计算所消耗的时间也就相对较多，但较近些年的血流仿真方法在效率上还是有很大提升的。表1是两种粒子数情况下3种算法的计算效率（渲染部分均采用Marching Cubes算法）。结果表明，本研究算法能够高效、快速地仿真皮肤表面的血液流动，降低了计算过程的复杂度，在保证真实感的前提下提高了血流模拟的实时性。

表1 两种粒子数情况下3种算法的计算效率

6 结论

流血现象广泛存在于我们的日常生活中，比如意外划伤出血、手术中皮肤表面切割出血、动脉出血、器官表面渗血等，因此对于血流的仿真就成为虚拟现实领域的一大热点。又因为对实时性和真实感的要求较高，血流仿真也是医学仿真领域的一大难点。本研究在SPH算法模拟流体的基础上加入形状约束算法来控制血粒子的运动规律，并通过调节插值系数来实现血流的实时仿真，可视化部分采用经典的Marching Cubes算法，最终在保证真实感的前提下实现了血流现象的实时仿真。未来工作可以考虑使用GPU对本研究算法进行加速，也可使用本研究算法模拟更复杂的场景。

[1] Deschamps T,Schwartz P．Vessel segmentation and blood flow simulation using level-sets and embedded boundary methods[C]．International Congress Series,Nederland:2004．

[2] Cebral JR,Lohner R．Efficient Simulation of blood flow past complex endovascular devices using an adaptive embedding technique[J]．IEEE Trans Med Imaging,2005,24(4):468-476．

[3] Liu XM,Hen CY．Bleeding simulation based particle system for surgical simulator[C]．Pacific-Asia Conference on Knowledge Engineering and Software Engineering,Shenzhen:2009．

[4] 杨礼波．虚拟手术系统中的流血模拟[D]．郑州:华北水利水电学院,2011．

[5] 温婵娟,欧嘉蔚,贾金原．GPU通用计算平台上的SPH流体模拟[J]．计算机辅助设计与图形学学报,2010,22(3):407-411．

[6] 施鹏,熊岳山,徐凯,等．虚拟肝脏手术中实时动态渗血效果模拟[J]．计算机应用,2013,33(10):2911-2913．

[7] Gerszewski D,Bhattacharya H,Bargteil AW．Point-based method for animating elastoplastic solids[C]．Association for Computing Machinery,New Orleans:2009．

[8] Chang YZ,Bao K,Zhu J,et al．High viscosity fluid simulation using particle-based method[C]．International Symposium on Virtual Reality Innovation,Singapore,2011．

[9] Muller M,Heidelberger B,Teschner M,et al．Meshless deformations based on shape matching[J]．ACM Trans Graph,2005,24(3):471-478．

[10] 万旺根,林继承,余小清,等．基于粒子系统和形状匹配的实时无网格变形仿真[J]．计算机应用,2008,28(12):3007-3009．

[11] 唐勇,陈静,吕梦雅．基于SPH及形状约束的粘弹性流体的实时模拟[J]．小型微型计算机系统,2013,34(11):2626-2629．

Real-Time Simulation of Blood Flow Based on SPH and Shape Constrain

CHEN Guo-dong, DANG Qi-qi, YE Dong-wen
School of Physics and Information Engineering, Fuzhou University, Fuzhou Fujian 350000, China

When Smoothed Particle Hydrodynamics (SPH) was deployed to simulate blood flow, the large amount of computation in solving the viscoelastic item of Navier-Stokes (N-S) equations made it dif fi cult to achieve real-time simulation of blood fl ow. In view of the problem, this paper proposed a realtime blood fl ow simulation method based on SPH and shape constrain. Firstly, the blood fl ow simulation model should be determined and the entire blood system was treated as a discrete particle system. Then, the fluid speed was acquired through solving the blood flow control equations with deployment of the SPH method. On the basis of SPH simulation of fl uid, the shape constraining method was used to control the movement of blood particles so as to obtain the solid speed. Finally, the actual speed of blood particles were obtained by linear interpolation of the fluid speed and solid speed so as to update and visualize the position of blood particles and realize real-time simulation of blood fl ow. The experimental results showed that the algorithm could simulate the viscosity effect of blood with deployment of shape constraining, which avoided the complex computation process of the traditional methods, achieved fast blood fl ow simulation results and met the real-time requirements.

smoothed particle hydrodynamics; viscous force; shape constraining; simulation for blood fl ow; real-time

TP391.41

A

10.3969/j.issn.1674-1633.2015.06.005

1674-1633(2015)06-0023-05

2015-03-01

福建省科技计划重点项目支持（2011H0027）；福建省自然科学基金项目支持（2013J05090）。

陈国栋，副研究员，主要研究方向为计算机图形学和虚拟现实技术。

通讯作者邮箱：1285056478@qq．com