一、卡尔曼滤波器要解决的问题

  首先说一下卡尔曼滤波器要解决的是哪一类问题,这类系统应该如何建模。这里说的是线性卡尔曼滤波器,顾名思意,那就是线性动态的离散系统。这类系统可以用如下两个方程来表示:

\[\begin{array}{l}
 x(n + 1) = {\bf{F}}(n + 1,n)x(n) + {v_1}(n) \\
 y(n) = {\bf{C}}(n)x(n) + {v_2}(n) \\
 \end{array}\]

  其中:  

  x(n)表示系统的状态

  F(n+1,n)为状态转移矩阵,表示状态随时间的变化规律。通俗的讲,从当前状态到下一个状态之间有什么关系。

  C(n)表示观测值与状态的关系

  y(n)表示状态的观测值

  v1表示系统过程的噪声

  v2表示观测过程中产生的噪声

  上面的两个方程中,第一个方程是过程方程,它表示系统状态x(n)随时间的更新过程。第二个方程为测量方程,表示状态x(n)与测量结果y(n)的关系。这里我们要先对这两个方程中的概念做下解释。

状态是对系统特征进行的一个抽象,由预测系统未来特性时所需要的、与系统过去行为有关的最少数据组成。

  这个概念不好理解吧!那么举个例子。相信不少朋友在网上看到过有人拿房间温度测量来讲述卡尔曼滤波原理。这里房间里真实的温度就是状态,它可以是一个参数,也可以是多个参数。那么,用温度计测出来的值,就是这里的观测值y(n)。再说一个例子,假如我们要对一个运动的物体进行跟踪,那么,物体的位移和速度完全可以表示这个运动物体所组成的系统的主要特征。这时的状态就可以用一个具有位移和速度两个特征的向量来表示。解释到这里,相信很多朋友已经正确理解了状态这个概念,它表示的是系统客观存在的真实特征

  再说一下系统状态与其观测值之间为什么有C(n)的存在,这里它表示的是观测值与状态的关系。再拿室内测度测量来举例子,室内客观真实温度(未知量)做为这个系统中的状态,用温度计来测量这个状态。测出来的温度就是我们的观测温度y(n)。这里很可能系统状态跟其观测值不是简单的加个测量噪声的关系。那么这种关系就可以用C(n)来建模。

  上面的两个方程,就是线性卡尔曼滤波器对要解决的系统的建模。这里再次对这两个方程表示的模型做出更加通俗的解释:

  过程方程:它讲述的是系统从一个状态到另一个状态应该随时间如何变化

  测量方程:它讲述的是当前状态与当前测量值之间的关系

  也就是说,我们如果对一个系统感兴趣,首先要找出这个系统的一个或者几个主要特征(状态),然后对这几个特征进行观测,并得到一组观测值,通过不断的观测来认识这个系统,利用仅有的观测数据找到最优的方式来求解过程方程和测量方程。这就是卡尔曼滤波器要解决的问题。为了便于理解,后面我们只讨论只有一个特征的状态所组成的观测系统。

 

二、投影定理与新息过程

  现在我们已经清楚了卡尔曼滤波器要解决什么样的问题,对于这样的问题应该如何建立模型。下面我们来看下具体用什么样的方法来最优的求解过程方程和测量方程。卡尔曼滤波器的推导过程有两种思路:

  (1)正交投影定理(有的地方也称为射影定理,这里后面简称为投影定理)

  (2)新息过程

  卡尔曼老先生最初介绍卡尔曼滤波器时是使用正交投影定理来推导的。后来Kailath又提出来基于新息过程的推导方法。这里打算两个都介绍一下。上面说到,当对系统观测时,会得到一组观测值y(n),那好,我们首先考虑一个估计的问题,更具体的说是如何利用仅有的观测数据来预测状态x(n)的问题。由回归分析的内容我们可以知道,由观测值y求随机变量x的线性最小二乘估计(线性最小方差估计)可以表示为:

