水动力学模型实时校正方法对比

刘开磊，姚 成，李致家，阚光远，包红军

(1.河海大学水文水资源学院，江苏 南京 210098;2.中国气象局公共气象服务中心，北京 100081)

水动力学模型模拟洪水演进具有比较成熟的理论，同时又是存在很多不确定性因素的领域，受降雨分布、下垫面、人为调蓄等影响，水流运动情况相当复杂，使洪水模拟精度难以从改进模型结构方面进行提高。应用实时校正方法实现模型状态与预报结果的实时校正，可以在保证预见期的前提下达到较高的预报精度。依据校正对象进行分类，与水动力学模型联合应用的实时校正方法大致分为以下3类:

第一类是对状态变量的校正。理论上讲，状态调整实时校正最完善的是Kalman滤波理论，其基于模型和模型输入变量的先验知识 通过估计模型先验状态 先验误差协方差与残余的增益 对模型状态进行后验估计。将河道汇流关系用状态方程、量测方程进行描述，给出噪声估计特征即可实现递归计算。譬如:王井泉等［1］尝试将水文水力学方法耦合以作流量与水位预报，并将模型线性化与Kalman滤波方法耦合，实现半自适应的Kalman滤波用于水位校正;葛守西等［2］选择不同的状态变量构建Kalman滤波实时校正模型;王船海等［3］提出了“多步交替校正”的方法，分别建立水位与流量的状态空间方程，进行交替滤波计算。

第二类方法是对预报结果的校正。将河道汇流关系以确定性模型描述，累积预报值与相应时刻实测值之间的误差，再以系统识别方法对误差序列建立模型和估算参数［4］。作为对确定性模型的补充描述，这类方法就是模拟加校正模型，是对模型结果的校正，其中具有代表性的方法是基于具有遗忘因子的最小二乘法的误差自回归实时校正。传统误差自回归方法在水动力学模型模拟中已获得普遍应用，该方法的最大特点是简单，并且不需要提前假定误差的来源［5］。另外采用基于KNN［6］的非参数回归分析实时校正，该方法不需要建立自回归方程，不需要计算相关系数，与传统自回归方法有较大区别。

第三类方法是以上两种校正方法的综合，既包含模型运算过程中对状态变量的更新，也包含对模拟结果的实时校正。其基本思想在于实现在模型状态变量校正之后计算出相应的水位流量校正值，然后应用Kalman滤波或者误差自回归校正法等对模拟结果进行再校正，本质上是前两类方法的结合，以下简称综合法。

文本采用代表性的实时校正方法:基于误差自回归的实时校正(以下统称误差自回归校正法)、基于Kalman滤波的实时校正(以下称Kalman滤波校正法)、基于KNN的实时校正(以下称KNN校正法)，以及Kalman滤波校正与KNN校正法结合的综合法，分别与水动力学模型联合应用，在试验流域建立基于水动力学模型的洪水预报及实时校正系统，以对比各种校正方法的稳定性能及实时校正能力，并对校正结果进行比较分析，为实时校正方法的应用提供参考。

1 试验河段选择

试验河段(图1)选用淮河吴家渡—小柳巷段，主河道全长152 km，区间有临淮关水位站，为方便参数调节，河段内从已有的540个实测断面中选择12个断面，平均间隔14 km左右，经验证发现，根据所选择断面与540个断面所建立的水动力学模型模拟结果相差不大，12个断面满足模型的模拟精度要求。

河段洪水均以上游吴家渡来水为主，区间无降雨洪水，也无旁侧支流加入，行蓄洪调度行为较少，适于建立无旁侧入流影响的水动力学模型，便于比较各种实时校正方法的性能。

