水文模型不确定性研究方法耦合分析

郑言峰

(青岛市水利勘测设计研究院有限公司，山东 济南 250013)

水文模型的不确定性评价已经成为水文学者研究的热点问题之一[1-4]。水文模型不确定性的来源是多方面的，总体上可以分为以下3类：实测资料(降雨、蒸发、径流等)的不确定性、模型结构的不确定性以及模型参数的不确定性。为了综合分析以上因素对模型不确定性的影响，Beven和 Binley[5]于1992年提出了通用似然不确定性估计方法，即GLUE方法。该方法将Monte Carlo随机取样技术与Bayesian理论结合起来，原理简单，易于操作。然而，GLUE方法在很多实际应用中的结果[6]表明预测区间的覆盖率远低于给定的置信水平。XIONG等[7]通过引入“模拟偏差曲线”对GLUE方法进行了改进，在一定程度上提高了区间覆盖率，但同时也增加了区间带宽。

数据同化方法将模型和观测这两种研究手段进行了结合[8]，从而实现了在观测资料更新的同时可以更新模型的计算结果。它的主要思想是将最新的观测数据融合到所使用的模型之中，通过不断更新的观测数据来减小状态变量模拟误差，继而使用更新之后的状态变量重新预报从而来达到提高预报精度的目的。作为数据同化方法的典型代表，集合卡尔曼滤波(Ensemble Kalman Filter，EnKF)已被广泛应用到海洋、陆面和大气数据同化领域之中[9]。EnKF方法融合了蒙特卡洛(Monte Carlo，MC)采样技术，可以分别考虑模型观测资料的不确定性以及模型结构的不确定性。2005年，Moradkhani等[10]在EnKF的基础上又加入了模型参数的更新(即dual EnKF)，用以描述模型参数的不确定性。

为了探索解决在提高区间覆盖率的同时降低预测区间平均带宽问题，本文将dual EnKF方法与GLUE方法2种水文模型不确定性研究方法进行了耦合，得到了dual E-G方法，并以新安江模型为例，研究了dual E-G方法在汉江上游3个子流域(即甲河流域、天河流域和灞河流域)中的适用性。该成果可为水文模型预报、洪水决策风险分析等方面的深入研究提供有益参考。

1 研究区概况

汉江又称汉水，是长江中下游最大的一条支流，发源于秦岭南麓，干流流经陕西和湖北两省，于武汉市汇入长江。汉江干流全长1 570多km，流域面积大约15.9万km2。干流丹江口以上为上游，流域面积大约9.5万km2。汉江流域河流密布，沟壑交织，分布着大小河流及山沟数以千计。汉江流域幅员广阔，光、热和水资源空间差异比较大，气候特征多变，其降水主要来源于东南、西南2股暖湿气流，因此降水变率比较大，是中国旱、涝灾害多发的地区之一。作为南水北调中线工程水源地的丹江口水库，不仅承担着湖北武汉、襄阳等23个县市1亿多人口及124万多hm2耕地的防洪任务，还承担着河南、河北、北京、天津等四省市的20多座城市的供水任务，意义重大。但是21世纪以来，汉江流域出现的干旱问题，严重影响了流域内的生态环境与供水工程。研究汉江流域上降雨径流模拟的不确定性，不仅可以为决策者提供充足的风险信息，还可以实现水资源的合理配置，从而保证南水北调中线工程的顺利实施。

因此，选择了汉江上游的3个子流域(即甲河流域、天河流域和灞河流域)作为研究对象。3个子流域均为湿润地区，且资料长度均为8年。前5年(1980—1984年，共1 827天)的流域实测平均日降雨、日蒸发和日径流作为率定期，后3年(1985—1987年，共1 095天)的流域实测平均日降雨、日蒸发和日径流作为检验期。流域水文资料的具体情况见表1，3个子流域在汉江中的位置及流域见图1。

表1 研究流域资料

图1 甲河、天河及灞河在汉江中的位置及其流域

2 模型简介及评价指标

2.1 新安江模型简介

