CN111125893B - 一种基于dem和流量汇集的非弥散水流路径模拟方法 - Google Patents

Abstract

本发明提供了一种基于DEM和流量汇集的非弥散水流路径模拟方法，属于数字地形分析技术领域。其技术方案为：通过根据高程从高到低逐渐确定各栅格流量汇集点位置，模拟流量从汇集点沿局部排水方向流动的过程，确定每一个DEM栅格单元的下游单元，最终连接出非弥散水流路径。本发明的有益效果为：本发明为基于DEM的水文模型提供了更好地模拟水流在地表的流动过程的方法。

Description

一种基于DEM和流量汇集的非弥散水流路径模拟方法

技术领域

本发明涉及数字地形分析技术领域，尤其涉及一种基于DEM和流量汇集的非弥散水流路径模拟方法。

背景技术

水流路径指的是水受重力作用在地表的流动轨迹，目前的分布式水文模型大多使用水流路径来估算水流从生成点到河道或流域出口的流动距离，进而获取山坡宽度、汇流时间等模型参数。由于当下的水文模型大多借助栅格型数字高程模型(Digital ElevationModel，DEM)描绘地表地形，数字化的水流路径也大多基于DEM提取得到。目前基于DEM的水流路径提取算法包含两个大类，即弥散型和非弥散型。其中，弥散型算法认为水流是发散的，发源于一点的水流会不断分离成多份，流向多个下游区域。受限于水流不断发散的特点，这类算法难以提供准确的流动距离和流域边界。因此，非弥散型算法在部分水文模型中必不可少。

基于DEM的非弥散型水流路径提取算法普遍给每一个DEM栅格单元指定一个下游单元。从起始的上游栅格起，逐步寻找下游单元，连接上下游相邻单元的中心点，就构成了一条完整的水流路径。因此，相邻单元间的局部水流路径被限制为8个方向之一。但是，真实环境下局部水流的排水方向应该属于0°～360°之间无穷个方向之一。除了早期的D8算法和Rho8算法只在指向8个相邻的栅格单元中心的方向中选取局部排水方向，目前已经有多种算法可以计算得到0°～360°之间任意角度的局部排水方向，例如Tarboton(1997)发表的Dinf算法以及Hooshyar等(2016)基于切面曲率对Dinf流向进行优化的算法。由于这些算法大多模拟水流自单元中心起的排水方向，但排水方向又可能不指向相邻栅格单元的单元中心，因此当模拟的水流到达相邻单元中心时，就没有了进一步指示流动轨迹的排水方向。因此，如何借助这类局部排水方向构建完整的非弥散水流路径值得研究。

Orlandini等(2003)提出了一种D8-LTD算法，该算法将Dinf流向调整至8方向之一，并将Dinf流向和最终方向之间的偏差向下游传递，辅助下游单元的局部流向校正过程，但是D8-LTD算法只以实现起始单元的路径最优为目标。对与下游单元，仅将其水流路径作为其上游单元的一部分，因而难以实现全局水流路径最优。

如何解决上述技术问题为本发明面临的课题。

发明内容

本发明的目的在于提供一种在DEM上将数值为0°～360°之间任意角度的局部排水方向拟合为连续水流路径的方法，使得得到的水流路径能大致曲折环绕真实的水流路径向下游沿伸，通过将上下游的水流不断汇集到一处，从上游到下游不断确定每个单元的下游单元，并将局部排水方向与下游单元中心点方向间的偏差传递至下游，在下游的水流聚集点的位置确定过程中考虑该偏差，发明了一种基于DEM和流量汇集的非弥散水流路径模拟方法。

本发明是通过如下措施实现的：一种基于DEM和流量汇集的非弥散水流路径模拟方法，其中，所述方法具体步骤如下：

步骤S1、加载DEM高程数据，将DEM数据读取为二维数组A，创建4个行列数与A相同的数组M、N、L、S，其中M、N分别用于存储每个栅格的水流汇集点与中心点间相对坐标的横坐标X，纵坐标Y的坐标值，L用于记录每个单元的上游集水面积，S用于记录下游单元的中心点方向，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，同时创建优先队列Q，Q能存储栅格单元的信息，在数组A中行号、列号、高程，并能根据栅格的高程将其从高到低排序，遇到栅格高程相同，则将后加入Q的单元排在先加入Q的单元的后方；

