参考文献
- https://blog.csdn.net/qq_26682225/article/details/72649858
- On-Manifold Preintegration for Real-TimeVisual-Inertial Odometry
- 【泡泡机器人原创专栏】IMU预积分总结与公式推导(一)~(四)
一、预积分的由来
来协调imu数据和图像之间的关系。通过重新参数化[^1]来把关键帧之间的IMU数据积分成相对运动的约束,避免初始条件变化 [^2]造成的重复积分, 预积分的值并不是真实的测量值.
预积分与因子图结合, 增量式因子图
二、相关理论
- 基于滤波的方法能够快速的推断当前的状态量,但是精度会随着线性化误差的累积而恶化。
- 基于优化的完全平滑方法精度高一些,但计算量很大。
- 固定滞后的平滑方法(滑动窗)是一个折中。
采用了structless model(边际化),这样可以不用延迟视觉信息的处理,同时可以多次线性化视觉测量模型。具体代码的实现整合在了gtsam4.0中
- 参加估计的pose number:
full smoothers: 所有的pose.
fixed-lag: 滑动窗口内的pose
filter: 最近的pose
- 代表不确定性的方法:
EKF: 协方差矩阵
信息滤波[^3]和smoothers: 信息矩阵
- 测量模型被线性化的次数:
standard EKF :1
smoothing appoach : 多次
fixed-lag和filter的方式相比对外点有更强的鲁棒性,fixed-lag在边际化的会使H矩阵变得稠密,为了保持稀疏性,需要丢掉一些测量点。此外两者都需要采用first estimate jacobian来线性化,来保持系统能观性不变,否则会引入错误的信息。
对于一个系统,Observability性质(能观性),决定了这个系统在进行状态估计时,哪些自由度是可以被估计出来的。并且其能观性是不受估计方法(Closed-form 方法、EKF、或者Nonlinear Optimization等等)改变的。能观性通过Observability Matrix(能观性矩阵)体现,系统Unobservable的状态维数是这个矩阵零空间的维数。
比如,单目纯视觉SLAM里,尺度和6DOF的绝对位姿——总共7DOF,都是无法被估计。可以通过固定某一帧来确定6个自由度,但是尺度无法确定。
再比如Visual-Inertial系统里,在运动激励充分(足够多轴向有足够大的加速度/角速度)的情况下,尺度、滚转/俯仰角是可以被估计的,只剩下绝对偏航角、绝对位置无法获得,也就是说对于Visual-Inertial系统在合适的运动模式下Unobservable的维数是4。磁力计可以确定偏航角,加上固定某一帧可以确定相对位置。
不同的线性化状态点计算得到的系数矩阵(雅克比矩阵)会导致Observability Matrix(能观性矩阵)的零空间维数出现不一致(inconsistency),导致自由度不一样出现偏差。为了修正这个问题,提出了First Estimate Jacobian。用相同的状态点进行线性化计算就不会有这个问题了,具体来说,位姿采用propagated值而非updated值,landmark使用第一次估计值。First Estimate就是这个直观含义,即采用estimation from the first time。
理解不够透彻:可参考文献:
[1] Kelly, J., & Sukhatme, G. S. (2011). Visual-inertial sensor fusion: Localization, mapping and sensor-to-sensor self-calibration. The International Journal of Robotics Research, 30(1), 56-79
[3] Strasdat, H., Montiel, J. M. M., & Davison, A. J. (2010). Scale drift-aware large scale monocular SLAM. Robotics: Science and Systems VI.
[4] Huang, G. P., Mourikis, A. I., & Roumeliotis, S. I. (2009). A first-estimates Jacobian EKF for improving SLAM consistency. In Experimental Robotics(pp. 373-382). Springer Berlin Heidelberg.作者:
[8] Hermann, R., & Krener, A. J. (1977). Nonlinear controllability and observability. IEEE Transactions on automatic control, 22(5), 728-740.
预积分理论的贡献
- 使得测量与绝对位姿解耦
- 利用重力使得俯仰和横滚变得可观
- 与视觉结合使得bias可观
三、初步
a.黎曼几何概念
SO3群定义:
S
O
(
3
)
=
{
R
∈
R
:
R
T
R
=
I
,
d
e
t
(
R
)
=
1
}
.
{ SO(3) =\lbrace {\bf R} \in {\Bbb R} :{\bf R^TR} = {\bf I}, det({\bf R}) = 1 \rbrace. }
SO(3)={R∈R:RTR=I,det(R)=1}.
他也代表着一个光滑的流形(manifold), 流形的切空间(tan)表示为李代数
S
O
(
3
)
\mathcal SO(3)
SO(3), 它与一个3维向量的3*3的反对称矩阵(向量加^表示)一致.
ω
∧
=
[
ω
1
ω
2
ω
3
]
∧
=
[
0
−
ω
3
ω
2
ω
3
0
−
ω
1
−
ω
2
ω
1
0
]
∈
S
O
(
3
)
.
{\omega ^\wedge = \begin{bmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{bmatrix}^\wedge = \begin{bmatrix} 0 && -\omega_3 && \omega_2 \\ \omega_3 && 0 && -\omega_1 \\ -\omega_2 && \omega_1 && 0 \end{bmatrix} \in \mathcal{SO} (3)}.
ω∧=⎣⎡ω1ω2ω3⎦⎤∧=⎣⎡0ω3−ω2−ω30ω1ω2−ω10⎦⎤∈SO(3).
他们之间的指数映射和对数映射如下:
E
x
p
:
R
3
→
S
O
(
3
)
;
ϕ
→
e
x
p
(
ϕ
∧
)
L
o
g
:
S
O
(
3
)
→
R
3
;
R
→
l
o
g
(
R
)
\begin{aligned} Exp \ &{:}\quad \Bbb{R}^3 &{\rightarrow} \quad & SO(3) &\ ; &\quad\bf{\phi} &\rightarrow \quad &exp(\phi^\wedge) \\ Log \ &{:}\quad SO(3) &{\rightarrow} \quad & \Bbb{R}^3 &\ ; &\quad\bf{R} &\rightarrow \quad &log( {\bf R} ) \end{aligned}
Exp Log :R3:SO(3)→→SO(3)R3 ; ;ϕR→→exp(ϕ∧)log(R)
E
x
p
Exp
Exp映射:
E
x
p
(
ϕ
)
=
I
+
s
i
n
(
∥
ϕ
∥
)
∥
ϕ
∥
ϕ
∧
+
1
−
c
o
s
(
∥
ϕ
∥
)
∥
ϕ
∥
2
(
ϕ
∧
)
2
Exp({\mathbf \phi})=I+\frac{sin(\|{\mathbf \phi}\|)}{\|{\mathbf \phi}\|}{\mathbf \phi}^\wedge+\frac{1-cos(\|{\mathbf \phi}\|)}{\|{\mathbf \phi}\|^2}({\bf \phi}^\wedge)^2
Exp(ϕ)=I+∥ϕ∥sin(∥ϕ∥)ϕ∧+∥ϕ∥21−cos(∥ϕ∥)(ϕ∧)2
L
o
g
Log
Log映射:
l
o
g
(
R
)
=
ϕ
⋅
(
R
−
R
T
)
2
s
i
n
(
ϕ
)
,
其
中
旋
转
角
ϕ
=
c
o
s
−
1
(
t
r
(
R
)
−
1
2
)
log(R)=\frac{\phi \cdot ({\bf R-R}^T)}{2sin(\phi)} ,\quad 其中 \quad旋转角 \phi=cos^{-1}(\frac{tr(R)-1}{2})
log(R)=2sin(ϕ)ϕ⋅(R−RT),其中旋转角ϕ=cos−1(2tr(R)−1)
李群的乘法并不符合正常指数的运算法则,即不等于李代数的和,因此当变化量是微小量时使用BCH公式近似模型得到一阶近似:
E
x
p
(
ϕ
+
δ
ϕ
)
≈
E
x
p
(
ϕ
)
E
x
p
(
J
r
(
ϕ
)
δ
ϕ
)
Exp({\bf\phi}\ +\ \delta{\bf\phi} ) \approx Exp({\bf\phi}) \, Exp(J_r(\phi)\delta{\bf\phi})
Exp(ϕ + δϕ)≈Exp(ϕ)Exp(Jr(ϕ)δϕ)
其中
J
r
(
ϕ
)
=
I
−
1
−
cos
(
∥
ϕ
∥
)
∥
ϕ
∥
2
ϕ
∧
+
∣
∣
ϕ
∣
∣
−
sin
(
∣
∣
ϕ
∣
∣
)
∣
∣
ϕ
∣
∣
3
(
ϕ
∧
)
2
J_r(\phi) = {\bf I} - \frac {1 - \cos(\parallel\phi\parallel)} {\parallel\phi\parallel ^2}\phi^\wedge \ + \ \frac {\mid\mid\phi\mid\mid -\sin(\mid\mid\phi\mid\mid )}{\mid\mid\phi\mid\mid ^3}(\phi^\wedge)^2
Jr(ϕ)=I−∥ϕ∥21−cos(∥ϕ∥)ϕ∧ + ∣∣ϕ∣∣3∣∣ϕ∣∣−sin(∣∣ϕ∣∣)(ϕ∧)2
同理反之,乘法变加法也成立。
图中,
J
r
(
ϕ
)
J_r(\phi)
Jr(ϕ)将切线空间中的累加增量与左侧的乘法增量相关联。
R
E
x
p
(
ϕ
)
R
T
=
e
x
p
(
R
ϕ
∧
R
T
)
=
E
x
p
(
R
ϕ
)
R \ Exp({\bf\phi}) \ R^T\ = \ exp(R{\bf\phi}^\wedge R^T) \ = \ Exp(R{\bf\phi})
R Exp(ϕ) RT = exp(Rϕ∧RT) = Exp(Rϕ)
⟺ E x p ( ϕ ) R = R E x p ( R T ϕ ) . \iff Exp({\bf\phi}) \ R \ = \ R Exp(R^T{\bf\phi}). ⟺Exp(ϕ) R = RExp(RTϕ).
注意: 这里的 E x p Exp Exp和 e x p exp exp是有区别的,大写的不需要^这个。
b.SO(3)中不确定性的描述
由前面知识可知:
R
~
=
R
E
x
p
(
ϵ
)
,
ϵ
∼
N
(
0
,
Σ
)
\tilde{R} = R\ Exp(\epsilon), \qquad \epsilon \sim N(0,\ \Sigma)
R~=R Exp(ϵ),ϵ∼N(0, Σ)
其中,
R
R
R是不包含噪声的,
ϵ
\epsilon
ϵ是一个正态分布的小扰动.
我们对高斯扰动(3维高斯公式)进行积分:
∫
R
3
p
(
ϵ
)
d
ϵ
=
∫
R
3
α
e
−
1
2
∥
ϵ
∥
Σ
2
d
ϵ
=
1
,
\int_{\R^3} p(\epsilon)\ d\epsilon = \int_{\R^3} \alpha e^{-\frac{1}{2}\parallel\epsilon\parallel^2_\Sigma }d \epsilon=1,
∫R3p(ϵ) dϵ=∫R3αe−21∥ϵ∥Σ2dϵ=1,
其中
α
=
1
/
(
2
π
)
3
d
e
t
(
Σ
)
\alpha = 1/ \sqrt{(2\pi)^3 det(\Sigma)}
α=1/(2π)3det(Σ),
∥
ϵ
∥
Σ
2
=
ϵ
T
Σ
−
1
ϵ
\parallel \epsilon\parallel^2_\Sigma = \epsilon^T\Sigma^{-1}\epsilon
∥ϵ∥Σ2=ϵTΣ−1ϵ是马氏距离的平方.
当
∥
ϵ
∥
<
π
\parallel\epsilon\parallel < \pi
∥ϵ∥<π时, 根据式(8)可以推出
ϵ
=
L
o
g
(
R
−
1
R
~
)
\epsilon = Log(R^{-1}\tilde{R})
ϵ=Log(R−1R~)代入式(9)
∫
S
O
(
3
)
β
(
R
~
)
e
−
1
2
∥
L
o
g
(
R
−
1
R
~
)
∥
Σ
2
d
R
~
=
1
,
\int_{SO(3)} \beta(\tilde{R}) \ e^{-\frac{1}{2}\parallel Log(R^{-1}\tilde{R})\parallel^2_\Sigma }d \tilde{R} = 1,
∫SO(3)β(R~) e−21∥Log(R−1R~)∥Σ2dR~=1,
β
(
R
~
)
\beta(\tilde{R})
β(R~)是一个归一化因子, 他等于
β
(
R
~
)
=
α
/
∣
d
e
t
(
J
(
R
~
)
)
∣
\beta(\tilde{R})=\alpha/\mid det(J(\tilde{R}))\mid
β(R~)=α/∣det(J(R~))∣, 由于积分变量的变化
J
(
R
~
)
≐
J
r
(
L
o
g
(
R
−
1
R
~
)
)
J(\tilde{R})\doteq J_r(Log(R^{-1}\tilde{R}))
J(R~)≐Jr(Log(R−1R~))
关于积分变量的推导:
δ
R
=
L
o
g
(
E
x
p
(
ϵ
+
δ
ϵ
)
E
x
p
(
ϵ
)
)
=
L
o
g
(
E
x
p
(
ϵ
)
E
x
p
(
J
r
δ
ϵ
)
E
x
p
(
ϵ
)
)
=
J
r
δ
ϵ
\begin{aligned} \delta R&= Log( \frac {Exp(\epsilon+\delta\epsilon)}{Exp(\epsilon)}) \\ &=Log(\frac {Exp(\epsilon) Exp(J_r\delta\epsilon)}{Exp(\epsilon)}) \\ &=J_r\delta\epsilon \end{aligned}
δR=Log(Exp(ϵ)Exp(ϵ+δϵ))=Log(Exp(ϵ)Exp(ϵ)Exp(Jrδϵ))=Jrδϵ
在SO(3)上的高斯扰动:
p
(
R
~
)
=
β
(
R
~
)
e
−
1
2
∥
L
o
g
(
R
−
1
R
~
)
∥
Σ
2
p(\tilde{R}) = \beta(\tilde{R}) \ e^{-\frac{1}{2}\parallel Log(R^{-1}\tilde{R})\parallel^2_\Sigma }
p(R~)=β(R~) e−21∥Log(R−1R~)∥Σ2
ϵ
\epsilon
ϵ是小的协方差矩阵则,
β
≈
α
\beta\approx\alpha
β≈α, 即
J
(
R
~
)
=
1
J(\tilde{R})=1
J(R~)=1, 并且
R
~
\tilde{R}
R~接近R, 求其最大似然(log):
L
(
R
)
=
1
2
∥
L
o
g
(
R
−
1
R
~
)
∥
Σ
2
+
c
o
n
s
t
=
1
2
∥
L
o
g
(
R
~
−
1
R
)
∥
Σ
2
+
c
o
n
s
t
{\mathcal L}(R) = \frac{1}{2}\parallel Log(R^{-1}\tilde{R})\parallel^2_\Sigma + const = \frac{1}{2}\parallel Log(\tilde{R}^{-1}{R})\parallel^2_\Sigma + const
L(R)=21∥Log(R−1R~)∥Σ2+const=21∥Log(R~−1R)∥Σ2+const
其中const是
L
o
g
(
β
)
Log(\beta)
Log(β).
他的几何意义是SO(3)上的距离的平方, 不确定性权重为协方差的逆 Σ − 1 \Sigma^{-1} Σ−1.
c.流形上的高斯牛顿法
欧几里得空间的高斯牛顿法是迭代优化目标函数的二次近似(通常非凸),求解中通常简化为一系列的线性方程。
当变量属于流形时
min
x
∈
M
f
(
x
)
\min_{x\in\mathcal{M}}\ f(x)
x∈Mmin f(x)
为了简化,我们只考虑一个变量时。
这里我们不能直接近似为 x x x的二次函数。因为,
- 会使参数增多(9×3个)会导致欠定。
- 这么做会使结果不再是流形 M \mathcal M M上的了,不满足SO(3)定义。
通常的做法是定义一个映射
R
x
\mathcal{R}_x
Rx, 它是由李代数(tan 空间)上的
δ
x
\delta x
δx到
x
∈
M
x \in \mathcal{M}
x∈M流形上附近值的一个映射.如下:
min
x
∈
M
f
(
x
)
⇒
min
δ
x
∈
R
n
f
(
R
(
δ
x
)
)
.
\min_{x\in\mathcal{M}}\ f(x) \Rightarrow \min_{\delta x \in \R^n} f(\mathcal{R}(\delta x)).
x∈Mmin f(x)⇒δx∈Rnminf(R(δx)).
这种重新参数化叫做lifting.。这样就可以在欧几里得空间(也就是tan空间)进行优化,粗略地说当前估计定义在tan空间,而附近的增量是欧几里得空间。
在高斯牛顿框架下,求解当前估计的代价函数平方,得到最小时的李代数
δ
x
∗
\delta x^*
δx∗(tan空间),最终更新流形上的估计值
x
^
←
R
x
^
(
δ
x
∗
)
\hat x \gets {\mathcal R}_{\hat x}(\delta x^*)
x^←Rx^(δx∗)
这种方法是误差状态模型的基础和统一推广,通常是用于航空航天文献中的滤波方法.
本文中使用以下定义这种retraction(缩放?)
R
R
(
ϕ
)
=
R
E
x
p
(
δ
ϕ
)
,
δ
ϕ
∈
R
3
{\mathcal R}_R(\phi)=\rm R \ Exp(\delta \phi), \quad \delta\phi \in\R^3
RR(ϕ)=R Exp(δϕ),δϕ∈R3
R T ( δ ϕ , δ p ) = ( R E x p ( δ ϕ ) , p + R δ p ) , [ δ ϕ δ p ] ∈ R 6 {\mathcal R}_T(\delta\phi,\delta {\bf p})=({\rm R \ Exp(\delta\phi)}, {\bf p}+{\rm R} \delta {\bf p} ), \quad \begin{bmatrix}\delta \phi &\delta{\bf p}\end{bmatrix} \in \R^6 RT(δϕ,δp)=(R Exp(δϕ),p+Rδp),[δϕδp]∈R6
四、最大后验视觉惯导状态估计
IMU的坐标系与机体坐标系一致,相机与机体系之间的变换关系已知.
系统状态变量为姿态,位置,速度,偏差:
x
i
≐
[
R
i
,
p
i
,
v
i
,
b
i
]
{\bf x}_i \doteq [{\bf R}_i,{\bf p}_i,{\bf v}_i,{\bf b}_i]
xi≐[Ri,pi,vi,bi]
K
k
{\mathcal K}_k
Kk表示到时间k所有关键帧的序列,我们估计所有关键帧的状态:
K
k
≐
{
X
i
}
i
∈
K
k
{\mathcal K}_k\doteq \lbrace {\rm X}_i \rbrace_{i\in{\mathcal K}_k}
Kk≐{Xi}i∈Kk
测量的关键帧图像
C
i
{\mathcal C}_i
Ci包含许多路标点
l
l
l, 因此有观测值
z
i
l
{\bf z}_{il}
zil.
I
i
j
{\mathcal I}_{ij}
Iij表示两关键帧之间的IMU测量值. 因此测量可以表示为:
Z
k
≐
{
C
i
,
I
i
j
}
(
i
,
j
)
∈
K
k
{\mathcal Z}_k \doteq \lbrace{\mathcal C}_i ,{\mathcal I}_{ij} \rbrace _{(i, j)\in{\mathcal K}_k}
Zk≐{Ci,Iij}(i,j)∈Kk
X
k
{\mathcal X}_k
Xk的后验概率:
p
(
X
k
∣
Z
k
)
∝
p
(
X
0
)
p
(
Z
k
∣
X
k
)
=
(
a
)
p
(
X
0
)
∏
(
i
,
j
)
∈
K
k
p
(
C
i
,
I
i
j
∣
X
k
)
=
(
b
)
p
(
X
0
)
∏
(
i
,
j
)
∈
K
k
p
(
I
i
j
∣
x
i
,
x
j
)
∏
i
∈
K
k
∏
l
∈
C
i
p
(
z
i
l
∣
x
i
)
p({\mathcal X}_k | {\mathcal Z _k} )\propto p({\mathcal X}_0)p({\mathcal Z}_k | {\mathcal X }_k ) \overset {(a)}= p({\mathcal X}_0) \prod _{(i,j)\in{\mathcal K_k}}p({\mathcal C_i, I_{ij}|}{\mathcal X_k}) \\ \overset {(b)}=p({\mathcal X_0}) \prod _{(i,j)\in{\mathcal K_k}}p({\mathcal I_{ij}|x_i, x_j}) \prod_{i \in {\mathcal K_k}} \prod_{l \in {\mathcal C_i}}p(z_{il}|x_i)
p(Xk∣Zk)∝p(X0)p(Zk∣Xk)=(a)p(X0)(i,j)∈Kk∏p(Ci,Iij∣Xk)=(b)p(X0)(i,j)∈Kk∏p(Iij∣xi,xj)i∈Kk∏l∈Ci∏p(zil∣xi)
使用因子图来求解.
使用-log形式进行MAP估计, 可以写成残差像的和, 通俗说就是当前状态下预测值与测量值之间的不匹配度.
X
k
∗
≐
a
r
g
min
X
k
−
log
e
p
(
X
k
∣
Z
k
)
=
arg
min
X
k
∥
r
0
∥
Σ
0
2
+
∑
(
i
,
j
)
∈
K
k
∥
r
I
i
,
j
∥
Σ
i
,
j
2
+
∑
i
∈
K
k
∑
l
∈
C
i
∥
r
C
i
l
∥
Σ
C
2
\begin{aligned} {\mathcal X}_k^* &\doteq arg \min _{{\mathcal X}_k} -\log_e p({\mathcal X}_k|{\mathcal Z}_k) \\ &=\arg \min _{{\mathcal X}_k}\|{\bf r}_0\|^2_{\Sigma_0}+\sum_{(i, j)\in{\mathcal K}_k}\|{\bf r}_{{\mathcal I}_{i,j}}\|^2_{\Sigma_{i,j}}+\sum_{i \in{\mathcal K}_k} \sum_{l\in{\mathcal C_i}}\|{\bf r}_{{\mathcal C}_{il}}\|^2_{\Sigma_{\mathcal C}} \end{aligned}
Xk∗≐argXkmin−logep(Xk∣Zk)=argXkmin∥r0∥Σ02+(i,j)∈Kk∑∥rIi,j∥Σi,j2+i∈Kk∑l∈Ci∑∥rCil∥ΣC2
五、IMU模型和运动积分
陀螺仪测量模型:
B
ω
~
W
B
(
t
)
=
B
ω
W
B
(
t
)
+
b
g
(
t
)
+
η
g
(
t
)
{_B}\tilde{\bf \omega}_{WB}(t)={_B}{\omega}_{WB}(t)+{\bf b}^{g}(t)+{\bf \eta}^g(t)
Bω~WB(t)=BωWB(t)+bg(t)+ηg(t)
加速度计测量模型:
B
a
~
(
t
)
=
R
W
B
T
(
t
)
(
W
a
(
t
)
−
W
g
)
+
b
a
(
t
)
+
η
a
(
t
)
_B \tilde {\bf a}(t)={\rm R_{WB}^{T}}(t)(_{W}{\bf a}(t)-_{W}{\bf g})+{\bf b}^a(t)+{\bf \eta}^a(t)
Ba~(t)=RWBT(t)(Wa(t)−Wg)+ba(t)+ηa(t)
其中下标B表示机体坐标系, W表示世界坐标系(静止的), WB表示B系到W系转换.
b
{\bf b}
b为缓慢变化的bias;
η
\eta
η表示白噪声.
运动模型:
R
W
B
˙
=
R
W
B
B
ω
W
B
∧
,
W
v
˙
=
W
a
,
W
p
˙
=
W
v
\dot{R_{\rm WB}}={R_{\rm WB}} _{\rm B}{\bf \omega}_{\rm WB}^{\wedge}, \quad _{\rm W}\dot{\bf v}=_{\rm W}{\bf a}, \quad _{\rm W}\dot{\bf p}=_{\rm W}{\bf v}
RWB˙=RWBBωWB∧,Wv˙=Wa,Wp˙=Wv
在时间段
[
t
,
t
+
Δ
t
]
[t,t+\Delta t]
[t,t+Δt]上使用欧拉积分
R
W
B
(
t
+
Δ
t
)
=
R
W
B
(
t
)
E
x
p
(
B
ω
W
B
(
t
)
Δ
t
)
W
v
(
t
+
Δ
t
)
=
W
v
(
t
)
+
W
a
(
t
)
Δ
t
W
p
(
t
+
Δ
t
)
=
W
p
(
t
)
+
W
v
(
t
)
Δ
t
+
1
2
W
a
(
t
)
Δ
t
2
(28)
\begin{aligned} &{\rm R_{WB}}(t+\Delta t)={\rm R_{WB}}(t){\rm Exp}(_{\rm B}{\bf \omega}_{\rm WB}(t)\Delta t) \\ &_{\rm W}{\bf v}(t+\Delta t)=_{\rm W}{\bf v}(t)+_{\rm W}{\bf a}(t)\Delta t \\ &_{\rm W}{\bf p}(t+\Delta t)=_{\rm W}{\bf p}(t)+_{\rm W}{\bf v}(t)\Delta t+\frac{1}{2} {_{\rm W}{\bf a}(t)}\Delta t^2 \end{aligned} \tag {28}
RWB(t+Δt)=RWB(t)Exp(BωWB(t)Δt)Wv(t+Δt)=Wv(t)+Wa(t)ΔtWp(t+Δt)=Wp(t)+Wv(t)Δt+21Wa(t)Δt2(28)
将(25)(26)的测量模型公式代入得到:
R
(
t
+
Δ
t
)
=
R
(
t
)
E
x
p
(
(
ω
~
(
t
)
−
b
g
(
t
)
−
η
g
d
(
t
)
)
Δ
t
)
v
(
t
+
Δ
t
)
=
v
(
t
)
+
g
Δ
t
+
R
(
t
)
(
a
~
(
t
)
−
b
a
(
t
)
−
η
a
d
(
t
)
)
Δ
t
p
(
t
+
Δ
t
)
=
p
(
t
)
+
v
(
t
)
Δ
t
+
1
2
g
Δ
t
2
+
1
2
R
(
t
)
(
a
~
(
t
)
−
b
a
(
t
)
−
η
a
d
(
t
)
)
Δ
t
2
(29)
\begin{aligned} &{\rm R}(t+\Delta t)={\rm R}(t){\rm Exp}((\tilde{\omega}(t)-{\bf b}^{g}(t)-{\bf \eta}^{gd}(t))\Delta t) \\ &{\bf v}(t+\Delta t)={\bf v}(t)+{\bf g}\Delta t +{\rm R}(t)(\tilde {\bf a}(t)-{\bf b}^a(t)-{\bf \eta}^{ad}(t))\Delta t \\ &{\bf p}(t+\Delta t)={\bf p}(t)+{\bf v}(t)\Delta t+\frac{1}{2} {\bf g}\Delta t^2+\frac{1}{2}{\rm R}(t) (\tilde{\bf a}(t)-{\bf b}^a(t)-{\bf \eta}^{ad}(t))\Delta t^2 \end{aligned}\tag {29}
R(t+Δt)=R(t)Exp((ω~(t)−bg(t)−ηgd(t))Δt)v(t+Δt)=v(t)+gΔt+R(t)(a~(t)−ba(t)−ηad(t))Δtp(t+Δt)=p(t)+v(t)Δt+21gΔt2+21R(t)(a~(t)−ba(t)−ηad(t))Δt2(29)
这里简化了符号的表示. 公式中我们假设
R
(
t
)
{\rm R}(t)
R(t)是恒定不变的, 这样会导致产生误差, 但当IMU频率较高时, 这个影响会大大减小. 当使用的IMU的频率较小时, 需要使用更高阶的数值积分方法.
其中的白噪声离散形式协方差为原噪声协方差的除以时间间隔 Δ t \Delta t Δt.
六、流形上的IMU预积分
IMU预积分的初衷,是希望借鉴纯视觉SLAM中图优化的思想,将帧与帧之间IMU相对测量信息转换为约束节点(载体位姿)的边参与到优化框架中。
- IMU预积分理论最大的贡献是对这些IMU相对测量进行处理,使得它与绝对位姿解耦(或者只需要线性运算就可以进行校正),从而大大提高优化速度。
- 这种优化架构还使得加计测量中不受待见的重力变成一个有利条件——重力的存在将使整个系统对绝对姿态(指相对水平地理坐标系的俯仰角和横滚角,不包括真航向)可观。要知道纯视觉VO或者SLAM是完全无法得到绝对姿态的。
- 两种传感器融合可以利用冗余测量(例如两种方式都可以求取相对位姿)来抑制累积误差。
- IMU和视觉这两种不同源的测量,也使得IMU的bias可观,从而可以在优化中被有效估计。
- 纯单目视觉缺乏绝对尺度的问题,也可以由惯性信息的引入而得以解决。
将两关键帧 i i i和 j j j之间的综合成一个复合的测量值, 叫preintegrated IMU measurement, 它约束连续关键帧之间的运动.
两个关键帧之间的关系如下:
R
j
=
R
i
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
k
g
−
η
k
g
d
)
Δ
t
)
,
v
j
=
v
i
+
g
Δ
t
i
j
+
∑
k
=
i
j
−
1
R
k
(
a
~
k
−
b
k
a
−
η
k
a
d
)
Δ
t
p
j
=
p
i
+
∑
k
=
i
j
−
1
[
v
k
Δ
t
+
1
2
g
Δ
t
2
+
1
2
R
k
(
a
~
k
−
b
k
a
−
η
k
a
d
)
Δ
t
2
]
(30)
\begin{aligned} &{\bf R}_j={\bf R}_i \prod_{k=i}^{j-1}{\rm Exp((\tilde \omega_k-b_k^g-\eta_k^{gd})\Delta t)}, \\ &{\bf v}_j={\bf v}_i+{\bf g}\Delta t_{ij}+\sum_{k=i}^{j-1}{\bf R}_k(\tilde{\bf a}_k-{\bf b}_k^a-\eta_k^{ad})\Delta t \\ &{\bf p}_j={\bf p}_i+\sum_{k=i}^{j-1}[{\bf v}_k \Delta t+\frac{1}{2}{\bf g}\Delta t^2+\frac{1}{2}{\bf R}_k(\tilde{\bf a}_k-{\bf b}_k^a-\eta_k^{ad})\Delta t^2] \end{aligned}\tag {30}
Rj=Rik=i∏j−1Exp((ω~k−bkg−ηkgd)Δt),vj=vi+gΔtij+k=i∑j−1Rk(a~k−bka−ηkad)Δtpj=pi+k=i∑j−1[vkΔt+21gΔt2+21Rk(a~k−bka−ηkad)Δt2](30)
从该公式(30)可以看出当初始值改变, 如
R
i
{\bf R}_i
Ri改变会导致每一个公式都需要重新计算. 为了解决这个缺陷, 给出不受线性化点影响的增量表达式:
Δ
R
i
j
=
R
i
T
R
j
=
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
k
g
−
η
k
g
d
)
Δ
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
k
a
−
η
k
a
d
)
Δ
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
k
a
−
η
k
a
d
)
Δ
t
2
]
(31)
\begin{aligned} &\Delta {\bf R}_{ij}={\bf R}_i^T{\bf R}_j= \prod_{k=i}^{j-1}{\rm Exp((\tilde \omega_k-b_k^g-\eta_k^{gd})\Delta t)}, \\ &\Delta {\bf v}_{ij}={\bf R}_i^T({\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij})=\sum_{k=i}^{j-1}\Delta{\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a-\eta_k^{ad})\Delta t \\ &\Delta {\bf p}_{ij} ={\bf R}_i^T({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2)\\ & \qquad \ =\sum_{k=i}^{j-1}[\Delta{\bf v}_{ik} \Delta t+\frac{1}{2}\Delta{\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a-\eta_k^{ad})\Delta t^2] \end{aligned}\tag {31}
ΔRij=RiTRj=k=i∏j−1Exp((ω~k−bkg−ηkgd)Δt),Δvij=RiT(vj−vi−gΔtij)=k=i∑j−1ΔRik(a~k−bka−ηkad)ΔtΔpij=RiT(pj−pi−viΔtij−21gΔtij2) =k=i∑j−1[ΔvikΔt+21ΔRik(a~k−bka−ηkad)Δt2](31)
最后的式子应用了等差数列求和公式推导, 证明如下:
这里的速度和位置 Δ \Delta Δ值并不能代表真正的物理变化, 但却使得它与 i i i时刻状态和重力独立, 能够直接根据惯性测量值计算出来. 这些值仍是Bias的函数, 假设Bias在两关键帧之间保持不变.
A. Preintegrated IMU Measurements
(1) Δ R i j \Delta {\bf R}_{ij} ΔRij项
公式中包含的测量噪声使得MAP(最大后验估计)无法应用, 我们对公式进行变换, 把噪声独立出来, 使得log似然更容易获得. 假设
i
i
i时刻bias已知且不变. 使用公式(6)和公式(9)将连乘展开后, 从第二个开始使用公式(8)进行两两合并, 之后再进行两两合并, 直到前面只有第一个Exp的连乘, 后面只有第二个Exp的复合连乘)_进行变换:
Δ
R
i
j
≃
e
q
(
6
)
∏
k
=
i
j
−
1
[
E
x
p
(
(
ω
~
k
−
b
i
g
)
Δ
t
)
E
x
p
(
−
J
r
k
η
k
g
d
Δ
t
)
]
=
e
q
(
9
)
Δ
R
~
i
j
∏
k
=
i
j
−
1
E
x
p
(
−
Δ
R
~
k
+
1
j
T
J
r
k
η
k
g
d
Δ
t
)
≐
Δ
R
~
i
j
E
x
p
(
−
δ
ϕ
i
j
)
(32)
\begin{aligned} \Delta {\bf R}_{ij} &\overset{eq(6)} \simeq \prod_{k=i}^{j-1}\ [{\rm Exp((\tilde \omega_k-b_i^g)\Delta t)}{\rm Exp}(-{\rm J}_r^k \eta_k^{gd} \Delta t)] \\ &\overset{eq(9)} = \Delta \tilde {\bf R}_{ij}\prod_{k=i}^{j-1}{\rm Exp}(-\Delta \tilde {\bf R}_{k+1j}^T {\rm J}_r^k \eta_k^{gd} \Delta t) \\ &\doteq \Delta \tilde {\bf R}_{ij}{\rm Exp}(-\delta \phi_{ij}) \end{aligned}\tag {32}
ΔRij≃eq(6)k=i∏j−1 [Exp((ω~k−big)Δt)Exp(−JrkηkgdΔt)]=eq(9)ΔR~ijk=i∏j−1Exp(−ΔR~k+1jTJrkηkgdΔt)≐ΔR~ijExp(−δϕij)(32)
其中:
J
r
k
≐
J
r
k
(
(
ω
~
k
−
b
i
g
)
Δ
t
)
Δ
R
~
i
j
=
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
i
g
)
Δ
t
)
(33)
{\rm J}_r^k \doteq {\rm J}_r^k((\tilde \omega_k-b_i^g)\Delta t) \\ \Delta \tilde {\bf R}_{ij}=\prod_{k=i}^{j-1}\ {\rm Exp}((\tilde \omega_k-b_i^g)\Delta t) \tag {33}
Jrk≐Jrk((ω~k−big)Δt)ΔR~ij=k=i∏j−1 Exp((ω~k−big)Δt)(33)
δ
ϕ
i
j
\delta \phi_{ij}
δϕij定义为噪声.
总结:
Δ
R
i
j
≐
Δ
R
~
i
j
E
x
p
(
−
δ
ϕ
i
j
)
Δ
R
~
i
j
=
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
i
g
)
Δ
t
)
E
x
p
(
−
δ
ϕ
i
j
)
=
∏
k
=
i
j
−
1
E
x
p
(
−
Δ
R
~
k
+
1
j
T
J
r
k
η
k
g
d
Δ
t
)
(34)
\Delta {\bf R}_{ij}\doteq \Delta \tilde {\bf R}_{ij}{\rm Exp}(- \delta \phi_{ij}) \\ \Delta \tilde {\bf R}_{ij}=\prod_{k=i}^{j-1}\ {\rm Exp}((\tilde \omega_k-b_i^g)\Delta t) \\ {\rm Exp}(- \delta \phi_{ij}) =\prod_{k=i}^{j-1}{\rm Exp}(-\Delta \tilde {\bf R}_{k+1j}^T {\rm J}_r^k \eta_k^{gd} \Delta t) \tag {34}
ΔRij≐ΔR~ijExp(−δϕij)ΔR~ij=k=i∏j−1 Exp((ω~k−big)Δt)Exp(−δϕij)=k=i∏j−1Exp(−ΔR~k+1jTJrkηkgdΔt)(34)
(2) Δ v i j \Delta {\bf v}_{ij} Δvij项
将式(32)代入(31)的速度预积分公式中, 对Exp进行一阶近似, 并且忽略高次微小量.得到:
Δ
v
i
j
=
∑
k
=
i
j
−
1
Δ
R
~
i
k
E
x
p
(
−
δ
ϕ
i
k
)
(
a
~
k
−
b
k
a
−
η
k
a
d
)
Δ
t
=
∑
k
=
i
j
−
1
[
Δ
R
~
i
k
(
I
−
δ
ϕ
i
k
∧
)
(
a
~
k
−
b
k
a
)
Δ
t
−
Δ
R
~
i
k
η
k
a
d
Δ
t
]
=
∑
k
=
i
j
−
1
[
Δ
R
~
i
k
(
a
~
k
−
b
k
a
)
Δ
t
]
+
∑
k
=
i
j
−
1
[
Δ
R
~
i
k
(
a
~
k
−
b
k
a
)
∧
δ
ϕ
i
k
Δ
t
−
Δ
R
~
i
k
η
k
a
d
Δ
t
]
(35)
\begin{aligned} \Delta {\bf v}_{ij} &= \sum_{k=i}^{j-1}\Delta \tilde {\bf R}_{ik}{\rm Exp}(- \delta \phi_{ik})(\tilde{\bf a}_k-{\bf b}_k^a-\eta_k^{ad})\Delta t \\ &= \sum_{k=i}^{j-1}[ \Delta \tilde {\bf R}_{ik}({\bf I}-\delta \phi_{ik}^\wedge)(\tilde{\bf a}_k-{\bf b}_k^a)\Delta t - \Delta \tilde {\bf R}_{ik}\eta_k^{ad}\Delta t ] \\ &=\sum_{k=i}^{j-1}[ \Delta \tilde {\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a)\Delta t ]+\sum_{k=i}^{j-1}[ \Delta \tilde {\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a)^\wedge\delta \phi_{ik}\Delta t- \Delta \tilde {\bf R}_{ik}\eta_k^{ad}\Delta t ] \end{aligned} \tag {35}
Δvij=k=i∑j−1ΔR~ikExp(−δϕik)(a~k−bka−ηkad)Δt=k=i∑j−1[ΔR~ik(I−δϕik∧)(a~k−bka)Δt−ΔR~ikηkadΔt]=k=i∑j−1[ΔR~ik(a~k−bka)Δt]+k=i∑j−1[ΔR~ik(a~k−bka)∧δϕikΔt−ΔR~ikηkadΔt](35)
总结:
Δ
v
i
j
≜
Δ
v
~
i
j
−
δ
v
i
j
Δ
v
~
i
j
=
∑
k
=
i
j
−
1
[
Δ
R
~
i
k
(
a
~
k
−
b
k
a
)
Δ
t
]
δ
v
i
j
=
∑
k
=
i
j
−
1
[
Δ
R
~
i
k
η
k
a
d
Δ
t
−
Δ
R
~
i
k
(
a
~
k
−
b
k
a
)
∧
δ
ϕ
i
k
Δ
t
]
(36)
\Delta {\bf v}_{ij} \triangleq \Delta \tilde {\bf v}_{ij}-\delta {\bf v}_{ij} \\ \Delta \tilde {\bf v}_{ij}=\sum_{k=i}^{j-1}[ \Delta \tilde {\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a)\Delta t ] \\ \delta {\bf v}_{ij} = \sum_{k=i}^{j-1}[ \Delta \tilde {\bf R}_{ik}\eta_k^{ad}\Delta t -\Delta \tilde {\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a)^\wedge\delta \phi_{ik}\Delta t ] \tag {36}
Δvij≜Δv~ij−δvijΔv~ij=k=i∑j−1[ΔR~ik(a~k−bka)Δt]δvij=k=i∑j−1[ΔR~ikηkadΔt−ΔR~ik(a~k−bka)∧δϕikΔt](36)
(3) Δ p i j \Delta{\bf p}_{ij} Δpij项
同上(34)和(36)式代入(31)中得到:
Δ
p
i
j
=
∑
k
=
i
j
−
1
[
(
Δ
v
~
i
k
−
δ
v
i
k
)
Δ
t
+
1
2
Δ
R
~
i
k
(
I
−
δ
ϕ
i
k
∧
)
(
a
~
k
−
b
k
a
)
Δ
t
2
−
1
2
Δ
R
~
i
k
η
k
a
d
Δ
t
2
]
=
∑
k
=
i
j
−
1
[
Δ
v
~
i
k
Δ
t
+
1
2
Δ
R
~
i
k
(
a
~
k
−
b
k
a
)
Δ
t
2
]
+
∑
k
=
i
j
−
1
[
1
2
Δ
R
~
i
k
(
a
~
k
−
b
k
a
)
∧
δ
ϕ
i
k
Δ
t
2
−
1
2
Δ
R
~
i
k
η
k
a
d
Δ
t
2
−
δ
v
i
k
Δ
t
]
(37)
\begin{aligned} \Delta {\bf p}_{ij}&=\sum_{k=i}^{j-1}[ ( \Delta \tilde {\bf v}_{ik}-\delta {\bf v}_{ik})\Delta t + \frac{1}{2} \Delta \tilde {\bf R}_{ik}({\bf I}-\delta \phi_{ik}^\wedge)(\tilde{\bf a}_k-{\bf b}_k^a)\Delta t^2 -\frac{1}{2} \Delta \tilde {\bf R}_{ik}\eta_k^{ad}\Delta t^2 ]\\ &=\sum_{k=i}^{j-1}[ \Delta \tilde {\bf v}_{ik}\Delta t+\frac{1}{2}\Delta \tilde {\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a)\Delta t^2 ] \\ & \quad + \sum_{k=i}^{j-1}[\frac{1}{2}\Delta \tilde {\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a)^\wedge \delta\phi_{ik}\Delta t^2 -\frac{1}{2} \Delta \tilde {\bf R}_{ik}\eta_k^{ad}\Delta t^2-\delta {\bf v}_{ik}\Delta t] \end{aligned} \tag{37}
Δpij=k=i∑j−1[(Δv~ik−δvik)Δt+21ΔR~ik(I−δϕik∧)(a~k−bka)Δt2−21ΔR~ikηkadΔt2]=k=i∑j−1[Δv~ikΔt+21ΔR~ik(a~k−bka)Δt2]+k=i∑j−1[21ΔR~ik(a~k−bka)∧δϕikΔt2−21ΔR~ikηkadΔt2−δvikΔt](37)
总结:
Δ
p
i
j
≜
Δ
p
~
i
j
−
δ
p
i
j
Δ
p
~
i
j
=
∑
k
=
i
j
−
1
[
Δ
v
~
i
k
Δ
t
+
1
2
Δ
R
~
i
k
(
a
~
k
−
b
k
a
)
Δ
t
2
]
δ
p
i
j
=
∑
k
=
i
j
−
1
[
1
2
Δ
R
~
i
k
η
k
a
d
Δ
t
2
−
1
2
Δ
R
~
i
k
(
a
~
k
−
b
k
a
)
∧
δ
ϕ
i
k
Δ
t
2
+
δ
v
i
k
Δ
t
]
(38)
\Delta {\bf p}_{ij} \triangleq \Delta \tilde {\bf p}_{ij}-\delta {\bf p}_{ij} \\ \Delta \tilde {\bf p}_{ij} =\sum_{k=i}^{j-1}[ \Delta \tilde {\bf v}_{ik}\Delta t+\frac{1}{2}\Delta \tilde {\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a)\Delta t^2 ] \\ \delta {\bf p}_{ij} = \sum_{k=i}^{j-1}[\frac{1}{2} \Delta \tilde {\bf R}_{ik}\eta_k^{ad}\Delta t^2 -\frac{1}{2}\Delta \tilde {\bf R}_{ik}(\tilde{\bf a}_k-{\bf b}_k^a)^\wedge \delta\phi_{ik}\Delta t^2 + \delta {\bf v}_{ik}\Delta t] \tag{38}
Δpij≜Δp~ij−δpijΔp~ij=k=i∑j−1[Δv~ikΔt+21ΔR~ik(a~k−bka)Δt2]δpij=k=i∑j−1[21ΔR~ikηkadΔt2−21ΔR~ik(a~k−bka)∧δϕikΔt2+δvikΔt](38)
公式(34)(36)(38)代入公式(31)得到预积分测量模型,公式(39)左侧为预积分测量值(由IMU测量值和Bias的猜测估计值组成), 右侧为真值+预积分测量噪声.
Δ
R
~
i
j
=
R
i
T
R
j
E
x
p
(
δ
ϕ
i
j
)
Δ
v
~
i
j
=
R
i
T
(
v
j
−
v
i
−
g
Δ
t
i
j
)
+
δ
v
i
j
Δ
p
~
i
j
=
R
i
T
(
p
j
−
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
+
δ
p
i
j
(39)
\Delta \tilde {\bf R}_{ij} ={\bf R}_i^{\rm T}{\bf R}_j{\rm Exp}(\delta\phi_{ij}) \\ \Delta \tilde{\bf v}_{ij}={\bf R}_i^T({\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij})+\delta{\bf v_{ij}} \\ \Delta \tilde{\bf p}_{ij}={\bf R}_i^T({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2)+\delta{\bf p}_{ij} \tag{39}
ΔR~ij=RiTRjExp(δϕij)Δv~ij=RiT(vj−vi−gΔtij)+δvijΔp~ij=RiT(pj−pi−viΔtij−21gΔtij2)+δpij(39)
B. Noise Propagation
希望使噪声符合高斯分布, 定义噪声向量:
η
i
j
Δ
≜
[
δ
ϕ
i
j
T
,
δ
v
i
j
T
,
δ
p
i
j
T
]
T
∼
N
(
0
9
×
1
,
Σ
i
j
)
(40)
{\bf \eta}_{ij}^\Delta \triangleq[\delta\phi_{ij}^{\rm T}, \ \delta{\bf v}_{ij}^{\rm T}, \ \delta{\bf p}_{ij}^{\rm T}]^{\rm T} \sim {\mathcal N}(0_{9\times1},\Sigma_{ij}) \tag{40}
ηijΔ≜[δϕijT, δvijT, δpijT]T∼N(09×1,Σij)(40)
对旋转量预积分噪声
δ
ϕ
i
j
\delta \phi _{ij}
δϕij两端取-Log得到:
δ
ϕ
i
j
=
−
L
o
g
(
Π
k
=
i
j
−
1
E
x
p
(
−
Δ
R
~
k
+
1
j
T
J
r
k
η
k
g
d
Δ
t
)
)
(41)
\delta \phi_{ij}=-{\rm Log}(\Pi_{k=i}^{j-1}{\rm Exp}(-\Delta \tilde{\bf R}_{k+1j}^{\rm T}{\rm J}^k_r\eta^{gd}_k\Delta t)) \tag{41}
δϕij=−Log(Πk=ij−1Exp(−ΔR~k+1jTJrkηkgdΔt))(41)
其中的
η
k
g
d
\eta_k^{gd}
ηkgd为微小量, 利用公式
L
o
g
(
E
x
p
(
ϕ
)
E
x
p
(
δ
ϕ
)
)
=
ϕ
+
J
r
−
1
(
ϕ
)
δ
ϕ
(42)
\rm Log( Exp(\phi)\ Exp(\delta \phi))=\phi+J_r^{-1}(\phi) \delta \phi \tag{42}
Log(Exp(ϕ) Exp(δϕ))=ϕ+Jr−1(ϕ)δϕ(42)
得到
δ
ϕ
i
j
≃
Σ
k
=
i
j
−
1
Δ
R
~
k
+
1
j
T
J
r
k
η
k
g
d
Δ
t
(43)
\delta \phi_{ij} \simeq \Sigma_{k=i}^{j-1}\Delta\tilde{\bf R}_{k+1j}^{\rm T}{\rm J}^k_r\boldsymbol{\eta}^{gd}_k\Delta t \tag{43}
δϕij≃Σk=ij−1ΔR~k+1jTJrkηkgdΔt(43)
证明:
由公式(43)可知, δ ϕ i j \delta \phi _{ij} δϕij是 η k g d \boldsymbol{\eta}^{gd}_k ηkgd线性组合, 因此也是零均值的高斯分布.
对于 δ v i j , δ p i j \delta{\bf v}_{ij}, \ \delta{\bf p}_{ij} δvij, δpij, 根据(36)(38)的增量项公式可知它们是 δ ϕ i j \delta \phi _{ij} δϕij和 η k a d \boldsymbol{\eta}^{ad}_k ηkad的线性组合, 因此也符合零均值的高斯分布.
B1. Iterative Noise Propagation
我们可以通过递推的形式求解
η
i
j
Δ
\boldsymbol{\eta}_{ij}^{\Delta}
ηijΔ和其协方差
Σ
i
j
\boldsymbol{\Sigma}_{ij}
Σij.
δ
ϕ
i
j
=
Δ
R
~
j
j
−
1
δ
ϕ
i
j
−
1
+
J
r
j
−
1
η
j
−
1
g
d
Δ
t
δ
v
i
j
=
δ
v
i
j
−
1
+
Δ
R
~
i
j
−
1
η
j
−
1
a
d
Δ
t
−
Δ
R
~
i
j
−
1
(
a
~
j
−
1
−
b
i
a
)
∧
δ
ϕ
i
j
−
1
Δ
t
δ
p
i
j
=
δ
p
i
j
−
1
+
δ
v
i
j
−
1
Δ
t
+
1
2
Δ
R
~
i
j
−
1
η
j
−
1
a
d
Δ
t
2
−
1
2
Δ
R
~
i
j
−
1
(
a
~
j
−
1
−
b
i
a
)
∧
δ
ϕ
i
j
−
1
Δ
t
2
(44)
\begin{aligned} \delta \phi_{ij}&=\Delta \tilde{\bf R}_{jj-1}\delta\phi_{ij-1}+{\bf J}_r^{j-1}\boldsymbol{\eta}_{j-1}^{gd}\Delta t \\ \delta{\bf v}_{ij}&= \delta {\bf v}_{ij-1}+ \Delta \tilde {\bf R}_{ij-1}\eta_{j-1}^{ad}\Delta t -\Delta \tilde {\bf R}_{ij-1}(\tilde{\bf a}_{j-1}-{\bf b}_i^a)^\wedge\delta \phi_{ij-1}\Delta t \\ \delta{\bf p}_{ij}&=\delta{\bf p}_{ij-1}+\delta {\bf v}_{ij-1}\Delta t+\frac{1}{2} \Delta \tilde {\bf R}_{ij-1}\boldsymbol{\eta}_{j-1}^{ad}\Delta t^2 -\frac{1}{2}\Delta \tilde {\bf R}_{ij-1}(\tilde{\bf a}_{j-1}-{\bf b}_i^a)^\wedge \delta\phi_{ij-1}\Delta t^2 \end{aligned} \tag{44}
δϕijδvijδpij=ΔR~jj−1δϕij−1+Jrj−1ηj−1gdΔt=δvij−1+ΔR~ij−1ηj−1adΔt−ΔR~ij−1(a~j−1−bia)∧δϕij−1Δt=δpij−1+δvij−1Δt+21ΔR~ij−1ηj−1adΔt2−21ΔR~ij−1(a~j−1−bia)∧δϕij−1Δt2(44)
其中第一个公式证明如下, 后两个公式即拆出一项即可.
定义IMU的测量噪声向量为
η
k
d
=
[
η
k
g
d
,
η
k
a
d
]
(45)
\boldsymbol{\eta}^d_{k}=[\boldsymbol{\eta}^{gd}_{k}, \ \boldsymbol{\eta}^{ad}_{k}] \tag{45}
ηkd=[ηkgd, ηkad](45)
将式(44)写成矩阵形式:
η
i
j
Δ
=
A
j
−
1
η
i
j
−
1
Δ
+
B
j
−
1
η
j
−
1
d
(46)
\boldsymbol{\eta}^\Delta_{ij}={\bf A}_{j-1}\boldsymbol{\eta}^\Delta_{ij-1}+{\bf B}_{j-1}\boldsymbol{\eta}^d_{j-1} \tag{46}
ηijΔ=Aj−1ηij−1Δ+Bj−1ηj−1d(46)
其中
A
j
−
1
=
[
Δ
R
~
j
j
−
1
0
0
−
Δ
R
~
i
j
−
1
(
a
~
j
−
1
−
b
i
a
)
∧
Δ
t
I
0
−
1
2
Δ
R
~
i
j
−
1
(
a
~
j
−
1
−
b
i
a
)
∧
Δ
t
2
Δ
t
I
I
]
(47)
{\bf A}_{j-1}=\begin{bmatrix} \Delta\tilde{\bf R}_{jj-1} &0 &0 \\ -\Delta \tilde {\bf R}_{ij-1}(\tilde{\bf a}_{j-1}-{\bf b}_i^a)^\wedge\Delta t & {\bf I} & 0 \\ -\frac{1}{2}\Delta \tilde {\bf R}_{ij-1}(\tilde{\bf a}_{j-1}-{\bf b}_i^a)^\wedge\Delta t^2 &\Delta t{\bf I} &{\bf I} \end{bmatrix} \tag{47}
Aj−1=⎣⎡ΔR~jj−1−ΔR~ij−1(a~j−1−bia)∧Δt−21ΔR~ij−1(a~j−1−bia)∧Δt20IΔtI00I⎦⎤(47)
B j − 1 = [ J r j − 1 Δ t 0 0 Δ R ~ i j − 1 Δ t 0 1 2 Δ R ~ i j − 1 Δ t 2 ] (48) {\bf B}_{j-1} = \begin{bmatrix} {\bf J}_r^{j-1}\Delta t & 0 \\ 0 &\Delta \tilde {\bf R}_{ij-1}\Delta t \\ 0& \frac{1}{2}\Delta \tilde {\bf R}_{ij-1}\Delta t^2 \end{bmatrix} \tag{48} Bj−1=⎣⎡Jrj−1Δt000ΔR~ij−1Δt21ΔR~ij−1Δt2⎦⎤(48)
据此可以求得预积分测量噪声的协方差矩阵:
Σ
i
j
=
A
j
−
1
Σ
i
j
−
1
A
j
−
1
T
+
B
j
−
1
Σ
η
B
j
−
1
T
(49)
\Sigma_{ij}={\bf A}_{j-1}\Sigma_{ij-1}{\bf A}_{j-1}^T+{\bf B}_{j-1}\Sigma_{\eta}{\bf B}^T_{j-1} \tag{49}
Σij=Aj−1Σij−1Aj−1T+Bj−1ΣηBj−1T(49)
C. Incorporating Bias Updates
之前的计算都是在bias恒定的假设下进行的,但这是不现实的。重新计算预积分值会耗费大量的资源,因此使用一阶展开式来计算在bias增加一个小量时的预积分值的更新。
Δ
R
~
i
j
(
b
^
i
g
)
≈
Δ
R
~
i
j
(
b
‾
i
g
)
⋅
E
x
p
(
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
b
i
g
)
Δ
v
~
i
j
(
b
^
i
g
,
b
^
i
a
)
≈
Δ
v
~
i
j
(
b
‾
i
g
,
b
‾
i
a
)
+
∂
Δ
v
‾
i
j
∂
b
‾
g
δ
b
i
g
+
∂
Δ
v
‾
i
j
∂
b
‾
a
δ
b
i
a
Δ
p
~
i
j
(
b
^
i
g
,
b
^
i
a
)
≈
Δ
v
~
i
j
(
b
‾
i
g
,
b
‾
i
a
)
+
∂
Δ
p
‾
i
j
∂
b
‾
g
δ
b
i
g
+
∂
Δ
p
‾
i
j
∂
b
‾
a
δ
b
i
a
(50)
\Delta \tilde {\bf R}_{ij} ({\bf \hat b}^g_i )\approx\Delta \tilde {\bf R}_{ij} (\overline {\bf b}^g_i) \cdot {\rm Exp} \left( \frac{\partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i \right) \\ \Delta \tilde {\bf v}_{ij} ({\bf \hat b}^g_i ,{\bf \hat b}^a_i) \approx\Delta \tilde {\bf v}_{ij} ({\bf \overline b}^g_i ,{\bf \overline b}^a_i)+\frac{\partial \Delta \overline {\bf v}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i+\frac{\partial \Delta \overline {\bf v}_{ij} }{\partial \overline {\bf b}^a} \delta{\bf b}^a_i \\ \Delta \tilde {\bf p}_{ij} ({\bf \hat b}^g_i ,{\bf \hat b}^a_i) \approx \Delta \tilde {\bf v}_{ij} ({\bf \overline b}^g_i ,{\bf \overline b}^a_i)+\frac{\partial \Delta \overline {\bf p}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i+\frac{\partial \Delta \overline {\bf p}_{ij} }{\partial \overline {\bf b}^a} \delta{\bf b}^a_i \tag{50}
ΔR~ij(b^ig)≈ΔR~ij(big)⋅Exp(∂bg∂ΔRijδbig)Δv~ij(b^ig,b^ia)≈Δv~ij(big,bia)+∂bg∂Δvijδbig+∂ba∂ΔvijδbiaΔp~ij(b^ig,b^ia)≈Δv~ij(big,bia)+∂bg∂Δpijδbig+∂ba∂Δpijδbia(50)
为了求得雅克比系数,进行如下定义:
Δ
R
^
i
j
=
Δ
R
~
i
j
(
b
^
i
g
)
,
Δ
R
‾
i
j
=
Δ
R
~
i
j
(
b
‾
i
g
)
Δ
v
^
i
j
=
Δ
v
~
i
j
(
b
^
i
g
)
,
Δ
v
‾
i
j
=
Δ
v
~
i
j
(
b
‾
i
g
)
Δ
p
^
i
j
=
Δ
p
~
i
j
(
b
^
i
g
)
,
Δ
p
‾
i
j
=
Δ
p
~
i
j
(
b
‾
i
g
)
(51)
\Delta \hat {\bf R}_{ij}=\Delta \tilde {\bf R}_{ij} ({\bf \hat b}^g_i ) , \ \Delta \overline {\bf R}_{ij}=\Delta \tilde {\bf R}_{ij} ({\bf \overline b}^g_i ) \\ \Delta \hat {\bf v}_{ij}=\Delta \tilde {\bf v}_{ij} ({\bf \hat b}^g_i ) , \ \Delta \overline {\bf v}_{ij}=\Delta \tilde {\bf v}_{ij} ({\bf \overline b}^g_i ) \\ \Delta \hat {\bf p}_{ij}=\Delta \tilde {\bf p}_{ij} ({\bf \hat b}^g_i ) , \ \Delta \overline {\bf p}_{ij}=\Delta \tilde {\bf p}_{ij} ({\bf \overline b}^g_i ) \tag{51}
ΔR^ij=ΔR~ij(b^ig), ΔRij=ΔR~ij(big)Δv^ij=Δv~ij(b^ig), Δvij=Δv~ij(big)Δp^ij=Δp~ij(b^ig), Δpij=Δp~ij(big)(51)
使用与式(32)和(43)相同的方法,推导得到
Δ
R
^
i
j
=
Δ
R
~
i
j
(
b
^
i
g
)
=
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
i
g
)
Δ
t
)
=
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
‾
i
g
−
δ
b
i
g
)
Δ
t
)
=
∏
k
=
i
j
−
1
E
x
p
(
(
ω
~
k
−
b
‾
i
g
)
Δ
t
)
E
x
p
(
−
J
r
k
δ
b
i
g
Δ
t
)
=
Δ
R
‾
i
j
∏
k
=
i
j
−
1
E
x
p
(
−
Δ
R
~
k
+
1
j
(
b
‾
i
)
T
J
r
k
δ
b
i
g
Δ
t
)
=
Δ
R
‾
i
j
E
x
p
(
∑
k
=
i
j
−
1
−
Δ
R
~
k
+
1
j
(
b
‾
i
)
T
J
r
k
δ
b
i
g
Δ
t
)
(52)
\begin{aligned} \Delta \hat {\bf R}_{ij}&=\Delta \tilde {\bf R}_{ij} ({\bf \hat b}^g_i ) \\ &=\prod_{k=i}^{j-1}\ {\rm Exp}((\tilde \omega_k-b_i^g)\Delta t) \\ &=\prod_{k=i}^{j-1}\ {\rm Exp}((\tilde \omega_k-\overline b_i^g - \delta b_i^g)\Delta t ) \\ &=\prod_{k=i}^{j-1}\ {\rm Exp}((\tilde \omega_k- \overline b_i^g)\Delta t) {\rm Exp}(-{\rm J}_r^k \delta{\bf b}^g_i\Delta t) \\ &=\Delta \overline {\bf R}_{ij} \prod_{k=i}^{j-1}{\rm Exp}(-\Delta \tilde {\bf R}_{k+1j}(\overline{\bf b}_i)^{\rm T} {\rm J}_r^k \delta{\bf b}^g_i\Delta t) \\ &=\Delta \overline {\bf R}_{ij}{\rm Exp}( \sum_{k=i}^{j-1} -\Delta \tilde {\bf R}_{k+1j}(\overline{\bf b}_i)^{\rm T} {\rm J}_r^k \delta{\bf b}^g_i\Delta t) \end{aligned} \tag{52}
ΔR^ij=ΔR~ij(b^ig)=k=i∏j−1 Exp((ω~k−big)Δt)=k=i∏j−1 Exp((ω~k−big−δbig)Δt)=k=i∏j−1 Exp((ω~k−big)Δt)Exp(−JrkδbigΔt)=ΔRijk=i∏j−1Exp(−ΔR~k+1j(bi)TJrkδbigΔt)=ΔRijExp(k=i∑j−1−ΔR~k+1j(bi)TJrkδbigΔt)(52)
将其代入
Δ
p
^
i
j
,
Δ
v
^
i
j
\Delta \hat {\bf p}_{ij}, \ \Delta \hat {\bf v}_{ij}
Δp^ij, Δv^ij之中,并对Exp进行一阶近似可推导出其它雅克比系数。
∂
Δ
R
‾
i
j
∂
b
‾
g
=
−
∑
k
=
i
j
−
1
[
Δ
R
~
k
+
1
j
(
b
‾
i
)
T
J
r
k
Δ
t
]
∂
Δ
v
‾
i
j
∂
b
‾
a
=
−
∑
k
=
i
j
−
1
Δ
R
‾
i
k
Δ
t
∂
Δ
v
‾
i
j
∂
b
‾
g
=
−
∑
k
=
i
j
−
1
Δ
R
‾
i
k
(
a
‾
k
−
b
‾
i
a
)
∧
∂
Δ
R
‾
i
k
∂
b
‾
g
Δ
t
∂
Δ
p
‾
i
j
∂
b
‾
a
=
∑
k
=
i
j
−
1
∂
Δ
v
‾
i
k
∂
b
‾
a
Δ
t
−
1
2
Δ
R
‾
i
k
Δ
t
2
∂
Δ
p
‾
i
j
∂
b
‾
g
=
∑
k
=
i
j
−
1
∂
Δ
p
‾
i
k
∂
b
‾
g
Δ
t
−
1
2
Δ
R
‾
i
k
(
a
‾
k
−
b
‾
i
a
)
∧
∂
Δ
R
‾
i
k
∂
b
‾
g
Δ
t
2
\begin{aligned} \frac{\partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} &=-\sum_{k=i}^{j-1}[ \Delta \tilde {\bf R}_{k+1j}(\overline{\bf b}_i)^{\rm T} {\rm J}_r^k \Delta t] \\ \frac{\partial \Delta \overline {\bf v}_{ij} }{\partial{\bf \overline b}^a} &= -\sum_{k=i}^{j-1}\Delta \overline {\bf R}_{ik} \Delta t \\ \frac{\partial \Delta \overline {\bf v}_{ij} }{\partial{\bf \overline b}^g} &= -\sum_{k=i}^{j-1}\Delta \overline {\bf R}_{ik}(\overline {\bf a}_k-\overline{\bf b}_i^a)^\wedge \frac{\partial \Delta \overline {\bf R}_{ik} }{\partial \overline {\bf b}^g}\Delta t \\ \frac{\partial \Delta \overline {\bf p}_{ij} }{\partial{\bf \overline b}^a} &= \sum_{k=i}^{j-1} \frac{\partial \Delta \overline {\bf v}_{ik} }{\partial \overline {\bf b}^a}\Delta t- \frac{1}{2}\Delta \overline {\bf R}_{ik} \Delta t^2 \\ \frac{\partial \Delta \overline {\bf p}_{ij} }{\partial{\bf \overline b}^g} &= \sum_{k=i}^{j-1} \frac{\partial \Delta \overline {\bf p}_{ik} }{\partial \overline {\bf b}^g}\Delta t- \frac{1}{2}\Delta \overline {\bf R}_{ik}(\overline {\bf a}_k-\overline{\bf b}_i^a)^\wedge \frac{\partial \Delta \overline {\bf R}_{ik} }{\partial \overline {\bf b}^g}\Delta t^2 \end{aligned}
∂bg∂ΔRij∂ba∂Δvij∂bg∂Δvij∂ba∂Δpij∂bg∂Δpij=−k=i∑j−1[ΔR~k+1j(bi)TJrkΔt]=−k=i∑j−1ΔRikΔt=−k=i∑j−1ΔRik(ak−bia)∧∂bg∂ΔRikΔt=k=i∑j−1∂ba∂ΔvikΔt−21ΔRikΔt2=k=i∑j−1∂bg∂ΔpikΔt−21ΔRik(ak−bia)∧∂bg∂ΔRikΔt2
D. Preintegrated IMU Factors
预积分测量模型使用了大量的一阶近似和高斯分布近似,因此存在残差。残差是预积分的计算值(使用非IMU方法,如视觉)与测量值之间的差。优化最终的目的就是通过调整(lifting)导航状态使残差(的加权范数)最小化。
根据公式(39)可以得到残差的定义:
r
I
i
j
≐
[
r
Δ
R
i
j
T
,
r
Δ
v
i
j
T
,
r
Δ
p
i
j
T
]
T
∈
R
9
(53)
{\bf r}_{\mathcal I_{ij}} \doteq [{\bf r}^{\rm T}_{\Delta {\bf R}_{ij}},\ {\bf r}^{\rm T}_{\Delta {\bf v}_{ij}}, \ {\bf r}^{\rm T}_{\Delta {\bf p}_{ij}}]^{\rm T} \in \R^9 \tag{53}
rIij≐[rΔRijT, rΔvijT, rΔpijT]T∈R9(53)
r Δ R i j ≐ L o g ( ( Δ R ~ i j ( b ‾ i g ) ⋅ E x p ( ∂ Δ R ‾ i j ∂ b ‾ g δ b i g ) ) T R i T R j ) = L o g [ ( Δ R ^ i j ) T Δ R i j ] r Δ v i j ≐ R i T ( v j − v i − g Δ t i j ) − [ Δ v ~ i j ( b ‾ i g , b ‾ i a ) + ∂ Δ v ‾ i j ∂ b ‾ g δ b i g + ∂ Δ v ‾ i j ∂ b ‾ a δ b i a ] = Δ v i j − Δ v ^ i j r Δ p i j ≐ R i T ( p j − p i − v i Δ t i j − 1 2 g Δ t i j 2 ) − [ Δ p ~ i j ( b ‾ i g , b ‾ i a ) + ∂ Δ p ‾ i j ∂ b ‾ g δ b i g + ∂ Δ p ‾ i j ∂ b ‾ a δ b i a ] = Δ p i j − Δ p ^ i j (54) \begin{aligned} {\bf r}_{\Delta {\bf R}_{ij}} &\doteq {\rm Log}\left( \left(\Delta \tilde {\bf R}_{ij} (\overline {\bf b}^g_i) \cdot {\rm Exp} \left( \frac{\partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i \right) \right)^{\rm T} {\bf R}_i^{\rm T}{\bf R}_j \right) \\ &= {\rm Log} [(\Delta \hat {\bf R}_{ij})^T\Delta {\bf R}_{ij}]\\ \\ {\bf r}_{\Delta {\bf v}_{ij}} &\doteq {\bf R}_i^T({\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij}) \\ &- \left[\Delta \tilde {\bf v}_{ij} ({\bf \overline b}^g_i ,{\bf \overline b}^a_i)+\frac{\partial \Delta \overline {\bf v}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i+\frac{\partial \Delta \overline {\bf v}_{ij} }{\partial \overline {\bf b}^a} \delta{\bf b}^a_i \right]\\ &= \Delta {\bf v}_{ij} -\Delta \hat {\bf v}_{ij}\\ \\ {\bf r}_{\Delta {\bf p}_{ij}} &\doteq {\bf R}_i^T({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2) \\ &- \left[\Delta \tilde {\bf p}_{ij} ({\bf \overline b}^g_i ,{\bf \overline b}^a_i)+\frac{\partial \Delta \overline {\bf p}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i+\frac{\partial \Delta \overline {\bf p}_{ij} }{\partial \overline {\bf b}^a} \delta{\bf b}^a_i \right] \\ &=\Delta {\bf p}_{ij} - \Delta \hat {\bf p}_{ij} \end{aligned} \tag{54} rΔRijrΔvijrΔpij≐Log((ΔR~ij(big)⋅Exp(∂bg∂ΔRijδbig))TRiTRj)=Log[(ΔR^ij)TΔRij]≐RiT(vj−vi−gΔtij)−[Δv~ij(big,bia)+∂bg∂Δvijδbig+∂ba∂Δvijδbia]=Δvij−Δv^ij≐RiT(pj−pi−viΔtij−21gΔtij2)−[Δp~ij(big,bia)+∂bg∂Δpijδbig+∂ba∂Δpijδbia]=Δpij−Δp^ij(54)
bias常常被当做状态量进行估计,使用估计其偏差的方式,因此全部的导航状态为 R i , p i , v i , R j , p j , v j , δ b i g , δ b i a {\bf R}_i,{\bf p}_i,{\bf v}_i,{\bf R}_j,{\bf p}_j,{\bf v}_j,{\bf \delta b}_i^g,{\bf \delta b}_i^a Ri,pi,vi,Rj,pj,vj,δbig,δbia。
注: lifting我理解的是,在待优化变量上增加一个微小量进行更新,通过调整这个微小量达到优化的目的。
D1. Jacobians of Residual Errors
(1). r Δ R i j {\bf r}_{\Delta {\bf R}_{ij}} rΔRij的雅克比
根据公式(54)的第一个式子可知:
p i , p j , v i , v j , δ b i a {\bf p}_i, \ {\bf p}_j, \ {\bf v}_i, \ {\bf v}_j, \ \delta{\bf b}_i^a pi, pj, vi, vj, δbia未出现在公式中,因此雅克比为0;
R
i
{\bf R}_i
Ri进行lifting计算雅克比:
r
Δ
R
i
j
(
R
i
E
x
p
(
δ
ϕ
i
)
)
=
L
o
g
[
(
Δ
R
^
i
j
)
T
(
R
i
E
x
p
(
δ
ϕ
i
)
)
T
R
j
]
(
展
开
转
置
并
将
装
置
用
负
号
代
替
)
=
L
o
g
[
(
Δ
R
^
i
j
)
T
E
x
p
(
−
δ
ϕ
i
)
R
i
T
R
j
]
(
使
用
A
d
j
o
i
n
t
性
质
)
=
L
o
g
[
(
Δ
R
^
i
j
)
T
R
i
T
R
j
E
x
p
(
−
R
j
T
R
i
δ
ϕ
i
)
]
(
前
边
是
原
残
差
公
式
)
=
L
o
g
[
E
x
p
(
r
Δ
R
i
j
(
R
i
)
)
⋅
E
x
p
(
−
R
j
T
R
i
δ
ϕ
i
)
]
(
使
用
B
C
H
近
似
公
式
)
≈
r
Δ
R
i
j
(
R
i
)
−
J
r
−
1
(
r
Δ
R
i
j
(
R
i
)
)
R
j
T
R
i
δ
ϕ
i
=
r
Δ
R
i
j
−
J
r
−
1
(
r
Δ
R
i
j
)
R
j
T
R
i
δ
ϕ
i
\begin{aligned} &{\bf r}_{\Delta{\bf R}_{ij}}({\bf R}_i{\rm Exp}(\delta\phi_i)) \\ &={\rm Log}\left[ (\Delta\hat{\bf R}_{ij})^T({\bf R}_i{\rm Exp}(\delta\phi_i))^T{\bf R}_j \right] {\text (展开转置并将装置用负号代替)} \\ &={\rm Log}\left[ (\Delta\hat{\bf R}_{ij})^T{\rm Exp}(-\delta\phi_i){\bf R}_i^T{\bf R}_j \right] {\text ( 使用Adjoint性质} )\\ &={\rm Log}\left[ (\Delta\hat{\bf R}_{ij})^T{\bf R}_i^T{\bf R}_j{\rm Exp}(-{\bf R}_j^T{\bf R}_i\delta\phi_i) \right] {\text (前边是原残差公式)} \\ &={\rm Log}\left [ {\rm Exp}\left({\bf r}_{\Delta{\bf R}_{ij}}({\bf R}_i) \right)\cdot {\rm Exp}(-{\bf R}_j^T{\bf R}_i\delta\phi_i) \right ] {\text(使用BCH近似公式)} \\ &\approx{\bf r}_{\Delta{\bf R}_{ij}}({\bf R}_i)-{\bf J}_r^{-1}\left({\bf r}_{\Delta{\bf R}_{ij}}({\bf R}_i)\right){\bf R}_j^T{\bf R}_i\delta\phi_i \\ &={\bf r}_{\Delta{\bf R}_{ij}}-{\bf J}_r^{-1}\left({\bf r}_{\Delta{\bf R}_{ij}}\right){\bf R}_j^T{\bf R}_i\delta\phi_i \end{aligned}
rΔRij(RiExp(δϕi))=Log[(ΔR^ij)T(RiExp(δϕi))TRj](展开转置并将装置用负号代替)=Log[(ΔR^ij)TExp(−δϕi)RiTRj](使用Adjoint性质)=Log[(ΔR^ij)TRiTRjExp(−RjTRiδϕi)](前边是原残差公式)=Log[Exp(rΔRij(Ri))⋅Exp(−RjTRiδϕi)](使用BCH近似公式)≈rΔRij(Ri)−Jr−1(rΔRij(Ri))RjTRiδϕi=rΔRij−Jr−1(rΔRij)RjTRiδϕi
R
j
{\bf R}_j
Rj:
r
Δ
R
i
j
(
R
i
E
x
p
(
δ
ϕ
i
)
)
=
L
o
g
[
(
Δ
R
^
i
j
)
T
R
i
T
R
j
E
x
p
(
δ
ϕ
i
)
]
(
前
边
原
公
式
,
并
利
用
B
C
H
公
式
)
=
r
Δ
R
i
j
(
R
j
)
+
J
r
−
1
(
r
Δ
R
i
j
(
R
j
)
)
δ
ϕ
j
=
r
Δ
R
i
j
+
J
r
−
1
(
r
Δ
R
i
j
)
δ
ϕ
j
\begin{aligned} &{\bf r}_{\Delta{\bf R}_{ij}}({\bf R}_i{\rm Exp}(\delta\phi_i)) \\ &={\rm Log}\left[ (\Delta\hat{\bf R}_{ij})^T{\bf R}_i^T{\bf R}_j{\rm Exp}(\delta\phi_i) \right] {\text (前边原公式,并利用BCH公式)} \\ &={\bf r}_{\Delta{\bf R}_{ij}}({\bf R}_j)+{\bf J}_r^{-1}({\bf r}_{\Delta{\bf R}_{ij}}({\bf R}_j)) \delta \phi_j \\ &={\bf r}_{\Delta{\bf R}_{ij}}+{\bf J}_r^{-1}({\bf r}_{\Delta{\bf R}_{ij}}) \delta \phi_j \end{aligned}
rΔRij(RiExp(δϕi))=Log[(ΔR^ij)TRiTRjExp(δϕi)](前边原公式,并利用BCH公式)=rΔRij(Rj)+Jr−1(rΔRij(Rj))δϕj=rΔRij+Jr−1(rΔRij)δϕj
$\delta{\bf b}_i^g $:
r
Δ
R
i
j
(
δ
b
i
g
+
δ
~
b
i
g
)
=
L
o
g
{
[
Δ
R
~
i
j
(
b
‾
i
g
)
⋅
E
x
p
(
∂
Δ
R
‾
i
j
∂
b
‾
g
(
δ
b
i
g
+
δ
~
b
i
g
)
]
T
R
i
T
R
j
}
(
将
E
x
p
使
用
B
C
H
近
似
公
式
展
开
)
≈
L
o
g
{
[
Δ
R
~
i
j
(
b
‾
i
g
)
⋅
E
x
p
(
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
b
i
g
)
E
x
p
(
J
r
(
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
b
i
g
)
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
~
b
i
g
)
]
T
Δ
R
i
j
}
(
第
二
个
E
x
p
之
前
是
预
积
分
增
量
)
=
L
o
g
{
[
Δ
R
^
i
j
⋅
E
x
p
(
J
r
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
~
b
i
g
)
]
T
Δ
R
i
j
}
(
展
开
转
置
部
分
,
合
并
第
二
部
分
)
=
L
o
g
[
E
x
p
(
−
J
r
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
~
b
i
g
)
E
x
p
(
r
Δ
R
i
j
(
δ
b
i
g
)
)
]
(
使
用
A
d
j
o
i
n
t
性
质
,
把
微
小
量
换
到
右
侧
为
了
J
r
)
=
L
o
g
[
E
x
p
(
r
Δ
R
i
j
(
δ
b
i
g
)
)
E
x
p
(
−
E
x
p
(
−
r
Δ
R
i
j
(
δ
b
i
g
)
)
J
r
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
~
b
i
g
)
]
(
使
用
B
C
H
近
似
公
式
,
这
里
有
两
个
J
r
注
意
区
分
)
=
r
Δ
R
i
j
−
J
r
−
1
(
r
Δ
R
i
j
)
⋅
E
x
p
(
−
r
Δ
R
i
j
)
⋅
J
r
(
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
b
i
g
)
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
~
b
i
g
\begin{aligned} &{\bf r}_{\Delta{\bf R}_{ij}}(\delta{\bf b}_i^g+\tilde\delta{\bf b}_i^g) \\ &={\rm Log}\left\{ \left [ \Delta \tilde {\bf R}_{ij} (\overline {\bf b}^g_i) \cdot {\rm Exp} \left( \frac{\partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} (\delta{\bf b}^g_i + \tilde \delta{\bf b}_i^g\right) \right]^T {\bf R}_i^T{\bf R}_j\right\} \\ &{\text (将Exp使用BCH近似公式展开) }\\ &\approx {\rm Log}\left\{ \left [ \Delta \tilde {\bf R}_{ij} (\overline {\bf b}^g_i) \cdot {\rm Exp} \left( \frac{\partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i \right) {\rm Exp} \left({\bf J}_r \left( \frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i\right)\frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \tilde \delta{\bf b}_i^g \right) \right]^T \Delta{\bf R}_{ij}\right\} \\ &{\text (第二个Exp之前是预积分增量) }\\ &= {\rm Log}\left\{ \left [ \Delta \hat {\bf R}_{ij} \cdot {\rm Exp} \left({\bf J}_r \frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \tilde \delta{\bf b}_i^g \right) \right]^T \Delta{\bf R}_{ij}\right\} \\ &{\text (展开转置部分,合并第二部分) }\\ &= {\rm Log}\left[ {\rm Exp} \left( -{\bf J}_r \frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \tilde \delta{\bf b}_i^g \right) {\rm Exp}({\bf r}_{\Delta{\bf R}_{ij}}(\delta{\bf b}_i^g)) \right] \\ &{\text (使用Adjoint性质,把微小量换到右侧为了J_r) } \\ &= {\rm Log}\left[ {\rm Exp}({\bf r}_{\Delta{\bf R}_{ij}}(\delta{\bf b}_i^g)) {\rm Exp} \left( - {\rm Exp}(-{\bf r}_{\Delta{\bf R}_{ij}}(\delta{\bf b}_i^g)){\bf J}_r \frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \tilde \delta{\bf b}_i^g \right) \right] \\ &{\text (使用BCH近似公式,这里有两个J_r注意区分) } \\ &={\bf r}_{\Delta{\bf R}_{ij}}-{\bf J}_r^{-1}({\bf r}_{\Delta{\bf R}_{ij}})\cdot{\rm Exp}(-{\bf r}_{\Delta{\bf R}_{ij}})\cdot {\bf J}_r \left( \frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i\right) \frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \tilde \delta{\bf b}_i^g \end{aligned}
rΔRij(δbig+δ~big)=Log{[ΔR~ij(big)⋅Exp(∂bg∂ΔRij(δbig+δ~big)]TRiTRj}(将Exp使用BCH近似公式展开)≈Log{[ΔR~ij(big)⋅Exp(∂bg∂ΔRijδbig)Exp(Jr(∂bg∂ΔRijδbig)∂bg∂ΔRijδ~big)]TΔRij}(第二个Exp之前是预积分增量)=Log{[ΔR^ij⋅Exp(Jr∂bg∂ΔRijδ~big)]TΔRij}(展开转置部分,合并第二部分)=Log[Exp(−Jr∂bg∂ΔRijδ~big)Exp(rΔRij(δbig))](使用Adjoint性质,把微小量换到右侧为了Jr)=Log[Exp(rΔRij(δbig))Exp(−Exp(−rΔRij(δbig))Jr∂bg∂ΔRijδ~big)](使用BCH近似公式,这里有两个Jr注意区分)=rΔRij−Jr−1(rΔRij)⋅Exp(−rΔRij)⋅Jr(∂bg∂ΔRijδbig)∂bg∂ΔRijδ~big
总结:
∂
r
Δ
R
i
j
∂
δ
ϕ
i
=
−
J
r
−
1
(
r
Δ
R
i
j
)
R
j
T
R
i
∂
r
Δ
R
i
j
∂
δ
ϕ
j
=
J
r
−
1
(
r
Δ
R
i
j
)
∂
r
Δ
R
i
j
∂
δ
p
i
=
0
∂
r
Δ
R
i
j
∂
δ
p
j
=
0
∂
r
Δ
R
i
j
∂
δ
v
i
=
0
∂
r
Δ
R
i
j
∂
δ
v
j
=
0
∂
r
Δ
R
i
j
∂
δ
~
b
i
a
=
0
∂
r
Δ
R
i
j
∂
δ
~
b
i
g
=
α
\begin{aligned} \frac{\partial {\bf r}_{\Delta{\bf R}_{ij}}}{\partial\delta\phi_i}&= -{\bf J}_r^{-1}\left({\bf r}_{\Delta{\bf R}_{ij}}\right){\bf R}_j^T{\bf R}_i & \frac{\partial {\bf r}_{\Delta{\bf R}_{ij}}}{\partial\delta\phi_j} &={\bf J}_r^{-1}({\bf r}_{\Delta{\bf R}_{ij}}) \\ \frac{\partial {\bf r}_{\Delta{\bf R}_{ij}}}{\partial\delta {\bf p}_i}&=0 & \frac{\partial {\bf r}_{\Delta{\bf R}_{ij}}}{\partial\delta {\bf p}_j} &=0 \\ \frac{\partial {\bf r}_{\Delta{\bf R}_{ij}}}{\partial\delta {\bf v}_i}&=0 &\frac{\partial {\bf r}_{\Delta{\bf R}_{ij}}}{\partial\delta {\bf v}_j}& = 0\\ \frac{\partial {\bf r}_{\Delta{\bf R}_{ij}}}{\partial \tilde\delta{\bf b}^a_i}&= 0 &\frac{\partial {\bf r}_{\Delta{\bf R}_{ij}}}{\partial\tilde\delta{\bf b}^g_i} &= \alpha \\ \end{aligned}
∂δϕi∂rΔRij∂δpi∂rΔRij∂δvi∂rΔRij∂δ~bia∂rΔRij=−Jr−1(rΔRij)RjTRi=0=0=0∂δϕj∂rΔRij∂δpj∂rΔRij∂δvj∂rΔRij∂δ~big∂rΔRij=Jr−1(rΔRij)=0=0=α
其中
α
=
J
r
−
1
(
r
Δ
R
i
j
)
⋅
E
x
p
(
−
r
Δ
R
i
j
)
⋅
J
r
(
∂
Δ
R
‾
i
j
∂
b
‾
g
δ
b
i
g
)
∂
Δ
R
‾
i
j
∂
b
‾
g
\alpha ={\bf J}_r^{-1}({\bf r}_{\Delta{\bf R}_{ij}})\cdot{\rm Exp}(-{\bf r}_{\Delta{\bf R}_{ij}})\cdot {\bf J}_r \left( \frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} \delta{\bf b}^g_i\right) \frac{ \partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g}
α=Jr−1(rΔRij)⋅Exp(−rΔRij)⋅Jr(∂bg∂ΔRijδbig)∂bg∂ΔRij;
(2). r Δ v i j {\bf r}_{\Delta {\bf v}_{ij}} rΔvij的雅克比
根据公式(54)的第二个式子可知:
${\bf R}_j, \ {\bf p}_i , \ {\bf p}_j $未出现在公式中雅克比为0;
δ b i a , δ b i g \delta{\bf b}_i^a, \ \delta{\bf b}_i^g δbia, δbig在公式中是线性的因此雅克比即为系数;
R
i
{\bf R}_i
Ri:
r
Δ
v
i
j
(
R
i
E
x
p
(
δ
ϕ
i
)
)
=
(
R
i
E
x
p
(
δ
ϕ
i
)
)
T
(
v
j
−
v
i
−
g
Δ
t
i
j
)
−
Δ
v
^
i
j
(
展
开
转
置
,
对
E
x
p
进
行
一
阶
近
似
)
≈
(
I
−
(
δ
ϕ
i
)
∧
)
⋅
R
i
T
⋅
(
v
j
−
v
i
−
g
Δ
t
i
j
)
−
Δ
v
^
i
j
(
展
开
第
一
个
括
号
,
凑
成
一
项
原
残
差
量
)
=
r
Δ
v
i
j
−
(
δ
ϕ
i
)
∧
⋅
R
i
T
⋅
(
v
j
−
v
i
−
g
Δ
t
i
j
)
(
利
用
a
∧
b
=
−
b
∧
a
推
导
)
=
r
Δ
v
i
j
+
[
R
i
T
⋅
(
v
j
−
v
i
−
g
Δ
t
i
j
)
]
∧
δ
ϕ
i
\begin{aligned} &{\bf r}_{\Delta{\bf v}_{ij}}({\bf R}_i{\rm Exp}(\delta\phi_i)) \\ &=({\bf R}_i{\rm Exp}(\delta\phi_i))^T({\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij})-\Delta\hat{ \bf v}_{ij} \\ &{\text (展开转置,对Exp进行一阶近似)}\\ &\approx\left({\bf I} -(\delta\phi_i)^\wedge \right)\cdot{\bf R}^T_i \cdot({\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij})-\Delta\hat{ \bf v}_{ij} \\ &{\text (展开第一个括号,凑成一项原残差量)} \\ &={\bf r}_{\Delta{\bf v}_{ij}} - (\delta\phi_i)^\wedge\cdot{\bf R}^T_i \cdot({\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij}) \\ &{\text (利用a^\wedge b=-b^\wedge a推导)} \\ &={\bf r}_{\Delta{\bf v}_{ij}} + \left[ {\bf R}^T_i \cdot({\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij}) \right]^\wedge\delta\phi_i \end{aligned}
rΔvij(RiExp(δϕi))=(RiExp(δϕi))T(vj−vi−gΔtij)−Δv^ij(展开转置,对Exp进行一阶近似)≈(I−(δϕi)∧)⋅RiT⋅(vj−vi−gΔtij)−Δv^ij(展开第一个括号,凑成一项原残差量)=rΔvij−(δϕi)∧⋅RiT⋅(vj−vi−gΔtij)(利用a∧b=−b∧a推导)=rΔvij+[RiT⋅(vj−vi−gΔtij)]∧δϕi
v
i
{\bf v}_i
vi:
r
Δ
v
i
j
(
v
i
+
δ
v
i
)
=
R
i
T
(
v
j
−
v
i
−
δ
v
i
−
g
Δ
t
i
j
)
−
Δ
v
^
i
j
=
r
Δ
v
i
j
(
v
i
)
−
R
i
T
δ
v
i
\begin{aligned} &{\bf r}_{\Delta{\bf v}_{ij}}({\bf v}_i+\delta{\bf v}_i) \\ &={\bf R}_i^T({\bf v}_j-{\bf v}_i-\delta{\bf v}_i-{\bf g}\Delta t_{ij})-\Delta\hat{ \bf v}_{ij} \\ &={\bf r}_{\Delta{\bf v}_{ij}}({\bf v}_i) -{\bf R}_i^T\delta{\bf v}_i \end{aligned}
rΔvij(vi+δvi)=RiT(vj−vi−δvi−gΔtij)−Δv^ij=rΔvij(vi)−RiTδvi
v
j
{\bf v}_j
vj:
r
Δ
v
i
j
(
v
j
+
δ
v
j
)
=
R
i
T
(
v
j
+
δ
v
j
−
v
i
−
g
Δ
t
i
j
)
−
Δ
v
^
i
j
=
r
Δ
v
i
j
(
v
j
)
+
R
i
T
δ
v
j
\begin{aligned} &{\bf r}_{\Delta{\bf v}_{ij}}({\bf v}_j+\delta{\bf v}_j) \\ &={\bf R}_i^T({\bf v}_j+\delta{\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij})-\Delta\hat{ \bf v}_{ij} \\ &={\bf r}_{\Delta{\bf v}_{ij}}({\bf v}_j) +{\bf R}_i^T\delta{\bf v}_j \end{aligned}
rΔvij(vj+δvj)=RiT(vj+δvj−vi−gΔtij)−Δv^ij=rΔvij(vj)+RiTδvj
总结:
∂
r
Δ
v
i
j
∂
δ
ϕ
i
=
[
R
i
T
⋅
(
v
j
−
v
i
−
g
Δ
t
i
j
)
]
∧
∂
r
Δ
v
i
j
∂
δ
ϕ
j
=
0
∂
r
Δ
v
i
j
∂
δ
p
i
=
0
∂
r
Δ
v
i
j
∂
δ
p
j
=
0
∂
r
Δ
v
i
j
∂
δ
v
i
=
−
R
T
∂
r
Δ
v
i
j
∂
δ
v
j
=
R
T
∂
r
Δ
v
i
j
∂
δ
~
b
i
a
=
−
∂
Δ
v
‾
i
j
∂
b
‾
i
a
∂
r
Δ
v
i
j
∂
δ
~
b
i
g
=
−
∂
Δ
v
‾
i
j
∂
b
‾
i
g
\begin{aligned} \frac{\partial {\bf r}_{\Delta{\bf v}_{ij}}}{\partial\delta\phi_i}&= \left[ {\bf R}^T_i \cdot({\bf v}_j-{\bf v}_i-{\bf g}\Delta t_{ij}) \right]^\wedge & \frac{\partial {\bf r}_{\Delta{\bf v}_{ij}}}{\partial\delta\phi_j} &= 0\\ \frac{\partial {\bf r}_{\Delta{\bf v}_{ij}}}{\partial\delta {\bf p}_i}&=0 & \frac{\partial {\bf r}_{\Delta{\bf v}_{ij}}}{\partial\delta {\bf p}_j} &=0 \\ \frac{\partial {\bf r}_{\Delta{\bf v}_{ij}}}{\partial\delta {\bf v}_i}&=-{\bf R}^T &\frac{\partial {\bf r}_{\Delta{\bf v}_{ij}}}{\partial\delta {\bf v}_j}& = {\bf R}^T \\ \frac{\partial {\bf r}_{\Delta{\bf v}_{ij}}}{\partial \tilde\delta{\bf b}^a_i}&= - \frac{\partial \Delta \overline {\bf v}_{ij} }{\partial{\bf \overline b}^a_i} &\frac{\partial {\bf r}_{\Delta{\bf v}_{ij}}}{\partial\tilde\delta{\bf b}^g_i} &= -\frac{\partial \Delta \overline {\bf v}_{ij} }{\partial{\bf \overline b}^g_i} \\ \end{aligned}
∂δϕi∂rΔvij∂δpi∂rΔvij∂δvi∂rΔvij∂δ~bia∂rΔvij=[RiT⋅(vj−vi−gΔtij)]∧=0=−RT=−∂bia∂Δvij∂δϕj∂rΔvij∂δpj∂rΔvij∂δvj∂rΔvij∂δ~big∂rΔvij=0=0=RT=−∂big∂Δvij
(3).
r
Δ
p
i
j
{\bf r}_{\Delta {\bf p}_{ij}}
rΔpij的雅克比
根据公式(54)的第三个式子可知:
${\bf R}_j , \ {\bf v}_j $未出现在公式中雅克比为0;
δ b i a , δ b i g \delta{\bf b}_i^a, \ \delta{\bf b}_i^g δbia, δbig在公式中是线性的因此雅克比即为系数;
R
i
{\bf R}_i
Ri:
r
Δ
p
i
j
(
R
i
E
x
p
(
δ
ϕ
i
)
)
=
(
R
i
E
x
p
(
δ
ϕ
i
)
)
T
(
p
j
−
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
−
Δ
p
^
i
j
(
展
开
转
置
项
,
对
E
x
p
进
行
一
阶
近
似
)
≈
(
I
−
(
δ
ϕ
i
)
∧
)
R
i
T
(
p
j
−
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
−
Δ
p
^
i
j
(
乘
开
,
并
且
跟
之
前
一
样
)
=
r
Δ
p
i
j
+
[
R
i
T
(
p
j
−
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
]
∧
⋅
δ
ϕ
i
\begin{aligned} &{\bf r}_{\Delta{\bf p}_{ij}}({\bf R}_i{\rm Exp}(\delta\phi_i)) \\ &=({\bf R}_i{\rm Exp}(\delta\phi_i)) ^T\left({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right)-\Delta\hat{\bf p}_{ij} \\ &{\text (展开转置项,对Exp进行一阶近似)} \\ &\approx \left( {\bf I}-(\delta\phi_i)^\wedge \right){\bf R}_i^T\left({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right)-\Delta\hat{\bf p}_{ij} \\ &{\text (乘开,并且跟之前一样)} \\ &={\bf r}_{\Delta{\bf p}_{ij}}+\left[ {\bf R}_i^T\left({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right) \right]^\wedge \cdot\delta\phi_i \end{aligned}
rΔpij(RiExp(δϕi))=(RiExp(δϕi))T(pj−pi−viΔtij−21gΔtij2)−Δp^ij(展开转置项,对Exp进行一阶近似)≈(I−(δϕi)∧)RiT(pj−pi−viΔtij−21gΔtij2)−Δp^ij(乘开,并且跟之前一样)=rΔpij+[RiT(pj−pi−viΔtij−21gΔtij2)]∧⋅δϕi
p
i
{\bf p}_i
pi:
r
Δ
p
i
j
(
p
i
+
R
i
⋅
δ
p
i
)
=
R
i
T
(
p
j
−
p
i
−
R
i
⋅
δ
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
−
Δ
p
^
i
j
=
R
i
T
(
p
j
−
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
−
Δ
p
^
i
j
−
I
⋅
δ
p
i
=
r
Δ
p
i
j
−
I
⋅
δ
p
i
\begin{aligned} &{\bf r}_{\Delta{\bf p}_{ij}}({\bf p}_i+{\bf R}_i \cdot\delta{\bf p}_i) \\ &={\bf R}_i ^T\left({\bf p}_j-{\bf p}_i-{\bf R}_i \cdot\delta{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right)-\Delta\hat{\bf p}_{ij} \\ &={\bf R}_i ^T\left({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right)-\Delta\hat{\bf p}_{ij} -{\bf I}\cdot\delta{\bf p}_i\\ &={\bf r}_{\Delta{\bf p}_{ij}}-{\bf I}\cdot\delta{\bf p}_i \end{aligned}
rΔpij(pi+Ri⋅δpi)=RiT(pj−pi−Ri⋅δpi−viΔtij−21gΔtij2)−Δp^ij=RiT(pj−pi−viΔtij−21gΔtij2)−Δp^ij−I⋅δpi=rΔpij−I⋅δpi
p
j
{\bf p}_j
pj:
r
Δ
p
i
j
(
p
j
+
R
j
⋅
δ
p
j
)
=
R
i
T
(
p
j
+
R
j
⋅
δ
p
j
−
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
−
Δ
p
^
i
j
=
R
i
T
(
p
j
−
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
−
Δ
p
^
i
j
+
R
i
T
R
j
⋅
δ
p
j
=
r
Δ
p
i
j
+
R
i
T
R
j
⋅
δ
p
j
\begin{aligned} &{\bf r}_{\Delta{\bf p}_{ij}}({\bf p}_j+{\bf R}_j \cdot\delta{\bf p}_j) \\ &={\bf R}_i ^T\left({\bf p}_j+{\bf R}_j \cdot\delta{\bf p}_j -{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right)-\Delta\hat{\bf p}_{ij} \\ &={\bf R}_i ^T\left({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right)-\Delta\hat{\bf p}_{ij} +{\bf R}_i ^T{\bf R}_j \cdot\delta{\bf p}_j\\ &={\bf r}_{\Delta{\bf p}_{ij}}+{\bf R}_i ^T{\bf R}_j \cdot\delta{\bf p}_j \end{aligned}
rΔpij(pj+Rj⋅δpj)=RiT(pj+Rj⋅δpj−pi−viΔtij−21gΔtij2)−Δp^ij=RiT(pj−pi−viΔtij−21gΔtij2)−Δp^ij+RiTRj⋅δpj=rΔpij+RiTRj⋅δpj
v
i
{\bf v}_i
vi:
r
Δ
p
i
j
(
v
i
+
δ
v
i
)
=
R
i
T
(
p
j
−
p
i
−
v
i
Δ
t
i
j
−
δ
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
−
Δ
p
^
i
j
=
r
Δ
p
i
j
−
R
i
T
Δ
t
i
j
⋅
δ
v
i
\begin{aligned} &{\bf r}_{\Delta{\bf p}_{ij}}({\bf v}_i+\delta{\bf v}_i) \\ &={\bf R}_i ^T\left({\bf p}_j -{\bf p}_i-{\bf v}_i \Delta t_{ij}-\delta{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right)-\Delta\hat{\bf p}_{ij} \\ &={\bf r}_{\Delta{\bf p}_{ij}}-{\bf R}_i ^T\Delta t_{ij}\cdot\delta{\bf v}_i \end{aligned}
rΔpij(vi+δvi)=RiT(pj−pi−viΔtij−δviΔtij−21gΔtij2)−Δp^ij=rΔpij−RiTΔtij⋅δvi
总结:
∂
r
Δ
p
i
j
∂
δ
ϕ
i
=
[
R
i
T
(
p
j
−
p
i
−
v
i
Δ
t
i
j
−
1
2
g
Δ
t
i
j
2
)
]
∧
∂
r
Δ
p
i
j
∂
δ
ϕ
j
=
0
∂
r
Δ
p
i
j
∂
δ
p
i
=
−
I
3
×
1
∂
r
Δ
p
i
j
∂
δ
p
j
=
R
i
T
R
j
∂
r
Δ
p
i
j
∂
δ
v
i
=
−
R
i
T
Δ
t
i
j
∂
r
Δ
p
i
j
∂
δ
v
j
=
0
∂
r
Δ
p
i
j
∂
δ
~
b
i
a
=
−
∂
Δ
p
‾
i
j
∂
b
‾
i
a
∂
r
Δ
p
i
j
∂
δ
~
b
i
g
=
−
∂
Δ
p
‾
i
j
∂
b
‾
i
g
\begin{aligned} \frac{\partial {\bf r}_{\Delta{\bf p}_{ij}}}{\partial\delta\phi_i}&= \left[ {\bf R}_i^T\left({\bf p}_j-{\bf p}_i-{\bf v}_i \Delta t_{ij}-\frac{1}{2}{\bf g}\Delta t_{ij}^2 \right) \right]^\wedge & \frac{\partial {\bf r}_{\Delta{\bf p}_{ij}}}{\partial\delta\phi_j} &= 0\\ \frac{\partial {\bf r}_{\Delta{\bf p}_{ij}}}{\partial\delta {\bf p}_i}&=-{\bf I}_{3\times1} & \frac{\partial {\bf r}_{\Delta{\bf p}_{ij}}}{\partial\delta {\bf p}_j} &={\bf R}_i ^T{\bf R}_j \\ \frac{\partial {\bf r}_{\Delta{\bf p}_{ij}}}{\partial\delta {\bf v}_i}&=-{\bf R}_i ^T\Delta t_{ij} &\frac{\partial {\bf r}_{\Delta{\bf p}_{ij}}}{\partial\delta {\bf v}_j}& =0 \\ \frac{\partial {\bf r}_{\Delta{\bf p}_{ij}}}{\partial \tilde\delta{\bf b}^a_i}&= - \frac{\partial \Delta \overline {\bf p}_{ij} }{\partial{\bf \overline b}^a_i} &\frac{\partial {\bf r}_{\Delta{\bf p}_{ij}}}{\partial\tilde\delta{\bf b}^g_i} &= -\frac{\partial \Delta \overline {\bf p}_{ij} }{\partial{\bf \overline b}^g_i} \\ \end{aligned}
∂δϕi∂rΔpij∂δpi∂rΔpij∂δvi∂rΔpij∂δ~bia∂rΔpij=[RiT(pj−pi−viΔtij−21gΔtij2)]∧=−I3×1=−RiTΔtij=−∂bia∂Δpij∂δϕj∂rΔpij∂δpj∂rΔpij∂δvj∂rΔpij∂δ~big∂rΔpij=0=RiTRj=0=−∂big∂Δpij
七、无结构化视觉因子
对于视觉测量引入无结构化模型(structureless model),关键就是线性化的消除MAP优化中的路标变量。对于视觉测量的残差有以下模型:
∑
i
∈
K
k
∑
l
∈
C
i
∥
r
C
i
l
∥
Σ
C
2
=
∑
l
=
1
L
∑
i
∈
X
(
l
)
∥
r
C
i
l
∥
Σ
C
2
r
C
i
l
=
z
i
l
−
π
(
R
i
,
p
i
,
ρ
l
)
(55)
\sum_{i \in{\mathcal K}_k} \sum_{l\in{\mathcal C_i}}\|{\bf r}_{{\mathcal C}_{il}}\|^2_{\Sigma_{\mathcal C}}=\sum_{l=1}^{L}\sum_{i\in {\mathcal X}(l)}\|{\bf r}_{{\mathcal C}_{il}}\|^2_{\Sigma_{\mathcal C}} \\ {\bf r}_{{\mathcal C}_{il}}={\bf z}_{il}-\pi({\bf R}_i,{\bf p}_i,\rho_l) \tag{55}
i∈Kk∑l∈Ci∑∥rCil∥ΣC2=l=1∑Li∈X(l)∑∥rCil∥ΣC2rCil=zil−π(Ri,pi,ρl)(55)
公式的含义是对于L个路标中的一个,在其可被观测到的关键帧上计算测量值与重投影之间的误差的二范数。
其中
ρ
l
\rho _l
ρl是路标的位置,这会对计算产生负影响,因此structureless模型就是要避免优化路标。对代价函数进行缩放,原残差公式变为:
∑
l
=
1
L
∑
i
∈
X
(
l
)
∥
z
i
l
−
π
ˇ
(
δ
ϕ
i
,
δ
p
i
,
δ
ρ
l
)
∥
Σ
C
2
(56)
\sum_{l=1}^{L}\sum_{i\in {\mathcal X}(l)}\| {\bf z}_{il}-\check \pi({\bf \delta \phi}_i,{\delta \bf p}_i,\delta\rho_l) \|^2_{\Sigma_{\mathcal C}} \tag{56}
l=1∑Li∈X(l)∑∥zil−πˇ(δϕi,δpi,δρl)∥ΣC2(56)
对其进行线性化, 并将第二个求和写入矩阵中:
∑
l
=
1
L
∥
F
l
δ
T
X
(
l
)
+
E
l
δ
ρ
l
−
b
l
∥
2
(57)
\sum_{l=1}^{L}\| {\bf F}_l\delta{\bf T}_{\mathcal X(l)}+{\bf E}_l\delta \rho_l-{\bf b}_l \|^2 \tag{57}
l=1∑L∥FlδTX(l)+Elδρl−bl∥2(57)
F
l
,
E
l
{\bf F}_l , \ {\bf E}_l
Fl, El是雅克比系数,
b
l
{\bf b}_l
bl是线性化点处线性化后的残差。
其他参数给定时,对于
δ
ρ
l
\delta \rho_l
δρl是一个二次函数求最小值的问题,最小值取在-b/2a处,因此有:
δ
ρ
l
=
−
(
E
l
T
E
l
)
−
1
E
l
T
(
F
l
δ
T
X
(
l
)
−
b
l
)
(58)
\delta \rho_l=-({\bf E}_l^{\rm T}{\bf E}_l)^{-1}{\bf E}_l^{\rm T}({\bf F}_l \delta{\bf T}_{\mathcal X(l)}-{\bf b}_l)\tag{58}
δρl=−(ElTEl)−1ElT(FlδTX(l)−bl)(58)
式(58)代入(57)消除
δ
ρ
l
\delta \rho_l
δρl:
∑
l
=
1
L
∥
(
I
−
E
l
(
E
l
T
E
l
)
−
1
E
l
T
)
(
F
l
δ
T
X
(
l
)
−
b
l
)
∥
2
(59)
\sum_{l=1}^L\|({\bf I}-{\bf E}_l({\bf E}_l^{\rm T}{\bf E}_l)^{-1}{\bf E}_l^{\rm T})({\bf F}_l \delta{\bf T}_{\mathcal X(l)}-{\bf b}_l)\|^2 \tag{59}
l=1∑L∥(I−El(ElTEl)−1ElT)(FlδTX(l)−bl)∥2(59)
I
−
E
l
(
E
l
T
E
l
)
−
1
E
l
T
{\bf I}-{\bf E}_l({\bf E}_l^{\rm T}{\bf E}_l)^{-1}{\bf E}_l^{\rm T}
I−El(ElTEl)−1ElT是与
E
l
{\bf E}_l
El垂直的。
这个方法在BA中叫做_Schur complement trick_,标准的做法是通过回代更新 ρ l \rho_l ρl的线性化点。相对的,我们使用三角化在位姿的线性化点更新路标的位置。使用这种方法可以把大量的位姿和路标变成更小的只包括位姿的L个因子。通过路标点被多个关键帧观测来保证连续性。同样的方法在MSCKF中也有使用,但是它只线性化和吸收测量值一次,直到同一个路标点的所有测量被观测到时才处理。这不适用于基于优化的本方案,它可以对新测量量进行多次线性化和增量更新。