CN105469061A - 地形特征线提取方法及装置 - Google Patents

Abstract

本发明实施例提供了一种地形特征线提取方法及装置，改善了现有技术中特征线提取方法无法根据特征显著程度对提取结果进行方便筛选的问题。该地形特征线提取方法，包括：对采集到的地形的数字高程模型DEM数据进行预处理；提取出多个特征点并计算各特征点的特征显著度；沿每个特征点的特征方向或与该特征点的特征方向相邻的方向进行特征延伸；得到特征图，并根据最小生成树算法生成特征森林；对特征森林进行枝干分解，提取出所有特征线；从各特征线中选取出显著度最高的预设数量的特征线作为地形的特征线集合。使用该地形特征线提取方法及装置，可以根据特征显著度对提取结果进行方便的筛选，提高了特征性提取的抗噪能力和特征点连通性。

Description

地形特征线提取方法及装置

技术邻域

本发明涉及地形处理技术邻域，具体而言，涉及一种地形特征线提取方法及装置。

背景技术

山脊线和山谷线是非常重要的地形特征线，描述了最主要的地形骨架结构，自动提取山脊线和山谷线对于地形地貌研究有着重要意义，近年来也被用于基于样本的地形合成。

从规则格网数字高程模型(DigitalElevationModel，DEM)中提取山脊线和山谷线一直是一个研究热点，主要方法包括地表流水模拟分析法、几何形态分析法和基于图像处理技术的方法。

其中，地表流水模拟分析法依据地表水流从高处流向低处的自然规律，通过计算汇水量生成连续河网，在此基础上细化河网水系得到山谷线，而将河网公共边界作为山脊线，此类方法提取结果整体性较好，但计算量大，效率较低，且对平原地区水流方向的判定比较困难；几何形态分析法通常利用山脊和山谷点为断面极大或极小值点，采用局部窗口对断面进行分析以确定特征点，此类方法简单高效，提取精度相对较高，但结果容易受窗口大小的影响，且对噪声敏感，特征线连续性较差；基于图像处理的方法将DEM看作灰度图，利用图像处理中的技术进行线型特征检测，如利用图像处理中的灰度形态学算子提取山脊线和山谷线，取得较好效果，但在平坦山脊和山谷处易产生断裂。

上述方法均不能根据特征强弱对提取结果进行方便的筛选控制，此外基于规则格网DEM的特征线提取方法总体上抗噪性差，且特征点连接成线时较为困难。

发明内容

有鉴于此，本发明实施例的目的在于提供一种地形特征线提取方法及装置，以改善现有技术中基于规则格网DEM的特征线提取方法不能根据特征强弱对提取结果进行方便的筛选控制，总体上抗噪性差，且特征点连接成线时较为困难的问题。

为了实现上述目的，本发明实施例采用的技术方案如下：

第一方面，本发明实施例提供了一种地形特征线提取方法，包括：

对采集到的地形的数字高程模型DEM数据进行预处理；

分别沿多个方向对所述DEM数据进行扫描，提取出多个特征点并计算各所述特征点的特征显著度；

沿每个特征点的特征方向或与该特征点的特征方向相邻的方向进行特征延伸；

将每个特征点作为顶点，将相邻特征点连线构成边得到特征图，并根据最小生成树算法生成特征森林；

对所述特征森林进行枝干分解，提取出所有特征线，将各所述特征线上的所有特征点的特征显著度之和作为各所述特征线的显著度；

从各所述特征线中选取出显著度最高的预设数量的特征线作为所述地形的特征线集合。

第二方面，本发明实施例提供了一种地形特征线提取装置，包括：

预处理单元，用于对采集到的地形的数字高程模型DEM数据进行预处理；

数据扫描单元，用于分别沿多个方向对经所述预处理单元预处理后的DEM数据进行扫描，提取出多个特征点并计算各所述特征点的特征显著度；

特征延伸单元，用于沿每个特征点的特征方向或与该特征点的特征方向相邻的方向进行特征延伸；

特征森林生成单元，用于将每个特征点作为顶点，将相邻特征点连线构成边得到特征图，并根据最小生成树算法生成特征森林；

特征线提取单元，用于对所述特征森林进行枝干分解，提取出所有特征线，将各所述特征线上的所有特征点的特征显著度之和作为各所述特征线的显著度；

特征线集合选取单元，用于从各所述特征线中选取出显著度最高的预设数量的特征线作为所述地形的特征线集合。

本发明实施例中所提供的地形特征线提取方法及装置，能够根据特征显著度对提取结果进行筛选控制，提取预设数量的最显著的山谷或山脊线作为地形的特征线集合，提取的特征线与实际地形相符，有良好的抗噪能力且能较好地处理平坦地形特征，连通性较好。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本邻域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1示出了本发明实施例所提供的一种方法流程示意图；

图2示出了本发明实施例所提供的一种地形断面示意图；

图3示出了本发明实施例所提供的一种水平方向断面扫描示意图；

图4示出了本发明实施例所提供的一种垂直方向断面扫描示意图；

图5示出了本发明实施例所提供的一种右斜45度方向断面扫描示意图；

图6示出了本发明实施例所提供的一种左斜45度方向断面扫描示意图；

图7示出了本发明实施例所提供的一种全局断面行数据处理流程示意图；

图8示出了本发明实施例所提供的一种提取关键点示意图；

图9示出了本发明实施例所提供的一种提取候选特征点示意图；

图10示出了本发明实施例所提供的一种提取特征点示意图；

图11示出了本发明实施例所提供的一种特征点属性及坡面角示意图；

图12示出了本发明实施例所提供的一种特征延伸示意图；

图13示出了本发明实施例所提供的一种枝干分解示意图；

图14示出了本发明实施例所提供的一种实验区域灰度晕渲图；

图15示出了本发明实施例所提供的一种特征延伸前的山脊特征点的示意图；

图16示出了本发明实施例所提供的一种特征延伸后的山脊特征点的示意图；

