0. 简介

占用地图是机器人系统中推理环境未知和已知区域的基本组成部分。《Occupancy Grid Mapping without Ray-Casting for High-resolution LiDAR Sensors》介绍了一种高分辨率LiDAR传感器的高效占用地图框架，称为D-Map。该框架引入了三个主要创新来解决占用地图的计算效率挑战。首先，我们使用深度图像来确定区域的占用状态，而不是传统的射线投射方法。其次，我们在基于树的地图结构上引入了一种高效的树上更新策略。这两种技术避免了对小单元的冗余访问，显著减少了需要更新的单元数量。第三，我们通过利用LiDAR传感器的低误报率，在每次更新时从地图中删除已知单元。这种方法不仅通过减小地图大小提高了框架的更新效率，还赋予了它一种有趣的递减特性，我们称之为D-Map。为了支持我们的设计，我们提供了深度图像投影的准确性和占用更新的时间复杂度的理论分析。此外，我们在公共和私有数据集上进行了广泛的基准实验，使用了各种LiDAR传感器。我们的框架在与其他最先进方法相比的同时，展示了卓越的效率，同时保持了可比较的映射精度和高内存效率。我们展示了D-Map在手持设备和携带高分辨率LiDAR的航空平台上进行实时占用地图的两个真实世界应用。此外，我们在GitHub上开源了D-Map的实现

图1. 我们提出的框架D-Map作为自主无人机探索古堡任务中的实时高分辨率占据地图模块。(a) 无人机采集的高保真度点云。(b) 场景的鸟瞰图。(c) 航空平台携带了一台128通道激光雷达(OS1-128)来执行探索任务。本文的附带视频可在Youtube上观看：youtu.be/m5QQPbkYYnA。

1. 主要贡献

在这项工作中，我们提出了一种新颖的映射框架，有效解决了基于激光雷达的占据地图中的上述问题。我们的贡献可以总结如下：

• 我们提出了一种基于深度图像投影的占据状态确定方法，以减轻传统射线投射技术中的计算负载。这种基于投影的方法可以对任意大小的单元进行占据状态的确定，从而在大规模环境中实现后续的高效更新。
• 我们提出了一种基于混合地图结构的新颖的树上更新策略，以提供计算和内存效率之间的卓越平衡。混合地图结构将未知空间存储在八叉树上，从而实现对大规模未知空间的内存高效表示，而占据空间则存储在哈希网格地图上。在效率方面，所提出的策略允许在八叉树上确定大单元的占据状态，从而避免对小单元进行不必要的更新，提高了效率。
• 我们利用激光雷达测量的低误报率，直接在每次更新时删除具有确定状态（即占据或空闲）的单元。这种方法使我们的地图结构具有递减性质，因此我们将我们的框架称为D-Map，提供更高的计算效率和更少的内存使用。
• 我们对所提出的占据状态确定方法的准确性进行了深入分析，并对D-Map中的更新和查询的时间复杂度进行了分析。具体而言，我们推导出了一个分析函数，用于量化与深度图像分辨率相关的准确性损失。对D-Map上的更新的时间复杂度分析为我们相对于依赖射线投射的最先进方法的卓越性能提供了理论支持。
• 我们在Github上提供了D-Map的实现：github.com/hku-mars/D-Map，以促进我们工作的可复现性和进一步的发展。

2. 综述

图2展示了提出的D-Map框架的概述。橙色块中勾勒的地图结构在第V-A节中进行了解释。它由两部分组成：占用地图和未知地图。绿色块表示占用更新策略的流程。在每次更新时，从传感器姿态处的点云生成深度图像（详见第IV-A节）。随后，在深度图像上构建2D段树，以实现高效的占用状态确定（详见第IV-B节和第IV-C节）。更新未知地图的整个过程在第V-B节中进行了描述，并总结如下。单元提取模块从八叉树上按照从最大到最小的尺寸收回未知单元，将它们投影到深度图像上，并确定它们的占用状态。确定为已知的单元直接从地图中移除，而未知的单元保留下来，而未确定的单元则被分割成更小的单元以进一步确定占用状态。这种由粗到细的过程可以直接对大单元进行更新，而无需查询每个小单元。此外，已知单元的移除赋予了我们的框架递减属性，从而在计算和内存方面提供了高效性。

图2. D-Map的框架概述。蓝色块显示了D-Map的输入，包括点云和相应的传感器里程计。
橙色块是D-Map的占据地图结构，由一个哈希网格地图和一个八叉树组成，用于维护占据空间和未知空间。占据更新策略显示在绿色块中，该策略从八叉树中提取感知区域内的单元格，并根据使用深度图像的占据状态确定方法进行操作。

3. 深度图像上的占用状态确定

本节描述了如何从传入的点云栅格化深度图像，并确定其占用状态。

3.1 深度图像栅格化

