面向视网膜脱离手术的硅油填充模拟

徐衍睿，班晓娟，王笑琨✉，王 宇，尹 豆，周 靖，黄厚斌，朱志鸿

1) 北京科技大学北京材料基因工程高精尖创新中心，北京 100083 2) 北京科技大学计算机与通信工程学院，北京 100083 3) 北京科技大学人工智能研究院，北京 100083 4) 北京科技大学材料领域知识工程北京市重点实验室，北京 100083 5) 中国人民解放军总医院海南医院，三亚 572013

孔源性视网膜脱离（Rhegmatogenous retinal detachment, RRD）是一种常见的玻璃体视网膜疾病，通常是由于视网膜萎缩变性或玻璃体牵引形成视网膜神经上皮全层裂孔，引起变性液化的玻璃体经裂孔进入视网膜下，从而导致视网膜神经上皮层与色素上皮层之间的分离[1]. 该疾病的主要治疗方式是通过玻璃体切割联合硅油眼内填充来实现视网膜复位.

在玻璃体切割联合硅油填充手术中，硅油需填充满整个玻璃体腔，顶压视网膜裂孔，防止水液进入视网膜下，有效修复视网膜. 但过量填充硅油会增加诸多并发症产生几率，如硅油乳化引起增生性玻璃体视网膜病变[2]，硅油进入前房导致角膜病变[3]，填充量过多引起瞳孔阻滞型青光眼[4]等.因此，掌握眼球内填充状态，控制最佳填充量，是目前手术过程迫切亟需解决的问题.

为此，本文采用物理模拟与三维可视化结合的方式，对面向孔源性视网膜脱离的硅油填充过程进行模拟，为硅油填充后的眼内状态分析和估测硅油填充量提供辅助参考.首先基于手术过程医学参考资料与眼球几何结构建立术中眼球三维模型，并进行离散粒子采样；然后，根据水与硅油的不同物理性质，对水-硅油两相流动及交互进行建模；最后，构建固液交互模型，实现两相液体（水和硅油）与固体边界（眼球）的交互仿真.

论文的其余部分组织如下：首先介绍相关领域的研究工作；其次对基于物理的视网膜脱离手术硅油填充模拟等主要内容进行阐述；接下来对所提方法进行验证，论述实验结果；最后归纳与总结.

1 相关工作

随着现代医疗技术的进步，孔源性视网膜脱离治疗越来越得到重视，基于玻璃体切割联合眼内硅油填充手术治疗孔源性视网膜脱离已在国际上广泛应用，手术设备和技术日益成熟. 作为常见的眼科致盲性疾病，孔源性视网膜脱离手术治疗的预后对于患者之后的生活质量至关重要，眼内硅油填充所带来的视网膜损伤、继发性青光眼等各类并发症也逐渐引起高度重视.为了提高手术质量和治疗效果，采用计算机模拟技术辅助视网膜手术分析，是有效辅助手术方案制定以及提高患者视力预后的重要手段.

1.1 视网膜脱离手术

发生孔源性视网膜脱离时，液化玻璃体经视网膜神经上皮全层裂孔进入视网膜下，导致视网膜神经上皮层与色素上皮层之间的分离，如图1所示.在上世纪60年代以前，受限于当时对玻璃体生理和病理的认识，以及精密仪器和设备的缺乏，在相当长一段时期内，玻璃体曾被认为是眼科手术的禁区. 直到20世纪80年代，随着显微外科技术和显微器械的发展，以及对玻璃体视网膜病变定义、命名和分类的进一步认识，玻璃体手术迎来了巨大的发展[5].

图1 孔源性视网膜脱离（a）与正常眼底照相（b）Fig.1 Photographs of rhegmatogenous retinal detachment (a) and normal fundus (b)

而20世纪90年代中期眼内填充物（包括空气、惰性气体、硅油和过氟化碳液体等）的发现和应用，进一步完善了玻璃体手术，为日益增加的玻璃体视网膜病变患者，提供了更好的疗效[6]. 其中硅油以其理化性质稳定、生物耐受性好等优势，现已作为一种安全有效的眼内填充剂被广泛应用[7].

在传统的玻璃体切割联合硅油填充手术中，硅油需填充满整个玻璃体腔（图2），以致引起诸多并发症. Kawaguchi[2]发现硅油乳化会引起增生性玻璃体视网膜病变、复发性视网膜脱离、眼内炎症、继发性青光眼等疾病，同时李敏等[8]使用重硅油填充治疗下方裂孔源性视网膜脱离伴严重增生性玻璃体视网膜病变(PVR)，并对手术效果及并发症进行了评价，Odrobina等[9]确定硅油的长期存留还会诱发并发性白内障、减少脉络膜厚度、引起视网膜变性，Roca等[10]尤其提出对黄斑部视网膜的顶压会严重影响视力预后.Moussa等[3]发现如果硅油进入前房会导致角膜病变（如角膜内皮失代偿、角膜带状变性），填充量过多会引起瞳孔阻滞型青光眼[4].

