卡尔曼滤波通俗介绍
易于理解的介绍,应该是属于文字逻辑,而不是公式逻辑
参考文献
建议按照顺序阅读
【1】【易懂教程】我是如何十分钟理解与推导贝叶斯滤波(Bayes Filter)算法?
【4】详解卡尔曼滤波器
卡尔曼滤波的作用
-
卡尔曼滤波用于优化我们感兴趣的量,当这些量无法直接测量但可以间接测量时。
-
用于估算系统状态,通过组合各种受噪音的传感器测量值
从贝叶斯滤波出发
本部分并不需要真正的了解贝叶斯滤波,只需要理解卡尔曼滤波和它的关系,更清晰的理解卡尔曼滤波
贝叶斯滤波的工作就是根据不断接收到的新信息和提供的一些已经知道的统计值,来不断更新概率,更新概率值的方法是根据概率论中的条件概率计算公式。
贝叶斯滤波公式如下:
P
(
x
t
∣
z
t
,
u
t
,
x
t
−
1
)
=
η
P
(
z
t
∣
x
t
)
P
(
x
t
∣
u
t
,
x
t
−
1
)
P(x_t|z_t,u_t,x_{t-1})=\eta P(z_t|x_t)P(x_t|u_t,x_{t-1})
P(xt∣zt,ut,xt−1)=ηP(zt∣xt)P(xt∣ut,xt−1)
在这个公式中,有两个状态估计假设
- 当前状态 x t x_t xt只与上个时刻状态 x t − 1 x_{t-1} xt−1以及控制命令 u t u_t ut有关(可以推算出 x t x_t xt的预估值)
- 当前的观测量 z t z_t zt只与当前状态 x t x_t xt有关(可以推算出 x t x_t xt的观测值)
贝叶斯滤波是一种思想,提出了如何利用观测值和预估值推算出 x t x_t xt的预测值的可信度。
如果你对观测值,预估值,预测值还有点疑惑,没关系,我们还没有正式介绍它们。
而卡尔曼滤波是对贝叶斯滤波的一种具体实现,具体实现的方式是
- 状态量是线性变化
- 变量和噪声都是正态分布
- 只要知道 μ \mu μ, σ 2 \sigma^2 σ2,则有关变量和噪声的所有信息都知道
- 正太分布经过线性变换依然是线性(刚好与状态线性变换假设对应)
正式进入卡尔曼滤波
还记得卡尔曼滤波的假设吗?所有变量和噪声都是正太分布。
那么显然,对于正太分布
x
∼
N
(
μ
,
σ
2
)
x\sim N(\mu,\sigma^2)
x∼N(μ,σ2),
x
x
x最有可能的值是
μ
\mu
μ,当然,这个数据并不一定靠谱,因为只是最有可能,而不是绝对,那么对于用
μ
\mu
μ代表
x
x
x就有一个可信度,它就是
σ
2
\sigma^2
σ2。
μ
可
能
值
σ
2
可
信
度
\mu \qquad可能值\\ \sigma^2 \qquad可信度
μ可能值σ2可信度
这里值的一提的是,我们往往用向量、矩阵去存储数据、描述状态、代表变化过程等,所以对于
x
∼
N
(
μ
,
σ
2
)
x\sim N(\mu,\sigma^2)
x∼N(μ,σ2),我们可能会有点小变化,比如
σ
2
\sigma^2
σ2变成了协方差矩阵
Σ
\Sigma
Σ。当然这些在概念理解时无伤大雅,只是在公式理解时会出现。
μ
→
可
能
值
Σ
可
信
度
\overrightarrow{\mu } \qquad可能值\\ \Sigma \qquad可信度
μ可能值Σ可信度
好的,我们开始吧。
预估过程
对于一个系统,我们现在知道了它上个时刻的状态
x
^
k
−
1
\hat{x}_{k-1}
x^k−1和对应的可信度
P
k
−
1
P_{k-1}
Pk−1
μ
→
可
能
值
x
^
k
−
1
Σ
可
信
度
P
k
−
1
\overrightarrow{\mu } \qquad可能值 \qquad \hat{x}_{k-1}\\ \Sigma \qquad可信度 \qquad P_{k-1}
μ可能值x^k−1Σ可信度Pk−1
这里值的一提的是, ^ \quad \hat{}\quad ^修饰的向量,表示我们已知(无论是以怎么样的方式得到的,如测量,预估等)。而其对应的可信度,往往是我们卡尔曼滤波需要调的参数
很显然的是,我们对状态的变化有一个数学模型的描述,这个数学模型能够由上一时刻的状态(
x
k
−
1
x_{k-1}
xk−1),推出(也叫预估)这一时刻的状态(
x
k
x_k
xk),如在最普通的卡尔曼滤波中,这个数学模型是线性的。
x
^
k
=
A
x
^
k
−
1
+
b
→
\hat{x}_k=A\hat{x}_{k-1}+\overrightarrow{b }
x^k=Ax^k−1+b
当然还有很多其他模型,如对于一个运动系统,可以是匀速直线运动,也可以是匀加速运动等,对于这个数学模型的预估过程(使用数学模型,由上一时刻的状态(
x
k
−
1
x_{k-1}
xk−1),推出这一时刻的状态(
x
k
x_k
xk)),我们用状态转移矩阵
F
k
F_k
Fk统一表示。
x
^
k
−
1
→
预
估
过
程
F
k
x
^
t
−
=
F
k
x
^
k
−
1
P
k
−
1
→
预
估
过
程
F
k
P
k
−
=
F
k
P
k
−
1
F
k
T
\hat{x}_{k-1}\xrightarrow{预估过程F_k}\hat{x}_t^{-}=F_k\hat{x}_{k-1}\\ P_{k-1}\xrightarrow{预估过程F_k}P_k^{-}=F_kP_{k-1}F_k^T
x^k−1预估过程Fkx^t−=Fkx^k−1Pk−1预估过程FkPk−=FkPk−1FkT
有对于 x ^ t − \hat{x}_t^{-} x^t−这样在右上角有 − \quad ^- \quad −的数据,表示预估值,代表着单纯由我们的数学模型推测出来的值,但并不是卡尔曼滤波真正得到的值,只是一个过程量。
对于这个预估过程很好理解,一切都是我们知道的!(牢记卡尔曼滤波的假设!),因此无论是状态向量还是可信度,他们的预估值都是完全在我们掌握之下的,因此无需添加一些量进行修正。
控制过程
对于一个系统,我们会进行一些控制使其状态发生变化,像开车时,我们一直匀速运动,那么预估过程也会是匀速模型,而我们突然踩油门,会有一个加速的控制,这些是我们数学模型没有掌握的,因此我们需要添加控制量对应的项。我们掌握了很多具体控制,如踩多少油门,这些具体控制用
B
^
k
\hat{B}_{k}
B^k表示控制过程,
u
^
t
\hat{u}_t
u^t表示控制量。
x
^
k
−
=
F
k
x
^
k
−
1
→
控
制
过
程
B
^
k
,
控
制
量
u
^
k
x
^
k
−
=
F
k
x
^
k
−
1
+
B
k
u
^
k
−
P
k
−
=
F
k
P
k
−
1
F
k
T
→
控
制
过
程
B
^
k
,
控
制
量
u
^
t
P
k
−
=
F
k
P
k
−
1
F
k
T
\hat{x}_k^{-}=F_k\hat{x}_{k-1}\xrightarrow{控制过程\hat{B}_{k},控制量\hat{u}_k}\hat{x}_k^{-}=F_k\hat{x}_{k-1}+B_k\hat{u}_k^{-}\\ P_k^{-}=F_kP_{k-1}F_k^T\xrightarrow{控制过程\hat{B}_{k},控制量\hat{u}_t}P_k^{-}=F_kP_{k-1}F_k^T
x^k−=Fkx^k−1控制过程B^k,控制量u^kx^k−=Fkx^k−1+Bku^k−Pk−=FkPk−1FkT控制过程B^k,控制量u^tPk−=FkPk−1FkT
而**还有很多我们没有掌握的“控制”**会对汽车的运行造成影响,如风速、车胎磨损、路面情况等,这些我们没有掌握但是造成影响的“控制”,我们称其为噪声(干扰),它们对我们状态量的改变未知(因为我们对它们没有掌握!),因此会对我们修正的预估值的可信度产生影响,这个影响称为
Q
k
Q_k
Qk。注意,我们的状态量只是最可能的量,因此不一定是准确的。(牢记卡尔曼滤波的假设!)
x
^
k
−
=
F
k
x
^
k
−
1
+
B
k
u
^
k
−
→
噪
声
Q
k
x
^
k
−
=
F
k
x
^
k
−
1
+
B
k
u
^
k
−
P
k
−
=
F
k
P
k
−
1
F
k
T
→
噪
声
Q
k
P
k
−
=
F
k
P
k
−
1
F
k
T
+
Q
k
\hat{x}_k^{-}=F_k\hat{x}_{k-1}+B_k\hat{u}_k^{-}\xrightarrow{噪声Q_k}\hat{x}_k^{-}=F_k\hat{x}_{k-1}+B_k\hat{u}_k^{-}\\ P_k^{-}=F_kP_{k-1}F_k^T\xrightarrow{噪声Q_k}P_k^{-}=F_kP_{k-1}F_k^T+Q_k
x^k−=Fkx^k−1+Bku^k−噪声Qkx^k−=Fkx^k−1+Bku^k−Pk−=FkPk−1FkT噪声QkPk−=FkPk−1FkT+Qk
更新过程
以上是我们完全根据上一个状态的预估,但实际上,系统的状态可能在这个时刻突然发生改变(注意与控制过程造成状态改变的区别)
- 控制过程造成的状态改变往往是平滑的,即使有噪声的影响,也不会改变整体曲线的平滑改变的趋势。
- 这里说的系统状态的改变,指的是陡峭尖锐的曲线变化。
举一个例子:原本我们正常的在路上驱车行使,刚开始是匀速行驶(这时我们可以给自己车建立一个数学模型了——线性的变换;实际上,绝大部分过程的数学模型都可设置成线性的,因为时间间隔足够小,以直代曲的思想),这时我们进入了预估过程;我们需要超车,于是踩油门,因此进入了控制过程;不幸的是在 t t t时刻,我们超车失败,撞在了装了30吨货物的大卡车上,速度瞬间变成0,发动机也报废,而因为这只是瞬间发生的事件,因此我们仍然在踩油门,即预估过程和控制过程一直在继续,但是实际上我们已经完全停下来了,但是我们的滤波没有即时跟上我们的状态变化。
我们希望我们的滤波能够及时的反映一些意外,因此我们需要传感器数据对卡尔曼滤波进行修正(更新),这是更新过程。
我们有很多传感器可以帮助我们及时跟进系统的状态变化,对我们的预估值进行更新。依然是关于汽车的例子,比如我们有发动机的温度计和GPS两个传感器,对发动机的温度( T T T)和汽车速度( v x , v y v_x,v_y vx,vy)进行观测(这时候 x t → \overrightarrow{x_t} xt是三维向量 ( T , v x , v y ) T (T,v_x,v_y)^T (T,vx,vy)T,这里上标T表示转置),我们会有一个2X3的矩阵,把
我们有m个传感器,会对n个状态值进行观测(
x
^
t
\hat{x}_t
x^t是n维向量),那我们就会有一个
m
×
n
m\times n
m×n的矩阵,表示由状态的预估值向传感器数据的预估值的转换矩阵,这个矩阵是
H
k
H_k
Hk,而传感器的预估值为
z
^
k
−
\hat{z}_k^-
z^k−,它的可信度是
R
k
−
R_k^-
Rk−。
x
^
k
−
→
由
状
态
量
到
传
感
器
数
据
的
推
算
过
程
H
k
z
^
k
−
=
H
k
x
^
k
−
P
k
−
→
H
k
R
k
−
=
H
k
P
k
H
k
T
\hat{x}_k^{-}\xrightarrow{由状态量到传感器数据的推算过程H_{k}}\hat{z}_k^{-}=H_k\hat{x}_{k}^-\\ P_k^{-}\xrightarrow{\qquad\qquad H_{k}\qquad\qquad}R_k^{-}=H_kP_{k}H_k^T
x^k−由状态量到传感器数据的推算过程Hkz^k−=Hkx^k−Pk−HkRk−=HkPkHkT
而我们还有
t
t
t时刻的传感器测量值为
z
→
k
\overrightarrow{z}_k
zk,它的可信度是
R
k
R_k
Rk。那么我们该怎么样用
z
→
k
\overrightarrow{z}_k
zk,
R
k
R_k
Rk更新
z
^
k
−
\hat{z}_k^-
z^k−,
R
k
−
R_k^-
Rk−呢?
我们的基本思路是 z → k \overrightarrow{z}_k zk相信一些, z ^ k − \hat{z}_k^- z^k−也相信一些,在 z → k \overrightarrow{z}_k zk和 z ^ k − \hat{z}_k^- z^k−之间取值,即取 z → k \overrightarrow{z}_k zk和 z ^ k − \hat{z}_k^- z^k−的加权平均数。有了这个想法后,有一个很简单的算法:相乘!(不要忘记我们的假设:都是高斯分布!)
N ( μ 1 , Σ 1 ) N(\mu_1,\Sigma_1) N(μ1,Σ1)与 N ( μ 2 , Σ 2 ) N(\mu_2,\Sigma_2) N(μ2,Σ2)相乘,结果还是高斯分布 N ( μ , Σ ) N(\mu,\Sigma) N(μ,Σ)!高斯万岁!
μ = μ 1 Σ 2 + μ 2 Σ 1 Σ 1 + Σ 2 Σ = Σ 1 Σ 2 Σ 1 + Σ 2 \mu=\frac{\mu_1\Sigma_2+\mu_2\Sigma_1}{\Sigma_1+\Sigma_2} \qquad \Sigma=\frac{\Sigma_1\Sigma_2}{\Sigma_1+\Sigma_2} μ=Σ1+Σ2μ1Σ2+μ2Σ1Σ=Σ1+Σ2Σ1Σ2
可以化简为
μ = μ 1 + k ( μ 2 − μ 1 ) Σ = Σ 1 − k Σ 1 \mu=\mu_1+k(\mu_2-\mu_1) \qquad \Sigma=\Sigma_1-k\Sigma_1 μ=μ1+k(μ2−μ1)Σ=Σ1−kΣ1
其中 k = Σ 1 Σ 1 + Σ 2 k=\frac{\Sigma_1}{\Sigma_1+\Sigma_2} k=Σ1+Σ2Σ1,这里有一个比较重要的问题是 Σ 1 , Σ 2 {\Sigma_1,\Sigma_2} Σ1,Σ2是协方差矩阵,没有“倒数”的概念,应该是“逆”,但是这样写可以在更直观的理解加权平均数的概念,所以仅仅在理解上,我们做以下的对应:
( Σ 1 + Σ 2 ) − 1 = 1 Σ 1 + Σ 2 (\Sigma_1+\Sigma_2)^{-1}=\frac{1}{\Sigma_1+\Sigma_2} (Σ1+Σ2)−1=Σ1+Σ21
我们相乘,可得
z
^
k
′
=
z
^
k
−
+
K
(
z
→
k
−
z
^
k
−
)
R
k
′
=
R
k
−
−
K
R
k
−
K
=
H
k
P
k
H
k
T
H
k
P
k
H
k
T
+
R
k
\hat{z}_k^{'}=\hat{z}_k^-+K(\overrightarrow{z}_k-\hat{z}_k^-) \qquad R_k^{'}=R_k^--KR_k^-\\ K= \frac{H_kP_kH_k^T}{H_kP_kH_k^T+R_k}
z^k′=z^k−+K(zk−z^k−)Rk′=Rk−−KRk−K=HkPkHkT+RkHkPkHkT
最后,我们需要推算出k时刻的最优值
z
^
k
′
=
z
^
k
−
+
K
(
z
→
k
−
z
^
k
−
)
→
由
传
感
器
数
据
到
状
态
量
的
推
算
过
程
H
k
−
x
^
k
=
H
k
−
1
z
^
k
′
=
x
^
k
−
+
K
′
(
z
→
k
−
H
k
−
1
x
^
k
′
)
R
k
′
=
R
k
−
−
K
R
k
−
→
H
k
−
P
k
=
H
k
−
1
R
k
′
H
k
T
−
1
=
P
k
−
−
K
′
H
k
P
k
−
K
′
=
P
k
H
k
T
H
k
P
k
H
k
T
+
R
k
\hat{z}_k^{'}=\hat{z}_k^-+K(\overrightarrow{z}_k-\hat{z}_k^-) \xrightarrow{由传感器数据到状态量的推算过程H_{k}^-}\hat{x}_k=H_k^{-1}\hat{z}_k^{'}=\hat{x}_k^-+K^{'}(\overrightarrow{z}_k-H_k^{-1}\hat{x}_k^{'})\\ R_k^{'}=R_k^--KR_k^-\xrightarrow{\qquad\qquad H_{k}^-\qquad\qquad}P_k=H_k^{-1}R_k^{'}{H_k^{T}}^{-1}=P_k^--K^{'}H_kP_k^{-}\\ K^{'}=\frac{P_kH_k^T}{H_kP_kH_k^T+R_k}
z^k′=z^k−+K(zk−z^k−)由传感器数据到状态量的推算过程Hk−x^k=Hk−1z^k′=x^k−+K′(zk−Hk−1x^k′)Rk′=Rk−−KRk−Hk−Pk=Hk−1Rk′HkT−1=Pk−−K′HkPk−K′=HkPkHkT+RkPkHkT
对于
x
^
k
−
\hat{x}_k^-
x^k−
x
^
k
=
(
1
−
K
′
H
k
)
x
^
k
−
+
K
′
H
k
x
→
k
\hat{x}_k=(1-K^{'}H_k)\hat{x}_k^-+K^{'}H_k\overrightarrow{x}_k\\
x^k=(1−K′Hk)x^k−+K′Hkxk
注意:这里有一步
z
→
k
=
H
K
x
→
k
\overrightarrow{z}_k=H_K\overrightarrow{x}_k
zk=HKxk的变换
所以,我们可以直观的看出:卡尔曼滤波是一种权重不断自我更新的加权平均数
复习
请给下面的图填写完整吧!(实际上是我的绘图软件没法打数学公式)