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

本发明专利技术提供了一种融合地表形态特征的坡长提取方法，包括：步骤一，DEM数据预处理：步骤二，根据正DEM计算D8流向和等高线曲率,步骤三，根据负DEM计算up_iGD8上坡向：步骤四，根据D8计算汇水面积并设置阈值确定沟道截断矩阵：步骤五，根据up_iGD8上坡向计算坡度并确定坡度截断矩阵：步骤六，根据up_iGD8上坡向计算上坡坡度线长度：步骤七，根据坡长计算公式计算得到坡长。坡长计算公式为：本发明专利技术的方法在连续的解空间中推导了坡长理论模型，可以得到一种直接的坡长计算方法，避免了对上坡汇水面积和有效等高线宽度的估计来提取坡长，有效降低了计算原理中潜在的误差风险。低了计算原理中潜在的误差风险。低了计算原理中潜在的误差风险。

【技术实现步骤摘要】
一种融合地表形态特征的坡长提取方法


[0001]本专利技术属于数字地形分析领域，涉及坡长提取，具体涉及一种融合地表形态特征的坡长提取方法。

技术介绍

[0002]通用土壤流失方程(USLE)和修正的通用土壤流失方程(RUSLE)因其高效性和健壮性被广泛应用于土壤侵蚀预测和评价。坡长因子(L
‑
factor)是USLE/RUSLE模型中的一个重要地形因子，L
‑
factor是坡长(Slope Length)除以标准长度(22.13m)后取坡度指数的一个比例参数，反映了一种水流携带泥沙的侵蚀强度。其中，坡长的基本定义是从径流源点到坡面泥沙发生沉积或到河网截断处的水平投影距离。
[0003]然而坡长因用户的不同理解而存在诸多的计算方法，在区域尺度中一般基于数字高程模型(DEM)进行提取，广泛研究的坡长提取包括：基于一维径流路径的直接计算方法，基于单位汇水面积的代替坡长计算法。然而，直接计算法通过计算流路的累计长度导致不能反映坡面的收敛、发散的变化，其结果与地表的形态变化吻合程度较差；代替计算法主要是通过估算上坡水流的汇水面积和所在出水口的有效等高线宽度，将两者做比按照单位汇水面积来代替坡长进行计算，这种方法在估算上坡汇水面积和有效等高线宽度时的误差较大，严重影响坡长的提取结果。前人研究中多集中于用不同方法计算上坡汇水面积来尝试提取到更合理的坡长，但上坡汇水面积和有效等高线宽度的两因素误差未曾得到有效的解决。
[0004]针对现有技术存在的不足，本专利技术的目的在于，提供一种融合地表形态特征的坡长提取方法，解决现有技术中的坡长提取方法的误差较大的技术问题。
[0005]为了解决上述技术问题，本专利技术采用如下技术方案予以实现：
[0006]一种融合地表形态特征的坡长提取方法，该方法采用的坡长计算公式为：
[0007][0008]式中：
[0009]λ表示坡长；
[0010]λ0表示局部最高点的坡长；
[0011]l表示坡度线长度；
[0012]k
c
表示等高线曲率。
[0013]本专利技术还具有如下技术特征：
[0014]该方法在DEM中提取坡长，具体包括以下步骤：
[0015]步骤一，DEM数据预处理：
[0016]步骤1.1，读取DEM数据：
[0017]Step101，创建demMap对象，读取DEM数据的头信息；
[0018]所述的头信息为栅格尺寸cellsize、栅格总行数rows和栅格总列数cols；
[0019]Step102，读取栅格矩阵并记录到浮点型的二维数组demData中；
[0020]步骤1.2，DEM数据缺失修复与填洼处理：
[0021]采用单环搜索方法对DEM数据进行洼地填充，遍历二维数组demData，在3
×
3的搜索窗口中，每次判断中心栅格数值是否均小于周围8个，如果满足该条件就将中心栅格赋值为8邻域中的最小值；如果中心栅格缺失数据，则将中心栅格赋值为8邻域中的最小值；
[0022]步骤二，根据正DEM计算D8流向和等高线曲率,：
[0023]步骤2.1，计算D8的流向：
[0024]经过步骤一的DEM数据预处理后的DEM即为正DEM；
[0025]遍历二维数组demData，在3
×
3的搜索窗口中，每次计算中间栅格与8邻域栅格的带权高程差，即中间栅格数值减去邻域栅格数值，做差后除以两个栅格的距离，其中坐标系方向距离为1，对角线方向距离为取带权高程差最大的邻域栅格的方向为中心栅格的流向；
[0026]D8的可能的方向包括：东、东南、南、西南、西、西北、北和东北这8个方向，分别记录方向代号为1、2、4、8、16、32、64和128；
[0027]步骤2.2，计算等高线曲率：
[0028]等高线曲率的计算方法采用三阶不带权差分方法计算，其计算公式为：
[0029][0030]式中：
[0031]k
c
表示等高线曲率；
[0032]d表示栅格尺寸cellsize
[0033]x和y分别表示形态学连续解空间中的地表曲面在二位平面上投影的横坐标和纵坐标；
[0034]z表示形态学连续解空间中的高程；
[0035]p表示z对x的一阶偏导；q表示z对y的一阶偏导；r表示Z对x的二阶偏导；s表示z对xy的二阶偏导；t表示z对y的二阶偏导；
[0036]步骤三，根据负DEM计算up_iGD8上坡向：
[0037]负DEM数据是将预处理后的正DEM进行数值的反操作得到，即让二维数组demData的每个元素减去全部数据中的最小值，就可得到负DEM二维数组；
[0038]在正DEM中计算每个栅格最可能的上坡方向如同在负DEM中计算最陡流出方向，即原DEM的up_iGD8上坡向，具体的计算过程与步骤2.1相同；
[0039]步骤四，根据D8计算汇水面积并设置阈值确定沟道截断矩阵：
[0040]在步骤2.1中得到了每个栅格的流向，根据此流向进行河网提取，所述的河网提取采用流域拓扑方法来计算，通过计算得到的汇水面积，设置一个河网阈值，即当汇水面积矩阵中的元素大于规定阈值后标记该元素的位置为沟道截断位置；依次将满足条件的元素保留在新的矩阵中，其数值记为1，其它元素记为0，最终生成沟道截断矩阵；
[0041]步骤五，根据up_iGD8上坡向计算坡度并确定坡度截断矩阵：
[0042]在步骤三中得到了每个栅格的up_iGD8上坡向，坡度是根据坡向的进一步计算，选择中心栅格的上坡向所指的邻域栅格，计算两者高程差Δh和水平距离Δd，然后取反正切，记为式VI：
[0043]slope＝180/π
·
arctan(Δh/Δd)
ꢀꢀꢀꢀ
式VI；
[0044]当坡向为坐标系方向时，水平距离Δd为1倍的栅格尺寸；当坡向为对角线方向时，水平距离Δd为倍的栅格尺寸；坡度截断的标准是：
[0045]如果当前栅格的坡度大于2.86
°
时，坡度变化率超过50％发生截断；
[0046]如果当前栅格的坡度小于等于2.86
°
时，坡度变化率超过70％发生截断；
[0047]将满足上述条件的栅格的位置记录下来，保留在新的矩阵中，其数值记为1，其它元素记为0，最终生成坡度截断矩阵；
[0048]步骤六，根据up_iGD8上坡向计算上坡坡度线长度：
[0049]采用局部坡向合并策略计算上坡坡度线长度；
[0050]所述的局部坡向合并策略的具体步骤如下：
[0051]Step601，从每个栅格开始建立辅助队列，按照上坡向所指的方向依次将邻域栅格加入队列，并将Δd进行累加；
[0052]Step602，根据Step601的方法，每累加3个步长，记录初始栅格和本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种融合地表形态特征的坡长提取方法，其特征在于，该方法采用的坡长计算公式为：式中：λ表示坡长；λ0表示局部最高点的坡长；l表示坡度线长度；k
c
表示等高线曲率。2.如权利要求1所述的融合地表形态特征的坡长提取方法，其特征在于，该方法在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表示栅格尺寸cellsizex和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/Δ...

【专利技术属性】
技术研发人员：张宏鸣，董良，葛晨宇，许伊昆，刘子涵，樊晓，
申请(专利权)人：西北农林科技大学，
类型：发明
国别省市：