一種硅油填充術(shù)眼內(nèi)硅油乳化過(guò)程模擬可視化方法

班曉娟 王佳敏王笑琨* 張雅斕徐衍睿宋重明黃厚斌 朱志鴻

①(北京材料基因工程高精尖創(chuàng)新中心(北京科技大學(xué)) 北京 100083)

②(計(jì)算機(jī)與通信工程學(xué)院、人工智能研究院(北京科技大學(xué)) 北京 100083)

③(材料領(lǐng)域知識(shí)工程北京市重點(diǎn)實(shí)驗(yàn)室(北京科技大學(xué)) 北京 100083)

④(中國(guó)人民解放軍總醫(yī)院海南醫(yī)院 三亞 572013)

1 引言

孔源性視網(wǎng)膜脫離(Rhegmatogenous Retinal Detachment， RRD)以視網(wǎng)膜裂孔形成為特征，是視網(wǎng)膜神經(jīng)上皮層和視網(wǎng)膜色素上皮層分離的一類(lèi)視網(wǎng)膜病變。玻璃體切割(Pars Plana Vitrectomy，PPV)聯(lián)合硅油填充手術(shù)是孔源性視網(wǎng)膜脫離的主要治療方式。硅油是一種安全有效的眼內(nèi)填充劑[1]，其理化性質(zhì)穩(wěn)定，沒(méi)有生物毒性相對(duì)安全，良好的透光性和屈光性不影響術(shù)后對(duì)眼底的觀(guān)察與治療，有利于視力恢復(fù)，強(qiáng)表面張力可以有效封閉視網(wǎng)膜裂孔，并防止眼內(nèi)低張力液體穿透裂孔造成視網(wǎng)膜2次脫離。但是，硅油長(zhǎng)期存在于眼內(nèi)會(huì)發(fā)生乳化，乳化部分表面張力降低。填充過(guò)少，可能導(dǎo)致剩余硅油不足以完整覆蓋視網(wǎng)膜裂孔，降低手術(shù)成功率；填充過(guò)多，乳化硅油進(jìn)入前房，會(huì)導(dǎo)致一系列眼部并發(fā)癥，如增生性玻璃體視網(wǎng)膜病變、復(fù)發(fā)性視網(wǎng)膜脫離、眼內(nèi)炎癥、繼發(fā)性青光眼等[2]。因此，考慮硅油乳化情況并適量填充硅油，預(yù)測(cè)硅油乳化狀態(tài)對(duì)手術(shù)預(yù)后的影響，是提高手術(shù)成功率、減少眼部并發(fā)癥發(fā)生的關(guān)鍵。

隨著醫(yī)學(xué)與計(jì)算科學(xué)的交叉學(xué)科發(fā)展，醫(yī)學(xué)可視化研究與計(jì)算機(jī)技術(shù)輔助醫(yī)療診斷已取得了長(zhǎng)足進(jìn)展。蔡軼珩等人[3]利用機(jī)器學(xué)習(xí)方法，將血管像素與非血管像素看作二分類(lèi)，實(shí)現(xiàn)從背景中迅速分割出視網(wǎng)膜血管。陳強(qiáng)等人[4]利用隨機(jī)森林分類(lèi)器，通過(guò)頻譜域光學(xué)相干層析技術(shù)對(duì)視網(wǎng)膜神經(jīng)纖維層進(jìn)行分割，以輔助診斷青光眼等疾病。最近徐衍睿等人[5]提出一種針對(duì)硅油填充手術(shù)的3維建模與可視化方案輔助醫(yī)生進(jìn)行手術(shù)流程規(guī)劃，但是該方法的相間壓力計(jì)算并不精確，并且沒(méi)有考慮硅油乳化的擴(kuò)散現(xiàn)象。本文利用醫(yī)學(xué)影像眼球3維模型，提出一種計(jì)算機(jī)圖形學(xué)基于物理的流體模擬方法，對(duì)眼內(nèi)硅油填充過(guò)程以及乳化過(guò)程進(jìn)行可視化仿真建模，從而輔助醫(yī)生判斷硅油填充量和取出時(shí)間，以減少眼部并發(fā)癥的發(fā)生，提升手術(shù)預(yù)后。

