基于GIS的高质量约束Delaunay三角网格剖分

第26卷 第5期2010年9月地理与地理信息科学

Geog ra phy and Geo-Infor matio n Science V ol.26 N o.5September 2010

收稿日期:2010-04-21; 修订日期:2010-06-10

基金项目:国家自然科学基金重点项目(50839001);国家自然科学基金项目(50874021、50779006);辽宁省高等学校科研项目计划(L20100321) 作者简介:赵晓东(1969-),男,博士,教授,从事GIS 在采矿、岩土工程应用和水动模型耦合的研究。E-mail:xdong.zhao@

基于GIS 的高质量约束Delaunay 三角网格剖分

赵晓东1

,晏小宝1

,沈永明2

,王 亮

2

(1.大连大学院士创业园中日地层环境科学研究中心,辽宁大连116622;2.大连理工大学海岸和近海工程国家重点实验室,辽宁大连116023)

摘要:在分析现有非结构化网格剖分算法的基础上,提出了一种GIS 支持下的改进分治算法实现约束Delaunay 三角网格剖分。该方法利用了GIS 的空间拓扑关系对算法输入数据进行预处理,基于三角形的统一数据结构实现了网格细化,对输出剖分网格进行准确的拓扑和约束条件的检查,并基于推进阵面算法思想,结合空间邻近拓扑关系实现了三角剖分节点和网格的重新编号,方便了实际问题中开边界条件的赋值,提高了计算效率。实例应用表明,该方法大大简化了数值模型非结构化网格剖分的前处理过程,集成了几种综合算法的优点,在保证原分治算法时间复杂度的基础上,提高了约束条件下Delauna y 三角网格生成的质量。关键词:网格剖分;GIS;约束Delaunay 三角剖分

中图分类号:P208 文献标识码:A 文章编号:1672-0504(2010)05-0024-05

0 引言

在应用有限元、有限差分和有限体积法对力学问题进行数值计算的前处理中,网格自动剖分占有重要的位置。目前,实现网格剖分的算法很多,而且不断有新的算法推出,其相关研究领域不仅针对数值模拟计算,对GIS 数据表达、地学分析、计算机视觉、表面目标重构等众多领域也是一项重要的应用技术[1-3]。网格剖分可分为结构化网格(Structured Gr id)和非结构化网格(Unstructured Grid)两种。非结构化网格易于控制网格单元大小、形状及网格点位置,可以实现合理分布网格的密度,提高计算精度,因此具有比结构化网格更大的灵活性和对复杂边界更强的适应性。生成二维非结构化网格的常用方法[4]有四叉树法(Quadtr ee)、Delaunay 三角剖分法(Delaunay T riangulation,DT)和推进阵面法(Ad -v ancing Fr ont M ethod,AFM)以及几种方法的综合和改进。由于Delaunay 三角网是Vor ono i 图的直线对偶图,具有空外接圆和最大最小角特性,还可以尽可能地避免病态三角形的出现,实现约束条件下的三角剖分,因此,在数值计算非结构化网格剖分和GIS 三角网(T IN)的数据格式表达中广泛使用。

目前常见的构建Delaunay 三角网的算法有:分治算法(Div ide -and -conquer )、逐点插入算法(Incre -m ental Insertio n)、生长算法和扫描线算法(Sw eep -line)[3,5-7],其中以分治算法效率最高,时间复杂度

为O(N logN)[5]。由于实际的不同需要,除几何约束条件外,对三角剖分的角度、大小和节点物理量的控制也都有相应约束,目前的剖分算法尚不能满足以上全部约束条件进行剖分。本文基于改进的分治

算法,提供了梯度变化和针对实际问题的单元和节点重编号,并利用GIS 强大的空间数据处理和分析功能对算法做数据前后处理,保证了高质量约束Delaunay 三角网格剖分的生成。

1 约束Delaunay 三角化方法

分治算法[5]的基本思路是把输入点集分割为数个较小的点集,在各个子点集内生成小三角网,再逐级合并的一个递归过程。子点集最终只有两点或者三点,形成一条边、共线边或者三角形。算法必须将点集连接为凸区域(凸壳),对凹区域和多连通区域会产生域外三角形,无法保证约束边的存在。约束Delaunay 三角化的分治算法的关键过程及实现方法如下。

1.1 点集划分

原分治算法只有Y 轴方向的分割,这里改进为交替分割[8],即点域X 轴方向的长度大于Y 轴方向的长度,则以X 轴方向对半分割点集,否则以Y 轴方向对半分割点集。该过程为递归调用的分割过程。

