基于参数分级优化的水文气象耦合预报研究

王秀杰，李丹丹，苑希民，徐浩田，齐喜玲

(1.天津大学水利工程仿真与安全国家重点实验室，天津 300072；2. 天津大学建筑工程学院，天津 300072)

山区中小流域在短历时强降雨下产汇流时间短，洪水成峰快[1-3]，同时流域自身调蓄洪水的能力较小，常常引发泥石流、滑坡等地质灾害，破坏力巨大，严重制约了我国经济社会的可持续发展．洪水预报作为一项重要的防洪非工程措施，可提前预测水库、河道等重要防洪控制点的水位、流量等信息，在防洪减灾中发挥的重要基础作用越来越明显．然而，山区中小流域由于缺乏准确的地面覆盖和气象驱动数据集以及成熟的模型参数转移技术，洪水预报面临着流域气象水文表征及流域响应表征皆具非平稳性的挑战[4]. 其中，模型参数的不确定性是中小流域流域响应非平稳性的重要表征．随着PUB(无资料流域水文预报)计划的实施[5]，模型参数在空间异质方面的推求取得一定进展，空间邻近法、回归分析法以及物理相似法等[6]参数区域化方法应运而生，但中小流域模型参数的时间异质性却很少被考虑．降雨的时空变异性导致了中小流域突发降雨过程中传输损失的不确定性[7]、模型参数的非平稳性[8-10]等问题，因此，仅通过率定一组水文模型参数来寻求一个流域径流形成的一般性或平均化规律[11-13]对全流域洪水进行预报必然存在较大误差．

中小流域水文气象表征的非平稳性主要体现在输入资料的短缺[14-15]．地面天气站点数据和卫星监测数据为流域水文规律的正确识别发挥了重要作用，但由于站点地理分布不均，某些数据未能观测到表征天气的全部特征或因受限于观测仪器性能(如噪声数据)等问题很可能造成预报的振荡与波动．同时，世界气象组织(WMO)数据显示的世界范围内水文数据观测网络的逐渐减少[16]也无疑对水文预报精度提出挑战，特别是受社会经济、技术及水资源战略地位影响的发展中国家．针对以上问题，气候再分析将过去的观测结果与模型结合起来以生成多个气候变量的一致时间序列，在统计意义上最真实地再现了当时的大气状态．同时，由于水文站网的限制及预见期延长的需求，目前基于卫星、雷达等遥感技术的定量降水估算和基于数值模式的定量降水预报(QPFs)在洪水预报领域不断推广应用，促进了水文气象耦合预报研究的不断深入[17]．交互式全球大集合预报系统的集合预报(TIGGE)将“单一”的确定性预报转变为多值概率预报，可定量估计天气预报的不确定性，为资料缺乏地区的洪水预报提供了气象数据来源途径[18-19].

基于以上问题，本文提出一种基于降雨时程分配指标(最大降雨强度和分配方差)的参数分级优化方法，同时引入ERA5 再分析数据和集合预报产品，对广西驮英水库上游洪水开展预报研究，为中小流域洪水预警预报和防御提供客观事实依据，对提高应对极端洪水的决策支持能力具有重要意义．

1 研究方法

1.1 基于降雨时程分配指标的参数分级优化方法

