卡尔曼滤波器及matlab代码

实用标准文档文案大全信息融合大作业——维纳最速下降法滤波器，卡尔曼滤波器设计及Matlab仿真时间：2010-12-6专业：信息工程班级：09030702学号：2007302171姓名：马志强实用标准文档文案大全1.滤波问题浅谈估计器或滤波器这一术语通常用来称呼一个系统，设计这样的系统是为了从含有噪声的数据中提取人们感兴趣的，接近规定质量的信息。由于这样一个宽目标，估计理论应用于诸如通信、雷达、声纳、导航、地震学、生物医学工程、金融工程等众多不同的领域。例如，考虑一个数字通信系统，其基本形式由发射机、信道和接收机连接组成。发射机的作用是把数字源(例如计算机)产生的0、1符号序列组成的消息信号变换成为适合于信道上传送的波形。而由于符号间干扰和噪声的存在，信道输出端收到的信号是含有噪声的或失真的发送信号。接收机的作用是，操作接收信号并把原消息信号的一个可靠估值传递给系统输出端的某个用户。随着通信系统复杂度的提高，对原消息信号的还原成为通信系统中最为重要的环节，而噪声是接收端需要排除的最主要的干扰，人们也设计出了针对各种不同条件应用的滤波器，其中最速下降算法是一种古老的最优化技术，而卡尔曼滤波器随着应用条件的精简成为了普适性的高效滤波器。2．维纳最速下降算法滤波器2.1最速下降算法的基本思想考虑一个代价函数𝐽(𝑤)，它是某个未知向量𝑤的连续可微分函数。函数𝐽(𝑤)将𝑤的元素映射为实数。这里，我们要寻找一个最优解𝑤。使它满足如下条件𝐽(𝒘0)≤𝐽(𝒘)(2.1)这也是无约束最优化的数学表示。特别适合于自适应滤波的一类无约束最优化算法基于局部迭代下降的算法：从某一初始猜想𝑤(0)出发，产生一系列权向量𝑤(1),𝑤(2),⋯，使得代价函数𝐽(𝑤)在算法的每一次迭代都是下降的，即𝐽(𝒘(𝑛+1))𝐽(𝒘(𝑛))其中𝑤(𝑛)是权向量的过去值，而𝑤(𝑛+1)是其更新值。我们希望算法最终收敛到最优值𝑤0。迭代下降的一种简单形式是最速下降法，该方法是沿最速下降方向连续调整权向量。为方便起见，我们将梯度向量实用标准文档文案大全表示为𝒈=𝜵𝑱(𝒘)=𝝏𝑱(𝒘)𝝏𝒘(2.2)因此，最速下降法可以表示为𝒘(𝑛+1)=𝒘(𝑛)−12𝜇𝒈(𝑛)(2.3)其中𝑛代表进程，𝜇是正常数，称为步长参数，1/2因子的引入是为了数学上处理方便。在从𝑛到𝑛+1的迭代中，权向量的调整量为𝛿𝒘(𝑛)=𝒘(𝑛+1)−𝒘(𝑛)=−12𝜇𝒈(𝑛)(2.4)为了证明最速下降算法满足式(2.1)，在𝒘(𝑛)处进行一阶泰勒展开，得到𝐽(𝒘(𝑛+1))≈𝐽(𝒘(𝑛))+𝒈𝐻(𝑛)𝛿𝒘(𝑛)(2.5)此式对于𝜇较小时是成立的。在式(2.4)中设𝒘为负值向量，因而梯度向量𝒈也为负值向量，所以使用埃尔米特转置。将式(2.4)用到式(2.5)中，得到𝐽(𝒘(𝑛+1))≃𝐽(𝒘(𝑛))−12𝜇‖𝒈(𝑛)‖2此式表明当𝜇为正数时，𝐽(𝒘(𝑛+1))𝐽(𝒘(𝑛))。因此，随着𝑛的增加，代价函数𝐽(𝑛)减小，当𝑛=∞时，代价函数趋于最小值J∞。2.2最速下降算法应用于维纳滤波器考虑一个横向滤波器，其抽头输入为𝑢(𝑛),𝑢(𝑛−1),⋯,𝑢(𝑛−𝑀+1)，对应的抽头权值为𝑤0(𝑛),𝑤1(𝑛),⋯,𝑤𝑀−1(𝑛)。抽头输入是来自零均值、相关矩阵为𝑅的广义平稳随机过程的抽样值。除了这些输入外，滤波器还要一个期望响应𝑑(𝑛)，以便为最优滤波提供一个参考。在时刻𝑛抽头输入向量表示为𝒖(𝑛)，滤波器输出端期望响应的估计值为𝑑̂(𝑛|𝑈𝑛),其中𝑈𝑛是由抽头输𝑢(𝑛),𝑢(𝑛−1),⋯,𝑢(𝑛−𝑀+1)所张成的空间。空过比较期望响应𝑑(𝑛)及其估计值，可以得到一个估计实用标准文档文案大全误差𝑒(𝑛)，即𝑒(𝑛)=𝑑(𝑛)−𝑑̂(𝑛|𝑈𝑛)=𝑑(𝑛)−𝒘𝐻(𝑛)𝒖(𝑛)(2.6)这里𝒘𝐻(𝑛)𝒖(𝑛)是抽头权向量𝒘(𝑛)与抽头输入向量𝒖(𝑛)的内积。𝒘(𝑛)可以进一步表示为𝒘(𝑛)=[𝑤0(𝑛),𝑤1(𝑛),⋯,𝑤𝑀−1(𝑛)]𝑇同样，抽头输入向量𝒖(𝑛)可表示为𝒖(𝑛)=[𝑢(𝑛),𝑢(𝑛−1),⋯,𝑢(𝑛−𝑀+1)]𝑇如果抽头输入向量𝒖(𝑛)和期望响应𝑑(𝑛)是联合平稳的，此时均方误差或者在时刻𝑛的代价函数𝐽(𝑛)是抽头权向量的二次函数，于是可以得到𝐽(𝑛)=𝜎𝑑2−𝒘𝐻(𝑛)𝒑−𝒑𝑯𝒘(𝑛)+𝒘𝐻(𝑛)𝑹𝒘(𝑛)(2.7)其中，𝜎𝑑2为目标函数𝑑(𝑛)的方差，𝒑抽头输入向量𝒖(𝑛)与期望响应𝑑(𝑛)的互相关向量，及𝑹为抽头输入向量𝒖(𝑛)的相关矩阵。从而梯度向量可以写为𝜵𝑱(𝑛)=[𝜕𝐽(𝑛)𝜕𝑎0(𝑛)+𝑗𝜕𝐽(𝑛)𝜕𝑏0(𝑛)𝜕𝐽(𝑛)𝜕𝑎1(𝑛)𝜕𝐽(𝑛)𝜕𝑎𝑀−1(𝑛)+⋮+𝑗𝜕𝐽(𝑛)𝜕𝑏1(𝑛)𝜕𝐽(𝑛)𝜕𝑏𝑀−1(𝑛)]=−2𝒑+2𝑹𝒘(𝑛)(2.8)其中在列向量中𝜕𝐽(𝑛)𝜕𝑎𝑘(𝑛)和𝜕𝐽(𝑛)𝜕𝑏𝑘(𝑛)分别是代价函数𝑱(𝑛)对应第𝑘个抽头权值𝑤𝑘(𝑛)的实部𝑎𝑘(𝑛)和虚部𝑏𝑘(𝑛)的偏导数。对最速下降算法应用而言，假设式(2.8)中相关矩阵𝑹和互相关向量𝒑已知，则对于给定的抽头权向量𝒘(𝑛+1)为𝒘(𝑛+1)=𝒘(𝑛)+𝜇[𝒑−𝑹𝒘(𝑛)(2.9)它描述了为那滤波中最速下降法的数学表达式。3.卡尔曼滤波器3.1卡尔曼滤波器的基本思想卡尔曼滤波器是用状态空间概念描述其数学公式的，另外新颖的特点是，实用标准文档文案大全他的解递归运算，可以不加修改地应用于平稳和非平稳环境。尤其是，其状态的每一次更新估计都由前一次估计和新的输入数据计算得到，因此只需存储前一次估计。除了不需要存储过去的所有观测数据外，卡尔曼滤波计算比直接根据滤波过程中每一步所有过去数据进行估值的方法都更加有效。𝑣1(𝑛)𝒙(𝑛+1)𝒙(𝑛)y(𝑛)𝒗𝟐(𝑛)图3.1线性动态离散时间系统的信号流图表示“状态”的概念是这种表示的基础。状态向量，简单地说状态，定义为数据的最小集合，这组数据足以唯一地描述系统的自然动态行为。换句话说，状态由预测系统未来特性时所素要的，与系统的过去行为有关的最少的数据组成。典型地，比较有代表性的情况是，状态𝑥(𝑛)是未知的。为了估计它，我们使用一组观测数据，在途中用向量𝒚(𝑛)表示。y(𝑛)成为观测向量或者简称观测值，并假设它是𝑁维的。在数学上，图3.1表示的信号流图隐含着一下两个方程：(1)过程方程𝒙(𝑛+1)=𝑭(𝑛+1,𝑛)𝒙(𝑛)+𝒗1(𝑛)(3.1)式中，M×1向量𝒗1(𝑛)表示噪声过程，可建模为零均值的白噪声过程，且其相关矩阵定义为𝑬[𝒗1(𝑛)𝒗𝟏𝑯(𝑘)]={𝑸𝟏(𝑛)𝑛=𝑘𝑶𝑛≠𝑘(2)测量方程𝒚(𝑛)=𝑪(𝑛)𝒙(𝑛)+𝒗𝟐(𝑛)(3.2)其中𝑪(𝑛)是已知的N×M测量矩阵。N×1向量𝒗𝟐(𝑛)称为测量噪声，建模+z−1𝐈+F(n+1,n)𝑪(𝑛)实用标准文档文案大全为零均值的白噪声过程，其相关矩阵为𝑬[𝒗2(𝑛)𝒗𝟐𝑯(𝑘)]={𝑸𝟐(𝑛)𝑛=𝑘𝑶𝑛≠𝑘(3.3)测量方程(3.2)确立了可测系统输出𝒚(𝑛)与状态𝒙(𝑛)之间的关系，如图3.1所示。3.2新息过程为了求解卡尔曼滤波问题，我们将应用基于新息过程的方法。根据之前所述，用向量𝒚̂(𝑛|𝑦𝑛−1)表示𝑛=1时刻到𝑛−1时刻所有观测数据过去值给定的情况下，你时刻观测数据𝒚(𝑛)的最小均方估计。过去的值用观测值𝒚(1),𝒚(2),⋯,𝒚(𝑛−1)表示，他们张成的向量空间用𝑦𝑛−1表示。从而可以定义新息过程如下：𝜶(𝑛)=𝒚(𝑛)−𝒚̂(𝑛|𝑦𝑛−1)(3.4)其中M×1向量𝜶(𝑛)表示观测数据𝒚(𝑛)的新息。3.3应用新息过程进行状态估计下面，我们根据信息过程导出状态𝒙(𝑖)的最小均方估计。根据推导，这个估计可以表示成为新息过程𝜶(1),𝜶(2),⋯,𝜶(𝑛)序列的线性组合，即𝒙̂(𝑖|𝑦𝑛)=∑𝑩𝑖(𝑘)𝜶(𝑘)𝑛𝑘=1(3.5)其中{𝑩𝑖(𝑘)}𝑘=1𝑛是一组待定的𝑀×𝑁矩阵。根据正交性原理，预测状态误差向量与新息过程正交，即𝐸[𝜺(𝑖,𝑛)𝜶𝐻(𝑚)]=𝐸{[𝒙(𝑖)−𝒙̂[𝑖|𝑦𝑛]]𝜶𝑯(𝑚)}=𝑶(3.6)将式(3.5)代入式(3.6)，并利用新息过程的正交性质，即得𝐸[𝒙(𝑖)𝜶𝑯(𝑚)]=𝑩𝒊(𝑚)𝐸[𝜶(𝑚)𝜶𝑯(𝑚)]=𝑩𝒊(𝑚)𝑹(𝑚)(3.7)因此，式(3.7)两边同时右乘逆矩阵𝑹−1(𝑚)，可得𝑩𝒊(𝑚)的表达式为𝑩𝒊(𝑚)=𝐸[𝒙(𝑖)𝜶𝑯(𝑚)]𝑹−1(𝑚)实用标准文档文案大全(3.8)最后，将式(3.8)带入式(3.5)，可得最小军方差估计𝒙̂(𝑖|𝑦𝑛)=∑𝐸[𝒙(𝑖)𝜶𝑯(𝑘)]𝑹−1(𝑘)𝜶(𝑘)𝑛𝑘=1=∑𝐸[𝒙(𝑖)𝜶𝑯(𝑘)]𝑹−1(𝑘)𝜶(𝑘)𝑛−1𝑘=1+𝐸[𝒙(𝑖)𝜶𝑯(𝑛)]𝑹−1(𝑛)𝜶(𝑛)(3.9)故对于𝑖=𝑛+1，有𝒙̂(𝑛+1|𝑦𝑛)=∑𝐸[𝒙(𝑛+1)𝜶𝑯(𝑘)]𝑹−1(𝑘)𝜶(𝑘)𝑛−1𝑘=1+𝐸[𝒙(𝑛+1)𝜶𝑯(𝑛)]𝑹−1(𝑛)𝜶(𝑛)(3.10)然而，𝑛+1时刻的状态𝒙(𝑛+1)与𝑛时刻的状态𝒙(𝑛)的关系式由式可以推导出对于0≤k≤n，有𝐸[𝒙(𝑛+1)𝜶𝑯(𝑘)]=𝐸{[𝑭(𝑛+1,𝑛)𝒙(𝑛)+𝒗1(𝑛)]𝜶𝐻(𝑘)}=𝑭(𝑛+1,𝑛)𝐸[𝒙(𝑛)𝜶𝑯(𝑘)](3.11)其中𝜶(𝑘)只与观测数据𝒚(1),𝒚(2),⋯,𝒚(𝑘)有关。因此可知，𝒗1(𝑛)与𝜶(𝑘)彼此正交(其中0≤k≤n)。利用式(3.11)以及当𝑖=𝑛时𝒙̂(𝑖|𝑦𝑛)的计算公式，可将式(3.10)右边的求和项改写为∑𝐸[𝒙(𝑛+1)𝜶𝑯(𝑘)]𝑹−1(𝑘)𝜶(𝑘)𝑛−1𝑘=1=𝐹(𝑛+1,𝑛)∑𝐸[𝒙(𝑛)𝜶𝑯(𝑘)]𝑹−1(𝑘)𝜶(𝑘)𝑛−1𝑘=1=𝑭(𝑛+1,𝑛)𝒙̂(𝑛|𝑦𝑛−1)(3.12)为了进一步讨论，引入如下基本定义。3.4卡尔曼增益定义M×N矩阵𝑮(𝑛)=𝐸[𝒙(𝑛+1)𝜶𝐻(𝑛)]𝑹−1(𝑘)(3.13)实用标准文档文案大全其中𝐸[𝒙(𝑛+1)𝜶𝐻(𝑛)]是状态向量𝒙(𝑛+1)和新息过程𝜶(𝑛)的互相关矩阵。利用这一定义和式(3.12)的结果，可以将式(3.10)简单重写为𝒙̂(𝑛+1|𝑦𝑛)=𝑭(𝑛+1,𝑛)𝒙̂(𝑛|𝑦𝑛−1)+𝑮(𝑛)𝜶(𝑛)(3.14)式(3.14)具有明确的物理意义。它标明：线性动态系统状态的最小均方估计𝒙̂(𝑛+1|𝑦𝑛)可以由前一个估计𝒙̂(𝑛|𝑦𝑛−1)求得。为了表示对卡尔曼开创性贡献的认可，将矩阵𝑮(𝑛)称为卡尔曼增益。现在剩下唯一要解决的问题是，怎样以一种便于计算的形式来表示卡尔曼增益