基于多源降水融合驱动的WRF-Hydro模型在中小河流洪水预报中的适用性

晁丽君，张 珂，陈新宇，王国庆

(1.河海大学水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；2.长江保护与绿色发展研究院，江苏 南京 210098；3.河海大学水文水资源学院，江苏 南京 210098；4.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，江苏 南京 210029；5.南京水利科学研究院水利部应对气候变化中心，江苏 南京 210029)

“降雨之后发生了什么”是水文学最基本的一个科学问题[1]。局部强降水容易引发洪水，尤其在地形复杂的山区型中小河流。我国有数以万计的中小河流，这些河流洪水频发、灾害严重，是当前洪水预报应重点关注的区域[2]。中小河流的洪水具有发生快、破坏性强、灾害严重的特点，因此其洪水预报应具有及时性、准确性和可靠性[2]等特点。随着对地观测技术、遥感技术以及一系列反演算法的发展，获取空间连续水文要素已成为现实[3-5]。降雨作为洪水预报过程中最重要的驱动要素，其精度、时效性和可靠性直接影响了洪水预报的精度与可靠性[6]。有关研究表明，降雨径流过程不确定性的70%～80%来源于降雨质量和时空变异性误差[7]，因此，准确获取高时空分辨率的降水数据对于中小河流的水文过程分析、洪涝灾害的预报预警具有重要的现实意义[8]。目前，星地多源降水融合已成为提高降水产品精度的一种常用方法[9]。

气候变化和人类活动改变了天然的产汇流规律，对中小河流的影响尤为突显，提高中小河流洪水预报的精度[10-11]势在必行。自20世纪30年代以来，作为洪水预报的重要手段[12]，用于洪水预报或预测的降雨径流模型层出不穷[7]。随着计算机和GIS技术的发展，美国、欧洲、中国等国家开发了一系列的洪水预报系统，如水利部水文局开发的中国洪水预报系统(NFFS)[13]、美国气象局开发的美国洪水预报系统(AHPS)、美国水文研究中心(HRC)开发的突发性洪水预警指南系统(CAFFG)、欧盟委员会联合研究中心(JRC)研发的欧洲洪水预警系统(EFAS)、欧洲中期天气预报中心(ECMWF)开发的集合预报系统(EPS)[14]等。近十几年来，洪水预报从大尺度向精细化、网格化、智能化预报预警方向发展[15-17]，美国研发的新一代国家水文模型(NWM)即WRF-Hydro模型框架，提高了区域的洪水预报能力[18-21]。本文通过星地多源降水融合获取高时空分辨率的降水数据，并驱动WRF-Hydro模型，探讨了WRF-Hydro模型在中小河流洪水预报中的适用性和合理性，为WRF-Hydro模型在中国区域的应用提供参考。

1 研究区域与数据简介

1.1 研究区域

以典型的山区中小河流子午河为研究区。子午河系汉江上游北岸的一级支流，地理位置为33°18′N～33°44′N，107°51′E～108°30′E(图1)。子午河属山溪性河流，发源于秦岭南麓，主要流经秦岭深山区，流域呈扇形。河流全长161 km，流域面积3 028 km2，河道平均比降0.544%。子午河由上游的椒溪河、蒲河、汶水河和下游的长安河汇合而成。主源汶水河发源于宁陕、周至、户县交界的秦岭南麓，由东北向西南流径宁陕县境内，在宁陕县与佛坪县交界处与浦河、椒溪河汇合后称之为子午河。子午河由北向南于石泉县三华石乡白沙渡附近入汉江。三河口(椒溪河、蒲河、汶水河交汇处)以上汶水河全长106 km，流域面积为1 094 km2，河道平均比降0.93%；椒溪河长70 km，流域面积596 km2，河床平均宽60 m，河道平均比降1.87%，年平均径流量2.665亿m3；蒲河长58 km，流域面积496 km2，河道平均比降0.266%，平均径流深430 mm，平均流量为5.1 m/s。

图1 子午河流地理位置、水系、子流域、站点分布 Fig.1 Map of the Ziwuhe watershed with location, river network, division of sub-watersheds, and locations of the rain gauge and hydrological stations

子午河流域地势北高南低，大部分由花岗岩、片麻岩等变质岩组成，山高谷深、沟壑纵横。秦岭主峰海拔2 965 m，为流域最高峰。流域内山地因受河流切割，多形成自南向北伸展的山梁，谷底切割较深，峡谷壁立，海拔在1 000 m左右，因此水能资源丰富。流域内森林资源丰富，植被良好，林木茂密，森林覆盖率达70%。

