CN112884223A - 基于多源卫星降水信息和产流约束校正的无资料地区洪水实时预报方法 - Google Patents

Abstract

本发明提供一种基于多源卫星降水信息和产流约束校正的无资料地区洪水实时预报方法，属于洪水实时预报技术领域。首先，拆分流域出口流量组成，获得校正对象。其次，采用多源卫星降水产品驱动水文模型后得到第一步获得的校正对象的多个序列，将多序列的最值形成的范围作为约束，进行产流约束校正和流量预报。最后，根据传统的评价指标评价每场洪水的预报表现，并进行参数率定。本发明通过将多源卫星降水信息转化为产流的上下限，从而对产流进行有约束的校正，可大幅提高卫星降雨产品驱动的水文模型的洪水预报精度，以支持无降雨资料地区洪水预报，并为多源卫星降雨产品的应用提供新思路。

Description

基于多源卫星降水信息和产流约束校正的无资料地区洪水实 时预报方法

技术领域

本发明属于洪水实时预报技术领域，涉及一种基于多源卫星降水信息和产流约束校正的无资料地区洪水实时预报方法。

背景技术

洪水实时预报可有效减轻洪水造成的损失。由于降雨是驱动水文模型的首要信息，因此精确的降雨信息在洪水实时预报中起着至关重要的作用。目前，基于地面站观测的降雨信息被广泛应用于洪水预报。但是，我国河流众多，受投资资金、地形气候以及设备维护等条件限制，目前仍有很多地区无降雨监测站点。近年来，随着卫星遥感技术的迅猛发展，定量降雨分析技术为解决无资料地区洪水预报提供了新的思路。目前，全球多机构通过融合多源卫星传感器发布了多种卫星降水产品，且每种产品包含多个版本，因此卫星降水产品版本繁多且各有优劣。

作为一种间接估测结果，卫星遥感估测降雨产品往往存在着系统误差和不确定性，其中近实时产品的误差尤为显著。即便如此，这仅有的近实时卫星降水产品对于无降雨资料地区是尤为可贵的，因为对于洪水实时预报而言，卫星降水产品的滞时长短直接影响洪水预报的预见期和时效性。大量卫星降水产品评价的研究表明，尽管近实时产品的精度明显低于后处理产品，但它们都能一定程度上描述降雨的特征和趋势。因此，亟需发明一项针对近实时卫星降水产品的洪水实时预报方法，从而充分利用近实时降水产品的有效信息。

许多研究在提高卫星降雨产品驱动的洪水预报精度方面做出尝试，其中最为常见的方法就是以卫星降雨产品为输入重率定水文模型，但由此得到的模型参数会较大程度地偏离参数的物理意义范围，导致模型参数不具备代表性，洪水预报结果存在较大不确定性。此外，进一步提高卫星降雨产品驱动的洪水预报精度的方法主要包含三类。第一类，在驱动水文模型之前，根据少量的或流域周边的地面站观测降雨与多源卫星降雨产品的误差统计关系，校正产品的误差，再以校正后的产品驱动水文模型进行洪水预报；第二类，利用多源卫星降雨产品分别直接驱动水文模型，再基于权重或贝叶斯平均等方法融合水文模型的模拟结果；第三类，在水文模型的模拟过程中利用实时观测的流量或土壤含水量更新校正洪水预报结果。多研究对比表明，由于流量可以被直接观测，利用实时流量校正预报结果在上述三类方法中表现最好。但利用观测流量直接校正预报流量大大缩短了洪水预报的预见期，因此利用观测流量间接校正产流或降雨成为了更优选择。此外，由于校正过程中缺乏约束，往往导致校正后的洪水过程跳动现象严重，洪水预报结果不稳定。但目前仍没有任何一种方法是基于多源卫星降水产品和产流约束校正来实时预报无资料地区的洪水。

发明内容

针对现有技术存在的问题，本发明提供一种基于多源卫星降水信息和产流约束校正的无资料地区洪水实时预报方法。

本发明采用的技术方案为：

一种基于多源卫星降水信息和产流约束校正的无资料地区洪水实时预报方法，包括以下步骤：

第一步，拆分流域出口流量组成，获得本发明的校正对象。

