CN105550691B - 基于尺度空间的自适应山谷山脊线提取方法及系统 - Google Patents

Abstract

本发明提供一种基于尺度空间的自适应山谷山脊线提取方法及系统，包括首先对原始DEM过采样作为处理的初始DEM，根据初始DEM建立第0组尺度空间；根据金字塔的总层数N，先挑选第0组尺度中空间尺度最大的DEM作为金字塔最底层，然后将第0层DEM进行降采样，并相应建立第1组尺度空间，再从中挑选尺度最大的DEM作为DEM金字塔的第1层，依次类推；从DEM金字塔的顶层开始，进行自适应多角度地形断面高程极值法提取，并经过后处理得到山脊线和山谷线，然后逐层精化提取结果。本发明与现有方法相比，既能兼顾山脊(谷)提取的整体趋势及细节变化，保证提取的精度，又能快速地得到提取结果，保证提取的效率。

Description

基于尺度空间的自适应山谷山脊线提取方法及系统

技术领域

本发明属于数字地形分析技术领域，涉及一种基于尺度空间的自适应山谷山脊线提取方法。

背景技术

数字高程模型(Digital Elevation Model，简称DEM)是对地球表面的一种离散化的数学表示，也是遥感和地理信息系统中进行三维空间数据处理和数字地形分析的核心数据之一。而山脊线和山谷线是可以从数字高程模型中提取的一种表示山地地形变化的重要分界线，在地形表示、数字地形分析、测绘及工程设计中都有着重要的应用。因此从DEM中自动高效地提取山脊线和山谷线一直都是数字地形分析领域的一个重要课题。

现有的提取方法从原理上可以分为整体算法和局部算法两类。整体算法主要有各种基于地形流水模拟的方法。整体算法具有较强的抗噪音的能力，但对于起伏较小的山脉，提取效果较差，并且计算量大，算法计算量随规则格网的数量的增加成平方关系增长。这些缺陷使得该方法用于山脊线和山谷线的提取有着许多不便之处。局部算法有地形断面高程极值法、曲面拟合的方法等，局部算法的主要特点是计算量小，速度快，但其不能够顾及地形的整体变化规律，使得提取结果也会有所遗漏，并且容易受到噪声影响，会给后续山脊线和山谷线的正确判定带来不便。这两类方法无法兼顾整体和细节，精度与效率，因此不能满足自动提取山脊线和山谷线的效率和精度要求。

发明内容

针对现有山脊(谷)线提取方法的缺点，本发明的目的是提供一种快速有效的从DEM中提取山脊(谷)的方法——一种基于尺度空间的自适应山谷山脊线提取方法，既能兼顾山脊(谷)提取的整体趋势及细节变化，保证提取的精度，又能快速地得到提取结果，保证提取的效率。

为实现上述目的，本发明的技术方案提供一种基于尺度空间的自适应山谷山脊线提取方法，包括以下步骤，

步骤1，尺度空间的生成，包括首先对原始DEM过采样作为处理的初始DEM，对初始DEM进行从低到高地不同程度的高斯平滑得到一组DEM的尺度空间，称为第0组尺度空间；

步骤2，DEM金字塔的建立，包括根据金字塔的总层数N，先挑选第0组尺度中空间尺度最大的DEM作为DEM金字塔的最底层，记为第0层，然后将第0层DEM进行降采样，并相应建立第1组尺度空间，再从中挑选尺度最大的DEM作为DEM金字塔的第1层，依次类推，最终建立N层DEM金字塔，最顶层记为第N-1层；

步骤3，从DEM金字塔的顶层开始，进行自适应多角度地形断面高程极值法提取，并经过后处理得到山脊线和山谷线，然后逐层精化提取结果；包括以下子步骤，

步骤3.1，对第N-1层DEM利用自适应多角度地形断面高程极值法和后处理进行提取，得到第N-1层的提取结果；令当前层标记i＝N-2；

步骤3.2，对第i层DEM利用自适应多角度地形断面高程极值法进行提取，将第i+1级的提取结果映射到第i级DEM上，并与第i级DEM上的自适应多角度地形断面高程极值法提取结果进行叠加，然后进行后处理，得到第i级DEM上的提取结果；

