Print

发布时间: 2017-11-16
摘要点击次数:
全文下载次数:
DOI: 10.11834/jig.170252
2017 | Volume 22 | Number 11




    地理信息技术    




  <<上一篇 




  下一篇>> 




显著度可控的DEM地形特征线提取
expand article info 邹昆1,2, 翁宏章2, 李文生1,2
1. 电子科技大学中山学院计算机学院, 中山 528402;
2. 电子科技大学计算机科学与工程学院, 成都 611731
收稿日期: 2017-06-21; 修回日期: 2017-08-20
基金项目: 国家自然科学基金项目（61502088）；广东省自然科学基金项目（2016A030313018）；广东省高等学校优秀青年教师培养计划项目（YQ2015241，YQ2015242）
第一作者简介: 邹昆(1980-), 男, 副教授, 2008年于华南理工大学获计算机应用技术专业博士学位, 主要研究方向为数字图像处理、纹理合成、地形合成。E-mail:cszoukun@foxmail.com.
中图法分类号: TP391
文献标识码: A
文章编号: 1006-8961(2017)11-1611-12

摘要

目的 基于数字高程模型（DEM）的地形山脊线和山谷线提取对地形模型简化、基于样本的地形合成和地形地貌研究有重要意义，针对许多传统算法无法对所提取特征线的显著度进行方便准确的控制，以及不支持环形特征线提取的问题，提出一种新的显著度可控的DEM地形特征线提取算法。方法 首先利用全局断面扫描算法提取特征点并计算各特征点的显著度，然后根据特征点的特征方向进行特征延伸以增强特征连通性，接着采用改进的Hilditch细线化算法对特征点集合进行细线化处理，之后为相邻特征点添加特征边，构成特征图，利用环路检测与破环算法检测特征图中的环路，并破除冗余小环路，最后根据分支显著度的相似度和分支方向一致性进行特征图分解，计算分解得到特征线的显著度并筛选得到最终特征线。结果 使用真实DEM数据提取最显著的若干条特征线，与现有的基于特征显著度的地形特征线提取算法进行对比，本文算法对特征图的分解能够更准确地提取主干特征线，而基于显著度的特征线筛选控制也更加准确合理。对提出的环路检测与破环算法进行实验验证，该算法能保留大的山脊线环路，破除小的冗余环路。结论 实验结果表明，本文算法能有效实现显著度可控的山脊线和山谷线自动提取，提取结果与人眼观察结果基本一致，同时能够支持含有环形特征的地形。

关键词

地形特征线提取; 特征显著度; 数字高程模型; 特征图分解; 特征线筛选

Saliency-controllable topographic feature line extraction from DEM
expand article info Zou Kun1,2, Weng Hongzhang2, Li Wensheng1,2
1. School of Computer Engineering, Zhongshan Institute, University of Electronic Science and Technology of China, Zhongshan 528402, China;
2. School of Computer Science and Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China
Supported by: National Natural Science Foundation of China (61502088)

Abstract

Objective Ridge lines and valley lines are the most important terrain feature lines, which delineate skeleton structures of terrain.The automatic extraction of them from digital elevation models (DEM) is of great significance in applications such as relief automated generalization, hydrological analysis and geographical information system, and can also provide matching features for terrain model simplification and sample-based terrain synthesis.For most traditional algorithms, feature lines can only be filtered based on their length by eliminating endmost short branches, and it is difficult to accurately control the saliency of extracted feature lines.Because feature lines are usually modeled as trees, and break-all-loop algorithms such as polygon breaking and minimum spanning tree are commonly used, ringlike features such as craters cannot be correctly extracted without notches.Therefore, a new algorithm for saliency-controllable topographic feature line extraction from DEM was proposed based on feature saliency. Method Our algorithm consisted of five steps.First, feature points were extracted by the global profile scan algorithm which can better eliminate noise and pseudo feature points, and feature saliency of each feature point was calculated according to its height, ambilateral height drops and slopes.Second, feature connectivity was enhanced by feature extension based on feature directions of feature points.Third, feature point set was thinned to be one-pixel-wide with an improved Hilditch algorithm, in which a neighborhood scan process was added when finding a point that can be deleted, and among all points satisfying remove conditions in the neighborhood, the one with the smallest feature saliency was chosen to be deleted.In such a way, thinning result had almost nothing to do with the global scanning sequence, and salient feature points could be preserved.Fourth, feature graph was constructed by taking feature points as nodes and adding edges to all neighboring ones.It was highly likely that many loops existed in the graph, and most of them were redundant small ones.In order to support extraction of ringlike features, we broke small loops while preserving big ones.A loop detecting and breaking algorithm was devised to handle all kinds of loops, such as nested loops, loops containing common sides and loops with no interior points.In the end, feature graph was decomposed into branches and feature lines were constructed by assembling consecutive branches, requiring that all branches inside a feature line had similar feature saliencies and directions of neighboring branches were consistent.Then, feature saliency of each feature line was calculated, based on which all feature lines were sorted in descending order.And by selecting a certain number of headmost feature lines in the sorted list, it was simple to implement feature saliency-based filtering. Result To validate the effectiveness of the proposed algorithm, a comparison with the existing feature saliency-based topographic feature line extraction algorithm was made by extracting some most significant feature lines from real DEM data.The proposed algorithm can extract feature line trunks more accurately and can effectively avoid overextension of trunks, which is attributed to the new feature graph decomposition strategy.In terms of saliency-based control, by simultaneously considering feature line length, average and maximum feature saliency of all feature points on the line, the new formula of feature saliency calculation for feature lines is more reasonable, and so the proposed algorithm can provide more accurate control over the saliency of extracted feature lines.To verify the proposed loop detecting and breaking algorithm, intermediate results before and after loop breaking were shown, and a synthesized DEM with a big crater was tested. Results show that the proposed algorithm can well preserve the large ridge line loop while breaking all small redundant ones.Because differences in thinning results for real DEM data are at pixel level and can hardly be perceived by eye, so a comparison was made based on a man-made example of 2-pixel-wide feature point band to exemplify the proposed improved thinning algorithm.And the improved algorithm deletes feature points with smaller feature saliency which means true features are more likely to be preserved. Conclusion The proposed algorithm can effectively extract ridge and valley lines with controllable saliency, which are consistent with human eye observation, and can support terrains with ringlike features.The main contributions include three points.First, by exploiting the proposed feature graph decomposition algorithm and applying the devised formula to calculate feature saliency of feature lines, decomposition results become more reasonable, and more precise control over the feature saliency of extracted feature lines can be provided.Second, a novel algorithm was proposed to detect all redundant small loops while preserving large ones, and it can handle all complex loop cases.At last, a further improvement upon the improved Hilditch algorithm was made based on feature saliency, so that more salient feature points will be retained and the resulting skeleton lines may be more accurate for feature point bands of 2 or 3 pixel wide.However, there are several shortcomings of the proposed algorithm.First, there are many parameters which are set mainly based on experience.Especially, the branch feature saliency dissimilarity threshold used in feature graph decomposition should be carefully tuned, which has a great influence on experimental results. Second, compared with some break-all-loop algorithms such as minimum spanning tree, the loop detecting and breaking algorithm is inefficient and it becomes a bottleneck of the whole algorithm.Future work includes improving the efficiency of the loop detecting and breaking algorithm as well as parameter optimization and automatic parameter setting.

