连续API产汇流模型的改进及应用

2018-01-04

水利水电快报 2017年12期

关键词：基流径流降雨

(长江水利委员会长江上游水文水资源勘测局，重庆 400020)

试验与研究

连续API产汇流模型的改进及应用

曹磊陈川建

(长江水利委员会长江上游水文水资源勘测局，重庆 400020)

结合长江流域上游水文特性和实际预报工作，采用连续的日消退系数计算方法、连续的地下水分割方式、自动计算降雨径流相关曲线(P～Pa～R曲线)和瞬时单位线计算方法，并结合降雨径流相关计算流程进行产汇流计算和评定，得到一套较完整的连续API降雨径流率定参数。通过该方法将API模型归纳成11个数字参数，方便利用智能优化算法进行优化率定计算，提高了API洪水预报模型的率定精度。提出的连续API模型的计算方法，为降雨径流相关模型方案的建立、率定和实际应用提供了一套新思路。

产汇流计算；连续API模型；蓄满产流；模型参数率定

API模型区间产汇流模型作为一种洪水预报模型在大多数流域得到广泛使用。该模型的参数率定涉及到场次洪水的选峰、地下水分割、净流量计算、时段单位线率定等多方面，率定过程较复杂，主要靠人工结合计算机交互技术，采用单场次洪水进行交互分析和数据预处理才能完成，特别是利用人工退水曲线的方式来分割场次净雨具有一定的随意性。随着地下水的分割[1]和还原[2-4]方法的使用，以及利用瞬时单位线取代时段单位线，逐渐与计算机技术结合后，提出了连续API模型并得到改进[5-6]，已在区间流计算方面得到推广和应用[7]，取得了一定的效果。然而，目前还是通过选择场次洪水人工交互分析得到连续API模型中退水曲线及场次洪水的分割，以及核心参数(P～Pa～R曲线)。同时，日消退系数K值，仍然靠人为来确定，这些都约束了连续API模型的率定和自动计算的应用，没有真正达到连续率定和连续计算的目的。本文以蓄满产流原理为基础，结合地下水分割和分布水文模型中的植被指数等，结合计算机技术提出了一套真正意义上连续率定、连续计算的连续API模型，包括自动进行地下水还原计算、P～Pa～R曲线率定、计算瞬时单位线和每月K值的计算过程。

1 模型介绍

API模型是前期降雨指数(Antecedent Precipitation Index)模型的缩写，其使用的表达式是传统的降雨径流相关图，故又称为降雨径流经验相关法。在20世纪40年代，R.K.林斯雷等提出了以季节、历时、前期降雨指数API为参数的五变数降雨径流相关图，配合时段单位线技术，成为一种有效的降雨洪水预报模型。1969年，西纳(Sittner)等提出了模型地下径流的建议，结合谢尔曼时段单位线即构成了可模拟流域降雨净雨过程的API模型。

API模型属于多输入、单输出静态的系统数学模型。模型建立之初，主要用于一次洪水径流量预报。在我国，模型经简化，形成了四变数关系线、三变数关系图(P～Pa～R)或两变数关系图(P+Pa～R)。简化后的API模型便于方案的率定和使用。其降雨产流函数关系可表示为：R=f(P，Pa，T，季节) 、R=f(P，Pa，T) 、R=f(P，Pa)或R=f(P，P+Pa)。其中，北方流域使用较多P+Pa～R产流方案，而南方、特别是长江流域多采用P～Pa～R三轴曲线。API模型的主要参数有：最大初损(IM)、前期降雨指数(Pa)、日消退系数(K)、降雨净流相关曲线(P～Pa～R)和流域时段单位线(UH)，计算时采用流域平均降雨来推算产汇流，参数的求解和含义可参考其他相应文献。

在模型使用过程中，水文预报工作者对该模型进行了不断改进[6]，包括Pa计算方式的改进，并将分块计算方式引入到了API模型中[6]，考虑到降雨分布不均匀情况下的产流计算；将API的计算方式与计算机结合提出API模型的连续计算方法[7]；将纳西(Nash)瞬时单位线(IUH)引入模型计算中，利用调蓄滞时(Ku)和串联水库个数(N)两个参数来确定流域瞬时单位线，再根据流域面积和计算时段长也可将瞬时单位线转化成时段单位线，使得单位线的率定更方便。

2 模型的不足与改进

不管是传统的API模型还是改进后的模型，在进行模型率定时，主要是靠人工结合计算机交互功能，经过场次样本选取、基流分割、径流计算、单位线推求、样本评定5个主要的复杂步骤来完成。特别是地下水分割时具有一定的随意性和不确定性(如图1)，使得模型的率定和实际使用过程中出现不一致的情况，这一定程度上制约了API模型的率定与使用，影响了模型计算精度。

