洪水概率预报评价指标研究*

2020-03-27蒋晓蕾梁忠民胡义明李彬权

湖泊科学 2020年2期

蒋晓蕾，梁忠民，胡义明，王军，李彬权

(河海大学水文水资源学院,南京 210098)

由于自然过程的复杂性和人类认识水平的局限性，使得单一预报结果不可避免的存在不确定性，因此发展起来了集合预报与概率预报[1]. 集合预报给出了预报量可能的结果集，是采用集合形式量化预报的不确定性. 而概率预报，顾名思义，估计了预报量(如流量)的概率分布，通过概率形式量化预报的不确定性. 一方面，集合预报仅提供了预报量有限的可能结果，而概率预报量化了预报量所有结果的可能性，因此，概率预报是集合预报的一种“升华”，较集合预报提供更为全面的信息. 另一方面，由于很难直接估计预报量的理论分布函数，因此，往往采用预报量有限的可能结果集来近似替代其分布函数，估计其预报不确定性，即集合预报可以视为概率预报的一种表现形式[2]. 本文将集合预报视作概率预报的一种表现形式，对洪水概率预报结果的评价体系开展研究.

就方法而言，目前的洪水概率预报方法可以大致分为两类途径：全要素耦合途径和总误差分析途径.

全要素耦合途径是指分别识别、量化洪水预报过程中各阶段的不确定性，并将其与确定性预报模型耦合，从而实现洪水概率预报. 这些不确定性一般包括输入(如降雨)[3-4]、模型结构[5]、模型参数和模型状态[6]的不确定性. 代表性方法包括贝叶斯总误差分析方法BATEA[7](BAyesian Total Error Analysis)、贝叶斯综合不确定性估计方法IBUNE[8](Integrated Bayesian UNcertainty Estimator)等. 这类方法虽然可以量化洪水预报过程中的各类不确定性要素，识别不确定性来源，但往往计算耗时，不利于生产实际应用. 总误差分析途径是指不识别洪水预报具体过程的不确定性，仅对最终预报结果进行分析，量化其不确定性，进而实现概率预报[9-10]. 常用的方法包括水文不确定性处理器HUP[11](Hydrologic Uncertainty Processor)、模型条件处理器MCP[12](Model Conditional Processor)、三维误差矩阵[13]等. 这一途径避免了处理洪水预报子过程的不确定性，往往计算省时，结构灵活，可以与各类确定性预报模型进行衔接，但是无法溯源预报的不确定性.

1 概率预报“精度-可靠度”联合评价指标体系

本节首先对现有的精度评价指标和可靠度评价指标进行分析和总结，在可靠度评价指标中，提出覆盖率判定系数用以评价概率预报的整体可靠度. 在此基础上，建立了洪水概率预报“精度-可靠度”联合评价指标体系.

1.1 精度评价指标

(1)洪峰相对误差REP(%)和洪量相对误差REV(%).

《水文情报预报规范》[20]规定：降雨径流预报以实测洪峰流量的20%作为许可误差. 因此，REP的许可误差为±20%.

将洪量平铺于流域面积，可以计算出相应的径流深，因此，径流深误差与洪量误差相等. 《水文情报预报规范》[20]规定：径流深预报以实测值的20%作为许可误差. 因此，可以认为REV的许可误差为±20%.

REP和REV均评价了预报结果的偏差(Bias，预报是否偏大/偏小于观测).

(2)洪峰滞时ET.

洪峰滞时指预测洪峰滞后于实测洪峰的时段长. 《水文情报预报规范》[20]规定：滞时绝对值的许可误差为预见期的30%，许可误差小于3 h或一个计算时段长时，则以3 h或一个计算时段长作为许可误差.

(3)确定性系数/纳什系数[21](Deterministic Coefficient/Nash-Sutcliffe Efficiency，NSE).

(1)

NSE评价了洪水过程的拟合效果，其取值范围是(-∞,1]，NSE越接近于1，表明洪水过程的拟合效果越好. 除了NSE以外，均方误差(Mean Squared Error，MSE)[22]也可以用来评估洪水过程的拟合情况：

