三角网数字地面模型及程序实现

1、三角网数字地面模型及程序实现袁功青（华南理工大学 土木与交通学院 ，广东 广州 510640）第一章 三角网数字地面模型基本原理及特点1.1 数字地面模型的含义数字地形模型（Digital Terrain Model，简称DTM或“数模”）是一个表示地形特征的、空间分布的、有规则的数字阵列，也就是地表单元平面位置及其地形属性的数字化信息的有序集合。它是地表2维地理空间位置和其相关的地表属性信息的数字化表现，可表示为：i=(i,i)=1, （式1-1）式中i是任一平面位置(i,i)的地表特有信息值，一般有：1) 基本地貌信息如高程、坡度、坡向等地貌因子；2) 自然地理环境信息如土壤、植被、气候、

2、地质分布等；3) 自然构筑物如河流、水系等和人工构筑物如公路、铁路、居民地等；4) 社会人文经济信息如人口分布、工农业产值等。根据不同的i值，其名称也稍有不同，如i为高程时，称为数字高程模型(Digital Height(Elevation) Model, DHM(DEM)；当i为土壤分布时，称为数字土壤模型。然而，不管怎样，都是对叠加在2维地理位置上的地表属性信息的数字描述，统称为数字地面模型，考虑到它的多样性，一般用DTMs表示。表达了平面位置(i,i)和i的空间相关关系，从这一意义上讲，DTMs是2.5维的而非3维。本质上讲，DTMs是对地形表面在地形采样数据基础上的表面重构。因此，对于

3、一个通用的DTMs系统，一般要经过数据采集、DTM建立和应用模型建立3个步骤，图1.1详细表达了DTMs系统的任务与建立步骤。地形表面信息解释遥感卫星数据航空摄影相片地形图数字化可视化基本应用地面测量数据交互修改 图1.1 DTMs系统的建立步骤与任务地表的采样数据是一个无序的数据集合，为便于计算机的管理及应用，需对该数据建立结构化的表达式，即隐式或显式的表达数据之间的拓扑关系，可用散点的形式和网格结构的形式来表达，目前主要是以网格结构的形式来表达。根据网格结构的不同，有规则网格(矩形、正三角形、正六边形网格等)和非规则网格(三角形、四边形网格等)。在规则网格中，矩形网格平面位置隐含于行列号中

4、，称为网格或rid，因网格结构简单而被普遍采用。三角形是平面图形中最简单的几何图形，因而在不规则的网格中经常使用，一般称为不规则三角网DTM或TIN。在应用中，数字地形模型一般又分为散点数模、规则格网数模（Regular Grid简称 RG）、不规则三角网数模（Triangulation Irregular Network，简称TIN），它们具有各自的特点、内插方法不适用范围。1.2 数字地面模型的种类及特点由于数模原始数据点的分布形式不同，数据采集的方式不同，以及数据处理、内插的方法不同和最后的输出格式不同等原因，数字地面模型的种类较多。根据数模中已知数据点的分布形式并考虑到数据输出格式及数

5、据处理方式，可将数字地面模型大致地分为规则数模、半规则数模和不规则数模三大类。各类数模的主要特点如下：1.2.1 规则数模 规则数模是指原始地形点之间均有固定的联系，如方格网数模、矩形格网数模和正三角形格网数模等。在格网之间待定点的高程，常采用局部多项式进行内插。由于每个已知点相对于周围已知点的位置是固定的，所以若按规则格网方式采集数据建立数模，量测地形点简单、客观，不需判读地形，易于实现数据采集的自动化、半自动化。由于格网是规则等距的，在计算机中只需储存各个格网节点的高程值，而平面坐标只需记录第一个格网节点即可，其余节点的平面坐标根据其格网管理信息很容易在计算机中确定和恢复，这可大为节省计算

6、机内存。这种只记录节点高程的数模，也称为数字高程模型(Digital Elevation Model 或 Digital Height Model，简称DEM或DHM)。该类数模的另一优点是输出形式简单、数据结构良好、便于应用，内插待定点高程时，检索与内插简单快速。这种数模的最大缺点是原始数据不能适应地形的变化，除十分均匀的地形外，已知点没有与地形特征点联系起来，易遗漏地形变化点。由于数据采集按规则格网方式进行，一旦间距给定，所有已知点平面位置就是固定的，从而导致地形变化大的地方地面信息不足，而地形均匀、平缓的区域冗余数据点太多这种精度不一的现象，此外，由于规则格网节点不能兼顾地形变化线和地形

7、特征点，格网中也就难于确定地面坡度的变化，从而导致高程内插精度降低。若要使规则格网数模更好地表示地形，则只有将格网间距缩小，这将导致原始数据的采集工作成倍增加。规则格网数模一般适用于地形较平缓和变化均匀的区域，以及用于搜索地形等高线、绘制地形全景透视和对内插速度要求极高的路线平面优化中内插地面线等方面。1.2.2 半规则数模半规则数模是指各原始数据点之间均有一定联系，如用地形断面或等高线串表示的数模。当以断面线高程表示地形时，任一串原始地形点均表示某一个特定的地形剖面。各断面之间的距离及断面上点与点之间的距离可以是固定的，也可随地形变化而定。这种数模相当于用一批密集的、相互关联的地形剖面叠合在

8、一起模拟地形，待定点的高程由相邻断面的已知点高程内插得到。如一种称为鱼骨式数模的数字地面模型就属于此类，这是在路线设计方案确定以后,沿路线纵、横断面采集地面点数据，用一批相互关联的横断面地形剖面叠合在一起而形成的数字地面模型。 沿等高线采集的一系列同一高程的X、Y坐标的地形高程信息(二维串线)，以及沿地形特征线、断裂线、地物、水系等各种信息采集的一系列X、Y、Z三维坐标信息(三维串线)所组成的数模，称串状数模。由于等高线具有在地形坡度大的地方密，在地形平缓处稀的特征，故地形点的分布与密度能很好地适应地形的变化。此外，由于各种断裂线、地物和水系信息可用三维串线在数模中表示，方便地进行处理，使程序

9、功能大为加强，这是串状数模最为突出的优点。英国的MOSS系统所采用的就是串状数字地面模型。 半规则数模能较好地适应地形变化，内插精度较高，但数据采集不能实现自动化，原始数据的分布与密度易受操作人的主观影响，建立数模过程中的程序处理较规则数模复杂。1.2.3 不规则数模 不规则数模其原始地形数据点之间无任何联系，点的分布是随机的，一般常采集地形特征点、变坡点、反坡点、山脊线、山谷线等处，常见的有散点数模、三角网数模等。 散点数模是将原始地形点看作一些随机分布的“离散点”，可认为点与点之间无任何联系。在这种已知点中直接进行内插其精度是不可靠的，因为不能确定由哪些点构成实际的地表面。所以内插时先得利

10、用已知点拟合一个局部的或区域的内插表面，然后再由该内插表面确定待定点的高程，这是这一类散点数模的共同特点。这种高程内插也称为移动曲面拟合法，在坡度变化均匀的地形条件下，这种拟合较符合地形的实际情况，可达到较高的内插精度。在坡度变化大，地形起伏急剧改变的区域，特别是地物、地形断裂线附近，这种曲面拟合则明显地难于很好地描述出地形局部形状，而造成内插点的高程误差偏大。所以，为提高散点数模的内插精度，使之满足工程设计需要，除原始数据点的分布必须符合地形变化，保持足够的密度外，程序中必须具备对地物、地形断裂线处理的功能。 从数模的精度和计算速度两方面来考察，散点数模不失为一种简单而有较的方法，具有很大的

11、实用价值，在国内外许多实用程序中广泛地得到应用。 不规则数模总的特点是：数据采集是随机的，一般都是取地形特征点，所以能较好地适应地形变化，内插精度较高。其缺点是采样需要人工判读地形，从而增加了数据采集的难度，此外构造数模较复杂，计算时间较长。由于该类数模优点较为明显，所以应用最为广泛。1.3 三角网数字地面模型基本原理三角网数字地面模型(Triangulated Irregular Network)TIN是用一系列的互不交又、互不重复的三角形逼近地形表面，与Grid相比，TIN在拟合地表上更灵活、更精确。但由于结构的不同，Grid许多成熟的技术并不能完全移植到TIN中。从结构上讲，TIN是一典

12、型的矢量数据结构。它主要通过节点(地表采样点)、三角形边和三角形面之间的关系来显式或隐式的表达底细功能散点的拓扑关系，因此设计一个高效的、结构紧凑的、维护方便的TIN的存储与组织结构对TIN的应用与维护是至关重要的。Tin的基本单元三角形的几何形状直接决定着TIN的应用质量。由于地形的自相关性，相互越接近的地形采样点，其之间的关联程度越大;同时，理论与实践均证明，狭长的三角形其插值精度比规则的三角形插值精度可信度要低。因此，在Tin中对三角形的几何形状有着严格的要求。一般有三条原则:l)尽量接近正三角形;2)保证最近的点形成三角形;3)三角形网络唯一。一个良好的数据结构和三角划分准则，必须由一

