CN105160658A - 一种基于子流域边界和流路特征的山脊线提取方法 - Google Patents

Abstract

本发明公开了一种基于子流域边界和流路特征的山脊线提取方法，以流域分水岭提取为基础，根据水流流路特征去除分水岭中的非山脊部分，以实现山脊线的准确提取。技术方案为：基于DEM进行填洼、流向计算和累计汇流量计算；设定汇流累积量阈值实现河道和子流域提取，获得子流域边界；基于流向和河流栅格数据计算每个栅格单元的入河流路长度；基于DEM和河流栅格数据计算每个栅格的流域相对高程HAND；提取入河流路长度和HAND值大于一定阈值的子流域边界，并进行去噪、细化和矢量化，最终得到山脊线。本发明突破了现有山脊线提取方法受DEM噪声影响大、提取的山脊线形状破碎等问题，在资源环境管理与地理国情普查等领域具有广阔的应用前景。

Description

一种基于子流域边界和流路特征的山脊线提取方法

技术领域

本发明涉及一种基于子流域边界和流路特征的山脊线提取方法，技术应用领域为资源环境管理与地理国情普查。

背景技术

由大自然长期作用所形成的山脊线在几何和物理两个方面均具有重要意义。在几何意义方面，它构成了地形起伏变化的分界线(骨架线)，因而在地形表达和DEM构建方面有重要作用；在物理意义方面，山脊线具有分水性，因而对地理过程分析和工程设计实施具有重要参考意义。

已有的山脊线提取方法可分为两类：基于局部地形几何特征的方法和基于整体地表流水模拟的方法。基于局部地形几何特征的方法包括等高线曲率最大判别法、等高线骨架化法和地形断面极值法。等高线曲率最大判别法和地形断面极值法均先提取山脊候选点，再根据一定的规则进行筛选得到山脊线。在提取山脊候选点时，等高线曲率最大判别法先计算等高线上各点的曲率值，然后找出局部曲率最大值点，之后由局部曲率最大值来得到山脊线。该方法的缺点在于由于地形噪声的存在,判定山脊线上的候选点时采用同一阈值,因而不能顾及区域地形变化的特殊规律,当阈值选择过大时,会遗漏一些山脊线和山谷线上的点；反之，则会增加候选点中的噪声。地形断面极值法通过计算地形断面上高程极值点的办法确定山脊候选点，由于该方法通常只采用两个正交方向上的地形断面，因此它会丢失某些方向的山脊点。等高线骨架化方法又称为中心轴化法，该方法先求取每条等高线弯曲部分的骨架，然后连接相邻等高线的骨架以得到山脊线，该方法将山脊线两侧的地形视为对称变化，这与多数地形的实际变化不相符合。

基于整体地表流水模拟的方法有三维地形表面流水数字模拟法和等高线垂线跟踪法。三维地形表面流水数字模拟法通过分析流水的运动情况提取山脊线，具体实现方法为：先根据流水方向计算每个点的汇流累积量，并根据汇流累积值及其变化规律确定合水线；然后，找出各合水线的汇水区域边界线(分水线)作为山脊线。等高线垂线跟踪法与三维地形表面流水数字模拟法基本原理相同，不同之处是后者通过汇流累积量确定合水线和分水线，而前者通过跟踪坡向的方式得到合水线和分水线。由于分水线与山脊线的概念不同，这类方法所确定的山脊线皆为闭合曲线，与实际山脊线存在较大差异。

综上所述，现有的山脊线提取方法受到算法原理、地形噪声等方面的影响，实用性受到很大限制。

发明内容

本发明要解决的技术问题是：克服现有山脊线提取方法受DEM噪声影响大、提取的山脊线形状破碎等问题，以流域分水岭提取为基础，提供一种基于子流域边界和流路特征的山脊线提取方法，根据水流流路特征去除分水岭中的非山脊部分，从而实现山脊线的准确提取。