(2)

理论上，NSE与MSE存在转化关系[22]：

(3)

由公式(2)和(3)可知，MSE的取值受到了流量级别的影响，不利于对比不同量级洪水过程的拟合效果；而NSE在一定程度上消除了流量量级的影响，因此，在实际应用中NSE比MSE应用更为广泛，更具说服力. 对NSE进行分解，可以将NSE分解为3个评价因子[22]：

NSE=2·α·r-α2-β2

(4)

(5)

(6)

(7)

式中，σpre为倾向值预报的标准差；σobs为实测流量的标准差；μpre为倾向值预报的均值；μobs为实测流量的均值；Covpre,obs为倾向值预报与实测流量的协方差；α为方差因子；β为均值因子；r为预报与实测的线性相关系数.

由公式(4)～(7)可知，NSE评估了倾向值预报与实测流量之间均值和标准差的差异，同时考虑了两者的线性相关程度. 当倾向值可以完美拟合实测流量过程时(完美预报)，两者的均值相等、标准差相等、线性相关程度最高，即α=1；β=0；r=1，此时NSE=1. 根据NSE的评价因子(β、α、r)，Gupta等[22]提出了Kling-Gupta Efficiency (KEG)评价指标：

(8)

(9)

(10)

G3=(r-1)2

(11)

式中, 变量含义与式(5)～(7)中变量相同.KGE的取值范围为(-∞,1]，KGE越接近于1，表明洪水过程的拟合效果越好.

类似于NSE，KGE也考虑了均值因子G2、方差因子G1和线性相关性因子G3共3个评价因子. 不同于NSE，KGE中因子对KGE的影响关系更直观：KGE与3个因子均为单调递减关系. 因此，通过计算KGE及其评价因子，可以提供预报流量更细致的精度评价.

(4)基准系数[23](Benchmark Efficiency，BE).

(12)

式中，Qb(i)为基准预报流量，m3/s；其他变量含义同前.

基准系数BE评价了模型的预报技巧(Forecast Skill，预报相对于某一标准是否具有较高的精度)，一般用于对比倾向值预报Qpre与基准预报Qb在洪水过程拟合中的相对好坏：BE=0说明Qpre与Qb在洪水过程拟合中表现相当；BE>0说明Qpre比Qb更优；BE<0说明Qpre比Qb拟合效果差.BE在形式上与NSE相似，可以认为NSE是以实测平均值作为基准预报的BE值.

1.2 可靠度评价指标

(1)区间离散度[24-26]：

(13)

式中，DI指某一置信度条件下预报区间的平均离散度；N为预报时刻数；D(i)指第i时刻预报区间的离散度(预报区间宽度与实测的比值)，具体计算如下：

(14)

式中，qu(i)、qd(i)指预报区间的上限和下限，m3/s；其他变量含义同前.

D(i)的计算过程消除了流量级别对区间宽度(b(i)=qu(i)-qd(i))的影响，因此DI的取值不受流量级别影响，可用以对比不同量级洪水的区间预报结果：离散度越小区间预报效果越好，但是有关DI阈值的研究仍较为匮乏. 我国《水文情报预报规范》[20]中规定洪峰预报允许误差为20%，借鉴这一概念，本文将洪峰处的区间离散度定义为洪峰离散度Dpeak，并确定其阈值如下：

将洪峰预报区间边界与实测洪峰相对距离的最大允许值定为0.2，即:

(15)

(16)

式中，qu，peak和qd,peak分别为洪峰预报区间的上下边界，m3/s；Qobs,peak为洪峰观测流量，m3/s.

(2)区间覆盖率[24-25](Containing Ratio,CR)：

(17)

(18)

理论上来说，CR值越接近区间置信度，概率预报结果越合理. 如置信度为90%的区间预报结果，其CR值越接近90%，该预报区间越合理.

(3)覆盖率判定系数(Containing Ratio Coefficient,CRC).