13、个高效的算法和程序实现。算法在具体应用中发挥的作用由算法本身的性能和实现它的程序质量共同决定，而程序的好坏在很大程度上依赖于算法的原理。对算法本身在理论上进行分析论证，寻求一种高效率、高精度、适用面广的算法是TIN建立中的重要一环。由上述的分析知，TIN的数据组织+三角划分准则+算法和程序构成了TIN的基本理论体系框架。图2表示了这一结构。三角划分准则三角划分准则三角划分准则三角划分准则数据组织结构图2 TIN模型的理论体系结构第二章 三角网生成算法目前国内外建立不规则三角形网的主要方法有三类：自动联结三角网法；径向扫描法（Radial Sweeping Algorithm简称 RSA）； 泰

14、森三角形法（Thiessen）。下面分别对这几种方法作一些介绍和评价。1 自动联结三角网法自动联结三角形网，要求建立在可能获得最佳三角形的条件下。这个条件就是在使用周围邻近的离散点组成三角形网时，应尽可能地确保每个三角形都是锐角三角形或三边的长度尽量相等，避免出现过大的钝角和过小的的锐角。其判别法则是：最小距离和法则与最短外接圆半径法则。最小距离和法则是指待插入点到基边两端点的距离和为最小，在用程序具体实现时，是用“该点对基边两端点张角最大或其余弦值最小”原则来实现的。最短外接圆半径法则是指三角形外接圆半径最小。从余弦定理公式：c2=r2+ r2-2 r2cos() （式2-1）可以看出：co

