FAST-LIO: A Fast, Robust LiDAR-Inertial Odometry Package by Tightly-Coupled Iterated Kalman Filter

Wei Xu and Fu Zhang

Abstract—This letter presents a computationally efficient and robust LiDAR-inertial odometry framework. We fuse LiDAR feature points with IMU data using a tightly-coupled iterated extended Kalman filter to allow robust navigation in fast-motion, noisy or cluttered environments where degeneration occurs. To lower the computation load in the presence of a large number of measurements, we present a new formula to compute the Kalman gain. The new formula has computation load depending on the state dimension instead of the measurement dimension. The proposed method and its implementation are tested in various indoor and outdoor environments. In all tests, our method produces reliable navigation results in real-time: running on a quadrotor onboard computer, it fuses more than 1200 effective feature points in a scan and completes all iterations of an iEKF step within 25 ms. Our codes are open-sourced on Github.¹

Index Terms—Aerial systems, localization, perception and autonomy, sensor fusion.

I. INTRODUCTION

IMULTANEOUS localization and mapping (SLAM) is a fundamental prerequisite of mobile robots, such as unmanned aerial vehicles (UAVs). Visual (-inertial) odometry (VO), such as Stereo VO [1], [2] and Monocular VO [3], [4] are commonly used on mobile robots due to their lightweight and low-cost. Although providing rich RGB information, visual solutions lack direct depth measurements and require much computation resources to reconstruct the 3D environment for trajectory planning. Moreover, they are very sensitive to lighting conditions. Light detection and ranging (LiDAR) sensors could overcome all these difficulties but have been too costly (and bulky) for small-scale mobile robots.

Solid-state LiDARs recently emerge as main trends in LiDAR developments, such as those based on microelectro-mechanical-system (MEMS) scanning [5] and rotating prisms [6]. These LiDARs are very cost-effective (in a cost range similar to global shutter cameras), lightweight (can be carried by a small-scale UAV), and of high performance (producing active and direct 3D measurements of long-range and high-accuracy).

These features make such LiDARs viable for UAVs, especially industrial UAVs, which need to acquire accurate 3D maps of the environments (e.g., aerial mapping) or may operate in cluttered environments with severe illumination variations (e.g., post-disaster search and inspection).

Despite the great potentiality, solid-state LiDARs bring new challenges to SLAM: 1) the feature points in LiDAR measurements are usually the geometrical structures (e.g., edges and planes) in the environments. When the UAV is operating in cluttered environments where no strong features are present, the LiDAR-based solution easily degenerates. This problem is more obvious when the LiDAR has a small FoV. 2) Due to the high-resolution along the scanning direction, a LiDAR scan usually contains many feature points (e.g., a few thousand). While these feature points are not adequate to reliably determine the pose in case of degeneration, tightly fusing such a large number of feature points to IMU measurements requires tremendous computational resources that are not affordable by the UAV onboard computer. 3) Since the LiDAR samples point sequentially with a few laser/receiver pairs, laser points in a scan are always sampled at different times, resulting in motion distortion that will significantly degrade a scan registration [7]. The constant rotations of UAV propellers and motors also introduce significant noises to the IMU measurements.

To make the LiDAR navigation viable for small-scale mobile robots such as UAVs, we propose the FAST-LIO, a computationally efficient and robust LiDAR-inertial odometry package. More specifically, our contributions are as follows: 1) To cope with fast-motion, noisy or cluttered environments where degeneration occurs, we adopt a tightly-coupled iterated Kalman filter to fuse LiDAR feature points with IMU measurements. We propose a formal back-propagation process to compensate for the motion distortion; 2) To lower the computation load caused by a large number of LiDAR feature points, we propose a new formula for computing the Kalman gain and prove its equivalence to the conventional Kalman gain formula. The new formula has a computation complexity depending on the state dimension instead of the measurement dimension. 3) We implement our formulations into a fast and robust LiDAR-inertial odometry software package. The system is able to run on a small-scale quadrotor onboard computer. 4) We conduct experiments in various indoor and outdoor environments and with actual UAV flight tests (Fig. 1) to validate the system’s robustness when fast motion or intense vibration noise exists.

The remaining letter is organized as follows: In Section II, we discuss relevant research works. We give an overview of the complete system pipeline and the details of each key component in Section III. The experiments are presented in Section IV, followed by conclusions in Section V.

Fig. 1. Our LiDAR-inertial navigation system runs on a Livox AVIA LiDAR² and a DJI Manifold 2-C onboard computer³, all on a customized small-scale quadrotor UAV (280 mm wheelbase). The RGB camera is not used in our algorithm, but only for visualization. Video is available at https://youtu.be/iYCY6T79oNU.

II. RELATED WORKS

Existing works on LiDAR SLAM are extensive. Here we limit our review to the most relevant works: LiDAR-only odometry and mapping, loosely-coupled and tightly-coupled LiDAR-Inertial fusion methods.

A. LiDAR Odometry and Mapping

Besl et al. [7] propose an iterated closest points (ICP) method for scan registration, which builds the basis for LiDAR odometry. ICP performs well for dense 3D scans. However, for sparse point clouds of LiDAR measurements, the exact point matching required by ICP rarely exists. To cope with this problem, Segal et al. [8] propose a generalized-ICP based on the point-to-plane distance. Then Zhang et al. [9] combine this ICP method with a point-to-edge distance and developed a LiDAR odometry and mapping (LOAM) framework. Thereafter, many variants of LOAM have been developed, such as LeGO-LOAM [10] and LOAM-Livox [11]. While these methods work well for structured environments and LiDARs of large FoV, they are very vulnerable to featureless environments or small FoV Li-DARs [11].

B. Loosely-Coupled LiDAR-Inertial Odometry

IMU measurements are commonly used to mitigate the problem of LiDAR degeneration in featureless environments. Loosely-coupled LiDAR-inertial odometry (LIO) methods typically process the LiDAR and IMU measurements separately and fuse their results later. For example, IMU-aided LOAM [9] takes the pose integrated from IMU measurements as the initial estimate for LiDAR scan registration. Zhen et al. [12] fuse the IMU measurements and the Gaussian Particle Filter output of

LiDAR measurements using the error-state EKF. Balazadegan et al. [13] add the IMU-gravity model to estimate the 6-DOF ego-motion to aid the LiDAR scan registration. Zuo et al. [14] use a Multi-State Constraint Kalman Filter (MSCKF) to fuse the scan registration results with IMU and visual measurements. A common procedure of the loosely-coupled approach is obtaining a pose measurement by registering a new scan and then fusing the pose measurement with IMU measurements. The separation between scan registration and data fusion reduces the computation load. However, it ignores the correlation between the system’s other states (e.g., velocity) and the pose of the new scan. Moreover, in the case of featureless environments, the scan registration could degenerate in certain directions and causes unreliable fusion in later stages.

C. Tightly-Coupled LiDAR-Inertial Odometry