图1 地下水不同分割示意

针对API模型的不足做了改进，提出如下连续API模型的率定和计算方法。

2.1 基流分割与还原

2.1.1 常规分割方法的局限性

在日常的API模型参数率定过程中，通常采用图1中列出的方法来进行基流分割，分割方式方法有一定的随意性，计算的净雨量和实际净雨量有时差异较大。在API模型方案参数的实际使用过程中，降雨未结束洪峰未出现时(预报根据时间后的过程)，无法确定直接径流停止点(B点)(平割法除外)，同时也很难确定计算的直接径流和地下水过程的总流量过程。不少使用者在直接径流的基础上加一个基流定值作为流域出口总出流，相当于不管采用哪种分割方式确定的方案都当成平割法进行处理，这种处理在一定程度上降低了预报的精度，计算流程也不是太合理。

2.1.2 数字滤波

滤波起源于通信，它是将信号中特定波段频率滤除得到有用信号的一种技术。内森(Nathan)和麦克马洪(McMahon)在1990年首次将数字滤波技术应用到基流分割中，该方法已成为近年应用最广泛的基流分割方法[8]。分割方程如下：

(1)

Qb2=Q2-Qd2

(2)

式中，Qd2和Qd1分别为滤波得到的第2时刻和第1时刻的地表径流，m3/s；Q2和Q1分别为第2时刻和第1时刻的径流，m3/s；Qb2为第2时刻的基流，m3/s；α为滤波系数。

这是一个递推公式，基于较长的系列资料进行计算时，特别是枯季开始时，Qd1比较小，可设置一个初值，利用公式(1)计算得到Qd2，再利用公式(2)计算得到第2时刻的基流，以此类推计算得到整个洪水过程的基流。这种方法在水文中得到不断的应用和推广，取得了较好的效果[9-11]。

2.1.3 地下水还原计算[3]

1969年美国天气局提出API模型中地下水出流采用马斯京根法线性水库(x=0)演算，模型中假定入渗到地下水库中的下渗量f与直接径流Q呈线性关系，即：

f=Z(Q-Qg)

(3)

地下水汇流用马法演算，可以得到如下公式：

Qg2=C0(f1+f2)+C2Qg1

(4)

将式(3)带入式(4)得到如下地下水推算公式：

(5)

式中，Qg1为上时段初地下水流量，m3/s；Q1、Q2为时段始末的直接径流量，m3/s；C0、C2为马法演算参数；Qg2为要计算的时段末地下水流量，m3/s；Z为地下水入流比例系数。

该方法基于直接径流进行地下水还原计算，避开了地下水的分割难题，在实际预报过程中采用公式(5)即可推导得出。

2.2 P～Pa～R参数计算改进

根据新安江三水源模型的蓄水容量曲线和降雨径流关系的诠释，蓄水容量曲线与降雨径流具有对应关系，故流域蓄水容量曲线的线形经确定后降雨径流关系曲线便随之确定[3]。蓄水容量曲线的线性在湿润地区广泛采用抛物线，因此,给定流域最大平均蓄水容量Wm(API模型里为IM)和抛物线的指数b，根据不同的前期蓄水量(在蓄满产流中为W0，在API模型中为Pa)，可对应生成一簇P～Pa～R曲线或一条P+Pa～R曲线，把复杂的关系曲线族归结于IM、b两个参数。

根据蓄满产流的原理，只要有降雨，根据蓄水容量分布曲线总有部分流域产流。但在实际工作中发现，在降雨较小的情况下，流域基本无产流或产流较小,连续API模型中设置了初始损失降雨Ps，使得在Pa较小时，曲线直接交Y轴与最大初损，Pa越大，初损逐步递减；有的流域土壤厚度较大，很难蓄满达到全流域产流的情况，在连续API模型中还设置了将整个曲线族按照曲线与Y轴交线进行旋转的参数β，得到最后相关线族。归纳起来，P～Pa～R曲线簇由流域最大缺水(IM)、蓄水容量曲线指数(b)、流域最大损失(Ps)和旋转角度(β)4个参数可推导出来。

2.3 K值的改进

早期的日消退系数K采用一个常值(通常取值0.85左右[4]),为了反映不同季节不同蒸散发，后来多采用每月一个值，在一定程度上更趋于合理，但相当于把原来单独的一个参数无形中变成了12个参数，增加了模型率定的难度。从K值反应的蒸散发的含义出发，力求算法简便。K.本旺(Keith Benven)[12-13]归纳了前人的研究成果，整理得到估算潜在蒸散发的季正弦曲线公式，公式如下：

(6)

另外，程根伟[14]等提出物候指数公式：

(7)

