【VIO】第3讲 基于滑动窗口算法的VIO

1.从高斯分布到信息矩阵

(1)SLAM 问题概率建模

具体参见:https://blog.csdn.net/ASUNAchan/article/details/124654207?spm=1001.2014.3001.5501

(2)高斯分布和协方差矩阵

​ 1)多元高斯分布

​ 零均值高斯分布:
p ( x ) = 1 Z e − 1 2 x T Σ − 1 x p(x) = \frac{1}{Z}e^{-\frac{1}{2}x^T\Sigma^{-1} x} p(x)=Z1e21xTΣ1x
​ 其中: Σ \Sigma Σ 为协方差矩阵,协方差矩阵的逆记作 Λ = Σ − 1 \Lambda = \Sigma^{-1} Λ=Σ1 Σ i j = E ( x i x j ) \Sigma_{ij} = E(x_ix_j) Σij=E(xixj) 为对应元素求期望。

​ 2)例1
x 2 = v 2 x 1 = w 1 x 2 + v 1 x 3 = w 3 x 2 + v 3 x_2 = v_2 \\ x_1 = w_1 x_2 + v_1 \\ x_3 = w_3 x_2 + v_3 x2=v2x1=w1x2+v1x3=w3x2+v3
​ 推导过程略

​ 结论:

​ 1 协方差逆矩阵中如果坐标为 ( i , j ) (i, j) (i,j) 的元素为 0,表示元素 i i i j j j 关于其他变量条件独立

​ 2 协方差中非对角元素 Σ i j > 0 \Sigma_{ij} > 0 Σij>0 表示两变量是正相关

​ 3 信息矩阵中非对角元素为负数,甚至为 0。 Λ 12 < 0 \Lambda_{12} < 0 Λ12<0 表示在变量 x 3 x_3 x3 发生的条件下,元素 x 1 x_1 x1 x 2 x_2 x2 正相关。

​ 3)例2
x 2 = w 1 x 1 + w 3 x 3 + v 2 x_2 = w_1 x_1 + w_3 x_3 + v_2 x2=w1x1+w3x3+v2

2.舒尔补应用:边际概率, 条件概率

(1)舒尔补的概念

​ 给定任意的矩阵块 M:
M = [ A B C D ] M = \left[ \begin{matrix} A & B \\ C & D \end{matrix} \right] M=[ACBD]
如果,矩阵块 D 是可逆的,则 A − B D − 1 C A−BD^{−1} C ABD1C 称之为 D 关于 M 的舒尔补。
[ I − B D − 1 0 I ] [ A B C D ] [ I 0 − D − 1 C I ] = [ A − B D − 1 C 0 0 D ] \left[ \begin{matrix} I & -BD^{-1} \\ 0 & I \end{matrix} \right] \left[ \begin{matrix} A & B \\ C & D \end{matrix} \right] \left[ \begin{matrix} I & 0 \\ -D^{-1}C & I \end{matrix} \right] =\left[ \begin{matrix} A-BD^{-1}C & 0 \\ 0 & D \end{matrix} \right] [I0BD1I][ACBD][ID1C0I]=[ABD1C00D]
如果,矩阵块 A 是可逆的,则 D − C A − 1 B D − CA^{−1} B DCA1B 称之为 A 关于 M 的舒尔补。
[ I 0 − C A − 1 I ] [ A B C D ] [ I − A − 1 B 0 I ] = [ A 0 0 D − C A − 1 B ] \left[ \begin{matrix} I & 0 \\ -CA^{-1} & I \end{matrix} \right] \left[ \begin{matrix} A & B \\ C & D \end{matrix} \right] \left[ \begin{matrix} I & -A^{-1}B \\ 0 & I \end{matrix} \right] =\left[ \begin{matrix} A & 0 \\ 0 & D-CA^{-1}B \end{matrix} \right] [ICA10I][ACBD][I0A1BI]=[A00DCA1B]
然后,从对角形恢复M矩阵
[ I B D − 1 0 I ] [ A − B D − 1 C 0 0 D ] [ I 0 D − 1 C I ] = [ A B C D ] \left[ \begin{matrix} I & BD^{-1} \\ 0 & I \end{matrix} \right] \left[ \begin{matrix} A-BD^{-1}C & 0 \\ 0 & D \end{matrix} \right] \left[ \begin{matrix} I & 0 \\ D^{-1}C & I \end{matrix} \right] =\left[ \begin{matrix} A & B \\ C & D \end{matrix} \right] [I0BD1I][ABD1C00D][ID1C0I]=[ACBD]

