地表坡向约束的高分辨率DEM地形特征增强提取方法

本发明的实施例一般涉及地理空间数据处理领域，并且更具体地，涉及一种地表坡向约束的高分辨率dem地形特征增强提取方法。

背景技术：

数字高程模型(digitalelevationmodel，dem)是表征地形特征的三维产品。随着地球观测技术手段的进步，大规模多分辨率dem被提供给公众使用。在一些表示地形的可视化方法中，曲率、坡度、山体阴影和坡向等参数被用作地形特征提取的参数。

传统方法更适用于“噪音小”且更为光滑的低分辨率的dem产品，但若想从高分辨率dem数据中获得地表特征的细节，传统方法则失去优势。这是因为高分辨率dem产品因为能展示更多的地表粗糙度的细节信息，进而形成“椒盐噪声”，从而带来负面影响。

技术实现要素：

根据本发明的实施例，提供了一种地表坡向约束的高分辨率dem地形特征增强提取方法，该方法包括：

dem生成栅格产品，根据所述栅格产品的栅格像素在其周围均匀分布的8个方向上的相邻像素，计算所述栅格像素的坡向值，将连续的坡向通过离散的栅格进行划分，得到地表坡向的多个值域区间；

设置地表坡向的坡向粒度，以所述坡向粒度的方向为约束条件，建立适用于高分辨率dem地貌特征提取的增强模型；

通过所述增强模型从米级或亚米级的高分辨率dem数据中提取地形特征。

进一步地，所述计算所述栅格像素的坡向值，包括：

a(x0，y0)＝57.29578*atan2(φ1-φ2)

其中，a(x0，y0)为栅格像素(x0，y0)的坡向值，φ1和φ2分别为第一角度和第二角度，通过下式计算得到：

其中，nw、n、ne、w、e、sw、s、se分别表示8个方向，即西北、北、东北、西、东、西南、南、东南。

进一步地，所述坡向粒度为：

其中，h为坡向粒度；n为离散的坡向产品的值域区间个数。

进一步地，所述建立适用于高分辨率dem地貌特征提取的增强模型，包括：

建立坡向约束的基于曲率的模糊分类提取模型；和/或

建立坡向约束的基于曲率的lbp提取模型；和/或

建立坡向约束的基于曲率的slic提取模型。

进一步地，所述建立坡向约束的基于曲率的模糊分类提取模型，包括：

计算任一像素点和其相邻的8个像素之间的坡向差；

判断两个相反方向的坡向差，如果其中一个坡向差大于坡向粒度，另一个坡向差小于坡向粒度，则将所述像素点标记为1；否则，将所述像素点标记为0；生成坡向标记地图；

建立坡向约束的基于曲率的模糊分类提取模型。

进一步地，所述坡向约束的基于曲率的模糊分类提取模型为：

其中，f(xo，y0)是坡向约束的基于曲率的模糊分类提取模型；c(xo，y0)是像素的曲率均值；θ1、θ2是所要提取特征要素各类别的阈值；ad+(xo，yo)和ad-(xo，yo)表示在两个相反方向上的坡向差。

进一步地，所述建立坡向约束的基于曲率的lbp提取模型，包括：

设置高分辨率dem中的3×3的像素窗口，通过三元模式将所述像素窗口的中心像素与其8个相邻像素的值进行比较，得到lbp提取模型：

其中，p0为像素窗口的中心像素值；pa为8个相邻像素中任一个的值；lbp(p0-pa)为坡向约束的基于曲率的lbp提取模型。

进一步地，所述建立坡向约束的基于曲率的slic提取模型，包括：

步骤1：初始化种子点；

步骤2：计算种子点的n*n邻域内的所有像素点的梯度值，将梯度值最小的像素点作为优化后的种子点；

步骤3：搜索优化后的每个种子点在搜索范围内的像素点；

步骤4：计算搜索范围内每个像素点到对应种子点的距离，得到距离度量；所述距离包括属性距离和空间距离；

步骤5：将到所述像素点的距离度量最小的种子点作为该像素点的聚类中心，为每个像素点分配对应的类标签；

步骤6：迭代步骤1～步骤5，直至误差收敛，得到坡向约束的基于曲率的slic提取模型。

进一步地，还包括：

建立标记表，用于对遍历过的像素点进行标记；所述标记表内的元素均为-1，按照从左到右，从上到下的顺序，将不连续的超像素以及尺寸小于尺寸阈值的超像素重新分配给邻近的超像素，并分配标签给遍历过的像素点进行标记，直至遍历所有像素点。