基于物理的流體模擬可視化方法可分為歐拉網(wǎng)格法、拉格朗日粒子法以及混合方法[6]。歐拉網(wǎng)格法將空間離散化為網(wǎng)格，可簡(jiǎn)化多相流場(chǎng)景相間動(dòng)量交換的計(jì)算過(guò)程，便于耦合數(shù)據(jù)驅(qū)動(dòng)方法以加快計(jì)算速度[7]，但容易產(chǎn)生數(shù)值耗散；拉格朗日粒子法將流體和剛體離散化為粒子，更適合模擬具有復(fù)雜流變性的液體[8]以及流體自由表面的精細(xì)動(dòng)態(tài)[9]。其中光滑粒子流體動(dòng)力學(xué)(Smoothed Particle Hydrodynamics， SPH)方法采用粒子積分近似[10]，能夠捕捉流體各方面動(dòng)態(tài)細(xì)節(jié)。Becker等人[11]提出基于Tait方程的弱可壓縮SPH(Weakly Compressible SPH， WCSPH)方法以逼近流體不可壓縮特性。隨后Solenthaler等人[12]提出預(yù)測(cè)-矯正SPH(Predictive-Corrective Incompressible SPH， PCISPH)方法，通過(guò)預(yù)測(cè)-矯正隱式迭代方案確定流場(chǎng)壓強(qiáng)，相較于WCSPH有了較大的效率提升。Ihmsen等人[13]提出隱式不可壓縮SPH(Implicit Incompressible SPH， IISPH)方法，進(jìn)一步提高了求解器的收斂速度與魯棒性。在此基礎(chǔ)上，Bender等人[14]組合無(wú)散度求解器和恒定密度求解器，提出了無(wú)散度SPH(Divergence-Free SPH， DFSPH)方法，進(jìn)一步提升模擬性能。對(duì)于SPH方法在預(yù)測(cè)-矯正過(guò)程中存在的數(shù)值耗散問(wèn)題，Wang等人[15]提出一種基于朗肯渦模型的湍流精細(xì)化方法，通過(guò)恢復(fù)由于粒子缺乏旋轉(zhuǎn)而損失的能量來(lái)增強(qiáng)表面細(xì)節(jié)。進(jìn)一步，Liu等人[16]提出了一種基于流函數(shù)的渦度恢復(fù)方法，將渦量場(chǎng)和速度場(chǎng)聯(lián)系起來(lái)，該方法可以增強(qiáng)現(xiàn)有的渦旋并且捕獲額外的湍流。

物理學(xué)研究中，“相”指不同物態(tài)或同一物態(tài)的不同物理性質(zhì)或力學(xué)狀態(tài)，多相流即系統(tǒng)中同時(shí)存在兩個(gè)或兩個(gè)以上的流體相態(tài)。硅油填充手術(shù)過(guò)程中，涉及硅油與水兩相液體的交互運(yùn)動(dòng)，需要考慮不同液體的物理性質(zhì)，以實(shí)現(xiàn)混溶、不混溶等復(fù)雜現(xiàn)象的模擬。與歐拉網(wǎng)格法相比，粒子表示方法能更加清晰地表示流體界面，更適合應(yīng)用于多相流模擬。Müller等人[17]采用擴(kuò)散方程，使用粒子方法實(shí)現(xiàn)了可混溶液體的擴(kuò)散模擬。為了進(jìn)一步豐富多相流交互模擬效果，Liu等人[18]引入體積分?jǐn)?shù)方案與SPH方法相結(jié)合。這些方案僅利用濃度差驅(qū)動(dòng)不同相之間的擴(kuò)散效應(yīng)，忽略了流體運(yùn)動(dòng)和作用力分布導(dǎo)致的多相流體的混合與分離。Ren等人[19]將漂移速度引入SPH公式，成功模擬了由于各相之間的相對(duì)運(yùn)動(dòng)而引起的混合和非混合現(xiàn)象。然而，該方法將質(zhì)心速度作為混合粒子的速度，這會(huì)造成非無(wú)散度速度場(chǎng)和粒子內(nèi)部膨脹，并且該方法與不可壓縮求解器不兼容。Yang等人[20]基于亥姆霍茲自由能，擴(kuò)展了多相流模擬方法，以靈活地捕捉各種混合與分離過(guò)程。Yan等人[21]進(jìn)一步提出多相流與剛體耦合方法，實(shí)現(xiàn)可變形體、顆粒材料、多種流體之間的相互作用效果模擬。Jiang等人[22]在體積加權(quán)混合速度的基礎(chǔ)上，建立了一種新的多相流無(wú)散度混合模型，通過(guò)自適應(yīng)單相不可壓縮求解器求解混合速度，并確保質(zhì)量守恒，大大提高了模擬的精度。但是該方法數(shù)密度的簡(jiǎn)單實(shí)現(xiàn)限制了可以處理的多相流密度比，并且引入了非守恒的粒子內(nèi)界面動(dòng)量源，使得模擬過(guò)程不穩(wěn)定。

