# Tightly-Coupled LiDAR-Visual-Inertial SLAM and Large-Scale Volumetric Occupancy Mapping

Simon Boche$^{1}$, Sebastián Barbas Laina$^{1}$, Stefan Leutenegger$^{1,2,3}$

Abstract— Autonomous navigation is one of the key requirements for every potential application of mobile robots in the real-world. Besides high-accuracy state estimation, a suitable and globally consistent representation of the 3D environment is indispensable. We present a fully tightly-coupled LiDAR-Visual-Inertial SLAM system and 3D mapping framework applying local submapping strategies to achieve scalability to large-scale environments. A novel and correspondence-free, inherently probabilistic, formulation of LiDAR residuals is introduced, expressed only in terms of the occupancy fields and its respective gradients. These residuals can be added to a factor graph optimisation problem, either as frame-to-map factors for the live estimates or as map-to-map factors aligning the submaps with respect to one another. Experimental validation demonstrates that the approach achieves state-of-the-art pose accuracy and furthermore produces globally consistent volumetric occupancy submaps which can be directly used in downstream tasks such as navigation or exploration.

### I. INTRODUCTION

Robust and accurate state estimation is one of the fundamental components for autonomous navigation of robotic systems. But positioning is only part of the problem. An accurate and especially globally consistent representation of the 3D environment that the robot is operating in, is also indispensable. Simultaneous Localisation and Mapping (SLAM) approaches fusing multiple sensor sources, such as stereo vision, Inertial Measurement Units (IMU) or Light Detection and Ranging (LiDAR) sensors, have proven to achieve accurate performance in localisation. Robustness is gained by fusing complementary sensors to overcome degrading scenarios for each of the individual sensors. Most state-of-the-art LiDAR-Visual-Inertial (LVI) SLAM systems are representing the 3D world in terms of features or point clouds, e.g. [^1], [^2], [^3], [^4]. Lately, also surfels have been adopted to LiDAR-based SLAM systems [^5]. While these representations may prove suitable for state estimation and surface reconstruction, they cannot be directly used for downstream tasks such as robotic path planning, where an explicit representation of observed free space is desirable.

That is why, lately, there has been research in coupling SLAM and volumetric mapping. To account for potential pose drift, which would lead to degrading map accuracy and

![Image](images/3655bd0e47f3f1a18658b85d9e6449b3ea8f73526a400772d53bb5e0779694a5.jpg)

<div style="text-align: center;">Fig. 1: Horizontal slices of the occupancy fields and 3D reconstructions from Sequence Exp21 of the HILTI SLAM Challenge [^8]. All submaps are overlayed. Meshes of the surfaces (brown) are extracted from the occupancy field as zero-crossings. The blue area denotes free space extracted for gravity-aligned slices through the submaps. The estimated trajectory is shown in green.</div>

consistency, submapping approaches have recently gained interest. The idea behind that is to divide the environment into several local submaps, which will stay accurate and consistent due to locally limited drift. Different strategies on how to create these submaps and to keep global consistency across submaps have been investigated, e.g. in [^6] or [^7].

Most existing methods decouple the problems of SLAM and volumetric mapping and rearrange the relative position and orientation of all submaps with respect to each other in a loosely-coupled optimisation problem. In contrast, we aim to formulate a whole unified tightly-coupled problem to keep our 3D map representations consistent at all times which makes them suitable for usage in robotic applications. In short, the key contributions of this work are the following.

- We present a tightly-coupled optimisation-based LVI-SLAM and occupancy mapping framework leveraging local submaps to ensure global consistency.

- We introduce a novel residual formulation for LiDAR-based error terms that directly operates on the occupancy values and field gradients. Both can be efficiently queried from the submaps without an expensive data association step as for example in Iterative Closest Point (ICP) approaches. These residuals can be added to the factor graph in two ways, either as factors for every live frame or as factors between submaps.

- We evaluate the approach quantitatively in terms of localisation accuracy and qualitatively regarding map quality. We demonstrate that our system yields globally consistent maps for different LiDAR sensors while achieving state-of-the-art localisation performance.

---

The remainder of this paper is structured as follows: after a brief summary of related works in Sec. II, Sec. III presents the problem statement and notations. Sec. IV sketches the overall structure of the system before the mapping and estimation approaches are described in Sec. V and VI. An experimental evaluation is given in Sec. VII.

### II. RELATED WORK

One of the first breakthroughs and still influential works in the field of LiDAR-based SLAM methods was LOAM [^9]. Its idea was the extraction of geometric features such as edges and planes, and matching them across frames. This basic concept has been widely adopted in many following works. While LOAM was also able to fuse IMU measurements in a loosely-coupled way, LIO-SAM [^10] formulates a tightly-coupled pose graph optimisation fusing edge and plane error residuals with IMU error residuals. Also filter-based approaches, such as FAST-LIO [^11] or FAST-LIO2 [^12] have successfully demonstrated high-accuracy localisation fusing LiDAR and IMU. While FAST-LIO still uses also plane and edge features, FAST-LIO2 achieves significant speed-up by directly operating on the raw points. [^13] establishes sliding-window bundle adjustment (BA) for LiDAR scans.

Recently, LiDAR-Visual-Inertial SLAM systems have become increasingly popular. A large amount of these approaches makes use of LOAM's idea of geometric feature extraction. LVI-SAM [^3] combines two subsystems, a Visual-Inertial (VI) and a LiDAR-Inertial (LI), in a tightly-coupled way to complement each other in challenging scenarios. The LI system builds a factor graph based on IMU preintegration errors and edge and plane residuals. VILENS [^1] also builds an optimisation problem consisting of visual, inertial, leg odometry, and LiDAR-based line and plane residuals. Another line of research uses Multi-State Constraint Kalman Filter (MSCKF [^14]) based approaches for tightly-coupled, filter-based fusion of visual, inertial and LiDAR measurements. Examples are [^15] and [^16].

