[发明专利]基于分水岭机制的自适应斜坡单元提取方法

权利要求书

1.基于分水岭机制的自适应斜坡单元提取方法，包括如下步骤：

步骤一：首先，对数字高程图进行预处理操作，生成反向高程图，并对两张高程图进行填洼，通过grass gis的python库函数计算得到高程图的坡向栅格图aspect、坡度栅格图slope以及流向栅格图dir，通过栅格计算器计算aspect的cos值和sin值即坡向的水平分量sinno和垂直分量cosno，后续计算坡向一致性用到，使用汇流累积算法分别计算正反高程图的汇流累积量su_acc和re_acc，汇流累积量算法如下：

输入流向栅格图dir

(1.1)初始化流入量矩阵NIDP为dir相同尺寸的全零矩阵，dir中nodata的部分，NIDP也置为nodata；

(1.2)使用NIDP-base算法，遍历dir计算NIDP，周围每有一个cell流入遍历到的cell，该cell的值加1；

(1.3)初始化汇流累积量矩阵FlowAccu，与dir尺寸相同的全1矩阵；

(1.4)遍历NIDP矩阵，找到NIDP值为0的cell作为起点，根据dir矩阵的流向，流到NIDP中不为1的位置，并且将该位置的NIDP的值减一，FlowAccu矩阵的值加上流出位置的值。

(1.5)将FlowAccu导出tif

步骤二：通过grass gis中的watershed分水岭算法，设定一个最小面积阈值(通过二分法和设定的斜坡单元密度来寻找一个最优值)，提取得到初始粗斜坡单元slu，大部分实际应用中不会关注平坡区域，因此提取出坡度小于五度的平坡区域mask，剔除掉mask中细碎的部分，裁剪掉slu中mask部分，如果有行政边界矢量图就叠加到slu上，对每个斜坡单元根据行政区划编号，如果没有行政边界矢量图则将每个斜坡单元的行政区划编号设为0，溶解因计算误差带来的细碎斜坡单元，最后对斜坡单元边界线进行平滑就可以得到斜坡单元的最初版本，这个版本的斜坡单元密度比目标密度要稍大，以便后续的细化溶解。

步骤三：提取斜坡单元以及线段特征，首先根据步骤一的特征栅格图提取斜坡单元对应的面特征，如面积、坡向、行政区划等，再构造斜坡单元的邻接关系表提取邻接单元交线，通过su_acc和re_acc提取交线的汇流累积量，将su_acc和re_acc中的值量化后取两者中的较大值，得到交线的汇流累积量特征。

步骤四：遍历斜坡单元，根据步骤三得到的面积特征，判断斜坡单元是否需要溶解，计算需要溶解的斜坡单元每条边的坡向一致性，并挑选出满足坡向一致性阈值的边，坡向一致性阈值计算如下：

其中th_consist是某个斜坡单元的阈值，k为阈值控制系数，越小越容易被溶解，area_norm为某个斜坡单元面积标准化的结果。再从这些边中筛选出与邻接斜坡单元行政区划一致的边，挑选出汇流累积量最大的边作为要溶解的边，存入fush_line中。

步骤五：输入fush_line，找到可以进行连续溶解的斜坡单元，减少迭代次数，left表示线段左边的斜坡单元，right为右边。

(5.1)初始化is_fush为-1，表示未被修改；

(5.2)遍历fush_line，跳过is_fush为1的线段；

(5.3)将线段left的邻居中未被修改的线段加入队列中；

(5.4)队列中的线段出队列，并把融合标记设为left的id，is_fush置为1，并将该线段的邻居线段加入队列，直到队列为空；

(5.5)遍历结束，更新fush_line；

更新fush_line可以使得一些连续的斜坡单元可以同时溶解，使得溶解迭代的次数减少，提高溶解的速率。

步骤六：将fush_line中线段的left和right的斜坡单元的溶解标记都设置为该线段的溶解标记，然后将相同溶解标记的斜坡单元进行溶解，计算溶解过后的斜坡单元的密度，如果密度满足目标密度或者达到最大迭代次数，将斜坡单元输出；否则返回步骤三，提取溶解过后的斜坡单元特征，再次进行溶解，需要注意的是，为了提高计算速度，没有变化的交线特征不需要重新提取。