Accurate Monocular Visual-inertial SLAM using a Map-assisted EKF Approach

TL;DR

This paper introduces a real-time monocular visual-inertial SLAM system that combines an EKF with bundle adjustment and loop closure to achieve accurate, consistent, and long-term localization on a standard CPU.

Contribution

It proposes a novel tightly-coupled visual-inertial SLAM approach that integrates EKF, bundle adjustment, and loop closure for improved accuracy and consistency.

Findings

Outperforms existing methods on public datasets

Provides robust long-term localization in real-world scenarios

Achieves real-time performance on standard CPU hardware

Abstract

This paper presents a novel tightly-coupled monocular visual-inertial Simultaneous Localization and Mapping algorithm, which provides accurate and robust localization within the globally consistent map in real time on a standard CPU. This is achieved by firstly performing the visual-inertial extended kalman filter(EKF) to provide motion estimate at a high rate. However the filter becomes inconsistent due to the well known linearization issues. So we perform a keyframe-based visual-inertial bundle adjustment to improve the consistency and accuracy of the system. In addition, a loop closure detection and correction module is also added to eliminate the accumulated drift when revisiting an area. Finally, the optimized motion estimates and map are fed back to the EKF-based visual-inertial odometry module, thus the inconsistency and estimation error of the EKF estimator are reduced. In this way, the system can continuously provide reliable motion estimates for the long-term operation. The performance of the algorithm is validated on public datasets and real-world experiments, which proves the superiority of the proposed algorithm.

Tables

Table 1. TABLE I: Translation RMSE of the trajectories estimated from the proposed method and EKF VIO on the EuRoC MAV dataset

Sequence	Our Method With Loop	Our Method Without Loop	EKF
V1 01 easy	0.080	0.080	0.087
V1 02 medium	0.043	0.099	0.170
V1 03 difficult	0.124	0.245	0.301
V2 01 easy	0.052	0.052	0.082
V2 02 medium	0.042	0.042	0.191
V2 03 difficult	0.074	0.275	0.368
MH 01 easy	0.021	0.021	0.175
MH 02 easy	0.071	0.071	0.277
MH 03 medium	0.061	0.061	0.307
MH 04 difficult	0.064	0.064	0.309
MH 05 difficult	0.048	0.056	0.529

Table 2. TABLE II: Translation RMSE of the trajectories estimated from different approaches on the EuRoC MAV dataset. The best results are given in bold.

sequence	Our Method With Loop	Our Method Without Loop	VINS-MONO With Loop	VINS-MONO Without Loop	ORB-VIN	OKVIS	ROVIO
V1 01 easy	0.080	0.080	0.081	0.088	$0.027$	0.089	1.412
V1 02 medium	0.043	0.099	0.042	0.068	$0.028$	0.141	0.160
V1 03 difficult	$0.124$	0.245	0.156	0.160	X	0.262	0.170
V2 01 easy	0.052	0.052	0.063	0.068	$0.032$	0.135	0.236
V2 02 medium	0.042	0.042	0.066	0.084	$0.041$	0.155	0.408
V2 03 difficult	$0.074$	0.275	0.157	0.159	$0.074$	0.279	0.213
MH 01 easy	$0.021$	$0.021$	0.098377	0.301814	0.075	0.309	0.354
MH 02 easy	$0.071$	$0.071$	0.152	0.249	0.084	0.293	0.594
MH 03 medium	$0.061$	$0.061$	0.080	0.173	0.087	0.310	0.310
MH 04 difficult	$0.064$	$0.064$	0.129	0.323	0.217	0.360	1.058
MH 05 difficult	$0.048$	0.063	0.077	0.257	0.082	0.404	1.241

Equations

$$\mathbf{\xi }=\mathrm{Log}(\mathrm{R})={\log (\mathrm{R})}^{\vee }$$

$$\mathrm{R}(\mathbf{\xi })=\mathrm{Exp}(\mathbf{\xi })=\exp ({\mathbf{\xi }}^{\wedge })$$

$$\mathrm{R}=\hat{\mathrm{R}}\oplus \delta \mathbf{\xi }=\hat{\mathrm{R}}\mathrm{Exp}(\delta \mathbf{\xi })$$

$$\xi =\hat{\xi }\oplus \delta \mathbf{\xi }=\mathrm{Log}(\mathrm{Exp}(\hat{\xi })\mathrm{Exp}(\delta \mathbf{\xi }))$$

$$\begin{split}\widetilde{\bm{\mathrm{z}}}_{kl}&=\bm{\mathrm{z}}_{kl}+\bm{\sigma}_{kl}\\ &=\pi(\bm{\mathrm{f}}^{C_{k}}_{l})+\bm{\sigma}_{kl}\\ &=\left[\begin{array}[]{c}f_{x}\frac{x^{C_{k}}_{l}}{z^{C_{k}}_{l}}+c_{x}\\ f_{y}\frac{y^{C_{k}}_{l}}{z^{C_{k}}_{l}}+c_{y}\end{array}\right]+\bm{\sigma}_{kl}\\ \end{split}$$

$${\mathrm{X}}_k={[{\mathrm{X}}_{B_k}^{\mathrm{T}}{\mathrm{X}}_{L_k}^{\mathrm{T}}]}^{\mathrm{T}}$$

$${\mathrm{X}}_{B_k}={[{{\xi }_{B_k}^W}^{\mathrm{T}}{{\mathrm{p}}_{B_k}^W}^{\mathrm{T}}{{\mathrm{v}}_{B_k}^W}^{\mathrm{T}}{\mathrm{b}}_{a_k}^{\mathrm{T}}{\mathrm{b}}_{g_k}^{\mathrm{T}}]}^{\mathrm{T}}$$

$$\delta {\mathrm{X}}_{B_k}={[\delta {\xi }_{B_k}^{\mathrm{T}}\delta {{\mathrm{p}}_{B_k}^W}^{\mathrm{T}}\delta {{\mathrm{v}}_{B_k}^W}^{\mathrm{T}}\delta {\mathrm{b}}_{a_k}^{\mathrm{T}}\delta {\mathrm{b}}_{g_k}^{\mathrm{T}}]}^{\mathrm{T}}$$

$${\mathrm{X}}_{L_k}={[{{\mathrm{f}}_1^W}^{\mathrm{T}}\cdots {{\mathrm{f}}_m^W}^{\mathrm{T}}]}^{\mathrm{T}}$$

$${\mathrm{f}}_l^W={[x_l{y}_l{z}_l{\theta }_l{\varphi }_l{\rho }_l]}^{\mathrm{T}}$$

$$\delta {\mathrm{X}}_k={[\delta {\mathrm{X}}_{B_k}^{\mathrm{T}}\delta {{\mathrm{f}}_1^W}^{\mathrm{T}}\cdots \delta {{\mathrm{f}}_m^W}^{\mathrm{T}}]}^{\mathrm{T}}$$

$$\begin{array}{rl}& {\xi }_{B_{k|k-1}}^W=\mathrm{Log}\left(\mathrm{Exp}({\xi }_{B_{k-1}}^W)\mathrm{Exp}\left(({\tilde{\omega }}_{k-1}-{\mathrm{b}}_{g_{k-1}}-{\mathrm{n}}_{gd})\Delta t\right)\right)\\ & {\mathrm{p}}_{B_{k|k-1}}^W={\mathrm{p}}_{B_{k-1}}^W+{\mathrm{v}}_{B_{k-1}}^W\Delta t\\ & {\mathrm{v}}_{B_{k|k-1}}^W={\mathrm{v}}_{B_{k-1}}^W+({\mathrm{R}}_{B_{k-1}}^W({\tilde{\mathrm{a}}}_{k-1}-{\mathrm{b}}_{a_{k-1}}-{\mathrm{n}}_{ad})+{\mathrm{g}}^W)\Delta t\\ & {\mathrm{b}}_{a_{k|k-1}}={\mathrm{b}}_{a_{k-1}}\\ & {\mathrm{b}}_{g_{k|k-1}}={\mathrm{b}}_{g_{k-1}}\end{array}$$

