陳國棟，黨琪琪，葉東文

福州大學(xué) 物理與信息工程學(xué)院，福建福州 350000

基于SPH及形狀約束的血流實(shí)時(shí)仿真

陳國棟，黨琪琪，葉東文

福州大學(xué) 物理與信息工程學(xué)院，福建福州 350000

針對采用光滑粒子流體動力學(xué)（SPH）方法仿真血液流動時(shí)，由于求解Navier-Stokes方程粘滯力項(xiàng)的運(yùn)算量大而很難實(shí)現(xiàn)實(shí)時(shí)性仿真的問題，本文提出了一種基于SPH及形狀約束的血流實(shí)時(shí)仿真方法：首先確定血流仿真模型并將整個(gè)血流仿真系統(tǒng)離散為一個(gè)粒子系統(tǒng)，其次通過SPH方法求解血流控制方程求得血液的流體速度，然后在SPH方法模擬流體的基礎(chǔ)上通過形狀約束算法控制血粒子的運(yùn)動規(guī)律求得血液的固體速度，最后對SPH方法求得的流體速度和形狀約束算法求得的固體速度進(jìn)行線性插值得到血粒子的實(shí)際速度，進(jìn)而更新血粒子的位置信息并進(jìn)行可視化處理，最終實(shí)現(xiàn)血液流動現(xiàn)象的實(shí)時(shí)性仿真。實(shí)驗(yàn)結(jié)果表明，該方法通過形狀約束算法可以模擬出血液的粘性效果，避免了傳統(tǒng)方法中的復(fù)雜計(jì)算過程，能夠快速地模擬出血流效果，滿足實(shí)時(shí)性要求。

光滑粒子流體動力學(xué)；粘滯力；形狀約束；血流仿真；實(shí)時(shí)性

0 前言

隨著計(jì)算機(jī)圖形處理技術(shù)的進(jìn)步，虛擬手術(shù)仿真系統(tǒng)成為外科醫(yī)生臨床實(shí)驗(yàn)和醫(yī)學(xué)院學(xué)生接受教育的一個(gè)重要平臺，醫(yī)生在進(jìn)行相應(yīng)的真人操作之前在虛擬環(huán)境中進(jìn)行反復(fù)的訓(xùn)練可提高操作的熟練程度，進(jìn)而提高手術(shù)成功的幾率。虛擬手術(shù)仿真系統(tǒng)是一個(gè)應(yīng)用前景廣闊的平臺，它可以在避免病人參與的前提下為學(xué)生提供適當(dāng)?shù)呐嘤?xùn)環(huán)境，提高醫(yī)學(xué)院學(xué)生的學(xué)習(xí)效率。

在手術(shù)中，外科醫(yī)生經(jīng)常會使用一些工具來執(zhí)行切口及剝離組織，很可能會導(dǎo)致出血，因此血流的模擬是虛擬手術(shù)仿真系統(tǒng)的重要組成部分。血液作為粘彈性流體所具有的特殊物理性質(zhì)使得血流的實(shí)時(shí)性模擬成為一大挑戰(zhàn)性難題。流體運(yùn)動規(guī)律由一組被稱為Navier-Stokes（N-S）方程組的偏微分方程來控制，這些方程的積分求解隨著時(shí)間的推移在計(jì)算上是非常復(fù)雜的，實(shí)時(shí)模擬是很困難的，很多計(jì)算機(jī)科學(xué)家已做出許多嘗試來加快其求解速率。

1 相關(guān)工作