步骤3.3，判断是否i＝0，是则输出结果，否则令i＝i-1，返回步骤3.2；

所述自适应多角度地形断面高程极值法，是利用一个窗口大小根据地形起伏度自适应变化的模板在多个方向对DEM进行剖面分析，找出各个断面上的极大值点，对DEM中任一格网点，若同时在多个断面中判断为极大值点，则当作山脊线的候选点；找出各个断面上的极小值点，对DEM中任一格网点，若同时在多个断面中判断为小值点，则当作山谷线的候选点。

而且，所述后处理包括形态学闭操作和Hilditch算法细化。

而且，所述窗口大小根据地形起伏度自适应变化，实现方式为首先将DEM均匀分块为公里格网，然后通过求各分块的地形起伏度，在进行剖面分析时，根据正在分析的格网点所在分块的地形起伏度及DEM分辨率计算该点的窗口的大小，

计算地形起伏度如下式，

R_m，n＝Max_m，n-Min_m，n

其中，R_m，n为地形起伏度，m，n分别为当前公里格网的行列编号，Max_m，n为当前公里格网内部高程最大值，Min_m，n为当前公里格网内部高程最小值；

计算窗口大小如下式，

其中，w_m，n为适宜窗口大小，L_m，n，W_m，n为DEM格网的实际长度和宽度，R_m，n为地形起伏度，GSD_Dem为DEM的分辨率，m，n分别为当前公里格网的行列编号。

而且，所述在多个方向对DEM进行剖面分析，包括对DEM中每个格网点在0°、45°、90°、135°四个方向上进行剖面分析。

而且，金字塔的总层数N按下式确定，

N＝max((log₂(min(width，height)))-o_min-3，1)

其中width为DEM的横向格网数目，height为DEM的纵向格网数目，o_min为预设参数。

本发明还提供一种基于尺度空间的自适应山谷山脊线提取系统，包括以下模块，

初始化模块，用于首先对原始DEM过采样作为处理的初始DEM，对初始DEM进行从低到高地不同程度的高斯平滑得到一组DEM的尺度空间，称为第0组尺度空间；

DEM金字塔建立模块，用于根据金字塔的总层数N，先挑选第0组尺度中空间尺度最大的DEM作为DEM金字塔的最底层，记为第0层，然后将第0层DEM进行降采样，并相应建立第1组尺度空间，再从中挑选尺度最大的DEM作为DEM金字塔的第1层，依次类推，最终建立N层DEM金字塔，最顶层记为第N-1层；

提取模块，用于从DEM金字塔的顶层开始，进行自适应多角度地形断面高程极值法提取，并经过后处理得到山脊线和山谷线，然后逐层精化提取结果；包括以下单元，

初始提取单元，用于对第N-1层DEM利用自适应多角度地形断面高程极值法和后处理进行提取，得到第N-1层的提取结果；令当前层标记i＝N-2；

叠加单元，用于对第i层DEM利用自适应多角度地形断面高程极值法进行提取，将第i+1级的提取结果映射到第i级DEM上，并与第i级DEM上的自适应多角度地形断面高程极值法提取结果进行叠加，然后进行后处理，得到第i级DEM上的提取结果；

迭代单元，用于判断是否i＝0，是则输出结果，否则令i＝i-1，命令叠加单元工作；

所述自适应多角度地形断面高程极值法，是利用一个窗口大小根据地形起伏度自适应变化的模板在多个方向对DEM进行剖面分析，找出各个断面上的极大值点，对DEM中任一格网点，若同时在多个断面中判断为极大值点，则当作山脊线的候选点；找出各个断面上的极小值点，对DEM中任一格网点，若同时在多个断面中判断为小值点，则当作山谷线的候选点。

而且，所述后处理包括形态学闭操作和Hilditch算法细化。

而且，所述窗口大小根据地形起伏度自适应变化，实现方式为首先将DEM均匀分块为公里格网，然后通过求各分块的地形起伏度，在进行剖面分析时，根据正在分析的格网点所在分块的地形起伏度及DEM分辨率计算该点的窗口的大小，

计算地形起伏度如下式，

R_m，n＝Max_m，n-Min_m，n

其中，R_m，n为地形起伏度，m，n分别为当前公里格网的行列编号，Max_m，n为当前公里格网内部高程最大值，Min_m，n为当前公里格网内部高程最小值；