$$\delta {\mathrm{X}}_{B_{k|k-1}}={\mathrm{\Phi }}_k\delta {\mathrm{X}}_{B_{k-1}}+{\mathrm{G}}_k{\mathrm{n}}_B$$

$$\begin{split}\bm{\mathrm{P}}_{k|k-1}&=\left[\begin{array}[]{cc}\bm{\mathrm{P}}_{B_{k|k-1}}&\bm{\mathrm{P}}_{BL_{k|k-1}}\\ \bm{\mathrm{P}}_{LB_{k|k-1}}&\bm{\mathrm{P}}_{L_{k|k-1}}\end{array}\right]\\ &=\left[\begin{array}[]{cc}\bm{\mathrm{\Phi}}_{k}\bm{\mathrm{P}}_{B_{k-1}}{\bm{\mathrm{\Phi}}_{k}}^{\mathrm{T}}+\bm{\mathrm{G}}_{k}\bm{\mathcal{Q}}{\bm{\mathrm{G}}_{k}}^{\mathrm{T}}&\bm{\mathrm{\Phi}}_{k}\bm{\mathrm{P}}_{BL_{{k-1}}}\\ \bm{\mathrm{P}}_{LB_{{k-1}}}{\bm{\mathrm{\Phi}}_{k}}^{T}&\bm{\mathrm{P}}_{L_{{k-1}}}\end{array}\right]\\ \end{split}$$

$$\bm{\mathrm{y}}_{l}^{W}=\left[\begin{array}[]{c}x_{l}\\ y_{l}\\ z_{l}\end{array}\right]+\frac{1}{\rho_{l}}\bm{\mathrm{m}}(\theta_{l},\phi_{l})$$

$$\bm{\mathrm{m}}(\theta_{l},\phi_{l})=\left[\begin{array}[]{c}\mathrm{cos}\phi_{l}\mathrm{sin}\theta_{l}\\ -\mathrm{sin}\phi_{l}\\ \mathrm{cos}\phi_{l}\mathrm{cos}\theta_{l}\end{array}\right]$$

$$\begin{array}{rl}& {\tilde{\mathrm{z}}}_{kl}={\mathrm{h}}_{kl}({\mathrm{X}}_{B_{k|k-1}},{\mathrm{f}}_l^W)+{\sigma }_{kl}\\ & =\pi \left({\mathrm{f}}_l^{C_k}\right)+{\sigma }_{kl}\end{array}$$

$$\begin{array}{rl}{\mathrm{r}}_{kl}& ={\tilde{\mathrm{z}}}_{kl}-{\mathrm{h}}_{kl}({\mathrm{X}}_{B_{k|k-1}},{({\mathrm{X}}_{L_k})}_l)\\ & \simeq {\mathrm{H}}_{B_{kl}}\delta {\mathrm{X}}_{B_{k|k-1}}+{\mathrm{H}}_{f_{kl}}\delta {\mathrm{f}}_l^W+{\sigma }_{kl}\end{array}$$

$${\mathrm{H}}_{kl}=[{\mathrm{H}}_{B_{kl}}0\cdots {\mathrm{H}}_{f_{kl}}0\cdots ]$$

$${\mathrm{K}}_k={\mathrm{P}}_{k|k-1}{\mathrm{H}}_k^{\mathrm{T}}{({\mathrm{H}}_k{\mathrm{P}}_{k|k-1}{\mathrm{H}}_k^{\mathrm{T}}+\mathrm{\Sigma })}^{-1}$$

$$\begin{array}{rl}{\mathrm{X}}_{k|k}& ={\mathrm{X}}_{k|k-1}\circ {\mathrm{K}}_k{\mathrm{r}}_k\\ {\mathrm{P}}_{k|k}& =(\mathrm{I}-{\mathrm{K}}_k{\mathrm{H}}_k){\mathrm{P}}_{k|k-1}\end{array}$$

$${[\breve{x_l}\breve{y_l}\breve{z_l}]}^{\mathrm{T}}={\mathrm{R}}_{B_{k|k}}^W{\mathrm{p}}_C^B+{\mathrm{p}}_{B_{k|k}}^W$$

$$\left[\begin{array}[]{c}\breve{\theta_{l}}\\ \breve{\phi_{l}}\end{array}\right]=\left[\begin{array}[]{c}\mathrm{arctan}(x_{kl}^{W},z_{kl}^{W})\\ \mathrm{arctan}(-y_{kl}^{W},\sqrt{{x_{kl}^{W}}^{2}+{z_{kl}^{W}}^{2}})\end{array}\right]$$

$$\bm{\mathrm{\tau}}_{{kl}}^{W}=\left[\begin{array}[]{c}x_{kl}^{W}\\ y_{kl}^{W}\\ z_{kl}^{W}\end{array}\right]=\bm{\mathrm{R}}^{W}_{B_{k|k}}\bm{\mathrm{R}}^{B}_{C}\left[\begin{array}[]{c}\frac{u_{kl}-c_{x}}{f_{x}}\\ \frac{v_{kl}-c_{y}}{f_{y}}\\ 1\end{array}\right]$$

$$\bm{\mathrm{P}}_{k|k}\leftarrow\bm{\mathrm{J}}\left[\begin{array}[]{cc}\bm{\mathrm{P}}_{k|k}&\bm{\mathrm{0}}_{(15+6m)\times 6}\\ \bm{\mathrm{0}}_{6\times(15+6m)}&\bm{\mathrm{\Sigma}}_{h\rho}\end{array}\right]{\bm{\mathrm{J}}}^{\mathrm{T}}$$

$$\mathcal{X}=\{{\mathrm{X}}_{B_i},\cdots ,{\mathrm{X}}_{B_j},{\mathrm{L}}_1^W,\cdots ,{\mathrm{L}}_n^W\}$$

$$\begin{array}{rl}f(\mathcal{X})=& \Vert {\mathrm{r}}_p{\Vert }_{{\mathbf{\Sigma }}_p}^2+\sum _{k\in \mathcal{K}}\rho \left(\Vert {\mathrm{r}}_{{\mathcal{I}}_{k-1k}}{\Vert }_{{\mathbf{\Sigma }}_{{\mathcal{I}}_{k-1k}}}^2\right)\\ & +\sum _{k\in \mathcal{K},l\in \mathcal{C}}\rho \left(\Vert {\mathrm{r}}_{{\mathcal{C}}_{kl}}{\Vert }_{{\mathbf{\Sigma }}_{{\mathcal{C}}_{kl}}}^2\right)\end{array}$$

$${\mathcal{X}}^{\ast }=\underset{\mathcal{X}\in \mathcal{M}}{\mathrm{argmin}}f(\mathcal{X})$$

$$\begin{array}{rl}& \Delta {\tilde{\mathrm{R}}}_{ij}=\prod _{k=i}^{j-1}\mathrm{Exp}\left(({\tilde{\omega }}_k-{\mathrm{b}}_{g_i})\Delta t\right)\\ & \Delta {\tilde{\mathrm{p}}}_{ij}=\sum _{k=i}^{j-1}\left(\Delta {\tilde{\mathrm{v}}}_{ik}\Delta t+\frac{1}{2}\Delta {\tilde{\mathrm{R}}}_{ik}({\tilde{\mathrm{a}}}_k-{\mathrm{b}}_{a_i})\Delta t^2\right)\\ & \Delta {\tilde{\mathrm{v}}}_{ij}=\sum _{k=i}^{j-1}\Delta {\tilde{\mathrm{R}}}_{ik}({\tilde{\mathrm{a}}}_k-{\mathrm{b}}_{a_i})\Delta t\end{array}$$

