基于栅格型新安江模型的中小河流精细化洪水预报

姚 成，李致家，张 珂，朱跃龙，刘志雨，黄迎春，龚珺夫，张锦堂，童冰星

(1.河海大学水文水资源学院，江苏 南京 210098；2.河海大学计算机与信息学院，江苏 南京 210098；3.水利部信息中心，北京 100053；4.安徽省水文局，安徽 合肥 230022)

我国中小河流众多，洪水频发，灾害严重，是当前防洪工作的重点薄弱环节。据统计，一般年份中小河流的水灾损失占全国水灾总损失的70%～80%，2010年以来水灾造成的人员死亡中有2/3以上发生在中小河流[1]。目前，我国中小河流的洪水预报预警整体水平不高，尚不能满足防洪需求。水文实测、遥感遥测、DEM、土壤植被等组成的水文大数据日益丰富，研究水文大数据驱动的精细化洪水预报已成为提高中小河流预报精度的必然途径[2]。

精细预报要求预报模型或方法能够精准描述中小河流产汇流过程、考虑人类活动影响，同时还能够提供高分辨率、多要素预报产品，以满足新时期影响预报与风险预警的业务需求[3]。大量研究表明，分布式水文模型是实现精细预报的重要支撑[4-11]，如，2015年美国国家海洋和大气管理局联合多个科研机构与政府部门，以分布式模型WRF-Hydro为基础，共同研发了美国国家水模型，美国天气预报中心在2016年的业务预报中，对该模型进行了初步运行，实现了覆盖美国大陆1 km网格上的连续计算，预报成果不仅包括美国270万条河流断面的流量与流速变化，也包括土壤含水量等变量的空间分布，进一步验证了分布式模型在精细化洪水预报中的有效性[3]；2017年，李致家等[12]以正交网格分布式模型为基础，提出了基于蓄满产流与超渗产流空间组合框架的精细化降雨径流模型，考虑参数的地理空间分布，建立了模型参数与下垫面特征间的定量关系，在湿润、半湿润半干旱流域均取得了良好的实际应用效果。

本文以栅格型新安江模型(Grid-Xin’anjiang model，GXM)为基础[13-17]，采用考虑中小河流下垫面条件复杂、人类活动影响大等特点的精细化产汇流模拟方法，并对参数及其空间分布进行客观估计，旨在进一步提高中小河流洪水预报精度，实现流域内任意网格单元流量、土壤含水量等不同水文要素的精细预报，为中小河流影响预报与风险预警业务提供更丰富的预报产品。

1 模 型 原 理

GXM是以正交网格为计算单元，以新安江模型为基础构建的分布式水文模型，该模型由植被冠层截留模块、河道降水模块、三层蒸散发模块、蓄满产流模块、分水源模块、坡地汇流模块与河道汇流模型组成。在产流与分水源计算时，新安江模型采用只有统计意义的概念性元件，即蓄水容量分布曲线，来考虑土壤含水量空间分布不均的问题。而GXM则是利用地形地貌以及土壤植被等下垫面分布特征，确定任意网格单元的张力水蓄水容量与自由水蓄水容量，进而对土壤含水量的空间分布进行具体描述。同时，为了更好地考虑中小河流下垫面条件复杂、人类活动影响大等特点，GXM采用了更具针对性的产汇流模拟方法。

1.1 沿程水流再分配过程模拟

GXM在进行逐网格汇流演算时，对沿程水流的再分配过程进行了模拟。即对于任意网格单元而言，当其土壤含水量未达到田间持水量之前，上游来水优先补充该网格的土壤缺水量，直至其蓄满为止。此外，若网格单元有河道存在，属于河道网格时，则该网格地表径流先按一定的比例汇入河道，然后再汇至下游网格。其中，网格单元上游来水的计算方法为

(1)

式中：Qsup,i——当前网格i的上游网格入流量；Qsout,j——上游第j个网格地表径流的出流量；fch,j——上游第j个网格地表径流出流量汇入河道的比例，对于坡地网格，该比例取0；m——与网格i相邻的上游网格总数。

1.2 库塘坝蓄泄影响的定量化模拟

中小河流内中小型水库、水塘和水坝等蓄水工程种类多、数量大，且大都无闸门控制、缺乏实测水文资料，如何在作业预报中考虑库塘坝蓄泄影响是中小河流洪水预报面临的重要问题。GXM首先根据库塘坝经纬度坐标，确定其所处的网格位置，建立库塘坝的拓扑结构；对于有实测资料的库塘坝，在模型汇流演算中，以库塘坝所处的网格为阻断节点，根据其实测资料进行蓄泄，然后汇流至下游网格；对于无资料库塘坝，引入蓄满率θ及蓄满率临界状态λC，通过设定蓄泄曲线的方式进行蓄泄演算[18]。图1为场次洪水过程中库塘坝时段蓄泄率随蓄满率的变化过程示意图。当θ≤λC时，库塘坝以拦蓄为主；反之，以泄流为主。对于一场洪水过程而言，库塘坝的初始蓄满率θ0低于λC，此时库塘坝仅拦蓄；随着降雨增加，θ增加至λC；若继续降雨，θ>λC，此时库塘坝表现为泄流。

