文章目录
- 问题例子example
- 求解思路
- 边缘化
- 边缘化原理
- 边缘化的实际步骤
- marg先验约束公式
- 先验约束公式1
- 先验约束公式2
- marg的问题及FEJ
- 实例分析:VINS-Mono中的滑动窗口策略
- 边缘化的代码实现(伪代码)
- 参考
本文简要介绍VIO常用的滑动窗口及边缘化的原理,以及存在问题。
边缘化其实思想并不复杂,能够引入一定先验信息的同时,也会带来一些弊端,本文以简单示例,附加简单清晰的原理讲解滑动窗口及边缘化的原理。
滑动窗口是一种动态管理优化问题规模的技术,旨在维持实时性。其核心思想是仅保留最近一段时间内的状态变量(如相机位姿、路标点、IMU参数等),移除旧状态以避免计算量无限增长。
-
窗口更新策略:
- 新增状态:当新数据(如图像、IMU测量)到达时,将新状态(如位姿 T k + 1 T_{k+1} Tk+1、路标点 p n e w p_{new} pnew)加入窗口。
- 移除旧状态:当窗口容量达到上限时,按特定规则(如时间顺序或信息量评估)选择旧状态进行边缘化。
-
窗口大小选择:
- 静态窗口:固定大小(如10帧),适合计算资源受限场景。
- 动态窗口:根据场景复杂度或运动速度调整大小(如运动剧烈时增大窗口)。
问题例子example
设一个观测系统,对应的图模型如下:
图中:
圆圈:表示顶点,需要估计的变量。即该图模型中,需要求解的
连线:表示关系边,表示各个顶点(待求解变量)之间的观测约束(残差)。
根据该图模型建立的最小二乘问题如下:
ξ = a r g m i n ξ 1 2 ∑ i ∥ r i ∥ Σ i 2 \boldsymbol{\xi}=\underset{\boldsymbol{\xi}}{\operatorname*{\operatorname*{\mathrm{argmin}}}}\frac{1}{2}\sum_i\|\mathbf{r}_i\|_{\boldsymbol{\Sigma}_i}^2 ξ=ξargmin21i∑∥ri∥Σi2
其中:
ξ = [ ξ 1 ξ 2 . . . ξ 6 ] , r = [ r 12 r 13 r 14 r 15 r 56 ] \boldsymbol{\xi}= \begin{bmatrix} \boldsymbol{\xi}_{1} \\ \boldsymbol{\xi}_{2} \\ ... \\ \boldsymbol{\xi}_{6} \end{bmatrix},\mathbf{r}= \begin{bmatrix} \mathbf{r}_{12} \\ \mathbf{r}_{13} \\ \mathbf{r}_{14} \\ \mathbf{r}_{15} \\ \mathbf{r}_{56} \end{bmatrix} ξ= ξ1ξ2...ξ6 ,r= r12r13r14r15r56
ξ \boldsymbol{\xi} ξ表示待求解的变量,对应图中的圆圈节点,图中共6个节点。
r \mathbf{r} r表示观测约束,下标表示产生联系的两个节点。图中共有5个观测。
求解思路
对于上述最小二乘问题,使用高斯牛顿法进行迭代求解时,其变量更新求解公式如下:
J ⊤ Σ − 1 J ⏟ H o r Λ δ ξ = − J ⊤ Σ − 1 r ⏟ b (1) \underbrace{\mathbf{J}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{J}}_{\mathrm{H~}or\boldsymbol{\Lambda}}\delta\boldsymbol{\xi}=\underbrace{-\mathbf{J}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{r}}_{\mathbf{b}}\tag{1} H orΛ J⊤Σ−1Jδξ=b −J⊤Σ−1r(1)
J \mathbf{J} J魏雅可比矩阵。雅可比矩阵行与观测数(边线约束)相关,列与求解变量数相关。
J = ∂ r ∂ ξ = [ ∂ r 12 ∂ ξ ∂ r 13 ∂ ξ ∂ r 14 ∂ ξ ∂ r 15 ∂ ξ ∂ r 56 ∂ ξ ] = [ J 1 J 2 J 3 J 4 J 5 ] , J ⊤ = [ J 1 ⊤ J 2 ⊤ J 3 ⊤ J 4 ⊤ J 5 ⊤ ] \mathbf{J}=\frac{\partial\mathbf{r}}{\partial\boldsymbol{\xi}}= \begin{bmatrix} \frac{\partial\mathbf{r}_{12}}{\partial\boldsymbol{\xi}} \\ \frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}} \\ \frac{\partial\mathbf{r}_{14}}{\partial\boldsymbol{\xi}} \\ \frac{\partial\mathbf{r}_{15}}{\partial\boldsymbol{\xi}} \\ \frac{\partial\mathbf{r}_{56}}{\partial\boldsymbol{\xi}} \end{bmatrix}= \begin{bmatrix} \mathbf{J}_{1} \\ \mathbf{J}_{2} \\ \mathbf{J}_{3} \\ \mathbf{J}_{4} \\ \mathbf{J}_{5} \end{bmatrix},\mathbf{J}^{\top}= \begin{bmatrix} \mathbf{J}_{1}^{\top} & \mathbf{J}_{2}^{\top} & \mathbf{J}_{3}^{\top} & \mathbf{J}_{4}^{\top} & \mathbf{J}_{5}^{\top} \end{bmatrix} J=∂ξ∂r= ∂ξ∂r12∂ξ∂r13∂ξ∂r14∂ξ∂r15∂ξ∂r56 = J1J2J3J4J5 ,J⊤=[J1⊤J2⊤J3⊤J4⊤J5⊤]
将上述雅可比矩阵带入公式 ( 1 ) (1) (1)可得参数更新的残差累加形式:
∑ i = 1 5 J i ⊤ Σ i − 1 J i δ ξ = − ∑ i = 1 5 J i ⊤ Σ i − 1 r i \sum_{i=1}^5\mathbf{J}_i^\top\boldsymbol{\Sigma}_i^{-1}\mathbf{J}_i\delta\boldsymbol{\xi}=-\sum_{i=1}^5\mathbf{J}_i^\top\boldsymbol{\Sigma}_i^{-1}\mathbf{r}_i i=1∑5Ji⊤Σi−1Jiδξ=−i=1∑5Ji⊤Σi−1ri
由于一个观测约束通常几个(通常为2个)状态量,如在某帧相机下对一个特征点的观测就是一个约束,这个约束涉及到的变量为一个特征点和一个相机位姿。从而雅可比矩阵通常是稀疏矩阵,无关的节点对应项为0。
如:
J 2 = ∂ r 13 ∂ ξ = [ ∂ r 13 ∂ ξ 1 0 ∂ r 13 ∂ ξ 3 0 0 0 ] H 2 = J 2 ⊤ Σ 2 − 1 J 2 = [ ( ∂ r 13 ∂ ξ 1 ) ⊤ Σ 2 − 1 ∂ r 13 ∂ ξ 1 0 ( ∂ r 13 ∂ ξ 1 ) ⊤ Σ 2 − 1 ∂ r 13 ∂ ξ 3 0 0 0 0 0 0 0 0 0 ( ∂ r 13 ∂ ξ 3 ) ⊤ Σ 2 − 1 ∂ r 13 ∂ ξ 1 0 ( ∂ r 13 ∂ ξ 3 ) ⊤ Σ 2 − 1 ∂ r 13 ∂ ξ 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ] \begin{aligned} & \mathbf{J}_2=\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}}= \begin{bmatrix} \frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_1} & \boldsymbol{0} & \frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_3} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \end{bmatrix} \\ & \boldsymbol{H}_2=\mathbf{J}_2^\top\boldsymbol{\Sigma}_2^{-1}\mathbf{J}_2= \begin{bmatrix} (\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_1})^\top\boldsymbol{\Sigma}_2^{-1}\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_1} & \boldsymbol{0} & (\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_1})^\top\boldsymbol{\Sigma}_2^{-1}\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_3} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\ (\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_3})^\top\boldsymbol{\Sigma}_2^{-1}\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_1} & \boldsymbol{0} & (\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_3})^\top\boldsymbol{\Sigma}_2^{-1}\frac{\partial\mathbf{r}_{13}}{\partial\boldsymbol{\xi}_3} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \end{bmatrix} \end{aligned} J2=∂ξ∂r13=[∂ξ1∂r130∂ξ3∂r13000]H2=J2⊤Σ2−1J2= (∂ξ1∂r13)⊤Σ2−1∂ξ1∂r130(∂ξ3∂r13)⊤Σ2−1∂ξ1∂r13000000000(∂ξ1∂r13)⊤Σ2−1∂ξ3∂r130(∂ξ3∂r13)⊤Σ2−1∂ξ3∂r13000000000000000000000
同理,可以得知 H 1 , H 3 , H 4 , H 5 \boldsymbol{H}_1,\boldsymbol{H}_3,\boldsymbol{H}_4,\boldsymbol{H}_5 H1,H3,H4,H5也是稀疏的。
进一步,将五个残差的信息矩阵加起来,得到 H = ∑ i = 1 5 J i ⊤ Σ i − 1 J i δ ξ \boldsymbol{H}=\sum_{i=1}^5\mathbf{J}_i^\top\boldsymbol{\Sigma}_i^{-1}\mathbf{J}_i\delta\boldsymbol{\xi} H=∑i=15Ji⊤Σi−1Jiδξ,对应的 b = ∑ i = 1 5 J i ⊤ Σ i − 1 r i b=\sum_{i=1}^5\mathbf{J}_i^\top\boldsymbol{\Sigma}_i^{-1}\mathbf{r}_i b=∑i=15Ji⊤Σi−1ri该矩阵的稀疏性如下图所示。
边缘化
滑动窗口方法会仅保留一定数量的最近帧(通常是当前帧和若干个前一帧),而将较旧的帧从优化中移除。这样做的目的是减少计算负担,保持实时性能,同时避免过度积累历史信息。
边缘化(Marginalization)是与滑动窗口密切相关的技术,主要用于处理窗口内移除的旧帧或者不再重要的路标点。在SLAM中,当一个变量被从滑动窗口中移除时,其状态信息不能完全丢弃,而是通过边缘化技术将其状态融合到其他变量中,以保留系统的长期一致性。边缘化通过数学手段对被移除的帧进行降维处理,减少对计算资源的需求,同时确保优化问题的解仍然尽可能准确。
边缘化原理
边缘化即将旧的变量移出窗口状态有两种理解方法,一种可以从概率的角度去分析,参考链接。另一种是基于Schur补的思想。参考链接。
以下从schur补来说明边缘化的原理:即通过消元法将需要marg的变量抹除。
设滑窗中的优化变量为 ξ = [ ξ m , ξ r ] T \xi = [\xi_m, \xi_r]^{\rm T} ξ=[ξm,ξr]T,对应的Hessian矩阵 H H H和雅可比矩阵 J J J如下
H = [ H m m H m r H r m H r r ] , b = [ J m J r ] H=\begin{bmatrix}H_{mm}&H_{mr}\\H_{rm}&H_{rr}\end{bmatrix},b=\begin{bmatrix}J_m\\J_r\end{bmatrix} H=[HmmHrmHmrHrr],b=[JmJr]
使用迭代法计算参数更新值 ξ = [ ξ m , ξ r ] T \xi = [\xi_m, \xi_r]^{\rm T} ξ=[ξm,ξr]T时,(牛顿法为例),计算公式如下:
Δ ξ = − H − 1 b \Delta \xi =-H^{-1}b Δξ=−H−1b
设 ξ m \xi_m ξm为待marg的变量,利用消元法消去上式中的 ξ m \xi_m ξm,则可以得到:
H p = H r r − H r m H m m − 1 H m r H_p = H_{rr}-H_{rm}H_{mm}^{-1}H_{mr} Hp=Hrr−HrmHmm−1Hmr b p = b r − H r m H m m − 1 b m b_p = b_r -H_{rm}H_{mm}^{-1}b_m bp=br−HrmHmm−1bm
在上述的问题例子中marg掉变量 ξ 1 \xi_1 ξ1,使用上述的marg公式的示意图如下:
SLAM中的Hessian矩阵通常比较稀疏,这是由于通常通常一帧的观测对应很多的观测点。Hessian矩阵中元素的0或非0表征者各个变量(元素)之间的联系与否。
看到经过marg边缘化之后,原先稀疏的Hessian矩阵变得稠密,从而表明。marg边缘化会使得原先没有联系的变量发生联系。
当有新的变量 ξ 7 \xi_7 ξ7进入系统,与变量 ξ 2 \xi_2 ξ2发生观测联系,如下示意图:
将新的观测与变量加入Hessian矩阵。首先将marg后得到的Hessian矩阵 H p H_p Hp进行扩容,而后加入新的变量及约束构成的Hessian矩阵 H 7 H_7 H7,如下示意图。
边缘化的实际步骤
在非线性优化中,边缘化需处理线性化点固定问题,具体步骤如下:
-
构建增量方程:
在当前线性化点(如 ξ 0 \xi_0 ξ0)处计算残差的雅可比矩阵 J J J 和 Hessian 矩阵 H = J T J H = J^T J H=JTJ。 -
选择边缘化变量:
确定要移除的变量 ξ m \xi_m ξm(如窗口中最旧的位姿 T 1 T_1 T1 及其关联的路标点)。 -
舒尔补操作:
对 H H H 和 b b b 进行分块,计算先验信息 H prior H_{\text{prior}} Hprior 和 b prior b_{\text{prior}} bprior。 -
构造先验残差:
将先验信息表示为虚拟残差项 r p r_{\text{p}} rp,其对应的优化问题为:
min Δ ξ r ∥ r p + J p ( ξ − ξ p ) ∥ 2 \min_{\Delta \xi_r} \left\| r_{\text{p}} + J_{\text{p}}(\xi-\xi_p) \right\|^2 Δξrmin∥rp+Jp(ξ−ξp)∥2 -
信息矩阵叠加:
边缘化后的先验信息作为新的约束加入后续优化中:
-
后续优化问题的 Hessian 矩阵为:
H total = H new + H p H_{\text{total}} = H_{\text{new}} + H_{\text{p}} Htotal=Hnew+Hp -
误差向量叠加:
b total = b new + b p b_{\text{total}} = b_{\text{new}} + b_{\text{p}} btotal=bnew+bp
marg先验约束公式
先验约束公式1
marg后得到的先验信息约束方程可以等价为如下公式:
e p = ∣ ∣ r p + J p ( ξ − ξ p ) ∣ ∣ (2) e_p = ||r_p+J_p(\xi-\xi_p)|| \tag2 ep=∣∣rp+Jp(ξ−ξp)∣∣(2)
注意:这里的 J p J_p Jp及 r p r_p rp等均是由 H p H_p Hp及 b p b_p bp中反构造出的常参数,通常在边缘化迭代中不需要计算,这里的公式是用于说明边缘化的数学约束意义。
H p = J p T J p H_p = J_p^TJ_p Hp=JpTJp, ξ p \xi_p ξp为在新的数据进入之前,优化计算得到的变量值。 r p r_p rp代指加入新数据之前,迭代优化最后剩余的残差。继而可得其对应的迭代公式: b p = J p T ( r p + J p ( ξ − ξ p ) ) , H p = J p T J p , Δ ξ p = − H p − 1 b p b_p=J_p^{\rm T}(r_p+J_p(\xi-\xi_p)), H_p = J_p^TJ_p,\ \ \Delta \xi_p = -H_p^{-1}b_p bp=JpT(rp+Jp(ξ−ξp)),Hp=JpTJp, Δξp=−Hp−1bp
marg后的第一次迭代中, ξ = ξ p , b p = J p T r p \xi =\xi_p,\ \ b_p=J_p^{\rm T}r_p ξ=ξp, bp=JpTrp。
通常在marg后会加入其他信息到hessian矩阵及 b b b中,综合求解新的变量更新值 Δ ξ \Delta \xi Δξ。注意这里 Δ ξ \Delta\xi Δξ表示综合新的信息算出的更新值,而 Δ ξ p \Delta\xi_p Δξp为使用先验约束算出的变量更新量。
在后续迭代时,变量更新公式: ξ ← ξ + Δ ξ \xi\leftarrow \xi+\Delta\xi ξ←ξ+Δξ。
在变量marg之后,marg得到的稠密先验矩阵 H p H_p Hp随着变量的变化不再能更新,而其偏置量 b p b_p bp可以使用如下公式进行更新:
将 b p k + 1 b_p^{k+1} bpk+1在 ξ p \xi_p ξp处求偏导:
b p k + 1 = b p k + ∂ b p k ∂ ξ ∣ ξ p δ ξ = b p + ∂ J p T ( r p + J p ( ξ − ξ p ) ) ∂ ξ ∣ ξ p δ ξ = b p + ∂ J p T r p + ∂ H p ( ξ − ξ p ) ∂ ξ δ ξ = b p + H p δ ξ = b p + H p ( ξ k + 1 − ξ p ) b_p^{k+1}=\left.b_p^k+\frac{\partial b_p^k}{\partial\xi}\right|_{\xi_p}\delta\xi\\ =b_p+\left.\frac{\partial J_p^{\rm T}(r_p+J_p(\xi-\xi_p))}{\partial\xi}\right|_{\xi_p}\delta\xi\\ =b_p+\frac{\partial J_p^{\rm T}r_p+\partial H_p(\xi-\xi_p)}{\partial\xi}\delta\xi\\ =b_p+H_p\delta\xi\\ =b_p+H_p(\xi^{k+1}-\xi_p) bpk+1=bpk+∂ξ∂bpk ξpδξ=bp+∂ξ∂JpT(rp+Jp(ξ−ξp)) ξpδξ=bp+∂ξ∂JpTrp+∂Hp(ξ−ξp)δξ=bp+Hpδξ=bp+Hp(ξk+1−ξp)
也可以使用推导得到偏置的更新公式。
b p k + 1 = J p T ( r p + J p ( ξ k + Δ ξ − ξ p ) ) = b p + H p ( ξ k + 1 − ξ p ) b^{k+1}_p=J_p^{\rm T}(r_p+J_p(\xi^k+\Delta\xi-\xi_p))=b_p+H_p(\xi^{k+1}-\xi_p) bpk+1=JpT(rp+Jp(ξk+Δξ−ξp))=bp+Hp(ξk+1−ξp)
进一步的先验更新公式:
Δ ξ p k + 1 = − H p − 1 ( b p + H p ( ξ k + 1 − ξ p ) ) = − H p − 1 b p − ( ξ k + 1 − ξ p ) ) \Delta\xi_p^{k+1}=-H_p^{-1}(b_p+H_p(\xi^{k+1}-\xi_p))\\ =-H_p^{-1}b_p-(\xi^{k+1}-\xi_p)) Δξpk+1=−Hp−1(bp+Hp(ξk+1−ξp))=−Hp−1bp−(ξk+1−ξp))
先验约束公式2
上述约束可以进一步转化为:
e p ′ = ∣ ∣ ξ − ( ξ p − H p − 1 b p ) ∣ ∣ (3) e_p'=||\xi -(\xi_p-H_p^{-1}b_p)||\tag3 ep′=∣∣ξ−(ξp−Hp−1bp)∣∣(3)
其中, b p = J p T r p b_p=J_p^{\rm T}r_p bp=JpTrp。此时,该约束方程对应的雅可比矩阵等等如下:
J p ′ = I , H p ′ = I , b p ′ = J p ′ T [ ξ − ( ξ p − H p − 1 b p ) ] = ξ − ( ξ p − H p − 1 b p ) J_p'=I, \ \ H_p'=I,\ \ b_p' =J_p'^{\rm T}[\xi -(\xi_p-H_p^{-1}b_p)]=\xi -(\xi_p-H_p^{-1}b_p) Jp′=I, Hp′=I, bp′=Jp′T[ξ−(ξp−Hp−1bp)]=ξ−(ξp−Hp−1bp)
同样地,marg后第一次迭代, ξ = ξ p , b p ′ = H p − 1 b p \xi =\xi_p,\ \ b_p'=H_p^{-1}b_p ξ=ξp, bp′=Hp−1bp
Δ ξ p = − I ⋅ b p ′ = − H p − 1 b p \Delta \xi_p =-I\cdot b_p'= -H_p^{-1}b_p Δξp=−I⋅bp′=−Hp−1bp
后续随着变量 ξ \xi ξ的更新,Hessian阵 H p ′ H_p' Hp′为单位矩阵,不参与优化。偏置 b p ′ b_p' bp′的更新公式如下:
Δ b p ′ = Δ ξ \Delta b_p'=\Delta\xi Δbp′=Δξ
进一步的先验更新公式:
Δ ξ p k + 1 ′ = − I ( ξ k + 1 − ( ξ p − H p − 1 b p ) ) = − H p − 1 b p − ( ξ k + 1 − ξ p ) ) \Delta\xi_p^{k+1'}=-I(\xi^{k+1}-(\xi_p-H_p^{-1}b_p))\\ =-H_p^{-1}b_p-(\xi^{k+1}-\xi_p)) Δξpk+1′=−I(ξk+1−(ξp−Hp−1bp))=−Hp−1bp−(ξk+1−ξp))
从而我们得到同样的变量更新结果。从而可知这两个先验约束方程(2)和(3)等价。
由(3)我们可以很清楚地看出,marg后的先验约束,即约束变量维持 ξ = ( ξ p + H p − 1 b p ) \xi =(\xi_p+H_p^{-1}b_p) ξ=(ξp+Hp−1bp),即将marg之前算得的变量值,作为一种约束。在加入后续新的观测时,约束变量在marg之前算出的结果附近。
marg的问题及FEJ
边缘化带来的两个挑战:
-
Hessian矩阵变稠密
首先。marg会将Hessian矩阵变的稠密,而信息矩阵的稀疏性,使得使用回代法求解变量更新值时会带来一定的便利,如方法SAM。 -
变量更新的线性化点不一致的问题
marg后得到的先验信息矩阵包含了已经被marg掉的变量节点与窗口现有节点之间的约束关系,当窗口现有节点更新时,原本先验信息矩阵与现有节点相关的约束也应该进行更新。
而,使用marg在迭代优化算法中,marg后得到Hessian矩阵在新的迭代过程中不能再被优化,即图中的紫蓝色部分固定。而新的观测对应信息矩阵,红色部分,会在后续的迭代过程中进一步被优化。从而出现部分变量线性化点不一致的情况,从而使得解空间发生变化。
如本例中将变量 ξ 1 \xi_1 ξ1边缘化,而窗口中仍存在与 ξ 1 \xi_1 ξ1相关的其他变量(本例中为 ξ 2 , ξ 3 , ξ 5 , ξ 5 \xi_2,\xi_3,\xi_5,\xi_5 ξ2,ξ3,ξ5,ξ5),如下图所示。
• 红色为被 marg 变量以及测量约束。
• 绿色为跟 marg 变量有关的保留变量。
• 蓝色为和 marg 变量无关联的变量。
当 ξ 2 , ξ 3 , ξ 5 , ξ 5 \xi_2,\xi_3,\xi_5,\xi_5 ξ2,ξ3,ξ5,ξ5这些变量被进一步优化时,图中原本它们与 ξ 1 \xi_1 ξ1构成的约束(图中红线)对应的数据对应的线性化点也应该进行更新,但是由于边缘化后,Hessian矩阵的更新变的非常复杂。而新的变量 ξ 7 \xi_7 ξ7带来的与 ξ 2 \xi_2 ξ2构建的约束能够随着 ξ 2 \xi_2 ξ2的更新而更新。从而使得先验对于 ξ 2 \xi_2 ξ2的求导位置与新的约束对于 ξ 2 \xi_2 ξ2的求导位置不一致。从而使得 ξ 2 \xi_2 ξ2的解偏离真正的解。
这点也可以,从约束方程(3)中进行具体理解,即先验约束即使在得到新的约束方程后,会约束变量保持先验得到的变量结果,而所有变量的估计是不可能不存在误差的,因而先验得到的变量值也存在误差,进而进一步影响最后的变量估计,这个误差是随着时间累积的,即使有更多的观测约束,也消弭不了。最终滑动窗口+边缘化会表现出随着时间累积而误差累积的情况。
解决方案:FEJ(First Estimate Jacobians)强制所有涉及 ξ m \xi_m ξm 的雅可比在首次线性化点 ξ 0 \xi_0 ξ0 处计算。但此方法需要 ξ 0 \xi_0 ξ0本身接近于真实值(或者与后续迭代过程中的实际变量值(本该的线性化点)差距不大),否则会影响迭代速度,甚至解算失误的问题。从而该方法也会由此引入新累积误差。
实例分析:VINS-Mono中的滑动窗口策略
- 窗口结构:包含10个关键帧及其关联的路标点、IMU预积分。
- 边缘化规则:
- 当新帧加入时,若窗口已满,边缘化最旧帧及其观测的路标点。
- 若最旧帧为关键帧,保留其观测的路标点;否则直接丢弃。
- FEJ应用:所有被边缘化的路标点的雅可比矩阵固定在首次观测时的位姿处计算。
边缘化的代码实现(伪代码)
// 滑动窗口管理
void SlidingWindow::addNewFrame(Frame& new_frame) {window.push_back(new_frame);if (window.size() > WINDOW_SIZE) {Frame old_frame = window.front();window.pop_front();marginalize(old_frame); // 执行边缘化}
}// 边缘化操作
void marginalize(Frame& old_frame) {// 1. 构建Hessian矩阵和误差向量MatrixXd H;VectorXd b;buildHessian(old_frame, H, b);// 2. 分块:保留变量ξ_r 和边缘化变量ξ_mMatrixXd H_rr = H.block(r_indices, r_indices);MatrixXd H_rm = H.block(r_indices, m_indices);MatrixXd H_mm = H.block(m_indices, m_indices);VectorXd b_r = b.segment(r_indices);VectorXd b_m = b.segment(m_indices);// 3. 舒尔补计算先验信息MatrixXd H_prior = H_rr - H_rm * H_mm.inverse() * H_rm.transpose();VectorXd b_prior = b_r - H_rm * H_mm.inverse() * b_m;// 4. 将先验信息添加到后续优化问题optimizer.addPriorConstraint(H_prior, b_prior);
}
- 滑动窗口:通过动态管理状态变量控制计算复杂度。
- 边缘化:将移除变量的影响转化为先验约束,避免信息丢失。
- 舒尔补:数学工具实现信息矩阵的压缩。
- FEJ:解决线性化点不一致,减少误差累积。
参考
《手写VIO 贺一家,高翔,崔华坤》