流域尺度的水文模型可以被表达为：

Q＝f(P，E，X，θ) (1)

其中，Q＝[Q₁，Q₂，Q₃，…，Q_m]^T、P＝[P₁，P₂，P₃，…，P_n]^T、E＝[E₁，E₂，E₃，…，E_m]^T和X＝[X₁，X₂，X₃，…，X_m]^T分别是流量、降雨、蒸散发和状态变量的时间序列，m 和n分别是洪水事件和降雨事件的持续时间。进一步地，X_i＝g₁(P，E，X₀，θ′)，其中P和E中元素个数为i；X₀是此次洪水事件初始状态变量，是从水文年的起始时间算起的；g₁(·)是水文模型的子函数。θ是水文模型的参数向量，θ′是θ的子集。

产流量作为一个可以综合反映函数(1)中自变量的中间变量，数值上等于净雨，表达式为：

R_i＝f₁(P_i，E_i，X_i，θ′) (2)

其中，R_i为i时段的产流量，f₁(·)是水文模型的产流模块。由此，式(1)也可以表达为：

Q＝Q(R，θ) (3)

其中，R＝[R₁，R₂，R₃，…，R_n]是产流量的时间序列。

当任一时段的降雨、蒸散发和流域状态变量驱动水文模型时，可以得到：

q_·，i＝f(P_i，E_i，X_i，θ) (4)

其中，q_·，i＝[q_1，i，q_2，i，q_3，i，...q_m，i]^T是水文模型对R_i的响应，即R_i在流域出口的流量过程。与从式(1)到式(3)的变形类似，式(4)可做如下表达：

q_·，i＝q(R_i，θ) (5)

其中，q_j，i＝q(R_i，θ，_j)；当i≥j时,q(R_i，θ，j)＝0。为校正R_i的误差，q_·，i是本发明的校正对象。至此，如式(6)所示，流量过程可以被分解为q_·，i。此分解过程如图1所示。

式(6)的矩阵形式为：

Q＝LA (7)

其中，A＝[1，1，…，1]^T，

第二步，采用多源卫星降水产品驱动水文模型后得到第一步获得的校正对象的多个序列，将多序列的最值形成的范围作为约束，进行产流约束校正和流量预报。

(1)约束

对于任一时段，有h种卫星降雨产品可供驱动水文模型，由此可得h种产流量，即

进而可得到h种水文模型对R_i的响应，即

也就是有h种L，即L⁽¹⁾，…，L^(h)。据此可获得产流约束校正所需的约束的上限

和下限(L)，其表达式如式(8-9)所示：

(2)校正

Q_o为流域出口断面的水文站观测流量，当有m′(m′≤n＜m)个流量观测值 Q_o时，也已知m′个时段的h种径流量，因此Q_o可以做如下表示：

Q_o＝L _(m′×m′)A+ΔL_(m′×m′)A′+ε (10) 其中，

和L _(m′×m′)分别是

和L的前m’阶子阵；A′＝[α₁，α₂，…，α_m′]^T是被约束的校正系数向量，其中0≤α_i≤1，α_i越大，说明给上限更大的权重，给下限更小的权重，α_i＝0或α_i＝1说明以R_i计算得到的q近似等于下限或上限；ε＝[e₁，e₂，…，e_m′]^T是随机误差向量，根据最小二乘法，A′满足

即：

A′＝(ΔL_(m′×m′)TΔL_(m′×m′))-1ΔL_(m′×m′) ^T(Q_o-L _(m′×m′)A) (12)

此外，α_i＜0或α_i＞1分别代表校正后的q低于下限或高于上限。

为了约束校正产流量，A′中的元素也被约束，即当α_i＞1，则α_i＝1；当α_i＜0，则α_i＝0。由于L _(m′×m′)和ΔL_(m′×m′)是严格下三角矩阵，A′中仅有前(m′-T)个元素有确定解，余下的T个元素不唯一，其中T是水文模型中流量对产流响应的滞时。

(3)预报

