一种考虑山区平原地貌差异的子流域单元划分方法

摘要

本发明公开了一种考虑山区平原地貌差异的子流域单元划分方法，包括计算流域内各栅格流向和汇流累积量，基于河网阈值提取模拟河网；设定平原单元长度，对平原区及其周边栅格进行N*N大小聚合，对每片N*N区域设置平原单元标记序号；从流域出口开始，以模拟河网栅格为主干，向上溯源遍历，根据遇到的栅格类型进行计算单元划分，对划分后的单元赋值唯一的子流域编号；确定相邻平原单元和山区单元之间的流向；根据流向关系，从河流出口单元开始溯源遍历，生成子流域上下游拓扑关系表。优点是：更好、更全面地刻画流域内不同地貌的地表、地下水流动机制，提高分布式水文模型的水文模拟精度，为分布式水文模型在大尺度流域的构建和应用提供理论和实践支撑。

说明书

技术领域

本发明涉及子流域单元划分技术领域，尤其涉及一种考虑山区平原地貌差异的子流域单元划分方法。

背景技术

在大尺度流域分布式水文模型构建应用中，采用子流域进行计算单元划分是常用的方法。相比于纯栅格单元划分，该方法具有能够保持各计算单元间的水文过程相对独立，且减少计算单元数量的优点。子流域划分方法假定地表汇流过程和地下汇流过程有着相同的分水岭，且其内地表、地下产流由分水岭向河道汇集，最终流出子流域。由于山区地表地下分水岭明显，因此比较适用于山区。但是，对平原区而言，由于地势平坦，各单元之间无明显的分水岭划分边界以及恒定的地表地下水流向，表现为地表水系散乱，地下水流动方向不固定，而且地表和地下各单元之间也存在比较活跃的水量交换过程，因此平原区是不适合采用子流域进行计算单元划分。对于大尺度流域而言，流域内既存在山区又存在平原区，如果只采用单一的子流域对整个流域进行计算单元划分，将无法准确描述流域内平原区水循环过程；如果只采用栅格单元进行划，则计算单元数量会变得非常多，会导致计算维度灾。

发明内容

本发明的目的在于提供一种考虑山区平原地貌差异的子流域单元划分方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明采用的技术方案如下：

一种考虑山区平原地貌差异的子流域单元划分方法，包括如下步骤，

S1、以DEM数据为基础，计算流域内各栅格流向和汇流累积量，并基于设定的河网阈值提取模拟河网；

S2、设定平原单元长度，对平原区及其周边栅格进行N*N大小聚合，并对每片N*N区域设置平原单元标记序号；

S3、从流域出口开始，以模拟河网栅格为主干，向上溯源遍历，根据遇到的栅格类型进行计算单元划分，并对划分后的各单元赋值一个唯一的子流域编号，该编号为从1开始的自然数；划分出的单元包括平原单元和山区单元；

S4、利用栅格流向信息，确定相邻平原单元和山区单元之间的流向；

S5、根据平原单元和山区单元之间的流向关系，从河流出口单元开始溯源遍历，生成子流域上下游拓扑关系表。

优选的，栅格流向采用D8算法计算得出，即将该栅格流向周边8个栅格中坡度最陡的栅格的流向作为该栅格的流向；栅格的汇流累积量即为上游所有流入当前栅格的栅格数量之和。

优选的，河网阈值用于定义各个栅格是属于河网栅格或坡面栅格，当栅格的汇流累积量大于河网阈值，则将该栅格定义为河网栅格，否则将该栅格定义为坡面栅格；所有的所述河网栅格构成模拟河网，所述河网阈值须保证提取的模拟河网在源头区与实际河网源头一致。

优选的，步骤S2具体包括如下内容，

S21、将平原单元长度设定为栅格长度的整数倍，以确保划分的平原单元和山区单元边界吻合；

S22、从栅格图层起始行列开始，逐行逐列按照步长进行累加遍历，每次对N行N列共N*N个相邻栅格进行检查，如果该区域内的平原栅格的个数大于或等于山区栅格的个数，则将该N*N区域标记为平原单元，并对区域内所有栅格设置相同的平原单元标记序号，所述序号为从1开始的自然数；若平原栅格的个数小于山区栅格的个数，则不做处理。