本文針對(duì)治療孔源性視網(wǎng)膜脫離的玻璃體切割聯(lián)合硅油填充術(shù)，提出一種眼內(nèi)硅油乳化過(guò)程模擬可視化方法，對(duì)硅油填充過(guò)程及乳化過(guò)程的多相流體動(dòng)態(tài)交互與互溶擴(kuò)散過(guò)程進(jìn)行模擬。為處理眼內(nèi)硅油-水兩種不同密度流體交互，提出體積不可壓縮多相流體計(jì)算框架進(jìn)行耦合模擬；結(jié)合內(nèi)聚力與曲率力描述宏觀(guān)兩相流體耦合中的表面張力作用；引入力平衡散體動(dòng)力學(xué)模型，描述硅油乳化過(guò)程中的互溶擴(kuò)散效應(yīng)。

2 基于SPH方法的流體模擬基礎(chǔ)

2.1 流體力學(xué)控制方程

物理上用于描述不可壓縮流體動(dòng)量守恒的納維-斯托克斯(Navier-Stocks， N-S)方程為

2.2 SPH離散化

3 可混溶多相流體耦合求解器

圖1 SPH數(shù)值近似示意圖

面向硅油填充術(shù)中的硅油-水多相流體交互問(wèn)題，表面張力所導(dǎo)致的交互界面曲率效果，以及硅油乳化效應(yīng)所產(chǎn)生的多相流體互溶擴(kuò)散運(yùn)動(dòng)，本節(jié)基于傳統(tǒng)用于單相流體模擬過(guò)程的SPH不可壓縮流體模擬方法，構(gòu)建了硅油-水耦合乳化可視化模擬框架，如圖2所示。首先，針對(duì)傳統(tǒng)SPH模擬方法在處理非均一靜態(tài)密度流場(chǎng)的數(shù)值計(jì)算誤差問(wèn)題，使用基于體積不可壓縮的多相流體壓強(qiáng)計(jì)算模型，增強(qiáng)數(shù)值計(jì)算準(zhǔn)確性；其次，采用內(nèi)聚力與曲率力模擬兩相流體間表面張力作用；再次，引入力平衡散體動(dòng)力學(xué)模型，模擬硅油-水兩相間因乳化產(chǎn)生的互溶擴(kuò)散效應(yīng)。

圖2 可混溶硅油-水耦合乳化模擬框架

3.1 體積不可壓縮的多相流體模擬

時(shí)下最先進(jìn)的無(wú)散度隱式不可壓縮方法DFSPH[14]通過(guò)流體壓縮狀態(tài)描述壓強(qiáng)數(shù)值。即根據(jù)流體密度壓縮狀態(tài)判定壓強(qiáng)大小，抵消其他項(xiàng)導(dǎo)致的壓縮，使流體整體不可壓縮。根據(jù)式(2)，有

3.2 相間表面張力建模

為模擬硅油液滴的強(qiáng)表面張力作用，本文引入單相流體間的內(nèi)聚力(cohesion)以及交互過(guò)程中令表面積最小化的曲率力(curvature)[24]。

對(duì)于內(nèi)聚力，可直接表現(xiàn)為同種粒子間的相互吸引作用，在考慮牛頓第二定律的基礎(chǔ)上，引入符合高斯曲線(xiàn)形狀的吸引規(guī)律

即為A=1情況下的SPH方法數(shù)值計(jì)算過(guò)程。位于流體表面的流體粒子由于缺乏鄰居，故會(huì)產(chǎn)生顏色場(chǎng)偏小問(wèn)題，通過(guò)式(5)可計(jì)算顏色場(chǎng)梯度，從而獲取界面法線(xiàn)值

3.3 硅油乳化力平衡散體動(dòng)力學(xué)建模