本发明的技术解决方案为：基于DEM(DigitalElevationModel)进行填洼、流向计算、累计汇流量计算，并设定汇流累积量阈值实现河道和子流域提取，获得子流域边界；基于流向和河流计算每个栅格单元的入河流路长度；基于DEM和河流计算每个栅格的流域相对高程HAND(HeightAbovetheNearestDrainage)；提取入河流路长度和HAND大于一定阈值的子流域边界，并进行去噪、细化和矢量化，最终得到完整、连续的山脊线。具体步骤如下：

(1)采用栅格DEM(DigitalElevationModel)为输入数据，进行填洼、水流方向计算(D8单流向算法)，得到流向栅格数据；

(2)根据流向数据进行累计汇流量计算，并根据用户指定的阈值提取河流，得到河流栅格数据；

(3)根据河流和流向栅格数据，进行子流域划分得到子流域栅格图层，并提取子流域边界；从子流域栅格图层提取子流域边界的方法为：遍历子流域栅格，判断某一栅格是否与其周围所有栅格具有相同的子流域编码，如果该条件为True，则该栅格在子流域内部，如果为False，则该栅格位于两个子流域的边界；

(4)根据河流和流向栅格图层(二者的空间范围和分辨率相同)，计算入河流路长度，计算方法为：遍历每个栅格单元，执行如下操作：首先，沿流线向下游追踪，将流路中的每个栅格单元加入堆栈，满足如下情况之一时当前追踪结束：①遇到河道栅格；②遇到流路长度已知的栅格；③超出栅格图层的边界；然后，逐个弹出当前流路堆栈中的栅格单元，按如下方式计算入河流路长：①若为河流栅格，流路长为0；②若为栅格图层的边界，入河流路长为该栅格自身的流路长，即直角边分别为水平流路长和相邻栅格单元高程差的直角三角形斜边(流向与坐标轴平行时，水平流路长为栅格单元边长；流向与坐标轴呈45度角时，水平流路长为乘以栅格单元边长)；③对于坡面栅格单元，入河流路长为当前栅格单元的流路长加上其下游单元的入河流路长；

(5)以河流、流向和DEM栅格图层(三者的空间范围和分辨率相同)为输入，遍历每个栅格单元，寻找距其最近的河道单元，以当前栅格单元的高程减去距其最近的河道单元的高程，得到流域相对高程(HAND)；具体计算方法为：遍历每个栅格单元，执行如下操作：首先，沿流线向下流追踪，将流路中的每个栅格单元加入堆栈，满足如下情况之一时当前追踪结束：①遇到河道栅格，记下该河道栅格的高程；②遇到已知最近河道单元的栅格，记下其最近河道栅格的高程；③超出栅格图层的边界，以流路中边界栅格的高程作为该流路上最近河道栅格的高程；然后，逐个弹出当前流路堆栈中的栅格单元，以栅格单元自身高程减去距其最近的河道栅格的高程，得到HAND值；

(6)提取入河流路长度和HAND值大于一定阈值的分水岭，作为备选山脊线。入河流路长和HAND阈值可以用户指定，也可以根据规则计算默认值，其计算方法为：判断输入的HAND阈值是否为空，如果为空，则计算流域范围内HAND图层的平均值，并将其作为HAND阈值；判断输入的入河流路长度阈值是否为空，如果为空，则判断分水岭上所有栅格单元入河流路长度的平均值，并将其作为流路长度阈值；

(7)对备选山脊线进行滤波，以去除噪声。滤波方法为：遍历备选山脊线栅格图层，判断滤波窗口是否完全包含了某一块孤立的备选山脊线(滤波窗口的边界上没有山脊线栅格单元)，如果是，则去除该窗口内的山脊线；滤波窗口大小默认为15*15个栅格单元，用户可根据需要进行修改；

(8)对滤波后的山脊线进行细化和矢量化，得到最终的山脊线；细化和矢量化操作可以调用ArcGIS、GRASS等成熟GIS软件实现。