[ I 0 C A − 1 I ] [ A 0 0 D − C A − 1 B ] [ I A − 1 B 0 I ] = [ A B C D ] \left[ \begin{matrix} I & 0 \\ CA^{-1} & I \end{matrix} \right] \left[ \begin{matrix} A & 0 \\ 0 & D-CA^{-1}B \end{matrix} \right] \left[ \begin{matrix} I & A^{-1}B \\ 0 & I \end{matrix} \right] =\left[ \begin{matrix} A & B \\ C & D \end{matrix} \right] [ICA10I][A00DCA1B][I0A1BI]=[ACBD]

所以M矩阵的逆矩阵 M − 1 M^{-1} M1 为:
[ A B C D ] − 1 = [ I 0 − D − 1 C I ] [ ( A − B D − 1 C ) − 1 0 0 D − 1 ] [ I − B D − 1 0 I ] \left[ \begin{matrix} A & B \\ C & D \end{matrix} \right]^{-1} =\left[ \begin{matrix} I & 0 \\ -D^{-1}C & I \end{matrix} \right] \left[ \begin{matrix} (A-BD^{-1}C)^{-1} & 0 \\ 0 & D^{-1} \end{matrix} \right] \left[ \begin{matrix} I & -BD^{-1} \\ 0 & I \end{matrix} \right] [ACBD]1=[ID1C0I][(ABD1C)100D1][I0BD1I]

[ A B C D ] − 1 = [ I − A − 1 B 0 I ] [ A − 1 0 0 ( D − C A − 1 B ) − 1 ] [ I 0 − C A − 1 I ] \left[ \begin{matrix} A & B \\ C & D \end{matrix} \right]^{-1} = \left[ \begin{matrix} I & -A^{-1}B \\ 0 & I \end{matrix} \right] \left[ \begin{matrix} A^{-1} & 0 \\ 0 & (D-CA^{-1}B)^{-1} \end{matrix} \right] \left[ \begin{matrix} I & 0 \\ -CA^{-1} & I \end{matrix} \right] [ACBD]1=[I0A1BI][A100(DCA1B)1][ICA10I]

(2)舒尔补应用于多元高斯分布

假设多元变量x服从高斯分布,且由两部分组成:
x = [ a , b ] T x = [a, b]^T x=[a,b]T
变量之间的协方差矩阵:
K = [ A C T C D ] K = \left[ \begin{matrix} A & C^T \\ C & D \end{matrix} \right] K=[ACCTD]
其中: A = c o v ( a , a ) A=cov(a,a) A=cov(a,a) D = c o v ( b , b ) D=cov(b, b) D=cov(b,b) C = c o v ( a , b ) C=cov(a, b) C=cov(a,b)

