最佳估计理论之卡尔曼滤波

1、管理系统分析专题报告1/30最佳状态估计理论卡尔曼滤波方法研究摘要：在任何系统中，为了对系统形成有效的控制，系统状态的准确把握显得尤为重要。我们可以通过一系列的手段对特定系统进行观测，以估计出系统的过去、现在、未来的状态，具体应用可分别表现为对过去状况的评估、当前状态的实时控制、趋势的准确预测等等。最优估计出系统状态过程中，实际的量测往往是存在诸如来自系统自身、测量工具等所带来的干扰，控制论中将这种干扰定义为噪声。如何去寻求滤除这种噪声干扰，便成为最佳系统状态估计首先必须解决的问题。Kalman滤波等一些滤波算法便因此应运而生，其作为一种最优估计理论与方法，由于它的实时递推、存储量小和简单易行

2、的特点，在工程应用中受到了重视，广泛应用于信号处理、控制、通信、航天、制导、目标跟踪、石油勘探、故障诊断、卫星测控、GPS定位、检测与估计及机器人等等领域。卡尔曼滤波随时间及研究的发展，已形成了多种多样的理论和应用的形式。本次的学习带着了解认识该滤波算法思想和数学思维的目的，只对一般卡尔曼滤波问题（最优预测与最优滤波）进行了基本的研究。据此本文主要内容安排包括：1、卡尔曼滤波的理论背景；2、理论基础（选取离散系统的最优预测问题作为代表）阐述；3、工程扩展应用状况；4、理论局限性等方面。关键字：卡尔曼滤波；最优估计；状态估计；控制2/30目录1绪论 . 41.1卡尔曼滤波理论研究及应用概述 .

3、41.2维纳滤波简述 . 51.卡尔曼滤波理论概述 . 72卡尔曼滤波理论基础 . 92.1卡尔曼滤波问题的提法 . 92.2离散系统卡尔曼最优预测基本方程的推导（举例） . 122.2.1求解基本过程推演. 132.2.2卡尔曼预测估计递推方程的计算步骤 . 192.2.3应用算例 . 203工程扩展应用举例 . 223.1卡尔曼滤波在飞机控制中的应用 . 223.2基于卡尔曼滤波方法的时用水量预测（定量描述略） . 264卡尔曼滤波局限性分析 . 284.1稳定性定理 . 284.2滤波的发散问题 . 284.3卡尔曼滤波理论的进一步发展 . 29参考文献. 303/301绪论1.1卡尔曼

4、滤波理论研究及应用概述在自动控制、航空与航天、通讯、导航和工业生产等领域中，越来越多地遇到“估计”问题。所谓“估计”，简单地说，就是从观测数据中提取信息。在飞行器导航中，要从带有随机干扰的观测数据中，估计出飞行器的位置，速度和加速度等运动状态变量，这就遇到状态变量的估计问题，这些变量都是随机过程。因此，“估计”的任务就是从带有随机误差的观测数据中估计出某些参数或某些状态变量，这些被估计参数或状态变量可统称为被称为被估量。我们希望估计出来的参数或状态变量越接近实际值越好。我们希望估计是最好的，因此提出最优估计问题。所谓最优估计，是指在某一确定的准则条件下，按照某种统计意义上来说，估计达到最优。为

5、了正确地解决参数估计和状态估计问题，首先要研究估计方法。最早的估计方法是高斯（Gauss,K.F.）于1795年在他的关于天体运动理论一书中提出的最小二乘法。最小二乘法没有考虑被估参数和观测数据的统计特性，因此这种方法不是最优估计方法。在192年费歇（Fisher,R.A.）提出了极大似然估计方法：从概率密度出发来考虑估计问题，对估计理论做出了重大贡献。对于随机过程的估计问题，到本世纪三十年代才积极开展起来。主要成果为1940年美国学者维纳（Wiener,N）所提出的在频域中设计统计最优滤波器的方法，这一方法称为维纳滤波。同时期，苏联学者哥尔莫郭洛夫（A.H.）提出并初次解决了离散平稳随机序列

