CN108399309A - 一种大尺度复杂地形区分布式水文模型的子流域划分方法 - Google Patents

Abstract

本发明公开一种大尺度复杂地形区分布式水文模型的子流域划分方法，包括：最大和最小集水面积阈值的确定、流域出水口的确定、多阈值虚拟河网融合与子流域划分。本发明可实现流域出水口的自动识别，高效快速。另外，初步解决了内流区河网提取失真问题，可准确描述所划分子流域间的拓扑关系。在准确拟合研究区范围的同时有效控制了子流域的数量，兼顾了模型模拟的精度和效率需求。该方法在大尺度、复杂地形区的子流域划分上适用性良好，相关结果可对宏观尺度区域分布式水文模型构建给予有力支持。

Description

一种大尺度复杂地形区分布式水文模型的子流域划分方法

技术领域

本发明涉及水利工程技术领域，尤其涉及子流域划分方法，具体为一种大尺度复杂地形区分布式水文模型的子流域划分方法。

背景技术

按照研究对象的地理特征，分布式水文模拟可分为两类：一是针对完整流域的流域分布式水文模拟，二是针对行政区或人工划定区的区域分布式水文模拟。当研究对象空间尺度较大时，为兼顾模拟的精度和效率，一般以GIS为平台，基于栅格DEM提取虚拟河网，将研究区划分为一定数量的子流域，以子流域作为基本模拟单元。如何科学、合理地进行的子流域划分对构建大尺度分布式水文模型至关重要。然而，现有的子流域划分方法仍存在一些不足。一方面，传统的子流域划分方法多以流域为服务对象，服务于区域，特别是大尺度区域的针对性研究较少。另一方面，目前的子流域划分方法并不能完全满足在大尺度复杂地形区域构建分布式水文模型的需求。大尺度复杂地形区的子流域划分往往面临诸如以下难题：（1）研究对象不是完整的流域，不具备“封闭”特性，外边界出水口数目不固定；（2）河流众多且相互独立，水系特征相差较大，内流水系和外流水系并存；（3）复杂的地形地貌导致水系空间分布的差异性较大；（4）高程变化剧烈的微地形区的存在造成研究区内无汇流联系的各流域之间集水面积相差悬殊。若使用传统的基于地表径流漫流模型的子流域划分方法，难以精准确定流域出水口，导致模拟范围与研究区域偏离过大，内流区的真实状况也难以准确反映。另外，单一的集水面积阈值无法兼顾大尺度复杂地形区分布式水文模型计算效率和模拟精度的双重需求。集水面积过大时，研究区模拟面积与实际面积无法良好吻合，随着集水面积的减小，划分的子流域个数会呈指数增加，大大加重模型计算负担。

发明内容

本发明所要解决的技术问题，在于克服现有子流域划分技术存在的缺陷，立足于DEM河网干支拓扑编码规则提出了一种针对大尺度复杂地形区，基于流域出水口自动识别和多阈值虚拟河网融合技术的子流域划分方法。该发明能够在准确把握研究区范围的基础上控制子流域数量，同时兼顾模型的计算精度和效率，为大尺度分布式水文模拟提供支撑。

本发明的目的是这样实现的：一种大尺度复杂地形区分布式水文模型的子流域划分方法，所述方法包括三个部分：包括最大和最小集水面积阈值的确定、流域出水口的确定、多阈值虚拟河网融合与子流域划分。步骤如下：

确定最大和最小集水面积阈值的具体步骤：

步骤一：基于研究区范围裁剪DEM数据，得到剪裁后的DEM数据A。

步骤二：利用实测河网对DEM数据A进行烧录处理，得到河道烧录后的DEM数据B；进行河道烧录的目的在于有效提取出同实测河网高拟合度的模拟河网。

步骤三：逐一排查内流区的实测河网，选择每条内流水系的末端作为该水系所在流域的潜在出口，将该潜在出口所在DEM栅格的高程值设置为NoData，得到设置后的DEM数据C。

步骤四：新DEM 数据C进行洼地填充处理，得到无洼地的DEM数据D。在虚拟河网提取时，由于洼地的存在，导致水流无法顺利流向出水口。所以需要对洼地进行处理，保证流域内所有栅格都能流向流域出水口。

