参考文献
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 中
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 来线性化,来保持系统能观性不变 ,否则会引入错误的信息。
SLAM中的First-Estimates 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. }
S O ( 3 ) = { R ∈ R : R T R = I , d e t ( R ) = 1 } .
他也代表着一个光滑的流形(manifold), 流形的切空间(tan)表示为李代数S O ( 3 ) \mathcal SO(3) S O ( 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 − ω 3 0 ω 1 ω 2 − ω 1 0 ⎦ ⎤ ∈ S O ( 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}
E x p L o g : R 3 : S O ( 3 ) → → S O ( 3 ) R 3 ; ; ϕ R → → e x p ( ϕ ∧ ) l o g ( R ) E x p Exp E x p 映射: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
E x p ( ϕ ) = I + ∥ ϕ ∥ s i n ( ∥ ϕ ∥ ) ϕ ∧ + ∥ ϕ ∥ 2 1 − c o s ( ∥ ϕ ∥ ) ( ϕ ∧ ) 2 L o g Log L o g 映射: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})
l o g ( R ) = 2 s i n ( ϕ ) ϕ ⋅ ( R − R T ) , 其 中 旋 转 角 ϕ = c o s − 1 ( 2 t r ( 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})
E x p ( ϕ + δ ϕ ) ≈ E x p ( ϕ ) E x p ( J r ( ϕ ) δ ϕ )
其中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
J r ( ϕ ) = I − ∥ ϕ ∥ 2 1 − cos ( ∥ ϕ ∥ ) ϕ ∧ + ∣ ∣ ϕ ∣ ∣ 3 ∣ ∣ ϕ ∣ ∣ − sin ( ∣ ∣ ϕ ∣ ∣ ) ( ϕ ∧ ) 2
同理反之,乘法变加法也成立。
图中,J r ( ϕ ) J_r(\phi) J r ( ϕ ) 将切线空间中的累加增量与左侧的乘法增量相关联。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 E x p ( ϕ ) R T = e x p ( R ϕ ∧ R T ) = E x p ( R ϕ )
⟺ E x p ( ϕ ) R = R E x p ( R T ϕ ) .
\iff Exp({\bf\phi}) \ R \ = \ R Exp(R^T{\bf\phi}).
⟺ E x p ( ϕ ) R = R E x p ( R T ϕ ) .
注意: 这里的E x p Exp E x p 和e x p exp e x p 是有区别的,大写的不需要^这个。
b.SO(3)中不确定性的描述
由前面知识可知:R ~ = R E x p ( ϵ ) , ϵ ∼ N ( 0 , Σ )
\tilde{R} = R\ Exp(\epsilon), \qquad \epsilon \sim N(0,\ \Sigma)
R ~ = R E x p ( ϵ ) , ϵ ∼ 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,
∫ R 3 p ( ϵ ) d ϵ = ∫ R 3 α e − 2 1 ∥ ϵ ∥ Σ 2 d ϵ = 1 ,
其中α = 1 / ( 2 π ) 3 d e t ( Σ ) \alpha = 1/ \sqrt{(2\pi)^3 det(\Sigma)} α = 1 / ( 2 π ) 3 d e t ( Σ ) , ∥ ϵ ∥ Σ 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}) ϵ = L o g ( R − 1 R ~ ) 代入式(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,
∫ S O ( 3 ) β ( R ~ ) e − 2 1 ∥ L o g ( R − 1 R ~ ) ∥ Σ 2 d R ~ = 1 , β ( R ~ ) \beta(\tilde{R}) β ( R ~ ) 是一个归一化因子, 他等于β ( R ~ ) = α / ∣ d e t ( J ( R ~ ) ) ∣ \beta(\tilde{R})=\alpha/\mid det(J(\tilde{R}))\mid β ( R ~ ) = α / ∣ d e t ( 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 ~ ) ≐ J r ( L o g ( R − 1 R ~ ) )
关于积分变量的推导 :δ 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 = L o g ( E x p ( ϵ ) E x p ( ϵ + δ ϵ ) ) = L o g ( E x p ( ϵ ) E x p ( ϵ ) E x p ( J r δ ϵ ) ) = J r δ ϵ
在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 − 2 1 ∥ L o g ( R − 1 R ~ ) ∥ Σ 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 ) = 2 1 ∥ L o g ( R − 1 R ~ ) ∥ Σ 2 + c o n s t = 2 1 ∥ L o g ( R ~ − 1 R ) ∥ Σ 2 + c o n s t
其中const是L o g ( β ) Log(\beta) L o g ( β ) .
他的几何意义是SO(3)上的距离的平方, 不确定性权重为协方差的逆Σ − 1 \Sigma^{-1} Σ − 1 .
c.流形上的高斯牛顿法
欧几里得空间的高斯牛顿法是迭代优化目标函数的二次近似(通常非凸),求解中通常简化为一系列的线性方程。
当变量属于流形时min x ∈ M f ( x )
\min_{x\in\mathcal{M}}\ f(x)
x ∈ M min f ( x )
为了简化,我们只考虑一个变量时。
这里我们不能直接近似为x x x 的二次函数。因为,
会使参数增多(9×3个)会导致欠定。
这么做会使结果不再是流形M \mathcal M M 上的了,不满足SO(3)定义。
通常的做法是定义一个映射R x \mathcal{R}_x R x , 它是由李代数(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 ∈ M min f ( x ) ⇒ δ x ∈ R n min f ( R ( δ x ) ) .
这种重新参数化叫做lifting .。这样就可以在欧几里得空间(也就是tan空间)进行优化,粗略地说当前估计定义在tan空间,而附近的增量是欧几里得空间。
在高斯牛顿框架下,求解当前估计的代价函数平方,得到最小时的李代数δ x ∗ \delta x^* δ x ∗ (tan空间),最终更新流形上的估计值x ^ ← R x ^ ( δ x ∗ )
\hat x \gets {\mathcal R}_{\hat x}(\delta x^*)
x ^ ← R x ^ ( δ 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
R R ( ϕ ) = R E x p ( δ ϕ ) , δ ϕ ∈ R 3
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
R T ( δ ϕ , δ p ) = ( R E x p ( δ ϕ ) , p + R δ p ) , [ δ ϕ δ p ] ∈ R 6
四、最大后验视觉惯导状态估计
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]
x i ≐ [ R i , p i , v i , b i ] K k {\mathcal K}_k K k 表示到时间k所有关键帧的序列,我们估计所有关键帧的状态:K k ≐ { X i } i ∈ K k
{\mathcal K}_k\doteq \lbrace {\rm X}_i \rbrace_{i\in{\mathcal K}_k}
K k ≐ { X i } i ∈ K k
测量的关键帧图像C i {\mathcal C}_i C i 包含许多路标点l l l , 因此有观测值z i l {\bf z}_{il} z i l . I i j {\mathcal I}_{ij} I i j 表示两关键帧之间的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}
Z k ≐ { C i , I i j } ( i , j ) ∈ K k X k {\mathcal X}_k X k 的后验概率: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 ( 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 )
使用因子图来求解.
使用-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}
X k ∗ ≐ a r g X k min − log e p ( X k ∣ Z k ) = arg X k min ∥ 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
五、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 ω ~ W B ( t ) = B ω W B ( t ) + b g ( 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)
B a ~ ( t ) = R W B T ( t ) ( W a ( t ) − W g ) + b a ( 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}
R W B ˙ = R W B B ω W B ∧ , W v ˙ = W a , W p ˙ = W v
在时间段[ 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}
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 + 2 1 W a ( t ) Δ t 2 ( 2 8 )
将(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 ) 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 + 2 1 g Δ t 2 + 2 1 R ( t ) ( a ~ ( t ) − b a ( t ) − η a d ( t ) ) Δ t 2 ( 2 9 )
这里简化了符号的表示. 公式中我们假设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}
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 + 2 1 g Δ t 2 + 2 1 R k ( a ~ k − b k a − η k a d ) Δ t 2 ] ( 3 0 )
从该公式(30)可以看出当初始值改变, 如R i {\bf R}_i R i 改变会导致每一个公式都需要重新计算. 为了解决这个缺陷, 给出不受线性化点影响的增量表达式:Δ 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}
Δ 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 − 2 1 g Δ t i j 2 ) = k = i ∑ j − 1 [ Δ v i k Δ t + 2 1 Δ R i k ( a ~ k − b k a − η k a d ) Δ t 2 ] ( 3 1 )
最后的式子应用了等差数列求和公式推导, 证明如下:
这里的速度和位置Δ \Delta Δ 值并不能代表真正的物理变化, 但却使得它与i i i 时刻状态和重力独立, 能够直接根据惯性测量值计算出来. 这些值仍是Bias的函数, 假设Bias在两关键帧之间保持不变.
A. Preintegrated IMU Measurements
(1)Δ R i j \Delta {\bf R}_{ij} Δ R i j 项
公式中包含的测量噪声使得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}
Δ 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 ) ( 3 2 )
其中: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}
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 ) ( 3 3 ) δ ϕ i j \delta \phi_{ij} δ ϕ i j 定义为噪声.
总结:Δ 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}
Δ 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 ) ( 3 4 )
(2)Δ v i j \Delta {\bf v}_{ij} Δ v i j 项
将式(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}
Δ 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 ] ( 3 5 )
总结:Δ 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}
Δ 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 ] ( 3 6 )
(3)Δ p i j \Delta{\bf p}_{ij} Δ p i j 项
同上(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}
Δ p i j = k = i ∑ j − 1 [ ( Δ v ~ i k − δ v i k ) Δ t + 2 1 Δ R ~ i k ( I − δ ϕ i k ∧ ) ( a ~ k − b k a ) Δ t 2 − 2 1 Δ R ~ i k η k a d Δ t 2 ] = k = i ∑ j − 1 [ Δ v ~ i k Δ t + 2 1 Δ R ~ i k ( a ~ k − b k a ) Δ t 2 ] + k = i ∑ j − 1 [ 2 1 Δ R ~ i k ( a ~ k − b k a ) ∧ δ ϕ i k Δ t 2 − 2 1 Δ R ~ i k η k a d Δ t 2 − δ v i k Δ t ] ( 3 7 )
总结:Δ 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}
Δ p i j ≜ Δ p ~ i j − δ p i j Δ p ~ i j = k = i ∑ j − 1 [ Δ v ~ i k Δ t + 2 1 Δ R ~ i k ( a ~ k − b k a ) Δ t 2 ] δ p i j = k = i ∑ j − 1 [ 2 1 Δ R ~ i k η k a d Δ t 2 − 2 1 Δ R ~ i k ( a ~ k − b k a ) ∧ δ ϕ i k Δ t 2 + δ v i k Δ t ] ( 3 8 )
公式(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 ~ 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 − 2 1 g Δ t i j 2 ) + δ p i j ( 3 9 )
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}
η i j Δ ≜ [ δ ϕ i j T , δ v i j T , δ p i j T ] T ∼ N ( 0 9 × 1 , Σ i j ) ( 4 0 )
对旋转量预积分噪声δ ϕ i j \delta \phi _{ij} δ ϕ i j 两端取-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}
δ ϕ i j = − L o g ( Π k = i j − 1 E x p ( − Δ R ~ k + 1 j T J r k η k g d Δ t ) ) ( 4 1 )
其中的η k g d \eta_k^{gd} η k g d 为微小量, 利用公式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}
L o g ( E x p ( ϕ ) E x p ( δ ϕ ) ) = ϕ + J r − 1 ( ϕ ) δ ϕ ( 4 2 )
得到δ ϕ 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}
δ ϕ i j ≃ Σ k = i j − 1 Δ R ~ k + 1 j T J r k η k g d Δ t ( 4 3 )
证明:
由公式(43)可知, δ ϕ i j \delta \phi _{ij} δ ϕ i j 是η k g d \boldsymbol{\eta}^{gd}_k η k g d 线性组合, 因此也是零均值的高斯分布.
对于δ v i j , δ p i j \delta{\bf v}_{ij}, \ \delta{\bf p}_{ij} δ v i j , δ p i j , 根据(36)(38)的增量项公式可知它们是δ ϕ i j \delta \phi _{ij} δ ϕ i j 和η k a d \boldsymbol{\eta}^{ad}_k η k a d 的线性组合, 因此也符合零均值的高斯分布.
B1. Iterative Noise Propagation
我们可以通过递推的形式求解η i j Δ \boldsymbol{\eta}_{ij}^{\Delta} η i j Δ 和其协方差Σ i j \boldsymbol{\Sigma}_{ij} Σ i j .δ ϕ 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}
δ ϕ i j δ v i j δ p i j = Δ R ~ j j − 1 δ ϕ i j − 1 + J r j − 1 η j − 1 g d Δ t = δ 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 − 1 + δ v i j − 1 Δ t + 2 1 Δ R ~ i j − 1 η j − 1 a d Δ t 2 − 2 1 Δ R ~ i j − 1 ( a ~ j − 1 − b i a ) ∧ δ ϕ i j − 1 Δ t 2 ( 4 4 )
其中第一个公式证明如下, 后两个公式即拆出一项即可.
定义IMU的测量噪声向量为η k d = [ η k g d , η k a d ] (45)
\boldsymbol{\eta}^d_{k}=[\boldsymbol{\eta}^{gd}_{k}, \ \boldsymbol{\eta}^{ad}_{k}] \tag{45}
η k d = [ η k g d , η k a d ] ( 4 5 )
将式(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}
η i j Δ = A j − 1 η i j − 1 Δ + B j − 1 η j − 1 d ( 4 6 )
其中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}
A j − 1 = ⎣ ⎡ Δ R ~ j j − 1 − Δ R ~ i j − 1 ( a ~ j − 1 − b i a ) ∧ Δ t − 2 1 Δ R ~ i j − 1 ( a ~ j − 1 − b i a ) ∧ Δ t 2 0 I Δ t I 0 0 I ⎦ ⎤ ( 4 7 )
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}
B j − 1 = ⎣ ⎡ J r j − 1 Δ t 0 0 0 Δ R ~ i j − 1 Δ t 2 1 Δ R ~ i j − 1 Δ t 2 ⎦ ⎤ ( 4 8 )
据此可以求得预积分测量噪声的协方差矩阵:Σ 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}
Σ i j = A j − 1 Σ i j − 1 A j − 1 T + B j − 1 Σ η B j − 1 T ( 4 9 )
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 ~ i j ( b ^ i g ) ≈ Δ R ~ i j ( b i g ) ⋅ E x p ( ∂ b g ∂ Δ R i j δ b i g ) Δ v ~ i j ( b ^ i g , b ^ i a ) ≈ Δ v ~ i j ( b i g , b i a ) + ∂ b g ∂ Δ v i j δ b i g + ∂ b a ∂ Δ v i j δ b i a Δ p ~ i j ( b ^ i g , b ^ i a ) ≈ Δ v ~ i j ( b i g , b i a ) + ∂ b g ∂ Δ p i j δ b i g + ∂ b a ∂ Δ p i j δ b i a ( 5 0 )
为了求得雅克比系数,进行如下定义:Δ 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 ^ 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 ) ( 5 1 )
使用与式(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 ^ 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 ) ( 5 2 )
将其代入Δ p ^ i j , Δ v ^ i j \Delta \hat {\bf p}_{ij}, \ \Delta \hat {\bf v}_{ij} Δ p ^ i j , Δ v ^ i j 之中,并对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}
∂ b g ∂ Δ R i j ∂ b a ∂ Δ v i j ∂ b g ∂ Δ v i j ∂ b a ∂ Δ p i j ∂ b g ∂ Δ p i j = − k = i ∑ j − 1 [ Δ R ~ k + 1 j ( b i ) T J r k Δ t ] = − k = i ∑ j − 1 Δ R i k Δ t = − k = i ∑ j − 1 Δ R i k ( a k − b i a ) ∧ ∂ b g ∂ Δ R i k Δ t = k = i ∑ j − 1 ∂ b a ∂ Δ v i k Δ t − 2 1 Δ R i k Δ t 2 = k = i ∑ j − 1 ∂ b g ∂ Δ p i k Δ t − 2 1 Δ R i k ( a k − b i a ) ∧ ∂ b g ∂ Δ R i k Δ t 2
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}
r I i j ≐ [ r Δ R i j T , r Δ v i j T , r Δ p i j T ] T ∈ R 9 ( 5 3 )
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 Δ R i j r Δ v i j r Δ p i j ≐ L o g ( ( Δ R ~ i j ( b i g ) ⋅ E x p ( ∂ b g ∂ Δ R i j δ b i g ) ) T R i T R j ) = L o g [ ( Δ R ^ i j ) T Δ R i j ] ≐ R i T ( v j − v i − g Δ t i j ) − [ Δ v ~ i j ( b i g , b i a ) + ∂ b g ∂ Δ v i j δ b i g + ∂ b a ∂ Δ v i j δ b i a ] = Δ v i j − Δ v ^ i j ≐ R i T ( p j − p i − v i Δ t i j − 2 1 g Δ t i j 2 ) − [ Δ p ~ i j ( b i g , b i a ) + ∂ b g ∂ Δ p i j δ b i g + ∂ b a ∂ Δ p i j δ b i a ] = Δ p i j − Δ p ^ i j ( 5 4 )
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 R i , p i , v i , R j , p j , v j , δ b i g , δ b i a 。
注 : lifting我理解的是,在待优化变量上增加一个微小量进行更新,通过调整这个微小量达到优化的目的。
D1. Jacobians of Residual Errors
(1). r Δ R i j {\bf r}_{\Delta {\bf R}_{ij}} r Δ R i j 的雅克比
根据公式(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 p i , p j , v i , v j , δ b i a 未出现在公式中,因此雅克比为0;
R i {\bf R}_i R i 进行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 Δ 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 R j {\bf R}_j R j :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 Δ 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
$\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 Δ R i j ( δ b i g + δ ~ b i g ) = L o g { [ Δ R ~ i j ( b i g ) ⋅ E x p ( ∂ b g ∂ Δ R i j ( δ 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 ( ∂ b g ∂ Δ R i j δ b i g ) E x p ( J r ( ∂ b g ∂ Δ R i j δ b i g ) ∂ b g ∂ Δ R i j δ ~ b i g ) ] T Δ R i j } ( 第 二 个 E x p 之 前 是 预 积 分 增 量 ) = L o g { [ Δ R ^ i j ⋅ E x p ( J r ∂ b g ∂ Δ R i j δ ~ b i g ) ] T Δ R i j } ( 展 开 转 置 部 分 , 合 并 第 二 部 分 ) = L o g [ E x p ( − J r ∂ b g ∂ Δ R i j δ ~ 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 ∂ b g ∂ Δ R i j δ ~ 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 ( ∂ b g ∂ Δ R i j δ b i g ) ∂ b g ∂ Δ R i j δ ~ b i g
总结:∂ 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 Δ R i j ∂ δ p i ∂ r Δ R i j ∂ δ v i ∂ r Δ R i j ∂ δ ~ b i a ∂ r Δ R i j = − J r − 1 ( r Δ R i j ) R j T R i = 0 = 0 = 0 ∂ δ ϕ j ∂ r Δ R i j ∂ δ p j ∂ r Δ R i j ∂ δ v j ∂ r Δ R i j ∂ δ ~ b i g ∂ r Δ R i j = J r − 1 ( r Δ R i j ) = 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} α = J r − 1 ( r Δ R i j ) ⋅ E x p ( − r Δ R i j ) ⋅ J r ( ∂ b g ∂ Δ R i j δ b i g ) ∂ b g ∂ Δ R i j ;
(2). r Δ v i j {\bf r}_{\Delta {\bf v}_{ij}} r Δ v i j 的雅克比
根据公式(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 δ b i a , δ b i g 在公式中是线性的因此雅克比即为系数;
R i {\bf R}_i R i :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 Δ 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 v i {\bf v}_i v i :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 Δ 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 v j {\bf v}_j v j :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 Δ 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
总结:∂ 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 Δ v i j ∂ δ p i ∂ r Δ v i j ∂ δ v i ∂ r Δ v i j ∂ δ ~ b i a ∂ r Δ v i j = [ R i T ⋅ ( v j − v i − g Δ t i j ) ] ∧ = 0 = − R T = − ∂ b i a ∂ Δ v i j ∂ δ ϕ j ∂ r Δ v i j ∂ δ p j ∂ r Δ v i j ∂ δ v j ∂ r Δ v i j ∂ δ ~ b i g ∂ r Δ v i j = 0 = 0 = R T = − ∂ b i g ∂ Δ v i j (3). r Δ p i j {\bf r}_{\Delta {\bf p}_{ij}} r Δ p i j 的雅克比
根据公式(54)的第三个式子可知:
${\bf R}_j , \ {\bf v}_j $未出现在公式中雅克比为0;
δ b i a , δ b i g \delta{\bf b}_i^a, \ \delta{\bf b}_i^g δ b i a , δ b i g 在公式中是线性的因此雅克比即为系数;
R i {\bf R}_i R i :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 Δ 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 − 2 1 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 − 2 1 g Δ t i j 2 ) − Δ p ^ i j ( 乘 开 , 并 且 跟 之 前 一 样 ) = r Δ p i j + [ R i T ( p j − p i − v i Δ t i j − 2 1 g Δ t i j 2 ) ] ∧ ⋅ δ ϕ i p i {\bf p}_i p i :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 Δ 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 − 2 1 g Δ t i j 2 ) − Δ p ^ i j = R i T ( p j − p i − v i Δ t i j − 2 1 g Δ t i j 2 ) − Δ p ^ i j − I ⋅ δ p i = r Δ p i j − I ⋅ δ p i p j {\bf p}_j p j :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 Δ 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 − 2 1 g Δ t i j 2 ) − Δ p ^ i j = R i T ( p j − p i − v i Δ t i j − 2 1 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 v i {\bf v}_i v i :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 Δ p i j ( v i + δ v i ) = R i T ( p j − p i − v i Δ t i j − δ v i Δ t i j − 2 1 g Δ t i j 2 ) − Δ p ^ i j = r Δ p i j − R i T Δ t i j ⋅ δ v i
总结:∂ 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 Δ p i j ∂ δ p i ∂ r Δ p i j ∂ δ v i ∂ r Δ p i j ∂ δ ~ b i a ∂ r Δ p i j = [ R i T ( p j − p i − v i Δ t i j − 2 1 g Δ t i j 2 ) ] ∧ = − I 3 × 1 = − R i T Δ t i j = − ∂ b i a ∂ Δ p i j ∂ δ ϕ j ∂ r Δ p i j ∂ δ p j ∂ r Δ p i j ∂ δ v j ∂ r Δ p i j ∂ δ ~ b i g ∂ r Δ p i j = 0 = R i T R j = 0 = − ∂ b i g ∂ Δ p i j
七、无结构化视觉因子
对于视觉测量引入无结构化模型(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 ∈ 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 ) ( 5 5 )
公式的含义是对于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 ∑ L i ∈ X ( l ) ∑ ∥ z i l − π ˇ ( δ ϕ i , δ p i , δ ρ l ) ∥ Σ C 2 ( 5 6 )
对其进行线性化, 并将第二个求和写入矩阵中:∑ 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 ∥ F l δ T X ( l ) + E l δ ρ l − b l ∥ 2 ( 5 7 ) F l , E l {\bf F}_l , \ {\bf E}_l F l , E l 是雅克比系数,b l {\bf b}_l b l 是线性化点处线性化后的残差。
其他参数给定时,对于δ ρ 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 = − ( E l T E l ) − 1 E l T ( F l δ T X ( l ) − b l ) ( 5 8 )
式(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 − E l ( E l T E l ) − 1 E l T ) ( F l δ T X ( l ) − b l ) ∥ 2 ( 5 9 ) 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 − E l ( E l T E l ) − 1 E l T 是与E l {\bf E}_l E l 垂直的。
这个方法在BA中叫做_Schur complement trick_,标准的做法是通过回代更新ρ l \rho_l ρ l 的线性化点。相对的,我们使用三角化在位姿的线性化点更新路标的位置。使用这种方法可以把大量的位姿和路标变成更小的只包括位姿的L个因子。通过路标点被多个关键帧观测来保证连续性。同样的方法在MSCKF中也有使用,但是它只线性化和吸收测量值一次,直到同一个路标点的所有测量被观测到时才处理。这不适用于基于优化的本方案,它可以对新测量量进行多次线性化和增量更新。