维纳滤波器和卡尔曼滤波器

1、 ( )( )( )x ns nw n (7-1)则输出( )y n为( )( )( )( ) ()my nx nh nh m x nm (7-2)()()(nwnsnx)( )(nsny )(nh图7-1 维纳滤波器的输入输出关系)( )()(nsnsne (7-3)22)( )()(nsnsEneE(7-4) 0n当0,h(n)0, 0)(nnh当0mm)h(m)x(n(n)s y(n)(7-5)202) )()()()(mmnxmhnsEneE(7-6)2 , 1 , 00)() )()()(20jjnxmnxmhnsEmopt(7-7)即0)()()()()(0jjnxmnxEmhjn

2、xnsEmopt(7-8)用相关函数R来表达上式，则得到维纳霍夫方程的离散形式：0)()()(0jmjRmhjRmxxoptxs从维纳霍夫方程中解出的h就是最小均方误差下的最佳h， 即)(nhopt 求到 )(nhopt，这时的均方误差为最小： 20min2) )()()()(moptmnxmhnsEneE )()()()()()()(2)(0002mroptoptmrnxrhmnxmhmnxmhnsnsE000)()()()()(2)0(mrxxoptoptmxsoptssrmRrhmhmRmhR由式（7-9）进一步化简得：0min2)()()0()(mxsoptssmRmhRneE(7-1

3、0)10)()()( )(Nmmnxmhnsny (7-11)2102) )()()()(NmmnxmhnsEneE (7-12) 12 , 1 , 00)() )()()(210NjjnxmnxmhnsENmopt(7-13)1, 1 , 0)()()()()(10NjjnxmnxEmhjnxnsENmopt .(7-14)1, 2 , 1 , 0)()()(10NjmjRmhjRNmxxoptxs (7-15)于是得到N个线性方程：)0() 1()2() 1 () 1()0() 1(1)2() 1()0() 1 () 1 ()0() 1 (1) 1() 1() 1 () 1 ()0()0(

4、)0(0 xxxxxxxsxxxxxxxsxxxxxxxsRNhNRhNRhNRNjNRNhRhRhRjNRNhRhRhRj写成矩阵形式有： ) 1() 1 ()0() 1() 1 ()0()0()2() 1()2()0() 1 () 1() 1 ()0(NRRRNhhhRNRNRNRRRNRRRxsxsxsxxxxxxxxxxxxxxxxxx(7-16)简化形式：RxxH=Rxs (7-17)式中，Hh(0) h(1) h(N-1)，是待求的单位脉冲响应；Rxs) 1(),1 (),0(NRRRxsxsxs，是互相关序列； Rxx )0()2() 1()2()0() 1 () 1() 1 (

5、)0(xxxxxxxxxxxxxxxxxxRNRNRNRRRNRRR，是自相关矩阵 只要Rxx是非奇异的，就可以求到H：H=Rxx1 Rxs (7-18)求得)(nhopt后，这时的均方误差为最小： 210min2) )()()()(NmoptmnxmhnsEneE )()()()()()()(2)(1010102NmNroptoptNmrnxrhmnxmhmnxmhnsnsE101010)()()()()(2)0(NmNrxxoptoptNmxsoptssrmRrhmhmRmhR 10min2)()()0()(NmxsoptssmRmhRneE (7-19)用有限长的)(nh来实现维纳滤波时

6、，当已知观测值的自相关和观测值与信号的互相关时就可以按照式（7-15）在时域里求解 )(nhopt但是当N比较大时，计算量很大，并 且涉及到求自相关矩阵的逆矩阵问题。 )(ns)(nw0)()(mRmRwssw则有)()()()()()()()(mRmnsnwmnsnsEmnsnxEmRssxs)()()()()()()(mRmRmnwmnsnwnsEmRwwssxx则式（7-15）和式（7-19）化为：1, 2 , 1 , 0)()()()(10NjmjRmjRmhjRwwNmssoptss(7-20) 10min2)()()0()(NmssoptssmRmhRneE (7-21)【例7-1

7、】已知图7-1中 )()()(nwnsnx且 )(ns与 )(nw统计独立，其中 )(ns的自相关序列为 mssmR6 . 0)()(nw， 是方差为1的单位白噪声，维纳滤波器来估计 )(ns，并求最小均方误差。 解依题意，已知信号的自相关和噪声的自相关为： mssmR6 . 0)()()(mmRww，代入式（7-20）得 ) 1 (2)0(6 . 06 . 01) 1 (6 . 0)0(210hhjhhj解得： )0(h0.451， ) 1 (h0.165。 将上述结果代入式（7-21），求得最小均方误差：45. 0) 1 (6 . 0)0(1)()()0()(10min2hhmRmhRne

8、Emssss若要进一步减小误差可以适当增加维纳滤波的阶数，但相应的计算量也会增加。 提出的白化的方法求解维纳霍夫方程，得到系统函数 H(z)。 由第三章的知识，我们知道随机信号都可以看成是由一白色噪声 w1(n) 激励一个物 理可实现 的系统 或模型的响应，如图7-2所示， 其中A(z) 表示系统 的传递函数。由于 x(n) = s (n) + w(n)， 在图7-2的基础上给出x(n)的信号模型 )(ns)(zA)(1nw 图7-3 )(nx的信号模型 )(ns)(ns)(zA)(1nw)(nw)(nx图7-4 维纳滤波器的输入信号模型)(nx)(zB)(1nw)()(2111mmRwww2

9、1w)(mRss )()()()()()()(11rkssrmnwraknwkaEmnsnsEmR)()()(11rkmRrakawwkr令 krl上式 lkwwwwkllkakalmRlmRlkaka)()()()()()(1111令 )()()()()(lalalkakalfl)()()()()()()()(111111mamamRmfmRlflmRmRwwwwlwwss对式（7-22）进行Z变换得到系统函数和相关函数的z变换之间的关系：)()()(121zAzAzRwss (7-23)同样，对图7-4进行z变换得)()()(121zBzBzRwxx(7-24)图7-4中利用卷积性质还可以

10、找到互相关函数之间的关系：)()()()()()(1mnsknwkbEmnsnxEmRkxs)()()()(11mbmRkmRkbswksw两边z变换得到)()()(11zBzRzRswxs (7-25)如果已知观测信号的自相关函数，求它的z变换， 然后找到该函数的成对零点、极点，取其中在单位圆内的那一半零点 、极点构成)(zB，另外在)(zB)(1zB单位圆外的零、极点构成，这样就保证了 是因果的，并且是最小相位系统。从图7-4可得)()(1)(1zXzBzW (7-26)由于系统函数)(zB的零点和极点都在单位圆内，)(1zB)(nx即是一个物理可实现的最小相位系统，则也是一个物理可实现的

11、最小相移网络函数。我们就可以利用式（7-26）对 )(nx)(1nw)(1zB进行白化，即把当作输入，当作输出，是系统传递函数。)(zH)(zHopt)(1zB)()()(zBzGzH (7-27)()()(nwnsnx)( )(nsny)(nh（a）)(nx)( )(nsny)(1zB)(zG)(1nw（b）图7-5 利用白化方法求解模型)(nx)(mRxx)(zRxx)()()(121zBzBzRwxx)(zB)()()(zBzGzH)(zB)(zHopt01)()()( )(mmnwmgnsny (7-28)均方误差为：2012) )()()()(mmnwmgnsEneE )()()()

