一种用于降水量估算的卡尔曼滤波改进方法与流程

本发明涉及气象科学，特别是降水量的估算方法。

背景技术：

目前，常规探空资料依然被作为我国对短时强降水等灾害性天气数值预报模式的初始资料。然而，虽然在20世纪60年代以后，我国不断加大自动气象观测站站点建设，加密自动气象站观测网，但这些传统气象探空资料在反映中小尺度降水的时间间隔和空间分布方面仍存在很大的困难，同时由于地形等因素的影响，这种雨量信息往往不能代表真实雨量的分布。

利用卡尔曼滤波能够实现对区域降水总量的估算，将传统的以毫米为计量单位的点状计量方式转换成以吨位计量的面状计量方式。通过对降水进行统计研究发现，自动气象站与人工观测站联合观测插值得到的面降雨总量与通过卡尔曼滤波法估算得到的面降雨总量，误差可以控制在15％以内。如何进一步的提高卡尔曼滤波估算的准确度和稳定性一直是本领域内研究人员的研究重点之一。

技术实现要素：

发明目的：为了克服现有技术中存在的不足，本发明提供一种用于降水量估算的卡尔曼滤波改进方法，用于解决现有的卡尔曼滤波估算面降雨总量的方法误差较大且稳定性差的技术问题。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种用于降水量估算的卡尔曼滤波改进方法，在采样频率不变的前提下，给定如下初始值，包括和P_k，这里k＝0，是根据第一次测量值设定的，P₀是该测量值的误差方差，P₀一般设为0.1然后顺序执行以下步骤：

步骤一、利用(1)式进行状态预测，即求解状态预测估计矩阵

(1)式中,

为第k次状态预测估计矩阵，应用在降水量估计中，即指不同时间的降水量的估值；

Y_k为中间计算值，且Y_k＝E_k-H_kK_k，这里E_k为单位矩阵，H_k为第k次系统观测矩阵；

K_k为第k次系统滤波增益矩阵；

如果采样频率确定，则Y_k和K_k各个元素均为常数；

U_k为第k次系统观测值，能够直接得到；

步骤二、利用(2)式进行误差方差预测，即求解预测误差方差矩阵

(2)式中，

P_1k为第k次一步预测误差方差矩阵；

Φ_k为第k次系统状态转移矩阵；为Φ_k的转置；

Γ_k为第k次系统噪声矩阵；为Γ_k的转置；

H_k-1为第k-1次系统观测矩阵；

K_k-1为第k-1次系统滤波增益矩阵；

Q_k-1为第k-1次系统协方差矩阵；

步骤三、利用(3)式确定系统滤波增益矩阵

(3)式中，

K_k为第k次系统滤波增益矩阵；Q_k为第k次系统协方差矩阵；H_k为第k次系统观测矩阵,为H_k的转置；

步骤四、利用(4)式求解估计误差方差矩阵

P_k＝(E_k-K_kH_k)P_1k (4)

(4)式中，

P_k为第k次估计误差方差矩阵；

步骤五、令k＝k+1，若满足k〉N，则算法结束；若满足k≤N，则返回步骤一；这里N是指总共的降水量测量次数，为预设值。

有益效果：

本发明提供一种用于降水量估算的卡尔曼滤波改进方法，使用过程中确定采样频率不变的情况下，通过一次计算之后，那么大量复杂的滤波增益矩阵和估计误差方差矩阵可以预先离线完成，可以大大减少运算量，随着k不断增加的循环计算过程中都可以直接应用第一次计算获得的系统滤波增益矩阵和估计误差方差矩阵代入计算相关参数，避免了以往的算法中反复迭代的过程；

本发明方法能够实际应用于多普勒天气雷达降水估算研究，探测的结果比较符合实际情况，与原有的卡尔曼滤波估算法相比，估算稳定性与估算精度更高，可以满足实际应用的要求。

附图说明

图1为本发明的卡尔曼滤波改进方法流程图；

图2为各降水场次相对误差对比图；

图3为卡尔曼滤波法和自动站获得的时雨量相关性；

