一、VIO残差函数的构建
1. 系统所需的状态变量
为了节约计算量采用滑动窗口形式的 Bundle Adjustment, 在
i
i
i 时刻, 滑动窗口内待优化的系统状态量定义如下:
X
=
[
x
n
,
x
n
+
1
,
…
,
x
n
+
N
,
λ
m
,
λ
m
+
1
,
…
,
λ
m
+
M
]
x
i
=
[
p
w
b
i
,
q
w
b
i
,
v
i
w
,
b
a
b
i
,
b
g
b
i
]
?
,
i
∈
[
n
,
n
+
N
]
\begin{aligned} &\mathcal{X}=\left[\mathbf{x}_{n}, \mathbf{x}_{n+1}, \ldots, \mathbf{x}_{n+N}, \lambda_{m}, \lambda_{m+1}, \ldots, \lambda_{m+M}\right] \\ &\mathbf{x}_{i}=\left[\mathbf{p}_{w b_{i}}, \mathbf{q}_{w b_{i}}, \mathbf{v}_{i}^{w}, \mathbf{b}_{a}^{b_{i}}, \mathbf{b}_{g}^{b_{i}}\right]^{\top}, i \in[n, n+N] \end{aligned}
?X=[xn?,xn+1?,…,xn+N?,λm?,λm+1?,…,λm+M?]xi?=[pwbi??,qwbi??,viw?,babi??,bgbi??]?,i∈[n,n+N]? 其中:
-
x
i
\mathbf{x}_{i}
xi? 包含
i
i
i 时刻 IMU 机体的在惯性坐标系中的位置,速度,姿态, 以及 IMU 机体坐标系中的加速度和角速度的偏置量估计。
-
n
,
m
n, m
n,m 分别是机体状态量, 路标在滑动窗口里的起始时刻。
-
N
N
N 滑动窗口中关键帧数量。
-
M
M
M 是被滑动窗口内所有关键帧观测到的路标数量。
2. 视觉重投影误差
2.1 视觉重投影误差
定义: 一个特征点在归一化相机坐标系下 的估计值与观测值的差。
r
c
=
[
x
z
?
u
y
z
?
v
]
\mathbf{r}_{c}=\left[\begin{array}{l} \frac{x}{z}-u \\ \frac{y}{z}-v \end{array}\right]
rc?=[zx??uzy??v?] 其中,待估计的状态量为特征点的三维空间坐标
(
x
,
y
,
z
)
?
(x, y, z)^{\top}
(x,y,z)?, 观测值
(
u
,
v
)
?
(u, v)^{\top}
(u,v)? 为特征在相机归一化平面的坐标。
2.2 逆深度参数化
特征点在归一化相机坐标系与在相机坐标系下的坐标关系为:
[
x
y
z
]
=
1
λ
[
u
v
1
]
\left[\begin{array}{l} x \\ y \\ z \end{array}\right]=\frac{1}{\lambda}\left[\begin{array}{l} u \\ v \\ 1 \end{array}\right]
???xyz????=λ1????uv1???? 其中
λ
=
1
/
z
\lambda=1 / z
λ=1/z 称为逆深度。
- 使用逆深度的缘故是相较于深度值其
更符合高斯分布的特性 ; - 深度值一般很大,会影响优化的节奏,使用
倒数更能保证数值的稳定性 。
2.3 VIO 中基于逆深度的重投影误差
特征点逆深度在第
i
i
i 帧中初始化得 到,在第
j
j
j 帧又被观测到,预测其 在第
j
j
j 中的坐标为:
[
x
c
j
y
c
j
z
c
j
1
]
=
T
b
c
?
1
T
w
b
j
?
1
T
w
b
i
T
b
c
[
1
λ
u
c
i
1
λ
v
c
i
1
λ
1
]
\left[\begin{array}{c}x_{c_{j}} \\ y_{c_{j}} \\ z_{c_{j}} \\ 1\end{array}\right]=\mathbf{T}_{b c}^{-1} \mathbf{T}_{w b_{j}}^{-1} \mathbf{T}_{w b_{i}} \mathbf{T}_{b c}\left[\begin{array}{c}\frac{1}{\lambda} u_{c_{i}} \\ \frac{1}{\lambda} v_{c_{i}} \\ \frac{1}{\lambda} \\ 1\end{array}\right]
?????xcj??ycj??zcj??1??????=Tbc?1?Twbj??1?Twbi??Tbc??????λ1?uci??λ1?vci??λ1?1?????? 其中:
T
=
[
R
t
0
T
1
]
\begin{aligned} T=\left[\begin{array}{ll} R & t \\ 0^{T} & 1 \end{array}\right] \end{aligned}
T=[R0T?t1?]?
- 将i帧中观测到的数据变换到相机坐标系,将
相机坐标系变换到body坐标系 ; - 将第i个body坐标系变换到世界坐标系;
- 将世界坐标系变换到第j个body坐标系;
- 将
body坐标系变换到相机坐标系 ,得到第j帧的预测值。
这期间相对于纯视觉多了相机坐标系变换到body坐标系 ,然后再由body坐标系变换回相机坐标系 的过程。
视觉重投影误差为:
r
c
=
[
x
c
j
z
c
j
?
u
c
j
y
c
j
z
c
j
?
v
c
j
]
\begin{aligned} &\mathbf{r}_{c}=\left[\begin{array}{l} \frac{x_{c_{j}}}{z_{c_{j}}}-u_{c_{j}} \\ \frac{y_{c_{j}}}{z_{c_{j}}}-v_{c_{j}} \end{array}\right] \\ \end{aligned}
?rc?=???zcj??xcj????ucj??zcj??ycj????vcj???????
3. 预积分模型由来及意义
3.1 为什么需要预积分?
IMU的测量值为
ω
,
b
\omega, b
ω,b ,则有
ω
~
b
=
ω
b
+
b
g
+
n
g
a
~
b
=
q
b
w
(
a
w
+
g
w
)
+
b
a
+
n
a
\begin{gathered} \tilde{\omega}^{b}=\omega^{b}+b^{g}+n^{g} \\ \tilde{a}^{b}=q_{b w}\left(a^{w}+g^{w}\right)+b^{a}+n^{a} \end{gathered}
ω~b=ωb+bg+nga~b=qbw?(aw+gw)+ba+na? PVQ对时间的导数可写成
p
˙
w
b
t
=
v
t
w
v
˙
t
w
=
a
t
w
q
˙
w
b
t
=
q
w
b
t
?
[
0
1
2
ω
b
t
]
\begin{gathered} \dot{p}_{w b_{t}}=v_{t}^{w} \\ \dot{v}_{t}^{w}=a_{t}^{w} \\ \dot{q}_{w b_{t}}=q_{w b_{t}} \otimes\left[\begin{array}{c} 0 \\ \frac{1}{2} \omega^{b_{t}} \end{array}\right] \end{gathered}
p˙?wbt??=vtw?v˙tw?=atw?q˙?wbt??=qwbt???[021?ωbt??]? 根据上面的导数关系,可以从第
i
\mathrm{i}
i 时刻的
P
V
Q
\mathrm{PVQ}
PVQ ,通过对
I
M
U
\mathrm{IMU}
IMU 的测量值进行积分,得到第时刻的 PVQ:
p
w
b
j
=
p
w
b
i
+
v
i
w
Δ
t
+
?
t
∈
[
i
,
j
]
(
q
w
b
t
a
b
t
?
g
w
)
δ
t
2
v
j
w
=
v
i
w
+
∫
t
∈
[
i
,
j
]
(
q
w
b
t
a
b
t
?
g
w
)
δ
t
q
w
b
j
=
∫
t
∈
[
i
,
j
]
q
w
b
t
?
[
0
1
2
ω
b
t
]
δ
t
\begin{aligned} &\mathbf{p}_{w b_{j}}=\mathbf{p}_{w b_{i}}+\mathbf{v}_{i}^{w} \Delta t+\iint_{t \in[i, j]}\left(\mathbf{q}_{w b_{t}} \mathbf{a}^{b_{t}}-\mathbf{g}^{w}\right) \delta t^{2} \\ &\mathbf{v}_{j}^{w}=\mathbf{v}_{i}^{w}+\int_{t \in[i, j]}\left(\mathbf{q}_{w b_{t}} \mathbf{a}^{b_{t}}-\mathbf{g}^{w}\right) \delta t \\ &\mathbf{q}_{w b_{j}}=\int_{t \in[i, j]} \mathbf{q}_{w b_{t}} \otimes\left[\begin{array}{c} 0 \\ \frac{1}{2} \boldsymbol{\omega}^{b_{t}} \end{array}\right] \delta t \end{aligned}
?pwbj??=pwbi??+viw?Δt+?t∈[i,j]?(qwbt??abt??gw)δt2vjw?=viw?+∫t∈[i,j]?(qwbt??abt??gw)δtqwbj??=∫t∈[i,j]?qwbt???[021?ωbt??]δt?
- 问题: 为什么需要预积分?
每次
q
ω
b
t
q_{\omega b_{t}}
qωbt?? 更新后,都需要重新进行积分,运算量大 ;我们知道在紧耦合中待优化的状态变量是
p
p
p,
q
q
q,
v
v
v在优化的一次次迭代中这些变量是在不停变化的。由公式可以看到当
b
i
b_i
bi?时刻的状态发生改变的时候,我们需要一次次的积分来求取
b
j
b_j
bj?时刻的状态,这个是很耗费计算资源的,为了节省计算资源和时间,我们希望能够避免积分部分的重复计算。
3.2 怎么预积分?
一个很简单的公式转换,就可以将积分模型转为预积分模型:
q
w
b
t
=
q
w
b
i
?
q
b
i
b
t
\mathbf{q}_{w b_{t}}=\mathbf{q}_{w b_{i}} \otimes \mathbf{q}_{b_{i} b_{t}}
qwbt??=qwbi???qbi?bt?? 那么,PVQ 积分公式中的积分项则变成相对于第i 时刻的姿态 ,而不是相对于世界坐标系的姿态(注意,积分的变化 ):
p
w
b
j
=
p
w
b
i
+
v
i
w
Δ
t
?
1
2
g
w
Δ
t
2
+
q
w
b
i
?
t
∈
[
i
,
j
]
(
q
b
i
b
t
a
b
t
)
δ
t
2
v
j
w
=
v
i
w
?
g
w
Δ
t
+
q
w
b
?
∫
t
∈
[
i
,
j
]
(
q
b
z
b
t
a
b
t
)
δ
t
q
w
b
j
=
q
w
b
i
∫
t
∈
[
i
,
j
]
q
b
i
b
t
?
[
0
1
2
ω
b
t
]
δ
t
\begin{aligned} &\mathbf{p}_{w b_{j}}=\mathbf{p}_{w b_{i}}+\mathbf{v}_{i}^{w} \Delta t-\frac{1}{2} \mathbf{g}^{w} \Delta t^{2}+\mathbf{q}_{w b_{i}} \iint_{t \in[i, j]}\left(\mathbf{q}_{b_{i} b_{t}} \mathbf{a}^{b_{t}}\right) \delta t^{2} \\ &\mathbf{v}_{j}^{w}=\mathbf{v}_{i}^{w}-\mathbf{g}^{w} \Delta t+\mathbf{q}_{w b_{\imath}} \int_{t \in[i, j]}\left(\mathbf{q}_{b_{z} b_{t}} \mathbf{a}^{b_{t}}\right) \delta t \\ &\mathbf{q}_{w b_{j}}=\mathbf{q}_{w b_{i}} \int_{t \in[i, j]} \mathbf{q}_{b_{i} b_{t}} \otimes\left[\begin{array}{c} 0 \\ \frac{1}{2} \omega^{b_{t}} \end{array}\right] \delta t \end{aligned}
?pwbj??=pwbi??+viw?Δt?21?gwΔt2+qwbi???t∈[i,j]?(qbi?bt??abt?)δt2vjw?=viw??gwΔt+qwb???∫t∈[i,j]?(qbz?bt??abt?)δtqwbj??=qwbi??∫t∈[i,j]?qbi?bt???[021?ωbt??]δt? 这里的预积分量仅仅跟IMU 测量值有关 ,它将一段时间内的IMU 数据直接积分起来就得到了预积分量 :
α
b
i
b
j
=
?
t
∈
[
i
,
j
]
(
q
b
i
b
t
a
b
t
)
δ
t
2
β
b
i
b
j
=
∫
t
∈
[
i
,
j
]
(
q
b
i
b
t
a
b
t
)
δ
t
q
b
i
b
j
=
∫
t
∈
[
i
,
j
]
q
b
i
b
t
?
[
0
1
2
ω
b
t
]
δ
t
\begin{aligned} \boldsymbol{\alpha}_{b_{i} b_{j}} &=\iint_{t \in[i, j]}\left(\mathbf{q}_{b_{i} b_{t}} \mathbf{a}^{b_{t}}\right) \delta t^{2} \\ \boldsymbol{\beta}_{b_{i} b_{j}} &=\int_{t \in[i, j]}\left(\mathbf{q}_{b_{i} b_{t}} \mathbf{a}^{b_{t}}\right) \delta t \\ \mathbf{q}_{b_{i} b_{j}} &=\int_{t \in[i, j]} \mathbf{q}_{b_{i} b_{t}} \otimes\left[\begin{array}{c} 0 \\ \frac{1}{2} \boldsymbol{\omega}^{b_{t}} \end{array}\right] \delta t \end{aligned}
αbi?bj??βbi?bj??qbi?bj???=?t∈[i,j]?(qbi?bt??abt?)δt2=∫t∈[i,j]?(qbi?bt??abt?)δt=∫t∈[i,j]?qbi?bt???[021?ωbt??]δt? 重新整理下PVQ 的积分公式,有:
[
p
w
b
j
v
j
w
q
w
b
j
b
j
a
b
j
g
]
=
[
p
w
b
i
+
v
i
w
Δ
t
?
1
2
g
w
Δ
t
2
+
q
w
b
i
α
b
i
b
j
v
i
w
?
g
w
Δ
t
+
q
w
b
i
β
b
i
b
j
q
w
b
q
q
b
i
b
j
b
i
a
b
i
g
]
\left[\begin{array}{c} \mathbf{p}_{w b_{j}} \\ \mathbf{v}_{j}^{w} \\ \mathbf{q}_{w b_{j}} \\ \mathbf{b}_{j}^{a} \\ \mathbf{b}_{j}^{g} \end{array}\right]=\left[\begin{array}{c} \mathbf{p}_{w b_{i}}+\mathbf{v}_{i}^{w} \Delta t-\frac{1}{2} \mathbf{g}^{w} \Delta t^{2}+\mathbf{q}_{w b_{i}} \alpha_{b_{i} b_{j}} \\ \mathbf{v}_{i}^{w}-\mathbf{g}^{w} \Delta t+\mathbf{q}_{w b_{i}} \boldsymbol{\beta}_{b_{i} b_{j}} \\ \mathbf{q}_{w b_{\mathbf{q}}} \mathbf{q}_{b_{i} b_{j}} \\ \mathbf{b}_{i}^{a} \\ \mathbf{b}_{i}^{g} \end{array}\right]
???????pwbj??vjw?qwbj??bja?bjg?????????=???????pwbi??+viw?Δt?21?gwΔt2+qwbi??αbi?bj??viw??gwΔt+qwbi??βbi?bj??qwbq??qbi?bj??bia?big?????????
3.3 预积分是什么?
预积分其实是加速度和角速度引起的
b
j
b_j
bj?时刻的状态量关于
b
i
b_i
bi?时刻的增量 。 所谓的预积分,就是预先积分好,之后需要的时候拿来用。我们知道IMU测量的是线性加速度和角速度。状态量中速度是关于加速度的积分,位置是关于加速度的二次积分,姿态是关于角速度的积分。因此,在当前坐标系下(因为IMU的测量值就是在当前坐标系下的)把IMU测量值积分好,当需要的是,只需要根据情况进行坐标系转换和加减了。 可以看到预积分与待优化的状态量(忽略零偏)无关,它是一个固定的值。
4. 预积分量方差的计算
预积分误差定义:一段时间内IMU 构建的预积分量作为测量值,对两时刻之间的状态量进行约束
[
r
p
r
q
r
v
r
b
a
r
b
g
]
15
×
1
=
[
q
b
i
w
(
p
w
b
j
?
p
w
b
i
?
v
i
w
Δ
t
+
1
2
g
w
Δ
t
2
)
?
α
b
i
b
j
2
[
q
b
j
b
i
?
(
q
b
i
w
?
q
w
b
j
)
]
x
y
z
q
b
i
w
(
v
j
w
?
v
i
w
+
g
w
Δ
t
)
?
β
b
i
b
j
b
j
a
?
b
i
a
b
j
g
?
b
i
g
]
\left[\begin{array}{c} \mathbf{r}_{p} \\ \mathbf{r}_{q} \\ \mathbf{r}_{v} \\ \mathbf{r}_{b a} \\ \mathbf{r}_{b g} \end{array}\right]_{15 \times 1}=\left[\begin{array}{c} \mathbf{q}_{b_{i} w}\left(\mathbf{p}_{w b_{j}}-\mathbf{p}_{w b_{i}}-\mathbf{v}_{i}^{w} \Delta t+\frac{1}{2} \mathbf{g}^{w} \Delta t^{2}\right)-\alpha_{b_{i} b_{j}} \\ 2\left[\mathbf{q}_{b_{j} b_{i}} \otimes\left(\mathbf{q}_{b_{i} w} \otimes \mathbf{q}_{w b_{j}}\right)\right]_{x y z} \\ \mathbf{q}_{b_{i} w}\left(\mathbf{v}_{j}^{w}-\mathbf{v}_{i}^{w}+\mathbf{g}^{w} \Delta t\right)-\beta_{b_{i} b_{j}} \\ \mathbf{b}_{j}^{a}-\mathbf{b}_{i}^{a} \\ \mathbf{b}_{j}^{g}-\mathbf{b}_{i}^{g} \end{array}\right]
???????rp?rq?rv?rba?rbg?????????15×1?=???????qbi?w?(pwbj???pwbi???viw?Δt+21?gwΔt2)?αbi?bj??2[qbj?bi???(qbi?w??qwbj??)]xyz?qbi?w?(vjw??viw?+gwΔt)?βbi?bj??bja??bia?bjg??big????????? 上面误差中位移,速度,偏置都是直接相减 得到。第二顶是关于四元数的旋转误差 ,其中 [.] xyz 表示只取四元数的虚部
(
x
,
y
,
z
)
(x, y, z)
(x,y,z) 组成的三维向量。
5. 预积分离散方法
关于预积分的计算,前面提到过有欧拉法与中值法。这里使用mid-point 方法 ,即两个相邻时刻k 到
k
+
1
k+1
k+1 的位姿是用两个时刻的测量值
a
,
w
a, w
a,w 的平均值来计算:
ω
=
1
2
(
(
ω
b
k
?
b
k
g
)
+
(
ω
b
k
+
1
?
b
k
g
)
)
q
b
i
b
k
+
1
=
q
b
i
b
k
?
[
1
1
2
ω
δ
t
]
a
=
1
2
(
q
b
i
b
k
(
a
b
k
?
b
k
a
)
+
q
b
i
b
k
+
1
(
a
b
k
+
1
?
b
k
a
)
)
α
b
i
b
k
+
1
=
α
b
i
b
k
+
β
b
i
b
k
δ
t
+
1
2
a
δ
t
2
β
b
i
b
k
+
1
=
β
b
i
b
k
+
a
δ
t
b
k
+
1
a
=
b
k
a
+
n
b
k
a
δ
t
b
k
+
1
g
=
b
k
g
+
n
b
k
g
δ
t
\begin{aligned} \boldsymbol{\omega} &=\frac{1}{2}\left(\left(\boldsymbol{\omega}^{b_{k}}-\mathbf{b}_{k}^{g}\right)+\left(\boldsymbol{\omega}^{b_{k+1}}-\mathbf{b}_{k}^{g}\right)\right) \\ \mathbf{q}_{b_{i} b_{k+1}} &=\mathbf{q}_{b_{i} b_{k}} \otimes\left[\begin{array}{c} 1 \\ \frac{1}{2} \boldsymbol{\omega} \delta t \end{array}\right] \\ \mathbf{a} &=\frac{1}{2}\left(\mathbf{q}_{b_{i} b_{k}}\left(\mathbf{a}^{b_{k}}-\mathbf{b}_{k}^{a}\right)+\mathbf{q}_{b_{i} b_{k+1}}\left(\mathbf{a}^{b_{k+1}}-\mathbf{b}_{k}^{a}\right)\right) \\ \alpha_{b_{i} b_{k+1}} &=\alpha_{b_{i} b_{k}}+\beta_{b_{i} b_{k}} \delta t+\frac{1}{2} \mathbf{a} \delta t^{2} \\ \beta_{b_{i} b_{k+1}} &=\beta_{b_{i} b_{k}}+\mathbf{a} \delta t \\ \mathbf{b}_{k+1}^{a} &=\mathbf{b}_{k}^{a}+\mathbf{n}_{\mathbf{b}_{k}^{a} \delta t} \\ \mathbf{b}_{k+1}^{g} &=\mathbf{b}_{k}^{g}+\mathbf{n}_{\mathbf{b}_{k}^{g} \delta t} \end{aligned}
ωqbi?bk+1??aαbi?bk+1??βbi?bk+1??bk+1a?bk+1g??=21?((ωbk??bkg?)+(ωbk+1??bkg?))=qbi?bk???[121?ωδt?]=21?(qbi?bk??(abk??bka?)+qbi?bk+1??(abk+1??bka?))=αbi?bk??+βbi?bk??δt+21?aδt2=βbi?bk??+aδt=bka?+nbka?δt?=bkg?+nbkg?δt??
6. 预积分量的方差
疑问?: 一个 IMU 数据作为测量值的噪声方差我们能够标定。现在,一段时间内多个 IMU 数据积分形成的预积分量的方差 呢? Covariance Propagation(协方差传播) 举例说明: 已知一个变量
y
=
A
x
,
x
∈
N
(
0
,
Σ
x
)
y=A x, x \in N\left(0, \Sigma_{x}\right)
y=Ax,x∈N(0,Σx?), 则有
Σ
y
=
A
Σ
x
A
?
\Sigma_{y}=A \Sigma_{x} A^{\top}
Σy?=AΣx?A?
Σ
y
=
E
(
(
A
x
)
(
A
x
)
?
)
=
E
(
A
x
x
?
A
?
)
=
A
Σ
x
A
?
\begin{aligned} \Sigma_{y}=& E\left((A x)(A x)^{\top}\right) \\ =& E\left(A x x^{\top} A^{\top}\right) \\ &=A \Sigma x A^{\top} \end{aligned}
Σy?==?E((Ax)(Ax)?)E(Axx?A?)=AΣxA?? 所以,要推导预积分量的协方差,我们需要知道 imu 噪声和预积分量之间的线性递推关系 。 假设已知了相邻时刻误差的线性传递方程:
η
i
k
=
F
k
?
1
η
i
k
?
1
+
G
k
?
1
n
k
?
1
\eta_{i k}=F_{k-1} \eta_{i k-1}+G_{k-1} n_{k-1}
ηik?=Fk?1?ηik?1?+Gk?1?nk?1? 比如: 状态量误差为
η
i
k
=
[
δ
θ
i
k
,
δ
v
i
k
,
δ
p
i
k
]
\eta_{i k}=\left[\delta \theta_{i k}, \delta v_{i k}, \delta p_{i k}\right]
ηik?=[δθik?,δvik?,δpik?], 测量唱声为
n
k
=
[
n
k
g
,
n
k
a
]
n_{k}=\left[n_{k}^{g}, n_{k}^{a}\right]
nk?=[nkg?,nka?] 。 误差的传递由两部分组成当前时刻的误差 传递给下一时刻,当前时刻测量噪声 传递给下一时刻。 一个有趣的例子: 综艺节目中常有传递信息的节目,前一个人根据上一个人的信息 + 自己的理解(测量)传递给下一个人,导致这个信息越传越错。 协方差矩阵可以通过递推计算得到:
Σ
i
k
=
F
k
?
1
Σ
i
k
?
1
F
k
?
1
?
+
G
k
?
1
Σ
n
G
k
?
1
T
\Sigma_{i k}=F_{k-1} \Sigma_{i k-1} F_{k-1}^{\top}+G_{k-1} \Sigma_{n} G_{k-1}^{T}
Σik?=Fk?1?Σik?1?Fk?1??+Gk?1?Σn?Gk?1T? 其中,
Σ
n
\Sigma_{n}
Σn? 是测量噪声的协方差矩阵, 方差从 i 时刻开始进行递推,
Σ
i
i
=
0
\Sigma_{i i}=0
Σii?=0 。
7. 状态误差线性递推公式的推导
通常对于状态量之间的递推关系是非线性的方程如
x
k
=
f
(
x
k
?
1
,
u
k
?
1
)
x_{k}=f\left(x_{k-1}, u_{k-1}\right)
xk?=f(xk?1?,uk?1?), 其中状态量为
x
,
u
x, u
x,u 为系统的输入量。
- 我们可以用两种方法来推导状态误差传递的线性递推关系:
- 一种是基于
一阶泰勒展开 的误差递推方程。 - 一种是基于
误差随时间变化 的递推方程。
7. 1基于泰勒展开的误差传送(应用于 EKF 的协方差 预测)
令状态量为
x
=
x
^
+
δ
x
x=\hat{x}+\delta x
x=x^+δx, 其中, 真值为
x
^
\hat{x}
x^, 误差为
δ
x
\delta x
δx 。另外, 输入量
u
u
u 的噪声为
n
n
n 。 非线性系统
x
k
=
f
(
x
k
?
1
,
u
k
?
1
)
x_{k}=f\left(x_{k-1}, u_{k-1}\right)
xk?=f(xk?1?,uk?1?) 的状态误差的线性递推关系如下:
δ
x
k
=
F
δ
x
k
?
1
+
G
n
k
?
1
\delta x_{k}=F \delta x_{k-1}+G n_{k-1}
δxk?=Fδxk?1?+Gnk?1? 其中,
F
\mathrm{F}
F 是状态量
x
k
x_{k}
xk? 对状态量
x
k
?
1
x_{k-1}
xk?1? 的雅克比矩阵,
G
\mathrm{G}
G 是状态量
x
k
x_{k}
xk? 对输入量
u
k
?
1
u_{k-1}
uk?1? 的雅克比矩阵。 证明:对非线性状态方程 进行一阶泰勒展开有:
x
k
=
f
(
x
k
?
1
,
u
k
?
1
)
x
^
k
+
δ
x
k
=
f
(
x
^
k
?
1
+
δ
x
k
?
1
,
u
^
k
?
1
+
n
k
?
1
)
x
^
k
+
δ
x
k
=
f
(
x
^
k
?
1
,
u
^
k
?
1
)
+
F
δ
x
k
?
1
+
G
n
k
?
1
\begin{gathered} x_{k}=f\left(x_{k-1}, u_{k-1}\right) \\ \hat{x}_{k}+\delta x_{k}=f\left(\hat{x}_{k-1}+\delta x_{k-1}, \hat{u}_{k-1}+n_{k-1}\right) \\ \hat{x}_{k}+\delta x_{k}=f\left(\hat{x}_{k-1}, \hat{u}_{k-1}\right)+F \delta x_{k-1}+G n_{k-1} \end{gathered}
xk?=f(xk?1?,uk?1?)x^k?+δxk?=f(x^k?1?+δxk?1?,u^k?1?+nk?1?)x^k?+δxk?=f(x^k?1?,u^k?1?)+Fδxk?1?+Gnk?1??
7.2 基于误差随时间变化的递推方程
如果我们能够推导状态误差随时间变化的导数关系, 比如:
δ
x
′
=
A
δ
x
+
B
n
\delta x^{\prime}=A \delta x+B n
δx′=Aδx+Bn 则误差状态的传递方程为:
δ
x
k
=
δ
x
k
?
1
+
δ
x
k
?
1
′
Δ
t
→
δ
x
k
=
(
I
+
A
Δ
t
)
δ
x
k
?
1
+
B
Δ
t
n
k
?
1
\begin{gathered} \delta x_{k}=\delta x_{k-1}+\delta x_{k-1}^{\prime} \Delta t \\ \rightarrow \delta x_{k}=(I+A \Delta t) \delta x_{k-1}+B \Delta t n_{k-1} \end{gathered}
δxk?=δxk?1?+δxk?1′?Δt→δxk?=(I+AΔt)δxk?1?+BΔtnk?1?? 这两种推导方式的可以看出有:
F
=
I
+
A
Δ
t
,
G
=
B
Δ
t
F=I+A \Delta t, G=B \Delta t
F=I+AΔt,G=BΔt
- 第一种方法不是很好么,为什么会随着去弄误差随时间的变化呢?
这是因为 VIO 系统中已经知道了状态的导数和状态之间的转移矩阵 。如:我们已经知道速度和状态量之间的关系:
v
˙
=
R
a
b
+
g
\dot{v}=R a^{b}+g
v˙=Rab+g 那我们就可以推导速度的误差和状态误差之间的关系,再每一项上都加上各自的误差 就有:
v
˙
+
δ
v
˙
=
R
(
I
+
[
δ
θ
]
×
)
(
a
b
+
δ
a
b
)
+
(
g
+
δ
g
)
δ
v
˙
=
R
δ
a
b
+
R
[
δ
θ
]
×
(
a
b
+
δ
a
b
)
+
δ
g
δ
v
˙
=
R
δ
a
b
?
R
[
a
b
]
×
δ
θ
+
δ
g
\begin{gathered} \dot{v}+\delta \dot{v}=R\left(I+[\delta \theta]_{\times}\right)\left(a^{b}+\delta a^{b}\right)+(g+\delta g) \\ \delta \dot{v}=R \delta a^{b}+R[\delta \theta]_{\times}\left(a^{b}+\delta a^{b}\right)+\delta g \\ \delta \dot{v}=R \delta a^{b}-R\left[a^{b}\right]_{\times} \delta \theta+\delta g \end{gathered}
v˙+δv˙=R(I+[δθ]×?)(ab+δab)+(g+δg)δv˙=Rδab+R[δθ]×?(ab+δab)+δgδv˙=Rδab?R[ab]×?δθ+δg? 由此就能依次类推,轻易写出整个
A
A
A 和
B
B
B 其他方程了。
8. 预积分的误差递推公式推导
首先回顾预积分的误差递推公式, 将测量噪声 也考虑进模型:
ω
=
1
2
(
(
ω
b
k
+
n
k
g
?
b
k
g
)
+
(
ω
b
k
+
1
+
n
k
+
1
g
?
b
k
g
)
)
q
b
i
b
k
+
1
=
q
b
i
b
k
?
[
1
1
2
ω
δ
t
]
a
=
1
2
(
q
b
i
b
k
(
a
b
k
+
n
k
a
?
b
k
a
)
+
q
b
i
b
k
+
1
(
a
b
k
+
1
+
n
k
+
1
a
?
b
k
a
)
)
α
b
i
b
k
+
1
=
α
b
i
b
k
+
β
b
i
b
k
δ
t
+
1
2
a
δ
t
2
β
b
i
b
k
+
1
=
β
b
i
b
k
+
a
δ
t
b
k
+
1
a
=
b
k
a
+
n
b
k
a
δ
t
b
k
+
1
g
=
b
k
g
+
n
b
k
g
δ
t
\begin{aligned} \boldsymbol{\omega} &=\frac{1}{2}\left(\left(\boldsymbol{\omega}^{b_{k}}+\mathbf{n}_{k}^{g}-\mathbf{b}_{k}^{g}\right)+\left(\boldsymbol{\omega}^{b_{k+1}}+\mathbf{n}_{k+1}^{g}-\mathbf{b}_{k}^{g}\right)\right) \\ \mathbf{q}_{b_{i} b_{k+1}} &=\mathbf{q}_{b_{i} b_{k}} \otimes\left[\begin{array}{c} 1 \\ \frac{1}{2} \boldsymbol{\omega} \delta t \end{array}\right] \\ \mathbf{a} &=\frac{1}{2}\left(\mathbf{q}_{b_{i} b_{k}}\left(\mathbf{a}^{b_{k}}+\mathbf{n}_{k}^{a}-\mathbf{b}_{k}^{a}\right)+\mathbf{q}_{b_{i} b_{k+1}}\left(\mathbf{a}^{b_{k+1}}+\mathbf{n}_{k+1}^{a}-\mathbf{b}_{k}^{a}\right)\right) \\ \boldsymbol{\alpha}_{b_{i} b_{k+1}} &=\boldsymbol{\alpha}_{b_{i} b_{k}}+\boldsymbol{\beta}_{b_{i} b_{k}} \delta t+\frac{1}{2} \mathbf{a} \delta t^{2} \\ \boldsymbol{\beta}_{b_{i} b_{k+1}} &=\boldsymbol{\beta}_{b_{i} b_{k}}+\mathbf{a} \delta t \\ \mathbf{b}_{k+1}^{a} &=\mathbf{b}_{k}^{a}+\mathbf{n}_{\mathbf{b}_{k}^{a}} \delta t \\ \mathbf{b}_{k+1}^{g} &=\mathbf{b}_{k}^{g}+\mathbf{n}_{\mathbf{b}_{k}^{g}} \delta t \end{aligned}
ωqbi?bk+1??aαbi?bk+1??βbi?bk+1??bk+1a?bk+1g??=21?((ωbk?+nkg??bkg?)+(ωbk+1?+nk+1g??bkg?))=qbi?bk???[121?ωδt?]=21?(qbi?bk??(abk?+nka??bka?)+qbi?bk+1??(abk+1?+nk+1a??bka?))=αbi?bk??+βbi?bk??δt+21?aδt2=βbi?bk??+aδt=bka?+nbka??δt=bkg?+nbkg??δt? 确定误差传递的状态量,噪声量 ,然后开始构建传递方程。 用前面一阶泰勒展开的推导方式
δ
x
k
=
F
δ
x
k
?
1
+
G
n
k
?
1
\delta x_{k}=F \delta x_{k-1}+G n_{k-1}
δxk?=Fδxk?1?+Gnk?1?, 我们?望能推导出如下的形式:
[
δ
α
b
k
+
1
b
k
+
1
′
δ
θ
b
k
+
1
b
k
+
1
′
δ
β
b
k
+
1
b
k
+
1
′
δ
b
k
+
1
a
δ
b
k
+
1
g
]
=
F
[
δ
α
b
k
b
k
′
δ
θ
b
k
b
k
′
δ
β
b
k
b
k
′
δ
b
k
a
δ
b
k
g
]
+
G
[
n
k
a
n
k
g
n
k
+
1
a
n
k
+
1
g
n
b
k
a
n
b
k
g
]
\left[\begin{array}{c} \delta \boldsymbol{\alpha}_{b_{k+1} b_{k+1}^{\prime}} \\ \delta \boldsymbol{\theta}_{b_{k+1} b_{k+1}^{\prime}} \\ \delta \boldsymbol{\beta}_{b_{k+1} b_{k+1}^{\prime}} \\ \delta \mathbf{b}_{k+1}^{a} \\ \delta \mathbf{b}_{k+1}^{g} \end{array}\right]=\mathbf{F}\left[\begin{array}{c} \delta \boldsymbol{\alpha}_{b_{k} b_{k}^{\prime}} \\ \delta \boldsymbol{\theta}_{b_{k} b_{k}^{\prime}} \\ \delta \boldsymbol{\beta}_{b_{k} b_{k}^{\prime}} \\ \delta \mathbf{b}_{k}^{a} \\ \delta \mathbf{b}_{k}^{g} \end{array}\right]+\mathbf{G}\left[\begin{array}{c} \mathbf{n}_{k}^{a} \\ \mathbf{n}_{k}^{g} \\ \mathbf{n}_{k+1}^{a} \\ \mathbf{n}_{k+1}^{g} \\ \mathbf{n}_{\mathbf{b}_{k}^{a}} \\ \mathbf{n}_{\mathbf{b}_{k}^{g}} \end{array}\right]
???????δαbk+1?bk+1′??δθbk+1?bk+1′??δβbk+1?bk+1′??δbk+1a?δbk+1g?????????=F???????δαbk?bk′??δθbk?bk′??δβbk?bk′??δbka?δbkg?????????+G?????????nka?nkg?nk+1a?nk+1g?nbka??nbkg????????????
F
,
G
\mathrm{F}, \mathrm{G}
F,G 为两个时刻间的协方差传递矩阵。 这里我们直接给出
F
,
G
F, G
F,G 的最终形式,后面会对部分项进行详细推导:
F
=
[
I
f
12
I
δ
t
?
1
4
(
q
b
i
b
k
+
q
b
i
b
k
+
1
)
δ
t
2
f
15
0
I
?
[
ω
]
×
0
0
?
I
δ
t
0
f
32
I
?
1
2
(
q
b
i
b
k
+
q
b
i
b
k
+
1
)
δ
t
f
35
0
0
0
I
0
0
0
0
0
I
]
\mathbf{F}=\left[\begin{array}{ccccc} \mathbf{I} & \mathbf{f}_{12} & \mathbf{I} \delta t & -\frac{1}{4}\left(\mathbf{q}_{b_{i} b_{k}}+\mathbf{q}_{b_{i} b_{k+1}}\right) \delta t^{2} & \mathbf{f}_{15} \\ \mathbf{0} & \mathbf{I}-[\boldsymbol{\omega}]_{\times} & \mathbf{0} & \mathbf{0} & -\mathbf{I} \delta t \\ \mathbf{0} & \mathbf{f}_{32} & \mathbf{I} & -\frac{1}{2}\left(\mathbf{q}_{b_{i} b_{k}}+\mathbf{q}_{b_{i} b_{k+1}}\right) \delta t & \mathbf{f}_{35} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I} \end{array}\right]
F=???????I0000?f12?I?[ω]×?f32?00?Iδt0I00??41?(qbi?bk??+qbi?bk+1??)δt20?21?(qbi?bk??+qbi?bk+1??)δtI0?f15??Iδtf35?0I????????
G
=
[
1
4
q
b
i
b
k
δ
t
2
g
12
1
4
q
b
i
b
k
+
1
δ
t
2
g
14
0
0
0
1
2
I
δ
t
0
1
2
I
δ
t
0
0
1
2
q
b
i
b
k
δ
t
g
32
1
2
q
b
i
b
k
+
1
δ
t
g
34
0
0
0
0
0
0
I
δ
t
0
0
0
0
0
0
I
δ
t
]
\mathbf{G}=\left[\begin{array}{cccccc} \frac{1}{4} \mathbf{q}_{b_{i} b_{k}} \delta t^{2} & \mathbf{g}_{12} & \frac{1}{4} \mathbf{q}_{b_{i} b_{k+1}} \delta t^{2} & \mathbf{g}_{14} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \frac{1}{2} \mathbf{I} \delta t & \mathbf{0} & \frac{1}{2} \mathbf{I} \delta t & \mathbf{0} & \mathbf{0} \\ \frac{1}{2} \mathbf{q}_{b_{i} b_{k}} \delta t & \mathbf{g}_{32} & \frac{1}{2} \mathbf{q}_{b_{i} b_{k+1}} \delta t & \mathbf{g}_{34} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I} \delta t & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I} \delta t \end{array}\right]
G=???????41?qbi?bk??δt2021?qbi?bk??δt00?g12?21?Iδtg32?00?41?qbi?bk+1??δt2021?qbi?bk+1??δt00?g14?21?Iδtg34?00?000Iδt0?0000Iδt????????
二、残差Jacobian的推导
1. 视觉重投影残差的Jacobian
视觉残差为: 对于第
i
\mathrm{i}
i 帧中的特征点, 它投影到第
j
\mathrm{j}
j 帧相机坐标系下的值为:
[
x
c
j
y
c
j
z
c
j
1
]
=
T
b
c
?
1
T
w
b
j
?
1
T
w
b
i
T
b
c
[
1
λ
u
c
i
1
λ
v
c
i
1
λ
1
]
\left[\begin{array}{c} x_{c_{j}} \\ y_{c_{j}} \\ z_{c_{j}} \\ 1 \end{array}\right]=\mathbf{T}_{b c}^{-1} \mathbf{T}_{w b_{j}}^{-1} \mathbf{T}_{w b_{i}} \mathbf{T}_{b c}\left[\begin{array}{c} \frac{1}{\lambda} u_{c_{i}} \\ \frac{1}{\lambda} v_{c_{i}} \\ \frac{1}{\lambda} \\ 1 \end{array}\right]
?????xcj??ycj??zcj??1??????=Tbc?1?Twbj??1?Twbi??Tbc??????λ1?uci??λ1?vci??λ1?1?????? 拆成三维坐标形式为:
f
c
j
=
[
x
c
j
y
c
j
z
c
j
]
=
R
b
c
?
R
w
b
j
?
R
w
b
i
R
b
c
1
λ
[
u
c
i
v
c
i
1
]
+
R
b
c
?
(
R
w
b
j
?
(
(
R
w
b
i
p
b
c
+
p
w
b
i
)
?
p
w
b
j
)
?
p
b
c
)
\begin{aligned} \mathbf{f}_{c_{j}}=\left[\begin{array}{c} x_{c_{j}} \\ y_{c_{j}} \\ z_{c_{j}} \end{array}\right] &=\mathbf{R}_{b c}^{\top} \mathbf{R}_{w b_{j}}^{\top} \mathbf{R}_{w b_{i}} \mathbf{R}_{b c} \frac{1}{\lambda}\left[\begin{array}{c} u_{c_{i}} \\ v_{c_{i}} \\ 1 \end{array}\right] \\ &+\mathbf{R}_{b c}^{\top}\left(\mathbf{R}_{w b_{j}}^{\top}\left(\left(\mathbf{R}_{w b_{i}} \mathbf{p}_{b c}+\mathbf{p}_{w b_{i}}\right)-\mathbf{p}_{w b_{j}}\right)-\mathbf{p}_{b c}\right) \end{aligned}
fcj??=???xcj??ycj??zcj???????=Rbc??Rwbj???Rwbi??Rbc?λ1????uci??vci??1????+Rbc??(Rwbj???((Rwbi??pbc?+pwbi??)?pwbj??)?pbc?)? 再推导各类 Jacobian 之前, 为了简化公式, 先定义如下变量:
f
b
i
=
R
b
c
f
c
i
+
p
b
c
f
w
=
R
w
b
i
f
b
i
+
p
w
b
i
f
b
j
=
R
w
b
j
?
(
f
w
?
p
w
b
j
)
\begin{aligned} \mathbf{f}_{b_{i}} &=\mathbf{R}_{b c} \mathbf{f}_{c_{i}}+\mathbf{p}_{b c} \\ \mathbf{f}_{w} &=\mathbf{R}_{w b_{i}} \mathbf{f}_{b_{i}}+\mathbf{p}_{w b_{i}} \\ \mathbf{f}_{b_{j}} &=\mathbf{R}_{w b_{j}}^{\top}\left(\mathbf{f}_{w}-\mathbf{p}_{w b_{j}}\right) \end{aligned}
fbi??fw?fbj???=Rbc?fci??+pbc?=Rwbi??fbi??+pwbi??=Rwbj???(fw??pwbj??)? Jacobian 为视觉误差对两个时刻的状态量,外参,以及逆深度求导:
J
=
[
?
r
c
?
[
δ
p
b
i
b
i
′
δ
θ
b
i
b
i
′
]
?
r
c
?
[
p
b
j
b
j
′
δ
θ
b
j
b
j
′
]
?
r
c
?
[
δ
p
c
c
′
δ
θ
c
c
′
]
?
r
c
?
δ
λ
]
\mathbf{J}=\left[\begin{array}{lll}\frac{\partial \mathbf{r}_{c}}{\partial\left[\begin{array}{l}\delta \mathbf{p}_{b_{i} b_{i}^{\prime}} \\ \delta \boldsymbol{\theta}_{b_{i} b_{i}^{\prime}}\end{array}\right]} \quad \frac{\partial \mathbf{r}_{c}}{\partial\left[\begin{array}{l}\mathbf{p}_{b_{j} b_{j}^{\prime}} \\ \delta \boldsymbol{\theta}_{b_{j} b_{j}^{\prime}}\end{array}\right]} \quad \frac{\partial \mathbf{r}_{c}}{\partial\left[\begin{array}{l}\delta \mathbf{p}_{c c^{\prime}} \\ \delta \boldsymbol{\theta}_{c c^{\prime}}\end{array}\right]} \quad \frac{\partial \mathbf{r}_{c}}{\partial \delta \lambda}\end{array}\right]
J=[?[δpbi?bi′??δθbi?bi′???]?rc???[pbj?bj′??δθbj?bj′???]?rc???[δpcc′?δθcc′??]?rc???δλ?rc???]
参考文献
VIO 中残差雅可比的推导
|