15、s()-1/ r2。其中：c代表基边的固定长度，r代表外接圆半径。由于圆心角2，即有 cos(2) -1/ r2，所以两种判别法则是保持一致的。 自动联接三角形网法的具体思路是：1）确定第一个三角形在散点域中任取一离散点A作为第一个三角形的第一顶点，找出距离该点最近的点B作为第一个三角形的第二顶点且构成基边AB，应用余弦定理找出和该基边对角为最大的点作为第一个三角形的第三顶点。值得注意的是：为了顾及地形特征，找第三点时还必须遵守 两条规则：不跨越地形特征线条件；不进入禁区（封闭区域）条件。2）三角形的扩展自动联结三角形网法中的三角形扩展是从第一个三角形的三条边依次开始的，应用余弦定理来找到第三

16、点。为了防止三角形被遗漏、重复或发生交叉，每形成一个新的三角形，都要与己形成的三角形作一些比较判别，如果重复或交叉，则该三角形不记录。自动联结三角形网法的优点是：原理明了，不考虑地形特征时编程简单。但其最大的缺点是：仅仅从最佳图形结构考虑，从而忽视了构网时的连续反复多次的三角形比较，使构网速度相当慢，导致数模的效率低下。而且一旦考虑地性线的介入，程序就变得复杂了。2 径向扫描法（Radial Sweeping Algorithm简称 RSA）径向扫描法（RSA）是Miratl和 Weingartan首次发表于一九八二年，其基本思想是分成以下两个阶段来完成不规则三角形网的建立。1）建立稀疏的辐射

17、三角网首先是选择离区域形心最近的点作为辐射中心，由这一点向其它地形点作射线并计算其相对于辐射中心的距离和方位，对这些离散点以方位作为关键词进行排序，若方位相同则按距离，再把这些排序好的点与它的后继点连一线段，形成一系列细而长的三角形。如果该点同前一点居于同一方位，则前一点、该点、后继点三点组成新的三角形。在处理过程中，用一个链表或动态数组将每一个点的信息存储起来，形成一个外边界。然后依次从链表或动态数组中连续取出三点，判别是否能形成内部三角形，若能，则记录该三角形，且把第二点从边界表中去掉。每当优点删去时，必须从头开始去进行点检验，最后得到一个凸形的闭合边界。2）优化初始三角网径向扫描法（RS

18、A）法在优化三角网过程中考虑了两个因素，一是图形结构；二是地形起伏状况。显然，就图形结构而言，等边三角形是最佳图形；就地形起伏情况而言，为使三角形网能很好地覆盖于地表，三角形边不应有“悬空”和“切割” 现象。具体到该法中，在考虑图形结构时，若相邻两个三角形组成凸四边形，应检查四边形两对角线之比，如果超过约定值，则以短边作为三角边；在考虑地形起伏清况时，则检查两对角线高差之比，一般选取高差大的两对角点连线作为三角形边。当二者出现矛盾时，视其影响而定。经过优化处理后，所得的三角网就比较合理了。径向扫描法（RSA）法思路简明，容易实现。但是，约定值若没有一个客观的衡量标准，导致网形不唯一，特别是引入

19、地性结构线后，难免会出现三角网的“悬空”或“切割”边，致使歪曲地形；同时三角网的优化工作也是一个相当复杂的过程，故该法在建立高精度的DTM时是不可取的。3 泰森（Thiessen）三角形法泰森多边形的概念是:将分布在平面区域上的一组离散点用直线分隔,使每个离散点都包含在一个多边形之内。进行分隔的规则是,每个多边形内只包含一个离散点,而且包含离散点Pi的多边形中,任意一点Q到Pi的距离都小于Q点到任一其他离散点Pj的距离。把每两个相邻的泰森多边形中的离散点用直线连结,生成三角形网,这样数模便完成。泰森模型的显著特点是:每个三角形的外接圆内不包含其他离散点,而且三角形的最小内角达到最大值,又称De

20、launay三角形和Delaunay三角形格网。泰森三角形法是基于图论学的表面三角剖分方法。该法数学理论体系严密且从图论角度出发，因而具有良好的图形结构和网形唯一性。但是有两个限制条件：散点域必须为凸域；散点域中无四点共圆。由于其优点明显，国内外对其研究也多，发明了许多算法。本文将在下一章就改进的Thiessen三角形法，即：Delaunay三角形法构网建立数模进行详细介绍。第三章 三角网数字地面模型程序实现在地学领域，经常需要处理大量分布于地域内的离散数据。为了合理利用一定地域内不均匀分布的离散地形点，1908年，G.Voronoi首先在数学上限定了每个离散点数据的有效作用范围，即其有效反映