$$\begin{array}{rl}& \Delta {\tilde{\mathrm{R}}}_{ij}({\mathrm{b}}_{g_i})=\Delta {\tilde{\mathrm{R}}}_{ij}({\bar{\mathrm{b}}}_{g_i})\mathrm{Exp}\left(\frac{\partial \Delta {\bar{\mathrm{R}}}_{ij}}{\partial {\mathrm{b}}_g}\delta {\mathrm{b}}_g\right)\\ & \Delta {\tilde{\mathrm{p}}}_{ij}({\mathrm{b}}_{g_i},{\mathrm{b}}_{a_i})=\Delta {\tilde{\mathrm{p}}}_{ij}({\bar{\mathrm{b}}}_{g_i},{\bar{\mathrm{b}}}_{a_i})+\frac{\partial \Delta {\bar{\mathrm{p}}}_{ij}}{\partial {\mathrm{b}}_g}\delta {\mathrm{b}}_g+\frac{\partial \Delta {\bar{\mathrm{p}}}_{ij}}{\partial {\mathrm{b}}_a}\delta {\mathrm{b}}_a\\ & \Delta {\tilde{\mathrm{v}}}_{ij}({\mathrm{b}}_{g_i},{\mathrm{b}}_{a_i})=\Delta {\tilde{\mathrm{v}}}_{ij}({\bar{\mathrm{b}}}_{g_i},{\bar{\mathrm{b}}}_{a_i})+\frac{\partial \Delta {\bar{\mathrm{v}}}_{ij}}{\partial {\mathrm{b}}_g}\delta {\mathrm{b}}_g+\frac{\partial \Delta {\bar{\mathrm{v}}}_{ij}}{\partial {\mathrm{b}}_a}\delta {\mathrm{b}}_a\end{array}$$

Full text

I Introduction

Concurrent motion estimation and map reconstruction by combining visual and inertial measurements has received significant interest in the field of Robotics and Computer Vision communities. This visual-inertial sensor suite can serve as an ideal alternative to GPS in environments where GPS is denied, since both sensors are small, lightweight, cheap enough and complementary. On the one hand, Visual SLAM can provide good tracking and rich map information in visually distinguishable environments. However due to the sensor limitation, visual simultaneous localization and mapping(SLAM) is sensitive to motion blur, occlusions and illumination changes. On the other hand, inertial sensors are able to provide self-motion information at high frequency, so inertial navigation is robust to aggressive motion and can provide absolute scale for the motion. Whereas the result of inertial navigation is noisy and diverges even in a few seconds. Therefore, fusing measurements from the inertial sensor to the visual SLAM in a tightly-coupled way, both the robustness and the accuracy of motion tracking can be dramatically improved. The advantage of visual-inertial sensor fusion is the most obvious in a monocular visual-inertial setup, because the scale of the motion estimation and map structure computed from monocular SLAM is ambiguous and liable to drift over time.

In this paper, we aim to build a system which enables the accurate and robust motion tracking within the accurately and consistently reconstructed map. Firstly, we employ the EKF-based visual-inertial odometry(VIO) to track the 3D motion of the IMU body frame in real time. Whereas due to the linearization errors, the estimator tends to be inconsistent, which results in large estimation errors and divergence. So we intend to perform bundle adjustment(BA) to reduce the inconsistency of the estimator. It is demonstrated in [1] that the BA techniques can achieve better accuracy than filtering techniques, because optimization methods can iteratively relinearize measurement equations to better deal with their nonlinearity. However, the optimization method leads to a high computational cost. In addition, increasing the number of feature correspondences and keyframes in the window of local BA will lead to significant increase in accuracy, whereas feature extraction and matching, optimization for the local BA are also time consuming. So if we increase the number of features, and perform the local BA after the EKF for each frame to increase the accuracy of the system, the system will incapable of real time operation. Thus, in order to ensure both real time operation and accurate global map reconstruction, we extract a vast number of features only for those selected keyframes, we use these features to construct global map and perform the local BA in a parallel thread. In this way, we can properly solve the loss of accuracy due to the linearization in real-time. Besides, if the system is unable to close loops, the drift of the estimated trajectory will accumulate without bound, even if the sensor is continuously revisiting the same place. Therefore, we add a loop closure module in a new parallel thread for reducing the accumulated drift when returning to an already mapped area.

In summary, we propose a tightly-coupled monocular visual-inertial SLAM(VISLAM) system, which is able to perform real-time, accurate, robust and long term localization and map reconstruction. Our approach operates in three parallel threads, one thread is used to perform the EKF VIO, and the other two threads, one for BA and the other one for loop closing, are used to compensate for the error growth of the estimated trajectory and construct an accurate and consistent global map. In VIO thread, since the computational cost of the EKF is quadratic in the number of features, for the real-time state estimation we extract the appropriate number of features for each frame, and extract a lot of features only in selected keyframes for performing local BA and loop closure. Finally, we design a feedback mechanism to increase both the consistency and accuracy of the EKF estimator, which is achieved by feeding back the optimized state and globally consistent map to update the state vector of the EKF VIO module.

In experiments, the results demonstrate the benefits of our system towards the EKF based VIO. We also compare to the state-of-the-art VIO and VISLAM approaches, and demonstrate the superior performance of our method.

Our monocular visual-inertial SLAM algorithm is shown in Fig.1. The system is complete and drift-free in large scale environments. The remainder of the paper is organized as follows. In Section II, we describe the relevant literature. Notations are given in Section III. In Section IV, EKF based VIO is presented and the Jacobian matrices of EKF are given in appendix. A tightly-coupled joint visual-inertial nonlinear optimization is introduced In Section V. Section VI introduce our tightly-coupled VISLAM approach. Experiments results are shown in Section VII. Finally, the paper is concluded in Section VIII.

II RELATED WORK

There are a vast amount of research towards visual SLAM problem, we refer to the review paper [2] for the progress made in the past few decades. In this section, we will discuss the most relevant works on monocular VIO and VISLAM system.

The fusion for visual and inertial measurements is usually divided into two classes. Loosely-coupled approaches, e.g. [3][4], process the visual and inertial measurements separately. Therefore the accumulated drift in vision module cannot be eliminated from the usability of inertial measurements, which leads the resulted estimate to be sub-optimal. Tightly-coupled ones interested in this work, e.g. [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19], perform VIO or VISLAM system by considering the tight interaction between visual and inertial measurements to optimally exploit the both sensing cues, thus achieve higher precision at the expensive of additional computational complexity. Besides, for tightly-coupled VIO/VISLAM solutions, two methodologies have been prevalent: filtering-based methods [5, 6, 7, 8, 9, 10, 11, 12, 13] and BA-based methods [14, 15, 16, 17, 18, 19, 20].

Historically, the monocular SLAM problem has been addressed with filtering method, which operates on the mean and covariance of the probabilistic distribution in a kalman filtering framework. Filtering based approaches require fewer computational resources due to the continuous marginalization of past state, however the system get slightly lower accuracy due to the linearization error. According to the way processing the measurement information, the recursive filtering approaches can be classified into two main categories: extended kalman filter (EKF) based methods [5, 6, 7, 8, 9] and sliding window filtering approaches [10, 11, 12, 13]. The state vector of EKF-based SLAM algorithms include both the pose of the platform and a set of feature positions, so as long as these features are continuously observed and contained in the state vector, the estimated pose relative to these features will not drift. However it have high computational complexity (quadratic in the number of features in the state vector), therefore only currently observable landmarks are tracked to ensure real-time operation. In contrast, sliding window filtering approaches maintain a sliding window of past camera poses in the state vector, and use the feature measurements to impose probabilistic constraints on these poses, therefore keep computational complexity only linear in the number of features by excluding point features from filter state vector. Generally, the VIO problem has four unobservable directions, but since the linearization errors, the system only have three unobservable directions, which renders the filter inconsistent. So papers [8, 9, 11, 12, 13] were proposed a series of methods, e.g. first-estimates Jacobian and constraint of system observability, to improve the consistency of the system. If the measurement models were linear, both methods yield the same result equal to the MAP estimate.

