1 引 言

地形特征线又叫地性线，包括山脊线和山谷线。它描述了地形的骨架结构，揭示了地貌形态的本质，在数字高程模型(DEM)生成、地貌综合、流域分析、水文分析、河系自动生成等研究中具有重要作用。目前，地形特征线的自动提取算法均源于两种原理：基于几何形态分析和基于地形表面流水分析^[1]，具体的方法前者有等高线曲率判别法、多因子特征提取法^{[2, 3]}、等高线骨架化法。根据数据源不同还可以将地形特征线提取方法分为基于等高线的地形特征提取^{[4, 5, 6, 7]}和基于DEM的地形特征提取^{[8, 9, 10, 11, 12]}两种，本文的方法属于前者。

与其他方法相比，利用Delaunay三角网(Delaunay triangulation，DT)模型和约束型Delaunay三角网(constrained Delaunay triangulation，CDT)^{[13, 14, 15, 16, 17]}来探测和表达曲线的弯曲特征能较好地顾及弯曲的对称性与层次性。然而，简单地将DT模型用于等高线弯曲检测，相邻曲线上的点会相互干扰，造成一些关键弯曲无法提取，或者提取出一些小于人们认知意义上的弯曲^[15]。为避免这种干扰文献^[15]提出了联合DT模型。但是，由于等高线的嵌套特性，基于单条等高线生成的三角网之间会存在大量的相交三角形，而且由于构成等高线的离散点分布疏密程度不同，在相邻等高线的影响下，难免会遗漏少量弯曲特征，进而导致局部地形特征线在连接时出现中断。为此文献^[16]预先对满足一定条件的等高线进行插值处理，再采用联合DT的方法提取地形特征线，从而解决大部分地形特征线中断的问题。但是这样获取的地形特征线很不平滑，局部呈现锯齿状，给地形特征线的评价增加了难度和不准确因素。

针对上述问题，本文提出了一种基于特征段CDT的地形特征线生成方法。

基于特征段CDT的地形特征线生成方法在对等高线划分特征段的基础上生成特征段CDT，然后利用特征段CDT获得局部地形特征线树，进而确定地形特征点，并利用特征段CDT进行地形特征线的自动追踪，获取地形特征线。

图 1 总体思路示意图 Fig. 1 Schematic diagram of the general idea

图选项

特征段又称凹凸段、特征弯曲，是指等高线上对应着分水地带(山脊)与汇水地带(山谷)的弧段。特征段划分一般采用矢量叉积法^{[2, 18]}。图 2为采用此法获得的特征段示意图。

图 2 凹凸点构成的特征段 Fig. 2 The feature segment constituted by concave and convex points

图选项

经过凹凸段划分获得的特征段是一些离散点的有序排列，而DT适合对大量不规则离散点构造Delaunay三角网，对少量点构造的三角网的规整性往往差强人意。不规则的三角网很容易丧失探测和表达弯曲特征的条件。本文改进传统等高线三角网的生成方法，提出了三角网生长法对特征段构建CDT。

三角网生长法先确定起始基线，然后通过一定的搜寻规则寻找第3点生长出第1个CDT，再以两条新边作为新的基线继续CDT的“生长”，直到该段上所有的点都参与构网。本文采用“最小角最大”原则作为寻找第3点的限制条件。需要注意的是凹段和凸段的弯曲方向相反，在实际操作中应有所区别。

基于特征段的CDT生成步骤为：

(1) 选取特征段第1个点P₁和第2个点P₂，连成有向线段P₁P₂作为初始基线。

(2) 沿基线固定一侧搜寻可能为第3点的候选点集C，寻找与基线两端点连线夹角最大点作为第3点，假设为M点，则有∠P₁MP₂最大。生成三角形P₁MP₂。在哪一侧搜寻取决于特征段的凹凸性，对于凹段，应沿基线P₁P₂方向右侧查找候选点集C；对于凸段，则沿基线P₁P₂方向左侧查找候选点集C。

(3) 分别以有向线段P₁M以及有向线段MP₂作为新基线，重复步骤(2)、(3)直到所有点处理完毕。

