參考文獻
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