Unlike the loosely-coupled methods, tightly-coupled LiDARinertial odometry methods typically fuse the raw feature points (instead of scan registration results) of LiDAR with IMU data. There are two main approaches to tightly-coupled LIO: optimization-based and filter-based. Geneva et al. [15] use a graph optimization with IMU pre-integration constrains [16] and plane constraints [17] from LiDAR feature points. Recently, Ye et al. [18] propose the LIOM package which uses a similar graph optimization but is based on edge and plane features. For filter-based methods, Bry [19] uses a Gaussian Particle Filter (GPF) to fuse the data of IMU and a planar 2D LiDAR. This method has also been used in the Boston Dynamics Atlas humanoid robot. Since the computation complexity of particle filter grows quickly with the number of feature points and the system dimension, Kalman filter and its variants are usually more preferred, such as extended Kalman filter [20], unscented Kalman filter [21], and iterated Kalman filter [22].

Our method falls into the tightly-coupled approach. We adopt an iterated extended Kalman filter similar to [22] to mitigate linearization errors. Kalman filter (and its variants) has a time complexity $O (m^{2})$ where m is the measurement dimension [23], this may lead to remarkably high computation load when dealing with a large number of LiDAR measurements. Naive downsampling would reduce the number of measurements but at the cost of information loss. [22] reduces the number of measurements by extracting and fitting ground planes similar to [10]. This, however, does not apply to aerial applications where the ground plane may not always present.

III. METHODOLOGY

A. Framework Overview

This letter will use the notations shown in Table I. The overview of our workflow is shown in Fig. 2. The LiDAR inputs are fed into the feature extraction module to obtain planar and edge features. Then the extracted features and IMU measurements are fed into our state estimation module for state estimation at 10 Hz − 50 Hz. The estimated pose then registers the feature points to the global frame and merges them with the feature points map built so far. The updated map is finally used to register further new points in the next step.

Fig. 2. System overview of FAST-LIO. (a): the overall pipeline; (b): the forward and backward propagation.

B. System Description

1) - / Operator: Let M be the manifold of dimension n in consideration $(e . g ., M = S O (3))$ . Since manifolds are locally homeomorphic to $R^{n}$ , we can establish a bijective mapping from a local neighborhood on M to its tangent space $R^{n}$ via two encapsulation operators - and [24]:

\begin{matrix} ⊞ : M \times R^{n} \to M; \overset{⊖ d \dag, M \times M \to R^{n}}{⊔ \dag, M} \\ M = S O (3) : R ⊞ r = R E x p (r); R_{1} ⊞ R_{2} = L o g (R_{2}^{T} R_{1}) \\ M = R^{n} : a ⊞ b = a + b; a ⊞ b = a - b \end{matrix}

where $\begin{array} { r } { \mathrm { E x p } ( \mathbf { r } ) { = } \mathbf { I } + \frac { \mathbf { r } } { \Vert \mathbf { r } \Vert } \sin ( \Vert \mathbf { r } \Vert ) { + } \frac { \mathbf { r } ^ { 2 } } { \Vert \mathbf { r } \Vert ^ { 2 } } ( 1 - \cos ( \Vert \mathbf { r } \Vert ) ) } \end{array} $i s t h e e x p o n e n t i a l m a p [24] a n d$ \mathrm { L o g ( \cdot ) } $i s i t s i n v e r s e m a p . F o r a c o m p o u n d m a n i f o l d$ \mathcal { M } = S O ( 3 ) \times \mathbb { R } ^ { n }$ we have:

Π_{a}^{\tilde{A}} \equiv [r_{b}^{\tilde{r}}] = [R \oplus r_{\tilde{c}}]; [R_{1}] \equiv [R_{2}] = [R_{1} ⊖ R_{2}]

From the above definition, it is easy to verify that

(x ⊞ u) ⊞ x = u; x ⊞ (y ⊞ x) = y; \forall x, y \in M, \forall u \in R^{n} .

2) Continuous Model: Assuming an IMU is rigidly attached to the LiDAR with a known extrinsic $^{I} T_{L} = {\bar{(}}^{I} {\bar{R}}_{L},^{I} p_{L})$ Taking the IMU frame (denoted as I) as the body frame of reference leads to a kinematic model:

\begin{matrix} Λ^{G} {\dot{p}}_{I} = Λ^{G} v_{I}, Λ^{G} {\dot{v}}_{I} = Λ^{G} R_{I} (a_{m} - b_{a} - n_{a}) + Λ^{G} g, Λ^{G} \dot{g} = 0 \\ Λ^{G} {\dot{R}}_{I} = Λ^{G} R_{I} {ω ω - b_{ω} - n_{ω}}_{Λ}, {\dot{b}}_{ω} = n_{b ω}, {\dot{b}}_{a} = n_{b a} \end{matrix}

where $Λ^{G} p_{I}, Λ^{G} R_{I}$ are the position and attitude of IMU in the global frame $(i . e .$ , the first IMU frame, denoted as $G), ω^{G} g$ is the unknown gravity vector in the global frame, $a_{m}$ and $ω_{m}$ are IMU measurements, $n_{a}$ and $n_{ω}$ are the white noise of IMU measurements, $b_{a}$ and $b_{ω}$ are the IMU bias modelled as the random walk process with Gaussian noises $n_{b a}$ and $n_{b ω},$ and the notation $⌊ a ⌋_{\land}$ denotes the skew-symmetric matrix of vector $a \in R^{3}$ that maps the cross product operation.

3) Discrete Model: Based on the - operation defined above, we can discretize the continuous model in (1) at the IMU sampling period $Δ t$ using a zero-order holder. The resultant discrete model is

x_{i + 1} = x_{i} ⊞ (Δ t f (x_{i}, u_{i}, w_{i}))

where i is the index of IMU measurements, the function $f,$ state $x,$ input u and noise w are defined as below:

\begin{array} { r l } & { M = S O ( 3 ) \times { \mathbf B } ^ { 1 5 } , \ \mathrm { d } \mathrm { m } (  { \mathcal { M } } ) = 1 8 } \\ & { \quad \times = [ { \mathbf R } _ { T } ^ { T } \quad G _ { { \mathbf B } } ^ { T } \quad { \mathbf G } _ { T } ^ { T } \quad { \mathbf b } _ { T } ^ { T } \quad { \mathbf b } _ { \omega } ^ { T } \quad { \mathbf b } _ { R } ^ { T } \quad G _ { { \mathbf B } } ^ { T } \quad { \mathbf G } _ { { \mathbf B } } ^ { T } ] ^ { T } \in \mathcal { M } } \\ & { \quad = [ \omega _ { m } ^ { T } \quad { \mathbf a } _ { m } ^ { T } ] ^ { T } , \ w \stackrel { \cdot } { = } [ { \mathbf n } _ { \omega } ^ { T } \quad \mathbf { n } _ { \omega } ^ { T } \quad \mathbf { n } _ { \mathbf { b } \omega } ^ { T } \quad \mathbf { n } _ { \mathbf { b } \omega } ^ { T } \quad \mathbf { n } _ { \mathbf { b } \omega } ^ { T } ] ^ { T } } \\ & { \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad } \\ & { \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad } \\ &  \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \mathrm { h } _ { \mathbf { b } \omega _ { i } } \quad \mathbf { n } _ { \mathbf { b } \omega } ^ { T } ) = (  \mathbf { G } _ { { \mathbf B } _ { i } } \quad ( \mathbf { n } _ { m _ { i } } - \mathbf { b } _ { \mathbf { a } _ { m } _ { i } } ) \quad + { \mathbf G } _ { { \mathbf B } _ { i } } | \quad \quad \quad \quad \quad \quad \end{array}\tag{3}

4) Preprocessing of LiDAR Measurements: LiDAR measurements are point coordinates in its local body frame. Since the raw LiDAR points are sampled at a very high rate (e.g., 200 kHz), it is usually not possible to process each new point once being received. A more practical approach is to accumulate these points for a certain time and process them all at once. In FAST-LIO, the minimum accumulation interval is set to 20 ms, leading to up to 50 $H z$ full state estimation (i.e., odometry output) and map update as shown in Fig. 2 (a). Such an accumulated set of points is called a scan, and the time for processing it is denoted as $t_{k}$ (see Fig. 2 (b)). From the raw points, we extract planar points with high local smoothness [9] and edge points with low local smoothness as in [11]. Assume the number of feature points is $m .$ , each is sampled at time $ρ_{j} \in (t_{k - 1}, t_{k}]$ and is denoted as $L_{j} {σ_{p}}_{f_{j}}$ , where $L_{j}$ is the LiDAR local frame at the time $ρ_{j}$ . During a LiDAR scan, there are also multiple IMU measurements, each sampled at time $τ_{i} \in [t_{k - 1}, t_{k}]$ with the respective state $x_{i}$ as in (2). Notice that the last LiDAR feature point is the end of a scan, i.e., $ρ_{m} = t_{k}$ , while the IMU measurements may not necessarily be aligned with the start or end of the scan.

C. State Estimation

To estimate the states in the state formulation (2), we use an iterated extended Kalman filter. Moreover, we characterize the estimation covariance in the tangent space of the state estimate as in [24], [25]. Assume the optimal state estimate of the last LiDAR scan at $t_{k - 1} i s {\bar{x}}_{k - 1}$ with covariance matrix ${\bar{P}}_{k - 1}$ . Then ${\bar{P}}_{k - 1}$ represents the covariance of the random error state vector defined below:

{\tilde{x}}_{k - 1} ≐ x_{k - 1} ⊡ u {\bar{x}}_{k - 1} = {[δ θ^{T} Σ^{G} {\tilde{p}}_{I}^{T} Σ^{G} {\tilde{v}}_{I}^{T} {\tilde{b}}_{ω}^{T} {\tilde{b}}_{a}^{T} Σ^{G} {\tilde{g}}^{T}]}^{T}

where $δ θ = L o g (G {\bar{R}}_{I}^{T G} R_{I})$ is the attitude error and the rests are standard additive errors (i.e., the error in the estimate x¯ of a quantity x is $\tilde{x} = x - \bar{x})$ . Intuitively, the attitude error $δ θ$ describes the (small) deviation between the true and the estimated attitude. The main advantage of this error definition is that it allows us to represent the attitude uncertainty by the $3 \times 3$ covariance matrix $E {δ θ δ θ^{T}}$ . Since the attitude has $3^{\circ}$ of freedom (DOF), this is a minimal representation.

1) Forward Propagation: The forward propagation is performed once receiving an IMU input (see Fig. 2). More specifically, the state is propagated following (2) by setting the process noise $w_{i}$ to zero:

{\hat{x}}_{i + 1} = {\hat{x}}_{i} ⊞ (Δ t f ({\hat{x}}_{i}, u_{i}, 0)); {\hat{x}}_{0} = {\bar{x}}_{k - 1} .

To propagate the covariance, we use the error state dynamic model obtained below:

\begin{matrix} {\tilde{x}}_{i + 1} = x_{i + 1} ⊞ {\hat{x}}_{i + 1} \\ = (x_{i} ⊞ Δ t f (x_{i}, u_{i}, w_{i})) ⊞ ({\hat{x}}_{i} ⊞ Δ t f ({\hat{x}}_{i}, u_{i}, 0)) \\ \overset{(23)}{≃} F_{\tilde{x}} {\tilde{x}}_{i} + F_{w} w_{i} . \end{matrix}

The matrix $F_{\tilde{x}}$ and $F_{w}$ in (5) is computed following the Appendix. V-A. The result is shown in (7) shown at the bottom of this page, where $\widehat { \pmb { \omega } } _ { i } = \pmb { \omega } _ { m _ { i } } - \widehat { \mathbf { b } } _ { \omega _ { i } } , \widehat { \mathbf { a } } _ { i } = \mathbf { a } _ { m _ { i } } - \widehat { \mathbf { b } } _ { \mathbf { a } } $a n d$ { \mathbf { A } ( \mathbf { u } ) ^ { - 1 } }$ follows the same definition in [26] as below:

\begin{matrix} {A (u)}^{- 1} = I - \frac{1}{2} ⌊ u ⌋_{\land} + (1 - α (‖ u ‖)) \frac{⌊ u ⌋_{\land}^{2}}{‖ u ‖^{2}} \\ α (m) = \frac{m}{2} c o t (\frac{m}{2}) = \frac{m}{2} \frac{c o s (m / 2)}{s i n (m / 2)} \end{matrix}

Denoting the covariance of white noises w as Q, then the propagated covariance ${\hat{P}}_{i}$ can be computed iteratively following the below equation.

\begin{matrix} {\hat{P}}_{i + 1} = F_{\tilde{x}} {\hat{P}}_{i} F_{\tilde{x}}^{T} + F_{w} Q F_{w}^{T}; {\hat{P}}_{0} = {\bar{P}}_{k - 1} . \end{matrix}

The propagation continues until reaching the end time of a new scan at $t_{k}$ where the propagated state and covariance are denoted as ${\hat{x}}_{k}, {\hat{P}}_{k}$ . Then ${\hat{P}}_{k}$ represents the covariance of the error between the ground-truth state $x_{k}$ and the state propagation ${\hat{x}}_{k} (i . e ., x_{k} \in {\hat{x}}_{k})$

2) Backward Propagation and Motion Compensation: When the points accumulation time interval is reached at time $t_{k}$ , the new scan of feature points should be fused with the propagated state ${\hat{x}}_{k}$ and covariance ${\hat{P}}_{k}$ to produce an optimal state update. However, although the new scan is at time $t_{k}$ , the feature points are measured at their respective sampling time $ρ_{j} \leq t_{k}$ (see Section. III-B4 and Fig. 2 (b)), causing a mismatch in the body frame of reference.