以图 3所示凸段为例，首先找到初始基线P₁P₂。本文采用矢量叉积法搜索位于基线左侧的点，对于第n个点P_n，如果P₁P₂×P₁P_n＞0(凹段叉积小于0)，则点P_n就位于有向线段P₁P₂的左侧。根据这一原理，P₃、P₄、…、P₁₅、P₁₆都位于P₁P₂的左侧，这时利用余弦定理找到与基线首尾点连线夹角α最大的点P₁₆。于是生成第一个Delaunay三角形P₁P₂P₁₆。然后分别以P₁P₁₆和P₁₆P₂为新的基线利用上述方法构建新三角形。另外需要注意的是如果上图为凹段，则其点的顺序跟凸段正好相反。

图 3 三角网生长法构建特征段CDT示意图 Fig. 3 Diagram about building CDT on feature segment using of triangulation growth method

图选项

使用上述方法生成的特征段CDT具有如下特点：

(1) 特征段CDT中三角形个数与特征段曲率成正比。

(2) 相邻等高线上属于同一山脊(或山谷)的特征段CDT具有形态相似性。

(3) 同一条等高线上的相邻特征段CDT在拐点处不相交。

特征段CDT的第1个特点表明特征段CDT与特征段曲率具有紧密的联系，这为利用特征段CDT描述特征段及通过特征段CDT获取地形特征点提供可靠性依据。第2个和第3个特点则为地形特征点匹配和地形特征线连接提供了一个简单易行的方法。

利用特征段CDT可以构建局部地形特征线树，局部地形特征线树相互连接构建完整的地形结构线树，可以用于地貌分析和水系流域分析中。

本文将基于每个特征段CDT构建的地形结构线称为局部地形特征线树。构建局部地形特征线树之前需对三角网中三角形进行分类。文献^[19]在基于Delaunay三角网提取多边形骨架线时将三角网中的三角形分为端点三角形、主干三角形和分支三角形三类，本文借鉴了这种分类方法。

构建局部地形特征线树的过程可以分为两个部分，第1个部分主要是处理初始三角形，第2个部分主要是处理初始三角形之外剩余三角形。初始三角形指包含初始边(由特征段首尾点相连接而构成的边)的三角形，因为其分类特点与其余三角形不同，故需单独处理。基于特征段三角网构建局部地形特征线树的步骤为：

(1) 在特征段CDT中找出初始三角形,不同弯曲形态中的初始三角形可归结为端点三角形和分支三角形两种类型。根据初始三角形的类型获取相应地形特征线片段，然后将其邻接三角形作为下一个要处理的三角形转入步骤(2)中。

(2) 将步骤(1)所得初始三角形的邻接三角形作为第1个三角形，利用Delaunay三角网中三角形的邻接关系依次获得剩余三角形，递归处理这些三角形。处理过程根据邻接三角形类型分为下面3种情况：① 如果为端点三角形，找到其邻接边，然后将邻接边所对顶点作为叶节点插入到当前节点下，结束递归，完成当前特征段局部地形特征线树的构建;② 如果为主干三角形，找到另外一条未被处理的邻接边，将该邻接边的中点作为当前节点的子节点插入到局部地形特征线树中，并标记当前主干三角形为已处理三角形，然后将未被处理的邻接三角形作为下一个待处理三角形进行递归处理;③ 如果为分支三角形，首先找到已被处理的邻接三角形所对应的邻接边，根据此边定位另外两条邻接边，分别将它们的中点作为当前节点的子节点插入到局部地形特征线树中，将它们对应的邻接三角形作为待处理三角形，并标记当前三角形为已处理三角形。

基于特征段CDT构建局部地形特征线树为特征点识别和地形特征线连接奠定了基础。局部地形特征线树的叶节点可以看做是相应特征段上的特征点。基于特征段的CDT虽然立足于特征段，但其探测出的特征点又不局限于特征段的个数，还与特征段形态有关。因此该方法克服了文献^[2]中提出的基于特征段的多因子特征提取法应用于闭合凸多边形等高线上的缺陷。

特征点指构成山脊线和山谷线的点。对于某些包含多个端点三角形的特征段CDT，基于其构建的局部地形特征线树会有多个叶节点，相应的特征段会具有多个特征点。然而一个特征段包含多个特征点不利于特征点的匹配(特征点匹配是指根据已知特征点找到属于同一地形特征线的下一个特征点)和特征线的连接。基于此，本方法对特征段进行再分。特征段再分的步骤为：

(1) 在原始特征段CDT中查找端点三角形的个数，如果小于2则不需再分，若大于等于2，进入步骤(2)。