SCS-CN 模型因其计算过程简单、所需参数较少，被广泛应用于缺资料地区的地表径流估算中[20].模型地表产流中仅有一个反映流域土地利用、土壤类型、前期土壤湿度及水文条件等特征的综合径流曲线数CN；地下汇流采用线性水库贮水结构进行计算，根据退水过程分析得到地下水库的蓄泄系数后，仅有平均下渗率f 需要率定．目前，国内外针对SCS-CN模型的研究主要集中在以下3 个方面[21]：①区域径流量估算；②模型结构改进；③模型参数率定与改进．其中，基于降雨-径流观测数据率定CN 值的研究较为广泛，具体表现为：①根据降雨-径流反推得到相应的潜在蓄水能力S 值与流域综合CN 值，计算出CN 值与降雨量之间的函数关系[22]，或进一步用平均值法、对数频率分布法、渐近线法和中值法等推算综合CN 值[23]；②用参数自动优选确定流域综合最优CN 值．研究表明地区不同，CN 的取值范围也不同：在黄土高原干旱半干旱气候地区，AMCⅡ条件下CN取值基本在70～80 范围内，在不同土壤湿润条件下CN 取值变化范围大，如北京密云石闸自然植被条件下的CN 取值在57.7～89.3 变化范围内[24]；南方地区的CN 取值总体比北方大，如湖北省松柏站区域在场次降雨下CN 值基本在90～97 范围内[25]；江西吉泰盆地红壤丘陵区的CN 值在AMCⅠ条件下的变化范围为86.57～91.48[26]．同样地，平均下渗率f 在不同场次降雨下也呈现不同特性，如章俊霞等[27]在研究影响南方红壤土入渗因素时，得出多场次降雨下平均下渗率主要集中在1～10 mm/h 范围内．参数CN 和f与降雨总量及降雨强度密切相关．由于降雨总量受降雨历时的影响，暴雨侵蚀能力的非平稳性将对参数的分级率定造成影响．因此，引入降雨时程分配指标，将降雨总量对参数的影响分解为最大降雨强度(mm/h)与降雨方差(降雨值与最大值间的方差(n 为降雨历时))两部分(见表1)，旨在通过考虑降雨时程分配引起的参数时间异质性的参数分级优化方法进一步提高洪水预报精度．

表1 参数分级优化参考表Tab.1 Reference table for hierarchical optimization of parameters

根据率定的不同CN 值利用式(1)～(3)推导流域的产流过程[28]．

式中：S为潜在蓄水能力，mm；Ia为初损，mm；P为降雨过程，mm．

同样地，根据率定的不同等级的平均下渗率f计算相应的地下汇流过程．

1.2 气候再分析数据集及修正

ERA5 是ECMWF 利用先进的建模技术和数据同化系统重新分析已有观测数据而得到的再分析数据集，包括每小时的大气、陆地表面和海况等参数.ERA5-Land 是一个9 km 分辨率的全球陆地表面数据集，与ERA5 相比，它以增强的分辨率提供了几十年来土地变量演变的一致视图，为缺资料地区的水文预报提供了更为精确的数据来源．与其他任何模拟一样，ERA5-Land 数据集提供的估算值具有一定程度的不确定性，数值模型只能或多或少准确地表示控制地球系统不同组成部分的实际物理过程，模型估算的不确定性会随着时间的推移而增加．若将气候模式结果直接应用于水文驱动，其偏差会对模拟产生很大影响．扰动法(平均态订正)[29]和分位数映射法(分布态订正[30])是气候模式结果订正常用的两种方法，其中，基于概率分布的分位数映射法在高时间分辨率上订正效果较好．分位数映射法中传递函数包括基于理论概率分布函数[31]和经验概率分布函数[32]，其中，基于非参数转换的经验概率分布由于不需要对原始数据做前提假设订正而得到广泛应用．

1.3 数值集合预报检验与评估方法

国际全球集合(the international grand global ensemble，TIGGE)是一个世界气象组织的项目，旨在提高一天到两周的恶劣天气事件预测的准确性[33].TIGGE 考虑到初始状态条件的扰动和数值模式中物理过程的不确定性，可捕获一系列现实的可能结果．集合预报数据采用Gridded Binary(GRIB) edition 2 格式进行存储，预报步长为6 h，预报最长时效可达16 d．该数据集由来自10 个全球NWPs(数值预报)中心的集合预报数据组成[34]，预报数据可在欧洲气象中心官网进行免费下载．

