一种基于实测径流数据的水文预报状态变量实时修正方法

摘要

本发明公开了一种基于实测径流数据的水文预报状态变量实时修正方法，该方法步骤如下：步骤1，收集和整理水文气象实测数据和数值降水预报数据；步骤2，根据步骤1中的水文气象实测数据构建流域水文模型，进行模型参数优化，确定模型最优参数组；步骤3，选定影响模型产汇流的重要状态变量，根据模型产汇流原理和步骤2中的最优参数组反推状态变量的计算公式，并依据实时观测径流计算当前时刻的状态变量；步骤4，利用数据同化方法将步骤3计算得到的状态变量作为模型初始条件输入水文模型，引入数值降水预报数据作为模型驱动进行洪水预报。该方法可提高流域水文模型对洪水的模拟和预报精度，有效延长洪水预见期。

权利要求

1.一种基于实测径流数据的水文预报状态变量实时修正方法，其特征在于，包括如下步骤：步骤1，收集和整理水文气象实测数据和数值降水预报数据；

步骤2，根据步骤1中的水文气象实测数据构建流域水文模型，进行模型参数优化，确定模型最优参数组；

步骤3，选定影响模型产汇流的重要状态变量，根据模型产汇流原理和步骤2中的最优参数组反推状态变量的计算公式，并依据实时观测径流计算当前时刻的状态变量；状态变量的反推方法为自下而上法，即根据径流数据和产汇流公式进行溯源反推；

步骤4，利用数据同化方法将步骤3计算得到的状态变量作为模型初始条件输入流域水文模型，引入步骤1中的数值降水预报数据作为模型驱动进行洪水预报。

2.如权利要求1所述的基于实测径流数据的水文预报状态变量实时修正方法，其特征在于，步骤1中所需要的水文气象实测数据包括逐日降水数据、蒸散发数据和流域出口径流数据；所述数值降水预报数据来自TIGGE数据集，包括欧洲中值天气预报中心ECMWF、美国环境预报中心NCEP、中国气象局CMA或加拿大气象中心CMC预报中心发布的数值降水预报数据。

3.如权利要求1所述的基于实测径流数据的水文预报状态变量实时修正方法，其特征在于，步骤1中所涉及的水文气象实测数据为流域面平均数据，经由泰森多边形方法计算。

4.如权利要求1所述的基于实测径流数据的水文预报状态变量实时修正方法，其特征在于，所述步骤2中的流域水文模型采用CR4J、HBV或新安江模型。

5.如权利要求1所述的基于实测径流数据的水文预报状态变量实时修正方法，其特征在于，步骤2中用于模型参数优化的方法为NSGA‑II、εNSGA‑II或SCE‑UA全局优化算法。

6.如权利要求1所述的基于实测径流数据的水文预报状态变量实时修正方法，其特征在于，步骤3中影响模型产汇流的重要状态变量包括模型中的土壤水储量。

7.如权利要求1所述的基于实测径流数据的水文预报状态变量实时修正方法，其特征在于，步骤4中的数据同化方法为直接同化法。

说明书

一种基于实测径流数据的水文预报状态变量实时修正方法

技术领域

[0001] 本发明涉及一种基于实测径流数据的水文预报状态变量实时修正方法，属于洪水预报领域，主要用于水利部门开展洪水预报预警工作。

背景技术

[0002] 由于城市化和气候变化的影响，洪水干旱等自然灾害的发生越来越频繁，严重威胁着国民经济和人民生命财产安全，利用实测径流数据实时修正水文模型状态变量，提高模型对洪水的模拟和预报能力，对缓解灾情有着重要意义。