计算窗口大小如下式，

其中，w_m，n为适宜窗口大小，L_m，n，W_m，n为DEM格网的实际长度和宽度，R_m，n为地形起伏度，GSD_Dem为DEM的分辨率，m，n分别为当前公里格网的行列编号。

而且，所述在多个方向对DEM进行剖面分析，包括对DEM中每个格网点在0°、45°、90°、135°四个方向上进行剖面分析。

而且，金字塔的总层数N按下式确定，

N＝max((log₂(min(width，height)))-o_min-3，1)

其中，width为DEM的横向格网数目，height为DEM的纵向格网数目，o_min为预设参数。

与现有的山脊线和山谷线的提取方法相比，本发明的有益效果是：

1)本发明与现有方法相比，克服了现有方法对地形噪声敏感的缺点，更能兼顾山脊(谷)提取的整体趋势及细节变化，更能够保证提取的精度。这是因为本发明采用的尺度空间变换，通过建立尺度空间与高斯金字塔，在多个尺度上对DEM数据进行处理，从粗到细，逐层精化。本发明中算法与局部算法的算法复杂度均为O(n)，而整体算法的复杂度为O(n²)，效率优于整体算法，同时本发明提取结果又避免了局部算法对噪声敏感的缺陷。在保证整体趋势提取效果的也避免了地形图中噪声的影响，同时保证了计算效率。

2)本发明采用的自适应多角度地形断面高程极值法，采用自动化计算剖面窗口大小，减少了人工干预的影响，自动化程度高，具有更好的自适应性，并从多个角度进行剖面分析，能够有效避免地形断面高程极值法因为剖面角度而产生的遗漏问题，使得提取结果更加完整。

3)本发明提取的山脊线和山谷线的算法在理论上正确，实际应用中可行，能够有效的从数字格网高程模型中提取到山脊线和山谷线，提取结果与实际地形变化相一致，具有较大的实际应用价值。

附图说明

图1为本发明实施例的流程图。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

采用真实的DEM数据对基于尺度空间及多角度地形断面的山谷山脊线提取方法进行试验，提取得到的山脊线和山谷线与地形起伏基本一致，证明了本发明在理论上是正确有效的，并且在实际应用中也是切实可行的。其效果见附图，可以看出本发明提取得到的山脊线和山谷线与实际地形在整体和细节上均有较好的一致性。

本发明的实施例是对浙江省某区域的5m分辨率，大小为1673x1909的DEM数据进行山脊线和山谷线的提取，参照图1，本发明实施例的具体步骤如下：

步骤a，尺度空间的生成：首先对原始DEM过采样作为新的初始DEM进行处理，对其进行从低到高地不同程度的高斯平滑得到一组DEM的尺度空间，称为第0组尺度空间。

将原始DEM过采样是为了防止因高斯平滑而导致的DEM细节信息的丢失，更好的保留原始DEM信息。

尺度空间理论的目的是模拟图像数据的多尺度特征，DEM也可以视为一种特殊的图像数据，DEM每个格网点的高程与图像中各像素的灰度值相对应。尺度空间方法的基本思想是：在视觉信息处理模型中引入一个被视为尺度的参数，通过连续变化尺度参数获得不同尺度下的视觉处理信息，然后综合这些信息以深入地挖掘图像的本质特征。尺度空间方法将传统的单尺度视觉信息处理技术纳入尺度不断变化的动态分析框架中，因此更容易获得图像的本质特征。通过尺度空间的建立也能够获得在尺度不断变化中得到从整体到细节的山脊线和山谷线，能够提高提取的完整性和准确性。高斯卷积核是实现尺度变换的唯一线性核，于是一幅二维图像(DEM)的尺度空间定义为：

其中G(x，y，σ)是尺度可变高斯函数，(x，y)是空间坐标，σ是尺度坐标，σ的大小决定图像(DEM)的平滑程度，大尺度对应图像(DEM)的概貌特征，小尺度对应图像(DEM)的细节特征。大的σ值对应粗糙尺度(低分辨率)，反之，小的σ值对应精细尺度(高分辨率)。