21、区域信息的范围，并定义了二维平面上的Voronoi图 (简称为V-图，也称为Thiessen图，Dirichlet图，Vigner-Seithz图)。1911年，荷兰的应用V-图进行了大区域内的平均降水量研究。1934年，B.Delaunay由V-图演化出了更易于分析应用的Delaunay三角网。从此，V-图和Delaunay三角网就成了被普遍接受和广泛采用的分析研究区域离散数据的有力工具。当然，它的应用不仅适用于地学，而且活跃于所有与2.5维分析有关的领域。3.1 Delaunay三角网的定义1）Delaunay三角化的理论基础一Voronoi图1908年，首先在数学上限定了每个离散点数据的

22、有效作用范围，即其有效反映区域信息的范围，并定义了二维平面上的Voronoi图(以下称为v-图)。1934年Boly由v图演化出了更易于分析应用的Delaunay三角网(以下称为D一三角网)。从此，V一图和D一三角网就成了被普遍接受和广泛采用的分析离散数据的有力工具。在TIN的生成中，V一图和D一三角网有着广泛的应用。Voronoi图是由不规则网格定义的最基础和有用的结构之一。作为一个有着广泛用途的几何结构，V一图可以应用于许多领域，并且经常以首先将它应用与某一特定领域的人的名字来命名它。V一图在地理方面的应用是从气象学家开始的。他用V一图去解释那些有着不均匀天气状况分布的地区，因此在此领域V

23、一图又称为Thiessen多边形。在图3.1中可以看到，V一图是通过切分一个中心点和它周围点之间的连线来定义的。切分线和连线之间的关系是互相垂直的。当我们对整个区域中的每一个点应用这种方法时，整个区域会被相邻的多边形所覆盖。图3.2是由一组自由分布的点组成的V一图。图3.1 Voronoi多边形的构建 图3.2 由一组自由分布的点组成的V一图Voronoi多边形有如下五点性质：l 对于平面区域V上给定的K个离散点，将S分成K个Voronoi多边形的分法是唯一的。相应地，作为伴生图形的Delaunay三角网唯一。l 任意两个Voronoi多边形不存在公共区域。相应地，作为伴生图形的Delauna

24、y三角网中没有三角形重复、交叉。l Voronoi多边形是一凸多边形。l 点集V中任意三点vi、vj、vk构成Delaunay三角形的充要条件是：过vi、vj、vk三点的圆内不含有V中的其余任何点。l 由点集V中 N个散点构成的 Delaunav三角形数是：Nt=2(N-1)-Nb，三角形的边数是：Ne=3(N-1)-Nb，其中Nb是V-图外围边界顶点个数。2） Delaunay三角网的定义作为Voronoi图的对偶，Delaunay三角网与Voronoi图紧密相关，它以第一位发现这一对称关系的科学家B.Delaunay的名字命名。Delaunay三角网的定义是:连接所有相邻的Voronoi多

25、边形的生长中心所形成的三角网称为Delannay三角网。图3.3 一个点集的Voronoi图(虚线)和它的Delaunay三角剖分(实线)。3.2 Delaunay三角网的特性Delaunay三角网有几个非常重要的性质，如下所示：(l) 空外接圆性质:设T是平面有限点集V的一个三角剖分。T中的一个三角形的外接圆如果不包含V中的其它点，我们就称它满足空外接圆性质。V的一个三角剖分:，当且仅当它的每一个三角形都满足空外接圆性质时，它才被称为一个Delaunay三角剖分(见图3.4)。图3.4 (a)不是Delaunay三角剖分 (b)是Delaunay三角剖分(2) 最大一最小角性质:设T是V的一

26、个三角剖分，e是T的一条边，Q是T中和e相邻的两个三角形组成的四边形。当交换Q中的对角线不增加6个内角中的最小角时，我们称边e满足最大一最小角性质。满足最大一最小角性质的边也被称为局部优化的边(见图3.5)。V的三角剖分T，当且仅当T的每一条边都局部优化时，才是Delaunay三角剖分。图3.5加深的边在 (a)中没有局部优化，在(b)中被局部优化(3) Delaunay三角网网形唯一性:Delaunay三角网是建立在控制点周围最紧密相连的双重格网，这种网形只与离散点的分布有关，而与建网的途径即离散点参加构网的先后顺序无关。所以，只要离散点平面位置相对固定，其网形都一样。由于几个性质，决定了D

27、elaunay三角网具有极大的应用价值。同时，它也是二维平面三角网中唯一的、最好的。3.3 Delaunay三角网构网方法在GIS中，2.5维的分析处理由DTM(数字地形模型)模型进行。由上节可知DTM主要由格网（rid）与不规则三角网 (TIN)两种数据格式来表示，而以后者更为重要。格网的优点是：充分表现了高程的细节变化，拓扑关系简单，算法实现容易，某些空间操纵及存储方便；而它的局限性是：占用巨大的存储空间，不规则的地面特征与规则的数据表示之间的不协调。不规则三角网的优点是：高效率的存储，数据结构简单，与不规则的地面特征和谐一致，可以表示线性特征和迭加任意形状的区域边界，易于更新，可适应各种

