卡尔曼滤波最完整公式推导

技术讨论 chengzi ⋅ 于 1个月前 ⋅ 150 阅读

作者丨东林钟声@知乎(已授权)
来源丨https://zhuanlan.zhihu.com/p/341440139​
编辑丨极市平台
卡尔曼滤波是一种利用线性系统状态方程,通过系统输入输出观测数据,对系统状态进行最优估计的算法。由于观测数据中包括系统中的噪声和干扰的影响,所以最优估计也可看作是滤波过程。上面一段话来自百度百科,其实最核心的意思就是卡尔曼滤波可以很好地从带有噪声的数据过程中估计状态。而且卡尔曼滤波也是阿波罗登月中使用的突破性技术之一,正好我国嫦娥五号采土归来,正好回顾一下卡尔曼滤波的技术细节,对公式推导做一个完整的整理。

1、定义

  • 状态方程:

$$
x(k+1) = F(k)x(k) + G(k)u(k) + v(k) \ u(k)为输入信号、x(k)为状态变量、v(k)表示过程噪声\ 其中 \text{cov}{v(k)} = Q(k)表示协方差矩阵
$$

  • 观测方程:

$$
z(k) = H(k)x(k) + \omega(k) \ z(k)为观测变量、\omega(k)表示观测噪声 \ 其中 \text{cov}{\omega(k)} = R(k)表示协方差矩阵
$$

  • 符号定义:

首先定义状态估计误差

$$
\tilde{x}(k+1|k) = x(k+1) - \hat{x}(k+1|k)\ \hat{x}(k+1|k)表示估计在第k时刻,对x(k+1)的估计
$$

理解卡尔曼公式推导的核心,是一定要理解这里的符号标记。 核心是在第k时刻,对k+1时刻的估计!所以后面都要记住,在第k时刻的时候,对k+1时刻的任何值,都只能是估计(预测未来值)。具体怎么估计的,如下:

$$
状态估计方程\quad\hat{x}(k+1|k) = F(k)\hat{x}(k|k) + G(k)u(k)\
$$

同理,我们也可以得到观测估计误差

$$
\tilde{e}(k+1)= z(k+1) - \hat{z}(k+1|k) \ 观测估计方程 \quad \hat{z}(k+1|k) = H(k+1)\hat{x}(k+1|k)
$$

基于上面的公式,我们定义两个重要的误差协方差矩阵

$$
状态误差协方差矩阵\quad P(k+1|k) = \text{cov}{\tilde{x}(k+1|k)} \ 观测误差协方差矩阵\quad S(k+1) = \text{cov}{\tilde{e}(k+1)}
$$

  • 最终目的

其实卡尔曼的最终目的,是得到一个基于误差能够不断修正迭代式估计表达式,其具体形式应该如下:

$$
\hat{x}(k+1|k+1) = \hat{x}(k+1|k) + W(k+1)\tilde{e}(k+1) \
$$

这个式子非常直观,就是基于误差去修正,怎么最优的去修正?就是我们怎么来算这个W,也就是卡尔曼增益(Kalman Gain)。

2、推导

为了得到不停的基于误差修正的一种计算模式,可能要能够得到一些关键的递推形式,如我之前发表过关于递推最小二乘的推导。首先我们推导关于状态误差协方差矩阵:

$$
\begin{aligned}P(k+1|k) &= \text{cov}{x(k+1) - \hat{x}(k+1|k)} \ &=\text{cov}{F(k)x(k) + G(k)u(k) + v(k) - F(k)\hat{x}(k|k) - G(k)u(k)} \ &=\text{cov}{F(k)(x(k) - \hat{x}(k|k)) + v(k)} \ &=F(k)\text{cov}{x(k) - \hat{x}(k|k)}F^T(k) + \text{cov}{v(k)} \ &= F(k)P(k|k)F^T(k) + Q(k) \end{aligned}
$$