To compensate the relative motion (i.e., motion distortion) between time $ρ_{j}$ and time $t_{k},$ , we propagate (2) backward as $\check { \mathbf { x } } _ { j - 1 } = \check { \mathbf { x } } _ { j } \boxplus ( - \Delta t \mathbf { f } ( \check { \mathbf { x } } _ { j } , \mathbf { u } _ { j } , \mathbf { 0 } ) ) $), s t a r t i n g f r o m z e r o p o s e a n d r e s t s t a t e s (e . g ., v e l o c i t y a n d b i a s) f r o m$ \widehat { \mathbf { x } } _ { k }$ . The backward propagation is performed at the frequency of feature point, which is usually much higher than the IMU rate. For all the feature points sampled between two IMU measurements, we use the left IMU measurement as the input in the back propagation. Furthermore, noticing that the last three block elements (corresponding to the gyro bias, accelerometer bias, and gravity) of $f (x_{j}, u_{j}, 0)$ (see (3)) are zeros, the back propagation can be reduced to:

\begin{matrix} I_{k} {\overset{ˇ}{p}}_{I_{j - 1}} =^{I_{k}} {\overset{ˇ}{p}}_{I_{j}} -^{I_{k}} {\overset{ˇ}{v}}_{I_{j}} Δ t, {s . f .}^{I_{k}} {\overset{ˇ}{p}}_{I_{m}} = 0; \\ I_{k} {\overset{ˇ}{v}}_{I_{j - 1}} =^{I_{k}} {\overset{ˇ}{v}}_{I_{j}} -^{I_{k}} {\overset{ˇ}{R}}_{I_{j}} (a_{m_{i - 1}} - {\hat{b}}_{a_{k}}) Δ t -^{I_{k}} {\hat{g}}_{k} Δ t, \\ {s . f .}^{I_{k}} {\overset{ˇ}{v}}_{I_{m}} =^{G} {\hat{R}}_{I_{k}}^{T}^{G} {\hat{v}}_{I_{k}},^{I_{k}} {\hat{g}}_{k} =^{G} {\hat{R}}_{I_{k}}^{T} {\hat{g}}_{k}; \\ I_{k} {\overset{ˇ}{R}}_{I_{j - 1}} =^{I_{k}} {\overset{ˇ}{R}}_{I_{j}} E x p (({\hat{b}}_{ω_{k}} - ω_{m_{i - 1}}) Δ t), {s . f .}^{I_{k}} R_{I_{m}} = I, \end{matrix}

where $ρ_{j - 1} \in [τ_{i - 1}, τ_{i})$ and s.f. means “starting from”.

The backward propagation will produce a relative pose between time $ρ_{j}$ and the scan-end time $\begin{matrix} t_{k} : & ^{I_{k}} {\overset{ˇ}{T}}_{I_{j}} = \end{matrix}$ $(^{I_{k}} {\overset{ˇ}{R}}_{I_{j}},^{I_{k}} {\overset{ˇ}{p}}_{I_{j}})$ . This relative pose enables us to project the local measurement ${L_{j}}_{p_{f_{j}}}$ to scan-end measurement ${\bar{L_{k}}}_{\bar{p_{f_{j}}}}$ as follows (see Fig. 2):

^{L_{k}} p_{f_{j}} =^{I} T_{L}^{- 1 I_{k}} {\overset{ˇ}{T}}_{I_{j}}^{I} T_{L}^{L_{j}} p_{f_{j}},

where $I_{T_{L}}$ is the known extrinsic (see Section. III-B2). Then the projected point $\boldsymbol { L } _ { k } _ { \mathbf { p } _ { f _ { j } } }$ is used to construct a residual in the following section.

3) Residual Computation: With the motion compensation in (10), we can view the scan of feature points ${{^{L_{k}}}_{p_{f_{j}}}}$ all sampled at the same time $t_{k}$ and use it to construct the residual. Assume the current iteration of the iterated Kalman filter is $κ,$ and the corresponding state estimate is ${\hat{x}}_{k}^{κ}$ . When $κ = 0, {\hat{x}}_{k}^{κ} = {\hat{x}}_{k}$ , the

\mathbf { F } _ { \widetilde { \mathbf { x } } } = [ \begin{array} { c c c c c c c } { \mathrm { E x p } ( - \widetilde { \omega } _ { i } \Delta t ) } & { \mathbf { 0 } } & { \mathbf { 0 } } & { - \mathbf { A } ( \widetilde { \omega } _ { i } \Delta t ) ^ { T } \Delta t } & { \mathbf { 0 } } & { \mathbf { 0 } } \\ { \mathbf { 0 } } & { \mathbf { I } \ \mathbf { L } \Delta t } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } \\ { - ^ { G } \widehat { \mathbf { R } } _ { I _ { i } }  \widehat { \mathbf { a } } _ { i }  _ { , N } \Delta t } & { \mathbf { 0 } } & { \mathbf { I } } & { \mathbf { 0 } } & { - ^ { G } \widehat { \mathbf { R } } _ { I _ { i } } \Delta t \mathbf { \Delta } \mathbf { I } \Delta t } \\ { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { I } } & { \mathbf { 0 } } & { \mathbf { 0 } } \\ { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { I } } & { \mathbf { 0 } } \\ { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { I } } \end{array} ] , \mathbf { F } _ { \mathbf { w } } = [ \begin{array} { c c c c c c } { - \mathbf { A } ( \widetilde { \omega } _ { i } \Delta t ) ^ { T } \Delta t } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } \\ { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { 0 } } \\ { \mathbf { 0 } } & { - ^ { G } \widehat { \mathbf { R } } _ { I _ { i } } \Delta t } & { \mathbf { 0 } } & { \mathbf { 0 } } \\ { \mathbf { 0 } } & { \mathbf { 0 } } & { \mathbf { L } \Delta t } &  \mathbf { 0 }  \end{array}\tag{7}

predicted state from the propagation in (4). Then, the feature points ${{^{L_{k}}}_{p_{f_{j}}}}$ can be transformed to the global frame as below:

^{G} {\hat{p}}_{f_{j}}^{κ} =^{G} {\hat{T}}_{I_{k}}^{κ}^{I} T_{L}^{L_{k}} p_{f_{j}}; j = 1, \dots, m .

For each LiDAR feature point, the closest plane or edge defined by its nearby feature points in the map is assumed to be where the point truly belongs to. That is, the residual is defined as the distance between the feature point’s estimated global frame coordinate $G_{{\hat{p}}_{f_{j}}} κ_{}$ and the nearest plane (or edge) in the map. Denoting $u_{j}$ the normal vector (or edge orientation) of the corresponding plane (or edge), on which lying a point $^{G} q_{j}$ then the residual $z_{j}^{κ}$ is computed as:

z_{j}^{κ} = G_{j} (G_{{\hat{p}}_{f_{j}}^{κ}} - G_{q_{j}})