图2 玻璃体切割联合硅油填充术示意图. （a）孔源性视网膜脱离；（b）硅油眼内填充Fig.2 Schematic diagram of vitrectomy combined with silicone oil tamponade: (a) rhegmatogenous retinal detachment; (b) silicone oil tamponade

近年来随着医疗技术发展，掌握眼球内及填充状态，以最少量的硅油达到最佳的治疗效果，是当前研究的热点问题. 计算机技术已在医疗辅助方向和图像识别与分析领域取得了极大的进展. 童何俊和付冬梅[11]提出的基于参考模型的视网膜特征量化方法提供了一系列适用于计算机判断分析视网膜状态的可量化特征；马博渊等[12]总结了深度学习算法在图像区域组织分割中的应用，分析了如何通过提升图像分割精度来提升三维组织重构准确性.而基于物理的流体模拟能够较好的仿真眼球内的情况和分析硅油与眼球间的物理属性，通过模拟眼球内压力、硅油表面张力、硅油比重等物理参数进一步探讨眼球硅油手术填充情况，可以尽可能地辅助医生决策，减少术后硅油并发症的发生.

1.2 基于物理的流体模拟

计算机图形学领域中，流体模拟技术在过去数十年一直备受关注，得到了长足发展. 传统流体模拟方法可分为欧拉方法与拉格朗日方法，拉格朗日流体模拟方法尤为在处理质量守恒及自由表面相关问题具有显著优势.其中光滑粒子动力学(Smoothed particle hydrodynamics, SPH)将连续计算空间进行离散化，成为互相作用的粒子，并保持模拟过程中的质量守恒，适用于流体运动过程中捕捉流体各方面细节. 下面将在不可压缩模拟，多相流模拟，固液交互模拟方向进行深入论述.

SPH方法最初采用的是理想气体方程，该方法在视觉上会有很强的压缩感，但也会使得结果不那么逼真[13].Becker和Teschner[14]采用Tait方程取代之前的理想气体方程并结合硬度很高的控制系数，称之为WCSPH（Weakly compressible SPH），其能够明显增加效果真实性，但是也会降低算法的效率.为在进一步提升模拟过程中流体不可压缩性的同时让仿真过程更为高效稳定，隐式动力学计算模型相继被提出.Solenthaler和Pajarola[15]提出的PCISPH（Predictive-corrective incompressible SPH），其能够设置全局最大密度波动，用预测-矫正的迭代过程实现不可压缩性.相比于预测矫正方法，隐式不可压缩SPH（Implicit incompressible SPH, IISPH）[16]通过构造适合隐式求解的压力泊松方程（Pressure poisson equation, PPE），并使用共轭梯度法或松弛雅各比方法求解线性方程组，实现了更为高效仿真算法. 相较于PCISPH，IISPH在迭代时间步长、收敛稳定性、计算效率上均有所提升，相同场景下整体仿真速度可实现3～5倍提升.

除对于流体本身的动力学计算外，流体仿真中还需考虑，流体与固体的耦合交互作用，设立合理的边界处理条件.Monaghan和Kajtar[17]改良了早期对于流体的惩罚力方法，通过引入三次样条函数建立了排斥力模型，提升了交互边界的稳定性.随后Becker和Tessendorf[18]在固液间构建了直接力模型，防止流体粒子穿过边界，进行了对于不同滑移条件的模拟，并将隐式求解方法运用到流固耦合处理中来，基于预测-矫正方案利用直接力矫正粒子的速度位置信息，实现了固液双向耦合效果.为进一步提升流体粒子匮乏以及不均匀分布时的固液耦合表现，Harada等[19]提出了基于固液粒子距离的密度权重函数，以弥补空缺邻居粒子导致的数值计算误差，并进一步解决了流体黏附问题. Akinci等[20]利用泊松盘采样方法对于刚体模型进行单层采样，并利用镜像方法计算刚体粒子对于流体粒子的数值贡献，通过保证系统动量守恒提升了模拟过程真实性. Macklin等[21]进一步统一固体和流体仿真算法，提出了基于位置的流体粒子动力学（Position based fluids, PBF）模型，通过统一粒子约束实现相互作用.