對(duì)于硅油和水之間表現(xiàn)為互溶擴(kuò)散現(xiàn)象的乳化效果，本文采取體積分?jǐn)?shù)方案，令同一離散粒子內(nèi)部同時(shí)存在多相流體。相較于Jiang等人[22]提出的無(wú)散度混合模型，本文做了兩個(gè)改進(jìn)，分別是粒子內(nèi)相交換過(guò)程和粒子間預(yù)測(cè)-矯正過(guò)程的解耦及隱式計(jì)算漂移速度，使得本文方法在滿(mǎn)足質(zhì)量守恒的同時(shí)滿(mǎn)足粒子內(nèi)動(dòng)量守恒，提高了可混溶多相流體交互的穩(wěn)定性，并且該方法可以便利地集成到各種不可壓縮流體求解器中。

對(duì)于多相流場(chǎng)景，在每個(gè)粒子內(nèi)部各相混合的狀態(tài)下，粒子的速度可以描述為

對(duì)于動(dòng)力學(xué)過(guò)程，本文引入力平衡方程模型[25]，將粒子間相互作用與粒子內(nèi)物相交換兩個(gè)過(guò)程解耦，顯著提高了隱式SPH方法的穩(wěn)定性，在粒子i中相k的速度變化率可寫(xiě)成

從而可獲取粒子間與粒子內(nèi)部多相流運(yùn)動(dòng)交互過(guò)程聯(lián)系。

圖3 具有混合相的粒子

4 實(shí)驗(yàn)

為驗(yàn)證本文方法的有效性，本節(jié)首先進(jìn)行硅油-水兩相流體耦合實(shí)驗(yàn)，并驗(yàn)證兩相流體耦合中的表面張力作用；然后開(kāi)展墨滴擴(kuò)散實(shí)驗(yàn)，模擬可混溶多相流相間交互效果；最后對(duì)眼球硅油填充進(jìn)行仿真，進(jìn)行術(shù)后硅油乳化的模擬。本文仿真算法基于C++編寫(xiě)，利用Eigen作為數(shù)學(xué)計(jì)算工具，并使用文獻(xiàn)[26]的表面重構(gòu)方法由粒子生成3D網(wǎng)格?？梢暬矫妫褂肙penGL作為流體粒子可視化工具，顯示實(shí)時(shí)仿真效果，并使用Blender進(jìn)行離線(xiàn)動(dòng)畫(huà)渲染。本文開(kāi)展實(shí)驗(yàn)所用的計(jì)算配置是一臺(tái)128 GB內(nèi)存工作站，處理器為2個(gè)16核CPU，主頻為2.30 GHz。

4.1 瑞利泰勒不穩(wěn)定性實(shí)驗(yàn)

本節(jié)通過(guò)實(shí)驗(yàn)驗(yàn)證體積不可壓縮SPH方法模擬兩相流體耦合的優(yōu)越性，以及表面張力在流體耦合中的作用。瑞利泰勒不穩(wěn)定性實(shí)驗(yàn)如圖4所示。圖4(a)—圖4(d)使用DFSPH方法進(jìn)行模擬，圖4(e)—圖4(h)使用本文方法，圖4(i)—圖4(l)使用本文方法結(jié)合表面張力作用，黃色相為硅油，藍(lán)色相為水，兩相流體不可混溶，密度比為0.7:1，硅油相因?yàn)槊芏容^小而上浮，從而形成對(duì)流效果。

圖4 瑞利泰勒不穩(wěn)定性實(shí)驗(yàn)

第1行是使用DFSPH方法進(jìn)行模擬得到的結(jié)果，第2行使用本文所提體積不可壓縮多相流體壓強(qiáng)計(jì)算模型，第3行在體積不可壓縮模型計(jì)算壓強(qiáng)的同時(shí)處理了表面張力的作用。對(duì)比3種方法不同時(shí)刻模擬結(jié)果可知，DFSPH密度不可壓縮計(jì)算方法由于相間壓強(qiáng)計(jì)算誤差產(chǎn)生了嚴(yán)重的滲透現(xiàn)象，本文體積不可壓縮計(jì)算方法相間壓強(qiáng)計(jì)算更加穩(wěn)定，可以保留更多的耦合細(xì)節(jié)；考慮表面張力作用，硅油相產(chǎn)生表面積最小化趨勢(shì)，可以有效緩解低密度比流體對(duì)流效果的紊亂狀態(tài)，避免產(chǎn)生硅油小滴。

圖5對(duì)比了DFSPH和本文方法在該場(chǎng)景下每幀壓強(qiáng)求解器的迭代次數(shù)。DFSPH方法由于沒(méi)有表面張力，兩相流體對(duì)流效果紊亂，易產(chǎn)生硅油小滴，相間壓強(qiáng)計(jì)算不穩(wěn)定，迭代次數(shù)周期性出現(xiàn)峰值；本文方法考慮了硅油相的表面張力作用，可以穩(wěn)定模擬兩相流體對(duì)流效果，在流體塊碰撞容器壁后迭代次數(shù)趨于穩(wěn)定。