6、的观测和外推问题。维纳滤波和哥尔莫郭洛夫滤波方法，局限于处理平稳随机过程，并只能提供稳态的最优估值。这些滤波方法在工程实践上也遇到许多困难，在实际应用上受到一定的限制。1960年卡尔曼（Kalman,R.E.）和布西（Bucy,R.S）提出了最优递推滤波方法，可用数字计算机来实现。卡尔曼滤波既适用于平稳随机过程，又适用于非平稳随机过程，因此卡尔曼滤波方法得到广泛的应用。本文要学习的卡尔曼滤波理论，正是源于他的博士论文和1960年发表的论文A New Approach to Linear Filtering and Prediction Problems（线性滤波与预测问题的新方法）。卡尔曼滤波

7、在工程实践中，特别在航空空间技术中迅速得到应用。例如在测4/30轨问题和惯性导航等方面都应用卡尔曼滤波。所谓测轨问题，就是在不同时刻，对飞行器进行观测，根据观测数据，应用卡尔曼滤波方法，估计出这个飞行器每时每刻的状态变量，如飞行器的位置、速度、加速度及阻力系数等物理量。以便对飞行器进行导航、制导和拦截。斯坦利.施密特(Stanley Schmidt)首次实现了卡尔曼滤波器。卡尔曼在NASA埃姆斯研究中心访问时，发现他的方法对于解决阿波罗计划的轨道预测很有用，后来阿波罗飞船的导航电脑便使用了这种滤波器。关于这种滤波器的论文由Swerling(1958)、Kalman(1960)与Kalman a

8、nd Bucy(1961)发表。在其它方面卡尔曼滤波也得到应用。目前,卡尔曼滤波已经有很多不同的实现。卡尔曼最初提出的形式现在一般称为简单卡尔曼滤波器。除此以外，还有施密特扩展滤波器、信息滤波器以及很多Bierman，Thornton开发的平方根滤波器的变种。也许最常见的卡尔曼滤波器是锁相环，它在收音机、计算机和几乎任何视频或通讯设备中广泛存在。当前对于卡尔曼滤波较前沿的研究与应用有很多，诸如机器人导航、控制；传感器数据融合；雷达系统以及导弹追踪；计算机图像处理，包括头脸识别、图像分割、图像边缘检测等。1.2维纳滤波简述为了实现最优估计，最佳线性滤波理论起源于40年代美国科学家Wiener和前

9、苏联科学家等人的研究工作，后人统称为维纳滤波理论。这是在信号和干扰都表示成有理谱密度的情况下，找出最优滤波器，使得实际输出与希望输出之间的均方误差最小。维纳滤波问题的关键是导出维纳霍夫积分方程，解这一积分方程可得最优滤波器的脉冲过渡函数，从脉冲过渡函数可得滤波器的传递函数。通常，解维纳霍夫积分方程是很困难的，即使对少数情况，能得解析解，但在工程上也往往难以实现。特别对于有限记忆或非常平稳过程，维纳滤波问题变得更为复杂了。维纳滤波问题的提法简单罗列如下：设系统的观测方程为z(t)=x(t)+v(t)；式中，x(t)为有用信号z(t)为观测信号5/30v(t)为观测误差。设x(t)、z(t)和v(

10、t)都是均值为零并具有各态历经性的平稳随机过程。(t)。维纳滤波的任务就是根据观测值z(t)估计x(t)，使x(t)估值接近于x设计出一个线性定常系统L，如图2-1所示，使得系统的输出y(t)与x(t)具有最小方差，即J=Ex(t)-y(t)=min 2 (1-3)(t)。 这样y(t)就作为x(t)的估值x如果系统的脉冲过渡函数为h(l)，则y(t)=0h(l)z(t-l)dl (1-4)y(t)是系统L根据输入信号z(t)在(-，t)上的全部过去值所给出的实际输出，如图2-2所示。y(t)是z(t-l)的线性函数(l0)。根据问题的性质，维纳滤波有下列三个条件：信号与噪声都是均值为零并具有