图1 无资料库塘坝蓄泄关系示意图Fig.1 Schematic of water storage and discharge curves of ungauged reservoirs and ponds

对于任意蓄水工程而言，模型假设其蓄泄曲线均可用幂函数表示：

(2)

式中：θt——t时段库塘坝的蓄满率；ωs,t——t时段库塘坝的蓄水率；ms——拦蓄曲线形状参数，反映库塘坝的蓄水能力，一般取值为0～1；ωd,t——t时段库塘坝的泄流率；md——泄流曲线形状参数，反映库塘坝的泄流能力，一般取值为0～3。θ0可通过模型的连续计算，利用土壤含水量和前期累积降雨量进行计算[18]；也可以通过Landsat、高分一号等卫星影像数据，采用最大类间方差迭代提取方法，对库塘坝水体面积进行提取[19]，再根据水体面积-蓄量关系进行计算。

若令泄流为正，则t时段的蓄泄量ΔQt为

(3)

式中：QSUM,t——t时段库塘坝所处网格对应的上游地表径流、壤中流、地下径流及河道径流的累积流量。

1.3 考虑河道特征的数字河道汇流演算方法

GXM在建立之初，采用传统的Muskingum法进行逐网格河道汇流演算。为了进一步适应不同河道特征，并考虑中小河流河道水力要素的时空变化，后又引入基于扩散波(diffusion wave，DW)与Muskingum-Cunge-Todini(MCT)[20]的数字河道汇流演算方法[21]。与Muskingum法相比，DW与MCT方法可以根据河道断面特征，更好地考虑中小河流河道流量、流速、水深等要素的时空变异性，进而更真实地模拟中小河流河道汇流过程。

GXM采用基于地形特征的中小河流河宽估算方法[22]对河道进行数字化，以提高河道网格单元的河宽估算精度，计算方法为

B=δfte+B0

(4)

式中：B——河道网格i的河宽；fte——河道网格单元地形因子；δ——河宽比例系数；B0——基础河宽。其中，δ与B0均可通过遥感影像图量测或断面实测资料进行确定，fte可通过河道网格单元对应的上游累积汇水面积与坡度原点矩进行计算[22]。

由于GXM采用了蓄满产流原理，因此模型理论上仅适用于湿润半湿润地区。若将其用于干旱半干旱地区，可采用蓄满产流与超渗产流的时空动态组合框架[12]。若将模型用于地下水超采地区时，如海河流域中小河流[23-24]，可在模型中增加基于“带门槛水库”的地下拦蓄模块[25]。

2 模型参数空间分布估计

分布式模型网格单元多达数千个，其参数在空间上呈不均匀分布，若只依赖于出口断面实测水文资料，难以保证参数确定和预报结果的合理性。解决该问题的主要途径是以参数物理意义为基础，采用理论推导与数据挖掘等技术，建立参数与地形、土壤、植被等下垫面特征的定量关系，对模型参数的空间分布进行客观估计，以降低参数对实测资料的依赖性，保证参数估计的合理性。

根据现有针对GXM参数及其空间分布估计的研究方法[17,26-27]，植被叶面指数、植被高度与坡面汇流曼宁糙率系数等参数，可以根据网格单元植被类型，利用已有研究结论和文献，直接进行赋值；网格单元张力水蓄水容量、自由水蓄水容量、壤中流与地下径流出流系数等参数，均可以通过其与地形指数、土壤类型、土层厚度等特征的定量关系进行估算；河道网格单元对应的河道宽度、河道长度等参数，可直接利用DEM、遥感影像等进行计算或测量；地表径流汇入河道比例参数可采用基于流向的面积比例法进行估算。在此基础上，再利用实测水文资料，对蓄满率、蓄满率临界状态等剩余参数进行率定或对估算参数进行微调，即可获得每个网格单元的模型参数组。

3 降雨场与下垫面空间分布提取

在GXM应用中，所需的降雨场主要包括实况降雨场和预报降雨场。其中，实况降雨场主要通过两种途径获取：一种是直接利用地面雨量站实时数据，采用泰森多边形、反距离平方等方法进行空间插值；另一种则是利用雷达数据、卫星数据与地面雨量站数据，采用混合地理加权回归和基于贝叶斯与机器学习的多源信息融合等方法[2]，以进一步提高实况降雨场的时空分辨率。对于预报降雨场而言，模型现阶段主要采用WRF模式以及EC、NCEP、GRAPES等数值预报产品[28]。