步骤五：采用D8算法，对无洼地的DEM 数据D求取水流方向矩阵，得到水流方向数据E。D8 算法，又称为最陡坡降算法，假设单个栅格中的水流只能流入与之相邻的8 个栅格中的一个。它用最陡坡度来确定水流的方向，即在3×3 的DEM 栅格上，计算中心栅格与各相邻栅格间的距离权落差（即栅格中心点落差除以栅格中心点之间的距离），取距离权落差最大的栅格为中心栅格的流出栅格。

步骤六：对水流方向数据E求取汇流累积量矩阵，得到汇流累积量数据F。汇流累积量是指流入当前栅格的所有上游栅格的个数，作为当前栅格的汇流累积量值，其值大小代表了上游汇水面积的大小。

步骤七：初设一个集水面积阈值G，对汇流累积量数据F进行重分类处理，生成单值的虚拟河网栅格数据H，使得大于该阈值的栅格值为1，小于该阈值的栅格值设定成NoData。将虚拟河网栅格转为矢量数据，计算虚拟河网密度I。

步骤八：设置不同的集水面积阈值，重复步骤七，统计虚拟河网密度I随集水面积阈值G的变化关系曲线，找出模拟河网密度I随阈值G的增大变化不再显著的点，选择该点对应的集水面积阈值作为本发明的最大集水面积阈值J。

步骤九：以1为步长，调小最大集水面积阈值J，重复步骤七，生成虚拟河网栅格数据H。

步骤十：根据虚拟河网栅格数据H，划分研究区的子流域K，计算研究区模拟范围的面积。

步骤十一：当研究区模拟范围吻合程度（研究区模拟面积与实际面积之比）小于99.9%时，重复步骤九和步骤十，直至研究区模拟范围吻合程度（研究区模拟面积与实际面积之比）达到或超过99.9%。选择此时对应的集水面积阈值作为本发明的最小集水面积阈值L。

确定流域出水口的具体步骤：

步骤一：基于水流方向数据E，按照D8流向编码，遍历研究区外边界以及内流区流域潜在出口栅格的周围8个栅格的流向，如果某栅格流向的下一个栅格无流向值，则将此栅格标记为“备选”流域出水口，生成“备选”流域出水口数据M。

步骤二：利用汇流累积量数据F，确定“备选”流域出水口数据所在栅格的汇流累积量，若得到的汇流累积量大于最小集水面积阈值L，则选择该出水口作为最终的流域出水口。遍历所有的“备选”流域出水口，生成流域出水口数据N。对每一个流域出水口设定唯一值。

多阈值虚拟河网融合与子流域划分的具体步骤：

步骤一：以流域出水口数据N所在栅格为起点O，根据水流流向栅格判断识别该栅格的上游河道汇入栅格，可能有一个或多个。若为一个，则将此上游河道栅格直接作为干流汇入栅格P，赋值为1。若为多个，比较其汇流累积数大小，选择其中最大值所在栅格作为干流汇入栅格P。

步骤二：将步骤一得到的干流汇入栅格P作为起点O，重复步骤一，继续向上游追溯，直至河道源头。最终提取所有的干流汇入栅格P，生成干流虚拟河网Q。

步骤三：利用最大集水面积阈值J和汇流累积量数据F，提取生成虚拟河网R。

步骤四：将干流虚拟河网Q和虚拟河网R融合，生成新的虚拟河网S。

步骤五：基于水流方向数据E和汇水交汇点（节点），对虚拟河网S进行分割，并对每一段河网设定唯一值，得到分割后的汇水网格数据T。

步骤六：基于分割后的汇水网络数据T和水流方向数据E，划分得到研究区的子流域U。

本发明产生的有益效果是：1）可实现流域出水口的自动识别，高效快速；2）初步解决了内流区河网提取失真问题，可准确描述所划分子流域间的拓扑关系；3）在准确拟合研究区范围的同时有效控制了子流域的数量，兼顾了模型模拟的精度和效率需求。该方法在大尺度、复杂地形区的子流域划分上适用性良好，相关结果可对宏观尺度区域分布式水文模型构建给予有力支持。

附图说明

下面结合附图和具体实施例，对本发明作进一步详细说明。

图1是本发明的实施例所述方法的流程图；

图2是本发明确定的我国内流区流域潜在出口分布图；

图3是本发明确定的全国虚拟河网密度随集水面积阈值的变化关系曲线图；

图4是本发明确定的全国模拟范围吻合程度和子流域个数随集水面积阈值的变化关系曲线图；

图5是本发明所述240km²集水面积阈值下我国沿海地区子流域划分结果局部示意图；

图6是本发明采用的D8 算法编码示意图；

图7是本发明确定的全国流域出水口分布图；