圖5 DFSPH和本文方法迭代次數(shù)對(duì)比

該實(shí)驗(yàn)證明了本文方法可有效模擬眼球內(nèi)硅油-水兩相耦合，并表現(xiàn)了宏觀(guān)兩相流體耦合過(guò)程中表面張力的作用。

4.2 墨滴實(shí)驗(yàn)

本實(shí)驗(yàn)?zāi)M紅色墨滴落入水中的擴(kuò)散過(guò)程，展示可混溶多相流體相間交互效果，以輔助硅油乳化過(guò)程可視化分析。紅色相為墨滴，透明相為水，密度比為5:1，墨滴墜落的過(guò)程中會(huì)產(chǎn)生相間擴(kuò)散效果，如圖6所示。

該實(shí)驗(yàn)表明本文的力平衡散體動(dòng)力學(xué)模型可有效模擬可混溶多相流體相間交互效果，以輔助實(shí)現(xiàn)硅油乳化過(guò)程可視量化分析。

4.3 硅油注入與乳化實(shí)驗(yàn)

本實(shí)驗(yàn)?zāi)M硅油填充術(shù)的硅油注入過(guò)程及術(shù)后的硅油乳化過(guò)程。該實(shí)驗(yàn)通過(guò)力平衡散體動(dòng)力學(xué)模型中的擴(kuò)散系數(shù)控制硅油乳化程度，并考慮了表面張力作用，表面張力系數(shù)設(shè)置為3。圖7為術(shù)后同一時(shí)刻眼內(nèi)腔硅油乳化不同程度的仿真結(jié)果，(a)圖擴(kuò)散系數(shù)為0.01，硅油乳化較輕，(b)圖擴(kuò)散系數(shù)為0.1，硅油乳化情況較為嚴(yán)重。

該實(shí)驗(yàn)驗(yàn)證了本文方法可以對(duì)硅油填充術(shù)后眼內(nèi)腔硅油乳化過(guò)程進(jìn)行有效模擬，以輔助術(shù)者判定所需合適的硅油注入量及術(shù)后硅油取出時(shí)間。

5 結(jié)束語(yǔ)

玻璃體切割聯(lián)合硅油填充術(shù)在治療孔源性視網(wǎng)膜脫離手術(shù)中日趨流行，然而硅油長(zhǎng)期存在于眼內(nèi)會(huì)發(fā)生乳化，進(jìn)而導(dǎo)致角膜帶狀變性、繼發(fā)性青光眼、并發(fā)性白內(nèi)障等嚴(yán)重的眼部并發(fā)癥，所以預(yù)測(cè)硅油乳化的程度并及時(shí)取出硅油，可以在保證最好的治療效果的同時(shí)減少并發(fā)癥的發(fā)生。本文提出一種硅油填充術(shù)后眼內(nèi)硅油乳化過(guò)程計(jì)算機(jī)模擬可視化方法，構(gòu)建了體積不可壓縮多相流體計(jì)算框架進(jìn)行眼內(nèi)腔硅油-水的耦合模擬，采用內(nèi)聚力和曲率力相結(jié)合處理表面張力作用，并引入力平衡散體動(dòng)力學(xué)模型，實(shí)現(xiàn)眼內(nèi)腔硅油乳化的仿真過(guò)程。實(shí)驗(yàn)結(jié)果表明，本文方法能夠較好地實(shí)現(xiàn)硅油-水穩(wěn)定耦合、表面張力、乳化擴(kuò)散等效果的模擬與可視化。

本文為預(yù)測(cè)術(shù)后眼內(nèi)腔硅油乳化狀態(tài)提供了一種有效方法，也給物理模擬輔助眼科學(xué)帶來(lái)了新的研究思路。但是，影響硅油乳化的因素有很多，例如硅油自身的黏度、純度、揮發(fā)性等，以及眼內(nèi)出血、眼球過(guò)度運(yùn)動(dòng)等眼內(nèi)因素，如何考慮實(shí)際情況合理設(shè)置硅油乳化的難易程度，精確地模擬眼內(nèi)腔硅油乳化過(guò)程，以輔助判定硅油取出時(shí)間，是下一步研究的重點(diǎn)與難點(diǎn)。