摘 要：在常规逐点插入算法的基础上，提出了一种改进的逐点插入构建Delaunay三角网的算法。引入散乱点集有序化、三角形单元分类的方法快速生成Delaunay三角网。
关键词：Delaunay三角网 逐点插入 地震层位模型 OpenGL

在三维地理信息系统中，地震层位模型的可视化研究有着广泛的应用^[1]。如在油藏勘探领域，通过重构地表特征层位、判断其走向可预测储层的位置。地震层位模型是构建地表层位表面空间位置与其相关属性信息的数字化表示，它是对地表层位在地震波采样数据基础上的表面重构^[2]。
由地震波采样得到并经过速度场解释之后的地震数据表现为三维空间的散乱数据点集。针对该类数据的特点，普遍采用基于三角网的建模方法构造层位模型。其中，Delaunay三角网具有良好的形态，在表达地质形态方面表现较为出色。
如何快速、高效地构建Delaunay三角网一直是众多学者研究和关注的焦点。迄今为止出现了不少成熟的算法，基本算法为逐点插入法、分割-合并法及三角网生长法等^[3]。其中三角网生长算法由于算法效率低，目前较少采用；分割-合并算法最为高效，但相对复杂，且深度递归，对内存要求高；逐点插入法实现简单、占用内存小，但其时间复杂度较高，效率低于分割-合并算法。
在油藏勘探过程中，有时须补测部分数据，这会引起地震数据的动态变化。因此需要对已存在的Delaunay三角网进行局部更新，按需插入和删除某些点。这就要求有一种能够快速构建三角网的方法。
考虑到工程应用的需要，本文在常规逐点插入算法的基础上，引入散乱点集有序化、三角形单元分类的方法快速构建Delaunay三角网，有效地解决了地震层位建模的工程问题。实测结果表明，改进后算法的时间复杂度大为降低，从而提高了算法的效率。
1 逐点插入算法流程
假设给定散乱数据点集P={V_i|i=1，……，n}，其中V_i的三维坐标为(x_i，y_i，z_i)，n为点的个数。经过Delaunay三角化之后，输出Delaunay三角网D={T_j|j=1，……，m}，其中Tj表示第j个三角形，m为三角形的个数。逐点插入法的基本思想为：在某一已存在的三角网中插入一点，遍历所有三角形，查找外接圆包含该点的三角形集合，然后用LOP(Local Optimization Procedure)法则优化^[4]，确保所建三角网符合Delaunay法则。以下用伪码描述逐点插入算法的流程。
SUBROUTINE：Delaunay Triangulate
输入：散乱数据点集P={V_i|i=1，……，n}
输出：Delaunay三角网D={T_j|j=1，……，m}
初始化三角形数组
确定包围三角形，把包围三角形的顶点添加到顶点数组的末尾
把包围三角形添加到三角形数组
FOR i=1 to n DO
对于顶点数组中的顶点Vi，初始化边缓冲区
FOR j=1 to m DO
对于三角形数组中的三角形Tj，计算其外接圆圆心和半径
IF Vi位于外接圆内 THEN
把Tj的三条边添加到边缓冲区，从三角形数组中删除Tj
ENDIF
ENDFOR
删除边缓冲区中所有重复指定的边，只保留闭合多边形的边
把由Vi和闭合多边形的边形成的所有新三角形添加到三角形数组
ENDFOR
从三角形数组中删除所有包含包围三角形顶点的三角形
从三角形数组中删除包围三角形
END
2 算法改进
2.1 自定义数据结构
合理的数据结构可以提高算法的执行效率。本文设计了点、边和三角形三种数据类型来存储空间数据点集，表达三类元素之间的拓扑关系。
(1)点数据类型用来存储三角网中所有的数据点坐标。
typedef struct
{
float x，y，z； //空间散乱点坐标
} POINT3D；
(2)边数据类型存储组成边的二个顶点的索引号，记录三角网中的边与点之间的拓扑关系。
typedef struct
{
long p1，p2； //两个顶点的索引号
} EDGE；
(3)三角形数据类型存储组成三角形的三个顶点的索引号，记录三角网中三角形与点、边之间的拓扑关系。
typedef struct
{
long p1，p2，p3； //三个顶点的索引号
} TRIANGLE；
可见，空间数据点集坐标存储在顶点数组中，边缓冲区和三角形数组存储对应点的索引号，这样既节省存储空间，又便于数据的直接访问。而且由于数组元素的随机访问速度最快，可提高搜索效率。
2.2 散乱点集有序化
常规逐点插入算法依次从存储散乱数据点集的数组中提取一点，然后遍历当前三角形数组中所有的三角形，判断其外接圆是否包含当前待插入点。这一步是点插入过程中最费时的一步，其计算效率取决于三角网的数据结构和搜索方法。本文在对散乱数据点集进行Delaunay三角化之前，首先对其预处理，使散乱点集有序化。具体处理过程如下：(1)遍历原始散乱数据点集，计算其X轴坐标范围(X_min，X_max)和Y轴坐标范围(Y_min，Y_max)，通过比较得出坐标分量变化较大的坐标轴；(2)以此坐标轴的坐标分量为比较依据，按照从小到大的顺序对散乱数据点集进行排序，得到有序的空间数据点集；(3)将排序之后的点集投影在XOY平面上，删除X、Y坐标分量相等的重复点，去除奇异性，保证惟一性。
2.3 三角形单元分类
每次在已生成的三角网中插入一点之后，一些不符合Delaunay准则的三角形被删除，同时又生成一批新的三角形。常规逐点插入算法对这些三角形不加以区分，每次插入操作都必须遍历所有三角形，计算其外接圆圆心和半径，判断插入点是否位于其外接圆内。而实际上，大部分三角形是不需要进行查找、判断的，这些冗余操作大大降低了算法执行效率，也是导致常规逐点插入算法时间复杂度高的主要原因。本文针对这一瓶颈，对逐点插入过程中生成的三角形单元分类为已完成三角化和未完成三角化二类并加以标识，从而只需要对那些标识为未完成的三角形进行操作，算法效率大为提高。
详细的优化步骤：(1)每产生一个新三角形，初始化为未完成三角化并赋以标识；(2)在已形成的三角网中插入一点时，循环判断当前三角形数组中的所有三角形的分类标识。若标识为已完成则跳过该三角形继续判断下一个三角形；否则进一步计算其外接圆圆心与半径，运用LOP法则进行优化，同时根据三角形单元分类规则更新该三角形的分类标识。
三角形单元分类规则是建立在散乱点集有序化基础上的。假设散乱点集有序化过程中以X轴的坐标分量为比较依据，则定义分类规则如下：
定义 已知插入点V_i(x_i，y_i，z_i)(i=1，2，3，……，n，n为点的个数)和待分类三角形T_j(p1_j，p2_j，p3_j)(j=1，2，3，……，m，m为三角形的个数)，其中T_j的外接圆圆心为C_j(cx_j，cy_j，cz_j)、半径为r_j。若x_i-cx_j≥r_j，则三角形Tj标识为已完成三角化，否则标识为未完成三角化。
分类规则的直观解释如图1所示。若T_j外接圆圆心C_j到经过插入点V_i且垂直于X轴的直线MN的距离大于其半径r_j的长度，则可保证所有X坐标分量不小于插入点V_i的点都位于T_j外接圆外。又因为散乱点集有序化之后按照X坐标分量以升序排列，故对于点V_i之后的插入点V_k(k> i)而言，永远都不会位于三角形T_j的外接圆内，从而T_j可标识为已完成三角化。

