CN113379828B - 一种融合地表形态特征的坡长提取方法 - Google Patents

Abstract

本发明提供了一种融合地表形态特征的坡长提取方法，包括：步骤一，DEM数据预处理：步骤二，根据正DEM计算D8流向和等高线曲率：步骤三，根据负DEM计算up_iGD8上坡向：步骤四，根据D8计算汇水面积并设置阈值确定沟道截断矩阵：步骤五，根据up_iGD8上坡向计算坡度并确定坡度截断矩阵：步骤六，根据up_iGD8上坡向计算上坡坡度线长度：步骤七，根据坡长计算公式计算得到坡长。坡长计算公式为：

本发明的方法在连续的解空间中推导了坡长理论模型，可以得到一种直接的坡长计算方法，避免了对上坡汇水面积和有效等高线宽度的估计来提取坡长，有效降低了计算原理中潜在的误差风险。

Description

一种融合地表形态特征的坡长提取方法

技术领域

本发明属于数字地形分析领域，涉及坡长提取，具体涉及一种融合地表形态特征的坡长提取方法。

背景技术

通用土壤流失方程(USLE)和修正的通用土壤流失方程(RUSLE)因其高效性和健壮性被广泛应用于土壤侵蚀预测和评价。坡长因子(L-factor)是USLE/RUSLE模型中的一个重要地形因子，L-factor是坡长(Slope Length)除以标准长度(22.13m)后取坡度指数的一个比例参数，反映了一种水流携带泥沙的侵蚀强度。其中，坡长的基本定义是从径流源点到坡面泥沙发生沉积或到河网截断处的水平投影距离。

然而坡长因用户的不同理解而存在诸多的计算方法，在区域尺度中一般基于数字高程模型(DEM)进行提取，广泛研究的坡长提取包括：基于一维径流路径的直接计算方法，基于单位汇水面积的代替坡长计算法。然而，直接计算法通过计算流路的累计长度导致不能反映坡面的收敛、发散的变化，其结果与地表的形态变化吻合程度较差；代替计算法主要是通过估算上坡水流的汇水面积和所在出水口的有效等高线宽度，将两者做比按照单位汇水面积来代替坡长进行计算，这种方法在估算上坡汇水面积和有效等高线宽度时的误差较大，严重影响坡长的提取结果。前人研究中多集中于用不同方法计算上坡汇水面积来尝试提取到更合理的坡长，但上坡汇水面积和有效等高线宽度的两因素误差未曾得到有效的解决。

发明内容

针对现有技术存在的不足，本发明的目的在于，提供一种融合地表形态特征的坡长提取方法，解决现有技术中的坡长提取方法的误差较大的技术问题。

为了解决上述技术问题，本发明采用如下技术方案予以实现：

一种融合地表形态特征的坡长提取方法，该方法采用的坡长计算公式为：

式中：

λ表示坡长；

λ₀表示局部最高点的坡长；

l表示坡度线长度；

k_c表示等高线曲率。

本发明还具有如下技术特征：

该方法在DEM中提取坡长，具体包括以下步骤：

步骤一，DEM数据预处理：

步骤1.1，读取DEM数据：

Step101，创建demMap对象，读取DEM数据的头信息；

所述的头信息为栅格尺寸cellsize、栅格总行数rows和栅格总列数cols；

Step102，读取栅格矩阵并记录到浮点型的二维数组demData中；

步骤1.2，DEM数据缺失修复与填洼处理：

采用单环搜索方法对DEM数据进行洼地填充，遍历二维数组demData，在3×3的搜索窗口中，每次判断中心栅格数值是否均小于周围8个，如果满足该条件就将中心栅格赋值为8邻域中的最小值；如果中心栅格缺失数据，则将中心栅格赋值为8邻域中的最小值；

步骤二，根据正DEM计算D8流向和等高线曲率：

步骤2.1，计算D8的流向：

经过步骤一的DEM数据预处理后的DEM即为正DEM；

