[发明专利]一种融合地表形态特征的坡长提取方法

权利要求书

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遍历完全部数据后得到一个汇水面积矩阵，记录了每个栅格元素的累计流量。