VINS-Mono论文知识点精炼(翻译)

论文地址: VINS-Mono: A Robust and Versatile Monocular Visual-Inertial State Estimator

摘要

  • VINS-Mono:一种鲁棒的和多功能的单目视觉惯性状态估计器
  • 鲁棒的估计器初始化
  • 紧耦合非线性优化方法融合IMU测量值和特征观测值,获取高精度视觉惯性里程计
  • 使用低功耗回环检测模块
  • 四自由度全局位姿图优化
  • 地图可保存、加载和拼接

关键字: 单目视觉惯性系统、状态估计、传感器融合、同时定位建图

引言

  • 单目
    优点: 尺寸小、低功耗、安装简易.
    缺点: 无法获取尺度信息
  • 单目+IMU
    优点: 获取尺度信息、部分姿态信息(roll, pitch)、提升跟踪性能(当光照变换、纹理缺失或运动模糊时)
    问题:(1). 需要严谨的初始化: 视觉惯性系统是一个高度非线性系统,初始化很具挑战性;
            (2). 启动有限制: 大多数视觉惯性系统需要在指定的地方启动,并且需要缓慢小心移动;
            (3). 漂移: 已经存在基于回环检测、重定位和全局优化方法消除漂移;
            (4). 地图的重用: 需要有能够对地图进行保存和加载的方法.
  • VINS-Mono特点
    (1). 鲁棒的初始化程序能够从未知的初始状态引导系统;
    (2). 紧耦合的基于优化的单目VIO,带有相机IMU外参标定和IMU偏置矫正;
    (3). 在线重定位和四自由度全局位姿图优化;
    (4). 位姿图重用: 保存、加载和融合多个局部位姿图.

概述

在这里插入图片描述
从上面的框架图中可以获得VINS-Mono的组成如下:

  • 测量的预处理
    (1) 图像特征提取和跟踪
    (2) 预积分两图像帧之间的IMU测量值
  • 初始化
    提供引导基于非线性优化VIO系统所需的值,包括位姿、速度、重力向量、陀螺仪偏置和三维特征位置
  • 基于优化的VIO
  • 重定位
  • 位姿图优化
    (1) 利用经过几何一致性验证后的重定位结果,执行全局位姿图优化,消除漂移
    (2) 位姿图的重用

测量的预处理

A.视觉预处理

  • 特征跟踪: 使用KLT稀疏光流算法
  • 特征检测: (1) 检测新特征,维护每张图像最小特征数量(100~300)
                    (2) 设定相邻特征之间的最小距离,保证特征分布均匀
  • 离群特征丢弃: 使用基于基础矩阵的RANSAC方法
  • 特征生成: 去畸变,然后投影到单位球面上
  • 关键帧生成: 满足足够视差或跟踪质量太差任一条件即可

