IMU驱动系统的误差状态运动学
误差状态卡尔曼滤波器简介
在误差状态滤波器的公式中,我们谈论真实状态、名义状态和误差状态值,其中真实状态被表示为名义状态和误差状态的适当组合(线性求和、四元数乘积或矩阵乘积)。这样做的想法是将名义状态视为大信号(可在非线性方式下积分),将误差状态视为小信号(因此可线性积分并适用于线性高斯滤波)。
误差状态滤波器可以这样解释。一方面,高频IMU数据 u m \textbf{u}_m um被集成到名义状态 x \textbf{x} x中。该名义状态不考虑噪声项 w \textbf{w} w和其他可能的模型缺陷。因此,它将累积误差。这些误差被收集到误差状态 δ x δ\textbf{x} δx中,并且通过误差状态卡尔曼滤波器(ESKF)进行估计,这一次包括所有的噪声和扰动。误差状态由小信号量组成,其演变函数由正确定义的(时变)线性动态系统确定,其动态、控制和测量矩阵从名义状态的值计算而来。与名义状态的集成并行,ESKF预测误差状态的高斯估计。它只进行预测,因为现在还没有其他可用的测量来纠正这些估计。滤波器校正是在到达IMU以外的信息(例如GPS、视觉等)时进行的,这些信息能够使误差可观测,并且通常比集成阶段低得多。这种校正提供了误差状态的后验高斯估计。此后,误差状态的平均值被注入到名义状态中,然后重置为零。误差状态的协方差矩阵被方便地更新以反映此重置。系统会一直这样运行下去。
连续时间的系统运动学
误差状态卡尔曼滤波中的所有变量:
误差状态的幅值 | 真实值 | 名义值 | 误差 | 组合 | 测量值 | 噪声 |
---|---|---|---|---|---|---|
状态变量集合 | x t \textbf{x}_t xt | x \textbf{x} x | δ x δ\textbf{x} δx | x t = x ⊕ δ x \mathbf x_t=\mathbf x\oplus\delta\mathbf x xt=x⊕δx | ||
位置 | p t \mathbf{p}_t pt | p \mathbf{p} p | δ p \delta\mathbf{p} δp | p t = p + δ p \mathbf{p}_t=\mathbf{p}+\delta\mathbf{p} pt=p+δp | ||
速度 | v t \mathbf{v}_t vt | v \mathbf{v} v | δ v \delta\mathbf{v} δv | v t = v + δ v \mathbf{v}_t=\mathbf{v}+\delta\mathbf{v} vt=v+δv | ||
四元数 | q t \mathbf{q}_t qt | q \mathbf{q} q | δ q \delta\mathbf{q} δq | q t = q ⊗ δ q \mathbf q_{t}=\mathbf q\otimes\delta\mathbf q qt=q⊗δq | ||
旋转矩阵 | R t \mathbf{R}_t Rt | R \mathbf{R} R | δ R \delta\mathbf{R} δR | R t = R δ R \mathbf{R}_{t}=\mathbf{R}\delta\mathbf{R} Rt=RδR | ||
角度矢量 | δ θ \delta\boldsymbol{\theta} δθ | δ q = e δ θ / 2 δ R = e [ δ θ ] × \begin{array}{c}\delta\mathbf{q}=e^{\delta\boldsymbol{\theta}/2}\\ \delta\mathbf{R}=e^{[\delta\boldsymbol{\theta}]_\times}\end{array} δq=eδθ/2δR=e[δθ]× | ||||
加速度计偏置 | a b t \mathbf{a}_{bt} abt | a b \mathbf{a}_{b} ab | δ a b δ\textbf{a}_{b} δab | a b t = a b + δ a b \mathbf{a}_{b t}=\mathbf{a}_{b}+\delta\mathbf{a}_{b} abt=ab+δab | a w \mathbf{a}_{w} aw | |
陀螺仪偏置 | ω b t \omega_{bt} ωbt | ω b \omega_{b} ωb | δ ω b δ\omega_{b} δωb | ω b t = ω b + δ ω b \omega_{b t}=\omega_{b}+\delta \omega_{b} ωbt=ωb+δωb | ω w \omega_{w} ωw | |
重量矢量 | g t \mathbf{g}_{t} gt | g \mathbf{g} g | δ g δ\textbf{g} δg | g t = g + δ g \mathbf{g}_{ t}=\mathbf{g}+\delta\mathbf{g} gt=g+δg | ||
加速度 | a t \mathbf{a}_{t} at | a m \mathbf{a}_{m} am | a n \mathbf{a}_{n} an | |||
角速度 | ω t \omega_{t} ωt | ω m \omega_{m} ωm | ω n \omega_{n} ωn |
真实状态运动学
真实值运动学方程:
p
˙
t
=
v
t
v
˙
t
=
a
t
q
˙
t
=
1
2
q
t
⊗
ω
t
a
˙
b
t
=
a
w
ω
˙
b
t
=
ω
w
g
˙
t
=
0
\begin{array}{l}\dot{\mathbf{p}}_t=\mathbf{v}_t\\ \dot{\mathbf{v}}_t=\mathbf{a}_t\\ \dot{\mathbf{q}}_t=\frac{1}{2}q_t\otimes\mathbf{\omega}_t\\ \dot{\mathbf{a}}_{bt}=\mathbf{a}_w\\ \dot{\mathbf{\omega}}_{bt}=\mathbf{\omega}_w\\ \dot{\mathbf{g}}_t=0\end{array}
p˙t=vtv˙t=atq˙t=21qt⊗ωta˙bt=awω˙bt=ωwg˙t=0
机身坐标系下的加速度传感器读数
a
m
\mathbf{a}_{m}
am和角速度传感器读数
w
m
w_{m}
wm被用来获取机器人在当前时刻的真实加速度
a
t
\mathbf{a}_{t}
at和角速率
w
t
w_{t}
wt。由于传感器本身存在噪声,因此
a
m
\mathbf{a}_{m}
am和
w
m
w_{m}
wm是带噪声的传感器读数。
a
m
=
R
t
⊤
(
a
t
−
g
t
)
+
a
b
t
+
a
n
ω
m
=
ω
t
+
ω
b
t
+
ω
n
\begin{array}{l}\mathbf{a}_m=\mathbf{R}_t^\top(\mathbf{a}_t-\mathbf{g}_t)+\mathbf{a}_{bt}+\mathbf{a}_n\\ \mathbf{\omega}_m=\mathbf{\omega}_t+\mathbf{\omega}_{bt}+\mathbf{\omega}_n\end{array}
am=Rt⊤(at−gt)+abt+anωm=ωt+ωbt+ωn
通过
R
t
≜
R
{
q
t
}
\mathbf{R}_{t}\triangleq\mathbf{R}\{\mathbf{q}_{t}\}
Rt≜R{qt}分离出真实值。
a
t
=
R
t
(
a
m
−
a
b
t
−
a
n
)
+
g
t
ω
t
=
ω
m
−
ω
b
t
−
ω
n
\begin{array}{l}\mathbf{a}_t=\mathbf{R}_t(\mathbf{a}_m-\mathbf{a}_{bt}-\mathbf{a}_n)+\mathbf{g}_t\\ \mathbf{\omega}_t=\mathbf{\omega}_m-\mathbf{\omega}_{bt}-\mathbf{\omega}_n\end{array}
at=Rt(am−abt−an)+gtωt=ωm−ωbt−ωn
代入
p
˙
t
=
v
t
v
t
=
R
t
(
a
m
−
a
b
t
−
a
n
)
+
g
t
q
˙
t
=
1
2
q
t
⊗
(
ω
m
−
ω
b
t
−
ω
n
)
a
˙
b
t
=
a
w
ω
˙
t
=
ω
w
g
˙
t
=
0
\begin{array}{l}\dot{\textbf{p}}_t=\textbf{v}_t\\ \textbf{v}_t=\textbf{R}_t(\textbf{a}_m-\textbf{a}_{bt}-\textbf{a}_n)+\textbf{g}_t\\ \dot{\textbf{q}}_t=\frac{1}{2}\textbf{q}_t\otimes(\omega_m-\omega_{bt}-\omega_n)\\ \dot{\textbf{a}}_{bt}=\textbf{a}_w\\ \dot{\omega}_t=\omega_w\\ \dot{\textbf{g}}_t=0\end{array}
p˙t=vtvt=Rt(am−abt−an)+gtq˙t=21qt⊗(ωm−ωbt−ωn)a˙bt=awω˙t=ωwg˙t=0
设
x
˙
t
=
f
t
(
x
t
,
u
,
w
)
{\dot{\mathbf{x}}}_{t}=f_{t}(\mathbf{x}_{t},\mathbf{u},\mathbf{w})
x˙t=ft(xt,u,w),该系统受IMU噪声读数
u
m
\mathbf{u}_{m}
um的控制,并受白高斯噪声
w
\textbf{w}
w的扰动
x
t
=
[
p
t
v
t
q
t
a
d
g
t
]
u
=
[
a
m
−
a
n
ω
m
−
ω
n
]
v
=
[
a
w
ω
w
]
\mathbf{x}_t=\begin{bmatrix}\mathbf{p}_t\\ \mathbf{v}_t\\ \mathbf{q}_t\\ \mathbf{a}_d\\ \mathbf{g}_t\\ \end{bmatrix}\quad\mathbf{u}=\begin{bmatrix}\mathbf{a}_m-\mathbf{a}_n\\ \mathbf{\omega}_m-\mathbf{\omega}_n\end{bmatrix}\quad\mathbf{v}=\begin{bmatrix}\mathbf{a}_w\\ \mathbf{\omega}_w\end{bmatrix}
xt=
ptvtqtadgt
u=[am−anωm−ωn]v=[awωw]
重力矢量 g t \mathbf{g}_{t} gt由滤波器估计。系统从一个固定且任意已知的初始方向开始,即 q t ( t = 0 ) = q 0 \mathbf{q}_t(t=0)=\mathbf{q}_0 qt(t=0)=q0,由于初始方向通常不在水平面上,因此初始重力向量通常是未知的,通常将初始方向设为 q 0 = ( 1 , 0 , 0 , 0 ) \mathbf{q}_0=(1,0,0,0) q0=(1,0,0,0)。因此, R 0 = R { q 0 } = I \mathbf{R}_0=\mathbf{R}\{\mathbf{q}_0\}=\mathbf{I} R0=R{q0}=I。还可以从系统中删除与重力相关的所有方程,并采用更传统的方法,即假设重力向量为 g = Δ ( 0 , 0 , − 9.8 x x ) \mathbf{g}\stackrel{\Delta}{=}\left(0,0,-9.8x x\right) g=Δ(0,0,−9.8xx),其中 x x xx xx是实验现场重力向量的适当小数位数,同时假设初始方向 q 0 \mathbf{q}_0 q0是不确定的。
名义状态运动学
名义状态运动学对应于没有噪声或扰动的建模系统。
p
˙
=
v
v
˙
=
R
(
a
m
−
a
b
)
+
g
q
˙
=
1
2
a
⊗
(
ω
m
−
ω
b
)
a
˙
b
=
0
ω
˙
b
=
0
g
˙
=
0
\begin{array}{l}\dot{\mathbf{p}}=\mathbf{v}\\ \dot{\mathbf{v}}=\mathbf{R}(\mathbf{a}_m-\mathbf{a}_b)+\mathbf{g}\\ \dot{\mathbf{q}}=\frac{1}{2}\mathbf{a}\otimes(\mathbf{\omega}_m-\mathbf{\omega}_b)\\ \dot{\mathbf{a}}_b=0\\ \dot{\mathbf{\omega}}_b=0\\ \dot{\mathbf{g}}=0\end{array}
p˙=vv˙=R(am−ab)+gq˙=21a⊗(ωm−ωb)a˙b=0ω˙b=0g˙=0
误差状态运动学
目标是确定误差状态的线性化动力学。
δ
˙
p
=
δ
v
δ
˙
v
=
−
R
[
a
m
−
a
b
]
×
δ
θ
−
R
δ
a
b
+
δ
g
−
R
a
n
δ
˙
θ
=
−
[
ω
m
−
ω
b
]
×
δ
θ
−
δ
ω
b
−
ω
n
δ
˙
a
b
=
a
w
δ
˙
ω
b
=
ω
w
δ
g
˙
=
0
\begin{aligned}&\dot{\delta}\textbf{p}=\delta\textbf{v}\\ &\dot{\delta}\textbf{v}=-\textbf{R}\left[\textbf{a}_m-\textbf{a}_b\right]_\times\delta\boldsymbol{\theta}-\textbf{R}\delta\textbf{a}_b+\delta\textbf{g}-\textbf{R}\textbf{a}_n\\ &\dot{\delta}{\boldsymbol{\theta}}=-\left[\omega_m-\omega_b\right]_\times\delta\boldsymbol{\theta}-\delta\omega_b-\omega_n\\ &\dot{\delta}{{\textbf{a}}}_b={\textbf{a}}_w\\ &\dot{\delta}{{\omega}}_b={\omega}_w\\ &{{\delta}}{\dot{\textbf{g}}}=0\end{aligned}
δ˙p=δvδ˙v=−R[am−ab]×δθ−Rδab+δg−Ranδ˙θ=−[ωm−ωb]×δθ−δωb−ωnδ˙ab=awδ˙ωb=ωwδg˙=0
其中,速度和方向误差的方程需要对角速度和重力矢量的误差方程进行一些操作,以获得线性化的动力学。
离散时间内的系统运动学
上述微分方程需要被积分成差分方程,以考虑离散的时间间隔∆t>0,积分方法可以有所不同。在某些情况下,可以使用精确的闭式解决方案。在其他情况下,可以使用不同精度的数值积分方法。请参见附录中有关积分方法的相关细节。
需要对以下子系统进行积分:
- 名义状态。
- 误差状态。
(a)确定性部分:状态动态和控制。
(b)随机部分:噪声和扰动。
对于名义状态和误差状态,需要对确定性部分和随机部分进行积分。确定性部分包括状态动态和控制,随机部分包括噪声和扰动。对于不同的子系统和积分方法,需要进行适当的选择和调整,以获得准确的结果。
名义状态的运动学
p ← p + v Δ t + 1 2 ( R ( a m − a b ) + g ) Δ t 2 v ← v + ( R ( a m − a b ) + g ) Δ t q ← q ⊗ q { ( ω m − ω b ) Δ t } a b ← a b g ← g \begin{array}{l}\textbf{p}\leftarrow\textbf{p}+\textbf{v}\Delta t+\frac{1}{2}(\textbf{R}(\textbf{a}_m-\textbf{a}_b)+\textbf{g})\Delta t^2\\ \textbf{v}\leftarrow\textbf{v}+(\textbf{R}(\textbf{a}_m-\textbf{a}_b)+\textbf{g})\Delta t\\ \textbf{q}\leftarrow\textbf{q}\otimes\textbf{q}\{(\omega_m-\omega_b)\Delta t\}\\ \textbf{a}_b\leftarrow\textbf{a}_b\\ \textbf{g}\leftarrow\textbf{g}\end{array} p←p+vΔt+21(R(am−ab)+g)Δt2v←v+(R(am−ab)+g)Δtq←q⊗q{(ωm−ωb)Δt}ab←abg←g
其中, x ← f ( x , ∙ ) x\leftarrow f(x,\bullet) x←f(x,∙)代表 x k + 1 = f ( x k , ∙ k ) x_{k+1}=f(x_k,\bullet_k) xk+1=f(xk,∙k)类型的时间更新, R ≜ R { q } \mathbf{R} \triangleq\mathbf{R}\{\mathbf{q}\} R≜R{q}是与当前名义方向 q \mathbf{q} q相关的旋转矩阵, q { v } \mathbf{q}\{v\} q{v}是与旋转 v v v相关的四元数。
误差状态的运动学
δ p ← δ p + δ v Δ t δ v ← δ v + ( − R [ a m − a b ] × δ θ − R δ a b + δ g ) Δ t + v i δ θ ← R ⊤ { ( ω m − ω b ) Δ t } δ θ − δ ω b Δ t + θ i δ a b ← δ a b + a i δ ω b ← δ ω b + ω i δ g ← δ g \begin{array}{l}\delta\textbf{p}\leftarrow\delta\textbf{p}+\delta\textbf{v}\Delta t\\ \delta\textbf{v}\leftarrow\delta\textbf{v}+(-\textbf{R}\left[\textbf{a}_m-\textbf{a}_b\right]\times\delta\boldsymbol{\theta}-\textbf{R}\delta\textbf{a}_b+\delta\textbf{g})\Delta t+\textbf{v}_\textbf{i}\\ \delta\boldsymbol{\theta}\leftarrow\textbf{R}^\top\left\{(\omega_m-\omega_b)\Delta t\right\}\delta\boldsymbol{\theta}-\delta\omega_b\Delta t+\boldsymbol{\theta}_\textbf{i}\\ \delta\textbf{a}_b\leftarrow\delta\textbf{a}_b+\textbf{a}_\textbf{i}\\ \delta \omega_b\leftarrow\delta \omega_b+\omega_\textbf{i}\\ \delta\textbf{g}\leftarrow\delta\textbf{g}\end{array} δp←δp+δvΔtδv←δv+(−R[am−ab]×δθ−Rδab+δg)Δt+viδθ←R⊤{(ωm−ωb)Δt}δθ−δωbΔt+θiδab←δab+aiδωb←δωb+ωiδg←δg
其中,,
v
i
\textbf{v}_\textbf{i}
vi、
θ
i
\boldsymbol{\theta}_\textbf{i}
θi、
a
i
\textbf{a}_\textbf{i}
ai和
ω
i
\omega_\textbf{i}
ωi是随机脉冲,作用于速度、方向和偏差估计上,由白高斯过程建模。它们的均值为零,协方差矩阵是通过将
a
n
\mathbf{a}_{n}
an、
ω
n
\omega_{n}
ωn、
a
ω
\mathbf{a}_{\omega}
aω和
ω
ω
\omega_{\omega}
ωω的协方差积分到步长时间
Δ
t
\Delta t
Δt上得到的。
V
i
=
σ
a
~
n
2
Δ
t
2
I
[
m
2
/
s
2
]
Θ
i
=
σ
ω
~
n
2
Δ
t
2
I
[
r
a
d
2
]
A
i
=
σ
a
w
2
Δ
t
I
[
m
2
/
s
4
]
Ω
i
=
σ
ω
w
2
Δ
t
I
[
r
a
d
2
/
s
2
]
\begin{aligned}\textbf{V}_\textbf{i}&=\sigma_{\tilde{a}_n}^2\Delta t^2\textbf{I}&\quad&[m^2/s^2]\\ \boldsymbol{\Theta}_\textbf{i}&=\sigma_{\tilde{\omega}_n}^2\Delta t^2\textbf{I}&\quad&[rad^2]\\ \textbf{A}_\textbf{i}&=\sigma_{{a}_w}^2\Delta t\textbf{I}&\quad&[m^2/s^4]\\ \boldsymbol{\Omega}_\textbf{i}&=\sigma_{{\omega}_w}^2\Delta t\textbf{I}&\quad&[rad^2/s^2]\\ \end{aligned}
ViΘiAiΩi=σa~n2Δt2I=σω~n2Δt2I=σaw2ΔtI=σωw2ΔtI[m2/s2][rad2][m2/s4][rad2/s2]
其中, σ a ~ n [ m 2 / s 2 ] \sigma_{\tilde{a}_n}[m^2/s^2] σa~n[m2/s2], σ ω ~ n [ r a d / s ] \sigma_{\tilde{\omega}_n}[rad/s] σω~n[rad/s], σ a w [ m / s 2 s ] \sigma_{{a}_w}\left[m/s^2\sqrt{s}\right] σaw[m/s2s]和 σ ω w [ r a d / s s ] \sigma_{{\omega}_w}\left[rad/s\sqrt{s}\right] σωw[rad/ss]将根据IMU数据表中的信息或实验测量结果来确定。
误差状态的雅可比矩阵和扰动矩阵
雅可比矩阵可以通过对前一节中的误差状态差分方程进行简单的推导得到。为了以简洁的形式书写这些方程,考虑名义状态向量
x
\mathbf{x}
x、误差状态向量
δ
x
\delta\mathbf{x}
δx、输入向量
u
m
\mathbf{u}_m
um和扰动脉冲向量
i
\mathbf{i}
i,如下所示
x
=
[
p
v
q
a
b
ω
b
g
]
,
δ
x
=
[
δ
p
δ
v
δ
θ
δ
a
b
δ
ω
b
δ
g
]
,
u
m
=
[
a
m
ω
m
]
,
i
=
[
v
i
θ
i
a
i
ω
i
]
\mathbf{x}=\begin{bmatrix}\mathbf{p}\\ \mathbf{v}\\ \mathbf{q}\\ \mathbf{a}_b\\\mathbf{\omega}_b\\ \mathbf{g}\end{bmatrix},\quad\delta\mathbf{x}=\begin{bmatrix}\delta\mathbf{p}\\ \delta\mathbf{v}\\\delta\boldsymbol{\theta}\\ \delta\mathbf{a}_b\\ \delta\mathbf{\omega}_b\\ \delta\mathbf{g}\end{bmatrix},\quad\mathbf{u}_m=\begin{bmatrix}\mathbf{a}_m\\ \mathbf{\omega}_m\end{bmatrix},\quad\mathbf{i}=\begin{bmatrix}\mathbf{v}_\mathbf{i}\\ \boldsymbol{\theta}_\mathbf{i} \\\mathbf{a}_\mathbf{i}\\ \mathbf{\omega}_\mathbf{i}\end{bmatrix}
x=
pvqabωbg
,δx=
δpδvδθδabδωbδg
,um=[amωm],i=
viθiaiωi
误差状态系统改写成:
δ
x
←
f
(
x
,
δ
x
,
u
m
,
i
)
=
F
x
(
x
,
u
m
)
⋅
δ
x
+
F
i
⋅
i
\delta\textbf{x}\leftarrow f(\textbf{x},\delta\textbf{x},\textbf{u}_m,\textbf{i})=\textbf{F}_\textbf{x}(\textbf{x},\textbf{u}_m)\cdot\delta\textbf{x}+\textbf{F}_\textbf{i}\cdot\textbf{i}
δx←f(x,δx,um,i)=Fx(x,um)⋅δx+Fi⋅i
ESKF的预测方程:
δ
x
^
←
F
x
(
x
,
u
m
)
⋅
δ
x
^
P
←
F
x
P
F
x
⊤
+
F
i
Q
i
F
i
⊤
\begin{array}{l}\hat{{\delta}{\mathbf{x}}}\leftarrow\mathbf{F}_{\mathbf{x}}(\mathbf{x},\mathbf{u}_m)\cdot\hat{{\delta}{\mathbf{x}}}\\ \mathbf{P}\leftarrow\mathbf{F}_{\mathbf{x}}\mathbf{P}\mathbf{F}_{\mathbf{x}}^{\top}+\mathbf{F}_{\mathbf{i}}\mathbf{Q}_{\mathbf{i}}\mathbf{F}_{\mathbf{i}}^{\top}\end{array}
δx^←Fx(x,um)⋅δx^P←FxPFx⊤+FiQiFi⊤
其中,
δ
x
∼
N
{
δ
x
^
,
P
}
\delta\mathbf{x}\sim\mathcal{N}\{\hat{\delta\mathbf{x}},\mathbf{P}\}
δx∼N{δx^,P};
F
x
\mathbf{F}_{\mathbf{x}}
Fx和
F
i
\mathbf{F}_{\mathbf{i}}
Fi是
f
(
)
f()
f()相对于误差和扰动矢量的雅可比矩阵;
Q
i
\mathbf{Q}_{\mathbf{i}}
Qi是扰动脉冲的协方差矩阵。
上述雅可比矩阵和协方差矩阵的表达方式详见下文。这里出现的所有与状态有关的值都是直接从名义状态中提取的。
F
x
=
∂
f
∂
δ
x
∣
x
,
u
m
=
[
I
I
Δ
t
0
0
0
0
0
I
−
R
[
a
m
−
a
b
]
×
Δ
t
−
R
Δ
t
0
I
Δ
t
0
0
R
⊤
{
(
ω
m
−
ω
b
)
Δ
t
}
0
−
I
Δ
t
0
0
0
0
I
0
0
0
0
0
0
I
0
0
0
0
0
0
I
]
\textbf{F}_\textbf{x}=\frac{\partial f}{\partial\delta\textbf{x}}\bigg|_{\textbf{x},\textbf{u}_m}=\begin{bmatrix}\textbf{I}&\textbf{I}\Delta t&0&0&0&0\\ 0&\textbf{I}&-\textbf{R}[\textbf{a}_m-\textbf{a}_b]_\times\Delta t&-\textbf{R}\Delta t&0&\textbf{I}\Delta t\\ 0&0&\textbf{R}^\top\{(\omega_m-\omega_b)\Delta t\}&0&-\textbf{I}\Delta t&0\\ 0&0&0&\textbf{I}&0&0\\ 0&0&0&0&\textbf{I}&0\\ 0&0&0&0&0&\textbf{I}\end{bmatrix}\\
Fx=∂δx∂f
x,um=
I00000IΔtI00000−R[am−ab]×ΔtR⊤{(ωm−ωb)Δt}0000−RΔt0I0000−IΔt0I00IΔt000I
F i = ∂ f ∂ i ∣ x , u m = [ 0 0 0 0 I 0 0 0 0 I 0 0 0 0 I 0 0 0 0 I 0 0 0 0 ] , Q i = [ V i 0 0 0 0 Θ i 0 0 0 0 A i 0 0 0 0 Ω i ] \textbf{F}_\textbf{i}=\frac{\partial f}{\partial\textbf{i}}\bigg|_{\textbf{x},\textbf{u}_m}=\begin{bmatrix}0&0&0&0\\\textbf{I}&0&0&0\\ 0&\textbf{I}&0&0\\ 0&0&\textbf{I}&0\\ 0&0&0&\textbf{I}\\ 0&0&0&0\end{bmatrix}\quad,\quad \textbf{Q}_\textbf{i}=\begin{bmatrix}\textbf{V}_\textbf{i}&0&0&0\\ 0&\boldsymbol{\Theta}_\textbf{i}&0&0\\ 0&0&\textbf{A}_\textbf{i}&0\\ 0&0&0&\boldsymbol{\Omega}_\textbf{i}\end{bmatrix} Fi=∂i∂f x,um= 0I000000I000000I000000I0 ,Qi= Vi0000Θi0000Ai0000Ωi
将IMU与互补传感器数据融合
当出现比IMU更多种类的信息,例如GPS或视觉信息时,会对ESKF进行纠正。在设计良好的系统中,这应该使IMU偏差可观测,并允许ESKF正确估计它们。有无数的可能性,最流行的是GPS + IMU,单眼视觉 + IMU和立体视觉 + IMU。
虽然IMU信息到目前为止用于向ESKF进行预测,但此其他信息用于纠正滤波器,从而观测IMU偏差误差。纠正包括三个步骤:
- 通过滤波器校正观测误差状态,
- 注入观测到的误差到名义状态中,并
- 重置误差状态。
通过滤波校正观测误差状态
假设像往常一样,传感器提供的信息取决于状态,例如:
y
=
h
(
x
t
)
+
v
\mathbf{y}=h(\mathbf{x}_{t})+v
y=h(xt)+v
其中
h
(
)
h()
h()是系统状态(真实状态)的一般非线性函数,
v
v
v是具有协方差
V
\textbf{V}
V的白高斯噪声,
v
∼
N
{
0
,
V
}
v\sim\mathcal{N}\{0,\mathbf{V}\}
v∼N{0,V}
滤波器是在估计错误状态,因此,滤波器的修正方程为:
K
=
PH
⊤
(
H
PH
⊤
+
V
)
−
1
δ
x
^
←
K
(
y
−
h
(
x
^
t
)
)
P
←
(
I
−
KH
)
P
\begin{aligned}\textbf{K}&=\textbf{PH}^\top(\textbf{H}\textbf{PH}^\top+\textbf{V})^{-1}\\ &\hat{{\delta}\textbf{x}}\leftarrow\textbf{K}(\textbf{y}-h(\hat{\textbf{x}}_t))\\ &\textbf{P}\leftarrow(\textbf{I}-\textbf{K}\textbf{H})\textbf{P}\end{aligned}
K=PH⊤(HPH⊤+V)−1δx^←K(y−h(x^t))P←(I−KH)P
要求雅可比矩阵
H
\textbf{H}
H相对于误差状态
δ
x
{\delta}\textbf{x}
δx定义,并在最佳真实状态估计
x
^
t
=
x
⊕
δ
x
^
{\hat{\mathbf{x}}}_{t}=\mathbf{x}\oplus\hat{\delta\mathbf{x}}
x^t=x⊕δx^处进行评估。由于误差状态平均值在此阶段为零(还没有观测到它),因此有
x
^
t
=
x
{\hat{\mathbf{x}}}_{t}=\mathbf{x}
x^t=x,并且可以使用名义误差
x
\textbf{x}
x作为评估点,从而得出以下结果:
H
≡
∂
h
∂
δ
x
∣
x
\mathbf{H}\equiv\left.\frac{\partial h}{\partial\delta\mathbf{x}}\right|_{\mathbf{x}}
H≡∂δx∂h
x
滤波器校正计算的雅可比矩阵
上述的雅可比矩阵可以通过多种方式计算。最具说明性的方法是利用链式法则,
H
≜
∂
h
∂
δ
x
∣
x
=
∂
h
∂
x
t
∣
x
∂
x
t
∂
δ
x
∣
x
=
H
x
X
δ
x
\textbf{H}\triangleq\left.\frac{\partial h}{\partial\delta\textbf{x}}\right|_{\textbf{x}}=\left.\frac{\partial h}{\partial\textbf{x}_t}\right|_{\textbf{x}}\left.\frac{\partial\textbf{x}_t}{\partial\delta\textbf{x}}\right|_{\textbf{x}}=\textbf{H}_{\textbf{x}}\textbf{X}_{\delta\textbf{x}}
H≜∂δx∂h
x=∂xt∂h
x∂δx∂xt
x=HxXδx
其中,
H
≜
∂
h
∂
δ
x
∣
x
\textbf{H}\triangleq\left.\frac{\partial h}{\partial\delta\textbf{x}}\right|_{\textbf{x}}
H≜∂δx∂h
x是
h
(
)
h()
h()相对于其自身参数的标准雅各布系数(即在常规EKF中使用的雅各布系数)。链式规则的第一部分取决于所使用的特定传感器的测量功能.
第二部分,
X
δ
x
≜
∂
x
t
∂
δ
x
∣
x
\textbf{X}_{\delta\textbf{x}}\triangleq\left.\frac{\partial\textbf{x}_t}{\partial\delta\textbf{x}}\right|_{\textbf{x}}
Xδx≜∂δx∂xt
x,是真实状态相对于误差状态的雅可比矩阵。这部分可以在这里导出,因为它只取决于状态的ESKF组成。有导数,
X
δ
x
=
[
∂
(
p
+
δ
p
)
∂
δ
p
∂
(
v
+
δ
v
)
∂
δ
v
0
∂
(
q
⊗
δ
q
)
∂
δ
θ
∂
(
a
b
+
δ
a
b
)
∂
δ
a
b
0
∂
(
ω
b
+
δ
ω
b
)
∂
δ
ω
b
∂
(
g
+
δ
g
)
∂
δ
g
]
\mathbf{X}_{\delta \mathbf{x}}=\left[\begin{array}{ccccc} \frac{\partial(\mathbf{p}+\delta \mathbf{p})}{\partial \delta \mathbf{p}} & & & & \\ & \frac{\partial(\mathbf{v}+\delta \mathbf{v})}{\partial \delta \mathbf{v}} & & & & 0 & \\ & & \frac{\partial(\mathbf{q} \otimes \delta \mathbf{q})}{\partial \delta \boldsymbol{\theta}} & & \\ & & & & \frac{\partial\left(\mathbf{a}_{b}+\delta \mathbf{a}_{b}\right)}{\partial \delta \mathbf{a}_{b}} & \\ &0 & & & & \frac{\partial\left(\boldsymbol{\omega}_{b}+\delta \boldsymbol{\omega}_{b}\right)}{\partial \delta \boldsymbol{\omega}_{b}} & \\ & & & & & &\frac{\partial(\mathbf{g}+\delta \mathbf{g})}{\partial \delta \mathbf{g}} \end{array}\right]
Xδx=
∂δp∂(p+δp)∂δv∂(v+δv)0∂δθ∂(q⊗δq)∂δab∂(ab+δab)0∂δωb∂(ωb+δωb)∂δg∂(g+δg)
这导致了除了
4
×
3
4×3
4×3四元数项
Q
δ
θ
=
∂
(
q
⊗
δ
q
)
/
∂
δ
θ
\textbf{Q}_{\delta\boldsymbol{\theta}}=\partial(\textbf{q}\otimes\delta\textbf{q})/\partial\delta\boldsymbol{\theta}
Qδθ=∂(q⊗δq)/∂δθ之外的所有恒等式
3
×
3
3×3
3×3块(例如,
∂
(
p
+
δ
p
)
∂
δ
p
=
I
3
\frac{\partial(\textbf{p}+\delta\textbf{p})}{\partial\delta\textbf{p}}=\textbf{I}_3
∂δp∂(p+δp)=I3)。因此,有这样的形式,
X
δ
x
≜
∂
x
t
∂
δ
x
∣
x
=
[
I
6
0
0
0
Q
δ
θ
0
0
0
I
9
]
\textbf{X}_{\delta\textbf{x}}\triangleq\left.\frac{\partial\textbf{x}_t}{\partial\delta\textbf{x}}\right|_{\textbf{x}}=\begin{bmatrix}\textbf{I}_6&0&0\\ 0&\textbf{Q}_{\delta\theta}&0\\ 0&0&\textbf{I}_9\end{bmatrix}
Xδx≜∂δx∂xt
x=
I6000Qδθ000I9
四元数项
Q
δ
θ
\textbf{Q}_{\delta\boldsymbol{\theta}}
Qδθ可推导为:
Q
δ
θ
=
1
2
[
−
q
x
−
q
y
−
q
z
q
w
−
q
z
q
y
q
z
q
w
−
q
x
−
q
y
q
x
q
w
]
\textbf{Q}_{\delta\boldsymbol{\theta}}=\frac{1}{2}\begin{bmatrix}-q_x&-q_y&-q_z\\ q_w&-q_z&q_y\\ q_z&q_w&-q_x\\ -q_y&q_x&q_w\end{bmatrix}
Qδθ=21
−qxqwqz−qy−qy−qzqwqx−qzqy−qxqw
将观测误差注入名义状态
在ESKF更新之后,使用适当的组合(求和或四元数乘积)将观测误差状态与名义状态进行更新。
x
←
x
⊕
δ
x
^
\textbf{x}\leftarrow\textbf{x}\oplus\hat{\delta\textbf{x}}
x←x⊕δx^
也就是,
p
←
p
+
δ
p
^
v
←
v
+
δ
v
^
q
←
q
⊗
q
{
δ
θ
^
}
a
b
←
a
b
+
δ
a
b
^
ω
b
←
ω
b
+
δ
ω
^
b
g
←
g
+
δ
g
^
\begin{array}{l} \mathbf{p} \leftarrow \mathbf{p}+\hat{\delta \mathbf{p}} \\ \mathbf{v} \leftarrow \mathbf{v}+\hat{\delta \mathbf{v}} \\ \mathbf{q} \leftarrow \mathbf{q} \otimes \mathbf{q}\{\hat{\delta \boldsymbol{\theta}}\} \\ \mathbf{a}_{b} \leftarrow \mathbf{a}_{b}+\hat{\delta {\mathbf{a}}_{b}} \\ \boldsymbol{\omega}_{b} \leftarrow \boldsymbol{\omega}_{b}+\delta \hat{\boldsymbol{\omega}}_{b} \\ \mathbf{g} \leftarrow \mathbf{g}+\hat{\delta \mathbf{g}} \end{array}
p←p+δp^v←v+δv^q←q⊗q{δθ^}ab←ab+δab^ωb←ωb+δω^bg←g+δg^
重置ESKF
在将误差注入到名义状态后,误差状态均值 δ x ^ \hat{\delta\textbf{x}} δx^将被重置。这在方向部分尤其重要,因为新的方向误差将相对于新名义状态的方向框架局部表示。为了使ESKF更新完整,误差的协方差需要根据这个修改进行更新。
称误差重置函数为
g
(
)
g()
g()。写法如下:
δ
x
←
g
(
δ
x
)
=
δ
x
⊖
δ
x
^
\delta\mathbf{x}\leftarrow g(\delta\mathbf{x})=\delta\mathbf{x}\ominus\hat{\delta\mathbf{x}}
δx←g(δx)=δx⊖δx^
其中,
⊖
\ominus
⊖代表
⊕
\oplus
⊕的组合的逆。因此,ESKF的误差重置操作为:
δ
x
^
←
0
P
←
G
PG
⊤
\begin{array}{l}\hat{\delta\text{x}}\leftarrow0\\ \text{P}\leftarrow\text{G}\text{PG}^\top\end{array}
δx^←0P←GPG⊤
其中G是雅可比矩阵,定义如下:
G
≜
∂
g
∂
δ
x
∣
δ
x
^
\textbf{G}\triangleq\frac{\partial g}{\partial\delta\textbf{x}}\bigg|_{\hat{\delta\textbf{x}}}
G≜∂δx∂g
δx^
类似于上述更新雅可比矩阵所发生的情况,这个雅可比矩阵在除了方向误差之外的所有对角块上都是相同的。在此给出完整的表达式,
∂
δ
θ
+
/
∂
δ
θ
=
I
−
[
1
2
δ
θ
^
}
]
×
\partial\delta\boldsymbol{\theta}^{+}/\partial\delta\boldsymbol{\theta}=\mathbf{I}-\left[{\frac{1}{2}}{\hat{\delta \boldsymbol{\theta}}\}}\right]_{\times}
∂δθ+/∂δθ=I−[21δθ^}]×,
G
=
[
I
6
0
0
0
I
−
[
1
2
δ
θ
^
}
]
×
0
0
0
I
9
]
\textbf{G}=\begin{bmatrix}\textbf{I}_6&0&0\\ 0&\textbf{I}-\begin{bmatrix}\frac{1}{2}\hat{\delta \boldsymbol{\theta}}\}\end{bmatrix}_\times&0\\ 0&0&\textbf{I}_9\end{bmatrix}
G=
I6000I−[21δθ^}]×000I9
在大多数情况下,误差项ˇδθ可以忽略,只产生一个雅可比矩阵
G
=
I
18
\textbf{G}=\textbf{I}_{18}
G=I18,从而导致微小的误差复位。
使用全局误差的ESKF
本节探讨在全局参考系下定义角度误差与迄今为止所使用的局部定义之间的差异。在全局参考系下定义角度误差δθ意味着在左侧进行组合,即:
q
t
=
δ
q
⊗
q
=
q
{
δ
θ
}
⊗
q
\mathbf{q}_{t}=\delta\mathbf{q}\otimes\mathbf{q}=\mathbf{q}\{\delta\mathbf{\boldsymbol{\theta}}\}\otimes\mathbf{q}
qt=δq⊗q=q{δθ}⊗q
这里,无论角度误差是全局定义还是局部定义,仍然保留角速率向量
ω
ω
ω的局部定义,即在连续时间中,
q
˙
=
1
2
q
⊗
ω
\dot{\mathbf{q}}=\frac{1}{2}\mathbf{q}\otimes\mathbf{\omega}
q˙=21q⊗ω,因此在离散时间中,
q
←
q
⊗
q
{
ω
Δ
t
}
\textbf{q}\leftarrow\textbf{q}\otimes\textbf{q}\{\omega\Delta t\}
q←q⊗q{ωΔt}。这是为了方便起见,因为陀螺仪提供的角速率测量值是以体坐标系为基准的,也就是局部坐标系。
连续时间中的系统运动学
真实状态和名义状态运动学
真实运动学和名义运动学不包含误差,它们的方程不变。
误差状态运动学
误差状态运动学方程
δ
p
˙
=
δ
v
δ
v
˙
=
−
[
R
(
a
m
−
b
a
b
)
]
×
δ
θ
−
R
δ
a
b
+
δ
g
−
Ra
n
δ
θ
˙
=
−
R
δ
ω
b
−
R
ω
n
δ
a
˙
b
=
a
w
δ
ω
˙
b
=
ω
w
δ
g
˙
=
0
\begin{aligned}&\dot{\delta\textbf{p}}=\delta\textbf{v}\\ &\dot{\delta\textbf{v}}=-\left[\textbf{R}(\textbf{a}_m-b\textbf{a}_b)\right]_\times\delta\boldsymbol{\theta}-\textbf{R}\delta\textbf{a}_b+\delta\textbf{g}-\textbf{Ra}_n\\ &\dot{\delta\boldsymbol{\theta}}=-\textbf{R}\delta\omega_b-\textbf{R}\omega_n\\ &\delta\dot{\textbf{a}}_b=\textbf{a}_w\\ &\delta{\dot{\omega}}_b=\omega_w\\ &\dot{\delta\textbf{g}}=0\end{aligned}
δp˙=δvδv˙=−[R(am−bab)]×δθ−Rδab+δg−Ranδθ˙=−Rδωb−Rωnδa˙b=awδω˙b=ωwδg˙=0
离散时间中的系统运动学
名义状态
名义状态方程不涉及误差,因此与角度误差局部定义情况下的方程相同。
误差状态
δ p ← δ p + δ v Δ t δ v ← δ v + ( − [ R ( a n − a b ) ] × δ θ − R δ a b + δ g ) Δ t + v i δ θ ← δ θ − R δ ω b Δ t + θ i δ a b ← δ a b + a i δ ω b ← δ ω b + ω i δ g ← δ g \begin{array}{l}\delta\textbf{p}\leftarrow\delta\textbf{p}+\delta\textbf{v}\Delta t\\ \delta\textbf{v}\leftarrow\delta\textbf{v}+(-\left[\textbf{R}(\textbf{a}_n-\textbf{a}_b)\right]_{\times}\delta\boldsymbol{\theta}-\textbf{R}\delta\textbf{a}_b+\delta\textbf{g})\Delta t+\textbf{v}_\textbf{i}\\ \delta\boldsymbol{\theta}\leftarrow\delta\boldsymbol{\theta}-\textbf{R}\delta\omega_b\Delta t+\boldsymbol{\theta}_\textbf{i}\\ \delta\textbf{a}_b\leftarrow\delta\textbf{a}_b+\textbf{a}_\textbf{i}\\ \delta\omega_b\leftarrow\delta\omega_b+\omega_\textbf{i}\\ \delta\textbf{g}\leftarrow\delta\textbf{g}\end{array} δp←δp+δvΔtδv←δv+(−[R(an−ab)]×δθ−Rδab+δg)Δt+viδθ←δθ−RδωbΔt+θiδab←δab+aiδωb←δωb+ωiδg←δg
误差状态雅可比矩阵和扰动矩阵
通过简单地推导上述方程,可以得到转移矩阵:
F
x
=
[
I
I
Δ
t
0
0
0
0
0
I
−
[
R
(
a
m
−
a
b
)
]
×
Δ
t
−
R
Δ
t
0
I
Δ
t
0
0
I
0
−
R
Δ
t
0
0
0
0
I
0
0
0
0
0
0
I
0
0
0
0
0
0
I
]
\mathbf{F}_{\mathbf{x}}=\begin{bmatrix}\mathbf{I}&\mathbf{I}\Delta t&0&0&0&0\\ 0&\mathbf{I}&-[\mathbf{R}(\mathbf{a}_{m}-\mathbf{a}_{b})]_{\times}\Delta t&-\mathbf{R}\Delta t&0&\mathbf{I}\Delta t\\ 0&0&\mathbf{I}&0&-\mathbf{R}\Delta t&0\\ 0&0&0&\mathbf{I}&0&0\\ 0&0&0&0&\mathbf{I}&0\\ 0&0&0&0&0&\mathbf{I}\end{bmatrix}
Fx=
I00000IΔtI00000−[R(am−ab)]×ΔtI0000−RΔt0I0000−RΔt0I00IΔt000I
扰动雅可比矩阵和扰动矩阵保持不变。
F
i
=
[
0
0
0
0
I
0
0
0
0
I
0
0
0
0
I
0
0
0
0
I
0
0
0
0
]
,
Q
i
=
[
V
i
0
0
0
0
Θ
i
0
0
0
0
A
i
0
0
0
0
Ω
i
]
\textbf{F}_\textbf{i}=\begin{bmatrix}0&0&0&0\\\textbf{I}&0&0&0\\ 0&\textbf{I}&0&0\\ 0&0&\textbf{I}&0\\ 0&0&0&\textbf{I}\\ 0&0&0&0\end{bmatrix}\quad,\quad \textbf{Q}_\textbf{i}=\begin{bmatrix}\textbf{V}_\textbf{i}&0&0&0\\ 0&\boldsymbol{\Theta}_\textbf{i}&0&0\\ 0&0&\textbf{A}_\textbf{i}&0\\ 0&0&0&\boldsymbol{\Omega}_\textbf{i}\end{bmatrix}
Fi=
0I000000I000000I000000I0
,Qi=
Vi0000Θi0000Ai0000Ωi
与互补传感器数据融合
当考虑全局角度误差时,涉及ESKF机制的融合方程仅略有变化。将通过ESKF校正中的误差状态观测、误差注入到名义状态以及重置步骤来修订这些变化。
误差状态观测
与局部误差定义相比,唯一的区别在于将方向与角度误差相关联的观测函数的雅可比矩阵块。
Q
δ
θ
≜
1
2
[
−
q
x
−
q
y
−
q
z
q
w
q
z
−
q
y
−
q
z
q
w
q
x
q
y
−
q
x
q
w
]
\mathbf{Q}_{\delta\theta}\triangleq\frac{1}{2}\begin{bmatrix}-q_x&-q_y&-q_z\\ q_w&q_z&-q_y\\ -q_z&q_w&q_x\\ q_y&-q_x&q_w\end{bmatrix}
Qδθ≜21
−qxqw−qzqy−qyqzqw−qx−qz−qyqxqw
将观测误差注入名义状态
p ← p + δ p v ← v + δ v q ← q { δ θ ^ } ⊗ q a b ← a b + δ a b ω b ← ω b + δ ω b g ← g + δ g \begin{array}{l}\textbf{p}\leftarrow\textbf{p}+\delta\textbf{p}\\ \textbf{v}\leftarrow\textbf{v}+\delta\textbf{v}\\ \textbf{q}\leftarrow\textbf{q}\{\hat{\delta \boldsymbol{\theta}}\}\otimes\textbf{q}\\ \textbf{a}_b\leftarrow\textbf{a}_b+\delta\textbf{a}_b\\ \omega_b\leftarrow\omega_b+\delta\omega_b\\ \textbf{g}\leftarrow\textbf{g}+\delta\textbf{g}\end{array} p←p+δpv←v+δvq←q{δθ^}⊗qab←ab+δabωb←ωb+δωbg←g+δg
重置ESKF
根据以下公式,ESKF误差均值被重置,协方差被更新:
δ
x
^
←
0
P
←
G
PG
⊤
\begin{array}{l}\hat{\delta\text{x}}\leftarrow0\\ \text{P}\leftarrow\text{G}\text{PG}^\top\end{array}
δx^←0P←GPG⊤
雅可比矩阵为:
G
=
[
I
6
0
0
0
I
+
[
1
2
δ
θ
^
}
]
×
0
0
0
I
9
]
\textbf{G}=\begin{bmatrix}\textbf{I}_6&0&0\\ 0&\textbf{I}+\begin{bmatrix}\frac{1}{2}\hat{\delta \boldsymbol{\theta}}\}\end{bmatrix}_\times&0\\ 0&0&\textbf{I}_9\end{bmatrix}
G=
I6000I+[21δθ^}]×000I9