零基础读懂“扩展卡尔曼滤波”——中篇

https://mp.weixin.qq.com/s/cr1u840FeBLalrmtJQs8TQ

本篇文章分上、中、下三篇,上篇从标准卡尔曼滤波开始,中篇加入更真实的系统模型,下篇从传感器的数据融合中实现扩展卡尔曼滤波。

8. 一个更真实的模型

现在重新回顾一下描述系统的状态方程和观测方程:
x_k=ax_{k-1}
z_k=x_k+v_k

x_k 表示系统的当前状态,x_{k-1}表示系统的前一个状态;a是一个常量,z_k是系统当前状态的观测值;v_k是系统观测噪声;

这两个方程已经适用于大多数系统,但仍然不够普适性;现在依然以飞机的飞行为例,我们并没有考虑到飞行员对飞机的操作和控制,飞行员操作控制杆向前或向后移动,对飞机输入控制量,最终对飞机产生控制。我们将这个控制量称为u_k,表示飞行员发送到飞机的控制信号的当前值,并且我们对这个控制量加上一个缩放因子b,因此整个状态方程变为:

x_k=ax_{k-1}+bu_k

一般而言,除了噪声以外的任何信号都可以通过常量进行一定比例的缩放,因此观测方程也可以这样表示:

z_k=cx_k+v_k

9. 调整估计值

通过以上的方程我们更真实的描述一个系统的状态和观测,由此预测和更新方程也需要做相应的调整:

预测:

\hat x_k=a\hat x_{k-1}+bu_k

更新:

g_k=p_kc/(cp_kc+r)
\hat x_k \gets\hat x_k+g_k(z_k-c\hat x_k)
p_k\gets(1-g_kc)p_k

仍然以飞机飞行为例,我们增加一个控制信号,表示飞行员稳步拉回控制杆,并提高飞机飞行高度,原始信号为蓝色表示,观测值以红色表示,绿色信号为卡尔曼滤波后的值。

image

10. 向系统中添加速度

现在回想一下飞机下降时飞行高度的原始方程:

altitude_{current}=0.98\times altitude_{previous}

以更通用的方程表示为:

x_k=ax_{k-1}

我们都知道飞机飞行的高度即海拔,可以看做是一种距离,距离总是和速度与时间相关的:

distance=velocity\times time

再进一步思考,当前距离和前一时刻的距离的关系:

distance_{curr}=distance_{prev}+velocity_{prev}\times(time_{curr}-time_{prev})

下标curr代表当前,prev代表前一个。也就是说,飞机现在的位置等于前一时刻的位置加上刚刚飞行过的垂直距离。如果按照固定的时间间隔进行计算,比如说1秒,或者1分钟等,那么可以将上式简化为:

distance_{curr}=distance_{prev}+velocity_{prev}\times timestep

貌似与我们的状态方程更接近了一步,但是还差一些,接下来我们引入向量和矩阵。

11. 向量和矩阵

通常情况下系统的状态并不是一个变量,比如飞机的高度和速度,那么我们就可以使用向量表示系统的状态:

x_k \equiv \left[ \begin{matrix} distance_k \\ velocity_k \end{matrix} \right]

符号\equiv表示一个定义,系统状态定义为一个距离和速度的矢量。 当我们用一个矩阵去乘以一个向量时得到的结果依然是一个向量,比如:

\left[ \begin{matrix} a & b \\ c & d \end{matrix} \right] \left[ \begin{matrix} x \\ y \end{matrix} \right] = \left[ \begin{matrix} ax+by \\ cx+dy \end{matrix} \right]

因此,我们将系统状态方程中的常量A定义为一个矩阵,因为系统的状态并非是一种:

A= \left[ \begin{matrix} 1 & timestep \\ 0 & 1 \end{matrix} \right]

再次回顾一下系统状态方程:

x_k=Ax_{k-1}

将系统状态矢量代入:

\left[ \begin{matrix} distance_k\\ velocity_k \end{matrix} \right] = \left[ \begin{matrix} 1 & timestep \\ 0 & 1 \end{matrix} \right] \left[ \begin{matrix} distance_{k-1}\\ velocity_{k-1} \end{matrix} \right] = \left[ \begin{matrix} 1\times distance_{k-1}+timestep\times velocity_{k-1}\\ 0\times distance_{k-1}+1\times velocity_{k-1} \end{matrix} \right] = \left[ \begin{matrix} distance_{k-1}+timestep\times velocity_{k-1}\\ velocity_{k-1} \end{matrix} \right]

从公式中可以清晰的看出距离等于前一刻的距离与速度和观测时间间隔的乘积之和,而速度和上一时刻的速度相等;这表明飞机是匀速上升或下降的,但现实中速度不可能做到匀速。当速度是变化的,即系统中还存在一个加速度变量,那么我们的系统状态方程可以表示为:

\left[ \begin{matrix} distance_k\\ velocity_k\\ acceleration_k \end{matrix} \right] = \left[ \begin{matrix} 1 & timestep & 0 \\ 0 & 1 & timestep \\ 0 & 0 & 1 \end{matrix} \right] \left[ \begin{matrix} distance_{k-1}\\ velocity_{k-1}\\ acceleration_{k-1} \end{matrix} \right]

12. 再看预测和更新

系统状态的修正公式:

x_k=Ax_{k-1}