[0003] 在实际流域中，由于诸多因素的影响，洪水预报存在较大的误差。其中，影响洪水预报精度的主要因素包括输入降水预报的误差、水文模型的误差、观测资料误差等，因此利用实时观测数据及时修正模型误差，可以显著弥补洪水预报的不足。目前大部分研究集中于对土壤含水量、降水预报以及输出流量的实时修正，考虑到水文预报是根据水文模型经由一系列的中间过程计算得到的，状态变量作为计算过程的中间值，其取值对模型计算有着直接的影响，初值不同，计算的结果相差很大。因此，对状态变量进行实时修正，可以起到更为直接有效的效果。

[0004] 一般来说，水文模型尤其是集总概念式水文模型的状态变量一般均具有一定的物理意义，原则上状态变量值都是可以测量的，但是在实际应用中往往难以直接测量，主要原因包括流域的不均匀性、水文模型状态变量的多样性等。目前状态变量的确定方式有两种方法，一种是在年初人为给定状态变量初值，另一种则是根据人工经验给定计算开始时刻的状态变量初值。但这两种方法均存在依靠经验、过度主观的缺陷，造成状态变量的给定存在较大误差。初值误差在模型系统内不断传输、繁衍，最终导致模型输出结果与实际相差甚远。

[0005] 现今国内外洪水预报的实时修正方法有很多，其中基于自回归的AR模型和卡尔曼滤波模型最为常用。AR模型结构简单、计算简便，但其修正效果依赖于前后时段径流误差的相关性，并且往往对洪峰的修正效果有限。卡尔曼滤波较为复杂，可容纳的信息量大，但实时洪水常常难以提供修正所需要的足够信息，大大限制了该方法的使用，导致其修正效果与AR模型效果类似。

发明内容

[0006] 为解决上述现有技术的不足，本发明的目的在于提供一种基于实测径流数据的水文预报状态变量实时修正方法，以实现洪水实时修正，提高洪水预报的精度并延长预见期。

[0007] 为实现上述目标，本发明通过如下技术方案实现：

[0008] 一种基于实测径流数据的水文预报状态变量实时修正方法，包括如下步骤：

[0009] 步骤1，收集和整理水文气象实测数据和数值降水预报数据；

[0010] 步骤2，根据步骤1中的水文实测气象数据构建流域水文模型，进行模型参数优化，确定模型最优参数组；

[0011] 步骤3，选定影响模型产汇流的重要状态变量，根据模型产汇流原理和步骤2中的最优参数组反推状态变量的计算公式，并依据实时观测径流计算当前时刻的状态变量；

[0012] 步骤4，利用数据同化方法将步骤3计算得到的状态变量作为模型初始条件输入水文模型，引入数值降水预报数据作为模型驱动进行洪水预报。

[0013] 上述技术方案中，进一步地，所述步骤1中的水文气象实测数据主要包括逐日降水数据、蒸发数据和流域出口径流数据等。所述集合降水预报数据来自TIGGE(THORPEX Interactive Grand Global Ensemble，全球交互式超级大集合)数据集，包括但不限于欧洲中值天气预报中心ECMWF(European Centre for Medium‑Range Weather Forecasts)、美国环境预报中心NCEP(National Centers for Environmental Prediction)、中国气象局CMA(China Meteorological Administration)、加拿大气象中心CMC(Canadian Meteorological Centre)等预报中心发布的中短期数值降水预报数据。

[0014] 进一步地，所述步骤1中的水文气象实测数据为流域面平均数据，经由泰森多边形方法计算。泰森多边形计算面雨量的包括以下步骤：

[0015] (1)在ArcGIS软件中输入流域轮廓图、雨量站(P1,P2,…Pi…,PN)，将流域划分为N个泰森多边形；

[0016] (2)计算第i个雨量站Pi所在的泰森多边形的面积，计算其与流域总面积的比值wi，作为雨量站Pi对流域面雨量的贡献值；

[0017] (3)流域面雨量的计算公式为：

[0018] 进一步地，所述步骤2中的流域水文模型根据实际情况进行选择，但主要考虑集总概念性水文模型，包括但不限于CR4J、HBV和新安江模型等。该类模型产汇流计算公式简单，更容易由实测径流数据反推状态变量。