对于预报而言，不仅A′中前(m′-T)个元素很重要，后T个元素也很重要，因为未来的流量过程受前m′个时段的所有的降雨的影响。考虑到A′中的后T个元素也是流量预报的关键，本发明令其都等于α_m′-T，即A′中第(m′-T)个元素。至此预报流量在(m′+1)和m之间的Qc时间序列可以用A′来更新，如式(13)所示：

Q_c(m-m′)＝L _{[(m-m′)×m′]}A+ΔL_{[(m-m′)×m′]}A′

其中，

和L _{[(m-m′)×m′]}分别代表

和L的m′+1～m行、1～m′列子矩阵； Q_c(m-m′)＝[Q_c(R，θ，m′+1)，Q_c(R，θ，m′+2)，…，Q_c(R，θ，m)]^T是Q_c的子向量；Q_c是校正后的预报流量的时间序列。随着m′的增加，Q_c(m-m′)持续更新，从而使得Q_c持续更新。一旦m′＞n，n个时段内所有的q都被校正了，Q_c也就不在更新。因此n次更新后可获得最终的Q_c，它即将用来率定参数和评价预报结果。

第三步，基于本校正和预报方法的参数重率定

(1)评价指标

本发明选用以下4个指标评价每场洪水的预报表现：纳什系数(NS)、洪峰相对误差(RPF)、峰现时间误差(EPT)和产流相对误差(RRD)。表达式如下：

EPT＝T_cp-T_op (16)

其中，NS由场次洪水计算得到；Q_cp和Q_op分别代表预报和观测的场次洪水的洪峰流量；T_cp和T_op分别代表预报和观测的场次洪水峰现时间；R_c和R_o分别代表预报和观测的场次洪水产流量。

(2)参数率定

θ的组成元素根据水文模型的选取而不同，并且需要优化算法进行率定。

采用NS_all，

和

作为优化目标来率定参数，其中，NS_all代表 NS由所有场次洪水计算得到，

表示RPF的绝对值的平均值，

表示EPT的绝对值的平均值，

表示RRD的绝对值的平均值。

本发明的效果和益处为：本发明通过将多源卫星降水信息转化为产流的上下限，从而对产流进行有约束的校正，可大幅提高卫星降雨产品驱动的水文模型的洪水预报精度，以支持无降雨资料地区洪水预报，并为多源卫星降雨产品的应用提供了新思路。

附图说明

图1是本发明涉及的流量过程分解，图中字母含义与公式中涉及字母的含义相同。

图2是本发明以T＝1和m’＝3为例的校正和预报运算。

图3是基于本发明的参数率定框架。

具体实施方式

本发明在以多源卫星降水信息驱动水文模型的基础上，提出一种基于多源卫星降水信息和产流约束校正的无资料地区洪水实时预报方法。

下面通过实施例，对本发明做进一步说明。

小二沟流域位于嫩江西岸，发源于内蒙古自治区大兴安岭地区，控制流域面积16761km²，汇流时间大于等于1天。流域整体呈西高东低，土壤肥沃，两岸开垦了大片农田，右岸甘南县境内建有查哈阳大型灌区，是黑龙江省最大的灌区，我国东北地区四大灌区之一。诺敏河径流补给的主要来源是降雨和季节性融雪，受温带季风气候影响，冬季寒冷干燥，夏季高温多雨，流域内降雨量年际变化较大，年内分配不均，汛期集中在6～9月，期间径流之和占全年径流 70％以上，易导致洪灾发生。以该区14场洪水为实例进行模型多场次率定和验证，其中9场洪水用于率定参数，5场洪水用于验证模型。选用4种近实时卫星降水产品，即IMERG-v06-Early(IMERG-E)、IMERG-v06-Late(IMERG-L)、 GSMaP-v06-NRT(GSMaP-N)以及GSMaP-v06-Gauge-NRT(GSMaP-GN)为实例驱动模型，表1为选用的降雨产品的基本信息。

表1选用的降水产品的基本信息

具体步骤如下：

(1)本次实例研究中选用新安江模型作为集总式水文模型，需首先根据小二沟流域的特性设置部分模型参数，如流域面积、计算步长等。本次实例研究中选用粒子群优化算法，需设置粒子群优化算法的运行参数，学习因子、惯性权重、迭代次数、搜索空间纬度等。