区别于单相流体，多相流体模拟需考虑具有不同性质流体之间的相互作用. Macklin等[22]将体积分数概念引入图形学领域，用来表示不同相在离散空间中的分布.SPH方法已与体积分数方案结合，以模拟多种流体交互运动. Ren等[23]通过计算相间漂移速度实现多相流体交互，随后Yan等[24]将Ren的方法扩展到包含固相的多相流模拟中.根据流体雷诺数特性，含有离散相的混合物需采用不同算法进行模拟. Lin[25]对流动缓慢的多孔介质流进行了模拟. Nielsen和Østerby[26]通过粒子间的黏性项在各相间进行动量交换，实现了对高雷诺数流体混合物的模拟.使用单相流体渲染方法作用于多相流体时，会导致相间交互产生误差缺陷.

采用基于物理的流体模拟技术可视化手术中硅油填充过程，可以辅助医生进行视网膜脱离手术治疗，帮助医生了解硅油表面张力以及其他物理参数，将为治疗孔源性视网膜脱离及其它玻璃体视网膜疾病提供更安全、更有效的新型手术治疗方法.

2 基于物理的视网膜脱离手术硅油填充模拟

为了多视角观测眼球结构和验证视网膜裂孔角度与剥离状态，需对玻璃体切除术后眼内腔中硅油-水两相液体的交互环境进行模拟. 本文采用基于物理的模拟技术，针对不同液体物理特性，对两相液体运动进行建模；分析硅油-水接触面处相间交互状态，对相间受力进行建模；考虑两相流相间的表面张力作用，对液相接触表面曲率进行计算建模. 本章主要内容分为两部分：水-硅油两相液体耦合及临界面精细化模拟和面向眼内腔环境的固液交互模拟.主要流程如图3所示.

图3 孔源性视网膜脱离治疗建模分析流程图Fig.3 Modeling and analysis flow chart of rhegmatogenous retinal detachment

2.1 水-硅油两相液体耦合及临界面精细化模拟

对眼内的水和硅油两相流体填充过程中不同性质流体交互进行分析，需构建基于粒子的离散化模型，用于模拟多相流体整体运动，综合分析多相流物理特性，建立多相流体运动模型.

2.1.1 纳维-斯托克斯方程与SPH离散化

纳维-斯托克斯方程（Navier-Stokes Euqation，N-S）为描述黏性不可压缩流体运动的物理方程，可以被表示为[27]：

在基于SPH方法的流体模拟中，模拟区域内任意物理场A在坐标x处的取值可以根据狄拉克函数性质，被近似为[28]:

其中，W为一归一化类高斯核函数，称为光滑核函数（Smoothing kernel function），其中h被称为支持半径（Supporting radius），本文选取三次样条函数（Cubic spline kernel)[28]作为光滑核函数，为对应的体积积分变量.SPH方法将经典物理学中的连续介质离散化为宏观质点，如图4所示.

图4 SPH方法下粒子数值近似示意图Fig.4 Numerical approximation of particles using smoothed particle hydrodynamics method

由于光滑核函数二阶导数计算会产生高阶数值误差不稳定现象，故的拉普拉斯算子近似可采用人工形式[28]:

而N-S方程中的梯度、拉普拉斯算子可由式（4）和式（5）进行计算，本文基于流体不可压缩特性，通过SPH近似密度与流体静态密度的压缩比，利用压强力抵消数值计算过程流体压缩性，使用隐式不可压缩SPH算法IISPH[16]进行流体模拟.

2.1.2 水-硅油两相液体耦合

为实现多相流体相间交互动态效果，需提取多相流体交互界面位置.本文采用颜色场方法[29]通过设定阈值提取边界，实现对于多相流的界面效果控制.

界面张力的定义如下：

2.1.3 液体表面张力建模

表面张力是展现流体微观特性最常见和最重要的物理性质，为了实现其效果，本文引入内聚力模型[27]，考虑分子间的吸引力和排斥力影响，当它们之间的距离低于阈值时，粒子彼此吸引，当它们之间的距离超过阈值时，粒子彼此排斥，可以表示如下:

式中，B是样条函数，本文参考[30]所使用样条函数，可表示为:

2.2 面向眼内腔环境的固液交互模拟

对眼内腔环境中固液交互进行建模，需建立离散化边界模型，实现基于粒子的稳定交互；根据眼球内壁物质特性，进行眼球内部硅油、水固液交互模拟仿真.

2.2.1 固液双向耦合密度计算