where $G_{j} = u_{j}^{T}$ for planar features and $G_{j} = ⌊ u_{j} ⌋_{\land}$ for edge features. The computation of the $u_{j}$ and the search of nearby points in the map, which define the corresponding plane or edge, is achieved by building a KD-tree of the points in the most recent map [11]. Moreover, we only consider residuals whose norm is below certain threshold (e.g., 0.5m). Residuals exceeding this threshold are either outliers or newly observed points.

4) Iterated State Update: To fuse the residual $z_{j}^{κ}$ computed in (12) with the state prediction ${\hat{x}}_{k}$ and covariance ${\hat{P}}_{k}$ propagated from the IMU data, we need to linearize the measurement model that relates the residual $z_{j}^{κ}$ to the ground-truth state $x_{k}$ and measurement noise. The measurement noise originates from the LiDAR ranging and beam-directing noise $L_{^{j} n_{f_{j}}}$ when measuring the point ${\bar{L}}_{j} {\bar{p}}_{f_{j}}$ . Removing this noise from the point measurement $L_{j} σ_{p_{f_{j}}}$ leads to the true point location

{\sp L}_{j} p_{f_{j}}^{g t} = {\sp L}_{j} p_{f_{j}} - {\sp L}_{j} n_{f_{j}} .

This true point, after projecting to the frame $L_{k}$ via (10) and then to the global frame with the ground-truth state $x_{k}$ (i.e, pose), should lie exactly on the plane (or edge) in the map. That is, plugging (13) into (10), then into (11), and further into (12) should result in zero. i.e.,

0 = h_{j} (x_{k}, Λ_{}^{L_{j}} n_{f_{j}}) = G_{j} (^{G} T_{I_{k}}^{I_{k}} {\hat{T}}_{I_{j}}^{I} T_{L} (^{L_{j}} p_{f_{j}} -^{L_{j}} n_{f_{j}}) -^{G} q_{j})

Approximating the above equation by its first order approximation made at ${\hat{x}}_{k}^{κ}$ leads to

\begin{matrix} 0 = h_{j} (x_{k}, ξ^{L_{j}} n_{f_{j}}) ≃ h_{j} ({\hat{x}}_{k}^{κ}, 0) + H_{j}^{κ} {\tilde{x}}_{k}^{κ} + v_{j} \\ = z_{j}^{κ} + H_{j}^{κ} {\tilde{x}}_{k}^{κ} + v_{j} \end{matrix}

where ${\tilde{x}}_{k}^{κ} = x_{k} ⊟ {\hat{x}}_{k}^{κ}$ (or equivalently $x_{k} = {\hat{x}}_{k}^{κ} E {\tilde{x}}_{k}^{κ})$ , H^κ is the Jacobin matrix of $\mathbf { h } _ { j } ( \widehat { \mathbf { x } } _ { k } ^ { \kappa } \boxed { \mathbb { H } } \widetilde { \mathbf { x } } _ { k } ^ { \kappa } , { \cal L } _ { j } \mathbf { n } _ { f _ { j } } ) $w i t h r e s p e c t t o$ \widetilde { \mathbf { x } } _ { k } ^ { \kappa } $e v a l u a t e d a t z e r o, a n d$ \mathbf { v } _ { j } \in \mathcal { N } ( \mathbf { 0 } , \mathbf { R } _ { j } ) $c o m e s f r o m t h e r a w m e a s u r e m e n t n o i s e$ L _ { j } { \bf \Pi } _ { \mathbf { n } _ { f _ { i } } }$

Notice that the prior distribution of $x_{k}$ obtained from the forward propagation in Section. III-C1 is for

\begin{matrix} x_{k} ⊟ {\hat{x}}_{k} = ({\hat{x}}_{k}^{κ} ⊞ {\tilde{x}}_{k}^{κ}) ⊞ {\hat{x}}_{k} = {\hat{x}}_{k}^{κ} ⊞ {\hat{x}}_{k} + J^{κ} {\tilde{x}}_{k}^{κ} \end{matrix}

where $J^{κ}$ is the partial differentiation of $\left( \widehat { \mathbf { x } } _ { k } ^ { \kappa } \boxplus \widetilde { \mathbf { x } } _ { k } ^ { \kappa } \right) \boxminus \widehat { \mathbf { x } } _ { k } $w i t h r e s p e c t t o$ \widetilde { \mathbf { x } } _ { k } ^ { \kappa }$ evaluated at zero:

J^{κ} = [\begin{matrix} A {(^{G} {\hat{R}}_{I_{k}}^{κ} ⊞^{G} {\hat{R}}_{I_{k}})}^{- T} & 0_{3 \times 15} \\ 0_{15 \times 3} & I_{15 \times 15} \end{matrix}]

where $A (\cdot)^{- 1}$ is defined in (6). For the first iteration (i.e., the case of extended Kalman filter), ${\hat{x}}_{k}^{κ} = {\hat{x}}_{k}$ , then $J^{κ} = I$

Combining the prior in (15) with the posteriori distribution from (14) yields the maximum a-posteriori estimate (MAP):

{min}_{{\tilde{x}}_{k}^{κ}} (‖ x_{k} ⊞ {\hat{x}}_{k} ‖_{{\hat{P}}_{k}^{- 1}}^{2} + \sum_{j = 1}^{m} ‖ z_{j}^{κ} + H_{j}^{κ} {\tilde{x}}_{k}^{κ} ‖_{R_{j}^{- 1}}^{2})

where $‖ x ‖_{M}^{2} = x^{T} M x$ . Substituting the linearization of the prior in (15) into (17) and optimizing the resultant quadratic cost leads to the standard iterated Kalman filter [22], which can be computed below (to simplify the notation, let $H = [H_{1}^{κ^{T}}, \dots, H_{m}^{{\bar{κ}}^{T}}]^{T}$ $R = d i a g (R_{1}, \cdot \cdot \cdot R_{m})$ , P₌ $(J^{κ})^{- 1} {\hat{P}}_{k} (J^{κ})^{- T}$ , and $z_{k}^{κ} = [z_{1}^{κ^{T}}, \dots, z_{m}^{κ^{T}}]^{T})$

\begin{matrix} K = {P H}^{T} (H {P H}^{T} + R)^{- 1}, \\ {\hat{x}}_{k}^{κ + 1} = {\hat{x}}_{k}^{κ} H (- K z_{k}^{κ} - (I - K H) (J^{κ})^{- 1} ({\hat{x}}_{k}^{κ} ⊞ {\hat{x}}_{k})) . \end{matrix}

The updated estimate ${\hat{x}}_{k}^{κ + 1}$ is then used to compute the residual in Section. III-C3 and repeat the process until convergence (i.e., $‖ {\hat{x}}_{k}^{κ + 1} ⊖ {\hat{x}}_{k}^{κ} ‖ < ϵ)$ . After convergence, the optimal state estimation and covariance is:

{\bar{x}}_{k} = {\hat{x}}_{k}^{κ + 1}, {\bar{P}}_{k} = (I - K H) P

