找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

查看: 22451|回复: 54

扩展卡尔曼滤波器(EKF):一个面向初学者的交互式教程

[复制链接]
发表于 2015-5-2 10:02:13 | 显示全部楼层 |阅读模式
本帖最后由 Bluesky 于 2015-5-2 10:07 编辑

来源:home.wlu.edu/~levys/kalman_tutorial/

翻译:BlueSky  - 爱无人机论坛


在学习一些开源飞控系统像OpenPilotPixhawk时,我经常会碰到一个叫扩展卡尔曼滤波器(EKF)的东西。 谷歌上搜索这个词汇可以找到很多不同的网页以及参考文献,但我发现大部分都艰涩难懂导致像我这样的新手无法学习。 [1] 于是后来我决定编写一个从最基本的原理开始讲解的EKF教程,用于教学和学习。这个教程假定你只有高中数学水平,并会根据教学需求来介绍一些像线性代数这样更加高级的概念,而不是像其它文献一样默认你已经懂了。我们将从一些简单的例子和标准(线性)卡尔曼滤波器开始,并慢慢去理解教程后部分的实际中的扩展卡尔曼滤波器的实现。


1: 一个单的例子

想象一个飞机将要准备着陆。尽管可能有很多事情我们需要担心,如空速,燃料等,但最明显需要关注的地方是飞机此时的海拔高度。作为一个简单地近似,我们可以把当前的高度认为是之前高度的一部分。比如,如果我们观察到每次飞机的高度降低2%,则可以认为当前时间飞机的高度为上一次时间高度的98%:

  altitudecurrent_time = 0.98 * altitudeprevious_time
(altitude:高度,current_time:当前时间,previous_time:上次时间)

类似于这样使用变量上一次的值来定义当前值的方程,工程师使用“递归”这个词来表示:为了计算当前值,我们必须回归到上一次的值。最后我们会回到一个最初的”基本情况“,类似于一个已知的起始高度。

尝试移动滑块,可以看到飞机的高度变化曲线随着不同百分比而改变。


[1] 有两个值得注意的例外 the Wikipedia page, 我从这里借用了一些东西.
这个教程也从OpenPilot和DIY Drones的社区成员那里得到了不少有益的帮助 (特别是Tim Wilkins, 修正了我解释里的一些错误).

回复

使用道具 举报

 楼主| 发表于 2015-5-2 10:24:15 | 显示全部楼层
本帖最后由 Bluesky 于 2015-5-2 11:00 编辑

3: 放到一起

此时我们有了两个描述飞机状态的方程:

8.jpg

这些方程非常容易理解,但他们不具有通用性,仍不足以描述那些除了我们这个“飞机-高度”样例之外的系统。为了使方程变得更具一般性,工程师们采用数学上约定俗成的符号像 9.jpg
,   10.jpg
, 来表示变量, 11.jpg
12.jpg
表示常量,下标 13.jpg
表示时间[3]。然后我们的方程变成如下形式:
6.jpg
代表系统的当前状态, 14.jpg
表示上一次的状态,是一个常量(在我们的例子里为0.98), 15.jpg
是我们对系统当前的观察值, 16.jpg
是当前的测量噪声。

卡尔曼滤波器这么流行的原因之一便是它可以让我们在给出观察值常量,以及测量噪声的整体数量 19.jpg
之后,便能得到对系统当前状态 17.jpg
的一个非常好的估计值

另外,我们应当考虑到,飞机的真实高度变化并非是一条平滑的曲线。任何一个曾坐过飞机的人都会告诉你,飞机在下降着陆时通常会遇到一定震荡。这个震荡被定义为高度上的噪声,可以视为另一种噪声信号:

18.jpg
altitude:高度,turbulence:震荡

一般而言:
7.jpg

20.jpg
被称为过程噪声,就像震荡,它是飞机降落过程中固有的一部分,而不是人工的观察和测量。为了将注意力放在其它主题上,我们将忽略这个过程噪声一小会,直到传感器融合那一章。[4]


[3] 为什么不用t 来表示时间?可能是因为这里时间被视为一个序列的离散步骤,这种情况下通常用k 表示。
[4] 我还忽略了状态方程中的控制信号,因为它对于大多数卡尔曼滤波器方程都是次要的,且当需要的时候可以很容易加上 。