Key words

topographic feature line extraction; feature saliency; digital elevation model; feature graph decomposition; feature line filtering

0 引言

山脊线和山谷线描述了地形的骨架结构，是最重要的地形特征线，其自动提取在自动地貌综合、水文分析、地理信息系统等方面有着非常重要的应用，同时也可为地形模型简化[1-2]和基于样本的地形合成[3-4]提供重要的特征依据。根据数据源的不同，可分为基于等高线的方法[5]和基于数字高程模型(DEM)的方法[6-18]，后者由于数据形态与图像相似以及应用的相关性，除测绘学界外也受到了图形图像领域学者的广泛关注。

现有的基于DEM的地形山脊线和山谷线提取算法可分为三大类：基于图像处理技术检测曲线结构的算法[6-9]、基于地表流水物理模拟的算法[10-11]以及基于地表几何形态分析的算法[12-18]。基于图像处理技术检测曲线结构的算法将检测一般图像中曲线结构的方法应用于DEM特征线检测，在线形特征不明显的山脊/山谷分支连接处易产生断裂。基于地表流水物理模拟的算法通过模拟地表水流过程，确定各点水流方向和汇水量，将汇水量大于阈值的特征点连接形成山谷线，将汇水区域的边界作为山脊线，此类算法提取的特征线整体性很好，但提取的山脊线为闭合曲线，不符合实际，此外较难判定平坦区域的水流方向。基于地表几何形态分析的算法利用山脊和山谷点为相应方向上的高度极值点来提取地形特征线，此类算法简单高效，提取精度较高，但容易受噪声干扰，特征连续性较差，此外提取结果会受到检测窗口大小的影响。

近年来，基于地表几何形态分析的算法取得了较大进展。Chen等人[12]将纵、横及两个对角线方向上的局部断面极值点作为特征点，然后通过细线化筛选特征点，最后通过剪枝等后处理得到最终特征线，该算法简单易行，但所提取特征线的连通性不佳。Huang等人[13]通过计算梯度方向来选择最优断面方向，进一步完善了地形断面极值法。Chang等人[14-15]提出了剖面识别和多边形拆解算法(PPA)，通过降低特征点的入选条件很好地保证了特征线的连续性，但其多边形拆解算法复杂度太高，此外由于多边形拆解会破开所有环路，导致PPA算法不支持环形山、火山口等环形地形特征的提取。Bangay等人[16]通过将多边形拆解算法替换为最小生成树算法，大幅提升了特征提取速度。Zhang等人[17]在PPA的基础上进行改进，提出了基于形态学和几何剖面识别的骨架特征提取算法，通过使用基于形态学的拆解，仅拆除不稳定环路，保留稳定的环路，从而支持环形地形特征的提取，同时也提高了效率。然而上述算法均存在以下两方面问题：一是不能够对所提取特征线的显著程度进行方便的筛选控制；二是在识别特征点时均使用局部窗口，信息的不全面导致易受噪声和伪特征影响，窗口大小的设定也对结果有较大影响。Zou等人[18]提出一种基于特征显著度的地形特征线提取算法，利用4个主要方向的全局断面扫描来识别特征点，通过最小生成树算法破开环路，最后基于特征显著度对特征线进行筛选控制。该算法较好地解决了上述两个问题，但其在将特征树分解成特征线时容易造成主干特征线的过度延伸或是提取不准确，同时其特征线显著度定义方式不是十分合理，造成对特征线显著度的控制不够准确，此外该算法也不支持环形地形特征的提取。

为了实现对所提取特征线的显著度进行更准确的控制，同时支持环形地形特征提取，本文在前期工作[18]的基础上提出一种改进的显著度可控的地形特征线提取算法。该算法采用了文献[18]中的基于全局断面扫描的特征点提取及特征延伸方法，但后续处理皆不相同。一方面，借鉴Zhang等人算法[17]的思想，将基于最小生成树的破环算法替换为细线化和基于区域填充的环路检测和破环算法，以支持环形地形特征提取。与文献[17]算法的区别在于：采用了不同的细线化算法，并在细线化时加入特征显著度的考虑，以使得特征位置更准确；此外采用新的环路检测与破环算法，使得能够同时处理所有环路，且能够应对嵌套环路等复杂环路情况。另一方面，采用了不同的特征图分解算法以及不同的特征线显著度计算方法，使得基于显著度的筛选控制更准确。实验结果表明，与现有算法相比，该算法对特征线显著度的控制更符合人眼观察实际，同时能够支持含有环形特征的地形。

1 本文算法

本文地形特征线提取算法的流程如图 1所示(以山脊线提取为例)。



图 1 本文算法流程
Fig. 1 The workflow of the proposed algorithm

算法详细步骤如下：

1) 特征点提取。采用全局断面扫描算法提取特征点，同时计算各特征点的特征显著度。