所以:
p ( a , b ) = p ( a ) p ( b ∣ a ) ∝ e x p ( − 1 2 [ a b ] [ A C T C D ] − 1 [ a b ] ) p(a,b)=p(a)p(b|a) \propto exp( -\frac{1}{2} \left[ \begin{matrix} a & b \end{matrix} \right] \left[ \begin{matrix} A & C^T \\ C & D \end{matrix} \right]^{-1} \left[ \begin{matrix} a \\ b \end{matrix} \right] ) p(a,b)=p(a)p(ba)exp(21[ab][ACCTD]1[ab])
使用舒尔补的公式对高斯分布分解:
p ( a , b ) ∝ e x p ( 1 2 [ ( a T A − 1 a ) + ( b − C A − 1 a ) T Δ A − 1 ( b − C A − 1 a ) ] ) ∝ e x p ( 1 2 a T A − 1 a )    e x p [ 1 2 ( b − C A − 1 a ) T Δ A − 1 ( b − C A − 1 a ) ] p(a,b)\\ \propto exp( \frac{1}{2}[(a^TA^{-1}a) + (b-CA^{-1}a)^T\Delta_A^{-1}(b-CA^{-1}a)] )\\ \propto exp( \frac{1}{2}a^TA^{-1}a)\; exp [\frac{1}{2}(b-CA^{-1}a)^T\Delta_A^{-1}(b-CA^{-1}a) ] p(a,b)exp(21[(aTA1a)+(bCA1a)TΔA1(bCA1a)])exp(21aTA1a)exp[21(bCA1a)TΔA1(bCA1a)]
其中:
p ( b ∣ a ) ∼ N ( C A − 1 a , D − C A − 1 B ) p ( a ) ∼ N ( 0 , A ) p(b|a)\sim N(CA^{-1}a, D−CA^{−1} B)\\ p(a) \sim N(0,A) p(ba)N(CA1a,DCA1B)p(a)N(0,A)
p ( a ) , p ( b ∣ a ) p(a),p(b|a) p(a),p(ba) 的信息矩阵:
在这里插入图片描述

由条件概率 p ( b ∣ a ) p(b|a) p(ba) 的协方差为 Δ A \Delta_A ΔA 及上述公式,得信息矩阵 Δ A − 1 = Λ b b \Delta_A^{-1} = \Lambda_{bb} ΔA1=Λbb

由边际概率 p ( a ) p(a) p(a) 的协方差为 A A A 及上述公式,得信息矩阵 A − 1 = Λ a a − Λ a b Λ b b − 1 Λ b a A^{-1} = \Lambda_{aa} - \Lambda_{ab}\Lambda_{bb}^{-1}\Lambda_{ba} A1=ΛaaΛabΛbb1Λba

(3)总结

边际概率对于协方差矩阵的操作是很容易的,但不好操作信息矩阵。条件概率恰好相反,对于信息矩阵容易操作,不好操作协方差矩阵。
在这里插入图片描述

3.滑动窗口算法

(1)例子

​ 有如下最小二乘系统,对应的图模型如有图所示


ξ = a r g    m i n ξ    1 2 ∑ i ∣ ∣ r i ∣ ∣ Σ i 2 \xi = arg\;min_{\xi}\;\frac{1}{2}\sum_i ||r_i||^2_{\Sigma_i} ξ=argminξ21i∣∣riΣi2
其中: ξ = [ ξ 1 , ξ 2 , ⋯   , ξ 6 ] \xi = [\xi_1, \xi_2, \cdots, \xi_6] ξ=[ξ1,ξ2,,ξ6]^T, r = [ r 12 , r 13 , r 14 , r 15 , r 56 ] T r=[r_{12},r_{13},r_{14},r_{15},r_{56}]^T r=[r12,r13,r14,r15,r56]T

针对上述最小二乘问题,对应高斯牛顿求解:
J T Σ − 1 J Δ ξ = − J T Σ − 1 r J^T\Sigma^{-1}J\Delta \xi = -J^T\Sigma^{-1}r JTΣ1JΔξ=JTΣ1r

所以可以将公式写成:
∑ i = 1 5 J i T Σ i − 1 J i Δ ξ = − ∑ i = 1 5 J i T Σ i − 1 r \sum_{i=1}^{5}J_i^T\Sigma_i^{-1}J_i\Delta \xi = -\sum_{i=1}^{5}J_i^T\Sigma_i^{-1}r i=15JiTΣi1JiΔξ=i=15JiTΣi1r

(2)信息矩阵的稀疏性

​ 由于每个残差只和某几个状态量有关,因此,雅克比矩阵求导时,无关项的雅克比为 0。比如:
J 2 = ∂ r 13 ∂ ξ = [ ∂ r 13 ∂ ξ 1 , 0 , ∂ r 13 ∂ ξ 3 , 0 , 0 , 0 ] J_2 = \frac{\partial r_{13}}{\partial \xi} = [\frac{\partial r_{13}}{\partial \xi_1}, 0, \frac{\partial r_{13}}{\partial \xi_3}, 0, 0, 0 ] J2=ξr13=[ξ1r13,0,ξ3r13,0,0,0]