5 结 论
本文在仔细研究常规逐点插入算法的基础上，对其进行了优化，使得改进之后的算法不但继承了原算法流程清晰、实现简单、占用内存较小、可局部动态更新三角网的优点，而且提高了执行效率。该算法运用于实际工程中构建地震层位模型，取得了良好的效果。
当前，包括地震层位建模在内的三维地理信息系统正处于从研究走向实用的过渡阶段。在实际应用中，有时还存在着约束三角剖分，也就是部分散乱点之间常常存在着某种约束关系，如地表模型中的断裂线等。在对此种数据集进行三角网生成时，生成的三角网应保持其原有的约束关系，即约束数据集下的三角剖分。因此，如何在约束条件下快速准确地生成Delaunay三角网，是下一步研究工作的重点。
参考文献
1 包世泰，夏斌，崔学军等.地质信息模型研究及其应用.大地构造与成矿学，2004；28(4)
2 胡金星，潘懋，马照亭等.高效构建Delaunay三角网数字地形模型算法研究.北京大学学报(自学科学版)，2003；39(5)
3 武晓波，王世新，肖春生.Delaunay三角网的生成算法研究.测绘学报，1999；28(1)
4 Lee D T，Schachter B J.Two Algorithms for Constructing a Delaunay Triangulation.Int J.of Computer and Information Sciences，1980；9(3)