28、分布密度的数据等；而它的局限性是：算法实现比较复杂和困难。在现今的GIS系统中，基本上均支持以上两种数据格式，以不规则三角网为主，格网为辅，并提供相互转换功能。历经二十余年的不懈努力，TIN的许多算法难关已被攻克。TIN的生成算法中，三种算法：分割-归并法、逐点插入法和逐步生长法最终被普遍接受和采用，而又以前二者更为流行。分割-归并法时间效率高，但占用内存空间较多；逐点插入法时间效率较差，占用内存空间较少。现在有人在这两种算法的基础上，提出了一种综合性能的算法合成算法，它的性能优于以上两种算法。下面就Delaunay三角网的三类建网方法：分治算法、逐点插入法、逐步生长法作进一步介绍和评价。3.

29、3.1 分治算法（Divide-and-Conquer）Shamos和Hoey提出了分治算法思想，并给出了一个生成V-图的分治算法。Lewis和Robinson于1978年提出了一种建立三角网的方法，以处理限定边界的等高线和有限元分析应用，他们将分治算法思想应用于生成Delaunay三角网，他们给出了一个“问题简化”算法，递归地分割点集，直至子集中只包含三个点而形成三角形，然后自下而上地逐级合并生成最终的三角网。以后Lee和Schachter又改进和完善了Lewis和Robinson的算法，提出了对平面点建立Delaunay三角网的分而治之算法。这种算法首先将数据排序，分成两个互不相交的子集，

30、在每一子集建立三角网后，将两个三角网合并以生成最终的狄洛尼三角网。Dwyer（1987）对分而治之算法作了一些改进。通过将数据分成垂直条块，进而用相交边界将条块再一次细分为区域，并应用带约束条件的Lawson LOP将对角线进行交换，Dwyer的算法能处理带约束条件的数据。Lee和Schachter算法的基本步骤是：1 把点集V以横坐标为主，纵坐标为辅按升序排序，然后递归地执行步骤26；2 把点集V分为近似相等的两个子集VL和VR；3 在VL和VR中生成三角网；4 用Lawson提出的局部优化算法（LOP）优化所生成的三角网,使之成为Delaunay三角网；5 找出连接VL和VR中两个凸壳的底

31、线和顶线；6 由底线至顶线合并VL和VR中两个三角网。以上步骤显示，分治算法的基本思路是使问题简化，把点集划分到足够小，使其易于生成三角网，然后把子集中的三角网合并生成最终的三角网，用局部优化算法（LOP）保证其成为Delaunay三角网。不同的实现方法可有不同的点集划分法、子三角网生成法及合并法。3.3.2 逐点插入法（Iteration）Lawson提出了用逐点插入法建立Delaunay三角网的算法思想，它是将未处理的点加入到已经存在的Delaunay三角网中，每次插入一个点，然后将Delaunay三角网重新定义。Lee和Schachter，Bowyer，Watson，Sloan，Mace

32、donio和Pareschi，Floriani和Puppo，Tsai先后进行了发展和完善。逐点插入算法的基本步骤是：1 定义一个包含所有数据点的初始多边形；2 在初始多边形中建立初始三角网；3 插入一个数据点P，在已经存在的三角网中找出包含P的三角形t，把P与t的三个顶点相连，生成三个新的三角形；4 用Lawson的LOP算法优化局部三角网；5 是否所有离散点都参与了构网？若是，构网完毕；否则，转到3。从上述步骤可以看出，逐点插入算法的思路非常简单，清晰。先在包含所有数据点的一个多边形中建立初始三角网，然后将余下的点逐一插入，用LOP算法确保其成为Delaunay三角网。同时，该算法占用内存小

33、，时间复杂度为O(N5/4)，运算速度较快，是一种比较理想的Delaunay三角网构网方法。各种实现方法的差别在于其初始多边形的不同以及建立初始三角网的方法不同。3.3.3 逐步生长法Green和Sibson首次实现了一个生成Dirichlet多边形图的生长算法。他们在1978年提出的递归方法通过顺序扫描数据点，计算平面Dirichlet多边形。在此过程中随着新点的加入，以递归的方式改变新点周围Dirichlet多边形的邻接关系。Brassel和Reif稍后也发表了类似的算法。这种算法从任意一点开始，在一维排序链表中以循环搜索方式查找此任意点的泰森邻域点，将建模区域划分为泰森多边形，然后把这一

34、点与其所有泰森邻域点连接以形成了Delaunay三角网，并由此生长，直至所有点都包含到三角网中。McCullagh和Ross通过把点集分块和排序改进了点搜索方法，减少了搜索时间。Maus也给出了一个非常相似的算法。逐步生长法的基本步骤是：1以任一点为起始点，找出与起始点最近的数据点相互连接形成Delaunay三角形的一条边作为基线；2按Delaunay三角网的判别法则(即它的两个基本性质)，找出与基线构成Delaunay三角形的第三点；3基线的两个端点与第三点相连，成为新的基线；4迭代以上两步直至所有基线都被处理。上述过程表明, 逐步生长法的思路是，先找出点集中相距最短的两点连接成为一条Del