Λ 2 = J 2 T Σ 2 − 1 J 2 \Lambda_2 = J_2^T\Sigma_2^{-1}J_2 Λ2=J2TΣ21J2

​ 同理,可以计算 $\Lambda_1 , \Lambda_3 , \Lambda_4 ,\Lambda_5 $,并且也是稀疏的。

​ 可以结合十四讲第9章以及第8章直接法部分来更好理解。

(3)基于边际概率的滑动窗口算法

​ 1 为什么 SLAM 需要滑动窗口算法?

​ 随着 VSLAM 系统不断往新环境探索,就会有新的相机姿态以及看到新的环境特征,最小二乘残差就会越来越多,信息矩阵越来越大,计算量将不断增加。

​ 为了保持优化变量的个数在一定范围内,需要使用滑动窗口算法动态增加或移除优化变量。

​ 2 滑动窗口算法大致流程

​ 增加新的变量进入最小二乘系统优化;如果变量数目达到了一定的维度,则移除老的变量。SLAM 系统不断循环前面两步。

​ 3 利用边际概率移除老的变量

​ 直接丢弃变量和对应的测量值,会损失信息。正确的做法是使用边际概率,将丢弃变量所携带的信息传递给剩余变量。

​ 如果是直接丢弃,信息矩阵如何变化?用边际概率来操作又会带来什么问题?

(4)例子

​ 如下图优化系统中,随着 x t + 1 x_{t+1} xt+1 的进入,变量 x t x_t xt 被移除。

​ marginalization 会使得信息矩阵变稠密!原先条件独立的变量,可能变得相关。

4.滑动窗口中的 FEJ 算法

(1)回顾例子

​ 1 如图所示,在 t ∈ [ 0 , k ] s t\in[0, k]s t[0,k]s 时刻, 系统中状态 量为 ξ i , i ∈ [ 1 , 6 ] \xi_i , i \in [1, 6] ξi,i[1,6]。第 k ′ k^′ k 时刻,加入新的观 测和状态量 ξ 7 \xi_7 ξ7

​ 2 在第 k 时刻,最小二乘优化完以后,marg 掉变量 ξ 1 \xi_1 ξ1。被 marg 的状态量记为 x m x_m xm, 剩余 的变量 ξ i , i ∈ [ 2 , 5 ] \xi_i , i ∈ [2, 5] ξi,i[2,5] 记为 x r x_r xr

​ 3 marg 发生以后, x m x_m xm 所有的变量以及对应的 测量将被丢弃。同时,这部分信息通过 marg 操作传递给了保留变量 x r x_r xr ,marg 变量的信息跟 ξ 6 \xi_6 ξ6 不相关。

​ 4 第 k ′ k^′ k 时刻,加入新的状态量 ξ 7 \xi_7 ξ7(记作 x n x_n xn) 以及对应的观测,开始新一轮最小二乘优化。

(2)marg 前后

​ 已知的是:
∑ i = 1 5 J i T Σ i − 1 J i Δ ξ = − ∑ i = 1 5 J i T Σ i − 1 r \sum_{i=1}^{5}J_i^T\Sigma_i^{-1}J_i\Delta \xi = -\sum_{i=1}^{5}J_i^T\Sigma_i^{-1}r i=15JiTΣi1JiΔξ=i=15JiTΣi1r
​ marg 前,变量 x m x_m xm 以及对应测量 S m S_m Sm 构建的最小二乘信息矩阵为:
在这里插入图片描述

​ marg 后,变量 x m x_m xm 的测量信息传递给了变量 x r x_r xr
b p ( k ) = b m r ( k ) − Λ r m ( k ) Λ m m ( k ) − 1 b m m ( k ) b_p(k) = b_{mr}(k) -\Lambda_{rm}(k)\Lambda_{mm}(k)^{-1}b_{mm}(k) bp(k)=bmr(k)Λrm(k)Λmm(k)1bmm(k)