To minimise drift in long-term scenarios, the latest research commonly applies the concept of submapping. This originates from early SLAM research, such as the Atlas framework [^17]. In this context, additional factors can be derived to align submaps and to eliminate drift. [^1] uses local point-cloud submaps and adds ICP odometry measurements into the factor graph. Wildcat [^5], a sliding-window optimisation-based LiDAR-Inertial Odometry system, achieves peak state-of-the-art robustness and accuracy by building local surfel submaps and aligning the submaps in a pose graph optimisation. Most of the aforementioned work uses feature-based or surface-based 3D representations. While they might be sufficient for high-accuracy localisation and 3D reconstruction, they are not suitable for navigation due to their lack of representing observed free space.

Submapping on volumetric maps has been addressed in various works. [^18] builds occupancy submaps based on OctoMap [^19] and aligns them by standard ICP registration. Voxgraph [^6] uses a Visual-Inertial Odometry to provide poses for integration into TSDF maps. Upon completion, ESDF fields are generated and submaps are aligned in a back-end pose graph optimisation using correspondence-free error terms based on ESDF values. New submaps are created at a fixed frequency. Wang et al. [^20] instead use occupancy maps as their 3D representation. Using Supereight2 [^21] as adaptive-resolution mapping approach, new submaps are spawned based on the distance travelled. Submaps are rearranged based on updates from the visual-inertial estimator. The follow-up work [^7] improved the submap creation by evaluating the point cloud overlaps of new scans and alignment of submaps is based on ICP. In this work, we will adopt the concept of submapping. In contrast to [^6], [^20], [^7], the global alignment of submaps is not decoupled from the estimator but provides direct feedback. We formulate correspondence-free residuals as in [^6] without the expensive need to extract ESDFs as we directly use the available occupancy information.

### III. PRELIMINARIES

### A. Problem Statement

In this work, we aim to build an accurate and globally consistent volumetric representation of the environment around a mobile robot, which can be used in downstream tasks, most prominently navigation. We do so by taking a Visual-Inertial SLAM System, OKVIS2 [^22], as our starting point, as well as an adaptation of the adaptive-resolution occupancy mapping framework Supereight2 [^21]. As VI-SLAM tends to be locally consistent but suffers from the accumulation of larger drift over time, we tackle the problem of degrading map quality by applying submapping strategies. Furthermore, to increase the accuracy of the estimator while simultaneously ensuring global consistency of overlapping maps, we integrate LiDAR constraints based on frame-to-map as well as map-to-map alignment factors into the state estimator in a tightly-coupled way.

### B. Notation

Throughout this work, the following notation will be used: coordinate frames are written as $\underline{\pmb{F}}_{A}$ and a vector expressed in this reference frame will be denoted as $_{A}\mathbf{r}$. The rigid body transformation which transforms points from a reference frame $\underline{\pmb{F}}_{B}$ to another reference frame $\underline{\pmb{F}}_{A}$ is given by $\pmb{T}_{AB} \in SE(3)$ and can be decomposed into a rotation matrix $\mathbf{C}_{AB} \in SO(3)$ and the translational component $_{A}\mathbf{r}_{B}$. We also denote the rotation $\mathbf{C}_{AB}$ with its unit quaternion form $\mathbf{q}_{AB}$. The most important reference frames that will be used are: a fixed world reference frame $\underline{\pmb{F}}_{W}$, the camera coordinate frames $\underline{\pmb{F}}_{C_i}$ for $i = 1 \ldots N$ cameras, the IMU sensor frame $\underline{\pmb{F}}_{S}$, the LiDAR sensor frame $\underline{\pmb{F}}_{L}$ and a map frame $\underline{\pmb{F}}_{M}$. Furthermore, $[\cdot]^{\times}$ represents the skew-symmetric matrix of a 3D vector.

### IV. SYSTEM OVERVIEW

In this work, we build a fully tightly-coupled system combining volumetric mapping and LiDAR-Visual-Inertial (LVI) SLAM. A high-level overview is sketched in Fig. 2. Inputs considered are IMU measurements, stereo camera frames and

---

![Image](images/ceed19be76f27b38a79cca8f5b690f716fd6575b4652103d3be1a452d1c10b38.jpg)

<div style="text-align: center;">Fig. 2: High-level structure of the proposed tightly-coupled LVI-SLAM and submapping framework. Motion compensation for LiDAR point clouds is performed to formulate live frame-to-map factors and to update the current active submap. Upon submap completion, the most overlapping previous submap is determined and map-to-map factors are fused into the optimisation problem.</div>

potentially unordered LiDAR point clouds. By the use of state estimates and IMU measurements integrated to a point's individual timestamp, motion compensation is performed to account for dynamic movement of the LiDAR sensor. After transforming the measurements to the respective map frame, the occupancy field of the current active submap is updated. Tight coupling between the submapping module and the LiDAR-Visual-Inertial estimator is achieved in two ways. Live frame-to-map factors are formulated for the state corresponding to the most recent stereo frame. Furthermore, upon completion of the active submap, the most overlapping submap is detected and map-to-map factors are added to the LVI state estimator. These factors as well as details regarding the submapping strategies will be given in the following sections.

### V. OCCUPANCY MAPPING AND SUBMAPPING

As a volumetric occupancy mapping approach, we adopted the octree-based multi-resolution volumetric mapping framework Supereight2 [^21] specifically for the usage with different kinds of LiDAR sensors.

### A. Occupancy Mapping for dynamical LiDAR sensors

In contrast to [^20], we do not use range image projections of LiDAR point clouds as input to our mapping system. Instead, to account for dynamically moving sensors, we added a ray-based integration interface for Supereight2 that

![Image](images/bbb9744dae2104ee4c363044a9bebe24148a095c676be719b86dbc2eb7a10e6e.jpg)

<div style="text-align: center;">Fig. 3: Inverse sensor model used in Supereight2 [^21] for two measurements (1m, 5m). The log-odds occupancy probability along a ray measurement is expressed as a function of the difference  $d_{r}$ between query points along the ray and the measured distance  $z_{r}$. The occupancy values are clipped in front of the surface at a minimum  $L_{min}$ reached at  $3\sigma$ where  $\sigma$ is a distance dependent uncertainty value. It grows linearly up to half the surface thickness  $\tau(z_{r})$. For more details, see [^21].</div>

