HEC-HMS水文建模系统原理·方法·应用
上QQ阅读APP看本书,新人免费读10天
设备和账号都新为新人

第9章 HEC-HMS模型校验

本章介绍HEC-HMS中的模型校验的概念和参数搜索优化方法。HEC-HMS模型校验就是使用实测的水文气象数据搜寻计算系统的参数值,当搜寻计算结果与实测径流吻合良好时确定模型的参数值。这种搜通常也被称为为优化。首先介绍了HEC-HMS模型校验的概念,之后总结了HEC-HMS模型校验方法,提出了模型校验的符合度指标,最后介绍了单变量梯度搜索算法和Nelder and Mead算法。

9.1 校验的定义

每一个包含在HEC-HMS中的模型都含有参数。使用模型进行径流估算或水文过程线演进时需要指定参数的值。前面的章节认识了这些参数并介绍了如何从集水区及河道的特性估算这些参数。例如,第6章介绍的动波直接径流模型中有一个反映地面粗糙度的参数N;这个参数可以从所了解的流域的土地利用情况估算。

但是,如在第2章所述,HEC-HMS中的一些模型含有无法从河道或集水区实测的特性预测的参数。斯奈德单位线模型的Cp就是一个例子;该参数没有直接的物理意义。同样,马斯京根演进模型参数x也是不能量测的;它只是一个简单的表示河道中上下游流量在计算时的相对重要性的权值。方程(8-18)提供了从河道特性估算x的方法,但是这只是一个近似,只适合有限的几种情况。

那么如何确定合适的参数值呢?如果可以得到实测的降水量和河流流量,那么可通过校验来确定合适的参数值。校验使用实测的水文气象数据搜寻系统的参数值,这些参数能使计算结果与实测径流吻合良好。这种搜寻常被称为优化。

9.2 HEC-HMS校验方法总结

在HEC-HMS中,是根据图9-1所示的过程系统地搜索最优参数值。该过程从数据收集开始。对于降雨-径流模型,所需数据是降水和时间序列。对于演进模型,需要演进河道实测的入流和出流量。表9-1和表9-2为数据的收集提供了某些提示。

表9-1 降雨径流模型校验的数据收集说明

图9-1 校验过程的概念图

表9-2 演进模型校验的数据收集说明

下一步是选择参数的初始估算值。如任何的搜索方法一样,初始估算值设置得好(搜索的开始),就会越快搜索到结果。在前面章节中对参数估算值的提示可能会很有用。

给定这些参数的初始估算值,使用实测的边界条件(降雨或上游流量),HEC-HMS的模型可以计算输出结果也可以计算集水区径流过程或河道出流过程线。

这时,HEC-HMS将计算的水文过程与实测的水文过程进行比较。例如,HEC-HMS计算图9-2中用点画线表示的水文过程线,并将它与用实线表示的实测的水文过程线径行比较。比较的目的是判断模型是否符合真实的水文系统。比较方法将在本章的后面介绍。

图9-2 计算与实测的水文过程线的比较

如果符合度得不到满足,HEC-HMS会系统地调整参数并进行迭代运算。调整参数的算法在本章的后面加以介绍。

当符合度得到满足后,HEC-HMS将报告最优的参数值。并假设这些参数值可以用于径流或演进计算,这也是洪水径流分析的目的。

9.3 符合度指标

为了比较计算过程线与实测过程线,HEC-HMS将计算符合度指标。包含在HEC-HMS中的用于搜索最优指标模型参数的算式也被称为目标函数。在HEC-HMS中,根据分析的需要,可以使用四个目标函数中的一个。这四种校验方法的目的是找出使目标函数的值最小的合理的参数值。可选的目标函数(如表9-3所示)是:

表9-3 HEC-HMS 校验的目标函数

Z为目标函数;NQ为计算的过程线纵坐标数目;q0i)为实测流量;qsi) 为用所选的模型参数计算的计算流量;q0(peak) 为实测峰值;q0(mean) 为实测流量的平均值;qs(peak)为计算的峰值。