In [1], it was shown that nonlinear optimization-based approaches provide better accuracy than filtering approaches by its capability to relinearize the state at each iteration, therefore avoiding integrated error from linearization, however it leads to higher computational demands. In following, we introduce several classic tightly-coupled BA-based VIO and VISLAM system. OKVIS [14] presented an approach to tightly integrate inertial measurements into keyframe-based visual SLAM, which makes the nonlinear cost function comprised of IMU error term with the landmark reprojection error term to be jointly optimized. Additionally, marginalization of old state is used to maintain a bounded-sized optimization window, therefore OKVIS has achieved increased accuracy and robustness in real-time operation. However, the system needs to repeatedly compute the IMU integration when the linearization point changes. To eliminate this repeated computation, Foster et al. provided a preintegration theory for inertial measurements in [15] that properly address the manifold structure of the rotation group based on the work of [21]. Then, the preintegrated IMU model and structureless visual model are seamlessly integrated in a fully probabilistic manner to build a much more computationally efficient optimization method for state estimation. Therefore, the system achieves better accuracy than Project Tango [22] by using SVO as front-end and the visual-inertial joint optimization in back-end. Tightly-coupled visual-inertial odometry methods mentioned above all lack the capability to close loops and reuse an already reconstructed map, due to the marginalization of past states to maintain a constant computational cost or the use of full smoothing. Thus, ORB-VISLAM [17] presented a real-time tightly-coupled monocular visual-inertial SLAM system, which enables the loop closure and the reuse of previously estimated 3D map. The system achieved higher accuracy than the fully direct, stereo visual-inertial odometry [18], because there is no drift accumulation for localization in already mapped areas. Recently, a novel real-time, tightly-coupled, sliding window optimization based versatile monocular visual-inertial odometry [19][20] was proposed, in which, the state of the system and a representation of the environment are estimated by local BA in one thread, and loops are closed in lightweight manner in parallel thread.

Both filtering-based methods and BA-based methods have their merits, so in this work, we tightly fuse both methods to achieve the best accuracy, robustness and efficiency.

III NOTATIONS

Throughout the paper, we denote the world reference frame as $(\cdot)^{W}$, and denote the IMU body frame and camera frame for the $k^{th}$ image as $(\cdot)^{B_{k}}$ and $(\cdot)^{C_{k}}$ respectively. In addition, we employ $\bm{\mathrm{R}}^{\mathcal{F}_{1}}_{\mathcal{F}_{2}}\in\bm{\mathrm{SO}}(3)$ to represent rotation from frame $\{\mathcal{F}_{2}\}$ to $\{\mathcal{F}_{1}\}$, ${\bm{\mathrm{p}}}^{\mathcal{F}_{1}}_{\mathcal{F}_{2}}\in\mathbb{R}^{3}$ and ${\bm{\mathrm{v}}}^{\mathcal{F}_{1}}_{\mathcal{F}_{2}}\in\mathbb{R}^{3}$ to describe the 3D position and velocity of frame $\{\mathcal{F}_{2}\}$ with respect to the frame $\{\mathcal{F}_{1}\}$. Besides, the rotation and translation between the rigidly mounted camera-IMU sensor are denoted as $\bm{\mathrm{R}}^{B}_{C}$ and $\bm{\mathrm{p}}^{B}_{C}$, which was computed from the calibration.

For the over-parameterized rotation matrix, a vector $\bm{\mathrm{\xi}}\in\mathbb{R}^{3}$ can be computed from the tangent space $\mathfrak{so}(3)$ of manifold $\bm{\mathrm{SO}}(3)$ to provide a minimal representation. The Lie algebra $\mathfrak{so}(3)$ and Lie group $\bm{\mathrm{SO}}(3)$ are related through the logarithm map and exponential map:

[TABLE]

[TABLE]

where ${(\cdot)}^{\wedge}$ operator maps a vector in $\mathbb{R}^{3}$ to a $3\times 3$ skew symmetric matrix, and ${(\cdot)}^{\vee}$ is the inverse operator. The formula for $\mathrm{log}(\cdot)$ and $\mathrm{exp}(\cdot)$ can be found in [23].

Furthermore, the uncertainty of rotation $\bm{\mathrm{R}}\in\bm{\mathrm{SO}}(3)$ and $\bm{\mathrm{\xi}}\in\mathbb{R}^{3}$ are described as:

[TABLE]

[TABLE]

where $\hat{\bm{\mathrm{R}}}$ and $\hat{\bm{\mathrm{\xi}}}$ are the mean estimate of $\bm{\mathrm{R}}$ and $\bm{\mathrm{\xi}}$ respectively, and $\delta\bm{\xi}\sim\mathcal{N}(\bm{0},\bm{\mathrm{\Sigma}})$ is a normally distributed perturbation.

For visual measurements, we consider a projection function $\pi:\mathbb{R}^{3}\to\Omega\subset\mathbb{R}^{2}$, which projects the $l^{th}$ map point expressed in the current camera frame $\bm{\mathrm{f}}^{C_{k}}_{l}={[x^{C_{k}}_{l}\ y^{C_{k}}_{l}\ z^{C_{k}}_{l}]}^{\mathrm{T}}\in\mathbb{R}^{3}$ onto 2D points on the image plane $\bm{\mathrm{z}}_{kl}={[u_{kl}\ v_{kl}]}^{\mathrm{T}}\in\Omega$:

[TABLE]

where $\widetilde{\bm{\mathrm{z}}}_{kl}$ is the corresponding feature measurement, and $\bm{\sigma}_{kl}$ is the $2\times 1$ measurement noise with covariance $\bm{\Sigma}_{\sigma_{kl}}$ associated to the feature scale. In addition, $f_{x}$, $f_{y}$ are focal length and $c_{x}$, $c_{y}$ are principle point, which are known from calibration.

IV VIO Description

In this section, we describe the EKF VIO system, which is based on the work of [24]. An overview of the algorithm is given in VIO part of Fig. 1. Inertial measurements are used to predict the motion movement in the prediction stage, then the state is updated using the matched visual features. In this way, we can estimate the state of the body frame efficiently.

IV-A Full State Vector

The state vector to be estimated comprises the IMU state $\bm{\mathrm{X}}_{B_{k}}$ and a set of landmark parameters ${\bm{\mathrm{X}}_{L_{k}}}$:

[TABLE]

The IMU state is formulated by the vector:

[TABLE]

where ${\bm{\mathrm{\xi}}^{W}_{B_{k}}}\in\mathbb{R}^{3}$ is the Lie algebra of orientation ${\bm{\mathrm{R}}^{W}_{B_{k}}}\in\bm{\mathrm{SO}}(3)$ from frame $\{B_{k}\}$ to $\{W\}$, ${\bm{\mathrm{p}}^{W}_{B_{k}}}\in\mathbb{R}^{3}$ and ${\bm{\mathrm{v}}^{W}_{B_{k}}}\in\mathbb{R}^{3}$ are the 3D position and velocity of frame $\{B_{k}\}$ with respect to $\{W\}$, as well as $\bm{\mathrm{b}}_{a}\in\mathbb{R}^{3}$ and $\bm{\mathrm{b}}_{g}\in\mathbb{R}^{3}$ are additive accelerometer and gyroscope biases respectively. Following (7), the IMU error state vector is defined as:

[TABLE]

where we use the standard additive error for the 3D position, velocity and biases, while for rotation, the error is defined as (4).

Assuming that $m$ features are included in the map at time-step $k$, then the coordinates of features are:

[TABLE]

The position of the $l^{th}$ landmark ${\bm{\mathrm{f}}_{l}^{W}}$ is paramterized in inverse depth coordinates as:

[TABLE]

