以NASA之名: 卡尔曼滤波器

'That's one small step for man,one giant leap for mankind.'Neil Alden Armstron

引言

二十世纪的阿波罗登月计划在人类历史上是浓墨重彩的一笔, 是人类科学发展极其重要的里程碑. 在此计划中, 阿姆斯特朗在月球上说出了上面的一句话,是对此计划最最恰当的注释. 说起来这个计划很''简单'': 送人到月球转一圈,然后再回来. 这么一个'简单'的计划实施起来得有多困难大家心中早有尺度.而卡尔曼滤波器(Kalman filter)就是其中的功臣.

荣耀骑士

卡尔曼滤波(器)是以其主要贡献者 Rudolf Emil Kalman 命名的. 这个滤波器做什么的呢?

现在,时间闪回到 1969年 七月, 当时 Armstrong 正坐在阿波罗11号中,飞向月球. 在此期间掌握飞船的位置及其速度是非常重要的.因为只要偏差超过一定阈值, 飞船的超高速就会将这偏差迅速放大,进而飞船就会偏离预定轨道,最终,极有可能飞向太空就再也回不来了.

那要如何得到飞船在某个时间点的位置与速度呢?

我们知道地面上的科学家在将阿波罗11号发射升空之前肯定会预先计算飞船每个点的位置,或者一定有一个计算精确的计算公式.但, 无论多么精确的公式都无法包含所有因素,不确定性是无法避免的. 而阿波罗上面肯定会有一套复杂的传感系统来测算飞船的位置与速度, 然而, 传感觉器必然地会有偏差. 也就是说我们会从两种不同的途径获得两套参数(位置与速度), 这两套参数一般是不同的,却又都是不准确的.

我们当然会想到:能不能综合地考虑这两套参数,从而获得一套靠谱的参数,来指导飞船航行? 那这两个参数如何综合考虑呢? 在此关键时刻,荣耀骑士 — 卡尔曼滤波器登场了.它的出现就是解决此等棘手问题的.

卡尔曼滤波器是如何工作的呢? 我们知道科学家的公式是不会错的,只不过外界因素干扰,才使其失准的.因此我们以其为基准,用传感器的参数来对基校准. 传感器的参数也有不确定性. 我们的目的是获得更精准的参数, 故需要在传感参数校正公式参数时,首要任务即是在当前所有信息面前最小化不确定性.这就是卡尔曼滤波器的作用. 至于它是如何最小化不确定性的,就要涉及公式了, 你确定要上吗? 如果是,请继续~

卡尔曼滤波器*

卡尔曼滤波(器)基于线性动态系统,是建立在线性模型上的带有白噪声的马氏链. 其与隐马尔可夫模型(HMM)很相像,最主要的不同点在其隐状态是连续的,而HMM是离散的.