1.绝对误差的和

目标函数将每一个计算的过程线的纵坐标与实测值进行比较,均匀地将权值赋予每一个值。这时,比较的指标就是纵坐标的差。但是,因为这个差值可正可负,简单求和可以使正负的差值相抵消。在水文模拟中,不希望出现正的或负的偏差,就像不希望高估或低估出现的情况一样。为此,该函数是差值的绝对值的总和。这样,这个函数间接地成为比较两个过程线的峰值、水量、峰值时间的符合度的一个度量。如果该函数值为零,则表示它们完全相符:即所有的计算的过程线的纵坐标完全等于实测值。但是这种情况很少见。

2.残差的平方和

这是模型校验常用的一个目标函数。该函数也比较所有的纵坐标,但用差的平方表示符合度。这样,10m3/s的偏差就比1m3/s的偏差差100倍。将偏差平方后可能会产生不必要的高估和低估。同样,该函数间接地成为比较两个过程线的峰值、水量、峰值时间的量值的一个度量。

3.峰值误差百分比

这个指标仅量测计算的峰值流量与实测的峰值之间的符合度。它将符合度量化为表示为差的绝对值的百分比。之所以这样处理是不希望低估和高估偏差。它不能反映水量或时间的误差。这个目标函数在需要对峰值流量和峰值水位进行限制的规划和设计中是一个合乎逻辑的选择。洪泛区的管理研究常就是这种情况,试图用流量和水位之间的关系来限制遭受淹没的区域的面积。

4.峰值加权均方根误差

该函数与计算程序HEC-1(USACE,1998)中的校验函数相同。它比较所有纵坐标,将偏差平方,并给平方后的偏差赋权值。分配到每一个纵坐标上的权值与纵坐标值成正比。比实测过程线的均值大的纵坐标分配的权值大于1.00,比均值小的赋予小于1.00的权值。峰值实测坐标被赋予的权值最大

除了符合度的数值度量,HEC-HMS也提供图表比较,这使得模型与实测水文系统的符合度的可视化。HEC-HMS显示的计算的水文过程线与图9-2所示的极为相似。除此之外,还可以显示散点图,如图9-3所示。该图是用每一个时间步长的计算值与同一时间的实测值绘出的图。通过检查这些图有助于识别由于参数选择而产生的模型偏差。图中的直线表示流量计算值与实测值相等:如果这些点在该线之上,表明使用所指定的参数的模型预测出的值正好与实测值相等。在直线上方的点表示模型预测值比实测值高。在直线之下方的点表示预测值小于实测值。如果所有的点都在该等值的直线上方,说明模型有偏差;因为它总是给出过高的估算值。同样,如果所有的点都在该直线的下方,这表明模型总是发生低估。如果点分布于该直线的上方和下方的数目相同,这表明该经过校验的模型发生高估或低估的可能性相同。

等值线周围点的散布情况也提供了模型符合度的一种表示。如果点的分布比较分散,说明模型不能很好与实测值相匹配,预测中的随机误差与流量值的大小有较大的关系。如果点的分散很小,说明模型参数符合较好。

图9-3 计算与观测流量散点图

HEC-HMS还计算并绘制流量计算值与实测值之间的残差的时间序列。图9-4是一个例子。该图表示估算误差在模拟期间的分布。检查该图有助于将注意力集中到估算时需要特别关注的参数上。例如,如果最大的残差集中在径流开始时,说明初始损失参数可能选择不当。

图9-4 流量残差图

9.4 各种搜索方法

如早前所述,HEC-HMS中的校验的目的是为了确定合理的模型参数,使用目标函数之一作为度量时,这些参数使计算的过程线与实测过程线之间的符合度最好。

如图9-1所示,这种搜索是一种试错法式的搜索。选择试算参数,使用模型计算结果,并计算结果的误差。如果误差不可接受,HEC-HMS就会改变试算参数并重复计算。根据单变量梯度搜索算法或者是Nelder and Mead单一搜索算法改变试算参数。