由于数值产品的不确定性会显著影响流量预测的准确性，因此，数值天气预报数据在应用前的检验与评估至关重要．目前NWPs 的预测评估方法可以分为两大类：①集成预报的概率评估，例如Brier Score(BS 评分)、Brier Skill Score(BSS 评分)[35]、相对作用特征(ROC)、Talagrand 分布等；②确定性预报的确定性评估方法，典型的有距平相关系数(ACC)、均方根误差(RSEM)、TS 评分等．其中，基于确定性评估的TS 评分及基于概率评估的BS、BSS 评分在气象数据可靠性评估中被广泛应用．

1.4 灰色关联分析法

灰色关联分析由于其对数据要求较低、原理简单，被广泛应用于进行各评价对象与参考对象之间关联度的优劣排序，量化两系统的因素动态发展态势的关联性大小[36]．下面选取降雨特征中的累积降雨量、最大降雨强度、降雨方差与流量峰值之间的关联性进行说明．为了量化不同预见期下累积降雨量、最大降雨强度及降雨方差对洪水量级的影响程度，建立流量、水位与最敏感因素之间的经验公式，便于流量的快速估算，对以上3 个因子进行灰色关联度分析(GRA)[37]．

1) 数据标准化处理

本文以流量值X0={X0(1), X0(2),…, X0( k)}( k =1,2,3,…, m)为参考序列，分别以不同预见期下的累积降雨量X1、最大降雨强度X2、降雨方差X3为比较序列．为消除量纲带来的数值范围影响，现对数列进行量纲归一化处理．具体表现为用每个数列第一个数 Xi(1)去除其他数 Xi( k )，再求各因子与参考数列的差值序列的绝对值

2) 比较数列与参考数列的关联分析

式中：ξi(k )为第i 比较序列中第k 个值与参考值的关联系数；为第i 比较序列与参考序列间的绝对差，Δi(max)和Δi(min)分别为最大值与最小值；分别代表两级最小值和两级最大值．

2 实例分析

驮英水库位于广西宁明县那堪乡垌中村上游约6 km 的珠江流域西江水系明江支流公安河上游河段．水库上游流域面积606 km2，研究区域经纬度范围 为21° 35′ 45′ ′N ～21 °55′ 5.9′ ′N ，107 °17′ 7 .1′ ′E ～107°34 ′3 .7 ′E ，属亚热带季风气候[38]．流域的土壤特性、水文分组、土地利用情况见图1．流域内83.87%面积的土壤属性为壤土，饱和导水率在1.8～18 mm/h 范围内，水文分组为C 组，剩余16.13%的流域面积水文分组为D 组．夏末秋季，受南海及太平洋台风影响，工程区多出现台风雨，强度大，雨量集中，具有典型的山区性河流特性．驮英水库未建前曾设有驮英水文站，目前已停测，现有原驮英水文站1987—1999年共计13 年水文气象资料．驮英水库上游建有九特、那驮、板固、叫弄、枯强、潭昔、小平7 个雨量站，现仅有2012—2019 年实测资料．降雨径流数据在时间上的不对应给模型的率定带来困难．现引入ERA5-Land 降水再分析逐小时数据集，水平分辨率为0.1°×0.1°．从原驮英水文站记录的1987—1999年历史流量数据选取多个流量量级的场次洪水，并在ERA5-Land 数据集中下载相应时间的栅格降水数据，通过算数平均法计算流域的面平均降雨量．利用分位数映射法订正的流域面降雨量和驮英水文站相应时间的实测流量分别对驮英水文站的实测降雨过程进行强度和分布的校正，最终使用校正后的降雨数据开展模型的率定和验证．同时，由于目前坝址以上流域尚未建设水文、水位测站，驮英水库流域洪水预报结果的校核数据则需根据实测水位数据(水位尺测量)与水位-库容关系、水位-下泄流量关系进行反演计算，入库流量反演计算步骤见图2．

图1 驮英水库上游流域土地利用和土壤类型分布Fig.1 Distribution diagram of land use and soil types in the upstream basin of Tuoying reservoir

2.1 ERA5-Land降水再分析数据适用性检验与修正

