MDSP第四章维纳滤波和卡尔曼滤波

第二章维纳滤波和卡尔曼滤波2.1引言2.2维纳滤波器的离散形式——时域解2.3离散维纳滤波器的z域解2.4维纳预测2.5卡尔曼(Kalman)滤波2.1引言在生产实践中，我们所观测到的信号都是受到噪声干扰的。如何最大限度地抑制噪声，并将有用信号分离出来，是信号处理中经常遇到的问题。换句话说，信号处理的目的就是要得到不受干扰影响的真正信号。相应的处理系统称为滤波器。这里，我们只考虑加性噪声的影响，即观测数据()xn是信号()sn与噪声()vn之和（如图2.1.1所示），即()()()xnsnvn(0.0.1)我们的目的是为了得到不含噪声的信号()sn，也称为期望信号，若滤波系统的单位脉冲响应为()hn（如图2.1.2所示），系统的期望输出用()dyn表示，()dyn应等于信号的真值()sn；系统的实际输出用()yn表示，()yn是()sn的逼近或估计，用公式表示为()()dynsn，x(n)s(n)v(n)图2.1.1观测信号的组成h(n)x(n)s(n)＋v(n)y(n)图2.1.2信号处理的一般模型ˆ()()ynsn。因此对信号()xn进行处理，可以看成是对期望信号的估计，这样可以将()hn看作是一个估计器，也就是说，信号处理的目的是要得到信号的一个最佳估计。那么，采用不同的最佳准则，估计得到的结果可能不同。所得到的估计，在通信中称为波形估计；在自动控制中，称为动态估计。假若已知(1),(2),,()xnxnxnm，要估计当前及以后时刻的信号值ˆ(),0snNN，这样的估计问题称为预测问题；若已知(),(1),,()xnxnxnm，要估计当前的信号值()sn，称为过滤或滤波；根据过去的观测值(1),(2),,()xnxnxnm，估计过去的信号值(),1snNN，称为平滑或内插。维纳(Wiener)滤波与卡尔曼(Kalman)滤波就是用来解决这样一类从噪声中提取信号的过滤或预测问题，并以估计的结果与信号真值之间的误差的均方值最小作为最佳准则。维纳滤波是在第二次世界大战期间，由于军事的需要由维纳提出的。1950年，伯特和香农给出了当信号的功率谱为有理谱时，由功率谱直接求取维纳滤波器传输函数的设计方法。维纳滤波器的求解，要求知道随机信号的统计分布规律（自相关函数或功率谱密度），得到的结果是封闭公式。采用谱分解的方法求解，简单易行，具有一定的工程实用价值，并且物理概念清楚，但不能实时处理；维纳滤波的最大缺点是仅适用于一维平稳随机信号。这是由于采用频域设计法所造成的，因此人们逐渐转向在时域内直接设计最佳滤波器的方法。早在20世纪40年代，开始有人用状态变量模型来研究随机过程，到60年代初，由于空间技术的发展，为了解决对非平稳、多输入输出随机序列的估计问题，卡尔曼提出了递推最优估计理论。它用状态空间法描述系统，即已知前一状态的估计值和最近一个观测数据，采取递推的算法估计当前的状态值。由于卡尔曼滤波采用递推法，适合于计算机处理，并且可以用来处理多维和非平稳随机信号，现已广泛应用于很多领域，并取得了很好的结果，但它只是维纳滤波的一种算法。卡尔曼滤波已经出现，就受到人们的很大重视，并在实践中不断丰富和完善，其中一个成功的应用是设计运载体的高精度组合导航系统[5]。前面讲到，维纳滤波和卡尔曼滤波都是采用最小方差估计，即采用最小均方误差准则（MMSE，MinimumMeanSquareError）。一般来讲，均方误差准则对大误差抑制能力强，而对小误差不敏感。下面我们首先介绍维纳滤波的时域和复频域的求解方法，然后介绍维纳预测、特别是纯预测问题的求解，最后介绍卡尔曼滤波的内容。2.2维纳滤波器的离散形式——时域解设计维纳滤波器的任务，实际上就是选择()hn，使其输出信号()yn与期望信号()dn误差的均方值为最小，即维纳滤波器是一个均方误差最小准则下的最小滤波器。维纳滤波最初是对连续时间信号以模拟滤波器的形式出现的，尔后才有离散形式，这两种形式解决问题的思路是基本一致的。在这里我们只讨论离散维纳滤波器，模拟维纳滤波器的设计以习题的形式出现，共大家练习。假设滤波系统()hn是一个线性时不变系统，它的单位脉冲响应和输入信号都是复函数，设定()()j(),0,1,2,hnanbnn(0.1.1)2.2.1维纳滤波器时域求解的方法根据线性系统的基本理论，并考虑到系统的因果性，可以得到滤波器的输出()yn,0()()()()()mynxnhnhmxnm(0.1.2)设期望信号为()dn，误差信号()en及其均方值2[|()|]Een分别为()()()()()endnynsnyn(0.1.3)2220[|()|][|()()|]()()()mEenEdnynEdnhmxnm(0.1.4)要使均方误差为最小，须满足2[|()|]0jEenh(0.1.5)这里，jh表示()hj；同理，可以用,jjab分别表示(),()ajbj。由于误差的均方值是一标量，因此（2.2.5）式是一个标量对复函数的求导问题，它等价于2222[|()|][|()|][|()|][|()|]0,0,1,2,jjjjjjjjabEenEenahbhEenEenjjab(0.1.6)记,0,1,2,jjjjjab(0.1.7)则(2.2.6)式可以写为2[|()|]0jEen(0.1.8)将(2.2.8)式展开**2**()()()()[|()|]()()()()jjjjjenenenenEenEenjenenjenabab(0.1.9)又根据(2.2.1)~(2.2.3)式()()jenxnja(0.1.10)()()jenjxnjb(0.1.11)**()()jenxnja(0.1.12)**()()jenjxnjb(0.1.13)将(2.2.10)～(2.2.13)式代入（2.2.9）式，得2*[|()|]2[()()]jEenExnjen(0.1.14)因此[()()]0,0,1,2,Exnjenj(0.1.15)上式说明，均方误差达到最小值的充要条件是误差信号与任一进入估计的输入信号正交，这就是通常所说的正交性原理。它的重要意义在于提供了一个数学方法，用以判断线性滤波系统是否工作于最佳状态。下面计算输出信号与误差信号的互相关函数0[()()][()()()]jEynenEhjxnjen*0()[()()]jhjExnjen(0.1.16)假定滤波器工作于最佳状态，滤波器的输出()optyn与期望信号()dn的误差为()opten，把(2.2.15)式代入上式，得到*opt[()()]0optEynen(0.1.17)这样，在滤波器工作于最佳状态时，期望信号()dn、估计值()yn与误差信号()en的几何关系如图2.2.1所示。图中，用()ne表示误差信号所构成的向量，()optne表示滤波器处于最佳工作状态时，误差信号所构成的向量。同理，()nd表示期望信号所构成的向量，()optny表示滤波器处于最佳工作状态时，输出信号所构成的向量。这样的向量表示方法在第三章介绍最小二乘滤波器时，还要用到。图2.2.1表明滤波器处于最佳工作状态时，估计值加上估计偏差等于期望信号，即optopt()()()nnndye(0.1.18)注意我们所研究的是随机信号，图2.2.1中各矢量的几何表示应理解为相应量的统计平均或者是数学期望。再从能量的角度来看，假定输入信号和期望信号都是零均值，应用正交性原理，则222dopt[||]Eoptye，因此在滤波器处于最佳状态时，估计值的能量总是小于等于期望信号的能量。2.2.2维纳—霍夫方程将（2.2.15）式展开，可以得到eopt(n)d(n)yopt(n)图2.2.1期望信号、估计值与误差信号的几何关系***0()()()()0mExnkdnhmxnm(0.1.19)将输入信号分配进去，得到*0()()(),0,1,2,dxxxmrkhmrmkk对上式两边取共轭，利用相关函数性质：()()yxxyrkrk(或()()xyyxrkrk)，得到0()()()()(),0,1,2,xdxxxxmrkhmrkmhkrkk(0.1.20)(2.2.20)式称为维纳—霍夫(Wiener－Hopf)方程。当()hn是一个长度为M的因果序列(即()hn是一个长度为M的FIR滤波器)时，维纳—霍夫方程表述为10()()()()(),0,1,2,Mxdxxxxmrkhmrkmhkrkk(0.1.21)把k的取值代入(2.2.21)式，得到1212120(0)(1)(1)(0)1(1)(0)(2)(1)1(1)(2)(0)(1)xxxxMxxxdxxxxMxxxdxxxxMxxxdkhrhrhrMrkhrhrhrMrkMhrMhrMhrrM(0.1.22)定义12Mhhhh(0)(1)(1)xdxdxdxdrrrMR(0)(1)(1)(1)(0)(2)(1)(2)(0)xxxxxxxxxxxxxxxxxxxxrrrMrrrMrMrMrR（2.2.22）式可以写成矩阵的形式，即xdxxRRh(0.1.23)对上式求逆，得到1optxxxdhRR(0.1.24)上式表明已知期望信号与观测数据的互相关函数及观测数据的自相关函数时，可以通过矩阵求逆运算，得到维纳滤波器的最佳解。同时可以看到，直接从时域求解因果的维纳滤波器，当选择的滤波器的长度M较大时，计算工作量很大，并且需要计算xxR的逆矩阵，从而要求的存贮量也很大。此外，在具体实现时，滤波器的长度是由实验来确定的，如果想通过增加长度提高逼近的精度，就需要在新M基础上重新进行计算。因此，从时域求解维纳滤波器，并不是一个有效的方法。2.2.3估假定所研究的信号都是零均值的，滤波器为FIR型，长度等于M，将(2.2.2)式和(2.2.3)式代入(2.2.4)式，可以得到1220[|()|][|()|]()[()()]MkEenEdnhkExnkdn111000()[()()]()()[()()]MMMkkihkExnkdnhkhiExnkxni(0.1.25)上式可以进一步化简得到111122*0000[|()|]()()()()()()()MMMMdxdxdxxkkkiEenhkrkhkrkhkhirki2dxdxdxx†††hRRhhRh2111()()dxdxxxdxxxdxxxxxd††RRRhRRRhRR(0.1.26)可以看出，均方误差与滤波器的单位脉冲响应是一个二次函数关系。由于单位脉冲响应()nh为M维向量，因此均方误差是一个超椭圆抛物形曲面，该曲面有极小点存在。当滤波器工作于最佳状态时，均方误差取得最小值。将(2.2.24)式，即1optxxxdhRR，代入(2.2.26)式，得到最小均方误差2212min[|()|]dxdxxxddxdoptEen††RRRRh(0.1.27)下面举一个具体的数值例子，说明维纳滤波器的求解方法。例2.2.1设2()()()ynxnvn，2()vn是一白噪声，方差220.1。期望信号1()xn的信号模型如图2.2.2(a)所示，其中白噪声1()vn的方差210.27，且00.8458b。()xn的信号模型如图2.2.2(b)所示，10.9458b。假定1()vn与2()vn、1()xn与()yn不相关，并