优选的，所述平原单元为N*N栅格型单元，山区单元为根据汇流路径确定的子流域型单元，所述平原单元包括单元内有河道的平原河道单元和单元内无河道的平原坡面单元；所述山区单元包括单元内有河道的山区河道单元和单元内无河道的山区坡面单元。

优选的，步骤S3中，对溯源遍历的栅格根据其属性分别使用不同的编码处理模式；若遍历起始栅格为山区河道栅格，则按照山区河道栅格编码模式处理；若为平原河道栅格，则按照平原栅格编码模式进行处理；直到河流源头结束遍历；

所述平原栅格编码模式为，

查询当前平原栅格的平原单元标记序号，将该平原单元周围(2N-1)*(2N-1)范围内平原单元标记序号等于该栅格平原单元标记序号的栅格赋值相同的子流域编号；若该平原单元中有河道栅格，则标记为平原河道单元，否则标记为平原坡面单元；

以流入该平原单元的栅格为起点，重新进行遍历搜索，若流入的栅格为平原栅格，则按照平原栅格编码模式进行处理；若流入的栅格是山区河道栅格，则按照山区河道栅格编码模式进行处理；若流入的栅格为山区坡面栅格，则按照山区坡面栅格编码模式进行处理；

所述山区河道栅格编码模式为，

对河网栅格进行溯源遍历，直到遇到人为分割点、分岔河道或平原单元标记栅格为止，将其遍历的河道栅格作为一个河段，并赋值一个子流域编号；以该河段内的河道栅格为起点，溯源追迹坡面栅格直到无上游流入栅格或者上游流入栅格为平原单元标记栅格为止，将所有汇入的坡面栅格赋值该河道栅格对应的子流域编号，所有的河道栅格和坡面栅格构成一个山区单元，并将该单元类型标记为山区河道单元；

若在山区河道单元上游有平原栅格流入，则以各流入平原栅格为起点按照平原栅格编码模式进行处理；

所述山区坡面栅格编码模式为，

只有当有山区坡面栅格流入平原单元栅格时才存在山区坡面栅格编码模式，此时需要查找到流入同一个平原单元的所有山区坡面栅格，并以此为起点溯源遍历上游流入的山区坡面栅格，直到无上游流入栅格或者上游流入栅格为平原单元标记栅格为止，将所遍历的栅格都划归为同一个山区单元，其子流域编号等于流入平原单元子流域编号加1，并将该单元类型标记为山区坡面单元。

若山区坡面栅格溯源遍历过程中，遇到上游流入栅格为平原单元而终止的情况，需要对遇到的平原栅格采用平原栅格编码模式进行处理。

优选的，步骤S4中，河道类型的单元按照河道栅格的流向确定单元间上下游关系；山区坡面单元的流向为出口栅格流向的栅格所在的平原单元；平原坡面单元按照平原坡面单元流向确定规则加以确定；

所述平原坡面单元流向确定规则为，确定平原坡面单元的高程及其周边单元的高程，利用平原坡面单元的高程及周边单元的高程计算平原坡面单元到周边单元的比降，则该平原坡面单元的流向是流向比降最陡的那个周边单元；所述比降即平原坡面单元的高程减去周边单元的高程后除以两者之间的距离。

优选的，所述平原单元的高程等于单元内所有栅格高程的平均值；所述山区单元的高程为山区单元内与该平原坡面单元直接相邻的各个栅格的高程的平均值。

优选的，平原单元和周边平原单元之间的距离采用形心点直线距离表示，上下左右相邻单元距离为1个平原单元长度，对角线相邻单元距离为

优选的，生成的子流域上下游拓扑关系表中包括当前子流域编码，下游子流域编码，上游子流域编码，河道坡面单元标记，平原山区单元标记；子流域编号等于0表示没有上游或下游子流域；河道坡面单元标记等于0表示坡面，等于1表示河道；平原山区单元标记等于1表示坡面，等于0表示山区。