2) 特征延伸。根据特征点的特征方向进行特征延伸，以增强特征连通性。

3) 细线化。对特征点集合进行细线化处理，使其所构成的特征点带为单像素宽。

4) 环路检测与破环。为相邻特征点添加特征边，构成特征图；采用环路检测与破环算法检测特征图中的所有环路，并破除小面积冗余环路(由于破环过程只是删除连接相邻特征点的边而未删除特征点，若未对DEM进行下采样，难以看出差别，故图 1中未显示破环结果)。

5) 特征图分解与特征线筛选。根据相邻分支之间显著度的相似度以及分支方向的一致性将特征图分解为一条条的特征线，计算所得到的特征线的显著度，并基于特征线显著度进行特征线筛选。

由于可通过将DEM数据取反然后用提取山脊线的方法提取山谷线，下面以山脊线提取为例进行算法描述。

1.1 特征点提取与特征延伸

采用与文献[18]相同的方法提取特征点和进行特征延伸：

1) 采用全局断面扫描算法提取山脊点。分别沿水平、垂直及两个斜45°方向对DEM数据进行逐行扫描，根据前后地形走势、高度差等信息提取出山脊点。由于利用了全局信息，更容易剔除噪声和伪特征点。在确定山脊点之后，利用扫描过程中获得的信息，根据高度、两侧高度落差和两侧坡度计算每一山脊点${{p_i}}$的特征显著度[18]

$ \begin{array}{l} \begin{array}{*{20}{l}} {S\left( {{p_i}} \right) = } \end{array}\\ {k_{\rm{h}}}\frac{{{h_i}}}{{{h_{{\rm{max}}}}}} + {\rm{ }}{k_{\rm{d}}}\frac{{|\Delta h_i^{\rm{l}}| + \Delta |h_i^{\rm{r}}|}}{{2\Delta {h_{{\rm{max}}}}}} + \\ \frac{{{k_{\rm{s}}}}}{2}\left( {\frac{{|\Delta h_i^{\rm{l}}|}}{{\sqrt {\Delta h{{_i^{\rm{l}}}^2} + {\rm{ }}d{{_i^{\rm{l}}}^2}} }} + \frac{{|\Delta h_i^{\rm{r}}|}}{{\sqrt {\Delta h{{_i^{\rm{r}}}^2} + {\rm{ }}d{{_i^{\rm{r}}}^2}} }}} \right) \end{array} $ (1)

式中，${{h_i}}$为${{p_i}}$的高度，Δ ${h_i^{\rm{l}}}$和Δ ${h_i^{\rm{r}}}$为${{p_i}}$在相应扫描方向上左右两边的高度落差，假设达到此高度落差的最近点为${p_i^{\rm{l}}}$和${p_i^{\rm{r}}}$，则${d_i^{\rm{l}}}$和${d_i^{\rm{r}}}$分别为${{p_i}}$与${p_i^{\rm{l}}}$和${p_i^{\rm{r}}}$的水平距离，${{h_{{\rm{max}}}}}$为最大高度，Δ${{h_{{\rm{max}}}}}$为最大高度落差，$k_\rm{h}$、$k_\rm{d}$和$k_\rm{s}$为权重系数。如果${{p_i}}$在多个方向的扫描中均被检测为特征点(即计算了多个$S$ (${{p_i}}$)值)，仅保留最大$S$(${{p_i}}$)值。

2) 特征延伸。对于每一山脊点尝试沿其特征方向及特征方向的相邻方向进行延伸(特征方向与断面扫描方向垂直)，如果当前考察点的高度不小于两侧相邻点高度，则将其纳入新增山脊点，新增山脊点的显著度为0(以便于后面筛选时将主要由新增山脊点构成的山脊线筛选掉)。

1.2 细线化

初步得到的山脊点中可能存在孤立点，另外由于沿多个方向进行断面扫描，在部分显著山脊点的两侧存在不同特征方向的显著度较弱的相邻山脊点，从而构成了局部位置上的山脊点宽带(通常为2 3像素宽)。孤立点需要剔除，而山脊点宽带需要进行细线化处理。

Chen等人[12]采用改进的Hilditch细线化算法来完成此任务(具体形式化描述可参见文献[12])：

1) 按扫描线顺序依次考察所有山脊点，将同时满足以下5个条件的点标记为待删除点(如无特别说明，在测试这些条件时将已标记的待删除点仍视为山脊点)：

(1) 当前点为山脊点；

(2) 当前点位于山脊点宽带边缘(4邻域内至少有1个非山脊点)；

(3) 当前点不是端点(4邻域内至少有2个山脊点)；

(4) 删除当前点后3×3邻域内的山脊点仍然是8连通的；

(5) 线宽为2的部分只消除一侧(对于8邻域内已标记的每一待删除点，将其视为非山脊点时均仍满足条件(4)，若8邻域内不存在已标记的待删除点，同样视为满足本条件)。

2) 将所有标记为待删除的点从山脊点集合中删除，如果未删除任何点则转第3) 步，否则转步骤1)。

3) 删除所有孤立山脊点。

然而由于该算法沿扫描线标记待删除点，会优先删除先扫描到的满足条件的山脊点，而并未考虑构成局部宽带的山脊点之间的显著度差异。

本文在其基础上进行改进，对步骤1) 进行修改：仍然按扫描线顺序考察每一山脊点$p$，但如果当前山脊点$p$满足上述5个删除条件，查看其5×5邻域内是否有其他山脊点满足删除条件，将邻域内所有满足删除条件的山脊点中特征显著度最小的山脊点标记为待删除点，继而考察扫描线顺序上$p$的下一个未被标记为待删除点的山脊点，其他处理与文献[12]一致。这样对于线宽为2或3的山脊线，会优先删除特征显著度小的山脊点。