为了确定占用状态，由激光雷达传感器捕获的点云被栅格化成当前传感器姿态下的深度图像。为确保状态确定的准确性，深度图像的分辨率应足够小，以便从地图到深度图像的单元格的投影区域大于一个像素。如图3所示，我们使用以下方程确定深度图像分辨率ψ_{map}相对于地图分辨率d和激光雷达的检测范围R：

然而，高分辨率地图将导致一个巨大尺寸的高分辨率深度图像，其中许多像素为空，这是由于点云数量远远小于深度图像的尺寸所致。为了解决这个问题，我们通过LiDAR角分辨率来限制深度图像的分辨率，这是激光雷达发射和接收的两个激光脉冲之间的最小角度。具体来说，我们定义标准深度图像分辨率ψ_I为

其中 ψ_{lidar} 是激光雷达的角分辨率。请注意，为了定义简单起见，我们不区分垂直和水平角分辨率。此外，当相邻点投影到同一像素时，我们保留最小深度值。

3.2 2-D Segment Tree

为了确定地图中单元格的占用状态，采用了一个两步过程，首先将单元格投影到深度图像上，然后比较投影深度与相应区域的最小和最大深度值。由于单元格在深度图像上的投影区域随单元格位置而变化，因此采用了2-D段树结构来加快对深度图像上最小和最大值的高效查询，具体如下所述。段树是一种完全平衡的二叉树，通过表示一组区间来高效提供范围查询[53]。图4描述了通过1-D段树查询最小值的过程。段树通过递归地将数组一分为二来构建，直到每个节点包含一个单独的元素。由于段树上的每个节点表示数组的一个区间，因此在构建过程中预处理区间中的值的总体信息，如最小值和最大值，以加速后续查询。在查询时，段树使用一部分节点（图4中的着色节点）检索查询区间的最小表示。结果是通过总结从检索节点中获取的信息，比直接查询更少的操作得到的。在1-D段树上的查询时间复杂度为O(log N)，其中N是离散数组中的元素数量，而直接查询的时间复杂度为O(N)。将1-D段树扩展为2-D结构的方法涉及构建“段树的段树”，如[54]中提出的。外层的段树通过行来分割2-D数组，对外部段树的每个节点，构建一个1-D内部段树来维护覆盖行的列信息。2-D段树上的查询首先搜索外部树，找到代表被查询区域覆盖的行的节点，然后遍历相应节点上的内部段树，以检索被覆盖的列。最后，从检索节点上存储的信息中总结出被查询区域的结果。在2-D数组大小为N×M的情况下，2-D段树上的查询时间复杂度为O(log N log M)，而直接查询的时间复杂度为O(N^2)。给定从点云光栅化的深度图像，我们构建了一个2-D段树来维护每个树节点上的最小和最大深度值，分别表示为dMin和dMax。此外，我们还跟踪了在每个节点覆盖区域内由点云占据的像素数量，表示为dSum。

图3. 该图说明了地图分辨率d、检测范围R和深度图分辨率ψ_{map}之间的空间关系。

图4. 该图示例了在一维段树上快速查询像素范围[2, 7]内最小值的示例。从段树的根节点开始，范围查询沿着树递归搜索，直到当前节点范围完全被查询范围覆盖为止，此时将返回在树构建过程中保存在节点上的节点范围的最小值。在这个示例中，范围[2, 7]导致四个节点分别表示范围[2, 2]，[3, 4]，[5, 6]和[7, 7]。范围的最小值可以高效地从这四个节点中获取，而不是计算数组中的六个元素。

3.3 占用状态确定(重点内容)

我们以五个单元格为例介绍了我们确定单元格占用状态的方法原理，如图5所示，编号从1到5。我们首先根据单元格是否完全位于激光雷达感知区域内对其进行分类。如图5(b)中的俯视图所示，Grid 1、Grid 2和Grid 3完全位于感知区域内，而Grid 4和Grid 5只有部分位于其中。在完全位于感知区域内的单元格中，Grid 3由于位于观察环境中所有物体的前方而被确定为已知；Grid 1由于位于物体后方而被确定为未知。Grid 2的占用状态仍未确定，因为其一部分位于物体前方，另一部分位于物体后方。对于部分位于感知区域内的单元格，Grid 5由于位于物体后方而被确定为未知。虽然位于物体前方，但Grid 4的占用状态由于其位置并非完全位于感知区域内而仍未确定。

根据上述原则，D-Map将完全位于或与激光雷达感知区域相交的单元格投影到最近点云所栅格化的当前深度图像上，该点云代表环境中的物体。通过比较它们的深度值来确定单元格与物体之间的相对位置。深度图上的占用状态确定过程如算法1所述，并在下文中进行解释。
考虑到一致性和计算效率，我们通过单元格的内切圆将其投影到深度图像上，如图5(b)和(c)所示。首先将单元格中心C表示为球坐标(r, θ, ϕ)（第2行）。描述单元格覆盖的深度范围为其相对于中心C和单元格大小L的最小深度BoxMin和最大深度BoxMax（第3行）。深度图上投影区域的边界（即图5(c)中的查询像素）如下获得。