采用Preissmann四点隐式差分格式对一维圣维南方程组进行离散化，建立水动力学模型，以上游吴家渡流量过程作为上边界(非实测期内采用上游水文学模型预报的流量过程作为上边界实现模拟预报)，下游小柳巷水位流量的多元回归关系作为下边界条件。区间的临淮关站点具有较长年限的水位实测资料，可实现对整个试验河段的一元三点不等距插值，其中临淮关站点水位(无插补误差)可以体现Kalman滤波校正模拟效果。

图1 试验河段示意图Fig.1 Sketch map of experimental river reach

2 实时校正理论与模型构建

2.1 误差自回归校正理论

基于误差自回归的实时校正理论认为误差序列存在序贯相关性，通过建立误差自回归方程可以较准确地估计当前预报误差，适用于任何预报方法［7］。

设{y(t)}是流量误差时间序列，建立自回归模型如下:

式中:ext——预见期;N——自回归的阶数;α0，α1，…，αN——系数;e(t)——方程的余项。

由于在洪水演进过程中流量与水位之间具有相关转换的关系，以流量误差为例，研究认为流量误差不仅与历史流量预报误差相关，也与相应时刻的水位预报误差z(t)相关，建立流量误差自回归模型如下:

在实际模拟过程中首先需要确定模型的阶数，并选用合适的方法进行系数估计。研究采用损失函数法［8］进行阶数确定，建立小柳巷站点流量误差平方和与N的相关曲线定阶，模拟流量通过基于多元回归分析的水位～流量关系转换为水位模拟值。

2.2 KNN校正法原理

KNN非参数估计技术是近几年来在气象数值预报使用中颇为成功的一种方法，它仅仅依靠已经积累的、包含系统潜在关系的大量数据对目标状态进行估计，不需要建立相关的模型，不需要对模拟过程的先验知识［6］。洪水演进过程中相近的河道状况及相似的天气条件下往往会产生相似的洪水过程和模拟误差，例如受变动回水、洪水涨落、溃口等因素影响，低水状态的水位～流量关系难以精确模拟，导致水动力学模型在低水状态的模拟预报精度下降;在模型模拟中无法精确考虑到区间产汇流过程，因此大洪水时的模拟洪峰比实际洪峰要小一些，模拟的峰后过程往往比实际过程更陡。类似现象的客观存在，使KNN校正法有可能在洪水预报实时校正中发挥作用［9］，该方法可直接根据训练样本数据建立基于KNN非参数的多输入、单输出模型，避免最小二乘法计算量较大、洪峰数据拟合较差的问题。此外，还可以避免Kalman滤波校正法需要将模型线性化，校正模型建立不易，且易受状态向量、观测向量选择的影响，受噪声协方差矩阵影响较大，校正效果不稳定问题。

KNN校正法是一种统计学习方法，基于KNN的实时校正方法由4个部分组成，第一部分用于存放全部历史误差数据，一般取模拟时段的前70%用于以后时段内水文要素的实时校正，第二部分是近邻样本数据对的向量化，本系统中基于误差自回归的方法，认为当前时刻的预报误差与前面多个时段的(NT)的预报误差成某种相关关系，pi+ext=f(pi-θ，…，pi-2，pi-1)。ΔT代表计算步长，本研究中取ΔT=1h。设定预见期为ext，需要构建m个误差向量，则可以从N+ext时刻开始至(m-1)+N+ext建立m个历史预报误差向量集。第三部分最近邻样本识别是近邻误差向量样本的对比与搜索算法的实现，本系统采用欧氏距离作为评价向量间相似相关性的指标，以快速排序算法找出与当前预报误差相关的向量最相似的K个预报误差向量。第四步误差估计，应用反距离权重法对K个预报误差向量所对应的预报误差加权作为对当前误差的估计(见图2)。

图2 KNN实时校正方法原理示意图Fig.2 SchematicdiagramofprincipleofKNNreal-timecorrectionmethod