步骤S2、忽略第一列、最后一列、第一行、最后一行，扫描数组A中剩下的单元，判断数组A的每一个单元是否是无效值，若判别一个单元是无效值，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若判别一个单元不是无效值，则扫描其周边8个相邻单元，若相邻单元中存在无效值单元，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若相邻单元中没有无效值单元，则在数组M、N中将第p行、第q列的数值赋为0，在L中将第p行、第q列的数值赋为1，并将A中的当前单元插入优先队列Q；

步骤S3、取出优先队列Q最前端的单元C，读取单元C在数组A中的行数a，列数b，对应高程z值，若矩阵S的第a行，第b列已被赋值，则直接进行步骤S8，否则从步骤S4开始进行；

步骤S4、计算单元C的局部排水方向γ；

步骤S5、分别从数组M、N中读取各自第a行，第b列的数值y_c和x_c，y_c和x_c是单元C的水流汇集点相对单元C中心点的相对坐标W(x_c，y_c)，以C为中心，依次连接其周围8个相邻单元的中心点构成矩形窗口，水流从C的水流汇集点W(x_c，y_c)出发沿γ方向流动到达矩形窗口边缘，到达点记为R；

步骤S6、获取点R所在的单元T，获取T的中心点相对C单元中心点的方向，将其记入数组S的第a行，第b列；

步骤S7、若单元T在数组M中对应的同行同列单元不是无效值，则计算点R相对单元T中心点的相对坐标(x₁，y₁)，分别从数组L读取C和T对应行列的单元的数值，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值；

步骤S8、判别优先队列Q是否为空队列，若Q不是空队列，则返回步骤S3，否则扫描数组S中的每一个单元的方向数值，若数值不是无效值，则根据该方向值连接每个单元和下游单元的中心点，可得到DEM上模拟出的水流路径。

作为本发明提供的一种基于DEM和流量汇集的非弥散水流路径模拟方法进一步优化方案，所述步骤S7中，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，

其中，坐标(x，y)由以下方程得到：

其中，上述方程中c₀和c₁分别为单元T和C的集水面积值。

为了更好地实现上述发明目的，本发明还提供一种基于DEM和流量汇集的非弥散水流路径模拟方法，其中，所述方法具体步骤如下：

步骤S1、加载DEM高程数据，将DEM数据读取为二维数组A，创建4个行列数与A相同的数组M、N、L、S，其中M、N分别用于存储每个栅格的水流汇集点与中心点间相对坐标的横坐标X，纵坐标Y的坐标值，L用于记录每个单元的上游集水面积，S用于记录下游单元的中心点方向，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，同时创建优先队列Q，Q能存储栅格单元的信息，在数组A中行号、列号、高程，并能根据栅格的高程将其从高到低排序，遇到栅格高程相同，则将后加入Q的单元排在先加入Q的单元的后方；

步骤S1中的无效值具体为：

由于DEM包含的有效区域会出现不规则形状，存储时存储为矩形，通过创建额外的栅格，将区域填充为包含有效区域的最小矩形区域，额外创建的栅格普遍被赋予一个固定值，该固定值被定义为无效值，无效值被定义为-9999；

步骤S2、忽略第一列、最后一列、第一行、最后一行，扫描数组A中剩下的单元，判断数组A的每一个单元是否是无效值，若判别一个单元是无效值，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若判别一个单元不是无效值，则扫描其周边8个相邻单元，若相邻单元中存在无效值单元，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若相邻单元中没有无效值单元，则在数组M、N中将第p行、第q列的数值赋为0，在L中将第p行、第q列的数值赋为1，并将A中的当前单元插入优先队列Q；

步骤S3、取出优先队列Q最前端的单元C，读取单元C在数组A中的行数a，列数b，对应高程z值，若矩阵S的第a行，第b列已被赋值，则直接进行步骤S8，否则从步骤S4开始进行；

步骤S4、计算单元C的局部排水方向γ；

具体内容为：所述局部排水方向γ的范围为0°-360°，且可以为范围内的任意角度，使用Dinf方法得到的局部排水方向，将DEM的正北方向定义为0°，并顺时针增大；

步骤S5、分别从数组M、N中读取各自第a行，第b列的数值y_c和x_c，y_c和x_c是单元C的水流汇集点相对单元C中心点的相对坐标W(x_c，y_c)，以C为中心，依次连接其周围8个相邻单元的中心点构成矩形窗口，水流从C的水流汇集点W(x_c，y_c)出发沿γ方向流动到达矩形窗口边缘，到达点记为R；