然后，我们准备以下信息来确定单元的占用状态。首先，计算投影区域覆盖的查询像素数量，即在ProjPixels（第5行）中计数。其次，通过查询二维段树（如第IV-B节所述）提供投影区域内占用像素dSum、最小深度dMin和最大深度dMax的数量（第6行）。第三，我们定义单元的观测完整度α，以确定当前深度图像是否完全观测到单元，该值通过将占用像素dSum除以深度图像上的查询像素ProjPixels计算得出（第7行）。计算得到的α量化了单元的投影区域是否已被当前深度图像中的点充分观测到。较大的α表明单元的大部分位于感知区域内，其投影区域中的大多数像素在深度图像上被LiDAR点积极占用。只有具有较大α值的单元才会更新其占用状态。

通过比较深度图像的深度范围（由最小深度dMin和最大深度dMax表示）与单元的深度范围（由最小深度BoxMin和最大深度BoxMax表示，如图6所示），获取单元的占用状态。首先，通过比较观测完整度α与完整度阈值ε，对单元进行分类。对于完全观测到的单元（即第9行中的α > ε），如果深度图像的dMax小于单元的BoxMin（第10∼12行），则确定占用状态为未知。这种情况表明单元位于环境中的物体后方，如图6(a)所示。当深度图像的dMin大于单元的BoxMax时，确定单元的占用状态为已知（第13∼15行），表明它位于物体前方，如图6(b)所示。否则，无法确定单元的占用状态，如图6(c)所示。对于未完全观测到的单元（即α ⩽ ε），如果深度图像的dMax小于单元的BoxMin（第20∼22行），则确定占用状态为未知，表明单元位于物体后方，如图6(d)所示。如果不满足此条件，则考虑由于观测不完整而无法确定单元的占用状态（第22∼24行），如图6(e)所示。

在某些情况下，单元的投影区域可能小于深度图像上的一个像素区域，如图7所示。这种情况发生在深度图像分辨率ψ_I由LiDAR的角分辨率计算得出时，被称为“单像素情况”。在这种情况下，回想一下，像素仅保存一个具有最小深度的LiDAR点，如第IV-A节所述，我们应进一步确定LiDAR点是否实际穿过了单元，即检查点是否位于单元的投影区域内。如果LiDAR点位于投影区域外，则单元的占用状态仍然未知。否则，我们使用与算法1类似的逻辑来确定占用状态。

3.4 深度图像分辨率分析

深度图像的分辨率是一个关键参数，对于平衡占用状态确定的准确性和效率起着至关重要的作用。低分辨率的深度图像更容易查询，但会导致原始点云信息的丢失增加。相反，高分辨率的深度图像提供更好的准确性，但需要更大的计算量，从而降低了效率。当考虑具有标准分辨率ψ_I的深度图像M_I时，可以通过放宽因子γ将分辨率放宽到ψ = γψ_I，以实现更高的效率，但牺牲一定的准确性。为了评估分辨率放宽因子γ导致的准确性损失，我们推导了一个描述保留准确性的函数f(γ)，作为因子γ的函数来描述从原始深度图像M_I中保留的准确性。
我们假设使用激光雷达扫描一个静止的大型开放环境，没有障碍物。当通过投影到具有标准分辨率ψ_I的深度图像MI来确定自由空间时，结果的占用地图应该在一个3D球形域内显示完全的自由空间，其体积表示为V_I。然而，当将深度图像分辨率放宽到ψ = γψ_I以获得放宽的深度图像M时，由于单像素情况，一些自由空间会被错误地确定为未知。我们将这种情况的体积表示为V。V_I和V的示意图在补充材料中的图1中呈现[55]。准确性函数f(γ)被定义为正确确定为自由空间的比率：

通过一系列模拟实验验证了准确性函数 f(γ)。这些实验涉及各种LiDAR检测范围和地图分辨率。具体而言，给定检测范围 R 和地图分辨率 d，会在半径为 R 的球面上生成一批覆盖检测范围内自由空间的点云。为了计算体积 V_I，我们计算分辨率为 d 的地图中被投影到原始深度图像 M_I上并乘以单位体积 d^3 的结果，以计算被识别为自由空间的单元格数量。类似地，通过将地图投影到放松的深度图像 M 上，我们得到了体积 V。准确性函数值 f(γ) 使用式(4)中的定义进行计算。在模拟实验中，我们将 LiDAR 的视场角 Φ 设置为 [−15^◦, 15^◦]。

模拟实验的结果见图8。理论模型很好地拟合了实验曲线。然而，当 γ 接近 γ_0 时，由于离散化计算覆盖区域的单元格数量，会出现轻微误差。这些误差在低分辨率和短检测范围的曲线中显而易见（例如，在 R = 10m 和 d = 0.5m 的设置中，即图8中的绿色实线）。