KNN校正法也是一种仿真预报员根据河道汇流(或降雨-径流等)原理与经验进行误差数据分析预报的方法，从实现过程看，这种方法是集数值预报、误差自回归、统计校正方法为一体的综合校正技术。KNN校正法原理清晰 建模计算简便 试验证明其对洪水预报校正能力较强 并且能在较长预见期内保证较高的精度2.3 Kalman滤波校正理论与模型构建

应用Kalman滤波校正法实现水动力学模型实时校正需要处理几个关键的问题［7］:(a)水动力模型的线性化;(b)状态向量、量测向量的选择;(c)误差协方差阵的结构与参数设定。其中如何设定误差协方差阵的参数是至今仍未解决的问题，本文采用Sage-Husa自适应滤波方法实时更新Kalman滤波校正法中模型噪声协方差阵的思路，经淮河流域2003—2009年多场洪水验证表明该方法较稳定，模拟效果较好，可以为解决协方差阵参数特征估计问题提供一定的借鉴。

2.3.1 Kalman滤波校正法与水动力学模型耦合的实时校正模型构建

描述水动力学方程的圣维南方程组包括连续方程和动量方程，离散化之后的方程组对于i河段的连续方程和动量方程则分别简化为［10］

式中:i——断面序列号，取值1 ～N;j——时间序列号，取值 1 ～t;C、D、E、F、G——对应断面的时段参数。在建模过程中只要给出初始时刻各断面的水位流量状况，就可以根据相应的参数公式由上一时刻的水流状态计算得出，因此只需已知河道初始流量及水位即可递推计算参数值。

由方程(3)(4)对N个断面的N-1个河段可写出2N-2个方程，加入上下边界条件(根据研究流域特点，上边界选用流量过程，下边界选用水位-流量关系)构成2N个线性方程，其中含有2N个未知量(N个断面的Zj和Qj)，该方程组有唯一解。根据式(3)(4)，方程组的矩阵形式为［3］

式中:H——系数矩阵，为时变矩阵，其元素值由相应的水动力学模型推算得出，但是对于当前时刻是已知的。

观察式(5)，设定状态向量 xt=［G0，D1，φ1，D2，φ2，…，DN-1，φN-1，GN］，观测向量 yt=［Q1，Z1，Q2，Z2，Q3，…，ZN-1，QN，ZN］，可以构建Kalman滤波方程，其基本方程状态空间方程及观测方程分别记为式(6)、式(7):

式中:ωt、νt——模型误差向量、观测误差向量。

与注塑部件生产和组装有关的工艺过程为最新一代能够处理图像视觉的传感器提供了广泛的应用空间。Sensopart传感器公司的Klaus Berdel指出：“通常注塑成型过程中最常见的问题是短射引致的模具不完全填充以及过度注射导致的毛刺和飞边。”利用预置好的工具对工件的表面和轮廓进行检测很容易发现部件的几何形状有这样或那样的偏差，此外视觉传感器还能够及时检测插入错误，例如在注塑问题的模制插头中发现未校准的金属销。

由以上方程建立的基于Kalman滤波的水动力学模型实时校正系统［7］可以实现对Kalman滤波方程状态向量X的实时校正，需要注意的是:(a)该状态向量区别于水动力学模型的状态(水位、流量)，校正后的X回代入式(3)(4)实现对模型状态的实时校正;(b)观测向量只能选择有实测资料的水位或者流量，但是由于河段站点较少，资料获取受限，在应用中采用一元三点不等距插值插补无实测断面的水位、流量，由于资料缺乏引起的观测误差体现在观测误差协方差矩阵R中。

综合法［11］本质上为两步实时校正的结合，第一步对模型状态校正，方法的选择上首推Kalman滤波校正法;第二步对模型结果实时校正，考虑到误差自回归校正法对洪水过程发生突变部位的跟踪校正能力较差，因此在方法上优先选择KNN校正法。两种方法前文已有叙述。

2.3.2 Kalman滤波模型参数设定

