扩展卡尔曼滤波器
前言
卡尔曼滤波器作为状态估计的一种重要方法,有着其使用范围的局限性。即,卡尔曼滤波器仅适用于线性系统中。
这是因为卡尔曼滤波器使用的前提是先验分布是高斯似然分布的共轭先验。这就要求先验分布也是高斯的,并且可以确保后验分布与先验分布是属于同一个分布族的,即同属于高斯分布。
根据高斯分布的性质,高斯分布在经过一个线性系统后,仍为高斯分布。但是高斯分布在经过一个非线性系统后,则无法确保其仍为高斯分布。直接使用线性卡尔曼滤波器的公式便不再可行。
因此,需要找到一种方法,使得在非线性系统下仍可以使用卡尔曼滤波器公式。
一、泰勒展开
扩展卡尔曼滤波器的主要思想是对非线性系统做线性化处理,转变为线性系统,并在此基础上实现预测和滤波。其使用的线性化处理方式是泰勒展开。
先回顾以下泰勒展开的公式。假设有一函数为
f
(
x
)
f(x)
f(x),在
x
0
x_{0}
x0处做泰勒展开后,可以得到
f
(
x
)
=
f
(
x
0
)
+
∂
f
∂
x
(
x
−
x
0
)
+
.
.
.
+
o
(
x
−
x
0
)
f(x) = f(x_{0})+\frac{\partial f}{\partial x}(x-x_{0}) + ... + o(x-x_{0})
f(x)=f(x0)+∂x∂f(x−x0)+...+o(x−x0)
其中,
o
(
x
−
x
0
)
o(x-x_{0})
o(x−x0)为高阶项。
对于多元函数,泰勒展开公式类似与一元函数的公式。假设有一个多元函数
f
(
x
,
y
)
f(x,y)
f(x,y),在
x
0
x_{0}
x0,
y
0
y_{0}
y0做泰勒展开后,可以得到
f
(
x
,
y
)
=
f
(
x
0
,
y
0
)
+
∂
f
∂
x
(
x
−
x
0
)
+
∂
f
∂
x
(
y
−
y
0
)
+
.
.
.
+
o
(
x
−
x
0
)
+
o
(
y
−
y
0
)
f(x,y) = f(x_{0},y_{0})+\frac{\partial f}{\partial x}(x-x_{0})+\frac{\partial f}{\partial x}(y-y_{0}) + ... + o(x-x_{0}) +o(y-y_{0})
f(x,y)=f(x0,y0)+∂x∂f(x−x0)+∂x∂f(y−y0)+...+o(x−x0)+o(y−y0)
二、扩展卡尔曼滤波器
1.推导
首先处理状态转移函数的线性化。
假设
k
+
1
k+1
k+1时刻的状态向量表示为
x
k
+
1
=
f
(
x
k
,
u
k
,
v
k
)
x_{k+1} = f(x_{k},u_{k},v_{k})
xk+1=f(xk,uk,vk)
其中,
f
(
⋅
)
f(\cdot)
f(⋅)为状态转移方程,
x
k
x_{k}
xk为
k
k
k时刻的状态向量,
u
k
u_{k}
uk为k时刻的机动项,
v
k
v_{k}
vk为k时刻的过程噪声。假设过程噪声服从均值为0的高斯分布,
v
k
∼
N
(
0
,
Q
)
v_{k}\sim N(0,Q)
vk∼N(0,Q)。由于在实际情况下,我们无法获取
k
k
k时刻的状态向量,因此使用
k
k
k时刻的状态向量估计值
x
^
k
∣
k
\hat x_{k|k}
x^k∣k代替。同时,为了推导方便,暂时不考虑机动项
u
k
u_{k}
uk,并令过程噪声
v
k
v_{k}
vk为0。由此得到下式
x
k
+
1
=
f
(
x
k
,
v
k
)
=
f
(
x
^
k
∣
k
,
0
)
+
∂
f
∂
x
k
∣
x
k
=
x
^
k
∣
k
(
x
k
−
x
^
k
∣
k
)
+
∂
f
∂
v
k
∣
v
k
=
0
(
v
k
−
0
)
+
.
.
.
+
o
(
x
k
−
x
^
k
∣
k
)
+
o
(
v
k
−
0
)
\begin{aligned} x_{k+1}&= f(x_{k},v_{k}) \\&= f(\hat x_{k|k},0)+\frac{\partial f}{\partial x_{k}|_{x_{k}=\hat x_{k|k}}}(x_{k} - \hat x_{k|k}) + \frac{\partial f}{\partial v_{k}|_{v_{k}=0}}(v_{k} - 0)+...+ o(x_{k} - \hat x_{k|k}) +o(v_{k} - 0) \end{aligned}
xk+1=f(xk,vk)=f(x^k∣k,0)+∂xk∣xk=x^k∣k∂f(xk−x^k∣k)+∂vk∣vk=0∂f(vk−0)+...+o(xk−x^k∣k)+o(vk−0)
令
∂
f
∂
x
k
∣
x
k
=
x
^
k
∣
k
=
F
\frac{\partial f}{\partial x_{k}|_{x_{k}=\hat x_{k|k}}} = F
∂xk∣xk=x^k∣k∂f=F,
∂
f
∂
v
k
∣
v
k
=
0
=
V
\frac{\partial f}{\partial v_{k}|_{v_{k}=0}} = V
∂vk∣vk=0∂f=V,并忽略高阶项,仅保留一阶项,则上式可简写为
x
k
+
1
≈
f
(
x
^
k
∣
k
,
0
)
+
F
(
x
k
−
x
^
k
∣
k
)
+
V
v
k
x_{k+1}\approx f(\hat x_{k|k},0) + F(x_{k} - \hat x_{k|k})+Vv_{k}
xk+1≈f(x^k∣k,0)+F(xk−x^k∣k)+Vvk
由上式可推导出预测状态向量的表达式
x
^
k
+
1
∣
k
≈
E
{
f
(
x
^
k
∣
k
,
0
)
+
F
(
x
k
−
x
^
k
∣
k
)
+
V
v
k
}
=
f
(
x
^
k
∣
k
,
0
)
\hat x_{k+1|k}\approx E\left \{ f(\hat x_{k|k},0) + F(x_{k} - \hat x_{k|k})+Vv_{k} \right \} = f(\hat x_{k|k},0)
x^k+1∣k≈E{f(x^k∣k,0)+F(xk−x^k∣k)+Vvk}=f(x^k∣k,0)
预测误差协方差矩阵的表达式为
P
k
+
1
∣
k
=
E
{
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
T
}
=
E
{
(
F
(
x
k
−
x
^
k
∣
k
)
+
V
v
k
)
(
F
(
x
k
−
x
^
k
∣
k
)
+
V
v
k
)
T
}
=
F
P
k
∣
k
F
T
+
V
Q
V
T
\begin{aligned} P_{k+1|k} &= E\left \{ (x_{k+1}-\hat x_{k+1|k})(x_{k+1}-\hat x_{k+1|k})^{T}\right \} \\ &= E\left \{ (F(x_{k} - \hat x_{k|k})+Vv_{k})(F(x_{k} - \hat x_{k|k})+Vv_{k})^{T}\right \} \\ &= FP_{k|k}F^{T} + VQV^{T} \end{aligned}
Pk+1∣k=E{(xk+1−x^k+1∣k)(xk+1−x^k+1∣k)T}=E{(F(xk−x^k∣k)+Vvk)(F(xk−x^k∣k)+Vvk)T}=FPk∣kFT+VQVT
同理,接下来再对观测函数做线性化。
假设
k
+
1
k+1
k+1时刻的观测向量表示为
z
k
+
1
=
h
(
x
k
+
1
,
w
k
+
1
)
z_{k+1} = h(x_{k+1},w_{k+1})
zk+1=h(xk+1,wk+1)
其中,
h
(
⋅
)
h(\cdot)
h(⋅)为观测函数,
x
k
x_{k}
xk为
k
+
1
k+1
k+1时刻的状态向量,
w
k
w_{k}
wk为k+1时刻的观测噪声。假设观测噪声服从均值为0的高斯分布,
w
k
∼
N
(
0
,
R
)
w_{k}\sim N(0,R)
wk∼N(0,R)。使用
k
+
1
k+1
k+1时刻的状态向量预测值
x
^
k
+
1
∣
k
\hat x_{k+1|k}
x^k+1∣k代替
x
k
+
1
x_{k+1}
xk+1,同时令
w
k
w_{k}
wk为0。由此得到下式
z
k
+
1
=
h
(
x
k
+
1
,
w
k
+
1
)
=
h
(
x
^
k
+
1
∣
k
,
0
)
+
∂
h
∂
x
k
+
1
∣
x
k
+
1
=
x
^
k
+
1
∣
k
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
+
∂
h
∂
w
k
+
1
∣
w
k
+
1
=
0
(
w
k
+
1
−
0
)
+
.
.
.
+
o
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
+
o
(
w
k
+
1
−
0
)
\begin{aligned} z_{k+1} &= h(x_{k+1},w_{k+1}) \\ &= h(\hat x_{k+1|k},0) + \frac{\partial h}{\partial x_{k+1}|_{x_{k+1}=\hat x_{k+1|k}}}(x_{k+1}-\hat x_{k+1|k})+\frac{\partial h}{\partial w_{k+1}|_{w_{k+1}=0}}(w_{k+1}-0) +...+o(x_{k+1}-\hat x_{k+1|k})+o(w_{k+1}-0) \end{aligned}
zk+1=h(xk+1,wk+1)=h(x^k+1∣k,0)+∂xk+1∣xk+1=x^k+1∣k∂h(xk+1−x^k+1∣k)+∂wk+1∣wk+1=0∂h(wk+1−0)+...+o(xk+1−x^k+1∣k)+o(wk+1−0)
令
∂
h
∂
x
k
+
1
∣
x
k
+
1
=
x
^
k
+
1
∣
k
=
H
\frac{\partial h}{\partial x_{k+1}|_{x_{k+1}=\hat x_{k+1|k}}} = H
∂xk+1∣xk+1=x^k+1∣k∂h=H,
∂
f
∂
w
k
∣
w
k
=
0
=
W
\frac{\partial f}{\partial w_{k}|_{w_{k}=0}} = W
∂wk∣wk=0∂f=W,并忽略高阶项,仅保留一阶项,则上式可简写为
z
k
+
1
≈
h
(
x
^
k
+
1
∣
k
,
0
)
+
H
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
+
W
w
k
+
1
z_{k+1}\approx h(\hat x_{k+1|k},0) + H(x_{k+1} - \hat x_{k+1|k})+Ww_{k+1}
zk+1≈h(x^k+1∣k,0)+H(xk+1−x^k+1∣k)+Wwk+1
由上式可推导出预测观测向量的表达式
z
^
k
+
1
≈
E
{
h
(
x
^
k
+
1
∣
k
,
0
)
+
H
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
+
W
w
k
+
1
}
=
h
(
x
^
k
+
1
∣
k
,
0
)
\hat z_{k+1}\approx E\left \{ h(\hat x_{k+1|k},0) + H(x_{k+1} - \hat x_{k+1|k})+Ww_{k+1} \right \} = h(\hat x_{k+1|k},0)
z^k+1≈E{h(x^k+1∣k,0)+H(xk+1−x^k+1∣k)+Wwk+1}=h(x^k+1∣k,0)
则
P
x
z
,
k
+
1
∣
k
=
E
{
(
x
k
+
1
−
x
^
k
+
1
)
(
z
k
+
1
−
z
^
k
+
1
∣
k
)
T
}
=
E
{
(
x
^
k
+
1
∣
k
−
z
k
+
1
)
(
H
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
+
W
w
k
+
1
)
T
}
=
P
k
+
1
∣
k
H
\begin{aligned} P_{xz,k+1|k} &= E\left \{ ( x_{k+1} - \hat x_{k+1})(z_{k+1} - \hat z_{k+1|k})^{T}\right \} \\ &= E\left \{ (\hat x_{k+1|k} - z_{k+1})(H(x_{k+1} - \hat x_{k+1|k})+Ww_{k+1})^{T}\right \} \\ &= P_{k+1|k}H \end{aligned}
Pxz,k+1∣k=E{(xk+1−x^k+1)(zk+1−z^k+1∣k)T}=E{(x^k+1∣k−zk+1)(H(xk+1−x^k+1∣k)+Wwk+1)T}=Pk+1∣kH
P
z
z
,
k
+
1
∣
k
=
E
{
(
z
k
+
1
−
z
^
k
+
1
∣
k
)
(
z
k
+
1
−
z
^
k
+
1
∣
k
)
T
}
=
E
{
(
H
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
+
W
w
k
+
1
)
(
H
(
x
k
+
1
−
x
^
k
+
1
∣
k
)
+
W
w
k
+
1
)
T
}
=
H
P
k
+
1
∣
k
H
T
+
W
R
k
+
1
W
T
\begin{aligned} P_{zz,k+1|k} &= E\left \{ (z_{k+1} - \hat z_{k+1|k})(z_{k+1} - \hat z_{k+1|k})^{T}\right \} \\ &= E\left \{ (H(x_{k+1} - \hat x_{k+1|k})+Ww_{k+1})(H(x_{k+1} - \hat x_{k+1|k})+Ww_{k+1})^{T}\right \} \\ &= HP_{k+1|k}H^{T} + WR_{k+1}W^{T} \end{aligned}
Pzz,k+1∣k=E{(zk+1−z^k+1∣k)(zk+1−z^k+1∣k)T}=E{(H(xk+1−x^k+1∣k)+Wwk+1)(H(xk+1−x^k+1∣k)+Wwk+1)T}=HPk+1∣kHT+WRk+1WT
可以得到卡尔曼增益为
K
=
P
x
z
P
z
z
−
1
=
P
k
+
1
∣
k
H
(
H
P
k
+
1
∣
k
H
T
+
W
R
k
+
1
W
T
)
−
1
K = P_{xz}P_{zz}^{-1} = P_{k+1|k}H(HP_{k+1|k}H^{T} + WR_{k+1}W^{T})^{-1}
K=PxzPzz−1=Pk+1∣kH(HPk+1∣kHT+WRk+1WT)−1
最后通过滤波步骤,得到后验的状态空间均值与方差
x
^
k
+
1
∣
k
+
1
=
x
^
k
+
1
∣
k
+
K
(
z
k
+
1
−
z
^
k
+
1
∣
k
)
\hat x_{k+1|k+1} = \hat x_{k+1|k} + K(z_{k+1} - \hat z_{k+1|k})
x^k+1∣k+1=x^k+1∣k+K(zk+1−z^k+1∣k)
P
k
+
1
∣
k
+
1
=
(
I
−
K
H
)
P
k
+
1
∣
k
(
I
−
K
H
)
T
+
K
R
k
+
1
K
T
P_{k+1|k+1}=(I-KH)P_{k+1|k}(I-KH)^{T}+KR_{k+1}K^{T}
Pk+1∣k+1=(I−KH)Pk+1∣k(I−KH)T+KRk+1KT
2.总结
扩展卡尔曼滤波器也是由五个公式构成的
预测步:
x
^
k
+
1
∣
k
=
f
(
x
^
k
∣
k
,
0
)
\hat x_{k+1|k}= f(\hat x_{k|k},0)
x^k+1∣k=f(x^k∣k,0)
P
k
+
1
∣
k
=
F
P
k
∣
k
F
T
+
V
Q
V
T
P_{k+1|k} = FP_{k|k}F^{T} + VQV^{T}
Pk+1∣k=FPk∣kFT+VQVT
卡尔曼增益计算:
K = P x z P z z − 1 = P k + 1 ∣ k H ( H P k + 1 ∣ k H T + W R k + 1 W T ) − 1 K = P_{xz}P_{zz}^{-1} = P_{k+1|k}H(HP_{k+1|k}H^{T} + WR_{k+1}W^{T})^{-1} K=PxzPzz−1=Pk+1∣kH(HPk+1∣kHT+WRk+1WT)−1
滤波步:
x
^
k
+
1
∣
k
+
1
=
x
^
k
+
1
∣
k
+
K
(
z
k
+
1
−
z
^
k
+
1
∣
k
)
\hat x_{k+1|k+1} = \hat x_{k+1|k} + K(z_{k+1} - \hat z_{k+1|k})
x^k+1∣k+1=x^k+1∣k+K(zk+1−z^k+1∣k)
P
k
+
1
∣
k
+
1
=
(
I
−
K
H
)
P
k
+
1
∣
k
(
I
−
K
H
)
T
+
K
R
k
+
1
K
T
P_{k+1|k+1}=(I-KH)P_{k+1|k}(I-KH)^{T}+KR_{k+1}K^{T}
Pk+1∣k+1=(I−KH)Pk+1∣k(I−KH)T+KRk+1KT