12、()()()(2)(0011012mrmrnwrgmnwmgmnwmgnsnsE000)()()()()(2)0(111mrwwmswssrmRrgmgmRmgR)()(2111mmRwww02202)()()(2)0()(11mwmswssmgmRmgRneE02202)(1)()()0(11111mswwmwswwssmRmRmgR均方误差最小也就是上式的中间一项最小，所以0,)()(211mmRmgwswopt (7-30)注意，这里的)(mg是因果的。对该式求z变换，得到211)()(wswoptzRzG (7-31)(1zRsw)(1mRsw表示对求单边z变换。所以维纳霍夫方程的系统

13、函数解表示为：)()()(zBzGzHoptopt)()(211zBzRwsw由式（7-25）上式可以表示为：)()()(zBzGzHoptopt)()(211zBzRwsw)()(/ )(211zBzBzRwxs因果的维纳滤波器的最小均方误差为： min2)(neE022)(1)0(11mswwssmRRmswwssmumRR)()(1)0(2211 (7-33)利用帕塞伐尔(Parseval)定理，上式可用z域来表示：min2)(neEcxsoptsszdzzRzHzRj)()()(211(7-34) min2)(neEcxsoptsszdzzRzHzRj)()()(211 (7-34)(

14、)()(nwnsnx)(ns)(nw)(nsmssmR8 . 0)()(nw)(ns【例7-2】已知图7-1中，且统计独立，其中的自相关序列为，是方差为1的单位白噪声，试，并求最小均方误差。与设计一个物理可实现的维纳滤波器来估计mssmR8 . 0)()()(mmRww0)(mRsw)()(mRmRssxs解依题意，已知，步骤1)()()(mRmRmRwwssxx求z变换25. 18 . 0 ,)8 . 01)(8 . 01 ()5 . 01)(5 . 01 (6 . 11)8 . 01)(8 . 01 (36. 0)(111zzzzzzzzRxx步骤2由于)()()(121zBzBzRwxx

15、，容易找到最小相位系统和白噪声方差zzzzB8 . 0 ,8 . 015 . 01)(1125. 1,8 . 015 . 01)(1zzzzB6 . 121w步骤3利用式（732） )(zHopt)()(/ )(211zBzBzRwxs)5 . 01)(8 . 01 (36. 0)5 . 01 (6 . 18 . 01111zzzz对括号里面求反变换，注意括号内的收敛域为28 . 0 z)5 . 01)(8 . 01 (36. 011zzZ) 1()2(6 . 0)()8 . 0(6 . 0nununn取因果部分，也就是第一项，所以118 . 0116 . 0)5 . 01)(8 . 01 (

16、36. 0zzz)(zHopt0,)5 . 0(375. 0)(nnhn步骤4最小均方误差为 min2)(neEcxsoptsszdzzRzHzRj)()()(211cdzzzzzj)5 . 0)(25. 1)(8 . 0()5 . 0625. 0(45. 021取单位圆为积分围线，有两个单位圆内的极点，0.8和0.5，求它们的留数和，所以min2)(neE375. 0)25. 15 . 0)(8 . 05 . 0()5 . 05 . 0*625. 0(45. 0)5 . 08 . 0)(25. 18 . 0()5 . 08 . 0*625. 0(45. 0)(nx) 1( nx)2( nx)

17、( )(nsny)( )(Nnsny0)(nyd)(ny)()()(nwnsnx)( )(Nnsny)(nh)()(Nnsnyd图7-6维纳预测器本节和上节一样着重讨论预测器的系统函数以及预测的均方误差，维纳预测器和维纳滤波器比较类似，因而分析方法也都可以借鉴前面的内容。 )(nh0, 0)(nnh当对于图7-6模型，设是物理可实现的，也即，则有是因果序列：0)()()( )(mmnxmhNnsny (7-35)202) )()()()(mmnxmhNnsEneE (7-36)要使得均方误差最小，则将上式对各)(mh，m0，1，求偏导，并且等于零，得：)(mh2 , 1 , 00)() )()