利用非参数转换RQUANT 建立水库上游7 个雨量站2012—2019 年间场次降雨平均值与对应ERA5-Land 逐小时再分析降水数据之间的传递函数，并将传递函数应用于1987—1999 年间ERA5-Land 场次降水分布与2012—2019 年间ERA5-Land降水分布相似(强度及分布)的降雨数据的偏差修正，利用修正后的ERA5-Land 降雨数据和逐小时流量数据对原驮英水文站的降雨数据分别进行降雨强度和降雨分布的校正，最终获得修正后的1987—1999 年用于模型率定和验证的降水数据．以2013 年第30号台风“海燕”(1979 年以来登陆广西省最强台风之一)形成的暴雨过程为例，通过分位数映射得到的偏差订正结果见图3，修正前后的误差泰勒图见图4(A点为修正前ERA5-Land 降水预报性能；B 点为修正后ERA5-Land 降水预报性能；Buoy 为实测数据).从图中可以看出，修正后的降水真实性能得到有效提高.

图2 入库流量反演计算流程Fig.2 Flow chart of inversion calculation process of discharge

图3 降水偏差订正结果Fig.3 Corrected results of precipitation deviation

图4 降水修正前后泰勒图评估Fig.4 Taylor chart assessment before and after correction of precipitation

2.2 参数分级优化率定结果

选取1987—1996 年订正后的降水和驮英站流量资料对模型中的CN 和f 两参数进行精细化率定，利用1997—1999 年数据进行模型验证．首先计算场次降雨的时程分配特征(最大降雨强度和分配方差)，并利用不同的降雨类型数据对模型参数进行率定和验证，结果显示不同的降雨类型率定结果不同．最大降雨强度对参数率定结果影响最为显著，说明了短历时强降雨过程对洪水产流机制的重要性；降雨分配方差对参数的影响使其在相邻值范围内取值．根据率定的不同结果选取对应的相邻降雨特征值区间之间的中间整数作为级别划分阈值．经率定、验证得到模型参数，结果见表2，其中CN 值范围在95～98，与湖北松柏站区域的CN 值范围(90～97)较接近[25]，参数f 的率定结果与南方红壤土入渗研究得出的1～10 mm/h 范围较接近[27]．由于研究区域地处湿润地区，降雨前期土壤的湿润程度在降雨初期很容易达到饱和，对径流产生的影响很小．因此，参数CN 值主要受下渗截流的影响，平均下渗率f 越大，相对的地表产流越小，导致CN 值越小．因此，可根据二者的对应关系由其中一个参数得到另一个参数．二者之间的具体关系见表3．

表2 参数CNi 与 fi 精细化率定表Tab.2 Precise calibration table for parameters CNi and fi

表3 f i与CN i值对应取值表Tab.3 Corresponding vaule table for fi and CNi

2.3 集合预报数据的检验与评估

选取2016—2018 年夏季(6—8 月)TIGGE 4 套集合预报产品(ECMWF、NCEP、UKMO、JMA)的集合日平均降雨数据作为评估数据．起报时间为每日0时，输出间隔24 h，预报时效为1～10 d，空间分辨率选择最高0.5° ×0.5 °；实测降水数据来自流域上游7个雨量站记录的2012—2018 年对应的日累积降雨量，时间标准采用世界时．由于研究区域基本分布在数据集预报范围0.5° ×0.5 °内，因此不需要对数据集中的降雨预报数据做降尺度处理，可直接将其作为研究区域的预报面雨量．利用TS 评分、BS 评分、BSS评分分别对ECMWF、NCEP、UKMO、JMA 4 个气象中心对该场降雨的集合预报产品进行可靠性检验．

为提供TS 评分与BS 评分过程中的分级可靠性检验依据，首先将累积降雨量按照《中短期天气预报质量检验办法》对降雨量的划分原则进行分类统计，具体等级划分标准见表4 所示．

表4 降水等级划分标准Tab.4 Classification criteria for precipitation

TS 评分与BS 评分的取值均在0～1 范围内，TS值越大代表预报效果越好，TS＝1 代表全部命中．相反，BS 值越小代表对降雨事件的预报效果越佳，BSS越大表示预报效果越好．以集合预报平均值作为预报结果的各产品TS 评分及BS 评分结果分别见表5和表6．

