自适应网格在水平集方法中的优化与应用

任宇佳， 马朝青,2

(1.烟台大学 计算机与控制工程学院 山东 烟台 264010；2.高端海洋工程装备智能技术烟台市重点实验室 山东 烟台 264010)

0 引言

水平集方法(level set method，LSM)是一种用于界面追踪和形状建模的数值技术。水平集方法的优点是可以在笛卡尔网格(Cartesian grid)上基于欧拉法(Eulerian approach)直接对演化中的曲线、曲面进行描述和数值计算,不必对曲线、曲面进行参数化[1-2]。水平集方法的另一个优点是可以方便地追踪物体的拓扑结构的改变。在水平集方法应用中，为了在控制算法时间复杂度的同时保证轮廓的分辨率和水平集方程的求解精度，通常在算法中引入自适应网格[3-4]。

自适应网格方法(adaptive mesh refinement，AMR)是求解偏微分方程的工具。AMR的核心思想是：在迭代过程中，当某些区域变化剧烈(例如大变形、激波面、触间断面和滑移面等)时，该区域网格自动实现分级细化。因为AMR网格生成过程简单快速、自动化程度高、易于实现，在众多领域被广泛应用[5-6]。

但是采用AMR的水平集方法在实现过程中仍存在难点，由于AMR是一种分级细化的不均匀网格，在不同细化级别之间的节点(T节点)处会出现邻居节点集合不完整的情况，影响偏微分方程的近似求解。为了解决这一问题，近几年AMR相关研究中提出了最长边二分法[7-8]、最新顶点二分法[9-10]和red-green细化[11-13]等方法。使用最长二分法细化后不能保证子网格的形状规律性。尽管在文献[14-15]中证明了二维三角网格经过最长边二分法后得到的子网格的形状规律性，但在三维中仍是一个不确定问题。后两种则可以保持初始网格中某一元素的所有子代的形状规律性，相较于red-green细化,最新顶点二分法属于复杂度较高的网格嵌套细化。事实上，red-green细化是对带有T节点的1-irregular网格常规细化的改进，其目的是实现细化网格的一致性。该细化规则首先在文献[11-12]提出，并在文献[16-18]中成功地扩展到高维的单纯网格。

水平集方程的求解是一个多次迭代求解的过程。引入AMR的水平集方法在结束每个时间步的计算后，需要重新确定细化和非细化网格，如果每次都进行网格的全局遍历会增加复杂度。针对这一问题，文献[19]提出了一种快速判断方法，即将遍历区域缩小到每个时间步计算所得曲面的临近网格，有效地减少了界面演化时的计算工作量。文献[20]开发了一种定位水平集，不需要在空间域中显式地找到界面位置的快速方法。该方法解决了如何将只在界面上给出的量扩展到界面邻域的重要问题。文献[21]将窄带法和快速行进近似法(the fast marching approximations)相结合，提出了一种用于研究生长过程中界面演化的算法，在重要的物理边界上建立了一个最高分辨率的可控有限宽度区域。

本文为了提高AMR水平集方法的迭代效率，在采用了局部水平集方法的基础上结合了窄带法[21]对需要标记的区域进行选择，优化了细化区域的选择机制。同时，采用基于red-green细化的改进算法优化了T节点的近似计算。最后，为验证所提出算法的效果，将该算法应用于基于水平集的血栓生长模型中，优化了模型中水平集求解部分，并利用该模型进行了血栓形成过程的仿真实验。

1 局部水平集方法基本原理

水平集方法[22]在解决轮廓演化问题时，把低维的计算问题上升到高一维，即把N维的轮廓利用N+1维的水平集函数进行描述，获得一个闭合的超曲面，水平集函数值为0的零水平集即为被描述的轮廓。这样，轮廓的运动就可以通过闭合超曲面的演化方程得到，仅需要确定演化过程中零水平集的位置就可对给定轮廓进行追踪。

局部水平集的基本思想是在零水平集周围设定一条覆盖零水平集的带状区域，区域内的水平集符合|φ|

在局部水平集方法[23]中，局部水平集方程被定义为

(1)

其中：φ(x,t)是局部水平集函数；截断函数C(φ)为

(2)

水平集的初始化仅在T0={x:|φ(x)<δ|}处进行。