(2)参数率定：参数率定过程如图3所示，优化算法生成初始的模型参数粒子种群后，以4种近实时卫星降水产品驱动本发明提出的方法。当进行第m′次更新时，以T＝1，m′＝3为例的校正和预报过程如图2所示，涉及的各向量和矩阵的前三阶用于校正过程，利用最小二乘法解得校正系数向量A′，其中仅α₁和α₂为确定解；预报阶段令α₃＝α₂，以图2涉及各向量的4～m阶进行正向计算，更新Q_c的4～m个元素，即m′+1～m个元素。率定期的每场洪水利用观测流量进行 n次产流约束校正，得到每场洪水经过校正后的预报流量序列Qc。利用Qc和 Qo得到的评价指标，评价基于此参数粒子的洪水预报的精度。优化算法再更新参数种群和粒子变化方向，最终得到水文模型参数的全局最优解。

(3)模型验证：基于率定得到的模型参数最优参数，以4种近实时卫星降水产品驱动本发明提出的方法。当进行第m′次更新时，以T＝1，m′＝3为例的校正和预报过程如图2所示，涉及的各向量和矩阵的前三阶用于校正过程，利用最小二乘法解得校正系数向量A′，其中仅α₁和α₂为确定解；预报阶段令α₃＝α₂，以图2涉及各向量的4～m阶进行正向计算，更新Q_c的4～m个元素，即m′+1～m 个元素。验证期的每场洪水利用观测流量进行n次产流约束校正，得到每场洪水经过校正后的预报流量序列Qc。利用Qc和Qo得到的4个评价指标，评价基于此参数粒子的洪水预报的精度。若评价结果满足应用精度需求，则在实际工程中利用此参数向量预报洪水，否则，重新率定参数。

小二沟流域基于上述4种近实时卫星降水信息和产流约束校正的洪水预报方法，在率定期和验证期的预报效果评价指标如表2所示。可见除率定期的峰现时间误差和验证期的产流相对误差本发明方法略劣于未校正方法，其余所有基于本发明的评价指标在率定期和验证期都显著优于未校正方法的评价指标。此外，基于本发明的预报效果评价指标在率定期和验证期的表现稳定，说明基于本发明方法率定得到的参数具有代表性，能够有效描述流域的雨洪转化传递规律。

表2本发明在小二沟流域的率定期和验证期的预报效果评价指标

以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制，应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些均属于本发明的保护范围。

Claims (2)

1.一种基于多源卫星降水信息和产流约束校正的无资料地区洪水实时预报方法，其特征在于，包括以下步骤：

第一步，拆分流域出口流量组成，获得校正对象；

流域尺度的水文模型表达为：

Q＝f(P，E，X，θ) (1)

其中，Q＝[Q₁，Q₂，Q₃，…，Q_m]^T、P＝[P₁，P₂，P₃，…，P_n]^T、E＝[E₁，E₂，E₃，…，E_m]^T和X＝[X₁，X₂，X₃，…，X_m]^T分别是流量、降雨、蒸散发和状态变量的时间序列，m和n分别是洪水事件和降雨事件的持续时间；X_i＝g₁(P，E，X₀，θ′)，其中P和E中元素个数为i；X₀是此次洪水事件初始状态变量，从水文年的起始时间算起的；g₁(·)是水文模型的子函数；θ是水文模型的参数向量，θ′是θ的子集；

产流量作为能够综合反映函数(1)中自变量的中间变量，数值上等于净雨，表达式为：

R_i＝f₁(R_i，E_i，X_i，θ′) (2)

其中，R_i为i时段的产流量，f₁(·)是水文模型的产流模块；由此，式(1)也可以表达为：

Q＝Q(R，θ) (3)

其中，R＝[R₁，R₂，R₃，…，R_n]是产流量的时间序列；

当任一时段的降雨、蒸散发和流域状态变量驱动水文模型时，得到：

q_·，i＝f(P_i，E_i，X_i，θ) (4)

其中，q_·，i＝[q_1，i，q_2，i，q_3，i，...，q_m，i]^T是水文模型对R_i的响应，即R_i在流域出口的流量过程；与从式(1)到式(3)的变形类似，式(4)做如下表达：