本发明的原理是：山脊线具有分水性，即：山脊线是分水岭的一部分，所以先提取分水岭，再去除其中不是山脊的部分，就能得到山脊线。分水岭提取可以通过先进行子流域划分，再提取子流域边界的方式进行；山脊线提取可根据山脊线上的流路特征进行，即：山脊线相对其最近河流的高度大于一定阈值，且从山脊线上的一点开始的流路长度也大于一定阈值，通过设定合适的阈值就可以实现从分水岭中提取山脊线的目的。

本发明与现有技术相比的优点在于：基于局部地形几何特征的方法所提取的离散地形特征点相互连接存在较大困难，甚至无法连接，而且此类方法受DEM噪声影响较大。而基于地表流水模拟的方法虽然提取结果整体性较好，但对格网DEM数据质量的要求较高，算法相对复杂，速度也较慢；此外，由于此类方法提取的闭合曲线与实际的山脊线有所差异，这会使山脊线提取的准确性受到影响。本发明克服现有山脊线提取方法受DEM噪声影响大、提取的山脊线形状破碎等问题，以流域分水岭提取为基础，根据流路特征去除分水岭中的非山脊部分，从而实现山脊线的准确提取。

附图说明

图1为本发明的实现流程图；

图2为入河流路长度计算流程图；

图3为相对高程HAND计算流程图；

图4为太湖流域丘陵区中田舍流域的山脊线提取结果。

具体实施方式

下面结合如图1所示的流程图，以太湖流域丘陵区中田舍流域(万km²)的山脊线提取为例，说明本发明的具体实施方法：

(1)选择太湖流域典型的低山丘陵小流域——中田舍流域为研究区，该流域位于江苏省溧阳市饮，流域干流长10.4km，面积约42km²。流域内存在两种不同类型的地貌，南部以山地为主，高程范围大致为180-530m，主要土壤类型为砾石土；北部以残丘为主，残丘的高程在100m左右，主要土壤类型为黄砂土。

(2)采用中田舍流域DEM(分辨率25m)为输入数据，利用TauDEM软件，进行填洼、水流方向计算(D8单流向算法)，得到流向栅格数据；并根据流向数据进行累计汇流量计算，以50个栅格为汇流阈值提取河流，得到河流栅格数据。

(3)根据河流和流向栅格数据，进行子流域划分得到子流域栅格图层，并通过Python语言编程实现子流域边界的提取，读写数据采用GDAL库。

(4)根据图2所示流程图，通过Python语言编程，根据河流和流向栅格图层(二者的空间范围和分辨率相同)，计算入河流路长度。

(5)根据图3所示流程图，以河流、流向和DEM栅格图层(三者的空间范围和分辨率相同)为输入，通过Python语言编程遍历每个栅格单元，寻找距其最近的河道单元，以当前栅格单元的高程减去距其最近的河道单元的高程，得到流域相对高程(HAND)。

(6)通过Python编程，提取入河流路长度和HAND值大于一定阈值的分水岭，作为备选山脊线；入河流路长度阈值设置为200m，HAND阈值采用程序自动计算的默认值。

(7)通过Python编程，对备选山脊线进行滤波，以去除噪声；滤波方法为：遍历备选山脊线栅格图层，判断滤波窗口是否完全包含了某一块孤立的备选山脊线(滤波窗口的边界上没有山脊线栅格单元)，如果是，则去除该窗口内的山脊线。滤波窗口大小采用默认值15*15个栅格单元。

(8)利用ArcGIS软件，对滤波后的山脊线进行细化和矢量化，得到山脊线。

山脊线的提取效果如图4所示，可以看出：与已有方法相比，使用本发明方法所提取的山脊线因为不受固定阈值和山脊点分布方向的限制，所以可以较全面的提取研究区范围内的所有山脊线，而且该方法还具有计算简便，准确性高等优点，在地形分析、流域模拟和土壤制图等方面具有广泛的应用前景。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (5)

1.一种基于子流域边界和流路特征的山脊线提取方法，其特征在于，具体实现步骤如下：