遍历二维数组demData，在3×3的搜索窗口中，每次计算中间栅格与8邻域栅格的带权高程差，即中间栅格数值减去邻域栅格数值，做差后除以两个栅格的距离，其中坐标系方向距离为1，对角线方向距离为

取带权高程差最大的邻域栅格的方向为中心栅格的流向；

D8的可能的方向包括：东、东南、南、西南、西、西北、北和东北这8个方向，分别记录方向代号为1、2、4、8、16、32、64和128；

步骤2.2，计算等高线曲率：

等高线曲率的计算方法采用三阶不带权差分方法计算，其计算公式为：

式中：

k_c表示等高线曲率；

d表示栅格尺寸cellsize

x和y分别表示形态学连续解空间中的地表曲面在二维平面上投影的横坐标和纵坐标；

z表示形态学连续解空间中的高程；

p表示z对x的一阶偏导；q表示z对y的一阶偏导；r表示z对x的二阶偏导；s表示z对xy的二阶偏导；t表示z对y的二阶偏导；

步骤三，根据负DEM计算up_iGD8上坡向：

负DEM数据是将预处理后的正DEM进行数值的反操作得到，即让二维数组demData的每个元素减去全部数据中的最小值，就可得到负DEM二维数组；

在正DEM中计算每个栅格最可能的上坡方向如同在负DEM中计算最陡流出方向，即原DEM的up_iGD8上坡向，具体的计算过程与步骤2.1相同；

步骤四，根据D8计算汇水面积并设置阈值确定沟道截断矩阵：

在步骤2.1中得到了每个栅格的流向，根据此流向进行河网提取，所述的河网提取采用流域拓扑方法来计算，通过计算得到的汇水面积，设置一个河网阈值，即当汇水面积矩阵中的元素大于规定阈值后标记该元素的位置为沟道截断位置；依次将满足条件的元素保留在新的矩阵中，其数值记为1，其它元素记为0，最终生成沟道截断矩阵；

步骤五，根据up_iGD8上坡向计算坡度并确定坡度截断矩阵：

在步骤三中得到了每个栅格的up_iGD8上坡向，坡度是根据坡向的进一步计算，选择中心栅格的上坡向所指的邻域栅格，计算两者高程差Δh和水平距离Δd，然后取反正切，记为式Ⅵ：

slope＝180/π·arctan(Δh/Δd) 式VI；

当坡向为坐标系方向时，水平距离Δd为1倍的栅格尺寸；当坡向为对角线方向时，水平距离Δd为

倍的栅格尺寸；坡度截断的标准是：

如果当前栅格的坡度大于2.86°时，坡度变化率超过50％发生截断；

如果当前栅格的坡度小于等于2.86°时，坡度变化率超过70％发生截断；

将满足上述条件的栅格的位置记录下来，保留在新的矩阵中，其数值记为1，其它元素记为0，最终生成坡度截断矩阵；

步骤六，根据up_iGD8上坡向计算上坡坡度线长度：

采用局部坡向合并策略计算上坡坡度线长度；

所述的局部坡向合并策略的具体步骤如下：

Step601，从每个栅格开始建立辅助队列，按照上坡向所指的方向依次将邻域栅格加入队列，并将Δd进行累加；

Step602，根据Step601的方法，每累加3个步长，记录初始栅格和第四个栅格的所在行列位置，用这两个栅格的欧氏距离代替之前的3次累计步长；

Step603，按照Step601和Step602全部计算完整个坡度线长度，一般在DEM的局部最高点处停止计算，局部最高点的特征是其上坡向不存在；

Step604，按照Step601至Step603的方法遍历完二维数组demData的全部数据；

步骤七，根据坡长计算公式计算得到坡长：

Step701，确定DEM数据中的局部最高点，即up_iGD8上坡向不存在的栅格位置，选择一个局部最高点建立一个搜索队列，并将其作为第一个元素放入队列，初始化当前栅格的上坡截断状态，即CUT＝0；