回复 支持 2 反对 0

使用道具 举报

 楼主| 发表于 2015-5-2 10:06:49 | 显示全部楼层
本帖最后由 Bluesky 于 2015-5-2 10:15 编辑

2: 噪声的处

当然,现实世界中会用一个传感器像GPS或气压计来测量高度值。这种传感器会输出具有不同精确度的值。[2] 如果传感器输出的是一个恒定的量,那我们可以对变量进行简单地加减来确定我们的高度。然而通常情况下,传感器的精度时刻在发生着不可预知的变化,使得传感器读出的是真实高度值的一个具有“噪声”的版本:

  observed_altitudecurrent_time = altitudecurrent_time + noisecurrent_time
(observed_altitude:观察的高度,altitude:真实高度,noise:噪声,current_time:当前时间)

尝试移动上方的滑块,可以看到噪声对观察的高度值的影响。噪声可表示为可观察的高度范围的百分比。


3.jpg

0% noise


4.jpg

6% noise

5.jpg

15% noise



[2]比如,Garmin公司公布其气压高度计读数精度为“适当校准下为10英尺”。打个比方,当高度计读到的值为1000英尺时,当前的真实高度可能是从990到1010英尺之间的任何一个数。


回复 支持 反对

使用道具 举报

 楼主| 发表于 2015-5-2 11:04:49 | 显示全部楼层
本帖最后由 Bluesky 于 2015-5-2 11:20 编辑

4: 状态估计

在这里(忽略过程噪声)是两个描述我们观察到的系统状态的方程:
21.jpg

既然我们的目的是想从观察值 35.jpg
得到状态 9.jpg
,我们可以将第二个方程重写为:
22.jpg

当然你很容易发现这里有个问题是我们不知道当前的测量噪声 16.jpg
,它是不可预测的。幸运的是,卡尔曼具有洞察力,我们可以依靠当前的观察值和上一次的状态估计值来估计当前的状态。工程师们通常在一个变量上面加一个“帽子”符号^ 来表示这个变量是一个估计值。因此 29.jpg
表示当前状态的估计。然后我们可以通过权衡上一次的状态和当前的观察值来表示当前的估计:
23.jpg

其中 31.jpg
代表“增益”,即权重值[5]。这个方程被红色高亮,因为它会在我们实现卡尔曼滤波器时被直接用到。
现在,这一切似乎看起来相当复杂了,但思考一下,当增益系数 32.jpg
取两个极限值时会发生什么事情。
33.jpg
,我们得到
24.jpg

换句话说,当增益为0时,观察值是无效的,我们得到这样一个方程:当前的状态等于上一次的状态
34.jpg
,我们得到
25.jpg

即当增益为1时,上一次的状态便无关紧要了,我们完全依靠当前的观察值来获得当前的状态估计。

当然,实际的增益值可能介于两个极限值之间。尝试移动下面的滑块来改变增益系数,可以看到它对当前状态估计的影响。
26.jpg

27.jpg
28.jpg



[5] 通常用变量k 来表示增益,因为它以卡尔曼增益Kalman gain)闻名。怀着对Rudolf Kalman的尊敬,我发现使用相同的字母k来表示变量和下标容易被混淆,所以我选择了字母g 代替。
回复 支持 反对

使用道具 举报

 楼主| 发表于 2015-5-2 11:23:02 | 显示全部楼层
本帖最后由 Bluesky 于 2015-5-5 23:54 编辑

5: 计算增

所以我们现在有了一个方程,可以实际用于计算当前状态的估计 29.jpg
,基于上一次的状态 42.jpg
,当前的观察值 15.jpg
和当前的增益 32.jpg
36.jpg

然而,我们应该如何计算增益系数?答案是:间接地从噪声计算。要记得,每次的观测值都是与特定的噪声相关联的:
37.jpg

我们不知道一次观测中的个体噪声值,但通常我们可以知道噪声的平均值:比如,一个公布了精度值的传感器通常会告诉我们它的输出噪声平均近似是多少,我们用 r来表示这个值。它没有下标,因为它跟时间无关,而是代表着一个传感器的性能。然后我们可以用 r来计算当前的增益 51.jpg
38.jpg