进一步地，所述计算搜索范围内每个像素点到对应种子点的距离，得到距离度量，包括：

其中，dz代表两个像素间的属性距离，dxy代表两个像素间的空间距离，dattr是指两个像素间的高程差，dx和dy分别是指两个像素在水平和垂直维度上的距离；

所述距离度量dtotal为：

其中，θ是指空间距离与属性距离的比率；n为分割后的超级像素的近似数。

应当理解，发明内容部分中所描述的内容并非旨在限定本发明的实施例的关键或重要特征，亦非用于限制本发明的范围。本发明的其它特征将通过以下的描述变得容易理解。

本发明能够针对高分辨率dem地形特征提取的需求，从高分辨率的地形图中降低噪声的负面影响，提高高分辨率dem的利用率，建立地表坡向约束的高分辨率dem地形特征提取模型，增强表征地貌特征的能力。

附图说明

结合附图并参考以下详细说明，本发明各实施例的上述和其他特征、优点及方面将变得更加明显。在附图中，相同或相似的附图标记表示相同或相似的元素，其中：

图1示出了根据本发明的地表坡向约束的高分辨率dem地形特征增强提取方法的流程图；

图2示出了本发明的坡向约束的基于曲率的模糊分类提取模型与现有技术中坡向约束的基于曲率的模糊分类提取方法的提取结果示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的全部其他实施例，都属于本发明保护的范围。

另外，本文中术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，a和/或b，可以表示：单独存在a，同时存在a和b，单独存在b这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。

本发明中，能够针对高分辨率dem地形特征提取的需求，从超高分辨率的地形图中降低噪声的负面影响，提高高分辨率dem的利用率，建立地表坡向约束的高分辨率dem地形特征提取模型，增强表征地貌特征的能力。

图1示出了本发明的地表坡向约束的高分辨率dem地形特征增强提取方法的流程图。

该方法包括：

s100、dem生成栅格产品，根据所述栅格产品的栅格像素在其周围均匀分布的8个方向上的相邻像素，计算所述栅格像素的坡向值，将连续的坡向通过离散的栅格进行划分，得到地表坡向的多个值域区间。

dem是digitalelevationmodel的缩写，表示数字高程模型，是通过有限的地形高程数据实现对地面地形的数字化模拟(即地形表面形态的数字化表达)，它是用一组有序数值阵列形式表示地面高程的一种实体地面模型。

通过dem生成栅格产品，并在由所述栅格产品的栅格数据中，选择一坐标点上的一个像素，例如坐标(xo，yo)上的一个像素；在所述像素周围选择均匀分布的8个方向，例如nw、n、ne、w、e、sw、s和se，分别表示西北、北、东北、西、东、西南、南、东南。所述8个方向均分360度，即每相邻两个方向之间的夹角为45度。在所述8个方向上选择栅格像素的相邻像素。

所述栅格像素的坡向值为：

a(x0，y0)＝57.29578*atan2(φ1-φ2)

其中，a(x0，y0)为栅格像素(x0，y0)的坡向值，φ1和φ2分别为第一角度和第二角度，通过下式计算得到：

其中，nw、n、ne、w、e、sw、s、se分别表示8个方向，即西北、北、东北、西、东、西南、南、东南。

离散的栅格为若干个，包括基于4个、8个或16个方向集的栅格。所述连续的坡向通过离散的栅格划分为若干个值域区间。

在上述实施例中，在西北、北、东北、西、东、西南、南、东南，8个方向上将连续的坡向产品划分为离散的栅格，则此离散坡向产品中的每个像素都被划分为8个值域区间，即：0°～45°、45°～90°、90°～135°、135°～180°、180°～225°、225°～270°、270°～315°、315°～360°。

又例如，在4个方向上将连续的坡向产品划分为离散的栅格，则此离散的坡向产品中的每个像素都被划分为4个阈值区间，即：北向区间对应于0°至45°和315°至360°的范围；东向区间为45°至135°；南向区间为135°至225°；西向区间为225°至315°。

值得注意的是，本发明中出现的高分辨率为垂直分辨率不小于720。

s200、设置地表坡向的坡向粒度，以所述坡向粒度的方向为约束条件，建立适用于高分辨率dem地貌特征提取的增强模型。

所述地表坡向的坡向粒度为：

其中，h为坡向粒度；n为离散的坡向产品的值域区间个数。

作为本发明的一种实施例，当地表坡向被划分为4个值域区间时，坡向粒度为90°；又例如，当地表坡向被划分为8个值域区间时，坡向粒度为45°；又例如，当地表坡向被划分为16个值域区间时，坡向粒度为22.5°。