18、()(20jjnxmnxmhNnsEmopt (7-37) 用相关函数R来表达上式：0)()()()()(0jjnxmnxEmhjnxNnsEmopt(7-38)0)()()(0jmjRmhjNRmxxoptxs (7-39)()(Nnsnyd)()()()(NmRNmnsnxEmRxsxyd由于，则，z变换得)()(zRzzRxsNxyd)()(11zRzzRxsNxyd (7-40)()(Nnsnyd)()(nsnyd借鉴维纳滤波器的结果类似给出维纳预测器的最佳传递函数，对应维纳预测器，对应维纳滤波器，故因果的预测器的传递函数为： )(zHopt)()(/ )(211zBzBzRwxyd)

19、()(/ )(211zBzBzRzwxsN (7-41)最小均方误差为 min2)(NneE022)(1)0(11mywwssmRRd (7-42)利用帕塞伐尔(Parseval)定理，上式可用z域来表示min2)(NneEcxyoptsszdzzRzHzRjd)()()(211cxsNoptsszdzzRzzHzRj)()()(211)()()(nwnsnx)(ns)(nw)(nsmssmR8 . 0)()(nw) 1( ns【例7-3】已知图7-6中，且与统计独立，其中的自相关序列为，是方差为1的单位白噪声，并求最小均方误差。试设计一个物理可实现的维纳预测器估计mssmR8 . 0)()(

20、)(mmRww0)(mRsw解依题意已知， )()(mRmRssxs)()()(mRmRmRwwssxx求z变换：25. 18 . 0 ,)8 . 01)(8 . 01 ()5 . 01)(5 . 01 (6 . 11)8 . 01)(8 . 01 (36. 0)(111zzzzzzzzRxx由于)()()(121zBzBzRwxx，容易找到最小相位系统和白噪声方差：zzzzB8 . 0 ,8 . 015 . 01)(1125. 1,8 . 015 . 01)(1zzzzB6 . 121w由式（7-41），N1，)(zHopt)()(/ )(211zBzBzzRwxs)()(/ )(211zB

21、zBzzRwss求z变换：25. 18 . 0 ,)8 . 01)(8 . 01 ()5 . 01)(5 . 01 (6 . 11)8 . 01)(8 . 01 (36. 0)(111zzzzzzzzRxx由于)()()(121zBzBzRwxx，容易找到最小相位系统和白噪声方差：zzzzB8 . 0 ,8 . 015 . 01)(1125. 1,8 . 015 . 01)(1zzzzB6 . 121w由式（7-41），N1，)(zHopt)()(/ )(211zBzBzzRwxs)()(/ )(211zBzBzzRwss)5 .01)(8 .01 (36.0)5 .01 (6 .18 .01

22、111zzzzz对括号里面求z反变换，注意括号内的收敛域为：28 . 0 z，)5 . 01)(8 . 01 (36. 011zzzZ) 1()2(2 . 1)()8 . 0(48. 0nununn取因果部分，也就是第一项，所以118.01148.0)5.01)(8.01(36.0zzzz)(zHopt11115 . 013 . 0)8 . 01 (48. 0)5 . 01 (6 . 18 . 01zzzz0,)5 . 0(3 . 0)(nnhn把上式写成差分方程形式有：)( 5 . 0)(3 . 0) 1( nsnxnsmin2)(neEcssoptsszdzzRzzHzRj)()()(21

23、11最小均方误差为：cdzzzj6 . 0)8 . 01)(5 . 0(36. 021)(nw)(Nns)()(nsnx)( )(Nnxny)(nh)()(Nnxnyd图7-7 N步纯预测器)()(nsnx)()()(121zBzBzRwxx)()()()()(121zBzBzRzRzRwxsssxx这时，用白化法来求解预测器的系统函数。因为，从而有： (7-44)(zHopt)()(/ )(211zBzBzRzwxsN)()()()(/ )()(211211zBzBzzBzBzBzBzNwwN将上式代入式（7-41）、（7-43）得：min2)(NneEcxsNoptsszdzzRzzHzR

24、j)()()(211cNNwzdzzBzzBzzBzBj)()()()(21121假设B(z）是b（n）的z变换，且b（n）是实序列，则上式可以利用帕塞伐尔定(Parseval)理进一步化简：min2)(NneE )()()()(221nnwNnbnuNnbnb (7-46)又因为B(z）是最小相位系统，一定是因果的，上式可以简化min2)(NneE102202022)( )()(11NnwnnwnbNnbnb (7-47)上式说明最小均方误差随着N的增加而增加，也即预测距离越远误差越大。 )()(nsnx)(nsmssmR8 . 0)()( NnsmssmR8 . 0)(【例7-4】已知图7

25、-7中，其中的自相关序列为，试设计一个物理可，并求最小均方误差。，则实现的维纳预测器来估计解依题意，已知25. 18 . 0 ,)8 . 01)(8 . 01 (36. 0)(1zzzzRss因为 )()()()()(121zBzBzRzRzRwxsssxx容易找到最小相位系统和白噪声方差:zzzB8 . 0 ,8 . 011)(1)(8 . 0)(nunbn25. 1,8 . 011)(1zzzB36. 021w利用式（7-45）:)(zHopt8 . 011)8 . 01 ()()(11zzzzBzBzNN)(8 . 08 . 01111NnuzzZNnN0n因为 ，只取的部分，有:)(8

26、 . 08 . 01111nuzzZNnN)8 . 0118 . 08 . 01111zzzNN)(zHoptNoptzH8 . 0)(回到z域有：，代入得:最小均方误差为：min2)(NneENNnnNnwnunb210210228 . 01)(8 . 036. 0)(1它说明当N越大，误差越大，当N0时，没有误差。把上述结果用模型表示如图7-8所示。)()(nsnx)(8 . 0)( )(nsNnsnyN)(8 . 0)(nnhNopt图7-8 例题7-3的纯预测模型0)()()( )(mmnxmhNnxnypmmnxmhnxny1)()()( )(这就是一步线性预测公式，常常用下列符合表

27、示pmpmmnxanx1)()( (7-48)式中p为阶数，)(mhapm。预测的均方误差为:) )()()(212pmpmmnxanxEneEpmplxxplpmpmxxpmxxmlRaamRaR111)( )( 2)0( (7-49)要使得均方误差最小，将上式右边对pma求偏导并且等于零，得到p个等式 pmaplmlRalRpmxxpmxx, 2 , 10)()(1 (7-50)最小均方误差：pmxxpmxxmRaRneE1min2)()0()( (7-51)式（7-50）就是YuleWalker（Y-W）方程，和第三章AR模型参数估计的方程一致，如何去求解该方程在第三章有详细介绍。把Yu

28、leWalker（Y-W）方程和维纳霍夫方程进行比较，维纳霍夫方程要估计的量是s(n)，Y-W方程要估计的量是x(n)本身，因而解维纳霍夫方程要已知x(n)、y(n)的互相关函数，实际中这个互相关函数往往是未知的，而解Y-W方程只需要知道观测信号的自相关函数。因此Y-W方程比W-H方程更具有实用价值。)( nsmssmR8 . 0)(mssxxmRmR8 . 0)()(64. 0)2(, 8 . 0) 1(, 1)0(xxxxxxRRR例7-5】已知图7-7中x(n)=s(n)，其中的自相关序列为的可实现的一步线性预测器，并求最小均方误差。解， ，试设计一个p2利用Y-W方程2 , 10)()

29、(21lmlRalRmxxpmxx，可以列出2个方程式08 . 064. 008 . 08 . 02121ppppaaaa0, 8 . 021ppaa) 1(8 . 0)( nxnx解得：，也即36. 064. 01)()0()(1min2pmxxpmxxmRaRneE结果和例（7-4）N=1时一致。),(),(),(ftWftSftXii),(),(),(1),(1ftWftSftXNftXNiiiNitwNtstxiNi, 2 , 1),(1)()(1对每次观测用短时傅立叶变换求时频表示（TFR）：对N次观测的时频表示（TFR）求平均：，样本平均的时频表示（TFR）为： (1)样本平均为：

30、),(1),(),(ftWNftSftX (2)从式（2）可以得到一个基于样本平均的简单时频平面后验维纳滤波器：),(1),(),(),(ftWNftSftSfth (3) ),(1),(),(1),(),(),(11NiiiNiiiftXIFNftWftSCOVNftWftSftX .(4),(),(),(),(),(),(ftXIFftWftSCOVftWftSftX (5),(ftW),(ftWi式中COV表示信号和噪声之间的方差，也就是考虑了信号和噪声并非相互独立；IF是干扰项；表示样本平均的噪声功率；表示样本噪声功率的平均。)(nx),2(),1(nxnx)(ns)(zH)(nh)

31、1( ns)(nx)( ns它是根据前一个估计值和最近一个观测值来估计信号的当前推方法进行估计的，因而卡尔曼滤波对信号的平稳性和时不变性不做要求。我们利用维纳滤波的模型引入到卡尔曼滤波的信号模型。，它是用状态方程和递) 1() 1()(1nwnasns)(ns)(1nw)(zA上式也就是一阶AR模型。在卡尔曼滤波中信号)(ns被称为是状态变量，用矢量的 S(k)S(k)形式表示为，在k时刻的状态用表示， 在k1时刻的状态用1)S(k 表示。 )(1nw(k)w1A(k)激励信号也用矢量表示为，激励和响应之间的关系用传递矩阵来表 示，它是由系统的)(zA有一定关系。有了这些假设后结构确定的，与我

32、们给出状态方程：1)(kw1)A(k)S(kS(k)1 (7-53)S(k)1)S(k 1)S(k 上式表示的含义就是在k时刻的状态可以由它的前一个时刻的状态来求得，即认为k1时刻以前的各状态都已记忆在状态中了 )(ns)(zA)(1nw)(ns)(zA)(1nw)(nw)(nx图7-11 维纳滤波的信号模型和观测信号模型)()()(nwnsnx)(nw卡尔曼滤波是根据系统的量测数据（即观测数据）对系统的运动进行估计的，所以除了状态方程之外，还需要量测方程。还是从维纳滤波的观测信号模型入手，图7-11的右图，观测数据和信号的关系为：，一般是均值为零的高斯X(k)S(k)白误差矢量，则量测矢量与

33、状态矢量w(k)S(k)X(k) (7-54)上式和维纳滤波的)()()(nwnsnx概念上是一致的，也就是说卡尔曼滤波的一维信号模型和维纳滤波的信号模型是一致的。把式(7-54)推广就得到更普遍的多维量测方程w(k)C(k)S(k)X(k) (7-55)C(k)X(k)S(k)上式中的称为量测矩阵，它的引入原因是，的维数不一定与状态矢量的维数相同，因为我们不一定能观测到所有需要的状态参数 量测矢量X(k)1mS(k)1nC(k)nmw(k)1m假如是的矢量，是的矢量，就是的矩阵，是的矢量。1)(kw1)A(k)S(kS(k)1w(k)C(k)S(k)X(k)S(k)C(k)1)A(k 1zw

34、(k)(k)w1X(k)1)S(k 图7-12 卡尔曼滤波的信号模型w(k)S(k)X(k)25. 18 . 0 ,)8 . 01)(8 . 01 (36. 0)(1zzzzRss)()(mmRwwA(k)C(k)【例7-6】设卡尔曼滤波中量测方程为，已知信号的自相关函数的z变换为，噪声的自相关函数为：，信号和噪声统计独立。求卡尔曼滤波和。信号模型中的)()()(121zAzAzRwss)8 . 01)(8 . 01 (36. 011zzzz解根据等式：可以求得：)()(8 . 01)(111zWzSzzzA)()(8 . 0) 1(1kwksks8 . 0A(k)w(k)S(k)X(k)C(

35、k)变换到时域得：因此 又因为，所以1。S(k)(k)S1)(kw1)A(k)S(kS(k)1 (7-56)w(k)C(k)S(k)X(k) (7-57)A(k)C(k)X(k)1)(kS(k)S上式中和是已知的，已知，现在的问题就是如何来求当前时刻。 是观测到的数据，也是已知的，假设信号的 上一个估计值 的估计值(k)w1w(k)S(k)(k)w1w(k)(k)S(k)X(k)S(k)X上两式中如果没有与，可以立即求得，估计问题的出现就是因为信号与噪声的与，用上两式和分别用和表示，得：叠加。假设暂不考虑得到的1)(kSA(k)(k)S (7-58)1)(kSC(k)A(k)(k)SC(k)(

36、k)X (7-59)X(k)(k)X必然，观测值和估计值之间有误差 ，它们之间的差(k)X称为新息（innovation）：(k)XX(k)(k)X (7-60)(k)w1w(k)显然，新息的产生是由于我们前面忽略了与所引起的，也就是说新息里面 (k)w1w(k)(k)XH(k)(k)w1S(k)包含了与的信息成分。因而我们用新息乘以一个修正矩阵，用它来代替式（7-56）的来对进行估计：(k)XH(k)1)(kSA(k)(k)S1)(kSC(k)A(k)H(K)X(k)1)(kSA(k) (7-61)S(k)X(k)(k)S(k)SC(k)A(k)1zX(k)1)(kSH(k)(k)X(k)X

37、由（7-56）（7-61）可以画出卡尔曼滤波对进行估计的递推模型，如图7-13所示，输入为观测值，输出为信号估计值。图7-13 卡尔曼滤波的一步递推法模型(k)SH(k)1)(kSC(k)A(k)w(k)(k)H(K)C(k)S1)(kSA(k)(k)S1)(kSC(k)A(k)w(k)1)(kw1)(kSA(k)H(K)C(k)1)(kSA(k)1H(k)w(k)1)(kw1)(kS(k)H(K)C(k)AH(k)C(k)1)I(kSA(k)1.(7-62)H(k)H(k)根据上式来求最小均方误差下的，然后把求到的代入（7-61）则可以得 到估计值(k)S。(k)SS(k)(k)S(k)设真

38、值和估计值之间的误差为：，误差是个矢量，因而均方误差是一个矩阵，用表示。把式（7-62）代入得：(k)SS(k)(k)SH(k)w(k)1)(kw1)(kS1)A(k)S(kH(K)C(k)I1.(7-63)均方误差矩阵：(k)S(k)SE(k) (7-64)表示对向量取共轭转置。为了计算方便，令(k)S(k)(S(k)SE(S(k)(k) (7-65) 找到和均方误差矩阵的关系：1)(kSA(k)1)(kw1)k1)(A(k)S(kSA(k)1)(kw1)E(A(k)S(k(k)111)(k1)w(kEwA(k)1)(kS1)1)(S(k(kS1)A(k)ES(k111)Q(k1)A(k)A

39、(k)(k (7-66)(k)w1w(k)j)Q(k)(k(j)(k)wEw11把式（7-63）代入式（7-64），并且利用条件：与都是零均值的高斯白噪声，且它们和 互不相关，协方差矩阵分别为j)R(k)(kEw(k)w(j)1)S(k)(kw111)(kS)(kw11w(k)与不相关；与及不相关。最后化简得：(k)S(k)SE(k)k)H(k)R(k)H(H(k)C(k)1)IQ(k1)A(k)A(k)(kH(K)C(k)I. (7-67)把式（7-66）代入（7-67）得(k)k)H(k)R(k)H(H(k)C(k)(k)IH(K)C(k)IR(k)H(k)(k)C(k)H(k)C(k)H

40、(k)(k)C(k)(k)H(K)C(k)(k)SSR(k)(k)C(k)C(k)(k)C(k)U令 ，代入上式化简： (k)H(k)H(k)SSUH(k)H(K)U(k)111)U(SH(k)S)U(SH(k)SU)U(SS(k) (7-68)01)U(SH(k)SH(k)上式第一项和第二项与修正矩阵无关，第三项是，于是可以求得最小均方误差下的修正矩阵为：H(k)半正定矩阵，要使得均方误差最小，则必须1R(k)(k)C(k)C(k)(k)C(k)H(k) (7-69)把上式代入(7-61)即可得均方误差最小条件下的(k)S递推公式。相应的式（7-68）的第三项为零，得最小均方误差为：(k)1

41、U)U(SS(k)(k)H(k)C(k)I(7-70)综上所述，得到卡尔曼滤波的一步递推公式：(k)1)Q(k1)A(k)A(k)(k (7-71)1R(k)(k)C(k)C(k)(k)C(k)H(k) (7-72)(k)(k)H(k)C(k)I (7-73)(k)S1)(kSC(k)A(k)H(K)X(k)1)(kSA(k) (7-74)(k)S(k)S(0)0(0SE)(S有了上面四个递推公式后我们就可以得到和。如果初始状态的统计特性已知，并且令)0(var()0()0()(0()0()0(SSSSSE且矩阵 R(k)C(k),Q(k),A(k),都是已知的，以及观测量X(k)也是已知的，

42、就能用递推 (k)S(k)0() 1 (计算法得到所有的和：将初始条件代入式（7-71）求得； ) 1 () 1 (H) 1 () 1 (H) 1 ()0(0SE)(S) 1 (H)(S 1将代入式（7-72）求得和代入式（7-73）求得；将初始条件和代入式（7-74）求得；依此类推。这样递推用计算机实现 ；将非常方便。和维纳滤波一样，卡尔曼滤波也可以推广到卡尔曼预测，推导过程和维纳滤波到维纳预测类似，也同样有纯卡尔曼预测，这里不再推导。w(k)S(k)X(k)25. 18 . 0 ,)8 . 01)(8 . 01 (36. 0)(1zzzzRss)()(mmRww1)0(, 0) 1(S(k

43、)S(k)(k)S(k)【例7-7】设卡尔曼滤波中量测方程为，已知信号的自相关函数的z变换为，噪声的自相关函数为，信号和噪声统计独立，已知，在k0时刻开始观测信号。试用卡尔曼滤波的公式求和，k0,1,2,3,4,5,6,7；以及和。稳态时的8 . 0A(k)C(k)36. 021wQ(k)1)(var(kwR(k)解由例7-6的结果知，1，把它们代入式(7-71) (7-74)得(k)36. 064. 01)(k（1）1(k)(k)H(k)1（2）(k)(k)H(k) 1 （3） (k)S1)(kSH(K)X(k)1)(kS8 . 08 . 0（4） 1(k)(k)H(k)1(k) (k)1/

44、(k)由于是一维情况，求逆，把（1）代入（2）、（3）式，消去，再把（2）和（3）联立，得到H(k)1.361)(k0.361)(k(k)64. 064. 0 （5） 初始条件为1)0(, 0) 1(S，k0开始观测，利用等式（4），（5）进行递推得：)0()0(H)0(0X)(Sk0，1.0000，1.0000，；) 1 () 1 (5 . 0)0(4 . 01XS)(S k1，0.5000，0.5000，) 1 (H)2()2(H)2(4048. 0) 1 (4762. 02XS)(Sk2，0.4048，0.4048， )3()3(H)3(3824. 0)2(4941. 03XS)(Sk3

45、，0.3824，0.3824， )4()4(H)4(3768. 0)3(4985. 04XS)(Sk4，0.3768，0.3768， )5()5(H)5(3755. 0)4(4996. 05XS)(Sk5，0.3755，0.3755， )6()6(H)6(3751. 0)5(4999. 06XS)(Sk6，0.3751，0.3751，；)7()7(H)7(3750. 0)6(5000. 07XS)(Sk7，0.3750，0.3750， 如果给定每个时刻的观察值就可以得到每一时刻的信号估计值，上面是递推过程，还没有达到稳态的情况。假设到了某一时刻k1，前后时刻的均方误差相等，也就是误差不再随着递推增加而下降，达到最小的均方误差了，即稳态情况，式（5）中的误差) 1()(kk，代入（5）式可以计算到稳态时的均方误差为 375. 0) 1()(kk即稳态时的修正矩阵375. 0)(kH，代入式（4）得稳态时的信号估计： (k)SX(k)1)(kS375. 05 . 0)(zH15 . 01375. 0z化到z域有：。)(zHopt15 . 018/3z0,)5 . 0(375. 0)(nnhn