图17示出了本发明实施例所提供的一种枝干裁剪结果示意图；

图18示出了本发明实施例所提供的另一种枝干裁剪结果示意图；

图19示出了本发明实施例所提供的另一种枝干裁剪结果示意图；

图20示出了本发明实施例所提供的一种处理后的山脊线示意图；

图21示出了本发明实施例所提供的一种处理后的山谷线示意图；

图22示出了本发明实施例所提供的一种装置的结构示意图；

图23示出了本发明实施例所提供的另一种装置的结构示意图。

具体实施方式

下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本邻域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

考虑到现有技术中基于规则格网DEM的特征线提取方法不能根据特征强弱对提取结果进行方便的筛选控制，总体上抗噪性差，且特征点连接成线时较为困难。本发明实施例提出了一种新的地形特征线提取方法，通过基于特征显著度的裁剪实现对提取结果特征显著程度的控制，通过全局断面扫描避免窗口参数选择，提高抗噪能力及对平坦地形特征的妥善处理能力，通过基于特征方向的特征延伸提高连通性。

如图1所示，该地形特征线提取方法，包括：

步骤S100：对采集到的地形的数字高程模型DEM数据进行预处理；

于本步骤中，优选采用高斯滤波与下采样对DEM数据进行预处理。采用下采样一方面是为了保证最终得到的是单像素宽的特征线，另一方面是为了降低后续步骤的计算量。下采样时，在原DEM中沿水平/垂直方向每n列/行取1个样本点，n值可根据DEM分辨率和提取精度要求选取。为了提高下采样后DEM的保真度，同时过滤细微噪声，优选在下采样前先进行高斯滤波，滤波器窗口大小可设为2n-1。

步骤S101：分别沿多个方向对所述DEM数据进行扫描，提取出多个特征点并计算各所述特征点的特征显著度；

本步骤是进行全局断面扫描和处理。断面极值法是一种简单高效的确定候选特征点的方法，通常将断面上的高程极大值点作为山脊线上的候选点，将极小值点作为山谷线上的候选点。现有技术中通常采用局部窗口进行断面分析，采用局部窗口分析时利用的信息有局限性，效果容易受窗口大小和噪声影响。如图2给出了一个DEM断面示意图，人眼易判断出A处为平原，B处为小平台，C、E两处为山脊，D处为山谷，D和F处的小突起为噪声，G处可能为山谷，但由于过于靠近边界且右方只有微小上升，无法确定。但如果使用局部窗口进行分析，如果窗口过小，则会丢失C处中间的山脊点，而D处和F处的小突起易被误判为山脊点，A的两侧及B处的平坦特征也可能被误判。造成这种差异的原因在于人眼是利用全局信息进行判断的，受此启发本步骤中采用基于全局断面扫描的特征点确定方法。为了尽量不遗漏特征，本步骤中优选分别沿水平(参阅图3)、垂直(参阅图4)及两个斜45度方向(参阅图5、图6)对DEM数据进行逐行扫描，并对每行数据进行处理，提取特征点并计算特征显著度。其中，两个斜45度方向的扫描优选忽略两边像素数少于3的行。而全局断面行数据处理流程如图7所示，具体描述如下。

步骤S200：分别沿水平、垂直和两个斜四十五度方向对所述DEM数据进行逐行扫描。

步骤S201：逐个访问扫描得到的每行DEM数据中的每个数据点p_i。

步骤S202，提取关键点：提取出每行DEM数据中的起始数据点和终止数据点，将所述起始数据点和终止数据点添加至关键点集合K；计算每行DEM数据中非起始数据点和终止数据点的其余数据点分别与相邻数据点的高程差值，根据所述高程差值判断各所述数据点的前后地形走势，将前后地形走势相反的数据点添加至所述关键点集合K；从位于平坦地势起点的数据点开始不断访问后续数据点直至高程差值不为0，将平坦地势起点至高程差值不为0的数据点的距离l作为平坦地势长度，判断距离l是否大于W_THRESH，如果是，则将位于平坦地势起点的所述数据点添加至所述关键点集合K，其中，W_THRESH为鉴别平坦线型特征和平坦区域特征的阈值；

于本步骤中，逐个访问每一数据点p_i，如果是首尾点则直接加入关键点集合K，否则通过计算数据点与相邻数据点的高程差值判断数据点前后的地形走势，有三种情况。

第一种，数据点前后地形走势同为上升或下降，则不做任何处理。

第二种，数据点前后地形走势刚好相反，则将此数据点加入关键点集合K。

第三种，数据点p_i位于平坦地势起点，则不断访问后续点直至高程差值不为0，从而获知平坦地势长度l及后续走势。如果l大于阈值W_THRESH，或是前后走势相反，则将数据点p_i加入关键点集合K。其中W_THRESH用于鉴别平坦线型特征(例如平顶山脊)和平坦区域特征(例如平原)。

由步骤S202得到的关键点主要是地形断面上的折点，但剔除了类似图2中B处的平坦地势点，此外将断面首尾点加入关键点集合K主要是为了方便后续分析。本步骤中针对图2提取关键点的结果如图8所示，其中，关键点后方的水平线条给出了对应关键点的长度(非平坦地势起点的长度为0)。

步骤S203，剔除噪声和伪特征，提取候选特征点：根据地形走势和所述关键点集合K中相邻关键点的高度差进行分析判断，去除所述关键点集合K中的噪声和伪特征，得到候选特征点集合C；

于本步骤中，优选对关键点集合K进行一次遍历，根据当前地形总体走势Trend以及相邻关键点的高度差进行分析判断，当前地形总体走势Trend主要分为三种：UNDEFINED(平缓)、ASCEND(上升)、DESCEND(下降)。将首个关键点直接加入候选特征点集合C，置Trend为UNDEFINED，然后逐个访问后续关键点。设前一候选特征点为c，前一关键点为k_i-1，噪声阈值为V_THRESH，对每一关键点k_i进行如下操作。

