Tightly Coupled 3D Lidar Inertial Odometry and Mapping

摘要
I. 介绍
II. 相关工作
III. 符号与说明
- A. 符号
- B. IMU运动
- - 1) 状态
  - 2) 运动模型
  - 3) 预积分
IV. 紧耦合Lidar-IMU里程计
- A. Lidar-IMU里程计概述
- B. 去畸变和特征提取
- C. 相对Lidar测量
- D.激光雷达扫描匹配
- E. 优化
V.使用旋转约束进行细化
VI. 实验
VII. 测试与分析

本文主要是对LIO-Mappingt论文的一些学习总结，具体的参考文献请参考原文，如有不对之处，希望大家提出宝贵的意见，共同进步~

摘要

运动估计是移动机器人的基本要求之一，通过传感器融合我们可以弥补单一传感器的缺陷并提供更可靠的运动估计。在这篇文章中，我们介绍一种Lidar-IMU紧耦合融合方法。提供通过联合最小化雷达和IMU测量的误差，雷达-IMU里程计（LIO）能在长时间的实验中，甚至在激光测量环境出现退化的情况下，也能保持较小的漂移。除此之外，为了获得更可靠的雷达位姿估计，此外，为了获得更可靠的位姿估计，提出了一种旋转约束的优化算法(LIO-mapping)，进一步将激光雷达位姿与全局地图对齐。实验结果表明，该方法能够在IMU高更新速率下，高精度地估计传感器的位姿，即使是在快速运动的条件下，或在特征不足的情况下。

I. 介绍

自我运动估计在许多导航任务中起着重要的作用，是自主机器人的关键问题之一。它提供了机器人的姿态信息，并可以为姿态控制提供即时反馈。此外，它与各种感知环境的传感器一起，为同步定位和建图(SLAM)提供关键信息。准确估计机器人的姿态有助于降低风险，并有助于成功的运动规划。

激光雷达传感器可以对周围环境进行距离测量，在机器人系统中得到了广泛的应用。典型的3D激光雷达可以在360度水平视场(FOV)下，以10Hz左右的频率感知周围环境。此外，作为一种主动传感器，它对光照具有不变性。高可靠性和高精度使激光雷达传感器成为姿态估计的热门选择。

尽管它有很多优点，但是激光雷达传感器并不完美，也有一些缺点。从激光雷达本身来看，其垂直分辨率较低，获得的稀疏点云提供的特征有限，使得特征跟踪成为一个棘手的问题。此外，安装在移动机器人上的激光雷达会产生运动畸变，直接影响传感精度。在现实场景中，有一些激光雷达退化的情况，在这种情况下，激光雷达接收的点云比较少。例如，3D激光雷达在狭窄的走廊环境中只接收到少量可用点。接收的点云主要来自侧墙，只有一小部分点云可以从天花板和地板上观察到。在这种情况下，匹配的激光雷达特征很容易导致不受约束的姿态估计。另一个缺点是它的更新率低。这限制了它在需要快速响应的任务上的应用，如机器人姿态的控制。

在本文中，我们提出了一种紧耦合的3D Lidar-IMU位姿估计算法来克服上述问题。激光雷达和IMU的测量数据用于联合优化。为了实现实时和一致的估计，对历史位姿采用滑动窗口和边缘化优化，然后进行旋转约束的优化。我们工作的主要贡献如下:

提出了一种紧耦合激光雷达- imu测程算法。它提供了实时准确的状态估计和较高的更新速度。
基于lidar-IMU里程计的先验，旋转约束精化方法进一步优化最终的姿态和生成的点云地图。即使在一些激光雷达退化的情况下，它也能确保一致和稳健的估计。
通过大量的室内和室外测试，验证了算法的有效性。它优于最先进的激光纯或松散耦合激光- imu算法。
源代码可在online1上获得。这是社区提供的第一个紧密耦合激光雷达和IMU融合的开源实现。

本文的其余部分组织如下：第二节对相关工作进行了回顾。第三节解释了符号和一些初步说明。所提出的里程计和精化方法分别在第四章和第五章。实现和测试见第六章和第七章。结论见第八章。

II. 相关工作

这一部分分别介绍了基于松耦合的Lidar-IMU以及基于紧耦合的Lidar-IMU融合算法的研究现状，并简单介绍了作者提出的融合算法。

