在求学的道路上还是一定要有求知和专注的精神,之前是涉猎的比较多,导致现在面广而不实,只能从头开始学习,近期在学习卡尔曼滤波,久闻其名,却不知道如何应用。这次根据看到的几篇博客,来自己记录一下,增强记忆,也便于各位朋友来指点:

一、什么是卡尔曼滤波器

 卡尔曼滤波器,是一种“optimal recursive data processing algorithm”方法,最优化自回归数据处理算法,在机器人导航、

控制、传感器信号融合以及导弹等方面得到了广泛应用,近年来在图像处理方面又得到了广泛应用。而我要学习它,则是因为近

段时间以来,接触到定位导航方面的事情。

二、卡尔曼滤波器算法

 首先我们引入一个线性离散控制过程的系统:

 X(K)=A*X(K-1)+B*U(K)+W(K)

 Z(K)=HX(K)+V(K)

上式中,X(K)是系统当前的状态,U(K)是K时刻系统的控制量,Z(K)是K时刻系统的测量值,A、B是系统参数,H是测量系统参数

W和V分别表示过程系统和测量系统的噪声,服从高斯分布,协方差分别为Q和R。

卡尔曼滤波器主要包含以下五个公式:

X(K|K-1)=A*X(K-1|K-1)+B*U(K)................................................................(1)

上式表示由上一时刻的状态预测当前时刻的状态

 P(K|K-1)=A*P(K-1|K-1)*A'+Q....................................................................(2)

 上式中,P(K|K-1)表示X(K|K-1)对应的协方差,P(K-1|K-1)代表X(K-1|K-1)对应的协方差。

X(K|K)=X(K|K-1)+Kg(K)*(Z(K)-H*X(K|K-1))..............................................(3)

上式中,X(K|K)表示由预测值和测量值得到的当前时刻的最优估计值,Kg(K)代表卡尔曼增益。

Kg(K)=P(K|K-1)*H'/[H*P(K|K-1)H'+R]........................................................(4)

 通过式(3)已经得到了当前状态的最优估计,但是为了使得系统能够不断迭代,需要更新得到X(K|K)对应的协方差:

 P(K|K)=(I-Kg(K)*H)*P(K|K-1).........................................................................(5)

至此,卡尔曼滤波器的算法已经完毕,网上还有很多温度或小车的例子,此处不再赘述。

三、Matlab实现

******************************************************************************


*****************************************************************************



对比四条曲线,会发现跟随性能不一样,是由于曲线值太小,噪声相对太大的缘故,将值变大,噪声相对就小了,跟随就好了

%第一条曲线 

 
 
 %t=0.1:0.1:6; 

 
 
 %X=t; 

 
 
 %第二条曲线 

 
 
 %t=1:1:60; 

 
 
 %X=t; 

 
 
 %第三条曲线 

 
 
 %t=1:1:60; 

 
 
 %X=sin(t); 

 
 
 %第四条曲线 

 
 
 t=1:1:60; 

 
 
 X=60* 
 sin 
 (t); 

 
 
  

 

    

 
 
 N=1:60; 

 
 
 plot 
 (N,X, 
 'r' 
 ); 

 
 
 hold  
 on; 

 

    

 
 
   

 
 
 %系统方程:X(k+1)=A*X(k)+w(k)  

 
 
 %观测方程:Z(k)=H*Z(k)+v(k)  

 
 
 %这个是生成高斯噪声的随机数 

 
 
 w= 
 randn 
 (1,60);    

 
 
 v= 
 randn 
 (1,60); 

 

    

 
 
 A=1; 

 
 
 H=1; 

 
 
 X_k(1)=0;                       
 %状态估计初值              

 
 
 P_kk(1)=0;                      
 %P(k/k) 

 
 
 P_k(1)=0;                       
 %P(k/k-1) 这个初始化不需要,就给你们看看变量的对应 

 
 
 Z_k(1)=X_k(1)+w(1) ;           
 %测量值 

 

    

 
 
 R=( 
 std 
 (v)).^2;  

 
 
 Q=( 
 std 
 (w)).^2;  

 

    

 
 
 Kg(1)=P_kk(1)*H 
 '/(H*P_kk(1)*H' 
 +R);         
 %卡尔曼增益 Kg 

 
 
 P_k(1)=A*P_kk(1)*A'+Q ;                    
 %方差预测   P_k/k-1 

 

    

 
 
 for  
 i 
 =2:60 

 
 
  
 X_k( 
 i 
 )=A*X_k( 
 i 
 -1)+Kg( 
 i 
 -1)*(Z_k( 
 i 
 -1)-H*A*X_k( 
 i 
 -1));  
 %这边就直接代入式(1),所以没出现 

 
 
  
 Kg( 
 i 
 )=P_k( 
 i 
 -1)*H 
 '\((H*P_k(i-1)*H' 
 +R)) ; 

 
 
  
 P_k( 
 i 
 )=A*P_kk( 
 i 
 -1)*A'+Q; 

 
 
  
 P_kk( 
 i 
 )=P_k( 
 i 
 -1)-Kg( 
 i 
 )*H*P_kk( 
 i 
 -1) ; 

 
 
  
 Z_k( 
 i 
 )=H*X( 
 i 
 )+w( 
 i 
 ) ; 

 
 
 end  

 
 
 n=1:60; 

 
 
 plot 
 (n,X_k); 

 ************************************************************************************************** 

********************************************************************************************************
clear
N=200;
w(1)=0;
w=randn(1,N)
x(1)=0;
a=1;
for k=2:N;
x(k)=a*x(k-1)+w(k-1);
end


V=randn(1,N);
q1=std(V);
Rvv=q1.^2;
q2=std(x);
Rxx=q2.^2; 
q3=std(w);
Rww=q3.^2;
c=0.2;
Y=c*x+V;

p(1)=0;
s(1)=0;
for t=2:N;
p1(t)=a.^2*p(t-1)+Rww;
b(t)=c*p1(t)/(c.^2*p1(t)+Rvv);
s(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1));
p(t)=p1(t)-c*b(t)*p1(t);
end

t=1:N;
plot(t,s,'r',t,Y,'g',t,x,'b');

函数先抄上,尚未验证,等装了matlab再回来补充验证结果。