式中，Day为时间日序列的日数，即从年初开始的天数，d；A、B分别为与植物有关系的系数；91为水文年与自然年的日差值。

根据以上理论和API模型的日消退系数的含义，将该模型的日消退系数公式定义为：

(8)

图2 日消退系数变化

2.4 API模型中Pa数据计算模式

最初的API模型中Pa值按照场次降雨进行计算。在连续API模型中，利用计算机连续计算时，对于连续分散的小雨的情况，很难确定降雨开始时间及Pa值，这也会给预报计算带来一定的误差。连续API模型沿用原API模型的日模型计算结构，进行每日滚动计算日Pa值。每日切换一个Pa值对应的一条P～Pa～R关系曲线，每时段采用累积降雨进行查线的方法计算。

2.5 改进后的结构

图3 连续API模型结构示意

该流程利用多个数据化参数，概括了连续API模型结构和率定、计算的顺序，可利用遗传算法等AI算法进行优化计算，得到较为精确的模型参数。

表1 连续API模型优化计算参数

注：采用优化计算方法自动率定不包括滤波系数α。

3 连续API模型的应用

连续API模型可以应用人工逐步率定，也可以利用优化算法进行自动优化计算。

3.1 人工交互率定

人工交互率定可以分为两个步骤完成。

(1) 利用人工进行交互率定时，可先通过调整滤波系数(α)，将整个过程分成地下径流和直接径流，再调整地下水马法参数(Kmas)和地下水入流比例系数(Z)，使得模拟的地下径流跟分割的地下径流过程吻合且具有一定的合理性。

3.2 优化计算

连续API模型11个参数(自动率定为10个参数)都已数字化，可以直接利用如遗传算法等优化计算方法来进行优化率定。优化计算的适应度函数除了利用预报过程的确定性系数外，还可以结合计算过程的洪峰高度和洪峰时间误差等系数统一进行约束，达到洪峰高度和长系列过程都能满足要求的最优参数。也可以先通过人工交互的方式调整好几个不敏感的参数并进行固定，其他剩下的参数再利用优化计算法进行优化计算，最终得到合理的优化模型参数。

3.3 应用实例

以长江上游某流域为例，其集水面积1 441 km2，采用雨洪对应较好的1990年洪水摘录资料，计算时段长Δt=1 h，利用遗传算法，设置50个基因个体，对该连续API模型进行优化计算。采用完全确定性系数[15]控制和多洪峰(峰时和峰量)结合确定性系数控制(三者权值分别为20%，20%和60%)两种方法分别进行优化计算，每种方式各计算10次，最后对计算结果进行分析。

分析得出，按照确定性系数控制条件每次计算得到的最低指标都在乙级以上，按照多目标权值控制计算模式计算结果适应度函数平均值有所提高，但整个长系列的确定性系数稍有减小，对于洪峰的高度控制较好。计算结果见图4，模型参数如表1所示。

图4 连续API模型优化计算结果

利用该表中优化后的参数，逐年对1990～2012年23 a长系列资料进行检验，方案合格率达到乙级的年份所占百分比在75%以上，效果良好。

4 展望与不足

连续API模型采用11个数字参数很好地模拟了API模型的率定和实际使用计算，也为采用如遗传算法等优化计算方法进行API模型参数的优化计算提供了可能，一定程度简化了API模型率定的过程。通过控制优化计算的控制条件，采用传统的“确定性系数”结合洪峰高度和洪峰时间的适应度函数，提高并优化了计算的精度。

但由于连续API模型对一些参数进行了概化，也存在以下一些不足之处。

(1) 瞬时单位线由调蓄滞时(Ku)和串联水库个数(N)决定，在整个汛期过程采用的都是一条单位线，不能对单位线进行分类和分组。

(2) 公式(3)中，下渗量f和直接径流成比列关系只是一个假设，实际上f的量值既与土壤下渗特性有关，还与包气带调蓄作用有关[3]。

(3) 该方法只能提供三参数和两参数的API模型计算，未能考虑四变数和五变数的模型。

(4) 连续计算后，计算结果为某个流域的整个汛期产汇流过程，不是单场次洪水样本，只能使用确定性系数法来评定。若使用该连续API模型还需要进行场次评定，也可在整个长系列计算完成后，选择每场降雨和对应的洪水过程作为样本，设定样本的开始时间和结束时间对每场洪水的预见期、变幅、净雨、洪峰和峰现时间等进行场次精度评定和统计。

(5) 由于受到优化计算算法的限制，每次得到的参数不一定为最优，计算完成后最好进行合理性分析，并进行多次计算，以选择精度高并且又合理的参数为最后结果，在多次计算结果差别不大的情况下也可取平均值，需要具体问题具体分析。