A problem with the commonly used Kalman gain form in (18) is that it requires to invert the matrix ${H P H}^{T} \bar{+}$ R which is in the dimension of the measurements. In practice, the number of LiDAR feature points are very large in number, inverting a matrix of this size is prohibitive. As such, existing works [22], [27] only use a small number of measurements. In this letter, we show that this limitation can be avoided. The intuition originates from (17) where the cost function is over the state, hence the solution should be calculated with complexity depending on the state dimension. In fact, if directly solving (17), we can obtain the same solution in (18) but with a new form of Kalman gain shown below:

K = (H^{T} R^{- 1} H + P^{- 1})^{- 1} H^{T} R^{- 1} .

We prove in Appendix B that the two forms of Kalman gains are indeed equivalent based on the matrix inverse lemma [28]. Since the LiDAR measurements are independent, the covariance matrix R is (block) diagonal and hence the new formula only requires to invert two matrices both in the dimension of state instead of measurements. The new formula greatly saves the computation as the state dimension is usually much lower than measurements in LIO (e.g., more than 1000 effective feature points in a scan for 10 Hz scan rate while the state dimension is only 18).

Algorithm 1: State Estimation.
Input:Last optimal estimation ${\bar{x}}_{k - 1}$ and ${\bar{P}}_{k - 1},$
IMU inputs $(a_{m}, ω_{m})$ in current scan;
LiDAR feature points ${L_{j}}_{p_{f_{j}}}$ in current scan. Forward
propagation to obtain state prediction ${\hat{x}}_{k}$ via (4) and
covariance prediction ${\hat{P}}_{k}$ via (8);
Backward propagation to obtain ${ L } _ { k } _ { \mathbf { p } _ { f _ { j } } }$ via (9), (10);
$κ = - 1, {\hat{x}}_{k}^{κ = 0} = {\hat{x}}_{k};$
repeat $| \widehat { \mathbf { x } } _ { k } ^ { \kappa + 1 } \boxdot { \mathbf { x } } \widehat { \mathbf { x } } _ { k } ^ { \kappa } | < \epsilon \kappa = \kappa + 1 ;$
Compute $J^{κ}$ via (16) and $P = (J^{κ})^{- 1} {\hat{P}}_{k} (J^{κ})^{- T};$
Compute residual $z_{j}^{κ}$ (12) and Jocobin H^κ (14);
Compute the state update ${\hat{x}}_{k}^{κ + 1}$ via (18) with the Kalman
gain K from $(20);$
$\bar { \mathbf { x } } _ { k } = \widehat { \mathbf { x } } _ { k } ^ { \kappa + 1 } ; \bar { \mathbf { P } } _ { k } = ( \mathbf { I } - \mathbf { K H } ) \mathbf { P } .$
Output:Current optimal estimation ${\bar{x}}_{k}$ and ${\bar{P}}_{k} .$

5) The Algorithm: Our state estimation is summarized in Algorithm 1.

D. Map Update

With the state update ${\bar{x}}_{k}$ (hence ${ } ^ { G } \bar { \mathbf { T } } _ { I _ { k } } = ( { } ^ { G } \bar { \mathbf { R } } _ { I _ { k } } , { } ^ { G } \bar { \mathbf { p } } _ { I _ { k } } ) ) $e a c h f e a t u r e p o i n t$ ( ^ { L _ { k } } \mathbf { p } _ { f _ { j } } ) $p r o j e c t e d t o t h e b o d y f r a m e$ L _ { k }$ (see (10)) is then transformed to the global frame via:

^{G} {\bar{p}}_{f_{j}} =^{G} {\bar{T}}_{I_{k}}^{I} T_{L}^{L_{k}} p_{f_{j}}; j = 1, \dots, m .

These features points are finally appended to the existing map containing feature points from all previous steps.

E. Initialization

To obtain a good initial estimate of the system state (e.g., gravity vector ${\overset{˘}{G}}_{g}$ , bias, and noise covariance) so to speedup the state estimator, initialization is required. In FAST-LIO, the initialization is simple: keeping the LiDAR static for several seconds (2 seconds for all the experiments in this letter), the collected data is then used to initialize the IMU bias and the gravity vector. If non-repetitive scanning is supported by the LiDAR (e.g., Livox AVIA), keeping static also allows the LiDAR to capture an initial high-resolution map that is beneficial for the subsequent navigation.

IV. EXPERIMENT RESULTS

A. Computational Complexity Experiments

In order to validate the computational efficiency of the proposed new formula for computing Kalman gains. We intentionally replace the computation of Kalman gains with the old formula in our system and compare their computation time under the same system pipeline and number of feature points. The results are shown in Table. II. It is obvious that the complexity of the new formula is much lower than the old one.

TABLE I
SOME IMPORTANT NOTATIONS

Symbols	Meaning
$t_{k}$	The scan-end time of the k-th LiDAR scan.
$τ_{i}$	The i-th IMU sample time in a LiDAR scan.
$ρ_{j}$	The j-th feature point's sample time in a LiDAR scan.
$\dot{I_{i}}, I_{j}, I_{k}$	The IMU body frame at the time $τ_{i}, ρ_{j}$ and $t_{k} .$
$L_{j}, {\dot{L}}_{k}$	The LiDAR body frame at the time $ρ_{j}$ and $t_{k} .$
$x, \hat{x}, \bar{x}$	The ground-true, propagated, and updated value of x. $\overset{―}{x} .$
$\tilde{x}$	The error between ground-true x and its estimation
${\hat{x}}^{κ}$	The κ-th update of x in the iterated Kalman filter.
$x_{i}, x_{j}, x_{k}$	The vector (e.g.,state) x at time $τ_{i}, ρ_{j}$ and $t_{k} .$
${\overset{ˇ}{x}}_{j}$	Estimate of $x_{j}$ relative to $x_{k}$ in the back propagation.

TABLE II

THE RUNNING TIMES OF TWO KALMAN GAIN FORMULAS

Feature Num.	307	717	998	1243	1453	1802
Old Formula (ms)	7.1	23.4	109.3	251	1219	1621
New Formula (ms)	0.07	0.11	0.25	0.37	0.59	1.16

Fig. 3. During the flight experiment, the UAV is automatically flying in a circle path with 1.8 m radius and 1.4 m height. The circle path is conducted repeatedly for 4 times with different periods (6-10 s). The yaw command of the UAV maintains constant during the flight. In the end, the UAV is manually controlled to land at the take-off point, which enables us to measure the drift.

B. UAV Flight Experiments