如果关键点k_i与候选特征点c高度相同，则视二者处于同一平坦地势，做以下处理：修改c的长度，如果长度大于W_THRESH，置Trend为UNDEFINED，转访问下一关键点。

如果关键点k_i的长度大于W_THRESH，置Trend为UNDEFINED，将关键点k_i加入候选特征点集合C，转访问下一关键点。

判断当前Trend的值，有三种情况：

第一种，Trend为UNDEFINED。如果关键点k_i与候选特征点c的高度差Δh的绝对值大于V_THRESH，将关键点k_i加入候选特征点集合C，并根据高度差Δh的符号修改Trend(Δh>0则置Trend为ASCEND，Δh<0则置Trend为DESCEND)。

第二种，Trend为ASCEND。如果关键点k_i的高度大于关键点k_i-1的高度，且也大于候选特征点c的高度，说明候选特征点c为伪特征点，将其替换为关键点k_i；如果关键点k_i的高度小于关键点k_i-1的高度，且候选特征点c与关键点k_i的高度差大于V_THRESH，则将关键点k_i加入候选特征点集合C，并置Trend为DESCEND。

第三种，Trend为DESCEND。如果关键点k_i的高度小于关键点k_i-1的高度，且也小于选特征点c的高度，说明选特征点c为伪特征点，将其替换为关键点k_i；如果关键点k_i的高度大于关键点k_i-1的高度，且关键点k_i与选特征点c的高度差大于V_THRESH，则将关键点k_i加入候选特征点集合C，并置Trend为ASCEND。

另，对于末尾关键点，可不用经历本步骤，直接加入候选特征点集合C。图9给出了提取的候选特征点，图8中的关键点k₅、k₈、k₉被当作噪声剔除，而关键点k₆由于和前一候选特征点c₄高度相同，合并到候选特征点c₄，延长了候选特征点c₄的长度。

步骤S204，提取特征点：将所述候选特征点集合C中长度小于W_THRESH且与两侧相邻候选特征点的高度差的绝对值大于V_THRESH的候选特征点加入特征点集合F中，其中，V_THRESH为噪声阈值；

于本步骤中，优选对候选特征点集合C进行一次遍历，将长度小于W_THRESH且与两侧相邻候选特征点的高度差绝对值均大于V_THRESH的候选特征点加入特征点集合F。不用判断两侧高度差符号相反是因为步骤S203中已保证这一点。图10给出了提取的特征点的示意图。

总体上看，地形特征点分为两类：山脊点和山谷点，本步骤提取的特征点既有山脊点也有山谷点，在一次包含水平、垂直及两个斜45度方向四个方向的逐行扫描(全局扫描)中可同时检测出来，特征线(山脊线或山谷线)是由一系列相邻的特征点构成的。

步骤S205，计算特征显著度：根据所述特征点集合F中每个特征点的高度、每个特征点与相邻候选特征点的高度差，以及每个特征点与相邻候选特征点的距离计算得到各所述特征点的特征显著度。

对每一特征点f_i已记录有以下属性供计算特征显著度：高度h_i，与相邻候选特征点的高度差Δh_li和Δh_ri，以及与相邻候选特征点的距离d_li和d_ri，其中，距离用像素单位即可，对于斜45度方向的断面，须乘以缩放因子。例如，图10中f₁的属性如图11所示。对于山脊点的显著度计算，考虑以下因素：

归一化高度：高度越高，山脊特征越显著。

平均归一化高度落差：两侧高度落差越大，山脊特征越显著。这里高度差取的是与相邻候选特征点的高度差，而不是与相邻特征点的高度差，因为平原上山脊的高度落差应计算山脊点与平原的海拔差值，例如图10中f₁的左侧落差Δh_l1是其与图9中c₂的高度差值，如图11所示。

平均归一化坡度：指沿断面扫描方向的特征点两侧坡度，坡度越大，山脊特征越显著。为了简化计算，采用坡面角的正弦值作为归一化坡度，例如图11中f₁两侧的坡面角分别为α_l1和α_r1。

由此得到山脊点的特征显著度的计算公式为：

$\begin{array}{c}Sig\left(f_i\right)=k_h\frac{h_i}{h_{\max }}+k_d\frac{\left|{\mathrm{\&Delta;h}}_{li}\right|+\left|{\mathrm{\&Delta;h}}_{ri}\right|}{2\left|{\mathrm{\&Delta;h}}_{\max }\right|}\\ +\frac{k_s}{2}(\frac{\left|{\mathrm{\&Delta;h}}_{li}\right|}{\sqrt{{{\mathrm{\&Delta;h}}_{li}}^2+{d_{li}}^2}}+\frac{\left|{\mathrm{\&Delta;h}}_{ri}\right|}{\sqrt{{{\mathrm{\&Delta;h}}_{ri}}^2+{d_{ri}}^2}})\end{array}.$

其中，h_i为特征点f_i的高度，Δh_li和Δh_ri为特征点f_i与相邻候选特征点的高度差，d_li和d_ri为特征点f_i与相邻候选特征点的距离，h_max为所述地形的DEM数据中的最大高度，h_min为所述地形的DEM数据中的最小高度，Δh_max为所述地形的DEM数据中的最大高度落差，取值为h_max-h_min，k_h为高度权重系数，k_d为高度落差权重系数，k_s为坡度权重系数。优选k_h为归一化高度权重系数，k_d为平均归一化高度落差权重系数，k_s为平均归一化坡度权重系数。

山谷点的显著度计算考虑的因素类似，区别仅在于归一化高度，此时高度越低，特征越显著，因此计算公式为：

$\begin{array}{c}Sig\left(f_i\right)=k_h\frac{h_{\max }-h_i}{h_{\max }}+k_d\frac{\left|{\mathrm{\&Delta;h}}_{li}\right|+\left|{\mathrm{\&Delta;h}}_{ri}\right|}{2\left|{\mathrm{\&Delta;h}}_{\max }\right|}\\ +\frac{k_s}{2}(\frac{\left|{\mathrm{\&Delta;h}}_{li}\right|}{\sqrt{{{\mathrm{\&Delta;h}}_{li}}^2+{d_{li}}^2}}+\frac{\left|{\mathrm{\&Delta;h}}_{ri}\right|}{\sqrt{{{\mathrm{\&Delta;h}}_{ri}}^2+{d_{ri}}^2}})\end{array}.$