fuses measurements on a per-ray-basis. We use fairly standard occupancy mapping in log-odds space with a simplified inverse sensor model that approximates the probabilistic nature of the sensor in terms of range or depth uncertainty and outlier probabilities. The log-odds occupancy $l^k(_{M}\mathbf{p})$ of a 3D point $M\mathbf{p}$ in a map frame $\mathbf{\mathcal{F}}_M$ at step k is given by

 $$ l^{k}\left(_{M}\mathbf{p}\right)=\log\frac{P_{\mathrm{o c c}}\left(_{M}\mathbf{p}|z_{r}\right)}{1-P_{\mathrm{o c c}}\left(_{M}\mathbf{p}|z_{r}\right)}, $$

where $z_{r}$ is the measured distance of a single ray and $l^{k}$ is following a piece-wise linear function along the ray as shown in Fig. 3.

Using clamped weights, we apply the same additive Bayesian updates as in [^21]:

 $$ \begin{aligned}\bar{L}^{k}\left(_{M}\mathbf{p}\right)&=\frac{\bar{L}^{k-1}\left(_{M}\mathbf{p}\right)w^{k-1}+l^{k}\left(_{M}\mathbf{p}\right)}{w^{k-1}+1}\\w_{k}&=\min\left\{w_{k-1}+1,w_{\max}\right\}\end{aligned} $$

The accumulated log-odds can still be preserved as

 $$ \begin{array}{r}{\bar{L}^{k}\left(_{M}\mathbf{p}\right)=\bar{L}^{k}\left(_{M}\mathbf{p}\right)w_{k}.}\end{array} $$

To keep the octree representation consistent and the resolution as coarse as possible, similar up-propagation and tree pruning strategies as in [^21] are applied.

### B. Submapping Strategy

Similar to [^7], our goal is to leverage the available sensor information as a decision criterion when to spawn a new submap: we evaluate per-frame the ratio between incoming measurements that correspond to already observed space and the total number of measurements. If this ratio falls below a certain threshold $\lambda_{\text{overlap}}$, a new submap will be created with the next visual keyframe generated as in [^22]. Every submap is anchored in the IMU frame $\mathcal{F}_S$ at the time step k of the corresponding keyframe. All submaps are stored together with their respective keyframe poses $\mathcal{T}_{WS^k}$.

### VI. LiDAR-VISUAL-INERTIAL ESTIMATOR

In the following, we present the extension of the OKVIS2 [^22] factor graph.

### A. Factor Graph and Optimisation Problem

The factor graph representation used in this work is shown in Fig. 4. On top of visual, inertial and pose graph factors, that are already used in OKVIS2, two different types of LiDAR factors are added to the graph. For every state in the real-time optimisation window, LiDAR factors with respect to the last completed submap are added. Furthermore, map-to-map LiDAR factors are added to constrain submaps with respect to each other by aggregating measurements between submaps. The state vector that is estimated remains:

 $$ \mathbf{x}=\left[_{W}\mathbf{r}_{S}^{T},\mathbf{q}_{W S}{}^{T},{}_{W}\mathbf{v}^{T},\mathbf{b}_{\mathrm{g}}^{T},\mathbf{b}_{\mathrm{a}}^{T}\right]^{T}, $$

where $w_{\mathbf{r}_S}$ and $\mathbf{q}_{WS}$ denote the position and orientation of the IMU sensor frame in the fixed world frame, and $w_V$ describes the velocity of the IMU with respect to the world frame. $\mathbf{b}_g$ and $\mathbf{b}_a$ stand for gyroscope and accelerometer biases, respectively.

---

![Image](images/436c948e9d3f8fa1c87f6f0e4524f2c769b0e7d7b2361ce6b5cd1a3e60b0a740.jpg)

<div style="text-align: center;">● (many) landmarks ⚫ keyframe pose (with observ.) ⚫ non-key pose (with observ.) ● submap keyframe ● posegraph pose (no observ.) ● speed/biases ⚫ many reprojection errors ⚫ IMU error ⚫ posegraph errors ⚫ frame-to-map errors ⚫ map-to-map errors ⚫ pose prior ⚫ speed/bias prior ⚫ estimated variable ⚫ fixed parameters ⚫ active error term ⚫ inactive error term</div>

<div style="text-align: center;">Fig. 4: Optimisation Factor Graph of the LiDAR-Visual-Inertial Estimator. Left: The real-time estimator connects set of current keyframe states and non-keyframe states by IMU errors and visual reprojection errors. It also shows the active submap keyframe and the last completed submap keyframe. For every state in the optimisation window, we can formulate live LiDAR factors between every live frame and the last completed submap. Right: OKVIS2 connects keyframe poses through relative pose errors at a later stage. In addition to the live LiDAR factors, measurements between frames can be aggregated and map-to-map LiDAR factors can be added to the factor graph. Every map will be connected to the previous submap and optionally an older submap if the geometric overlap surpasses a threshold.</div>

### B. Error Residuals

Here, for convenience, a summary of the original OKVIS2 factors is given before deriving the new LiDAR factors.

1) Visual Reprojection Errors: OKVIS2 [^22] uses 2D reprojection errors $\mathbf{e}_{\mathrm{r}}^{i,j,k}$ of the j-th landmark in the frame of the i-th camera at a timestamp k and its corresponding observation $\tilde{\mathbf{z}}_{\mathrm{r}}^{i,j,k}$ in the image is given by

 $$ \mathbf{e}_{\mathrm{r}}^{i,j,k}=\tilde{\mathbf{z}}_{\mathrm{r}}^{i,j,k}-\mathbf{h}\left(\mathbf{T}_{S C_{i}}^{-1}\mathbf{T}_{S^{k}W^{}W}\mathbf{l}^{j}\right). $$

Hereby, h ( $\cdot$) denotes the camera projection.

