一、状态估计概念
由于位姿和路标点都是待估计的变量,我们改变一下记号,令 x k {\boldsymbol{x}_{k}} xk为 k k k时刻的所有未知量。它包含了当前时刻的相机位姿与 m m m个路标点。在这种记号的意义下写成:
x k ≜ { x k , y 1 , . . . , y m } (1) {{x}_{k}}\triangleq \left\{ \boldsymbol{{x}_{k}},\text{ }\boldsymbol{{y}_{1}},\text{ }.\text{ }.\text{ }.\text{ },\text{ }{\boldsymbol{y}_{m}} \right\}\tag{1} xk≜{xk, y1, . . . , ym}(1)
把 k k k时刻的所有观测记作 z k {\boldsymbol{z}_{k}} zk:
{ x k = f ( x k − 1 , u k ) + w k z k = h ( x k ) + v k k = 1 , … , N (2) \left\{\begin{array}{l} \boldsymbol{x}_k=f\left(\boldsymbol{x}_{k-1}, \boldsymbol{u}_k\right)+\boldsymbol{w}_k \\ \boldsymbol{z}_k=h\left(\boldsymbol{x}_k\right)+\boldsymbol{v}_k \end{array} \quad k=1, \ldots, N\right. \tag{2} {xk=f(xk−1,uk)+wkzk=h(xk)+vkk=1,…,N(2)
现在考虑第 k k k时刻的情况。我们希望用过去0到 k k k中的数据,来估计现在的状态分布:
P ( x k ∣ x 0 , u 1 : k , z 1 : k ) (3) P\left( \boldsymbol{{x}_{k}}|{\boldsymbol{x}_{0}},{\boldsymbol{u}_{1:k}},{\boldsymbol{z}_{1:k}} \right)\tag{3} P(xk∣x0,u1:k,z1:k)(3)
下标 0 : k 0:k 0:k 表示从 0 时刻到 k k k时刻的所有数据。请注意 z k {{z}_{k}} zk来表达所有在 k k k时刻的观测数据,注意它可能不止一个,只是这种记法更加方便。 下面我们来看如何对状态进行估计。按照 Bayes 法则,把 z k { \boldsymbol{z}_{k}} zk与 x k { \boldsymbol{x}_{k}} xk交换位置,有:
P
(
x
k
∣
x
0
,
u
1
:
k
,
z
1
:
k
)
∝
P
(
z
k
∣
x
k
)
P
(
x
k
∣
x
0
,
u
1
:
k
,
z
1
:
k
−
1
)
(4)
P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k}\right) \propto P\left(\boldsymbol{z}_k \mid \boldsymbol{x}_k\right) P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k-1}\right)\tag{4}
P(xk∣x0,u1:k,z1:k)∝P(zk∣xk)P(xk∣x0,u1:k,z1:k−1)(4)
这里第一项称为似然,第二项称为先验。似然由观测方程给
定,而先验部分,我们要明白当前状态
x
k
{{x}_{k}}
xk是基于过去所有的状态估计得来的。按照
x
k
−
1
{{x}_{k-1}}
xk−1为条件概率展开:
P ( x k ∣ x 0 , u 1 : k , z 1 : k − 1 ) = ∫ P ( x k ∣ x k − 1 , x 0 , u 1 : k , z 1 : k − 1 ) P ( x k − 1 ∣ x 0 , u 1 : k , z 1 : k − 1 ) d x k − 1 (5) P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k-1}\right)=\int P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_{k-1}, \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k-1}\right) P\left(\boldsymbol{x}_{k-1} \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k-1}\right) \mathrm{d} \boldsymbol{x}_{k-1}\tag{5} P(xk∣x0,u1:k,z1:k−1)=∫P(xk∣xk−1,x0,u1:k,z1:k−1)P(xk−1∣x0,u1:k,z1:k−1)dxk−1(5)
如果我们考虑更久之前的状态,也可以继续对此式进行展开,但现在我们假设了马尔可夫性,所以只关心 k k k时刻和 k − 1 k-1 k−1时刻的情况。
x 0 ⊕ u 1 → x 1 {{x}_{0}}\oplus {{u}_{1}}\to {{x}_{1}} x0⊕u1→x1
x 1 ⊕ u 2 → x 2 {{x}_{1}}\oplus {{u}_{2}}\to {{x}_{2}} x1⊕u2→x2
由于假设了马尔可夫性,当前时刻状态只和上一个时刻的状态
x
k
−
1
{\boldsymbol{x}_{k-1}}
xk−1和输入
u
k
{\boldsymbol{u}_{k}}
uk有关,式(5)中等式右侧第一部分可进一步简化
P
(
x
k
∣
x
k
−
1
,
x
0
,
u
1
:
k
,
z
1
:
k
−
1
)
=
P
(
x
k
∣
x
k
−
1
,
u
k
)
(6)
P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_{k-1}, \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k-1}\right)=P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_{k-1}, \boldsymbol{u}_k\right)\tag{6}
P(xk∣xk−1,x0,u1:k,z1:k−1)=P(xk∣xk−1,uk)(6)
第二部分可简化为:
P
(
x
k
−
1
∣
x
0
,
u
1
:
k
,
z
1
:
k
−
1
)
=
P
(
x
k
−
1
∣
x
0
,
u
1
:
k
−
1
,
z
1
:
k
−
1
)
(7)
P\left(\boldsymbol{x}_{k-1} \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k-1}\right)=P\left(\boldsymbol{x}_{k-1} \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k-1}, \boldsymbol{z}_{1: k-1}\right)\tag{7}
P(xk−1∣x0,u1:k,z1:k−1)=P(xk−1∣x0,u1:k−1,z1:k−1)(7)
这是考虑到 k k k时刻的输入量 u k {\boldsymbol{u}_{k}} uk与 k − 1 k-1 k−1时刻的状态 x k − 1 {\boldsymbol{x}_{k-1}} xk−1无关,所以我们把 u k {\boldsymbol{u}_{k}} uk拿掉。
这一项实际是 k − 1 k-1 k−1时刻的 x k − 1 {\boldsymbol{x}_{k-1}} xk−1状态分布。直觉上看,就是用 x 0 {\boldsymbol{x}_{0}} x0、一系列的动作输入 u 1 : k {\boldsymbol{u}_{1:k}} u1:k和一系列观测 z 1 : k − 1 {\boldsymbol{z}_{1:k-1}} z1:k−1估计得到的 x k − 1 {\boldsymbol{x}_{k-1}} xk−1。
这样就可以递推:
x 0 ⊕ u 1 → x 1 {{x}_{0}}\oplus {{u}_{1}}\to {{x}_{1}} x0⊕u1→x1
x 1 ⊕ u 2 → x 2 {{x}_{1}}\oplus {{u}_{2}}\to {{x}_{2}} x1⊕u2→x2
于是, ( 6 ) (6) (6)和 ( 7 ) (7) (7)的到 ( 5 ) (5) (5)的过程说明了,“如何把 k − 1 k-1 k−1时刻的状态分布推导至 k k k时刻”这样一件事(也就是后面要说的预测部分)。也就是说,在程序运行期间,我们只要维护一个状态量,对它进行不断地迭代和更新即可。进一步,如果假设状态量 x \boldsymbol{x} x服从高斯分布,因为高斯分布可以由一个均值和一个方差唯一确定,那我们只需考虑维护状态量的均值和协方差即可。
如果我们假设系统为线性系统,也就是该系统可以由下面的线性方程来描述:
{
x
k
=
A
k
x
k
−
1
+
u
k
+
w
k
z
k
=
C
k
x
k
+
v
k
k
=
1
,
…
,
N
(8)
\left\{\begin{array}{l} \boldsymbol{x}_k=\boldsymbol{A}_k \boldsymbol{x}_{k-1}+\boldsymbol{u}_k+\boldsymbol{w}_k \\ \boldsymbol{z}_k=\boldsymbol{C}_k \boldsymbol{x}_k+\boldsymbol{v}_k \end{array} \quad k=1, \ldots, N\right. \tag{8}
{xk=Akxk−1+uk+wkzk=Ckxk+vkk=1,…,N(8)
A k \boldsymbol{A}_k Ak就是状态转移矩阵,有的也会为输入增加控制矩阵 B B B
{ x k = A k x k − 1 + B k u k + w k z k = C k x k + v k k = 1 , … , N (9) \left\{\begin{array}{l} \boldsymbol{x}_k=\boldsymbol{A}_k \boldsymbol{x}_{k-1}+\boldsymbol{B}_k \boldsymbol{u}_k+\boldsymbol{w}_k \\ \boldsymbol{z}_k=\boldsymbol{C}_k \boldsymbol{x}_k+\boldsymbol{v}_k \end{array} \quad k=1, \ldots, N\right. \tag{9} {xk=Akxk−1+Bkuk+wkzk=Ckxk+vkk=1,…,N(9)
假设所有的状态和噪声均满足高斯分布。记这里的噪声服从零均值高斯分布:
w k ∼ N ( 0 , R ) . v k ∼ N ( 0 , Q ) (10) {w_k} \sim N({\bf{0}},R).\quad {v_k} \sim N({\bf{0}},Q)\tag{10} wk∼N(0,R).vk∼N(0,Q)(10)
为了简洁我省略了 R R R和 Q Q Q的下标。现在,利用马尔可夫性,假设我们知道了 k − 1 {k-1} k−1时刻的后验(在 k − 1 {k-1} k−1时刻看来)状态估计: x ^ k − 1 \hat{\boldsymbol{x}}_{k−1} x^k−1 和它的协方差 P ^ k − 1 \hat{\boldsymbol{P}}_{k−1} P^k−1,现在要根据 k − 1 {k-1} k−1 时刻的输入和观测数据,确定 xk 的后验分布。
为区分推导中的先验和后验,我们在记号上作一点区别:以上帽子 x ^ k \hat{\boldsymbol{x}}_{k} x^k表示后验,以下帽子 x ˇ k \check{\boldsymbol{x}}_{k} xˇk表示先验分布:
二、推导
总体的思想:
先验就是预测值,不准的值。后验就是修正过的值,精确值。
有了上一时刻的后验(先别管怎么得到的,反正就是有了)。通过运动方程可以大概地预测出下一时刻的值,这个值不准,因为通过运动的性质进行估计会有误差,所以叫做先验值。同时我们还有观测,可以用观测到的信息修正这个先验值,修正之后得到了更精确的值叫做后验。有了后验,那么下一时刻就可以重复上面的过程,用运动方程预测,用观测方程修正
1、通过运动方程确定先验
通过运动方程确定 x k \boldsymbol{x}_{k} xk的先验分布。根据高斯分布的性质,显然有:
P ( x k ∣ x 0 , u 1 : k , z 1 : k − 1 ) = N ( A k x ^ k − 1 + u k , A k P ^ k − 1 A k T + R ) (11) P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k-1}\right)=N\left(\boldsymbol{A}_k \hat{\boldsymbol{x}}_{k-1}+\boldsymbol{u}_k, \boldsymbol{A}_k \hat{\boldsymbol{P}}_{k-1} \boldsymbol{A}_k^T+\boldsymbol{R}\right) \tag{11} P(xk∣x0,u1:k,z1:k−1)=N(Akx^k−1+uk,AkP^k−1AkT+R)(11)
这一步称为预测,它显示了如何从上一个时刻的状态,根据输入信息(但是有噪声),推断当前时刻的状态分布,也就是根据上一时刻的后验 x ^ k − 1 \hat{\boldsymbol{x}}_{k−1} x^k−1 和它的协方差 P ^ k − 1 \hat{\boldsymbol{P}}_{k−1} P^k−1,使用运动方程来预测下一时刻的先验 x ˇ k \check{\boldsymbol{x}}_{k} xˇk 和它的协方差$\check{\boldsymbol{P}}_{k}
注意这里的 P ( x k ∣ x 0 , u 1 : k , z 1 : k − 1 ) P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{1: k-1}\right) P(xk∣x0,u1:k,z1:k−1)使用了 k k k时刻的动作和 k − 1 k-1 k−1时刻以及之前的观测,并没有利用到 k {k} k时刻的观测,所以是不准的,是个先验的概率。
这个推导是通过线性高斯系统的性质得到的。
记先验的概率分布为:
x
ˇ
k
=
A
k
x
^
k
−
1
+
u
k
,
P
ˇ
k
=
A
k
P
^
k
−
1
A
k
T
+
R
(12)
\check{\boldsymbol{x}}_k=\boldsymbol{A}_k \hat{\boldsymbol{x}}_{k-1}+\boldsymbol{u}_k, \quad \check{\boldsymbol{P}}_k=\boldsymbol{A}_k \hat{\boldsymbol{P}}_{k-1} \boldsymbol{A}_k^{\mathrm{T}}+\boldsymbol{R} \tag{12}
xˇk=Akx^k−1+uk,Pˇk=AkP^k−1AkT+R(12)
2、通过观测方程修正先验,得到后验
由观测方程,我们可以计算在某个状态下,应该产生怎样的观测数据:
P ( z k ∣ x k ) = N ( C k x k , Q ) (13) P\left(\boldsymbol{z}_k \mid \boldsymbol{x}_k\right)=N\left(\boldsymbol{C}_k \boldsymbol{x}_k, \boldsymbol{Q}\right) \tag{13} P(zk∣xk)=N(Ckxk,Q)(13)
后验概率也就是式 ( 13 ) (13) (13)和式 ( 11 ) (11) (11)的乘积,也就是由式 ( 4 ) (4) (4)给出的贝叶斯公式。直接相乘非常复杂,但是我们可以利用高斯分布和高斯分布的乘积仍然是一个高斯分布这一性质,先假设结果符合的高斯分布 x k ∼ N ( x ^ k , P ^ k ) {x_k} \sim N(\hat{\boldsymbol{x}}_{k},\hat{\boldsymbol{P}}_{k}) xk∼N(x^k,P^k) 那么:
N ( x ^ k , P ^ k ) = η N ( C k x k , Q ) ⋅ N ( x ˇ k , P ˇ k ) (14) N\left(\hat{\boldsymbol{x}}_k, \hat{\boldsymbol{P}}_k\right)=\eta N\left(\boldsymbol{C}_k \boldsymbol{x}_k, \boldsymbol{Q}\right) \cdot N\left(\check{\boldsymbol{x}}_k, \check{\boldsymbol{P}}_k\right) \tag{14} N(x^k,P^k)=ηN(Ckxk,Q)⋅N(xˇk,Pˇk)(14)
把 ( 15 ) (15) (15)左右两侧取 l o g ( ∗ ) log(*) log(∗),系数项变为加上的常数项(与 x k \boldsymbol{x}_k xk无关),指数项相乘变为指数相加,那么指数部分相等:
( x k − x ^ k ) T P ^ k − 1 ( x k − x ^ k ) = ( z k − C k x k ) T Q − 1 ( z k − C k x k ) + ( x k − x ˇ k ) T P ˇ k − 1 ( x k − x ˇ k ) (15) \left(\boldsymbol{x}_k-\hat{\boldsymbol{x}}_k\right)^{\mathrm{T}} \hat{\boldsymbol{P}}_k^{-1}\left(\boldsymbol{x}_k-\hat{\boldsymbol{x}}_k\right)=\left(\boldsymbol{z}_k-\boldsymbol{C}_k \boldsymbol{x}_k\right)^{\mathrm{T}} \boldsymbol{Q}^{-1}\left(\boldsymbol{z}_k-\boldsymbol{C}_k \boldsymbol{x}_k\right)+\left(\boldsymbol{x}_k-\check{\boldsymbol{x}}_k\right)^{\mathrm{T}} \check{\boldsymbol{P}}_k^{-1}\left(\boldsymbol{x}_k-\check{\boldsymbol{x}}_k\right) \tag{15} (xk−x^k)TP^k−1(xk−x^k)=(zk−Ckxk)TQ−1(zk−Ckxk)+(xk−xˇk)TPˇk−1(xk−xˇk)(15)
为了求左侧的 x ^ k , P ^ k \hat{\boldsymbol{x}}_{k},\hat{\boldsymbol{P}}_{k} x^k,P^k,我们把两边展开,并比较 x ^ k \hat{\boldsymbol{x}}_{k} x^k的二次和一次系数。对于二次系数,有:
P ^ k − 1 = C k T Q − 1 C k + P ˇ k − 1 (16) \hat{\boldsymbol{P}}_k^{-1}=\boldsymbol{C}_k^T \boldsymbol{Q}^{-1} \boldsymbol{C}_k+\check{\boldsymbol{P}}_k^{-1} \tag{16} P^k−1=CkTQ−1Ck+Pˇk−1(16)
这里定义一个中间变量 K \boldsymbol{K} K,其实也就是卡尔曼增益(这一步推导需要用的SMW恒等式):
K = P ^ k C k T Q − 1 = P ˇ k C k T ( C k P ˇ k C k T + Q ) − 1 (17) \boldsymbol{K} =\hat{\boldsymbol{P}}_k \boldsymbol{C}_k^T \boldsymbol{Q}^{-1} =\check{\boldsymbol{P}}_k \boldsymbol{C}_k^T\left(\boldsymbol{C}_k \check{\boldsymbol{P}}_k \boldsymbol{C}_k^T+\boldsymbol{Q}\right)^{-1} \tag{17} K=P^kCkTQ−1=PˇkCkT(CkPˇkCkT+Q)−1(17)
在式 ( 16 ) (16) (16)左右两侧乘 P ^ k \hat{\boldsymbol{P}}_{k} P^k 有:
I = P ^ k C k T Q − 1 C k + P ^ k P ˇ k − 1 = K C k + P ^ k P ˇ k − 1 (18) \boldsymbol{I}=\hat{\boldsymbol{P}}_k \boldsymbol{C}_k^T \boldsymbol{Q}^{-1} \boldsymbol{C}_k+\hat{\boldsymbol{P}}_k \check{\boldsymbol{P}}_k^{-1}=\boldsymbol{K} \boldsymbol{C}_k+\hat{\boldsymbol{P}}_k \check{\boldsymbol{P}}_k^{-1} \tag{18} I=P^kCkTQ−1Ck+P^kPˇk−1=KCk+P^kPˇk−1(18)
于是有后验协方差 P ^ k \hat{\boldsymbol{P}}_k P^k:
P ^ k = ( I − K C k ) P ˇ k (19) \hat{\boldsymbol{P}}_k=\left(\boldsymbol{I}-\boldsymbol{K} \boldsymbol{C}_k\right) \check{\boldsymbol{P}}_k \tag{19} P^k=(I−KCk)Pˇk(19)
然后再算后验的均值比较一次项的系数,有:
− 2 x ^ k T P ^ k − 1 x k = − 2 z k T Q − 1 C k x k − 2 x ˇ k T P ˇ k − 1 x k (20) -2 \hat{\boldsymbol{x}}_k^T \hat{\boldsymbol{P}}_k^{-1} \boldsymbol{x}_k=-2 \boldsymbol{z}_k^T \boldsymbol{Q}^{-1} \boldsymbol{C}_k \boldsymbol{x}_k-2 \check{\boldsymbol{x}}_k^T \check{\boldsymbol{P}}_k^{-1} \boldsymbol{x}_k \tag{20} −2x^kTP^k−1xk=−2zkTQ−1Ckxk−2xˇkTPˇk−1xk(20)
利用协方差矩阵对称并且对称矩阵的逆矩阵也对称的性质整理得:
P
^
k
−
1
x
^
k
=
C
k
T
Q
−
1
z
k
+
P
ˇ
k
−
1
x
ˇ
k
(21)
\hat{\boldsymbol{P}}_k^{-1} \hat{\boldsymbol{x}}_k=\boldsymbol{C}_k^T \boldsymbol{Q}^{-1} \boldsymbol{z}_k+\check{\boldsymbol{P}}_k^{-1} \check{\boldsymbol{x}}_k \tag{21}
P^k−1x^k=CkTQ−1zk+Pˇk−1xˇk(21)
两侧乘以
P
^
k
\hat{\boldsymbol{P}}_{k}
P^k并代入式(17),得到了后验均值:
x
^
k
=
P
^
k
C
k
T
Q
−
1
z
k
+
P
^
k
P
ˇ
k
−
1
x
ˇ
k
=
K
z
k
+
(
I
−
K
C
k
)
x
ˇ
k
=
x
ˇ
k
+
K
(
z
k
−
C
k
x
ˇ
k
)
(23)
\begin{aligned} \hat{\boldsymbol{x}}_k & =\hat{\boldsymbol{P}}_k \boldsymbol{C}_k^T \boldsymbol{Q}^{-1} \boldsymbol{z}_k+\hat{\boldsymbol{P}}_k \check{\boldsymbol{P}}_k^{-1} \check{\boldsymbol{x}}_k \\ & =\boldsymbol{K} \boldsymbol{z}_k+\left(\boldsymbol{I}-\boldsymbol{K} \boldsymbol{C}_k\right) \check{\boldsymbol{x}}_k=\check{\boldsymbol{x}}_k+\boldsymbol{K}\left(\boldsymbol{z}_k-\boldsymbol{C}_k \check{\boldsymbol{x}}_k\right) \end{aligned} \tag{23}
x^k=P^kCkTQ−1zk+P^kPˇk−1xˇk=Kzk+(I−KCk)xˇk=xˇk+K(zk−Ckxˇk)(23)
三、总结
总而言之,上面的两个步骤可以归纳为“预测”(Predict)和“更新”(Update)两个步骤:
1、预测
计算先验概率分布(均值和方差)
x
ˇ
k
=
A
k
x
^
k
−
1
+
u
k
,
P
ˇ
k
=
A
k
P
^
k
−
1
A
k
T
+
R
(24)
\check{\boldsymbol{x}}_k=\boldsymbol{A}_k \hat{\boldsymbol{x}}_{k-1}+\boldsymbol{u}_k, \quad \check{\boldsymbol{P}}_k=\boldsymbol{A}_k \hat{\boldsymbol{P}}_{k-1} \boldsymbol{A}_k^T+\boldsymbol{R} \tag{24}
xˇk=Akx^k−1+uk,Pˇk=AkP^k−1AkT+R(24)
2、更新
(1)计算卡尔曼增益
K
k
=
P
ˇ
k
C
k
T
(
C
k
P
ˇ
k
C
k
T
+
Q
)
−
1
(25)
\boldsymbol{K}_k=\check{\boldsymbol{P}}_k \boldsymbol{C}_k^T\left(\boldsymbol{C}_k \check{\boldsymbol{P}}_k \boldsymbol{C}_k^T+\boldsymbol{Q}\right)^{-1} \tag{25}
Kk=PˇkCkT(CkPˇkCkT+Q)−1(25)
(2)计算后验概率分布(均值和方差)
x
^
k
=
K
k
z
k
+
(
I
−
K
k
C
k
)
x
ˇ
k
=
x
ˇ
k
+
K
k
(
z
k
−
C
k
x
ˇ
k
)
(26)
\hat{\boldsymbol{x}}_k =\boldsymbol{K}_k \boldsymbol{z}_k+\left(\boldsymbol{I}-\boldsymbol{K}_k \boldsymbol{C}_k\right) \check{\boldsymbol{x}}_k =\check{\boldsymbol{x}}_k+\boldsymbol{K}_k\left(\boldsymbol{z}_k-\boldsymbol{C}_k \check{\boldsymbol{x}}_k\right) \tag{26}
x^k=Kkzk+(I−KkCk)xˇk=xˇk+Kk(zk−Ckxˇk)(26)
P ^ k = ( I − K k C k ) P ˇ k (27) \hat{\boldsymbol{P}}_k=\left(\boldsymbol{I}-\boldsymbol{K}_k \boldsymbol{C}_k\right) \check{\boldsymbol{P}}_k \tag{27} P^k=(I−KkCk)Pˇk(27)
式$(17)可能会让人有些迷惑,我们使用
P
^
k
\hat{\boldsymbol{P}}_{k}
P^k表示
K
\boldsymbol{K}
K,然后又用
K
\boldsymbol{K}
K 来表示
P
^
k
\hat{\boldsymbol{P}}_{k}
P^k ,
有些循环定义的意思。但是注意,其实式$(16)已经可以算出来
P
^
k
\hat{\boldsymbol{P}}_{k}
P^k 了,然后我们用它计算了
K
\boldsymbol{K}
K ,式
(
17
)
(17)
(17)后面的那个等式就是计算出来的不依靠
P
^
k
\hat{\boldsymbol{P}}_{k}
P^k 表示的最终形式
四、通俗解释
x ^ k = K k z k + ( I − K k C k ) x ˇ k \hat{\boldsymbol{x}}_k =\boldsymbol{K}_k \boldsymbol{z}_k+\left(\boldsymbol{I}-\boldsymbol{K}_k \boldsymbol{C}_k\right) \check{\boldsymbol{x}}_k x^k=Kkzk+(I−KkCk)xˇk
当 C k = 1 \boldsymbol{C}_k=1 Ck=1时也就是观测和状态是同一种数据,比如观测是GPS测量到的车的位置,输入是车的里程计。要估计的量是车的位置。
卡尔曼增益其实就是用来对观测值(GPS测量)和先验值(里程计)做了一个加权平均,只不过这个卡尔曼增益是动态刷新的。
如果更相信观测,那么增大
K
k
{\boldsymbol{K}}_k
Kk 如果更相信运动预测,就要减小
K
k
{\boldsymbol{K}}_k
Kk
五、实际使用(套公式,参数初始化和整定)
(1)选择状态量,观测量,构建运动方程和观测方程
(2)初始化参数
x ^ k \hat{\boldsymbol{x}}_k x^k可以设置为 0 \bf{0} 0 ,因为不知道初始的 x ^ k \hat{\boldsymbol{x}}_k x^k,当然直到的话也可以设置为其他的值
P ^ k − 1 \hat{\boldsymbol{P}}_{k-1} P^k−1可以设置为 1 \bf{1} 1,这个值其实关系不大,因为很快就可以收敛。一般取个小值方便收敛,但是注意,一定不能取0
R \boldsymbol{R} R 和 Q \boldsymbol{Q} Q 的取值根据实际情况而定。
将式 ( 24 ) 带入 (24)带入 (24)带入(25),不难发现:如果更相信观测,那么增大 K k {\boldsymbol{K}}_k Kk ,也就是增大运动预测的噪声(过程噪声)的协方差 R \boldsymbol{R} R降低观测噪声协方差 Q \boldsymbol{Q} Q,如果更相信运动的预测,那么减小 K k {\boldsymbol{K}}_k Kk ,也就是增大 Q \boldsymbol{Q} Q降低 R \boldsymbol{R} R
这也是符合直觉的,如果相信某个传感器准确,也就是这个传感器数据的方差小。
(3)套公式
x
ˇ
k
=
A
k
x
^
k
−
1
+
u
k
,
P
ˇ
k
=
A
k
P
^
k
−
1
A
k
T
+
R
\check{\boldsymbol{x}}_k=\boldsymbol{A}_k \hat{\boldsymbol{x}}_{k-1}+\boldsymbol{u}_k, \quad \check{\boldsymbol{P}}_k=\boldsymbol{A}_k \hat{\boldsymbol{P}}_{k-1} \boldsymbol{A}_k^T+\boldsymbol{R}
xˇk=Akx^k−1+uk,Pˇk=AkP^k−1AkT+R
K
=
P
ˇ
k
C
k
T
(
C
k
P
ˇ
k
C
k
T
+
Q
)
−
1
\boldsymbol{K}=\check{\boldsymbol{P}}_k \boldsymbol{C}_k^T\left(\boldsymbol{C}_k \check{\boldsymbol{P}}_k \boldsymbol{C}_k^T+\boldsymbol{Q}\right)^{-1}
K=PˇkCkT(CkPˇkCkT+Q)−1
x
^
k
=
K
z
k
+
(
I
−
K
C
k
)
x
ˇ
k
=
x
ˇ
k
+
K
(
z
k
−
C
k
x
ˇ
k
)
\hat{\boldsymbol{x}}_k =\boldsymbol{K} \boldsymbol{z}_k+\left(\boldsymbol{I}-\boldsymbol{K} \boldsymbol{C}_k\right) \check{\boldsymbol{x}}_k =\check{\boldsymbol{x}}_k+\boldsymbol{K}\left(\boldsymbol{z}_k-\boldsymbol{C}_k \check{\boldsymbol{x}}_k\right)
x^k=Kzk+(I−KCk)xˇk=xˇk+K(zk−Ckxˇk)
P
^
k
=
(
I
−
K
C
k
)
P
ˇ
k
\hat{\boldsymbol{P}}_k=\left(\boldsymbol{I}-\boldsymbol{K} \boldsymbol{C}_k\right) \check{\boldsymbol{P}}_k
P^k=(I−KCk)Pˇk
非线性系统和拓展卡尔曼滤波器EKF
经典卡尔曼滤波器的局限性
经典卡尔曼滤波器只适用于线性系统(式 ( 12 ) (12) (12)),当时常见的系统(比如slam)都是非线性的
{ x k = f ( x ^ k − 1 , u k ) + w k z k = h ( x ^ k ) + v k k = 1 , … , N (8) \left\{\begin{array}{l} \boldsymbol{x}_k = f\left(\hat{\boldsymbol{x}}_{k-1}, \boldsymbol{u}_k\right)+\boldsymbol{w}_k \\ \boldsymbol{z}_k= h\left(\hat{\boldsymbol{x}}_{k}\right)+\boldsymbol{v}_k \end{array} \quad k=1, \ldots, N\right. \tag{8} {xk=f(x^k−1,uk)+wkzk=h(x^k)+vkk=1,…,N(8)
其中 f ( ) f() f()和 h ( ) h() h()都是非线性函数
通常的做法是,在某个点附近考虑运动方程以及观测方程的一阶泰勒展开,只保留一阶项,即线性的部分,然后按照线性系统进行推导。令 k − 1 k-1 k−1时刻的均值与协方差矩阵为 x ^ k − 1 \hat{\boldsymbol{x}}_{k-1} x^k−1, P ^ k − 1 \hat{\boldsymbol{P}}_{k-1} P^k−1。在 k k k 时刻,我们把运动方程和观测方程,在 x ^ k − 1 \hat{\boldsymbol{x}}_{k-1} x^k−1, P ^ k − 1 \hat{\boldsymbol{P}}_{k-1} P^k−1 处进行线性化(相当于一阶泰勒展开),有:
x k ≈ f ( x ^ k − 1 , u k ) + ∂ f ∂ x k − 1 ∣ x ^ k − 1 ( x k − 1 − x ^ k − 1 ) + w k \boldsymbol{x}_k \approx f\left(\hat{\boldsymbol{x}}_{k-1}, \boldsymbol{u}_k\right)+\left.\frac{\partial f}{\partial \boldsymbol{x}_{k-1}}\right|_{\hat{\boldsymbol{x}}_{k-1}}\left(\boldsymbol{x}_{k-1}-\hat{\boldsymbol{x}}_{k-1}\right)+\boldsymbol{w}_k xk≈f(x^k−1,uk)+∂xk−1∂f x^k−1(xk−1−x^k−1)+wk
记这里的偏导数为:
F
=
∂
f
∂
x
k
−
1
∣
x
^
k
−
1
\boldsymbol{F}=\left.\frac{\partial f}{\partial \boldsymbol{x}_{k-1}}\right|_{\hat{\boldsymbol{x}}_{k-1}}
F=∂xk−1∂f
x^k−1
z k ≈ h ( x ˇ k ) + ∂ h ∂ x k ∣ x ˇ k ( x k − x ^ k ) + n k \boldsymbol{z}_k \approx h\left(\check{\boldsymbol{x}}_k\right)+\left.\frac{\partial h}{\partial \boldsymbol{x}_k}\right|_{\check{\boldsymbol{x}}_k}\left(\boldsymbol{x}_k-\hat{\boldsymbol{x}}_k\right)+\boldsymbol{n}_k zk≈h(xˇk)+∂xk∂h xˇk(xk−x^k)+nk
记这里的偏导数为:
H
=
∂
h
∂
x
k
∣
x
ˇ
k
\boldsymbol{H}=\left.\frac{\partial h}{\partial \boldsymbol{x}_k}\right|_{\check{\boldsymbol{x}}_k}
H=∂xk∂h
xˇk
那么,在预测步骤中,根据运动方程有:
P ( x k ∣ x 0 , u 1 : k , z 0 : k − 1 ) = N ( f ( x ^ k − 1 , u k ) , F P ^ k − 1 F T + R k ) P\left(\boldsymbol{x}_k \mid \boldsymbol{x}_0, \boldsymbol{u}_{1: k}, \boldsymbol{z}_{0: k-1}\right)=N\left(f\left(\hat{\boldsymbol{x}}_{k-1}, \boldsymbol{u}_k\right), \boldsymbol{F} \hat{\boldsymbol{P}}_{k-1} \boldsymbol{F}^{\mathrm{T}}+\boldsymbol{R}_k\right) P(xk∣x0,u1:k,z0:k−1)=N(f(x^k−1,uk),FP^k−1FT+Rk)
这些推导和卡尔曼滤波是十分相似的。为方便表述,记这里先验和协方差的均值为
x ˇ k = f ( x ^ k − 1 , u k ) , P ˇ k = F P ^ k F T + R k \check{\boldsymbol{x}}_k=f\left(\hat{\boldsymbol{x}}_{k-1}, \boldsymbol{u}_k\right), \quad \check{\boldsymbol{P}}_k=\boldsymbol{F} \hat{\boldsymbol{P}}_k \boldsymbol{F}^T+\boldsymbol{R}_k xˇk=f(x^k−1,uk),Pˇk=FP^kFT+Rk
然后,考虑在观测中,我们有:
P ( z k ∣ x k ) = N ( h ( x ˇ k ) + H ( x k − x ˇ k ) , Q k ) P\left(\boldsymbol{z}_k \mid \boldsymbol{x}_k\right)=N\left(h\left(\check{\boldsymbol{x}}_k\right)+\boldsymbol{H}\left(\boldsymbol{x}_k-\check{\boldsymbol{x}}_k\right), \boldsymbol{Q}_k\right) P(zk∣xk)=N(h(xˇk)+H(xk−xˇk),Qk)
最后,根据最开始的 Bayes 展开式,可以推导出 xk 的后验概率形式。我们略去中间的推导过程,只介绍其结果。读者可以仿照着卡尔曼滤波器的方式,推导 EKF 的预测与更新方程。简而言之,我们会先定义一个卡尔曼增益:
K
k
=
P
ˇ
k
H
T
(
H
P
ˇ
k
H
T
+
Q
k
)
−
1
\boldsymbol{K}_k=\check{\boldsymbol{P}}_k \boldsymbol{H}^{\mathrm{T}}\left(\boldsymbol{H} \check{\boldsymbol{P}}_k \boldsymbol{H}^{\mathrm{T}}+\boldsymbol{Q}_k\right)^{-1}
Kk=PˇkHT(HPˇkHT+Qk)−1
在卡尔曼增益的基础上,后验概率的形式为:
x
^
k
=
x
ˇ
k
+
K
k
(
z
k
−
h
(
x
ˇ
k
)
)
,
P
^
k
=
(
I
−
K
k
H
)
P
ˇ
k
\hat{\boldsymbol{x}}_k=\check{\boldsymbol{x}}_k+\boldsymbol{K}_k\left(\boldsymbol{z}_k-h\left(\check{\boldsymbol{x}}_k\right)\right), \hat{\boldsymbol{P}}_k=\left(\boldsymbol{I}-\boldsymbol{K}_k \boldsymbol{H}\right) \check{\boldsymbol{P}}_k
x^k=xˇk+Kk(zk−h(xˇk)),P^k=(I−KkH)Pˇk
卡尔曼滤波器给出了在线性化之后,状态变量分布的变化过程。在线性系统和高斯噪声下,卡尔曼滤波器给出了无偏最优估计。而在 SLAM 这种非线性的情况下,它给出了单次线性近似下最大后验估计(MAP)
实际应用1 用陀螺仪数据估计角度
未完待续……敬请关注
实际应用2 融合陀螺仪数据和加速度计估计角度
未完待续……敬请关注