基于径流-蒸散发-蓄水量协同率定的平原区水文模拟研究

田益民, 鲍振鑫, 宋晓猛, 莫昱晨, 王国庆, 刘翠善

(1.贵州省水利水电勘测设计研究院有限公司,贵州 贵阳 550001; 2.中国矿业大学 资源与地球科学学院,江苏 徐州 221008; 3.南京水利科学研究院 水文水资源与水利工程科学国家重点实验室,江苏 南京 210029;4.水利部应对气候变化研究中心,江苏 南京 210029)

平原区是人口、城市、耕地等的密集区,水文模拟及预报对区域防汛抗旱、水资源规划与管理、农业生产等具有重要的支撑作用。以山坡水文学为基础构建的水文模型是对流域水文循环过程的一种数学概化,是认识水文基本规律,重现和预测流域水文过程的重要工具[1-3]。水文模型参数对水文模拟精度有重要的影响,一般利用流域出口断面的实测径流资料优化率定求解。但是,平原区无统一的流域出口断面,确定模型参数以模拟水文过程是水文模拟领域的关键且薄弱环节,是当前亟须解决的一项国际前沿和热点问题[4-7]。随着对地立体观测技术的不断进步,水文学家逐渐构建了全球蒸散发和陆地水储量等水文要素数据集,为水文模型参数率定与模拟精度评价提供了新的资料[8-12]。随着对水文基本规律认识的不断深入和科学计算能力的快速提高,水文模型的结构越来越复杂,基于物理机制分布式模型的应用越来越广泛[13-14],这些都为平原区水文模拟的发展提供了新方向[15]。

水文学家关于模型参数率定开展了大量研究,提出了单纯形法、Rosenbrock法、SCE-UA法等优化求解方法[16-22]。平原区无封闭流域出口断面,以往研究主要基于区域地下水位和土壤含水量等实测数据,率定陆面水文模型参数模拟垂向的降水和入渗过程,该方法侧重于评估区域水循环储量的模拟精度,对径流和蒸散发等通量的评估较为薄弱[23-24]。山丘区的水文模型参数率定重点将径流模拟的优劣作为主要评价目标,未考虑蒸散发和蓄水量等水文要素。但是,大部分降水消耗于蒸散发,形成径流的比例较少[25],在年尺度上,天然条件下区域水储量的变化不显著[26]。例如,黄河、海河等流域的多年平均降水径流系数均小于0.4[27-28]。由于模型参数对水分在流域中存储和耗散的分配过程具有不同的影响,仅根据某个水文要素优选的模型参数用以模拟其他水文要素时会产生较大的误差。PENG B等[29]在率定径流的基础上,将Grace卫星反演的流域水储量数据纳入模型校验标准,提高了水文模拟的精度。由于模型参数的相互作用、参数优化空间的高维脊、水文要素对模型参数的敏感性差异等问题,模型参数存在显著的“异参同效”现象,优化算法无法得到模型参数的唯一最优解,从而引起水文模拟的不确定性[30-35]。因此,对于径流模拟结果较好的模型参数组,在模拟蒸散发量和蓄水量时可能会产生较大的误差。对于地下水位和土壤含水量模拟结果都较好的模型参数组,在模拟径流和蒸散发量时也可能会产生较大的误差。如何同时考虑径流、蒸散发、蓄水量等水文循环通量和储量以提高模拟精度,缩减模型参数的待选样本空间,减弱参数率定的不确定性,是水文模拟研究的重要突破方向。

本文以黄淮海平原为研究对象,基于区域水循环通量和储量的实测与遥感反演数据,利用分布式可变下渗容量(Variable Infiltration Capacity,VIC)模型,构建基于径流、蒸散发和蓄水量的协同率定模型参数率定方法;采用随机森林等人工智能算法,分析模型参数的敏感性;利用广义似然不确定性估计(Generalised Likelihood Uncertainty Estimation,GLUE)方法,评估协同率定框架下的径流模拟不确定性,进而推动水文模拟技术的进步。

1 研究区域和数据资料

1.1 研究区域

黄淮海平原是我国三大平原之一,如图1所示,位于东经113.0°～121.5°、北纬30.0°～40.5°,涵盖北京、天津、河北、山东、河南、安徽及江苏等地区。