11、各态历经性的平稳随机过程；滤波器是一个物理可实现的线性定常系统。当lt称为预测（或外推）问题第二类：t1=t，称为滤波问题；第三类：t1k，称为预测（或外推）问题；第二类：j=k，称为滤波问题；第三类：jk称为平滑（或内插）问题。本文这里只讨论连续系统和离散系统的最优预测和最优滤波问题。并通过简述离散系统卡尔曼滤波最优预测基本方程的推导来说明卡尔曼滤波的具体过程，其余类型过程基本相似，在此只描述它们的推演处理方法及思想，详细过程略去。2.2离散系统卡尔曼最优预测基本方程的推导（举例）在推导卡尔曼预测基本方程时，为了简便起见，先不考虑控制信号的作用，这样，离散系统的差分方程(2-11)变成x(k

12、+1)=F(k+1,k)x(k)+G(k+1,k)w(k)观测方程仍为式(2-2)： (2-16)z(k)=H(k)x(k)+v(k)式中，w(k)和v(k)都均值为零的白噪声序列，w(k)和v(k)相互独立。 在采样间隔内w(k)和v(k)为常值，其统计特性如式(2-13)所示，即Ew(k)=Ev(k)=0Ew(k)wT(j)=QdkkjEv(k)vT(j)=RdkkjEw(k)vT(j)=0 状态向量的初值x(0)，其统计特性是给定的，即Ex(0)=m0TEx(0)-mxt-m=P0 000()给出观测序列z(0)，z(1)，L，z(k)，要求找出x(k+1)的线性最优估计(k+1/k)，

13、使得估计误差x%(k+1/k)=x(k+1)-x(k+1/k)的方差为最小，x即12/30T(k+!/k)Exk+1-xk+1/k()()x(k+1)-x=min (k+1/k)是z(0)，z(1)，要求估值xL，z(k)的线性函数，并且要求估计是无偏的，即(k+1/k)=EExx(k+1)2.2.1求解基本过程推演下面推导卡尔曼预测基本公式。推导的方法有几种，比较简易的方法是利用正交定理，用数学归纳法推导卡尔曼估计的基本递推估计公式。当获得观测值z(0)，z(1)，L，z(k-1)之后，假定已经得到状态向量x(k)的(k/k-1)。当还未获得k时刻的新观测值z(k)时，根据一个最优线性预测估

14、计x已有的观测值，可得k+1时刻的系统状态向量x(k+1)的两步预测估值(k+1/k-1)： x(k+1/k-1)=F(k+1,k)x(k/k-1) x (2-17)(k/k-1)是x(k)的一步最优线性估计，x(k+1/k-1)也是x(k+1)的由于x最优线性预测估计，这可用正交定理来证明。由式(2-16)减式(2-17)，可得%(k+1/k-1)=x(k+1)-x(k+1/k-1)x(k/k-1)=F(k+1,k)x(k)-x+G(k+1,k)w(k)%(k+1/k-1)=F(k+1,k)x%(k/k-1)+G(k+1,k)w(k) x%(k/k-1)=x(k)-x(k/k-1) 式中，x

15、 (2-18)(k/k-1)是x(k)的最优线性预测估值，%(k/k-1)因为x根据正交定理，估计误差x%(k/k-1)的线性变换必须正交于z(0)，z(1)，L，z(k-1)，所以x%(k/k-1)也必须正交于z(0)，F(k+1,k)xL，z(k-1)。式(2-18)中的w(k)是均值为零的白噪声序列，与z(0)，L，z(k-1)相互独立，因此w(k)正交于(k+1/k-1)是x(k+1)的最优两z(0)，L，z(k-1)。所以在没有获得之前z(k)，x步线性预测。13/30(k+1/k-1)来得到x(k+1)的在新观测值z(k)获取之后，可通过修正两步估值x一步预测估值x(k+1/k)。

16、设z(k)的预测估值(k/k-1) z(k/k-1)=H(k)x由式(2-2) (2-19)z(k)=H(k)x(k)+v(k)减式(2-19)，得z(k)的预测估计误差为%(k/k-1)=z(k)-z(k/k-1)z(k/k-1)=H(k)x(k)-x+v(k)%(k/k-1)+v(k)=H(k)x%(k/k-1)的原因有两个：造成z对k时刻状态向量x(k)的预测估计有误差；附加了v(k)时刻的观测噪声干扰。%(k/k-1)包含修正x(k+1/k-1)的新信息。这样，在获得之后，在两步估显然z(k+1/k-1)的基础上，用z%(k/k-1)去修正，便可得到k+1时刻状态向量值x(k+1/k)