Step702，采用层序搜索的方法进行当前栅格的8邻域遍历，根据up_iGD8上坡向确定能指向当前栅格的邻域，使其入队列，并按照式Ⅶ进行无截断的坡长计算；

Step703，对队列的首元素判断其位置是否和截断矩阵中的标记元素对应，如果队首元素在截断位置上，则比较当前栅格无截断的坡长值和上坡截断状态CUT，取两者的最大值以更新CUT；如果队首元素不在截断位置上，则不更新截断状态CUT的值；

Step704，用无截断的坡长减去CUT值以更新邻域栅格的真实坡长值；

Step705，将邻域的截断状态变量继承为当前栅格的上坡截断状态CUT值；

Step706，按照Step702至Step705，依次将邻域栅格能指向队首元素的情况遍历完，然后将当前栅格移除，更新队列的首元素；

Step707，按照Step702至Step706，依次处理完DEM中的一个局部最高点的搜索；

Step708，更换另一个局部最高点，按照Step701至Step707完成全局数据的搜索；

Step709，结束全局数据的搜索，导出截断后的坡长结果矩阵。

具体的，步骤四中，所述的流域拓扑方法具体步骤如下：

Step401，建立一个辅助队列，用来记录径流路径；

Step402，遍历二维数组demData，将每个栅格入队列，设置每个栅格的初始流量为1；

Step403，然后沿着当前栅格的流向向邻域的栅格延伸该路径，使邻域栅格依次入队，当前栅格依次出队，每次给邻域栅格的流量累加当前栅格的流量；

Step404，按照Step402和Step403遍历完全部数据后得到一个汇水面积矩阵，记录了每个栅格元素的累计流量。

本发明与现有技术相比，具有如下技术效果：

(Ⅰ)本发明的方法在连续的解空间中推导了坡长理论模型，可以得到一种直接的坡长计算方法，避免了对上坡汇水面积和有效等高线宽度的估计来提取坡长，有效降低了计算原理中潜在的误差风险。

(Ⅱ)本发明的方法根据理论模型的解可在DEM中实现，完整的将各个计算部分集成以便于GIS用户使用。从最终结果来看：提取到的坡长空间分布与地表形态结构具有较高的吻合程度，坡长提取结果可接受。

(Ⅲ)本发明的方法综合考虑了坡面二维形态学流域面积和坡长的截断条件，充分满足USLE/RUSLE中的定义要求，提取的坡长作为相关输入参数能够直接和USLE/RUSLE模型耦合。

(Ⅳ)与现有方法相比，本发明能有效降低坡长提取结果的误差，尤其是在水流强发散的坡面和截断位置，本发明提取的坡长分布结果能更好的与地表形态相吻合。本发明是基于形态学流域面积提取坡面地形因子的有益尝试，也为利用通用土壤流失方程进行土壤侵蚀评价与预测提供了技术支持。

附图说明

图1是坡长理论模型的示意图。

图2是基于格网DEM的坡长提取流程图。

图3是等高线曲率计算方法以及对应栅格位置图。

图4是上坡坡度线的计算方法示意图。

图5是地形图和提取的坡长的空间分布图。

图6是提取的坡长的累计频率曲线。

以下结合实施例对本发明的具体内容作进一步详细解释说明。

具体实施方式

基于现有技术中的情况，如何避免估算上坡汇水面积和有效等高线宽度的两因素计算误差，从而提取二维坡面形态学的坡长的是本申请的研究创新。基于连续解空间中坡长理论模型的解，结合坡面单位汇水面积的微分方程，融合坡度线、等高线曲率等地表形态特征，有效的避免了对上坡汇水面积和有效等高线宽度的两因素估计误差，能够准确地对坡长进行计算，对于区域尺度中提取坡长应用于土壤侵蚀预测起到重要的促进作用。