具体内容为：以C为中心，设其周边8个相邻单元如图3中编号，横、纵(x、y)坐标正方向如图3所示，依次连接其周围8个相邻单元的中心点构成矩形窗口，如图3中以B₂、B₄、B₆、B₈单元中点为顶点的窗口，让水流从C的水流汇集点W(如图3中(a)中的W₁，图3中(b)的W₂)出发沿γ方向(如图3中(a)、3中(b)中黑色箭头方向)流动到达矩形窗口边缘，到达点记为R(如图3中(a)中的R₁，图3中(b)中的R₂)，此时点R相对单元C中点的相对坐标(x_R，y_R)由以下方程获取：

步骤S6、获取点R所在的单元T，获取T的中心点相对C单元中心点的方向，将其记入数组S的第a行，第b列；

具体内容为：

由于相邻单元T的中心点只能位于单元C的中心点的8个方向之一，对这8个方向直接以数字1～8编号，与图3中相邻单元的下标编码相同，如图4所示，以数字1～8分别表示正北、东北、正东、东南、正南、西南、正西、西北8个方向；若根据步骤S5得到的到达点R相对单元C中心点的相对坐标为(x_R，y_R)，则记入数组S的方向编号n为：

步骤S7、若单元T在数组M中对应的同行同列单元不是无效值，则计算点R相对单元T中心点的相对坐标(x₁，y₁)，分别从数组L读取C和T对应行列的单元的数值，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值；

具体内容为：

若步骤S5得到的R相对单元C中心点的相对坐标为(x_R，y_R)，步骤S6中读取的方向编号为n，则点R相对单元T中心点的相对坐标(x₁，y₁)可以用下面的方程获取：

如图3所示，x₁和y₁中最多只有1个不为0，随后从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，坐标(x，y)由以下方程获取：

上述方程中c₀和c₁分别为单元T和C的集水面积值，将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值(c₀+c₁)；

步骤S8、判别优先队列Q是否为空队列，若Q不是空队列，则返回步骤S3，否则扫描数组S中的每一个单元的方向数值，若数值不是无效值，则根据该方向值连接每个单元和下游单元的中心点，可得到DEM上模拟出的水流路径。

作为本发明提供的一种基于DEM和流量汇集的非弥散水流路径模拟方法进一步优化方案，所述步骤S7中，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，

其中，坐标(x，y)由以下方程得到：

其中，上述方程中c₀和c₁分别为单元T和C的集水面积值。

本发明的有益效果为：本发明将DEM栅格单元按照高程值从高到低排序，并从高到低给各栅格指定一个下游单元，此过程中根据上游来水确定各栅格的水流汇集点，并借助局部排水方向确定下游单元，通过依次连接上下游单元的栅格中心得到最终的水流路径。相较于传统的水流路径模拟方法，本发明的方法具有所确定水流路径更为准确的优势。因此，本发明为基于DEM的水文模型提供了高精度的水流过程模拟方法。

附图说明

图1为本发明的整体流程图。

图2为本发明DEM上方向角度设置，以正北方向为0°，并顺时针增大的示意图。

图3为本发明从流量汇集点沿局部排水方向流动示意图；

其中，对于单元C，其内流量的可流动范围是单元B₂、B₄、B₆、B₈中心点连接成的窗口，

图3中(a)为流量聚集点在中心单元C的中心，且到达窗口北部边缘的情况；

图3中(b)为流量聚集点不在中心单元C的中心，且到达窗口东部边缘的情况；

图4为本发明DEM上8个方向的编号对应示意图；以数字1-8分别表示正北，东北，正东，东南，正南，西南，正西，西北8个方向。

图5为本发明实施例1使用的DEM示意图，该DEM中第A行，第5列为无效值单元。

图6为本发明模拟实施例1中DEM上水流路径的过程示意图；

其中，图6中(a)中的A数组即为DEM数组；

图6中(b)为根据步骤S2完成排序的优先队列Q；

图6中(c)为对DEM中高程最高的C3单元进行处理时考虑到的的3×3窗口；

图6中(d)为本发明操作结束后的数组M存储的数值；

图6中(e)为本发明操作结束后的数组N存储的数值；

图6中(f)为本发明操作结束后的数组L存储的数值；

图6中(g)为本发明操作结束后的数组S存储的数值；

图6中(h)为使用本发明得到的例1中DEM上的水流路径的可视化效果。