17、，即 x(k+1)的一步预测估计x(k+1/k)=x(k+1/k-1)+K(k)z%(k/k-1) x或(k+1/k)=F(k+1,k)x(k/k-1)+K(k)(k/k-1)xz(k)-H(k)x(2-20) 式中K(k)是待定矩阵，称为最优增益矩阵或加权矩阵。 式(2-20)可改写成(k+1/k)=F(k+1,k)x(k/k-1)x%(k/k-1)+K(k)v(k)+K(k)H(k)x (2-21)k+1时刻系统状态方程为x(k+1)=F(k+1,k)x(k)+G(k+1,k)w(k) (2-22)由式(2-22)减式(2-21)得x(k+1)的估计误差为14/30%(k+1/k)=%(k

18、/k-1)xF(k+1/k)-K(k)H(k)x+G(k+1/k)w(k)-K(k)v(k) (2-23) %(k/k-1)、w(k)、v(k)均分别正交于观察式(2-23)右边，x%(k+1/k)()正交于z(0)，z(1)，z(0)，z(1)，L，z(k-1)，因此，xL，z(k-1)。%(k+1/k)与z(k)正交，则x(k+1/k)就是x(k+1)的一步最优线性预测估值。若x%(k+1/k)与z(k)的正交条件： 因此，利用x%(k+1/k)zT(k)=0 Ex (2-24)来确定最优增益矩阵K(k)。把式(2-23)和式(2-2)代入式(2-24)，得%(k/k-1)+G(k+1/k

19、)w(k)-K(k)v(k)E【F(k+1/k)-K(k)H(k)xT(k/k-1)+H(k)x%(k/k-1)+v(k)】H(k)x=0(k/k-1)、w(k)和v(k)相互间都是正交的，因此上式可简考虑到(k/k-1)、x化为%(k/k-1)x%T(k/k-1)HT(k)-K(k)v(k)vT(k)=0EF(k+1/k)-K(k)H(k)x即Tk/k-1HTk-KkEvkvTk=0x%Fk+1/k-KkHkEk/k-1x()()()()()()()()()%(k/k-1)x%T(k/k-1)=P(k/k-1)，又有Ev(k)vT(k)=R，代入上式后设Exk可得最优增益矩阵为-1TT K(

20、k)=F(k+1,k)P(k/k-1)H(k)H(k)P(k/k-1)H(k)+R(k)(2-25)现在需确定估计误差方差阵P(k+1/k)的递推关系式。 根据估计误差方差阵的定义，有%(k+1/k)x%T(k+1/k) P(k+1/k)=Ex将式(2-23)代入上式得15/30%T(k/k-1)P(k+1/k)=E【F(k+1,k)-K(k)H(k)xT%(k/k-1)+G(k+1,k)w(k)-K(k)v(k)Fk+1,k-KkHk()()()xT+G(k+1/k)w(k)-K(k)v(k)】%(k/k-1)、w(k)和v(k)相互间都正交，可得： 考虑到xP(k+1/k)=F(k+1,k

21、)-K(k)H(k)P(k/k-1)F(k+1,k)-K(k)H(k)T+K(k)RKT(k)+G(k+1,k)QGT(k+1/k)kk式中Q=Ew(k)wT(k) k将上式展开，经整理后得P(k+1/k)=F(k+1,k)P(k/k-1)FT(k+1/k)-K(k)H(k)P(k/k-1)FT(k+1/k)-F(k+1,k)P(k/k-1)HT(k)KT(k)+K(k)H(k)P(k/k-1)HT(k)KT(k)+K(k)RKT(k)k+G(k+1,k)QGT(k+1/k)k上式右边第四项与第五项之和为 (2-26)K(k)H(k)P(k/k-1)HT(k)KT(k)+K(k)RKT(k)k

22、=K(k)H(k)P(k/k-1)HT(k)+RKT(k)k-1TT=F(k+1,k)P(k/k-1)H(k)H(k)P(k/k-1)H(k)+RkH(k)P(k/k-1)HT(k)+RKT(k)k=Fk+1,kPk/k-1HTkKTk()()()()显然，式(2-26)右边第四项与第五项之和在数值上等于第三项，但符号相反。这 样，式(2-25)右边第三、第四和第五项之和为零，所以P(k-1/k)=F(k+1/k)P(k/k-1)FT(k+1,k)-K(k)H(k)P(k/k-1)FT(k+1/k)+G(k+1/k)QGT(k+1/k)k16/30把式(2-25)的K(k)代入上式，可得估计误