在进行下垫面空间分布提取时，主要采用DEM高程、植被覆盖/利用、土壤类型与水利普查等数据，数据源包括SRTM数据库、UMD-LULC数据库、HWSD数据库等。同时，利用DEM高程数据，可进一步提取得到流向、水系、坡度、坡向、流径长度和地形指数等地形地貌特征。在此基础上，也可以利用每个网格单元的流向和上游累积汇水面积，获得逐网格汇流演算方法所需的网格单元演算次序[15]。

4 实 例 应 用

GXM自建立以来，已成功应用于江苏、浙江、安徽、陕西等不同省份多个湿润半湿润中小河流的洪水模拟及业务预报，均取得了良好的应用效果[2,13-17]。本文以新安江源头横江为例，开展GXM的应用分析。

4.1 研究区域概况

横江为典型的湿润区中小河流，位于皖南山区，其主要控制站屯溪站以上流域面积为2 670 km2，流域内植被丰富，森林覆盖率达70%以上。屯溪流域降水主要集中在汛期，汛期内降水约占全年总降水量的60%以上，暴雨类型以锋面型暴雨和低压型暴雨为主。屯溪流域属山区性河流，地势西高东低，高差变化悬殊，河道坡度大，洪水陡涨陡落，流速大、历时短。流域内中型水库1座，库容为2 900万m3。流域上游已建小型水库总库容2 200万m3，塘坝总库容2 555万m3，各小型水库与塘坝基本无调度运用方案。屯溪流域水系及中小水库塘坝分布见图2。

图2 屯溪流域水系及中小水库与塘坝分布Fig.2 Water system and spatial distribution of small and medium-sized reservoirs and ponds in Tunxi Basin

4.2 结果分析

GXM在屯溪流域应用时，网格大小采用1 km×1 km，即可以更好地与土壤、植被等数据的网格大小统一，也可以在保证应用精度的前提下，提高模型运算效率[29]。本文利用屯溪站1982—2020年日资料及同期66场洪水过程资料，对模型参数及其空间分布进行了估计、率定与检验。其中，1982—2008年资料用于参数的估计与率定，剩余资料用于参数的检验。结果表明，GXM在屯溪流域应用的径流深合格率为90.9%，洪峰合格率为83.3%，峰现时差合格率为90.9%，确定性系数均值为0.92。图3为屯溪站及休宁站2019年、2020年洪水过程预报结果，休宁站为屯溪站的嵌套站。可以看出，GXM不仅可以实现对屯溪断面出流过程的高精度模拟预报，同时可以在不修正参数及其空间分布情况下，直接给出上游嵌套断面休宁的高精度结果，也进一步验证了参数空间分布估算的合理性与准确性。

图3 屯溪站及休宁站2019年、2020年洪水过程预报结果Fig.3 Forecast results of flood hydrographs in 2019 and 2020 for Tunxi and Xiuning stations

除了预报断面的洪水过程，GXM还可以提供流量、水位、流速等不同水文要素的时空分布。图4为流量、缺水量、库塘坝蓄量变化量以及流速空间分布的模拟结果，验证了GXM对于不同水文要素时空分布预报的能力，相关预报产品可为进一步开展中小河流影响预报与风险预警业务提供支撑。

图4 流量、缺水量、库塘坝蓄量变化量以及流速的空间分布模拟结果Fig.4 Simulation results of spatial distribution of discharge,water shortage,storage change of reservoir and pond,and flow velocity

5 结 语

本文以GXM为基础，针对中小河流下垫面条件复杂、中小水库与塘坝数量大等特点，采用沿程水流再分配方法处理网格单元之间以及网格与河道之间的水流交换，利用库塘坝蓄泄影响的定量化模拟方法反映库塘坝蓄泄对产汇流过程的影响，采用数字河道汇流演算方法模拟中小河流河道特征，并为不同条件下进行逐网格汇流演算提供了Muskingum、DW、MCT等多种演算方法。此外，针对分布式水文模型参数空间分布估算难题，GXM以参数与地形地貌、土壤植被等下垫面特征的定量关系挖掘为基础，对模型参数的空间分布实现客观估计，降低参数对实测水文资料的依赖性。以安徽横江屯溪流域为例，对GXM精细预报效果进行验证，结果表明，GXM不仅可以实现对屯溪出口断面洪水过程的高精度预报，也可以实现对流域内部嵌套断面洪水过程的高精度预报；同时还可以获取流量、土壤含水量、流速、水深等不同要素的时空分布，可以为中小河流洪水防控提供丰富的预报产品。