2) IMU Errors: For the formulation of the IMU residuals, OKVIS2 adopts the IMU preintegration approach in [^23]. Between time steps k and n, the IMU error is:

 $$ \mathbf{e}_{\mathrm{s}}^{k}=\tilde{\mathbf{x}}^{n}\left(\mathbf{x}^{k},\tilde{\mathbf{z}}_{\mathrm{s}}^{k,n}\right) 囗 \mathbf{x}^{n}, $$

where $\tilde{\mathbf{x}}^{n}$ is the predicted state at an arbitrary time n as a function of the current state $\mathbf{x}^{k}$ and IMU measurements $\tilde{\mathbf{z}}_{\mathrm{s}}^{k,n}$. The $\boxminus$ performs regular subtraction except for the quaternion (see [^22]).

3) Relative Pose Errors: Additionally, relative pose errors $\mathbf{e}_{p}^{r,c}$ between time steps r and c are given by

 $$ \mathbf{e}_{\mathbf{p}}^{r,c}=\mathbf{e}_{\mathbf{p},0}^{r,c}+\left[\begin{matrix}{s^{r}\mathbf{r}_{S^{c}}-s^{r}\tilde{\mathbf{r}}_{S^{c}}}\\ {\mathbf{q}_{S^{r}S^{c}}\boxminus\tilde{\mathbf{q}}_{S^{r}S^{c}}}\\ \end{matrix}\right]. $$

with $S^r\tilde{\mathbf{r}}_S^c$ and $({\tilde{\mathbf{q}}_S}^{r_S})^{c}$ being nominal relative position and orientations expressed in the IMU frame $\underline{\mathscr{F}}_S$. We refer the reader to [^22] for a detailed derivation of the constant $\mathbf{e}_{p,0}^{r,c}$ as well as error Jacobians and error weights $\mathbf{W}_{p}^{r,c}$.

4) Submap-based LiDAR Errors: As stated in the previous section, there are two types of LiDAR residuals. In both cases, the error residuals are formulated based on two states. On the one hand, for live residuals, these two states are given by the current frame and the keyframe associated to the last completed submap. On the other hand, for map-to-map residuals, both states will be submap anchor frames. Regarding the formulation of the error residuals, however, they do not differ.

Given a completed submap anchored at keyframe pose $T_{WS^a}$ and a point cloud $\mathcal{P}_b$ associated to another state $T_{WS^b}$, we formulate the residual for every $S^b \mathbf{p} \in \mathcal{P}_b$ as:

 $$ e_{1}^{a,b}\left(_{S^{a}}\mathbf{p}\right)=\frac{d}{\sigma}=\frac{L\left(_{S^{a}}\mathbf{p}\right)}{\sqrt{\frac{L_{\min}^{2}}{9}+\sigma_{z}^{2}\left|\nabla L\left(_{S^{a}}\mathbf{p}\right)\right|^{2}}}, $$

where $S^a \mathbf{p} = T_{S^a S^b S^b} \mathbf{p}$ with $T_{S^a S^b} = T_{WS^a}^{-1} T_{WS^b}$. The idea here is that every measured point should be on a surface in the 3D map; and the distance d of the point from the nearest surface can be extrapolated from the occupancy value $L(\cdot)$ and the occupancy gradient $\nabla L(\cdot)$ assuming a linear behavior as in the sensor model for mapping. From the model (Fig. 3), we can derive the distance d and the map uncertainty as:

 $$ d=\frac{L}{|\nabla L|},\qquad\sigma_{\mathrm{m a p}}=\frac{L_{\mathrm{m i n}}}{3|\nabla L|}. $$

 $L_{min}$ is a configuration parameter and denotes the saturation minimum log-odds occupancy value. With the sensor-specific measurement uncertainty $\sigma_{z}$, we can formulate the weighted residual in Eqn. (8) using the total uncertainty

 $$ \sigma=\sqrt{\sigma_{\mathrm{m a p}}^{2}+\sigma_{z}^{2}}. $$

Furthermore, we want to derive analytic Jacobians, whereby using a perturbation $\delta\chi_T = [\delta\mathbf{r}, \delta\boldsymbol{\alpha}]$ for poses around linearisation points $\overline{\mathbf{r}}$ and $\overline{\mathbf{C}}$:

 $$ \begin{aligned}&\mathbf{r}=\bar{\mathbf{r}}+\delta\mathbf{r},\\&\mathbf{C}=\operatorname{Exp}\left(\delta\boldsymbol{\alpha}\right)\overline{\mathbf{C}}.\\ \end{aligned} $$

We further define the error state as $\delta\chi = [\delta\chi_{T_{WS}}, \delta\chi_{sb}]$, with $\delta\chi_{sb}$ denoting an additive perturbation of the speed and biases. With the LiDAR residuals only depending on pose states, but not on speed and biases, we can compute Jacobians with respect to the poses of frames $\overrightarrow{F}_{S^{a}}$ and $\overrightarrow{F}_{S^{b}}$ leveraging the chain rule:

 $$ \frac{\partial e_{1}^{a,b}}{\delta\mathbf{x}_{T_{W S^{m}}}}=\frac{\partial e_{1}^{a,b}}{\partial_{S^{a}}\mathbf{p}}\frac{\partial_{S^{a}}\mathbf{p}}{\delta\mathbf{x}_{T_{W S^{m}}}}, $$

---

with $m \in \{a, b\}$. The individual Jacobians in Eq. (12) are:

 $$ \begin{align*}\frac{\partial e^{a,b}_{1}}{\partial_{S^{a}}\mathbf{p}}&=\frac{\nabla L}{\sqrt{\frac{L^{2}_{\min}}{9}+|\nabla L|^{2}\sigma^{2}_{z}}}\\\frac{\partial_{S^{a}}\mathbf{p}}{\delta\mathcal{X}_{T_{WS^{a}}}}&=\mathbf{C}_{S^{a}W}\begin{bmatrix}-\mathbf{I}_{3}&\left[\mathbf{C}_{WS^{b}S^{b}}\mathbf{p}+_{W}\mathbf{r}_{S^{b}}-_{W}\mathbf{r}_{S^{a}}\right]^{\times}\end{bmatrix}\\\frac{\partial_{S^{a}}\mathbf{p}}{\delta\mathcal{X}_{T_{WS^{b}}}}&=\mathbf{C}_{S^{a}W}\begin{bmatrix}\mathbf{I}_{3}&-\left[\mathbf{C}_{WS^{b}S^{b}}\mathbf{p}\right]^{\times}\end{bmatrix}.\tag{13}\end{align*} $$

