|
本节内容需要各位拥有一定的线性代数的计算基础。（最好会用matlab） 根据前面的描述，我们把房间看成一个系统，然后对这个系统建模。当然，我们见的模型不需要非常地精确。我们所知道的这个房间的温度是跟前一时刻的温度相同的，所以A=1。没有控制量，所以U(k)=0。因此得出：
X(k|k-1)=X(k-1|k-1) ……….. (6)
式子（2）可以改成：
P(k|k-1)=P(k-1|k-1) +Q ……… (7)
因为测量的值是温度计的，跟温度直接对应，所以H=1。式子3，4，5可以改成以下：
X(k|k)= X(k|k-1)+Kg(k) (Z(k)-X(k|k-1)) ……… (8)
Kg(k)= P(k|k-1) / (P(k|k-1) + R) ……… (9)
P(k|k)=（1-Kg(k)）P(k|k-1) ……… (10)
现在我们模拟一组测量值作为输入。假设房间的真实温度为25度，我模拟了200个测量值，这些测量值的平均值为25度，但是加入了标准偏差为几度的高斯白噪声。
为了令卡尔曼滤波器开始工作，我们需要告诉卡尔曼两个零时刻的初始值，是X(0|0)和P(0|0)。他们的值不用太在意，随便给一个就可以了，因为随着卡尔曼的工作，X会逐渐的收敛。但是对于P，一般不要取0，因为这样可能会令卡尔曼完全相信你给定的X(0|0)是系统最优的，从而使算法不能收敛。我选了 X(0|0)=1度，P(0|0)=10。
该系统的真实温度为25度。
卡尔曼滤波是以最小均方误差为估计的最佳准则，来寻求一套递推估计的算法，其基本思想是：采用信号与噪声的状态空间模型，利用前一时刻地估计值和现时刻的观测值来更新对状态变量的估计，求出现时刻的估计值。 现设线性时变系统的离散状态防城和观测方程为： X(k) = F(k,k-1)·X(k-1)+T(k,k-1)·U(k-1) Y(k) = H(k)·X(k)+N(k) 其中: X(k)和Y(k)分别是k时刻的状态矢量和观测矢量 F(k,k-1)为状态转移矩阵 U(k)为k时刻动态噪声 T(k,k-1)为系统控制矩阵 H(k)为k时刻观测矩阵 N(k)为k时刻观测噪声
则卡尔曼滤波的算法流程为： 预估计X(k)^= F(k,k-1)·X(k-1)
1. 计算预估计协方差矩阵
C(k)^=F(k,k-1)×C(k)×F(k,k-1)'+T(k,k-1)×Q(k)×T(k,k-1)'
Q(k) = U(k)×U(k)'
2. 计算卡尔曼增益矩阵
K(k) = C(k)^×H(k)'×[H(k)×C(k)^×H(k)'+R(k)]^(-1)
R(k) = N(k)×N(k)'
3. 更新估计
X(k)~=X(k)^+K(k)×[Y(k)-H(k)×X(k)^]
4. 计算更新后估计协防差矩阵
C(k)~ = [I-K(k)×H(k)]×C(k)^×[I-K(k)×H(k)]'+K(k)×R(k)×K(k)'
5. X(k+1) = X(k)~
C(k+1) = C(k)~
重复以上步骤
其c语言实现代码如下： - #include "stdlib.h"
- #include "rinv.c"
- int lman(n,m,k,f,q,r,h,y,x,p,g)
- int n,m,k;
- double f[],q[],r[],h[],y[],x[],p[],g[];
- {
- int i,j,kk,ii,l,jj,js;
- double *e,*a,*b;
- e=malloc(m*m*sizeof(double));
- l=m;
- if (l<n) l=n;
- a=malloc(l*l*sizeof(double));
- b=malloc(l*l*sizeof(double));
- for (i=0; i<=n-1; i++)
- for (j=0; j<=n-1; j++)
- { ii=i*l+j; a[ii]=0.0;
- for (kk=0; kk<=n-1; kk++)
- a[ii]=a[ii]+p[i*n+kk]*f[j*n+kk];
- }
- for (i=0; i<=n-1; i++)
- for (j=0; j<=n-1; j++)
- { ii=i*n+j; p[ii]=q[ii];
- for (kk=0; kk<=n-1; kk++)
- p[ii]=p[ii]+f[i*n+kk]*a[kk*l+j];
- }
- for (ii=2; ii<=k; ii++)
- { for (i=0; i<=n-1; i++)
- for (j=0; j<=m-1; j++)
- { jj=i*l+j; a[jj]=0.0;
- for (kk=0; kk<=n-1; kk++)
- a[jj]=a[jj]+p[i*n+kk]*h[j*n+kk];
- }
- for (i=0; i<=m-1; i++)
- for (j=0; j<=m-1; j++)
- { jj=i*m+j; e[jj]=r[jj];
- for (kk=0; kk<=n-1; kk++)
- e[jj]=e[jj]+h[i*n+kk]*a[kk*l+j];
- }
- js=rinv(e,m);
- if (js==0)
- { free(e); free(a); free(b); return(js);}
- for (i=0; i<=n-1; i++)
- for (j=0; j<=m-1; j++)
- { jj=i*m+j; g[jj]=0.0;
- for (kk=0; kk<=m-1; kk++)
- g[jj]=g[jj]+a[i*l+kk]*e[j*m+kk];
- }
- for (i=0; i<=n-1; i++)
- { jj=(ii-1)*n+i; x[jj]=0.0;
- for (j=0; j<=n-1; j++)
- x[jj]=x[jj]+f[i*n+j]*x[(ii-2)*n+j];
- }
- for (i=0; i<=m-1; i++)
- { jj=i*l; b[jj]=y[(ii-1)*m+i];
- for (j=0; j<=n-1; j++)
- b[jj]=b[jj]-h[i*n+j]*x[(ii-1)*n+j];
- }
- for (i=0; i<=n-1; i++)
- { jj=(ii-1)*n+i;
- for (j=0; j<=m-1; j++)
- x[jj]=x[jj]+g[i*m+j]*b[j*l];
- }
- if (ii<k)
- { for (i=0; i<=n-1; i++)
- for (j=0; j<=n-1; j++)
- { jj=i*l+j; a[jj]=0.0;
- for (kk=0; kk<=m-1; kk++)
- a[jj]=a[jj]-g[i*m+kk]*h[kk*n+j];
- if (i==j) a[jj]=1.0+a[jj];
- }
- for (i=0; i<=n-1; i++)
- for (j=0; j<=n-1; j++)
- { jj=i*l+j; b[jj]=0.0;
卡尔曼滤波的内容简单介绍到这里，比较困难，对于想搞自动控制和四元数还有姿态控制（包括姿态融合）的同志们一定要弄懂，卡尔曼核心就是观测状态和实际状态的问题。希望各位能看明白。如果写的不好，希望指正。
如果各位想要真正弄懂卡尔曼滤波，还需要参考其他的书籍和相关论文，因为这种算法可以详细的写成一本书。
| 
/3 
发表于 2015-9-15 14:46:43
沉默卡
照妖镜