表5 ECMWF、NCEP、UKMO、JMA 针对各级降雨预报的TS评分Tab.5 Threat scores of rainfall forecast at several levels for ECMWF，NCEP，UKMO，and JMA

表6 ECMWF、NCEP、UKMO、JMA 针对各级降雨预报的BS评分Tab.6 BS scores of rainfall forecast at several levels for ECMWF，NCEP，UKMO，and JMA

图5 预报产品综合评价性能Fig.5 Comprehensive evaluation chart of forecast products

图6 2016—2018年夏季各模式日降雨数据的Brier技巧评分Fig.6 Brier skill scores for daily rainfall data in summer of each model from 2016 to 2018

表5 定量揭示了4 种预报产品在不同降雨等级上的预报效果均表现出随着降雨量级的增大命中率逐渐减小．具体而言，在中雨预报范围内JMA 表现更为突出，NCEP 在大雨级别上的预报效果相比其他产品结果较差．表6 中BS 评分结果说明NCEP 对中雨、JMA 对大雨的捕捉概率优于其他产品，暴雨级别上ECMWF 和NCEP 的预报能力更突出．为直观考量4 个预报产品的综合预报性能，现计算BS 值的相反值(1-BS)，并分别赋予两评价指标0.5 的权重综合计算4 个数据中心的预报性能，计算结果见图5．从图5 中可以看出4 个气象中心在小雨预报上均表现最好，其次为暴雨、中雨和大雨；4 种预报产品的预报性能在不同降雨等级上均有不同程度的差异；NCEP、ECMWF 的预报效果分别以综合性能0.50 和0.49 在整体预报上更胜一筹，其次为UKMO、JMA.为进一步比较各集合预报在预报时效范围内日尺度上的预报性能，利用BSS 评分对预报产品1～10 d 预报时效内的预报结果进行评估．从图6 中可以看出各模式的预报性能随着预报时效的延长而逐渐变差.1～6 d 预报中NCEP 的预报效果最好，ECMWF 预报效果与之相差甚微，其次为JMA、UKMO．7～10 d的预报中，ECMWF 预报效果最好，NCEP、UKMO、JMA 次之．综上，ECMWF 与NCEP 的预报效果始终最好，与图5 的评估结果一致．6 d 内的短期预报可选用NCEP 预报中心的集合预报进行驮英水库流域的洪水预报．根据NCEP 逐6 h 的数据产品及研究区域的典型雨型，将每6 h 的累积降雨按照计算的比例进行逐小时降雨量分配得到最终的降雨预报结果，见图7．从图7 中可以看出，实测降雨的雨型分布与研究区域的设计雨型分布大体相同，预报过程与实测过程总体趋势与峰现时间均表现一致，峰值误差为14.06%，纳什系数为0.87，并且在逐小时不同降雨量级的预报效果上均表现良好．

图7 NCEP预报产品与实测数据对比Fig.7 Comparison between NCEP forecast products and measured data

2.4 洪水预报结果分析

2019 年8 月1 日22：00 第7 号台风“韦帕”(热带风暴级)进入广西省北部湾后沿海岸线西行．2 日21：20 前后在防城港市沿海再次登陆．受其影响，驮英水库上游流域产生大范围降雨，局部地区达到暴雨级别，且持续时间较长．根据NCEP 集合预报产品最终确定此次过程的降雨历时分布见图8．由于工程区尚未建水文站，入库流量可由实测水位数据进行反演复核．根据图8 的降雨分布确定参数取值f＝9，CN＝95，预报计算结果见图9．从图中可以看出计算结果与复核结果的总体趋势与峰现时间均表现一致，流量峰值与降雨峰值峰现时间相差6 h，预报洪峰流量为1 531 m3/s，复核流量峰值为1 416 m3/s，峰值误差为8.12%；计算结果与复核结果的平均相对误差Re＝13.09%，纳什系数NSE＝0.86，模型计算结果的相对误差均小于20%，满足《水文情报预报规范》中的误差要求．