23、差方差阵的递推关系式为P(k+1/k)=F(k+1,k)P(k/k-1)FT(k+1/k)-1TT(2-F(k+1,k)P(k/k-1)H(k)H(k)P(k/k-1)H(k)-RkH(k)P(k/k-1)FT(k+1,k)+G(k+1,k)QGT(k+1,k)k27)方程(2-20)、(2-25)和(2-27)构成一组完整的最优线性估计方程，现综合如下： 最优预测估计方程为式(2-20)，即(k+1/k)=F(k+1,k)x(k/k-1)+K(k)(k/k-1)xz(k)-H(k)x最优增益矩阵方程为(2-25)，即-1TT K(k)=F(k+1,k)P(k/k-1)H(k)H(k)P(k/

24、k-1)H(k)+Rk估计误差方差阵的递推方程为式(2-27)，即P(k+1/k)=F(k+1,k)P(k/k-1)FT(k+1/k)-1TT -F(k+1,k)P(k/k-1)H(k)H(k)P(k/k-1)H(k)-RkH(k)P(k/k-1)FT(k+1,k)+G(k+1,k)QGT(k+1,k)k从式(2-27)可看出，估计误差方差阵与Q和R有关，而与观测值z(k)无关。因kk此，可事先估计误差方差阵P(k+1/k)，同时也可算出最优增益矩阵K(k)。 按式(2-16)和式(2-2)作出系统模型方块图，如图2-3所示。图2-4表示由方程(2-21)所描述的卡尔曼最优预测估计方块图。17

25、/30从图2-4可看出，最优预测估计由三部分组成：估值为观测值的线性函数；最优增益矩K(k)；单位负反馈回路。现在需要验证第二节时最优估计提出的三项要求：估值为观测值的线性函数；估计值的方差为最小；估值是无偏的，即(k+1/k)Ex=Ex(k+1)在上面推导预测估计方程时，是按照和两项要求推导的，现需要说明估值是无偏的问题。对式(2-16)的两端取数学期望，考虑到Ew(k)=0，可得Ex(k+1)=F(k+1,k)Ex(k) (2-28)再对式(2-21)的两端取数学期望，考虑到Ev(k)=0，可得(k+1/k)(k/k-1)%(k/k-1)Ex=F(k+1,k)Ex+K(k)H(k)Ex(k

26、/k-1)(k/k-1)=F(k+1,k)Ex+K(k)H(k)Ex(k)-x将式(2-28)减式(2-29)，得到 (2-29)(k+1/k)(k/k-1)Ex(k+1)-x=F(k+1,k)Ex(k)-x(k/k-1)-K(k)H(k)Ex(k)-x如果初始条件为 (2-30)(0/0)Ex=Ex(0)=m0或18/300/0x则 (-)=Ex(0)=m0()%0/0=E0/0=0 Exx(0)-Ex-根据式(2-30)的递推关系，可得(1/0)Ex(1)-x=0(2/1)Ex(2)-x=0M(k+1/k)Ex(k+1)-x=0()因此(k+1/k)Ex=Ex(k+1)k=0,1,2L0/