图 2给出了改进前后的细线化算法针对一个双像素宽特征点带示例的效果对比，其中空白格代表非特征点，方格中央非0数字表示特征点的特征显著度，0表示细线化算法删除的特征点，方格右下角的数字表示在扫描过程中被标记为待删除点的顺序。经过一次扫描，Chen算法[12]将优先扫描到的满足删除条件的6个特征点删除了(剩下2个点不满足步骤1) 中条件(5))，而本文算法则在一次扫描过程中删除了特征显著度较低的4个点(依次为考察显著度为10、5、7、6的点时标记为待删除点的)。此外，由于多像素宽特征点带两端的像素并非“端点”(满足条件(3))，Chen算法[12]会缩短此类特征点带的长度，而本文算法可一定程度减少此现象的发生。实际上，对于双像素宽特征线，3×3邻域已足够，采用5×5邻域主要是为了使算法也适用于3像素宽的特征点带。



图 2 双像素宽特征点带细线化对比
Fig. 2 Comparison of thinning for a 2-pixel-wide feature point band ((a) an example of 2-pixel-wide feature point band; (b) Chen[12] result; (c) ours)

1.3 环路检测与破环

假定DEM为$M$行$N$列，在细线化完成后，得到$M$行$N$列的特征点标记图$\mathit{\boldsymbol{L}}$，其第$i$行第$j$列元素为

$ L\left( {i, j} \right) = L\left( {{p_{i, j}}} \right) = \left\{ \begin{array}{l} 0{\rm{ }}\;\;\;{p_{i, j}}为特征点\\ 1{\rm{ }}\;\;\;{p_{ij}}非特征点 \end{array} \right. $

根据$\mathit{\boldsymbol{L}}$对相邻特征点进行全连接操作，形成特征边，如果出现边交叉的情况，剔除显著度较低的边(特征边的显著度定义为构成该边的特征点显著度之和)。此时得到的特征图中通常存在环路，其中绝大多数为面积非常小的环路，尤其是小三角形，这些小环路应该被破除，而由于环形山脊等环形特征在自然界及合成地形中广泛存在，需要保留较大的环路。Zhang等人[17]分开处理冗余小三角形和其他环路，本文提出一种新的环路检测与破环算法来同时处理所有环路，且能够应对各种复杂环路情况，例如嵌套环路、共边环路、共边嵌套环路等。

算法思想是：将特征点标记图$\mathit{\boldsymbol{L}}$放大至3倍大小，使得任何环路均存在内部元素，然后检测由特征线环路划分得到的多个4连通区域，为不同区域赋予不同的区域标记，从而得到区域标记图$\mathit{\boldsymbol{R}}$ (图 3给出了一个示例)；在此基础上，根据特征边左右两侧元素的区域标记来区分特征边的类型(非环路边、一般环路边或环路公共边)，并检测环路；最后根据环路面积大小选择保留或破除环路，在破环时，通过分析构成环路的特征边的情况来获取环路类型信息(简单环路、内部嵌套环路或是与其他环路有公共边的环路)，并进行有针对性地破环。



图 3 复杂环路示例
Fig. 3 An example of complicated loops ((a) feature point labels $\mathit{\boldsymbol{L}}$; (b) region labels $\mathit{\boldsymbol{R}}$)

环路检测与破环算法采用如下符号约定：

1) $R$ ($i$，$j$)为3$M$行3$N$列的区域标记图$\mathit{\boldsymbol{R}}$中第$i$行第$j$列的元素，即

$ R\left( {i, j} \right) = R\left( {{p_{i, j}}} \right) = \left\{ \begin{array}{l} 0{\rm{ }}\;\;\;\;{p_{i, j}}位于特征线上\\ id > 0{\rm{ }}\;\;\;\;{p_{i, j}}在标记为id的区域内\\ - 1{\rm{ }}\;\;\;\;{p_{i, j}}状态不确定 \end{array} \right. $

2) $\mathit{\boldsymbol{E}}$为特征边集合，对于集合中的特征边$e$，$p$$_e$ $^{{\rm{left}}}$和$p$$_e$$^{{\rm{right}}}$表示$e$在$\mathit{\boldsymbol{R}}$中对应位置左右两侧的元素，如果$\mathit{\boldsymbol{R}}$ ($p$$_e$$^{{\rm{left}}}$)≠ $\mathit{\boldsymbol{R}}$ ($p$$_e$$^{{\rm{right}}}$)，则$e$为环路边；${v_e}$为访问标识，即

$ \begin{array}{*{20}{l}} {{v_e} = \left\{ \begin{array}{l} 0{\rm{ }}\;\;\;\;e为非环路边\\ 1{\rm{ }}\;\;\;\;e为一般环路边\\ 2{\rm{ }}\;\;\;\;e为环路公共边 \end{array} \right.} \end{array} $

3) $\mathit{\boldsymbol{A}}$ [$id$]表示$\mathit{\boldsymbol{R}}$中标记为$id$的区域的面积(内部元素个数)。

环路检测与破环算法步骤描述如下：

1) 构建3$M$行3$N$列的区域标记图$\mathit{\boldsymbol{R}}$，将其所有元素值均初始化为-1；

2) 遍历特征边集合$\mathit{\boldsymbol{E}}$，对于其中每一条特征边$e$，假定其两个端点的坐标分别为($i$$_1$，$j$$_1$)和($i$$_2$，$j$$_2$)，其在$\mathit{\boldsymbol{R}}$中对应的坐标点为(3$i$$_1$，3$j$$_1$)和(3$i$$_2$，3$j$$_2$)，将$R$ (3$i$$_1$，3$j$$_1$)和$R$ (3$i$$_2$，3$j$$_2$)及$\mathit{\boldsymbol{R}}$中二者连线上的所有元素的区域标记置为0；

3) 初始化起始坐标$p$$_{{\rm{start}}}$为(0，0)，$id$=1，$\mathit{\boldsymbol{A}}$ [$id$]=0；

4) 从$p$$_{{\rm{start}}}$开始沿扫描线顺序遍历$\mathit{\boldsymbol{R}}$，寻找$R$ ($p$)=-1的元素$p$，找到则令$p$$_{{\rm{start}}}$=$p$，然后转步骤5)，找不到则转步骤6)；