公式中的i是特征点的编号，可根据Δh_li或Δh_ri的符号判定特征点f_i是山脊点还是山谷点。对于长度l大于0的平坦地势特征点f，优选将其后的l个点均作为特征点，并将它们的特征显著度设为与f相同，但将处于中间位置的特征点显著度加上一个微小增量，从而引导最终生成的单像素特征线的位置处于平顶山脊或平底山谷的中央。由于需要进行四次全局扫描，对于每个DEM样本点，优选保存其特征显著度最大值，并记录对应的扫描方向。

步骤S102：沿每个特征点的特征方向或与该特征点的特征方向相邻的方向进行特征延伸；

仅采用上述步骤可能会有遗漏，无法提取出所有的特征点，造成特征线的断裂，即一条特征线可能由于中间某些特征点没有提取出来从而被分隔成了几条短的特征线，特征延伸的目的就是根据已经检测出的特征点，尝试沿特征方向延伸去修补这些断裂。如果一条显著的特征线因为断裂被分成了几段，那后续操作是将其作为几条不同的特征线处理的，其特征显著度也会大大降低，从而可能会在基于特征显著度的筛选中被筛掉，即使不被筛掉，提取的结果也是与实际地形不符的，因而需进行特征延伸。

虽然沿四个方向进行了全局断面扫描，但仍可能存在局部特征断裂，尤其是在显著特征线与其分支特征线的连接处附近，例如主脉与支脉的连接处附近，这部分区域地势起伏相对平缓，特征点容易被当作噪声过滤掉。对于每一个特征点尝试沿该点的特征方向进行延伸，特征方向与断面扫描方向垂直。如图12所示，对于特征点p₀(假设为山脊点)尝试沿左上和右下两个方向延伸，以沿左上方向延伸为例，需进行以下操作。

考察p₅、p₆和p₇三个特征点，如果其中有一个点已为特征点，则无需沿此方向延伸，否则继续；

依次考查p₆、p₅和p₇，在特征点p₀的8邻域内判断其是否大于等于相邻点高度，即为局部最大值，如果是则将p_i设为特征点，并将其显著度设为0，特征方向设为p₀p_i连线方向；

将新增特征点p_i设为待延伸点，并根据其特征方向尝试进行进一步延伸。

对于山谷点的延伸，将上述步骤中的局部最大值判断改为局部最小值判断即可。

由图12可知，特征点p₀的特征方向所对应的延伸方向有多个。例如:若特征方向是水平方向，延伸方向有向左延伸和向右延伸两个方向，在向左延伸时，要考虑左方、左上和左下三个方向的相邻点，在向右延伸时，要考虑右方、右上和右下三个方向的相邻点，在向左延伸时，如果左边、左上和左下有一个点已经是特征点，则无需向左延伸了，但还要尝试向右延伸，向右延伸的判断方式与向左延伸类似，且特征点p₀的特征方向所对应的多个延伸方向的判断方式均类似，因而在此不作重复说明。

为了更好地保证连通性，优选延伸操作在特征点的筛选之前进行，而由于延伸得到的特征点的显著度为0，在后面基于显著度的裁剪中很容易将“过度延伸”的特征线，即主要由延伸特征点构成的特征线裁剪掉。

步骤S103：将每个特征点作为顶点，将相邻特征点连线构成边得到特征图，并根据最小生成树算法生成特征森林；

在得到特征点后，构造特征图，将特征点作为图的顶点，将相邻特征点连线构成边，边的权值赋为两个特征点的显著度之和，如出现两条边交叉的情况，则剔除权值低的边。

接着需要打破特征线环路，生成特征森林，本实施例中优选利用Kruskal最小生成树算法生成特征森林。

步骤S104：对所述特征森林进行枝干分解，提取出所有特征线，将各所述特征线上的所有特征点的特征显著度之和作为各所述特征线的显著度；

枝干分解即从特征森林中提取特征线，图13给出了枝干分解的示意图，其中特征线越粗表示其显著度越高，分解结果是p₀p₁、p₂p₃、p₄p₇、p₆p₅四条特征线。分解方法是遍历所有叶子节点(度数为1的特征点)，对每一叶子节点p_i做如下处理：

如p_i属于一棵未访问过的特征树，通过两次深度优先遍历找到该树的最长路径，即显著度最大路径，作为该特征树的主干，对应主脉或干流。

如p_i不在主干上，则从其开始进行路径搜索，直至遇到分支节点p_j(度数大于2)，提取特征线p_ip_j，并对p_j进行考察。如果p_j不在主干上，且与p_j相连的所有叶子分支(分支一端为叶子节点)均已提取，则选取以p_j为端点的显著度最大的叶子分支进行延伸，直至遇到下一分支节点p_k，并对p_k进行与p_j相同的考察和处理，可以得出，上述过程是重复迭代进行的。

对于分解得到的所有特征线，将特征线上所有特征点的显著度之和作为特征线的显著度。

步骤S105，枝干裁剪：从各所述特征线中选取出显著度最高的预设数量的特征线作为所述地形的特征线集合。

特征森林中通常存在较多干扰特征线，现有技术中大都通过剔除长度较短的特征线段来消除干扰，然而长度并不能客观地反映特征的强弱，显著特征线的两端也可能存在较短的显著特征段，因此本发明实施例中采用基于特征显著度的特征线裁剪剔除干扰特征线。优选根据显著度对特征线进行降序排序，然后可根据用户指定的阈值进行特征线的裁剪，阈值可以根据要提取的最显著特征线的条数确定，或是根据最大显著度的百分比确定。

将特征线上特征点的坐标乘以下采样使用的n值即可得到原图中的特征位置。为了提高提取特征线的质量，在上述步骤基础上，优选还进行以下后续处理。

