卡尔曼滤波算法总结.docVIP

2015.12.12 void Kalman_Filter(float Gyro,float Accel) { Angle+=(Gyro - Q_bias) * dt; Pdot[0]=Q_angle - PP[0][1] - PP[1][0]; Pdot[1]= - PP[1][1]; Pdot[2]= - PP[1][1]; Pdot[3]=Q_gyro; PP[0][0] += Pdot[0] * dt; PP[0][1] += Pdot[1] * dt; PP[1][0] += Pdot[2] * dt; PP[1][1] += Pdot[3] * dt; Angle_err = Accel - Angle; PCt_0 = C_0 * PP[0][0]; PCt_1 = C_0 * PP[1][0]; E = R_angle + C_0 * PCt_0; K_0 = PCt_0 / E; K_1 = PCt_1 / E; t_0 = PCt_0; t_1 = C_0 * PP[0][1]; PP[0][0] -= K_0 * t_0; PP[0][1] -= K_0 * t_1; PP[1][0] -= K_1 * t_0; PP[1][1] -= K_1 * t_1; Angle += K_0 * Angle_err; Q_bias += K_1 * Angle_err; Gyro_x = Gyro - Q_bias; } 首先是卡尔曼滤波的5个方程： 5个式子比较抽象，现在直接用实例来说： 一、卡尔曼滤波第一个式子 对于角度来说，我们认为此时的角度可以近似认为是上一时刻的角度值加上上一时刻陀螺仪测得的角加速度值乘以时间，因为，角度微分等于时间的微分乘以角速度。但是陀螺仪有个静态漂移（而且还是变化的），静态漂移就是静止了没有角速度然后陀螺仪也会输出一个值，这个值肯定是没有意义的，计算时要把它减去。 由此我们得到了当前角度的预测值Angle Angle=Angle+(Gyro - Q_bias) * dt; 其中等号左边Angle为此时的角度，等号右边Angle为上一时刻的角度，Gyro 为陀螺仪测的角速度的值，dt是两次滤波之间的时间间隔，我们的运行周期是4ms或者6ms。 同时 Q_bias也是一个变化的量。 但是就预测来说认为现在的漂移跟上一时刻是相同的，即 Q_bias=Q_bias 将上面两个式子写成矩阵的形式 得到上式，这个式子对应于卡尔曼滤波的第一个式子 为2维列向量，A为2维方阵，为2维列向量，B 为2维列向量，为 二卡尔曼滤波第二个式子 接着是预测方差阵的预测值，这里首先要给出两个值，一个是漂移的噪声，一个是角度值的噪声，（所谓噪声就是数据的方差值） 的协方差矩阵，即， 因为漂移噪声和角度噪声是相互独立的，则。 又由性质可知即方差，所以得到的矩阵如下 ，这里的两个方差值是开始就给出的常数 程序中的定义如下float Q_angle=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)设为,第一式已知A为 则计算A P(k-1|k-1) A’+Q（就是个矩阵乘法和加法，算算吧）结果如下 很小为了计算简便忽略不计。 于是得到 a,b,c,d分别和矩阵的P[0][0],P[0][1],P[1][0],P[1][1] 计算过程转化为如下程序，代换即可 Pdot[0]=Q_angle - PP[0][1] - PP[1][0]; Pdot[1]= - PP[1][1]; Pdot[2]= - PP[1][1];/ Pdot[3]=Q_gyro; PP[0][0] += Pdot[0] * dt; PP[0][1] += Pdot[1] * dt; PP[1][0] += Pdot[2] * dt; PP[1][1] += Pdot[3] * dt; 三，这里是卡尔曼滤波的第三个式子 Kg(k)= P(k|k-1) H’ / (H P(k|k-1) H’ + R) ……… (3)//计算卡尔曼增益 即计算卡尔曼增益，这是个二维向量设为，这里的 = 为由此kg= P(K|K-1)+R，这里又有一个常数R，程序中的定义如下 float R_angle=0