一种基于等高线的地形特征线提取方法

张 尧，樊 红，李玉娥

1.四川省基础地理信息中心，四川 成都 610041；2.武汉大学 测绘遥感信息工程国家重点实验室，湖北 武汉 430079；3.四川智图信息技术有限公司，四川 成都 610041

1 引 言

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

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

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

2 基本思路

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

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

3 具体方法

3.1 特征段划分

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

3.2 基于特征段构建约束Delaunay三角网

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

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

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

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

（1）选取特征段第1个点P1和第2个点P2，连成有向线段P1P2作为初始基线。

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

（3）分别以有向线段P1M以及有向线段MP2作为新基线，重复步骤（2）、（3）直到所有点处理完毕。

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

图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可以构建局部地形特征线树，局部地形特征线树相互连接构建完整的地形结构线树，可以用于地貌分析和水系流域分析中。

3.3 基于特征段三角网构建局部地形特征线树

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

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

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

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

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

3.4 特征点的确定

特征点指构成山脊线和山谷线的点。对于某些包含多个端点三角形的特征段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

3.5 地形特征线自动追踪

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

3.5.1 特征点匹配原则

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

3.5.2 地形特征线追踪与连接

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

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

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

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

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

（1）从等高线树的根结点N0开始，找出根节点等高线上的特征段集C0。

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

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

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

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

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

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

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

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

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

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

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

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

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

4 结 论

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

［1］ WU Yanlan.Research on Map Algebra Model and Methods of 3DRelief Generalization［D］.Wuhan：Wuhan University，2004：8－17.（吴艳兰.地貌三维综合的地图代数模型和方法研究［D］.武汉：武汉大学，2004：8－17.）

［2］ ZHANG Linlin.Research on Structured Method of Relief Generalization［D］.Zhengzhou：The Information Engineering University，2005：37－41.（张琳琳.地貌自动综合的结构化方法研究［D］.郑州：信息工程大学，2005：37－41.）

［3］ CHEN Haiyan，WAN Gang.Research on the Algorithm of Automatic Building Topgraphic Structures Using the Contour Data［J］.Bulletin of Surveying and Mapping，2003（3）：21－23.（陈海燕，万刚.利用等高线数据自动生成地性结构线的算法研究［J］.测绘通报，2003（3）：21－23.）

［4］ MAO Kebiao，CHEN Xiangdong.A Research on Constructing Terrain Features Automatically［J］.Science of Surveying and Mapping，1995（3）：12－18.（毛可标，陈向东.地形结构线自动生成方法研究［J］.测绘科学，1995（3）：12－18.）

［5］ SONG Yonghak，SHAN Jie.An Adaptive Approach to Topographic Feature Extraction from Digital Terrain Models［J］.Photogrammetric Engineering and Remote Sensing，2009，75（3）：281－290.

［6］ LIU Huijie，JIN Hailiang，MIAO Baoliang.An Algorithm for Extracting Terrain Structure Lines Based on Contour Data［C］∥Proceedings of SPIE－The International Society for Optical Engineering.Beijing：［s.n.］，2010.

［7］ CHEN X D.Automatic Building Topographic Structures Lines Using the Digitized Contour Data［C］∥Proceedings of SPIE－The International Society for Optical Engineering.Beijing：［s.n.］，2010.

［8］ KWEON I S，KANADE T.Extracting Topographic Terrian Features from Elevation Maps［J］.Computer Vision，Graphics，and Image Processing：Image Understanding，1994，59（2）：171－182.

［9］ ZHU Qing，ZHAO Jie，ZHONG Zheng.The Extraction of Topographic Patterns Based on Regular Grid DEMs［J］.Acta Geodaetica et Cartographica Sinica，2004，33（1）：77－82.（朱庆，赵杰，钟正.基于规则格网DEM的地形特征提取算法［J］.测绘学报，2004，33（1）：77－82.）

［10］ KWEON I，KANADE T.Extracting Topographic Terrain Features from Elevation Maps［J］.CVGIP：Image Understanding，1994，59（2）：171－182.

［11］ BARKOWSKY T，LATECKI L J，RICHTER K F.Schematizing Maps：Simplification of Geographic Shape by Discrete Curve Evolution.Spatial Cognition II－Integrating Abstract Theories，Empirical Studies，Formal Models，and Practical Applications［M］.Berlin：Springer，2000.

［12］ PRICE CURTIS V，WOLOCK DAVID M，AYERS MARK A.Extraction of Terrain Features from Digital Elevation Models［C］∥Proceedings of 1989National Conference Hydraulic.London：［s.n.］，1989：845－850.

［13］ AI Tinghua，GUO Renzhong，LIU Yaolin.A Binary Tree Representation of Curve Hierarchical Structure in Depth［J］.Acta Geodaetica et Cartographica Sinica，2001，30（4）：343－348.（艾廷华，郭仁忠，刘耀林.曲线弯曲深度层次结 构 的 二 叉 树 表 达 ［J］.测 绘 学 报，2001，30（4）：343－348.）

［14］ GUO Qingsheng，YANG Zuqiao，FENG Ke.Extracting Topographic Characteristic Line from Contours［J］.Geomatics and Information Science of Wuhan University，2008，33（3）：253－256.（郭庆胜，杨族桥，冯科.基于等高线提取地形特征线的研究［J］.武汉大学学报：信息科学版，2008，33（3）：253－256.）

［15］ WU Fan，SU Weimin，YANG Yingwei，et al.Extracting Terrain Features from Contour Maps Based on United－Delaunay－Triangulaltion Model ［J］.Journal of China University of Minig and Technology，2007，36（2）：172－175.（吴凡，粟卫民，杨英伟，等.基于联合Delaunay三角网的等高线地形特征提取研究［J］.中国矿业大学学报，2007，36（2）：172－175.）

［16］ YANG Yingwei，LUO Juan.Research on Contour Group Generalization Based on United Delaunay Triangulation Model［J］.Geospatial Information，2010，8（2）：1－9.（杨英伟，罗娟.联合Delaunay三角网的等高线群综合研究［J］.地理空间信息，2010，8（2）：1－9.）

［17］ AI Tinghua，ZHU Guorui，ZHANG Genshou.Extraction of Landform Features and Organization of Valley Tree Structure Based on Delaunay Triangulation Mode［J］.Journal of Remote Sensing，2003，7（4）：292－299.（艾廷华，祝国瑞，张根寿.基于Delaunay三角网模型的等高线地形特征提取及谷地树结构化组织［J］.遥感学报，2003，7（4）：292－299.）

［18］ LIU Ying.The Representation Identification and Generalization of Spatial Graphics［D］.Zhengzhou：The Information Engineering University，2005：100－102.（刘颖.空间图形的表达、识别与综合［D］.郑州：信息工程大学，2005：100－102.）

［19］ FAN Hong.On Automatic Map Labeling［M］.Beijing：Publishing House of Surveying and Mapping，2004：81－82.（樊红.地图注记自动配置的研究［M］.北京：测绘出版社，2004：81－82.）