图7为本发明实施例2使用的DEM使用不同方法模拟的水流路径示意图；

其中，图7中(a)为使用的倾斜平面的DEM高程图；

图7中(b)为使用的倾斜凹面的DEM高程图；

图7中(c)和图7中(d)分别为使用经典D8方法得到的平面、凹面上的水流路径；

图7中(e)和图7中(f)分别为使用D8-LTD方法得到的平面、凹面上的水流路径；

图7中(g)和图7中(h)分别为使用本发明得到的平面、凹面上的水流路径；

在图7中(c)至图7中(h)中，每幅图内从左下向右上的5条直线为间距为8m的等高线，其余的折线均为还原的水流路径。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。当然，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

需要说明的是，在不冲突的情况下，本发明创造中的实施例及实施例中的特征可以相互组合。

在本发明创造的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明创造和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明创造的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明创造的描述中，除非另有说明，“多个”的含义是两个或两个以上。

在本发明创造的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明创造中的具体含义。

实施例1

参见图1至图7，本发明提供其技术方案为，为了更好地实现上述发明目的，本发明还提供一种基于DEM和流量汇集的非弥散水流路径模拟方法，其中，所述方法具体步骤如下：

步骤S1、加载DEM高程数据，将DEM数据读取为二维数组A，创建4个行列数与A相同的数组M、N、L、S，其中M、N分别用于存储每个栅格的水流汇集点与中心点间相对坐标的横坐标X，纵坐标Y的坐标值，L用于记录每个单元的上游集水面积，S用于记录下游单元的中心点方向，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，同时创建优先队列Q，Q能存储栅格单元的信息，在数组A中行号、列号、高程，并能根据栅格的高程将其从高到低排序，遇到栅格高程相同，则将后加入Q的单元排在先加入Q的单元的后方；

步骤S1中的无效值具体为：

由于DEM包含的有效区域会出现不规则形状，存储时存储为矩形，通过创建额外的栅格，将区域填充为包含有效区域的最小矩形区域，额外创建的栅格普遍被赋予一个固定值，该固定值被定义为无效值，无效值被定义为-9999；

步骤S2、忽略第一列、最后一列、第一行、最后一行，扫描数组A中剩下的单元，判断数组A的每一个单元是否是无效值，若判别一个单元是无效值，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若判别一个单元不是无效值，则扫描其周边8个相邻单元，若相邻单元中存在无效值单元，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若相邻单元中没有无效值单元，则在数组M、N中将第p行、第q列的数值赋为0，在L中将第p行、第q列的数值赋为1，并将A中的当前单元插入优先队列Q；

步骤S3、取出优先队列Q最前端的单元C，读取单元C在数组A中的行数a，列数b，对应高程z值，若矩阵S的第a行，第b列已被赋值，则直接进行步骤S8，否则从步骤S4开始进行；

步骤S4、计算单元C的局部排水方向γ；

具体内容为：所述局部排水方向γ的范围为0°-360°，且可以为范围内的任意角度，使用Dinf方法得到的局部排水方向，将DEM的正北方向定义为0°，并顺时针增大；

步骤S5、分别从数组M、N中读取各自第a行，第b列的数值y_c和x_c，y_c和x_c是单元C的水流汇集点相对单元C中心点的相对坐标W(x_c，y_c)，以C为中心，依次连接其周围8个相邻单元的中心点构成矩形窗口，水流从C的水流汇集点W(x_c，y_c)出发沿γ方向流动到达矩形窗口边缘，到达点记为R；

具体内容为：以C为中心，设其周边8个相邻单元如图3中编号，横、纵(x、y)坐标正方向如图3所示，依次连接其周围8个相邻单元的中心点构成矩形窗口，如图3中以B₂、B₄、B₆、B₈单元中点为顶点的窗口，让水流从C的水流汇集点W(如图3中(a)中的W₁，图3中(b)中的W₂)出发沿γ方向(如图3中(a)、3中(b)中黑色箭头方向)流动到达矩形窗口边缘，到达点记为R(如图3中(a)中的R₁，图3中(b)中的R₂)，此时点R相对单元C中点的相对坐标(x_R，y_R)由以下方程获取：

步骤S6、获取点R所在的单元T，获取T的中心点相对C单元中心点的方向，将其记入数组S的第a行，第b列；

具体内容为：