受其他视觉惯性工作[10]，[11]的启发，我们设计了紧耦合Lidar-IMU融合的方法。我们“预积分”并使用原始IMU观测和激光雷达观测来优化整个系统的状态，这可以在激光退化或运动快速的情况下工作。据我们所知，我们的算法是少数适用于复杂3D环境的3D Lidar-IMUu融合算法之一。

III. 符号与说明

A. 符号

我们将3D激光雷达传感器捕获的每一条激光测量线表示为scan CCC，而sweep SSS表示包含一次测量中所有scan。例如，16行3D激光雷达在一次扫描中包含16次scan。

在以下部分中，我们将变换矩阵表示为Tba∈SE(3)T_b^a \in SE(3)Tba∈SE(3) ，表示将Fb\mathcal{F_b}Fb坐标系下的点 xb∈Rx^b \in \mathbb{R}xb∈R转换到Fa\mathcal{F_a}Fa坐标系下。Tbaˉ\bar{T_b^a}Tbaˉ表示IMU预测的转换矩阵，Rba∈SO(3)R_b^a \in SO(3)Rba∈SO(3)和pbap_b^apba是TbaT_b^aTba的旋转矩阵和平移向量。同时与RbaR_b^aRba对应的是用四元数qbaq_b^aqba表示旋转，⊗\otimes⊗表示四元数乘法。我们使用ak^\hat{a_k}ak^和wk^\hat{w_k}wk^表示IMU在kkk时刻的原始测量数据。Fa\boldsymbol{F_a}Fa表示在初始坐标系Fa\mathcal{F_a}Fa下提取的特征点，可以通过通过转换矩阵转换到Fb\mathcal{F_b}Fb坐标系下，表示为Fab\boldsymbol{F_a^b}Fab。

B. IMU运动

1) 状态

FBi\mathcal{F_{B_{i}}}FBi和FLi\mathcal{F_{L_{i}}}FLi分别表示IMU和Lidar中心的参考坐标系，与此同时，我们获得每一个sweep SiS_iSi的时间戳iii。我们需要估计的IMU状态XBiWX_{B_i}^WXBiW是在世界坐标系FW\mathcal{F_W}FW下的，Lidar和IMU之间的外参可以表述为TBLT_B^LTBL,IMU在iii时刻的状态和外参可以表述为：

其中pBiWp^W_{B_i}pBiW、vBiWv^W_{B_i}vBiW和qBiWq_{B_i}^WqBiW分别表示位置、速度以及IMU坐标系到世界坐标系的转换。而bab_aba和bgb_gbg分别表示加速度计和陀螺仪的零偏。

2) 运动模型

通过IMU中加速度计和陀螺仪的输入，我们可以通过公式2所示的离散化方程将IMU的之前的状态XBiWX_{B_i}^WXBiW更新到当前的状态XBjWX_{B_j}^WXBjW,其中Δt\Delta tΔt表示两次连续IMU测量的时间间隔。同时我们将激光雷达第k=ik=ik=i帧和第k=jk=jk=j帧中间的所有IMU测量值进行积分。我们不严谨的使用k=j−1k=j-1k=j−1表示在第k=jk=jk=j帧之前的IMU测量值。

其中gWg^WgW是世界坐标系下的重力向量，为了使用方便，我们速记(⋅)i≐(⋅)BiW，Δtij=∑k=ij−1(\cdot )_i \doteq (\cdot)^W_{B_i}，\Delta t_{ij} = \sum^{j-1}_{k=i}(⋅)i≐(⋅)BiW，Δtij=∑k=ij−1，以及将∏j−1k=i\prod{j-1}_{k=i}∏j−1k=i表示为四元数的连续相乘。

3) 预积分

在iii时刻和jjj时刻之间的运动，我们可以通过预积分来表示：zji={Δpij,Δvij,Δqij}z^i_j=\{\Delta p_{ij} , \Delta v_{ij} , \Delta q_{ij}\}zji={Δpij,Δvij,Δqij},其中协方差矩阵CBjBiC_{B_j}^{B_i}CBjBi表示状态误差模型。

IV. 紧耦合Lidar-IMU里程计

为了进行有效的运动估计，许多雷达建图的工作如[7]、[1]和[13]将任务分为两个部分：里程计和建图。受这些工作的启发，我们提出的系统由两个并行的部分组成。在Sec. IV介绍了第一部分，一个紧耦合的Lidar-IMU里程计，优化滑动窗口内的所有状态。在Sec. V介绍了第二部分：细化旋转约束（建立全局一致地图的过程），使用已经优化的位姿和重力约束中的信息将激光雷达sweep与全局地图对齐。

