基于多流域同步率定的HBV模型参数移植方法研究

2021-03-05 08:18盛奕华衣学军黄迎春李京兵
中国农村水利水电 2021年2期
关键词:水文平均值流域

盛奕华,衣学军,黄迎春,李京兵,姚 成,李 昊

(1.河海大学水文水资源学院,南京 210098;2.烟台市水文局,山东 烟台 264000;3.安徽省水文局,合肥 230022;4.烟台市城市水源工程运行维护中心,山东 烟台 264000)

流域水文模型作为一种研究流域内复杂水文现象的有效工具,在洪水预报、水资源规划和管理等工作中发挥着重要作用[1]。水文模型的参数主要依赖于水文气象观测资料,通过率定的方式来确定,因此无资料地区水文模型参数的确定问题是模型应用的主要难点。2003年国际水文科学协会(IAHS)就此提出了无资料流域水文预测的PUB计划(Prediction in Ungauged Basins)[2-4],该计划通过10多年的研究取得了丰硕成果,形成了无资料地区水文模拟和洪水预报的方法集。目前国内外常用的无资料地区参数移植方法主要基于区域化方法,常用的方法有:多元回归法[5]、空间临近法[5-7]、流域物理特性相似法[5,8,9]以及全局平均法[5,10]。

水文降雨-径流过程容易受到流域气象条件的影响,研究表明同一流域在干旱、湿润等不同气象条件下,水文响应可能存在一定的差异性。水文模型参数的率定依赖于所选取的历史数据,若率定期的数据代表性不足,则容易出现所率定的模型参数与率定期特殊的气象条件[11]“过度拟合”的现象。这种气象条件的不确定性增加了参数的不确定性[12],使得在无资料地区模型参数移植的稳定性不高。目前的区域化方法若只考虑两个流域之间的参数移植,难以有效解决由于气象因素导致的水文模型“过度拟合”问题,而全局平均法[13]虽然考虑了多个流域,但该方法只是将模型参数的平均值直接进行移植,取得的模拟效果较差。

基于以上问题,本文开展基于多目标优化函数[14-16]的多流域HBV模型同步率定方法的研究,以降低上述参数不确定性的影响,通过获取多个流域能够同时适用的共享参数,提高共享参数在无资料地区的可移植性。本文以概念性水文模型HBV模型[17]为例,对安徽屯溪流域进行模型率定和参数移植,通过与单个流域之间的模型参数移植效果进行比较,分析基于多流域同步率定方法的参数移植效果。

1 研究区域与数据

选取安徽省南部山区的屯溪流域作为本文的研究区,该流域临近中国东南沿海,居新安江上游,处于亚热带季风气候区,季节分明,气候温和,多年平均气温约17 ℃。其多年平均降雨量约为1 800 mm,降雨主要集中在6至9月的汛期。流域面积为2 692 km2,地势西北高东南低,在屯溪水文站上游还嵌套有6个水文站,分别是榆村水文站、流口水文站、新亭水文站、呈村水文站、万安水文站和月潭水文站(图1)。

选取屯溪流域及其嵌套子流域的径流资料用于HBV模型参数的率定和验证,各流域的地貌特征及资料系列见表1。用于日径流模拟的资料序列较长,其中丰、枯水年份兼具,具有一定的代表性。表中流域地貌特征等数据是基于数字高程模型(DEM)提取得到的。

表1 屯溪流域及其嵌套子流域概况Tab.1 Characteristics of Tunxi catchment and the nested sub-catchment

2 模型与方法

2.1 HBV模型

HBV(Hydrologiska Byärns Vattenbalansavdelning)模型是由瑞典气象水文研究所(SMHI)开发的概念性水文模型,该模型由于结构简单,需要通过历史资料率定的参数较少且应用简便,在瑞典、美国等国家得到了广泛应用,取得了良好的应用效果[17]。该模型依据水量平衡原理进行建模,首先通过度日因子法计算融积雪量,并通过实际土壤蓄水量函数进行地下水补给量和流域实际蒸散发量的计算,通过线性水库方程将径流分为地表径流、壤中流和地下径流3种,最后基于单位线方法的三角权重函数进行出口流量过程的演算[18,19]。HBV模型结构主要包含融雪计算、土壤水分计算和流域产汇流3个部分[20]。由于本文研究的区域位于亚热带地区,降雪比重低,因此在模型中不考虑融雪计算部分。

本研究所构建的HBV模型需要率定的参数有7个:FC、Beta、L、K0、K1、KD和K2,各个参数的物理意义[21,22]及取值范围如表2所示。

表2 参数取值范围Tab.2 Description and threshold values of the calibrated parameters

2.2 基于ROPE方法的模型参数优选

ROPE(Robust Parameter Estimation algorithm)[14,15]是由Bardossy和Singh提出的一种高效稳健的模型参数自动优化算法。该方法基于深度函数的概念,把模型参数组作为一个多维度空间,通过蒙特卡洛采样不断迭代生成高维空间内部的新参数组,从而获得预定数量的、具有稳定模拟效果的参数组。ROPE算法采用的是半空间深度函数判断新参数组的深度值。半空间深度函数可计算数据点相对于数据集的半空间深度(HD),其定义如下:

