基于Priority‑Flood的内流流域提取方法

摘要

本发明公开了一种基于Priority‑Flood的内流流域提取方法，采用先漫水示踪法，通过模拟水流注满洼地并溢出的过程，获取洼地内正确的汇流方向并同时获得水溢出每个洼地时的“出水口点”以及该洼地的面积；然后根据洼地面积大小筛选符合要求的洼地；根据流向数据，将流向出水口点的所有栅格归并为一个区域，即完成内流流域的划分与提取。本发明适用于基于DEM进行内流区的流域提取，无需填充洼地，且只需对所有栅格遍历一次，能够自动化、高效的提取内流流域，为下一步内流流域的地形地貌和水文分析奠定重要的基础。

说明书

技术领域

本发明涉及一种内流流域提取方法，特别是一种基于Priority-Flood的内流流域提取方法。

背景技术

内流区指的是它的集水区与外部海洋环境没有水文联系的流域，内流区域占地球表面的20％，分布范围广泛，是地球地貌、地表生态环境的重要组成部分。针对内流流域自身的地形、水文特点，进行流域单元的划分和提取，是研究内流区及内流流域的水文过程、地貌发育、生物演化等特征重要基础。

基于数字高程模型(DEM)的水文分析方法是目前流域划分提取和研究分析的常用方法。其中流向模型算法是首要步骤且是最重要的基础部分，目前有单流向算法和多流向算法。D8(Deterministic eight-neighbors)作为典型的单流向算法，将水流方向可能性抽象为平面上的八个方向，限制了对水流实际上无穷多种运动方向的描述精度。后续也提出来诸如Rho8、Rho4、DEMON和Dinf等诸多改进方法，但单流向算法将所有水量汇入统一下游单元的设计必然导致它不适合模拟在地势平坦的坡面上水流散漫流动的情况，在平面区域会出现大量不合理的平行河网。为解决该问题，人们提出来多流向算法，可以使中心格网单元有多个流出单元，并按一定的权重分配水量。但合理参数值的选取成为难以确定的问题。不管基于何种流向算法，基于DEM的流域划分提取方法通常包括三个步骤：计算流向—提取河网—划分流域。

在现有的流向算法中，计算流向的必要前提步骤是填洼，因为洼地的存在会令水流无法流出，导致接下来的河网及流域提取失败。但是针对于内流区，此种方法显然是错误的。内流区本身就是作为洼地存在，这也是内流区与外界没有水文联系的原因。按常规的方法填充洼地之后，丧失了整个内流区的地形地貌信息，整个内流区作为一个平地存在，后续提取汇流网络时会产生大量的平行河网，所划分和提取内流流域的结果也是错误的。这也就导致了现有的方法无法用于内流区流域的划分与提取。

发明内容

本发明所要解决的技术问题是提供一种基于Priority-Flood的内流流域提取方法，它能够自动化、高效的提取内流流域。

为解决上述技术问题，本发明所采用的技术方案是：

一种基于Priority-Flood的内流流域提取方法，其特征在于包含以下步骤：

步骤一：定义存储数据及判断栅格访问顺序的三个队列及流向栅格：包括平地队列、洼地队列、坡地队列，以及用来存储流向结果的栅格数据；

步骤二：将所有边界栅格压入坡地队列中：首先选取DEM数据的边界栅格，然后根据制定的边界栅格流向规则计算每一个栅格的流向，并将流向栅格对应位置赋以相应的流向编码，然后将边界栅格标记为已处理，最后将所有栅格压入坡地队列中；

步骤三：循环处理三个队列中的元素，直到三个队列都为空：根据元素选取顺序决策规则，从三个队列中选择弹出一个元素进行处理：循环处理弹出栅格所有未处理的邻域栅格；

步骤四：根据高程判断待处理栅格所属类型：通过比较待处理栅格高程与溢出高程的大小，确定邻域栅格所属的队列；

步骤五：根据待处理栅格类型选择不同的流向计算方法：通过比较待处理栅格高程与中心栅格高程的大小，根据不同的情况选择相应的流向计算方法，然后将流向栅格对应位置赋以相应的流向编码；

步骤六：判断待处理栅格点是否为洼地出水口点：确立洼地出水口的特征，判断待处理栅格是否为洼地出水口点，若是，标记此栅格点；更新洼地面积：判断待处理栅格是否为洼地，若是，更新其所属洼地面积大小，然后将其压入所属的队列中，并将待处理栅格标记为已处理；

步骤七：筛选洼地，归并流域：选择合适的洼地面积阈值，筛选出符合要求的洼地及其洼地出水口点，然后递归处理将所有流向出水口点的栅格标记为同一流域。

进一步地，所述步骤一中三个栅格要素队列包括平地队列为普通的先进先出队列、洼地队列也为普通的先进先出队列和坡地队列为最小优先队列。