实施例中，首先对DEM数据利用双线性内插进行2倍过采样，得到新的DEM的数据。这样可以在更细的层次上进行高斯平滑，更好的保留原始信息，提高提取精度。这里将DEM视为一种特殊的影像，DEM格网点高程对应于影像的每个像素灰度值。双线性内插值算法描述如下：

对于一个目的像素，设置坐标通过反向变换得到的浮点坐标为(i+u，j+v)(其中i、j均为浮点坐标的整数部分，u、v为浮点坐标的小数部分，是取值[0，1)区间的浮点数)，则这个像素的值f(i+u，j+v)可由原DEM中坐标为(i，j)、(i+1，j)、(i，j+1)、(i+1，j+1)所对应的周围四个像素的值决定，即：

f(i+u，j+v)＝(1-u)(1-v)f(i，j)+(1-u)vf(I，j+1) (2)

+u(1-v)f(i+1，j)+uvf(i+1，j+1)

其中f(i，j)表示源DEM(i，j)处的格网点高程值，以此类推。

然后对其进行从低到高的不同程度S次高斯平滑得到第0组尺度空间，具体实施时本领域技术人员可自行预设S的取值。实施例中S取3，计算方法见公式(1)，尺度可变高斯函数具体公式如下：

其中，

其中G(x，y，σ)是尺度可变高斯函数，I(x，y)代表DEM在空间坐标(x，y)处的高程值，(x，y)代表空间坐标，σ是尺度坐标，S为高斯滤波次数，σ₀为初始值。

对于第0组尺度空间，3次高斯平滑σ的大小分别为2σ₀，后续操作中每层金字塔都要建立相应的尺度空间，在步骤b中将详细介绍σ取值方法。

步骤b，DEM金字塔的建立：选择尺度空间中一定尺度的DEM降采样并重复建立尺度空间，选择一定尺度的DEM进行降采样得到上层DEM，最终建立N层DEM金字塔。

具体可为挑选第0组尺度空间尺度最大(σ值最大)的DEM作为DEM金字塔的最底层即第0层，然后将第0层DEM进行降采样，并相应建立第1组尺度空间，再从中挑选尺度最大的DEM作为DEM金字塔的第1层，后面每层DEM都是下一层DEM进行降采样并建立尺度空间选择尺度最大的DEM的结果，最终建立N层DEM金字塔，最顶层记为第N-1层。建立DEM金字塔的目的有两个，一是为提高山脊(谷)线提取的效率；二是减少了DEM中地形噪声的干扰，有利于提取山脊(谷)线的主体部分。

实现过程如下：

首先挑选第0组尺度空间中尺度最大即高斯平滑的σ值最大的DEM作为DEM金字塔的最底层即第0层，然后将第0层DEM利用双线性内插方法进行2倍降采样(和2倍过采样时相反过程，内插方法见步骤a)，之后利用尺度变换建立第1组尺度空间，再从中挑选尺度最大的DEM作为DEM金字塔的第1层，后面每层DEM都是下一层DEM进行2倍降采样并建立尺度空间后选择尺度最大的DEM的结果。可采用以下迭代过程实现：

1)初始化t＝0，挑选第0组尺度空间中尺度最大的DEM作为DEM金字塔的最底层；

2)对第t层DEM利用双线性内插方法进行2倍降采样，利用尺度变换建立第t+1组尺度空间，从中挑选尺度最大的DEM作为DEM金字塔的第t+1层；

3)判断是否t＝N-1，是则输出结果，否则令t＝t+1，返回2)。

具体实施时，本领域技术人员可自行预设金字塔的总层数N。考虑到当DEM格网小于一定值(一般为256)，可以不必继续重采样，建议金字塔的总层数N计算方法如下：

N＝max((log₂(min(width，height)))-o_min-3，1) (5)

其中width为DEM的横向格网数目，height为DEM的纵向格网数目，o_min的取值也可由本领域技术人员预先设定，实施例中o_min＝5。

在建立DEM金字塔过程中，生成的每组尺度空间的σ值的确定方法如下：

