Kalman滤波器的基本原理及仿真

Kalman滤波器的基本原理及仿真摘要：Kalman滤波是对线性最小均方误差滤波的另一种处理方法，实际是维纳滤波的一种递推算法。它采用的递推算法利用了前一时刻的估计值和新的观测值，大大提高了处理的实时性，同时也能自动跟踪随机信号统计特性的非平稳变化，对于解决很大部分的问题，他是最优，效率最高甚至是最有用的，因此得到了广泛的应用。Kalman滤波的应用包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。关键字：Kalman滤波 线性最小均方误差滤波 估计值 观测值一、Kalman滤波器的提出Kalman滤波器是源于匈牙利数学家Rudolf Emil Kalman的博士论文和1960年发表的论文A New Approach to Linear Filtering and Prediction Problems（线性滤波与预测问题的新方法）。在信号处理，通信和现代控制系统中，需要对一个随机动态系统的状态进行估计，由一个测量装置对系统状态进行测量，通过记录的测量值对状态进行最优估计例如：对一个一阶AR模型 的输出状态进行估计。观测方程是 是测量引入的白噪声，通过各的值估计。这类最优估计问题称为卡尔曼滤波。二、Kalman滤波器的基本思想 利用观测数据对状态变量的预测估计进行修正，以得到状态变量的最优估计，即 最优估计预测估计修正 三、Kalman滤波器的特点（1）算法是递推的，时域内设计滤波器，适用于多维随机过程的估计；（2）用递推法计算，不需要知道全部过去的值。用状态方程描述状态变量的动态变化规律，因此，信号可以是平稳的，也可以是非平稳的；（3）误差准则仍为均方误差最小准则。四、有关Kalman滤波器的探索设我们要研究的对象是一个房间的温度。根据经验判断，这个房间的温度是恒定的，也就是下一分钟的温度等于现在这一分钟的温度（假设我们用一分钟来做时间单位）。假设你对你的经验不是100%的相信，可能会有上下偏差几度。我们把这些偏差看成是高斯白噪声（White Gaussian Noise），也就是这些偏差跟前后时间是没有关系的而且符合高斯分配（Gaussian Distribution）。另外，我们在房间里放一个温度计，但是这个温度计也不准确的，测量值会比实际值偏差。我们也把这些偏差看成是高斯白噪声。 在对于某一分钟我们有两个有关于该房间的温度值：你根据经验的预测值（系统的预测值）和温度计的值（测量值）。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。 现在要估算k时刻的是实际温度值。首先你要根据k-1时刻的温度值，来预测k时刻的温度。因为你相信温度是恒定的，所以你会得到k时刻的温度预测值是跟k-1时刻一样的，假设是23度，同时该值的高斯噪声的偏差是5度（5是这样得到的：如果k-1时刻估算出的最优温度值的偏差是3，你对自己预测的不确定度是4度，他们平方相加再开方，就是5）。然后，你从温度计那里得到了k时刻的温度值，假设是25度，同时该值的偏差是4度。由于我们用于估算k时刻的实际温度有两个温度值，分别是23度和25度。究竟实际温度是多少呢？相信自己还是相信温度计呢？究竟相信谁多一点，我们可以用他们的covariance来判断。因为Kg2=52/(52+42)，所以Kg=0.78，我们可以估算出k时刻的实际温度值是：23+0.78*(25-23)=24.56度。可以看出，因为温度计的covariance比较小（比较相信温度计），所以估算出的最优温度值偏向温度计的值。 我们已经得到k时刻的最优温度值了，下一步就是要进入k+1时刻，进行新的最优估算。到现在为止，好像还没看到什么自回归的东西出现。对了，在进入k+1时刻之前，我们还要算出k时刻那个最优值（24.56度）的偏差。算法如下：(1-Kg)*52)0.5=2.35。这里的5就是上面的k时刻你预测的那个23度温度值的偏差，得出的2.35就是进入k+1时刻以后k时刻估算出的最优温度值的偏差（对应于上面的3）。 就是这样，卡尔曼滤波器就不断的把covariance递归，从而估算出最优的温度值。他运行的很快，而且它只保留了上一时刻的covariance。上面的Kg，就是卡尔曼增益（Kalman Gain）。他可以随不同的时刻而改变他自己的值五、Kalman滤波器的的基本理论 5.1、标量（一阶递归模型）卡尔曼滤波器 信号模型和测量模型 为使卡尔曼滤波过程的物理意义明确，采用下列符号：（1）用 代替，表示用n时刻及以前所有数据对所作的最佳线性估计。（2）用代替 ，表示用n-1时刻及以前所有数据对所作的最佳线性估计。 图一对于图一所示的系统，他的迭代形式如下：这是Kalman滤波器的标准形式。5.1.1、物理意义一步预测：是利用n-1时刻及以前的数据对 的预测。第二步预测：是对的最佳估计。新息（Innovation）: 是对 的预测误差，代表中所含的无法预测的信息。选择适当的系数 对新息进行加权，作为对预测值的修正值，修正后得到对信号的最佳估计为：信号n时刻的最佳估计为n时刻的一步预测与新息的加权 和。不同时间的最佳加权系数是不同的，相应的均方误差最小，即现在来求取最佳修正加权系数：滤波误差及功率：令需使误差 和 新息正交。设信号的一步预测误差为：则：得到卡尔曼增益：上式表明：最佳加权系数 ，因预测越不准确，利用 进行的修正就越多。一步预测误差功率：一步预测误差方程： 图二Kalman滤波的结构框图5.2、矢量Kalman滤波器实际应用中常需要根据观测数据同时估计若干个信号，或者估计一个高阶自回归过程，例：对象：q阶AR过程：模型方程和量测方程：若同时估计若干个信号：信号矢量，噪声矢量，参数矩阵如下：5.2.1矢量卡尔曼滤波器的计算公式滤波方程：预测误差方程：增益方程：滤波误差方程：六、卡尔曼滤波器算法流程：七、Kalman滤波器的matlab实现算法7.1、Kalman滤波器的设计离散系统： xn+1=Axn+B(un+wn) (7.1) yn=Cxn (7.2) 其中, wn是在输入端加入的高斯噪声。 状态矩阵参数分别为A = 1.1269-0.49400.1129 1.0000 0 0 0 1.0000 0;B = -0.3832 0.5919 0.5191;C = 1 0 0;我们的目标是设计Kalman滤波器, 在给定输入un和带噪输出测量值 yvn=Cxn+vn的情况下估计系统的输出。 其中, vn是高斯白噪声。 1） 离散Kalman滤波器上述问题的稳态Kalman滤波器方程如下： 测量值修正计算 (7.4)(7.3) 2） 稳态设计 我们可以通过kalman函数设计上述稳态滤波器。 首先定义带噪声的系统模型： xn+1=Axn+Bun+Bwn (状态方程) yn=Cxn (测量方程) 具体的程序代码如下： Plant = ss (A, B B, C, 0, -1, inputname, u w, outputname, y);假设Q=R=1, 下面可以设计离散Kalman滤波器： Q = 1; R = 1; kalmf, L, P, M = kalman(Plant, Q, R); 图三 Kalman滤波器滤波器的功能是在已知输入噪声方差的条件下尽可能消除输出信号中的噪声影响。图四显示了滤波前后的不同输出信号。 下面用程序来比较滤波后输出信号与系统实际信号相对理想输出的误差。 a = A; b = B B 0*B; c = C;C; d = 0 0 0;0 0 1;P = ss(a,b,c,d,-1, inputname, u w v,outputname, yyv);sys = parallel(P, kalmf, 1, 1, , ) % 创建并联系统% 将系统输出yv正反馈到滤波器的输入端, 形成闭环系统SimModel = feedback(sys, 1, 4, 2, 1)% 从I/O列表中删除yvSimModel = SimModel(1 3, 1 2 3)% 产生高斯噪声信号t = 0:100;u = sin(t/5);n = length(t)randn(seed, 0)图四 滤波前后输出信号的比较 w = sqrt(Q)*randn(n, 1);v = sqrt(R)*randn(n, 1);% 系统仿真out, x = lsim(SimModel, w, v, u);y = out(:, 1); % 系统真实（理想）输出响应ye = out(:, 2);% 滤波后的系统输出yv = y + v;% 系统输出的测量值% 比较结果subplot(211), plot(t, y, -, t, ye, -), xlabel(No. of samples), ylabel(Output)title(Kalman filter response)subplot(212), plot(t, y-yv, -., t, y-ye, -), xlabel(No. of samples), ylabel(Error)比较结果如图五所示。图中上面的图形显示的是真实响应y（虚线）和滤波后的输出ye（实线）， 下面的图形比较测量误差（虚线）与估计误差（实线）。 该图表明, 滤波器最大程度地消除了系统输出中的噪声影响。 这可以通过计算误差的协方差进行验证。 MeasErr = y-yv;MeasErrCov = sum(MeasErr.*MeasErr)/length(MeasErr);EstErr = y-ye;EstErrCov = sum(EstErr.*EstErr)/length(EstErr);滤波前误差（测量误差）的协方差为MeasErrCovMeasErrCov = 1.1138而滤波后的误差（估计误差）的协方差仅为 EstErrCovEstErrCov = 0.27227.2完整程序A = 1.1269 -0.4940 0.1129 1.0000 0 0 0 1.0000 0; B = -0.3832 0.5919 0.5191; C = 1 0 0; Plant = ss(A,B B,C,0,-1,inputname,u w,outputname,y);Q = 1; R = 1; kalmf,L,P,M = kalman(Plant,Q,R); kalmf = kalmf(1, :);a = A;b = B B 0*B; c = C;C;d = 0 0 0;0 0 1;P = ss(a,b,c,d,-1,inputname,u w v,outputname,y yv); sys = parallel(P,kalmf,1,1,) ;% Close loop around input #4 and output #2 SimModel = feedback(sys,1,4,2,1) ;% Delete yv from I/O list SimModel = SimModel(1 3,1 2 3) ;SimModelt = 0:100; u = sin(t/5); n = length(t) ;randn(seed,0) ;w = sqrt(Q)*randn(n,1); v = sqrt(R)*randn(n,1); out,x = lsim(SimModel,w,v,u); y = out(:,1); % true response ye = out(:,2); % filtered response yv = y + v; % measured response figure(1); plot(t,y,-,t,ye,-),xlabel(No. of samples), ylabel(Output) title(Kalman filter response) figure(2);plot(t,(y-yv)/y,-,t,(y-ye)/y,-), xlabel(No. of samples), ylabel(Error) M