9.4.1 单变量梯度搜索算法

HEC-HMS中的单变量搜索算法可对参数估算值进行连续的修正。也就是说,如果xk代表目标函数fxk)第k回试算的参数估算值,那么在k+1回定义一个新的估算值xk+1

式中:Δxk为参数的改正量。

HEC-HMS中该搜索的目的是选择Δxk使其估算值逐渐接近可以产生目标函数最小值的参数值。一次改正一般不能达到最小值,因此方程会被反复地使用。

HEC-HMS中使用的梯度法的基础是牛顿法。牛顿法使用下面的方法定义Δxk

(1)用泰勒级数将目标函数近似为:

式中:fxk+1)为k次循环计算时的目标函数;df(·)/dx和d2f(·)/dx2分别为目标函数的一阶和二阶导数。

(2)最好选择xk+1使fxk+1)最小。这时fxk+1)的导数应等于零。为此,求出方程(9-2)的导数并使之等于零,同时省去高次项。得到:

该方程重新整理后与方程(9-1)结合得到

HEC-HMS使用了一种数值计算方法计算每一个循环k时的导数df(·)/dx和d2f(·)/dx2。计算如下:

(1)定义=0.99xk=0.98xk为与xk相邻的两个替代参数,并计算每一个参数的目标函数。

(2)计算它们之间的差值,得到

(3)导数df(·)/dx被近似为Δ1,导数d2f(·)/dx2被近似为Δ21。严格地讲,当这些近似表达式被代入方程(9-4)时,就得到了牛顿法中的改正值Δxk

正如在HEC-HMS中实现的那样,根据HEC研究人员的模型校验经验对该改正值进行了修正。即用下式进行修正:

其中C如表9-4所示。

除了这个修正,HEC-HMS对每一个xk+1进行测试,以确定是否满足fxk+1)<fxk)。如果不是,就定义新的公式试算:

如果fxk+2)>fxk),搜索结束,因为已不需要修改。

如果通过校验要找的不只是一个参数,那么通过令其他参数为常数,然后对每一个参数先后使用这一搜索方法。例如,欲求出斯奈德的Cptp,在初次估算时HEC-HMS保持tp不变并首先调试Cp。接着,保持新的调整过的Cp值不变,算法调整tp。这种先后的调整重复了四次。之后,该算法评价这些参数的最后一次调整,以确定能对目标函数产生最大减小的参数值。使用此处定义的方法调整这个参数。该过程一直持续,直到附加的调整使目标函数的减小值小于1%。

9.4.2 Nelder and Mead算法

Nelder and Mead算法不使用目标函数的导数进行最优参数值的搜索。相反,该算法依靠一种更简单的搜索。在搜索中,用下面的方法选择参数估算值,即根据前一步的试算获得的结果识别好的估算值并排除差的估算值,从用好的估算值建立的模式生成更好的估算值。

Nelder and Mead搜索法使用一个单纯形,即一组替代的参数值。对于一个有n个参数的模型,该单纯形有着n+1个不同集合的参数。例如,如果该模型有两个参数,这两个参数的每一个参数的由三个估算值组成的集合包含在这个单纯形中。

几何上看,可将n个模型参数视为空间的维数,单纯形作为n维空间中的多面体,每一个参数集合作为该多面体的n+1个顶点。在双参数模型情况,单纯形就是一个二维空间的三角形,如图9-5所示。

表9-4 HEC-HMS单变量梯度搜索法改正系数

图9-5 二参数模型的初始单纯形

Nelder and Mead算法改进单纯形以找出一个顶点,在这个顶点上目标函数的只是最小的。为此,该方法进行下面的运算:

(1)比较。

改进的第一步是找出产生最坏(最大)的目标函数值的单存形的顶点,以及产生最好的(最小)的目标函数值的单存形的顶点。在图9-6中,分别用WB标记表示。