本发明基于坡面的二维形态学流域面积原理，综合考虑了坡面水流的二维过程和坡长截断的条件，在形态学连续解空间中建立坡长理论模型。根据模型的解设计并实现的一种融合地形表面特征(坡度线，等高线曲率，径流截断因子)的坡长提取方法。

所述的坡长理论模型的构建过程为：

坡长理论模型的基本假设是：地形表面为纯形态学曲面，曲面一阶导数和二阶导数连续，不考虑地表溶质的下渗等复杂过程。

若用一个隐式方程式Ⅰ表示地表曲面，沿着垂直于等高线方向的线积分过程就是坡度线的计算方法。

z＝f(x，y) 式I；

式中：

z表示形态学连续解空间中的高程；

x和y分别表示形态学连续解空间中的地表曲面在二维平面上投影的横坐标和纵坐标。

在形态学连续解空间中，定义曲面上有起始点P₁(x₁,y₁)和结束点P₂(x₂,y₂)，则两点间的坡度线长度l为式Ⅱ。

式中：

l表示坡度线长度；

f_x表示z对x的一阶偏导；

f_y表示z对y的一阶偏导。

根据坡面上二维形态学流域的概念，结合坡长和L-factor的定义，坡长是一定坡面的投影面积对坡宽的比值，在连续解空间中，坡长λ的计算式如式Ⅲ所示：

式中：

λ表示坡长；

A表示上坡汇水面积；

w表示有效等高线宽度(即，坡宽)。

通过分析推导式Ⅱ和式Ⅲ的非线性关系，可以得到一个λ对l的偏微分方程如式Ⅳ所示；

式中：

k_c表示等高线曲率。

通过式Ⅳ的计算，能有效避免估计上坡贡献面积和有效等高线宽度。

上述过程主要考虑了坡长的中间过程，对于定义中的坡长截断条件的判断，主要分为两种情况，即：

坡度变化在一定条件下形成的坡面泥沙淤积饱和的情况称之为坡度截断；

径流汇集形成河道对坡长的截断称之为沟道截断。

图1是坡长理论模型的示意图，只要能确定坡长的两种截断方式得到截断点的位置，根据式式Ⅱ和式Ⅳ就可以对坡长进行提取。

结合以上坡长理论模型，可在格网DEM数据中对其进行离散化求解。DEM是将连续的地表离散化为规则正方形的栅格，栅格尺寸为DEM分辨率，栅格的行数×列数为DEM数据量。

以下给出本发明的具体实施例，需要说明的是本发明并不局限于以下具体实施例，凡在本申请技术方案基础上做的等同变换均落入本发明的保护范围。

实施例1：

遵从上述技术方案，如图2所示，本实施例给出一种融合地表形态特征的坡长提取方法，该方法在DEM中提取坡长，具体包括以下步骤：

步骤一，DEM数据预处理：

步骤1.1，读取DEM数据：

Step101，创建demMap对象，读取DEM数据的头信息；

所述的头信息为栅格尺寸cellsize、栅格总行数rows和栅格总列数cols；

Step102，读取栅格矩阵并记录到浮点型的二维数组demData中；

步骤1.2，DEM数据缺失修复与填洼处理：

针对原始DEM数据中可能存在的数据问题进行预处理，保证后续计算的有效性。

采用单环搜索方法对DEM数据进行洼地填充，遍历二维数组demData，在3×3的搜索窗口中，每次判断中心栅格数值是否均小于周围8个，如果满足该条件就将中心栅格赋值为8邻域中的最小值。如果中心栅格缺失数据，则将中心栅格赋值为8邻域中的最小值；

步骤二，根据正DEM计算D8流向和等高线曲率：

步骤2.1，计算D8的流向：

D8流向的原理是找到每个栅格最可能流向8个邻域的方向。具体方法是：

经过步骤一的DEM数据预处理后的DEM即为正DEM；