35、aunay边，然后按Delaunay三角网的判别法则找出包含此边的Delaunay三角形的另一端点，依次处理所有新生成的边，直至最终完成。该算法时间复杂度为0（N3/2），时间效率较低，现已基本淘汰。各种不同的实现方法多在搜寻“第三点”上做文章。Tsai为比较算法性能，给出了一张各种算法的时间复杂度对照表，如表1所示。表1几种Delaunay三角网生成算法的时间复杂度算 法一般情况最坏情况Lawson(1977)2O(N4/3)O(N2)Green和Sibson(1978)3O(N3/2)O(N2)Lewis和and Robinson(1978)1O(N logN)O(N2)Brassel和R

36、e if(1979)3O(N3/2)O(N2)MaCullagh和Ross(1980)3O(N3/2)O(N2)Lee和Schachlter(1980)1O(N logN)O(N logN)Lee和Schachlter(1980)2O(N3/2)O(N2)Bowyer(1981)2O(N3/2)O(N2)Watson(1981)2O(N3/2)O(N2)Mirante和Weigarten(1982)3O(N3/2)O(N2)Sloan(1987)2O(N5/4)O(N2)Dwyer(1987)1O(N logN)O(N logN)Chew(1989)1O(N logN)O(N logN)Mac

37、edonio和Pareschi(1991)2O(N3/2)O(N2)上标说明：1. 分割-归并法；2. 逐点插入法；3. 逐步生长法。表中N为数据点数，O(f(N)表示算法的时间复杂度，它以算法中频度最大的语句频度f(N)来度量。3.3.4 一种新的合成算法由上面的介绍我们可以看到，目前采用较多的前两类算法各具优势又有局限，同时，它们又具有明显的互补性。分治算法时间性能好，空间性能差；逐点插入法空间性能好，时间性能差。我们评价一个算法的优劣是看它对时间和空间的消耗，即时空性能的综合表现。因此，就产生了一种合成算法即把以上两种算法结合起来，取长补短，来提高算法的性能。这种结合的具体做法是：以分治

38、算法为主体，当递归分割数据集的过程进行到子集中的数据量小于一个预定值分割阈值时终止，然后用逐点插入法在子集中生成子三角网。这种算法就是合成算法。它的流程图见图5。开 始插入法插 入 法（TL/TR）排 序寻找上下切线模块HULLNL>Nd?NR>Nd?合成算法TL合成算法TR分割点集VNv>Nd合并模块MERGE结束NYNNYY图3.6合成算法主程序流程图其中表示数据集：是的数据量；是分割阈值；,分别表示两个子集的数据量；,分别表示在子集中建立的两个子三角网。图3.6表明，为便于分割，首先按升序以横坐标为主，纵坐标为辅把数据排序。然后比较与，如>则分割数据集，否则用逐点

39、插入法生成三角网。合成算法包含3个模块：逐点插入法模块、寻找上下切线模块以及合并模块。1逐点插入法模块逐点插入法有多种实现方法，该算法选用acedonio和areschi提出的方法。在这个方法中，把数据中的凸壳作为超级多边形。凸壳就是由数据中的凸集顺序相连形成的多边形。逐点插入法模块由以下4个步骤完成。(1)找出凸壳。这一步用arkin 在 1991年提出的算法来做。为提高搜索效率，先将数据集划分为矩形栅格。设每个栅格单元内包含的数据为4，则栅格的行列数为/4。因具有,+,-的最大值及最小值的点是凸集中的点，把它们按逆时针方向顺序相连作为初始凸壳。然后用递归函数convex (,)完成整个凸壳

40、。,是凸壳上的相邻点。convex (,)找出凸壳上,间的另一点。如果不存在,则返回空。是与相交的栅格中位于右侧且与距离最大的点，或当右侧无点，但上有点，且位于,之间的点。(2)建立初始三角网。以凸壳上值最小的点为出发点，按序与其余点相连。(3)插入点。先找出包含要插入的点的三角形，然后连接与的3个顶点。(4)优化。用Delaunay三角网的空心圆性质考查新生成的三角形，如不满足，则调换与相邻三角形所组成的四边形中相连的对角线。向相邻三角形扩展此过程直到满足条件为止。后两步循环执行直到插入全部数据。2寻找上下切线模块这个模块找出连接左右两个子三角网的凸壳和的上切线和下切线。设是中值最大的点；是

41、中值最小的点；是上逆时针方向的邻点；是上顺时针方向的邻点。寻找和的下切线过程由两个循环完成。(1) 如位于右侧，则成为新的，逆时针方向的邻点成为新的，否则退出。(2) 如位于右侧，则成为新的，逆时针方向的邻点成为新的，否则退出。用类似的方法可以生成和的上切线。3合并模块子三角网和的合并由一个循环完成。以和的下切线为起始基线，分别在和中找出与它构成Delaunay三角形的第3个顶点1和1，在这两个三角形中选择外接圆半径小的连接到三角网中。新生成的边作为新的基线继续这个过程，当基线成为和的上切线时终止。在寻找上下切线模块和合并模块中，要频繁用到两个函数pred (,)和succ (,)。pred

42、(,)是在包含公共端点的一簇线段中，得到与边顺时针方向邻边的另一个端点。succ (,)相反，是得到与边逆时针方向邻边的另一个端点。3.4 局部优化过程（LOP-Local Optimization Procedure）Lawson提出的局部优化过程是所有生成Delaunay三角网算法都要用到的关键过程，其理论依据就是Delaunay三角网的空外接圆性质。前面所讲到的采用逐步插入法插入一个点P时，局部三角网优化过程是：三角形入栈出栈Delaunay三角化入栈出栈，直到栈空为止，其步骤如下（示意图如图3.7所示）： 图 3.7 局部优化过程（1）、建立堆栈A和B；（2）、连接点P与其所在三角形（

43、t4）的三顶点，构成三个三角形（t4、t5和t6）；(3)、初始化堆栈。把以P为顶点的三角形（t4、t5和t6）放入堆栈A中，把与P相对的三角形（t3、t1和t2）放入堆栈B中，且必须保证两组三角形中处于同一位置的两个三角形相邻；(4)、出栈。检查A、B（其实A、B是同步的）是否为空，若是，则过程完毕；否则分别从栈A、B中取出一个三角形（t4、t3）；(5)、Delaunay三角化。检查P是否在从栈B中取出的三角形（t3）的外接圆内（不包括刚好在圆上）。若不是，转4；否则，由这两个三角形形成凸四边形，交换其对角线形成新的三角形且赋与相同的编号（t4、t3）；(6)、入栈。把刚形成的三角形（t4

44、、t3）放入栈A，把新的与P相对的两个三角形（t7、t8）放入栈B，且注意次序。再转4。3.5 初始多边形及三角形的形成Delaunay三角网的外边界是凸多边形，采用逐点插入法首先得建立初始多边形，即“凸壳”及初始三角网。生成凸壳的方法、凸壳的形状多种多样，相应地生成的初始三角形也不一样。3.5.1 构造三角形作为凸壳这里介绍一种采用三角形作为初始多边形的方法，其方法大致如下：1计算dmax的值。dx=xmax-xmin，dy=ymax-ymin，dmax=max（dx，dy）2标准化坐标值。目的是使点的x、y值限定在01的范围内。x=（x-xmin）/dmax，y=（y-ymin）/dmax

45、3构造凸壳。有关文献中非常形象地采用了术语“supertriangle”，意思是“特大三角形”。方法是在原有点集中增加三个点，用标准化坐标表示分别为（-100，-100）、（100，-100）和（0，100）。这三点构成的三角形包容了整个点集，既是初始多边形，又是初始三角网。其实为构造凸壳增加的三点可任意取值，条件是包容整个点集。该法的缺点是标准化坐标会影响数据精度。3.5.2 找边界点构造凸壳该方法是首先定义初始凸壳，用程序递规找出其余边界点，构成凸壳。然后，用凸壳上某一点为固定点，与其他点连直线构成初始三角网。最后借助Lawson的LOP算法优化初始三角网成为新的初始三角网。该算法麻烦，比

46、较费时。3.5.3 构造四边形作为凸壳如果采用四边形作为初始凸多边形，其过程如下：1计算散点集的向最小值xmin、最大值xmax和y向最小值ymin、最大值ymax，且进行修正。xmin=xmin-TOL；xmax=xmax+TOL；ymin=ymin-TOL；ymax=ymax+TOL；2构造凸多边形构造四个虚点并增加到散点集中：A（xmin，ymin，0），B（xmax，ymin，0），C（xmax，ymax，0）和D（xmin，ymax，0）。3构造初始三角网首先按逆时针顺序构造两个三角形：一个由A、B、C三点构成，三角形编号为0；另一个由A、C、D构成，三角形编号为1。然后利用LOP算

47、法优化初始三角网。3.6 点所在三角形的判别方法采用逐点插入法时，程序首先需要找到包含该点的三角形。在空间分析学中，将点元素与面元素的空间拓扑计算归结为著名的“点在多边形中（Point-in-polygon）”的识别问题。三角形是一个凸多边形，判别方法更多，目前主要有以下几种：面积法、交点法、方向角法和循环法。3.6.1 面积法基本原理：根据三角形之间的面积关系来判定。如图3.8（a）所示，把待判定点P与三角形t的三个顶点A、B、C直线相连，形成四个三角形，由数学上的“图解法”可得出如下结论：1若SABC=SABP+SBCP+SCAP，则点P在三角形t内；2若SABCSABP+SBCP+SCA

48、P，则点P在三角形t外；评述：该法搜索三角形的量大，计算量也大，因而速度也慢。（a）面积法 （b）交点法 （c）循环法图3.8 包含点P的三角形的判定方法3.6.2 交点法 有的资料上称作“垂线法”或“半线理论”，其基本原理：计算并根据过点的某一射线与三角形相交的交点的分布情况来判定。如图3.8（b）所示，由点P向右作平行于x轴或向下作平行于y轴的射线，（不用计算）判定该射线与三角形t的边的交点个数，并根据交点个数得出如下结论：1)、若交点数为奇数，则点P在三角形t内；2)、若交点数为偶数，则点P在三角形t外；需注意的是：在统计交点数时，当交点为三角形t的一顶点时，采取“舍下取上”的原则，即：