B.IMU预积分

  • IMU噪声和偏置
    陀螺仪和加速度计测量值:
    a ^ t = a t + b a t + R w t g w + n a ω ^ t = ω t + b w t + n w (1) \begin{aligned} &\hat{\mathbf{a}}_{t} = \mathbf{a}_{t} + \mathbf{b}_{a_{t}} + \mathbf{R}_{w}^{t}\mathbf{g}^{w} + \mathbf{n}_{a} \\ &\hat{\boldsymbol{\omega}}_{t} = \boldsymbol{\omega}_{t} + \mathbf{b}_{w_{t}} + \mathbf{n}_{w} \tag{1} \end{aligned} a^t=at+bat+Rwtgw+naω^t=ωt+bwt+nw(1)
    陀螺仪和加速度测量噪声为高斯白噪声: n a ∼ N ( 0 , σ a 2 ) , n w ∼ ( 0 , σ w 2 ) \mathbf{n}_{a}\sim\mathbf{N}(0,\boldsymbol{\sigma}_{a}^{2}),\mathbf{n}_{w}\sim(0,\boldsymbol{\sigma}_{w}^{2}) naN(0,σa2),nw(0,σw2)
    陀螺仪和加速度计被建模成随机游走模型,满足: b ˙ a t = n b a , b ˙ w t = n b w \dot{\mathbf{b}}_{a_{t}} = \mathbf{n}_{b_{a}},\dot{\mathbf{b}}_{w_{t}} = \mathbf{n}_{b_{w}} b˙at=nba,b˙wt=nbw,其中 n b a ∼ N ( 0 , σ b a 2 ) , n b w ∼ ( 0 , σ b w 2 ) \mathbf{n}_{b_{a}}\sim\mathbf{N}(0,\boldsymbol{\sigma}_{b_{a}}^{2}),\mathbf{n}_{b_{w}}\sim(0,\boldsymbol{\sigma}_{b_{w}}^{2}) nbaN(0,σba2),nbw(0,σbw2).
  • 预积分(PVQ增量)
    α b k + 1 b k = ∬ t ∈ [ t k , t k + 1 ] R t b k ( a ^ t − b a t ) d t 2 β b k + 1 b k = ∫ t ∈ [ t k , t k + 1 ] R t b k ( a ^ t − b a t ) d t γ b k + 1 b k = ∫ t ∈ [ t k , t k + 1 ] 1 2 Ω ( ω ^ t − b w t ) γ t b k d t (3) \begin{aligned} &\boldsymbol{\alpha}_{b_{k+1}}^{b_{k}} = \iint_{t\in[t_{k},t_{k+1}]}\mathbf{R}_{t}^{b_{k}}(\hat{\mathbf{a}}_{t} - \mathbf{b}_{a_{t}})dt^{2} \\ &\boldsymbol{\beta}_{b_{k+1}}^{b_{k}} = \int_{t\in[t_{k},t_{k+1}]}\mathbf{R}_{t}^{b_{k}}(\hat{\mathbf{a}}_{t} - \mathbf{b}_{a_{t}})dt \\ &\boldsymbol{\gamma}_{b_{k+1}^{b_{k}}} = \int_{t\in[t_{k},t_{k+1}]}\frac{1}{2}\boldsymbol{\Omega}(\hat{\boldsymbol{\omega}}_{t}-\mathbf{b}_{w_{t}})\boldsymbol{\gamma}_{t}^{b_{k}} dt \tag{3} \end{aligned} αbk+1bk=t[tk,tk+1]Rtbk(a^tbat)dt2βbk+1bk=t[tk,tk+1]Rtbk(a^tbat)dtγbk+1bk=t[tk,tk+1]21Ω(ω^tbwt)γtbkdt(3)
    其中:
    Ω ( ω ) = [ − ⌊ ω ⌋ × ω − ω T 0 ] , ⌊ ω ⌋ × = [ 0 − ω z ω y ω z 0 − ω x − ω y ω x 0 ] (4) \boldsymbol{\Omega}(\boldsymbol{\omega}) = \begin{bmatrix} -\left \lfloor \boldsymbol{\omega} \right \rfloor_{\times} & \boldsymbol{\omega} \\ -\boldsymbol{\omega}^{T} & 0 \end{bmatrix}, \left\lfloor \boldsymbol{\omega} \right\rfloor_{\times}= \begin{bmatrix} 0 & -\boldsymbol{\omega}_{z} & \boldsymbol{\omega}_{y} \\ \boldsymbol{\omega}_{z} & 0 & -\boldsymbol{\omega}_{x} \\ -\boldsymbol{\omega}_{y} & \boldsymbol{\omega}_{x} & 0 \tag{4} \end{bmatrix} Ω(ω)=[ω×ωTω0],ω×=0ωzωyωz0ωxωyωx0(4)
    从预积分增量的表达式中可以看出预积分增量是用 b k b_{k} bk作为参考坐标系,因此其它时刻的状态发生改变时不需要对PVQ增量进行重新积分.但表达式显示该积分增量与陀螺仪和加速度计偏置有关,也就意味着如果在优化过程中出现陀螺仪和加速度计的偏置发生改变,需要在优化的过程中不断进行上述积分过程,这是个巨大的硬件资源消耗,因此下面提出了一种近似方法.
  • 偏置矫正
    如果偏置 b a t , b w t \mathbf{b}_{a_{t}},\mathbf{b}_{w_{t}} bat,bwt发生小的改变,则按下式进行矫正:
    α b k + 1 b k = α ^ b k + 1 b k + J b a α δ b a k + J b w α δ b w β b k + 1 b k = β ^ b k + 1 b k + J b a β δ b a k + J b w β δ b w γ b k + 1 b k = γ ^ b k + 1 b k ⊗ [ 1 1 2 J b w γ δ b w ] (5) \begin{aligned} &\boldsymbol{\alpha}_{b_{k+1}}^{b_{k}} = \hat{\boldsymbol{\alpha}}_{b_{k+1}}^{b_{k}} + \mathbf{J}_{b_{a}}^{\alpha}\delta\mathbf{b}_{a_{k}} + \mathbf{J}_{b_{w}}^{\alpha}\delta\mathbf{b}_{w} \\ &\boldsymbol{\beta}_{b_{k+1}}^{b_{k}} = \hat{\boldsymbol{\beta}}_{b_{k+1}}^{b_{k}} + \mathbf{J}_{b_{a}}^{\beta}\delta\mathbf{b}_{a_{k}} + \mathbf{J}_{b_{w}}^{\beta}\delta\mathbf{b}_{w} \\ &\boldsymbol{\gamma}_{b_{k+1}^{b_{k}}} = \hat{\boldsymbol{\gamma}}_{b_{k+1}^{b_{k}}} \otimes \begin{bmatrix} 1\\ \frac{1}{2}\mathbf{J}_{b_{w}}^{\gamma}\delta\mathbf{b}_{w} \end{bmatrix} \tag{5} \end{aligned} αbk+1bk=α^bk+1bk+Jbaαδbak+Jbwαδbwβbk+1bk=β^bk+1bk+Jbaβδbak+Jbwβδbwγbk+1bk=γ^bk+1bk[121Jbwγδbw](5)
    否则的话,由于偏置有明显的变化,我们需要重新利用的新的偏置对PVQ增量预积分公式进行重新积分.该策略能够有效的避免了在优化过程中频繁的进行积分,降低了对硬件资源的使用.

初始化

A.滑窗内纯视觉SfM