遍历二维数组demData，在3×3的搜索窗口中，每次计算中间栅格与8邻域栅格的带权高程差，即中间栅格数值减去邻域栅格数值，做差后除以两个栅格的距离，其中坐标系方向距离为1，对角线方向距离为

取带权高程差最大的邻域栅格的方向为中心栅格的流向；

D8的可能的方向包括：东、东南、南、西南、西、西北、北和东北这8个方向，分别记录方向代号为1、2、4、8、16、32、64和128；

步骤2.2，计算等高线曲率：

等高线曲率的计算方法采用三阶不带权差分方法计算，其计算公式为：

式中：

k_c表示等高线曲率；

d表示栅格尺寸cellsize

x和y分别表示形态学连续解空间中的地表曲面在二维平面上投影的横坐标和纵坐标；

z表示形态学连续解空间中的高程；

p表示z对x的一阶偏导；q表示z对y的一阶偏导；r表示z对x的二阶偏导；s表示z对xy的二阶偏导；t表示z对y的二阶偏导；

式Ⅴ中一阶偏导和二阶偏导对应的栅格的位置如图3所示。

步骤三，根据负DEM计算up_iGD8上坡向：

负DEM数据是将预处理后的正DEM进行数值的反操作得到，即让二维数组demData的每个元素减去全部数据中的最小值，就可得到负DEM二维数组；此操作发目的是计算每个栅格的最陡上坡方向，以得到目标中的上坡坡度线长度。

在正DEM中计算每个栅格最可能的上坡方向如同在负DEM中计算最陡流出方向，即原DEM的up_iGD8上坡向，具体的计算过程与步骤2.1相同；

步骤四，根据D8计算汇水面积并设置阈值确定沟道截断矩阵：

在步骤2.1中得到了每个栅格的流向，根据此流向进行河网提取，所述的河网提取采用流域拓扑方法来计算，通过计算得到的汇水面积，设置一个河网阈值，即当汇水面积矩阵中的元素大于规定阈值后标记该元素的位置为沟道截断位置；依次将满足条件的元素保留在新的矩阵中，其数值记为1，其它元素记为0，最终生成沟道截断矩阵；

步骤四中，所述的流域拓扑方法具体步骤如下：

Step401，建立一个辅助队列，用来记录径流路径；

Step402，遍历二维数组demData，将每个栅格入队列，设置每个栅格的初始流量为1；

Step403，然后沿着当前栅格的流向向邻域的栅格延伸该路径，使邻域栅格依次入队，当前栅格依次出队，每次给邻域栅格的流量累加当前栅格的流量；

Step404，按照Step402和Step403遍历完全部数据后得到一个汇水面积矩阵，记录了每个栅格元素的累计流量；

步骤五，根据up_iGD8上坡向计算坡度并确定坡度截断矩阵：

在步骤三中得到了每个栅格的up_iGD8上坡向，坡度是根据坡向的进一步计算，选择中心栅格的上坡向所指的邻域栅格，计算两者高程差Δh和水平距离Δd，然后取反正切，记为式Ⅵ：

slope＝180/π·arctan(Δh/Δd) 式VI；

当坡向为坐标系方向时，水平距离Δd为1倍的栅格尺寸；当坡向为对角线方向时，水平距离Δd为

倍的栅格尺寸；坡度截断的标准是：

如果当前栅格的坡度大于2.86°时，坡度变化率超过50％发生截断；

如果当前栅格的坡度小于等于2.86°时，坡度变化率超过70％发生截断；

将满足上述条件的栅格的位置记录下来，保留在新的矩阵中，其数值记为1，其它元素记为0，最终生成坡度截断矩阵；

步骤六，根据up_iGD8上坡向计算上坡坡度线长度：

在步骤三中得到了每个栅格的up_iGD8上坡向，沿着上坡向进行累加能够得到上坡坡度线的累计长度。但是由于上坡向的方向在离散的规则格网DEM中只有8个方向，按照局部上坡向进行累加得到的有锯齿状的曲折，其累计长度必定造成不确定的误差。因此本方法在此处设计了一个局部坡向合并策略来减少坡度线长度的误差，示意图如图4所示。