\[\hat x = E[x] + Cov(x,y)D{(y)^{ - 1}}\{ y - E[y]\} \]

  这时,随机变量x与其估计的差值跟观测值y正交(不相关),我们就称x的估计值为x在y的上投影,记为:

\[\hat x = proj[x|y] = proj[x|y(1),y(2),...,y(k)]\]

  其几何意义表示如下:

Kalman Filter卡尔曼滤波 java实现 详解卡尔曼滤波_卡尔曼滤波

新息,它表示当前观测值与其估计值的误差。可以表示为下面的方式:

\[\alpha (n) = y(n) - proj[y(n)|y(1),y(2),...,y(n - 1)] = y(n) - \hat y(n|{Y_{n - 1}})\]

  新息几何意义表示如下:

Kalman Filter卡尔曼滤波 java实现 详解卡尔曼滤波_卡尔曼滤波_02

  再来看新息产生的过程,我们得到一个观测值,就能同时得到一个对其对应的误差。反之也是成立的,也就是说观测值与新息,存在着一一对应的关系,那如果我们想要得到最终的系统状态,也新息来代替观测值来估计也是可以的,即:

\[\{ y(1),y(2),...,y(n)\}  \mathbin{\lower.3ex\hbox{$\buildrel\textstyle\rightarrow\over
{\smash{\leftarrow}\vphantom{_{\vbox to.5ex{\vss}}}}$}} \{ \alpha (1),\alpha (2),...,\alpha (n)\} \]

  因此

\[proj[x|y(1),y(2),...,y(n)] = proj[x|\alpha (1),\alpha (2),...,\alpha (n)]\]

  写成递推的形式为:

\[\begin{array}{l}
 proj[x|y(1),y(2),...,y(n)] = proj[x|\alpha (1),\alpha (2),...,\alpha (n)] \\
  = proj[x|\alpha ] = E[x] + \sum\limits_{i = 1}^{n - 1} {E\{ x} {\alpha ^T}(i)\} E{\{ \alpha (i){\alpha ^T}(i)\} ^{ - 1}}\alpha (i) + E\{ x{\alpha ^T}(n)\} E{\{ \alpha (n){\alpha ^T}(n)\} ^{ - 1}}\alpha (n) \\
  = proj[x|\alpha (1),\alpha (2),...,\alpha (n - 1)] + E\{ x{\alpha ^T}(n)\} E{\{ \alpha (n){\alpha ^T}(n)\} ^{ - 1}}\alpha (n) \\
  = proj[x|y(1),y(2),...,y(n - 1)] + E\{ x{\alpha ^T}(n)\} E{\{ \alpha (n){\alpha ^T}(n)\} ^{ - 1}}\alpha (n) \\
 \end{array}\]

   看来,利用投影的思想,我们是可以递归的表达出来状态的预测过程的。也可以说,在几何上,kalman滤波器可以看做状态在由观测生成的线性空间上的投影。但如果只分析到这里,估计很多朋友可能还不能真正的理解卡尔曼滤波器的思想。那么我们再以新息过程来分析下如何来求解过系统状态。

 

  认真考虑下就会发现,整个过程有两个误差。一个是观测过程中的预测误差:新息;一个是状态更新过程中的预测误差(一直隐含在投影过程中没有被注意到),我们来看下它们之间的关系。

  由过去所有的观测值y(1),...,y(n-1)和系统的测量方程,可以得到当前观测值y(n)的最小均方估计为

\[\hat y(n|{Y_{n - 1}}) = {\bf{C}}(n)\hat x(n|{{\bf{Y}}_{n - 1}})\]

  于是,新息过程还可以表示为:

\[\begin{array}{*{20}{c}}
   {\alpha (n) = y(n) - \hat y(n|{{\bf{Y}}_{n - 1}})}  \\
   { = y(n) - {\bf{C}}(n)\hat x(n|{{\bf{Y}}_{n - 1}})}  \\
   { = {\bf{C}}(n)\varepsilon (n,n - 1) + {v_2}(n)}  \\
\end{array}\]

  这里ε(n,n-1)是状态x(n)和其一步预测值之差

  \[\varepsilon (n,n - 1) = x(n) - \hat x(n|{{\bf{Y}}_{n - 1}})\]

  如果我们定义