图8 基于NCEP集合预报的“韦帕”台风降雨过程分布Fig.8 Rainfall distribution of typhoon “Weipa”based on NCEP ensemble forecast

图9 计算结果与复核结果对比Fig.9 Comparison between computational results and review results

2.5 洪水敏感性分析与流量估算

对于中小流域，降水的时间分辨率比空间分辨率更敏感[39]．因此，评估降雨变化过程如何影响预测流量过程和引起的不确定性问题不容忽视，尤其气候变化可能会导致极端事件的频率增加[40]时．洪水量级对累积降雨量W、最大降雨强度P以及降雨分配(方差) S2较为敏感，利用灰色关联分析法分析3 个主要影响因子与洪水量级的关联度，结果见表7．可以看出，3 个影响因子与流量峰值的关联度最大值分别为0.917、0.899、0.797，最小值分别为 0.905、0.879、0.770，说明洪水量级与3 个影响因子有很强的相关性；不同预见期下均表现为累积降雨量对洪水的影响最为显著(最大关联度为 0.917，最小关联度为0.905)，其次为最大降雨强度(最大关联度为0.899，最小关联度为0.879)；3 个预见期下累积降雨量、最大降雨强度与洪水量级的关联度差别均很小，从而说明了短历时强降雨过程对洪水的产流机制至关重要；同时，3 种预见期下相应影响因子的关联性值相差甚微，说明3 个影响因子与洪水量级的关联性对预见期的敏感性很低．

为了解不同预见期(6 h、12 h、24 h)下的累积降雨量与洪峰流量及坝前水位之间的关系，将累积降雨量与流量及坝前水位数据绘制在坐标图中，并用曲线拟合，结果显示洪峰流量、水位与累积降雨量的经验公式皆可用不同形式的幂函数进行表示，见图10．根据经验公式可依据累积降雨量实现相应的流量与坝前水位的快速动态估算，可为应急抢险争取了更多宝贵时间．

表7 灰色关联分析结果Tab.7 Results of grey relational analysis

图10 不同预见期下累积降雨量与洪峰流量、坝前水位的关系Fig.10 Relationship between cumulative rainfall ，peak flow and water level in front of the dam under different forecast periods

3 结 论

本文针对模型参数非平稳性导致的仅率定一组平均水文参数造成预报误差问题，提出了基于降雨时程分配指标(最大降雨强度和方差)的参数分级优化方法，同时引入降雨再分析数据(ERA5-Land)和集合预报产品进行洪水预报，并利用灰色关联度量化降雨过程中的累积降雨量、最大降雨强度及分配方差等不确定性对径流产生的影响，建立流量、水位与最敏感因素之间的经验公式．主要结论如下．

(1) 洪水预报模型中引入最大降雨强度及降雨方差的参数分级优化方法综合考虑了降雨的分配引起的产流规律的变化，将其用于实例分析中的平均百分比误差为13.09%，纳什系数为0.86．可见该方法明显减少了参数不确定性对水文过程的影响，洪水预报精度得到显著提高．

(2) TIGGE 集合预报对资料短缺地区的气象水文预报提供了重要的气象数据来源，其提供的多值概率预报定量估计了天气预报的不确定性，进一步保证了预见期延长下的流量预测精度．但TIGGE 中的集合预报产品对各地区的适用效果各不相同，在研究区域的洪水预报过程中，通过多种产品对比分析选取预报性能最好的集合预报产品，从模型输入角度进一步减小洪水预报的不确定性．同时，气候再分析数据集也为缺资料地区的洪水预报提供了可靠的数据来源.

(3) 在累积降雨量、最大降雨强度及分配方差3个影响因子中，洪水对累积降雨量的敏感性最为显著．同时，不同预见期的洪峰流量、水位与累积降雨量之间存在确定的幂律关系，通过这些关系可快速预估水库流量与相应水位，为水库进行实时洪水调度提供可靠依据．