图8是本发明确定的海南岛虚拟河网融合示意图；

图9是本发明得到的全国子流域划分结果图；

图10是本发明得到的全国子流域划分结果局部示意图。

具体实施方式

下面结合附图和具体实施例，对本发明做进一步详细说明：

我国是一个典型的大尺度复杂地形区，具有以下几个显著特征：1）幅员面积大，南北跨越约50个纬度，东西横跨61个经度；2）河湖众多，内流水系和外流水系并存；3）海岸线长，近海地区地形起伏变化剧烈。以我国为例展开，具有较好代表性。

所用DEM数据为美国太空总署（NASA）和国防部国家测绘局（NIMA）联合测量、生产的90m×90m的SRTM栅格数据，将其重采样至1km×1km的DEM栅格数据作为本实例展开的基础数据。同时，全国1:400万实测水系和全国范围面矢量地图数据用于修正DEM数据。实例以ArcGIS10.2软件为平台，编程计算依靠Fortran语言实现。

第一步：确定全国子流域划分的最大和最小集水面积阈值

（1）DEM数据裁剪。利用Extract by Mask 工具，基于全国范围面矢量地图数据裁剪得到全国DEM数据A。

（2）河道烧录处理。利用Polyline to Raster工具将实测水系线矢量数据转为栅格数据，然后利用Raster Calculator工具将全国DEM数据A中实测水系所在栅格高程降低1000m，其余栅格高程保持不变，得到河道烧录处理后的全国DEM数据B。

（3）排查我国范围内的内流水系，在所有四级及以上河流和部分五级河流的末端人工设定流域潜在出口，共计73个（图2）。将全国DEM数据B中流域潜在出口所在栅格高程设置为空，得到新的全国DEM数据C。

（4）填洼处理。利用Fill 工具对全国DEM数据C中的洼地进行填充处理，得到无洼地的全国DEM数据D。

（5）水流方向计算。基于全国DEM数据D，利用Flow Direction 工具计算水流方向，得到全国水流方向数据E。

（6）汇流累积量计算。基于全国水流方向数据E，利用Flow Accumulation 工具，计算得到全国汇流累积量数据F。

（7）全国最大集水面积阈值确定。在10~1000km²的范围内选取若干集水面积阈值，生成单值的全国虚拟河网栅格数据H，计算所对应的全国虚拟河网密度I，绘制两者变化关系曲线（图3）。可以发现，阈值在200~300km²范围内曲线趋于稳定。通过进一步排查，最终确定以240km²作为全国范围最大集水面积阈值J。此时得到的全国模拟范围吻合程度（全国模拟面积与实际面积之比）为98.51%，沿海岸线大量高程起伏剧烈的微地形区未能完成子流域划分（图4和图5）。

（8）全国最小集水面积阈值确定。不断调小全国集水面积阈值，当全国集水面积阈值降低至4km²时，全国模拟范围吻合程度（全国模拟面积与实际面积之比）达到99.92%。选择以4km²作为全国范围最小集水面积阈值L，此时划分的子流域个数增至100.26万，严重增加模型计算负荷（图4）。

第二步：确定全国流域出水口

（1）基于全国水流方向数据E，按照D8流向编码（图6），遍历全国边界以及内流区流域潜在出口栅格的周围8个栅格的流向，如果某栅格流向的下一个栅格无流向值，则将此栅格标记为“备选”流域出水口，生成全国“备选”流域出水口数据M。

（2）利用全国汇流累积量数据F，确定全国“备选”流域出水口数据所在栅格的汇流累积量，若得到的汇流累积量大于4km²，则选择该出水口作为最终的流域出水口。遍历全国所有的“备选”流域出水口，提取生成全国流域出水口数据N。全国范围内共确定4121个流域出水口，分布于海岸线、边境线和内流区（图7）。

第三步：多阈值虚拟河网融合与子流域划分

（1）以全国4121个流域出水口为起点逆流而上，通过全国水流流向矩阵和水流累积矩阵向上游追迹流域出水口所在河网干流，赋值为1，生成全国干流虚拟河网Q。

（2）以240km²集水面积阈值，利用Con工具提取全国虚拟河网R。

（3）利用Raster Calculator工具将全国干流虚拟河网Q和全国虚拟河网R融合，生成新的全国虚拟河网S。以海南岛进行局部示意，见图8。

（4）利用Stream Link 工具，基于全国虚拟河网S和水流方向数据E，对全国虚拟河网S进行分割，得到分割后的全国汇水网格数据T。