[0019] 进一步地，所述步骤2中用于模型参数优化的方法为全局优化算法，包括但不限于NSGA‑II、εNSGA‑II、SCE‑UA等。

[0020] 进一步地，所述步骤3中影响模型产汇流的重要状态变量主要包括模型中的土壤水储量，通过文献阅读确定。包括但不限于GR4J模型中的产流水储量S和汇流水储量R、HBV模型中的上层水储量SW、下层水储量GW等。

[0021] 进一步地，所述步骤3中状态变量的反推方法为经验自下而上法，即根据径流数据和产汇流公式进行溯源反推。

[0022] 进一步地，所述步骤4中的数据同化方法为直接同化法，将反推状态变量直接输入模型作为模型当前时刻的状态变量取值。

[0023] 进一步地，所述步骤4中的洪水预报为逐日洪水滚动预报。

[0024] 通过采用上述技术手段，本发明的有益效果为：

[0025] (1)本发明方法提出的基于实测径流的状态变量反推方法能够解决当前状态变量修正主要基于经验以及状态变量难以实测获取的问题，本发明方法简便，对数据需求小，效果显著，是解决水文预报初始误差的重要手段。

[0026] (2)本发明直接利用径流数据实现洪水预报状态变量的实时修正，能够显著减小水文模型的初始误差，大幅度提高模型对洪水的模拟和预报精度。本发明方法简单有效，具有较高的实用价值，在业务预报中具有较高的应用前景。

附图说明

[0027] 图1为本发明的一种基于实测径流数据的水文预报状态变量实时修正方法整体实现流程图。

[0028] 图2是本发明的HBV模型的原理图。

[0029] 图3是本发明的实时修正方法对模型模拟效果的影响。

[0030] 图4是经本发明实时修正的洪水预报与未修正预报的CRPS指标图。

[0031] 图5是未实时修正的洪水预报在预见期为1天、2天、3天的表现。

[0032] 图6是经本发明实时修正的洪水预报在预见期为1天、2天、3天的表现。

具体实施方式

[0033] 下面结合附图和具体实施方式对本发明做进一步说明。采用HBV模型作为案例进行具体实施，包括如下步骤：

[0034] 步骤1，流域地理信息、水文气象实测数据准备

[0035] 使用ArcGIS对目标流域的地理信息数据进行提取和处理，获取流域轮廓以及流域集水区面积。

[0036] 流域潜在蒸发根据彭曼公式基于日最高温、日最低温、风速、相对湿度以及日照时长等数据计算，具体计算过程如下：

[0037]

[0038] 其中PET为潜在蒸散发(mm/day)；Δ是饱和水汽压‑温度曲线斜率(kPa/℃)；Rn和G2
分别是净辐射和地表热通量(MJ/m .day)；γ是湿度计常数(kPa/℃)；T、u2、es和ea分别为距地面2米处的平均温度(℃)、平均风速(m/s)、饱和水汽压(kPa)和实际水汽压(kPa)。

[0039] 面雨量根据泰森多边形法计算流域平均，包括以下步骤：

[0040] 1)在ArcGIS软件中输入流域轮廓图、雨量站(P1,P2,…Pi…,PN)，将流域划分为N个泰森多边形；

[0041] 2)计算第i个雨量站Pi所在的泰森多边形的面积，计算其与流域总面积的比值wi，作为雨量站Pi对流域面雨量的贡献值；

[0042] 3)流域面雨量的计算公式为：

[0043] 步骤2，水文模型参数优化

[0044] 根据以下目标函数NSE(Nash‑Sutcliff Efficiency，确定性系数)利用全局优化算法SCE‑UA来率定水文模型：

[0045]

[0046] 式中N为日径流数量；Qobs和Qsim分别为实测日流量和模拟日流量； 是实测径流的平均值。

[0047] HBV模型待优化参数及其范围如表1‑1所示。