40.jpg
是一个通过递归计算得到的预测误差: [6]
39.jpg

在继续之前,先让我们与状态估计公式一起思考一下上面两个公式的含义

假如说误差 41.jpg
在我们的上一次预测中为0,则我们当前的增益 50.jpg
则为0/(0+r)=0那我们的下一次状态估计将等同于当前的状态估计。这是有道理的,因为如果我们的预测是准确的预测误差为0),便不需要调整我们当前的状态估计值

另一种极端情况,如预测误差为1,则此时增益为1/(1+r)如果 r 为0 — 即,我们的系统里的噪声小到可以忽略不计 —  则增益为1,然后新的状态估计 17.jpg
将受到观测值的强烈影响。

但随着 r 增大,增益可以变得任意小。换言之,当系统噪声过大时,坏的预测将被忽略。

第三个公式呢?递归上一次的值 47.jpg
和当前增益 48‘.jpg
来计算预测误差 49.jpg
?同样的,设增益为极值时有助于我们思考发生了什么:当 43.jpg
, 有 44.jpg
。所以,类似于状态估计,增益为0意味着预测误差没有变化。而当 45.jpg
时,有 46.jpg
。因此,最大增益值对应着零预测误差,即仅使用当前的观测值来更新当前的状态。

[6]技术性地来说r实际上噪声信号的方差;即,个体噪声值与其平均值之差的平方的和的平均数。如果这个噪声方差随着时间而改变,卡尔曼滤波器将会运行得非常好,但在大多数应用中它可以被假定为常量。同样的道理, 13.jpg
步骤k中估计过程的协方差;它是我们的预测的均方误差的平均值。事实上,正如Tim Wikin向我指出的,状态是随机的变量/矢量(一个随机过程的瞬时值),它从来没有一个“对的”值可言!估计仅仅是描述状态的过程模型的最具可能性的值。


回复 支持 反对

使用道具 举报

 楼主| 发表于 2015-5-2 11:44:17 | 显示全部楼层
本帖最后由 Bluesky 于 2015-5-24 13:38 编辑

6.预测和更新

我们差不多准备好运行卡尔曼滤波器并看到一些结果了。首先,可能你想知道最初的状态方程中的常量a 发生了什么:
5.JPG
它似乎在我们的状态估计方程中已经消失了:
6.JPG
答案是,基于不同的信息这两个方程我们都要分别用于状态估计。最初的方程代表着一个表示“状态应该成为什么样”的预测,而第二个方程代表着一个对这个预测的更新,基于一个观测。 [7] 所以我们重写一下最初的方程,用x 上添加一个帽子符号,表示这是一个估计:
7.JPG
最后,我们还要在误差的预测中使用常量 a  :[8]
8.JPG
这两个红色标注的方程代表着卡尔曼滤波器的预测阶段。这个循环——预测/更新,预测/更新,...将会重复我们所需要的次数。

[7] 从技术上来讲,第一个估计称为先验估计,第二个为后验估计,且大部分文献会引入一些额外的上标或下标加以区别。因为我想让事情变得简单(可以让你用你最喜欢的编程语言来写代码!),所以我避免让符号变得更复杂。

[8] Zichao Zhang 友好地向我指出,我们乘以a两次,是因为预测误差 9.JPG
本身是一个均方误差;因此,它是被状态值
10.JPG
的相关系数的平方缩放。用
11.JPG
来代表预测误差而不是 12.JPG
的原因将会在第十二章说明。
回复 支持 反对

使用道具 举报

 楼主| 发表于 2015-5-2 11:44:26 | 显示全部楼层
占位编辑
回复

使用道具 举报

 楼主| 发表于 2015-5-2 11:44:33 | 显示全部楼层
占位编辑
回复

使用道具 举报

 楼主| 发表于 2015-5-2 11:44:41 | 显示全部楼层
占位编辑
回复

使用道具 举报

 楼主| 发表于 2015-5-2 11:44:46 | 显示全部楼层
占位编辑
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表