1、卡尔曼滤波方法应用南京信息工程大学气象台一、滤波的气象意义二、卡尔曼滤波方法三、递推滤波系统的参数计算方法四、递推系统制作预报的业务流程五、应用中的若干问题讨论六、应用步骤Rudolf Emil Kalman 匈牙利数学家 BS&MS at MIT(Massachusetts Institute of Technology)PhD at Columbia 1960年发表的论文A New Approach to Linear Filtering and Prediction Problems（线性滤波与预测问题的新方法）Signal Processing 数字滤波：通过一种算法排除可能的随机干扰
2、，提高检测精度的一种手段 线性系统 f(A+B)=f(A)+f(B)数学方法处理 噪声信号输入-尽可能少噪声输出 数值预报产品的释用技术方法：1、人的经验为主的定性方法天气 学方法）。2 2、客观定量方法（统计学方法、动 力释用方法、神经元网络）。MOS方法是被广泛释用的数值产品方法，是以数值产品历史资料为基础建立MOS方程的，资料年限太短（不足一年），方程统计特性差，资料年限长（2-3年），方程统计特性好，但在积累资料及用MOS方程作预报时不能改进及更新模式。在数值预报迅速发展的今天显然是不可能的。MOS方法示意图解决途径如下：1、根据新模式的统计特征，对MOS方程进行订正。2、用新模式重新
3、对2-3年的历史样本进行计算，以积累数值产品历史资料。3、只需少量的数值产品历史资料，建立能适应数 值模式变化的统计模型，这种方法越来越得到 广大气象工作者的重视，卡尔曼滤波方法就具 有这种特点。目前，我国数值预报发展迅速，数值模式更新快，广大台站积累足够供建立MOS方程使用的数值产品历史资料比较困难，因此，卡尔曼滤波方法在我国天气预报中有广泛的应用前景。不需要保存全部历史资料数据，可借助于前时刻的滤波结果，递推出现时刻的状态估计量，大大减少了存储量和计算量。一般为具有线性变化特征的连续性变量。飞行 潜艇导航 导弹弹道计算 (1969年的APPOLO)气象业务预报(1987年)(应用成功的主要
4、是北欧国家，如芬兰、瑞士、丹麦等)最高最低气温预报。(1992年日本制作56个站)Application 视频跟踪 一、滤波的气象意义 在实际问题中，常常遇到所获得的信息混杂着其它噪音，希望排除无用的干扰而能最佳估计出有用的信息，滤波是处理这类实际问题的重要方法。预报员每天用各种方法制作天气要素预报，可以得到带有误差的预报值时间序列，造成预报误差的原因很多，我们试图订正它。根据滤波的基本思想，卡尔曼滤波可以用于处理一系列带有误差的预报值而得到它的最佳估算值，这对提高预报精度具有重要的现实意义。卡尔曼滤波方法通过利用前一时刻卡尔曼滤波方法通过利用前一时刻预报误差反馈到原来的预报方程，及预报误差反
5、馈到原来的预报方程，及时修正预报方程系数，以此提高下一时修正预报方程系数，以此提高下一时刻的预报精度，这是卡尔曼滤波方时刻的预报精度，这是卡尔曼滤波方法用于天气预报的气象意义。法用于天气预报的气象意义。而MOS方程一旦建立之后，在制作预报过程中，预报误差不能反馈到MOS方程中，更不能修正方程系数，这就是这两种方法的重要区别之一。预测值(假设当前值仅跟上一时刻有关)但变化中可能有噪声 观测值 读数会有误差 两种噪声相互无关 根据连续的观测值来推算实际变量值的变化 KF是根据上一状态的估计值和当前状态的观测值推出当前状态的估计值的滤波方法S(t)=f(S(t-1),O(t)它是用状态方程和递推方法
6、进行估计的，因而卡尔曼滤波对信号的平稳性和时不变性不做要求维纳滤波：使用全部观测值保证平稳性卡尔曼滤波方法示意图卡尔曼滤波方法示意图 二、卡尔曼滤波方法二、卡尔曼滤波方法 递推滤波可用于解决如何利用前一前一时刻预报误差时刻预报误差来及时修正修正预报方程系系数数这一问题。滤波对象假定是离散时离散时间线性间线性动态系统，并认为天气预报对象是具有这种特征的动态系统，可用以下两组方程来描述：Y Yt t=X=Xt tt t+t t （1 1）t t=t-1t-1+t-1t-1 （2 2）(1)(1)式为预报方程，t t为量测噪声，是n n维随机向量；Yt是n n维 量 测 变 量（预 报 量），可 用
7、 下 式 表 示：Yt=y y1 1,y,y2 2,y,yn n tT，X Xt t是n nm m维的预报因子矩阵，t t是m m维回归系数。在递推滤波方法中，将t t作为状态向量，它是变化的，用状态方程(2)(2)式来描述其变化。(2)(2)式中t-1t-1是动态噪声。ttnmnnmmtTtntTVpRHIvTTxxxxxxxxxXyyyY8501000500850850100021222211121121),(月降水月均北极涡Tt，风）温度，月平均气温，(动态噪声t-1t-1与量测噪声t t都是随机向量，并假定二者互不相关、均值为零、方差分别为W W 和V V 的白噪声。通常用 Y Yt
8、t=X=Xt tt t+t t (1)(1)t t=t-1t-1+t-1t-1 (2)(2)两方程来描述离散时间的线性动态系统。具有这种特征的天气预报对象所关心的是它的状态向关心的是它的状态向量的变化量的变化。根据上述对t-1t-1和t t的假定，运用广义最小二乘法，可以得到一组递推滤波公式，这一组公式组成了递推滤波系统。1.1.运用运用t-1t-1时刻预测时刻预测t t Y Yt=X Xtt-1 2.t-12.t-1时刻最优估计值误差方差时刻最优估计值误差方差 C Ct-1,t,t时刻预测值方差时刻预测值方差 R Rt R Rt=C=Ct-1+W W 3.Kalman3.Kalman增益增益
9、 KG=KG=R Rt X XtT/(/(X Xt R Rt X XtT+V)=+V)=At 4.4.估计估计t t时刻系数时刻系数t t=t-1+A+At(Y Yt)A At=R Rt X XtT t-1 t=X Xt R Rt X XtT+V5.Kalman5.Kalman增益更新得增益更新得t t时刻最优估计值误差方差时刻最优估计值误差方差 C Ct C Ct=(I-KG=(I-KG Xt)R Rt=R Rt-A At tA AtTx Y Yt=X Xtt-1 R Rt=C=Ct-1+W W t=X Xt R Rt X XtT A At=R Rt X XtT t-1 t=t-1 +A A
10、t(Y Yt)C Ct=R Rt-A At tA AtT上述六个公式组成的递推滤波系统体现了卡尔曼滤波的基本思想。每加进一次新的量测(，X Xt),只需利用已算出的前一次滤波值t-1 和滤波误差方差阵C Ct-1，便可算出新的状态滤波值t和新的滤波误差方差阵C Ct，就能通过公式得到t+1时刻的预报值。这样不论预报次数如何增加，不需要存储大量历史的量测数据，大大减少了计算机的存贮，而且只进行矩阵的加、减、乘和求逆，通常计算量不大，从而满足了应用滤波的实时性要求。这就是卡尔曼滤波方法的优点。三、递推滤波系统的参数计算方法三、递推滤波系统的参数计算方法 分析上面的一组递推公式可以得知，t t，C
11、Ct，W W，V V是重要参数，在确定这四个参数的基础上，利用数值模式提供的预报因子X Xt t、前一次预报量及其观测值，才能通过更新预报方程系数制作预报，因此，必须研究这四个参数的计算方法。1 1、递推系统递推系统参数初值参数初值的计算方法的计算方法 要反复运算上述六个公式来实现递推过程，必须首先确定初值0 0，C C0。我们通常采用以下客观方法：0 0 的确定。C C0的确定。2 2、递推系统参数递推系统参数W，V的计算方法的计算方法 W W、V V 分别是动态噪声和量测噪声的方差阵，可以假定随机扰动的特性不随时间变化，但是，必须在应用上述递推系统之前确定。W的确定：根据白噪音的假定，W的
12、非对 角线元素均为零。mwwwW0000021 30)3.0(00030)6.0(00030)6.0()(000)(000)(222232221TTTW可以用可以用的变化来估算的变化来估算W W 值值 V V的确定。根据白噪音的假定，V的非对角线元素均为零：nVVVV0000021 利用样本资料对预报量Y的n分量(y y1 1,y y2 2,y,yn n)建立回归方程后，可以求出n个残差(q q1,q,q2 2,q,qn n)，从回归分析得知：21)(itktitiyyq q q1 1/(k-m-1),q(k-m-1),q2 2/(k-m-1),.q(k-m-1),.qn n/(k-m-1)(
13、k-m-1)分别为v v1 1,v,v2 2,v,vn n的无偏估计值，其中k是样本容量，m是因子个数，必须km+1，因此有:1000100121mkqmkqmkqVn 我们只要用少量（2个月）的量测（X Xt，Y Yt）样本资料，就能得到这四个递推系统参数0，C C0，W，V。3、递推过程中的参数计算方法 系数的更新原理是在已知前一时刻（t-1）的系数t-1的基础上加上订正项，获取订正项构成了递推的主要过程。除了预报误差对方程系数更新有重要影响外，预报因子质量也是最重要的因素之一。应用递推系统的过程是每增加一次新的量测Xt和Yt时，利用W、V，前一次的系数t-1及其误差C Ct-1就可推算下
14、一时刻的t和C Ct，同时又作了要素预报，如此反复循环进行。四、递推系统制作预报的业务流程四、递推系统制作预报的业务流程 卡尔曼滤波系统适用于制作温度、湿度和风等连续性预报量的预报，为预报员提供这类 客观指导预报产品。因此，在计算机上建立了用递推方法自动制作上述预报的自动化业务流程。1、流程结构、流程结构 数值产品格点值、站点的天气要素观测值 预报量的预报值文件、预报方程系数文 件、预报方程系数误差的方差文件 （这些文件的内容随递推过程不断更新）流程 输入两类基本数据文件，就可依次递推系统中各个参数，得到系统本身生成的数据文件，作为下一时刻运行递推系统的输入信息。随机误差方差在递推起始被确定后
15、，不再随递推过程改变。2、流程特点、流程特点 整个流程计算量不大，存储空间也小，一般气象台站配备微机就能应用。五、应用中的若干问题讨论五、应用中的若干问题讨论1、预报对象的选择预报对象的选择 预报对象最好选择具有线性变化特征的连续性变量，如温度、湿度、风等。2、预报因子的选择预报因子的选择 预报因子与预报对象之间相关程度高而且预报因子要具有较高的精度，预报因子的个数不宜过多，一般不超过4个。3、递推滤波的时间间隔、递推滤波的时间间隔 递推滤波的时间间隔不宜长，一般在短时或短期预报中应用卡尔曼滤波方法优于中期预报。4、预报精度、预报精度 选择好的预报因子是至关重要的。5、预报滞后现象、预报滞后现
16、象 预报值的变化滞后于观测实况的变化，尤其在预报对象发生剧烈变化时比较明显，要克服这一现象有待进一步研究。北京地区北京地区19891989年年1111月月1212月月9595年年1212月至月至9696年年4 4月侯平均温度实况与预报对照图月侯平均温度实况与预报对照图 a.a.北京北京 b.b.福州福州 c.c.成都成都 六、应用步骤六、应用步骤初选因子 选取一些与预报量相关高，影响大的因子 收集资料 收集所选因子的两个月的数值产品资料 收集预报量的站点实测资料 因子筛选 建立逐步回归方程，筛选因子，因子数最好不超过4个 确定递推系统参数0 C C0 W V 利用所收集的两个月的数值产品资料建立 两个回归方程，求取回归系数，即0、T C C0是0的误差方差阵，C C0取为零方阵 wj=(j)T-(j)02/T j=1,2,m vj=q qj j/(k-m-1)/(k-m-1)将上述参数、资料输入计算流程结束！卡尔曼滤波系统制作天气要素预报的业务流程卡尔曼滤波系统制作天气要素预报的业务流程 建立数据文件建立数据文件卡尔曼滤波系统制作天气要素预报的业务流程卡尔曼滤波系统制作天气要素预报的业务流程 递推系统计算流程递推系统计算流程返回返回
