1.卡尔曼滤波是什么?
为什么要对事物状态进行估计?由于我们无法准确知道物体的当前状态,为了获得事物状态我们需要测量。但是测量值并不是准确的,总会存在噪声。卡尔曼滤波是一种结合预测(先验分布)和测量更新(似然)的状态估计算法。
预测模块就是对物体的运动建立运动模型,本文栗子中的行人状态估计我们采用恒速度模型(CV),通过对上一时刻的最优估计进行运动模型转换,得到当前时刻的估计状态,以及预测误差。
更新模块是根据当前时刻的测量值和预测状态得到当前时刻的状态最优估计。
2.卡尔曼公式
Predict:
X
k
=
F
X
^
k
−
1
+
B
U
k
+
w
k
{{X}_{k}}=F{{\overset{\hat{\ }}{\mathop{X}}\,}_{k-1}}+B{{U}_{k}}+{{w}_{k}}
Xk=FX ^k−1+BUk+wk
P
k
=
F
P
k
−
1
F
+
Q
k
{{P}_{k}}=F{{P}_{k-1}}F+{{Q}_{k}}
Pk=FPk−1F+Qk
其中
X
k
{{X}_{k}}
Xk为k时刻的状态预测值,F为状态转移矩阵,
X
^
k
−
1
{{\overset{\hat{\ }}{\mathop{X}}\,}_{k-1}}
X ^k−1为k-1时刻的最优估计。
B
U
k
B{{U}_{k}}
BUk为外部输入,
w
k
{{w}_{k}}
wk为过程激励噪声。P为预测误差,
P
k
{{P}_{k}}
Pk表示为k次先验估计协方差矩阵,
P
k
−
1
{{P}_{k-1}}
Pk−1为k-1次后验估计协方差矩阵。
Q
k
{{Q}_{k}}
Qk为过程激励噪声协方差矩阵。
Update:
K
K
=
P
K
H
T
(
H
P
K
H
T
+
R
)
−
1
{{K}_{K}}={{P}_{K}}{{H}^{T}}{{(H{{P}_{K}}{{H}^{T}}+R)}^{-1}}
KK=PKHT(HPKHT+R)−1
X
K
^
=
X
K
+
K
K
(
Z
K
−
H
X
K
)
\overset{\hat{\ }}{\mathop{{{X}_{K}}}}\,={{X}_{K}}+{{K}_{K}}({{Z}_{K}}-H{{X}_{K}})
XK ^=XK+KK(ZK−HXK)
P
K
=
(
1
−
K
K
H
)
P
K
{{P}_{K}}=(1-{{K}_{K}}H){{P}_{K}}
PK=(1−KKH)PK
其中K为卡尔曼增益,H为观测矩阵,R为测量噪声协方差矩阵,
Z
K
{{Z}_{K}}
ZK为k时刻测量值。
3.以行人状态估计举个栗子
①predict:
先对行人状态建立一个运动模型,我们在此选择恒速度运动模型,
X
=
(
p
x
,
p
y
,
v
x
,
v
y
)
T
X={{({{p}_{x}},{{p}_{y}},{{v}_{x}},{{v}_{y}})}^{T}}
X=(px,py,vx,vy)T,此时
B
U
k
B{{U}_{k}}
BUk为0,过程噪声
w
k
{{w}_{k}}
wk其实是行人加减速引起的。公式如下:
X
K
=
[
1
0
△
t
0
0
1
0
△
t
0
0
1
0
0
0
0
1
]
[
p
x
p
y
v
x
v
y
]
k
−
1
+
[
1
2
a
x
△
t
2
1
2
a
y
△
t
2
a
x
△
t
a
y
△
t
]
k
−
1
{{X}_{K}}=\left[ \begin{matrix} 1 & 0 & \vartriangle t & 0 \\ 0 & 1 & 0 & \vartriangle t \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right]{{\left[ \begin{matrix} {{p}_{x}} \\ {{p}_{y}} \\ {{v}_{x}} \\ {{v}_{y}} \\ \end{matrix} \right]}_{k-1}}+{{\left[ \begin{matrix} \tfrac{1}{2}{{a}_{x}}\vartriangle {{t}^{2}} \\ \tfrac{1}{2}{{a}_{y}}\vartriangle {{t}^{2}} \\ {{a}_{x}}\vartriangle t \\ {{a}_{y}}\vartriangle t \\ \end{matrix} \right]}_{k-1}}
XK=⎣⎢⎢⎡10000100△t0100△t01⎦⎥⎥⎤⎣⎢⎢⎡pxpyvxvy⎦⎥⎥⎤k−1+⎣⎢⎢⎡21ax△t221ay△t2ax△tay△t⎦⎥⎥⎤k−1
P
k
=
F
P
k
−
1
F
+
Q
k
{{P}_{k}}=F{{P}_{k-1}}F+{{Q}_{k}}
Pk=FPk−1F+Qk
其中P第一次应该初始化为???由于w(过程激励噪声)是随机代入的,所以为高斯白噪声,满足
p
(
w
)
~
N
(
0
,
Q
)
p(w)\tilde{\ }N(0,Q)
p(w) ~N(0,Q),Q即过程噪声协方差矩阵,表示如下:
Q
=
[
σ
p
x
2
σ
p
x
p
y
σ
p
x
v
x
σ
p
x
v
y
σ
p
y
p
x
σ
p
y
2
σ
p
y
v
x
σ
p
y
v
y
σ
v
x
p
x
σ
v
x
p
y
σ
v
x
2
σ
v
x
v
y
σ
v
y
p
x
σ
v
y
p
y
σ
v
y
v
x
σ
v
y
2
]
Q=\left[ \begin{matrix} \sigma _{{{p}_{x}}}^{2} & {{\sigma }_{{{p}_{x}}{{p}_{y}}}} & {{\sigma }_{{{p}_{x}}{{v}_{x}}}} & {{\sigma }_{{{p}_{x}}{{v}_{y}}}} \\ {{\sigma }_{{{p}_{y}}{{p}_{x}}}} & \sigma _{{{p}_{y}}}^{2} & {{\sigma }_{{{p}_{y}}{{v}_{x}}}} & {{\sigma }_{{{p}_{y}}{{v}_{y}}}} \\ {{\sigma }_{{{v}_{x}}{{p}_{x}}}} & {{\sigma }_{{{v}_{x}}{{p}_{y}}}} & \sigma _{{{v}_{x}}}^{2} & {{\sigma }_{{{v}_{x}}{{v}_{y}}}} \\ {{\sigma }_{{{v}_{y}}{{p}_{x}}}} & {{\sigma }_{{{v}_{y}}{{p}_{y}}}} & {{\sigma }_{{{v}_{y}}{{v}_{x}}}} & \sigma _{{{v}_{y}}}^{2} \\ \end{matrix} \right]
Q=⎣⎢⎢⎡σpx2σpypxσvxpxσvypxσpxpyσpy2σvxpyσvypyσpxvxσpyvxσvx2σvyvxσpxvyσpyvyσvxvyσvy2⎦⎥⎥⎤
Q表示的是噪声
[
1
2
a
x
△
t
2
1
2
a
y
△
t
2
a
x
△
t
a
y
△
t
]
T
{{\left[ \begin{matrix} \tfrac{1}{2}{{a}_{x}}\vartriangle {{t}^{2}} & \tfrac{1}{2}{{a}_{y}}\vartriangle {{t}^{2}} & {{a}_{x}}\vartriangle t & {{a}_{y}}\vartriangle t \\ \end{matrix} \right]}^{T}}
[21ax△t221ay△t2ax△tay△t]T各个变量之间的协方差,假设
p
(
a
x
)
~
N
(
0
,
a
)
p(ax)\tilde{\ }N(0,a)
p(ax) ~N(0,a),
p
(
a
y
)
~
N
(
0
,
b
)
p(ay)\tilde{\ }N(0,b)
p(ay) ~N(0,b),则此时Q表示如下:
Q
=
[
1
4
a
x
△
t
4
0
1
2
a
x
△
t
3
0
0
1
4
a
y
△
t
4
0
1
2
a
y
△
t
3
1
2
a
x
△
t
3
0
a
x
△
t
2
0
0
1
2
a
y
△
t
3
0
a
y
△
t
2
]
Q=\left[ \begin{matrix} \tfrac{1}{\text{4}}{{a}_{x}}\vartriangle {{t}^{\text{4}}} & \text{0} & \tfrac{1}{\text{2}}{{a}_{x}}\vartriangle {{t}^{\text{3}}} & \text{0} \\ \text{0} & \tfrac{1}{\text{4}}{{a}_{y}}\vartriangle {{t}^{\text{4}}} & 0 & \tfrac{1}{\text{2}}{{a}_{y}}\vartriangle {{t}^{\text{3}}} \\ \tfrac{1}{\text{2}}{{a}_{x}}\vartriangle {{t}^{\text{3}}} & 0 & {{a}_{x}}\vartriangle {{t}^{\text{2}}} & 0 \\ 0 & \tfrac{1}{\text{2}}{{a}_{y}}\vartriangle {{t}^{\text{3}}} & 0 & {{a}_{y}}\vartriangle {{t}^{\text{2}}} \\ \end{matrix} \right]
Q=⎣⎢⎢⎡41ax△t4021ax△t30041ay△t4021ay△t321ax△t30ax△t20021ay△t30ay△t2⎦⎥⎥⎤
②update
假设此时测量值来自于里程计,即测量设备是底层光电编码器,则
v
x
{{v}_{x}}
vx
v
y
{{v}_{y}}
vy是可观测的,观测矩阵H表示如下:
H
=
[
0
0
1
0
0
0
0
1
]
H=\left[ \begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right]
H=[00001001]
注意:状态观测矩阵,与建立的运动模型和传感器都有关系,在该例子中CV模型下,如果传感器是里程计,则观测值为速度;如果传感器是激光雷达则观测值应该为位置xy,如果是毫米波雷达的话,由于毫米波雷达的输出结果是极坐标下的结果不是x,y,vx,vy,转换公式为非线性的,所以观测矩阵应该计算其雅可比矩阵。
测量噪声R表示如下:
R
=
[
σ
V
X
2
0
0
σ
V
Y
2
]
R=\left[ \begin{matrix} \sigma _{{{V}_{X}}}^{2} & 0 \\ 0 & \sigma _{{{V}_{Y}}}^{2} \\ \end{matrix} \right]
R=[σVX200σVY2]
测量噪声是测量传感器的性能指标,一般可由生产厂家提供。
得到卡尔曼增益后,即可得到当前状态的最有估计。并更新误差协方差矩阵。
到此,卡尔曼滤波过程就结束了。
本文参考大神博客:https://blog.csdn.net/AdamShan/article/details/78248421
github代码推荐:https://github.com/JunshengFu/tracking-with-Extended-Kalman-Filter