一種硅油填充術眼內硅油乳化過程模擬可視化方法

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

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

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

③(材料領域知識工程北京市重點實驗室(北京科技大學) 北京 100083)

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

1 引言

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

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

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

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

本文針對治療孔源性視網膜脫離的玻璃體切割聯合硅油填充術，提出一種眼內硅油乳化過程模擬可視化方法，對硅油填充過程及乳化過程的多相流體動態交互與互溶擴散過程進行模擬。為處理眼內硅油-水兩種不同密度流體交互，提出體積不可壓縮多相流體計算框架進行耦合模擬；結合內聚力與曲率力描述宏觀兩相流體耦合中的表面張力作用；引入力平衡散體動力學模型，描述硅油乳化過程中的互溶擴散效應。

2 基于SPH方法的流體模擬基礎

2.1 流體力學控制方程

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

2.2 SPH離散化

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

圖1 SPH數值近似示意圖

面向硅油填充術中的硅油-水多相流體交互問題，表面張力所導致的交互界面曲率效果，以及硅油乳化效應所產生的多相流體互溶擴散運動，本節基于傳統用于單相流體模擬過程的SPH不可壓縮流體模擬方法，構建了硅油-水耦合乳化可視化模擬框架，如圖2所示。首先，針對傳統SPH模擬方法在處理非均一靜態密度流場的數值計算誤差問題，使用基于體積不可壓縮的多相流體壓強計算模型，增強數值計算準確性；其次，采用內聚力與曲率力模擬兩相流體間表面張力作用；再次，引入力平衡散體動力學模型，模擬硅油-水兩相間因乳化產生的互溶擴散效應。

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

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

時下最先進的無散度隱式不可壓縮方法DFSPH[14]通過流體壓縮狀態描述壓強數值。即根據流體密度壓縮狀態判定壓強大小，抵消其他項導致的壓縮，使流體整體不可壓縮。根據式(2)，有

3.2 相間表面張力建模

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

對于內聚力，可直接表現為同種粒子間的相互吸引作用，在考慮牛頓第二定律的基礎上，引入符合高斯曲線形狀的吸引規律

即為A=1情況下的SPH方法數值計算過程。位于流體表面的流體粒子由于缺乏鄰居，故會產生顏色場偏小問題，通過式(5)可計算顏色場梯度，從而獲取界面法線值

3.3 硅油乳化力平衡散體動力學建模

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

對于多相流場景，在每個粒子內部各相混合的狀態下，粒子的速度可以描述為

對于動力學過程，本文引入力平衡方程模型[25]，將粒子間相互作用與粒子內物相交換兩個過程解耦，顯著提高了隱式SPH方法的穩定性，在粒子i中相k的速度變化率可寫成

從而可獲取粒子間與粒子內部多相流運動交互過程聯系。

圖3 具有混合相的粒子

4 實驗

為驗證本文方法的有效性，本節首先進行硅油-水兩相流體耦合實驗，并驗證兩相流體耦合中的表面張力作用；然后開展墨滴擴散實驗，模擬可混溶多相流相間交互效果；最后對眼球硅油填充進行仿真，進行術后硅油乳化的模擬。本文仿真算法基于C++編寫，利用Eigen作為數學計算工具，并使用文獻[26]的表面重構方法由粒子生成3D網格?？梢暬矫?，使用OpenGL作為流體粒子可視化工具，顯示實時仿真效果，并使用Blender進行離線動畫渲染。本文開展實驗所用的計算配置是一臺128 GB內存工作站，處理器為2個16核CPU，主頻為2.30 GHz。

4.1 瑞利泰勒不穩定性實驗

本節通過實驗驗證體積不可壓縮SPH方法模擬兩相流體耦合的優越性，以及表面張力在流體耦合中的作用。瑞利泰勒不穩定性實驗如圖4所示。圖4(a)—圖4(d)使用DFSPH方法進行模擬，圖4(e)—圖4(h)使用本文方法，圖4(i)—圖4(l)使用本文方法結合表面張力作用，黃色相為硅油，藍色相為水，兩相流體不可混溶，密度比為0.7:1，硅油相因為密度較小而上浮，從而形成對流效果。

圖4 瑞利泰勒不穩定性實驗

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

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

圖5 DFSPH和本文方法迭代次數對比

該實驗證明了本文方法可有效模擬眼球內硅油-水兩相耦合，并表現了宏觀兩相流體耦合過程中表面張力的作用。

4.2 墨滴實驗

本實驗模擬紅色墨滴落入水中的擴散過程，展示可混溶多相流體相間交互效果，以輔助硅油乳化過程可視化分析。紅色相為墨滴，透明相為水，密度比為5:1，墨滴墜落的過程中會產生相間擴散效果，如圖6所示。

該實驗表明本文的力平衡散體動力學模型可有效模擬可混溶多相流體相間交互效果，以輔助實現硅油乳化過程可視量化分析。

4.3 硅油注入與乳化實驗

本實驗模擬硅油填充術的硅油注入過程及術后的硅油乳化過程。該實驗通過力平衡散體動力學模型中的擴散系數控制硅油乳化程度，并考慮了表面張力作用，表面張力系數設置為3。圖7為術后同一時刻眼內腔硅油乳化不同程度的仿真結果，(a)圖擴散系數為0.01，硅油乳化較輕，(b)圖擴散系數為0.1，硅油乳化情況較為嚴重。

該實驗驗證了本文方法可以對硅油填充術后眼內腔硅油乳化過程進行有效模擬，以輔助術者判定所需合適的硅油注入量及術后硅油取出時間。

5 結束語

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

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