Note that the overall idea of this optimisation resembles a point-to-plane ICP. In contrast to the usual point-to-plane ICP, the expensive step of data association and normal computation can be omitted as occupancy values and gradients can be directly retrieved from the occupancy field. Furthermore, outliers are implicitly covered as they will most likely lie in either free or unobserved space with zero or invalid gradients.

5) Optimisation Problem: All of the aforementioned factors are combined in the overall minimisation objective:

 $$ \begin{aligned}c\left(\mathbf{x}\right)&=\frac{1}{2}\sum_{i}\sum_{k\in\mathcal{K}}\sum_{j\in\mathcal{J}(i,k)}\rho\left(\mathbf{e}_{\mathrm{r}}^{i,j,k^{T}}\mathbf{W}_{\mathrm{r}}\mathbf{e}_{\mathrm{r}}^{i,j,k}\right)\\&+\frac{1}{2}\sum_{k\in\mathcal{P}\cup\mathcal{K}\backslash f}\mathbf{e}_{\mathrm{s}}^{k^{T}}\mathbf{W}_{\mathrm{s}}^{k}\mathbf{e}_{\mathrm{s}}^{k}+\frac{1}{2}\sum_{r\in\mathcal{P}}\sum_{c\in\mathcal{C}(r)}\mathbf{e}_{\mathrm{p}}^{r,c^{T}}\mathbf{W}_{\mathrm{r}}^{r,c}\mathbf{e}_{\mathrm{p}}^{r,c}\\&+\frac{1}{2}\sum_{k\in\mathcal{K}}\sum_{\mathbf{p}\in\mathcal{L}_{k}}e_{1}^{C,k^{2}}+\frac{1}{2}\sum_{b\in\mathcal{M}}\sum_{a\in\mathcal{A}_{b}}\sum_{\mathbf{p}\in\mathcal{L}_{b}}e_{1}^{a,b^{2}}.\quad(14)\end{aligned} $$

Here, the set $\mathcal{K}$ contains the most recent frames as well as keyframes with observations of visible landmarks in $\mathcal{J}(i,k)$. $\mathcal{P}$ contains all pose graph frames and f denotes the most current frame. $\mathcal{C}(r) \subset \mathcal{P}$ is the set of all pose graph frames connected to a frame r. Furthermore, the set $\mathcal{L}_k$ denotes the set of all LiDAR measurements associated to a frame k. $\mathcal{M}$ is the set of all past submaps, and $\mathcal{A}_b$ the set of all submap frames connected to a submap frame $\boldsymbol{T}_{WS^b}$ via map-to-map residuals. C denotes the last completed submap frame.

### VII. EXPERIMENTAL RESULTS

We will quantitatively evaluate the accuracy of the proposed tightly-coupled LiDAR-Visual-Inertial SLAM system in Section VII-B. In section VII-C, we will qualitatively demonstrate that the proposed system yields globally consistent submaps that can be used in robotic applications. Finally, in Section VII-D, we will briefly analyse the computational performance of the proposed system.

### A. Implementation Details

All experiments were performed on an Intel i7-11700K CPU with 3.6 GHz and 64 GB RAM. Google's non-linear least squares optimisation framework Ceres [^25] has been used with an implementation of analytical Jacobians. As a threshold $\lambda_{overlap}$ for new submap creation, a value of 0.4 was chosen. For every live state, we added 100 frame-to-map residuals. Upon completion of submaps, 1000 map-to-map residuals were added. For that, points were randomly sampled from the aggregated point clouds between the two corresponding frames. Submaps have a maximum dimension of 15.36 m and a finest resolution of 3 cm.

### B. HILTI Slam Challenge

We evaluated our approach on the HILTI 2022 SLAM Challenge benchmark [^8]. The Hilti-Oxford dataset was collected using a handheld sensor device consisting of an Alphasense five-camera module, an IMU and a 32 beam Hesai PandarXT-32 LiDAR sensor. With the available calibration datasets, we ran our own sensor calibration of the multicamera system and IMU-camera extrinsics using Kalibr [^26]. As we experienced that the performance on the different sequences dependent significantly on the calibration parameters that were used, we also performed online calibration of the IMU-camera extrinsics. The challenge includes 8 different sequences with varying levels of difficulty. Difficulty levels are based on challenging environments (e.g. bad lighting) or aggressive motions of the sensor suite. Out of the 8 sequences, 3 were taken on a construction site (Exp01, Exp02 and Exp03), 1 was taken in a long corridor of an office building (Exp07) and the remaining 4 sequences were recorded in and around the Oxford Sheldonian Theatre (Exp09, Exp11, Exp15 and Exp21). After the challenge, 3 additional sequences from the construction site (Exp04, Exp05 and Exp06) were added to the automatic evaluation system. For evaluation, sparse ground-truth with millimeter accuracy is provided. An automatic evaluation system is available that uses a score based on the Absolute Trajectory Error (ATE) as an evaluation metric. First, the estimated trajectories are aligned with the sparse ground-truth using SE(3) Umeyama alignment. For every ground-truth control point, the corresponding estimate is retrieved and awarded a score ranging from 0 to 10 (10 for errors below 1 cm, 0 for errors above 10 cm, for more details see [^8]). The total score $S_{j}$ for a dataset j with N control points is computed as

 $$ S_{j}=\left(\frac{1}{10N}\sum_{i=0}^{N}s_{i}\right)\times100. $$

A full score of 100 would correspond to a localisation error below 1 cm across the whole trajectory.

Table I shows the final scores of our approach on all challenge sequences with identical parameters. We report the causal as well as non-causal estimates including a final Bundle Adjustment. We compare it to the Visual-Inertial SLAM (including loop closures) performance of OKVIS2 [^22] as a baseline. The VI baseline used all 5 cameras and the same visual frontend parameters as in the LVI case. Furthermore, we compare it to other state-of-the-art methods that have been published an evaluated on the benchmark including the – at the time of the challenge – winning algorithm Wildcat [^5].