图1 黄淮海平原

黄淮海平原跨越10个水资源二级区,面积约为31万km2,占全国国土面积的3.10%。黄淮海平原地势平缓,是我国重要的农业生产基地,其中耕地面积和粮食播种面积分别约占全国的24.96%和33.36%[36]。主要作物为小麦和玉米,农业种植制度以冬小麦与夏玉米轮作为主。该区域位于温带大陆季风气候带,具有雨热同期的特点,降水四季分明,年降雨量为500～900 mm,年降水量不均匀,多集中于夏季。黄淮海平原水资源短缺问题突出,每亩水资源量仅为全国平均水平的12.50%。

1.2 数据资料

本研究的计算时段为2003年1月—2016年12月,采用的数据包含最高气温、最低气温、降水、风速、径流、蒸散发、蓄水量、植被覆盖、土壤质地、数字高程等要素。区域内的逐日最高气温、最低气温、降水、风速等气象要素数据源自于中国气象数据共享网(https://data.cma.cn/)。平原区无封闭的流域出口,受上游来水和人类取用水等多种因素影响,水文站的实测径流不能反映研究区域的天然径流量,无法用于水文模型的率定和校验。平原区天然径流还原是一项繁巨的工作,目前可利用的相对权威的数据主要为水利部主持的全国水资源调查评价数据和1997年以来发布的《中国水资源公报》的相关数据。该数据描述的是各个水资源分区的本地天然年径流量,考虑了取用水和相关水利工程设施对径流量的影响,可用来评价水文模型模拟年径流量的精度,用其来表征黄淮海平原及其内部10个水资源二级区年尺度的天然径流量。蒸散发数据由英国布里斯托大学地理科学学院水文气象系发布的全球Gleam v3.3遥感数据产品估算[37]。该产品由多种卫星遥感数据同化得到,空间分辨率为0.25°×0.25°,时间分辨率为月,包括a、b两套数据集,本研究使用的是Gleam v3.3b产品。相关研究表明,Gleam v3.3蒸散发产品的精度较高,在黄淮海流域具有较好的适用性[38]。在多年水循环条件下,流域的蓄水量一般呈现稳定状态,本研究利用实测降水和还原的天然径流数据,在年尺度上对10个水资源二级区的Gleam蒸散发数据做偏差校正,对每年蒸发的月数据都乘以校正后相对应的比例因子,以更好地反映黄淮海平原的实际蒸散发情况。区域蓄水量由Grace卫星反演的逐月陆地水储量数据估算,数据空间分辨率为1°×1°,数据来源于JPL(https://grace.jpl.nasa.gov/data/get-data/jpl-global-mascons/)[39]。植被覆盖数据采用的是马里兰大学发布的全球1 km×1 km土地覆盖数据集。土壤数据源自寒区旱区科学中心(http://westdc.westgis.ac.cn/)基于世界土壤数据库整理的中国土壤数据集,分辨率为5 km×5 km。数字高程数据源自地理空间数据云(http://www.gscloud.cn/)SRTMDEM数字高程数据,分辨率为90 m×90 m。

研究建立的VIC模型网格空间分辨率为50 km×50 km,对于不同的数据类型有不同的处理原则。气象驱动数据:由45个气象站点的降水、风速、最高气温、最低气温4个气象要素根据空间插值的方法获得网格型逐日气象数据。植被数据:统计每个网格内所有的土地覆盖类型,计算每个网格内每种土地覆盖类型的占比,同时加上每种土地覆盖类型的逐月植被叶面积指数等理化性质。土壤数据:统计每个网格中该土壤占比最大的土壤作为网格的土壤类型,同时附上代表该土壤类型的理化性质。

2 研究方法

2.1 VIC模型简介

VIC模型是美国华盛顿大学、加利福尼亚大学伯克利分校以及普林斯顿大学共同研发的大尺度水文模型,是一个基于空间网格化的半分布式水文模型[40-41]。VIC模型可同时进行陆-气间能量平衡和水量平衡的模拟,弥补了传统水文模型不能对热量过程模拟的缺陷。同时VIC模型能够描述地表主要的水文气象过程,如:地表径流、基流、土壤层蒸发、植被蒸散发、地表截留蒸发等。它考虑了蓄满产流和超渗产流两种机制,分析了土壤性质的次网格化非均匀性对产流的影响,用下渗能力面积分配曲线和土壤水容量面积分配曲线表示土壤不均匀性对产流的影响。VIC模型将土壤分为3层,其中上两层产生地表径流,下层产生基流。模型采用半分布式降雨-径流ARNO模型方案模拟基流过程,在某一阈值下基流是线性消退的,当土壤含水量高于这一阈值时,基流是非线性的。

VIC模型有6个参数需要率定,包括可变下渗曲线指数(B)、非线性基流开始时基流值与最大基流的比值(DS)、最大基流流速(Dm)、下层土壤最大含水量比例系数(WS)、第2层土壤厚度(d2)和第3层土壤厚度(d3),这些参数是影响VIC模型陆面模式中产流的参数。将黄淮海平原划分成253个网格,如图1所示,在10个水资源二级区中分别率定模型参数。其中2003—2012年的资料用于模型率定,2013—2016年的资料用于模型检验。选用均方根误差(RMSE)、相关系数(r)以及相对偏差(BIAS)3个指标评价水文模拟的效果[42]。由于本研究的时间尺度随评价的水文要素的不同而不同,相关评价指标主要为不同水文要素评价模型模拟的误差。构建VIC模型的时间精度与输入的气象数据的一致,时间尺度为日。由于平原区还原的天然径流量的时间尺度为年,则统计VIC模型陆面模式输出的逐年总径流与还原的天然径流做对比。Glean v3.3b数据的时间尺度为月,则统计VIC模型输出的逐月总蒸发和Glean v3.3b的蒸发做对比。模型的蓄水量由VIC模型3层土壤含水量之和来表征,Grace卫星水储量作为观测的蓄水量时间尺度为月,则统计模型逐月总土壤含水量的平均值与Grace卫星水储量做比较。在月、年等长时间尺度上,从区域水量平衡的角度评价径流总量、蒸散发量和蓄水量的模拟精度,成果可为水资源评价及相关模型参数率定提供一定参考。

2.2 基于协同率定的水文模拟框架

天然水循环过程中,降水主要影响径流、蒸散发和蓄水量的变化,其水量平衡方程为:

P=R+E+ΔS。

(1)

式中P、R、E、ΔS分别为降水量、径流深、蒸散发量、蓄水量变化。

完整的水文模拟应同时考虑径流、蒸散发和蓄水量3个要素的模拟精度。根据实测和遥感反演同化数据,构建了基于径流、蒸散发和蓄水量的水文模型协同参数率定方法与模拟框架。为了与传统方法对比,将以往只将径流模拟优劣作为模型参数率定评价指标的方法称为单一径流率定法,将本研究同时考虑径流、蒸散发和蓄水量3个要素模拟精度的方法称为协同率定法。径流模拟精度由VIC模型模拟的总径流量与还原的天然径流量对比评价。蒸散发模拟精度由VIC模型模拟的各组分蒸散发之和与校正的Glean v3.3b产品对比评价。对于蓄水量模拟精度,由于Grace卫星反演的是水储量与卫星发射时刻区域水储量的相对值,而不是某一时刻水储量的绝对值。考虑到卫星反演水储量和模型模拟值的量级和单位不统一,在评价蓄水量模拟精度时,需对数据作标准化处理,利用相关系数来评价。综合考虑以上3个水文要素的评价指标,首先利用蒙特卡罗法生成多组参数,将径流评估的均方根误差、蒸散发评估的均方根误差和蓄水量相关系数进行归一化处理,得到每组3种评价指标的权重;然后,选择归一化之后的径流均方根误差加上蒸散发均方根误差减去蓄水量相关系数,得到的最小参数组作为模型最后的目标参数:

Omin=ZSRMSER+ZSRMSEE-ZSRS。

(2)

式中:Omin为模型最优参数率定目标;ZSRMSER为评估径流的均方根误差归一化后的指标;ZSRMSEE为评估蒸散发量的均方根误差归一化后的指标;ZSRS为评估蓄水量相关系数归一化后的指标

2.3 模型参数敏感性与不确定性分析

模型参数的取值对水文模拟精度有显著的影响。利用基于随机森林算法构建的模型参数与径流模拟均方根误差的统计模型来评估VIC模型的6个参数的敏感性。随机森林算法是在Bagging集成学习理论基础上发展起来的一种机器学习方法,被广泛应用于回归问题和分类问题,也可通过分析若干自变量对因变量的影响程度来研究特征选择问题[43-44]。该方法主要采用重抽样方法,从原始样本中随机重抽取多个样本,对每个样本利用决策树方法建立模型,然后将所有样本的计算结果集合为最终结果。分析VIC模型6个参数的敏感性时,首先,根据参数的物理意义和取值范围,利用蒙特卡罗方法生成多个参数组样本;然后,分别利用每一组参数,模拟径流过程并统计其均方根误差,再将径流均方根误差作为因变量,模型参数作为自变量,利用随机森林算法构建回归模型;最后,利用随机森林算法中各个自变量的重要程度指标表征VIC模型6个参数的敏感性。

GLUE法是一种基于异参同效和贝叶斯理论的不确定性分析方法[45-46]。该方法是在敏感性分析基础上演变而来的一种评价不同模型参数取值引起的径流模拟不确定性的估计方法。利用GLUE方法评估传统单一径流率定法和协同率定法模拟径流的不确定性,对比分析协同率定能否减弱径流模拟的不确定性。根据水文要素模拟评价指标,设定参数筛选阈值,根据模型参数敏感性分析中的蒙特卡罗随机模拟结果,选取有效的参数样本。以往只将径流作为单一评价指标,而本研究同时考虑径流、蒸散发和蓄水量等多个评价指标。首先,选取满足径流和蒸散发均方根误差按从小到大排列的前20%的参数组和蓄水量相关系数的前20%的参数组[35]作为有效参数样本;然后,将水文要素的评价指标作为每个参数样本的权重,并对其归一化后作为每个参数样本的概率;最后,将每个时刻所有有效参数样本模拟的水文要素值从小到大排序,统计其累积概率,分析水文模拟的置信区间,评估其不确定性。

3 结果与分析

3.1 模型参数敏感性分析

3.1.1 模型参数与径流模拟评价指标相关性

对于黄淮海平原及其10个水资源二级区,VIC模型参数与径流模拟3个评价指标的相关性如图2所示。由于模型率定参数与评价指标之间存在非线性关系,当率定参数超过最优解时,指标评价效果会随之降低。为了显示出率定参数与评价指标的相关性,图2中的模型率定参数是相对于协同率定法估算最优解的绝对偏差。结果显示:d2与模型3个评价指标的相关性较显著,与均方根误差和相对偏差的相关性最明显,其中大部分研究区的d2与均方根误差的相关性系数超过0.50,部分区域的相关性系数超过0.70;d2与相对偏差的相关性系数基本稳定在0.65以上,与相关系数的相关性系数基本在0.40上下波动。可变下渗曲线指数B与3个评价指标呈中低度相关性,其相关系数的绝对值为0.1～0.4。其余4个参数与模拟径流评价指标的相关性较弱,相关系数绝对值基本都小于0.1。总体而言,VIC模型参数d2对径流模拟结果具有重要的影响,其次是参数B,径流模拟结果对其余4个参数的敏感性较弱。

图2 模型参数与径流评价指标的相关性

将模拟径流的均方根误差作为因变量,6个模型参数作为自变量,基于随机森林算法,构建模型参数与模拟均方根误差的回归模型。利用随机森林算法中各个自变量的重要程度指标表征VIC模型6个参数的敏感性,其结果如图3所示。

图3 模型参数重要程度指标

由图3可知,d2最重要(重要程度指标为84%),其次是B(重要程度指标为80%),参数DS、Dm、WS、d3的重要性较小,其重要程度指标分别为3.0%、4.0%、0.9%、0.1%。在VIC模型中,由于d2表征网格的蓄水能力,其直接影响下渗能力,进而控制地表径流产流量的大小,所以模型对d2的敏感性最高。可变下渗曲线指数B反映区域下渗能力空间分布的不均匀性,影响模型的产流过程,所以模拟结果对参数B也有一定的敏感性。而其余4个参数主要控制基流的产生,对地表径流的影响较弱。该结果与孙若辰[47]得出的相关结论相符。

3.1.2 模型参数优选空间分析

VIC模型对6个参数中d2最为敏感。以受人类活动影响较小的淮河上游王家坝以上的水资源二级区作为研究对象,分析d2的不同取值对径流量、蒸散发量和蓄水量模拟评价指标的影响,结果如图4所示。由图4可见:只考虑径流评价指标时,当均方根误差较小时,d2的取值范围为0.15～0.40 m;只考虑蒸散发时,d2在0.35 m附近;蓄水量较难准确模拟,其参数的取值范围较大。此外,分析了每组参数样本对径流量、蒸散发量和蓄水量的协同影响。以图4(a)和4(b)为例,当d2小于0.35 m时,虽然径流量、蒸散发量和蓄水量独立评估时都能获得较好的模拟效果,但是没有一个参数样本能使3个水文要素同时实现较好的模拟效果。也就是说,当d2小于0.35 m时,对于VIC模型的其余5个参数与径流量、蒸散发量和蓄水量的参数最优解不一致。因此,同时考虑3个水文要素的模拟效果能够缩小模型参数的优选空间,从而减少模型参数优选的不确定性。

图4 径流、蒸散发、蓄水量协同率定下的参数d2优选空间

3.2 径流量、蒸散发量、蓄水量模拟

基于构建的径流-蒸散发-蓄水量的协同率定方法,确定的模型参数率定结果见表1。利用VIC模型模拟的2003—2016年黄淮海平原径流量、蒸散发量和蓄水量过程,及其与传统单一径流率定法的对比结果如图5所示,相关评价指标见表2。其中,2003—2012年数据用于模型参数率定,2013—2016年数据用于模型校验。结果表明,VIC模型能够较好地模拟黄淮海平原水文过程。

表1 黄淮海平原协同率定参数率定结果

表2 黄淮海平原径流、蒸散发和蓄水量模拟评价结果

图5 黄淮海平原径流、蒸散发和蓄水量过程模拟

从径流和蒸散发的模拟情况看,无论是率定期还是检验期,基于3种水文要素的协同率定方法的VIC模型都能较好地模拟研究区径流与蒸散发的变化过程,其相关系数均在0.90以上,相对偏差都较小。对于蓄水量,协同率定方法在率定期的相关系数超过0.60,但是检验期的模拟精度较差,相关系数仅为0.11。在天然状态下,流域多年水储量的补排情况呈稳定状态[48]。但是,受人类活动的影响,近几十年,黄淮海平原地下水超采严重,严重干扰了区域水量平衡状态,流域蓄水量损耗量持续大于补给量,地下水位持续下降,蓄水量显著减少[49-50]。由于VIC模型未考虑人类活动对区域水平衡的影响,利用2003—2012年资料率定参数模拟的2013—2016年蓄水量比实测值偏大。单一径流率定法只能较好地模拟径流过程,而蒸散发和蓄水量的模拟效果较差。两种参数率定方法相比,模拟径流的精度相似,但是协同率定方法模拟的蒸散发和蓄水量精度高于单一径流率定法的模拟精度。例如:在率定期,协同率定方法模拟的径流量相对偏差、蒸散发相对偏差和蓄水量相关系数分别为-1.11%,-0.47%和0.61,而单一径流率定法计算的3个指标分别为0.66%、-1.86%和0.44;在检验期,协同率定方法模拟的3种水文要素结果都要好于单一径流率定法的。因此,基于径流-蒸散发-蓄水量的协同率定方法在不降低径流模拟精度的同时,能够更准确地模拟蒸散发和蓄水量过程,适用于水文循环多要素系统模拟。

为更好地反映黄淮海平原水文过程的空间差异性,分别在10个水资源二级区,利用协同率定方法率定VIC模型参数,验证径流量、蒸散发量和蓄水量等水文要素的模拟性能,评价结果如图6所示。10个水资源二级区独立率定模拟参数的模拟结果与整个黄淮海平原统一率定的结果相似。从径流量和蒸散发量的评价指标来看,模型在率定期和检验期的模拟效果都较好,相关系数都在0.90左右,相对偏差基本都小于10%;从蓄水量的模拟结果来看,率定期的相关系数为0.60左右,但是检验期的相关系数为0.40左右。这是因为黄淮海平原水储量处于长期下降的状态[49],导致在检验期蓄水量比实测值偏大,其相关性减弱。总体上VIC模型能够较好地模拟黄淮海平原10个水资源二级区的径流量、蒸散发量和蓄水量等水文要素,适用于该区域的水文模拟。

图6 黄淮海平原10个水资源二级区水文要素模拟评价指标

3.3 水文模拟不确定性评估

3.3.1 协同率定参数优选

由于模型参数之间的相互作用,水文模型存在异参同效现象,即:不同的模型参数组能够得到相似的径流模拟效果,但有可能模拟的蒸散发和蓄水量差别较大。本研究分析了径流均方根误差、蒸散发均方根误差和蓄水量相关系数的关系,结果如图7所示。从图7(a)可以看出,满足径流模拟较好的参数有很多组,但大部分参数对蒸散发和蓄水量的模拟较差,只有少部分能够较好地同时模拟这3个水文要素。同时,进一步分析了仅考虑径流均方根误差按由小到大排名前10的参数模拟蒸散发和蓄水量的情况,结果如图7(b)所示。从图7(b)可看出,虽然前10组参数样本的径流均方根误差都很低,但其蒸散发均方根误差和蓄水量却不一定都模拟得好。例如:径流模拟排名第一的参数组中,径流均方根误差取值最小,蒸散发均方根误差为15.4,而蓄水量相关系数只有0.4;径流模拟排名第三的参数组中(其径流模拟的均方根误差与第一组的相比变化很小),蒸散发均方根误差仅为8.0,蓄水量相关系数为0.6。因此,协同率定的参数率定方法能够减小模型参数的可能优选空间,进而缩小径流模拟的置信区间,在几乎不降低径流模拟精度的情况下,大幅度提高蒸散发和蓄水量的模拟精度。

3.3.2 径流模拟的不确定性

利用GLUE方法评估了单一径流率定和协同率定两种参数率定方法模拟径流的不确定性,结果见图8和表3。两种参数率定方法模拟径流的90%置信区间都能覆盖观测径流。相比而言,协同率定方法的90%置信区间比单一径流率定的90%置信区间平均减少26%。协同率定法模拟径流的50%估计值和观测值的拟合程度较好,其相对偏差以及均方根误差的评价效果均好于单一径流率定的结果。可见,通过多目标的参数率定,减小了径流模拟的不确定性。

表3 两种参数率定方法模拟径流的不确定性对比

图8 参数率定方法在黄淮海平原模拟径流的不确定性

4 结论

根据水资源评价还原的天然径流、遥感同化的蒸散发和Grace卫星反演的水储量等区域水循环通量和储量多源要素数据集,构建了基于径流、蒸散发和蓄水量的协同率定方法,评估了VIC模型在黄淮海平原水文模拟中的适用性,在模型参数敏感性分析的基础上,研究了协同率定方法对参数优选空间、异参同效现象和径流模拟不确定性等方面的影响,主要结论如下:

1)VIC模型能够较好地模拟黄淮海平原的水文要素,模拟结果与实测值相比,模拟径流和蒸散发在率定期的相关系数大于0.9,蓄水量的相关系数大于0.6。

2)对于VIC模型的6个参数,d2对模拟结果的影响最显著,其次是B,模拟结果对其余参数的敏感性较弱。

3)与单一径流率定方法相比,基于径流-蒸散发-蓄水量的协同率定方法缩小了模型参数的优选空间,其径流模拟的90%置信区间比单一径流率定的90%置信区间大幅度减小,在几乎不降低径流模拟精度的情况下大幅度提高了蒸散发和蓄水量的模拟精度。

本研究提出的协同参数率定方法可应用于其他平原区或者山丘区进一步检验其效果和适用性。随着更多水文要素数据集的不断完善,可以在此基础上,深入研究更多率定目标对水文模型不确定性的影响,更好地评估水文模型对水循环系统各组分和分量的模拟性能,为流域水文模拟技术进步提供更有力的科学支撑。