由于下采样丢失了部分信息，特征点可能并非断面上的极大值点，因此对其进行位置修正：沿与特征方向垂直的断面方向，以(n-1)/2为半径搜索高度最大值点代替当前特征点。

部分特征线的末尾存在不自然的短钩特征，可通过计算曲率和到尾端的距离进行判断并移除。

为了使特征线更加平滑和自然，对特征线进行平滑处理，即用相邻特征点的坐标均值作为当前特征点的坐标。

本发明实施例中所提供的地形特征线提取方法，能够根据特征显著度对提取结果进行筛选控制，提取预设数量的最显著的山谷或山脊线作为地形的特征线集合，提取的特征线与实际地形相符，有良好的抗噪能力且能较好地处理平坦地形特征，连通性较好。

为了更好地体现本发明实施例的实现效果，进行了以下实验。

选用某山区30米DEM数据验证本方法的有效性，如图14所示。其含有丰富的山脊、山谷特征以及噪声，在左方主要山脊交汇处存在较宽的平顶山脊特征。实验中参数选择如下：下采样的n值取3，但在生成图15和图16时为了使结果显示更加清晰，未进行下采样；W_THRESH取值为15，V_THRESH取地形最大高度差的5％；在计算特征显著度时，权重系数k_h、k_d和k_s分别取3、2、1。

图15在图14基础上叠加显示了提取的山脊特征点，可以看到，特征线存在断裂情况，断裂位置一般位于主山脊与分支山脊的连接处；此外不少山脊线并非单像素的，这是因为在部分显著特征点的两侧，存在不同特征方向的相邻特征，这是沿多个方向进行断面扫描的结果，两侧特征点的特征显著度相对较弱，在生成特征森林时会形成主要山脊旁的短枝，参阅图17。在后续的枝干裁剪中很容易被剔除；图16在图14基础上叠加显示了特征延伸后的特征点，可以看到特征延伸大幅增强了特征线的连通性，有部分特征线主要由延伸的特征点构成，由于这类特征点显著度为0，这些过度延伸的特征线在枝干裁剪中同样易被剔除。

图17～19给出了针对山脊线提取的三种枝干裁剪的结果，分别根据特征显著度取前3、前100和全部特征线(即无裁剪)。由于特征线的显著度是由特征边的显著度累加而来，实际上考虑了高度、高度落差、坡度及长度四方面因素，从图中可以看到裁剪结果与人眼观测到的地形特征显著程度基本一致。然而，由于在进行枝干分解时追求主干显著度最大化，可能将主干两端的细小分支也加到主干上，从而使这些细小分支在枝干裁剪中无法裁剪掉，例如图19中左下角特征线的长钩状特征。可通过在进行枝干分解时考虑特征线段间的特征相似性来解决此问题。

图20和图21分别给出了提取的最终山脊线和山谷线，在进行枝干裁剪时分别取前100和前50条特征线。可看到提取的特征线与实际地形是相符的，噪声和平顶特征得到了较好的处理。图21中右上和右方有两处特征线穿过了山脊，是因为这两处存在鞍点。

传统的山脊线和山谷线提取方法不能根据特征强弱对提取结果的精细程度进行方便的控制，为此本发明实施例提出了一种新的规则格网DEM地形特征线自动提取方法。该方法利用高度、高度落差和坡度信息计算特征点的显著度，通过枝干分解以及基于特征显著度的枝干裁剪实现对提取结果的控制；此外通过全局断面扫描增强特征点识别的准确性与抗噪性，以及更好地处理平坦特征。实验结果证实了该方法的有效性，在可对提取结果的特征等级提供方便有效控制的同时，提取的特征线与人工识别结果有较强一致性。

实施例2

如图22所示，本发明实施例提供了一种地形特征线提取装置，包括：预处理单元300，用于对采集到的地形的数字高程模型DEM数据进行预处理；数据扫描单元301，用于分别沿多个方向对经所述预处理单元300预处理后的DEM数据进行扫描，提取出多个特征点并计算各所述特征点的特征显著度；特征延伸单元302，用于沿每个特征点的特征方向或与该特征点的特征方向相邻的方向进行特征延伸；特征森林生成单元303，用于将每个特征点作为顶点，将相邻特征点连线构成边得到特征图，并根据最小生成树算法生成特征森林；特征线提取单元304，用于对所述特征森林进行枝干分解，提取出所有特征线，将各所述特征线上的所有特征点的特征显著度之和作为各所述特征线的显著度；特征线集合选取单元305，用于从各所述特征线中选取出显著度最高的预设数量的特征线作为所述地形的特征线集合。

如图23所示，所述数据扫描单元301，包括：逐行扫描子单元400，用于分别沿水平、垂直和两个斜四十五度方向对所述DEM数据进行逐行扫描；数据点扫描子单元401，用于逐个访问扫描得到的每行DEM数据中的每个数据点p_i；关键点集合提取子单元402，用于提取出每行DEM数据中的起始数据点和终止数据点，将所述起始数据点和终止数据点添加至关键点集合K；计算每行DEM数据中非起始数据点和终止数据点的其余数据点分别与相邻数据点的高程差值，根据所述高程差值判断各所述数据点的前后地形走势，将前后地形走势相反的数据点添加至所述关键点集合K；从位于平坦地势起点的数据点开始不断访问后续数据点直至高程差值不为0，将平坦地势起点至高程差值不为0的数据点的距离l作为平坦地势长度，判断距离l是否大于W_THRESH，如果是，则将位于平坦地势起点的所述数据点添加至所述关键点集合K，其中，W_THRESH为鉴别平坦线型特征和平坦区域特征的阈值；候选特征点集合选取子单元403，用于根据地形走势和所述关键点集合K中相邻关键点的高度差进行分析判断，去除所述关键点集合K中的噪声和伪特征，得到候选特征点集合C；特征点集合确认子单元404，用于将所述候选特征点集合C中长度小于W_THRESH且与两侧相邻候选特征点的高度差的绝对值大于V_THRESH的候选特征点加入特征点集合F中，其中，V_THRESH为噪声阈值；特征显著度计算子单元405，用于根据所述特征点集合F中每个特征点的高度h_i、每个特征点与相邻候选特征点的高度差h_li、h_ri，以及每个特征点与相邻候选特征点的距离d_li、d_ri计算得到各所述特征点的特征显著度。

