CN106971034A - 一种基于集合卡尔曼滤波的对无资料地区径流量推求方法 - Google Patents

Abstract

本发明公开了一种基于集合卡尔曼滤波的对无资料地区径流量推求方法，涉及水文模拟技术领域。所述方法：计算全流域的优化参数值；计算子流域A的出口点C点的估计流量及其分布规律；计算变量和参数的集合扰动值；计算下游子流域B的出口点即水文站D的径流量预报集合；对水文站D的实测径流值进行高斯扰动，并将得到高斯扰动值作为观测数据集合；融入观测数据集合进行同化更新；最终得到C点径流量Q_C的最优估计值。本发明简化了滤波同化过程，较好的提高了同化精度；同时也避免了常规方法对无资料地区求径流量的缺陷。

Description

一种基于集合卡尔曼滤波的对无资料地区径流量推求方法

技术领域

本发明涉及水文模拟技术领域，尤其涉及一种基于集合卡尔曼滤波的对无资料地区径流量推求方法。

背景技术

自然界中有不少区域因雨量或水文站早期建设少、历史径流资料序列短，无法满足预报需求；同时，因人类活动等影响，导致流域产汇流特性发生较大变化，历史资料无法反映当前来水特性，历史资料不可用。因此对无资料地区的研究也是亟待解决的。

目前，对无资料地区径流预报常用的方法为区域化方法，即通过某种途径，利用有资料流域的模型参数推求无资料流域的模型参数，从而对无资料流域进行预报。常用的区域化方法有属性相似法、空间相近法和回归法。其中，属性相似法指找出与研究流域属性如地形、土壤和气候等相似的流域，并把其参数作为研究流域的参数，属性相似法研究根据为同一区域的物理和气候属性相对一致，因此相邻流域的水文行为相似。空间相近法指找出与研究流域(如无资料流域)距离上相近的一个或者多个有资料流域，并把有资料流域参数作为研究流域的参数。回归法指根据有资料流域的模型参数和流域属性，建立二者之间的多元回归方程，从而利用无资料流域的流域属性推求其模型参数。

但是，现有技术存在如下缺点：

一、未考虑不同参数的区别：目前，无资料地区参数确定通常采用移用或回归分析的方法。现有，确定移用时通常先找到参照流域，而后所有参数统一移用或回归分析。可是实际上，水文模型中每个参数代表的物理意义不同，因每个流域都有其独特性，很难找到与模型中所有参数的相关特征均相似的流域，虽现今已有研究仅依据一种或几种相关特征寻找参照流域，而所有参数统一移用或采用同一些流域，但这种方式很明显无法反映参数的独特性。

二、相似流域的选取缺乏全面性和客观性：现有研究通常选用一种或几种地形、地貌特征参数的综合指数如地形指数作为属性判别寻找相似流域，而后进行参数的移用或分析。因为不同水文模型其参数的物理意义不同，仅仅依据一种或几种流域特征参数进行判别是不够全面和客观的，很难反映真正的流域相似，并且采用所有流域特征进行分析存在信息量大、分类复杂以及无法判断的问题。

三、主观判别和效率低：现有技术均采用人工判别处理的方式，解决无资料地区的洪水预报问题，此技术难以大量推广，并且其中包含人为的主观性和效率低下的问题，寻找其他有效简便的方法是目前急需解决的。

发明内容

本发明的目的在于提供一种基于集合卡尔曼滤波的对无资料地区径流量推求方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明所述基于集合卡尔曼滤波的对无资料地区径流量推求方法，所述方法包括：

S1，计算全流域的优化参数值

全流域中存在至少一个无水文资料子流域和至少一个有水文资料子流域，所述无水文资料的子流域的上游流域和下游流域均为全流域中子流域；

对全流域进行参数优化，得到参数最优值，选取任意两个参数的最优值进行数据同化处理；两个参数分别为参数β和参数γ；

S2，计算子流域A的出口的估计流量及其分布规律

设定无水文资料子流域为子流域A，子流域A的出口点为C点，子流域A的下游流域为有水文资料子流域B，子流域B中的出口点为水文站D；

子流域B出口对应的径流量是已知的，记为Q_D；

然后，通过公式(1)计算子流域A的出口点C点的初始估计流量Q_C′；

Q_C′＝Q_D×(S_A/S_A+B) (1)；

其中，S_A为子流域A的面积，S_A+B为子流域A和子流域B的面积之和；

S3，计算变量和参数的集合扰动值