对于研究水文模型的不确定性而言，概念性水文模型优于系统模型和物理模型。概念性模型主要利用一些简单的物理概念和经验公式，来近似描述水流在流域的运动状态。系统模型主要是将所研究的流域或区间看作一种动力系统，利用输入与输出资料，建立某种数学关系，这种模型只关心模拟结果的精度，并不考虑输入与输出之间的物理关系。物理模型主要根据水流连续方程和动量方程来求解水流在流域的时间和空间变化规律，虽然模型物理意义更加明确，但其对资料要求过于严格，一般流域很难满足。

新安江(三水源)模型是河海大学(原华东水利学院)赵人俊等在1973年提出来的降雨径流模型，简称新安江模型。由于新安江模型采用蓄满产流进行产流计算，因此主要适用于湿润、半湿润地区。近几十年，新安江模型不断改进，已成为中国应用较为广泛的一个流域水文模型。

新安江模型作为国内常用的一种概念性水文模型，大部分模型参数都具有明显的物理意义，不仅模拟精度高，而且对水文资料的要求比较低，对于日模型来说，只需要流域日蒸发资料、日降雨资料以及日径流资料即可。新安江模型共有15个模型参数。①K：流域蒸散发能力Ep与蒸发皿实测值EM之比，简称蒸散发折算系数；②IMP：一般将整个流域面积分为透水面积与不透水面积，IMP是指全流域面积中不透水面积所占的比例；③B：蓄水容量曲线的指数；④WUM：流域上层蓄水容量(mm)；⑤WLM：流域下层蓄水容量(mm)；⑥WDM：流域深层蓄水容量(mm)；⑦C：深层蒸散发系数；⑧SM：自由水蓄水容量(mm)，反映表层土蓄水能力的大小，直接决定着地面径流和地下径流所占的比重；⑨EX：自由水蓄水容量曲线的指数；⑩KG：地下水出流系数；KI：壤中流出流系数，主要取决于表层土的渗透性；CG：地下径流的消退系数，反映地下径流的消退规律；CI：壤中流的消退系数，反映壤中流的消退规律；nn：反映流域调蓄能力的参数，相当于线性水库的个数或水库的调节次数；KK：线性水库的蓄泄系数，相当于流域汇流时间的参数，具有时间因次。

2.2 似然函数的选择

在对水文模型进行率定时，需要定义一个评价模拟实测流量好坏的标准，即似然函数。似然函数必须满足一定的特点：当模拟系列与实测系列没有相似性时它的值应为0，并且相似性越高，它的值应越大。不同的似然函数率定得到的参数也会有所不同，有的似然函数对高水的模拟效果较好，有的似然函数对低水的模拟效果较好，而纳西效率系数因对整体模拟效果较好而被广泛应用于水文模型参数估计研究中。因此，采用纳西效率系数R2作为似然函数，表达见式(1)：

(1)

当R2的值等于1时，模拟流量序列完全重现了相应的实测流量序列，此时的模拟效果最好；当R2的值等于0时，模拟流量序列对同期观测的拟合程度与实测流量序列长时间尺度上的平均值对实测流量序列的拟合程度相当；R2的值越小，模型模拟的序列对相应观测的拟合程度越差[11]。

2.3 预测区间评价指标

为了反映水文模型模拟和预报实测流量的不确定性，每个实测流量的预测结果不应该只是单纯的点估计，而应该是给定置信水平下的区间估计。XIONG等[12]从区间覆盖率、区间带宽和区间不对称性3个方面详细介绍了7个用于评价预测区间好坏的指标。因此，从3个方面中各选择一个比较常用的指标(区间覆盖率、平均带宽和平均不对称程度)来评价预测区间的性质。它们的定义如下。

a)区间覆盖率(CR)。区间覆盖率是指实测流量落入预报区间的个数与实测流量数据总个数的比值。自从Beven提出GLUE方法以后，区间覆盖率就被广泛应用于评价预测区间的好坏。CR越大，表示预测区间覆盖率越高，CR最为理想的值为1，这意味着预测区间覆盖了所有的实测流量数据。

b)平均带宽(B)。平均带宽是所有流量的预测区间上下界差值的平均值，单位：m3/s。

(2)

QU(j)、QL(j)分别表示第j时刻的流量预测区间上下界。对于指定的置信水平，在保证有较高的覆盖率前提下，预报区间的平均带宽越窄越好，此时的预测区间的不确定性也较小。