Kalman滤波模型中的Q和R矩阵是滤波器设计的关键参数，参数敏感性较强，对模型校正能力影响较大，然而其矩阵形式与矩阵元素取值相对难以确定。国内外对其参数形式与取值的讨论较多［12］，没有绝对适合的方法来确定噪声矩阵的形式与取值方法。

水动力学模型中有实测资料的站点远少于需要模拟的站点，因此理论上Q和R的特征无法全部估计出来，由于模型以及Q和R需要事先给出，而且Q和R设定不当还可能导致滤波发散，于是引入实时滤波技术，应用Sage-Husa自适应滤波方法实现对系统模型与观测噪声协方差矩阵的实时更新。

Sage-Husa自适应滤波方法自1969年提出并获得应用之后也获得很多的发展变化，然而其本质上仍然是对稳态系统的噪声估计，对于明显非稳态的典型水文过程其应用受到局限，试验研究发现Sage-Husa自适应滤波可以实现对水动力学洪水实时校正模型中q和Q的准确估计，加速校正模型的收敛，而对r和R则不能正确估计，甚至会导致校正模型发散，于是建立的Sage-Husa自适应滤波校正方法如下:

其中

式中:K——增益矩阵;z——对应时刻系统的观测量。

从式(8)可以看到当t很大时(t-1)/t趋近于1，1/t趋近于零，参数估计取决于上一时刻的参数，也就是参数估计逐步进入稳态，式(8)起不到实时校正的作用。

3 实时校正方法对比分析

模拟采用2003—2009年共5场洪水，以河道末端面(小柳巷)的实时校正结果作为指标进行多方法的模拟验证。Kalman滤波校正法在预见期内应用外推方法求得对应时刻的量测向量值。根据相关文献［13］，采用洪峰相对误差、确定性系数以体现各方法的模拟结果，模拟结果对比见表1。

表1 不同校正方法在不同预见期内的校正性能指标对比Table 1 Comparison of correction performance indices of different correction methods during different forecast periods

a.从以上4种校正方法校正结果的对比可以发现，除2003年、2007年2年中洪峰相对误差稍大，在其余的2005年、2006年、2009年洪水的8ΔT预见期内模拟预报结果，其确定性系数都能达到0.950以上，洪峰相对误差能保证在5.0%以下，说明文中所用到的误差自回归校正法、KNN校正法、Kalman滤波校正法以及综合法均能在实时洪水预报校正中取得良好的效果，足够为洪水预报调度与防洪决策提供理论数据支持。

b.2003年这场洪水的校正过程中误差自回归校正法在洪峰部分误差较大，从表1可见该方法洪峰相对误差较大，在洪峰部位的模拟校正结果跳变，说明误差自回归方法在洪峰或者洪水状况突变部分跟踪能力不够，从而导致对预报误差的预测偏离实际情况。究其原因，经分析后认为2003年、2007年汛期洪水量级较大，为实现错峰削洪，上游邱家湖、荆山湖等行蓄洪区使用频率增大，人为干预河道洪水演进过程的行为较多，尤其在洪峰前后行蓄洪区频繁使用，影响了河道洪水的运动状态，时间顺序的误差线性相关关系不再特别明显，影响了误差自回归方法在洪峰部位的校正效果。2003年、2007年各预见期内校正结果的变化也说明了这一点。

4 结 语

水力学往往被应用于传统水文学方法无法解决或解决不好的问题［1］，然而由于其模型本身的复杂性，限制了其预报精度的进一步提高［14］。本文对具有代表性的误差自回归校正法、KNN校正法、Kalman滤波校正法以及综合法进行讨论分析，得出以下结论:

a.对于类似2005年、2006年、2009年的中小型洪水过程，区间行蓄洪调度行为较少，4种方法均对预报结果的精度提升明显，能够较好地满足预报校正的需要。