除了空间分辨率和地图的尺寸会影响地形地貌的特征表现，地表坡向的尺度划分也会影响地形表现。故，建立坡向粒度以表征坡向栅格产品中两个像素之间可以表示的最小方向差异。不同坡向粒度的条件下所表示的地表特征不同。通过增加方向粒度，可以获得更多的地形特征和地表细节。

所述建立适用于高分辨率dem地貌特征提取的增强模型，包括：

建立坡向约束的基于曲率的模糊分类提取模型；和/或

建立坡向约束的基于曲率的lbp提取模型；和/或

建立坡向约束的基于曲率的slic提取模型。

所述建立坡向约束的基于曲率的模糊分类提取模型，包括：

首先，计算任一像素点和其相邻的8个像素之间的坡向差。

作为本发明的一种实施例，从坡向栅格中选定一像素点(x0，y0)，预设该像素点的坡向值a(x0，y0)；计算像素点(x0，y0)和其相邻的8个像素之间的坡向差，用addir表示，其中，dir表示8个方向的索引值，例如，像素点(x0，y0)和其相邻的东向像素点的坡向差，表示为adeast。

其次，判断两个相反方向的坡向差，如果其中一个坡向差大于坡向粒度，另一个坡向差小于坡向粒度，则将所述像素点标记为1；否则，将所述像素点标记为0；生成坡向标记地图。

作为本发明的一种实施例，在与像素点(x0，y0)相邻的8个像素中存在4对互为相反方向的像素，分别为南和北、东和西、东北和西南、东南和西北。8个像素的坡向粒度为22.5°。

在上述互为相反方向的像素对中，如果其中一个像素对的坡向差小于坡向粒度，例如东和西方向上的adeast和adwest，满足一个大于22.5°，且另一个小于22.5°的条件，则该像素(x0，y0)标记为1；例如，adeast且adwest小于22.5，则a(x0，y0)标记为1。

如果不满足上述条件，即互为相反方向的像素对中，两个像素的坡向差均小于或均大于坡向粒度，则将像素(x0，y0)标记为0。

依据上述标记规则，对坡向栅格中的每个栅格像素进行标记0或1，生成坡向标记地图。

最后，建立坡向约束的基于曲率的模糊分类提取模型。

其中，f(xo，y0)是坡向约束的基于曲率的模糊分类提取模型；c(xo，y0)是像素的曲率均值；θ1、θ2是所要提取特征要素各类别的阈值；ad+(xo，yo)和ad-(xo，yo)表示在两个相反方向上的坡向差。

如图2所示为本发明的坡向约束的基于曲率的模糊分类提取模型与现有技术中坡向约束的基于曲率的模糊分类提取方法的提取结果示意图；其中图2a-1和图2a-2是利用现有技术中坡向约束的基于曲率的模糊分类提取方法的提取的火山口边界结果，而图2b-1和图2b-2是利用本发明的改进的坡向约束的基于曲率的模糊分类提取模型进行提取的火山口边界示意图，能够看出，在中低分辨率dem的脊线和山谷线提取中，平均曲率通常提取出正确的结果。然而，高分辨率dem中的地表细节将导致仅用平均曲率会提取过多非脊线或山谷线。

综上，面对高分辨率曲率图上复杂的地形，通过坡向值，可以更容易地将这两个火山口和山谷线的边缘与背景区分开来。因此，将坡向和平均曲率相结合的方法可以处理高分辨率dem中的噪声和接近脊线和山谷线的特征，有效地提高基于曲率的性能。

lbp(localbinarypatterns)又叫做局部二值模式，是一种有效的纹理描述算子，度量和提取图像局部的纹理信息，对光照具有不变性。本发明对现有的lbp算子进行了改进，建立了坡向约束的基于曲率的lbp提取模型。

所述建立坡向约束的基于曲率的lbp提取模型，过程包括：

首先，设置高分辨率dem中的3×3的像素窗口；

其次，通过三元模式将所述像素窗口的中心像素与其8个相邻像素的值进行比较，得到lbp提取模型。

所述三元模式即1，0，-1，通过1，0，-1的三元模式来表示地形特征。

所述lbp提取模型为：

其中，p0为像素窗口的中心像素值；pa为8个相邻像素中任一个的值；lbp(p0-pa)为坡向约束的基于曲率的lbp提取模型。

