流形上的預積分


參考文獻

  • https://blog.csdn.net/qq_26682225/article/details/72649858
  • On-Manifold Preintegration for Real-TimeVisual-Inertial Odometry
  • 【泡泡機器人原創專欄】IMU預積分總結與公式推導(一)~(四)

一、預積分的由來

[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-w713fQFN-1584383006168)(../img/預積分理論/rec-18-0417-1014.png)]

來協調imu數據和圖像之間的關係。通過重新參數化[^1]來把關鍵幀之間的IMU數據積分成相對運動的約束,避免初始條件變化 [^2]造成的重複積分, 預積分的值並不是真實的測量值.


預積分與因子圖結合, 增量式因子圖

二、相關理論

  • 基於濾波的方法能夠快速的推斷當前的狀態量,但是精度會隨着線性化誤差的累積而惡化。
  • 基於優化的完全平滑方法精度高一些,但計算量很大。
  • 固定滯後的平滑方法(滑動窗)是一個折中。

採用了structless model(邊際化),這樣可以不用延遲視覺信息的處理,同時可以多次線性化視覺測量模型。具體代碼的實現整合在了gtsam4.0


  • 參加估計的pose number:

full smoothers: 所有的pose.

fixed-lag: 滑動窗口內的pose

filter: 最近的pose

  • 代表不確定性的方法:

EKF: 協方差矩陣

信息濾波[^3]和smoothers: 信息矩陣

  • 測量模型被線性化的次數:

standard EKF :1

smoothing appoach : 多次


fixed-lag和filter的方式相比對外點有更強的魯棒性,fixed-lag在邊際化的會使H矩陣變得稠密,爲了保持稀疏性,需要丟掉一些測量點。此外兩者都需要採用first estimate jacobian來線性化,來保持系統能觀性不變,否則會引入錯誤的信息。

  1. 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. }
他也代表着一個光滑的流形(manifold), 流形的切空間(tan)表示爲李代數 S O ( 3 ) \mathcal SO(3) , 它與一個3維向量的3*3的反對稱矩陣(向量加^表示)一致.
ω = [ ω 1 ω 2 ω 3 ] = [ 0 ω 3 ω 2 ω 3 0 ω 1 ω 2 ω 1 0 ] S O ( 3 ) . {\omega ^\wedge = \begin{bmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{bmatrix}^\wedge = \begin{bmatrix} 0 && -\omega_3 && \omega_2 \\ \omega_3 && 0 && -\omega_1 \\ -\omega_2 && \omega_1 && 0 \end{bmatrix} \in \mathcal{SO} (3)}.
他們之間的指數映射和對數映射如下:
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 Exp 映射:
E x p ( ϕ ) = I + s i n ( ϕ ) ϕ ϕ + 1 c o s ( ϕ ) ϕ 2 ( ϕ ) 2 Exp({\mathbf \phi})=I+\frac{sin(\|{\mathbf \phi}\|)}{\|{\mathbf \phi}\|}{\mathbf \phi}^\wedge+\frac{1-cos(\|{\mathbf \phi}\|)}{\|{\mathbf \phi}\|^2}({\bf \phi}^\wedge)^2
L o g Log 映射:
l o g ( R ) = ϕ ( R R T ) 2 s i n ( ϕ ) , ϕ = c o s 1 ( t r ( R ) 1 2 ) log(R)=\frac{\phi \cdot ({\bf R-R}^T)}{2sin(\phi)} ,\quad 其中 \quad旋轉角 \phi=cos^{-1}(\frac{tr(R)-1}{2})
李羣的乘法並不符合正常指數的運算法則,即不等於李代數的和,因此當變化量是微小量時使用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})
其中
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
同理反之,乘法變加法也成立。

[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-FJyvRPcw-1584383006172)(../img/預積分理論/rec-18-0417-2014.png)]

圖中, J r ( ϕ ) J_r(\phi) 將切線空間中的累加增量與左側的乘法增量相關聯。
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})

      E x p ( ϕ )   R   =   R E x p ( R T ϕ ) . \iff Exp({\bf\phi}) \ R \ = \ R Exp(R^T{\bf\phi}).