b.对于类似2003年、2007年在上游或区间内部有较频繁行蓄洪调度行为的洪水过程，误差自回归实时校正方法对洪峰部位校正结果出现跳变，说明该方法已经不再可靠，应考虑采用其他3种校正方法。

c.基于Kalman滤波的实时校正方法可实现全断面水文要素的校正，相对于其他方法，其校正的范围更广，然而其噪声特征的估计是一个难点问题，较大程度地限制了该方法在水文领域的应用效果，文中提到的对各协方差矩阵分类进行估计的方法可以为噪声估计提供参考。

d.Kalman滤波校正与KNN校正法相结合的综合法能够较好地避免由于滤波参数估计不准确造成的误差放大，模拟校正的结果显示其模拟稳定性相对于Kalman滤波校正法以及KNN校正法都有提升，模拟精度较高，因此可以作为水动力学实时校正方法加以推广。

综上所述，在实际的应用中可以优先考虑综合法或KNN校正法、Kalman滤波校正法，并针对具体的河道优选校正模型。

［1］王井泉，李致家.卡尔曼半自适应滤波水位实时预报模型研究［J］.人民长江，2000，31(1):28-30.(WANG Jingquan，LI Zhijia.Semi-adaptive Kalman wave filter updating model for real time forecasting of water level［J］.Yangtze River，2000，31(1):28-30.(in Chinese))

［2］葛守西，程海云，李玉荣，等.水动力学模型卡尔曼滤波实时校正技术［J］.水利学报，2005，36(6):687-693.(GE Shouxi，CHENG Haiyun，LI Yurong，et al.Real time updating of hydrodynamic model by using Kalman filter［J］.Journal of Hydraulic Engineering，2005，36(6):687-693.(in Chinese))

［3］王船海，白耀玲.基于卡尔曼滤波技术的河道汇流实时校正［J］.河海大学学报:自然科学版，2007，35(2):181-185.(WANG Chuanhai，BAI Yaoling.Real-time correction of flow concentration of river channels based on Kalman filter technique［J］.Journal of Hohai University:Natural Sciences，2007，35(2):181-185.(in Chinese))

［4］ REFSGAARD J C.Validation and intercomparison of different updating procedures for real-time forecasting［J］.Nordic Hydrology，1997，28:65-84.

［5］杨敏.水文时间序列相似性模型的研究与应用［D］.南京:河海大学，2002.

［6］KARLSSON M，YAKOWITZ S.Nearest-neighbor methods for nonparametric rainfall-runoff forecasting［J］.Water Resource Research，1987，23(7):1300-1308.

［7］葛守西.现代洪水预报技术［M］.北京:中国水利水电出版社，1999:72-129.

［8］李致家，孔凡哲，王栋，等.现代水文模拟与预报技术［M］.南京:河海大学出版社，2010:125-127.

［9］DAVID R M.HandbooKof Hydrology［M］.New York:McGraw-Hill，1993:960-986.

［10］夏自强.多汊道串联感潮河段洪水演算方法探讨［D］.南京:河海大学，1988.

［11］CLOKE H L，PAPPERGER F.Ensemble flood forecasting:a review［J］.Journal of Hydrology，2009，375(3/4):613-626.

［12］HENRIKM，CLAUS S.Adaptive state updating in real-time river flow forecasting:a combined filtering and error forecasting procedure［J］.Journal of Hydrology，2005，308(1/4):302-312.

［13］VAN S N，RONSYN J，WILLEMS P.A non-parametric data-based approach for probabilistic flood forecasting in support of uncertainty communication［J］.Environmental Modeling＆ Software，2012，33:92-105.

［14］赖锡军，李纪人.洪水淹没范围数据与动力学模型的融合［J］.河海大学学报:自然科学版，2009，37(5):529-533.(LAI Xijun，LI Jiren.Assimilation of flood extent data into hydrodynamic model［J］.Journal of Hohai University:Natural Sciences，2009，37(5):529-533.(in Chinese))