时间序列八: 以NASA之名: 卡尔曼滤波器

以NASA之名: 卡尔曼滤波器

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

目录

• 以NASA之名: 卡尔曼滤波器
  • 引言
  • 荣耀骑士
  • 卡尔曼滤波器*
  • 参考文献:

引言

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

荣耀骑士

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

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

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

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

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

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

卡尔曼滤波器*

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

模型:

[left{ egin{array} x_{k} = A_kx_{k-1}+B_k u_k + w_k\ y_k = C_k x_k + v_k end{array} ight. ]

其中所有变量均为矩阵(或向量), (x_k): k 时刻的状态向量; (A_k): 状态转移矩阵; (B_k): 控制矩阵,(u_k): 控制向量; 一般,将此控制部分当作确定部分而不包含在模型中,因此上式可简化为:

[left{ egin{array} x_{k} = A_kx_{k-1} + w_k\ y_k = C_k x_k + v_k end{array} ight. ]

(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):

[ ilde{x}_k = x_k - hat{x}_k\ min E[ ilde{x}_k ilde{x}_k^T] ightarrow 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) 现在称为状态校正后的估计值.

观测向量的误差(也即新息):

[ ilde{y}_k = y_k - hat{y}_k ]

用观测值来校正状态变量:

[egin{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)( ilde{y}_k)即是在 k 时刻之前没有但在 t 时刻产生的新的新信息. 用代数的语言即是 新息与过去的信息是正交的:

[E[ ilde{y}_ky_j^T] = oldsymbol{O}, qquad 1le j<k ]

而且可以容易的推论出, 新息之间也是正交的:

[E[ ilde{y}_i ilde{y}_j^T] =oldsymbol{O},qquad i e j, andquad 1le i,jle k ]

从而可知, 新息( ilde{y}_k)是一个与k 之前时刻的数据不相关,且有白噪声性质的随机过程.

校正后的状态变量估计误差及其均方值 ((P_k)) 还有未经校正的状态变量估计误差的均方值 ((P_k')):

[egin{array}\ ilde{x}_k &=& x_k - hat{x}_k\ P_k &=& E[ ilde{x}_k ilde{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 ightarrow H_k ightarrow hat{x}_k ]

[egin{array}\ ilde{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{egin{array}\ hat{x}_k &=& A_k hat{x}_{k-1} + H_k(y_k - C_kA_khat{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} ight. ]

这就是卡尔曼滤波器的递推公式了.

参考文献:

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.