目前，醫(yī)學(xué)仿真領(lǐng)域的許多國內(nèi)外學(xué)者對于血流的實(shí)時(shí)仿真已做了大量的工作。Deschamps等[1]使用嵌入式邊界法對血管中的血流進(jìn)行了模擬，取得了一定的仿真效果，但實(shí)時(shí)性和真實(shí)感都有待提高。Cebral等[2]應(yīng)用計(jì)算流體動力學(xué)對動脈血管中的血液流動模型進(jìn)行了研究。Liu等[3]基于粒子系統(tǒng)模擬了動脈血流噴出以及血液在空中自由飛濺的效果，其算法的主要特點(diǎn)是視覺上的逼真性，但在一定程度上忽略了血液精確的物理特性，如果將血液的復(fù)雜粘性考慮在內(nèi)，則很難達(dá)到實(shí)時(shí)性要求。楊禮波等[4]在GPU的CUDA運(yùn)算平臺上采用光滑粒子流體動力學(xué)（Smoothed Particle Hydrodynamics, SPH)）方法對N-S方程進(jìn)行求解[5]，實(shí)現(xiàn)了血液流動的仿真效果且實(shí)時(shí)性有了一定的提高，但是N-S方程中粘滯力項(xiàng)的復(fù)雜求解過程使得其實(shí)時(shí)性的提升受到了一定程度的限制。施鵬等[6]通過對傳統(tǒng)拉格朗日粒子法模擬血流的模型進(jìn)行簡化實(shí)現(xiàn)了動態(tài)模擬滲血的實(shí)時(shí)性要求，但真實(shí)感有待進(jìn)一步提升。

無網(wǎng)格法是一種多功能的仿真方法，既可用于彈性固體的仿真，也可用于粘性流體的仿真。Gerszewski等[7]將無網(wǎng)格法用于仿真彈塑性固體，通過添加彈性力項(xiàng)而實(shí)現(xiàn)了彈塑性固體的動畫效果。Chang等[8]使用基于粒子的方法對高粘彈性流體進(jìn)行了實(shí)時(shí)仿真，并取得了較好的仿真效果。形狀約束算法是一種基于無網(wǎng)格的方法，可試圖通過恢復(fù)物體的原有形狀而使物體具有粘性特性。

血液本身具有流體的屬性，又因?yàn)檠赫承缘拇嬖谑蛊湟簿哂辛祟愃茝椥怨腆w的屬性，因此血液屬于一種粘彈性流體?；谝陨媳尘?，本文提出了一種簡單、快速的仿真血液流動的方法：首先使用SPH方法求解血液流體控制方程獲得血粒子的流體速度，其次在SPH方法模擬血液流動的基礎(chǔ)上通過形狀約束算法來控制血粒子的運(yùn)動規(guī)律從而獲得血粒子的固體速度，然后通過插值系數(shù)將兩種算法的速度進(jìn)行耦合得到血粒子的真實(shí)速度，最后更新血粒子的位置信息并進(jìn)行可視化處理。本研究首先確定血流仿真模型并進(jìn)行SPH方法求解，然后將形狀約束算法用于血流仿真系統(tǒng)求解血粒子的實(shí)際速度，最后進(jìn)行實(shí)驗(yàn)來驗(yàn)證該算法的可行性?？梢暬糠植捎媒?jīng)典的Marching Cubes等值面構(gòu)造算法。

2 血液流體與固體速度

2.1 血液流體速度

正常生理?xiàng)l件下，血液是不可壓縮的，且流動緩慢，因此將其作為不可壓牛頓流體來進(jìn)行模擬。其控制方程組如下所示：

對于壓力項(xiàng)，核函數(shù)取如下形式：

對于粘滯力項(xiàng)，本研究使用形狀約束算法來實(shí)現(xiàn)血粒子之間的粘性效果，在此不予考慮。

對于重力項(xiàng)，核函數(shù)取如下形式：

其中h為光滑長度。

由牛頓第二定律可知：

由式(5)可以求得血粒子的加速度，進(jìn)而計(jì)算出血粒子的流體速度：

2.2 血液固體速度

2.2.1 形狀匹配思想

形狀匹配的整體思想[9]相對比較簡單，其概述圖形見圖1。將形狀匹配的基本思想轉(zhuǎn)換為數(shù)學(xué)計(jì)算，闡述為兩個(gè)平移和一個(gè)旋轉(zhuǎn)：首先平移原始對象到坐標(biāo)零點(diǎn)，然后對其進(jìn)行旋轉(zhuǎn)，最后再平移它到目標(biāo)點(diǎn)（圖2）。

圖1 粒子形狀匹配思想概述圖形

圖2 形狀匹配旋轉(zhuǎn)、平移示意圖

2.2.2 形狀匹配算法