区间覆盖率CR只能考虑某一个置信度(如90%)预报区间包含实测的情况，并不能代表整个概率预报结果的合理性，因此，本文提出覆盖率判定系数CRC评价指标，对所有置信度区间预报结果做出整体评估.

根据区间覆盖率的含义，当概率预报为完美预报时，对所有的置信度X=0,…,100%，其对应的覆盖率CR，应满足CR=X. 点绘(X,CR)时，所有的点据应位于45°线上. 然而，实际中不可能做到完美预报，(X,CR)点往往散落在45°线附近，此时，可以计算这些点据与45°线的确定性系数R(即CRC)，来判断概率预报的合理性，覆盖率判定系数计算公式如下：

(19)

CRC的取值范围是(-∞，1]，CRC=1为完美概率预报. 当CRC大于某一阈值时，线性关系CR=X显著，概率预报结果合理，而这一阈值随实际研究问题不同而变化，应通过统计检验来确定. 在水文实际应用中，一般认为相关系数r≥0.8时变量存在较显著线性相关性，根据r与确定性系数R的关系(r2=R)，本文认为在不便于采用统计检验方式确定CRC阈值时，可以将0.64作为CRC的阈值，认为CRC≥0.64时45°线可以较好的拟合(X,CR)点据，即概率预报结果合理.

由CRC计算公式可知，CRC同时考虑了所有置信度条件下的区间预报结果，可以衡量概率预报的整体合理性.

就思路而言，CRC与Q-Q(Quantiles-Quantiles)图方法[27]、概率图(Probability Plot)法[28]类似，三者的45°线均代表了完美预报. 然而，Q-Q图方法评价的是各个分位数的合理性，概率图法评价了累计概率的合理性，CRC考虑的是区间覆盖率的合理性，三者的关注点不同. 考虑到实际洪水概率预报中多关注区间预报的覆盖率，因此本文认为CRC方法较Q-Q图方法、概率图方法更适用于实际生产应用.

(4)单位区间离散化系数.

一般而言，区间宽度越大，离散度越大，区间覆盖率越大，反之亦然，区间离散度和覆盖率在某种程度上是一组矛盾的评价指标. 为了解决这一矛盾，提出了单位区间离散化系数[17, 29](Percentage of observations bracketed by the Unit Confidence Interval，PUCI)，用以描述“平均相对区间宽度”所包含的实测点据比例.

(20)

式中变量含义同前. 理论上，PUCI的取值范围是[0,∞).PUCI值越大，表明该置信度的区间预报结果合理性越差.

对不同置信度预报区间的PUCI值进行平均，可以得到平均单位区间离散化系数(ACI)：

(21)

式中变量含义同前. 通过计算ACI值，可以估计区间预报的整体合理性：ACI越小，概率预报结果越合理. 然而，由于ACI的允许误差很难界定，因此不建议采用ACI对概率预报结果合理性进行评估.

(5)连续排位评分(Continuous Ranked Probability Score,CRPS)和连续排位技巧评分(Continuous Ranked Probability Skill Score,CRPSS)

CRPS[30-31]评估了概率预报的预报能力. 对任一时刻，CRPS计算方式如下：

(22)

CRPS值[31-32]越小，说明预测分布在观测值附近越呈现“尖瘦”形态，预测分布的“准确性”越高. 但目前研究尚未给出CRPS的合理阈值，往往采用CRPS比较多个预报的相对预报能力[33]，进而提出了CRPSS[19]，其计算公式如下：

(23)

式中，CRPSfore为待评价概率预报的CRPS值；CRPSref为参考预报的CRPS值，对于确定性预报，其CRPS值即为平均绝对误差[30].

连续排位技巧评分CRPSS是一种技巧(Skill)评分指标，取值范围是(-∞，1]. 当CRPSS>0时，概率预报较参考预报具有更高的预报能力；当CRPSS<0时，待评价概率预报的“准确性”不如参考预报；当CRPSS=0时，待评价概率预报与参考预报预报能力相当.

(6)区间对称度和对称比

对称度[24]Λ1评估了某一置信度预报区间关于实测流量的对称性，其计算公式如下：

(24)