c)平均不对称程度(S)。

(3)

理想的预测区间是不仅能覆盖大部分的实测流量，而且预测区间的平均带宽值较小，并且预测下界与实测流量之间的距离近似等于预测上界与实测流量之间的距离。

3 不确定性研究方法

3.1 GLUE方法

GLUE方法中的一个很重要的观点是：导致模型模拟结果好坏的不是模型单个参数，而是模型的参数组合。目前GLUE方法已成为水文模型不确定研究的主流方法[1-3,13-14]。GLUE方法的主要思路如下：选定模型及似然函数，并设定似然阈值；在预先设定的参数分布(一般假定为均匀分布)取值空间内，利用Monte Carlo随机取样获得一定数量的有效参数组(即似然函数值高于似然阈值的参数组)；在给定的置信水平下，求出每个实测流量的预测区间。

GLUE方法共有3个参数，即似然函数阈值TH、有效参数组的个数M、置信水平α。每个参数的取值如下：TH=0.5，M=2000，α=0.8。

3.2 dual EnKF方法

在EnKF方法中假定水文模型的参数不随时间变化，模型的参数可以通过模型率定获得。但是由于模型驱动数据和观测数据存在一定误差，不能保证模型的模拟效果不会随着时间的变化而有所差异。因此，有必要将模型参数的更新纳入预报的过程之中[10,15-16]。模型参数的更新过程类似于状态变量的更新过程。

dual EnKF方法的具体步骤如下。

步骤一选择水文模型，准备必要的驱动数据(流域平均日降雨资料P、流域平均日蒸发资料E)和流域实测日流量资料Q，确定一组有效参数组θ0。

步骤二确定集合成员个数N、日蒸发的测量误差超参数γE、日降雨的测量误差超参数γP和日径流的测量误差超参数γQ。

步骤三分别对日蒸发、日降雨和日径流资料加上一系列高斯白噪声扰动。

(4)

(5)

(6)

步骤四计算每个成员的状态变量先验值。

(7)

t=1时，

步骤五计算每个成员的模拟流量。

(8)

步骤六计算每个模型参数的增益。

(9)

步骤七计算每个成员t+1时刻的模型参数更新值。

(10)

步骤八使用更新之后的模型参数重新计算状态变量值。

(11)

步骤九使用更新之后的模型参数重新计算模拟流量值。

(12)

步骤十计算每个状态变量的增益。

(13)

步骤十一计算每个成员的状态变量后验值。

(14)

dual EnKF方法共有4个参数，即集合成员个数N和3个超参数γE、γP、γQ。随着集合成员个数N的增加，集合成员的代表性逐渐增强、分析误差逐渐减小，但却加重了计算负担。因此，需要选择一个合适的集合数既能满足集合成员代表性的要求，又要考虑计算负担问题[17]。超参数的大小直接影响着驱动数据和观测数据误差的大小。Moradkhani等[10]通过研究发现当N=50、γE=γP=γQ=0.1时，模拟效果较好。因此，本文中的4个参数的取值与Moradkhani的研究结果保持一致。

3.3 dual E-G方法

目前大多数基于GLUE方法的水文模型不确定性研究结果表明：随着覆盖率的增高，平均带宽反而增大[7,12,14]，无法提供有效的决策依据。

为了探索解决在提高区间覆盖率的同时降低预测区间平均带宽问题，将dual EnKF方法与GLUE方法2种水文模型不确定性研究方法进行了耦合，得到了dual E-G方法。具体耦合方法如下：①首先按照GLUE方法产生似然函数大于阈值TH的M组有效参数组；②将每组模型参数产生的服从均值为参数本身，方差为参数的0.1倍的N个参数组作为每个成员的参数初始值；③然后按照dual EnKF方法进行模型参数更新；④最后将得到的N个模拟流量系列的均值作为该组参数的最优模拟系列，则M组有效参数组便得到M个模拟流量系列；⑤将M个模拟流量系列同一时刻的所有模拟流量Q(i,j)分别从小到大排序；⑥对于给定的置信水平α(0<α<1)，分别计算第j时刻实测流量的下界预测值QL(j)和上界预测值QU(j)。

4 实例计算与分析

