快速Delaunay三角剖分及其在地震层位建模中的应用
摘 要：在常规逐点插入算法的基础上，提出了一种改进的逐点插入构建Delaunay三角网的算法。引入散乱点集有序化、三角形单元分类的方法快速生成Delaunay三角网。
关键词：Delaunay三角网 逐点插入 地震层位模型 OpenGL
在三维地理信息系统中，地震层位模型的可视化研究有着广泛的应用[1]。如在油藏勘探领域，通过重构地表特征层位、判断其走向可预测储层的位置。地震层位模型是构建地表层位表面空间位置与其相关属性信息的数字化表示，它是对地表层位在地震波采样数据基础上的表面重构[2]。
由地震波采样得到并经过速度场解释之后的地震数据表现为三维空间的散乱数据点集。针对该类数据的特点，普遍采用基于三角网的建模方法构造层位模型。其中，Delaunay三角网具有良好的形态，在表达地质形态方面表现较为出色。
如何快速、高效地构建Delaunay三角网一直是众多学者研究和关注的焦点。迄今为止出现了不少成熟的算法，基本算法为逐点插入法、分割-合并法及三角网生长法等[3]。其中三角网生长算法由于算法效率低，目前较少采用；分割-合并算法最为高效，但相对复杂，且深度递归，对内存要求高；逐点插入法实现简单、占用内存小，但其时间复杂度较高，效率低于分割-合并算法。
在油藏勘探过程中，有时须补测部分数据，这会引起地震数据的动态变化。因此需要对已存在的Delaunay三角网进行局部更新，按需插入和删除某些点。这就要求有一种能够快速构建三角网的方法。
考虑到工程应用的需要，本文在常规逐点插入算法的基础上，引入散乱点集有序化、三角形单元分类的方法快速构建Delaunay三角网，有效地解决了地震层位建模的工程问题。实测结果表明，改进后算法的时间复杂度大为降低，从而提高了算法的效率。
1 逐点插入算法流程
假设给定散乱数据点集P={Vi|i=1，……，n}，其中Vi的三维坐标为(xi，yi，zi)，n为点的个数。经过Delaunay三角化之后，输出Delaunay三角网D={Tj|j=1，……，m}，其中Tj表示第j个三角形，m为三角形的个数。逐点插入法的基本思想为：在某一已存在的三角网中插入一点，遍历所有三角形，查找外接圆包含该点的三角形集合，然后用LOP(Local Optimization Procedure)法则优化[4]，确保所建三角网符合Delaunay法则。以下用伪码描述逐点插入算法的流程。
SUBROUTINE：Delaunay Triangulate
输入：散乱数据点集P={Vi|i=1，……，n}
输出：Delaunay三角网D={Tj|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 自定义数据结构