本发明的有益效果是：1、以DEM为基础，考虑山区平原区地貌差异，进行子流域划分和编码，能够保证划分得到的子流域基本符合存在平原区的流域内各单元间的上下游关系，且能根据需求将平原区内的小栅格聚合为一个大的平原单元进行子流域编码，保证子流域范围划分的正确性。2、能够更好、更全面地刻画流域内不同地貌的地表、地下水流动机制，提高分布式水文模型的水文模拟精度，为分布式水文模型在大尺度流域的构建和应用提供理论和实践支撑。

附图说明

图1是本发明实施例中子流域划分方法的流程示意图；

图2是本发明实施例中平原区的待划分流域的示意图；

图3是本发明实施例中平原范围及平原标记序号的示意图；

图4是本发明实施例中平原栅格编码模式的流程示意图；

图5是本发明实施例中山区河道栅格编码模式的流程示意图；

图6是本发明实施例中山区坡面栅格编码模式的流程示意图；

图7是本发明实施例中单元间流向示意图；

图8是本发明实施例中子流域划分结果示意图；

图9是本发明实施例中子流域上下游拓扑关系表示意图。

具体实施方式

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

如图1所示，本发明提供了一种考虑山区平原地貌差异的子流域单元划分方法，包括如下步骤，

S1、以DEM数据为基础，计算流域内各栅格流向和汇流累积量，并基于设定的河网阈值提取模拟河网；

S2、设定平原单元长度，对平原区及其周边栅格进行N*N大小聚合，并对每片N*N区域设置平原单元标记序号；

S3、从流域出口开始，以模拟河网栅格为主干，向上溯源遍历，根据遇到的栅格类型进行计算单元划分，并对划分后的各单元赋值一个唯一的子流域编号，该编号为从1开始的自然数；划分出的单元包括平原单元和山区单元；

S4、利用栅格流向信息，确定相邻平原单元和山区单元之间的流向；

S5、根据平原单元和山区单元之间的流向关系，从河流出口单元开始溯源遍历，生成子流域上下游拓扑关系表。

本实施例中，本发明的子流域划分方法具体包括以上五个步骤，下面分别针对这五个步骤的具体内容进行解释。

一、提取模拟河网

该部分即为步骤S1的内容，步骤S1的具体内容为以DEM数据为基础，采用D8算法计算流域内各栅格的流向，并基于流向计算各栅格汇流累积数设定河道阈值，将汇流累积数大于阈值的栅格设为河网栅格，从而提取出模拟河网。

采用D8算法计算流域内各栅格的流向的具体过程为，将该栅格流向周边8个栅格中坡度最陡的栅格的流向作为该栅格的流向；栅格的汇流累积量计算过程为，上游所有流入当前栅格的栅格数量之和。

本实施例中，河网阈值用于定义各个栅格是属于河网栅格或坡面栅格，当栅格的汇流累积量大于河网阈值，则将该栅格定义为河网栅格，否则将该栅格定义为坡面栅格；所有的所述河网栅格构成模拟河网。

河网阈值的具体取值由用户进行确定，所述河网阈值须保证提取的模拟河网在源头区与实际河网源头一致；河网阈值的取值过程为，通过实验从小逐渐加大，直到提取的模拟河网在源头区同实际河网源头相一致。

二、设置平原单元标记序号

该部分即为步骤S2的内容，步骤S2包括如下两个步骤，如图2和图3所示；

S21、将平原单元长度设定为栅格长度的整数倍，以确保划分的平原单元和山区单元边界吻合；

平原单元长度为DEM栅格长度的整数倍(N倍)，确保划分的平原单元和山区子流域单元边界相吻合；单元长度越大，平原单元个数越少，单元长度越小，平原单元个数就越多。平原单元长度值具体由用户研究确定，一般推荐为3-10个DEM栅格长度。