where ${(x_{l},y_{l},z_{l})}^{\mathrm{T}}$ is the camera position, in which the $l^{th}$ landmark was firstly observed. $\theta_{l}$ and $\phi_{l}$ are the azimuth and elevation angle defining unit ray(expressed in the world frame) that goes from the camera center ${(x_{l},y_{l},z_{l})}^{\mathrm{T}}$ to the $l^{th}$ landmark, and ${\rho_{l}}$ is its inverse depth along the unit ray.

Therefore the EKF error state vector is expressed as:

[TABLE]

where we use the standard additive error for the landmark.

IV-B IMU Propagation Model

Different to many other visual-inertial methods, we define the state propagation model directly in discrete-time to make the derivatives required for the EKF prediction are available in close-form. Using the measured acceleration $\widetilde{\bm{\mathrm{a}}}_{k-1}$ and angular velocity $\widetilde{\bm{\mathrm{\omega}}}_{k-1}$ obtained from IMU, the discrete-time propagation model $\bm{\mathrm{X}}_{B_{k|k-1}}=\bm{\mathrm{f}}_{k}(\bm{\mathrm{X}}_{B_{k-1}})$ is:

[TABLE]

where $\bm{\mathrm{R}}^{W}_{B_{k-1}}=\mathrm{Exp}({\bm{\mathrm{\xi}}^{W}_{B_{k-1}}})$, $\bm{\mathrm{n}}_{gd}\sim\mathcal{N}(\bm{\mathrm{0}},{\bm{\Sigma}_{g}}/{\Delta t})$ and $\bm{\mathrm{n}}_{ad}\sim\mathcal{N}(\bm{\mathrm{0}},\bm{\Sigma}_{a}/\Delta t)$ are discrete-time white Gaussian noise for inertial measurements, and $\bm{\mathrm{g}}^{W}$ is the constant gravity. In this work, we ignore the slow random walk of the inertial biases, so the biases are considered fixed and estimated as part of the state. The linearized discrete-time IMU error state propagation model is represented as:

[TABLE]

where $\bm{\mathrm{n}}_{B}={[{\bm{\mathrm{n}}_{ad}}^{\mathrm{T}}\ \ {\bm{\mathrm{n}}_{gd}}^{\mathrm{T}}]}^{\mathrm{T}}$ is the system noise with covariance \bm{\mathcal{Q}}=\left[\begin{array}[]{cc}\bm{\Sigma}_{a}/\Delta t&\bm{\mathrm{0}}_{3\times 3}\\ \bm{\mathrm{0}}_{3\times 3}&\bm{\Sigma}_{g}/\Delta t\end{array}\right]. The matrices $\bm{\mathrm{\Phi}}_{k}$ and $\bm{\mathrm{G}}_{k}$ in (13) are Jacobians of $\bm{\mathrm{f}}_{k}(\cdot)$ with respect to IMU state and the system noise, these representation can be found in appendix A.

Therefore, the covariance matrix is propagated as follows:

[TABLE]

IV-C Measurement Model

The inverse depth parameterization is used for features to (1) enhance the degree of linearity for measurement equations, and (2) better deal with the low parallax features. The Euclidean XYZ coordinates of the $l^{th}$ feature in world frame $\bm{\mathrm{y}}_{l}^{W}$ can be transformed from its inverse depth representation $\bm{\mathrm{f}}_{l}^{W}$ as:

[TABLE]

[TABLE]

Thus the measurement model representing the projection of the $l^{th}$ landmark to the ${k}^{th}$ image is:

[TABLE]

where \bm{\mathrm{f}}_{l}^{C_{k}}=\bm{\mathrm{R}}_{W}^{C_{k|k-1}}\left(\rho_{l}\left(\left[\begin{array}[]{c}x_{l}\\ y_{l}\\ z_{l}\end{array}\right]-\bm{\mathrm{p}}^{W}_{C_{k|k-1}}\right)+\bm{\mathrm{m}}(\theta_{l},\phi_{l})\right), $\bm{\mathrm{R}}_{W}^{C_{k|k-1}}={(\bm{\mathrm{R}}^{W}_{B_{k|k-1}}\bm{\mathrm{R}}^{B}_{C})}^{\mathrm{T}}$ and $\bm{\mathrm{p}}^{W}_{C_{k|k-1}}={\bm{\mathrm{p}}^{W}_{B_{k|k-1}}}+\bm{\mathrm{R}}^{W}_{B_{k|k-1}}{\bm{\mathrm{p}}^{B}_{C}}$. From the measurement model, we compute the reprojection error and its linearized approximation as:

[TABLE]

where the matrices $\bm{\mathrm{H}}_{B_{kl}}$ and $\bm{\mathrm{H}}_{f_{kl}}$ are derivatives of the measurement model with respect to the IMU state estimate and the $l^{th}$ feature position respectively, their expressions are given in appendix B.

Therefore, we can obtain the measurement Jacobian matrix as:

[TABLE]

IV-D Filter Update

To perform an update of the estimated state, we stack the $m$ individual measurement residual $\bm{\mathrm{r}}_{kl}$ at time-step $k$ together to form a single $2m\times 1$ residual vector $\bm{\mathrm{r}}_{k}={\left[\bm{\mathrm{r}}_{k1}^{\mathrm{T}}\cdots\bm{\mathrm{r}}_{kl}^{\mathrm{T}}\cdots\bm{\mathrm{r}}_{km}^{\mathrm{T}}\right]}^{\mathrm{T}}$. In the same way, the measurement Jacobians are also combined to a single $2m\times n$ measurement matrix $\bm{\mathrm{H}}_{k}={\left[\bm{\mathrm{H}}_{k1}^{T}\cdots\bm{\mathrm{H}}_{kl}^{\mathrm{T}}\cdots\bm{\mathrm{H}}_{km}^{\mathrm{T}}\right]}^{\mathrm{T}}$. Then the kalman gain is computed as:

[TABLE]

where $\bm{\mathrm{\Sigma}}$ is the stacked $2m\times 2m$ measurement uncertainty. Finally, the full state and covariance are updated by:

[TABLE]

where $\circ$ operator is equal to the $\oplus$ operator in (4) for the orientation and the addition of vector for other state.

Based on the predicted IMU state, features used to update the state are searched with the optical flow method. Prior to using each feature’s measurement to update, for each matched feature, the Mahalanobis distance $d=\bm{\mathrm{r}}_{kl}^{\mathrm{T}}{(\bm{\mathrm{H}}_{kl}\bm{\mathrm{P}}_{k|k-1}\bm{\mathrm{H}}_{kl}^{\mathrm{T}}+\bm{\Sigma}_{\sigma_{kl}})}^{-1}\bm{\mathrm{r}}_{kl}$ is firstly computed to reject outliers. If $d$ is smaller than a threshold given by the 95-th percentile of the $\chi^{2}$ distribution, the feature is accepted as an inlier, and used for the filter update. We perform the filter update by combining 1-point RANSAC method as in [25] to find reliable inliers.

IV-E State Augmentation

Once a new feature $l$ is needed at time-step $k$, the initial position $\bm{\mathrm{f}}_{l}^{W}={[\breve{x_{l}}\ \breve{y_{l}}\ \breve{z_{l}}\ \breve{\theta_{l}}\ \breve{\phi_{l}}\ \breve{\rho_{l}}]}^{\mathrm{T}}$ of the new feature is computed as follows. The camera position is computed by:

[TABLE]

Besides, from the observation ${[u_{kl}\ v_{kl}]}^{\mathrm{T}}$ of the new feature in the image, the angles $\breve{\theta_{l}}$ and $\breve{\phi_{l}}$ defining its direction are calculated as:

[TABLE]

[TABLE]

The initial value for ${\breve{\rho}_{i}}$ and its standard deviation $\sigma_{\rho}$ are set as in [26]. Then the initial position for new feature is appended to the state vector, and the state covariance matrix is also augmented accordingly:

[TABLE]