模型:
\[
\left\{
\begin{array}\
x_{k} = A_kx_{k-1}+B_k u_k + w_k\\
y_k = C_k x_k + v_k
\end{array}
\right.
\]
其中所有变量均为矩阵(或向量), \(x_k\): k 时刻的状态向量; \(A_k\): 状态转移矩阵; \(B_k\): 控制矩阵,\(u_k\): 控制向量; 一般,将此控制部分当作确定部分而不包含在模型中,因此上式可简化为:
\[
\left\{
\begin{array}\
x_{k} = A_kx_{k-1} + w_k\\
y_k = C_k x_k + v_k
\end{array}
\right.
\]
\(y_k\) : k时刻的观测向量; \(C_k\) : 观测矩阵; \(w_k, v_k\)均为零均值白噪声,方差分别为 \(Q_k, R_k\), 其互不相关,并与初始状态亦不相关:
\[
w_k : E[w_k] = 0,\sigma_k^2 = Q_k;\\
v_k: E[v_k] = 0,\sigma_k^2 = R_k;\\
cov[x_0, w_k] = 0;\\
cov(x_0,v_k) = 0;\\
E[w_kv_k^T] = 0.\\
\]
卡尔曼滤波器即是利用递推方法及状态方程(式1或式2) 寻找最小均方误差状态变量 \(x_k\) 的估计值 \(\hat{x}_k\):
\[
\tilde{x}_k = x_k - \hat{x}_k\\
\min E[\tilde{x}_k\tilde{x}_k^T] \rightarrow \hat{x}_k
\]
当不考虑白噪声时:
\[
\hat{x}'_k = A_k \hat{x}_{k-1}\\
\hat{y}_k = C_k \hat{x}'_k = C_k A_k \hat{x}_{k-1}
\]
其中 \(\hat{x}'_k\) 为 k 时刻状态的估计值(由上一时刻状态经转移矩阵得到), \(\hat{y}_k\) 为测量值的估计值(由状态估计值推出), 而 \(\hat{x}_k\) 现在称为状态校正后的估计值.

观测向量的误差(也即新息):
\[
\tilde{y}_k = y_k - \hat{y}_k
\]
用观测值来校正状态变量:
\[
\begin{array}
\\
\hat{x}_k &=&A_k \hat{x}_{k-1} + H_k (y_k - \hat{y}_k)\\
&=&A_k \hat{x}_{k-1} + H_k (y_k - \ C_k A_k \hat{x}_{k-1})
\end{array}
\]
其中, \(H_k\)称为增益矩阵,即为一新息加权矩阵.

新息(Innovation)\(\tilde{y}_k\)即是在 k 时刻之前没有但在 t 时刻产生的新的新信息. 用代数的语言即是 新息与过去的信息是正交的:
\[
E[\tilde{y}_ky_j^T] = \boldsymbol{O}, \qquad 1\le j<k
\]
而且可以容易的推论出, 新息之间也是正交的:
\[
E[\tilde{y}_i\tilde{y}_j^T] =\boldsymbol{O},\qquad i\ne j, and\quad 1\le i,j\le k
\]
从而可知, 新息\(\tilde{y}_k\)是一个与k 之前时刻的数据不相关,且有白噪声性质的随机过程.

校正后的状态变量估计误差及其均方值 (\(P_k\)) 还有未经校正的状态变量估计误差的均方值 (\(P_k'\)):
\[
\begin{array}\\
\tilde{x}_k &=& x_k - \hat{x}_k\\
P_k &=& E[\tilde{x}_k\tilde{x}_k^T] = E[(x_k - \hat{x}_k)(x_k - \hat{x}_k)^T] \\
P_k' &=& E[(x_k - \hat{x}_k')(x_k - \hat{x}_k')^T]\\
\end{array}
\]
卡尔曼滤波器即是最小化 \(P_k\)(不确定性) 的方式来得到 \(H_k\) ,最终得到 \(\hat{x}_k\) :
\[
\min P_k \rightarrow H_k \rightarrow \hat{x}_k
\]

\[
\begin{array}\\
\tilde{x}_k &=& x_k - \hat{x}_k\\
&=& A_kx_{k-1} + w_k -(A_k \hat{x}_{k-1} + H_k (y_k - \ C_k A_k \hat{x}_{k-1}))\\
&=& A_kx_{k-1} + w_k\\
&&-(A_k \hat{x}_{k-1} + H_k (\ C_k (A_k x_{k-1}+w_k)+ v_k - \ C_k A_k \hat{x}_{k-1}))\\
&=& (A_k -H_kC_kA_k)(x_{k-1} - \hat{x}_{k-1}) + (I - H_kC_k)w_k - H_k v_k\\
\end{array}
\]

偷个懒吧, 接着就是'暴力'推导, 没有什么技术含量, 只看你有没有仔细,否则结果一定会出来.

经过几张\(A_4\) 纸的推导, 你可以得到以下几个重要公式:
\[
\left\{\begin{array}\\
\hat{x}_k &=& A_k \hat{x}_{k-1} + H_k(y_k - C_kA_k\hat{x}_{k-1})\\
H_k &=& P_k'C_k^T(C_kP'_kC_k^T+R_k)^{-1}\\
P_k' &=& A_kP_{k-1}A_k^T + Q_{k-1}\\
P_k &=& (I - H_kC_k)P_k'
\end{array}\right.
\]
这就是卡尔曼滤波器的递推公式了.

参考文献:

  1. 现代信号处理(第二版), 张贤达, 2002, 清华大学出版社.
  2. How a Kalman filter works, in pictures, 2015, bzarg.com.
  3. Kalman filter with Real-Time Applications, 2009,Springer.
  4. Understanding the Basis of the Kalman filter via a simple and intuitive derivation.

时间序列八: 以NASA之名: 卡尔曼滤波器的更多相关文章

  1. 对Kalman(卡尔曼)滤波器的理解

    1.简单介绍(Brief Introduction) 在学习卡尔曼滤波器之前,首先看看为什么叫"卡尔曼". 跟其它著名的理论(比如傅立叶变换.泰勒级数等等)一样.卡尔曼也是一个人的 ...

  2. 对Kalman(卡尔曼)滤波器的理解@@zz

    1.简介(Brief Introduction) 在学习卡尔曼滤波器之前,首先看看为什么叫“卡尔曼”.跟其他著名的理论(例如傅立叶变换,泰勒级数等等)一样,卡尔曼也是一个人的名字,而跟他们不同的是,他 ...

  3. opencv 卡尔曼滤波器例子,自己修改过

    一.卡尔曼滤波器的理论解释 http://blog.csdn.net/lindazhou2005/article/details/1534234(推荐) 二.代码中一些随机数设置函数,在opencv中 ...

  4. 卡尔曼滤波器 Kalman Filter (转载)

    在学习卡尔曼滤波器之前,首先看看为什么叫“卡尔曼”.跟其他著名的理论(例如傅立叶变换,泰勒级数等等)一样,卡尔曼也是一个人的名字,而跟他们不同的是,他是个现代人! 卡 尔曼全名Rudolf Emil ...

  5. [转载]卡尔曼滤波器及其基于opencv的实现

    卡尔曼滤波器及其基于opencv的实现 源地址:http://hi.baidu.com/superkiki1989/item/029f65013a128cd91ff0461b 这个是维基百科中的链接, ...

  6. kalman filter卡尔曼滤波器- 数学推导和原理理解-----网上讲的比较好的kalman filter和整理、将预测值和观测值融和

    = 参考/转自: 1 ---https://blog.csdn.net/u010720661/article/details/63253509 2----http://www.bzarg.com/p/ ...

  7. 卡尔曼滤波器实例:跟踪自由落体运动:设计与Matlab实现

    [首发:cnblogs    作者:byeyear    Email:byeyear@hotmail.com] 本文所用实例来自于以下书籍: Fundamentals of Kalman Filter ...

  8. [Math]理解卡尔曼滤波器 (Understanding Kalman Filter) zz

    1. 卡尔曼滤波器介绍 卡尔曼滤波器的介绍, 见 Wiki 这篇文章主要是翻译了 Understanding the Basis of the Kalman Filter Via a Simple a ...

  9. [Math]理解卡尔曼滤波器 (Understanding Kalman Filter)

    1. 卡尔曼滤波器介绍 卡尔曼滤波器的介绍, 见 Wiki 这篇文章主要是翻译了 Understanding the Basis of the Kalman Filter Via a Simple a ...

随机推荐

  1. 自己做的一个小demo

    上图: 主段代码: <script type="text/javascript"> var getRandomColor = function(){ return (f ...

  2. Azure Stack如何解决混合云的种种挑战

    微软希望能够通过Azure Stack来帮助企业连接他们的私有云和公共云.但这仍然是一项进行中的工作. 大多数企业都不愿意将所有IT运营都放到公有云中.相反,他们希望可以灵活的在这两个共享的基础架构即 ...

  3. SVN与CVS的区别大全(转载)

    本节讲解SVN与CVS的区别,主要包括是否更好的冲突标识与处理,是否有更多的本地/离线操作以及元数据管理问题. 更好的冲突标识与处理     通过是否进行更好的冲突标识与处理看SVN与CVS的区别:C ...

  4. 【BZOJ】1070: [SCOI2007]修车

    1070: [SCOI2007]修车 Description 同 一时刻有N位车主带着他们的爱车来到了汽车维修中心.维修中心共有M位技术人员,不同的技术人员对不同的车进行维修所用的时间是不同的.现在需 ...

  5. Jquery Mobile 记录

    使用的是C#语言,.Net+Jquery Mobile 框架开发 1.使用水平组切换操作 <fieldset id="Tfdset1" data-role="con ...

  6. 流动python - 一个极简主义event制

    event至少该系统的核心,以满足: 1.存储容器事件,可以被添加到事件来删除 2.触发事件fire 守则. class Event(list): def __call__(self, *args, ...

  7. Python爬虫:学爬虫前得了解的事儿

    这是关于Python的第14篇文章,主要介绍下爬虫的原理. 提到爬虫,我们就不得不说起网页,因为我们编写的爬虫实际上是针对网页进行设计的.解析网页和抓取这些数据是爬虫所做的事情. 对于大部分网页来讲, ...

  8. Emulator: glTexImage2D: got err pre :( 0x502 internal 0x1908 format 0x1908 type 0x1401

    Go to Tools > AVD Manager > Virtual device configuration > Show advanced settings > Boot ...

  9. Java笔记 #02# 带资源的try语句

    索引 普通的 try.java 带资源的 try.java 当资源为 null 的情况 可以参考的文档与资料 / test.txt 待读取的内容 hello. / 普通的 try.java 读取 te ...

  10. P2831 愤怒的小鸟(状压dp)

    P2831 愤怒的小鸟 我们先预处理出每个猪两两之间(设为$u,v$)和原点三点确定的抛物线(当两只猪横坐标相等时显然无解) 处理出$u,v$确定的抛物线一共可以经过多少点,记为$lines[u][v ...