注意: 這裏的 E x p Exp e x p exp 是有區別的,大寫的不需要^這個。


b.SO(3)中不確定性的描述

由前面知識可知:
R ~ = R   E x p ( ϵ ) , ϵ N ( 0 ,   Σ ) \tilde{R} = R\ Exp(\epsilon), \qquad \epsilon \sim N(0,\ \Sigma)
其中, 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,
其中 α = 1 / ( 2 π ) 3 d e t ( Σ ) \alpha = 1/ \sqrt{(2\pi)^3 det(\Sigma)} , ϵ Σ 2 = ϵ T Σ 1 ϵ \parallel \epsilon\parallel^2_\Sigma = \epsilon^T\Sigma^{-1}\epsilon 是馬氏距離的平方.

ϵ < π \parallel\epsilon\parallel < \pi 時, 根據式(8)可以推出 ϵ = L o g ( R 1 R ~ ) \epsilon = Log(R^{-1}\tilde{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,
β ( R ~ ) \beta(\tilde{R}) 是一個歸一化因子, 他等於 β ( R ~ ) = α / d e t ( J ( R ~ ) ) \beta(\tilde{R})=\alpha/\mid det(J(\tilde{R}))\mid , 由於積分變量的變化 J ( R ~ ) J r ( L o g ( R 1 R ~ ) ) J(\tilde{R})\doteq J_r(Log(R^{-1}\tilde{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}
在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 }
ϵ \epsilon 是小的協方差矩陣則, β α \beta\approx\alpha , 即 J ( R ~ ) = 1 J(\tilde{R})=1 , 並且 R ~ \tilde{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
其中const是 L o g ( β ) Log(\beta) .

他的幾何意義是SO(3)上的距離的平方, 不確定性權重爲協方差的逆 Σ 1 \Sigma^{-1} .

c.流形上的高斯牛頓法

歐幾里得空間的高斯牛頓法是迭代優化目標函數的二次近似(通常非凸),求解中通常簡化爲一系列的線性方程。

當變量屬於流形時
min x M   f ( x ) \min_{x\in\mathcal{M}}\ f(x)
爲了簡化,我們只考慮一個變量時。

這裏我們不能直接近似爲 x x 的二次函數。因爲,

  • 會使參數增多(9×3個)會導致欠定。
  • 這麼做會使結果不再是流形 M \mathcal M 上的了,不滿足SO(3)定義。

通常的做法是定義一個映射 R x \mathcal{R}_x , 它是由李代數(tan 空間)上的 δ x \delta x x M x \in \mathcal{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)).
這種重新參數化叫做lifting.。這樣就可以在歐幾里得空間(也就是tan空間)進行優化,粗略地說當前估計定義在tan空間,而附近的增量是歐幾里得空間。

在高斯牛頓框架下,求解當前估計的代價函數平方,得到最小時的李代數 δ x \delta x^* (tan空間),最終更新流形上的估計值
x ^ R x ^ ( δ x ) \hat x \gets {\mathcal R}_{\hat x}(\delta 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 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


四、最大後驗視覺慣導狀態估計

IMU的座標系與機體座標系一致,相機與機體系之間的變換關係已知.

[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-0hwMfUZZ-1584383006176)(../img/預積分理論/rec-18-0704-1743.png)]

系統狀態變量爲姿態,位置,速度,偏差:
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]
K k {\mathcal K}_k 表示到時間k所有關鍵幀的序列,我們估計所有關鍵幀的狀態:
K k { X i } i K k {\mathcal K}_k\doteq \lbrace {\rm X}_i \rbrace_{i\in{\mathcal K}_k}
測量的關鍵幀圖像 C i {\mathcal C}_i 包含許多路標點 l l , 因此有觀測值 z i l {\bf z}_{il} . I i j {\mathcal I}_{ij} 表示兩關鍵幀之間的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}
X k {\mathcal 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)
使用因子圖來求解.

使用-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}


五、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 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表示機體座標系, W表示世界座標系(靜止的), WB表示B繫到W系轉換. b {\bf 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}