where \bm{\mathrm{\Sigma}}_{h\rho}=\left[\begin{array}[]{cc}\bm{\Sigma}_{\sigma_{kl}}&\bm{\mathrm{0}}_{2\times 1}\\ \bm{\mathrm{0}}_{1\times 2}&{\sigma_{\rho}}^{2}\end{array}\right] denotes the uncertainty of the visual measurements and the initial inverse depth, for the Jacobian $\bm{\mathrm{J}}$ we refer the reader to appendix C.

V VISUAL-INERTIAL BUNDLE ADJUSTMENT

Once a frame processed by EKF based VIO is selected as a keyframe, we apply a nonlinear optimization in a sliding window to improve the accuracy of the estimated state. In this section, we combine the visual and inertial measurements in an unified formulation.

V-A Bundle Adjustment Representation

We formulate a joint optimization problem to optimally estimate the full state in a sliding window using all the available inertial and visual measurements. Full state in sliding window contain a set of successive keyframes from $i$ to $j$ and $n$ landmarks visible by the keyframes in sliding window, which is denoted as:

[TABLE]

where $\bm{\mathrm{L}}_{l}^{W}$ is the position of the 3D map point expressed in Euclidean XYZ coordinates. We denote the keyframes and visual measurements in sliding window as $\mathcal{K}$ and $\mathcal{C}$ respectively. Then the energy function that we want to minimize is given by:

[TABLE]

where $\bm{\mathrm{r}}_{p}$, $\bm{\mathrm{r}}_{\mathcal{I}_{k-1k}}$, $\bm{\mathrm{r}}_{\mathcal{C}_{kl}}$ are prior error, temporal IMU error and reprojection error respectively, as well as $\bm{\Sigma}_{p}$, ${\bm{\Sigma}}_{\mathcal{I}_{k-1k}}$, ${\bm{\Sigma}}_{\mathcal{C}_{kl}}$ are the corresponding covariance matrices, and $\rho$ is the Huber robust cost function. The optimization problem can be interpreted as a factor graph shown in Fig. 2.

Therefore the best estimate for variable $\bm{\mathcal{X}}$ can be obtained by minimizing the objective function on manifold:

[TABLE]

Detailed IMU and visual residuals are provided in the following subsections. The least squares problem are solved by Gauss-Newton method implemented in g2o[27] or ceres solver.

V-B Inertial Measurement Model

IMU measurements arrive at a much higher frequency than the visual measurements. So in order to avoid the frequent integration whenever the linearization point changes, we adopt the IMU preintegration approach proposed in [15]. The IMU preintegraton is independent of the initial conditions, and can incorporate the change of IMU biases. The concept was firstly proposed in [21].

We integrate all the IMU measurements $\{\widetilde{\bm{\mathrm{a}}}_{k},\ \widetilde{\bm{\mathrm{\omega}}}_{k}\}$ between keyframes i and j to compute the IMU preintegration $\Delta{\widetilde{\bm{\mathrm{I}}}_{ij}}=[\Delta\widetilde{\bm{\mathrm{R}}}_{ij},\Delta\widetilde{\bm{\mathrm{p}}}_{ij},\Delta\widetilde{\bm{\mathrm{v}}}_{ij}]$ on manifold as:

[TABLE]

Furthermore, given a bias update $\delta\bm{\mathrm{b}}$ and using the first-order expansion, the preintegrated IMU measurement can be updated as:

[TABLE]

where Jacobians $\frac{\partial\Delta(\cdot)}{\partial\bm{\mathrm{b}}_{\cdot}}$ describe how a change in the bias estimate effects the preintegrated IMU measurements. The derivation of the Jacobians can be found in [15]. From geometric constraints, we get the IMU measurement as:

[TABLE]

where ${[\delta\bm{\mathrm{\xi}}_{ij}^{\mathrm{T}}\ \delta\bm{\mathrm{t}}_{ij}^{\mathrm{T}}\ \delta\bm{\mathrm{v}}_{ij}^{\mathrm{T}}]}^{\mathrm{T}}\in\mathcal{N}(\bm{\mathrm{0}},\bm{\Sigma}_{\mathcal{I}_{ij}})$ is zero-mean white Gaussian noise. Given the above measurement model, the IMU preintegration residual ${\bm{\mathrm{r}}_{\Delta_{ij}}}={[\bm{\mathrm{r}}_{\Delta\bm{\mathrm{R}}_{ij}}^{\mathrm{T}}\ \bm{\mathrm{r}}_{\Delta\bm{\mathrm{p}}_{ij}}^{\mathrm{T}}\ \bm{\mathrm{r}}_{\Delta\bm{\mathrm{v}}_{ij}}^{\mathrm{T}}]}^{\mathrm{T}}\in\mathbb{R}^{9}$ is:

[TABLE]

The corresponding covariance matrix ${\bm{\Sigma}}_{\Delta_{ij}}$ can be calculated by incrementally propagating the preintegration noise from keyframe $i$ to $j$. The detailed derivatives about Jaocbians and the uncertainty propagation on manifold can refer to paper [15].

V-C Bias Model

Biases are slowly time-varying, so for the biases between consecutive keyframes i and j, we have:

[TABLE]

where $\bm{\mathrm{\eta}}_{b_{gd}}$ and $\bm{\mathrm{\eta}}_{b_{ad}}$ are the discretized bias random walk with covariance $\bm{\Sigma}_{b_{gd}}$ and $\bm{\Sigma}_{b_{ad}}$. Therefore, we express the bias error ${\bm{\mathrm{r}}_{b}=[\bm{\mathrm{r}}_{g}^{\mathrm{T}}\ \bm{\mathrm{r}}_{a}^{\mathrm{T}}}]\in\mathbb{R}^{6}$ as:

[TABLE]

V-D Visual Measurement Model

Through the measurement model in (LABEL:visualmm), the reprojection residual $\bm{\mathrm{r}}_{\mathcal{C}_{kl}}\in\mathbb{R}^{2}$ for the $l^{th}$ map point seen by the $k^{th}$ keyframe is:

[TABLE]

The corresponding covariance matrix ${{\bm{\Sigma}}_{\mathcal{C}_{kl}}}$ is equal to $\bm{\Sigma}_{\sigma_{kl}}$.

V-E Error Term Representation

In this section, we give detailed representation for the IMU error term $\|{\bm{\mathrm{r}}_{\mathcal{I}_{ij}}}\|^{2}_{{\bm{\Sigma}}_{\mathcal{I}_{ij}}}$ and the reprojection error term $\|{\bm{\mathrm{r}}_{\mathcal{C}_{kl}}}\|^{2}_{{\bm{\Sigma}}_{\mathcal{C}_{kl}}}$ in (27). Given the inertial measurement model in Section V-B and bias model in Section V-C, the IMU error term is:

[TABLE]

where {\bm{\Sigma}}_{b_{d}}=\left[\begin{array}[]{cc}\bm{\Sigma}_{b_{gd}}&\bm{\mathrm{0}}\\ \bm{\mathrm{0}}&\bm{\Sigma}_{b_{ad}}\end{array}\right].

In addition, given the visual measurement model in Section V-D, the reprojection error term is:

[TABLE]

VI tightly-coupled monocular VISLAM

In this section, we introduce our tightly-coupled visual-inertial monocular SLAM approach, which combines the EKF-based VIO front-end with the BA and loop closure back-end to provide the accurate and robust state estimation. An overview of our system is shown in Fig.1. When new inertial sensor is used, we firstly perform the EKF VIO alone with initial bias value of zero to obtain a good initial bias estimates before starting the all system, then set it as the initial bias value of our system.

VI-A Map Initialization

The map initialization is in charge of constructing an initial set of map points for the subsequent nonlinear optimization and loop closure. The initial map is created according to the estimated state from the EKF VIO. Since the accuracy of the initially created map will pose big influence on the accuracy of the whole system, we create initial map after running the EKF system about 10 seconds to converge.