采用局部坡向合并策略计算上坡坡度线长度；

所述的局部坡向合并策略的具体步骤如下：

Step601，从每个栅格开始建立辅助队列，按照上坡向所指的方向依次将邻域栅格加入队列，并将Δd进行累加；

Step602，根据Step601的方法，每累加3个步长，记录初始栅格和第四个栅格的所在行列位置，用这两个栅格的欧氏距离代替之前的3次累计步长；

如图4中，局部坡向AB、BC和CD被合并为AD，对应的累计长度

被替换为

Step603，按照Step601和Step602全部计算完整个坡度线长度，一般在DEM的局部最高点处停止计算，局部最高点的特征是其上坡向不存在；

如图4中，整个上坡坡度线路径的累计长度是

Step604，按照Step601至Step603的方法遍历完二维数组demData的全部数据；

进行局部坡向合并的策略，将部分线段合并，有效减少曲折的折线的累计长度，而不影响直线路径的累计长度，以减少坡度线长度的误差。

步骤七，根据坡长计算公式计算得到坡长：

坡长λ对坡度线l的偏微分方程式Ⅳ是一种非线性关系，反映了随地表等高线曲率变化下的平均坡长的结果，根据坡长理论模型的解进行数值近似计算得到坡长。

坡长理论模型的近似解，该方法采用的坡长计算公式即如下：

式中：

λ表示坡长；

λ₀表示局部最高点的坡长；

l表示坡度线长度；

k_c表示等高线曲率；

结合式7和前面步骤所求的中间结果，设计一个完整的迭代计算方法来提取坡长，包含以下步骤：

Step701，确定DEM数据中的局部最高点，即up_iGD8上坡向不存在的栅格位置，选择一个局部最高点建立一个搜索队列，并将其作为第一个元素放入队列，初始化当前栅格的上坡截断状态，即CUT＝0；

Step702，采用层序搜索的方法进行当前栅格的8邻域遍历，根据up_iGD8上坡向确定能指向当前栅格的邻域，使其入队列，并按照式Ⅶ进行无截断的坡长计算；

Step703，对队列的首元素判断其位置是否和截断矩阵中的标记元素对应，如果队首元素在截断位置上，则比较当前栅格无截断的坡长值和上坡截断状态CUT，取两者的最大值以更新CUT；如果队首元素不在截断位置上，则不更新截断状态CUT的值；

Step704，用无截断的坡长减去CUT值以更新邻域栅格的真实坡长值；

Step705，将邻域的截断状态变量继承为当前栅格的上坡截断状态CUT值；

Step706，按照Step702至Step705，依次将邻域栅格能指向队首元素的情况遍历完，然后将当前栅格移除，更新队列的首元素；

Step707，按照Step702至Step706，依次处理完DEM中的一个局部最高点的搜索；

Step708，更换另一个局部最高点，按照Step701至Step707完成全局数据的搜索。

Step709，结束全局数据的搜索，导出截断后的坡长结果矩阵。

应用例：

将上述实施例中的融合地表形态特征的坡长提取方法在一个真实地形数据中进行了实验。

(1)实验数据：

陕西省榆林市县南沟流域。DEM数据分辨率2.5m，海拔高度1014～1434m。

(2)实验方案：

根据实施例中的具体过程，按照图2所示的流程进行基于DEM的坡长提取。结合实际的地形数据，其中坡长的沟道截断的汇水面积阈值设置为2000m²，坡度截断的阈值为：坡度≤2.86°，发生截断的坡度变化阈值为70％；坡度＞2.86°，发生截断的坡度变化阈值为50％。

(3)实验结果与分析：

使用本发明方法提取的坡长的空间分布如图5所示，对比地形的细节，可以发现提取的坡长结果和地形吻合度较高，充分反映了表面凹凸的变化。