由于相邻单元T的中心点只能位于单元C的中心点的8个方向之一，对这8个方向直接以数字1～8编号，与图3中相邻单元的下标编码相同，如图4所示，以数字1～8分别表示正北、东北、正东、东南、正南、西南、正西、西北8个方向；若根据步骤S5得到的到达点R相对单元C中心点的相对坐标为(x_R，y_R)，则记入数组S的方向编号n为：

步骤S7、若单元T在数组M中对应的同行同列单元不是无效值，则计算点R相对单元T中心点的相对坐标(x₁，y₁)，分别从数组L读取C和T对应行列的单元的数值，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值；

具体内容为：

若步骤S5得到的R相对单元C中心点的相对坐标为(x_R，y_R)，步骤S6中读取的方向编号为n，则点R相对单元T中心点的相对坐标(x₁，y₁)可以用下面的方程获取：

如图3所示，x₁和y₁中最多只有1个不为0，随后从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，坐标(x，y)由以下方程获取：

上述方程中c₀和c₁分别为单元T和C的集水面积值，将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值(c₀+c₁)；

步骤S8、判别优先队列Q是否为空队列，若Q不是空队列，则返回步骤S3，否则扫描数组S中的每一个单元的方向数值，若数值不是无效值，则根据该方向值连接每个单元和下游单元的中心点，可得到DEM上模拟出的水流路径。

作为本发明提供的一种基于DEM和流量汇集的非弥散水流路径模拟方法进一步优化方案，所述步骤S7中，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，

其中，坐标(x，y)由以下方程得到：

其中，上述方程中c₀和c₁分别为单元T和C的集水面积值。

为了更清楚地说明本发明的技术方案，本发明还提供了以下具体实例进行验证：

具体实例1：

以图5所示的5×5规格的DEM为例，本例中的无效值设为常用的-9999，并以Dinf方法提供局部排水方向，处理流程与结果如图6所示；具体步骤如下：

步骤S1、将DEM数据读取为数组A(图6中(a))，创建数组M(图6中(d))、N(图6中(e))、L(图6中(f))、S(图6中(g))，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，并创建按单元从高到低排序的优先队列Q(图6中(b))；

步骤S2、除第一列、最后一列、第一行、最后一行外，扫描数组A中剩下的单元，判断数组A的每一个单元是否是无效值-9999，该步骤没有找到无效值单元，随后判别该单元周边8个相邻单元是否存在无效值单元，此处发现A中单元B4与无效值单元A5相邻，则依照当前单元在数组A中的行列数，即第B行，第4列，在数组M、N、L、S中均将第B行，第4列的数值赋为无效值-9999，同时将A中除第一列、最后一列、第一行、最后一行、以及B4外的其他单元加入Q(图6中(b))；

步骤S3、取出优先队列Q最前端的单元，即第C行、第3列、高程为7的单元；

步骤S4、使用Dinf方法计算得到DEM中C3单元的局部排水方向γ为198.43°；

步骤S5、分别从数组M、N中读取各自第C行、第3列的数值y_c和x_c，二者数值均为0，二者就是单元C的水流汇集点W相对单元C中心点的相对坐标(0，0)，让水流从C的水流汇集点W出发沿γ方向流动到达矩形窗口边缘点R，由方程(3)和(4)得到此时点R相对单元C中点的相对坐标(x_R，y_R)为(-0.33，-1)；

步骤S6、获取点R所在的单元T，获取T的中心点相对C单元中心点的方向，由方程(5)得到该方向值为5，也即此处点R所在的单元T即A中的单元D3，将方向值记入数组S的第C行、第3列(图6中(g))；

步骤S7、若单元D3在数组M中对应的单元，即第D行、第3列不是无效值-9999，计算点R相对单元T中心点的相对坐标(x₁，y₁)，由方程(6)和(7)得到该相对坐标为(-0.33，0)，再分别从数组L读取C3和D3对应行列的单元的数值，即单元C和T的集水面积值，此处单元C和T的集水面积值均为1，并从数组M、N读取单元D3对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，即(0，0)，使用方程(8)和(9)得到单元D3的新流量汇集点的相对坐标(-0.165，0)，并将T在数组M(图6中(d))、N(图6中(e))中的对应值，即第D行、第3列分别修改为-0.165和0，再将D3在数组L(图6中(f))中的相同行列的对应值修改为当前单元C3和D3的集水面积的累积值2；

