第1页 / 共7页
第2页 / 共7页
第3页 / 共7页
第4页 / 共7页
第5页 / 共7页
第6页 / 共7页
第7页 / 共7页
卡尔曼滤波算法入门以及C实现.doc
卡尔曼滤波算法及 C 语言实现
摘要：本文着重讨论了卡尔曼滤波器的原理，典型算法以及应用领域。清晰地阐述了
kalman filter 在信息估计方面的最优性能。着重介绍简单 kalman filter algorithm 的编程，使
用 kalman filter 的经典 5 个体现最优化递归公式来编程。通过 c 语言编写程序实现 kalman
filter 的最优估计能力。
关键词：kalman filter；最优估计；C 语言
1 引言
Kalman Filter 是一个高效的递归滤波器，它可以实现从一系列的噪声测量中，估计动态
系统的状态。起源于 Rudolf Emil Kalman 在 1960 年的博士论文和发表的论文《A New
Approach to Linear Eiltering and Prediction Problems》（《线性滤波与预测问题的新方法》）。并
且最先在阿波罗登月计划轨迹预测上应用成功，此后 kalman filter 取得重大发展和完善。它
的广泛应用已经超过 30 年，包括机器人导航，控制。传感器数据融合甚至在军事方面的雷
达系统以及导弹追踪等等，近年来更被广泛应用于计算机图像处理，例如头脸识别，图像分
割，图像边缘检测等等。
2 kalman filter 最 优 化 递 归 估 计
Kalman filter 是一个“optimal recursive data processing algorithm（最优化递归数据处理
方法）”。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的方法。而 kalman filter
最为核心的内容是体现它最优化估计和递归特点的 5 条公式。举一个例子来详细说明 5 条公
式的物理意义。
假设我们要研究的对象是某一个房间的温度信号。对于室温来说，一分钟内或一小段时
间内的值是基本上不变的或者变化范围很小。也就是说 1t 时刻的温度 1T 和 2t 时刻的温度 2T
基本不变，即
2 T
T 。在这个过程中，因为毕竟温度还是有所改变的，设有几度的偏差。
1
我们把这几度的偏差看成是高斯白噪声 )(tw ，也就是说
[
twE
)]
(
0
，
2
twD
)]
(
[
。除此
之外我们在用一个温度计来实时测量房间的温度值 Z ，但由于量具本身的误差，所测得的
温度值也是不准确的，也会和实际值偏差几度，把这几度的偏差看成是测量噪声 )(tv 。即满
足
([
tvE
)]
0
，
2
tvD
1
([
)]
。
此时我们对于这个房间的温度就得到了两个数据。一个是你根据经验得到的经验值
T ，一个是从温度计上得到的测量值 Z ，以及各自引入的高斯白噪声。下面就具体讲
2 T
1
1
解 kalman filter 来估计房间温度的原理与步骤。
要估计 K 时刻的实际温度值，首先要根据 K-1 时刻的温度值预测 K 时刻的温度，按照之
前 我 们 讨 论 的
2 T
T ， 若 k-1 时 刻 的 温 度 值 是
1
Tk
23
1
0
C
， 那 么 预 测 此 时 的
T
k
T
1
k
23
0
C
，假如该值的噪声是
05)(
kw
C
，5°是这样得到的，若果 k-1 时刻估算
出的最优温度值的噪声是
)(
kw
' 3
0
C
，预测的噪声是
)(
kv
' 4
0
C
，所以总体的噪声为
)(
kw
2'
)(
kw
)(
kv
2'
0
5
C
。此时再从温度计上得到 K 时刻的温度值为
T
kz
25
o
C
，设
该测量值的噪声是 C04 。
现在发现问题了，在 k 时刻我们就有了两个温度值
Tk
023
C
和
T
kz
25
o
C
，要信那个
呢，简单的求平均已经不能满足精度的要求了。我们可以用他们的协方差 covariance 来判断。
协方差本身就能体现两个信号的相关性，通过它就能判断到底真值更逼近于预测值还是测量
值。引入 kalman gain（ gk ），有公式计算 gk ，
2
kg
)(
kw
2
/(
)(
kw
2
2
'
5/(5))(
kv
2
2
2
)4
……（1）
所以 gk =0.78。我们可以估算出 K 时刻的实际温度值是，
TT
k
k
g
(
T
kz
T
k
)
23
78.0
25(
)23
56.24
o
C
……（2）
可以看出这个值接近于温度计测量到的值，所以估算出的最优温度值偏向温度计的值。
这时我们已经得到了 K 时刻的最优温度值，接下来估计 K+1 时刻的最优温度值。既然
kalman filter 是一个最优化的递归处理方法，那么递归就体现在该算法的一个核心参数 gk
上，由公式（1） gk 的算法可知每次计算时的 gk 是不一样的。这样我们要估计 K+1 时刻的
最优温度值，就得先算出 K 时刻的 gk ，然后才能利用公式（2）估计 K+1 时刻的最优温度值。
由此可以看出我们只需知道初始时刻的值和它所对应的协方差以及测量值，就可以进行
kalman 估计了。
3 Kalman Filter Algorithm
首先以一个离散控制过程为例讨论 kalman filter algorithm。该系统可用一个线性微分方
程来描述。
)(
kX
(
kXA
)1
)(
kWkUB
)(
……（3）
)(
kZ
)(
kXH
)(
kV
……（4）
2
（3）式和（4）式中，
)(kX 是 K 时刻的系统状态，
)(kU 是 K 时刻对系统的控制量，
A 和 B 是系统参数，对于多模型系统，它们为矩阵。
)(kZ 是 K 时刻的测量值，H 是测量系
统的参数，对于多测量系统，H 为矩阵。
)(kW 和 )(kV 分别表示系统和测量过程中的噪声，
使用 kalman filter 估计时，我们认为噪声满足高斯白噪声模型，设
)(kW 和 )(kV 的 covariance
分别为 Q 和 R。
讨论 kalman filter algorithm 的 5 个经典核心公式。
第一步，预测现在的状态：
|
kkX
(
)1
(
kXA
|1
k
)1
)(
kUB
……（5）
)1
是利用上一状态预测的结果，
(
kX
|1
k
)1
是上一时刻的最优
式（5）中
kkX
(
|
预测值，
)(kU 为现在状态的控制量，如果没有，可以为 0。
经过公式（5）后系统结果已经更新了，对应于
用 P 表示 covariance，
kkX
(
|
)1
的 covariance 还没有更新，
|
(
kkP
)1
(
kPA
|1
k
)1
T
QA
……（6）
式（6）中
(
kkP
|
)1
是
kkX
(
|
)1
对应的 covariance，
(
kP
|1
k
)1
是
(
kX
|1
k
)1
对应的 covariance， TA 是 A 的转置矩阵。Q 是系统的噪声，（5）和（6）式便是 kalman filter
中的前两个公式。对系统的预测。有了系统的预测，接下来就要参考测量值进行估计了。
|
kkX
(
)
|
kkX
(
)1
kk
g
()(
|
kkXHkZ
)(
(
))1
……（7）
由上面分析可知为了实现递归，每次的 gk 都是实时更新的。
)(
kk
g
(
KkP
|
)1
H
T
/(
|
kkPH
(
)1
H
T
R
)
……（8）
|
(
kkP
)
1(
|
kkPHkk
g
)(
(
)
)1
……（9）
这样每次
|
(
kkP
)
和
)(kkg 都需要前一时刻的值来更新，递归的估计下去。（5）~（9）
式便是 kalman filter algorithm 的五条核心公式。
4 利 用 C 语 言 编 程 实 现 Kalman Filter Algorithm
要求是给定一个固定量，然后由测量值来使用 kalman filter 估计系统真实值。
为了编程简单，我将（5）式中的 A=1，
)(kU
=0，（5）式改写为下面的形式，
3
|
kkX
(
)1
(
kX
|1
k
)1
……（10）
式（6）改写为，
|
(
kkP
)1
(
kP
|1
k
)1
Q
……（11）
再令 H=1，式（7），（8），（9）可改写为，
|
kkX
(
)
|
kkX
(
)1
kk
g
()(
)(
kZ
|
kkX
(
))1
……（12）
)(
kkg
(
KkP
|
/()1
|
(
kkP
)1
R
)
(
|
kkP
)
1(
(
kk
g
))
(
|
kkP
)1
……（13）
……（14）
使用 C 语言编程实现（核心算法）。
//x_last=x(k-1|k-1),x_mid=x(k|k-1)
//p_mid=p(k|k-1),p_last=p(k-1|k-1),Q=噪声
x_mid=x_last;
p_mid=p_last+Q;
kg=p_mid/(p_mid+R); //kg 为 kalman filter，R 为噪声
z_measure=z_real+frand()*0.03;//测量值
x_now=x_mid+kg*(z_measure-x_mid);//估计出的最优值
p_now=(1-kg)*p_mid;//最优值对应的 covariance
p_last = p_now;
x_last = x_now;
//更新 covariance 值
//更新系统状态值
5 算 法 测 试
为了测试 kalman filter algorithm，我设计了一个简单实验，来验证 kalman filter 的优良
性。程序中给定一个初值，然后给定一组测量值，验证 kalman filter 估值的准确性。
根据 kalman filter algorithm，我们需要给定系统初始值 x_last，系统噪声 Q 和测量噪声 R，
以及初始值所对应的协方差 P_last。为了验证优劣性，还需要给定真实值 z_real 来计算 kalman
filter 误差 error_kalman 以及测量误差 error_measure 以及它们在有限次的计算中的累积误
差，累积 kalman 误差 sumerror_kalman 和累积测量误差 sumerror_measure。
实验中给定 x_last=0，p_last=0,Q=0.018,R=0.0542.实验中可以通过适当改变 Q 和 R 来获
得更好的估计结果。也可以改变 p_last 和 x_last 的值，由于 kalman filter 是对协方差的递归
算法来估计信号数据的，所以 p_last 对算法结果的影响很大，图 3 就说明了这一情况，由于
在初始时就有协方差，所以在运行过程中算法累积误差相比初始时没有误差的就比较大。
给定值为 z_real=0.56 时运行结果如图 1 所示：
4
给定值 z_real=3.32 时的运行结果如图 2
图 1 真值为 0.56 的运行结果
图 2 真值为 3.32 的运行结果
图 3 为 Q，R 不变，p_last=0.02，x_last=0，z_real=0.56 时的测试结果。通过和前两次结
果比较发现 p_last 对估计结果影响较大，分析可知这种现象是符合 kalman filter 的，通过改
变 Q 和 R 的值也能改善算法的性能，但是实际操作中我们并不能控制这两个量。
5
图 3 改变 p_last 的测试结果
6 结论
本文通过对 kalman filter algorithm 的深入探讨，对 kalman filter 有了更深刻的认识，理
解了核心的 5 条公式的物理意义，以及 kalman filter 的思想，并通过理解算法编程实践，验
证了 kalman filter 在数据处理方面的优良性能。
并且通过实验结果分析了 kalman filter algorithm 的本质对每次估计产生的协方差递归结
合当前测量值来估计系统当前的最佳状态。如要改善算法的性能就必须要尽可能的减小系统
噪声和测量噪声，优化程序，减小估计的协方差。
参 考 文 献
[1]谭浩强.C 程序设计(第三版)[M].北京:清华大学出版社,2005,91~130.
[2]崔平远,黄晓瑞.基于联合卡尔曼滤波的多传感器信息融合算法及其应用[J].电机与
控制学报,2001,9(5): 204-207.
[3]党宏社,韩崇昭,段战胜.基于多卡尔曼滤波器的自适应传感器融合[J].系统工程与
电子技术,2004,5(26):311-313.
[4]文贡坚,王润生. 一种稳健的直线提取算法[J].软件学报,2001,11(11):1660-1665.
6
附录：源程序
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
double frand()
{
return 2*((rand()/(double)RAND_MAX) - 0.5);
//随机噪声}
void main()
{
float x_last=0;
float p_last=0.02;
float Q=0.018;
float R=0.542;
float kg;
float x_mid;
float x_now;
float p_mid;
float p_now;
float z_real=0.56;//0.56
float z_measure;
float sumerror_kalman=0;
float sumerror_measure=0;
int i;
x_last=z_real+frand()*0.03;
x_mid=x_last;
for(i=0;i<20;i++)
{ x_mid=x_last;
//x_last=x(k-1|k-1),x_mid=x(k|k-1)
p_mid=p_last+Q; //p_mid=p(k|k-1),p_last=p(k-1|k-1),Q=噪声
kg=p_mid/(p_mid+R); //kg 为 kalman filter，R 为噪声
z_measure=z_real+frand()*0.03;//测量值
x_now=x_mid+kg*(z_measure-x_mid);//估计出的最优值
p_now=(1-kg)*p_mid;//最优值对应的 covariance
printf("Real
position: %6.3f \n",z_real); //显示真值
printf("Mesaured position: %6.3f [diff:%.3f]\n",z_measure,fabs(z_real-z_measure));
//显示测量值以及真值与测量值之间的误差
printf("Kalman
position: %6.3f [diff:%.3f]\n",x_now,fabs(z_real - x_now)); //显示
kalman 估计值以及真值和卡尔曼估计值的误差
sumerror_kalman += fabs(z_real - x_now);
//kalman 估计值的累积误差
sumerror_measure += fabs(z_real-z_measure); //真值与测量值的累积误差
p_last = p_now; //更新 covariance 值
x_last = x_now; //更新系统状态值
}
printf("总体测量误差
: %f\n",sumerror_measure); //输出测量累积误差
printf("总体卡尔曼滤波误差: %f\n",sumerror_kalman);
//输出 kalman 累积误差
printf("卡尔曼误差所占比例: %d%% \n",100-(int)((sumerror_kalman/sumerror_measure)*100));
}
7
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