（5）利用Watershed 工具，基于分割后的全国汇水网络数据T和水流方向数据E，划分得到全国子流域U。全国共划分25174个子流域，子流域的平均面积为381km²，最大面积为5512km²，最小面积为10km²，全国模拟面积与实际面积之比达到99.92%，效果良好（图9和图10）。通过本发明，既明显提高了模型模拟精度，又极大减小了子流域个数，避免了计算灾难。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案（比如各种公式的运用、步骤的先后顺序等）进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。

Claims (1)

1.一种大尺度复杂地形区分布式水文模型的子流域划分方法，其特征在于，所述子流域划分方法包括三个部分：确定最大和最小集水面积阈值、确定流域出水口、多阈值虚拟河网融合与子流域划分；其具体步骤是；

确定最大和最小集水面积阈值：

步骤一：基于研究区范围裁剪DEM数据，得到裁剪后的DEM数据A；

步骤二：利用实测河网对DEM数据A进行烧录处理，得到河道烧录后的DEM数据B；

步骤三：逐一排查内流区的实测河网，选择每条内流水系的末端作为该水系所在流域的潜在出口，将该潜在出口所在DEM栅格的高程值设置为NoData，得到设置后的DEM数据C；

步骤四：对DEM 数据C进行洼地填充处理，得到无洼地的DEM数据D；

步骤五：采用D8算法，对无洼地的DEM 数据D求取水流方向矩阵，得到水流方向数据E；

步骤六：对水流方向数据E求取汇流累积量矩阵，得到汇流累积量数据F；

步骤七：初设一个集水面积阈值G，对汇流累积量数据F进行重分类处理，生成单值的虚拟河网栅格数据H，使得大于该阈值的栅格值为1，小于该阈值的栅格值设定成NoData；将虚拟河网栅格转为矢量数据，计算虚拟河网密度I；

步骤八：设置不同的集水面积阈值，重复步骤七，统计虚拟河网密度I随集水面积阈值G的变化关系曲线，找出虚拟河网密度I随集水面积阈值G的增大变化不再显著的点，选择该点对应的集水面积阈值作为最大集水面积阈值J；

步骤九：调小最大集水面积阈值J，重复步骤七，生成虚拟河网栅格数据H；

步骤十：根据虚拟河网栅格数据H，划分研究区的子流域K，计算研究区模拟范围的面积；

步骤十一：当研究区模拟范围吻合程度小于99.9%时，重复步骤九和步骤十，直至研究区模拟范围吻合程度达到或超过99.9%，选择此时对应的集水面积阈值作为最小集水面积阈值L，所述研究区模拟范围吻合程度为研究区模拟面积与实际面积之比；

确定流域出水口的具体步骤：

步骤一：基于水流方向数据E，按照D8流向编码，遍历研究区外边界以及内流区流域潜在出口栅格的周围8个栅格的流向，如果某栅格流向的下一个栅格无流向值，则将此栅格标记为“备选”流域出水口，生成“备选”流域出水口数据M；

步骤二：利用汇流累积量数据F，确定“备选”流域出水口数据所在栅格的汇流累积量，若得到的汇流累积量大于最小集水面积阈值L，则选择该出水口作为最终的流域出水口；遍历所有的“备选”流域出水口，生成流域出水口数据N；对每一个流域出水口设定唯一值；

多阈值虚拟河网融合与子流域划分的具体步骤：

步骤一：以流域出水口数据N所在栅格为起点O，根据水流流向栅格判断识别该栅格的上游河道汇入栅格，可能有一个或多个；若为一个，则将此上游河道栅格直接作为干流汇入栅格P，赋值为1；若为多个，比较其汇流累积数大小，选择其中最大值所在栅格作为干流汇入栅格P；

步骤二：将步骤一得到的干流汇入栅格P作为起点O，重复步骤一，继续向上游追溯，直至河道源头；最终提取所有的干流汇入栅格P，生成干流虚拟河网Q；

步骤三：利用最大集水面积阈值J汇流累积量数据F，提取生成虚拟河网R；

步骤四：将干流虚拟河网Q和虚拟河网R融合，生成新的虚拟河网S；

步骤五：基于水流方向数据E和汇水交汇点，对虚拟河网S进行分割，并对每一段河网设定唯一值，得到分割后的汇水网格数据T；

步骤六：基于分割后的汇水网络数据T和水流方向数据E，划分得到研究区的子流域U。