5) 从$p$$_{{\rm{start}}}$开始，使用基于4邻域的深度优先搜索遍历区域标记与$p$$_{{\rm{start}}}$相同的连通区域，每遍历到一个符合条件的元素$p$，设$R$ ($p$)=$id$，并执行$\mathit{\boldsymbol{A}}$ [$id$]++，遍历结束后设$id$ =$id$+1，$A$ [$id$]=0，然后转步骤4)；

6) 遍历特征边集合$\mathit{\boldsymbol{E}}$，针对每一条特征边$e$，根据其方向获取其在$\mathit{\boldsymbol{R}}$中左右两侧的元素$p$$_e$$^{{\rm{left}}}$和$p$$_e$$^{{\rm{right}}}$，并考查以下条件：(1) $R$ ($p$$_e$$^{{\rm{left}}}$)≠ $R$ ($p$$_e$$^{{\rm{right}}}$)；(2) $R$ ($p$$_e$$^{{\rm{left}}}$) > 1且$R$ ($p$$_e$$^{{\rm{right}}}$) > 1；(3) $R$ ($p$$_e$$^{{\rm{left}}}$) > 1或$R$ ($p$$_e$$^{{\rm{right}}}$) > 1，若同时满足条件(1) 和(2)，则${v_e}$=2；若满足条件(1) 和(3) 但不满足条件(2)，则${v_e}$ =1；其他情况则${v_e}$=0；

7) 执行$id$=$id$-1，若$id$≤1，算法结束；

8) 若$\mathit{\boldsymbol{A}}$ [$id$]大于指定阈值，则不破除围绕$id$对应区域的环路，转步骤7)；

9) 环路检测。遍历特征边集合$\mathit{\boldsymbol{E}}$，寻找所有满足以下两个条件的特征边$e$：(1)${v_e}$ > 0；(2)$p$$_e$$^{{\rm{left}}}$== $id$或$p$$_e$$^{{\rm{right}}}$== $id$，将$e$加入环路边集合$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$，并执行${v_e}$--；

10) 破环。统计$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$中满足访问标识${v_e}$> 0的特征边数量$count$，如果$count$等于$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$中元素数量，意味着$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$构成的环路为内部嵌套环路，但无公共边，此时将$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$中所有特征边的${v_e}$置0，同时将$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$中显著度最小的特征边从集合$E$中删除；如果$count$小于$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$中元素数量，意味着$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$构成的环路为简单环路或是与其他环路有公共边的环路，此时从$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$中选择${v_e}$为0(即非公共边)且显著度最小的特征边从集合$\mathit{\boldsymbol{E}}$中删除。完成破环后清空$\mathit{\boldsymbol{E}}$$_{{\rm{loop}}}$并转第7) 步。

上述算法利用了DEM最外围无特征点的事实(由特征点检测算法所决定)，因此标记为1的区域一定是外围的非环内区域，从而可按步骤6) 中方法判断环路边的类型；此外，从$id$大的区域开始进行破环，是因为按照步骤4)5) 进行区域标记的方法，对于嵌套环路(无公共边)，越靠内的环路$id$越大，从内部环路开始破除，能够避开环形区域(例如图 3(b)中标记为2的区域)的检测以及同时破除双环的复杂性。图 3给出了一个复杂环路的示例，其中包含有无公共边嵌套环路、有公共边嵌套环路以及有公共边非嵌套环路，按照本文算法均可予以有效破除。

1.4 特征图分解与特征线筛选

在破开小的冗余环路后，得到的是由特征边组成的可能带有较大环路的特征图，其中还存在大量显著度低的特征线，因此还要对特征线进行筛选。

现有方法[12, 14-15, 17]通常根据特征线的长度进行筛选，剔除短小的末端特征线分支，然而仅根据长度并不足以评判特征线的重要性，短小的末端分支也可能是构成显著的主山脊线的一部分。Zou等人[18]提出了一种基于特征显著度的筛选方法，将特征线上特征点的显著度之和作为特征线的显著度，以显著度最大化原则从特征树中搜寻主干，继而将特征树分解为一条条的特征线，然后根据显著度对特征线排序，保留最显著的若干条特征线。由于在搜寻主干时追求显著度最大化，该方法可能造成主干的过度延伸，而特征线显著度的定义方式决定了由弱特征点构成的较长特征线的显著度可能超过由强特征点构成的较短特征线，这样会造成主干提取不准确，以及特征线的显著度排序可能与实际情况不符；再者，该方法并不适用于带环路的特征图的分解。

针对以上问题，提出了一种新的特征图分解与特征线筛选算法。算法中采用了以下定义：

1) 特征点显著度。使用式(1) 计算得到的特征点的显著度值。

2) 分支。一条特征线段上除去两个端点后剩下点的度均为2，且两个端点均满足度大于2或度为1，则称该条线段为分支；对于由度为2的点构成的环路，也称其为一分支，其两个端点相同。

3) 分支显著度。分支上所有特征点的显著度之和。

4) 分支平均显著度。分支上所有特征点的显著度的平均值。

5) 分支走向。连接分支两端点所构成线段的方向。

6) 最显著分支。在特征图的一个连通分量中，拥有最大显著度的特征点所在的分支，若特征图中有多个连通分量，则会有多个最显著分支。

特征图分解的目标是将特征图分解为一条条特征线，为了解决过度延伸的问题，要求分解得到的特征线内部各分支的平均显著度相近。分解算法的思想是将特征图分解为分支集合，然后利用最大显著度点往往出现在主干特征线上这一事实，从最显著分支开始，根据分支走向及分支平均显著度的相似度来连接其他相邻分支，提取主干特征线，然后从主干特征线的各分叉点开始类似地通过连接相邻分支来提取其他特征线。算法步骤具体描述如下：

1) 扫描特征图，对每一未访问过的特征点。遍历其所在的连通分量，提取该连通分量中的所有分支，在此过程中计算各分支的显著度、长度以及平均显著度，并找到该连通分量的最显著分支。