Λ p ( k ) = Λ r r ( k ) − Λ r m ( k ) Λ m m − 1 ( k ) Λ m r ( k ) \Lambda_p(k) = \Lambda_{rr}(k) -\Lambda_{rm}(k)\Lambda_{mm}^{-1}(k)\Lambda_{mr}(k) Λp(k)=Λrr(k)Λrm(k)Λmm1(k)Λmr(k)

​ 下标 p 表示 prior. 即这些信息将构建一个关于 x r x_r xr 的先验信息。

​ 我们可以从 b p ( k ) b_p(k) bp(k), Λ p ( k ) \Lambda_p(k) Λp(k) 中反解出一个残差 r p ( k ) r_p(k) rp(k) 和对应的雅克比矩 阵 J p ( k ) J_p(k) Jp(k). 需要注意的是,随着变量 x r ( k ) x_r(k) xr(k) 的后续不断优化变化,残差 r p ( k ) r_p(k) rp(k) 或者 b p ( k ) b_p(k) bp(k) 也将跟着变化,但雅克比 J p ( k ) J_p(k) Jp(k) 则固定不变了。???

(3)新测量信息和旧测量信息构建新的系统

​ 在 k ′ k^′ k​ 时刻,新残差 r 27 r_{27} r27​ 和先验信息 b p ( k ) b_p(k) bp(k)​, Λ p ( k ) \Lambda_p(k) Λp(k)​ 以及残差 r 56 r_{56} r56​ 构建新的最小二乘问题:
b ( k ′ ) = Π T b p ( k ) − ∑ ( i , j ) ∈ S a ( k ′ ) J i j ( k ′ ) T Σ i j − 1 r i j ( k ′ ) b(k^\prime) = \Pi^T b_{p}(k) -\sum_{(i,j)\in S_a(k^\prime)} J_{ij}(k^\prime)^T\Sigma_{ij}^{-1}r_{ij}(k^\prime) b(k)=ΠTbp(k)(i,j)Sa(k)Jij(k)TΣij1rij(k)

Λ ( k ′ ) = Π T Λ p ( k ) Π − ∑ ( i , j ) ∈ S a ( k ′ ) J i j ( k ′ ) T Σ i j − 1 J i j ( k ′ ) \Lambda(k^\prime) = \Pi^T\Lambda_p(k)\Pi -\sum_{(i,j)\in S_a(k^\prime)} J_{ij}(k^\prime)^T\Sigma_{ij}^{-1}J_{ij}(k^\prime) Λ(k)=ΠTΛp(k)Π(i,j)Sa(k)Jij(k)TΣij1Jij(k)

​ 其中: Π = [ I d i m    x r      0 ] \Pi=[I_{dim\;x_r}\;\;0] Π=[Idimxr0] 将矩阵维度进行扩容, S a S_a Sa 用来表示除被 marg 掉的测量以外的其他测量,如 r 56 , r 27 r_{56}, r_{27} r56,r27

​ 出现的问题:

​ 由于被 marg 的变量以及对应的测量已被丢弃,先验信息 Λ p ( k ) \Lambda_p(k) Λp(k) 中关于 x r x_r xr 的雅克比在后续求解中没法更新。

​ 但 x r x_r xr 中的部分变量还跟其他残差有联系, 如 ξ 2 , ξ 5 \xi_2, \xi_5 ξ2,ξ5。这些残差如 r 27 r_{27} r27 ξ 2 \xi_2 ξ2 的雅克比会随着 ξ 2 \xi_2 ξ2 的迭代更新而不断在最新的线性化点处计算。

(4)滑动窗口算法导致的问题

​ 滑动窗口算法优化的时候,信息矩阵如上述公式变成了两部分,且这两部分计算雅克比时的线性化点(个人理解:时间变了,变量发生变化)不同。这可能会导致信息矩阵的零空间发生变化,从而在求解时引入错误信息。

​ 滑动窗口算法中,对于同一个变量,不同残差对其计算雅克比矩阵时线性化点可能不一致,导致信息矩阵可以分成两部分,相当于在信息矩阵中多加了一些信息,使得其零空间出现了变化。

​ 解决办法:First Estimated Jacobian

​ FEJ 算法:不同残差对同一个状态求雅克比时,线性化点必须一致。 这样就能避免零空间退化而使得不可观变量变得可观。

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值