(1)采用栅格DEM为输入数据，利用D8单流向算法进行填洼和水流方向计算，得到流向栅格数据；

(2)根据流向栅格数据进行累计汇流量计算，并根据用户指定的阈值提取河流，得到河流栅格数据；

(3)根据河流和流向栅格数据，进行子流域划分得到子流域栅格图层，并提取子流域边界；从子流域栅格图层提取子流域边界的方法为：遍历子流域栅格，判断某一栅格是否与其周围所有栅格具有相同的子流域编码，如果该条件为True，则该栅格在子流域内部，如果为False，则该栅格位于两个子流域的边界；

(4)根据河流栅格图层和流向栅格图层，计算入河流路长度，河流栅格图层和流向栅格图层的空间范围和分辨率相同；

(5)以河流栅格图层、流向栅格图层和DEM栅格图层为输入，遍历每个栅格单元，寻找距其最近的河道单元，以当前栅格单元的高程减去距其最近的河道单元的高程，得到流域相对高程HAND，其中河流栅格图层、流向栅格图层和DEM栅格图层的空间范围和分辨率相同；

(6)提取入河流路长度和相对高程HAND值大于一定阈值的分水岭，作为备选山脊线；入河流路长和相对高程HAND阈值可以由用户指定，也可以根据规则计算默认值；

(7)对备选山脊线进行滤波，以去除噪声；

(8)对滤波后的山脊线进行细化和矢量化，得到最终的山脊线。

2.根据权利要求1所述的一种基于子流域边界和流路特征的山脊线提取方法，其特征在于，所述步骤(4)计算入河流路长度的方法为：遍历每个栅格单元，执行如下操作：首先，沿流线向下游追踪，将流路中的每个栅格单元加入堆栈，满足如下情况之一时当前追踪结束：①遇到河道栅格；②遇到流路长度已知的栅格；③超出栅格图层的边界；然后，逐个弹出当前流路堆栈中的栅格单元，按如下方式计算入河流路长：①若为河流栅格，流路长为0；②若为栅格图层的边界，入河流路长为该栅格自身的流路长，即直角边分别为水平流路长和相邻栅格单元高程差的直角三角形斜边，当流向与坐标轴平行时，水平流路长为栅格单元边长；当流向与坐标轴呈45度角时，水平流路长为乘以栅格单元边长；③对于坡面栅格单元，入河流路长为当前栅格单元的流路长加上其下游单元的入河流路长。

3.根据权利要求1所述的一种基于子流域边界和流路特征的山脊线提取方法，其特征在于，所述步骤(5)的相对高程HAND的具体计算方法为：遍历每个栅格单元，执行如下操作：首先，沿流线向下流追踪，将流路中的每个栅格单元加入堆栈，满足如下情况之一时当前追踪结束：①遇到河道栅格，记下该河道栅格的高程；②遇到已知最近河道单元的栅格，记下其最近河道栅格的高程；③超出栅格图层的边界，以流路中边界栅格的高程作为该流路上最近河道栅格的高程；然后，逐个弹出当前流路堆栈中的栅格单元，以栅格单元自身高程减去距其最近的河道栅格的高程，得到HAND值。

4.权利要求1所述的一种基于子流域边界和流路特征的山脊线提取方法，其特征在于，所述步骤(6)的具体计算方法为：判断输入的相对高程HAND阈值是否为空，如果为空，则计算流域范围内HAND图层的平均值，并将其作为相对高程HAND阈值；判断输入的入河流路长度阈值是否为空，如果为空，则判断分水岭上所有栅格单元入河流路长度的平均值，并将其作为流路长度阈值。

5.权利要求1所述的一种基于子流域边界和流路特征的山脊线提取方法，其特征在于，所述步骤(7)的滤波方法为：遍历备选山脊线栅格图层，判断滤波窗口是否完全包含了某一块孤立的备选山脊线，其中滤波窗口的边界上没有山脊线栅格单元，如果是，则去除该窗口内的山脊线。