图4为卡尔曼滤波改进法和自动站获得的时雨量相关性。

具体实施方式

下面结合附图对本发明作更进一步的说明。

利用本发明进行广州市6次降水模拟与其他算法对比试验。

选取2013年3月间6场降水过程进行估算方法的实验检验。6次降雨过程时间为别为3月1日、3月20日、3月23日、3月26日、3月28日、3月31日，共计35小时累计降雨时次。本发明试验采用的是基于广州市域雷达半径110km之内，雷达回波资料取自投入实际业务运行的广州市CINRAD/SA多普勒天气雷达降水模式中VCP21采集到的每6min一次的雷达基本反射率数据。

步骤一、对雷达数据质量控制和数据预处理。

(1)根据现有经验，雷达的降水回波强度低于20dBZ的，有可能是降水，也有可能是由于雾或其它悬浮在空气中但与水分子无关的颗粒造成的。这种回波对于业务需求没有意义，因此将其视为无效降水数据。

(2)考虑到距离衰减对雷达回波的作用，当距离越远时雷达探测能力越弱。所以，如果最小阈值取值过小，就有可能将离雷达较远处的一些有用的降水回波剔除掉。综合上述因素，本实验取最小阈值为15dBZ。同时，由于冰水混合物相比于水滴的反射能力更强，当降水回波中有冰水混合物时，会造成反射率因子显著增大。故本实验对降水回波强度大于55dBZ的值，简单的以55dBZ代替，以防大范围降水过程中有少数小冰粒进入从而污染初值场，即全部初始赋值。

在实际应用中，首先将多普勒天气雷达回波强度产品数据所组成的基本产品数据dBZ，利用下述公式计算获得

上式中，

Z₀为晴空时的基本反射率；

Z为反射率因子；

然后，利用公式Z＝200I^1.6，将Z转化为降水强度I。

步骤二、自动气象站雨量数据质量控制。

质量控制时，对非常规的降水情况予以排除，包括如下2个方面：

首先，小时雨量超过1000mm的数据不参与统计；

其次，由于每小时小于5mm的降水属于弱降水，同时自动气象站雨量计本身误差较大，如果以此作为统计依据则参考性较差，所以每小时雨量小于5mm的数据本实验也不参与统计。

步骤三、利用几种方法计算出降水数据，并分别与实测值作对比。

分别运用传统的卡尔曼滤波法和本发明方法，对文中所述35个时次的多普勒天气雷达降水资料进行估算，并与同一降水场次相应自动气象站雨量计记录的实测值进行对比。为考查本发明方法估算降水的效果，本实验利用绝对误差、相对误差、估算结果作为降水估算精度的比较因子，并列于表1中。

表1

由表1可知，传统的卡尔曼滤波法除在2013年3月30日的降水估算量比实测值大以外，其他时段内估算得到的降水量与同时段内相应站点自动气象站雨量计观测得到的实测值相比都出现不同程度偏小的情况。同时，传统的卡尔曼滤波法得到的平均估算结果也比实测值要小。而本发明的卡尔曼滤波改进方法估算得到的各场次估算结果都比实测值稍大，同时平均估算结果也比实测值要高。可见，在算法性能稳定性方面，本发明的卡尔曼滤波改进方法比传统的卡尔曼滤波法更加稳定。

从两种算法平均估算结果的值来看，传统的卡尔曼滤波法的平均估算结果与实测值相比相对误差为2.5％。本发明的卡尔曼滤波改进方法的平均估算结果与实测值相比相对误差为2.4％。从图3、图4中可以看出，传统的卡尔曼滤波和自动气象站雨量计获得的时雨量相关性，与改进后的卡尔曼滤波和自动气象站雨量计获得的时雨量相关性两者的R²分别是0.9797和0.9924，差异效果比较明显。同时，由图中可以看出，改进后的卡尔曼滤波法与传统的卡尔曼滤波法的估算效果相比，点距比较集中，拟合度高，相关性更强。这说明了本发明的卡尔曼滤波改进方法相比传统的卡尔曼滤波法估算精度更高。

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