2 Red-green混合细化方法

本节将对现有的red-green及改进后的red-green分别进行细化描述。

2.1 传统的red-green细化

Red-green细化最初由文献[11]提出，是一种基于三角网格的局部细化方法。它的实质是用两种不同的策略对单元格进行细化，从而减少T节点。随着三维网格的发展，这种red-green技术也适用于三维四面体网格[16-18]。二维三角网格分割的red细化策略和green细化策略如图1所示。从图1(a)可知二维的red细化(也称为正则细化)是通过连接边的三个中点将一个三角形分成四个子三角形(称为红色元素)。从图1(b)可知green细化(也称为不规则求精)是连接顶点和与之相对的边的中点，通过green细化分割获得的两个子三角形称为绿色元素，两条分割边P3P1和P1P4称为绿色边(或不规则边)。元素T=T1∪T2∪T3∪T4的red细化为四个全等子域T1，T2，…，T4；元素T=T1∪T2的green细化为两个子三角形T1和T2。

图1 传统的red-green细化示意图Figure 1 The diagram of traditional red-green refinement

Red细化策略的优势在于，可以保证细化得到的子网格的形状与原始网格比较接近，但是会产生较多的T节点，在计算时增加难度；Green细化策略可以减少T节点的产生，但是不能保证子网格的形状。因此，需要将两种策略结合使用。

2.2 无T节点四边形网格的局部细化

原始的四边形网格细化一直使用四等分作为red细化策略。在实际应用中，如果每个边上最多允许一个T节点存在，red细化的对分完全可以满足四边形网格上的局部细化。但是在实际的计算中，T节点数量的增加会导致计算难度的升高。在数值算法中，为了保持离散解的连续性，必须对T节点施加约束。为了尽量减少求解过程中对T节点特殊处理，Bank等[11]给出了green细化策略，它们直接将T节点与四边形顶点相连。由于T节点出现在边的中间，在将T节点与不相邻顶点相连的情况下，即使初始网格为四边形，三角形仍将出现在细化网格中，生成的三角形和四边形的混合网格需要较为复杂的数据结构进行描述。因此，将red细化和green细化结合后可以形成新的混合细化，首先构造合适的red细化策略，保证细化网格与原始网格的一致性，然后辅助green细化策略进一步减少T节点的数量。同时，将混合细化方法扩展到四边形以消除T节点。

如图2所示，在二维混合细化算法中，原始网格首先通过三等分单元格K的边，将其分割为9个子单元格K1,K2,…,K9进行red细化。内部节点P5、P6、P7、P8相对原始单元格顶点P1、P2、P3、P4的位置如下：

图2 Red细化Figure 2 Red refinement

P5=(4P1+2P2+P3+2P4)/9；

(3)

P6=(2P1+P2+2P3+4P4)/9；

(4)

P7=(2P1+4P2+2P3+P4)/9；

(5)

P8=(P1+2P2+4P3+2P4)/9，

(6)

其中Pi=(xi,yi)，i=1,2,…,8为坐标。

如图3所示，处于细化边缘的非细化单元格根据邻接网格的细化情况可以分为四种模式：三侧细化；相对两侧细化；相邻两侧细化和单侧细化，分别对应产生6个、4个、4个、2个T节点(实心圆)。消除这些T节点的方案之一是该非细化网格进行同样的red细化(虚线)，但会产生新的T节点，数量分别是2个、4个、4个、6个(空心圆)。该方案仅对三侧细化模式具有减少T节点的效果。根据文献[11]提出的“Ke-1 Neighbor规则”，后三种情况需要采用green细化(实线)来消除T节点。

图3 边缘单元格的不同细化情况及相应的green细化模式Figure 3 The different adjacent refinement of edge cells and the corresponding green refinement mode

2.3 水平集方法中的笛卡尔细化网格

(7)

该方法求解水平集方程φi时需要使用{φi-3,φi-2,φi-1,φi,φi+1,φi+2,φi+3}。根据求解方法的精度不同，求解过程中使用的相邻节点水平集数量也不同。对于非均匀细化笛卡尔网格，两级细化分界处的单元格的T节点处，由于相邻数据不连续，需要通过与其相连节点的φ值求得近似相邻数据，用于水平集方程求解。

