前言
预积分的推导过程比较多,所以这里只记录关键结论。
其实这些公式不太好记忆,因为预积分推导过程的想法来源很巧妙,无法看出物理意义。
预积分定义式(必须记住)
一切推导的来源:
最好记忆的旋转相对运动量:
△
R
i
j
≐
R
i
T
R
j
=
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
g
,
k
−
η
g
d
,
k
)
△
t
)
\triangle R_{ij} \doteq R_{i}^{T}R_{j} = \prod_{k=i}^{j-1}Exp((\tilde{\omega}_{k}-b_{g,k}-\eta_{gd,k})\triangle t)
△Rij≐RiTRj=k=i∏j−1Exp((ω~k−bg,k−ηgd,k)△t)
△
v
i
j
≐
R
i
T
(
v
j
−
v
i
−
g
△
t
i
j
)
=
∑
k
=
i
j
−
1
△
R
i
k
(
a
~
k
−
b
a
,
k
−
η
a
d
,
k
)
△
t
\triangle v_{ij} \doteq R_{i}^{T}(v_{j} - v_{i} -g\triangle t_{ij}) = \sum_{k=i}^{j-1}\triangle R_{ik}(\tilde{a}_{k}-b_{a,k}-\eta_{ad, k})\triangle t
△vij≐RiT(vj−vi−g△tij)=k=i∑j−1△Rik(a~k−ba,k−ηad,k)△t
△
p
i
j
≐
R
i
T
(
p
j
−
p
i
−
v
i
△
t
i
j
−
1
2
g
△
t
i
j
2
)
=
∑
k
=
i
j
−
1
[
△
v
i
k
△
t
+
1
2
△
R
i
k
(
a
~
k
−
b
a
,
k
−
η
a
d
,
k
)
△
t
2
]
\triangle p_{ij} \doteq R_{i}^{T} \left ( p_{j} - p_{i} - v_{i}\triangle t_{ij}-\frac{1}{2}g\triangle t_{ij}^{2} \right ) = \sum_{k=i}^{j-1} \left[ \triangle v_{ik} \triangle t +\frac{1}{2}\triangle R_{ik} \left( \tilde{a}_{k} - b_{a,k} -\eta_{ad, k} \right) \triangle t^{2} \right]
△pij≐RiT(pj−pi−vi△tij−21g△tij2)=k=i∑j−1[△vik△t+21△Rik(a~k−ba,k−ηad,k)△t2]
η g d , k \eta_{gd,k} ηgd,k和 η a d , k \eta_{ad,k} ηad,k是离散化后的随机游走噪声(即建模为零均值白噪声高斯过程的测量噪声,因为也是角度或速度的随机游走所以才会这样描述)。
预积分的意义
IMU本身的数学模型中包括测量方程和噪声模型,再加上IMU运动学方程(本质上是微分方程组),可以推导出IMU在 △ t \triangle t △t时间段内的直接积分公式,在滤波方法中可以用直接积分作为预测,但是在优化方法中,我们并不希望优化过程随着IMU数据进行调用,这样做太浪费计算资源,我们更希望将IMU数据测量值组合在一起处理,所以累积的IMU观测就是所谓的IMU预积分。
直接积分的缺点是,它描述的过程和状态量有关,如果对 i i i时刻的状态进行优化,那么随后的所有IMU时刻的状态也会变化,整个积分要重新计算。
而通过预积分的定义可以看到,预积分定义最右侧公式与状态量无关(除了零偏,零偏发生变化时有单独的更新方法),所以当状态量发生变化时,也无需重新计算预积分量,更不必计算各个IMU时刻的状态量。此外预积分是累加或累乘的形式,给计算层面带来很大便利。
预积分的测量模型和噪声模型
推导过程假设了预积分过程中**零偏固定不变(所以下面公式中是 b g , i b_{g,i} bg,i或 b a , i b_{a,i} ba,i)**进行处理,实际零偏发生变化时预积分有单独的更新方式。
旋转部分
预积分旋转观测量(记住):
△
R
~
i
j
≐
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
g
,
i
)
△
t
)
\triangle \tilde{R}_{ij} \doteq \prod_{k=i}^{j-1}Exp((\tilde{\omega}_{k} - b_{g,i})\triangle t)
△R~ij≐k=i∏j−1Exp((ω~k−bg,i)△t)
旋转噪声项:
E
x
p
(
−
δ
ϕ
i
j
)
≐
∏
k
=
i
j
−
1
E
x
p
(
−
△
R
~
k
+
1
,
j
T
J
r
,
k
η
g
d
,
k
△
t
)
Exp(-\delta\phi_{ij}) \doteq \prod_{k=i}^{j-1}Exp(-\triangle\tilde{R}_{k+1,j}^{T}J_{r,k}\eta_{gd,k}\triangle t)
Exp(−δϕij)≐k=i∏j−1Exp(−△R~k+1,jTJr,kηgd,k△t)
J
r
,
k
J_{r,k}
Jr,k由BCH公式引入,和
ω
k
\omega_{k}
ωk、
b
g
,
k
b_{g,k}
bg,k有关(假设条件使得
b
g
,
k
=
b
g
,
i
b_{g,k} = b_{g,i}
bg,k=bg,i)。随机变量
δ
ϕ
i
j
\delta\phi_{ij}
δϕij只和随机变量
η
g
d
\eta_{gd}
ηgd有关,
η
g
d
\eta_{gd}
ηgd是均值为
0
0
0的白噪声,所以
δ
ϕ
i
j
\delta\phi_{ij}
δϕij均值为
0
0
0。
然后有:
△
R
i
j
=
△
R
~
i
j
E
x
p
(
−
δ
ϕ
i
j
)
\triangle R_{ij} = \triangle \tilde{R}_{ij} Exp(-\delta\phi_{ij})
△Rij=△R~ijExp(−δϕij)
旋转噪声项递推更新公式:
δ
ϕ
i
j
=
△
R
~
j
−
1
T
δ
ϕ
i
,
j
−
1
+
J
r
,
j
−
1
η
g
d
,
j
−
1
△
t
\delta\phi_{ij} = \triangle \tilde{R}_{j-1}^{T}\delta\phi_{i,j-1} + J_{r,j-1}\eta_{gd,j-1}\triangle t
δϕij=△R~j−1Tδϕi,j−1+Jr,j−1ηgd,j−1△t
旋转噪声项的协方差递推公式:
设
j
−
1
j-1
j−1时刻
δ
ϕ
i
,
j
−
1
\delta\phi_{i,j-1}
δϕi,j−1的协方差为
Σ
j
−
1
\Sigma_{j-1}
Σj−1,高斯白噪声
η
g
d
\eta_{gd}
ηgd的协方差为
Σ
η
g
d
\Sigma_{\eta_{gd}}
Σηgd(一般假设随机游走过程中高斯白噪声协方差保持不变,所以不区分
η
g
d
,
k
\eta_{gd,k}
ηgd,k)。
Σ
j
=
△
R
~
j
−
1
,
j
T
Σ
j
−
1
△
R
~
j
−
1
,
j
+
J
r
,
j
−
1
Σ
η
g
d
J
r
,
j
−
1
T
△
t
2
\Sigma_{j} = \triangle\tilde{R}_{j-1,j}^{T} \Sigma_{j-1} \triangle\tilde{R}_{j-1,j} + J_{r,j-1}\Sigma_{\eta_{gd}} J_{r,j-1}^{T} \triangle t^{2}
Σj=△R~j−1,jTΣj−1△R~j−1,j+Jr,j−1ΣηgdJr,j−1T△t2
速度部分
预积分速度观测量(记住):
△
v
~
i
j
≐
∑
k
=
i
j
−
1
△
R
~
i
k
(
a
~
k
−
b
a
,
i
)
△
t
\triangle \tilde{v}_{ij} \doteq \sum_{k=i}^{j-1}\triangle \tilde{R}_{ik}(\tilde{a}_{k}-b_{a,i})\triangle t
△v~ij≐k=i∑j−1△R~ik(a~k−ba,i)△t
速度噪声项:
−
δ
v
i
j
≐
∑
k
=
i
j
−
1
[
△
R
~
i
k
(
a
~
k
−
b
a
,
i
)
∧
δ
ϕ
i
k
△
t
−
△
R
~
i
k
η
a
d
,
k
△
t
]
-\delta v_{ij} \doteq \sum_{k=i}^{j-1}\left[ \triangle\tilde{R}_{ik}(\tilde{a}_{k}-b_{a,i})^{\wedge}\delta\phi_{ik}\triangle t - \triangle\tilde{R}_{ik}\eta_{ad,k}\triangle t \right]
−δvij≐k=i∑j−1[△R~ik(a~k−ba,i)∧δϕik△t−△R~ikηad,k△t]
然后有:
△
v
i
j
=
△
v
~
i
j
−
δ
v
i
j
\triangle v_{ij} = \triangle \tilde{v}_{ij} - \delta v_{ij}
△vij=△v~ij−δvij
速度噪声项递推公式:
δ
v
i
j
=
δ
v
i
,
j
−
1
−
△
R
~
i
,
j
−
1
(
a
~
j
−
1
−
b
a
,
i
)
∧
δ
ϕ
i
,
j
−
1
△
t
+
△
R
~
i
,
j
−
1
η
a
d
,
j
−
1
△
t
\delta v_{ij} = \delta v_{i,j-1} - \triangle\tilde{R}_{i,j-1}(\tilde{a}_{j-1}-b_{a,i})^{\wedge}\delta\phi_{i,j-1}\triangle t + \triangle\tilde{R}_{i,j-1}\eta_{ad,j-1}\triangle t
δvij=δvi,j−1−△R~i,j−1(a~j−1−ba,i)∧δϕi,j−1△t+△R~i,j−1ηad,j−1△t
位置(位移)部分
预积分位置(位移)观测量(记住):
△
p
~
i
j
≐
∑
k
=
i
j
−
1
[
(
△
v
~
i
k
△
t
)
+
1
2
△
R
~
i
k
(
a
~
k
−
b
a
,
i
)
△
t
2
]
\triangle \tilde{p}_{ij} \doteq \sum_{k=i}^{j-1}\left[ (\triangle\tilde{v}_{ik}\triangle t)+\frac{1}{2}\triangle \tilde{R}_{ik}(\tilde{a}_{k}-b_{a,i})\triangle t^{2}\right]
△p~ij≐k=i∑j−1[(△v~ik△t)+21△R~ik(a~k−ba,i)△t2]
位置(位移)噪声项:
−
δ
p
i
j
≐
∑
k
=
i
j
−
1
[
−
δ
v
i
k
△
t
+
1
2
△
R
~
i
k
(
a
~
k
−
b
a
,
i
)
∧
δ
ϕ
i
k
△
t
2
−
1
2
△
R
~
i
k
η
a
d
,
k
△
t
2
]
-\delta p_{ij} \doteq \sum_{k=i}^{j-1}\left[ -\delta v_{ik}\triangle t + \frac{1}{2}\triangle\tilde{R}_{ik}(\tilde{a}_{k}-b_{a,i})^{\wedge} \delta\phi_{ik}\triangle t^{2} - \frac{1}{2}\triangle\tilde{R}_{ik}\eta_{ad,k}\triangle t^{2} \right]
−δpij≐k=i∑j−1[−δvik△t+21△R~ik(a~k−ba,i)∧δϕik△t2−21△R~ikηad,k△t2]
然后有:
△
p
i
j
=
△
p
~
i
j
−
δ
p
i
j
\triangle p_{ij} = \triangle\tilde{p}_{ij} - \delta p_{ij}
△pij=△p~ij−δpij
位置(位移)噪声项递推公式:
δ
p
i
j
=
δ
p
i
,
j
−
1
+
δ
v
i
,
j
−
1
△
t
−
1
2
△
R
~
i
,
j
−
1
(
a
~
j
−
1
−
b
a
,
i
)
∧
δ
ϕ
i
,
j
−
1
△
t
2
+
1
2
△
R
~
i
,
j
−
1
η
a
d
,
j
−
1
△
t
2
\delta p_{ij} = \delta p_{i,j-1} + \delta v_{i,j-1}\triangle t - \frac{1}{2}\triangle\tilde{R}_{i,j-1}(\tilde{a}_{j-1}-b_{a,i})^{\wedge}\delta\phi_{i,j-1}\triangle t^{2} +\frac{1}{2} \triangle\tilde{R}_{i,j-1} \eta_{ad,j-1}\triangle t^{2}
δpij=δpi,j−1+δvi,j−1△t−21△R~i,j−1(a~j−1−ba,i)∧δϕi,j−1△t2+21△R~i,j−1ηad,j−1△t2
重要说明(重要!!!)
1、首先需要记住的公式均是预积分的定义式和观测量的定义式,其实二者形式几乎是相同的,只是观测量把定义式中的随机变量部分分离出去形成了噪声项。
2、说明一下本人对于预积分应用的理解,类比其他SLAM问题中的后端优化问题模型。可以看到预积分的公式可以写成“观测量=状态量推导的预测量+噪声”的形式,比如
△
R
~
i
j
=
△
E
x
p
(
−
δ
ϕ
i
j
)
\triangle \tilde{R}_{ij} = \triangle Exp(-\delta\phi_{ij})
△R~ij=△Exp(−δϕij),这个公式左侧的观测量由传感器观测值利用观测量定义式计算,右侧的预测量由本文第一部分中预积分定义式的中间公式计算,这样观察会发现,IMU预积分的数学模型和其他SLAM后端优化模型完全一样,求观测值-预测值的最小二乘问题其实就是求预积分的最大后验估计,而噪声项的协方差就是通过高斯分布的概率密度函数引入到优化问题中,指导着残差在求解优化问题中的比重。如果不采用完全准确的协方差,该模型将不会严格匹配预积分的最大后验估计,但是仍然具有最小二乘的意义,只不过残差在求解优化问题中的权重是自定义决定的。
预积分对于零偏的更新
前面的推导过程中,在观测量的定义部分假设了零偏不变,因此考虑零偏变化时观测量和零偏的关系。思路就是假设预积分观测和零偏之间是线性关系,利用泰勒展开保留一阶项(旋转部分的更新方式稍有不同)。重点在于雅可比矩阵的计算,最后还要把雅可比矩阵写成递推关系便于计算。
实际应用中,当优化后零偏发生变化时,利用递推公式先求雅可比矩阵,然后更新到预积分观测量上去。
公式太多了就不敲了(应该也不用记住)。
预积分残差对状态量的雅可比矩阵
这个公式也不敲了,需要的时候查书就可以了。
总结
在实际应用中,如果要进行IMU预积分,当一帧IMU信息到来时需要做:
1、计算三个预积分观测量
2、计算三个噪声项的协方差矩阵,作为后续图优化的信息矩阵
3、计算预积分观测量对零偏的雅可比矩阵,一共五个。
结束预积分计算后,可以把这些结果取出并在优化过程中使用。