1.2 约束特征线重构

在Delaunay 三角剖分上,采用逐次加入约束特

征线的方法进行约束重构,可以证明该重构方法是收敛的[4],具体步骤如下:1)对于任一条约束直线,查找与之相交的三角形,即该约束线穿越的所有三角形组成的影响域;2)删除该影响域内的所有三角形,形成一个多边形空腔,称为恢复域;3)以该约束直线为界将恢复域一分为二,形成两个多边形子域;4)对这两个多边形子域进行Delaunay 三角剖分,得到包含约束直线的影响域内的三角剖分;5)重复上述步骤,直至所有的约束线都被加入。1.3 基于三角形的统一数据结构

Delaunay 三角剖分涉及三角形和线段两类几何数据,如果将线段表示为退化的三角形势必增加算法特殊情况判断的冗余。如图1a 所示,利用/影子0层表示完整的三角形结构,凸壳上的每一线段有一个顶点在无穷远处(算法中采用空指针)的/影子0三角形,特征线段有两个/影子0三角形。图1b 中下端虚线三角形为子网合并中新建的/影子0三角形,阴影的三角形采用局部LOP 优化算法[7],保证三角形的Delaunay

特性。

图1 分治算法的三角形数据结构和子网合并[9]

Fig.1 Data s tructure of triangle and halfway through a merge

step of the dividing and conqu ering algorithm

1.4 网格细化修正

剖分细化是在网格中插入顶点细化网格,使得

网格最大和最小角满足约束。本文采用Shew chuk 的细化算法[9,10],它是对Ruppert 和Chew 提出算法[11,12]的综合改进算法,顶点插入使用LOP 算法[7],保证Delaunay 特性,遵循两个规则:

规则1:一条线段的直径圆是唯一且最小的包含该线段的圆。如图2所示,如果有一点位于直径园内,则称该线段被侵入。任何被侵入线段在其中点

位置插入顶点进行分割。

图2 线段被递归细分直至不被侵入为止[9]

Fig.2 Segments are split recursively until no

segments are encroached

规则2:一个三角形的外接圆是唯一经过三角形3个顶点的圆。如果三角形的一个角度过小或面积过大的情况下满足了约束条件,则称该三角形为坏三角形。如图3所示,以外心顶点(外接圆的中心)作为插入点进行分割。

图3 每个坏三角形线在外心顶点插入分割[9]

Fig.3 Each bad triangle is split by inserting a vertex

at its circumcenter

剖分细化具体算法如下:1)取约束线段集合中的每一条约束线段,检查它是否被其他点侵入,是则转第2步;否则转第3步。2)按规则1将约束线段分割,将新的约束线段加入约束线段集合,转第1步。3)取剖分三角形集合的三角形进行检查,全部合格,转第5步;否则,对于不满足约束条件的三角形按规则2插入点。4)如果插入点不侵入约束线段,将该点插入网格中,转第3步;否则,转第2步。5)算法收敛停止。

2 GIS 支持下约束Delaunay 三角网格剖分

高质量Delaunay 三角网格意味着必须满足一定的约束条件,而约束条件是根据工程需要或者数值模拟的条件给定的,如流体动力学问题一般在边

界、结构物和扩散源附近采取网格加密方法,这些都要对输入点集进行预处理;输出的剖分三角形内角不能太大也不能太小,面积梯度变化适中。因此,要实现本文提出的非结构网格生成算法,需要提取计算区域边界点的坐标,包括外边界和内边界上各节点光滑、间距的处理、内部约束特征线的设置等。如图4所示,这些预处理都可以在GIS 的支持下完成,同时可以实现空间数据的管理,保证数据的拓扑准确性和剖分三角网格的正确性。2.1 边界点的输入

剖分算法输入的为点集的坐标值,非实际的拓扑数据,而且边界点组织直接影响到剖分输出的结果,需要在输入前进行必要的预处理。如图4所示,按照拓扑关系,将输入数据划分为Poly gon 、Poly line 和Point 3种矢量数据,其涵盖了流体动力学问题中的海岸边界、岛屿、开边界和水深等高线等所有输入

数据。其中,海岸边界采用Do ug las -Peucker 算法做简化处理后,根据实际计算需要,对边界和岛屿边界

页

25第第5期

赵晓东等:基于G IS 的高质量约束Delaunay 三角网格剖分