依據(jù)形狀匹配思想，將血流模型仿真為一個(gè)簡單的粒子系統(tǒng)，這是本研究在2.1節(jié)中已經(jīng)完成的任務(wù)，然后以2.1節(jié)中所取的光滑半徑h為單位長度將整個(gè)血粒子系統(tǒng)劃分為多個(gè)相互重疊的微團(tuán)（區(qū)域），形狀約束算法能夠維持血液微團(tuán)的原有形狀，以實(shí)現(xiàn)血液的粘彈性效果（圖3），接下來的任務(wù)便是求解最佳的平移矢量和旋轉(zhuǎn)矩陣。

圖3 形狀約束算法示意圖

其中，ωi為加權(quán)因子，本研究選擇ωi=mi。物體的旋轉(zhuǎn)一般都是圍繞它的質(zhì)心進(jìn)行的，因此最佳的平移矢量為初始的形狀質(zhì)心到實(shí)際形狀質(zhì)心的平移矢量。質(zhì)心計(jì)算公式如下：

其中，t0為原始形狀的質(zhì)心，t為實(shí)際形狀的質(zhì)心。定義初始的相對位置為實(shí)際的相對位置為，接下來通過找到最優(yōu)的線性變換A來簡化最優(yōu)旋轉(zhuǎn)矩陣R的求解。公式(7)可以轉(zhuǎn)變?yōu)橐罁?jù)最小二乘法則有最優(yōu)變換矩陣A：

2.2.3 固體速度

在得到了旋轉(zhuǎn)矩陣R以及平移矢量t與t0之后，粒子目標(biāo)位置gi可通過下式求得：

其中，Nr為粒子i所屬區(qū)域的個(gè)數(shù)。如圖4所示，粒子b既被劃分到區(qū)域a又被劃分到區(qū)域c，所以有Nr=2。

圖4 粒子作用區(qū)域示意圖

基于計(jì)算出的目標(biāo)位置（圖5），一個(gè)試圖恢復(fù)物體原有形狀的額外速度被引入，如下式所示：

其中，fext為粒子所受到的外力（本文僅包含重力）。Ui為粒子i所屬區(qū)域的集合。

圖5 基于目標(biāo)位置的形狀恢復(fù)示意圖

3 血流速度求解

由2.1節(jié)和2.2節(jié)可計(jì)算出兩種算法下血粒子的速度，然后通過線性插值[11]來求得血粒子的實(shí)際速度。所謂線性插值即通過參數(shù)ω和(1-ω)賦予兩個(gè)速度值不同的權(quán)重且對兩個(gè)速度值進(jìn)行線性疊加，如式(12)所示：當(dāng)ω=0時(shí)，速度部分只有一項(xiàng)，即為形狀約束算法求得的速度，該速度試圖恢復(fù)血液的原始形狀而使血液保持固體的屬性；當(dāng)1-ω=0（即ω=1）時(shí)，速度部分也只有一項(xiàng)，即為不考慮粘滯力情況下SPH算法求得的速度，該速度因?yàn)椴豢紤]粘滯性而具有理想水的屬性；當(dāng)ω的取值為0～1時(shí)，線性疊加后的速度具有理想水和彈性固體之間的屬性，且ω的值越大越接近理想水的屬性，反之越接近彈性固體的屬性。血液的屬性處于理想流體水與彈性固體之間，是一種粘彈性流體，所以該算法對于血流仿真具有可行性。

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

由方程組(1)中的動量守恒方程可以看出，在用SPH算法仿真血流時(shí)，粘滯力主要取決于血液粘滯系數(shù)μ。在本研究的算法中，則是通過調(diào)整插值系數(shù)ω進(jìn)行仿真實(shí)驗(yàn)來確定血液粘性所對應(yīng)的近似ω值。

4 算法思路

