一种基于地表凹度指数的沟谷网络自动提取方法

1.本发明涉及地形分析领域，具体是一种基于地表凹度指数的沟谷网络自动提取方法。


背景技术：

2.沟谷是指地球表面由流水冲击侵蚀形成的一种狭窄凹地，在流水侵蚀严重的地区，不同级别的沟谷镶嵌、互联形成一个地区的沟谷网络。在流水侵蚀严重的区域，例如我国黄土高原地区，准确提取沟谷网络对农业生产规划、水土保持、防沙治沙等工作，以及地貌发育演化等科学研究都具有重要意义。
3.当前沟谷网络的提取主要方式是基于dem数据，利用地表径流漫流模型，模拟地表径流在地表的流动来产生水系，即沟谷网络。其基本过程是先对dem数据进行填洼处理，消除dem数据中的局部洼地，然后依据高程差最大的原则确定每一个栅格单元上水流的方向，从而模拟出水流在一个区域表面上的汇流路线；然后，基于这个汇流路线计算每个栅格单元所经过的汇水单元的个数，统计得出的总数称为该栅格单元的汇流累积量；一个栅格的汇流累积量达到一定数值范围后就可以认为该栅格成为水系或沟谷的一部分，即通过判断栅格的汇流累积量值来确定沟谷的位置。该方法依据水文学汇流概念判别水流路径并模拟地表径流，有一定的模型基础，同时算法简单，能够直接产生连续的流路，是当前普遍采用的沟谷网络提取方法。
4.上述方法的一个关键步骤就是确定一个汇流累积量阈值，该阈值决定了一条沟谷水系的起始位置，也决定了一个区域的沟谷网络的数量和分布。因此关于沟谷水系提取的主要研究就集中在如何确定这个汇流累积量阈值。综合分析已有的研究，目前采用较多的方法有两种：第一种是实验法，即采用不同阈值提取沟谷网络，然后与已有的水系图、遥感影像等数据叠加比对，确定最适宜的汇流累积量阈值；第二种是寻找突变点法，主要思路是计算不同汇流累积量阈值下的沟谷网络，然后统计其总长度，绘制不同阈值参数与沟谷网络总长度之间的散点图，将散点图拟合成一条曲线，寻找曲线的突变点作为最佳阈值，而寻找突变点的典型方法为均值变点分析法。
5.以上两种方法存在明显的缺陷：第一，工作量较大，两种方法都需要设置不同阈值参数进行沟谷网络提取，同时后续的结果比对、突变点计算等也比较耗时；第二，当前在提取沟谷网络时所有方法都试图确定一个最优的汇流累积量阈值，而忽视了该阈值的可变性，即便在同一个区域，不同地表的侵蚀程度存在差异，沟谷起始位置对应的汇流累积量阈值参数也可能时是不同的。因此如果固定一个阈值参数，即便确定一个适宜的参数值，所提取的沟谷网络也存在较大误差。
6.为解决以上技术问题，本发明从沟谷的定义出发，利用其在地表中呈现的“狭窄凹地”特征，构建了一个基于dem数据计算的地表凹度指数，通过该指数首先确定沟谷在区域地表的起始位置，然后结合地表径流漫流模型，从沟谷起始位置，顺着汇流路径直到流域出水口。这种方式在提取沟谷网络过程中不再固定单一的汇流累积量阈值参数，而是从地表
形态的角度先寻找确定沟谷起始位置，然后结合地表径流漫流模型也能得到连续完整的沟谷网络。本发明方法易于操作，所提取的沟谷网络与实际地表中的沟谷网络更加吻合，在农业规划、水土保持规划等实践项目和地貌学相关科学研究领域具有重要的应用价值。


技术实现要素：