27、0=E所以，只要初始估计选为Exx(0)，所得估计是无偏的。 -2.2.2卡尔曼预测估计递推方程的计算步骤卡尔曼预测估计递推方程的计算步骤如下：在t时刻给定初值： 0估值误差方差阵初值：P0/0()0/0x(-(0)=E)=xx(0)=m0 (T=P(0)=Ex0-x0x0-x0()()()() -根据公式(2-25)计算t时刻最优增益阵K(0)： 0-1TT K(0)=F(1，0)P(0)H(0)H(0)P(0)H(0)+R0(1/0)： 根据公式(2-20)计算x(1)的最优估值x(1/0)=F(1，(0)+K(0)(0)x0)xz(0)-H(0)x根据公式(2-27)计算P(1/0)：-

28、1TTTP(1/0)=F(1，0)P(0)F(1，0)-F(1，0)P(0)H(0)H(0)P(0)H(0)+R0H(0)P(0)FT(1，0)+G(1，0)QG(1，0)019/30根据已知的P(1/0)计算t时刻的K(1)。 1(2/1)。 根据K(1)计算x(2)的估值x重复上述递推计算步骤，可得(3/2)，(k/k-1)，P(k/k-1),K(k)，x(k+1/k) P(2/1),K(2),xL，x2.2.3应用算例例2-1设系统状态方程和观测方程为x(k+1)=0.5x(k)+w(t)z(k)=x(k)+v(k)w(k)和v(k)都是均值为零的白噪声序列，且不相关，其统计特性如下：E

29、w(k)=0,Ev(k)=0Ew(k)w(j)=1d初值Ex(0)=m0=0,kj,Ev(k)v(j)=2d kjP0/0(-)=1观测值z(0)=0，z(1)=4，z(2)=2试求x(k)的最优预测估值。解:与前面的有关方程对照，可得F(k+1,k)=0.5，G(k+1,k)=1H(k)=1，Q(k)=1，Rk=2k最优预测估值方程为(k+1/k)=0.5x(k/k-1)+K(k)(k/k-1)xz(k)-x最优增益矩阵为K(k)=0.5P(k/k-1)P(k/k-1)+2估值误差方差阵递推方程为 -12-1P(k+1/k)=0.25P(k/k-1)-0.25Pk/k-1Pk/k-1+2()

30、()+1(1/0)，x(2/1)和x(3/2) 下面计算x20/30P0/0(-)=116K(0)=0.51(1+2)-1=1P(1/0)=0.25-0.25+1=1.1673K1=0.51.1671.167+2-1=0.1842 ()()P(2/1)=0.251.167-0.251.167(1.167+2)-1+1=1.1842K(2)=0.51.1842(1.1842+2)-1=0.18590/0取x(-)的初值0/0x(-)=Ex(0)=m0=01(1/0)=0.50+0=0x6(2/1)=0.50+0.1842(4-0)=0.7368x(3/2)=0.50.7368+0.1859(2-

31、0.7368)=0.6032x2.2.4其它系统估计问题基本方程1.连续系统卡尔曼最优滤波2.离散系统卡尔曼最优滤波3.噪声为有色噪声情况4.系统噪声与观测噪声相关情况与离散系统卡尔曼最优预测推导过程和基本思想相似，最终得出卡尔曼滤波的三个重要方程。由于时间限制，这里不再详细阐述。21/303工程扩展应用举例3.1卡尔曼滤波在飞机控制中的应用卡尔曼滤波在航空和航天中，以及在其他工程领域中都得到应用。对于每一个具体应用问题，都要求深入了解问题的物理实质和工程实际问题，因此比较复杂。这里仅对卡尔曼滤波在飞机自动驾驶中的经典应用进行了探讨。图3-1自动驾驶仪（附注：自动驾驶仪是按技术要求自动控制飞行

32、器轨迹的调节设备，其作用主要是保持飞机姿态和辅助驾驶员操纵飞机。对无人驾驶飞机，它将与其他导航设备配合完成规定的飞行任务。导弹上的自动驾驶仪起稳定导弹姿态的作用，故称导弹姿态控制系统。自动驾驶仪是模仿驾驶员的动作驾驶飞机的。它由敏感元件、计算机和伺服机构组成。当飞机偏离原有姿态时，敏感元件检测变化，计算机算出修正舵偏量，伺服机构将舵面操纵到所需位置。）22/30图3-2飞机纵向运动参量关系图设飞机在垂直平面内运动（纵向运动）各参量的关系如图3-2所示。 设O为飞机重心，OX1为飞机纵轴，v为飞机速度向量。OX1轴与地平面的夹角J称为俯仰角。速度向量v与地平面的夹角q称为航迹角。a称为攻角，d为

33、飞机操纵面的偏转角。J、q和a成下列关系： rrJ=q+a(3-1)飞机的纵向运动方程为&+aJ&+aa=ad(3-2) J123&-a&+a4a=0(3-3) J飞机在飞行过程中往往受到扰动气流的作用，可把这一扰动作用归结为在攻角a中出现了附加的随机干扰项w1，则方程(3-2)和式(3-3)变成&+aJ&J1+a2a=a3d-a2w1(3-4)&-a&-a4a=a1w1(3-5) J把上述二式写成状态方程，设状态变量为&，x=J(3-6) x1=a，x2=J3则状态方程&1=-a4x1+x2-a1w1(3-7) x&2=-a2x1-a1x2+a3d-a2w1(3-8) x&3=x2(3-9)

34、 x23/30把上述方程组写成矩阵形式：&=Ax+Bu-Tw(3-10) x式中x=x1，x2，x3，u=d，w=w1T-a1A=-a201-a1100a 0，B=300&，由于有测量误在飞机上用二自由度测速陀螺仪测量飞机的俯仰角速度J差，测速陀螺的输出为&+v=x+n z1=J121(3-11)&，由于有测量误差，三自由度在飞机上用三自由度陀螺测量飞机的俯仰角J陀螺的输出为z2=J+v2=x3+v2(3-12)把式(3-11)和式(3-12)联合起来得到矩阵形式z=Hx+v 式中(3-13)z1010v1z=，H=，v=v z00122设w1，v1，v2都是均值为零的白噪声，且互不相关，即E

35、w1=Ev1=Ev2=02Ewtwt=qd(t-t) ()()111r12Ev(t)v(t)=0T0d(t-t) 2r2TEwtv()(t)1=0&和q的估值。 求a、J解根据前面所述的连续系统卡尔曼滤波基本方程（）三式，可得最优滤波方程为24/30&(t/t)x1-a1&(t/t)=-ax22&0xt/t)3(1-a111(t/t)00xk112(t/t)+a3u+k210x(t/t)00xk3131(t/t)xk12z1010k22-xt/t)z0012(2k32xt/t()3(3-14)最优增益矩阵为k11K=k21k31k12p11pk22=21k32p31p12p22p32p1300

36、-210r1p230p3301p12r1-20-2=pr221r2-2p32r1-2p13r2-2p23r2-2p33r2-2(3-15)&=0时，可得代数黎卡提方程： 当滤波达到稳态时，即当P-a110p11p12p13p11-a-a0p+ppp2121222321100p31p32p33p31-a1p11p12q2-a-a0-p+-a21p2222110p31p32p11010001p21p31p12p22p32p13p23=0p33p13-a1-a201-a1ap22p231p32p33000p1300-r1-20p2310-20-r201p33p12(3-16)P为对称矩阵，因此p12

37、=p21，p13=p31，p32=p23设a1=1.76秒，a2=4.2秒，a3=7.4秒，a4=0.77秒q1=1.2度，r1=0.85度/秒，r2=0.85度求解式(3-16)得p11=0.86，p22=2.78/秒，p33=0.48 p31=p13=0.41秒，p21=p12=0.58，p23=p32=0.17求解式(3-15)得k1=0.81秒，k12=0.57，k21=3.8625/30k22=0.24秒，k21=0.23秒，k32=0.67&和q的估值。 把求得的K代入式(3-14)，可得a、J&=x&=x1(t/t)，J2(t/t)，J=x3(t/t)a3.2基于卡尔曼滤波方法的

38、时用水量预测（定量描述略）杭州位于钱塘江畔,是一个人口密集城市,其设计供水量能力为140万t/d。本文以杭州市的时用水量为例,分别采用卡尔曼滤波模型以及卡尔曼滤波季节模型进行预测,对结果进行比较分析(见二图),研究上述2种方法的利弊。从图中可以看出,对于1h预测,非季节模型误差有较大波动,改进后的模型,即26/30卡尔曼滤波季节模型相对较为平稳。总体来说季节模型预测精度相对较高。评价模型预测精度,还需查看其预测误差分布特征7。图3为卡尔曼滤波常规模型与季节模型预测误差分布比较,采用了1个月数据,用2种方法进行了预测,从图中可以看出,2种方法误差都主要集中在-0.1至0.1区域,呈对称分布。从1

39、个月所有数据来看,误差大致呈正态分布,季节模型较常规模型趋势图更为尖锐,表明季节模型比常规模型预测效果更为理想。综合上图和下图,比较得出,卡尔曼滤波季节模型具有更高的预测精度。采用基于卡尔曼滤波方法预测模型,进行时用水量预测。鉴于城市时用水量具有周期性、趋势性及随机扰动性,引入季节因子,形成卡尔曼滤波方法季节模型。以杭州市实测数据为例,分别采用卡尔曼滤波模型和卡尔曼滤波季节模型分别对该市时用水量进行预测,预测结果表明,卡尔曼滤波方法季节预测模型精度相对较高。从计算模型看,季节预测模型要求低,但预测精度较高,具有一定实用性,能够满足供水系统优化调度的要求。27/304卡尔曼滤波局限性分析4.1稳定性定理前面比较详细地讨论了线性系统卡尔曼滤波基本方程的推导和递推计算步骤。在估值计算时，需要利用一连串的观测数据，按照滤波基本方程进步递推计(t1/t)。这里，需要确切知道x(t)的初值和算，可得状态向量x(k)的最优估值x估计误差方差阵的初值P(0/0)。但在许多实际问题中，往往不可能确切知道初值x(0)和P(0/0)，甚至根本不知道这些初值。为了进步滤波计算，只能假定初