2) 从每个最显著分支开始，提取其所属特征图连通分量的主干特征线。将最显著分支设为当前分支，往两端遍历与其相邻的未访问分支，针对与当前分支走向一致的分支(两分支走向偏转角不超过90°)，利用式(2) 计算其平均显著度与最显著分支的平均显著度之间的差异，找到差异最小的相邻分支，如果此差异大于指定阈值，则终止相应方向的遍历，否则将该相邻分支纳入主干特征线并设为当前分支，继而用同样的方法继续遍历。当延最显著分支两个端点方向的遍历都终止时，主干特征线提取完毕。在提取主干特征线过程中，将遍历到的度大于2的分支端点都加入集合$\mathit{\boldsymbol{BP}}$中。

$ dif = \frac{{|Av{e_{{\rm{MB}}}} - Av{e_{{\rm{CB}}}}|}}{{Av{e_{{\rm{MB}}}}}} $ (2)

式中，$dif$为当前分支的平均显著度${Av{e_{{\rm{CB}}}}}$与最显著分支的平均显著度${Av{e_{{\rm{MB}}}}}$之间的差异。

3) 若$\mathit{\boldsymbol{BP}}$为空，则算法结束；否则，从中取出一分支点$p$。

4) 选取以$p$为端点的一未访问分支作为参考分支，沿该分支的另一端点方向遍历相邻分支，采用与步骤2) 类似的方法，根据分支走向及与参考分支之间的平均显著度差异来选取相邻分支，纳入到当前提取的特征线中，并将遍历到的度大于2的分支端点都加入到集合$\mathit{\boldsymbol{BP}}$中。在结束当前特征线的提取后，考察是否有其他以$p$为端点的未访问分支，有则将其作为参考分支继续提取特征线，没有则转第3) 步。

在将特征图分解为一系列特征线后(假定共有$n$ 条特征线)，计算各特征线的显著度，将特征线按照显著度从高到底排序，保留显著度最大的前$m$条特征线作为最终提取结果($m$可由用户事先指定，或是在排序后由用户交互式地指定)。其中特征线${{L_i}}$的显著度

$ \begin{array}{l} S\left( {{L_i}} \right) = \\ \left( {\frac{{Av{e_i} - Av{e_{{\rm{min}}}}}}{{Av{e_{{\rm{max}}}} - Av{e_{{\rm{min}}}}}}} \right){^{{k_{\rm{a}}}}} \times \left( {\frac{{Le{n_i} - Le{n_{{\rm{min}}}}}}{{Le{n_{{\rm{max}}}} - Le{n_{{\rm{min}}}}}}} \right){^{{k_{\rm{l}}}}} \times \\ \left( {\frac{{Ma{x_i} - Ma{x_{{\rm{min}}}}}}{{Ma{x_{{\rm{max}}}} - Ma{x_{{\rm{min}}}}}}} \right){^{{k_{\rm{m}}}}} \end{array} $ (3)

式中，$k$$_{\rm{a}}$、$k$$_{\rm{l}}$、$k$$_{\rm{m}}$为权重系数，${Av{e_i}}$为${{L_i}}$上所有特征点的平均显著度，${Av{e_{{\rm{min}}}}}$ =min{${Av{e_i}}$|$i$∈[1，$n$]}，${Av{e_{{\rm{max}}}}}$ = max{${Av{e_i}}$|$i$∈[1，$n$]}，${Le{n_i}}$为构成${{L_i}}$的特征点个数，也即${{L_i}}$的长度，而${Le{n_{{\rm{min}}}}}$和${Le{n_{{\rm{max}}}}}$分别为所有特征线长度中的最小和最大值，${Ma{x_i}}$为${{L_i}}$上特征点显著度的最大值，${Ma{x_{{\rm{min}}}}}$ =min{${Ma{x_i}}$ |$i$∈[1，$n$]}，${Ma{x_{{\rm{max}}}}}$ =max{${Ma{x_i}}$|$i$∈[1，$n$]}。

2 实验结果与分析

为了验证本文算法的有效性，选用了某山区30 m精度DEM数据(如图 4所示)进行测试。由于本文算法实现的是显著度可控的DEM地形特征线提取，因此选择与文献[18]基于特征显著度的地形特征线提取算法进度对比测试。



图 4 用于对比实验的某山区DEM的灰度图
Fig. 4 Grayscale images of two real DEMs for experiment
((a) 480×480 pixels DEM; (b) 512×512 pixels DEM)

在实验中，文献[18]算法为了保证得到单像素宽的特征线，须进行下采样，实验中每3行/列取1行/列数据，而本文算法采用了细线化来保证得到单像素宽的特征线，所以无需进行下采样(在实验中，针对图 4(a)地形未进行下采样，针对图 4(b)进行了下采样，采样率与文献[18]算法一致)；本文算法与文献[18]算法采用相同的特征点提取算法，提取特征点时用于鉴别平顶山脊或平底山谷与其他平坦地形特征(例如高原、平原和盆地)的宽度阈值取值为15，用于鉴别噪声的高度阈值取值为该地形图最大高度差的3 %；特征点显著度的计算方法也相同(按式(1) 计算)，权重系数$k$h、$k$d和$k$s分别取3、2、1。本文算法所独有的参数设置如下：在进行破环时保留环内面积大于105个像素的环路；特征图分解算法步骤2) 中的平均显著度差异阈值通常设为60 % 80 %之间；按式(3) 计算特征线显著度时权重系数$k$$_{\rm{a}}$、$k$$_{\rm{l}}$、$k$$_{\rm{m}}$均取1。本文算法提取山谷线的方法是将DEM中地形高度取反后，再进行山脊线提取。

实验结果对比如图 5-图 8所示，其中最显著的5条特征线用不同颜色标识出来，对应数字表示特征线的显著度，在提取的前40条最显著山脊线结果中，用蓝圈标出了部分在对比中有明显缺陷之处。



图 5 最显著山脊线提取结果对比1
Fig. 5 Comparison one of extracting the most prominent ridge lines
((a) reference[18] result (5 ridge lines); (b) ours (5 ridge lines); (c) reference[18] result (40 ridge lines); (d) ours (40 ridge lines))

