集合卡尔曼滤波在河道洪水预报中的应用.doc

基于集合Kalman滤波的河道洪水预报研究岳延兵1,2，李致家1，李振兴2 (1.河海大学水文水资源学院、江苏 南京 210098；2.山西水利职业技术学院、山西 运城 044004)摘 要：基于集合Kalman滤波的数据同化技术应用逐渐成熟，研究了集合Kalman滤波法与最小二乘法对马斯京根模型参数率定的河道洪水预报技术，利用长水白马寺实测洪水进行了对比检验，讨论了集合Kalman滤波的多时段预报洪水过程及其特点。试验结果表明，应用集合Kalman滤波技术优于最小二乘法的预报效果，可有效延长洪水预报的时段和提高预报精度。关键词：集合Kalman滤波；数据同化；信息融合；水文预测1河道洪水预报研究简述河道洪水预报，是依据洪水波在河道中的运动、传播规律为基础，由河道上断面的水位、流量值预报下游断面水位、流量值的方法1-2。常见的河道洪水预报方法有两类：相应水位（流量）法和洪水演算法。按照描述洪水传播规律形式的不同，洪水演算法可分为水文学方法、水力学方法和基于系统理论的方法。基于系统理论的方法直接描述河道系统输入与输出之间的因果关系，不依赖于河道地形资料。基于系统理论的河道洪水演算方法具体可分为线性系统模型、时间序列模型、神经网络模型等2。目前利用小波分析、遗传算法、神经网络与混沌分析等进行信息融合的预测方法也逐渐发展起来3。2集合Kalman滤波同化技术Kalman滤波不假设背景场的误差是定常的，而是利用协方差矩阵预报方程来计算4。其原理是先进的，但是在实际应用中却因为计算量和存储量太大，在可预见的将来还无法应用到像最优插值方法现在应用到的那种规模的数值模式和观测系统。Kalman光滑是Kalman滤波的一种完善，它可以利用后一时刻的资料来改进前一时刻的估计。Kalman光滑分为固定区间光滑和固定时滞光滑两种，其计算量比Kalman滤波更大。目前都在研究Kalman滤波的简化形式，比如基于粗空间分辨率的模式来构造协方差矩阵预报方程。集合Kalman滤波(Evensen, 1994)是针对Kalman滤波中的协方差矩阵预报模式有时会出现计算不稳定等问题而提出的，其主要思想是抛弃协方差矩阵预报模式，而是直接利用Monte-Carol方法多次积分的模式，从而得到背景场的误差协方差矩阵，这个方法的计算存储量比Kalman滤波小很多5。集合Kalman(EnKF)是一种纯的统计Monte-Carlo罗方法，模式状态的集合在状态空间中随时间演变，集合的均值为模式预报值的最优估计，集合的发散分布代表了误差方差。在测量时刻，每一个测量值被另一个集合所取代，集合的期望值就是测量值的最优估计，集合的方差反映了测量误差。EnKF同化方法包含预报部分和分析部分6-8。预报部分：对于一个水文计算模型，模式将分析阶段得到的分析集合作为初始场向前预报至下一个观测资料时刻（即数据同化时刻t）： （1）式中：X为状态变量；M为预报模式。分析部分：在预报部分得到模型预测集合，用集合平均代替真实值，同化观测资料求得分析集合平均解。河道洪水流量是一个随机动态过程，测量流量值是含误差的不定值。集合Kalman滤波参数估计方法通过对观测值施加观测噪声，可以模拟现场观测值的不确定性。对于集合Kalman滤波而言，加入一定的噪声还可以提高滤波器的稳定性，保证估计量的标准差不会急剧下降6。采用Monte Carlo模拟高斯白噪声（均值为0，标准差为观测值乘以噪声比例系数），流量观测向量变为一组经过扰动的观测集合，其第i个样本可表示为： （2）式中：H为观测算子；e为观测误差。根据实测和预测值，结合预测误差矩阵，就可以直接更新模型参数： （3） （4） （5） （6）式中：K为卡尔曼增益矩阵；O为观测误差协方差矩阵；B为背景场误差协方差矩阵。集合Kalman滤波进行数据同化的计算流程如图1所示。预报部分变量的集合预报 分析部分1 预报误差的分析 2 计算增益矩阵 3 对变量的分析 图1集合Kalman滤波的计算流程3 最小二乘法率定马斯京根模型参数与检验为了便于对比研究，首先利用最小二乘法率定马斯京根模型参数并检验预报效果。马斯京根模型的槽蓄方程式和流量演算方程式如下9： （7） （8）式中：I为上游输入流量；O为下游输出流量；K为蓄量流量关系曲线的坡度；X为流量比重系数；C为流量演算方程系数。利用长水白马寺实测洪水数据，应用最小二乘法对马斯京根模型参数率定结果：X=0.35，K=11.86 ，rt=12。计算得：C0=0.134875，C1=0.740462，C2=0.124663。对白马寺站进行了1975年-1989年共4场模拟洪水预报，结果见表1、图2所示。表1 白马寺4场检验洪水的模拟洪水预报特征值洪号开始时间结束时间预报洪峰流量(m3)实测洪峰流量(m3)峰值绝对误差(m3)平均绝对误差(m3)确定性系数19750805141975-8-5 14:001975-8-14 8:001539192038176.20.932719810906001981-9-6 0:001981-9-14 10:001034111076360.955019831004131983-10-4 13:001983-10-13 7:00189218801261.810.974719890709101989-7-9 10:001989-7-17 23:00122413209637.350.9772 图2 1975年1989年共4场检验洪水实测与预报流量过程对比4 集合Kalman滤波同化洪水预报模型与检验对白马寺站进行了1975年-1989年共4场逐时段预见期的模拟洪水预报，应用集合Kalman滤波模型预报结果见表2 、图3所示。表2 白马寺4场检验洪水的逐时段模拟洪水预报特征值洪号开始时间结束时间预报洪峰流量(m3)实测洪峰流量(m3)峰值绝对误差(m3)平均绝对误差(m3)确定性系数19750805141975-8-5 14:001975-8-14 8:0019771920577.540.997119810906001981-9-6 0:001981-9-14 10:001117111073.730.998219831004131983-10-4 13:001983-10-13 7:0019131880338.410.998319890709101989-7-9 10:001989-7-17 23:0013581320384.950.9985图3 1975年-1989年共4场检验洪水实测与预报流量过程对比计算结果显示，应用集合Kalman模型预报结果优于最小二乘法的预报结果，试验效果较好。为了进一步研究延长预见期后的情况，下面应用集合Kalman滤模型在逐时段rt基础上进行2rt时段模拟洪水预报。计算预报结果见表3 、图4所示。表3 白马寺4场检验洪水的2rt时段模拟洪水预报特征值洪号开始时间结束时间预报洪峰流量(m3)实测洪峰流量(m3)峰值绝对误差(m3)平均绝对误差(m3)确定性系数19750805141975-8-5 14:001975-8-14 8:002093192017318.110.988819810906001981-9-6 0:001981-9-14 10:0011511110419.740.991619831004131983-10-4 13:001983-10-13 7:00191918803919.980.993719890709101989-7-9 10:001989-7-17 23:00141313209315.500.9889图4 1975年1989年共4场检验洪水实测与rt、2rt时段预报流量对比计算结果显示，延长预见期后，应用集合Kalman滤波模型预报结果仍优于应用最小二乘法的预报结果。集合Kalman滤波模型预报结果之间对比：预见期2rt时段比rt时段流量预报的峰值绝对误差略大，确定性系数降低，但都满足要求，总体效果较好。预见期2rt时段模拟洪水过程比rt时段模拟洪水过程更敏感，这是误差放大的结果。为了进一步研究更长预见期多倍rt时段后的预报效果，下面计算了试验的4场检验洪水预见期分别为18rt时段预报特征值并绘制了关系曲线。表4 预见期18rt时段的白马寺4场检验洪水预报特征值 N rt洪号123456781975080514峰值绝对误差(m3)5717330910361285-94-10-823平均绝对误差(m3)7.5418.1142.1786.07132.8155.9197.3228.7确定性系数0.99710.98880.95810.82400.65910.53180.36710.24831981090600峰值绝对误差(m3)741141282559-13272142平均绝对误差(m3)3.739.7420.4640.0264.4483.5119138.7确定性系数0.99820.99160.97320.92390.82860.79840.62830.41821983100413峰值绝对误差(m3)333963372253599102522平均绝对误差(m3)8.4119.9841.8476.27112.6150.8189.4226.8确定性系数0.99830.99370.97790.93410.87980.81870.75190.66051989070910峰值绝对误差(m3)3893265667412-10079-110平均绝对误差(m3)4.9515.5030.2252.5480.0792.03115.9129.1确定性系数0.99850.98890.95930.90180.82920.80400.75680.7225图5 1975年1989年共4场检验洪水的预见期与预报特征值关系计算结果显示，随着预见期的多倍延长，洪水流量预报的峰值绝对误差逐渐放大并在5rt时段后出现波动，峰值平均绝对误差也逐渐放大，确定性系数逐渐降低。根据水文资料情况，rt的单位可分别选择小时、日、周、月等，以满足不同水文预报的需要。集合Kalman滤波数据同化技术在水文预报领域的应用还需进一步研究。5 结 论集合Kalman滤波数据同化方案可以融合来自多源的信息，试验显示较最小二乘率定马斯京根参数法预报效果好。集合Kalman滤波数据同化预报洪水过程随rt时段成倍增长而敏感，且预报的峰值绝对误差逐渐偏大，在确定性系数满意的前提下预报时段可有效延长。本文仅利用了水文观测数据本身进行同化试验来验证算法，今后可以利用卫星、雷达等各种间接观测信息进行信息融合，不断修正模型的预测轨迹，提高模型估计的精度。目前建立的数据同化方案多以研究单一变量为主，并开展理想和单点试验作为验证，数据同化算法的研究和发展尚需不断完善和测试。参考文献：1 Anderson M G, Burt T P. Hydrological ForecastingM. New York, John Wiley & Sons, 1985.2 芮孝芳. 洪水预报理论的新进展及现行方法的适用性J. 水利水电科技进展. 2001, 21(5): 14. Rui Xiao-fang . Advances in Theories of Flood Prediction and Suitability of Existing Methods j. Advances In Science and Technology of Water Resources. 2001, 21(5): 14.3 岳延兵、李致家、范敏.基于信息融合的水文预测技术研究J.水资源与水工程学报.2009,10(5):91-95. Yue Yan-bing, Li Zhi-jia, Fan Min. Hydrological Prediction Technology Research Based on the Information Integrationj. Journal of Water Resources and Water Engineering.2009,10(5):91-95.4 高拴柱 . 集合Kalman滤波资料同化技术及研究现状J .气象. 2005（6）38. Gao Shuanzhu . Review on Ensemble Kalman Filter Data Assimilation j. Meteorological Monthly. 2005（6）38.5 Jurgen Komma,Gunter Bloschl,Christian Reszler . Soil moisture updating by Ensemble Kalman Filtering in real-time flood forecastingJ. Journal of Hydrology. 2008（357）,228242.6Burgers G, van Leeuwen PJ, Evensen G. Analysis scheme in the ensemble Kalman FilterJ. Mon Weather Rev,1998（126）1719-1724.7Reichle RH, McLaughlin D, Entekhabi D., Hydrologic data assimilation with the ensemble Kalman filterJ.Monthly Weather Rev, 2002（130）103-14.8Evensen G. The Ensemble Kalman Filter:Theoretical formulation and practical implementationJ. Ocean Dynamics, 2003（53）343-367.9 包为民. 水文预报M. 北京：中国水利水电出版社.2006. Bao Weimin. Hydrological forecasting M.Beijing: China Water Power Press.2006Study of River Channel Flood Forecasting by Ensemble Kalman FilterYue Yanbing 1，2， LI Zhijia 1 ，li Zhengxi 2(1. College of Hydrology and Water Resources, Hohai University, Nanjing 210098;2. Shanxi Conservancy Technical College, Yuncheng 044004 ) Abstract: Techno