首先对三维眼球模型进行点云采样，本文使用泊松盘采样算法实现对于模型的离散化；将相对均匀离散点赋予刚体边界粒子，用于与流体粒子进行交互，图5所示，使其能够应对不同形状的眼球特征，包括由一层或一排边界粒子组成低维度刚体表达；采用Shepherd核方法对于欠采样与粒子采样非均一区域进行数值加权，缓解SPH数值近似过程中邻居粒子过密或过疏导致的数值不稳定以及颗粒黏附伪影问题.

图5 边界处理示意图Fig.5 Schematic diagram depicting boundary handling

对Shepherd核方法，其机制在于，通过考虑邻居刚体粒子相对位置与数量（b ∈ j），显式地计算点云离散化后的每一刚体边界粒子所代表的刚体体积：

Shepherd核方法对于密集采样的区域，边界粒子的体积会变小；对于稀疏采样的区域，边界粒子的体积会变大. 根据镜像原理，可令边界刚体粒子具有同流体粒子相同静态密度，则边界刚体粒子质量只同其体积与流体粒子差异相关，对于流体粒子i ，刚体粒子的质量可表达为：

2.2.2 边界流体压力

在实践中，从液体施加到刚体某些区域的压力对附近的流体粒子没有运动学影响.因此可将刚体边界粒子施 加到流体粒子的压力[16]记为：

该过程仅使用流体粒子的密度和压力，并遵循牛顿第二定律，使从流体粒子到刚体边界粒子压强力为式（16）取反：

边界力的大小基于流体粒子的压强力，该力随着流体粒子接近眼球刚体边界粒子而增加.由于靠近刚体边界的流体粒子的压力将导致其远离边界，且流体之间的压强作用亦会反映在边界压强上，因此本方法无需使用额外的力或位置校正，便可消除流体黏附伪影，并防止流体粒子渗透出眼球边界.

2.2.3 流体对固体的黏附力

为表现流体对于固体的黏附效应，参考式（11），将流体与边界粒子之间的黏附力计算为[27]：

3 实验

为验证本文方法的有效性，本节首先进行硅油-水两相流交互实验以及表面张力对比实验，随后对眼球硅油填充进行仿真，最后进行硅油填充术医疗场景下注硅油排水的模拟场景.本文仿真算法基于C++编写，并使用OpenMP进行多线程并行计算，利用Eigen作为数学计算工具.可视化方面，使用OpenGL作为流体粒子可视化工具，显示实时仿真效果，并使用Blender进行离线动画渲染.

3.1 水与硅油的交互模拟验证

本节通过实验验证两相流交互模型有效性.图6中为两相流溃坝实验，硅油和水分别为两相流体，两种流体在模拟开始时因重力作用而下落，随后接触并逐渐混合.该过程中由于硅油密度略低于水，会因压强力与界面力的作用而不断上浮；而拥有较大密度的水相则会不断下沉.最终两相流体实现完全分层，并具有明显交界面.

图6 两相流溃坝实验. （a～d）流体运动过程粒子状态；（e～h）渲染后效果Fig.6 Dam break experiment of two-phase flow: (a-d) particle state of fluid motion; (e-h) post-render effect

该实验表现了眼球内部两相流体交互运动过程中受重力与表面张力影响而产生的分层作用，证明了本文方法可有效模拟硅油-水两相流体的动力学交互过程.

3.2 流体表面张力实验

本实验通过令水块跌落入容器并观察液体形态，验证本文表面张力方法有效性.图7展示了水块掉落在水平平板上后形成的溅落水面情况. 水块掉落后先水平平铺，并散至四周，最终静止.

图7 不同张力系数下的水块冲击表现结果. （a～d）第22帧时 α =0,0.1,0.5,0.8的效果；（e～h）第69帧时 α =0,0.1,0.5,0.8的效果；（i～l）第503帧时（静止后）α =0,0.1,0.5,0.8的效果Fig.7 Impact performance of water blocks with different tension coefficients: (a-d) effect at frame 22 when α =0,0.1,0.5,0.8; (e-h) effect at frame 69 when α =0,0.1,0.5,0.8; (i-l) effect at frame 503 (after rest) when α=0,0.1,0.5,0.8

当张力系数设置为0时（无表面张力作用），可以观察到液体呈均匀分布，覆盖绝大部分底板，未出现聚集现象，无高度突起. 在张力系数逐渐增大过程中，以第22帧为例，由于表面张力的作用，液体溅开的范围逐渐变小，液体内聚力同时增大，使得水珠趋向于聚集效果.由第503帧（静止）可以看到，因为系数增大，底板上液体覆盖面积也越来越小，侧面说明了表面张力是影响液体聚集增加高度的主要原因.

