一种基于数字高程模型的坡形提取与识别方法

摘要

本发明是一种基于数字高程模型的坡形提取与识别方法，目的是提供一种利用数字高程模型作为数据源提取各种坡形分类信息的方法。针对如今以数字高程模型为数据源提取坡形过程中的各种困难，提出了在划分斜坡单元的基础上，利用做等高线垂线的方法找到斜坡的最陡路径，以此得到的纵剖面线可以较好的反映出斜坡的真实坡形。之后利用三次函数拟合的方法，分为多种情况对坡形进行了识别，将主观过程转化为了数学计算，既能保证贴近真实情况还能很好的把握住整体坡形。本发明较好地克服了坡形提取与识别过程中不确定性因素影响和难以客观量化的困难，提出了一套切实可行的提取与识别方法，具有一定的实际应用价值。

说明书

技术领域

本发明涉及地理信息领域。

背景技术

坡形是一种重要的地形特征，也是影响滑坡发育的决定性因素， 因此无论是在滑坡预测还是在土地评价中，都会涉及到坡形的提取。 而如今最常用的地形信息获取数据源就是数字高程模型(以下简称 DEM)，它是地表形态高程属性的数字化表达。

目前，随着GIS技术的飞速发展，利用DEM作为信息源，提取 各种地形定量因子，例如坡度、坡向、地表曲率、地形起伏度、地 表粗糙度等，已经成为十分成熟的技术。然而，坡形的提取却由于 其无法准确定量化计算且存在大量不确定性因素等原因一直被研究 者所回避。

发明内容

本发明的目的是提供一种利用DEM作为数据源提取各种坡形分 类信息的方法。

为实现本发明目的，其具体的技术方案为：

一种基于数字高程模型的坡形提取与识别方法，包括如下步骤：

步骤1.将所测区域的DEM利用地表水文分析提取出的汇水线与 分水线(即山谷和山脊线)所形成的集水流域，划分为若干个斜坡 单元；

步骤2.对于每一个斜坡单元，从单元内所有高程数据点中找出 高程最大的一个点，当高程最大的点不止一个时可任取一个；

步骤3.在每一个斜坡单元内，利用单元内所有数据点的高程值 做出斜坡在单元范围内的所有等高线，数据点不足时需要进行插值 处理，等高线越密最终所得到的坡形结果越准确；

步骤4.以步骤2中得出的高程最大点为起点，从高到低依次找 出下一条等高线上的一个点，使之与找出的上一点距离最近,所有找 到点的连线可称为斜坡的一条最陡路径；

步骤5.将步骤4中找到的每个点放置于直角坐标系中，以每个 点的高程值作为纵坐标，每两个相邻点的横坐标间隔为它们之间的 水平距离差。连接直角坐标系中各点即可得到可以合理表达坡形的 斜坡纵剖面曲线；

步骤6.将步骤5中直角坐标系内的所有点拟合成一个三次函 数曲线的一部分，并找出这个三次函数的极值点与拐点；

步骤7.如果步骤6所拟合出的曲线上靠近两端点的某一指定范 围内存在极值点，则舍去极值点与端点之间所有的点；

步骤8.经步骤7处理后的曲线上，靠近中点的某一指定范围内 如果存在拐点，则以拐点为分界点，判断：(1)如果拐点以上部分的 两个端点纵坐标平均值与中点纵坐标值之差大于某一特定参数R(R 为正)，拐点以下部分的中点纵坐标值与两个端点纵坐标平均值之 差大于参数R，则识别结果为上凹下凸形坡；(2)如果拐点以上部分 的两个端点纵坐标平均值与中点纵坐标值之差大于参数R，拐点以下 部分的中点纵坐标值与两个端点纵坐标平均值之差不大于参数R，则 识别结果为凹形坡；(3)如果拐点以上部分的中点纵坐标值与两个端 点纵坐标平均值之差大于参数R，拐点以下部分的两个端点纵坐标平 均值与中点纵坐标值之差大于参数R，则识别结果为上凸下凹形坡； (4)如果拐点以上部分的中点纵坐标值与两个端点纵坐标平均值之 差大于参数R，拐点以下部分的两个端点纵坐标平均值与中点纵坐标 值之差不大于参数R，则识别结果为凸形坡；(5)如果拐点以上部分 的中点纵坐标值与两个端点纵坐标平均值之差的绝对值不大于参数 R，拐点以下部分的两个端点纵坐标平均值与中点纵坐标值之差大于 参数R，则识别结果为凹形坡；(6)如果拐点以上部分的中点纵坐标 值与两个端点纵坐标平均值之差的绝对值不大于参数R，拐点以下部 分的中点纵坐标值与两个端点纵坐标平均值之差大于参数R，则识别 结果为凸形坡；(7)如果拐点以上部分的中点纵坐标值与两个端点纵 坐标平均值之差的绝对值不大于参数R，拐点以下部分的中点纵坐标 值与两个端点纵坐标平均值之差的绝对值也不大于参数R，则识别结 果为平坡；

步骤9.经步骤7处理后的曲线上，靠近中点的某一指定范围内 如果不存在拐点，则判断：(1)如果两个端点纵坐标平均值与中点纵 坐标值之差大于某一特定参数r(r为正)，则识别结果为凹形坡； (2)如果中点纵坐标值与两个端点纵坐标平均值之差大于参数r，则 识别结果为凸形坡；(3)如果中点纵坐标值与两个端点纵坐标平均值 之差的绝对值不大于参数r，则识别结果为平坡。

本发明的有益效果在于：

(1)在划分了斜坡单元的基础上提取坡形，而非网格单元，使得 最终得出的坡形分布结果更加准确，避免了一个单元内包含多个斜 坡的情况；