[0048] 表1‑1HBV待优化参数范围

[0049]

[0050]

[0051] 复合形交叉进化算法SCE‑UA是在单纯形法的基础上融合自然界生物竞争进化理论和基因算法基本原理等概念综合而成的一种全局优化算法，在解决全局优化问题方面更加有效、灵活、稳健，被认为是非线性复杂模型参数优选最有效的方法之一。算法步骤如下：

[0052] 1)初始化。假设待优化问题为n维问题(n为待优化的参数数目)，参与进化的复合型个数为p，每个复合型所包含的顶点数为m＝2n+1，计算样本点数目为s＝pm。p的取值需要单独考虑，参考文献选取p＝5；

[0053] 2)产生样本点。在可行空间内(参数上下限之间)随机生成s个样本点记为x1,x2…xs，计算每个样本点对应的目标函数值记为f1,f2…fs；

[0054] 3)样本点群排序。将得到的s个样本点的目标函数值按升序排列，记为fi，与之对应的xi记为(xi，fi)，并储存到数组D中；

[0055] 4)进行复合形计算。将数组D分成p个复合形记为A1,A2…Ap，每个复合形包含m个样本点。依据竞争的复合形演化算法CCE分别对各复合形进行演化。然后将A1,A2…Ap代替到D中，并对D按目标函数的升序进行排序；

[0056] 5)收敛性判断。若达到收敛条件(本案例选取条件为：目标函数值的改进程度已经不能达到0.01％)则停止，否则回到第4)步。

[0057] 步骤3，根据实测径流反推模型状态变量

[0058] 基于文献阅读以及HBV模型的产汇流原理(图2)，不难看出HBV模型存在两个重要的中间状态变量，即上层土壤储水量SW，下层土壤储水量GW。模型产汇流公式如下：

[0059] Qf(t)＝Kf*SW(t)1+∝

[0060] Qs(t)＝Ks*GW(t)

[0061]

[0062] 其中，Qf和Qs分别代表地表径流和地下径流，Area为流域面积。

[0063] 当已知某一时刻的实测径流Qobs(t0),可以根据当前时刻模型的模拟地表径流Qf,sim(t)与模拟总径流Qsim(t)的比值推估当前时刻的状态变量SW和GW：

[0064]

[0065]

[0066]

[0067]

[0068] 步骤4，洪水预报与评估

[0069] 将SW(t0)和GW(t0)代入HBV并作为t0时刻模型的初始条件，输入气象数据进行洪水模拟和预报。洪水的评估指标采用CRPS(Continuous Ranked Probability Score)，该指标量化了集合预报连续概率分布和实测值之间的一致程度，是最为常用的集合洪水预报评估指标。CRPS越小表示预报误差越好，其计算公式如下：

[0070]

[0071] 式中，N是样本数量，x为径流预报，xa是实测径流；F是径流的分布函数，H是著名的Heaviside函数：

[0072]

[0073] 图3表明经本发明实时修正后HBV模型的径流模拟能力明显提升，模拟径流与实测径流的水文过程线符合得更好，变化完全一致。

[0074] 图4显示本发明实时修正径流预报与实测径流计算得到的CRPS指标比未实时修正预报径流的CRPS低，说明本发明能够显著降低预报误差，提高预报精度。特别是预见期为1～3天时，本发明实时修正预报的CRPS不到未修正预报CRPS的50％，对径流预报的改善效果显著。

[0075] 图5和图6对比可知，在预见期为1天时，本发明可以将洪水预报有效提升至与实测径流基本一致的水平；预见期为2～3天时，本发明能够显著提升涨水段和退水段的洪水预报。

[0076] 以上所述仅是本发明的实例实施而已，并不用于限制本发明，对于本领域的研究者来说，本发明可以有各种更改和变化。凡是在本发明的权利要求限定范围内，所做的任何修改、等同替换、改进等，均应在本发明的保护范围之内。