设x是d维空间中的一点,F是d维空间上的一个分布函数,那么x相对于F的半空间深度就是包含x的半空间上的最小概率。

HD(x,F)=inf{F(H),x∈H},x∈Rd

(1)

式中:H是包含x的封闭半空间;inf为取下界。

根据半空间深度函数的性质,深度值大于1的点是位于d维空间的分布函数内部的,因此通过随机采样和深度值计算,可以生成d维空间内部的新点距。

ROPE方法的主要步骤如下:

(1)选取需要进行率定的模型参数,根据参考文献等资料[2,14-16],确定模型参数的取值范围;

(2)根据参数取值范围,通过蒙特卡洛随机采样方法生成n组参数,记为Xn(本研究中n=10 000);

(3)确定模型的模拟效率评判准则,对所有的参数组Xn进行水文模拟;

2.3 目标函数

2.3.1 单一目标优化

对于流域自身的参数率定,采用单一目标优化的方法,选取纳什效率系数(NS)为评价指标,本研究中将编号为i的流域中参数组θ的目标函数表示为Oi(θ),其计算公式如下:

(2)

NS用于描述洪水过程模拟序列与实测序列之间偏差的大小,NS值越接近1,表明模拟序列偏离实测序列的程度越小,模拟精度越高。

2.3.2 多目标优化

对于多个流域同步率定,需要综合考虑所有参与率定流域的模拟效果。本文通过构建基于折衷规划算法的多目标优化函数来实现多个流域水文模型的同时率定,获取能够适用于所有参与率定流域的参数组。由于每个流域都有自身目标函数的理想解(即单一目标优化时的最优解),而目标函数之间存在一定的冲突,所以考虑折衷规划的思路,寻求与所有理想解距离最小的折衷解。由于本研究中的目标函数具有相同的量度,于是简化为如下目标函数:

(3)

3 结果与分析

3.1 单个流域参数率定与移植

基于上述ROPE参数率定优选算法,选择单一目标函数,采用HBV模型对7个研究流域分别进行参数率定,并将单个流域的率定参数移用到其他6个流域,从而分析单个流域参数率定和移植的效果。

通过ROPE方法迭代,每个流域均得到了10 000组模拟结果非常接近的最佳参数组,由于这10 000组参数的模拟结果非常接近,为简单起见,

采用10 000组参数的模拟均值作为每个流域的模拟结果。图2展示了HBV模型在不同流域率定期与验证期的10 000组参数NS效率系数的平均值。可以看出模拟结果NS平均值均高于0.8,说明基于ROPE方法进行参数优选的HBV模型能很好地用于屯溪嵌套流域的洪水模拟,精度较高。

图2 单个流域率定参数的模拟NS均值Fig.2 Model performance for single catchment calibration and validation

图3显示了单个流域率定的10 000组参数移植到其他6个流域的模拟结果平均值。由图3可以看出,单个流域率定参数结果均优于移用其他流域参数的模拟结果。此外,不同流域的参数移植效果不同,大部分流域通过移植其他流域参数,能够获取理想的模拟结果(WA,CC,TX,LK,YT)。同时,部分移植参数的模拟结果较差(XT,YC),表明这些流域率定的参数不适合移植到其他流域。

图3 移植参数的模拟效果NS平均值Fig.3 Model performance by transferring parameters from other catchments

本文通过绘制流域参数率定和移植的结果矩阵图(图4),以进一步分析流域之间的参数移植效果。从矩阵图中可以明显看出,来自不同流域的移植参数的模拟效果值NS出现较大的变化。同时可以看出,该参数移植效果矩阵是不对称的,即存在两个流域A和B,当A流域率定的参数用于B流域时,可以取得理想的模拟结果,而B流域率定的参数却不适用于A流域的模拟。以YT流域和YC流域的NS指标为例,使用YC流域率定得到的参数在YT流域模拟,模拟效果NS平均值比YT流域自身率定NS值低0.2,移植效果差;而使用YT流域率定得到的参数在YC流域模拟,模拟效果NS平均值比YC流域自身率定NS值低0.07,移植效果合理。这种不对称现象,进一步印证了“过度拟合”现象的存在。即率定期内某些主要的径流过程可能在一个流域内出现频繁,参数在率定过程中体现了该过程的影响,使得参数的可移植性增强;而在另一个流域则可能历史洪水信息不够完整,参数对气象条件“过度拟合”,从而导致参数从一个流域到另一个流域的移植性较弱。另外考虑到由于屯溪流域及其嵌套子流域间的坡度值存在较大差异,对于汇流过程影响较大,造成流域的调蓄作用强弱不一,从而参数移植效果良莠不一。

图4 参数移植的NS平均值矩阵Fig.4 Color code matrix for parameter transfer

3.2 多流域同步率定