在時間段 [ t , t + Δ t ] [t,t+\Delta 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}
將(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 ) {\rm R}(t) 是恆定不變的, 這樣會導致產生誤差, 但當IMU頻率較高時, 這個影響會大大減小. 當使用的IMU的頻率較小時, 需要使用更高階的數值積分方法.

其中的白噪聲離散形式協方差爲原噪聲協方差的除以時間間隔 Δ t \Delta t .


六、流形上的IMU預積分

IMU預積分的初衷,是希望借鑑純視覺SLAM中圖優化的思想,將幀與幀之間IMU相對測量信息轉換爲約束節點(載體位姿)的邊參與到優化框架中。

  • IMU預積分理論最大的貢獻是對這些IMU相對測量進行處理,使得它與絕對位姿解耦(或者只需要線性運算就可以進行校正),從而大大提高優化速度。
  • 這種優化架構還使得加計測量中不受待見的重力變成一個有利條件——重力的存在將使整個系統對絕對姿態(指相對水平地理座標系的俯仰角和橫滾角,不包括真航向)可觀。要知道純視覺VO或者SLAM是完全無法得到絕對姿態的。
  • 兩種傳感器融合可以利用冗餘測量(例如兩種方式都可以求取相對位姿)來抑制累積誤差。
  • IMU和視覺這兩種不同源的測量,也使得IMU的bias可觀,從而可以在優化中被有效估計。
  • 純單目視覺缺乏絕對尺度的問題,也可以由慣性信息的引入而得以解決。

將兩關鍵幀 i i 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}
從該公式(30)可以看出當初始值改變, 如 R i {\bf 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}
最後的式子應用了等差數列求和公式推導, 證明如下:

[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-0g4CiC37-1584383006179)(../img/預積分理論/IMU預積分證明.jpeg)]

這裏的速度和位置 Δ \Delta 值並不能代表真正的物理變化, 但卻使得它與 i i 時刻狀態和重力獨立, 能夠直接根據慣性測量值計算出來. 這些值仍是Bias的函數, 假設Bias在兩關鍵幀之間保持不變.

A. Preintegrated IMU Measurements

(1) Δ R i j \Delta {\bf R}_{ij}

公式中包含的測量噪聲使得MAP(最大後驗估計)無法應用, 我們對公式進行變換, 把噪聲獨立出來, 使得log似然更容易獲得. 假設 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}
其中:
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}
δ ϕ i j \delta \phi_{ij} 定義爲噪聲.

總結:
Δ R i j Δ R ~ i j E x p ( δ ϕ i j ) Δ R ~ i j = k = i j 1   E x p ( ( ω ~ k b i g ) Δ t ) E x p ( δ ϕ i j ) = k = i j 1 E x p ( Δ R ~ k + 1 j T J r k η k g d Δ t ) (34) \Delta {\bf R}_{ij}\doteq \Delta \tilde {\bf R}_{ij}{\rm Exp}(- \delta \phi_{ij}) \\ \Delta \tilde {\bf R}_{ij}=\prod_{k=i}^{j-1}\ {\rm Exp}((\tilde \omega_k-b_i^g)\Delta t) \\ {\rm Exp}(- \delta \phi_{ij}) =\prod_{k=i}^{j-1}{\rm Exp}(-\Delta \tilde {\bf R}_{k+1j}^T {\rm J}_r^k \eta_k^{gd} \Delta t) \tag {34}

(2) Δ v i j \Delta {\bf v}_{ij}

將式(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 Δ 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}

(3) Δ p i j \Delta{\bf p}_{ij}

同上(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 Δ 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}
公式(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}

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 \delta \phi _{ij} 兩端取-Log得到:
δ ϕ i j = L o g ( Π k = i j 1 E x p ( Δ R ~ k + 1 j T J r k η k g d Δ t ) ) (41) \delta \phi_{ij}=-{\rm Log}(\Pi_{k=i}^{j-1}{\rm Exp}(-\Delta \tilde{\bf R}_{k+1j}^{\rm T}{\rm J}^k_r\eta^{gd}_k\Delta t)) \tag{41}
其中的 η k g d \eta_k^{gd} 爲微小量, 利用公式
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}
得到
δ ϕ 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}
證明:

[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-smJg9UE6-1584383006184)(../img/預積分理論/rec-18-0808-0948.png)]

由公式(43)可知, δ ϕ i j \delta \phi _{ij} η k g d \boldsymbol{\eta}^{gd}_k 線性組合, 因此也是零均值的高斯分佈.

對於 δ v i j ,   δ p i j \delta{\bf v}_{ij}, \ \delta{\bf p}_{ij} , 根據(36)(38)的增量項公式可知它們是 δ ϕ i j \delta \phi _{ij} η k a d \boldsymbol{\eta}^{ad}_k 的線性組合, 因此也符合零均值的高斯分佈.

B1. Iterative Noise Propagation

我們可以通過遞推的形式求解 η i j Δ \boldsymbol{\eta}_{ij}^{\Delta} 和其協方差 Σ i j \boldsymbol{\Sigma}_{ij} .
δ ϕ i j = Δ R ~ j j 1 δ ϕ i j 1 + J r j 1 η j 1 g d Δ t δ v i j = δ v i j 1 + Δ R ~ i j 1 η j 1 a d Δ t Δ R ~ i j 1 ( a ~ j 1 b i a ) δ ϕ i j 1 Δ t δ p i j = δ p i j 1 + δ v i j 1 Δ t + 1 2 Δ R ~ i j 1 η j 1 a d Δ t 2 1 2 Δ R ~ i j 1 ( a ~ j 1 b i a ) δ ϕ i j 1 Δ t 2 (44) \begin{aligned} \delta \phi_{ij}&=\Delta \tilde{\bf R}_{jj-1}\delta\phi_{ij-1}+{\bf J}_r^{j-1}\boldsymbol{\eta}_{j-1}^{gd}\Delta t \\ \delta{\bf v}_{ij}&= \delta {\bf v}_{ij-1}+ \Delta \tilde {\bf R}_{ij-1}\eta_{j-1}^{ad}\Delta t -\Delta \tilde {\bf R}_{ij-1}(\tilde{\bf a}_{j-1}-{\bf b}_i^a)^\wedge\delta \phi_{ij-1}\Delta t \\ \delta{\bf p}_{ij}&=\delta{\bf p}_{ij-1}+\delta {\bf v}_{ij-1}\Delta t+\frac{1}{2} \Delta \tilde {\bf R}_{ij-1}\boldsymbol{\eta}_{j-1}^{ad}\Delta t^2 -\frac{1}{2}\Delta \tilde {\bf R}_{ij-1}(\tilde{\bf a}_{j-1}-{\bf b}_i^a)^\wedge \delta\phi_{ij-1}\Delta t^2 \end{aligned} \tag{44}
其中第一個公式證明如下, 後兩個公式即拆出一項即可.

[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-uFBSa7dN-1584383006185)(../img/預積分理論/rec-18-0808-1127.png)]

定義IMU的測量噪聲向量爲
η k d = [ η k g d ,   η k a d ] (45) \boldsymbol{\eta}^d_{k}=[\boldsymbol{\eta}^{gd}_{k}, \ \boldsymbol{\eta}^{ad}_{k}] \tag{45}
將式(44)寫成矩陣形式:
η i j Δ = A j 1 η i j 1 Δ + B j 1 η j 1 d (46) \boldsymbol{\eta}^\Delta_{ij}={\bf A}_{j-1}\boldsymbol{\eta}^\Delta_{ij-1}+{\bf B}_{j-1}\boldsymbol{\eta}^d_{j-1} \tag{46}
其中
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}

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}