将初始估计流量Q_C′的均值和方差作为C点径流量Q_C初始设定的符合高斯分布的均值与方差；以参数β和参数γ的最优值为参考，对C点径流量Q_C、参数β和参数γ分别进行高斯扰动，分别得到变量、参数β和参数γ的高斯扰动值集合；

S4，计算下游有水文资料水文站D的径流量预报集合

在实测降水和蒸发数据的基础上，逐天运行水文模型，将S3中得到的C点径流量Q_C、参数β和参数γ的高斯扰动集合值逐一代入水文模型，并得到每个时间步长上的水文站D的径流量预报集合，每个时间步长上的水文站D的径流量预报集合数量大于等于1；

S5，对水文站D实测径流值进行高斯扰动，并将得到的高斯扰动集合作为观测数据集合；

S6，融入观测数据集合

在已知观测算子的条件下，计算卡尔曼增益，将每个所述水文站D的径流量预报集合融入观测数据集合进行更新；

S7，重复进行S4至S6，直到水文模型计算结束，最终得到C点径流量QC的最优估计值，所述最优估计值即为上游无资料子流域A出口点C点的径流量值。

优选地，S1中，对全流域进行参数优化采用的优化算法为所述粒子群算法。

优选地，所述水文模型为新安江模型。

优选地，步骤S3中，所述高斯扰动，按照公式(2)实现：

Ea＝ea+ε (2)；

其中，Ea表示增加高斯扰动后的模型状态变量或者观测值；ea表示初始模型状态变量或者观测值；ε为符合高斯分布的扰动值，即ε～N(0，1)。

优选地，步骤S6中，计算卡尔曼增益，具体按照下述计算：

K_t＝P_tH^T(HP_tH^T+R_t)^-1 (9)；

其中，m表示集合样本的总数；表示t时刻m个径流量预报集合的预报均值；表示t时刻水文模型第i个径流量预报集合中集合成员的一个预报状态向量；表示t时刻每一个径流量预报集合的集合成员与的扰动；y_i,t表示t时刻第i个集合成员的观测数据向量；表示t时刻m个观测数据集合的观测均值；Δy_i,t表示t时刻每一个观测数据的集合成员与的扰动；P_t表示t时刻的径流量预报集合的误差协方差；R_t表示t时刻观测数据集合的误差协方差矩阵；H为观测算子；K_t为卡尔曼增益，表示观测数据的权重；T表示观测算子的转置。

优选地，步骤S6中，对每一个集合进行更新，采用公式(10)进行更新：

其中，表示t时刻状态下的水文站D的径流量的最优估计集合，的均值即为t时刻状态下的最优估计值，表示t时刻的水文站D的径流量预报集合，y_t表示t时刻的观测数据集合，K_t为卡尔曼增益，表示观测数据的权重。

本发明的有益效果是：

本发明利用集合卡尔曼滤波推求无资料地区径流量的研究，为无资料地区的研究展开了新思路和新方法。

集合卡尔曼滤波的基本思想即是利用Monte Carlo方法设计预测状态的一个集合，该集合的平均可作为最佳估计，该集合的样本协方差即作为预测误差协方差的近似，该集合通过不断向前滤波，每个样本分别更新分析变量，而对变量的最佳估计即为各更新分析变量的样本平均。本发明提出的基于集合卡尔曼滤波的对无资料地区径流量推求方法，基于集合统计的思想，无须对非线性系统进行线性化，避免了Jacobian矩阵的繁冗计算，更简化了滤波同化过程。

本发明利用一种扩展模型状态变量方法，即将模型参数和变量置于联合的向量中，将上游无资料地区径流量看作变量，利用下游已知流域的参数并将其径流量作为观测值，同时实现变量与参数的同化，其中变量同化的均值即为无资料地区的径流量。基于集合统计的思想，无须对非线性系统进行线性化，避免了Jacobian矩阵的繁冗计算，并且同化过程中同时考虑输入降雨误差、模型误差以及观测误差，较好的提高了同化精度。同时也避免了常规方法对无资料地区求径流量的缺陷，如相似流域选取缺乏全面、客观性，现有技术未考虑不同参数的区别，主观判别、效率低等，该方法将会对无资料地区预报完善洪水预报理论、提高预报精度、水库大坝建设和调度设计等具有理论价值和指导意义。

附图说明

图1是全流域子流域分布及水文站位置简图；

图2是基于集合卡尔曼滤波的对无资料地区径流量推求方法的流程示意图；

图3是参数优化后出口点水文站D的径流模拟值与真实径流值对比示意图；

图4是基于集合卡尔曼滤波下的同化估计上游无资料地区的径流量与真实径流量对比示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