我们为了限制计算复杂性,我们只维护了滑窗内的数据.首先我们检查最新帧与先前所有帧之间的相关性,如果在最新帧和滑窗内其它的任意一帧之间存在稳定的特征跟踪(超过30个特征点)并且有足够的视差(超过20个像素),则使用5点法恢复这两帧之间的相对旋转和缩放过的移动.然后,我们随机设置一个尺度因子,三角化这两帧之间的所有特征点.基于这些三角化过的特征点,使用PnP方法估计窗口内所有其它帧的位姿.最后,使用一个全局BA方法最小化所有观察到的特征的总的重投影误差.设置第一个相机帧坐标系 ( ⋅ ) c 0 (\cdot)^{c_{0}} ()c0作为SfM的参考坐标系,则所有其它帧的位姿定义为 ( p ˉ c k c 0 , q c k c 0 ) (\mathbf{\bar{p}}_{c_{k}}^{c_{0}},\mathbf{q}_{c_{k}}^{c_{0}}) (pˉckc0,qckc0).给定相机与IMU之间的外参 ( p c b , q c b ) (\mathbf{p}_{c}^{b},\mathbf{q}_{c}^{b}) (pcb,qcb),我们可以将相机位姿从相机坐标系转换到基座标系(IMU)下:
q b k c 0 = q c k c 0 ⊗ ( q c b ) − 1 s p b k c 0 = s p c k c 0 − R b k c 0 p c b (6) \begin{aligned} \mathbf{q}_{b_{k}}^{c_{0}} &= \mathbf{q}_{c_{k}}^{c_{0}}\otimes(\mathbf{q}_{c}^{b})^{-1} \\ s\mathbf{p}_{b_{k}}^{c_{0}} & = s\mathbf{p}_{c_{k}}^{c_{0}} - \mathbf{R}_{b_{k}}^{c_{0}}\mathbf{p}_{c}^{b} \tag{6} \end{aligned} qbkc0spbkc0=qckc0(qcb)1=spckc0Rbkc0pcb(6)
式中 s s s为尺度因子(后面会计算), R b k c 0 \mathbf{R}_{b_{k}}^{c_{0}} Rbkc0 q b k c 0 \mathbf{q}_{b_{k}}^{c_{0}} qbkc0转换成的旋转矩阵.

B.视觉惯性对齐

