本发明涉及一种基于集合卡尔曼滤波的洪水预报方法,属于水文学领域。
背景技术：
洪水预报在防洪非工程措施中是一种重要内容，为防洪减灾、水资源利用与保护、水利工程建设与调度、工农业安全生产等服务，因此洪水预报的准确性和及时性有着重要意义。实时洪水预报是一种在联机水情测报系统中，使用实时雨、水情及其它有关水文气象信息作为洪水预报模型输人，并不断根据新信息校正或改善原有模型参数，力争预报结果逐步逼近真值的洪水预报。
然而利用水文模型进行预报时，误差总是存在的。预报误差包括模型自身的误差及输入资料的误差、模型初始状态误差，对预报值进行误差实时校正是十分必要的。多年来，水文学者对此进行了大量尝试，如ar模型、最小二乘法、系统动态响应曲线、卡尔曼滤波等等。但其中主要用于修正流域中间状态量的是卡尔曼滤波。
20世纪70年代初到80年代末，卡尔曼滤波技术被广泛应用于水文学、水动力学以及水质等方面的研究中。然而卡尔曼滤波在应用时会受制于线性系统的前提假设,该算法需要线性的模型算子和观测算子,对于复杂的非线性系统(如新安江模型)，给出线性的模型算子是非常困难的。基于蒙特卡洛法生成样本集合的集合卡尔曼方法避免了线性化系统这一步，同时在更新系统状态时不再需要事先确定模型协方差，加之其在计算机上易于实现而迅速在众多水文模型中得到应用。为了探索如何将集合卡尔曼与水文模型结合起来并得到最优的预报结果，许多水文学者做了大量的尝试和分析。在洪水预报中，中国应用最多、最广的水文模型就是新安江模型；现有技术中，将新安江模型与集合卡尔曼滤波相结合，在修正过程中模型预测集合可能出现集体偏大或集体偏小的情况，这造成最后的修正效果达不到理想程度。
技术实现要素：
发明目的：为了提高洪水预报的精度，本发明的目的在于提出一种基于集合卡尔曼滤波的洪水预报方法。
本发明技术方案如下。
一种基于集合卡尔曼滤波的洪水预报方法，包括以下步骤：
s1，根据流域的雨量站个数n,选取集合卡尔曼滤波器的系统状态量和观测量：系统状态量包括土壤含水量w、自由水蓄量s、产流面积fr、地表径流qs、壤中流qss、地下径流qs和模型计算流量值，观测量为流域出口实测流量值；
s2,设样本总数为n：对降雨量p、蒸发量ep、土壤含水量w、自由水蓄量s、产流面积fr、地表径流qs、壤中流qss和地下径流qg进行服从高斯分布的噪声干扰，分别形成由n个样本组成的集合；
s3,对系统状态量进行预报，得到水文模型预测集合；
通过观测转移矩阵h建立系统状态量和实测流量值之间的关系：
xt,i＝φ[xt-1,i,wt-1,i](1)
zt,i＝h(xt,i,vt,i)(2)
式中，xt,i表示t时刻第i个样本的水文模型预测集合，即t时刻，第i个样本的状态变量；xt,i由第t-1时刻第i个样本的状态变量xt-1,i及系统噪声wt-1,i推算，φ[]表示状态变量xt-1,i及系统噪声wt-1,i推算xt,i的推算函数；zt,i表示t时刻第i个样本的观测值，观测值zt,i通过观测转移矩阵h与状态变量xt,i和观测噪声vt,i建立关系；
s4,按比例选取水文模型预测集合xt,i中与观测值最接近的部分样本为真值，样本选取好后，选取的样本数目为m,用选取样本平均值
其中，et为模型预测集合与真实值之间的偏差(即预测的系统状态量和实测流量值之间的差值)，ett为et的转置矩阵；
s5，根据实测流量值的方差r，结合预测协方差矩阵pt，得到卡尔曼增益矩阵kt，并更新得到分析值xt,i，形成模型分析集合：
kt＝ptht(hptht+r)-1(6)
xt,i＝xt,i+kt[zt,i-hxt,i](7)
xt,i作为集合中n个样本更新后的系统状态量；
s6,对xt,i的样本集合取均值，得到分析值均值；判断t<l,如成立则重复步骤s3-s5,直至t<l等式不成立，则跳出循环，得到实时洪水预报的计算流量，其中t代表时段时刻；l代表结束时刻。
土壤含水量w、自由水蓄量s、产流面积fr、地表径流qs、壤中流qss、地下径流qg在初始时刻t＝1时的干扰表示初值存在不确定性。
较优地，每一时段对于降雨量p和蒸发量ep的噪声干扰表示在测量时存在不确定性；
初值存在不确定性和测量时存在不确定性通过影响每一时刻的模型中间状态量影响模型预报的出口断面流量，即模型计算流量q。
较优地，模型计算流量q是由水文模型计算得到的，是所有不确定性综合后反映到出口断面流量的不确定性；对实测流量值q0干扰成n个样本。
系统噪声和观测噪声皆为高斯白噪声。
较优地，步骤4中选取比例为10％-20％。
本发明的有益效果包括：
本发明提出了一种基于集合卡尔曼滤波的洪水预报方法，主要包括对选取流域中间状态量、降雨量、蒸发量进行干扰形成样本集合，通过新安江模型预报得到水文模型预测集合，由选取水文模型预测集合最接近实测值的一部分集合的平均值代替真实值计算得到预测协方差矩阵，然后根据实测值的方差结合模型预测集合得到卡尔曼增益矩阵，同时更新预报值，得到模型分析值。在将新安江水文模型与集合卡尔曼滤波相结合中根据流域出口断面计算流量与实测值之间的误差，得到流域出口断面计算流量的分析值同时对流域中间状态量进行修正时模型预测集合可能出现集体偏大或集体偏小的情况，这造成最后的修正效果达不到理想程度。因此，本申请同时对集合卡尔曼滤波进行一定的改进，使其修正效果更好，达到更优的洪水预报结果，提高洪水预报精度。
本发明公开一种基于集合卡尔曼滤波的洪水预报方法，对集合卡尔曼滤波进行一定的改进，使其修正效果更好，达到更优的洪水预报结果，提高洪水预报的精度；与现有技术误差方法相比，本发明可以应用于非线性系统，修正系统状态量，并准确预报之后的流域中间状态量和流量。
本申请在集合卡尔曼滤波与模型结合使用过程中，当某一时刻模型预测集合出现集体偏大或集体偏小时，修正的更精确，同时使之后的预报过程更精确。
附图说明
下面结合附图和实施例对本发明进一步说明；
图1是本发明改进后的集合卡尔曼滤波结构图。
图2是本发明与现有滤波方法的一次洪水预报的结果对比图。
具体实施方式
下面结合附图并通过具体实施例对本发明作进一步详述，以下实施例只是描述性的，不是限定性的，不能以此限定本发明的保护范围。
为了使本发明的技术手段、创作特征、工作流程、使用方法达成目的与功效，且为了使该评价方法易于明白了解，下面结合具体实施例，进一步阐述本发明。
如图1所示，一种基于集合卡尔曼滤波的洪水预报方方法：将改进集合卡尔曼滤波应用于新安江模型上，通过流域出口断面计算流量与实测值之间的误差，得到流域出口断面计算流量的分析值，同时对新安江模型土壤含水量、自由水蓄量、产流面积、地表径流、壤中流、地下径流进行修正，逐步降低误差的积累，提高洪水预报的精度。包括以下步骤：
s1，根据流域的雨量站个数n,选取集合卡尔曼滤波器的系统状态量和观测量：系统状态量包括土壤含水量、自由水蓄量、产流面积、地表径流、壤中流、地下径流和模型计算流量值，观测量为流域出口实测流量值；
s2,设样本总数为n：对降雨量p、蒸发量ep、土壤含水量w、自由水蓄量s、产流面积fr、地表径流qs、壤中流qss和地下径流qg进行服从高斯分布的噪声干扰，分别形成由n个样本组成的集合；
其中土壤含水量w、自由水蓄量s、产流面积fr、地表径流qs、壤中流qss、地下径流qg在初始时刻t＝1时的干扰表示初值存在一定程度不确定性，而每一时段对于降雨量p和蒸发量ep的噪声干扰是因其在测量时存在一定误差，而这些不确定性会通过影响每一时刻的新安江模型中间状态量，影响模型预报的出口断面流量即模型计算流量q。
模型计算流量q是由模型计算得到的，是所有不确定性综合后反映到出口断面流量的不确定性；同样的，对实测流量值q0干扰成n个样本；
s3,对系统状态量进行预报，得到水文模型预测集合；通过观测转移矩阵h建立系统状态量和实测流量值之间的关系：
xt,i＝φ[xt-1,i,wt-1,i](1)
zt,i＝h(xt,i,vt,i)(2)
式中，xt,i表示t时刻第i个样本的水文模型预测集合，即t时刻第i个样本的系统状态量，xt,i由第t-1时刻第i个样本的系统状态量xt-1,i及系统噪声wt-1,i推算；zt,i表示t时刻第i个样本的观测值，观测值zt,i通过h矩阵与系统状态量xt,i和观测噪声vt,i建立关系；系统噪声和观测噪声皆为高斯白噪声。
s4,按比例选取水文模型预测集合xt,i中与观测值最接近的部分样本为真值，选取比例根据实际情况可为10％或者20％，但需考虑当样本数比例过大，效果没有比例小的好且起不到改进效果，当样本比例选取过小，会造成每次模拟效果不同。当样本选取好后，选取的样本数目为m,用选取样本平均值
其中，et为模型预测集合与真实值之间的偏差，ett为et的转置矩阵；
s5，根据实测流量值的方差r，结合预测协方差矩阵pt，得到卡尔曼增益矩阵kt，并更新得到分析值xt,i，形成模型分析集合：
kt＝ptht(hptht+r)-1(6)
xt,i＝xt,i+kt[zt,i-hxt,i](7)
xt,i作为集合中n个样本更新后的系统状态量。
s6,对xt,i的样本集合取均值，得到分析值均值。判断t<l,如成立则重复s3-s5,直至等式不成立，则跳出循环，得到实时洪水预报的计算流量，t代表时段时刻；l代表结束时刻。
图2是本发明与现有滤波方法的一次洪水预报的结果对比图，对比本申请能够更准确进行洪水预测。
本领域内的技术人员可以对本发明进行改动或变型的设计但不脱离本发明的思想和范围。因此，如果本发明的这些修改和变型属于本发明权利要求及其等同的技术范围之内，则本发明也意图包含这些改动和变型在内。