实施例

本实施例所述基于集合卡尔曼滤波的对无资料地区径流量推求方法，所述方法包括：

S1，计算全流域的优化参数值

全流域中存在至少一个无水文资料子流域和一个有水文资料子流域，所述无水文资料的子流域的上游流域和下游流域均为全流域中子流域；

对全流域进行参数优化，得到参数最优值，选取任意两个参数的最优值进行数据同化处理；两个参数分别为参数β和参数γ；

S2，计算子流域A的出口的估计流量及其分布规律，所述分布规律包括估计流量的均值和方差

设定无水文资料子流域为子流域A，子流域A的出口为C，子流域A的下游流域为有水文资料子流域B，子流域B中的水文站为水文站D；

子流域B出口对应的径流量是已知的，记为Q_D；

然后，通过公式(1)计算子流域A的出口的初始估计流量Q_C′；

Q_C′＝Q_D×(S_A/S_A+B) (1)；

其中，S_A为子流域A的面积，S_A+B为子流域A和子流域B的面积之和；

S3，计算变量和参数的集合扰动值

将初始估计流量Q_C′的均值和方差作为C点径流量Q_C初始设定的符合高斯分布的均值与方差；以参数β和参数γ的最优值为参考，对C点径流量Q_C、参数β和参数γ分别进行高斯扰动，分别得到变量、参数β和参数γ的高斯扰动值集合；

S4，计算下游有水文资料水文站D的径流量预报集合

在实测降水和蒸发数据的基础上，逐天运行水文模型，将S3中得到的C点径流量Q_C、参数β和参数γ的高斯扰动值集合中的值逐一代入水文模型，并得到每个时间步长上的水文站D的径流量预报集合，每个时间步长上的水文站D的径流量预报集合数量大于等于1；

S5，对水文站D的实测径流值进行高斯扰动，并将得到高斯扰动集合作为观测数据集合；

S6，融入观测数据集合

在已知观测算子的条件下，计算卡尔曼增益，将每个所述水文站D的径流量预报集合融入观测数据集合进行更新；

S7，重复进行S4至S6，直到水文模型计算结束，最终得到C点径流量Q_C的最优估计值即为上游无资料子流域A出口点C点的径流量值。

更详细的解释说明为：

(一)S1中，对全流域进行参数优化采用的优化算法为所述粒子群算法。通过迭代的形式求得最优解。此次实施例选取种群规模为200，种群进化迭代次数为3000。

(二)所述水文模型为新安江模型。新安江模型适用性较广。本实施例中的两个参数为IMP(流域不透水面积)、KG(地下水出流系数)。

在本实施例中利用一种扩展模型状态变量方法，即将模型参数和变量置于联合的向量中，此次我们的实施例选取的流域，上游A无资料地区为雅砻江流域的甘孜水文站(假设其径流量未知，最后将真实径流值与同化估计值作对比，来验证同化结果的好坏)，下游B地区为雅江。将上游A无资料地区甘孜站的径流量看作变量，记作q，选取下游已知流域的两个参数，则状态变量的形式为X＝(q，Q，IMP，KG)，从实际理论中得知甘孜站的径流量q加上B流域的产流记为R_B，将总产流(q+R_B)代入新安江模型的汇流模块得到下游雅江站的径流量Q。同时，初始给定状态变量符合高斯分布的均值与方差。其中IMP～N(0.3，0.2²)，KG～N(0.4，0.2²)。q的均值和方差可根据上下游面积比例法求得，因为下游B雅江的径流量D是已知的，则得到上游A无资料地区甘孜站的初始估计流量q为D×(A/(A+B))，并由此推求q的分布规律(均值、方差)，即得q～N(390，10²)。本次发明中还考虑了模型误差以及观测误差，并用缩放因子的形式展现在程序中，更简便实用，模型及观测误差同样满足高斯分布，即服从N(0，0.01²)，N(0，0.1²)来表示模型和观测的不确定性。

步骤S2，实施例选取了三年日降雨蒸发数据，共计1095天，第一年作为预热期，第二年与第三年进行同化。逐天运行新安江模型，利用步骤S1得到的参数最优值，代入新安江模型中，得到模拟径流值，并将实际测量的径流量作为观测值，为之后进行更新同化模拟径流量数据。其中图3为参数优化后全流域的径流模拟值与真实径流值对比图。

(三)步骤S3中，所述高斯扰动，按照公式(2)实现：

Ea＝ea+ε (2)；

其中，Ea表示增加高斯扰动后的模型状态变量或者观测值；ea表示初始模型状态变量或者观测值；ε为符合高斯分布的扰动值，即ε～N(0，1)。

