卡尔曼滤波算法总结

1、2015.12.12void Kalman_Filter(float Gyro,float Accel) Angle+=(Gyro - Q_bias) * dt;Pdot0=Q_angle - PP01 - PP10;Pdot1= - PP11;Pdot2= - PP11;Pdot3=Q_gyro;PP00 += Pdot0 * dt;PP01 += Pdot1 * dt;PP10 += Pdot2 * dt;PP11 += Pdot3 * dt;Angle_err = Accel - Angle;PCt_0 = C_0 * PP00;PCt_1 = C_0 * PP10;E = R_angl

2、e + C_0 * PCt_0;K_0 = PCt_0 / E;K_1 = PCt_1 / E;t_0 = PCt_0; t_1 = C_0 * PP01;PP00 -= K_0 * t_0;PP01 -= K_0 * t_1;PP10 -= K_1 * t_0;PP11 -= K_1 * t_1;Angle+= K_0 * Angle_err;Q_bias += K_1 * Angle_err; Gyro_x = Gyro - Q_bias;首先是卡尔曼滤波的5个方程:P(k | k 1)H 7HP(k|k 1)H ' R)X(k|k 1) Kg (k)(Z(k) HX (k|k(1

3、) 先验估计(2) 协方差矩阵的预测(3) 计算卡尔曼增益1)(4)进行修正(5)更新协方差阵X(k|k 1) AX(k 1|k 1) Bu(k) P(k|k 1) AP(k 1| k 1)A' Q Kg(k)X(k|k)P(k|k) (I Kg(k) H) P(k |k 1)5个式子比较抽象，现在直接用实例来说:一、卡尔曼滤波第一个式子对于角度来说，我们认为此时的角度可以近似认为是上一时刻的角度值加上 上一时刻陀螺仪测得的角加速度值乘以时间，因为 d dt ，角度微分等于时 间的微分乘以角速度。但是陀螺仪有个静态漂移(而且还是变化的)，静态漂移就是静止了没有角速度然后陀螺仪也会输出一

4、个值，这个值肯定是没有意义的， 计算时要把它减去。由此我们得到了当前角度的预测值AngleAn gle=A ngle+(Gyro - Q_bias) * dt;其中等号左边Angle为此时的角度，等号右边 Angle为上一时刻的角度，Gyro 为陀螺仪测的角速度的值，dt是两次滤波之间的时间间隔，我们的运行周期是 4ms或者6ms同时Q_bias也是一个变化的量。但是就预测来说认为现在的漂移跟上一时刻是相同的，即Q_bias=Q_bias将上面两个式子写成矩阵的形式An gledt AngleQ_bias1 Q_biasdtGyro0得到上式，这个式子对应于卡尔曼滤波的第一个式子X(k |k

5、1) AX(k 1| k 1) Bu(k)X k|k 1为2维列向量QX，A为2维方阵dtX k-1|k 1 为 2 维列向量AngleQ_biasB为2维列向量dt0u k 为 Gyro卡尔曼滤波第二个式子接着是预测方差阵的预测值，这里首先要给出两个值，一个是漂移的噪声, 一个是角度值的噪声，(所谓噪声就是数据的方差值)P(k|k 1) AP(k 1| k 1)A' Q这里的Q为向量An gle的协方差矩阵，即Q_biascov(Angle,Angle) cov(Q_bias,Angle) cov(A ngle,Q_bias)cov(Q_bias,Q_biaS)因为漂移噪声和角度噪声

6、是相互独立的，则cov(Angle,Q_bias) 0。又由性质可知cov(x，x) D(x)即方差，所以得到的矩阵如下D(A ngle)00D(Q_bias)，这里的两个方差值是开始就给出的常数程序中的定义如下float Q_an gle=0.001;float Q_gyro=0.003;接着是这一部分 A P(k-1|k- 1) A ',其中的(P( k-1)|P(k-1)为上一时刻的 预测方差阵卡尔曼滤波的目标就是要让这个预测方差阵最小。其中 P(k-1|k-1)设为1式已知A为0dt1则计算A P(k-1|k- 1) A ' +Q (就是个矩阵乘法和加法，算算吧)结果如

7、下D(Angle) b d dta c dt b dt d.(dt)2c d dt2d.(dt)很小为了计算简便忽略不计。于是得到a c dt b dt D(Angle) b d dtc d dtda,b,c,d 分别和矩阵的 P00,P01,P10,P11计算过程转化为如下程序，代换即可Pdot0=Q_a ngle - PP01 - PP10;Pdot1= - PP11;Pdot2= - PP11;/Pdot3=Q_gyro;PP00 += Pdot0 * dt;PP01 += Pdot1 * dt;PP10 += Pdot2 * dt;PP11 += Pdot3 * dt;三，这里是卡尔曼

8、滤波的第三个式子计算卡尔曼增益1 0为由此kg=Kg(k)= P(k|k- 1) H ' / (H P(k |k- 1) H ' + R) (3)/即计算卡尔曼增益，这是个二维向量设为k1 ,这里的P(K|K-1)+R，这里又有一个常数R,程序中的定义如下 float R_an gle=0.5;这个指的是角度测量噪声值，则式子的分母=P00+R_a ngle即程序中的PCt_0 = C_0 * PP00;PCt_1 = C_0 * PP10;E = R_an gle + C_0 * PCt_0; 分子 卩°°P10K。于是求出K1K_0 = PCt_0 /

9、E;K_1 = PCt_1 / E;四，用误差还有卡尔曼增益来修正X(k|k)= X(k|k-1)+Kg(k) (Z(k) - H X(k|k-1) 通过卡尔曼增益进行修正这个矩阵带进去就行了 Z( k)=Accel 注意这个是加速度计算出来的角度Angle_err = Accel - Angle; 对应程序如下Angle += K_0 * Angle_err;Q_bias += K_1 * Angle_err;同时为了 PID 控制还有下次的使用把角速度算出来了Gyro_x = Gyro - Q_bias;五，最后一步对矩阵 P 进行更新，因为下一次滤波时要用到PP00 -= K_0 * t_0;PP01 -= K_0 * t_1;PP10 -= K_1 * t_0;PP11 -= K_1 * t_1;P(k|k)=(I-Kg(k) H ) P(k|k- 1) (5)/跟预测方差阵这个很简单，矩阵带进去算就行了六，总结 卡尔曼滤波一共只需要给很少的初始值量，float Q_angle=0.001;float