λ1(i)=|h(i)-0.5|

(25)

(26)

式中变量含义同前. 由Λ1计算公式可知，当i时刻的预报区间关于实测流量对称时，h(i)=0.5，λ1(i)=0.Λ1值越小，预报区间关于实测对称性越强. 一般的，要求Λ1<0.5. 此外，还可以采用平均对称度Λ2对区间对称性进行评价：

(27)

(28)

式中变量含义同前.Λ2值越小，预报区间关于实测对称性越强，当Λ2=0时，预报区间关于实测流量完全对称. 一般的，要求Λ2<1.

区间对称比[15]Λ3采用实测值大于区间上界的个数与小于区间下界个数的比值来衡量预报区间关于实测流量的对称性：

(29)

(30)

(31)

由于区间预报结果往往并不要求对称，同时，区间覆盖率与离散度已经在某种程度上反映了预报区间相对实测流量的对称性. 因此，对称性评价指标(对称度和对称比)在实际生产中的应用较少.

1.3 建议的“精度-可靠度”联合评价指标

由于洪水概率预报结果一般可以分为倾向值预报和区间预报两个部分，因此，应采用“精度-可靠度”联合评价方式对洪水概率预报结果进行评估，即分别评价倾向值预报的准确性(精度评价)和区间预报的合理性(可靠度评价)，具体评价指标如下：

在精度评价时，根据前文的介绍与分析，参考《水文情报预报规范》[20]的相关规定，建议采用：纳什系数NSE评价洪水过程，洪峰误差REF、洪峰滞时ET评价洪峰拟合情况，洪量误差REV评估水量平衡情况，基准系数BE比较不同定值预报(倾向值/确定性)对洪水过程的拟合情况. 当需要对定值预报过程进行更为细致的评价时，可以采用KGE及其3个评价因子G1、G2、G3对洪水过程的拟合效果进行评价.

在“可靠度”评价时，常关注置信度为90%的区间预报结果，因此，建议采用：洪峰处区间离散度Dpeak、单位区间离散化系数PUCI和覆盖率CR对该区间预报结果进行评价. 由于90%预报区间仅是众多预报区间中的一个，并不能代表其他置信度的预报区间结果，因此，应采用覆盖率判定系数CRC给出概率预报结果合理性的整体评估. 此外，可以采用连续排位技巧评分CRPSS对比多个概率预报结果的相对预报能力.

相关评价指标的具体计算公式见前文，此处不再赘述.

2 洪水概率预报评价指标应用

以淮河王家坝区间流域为例，在新安江模型预报基础上，采用水文不确定性处理器(Hydrologic Uncertainty Processor，HUP)[11]进行洪水概率预报，并采用“精度-可靠度”联合评价指标体系对洪水概率预报结果进行评估.

2.1 研究区概况

王家坝区间流域集水面积为7110 km2. 区域内年平均降水量800～1200 mm，降水年际变化大，时空分布不均匀. 洪水主要来源为上断面(息县、潢川、班台)河道来水和区间降雨产生的径流，区间降雨的汇流时间约72 h(雨峰至洪峰出现时间). 研究区域示意图见图1.

图1 研究区域示意图Fig.1 The sketch of research area

2.2 确定性预报

新安江模型是河海大学赵人俊[34]教授设计的国内第一个完整的流域水文模型. 三水源新安江模型蒸散发计算采用三层模型；产流计算采用蓄满产流模型；采用自由水蓄水库结构将总径流划分为地表径流、壤中流和地下径流3种；流域汇流计算采用线性水库；河网汇流采用滞洪演算法；河道汇流采用马斯京根分段连续演算法. 三水源新安江模型共有17个参数，包括7个产流参数、10个汇流参数.

以王家坝区间流域洪水为研究对象，上断面河道来水采用马斯京根法进行汇流计算. 选用1990-2013年共24年的资料进行日模型率定，对其间的12场洪水资料进行次洪模型计算(其中8场洪水用于率定，4场洪水用于验证)，本文以计算时段为3 h为例开展相关研究，新安江模型参数见表1，确定性预报精度统计见表2.