进一步地，所述步骤二中，边界栅格流向规则为：按照D8编码顺序遍历该边界栅格的邻域栅格，找到第一个在边界外的或者值为NoData的栅格，边界栅格即流向该栅格。

进一步地，所述步骤三中，元素选取决策顺序规则为：首先判断洼地队列是否为空，若不为空，弹出洼地队列队首元素作为中心栅格；否则，继续判断平地队列是否为空，若平地队列不为空，则弹出平地队列队首元素作为中心栅格；否则，则继续判断坡地队列是否为空，若坡地队列不为空，弹出坡地队列队首元素作为中心栅格，并将溢出高程设为该中心栅格高程值。

进一步地，所述步骤四中，根据高程判断待处理栅格所属类型的方法为：若待处理栅格高程小于溢出高程，则待处理栅格属于洼地；若待处理栅格高程等于溢出高程，则待处理栅格属于平地；若待处理栅格高程大于溢出高程，则待处理栅格属于坡地。

进一步地，所述步骤五中，选择不同的流向计算方法为：若待处理栅格高程值等于中心栅格高程值，则待处理栅格流向中心栅格；若待处理栅格高程值大于中心栅格，则找出待处理栅格邻域中所有已处理栅格，在这些栅格中找到高程值比待处理栅格高程值小且坡降最大的栅格，流向该栅格；若待处理栅格小于中心栅格，则找出待处理栅格邻域中所有已处理栅格，在这些栅格中找到高程值比待处理栅格高程值大且坡降最小的栅格，流向该栅格。

进一步地，所述步骤六中，判断待处理栅格点是否为洼地出水口点和更新洼地面积的方法为：当洼地队列为空时，第一个比溢出高程小的点即为新的洼地出水口点；更新洼地面积的方法：根据第四步中判断的待处理栅格所属的类型，若待处理栅格点属于洼地且不是洼地出水口点，则将当前洼地面积大小加一，并将其压入洼地队列中；若属于平地，则将其压入平地队列；若属于坡地，则将其压入坡地队列中。

进一步地，所述步骤七中，筛选洼地及归并流域的方法为：通过合适的面积阈值筛选洼地并寻找出水口点，并将流向出水口点的栅格压入队列中，弹出队列中元素，将流向该元素的栅格继续压入队列中，循环处理，直到队列为空。

本发明与现有技术相比，具有以下优点和效果：

1、本发明采用先漫水示踪法，通过模拟水流注满洼地并溢出的过程，获取洼地内正确的汇流方向并同时获得水溢出每个洼地时的“出水口点”以及该洼地的面积；然后根据洼地面积大小筛选符合要求的洼地；根据流向数据，将流向出水口点的所有栅格归并为一个区域，即完成内流流域的划分与提取。

2、本发明适用于基于DEM进行内流区的流域提取，无需填充洼地，且只需对所有栅格遍历一次，能够自动化、高效的提取内流流域，为下一步内流流域的地形地貌和水文分析奠定重要的基础。

3、改变了传统利用DEM计算流向时必须进行洼地填充的方法，使得DEM水文分析扩展到内流区域，打破了传统方法只能用于外流区域的局限；为进一步分析内流区流域的水文、地貌特征等奠定了基础。

附图说明

图1是本发明的基于Priority-Flood的内流流域提取方法的流程框图。

图2是本发明的流向算法流程图。

图3是本发明的洼地筛选及流向归并流程图。

图4是本发明的单流向D8流向编码示意图。

图5是本发明的优先漫水算法(Priority-Flood)的示意图。

图6是本发明的栅格类型示意图(深色栅格表示边界，外围浅色栅格表示NoData)。

图7是本发明的栅格流向计算示意图。

图8是本发明的整体流向示意图。

图9是本发明的洼地及内流流域示意图。

图10是本发明的实施例的青藏高原内流区内流流域划分图。

具体实施方式

下面结合附图并通过实施例对本发明作进一步的详细说明，以下实施例是对本发明的解释而本发明并不局限于以下实施例。

如图所示，本发明的基于Priority-Flood的内流流域提取方法，包含以下步骤：

步骤一：各种变量的定义及初始化。首先，读取原始DEM数据，定义一个与原始DEM大小相同的栅格存储流向数据。其次，定义存储栅格单元元数据的结构体，包括栅格行号、栅格列号及栅格高程值。该结构体即为之后三个队列中所要存储的数据。接着定义平地队列，用来存储标记为平地类型的栅格；定义坡地队列，用来存储标记为坡地类型的栅格；定义洼地队列，用来存储标记为洼地类型的栅格。其中，平地队列和洼地队列为普通的先进先出队列，坡地队列为最小优先队列(即属性值最小的元素拥有最高的优先级，会最优先被弹出处理)。接着定义洼地出水口标识为False，定义溢出高程为0。