获取多个流域能够共同适用的模型参数,是降低参数不确定性、提高参数的可移植性的关键。本文通过折衷规划方法构建多个流域同时率定的多目标优化函数,并利用ROPE参数优选方法,对所有7个流域的模型参数进行同步率定。通过ROPE方法迭代,优选出了在7个流域均具有较好模拟效果的10 000 组参数,这10 000 组参数在7个流域模拟结果的NS最大值与最小值的波动都很小。7个流域通过同时率定获得的NS效率系数平均值分别为0.89、0.87、0.87、0.82、0.78、0.87、0.83。图5显示了7个研究流域分别在率定期与验证期进行单个流域率定的NS平均值,和7个流域的同步率定的参数模拟NS平均值,以及单个流域在率定期之间参数移植的模拟结果。由图5可以看出,对于率定期来说,同步率定的参数模拟效果略差于流域自身率定的效果,但二者差距并不大,在部分流域(LK)二者的模拟效果相当。而将同步率定的参数模拟效果,与单个流域率定参数的移植效果对比,可以明显发现,多个流域同步率定获取的参数组模拟效果良好,优于大多数单个流域率定参数的移植效果。验证期的模拟结果验证了这一结论,在绝大多数流域移植参数模拟结果确定性系数与自身接近,只在YC流域下降略显明显,分析其原因可能是YC流域面积较其他流域明显为小,流域的调蓄作用较弱,造成移植参数模拟效果不好。综合以上说明同步率定获得的参数组较直接移植更为稳定,在一定程度上能够降低“过度拟合”的影响。

图5 同步率定的模拟效果NS平均值Fig.5 Model performance for simultaneous calibration and validation

3.3 基于留一法的参数移植

在上一节的同步率定中,所有流域都参与了模型率定,在本节中,采用留一法以进一步验证同时率定情况下参数移植效果。即对于所选的7个流域,每次保留一个流域假定为无资料流域,利用剩余6个流域进行同时率定,并将同时率定获得的参数移植到该无资料流域。图6为其在率定期与验证期的模拟结果NS平均值。

图6 留一法同步率定的模拟效果NS平均值Fig.6 Model performance for leave one out calibration and validation

从图6中可以看出,对于率定期来说,与上一节7个流域的同步率定的模拟结果NS平均值相比,留一法多流域同步率定的模拟结果NS平均值因未使用目标流域的数据,模拟结果略差于采用所有流域进行同步率定的结果,但与单个流域之间的参数移植结果NS平均值相比,留一法多流域同步率定的参数移植到TX、WA、CC和LK流域的结果具有明显优势。同时注意到,留一法多流域同步率定的模拟结果NS平均值与流域自身率定的NS平均值相差不大,甚至在部分流域(LK)二者模拟效果相当。验证期的模拟结果同样满足该规律,注意到与上节同步率定相比,留一法移植后的参数模拟结果NS值下降非常微小,移植参数具有较高的确定性。该结果表明,在无资料地区洪水模拟中,可用多流域同步率定法进行参数移植,以便获得较为稳定的移植效果,该方法在一定程度上降低了模型参数的不确定性,这将有助于提高无资料地区洪水预报的精度。

4 结 论

(1)采用基于深度函数的ROPE参数优选方法,HBV模型能很好地用于屯溪嵌套流域的洪水过程模拟,可以取得较好的模拟精度。该参数优选方法能够获取具有良好稳定性的模型参数,是水文预报参数选择的一种新途径。

(2)对于屯溪及其嵌套的6个子流域,单个流域之间的HBV模型参数移植结果参差不一,应考虑到观测资料对于参数不确定性的影响,不宜直接进行参数移植。

(3)多流域同步率定能有效降低HBV模型参数不确定性的影响,并获得较优的移植效果,提高无资料地区参数移植稳定性和洪水预报的参数确定性。水文模型不确定性来源众多,而本文仅讨论了针对模型参数的不确定性,对于自然不确定性、数据不确定性以及模型结构的不确定性还需加以研究,多流域同步率定能否普遍适用于水文模型的参数移植尚需进一步验证。同时,在多流域同步率定方法中还应进一步加强对水文相似性、流域个数等因素的考虑。本文所进行的是逐日资料的模拟,因此在模拟中主要关注了长期的水量平衡,所以未对峰现时间和洪峰模拟结果进行比对。在中小河流洪水预报中,小时尺度洪水的模拟更为重要,同步率定方法在次洪模型中的应用及洪水主要特征值的分析是我们下一步研究的重点。

猜你喜欢
水文平均值流域
发展水文经济 增强水文活力
昌江流域9次致洪大暴雨的空间分布与天气系统分析
浅谈水文档案的价值和开发利用
巧用1mol物质作标准 快速确定混合物组成
河南省小流域综合治理调查
江西省水文文化建设的思考
称“子流域”,还是称“亚流域”?
变力做功时运用F=F1+F2/2的条件
平面图形中构造调和平均值几例
2007/2008年度桑蚕干茧质量分析报告