本研究算法主要包括以下幾個(gè)步驟：首先是算法的初始化，將整個(gè)血流仿真系統(tǒng)初始化為一個(gè)粒子系統(tǒng)，并初始化血粒子的參數(shù)信息（比如血粒子的臨床參數(shù)信息，初始位置、速度以及加速度信息等）。其次是使用SPH算法求解血流控制方程進(jìn)而求解血流流體速度的部分，在該過程中要循環(huán)計(jì)算每一血粒子的密度、內(nèi)部壓力以及自身重力來求得粒子所受的合力，根據(jù)式(5)計(jì)算得出血粒子的加速度，再由式(6)計(jì)算得出血粒子的流體速度；然后利用形狀約束算法求解血流固體速度的部分，在該過程中要循環(huán)計(jì)算旋轉(zhuǎn)矩陣R和平移矢量t，根據(jù)式(10)求解血粒子的目標(biāo)位置，進(jìn)而根據(jù)式(11)求得血粒子的固體速度。最后根據(jù)式(12)得出血粒子的實(shí)際速度并更新血粒子的位置信息，使用Marching Cubes算法可視化處理后即可得出粘彈性流體仿真效果。具體步驟如下：

（1）初始化仿真環(huán)境：初始化血液參數(shù)信息，初始化血粒子的原始位置信息、原始速度以及原始加速度信息。

（2）確定血流的牛頓流體模型，通過SPH方法將連續(xù)的流體動力學(xué)方程轉(zhuǎn)變?yōu)殡x散形式，并求解離散方程得到血流的流體速度。

（3）利用形狀約束算法控制血粒子的運(yùn)動規(guī)律，求解血流的固體速度。

（4）通過線性插值系數(shù)ω得到流體的實(shí)際速度ui。

（5）根據(jù)血粒子的速度信息更新血粒子的位置信息。

（6）使用Marching Cubes算法實(shí)現(xiàn)可視化效果。

（7）調(diào)整ω進(jìn)行重復(fù)的仿真實(shí)驗(yàn)，直到得到類似血液的仿真結(jié)果。

整體算法流程，見圖6。

圖6 整體算法流程圖

5 實(shí)驗(yàn)結(jié)果與分析

5.1 實(shí)驗(yàn)環(huán)境

本研究實(shí)驗(yàn)的硬件環(huán)境：CPU雙核1.60 GHz，RAM 4 GB，顯卡NVIDIA GeForce GT 620M。

本研究實(shí)驗(yàn)的軟件環(huán)境：操作系統(tǒng)為Windows系統(tǒng)（64位），開發(fā)環(huán)境為Microsoft Visual Studio 2012，算法實(shí)現(xiàn)的語言采用C++，實(shí)驗(yàn)結(jié)果渲染采用OpenGL來實(shí)現(xiàn)。

5.2 實(shí)驗(yàn)結(jié)果與分析

在仿真實(shí)驗(yàn)中，SPH算法中支持域取單倍的光滑長度，形狀約束算法區(qū)域劃分長度也取單倍的光滑長度。為計(jì)算方便，血粒子質(zhì)量mi取值為1.0，粒子信息更新的時(shí)間步長Δt取為0.02 s。為了較好地進(jìn)行實(shí)時(shí)性的對比，實(shí)驗(yàn)中所提到的方法均用來仿真皮膚表面切割后的血液流動效果。實(shí)驗(yàn)結(jié)果中坐標(biāo)系O-XYZ為世界坐標(biāo)系，即虛擬的手術(shù)環(huán)境，XOZ平面為模擬的皮膚表面。皮膚表面切口效果如圖7所示；血流仿真結(jié)果如圖8～11所示。

圖7 皮膚表面切口效果示意圖

圖8 本研究算法血流仿真效果（粒子總數(shù)1000，ω=0．8）

圖9 傳統(tǒng)SPH算法血流仿真效果（粒子總數(shù)1000，血流粘滯系數(shù)μ=4 mPa · s）

圖10 算法血流仿真效果[2]（粒子總數(shù)1000，血流粘滯系數(shù)μ=4 mPa · s）

圖11 本研究算法血流仿真效果（粒子總數(shù)2000，ω=0．8）