步骤二：将所有边界栅格压入坡地队列中。首先，循环遍历所有栅格，对于所有的边界栅格(即栅格有效数据范围的边界，如果一个栅格的邻域栅格中存在空值或无数据，即判断该栅格为边界栅格)，按照边界栅格流向算法，计算该栅格的流向，修改对应位置的流向栅格的值。然后将该栅格压入坡地队列中。以上流向均采用D8流向编码。计算边界栅格流向的算法为：按照D8编码顺序遍历一个边界栅格的邻域，栅格流向第一个处于数据范围外或值为NoData的栅格。D8算法计算坡降的计算公式如下：

式中，Slp_c为待求中心栅格的最大等高线加权高差邻域单元的坡度，h_c为中心栅格单元的高程值，h_i为第i号邻域单元的高程值，L_i为第i号邻域单元的等高线加权值。Li的定义如下：

步骤三：处理三个队列元素。首先，判断三个队列是否都为空，若是，则流向计算结束，跳转到第七步；否则继续完成本步计算和接下来第四、五、六步。依次按照洼地栅格、平地栅格和坡地栅格的顺序，找到第一个不为空的队列，并将弹出队列中元素，将其称之为中心栅格。如果弹出的是坡地队列，还要设置溢出高程为中心栅格高程。如果洼地队列为空，需将洼地出水口标识设置为False，以便开始标记下一块洼地出水口。

步骤四：判断中心栅格邻域栅格所属类型。在第三步确定中心高程后，遍历处理每一个邻域栅格(称为待处理栅格)，若待处理栅格高程小于溢出高程，则待处理栅格属于洼地，在完成第五、六步计算后，将其压入洼地队列；若待处理栅格高程等于溢出高程，则待处理栅格属于平地，在完成第五、六步计算后，将其压入平地队列；若待处理栅格高程大于溢出高程，则待处理栅格属于坡地，在完成第五、六步处理后，将其压入坡地队列。

步骤五：计算待处理栅格的流向。若待处理栅格高程等于中心栅格高程，则待处理栅格流向中心栅格即可；若待处理栅格高程小于中心栅格高程，则先找到待处理栅格邻域中已处理的栅格，再从这些栅格中找到高程值大于待处理栅格且坡降最小的栅格，流向该栅格即可；若待处理栅格高程大于中心栅格高程，则先找到待处理栅格邻域中已处理的栅格，再从这些栅格中找到高程值小于待处理栅格且坡降最大的栅格，流向该栅格即可。

步骤六：判断出水口点并更新洼地面积。若待处理栅格点的高程小于中心栅格点，且洼地标识为False，则将该待处理栅格点标记为洼地出水口点，修改洼地标识为True，同时新建一条洼地记录，包括洼地出水口点的坐标位置、高程信息和洼地面积大小。若待处理栅格高程小于中心栅格高程但不是洼地出水口点，则该待处理栅格点属于洼地，将上一条洼地记录的面积加一。

步骤七：筛选洼地及归并流域。根据实地情况确定合适洼地面积阈值，选择符合要求的洼地，生成一个洼地列表。从列表中取出一块洼地，获得出水口点栅格，将其加入一个临时队列中，并赋以一个标记值(代表某一流域)。从临时队列中取出一个栅格元素，判断该栅格所有邻域栅格，如果邻域栅格流向该栅格且邻域栅格不属于洼地列表中洼地的出水口点，则将该邻域栅格标记为同一流域，并加入临时队列中。循环处理临时队列，直到临时队列为空。

下面以青藏高原的内流区(面积)为例对本发明进行说明。

1)获取数据。本例中使用分辨率90米的SRTM数据。根据青藏高原内流区范围截取相应的DEM数据。

2)队列及变量初始化。读入DEM数据，新建与DEM大小相等的流向数据，初始化相关的队列、变量等。

3)数据边界压入队列。识别数据边界，将边界存储到优先队列中。

4)队列数据处理。根据选择顺序选取不同的队列元素，处理元素的邻域栅格，计算各栅格的流向，判断洼地出水口点，计算洼地面积，并将其压入所属的队列中，循环处理直到三个队列都为空。最后，获得栅格流向，识别除青藏高原内流区内的所有洼地，以及每一块洼地的面积大小和相应的出水口点位置。

5)确定面积阈值并筛选。根据青藏高原样区的实际情况及研究目的，选择合适的洼地面积阈值，并根据阈值的大小筛选识别出来的洼地。

6)归并流域。根据所识别出来的洼地出水口的位置，将所有汇入该点的流向栅格归并为同一流域，合并小流域。最后得到青藏高原内流流域划分提取结果。

本说明书中所描述的以上内容仅仅是对本发明所作的举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种修改或补充或采用类似的方式替代，只要不偏离本发明说明书的内容或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。