據此可以求得預積分測量噪聲的協方差矩陣:
Σ 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}

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 = Δ 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}
使用與式(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}
將其代入 Δ p ^ i j ,   Δ v ^ i j \Delta \hat {\bf p}_{ij}, \ \Delta \hat {\bf v}_{ij} 之中,並對Exp進行一階近似可推導出其它雅克比係數。
Δ R i j b g = k = i j 1 [ Δ R ~ k + 1 j ( b i ) T J r k Δ t ] Δ v i j b a = k = i j 1 Δ R i k Δ t Δ v i j b g = k = i j 1 Δ R i k ( a k b i a ) Δ R i k b g Δ t Δ p i j b a = k = i j 1 Δ v i k b a Δ t 1 2 Δ R i k Δ t 2 Δ p i j b g = k = i j 1 Δ p i k b g Δ t 1 2 Δ R i k ( a k b i a ) Δ R i k b g Δ t 2 \begin{aligned} \frac{\partial \Delta \overline {\bf R}_{ij} }{\partial{\bf \overline b}^g} &=-\sum_{k=i}^{j-1}[ \Delta \tilde {\bf R}_{k+1j}(\overline{\bf b}_i)^{\rm T} {\rm J}_r^k \Delta t] \\ \frac{\partial \Delta \overline {\bf v}_{ij} }{\partial{\bf \overline b}^a} &= -\sum_{k=i}^{j-1}\Delta \overline {\bf R}_{ik} \Delta t \\ \frac{\partial \Delta \overline {\bf v}_{ij} }{\partial{\bf \overline b}^g} &= -\sum_{k=i}^{j-1}\Delta \overline {\bf R}_{ik}(\overline {\bf a}_k-\overline{\bf b}_i^a)^\wedge \frac{\partial \Delta \overline {\bf R}_{ik} }{\partial \overline {\bf b}^g}\Delta t \\ \frac{\partial \Delta \overline {\bf p}_{ij} }{\partial{\bf \overline b}^a} &= \sum_{k=i}^{j-1} \frac{\partial \Delta \overline {\bf v}_{ik} }{\partial \overline {\bf b}^a}\Delta t- \frac{1}{2}\Delta \overline {\bf R}_{ik} \Delta t^2 \\ \frac{\partial \Delta \overline {\bf p}_{ij} }{\partial{\bf \overline b}^g} &= \sum_{k=i}^{j-1} \frac{\partial \Delta \overline {\bf p}_{ik} }{\partial \overline {\bf b}^g}\Delta t- \frac{1}{2}\Delta \overline {\bf R}_{ik}(\overline {\bf a}_k-\overline{\bf b}_i^a)^\wedge \frac{\partial \Delta \overline {\bf R}_{ik} }{\partial \overline {\bf b}^g}\Delta t^2 \end{aligned}

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 Δ 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}

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

[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-XeYKcDM4-1584383006186)(../img/預積分理論/rec-18-0814-2236.png)]

: lifting我理解的是,在待優化變量上增加一個微小量進行更新,通過調整這個微小量達到優化的目的。

D1. Jacobians of Residual Errors

(1). r Δ R i j {\bf r}_{\Delta {\bf R}_{ij}} 的雅克比

根據公式(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 未出現在公式中,因此雅克比爲0;

R i {\bf 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 j {\bf 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}
$\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 δ ϕ 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}
其中 α = 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} ;

(2). r Δ v i j {\bf r}_{\Delta {\bf v}_{ij}} 的雅克比

根據公式(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 在公式中是線性的因此雅克比即爲係數;

R i {\bf 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}
v i {\bf 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}
v j {\bf 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 δ ϕ 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}
(3). r Δ p i j {\bf r}_{\Delta {\bf p}_{ij}} 的雅克比

根據公式(54)的第三個式子可知:

${\bf R}_j , \ {\bf v}_j $未出現在公式中雅克比爲0;

δ b i a ,   δ b i g \delta{\bf b}_i^a, \ \delta{\bf b}_i^g 在公式中是線性的因此雅克比即爲係數;

R i {\bf 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}
p i {\bf 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}
p j {\bf 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}
v i {\bf 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 δ ϕ 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}

七、無結構化視覺因子

對於視覺測量引入無結構化模型(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}
公式的含義是對於L個路標中的一個,在其可被觀測到的關鍵幀上計算測量值與重投影之間的誤差的二範數。

其中 ρ l \rho _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 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}

相關文章
相關標籤/搜索