文献[18]算法将特征点的显著度之和作为特征线的显著度，同时在进行特征树分解时追求主干的显著度最大化，造成了主干的过度延伸(例如图 5(a)中的绿色特征线以及图 6(a)、图 7(a)、图 8(a)中的红色特征线等等)以及主干提取不准确(例如图 5(a)中的黄色和蓝色特征线，图 6(a)中的黄色特征线，以及图 8(a)中的黄、绿、蓝特征线等等)，此外，也使得特征线的显著度排序与人眼观察有较大出入(例如图 5(a)中的紫色特征线显著程度明显不如图 5(b)中的蓝色特征线和其他4条特征线，不应出现在最显著的前5条特征线中)。



图 6 最显著山谷线提取结果对比1
Fig. 6 Comparison one of extracting the most prominent valley lines
((a) reference[18] result (5 valley lines); (b) ours (5 valley lines); (c) reference[18] result (40 valley lines); (d) ours (40 valley lines))



图 7 最显著山脊线提取结果对比2
Fig. 7 Comparison two of extracting the most prominent ridge lines
((a) reference[18] result (5 ridge lines); (b) ours (5 ridge lines); (c) reference[18] result (40 ridge lines); (d) ours (40 ridge lines))



图 8 最显著山谷线提取结果对比2
Fig. 8 Comparison two of extracting the most prominent valley lines
((a) reference[18] result (5 valley lines); (b) ours (5 valley lines); (c) reference[18] result (40 valley lines); (d) ours (40 valley lines))

而本文算法在进行特征图分解时将分支平均显著度相近的分支合并构成特征线，较好地避免了特征线的过度延伸，此外特征显著度的计算(按式(3) 计算)同时考虑了特征线的长度、线上特征点的平均显著度和最大显著度三方面因素，并将显著度表示为三者归一化后的乘积的形式，这样如果某特征线在三方面中有任一方面排名较后，总体排名都不会靠前。从最显著的前5条特征线提取结果对比中可以发现，本文算法对特征图的分解以及特征线的显著度排序更加合理，更接近人眼观察结果。

本文算法与文献[18]算法采用了同样的特征点提取和特征延伸方法，因此前40条最显著特征线提取结果看起来比较相似，但实际上有以下不同：1) 虽然部分区域的特征线提取结果整体看起来非常相近，但实际上是由不同的特征线组成的(即特征线的划分不一样)，与前5条特征线提取结果结合起来看则容易发现这一点。2) 由于本文算法在将特征图分解为特征线时考虑了构成特征线的分支之间方向和显著度的差异，避免了特征线方向在分支节点处发生急转(例如图 5(c)中上方篮圈标出的情况)，同时一些显著度较弱的分支在提取结果中被筛选掉了(例如图 7(c)和图 8(c)中上方中间蓝圈标出的分支以及一些特征线的末端分支)，如果降低特征图分解算法第2) 步中的平均显著度差异阈值，将筛选掉更多的显著度相对较弱的分支，但此阈值过小的话会影响特征线的整体性，即出现一条特征线被分解成多段的情况(图 5(d)中蓝圈标出的断裂即为此种情况，恰逢断裂处对应分支显著度较弱，被筛选掉了)。3) 由于在对图 4(a)DEM进行特征线提取时，本文算法未进行下采样，所提取的特征点与文献[18]算法对比有细微差异，本文算法支持在原始分辨率下进行特征提取，可避免下采样过程中带来的信息丢失。4) 由于特征线显著度计算方法的不同，入选前40条最显著山脊线的特征线也有差异。

为了验证本文环路检测与破环算法对小型环路的破除效果，选用图 1中的地形样本进行实验展示，实验中对地形样本进行了下采样处理，从而能够更清楚地看到破环效果。图 9给出了破环前后的中间结果，其中用黄色标出了破环前的环路及破环后的结果，由于这些环路均为小环路，需要全部破除，通过对比图 9(a)(b)可以看到，本文算法对小型环路能够实现稳定、完全的破除。



图 9 特征线提取中破环效果展示
Fig. 9 Illustration of loop breaking in feature line extraction
((a) before loop breaking; (b) after loop breaking)

由于含有大环路的DEM较难寻找，本文选取了如图 10(a)所示的含有环形山的合成地形DEM进行实验，来验证本文算法对大的环形特征线提取的支持，为了保证大型环路信息的准确性，未对该图进行下采样。图 10(b)显示了提取结果(仅保留了最显著的1条特征线)，可以看到，本文算法能够有效地保留大型环路。



图 10 环形山脊线提取
Fig. 10 Ringlike ridge line extraction
((a) synthesized DEM with a crater; (b) extraction result)

3 结论

在现有的基于特征显著度的地形特征线提取算法[18]的基础上，提出一种新的显著度可控的DEM地形特征线提取算法。第一，通过采用基于分支平均显著度相似性和分支方向一致性的特征图分解算法，以及在进行特征线显著度计算时综合考虑特征线长度、线上特征点平均显著度和最大显著度，提高了分解结果的合理性和基于显著度的特征线筛选的准确性；第二，提出一种新的环路检测和破环算法，使得能够处理各种复杂环路情况，在保留大环路的同时破除所有小环路；第三，对改进的Hilditch细线化算法[12]进一步改进，使得细线化结果与扫描顺序基本无关，且总是保留显著度较大的特征点。实验结果表明，与现有算法相比，该算法能够实现更准确的显著度控制并支持环形特征线的提取，提取结果与人眼观察结果基本一致。

然而，本文算法也存在以下不足：1) 算法参数较多，目前主要依据经验进行设置，其中特征图分解中使用的分支平均显著度差异阈值对结果影响较大，设置过小会导致特征线的过度划分，影响提取结果的整体性和连通性，而设置过大则无法剔除一些显著度较弱的分支；2) 环路检测与破环算法需要将特征点标记图放大3倍进行处理，效率不高(对于实验中所用的512×512左右的样本，在i7 2.6 GHz CPU和8 GB内存机器上执行时间在550 700 ms)，成为整个算法效率的瓶颈。因此，进一步优化参数选取，尝试实现关键参数的自动化设置，以及提高环路检测与破环算法的效率是下一步的研究目标。