其中，特征显著度计算子单元405，具体用于：

根据公式 $\begin{array}{c}Sig\left(f_i\right)=k_h\frac{h_i}{h_{\max }}+k_d\frac{\left|{\mathrm{\&Delta;h}}_{li}\right|+\left|{\mathrm{\&Delta;h}}_{ri}\right|}{2\left|{\mathrm{\&Delta;h}}_{\max }\right|}\\ +\frac{k_s}{2}(\frac{\left|{\mathrm{\&Delta;h}}_{li}\right|}{\sqrt{{{\mathrm{\&Delta;h}}_{li}}^2+{d_{li}}^2}}+\frac{\left|{\mathrm{\&Delta;h}}_{ri}\right|}{\sqrt{{{\mathrm{\&Delta;h}}_{ri}}^2+{d_{ri}}^2}})\end{array}$ 计算得到山脊特征点的特征显著度；

根据公式 $\begin{array}{c}Sig\left(f_i\right)=k_h\frac{h_{\max }-h_i}{h_{\max }}+k_d\frac{\left|{\mathrm{\&Delta;h}}_{li}\right|+\left|{\mathrm{\&Delta;h}}_{ri}\right|}{2\left|{\mathrm{\&Delta;h}}_{\max }\right|}\\ +\frac{k_s}{2}(\frac{\left|{\mathrm{\&Delta;h}}_{li}\right|}{\sqrt{{{\mathrm{\&Delta;h}}_{li}}^2+{d_{li}}^2}}+\frac{\left|{\mathrm{\&Delta;h}}_{ri}\right|}{\sqrt{{{\mathrm{\&Delta;h}}_{ri}}^2+{d_{ri}}^2}})\end{array}$ 计算得到山谷特征点的特征显著度；

其中，h_i为特征点f_i的高度，Δh_li和Δh_ri为特征点f_i与相邻候选特征点的高度差，d_li和d_ri为特征点f_i与相邻候选特征点的距离，h_max为所述地形的DEM数据中的最大高度，h_min为所述地形的DEM数据中的最小高度，Δh_max为所述地形的DEM数据中的最大高度落差，取值为h_max-h_min，k_h为高度权重系数，k_d为高度落差权重系数，k_s为坡度权重系数。

特征延伸单元302，具体用于：在山脊特征点p₀的邻域内，判断所述特征点p₀的特征方向及与所述特征点p₀的特征方向相邻的方向上与所述特征点p₀相邻的三个点中是否存在已为特征点的点，如果其中有一个点已为特征点，则无需沿所述特征点p₀的所述特征方向延伸，否则继续延伸；在所述特征点p₀的邻域内依次判断与所述特征点p₀相邻的三个点中是否存在某个点的高度大于等于相邻点高度，如果是，则将该点设为山脊特征点，并将该点的特征显著度设为0，特征方向设为该点与所述特征点p₀连线方向，并将该点设为待延伸点，根据其特征方向进行进一步延伸；

在山谷特征点p₁的邻域内，判断所述特征点p₁的特征方向及与所述特征点p₁的特征方向相邻的方向上与所述特征点p₁相邻的三个点中是否存在已为特征点的点，如果其中有一个点已为特征点，则无需沿所述特征点p₁的所述特征方向延伸，否则继续延伸；在所述特征点p₁的邻域内依次判断与所述特征点p₁相邻的三个点中是否存在某个点的高度小于相邻点高度，如果是，则将该点设为山谷特征点，并将该点的特征显著度设为0，特征方向设为该点与所述特征点p₁连线方向，并将该点设为待延伸点，根据其特征方向进行进一步延伸。

本发明实施例所提供的装置，其实现原理及产生的技术效果和前述实施例1相同，为简要描述，装置实施例部分未提及之处，可参考前述实施例1中相应内容。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术邻域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (10)

1.一种地形特征线提取方法，其特征在于，包括：

对采集到的地形的数字高程模型DEM数据进行预处理；

分别沿多个方向对所述DEM数据进行扫描，提取出多个特征点并计算各所述特征点的特征显著度；

沿每个特征点的特征方向或与该特征点的特征方向相邻的方向进行特征延伸；

将每个特征点作为顶点，将相邻特征点连线构成边得到特征图，并根据最小生成树算法生成特征森林；

对所述特征森林进行枝干分解，提取出所有特征线，将各所述特征线上的所有特征点的特征显著度之和作为各所述特征线的显著度；

从各所述特征线中选取出显著度最高的预设数量的特征线作为所述地形的特征线集合。

2.根据权利要求1所述的方法，其特征在于，所述分别沿多个方向对所述DEM数据进行扫描，提取出多个特征点并计算各所述特征点的特征显著度，包括：

分别沿水平、垂直和两个斜四十五度方向对所述DEM数据进行逐行扫描；

逐个访问扫描得到的每行DEM数据中的每个数据点p_i；

提取出每行DEM数据中的起始数据点和终止数据点，将所述起始数据点和终止数据点添加至关键点集合K；计算每行DEM数据中非起始数据点和终止数据点的其余数据点分别与相邻数据点的高程差值，根据所述高程差值判断各所述数据点的前后地形走势，将前后地形走势相反的数据点添加至所述关键点集合K；从位于平坦地势起点的数据点开始不断访问后续数据点直至高程差值不为0，将平坦地势起点至高程差值不为0的数据点的距离l作为平坦地势长度，判断距离l是否大于W_THRESH，如果是，则将位于平坦地势起点的所述数据点添加至所述关键点集合K，其中，W_THRESH为鉴别平坦线型特征和平坦区域特征的阈值；