步骤S8、判别优先队列Q是否为空队列，此时优先队列不为空，于是返回步骤S3，取出新的头部单元D4继续计算，如此循环直至为空，最终数组M、N、L、S的赋值如图中6(d)-6中(g)所示，此时否则扫描数组S，即图6中(g)中的每一个单元的方向数值，若数值不是无效值-9999，则根据该方向值连接每个单元和下游单元的中心点，即可得到DEM上模拟出的水流路径，如图6中(f)。

具体实例2：

以图7中(a)、图7中(b)中的平面和凹面DEM为例，两个DEM的尺寸均为51×51，高程分布已在图中标出，分别使用传统的D8方法、D8-LTD方法以及本发明获取两个面上的非弥散水流路径，可视化的水流路径如图7中(c)-7中(h)所示，图7中的黑色线条为水流路径，灰色线条为等高线，理想的水流路径应该垂直于等高线流动，可以发现D8方法得到的水流路径非常直且不垂直于等高线(图7中(c)，图7中(d))，虽然D8-LTD在平面上应用效果良好，水流路径能大致垂直于等高线(图7中(e))，但是该方法在凹面上得到的水流路径出现了紊乱交叉的现象(图7中(f))，本发明成功弥补了已有方法的不足，在两个DEM面上展现了良好的应用效果(图7中(g)、图7中(h))。

本发明通过设置水流汇集点，使水流从水流汇集点出发沿局部排水方向流动确定下游单元，相较于传统方法，使用本发明得到的水流路径精度更高。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于DEM和流量汇集的非弥散水流路径模拟方法，其特征在于，所述方法具体步骤如下：

步骤S1、加载DEM高程数据，将DEM数据读取为二维数组A，创建4个行列数与A相同的数组M、N、L、S，其中M、N分别用于存储每个栅格的水流汇集点与中心点间相对坐标的横坐标X，纵坐标Y的坐标值，L用于记录每个单元的上游集水面积，S用于记录下游单元的中心点方向，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，同时创建优先队列Q，Q能存储栅格单元的信息，在数组A中行号、列号、高程，并能根据栅格的高程将其从高到低排序，遇到栅格高程相同，则将后加入Q的单元排在先加入Q的单元的后方；

步骤S2、忽略第一列、最后一列、第一行、最后一行，扫描数组A中剩下的单元，判断数组A的每一个单元是否是无效值，若判别一个单元是无效值，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若判别一个单元不是无效值，则扫描其周边8个相邻单元，若相邻单元中存在无效值单元，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若相邻单元中没有无效值单元，则在数组M、N中将第p行、第q列的数值赋为0，在L中将第p行、第q列的数值赋为1，并将A中的当前单元插入优先队列Q；

步骤S3、取出优先队列Q最前端的单元C，读取单元C在数组A中的行数a，列数b，对应高程z值，若矩阵S的第a行，第b列已被赋值，则直接进行步骤S8，否则从步骤S4开始进行；

步骤S4、计算单元C的局部排水方向γ；

步骤S5、分别从数组M、N中读取各自第a行，第b列的数值y_c和x_c，y_c和x_c是单元C的水流汇集点相对单元C中心点的相对坐标W(x_c，y_c)，以C为中心，依次连接其周围8个相邻单元的中心点构成矩形窗口，水流从C的水流汇集点W(x_c，y_c)出发沿γ方向流动到达矩形窗口边缘，到达点记为R；

步骤S6、获取点R所在的单元T，获取T的中心点相对C单元中心点的方向，将其记入数组S的第a行，第b列；

步骤S7、若单元T在数组M中对应的同行同列单元不是无效值，则计算点R相对单元T中心点的相对坐标(x₁，y₁)，分别从数组L读取C和T对应行列的单元的数值，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值；

步骤S8、判别优先队列Q是否为空队列，若Q不是空队列，则返回步骤S3，否则扫描数组S中的每一个单元的方向数值，若数值不是无效值，则根据该方向值连接每个单元和下游单元的中心点，得到DEM上模拟出的水流路径。

2.根据权利要求1所述的基于DEM和流量汇集的非弥散水流路径模拟方法，其特征在于，所述步骤S7中，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，

其中，坐标(x，y)由以下方程得到：

其中，上述方程中c₀和c₁分别为单元T和C的集水面积值。

3.一种如权利要求1或2所述的基于DEM和流量汇集的非弥散水流路径模拟方法，其特征在于，所述方法具体步骤如下：