\[\begin{array}{l}
 {\bf{R}}(n) = E[\alpha (n){\alpha ^H}(n)] \\
 {\bf{K}}(n,n - 1) = E[\varepsilon (n,n - 1){\varepsilon ^H}(n,n - 1)] \\
 {{\bf{Q}}_2}(n) = {v_2}(n)v_2^H(n) \\
 \end{array}\]

  这里R(n)为新息误差相关矩阵,代表观测信息的质量优劣、K(n,n-1)为状态误差相关矩阵,代表状态时间更新的质量优劣、Q2(n)表示测量噪声向量的相关矩阵。再由观测方向可得:

\[{\bf{R}}(n) = {\bf{C}}(n){\bf{K}}(n,n - 1){{\bf{C}}^H}(n) + {{\bf{Q}}_2}(n)\]

  这就是状态预测误差与新息之间的二阶统计关系。因此,只要我们尽可能的优化时间更新的质量,同时提高观测信息的质量,那么,我们一定能估计出想要的系统状态。好,问题反过来,再回到文初提出来的问题,我们有了一组观测值,又得到了与比对应的新息序列。如何好好利用它们估计现我们想要的系统状态。下面就要好好说下:

  根据观测值与新息之间的一一对应关系,如果我们以新息过程来估计状的最小均方估计,则这个估计可以表示为新息过程序列的线性组合。即:

\[\hat x(i|{{\bf{Y}}_n}) = \sum\limits_{k = 1}^n {{{\bf{B}}_i}(k)\alpha (k)} \]

  这里B为待定系数矩阵,那么,根据线性滤波的正交原理,预测状态误差向量与新息过程正交,即:

\[E[\varepsilon (i,n){\alpha ^H}(m)] = E\{ [x(i) - x(i|{{\bf{Y}}_n})]{\alpha ^H}(m)\}  = 0\]

  再将新息过程导出来状态x(i)的最小均方估计代入,可得

\[E[x(i){\alpha ^H}(m)] = {{\bf{B}}_i}(m)E[\alpha (m){\alpha ^H}(m)] = {{\bf{B}}_i}(m){\bf{R}}(m)\]

  这里的R(m)为新息过程的相关矩阵。上式两边同时乘以R(m)的逆矩阵,可得待定系数矩阵B的表达式为

\[{{\bf{B}}_i}(m) = E[x(i){\alpha ^H}(m)]{{\bf{R}}^{ - 1}}(m)\]

  再将B代入新息过程导出来状态x(i)的最小均方估计公式中,可以得到状态x(i)的最小均方误差估计为:

\[\begin{array}{l}
 \hat x(i|{{\bf{Y}}_n}) = \sum\limits_{k = 1}^n {E[x(i){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k)}  \\
  = \sum\limits_{k = 1}^{n - 1} {E[x(i){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k)}  + E[x(i){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k) \\
 \end{array}\]

  对于i=n+1,有

\[\hat x(n + 1|{{\bf{Y}}_n}) = \sum\limits_{k = 1}^{n - 1} {E[x(n + 1){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k)}  + E[x(n + 1){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k)\]

   我们看到,这个结果与投影定理推导的结果的含义是相同的,只是使用的工具和准则不同。

 

  最后,我们再回想下观测数据与新息、状态之间的关系,以及根据这些关系导出的结果。不难得出结论,卡尔曼滤波器的出发点是尽可能的好好利用观测值,用尽可能少的计算量来估计系统状态。呵呵,是不是感觉数字的游戏还是很好玩的,并不一定很枯燥。好的,今天就分析到这里吧,后续将会使用这次分析的结果,逐渐对卡尔曼滤波器及其使用做一个说明。这里的分析内容来源于我在学习卡尔曼滤波器过程中的一些思考(参阅书籍:自适应滤波器原理、卡尔曼滤波原理与应用),水平有限,如有错误之处,还请各位大牛不吝赐教。