It can be seen that adding LiDAR factors increases the accuracy of OKVIS2 significantly, in the causal as well as the final optimised evaluation. Furthermore, the proposed approach achieves state-of-the-art localisation accuracy. At the moment, a final score of 427.65 achieves rank 7 out of

---

<table border=1 style='margin: auto; word-wrap: break-word;'><tr><td rowspan="2">Approach</td><td colspan="4">Sensors</td><td colspan="8">Sequence</td><td style='text-align: center; word-wrap: break-word;'>Score</td><td style='text-align: center; word-wrap: break-word;'></td></tr><tr><td style='text-align: center; word-wrap: break-word;'>C</td><td style='text-align: center; word-wrap: break-word;'>L</td><td style='text-align: center; word-wrap: break-word;'>V</td><td style='text-align: center; word-wrap: break-word;'>exp01(e)</td><td style='text-align: center; word-wrap: break-word;'>exp02(m)</td><td style='text-align: center; word-wrap: break-word;'>exp03(h)</td><td style='text-align: center; word-wrap: break-word;'>exp04(e)</td><td style='text-align: center; word-wrap: break-word;'>exp05(e)</td><td style='text-align: center; word-wrap: break-word;'>exp06(m)</td><td style='text-align: center; word-wrap: break-word;'>exp07(m)</td><td style='text-align: center; word-wrap: break-word;'>exp09(h)</td><td style='text-align: center; word-wrap: break-word;'>exp11(m)</td><td style='text-align: center; word-wrap: break-word;'>exp15(h)</td><td style='text-align: center; word-wrap: break-word;'>exp21(e)</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>OKVIS2 [^22]</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>8.46</td><td style='text-align: center; word-wrap: break-word;'>8.18</td><td style='text-align: center; word-wrap: break-word;'>0.00</td><td style='text-align: center; word-wrap: break-word;'>17.14</td><td style='text-align: center; word-wrap: break-word;'>23.33</td><td style='text-align: center; word-wrap: break-word;'>14.29</td><td style='text-align: center; word-wrap: break-word;'>0.00</td><td style='text-align: center; word-wrap: break-word;'>0.00</td><td style='text-align: center; word-wrap: break-word;'>0.00</td><td style='text-align: center; word-wrap: break-word;'>0.00</td><td style='text-align: center; word-wrap: break-word;'>16.64</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>OKVIS2 [^22]</td><td style='text-align: center; word-wrap: break-word;'></td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>30.77</td><td style='text-align: center; word-wrap: break-word;'>20.00</td><td style='text-align: center; word-wrap: break-word;'>0.00</td><td style='text-align: center; word-wrap: break-word;'>47.14</td><td style='text-align: center; word-wrap: break-word;'>56.67</td><td style='text-align: center; word-wrap: break-word;'>34.28</td><td style='text-align: center; word-wrap: break-word;'>31.67</td><td style='text-align: center; word-wrap: break-word;'>13.75</td><td style='text-align: center; word-wrap: break-word;'>28.00</td><td style='text-align: center; word-wrap: break-word;'>12.22</td><td style='text-align: center; word-wrap: break-word;'>136.41</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>VILENS [^1]</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>69.23</td><td style='text-align: center; word-wrap: break-word;'>49.09</td><td style='text-align: center; word-wrap: break-word;'>40.00</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>16.67</td><td style='text-align: center; word-wrap: break-word;'>5.00</td><td style='text-align: center; word-wrap: break-word;'>68.00</td><td style='text-align: center; word-wrap: break-word;'>17.78</td><td style='text-align: center; word-wrap: break-word;'>325.77</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>HKU $^{{*}}$</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>84.62</td><td style='text-align: center; word-wrap: break-word;'>70.00</td><td style='text-align: center; word-wrap: break-word;'>44.71</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>31.67</td><td style='text-align: center; word-wrap: break-word;'>7.50</td><td style='text-align: center; word-wrap: break-word;'>92.00</td><td style='text-align: center; word-wrap: break-word;'>24.44</td><td style='text-align: center; word-wrap: break-word;'>48.00</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>Wildcat [^5]</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>84.62</td><td style='text-align: center; word-wrap: break-word;'>84.55</td><td style='text-align: center; word-wrap: break-word;'>90.59</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>45.00</td><td style='text-align: center; word-wrap: break-word;'>44.38</td><td style='text-align: center; word-wrap: break-word;'>84.00</td><td style='text-align: center; word-wrap: break-word;'>46.67</td><td style='text-align: center; word-wrap: break-word;'>84.00</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>Ours</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>43.08</td><td style='text-align: center; word-wrap: break-word;'>17.27</td><td style='text-align: center; word-wrap: break-word;'>32.35</td><td style='text-align: center; word-wrap: break-word;'>57.14</td><td style='text-align: center; word-wrap: break-word;'>61.67</td><td style='text-align: center; word-wrap: break-word;'>17.14</td><td style='text-align: center; word-wrap: break-word;'>0.00</td><td style='text-align: center; word-wrap: break-word;'>16.25</td><td style='text-align: center; word-wrap: break-word;'>32.00</td><td style='text-align: center; word-wrap: break-word;'>51.11</td><td style='text-align: center; word-wrap: break-word;'>62.00</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>Ours</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>✓</td><td style='text-align: center; word-wrap: break-word;'>62.31</td><td style='text-align: center; word-wrap: break-word;'>39.09</td><td style='text-align: center; word-wrap: break-word;'>42.35</td><td style='text-align: center; word-wrap: break-word;'>68.57</td><td style='text-align: center; word-wrap: break-word;'>86.67</td><td style='text-align: center; word-wrap: break-word;'>47.14</td><td style='text-align: center; word-wrap: break-word;'>13.33</td><td style='text-align: center; word-wrap: break-word;'>38.12</td><td style='text-align: center; word-wrap: break-word;'>92.00</td><td style='text-align: center; word-wrap: break-word;'>64.44</td><td style='text-align: center; word-wrap: break-word;'>76.00</td></tr></table>