(2)利用等高线的原理找到斜坡的最陡路径，以此获得纵剖面曲 线，而非以直线路径获得纵剖面曲线，这样可以比较合理的反映出 斜坡的真实坡形；

(3)运用三次函数拟合的方法识别坡形，可以在保证贴近真实情 况的同时，较好地将主观过程转化为数学计算，并且可以较好的把 握住斜坡的整体坡形，而不被小范围内的坡形变化所影响；

(4)本发明可以完整地识别出出各种坡形情况，无论简单的凹凸 平坡还是复杂的“S”形坡。还可以根据需要设定参数来主观控制识 别过程，这些都使得最终识别结果更加准确。

附图说明

图1为本发明的流程图。

图2为本发明的计算机操作流程图。

图3为本发明提供的方法对某斜坡生成的最陡路径示意图。

图4为本发明提供的方法对某斜坡生成的纵剖面线及坡形识别 结果示意图。

图5为本发明提供的方法在某村镇提取的坡形分类示意图。

具体实施方式

下面结合附图和实施例对本发明的具体实施方法作进一 步说明，但不应该理解为本发明范围仅限于下述实施例。在不脱离 本发明上述方法的情况下，根据本领域普通技术知识和惯用手段， 做出各种替换和变更，均应包括在本发明的保护范围内。

如图2所示，为利用本发明的实施方法编织出的计算机程序流 程图，具体实施步骤包括：

(1)利用GIS软件中的水文分析工具，分别提取出正地形和反地 形的汇水线与分水线，然后把生成的集水流域与反向集水流域融合， 再经后期修改不合理的单元，最终得到所测区域的斜坡单元划分结 果；

(2)将GIS软件中DEM的高程值与xy坐标值均提取至鱼网点， 同时将每个鱼网点所处斜坡单元的编号一并提取至鱼网点，最后将 鱼网点数据导出；

(3)由于GIS软件无法处理较为复杂的计算过程，因此需要借助 二次开发编程实现，可以选择Python、C等各种计算机语言。本实 施例选择应用Matlab计算软件。将导出的数据通过read函数导入 至Matlab；

(4)应用循环语句依次提取出每个单元编号下所有的鱼网点数 据，之后的步骤5到12均是以一个单元内的点为操作对象进行的；

(5)找出单元内所有点中高程值最大的点，当找到多个时可任取 一个作为最陡路径的起点；

(6)利用Matlab中的griddata函数，将单元内的若干个散点差 值计算形成一个光滑的曲面，并利用contour3函数做出曲面的等高 线，等高线条数可以任意指定，等高线越密最终所得到的坡形结果 越准确。同时删除掉在实际高程范围以外的等高线，这些多余的等 高线是由于差值形成曲面时单元的边界被扩大所导致的；

(7)由于Matlab中的contour3函数在生成等高线后所反馈回来 的等高线信息无法直接使用，其中错乱的包含了等高线的条数、每 条等高线的高程值、每条等高线上各点的xy坐标值等各种数据，因 此需要单独提取并整理出对我们有用的信息，即每条等高线上点的 xy坐标值；

(8)从最高点出发，从高到低依次找出下一条等高线上的一个 点，使之与找出的上一点距离最近。之后运用scatter3函数画出单 元内点的空间分布，并在同一幅图中运用plot3函数画出所有找到 的点的连线；

(9)建立直角坐标系xy平面，取水平方向为x轴，竖直方向为y 轴，以步骤8找到的各点到y轴的水平距离作为横坐标，纵坐标则 取每个点的高程值，可以将最高点的横坐标定为0，之后计算中每相 邻两点间的水平距离，作为直角坐标系中每两个相邻点的横坐标间 距，从而绘制出纵剖面曲线；

(10)利用polyfit函数将步骤9直角坐标系中的所有点拟合成 一条三次函数曲线，并通过反馈回来的三次函数的系数计算出极值 点坐标和拐点坐标；

(11)当曲线两端含有极值点时，需要删除掉极值点与端点之间 的部分，这是为了避免局部的凹凸变化影响整体坡形的识别。其中 “两端”的定义可以通过对参数赋值来进行控制，如取曲线前后各 四分之一为“两端”范围；

(12)判断曲线中部是否存在拐点，其中“中部”的定义同样可 以通过对参数赋值来进行控制，如取曲线四分之二到四分之三的范 围为“中部”。如果存在拐点，则分别判断拐点以上部分和拐点以 下部分的凹凸性，上半部分凹下半部分凸的为上凹下凸形坡，上半 部分凹下半部分平的为凹形坡，由于对称性不可能存在上半部分凹 下半部分凹，上半部分凸下半部分凹的为上凸下凹形坡，上半部分 凸下半部分平的为凸形坡，由于对称性不可能存在上半部分凸下半 部分凸，上半部分平下半部分凸的为凸形坡，上半部分平下半部分 凹的为凹形坡，上半部分平下半部分平的为平坡。如果不存在拐点 则直接判断整条曲线的凹凸性。凹凸性的判断是通过比较曲线两端 点纵坐标平均值与曲线中点纵坐标值得到的，比较的过程中需要指 定一个阈值，两个值的差别在阈值范围内的视为平坡；

(13)对每一种坡形赋相应的数值，经过多次循环后得到每一个 单元的坡形值，将这些坡形值利用属性表连接导回到GIS软件中， 绘制出最终的所测区域坡形分布图。

取某一山区村镇作为实施例，利用上述程序，通过所测区域的 DEM数据得出了每个单元内斜坡的最陡路径以及纵剖面线，图3与图 4为其中一个单元的结果，并得到了各单元内斜坡的坡形识别结果， 最终得到该实施例地区的坡形分布图，如图5所示。