由實(shí)驗(yàn)結(jié)果可知，當(dāng)ω=0.8時(shí)，仿真效果最近似為血流仿真效果。對比圖8、圖9和圖10，可以得到如下結(jié)論：在相同的條件下（濕度、溫度、血流仿真模型取相同時(shí)），相同切口大?。捶抡鏁r(shí)所采用的粒子總數(shù)相同）時(shí)，本研究所采用方法的實(shí)時(shí)性與傳統(tǒng)SPH方法相比有很大的提升，與文獻(xiàn)[2]中的方法相比仿真效率也有很大的優(yōu)勢，幀率近似為傳統(tǒng)SPH方法的3倍、文獻(xiàn)[2]中方法的1.5倍。對比圖8和圖11，可以得到如下結(jié)論：在相同的條件下（濕度、溫度、血流仿真模型取相同值），切口大?。捶抡鏁r(shí)所采用的粒子總數(shù)不相同）不同時(shí)，切口越大（即仿真時(shí)所采用的粒子總數(shù)越多），實(shí)時(shí)性越低，仿真效率越低，但連續(xù)性較強(qiáng)，真實(shí)感較好。

本研究算法的實(shí)時(shí)性主要取決于血粒子目標(biāo)位置的計(jì)算，尤其是當(dāng)血粒子數(shù)目比較多的時(shí)候，血粒子位置計(jì)算所消耗的時(shí)間也就相對較多，但較近些年的血流仿真方法在效率上還是有很大提升的。表1是兩種粒子數(shù)情況下3種算法的計(jì)算效率（渲染部分均采用Marching Cubes算法）。結(jié)果表明，本研究算法能夠高效、快速地仿真皮膚表面的血液流動，降低了計(jì)算過程的復(fù)雜度，在保證真實(shí)感的前提下提高了血流模擬的實(shí)時(shí)性。

表1 兩種粒子數(shù)情況下3種算法的計(jì)算效率

6 結(jié)論

流血現(xiàn)象廣泛存在于我們的日常生活中，比如意外劃傷出血、手術(shù)中皮膚表面切割出血、動脈出血、器官表面滲血等，因此對于血流的仿真就成為虛擬現(xiàn)實(shí)領(lǐng)域的一大熱點(diǎn)。又因?yàn)閷?shí)時(shí)性和真實(shí)感的要求較高，血流仿真也是醫(yī)學(xué)仿真領(lǐng)域的一大難點(diǎn)。本研究在SPH算法模擬流體的基礎(chǔ)上加入形狀約束算法來控制血粒子的運(yùn)動規(guī)律，并通過調(diào)節(jié)插值系數(shù)來實(shí)現(xiàn)血流的實(shí)時(shí)仿真，可視化部分采用經(jīng)典的Marching Cubes算法，最終在保證真實(shí)感的前提下實(shí)現(xiàn)了血流現(xiàn)象的實(shí)時(shí)仿真。未來工作可以考慮使用GPU對本研究算法進(jìn)行加速，也可使用本研究算法模擬更復(fù)雜的場景。

[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] 楊禮波．虛擬手術(shù)系統(tǒng)中的流血模擬[D]．鄭州:華北水利水電學(xué)院,2011．

[5] 溫嬋娟,歐嘉蔚,賈金原．GPU通用計(jì)算平臺上的SPH流體模擬[J]．計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2010,22(3):407-411．

[6] 施鵬,熊岳山,徐凱,等．虛擬肝臟手術(shù)中實(shí)時(shí)動態(tài)滲血效果模擬[J]．計(jì)算機(jī)應(yīng)用,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] 萬旺根,林繼承,余小清,等．基于粒子系統(tǒng)和形狀匹配的實(shí)時(shí)無網(wǎng)格變形仿真[J]．計(jì)算機(jī)應(yīng)用,2008,28(12):3007-3009．

[11] 唐勇,陳靜,呂夢雅．基于SPH及形狀約束的粘彈性流體的實(shí)時(shí)模擬[J]．小型微型計(jì)算機(jī)系統(tǒng),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

福建省科技計(jì)劃重點(diǎn)項(xiàng)目支持（2011H0027）；福建省自然科學(xué)基金項(xiàng)目支持（2013J05090）。

陳國棟，副研究員，主要研究方向?yàn)橛?jì)算機(jī)圖形學(xué)和虛擬現(xiàn)實(shí)技術(shù)。

通訊作者郵箱：1285056478@qq．com