S22、从栅格图层起始行列开始，逐行逐列按照步长进行累加遍历，从(1,1)栅格开始检索，其下一行检索栅格为(N+1,1)、下一列检索栅格为(1，N+1)，依次类推；每次对N行N列共N*N个相邻栅格进行检查，如果该区域内的平原栅格的个数大于或等于山区栅格的个数，则将该N*N区域标记为平原单元，并对区域内所有栅格设置相同的平原单元标记序号，所述序号为从1开始的自然数；若平原栅格的个数小于山区栅格的个数，则不做处理。

步骤S2处理结束后，能够得到一个平原单元标记图层，图层中被认定为平原单元的栅格具有大于0的序号，被认定为山区的则无数据。该平原单元标记图层主要用于确定哪些栅格属于相同的平原单元。

具体的，步骤S2中可以设定平原单元长度为5个Dem栅格长度,从(1,1)栅格开始按5个栅格步长进行逐行逐列遍历，对平原区单元进行标记。如果被检索的5*5范围内的平原栅格个数大于等于山区栅格栅格个数，则将这25个栅格都赋值一个平原单元标记序号，如果小于山区栅格个数，则赋值为0。平原单元标记序号为从1开始依次增加的自然数。

三、对各个单元进行子流域编号

该步骤对应步骤S3，步骤S3中所述平原单元为步骤S2确定的N*N栅格型单元，山区单元为根据汇流路径确定的子流域型单元，如果单元内有河道则成为河道型单元，否侧成为坡面型单元。平原单元和山区单元交叉共4个小类别单元；具体的，所述平原单元包括单元内有河道的平原河道单元和单元内无河道的平原坡面单元；所述山区单元包括单元内有河道的山区河道单元和单元内无河道的山区坡面单元。

如图4至图6所示，本实施例中，步骤S3的具体执行过程为，对溯源遍历的栅格根据其属性分别使用不同的编码处理模式；若遍历起始栅格为山区河道栅格，则按照山区河道栅格编码模式处理；若为平原河道栅格，则按照平原栅格编码模式进行处理；直到河流源头结束遍历；

所述平原栅格编码模式为，

查询当前平原栅格的平原单元标记序号，将该平原单元周围(2N-1)*(2N-1)范围内平原单元标记序号等于该栅格平原单元标记序号的栅格赋值相同的子流域编号；若该平原单元中有河道栅格，则标记为平原河道单元，否则标记为平原坡面单元；

以流入该平原单元的栅格为起点，重新进行遍历搜索，若流入的栅格为平原栅格，则按照平原栅格编码模式进行处理；若流入的栅格是山区河道栅格，则按照山区河道栅格编码模式进行处理；若流入的栅格为山区坡面栅格，则按照山区坡面栅格编码模式进行处理；

所述山区河道栅格编码模式为，

对河网栅格进行溯源遍历，直到遇到人为分割点(水库、水文站位置点等)、分岔河道或平原单元标记栅格为止，将其遍历的河道栅格作为一个河段，并赋值一个子流域编号；以该河段内的河道栅格为起点，溯源追迹坡面栅格直到无上游流入栅格或者上游流入栅格为平原单元标记栅格为止，将所有汇入的坡面栅格赋值该河道栅格对应的子流域编号，所有的河道栅格和坡面栅格构成一个山区单元，并将该单元类型标记为山区河道单元；

若在山区河道单元上游有平原栅格流入，则以各流入平原栅格为起点按照平原栅格编码模式进行处理；

所述山区坡面栅格编码模式为，

只有当有山区坡面栅格流入平原单元栅格时才存在山区坡面栅格编码模式，此时需要查找到流入同一个平原单元的所有山区坡面栅格，并以此为起点溯源遍历上游流入的山区坡面栅格，直到无上游流入栅格或者上游流入栅格为平原单元标记栅格为止，将所遍历的栅格都划归为同一个山区单元，其子流域编号等于流入平原单元子流域编号加1，并将该单元类型标记为山区坡面单元。

若山区坡面栅格溯源遍历过程中，遇到上游流入栅格为平原单元而终止的情况，需要对遇到的平原栅格采用平原栅格编码模式进行处理。

具体的，步骤S3可以按照如下过程执行，从流域出水口开始，以模拟河网栅格为主干向上溯源遍历，根据遍历栅格属性采用不同的编码处理模式进行计算单元划分。

