# GRIL-Calib: Targetless Ground Robot IMU-LiDAR Extrinsic Calibration Method
using Ground Plane Motion Constraints

TaeYoung Kim1, Gyuhyeon Pak2 and Euntai Kim2,∗ ∗ is that the corresponding
author.1TaeYoung Kim is the Department of Vehicle Convergence Engineering,
Yonsei University, Seoul 03722, South Korea<EMAIL_ADDRESS>Pak
and Euntai Kim are with the Department of Electrical and Electronic
Engineering, Yonsei University, Seoul 03722, South Korea {ghpak,
<EMAIL_ADDRESS>
###### Abstract

Targetless IMU-LiDAR extrinsic calibration methods are gaining significant
attention as the importance of the IMU-LiDAR fusion system increases. Notably,
existing calibration methods derive calibration parameters under the
assumption that the methods require full motion in all axes. When IMU and
LiDAR are mounted on a ground robot the motion of which is restricted to
planar motion, existing calibration methods are likely to exhibit degraded
performance. To address this issue, we present GRIL-Calib: a novel targetless
Ground Robot IMU-LiDAR Calibration method. Our proposed method leverages
ground information to compensate for the lack of unrestricted full motion.
First, we propose LiDAR Odometry (LO) using ground plane residuals to enhance
calibration accuracy. Second, we propose the Ground Plane Motion (GPM)
constraint and incorporate it into the optimization for calibration, enabling
the determination of full 6-DoF extrinsic parameters, including theoretically
unobservable direction. Finally, unlike baseline methods, we formulate the
calibration not as sequential two optimizations but as a single optimization
(SO) problem, solving all calibration parameters simultaneously and improving
accuracy. We validate our GRIL-Calib by applying it to three public real-world
datasets and comparing its performance with that of existing state-of-the-art
methods in terms of accuracy and robustness. Our code is available at
https://github.com/Taeyoung96/GRIL-Calib.

###### Index Terms:

Calibration and identification, Sensor Fusion, Optimization, Mobile Robot.

##  I INTRODUCTION

LiDAR (Light Detection And Ranging) acquires more accurate 3D information than
other sensors, making it broadly applicable to a variety of robotic
applications such as ego-motion estimation [1], localization [2] and mapping
[3] system. However, LiDAR has the drawback that when it is applied to fast
motion, it is vulnerable to motion distortion. Recently, to overcome the
drawback, LiDAR has been used in conjunction with IMU (Inertial Measurement
Unit) [4], [5], [6]. The IMU effectively compensates for the motion distortion
of the LiDAR, as it measures acceleration and angular velocity at a high
frequency. Obviously, accurate extrinsic calibration is a fundamental
requirement in the IMU-LiDAR fusion system. Extrinsic calibration refers to
the process of establishing the spatial transformations between IMU and LiDAR,
thereby aligning the measurements from the two sensors. If the extrinsic
calibration is not accurate enough, the fusion of IMU and LiDAR becomes less
effective, leading to significant degradation in overall performance [7].

Figure 1:  Illustration of a ground robot equipped IMU ($I$) sensor and LiDAR
($L$) sensor. We leverage the IMU height from the ground ($\hat{d}_{I}$) as
minimal prior knowledge and the LiDAR ground segmentation points
($\mathcal{G}$) for IMU-LiDAR extrinsic calibration.

Therefore, a lot of research has been conducted regarding the extrinsic
calibration between IMU and LiDAR, and the research can be divided into two
directions: target-based and targetless approaches. Target-based approach uses
known, predefined objects such as cylinders [8] to estimate the spatial
transformations between IMU and LiDAR. Target-based approach is relatively
accurate, but it has the shortcoming that predetermined targets are required.
On the other hand, the targetless approach does not use any predefined objects
for calibration, receiving much attention due to ease of use and
repeatability. Instead, the targetless method requires comprehensive
excitation motion of both sensors to activate all Degrees of Freedom (DoF).
Achieving this motion is straightforward for handheld devices with integrated
IMU and LiDAR sensors, as they can easily perform a full range of motions.
However, for ground robots, it becomes a significant challenge because their
motion is limited to planar motion. Furthermore, from a theoretical
perspective, the lack of full motion will make the rank of Fisher information
matrix (FIM) less than its dimension, and the degenerate or rank-deficient FIM
will cause the estimated calibration parameters to have high uncertainty or
not to become unique. To address the issue, some papers perform observability
analysis [9], [10], [11], but the analysis only identifies unobservable
directions and does not directly estimate calibration parameters.

In this paper, we present GRIL-Calib: novel targetless Ground Robot IMU-LiDAR
extrinsic Calibration method. To address the aforementioned degenerate
problem, we utilize a ground plane motion (GPM) constraint. The GPM provides
geometric information about the robot’s motion and allows the full 6-DoF
extrinsic parameters even when some of the extrinsic parameters are
unobservable. The extrinsic calibration problem between IMU and LiDAR for
ground robots is depicted in Fig 1. Similar constraints related to the ground
have been proposed for LiDAR odometry (LO) [12] or visual-inertial odometry
(VIO) [13]. However, we propose a ground plane motion (GPM) constraint
specifically focused on the extrinsic calibration of the IMU and LiDAR. The
proposed method leverages a motion-based approach and integrates ground
observation, thereby offering a user-friendly and accurate solution for
extrinsic calibration of the IMU and LiDAR mounted on the ground robot. The
key contributions of this paper are summarized as follows:

  * •

We propose a novel extrinsic calibration method of IMU-LiDAR fusion system on
a ground robot and thus does not need to satisfy full excitation motion.

  * •

To improve the performance of the calibration, we proposed a novel LiDAR
Odometry (LO) that exploits ground plane observations. The ground plane
residuals are incorporated into an Iterative Error State Kalman Filter (IESKF)
to reduce drift around the z-axis.

  * •

We propose Ground Plane Motion (GPM) constraints and apply the constraint to
the nonlinear optimization for calibration. The GPM enables us to address
parameters including unobservable direction and also improve calibration
accuracy.

  * •

In the nonlinear optimization for calibration, we formulate the calibration of
the translation and rotation parameters as a Single Optimization (SO) problem
solve parameters simultaneously. The baseline methods, however, calibrate the
two parameters separately and sequentially, degrading the calibration
accuracy.

  * •

We validate our GRIL-Calib on M2DGR [14], Hilti [15], and S3E [16] datasets
where it achieves excellent results compared to state-of-the-art targetless
methods. We have open-sourced111https://github.com/Taeyoung96/GRIL-Calib our
algorithm to benefit the broader robotics community.

The rest part of this paper is organized as follows: Section II presents a
comprehensive review of recent IMU-LiDAR calibration methods. In Section III,
we present the system overview of GRIL-Calib and describe each key part in
detail at Section IV and Section V. Section VI describes the experiments
performed on various real-world datasets and the analysis of the results.
Finally, Section VII provides conclusions and discusses future work.

##  II Related Work

In early IMU-LiDAR calibration research, an arbitrary target is required for
calibration [8], [17]. However, the target-based method is not practical, so
targetless calibration methods have been studied more recently. Subodh et al.
[18] formulate a motion-based calibration based on the Extended Kalman Filter
(EKF). Additionally, Lv et al. [19] utilizes a continuous-time trajectory
formulation based on B-Spline. Also, OA-LICalib [11] which extends their
previous works [19], proposes two observability-aware modules designed to
address the degenerate case. Unfortunately, the proposed observability-aware
module is designed to only reject unnecessary motion. FAST-LIO2 [4] is an
accurate LiDAR-inertial odometry (LIO) algorithm based on an Iterated error
state Kalman filter (IESKF) which include extrinsic calibration parameters.
Due to the strong nonlinearity, the performance of the LIO system is highly
dependent on accurate initial states. To get an accurate estimate of the
initial states, Zhu et al. [20] propose the optimization-based algorithm that
aligns angular velocity and acceleration. Existing approaches [18], [11], [20]
require sufficient motion in all axes of the sensor to obtain accurate
results. These methods are practical when applied to handheld devices that can
move freely in various directions. However, existing methods encounter
challenges when applied to ground robots where motion is limited to planar
motion. In contrast, our method utilizes ground observation that can easily be
observed by a ground robot to obtain accurate 6-DoF extrinsic parameters, even
in the absence of sufficient motion in all axes. Recently, an extended Kalman
filter (EKF) based IMU-LiDAR calibration algorithm utilizing a GNSS
observation model in planar motion has been proposed [7]. However, it requires
an additional GNSS sensor, while our proposed algorithm calibrates IMU-LiDAR
without additional sensors.

##  III System Overview

The overall system of GRIL-Calib comprises two primary stages: pre-processing
and optimization, as shown in Fig. 2. Our framework is motivated by [20]. The
first stage is pre-processing, which generates additional information that is
critical for the following optimization stage. Specifically, the pre-
processing stage performs the following key tasks:

  * •

Ground Observation via ground segmentation, where ground points are extracted
from LiDAR data, facilitating the estimation of the LiDAR sensor’s relative
orientation to the ground.

  * •

Lidar Odometry (LO), which enhances LO performance by utilizing the extracted
ground plane and incorporating a proposed ground plane residual.

  * •

IMU Processing, involving the pre-processing of IMU measurements and the
calculation of the IMU’s relative orientation to the ground using the Madgwick
filter [21].

  * •

Data Synchronization, aimed at determining the approximate time offset between
IMU and LiDAR data using cross-correlation methods. Notably, this process
aligns with the approach proposed in [20].

The second stage is optimization, which aims to optimize the extrinsic
calibration parameters by utilizing the data acquired in the preprocessing
stage. In this optimization stage, we apply three types of constraints:
unified temporal-spatial constraints based on angular velocity and linear
acceleration, as well as proposed Ground Plane Motion (GPM) constraints. The
GPM constraint ensures that the robot’s motion remains on the ground and is
integrated into the objective function to achieve precise extrinsic
calibration while limiting planar motion. In contrast [20], where each type of
objective function is optimized individually and sequentially, we
simultaneously optimize all three types of constraints by assigning different
values of weights to each objective. This approach improves performance
compared to the method in [20].

Throughout the paper, we use the following notations: LiDAR frame denotes
$\left\\{L\right\\}$, IMU frame denotes $\left\\{I\right\\}$,
$\mathbf{R}^{I}_{L}\in SO(3)$ is the rotation matrix from IMU frame to LiDAR
frame. $\mathbf{p}^{I}_{L}\in\mathbb{R}^{3}$ is the translation vector from
the IMU frame to the LiDAR frame. $t^{I}_{L}\in\mathbb{R}$ is the time offset
between IMU and LiDAR. $\mathbf{b}_{a},\mathbf{b}_{g}\in\mathbb{R}^{3}$ means
the acceleration bias and angular velocity bias of the IMU, respectively.

Figure 2: The overview of proposed GRIL-Calib system. During the optimization
stage (Green Box), we use 3 constraints. See Section V for more details.

##  IV Preprocessing

In the preprocessing stage, we obtain some parameters and the sensor
measurements which will be used in the subsequent the optimization stage.

###  IV-A Ground Observation via LiDAR Ground Segmentation

First, the ground is segmented from the LiDAR input. LiDAR ground segmentation
aims to obtain the rotation matrix $\mathbf{R}^{G}_{L}$ of the LiDAR with
respect to the ground. When the current LiDAR scan $\mathcal{P}$ is presented,
we extract a set of ground points $\mathcal{G}$ using Patchwork++ [22], a
state-of-the-art (SOTA) ground segmentation algorithm. Then, we apply the
Principal Component Analysis (PCA) to the covariance matrix
$\mathcal{C}_{\mathcal{G}}$ of the ground points $\mathcal{G}$ identifying a
normal vector $\mathbf{n}_{L}=(a_{L},b_{L},c_{L})^{T}$ of the ground plane

$a_{L}x+b_{L}y+c_{L}z+d_{L}=0$ (1)

where $\mathbf{n}_{L}$ implies a normal vector to the ground in the LiDAR
coordinate $\left\\{L\right\\}$. We also define a normal vector
$\mathbf{n}_{G}=(0,0,1)^{T}$ with respect to the ground coordinate
$\left\\{G\right\\}$. From the two normal vectors $\mathbf{n}_{L}$ and
$\mathbf{n}_{G}$ , we obtain a rotation matrix $\mathbf{R}^{G}_{L}$ from
$\left\\{G\right\\}$ frame to $\left\\{L\right\\}$ frame by

$\mathbf{R}_{L}^{G}=\exp({{\left[{{\phi}_{GL}}{{\mathbf{n}}_{GL}}\right]}_{\times}})$
(2)

where
${{\phi}_{GL}}=\arccos({{\mathbf{n}}_{G}}^{T}{{\mathbf{n}}_{L}})\in\mathbb{R}$
and
${{\mathbf{n}}_{GL}}={{\mathbf{n}}_{G}}\times{{\mathbf{n}}_{L}}\in{{\mathbb{R}}^{3}}$
denote the rotation angle and vector between $\left\\{G\right\\}$ and
$\left\\{L\right\\}$, respectively. $\left[\cdot\right]_{\times}$ denotes the
skew-symmetric matrix, and $\exp\left(\cdot\right)$ denotes the exponential
map.

###  IV-B LiDAR Odometry Utilizing the Ground Plane Residual

We use Iterated Error state Kalman Filter (IESKF) to obtain LiDAR Odometry
(LO) with a similar procedure to [20]. However, we improved the performance of
LO by exploiting the assumption that ground robot is constrained to movement
on a 2D plane. When we define the state vector $\mathbf{x_{k}}$ of the LiDAR
by

${{\mathbf{x}}_{k}}={{\left[\log\left({{\mathbf{R}}_{{{L}_{k}}}}\right),{{\mathbf{p}}_{{{L}_{k}}}},{{\boldsymbol{\omega}}_{{{L}_{k}}}},{{\mathbf{v}}_{{{L}_{k}}}}\right]}^{T}}$
(3)

and its motion as follows,

$\mathbf{x}_{k+1}=\mathbf{x}_{k}\boxplus\left(\Delta
t\cdot\mathbf{F}\left(\mathbf{x}_{k},\mathbf{w}_{k}\right)\right)$ (4)

where
$\mathbf{F}({{\mathbf{x}}_{\mathbf{k}}},{{\mathbf{w}}_{\mathbf{k}}})={{\left[{{\boldsymbol{\omega}}_{{{L}_{k}}}},{{\mathbf{v}}_{{{L}_{k}}}},\boldsymbol{{{\varepsilon}}_{\boldsymbol{\omega}}},{\boldsymbol{{\varepsilon}}_{\mathbf{v}}}\right]}^{T}}$.
$\mathbf{R}_{L_{k}}$represents the attitude from the initial LiDAR frame
$\left\\{{{L}_{0}}\right\\}$ to the LiDAR frame $\left\\{{{L}_{k}}\right\\}$
at current time $k$, that is
${{\mathbf{R}}_{{{L}_{k}}}}=\mathbf{R}_{{{L}_{k}}}^{{{L}_{0}}}$,
$\log\left(\cdot\right)$ denotes the logarithmic map.
${{\mathbf{p}}_{{{L}_{k}}}}=\mathbf{p}_{{{L}_{k}}}^{{{L}_{0}}}$ represents the
position of the LiDAR at current time $k$.
${{\boldsymbol{\omega}}_{{{L}_{k}}}}$ and $\mathbf{v}_{L_{k}}$ represent the
angular and linear velocities of the LiDAR at time $k$, respectively.
${{\mathbf{w}}_{k}}=\left[{\boldsymbol{{\varepsilon}}_{\boldsymbol{\omega}}},{\boldsymbol{{\varepsilon}}_{\mathbf{v}}}\right]\sim\mathcal{N}\left(\mathbf{0},{{\mathbf{Q}}_{k}}\right)$
are a random walk noise represented by zero-mean Gaussian. $\Delta t$ is the
time interval between two consecutive scans, and the notation $\boxplus$
represents the encapsulated “boxplus” operations on the state manifold [23].

The prediction step of the IESKF consists of the prediction of the state and
the propagation of the covariance as follows,

$\displaystyle\hat{\mathbf{x}}_{k+1}$
$\displaystyle=\bar{\mathbf{x}}_{k}\boxplus\left(\Delta
t\cdot\mathbf{F}\left(\bar{\mathbf{x}}_{k},\mathbf{0}\right)\right)$ (5)
$\displaystyle\hat{\mathbf{P}}_{k+1}$
$\displaystyle=\mathbf{F}_{\tilde{\mathbf{x}}}\bar{\mathbf{P}}_{k}\mathbf{F}_{\tilde{\mathbf{x}}}^{T}+\mathbf{F}_{\mathbf{w}}\mathbf{Q}_{k}\mathbf{F}_{\mathbf{w}}^{T}$

where $\hat{(\cdot)}$, $\tilde{(\cdot)}$, and $\bar{(\cdot)}$ denote the
predicted, error and optimal states, respectively; $\mathbf{P}_{k}$ and
$\mathbf{Q}_{k}$ are the covariance matrices for error state
$\tilde{\mathbf{x}}_{k}$ and noise $\mathbf{w}_{k}$, respectilvey.
$\mathbf{F}_{\tilde{\mathbf{x}}}$ and $\mathbf{F}_{\mathbf{w}}$ are the
derivative of $(\mathbf{x}_{k+1}\boxminus\hat{\mathbf{x}}_{k+1})$ with respect
to $\tilde{\mathbf{x}}$ and $\mathbf{w}$ evaluated at zero, respectively. The
notation $\boxminus$ represents the encapsulated “boxminus” operations on the
state manifold [23].

The correction step of IESKF is conducted by a Maximum A Posteriori (MAP)
problem as in [20], where the prior distribution is obtained through state
prediction, and the likelihood is obtained through the first-order
approximation of point-to-plane measurement model [4].

$\mathcal{P}{{}_{k}}=\left\\{p_{j}^{{{L}_{k}}}\;|\;j=1,\cdots,{{N}_{k}}\right\\}$
is given at time $k$, the measurement model [4] is constructed based on the
idea that a current point $p_{j}^{{{L}_{k}}}$ acquired at time $k$ should
remain on a corresponding neighboring plane ${{\pi}_{j}}$ in the map, as
described below:

$0={{\mathbf{h}}_{j}}\left({{\mathbf{x}}_{k}},p_{j}^{{{L}_{k}}}+\boldsymbol{\varepsilon}_{j}^{{{L}_{k}}}\right)={{\left(\mathbf{u}_{j}^{{{L}_{0}}}\right)}^{T}}\left(\mathbf{T}_{{{L}_{k}}}^{{{L}_{0}}}\left(p_{j}^{{{L}_{k}}}+\boldsymbol{\varepsilon}_{j}^{{{L}_{k}}}\right)-\mathbf{q}_{j}^{{{L}_{0}}}\right)$
(6)

where $\mathbf{u}_{j}^{{{L}_{0}}}$ is the unit normal vector of the nearest
corresponding plane defined w.r.t $\left\\{{{L}_{0}}\right\\}$;
$\mathbf{T}_{{{L}_{k}}}^{{{L}_{0}}}$ represent a LiDAR pose including
($\mathbf{R}_{L_{k}}$,$\mathbf{p}_{L_{k}}$), and the point
$\mathbf{q}_{j}^{{{L}_{0}}}$ is a LiDAR point located on the plane within the
map. By applying the first-order approximation to (6), we obtain

$\displaystyle 0$
$\displaystyle={{\mathbf{h}}_{j}}\left({{\mathbf{x}}_{k}},\;p_{j}^{{{L}_{k}}}+\boldsymbol{\varepsilon}_{j}^{{{L}_{k}}}\right)$
(7)
$\displaystyle\simeq{{\mathbf{h}}_{j}}\left(\hat{\mathbf{x}}_{k},\;p_{j}^{{{L}_{k}}}+0\right)+\mathbf{H}_{j}\tilde{\mathbf{x}}_{k}+{{\mathbf{v}}_{j}}=\mathbf{z}_{j}+\mathbf{H}_{j}\tilde{\mathbf{x}}_{k}+{{\mathbf{v}}_{j}}$

where
$\mathbf{z}_{j}={{\mathbf{h}}_{j}}\left(\hat{\mathbf{x}}_{k},p_{j}^{{{L}_{k}}}+0\right)$,
$\mathbf{H}_{j}={{\left.\frac{\partial{{\mathbf{h}}_{j}}}{\partial{{\mathbf{x}}_{k}}}\right|}_{{{\mathbf{x}}_{k}}={{{\mathbf{\tilde{x}}}}_{k}}}}$,
and
${{\mathbf{v}}_{j}}\in\mathcal{N}\left(\mathbf{0},{{\mathbf{R}}_{j}}\right)$
is the disturbance which comes from the raw measurement noise
$\boldsymbol{\varepsilon}_{j}^{{{L}_{k}}}$.

Different from [20], however, we enhance the LO by exploiting the planar
constraint which effectively reduces drift in the LO. Specifically, we further
use the constraint that the roll and pitch of the rotation matrix
${{\mathbf{R}}_{{{L}_{k}}}}$ and the z-component of the translation vector
$\mathbf{p}_{L_{k}}$ should be consistent because the motion of the robot is
restricted to a 2D space. Thus, when we are given a normal vector
${{\mathbf{n}}_{{{L}_{k}}}}$ of the LiDAR ground points ${{\mathcal{G}}_{k}}$
at time $k$, the vector ${{\mathbf{n}}_{{{L}_{k}}}}$ should be orthogonal to
the x-y plane of the ground frame $\left\\{G\right\\}$.

${{[\mathbf{R}_{L}^{G}{{\mathbf{R}}_{{{L}_{k}}}}({{\mathbf{n}}_{{{L}_{k}}}}+{{\boldsymbol{\varepsilon}}_{{{\mathbf{n}}_{k}}}})]}_{1,2}}=0$
(8)

where ${{\boldsymbol{\varepsilon}}_{{{\mathbf{n}}_{k}}}}$ is the observation
noise to ${{\mathbf{n}}_{{{L}_{k}}}}$. $[\cdot]_{1,2}$ denotes the notation
that takes only the first and second elements in a ($3\times 1$) vector is
used. Further, the displacement of the z-component of
${{\mathbf{p}}_{{{L}_{k}}}}$ should be zero respect to ground frame
$\left\\{G\right\\}$, that is,

$\mathbf{e}_{3}^{T}\mathbf{R}{{{}_{L}^{G}}}{{\mathbf{p}}_{{{L}_{k}}}}=0.$ (9)

By combining (8) and (9), we define a ground plane measurement model
$\mathbf{h}_{\pi}$ given by

$0={{\mathbf{h}}_{\pi}}\left({{\mathbf{x}}_{k}},{{\mathbf{n}}_{{{L}_{k}}}}+{{\boldsymbol{\varepsilon}}_{{{\mathbf{n}}_{k}}}}\right)=\left(\begin{matrix}{{\left[\mathbf{R}_{L}^{G}{{\mathbf{R}}_{{{L}_{k}}}}\left({{\mathbf{n}}_{{{L}_{k}}}}+{{\mathbf{\varepsilon}}_{{{\mathbf{n}}_{k}}}}\right)\right]}_{1,2}}\\\
\mathbf{e}_{3}^{T}\mathbf{R}_{L}^{G}{{\mathbf{p}}_{{{L}_{k}}}}\\\
\end{matrix}\right)$ (10)

Using (10), by applying the first-order approximation, we derive the following
equation,

$\displaystyle 0=$
$\displaystyle\;{{\mathbf{h}}_{\pi}}\left({{\mathbf{x}}_{k}},{{\mathbf{n}}_{{{L}_{k}}}}+{{\boldsymbol{\varepsilon}}_{{{\mathbf{n}}_{k}}}}\right)$
(11) $\displaystyle\simeq$
$\displaystyle\;{{\mathbf{h}}_{\pi}}\left(\hat{\mathbf{x}}_{k},{{\mathbf{n}}_{{{L}_{k}}}}+\mathbf{0}\right)+\mathbf{H}_{\pi}\tilde{\mathbf{x}}_{k}+{{\mathbf{v}}_{\pi}}=\mathbf{z}_{\pi}+\mathbf{H}_{\pi}\tilde{\mathbf{x}}_{k}+{{\mathbf{v}}_{\pi}}$

where
$\mathbf{z}_{\pi}={{\mathbf{h}}_{\pi}}\left(\hat{\mathbf{x}}_{k},{{\mathbf{n}}_{{{L}_{k}}}}+\mathbf{0}\right)$,
$\mathbf{H}_{\pi}={{\left.\frac{\partial{{\mathbf{h}}_{\pi}}}{\partial{{\mathbf{x}}_{k}}}\right|}_{{{\mathbf{x}}_{k}}={{{\mathbf{\tilde{x}}}}_{k}}}}$
,and $\mathbf{v_{\pi}}\sim\mathcal{N}\left(\mathbf{0},\mathbf{S}_{k}\right)$.

By integrating the two likelihoods given by $\mathbf{h}_{j}$ in (7) and our
proposed $\mathbf{h}_{\pi}$ in (10) apply to MAP estimate, we obtain
${{\tilde{x}}_{k}}$ by

$\min_{\tilde{\mathbf{x}}_{k}^{\tau}}\left(\left\|\mathbf{x}_{k}\boxminus\hat{\mathbf{x}}_{k}\right\|_{\hat{\mathbf{P}}_{k}}^{2}+\sum_{j=1}^{m}\left\|\mathbf{z}_{j}^{\tau}+\mathbf{H}_{j}^{\tau}\tilde{\mathbf{x}}_{k}^{\tau}\right\|_{\mathbf{R}_{j}}^{2}+\left\|\mathbf{z}_{\pi}^{\tau}+\mathbf{H}_{\pi}^{\tau}\tilde{\mathbf{x}}_{k}^{\tau}\right\|_{\mathbf{S}_{k}}^{2}\right)$
(12)

where $\tau$ means the number of iterate times. $m$ is the number of LiDAR
points. When $\tilde{\mathbf{x}}_{k}$ decreases below a certain threshold, the
state is considered to be optimal. Based on the optimal state
$\bar{\mathbf{x}}_{k}$ we return,

$\mathbb{L}_{k}=\left\\{\boldsymbol{\omega}_{L_{k}},\mathbf{v}_{L_{k}},\boldsymbol{\alpha}_{L_{k}},\mathbf{a}_{L_{k}}\right\\}$
(13)

In (13), ${{\boldsymbol{\omega}}_{{{L}_{k}}}}$ and
${{\mathbf{v}}_{{{L}_{k}}}}$ are obtained from $\bar{\mathbf{x}}_{k}$, and the
angular acceleration $\boldsymbol{\alpha}_{L_{k}}$ and the linear acceleration
$\mathbf{a}_{L_{k}}$ of the LiDAR are computed by applying non-causal central
difference [24] to ${{\boldsymbol{\omega}}_{{{L}_{k}}}}$ and
${{\mathbf{v}}_{{{L}_{k}}}}$, respectively.

Figure 3: Visualization of timestamp for each IMU measurement and LiDAR
Odometry.

###  IV-C IMU Processing and Synchronization with LiDAR Odometry

In general, the IMU measurement model is defined as follows.

$\hat{\boldsymbol{\omega}}_{i}=\boldsymbol{\omega}_{i}+\mathbf{b}_{g}+\boldsymbol{\varepsilon}_{g},\quad\hat{\mathbf{a}}_{i}=\mathbf{a}_{i}+\mathbf{R}^{w}_{I_{i}}\mathbf{g}^{w}+\mathbf{b}_{a}+\boldsymbol{\varepsilon}_{a}$
(14)

The raw gyroscope and acceleration measurements at time $i$ denote as
$\hat{\boldsymbol{\omega}}_{i}$ and $\hat{\mathbf{a}}_{i}$. $\mathbf{b}_{g}$
is the gyroscope bias and $\boldsymbol{\varepsilon}_{g}$ is the gaussian white
noise for the gyroscope measurement. Also, $\mathbf{b}_{a}$ is the
acceleration bias and $\mathbf{g}^{w}$ is the gravity vector relative to world
frame $w$, and $\boldsymbol{\varepsilon}_{a}$ is a Gaussian white noise for
the acceleration measurement. Here, it should be noted that the IMU
measurements are much faster than LO, we use two different time indices $i$
and $k$ for IMU and LO, respectively, as shown in Fig 3.

Similar to [20], we remove the noise $\boldsymbol{\varepsilon}_{g}$ and
$\boldsymbol{\varepsilon}_{a}$ from the IMU measurements with a zero-phase
filter [25] in 14. Also, we down-sample to get the IMU measurement at
timestamp $k$ from the LiDAR.

$\boldsymbol{\bar{\omega}}_{I_{k}}=\boldsymbol{\omega}_{I_{k}}+\mathbf{b}_{g},\quad\bar{\mathbf{a}}_{I_{k}}=\mathbf{a}_{I_{k}}+\mathbf{R}^{w}_{I_{k}}\cdot\mathbf{g}^{w}+\mathbf{b}_{a}$
(15)

To achieve more precise synchronization between two sensors, we employ a
coarse-to-fine approach to determine the time offset, denoted as $t^{I}_{L}$.
As outlined in [20], $t^{I}_{L}$ is the sum of the coarse estimate $t_{c}$ and
the fine estimate $t_{f}$, thereby defined as $t^{I}_{L}=t_{c}+t_{f}$. Similar
to [20], we determine $t_{c}$ by using the cross-correlation method given by

$t_{c}^{*}=\underset{t_{c}}{\arg\max}\sum_{k=1}^{n}\left\|\bar{\boldsymbol{\omega}}_{I_{k+{t}_{c}}}\right\|\cdot\left\|\boldsymbol{\omega}_{L_{k}}\right\|$
(16)

In the following section, we derive the remaining component, $t_{f}$, as the
outcome of the optimization process. Further, by applying the IMU measurements
${{\boldsymbol{\omega}}_{{{I}_{k}}}}$ and ${{\mathbf{a}}_{{{I}_{k}}}}$ to the
Madgwick filter [21], we estimate the rotation matrix $\mathbf{R}^{G}_{I}$
from the ground frame $\left\\{G\right\\}$ to the IMU frame
$\left\\{I\right\\}$. Finally, as a result from the IMU processing and
synchronization procedure, we return

${{\mathbb{I}}_{k}}=\left\\{{{{\boldsymbol{\bar{\omega}}}}_{{{I}_{k}}}},{{{\bar{\boldsymbol{\alpha}}}}_{{{I}_{k}}}},{{{\mathbf{\bar{a}}}}_{{{I}_{k}}}},\mathbf{R}_{I}^{G}\right\\}$
(17)

##  V EXTRINSIC CALIBRATION VIA SINGLE OPTIMIZATION (SO)

Using the two sets of measurements ${{\mathbb{L}}_{k}}$ in (13) and
${{\mathbb{I}}_{k}}$ in (17) obtained from Section IV, we formulate the IMU-
LiDAR extrinsic calibration as an optimization problem. The calibration
parameter $\mathbf{x}_{cal}$ used in the calibration is given by,

$\mathbf{x}_{cal}=\left[\mathbf{R}^{I}_{L},\mathbf{p}^{I}_{L},t_{f},\mathbf{b}_{g},\mathbf{b}_{a}\right]$
(18)

###  V-A Unified Temporal-Spatial Constraints

We jointly utilize the two constraints proposed in [20]. The first constraint
uses the angular velocity measurements ${{\bar{\mathbf{\omega}}}_{{{I}_{k}}}}$
and ${{\mathbf{\omega}}_{{{L}_{k}}}}$ from IMU measurements and LO,
respectively.

$\bar{\boldsymbol{\omega}}_{I_{k^{\prime}}}+t_{f}\cdot\boldsymbol{\bar{\alpha}}_{I_{k^{\prime}}}=(\mathbf{R}^{I}_{L})^{T}\boldsymbol{\omega}_{L_{k}}+\mathbf{b}_{g}$
(19)

by compensating the time offset $k^{\prime}=k+t_{c}^{*}$ between two sensors,
where $t_{c}^{*}$ is given in 16 and assuming the IMU angular acceleration
$\boldsymbol{\bar{\alpha}}_{I_{k^{\prime}}}$ is constant over a small period
of time $t_{f}$. This is described as follows, according to (19)

$\mathbf{r}_{w}(\mathbf{R}^{I}_{L},\mathbf{b}_{g},\mathbf{t}_{f})=(\mathbf{R}^{I}_{L})^{\mathbf{T}}\boldsymbol{\omega}_{L_{k}}+\mathbf{b}_{g}-\bar{\boldsymbol{\omega}}_{I_{k^{\prime}}}-t_{f}\cdot\boldsymbol{\bar{\alpha}}_{I_{k^{\prime}}}$
(20)

The second constraint uses the acceleration measurements
${{{\mathbf{\bar{a}}}}_{{{I}_{k}}}}$ and ${{\mathbf{a}}_{{{L}_{k}}}}$ from IMU
and LO, respectively. As in [26], the accelerations
${{\mathbf{a}}_{{{I}_{k}}}}$ and ${{\mathbf{a}}_{{{L}_{k}}}}$of the two
sensors are related to each other by

$\mathbf{R}^{I}_{L}{\mathbf{\bar{a}}}_{I}=\mathbf{a}_{L}+\\{([\boldsymbol{\omega}_{L}]_{\times}^{2}+[\boldsymbol{\alpha}_{L}]_{\times})\cdot\mathbf{p}^{I}_{L}\\}$
(21)

where $\boldsymbol{\alpha}_{L}$ is obtained from ${{\mathbb{L}}_{k}}$. the
second constraint $\mathbf{r}_{a}$ is represented as follows.

$\displaystyle\mathbf{r}_{a}(\mathbf{R}^{I}_{L},\mathbf{p}^{I}_{L},\mathbf{b}_{a})$
$\displaystyle=\mathbf{R}^{I}_{L}(\mathbf{\bar{a}}_{I_{k^{\prime}}}-\mathbf{b}_{a})-\mathbf{a}_{L_{k}}$
(22)
$\displaystyle\quad-\\{([\boldsymbol{\omega}_{L_{k}}]_{\times}^{2}+[\boldsymbol{\alpha}_{L_{k}}]_{\times})\cdot\mathbf{p}^{I}_{L}\\}$

###  V-B Ground Plane Motion Constraints

In this subsection, we describe a novel constraint named Ground Plane Motion
(GPM). Our proposed constraint is slightly similar to the constraint proposed
in [13], where the plane motion information was leveraged to improve the
performance of VINS. However, the constraint in [13] was formulated for the
estimation of sensor motion, while our constraint is specifically designed for
extrinsic calibration between IMU and LiDAR. Specifically, our GPM is
parameterized in terms of $\mathbf{R}_{L}^{I}$ and $\mathbf{p}_{L}^{I}$. In
the formulation of the GPM, we directly measure the height of the IMU from the
ground $\hat{d_{I}}\in\mathbb{R}$ and use it as prior knowledge. Further, we
exploit the fact that the z-component of $\mathbf{p}{{{}_{L}^{I}}}$ could not
be observed because it is in the direction perpendicular to the planar motion.
Using the geometric relationship of the rotation parameters between the ground
and each sensor shown in Fig 4, we obtain

${{[\mathbf{R}_{I}^{G}\mathbf{R}_{L}^{I}\mathbf{R}_{G}^{L}\mathbf{e}_{3}]}_{1,2}}=0$
(23)
$\mathbf{e}_{3}^{T}{{(\mathbf{R}_{I}^{L}\mathbf{R}_{G}^{I}{{\mathbf{d}}_{I}}-\mathbf{R}_{G}^{L}{{\mathbf{d}}_{L}})}}=\mathbf{e}_{3}^{T}\mathbf{p}{{{}_{L}^{I}}}$
(24)

where ${{\mathbf{d}}_{L}}={{d}_{L}}\mathbf{e}_{3}$, and ${{d}_{L}}$ is the
height of the LiDAR from the ground and it is obtained from (1).
${{\mathbf{d}}_{I}}=\hat{{{d}_{I}}}\mathbf{e}_{3}$, and $\hat{{{d}_{I}}}$
denotes the height of the IMU from the ground
${{(\mathbf{R}_{L}^{G})}^{T}}=\mathbf{R}_{G}^{L}$ and
${{(\mathbf{R}_{I}^{G})}^{T}}=\mathbf{R}_{G}^{I}$, and they are obtained from
(1) and (17), respectively. The proposed Ground Plane Motion (GPM) constraint
$\mathbf{r}_{g}$ is represented as

${{\mathbf{r}}_{g}}(\mathbf{R}_{L}^{I},\mathbf{p}_{L}^{I})=\left(\begin{matrix}{{[\mathbf{R}_{I}^{G}\mathbf{R}_{L}^{I}\mathbf{R}_{G}^{L}\mathbf{e}_{3}]}_{1,2}}\\\
\mathbf{e}_{3}^{T}{{(\mathbf{R}_{I}^{L}\mathbf{R}_{G}^{I}{{\mathbf{d}}_{I}}-\mathbf{R}_{G}^{L}{{\mathbf{d}}_{L}}-\mathbf{p}_{L}^{I})}}\\\
\end{matrix}\right)$ (25)

Unlike [7], which performs the calibration while ignoring unobservable
direction, our Gril-Calib conducts extrinsic calibration including
unobservable direction even under the planar motion by utilizing Ground Plane
Motion (GPM) constraints.

Figure 4: Illustration of the geometric rotation relationship between ground
and each sensor. This robot platform is used in the M2DGR dataset.

###  V-C Single Optimization (SO) Process

The method proposed in [20] uses two distinct constraints, $\mathbf{r}_{w}$
and $\mathbf{r}_{a}$, separately over the two-stage optimization process.
First, $\mathbf{R}^{I}_{L}$ is obtained by minimizing the residual
$\mathbf{r}_{w}$, and then $\mathbf{p}^{I}_{L}$ is determined by minimizing
the residual $\mathbf{r}_{a}$ in [20]. However, this approach overlooks the
interaction and correlation between these constraints. Specifically, the
residual $\mathbf{r}_{a}$ depends on not only $\mathbf{p}^{I}_{L}$ but also
$\mathbf{R}^{I}_{L}$. simultaneously. In our Gril-Calib, we address this
ineffectiveness by reformulating the calibration as a Single Optimization
(SO). Our optimization problem could be depicted as follows,

$\begin{gathered}\mathbf{x}^{*}_{cal}=\underset{\mathbf{x}}{\arg\min}\sum\\{\rho_{w}\lVert\mathbf{r}_{w}(\mathbf{R}^{I}_{L},\mathbf{b}_{g},\mathbf{t}_{f})\rVert^{2}_{\Sigma_{\mathbf{r}_{w}}}+\\\
\rho_{a}\lVert\mathbf{r}_{a}(\mathbf{R}^{I}_{L},\mathbf{p}^{I}_{L},\mathbf{b}_{a})\rVert^{2}_{\Sigma_{\mathbf{r}_{a}}}+\rho_{g}\lVert\mathbf{r}_{g}(\mathbf{R}^{I}_{L},\mathbf{p}^{I}_{L})\rVert^{2}_{\Sigma_{\mathbf{r}_{g}}}\\}\end{gathered}$
(26)

and minimize the two residuals $\mathbf{r}_{w}$ and $\mathbf{r}_{a}$ along
with our proposed residual $\mathbf{r}_{g}$ simultaneously, where
$\Sigma_{\mathbf{r}_{w}},\Sigma_{\mathbf{r}_{a}},\Sigma_{\mathbf{r}_{g}}$
denotes the Mahalanobis distance for each covariance matrix. $\rho_{x}$, where
$x=w,a,g$ refers to the scaled Cauchy kernel function assigned to each
residual.

##  VI Experiments

###  VI-A Experiments Setup

TABLE I: RMSE result compared to existing methods on M2DGR dataset

M2DGR [14] |  | Gate 01
---
[0,140]
| Rotation 02
---
[0, 120]
| Hall 04
---
[0, 160]
| Lift 04
---
[0, 80]
| Street 08
---
[20, 130]
Rotation (°) | Translation (m) | Rotation (°) | Translation (m) | Rotation (°) | Translation (m) | Rotation (°) | Translation (m) | Rotation (°) | Translation (m)
ICL [18] | 29.549 | 1.210 | 11.581 | 0.839 | 48.628 | 7.517 | 9.998 | 0.390 | 7.205 | 0.847
FAST-LIO2 [4] | 3.388 | 0.392 | 3.686 | 0.394 | 3.325 | 0.387 | 3.667 | 0.401 | 3.335 | 0.385
OA-LICalib [11] | 0.642 | 0.237 | 0.922 | 0.208 | 0.851 | 1.742 | 0.902 | 0.021 | 0.770 | 0.241
LI-Init [20] | 1.120 | 0.383 | 1.305 | 0.311 | 2.108 | 0.236 | 0.904 | 0.263 | 0.784 | 0.223
GRIL-Calib (Proposed) | 0.376 | 0.045 | 0.818 | 0.033 | 0.472 | 0.027 | 0.154 | 0.066 | 0.328 | 0.043

  * •

The best-performing algorithms are shown in bold.

Our implementation is based on C++, using ROS (Robot Operating System) on
Ubuntu 18.04. We perform our experiments using a 4-core Intel® Core™ i5-6600
CPU and Titan XP GPU. We use Ceres-solver 222http://ceres-solver.org/ to
compute the non-linear optimization. We evaluated the algorithms using three
public datasets: M2DGR [14], Hilti [15], and S3E [16], which cover both indoor
and outdoor environments. Following the evaluation scheme used in [20], we
only evaluate spatial extrinsic calibration parameters in terms of RMSE (Root
Mean Squared Error), skipping the evaluation on IMU LiDAR time offset
$t^{I}_{L}$. To show the effectiveness of our Gril-Calib, we compare our
method with other recent IMU-LiDAR calibration methods which are named ILC
[18], FAST-LIO2 [4], OA-LICalib [11], LI-Init [20]. All algorithms are
reproduced in the same environment.

###  VI-B Extrinsic Calibration Performance on the M2DGR

The calibration result on the M2DGR dataset is summarized in Table I. M2DGR
uses Velodyne VLP-32C as the LiDAR and Handsfree A9 as the IMU. For a
consistent comparison, the initial rotation and translation parameters are set
to be the same for all algorithms. The initial rotation parameters are set to
$\mathbf{R}^{I}_{L}=(-5.0,-5.0,5.0)$ with the unit in degrees. The initial
translation parameters deviated by 0.402 meters from the ground truth based on
RMSE, represented as $\mathbf{p}^{I}_{L}=(0.6,0.45,0.6)$. We initialized all
remaining state vectors $\mathbf{t}_{f},\mathbf{b}_{g},\mathbf{b}_{a}$ to
zeros. Data was collected from both indoor (Hall 04 and Lift 04) and outdoor
(Gate 01, Rotation 02, and Street 08) environments, with partial
visualizations provided in Fig 5. The calibration result on the M2DGR is
summarized in Table I. The table shows that our Gril-Calib significantly
outperforms the other competing methods, specifically reducing RMSE by 75% for
rotations and 79% for translations compared to the baseline [20].

Figure 5: Trajectory visualization for Street 08 acquired outdoors (Left), and
Hall 04 acquired indoors (Right). We evaluated the algorithm with data
collected from a variety of environments.

(a) IMU-LiDAR extrinsic rotation parameters results

(b) IMU-LiDAR extrinsic translation parameter results

Figure 6: Visualization for each of the extrinsic parameters. The dotted line
represents the ground truth. Closer to the dotted line represents the best
estimation.

The extrinsic parameters of each calibration method are visualized in Fig. 6.
Fig. 6(a) and Fig, 6(b) are the results of the rotation and translation
matrices, respectively. In both figures, the dotted lines denote the ground
truth of the extrinsic parameters. Thus, the closer to the dotted line, the
better the calibration method is. ICL [18] is excluded from the figure because
its error is significantly larger than that of the other methods.

###  VI-C Extrinsic Calibration Performance on the Others

S3E and HILTI data are used as additional datasets, and their results are
summarized in Table II and Table III, respectively. Both of the datasets are
captured indoors. The S3E dataset uses a Velodyne VLP-16 and Xsens MTi-30 IMU,
while the HILTI dataset uses an Ouster OS0-64 LiDAR and AlphaSense IMU. As
shown in Tables II and III, our algorithm outperforms other competing methods
as in M2DGR, and it shows the reliable accuracy of our Gril-Calib across
various robot platforms. We also observe that OA-LICalib fails in calibration
due to a memory computation issue in the Baseline 3 sequence in Table III. The
reason might be that OA-LICalib uses surfel maps for calibration and the
computational cost rapidly increases with the number of surfels. In contrast,
our algorithm works smoothly without any memory issues. The two tables along
with Table I show the effectiveness and versatility of our approach in various
real-world scenarios.

TABLE II: RMSE result compared to existing methods on S3E dataset S3E [16] |  |  Lab 01 $(\beta)$
---
[65, 210]
|  Lab 02 $(\beta)$
---
[220, 340]
Rotation (°) | Translation (m) | Rotation (°) | Translation (m)
ICL [18] | 6.430 | 0.288 | 9.817 | 0.813
FAST-LIO2 [4] | 3.489 | 0.303 | 3.351 | 0.300
OA-LICalib [11] | 0.339 | 0.194 | 18.512 | 0.321
LI-Init [20] | 1.238 | 0.116 | 1.657 | 0.066
GRIL-Calib (Proposed) | 0.935 | 0.036 | 0.847 | 0.027
TABLE III: RMSE result compared to existing methods on HILTI dataset

HILTI [15] |  | Basement 3
---
[85, 315]
| Basement 4
---
[75, 225]
Rotation (°) | Translation (m) | Rotation (°) | Translation (m)
ICL [18] | 13.308 | 4.555 | 44.294 | 4.052
FAST-LIO2 [4] | 1.964 | 0.196 | 1.892 | 0.188
OA-LICalib [11] | $\times$ | $\times$ | 2.368 | 0.185
LI-Init [20] | 0.374 | 0.092 | 2.584 | 0.194
GRIL-Calib (Proposed) | 0.716 | 0.018 | 0.589 | 0.047

  * •

$\times$ means that the algorithm failed.

###  VI-D Robustness to the Initial Guess

The extrinsic calibration usually requires initial parameter values, and it
significantly affects the performance [27]. We conduct some experiments while
varying initial guesses to see how robust the competing calibration methods
are to the initial guesses. We use Gate 01 sequence in M2DGR. The results are
shown in Fig 8. Fig 8(a) shows the change of the calibration error while
fixing the initial $\mathbf{p}_{L}^{I}$ but changing initial
$\mathbf{R}_{L}^{I}$. Fig 8(b) is the result when $\mathbf{R}_{L}^{I}$ is
fixed but $\mathbf{p}_{L}^{I}$ is varied. Interestingly, OA-LICalib showed
relatively good performance in the previous three datasets, but as the initial
guess increases, its accuracy rapidly degrades. Instead, LI-Init demonstrates
robust performance against the initial guess, but it is less accurate than our
Gril-Calib due to the limited planar motion.

Figure 7: Evaluation of LiDAR-inertial odometry using different calibration
algorithms.

(a) RMSE result of IMU-LiDAR extrinsic rotation parameters
$\mathbf{R}^{I}_{L}$

(b) RMSE result of IMU-LiDAR extrinsic translation parameters
$\mathbf{p}^{I}_{L}$

Figure 8: Calibration performance with changes initial guess. TABLE IV: Ablation study for proposed method | GPM | SO |  LO w/ GP | RMSE
---|---|---|---|---
| Rotation (°) | Translation (m)
|  LI-Init [20]
---
(Baseline)
|  |  | 1.710 | 0.300
GRIL-Calib | ✓ |  |  | 1.421 | 0.088
(Proposed) | ✓ | ✓ |  | 0.685 | 0.054
✓ | ✓ | ✓ | 0.430 | 0.041

###  VI-E Application to IMU-LiDAR Fusion System

We compare the competing extrinsic calibration methods by applying the
respective results to IMU-LiDAR fusion system and evaluating the corresponding
odometry accuracy. In the experiment, we utilize Faster-LIO [6] as the IMU-
LiDAR fusion system. The results in Table I used extrinsic parameters of each
method. The experiment uses the hall 02 sequence since the sequence was not
used in Subsection VI-B. The Absolute Pose Error (APE) is depicted from
different calibration methods in Fig 7. We observe that the odometry obtained
by our Gril-Calib outperforms those of other methods. In particular, our Gril-
Calib achieves 8.6% improvement from LI-Init and 9.6% improvement from OA-
LICalib. This experiment shows the importance of accurate calibration when
applied to IMU-LiDAR fusion systems.

Figure 9: Comparative visualization of the z-component of a trajectory using
Rotation 02 sequence in M2DGR.

###  VI-F Ablation Study

We ablate three key modules comparing the result with that of our baseline
[20]. The results are summarized in Table IV. On the first row of the table,
‘GPM’, ‘SO’, and ‘LO w/ GP’ denote the three contributions made in our paper.
That is, they imply ‘calibration using GPM constraints’, ‘calibration via
single optimization’, and ‘LO with proposed ground plane residual’,
respectively. Notably, SO observed a 34% decrease in RMSE for translation
calibration. This reduction highlights that $\mathbf{R}^{I}_{L}$ and
$\mathbf{p}^{I}_{L}$ are highly correlated and the separate optimization which
ignores the correlation would lead to degraded accuracy. Our LO method with
the proposed ground plane residual, effectively minimizes z-axis drift in the
LO, as shown in Fig 9. This drift reduction leads to more stable and accurate
measurements, thereby enhancing the overall precision of our calibration
approach.

##  VII CONCLUSIONS

In this paper, we present GRIL-Calib: a novel targetless IMU-LiDAR extrinsic
calibration framework that mainly focuses on the ground robot. Our framework
aimed to address the challenge of obtaining accurate 6-DoF extrinsic
parameters including unobservable direction when the robot’s motion is limited
to a planar motion. To achieve this, we proposed the Ground Plane Motion (GPM)
constraint to utilize ground observation, which is easily accessible for the
ground robots. We also reformulate the calibration as a Single Optimization
(SO) problem to take advantage of the correlation of the residuals.
Additionally, we proposed ground plane residual to reduce the z-drift of the
LiDAR odometry, thereby improving the accuracy of the calibration. Finally, we
demonstrated the effectiveness and reliability of our GRIL-Calib by applying
it to three public datasets and varying some initial calibration parameters.
We acknowledge that our algorithm’s performance might be affected by uneven
ground. Future works should focus on calibration for such terrains.

## References

  * [1] I. Vizzo, T. Guadagnino, B. Mersch, L. Wiesmann, J. Behley, and C. Stachniss, “Kiss-icp: In defense of point-to-point icp–simple, accurate, and robust registration if done the right way,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 2, pp. 1029–1036, 2023.
  * [2] J. Wang, M. Xu, G. Zhao, and Z. Chen, “3d lidar localization based on novel nonlinear optimization method for autonomous ground robot,” _IEEE Transactions on Industrial Electronics_ , 2023.
  * [3] J. Zhang and S. Singh, “Loam: Lidar odometry and mapping in real-time.” in _Robotics: Science and systems_ , vol. 2, no. 9. Berkeley, CA, 2014, pp. 1–9.
  * [4] W. Xu, Y. Cai, D. He, J. Lin, and F. Zhang, “Fast-lio2: Fast direct lidar-inertial odometry,” _IEEE Transactions on Robotics_ , vol. 38, no. 4, pp. 2053–2073, 2022.
  * [5] T. Shan, B. Englot, D. Meyers, W. Wang, C. Ratti, and D. Rus, “Lio-sam: Tightly-coupled lidar inertial odometry via smoothing and mapping,” in _2020 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2020, pp. 5135–5142.
  * [6] C. Bai, T. Xiao, Y. Chen, H. Wang, F. Zhang, and X. Gao, “Faster-lio: Lightweight tightly coupled lidar-inertial odometry using parallel sparse incremental voxels,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 4861–4868, 2022.
  * [7] S. Li, X. Li, Y. Zhou, and C. Xia, “Targetless extrinsic calibration of lidar–imu system using raw gnss observations for vehicle applications,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 72, pp. 1–11, 2023.
  * [8] W. Liu and Y. Li, “Error modeling and extrinsic–intrinsic calibration for lidar-imu system based on cone-cylinder features,” _Robotics and Autonomous Systems_ , vol. 114, pp. 124–133, 2019.
  * [9] Y. Yang, P. Geneva, K. Eckenhoff, and G. Huang, “Degenerate motion analysis for aided ins with online spatial and temporal sensor calibration,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 2070–2077, 2019.
  * [10] Y. Yang, P. Geneva, X. Zuo, and G. Huang, “Online self-calibration for visual-inertial navigation: Models, analysis, and degeneracy,” _IEEE Transactions on Robotics_ , 2023.
  * [11] J. Lv, X. Zuo, K. Hu, J. Xu, G. Huang, and Y. Liu, “Observability-aware intrinsic and extrinsic calibration of lidar-imu systems,” _IEEE Transactions on Robotics_ , vol. 38, no. 6, pp. 3734–3753, 2022.
  * [12] X. Wei, J. Lv, J. Sun, E. Dong, and S. Pu, “Gclo: Ground constrained lidar odometry with low-drifts for gps-denied indoor environments,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 2229–2235.
  * [13] H. Li and J. Stueckler, “Visual-inertial odometry with online calibration of velocity-control based kinematic motion models,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 6415–6422, 2022.
  * [14] J. Yin, A. Li, T. Li, W. Yu, and D. Zou, “M2dgr: A multi-sensor and multi-scenario slam dataset for ground robots,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 2266–2273, 2021.
  * [15] M. Helmberger, K. Morin, B. Berner, N. Kumar, G. Cioffi, and D. Scaramuzza, “The hilti slam challenge dataset,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 7518–7525, 2022.
  * [16] D. Feng, Y. Qi, S. Zhong, Z. Chen, Y. Jiao, Q. Chen, T. Jiang, and H. Chen, “S3e: A large-scale multimodal dataset for collaborative slam,” _arXiv preprint arXiv:2210.13723_ , 2022.
  * [17] C. Le Gentil, T. Vidal-Calleja, and S. Huang, “3d lidar-imu calibration based on upsampled preintegrated measurements for motion distortion correction,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 2149–2155.
  * [18] S. Mishra, G. Pandey, and S. 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)_. IEEE, 2021, pp. 1–7.
  * [19] J. Lv, J. Xu, K. Hu, Y. Liu, and X. Zuo, “Targetless calibration of lidar-imu system based on continuous-time batch estimation,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 9968–9975.
  * [20] F. Zhu, Y. Ren, and F. Zhang, “Robust real-time lidar-inertial initialization,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 3948–3955.
  * [21] S. Madgwick _et al._ , “An efficient orientation filter for inertial and inertial/magnetic sensor arrays,” _Report x-io and University of Bristol (UK)_ , vol. 25, pp. 113–118, 2010.
  * [22] S. Lee, H. Lim, and H. Myung, “Patchwork++: Fast and robust ground segmentation solving partial under-segmentation using 3d point cloud,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 13 276–13 283.
  * [23] C. Hertzberg, R. Wagner, U. Frese, and L. Schröder, “Integrating generic sensor fusion algorithms with sound state representations through encapsulation of manifolds,” _Information Fusion_ , vol. 14, no. 1, pp. 57–77, 2013.
  * [24] D. Baleanu, O. Defterli, and O. P. Agrawal, “A central difference numerical scheme for fractional optimal control problems,” _Journal of Vibration and Control_ , vol. 15, no. 4, pp. 583–597, 2009.
  * [25] F. Gustafsson, “Determining the initial states in forward-backward filtering,” _IEEE Transactions on signal processing_ , vol. 44, no. 4, pp. 988–992, 1996.
  * [26] W. Xu, D. He, Y. Cai, and F. Zhang, “Robots’ state estimation and observability analysis based on statistical motion models,” _IEEE Transactions on Control Systems Technology_ , vol. 30, no. 5, pp. 2030–2045, 2022.
  * [27] S. Xu, J. S. Willners, Z. Hong, K. Zhang, Y. R. Petillot, and S. Wang, “Observability-aware active extrinsic calibration of multiple sensors,” in _2023 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2023, pp. 2091–2097.