49、当某线两端点的y坐标均小于等于yP时，不计交点；当某线两端点的y坐标只要一个大于yP时，就计交点。该方法在明确判定指定的三角形时还比较简单，且对凹凸多边形均适用，特别在基于闭合特征线的网形调整时非常实用。但在未明确对象时，搜索的三角形数会很多，计算量很大，导致速度慢。3.6.3 方向角法有的资料上也称为“夹角和法则”。基本原理是待定点P与三角形t的三个顶点连线段，按顶点次序各线段的夹角值代数和为一常数，如图3.9所示，其中各夹角规定逆时针为正，顺时针为负。（a）P在t内 （b）P在t外图3.9 包含点P的三角形的判定方法（二）方向角法由图可以得出以下结论：若方向角之和等于360º，则

50、点P在三角形t内；若方向角之和等于0º，则点P在三角形t外；该方向交的计算复杂，运行速度慢，且不便于识别三角形边界的点。3.6.4 循环法 循环法的基本原理是根据三角形的数据组织结构中顶点是按顺时针（clockwise）还是逆时针（anti-clockwise）排列及点P在三角形边的左右来判定。如图3.10所示，结论如下：当三角形顶点为逆时针排列时，若点P都在三角形三边的左侧，则点P在三角形t内；换个角度来看，若点P在某一边的右侧，则点P有可能在该边右侧的三角形内，使得搜寻三角形有某种方向性，如图3.10所示；图3.10 循环法定位点P所处三角形循环法只适用于凸多边形，对TIN是适宜