这里我们得到了P(k|k)P(k+1|k) 的一个递推形式,注意这里P(k|k) 的标记形式是为了更方便后面的推导,为什么这么说,其实是为了得到

$$
P(k|k) \rightarrow P(k+1|k) \rightarrow P(k+1|k+1)\
$$

这样一种能够不停向前迭代的形式。带着这个思考继续下面的推导,就会明白卡尔曼的巧妙之处。

我们先总结下,这里得到了P的一个递推形式:

$$
P(k+1|k) = F(k)P(k|k)F^T(k) + Q(k) \
$$

同理,也是可以得到观测误差协方差的一个递推形式:

$$
S(k+1) = H(k+1)P(k+1|k)H^T(k+1) + R(k+1) \
$$

会发现其实S(k+1)S(k) 没啥直接的关系,这就是为什么不把S(k+1) 写成像P那样 S(K+1|K) 的形式,因为没有必要,只有P是需要这样来写的。

这里再次思考,我们的目的是求一个最优的W,就是卡尔曼增益,具体怎么来最优?这里就用到了状态估计误差,就是说,我们的目的是让状态误差的平方和最小,这里可以使用P的迹。具体推导如下:

$$
\begin{aligned}P(k+1|k+1) &= \text{cov}{x(k+1) - \hat{x}(k+1|k+1)} \ &= \text{cov}{x(k+1) - (\hat{x}(k+1|k) + W(k+1)\tilde{e}(k+1)} \ &= \text{cov}{x(k+1) - (\hat{x}(k+1|k) + W(k+1)(z(k+1) - \hat{z}(k+1|k))}\ &=\text{cov}{x(k+1) - (\hat{x}(k+1|k) + W(k+1)[H(k+1)x(k+1) + \omega(k+1)- H(k+1)\hat{x}(k+1|k)]} \ &=\text{cov}\left[(I-W(k+1)H(k+1))(x(k+1)-\hat{x}(k+1|k)) - W(k+1)\omega(k)\right] \ &= (I-W{k+1}H{k+1})P(k+1|k)(I-W{k+1}H{k+1})^T + W{k+1}R(k+1)W^T{k+1} \end{aligned}
$$

然后使用P的迹对W求导,令其为0得到最优的W:

$$
\frac{\partial\ {trace(P(k+1|k+1))}}{W{k+1}} = -2(H{k+1}P(k+1|k))^T + 2W{k+1}S{k+1} = 0
$$

得到:

$$
W{k+1} = P(k+1|k)H^T{k+1}S^{-1}_{k+1}\
$$

然后再把W带回P(k+1|k+1)的式子中:

$$
P(k+1|k+1) = (I - W{k+1}H_{k+1})P(k+1|k)\
$$

到这里,推导完毕,得到了W的最优表达式,同时,也解决了我们上面提出的问题,一个完整的递推链:

$$
P(k|k) \rightarrow P(k+1|k) \rightarrow P(k+1|k+1)\
$$

建立完毕。那么通过这一套完整的递推链,给定P的一个初始估计P(0|0) 就可以按照下面的链来进行状态估计:

$$
P(0|0) \rightarrow P(1|0) \rightarrow P(1|1) \rightarrow P(2|1) \rightarrow \cdots\
$$

到这里卡尔曼滤波的公式推导完成了,这个版本是我在看过很多其他资料,反复提炼之后,标记最友好,推导最友好的一个版本,理解卡尔曼滤波一定要首先理解符号的定义,特别是下标,以及怎么样构建递推链。在掌握之后可以自行通过上面的推导总结出,其他教程反复提到的黄金五条公式,其实这些在我看来不是重点,重点是理解其本质原理。

3、资料

其实要弄懂卡尔曼滤波还真不简单,需要一些基本矩阵、统计、自控原理以及矩阵求导的知识,特别是对迹求导,这里推荐Matrix CookBook,里面详细列举了各种常用的求导方式。

小迷离

成为第一个点赞的人吧 :bowtie:
回复数量: 0
暂无回复~
您需要登陆以后才能留下评论!