x为一个向量,A是一个矩阵;大家可能还记得,这个方程的原始形式是这样的:

x_k=ax_{k-1}+bu_k

其中u_k是一个控制信号,b是一个比例因子,系统的观测方程为:

z_k=cx_k+v_k

其中z_k是观测量,v_k是由于传感器的不准确所带来的噪声。那么我们如何修正这些方程以适应新的向量矩阵形式呢?是的,线性代数让这些变得尤其简单。我们使用大写的字母写出缩放因子bc,使他们成为矩阵,而不是单一的标量值:

x_k=Ax_{k-1}+Bu_k
z_k=Cx_k+v_k

然后所有的变量包括状态、观测、噪声、控制等都被认为是向量,回顾一下我们前面的系统预测和更新方程:

预测:

\hat x_k=a\hat x_{k-1}+bu_k
p_k=ap_{k-1}a

更新:

g_k=p_kc/(cp_kc+r)
\hat x_k \gets \hat x_k+g_k(z_k-c\hat x_k)
p_k\gets(1-g_kc)p_k

大家可能会想,将这里所有的缩放因子都用大写表示,然后就可以了。然而事实并非如此,矩阵的乘法没有那么简单,这就是前面提到的为什么是ap_{k-1}a而不是a^2p_{k-1}A\times B有时候不等于B\times A,有时候需要转置矩阵,通过在矩阵的旁边加上一个上标^T来表示转置矩阵。转置矩阵是通过把每一行变成一列,每一列变成一行来实现的。看个实际例子:

\left[ \begin{matrix} a & b \\ c & d \end{matrix} \right]^T = \left[ \begin{matrix} a & c \\ b & d \end{matrix} \right]

由此,我们可以重写我们的预测方程: (关于卡尔曼滤波五个公式的由来,我会专门写一篇文章进行推导与分析)

\hat x_k=A\hat x_{k-1}+Bu_k
P_k=AP_{k-1}A^T

我们已经知道A为何是一个矩阵,但是P_k为何也是一个矩阵呢?

我们的第二个更新方程可以表示为:

\hat x_k \gets \hat x_k+g_k(z_k-C\hat x_k)

再看增益方程的原始形式:

g_k=p_kc/(cp_kc+r)

它包含一个除法,我们瞬间想到了除法可以用乘法代替:

a\times a^{-1}=1

那么上面的方程可以写为:

g_k=p_kc(cp_kc+r)^{-1}

我们现在把这些常量替换成矩阵:

G_k=P_kC^T(CP_kC^T+R)^{-1}
P_k=(I-G_kC)P_k

如何计算(CP_kC^T+R)^{-1}呢,求逆矩阵吗?是的,一个矩阵乘以它的逆矩阵后结果为一个单位矩阵

AA^{-1}=I

单位矩阵的定义如下,我们以一个 3x3 的矩阵举例:

\left[ \begin{matrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \right]

13. 传感器融合简单介绍

现在我们的卡尔曼滤波方程完全以矩阵形式表示,

模型:
x_k=Ax_{k-1}+Bu_k
z_k=Cx_k+v_k

预测:

\hat x_k=A\hat x_{k-1}+Bu_k
P_k=AP_{k-1}A^T

更新:

G_k=P_kC^T(CP_kC^T+R)^{-1}
\hat x_k=\hat x_{k}+G_k(z_k-C\hat x_k)
P_k\gets(I-G_kC)P_k

如果想在状态变量中添加以下额外的项似乎是一项很艰巨的任务。卡尔曼滤波最有价值的地方就是用在传感器的融合中。

回到飞机飞行的例子中,飞行员可以访问更多的信息,不仅仅是飞机的高度,仪表上可以显示飞机的水平速度、航向、经纬度、室外温湿度等信息。想象一下,假如一架飞机只有三个传感器,每一个传感器都有给定的状态,一个用于测量高度的气压计,一个用于航向的电子罗盘,一个测量空速的空速指示器等。

假设这些传感器是完全准确的,即没有任何噪声。那么我们的观测方程为:

z_k=Cx_{k-1}+v_k

将会变为:

\left[ \begin{matrix} barometer_k \\ compass_k \\ pitot_k \end{matrix} \right] = \left[ \begin{matrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \right] \left[ \begin{matrix} altitude_{k-1}\\ heading_{k-1} \\ airspeed_{k-1} \end{matrix} \right]

假如我们向系统中加入了一个GPS,同样用于测量高度,那么我们的观测方程将会变为:

\left[ \begin{matrix} barometer_k \\ compass_k \\ pitot_k \\ gps_k \end{matrix} \right] = \left[ \begin{matrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{matrix} \right] \left[ \begin{matrix} altitude_{k-1}\\ heading_{k-1} \\ airspeed_{k-1} \end{matrix} \right]

因此可以将多个传感器的读数结合起来,从而推断出某个状态的信息。就比如我们看病一样,对于重要的疾病总是寻求多个医生的诊断,对于一些重要的事情,最好有多个信息来源。

本篇先介绍这些,下一篇会以一个真实的传感器融合实例开始,并最终实现扩展卡尔曼滤波。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,332评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,508评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,812评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,607评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,728评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,919评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,071评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,802评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,256评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,576评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,712评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,389评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,032评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,026评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,473评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,606评论 2 350

推荐阅读更多精彩内容