为了消除T节点带来的附加计算，我们将改进后的red-green细化与窄带法相结合。首先对窄带区域内部的节点进行red细化，并对靠近窄带边界的节点进行green细化，以消除T节点。

如图4所示，深色区域为red细化网格，浅色为green细化网格，判断原始网格的“分类”：

1)与red细化相邻的原始单元格称为红色元素(图4(a))；

2)与green细化相邻的原始单元格称为绿色元素(图4(b))；

3)同时与red细化网格和green细化网格相邻的原始网格，也将其看作红元素(图4(c))。

图4中red细化网格和红色元素为深色，green细化网格和绿色元素为浅色。然后，忽略绿色元素，仅针对每一个红色元素进行如下操作：

图4 网格细化过程示例图Figure 4 Example diagram of mesh refinement process

1)统计红色元素边上的T节点数n；

2)将红色元素进行“预red细化”，统计细化后新增T节点数m；

3)若m/n<1，对该红色元素实施red细化，否则实施相应的green细化。

3 基于水平集的血栓形成计算模型

血栓[25]是血液循环中血管表面或心脏内壁上形成的凝块或沉积物。当血管内膜受损时，内皮下基质暴露，损伤附近的血小板被激活。随着血小板的粘附和聚集，形成以血小板为主要成分的初期血栓。根据血栓形成的原理，我们设计了初期血栓的混合模型，如图5所示。血栓生长模型包括血管、破损的血管内膜、血小板、血浆和初期血栓，在设计成管状的血管中，血浆自左向右流动，带着血小板，经过破损的血管内膜这一区域的血小板被凝血因子激活，然后成为初期血栓的一部分。

图5 血栓生长模型示意图Figure 5 Schematic diagram of thrombus growth model

血栓生长计算模型包括血浆、血小板和血栓三部分。如图6所示三个子模型的结构和关系：血浆子模型提供血流速度场并输送血小板；血小板模型是用来模拟血小板粘附和聚集的过程；血小板通过血栓子模型合并到血栓中。我们将初期血栓的表面定义为零水平集，并通过水平集方法追踪界面位置的变化，窄带的宽度值与血小板的直径成正比，确保新粘附上的血小板可以被包含在窄带内。

图6 血栓生长计算模型结构图Figure 6 Structure diagram of thrombus growth calculation model

血浆子模型是用来模拟血液流动的子模型，提供了血液流速，并且负责运输血小板。我们假设血浆是不可压缩的，通过求解Navier-Stokes方程获得血流速度。血小板是初期血栓的主要组成部分。血小板模型用于模拟血小板的激活、粘附和聚集反应。分别计算血小板在血流冲击力、破损内膜粘附力和初期血栓的聚集力作用下的移动速度，并使血小板在该速度下移动。

模拟血栓生长过程的子模型中，我们将血栓的表面定义为由零水平集描述的自由移动界面，水平集φ由符号距离公式计算得出：血栓内部φ<0；血栓外部φ>0；血栓表面φ=0。将破损表面定义为初期血栓的初始表面。在血栓生长过程中，发生粘附和聚集的血小板被定义为组成初期血栓的一部分，因此，当血小板模型中判定某一血小板发生粘附或聚集后，在该血小板所在位置定义界面移动速度，使血栓表面向外扩张产生扩张的表面通过水平集方法追踪界面位置的变化。在追踪过程中，定义窄带的宽度值不小于血小板的直径，以确保新捕获的血小板完全包含在窄带内。

4 实验结果与分析

为了验证改进的red-green细化策略在局部水平集中的有效性，我们建立了不同分辨率下的轮廓和窄带，并对细化数据进行统计并加以分析。在边长为5 mm的正方形网格中，建立半径为15 mm的圆形轮廓，窄带宽度为0.2 mm。网格的分辨率从25*25开始，以25为间隔递增至200*200。网格中red细化单元格数、red细化后T节点数和追加green细化单元格数的变化曲线如图7所示。随着分辨率的提高，被细化的网格数增长幅度非常大。由于圆形轮廓较为简单，T节点数目虽然增长，但幅度相对较小。消除这些red细化出现的T节点所追加的green细化单元格数目的增长幅度小于T节点的增长幅度，说明green细化在有效消除T节点时并不会带来更多的计算量。