其中σ₀为初始值，具体实施时，本领域技术人员可自行预设取值；o为金字塔中的第几层，S为一组尺度空间内影像的个数(和高斯滤波次数一致)，s用于标识一组尺度空间中的第几张图像。则最终金字塔中从下到上每层的尺度分别为2^o+1σ₀。实际应用中可选择直接生成最终尺度的DEM。这里在选定第0层金字塔的DEM后，直接对其进行双线性内插2倍降采样，并一步进行σ＝σ(o)的高斯平滑得到上一层金字塔DEM，再重复进行降采样，高斯平滑直到得到最终的DEM金字塔。σ(o)的确定方法如下：

σ(o)＝2^o+1σ₀ (7)

其中σ₀为初始值，o为金字塔中的第几层。

步骤c，从DEM金字塔的顶层开始，进行自适应多角度地形断面高程极值法提取：对DEM利用自适应多角度地形断面高程极值法进行山脊线和山谷线的提取，并经过一系列后处理将得到山脊线和山谷线。

本发明在现有地形断面高程极值法基础上进行自适应和多角度方面的改进，故命名为自适应多角度地形断面高程极值法。自适应多角度地形断面高程极值法是利用一个窗口大小可变的模板在多个方向对DEM进行剖面分析，找出各个断面上的极大(小)值点，对DEM中任一格网点，若同时在多个断面中判断为极大(小)值点，则将它当作山脊(谷)线的候选点。然后可以进行后处理，包括进行对提取结果进行形态学闭操作，最后利用Hilditch算法细化得到山脊(谷)线。

在提取中采用多个方向对DEM进行剖面分析是因为当剖面方向与山脊(谷)延伸方向成90°夹角时，判断的效果最好，当剖面方向与山脊(谷)延伸方向成0°夹角时，剖面分析法失效，而DEM中山脊(谷)的延伸方向可能是任意方向，进行剖面角度选择时应保证角度选择的丰富性、均匀性、对称性，这样在进行多角度的剖面分析才能够有效避免提取中因为剖面角度问题而产生的遗漏现象，才能够保证提取结果的完整性。

在进行自适应多角度地形断面高程极值法提取之前首先要确定进行剖面分析的窗口大小。剖面分析的适宜窗口大小与DEM分辨率及区域的起伏度有关，地形起伏越剧烈，剖面分析窗口应越小，DEM的分辨率越低，即每个DEM格网代表的实际面积越大，窗口应越小，为更好的反应地形起伏程度，进行剖面分析时窗口的大小在本发明中可以通过一种自动化方式计算得到，本发明中首先将DEM均匀分块，然后通过求各分块的地形起伏度，在进行剖面分析时，根据正在分析的格网点所在分块的地形起伏度及DEM分辨率计算该点适宜的分析窗口的大小。

在本实施例中采用划分公里格网计算地形起伏度的方法：

首先将DEM分为公里格网(即格网尺寸为1×1公里)，然后对每个公里格网内的DEM分别计算地形起伏度，并求得地形起伏度，公式如下：

R_m，n＝Max_m，n-Min_m，n (8)

其中R_m，n为地形起伏度，m，n分别为当前公里格网在划分的格网中的行列编号，Max_m，n为当前公里格网内部高程最大值，Min_m，n为当前公里格网内部高程最小值。

继而求得适宜窗口大小为：

其中w_m，n为适宜窗口大小，L_m，n，W_m，n为DEM格网的实际长度和宽度(实施例为1公里)，R_m，n为地形起伏度，GSD_Dem为DEM的分辨率，m，n分别为当前公里格网在划分的整个公里格网中的行列编号。因此，实际上一个公里格网内所有DEM格网点可以采用同样的窗口大小。适宜窗口大小必须是大于1的奇数，如果窗口大小是3则很容易受到噪声影响，5则能一定程度上克服噪声，所以实施例定义5为最小的适宜窗口大小。

然后对DEM金字塔最顶层的DEM的每个格网点在0°、45°、90°、135°四个方向上进行剖面分析，若其在两个或以上是剖面中判断为极大(小)值点，且与当前剖面的最小(大)值的差值达到5m以上，则将它当作山脊(谷)线的候选点，然后对提取结果进行二值化生成二值图，候选点所在像素的值为1，其他地方值为0，之后进行形态学闭操作，最后利用Hilditch算法细化并跟踪提取得到矢量山脊(谷)线。

闭操作和Hilditch细化算法为已有技术，为便于本领域技术人员实施应用到本发明起见，以下进行介绍。