在步骤S3中，根据初始设定的符合高斯分布的均值与方差，得到变量及两参数的集合高斯扰动值，并对高斯扰动值进行范围设定，保证同化的精确。

(四)步骤S4中，逐天运行新安江模型，每运行一天即产生100个随机集合，将步骤S3得到的变量和参数的集合扰动值以及实测降水与蒸发数据代入水文模型中，计算得到下游已知站点的径流量预报集合。考虑模型误差，使用者还可对径流量预报集合进行高斯扰动，得到径流量预报的扰动集合。

(五)步骤S6中，考虑符合高斯分布的观测误差，又已知观测算子，此次观测算子H＝(0，1，0，0)^T，计算卡尔曼增益，将上一步骤(5)得到的径流量扰动集合融入观测值，对每一个集合进行更新，得到C点径流量Q_C的最优估计值即为上游无资料子流域A出口点即甘孜站的径流量值。

具体按照下述计算：

K_t＝P_tH^T(HP_tH^T+R_t)^-1 (9)；

其中，m表示集合样本的总数；表示t时刻m个径流量预报集合的预报均值；表示t时刻水文模型第i个径流量预报集合中集合成员的一个预报状态向量；表示t时刻每一个径流量预报集合的集合成员与的扰动；y_i,t表示t时刻第i个集合成员的观测数据向量；表示t时刻m个观测数据集合的观测均值；Δy_i,t表示t时刻每一个观测数据的集合成员与的扰动；P_t表示t时刻的径流量预报集合的误差协方差；R_t表示t时刻观测数据集合的误差协方差矩阵；H为观测算子；K_t为卡尔曼增益，表示观测数据的权重。T表示观测算子的转置；

同时，步骤S6中，对每一个集合进行更新，采用公式(10)进行更新：

其中，表示t时刻状态下的子流域B出口点水文站D的径流量的最优估计集合，的均值即为t时刻状态下的最优估计，表示t时刻的子流域B出口点水文站D的径流量预报集合，y_t表示t时刻的观测数据集合，K_t为卡尔曼增益，表示观测数据的权重。

本实施例中的时间步长为2年，共730天，每计算一天都会算出100个集合数的径流量，更新时对100个集合分别进行更新，得到的均值就为该天的最优估计值。同样依次类推，下一天同样计算出100个集合。同化更新完一天的数据后，再重复运行下一天，产生100个预报集合，再进行更新求得变量最优估计值，以此类推，直到730天结束。在同化过程中，每个集合成员均为独立进行的。

经过整个同化更新过程后，可以得到图4。其中图4为基于集合卡尔曼滤波下的同化估计上游无资料地区的径流量与真实径流量对比图。

由图4基于集合卡尔曼滤波下的求解的同化估计上游无资料地区的径流量与真实径流量对比图，可以看出同化估计后的径流值基本与真实值吻合一致，其两者的相对误差均小于1，验证了同化估计的可行性。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：

本发明所述基于集合卡尔曼滤波的对无资料地区径流量推求方法，利用一种扩展模型状态变量方法，即将模型参数和变量置于联合的向量中，将上游无资料地区径流量看作变量，利用下游已知流域的参数并将其径流量作为观测值，同时实现变量与参数的同化，其中变量同化的均值即为无资料地区的径流量。基于集合统计的思想，无须对非线性系统进行线性化，避免了Jacobian矩阵的繁冗计算，并且同化过程中同时考虑输入降雨误差、模型误差以及观测误差，较好的提高了同化精度。同时也避免了常规方法对无资料地区求径流量的缺陷，如相似流域选取缺乏全面、客观性，现有技术未考虑不同参数的区别，主观判别、效率低等。

本发明所述方法将会对无资料地区预报完善洪水预报理论、提高预报精度、水库大坝建设和调度设计等具有理论价值和指导意义。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (6)

1.一种基于集合卡尔曼滤波的对无资料地区径流量推求方法，其特征在于，所述方法包括：

S1，计算全流域的优化参数值

全流域中存在至少一个无水文资料子流域和至少一个有水文资料子流域，所述无水文资料的子流域的上游流域和下游流域均为全流域中子流域；

对全流域进行参数优化，得到参数最优值，选取任意两个参数的最优值进行数据同化处理；两个参数分别为参数β和参数γ；

S2，计算子流域A的出口的估计流量及其分布规律

设定无水文资料子流域为子流域A，子流域A的出口点为C点，子流域A的下游流域为有水文资料子流域B，子流域B中的出口点为水文站D；