Firstly, we extract ORB features in the current frame $k$, and search for feature matches with the reference frame $r$. If sufficient feature correspondences are found, we perform the next step, else set the current frame as reference frame. The second step is to check the parallax of each correspondence, and pick out a set of feature matches $\mathcal{F}$ which have sufficient parallax. When the size of $\mathcal{F}$ is greater than a threshold, using the estimated pose from the EKF based VIO, we triangulate the matched features $\mathcal{F}$. Then, if enough map points are successfully created, a full BA combining reprojection error term and temporal IMU error term is applied to refine the initial map. Finally, using the optimized state and map to correct the EKF state according to the feedback mechanism described in Section VI-D.

VI-B FRONT-END

VI-B1 State Estimation

We firstly perform the VIO described in Section IV to estimate the state of current frame $\bm{\mathrm{X}}_{B_{k}}=\{{\bm{\mathrm{\xi}}^{W}_{B_{k}}}\ {\bm{\mathrm{p}}^{W}_{B_{k}}}\ {\bm{\mathrm{v}}^{W}_{B_{k}}}\ {\bm{\mathrm{b}}_{a_{k}}}\ {\bm{\mathrm{b}}_{g_{k}}}\}$. In order to ensure the real-time performance, we only remain those features visible in current frame to limit the number of keypoints in the state vector. In our work, we maintain 50 features in the EKF state. Finally, if the map is updated in the back-end, the state of the current frame will be corrected according to the feedback mechanism described in Section VI-D. In this way, the front-end can always provide reliable state estimates even if we run for long periods of time.

VI-B2 Keyframe Selection

After the state of a frame is estimated by the VIO, we adopt three criteria to determine whether this frame is a keyframe. (1) The time interval from last keyframe is beyond a certain threshold. This criteria ensures the accuracy of the system. Because the IMU just provide valuable information in short-time, if the time interval from the last keyframe is too long, the IMU constraint between two keyframes will become inaccurate. (2) The nonlinear optimization in back-end is finished. This criteria makes as many keyframes as possible to enhance the motion tracking accuracy. (3) The rotation angle from last keyframe is beyond a certain threshold. This criteria ensures the reconstruction of the globally consistent map.

Once a frame is selected as a keyframe, we extract the ORB features for the new keyframe and trigger the back-end to make the pose estimation more accurate.

VI-C BACK-END

Once a new keyframe is inserted, the nonlinear optimization described in Section V is performed to optimization the local map in a parallel thread. After the local BA is finished, some redundant keyframes will be culled to make the factor graph more concise. In loop closure thread, place recognition is performed. Once a loop is detected, a $\bm{\mathrm{Sim}}(3)$ optimization and a full BA is performed to eliminate the accumulated drift. We refer the interested readers to papers [17][28] for more details.

VI-D FEEDBACK MECHANISM

The EKF-based VIO can estimate the frame state efficiently. However since the accumulation of the linearization errors and the absence of the loop closure, the error of the state estimate will accumulate as time goes on. If the state provided in the front-end drift too much, the local BA will be hard to find best estimates. Therefore in order to constrain the error of the state estimated from VIO, we provide following feedback mechanism. The feedback mechanism is invoked whenever the map in the back-end is updated and contains two steps:

VI-D1 State Correction

After the BA and loop closure are performed in the back-end, the state estimation of the last keyframe $i$ is accurate enough. The key observation of the state correction is to improve the state estimation of the current frame $k$ by leveraging the optimized state of the last keyframe $i$. Therefore, given the estimated state $\bm{\mathrm{X}}_{B_{k}}$ from VIO, we optimize the state of the current frame by performing the nonlinear optimization as shown in Fig. 3, that is minimizing the following objective function:

[TABLE]

where $F_{k}$ denotes the features matched with the map points in current frame $k$. The form of $\|{\bm{\mathrm{r}}_{\mathcal{I}_{ik}}}\|^{2}_{{\bm{\Sigma}}_{\mathcal{I}_{ik}}}$ and $\|{\bm{\mathrm{r}}_{\mathcal{C}_{kl}}}\|^{2}_{{\bm{\Sigma}}_{\mathcal{C}_{kl}}}$ is the same as (36) and (37) respectively. Then the optimized state $\bm{\mathrm{X}}_{B_{k}}^{\ast}$ and its covariance matrix $\bm{\mathrm{\Upsilon}}_{B_{k}}$ obtained from the optimization are used to update the EKF state of the current frame as:

[TABLE]

VI-D2 EKF Map Correction

It is well known that the estimated pose relative to the map in the state vector is not drift. Therefore, as long as the position of map points in the EKF state vectorF is consistent with the optimized global map, the accuracy of the estimated state will accordingly increase. Therefore after the IMU state in EKF state vector is updated by state correction, we will add features in the optimized consistent map to the state vector. We denote $F_{k}^{{}^{\prime}}$ as a set of features matched with the optimized map in current frame $k$, the outliers are removed based on the optimized state. If we need to add $n$ new features to the EKF state, we randomly select n features in $F_{k}^{{}^{\prime}}$, then compute their initial position and add it to the filter state vector. For selected new feature $l$, the initial position ${\bm{\mathrm{p}}_{l}^{W}}^{\ast}={[x_{l}^{\ast}\ y_{l}^{\ast}\ z_{l}^{\ast}\ \theta_{l}^{\ast}\ \phi_{l}^{\ast}\ \rho_{l}^{\ast}]}^{T}$ is set as follows. ${[x_{l}^{\ast}\ y_{l}^{\ast}\ z_{l}^{\ast}]}^{T}$ is computed as (22) and ${[\theta_{l}^{\ast}\ \phi_{l}^{\ast}]}^{T}$ is computed as (23)(24) using the updated state $\bm{\mathrm{X}}_{{k|k}}^{\ast}$ from (39) and the observation of the selected new features in the current frame. In addition, for computing $\rho_{l}^{\ast}$, we firstly transform the map point in world frame $\bm{\mathrm{L}}_{l}^{W}$ to the current camera frame $\bm{\mathrm{L}}_{l}^{C_{k}}$:

[TABLE]

Then the initial inverse depth is obtained by $\rho_{l}^{\ast}=\frac{1}{\|\bm{\mathrm{L}}_{l}^{C_{k}}\|}$, and its variance is set as follows:

[TABLE]

where $\bm{\mathrm{J}}_{L}=-\frac{1}{\|\bm{\mathrm{L}}_{l}^{C_{k}}\|^{3}}{\bm{\mathrm{L}}_{l}^{C_{k}}}^{T}{({\bm{\mathrm{R}}^{W}_{B_{k}}}^{\ast}\bm{\mathrm{R}}^{B}_{C})}^{T}$ and $\bm{\mathrm{J}}_{Rt}=-\frac{1}{\|\bm{\mathrm{L}}_{l}^{C_{k}}\|^{3}}{\bm{\mathrm{L}}_{l}^{C_{k}}}^{T}\left[{\bm{\mathrm{R}}^{B}_{C}}^{T}\left({\bm{\mathrm{R}}^{W}_{B_{k}}}^{\ast T}(\bm{\mathrm{L}}_{l}^{W}-{\bm{\mathrm{p}}^{W}_{B_{k}}}^{\ast})\right)^{\wedge}\ -{\bm{\mathrm{R}}^{B}_{C}}^{T}\right]$ are the Jacobians of $\rho_{l}^{\ast}$ with respect to $\bm{\mathrm{L}}_{l}^{W}$ and $[{\bm{\mathrm{R}}^{W}_{B_{k}}}^{\ast}\ {\bm{\mathrm{p}}^{W}_{B_{k}}}^{\ast}]$ on manifold. Besides, $\bm{\Sigma}_{L}$ and $\bm{\Sigma}_{Rt}$ are the covariance matrix for $\bm{\mathrm{L}}_{l}^{W}$ and $[{\bm{\mathrm{R}}^{W}_{B_{k}}}^{\ast}\ {\bm{\mathrm{p}}^{W}_{B_{k}}}^{\ast}]$, which is computed from the BA in the back-end.