(2) 对每一个端点三角形，将其作为新段CDT的端点三角形，寻找它的邻接三角形。如果邻接三角形为空、初始三角形或分支三角形三种中的任一种，则将端点三角形的三个顶点作为新段的点集，此处新段完成；如果端点三角形的邻接三角形不为空，且不为初始三角形或分支三角形，则将其邻接三角形的第3个顶点添加到新段的点集中，并将邻接三角形添加到新段的特征段三角网中，再寻找邻接三角形的邻接三角形，依次递归直到邻接三角形为空、初始三角形或分支三角形中的任一种。

(3) 将再分的新特征段添加进特征段集合中，并在特征段集合中删除原始特征段，完成本条等高线上所有特征段的再分。

图 4是根据上述方法再分特征段后获得的新特征段与原特征段三角网及特征点对比图(这里特征点采用多因子特征法提取)。通过对比可以看出，对于闭合的凸多边形等高线，特征段再分后的特征点提取效果要明显优于特征段再分之前。

图 4 特征段再分前、后特征点提取效果对比 Fig. 4 Comparison of feature point extraction before and after dividing feature segment

图选项

本文直接将特征段局部地形特征线树的叶节点作为特征点。在对特征段进行再分之后，每个特征段的局部地形特征线树有且仅有一个叶节点，即为本特征段上的特征点。该方法充分利用了特征段CDT对弯曲特征的探测功能，与多因子特征提取法相比具有一定的优越性。

图 5为在一幅等高线图上采用两种方法获取特征点的比较。通过对比可以发现，本文方法除了不受特征段个数影响外还有另外两个优势：第一，提取的细节比较详细和丰富，可以提取出比较细微的地貌结构上的特征点，如图中灰色椭圆标识的地方，这些特征点所在的特征段都是要综合的对象；第二，很少受特征段对称性的影响，反观多因子提取法在非对称性特征段上提取的特征点往往不能精确反映原地形特征，尤其是靠近图幅边缘的那些被截断的特征段，如图中黑色椭圆标识的地方。

图 5 两种方法提取的特征点对比 Fig. 5 Comparison of feature points extracted by two different methods

图选项

地形特征线追踪是指通过特征点匹配确定地形特征线上的点。特征点匹配是指利用地形特征线上的已知特征点判断待判点是否为此地形特征线上的点。目前的研究者多采用联合多个匹配因子进行特征点匹配。匹配因子越多，地形特征线方向的判断就越精确，但是随之也增加了匹配的复杂性。本文结合前文生成的特征段CDT使用新的匹配因子--特征段三角网取得了比较理想的效果。

相交原则。如果参考特征段CDT有且只有一个待判特征段与之相交，称之为第1类相交，则待判特征段即为目标特征段，目标特征段上的特征点即为匹配点，如图 6中b圆内的两个特征段，右边为参考特征段，其CDT与左边等高线上的一个待判特征段相交。如果相交的待判特征段不止一个，称之为第2类相交(见图 6中c圆内所示特征段)，若无相交待判特征段(见图 6中a圆内所示特征段)则使用距离和角度原则。

图 6 相交原则示意图 Fig. 6 Diagram of the intersection principle

图选项

距离原则。取距离参考特征点最近的待判特征点为目标特征点。使用距离原则须同时考虑角度原则。

角度原则。参考特征段CDT初始三角形外边两端点分别与待判特征点连线的夹角不小于某一阈值。经过试验与比较本方法采用阈值为5°时，可以最大限度地剔除不合理特征点，保留真实特征点。

不跨越原则。参考特征点与匹配点的连线不跨越等高线及已有地形特征线。

上述原则中，相交原则和距离原则用于判断待判特征点是否为目标特征点，而角度原则和不跨越原则可以剔除不合理的目标特征点。利用上述原则进行特征点匹配的步骤为：

(1) 在待判特征段集合中查找是否有与参考特征段相交的待判特征段，若有则使用相交原则。对于第1类相交直接将目标特征段的特征点作为与参考特征点相匹配的点；对于第2类相交则将与参考特征段CDT相交的所有待判特征段构成一个集合，在该集合中利用距离原则查找匹配特征点。

(2) 对不适用相交原则的情况，即没有与参考特征段CDT相交的待判特征段或属于第2类相交时，使用距离原则。找出距参考特征点最近的待判特征点作为目标特征点，然后对目标特征点使用角度原则，如果它们构成的角度小于阈值则该参考特征点无匹配的特征点，如果构成的角度不小于阈值则目标特征点即为匹配点。