闭操作为形态学处理的常用的方式之一，相当于膨胀一次再腐蚀一次。闭操作可使轮廓线更光滑，消除狭窄的间断和细长的鸿沟，消除小的空洞，并填补轮廓线中的断裂。本发明从采用闭操作可以填补提取得到的山脊(谷)线中存在的小断裂，使得山脊(谷)线提取结果更加完整。闭操作的定义为：

其中A为二值图，B为结构元素，实施例取3×3的矩形模板，具体实施时本领域技术人员也可以设定为其他尺寸。

Hilditch细化算法是经典的二值图像细化算法，所谓的细化就是经过层层的剥离，从原来的图中去掉一些点，但仍要保持原来的形状，直到得到图像的骨架即为图像中形状的中轴线。Hilditch细化算法的步骤为：

对图像从左向右从上向下迭代每个像素，是为一个迭代周期。在每个迭代周期中，对于每一个像素p，如果它同时满足6个条件，则标记它。在当前迭代周期结束时，则把所有标记的像素的值设为背景值。如果某次迭代周期中不存在标记点(即满足6个条件的像素)，则算法结束。假设背景值为0，前景值为1，则：

6个条件为：

(I)p为1，即p不是背景；

(2)x1，x3，x5，x7不全部为1(否则把p标记删除，图像空心了)；

(3)x1～x8中，至少有2个为1(若只有1个为1，则是线段的端点。若没有为1的，则为孤立点)；

(4)p的8连通联结数为1；

(5)假设x3已经标记删除，那么当x3为0时，p的8联通联结数为1；

(6)假设x5已经标记删除，那么当x5为0时，p的8联通联结数为1。

其中，(5)、(6)是对(4)在特殊情况下的补充，如果满足(5)、(6)，也可认为p的8联通连接数为1。

联结数指在像素p的3×3邻域中，和p连接的图形分量的个数，如下表所示模板：

x4	x3	x2
			x5	p	x1
x6	x7	x8

8连通联结数计算公式是：

其中，x为如上表3×3二值图模板中的值，二值图中的值为0或1。

步骤d，逐层精化提取结果：因为步骤c是对N层DEM金字塔中的第N-1级DEM利用自适应多角度地形断面高程极值法进行提取，本步骤利用自适应多角度地形断面高程极值法对第N-2级DEM进行提取，并将第N-1级的提取结果映射到第N-2级DEM上，即按照地理位置映射，如果N-2级某个格网点对应N-1级上同样位置的格网点被标记为山脊(谷)线，则N-2级上该格网点也标记为山脊(谷)线，并与第N-2级DEM上的提取结果进行叠加，经过形态学闭操作，Hilditch算法细化等后处理操作(可参见步骤c)得到新的山脊(谷)线，作为第N-2级DEM上的提取结果。以此类推，直到得到第0级DEM的上的山脊线和山谷线，即为最终提取结果。

具体实施时，步骤c、d的实现可以设计流程为：

1)对第N-1层DEM利用自适应多角度地形断面高程极值法和后处理进行提取，得到第N-1层的提取结果；令当前层标记i＝N-2，

2)对当前层DEM(第i层)利用自适应多角度地形断面高程极值法进行提取，将第i+1层的提取结果映射到第i级DEM上，并与第i级DEM上利用自适应多角度地形断面高程极值法的提取结果进行叠加，然后进行后处理，得到第i级DEM上最终的提取结果；

3)判断是否i＝0，是则输出结果，即第0级DEM上最终的山脊线和山谷线提取结果；否则令i＝i-1，返回2)。

具体实施时，以上流程可采用计算机软件技术实现自动运行，也可以采用模块化方式提供相应系统。本发明实施例还提供一种基于尺度空间的自适应山谷山脊线提取系统，包括以下模块，

初始化模块，用于首先对原始DEM过采样作为处理的初始DEM，对初始DEM进行从低到高地不同程度的高斯平滑得到一组DEM的尺度空间，称为第0组尺度空间；