由表2可知：王家坝区间流域新安江模型预报的洪峰误差在4%～31%，最小洪峰误差为4.07%，有两场洪水洪峰误差超过20%，合格率为83%. 12场洪水的洪量误差均在15%以内；洪峰滞时在3个时段(9 h)以内；确定性系数均大于0.75. 综合而言，新安江模型预报洪水精度尚可.

表1 王家坝区间流域新安江模型参数

表2 新安江模型精度统计

2.3 概率预报模型

HUP模型[11]由Krzysztofowicz 和 Herr(2001年)提出，已被广泛用于洪水概率预报研究中[35-38]，本文仅对该模型进行简单介绍，其他细节可参考文献[11,39-40].

HUP模型中，首先采用正态分位数转换技术[41](Normal Quantile Transform)将确定性预报(新安江模型预报)结果与实测结果转换至正态空间. 然后在正态空间进行线性似然函数假定，并根据贝叶斯理论推求预报量的后验分布函数. 最后通过正态分位数逆转换过程，将正态空间的预测分布函数转换至原始空间，进而实现概率预报. 预报量hn的后验概率密度函数可以表示为：

(32)

表3 HUP模型相关参数

式中,sn为n时刻的预报量；n代表预见期，本文预见期为3 h(与确定性预报一致)；h0为起报时刻实测流量值；An、Bn、Dn、Tn为模型参数，结果见表3；Ω(·)表示正态分位数转换，Ω-1(·)表示正态分位数的逆转换过程；Γ(·)表示新安江模型预报或实测流量的边际分布，一般采用三参数对数威布尔分布，对数威布尔分布函数形式如下：

(33)

式中,x为新安江模型预报流量或实测流量；a、b、c为对数威布尔分布的3个参数. 采用上述8场洪水率定HUP模型相关参数，参数率定结果见表3，边际分布拟合情况如图2所示.

图2 边际分布拟合图Fig.2 Fitting graph of marginal distribution

由图2可知，对数威布尔分布理论概率曲线可以很好地拟合流量的经验累积概率. 计算对数威布尔分布与经验概率的平均绝对值误差：h0、h1和s1的平均绝对值误差均不超过0.02，因此，概率预报模型的边际分布拟合效果较好.

2.4 概率预报评价

在精度评价中，采用洪峰误差REF、洪量误差REV、洪峰滞时ET对HUP模型的倾向值预报(Q50预报值和均值预报)进行评估，并采用基准系数BE对比倾向值预报与新安江模型预报结果，相关数据见表4.

由表4可知：倾向值预报的REF在6%以内，小于新安江模型洪峰误差；倾向值预报的REV在7%以内；洪峰滞时ET在1个时段(3 h)以内. 通过倾向值预报的REF、REV和ET比较结果可知，概率预报的Q50和均值预报均比新安江模型预报精度更高. 此外，倾向值预报的BF均大于0，说明Q50预报和均值预报的洪水过程拟合效果比新安江模型更好. 为了更细致的对比倾向值预报的精度，分别计算了新安江模型和倾向值预报(Q50和均值)的KGE及其3个评价因子取值，见表5.

表4 洪水概率预报精度评价

由表5可知：概率预报倾向值预报(Q50和均值)的KGE均比新安江模型大，说明Q50和均值对实际流量过程的拟合精度较新安江模型预报更高. 此外，新安江模型的3个评价因子取值较小，说明新安江模型预报在方差、均值和线性相关性方面的模拟精度较高，而倾向值预报(Q50和均值)的3个评价因子几乎都小于新安江预报，说明倾向值预报在均值、方差和线性相关性的模拟精度都大于新安江模型.

在可靠度评价中，采用洪峰处区间离散度Dpeak、覆盖率CR和单位区间离散化系数PUCI对置信度为90%的区间预报结果进行评价，并采用覆盖率判定系数CRC对概率预报结果进行整体评估. 概率预报可靠度评价结果见表6. 图3点绘了验证期4场洪水置信度分别为90%、85%、……、10%时，对应的区间覆盖率.