# Observability-Aware Intrinsic and Extrinsic Calibration of LiDAR-IMU Systems

Jiajun Lv1,†, Xingxing Zuo1,†, Kewei Hu1, Jinhong Xu1, Guoquan Huang2, and
Yong Liu1,3 1The authors are with the Institute of Cyber-System and Control,
Zhejiang University, Hangzhou, China. Email: {lvjiajun314, xingxingzuo,
21932133<EMAIL_ADDRESS><EMAIL_ADDRESS>(Yong Liu is the
corresponding author.)2The author is with the Department of Mechanical
Engineering, University of Delaware, Newark, DE 19716, USA. Email:
<EMAIL_ADDRESS>Institute of Zhejiang University, Huzhou, 313000,
China. †The co-first authors have equal contributions.

###### Abstract

Accurate and reliable sensor calibration is essential to fuse LiDAR and
inertial measurements, which are usually available in robotic applications. In
this paper, we propose a novel LiDAR-IMU calibration method within the
continuous-time batch-optimization framework, where the intrinsics of both
sensors and the spatial-temporal extrinsics between sensors are calibrated
without using calibration infrastructure such as fiducial tags. Compared to
discrete-time approaches, the continuous-time formulation has natural
advantages for fusing high rate measurements from LiDAR and IMU sensors. To
improve efficiency and address degenerate motions, two observability-aware
modules are leveraged: (i) The information-theoretic data selection policy
selects only the most informative segments for calibration during data
collection, which significantly improves the calibration efficiency by
processing only the selected informative segments. (ii) The observability-
aware state update mechanism in nonlinear least-squares optimization updates
only the identifiable directions in the state space with truncated singular
value decomposition (TSVD), which enables accurate calibration results even
under degenerate cases where informative data segments are not available. The
proposed LiDAR-IMU calibration approach has been validated extensively in both
simulated and real-world experiments with different robot platforms,
demonstrating its high accuracy and repeatability in commonly-seen human-made
environments. We also open source our codebase to benefit the research
community: https://github.com/APRIL-ZJU/OA-LICalib.

###### Index Terms:

Sensor calibration, system observability, information matrix, continuous-time
representation, optimization.

##  I Introduction

The LiDAR-IMU (LI) system has prevailed in recent years [1, 2, 3, 4, 5, 6, 7],
as an increasing number of robotic applications require accurate and robust
LI-navigation solutions. Accurate LiDAR-IMU extrinsic calibration is a
prerequisite of LI navigation. Manually measuring the spatial extrinsic
parameters (relative translation and rotation) between the LiDAR and IMU might
be inaccurate or impractical, demanding an easy-to-use, automatic calibration
approach. Besides the extrinsic calibration, calibrating the intrinsics of
both sensors is equally important. IMU intrinsics [8, 9] modeling scale
imperfection and axis misalignment between the gyroscope and accelerometer,
are recommended to calibrate offline [10]. Similarly, LiDAR intrinsics [11,
12] such as each beam’s horizontal and vertical angles and range offset, have
been shown to be better calibrated for accurate registration and mapping
rather than using factory default. Therefore, in this paper, we perform both
intrinsic and extrinsic calibration for LI systems.

Many calibration methods rely on infrastructure (e.g., fiducial targets, or
turntables) to aid calibration [13, 14, 15], which however increases the
technical barrier for the widespread deployment of LI systems. LI sensors used
in practical applications may not be perfectly time synchronized,
necessitating the temporal calibration to compute the time offset between the
two sensors. While discrete-time representation of states is most commonly
used in calibration, its discretization introduces approximation when fusing
asynchronous measurements and incurs a higher computational cost when
processing high-rate data. Building upon our prior work [16], in this work, we
develop a continuous-time LI calibration method that computes both intrinsic
and extrinsic (spatial and temporal) LI sensing parameters without using
dedicated infrastructure.

It is well understood that (infrastructure-free) sensor calibration can be
heavily affected by the environment and motion because of observability (e.g.,
see [4, 10, 17]). In order to make the calibration tool easy-to-use for an end
user, it is desirable to enable our LI calibration method to automatically
address observability issues. To this end, we first perform observability-
aware automatic selection of informative data for offline calibration; that
is, we select only the most informative trajectory segments (rather than an
uninformative long-session trajectory with weakly-excited motions). This is
primarily due to the fact that if a calibration-data-collection trajectory is
not fully excited, some calibration parameters may become unobservable and
thus, cannot be computed or result in inconsistent results [10, 17]. However,
how to collect informative data for calibration is not trivial for non-expert
users. To alleviate this technology barrier, the proposed observability-aware
data selection will automatically select informative trajectory segments,
which in turn significantly reduces the computational cost as less data is
needed for accurate calibration. On the other hand, if degenerate or weakly-
excited motions inevitably occur – for example, in the autonomous-driving
case, straight-line translation without rotation, and single-axis rotation –
we also proposed to perform observability-aware state update during
optimization so that no spurious information would influx into unobservable
parameters. As a result, the proposed observability-aware calibration method,
termed as OA-LICalib, is able to perform efficient and robust LI calibration
with challenging motions and scenes.

In particular, the main contributions of this work include:

  * •

The proposed OA-LICalib is among the first to perform full LI calibration that
is able to calibrate all sensing parameters including the spatial (i.e.,
rigid-body transformation between the sensors) and temporal (i.e., time offset
between the sensors) extrinsics, and the intrinsics of LiDAR and IMU.

  * •

We propose the information-theoretic metrics to select informative trajectory
segments for the LI calibration, instead of naively processing all available
data regardless its observability, which significantly reduces the
computational cost.

  * •

We develop the observability-aware state update during optimization by
updating only the identifiable directions of state space under degenerate
motions (e.g., planar motion and constant velocity).

  * •

Extensive experiments in simulations and real world are conducted to examine
the feasibility and accuracy of the proposed OA-LICalib method, whose codebase
is also open sourced to broadly benefit the robotics community.

The remainder of this paper is organized as follows: After reviewing the
literature in Section II, we provide the background materials of continuous-
time batch optimization and present the sensor models in Section III. The
problem formulation is given in Section IV, while the proposed observability-
aware calibration method is presented in Section V. Simulated and real-world
experiments are carried out in Section VI and Section VII, respectively.

Finally, the paper is concluded in Section VIII.

##  II Related Work

Sensor calibration has a rich literature. Instead of providing a comprehensive
literature review, in this section, we only review the closely related work of
LI calibration.

###  II-A Spatial-temporal LI Calibration

In order to calibrate the rigidly-connected LiDAR and IMU, Geiger et al. [18]
proposed a motion-based calibration approach that performs the extrinsic
calibration by hand-eye calibration [19]. However, their approach expects each
sensor’s trajectory to be estimated independently and accurately, which is
difficult for the consumer-grade IMU. Gentil et al. [20] formulated the LI
calibration as a factor graph optimization problem. Gaussian progress (GP)
regression [21] was adopted to up-sample IMU measurements for pose
interpolation at every time instants of capturing LiDAR points. Although the
IMU data is continuously modeled based on GP, only the states at some specific
time instants are optimized in the discrete-time factor graph, which may
reduce accuracy. Recently, Mishra et al. [22] presented an EKF-based LI
extrinsic calibration method integrated into the visual-inertial navigation
framework OpenVINS [23].

Online extrinsic LI calibration approaches are also available. In [24, 4],
LiDAR, IMU, and camera are tightly fused with a lightweight EKF, where the
extrinsic parameters between the sensors are estimated online. Ye et al. [3]
presented a tightly coupled LI odometry with online extrinsic parameters
estimation, and adopted a linear motion model to eliminate motion distortion
during the sweep. Qiu et al. [25] developed an IMU-centric temporal offset and
extrinsic rotation calibration method for the LI system based on the motion
correlation. Recently, Xu et al. [26] proposed FAST-LIO2, which is a tightly-
coupled LiDAR-inertial odometry with online extrinsic calibration in an
iterated-EKF framework. Note that online calibration often assumes reasonable
initial values in order to guarantee convergence and may have extra
unobservable directions because of including calibration parameters into the
state vector [27, 4].

Continuous-time batch optimization with temporal basis functions is also
widely studied in calibration problems. Furgale et al. [28] detailed the
derivation and realization for a full SLAM problem based on B-spline basis
functions and evaluated the proposed framework within a camera-IMU calibration
problem, which was further extended to support both temporal and spatial
calibration [29]. Rehder et al. [30] adopted a similar framework for
calibrating the extrinsics between a camera-IMU system and a single-beam LiDAR
in two steps: the camera-IMU system is first calibrated with a chessboard, and
then the single-beam LiDAR is calibrated with respect to the camera-IMU
system.

###  II-B Intrinsic Calibration

Li et al. [8] proposed an online state estimation framework using a multi-
state-constraint Kalman filter (MSCKF) which simultaneously estimates
extrinsic parameters of a visual-inertial system, time offset, and intrinsics
of IMU and camera. Rehder et al. [13] extended the previous work [29] to
support the calibration of extrinsics and intrinsics of multiple IMUs in the
visual-inertial system. Yang et al. [10] investigated the necessity of online
IMU intrinsic calibration in the visual-inertial system and performed
observability analysis for different IMU intrinsic models.

As the closest to our work, Liu et al. [14] proposed to calibrate a LI
system’s spatial extrinsics and intrinsics by discrete-time bundle adjustment
optimization. A specially designed target composed of a cone and a cylinder is
required to determine the intrinsics of 3D LiDAR and formulate the constraints
involved with the rigid transformation between LiDAR and IMU. A high-precision
3-axis-adjustable turntable platform with accurate attitude readings is also
necessary to calibrate the IMU intrinsics. A close follow-up work [15]
simplifies the method [14] by removing the dedicated target requirement and
providing the concept to calibrate with the multiple geometric features
naturally existing in human-made environments. Point/sphere, line/cylinder,
and plane features are extracted and leveraged for providing valid
constraints. However, reliable and automatic geometric features extraction
from real-world scenes can be error-prone and susceptible; the extra facility
of a high-precision turntable platform remains necessary for IMU intrinsic
calibration.

###  II-C Observability Awareness in Calibration

It is known that we can assess the information content of trajectory segments
for camera and IMU calibration [9]. As such, Maye et al. [31] proposed
information-theoretic metrics to select informative data and discard redundant
data for calibration. A truncated QR decomposition of the Fisher information
matrix is employed to update estimates in only observable directions. Zhang et
al. [32] analyzed the geometric structures in optimized problem for point
cloud registration and determined well-conditioned directions to solve them
partially. Schneider et al. [9] developed an online self-calibration method in
discrete-time batch optimization for a visual-inertial system with the
extended capability of calibrating IMU intrinsics. Multiple types of
information metrics, including the trace, the determinant, and the maximal
eigen value of the covariance matrix, are investigated for informative segment
selection. The selected informative segments are fully excited with general
motion, making the calibrated parameters observable and solvable. Jiao et al.
[33] designed a sliding window-based multi-LiDAR odometry system with the
capability of online extrinsic calibration between multiple LiDARs, and the
singular values of Hessian matrix are leveraged to examine the convergence of
rotational extrinsic calibration.

###  II-D Extension of Our Previous Conference Publication [16]

While this paper is evolved from our prior work on continuous-time extrinsic
LI calibration of [16], there are significant contributions differentiating
this work from [16] and others. To the best of our knowledge, the proposed OA-
LICalib is among the first to fully calibrate all the sensing parameters of
the LI sensors, including not only the spatial extrinsics as in [16], but the
temporal extrinsics (time offset) as well as the intrinsic parameters of the
LiDAR and IMU. As it is often required to fully excite the sensor platform
during data collection in order for calibration results to converge, which
however might not be possible for some commonly-seen applications such as
self-driving cars [4], in this paper, we have proposed two observability-aware
strategies to address this issue to make the proposed OA-LICalib method more
easy-to-use for non-expert end users. Specifically, (i) the information-
theoretic based data selection policy selects the most informative trajectory
segments for calibration, which eases the effort in calibration data
collection and improves the calibration efficiency. (ii) The observability-
aware state update scheme only updates the states (including calibration
parameters) lying along the observable directions while preventing an influx
of spurious information in the unobservable directions.

##  III Problem Formulation

Figure 1: Illustration of different frames of reference used in the paper. The
trajectory of LiDAR can be represented with the IMU trajectory, the 6-DoF
rigid transformation ${}^{I}_{L}\mathbf{T}$, and the time offset $t_{c}$
between LiDAR and IMU, by
${}^{M}_{L}\mathbf{T}(\tau)=\left({}^{G}_{I}\mathbf{T}(\tau_{0}+t_{c}){}^{I}_{L}\mathbf{T}\right)^{\top}{}^{G}_{I}\mathbf{T}(\tau+t_{c})\
{}^{I}_{L}\mathbf{T}$. Figure 2: The intrinsics of an individual laser
comprising a multi-beam 3D LiDAR. The left figure depicts the ideal model,
while the three figures on the right illustrate the intrinsics arising from
correction factors and offsets.

We first present our notations used throughout the paper. We denote the 6-DoF
rigid transformation by
${}^{B}_{A}{\mathbf{T}}\in\mathit{SE}(3)\subset\mathbb{R}^{4\times 4}$, which
transforms the point ${}^{A}{\mathbf{p}}\in\mathbb{R}^{3}$ in frame $\\{A\\}$
into frame $\\{B\\}$.
${}^{B}_{A}{\mathbf{T}}=\begin{bmatrix}{}^{B}_{A}{\mathbf{R}}&{}^{B}{\mathbf{p}}_{A}\\\
\mathbf{0}&1\end{bmatrix}$ consists of rotational part
${}^{B}_{A}{\mathbf{R}}\in\mathit{SO}(3)$ and translational part
${}^{B}{\mathbf{p}}_{A}\in\mathbb{R}^{3}$. For simplicity, we omit the
homogeneous conversion in the rigid transformation by
${}^{B}{\mathbf{p}}={}^{B}_{A}{\mathbf{T}}{}^{A}{\mathbf{p}}$.
${}^{B}_{A}\bar{q}$ is the quaternion of rotation corresponding to the
rotation matrix ${}^{B}_{A}\mathbf{R}$. In our formulation, IMU reference
frame $\\{I\\}$ is assumed to be rigidly connected to the LiDAR frame
$\\{L\\}$. The map frame $\\{M\\}$ of LiDAR point cloud is determined by the
first LiDAR frame $\\{L_{0}\\}$ when the calibration starts. Similarly, the
global frame $\\{G\\}$ is determined as the first IMU frame $\\{I_{0}\\}$. The
spatial extrinsics, including the relative rotation and translation from the
LiDAR frame to the IMU frame, are denoted by ${}^{I}_{L}\bar{q}$ and
${}^{I}\mathbf{p}_{L}$, respectively. The temporal extrinsic parameter
modeling the time offset between LiDAR and IMU is $t_{c}$. The LiDAR scan
labeled with timestamp $\tau_{L}$ corresponds to IMU time instant
$t_{I}=\tau_{L}+t_{c}$. Figure 1 visualizes the different frames of reference
used in this paper.

The extrinsic and intrinsic LI calibration considered in this paper aims to
estimate the following states:

$\displaystyle\mathcal{X}=\\{\mathbf{x}_{p},\mathbf{x}_{q},\mathbf{x}_{Is},\mathbf{x}_{I},\mathbf{x}_{L},{}^{I}_{L}\bar{q},{}^{I}\mathbf{p}_{L},t_{c}\\}$
(1)

which includes the control points of B-splines
$\\{\mathbf{x}_{p},\mathbf{x}_{q}\\}$ to represent the continuous-time
trajectory, the IMU navigation states $\mathbf{x}_{Is}$, the IMU intrinsics
$\mathbf{x}_{I}$, the LiDAR intrinsics $\mathbf{x}_{L}$, and the LiDAR-IMU
spatiotemporal extrinsic parameters
$\\{{}^{I}_{L}\bar{q},{}^{I}\mathbf{p}_{L},t_{c}\\}$. Based on the LI sensor
data collected during the calibration process, we formulate the following
nonlinear least-squares (NLS) problem:

$\displaystyle\hat{\mathcal{X}}=\underset{\mathcal{X}}{\operatorname{argmin}}\;r,\;\;r=r_{I}+r_{L}\text{,}$
(2)

where $r_{I}$ is the IMU measurement residual cost function which will be
derived in detail in (12) and $r_{L}$ is the LiDAR measurement residual cost
function as derived in (22).

###  III-A Continuous-time Trajectory Representation

We employ B-spline to parameterize trajectory as it provides closed-form
analytic derivatives [34], enabling the effortless fusion of high-frequency
measurements for state estimation. B-spline also has the good property of
being locally controllable, which means the update of a single control point
only impacts certain consecutive segments of the spline [28]. This trait
yields a sparse system with a limited number of control points. To
parameterize the continuous-time 6-DoF trajectory, we leverage the split
representation [35] which represents the 3D translation and the 3D rotation
with uniform B-splines separately. Both the translation and the rotation
splines are parameterized with cumulative form. To be specific, the
translation $\mathbf{p}(t)$ of $d$ degrees over time $t\in[t_{i},t_{i+1})$ is
controlled by the temporally uniformly distributed translational control
points $\mathbf{p}_{i}$, $\mathbf{p}_{i+1}$, $\dots$, $\mathbf{p}_{i+d}$, and
the matrix format could be expressed:

$\displaystyle\mathbf{p}(t)=\mathbf{p}_{i}+\sum_{j=1}^{d}\mathbf{u}^{\top}{\tilde{\mathbf{M}}}^{(d+1)}_{(j)}\left(\mathbf{p}_{i+j}-\mathbf{p}_{i+j-1}\right)$
(3)

where $\mathbf{u}^{\top}=\begin{bmatrix}1&u&\dots&u^{d}\end{bmatrix}$ and
$u=(t-t_{i})/(t_{i+1}-t_{i})$. $\tilde{\mathbf{M}}^{(4)}_{(j)}$ is the j-th
column of the cumulative spline matrix $\tilde{\mathbf{M}}^{(4)}$ which only
depends on the corresponding degree of uniform B-spline. In this paper, cubic
($d=3$) B-spline is employed, thus the corresponding cumulative spline matrix
is

$\displaystyle\tilde{\mathbf{M}}^{(4)}=\frac{1}{6}\begin{bmatrix}6&5&1&0\\\
0&3&3&0\\\ 0&-3&3&0\\\ 0&1&-2&1\\\ \end{bmatrix}\text{.}$ (4)

Besides, we adopt the cumulative representation of B-spline to parameterize
the rotation on $\mathit{SO}(3)$ like [36, 35, 37], and quaternions are served
as the rotational control points of B-splines [38]:

$\displaystyle\bar{q}(t)=\bar{q}_{i}\otimes\prod_{j=1}^{d}exp\left(\mathbf{u}^{\top}{\tilde{\mathbf{M}}}^{(4)}_{(j)}log(\bar{q}_{i+j-1}^{-1}\otimes\bar{q}_{i+j})\right)$
(5)

where $\otimes$ denotes quaternion multiplication. $\bar{q}_{i}$ denotes the
rotational control point. $exp\left(\cdot\right)$ is the operation that
mapping Lie algebra elements to $\mathit{S}^{3}$, and $log\left(\cdot\right)$
is its inverse operation [38].

In this calibration system, the continuous 6-DOF poses of IMU
${}_{I}^{G}\bar{q}(t),{}^{G}{\mathbf{p}}_{I}(t)$ are represented in the global
frame $\\{G\\}$ by splines in the format of (3) and (5) with parameters of
control points $\mathbf{x}_{p}$ and $\mathbf{x}_{q}$ for 3D translation and
rotation, respectively. The derivatives of the splines with respect to time
can be easily computed [37], which lead to linear accelerations
${}^{I}\mathbf{a}(t)$ and angular velocities ${}^{I}\bm{\omega}(t)$ in the
local IMU reference frame:

$\displaystyle{}^{I}\mathbf{a}(t)$
$\displaystyle={}_{I}^{G}\mathbf{R}^{\top}(t)\left({}^{G}\ddot{\mathbf{p}}_{I}(t)-{}^{G}\mathbf{g}\right)$
(6) $\displaystyle{}^{I}\bm{\omega}(t)$
$\displaystyle={}_{I}^{G}\mathbf{R}^{\top}(t){}_{I}^{G}\dot{\mathbf{R}}(t)$
(7)

where ${}^{G}\mathbf{g}\in\mathbb{R}^{3}$ denotes the gravity vector in global
frame. ${}^{G}\mathbf{g}$ has only two degrees of freedom since its norm is
assumed to be a constant $\|{}^{G}\mathbf{g}\|\simeq 9.8$ in the system.

###  III-B IMU Model

IMU comprises a 3-axis gyroscope and a 3-axis accelerometer. Inspired by [8,
13, 10], the IMU measurements are modeled as:

$\displaystyle{}^{\omega}\bm{\omega}$
$\displaystyle=\mathbf{S}_{\omega}\mathbf{M}_{\omega}\
{}^{\omega}_{I}\mathbf{R}{}^{I}\bm{\omega}(t)+\mathbf{A}_{\omega}{}^{I}\mathbf{a}(t)+\mathbf{b}_{\omega}+\mathbf{n}_{\omega}$
$\displaystyle=h_{I_{\omega}}(t,\mathbf{x}_{q},^{G}\mathbf{g},\mathbf{S}_{\omega},\mathbf{M}_{\omega},{}^{\omega}_{I}\mathbf{R},\mathbf{A}_{\omega},\mathbf{b}_{\omega})$
(8) $\displaystyle{}^{a}\mathbf{a}$
$\displaystyle=\mathbf{S}_{a}\mathbf{M}_{a}\
{}^{a}_{I}\mathbf{R}{}^{I}\mathbf{a}(t)+\mathbf{b}_{a}+\mathbf{n}_{a}$
$\displaystyle=h_{I_{a}}(t,\mathbf{x}_{p},\mathbf{x}_{q},^{G}\mathbf{g},\mathbf{S}_{a},\mathbf{M}_{a},{}^{a}_{I}\mathbf{R},\mathbf{b}_{a})$
(9)

where $\mathbf{x}_{p}$ and $\mathbf{x}_{q}$ are the translational and
rotational control points of B-splines, respectively.
$\mathbf{b}_{\omega},\mathbf{b}_{a}$ are the biases of gyroscope and
accelerometer, which are assumed to be under Gaussian random walk.
$\mathbf{n}_{\omega}$ and $\mathbf{n}_{a}$ are the zero-mean white Gaussian
noise with covariance matrix ${\Sigma_{\omega}}$ and ${\Sigma_{a}}$,
respectively, which can be obtained from IMU noise model without approximation
or propagation. $\mathbf{S}_{\omega/a}$ is a diagonal matrix modeling the
scale imperfection as

$\displaystyle\mathbf{S}_{\omega}=\begin{bmatrix}S_{w1}&0&0\\\ 0&S_{w2}&0\\\
0&0&S_{w3}\\\ \end{bmatrix}\text{,}\
\mathbf{S}_{a}=\begin{bmatrix}S_{a1}&0&0\\\ 0&S_{a2}&0\\\ 0&0&S_{a3}\\\
\end{bmatrix}\text{.}$ (10)

$\mathbf{M}_{\omega/a}$ accounting for axis misalignment is an upper-
triangular matrix with diagonal elements as identity

$\displaystyle\mathbf{M}_{\omega}=\begin{bmatrix}1&M_{w1}&M_{w2}\\\
0&1&M_{w3}\\\ 0&0&1\\\ \end{bmatrix}\text{,}\
\mathbf{M}_{a}=\begin{bmatrix}1&M_{a1}&M_{a2}\\\ 0&1&M_{a3}\\\ 0&0&1\\\
\end{bmatrix}\text{.}$ (11)

$\mathbf{A}_{\omega}$ is a full $3\times 3$ matrix accounting for the
acceleration dependence (g-sensitivity) of the measurements.

Since the gyroscope $\\{\omega\\}$ and accelerometer $\\{a\\}$ inside the IMU
are individual sensors, there is a misalignment between them and the base IMU
frame $\\{I\\}$. ${}^{\omega/a}_{I}\mathbf{R}$ compensates for the rotational
misalignment, and the translational misalignment is omitted for calibration
since it is close to zero for single-chip MEMS sensors [39, 9, 10]. We can
only calibrate ${}^{\omega}_{I}\mathbf{R}$ or ${}^{a}_{I}\mathbf{R}$ since the
base IMU frame $\\{I\\}$ is aligned with either $\\{\omega\\}$ or $\\{a\\}$,
and calibrating both makes the system unobservable [10]. In this paper, we
choose to calibrate ${}^{\omega}_{I}\mathbf{R}$ only. Additionally, we ignore
the calibration of gravity sensitivity in analogy to [10].

Usually, IMU measurements are given at discrete time instants. The
measurements at timestamp $t_{k}$ are denoted as
${}^{\omega_{k}}\bm{\omega}_{m},{}^{a_{k}}\mathbf{a}_{m}$. Since the
trajectory is formulated in continuous time, we can readily get
${}^{I}\mathbf{a}(t_{k}),{}^{I}\bm{\omega}(t_{k})$ from (6) and (7) by
computing differentials of the trajectory. The IMU induces cost functions,
$r_{I}$, related to estimated variables including the control points of the
trajectory $\\{\mathbf{x}_{p},\mathbf{x}_{q}\\}$, the IMU navigation states
$\mathbf{x}_{Is}=\\{^{G}\mathbf{g},\mathbf{b}_{\omega},\mathbf{b}_{a}\\}$, and
the IMU intrinsics
$\mathbf{x}_{I}=\\{\mathbf{S}_{\omega},\mathbf{M}_{\omega},{}^{\omega}_{I}\mathbf{R},\mathbf{S}_{a},\mathbf{M}_{a}\\}$
are formulated as:

$\displaystyle r_{I}=r_{\omega}+r_{a}\text{,}\quad\text{where}$ (12)
$\displaystyle
r_{\omega}=\sum_{k}\frac{1}{2}\|{}^{\omega_{k}}\bm{\omega}_{m}\\!-\\!h_{I_{\omega}}(t_{k},\mathbf{x}_{q},\mathbf{x}_{I},\mathbf{x}_{Is})\|^{2}_{\Sigma_{\omega}}\text{,}$
(13) $\displaystyle
r_{a}=\frac{1}{2}\|{}^{a_{k}}\mathbf{a}_{m}\\!-\\!h_{I_{a}}(t_{k},\mathbf{x}_{p},\mathbf{x}_{q},\mathbf{x}_{I},\mathbf{x}_{Is})\|^{2}_{\Sigma_{a}}\text{,}$
(14)

$\|\mathbf{e}\|^{2}_{\Sigma}=\mathbf{e}^{\top}\Sigma^{-1}\mathbf{e}$ denotes
the squared energy norm weighted by the inverse of the covariance matrix
$\Sigma$, and $r_{\omega},r_{a}$ are residuals of angular velocity and linear
acceleration according to (8) and (9), respectively.

###  III-C 3D LiDAR Model

A mechanical 3D LiDAR with multiple laser beams [40] measures the ranges given
by individual laser heads pointing to different elevation angles. A LiDAR scan
is captured when the rigidly connected laser heads rotate around a mechanical
central axis. A laser head $i$ with elevation angle $\phi_{i}$ at LiDAR time
instant $\tau_{k}$ (with azimuth angle $\theta_{ik}$) gets a range measurement
$\rho_{ik}$, then the LiDAR point measurement transferred from the spherical
coordinate system can be denoted by:

$\displaystyle{}^{L_{k}}\mathbf{p}_{ik}=\begin{bmatrix}{}^{L_{k}}x_{ik}\\\
{}^{L_{k}}y_{ik}\\\
{}^{L_{k}}z_{ik}\end{bmatrix}=\begin{bmatrix}\rho_{ik}\cos\phi_{i}\cos\theta_{ik}\\\
\rho_{ik}\cos\phi_{i}\sin\theta_{ik}\\\
\rho_{ik}\sin\phi_{i}\end{bmatrix}\text{.}$ (15)

However, the above ideal model does not hold in practical manufacture since
spatial offsets exist between laser heads, and non-negligible errors exist in
both the range measurement and the azimuth angle $\theta_{ik}$. As shown in
Fig. 2, the intrinsics of an individual laser head $i$ consist of:

  * •

Elevation angle correction factor $\delta\phi_{i}$ and azimuth angle
correction factor $\delta\theta_{i}$.

  * •

Vertical spatial offset $V_{i}$ and horizontal spatial offset $H_{i}$.

  * •

Scale factor of range measurement $s_{i}$ and range measurement offset
$\delta\rho_{i}$.

Given the intrinsics, the zero-mean Gaussian measurement noise $n_{\rho,ik}$
over the range measurement, and defining the following factors:

$\left\\{\begin{array}[]{c}\bar{\phi}_{i}=\phi_{i}+\delta\phi_{i}\\\
\bar{\theta}_{ik}=\theta_{ik}+\delta\theta_{i}\\\
\bar{\rho}_{ik}=s_{i}\rho_{ik}+\delta\rho_{i}+n_{\rho,ik}\end{array}\right.$
(16)

The position of LiDAR point measurement is given by:

$\displaystyle{}^{L_{k}}\mathbf{p}_{ik}=\begin{bmatrix}\bar{\rho}_{ik}\cos\bar{\phi}_{i}\cos\bar{\theta}_{ik}+H_{i}\sin\bar{\theta}_{ik}\\\
\bar{\rho}_{ik}\cos\bar{\phi}_{i}\sin\bar{\theta}_{ik}+H_{i}\cos\bar{\theta}_{ik}\\\
\bar{\rho}_{ik}\sin\bar{\phi}_{i}+V_{i}\end{bmatrix}\text{.}$ (17)

The intrinsics of a 3D LiDAR with $l$ beams are denoted by

$\displaystyle\mathbf{x}_{L}=\\{\delta\phi_{i},\delta\theta_{i},V_{i},H_{i},s_{i},\delta\rho_{i}\\}_{|_{i=0,1,\cdots,l-1}}\text{.}$
(18)

For a LiDAR point measurement ${}^{L_{k}}\mathbf{p}_{ik}$ with an associated
3D plane with closest point parameterization [41],
${}^{M}\bm{\pi}_{j}={}^{M}d_{\pi,j}{}^{M}\mathbf{n}_{\pi,j}$, where
${}^{M}d_{\pi,j}$ and ${}^{M}\mathbf{n}_{\pi,j}$ denote the distance of the
plane to origin and unit normal vector expressed in map frame, respectively.
The point to plane distance is given by:

$\displaystyle{}^{M}\mathbf{p}_{ik}$
$\displaystyle={}^{M}_{L_{k}}\mathbf{R}(\tau_{k}){}^{L_{k}}\mathbf{p}_{ik}+{}^{M}\mathbf{p}_{L_{k}}(\tau_{k})$
(19) $\displaystyle\mathbf{z}_{ijk}$
$\displaystyle={}^{M}\mathbf{n}^{\top}_{\pi,j}{}^{M}\mathbf{p}_{ik}+{}^{M}d_{\pi,j}$
$\displaystyle=h_{L}(\tau_{k},\mathbf{x}_{p},\mathbf{x}_{q},\mathbf{x}_{L},{}^{I}_{L}\bar{q},{}^{I}\mathbf{p}_{L},t_{c})\text{.}$
(20)

The distance measurement above is related to variables including the control
points of trajectory $\\{\mathbf{x}_{p},\mathbf{x}_{q}\\}$, LiDAR intrinsics
$\mathbf{x}_{L}$, and spatial-temporal parameters
$\\{{}^{I}_{L}\bar{q},{}^{I}\mathbf{p}_{L},t_{c}\\}$. The 3D plane is
extracted from LiDAR point cloud map, detailed in Sec.IV-B, and assumed free
of noise when computing $\mathbf{z}_{ijk}$. Thus, ${\Sigma_{z_{ijk}}}$, the
covariance of the distance measurement $\mathbf{z}_{ijk}$, can be easily
propagated from the Gaussian noise $n_{\rho,ik}$ by first-order approximation.
The relative pose
$\\{{}^{M}_{L_{k}}\mathbf{R}(\tau_{k}),{}^{M}\mathbf{p}_{L_{k}}(\tau_{k})\\}$
in (20) is involved with the following equations from parameterized IMU
continuous-time trajectory:

$\displaystyle{}^{M}_{L_{k}}\mathbf{T}(\tau_{k})=\left({}^{G}_{I}\mathbf{T}(\tau_{0}+t_{c}){}^{I}_{L}\mathbf{T}\right)^{\top}{}^{G}_{I}\mathbf{T}(\tau_{k}+t_{c})\
{}^{I}_{L}\mathbf{T}\text{.}$ (21)

It is easy to find that the time offset $t_{c}$ between LiDAR and IMU is also
involved in the process of fetching pose ${}^{M}_{L_{k}}\mathbf{T}(\tau_{k})$
from continuous-time trajectory parameterized in the IMU time axis. In
practice, the LiDAR point measurement at timestamp $\tau_{k}$ is denoted as
$\mathbf{z}_{m,ijk}$. By minimizing the point to plane distance, we have the
following cost functions from LiDAR points associated with 3D planes:

$\displaystyle
r_{L}\\!=\\!\scalebox{0.9}{\mbox{$\displaystyle\frac{1}{2}\sum_{i}\sum_{j}\sum_{k}\|\mathbf{z}_{m,ijk}-h_{L}(\tau_{k},\mathbf{x}_{p},\mathbf{x}_{q},\mathbf{x}_{L},{}^{I}_{L}\bar{q},{}^{I}\mathbf{p}_{L},t_{c})\|^{2}_{\Sigma_{z_{m,ijk}}}$}}$
(22)

where $i,j,k$ denote the indexes of LiDAR beams, 3D planes and points in a
scanned ring by a laser beam, respectively.

Figure 3: The pipeline of proposed LiDAR-IMU calibration method, which allows
leveraging all the raw measurements from IMU and LiDAR sensor in a continuous-
time batch optimization framework. Details are provided in Section IV and
Section V.

##  IV Continuous-time Full LI Calibration

In this section, we present in detail the proposed continuous-time full LI
calibration method. Before diving into details, we first provide an overview
of the overall architecture as shown in Fig. 3. Specifically, we first
initialize the control points of rotation B-splines with the IMU angular
velocities, and the extrinsic rotation ${}^{I}_{L}\bar{q}$ by aligning the IMU
rotations with the LiDAR odometry (rotations) computed with NDT registration
[42]. By reducing the rotational distortion effect in LiDAR scans, we re-
compute the LiDAR odometry to improve its accuracy. We then initialize the
LiDAR surfels by fusing LiDAR scans with the computed LiDAR odometry and find
the point-to-surfel correspondences. With the LI measurements, we perform
batch optimization to estimate the considered states including B-splines
control points and spatial-temporal extrinsics, whose optimal estimates are
used to refine the surfels map and point-to-plane data association.

###  IV-A State Initialization

Initial guesses of the calibrated parameters are needed for the iterative non-
linear least-squares optimization. For initialization, we set IMU navigation
states as
$\hat{\mathbf{x}}_{{Is}_{0}}=\\{\|^{G}\mathbf{g}\|\mathbf{e}_{3},\mathbf{0}_{3},\mathbf{0}_{3}\\}$
and IMU intrinsics as
$\hat{\mathbf{x}}_{I_{0}}=\\{\mathbf{I}_{3},\mathbf{I}_{3},\mathbf{I}_{3},\mathbf{I}_{3},\mathbf{I}_{3}\\}$,
where $\mathbf{e}_{3}=[0,0,1]^{\top}$, $\mathbf{I}_{3}$ and $\mathbf{0}_{3}$
represent the $3\times 3$ identity matrix and zero matrix, respectively. As
for the initial values of LiDAR intrinsics, $\hat{\mathbf{x}}_{L_{0}}$ is
initialized as $\\{0,0,0,0,1,0\\}_{i}{|_{i=1,\cdots,l}}$. The time offset
$t_{c}$ between LiDAR and IMU is initialized as zero, which means the
timestamp $t_{k}$ of IMU is equal to the timestamp $\tau_{k}$ of LiDAR at the
beginning.

The initialization of control points and extrinsic parameters are described in
detail below. Prior to detailing the initialization process, we clarify the
inputs for initialization including discrete poses at about 10 Hz from LiDAR
odometry, angular velocities and linear accelerates at typically 400 Hz from
IMU. We adopt NDT registration to implement simple scan-to-map LiDAR odometry,
sufficient for sequences with a few seconds, and build a global map consists
of keyframed LiDAR scans determined according to the pose changes.

####  IV-A1 Initialize Control Points of Continuous-time Trajectory

Before detailing the concrete initialization method, we emphasize that the
spline trajectory is employed to represent the IMU trajectory rather than
LiDAR’s. If we model the continuous-time trajectory in LiDAR sensor, the
predicted linear acceleration and angular velocity in IMU frame need to be
transformed from LiDAR frame as

$\displaystyle{}^{I}{\mathbf{a}}(t)$
$\displaystyle={}^{I}_{L}{\mathbf{R}}({}^{M}_{L}{\mathbf{R}}(t)^{\top}\left({}^{M}{\mathbf{a}}(t)-{}^{M}{\mathbf{g}}\right)$
$\displaystyle+\lfloor{{}^{L}\dot{\bm{\omega}}(t)}\rfloor{}^{I}{\mathbf{p}}_{L}+\lfloor{{}^{L}{\bm{\omega}(t)}}\rfloor^{2}\
{}^{I}{\mathbf{p}}_{L})$ (23) $\displaystyle{}^{I}\bm{\omega}(t)$
$\displaystyle={}^{I}_{L}{\mathbf{R}}{}^{L}\bm{\omega}(t)\text{,}$ (24)

where $\lfloor\cdot\rfloor$ represents the skew-symmetric matrix of a 3D
vector. Note that angular accelerations ${}^{L}\dot{\bm{\omega}}(t)$ are
required to be computed, which makes the system complicated. Therefore, we
choose to parameterize the trajectory in IMU frame instead of LiDAR frame.

Given the raw angular velocity measurements from the IMU sensor, a series of
rotational control points $\mathbf{x}_{q}$ of the rotational trajectory can be
initialized by solving the following least-square problem:

$\displaystyle\hat{\mathbf{x}}_{q}\\!$
$\displaystyle=\\!\underset{\mathbf{x}_{q}}{\operatorname{argmin}}\;r_{\omega}$
$\displaystyle=\\!\underset{\mathbf{x}_{q}}{\operatorname{argmin}}\;\sum_{k}\frac{1}{2}\|{}^{\omega_{k}}\bm{\omega}_{m}\\!-\\!h_{I_{\omega}}(t_{k},\mathbf{x}_{q},\mathbf{x}_{I},\mathbf{x}_{Is})\|^{2}_{\Sigma_{\omega}}\text{,}$
(25)

where $r_{\omega}$ is defined in (13) and here we ignore the intrinsics of IMU
and the angular velocity bias of IMU. It is important to note that
${}^{G}_{I}\bar{q}(t_{0})$ is fixed to the identity quaternion during the
optimization. In (IV-A1), we try to initialize the rotational trajectory by
the raw IMU measurements rather than the integrated IMU poses, since the
latter are inaccurate and always affected by drifting IMU biases and noises.

For the initialization of translational control points $\mathbf{x}_{p}$, we
first utilize the LiDAR odometry to get descrete-time LiDAR poses. Then we
formulate and solve the following linear problem:

$\displaystyle\hat{\mathbf{x}}_{p}\\!=\\!\underset{\mathbf{x}_{p}}{\operatorname{argmin}}\sum_{k}\|{}^{G}_{M}\mathbf{T}{}^{M}\hat{\bm{p}}_{L_{k}}\\!-\\!{}^{G}_{I}\hat{\mathbf{R}}(t_{k},\hat{\mathbf{x}}_{q}){}^{I}\bm{p}_{L}\\!-\\!{}^{G}\bm{p}_{I}(t_{k},\mathbf{x}_{p})\|\text{,}$
(26)

where ${}^{M}\hat{\bm{p}}_{L_{k}}$ is the estimated position of scan at
timestamp $\tau_{k}$ from the LiDAR odometry;
${}^{G}_{I}\hat{\mathbf{R}}(t_{k},\hat{\mathbf{x}}_{q})$ controlled by
$\hat{\mathbf{x}}_{q}$ has been initialized in (IV-A1) and
${}^{G}\bm{p}_{I}(t_{k},\mathbf{x}_{p})$ controlled by $\mathbf{x}_{p}$ is the
translational trajectory to be initialized. Considering map frame and global
frame are determined as the first measurements from LiDAR and IMU,
respectively, ${}^{G}_{M}\mathbf{T}$ is equivalent to
${}^{I_{0}}_{L_{0}}\mathbf{T}$, that is ${}^{I}_{L}\mathbf{T}$. The initial
values of ${}^{I}_{L}\mathbf{T}$ computed in Sec. IV-A2 are used here. This
initialization method of control points is relatively coarse, but it can
effectively improve the convergence speed of the optimization problem,
compared to leave $\mathbf{x}_{p}$ uninitialized as zeros.

####  IV-A2 Extrinsic Initialization

Inspired by [43], we initialize the extrinsic rotation by aligning two
rotation sequences from LiDAR and IMU. With discrete LiDAR poses from LiDAR
odometry, it is easy to get the relative rotation between two consecutive
LiDAR scans, ${}^{L_{k}}_{L_{k+1}}\bar{q}$. Besides, the relative rotation
between time interval $[t_{k},t_{k+1}]$ in the IMU frame can also be obtained
from the initialized rotational trajectory as
${}_{I_{k+1}}^{I_{k}}\bar{q}={}_{I}^{G}\bar{q}^{-1}(t_{k})\otimes{}_{I}^{G}\bar{q}(t_{k+1})$.
The relative rotations at any $k$ from two sensors should satisfy the
following equation:

$\displaystyle{}^{I_{k}}_{I_{k+1}}\bar{q}\otimes{}_{L}^{I}\bar{q}={}_{L}^{I}\bar{q}\otimes{}^{L_{k}}_{L_{k+1}}\bar{q}\text{.}$
(27)

The above equation can be transferred into another equivalent representation
[44]:

$\displaystyle\mathcal{L}({}^{I_{k}}_{I_{k+1}}\bar{q}){}^{I}_{L}\bar{q}=\mathcal{R}({}^{L_{k}}_{L_{k+1}}\bar{q}){}^{I}_{L}\bar{q}$
$\displaystyle\Rightarrow$
$\displaystyle\left(\mathcal{L}({}^{I_{k}}_{I_{k+1}}\bar{q})-\mathcal{R}({}^{L_{k}}_{L_{k+1}}\bar{q})\right){}^{I}_{L}\bar{q}=\mathbf{0}\text{.}$
(28)

For a unit quaternion
$\bar{q}=\begin{bmatrix}\mathbf{q}_{v}&q_{w}\end{bmatrix}^{\top}$ with
imaginary part $\mathbf{q}_{v}$ and real part $q_{w}$, $\mathcal{L}(\bar{q})$
and $\mathcal{R}(\bar{q})$ are multiplication matrices defined as follow:

$\displaystyle\mathcal{L}({\bar{q}})$
$\displaystyle=\begin{bmatrix}q_{w}\mathbf{I}_{3}-\lfloor\mathbf{q}_{v}\rfloor&\mathbf{q}_{v}\\\
-\mathbf{q}^{T}_{v}&q_{w}\end{bmatrix}$ $\displaystyle\mathcal{R}({\bar{q}})$
$\displaystyle=\begin{bmatrix}q_{w}\mathbf{I}_{3}+\lfloor\mathbf{q}_{v}\rfloor&\mathbf{q}_{v}\\\
-\mathbf{q}^{T}_{v}&q_{w}\end{bmatrix}\text{.}$

Equation (28) is in a form of
${}^{k}_{k+1}\mathbf{A}{}^{I}_{L}\bar{q}=\mathbf{0}$. By stacking a series of
relative rotation measurements ${}^{k}_{k+1}\mathbf{A}$ row by row into a big
matrix $\mathbf{A}$, we have $\mathbf{A}{}^{I}_{L}\bar{q}=\mathbf{0}$, then
${}^{I}_{L}\hat{\bar{q}}$ can be computed as the right unit singular vector
corresponding to the smallest singular value of $\mathbf{A}$.

For the initialization of extrinsic translation ${}^{I}\bm{p}_{L}$, the
initial values of ${}^{I}\bm{p}_{L}$ are simply set to zeros, and the proposed
method is able to converge to reasonable values in practical experiments when
enough motion excitation exists.

(a) The surfels map in the first iteration.

(b) The surfels map in the second iteration.

Figure 4: Typical surfels map in an indoor scenario. The initial surfels map
suffers from motion distortion with the initialized LiDAR poses. With the
estimated states after just one iteration of batch optimization, the quality
of refined surfels map becomes significantly better. Points in different
surfels are colorized differently.

###  IV-B Data Association of LiDAR Measurements

With the estimated rotational extrinsics ${}^{I}_{L}\hat{\bar{q}}$, we can
transfer the rotational continuous-time curve in IMU frame into LiDAR frame
and remove the rotational motion distortion in LiDAR scans. Subsequently, we
perform the NDT-based LiDAR odometry again with the initial rotational guess
retrieved from the rotational continuous-time curve transformed. Since LiDAR
scans are partially undistorted and better initial guesses are provided, NDT-
based LiDAR odometry can generate more accurate estimated poses. With the
odometry results, a LiDAR point cloud map can be constructed by assembling
multiple LiDAR scans into an identical frame. We further uniformly discretize
the LiDAR point cloud map into 3D cells and compute the point distribution in
every 3D cell to identify surfels based on a plane-likeness coefficient [45]:

$\mathcal{P}=2\frac{\lambda_{1}-\lambda_{0}}{\lambda_{0}+\lambda_{1}+\lambda_{2}}\text{,}$
(29)

In the above equation, $\lambda_{0}\leq\lambda_{1}\leq\lambda_{2}$ are the
point distribution eigen values in a 3D cell. For the cell with planar point
distribution, $\mathcal{P}$ will be close to 1. The cell with plane-likeness
coefficient $\mathcal{P}$ larger than a threshold is determined as a surfel
cell, and the points inside this cell are used to fitting a plane. The
estimated plane is parameterized by the closest point as mentioned in Sec.
III-C. Some LiDAR points are associated with surfels if the point-to-plane
distance is less than a threshold, whereas the others are discarded. As
mentioned in Sec. III-C, the plane parameters are assumed to be perfect and
keep constant when we perform least-squares optimization, so are the surfels
map and the data associations. This is in a fashion of
Expectation–maximization algorithm. After the optimization, the surfels map,
the data associations, and the plane parameters will be updated with updated
estimated states, which will be detailed in Sec. IV-C. Additionally, two
practical tricks are engaged to repress the effect of outlier data
associations on optimization: firstly, only the point and plane correspondence
with a distance below 5 centimeters will be taken into the continuous-time
batch optimization; secondly, we rely on the M-estimation [46] idea and apply
the Huber robust kernel to deal with outliers.

###  IV-C Refinement

After the continuous-time batch optimization, the estimated states including
the extrinsics become more accurate. Thus, we leverage the current best
estimates to remove the motion distortion in the raw LiDAR scans, rebuild the
LiDAR surfels map, and update the point-to-surfel data associations. Note that
results from NDT-based LiDAR odometry are only utilized at the very beginning
for initializing the LiDAR poses and LiDAR map. The refinement step with
iterative batch-optimization will generate more accurate estimations and
reliable data association. The typical LiDAR surfels maps in the first and
second iterations of the batch optimization are shown in Fig. 4. The quality
of the map can be improved in a significant margin after one iteration in the
batch optimization.

As for the intrinsic calibration, a 16-beam LiDAR contains 96 parameters (6
parameters per beam), and some parameters of the first laser are fixed as a
reference during the calibration process. Only the distance offset and
distance scale are estimated for the first laser. Additionally, there are 15
parameters of IMU intrinsic parameters. Considering the large number of
intrinsic parameters to be estimated, to prevent the system from falling into
local optimum and reduce the complexity of the problem, we optimize the
spatial-temporal parameters only and fix the intrinsics at the initial values
firstly, and allow optimizing the intrinsic parameters a few iterations later
when the estimates of the trajectory and extrinsic parameters converge.
Furthermore, we correct the raw measurements of IMU and LiDAR based on the
estimated intrinsic parameters after they have gradually converged, instead of
using the poor intrinsic estimates to correct the measurements at the
beginning, which may hurt the calibration results. Throughout the experiments,
we start to calibrate the intrinsic parameters after 2 iterations of batch
optimization and correct raw measurements based on the estimated intrinsic
parameters after 11 iterations.

##  V Observability-Aware LI Calibration

In this section, we describe in detail the observability-aware strategies used
in the proposed LI calibration in order to improve its robustness to
degenerate motions and thus make it more easy-to-use for practitioners.

###  V-A Information-Theoretic Data Selection

Note that some data segments collected without sufficient motion excitation or
scene constraints may degrade the calibration [47], and blindly processing all
the available data may result in unnecessarily high computational cost. As
such, the proposed method seeks to automatically select the most informative
data segments for calibration. In particular, to facilitate the ensuing
derivations, assuming the white Gaussian noise and independent measurements,
the NLS problem (2) is equivalent to the maximum likelihood estimation (MLE)
problem and can be re-formulated as follows (see [48, 49]):

$\displaystyle\hat{\mathbf{x}}=\underset{\mathbf{x}}{\operatorname{argmin}}\sum_{i}\frac{1}{2}\|\mathbf{e}_{i}(\mathbf{x})\|^{2}_{\Sigma_{i}}$
(30)

where $\mathbf{e}_{i}$ denotes the residual induced from the $i$-th
measurement with covariance ${\Sigma_{i}}$. An iterative method such as
Gaussian-Newton and Levenberg-Marquard is often used to solve this problem.
That is, at each iteration, by linearizing the above nonlinear system at the
current estimate $\hat{\mathbf{x}}$,
$\mathbf{e}_{i}(\mathbf{x})=\mathbf{e}_{i}(\hat{\mathbf{x}})+\mathbf{J}_{i}\delta\mathbf{x}$,
with the Jacobian matrix $\mathbf{J}_{i}$, we have the following linear least-
squares problem with respect to the update increment (or error state estimate)
$\delta\hat{\mathbf{x}}$:

$\displaystyle\delta\hat{\mathbf{x}}=\underset{\delta\mathbf{x}}{\operatorname{argmin}}\sum_{i}\frac{1}{2}\|\mathbf{e}_{i}(\hat{\mathbf{x}})+\mathbf{J}_{i}\delta\mathbf{x}\|^{2}_{\Sigma_{i}}$
(31)

The optimal solution is given by the following normal equation:

$\displaystyle\underbrace{\left(\sum_{i}\mathbf{J}^{\top}_{i}{\Sigma^{-1}_{i}}\mathbf{J}_{i}\right)}_{\mathbf{A}}\delta\mathbf{x}$
$\displaystyle=\underbrace{-\sum_{i}\mathbf{J}_{i}^{\top}{\Sigma^{-1}_{i}}\mathbf{e}_{i}(\hat{\mathbf{x}})}_{\mathbf{b}}$
(32)

Once $\delta\hat{\mathbf{x}}$ is computed, we update the estimate as
$\hat{\mathbf{x}}\leftarrow\hat{\mathbf{x}}\boxplus\delta\mathbf{x}$, where
$\boxplus$ denotes the addition operation on manifold.

Note that $\mathbf{A}_{k\times k}$ in (32) is the Fisher information matrix of
the MLE (or NLS) problem (30) and encapsulates all the information contained
in the measurements. It is also important to note that $\mathbf{A}$ might be
rank deficient, for example, due to the lack of motion excitation and/or
environmental constraints (see [4]). It thus is critical to collect
informative data ensuring full-rank $\mathbf{A}$ (or observability [50]). To
this end, we first partition the entire collected data into multiple segments
with constant-time length. We then perform Singular Value Decomposition (SVD)
of $\mathbf{A}$ for each segment:

$\displaystyle\mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{U}^{\top}$ (33)

where
$\mathbf{S}=\operatorname{diag}(\sigma_{1},\sigma_{2},\cdots,\sigma_{k})$ is
the diagonal matrix of singular values in decreasing order, and
$\mathbf{U}=\begin{bmatrix}\mathbf{u}_{1},\mathbf{u}_{2},\cdots,\mathbf{u}_{k}\end{bmatrix}$
is an orthogonal matrix. As in [50], we employ the information metric (or
observability index) of an individual segment based on the minimum singular
value – that is, a larger one implies a more informative segment. With that,
for a specific sequence of collected data, we first evaluate the information
metric of each segment and then select the most informative ones for LI
calibration.

###  V-B Observability-Aware Update

While we have tried to collect the most informative data for calibration, it
still might be the case that the selected data cannot guarantee observability
and thus some states cannot be estimated. For example, when a ground vehicle
moves on planar surfaces in an urban environment, the lack of pitch and roll
rotation and translation on the vertical direction often leads to part of
states unobservable. As shown in [4], one-axis rotation is a so-called
degenerate motion causing the translation of extrinsics ${}^{I}\mathbf{p}_{L}$
partially unobservable. If this unobservability issue is not properly
addressed, the calibration performance would be largely degraded.

To make the proposed LI calibration robust to (weak) unobservability, we
explicitly enforce observability constraints during each iteration; that is,
only updating the states lying on the observable directions. To this end, we
employ the truncated SVD (TSVD) [51, 52] to perform low-rank approximation of
the information matrix, i.e., selecting the singular values
$(\sigma_{1},\sigma_{2},\cdots,\sigma_{l})$ over an information threshold
$\epsilon$ that is a design choice. With this observability-assurance
information matrix and SVD-based matrix pseudo-inverse, we have the following
solution [see (32)]:

$\displaystyle\delta\hat{\mathbf{x}}=\sum^{l}_{i=1}\frac{\mathbf{u}^{\top}_{i}\mathbf{b}\mathbf{u}_{i}}{\sigma_{i}}.$
(34)

As a result, as compared to solving for $\delta\hat{\mathbf{x}}$ by using all
the singular values of the full information matrix – which may erroneously
introduce observable directions due to numerical issues – we here explicitly
avoid updating the unobservable (or weakly observable) parameters, thus
conservatively ensuring observability.

##  VI Simulation Results

To validate the performance of the proposed continuous-time batch optimization
based LI calibration method, we conduct a series of simulations with ground
truth of all estimated parameters to examine the accuracy of intrinsic and
extrinsic calibration results, and also numerically study the sensitivity of
the system to the time offset in simulations.

Figure 5: The simulated environmental setup filled with planes in different
colors. The colored points indicate a simulated 3D LiDAR scan. Left: The
sinusoidal trajectory for testing calibration in full excitation motion.
Right: Figure-8-shape trajectory for testing calibration in degenerate motion.
Table I: Standard deviations for generating random LiDAR and IMU intrinsics in
simulation experiments.

Parameter Value Parameter Value Elevation angle correction $\delta\phi$ (∘)
0.1 Vertical spatial offset $V$ ($m$) 0.005 Azimuth angle correction
$\delta\theta$ (∘) 0.05 Horizontal spatial offset $H$ ($m$) 0.005 Range scale
factor$s$ (%) 0.005 Range offset $\delta\rho$ ($m$) 0.01 Gyro. Non-linearity
(%) 0.01 Accel. Non-linearity (%) 0.1 Gyro. Non-orthogonality (∘) 0.05 Accel.
Non-orthogonality (∘) 0.05

For the results presented below, we simulate an Xsens-100
IMU111https://www.xsens.com/products/mti-100-series/ and a 16-beam 3D LiDAR
VLP-16222https://velodyneLiDAR.com/vlp-16.html. The synthetic noises in
measurement are consistent with the real-world sensors. IMU measurements are
reported at 400Hz. LiDAR scans are received at 10Hz and the laser rays cover
360∘ horizontally and $\pm$15∘ vertically. The motion distortion effect is
also simulated to match reality. The intrinsic parameters of IMU and LiDAR are
randomly generated from normal distributions with standard deviations provided
in Table I. The extrinsic parameters from LiDAR to IMU are set as
$[0.3,0.15,0.05]^{\top}$ meter in translation and $[1,2,5]^{\top}$ degree in
rotation in all simulation experiments. Note that the extrinsic translation
has a minimum component of 0.05 meters and a maximum component of 0.3 meters
to verify the calibration accuracy of the proposed method for different
magnitudes of the translation.

For the simulated trajectories, a sinusoidal trajectory and a figure-8-shape
trajectory are introduced. The sinusoidal trajectory is formulated as:

$\left\\{\begin{array}[]{l}p_{x}(t)=2.0\cdot\cos(\frac{\pi}{5}\cdot t)+5\\\
p_{y}(t)=1.5\cdot\sin(\frac{\pi}{5}\cdot t)+5\\\
p_{z}(t)=0.8\cdot\cos(\frac{4\pi}{5}\cdot t)+5\\\ r_{x}(t)=0.4\cdot\cos(t)\\\
r_{y}(t)=0.6\cdot\sin(t)\\\ r_{z}(t)=0.7\cdot t\end{array}\right.$ (35)

and the trajectory of figure-8-shape is as follows:

$\left\\{\begin{array}[]{l}p_{x}(t)=2.0\cdot\cos(\frac{\pi}{5}\cdot t)\\\
p_{y}(t)=1.5\cdot\sin(\frac{\pi}{5}\cdot t)\cdot\cos(\frac{\pi}{5}\cdot
t)+5\\\ p_{z}(t)=2\\\ r_{x}(t)=r_{y}(t)=0\\\ r_{z}(t)=0.4\cdot\sin(t)\\\
\end{array}\right.$ (36)

where $t\in[0,10]$, $\mathbf{p}(t)=[p_{x}(t)\ \ p_{y}(t)\ \ p_{z}(t)]^{\top}$
is the translational trajectory in meter and
${\mathbf{r}}_{\text{euler}}(t)=[r_{z}(t)\ \ r_{y}(t)\ \ r_{z}(t)]^{\top}$
represents the euler angles in rad. The sinusoidal trajectory is with fully-
excited motions, while the figure-8-shape trajectory is a simulation of planar
motion, where the robot moves on the $x$-$y$ plane, and orientations change
along $z$-axis only. For both simulated trajectories, it is straightforward to
derive the linear acceleration and angular velocity at any time in order to
generate synthetic IMU readings. Note that in each experiment, we perform 10
Monte Carlo experiments with the same time span of 10 seconds.

The simulated environment is a room constructed with many planes and a volume
of $12\times 10\times 10\ m^{3}$ as shown in Fig. 5, where the two types of
trajectories are also plotted. The sinusoidal trajectory evaluates the
intrinsic and extrinsic calibration in full motion excitation while the
figure-8-shape trajectory for analyzing the extrinsic calibration performance
in degenerate motion.

To evaluate the impact of LiDAR intrinsic calibration on the accuracy of map
reconstruction, we introduce the mean map entropy (MME) metric [53], which
evaluates the sharpness of a map. The entropy $h$ for a map point
$\mathbf{p}_{k}$ is calculated by

$\displaystyle h(\mathbf{p}_{k})=\frac{1}{2}\ln{|2\pi
e\Sigma(\mathbf{p_{k}})|}$ (37)

where $\Sigma(\mathbf{p_{k}})$ is the sample covariance of map points within a
local radius $r$ around $\mathbf{p_{k}}$. Local radius $r$ is 0.3 meters in
our experimental evaluations. The MME $H(Q)$ is averaged over all map points
$Q$:

$\displaystyle H(Q)=\frac{1}{n}\sum_{k=1}^{n}h(\mathbf{p_{k}})\text{.}$ (38)

The smaller the MME value the sharper the map is (i.e., the better map
reconstruction).

Table II: Results of 10 Monte-Carlo simulations with the sinusoidal trajectory: mean and standard deviation of the extrinsic calibration results with and without calibrating the intrinsics of LiDAR and IMU. Note that the best results are marked in bold. Parameter | Gt | W/ Calib. Intr. | W/O Calib. Intr.
---|---|---|---
$p_{x}$ (cm) | 30 | 29.81 |  $\pm$ | 0.54 | 35.15 |  $\pm$ | 0.51
$p_{y}$ (cm) | 15 | 14.86 |  $\pm$ | 0.25 | 22.14 |  $\pm$ | 0.73
$p_{z}$ (cm) | 5 | 4.68 |  $\pm$ | 0.28 | 1.67 |  $\pm$ | 0.14
roll (deg) | 1 | 1.16 |  $\pm$ | 0.02 | 1.97 |  $\pm$ | 0.02
pitch (deg) | 2 | 1.91 |  $\pm$ | 0.02 | 2.41 |  $\pm$ | 0.03
yaw (deg) | 5 | 5.02 |  $\pm$ | 0.05 | 5.44 |  $\pm$ | 0.05
Table III: Results of simulations with the sinusoidal trajectory: time offset
calibrations.

| Time offset (ms)
---|---
True value | 1 | 2 | 3 | 5 | 8 | 12 | 21
Estimated | 1.01 | 1.86 | 2.77 | 4.79 | 7.87 | 12.94 | 20.63
Error | 0.01 | -0.14 | -0.23 | -0.21 | -0.13 | -0.06 | -0.37

Figure 6: Results of 10 Monte-Carlo simulations with the sinusoidal
trajectory: the IMU intrinsic parameters change with the increased number of
iterations. The ground truth is denoted with the red dashed line. Each row
from top to bottom represents scale and misalignment parameters of gyroscope,
and scale and misalignment parameters of accelerometer, and euler angles of
${}^{\omega}_{I}\mathbf{R}$ in (8), respectively.  Figure 7: Results of 10
Monte-Carlo simulations with the sinusoidal trajectory: the intrinsic
calibration results of a simulated 16-beam 3D LiDAR. For the first laser, only
the distance offset and distance scale are estimated while the other
parameters are fixed as a reference during the calibration process. The blue
dots are the ground truth and red lines are the mean of each parameter.

###  VI-A Calibration With Full Motion Excitation

####  VI-A1 Spatial-Temporal Calibration

The extrinsic calibration results with and without calibrating the intrinsic
parameters of LiDAR and IMU in 10 simulated sequences are shown in Table II.
We can find that with the awareness of intrinsic parameters, the proposed
method OA-LICalib achieves an accuracy of 0.40 cm in translation and 0.18
degree in rotation, and the calibration errors extremely increase without
intrinsic calibrations of LiDAR and IMU. It also confirms OA-LICalib’s
capability of calibrating extrinsic parameters with components at different
magnitudes, e.g., translation varies from 5cm to 30 cm, and rotation varies
from 1 degree to 5 degrees in this experiment. It is also interesting to note
that the magnitudes of the rotational and translational trajectories in (35)
affect the accuracy of the extrinsic calibration results. In the $x$-axis
translation and $z$-axis rotation, OA-LICalib gives calibrated parameters with
much lower relative error. This phenomenon agrees with the observability
analysis in [54] that calibration trajectory with strong motion excitation,
i.e., large accelerations and angular velocities, can lead to more accurate
calibrations of extrinsic parameters than weal excitation.

Furthermore, the temporal calibration results in 7 trails are shown in Table
III. Note that time offsets vary from 1 ms to 21 ms, are manually added to one
of the simulated sequences and calibrated. For the precision of temporal
calibration, the proposed method can achieve millisecond-level accuracy.

Table IV: Results of 10 Monte-Carlo simulations with the sinusoidal trajectory: mean and error of IMU intrinsic parameters after 14 iterations. Parameter | Gt | Estimated | Error
---|---|---|---
Sw0 | 0.99738 | 0.99714 | $\pm$ | 0.00020 | -0.00024
Sw1 | 0.99885 | 0.99882 | $\pm$ | 0.00031 | -0.00003
Sw2 | 1.00329 | 1.00169 | $\pm$ | 0.00132 | -0.00160
Sa0 | 0.98487 | 0.98374 | $\pm$ | 0.00071 | -0.00113
Sa1 | 1.01851 | 1.01780 | $\pm$ | 0.00164 | -0.00071
Sa2 | 0.97975 | 0.97936 | $\pm$ | 0.00037 | -0.00039
Mw0 | -0.00127 | -0.00117 | $\pm$ | 0.00015 | 0.00010
Mw1 | -0.00279 | -0.00388 | $\pm$ | 0.00193 | -0.00109
Mw2 | 0.00254 | 0.00577 | $\pm$ | 0.00254 | 0.00323
Ma0 | 0.01607 | 0.01553 | $\pm$ | 0.00129 | -0.00054
Ma1 | 0.02197 | 0.02238 | $\pm$ | 0.00049 | 0.00041
Ma2 | -0.00825 | -0.00825 | $\pm$ | 0.00052 | 0.00001
I.roll (deg) | 0.80000 | 0.81411 | $\pm$ | 0.01924 | 0.01411
I.pitch (deg) | 0.50000 | 0.48720 | $\pm$ | 0.02066 | -0.01280
I.yaw (deg) | 0.30000 | 0.31130 | $\pm$ | 0.03401 | 0.01130
Table V: Results of 10 Monte-Carlo simulations with the sinusoidal trajectory: the averaged absolute errors of the estimated LiDAR intrinsic parameters after 14 iterations.  Parameter | Error | Parameter | Error
---|---|---|---
$s\ (\%)$ | 0.0058 | $\pm$ | 0.0023 | $\delta\rho\ (mm)$ | 0.3048 | $\pm$ | 0.2167
$\delta\phi\ (deg)$ | 0.0075 | $\pm$ | 0.0081 | $V\ (mm)$ | 0.6383 | $\pm$ | 0.7838
$\delta\theta\ (deg)$ | 0.0029 | $\pm$ | 0.0019 | $H\ (mm)$ | 0.2111 | $\pm$ | 0.1106
Figure 8: The MME evaluation of the LiDAR map quality along with the number of
iterations in a typical simulated experiment. The evaluated point cloud map is
initially assembled with the poses from NDT-based LiDAR odometry, then
assembled with iteratively refined poses and intrinsics from the continuous-
time batch optimization. Figure 9: Simulation to examine observability
awareness: three mounting setups with different relative orientations between
the robot and the LiDAR-IMU sensor rig are tested (LiDAR and IMU are rigidly
connected, while LiDAR is not shown in this figure). From left to right: case
A, case B, case C. Among these three cases, the robots move in a shared
trajectory of figure-8 shape as (36), while the unobservable extrinsic
directions in local IMU frames are different.  Table VI: Results of simulation
with the figure-8 trajectory: the ground truth and the automatically-detected
unobservable extrinsic directions under three mounting cases are shown in Fig.
9. The head part corresponds to extrinsic rotation and the tail part
corresponds to extrinsic translation.

Case |  | Degenerate direction
---|---|---
| Truth | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 1.00000
Case A | Detected | -0.00001 | 0.00000 | 0.00005 | 0.00165 | -0.00110 | 1.00000
| Truth | 0.00000 | 0.00000 | 0.00000 | 0.50000 | 0.00000 | 0.86603
Case B | Detected | 0.00000 | 0.00000 | -0.00001 | 0.50117 | -0.00023 | 0.86535
| Truth | 0.00000 | 0.00000 | 0.00000 | 0.50000 | 0.43301 | 0.75000
Case C | Detected | 0.00006 | 0.00004 | 0.00007 | 0.50070 | 0.43296 | 0.74957

Table VII: Results of 10 Monte-Carlo simulations with the figure-8 trajectory:
mean of the estimated extrinsic parameters in three mounting cases with or
without observability-aware calibration.

|  |  ${}^{I}\mathbf{p}_{L}$ (cm) |  ${}^{I}_{L}\bar{q}$ in Euler (deg) | RMSE
---|---|---|---|---
Seq. | Method | $p_{x}$ | $p_{y}$ | $p_{z}$ | roll | pitch | yaw | Pos (cm) | Rot (deg)
Groundtruth | 30.00 | 15.00 | 5.00 | 1.00 | 2.00 | 5.00 |  |
| W/O Obser. | 27.24 | 11.49 | 7.17 | 0.85 | 2.40 | 4.83 | 2.87 | 0.27
Case A | W Obser. | 30.97 | 12.09 | 8.05 | 0.95 | 2.44 | 4.89 | 2.50 | 0.26
| W/O Obser. | 26.50 | 11.53 | 6.88 | 0.76 | 2.27 | 4.88 | 3.04 | 0.22
Case B | W Obser. | 29.73 | 11.94 | 5.69 | 1.01 | 2.27 | 4.81 | 1.82 | 0.19
| W/O Obser. | 26.32 | 11.20 | 6.03 | 0.08 | 2.28 | 4.86 | 3.11 | 0.56
Case C | W Obser. | 29.67 | 10.42 | 5.56 | 0.18 | 2.28 | 4.88 | 2.67 | 0.50

####  VI-A2 Intrinsic Calibration

Figure 6 illustrates the convergence process of IMU intrinsic parameters along
with iterations and visualizes the distribution of calibration results for 10
simulated trials. Note that I.roll, I.pitch and I.yaw are the corresponding
Euler angles of ${}^{\omega}_{I}\mathbf{R}$ in (8). The intrinsic parameters
remain unchanged during the first iteration and start to be optimized in the
second iteration, and the results gradually stabilize after about 8
iterations. Since the 11th iteration, the raw measurements are corrected based
on the estimated intrinsic parameters (as mentioned in Sec. IV-C). We can see
that the estimated intrinsic parameters finally converge to near the true
value, and the final estimated intrinsic results are summarized in Table IV.
The experiment indicates that the proposed method has high accuracy and
repeatability in intrinsic calibration.

The final calibration results of the simulated 16-beam LiDAR intrinsic
parameters in 10 trails are shown in Fig. 7. Besides, the averaged absolute
errors of the estimated LiDAR intrinsic parameters in 6 types are summarized
in Table V. The proposed method can identify six types of intrinsic parameters
of multi-beam LiDAR and estimate accurate results. We further evaluate the
quality of LiDAR point cloud map which is affected by both the estimated
trajectory and the calibrated parameters. Figure 8 shows the MME metric in a
typical simulated experiment. The evaluated point clouds consisting of around
$215k$ points are assembled by the poses from NDT-based LiDAR odometry
firstly, and refined poses from the continuous-time batch optimization. The
map becomes more accurate with increased iterations due to the batch
optimization over the continuous-time trajectory and calibrated parameters.
Notably, a significant improvement in map accuracy after the 11th iteration
due to the correction of the raw LiDAR points based on the current best
estimation of LiDAR intrinsics. It is obvious that intrinsic parameters
calibration significantly affects both the LiDAR map quality and the extrinsic
calibration. Therefore, we advocate calibrating the intrinsic parameters to
enable highly accurate LIDAR localization and mapping.

###  VI-B Calibration Under Degenerate Motion

In the same simulation environment with figure-8-shape trajectory, shown in
the right picture of Fig. 5, we investigate the extrinsic calibration accuracy
of OA-LICalib in the degenerate case. With the trajectory in (36), where
rotation happens along the $z$-axis only, the LiDAR-IMU extrinsic translation
along the rotation axis is unobservable according to the analysis in [4]. In
order to examine whether the proposed OA-LICalib is able to automatically
identify the unobservable directions, we also simulate three different
mounting cases (A, B, C) as shown in Fig. 9, where the LiDAR-IMU sensor rig is
rigidly attached to the robot platform at three different relative
orientations. Among the three cases, the robots move in the same trajectory
(as Eq. (36)), thus the unobservable directions of the extrinsic translation
are the identical vertical axis in the robot coordinate frame. However, due to
the different relative rotations between robot platform and the LiDAR-IMU
sensor rigs, the unobservable directions of the extrinsic translation,
${}^{I}\mathbf{p}_{L}$, are different when represented in local IMU frame.

Table VI summarizes the unobservable direction of the three mounting cases in
Fig. 9. For case A where the IMU $z$-axis is aligned with the vertical axis,
the unobservable direction is purely in the $z$-axis of extrinsic translation,
${}^{I}\mathbf{p}_{L}$. For a typical trial, OA-LICalib identifies the
singular value of the information matrix corresponding to extrinsics as
$[41760,\ 22193,\ 3525,\ 491,\ 295,\ 0.000]^{\top}$, and it’s obvious that
there exists one unobservable direction. The eigen vectors corresponding to
the minimum singular value (deemed as the detected unobservable direction) in
the three cases are also summarized in Table VI. We can see that the detected
unobservable directions are consistent with the truths. To verify the
effectiveness of the observability scheme in Section. V-B, we conduct the
ablation experiments on extrinsic calibration by starting from perturbed
initial extrinsic parameters (with perturbations of around 3 degrees and 3
centimeters on rotations and translations, respectively). Table VII shows the
extrinsic calibration results
($\\{{}^{I}_{L}\bar{q},\,{}^{I}\mathbf{p}_{L}\\}$) of 10 Monte-Carlo trials,
and the Root Mean Square Error (RMSE) is calculated by comparing the mean
calibration results to the ground truths. OA-LICalib with observability
awareness keeps the unobservable directions of extrinsic translation around
the initial values, and outperforms the method without observability
awareness.

##  VII Real-world Experiments

Table VIII: The details of sequences used in real-world experiments.

Dataset | Sequence | Config | Duration | Env. |  | Ave. Angular
---
Velocity (rad/s)
Motion Type | Excitation Type
|  | Vicon-
---
(01,02,…,10)
Hand-held | 15 | Indoor | 0.983 | Random Shaking | Full
|  | Stairs-
---
(01,02)
Hand-held | 60 | Indoor | 1.118 | Random-Shaking | Full
|  | Lab-
---
(01,02)
Hand-held | 60 | Indoor | 1.215 | Random Shaking | Full
| SKG-01 | On-Vehicle | 120 | Outdoor | 0.407 | Planar Motion | Degenerate
Self collection | YQ-01 | Hand held | 271 | Outdoor | 0.760 |  | Casual Walking &
---
occasional Random Shaking
Occasionally full
LIO-SAM [5] | Casual-walk | Hand held | 655 | Outdoor | 0.521 |  | Casual Walking &
---
occasional Random Shaking
Occasionally full
| IC_Office | Hand held | 200 | Indoor | 0.348 | Casual Walking | Occasionally full
| Office_Mitte | Hand held | 264 | Indoor | 0.350 | Casual-Walking | Occasionally full
Hilti SLAM challenge [55] | Basement1 | On-Vehicle | 330 | Indoor | 0.116 | Planar Motion | Degenerate
1 Corresponding to Basement_3 sequence in Hilti dataset. |

In the real-world experiments, we validate the proposed method on multiple our
own datasets as well as two public datasets. Table VIII summarizes the
sequences in terms of collection configuration, duration, environment, motion
type, etc. For the motion type:

  * •

Random Shaking: We hold the sensor by hand and shake the sensor suite
aggressively both in rotations and translations;

  * •

Casual Waking: We hold the sensor and walk around;

  * •

Planar Motion: The sensor rig is mounted on moving wheeled robots, see Fig. 10
(right).

The self-collected datasets are collected in both indoor and outdoor scenarios
by the self-assembled sensor rig shown in the left of Fig. 10, which consists
of 3 Xsens-100 IMUs sampled at 400 Hz and a Velodyne VLP-16 LiDAR sampled at
10 Hz. The Casual-walk sequence 333Available at http://www.udel.edu/009352 is
from LIO-SAM [5] and is collected by a Microstrain 3DM-GX5-25 IMU sampled at
500 Hz and a Velodyne VLP-16 LiDAR sampled at 10 Hz. The Hilti SLAM challenge
dataset [55] is intended for assessing the accuracy of state estimation in
challenging environments and it contains a number of visual, LiDAR and
inertial sensors. In the experiments, the data from Bosch BMI085 IMU sampled
at 200 Hz and the Ouster OS0-64 LiDAR sampled at 10 Hz are utilized for
evaluation.

The proposed method is extensively tested on the Stairs and Lab sequences to
examine the accuracy and repeatability of calibrated results. Due to the
absence of ground-truth extrinsic transformation between LiDAR and IMU in
real-world experiments, the relative poses of three components IMU inferred
from CAD assembly drawings are introduced as references. Additionally, the
proposed observability-aware calibration method incorporating the informative
data selection and observability-aware state update mechanism is thoroughly
evaluated. On the one hand, we examine the feasibility of the proposed method
in degenerate motion cases over the SKG-01 sequence and Basement sequence,
which are collected on the vehicles under planar motion. On the other hand, we
also analyze the proposed informative segment selection algorithm in the
YQ-01, Casual-walk, IC_Office, Office_Mitte sequences, among which the former
two sequences are outdoors while the latter two are indoors. Note that only
one of the three IMUs (Fig. 10), IMU1, is utilized to evaluate the SKG-01 and
YQ-01 sequences.

We also compare the proposed method with the open-source LiDAR-IMU calibration
toolkit Imu_Lidar_Calibration [22] (abbreviated as ILC) and our previous
method LI-Calib [16]. Notably, LI-Calib [16] does not have the observability-
aware scheme described in Sec. V or the explicit initialization of control
points (Sec. IV-A1), compared to the proposed OA-LICalib. In addition, two
online calibration methods, LIO-Mapping [3] (abbreviated as LIOM) and FAST-
LIO2 [26], are also included for the comparisons. LIOM and FAST-LIO2 are
LiDAR-IMU odometry and simultaneously calibrate the extrinsic parameters
between LiDAR and IMU. In practice, online calibration methods estimate the
whole trajectory of the entire calibration sequence, and the final converged
values of the extrinsic parameters are deemed the calibration results, as are
the ILC method. Since the ILC method particularly requires the sensor to keep
static for seconds to initialize the system, we skip some data segments at the
beginning in order to satisfy the static initialization requirements in some
sequences. Notably, under general motion scenarios, the proposed OA-LICalib
and our previous method LI-Calib are initialized from the identity extrinsic
rotation matrix and zero extrinsic translations without bells and whistles.
However, for the initial values of the extrinsics for these compared methods,
including ILC, LIOM, and FAST-LIO2, the extrinsic rotation parameters are
initialized by the reference parameters provided by the dataset, while the
translation extrinsic parameters are initialized as zero. These compared
methods fail to converge when starting from the identity extrinsic rotation
matrix due to the lack of an extrinsic initialization module.

Figure 10: Left: The self-assembled sensor rig used in real-world experiments.
Right: The ground wheeled robot for collecting SKG-01 sequence.

Figure 11: Top: Estimated continuous-time trajectories aligned with the
ground-truth. The color of the trajectory indicates ATE. Bottom: The
corresponding fitting results of the top-left trajectory. Only the $x$-axis is
plotted. The fitting errors almost have the same order of magnitude as the IMU
measurement noises.

###  VII-A Trajectory Representation Capability

In the Vicon sequences, ground truth poses are provided by the Vicon system at
100 Hz, with a duration of about 15 seconds to evaluate the accuracy of
estimated trajectories and the representation capability of the splines.
Comparing the estimated trajectories with the motion-capture system’s reported
trajectories, the average absolute trajectory error (ATE) [56] overall ten
sequences are 0.0183 m. Figure 11 shows two typical estimated trajectories
aligned with the ground truth and also illustrates fitting results. The
estimated B-Spline trajectory fits well with the ground-truth trajectory and
the acceleration measurements and angular velocity measurements, indicating
the B-Spline’s high representation capability.

Figure 12: Extrinisc translation calibration results from OA-LICalib over Stairs and Lab sequences. Top: Top view of the calibration results. The three black rectangular boxes represent the actual mounting positions of three IMUs, respectively. Bottom: Front view of the calibration results. Figure 13: Extrinsic Rotation calibration results from OA-LICalib over Stairs and Lab sequences. Table IX: The mean and standard deviation of the calibration results from different methods on Stairs and Lab sequences. The units for translation and rotation are in ‘cm’ and ‘degree’, respectively. The non-converged results of ILC and LIOM are excluded from the statistics.  |  |  ${}^{L}\mathbf{p}_{I}$ (cm) |  ${}^{L}_{I}\bar{q}$ in Euler (deg) | RMSE
---|---|---|---|---
Method | IMU | $p_{x}$ | $p_{y}$ | $p_{z}$ | roll | pitch | yaw | Pos (cm) | Rot (deg)
| IMU1 | 0 | -10 | -13 | 0 | 0 | -90 | - | -
| IMU2 | 10 | 0 | -13 | 0 | 0 | -90 | - | -
CAD Reference | IMU3 | -7 | 8 | -13 | 0 | 0 | 135 | - | -
| IMU1 | 9.61 | $\pm$ | 0.79 | -13.00 | $\pm$ | 4.78 | -7.34 | $\pm$ | 6.75 | 7.78 | $\pm$ | 0.47 | -4.95 | $\pm$ | 5.83 | -95.82 | $\pm$ | 2.18 | 11.65 | 5.04
| IMU2 | 19.65 | $\pm$ | 5.63 | -1.69 | $\pm$ | 3.36 | -4.44 | $\pm$ | 3.54 | 9.73 | $\pm$ | 1.40 | -3.48 | $\pm$ | 3.35 | -93.59 | $\pm$ | 1.59 | 6.67 | 6.30
ILC | IMU3 | 7.74 | $\pm$ | 2.13 | 7.20 | $\pm$ | 2.33 | -10.95 | $\pm$ | 0.49 | -6.63 | $\pm$ | 2.04 | 7.69 | $\pm$ | 2.05 | 127.47 | $\pm$ | 1.35 | 7.51 | 6.32
| IMU1 | -0.74 | $\pm$ | 3.77 | -12.87 | $\pm$ | 2.40 | -2.91 | $\pm$ | 4.31 | 9.64 | $\pm$ | 1.19 | 0.09 | $\pm$ | 0.25 | -93.31 | $\pm$ | 2.37 | 8.60 | 7.30
| IMU2 | 7.70 | $\pm$ | 6.42 | -5.85 | $\pm$ | 2.63 | -2.68 | $\pm$ | 2.71 | 9.82 | $\pm$ | 1.71 | 0.19 | $\pm$ | 0.29 | -94.06 | $\pm$ | 2.50 | 6.98 | 6.13
LIOM | IMU3 | -2.64 | $\pm$ | 2.54 | 11.21 | $\pm$ | 7.15 | 6.43 | $\pm$ | 3.96 | 7.89 | $\pm$ | 4.09 | 0.32 | $\pm$ | 0.27 | 131.30 | $\pm$ | 2.97 | 11.65 | 5.04
| IMU1 | -1.56 | $\pm$ | 1.28 | -10.56 | $\pm$ | 3.24 | -12.94 | $\pm$ | 2.78 | -0.52 | $\pm$ | 1.65 | 1.49 | $\pm$ | 1.06 | -91.00 | $\pm$ | 0.92 | 0.96 | 1.08
| IMU2 | 6.76 | $\pm$ | 1.89 | -3.08 | $\pm$ | 2.40 | -12.39 | $\pm$ | 1.96 | -1.05 | $\pm$ | 2.00 | 1.45 | $\pm$ | 0.89 | -91.34 | $\pm$ | 1.32 | 2.61 | 1.29
FAST-LIO2 | IMU3 | -7.84 | $\pm$ | 1.58 | 4.48 | $\pm$ | 2.72 | -10.41 | $\pm$ | 0.85 | -0.34 | $\pm$ | 1.03 | -0.7 | $\pm$ | 2.13 | 132.96 | $\pm$ | 0.78 | 2.57 | 1.26
| IMU1 | -0.97 | $\pm$ | 2.37 | -9.80 | $\pm$ | 1.49 | -12.98 | $\pm$ | 1.41 | 0.47 | $\pm$ | 0.40 | 0.18 | $\pm$ | 0.87 | -90.90 | $\pm$ | 0.29 | 0.57 | 0.60
| IMU2 | 9.64 | $\pm$ | 2.09 | -1.08 | $\pm$ | 1.10 | -12.47 | $\pm$ | 1.25 | 0.36 | $\pm$ | 0.39 | 0.42 | $\pm$ | 0.83 | -91.56 | $\pm$ | 0.27 | 0.73 | 0.96
LI-Calib | IMU3 | -6.73 | $\pm$ | 2.03 | 7.09 | $\pm$ | 1.33 | -12.66 | $\pm$ | 1.11 | -0.56 | $\pm$ | 0.60 | 0.89 | $\pm$ | 0.81 | 133.14 | $\pm$ | 0.31 | 0.58 | 1.23
| IMU1 | -0.22 | $\pm$ | 0.50 | -9.92 | $\pm$ | 0.31 | -13.48 | $\pm$ | 0.44 | 0.23 | $\pm$ | 0.04 | 0.21 | $\pm$ | 0.06 | -90.89 | $\pm$ | 0.07 | 0.31 | 0.54
| IMU2 | 10.16 | $\pm$ | 0.36 | -1.00 | $\pm$ | 0.31 | -12.83 | $\pm$ | 0.54 | 0.16 | $\pm$ | 0.07 | 0.46 | $\pm$ | 0.17 | -91.53 | $\pm$ | 0.05 | 0.59 | 0.93
OA-LICalib | IMU3 | -7.19 | $\pm$ | 0.53 | 7.46 | $\pm$ | 0.34 | -12.89 | $\pm$ | 0.61 | -0.43 | $\pm$ | 0.04 | 0.68 | $\pm$ | 0.05 | 133.14 | $\pm$ | 0.09 | 0.33 | 1.17
Table X: The OA-LICalib’s time consumption of main stages. Two segments from
Lab-01 and SKG-01 sequence are evaluated, respectively. The former segment is
under fully-excited motion, while the latter is under degeneration motion.

Step |  | Fully-excited Seq.
---
[0,15] of Lab-01
| Degenerate Seq.
---
[0,15] of SKG-01
Count |  | Processing
---
/Step(s)
Total(s) | Count |  | Processing
---
/Step(s)
Total(s)
Extri. Init. & Rerun Odom. | 1 | 15.73 | 15.73 | 1 | 12.25 | 12.25
Data Association | 1 | 1.80 | 1.80 | 1 | 3.31 | 3.31
Init. Traj. & First Opt. | 1 | 4.47 | 4.47 | 1 | 6.54 | 6.54
Refinement | 5 | 6.73 | 33.64 | 5 | 8.03 | 40.17
TSVD | 0 | - | - | 5 | 0.03 | 0.13
Total | - | - | 55.65 | - | - | 64.40

###  VII-B Calibration with Full Motion Excitation

We gathered four sequences of a 60-second duration with full excitation in
structured environments, Stairs-01 and Stairs-02 collected at the stairway,
Lab-01 and Lab-02 in the laboratory, verifying the calibration accuracy of the
full motion excitation data. Each sequence is split into four segments of 15
seconds and calibrated independently for calibration test. The extrinsic
calibration results of three IMUs with respect to the LiDAR by the proposed
OA-LICalib are shown in Fig. 12 and Fig. 13, where the black rectangles
indicate the CAD sketch of three IMUs. We can intuitively find that the
translational calibration results are distributed within rectangles, with the
distribution range in translation at around 2cm. The distribution range of
rotation concentrates within 0.2 degrees for each IMU. These results
demonstrate the proposed method’s high repeatability and reliability in
practical deployments.

Additionally, the quantitative comparison results are summarized in Table IX,
where RMSE is calculated by comparing it to the CAD reference. Note that LIOM
fails to converge on the Stairs sequences and ILC does not converge on the Lab
sequences in our experiments. Failing to convergence means rotational RMSE is
greater than 30 degrees. The non-converged results are excluded from the
statistics in Table IX. Although the discrete-time methods, ILC, LIOM and
FAST-LIO are initialized from the reference extrinsic rotation, the
continuous-time based methods LI-Calib and OA-LICalib show apparent advantages
over all of them, while starting from the identity extrinsic rotation matrix.
In addition, compared to LI-Calib, OA-LICalib takes more effort on the
initialization of control points (Sec. IV-A1), which benefits the convergence
and improves the accuracy of calibrations. It is worthwhile to mention that
the calibration results of OA-LICalib are with high calibration accuracy yet
low standard deviations, which demonstrate its great reliability. To give an
impression of the efficiency, the proposed OA-LICalib’s time consumption of
main stages in two typical 15s segment over the sequence Lab-01 (under fully-
excited motion) and sequence SKG-01 (under degenerate motion) are summarized
in Table X. Basically, the proposed method takes about only one minute to
generate the highly-accurate calibration results automatically on a 15s
calibration segment.

Figure 14: Evaluation in YQ-01 sequence. Top: The informativeness of every
individual segment of YQ-01 sequence, which has a total trajectory length of
363 meters; the corresponding three-axis angular velocity curves. Bottom: The
whole trajectory is plotted in dashed line and the selected segments in bright
colors. The trajectory direction is counterclockwise.  Table XI: The estimated
spatial parameters $\\{{}^{I}_{L}\bar{q},{}^{I}\mathbf{p}_{L}\\}$ and time
offset $t_{c}$ using one segment (Segment-$n$) or multi-segment ($n$ Segments)
over YQ-01 sequence.

|  ${}^{I}\mathbf{p}_{L}$ (cm) |  ${}^{I}_{L}\bar{q}$ in Euler (deg) |  | RMSE
---|---|---|---|---
Data | $p_{x}$ | $p_{y}$ | $p_{z}$ | roll | pitch | yaw |  Time Offset (ms) | Pos (cm) | Rot (deg)
CAD Reference | -10 | 0 | 13 | 180 | -180 | -90 | - |  |
Segment-1 | -11.0 | 0.3 | 14.6 | 179.75 | -179.76 | -89.25 | 7.4 | 1.10 | 0.48
Segment-2 | -11.0 | -0.2 | 14.0 | 179.77 | -179.81 | -89.23 | 7.6 | 0.82 | 0.48
Segment-3 | -11.2 | 0.7 | 13.5 | 179.77 | -179.8 | -89.18 | 7.4 | 0.85 | 0.51
Segment-4 | -11.0 | 0.0 | 14.7 | 179.73 | -179.77 | -89.2 | 7.8 | 1.14 | 0.51
Segment-5 | -9.7 | 0.3 | 12.8 | 179.71 | -179.81 | -89.18 | 7.4 | 0.27 | 0.51
Segment-6 | -10.6 | -0.2 | 15.1 | 179.78 | -179.82 | -89.19 | 7.9 | 1.27 | 0.50
2 Segments | -11.3 | 0.2 | 14.4 | 179.76 | -179.77 | -89.25 | - | 1.11 | 0.47
3 Segments | -11.3 | 0.3 | 14.2 | 179.76 | -179.77 | -89.24 | - | 1.04 | 0.48
4 Segments | -11.2 | 0.2 | 14.3 | 179.75 | -179.77 | -89.23 | - | 1.03 | 0.49
5 Segments | -11.1 | 0.2 | 14.1 | 179.75 | -179.79 | -89.23 | - | 0.91 | 0.48
6 Segments | -11.1 | 0.2 | 14.3 | 179.76 | -179.8 | -89.23 | - | 0.99 | 0.48

Figure 15: Top: The informativeness of every individual segment of IC-Office
sequence. Bottom: the corresponding linear acceleration curves and angular
velocity curves. Table XII: Comparison results of different methods on Hilti
dataset. LIOM fails to give reasonable results, thus not listed in the table.

|  |  | ${}^{I}\mathbf{p}_{L}$ (cm) | ${}^{I}_{L}\bar{q}$ in Euler (deg) |  | RMSE
---|---|---|---|---|---|---
Method | Sequence | Data | $p_{x}$ | $p_{y}$ | $p_{z}$ | roll | pitch | yaw |  Time Offset (ms) | Pos (cm) | Rot (deg)
Reference |  | 1.00 | -0.66 | 9.47 | 0.02 | -0.37 | -0.14 |  $\sim$1ms |  |
| Office_Mitte | All data | 0.43 | -1.02 | 6.58 | 0.05 | -0.02 | -0.35 | - | 1.71 | 0.24
FAST-LIO2 | IC_Office | All data | 0.69 | -1.21 | 5.02 | 0.03 | 0.12 | -0.59 | - | 2.59 | 0.39
| Office_Mitte | All data | -0.02 | -1.09 | 11.95 | 1.54 | -0.50 | 0.12 | -3.86 | 1.57 | 0.89
ILC | IC_Office | All data | -0.24 | 3.65 | 10.95 | -0.01 | 0.97 | 5.16 | -11.90 | 2.73 | 3.15
|  | Segment-1 | 0.20 | -0.70 | 9.10 | 0.12 | -0.07 | -0.26 | 1.62 | 0.51 | 0.20
|  | Segment-2 | 0.71 | 0.86 | 5.57 | 0.25 | -0.31 | 0.45 | 1.76 | 2.42 | 0.37
|  | Segment-3 | -0.44 | -0.66 | 6.96 | 0.04 | -0.20 | -0.08 | 1.65 | 1.67 | 0.10
|  | Segment-4 | -2.64 | -1.74 | 9.73 | 0.00 | 0.10 | 0.09 | 1.53 | 2.20 | 0.30
|  | 2-Segment | -0.15 | 0.16 | 7.21 | 0.16 | -0.14 | -0.11 | - | 1.54 | 0.16
|  | 3-Segment | -0.22 | 0.25 | 7.72 | 0.12 | -0.15 | -0.15 | - | 1.34 | 0.14
| Office_Mitte | 4-Segment | -0.20 | 0.02 | 7.69 | 0.11 | -0.13 | -0.18 | - | 1.30 | 0.15
|  | Segment-1 | 1.16 | -2.30 | 6.33 | -0.05 | -0.20 | -0.43 | 1.55 | 2.05 | 0.20
|  | Segment-2 | -1.49 | 0.58 | 8.81 | -0.03 | -0.42 | -0.02 | 0.96 | 1.65 | 0.08
|  | Segment-3 | 0.77 | 2.17 | 5.67 | 0.10 | -0.31 | -0.09 | 1.30 | 2.74 | 0.06
|  | Segment-4 | -0.17 | 0.17 | 7.34 | -0.20 | -0.25 | -0.07 | 1.43 | 1.48 | 0.15
|  | 2-Segment | 0.48 | -0.86 | 7.37 | -0.03 | -0.27 | -0.45 | - | 1.25 | 0.19
|  | 3-Segment | 0.64 | 0.24 | 6.76 | 0.01 | -0.27 | -0.27 | - | 1.66 | 0.09
OA-LICalib | IC_Office | 4-Segment | 0.53 | 0.30 | 6.91 | -0.02 | -0.27 | -0.23 | - | 1.60 | 0.08

Table XIII: Multiple-segment calibration using non-informative segment in IC-
Office sequences. The segments [60,75] and [105,120] are identified as
informative segments and segment [15,30] is identified as non-informative
segment, as shown in Fig. 15.

| ${}^{I}\mathbf{p}_{L}$ (cm) | ${}^{I}_{L}\bar{q}$ in Euler (deg) |  | RMSE
---|---|---|---|---
Data | $p_{x}$ | $p_{y}$ | $p_{z}$ | roll | pitch | yaw |  Time Offset (ms) | Pos (cm) | Rot (deg)
Reference | 1.00 | -0.66 | 9.47 | 0.02 | -0.37 | -0.14 |  $\sim$1ms |  |
[60,75] | 1.16 | -2.30 | 6.33 | -0.05 | -0.20 | -0.43 | 1.55 | 2.05 | 0.20
[105,120] | -1.49 | 0.58 | 8.81 | -0.03 | -0.42 | -0.02 | 0.96 | 1.65 | 0.08
[15,30] | 0.68 | -2.84 | 1.47 | 0.09 | -0.45 | 0.01 | 0.66 | 4.79 | 0.11
[60,75]+[105,120] | 0.48 | -0.86 | 7.37 | -0.03 | -0.27 | -0.45 | - | 1.25 | 0.19
[60,75]+[15,30] | 0.69 | -2.75 | 3.91 | -0.01 | -0.33 | -0.42 | - | 3.43 | 0.16

###  VII-C Calibration with Informative Segments

There is a relatively common situation in practice where some datasets provide
a large number of various data sequences yet without accurate extrinsic
parameters. There’s no possibility of collecting new sequences with the
sensors used in the existing datasets. Hence, if users want to get more
accurate extrinsic parameters, they need to complete the calibration based on
the existing data sequences. Special care should be taken when identifying
informative sequences for accurate calibration.

Regarding the implementation details about multi-segment calibration, we
firstly perform single segment calibration and identify informative segments;
afterward, with all the measurements from selected informative segments, we
conduct a total continuous-time batch optimization to simultaneously estimate
each segment’s individual states (including the continuous-time trajectory,
time offset $t_{c}$, IMU navigation states $\mathbf{x}_{Is}$), and the shared
spatial extrinsic parameters $\\{{}^{I}_{L}\bar{q},{}^{I}\mathbf{p}_{L}\\}$.
We evaluate the proposed informativeness-aware segment selection algorithm in
the self-collected sequence YQ-01 and public-available sequence Casual-walk
and Office_Mitte, IC_Office from the Hilti dataset. Like the Casual-walk
sequence collection process, we walk in a circle profile with the occasional
shake of the sensor to collect YQ-01 sequence. Again, we split sequences into
segments with a duration of 15-second for each segment.

####  VII-C1 YQ-01 Sequence

Figure 14 details the shape of the trajectory and the corresponding angular
velocity curves of the YQ-01 sequence. The top figure in Fig. 14 illustrates
each segment’s information metrics and angular velocities. The segments with
minimal singular value (see Sec. V-A) above 100 (black dash line) are selected
as the informative segments and the relative location of these segments are
shown in bright colors at the bottom of Fig. 14. It is interesting to see the
amplitudes of angular velocity can indistinctly reflect the information
metrics. The extrinsic calibration results with data of one segment
(Segment-$n$) and multi-segment ($n$ Segments) are summarized in Table XI.
Note that “$n$ Segments” means we take the first $n$ informative segments,
e.g., “3 Segments” represents Segment-1, Segment-2 and Segment-3 are jointly
utilized for calibration. The evaluation results show that the selected
single-segment calibration is able to generate reasonable calibration results,
which demonstrate the effectiveness of the informative segment selection.
Besides, multi-segment calibration is generally stable compared to single-
segment calibration, and four informative segments are enough to get a
reliable estimation of the extrinsics.

####  VII-C2 Hilti Sequences

Furthermore, we test the multi-segment calibration on two handheld sequences
of the Hilti dataset, the comparison results are reported in Table XII. Figure
15 reports the information metrics of every individual segment, linear
acceleration curves, and angular velocity curves in the IC_Office sequence,
and we pick out the four most informative segments for multi-segment
calibration. The results in Table XII suggest multi-segment calibration is
more stable and more accurate than single-segment calibration. The proposed
OA-LICalib with multi-segment calibration also significantly outperforms ILC
and FAST-LIO2.

Note that using minimal singular value to measure the informativeness of each
individual segment within a single sequence is feasible, however it is not
meaningful to compare the minimal singular values of different sequences under
significantly different scenarios. Minimal singular value is not a general
information metric [9], and it is affected by the sensory observations,
trajectory profile, surrounding environment, etc. This is reflected by the
different amplitudes of minimal singular values over different sequences, as
shown in Fig. 14 and Fig. 15. On the positive side, a minimal-singular-value-
based information metric is sufficient for informative segments selection from
a given sequence and leads to reasonable calibration results.

We further investigate the performance of non-informative segments and their
effect on multi-segment calibration results. Table XIII shows the calibration
results using multiple segments in IC-Office sequences. The segments [60, 75]
and [105, 120] are identified as informative segments, while the segment [15,
30] is a non-informative segment automatically by OA-LICalib, as shown in Fig.
15. The results in Table XIII suggest that the non-informative segment
generates a much bigger calibration error. When the non-informative segments
are included for multi-segment calibration, the calibration results may be
degraded compared to calibration using only informative segments. The multi-
segment calibration result is more like a weighted average of the calibration
results from their component segments, thus the calibration result is more
stable than single-segment calibration.

####  VII-C3 Casual-walk Sequence

The Casual-walk sequence is one of the common cases where the dataset is
publicly available but without providing carefully calibrated extrinsic
parameters. The rough extrinsic parameters are $[0,0,0]^{\top}$ in translation
and identity matrix in rotation. We try to conduct multi-segment calibration
on this sequence. This sequence was collected by randomly walking around the
campus and occasionally shaking the sensor vigorously, meaning that the
sequence is coupled with a number of segments suitable for calibration. We
perform multi-segment data calibration using only half of the sequence which
is divided into 20-second segments, and the results are shown in Table XIV.
Since there are no reliable reference values for the dataset’s extrinsics or
ground truth trajectory, we can only evaluate the quality of the LiDAR point
cloud map, which implicitly reflects the calibration accuracy. The MMEs of the
maps are -2.141, -2.163, and -2.185 for the rough extrinsics, the calibrated
spatial extrinsics, and the calibrated spatial-temporal extrinsics,
respectively. The map quality obtained from the calibrated spatial-temporal
extrinsics is the best, which verifies the effectiveness and accuracy of the
multi-segment calibration.

Table XIV: Evaluation of calibration results using single or multi informative
segments over the Casual-walk sequence.

|  ${}^{I}\mathbf{p}_{L}$ (cm) |  ${}^{I}_{L}\bar{q}$ in Euler (deg) |
---|---|---|---
Data | $p_{x}$ | $p_{y}$ | $p_{z}$ | roll | pitch | yaw |  Time Offset (ms)
Segment-1 | -0.9 | 0.4 | -4.7 | 178.88 | 179.83 | 178.52 | 8.3
Segment-2 | -0.8 | 0.7 | -2.9 | 178.92 | 179.81 | 178.53 | 8.1
Segment-3 | -1.6 | -0.1 | -3.2 | 178.98 | 179.90 | 178.50 | 9.1
Segment-4 | -0.2 | 0.0 | -1.2 | 178.91 | 179.77 | 178.72 | 8.5
2 Segments | -1.0 | 0.5 | -3.8 | 178.89 | 179.81 | 178.52 | -
3 Segments | -1.3 | 0.4 | -3.7 | 178.91 | 179.83 | 178.52 | -
4 Segments | -1.1 | 0.7 | -3.3 | 178.91 | 179.83 | 178.54 | -

Figure 16: The details about the SKG-01 sequence. Top: The environment for collecting SKG-01 sequence and the actual figure-8-shape trajectory. Middle: Concrete trajectories for each 15-second segment, the red stars indicate the start and end position of the trajectory, respectively. Bottom: Corresponding three-axis angular velocity curves. Time | Singular values in a decreasing order
---|---
[0,15] | 5545 | 4659 | 3664 | 958 | 282 | 0.010
[15,30] | 33595 | 16045 | 6844 | 930 | 724 | 0.016
[30,45] | 27126 | 10997 | 7585 | 1915 | 678 | 0.016
[45,60] | 47269 | 25342 | 5606 | 726 | 365 | 0.009
[60,75] | 9332 | 9186 | 3634 | 929 | 452 | 0.008
[75,90] | 60975 | 15907 | 5675 | 563 | 234 | 0.004
[90,105] | 18542 | 7932 | 760 | 290 | 125 | 0.006
[105,120] | 20997 | 10650 | 3940 | 811 | 8 | 0.012
Table XV: The singular value of each segment in SKG-01 sequence identified by
OA-LICalib. The secondary minimum singular value of the segments [90,105] and
[105,120] are relatively small.  Figure 17: Comparison of the estimated
${}^{I}\mathbf{p}_{L}$ under degenerate case over SKG-01 sequence: using OA-
LICalib with and without observability-aware calibration. References are from
the CAD sketch. Figure 18: Comparison of the estimated ${}^{I}\mathbf{p}_{L}$
under degenerate case over Basement sequence: using OA-LICalib with and
without observability awareness. The beginning of the sequence is skipped due
to less motion.  Table XVI: Comparison of the calibration results under
degenerate case over SKG-01 and Basement sequences.

|  |  |  ${}^{I}\mathbf{p}_{L}$ (cm) |  ${}^{I}_{L}\bar{q}$ in Euler (deg) | RMSE
---|---|---|---|---|---
Seq. | Method | Data | $p_{x}$ | $p_{y}$ | $p_{z}$ | roll | pitch | yaw | Pos (cm) | Rot (deg)
| Reference |  | -10.00 | 0.00 | 13.00 | 0.00 | 0.00 | 90.00 |  |
| FAST-LIO2 | [0,120] | 4.24 | 2.09 | 11.79 | 0.59 | 0.97 | 90.77 | 8.34 | 0.79
| ILC | [0,120] | -10.50 | -0.60 | 15.35 | -0.51 | 0.13 | 91.13 | 1.43 | 0.72
| W/O Obser. | [0,15] | -10.85 | -2.50 | 27.39 | -0.16 | 0.29 | 90.81 | 8.45 | 0.50
| W/O Obser. | [15,30] | -10.30 | -0.82 | 18.39 | -0.12 | 0.31 | 90.64 | 3.15 | 0.42
| W/O Obser. | [30,45] | -10.24 | -0.56 | 16.41 | -0.19 | 0.28 | 91.01 | 2.00 | 0.61
| W/O Obser. | [45,60] | -12.10 | -0.48 | 22.39 | -0.10 | 0.26 | 90.56 | 5.56 | 0.36
| W/ Obser. | [0,15] | -10.88 | 0.68 | 10.78 | -0.16 | 0.29 | 90.80 | 1.43 | 0.50
| W/ Obser. | [15,30] | -10.34 | 0.61 | 10.84 | -0.12 | 0.31 | 90.64 | 1.31 | 0.41
| W/ Obser. | [30,45] | -10.29 | 0.51 | 10.86 | -0.19 | 0.28 | 91.01 | 1.28 | 0.61
SKG-01 | W/ Obser. | [45,60] | -12.15 | 1.69 | 11.04 | -0.10 | 0.26 | 90.56 | 1.94 | 0.36
|  |  |  ${}^{I}\mathbf{p}_{L}$ (cm) |  ${}^{I}_{L}\bar{q}$ in Euler (deg) | RMSE
Seq. | Method | Data | $p_{x}$ | $p_{y}$ | $p_{z}$ | roll | pitch | yaw | Pos (cm) | Rot (deg)
| Reference |  | 1.00 | -0.66 | 9.47 | 0.02 | -0.37 | -0.14 |  |
| FAST-LIO2 | [85,325] | 0.96 | 1.53 | 7.33 | 1.22 | 0.91 | 1.51 | 1.77 | 1.39
| ILC | [85,325] | 11.05 | 10.67 | -46.97 | 4.42 | 5.38 | 9.01 | 33.74 | 6.74
| W/O Obser. | [85,100] | 0.34 | 0.27 | 26.72 | 0.16 | -0.34 | 0.27 | 9.98 | 0.25
| W/O Obser. | [100,115] | -0.04 | 1.53 | 7.28 | 0.31 | -0.46 | -0.29 | 1.88 | 0.19
| W/O Obser. | [115,130] | -2.78 | -1.65 | 64.22 | 0.10 | -0.41 | 1.25 | 31.69 | 0.80
| W/O Obser. | [130,145] | -0.70 | -0.08 | -1.99 | 0.19 | -0.22 | -0.12 | 6.70 | 0.13
| W/ Obser. | [85,100] | 0.21 | 0.25 | 7.52 | 0.16 | -0.34 | 0.26 | 1.32 | 0.24
| W/ Obser. | [100,115] | 0.06 | 1.37 | 7.47 | 0.31 | -0.45 | -0.29 | 1.73 | 0.19
| W/ Obser. | [115,130] | -3.32 | -1.71 | 7.56 | 0.10 | -0.41 | 1.30 | 2.79 | 0.83
Basement | W/ Obser. | [130,145] | -0.68 | -0.06 | 7.51 | 0.19 | -0.22 | -0.12 | 1.53 | 0.13

###  VII-D Calibration in Degenerate Case

In practical applications, the on-vehicle situation is a common case with
motion degeneration and there are high demands for accurate calibrations in
this scenario. In this section, we also test the extrinsic calibration in on-
vehicle conditions over the self-collected SKG-01 sequence and the publicly-
available Basement sequence. Regarding the initialization under degenerate
cases, since the compared methods ILC, LIOM, and FAST-LIO2 fail to be
initialized with the identity extrinsic rotation matrix, all the methods are
initialized with the perturbed extrinsics, with rotational perturbations of 2
degrees per axis and translational perturbations of 2 cm per axis. In general,
all the methods start from the same initial extrinsic parameters for fair
comparisons.

####  VII-D1 SKG-01 Sequence

To ensure the LiDAR sensor has sufficiently large visibility of the ground for
localization, we mount the sensor suite at a slight incline downwards as shown
in Fig. 10. We collect the SKG-01 sequence with a duration of 120 seconds,
while the ground robot moves on a planar ground surface. This sequence is
split into six 15-second segments, which are individually used for
calibration. Figure 16 illustrates the details of SKG-01 sequence and Fig. 17
displays translation calibration results of OA-LICalib with and without
observability awareness (Section V-B). Table XV reports the singular value of
each segment in SKG-01 sequence identified by OA-LICalib. For the first 6
segments, the eigen vector corresponding to the minimum singular value is
almost the same that is $[0.000,\ 0.000,\ 0.000,\ 0.007,\ -0.187,\ 0.982]$.
Note that the head part of the eigen vector corresponds to extrinsic rotation
and the tail part corresponds to extrinsic translation, which is determined by
the stack order of parameters $\delta\mathbf{x}$ in (32). The identified
degradation direction by the eigen vector is close to the $z$-axis, which is
reasonable and coincides with the fact that a certain angular offset exists
between the sensor frame and the ground plane. As shown in Fig. 17,
translation results in the $z$-axis of OA-LICalib without observability
awareness has a large deviation from the reference value. With the awareness
of observability, OA-LICalib is able to keep the main degenerate axis,
$z$-axis, as a prior and eliminate the coupled errors in the $x$-axis,
resulting in more accurate calibrated extrinsic parameters. However, for
segments at time intervals [90, 105] and [105,120], OA-LICalib estimates poor
results which is most likely due to insufficient excitation. From the angular
velocities curves of the SKG-01 sequence shown in Fig. 16, we can find that
the magnitudes of angular velocity are relatively small for these two
segments, leading to a deterioration of the observability, which is also
reflected by the singular values shown in Table XV.

####  VII-D2 Basement Sequence

We further test in the on-vehicle Basement sequence of the Hilti dataset. The
beginning of the Basement is skipped due to less motion. Fig. 18 reports the
calibration results. OA-LICalib with observability awareness automatically
identifies the principal non-observable direction as the $z$-axis of extrinsic
translation and succeeds in keeping it to the initial value. Besides, we
notice that the errors of final calibration results in the observable $x$-axis
and $y$-axis of translations are larger than those on the SKG-01 sequence.
Because the Hilti dataset is not specifically collected for calibration, and
the motion excitation in this sequence is insufficient, which is also
reflected by the small magnitude of averaged angular velocity shown in Table
VIII.

Table XVI reports the calibration results for the different methods. Again,
LIOM fails to give reasonable results thus are excluded from the table. The
main degenerate direction of SKG-01 and Basement sequences is the extrinsic
translation along the $z$-axis. Without the observability-aware module, our
method, FAST-LIO2, and ILC fail to give good calibrations along the $z$-axis
of extrinsic translation. With observability awareness, OA-LICalib is able to
generate accurate calibration results in on-vehicle conditions.

###  VII-E Remarks

It is clear from the previous simulation and experimental results that fully-
excited motions, whenever possible, are generally recommended for accurate LI
calibration, which agrees with the calibration literature [10, 54, 28]. For
example, the proposed OA-LICalib is able to provide accurate calibration
results from any segment with fully-excited motion in sequence Lab and Stairs,
as shown in Fig. 12 and Fig. 13 and Table IX. Among the sequences with
occasional fully-excited motion, e.g., YQ-01, IC_Office and Office_Mitte
sequences, the automatically-selected informative segments can also generate
reasonable calibration results as summarized in Table XI, XII. In contrast,
the OA-LICalib has larger errors on segments under degenerate planar motions,
such as [90,105] of the SKG-01 sequence (see Fig. 17) and [115,130] of the
Basement sequence (see Fig. 18), which are the segments with small angular
velocities (see Fig. 16).

While the fully-exited motions are recommended, they might not be feasible in
practice and it thus becomes necessary to address degeneracy. The proposed OA-
LICalib with observability-awareness is able to select informative segments
based on the singular values of the information matrix and tackle the
degeneracy via TSVD. Notably, the information metric of singular values is
affected by the sensory observations, motion profiles and environments. Some
other information metrics, for example, the trace, determinant, and maximal
eigenvalue of the information matrix can also be used [9].

It should be noted that the proposed OA-LICalib relies on the point-to-surfel
constraints, which well constrain the LI calibration problem in structured
environments. This implies that planar patches (surfels) exist in the
environments and can be identified. From our experience, no matter it is
indoors or outdoors, the OA-LICalib usually works well when there are enough
man-made structures (e.g., walls, buildings, staircases) within 50 meters and
the calibration data collection with fully-excited motion lasts for 10-20
seconds. However, in some cluttered environments without clear planar regions
or some scenarios with fairly sparse structures, the proposed approach would
have degraded performance. For example, it may fail to provide reasonable
results on some outdoor sequences of the Hilti dataset [55], in which the
scenarios are nearly empty without clear human-made structures. To address
this issue, resulting from insufficient point-to-surfel constraints from LiDAR
data, some ideas from advanced point cloud registration methods [57] could be
leveraged, such as Sdrsac [58], Teaser++ [59, 60], OPRANSAC [61], CVO [62],
and deep-learning-based methods [63].

##  VIII Conclusions and Future Work

In this paper, we have developed an observability-aware targetless LiDAR-IMU
calibration method, termed OA-LICalib, within the continuous-time batch
optimization framework. The proposed OA-LICalib calibrates not only the
spatial-temporal extrinsic parameters but also the intrinsic parameters of
both IMU and LiDAR sensors, while enforcing observability constraints during
update to address possible degeneracy. Specifically, two observability-aware
strategies are employed: (i) the informative data segment selection, and (ii)
the observability-aware update in back-end optimization. The former selects
the most informative data segments automatically for calibration among the
long-session data sequence, which lowers the computational consumption for
calibration as well as helps non-expert end users. The latter addresses the
degenerate motions via truncated SVD (TSVD) that is used to update only the
observable directions of the sensing parameters. The proposed method has been
extensively validated on both Monte-Carlo simulations and real-world
experiments, showing that the proposed OA-LICalib is able to provide accurate
spatial-temporal extrinsic and intrinsic calibration with high repeatability,
even under certain degenerate cases. For future work, the active calibration
guiding users to collect informative data deserves to be investigated.

## Acknowledgments

This work is supported by a Grant from the National Natural Science Foundation
of China (No. U21A20484) and Huang was partially supported by the University
of Delaware College of Engineering.

## References

  * [1] Ji Zhang and Sanjiv Singh  “LOAM: Lidar Odometry and Mapping in Real-time.”  In _Robotics: Science and Systems_ 2, 2014, pp. 9
  * [2] Tixiao Shan and Brendan Englot  “LeGO-lOAM: Lightweight and ground-optimized lidar odometry and mapping on variable terrain”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2018, pp. 4758–4765  IEEE
  * [3] Haoyang Ye, Yuying Chen and Ming Liu  “Tightly coupled 3d lidar inertial odometry and mapping”  In _2019 International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 3144–3150  IEEE
  * [4] Xingxing Zuo, Yulin Yang, Patrick Geneva, Jiajun Lv, Yong Liu, Guoquan Huang and Marc Pollefeys  “Lic-fusion 2.0: Lidar-inertial-camera odometry with sliding-window plane-feature tracking”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2020, pp. 5112–5119  IEEE
  * [5] Tixiao Shan, Brendan Englot, Drew Meyers, Wei Wang, Carlo Ratti and Daniela Rus  “Lio-sam: Tightly-coupled lidar inertial odometry via smoothing and mapping”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2020, pp. 5135–5142  IEEE
  * [6] Jiarong Lin, Xiyuan Liu and Fu Zhang  “A decentralized framework for simultaneous calibration, localization and mapping with multiple LiDARs”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2020, pp. 123–130  IEEE
  * [7] Jiajun Lv, Kewei Hu, Jinhong Xu, Yong Liu, Xiushui Ma and Xingxing Zuo  “CLINS: Continuous-Time Trajectory Estimation for LiDAR-Inertial System”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2021, pp. 6657–6663  IEEE
  * [8] Mingyang Li, Hongsheng Yu, Xing Zheng and Anastasios I Mourikis  “High-fidelity sensor modeling and self-calibration in vision-aided inertial navigation”  In _2014 IEEE International Conference on Robotics and Automation (ICRA)_ , 2014, pp. 409–416  IEEE
  * [9] Thomas Schneider, Mingyang Li, Cesar Cadena, Juan Nieto and Roland Siegwart  “Observability-aware self-calibration of visual and inertial sensors for ego-motion estimation”  In _IEEE Sensors Journal_ 19.10 IEEE, 2019, pp. 3846–3860
  * [10] Yulin Yang, Patrick Geneva, Xingxing Zuo and Guoquan Huang  “Online IMU Intrinsic Calibration: Is It Necessary?”  In _Proc. of the Robotics: Science and Systems, Corvallis, Oregon_ , 2020, pp. 716–20
  * [11] Gerardo Atanacio-Jiménez, José-Joel González-Barbosa, Juan B Hurtado-Ramos, Francisco J Ornelas-Rodríguez, Hugo Jiménez-Hernández, Teresa García-Ramirez and Ricardo González-Barbosa  “Lidar velodyne hdl-64e calibration using pattern planes”  In _International Journal of Advanced Robotic Systems_ 8.5 SAGE Publications Sage UK: London, England, 2011, pp. 59
  * [12] Jesse Levinson and Sebastian Thrun  “Unsupervised calibration for multi-beam lasers”  In _Experimental Robotics_ , 2014, pp. 179–193  Springer
  * [13] Joern Rehder, Janosch Nikolic, Thomas Schneider, Timo Hinzmann and Roland Siegwart  “Extending kalibr: Calibrating the extrinsics of multiple IMUs and of individual axes”  In _Proc. International Conference on Robotics and Automation_ , 2016, pp. 4304–4311  IEEE
  * [14] WI Liu and Yunwang Li  “Error modeling and extrinsic–intrinsic calibration for LiDAR-IMU system based on cone-cylinder features”  In _Robotics and Autonomous Systems_ 114 Elsevier, 2019, pp. 124–133
  * [15] Wanli Liu, Zhixiong Li, Reza Malekian, Miguel Sotelo, Zhenjun Ma and Weihua Li  “A Novel Multi-Feature based On-Site Calibration Method for LiDAR-IMU Systems”  In _IEEE Transactions on Industrial Electronics_ IEEE, 2019
  * [16] Jiajun Lv, Jinhong Xu, Kewei Hu, Yong Liu and Xingxing Zuo  “Targetless calibration of lidar-imu system based on continuous-time batch estimation”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2020, pp. 9968–9975  IEEE
  * [17] Faraz M Mirzaei, Dimitrios G Kottas and Stergios I Roumeliotis  “3D LIDAR-camera intrinsic and extrinsic calibration: Identifiability and analytical least-squares-based initialization”  In _The International Journal of Robotics Research_ 31.4 SAGE Publications Sage UK: London, England, 2012, pp. 452–467
  * [18] Andreas Geiger, Philip Lenz, Christoph Stiller and Raquel Urtasun  “Vision meets robotics: The KITTI dataset”  In _The International Journal of Robotics Research_ 32.11 Sage Publications Sage UK: London, England, 2013, pp. 1231–1237
  * [19] Radu Horaud and Fadi Dornaika  “Hand-eye calibration”  In _The International Journal of Robotics Research_ 14.3 Sage Publications Sage CA: Thousand Oaks, CA, 1995, pp. 195–210
  * [20] Cedric Le Gentil, Teresa Vidal-Calleja and Shoudong Huang  “3D Lidar-IMU Calibration based on Upsampled Preintegrated Measurements for Motion Distortion Correction”  In _Proc. International Conference on Robotics and Automation_ , 2018, pp. 2149–2155  IEEE
  * [21] Christopher KI Williams and Carl Edward Rasmussen  “Gaussian processes for machine learning”  MIT press Cambridge, MA, 2006
  * [22] Subodh Mishra, Gaurav Pandey and Srikanth Saripalli  “Target-free Extrinsic Calibration of a 3D-Lidar and an IMU”  In _2021 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems (MFI)_ , 2021, pp. 1–7  IEEE
  * [23] Patrick Geneva, Kevin Eckenhoff, Woosik Lee, Yulin Yang and Guoquan Huang  “OpenVINS: A Research Platform for Visual-Inertial Estimation”  In _Proc. of the IEEE International Conference on Robotics and Automation_ , 2020  URL: \url{https://github.com/rpng/open_vins}
  * [24] Xingxing Zuo, Patrick Geneva, Woosik Lee, Yong Liu and Guoquan Huang  “LIC-Fusion: LiDAR-Inertial-Camera Odometry”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2019, pp. 5848–5854  DOI: 10.1109/IROS40897.2019.8967746
  * [25] Kejie Qiu, Tong Qin, Jie Pan, Siqi Liu and Shaojie Shen  “Real-time temporal and rotational calibration of heterogeneous sensors using motion correlation analysis”  In _IEEE Transactions on Robotics_ 37.2 IEEE, 2020, pp. 587–602
  * [26] Wei Xu, Yixi Cai, Dongjiao He, Jiarong Lin and Fu Zhang  “Fast-lio2: Fast direct lidar-inertial odometry”  In _arXiv preprint arXiv:2107.06829_ , 2021
  * [27] Joern Rehder, Roland Siegwart and Paul Furgale  “A general approach to spatiotemporal calibration in multisensor systems”  In _IEEE Transactions on Robotics_ 32.2 IEEE, 2016, pp. 383–398
  * [28] Paul Furgale, Timothy D Barfoot and Gabe Sibley  “Continuous-time batch estimation using temporal basis functions”  In _Proc. International Conference on Robotics and Automation_ , 2012, pp. 2088–2095  IEEE
  * [29] Paul Furgale, Joern Rehder and Roland Siegwart  “Unified temporal and spatial calibration for multi-sensor systems”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2013, pp. 1280–1286  IEEE
  * [30] Joern Rehder, Paul Beardsley, Roland Siegwart and Paul Furgale  “Spatio-temporal laser to visual/inertial calibration with applications to hand-held, large scale scanning”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2014, pp. 459–465  IEEE
  * [31] Jérôme Maye, Paul Furgale and Roland Siegwart  “Self-supervised calibration for robotic systems”  In _2013 IEEE Intelligent Vehicles Symposium (IV)_ , 2013, pp. 473–480  IEEE
  * [32] Ji Zhang, Michael Kaess and Sanjiv Singh  “On degeneracy of optimization-based state estimation problems”  In _2016 IEEE International Conference on Robotics and Automation (ICRA)_ , 2016, pp. 809–816  IEEE
  * [33] Jianhao Jiao, Haoyang Ye, Yilong Zhu and Ming Liu  “Robust odometry and mapping for multi-lidar systems with online extrinsic calibration”  In _IEEE Transactions on Robotics_ IEEE, 2021
  * [34] Christiane Sommer, Vladyslav Usenko, David Schubert, Nikolaus Demmel and Daniel Cremers  “Efficient Derivative Computation for Cumulative B-Splines on Lie Groups”  In _arXiv preprint arXiv:1911.08860_ , 2019
  * [35] Adrian Haarbach, Tolga Birdal and Slobodan Ilic  “Survey of higher order rigid body motion interpolation methods for keyframe animation and continuous-time trajectory estimation”  In _2018 International Conference on 3D Vision (3DV)_ , 2018, pp. 381–389  IEEE
  * [36] Hannes Ovrén and Per-Erik Forssén  “Trajectory representation and landmark projection for continuous-time structure from motion”  In _The International Journal of Robotics Research_ 38.6 SAGE Publications Sage UK: London, England, 2019, pp. 686–701
  * [37] Christiane Sommer, Vladyslav Usenko, David Schubert, Nikolaus Demmel and Daniel Cremers  “Efficient Derivative Computation for Cumulative B-Splines on Lie Groups”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 11148–11156
  * [38] Myoung-Jun Kim, Myung-Soo Kim and Sung Yong Shin  “A general construction scheme for unit quaternion curves with simple high order derivatives”  In _Proceedings of the 22nd annual conference on Computer graphics and interactive techniques_ , 1995, pp. 369–376
  * [39] Dave Zachariah and Magnus Jansson  “Joint calibration of an inertial measurement unit and coordinate transformation parameters using a monocular camera”  In _2010 International Conference on Indoor Positioning and Indoor Navigation_ , 2010, pp. 1–7  IEEE
  * [40] Craig Glennie and Derek D Lichti  “Static calibration and analysis of the Velodyne HDL-64E S2 for high accuracy mobile scanning”  In _Remote sensing_ 2.6 Molecular Diversity Preservation International, 2010, pp. 1610–1624
  * [41] Patrick Geneva, Kevin Eckenhoff, Yulin Yang and Guoquan Huang  “LIPS: Lidar-inertial 3d plane slam”  In _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2018, pp. 123–130  IEEE
  * [42] Martin Magnusson  “The three-dimensional normal-distributions transform: an efficient representation for registration, surface analysis, and loop detection”, 2009
  * [43] Zhenfei Yang and Shaojie Shen  “Monocular visual–inertial state estimation with online initialization and camera–IMU extrinsic calibration”  In _IEEE Transactions on Automation Science and Engineering_ 14.1 IEEE, 2016, pp. 39–51
  * [44] Nikolas Trawny and Stergios I Roumeliotis  “Indirect Kalman filter for 3D attitude estimation”  In _University of Minnesota, Dept. of Comp. Sci. & Eng., Tech. Rep_ 2 Citeseer, 2005, pp. 2005
  * [45] Michael Bosse and Robert Zlot  “Continuous 3D scan-matching with a spinning 2D laser”  In _Proc. International Conference on Robotics and Automation_ , 2009, pp. 4312–4319  IEEE
  * [46] DQF De Menezes, Diego M Prata, Argimiro R Secchi and José Carlos Pinto  “A review on robust M-estimators for regression analysis”  In _Computers & Chemical Engineering_ Elsevier, 2021, pp. 107254
  * [47] Karol Hausman, James Preiss, Gaurav S Sukhatme and Stephan Weiss  “Observability-aware trajectory optimization for self-calibration with application to uavs”  In _IEEE Robotics and Automation Letters_ 2.3 IEEE, 2017, pp. 1770–1777
  * [48] Michael Kaess, Hordur Johannsson, Richard Roberts, Viorela Ila, John J Leonard and Frank Dellaert  “iSAM2: Incremental smoothing and mapping using the Bayes tree”  In _The International Journal of Robotics Research_ 31.2 Sage Publications Sage UK: London, England, 2012, pp. 216–235
  * [49] Kevin Eckenhoff, Patrick Geneva and Guoquan Huang  “Closed-form preintegration methods for graph-based visual–inertial navigation”  In _The International Journal of Robotics Research_ 38.5 SAGE Publications Sage UK: London, England, 2019, pp. 563–586
  * [50] C. Jauffret  “Observability and fisher information matrix in nonlinear regression”  In _Aerospace and Electronic Systems, IEEE Transactions on_ 43.2, 2007, pp. 756 –759  DOI: 10.1109/TAES.2007.4285368
  * [51] Per Christian Hansen  “The truncatedsvd as a method for regularization”  In _BIT Numerical Mathematics_ 27.4 Springer, 1987, pp. 534–553
  * [52] Jérôme Maye, Hannes Sommer, Gabriel Agamennoni, Roland Siegwart and Paul Furgale  “Online self-calibration for robotic systems”  In _The International Journal of Robotics Research_ 35.4 SAGE Publications Sage UK: London, England, 2016, pp. 357–380
  * [53] David Droeschel, Jörg Stückler and Sven Behnke  “Local multi-resolution representation for 6D motion estimation and mapping with a continuously rotating 3D laser scanner”  In _2014 IEEE International Conference on Robotics and Automation (ICRA)_ , 2014, pp. 5221–5226  IEEE
  * [54] Jonathan Kelly and Gaurav S Sukhatme  “Visual-inertial sensor fusion: Localization, mapping and sensor-to-sensor self-calibration”  In _The International Journal of Robotics Research_ 30.1 SAGE Publications Sage UK: London, England, 2011, pp. 56–79
  * [55] Michael Helmberger, Kristian Morin, Nitish Kumar, Danwei Wang, Yufeng Yue, Giovanni Cioffi and Davide Scaramuzza  “The Hilti SLAM Challenge Dataset”  In _arXiv preprint arXiv:2109.11316_ , 2021
  * [56] Michael Grupp  “evo: Python package for the evaluation of odometry and SLAM.”, https://github.com/MichaelGrupp/evo, 2017
  * [57] Xiaoshui Huang, Guofeng Mei, Jian Zhang and Rana Abbas  “A comprehensive survey on point cloud registration”  In _arXiv preprint arXiv:2103.02690_ , 2021
  * [58] Huu M Le, Thanh-Toan Do, Tuan Hoang and Ngai-Man Cheung  “Sdrsac: Semidefinite-based randomized approach for robust point cloud registration without correspondences”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 124–133
  * [59] Heng Yang and Luca Carlone  “A polynomial-time solution for robust registration with extreme outlier rates”  In _arXiv preprint arXiv:1903.08588_ , 2019
  * [60] Heng Yang, Jingnan Shi and Luca Carlone  “Teaser: Fast and certifiable point cloud registration”  In _IEEE Transactions on Robotics_ 37.2 IEEE, 2020, pp. 314–333
  * [61] Jiayuan Li, Qingwu Hu and Mingyao Ai  “Point cloud registration based on one-point ransac and scale-annealing biweight estimation”  In _IEEE Transactions on Geoscience and Remote Sensing_ IEEE, 2021
  * [62] Ray Zhang, Tzu-Yuan Lin, Chien Erh Lin, Steven A Parkison, William Clark, Jessy W Grizzle, Ryan M Eustice and Maani Ghaffari  “A new framework for registration of semantic point clouds from stereo and RGB-D cameras”  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 12214–12221  IEEE
  * [63] Shengyu Huang, Zan Gojcic, Mikhail Usvyatsov, Andreas Wieser and Konrad Schindler  “PREDATOR: Registration of 3D Point Clouds with Low Overlap”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 4267–4276