A. Lidar-IMU里程计概述

图1：激光雷达-IMU里程计框架。讨论每一块的小节在循环号之后给出

图2：滑动窗口。局部地图由第jjj帧之前且从ooo开始历史点云组成。优化窗口包含ppp和jjj之间的Lidar点云和距离变化，与[1]中一样；即选择面特征和线特征点

图1展示了我们提出的Lidar-IMU里程计的一个简单的概述，利用先前估计的状态，我们可以使用当前雷达的原始输入SjS_jSj已经历史时刻iii到当前时刻jjj的IMU原始输入IijI_{ij}Iij以进行一次新的状态优化。里程计估计步骤如下：

在SjS_jSj到达之前，IMU状态通过公式(2)迭代进行更新
同时，这些输入的Δpij\Delta p_{ij}Δpij、Δvij\Delta v_{ij}Δvij和Δqij\Delta q_{ij}Δqij经过预积分被用于联合优化。
当接收到一帧激光数据SjS_jSj时,对原始数据进行去畸变处理，以获得去畸变激光雷达帧Sjˉ\bar{S_j}Sjˉ（Sec. IV -B）。
接下来，进行特征提取减少数据维度和提取更重要的特征点FLi\boldsymbol{F_{L_i}}FLi（Sec. IV -B）。
根据之前相应的优化状态TBo,iT_{B_{o,i}}TBo,i和TBLT_B^LTBL（Sec. IV -c)，将滑动窗口下的历史雷达特征点云FLo,i\boldsymbol{F_{L_{o,i}}}FLo,i被合并成一个局部子图MLo,iLpM_{L_{o,i}}^{L_p}MLo,iLp。
通过预测的Lidar位姿Fj\boldsymbol{F_j}Fj，可以找到相对测量关系mLp+1,jm_{L_{p+1,j}}mLp+1,j(Sec. IV -C)。
最后一部是联合非线性优化，根据雷达测量的相对关系和IMU预积分可以获得滑动窗口内的状态的最大后验估计。优化结果用于更新步骤1）中的预测状态，避免IMU传播的漂移。

B. 去畸变和特征提取