在高分辨率dem中获得地表细节信息时，可以通过高程差很好地区分山脊、山肩或火山口边缘。在理想情况下，火山口边缘(山脊线)与肩部相比有明显的高差。然而，面对高分辨率dem中复杂的地形，基于高程的lbp无法预先定义精确的模式模板来表示复杂地形。

相比之下，尺度为1m的dem中的地表要复杂得多，并且与尺度为30m的dem中的地表差异显著。这说明在高分辨率dem上精确地表示复杂的地表高程变化对lbp是一个挑战。

当使用不同的dem导数时，像lbp这样的像素级方法的性能差异很大。与基于高程的lbp方法相比，结合坡向差异的lbp可以从高分辨率dem中提取更多的脊线。

slic(simplelineariterativeclustering)为简单的线性迭代聚类，是一种图像分割算法。本发明对现有的slic算法进行了改进，建立了坡向约束的基于曲率的slic提取模型。

所述建立坡向约束的基于曲率的slic提取模型，过程包括：

步骤1：初始化种子点。

按照设定的超像素个数，在图像内均匀的分配种子点。

作为本发明的一种实施例，设图片共有n个像素点，预分割为k个相同尺寸的超像素，则每个超像素的大小为n/k，则相邻种子点的距离(步长)近似为s＝sqrt(n/k)。

步骤2：计算种子点的n*n邻域内的所有像素点的梯度值，将梯度值最小的像素点作为优化后的种子点。

作为本发明的一种实施例，取n为3，则计算种子点在3*3邻域内的所有像素点的梯度值，得到邻域内梯度值最小的像素值，作为优化后的种子点。

步骤3：搜索优化后的每个种子点搜索范围内的像素点。

slic的搜索范围限制为2s*2s，s为相邻种子点的距离。即搜索范围是以每个种子点为中心的2s*2s的矩形范围。

步骤4：计算邻域内每个像素点到对应种子点的距离，得到距离度量；所述距离包括属性距离和空间距离；

对于每个搜索到的像素点，分别计算其到该种子点的距离，即计算两个像素间的属性距离dz和空间距离dxy：

其中，dz代表两个像素间的属性距离(例如高程、曲率或坡向)，dxy代表两个像素间的空间距离，dattr是指两个像素间的高程差，dx和dy分别是指两个像素在水平和垂直维度上的距离；

所述距离度量dtotal为：

其中，θ是指空间距离与属性距离的比率；θ值越大，则会使得结果包含更大尺寸的超级像素；n为分割后的超级像素的近似数。

步骤5：将到所述像素点的距离度量最小的种子点作为该像素点的聚类中心，为每个像素点分配对应的类标签。

当一个像素点被多个种子点搜索到时，通过计算能够得出该像素点与每个种子点之间的距离度量。选择其中距离度量最小的对应种子点作为该像素点的聚类中心。并将对应聚类中心的类标签分配给该像素点。

步骤6：迭代步骤1～步骤5，直至误差收敛，得到坡向约束的基于曲率的slic提取模型。

作为本发明的一种实施例，在上述步骤中，可以通过建立标记表增强连通性。具体包括：

建立标记表，用于对遍历过的像素点进行标记；所述标记表内的元素均为-1，按照从左到右，从上到下的顺序，将不连续的超像素以及尺寸小于尺寸阈值的超像素重新分配给邻近的超像素，并分配标签给遍历过的像素点进行标记，直至遍历所有像素点。尺寸阈值依实际情况进行选取。

基于高程和基于曲率的分割方法不能用以提取山脊线，因为其主要用于检测具有显著高程变化的边缘地块，而这些地貌大部分是山肩而不是山脊。坡向约束的基于曲率的slic提取模型相较于传统的基于高程和基于曲率的提取方法能够更好地识别山脊线，提升表征地貌的能力。

s300、通过所述增强模型从米级或亚米级的高分辨率dem数据中提取地形特征。

将待分析的米级或亚米级高分辨率dem数据作为数据源，选择上述坡向约束的基于曲率的模糊分类提取模型、坡向约束的基于曲率的lbp提取模型或坡向约束的基于曲率的slic提取模型作为分析模型，输出地形特征表现。

根据本发明的实施例，能够针对高分辨率dem地形特征提取的需求，从高分辨率的地形图中降低噪声的负面影响，提高高分辨率dem的利用率，建立地表坡向约束的高分辨率dem地形特征提取模型，增强表征地貌特征的能力。

需要说明的是，对于前述的各方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本发明并不受所描述的动作顺序的限制，因为依据本发明，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于可选实施例，所涉及的动作和模块并不一定是本发明所必须的。