VINS中的滑窗优化策略,将滑出窗外的帧与滑窗内的帧的约束使用边缘化的形式保存为先验误差进行后续非线性优化,以保留约束信息。本文对具体的方案进行记录。
紧耦合VIO处理的两种误差分别为IMU误差与图像误差,采用LM、GN、Dogleg等方式迭代求解待优化参数增量,以求得最大似然估计,最小化非线性化误差。
迭代优化的核心即求解增量方程:
H
δ
x
=
b
H\delta x=b
Hδx=b 其中
H
=
J
T
J
,
b
=
J
T
e
,
e
=
(
J
T
)
+
b
H=J^TJ,b=J^Te,e=(J^T)^+b
H=JTJ,b=JTe,e=(JT)+b.此处的
(
)
+
()^+
()+为伪逆,
e
e
e为残差,
J
J
J为残差对于状态的Jacobian。
由于引入滑窗后,对于滑到窗口之外的帧不再优化其参数,但其与滑窗内的数据仍然存在约束,直接丢掉这些窗外帧会造成约束信息丢失,所以要将其封装成先验信息,加入到大的非线性优化问题中作为一部分误差,这一步即为边缘化。假设要边缘化的状态为
x
2
x_2
x2,要保留的状态为
x
1
x_1
x1,对于上述增量方程,变为:
[
H
11
H
12
H
21
H
22
]
[
δ
x
1
δ
x
2
]
=
[
b
1
b
2
]
\left[ \begin{matrix} H_{11} & H_{12} \\ H_{21} & H_{22} \end{matrix} \right] \left[ \begin{matrix} \delta x_1\\ \delta x_2\end{matrix} \right]=\left[ \begin{matrix} b_1\\ b_2\end{matrix} \right]
[H11H21H12H22][δx1δx2]=[b1b2] 边缘化的方法为舒尔补:
[
I
−
H
12
H
22
−
1
0
I
]
[
H
11
H
12
H
21
H
22
]
[
δ
x
1
δ
x
2
]
=
[
I
−
H
12
H
22
−
1
0
I
]
[
b
1
b
2
]
[
H
11
−
H
12
H
22
−
1
H
21
0
H
21
H
22
]
[
δ
x
1
δ
x
2
]
=
[
b
1
−
H
12
H
22
−
1
b
2
b
2
]
\left[ \begin{matrix} I & -H_{12}H_{22}^{-1} \\ 0 & I \end{matrix} \right]\left[ \begin{matrix} H_{11} & H_{12} \\ H_{21} & H_{22} \end{matrix} \right] \left[ \begin{matrix} \delta x_1\\ \delta x_2\end{matrix} \right]=\left[ \begin{matrix} I & -H_{12}H_{22}^{-1} \\ 0 & I \end{matrix} \right]\left[ \begin{matrix} b_1\\ b_2\end{matrix} \right] \\ \left[ \begin{matrix} H_{11}-H_{12}H_{22}^{-1}H_{21} & 0 \\ H_{21} & H_{22} \end{matrix} \right]\left[ \begin{matrix} \delta x_1\\ \delta x_2\end{matrix} \right]=\left[ \begin{matrix} b_{1}-H_{12}H_{22}^{-1}b_{2}\\ b_2\end{matrix} \right]
[I0−H12H22−1I][H11H21H12H22][δx1δx2]=[I0−H12H22−1I][b1b2][H11−H12H22−1H21H210H22][δx1δx2]=[b1−H12H22−1b2b2] 这样只需要求解第一个方程,即
(
H
11
−
H
12
H
22
−
1
H
21
)
δ
x
1
=
b
1
−
H
12
H
22
−
1
b
2
,
(H_{11}-H_{12}H_{22}^{-1}H_{21})\delta x_1= b_{1}-H_{12}H_{22}^{-1}b_{2},
(H11−H12H22−1H21)δx1=b1−H12H22−1b2, 这样既没有损失约束信息,又不需要求解
δ
x
2
\delta x_2
δx2,把
H
11
−
H
12
H
22
−
1
H
21
H_{11}-H_{12}H_{22}^{-1}H_{21}
H11−H12H22−1H21记为
H
0
∗
H_0^*
H0∗,
b
1
−
H
12
H
22
−
1
b
2
b_{1}-H_{12}H_{22}^{-1}b_{2}
b1−H12H22−1b2记为
b
0
∗
b_0^*
b0∗。
上式为边缘化后待求解的增量方程,我们要根据这个增量方程恢复出边缘化封装的先验误差的具体形式
e
p
e_p
ep
随着迭代的过程,状态量被不断更新。但在边缘化时,被边缘化的状态值固定,舒尔补时使用的雅克比即为当时泰勒展开使用该固定的值(x线性化点)求得的雅克比,都需要固定。在VINS中,线性化点处的参数值x0被保存为keep_block_data。此即为EFJ,更多解释详看滑窗优化、边缘化、舒尔补、FEJ及fill-in问题。但在上述方程中,由于
b
=
J
T
e
b=J^Te
b=JTe,其中的
e
e
e随着状态的更新而变化,因此
b
b
b也需随之变化。
状态的更新(此处
δ
x
\delta x
δx表示当前参数值相对
x
0
x_0
x0的变化):
x
′
=
E
x
p
(
δ
x
)
⊗
x
0
x'=Exp(\delta x)\otimes x_0
x′=Exp(δx)⊗x0 此时b的更新可以表示为
b
=
b
0
+
∂
b
∂
x
∣
x
0
δ
x
=
b
0
+
J
T
∂
e
∂
x
∣
x
0
δ
x
=
b
0
+
J
T
J
δ
x
=
b
0
+
H
δ
x
b=b_0+\frac{\partial b}{\partial x}|_{x_0}\delta x\\=b_0+J^T\frac{\partial e}{\partial x}|_{x_0}\delta x\\=b_0+J^TJ\delta x=b_0+H\delta x
b=b0+∂x∂b∣x0δx=b0+JT∂x∂e∣x0δx=b0+JTJδx=b0+Hδx 有
[
b
1
b
2
]
=
[
b
1
,
0
b
2
,
0
]
+
[
H
11
H
12
H
21
H
22
]
[
δ
x
1
δ
x
2
]
b
1
=
b
1
,
0
+
H
11
δ
x
1
+
H
12
δ
x
2
b
2
=
b
2
,
0
+
H
21
δ
x
1
+
H
22
δ
x
2
\left[ \begin{matrix} b_1\\b_2\end{matrix} \right]=\left[ \begin{matrix} b_{1,0}\\b_{2,0}\end{matrix} \right]+\left[ \begin{matrix} H_{11} & H_{12} \\ H_{21} & H_{22} \end{matrix} \right] \left[ \begin{matrix} \delta x_1\\ \delta x_2\end{matrix} \right] \\ b_1=b_{1,0}+H_{11} \delta x_1+ H_{12} \delta x_2\\ b_2=b_{2,0}+H_{21} \delta x_1+ H_{22} \delta x_2
[b1b2]=[b1,0b2,0]+[H11H21H12H22][δx1δx2]b1=b1,0+H11δx1+H12δx2b2=b2,0+H21δx1+H22δx2 联立
b
∗
=
b
1
−
H
12
H
22
−
1
b
2
b^*=b_{1}-H_{12}H_{22}^{-1}b_{2}
b∗=b1−H12H22−1b2 有
b
∗
=
b
1
−
H
12
H
22
−
1
b
2
=
(
b
1
,
0
+
H
11
δ
x
1
+
H
12
δ
x
2
)
−
H
12
H
22
−
1
(
b
2
,
0
+
H
21
δ
x
1
+
H
22
δ
x
2
)
=
b
1
,
0
−
H
12
H
22
−
1
b
2
,
0
+
(
H
11
−
H
12
H
22
−
1
H
21
)
δ
x
1
=
b
0
∗
+
H
0
∗
δ
x
1
b^*=b_{1}-H_{12}H_{22}^{-1}b_{2}\\=(b_{1,0}+H_{11} \delta x_1+ H_{12} \delta x_2)-H_{12}H_{22}^{-1}(b_{2,0}+H_{21} \delta x_1+ H_{22} \delta x_2)\\ =b_{1,0}-H_{12}H_{22}^{-1}b_{2,0}+(H_{11}-H_{12}H_{22}^{-1}H_{21})\delta x_1\\ =b_0^*+H_0^*\delta x_1
b∗=b1−H12H22−1b2=(b1,0+H11δx1+H12δx2)−H12H22−1(b2,0+H21δx1+H22δx2)=b1,0−H12H22−1b2,0+(H11−H12H22−1H21)δx1=b0∗+H0∗δx1 由上述推导,原增量方程变为
H
0
∗
δ
x
=
J
l
T
J
l
δ
x
=
b
∗
=
b
0
∗
+
H
0
∗
d
x
=
b
0
∗
+
J
l
T
J
l
d
x
=
J
l
T
(
J
l
T
)
+
b
0
∗
+
J
l
T
J
l
d
x
=
J
l
T
(
(
J
l
T
)
+
b
0
∗
+
J
l
d
x
)
H_0^*\delta x=J_l^TJ_l\delta x=b^*=b_0^*+H_0^*dx\\ =b_0^*+J_l^TJ_ldx=J_l^T(J_l^T)^+b_0^*+J_l^TJ_ldx=J_l^T((J_l^T)^+b_0^*+J_ldx)
H0∗δx=JlTJlδx=b∗=b0∗+H0∗dx=b0∗+JlTJldx=JlT(JlT)+b0∗+JlTJldx=JlT((JlT)+b0∗+Jldx) 其中
J
l
J_l
Jl为
H
0
∗
H_0^*
H0∗分解出的等价雅克比,
δ
x
\delta x
δx为每次参数更新增量,
d
x
dx
dx为当前参数值与线性化点
x
0
x_0
x0的差值。
因此,边缘化后等价的先验误差
e
p
=
(
J
l
T
)
+
b
∗
=
(
J
l
T
)
+
b
0
∗
+
J
l
d
x
e_p=(J_l^T)^+b^*\\=(J_l^T)^+b_0^*+J_ldx
ep=(JlT)+b∗=(JlT)+b0∗+Jldx至此,边缘化的步骤可以总结为:
- 舒尔补求解 H 0 ∗ = H 11 − H 12 H 22 − 1 H 21 , b 0 ∗ = b 1 , 0 − H 12 H 22 − 1 b 2 , 0 H_0^*=H_{11}-H_{12}H_{22}^{-1}H_{21},b_0^*=b_{1,0}-H_{12}H_{22}^{-1}b_{2,0} H0∗=H11−H12H22−1H21,b0∗=b1,0−H12H22−1b2,0
- 分解 H 0 ∗ H_0^* H0∗求解 J l 、 J l T 、 ( J l T ) + J_l、J_l^T、(J_l^T)^+ Jl、JlT、(JlT)+,求解 e 0 = ( J l T ) + b 0 ∗ e_0=(J_l^T)^+b_0^* e0=(JlT)+b0∗,保存 J l , e 0 J_l,e_0 Jl,e0
- 求解 δ x \delta x δx由 H 0 ∗ δ x = b 0 ∗ H_0^*\delta x=b_0^* H0∗δx=b0∗(ceres完成)
- 使用
δ
x
\delta x
δx更新
x
x
x,计算当前状态
x
x
x与
x
0
x_0
x0的差值
d
x
dx
dx,更新
e
e
e:
e
p
=
e
0
+
J
l
d
x
e_p=e_0+J_ldx
ep=e0+Jldx(costfunction定义residuals,
同时定义costfunction需要的jacobians,把 J l J_l Jl由localsize转为globalsize)
ceres自动更新 b ∗ b^* b∗: b ∗ = J l T e p b^*=J_l^Te_p b∗=JlTep: - 求解 δ x \delta x δx由 H 0 ∗ δ x = b ∗ H_0^*\delta x=b^* H0∗δx=b∗(ceres完成)
- 循环4,5
附录:
使用H分解J的方式如下:
H
=
U
S
U
T
=
J
T
J
H=USU^T=J^TJ
H=USUT=JTJ 对H进行奇异值分解,由于H为半正定实矩阵,S也即为特征值矩阵,这样
J
=
S
1
2
U
T
(
J
T
)
+
=
S
−
1
2
U
T
J=S^\frac12U^T\\(J^T)^+=S^{-\frac12}U^T
J=S21UT(JT)+=S−21UT