根据地形走势和所述关键点集合K中相邻关键点的高度差进行分析判断，去除所述关键点集合K中的噪声和伪特征，得到候选特征点集合C；

将所述候选特征点集合C中长度小于W_THRESH且与两侧相邻候选特征点的高度差的绝对值大于V_THRESH的候选特征点加入特征点集合F中，其中，V_THRESH为噪声阈值；

根据所述特征点集合F中每个特征点的高度、每个特征点与相邻候选特征点的高度差，以及每个特征点与相邻候选特征点的距离计算得到各所述特征点的特征显著度。

3.根据权利要求2所述的方法，其特征在于，所述根据所述特征点集合F中每个特征点的高度、每个特征点与相邻候选特征点的高度差，以及每个特征点与相邻候选特征点的距离计算得到各所述特征点的特征显著度，包括：

根据公式 $\begin{array}{c}Sig\left(f_i\right)=k_h\frac{h_i}{h_{\max }}+k_d\frac{\left|{\mathrm{\&Delta;h}}_{li}\mathrm{|+|}{\mathrm{\&Delta;h}}_{ri}\right|}{2\left|{\mathrm{\&Delta;h}}_{\max }\right|}\\ +\frac{k_s}{2}\left(\frac{|{\mathrm{\&Delta;h}}_{li}|}{\sqrt{{{\mathrm{\&Delta;h}}_{li}}^2+{d_{li}}^2}}+\frac{|{\mathrm{\&Delta;h}}_{ri}|}{\sqrt{{{\mathrm{\&Delta;h}}_{ri}}^2+{d_{ri}}^2}}\right)\end{array}$ 计算得到山脊特征点的特征显著度；

根据公式 $\begin{array}{c}Sig\left(f_i\right)=k_h\frac{h_{\max }-h_i}{h_{\max }}+k_d\frac{\left|{\mathrm{\&Delta;h}}_{li}\mathrm{|+|}{\mathrm{\&Delta;h}}_{ri}\right|}{2\left|{\mathrm{\&Delta;h}}_{\max }\right|}\\ +\frac{k_s}{2}\left(\frac{|{\mathrm{\&Delta;h}}_{li}|}{\sqrt{{{\mathrm{\&Delta;h}}_{li}}^2+{d_{li}}^2}}+\frac{|{\mathrm{\&Delta;h}}_{ri}|}{\sqrt{{{\mathrm{\&Delta;h}}_{ri}}^2+{d_{ri}}^2}}\right)\end{array}$ 计算得到山谷特征点的特征显著度；

其中，h_i为特征点f_i的高度，Δh_li和Δh_ri为特征点f_i与相邻候选特征点的高度差，d_li和d_ri为特征点f_i与相邻候选特征点的距离，h_max为所述地形的DEM数据中的最大高度，h_min为所述地形的DEM数据中的最小高度，Δh_max为所述地形的DEM数据中的最大高度落差，取值为h_max-h_min，k_h为高度权重系数，k_d为高度落差权重系数，k_s为坡度权重系数。

4.根据权利要求1所述的方法，其特征在于，所述沿每个特征点的特征方向或与该特征点的特征方向相邻的方向进行特征延伸，包括：

在山脊特征点p₀的邻域内，判断所述特征点p₀的特征方向及与所述特征点p₀的特征方向相邻的方向上与所述特征点p₀相邻的三个点中是否存在已为特征点的点，如果其中有一个点已为特征点，则无需沿所述特征点p₀的所述特征方向延伸，否则继续延伸；在所述特征点p₀的邻域内依次判断与所述特征点p₀相邻的三个点中是否存在某个点的高度大于等于相邻点高度，如果是，则将该点设为山脊特征点，并将该点的特征显著度设为0，特征方向设为该点与所述特征点p₀连线方向，并将该点设为待延伸点，根据其特征方向进行进一步延伸；

在山谷特征点p₁的邻域内，判断所述特征点p₁的特征方向及与所述特征点p₁的特征方向相邻的方向上与所述特征点p₁相邻的三个点中是否存在已为特征点的点，如果其中有一个点已为特征点，则无需沿所述特征点p₁的所述特征方向延伸，否则继续延伸；在所述特征点p₁的邻域内依次判断与所述特征点p₁相邻的三个点中是否存在某个点的高度小于相邻点高度，如果是，则将该点设为山谷特征点，并将该点的特征显著度设为0，特征方向设为该点与所述特征点p₁连线方向，并将该点设为待延伸点，根据其特征方向进行进一步延伸。

5.根据权利要求1所述的方法，其特征在于，所述对所述特征森林进行枝干分解，包括：

遍历所述特征森林中所有度数为1的特征点p_i；

分别判断每个特征点p_i是否属于一个未访问过的特征树，通过深度优先遍历找到所述特征树中显著度最大的路径作为最长路径，将所述最长路径作为所述特征树的主干；

判断各所述特征点p_i是否在所述主干上，如果不在所述主干上，则从所述特征点p_i开始进行路径搜索，直至遇到度数大于2的分支节点p_j，提取特征线p_ip_j，如果所述分支节点p_j不在所述主干上，且与所述分支节点p_j相连的分支一端为叶子节点的所有分支另一端均为已提取特征线上的特征点，则选取以所述分支节点p_j为端点的分支一端为叶子节点的显著度最大的分支进行延伸，直至遇到下一分支节点p_k。

6.根据权利要求1～5任意一项所述的方法，其特征在于，所述方法还包括：移除所述特征线集合中的各所述特征线末尾的勾状特征，并对各所述特征线进行平滑处理。

7.一种地形特征线提取装置，其特征在于，包括：

预处理单元，用于对采集到的地形的数字高程模型DEM数据进行预处理；

数据扫描单元，用于分别沿多个方向对经所述预处理单元预处理后的DEM数据进行扫描，提取出多个特征点并计算各所述特征点的特征显著度；