In order to validate the robustness and computational efficiency of FAST-LIO in actual mobile robots, we build a smallscale quadrotor that can carry a Livox Avia LiDAR with 70^◦ FoV and a DJI Manifold 2-C onboard computer with a 1.8 GHz quad-core Intel i7-8550 U CPU and 8 GB RAM, as shown in Fig. 1. The UAV has only 280 mm wheelbase, and the LiDAR is directly installed on the airframe. The LiDAR-inertial odometry is sent to the flight controller tracking a circle trajectory (Fig. 3). The actual flight experiments show that FAST-LIO can achieve real-time and stable odometry output and mapping in a maximum of 50 Hz for the indoor environment. The flight trajectory and mapping results of 50 $H z$ frame rate indoor experiment are shown in Fig. 3. The average number of effective feature points and running time is 270 and 6.7 ms, respectively, the drift is smaller than 0.3% (0.08 m drift over 32 m trajectory). The flight video can be found at https://youtu.be/iYCY6T79oNU.

C. Indoor Experiments

Then we test FAST-LIO in a challenging indoor environment with large rotation speeds. In order to generate large rotation, the sensor suite is held on hands and undergoes quickly shaking. Fig. 4 shows the angular velocity and acceleration during the experiment. It is seen that the angular velocity often exceeds 100 deg/s. A state of the art implementation of LOAM on Livox LiDARs⁴ [11] and LOAM with IMU⁵ [9] are also tested as comparisons when the feature extraction are replaced with the one of FAST-LIO. The results show that FAST-LIO can output odometry faster and more stable than others, as shown in Fig. 5 and Table. III . It should be noted that the LOAM+IMU is a loosely-coupled method, hence results in inconsistent mapping. To further verify the mapping result, we perform a second experiment in the same environment but with a much slower motion. The map built by FAST-LIO is shown in the lower-right figure of Fig. 4. Since the two experiments have non-identical movements, it leads to slight visual differences at places occlusions occur. The rest mapping results are very close.

Fig. 4. The angular velocity and acceleration in the indoor experiments.

Fig. 5. The Mapping results of different LIO packages in an indoor environment with large rotation speed.

TABLE III
COMPARISON OF PROCESSING TIME FOR A LIDAR SCAN AT 10 HZ

Packages	Num. of effective features	Running time
LOAM	1107	59 ms
LOAM+IMU	1107	44 ms
FAST-LIO	1430	23 ms

Fig. 6. Mapping results of the Main Building, University of Hong Kong. The LiDAR platform, the same one in Fig. 1, is handheld randomly walking to scan the building. In order to show the drift, the experiment are started and ended at the same place.

Fig. 7. Comparison between LINS [22] and FAST-LIO. The data is from [22] and is collected by a Velodyne VLP-16 LiDAR and an Xsens MTiG-710 IMU, the green straight line in the center is the odometry output.

D. Outdoor Experiments

Here we show the performance of FAST-LIO in outdoor environments. Fig. 6 shows the mapping results (displaying all raw points) of the Main Building in the University of Hong Kong. The sensor suite is handheld during the data collection and returned to the starting position after traveling around 140m. The drift in this experiment is smaller than 0.05% (0.07 m drift over 140 m trajectory). The scan rate is set to 10 Hz in this experiment, and the average processing time of a scan is 25 ms with average 1497 effective feature points.

Further, we compare FAST-LIO with LINS⁶ [22]. To make a fair comparison, we use the dataset from LINS [22], which is a seaport area data collected by a Velodyne VLP-16 and an Xsens MTiG-710 IMU.⁶ The results show that the FAST-LIO can achieve better mapping accuracy (see Fig. 7) and only consumes 7.3 ms processing time in average while LINS takes 34.5 ms in average, both running at 10 Hz. It should be noted that since the EKF formula in LINS package has high computational complexity (see Section. III-C4)), it down samples the feature points to average 147 points in a scan (while 784 in a scan for FAST-LIO). This leads to degraded mapping accuracy for

LINS. The result in Fig. 7 shows all the feature points (before down-sample) of FAST-LIO and LINS. All the experiments are conducted on the DJI Manifold2 onboard computer.

V. CONCLUSION

This letter proposed FAST-LIO, a computationally efficient and robust LiDAR-inertial odometry framework by a tightly-coupled iterated Kalman filter. We used the forward and backward propagation to predict the states and compensate for the motion in a LiDAR scan. Besides, we proved and implemented an equivalent formula that can achieve much lower complexity for the Kalman gain computation. FAST-LIO was tested in the UAV flight experiment, challenging indoor environment with large rotation speed and outdoor environment. In all tests, our method produced precise, real-time, and reliable navigation results.

APPENDIX

A. Computation of $F_{\tilde{x}}$ and $F_{w}$

Recall $\begin{matrix} x_{i} = {\hat{x}}_{i} ⊞ {\tilde{x}}_{i} . \end{matrix}$ , denote $\mathbf { g } ( \widetilde { \mathbf { x } } _ { i } , \mathbf { w } _ { i } ) = \mathbf { f } ( \mathbf { x } _ { i } , \mathbf { u } _ { i } , \mathbf { w } _ { i } ) \Delta t = $ $ \mathbf { f } \left( \widehat { \mathbf { x } } _ { i } \right) $-$ \widetilde { \mathbf { x } } _ { i } , \mathbf { u } _ { i } , \mathbf { w } _ { i } ) \Delta t$ . Then the error state model (5) is rewriten as:

{\tilde{x}}_{i + 1} = \underset{G ({\tilde{x}}_{i}, H {\tilde{x}}_{i}) H g ({\tilde{x}}_{i}, w_{i})) \boxtimes ({\hat{x}}_{i} ⊞ g (0, 0))}{G ({\tilde{x}}_{i}, g ({\tilde{x}}_{i}, w_{i}))}

Following the chain rule of partial differention, the matrix $F_{\tilde{x}}$ and $F_{w}$ in (5) are computed as below.

\begin{matrix} F_{\tilde{x}} = {(\frac{\partial G ({\tilde{x}}_{i}, g (0, 0))}{\partial {\tilde{x}}_{i}} + \frac{\partial G (0, g ({\tilde{x}}_{i}, 0))}{\partial g ({\tilde{x}}_{i}, 0)} \frac{\partial g ({\tilde{x}}_{i}, 0)}{\partial {\tilde{x}}_{i}}) |}_{{\tilde{x}}_{i} = 0} \end{matrix}

\begin{matrix} F_{w} = {(\frac{\partial G (0, g (0, w_{i}))}{\partial g (0, w_{i})} \frac{\partial g (0, w_{i})}{\partial w_{i}}) |}_{w_{i} = 0} \end{matrix}

B. Equivalent Kalman Gain Formula

Based on the matrix inverse lemma [28], we can get:

{(P^{- 1} + H^{T} R^{- 1} H)}^{- 1} = P - P H^{T} {(H P H^{T} + R)}^{- 1} H P

Substituting above into (20), we can get:

\begin{matrix} K = {(H^{T} R^{- 1} H + P^{- 1})}^{- 1} H^{T} R^{- 1} \\ = P H^{T} R^{- 1} - P H^{T} {(H P H^{T} + R)}^{- 1} H P H^{T} R^{- 1} \end{matrix}