为了分析耦合效果，以新安江模型为例，模型参数取值范围见表2。对于给定的置信水平α(0<α<1)，采用3个评价指标(覆盖率、平均带宽和平均不对称性)评价了dual E-G方法在3个流域(甲河、天河、灞河)的预测区间性质，并与GLUE方法进行了对比，见表3、4。由表4可知，相比于GLUE方法，dual E-G方法在甲河、天河和灞河流域中计算得到的率定期的区间覆盖率分别提高了28%、41%和29%，平均带宽分别降低了8%、18%和6%，平均不对称程度分别降低了62%、75%和44%；检验期的区间覆盖率分别提高了31%、53%和32%，平均带宽分别降低了4%、1%和2%，平均不对称程度分别降低了67%、80%和59%。综上所述，dual E-G方法可以同时改善3种不确定性指标。

表2 新安江模型参数取值范围

表3 GLUE方法和dual E-G方法在不同流域上的预测区间性质

表4 dual E-G方法相比于GLUE方法的改善程度 %

为了便于比较不同方法计算出来的预测区间的性质，本文选用了甲河流域率定期中的1984年和检验期中的1985年的预测区间图来比较GLUE方法和dual E-G方法预测区间的优良性，见图2。可以看出，与GLUE方法相比，EnKF-G方法率定期模拟洪峰能力提高，检验期的低水部分模拟效果更好一些，且预测区间改善比较明显，不仅提高了区间覆盖率、降低了区间平均不对称程度，还在一定程度上降低了平均带宽。

a)GLUE方法率定期

d)dual E-G方法检验期

为了说明dual E-G方法降低不确定性的原理，以GLUE方法筛选出的新安江模型的第1组有效参数为例进行分析。对于GLUE方法，直接采用第1组有效参数计算的模拟流量系列计算纳西效率系数；对于dual E-G方法，将由第1组有效参数生成的服从均值为模型参数本身，方差为模型参数的0.1倍的50个参数组作为每个成员的初始参数组，然后再按照dual EnKF方法进行参数向量更新，从而计算出50个成员的模拟流量系列，将集合均值作为各时刻流量的最优预报值，组成最优模拟系列，进而计算纳西效率系数，结果见表5。由表5可知，dual E-G方法模拟精度明显高于GLUE方法，甲河流域率定期纳西效率系数提高了0.26，达到了0.89；灞河流域检验期纳西效率系数提高了一倍以上，达到0.62。

表5 第1组有效参数模拟系列纳西效率系数对比

GLUE方法能够充分考虑模型参数的空间不确定性，筛选出大于似然阈值的有效参数组，但是由于未能对参数进行实时更新，随着时间的变化，模型参数的时间不确定性逐渐累积，导致模拟精度有所降低；dual EnKF方法能够充分考虑模型参数的时间不确定性，对于给定位置的有效参数组能够通过实时数据及误差扰动进行不断更新，以达到最优的模拟效果，但是无法体现下垫面变化导致的参数空间不确定性影响。dual E-G方法将GLUE方法考虑的模型参数空间不确定性与dual EnKF方法考虑的模型参数时间不确定性进行了充分结合，使得模拟效果更加接近实际，从而表现出更高的模拟精度和更优的预测区间。

5 结语

单纯采用GLUE方法进行水文模型不确定性分析时不可避免地会产生以增加平均带宽为代价来提高区间覆盖率或以降低区间覆盖率为代价来减小平均带宽的问题，主要原因是由于“异参同效现象”的存在导致最优参数很难获得。dual EnKF方法可以根据每个成员的状态变量先验值，通过迭代更新使得每组参数都能够尽可能达到最优，从而提高每个成员的模拟效果。耦合后的dual E-G方法，能够充分考虑模型参数的空间不确定性和时间不确定性，可以使得M组有效参数组都能够尽可能达到最优，不仅可以提高模型预报精度，还可以改善预测区间性质，即可以同时改善区间覆盖率、平均带宽和平均不对称程度。

以新安江模型为例，研究了dual E-G方法在汉江3个子流域上的适用性。该成果可为水文模型预报、洪水决策风险分析等方面的深入研究提供有益参考。然而水文模型的不确定性是一个非常复杂的系统问题，dual E-G方法在不同模型、不同流域上的适用性仍有待进一步研究。