子流域B出口对应的径流量是已知的，记为Q_D；

然后，通过公式(1)计算子流域A的出口点C点的初始估计流量Q_C′；

Q_C′＝Q_D×(S_A/S_A+B) (1)；

其中，S_A为子流域A的面积，S_A+B为子流域A和子流域B的面积之和；

S3，计算变量和参数的集合扰动值

将初始估计流量Q_C′的均值和方差作为C点径流量Q_C初始设定的符合高斯分布的均值与方差；以参数β和参数γ的最优值为参考，对C点径流量Q_C、参数β和参数γ分别进行高斯扰动，分别得到变量、参数β和参数γ的高斯扰动值集合；

S4，计算下游有水文资料水文站D的径流量预报集合

在实测降水和蒸发数据的基础上，逐天运行水文模型，将S3中得到的C点径流量Q_C、参数β和参数γ的高斯扰动集合值逐一代入水文模型，并得到每个时间步长上的水文站D的径流量预报集合，每个时间步长上的水文站D的径流量预报集合数量大于等于1；

S5，对水文站D的实测径流值进行高斯扰动，并将得到的高斯扰动集合作为观测数据集合；

S6，融入观测数据集合

在已知观测算子的条件下，计算卡尔曼增益，将每个所述水文站D的径流量预报集合融入观测数据集合进行更新；

S7，重复进行S4至S6，直到水文模型计算结束，最终得到C点径流量Q_C的最优估计值，所述最优估计值即为上游无资料子流域A出口点的径流量值。

2.根据权利要求1所述方法，其特征在于，S1中，对全流域进行参数优化采用的优化算法为所述粒子群算法。

3.根据权利要求1所述方法，其特征在于，所述水文模型为新安江模型。

4.根据权利要求1所述方法，其特征在于，步骤S3中，所述高斯扰动，按照公式(2)实现：

Ea＝ea+ε (2)；

其中，Ea表示增加高斯扰动后的模型状态变量或者观测值；ea表示初始模型状态变量或者观测值；ε为符合高斯分布的扰动值，即ε～N(0，1)。

5.根据权利要求1所述方法，其特征在于，步骤S6中，计算卡尔曼增益，具体按照下述计算：

${\stackrel{\‾}{\hat{x}}}_{i,t}=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{\hat{x}}_{i,t}---\left(3\right);$

$\mathrm{\Δ}{\hat{x}}_{i,t}={\hat{x}}_{i,t}-{\stackrel{\‾}{\hat{x}}}_{i,t}---\left(4\right);$

${\stackrel{\‾}{y}}_{i,t}=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}{y}_{i,t}---\left(5\right);$

${\mathrm{\Δy}}_{i,t}=y_{i,t}-{\stackrel{\‾}{y}}_{i,t}---\left(6\right);$

$P_t=\frac{1}{m-1}\[\mathrm{\Δ}{\hat{x}}_{i,t}{\left(\mathrm{\Δ}{\hat{x}}_{i,t}\right)}^T\]---\left(7\right);$

$R_t=\frac{1}{m-1}\[{\mathrm{\Δy}}_{i,t}{\left({\mathrm{\Δy}}_{i,t}\right)}^T\]---\left(8\right);$

K_t＝P_tH^T(HP_tH^T+R_t)^-1 (9)；

其中，m表示集合样本的总数；表示t时刻m个径流量预报集合的预报均值；表示t时刻水文模型第i个径流量预报集合中集合成员的一个预报状态向量；表示t时刻每一个径流量预报集合的集合成员与的扰动；y_i,t表示t时刻第i个集合成员的观测数据向量；表示t时刻m个观测数据集合的观测均值；Δy_i,t表示t时刻每一个观测数据的集合成员与的扰动；P_t表示t时刻的径流量预报集合的误差协方差；R_t表示t时刻观测数据集合的误差协方差矩阵；H为观测算子；K_t为卡尔曼增益，表示观测数据的权重；T表示观测算子的转置。

6.根据权利要求1所述方法，其特征在于，步骤S6中，对每一个集合进行更新，采用公式(10)进行更新：

${\tilde{\hat{x}}}_t={\hat{x}}_t+K_t\[y_t-H\left({\hat{x}}_t\right)\]---\left(10\right);$

其中，表示t时刻状态下的水文站D的径流量的最优估计集合，的均值即为t时刻状态下的最优估计值，表示t时刻的水文站D的径流量预报集合，y_t表示t时刻的观测数据集合，K_t为卡尔曼增益，表示观测数据的权重。