3D激光雷达内部有旋转机构，可接收整圈数据。当3D激光雷达移动时，其原始数据SjS_jSjj会发生运动失真，这使得扫描中的点与真实位置不同。为了解决这个问题，我们使用IMU传播预测激光雷达运动TˉLj′Lj\bar{T}_{L_j'}^{L_j}TˉLj′Lj，并假设一个完整的sweep中运动符合线性模型。然后，每个点x(t)∈Sj⊂Rx(t) \in S_j \subset \mathbb{R}x(t)∈Sj⊂R通过TˉLj′Lj\bar{T}^{L_j}_{L_{j'}}TˉLj′Lj进行线性插值校正获得校正后的点云帧Sjˉ\bar{S_j}Sjˉ，其中t∈(tj′,tj]t \in (t_{j'} , t_j]t∈(tj′,tj]是扫描点的时间戳，tj′t_{j'}tj′和tjt_jtj分别是扫描开始和结束的时间戳。

为了提高计算效率，需要进行激光雷达特征提取。这里，我们只对平面或边缘上最相似的点感兴趣[1]，[7]，因为这些点可以从激光雷达扫描sweep的每根线束scan中提取。这些在点云Sjˉ\bar{S_j}Sjˉ中特征点FLj\boldsymbol{F_{L_j}}FLj通过曲率进行选择。

C. 相对Lidar测量

通过IMU和另一个传感器的融合，该传感器能够提供传感器对的相对姿态，估计当前状态。如果固定第一个第一个参考帧，XBWX^W_BXBW和TBLT^L_BTBL也是局部可观的。为了适当的融合来自IMU的预积分值，我们提出了使用两帧之间的雷达相对测量来约束雷达位姿，如算法1所示。在找到对应点之前，我们构建一个局部地图，因为单次扫描中的点不够密集，无法计算准确的对应。

为了简化，我们将预测变换矩阵TˉLjLp\bar{T}_{L_j}^{L_p}TˉLjLp表示为TLjLpT_{L_j}^{L_p}TLjLp，FLjLp\boldsymbol{F_{L_j}^{L_p}}FLjLp通过TˉLjLp\bar{T}_{L_j}^{L_p}TˉLjLp转换而来

局部地图包含来自离散时间戳NmN_mNm {o,....,p,.....,i}\{o,....,p,.....,i \}{o,....,p,.....,i}的激光雷达特征点。其中ooo，ppp和iii分别是窗口内第一帧、中心位置和最后一帧雷达点云的时间戳，如图2所示。局部地图MLo,iLpM^{L_p}_{L^{o,i}}MLo,iLp是建立在pivot lidar帧的坐标系下，通过已经优化好的TLγLpT^{L_p}_{L_\gamma}TLγLp将历史帧的点云转换成FLγLp\boldsymbol{F_{L_\gamma}^{L_p}}FLγLp，γ∈{o,......,i}\gamma \in \{ o , ...... ,i\}γ∈{o,......,i}。待估计的状态是时间戳NsN_sNs中的 {p+1,......,i,j}\{p+1,......,i,j\}{p+1,......,i,j}处的状态，其中p+1p+1p+1和jjj分别是靠近主节点pivot lidar帧时间戳以及窗口中的当前激光雷达点云帧时间戳。

通过建立的局部地图，我们可以找到MLo,iLpM^{L_p}_{L^{o,i}}MLo,iLp和FLα\boldsymbol{F}_{L_\alpha}FLα，α∈{p+1,......,j}\alpha \in \{p+1,......,j\}α∈{p+1,......,j}。我们将这种对应定义为相对激光雷达测量，因为它们是相对于pivot 位姿的，并且pivot 位姿会随着滑动窗口的改变而改变。Sec. IV -B 部分中，我们在FLα\mathcal{F}_{L_{\alpha}}FLα坐标系下提取的原始特征是平面点或者边缘点。在实际应用中，我们发现边缘点不能改善Lidar-IMU里程计的结果。因此，在下文中，我们只讨论平面特征。对变换后的每个特征点xLp∈FLαLp\boldsymbol{x}^{L_p} \in \boldsymbol{F}_{L_\alpha}^{L_p}xLp∈FLαLp使用KNN寻找MLo,iLpM^{L_p}_{L^{o,i}}MLo,iLp中的KKK个近似点π(xLp)\pi ({x}^{L_p})π(xLp)，然后对平面点进行拟合，将这些邻接点拟合到FLp\mathcal{F}_{L_{p}}FLp的一个平面上。平面点的系数可以用线性方程ωTx′+d=0,x′∈π(xLp)\omega^{T} x'+d=0, x' \in \pi({x}^{L_p})ωTx′+d=0,x′∈π(xLp)来求解，其中ω\omegaω是平面法线方向，ddd是到FLp\mathcal{F}_{L_{p}}FLp原点的距离。我们表示每个平面特征点x∈FLαx \in \mathcal{F}_{L_{\alpha}}x∈FLα的m=[x，ω，d]∈mLαm=[x，ω，d] \in m_{L_\alpha}m=[x，ω，d]∈mLα作为相对激光雷达测量之一。值得一提的是，在每个相对激光雷达测量m∈mLαm \in m_{L_\alpha}m∈mLα中，xxx在FLα\mathcal{F}_{L_{\alpha}}FLα中定义，并且ω\omegaω和ddd在FLp\mathcal{F}_{L_{p}}FLp中定义。

D.激光雷达扫描匹配

Lidar的相对测量关系可以提供pivot雷达帧位姿和接下来的雷达帧位姿的相对约束。我们的方法优化了滑动窗口内所有位姿，包括第一个位姿TLpWT^W_{L_p}TLpW,即FLp\mathcal{F}_{L_{p}}FLp未固定。因此，Lidar的损失函数中每一项都包括两帧Lidar点云的位姿，TLpWT^W_{L_{p}}TLpW和TLαWT^W_{L_{\alpha}}TLαW，α∈{p+1,.....j}\alpha \in \{p+1,.....j\}α∈{p+1,.....j}。优化pivot位姿将有助于更好地减小预积分误差，并确保传感器对与重力对准。我们估计的状态是IMU的状态，因此需要引入外部参数来表示IMU状态对激光雷达的约束。可以定义窗口中从后面激光雷达姿态到pivot姿态的相对变换：

利用前面的对应关系，每个相对激光雷达测量值m=[x，ω，d]∈mLα,α∈{p+1,......,j}m=[x，ω，d] \in m_{L_\alpha} , \alpha \in \{p+1 , ...... ,j \}m=[x，ω，d]∈mLα,α∈{p+1,......,j}的残差可以表示为点到面的距离:

E. 优化

为了获得优化的状态，采用了固定滞后平滑器和边缘化方法。固定滞后平滑器将NsN_sNs个IMU状态保持在滑动窗口中，从XBpWX_{B_p}^WXBpW到XBjWX_{B_j}^WXBjW，如图2所示。滑动窗口有助于限制计算量。当新的测量约束到来时，更平滑的将包括新状态并将窗口中最旧的状态边缘化。将优化的的所有状态为：

然后，最小化具有马氏范数的下面的代价函数，以获得状态XXX的MAP估计:

其中，rp(X)\boldsymbol{r}_{\mathcal p} (X)rp(X)是边缘化后的先验约束，rL(m,X)\boldsymbol{r}_{\mathcal L} (m,X)rL(m,X)是相对雷达约束的残差,rB(zβ+1β)\boldsymbol{r}_{\mathcal B} (z^{\beta}_{\beta +1})rB(zβ+1β)是IMU约束的残差，非线性最小二乘形式的代价函数可以用HθX=−bH\theta X=−bHθX=−b的高斯-牛顿算法来求解。我们使用Ceres Solver[15]来求解该问题。

激光雷达约束rL(m,X)\boldsymbol{r}_{\mathcal L} (m,X)rL(m,X)可以从公式(4)每次相对激光雷达测量导出。协方差矩阵CLαmC^m_{L_\alpha}CLαm由激光雷达精度决定。与文献[10]类似，IMU约束rB(zβ+1β)\boldsymbol{r}_{\mathcal B} (z^{\beta}_{\beta +1})rB(zβ+1β)可以由状态和IMU预积分得到，∣∣rp(X)∣∣2=bpTHp+bp||\boldsymbol{r}_{\mathcal p} (X)||^2 =b^T_{\mathcal p} H^+_\mathcal p b_\mathcal p∣∣rp(X)∣∣2=bpTHp+bp可以由舒尔补得到(详见[12])。

V.使用旋转约束进行细化

将特征点注册到全局地图(而不是局部地图)可以将激光雷达位姿约束到一致的世界坐标系FW\mathcal F_WFW。我们的细化方法使用了Sec. IV中的相对雷达测量mLm_LmL。因为全局地图是细化方法的附带产品，我们也称之为建图方法。将最新的激光雷达特征点与全局地图对齐的损失函数可定义为：

其中T=TLWT=T^W_LT=TLW是最新的待估计激光雷达姿态，而Lidar相对测量值mmm中的特征点xxx定义在FL\mathcal F_LFL中。系数ω\omegaω，ddd定义在FW\mathcal F_WFW。然后，我们可以使用类似的高斯-牛顿法来最小化CMC_\mathcal MCM。通过残差CMC_\mathcal MCM和Jacobian 矩阵JpCJ^C_pJpC和JθCJ^C_\thetaJθC进行优化，其中θ\thetaθ是相应四元数qqq的误差状态，但是由于累积的旋转误差和长期的操作，合并后的全局图不能与重力精确对准。这可能会导致进一步建图错误地与倾斜的地图对齐。根据[16]中使用SE(2)约束优化SE(3)的建议，我们提出了一种约束建图策略。该策略利用激光雷达IMU里程计的旋转约束，确保最终地图始终与重力对齐。图3示出了旋转约束建图的结构。

图3：旋转约束建图。里程计姿势首先用作全局点云注册的先验。然后将里程计中的旋转分量作为虚拟测量值应用于约束优化。

鉴于沿z轴的方向具有更高的不确定性，且方向的其他两个自由度更接近真实值，我们可以通过修改旋转的雅可比矩阵来约束损失函数，如[12]中的详细推导所示，
其中(⋅)˘\breve{(\cdot)}(⋅)˘表示上一次迭代的状态，Ω˘z\breve{\Omega }_zΩ˘z姿态相对于坐标系FWF_WFW的旋转信息矩阵的近似值，ϵx\epsilon _xϵx和ϵy\epsilon _yϵy可以通过x、y轴与z轴的信息比得到。

然后，我们使用JpCJ^C_pJpC和JθzCJ^C_{\theta_z}JθzC作为雅可比矩阵，这是优化步骤所需要的。激光雷达姿态的增量可表示为δθz\delta \theta_zδθz和δp\delta pδp，这可以实现更新激光雷达状态p~\tilde{p}p~和q~\tilde{q}q~。

VI. 实验

未完待续~

VII. 测试与分析