7.本发明针对上述问题，设计了一种能有效反映地表下凹特征的地表凹度指数，并基于规则格网dem数据实现该指数的计算，然后借助该指数实现沟谷网络的快速、准确提取。该方法包含基于dem数据的地表凹度指数计算与地表凹度指数辅助下的沟谷网络提取两个关键技术环节，能够有效解决当前沟谷网络提取中汇流累积量阈值参数难以有效确定的技术难题。
8.本发明采用的技术方案如下：
9.基于地表凹度指数的沟谷网络自动提取方法，包括如下步骤：
10.s1，输入dem数据，通过洼地填充、流向计算、汇流累积量计算等过程计算得到每个栅格的汇流累积量值，并依据汇流累积量值是否大于指定参数值对汇流累积量数据进行重分类，若大于，则取值为1，否则取值为0，得到研究区的第一个二值化栅格数据；
11.s2，基于dem数据，计算研究区的地形凹度指数，并依据地形凹度指数是否大于指定参数值对凹度指数进行重分类，若大于，则取值为1，否则取值为0，得到研究区的第二个二值化栅格数据；
12.s3，将步骤s1得到的研究区第一个二值化数据和步骤s2得到的第二个二值化栅格数据进行叠加乘法运算，得到研究区第三个二值化数据；
13.s4，将步骤s3得到的二值化数据，选择值为1的栅格，将其转换为矢量线格式，再提取每一条矢量线的两个端点；
14.s5，以步骤s4得到的矢量线端点，获取端点上的汇流累积量值，通过窗口分析的方式，搜寻汇流累积量值增加的栅格，直至流域出水口，最终得到研究区的沟谷网络。
15.进一步地，所述步骤s2的具体过程为：
16.s21，对研究区dem数据，读取dem数据格网分辨率h、行数m、列数n等基本信息；
17.s22，从第一个栅格开始，判断当前栅格的位置，如果当前栅格位于边界上，即其行号为1或m，或者列数为1或n，则对该栅格位置的地表凹度指数赋值为0，否则，执行s23步骤；
18.s23，当前栅格不在边界上，则以当前栅格为中心，与其左边、右边、上边、下边、左上方、右下方、右上方、左下方共计9个栅格组成一个3
×
3栅格窗口；
19.s24，根据窗口内不同栅格的空间位置关系，以及格网分辨率h，分别计算中心栅格与其左边和右边栅格，中心栅格与其上边和下边栅格，中心栅格与其左上方和右下方栅格，中心栅格与其右上方和左下方栅格的空间夹角的余弦值；
20.s25，比较步骤s24计算得到的四个余弦值，将最大值加上1作为中心栅格的地表凹度指数；
21.s26，重复步骤s22—s25，遍历所有栅格，完成每一个栅格地表凹度指数的计算；
22.s27，设置阈值参数，对步骤s21—s26计算得到的地表凹度指数进行判断，如果某栅格地表凹度指数大于该阈值参数，则将该栅格赋值为1，否则赋值为0，得到二值化数据。
23.进一步地，所述步骤s5的具体过程如下：
24.s51，定义一个与研究数据栅格分辨率、行列数等基本信息一致的栅格数据b，并赋初值0；
25.s52，读入步骤s4得到的点数据，获取每个点对应空间位置上的汇流累积量值，然后对对每一个点，依次执行步骤s53—s55；
26.s53，读取一个点，分别确定汇流累积量栅格数据a和步骤s51定义的栅格数据b中与该点相同位置的栅格单元，并将栅格数据b中对应的栅格单元标记为1；
27.s54，以步骤s53中确定的栅格数据a中栅格单元为当前处理栅格，在其周围(即左边、右边、上边、下边、左上方、右上方、左下方和右下方)栅格中搜寻汇流累积量值大于当前处理栅格单元汇流累积量值得栅格单元，若有多个，则选择最大的栅格作为当前处理栅格单元，同时将栅格数据b中与之对应的栅格单元标记为1；
28.s55，重复执行步骤s54，在执行过程中栅格数据b中的栅格单元可重复标记，直到栅格数据a中找不到汇流累积量值大于当前处理栅格单元汇流累积量值得栅格单元，说明已经到达流域出水口，则选择其他尚未处理的端点，重复执行步骤s53—s55；
29.s56，所有点处理完成后，栅格数据b中标记为1的栅格即组成了研究区的沟谷网络，将其转为实力数据输出得到研究区的沟谷网络数据。
30.本发明的技术特点及有益效果：
31.(1)本发明提出的地表凹度指数能够准确描述地形的下凹形态特征，该指数定义明确，且基于格网dem数据容易计算，可以为地表形态分析、数字地形分析等科研领域提供技术参考和方法借鉴。
32.(2)本发明提出的以地表凹度指数为辅助的沟谷网络提取流程，一方面结合了经典技术流程中的基本分析过程，另一方面通过计算地表的下凹形态来确定沟谷的起始位置，解决了经典方法中汇流累积量阈值参数难以确定的问题。
33.(3)本发明提出的以地表凹度指数为辅助的沟谷网络提取流程，不再以固定的汇流累积量阈值参数作为确定沟谷的唯一标准，而是将沟谷的形态特征融入进来，显著提高了地形变化复杂区域提取的沟谷网络的准确度。
附图说明
34.图1为本发明方法的流程示意图。
35.图2为本发明方法中在栅格单元之间构建的局部直角坐标系示意图。
36.图3为本发明方法中栅格单元中心点组成的向量示意图。
37.图4为本发明方法中寻找沟谷栅格过程示意图。
38.图5为本发明方法中沟谷栅格单元标记示意图。
39.图6为本发明方法中生成的沟谷示意图。
40.图7为本发明实施例中实验区计算的地表凹度指数示意图。
41.图8为本发明实施例中实验区计算的沟谷起始备选点示意图。
42.图9为本发明实施例中实验区提取的沟谷网络示意图。
43.图10为本发明实施例中实验区根据沟谷网络提取的沟谷起始点示意图。
44.图11为本发明实施例中实验区对比方法提取的沟谷网络示意图。
具体实施方式
45.以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
46.如图1所示，本发明基于地表凹度指数的沟谷网络自动提取方法的具体步骤如下：
47.步骤1：读入研究区dem数据，依次实施洼地填充、流向计算以及汇流累积量计算等过程，得到每一个栅格的汇流累积量值；
48.步骤2：以10为参数，判断步骤得到的汇流累积量栅格数据，如果栅格单元汇流累积量值大于10，则该栅格赋值为1，否则赋值为0，得到第1个二值化栅格数据；
49.步骤3：读入研究区dem数据，依次计算每个栅格单元的地表凹度指数，其中若待计算栅格位于边界上，即其处在第1行，或者最后一行，或第1列，或者最后一列，则将该栅格的地表凹度指数赋值为0，否则按照步骤4计算；
50.步骤4：对处在非边界位置上的栅格单元，以当前处理栅格为中心，与其周围8个栅格组成一个3
×
3分析窗口，以中心栅格中心点为原点，分别以过原点的水平方向和垂直方向建立局部直角坐标系(见图2)；
51.步骤5：基于步骤4构建的局部直角坐标系，以坐标系原点为中心点，分别生成中心点与其周围8个点的向量(见图3)，图中为展示向量，将两个坐标轴省去。然后分别计算向量p1o与p8o，p2o与p7o，p3o与p6o，p4o与p5o之间夹角的余弦值，计算公式为：
[0052][0053]
式中，θ表示两个向量之间的夹角，(x1,y1,z1)和(x2,y2,z2)分别是两个向量的空间坐标。
[0054]
计算完成后，取四个夹角余弦值的最大值，再加上1作为当前处理栅格单元的地表凹度指数值，即:
[0055]
con＝max(cosθ1 cosθ2 cosθ3 cosθ4) 1
ꢀꢀꢀꢀꢀꢀ
(公式2)
[0056]
步骤6：重复步骤4—5，完成研究区地表凹度指数的计算；
[0057]
步骤7：对计算的研究区内的地表凹度指数，以0.8为临界参数，将地表凹度指数大于0.8的赋值为1，将地表凹度指数小于或等于0.8的赋值为0，得到第2个二值化栅格数据；
[0058]
步骤8：两个二值化栅格数据做乘法叠加运算，则两个二值化栅格数据在值均为1的情况下叠加运算的结果为1，否则为0，得到第3个二值化栅格数据；
[0059]
步骤9：在第3个二值化栅格数据中选择值为1的栅格，将连续分布的栅格处于两端的栅格转为矢量点数据；
[0060]
步骤10：对步骤9中每一条连续栅格所获取的两个端点数据，分别获取点位置上的汇流累积量值，并将累积量值大的那个点删除；
[0061]
步骤11：生成结果栅格数据，赋初始值为0；
[0062]
步骤12：对步骤10得到的点数据，对每一个点依次处理，计算该点在栅格数据中的位置，确定汇流累积量栅格数据和结果栅格数据中的对应栅格单元，并将结果栅格单元赋
值为1；
[0063]
步骤13：以汇流累积量栅格数据中确定的当前处理栅格为中心，在其周边寻找汇流累积量值大于当前栅格汇流累积量值的栅格，将其作为新的当前处理栅格，并将结果栅格数据中位置对应的栅格数据标记为1(见图4、图5)；
[0064]
步骤14：不断重复步骤13，直到找不到汇流累积量值大于当前栅格汇流累积量值的栅格，说明到达流域出水口，一条汇水路径搜寻完成；
[0065]
步骤15：重复步骤12-14，所有点数据处理完成后，结果栅格数据中标记为1的即沟谷所在的位置(见图6)；
[0066]
步骤16：将结果数据中值为1的连续栅格转为矢量线，得到研究区域的沟谷网络，并输出。
[0067]
实施例
[0068]
为了验证本发明方法在提取沟谷网络时的可用性，实施例选择陕北黄土高于某小流域为实验区，该实验区是典型的流水侵蚀地形，沟谷网络系统发育十分完善，可作为不同方法沟谷网络提取结果差异性分析的理想实验区域。
[0069]
步骤1：对示例区域，经过洼地填充、流向计算、汇流累积量计算等过程，对计算的汇流累积量以10为临界参数，大于10的赋值为1，否则赋值为0，得到示例区域的第一个二值化栅格数据；
[0070]
步骤2：对示例区域，计算地表凹度指数(结果见图7)，然后以0.8为临界参数，大于0.8赋值为1，否则赋值为0，得到示例区域的第二个二值化栅格数据；
[0071]
步骤3：将步骤1得到的二值化栅格数据和步骤2得到的二值化栅格数据做乘法运算，得到示例区域第三个二值化栅格数据；
[0072]
步骤4：将步骤3得到的二值化栅格数据转成矢量线，再取线段两端点，获取端点上的汇流累积量值，并将同一个线段转成的端点中汇流累积量值较大的那个删除，得到可用点数据(图8)；
[0073]
步骤5：基于步骤4得到的点数据，以点为起始位置，采用搜寻汇流累积量值增加的相邻栅格的形式，获取汇水路径线；
[0074]
步骤6：将搜寻完成的汇水路径线上的栅格转为矢量数据，得到示例区域的沟谷网络结果(图9)，然后基于沟谷网络提取示例区域的沟谷起始点(图10)；
[0075]
步骤7：对示例区域，采用对比方法(以均值变点分析法确定阈值的方法)得到的沟谷网络结果(见图11)。
[0076]
通过以上对比分析，基于均值变点分析法计算的汇流累积量阈值提取的沟谷存在两个明显问题：第一是出现了大量伪沟谷，主要出现在小流域上游区域，该区域沟壑发育还不够完善，平坦地形较多，容易形成伪沟谷；第二个问题则是有很多长度较短的沟谷没有提取出来，这是因为这些沟谷集水面积较小，远没有达到确定的汇流累积量阈值，因此无法得到有效提取。而本发明方法提取的沟谷第一个问题完全不存在，因为伪沟谷主要分布在地形平坦地区，这些地区的地表凹度指数近似为0，所以本发明方法提取的沟谷可以有效避免出现伪沟谷；第二个问题对于本发明方法也会部分存在，主要是一些沟谷发育还不充分，或者处在沟谷发育的初期，地表下凹形态不明显，因而本发明方法无法有效提取，但其无法识别的沟谷的数量显著小于汇流累积量阈值法。因此，综合以上分析可以本发明方法提取沟
谷网络结果更准确。