参考文献

  • [1] Zhang H J, Sun J G, Lyu Y H, et al. A new simplification method for terrain model based on divergence function[J]. Chinese Journal of Computers, 2009, 32(5): 962–973. [张慧杰, 孙吉贵, 吕英华, 等. 一种新的基于发散度函数的地形模型简化方法[J]. 计算机学报, 2009, 32(5): 962–973. ] [DOI:10.3724/SP.J.1016.2009.00962]
  • [2] Zhang H J, Lü Y H, Liu S H. A terrain model simplification method based on adaptive areas division[J]. Journal of Computer Research and Development, 2010, 47(1): 53–61. [张慧杰, 吕英华, 刘淑华. 一种基于自适应区域分割的地形模型简化方法[J]. 计算机研究与发展, 2010, 47(1): 53–61. ]
  • [3] Zhou H, Sun J, Turk G, et al. Terrain synthesis from digital elevation models[J]. IEEE Transactions on Visualization and Computer Graphics, 2007, 13(4): 834–848. [DOI:10.1109/TVCG.2007.1027]
  • [4] Tasse F P, Gain J, Marais P. Enhanced texture-based terrain synthesis on graphics hardware[J]. Computer Graphics Forum, 2012, 31(6): 1959–1972. [DOI:10.1111/j.1467-8659.2012.03076.x]
  • [5] Xiong H J, Li X J. A new method to extract terrain feature lines[J]. Geomatics and Information Science of Wuhan University, 2015, 40(4): 498–502, 515. [熊汉江, 李秀娟. 一种提取山脊线和山谷线的新方法[J]. 武汉大学学报·信息科学版, 2015, 40(4): 498–502, 515. ] [DOI:10.13203/j.whugis20140004]
  • [6] Liu L, Zhang D, You J. Detecting wide lines using isotropic nonlinear filtering[J]. IEEE Transactions on Image Processing, 2007, 16(6): 1584–1595. [DOI:10.1109/TIP.2007.894288]
  • [7] Li S X, Chang H X, Zhu C F. Fast curvilinear structure extraction and delineation using density estimation[J]. Computer Vision and Image Understanding, 2009, 113(6): 763–775. [DOI:10.1016/j.cviu.2009.01.003]
  • [8] Liu Z J, Hu B G. High resolution DEM topographic feature line extraction algorithm using GPU[J]. Journal of Image and Graphics, 2012, 17(2): 249–255. [刘洲俊, 胡包钢. GPU加速的高分辨率DEM图像地形特征线提取算法[J]. 中国图象图形学报, 2012, 17(2): 249–255. ] [DOI:10.11834/jig.20120214]
  • [9] Kong Y P, Fang L, Jiang Y L, et al. A new method of extracting terrain feature lines by morphology[J]. Geomatics and Information Science of Wuhan University, 2012, 37(8): 996–999. [孔月萍, 方莉, 江永林, 等. 提取地形特征线的形态学新方法[J]. 武汉大学学报·信息科学版, 2012, 37(8): 996–999. ]
  • [10] Tribe A. Automated recognition of valley lines and drainage networks from grid digital elevation models:a review and a new method[J]. Journal of Hydrology, 1992, 139(1-4): 263–293. [DOI:10.1016/0022-1694(92)90206-B]
  • [11] Zhu Q, Zhao J, Zhong Z, et al. 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. ] [DOI:10.3321/j.issn:1001-1595.2004.01.014]
  • [12] Chen Y L, Liu D Y. A new method for automatic extraction of ridge and valley axes from DEM[J]. Journal of Image and Graphics, 2001, 6(12): 1230–1234. [陈永良, 刘大有. 一种新的山脊线和山谷线自动提取方法[J]. 中国图象图形学报, 2001, 6(12): 1230–1234. ] [DOI:10.11834/jig.2001012260]
  • [13] Huang P Z, Liu Z H. Extraction of ridge and valley from DEM based on gradient[J]. Geomatics and Information Science of Wuhan University, 2005, 30(5): 396–399. [黄培之, 刘泽慧. 基于地形梯度方向的山脊线和山谷线的提取[J]. 武汉大学学报·信息科学版, 2005, 30(5): 396–399. ] [DOI:10.3321/j.issn:1671-8860.2005.05.005]
  • [14] Chang Y C, Song G S, Hsu S K. Automatic extraction of ridge and valley axes using the profile recognition and polygon-breaking algorithm[J]. Computers & Geosciences, 1998, 24(1): 83–93. [DOI:10.1016/S0098-3004(97)00078-2]
  • [15] Chang Y C, Sinha G. A visual basic program for ridge axis picking on DEM data using the profile-recognition and polygon-breaking algorithm[J]. Computers & Geosciences, 2007, 33(2): 229–237. [DOI:10.1016/j.cageo.2006.06.007]
  • [16] Bangay S, de Bruyn D, Glass K.Minimum spanning trees for valley and ridge characterization in digital elevation maps[C]//Proceedings of the 7th International Conference on Computer Graphics, Virtual Reality, Visualisation and Interaction in Africa.Franschhoek, South Africa:ACM, 2010:73-82.[DOI:10.1145/1811158.1811171]
  • [17] Zhang H J, Liu Y X, Ma Z Q, et al. A terrain skeleton feature extraction method based on morphological encoding[J]. Journal of Computer Research and Development, 2015, 52(6): 1409–1423. [张慧杰, 刘亚鑫, 马志强, 等. 一种基于形态学编码的地形骨架特征提取方法[J]. 计算机研究与发展, 2015, 52(6): 1409–1423. ] [DOI:10.7544/issn1000-1239.2015.20131422]
  • [18] Zou K, Wo Y, Xu X.A feature significance-based method to extract terrain feature lines[J/OL]. Geomatics and Information Science of Wuhan University, (2016-11-01)[2017-02-01]. http://www.cnki.net/kcms/detail/42.1676.TN.20161101.1145.003.html. [邹昆, 沃焱, 徐翔. 利用特征显著度提取地形特征线的方法[J/OL]. 武汉大学学报·信息科学版, (2016-11-01)[2017-02-01]. http://www.cnki.net/kcms/detail/42.1676.TN.20161101.1145.003.html.[DOI:10.13203/j.whugis20150373]]