VII EXPERIMENTS

We make a complete evaluation of the proposed algorithm qualitatively and quantitatively on the EuRoC dataset[29]. The dataset contains 11 data sequences, which was recorded from a flying MAV in two different 30$m^{2}$ indoor rooms and a $300m^{2}$ industrial environment. Depending on the illumination, texture and motion dynamics, the data sequences are classified as easy, medium and difficult levels. The dataset provides synchronized global shutter WVGA stereo images at 20Hz, IMU measurements at 200Hz and ground truth state at 200Hz. We only use images from the left camera. Firstly, we evaluate the proposed algorithm qualitatively and quantitatively on the EuRoC dataset to show the accuracy of our system. Then we compare our method with other state-of-the-art approaches on the EuRoC dataset. Finally, the performance of our algorithm is validated again by indoor real-world experiments using the sensor of Intel RealSense ZR300. The experiments are performed on a laptop with Intel Core i5 2.2GHz CPU and an 8GB RAM.

VII-A Algorithm Evaluation

We successfully perform our algorithm on all 11 sequences of EuRoC dataset in real-time. Fig. 4 and Fig. 6 show the comparisons of the ground truth, the trajectory estimated by our algorithm and EKF VIO on V1 02 medium and MH 02 easy sequence respectively. The corresponding x,y,z translation error versus time is shown in Fig. 5 and Fig. 7. The estimated trajectories are aligned with the ground truth using the method of Horn[30]. As evident, the translation error of our approach is smaller than the error of EKF VIO, thus proving the superiority of our algorithm towards the EKF VIO method.

For quantitative analysis, table I shows the translation Root Mean Square Error(RMSE) of the estimated trajectory for each sequence, as proposed in [31]. The proposed method has achieved the average translation RMSE of 0.082m, 0.056m and 0.053m for V1, V2 and MH sequences with respect to 0.186m, 0.213m, 0.319m of EKF VIO system, which illustrates that our method reduced the error of 55%, 73%, 83% for V1, V2 and MH sequences. From the third and fourth columns of the table I, we can know that adding BA and feedback mechanism to EKF VIO system, the RMSE of sequences are much reduced, this is since (1) BA is able to relinearize measurement models to properly deal with the nonlinearity of the system, (2) increasing the number of feature matches in the window of local BA can greatly improve the accuracy, and (3) feedback of the optimized state and map to VIO can improve the consistency and accuracy of the EKF system. However, the advantage of BA and feedback mechnism is not well demonstrated in V1 02 medium, V1 03 difficult and V2 03 difficult sequences, because in which fast rotation and low texture are frequently happened. Both fast rotation and low texture result in fewer feature correspondences, and thereby make less constraints in local BA, thus the accuracy of the system is not improved greatly. Therefore, in these sequences, adding a loop closure achieved much better accuracy by eliminating the accumulated error when revisiting an already mapped area.

While adding BA, loop closure and feedback mechanism improves accuracy of the system, it increases the computational cost of the system due to the need of (1) extracting a lot of features for BA and loop closure, and (2) pose optimization and EKF update of the feedback mechanism. The computational increase only occurs in selected keyframes, so the proposed method will not incur too much computational cost. Our algorithm requires approximately $27$ msec for processing each image, so it can still run in real time, at about 40Hz.

VII-B Comparison to State-of-the-art Algorithms

We compare the proposed method with the state-of-the-art VINS-MONO [19][20], OKVIS [14], ROVIO [32] and ORB-VISLAM[17] method. VINS-MONO, OKVIS and ROVIO are open-source and contain the default parameters for the EuRoC dataset, for fair comparison, we only use left image. ORB-VISLAM shows its results in EuRoC dataset, which allowing for a direct comparison.

A comparison of the translation RMSE of the estimated trajectories on EuRoC dataset are shown in table II, X means the concerned method fails to run in the sequence. From these results, we can draw following conclusions. ORB-VISLAM and our algorithm have obtained the best accuracy, this is because local BA in both methods is performed in a parallel thread, so more feature correspondences are used in local BA. Besides, both methods can close loop to eliminate the accumulated error. Howerver, ORB-VISLAM fails to track the V1 03 difficult sequence. In comparison, VINS-MONO and OKVIS perform BA in the thread of tracking, so the number of features contained in local BA must be limited to ensure the real-time performance, therefore leading to slightly worse accuracy. VINS-MONO with loop can achieve better accuracy than OKVIS, due to its capability to close loop. In addition, ROVIO is an EKF method and not able to close loop, its the linearization error and the error accumulation make the method obtain diminished localization error. However, since it is a direct method, it can achieve minor drift in fast motion sequences of V1 03 difficult and V2 03 difficult.

For sequence MH 04 difficult, the estimated trajectories are shown in Fig. 8, and translation errors versus time are shown in Fig. 9. In this sequence, our loop closure is not triggered, however in the error plot, our method still achieved smallest translation error, which can prove the superiority of our algorithm again.

VII-C Indoor Real-world Experiment

We perform the indoor experiment in an 60$m^{2}$ office environment using the monocular-inertial Realsense ZR 300 sensor suite that provides images at the frequency of 20 Hz and IMU measurements at 200 Hz. As shown in the accompanying video111https://youtu.be/5_G8jUOjtN0, we hold the sensor suite by hand and walk in normal pace in the office, and we starts and ends at the same location. Fig. 10 shows the estimated trajectory, from which we can know there is no noticeable drifts occurred when we circle indoor. The end-to-end error is 0.055m with respect to the total length of 82m, it is just the 0.067% of the total trajectory length.

VIII Conclusion and future work

In this paper, we have presented a tightly-coupled monocular VISLAM system, which robustly tracks camera motion by EKF VIO, and perform non-linear optimization and loop closure to solve the linearization issues of the EKF system and eliminate the accumulated error. We also proposed a feedback mechanism to directly improve the consistency and accuracy of the EKF system. Therefore, our algorithm has achieved high accuracy, the performance of the proposed method is validated through experiments.

Point feature-based monocular VISLAM is prone to fail in poorly textured scenes or motion blurred images. Therefore in the future, we inted to deal with these specific situations for better accuracy and robustness. We also aim to build dense map to assist the understanding of the environment.

Acknowledgment

This paper is supported by National Science Foundation of China[grant number 61375081]; a special fund project of Harbin science and technology innovation talents research [grant number RC2013XK010002].

-A Matrices $\bm{\mathrm{\Phi}}_{k}$ and $\bm{\mathrm{G}}_{k}$

The Jacobian matrix $\bm{\mathrm{\Phi}}_{k}$ in (13) is:

[TABLE]

where $\bm{\mathrm{I}}_{3\times 3}$ is the $3\times 3$ identity matrix, and

[TABLE]

In addition, the Jacobian matrix $\bm{\mathrm{G}}_{k}$ in (13) is:

[TABLE]

with $\bm{\mathrm{\Phi}}_{\xi n_{g}}=\bm{\mathrm{\Phi}}_{\xi b_{g}}$ and $\bm{\mathrm{\Phi}}_{vn_{a}}=\bm{\mathrm{\Phi}}_{vb_{a}}$.

-B Measurement matrices

The Jacobian of the measurement model with respect to the IMU state in (18) is represented as:

[TABLE]

with:

[TABLE]

where $\bm{\mathrm{f}}_{l}^{C_{k}}={[x_{l}^{C_{k}}\ y_{l}^{C_{k}}\ z_{l}^{C_{k}}]}^{T}$. In addition, the Jacobian of the measurement model with respect to the $l^{th}$ feature position in (18) is:

[TABLE]

with:

[TABLE]

-C State Augmentation Jacobian

The Jacobian matrix used to augment the covariance matrix of the state vector is:

[TABLE]

with:

[TABLE]

where:

[TABLE]

and

[TABLE]