<div style="text-align: center;">TABLE I: Evaluation scores on HILTI22 SLAM Challenge. Bold: best score, underlined: second best; C: causal evaluation (definition in [^22]); L, V, I: LiDAR, Visual and Inertial measurements are considered; (e/m/h) classifies the difficulty level of the sequence: easy, medium or hard; *based on [^24]. Sequences exp04, exp05 and exp06 greyed out as they are not included in the final score.</div>

53 submissions on the challenge leaderboard. Note that also VILENS and Wildcat performed some form of offline batch or posegraph optimisation for the reported scores. For a better understanding, it can be stated that the reported scores here denote a mean ATE of 1 cm (or even below) to 3 cm in the final estimates. Only for exp07 where a long corridor with a lack of characteristic visual features or geometric features does not achieve the same performance and results in a mean ATE of approximately 7 cm.

### C. Mapping

1) HILTI Dataset: Fig. 1 shows an overlay of all submaps created on Exp21. Meshes of the surfaces and a slice through the occupancy fields at a height of 0.2 m are shown. It qualitatively demonstrates that we are able to provide accurate 3D reconstructions and sensible occupancy fields immediately usable for robot navigation without discernible inconsistencies between submaps.

2) Real-World Example: To qualitatively demonstrate the approach's versatility to a variety of different LiDAR sensor types, we additionally processed a dataset recorded on a Leica BLK2Fly drone. This dataset has been recorded in the outdoor area around an office building. The drone is equipped with five cameras and also five IMUs, out of which we only use one for state estimation. The LiDAR sensor is a single-beam dual-axis spinning LiDAR sensor. Our mapping approach avoids projection into very sparse range images. Fig. 5 shows a reconstruction as well as the trajectory and a slice through the occupancy field.

### D. Timings

We further show timings for map integration and graph optimisation of Exp15 in Table II. LVI Setup 1 here denotes

![Image](images/87b65823cd08f4615b4ab214ea67939a72446ef556b45a3fe2cb5e9ffb41ebbe.jpg)

<div style="text-align: center;">Fig. 5: Outdoor Example: Reconstruction (brown), horizontal slice of the free space (blue) and estimated trajectory (green).</div>

the original setup as it was used for the evaluation of all sequences in I. In Setup 2, we reduced the number of LiDAR residuals to 50 for live frame-to-map and 500 for map-to-map factors. Furthermore, the input point clouds were downsampled with a factor of 3. The results show that the rate of the optimisation can be increased from approximately 10 Hz to almost 15 Hz. Also the map integration can be sped up significantly by downsampling input point clouds for mapping. Instead of running at 10 Hz, we can run integration at 25 Hz. Note that this speed-up only leads to a very minor degradation in terms of localisation accuracy.

<table border=1 style='margin: auto; word-wrap: break-word;'><tr><td style='text-align: center; word-wrap: break-word;'>Function</td><td style='text-align: center; word-wrap: break-word;'>VI only</td><td style='text-align: center; word-wrap: break-word;'>LVI Setup 1</td><td style='text-align: center; word-wrap: break-word;'>LVI Setup 2</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>Optimisation</td><td style='text-align: center; word-wrap: break-word;'>49.0 $\pm$ 16.6</td><td style='text-align: center; word-wrap: break-word;'>98.9 $\pm$ 30.9</td><td style='text-align: center; word-wrap: break-word;'>79.3 $\pm$ 25.0</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>Batch Integration</td><td style='text-align: center; word-wrap: break-word;'>-</td><td style='text-align: center; word-wrap: break-word;'>98.4 $\pm$ 83.0</td><td style='text-align: center; word-wrap: break-word;'>38.2 $\pm$ 34.5</td></tr><tr><td style='text-align: center; word-wrap: break-word;'>ATE (Causal) [cm]</td><td style='text-align: center; word-wrap: break-word;'>23.5</td><td style='text-align: center; word-wrap: break-word;'>2.5</td><td style='text-align: center; word-wrap: break-word;'>3.8</td></tr></table>

<div style="text-align: center;">TABLE II: Per-frame average timings in ms and standard deviations for optimisation and map integration on Exp15.</div>

### VIII. CONCLUSION

Autonomous navigation in the real-world requires high-accuracy localisation, but also accurate and more importantly a globally consistent 3D representation of the environment. In this paper, we propose a fully tightly-coupled LiDAR-Visual-Inertial system that addresses both of these tasks simultaneously. Results from the state estimator are used to integrate incoming LiDAR measurements into local submaps. A novel, correspondence-free residual formulation has been introduced that only uses occupancy field values and gradients which can be efficiently queried. In a tightly-coupled approach, these LiDAR factors can be added to the optimisation problem either as live frame-to-map constraints or as map-to-map constraints across larger time spans. The proposed system has proven to achieve state-of-the-art performance in localisation while yielding consistent occupancy submaps. In future work, we are planning to improve on several aspects. These include amongst others a more realistic uncertainty model for the LiDAR sensor or an informed way of downsampling incoming LiDAR scans leveraging e.g. geometric characteristics. Another challenge, that we will address in the future is the robustness to challenging scenarios in which the visual frontend is prone to tracking failures. Finally, we are also planning to extend the current work to a full exploration framework navigating through submaps in the real-world.

---

## REFERENCES

[^1]: D. Wisth, M. Camurri, and M. Fallon, "Vilens: Visual, inertial, lidar, and leg odometry for all-terrain legged robots," IEEE Transactions on Robotics, vol. 39, no. 1, pp. 309–326, 2022.

[^2]: J. Lin and F. Zhang, "R 3 live: A robust, real-time, rgb-colored, lidar-inertial-visual tightly-coupled state estimation and mapping package," in 2022 International Conference on Robotics and Automation (ICRA). IEEE, 2022, pp. 10672–10678.