此外，实验中将本发明与现有方向进行了比较。常见的四种基于流向算法使用单位汇水面积代替计算方法来提取坡长，这些提取方法主要侧重于估算上坡汇水面积，进而除以每个栅格的有效等高线宽度，往往存在两因素的估算误差，并且不同的流向算法得到的坡长空间分布也存在巨大差异。而本发明不再延续此类方法，而是根据偏微分方程的解直接进行单位汇水面积的计算进而提取坡长。

对于结果的数值统计来看(如图6所示)，本发明95％的坡长值小于189m，在强发散区域，坡长的提取结果明显小于其它方法，和坡面的侵蚀原理相一致，在坡长的截断位置，相关截断条件的集成使得坡长的值得到有效控制。相对于其它方法，新的发明的方法提取的坡长的数据统计结果也更为合理。

Claims (2)

1.一种融合地表形态特征的坡长提取方法，其特征在于，该方法采用的坡长计算公式为：

式中：

λ表示坡长；

λ₀表示局部最高点的坡长；

l表示坡度线长度；

k_c表示等高线曲率；

该方法在DEM中提取坡长，具体包括以下步骤：

步骤一，DEM数据预处理：

步骤1.1，读取DEM数据：

Step101，创建demMap对象，读取DEM数据的头信息；

所述的头信息为栅格尺寸cellsize、栅格总行数rows和栅格总列数cols；

Step102，读取栅格矩阵并记录到浮点型的二维数组demData中；

步骤1.2，DEM数据缺失修复与填洼处理：

采用单环搜索方法对DEM数据进行洼地填充，遍历二维数组demData，在3×3的搜索窗口中，每次判断中心栅格数值是否均小于周围8个，如果满足该条件就将中心栅格赋值为8邻域中的最小值；如果中心栅格缺失数据，则将中心栅格赋值为8邻域中的最小值；

步骤二，根据正DEM计算D8流向和等高线曲率：

步骤2.1，计算D8的流向：

经过步骤一的DEM数据预处理后的DEM即为正DEM；

遍历二维数组demData，在3×3的搜索窗口中，每次计算中间栅格与8邻域栅格的带权高程差，即中间栅格数值减去邻域栅格数值，做差后除以两个栅格的距离，其中坐标系方向距离为1，对角线方向距离为

取带权高程差最大的邻域栅格的方向为中心栅格的流向；

D8的可能的方向包括：东、东南、南、西南、西、西北、北和东北这8个方向，分别记录方向代号为1、2、4、8、16、32、64和128；

步骤2.2，计算等高线曲率：

等高线曲率的计算方法采用三阶不带权差分方法计算，其计算公式为：

式中：

k_c表示等高线曲率；

d表示栅格尺寸cellsize

x和y分别表示形态学连续解空间中的地表曲面在二维平面上投影的横坐标和纵坐标；

z表示形态学连续解空间中的高程；

p表示z对x的一阶偏导；q表示z对y的一阶偏导；r表示z对x的二阶偏导；s表示z对xy的二阶偏导；t表示z对y的二阶偏导；

步骤三，根据负DEM计算up_iGD8上坡向：

负DEM数据是将预处理后的正DEM进行数值的反操作得到，即让二维数组demData的每个元素减去全部数据中的最小值，就可得到负DEM二维数组；

在正DEM中计算每个栅格最可能的上坡方向如同在负DEM中计算最陡流出方向，即原DEM的up_iGD8上坡向，具体的计算过程与步骤2.1相同；

步骤四，根据D8计算汇水面积并设置阈值确定沟道截断矩阵：

在步骤2.1中得到了每个栅格的流向，根据此流向进行河网提取，所述的河网提取采用流域拓扑方法来计算，通过计算得到的汇水面积，设置一个河网阈值，即当汇水面积矩阵中的元素大于规定阈值后标记该元素的位置为沟道截断位置；依次将满足条件的元素保留在新的矩阵中，其数值记为1，其它元素记为0，最终生成沟道截断矩阵；

