技术领域
[0001] 本
发明属于
水文
地貌学领域，特别涉及一种基于湖泊汇流关系的内流湖盆区流域自动划分方法。
背景技术
[0002] 流域是天然的地理单元，流域的空间信息广泛应用于水文建模(任立良和刘新良，2000；Liu et al.,2014)，地貌演化(祝士杰等，2013)、
土壤侵蚀(Vente et al.,2013)等众多领域。当前，
对流域信息的提取主要基于数字高程模型，通过开发相应的流域地形分析
算法，实现对流域边界的自动划分。
[0003] 流域划分的准确性主要依赖于数字高程模型的
精度和流域提取算法的性能。随着航空/天数字摄影测量、激光点
云、雷达测绘等技术的快速发展，地形数据的获取能
力有了极大的提升。当前在全球尺度上广泛应用的数据高程模型数据主要有SRTM DEM、ASTER GDEM、ALOSAW3D30等，这些数据也成为流域划分中的重要
基础数据(Yamazaki et al.,2017；Yue et al.,2017)。在算法方面，目前流域划分形成了较为固定的处理流程，主要包括了洼地填充、水流方向计算、汇流累计计算、河网提取及流域标识等步骤(O’Callaghan and Mark,1984；Tarboton et al.,1991；Freitas et al.,2016)。该算法流程主要针对外流区汇流特征设计，其基本假设在于，研究区内全部栅格产生的汇流均能够通过出水口栅格或边界栅格流出，这种思路能够确保所提取的河网信息具有较好的结构性特征，同时通过设定不同汇流面积
阈值，可以提取出具有层次嵌套结构的流域边界。
[0004] 虽然现有流域划分方法在外流区取得了较好结果，但并不适合于内流湖盆区的地表汇流特征。一个典型的内流湖盆区是由若干独立的内流盆地组合而成，不同内流盆地间的汇流网络是互不相连的，因此传统流域提取算法无法刻画出内流盆地的流域地形特征(Khan et al.,2014)。当前，内流湖盆区以其独特的自然环境得到学者们的广泛关注，典型代表是青藏高原内流区。学者们主要使用HydroSHEDS(http://hydrosheds.cr.usgs.gov)全球尺度流域数据集作为流域边界的参考数据。但该数据一方面在流域边界的表达上精度较差，另一方面也被指出存在着较多错误划分的流域边界，这些问题严重影响了内流湖盆区相关地学研究。部分学者也注意到内流湖盆区流域划分产品所存在的问题，并探讨了融合地形图、影像等多源数据的内流湖盆区流域划分结果的修正方法(Gao et al.,2018)。但这种思路对数据要求较高，很难推广到数据资料较为缺乏区域，同时因为需要人工判断的大量介入，处理的效率低，不适合于大区域尺度的研究。在此背景下，开展内流湖盆区流域自动划分方法的研究，可为内流湖盆区相关科学问题的探索提供重要的数据保障和方法
支撑，具有重要的科学意义。
[0005] 参考文献：
[0006] [1]任立良,刘新仁.数字高程模型信息提取与数字水文模型研究进展[J].水科学进展,2000,11(4):463-469.
[0007] [2]Liu J,Zhu A X,Liu Y,et al.,2014.A layered approach to parallel computing for spatially distributedhydrological modeling.Environmental Modelling&Software.51,221-227.
[0008] [3]祝士杰,汤国安,李发源,等.基于DEM的黄土高原面积高程积分研究[J].地理学报,2013,68(7):921-932.
[0009] [4]Vente J D,Poesen J,Verstraeten G,et al.,2013.Predicting soil erosion and sediment yield at regional scales:Where do we stand？.Earth Science Reviews,127(2),16-29.
[0010] [5]Yamazaki,D.,Ikeshima,D.,Tawatari,R.,et al.,2017.A high accuracy map of global terrain elevations.Geophysical Research Letters,322(22),109-127.
[0011] [6]O’Callaghan,J.F.,Mark,D.M.,1984.The extraction of drainage networks from digital elevation data.Comput.Vision,Graph.,Image Process.28,323–344.
[0012] [7]Tarboton,D.G.,Bras,R.L.,Rodriguez-Iturbe,I.,1991.On the extraction of channel networks from digital elevation data.Hydrological Processes,5(1),81-100.
[0013] [8]Freitas H R D A,Freitas C D C,Rosim S,et al.,2016.Drainage networks and watersheds delineation derived from TIN-based digital elevation models.Computers&Geosciences.92(C):21-37.
[0014] [9]Gao,Y.,Wang,W.,Yao,T.,et al.,2018.Hydrological network and classification of lakes on the ThirdPole.Journal ofHydrology,560,582-594.[0015] [10]Khan A,Richards K S,Parker G T,et al.,2014.How large is the Upper Indus Basin？The pitfalls ofauto-delineation using DEMs.Journal ofHydrology.509(4),442-453.
发明内容
[0016] 为解决现有内流湖盆区流域划分数据产品存在较多错误、缺乏自动化方法、无法准确提取出流域边界及水文地貌参数的问题，本发明提出了一种基于湖泊汇流关系的内流湖盆区流域划分方法，采用可公开获取的数字高程模型数据和湖泊分布数据，实现对内流湖盆区子流域的自动划分。
[0017] 为实现上述技术目的，本发明采用如下技术方案：
[0018] 基于湖泊汇流关系的内流湖盆区流域自动划分方法，包括以下过程：
[0019] 步骤1、基于数字高程模型与湖泊分布数据，计算水流方向，获取水流方向矩阵；
[0020] 步骤2、基于水流方向矩阵提取水文结构参数，包括径流源点、汇流网络、汇流网络与湖泊交点，并基于水流方向矩阵和汇流网络生成初始流域分割数据；
[0021] 步骤3、对于通过汇流网络相连的每一组湖泊，基于湖面高程和湖泊间的高程剖面，判断湖泊间汇流关系；
[0022] 步骤4、基于湖泊间汇流关系将湖泊划分为封闭型湖泊和过水型湖泊两类；并确定每一个过水型湖泊所对应的封闭型湖泊；
[0023] 步骤5，对于封闭型湖泊，基于径
流线追踪法结合峰值点控制，对步骤2中生成的初始流域分割数据进行归属判定，得到与封闭型湖泊相对应的流域划分结果。
[0024] 本发明所述的数字高程模型可选用SRTM-1/3GDEM、ALOS AW3D30或TanDEM等公开地形数据；湖泊分布数据可采用Pekel等公布的全球
水体数据集。
[0025] 作为本发明的进一步改进，对于洼地区域，采用Priority-Flood算法进行洼地流向的追踪；对于平地区域，采用G&M平地流向处理算法；其余区域采用D8算法进行流向计算。
[0026] 作为本发明的进一步改进，所述步骤2中，在水流方向矩阵的基础上，计算每个栅格的汇流累积值，同时根据预设的汇流累积阈值确定河网栅格，形成汇流网络矩阵；结合汇流网络矩阵和水流方向矩阵确定径流源点栅格；结合汇流网络矩阵和
栅格化后的湖泊数据，确定汇流网络与湖泊的交点；基于水流方向矩阵和汇流网络矩阵，生成河流连接矩阵，以实现对每段河流的编码，在此基础上，生成初始的流域分割数据。
[0027] 进一步的，径流源点、汇流网络和初始流域分割数据均采用相同的径流阈值确定；径流阈值优选为3-5平方公里。
[0028] 作为本发明的进一步改进，对于通过汇流网络相连的一组湖泊，生成湖泊间汇流网络所对应的高程剖面，根据剖面的整体趋势，判断湖泊间是否有地形阻隔；若有地形阻隔，则两湖泊不相连，若没有地形阻隔，则根据湖面高程判断两湖泊间的汇流方向，本发明中根据比较相连湖泊的湖面高程中位数判断湖泊间的汇流方向。
[0029] 其中，基于高程剖面的汇流关系确认方式具体为：
[0030] 计算高程剖面线上每个栅格的高程导数，当两个连续栅格的高程导数由正变负，则后一个栅格判定为局部峰值点；反之则为局部谷底点，谷底点和峰值点之间的栅格为高程上升区间；
[0031] 计算每段高程上升区间与整个剖面的长度比值和幅度比值，当长度比值或幅度比值≤阈值时，则认为该高程剖面起点湖泊将汇入高程剖面的终点湖泊；当长度比值或幅度比值＞阈值时，则认为该高程剖面两端湖泊不连通，同时需要记录该高程剖面的局部峰值点；
[0032] 进一步的，所述步骤3中，对于空间距离小于1公里的湖泊，直接采用湖泊的湖面高程判断湖泊间的汇流方向，默认两者间不存在地形阻隔，通过比较湖泊湖面高程，判断出湖泊间的汇流方向；反之则生成高程剖面，根据剖面的整体趋势，判断湖泊间是否有地形阻隔；进一步的，每一个栅格高程值采用以其为中心的5*5窗口的中位数，生成高程剖面线后进一步通过高斯滤波消除异常值干扰。
[0033] 作为本发明的进一步改进，所述步骤5中，径流线追踪的起点为径流源点和径流线与湖泊交点；依次处理全部径流源点和径流交点，沿汇流网络追踪直至遇到局部峰值点或湖泊区域，如果追踪终点是局部峰值点，则将全部经过的子流域暂不赋值；如果追踪的终点是封闭型湖泊，则将全部经过的子流域赋值为终点湖泊的编号；如果终点湖泊是过水型湖泊，则将全部经过的子流域赋值为该过水型湖泊对应的封闭型湖泊的编号；最后从局部峰值点开始沿汇流网络正向(顺着水流方向)和逆向(逆着水流方向)追踪直至遇到湖泊，如果追踪的终点是封闭型湖泊，则将全部经过的子流域赋值为终点湖泊的编号；如果终点湖泊是过水型湖泊，则将全部经过的子流域赋值为该过水型湖泊对应的封闭型湖泊的编号。
[0034] 其中，所述局部峰值点指，湖泊间汇流网络对应的高程剖面线上两个连续栅格的高程导数由正变负，则后一个栅格判定为局部峰值点。
[0035] 当完成全部初始子流域归属湖泊的划分后，采用ArcGIS
软件中Dissolve工具，可实现子流域的归并，即确定出每一个封闭型湖泊的流域范围。
[0036] 作为本发明的进一步改进，还包括步骤6，对于步骤5划分结果中流域面积超过预期面积的流域，设定集水区面积阈值，对其内部进行进一步划分，集水区面积阈值根据研究或应用需要设定，以满足使用需要。
[0037] 本发明具有以下两点优势：
[0038] (1)本发明提出了适用于内流湖盆区的流域自动划分方法，具有较高的划分精度和处理效率，为内流湖盆区的地貌、水文、生态等研究提供了重要的基础数据和技术支撑。
[0039] (2)本发明对数据的要求较低，基于可公开获取的数字高程模型和湖泊分布数据即可完成流域划分。同时当采用高
分辨率地形数据时，该方法也可充分挖掘高
质量地形数据的优势，提取出与真实地表更加吻合的流域边界。
附图说明
[0042] 图3本发明实施例中采用不同流向算法得到的汇流网络：(a)传统D8算法；(b)本发明中综合多种流向算法结果。
[0043] 图4本发明实施例中所得到水文结构特征提取结果：(a)全区域；(b)局部区域放大。
[0044] 图5本发明实施例中采用的基于湖泊间水位判断汇流关系示意图。
[0045] 图6本发明实施例中采用的基于径流高程剖面判断湖泊汇流关系示意图。
[0046] 图7本发明实施例中得到的湖泊类型划分结果。
[0047] 图8本发明实施例中采用的径流线追踪方法得到封闭型湖泊流域范围示意图。
[0048] 图9本发明实施例中得到的流域划分结果：(a)基于ALOS AW3D30数据；(b)HydroSHEDS数据提取结果。
具体实施方式
[0049] 下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施案例用于说明本发明，但不用来限制本发明的范围。
[0050] 本
申请的实施例以青藏高原内流区的部分区域为研究区。青藏高原被誉为“世界
屋脊”和“亚洲水塔”，该区域地形复杂、
气候独特、湖泊
冰川广泛分布，是当前全球变化研究的核心区域。青藏高原内流区主要分布在青藏高原的西北部，主要指羌塘高原和柴达木盆地等区域，该区域分布有大量低洼封闭湖盆，河流汇入这些洼地并形成为数众多的咸水湖。据统计，该区域大于1平方公里的湖泊约1200个。实施例选择了纳木错和色林错附近区域，总面积共计119,055km2。实验数据选择了ALOSAW3D30，数据分辨率为30米。
[0051] 如图2所示，是实施例的流程图，本实施例包括以下步骤：
[0052] 步骤一、综合多种流向算法的水流方向计算。如图3所示，基于传统D8算法在该区域所生成的汇流网络有明显的平行河网现象，特别是湖泊等平坦区域。而本实施例中通过对洼地、平地和其他起伏区域分别采用Priority-Flood、G&M算法和D8算法，所得到的汇流水网络较好的表达了该区域的汇流结构特征，消除了平行河网问题，为后续分析创造条件。
[0053] 步骤二、在流向计算的基础上，进一步生成了该区域水文结构特征。如图4所示，水文结构特征包含了径流源点、汇流网络、汇流网络与湖泊交点，以及初始子流域划分结果。
[0054] 步骤三、湖泊间汇流关系判断。对于如图5中所示空间距离较近湖泊(一般指＜1公里)，本实施例通过比较两个湖泊的湖面高程的中位数进行判断，湖面高程较高的湖泊汇入湖面高程较低的湖泊。对于空间距离较远的湖泊间汇流关系判断，首先记录下两个湖泊间的高程剖面信息，并运用高斯滤波消除噪声点；其次通过计算各坐标点的高程导数，确定出高程的上升区间；最后将每段高程上升区间与整个剖面的幅度比值及长度比值与阈值比较，判断出该上升区间是否能表征真实的地表阻隔，如果幅度比值或长度比值＞阈值，则可认为则该高程上升区间是由真实地形阻隔引起，如果均≤阈值，则判断为误差。如图6a中，两湖泊间虽然存在一个高程上升区间，但是其与整个剖面的幅度比值和长度比值均极小，可认为是数据误差造成的，因此两湖泊间可判断为连通状态。图6b中显示的是典型的不连通情况，高程剖面存在一个非常明显的上升区间，表明两湖泊间被山地阻隔。也有部分湖泊间是否连通并不明确，如图6c所示，根据高程剖面，两湖泊间存在有局部峰值点的可能性，这种情况下，阈值大小对判断结果有直接影响。在本实施例中，幅度比值阈值和长度比值阈值均设置为0.15。
[0055] 步骤四、湖泊类型划分。根据步骤三的判断结果，将只有汇入没有流出的湖泊定义为封闭型湖泊；将既有汇入又有流出的湖泊定义为过水型湖泊。如图7所示，该区域共有封闭型湖泊46个，过水型湖泊95个。
[0056] 步骤五，在湖泊类型划分的基础上，结合初始采用3km2面积阈值的子流域划分结果，确定出每一个封闭湖泊的范围。如图8所示，从径流源点沿着汇流网络进行追踪，直至遇到湖泊，由于该湖泊为封闭型湖泊，则将追踪经过的全部子流域赋值为该湖泊的编号。追踪起点还可以是汇流网络与湖泊的交点，这种情况下，当路径存在局部峰值点，即径流线两端的湖泊已经判断出是不连通的，则追踪过程中所经过的子流域暂不赋值，后续以局部峰值点为起点顺着水流方向和逆着水流方向分别追踪直至遇到湖泊，将所经过的全部子流域均赋值为终点封闭型湖泊的编号(图8b)。最后通过矢量数据的融合处理，可得到每一个封闭型湖泊对应的流域范围(图8c)。该区域全部处理后的结果如图9a所示，对比图9b中已有的HydroSHEDS数据可发现，本方法提取出的结果在部分子流域有明显不同，如本方法提取出的色林错流域面积约为45454.52km2远大于HydroSHEDS数据中28978.40km2的面积。而HydroSHEDS在该流域的划分结果已被证明具有
明显错误。
[0057] 步骤六、对于图9的第一次层级划分结果，如果部分子流域的面积超过实际应用需要，可在第一层次划分结果的基础上，通过设定集水区面积阈值，进行子流域的细化。
