信号处理Chapter 3 kalman filter

第三章卡尔曼滤波（TheKalmanfiltering）第一节卡尔曼滤波信号模型第二节卡尔曼滤波方法第三节卡尔曼滤波的应用1、信号模型状态方程和量测方程•维纳滤波的模型：信号可以认为是由白噪声激励一个线性系统的响应，假设响应和激励的时域关系可以用下式表示：(6-52)上式也就是一阶AR模型。)(ns)(1nw)(zA)1()1()(1nwnasns•在卡尔曼滤波中信号被称为是状态变量，用矢量的形式表示为，激励信号也用矢量表示为，激励和响应之间的关系用传递矩阵来表示，得出状态方程：(6-53)上式表示的含义就是在k时刻的状态可以由它的前一个时刻的状态来求得，即认为k－1时刻以前的各状态都已记忆在状态中了)(nsS(k))(1nw(k)w1A(k)1)(kw1)A(k)S(kS(k)1S(k)1)S(k1)S(k卡尔曼滤波是根据系统的量测数据（即观测数据）对系统的运动进行估计的，所以除了状态方程之外，还需要量测方程。在卡尔曼滤波中，用表示量测到的信号矢量序列，表示量测时引入的误差矢量，则量测矢量与状态矢量之间的关系可以写成(6-54)w(k)S(k)X(k)上式和维纳滤波的概念上是一致的，也就是说卡尔曼滤波的一维信号模型和维纳滤波的信号模型是一致的。把式(6-55)推广就得到更普遍的多维量测方程(6-55)上式中的称为量测矩阵，它的引入原因是，量测矢量的维数不一定与状态矢量的维数相同，因为我们不一定能观测到所有需要的状态参数。w(k)S(k)X(k)w(k)C(k)S(k)X(k)信号模型根据状态方程和量测方程，卡尔曼滤波的信号模型，如图6.12所示。•图6.12卡尔曼滤波的信号模型1)(kw1)A(k)S(kS(k)1w(k)C(k)S(k)X(k)S(k)C(k)1)A(k1zw(k)(k)w1X(k)1)S(k2、卡尔曼滤波方法（ThemethodofKalmanfiltering）卡尔曼滤波的一步递推法模型把状态方程和量测方程重新给出：(6-56)(6-57)假设信号的上一个估计值已知，现在的问题就是如何来求当前时刻的估计值。1)(kw1)A(k)S(kS(k)1w(k)C(k)S(k)X(k)1)(kSˆ(k)Sˆ•用上两式得到的和分别用和表示，得：(6-58)(6-59)•必然，观测值和估计值之间有误差，它们之间的差称为新息（innovation）：(6-60)显然，新息的产生是由于我们前面忽略了与所引起的1)(kSA(k)(k)Sˆˆ1)(kSC(k)A(k)(k)SC(k)(k)XˆˆˆX(k)(k)Xˆ(k)X~(k)XX(k)(k)Xˆ~(k)w1w(k)•用新息乘以一个修正矩阵，用它来代替式（6－56）的来对进行估计：(6-61)由（6－56）～（6－61）可以画出卡尔曼滤波对进行估计的递推模型，如图6.13所示(k)X~H(k)(k)w1S(k)(k)XH(k)1)(kSA(k)(k)S~ˆˆ1)](kSC(k)A(k)H(K)[X(k)1)(kSA(k)ˆˆS(k)•输入为观测值，输出为信号估计值。图1卡尔曼滤波的一步递推法模型X(k)(k)Sˆ(k)SˆC(k)A(k)1zX(k)1)(kSˆH(k)(k)Xˆ(k)X~卡尔曼滤波的递推公式从图1容易看出，要估计出就必须要先找到最小均方误差下的修正矩阵，结合式（6－61）、（6－56）、（6－57）得：(6-62)根据上式来求最小均方误差下的，然后把求到的代入（6－61）则可以得到估计值。(k)SˆH(k)1)](kSC(k)A(k)w(k)(k)H(K)[C(k)S1)(kSA(k)(k)Sˆˆˆ1)](kSC(k)A(k)w(k)1)](kw1)(kSA(k)H(K)[C(k)[1)(kSA(k)1ˆˆˆH(k)w(k)1)](kw1)(kS(k)H(K)C(k)[AH(k)C(k)]1)[I(kSA(k)1ˆˆH(k)H(k)(k)Sˆ•设真值和估计值之间的误差为：误差是个矢量，因而均方误差是一个矩阵，用表示。把式（6－62）代入得(6-63)均方误差矩阵：(6-64)表示对向量取共轭转置。(k)SS(k)(k)Sˆ~ε(k)(k)SS(k)(k)Sˆ~H(k)w(k)1)](kw1)](kS1)A(k)[S(kH(K)C(k)][[I1ˆ](k)S(k)SE[ε(k)τ~~为了计算方便，令(6-65)找到和均方误差矩阵的关系：(6-66)把式（6－63）代入式（6－64），最后化简得：](k))S(k))(S(k)SE[(S(k)(k)ετˆˆ]1))(kSA(k)1)(kw1)k1))(A(k)S((kSA(k)1)(kw1)E[(A(k)S(k(k)ετ11ˆˆ]1)(k1)w(kE[w]A(k)1))(kS1)1))(S(k(kS1)A(k)E[S(kτ11ττˆˆ1)Q(k1)A(k)A(k)ε(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•把式（6－66）代入（6－67）得令，代入上式化简：(6-68)要使得均方误差最小，则必须ε(k)ττk)H(k)R(k)H(H(k)C(k)](k)[IεH(K)C(k)][IττττR(k)]H(k)(k)C(k)εH(k)[C(k)H(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)ετ1τ1ττ1τ])U(S][H(k)S)U(S[H(k)SU)U(SS(k)ε01τ)U(SH(k)S求得最小均方误差下的修正矩阵为：(6-69)把上式代入(6-61)即可得均方误差最小条件下的递推公式。最小均方误差为：(6-70)1ττR(k)](k)C(k)ε[C(k)(k)C(k)εH(k)(k)Sˆε(k)τ1τU)U(SS(k)ε(k)εH(k)C(k)][I•综上所述，得到卡尔曼滤波的一步递推公式：(6-71)(6-72)(6-73)(6-74)(k)ε1)Q(k1)A(k)A(k)ε(kτ1ττR(k)](k)C(k)ε[C(k)(k)C(k)εH(k)ε(k)(k)εH(k)C(k)][I(k)Sˆ1)](kSC(k)A(k)H(K)[X(k)1)(kSA(k)ˆˆ【例】设卡尔曼滤波中量测方程为已知信号的自相关函数的z变换为，噪声的自相关函数为，信号和噪声统计独立，已知在k＝0时刻开始观测信号。试用卡尔曼滤波的公式求和，k＝0,1,2,3,4,5,6,7；以及稳态时的和。解：由例6-6的结果知，w(k)S(k)X(k)25.18.0,)8.01)(8.01(36.0)(1zzzzRss)()(mmRww1)0(,0)1(ˆS(k)Sˆε(k)(k)Sˆε(k)8.0A(k)1C(k)36.021wQ(k)1))(var(kwR(k)•把上式代入式(6-71)～(6-74)得（1）（2）（3）（4）求逆把（1）代入（2）、（3）式，消去，再把（2）和（3）联立，得到（5）(k)ε36.064.01)ε(k1](k)ε(k)[εH(k)1ε(k)(k)εH(k)][1(k)Sˆ1)](kSH(K)[X(k)1)(kSˆ8.0ˆ8.01](k)ε(k)[εH(k)1](k)ε[(k)ε1/(k)εH(k)1.361)ε(k0.361)ε(kε(k)64.064.0初始条件为，k＝0开始观测，利用等式（4），（5）进行递推得：•k＝0，1.0000，1.0000，•k＝1，0.5000，0.5000，•k＝2，0.4048，0.4048，•k＝3，0.3824，0.3824，•k＝4，0.3768，0.3768，•k＝5，0.3755，0.3755，•k＝6，0.3751，0.3751，•k＝7，0.3750，0.3750，•上面是递推过程，还没有达到稳态的情况。1)0(,0)1(ˆS)0()0(H)0(0ˆX)(S)1(）（1H)1(5.0)0(ˆ4.01ˆXS)(S)2()2(H)2(4048.0)1(ˆ4762.02ˆXS)(S)3()3(H)3(3824.0)2(ˆ4941.03ˆXS)(S)4()4(H)4(3768.0)3(ˆ4985.04ˆXS)(S)5()5(H)5(3755.0)4(ˆ4996.05ˆXS)(S)6()6(H)6(3751.0)5(ˆ4999.06ˆXS)(S)7()7(H)7(3750.0)6(ˆ5000.07ˆXS)(S•假设到了某一时刻k－1，前后时刻的均方误差相等，也就是误差不再随着递推增加而下降，达到最小的均方误差了，即稳态情况，式（5）中的误差代入（5）式可以计算到稳态时的均方误差为：即稳态时的修正矩阵，代入式4得稳态时的信号估计：化到z域有：。)1()(kk375.0)1()(kk375.0)(kHX(k)1)(kS375.0ˆ5.0(k)Sˆ)(zH15.01375.0z３卡尔曼滤波器的应用(ApplicationKalmanfilter)【例】已知条件和例6－2一样，状态方程和测量方程为：其中，信号和噪声统计独立。求卡尔曼滤波器的稳态和。1)(kw1)A(k)S(kS(k)1w(k)C(k)S(k)X(k)8.0A1C36.021wQ(k)1))(var(kwR(k)H(k)ε(k)•解：根据函数调用sys＝ss（A，B，C，D，1），得到离散卡尔曼状态模型，采样周期这里设为1。A，C已知，由于函数调用中是设计了两个观测信号的，我们这里只有一个观测信号，所以B取[01]，后一个1表示噪声的系数。D取0。实际的语句如下：•sys=ss(A,B,C,D,1)•然后调用函数[S，L，，H，]＝kalman（sys，Q，R），设计离散卡尔曼滤波器。实际语句和计算结果如下：•[s，l，，h，]=kalman(sys,0.36,1)•l=0.3000•=0.6000•h=0.3750•=0.3750•这里省略了输出的S，它表示的信息是达到稳态后系统状态模型，H和表示系统稳态的最终值)1(1kw•有了修正矩阵和均方误差，代入式（6－74）就可以根据观测信号得到卡尔曼滤波的估计值了。•从上面例题知道，只要确定了状态模型，就可以调用函数很快设计出卡尔曼滤波器，下面来看看卡尔曼滤波器在生物医学信号中的应用。•在生物医学信号处理中脑电图的肌电伪迹和其它噪声的消除，以及诱发电位的提取都有研究者尝试用卡尔曼滤波器来处理。本节介绍卡尔曼滤波器在诱发电位提取中的应用，方法如下：1.自发电位模型（EEG）和诱发电位（EP）模型的建立。•如图6.14所示，EEG信号通过用AR模型建立，激励是白噪声，EP信号的激励是单位脉冲序列，用等式表示如下：阶AR模型d表示从该时刻开始有单位脉冲刺激。pnwineaneEEGpii),()()(:1qiimiiidndinscnsEP01)()()(:图2EEG和EP模型•从图2知道，观测信号是EEG和EP的线性相加，用表示第i次刺激后测量的信号，对M次测量平均得：•叠加平均后的信号长度为N。利用先验知识建立好图6.14