思路: 对齐视觉SfM结构和IMU预积分结果.

  • 陀螺仪偏置标定
    通过SfM可获得第k帧和第k+1帧姿态信息: q b k c 0 , q b k + 1 c 0 \mathbf{q}_{b_{k}}^{c_{0}},\mathbf{q}_{b_{k+1}}^{c_{0}} qbkc0,qbk+1c0
    通过预积分可获第k帧和第k+1帧之间姿态信息的变化量: γ b k + 1 b k \boldsymbol{\gamma}_{b_{k+1}}^{b_{k}} γbk+1bk
    于是可通过上述信息建立相邻两帧之间的约束,构建残差方程:
    min ⁡ δ b w ∑ k ∈ B ∥ q b k + 1 c 0 − 1 ⊗ q b k c 0 ⊗ γ b k + 1 b k ∥ 2 γ b k + 1 b k ≈ γ ^ b k + 1 b k ⊗ [ 1 1 2 J b w γ δ b w ] (7) \begin{aligned} &\min_{\delta b_{w}} \sum_{k\in B} \left\| {\mathbf{q}_{b_{k+1}}^{c_{0}}}^{-1}\otimes\mathbf{q}_{b_{k}}^{c_{0}}\otimes\boldsymbol{\gamma}_{b_{k+1}}^{b_{k}} \right\|^{2} \\ &\boldsymbol{\gamma}_{b_{k+1}}^{b_{k}} \approx \hat{\boldsymbol{\gamma}}_{b_{k+1}}^{b_{k}} \otimes \begin{bmatrix} 1 \\ \frac{1}{2}\mathbf{J}^{\gamma}_{b_{w}}\delta \mathbf{b}_{w} \end{bmatrix} \end{aligned} \tag{7} δbwminkBqbk+1c01qbkc0γbk+1bk2γbk+1bkγ^bk+1bk[121Jbwγδbw](7)
    式中 B B B为滑动窗口内的所有帧索引集合.通过最小化公式7,即可获得陀螺仪偏置的增量,即实现陀螺仪偏置标定.在标定完陀螺仪偏置之后,需要重新利用新的偏置对IMU测量值进行预积分.
  • 速度、重力向量和尺度因子的初始化
    待初始化变量: χ I = [ v b 0 b 0 , v b 1 b 1 , . . . , v b n b n , g c 0 , s ] \chi_{I} = [\mathbf{v}_{b_{0}}^{b_{0}}, \mathbf{v}_{b_{1}}^{b_{1}}, ..., \mathbf{v}_{b_{n}}^{b_{n}}, \mathbf{g}^{c_{0}}, s] χI=[vb0b0,vb1b1,...,vbnbn,gc0,s]
    其中 v b k b k \mathbf{v}_{b_{k}}^{b_{k}} vbkbk是采集第 k k k张图像时基体坐标系下的速度, g c 0 \mathbf{g}^{c_{0}} gc0是在 c 0 c_{0} c0坐标系下的重力向量, s s s是单目SfM无法获取的尺度因子.
    考虑滑动窗口内的相邻两帧之间的关系,可建立如下PV增量和两帧状态的关系:
    α b k + 1 b k = R c 0 b k ( s ( p ˉ b k + 1 c 0 − p ˉ b k c 0 ) + 1 2 g c 0 Δ t k 2 − R b k c 0 v b k b k Δ t k ) β b k + 1 b k = R c 0 b k ( R b k + 1 c 0 v b k + 1 b k + 1 + g c 0 Δ t − R b k c 0 v b k b k ) (9) \begin{aligned} &\boldsymbol{\alpha}_{b_{k+1}}^{b_{k}} = \mathbf{R}_{c_{0}}^{b_{k}}(s(\bar{\mathbf{p}}_{b_{k+1}}^{c_{0}}-\bar{\mathbf{p}}_{b_{k}}^{c_{0}}) + \frac{1}{2}\mathbf{g}^{c_{0}}\Delta t_{k}^{2} - \mathbf{R}_{b_{k}}^{c_{0}}\mathbf{v}_{b_{k}}^{b_{k}}\Delta t_{k}) \\ &\boldsymbol{\beta}_{b_{k+1}}^{b_{k}} = \mathbf{R}_{c_{0}}^{b_{k}}(\mathbf{R}_{b_{k+1}}^{c_{0}}\mathbf{v}_{b_{k+1}}^{b_{k+1}}+\mathbf{g}^{c_{0}}\Delta t - \mathbf{R}_{b_{k}}^{c_{0}}\mathbf{v}_{b_{k}}^{b_{k}}) \tag{9} \end{aligned} αbk+1bk=Rc0bk(s(pˉbk+1c0pˉbkc0)+21gc0Δtk2Rbkc0vbkbkΔtk)βbk+1bk=Rc0bk(Rbk+1c0vbk+1bk+1+gc0ΔtRbkc0vbkbk)(9)
    将公式6带入9中,然后建立如下线性测量模型:
    z ^ b k + 1 b k = [ α ^ b k + 1 b k − p c b + R c 0 b k R b k + 1 c 0 p c b β ^ b k + 1 b k ] = H b k + 1 b k χ I + n b k + 1 b k (10) \hat{\mathbf{z}}_{b_{k+1}}^{b_{k}}= \begin{bmatrix} \hat{\boldsymbol{\alpha}}_{b_{k+1}}^{b_{k}} - \mathbf{p}_{c}^{b}+\mathbf{R}_{c_{0}}^{b_{k}}\mathbf{R}_{b_{k+1}}^{c_{0}}\mathbf{p}_{c}^{b} \\ \hat{\boldsymbol{\beta}}_{b_{k+1}}^{b_{k}} \end{bmatrix} = \mathbf{H}_{b_{k+1}}^{b_{k}}\chi_{I}+\mathbf{n}_{b_{k+1}}^{b_{k}} \tag{10} z^bk+1bk=[α^bk+1bkpcb+Rc0bkRbk+1c0pcbβ^bk+1bk]=Hbk+1bkχI+nbk+1bk(10)
    其中:
    H b k + 1 b k = [ − I Δ t 0 1 2 R c 0 b k Δ t k 2 R c 0 b k ( p ˉ c k + 1 c 0 − p ˉ c k c 0 ) − I R c 0 b k R b k + 1 c 0 R c 0 b k Δ t k 0 ] (11) \mathbf{H}_{b_{k+1}}^{b_{k}}= \begin{bmatrix} -\mathbf{I}\Delta t & \mathbf{0} & \frac{1}{2}\mathbf{R}_{c_{0}}^{b_{k}}\Delta t_{k}^{2} & \mathbf{R}_{c_{0}}^{b_{k}}(\bar{\mathbf{p}}_{c_{k+1}}^{c_{0}}-\mathbf{\bar{\mathbf{p}}}_{c_{k}}^{c_{0}}) \\ -\mathbf{I} & \mathbf{R}_{c_{0}}^{b_{k}}\mathbf{R}_{b_{k+1}}^{c_{0}} & \mathbf{R}_{c_{0}}^{b_{k}}\Delta t_{k} & \mathbf{0} \end{bmatrix} \tag{11} Hbk+1bk=[IΔtI0Rc0bkRbk+1c021Rc0bkΔtk2Rc0bkΔtkRc0bk(pˉck+1c0pˉckc0)0](11)
    式中 R b k c 0 , R b k + 1 c 0 , p ˉ c k c 0 , p ˉ c k + 1 c 0 \mathbf{R}_{b_{k}}^{c_{0}},\mathbf{R}_{b_{k+1}}^{c_{0}},\bar{\mathbf{p}}_{c_{k}}^{c_{0}},\bar{\mathbf{p}}_{c_{k+1}}^{c_{0}} Rbkc0,Rbk+1c0,pˉckc0,pˉck+1c0可通过SfM过程获得, Δ t \Delta t Δt是相邻两帧之间的间隔.因此可通过构建如下最小二乘问题:
    min ⁡ χ I ∑ ∥ z ^ b k + 1 b k − H b k + 1 b k χ I ∥ 2 (12) \min_{\chi_{I}}\sum\left\| \hat{\mathbf{z}}_{b_{k+1}}^{b_{k}} - \mathbf{H}_{b_{k+1}}^{b_{k}}\chi_{I}\right\|^{2} \tag{12} χIminz^bk+1bkHbk+1bkχI2(12)
    来优化求解出窗口中每一帧的速度、重力向量和尺度因子.
  • 重力细化
    大多数情况下重力向量的大小是已知的,因此重力向量只剩下2个自由度了.本文将重力向量按下图进行定义, g = g ( g ^ ⃗ + δ g ) , δ g = w 1 b 1 + w 2 b 2 \mathbf{g}=g(\vec{\hat{\mathbf{g}}}+\delta \mathbf{g}), \delta \mathbf{g}=w_{1}\mathbf{b}_{1}+w_{2}\mathbf{b}_{2} g=g(g^ +δg),δg=w1b1+w2b2,其中 g g g是已知的重力大小, g ⃗ \vec{g} g 是代表重力方向的单位向量, b 1 , b 2 \mathbf{b}_{1},\mathbf{b}_{2} b1,b2是重力向量切平面上的两个互相垂直的基向量.可看出该定义保证了重力向量在大小已知的情况下只有两个自由度.
    在这里插入图片描述

我们将新的重力向量 g \mathbf{g} g带入公司9,然后重新进行公式12的优化求解,即可对重力向量进行细化.

紧耦合单目VIO

在初始化完成之后,我们进行基于滑窗的紧耦合单目VIO算法获取高精度和鲁棒性的状态估计.