Now note that HPH $R^{- 1} = (H P H^{T} + R) R^{- 1} - I$ . Substituting it into above, we can get the standard Kalman gain formula in (18), as shown below.

\begin{matrix} K = P H^{T} R^{- 1} - P H^{T} R^{- 1} + P H^{T} {(H P H^{T} + R)}^{- 1} \\ = P H^{T} {(H P H^{T} + R)}^{- 1} . \end{matrix}

REFERENCES

[1] K. Sun et al., “Robust stereo visual inertial odometry for fast autonomous flight,” IEEE Robot. Automat. Lett., vol. 3, no. 2, pp. 965–972, Apr. 2018.

[2] A. Howard, “Real-time stereo visual odometry for autonomous ground vehicles,” in Proc. IEEE/RSJInt. Conf. Intell. Robots Syst., 2008, pp. 3946– 3952.

[3] T. Qin, P. Li, and S. Shen, “VINS-MONO: A robust and versatile monocular visual-inertial state estimator,” IEEE Trans. Robot., vol. 34, no. 4, pp. 1004–1020, Aug. 2018.

[4] C. Forster, M. Pizzoli, and D. Scaramuzza, “SVO: Fast semi-direct monocular visual odometry,” in Proc. IEEE Int. Conf. Robot. Automat. , 2014, pp. 15–22.

[5] D. Wang, C. Watkins, and H. Xie, “Mems mirrors for lidar: A review,” Micromachines, vol. 11, no. 5, pp. 456–479, 2020.

[6] Z. Liu, F. Zhang, and X. Hong, “Low-cost retina-like robotic lidars based on incommensurable scanning,” IEEE/ASME Trans. Mechatronics, pp. 1–1, 2021, early access, doi: 10.1109/TMECH.2021.3058173.

[7] P. J. Besl and N. D. McKay, “Method for registration of 3D shapes,” in Sensor fusion IV: control paradigms and data structures. Int. Soc. Opt. Photonics, vol. 1611, 1992, pp. 586–606.

[8] A. Segal, D. Haehnel, and S. Thrun, “Generalized-ICP,” in Proc. Robot.: Sci. Syst., Seattle, WA, vol. 2, no. 4, p. 435, 2009.

[9] J. Zhang and S. Singh, “LOAM: Lidar odometry and mapping in real-time,” in Robot.: Sci. Syst., vol. 2, no. 9, 2014.

[10] T. Shan and B. Englot, “LEGO-LOAM: Lightweight and groundoptimized lidar odometry and mapping on variable terrain,” in Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst., 2018, pp. 4758–4765.

[11] J. Lin and F. Zhang, “Loam livox: A fast, robust, high-precision lidar odometry and mapping package for LiDARs of small FoV,” in Proc. IEEE Int. Conf. Robot. Automat., 2020, pp. 3126–3131.

[12] W. Zhen, S. Zeng, and S. Soberer, “Robust localization and localizability estimation with a rotating laser scanner,” in Proc. IEEE Int. Conf. Robot. Automat., 2017, pp. 6240–6245.

[13] Y. Balazadegan Sarvrood, S. Hosseinyalamdary, and Y. Gao, “Visual-lidar odometry aided by reduced IMU,” ISPRS Int. J. Geo-inf., vol. 5, no. 1, p. 3, 2016.

[14] X. Zuo, P. Geneva, W. Lee, Y. Liu, and G. Huang, “Lic-fusion: Lidarinertial-camera odometry,” in Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst., 2019, pp. 5848–5854.

[15] P. Geneva, K. Eckenhoff, Y. Yang, and G. Huang, “Lips: Lidar-inertial 3 d plane slam,” in Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst., 2018, pp. 123–130.

[16] C. Forster, L. Carlone, F. Dellaert, and D. Scaramuzza, “On-manifold preintegration for real-time visual-inertial odometry,” IEEE Trans. Robot., vol. 33, no. 1, pp. 1–21, Feb. 2017.

[17] M. Hsiao, E. Westman, and M. Kaess, “Dense planar-inertial SLAM with structural constraints,” in Proc. IEEE Int. Conf. Robot. Automat., 2018, pp. 6521–6528.

[18] H. Ye, Y. Chen, and M. Liu, “Tightly coupled 3 d lidar inertial odometry and mapping,” in Proc. Int. Conf. Robot. Automat., 2019, pp. 3144–3150.

[19] A. Bry, A. Bachrach, and N. Roy, “State estimation for aggressive flight in gps-denied environments using onboard sensing,” in Proc. IEEE Int. Conf. Robot. Automat., 2012, pp. 1–8.

[20] J. A. Hesch, F. M. Mirzaei, G. L. Mariottini, and S. I. Roumeliotis, “A laser-aided inertial navigation system (l-ins) for human localization in unknown indoor environments,” in Proc. IEEE Int. Conf. Robot. Automat., 2010, pp. 5376–5382.

[21] Z. Cheng et al., “Practical phase unwrapping of interferometric fringes based on unscented Kalman filter technique,” Opt. Exp., vol. 23, no. 25, pp. 32337–32349, 2015.

[22] C. Qin, H. Ye, C. E. Pranata, J. Han, S. Zhang, and M. Liu, “Lins: A lidar-inertial state estimator for robust and efficient navigation,” in Proc. IEEE Int. Conf. Robot. Automat., 2020, pp. 8899–8906.

[23] M. Raitoharju and R. Piché, “On computational complexity reduction methods for kalman filter extensions,” IEEE Aerosp. Electron. Syst. Mag., vol. 34, no. 10, pp. 2–19, Oct. 2019.

[24] C. Hertzberg, R. Wagner, U. Frese, and L. Schröder, “Integrating generic sensor fusion algorithms with sound state representations through encapsulation of manifolds,” Inf. Fusion, vol. 14, no. 1, pp. 57–77, 2013.

[25] W. Xu, D. He, Y. Cai, and F. Zhang, “Robots state estimation and observability analysis based on statistical motion models,” 2020, arXiv:2010.05957.

[26] F. Bullo and R. M. Murray, “Proportional derivative (PD) control on the euclidean group,” in Proc. Eur. Control Conf., vol. 2, 1995, pp. 1091–1097.

[27] M. Bloesch, M. Burri, S. Omari, M. Hutter, and R. Siegwart, “Iterated extended kalman filter based visual-inertial odometry using direct photometric feedback,” Int. J. Robot. Res., vol. 36, no. 10, pp. 1053–1072, 2017.

[28] N. J. Higham, Accuracy and Stability ofNumerical Algorithms. Philadelphia, PA, USA: SIAM, 2002.

Authorized licensed use limited to: Jiangnan University. Downloaded on September 16,2026 at 03:27:52 UTC from IEEE Xplore. Restrictions apply.

本页由 md/2021_FAST-LIO__A_Fast__Robust_LiDAR-Inertial_Odometry_Package.md 预渲染生成。公式以 MathML 呈现，浏览器原生渲染，不依赖外部脚本或字体 CDN。