步骤S1、加载DEM高程数据，将DEM数据读取为二维数组A，创建4个行列数与A相同的数组M、N、L、S，其中M、N分别用于存储每个栅格的水流汇集点与中心点间相对坐标的横坐标X，纵坐标Y的坐标值，L用于记录每个单元的上游集水面积，S用于记录下游单元的中心点方向，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，同时创建优先队列Q，Q能存储栅格单元的信息，在数组A中行号、列号、高程，并能根据栅格的高程将其从高到低排序，遇到栅格高程相同，则将后加入Q的单元排在先加入Q的单元的后方；

步骤S1中的无效值具体为：

由于DEM包含的有效区域会出现不规则形状，存储时存储为矩形，通过创建额外的栅格，将区域填充为包含有效区域的最小矩形区域，额外创建的栅格普遍被赋予一个固定值，该固定值被定义为无效值，无效值被定义为-9999；

步骤S2、忽略第一列、最后一列、第一行、最后一行，扫描数组A中剩下的单元，判断数组A的每一个单元是否是无效值，若判别一个单元是无效值，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若判别一个单元不是无效值，则扫描其周边8个相邻单元，若相邻单元中存在无效值单元，则依照当前单元在数组A中的行列数，第p行，第q列，在数组M、N、L、S中均将第p行、第q列的数值赋为无效值，若相邻单元中没有无效值单元，则在数组M、N中将第p行、第q列的数值赋为0，在L中将第p行、第q列的数值赋为1，并将A中的当前单元插入优先队列Q；

步骤S3、取出优先队列Q最前端的单元C，读取单元C在数组A中的行数a，列数b，对应高程z值，若矩阵S的第a行，第b列已被赋值，则直接进行步骤S8，否则从步骤S4开始进行；

步骤S4、计算单元C的局部排水方向γ；

具体内容为：所述局部排水方向γ的范围为0°-360°，且为范围内的任意角度，使用Dinf方法得到的局部排水方向，将DEM的正北方向定义为0°，并顺时针增大；

步骤S5、分别从数组M、N中读取各自第a行，第b列的数值y_c和x_c，y_c和x_c是单元C的水流汇集点相对单元C中心点的相对坐标W(x_c，y_c)，以C为中心，依次连接其周围8个相邻单元的中心点构成矩形窗口，水流从C的水流汇集点W(x_c，y_c)出发沿γ方向流动到达矩形窗口边缘，到达点记为R；

具体内容为：以C为中心，设其周边8个相邻单元，横、纵(x、y)坐标正方向，依次连接其周围8个相邻单元的中心点构成矩形窗口，以B₂、B₄、B₆、B₈单元中点为顶点的窗口，让水流从C的水流汇集点W出发沿γ方向流动到达矩形窗口边缘，到达点记为R，此时点R相对单元C中点的相对坐标(x_R，y_R)由以下方程获取：

步骤S6、获取点R所在的单元T，获取T的中心点相对C单元中心点的方向，将其记入数组S的第a行，第b列；

具体内容为：

由于相邻单元T的中心点只能位于单元C的中心点的8个方向之一，对这8个方向直接以数字1～8编号，与相邻单元的下标编码相同，以数字1～8分别表示正北、东北、正东、东南、正南、西南、正西、西北8个方向；若根据步骤S5得到的到达点R相对单元C中心点的相对坐标为(x_R，y_R)，则记入数组S的方向编号n为：

步骤S7、若单元T在数组M中对应的同行同列单元不是无效值，则计算点R相对单元T中心点的相对坐标(x₁，y₁)，分别从数组L读取C和T对应行列的单元的数值，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值；

具体内容为：

若步骤S5得到的R相对单元C中心点的相对坐标为(x_R，y_R)，步骤S6中读取的方向编号为n，则点R相对单元T中心点的相对坐标(x₁，y₁)用下面的方程获取：

x₁和y₁中最多只有1个不为0，随后从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x₀，y₀)，以集水面积值作为权重，得到点(x₀，y₀)和点(x₁，y₁)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，坐标(x，y)由以下方程获取：

上述方程中c₀和c₁分别为单元T和C的集水面积值，将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值(c₀+c₁)；

步骤S8、判别优先队列Q是否为空队列，若Q不是空队列，则返回步骤S3，否则扫描数组S中的每一个单元的方向数值，若数值不是无效值，则根据该方向值连接每个单元和下游单元的中心点，得到DEM上模拟出的水流路径。

Images