图7 不同分辨率下细化单元格数及T节点数变化曲线Figure 7 Variation curve of refined cell number and T-node number under different resolutions

同时，针对静脉，我们使用第3节介绍的计算模型分别在无狭窄的直血管和有部分平滑狭窄的血管中对血栓的形成过程进行了仿真实验，各项血管参数分别为：血流速度1 m/s；血液黏度3.5×10-4Pa·s；血管直径4 mm；血管长度12 mm；破损区域面积1.57 mm×2 mm；狭窄程度25%；狭窄长度2 mm。并利用视觉化工具函式库(Visualization Toolkit，VTK)将仿真实验的结果进行可视化，方便观察血栓的形成过程，图8为两种血管的形态及破损血管内壁的位置的可视化图像。破损内壁初始为光滑表面，随着血小板被破损血管壁捕获成为血栓，表面不断扩展，形成增多并叠加的颗粒状突起。

图8 仿真血管Figure 8 Simulated blood vessels

在无狭窄血管和有狭窄血管中将血流速度都分别设置为三种情况：2 mm/s，4 m/s和6 m/s(对应雷诺数Re分别为0.01、0.02和0.03)，分别进行5次时长为50 ms的血栓生长仿真实验，将平均结果绘制生长曲线，如图9所示。狭窄血管中的生长速度总体大于无狭窄血管，两种血管中的血栓生长均随血流速度增加变快。将该结果与Kamada等的静脉血栓仿真[26]结果对比，血栓生长速率以及该速率与雷诺数的关系是一致的。不同之处在于，Kamada实验中考虑了血栓脱落，生长曲线存在中断，讨论脱落栓子阻塞血管的危险性，本文期望观察血栓在出现位置处的阻塞情况，未考虑血栓脱落，因此血栓在整个仿真时间内持续上升。

图9 不同血流速度下血栓生长曲线图Figure 9 Thrombus growth curve under different blood flow velocity

根据仿真过程中网格数据的比较，如图10所示，无论是否存在狭窄，如果在目标区域中使用均匀细化，即在覆盖全部血栓的区域建立六面体网格，细化网格的数量和增长率远高于基于窄带方法的局部细化。使用局部细化后，局部细化网格中仍存在约10%的T节点，求解演化方程时需要分别计算。引入green细化之后，尽管细化网格的总数增加了约3%，但green细化网格的数量仍小于T节点的数量。同时，在统计增长率中，发现尽管在狭窄血管的模拟中细化网格的总数大于非狭窄血管的数目，但其增长率小于非狭窄血管的增长率。这表明由于血管几何结构的复杂性，细化网格的数量仅在狭窄初期才大，但是在血栓生长过程中细化网格的数量却缓慢增加。由于窄带方法的局部细化是基于零水平集界面，该模型中的零水平集是血栓表面，细化网格也是血栓表面的增长。因此，在狭窄的情况下，细化网格的缓慢生长速度表明由于血管狭窄所导致的血小板颗粒的粘附和聚集范围更加集中。

图10 仿真实验中随血栓生长的网格数据比较Figure 10 Comparison of grid data with thrombus growth in simulation experiment

5 结论

本研究中，我们将改进后的red-green细化策略融合到局部水平集中，有效地提高了笛卡尔网格的运算效率。一方面，red-green细化的使用消除了T节点，避免了由于找不到邻居节点导致的计算难度的增加，提高了计算效率；另一方面，窄带法的应用，使得每个时间步后，只需要对窄带内的点进行是否需要细化或粗化的判断，定义在零水平集上的速度函数也只需要扩展到窄带内的部分，降低水平集方法的计算量和存储空间。我们先是使用简单的圆形轮廓对其进行实验，又将该方法应用与血栓生长过程的建模和仿真中。实验结果表明改进后的算法在网格分辨率增加时能够更好地减少计算量，并能够有效地提高计算效率。该方法对于血栓生长过程的建模与仿真有利于进一步分析血栓形成的相关影响因素。本文研究中，在选择细化区域时仅考虑了与零水平集的距离，未充分考虑轮廓形态的复杂程度。我们下一阶段的工作将继续针对自适应细化网格进行研究，将能够代表轮廓几何特征的曲率等元素加入判断细化区域的特征参数中。