[^3]: T. Shan, B. Englot, C. Ratti, and D. Rus, "Lvi-sam: Tightly-coupled lidar-visual-inertial odometry via smoothing and mapping," in 2021 IEEE international conference on robotics and automation (ICRA). IEEE, 2021, pp. 5692–5698.

[^4]: P. Dellenbach, J.-E. Deschaud, B. Jacquet, and F. Goulette, "Cticp: Real-time elastic lidar odometry with loop closure," in 2022 International Conference on Robotics and Automation (ICRA). IEEE, 2022, pp. 5580–5586.

[^5]: M. Ramezani, K. Khosoussi, G. Catt, P. Moghadam, J. Williams, P. Borges, F. Pauling, and N. Kottege, "Wildcat: Online continuous-time 3d lidar-inertial slam," arXiv preprint arXiv:2205.12595, 2022.

[^6]: V. Reijgwart, A. Millane, H. Oleynikova, R. Siegwart, C. Cadena, and J. Nieto, “Voxgraph: Globally consistent, volumetric mapping using signed distance function submaps,” IEEE Robotics and Automation Letters, vol. 5, no. 1, pp. 227–234, 2019.

[^7]: Y. Wang, M. Ramezani, M. Mattamala, S. T. Digumarti, and M. Fallon, "Strategies for large scale elastic and semantic lidar reconstruction," Robotics and Autonomous Systems, vol. 155, p. 104185, 2022.

[^8]: L. Zhang, M. Helmberger, L. F. T. Fu, D. Wisth, M. Camurri, D. Scaramuzza, and M. Fallon, “Hilti-oxford dataset: A millimeter-accurate benchmark for simultaneous localization and mapping,” IEEE Robotics and Automation Letters, vol. 8, no. 1, pp. 408–415, 2022.

[^9]: J. Zhang and S. Singh, “Loam: Lidar odometry and mapping in real-time.” in Robotics: Science and systems, vol. 2, no. 9. Berkeley, CA, 2014, pp. 1–9.

[^10]: T. Shan, B. Englot, D. Meyers, W. Wang, C. Ratti, and D. Rus, "Lio-sam: Tightly-coupled lidar inertial odometry via smoothing and mapping," in 2020 IEEE/RSJ international conference on intelligent robots and systems (IROS). IEEE, 2020, pp. 5135–5142.

[^11]: W. Xu and F. Zhang, “Fast-lio: A fast, robust lidar-inertial odometry package by tightly-coupled iterated kalman filter,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 3317–3324, 2021.

[^12]: W. Xu, Y. Cai, D. He, J. Lin, and F. Zhang, "Fast-lio2: Fast direct lidar-inertial odometry," IEEE Transactions on Robotics, vol. 38, no. 4, pp. 2053–2073, 2022.

[^13]: Z. Liu and F. Zhang, “Balm: Bundle adjustment for lidar mapping,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 3184–3191, 2021.

[^14]: A. I. Mourikis and S. I. Roumeliotis, “A multi-state constraint kalman filter for vision-aided inertial navigation,” in Proceedings 2007 IEEE international conference on robotics and automation. IEEE, 2007, pp. 3565–3572.

[^15]: X. Zuo, P. Geneva, W. Lee, Y. Liu, and G. Huang, "Lic-fusion: Lidar inertial-camera odometry," in 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). IEEE, 2019, pp. 5848–5854.

[^16]: X. Zuo, Y. Yang, P. Geneva, J. Lv, Y. Liu, G. Huang, and M. Pollefeys, "Lic-fusion 2.0: Lidar-inertial-camera odometry with sliding-window plane-feature tracking," in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). IEEE, 2020, pp. 5112–5119.

[^17]: M. Bosse, P. Newman, J. Leonard, M. Soika, W. Feiten, and S. Teller, “An atlas framework for scalable mapping,” in 2003 IEEE International Conference on Robotics and Automation (Cat. No. 03CH37422), vol. 2. IEEE, 2003, pp. 1899–1906.

[^18]: B.-J. Ho, P. Sodhi, P. Teixeira, M. Hsiao, T. Kusnur, and M. Kaess, “Virtual occupancy grid map for submap-based pose graph slam and planning in 3d environments,” in 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). IEEE, 2018, pp. 2175–2182.

[^19]: A. Hornung, K. M. Wurm, M. Bennewitz, C. Stachniss, and W. Burgard, “Octomap: An efficient probabilistic 3d mapping framework based on octrees,” Autonomous robots, vol. 34, pp. 189–206, 2013.

[^20]: Y. Wang, N. Funk, M. Ramezani, S. Papatheodorou, M. Popović, M. Camurri, S. Leutenegger, and M. Fallon, “Elastic and efficient lidar reconstruction for large-scale exploration tasks,” in 2021 IEEE

International Conference on Robotics and Automation (ICRA). IEEE, 2021, pp. 5035–5041.

[^21]: N. Funk, J. Tarrio, S. Papatheodorou, M. Popović, P. F. Alcantarilla, and S. Leutenegger, “Multi-resolution 3d mapping with explicit free space representation for fast and accurate mobile robot motion planning,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 3553–3560, 2021.

[^22]: S. Leutenegger, “Okvis2: Realtime scalable visual-inertial slam with loop closure,” arXiv preprint arXiv:2202.09199, 2022.

[^23]: C. Forster, L. Carlone, F. Dellaert, and D. Scaramuzza, “On-manifold preintegration for real-time visual-inertial odometry,” IEEE Transactions on Robotics, vol. 33, no. 1, pp. 1–21, 2016.

[^24]: C. Zheng, Q. Zhu, W. Xu, X. Liu, Q. Guo, and F. Zhang, "Fast-livo: Fast and tightly-coupled sparse-direct lidar-inertial-visual odometry," 2022.

[^25]: S. Agarwal, K. Mierle, and T. C. S. Team, "Ceres Solver," 3 2022. [Online]. Available: <https://github.com/ceres-solver/ceres-solver>

[^26]: J. Maye, P. Furgale, and R. Siegwart, “Self-supervised calibration for robotic systems,” in 2013 IEEE Intelligent Vehicles Symposium (IV). IEEE, 2013, pp. 473–480.