q_·，i＝q(R_i，θ) (5)

其中，q_j，i＝q(R_i，θ，j)；当i≥j时,q(R_i，θ，j)＝0；为校正R_i的误差，q_·，i是校正对象；至此，如式(6)所示，流量过程被分解为q_·，i，则：

第二步，采用多源卫星降水产品驱动水文模型后得到第一步获得的校正对象的多个序列，将多序列的最值形成的范围作为约束，进行产流约束校正和流量预报；

(1)约束

对于任一时段，有h种卫星降雨产品可供驱动水文模型，由此可得h种产流量，即

进而得到h种水文模型对R_i的响应，即

也就是有h种L，即L⁽¹⁾，…，L^(h)；据此获得产流约束校正所需的约束的上限

和下限(L)，其表达式如式(8-9)所示：

(2)校正

Q_o为流域出口断面的水文站观测流量，当有m′(m′≤n＜m)个流量观测值Q_o时，也已知m′个时段的h种径流量，因此Q_o做如下表示：

Q_o＝L _(m′×m′)A+ΔL_(m′×m′)A′+ε (10)

其中，

和L _(m′×m′)分别是

和L的前m’阶子阵；A′＝[α₁，α₂，…，α_m′]^T是被约束的校正系数向量，其中0≤α_i≤1，α_i＝0或α_i＝1说明以R_i计算得到的q近似等于下限或上限；ε＝[e₁，e₂，…，e_m′]^T是随机误差向量，根据最小二乘法，A′满足

即：

A′＝(ΔL_(m′×m′) ^TΔL_(m′×m′))^-1ΔL_(m′×m′) ^T(Q_o-L _(m′×m′)A) (12)

此外，α_i＜0或α_i＞1分别代表校正后的q低于下限或高于上限；

为约束校正产流量，A′中的元素也需要被约束，即当α_i＞1，则α_i＝1；当α_i＜0，则α_i＝0；A′中仅有前(m′-T)个元素有确定解，余下的T个元素不唯一，其中T是水文模型中流量对产流响应的滞时；

(3)预报

未来的流量过程受前m′个时段的所有的降雨的影响，令A′中的后T个元素都等于α_m′-T，即A′中第(m′-T)个元素；至此预报流量在(m′+1)和m之间的Qc时间序列可以用A′来更新，如式(13)所示：

其中，

和L _{[(m-m′)×m′]}分别代表

和L的m′+1～m行、1～m′列子矩阵；Q_c(m-m′)＝[Q_c(R，θ，m′+1)，Q_c(R，θ，m′+2)，…，Q_c(R，θ，m)]^T是Q_c的子向量；Q_c是校正后的预报流量的时间序列；

随着m′的增加，Q_c(m-m′)持续更新，从而使得Q_c持续更新；一旦m′＞n，n个时段内所有的q都被校正了，Q_c不在更新；因此n次更新后可获得最终的Q_c，用来率定参数和评价预报结果；

第三步，基于本校正和预报方法的参数重率定

(1)选用评价指标评价每场洪水的预报表现；所述评价指标包括纳什系数NS、洪峰相对误差RPF、峰现时间误差EPT和产流相对误差RRD

(2)参数率定

θ的组成元素根据水文模型的选取而不同，并且需要优化算法进行率定；

采用NS_all，

和

作为优化目标来率定参数，其中，NS_all代表NS由所有场次洪水计算得到，

表示RPF的绝对值的平均值，

表示EPT的绝对值的平均值，

表示RRD的绝对值的平均值。

2.根据权利要求1所述的一种基于多源卫星降水信息和产流约束校正的无资料地区洪水实时预报方法，其特征在于，第三步步骤(1)中评价指标纳什系数(NS)、洪峰相对误差(RPF)、峰现时间误差(EPT)、产流相对误差(RRD)的表达式如下：

EPT＝T_cp-T_op (16)

其中，NS由场次洪水计算得到；Q_cp和Q_op分别代表预报和观测的场次洪水的洪峰流量；T_cp和T_op分别代表预报和观测的场次洪水峰现时间；R_c和R_o分别代表预报和观测的场次洪水产流量。

Images