DEM金字塔建立模块，用于根据金字塔的总层数N，先挑选第0组尺度中空间尺度最大的DEM作为DEM金字塔的最底层，记为第0层，然后将第0层DEM进行降采样，并相应建立第1组尺度空间，再从中挑选尺度最大的DEM作为DEM金字塔的第1层，依次类推，最终建立N层DEM金字塔，最顶层记为第N-1层；

提取模块，用于从DEM金字塔的顶层开始，进行自适应多角度地形断面高程极值法提取，并经过后处理得到山脊线和山谷线，然后逐层精化提取结果；包括以下单元，

初始提取单元，用于对第N-1层DEM利用自适应多角度地形断面高程极值法和后处理进行提取，得到第N-1层的提取结果；令当前层标记i＝N-2；

叠加单元，用于对第i层DEM利用自适应多角度地形断面高程极值法进行提取，将第i+1级的提取结果映射到第i级DEM上，并与第i级DEM上的自适应多角度地形断面高程极值法提取结果进行叠加，然后进行后处理，得到第i级DEM上的提取结果；

迭代单元，用于判断是否i＝0，是则输出结果，否则令i＝i-1，命令叠加单元工作；

所述自适应多角度地形断面高程极值法，是利用一个窗口大小根据地形起伏度自适应变化的模板在多个方向对DEM进行剖面分析，找出各个断面上的极大值点，对DEM中任一格网点，若同时在多个断面中判断为极大值点，则当作山脊线的候选点；找出各个断面上的极小值点，对DEM中任一格网点，若同时在多个断面中判断为小值点，则当作山谷线的候选点。

各模块具体实现可参见相应步骤，本发明不予赘述。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (8)

1.一种基于尺度空间的自适应山谷山脊线提取方法，其特征在于：包括以下步骤，

步骤1，尺度空间的生成，包括首先对原始DEM过采样作为处理的初始DEM，对初始DEM进行从低到高地不同程度的高斯平滑得到一组DEM的尺度空间，称为第0组尺度空间；

步骤2，DEM金字塔的建立，包括根据金字塔的总层数N，先挑选第0组尺度中空间尺度最大的DEM作为DEM金字塔的最底层，记为第0层，然后将第0层DEM进行降采样，并相应建立第1组尺度空间，再从中挑选尺度最大的DEM作为DEM金字塔的第1层，依次类推，最终建立N层DEM金字塔，最顶层记为第N-1层；

步骤3，从DEM金字塔的顶层开始，进行自适应多角度地形断面高程极值法提取，并经过后处理得到山脊线和山谷线，然后逐层精化提取结果；包括以下子步骤，

步骤3.1，对第N-1层DEM利用自适应多角度地形断面高程极值法和后处理进行提取，得到第N-1层的提取结果；令当前层标记i＝N-2；

步骤3.2，对第i层DEM利用自适应多角度地形断面高程极值法进行提取，将第i+1级的提取结果映射到第i级DEM上，并与第i级DEM上的自适应多角度地形断面高程极值法提取结果进行叠加，然后进行后处理，得到第i级DEM上的提取结果；

步骤3.3，判断是否i＝0，是则输出结果，否则令i＝i-1，返回步骤3.2；

所述自适应多角度地形断面高程极值法，是利用一个窗口大小根据地形起伏度自适应变化的模板在多个方向对DEM进行剖面分析，找出各个断面上的极大值点，对DEM中任一格网点，若同时在多个断面中判断为极大值点，则当作山脊线的候选点；找出各个断面上的极小值点，对DEM中任一格网点，若同时在多个断面中判断为小值点，则当作山谷线的候选点；

所述窗口大小根据地形起伏度自适应变化，实现方式为首先将DEM均匀分块为公里格网，然后通过求各分块的地形起伏度，在进行剖面分析时，根据正在分析的格网点所在分块的地形起伏度及DEM分辨率计算该点的窗口的大小，

计算地形起伏度如下式，

R_m,n＝Max_m,n-Min_m,n

其中，R_m,n为地形起伏度，m,n分别为当前公里格网的行列编号，Max_m,n为当前公里格网内部高程最大值，Min_m,n为当前公里格网内部高程最小值；

计算窗口大小如下式，

其中，w_m,n为适宜窗口大小，L_m,n,W_m,n为DEM格网的实际长度和宽度，R_m,n为地形起伏度，GSD_Dem为DEM的分辨率，m,n分别为当前公里格网的行列编号。