A. 公式

  • 滑窗内完整状态向量定义
    χ = [ x 0 , x 1 , . . . x n , x c b , λ 0 , λ 1 , . . . λ m ] x k = [ p b k w , v b k w , q b k w , b a , b g ] , k ∈ [ 0 , n ] x c b = [ p c b , q c b ] (13) \begin{aligned} \boldsymbol{\chi} &= [\mathbf{x}_{0}, \mathbf{x}_{1},...\mathbf{x}_{n},\mathbf{x}_{c}^{b},\lambda_{0},\lambda_{1},...\lambda_{m}] \\ \mathbf{x}_{k} &= [\mathbf{p}_{b_{k}}^{w}, \mathbf{v}_{b_{k}}^{w}, \mathbf{q}_{b_{k}}^{w}, \mathbf{b}_{a},\mathbf{b_{g}}], k\in[0,n]\\ \mathbf{x}_{c}^{b} &= [\mathbf{p}_{c}^{b},\mathbf{q}_{c}^{b}] \end{aligned} \tag{13} χxkxcb=[x0,x1,...xn,xcb,λ0,λ1,...λm]=[pbkw,vbkw,qbkw,ba,bg],k[0,n]=[pcb,qcb](13)
    其中 x k \mathbf{x}_{k} xk是采集第 k k k张图像时的IMU状态,由PVQ状态量和IMU偏置组成. λ l \lambda_{l} λl是第 l l l个特征的逆深度(此深度是特征在第一次出现时定义的深度). x c b \mathbf{x}_{c}^{b} xcb是相机与IMU之间的外参.
  • 视觉惯性BA公式
    min ⁡ χ { ∥ r p − H p χ ∥ 2 + ∑ k ∈ B ∥ r B ( z ^ b k + 1 b k , χ ) ∥ P b k + 1 b k 2 + ∑ ( l , j ) ∈ C ρ ( ∥ r C ( z ^ l c j , χ ) ∥ P l c j 2 ) } (14) \min_{\boldsymbol{\chi}} \left\{ {\left\| \mathbf{r}_{p}-\mathbf{H}_{p}\boldsymbol{\chi} \right\|}^{2} + {\sum_{k\in B}\left\| \mathbf{r}_{B}(\hat{\mathbf{z}}_{b_{k+1}}^{b_{k}},\boldsymbol{\chi}) \right\|}_{\mathbf{P}_{b_{k+1}}^{b_{k}}}^{2} + \sum_{(l,j)\in C}\rho(\left\| \mathbf{r}_{C}(\hat{\mathbf{z}}_{l}^{c_{j}},\boldsymbol{\chi}) \right\|^{2}_{\mathbf{P}_{l}^{c_{j}}}) \right\} \tag{14} χminrpHpχ2+kBrB(z^bk+1bk,χ)Pbk+1bk2+(l,j)Cρ(rC(z^lcj,χ)Plcj2)(14)
    式中 ρ ( . ) \rho(.) ρ(.)为Huber norm损失函数,为了抑制噪声的影响,定义如下:
    ρ ( s ) = { s 2 s − 1 s ≤ 1 s > 1 (15) \rho(s) = \left\{\begin{matrix} s \\ 2\sqrt{s}-1 \end{matrix}\right. \begin{matrix} s\leq1 \\ s>1 \end{matrix} \tag{15} ρ(s)={s2s 1s1s>1(15)
    式中 r B , r C \mathbf{r}_{B},\mathbf{r}_{C} rB,rC分别是对于IMU和视觉测量的残差项,后面会详细介绍. [ r p , H p ] [\mathbf{r}_{p},\mathbf{H}_{p}] [rp,Hp]是来自于边缘化过程中的先验信息,后面也有介绍.该视觉惯性BA公式的求解是个最小二乘优化问题,使用Ceres-solver优化库进行求解.

B. IMU测量残差

以滑窗内两个连续帧 b k , b k + 1 b_{k},b_{k+1} bk,bk+1进行说明,可对IMU测量残差进行如下定义:
r B = [ δ α b k + 1 b k δ β b k + 1 b k δ θ b k + 1 b k δ b a δ b ] = [ R w b k ( p b k + 1 w − p b k w + 1 2 g w Δ t 2 − v b k w Δ t ) − α ^ b k + 1 b k R w b k ( v b k + 1 w + g w Δ t − v b k w ) − β ^ b k + 1 b k 2 [ q b k w − 1 ⊗ q b k + 1 w ⊗ ( γ b k + 1 b k ) − 1 ] x y z b a b k + 1 − b a b k b w b k + 1 − b w b k ] (16) \mathbf{r}_{B} = \begin{bmatrix} \delta \boldsymbol{\alpha}_{b_{k+1}}^{b_{k}} \\ \delta \boldsymbol{\beta}_{b_{k+1}}^{b_{k}} \\ \delta \boldsymbol{\theta}_{b_{k+1}}^{b_{k}} \\ \delta \mathbf{b}_{a} \\ \delta \mathbf{b}_{} \end{bmatrix}= \begin{bmatrix} \mathbf{R}_{w}^{b_{k}}(\mathbf{p}_{b_{k+1}}^{w} - \mathbf{p}_{b_{k}}^{w} + \frac{1}{2}\mathbf{g}^{w}\Delta t^{2} - \mathbf{v}_{b_{k}}^{w}\Delta t) - \hat{\boldsymbol{\alpha}}_{b_{k+1}}^{b_{k}} \\ \mathbf{R}^{b_{k}}_{w}(\mathbf{v}_{b_{k+1}}^{w} + \mathbf{g}^{w}\Delta t - \mathbf{v}_{b_{k}}^{w}) - \hat{\boldsymbol{\beta}}_{b_{k+1}}^{b_{k}} \\ 2[{\mathbf{q}_{b_{k}}^{w}}^{-1}\otimes\mathbf{q}_{b_{k+1}}^{w}\otimes(\boldsymbol{\gamma}_{b_{k+1}}^{b_{k}})^{-1}]_{xyz} \\ \mathbf{b}_{a_{b_{k+1}}} - \mathbf{b}_{a_{b_{k}}} \\ \mathbf{b}_{w_{b_{k+1}}} - \mathbf{b}_{w_{b_{k}}} \end{bmatrix} \tag{16} rB=δαbk+1bkδβbk+1bkδθbk+1bkδbaδb=Rwbk(pbk+1wpbkw+21gwΔt2vbkwΔt)α^bk+1bkRwbk(vbk+1w+gwΔtvbkw)β^bk+1bk2[qbkw1qbk+1w(γbk+1bk)1]xyzbabk+1babkbwbk+1bwbk(16)
式中 [ . ] x y z [.]_{xyz} [.]xyz指代的是将四元素转换成轴角表示, [ α ^ b k + 1 b k , β ^ b k + 1 b k , γ ^ b k + 1 b k ] [\hat{\boldsymbol{\alpha}}_{b_{k+1}}^{b_{k}},\hat{\boldsymbol{\beta}}_{b_{k+1}}^{b_{k}},\hat{\boldsymbol{\gamma}}_{b_{k+1}}^{b_{k}}] [α^bk+1bk,β^bk+1bk,γ^bk+1bk]是两个相邻图像帧之间的IMU测量值进行预积分的结果.

C.视觉测量残差

  • 残差理解
    与传统针孔相机定义重投影误差在一个常见的图像平面上不一样,我们定义相机测量残差在一个单位球上面.原因是因为几乎所有类型的摄像机(包括广角,鱼眼或全向摄像机)的光学元件都可以建模为连接单位球表面的单位射线。
  • 定义
    假设第 l l l个特征首先是在第 i i i张图像中被看到,则在第 j j j张图像中观察到特征 l l l的测量残差可定义如下:
    r C ( z ^ l c j , χ ) = [ b 1 , b 2 ] T ⋅ ( ϱ ˉ ^ l c j − ϱ l c j ∥ ϱ l c j ∥ ) ϱ ˉ ^ l c j = π c − 1 ( [ u ^ l c j v ^ l c j ] ) ϱ l c j = R b c ( R w b j ( R b i w ( R c b 1 λ l π c − 1 ( [ u ^ l c i v ^ l c i ] ) + p c b ) + p b i w − p b j w ) − p c b ) (17) \begin{aligned} & \mathbf{r}_{C}(\hat{\mathbf{z}}_{l}^{c_{j}},\boldsymbol{\chi}) = [\mathbf{b}_{1},\mathbf{b}_{2}]^{T}\cdot(\hat{\bar{\boldsymbol{\varrho}}}^{c_{j}}_{l}-\frac{\boldsymbol{\varrho}_{l}^{c_{j}}}{\left\| \boldsymbol{\varrho}_{l}^{c_{j}} \right\|}) \\ & \hat{\bar{\boldsymbol{\varrho}}}_{l}^{c_{j}} = \pi_{c}^{-1}(\begin{bmatrix} \hat{u}_{l}^{c_{j}} \\ \hat{v}_{l}^{c_{j}} \end{bmatrix}) \\ & \boldsymbol{\varrho}_{l}^{c_{j}} = \mathbf{R}_{b}^{c}(\mathbf{R}_{w}^{b_{j}}(\mathbf{R}_{b_{i}}^{w}(\mathbf{R}_{c}^{b}\frac{1}{\lambda}_{l}\pi_{c}^{-1}(\begin{bmatrix} \hat{u}_{l}^{c_{i}} \\ \hat{v}_{l}^{c_{i}} \end{bmatrix})+\mathbf{p}_{c}^{b})+\mathbf{p}_{b_{i}}^{w}-\mathbf{p}_{b_{j}}^{w})-\mathbf{p}_{c}^{b}) \end{aligned} \tag{17} rC(z^lcj,χ)=[b1,b2]T(ϱˉ^lcjϱlcjϱlcj)ϱˉ^lcj=πc1([u^lcjv^lcj])ϱlcj=Rbc(Rwbj(Rbiw(Rcbλ1lπc1([u^lciv^lci])+pcb)+pbiwpbjw)pcb)(17)
    式中 [ u ^ l c i , v ^ l c i ] [\hat{u}_{l}^{c_{i}},\hat{v}_{l}^{c_{i}}] [u^lci,v^lci]是第 l l l个特征在第 i i i张图像中第一次被观察到的测量值, [ u ^ l c j , v ^ l c j ] [\hat{u}_{l}^{c_{j}},\hat{v}_{l}^{c_{j}}] [u^lcj,v^lcj]是第 l l l个特征在第 j j j张图像中被观察到的测量值, π c − 1 \pi_{c}^{-1} πc1是反投影函数(将像素值通过内参投影到一个单位向量).由于视觉残差的自由度为2,于是我们按下图所示方式将残差向量投影到以 [ b 1 , b 2 ] [\mathbf{b}_{1},\mathbf{b}_{2}] [b1,b2]为基向量的切平面上,最后生成公式17所示的残差定义.
    在这里插入图片描述

D.边缘化

论文中写的不够详细,不如直接看大神贺一家的博客https://blog.csdn.net/heyijia0327/article/details/52822104

重定位

我们使用滑窗策略和边缘化策略方法限制了算法的复杂性,但会引入累积漂移问题.为了消除漂移,提出了一个融合了单目VIO的紧耦合重定位模块.该模块首先进行回环检测,判断是否观测到已经观测过的场景,然后建立回环帧与当前帧之间的特征关系.随后,这些特征关系将会被融入单目VIO模块,最后能够利用最小资源实现无漂移状态估计.多特征的多观测被直接用于重定位,能够提高精度和更好的状态估计顺滑性.上述过程如下图所示.
在这里插入图片描述

A.回环检测

我们利用前沿的词袋位置识别方法DBoW2进行回环检测.除了被用于单目VIO的角点特征外,另外500个角点也被检测,同时利用BRIEF描述器对角点进行描述.额外的角点特征被使用可以在回环检测时实现更好的召回率.描述子将会被作为单词在视觉数据集里面进行检索.DBoW2返回回环检测到的并经过时间和空间一致性检查的候选帧.我们维护所有的用于特征检索的描述子,但将原始图像丢弃,从而节约内存资源.

B.特征检索

当一个回环检测到,可以通过检索局部滑窗内的帧和回环候选帧之间的特征一致性建立两者之间的联系,其中一致性的检索是通过特征描述子匹配进行的.描述子匹配可能会引入误匹配,于是我们使用了两步几何离群值丢弃策略,减少误匹配.

  • 2-D-2-D: 基于基础矩阵的RANSAC方法去除误匹配
  • 3-D-2-D: 基于PnP方法的RANSAC方法去除误匹配

C.紧耦合重定位

重定位过程有效的对齐了当前滑窗和先前的位姿.在重定位过程中,我们将回环帧的位姿置为静态的,将滑窗内的所有IMU测量、局部视觉测量和检索到的特征一起进行优化.我们同公式(17)一样可以轻松的写出视觉测量模型残差,其中唯一的差异是回环帧的位姿 ( q ^ v w , p ^ v w ) (\hat{\mathbf{q}}_{v}^{w},\hat{\mathbf{p}}_{v}^{w}) (q^vw,p^vw)是直接来自于位姿图或过去的里程计输出,而且是静态的.最后,我们只需要在公式(14)的非线性代价函数中添加回环误差项构成新的误差函数,具体函数如下:
min ⁡ χ { ∥ r p − H p χ ∥ 2 + ∑ k ∈ B ∥ r B ( z ^ b k + 1 b k , χ ) ∥ P b k + 1 b k 2 + ∑ ( l , j ) ∈ C ρ ( ∥ r C ( z ^ l c j , χ ) ∥ P l c j 2 ) + ∑ ( l , v ) ∈ L ρ ( ∥ r C ( z ^ l v , χ , q ^ v w , p ^ v w ) ∥ P l c v 2 ) } \min_{\boldsymbol{\chi}} \left\{ {\left\| \mathbf{r}_{p}-\mathbf{H}_{p}\boldsymbol{\chi} \right\|}^{2} + {\sum_{k\in B}\left\| \mathbf{r}_{B}(\hat{\mathbf{z}}_{b_{k+1}}^{b_{k}},\boldsymbol{\chi}) \right\|}_{\mathbf{P}_{b_{k+1}}^{b_{k}}}^{2} + \sum_{(l,j)\in C}\rho(\left\| \mathbf{r}_{C}(\hat{\mathbf{z}}_{l}^{c_{j}},\boldsymbol{\chi}) \right\|^{2}_{\mathbf{P}_{l}^{c_{j}}}) + \sum_{(l,v)\in L}\rho(\left\| \mathbf{r}_{C}(\hat{\mathbf{z}}_{l}^{v},\boldsymbol{\chi},\hat{\mathbf{q}}_{v}^{w},\hat{\mathbf{p}}_{v}^{w}) \right\|^{2}_{\mathbf{P}_{l}^{c_{v}}}) \right\} χminrpHpχ2+kBrB(z^bk+1bk,χ)Pbk+1bk2+(l,j)Cρ(rC(z^lcj,χ)Plcj2)+(l,v)Lρ(rC(z^lv,χ,q^vw,p^vw)Plcv2)
式中 L L L是在回环帧中检索到的特征的集合. ( l , v ) (l,v) (l,v)指的是第 l l l个特征被第 v v v个回环帧观测到.注意,虽然代价函数与公式(14)有轻微的不同,但状态的维度并没有发生变化,因为回环帧的位姿是被置为静态的.当与当前滑窗建立了多个回环,我们优化将会使用所有的相关帧的相关特征.这将会给重定位过程引入多视图约束,结果将会更精确更顺滑.

全局位姿图优化和地图重用

A. 四自由度累积漂移方向

得益于重力的惯性测量,在世界坐标系中roll和pitch是绝对状态,而 x , y , z , y a w x,y,z,yaw x,y,z,yaw是相对与参考系的相对状态.因此累积漂移只发生在 x , y , z , y a w x,y,z,yaw x,y,z,yaw四个状态.为了充分利用这个有效信息来有效的矫正漂移,我们固定 r o l l , p i t c h roll,pitch roll,pitch值,只执行4自由度的位姿图优化

B. 添加关键帧到位姿图中

在VIO过程之后,关键帧将会被添加到位姿图中.每个关键帧作为位姿图中的一个节点,并且通过两种类型的边与其它节点相连.两种边类型如下:

  • 顺序边
    一个关键帧与它的前几帧建立几个顺序边.一个顺序边表示两个关键帧之间的相对转换,是直接通过VIO获得的.考虑关键帧 i i i和它的一个先前帧 j j j,顺序边只包含相对位置 p ^ i j i \hat{\mathbf{p}}_{ij}^{i} p^iji y a w yaw yaw ψ ^ i j \hat{\psi}_{ij} ψ^ij.
    p ^ i j i = R ^ i ( p ^ j w − p ^ i w ) ψ ^ i j = ψ ^ j − ψ ^ j (19) \begin{aligned} \hat{p}_{ij}^{i} &= \hat{\mathbf{R}}_{i}(\hat{\mathbf{p}}_{j}^{w} - \hat{\mathbf{p}}_{i}^{w}) \\ \hat{\psi}_{ij} &= \hat{\psi}_{j}-\hat{\psi}_{j} \end{aligned} \tag{19} p^ijiψ^ij=R^i(p^jwp^iw)=ψ^jψ^j(19)
  • 回环边
    如果关键帧有一个回环连接,它通过一个回环边连接回环帧.类似的,回环边只包含4自由度相对位姿变换,同公式(19)一致.回环边的值是通过重定位获得的.

c. 4自由度位姿图优化

我们定义第 i , j i,j i,j帧之间的边的残差如下:
r ( p i w , ψ i , p j w , ψ j ) = [ R ( ϕ ^ i , θ ^ i , ψ j ) − 1 ( p j w − p j w − p ^ i j i ) ψ j − ψ i − ψ ^ i j ] (20) \mathbf{r}(\mathbf{p}_{i}^{w},\psi_{i},\mathbf{p}_{j}^{w},\psi_{j})= \begin{bmatrix} \mathbf{R}(\hat{\phi}_{i},\hat{\theta}_{i},\psi_{j})^{-1}(\mathbf{p}_{j}^{w}-\mathbf{p}_{j}^{w}-\hat{\mathbf{p}}_{ij}^{i}) \\ \psi_{j}-\psi_{i}-\hat{\psi}_{ij} \end{bmatrix} \tag{20} r(piw,ψi,pjw,ψj)=[R(ϕ^i,θ^i,ψj)1(pjwpjwp^iji)ψjψiψ^ij](20)
式中 ϕ ^ i , θ ^ i \hat{\phi}_{i},\hat{\theta}_{i} ϕ^i,θ^i是固定的 r o l l , p i t c h roll,pitch roll,pitch角度估计值,是直接通过VIO获得的.
整个包含顺序边和回环边的位姿图是通过最小化下式代价函数获得的:
min ⁡ p , ψ { ∑ ( i , j ) ∈ S ∥ r i , j ∥ 2 + ∑ ( i , j ) ∈ L ρ ( ∥ r i , j ∥ 2 ) } (21) \min_{\mathbf{p},\psi} \left\{ \sum_{(i,j)\in S}\left\| \mathbf{r}_{i,j} \right\|^{2}+\sum_{(i,j)\in L}\rho(\left\| \mathbf{r}_{i,j} \right\|^{2}) \right\} \tag{21} p,ψmin(i,j)Sri,j2+(i,j)Lρ(ri,j2)(21)
式中 S S S是所有顺序边的集合, L L L是所有回环边的集合.虽然紧耦合重定位已经帮我们消除了错误的回环,我们任然添加Huber norm ρ ( ⋅ ) \rho(\cdot) ρ()进一步减少任意错误回环的影响.相反的,我们对顺序边使用任何损失函数,因为这些边是直接通过VIO获取的,已经包含该了足够多的离群值丢弃策略.
位姿图优化和重定位同步运行在两个不同的线程中,这保证了对最新优化的位姿图的快速使用.同样额,虽然当前位姿图优化没有完成,重定位任然可以使用已经存在的位姿图进行.

D. 位姿图融合

位姿图不仅能够优化当前的地图,而且也能够融合当前的地图和先前的地图.如果我们已经加载了一个先前构建的地图,并且在两幅地图中检测到了回环,我们就可以融合它们到一起.因为所有的边是相对约束,位姿图通过回环联系自动融合两幅地图.

E. 位姿图保存

我们的位姿图结构是非常简单的.我们只需要保存节点、边和每个关键帧的描述子.原始图像是丢弃的,为了减少内存的消耗.具体一点,我们保存的第 i i i关键帧是:
[ i , p ^ i w , q ^ i w , v , p ^ i v i , ψ ^ i v , D ( u , v , d e s ) ] (22) [i,\hat{\mathbf{p}}_{i}^{w},\hat{\mathbf{q}}_{i}^{w},v,\hat{\mathbf{p} }_{iv}^{i},\hat{\psi}_{iv},\mathbf{D}(u,v,des)] \tag{22} [i,p^iw,q^iw,v,p^ivi,ψ^iv,D(u,v,des)](22)
式中 i i i是帧索引, p ^ i w , q ^ i w \hat{\mathbf{p}}_{i}^{w},\hat{\mathbf{q}}_{i}^{w} p^iw,q^iw是关键帧的位置和姿态,通过VIO获得.如果这一帧有回环帧, v v v是回环帧的索引. p ^ i v i , ψ ^ i v \hat{\mathbf{p} }_{iv}^{i},\hat{\psi}_{iv} p^ivi,ψ^iv是这两帧之间的相对位置和 y a w yaw yaw角,通过重定位获得. D ( u , v , d e s ) \mathbf{D}(u,v,des) D(u,v,des)是特征集合,每个特征包含2自由度位置和它的描述子.

F. 位姿图加载

我们使用同保存格式相同的格式加载关键帧.每个关键帧是位姿图中的一个节点.节点的初始位姿是 p ^ i w , q ^ i w \hat{\mathbf{p}}_{i}^{w},\hat{\mathbf{q}}_{i}^{w} p^iw,q^iw.回环边直接通过回环信息 p ^ i v i , ψ ^ i v \hat{\mathbf{p} }_{iv}^{i},\hat{\psi}_{iv} p^ivi,ψ^iv建立.每个关键帧与它的相邻关键帧建立顺序边.在加载完位姿图之后,我们执行全局4自由度的位姿图优化.位姿图保存和加载的时间与位姿图的尺寸成正比.

  • 3
    点赞
  • 12
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值