技术领域
[0001] 本
发明涉及数字地形分析技术，具体涉及一种基于DEM的平地河网水流方向估算方法。
背景技术
[0002] 平地区域流向确定是数字地形分析与
地貌计量学中的提取河网的一个难题，尤其在地形因子提取中，长期以来缺少解决该问题有效的估算方法，进而限制地形因子相应参数的获得，这同时也限制了更大空间范围的
土壤侵蚀模型的应用。
[0003] 迄今为止，学者们已提出了不少从DEM中提取河网数据的
算法，可分为:1)基于局部地表特征的方法，即局部方法。最典型的是移动窗口算法，即使用一个矩形窗口对DEM格网阵列进行扫描，在窗口寻找凹形或“V”形区域，认为这些区域底部的单元格是沟谷所在
位置，但这种方法无法获取连续的河网；2)通过全局地表信息，建立一致的地表特征结构以获取地貌特征的方法，即全局方法。典型的有坡面径向水流算法，利用水流特征，能产生连续河网，通过计算汇水量实现河网的提取；还有将上述两种方法结合起来的方法。
[0004] 上述诸方法，主要参考平地周围高程，获得的流向几乎一致，提取河网时，易产生平行的
流线，无法提取清晰而独立的河网，而对平地采用逐级抬高增加高程来构造微坡面的方法，使得平地流向处理效率不高，并且DEM的微小
变形会使平地周围出现洼地，对结果产生较大影响。
发明内容
[0005] 针对上述
现有技术中存在的问题，本发明的目的在于，提供一种平地河网水流方向估算方法，以单流向算法为
基础，较准确地计算平行流域的流向，提高基于
地理信息系统的区域土壤侵蚀研究、水土保持环境效应评价、流域水文分析等的应用效率，满足
水文学、地理学的工作需要。
[0006] 为了实现上述任务，本发明采用以下技术方案：
[0007] 一种基于DEM的平地河网水流方向估算方法，包括以下步骤：
[0008] 步骤一，载入DEM数据，找出非平地区域的流向，根据已有的流向计算出汇水面积；
[0009] 步骤二，找出所有平地汇水出口栅格，具体为：
[0010] 遍历DEM中所有栅格数据，对每一个栅格进行以下处理：
[0011] 如果当前栅格有流向，则判断当前栅格周围的8个栅格中有没有流向、高程均与当前栅格相等的栅格a，如果有则将该栅格a加入到一个数组A中，则数组A中的栅格即为所有的平地汇水出口栅格；
[0012] 步骤三，对于每一
块平地区域，找出该平地区域在数组A中对应的所有平地汇水出口栅格，对这些栅格进行排序，选择高程最低的一个栅格作为该平地区域的汇水出口栅格b，将每一块平地区域筛选出的汇水出口栅格b放入一个数组B中；
[0013] 步骤四，确定每一块平地区域的流向并完善汇水
[0014] (1)从数组B中取出一个栅格b0，对于栅格b0对应的平地区域，找出这块平地区域中所有没有流向的栅格放在一个数组C中；
[0015] (2)将数组C中的栅格按照汇水量的由大到小进行排序，排序后选择汇水量最大的n个栅格，即为干流和支流的个数；
[0016] (3)对数组C中选出的n个栅格采用最
短路径算法进行处理，得到该平地区域的干流和支流；
[0017] (4)将数组C中除了选出n个栅格之外的其余栅格，采用就近原则算法进行处理，给平地区域中尚未分配流向的栅格进行流向的分配。
[0018] 进一步地，所述的最短路径算法步骤包括：
[0019] (3-1)从数组C中选出的n个栅格中汇水量最大的栅格开始，将该栅格作为平地入口点栅格，遍历出所有走一步能达到的位置，并将这些可以到达的位置放入一个平地流经队列中；
[0020] (3-2)对于平地流经队列中的每一个位置元素，依次遍历每一个元素可到达的位置，将这些可到达的位置不重复地记录到平地流经队列中，并记录这些可到达位置的父
节点位置；
[0021] (3-3)当遇到平地区域的汇水出口栅格时，停止寻找当前元素可到达的位置，并从平地出口栅格开始回溯，寻找每个当前位置的
父节点，并赋予父节点的水流方向和累积流量，直到回溯到平地入口点栅格，得到平地区域的干流；
[0022] (3-4)找到平地区域的干流后，从数组C中选出的n个栅格中选择另外的栅格作为平地入口点栅格，按照3-1至3-3相同的方法寻找平地区域中的支流，并将这些支流汇聚在离自身最近的流路上。
[0023] 进一步地，所述的就近原则算法的步骤包括：
[0024] (4-1)设置Rank数组，存放平地区域中已经有流向的栅格；
[0025] (4-2)遍历Rank数组，记当前遍历到的栅格为d，如果栅格d邻接的8个栅格中有在数组C中且没有流向的栅格e，则使栅格e流向d并将栅格e加入到Rank数组中；
[0026] (4-3)当Rank数组中的元素个数等于数组C中元素个数时停止遍历，则平地区域中所有栅格的流向都分配完毕。
[0027] 本发明具有以下技术特点：
[0028] 1.从本算法设计来看，本算法以平地栅格的初始累计流量为出发点，按照最短路径和就近原则的思想，以方向加权的遍历顺序分别确定平地栅格的流向。在最短路径方法中，规定方向加权可确保路径最短，从而使得程序运行速率较快。本算法在计算流向过程中也完成了汇水面积的计算，因此该方法适合用于流向及汇水面积同时计算，因此
软件工具生成的结果为流向和汇水两个文件。考虑效率问题，单独看待计算平地水流时间，该算法快于Barnes算法，但不适合单独计算水流方向。
[0029] 2.从本算法运行效果来看，本算法和Barnes算法以及J&D算法进行比较，本算法可以明显的看到平地中的主流，以及其他支流汇聚到主流或者附近河流之上，也解决了平地当中一部分水流方向散乱乱流的问题；Barnes算法计算平地流向和本算法的结果基本一致，但仍存在少数的平行流域，并且出现若干无流向点，导致形成不连续河网；J&D算法计算平地流向生产汇水面积也存在平行流问题。在运行时间上，Barnes算法通过增加DEM高程，不断重复计算和
迭代，因此计算量很大，运行时间较长；由于本算法不对DEM数据进行
修改，不会对地形造成任何微小的变化，且效率较高。
[0030] 3.本文在前期的研究基础上，通过与J&D方法、Barnes等人方法的对比，指出了它们对平地流向处理的不足，本发明直接对平地进行流向判断，不需对DEM进行修改，提出以平地栅格的初始汇流累积量为参考，采用最短路径方法和就近原则方法来计算平地的水流方向。该方法适合提取流向和汇水面积，结果合理，效率较高，具有较好的运行效果和可行性，是对数字地形分析方法的有益尝试。
附图说明
[0032] 图2为最短路径算法示意图；
[0033] 图3为就近原则算法示意图；
[0034] 图4(a)为本方案算法的平地流向示意图；
[0035] 图4(b)为Barnes算法平地流向示意图；
[0036] 图4(c)为本方案算法汇水面积示意图；
[0037] 图4(d)为Barnes算法汇水面积示意图；
[0038] 图5为不同算法下平地汇水面积示意图；
具体实施方式
[0039] 平地流向处理是水文应用中的一个重要环节，其目的在于合理确定水流方向，得到
精度较高的河网水系及数字流域特征，为水文、洪水、侵蚀等领域模拟提供基础数据。基于数字高程模型(DEM)的水文分析中，如河网水系提取、区域侵蚀制图等需要将栅格赋予水流方向。最棘手的问题就是平地区域的处理问题，前人的算法当中多以修正DEM值来确定平地水流方向，在合理性和效率方面还有待进一步的提高。本发明借鉴前人研究的基础上，根据平地处理的基本特征，设计并实现了平地水流方向处理的算法，以平地的初始汇流累积量为参考，提出最短路径算法和就近原则算法，并对相关结果进行了全面细致的分析和对比，具体方法如下：
[0040] 步骤一，载入DEM数据，根据水向低处流的原理，找出非平地区域的流向，根据已有的流向计算出汇水面积；这里由于平地区域没有流向，所以汇水面积是不完整的；每个栅格的汇水初始值记为单位1；
[0041] 步骤二，找出所有平地汇水出口栅格，具体为：
[0042] 遍历DEM中所有栅格数据，对每一个栅格进行以下处理：
[0043] 如果当前栅格有流向，则判断当前栅格周围的8个栅格(东，南，西，北，东南，西南，西北，东北)中有没有流向、高程均与当前栅格相等的栅格a，如果有则将该栅格a加入到一个数组A中，则数组A中的栅格即为所有的平地汇水出口栅格；如果当前栅格没有流向，则跳过继续遍历下一个栅格；这里的汇水出口栅格是指在一块平地区域中汇水最终流向平地外的栅格；对数组A中的栅格进行排序；
[0044] 步骤三，由于整个DEM数据中可能存在多块平地区域，因此数组A中集中了所有平地区域的汇水出口栅格。很显然在步骤二找到的平地汇水出口栅格中，可能一块平地区域会有好几个汇水出口栅格，我们每一块平地只保留一个出口栅格。由于前面的排序过程，使一块平地区域的所有汇水出口栅格都在数组A中比集中(图中或许有不同的平地高程相同所以一块平地的所有汇水出口栅格不一定是连续的)，所以处理时是比较好处理的：对于每一块平地区域，找出该平地区域在数组A中对应的所有平地汇水出口栅格，对这些栅格进行排序，选择高程最低的一个栅格作为该平地区域的汇水出口栅格b，将每一块平地区域筛选出的汇水出口栅格b放入一个数组B中；
[0045] 例如数组A中的元素为[q1,q2,q3...w1,w2...e1,e2,e3,e4...]，其中第一块平地区域对应的汇水出口栅格为q1,q2,q3，第二块平地区域对应的汇水出口栅格为w1,w2等等。步骤二中的筛选是先将所有符合条件的栅格选出来放在数组A中，然后每一块平地区域仅从其在数组A中对应的多个汇水出口栅格中选出一个作为这块平地区域的汇水出口栅格；如第一块平地区域对应的栅格按高程由大到小排序为q3,q1,q2，则选择q2作为该平地区域的汇水出口栅格b，并将这些数据放在一个数组B中，供后续处理。由于在步骤二的汇水过程中平地区域是没有流向的所以得到的汇水必然会在平地区域断流，在后续的平地流向确定中以汇水量的大小进行确定干流和支流，有了主流和支流有利于平地中水的汇聚(汇聚有效的避免了平行流)，因此我们只有处理了高处平地的流向并完成了平地的汇水及汇水留到下一块平(没有流向不能继续流了)的或流出界外，才能处理较低的平地。
[0046] 步骤四，确定每一块平地区域的流向并完善汇水
[0047] 对每一块平地区域，均按照下述过程进行干流、支流的确定：
[0048] (1)经过前面的步骤筛选，数组B中当前的每一个栅格仅对应一块平地区域；从数组B中取出一个栅格b0，该栅格当前的汇水量记为h0，对于栅格b0对应平地区域，找出这块平地区域中所有没有流向的栅格放在一个数组C中；C中只有一个栅格有流向，就是这块平地区域出水口栅格；
[0049] (2)将数组C中的栅格按照汇水量的由大到小进行排序，排序后选择汇水量最大的n个栅格，即为数组C的前n个元素，就是干流和支流的个数；
[0050] 在实际操作时，最多从C中选择4个栅格作为干流和支流上的栅格点，这4个栅格必然会是平地边缘的栅格，汇水量越大表示外部汇聚的水越多，则间接地说明的这个地点离水源头越远，因此选择汇水量最大的栅格作为干流上的栅格点；这里选用n＝4是考虑一般情况下平地区域面积不会太大，从实际的地形中我们知道一个汇水的交汇处一般不会超过4条较为明显的直流，因此这里选择了最多4个。
[0051] (3)对数组C中选出的n个栅格采用最短路径算法进行处理，得到该平地区域的干流和支流；最短路径算法即从C中选取一个元素作为起点，在C中找到有流向的栅格，从起点沿着这条流路流到有流向的栅格(第一个点必然会以最短路径流向该平地汇水出口栅格，从而形成主流，第二次的最短路径则是汇聚到主流上，第三条流向主流或是第二大主流，第四条会汇聚到前面的几条流径上)，将每一个有了流向的栅格加入到Rank数组中，Rank数组中第一个元素是平地出口点栅格，第二个则是与平地出口点栅格相邻且在干流上的栅格，后面的也是相当于平地修复后这块平地区域中流量从大到小的顺序排列。图2为寻找一条最短路径的过程，数字0表示平地入口点栅格；实心点表示平地出口点；小数点左边的数字表示已遍历过的节点序号；小数点之后的数字表示寻找方向的优先级，直边方向的优先级依次为1、2、3、4，斜边方向的优先级依次为5、6、7、8。
[0052] 最短路径算法可参考如下：
[0053] (3-1)从数组C中选出的n个栅格中汇水量最大的栅格开始，将该栅格作为平地入口点栅格，遍历出所有走一步能达到的位置，并将这些可以到达的位置放入一个平地流经队列中；
[0054] (3-2)对于平地流经队列中的每一个位置元素，依次遍历每一个元素可到达的位置，将这些可到达的位置不重复地记录到平地流经队列中，并记录这些可到达位置的父节点位置；
[0055] (3-3)当遇到平地区域的汇水出口栅格时，停止寻找当前元素可到达的位置，并从平地出口栅格开始回溯，寻找每个当前位置的父节点，并赋予父节点的水流方向和累积流量，直到回溯到平地入口点栅格，得到平地区域的干流；
[0056] (3-4)找到平地区域的干流后，从数组C中选出的n个栅格中选择另外的栅格作为平地入口点栅格，按照3-1至3-3相同的方法寻找平地区域中的支流，并将这些支流汇聚在离自身最近的流路上。
[0057] (4)将数组C中除了选出n个栅格之外的其余栅格，采用就近原则算法进行处理，给平地区域中尚未分配流向的栅格进行流向的分配。图3是一小块平的最短路径和就近原则后的流向结果，平地中深黑色的是用最短路径计算的，浅黑色的是有就近原则算的，数字表示就近原则中找到的先后顺序。
[0058] 就近原则算法的步骤包括：
[0059] (4-1)设置Rank数组，存放平地区域中已经有流向的栅格；
[0060] (4-2)遍历Rank数组，记当前遍历到的栅格为d，如果栅格d邻接的8个栅格中有在数组C中且没有流向的栅格e，则使栅格e流向d并将栅格e加入到Rank数组中；
[0061] (4-3)当Rank数组中的元素个数等于数组C中元素个数时停止遍历(此时也就是这块平地区域中的每一个栅格都有了流向)，则平地区域中所有栅格的流向都分配完毕。
[0062] 由于前面平地没有流向导致汇水积在了平地区域，在经过平地区域汇水面积计算后，原来的平地区域汇水都积在了平地汇水出口栅格b0处，设此时b0处的汇水量为h1，则汇水增量h2＝h1-h0；根据已有的流向将h2流下去，直到流出界外或没有流向了(又遇到了平地)。
[0063] 算法的分析与比较：
[0064] 图4为不同算法下平地流向及汇水面积示意图，其中图4(a)为本方案算法的平地流向示意图；图4(b)为Barnes算法平地流向示意图；图4(c)为本方案算法汇水面积示意图；图4(d)为Barnes算法汇水面积示意图。
[0065] 这四幅图都是韭园沟10DEM数据中的一处平地的处理结果图，通过这几幅图可以很显然的看出，Barnes算法在处理平的流向的是后没有考虑汇聚问题，导致图四中的汇水出现(看到两条支流已经流到一起了，然后又分开了，这显然不符合现实情况)。
[0066] 图5(a)为Barnes算法下局部地区流向示意图，图中数字代表：1:西2:西北3:北4:东北5:东6:东南7:南8:西南-1:无方向；通过图可以很明显的看到该算法的处理结果中出现了没有流向的栅格，这样汇水计算以及后续的工作就无法完成了。
[0067] 图5(b)为Barnes算法的汇水面积示意图；由于平地中有栅格没有流向所以汇水断断续续没有连成一个完整的水系。
[0068] 图5(c)为本方案算法汇水面积示意图；可以很明显的看到本算法的汇水结果是完整的，没有出现断的地方；
[0069] 图5(a)和图5(b)中右下
角的
颜色标示中，颜色标示汇水面积的大小，颜色越黑表示值越大，颜色越浅表示值越小。
[0070] 通过图4和图5可以看出本算法解决了现有算法中水到了汇聚的地方却没有汇聚的问题，也解决了现有算法结果中出现无流向的问题。