51、的。该法搜寻方向明确，提高了查找速度；思路明确，公式简单，计算量不大。3.7 数据组织结构的原则数据组织结构是DTM的基石，它的功能强弱将左右着整个数模系统的功能和性能。数据组织结构原则如下：数据结构的设计必须服务于系统功能的设计。它要为系统提供完整的数据信息，以利于各种应用；在满足系统功能的前提下，数据组织结构必须尽量简单，仅仅占用最小的存贮空间；快速查询功能。数据查询是DTM的各项功能中使用最为频繁的，但是DTM的数据量极大且杂，遍历数据是不可取的。快速查询意味着要付出一定的空间代价，须建立各种辅助索引，以换取查询速度的提高；动态拓扑结构。早期拓扑数据只是用于检查图形正确性的冗余数据，现在

52、拓扑查询之需已成为DTM的必备成分，原先拓扑数据结构是静态的，对任何图形的局部修改都需重建全局的拓扑关系，动态拓扑结构则能实时地响应局部拓扑关系的修改，而无需全局重新建立，提高了系统的速度。3.8 三角网数模的排序与检索目前，带状数字地面模型的数据检索有三级检索和二级检索机制两种。三级检索是把整个路段作为一个整体模型，建立“子区子块点、三角形链表”的检索机制，子区是沿路线划分建立的，如图3.11所示。图3.11 三级索引的子区划分有时候我们一方面整个路段的数据量很大，需很大计算机空间；另一方面，考虑到路线设计一般是分段设计的，没有必要这样做。也可以采用二级索引机制，以文件的形式分成几个子区，各

53、子区用均匀（正方形）格网建立“子块点、三角形链表”的索引。子区内数据组织如下所述。1计算子区地形散点的包容盒。Xmin=Xmin-TOL；Xmax=Xmax+TOL；Ymin=Ymin-TOL；Ymax=Ymax+TOL；其中：TOL为修正值，解决点落在格网线上的问题。2计算格网边长大小。经检验每个格网中约放n 个散点时，速度较好。则格网边长size为：size （式3-1）计算x、y轴向各自的格网数x_res,y_res。x_res int() （式3-2）y_res int() （式3-3）则总的格网数为：allres x_res×y_res （式3-4）3.格网排序。为了减小空

54、间及加快搜寻速度，可以对格网的组织由二维转换为一维连续排序。排序过程见图3.12所示。图3.12 格网排序过程4确定点P（x，y）所在的格网号。i int() （式3-5）j int() （式3-6）根据i、j计算格网号m_nBin:(i) 若i为奇数，则有：m_nBin=(i-1)*y_res+j (ii)若为偶数，则有：m_nBin=i*y_res-j+1 3.9 顾及地形特征线的网形调整地形特征线，包括山脊线、山谷线、陡坎等地形线及河流、房屋等地物轮廓线。在把这些特征线数据点与一般地形点一起完成初级网的建立之后，再用其信息进行网形调整，将使建立高精度的数字地面模型成为可能。3.9.1 基

55、本原理（1）、多边形的次三角剖分对一任意多边形的次三角剖分，虽然各次的三角形形状不同，但是其数目相等，不含边界的内部边数也相等。如图3.13所示。图3.13 多边形的N次三角剖分（2）、凹凸点判别原则对于按给定方向（顺时针或逆时针）组成的多边形01i-1ii+10，其中i相对于i-1，i+1的凹凸性为：（参见下图3.14）图3.14 凹凸点判别法则det （式3-7）当i-1、i、i+1为顺时针排列时：det0 i为凸点det0 i为凹点det0 i-1、i、i+1 三点共线当i-1、i、i+1为逆时针排列时：det0 i为凹点det0 i为凸点det0 i-1、i、i+1 三点共线3.9.2 对地形特征线的预处理（1）、加密地形特征线上的点一条地形特征线上的数据采集是在转折处（水平和竖向）发生的