如果遍历开始的栅格是山区河道栅格，调用山区河道栅格编码模式，以该栅格为起点继续向上溯源遍历查找河网栅格，直至遇到人为分割点、分岔河道、平原单元标记栅格结束，将其遍历的河道栅格作为一个河段，赋值一个子流域编号，然后以每个河道栅格为起点，溯源追迹坡面栅格直至无上游流入栅格或者上游流入栅格为平原单元标记栅格结束，将所有坡面栅格赋值该子流域编号，并将该单元类型标记为山区河道单元，结束一个河段后，开始下一河段的溯源查找。

如果遍历开始的栅格是平原单元标记栅格，则调用平原栅格编码模式，将周围(2N-1)*(2N-1)范围内平原单元标记序号等于该栅格平原单元标记序号的栅格赋值相同的子流域编号(该编号为前面子流域编号+1)，如果该平原单元中有河道栅格，则标记为平原河道单元，否则标记为平原坡面单元；然后以流入该平原单元的栅格为起点，重新遍历搜索，如果流入栅格是平原栅格，则递归调用平原栅格编码模式进行处理；如果流入栅格是山区河道栅格，则递归调用山区河道栅格编码模式进行处理；如果流入栅格是山区坡面栅格，则递归调用山区坡面栅格编码模式进行处理。

四、确定相邻平原单元和山区单元之间的流向

该步骤对应步骤S4，如图7所示，步骤S4中河道类型的单元按照河道栅格的流向确定单元间上下游关系；山区坡面单元的流向为出口栅格流向的栅格所在的平原单元；平原坡面单元按照平原坡面单元流向确定规则加以确定；

所述平原坡面单元流向确定规则为，确定平原坡面单元的高程及其周边单元的高程，利用平原坡面单元的高程及周边单元的高程计算平原坡面单元到周边单元的比降，则该平原坡面单元的流向是流向比降最陡的那个周边单元；所述比降即平原坡面单元的高程减去周边单元的高程后除以两者之间的距离。

本实施例中，所述平原单元的高程等于单元内所有栅格高程的平均值；所述山区单元的高程为山区单元内与该平原坡面单元直接相邻的各个栅格的高程的平均值。

本实施例中，平原单元和周边平原单元之间的距离采用形心点直线距离表示，上下左右相邻单元距离为1个平原单元长度，对角线相邻单元距离为

具体的，步骤S4可以按照如下过程执行，依据单元类型确定所有单元间的流向。河道类型的单元依据单元内河网水系拓扑关系确定单元间的流向；山区坡面单元则流向平原栅格编码模式中确定的平原单元；对于平原坡面单元，先确定平原坡面单元及其周边单元的高程，其次计算平原坡面单元到周边计算单元的比降，最后确定平原坡面单元的流向为从平原坡面单元到比降最低单元的方向，其中平原单元高程等于单元内的DEM高程平均值，平原单元相邻的山区单元高程等于山区单元内和该平原坡面单元直接相邻的几个栅格高程的平均值。

五、生成子流域上下游拓扑关系表

该步骤对应步骤S5，如图8和图9所示，步骤S5中生成的子流域上下游拓扑关系表中包括当前子流域编码，下游子流域编码，上游子流域编码(复数个)，河道坡面单元标记，平原山区单元标记；其中，子流域编号等于0表示没有上游或下游子流域；河道坡面单元标记等于0表示坡面，等于1表示河道；平原山区单元标记等于1表示坡面，等于0表示山区。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：

本发明提供了一种考虑山区平原地貌差异的子流域单元划分方法，该方法以DEM为基础，考虑山区平原区地貌差异，进行子流域划分和编码，能够保证划分得到的子流域基本符合存在平原区的流域内各单元间的上下游关系，且能根据需求将平原区内的小栅格聚合为一个大的平原单元进行子流域编码，保证子流域范围划分的正确性。该方法能够更好、更全面地刻画流域内不同地貌的地表、地下水流动机制，提高分布式水文模型的水文模拟精度，为分布式水文模型在大尺度流域的构建和应用提供理论和实践支撑。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。