步骤五，根据up_iGD8上坡向计算坡度并确定坡度截断矩阵：

在步骤三中得到了每个栅格的up_iGD8上坡向，坡度是根据坡向的进一步计算，选择中心栅格的上坡向所指的邻域栅格，计算两者高程差Δh和水平距离Δd，然后取反正切，记为式Ⅵ：

slope＝180/π·arctan(Δh/Δd) 式VI；

当坡向为坐标系方向时，水平距离Δd为1倍的栅格尺寸；当坡向为对角线方向时，水平距离Δd为

倍的栅格尺寸；坡度截断的标准是：

如果当前栅格的坡度大于2.86°时，坡度变化率超过50％发生截断；

如果当前栅格的坡度小于等于2.86°时，坡度变化率超过70％发生截断；

将满足上述条件的栅格的位置记录下来，保留在新的矩阵中，其数值记为1，其它元素记为0，最终生成坡度截断矩阵；

步骤六，根据up_iGD8上坡向计算上坡坡度线长度：

采用局部坡向合并策略计算上坡坡度线长度；

所述的局部坡向合并策略的具体步骤如下：

Step601，从每个栅格开始建立辅助队列，按照上坡向所指的方向依次将邻域栅格加入队列，并将Δd进行累加；

Step602，根据Step601的方法，每累加3个步长，记录初始栅格和第四个栅格的所在行列位置，用这两个栅格的欧氏距离代替之前的3次累计步长；

Step603，按照Step601和Step602全部计算完整个坡度线长度，一般在DEM的局部最高点处停止计算，局部最高点的特征是其上坡向不存在；

Step604，按照Step601至Step603的方法遍历完二维数组demData的全部数据；

步骤七，根据坡长计算公式计算得到坡长：

Step701，确定DEM数据中的局部最高点，即up_iGD8上坡向不存在的栅格位置，选择一个局部最高点建立一个搜索队列，并将其作为第一个元素放入队列，初始化当前栅格的上坡截断状态，即CUT＝0；

Step702，采用层序搜索的方法进行当前栅格的8邻域遍历，根据up_iGD8上坡向确定能指向当前栅格的邻域，使其入队列，并按照式Ⅶ进行无截断的坡长计算；

Step703，对队列的首元素判断其位置是否和截断矩阵中的标记元素对应，如果队首元素在截断位置上，则比较当前栅格无截断的坡长值和上坡截断状态CUT，取两者的最大值以更新CUT；如果队首元素不在截断位置上，则不更新截断状态CUT的值；

Step704，用无截断的坡长减去CUT值以更新邻域栅格的真实坡长值；

Step705，将邻域的截断状态变量继承为当前栅格的上坡截断状态CUT值；

Step706，按照Step702至Step705，依次将邻域栅格能指向队首元素的情况遍历完，然后将当前栅格移除，更新队列的首元素；

Step707，按照Step702至Step706，依次处理完DEM中的一个局部最高点的搜索；

Step708，更换另一个局部最高点，按照Step701至Step707完成全局数据的搜索；

Step709，结束全局数据的搜索，导出截断后的坡长结果矩阵。

2.如权利要求1所述的融合地表形态特征的坡长提取方法，其特征在于，步骤四中，所述的流域拓扑方法具体步骤如下：

Step401，建立一个辅助队列，用来记录径流路径；

Step402，遍历二维数组demData，将每个栅格入队列，设置每个栅格的初始流量为1；

Step403，然后沿着当前栅格的流向向邻域的栅格延伸该路径，使邻域栅格依次入队，当前栅格依次出队，每次给邻域栅格的流量累加当前栅格的流量；

Step404，按照Step402和Step403遍历完全部数据后得到一个汇水面积矩阵，记录了每个栅格元素的累计流量。

Images