一种基于DEM和流量汇集的非弥散水流路径模拟方法转让专利
申请号 : CN201911277951.2
文献号 : CN111125893B
文献日 : 2022-11-04
发明人 : 吴鹏飞 , 刘金涛 , 姚杰夫 , 费俊源 , 刘杨洋
申请人 : 河海大学
摘要 :
权利要求 :
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列的数值yc和xc，yc和xc是单元C的水流汇集点相对单元C中心点的相对坐标W(xc，yc)，以C为中心，依次连接其周围8个相邻单元的中心点构成矩形窗口，水流从C的水流汇集点W(xc，yc)出发沿γ方向流动到达矩形窗口边缘，到达点记为R；
步骤S6、获取点R所在的单元T，获取T的中心点相对C单元中心点的方向，将其记入数组S的第a行，第b列；
步骤S7、若单元T在数组M中对应的同行同列单元不是无效值，则计算点R相对单元T中心点的相对坐标(x1，y1)，分别从数组L读取C和T对应行列的单元的数值，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x0，y0)，以集水面积值作为权重，得到点(x0，y0)和点(x1，y1)的重心，作为单元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对应单元的数值作为当前流量汇集点的相对坐标(x0，y0)，以集水面积值作为权重，得到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，其中，坐标(x，y)由以下方程得到：
其中，上述方程中c0和c1分别为单元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列的数值yc和xc，yc和xc是单元C的水流汇集点相对单元C中心点的相对坐标W(xc，yc)，以C为中心，依次连接其周围8个相邻单元的中心点构成矩形窗口，水流从C的水流汇集点W(xc，yc)出发沿γ方向流动到达矩形窗口边缘，到达点记为R；
具体内容为：以C为中心，设其周边8个相邻单元，横、纵(x、y)坐标正方向，依次连接其周围8个相邻单元的中心点构成矩形窗口，以B2、B4、B6、B8单元中点为顶点的窗口，让水流从C的水流汇集点W出发沿γ方向流动到达矩形窗口边缘，到达点记为R，此时点R相对单元C中点的相对坐标(xR，yR)由以下方程获取：步骤S6、获取点R所在的单元T，获取T的中心点相对C单元中心点的方向，将其记入数组S的第a行，第b列；
具体内容为：
由于相邻单元T的中心点只能位于单元C的中心点的8个方向之一，对这8个方向直接以数字1～8编号，与相邻单元的下标编码相同，以数字1～8分别表示正北、东北、正东、东南、正南、西南、正西、西北8个方向；若根据步骤S5得到的到达点R相对单元C中心点的相对坐标为(xR，yR)，则记入数组S的方向编号n为：步骤S7、若单元T在数组M中对应的同行同列单元不是无效值，则计算点R相对单元T中心点的相对坐标(x1，y1)，分别从数组L读取C和T对应行列的单元的数值，从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x0，y0)，以集水面积值作为权重，得到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值；
具体内容为：
若步骤S5得到的R相对单元C中心点的相对坐标为(xR，yR)，步骤S6中读取的方向编号为n，则点R相对单元T中心点的相对坐标(x1，y1)用下面的方程获取：x1和y1中最多只有1个不为0，随后从数组M、N读取单元T对应单元的数值作为当前流量汇集点的相对坐标(x0，y0)，以集水面积值作为权重，得到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，坐标(x，y)由以下方程获取：上述方程中c0和c1分别为单元T和C的集水面积值，将T在数组L中的相同行列的对应值修改为当前单元C和T的集水面积的累积值(c0+c1)；
步骤S8、判别优先队列Q是否为空队列，若Q不是空队列，则返回步骤S3，否则扫描数组S中的每一个单元的方向数值，若数值不是无效值，则根据该方向值连接每个单元和下游单元的中心点，得到DEM上模拟出的水流路径。
说明书 :
一种基于DEM和流量汇集的非弥散水流路径模拟方法
技术领域
背景技术
时间等模型参数。由于当下的水文模型大多借助栅格型数字高程模型(Digital Elevation
Model，DEM)描绘地表地形，数字化的水流路径也大多基于DEM提取得到。目前基于DEM的水
流路径提取算法包含两个大类，即弥散型和非弥散型。其中，弥散型算法认为水流是发散
的，发源于一点的水流会不断分离成多份，流向多个下游区域。受限于水流不断发散的特
点，这类算法难以提供准确的流动距离和流域边界。因此，非弥散型算法在部分水文模型中
必不可少。
一条完整的水流路径。因此，相邻单元间的局部水流路径被限制为8个方向之一。但是，真实
环境下局部水流的排水方向应该属于0°～360°之间无穷个方向之一。除了早期的D8算法和
Rho8算法只在指向8个相邻的栅格单元中心的方向中选取局部排水方向，目前已经有多种
算法可以计算得到0°～360°之间任意角度的局部排水方向，例如Tarboton(1997)发表的
Dinf算法以及Hooshyar等(2016)基于切面曲率对Dinf流向进行优化的算法。由于这些算法
大多模拟水流自单元中心起的排水方向，但排水方向又可能不指向相邻栅格单元的单元中
心，因此当模拟的水流到达相邻单元中心时，就没有了进一步指示流动轨迹的排水方向。因
此，如何借助这类局部排水方向构建完整的非弥散水流路径值得研究。
但是D8‑LTD算法只以实现起始单元的路径最优为目标。对与下游单元，仅将其水流路径作
为其上游单元的一部分，因而难以实现全局水流路径最优。
发明内容
下游沿伸，通过将上下游的水流不断汇集到一处，从上游到下游不断确定每个单元的下游
单元，并将局部排水方向与下游单元中心点方向间的偏差传递至下游，在下游的水流聚集
点的位置确定过程中考虑该偏差，发明了一种基于DEM和流量汇集的非弥散水流路径模拟
方法。
标X，纵坐标Y的坐标值，L用于记录每个单元的上游集水面积，S用于记录下游单元的中心点
方向，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，同时创建优先队列Q，Q能存储栅格单元的信息，在数组A中行号、列号、高程，并能根据栅格的高
程将其从高到低排序，遇到栅格高程相同，则将后加入Q的单元排在先加入Q的单元的后方；
的行列数，第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；
元的中心点构成矩形窗口，水流从C的水流汇集点W(xc，yc)出发沿γ方向流动到达矩形窗口
边缘，到达点记为R；
单元T对应单元的数值作为当前流量汇集点的相对坐标(x0，y0)，以集水面积值作为权重，得
到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数
组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单
元C和T的集水面积的累积值；
游单元的中心点，可得到DEM上模拟出的水流路径。
坐标(x0，y0)，以集水面积值作为权重，得到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，
标X，纵坐标Y的坐标值，L用于记录每个单元的上游集水面积，S用于记录下游单元的中心点
方向，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，同时创建优先队列Q，Q能存储栅格单元的信息，在数组A中行号、列号、高程，并能根据栅格的高
程将其从高到低排序，遇到栅格高程相同，则将后加入Q的单元排在先加入Q的单元的后方；
值，该固定值被定义为无效值，无效值被定义为‑9999；
的行列数，第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；
元的中心点构成矩形窗口，水流从C的水流汇集点W(xc，yc)出发沿γ方向流动到达矩形窗口
边缘，到达点记为R；
发沿γ方向(如图3中(a)、3中(b)中黑色箭头方向)流动到达矩形窗口边缘，到达点记为R
(如图3中(a)中的R1，图3中(b)中的R2)，此时点R相对单元C中点的相对坐标(xR，yR)由以下
方程获取：
北、东北、正东、东南、正南、西南、正西、西北8个方向；若根据步骤S5得到的到达点R相对单元C中心点的相对坐标为(xR，yR)，则记入数组S的方向编号n为：
单元T对应单元的数值作为当前流量汇集点的相对坐标(x0，y0)，以集水面积值作为权重，得
到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数
组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单
元C和T的集水面积的累积值；
(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，坐标(x，y)由以下方程获取：
游单元的中心点，可得到DEM上模拟出的水流路径。
坐标(x0，y0)，以集水面积值作为权重，得到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，
局部排水方向确定下游单元，通过依次连接上下游单元的栅格中心得到最终的水流路径。
相较于传统的水流路径模拟方法，本发明的方法具有所确定水流路径更为准确的优势。因
此，本发明为基于DEM的水文模型提供了高精度的水流过程模拟方法。
附图说明
具体实施方式
限定本发明。
为对本发明创造的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明创造的描述中，除非另
有说明，“多个”的含义是两个或两个以上。
可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上
述术语在本发明创造中的具体含义。
下：
标X，纵坐标Y的坐标值，L用于记录每个单元的上游集水面积，S用于记录下游单元的中心点
方向，将M、N、L、S中的第一列、最后一列、第一行、最后一行的所有数值都赋为无效值，同时创建优先队列Q，Q能存储栅格单元的信息，在数组A中行号、列号、高程，并能根据栅格的高
程将其从高到低排序，遇到栅格高程相同，则将后加入Q的单元排在先加入Q的单元的后方；
值，该固定值被定义为无效值，无效值被定义为‑9999；
的行列数，第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；
元的中心点构成矩形窗口，水流从C的水流汇集点W(xc，yc)出发沿γ方向流动到达矩形窗口
边缘，到达点记为R；
出发沿γ方向(如图3中(a)、3中(b)中黑色箭头方向)流动到达矩形窗口边缘，到达点记为R
(如图3中(a)中的R1，图3中(b)中的R2)，此时点R相对单元C中点的相对坐标(xR，yR)由以下
方程获取：
北、东北、正东、东南、正南、西南、正西、西北8个方向；若根据步骤S5得到的到达点R相对单元C中心点的相对坐标为(xR，yR)，则记入数组S的方向编号n为：
单元T对应单元的数值作为当前流量汇集点的相对坐标(x0，y0)，以集水面积值作为权重，得
到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，并将T在数
组M、N中的对应值分别修改为x和y，随后将T在数组L中的相同行列的对应值修改为当前单
元C和T的集水面积的累积值；
(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，坐标(x，y)由以下方程获取：
游单元的中心点，可得到DEM上模拟出的水流路径。
坐标(x0，y0)，以集水面积值作为权重，得到点(x0，y0)和点(x1，y1)的重心，作为单元T的新流量汇集点的相对坐标(x，y)，
边8个相邻单元是否存在无效值单元，此处发现A中单元B4与无效值单元A5相邻，则依照当
前单元在数组A中的行列数，即第B行，第4列，在数组M、N、L、S中均将第B行，第4列的数值赋为无效值‑9999，同时将A中除第一列、最后一列、第一行、最后一行、以及B4外的其他单元加入Q(图6中(b))；
出发沿γ方向流动到达矩形窗口边缘点R，由方程(3)和(4)得到此时点R相对单元C中点的
相对坐标(xR，yR)为(‑0.33，‑1)；
行、第3列(图6中(g))；
的集水面积值均为1，并从数组M、N读取单元D3对应单元的数值作为当前流量汇集点的相对
坐标(x0，y0)，即(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；
值‑9999，则根据该方向值连接每个单元和下游单元的中心点，即可得到DEM上模拟出的水
流路径，如图6中(f)。
散水流路径，可视化的水流路径如图7中(c)‑7中(h)所示，图7中的黑色线条为水流路径，灰
色线条为等高线，理想的水流路径应该垂直于等高线流动，可以发现D8方法得到的水流路
径非常直且不垂直于等高线(图7中(c)，图7中(d))，虽然D8‑LTD在平面上应用效果良好，水
流路径能大致垂直于等高线(图7中(e))，但是该方法在凹面上得到的水流路径出现了紊乱
交叉的现象(图7中(f))，本发明成功弥补了已有方法的不足，在两个DEM面上展现了良好的
应用效果(图7中(g)、图7中(h))。