特征延伸单元，用于沿每个特征点的特征方向或与该特征点的特征方向相邻的方向进行特征延伸；

特征森林生成单元，用于将每个特征点作为顶点，将相邻特征点连线构成边得到特征图，并根据最小生成树算法生成特征森林；

特征线提取单元，用于对所述特征森林进行枝干分解，提取出所有特征线，将各所述特征线上的所有特征点的特征显著度之和作为各所述特征线的显著度；

特征线集合选取单元，用于从各所述特征线中选取出显著度最高的预设数量的特征线作为所述地形的特征线集合。

8.根据权利要求7所述的装置，其特征在于，所述数据扫描单元，包括：

逐行扫描子单元，用于分别沿水平、垂直和两个斜四十五度方向对所述DEM数据进行逐行扫描；

数据点扫描子单元，用于逐个访问扫描得到的每行DEM数据中的每个数据点p_i；

关键点集合提取子单元，用于提取出每行DEM数据中的起始数据点和终止数据点，将所述起始数据点和终止数据点添加至关键点集合K；计算每行DEM数据中非起始数据点和终止数据点的其余数据点分别与相邻数据点的高程差值，根据所述高程差值判断各所述数据点的前后地形走势，将前后地形走势相反的数据点添加至所述关键点集合K；从位于平坦地势起点的数据点开始不断访问后续数据点直至高程差值不为0，将平坦地势起点至高程差值不为0的数据点的距离l作为平坦地势长度，判断距离l是否大于W_THRESH，如果是，则将位于平坦地势起点的所述数据点添加至所述关键点集合K，其中，W_THRESH为鉴别平坦线型特征和平坦区域特征的阈值；

候选特征点集合选取子单元，用于根据地形走势和所述关键点集合K中相邻关键点的高度差进行分析判断，去除所述关键点集合K中的噪声和伪特征，得到候选特征点集合C；

特征点集合确认子单元，用于将所述候选特征点集合C中长度小于W_THRESH且与两侧相邻候选特征点的高度差的绝对值大于V_THRESH的候选特征点加入特征点集合F中，其中，V_THRESH为噪声阈值；

特征显著度计算子单元，用于根据所述特征点集合F中每个特征点的高度、每个特征点与相邻候选特征点的高度差，以及每个特征点与相邻候选特征点的距离计算得到各所述特征点的特征显著度。

9.根据权利要求8所述的装置，其特征在于，所述特征显著度计算子单元，具体用于：

根据公式 $\begin{array}{c}Sig\left(f_i\right)=k_h\frac{h_i}{h_{\max }}+k_d\frac{\left|{\mathrm{\&Delta;h}}_{li}\mathrm{|+|}{\mathrm{\&Delta;h}}_{ri}\right|}{2\left|{\mathrm{\&Delta;h}}_{\max }\right|}\\ +\frac{k_s}{2}\left(\frac{|{\mathrm{\&Delta;h}}_{li}|}{\sqrt{{{\mathrm{\&Delta;h}}_{li}}^2+{d_{li}}^2}}+\frac{|{\mathrm{\&Delta;h}}_{ri}|}{\sqrt{{{\mathrm{\&Delta;h}}_{ri}}^2+{d_{ri}}^2}}\right)\end{array}$ 计算得到山脊特征点的特征显著度；

根据公式 $\begin{array}{c}Sig\left(f_i\right)=k_h\frac{h_{\max }-h_i}{h_{\max }}+k_d\frac{\left|{\mathrm{\&Delta;h}}_{li}\mathrm{|+|}{\mathrm{\&Delta;h}}_{ri}\right|}{2\left|{\mathrm{\&Delta;h}}_{\max }\right|}\\ +\frac{k_s}{2}\left(\frac{|{\mathrm{\&Delta;h}}_{li}|}{\sqrt{{{\mathrm{\&Delta;h}}_{li}}^2+{d_{li}}^2}}+\frac{|{\mathrm{\&Delta;h}}_{ri}|}{\sqrt{{{\mathrm{\&Delta;h}}_{ri}}^2+{d_{ri}}^2}}\right)\end{array}$ 计算得到山谷特征点的特征显著度；

其中，h_i为特征点f_i的高度，Δh_li和Δh_ri为特征点f_i与相邻候选特征点的高度差，d_li和d_ri为特征点f_i与相邻候选特征点的距离，h_max为所述地形的DEM数据中的最大高度，h_min为所述地形的DEM数据中的最小高度，Δh_max为所述地形的DEM数据中的最大高度落差，取值为h_max-h_min，k_h为高度权重系数，k_d为高度落差权重系数，k_s为坡度权重系数。

10.根据权利要求9所述的装置，其特征在于，所述特征延伸单元，具体用于：

在山脊特征点p₀的邻域内，判断所述特征点p₀的特征方向及与所述特征点p₀的特征方向相邻的方向上与所述特征点p₀相邻的三个点中是否存在已为特征点的点，如果其中有一个点已为特征点，则无需沿所述特征点p₀的所述特征方向延伸，否则继续延伸；在所述特征点p₀的邻域内依次判断与所述特征点p₀相邻的三个点中是否存在某个点的高度大于等于相邻点高度，如果是，则将该点设为山脊特征点，并将该点的特征显著度设为0，特征方向设为该点与所述特征点p₀连线方向，并将该点设为待延伸点，根据其特征方向进行进一步延伸；

在山谷特征点p₁的邻域内，判断所述特征点p₁的特征方向及与所述特征点p₁的特征方向相邻的方向上与所述特征点p₁相邻的三个点中是否存在已为特征点的点，如果其中有一个点已为特征点，则无需沿所述特征点p₁的所述特征方向延伸，否则继续延伸；在所述特征点p₁的邻域内依次判断与所述特征点p₁相邻的三个点中是否存在某个点的高度小于相邻点高度，如果是，则将该点设为山谷特征点，并将该点的特征显著度设为0，特征方向设为该点与所述特征点p₁连线方向，并将该点设为待延伸点，根据其特征方向进行进一步延伸。