(2)对称。下一步是找出所有顶点的质心,不包括W点;该质心用图9-6中用标记C表示。算法接着从W定义一条通过该质心的直线,沿着该直线对称距离WC定义一个新的顶点R,如图9-6。

图9-6 单纯形的反射

图9-7 单纯形的扩展

(3)扩展。

如果用点R代表的该参数集合比该最好的顶点要好,或同样好,算法沿同一方向进一步扩展该单纯形,如图9-7所示。这就在图中定义了一个扩展的顶点E。如果该扩展的顶点比做好的还要好,就用这个扩展的顶点代替单纯形最差的顶点。如果扩展的顶点不必最好的好,那么就用反射的顶点代替最差的顶点。

(4)收缩。

如果反色顶点比最好顶点差,但是比某些其他的顶点要好(不包括最差的),通过替换最坏的顶点收缩该单纯形。如果反射顶点不比其他任何点好,排除最差点,收缩单纯形。这一过程如图9-8所示。为此,最差点沿着直线向质心移动。如果该收缩点的目标函数较好,就用这个点替代最差点。

图9-8 单纯形的收缩

图9-9 单纯形的减小

(5)减小。

如果收缩顶点没有改进,则单纯形通过向最佳点移动所有的顶点而被减小。这就会产生新的顶点R1R2,如图9-9所示。

当满足下面的准则时,Nelder and Mead的搜索就会结束:

(1)

式中:n为参数个数;j为顶点下标;c为质心点下标;zjzc分别对应顶点jc的目标函数。

(2)试算次数达到参数个数的50倍。

当搜索结束由最好顶点代表的参数就会被报告为最优参数值。

9.5 搜索的约束条件

系统工程师们将找出模型的最优参数值的数学问题称为约束优化问题。也就是说可行的范围和可接受的参数(被称为决策变量)是受到限制的。例如,马斯京根x参数小于0.0或大于0.5是不能接受的,不管结果如何一致。这样。在此范围外的搜索是没有必要的,在此范围以外找到的之都是不可接受的。对x的限制以及其他列于表9-5中的限制在搜索是要加以考虑。

在使用单一变量梯度或是Nelder and Mead算法的搜索中,HEC-HMS检查每一次试算以保证试算的值在可行的范围之内,如果不在这个范围,在继续计算前,HEC-HMS将增加试算值使其达到最小或减少试算值使之达到最大。

除了这些不能违反的约束,HEC-HMS也将考虑用户指定的软件约束条件。这些约束定义了对参数的限制。例如,常损失速率值缺省的可行性范围是0~300mm/h。但是,对于一个由密实的黏性土组成的排水区,这个速率可能小于15mm/h——使用比它大的值将受到质疑。需要指定15mm/h作为所需的土壤约束条件。这样,搜索在此范围外产生一个候选的参数时,目标函数就会被乘以一个罚因子。该罚因子的定义是:

式中:xi为估算的参数值ici为参数i最大或最小值;n为参数数量。该系数“规劝”搜索算法在软件约束范围内选择参数值。例如,如果均匀损失速率的搜索导致300mm/h的值,而软件指定的约束是15mm/h,目标函数的值将被乘以2(300-15+1)=572。即使符合度很好,这样的罚因子将使搜索算法从这个值移开而趋向15 mm/h。

表9-5 校验参数的约束条件

续表

参考文献

[1] Diskin,M.H.and Simon,E.(1977).“A procedure for the selection of objective functions for hydrologic simulation models.”Journal of Hydrology,34,129-149.

[2] Stephenson,D.(1979).“Direct optimization of Muskingum routing coefficients.”Journal of Hydrology,41,161-165.

[3] US Army Corps of Engineers,USACE (1998).HEC-1 flood hydrograph package user's manual.Hydrologic Engineering Center,Davis,CA.

[4] USACE (2000).HEC-HMS hydrologic modeling system user's manual.Hydrologic Engineering Center,Davis,CA.