(3) 对上述两步确定的匹配点使用不跨越原则。如果参考特征点与匹配点的连线(指参考特征点与匹配点之间的部分，不包括端点)与任一条等高线或已有的地形特征线相交则匹配点为伪匹配点，该参考特征点无匹配点。在查看参考特征点与匹配点的连线是否有与之相交的等高线时只查看参考特征点与匹配点所在的等高线即可。

等高线树是对等高线的结构化，它完美地表达了等高线的层次结构。将地形特征线的追踪与连接建立在等高线树的基础上可以提高追踪的效率。利用等高线树辅助地形特征线追踪须遵循下述规则：

(1) 山脊线的追踪适合从等高线树的根节点，即高程最低的等高线开始沿树干往下追踪，父节点等高线上的山脊点与子节点等高线上的山脊点具有多对一的关系。

(2) 山谷线的追踪适合从等高线树的叶节点，即高程最高(相对于树的每个分支来说)的等高线开始沿分支往上追踪，子节点等高线上的山谷点与父节点等高线上的山谷点具有多对一的关系。

(3) 无论是山脊线还是山谷线的追踪都必须沿树的每个分支往下或往上进行，不能跨越分支，从一个枝杈追踪到另一个枝杈。

明确上述规则之后，根据上节提出的特征点匹配原则进行地形特征线的追踪，主要步骤如下(以山脊线为例)：

(1) 从等高线树的根结点N₀开始，找出根节点等高线上的特征段集C₀。

(2) 从C₀中取出一个未处理过的凸段(凸段对应山脊)S，将S上的特征点P作为山脊线L的首点。在N₀的子节点N中寻找与S特征点相匹配(使用上节提出的匹配原则和步骤进行匹配)的凸段S′，如果S′已被处理过一次，即特征点P′也属于另一条特征线则将其特征点P′存进L中，并结束山脊线L的追踪。如果S′没有被处理过则将其特征点P′存进L中，并将S′标记为已处理。然后将N看作N₀，N的子节点看作N重复上述步骤，直到N₀为叶节点或无满足距离原则、角度原则和不跨越原则的特征段(点)时结束地形特征线L的追踪。如果N₀有多个子节点，则将所有子节点的特征段集合合并成一个更大的集合，在合并后的集合中利用上述方法寻找匹配特征段(点)。

(3) 从C₀中取下一个未处理的凸段C₁按步骤(2)进行处理。重复步骤(3)直到C₀中没有未被处理的凸段为止，以等高线根结点N₀开始的所有山脊线提取完毕。

(4) 获得N₀的子节点N，找出子节点N所表示的等高线上的特征段集C₁，将C₁看作C₀按步骤(2)、(3)步进行处理。将N看作N₀重复步骤(4)直到N₀无子节点为止，结束所有山脊线的提取。

谷底线的追踪与连接和山脊线相似，只是谷底线的追踪需从等高线树的叶节点开始，这里不再赘述。图 7为地形特征线追踪和连接的效果图。

图 7 地形特征线效果图 Fig. 7 Diagram of the terrain structure line

图选项

采用本方法获取的地形特征线具有以下几个优点：

(1) 能有效地提取地形特征，提取的地形特征线与实际地形的符合程度较高。

(2) 通过特征段CDT构建局部地形特征线树，进而获取地形特征点，并利用CDT作为特征点匹配的因子，这一过程基本不受地形复杂程度的影响。

(3) 提取的地形特征线比较完整，不需要进行二次连接。

(4) 提取的局部区域地形特征线自动构成树结构，方便地形特征线结构化。

(5) 提取的地形特征线主体部分比较平滑，有利于地形特征线的评价。

(6) 能提取出比较丰富的细节，较好地保留了细部，为地形图制图综合提供了详实的数据基础。

(7) 自动化程度较高，仅需输入角度阈值等参数便可由计算机自动完成地形特征线的追踪。

本文提出的基于特征段CDT的地形特征线提取方法利用CDT对等高线弯曲特征的探测功能，较好地解决了已有方法在特征点识别中的不足，根据特征段CDT追踪和连接地形特征线也取得了一定的效果，并通过试验证明了该方法用于地形特征线提取的合理性和准确性。利用本方法提取出的地形特征线可以直接进行地貌识别、流域分析，以及地貌结构化综合。该方法在角度原则阈值的选取上存在一些不足，即没有一个确切的方法，仅通过多次对比试验以及个人经验进行确定，还有待改进。