2.根据权利要求1所示基于尺度空间的自适应山谷山脊线提取方法，其特征在于：所述后处理包括形态学闭操作和Hilditch算法细化。

3.根据权利要求1所示基于尺度空间的自适应山谷山脊线提取方法，其特征在于：所述在多个方向对DEM进行剖面分析，包括对DEM中每个格网点在0°、45°、90゜、135°四个方向上进行剖面分析。

4.根据权利要求1或2或3所示基于尺度空间的自适应山谷山脊线提取方法，其特征在于：金字塔的总层数N按下式确定，

N＝max((log₂(min(width,height)))-o_min-3,1)

其中width为DEM的横向格网数目，height为DEM的纵向格网数目，o_min为预设参数。

5.一种基于尺度空间的自适应山谷山脊线提取系统，其特征在于：包括以下模块，

初始化模块，用于首先对原始DEM过采样作为处理的初始DEM，对初始DEM进行从低到高地不同程度的高斯平滑得到一组DEM的尺度空间，称为第0组尺度空间；

DEM金字塔建立模块，用于根据金字塔的总层数N，先挑选第0组尺度中空间尺度最大的DEM作为DEM金字塔的最底层，记为第0层，然后将第0层DEM进行降采样，并相应建立第1组尺度空间，再从中挑选尺度最大的DEM作为DEM金字塔的第1层，依次类推，最终建立N层DEM金字塔，最顶层记为第N-1层；

提取模块，用于从DEM金字塔的顶层开始，进行自适应多角度地形断面高程极值法提取，并经过后处理得到山脊线和山谷线，然后逐层精化提取结果；包括以下单元，

初始提取单元，用于对第N-1层DEM利用自适应多角度地形断面高程极值法和后处理进行提取，得到第N-1层的提取结果；令当前层标记i＝N-2；

叠加单元，用于对第i层DEM利用自适应多角度地形断面高程极值法进行提取，将第i+1级的提取结果映射到第i级DEM上，并与第i级DEM上的自适应多角度地形断面高程极值法提取结果进行叠加，然后进行后处理，得到第i级DEM上的提取结果；

迭代单元，用于判断是否i＝0，是则输出结果，否则令i＝i-1，命令叠加单元工作；

所述自适应多角度地形断面高程极值法，是利用一个窗口大小根据地形起伏度自适应变化的模板在多个方向对DEM进行剖面分析，找出各个断面上的极大值点，对DEM中任一格网点，若同时在多个断面中判断为极大值点，则当作山脊线的候选点；找出各个断面上的极小值点，对DEM中任一格网点，若同时在多个断面中判断为小值点，则当作山谷线的候选点；

所述窗口大小根据地形起伏度自适应变化，实现方式为首先将DEM均匀分块为公里格网，然后通过求各分块的地形起伏度，在进行剖面分析时，根据正在分析的格网点所在分块的地形起伏度及DEM分辨率计算该点的窗口的大小，

计算地形起伏度如下式，

R_m,n＝Max_m,n-Min_m,n

其中，R_m,n为地形起伏度，m,n分别为当前公里格网的行列编号，Max_m,n为当前公里格网内部高程最大值，Min_m,n为当前公里格网内部高程最小值；

计算窗口大小如下式，

其中，w_m,n为适宜窗口大小，L_m,n,W_m,n为DEM格网的实际长度和宽度，R_m,n为地形起伏度，GSD_Dem为DEM的分辨率，m,n分别为当前公里格网的行列编号。

6.根据权利要求5所示基于尺度空间的自适应山谷山脊线提取系统，其特征在于：所述后处理包括形态学闭操作和Hilditch算法细化。

7.根据权利要求5所示基于尺度空间的自适应山谷山脊线提取系统，其特征在于：所述在多个方向对DEM进行剖面分析，包括对DEM中每个格网点在0°、45°、90゜、135°四个方向上进行剖面分析。

8.根据权利要求5或6或7所示基于尺度空间的自适应山谷山脊线提取系统，其特征在于：金字塔的总层数N按下式确定，

N＝max((log₂(min(width,height)))-o_min-3,1)

其中，width为DEM的横向格网数目，height为DEM的纵向格网数目，o_min为预设参数。