子午河流域属于山地温带气候，其气候温和、雨量充沛。由于秦岭的阻断作用，夏季来自南方和东南方的潮湿空气携带大量的水汽，导致灾害天气频繁。子午河流域的降雨具有明显的季节性，夏季降水量占全年降水量的46%，而冬季降水量仅占全年的3%。多年平均降水量、蒸发量、径流量分别为：900 mm/a、500 mm/a、400 mm/a。降雨季节的分布不均，加上陡峭的地势，容易引发山洪和季节性洪水。

1.2 数据来源

研究数据分为气象驱动数据、实测降水和流量数据、静态下垫面数据。气象驱动数据作为WRF-Hydro模型的主要输入，包括降雨、温度、风速、气压、湿度、长波辐射、短波辐射。降水数据采用由CMORPH(CPC MORPHing technique)卫星降水与站点降水产生的融合降水数据，融合降水的时间和空间分辨率分别为1 h和1 km。 卫星降水数据源为CMORPH降水产品，来自美国气候预测中心(Climate Prediction Center，CPC)全球降水数据，其最小空间分辨率为8 km，时间分辨率最高可达0.5 h，空间覆盖范围为南北纬60°。

其他气象驱动数据来自全球陆面数据同化系统(Global Land Data Assimilation System，GLDAS)的产品(https：//disc.gsfc.nasa.gov/datasets?keywords=GLDAS)，在全球范围内提取流域范围内的数据，并通过空间上双线性插值和时间上线性内插，得到空间分辨率为1 km和时间分辨率为1 h的除降水之外的其他气象驱动数据。

2011年之前子午河内仅有5个雨量站(佛坪、四亩地、新厂街、钢铁、筒车湾)和2个水文站(大河坝、两河口)；2012年国家实施中小河流建设之后，共有16个雨量站、2个水文站。2015年之后在大河坝断面下游2 km处兴建了一座水电站。本文选用2013—2014年子午河流域18个站点的时段降雨和1个水文站(两河口)的时段流量用于WRF-Hydro模型计算结果和融合降水的验证，既具有代表性，又不受人类活动的影响。为了有效评估融合降水对中小河流洪水预报结果的影响，选用两河口水文站2013年和2014年共5场洪水(表1)进行模型计算。

表1 洪水场次和洪水起始时间

静态下垫面数据包括地形数据、土壤类型、植被类型、植被覆盖度。地形数据来自美国地质调查局(USGS)镜像网站(http：//www.hydrosheds.org)基于航天飞机雷达的多尺度水文数据和地图(HydroSHEDS)，包括1 km的DEM(图2(a))、流向、汇流累积量等。土壤和植被是WRF-Hydro模型中重要的静态输入数据，直接关系两个参数文件(VEGPATM.TBL和SOILPARM.TBL)的取值。土壤和植被数据来自WRF天气模式预处理WPS(WRF Preprocessing System)系统的静态数据库(WPS Geographical Static Data)，或WPS Geographical Input Data数据集。子午河流域内植被类型为常绿阔叶林、落叶针叶林、落叶阔叶林、稀疏草原(图2(b))，土壤类型分别为：黏土、砂土、砂壤土、壤砂土(图2(c))。植被覆盖度由MODIS的植被归一化指数(NDVI)计算得到，该流域植被覆盖度的范围为0.06～0.20(图2(d))。

图2 流域高程、土壤和植被空间分布Fig.2 Spatial distributions of watershed elevation, soil and vegetation

2 融合算法与WRF-Hydro模型

2.1 混合地理加权回归截尾函数(MGWR-BI)多源降水融合算法

基于经典的二源降水融合方法，建立融合降水、卫星降水和站点降水方程：

Pm=Ps+f(Pb)

(1)

其中Pb=Pg-Ps

式中：Pm为融合降水，mm/h；Ps为卫星降水，mm/h；Pg为站点降水；Pb为Pg减去Ps得到的偏差，mm/h。

由于CMORPH卫星降水的空间分辨率为8 km，时间分辨率为0.5 h，①将0.5 h的CMORPH卫星降水累加到1 h，并采用双线性插值法将空间分辨率为8 km的原始卫星降水数据(Po)降尺度到1 km(Ps)；②计算Pg与Ps之间的Pb，并运用MGWR-BI多源降水融合算法来估算每个网格的Pb；③根据式(1)将卫星降尺度降水扣除偏差得到融合降水数据，其时空分辨率为1 h、1 km，将其作为WRF-hydro模型的驱动数据[22]。

本文开发了一种基于MGWR-BI多源降水融合算法的“三步式”融合方法，每个时间步长的降水融合均通过该方法来实现。混合地理加权回归模型(MGWR)的方程可简写为

(2)

式中：Pb,i为第i个网格的降水偏差；n1为局部变量的个数，本文取值为5，局部变量主要有DEM、坡度、坡向、地表粗糙度、到海岸线的最短距离；l为全局变量；n2为全局变量的个数，本文取值为1；xik为第i个网格第k个局部变量；aik为第i个网格第k个局部变量对应的回归系数；xil为第i个网格第l个全局变量；ail为第i个网格第l个全局变量对应的回归系数；εi为随机误差。

MGWR方程求解的核心是回归系数aik、ail的计算，回归系数通过空间权重公式ωi来计算，由不同的空间权函数确定。常用的空间权函数主要有阈值法(TH)、距离反比法(IDS)、高斯函数法(GAU)、截尾函数法(BI)。本文选用截尾函数法(BI)进行空间权重的计算：

(3)

式中：dij为网格i到站点j之间的距离；b为描述权重与距离之间函数关系的非负衰减参数，依据AIC法则来确定。

2.2 WRF-Hydro模型

WRF-Hydro模型是由美国国家气象研究中心(NCAR)研发的[23-24]，基于WRF天气模式开发的陆气耦合模型系统[25-26]，主要包含大气模式、陆面模式(陆气参数化方案)、地下水模块、坡面河道汇流模块、水量管理方案5个部分[10]。

WRF-Hydro模型包含了降雨径流所涵盖的所有物理过程，产流由Noah-MP陆面模型计算，汇流部分的计算分为壤中流、坡面汇流、河道演算。土壤水分侧向流动按照壤中流的路径，饱和网格单元的下渗量由陆面模型计算得到；坡面汇流的计算采用完全不稳定、空间显式的扩散波方程。河道汇流的计算有马斯京根(Muskingum)方法、马斯京根-康吉(Muskingum-Cunge)方法、扩散波方法等计算方式可供选择。

WRF-Hydro模型能够根据不同的需求进行添加或改进，如耦合天气模式、陆面模式、不同水文模块等，因此是一个灵活的模型系统或称为可扩展的模型框架，能够实现网格与矢量的并行计算。计算效率的提高使得WRF-Hydro模型能够用于大规模的洪水预报[10]，不仅可以独立模拟用于多尺度洪水预报，还可以用于水资源管理、水质预测等多个方面[12]。目前，该模型在中国区域的应用研究较少，其在中国地区的研究具有更大的潜能。

WRF-Hydro模型输入数据的时空分辨率分别为1 h和1 km，同时模型输出数据的分辨率也设置为1 h和1 km，模型的时间计算步长设置为60 s。在计算过程中，WRF-Hydro模型首先通过陆面模型根据水量平衡进行产流计算，其次进行汇流计算。在产流计算过程中，陆面模型的参数化方案在WRF-hydro模型中已设置，仅根据需要选择参数化方法，无需进行参数率定。汇流计算涉及的参数可进行率定，也可通过WRFHydro_GIS_Pre-Processpr获取，因此本文未进行参数率定。

采用站点降水数据对融合降水进行精度验证，评价指标包括偏差(BIAS)、均方根误差(root mean square error，RMSE)、相关系数(R)。BIAS、RMSE用来描述站点降水数据与融合降水之间的误差与偏差，R用于描述站点降水数据与融合降水之间的拟合度。采用洪峰相对误差绝对值(absolute relative error of the flood peak，Dp)、洪量相对误差(Relative error of the flood volume，Dv)、确定性系数(Deterministic coefficient，Dc)3个常见的洪水预报评价指标评估WRF-Hydro模型的计算结果。

3 结 果 分 析

3.1 融合结果分析

由图3可知，5场洪水融合降水的BIAS范围为-0.65～0.35 mm/h，卫星降尺度降水的BIAS范围为-1.99～1.24 mm/h，融合降水的BIAS范围较小；除了2014083002场次洪水外，其他4场洪水融合降水BIAS的绝对值小于卫星降尺度降水BIAS的绝对值。5场洪水卫星降尺度降水的RMSE均大于融合降水的RMSE，二者的范围分别为3.29～7.85 mm/h、1.61～6.53 mm/h。5场洪水融合降水的R均大于0.65，卫星降尺度降水的R均小于0.15。通过偏差、均方根误差、相关系数指标的对比可见，相对于卫星降尺度降水，融合降水的精度有了较大的提高。

图3 各场次洪水实测站点降水与卫星降尺度降水、融合降水散点密度Fig.3 Scatter density maps for precipitation of observation station, downscaling precipitation of satellite and merged precipitation

3.2 WRF-Hydro模型计算结果对比分析

由于WRF-Hydro模型所有的输入数据都是网格化的数据，站点数据若作为输入需要进行插值，在插值过程中会产生误差，对模型计算结果造成一定的影响，导致精度降低。通过对融合降水数据的验证与分析可知，卫星降水与站点降水通过融合之后得到的融合降水数据精度有了较大幅度的提高。本文仅对原始卫星降尺度降水和融合降水数据驱动下WRF-Hydro模型计算的流量与实测流量进行对比。

两河口水文站5场洪水WRF-Hydro模型模拟的洪峰流量均小于实测洪峰流量，卫星降尺度降水模拟的洪峰流量更小(图4)。卫星降尺度降水模拟的洪峰相对误差绝对值均大于0.42，融合降水模拟的洪峰相对误差绝对值均小于0.15，融合降水数据驱动WRF-Hydro模型显著提高了洪峰模拟的精度(表2)。卫星降尺度降水模拟的洪量相对误差在-0.64～-0.47之间，融合降水模拟的洪量误差在-0.3～0.13之间，每场洪水的洪量相对误差均有所减小(表2)。相对于黑色的过程线，5场洪水红色的流量过程线更接近灰色填充区，说明融合降水作为驱动数据时，WRF-Hydro模型模拟的流量过程与实测流量过程拟合得更好(图4)。由表2可知，每场洪水融合降水数据驱动WRF-Hydro模型计算的确定性系数均大于卫星降尺度降水数据驱动下模型计算的确定性系数。通过与实测流量过程线对比，融合降水模拟精度高于卫星降尺度降水，WRF-Hydro模型在融合降水作为输入的条件下，模型计算精度的提高间接证明了融合方法的有效性。

表2 5场洪水模拟结果

图4 两河口水文站流量过程线对比Fig.4 Comparison of observed hydrographs and modeled hydrographs at Lianghekou station

3.3 WRF-Hydro模型验证

为了进一步评估WRF-Hydro模型在空间上的预报能力及空间计算的合理性，对比分析了在融合降水数据驱动下，子午河流域5场洪水19条水系洪峰流量的空间分布和峰现时间的空间分布(图5)。

子午河流域19条水系可分为一级水系(1、2、3、4、5、6、7、8、15、16)、二级水系(9、10、11、12、18)、三级水系(13、14、17)和出口断面(19)。图5(a)和图5(d)两场洪水的量级较小，两场洪水一级支流的洪峰流量均小于75 m3/s；图5(b)(c)(e)3场洪水的洪峰流量在第19号水系均大于1 200 m3/s，3场洪水一级支流的洪峰流量均小于220 m3/s。5场洪水二级支流的洪峰流量均大于对应的一级支流的洪峰流量，三级支流的洪峰流量均大于对应的二级支流的洪峰流量，并在流域出口断面达到最大值，分别为407.3 m3/s、1 287.7 m3/s、1 377.0 m3/s、361.6 m3/s、1 236.3 m3/s(图5)。出口断面的洪峰流量大于三级支流，三级支流大于二级支流，二级支流大于一级支流，洪峰流量的空间分布证明了WRF-Hydro模型计算结果的合理性。

图5 5场洪水每条水系洪峰流量和峰现时间空间分布Fig.5 Spatial distribution of simulated discharge and flood peak time for five flood events

5场洪水的峰现时间分别为2013-05-29T20：00、2013-07-18T16：00、2013-07-22T21：00、2014-08-31T03：00、2014-09-10T21：00，历时分别为184 h、160 h、131 h、177 h、469 h(表1)。将5场洪水的历时根据时间进行编号，5场洪水的模拟洪峰分别出现在第37个小时、第81个小时、第20个小时、第26个小时、第93个小时。

5场洪水一级水系的峰现时间分别为28～30 h、74～77 h、6～13 h、10～15 h、51～79 h，二级水系峰现时间分别为31～33 h、77～80 h、14～17 h、15～20 h、82～87 h，三级水系峰现时刻分别为35～36 h、80 h、17～19 h、23～24 h、87～92 h，出口断面峰现时间分别为37 h、81 h、20 h、26 h、93 h。此外，2014090701场次洪水的峰现时间从一级水系到三级水系变化较大，可能是由于该场洪水对应了多个降水事件，出现了多个洪峰(图5(j))。5场洪水的峰现时间从上游到下游呈现出一致性规律，峰现时间依次为一级水系、二级水系、三级水系、出口断面，即上游的水系先出现洪峰，下游的水系后出现洪峰，流域出口断面最后出现洪峰(图5)。

由5场洪水洪峰流量和峰现时间的空间分布可知，洪峰流量和峰现时间从一级水系到三级水系逐渐增加，在流域出口断面达到最大值，流域出口断面的流量主要来自椒溪河的贡献(图1和图5)；5场洪水的峰现时间从上游到干流再到下游逐渐增加，流域出口断面的洪峰最晚出现。以上充分说明了WRF-Hydro模型在产流、汇流计算方面的合理性。

5场洪水降水事件中，2013052808场次和2014083003场次洪水事件的融合降水总量相对较小，小于100 mm的占流域面积的75%以上(图6(a)和(d))；平均地表径流相对最小，小于25 mm的占流域面积的79%以上(图7(a)和(d))；降水总量和平均地表径流从侧面证明了这两场洪水事件的流量较小，实测流量均小于500 m3/s(图4(a)和(d))。2014090701场次洪水事件整体上融合降水总量最大(图6(e))，这是因为该场次洪水事件历时最长且包含了多个降水事件(图4)；融合降水总量在200～460 mm之间，大于400 mm、300～400 mm、200～300 mm的分别占流域总面积的24.76%、42.42%、32.81%；平均地表径流在20～80 mm之间，大于50 mm的占流域总面积的10.41%，平均地表径流较大的区域较少(图7)。2013071508场次和2013072202场次洪水的融合降水总量和平均地表径流均较大(图6(b)(c)、图7(b)(c))，其洪峰流量均超过了1 000 m3/s(图4(b)(c))。总之，对于每场洪水事件而言，融合降水总量与平均地表径流的空间分布高度一致。即融合降水总量较大的网格其平均地表径流也较大；融合降水总量较小的网格，其平均地表径流也较小；从侧面证明了模型在产流计算方面的合理性和正确性。

图6 各场次洪水融合降水总量空间分布Fig.6 Spatial distribution of total merged precipitation for each flood event

图7 各场次洪水降水模拟平均地表径流空间分布Fig.7 Spatial distribution of simulated average surface runoff from the precipitation of each flood event

4 结 语

本文基于混合地理加权回归截尾函数算法(MGWR-BI)开发了一种“三步式”降水融合方法，将站点实测降水与CMORPH卫星降尺度降水进行融合，以卫星降尺度降水和融合降水数据作为驱动数据，探讨了融合降水数据对中小河流洪水预报结果的影响。基于WRF-Hydro模型对子午河流域5场洪水进行模拟计算，通过流量过程对比和空间分析表明：WRF-Hydro模型在中小河流洪水预报方面具有合理性和适用性。

MGWR-BI多源降水融合算法实现了卫星降水数据的降尺度与融合，并且每场洪水对应的融合降水数据精度均高于原始卫星降尺度降水，表明MGWR-BI多源降水融合算法具有一定的稳健性，并且证明数据融合技术能够提高降水的质量和时空分辨率。

全流域19条水系5场洪水的洪峰流量和峰现时间的空间分布规律相似，其从大到小排列顺序为：流域出口、三级水系、二级水系、一级水系，证明WRF-Hydro模型在产流和汇流计算方面的合理性，同时说明该模型在空间计算方面具有无偏性的特点。WRF-Hydro模型能够预报某一流域尺度上任意一条水系的流量，将网格化的产流与矢量路径的汇流相结合，极大地提高了计算效率，这一特点体现了其在实时洪水预报中的应用潜力。

数据融合技术不仅能够获取高精度的降水数据，并且能够提高中小河流洪水预报的精度。因此，将数据融合技术引入WRF-Hydro模型，不仅为中小河流洪水预报提供了新的途径，而且为无资料地区的中小河流洪水预报研究提供了新思路，对中小河流的防灾减灾具有重要的实用价值和经济社会效益。