表面张力实验不同参数情况下液体表面积总和随时间变化折线如图8所示. 纵轴代表表面积，横轴代表时间帧. 从0～100帧来看，由于液体撞击底板导致表面积数值急剧上升，随后由于液体回滚表面积短时间缩小.在没有表面张力的作用下()，在该时间段表面积数值未出现明显震荡，而在表面张力作用的情况下(α =0.1,0.5,0.8)，表面积数值出现了明显的类似“”形振荡，表明在此期间表面张力在克服其他外力、压力的作用.

图8 不同张力系数下液体表面积变化率Fig.8 Rate of liquid surface area change under different tension coefficients

3.3 眼球内水和硅油两相交互实验

本节进行了眼球内部的水和硅油两相交互模拟实验，如图9所示，在3组实验中白色表示硅油，透明色表示水. 其中图 9（a）和（d）为两相均没有表面张力作用，图9（b）和（e）为只有硅油具有表面张力作用， α =0.85，图 9（c）和（f）为两相均有表面张力作用，其中硅油相表面张力 α =0.85，水相表面张力.可以发现，在没有表面张力的作用下，图9（a）中硅油上浮，但是出现聚集度较低的现象，相比于图 9（b）和（c），白色硅油呈现分散状态.由图9（b）可以发现，在水没有表面张力的作用下，水与硅油出现分界不明显现象.在两相均设置适当表面张力的情况下，图9（c）中两相界线明显，硅油呈现稳定的聚集状态，且有突起现象；同时，下方的水的聚集情况更加合理.

图9 眼球内两相液体交互. （a～c）第266帧两相均无表面张力、只有硅油具有表面张力、两相均有表面张力时的交互情况；（d～f）第376帧两相均无表面张力、只有硅油具有表面张力、两相均有表面张力时的交互情况Fig.9 Two-phase liquid interaction in the eyeball: (a-c) interaction effect of the two phases without surface tension, only the surface tension of silicone oil, and the surface tension of both phases in frame 266; (d-f) interaction effect of the two phases without surface tension, only the surface tension of silicone oil, and the surface tension of both phases in frame 376

3.4 眼内腔注入硅油流程模拟

为验证文中方法对于硅油填充术流程模拟的有效性，开展手术中硅油填充过程模拟实验，如图10所示，展示了眼内腔中注入硅油（白色）与排出水分（透明液体）的整个流程.在眼球内部完成玻璃体切割后，眼内腔充满水分，之后需要通过手术中对眼球置入的导管在注入硅油的同时排出较重的水.

图10 硅油填充手术流程模拟：硅油顺导管注入并排出眼内水的流程Fig.10 Simulation of silicone oil tamponade: flow of silicone oil injection and discharge of water from the hole along the guide tube

实验场景设置上，共向眼球内部置入两根导管，其中一根用于注入硅油，另一根用于排出水分，基于重力与大气压强以及相间作用下共同实现排水注油效果.

在起始状态中，眼球内部只存在透明的水分，插入较深导管可排除一定水分平衡内外压强.随后通过较浅导管注入硅油.由图可知，文章方法可有效模拟两导管形成的连通器效应，在诸如硅油过程中顺利将水排出，并全流程保持硅油在上，水在下.此外，由于文中表面张力方法作用，两流相间交界面曲率亦十分明显，可以真实地还原注入过程中两相流体交互状态.

4 结论

通过玻璃体切除联合硅油填充手术治疗孔源性视网膜脱离已经非常普遍，而硅油填充虽然可以用于贴附视网膜裂孔，但也会带来视网膜损伤和青光眼等并发症，所以在满足有效修复视网膜脱离的前提下，尽可能减少眼内硅油填充使用量，可有效提升治疗效果.但确切硅油填充量难以预测，术中不易掌握眼球内填充状态，为此本文采用基于物理的流体模拟方法对眼球内硅油填充过程进行模拟.根据实际医学治疗过程建立眼球三维模型并进行粒子采样，基于水与硅油的不同物理性质对水-硅油两相流动及交互进行建模，并构建固液交互模型，实现两相液体（水和硅油）与固体边界（眼内腔）的交互仿真.实验结果表明，本文方法能够较好地呈现眼球内多相运动、表面张力、边界处理等效果.

本文为了解硅油填充后的眼内状态提供了一种有效的方式，也给物理模拟辅助眼科学带来了新的研究思路.但是，由于眼球内结构的复杂性，物理模拟结果直接应用于孔源性视网膜脱离手术还存在一定的挑战，例如如何进行更加精细化的建模，如何根据模拟结果准确估测硅油填充量等，是下一步研究的重点与难点.