A Discrete-Time Attitude Observer on SO(3) for Vision and GPS Fusion

TL;DR

This paper introduces a discrete-time geometric attitude observer that fuses monocular vision and GPS data to improve orientation estimation accuracy and stability, effectively reducing drift and directly estimating North-East-Down frame orientation.

Contribution

It presents a novel attitude observer with rigorous stability analysis that fuses vision and GPS data, ensuring exponential convergence and robustness to large initial errors.

Findings

Observer effectively reduces visual odometry drift

Successfully estimates North-East-Down orientation

Converges exponentially to true attitude

Abstract

This paper proposes a discrete-time geometric attitude observer for fusing monocular vision with GPS velocity measurements. The observer takes the relative transformations obtained from processing monocular images with any visual odometry algorithm and fuses them with GPS velocity measurements. The objectives of this sensor fusion are twofold; first to mitigate the inherent drift of the attitude estimates of the visual odometry, and second, to estimate the orientation directly with respect to the North-East-Down frame. A key contribution of the paper is to present a rigorous stability analysis showing that the attitude estimates of the observer converge exponentially to the true attitude and to provide a lower bound for the convergence rate of the observer. Through experimental studies, we demonstrate that the observer effectively compensates for the inherent drift of the pure monocular vision based attitude estimation and is able to recover the North-East-Down orientation even if it is initialized with a very large attitude error.

Figures

Equations

$$R_{C_{k+1}}^{C_0}$$

$$p_{C_{k+1}}^{C_k}=\frac{1}{d}{R}_k^{\top }(p_{k+1}^A-p_k^A)$$

$$p_{k+1}^A-p_k^A={\int }_{t_k}^{t_{k+1}}{v}^A(s)\text{d}s\approx (t_{k+1}-t_k){\bar{v}}_k^A.$$

$${\bar{p}}_k^A$$

$${\bar{p}}_k^C$$

$${\bar{p}}_k^C=R_k^{\top }{\bar{p}}_k^A.$$

$$R_{k+1}=R_k{R}_{C_{k+1}}^{C_k},R_0=R_{C_0}^A.$$

$${\hat{R}}_{k+1}={\hat{R}}_k\exp (-{\hat{R}}_k^{\top }{\text{grad}}_1\varphi ({\hat{R}}_k,{\bar{p}}_k^C))R_{C_{k+1}}^{C_k},$$

$$\varphi ({\hat{R}}_k,{\bar{p}}_k^C)=({\hat{R}}_k{\bar{p}}_k^C-{\bar{p}}_k^A{)}^{\top }L({\hat{R}}_k{\bar{p}}_k^C-{\bar{p}}_k^A)$$

$$\displaystyle\text{grad}_{1}\!\phi(\hat{R}_{k},\bar{p}^{C}_{k})=-\Big{(}\!\!\big{(}L(\hat{R}_{k}\bar{p}^{C}_{k}-\bar{p}^{A}_{k})\!\big{)}_{\!\times}\hat{R}_{k}\bar{p}^{C}_{k}\Big{)}_{\!\times}\hat{R}_{k}.$$

$$\displaystyle\hat{R}_{k+\!1}\!\!=\!\exp\!\Big{(}\!\!\Big{(}\!\!\big{(}L(\hat{R}_{k}\bar{p}^{C}_{k}-\bar{p}^{A}_{k})\!\big{)}_{\!\times}\hat{R}_{k}\bar{p}^{C}_{k}\Big{)}_{\!\times}\!\Big{)}\hat{R}_{k}R^{{C}_{k}}_{{C}_{k+\!1}},~{}\hat{R}_{0}\!\!=\!\!I_{3}.$$

$$E_k={\hat{R}}_k{R}_k^{\top }.$$

$$E_{k+1}$$

$$\displaystyle=\exp\!\Big{(}\Big{(}\!\!\big{(}L(E_{k}\bar{p}^{A}_{k}-\bar{p}^{A}_{k})\!\big{)}_{\!\times}E_{k}\bar{p}^{A}_{k}\Big{)}_{\!\times}\Big{)}E_{k}.$$

$$\frac{1}{T+1}\sum _{i=k}^{k+T}(I_3-{\bar{p}}_i^A{\bar{p}}_i^A{}^{\top })\ge \beta I_3,$$

$$\Vert I_3-E_k{\Vert }_F\le \Vert I_3-E_0{\Vert }_F{\alpha }^k,$$

$${ϵ}_{k+1}={ϵ}_k-l{P}_k{ϵ}_k,$$

$$P_k:=-{\bar{p}}_k^A{}_{\times }{\bar{p}}_k^A{}_{\times }=I_3-{\bar{p}}_k^A{\bar{p}}_k^A{}^{\top }$$

$${\mathcal{L}}_k={ϵ}_k^{\top }{ϵ}_k.$$

$${\mathcal{L}}_{k+1}-{\mathcal{L}}_k$$

$$=-\gamma {ϵ}_k^{\top }{P}_k{ϵ}_k.$$

$${\mathcal{L}}_{k+1+T}-{\mathcal{L}}_k=-\gamma \sum _{i=k}^{k+T}{ϵ}_i^{\top }{P}_i{ϵ}_i$$

$$\sum _{i=k}^{k+T}{ϵ}_i^{\top }{P}_i{ϵ}_i\ge \frac{1}{2}\sum _{i=k}^{k+T}{ϵ}_k^{\top }{P}_i{ϵ}_k-\sum _{i=k}^{k+T}({ϵ}_i^{\top }-{ϵ}_k^{\top })P_i({ϵ}_i-{ϵ}_k).$$

$$\sum _{i=k}^{k+T}({ϵ}_i^{\top }-{ϵ}_k^{\top })P_i({ϵ}_i-{ϵ}_k)\le \frac{l^{2}T(T+1)}{2}\sum _{i=k}^{k+T}{ϵ}_i^{\top }{P}_i{ϵ}_i$$

$$\sum _{i=k}^{k+T}{ϵ}_k^{\top }{P}_i{ϵ}_k\ge \beta (T+1){ϵ}_k^{\top }{ϵ}_k=\beta (T+1){\mathcal{L}}_k$$

$$\sum _{i=k}^{k+T}{ϵ}_i^{\top }{P}_i{ϵ}_i\ge \frac{\beta (T+1)}{2}{\mathcal{L}}_k-\frac{l^{2}T(T+1)}{2}\sum _{i=k}^{k+T}{ϵ}_i^{\top }{P}_i{ϵ}_i,$$

$${\mathcal{L}}_{k+1+T}\le \bar{\alpha }{\mathcal{L}}_k.$$

$${\mathcal{L}}_k\le {\mathcal{L}}_{n(1+T)}\le {\bar{\alpha }}^n{\mathcal{L}}_0,\forall k\ge n(1+T)$$

$$\Vert {ϵ}_k\Vert \le \Vert {ϵ}_0\Vert {\alpha }^k$$

$$p_{C_{k+1}}^{C_0}-p_{C_k}^{C_0}$$

$$p_{C_{k+1}}^{C_0}-p_{C_k}^{C_0}={R_{C_0}^A}^{\top }(p_{C_{k+1}}^A-p_{C_k}^A).$$

$${ϵ}_i-{ϵ}_k=\sum _{j=k}^{i-1}{ϵ}_{j+1}-{ϵ}_j=-l\sum _{j=k}^{i-1}{P}_j{ϵ}_j.$$

Full text

I Introduction

Attitude estimation is the problem of determining the orientation of a vehicle with respect to a known frame of reference. Due to its broad application in robotics, navigation, and control of mechanical systems, this problem has been extensively studied in the past decades [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. High-end IMUs are commonly used for pose estimation in applications, such as mining, that require long term accurate localisation. However, high-end IMUs are usually expensive, heavy, and power hungry. Cameras are attractive alternatives for pose estimation due to their low cost, small size, and high accuracy. Lots of elegant visual odometry algorithms have been developed in the literature and have been successfully used for localisation of moving robots (see e.g. [11, 12, 13, 14, 15] and the references therein). Monocular visual odometry algorithms rely on matching consecutive images to estimate the relative transformation of camera between those images. This relative transformation is concatenated over time to compute the absolute pose of the camera with respect to the first camera frame. Despite their popularity, there are two potential problems with pure visual odometry in applications that require long term accurate positioning; drift over time, and lack of a global known reference frame.

The concatenated estimate of camera motion based purely on the visual measurements inevitably drifts over time, due to uncertainties in camera intrinsic parameters, aggregation or numerical errors, and measurement noise. Although attenuated, the drift still exists even if graph based bundle adjustment techniques are used to optimally estimate the relative motion over tens of images [16]. One way of mitigating the drift is to perform loop closure [17]. However, not all robot paths necessarily contain loops in practice. Moreover, effective loop closure might be difficult in some applications, such as mapping of mine sites, where there is perceptual aliasing in the environment (i.e. several distinct places look similar), or where large moving vehicles, dust, wind, etc. change the appearance of the same location [17]. Also in real time applications, even if effective loop closure could be done, the estimate of robot’s path still drifts between two consecutive occurrence of the loop closure.

The lack of a global known reference frame is another problem in pure vision based navigation. Many control objectives require knowledge of the robot’s attitude with respect to a fixed known frame of reference (e.g. the North-East-Down (NED) frame) to navigate the robot toward a target point in the environment. Vision based orientation estimates can not be computed with respect to the NED frame as camera alone does not have any information about the NED frame. Although the concatenation of relative transformations yields an estimate of the orientation with respect to the initial camera frame, the transformation of the initial camera frame with respect to the NED frame is still unknown.

Fusing GPS measurements with vision helps mitigating both of the problems discussed above. GPS velocity measurements contain the heading information, that help compensating for the inherent drift of visual odometry estimates. Also, GPS measures the motion directly with respect to a known reference, enabling estimation of the attitude directly with respect to the NED frame. Fusing vision with GPS measurements is particularly of interest in open cut mine sites where RTK GPS is employed to provide centimetre-level position accuracy, albeit with the aid of reference stations. GPS measurements can be quite noisy and unreliable in deep mine pits, hence require fusion with high bandwidth sensors such as IMUs in Inertial Navigation Systems (INS). Fusing the low bandwidth GPS measurements with the high bandwidth visual measurements is a promising solution for long term localisation with a comparable level of accuracy to a high-end INS at a fraction of the cost (see section VII).

Recursive stochastic filtering techniques, such as (extended) Kalman filtering, are very popular in sensor fusion for navigation [18, 1]. Nevertheless, Kalman filters usually require careful tuning, lack rigorous proof of stability when applied to non-linear systems, and are prone to divergence if not initialised carefully [19, 20, 3, 21, 22, 23]. To address these issues, geometric nonlinear observer design methods on Lie groups have been developed in the past decade [24, 21, 25, 26]. These methods do not provide as rich information as stochastic methods due to the lack of the state covariance. Nevertheless, continuous-time geometric nonlinear attitude observers developed based on these methods are simple to tune, computationally cheap, and have demonstrated good performance in applications where high rate IMUs are present [3, 4, 5, 6, 7, 8, 9, 10]. These methods, however, can not be directly applied to the vision-GPS fusion problem as the underlying kinematics of the this problem is naturally discrete-time due to the low rate of camera in most applications.

Elegant methods for designing discrete-time observers on Lie groups are developed in [27, 28]. Although suitable for satellite applications, the observer of [28] relies on attitude dynamics and the knowledge of external forces applied to the rigid body, hence is not applicable to the problem considered here. The observer of [27, equ. (18)] assumes availability of at least two non-collinear vectorial measurements and provides a novel stability analysis. As we will explain in this paper, however, GPS velocity measurements provide only a single vectorial measurement which prevents application of the stability theory of [27] to the particular problem considered here. To the best of our knowledge, there is no discrete-time geometric attitude observer in the literature concerning the problem of fusing vision and GPS velocity measurements (without the aid of IMU).

In this paper, we propose a discrete-time geometric observer on $\text{SO}(3)$ for fusing monocular visual odometry estimates with GPS velocity measurements. The structure of the proposed observer is inspired by [27] and the proof of stability provided here is an extension to that work based on inspirations from [22, 29, 5, 9]. We show that the estimates of relative translation of camera (provided by a visual odometry algorithm) can be modelled as a vectorial measurement whose associated reference direction is the GPS velocity measurement. Although this is the sole vectorial measurement available in this application, we show that the attitude can be asymptotically estimated if the direction of the velocity of vehicle varies sufficiently with time (see Assumption 1). A key contribution of this paper is to provide a rigorous stability analysis yielding a deep understanding of the underlying fundamental performance limits of the vision and GPS fusion problem in terms of the convergence rate of the observer. Our fusion method is independent of the choice of the visual odometry algorithm, giving the users an extra freedom to apply any algorithm that suites their particular application. We demonstrate the application of the developed theory by experimental studies showing excellent performance of the proposed observer in compensating for the drift in the visual odometry and recovering the orientation of vehicle with respect to the NED frame.

The structure of paper is as follows. In Sections II and III formulate, for the first time, the attitude estimation problem using vision and GPS velocity fusion as a discrete-time geometric observer design problem on $\text{SO}(3)$with an associated vectorial measurement belonging to the unit sphere. The attitude observer is proposed in Section IV and the stability analysis is performed in Section V (see Theorem 1). We provide simulation studies as well as experimental results in Sections VI and VII, respectively, and we finish the paper with conclusions in Section VIII.

II Problem formulation

Consider a vehicle equipped with a vision system and a GPS receiver moving in an environment with an unknown map. Assume that the GPS receiver and the camera are rigidly connected to the vehicle. Denote the North-East-Down (NED) reference frame (i.e. inertial frame) by $\{A\}$ and suppose that the reference frame $\{C\}$ is attached to the camera and moves with it. For simplicity, we refer to the camera reference frame when the $k$-th image was captured as the $k$-th camera frame. Denote the rotation of the $i$-th camera frame with respect to (wrt.) the $j$-th camera frame by $R^{{C}_{j}}_{{C}_{i}}\in\text{SO}(3)$. By matching corresponding feature points in two consecutive images $k$ and $k+1$, visual odometry algorithms are capable of computing the relative rotation between the camera frames of those two images [11, 12] We assume that $R^{{C}_{k}}_{{C}_{k+1}}$ is provided by a visual odometry algorithm, and we treat this variable as a known measurement. One can obtain the rotation of the $(k+1)$-th camera frame wrt. the initial camera frame $C_{0}$ by concatenating previous relative rotations, as follows.

[TABLE]

Visual odometry techniques are also able to compute the relative translation of the consecutive camera frames. For monocular cameras, however, such a relative translation is valid only up to a constant scale since the depth can not be measured in the monocular setup [30]. Denote the translation of the $i$-th camera frame wrt. the $j$-th camera frame (and expressed in the $j$-th camera frame) by $p^{{C}_{j}}_{{C}_{i}}\in\mathds{R}^{3}$. Assuming that the camera and GPS receiver are placed sufficiently close together111It is possible to extend the derivations of this paper to the case where the camera and the GPS receiver are not placed close together. We consider that as a possible future work., it is shown in Appendix A that

[TABLE]

where $p^{A}_{k+1}\in\mathds{R}^{3}$ is the $k$-th camera position (wrt. and expressed in the NED frame), $p^{{C}_{k}}_{{C}_{k+1}}$ is the relative translation of the $k$-th camera frame wrt. the $(k+1)$-th camera frames (expressed in the $k$-th camera frame) computed by the visual odometry algorithm and treated here as a known measurement, and $R_{k}\in\text{SO}(3)$ is the rotation of the $k$-th camera frame wrt. the NED frame. This rotation is related to $R^{{C}_{0}}_{{C}_{k}}$ via $R_{k}=R^{A}_{{C}_{0}}R^{{C}_{0}}_{{C}_{k}}$ where $R^{A}_{{C}_{0}}\in\text{SO}(3)$ is the unknown initial rotation of the camera frame wrt. the NED frame. One can use the GPS position measurements to compute $p^{A}_{k+1}-p^{A}_{k}$ in (2), nevertheless, numerical differencing the noisy GPS position measurements further increases the level of noise. For this reason, it is recommended to use the GPS velocity measurements to approximate $p^{A}_{k+1}$-$p^{A}_{k}$ as follows222Velocity measurements are computed in many GPS chips from the Doppler shift of the pseudo-range signals, yielding far more accurate velocity measurements compared to numerically differencing the range based position measurements. This accuracy is further enhanced if time-differenced carrier phase (TDCP) based GPS is used..

[TABLE]

If the vehicle’s velocity is approximately constant from $t_{k}$ to $t_{k+1}$, we simply use $\bar{v}_{k}^{A}=v_{k}^{A}\in\mathds{R}^{3}$ where $v_{k}^{A}$ is the GPS velocity measurement at time $t_{k}$. Otherwise, if the velocity changes linearly with time from $t_{k}$ to $t_{k+1}$, we use the approximation $\bar{v}_{k}^{A}=0.5(v_{k+1}^{A}+v_{k}^{A})$.

The aim of this paper is to estimate $R_{k}$ using the discrete-time dynamics model (1) and the output measurement model (2). That is to fuse the outputs of visual odometry algorithm $R^{{C}_{k}}_{{C}_{k+1}}$ and $p^{{C}_{k}}_{{C}_{k+1}}$ with the GPS velocity measurement $v_{k}$ (or equivalently the difference of consecutive position measurements if the velocity is unavailable). There are several challenges in such an estimation problem, some of which have been listed below.

• •

As opposed to classical linear estimation problem where the state space is a vector space, here the underlying system (1) evolves on the manifold $\text{SO}(3)$.

• •

The underlying dynamics is discrete-time, which prevents direct application of previous continuous-time observer design techniques on Lie groups [24, 21, 25, 26].

• •

The scale $d$ of the translation measurements provided by the visual odometry algorithm is unknown.

• •

The GPS measurements do not immediately provide information about the camera rotation. More specifically, the GPS measurements are related to the camera orientation via (2) where the initial camera orientation $R^{A}_{{C}_{0}}$ is unknown. Also, since (2) only provides a single $3$D vectorial measurement, the rotation is not instantaneously observable.

Note that we do not consider stochastic measurement noise in our modeling as our design framework is based on the deterministic geometric setup [24, 21, 25, 26].

III Vision-GPS fusion as a geometric observer design problem

In this section, we transform the underlying system kinematics (1) and the output model (2) to a form more suitable for observer design purpose. We normalize the measurements involved in (2) to eliminate the unknown scale $d$. Defining

[TABLE]

we have

[TABLE]

The variables $\bar{p}^{C}_{k}$ and $\bar{p}^{A}_{k}$ are unit vectors expressed in the camera frame and the NED frame, respectively. The normalization is particularly important since the scale of the visual odometry estimates drifts over time. The normalization not only isolates the effect of the scale drift, but also eliminates the need for knowledge of the velocity time stamp $t_{k}$. Equation (6) formulates the output measurement model of the system as a right action of $\text{SO}(3)$(the state space) on $\text{S}^{2}$ (the output manifold) where $\bar{p}^{A}_{k}$ represent a time-varying reference output in this context [24, 31, 22, 29]. This model directly relates the measurements $\bar{p}^{C}_{k}$ and $\bar{p}^{A}_{k}$ to the unknown state $R_{k}$ and is used in the observer design.

Multiplying the sides of (1) by the constant rotation $R^{A}_{{C}_{0}}$ from the left and recalling $R_{k}=R^{A}_{{C}_{0}}R^{{C}_{0}}_{{C}_{k}}$, we have

[TABLE]

Now, the observer design problem is simplified to estimate $R_{k}$ governing the dynamics (7) using the measurement of $R^{{C}_{k}}_{{C}_{k+1}}$ (as the input of the dynamics (7)) and the vectorial measurements $\bar{p}^{C}_{k}$ and $\bar{p}^{A}_{k}$ modeled via (6). The main advantage of formulating the problem as above compared to Section II is that the initial camera rotation $R^{A}_{{C}_{0}}$ now appear as an unknown initial condition of the estimation problem rather than as unknown variables in the system equation.

IV Observer design on $\text{SO}(3)$

The underlying kinematics of the system (7) is left-invariant and the considered output model (6) is right-invariant (in the sense of [24, 31]). Inspired by the novel work of [27], we consider the following structure for our observer.

[TABLE]

with $\hat{R}_{0}=I_{3}$, where $\phi:\text{SO}(3)\times\text{S}^{2}\to\mathds{R}^{+}$ is a right-invariant cost function (in the sense of [24]) and $\text{grad}_{1}\!\phi(\hat{R}_{k},\bar{p}^{C}_{k})\in T_{\hat{R}_{k}}\text{SO}(3)$ is the gradient of $\phi(.,\bar{p}^{C}_{k}):\text{SO}(3)\to\mathds{R}^{{}^{+}}$ wrt. a right-invariant Riemannian metric (see [24, Proposition 6.1.]). The motivation behind the observer structure (8) is that it generates an autonomous gradient-like estimation error dynamics if the reference output is constant [27]. We naturally choose the simple cost function

[TABLE]

where $L\in\mathds{R}^{3\times 3}$ is a constant positive definite matrix. Observe that $\hat{R}_{k}=R_{k}$ implies $\phi=0$. Employing standard differential geometry derivations, the gradient of $\phi$ wrt. the induced Riemannian metric $\langle\Omega_{1}R,\Omega_{2}R\rangle_{R}=\text{tr}(\Omega_{2}^{\top}\Omega_{1})$ is given by [32, Section VI]

[TABLE]

where $(.)_{\times}$ denotes the isomorphism from $\mathds{R}^{3}$ to $\mathfrak{so}(3)$ such that $a_{\times}b$ yields the vector cross product $a\times b$ for any $a,b\in\mathds{R}^{3}$. Using (8) and (10), and employing the property $\exp(\hat{R}^{\top}a_{\times}\hat{R})=\hat{R}^{\top}\exp(a_{\times})\hat{R}$, we propose the observer

[TABLE]

Remark 1

The observer (11) assumes GPS measurements are synchronized with images. In practice, the sampling rate of GPS is usually lower than that of the camera. In this case, one can introduce the auxiliary state $\hat{R}_{k}^{+}$ and break the observer (11) into two stages; prediction and update. In the prediction stage, the auxiliary state is updated by $\hat{R}_{k+1}^{+}=\hat{R}_{k}R^{{C}_{k}}_{{C}_{k+\!1}}$ each time a new relative orientation $R^{{C}_{k}}_{{C}_{k+\!1}}$ is provided by visual odometry. In the update stage, the state $\hat{R}$ is updated via \hat{R}_{k+\!1}=\exp\!\Big{(}\!\!\Big{(}\!\!\big{(}L(\hat{R}_{k}^{+}\bar{p}^{C}_{k}-\bar{p}^{A}_{k})\!\big{)}_{\!\times}\hat{R}_{k}^{+}\bar{p}^{C}_{k}\Big{)}_{\!\times}\!\Big{)}\hat{R}_{k}^{+} when GPS provides a new measurement. It is possible to show that the stability results of Section V still hold when the two-stage implementation of the observer is employed.

V Stability analysis

Although the structure of the proposed observer (11) is inspired by the novel work of [27], the stability analysis presented in [27, Proposition 4] does not apply here as the cost function (9) does not have a unique global minimum. Instead, we provide a new stability analysis by extending the results of [22, 29]. Consider the attitude estimation error.

[TABLE]

To prove $\hat{R}_{k}\to R_{k}$, we should show $E_{k}\to I_{3}$. Using (7), (11), (6), and (12), the dynamics of $E_{k}$ are given by

[TABLE]

One can verify that $E_{k}=I_{3}$ is an equilibrium point of the error dynamics (13). If $\bar{p}^{A}_{k}$ were constant and if we had two or more vector measurements, the elegant stability analysis of [27, Proposition 4] would be applied. However, only one vector measurement is available here. Since $\bar{p}^{A}_{k}$ is time-varying, the stability analysis is tedious even in the continuous-time case [22, 29]. Inspired by the continuous-time case, here we propose a novel stability analysis for the discrete-time case. The following assumption formally formulates the time-varying property of $\bar{p}^{A}_{k}$.

Assumption 1** (Persistency of excitation)**

There exist a positive integer $T$ and a constant $\beta>0$ such that

[TABLE]

for all non-negative integers $k$. $\square$

Assumption 1 could be written in other equivalent forms, see [22, Lemma 1]. This assumption is satisfied if $\bar{p}^{A}_{k}$ varies with time, i.e. if the vehicle’s direction of motion is not constant all the time (see [22, Lemma 1]). The only case where this assumption does not hold is when the vehicle’s trajectory is a perfectly straight line. It is worth mentioning that in the closely related problem of velocity-aided attitude estimation where an IMU is present instead of a camera, a similar condition as (14) is required unless an extra sensor (such as magnetometer) is available [33]. The more the direction of the motion of the vehicle changes, the larger $\beta$ gets. Nevertheless, since $\bar{p}^{A}_{k}$ is a unit vector, the largest singular value of the matrix $I_{3}-\bar{p}^{A}_{i}\bar{p}^{A}_{i}{}^{\top}$ is one and we have $\beta<1$ [22]. The following theorem summarizes the stability of the error dynamics (13). To simplify presentation of the results, we consider the observer gain $L=lI_{3}$ where $l$ is a scalar. It is straight-forward to extend the results to a more general positive definite gain matrix.

Theorem 1

Consider the observer (11) for the system (7) with the output (4)-(6). Suppose that Assumption 1 holds. Assume moreover that $L=lI_{3}$ where $0<l<2$. Then, the error dynamics (13) is locally exponentially stable to $E_{k}=I_{3}$. Moreover, there exists a compact region around $E_{k}=I$ in which

[TABLE]

where 0<\alpha:=\Big{(}1-\frac{\beta(T+1)(2l-l^{2})}{2+l^{2}T(T+1)}\Big{)}^{\frac{1}{2(1+T)}}<1 indicates a lower bound for the convergence rate and $\|I_{3}-E_{k}\|_{F}^{2}=\text{tr}((I_{3}-E_{k})^{\top}(I_{3}-E_{k}))=2\text{tr}(I_{3}-E_{k})$ is the induced Frobenius norm on $\text{SO}(3)$. $\square$

Proof of Theorem 1: the proof is mainly inspired by [22, Section V] which investigates the continuous-time scenario. There are subtle technical differences in the discrete-time case that prevent direct application of the continuous-time case. For this reason, here we provide a stability proof for the discrete-time case in detail. The first order approximation of the error $E_{k}$ around the identity is given by $E_{k}\approx I+{\epsilon_{k}}_{\times}$ where $\epsilon_{k}\in\mathds{R}^{3}$. Substituting for this approximation into the error dynamics (13) and resorting to Appendix B, we obtain the linearized error dynamics

[TABLE]

where

[TABLE]

is a time-varying positive semi-definite matrix whose rank is two. The Assumption 1 imposes a persistency of excitation condition on $P_{k}$ which can increase the rank of its summation over time to three. Consider the following Lyapunov candidate.

[TABLE]

Using (16) and noting that $P_{k}^{\top}=P_{k}$ and $P_{k}^{2}=P_{k}$ (due to (17) and the fact that $\bar{p}^{A}_{k}$ has a unit norm), we have

[TABLE]

where $\gamma:=2l-l^{2}$. Since $0<l<2$, we have $0<\gamma<1$. Since $P_{k}$ is positive semi-definite, we have $\mathcal{L}_{k+1}-\mathcal{L}_{k}\leq 0$ which indicates that the Lyapunov candidate is non-increasing along the system trajectories. Summing the sides of (20) from $k$ to $k+T$ yields

[TABLE]

Since $P_{i}\geq 0$, there exist $N_{i}\in\mathds{R}^{2\times 3}$ (with rank 2) such that $P_{i}=N_{i}^{\top}N_{i}$. Choose $a=N_{i}\epsilon_{k}$ and $b=N_{i}(\epsilon_{i}-\epsilon_{k})$ and employ the property $(a+b)^{\top}(a+b)\geq\frac{1}{2}a^{\top}a-b^{\top}b$ to obtain

[TABLE]

Using the derivations of Appendix C, we have

[TABLE]

Also, (14) and (18) yield

[TABLE]

Substituting for (23) and (24) into (22), we obtain

[TABLE]

which implies $\sum_{i=k}^{k+T}\epsilon_{i}^{\top}P_{i}\epsilon_{i}\geq=\frac{\beta(T+1)}{2+l^{2}T(T+1)}\mathcal{L}_{k}$. Substituting this into (21) yields

[TABLE]

where $\bar{\alpha}:=1-\frac{\beta(T+1)(2l-l^{2})}{2+l^{2}T(T+1)}$. Recalling that $0<l<2$ and $T\geq 1$, it is straight-forward to show that $1-\frac{\sqrt{2}}{2}\beta<\bar{\alpha}<1$. Employing (17) and recalling that $\bar{p}^{A}_{k}$ is a unit vector, we conclude that the largest singular value of $P_{k}$ is one. This, together with (14) implies that $0<\beta\leq 1$. Hence, $0<\frac{\sqrt{2}}{2}<\bar{\alpha}<1$. Taking $k=n(1+T)$, where $n$ is a non-negative integer, and successively employing (26) yields

[TABLE]

Hence, $\mathcal{L}_{k}$ decays exponentially to zero as $n$ goes to infinity. Substituting $\mathcal{L}_{k}$ for (18) and considering $n=\frac{k}{1+T}$ yields

[TABLE]

where $\alpha$ is defined in the statement of Theorem 1. Consequently, the equilibrium point $\epsilon_{k}=0$ of the linearized error dynamics (16) is uniformly globally exponentially stable [34]. Hence, by extending the results of [35, Theorem 4.15] to the discrete-time case according to the guidelines of [34], one concludes that the equilibrium point $E_{k}=I$ of the error dynamics (13) is uniformly locally exponentially stable. The exponential bound (15) follows from (28) directly. $\blacksquare$

By Theorem 1, the convergence rate of the observer depends on the value of $\beta$ which in turn depends on the motion of the vehicle (equ. (14)). Having an analytical solution for the convergence rate (provide by Theorem 1) makes it possible to obtain the value of the observer gain $l$ that yields the fastest convergence rate. This can be done by considering a moving horizon of time and numerically computing $T$ and $\beta$ satisfying (14) for that horizon. Given $\beta$ and $T$, one can analytical compute the value of the observer gain $l$ that maximizes the convergence rate $\alpha$. This yields a time-varying gain as $T$ and $\beta$ should be updates as the time horizon changes. It is straight-forward to extend the stability results of Theorem 1 to a time-varying gain.

VI Simulation results

Theorem 1 proves local stability of the error dynamics around the identity, but it does not provide an estimate of the domain of attraction. As the initial orientation of the camera wrt. to the NED frame is unknown, one cannot always choose the initial condition of the observer (11) close to the actual value of $R^{A}_{{C}_{0}}$. The aim of this section is to examine the domain of attraction of the error dynamics (13). Consider a camera moving with a constant velocity of $2\pi\approx 6.3$ (m/s) in a circular trajectory with the radius of $50$ (m) in the North-East plane. The NED frame is set to be at the centre of the circle. The camera starts its motion from the far North part of the path and moves clockwise to East. Assuming the sampling interval of $0.1$ seconds, the corresponding relative motions $R^{{C}_{k}}_{{C}_{k+1}}$ and $p^{{C}_{k}}_{{C}_{k+1}}$ are computed to model the output of visual odometry and passed to the observer (11) together with the linear velocity of the camera wrt. the NED frame. We neither model camera images nor we add noise to the measurements as the aim of this section is to only investigate the effect of initial condition on the attitude estimation error $\angle E_{k}=\text{acos}(1-\frac{1}{4}\|I-E_{k}\|_{F}^{2})$333The error $\angle E_{k}$ corresponds to the angle of rotation in the angle-axis decomposition of $E_{k}$ which models how far $\hat{R}_{k}$ is from $R_{k}$ [25, 22].. We perform $20$ simulations where in each simulation the observer is initialised with a random rotation corresponding to a random axis of rotation in a unit spare and a random angle of rotation between [math] and $179$ degrees. Fig. 1 shows the results of the simulations, demonstrating that the attitude estimation error converges to zero even when a very large initial estimation error is chosen.

VII Experimental results

In this section, we demonstrate the performance of the proposed observer (11) with real world experiments. The experiments have been done in collaboration with Maptek Pty Ltd444http://www.maptek.com/. The experimental setup consists of Maptek’s I-Site Drive system 555http://www.maptek.com/products/i-site/i-site_drive.html together with a Pointgrey Grasshopper3 GS3-U3-23S6C-C camera666https://www.ptgrey.com/grasshopper3-23-mp-color-usb3-vision-sony-pregius-imx174 mounted on a truck (see Fig 3). The I-Site Drive includes the NovAtel SPAN-CPT777http://www.novatel.com/products/span-gnss-inertial-systems/span-combined-systems/span-cpt/ high precision Inertial Navigation System (INS) consist of a RTK GPS and an inertial measurement unit, providing accurate estimates of the 6DoF transformation of the I-Site Drive system wrt. the inertial frame which is used as the ground truth $(R_{k},p_{k})$ in our experiments. The extrinsic camera to I-Site Drive calibration transformation have been carefully computed using a hand-eye calibration method that enables computing ground truth transformation of the camera wrt. the NED frame [36]. The camera is configured to capture $1920\times 1200$ pixels images at approximately $22$ fps. The camera is hardware synced to the I-Site Drive system to provide precise time stamping of measurements for accurate evaluation of the results. In our experiments, we use ORB-SLAM [13] as the visual odometry algorithm. ORB-SLAM is a state-of-the-art large scale simultaneous localisation and mapping software. Direct outputs of ORB-SLAM are estimates of the translation of the $k$-th camera frame wrt. the initial camera frame. Given two consecutive transformation estimates $(\check{R}^{{C}_{0}}_{{C}_{k}},\check{p}^{{C}_{0}}_{{C}_{k}})$ and $(\check{R}^{{C}_{0}}_{{C}_{k+1}},\check{p}^{{C}_{0}}_{{C}_{k+1}})$ as direct outputs of ORB-SLAM, we compute $R^{{C}_{k}}_{{C}_{k+1}}=\check{R}^{{C}_{0}}_{{C}_{k}}{\!}^{\top}\check{R}^{{C}_{0}}_{{C}_{k+1}}$ and $p^{{C}_{k}}_{{C}_{k+1}}=\check{R}^{{C}_{0}}_{{C}_{k}}{}^{\top}(\check{p}^{{C}_{0}}_{{C}_{k+1}}-\check{p}^{{C}_{0}}_{{C}_{k}})$ (recalling (1) and (29)) and we use them in the observer (11). For our experiments, we disable the loop closure thread of ORB-SLAM to observe the inherent drift of the pure visual odometry888We recall that effective loop closure is not possible in some mining applications due to the restrictions of vehicle’s path (not necessarily containing loops), high perceptual aliasing, large moving vehicles, and change in the appearance of the environment due to dust, wind, etc..

After careful calibration, the camera intrinsic parameters are estimated as $(1061.01,1062.22)$ for the focal length and $(951.28,604.86)$ for the principal point. The truck follows the path shown in Fig. 3 during which $4000$ images are captured. Since raw velocity measurement log was not available in our experiment, we use the numerical difference of the position estimates of the INS system as estimates of the linear velocity in our observer. This yields noisy estimates of the linear velocity, but is good enough for the purpose of demonstration in this section. We choose $L=0.004I_{3}$ and we initialise the observer (11) with $\hat{R}_{0}=I_{3}$ since the vision-GPS system does not initially have any information about the rotation of the camera wrt. the NED frame. Fig. 4 shows that despite the very large initial attitude estimation error, $\angle E_{k}$ converges to a very close neighborhood of zero. The observer is capable of accumulating the information of velocity measurements over time to asymptomatically estimate the camera’s orientation wrt. the NED frame. Fig. 5 shows the Euler angle decomposition of the ground truth orientation $R_{k}$ versus its estimate $\hat{R}_{k}$. The Euler angles estimates follow the ground truth angles very closely, demonstrating excellent performance of the observer. The covariance of the INS attitude estimate is about $0.2$ degrees in our experiment. Fig. 4 and Fig. 5 show that the vision-GPS system is capable of providing the level of accuracy close to a high-end INS system, albeit at a fraction of its cost.

As explained in Section II, one purpose of fusing GPS measurements with vision is to mitigates the inherent drift due to accumulation of the visual odometry error. In order to show this effect, we use the ground truth initial orientation $R^{A}_{{C}_{0}}$ to transform the direct orientation estimates of ORB-SLAM $\check{R}^{{C}_{0}}_{{C}_{k}}$ to the NED frame via $\check{R}_{k}:=R^{A}_{{C}_{0}}\check{R}^{{C}_{0}}_{{C}_{k}}$. This gives the estimates of $R_{k}$ with pure visual odometry if the initial attitude was known. We also initialise the observer (11) with $\hat{R}_{0}=R^{A}_{{C}_{0}}$ and obtain the attitude estimates $\hat{R}_{k}$. Fig. 6(a) shows that the attitude estimation error of the pure visual odometry $\angle(\check{R}_{k}R_{k}^{\top})$ (red curve) drifts over time while the proposed observer (11) successfully compensates for this drift by fusing the velocity such that the magnitude of the attitude error $\angle E_{k}$ (blue curve) remains bounded for all times. All of the results presented so far are obtained by processing the full resolution images. The visual odometry algorithm might not yield real time performance with full resolution of the camera. Here, we resized the images to $960\times 600$ and perform the same experiment as Fig. 6(a). Fig. 6(b) compares the attitude estimation error of the pure visual odometry versus that of the observer (11). Compared to Fig. 6(a), the accumulation of the visual odometry error is significantly higher while the observer error remains almost unchanged. Lastly, we point out that the presented results so far are obtained after very careful calibration of the camera intrinsic parameters. Even slight calibration error may significantly increase the drift of the pure visual odometry results while such drifts are compensated when the observer (11) is employed. To demonstrate this, we add $1$ percent error to the camera intrinsic parameters used in ORB-SLAM and reprocess the full resolution images. Fig. 6(c) compares the attitude estimation of the pure visual odometry versus the observer (11) both initialised with the ground truth orientation. Compared to Fig. 6(a), the attitude estimation error of the pure visual odometry significantly increases while the observer (11) successfully compensates for such a drift.

VIII Conclusion

This paper formulates the vision and GPS velocity fusion as an observer design problem on the Lie group $\text{SO}(3)$with a single vectorial output measurement whose associated reference output is time-varying. Inspired by [27, 22, 29, 5, 9], we propose a discrete-time attitude observer and provide a rigorous stability analysis showing the exponential convergence of the attitude estimation error to zero. We present experimental studies demonstrating that the observer is effectively capable of mitigating the inherent drift of pure visual odometry estimates. The observer also enables direct attitude estimation wrt. the NED frame, making it ideal for applications involving navigation and control of robots in the inertial frame. Where depth (i.e. scale) is available through stereo vision, extending the proposed observer to a estimate the pose on SE(3) is a potential future work.

Appendix A

GPS measurement model (2): similar to (1), one can concatenate the camera positions $p^{{C}_{0}}_{{C}_{k+1}}$ to obtain

[TABLE]

On the other hand, we have

[TABLE]

where ${R^{A}_{{C}_{0}}}^{\!\!\!\!\top}$ transforms the vector $p^{A}_{{C}_{k+1}}-p^{A}_{{C}_{k}}$ from the reference frame $\{A\}$ to the frame $\{C_{0}\}$. Substituting for (30) into (29) and noting that GPS positions are related to $p^{A}_{{C}_{k+1}}-p^{A}_{{C}_{k}}$ via $p^{A}_{k+1}-p^{A}_{k}=d(p^{A}_{{C}_{k+1}}-p^{A}_{{C}_{k}})$, yields (2).

Appendix B

Linearizion of error dynamics (13): using $E_{k}\approx I_{3}+{\epsilon_{k}}_{\times}$, we have \big{(}L(E_{k}\bar{p}^{A}_{k}-\bar{p}^{A}_{k})\!\big{)}_{\!\times}E_{k}\bar{p}^{A}_{k}\approx\big{(}L((I_{3}+{\epsilon_{k}}_{\times})\bar{p}^{A}_{k}-\bar{p}^{A}_{k})\!\big{)}_{\!\times}(I_{3}+{\epsilon_{k}}_{\times})\bar{p}^{A}_{k}\approx\big{(}L{\epsilon_{k}}_{\times}\bar{p}^{A}_{k}\big{)}_{\!\times}(I_{3}+{\epsilon_{k}}_{\times})\bar{p}^{A}_{k}\approx\big{(}L{\epsilon_{k}}_{\times}\bar{p}^{A}_{k}\big{)}_{\!\times}\bar{p}^{A}_{k} where we ignore the terms that are of order two or higher on $\epsilon_{k}$. Using the Taylor series expansion of the exponential map, we have \exp\!\Big{(}\Big{(}\!\!\big{(}L(E_{k}\bar{p}^{A}_{k}-\bar{p}^{A}_{k})\!\big{)}_{\!\times}E_{k}\bar{p}^{A}_{k}\Big{)}_{\!\times}\Big{)}E_{k}\approx(I_{3}+\big{(}(L{\epsilon_{k}}_{\times}\bar{p}^{A}_{k})_{\!\times}\bar{p}^{A}_{k}\big{)}_{\times}+\text{HOT})(I_{3}+{\epsilon_{k}}_{\times})\approx I_{3}+{\epsilon_{k}}_{\times}+\big{(}(L{\epsilon_{k}}_{\times}\bar{p}^{A}_{k})_{\!\times}\bar{p}^{A}_{k}\big{)}_{\times}+\text{HOT}. Substituting this into (13) and considering $E_{k+1}\approx I_{3}+{\epsilon_{k+1}}_{\times}$, we have I_{3}+{\epsilon_{k+1}}_{\times}=I_{3}+{\epsilon_{k}}_{\times}+\big{(}(L{\epsilon_{k}}_{\times}\bar{p}^{A}_{k})_{\!\times}\bar{p}^{A}_{k}\big{)}_{\times}. Canceling $I_{3}$ from the sides and taking the inverse of the operator $(.)_{\times}$ from the sides, we obtain $\epsilon_{k+1}=\epsilon_{k}+(L{\epsilon_{k}}_{\times}\bar{p}^{A}_{k})_{\!\times}\bar{p}^{A}_{k}$. Observing $(L{\epsilon_{k}}_{\times}\bar{p}^{A}_{k})_{\!\times}\bar{p}^{A}_{k}=\bar{p}^{A}_{k}{}_{\times}L\bar{p}^{A}_{k}{}_{\times}\epsilon_{k}$, noting $L=lI_{3}$, and using the property $a_{\times}b_{\times}=ba^{\top}-a^{\top}bI_{3}$, yields (16).

Appendix C

Proof of (23): using (16), we have

[TABLE]

Since the maximum singular value of $P_{i}$ is one, and using (31) we have \sum_{i=k}^{k+T}(\epsilon_{i}^{\top}-\epsilon_{k}^{\top})P_{i}(\epsilon_{i}-\epsilon_{k})\leq\sum_{i=k}^{k+T}(\epsilon_{i}^{\top}-\epsilon_{k}^{\top})(\epsilon_{i}-\epsilon_{k})=l^{2}\sum_{i=k}^{k+T}\big{(}\sum_{j=k}^{i-1}\epsilon_{j}^{\top}P_{j}\big{)}\big{(}\sum_{j=k}^{i-1}P_{j}\epsilon_{j}\big{)}. Using the Cauchy–-Schwarz inequality [37] and noting $P_{j}P_{j}=P_{j}$, we obtain $\sum_{i=k}^{k+T}(\epsilon_{i}^{\top}-\epsilon_{k}^{\top})P_{i}(\epsilon_{i}-\epsilon_{k})\leq l^{2}\sum_{i=k}^{k+T}(i-k)\sum_{j=k}^{i-1}\epsilon_{j}^{\top}P_{j}P_{j}\epsilon_{j}=l^{2}\sum_{i=k+1}^{k+T}(i-k)\sum_{j=k}^{i-1}\epsilon_{j}^{\top}P_{j}\epsilon_{j},$ where we used the fact that the value of the first summand evaluated at $i=k$ is zero and hence we can start the first summation from $i=k+1$. Changing the order of summations we obtain

[TABLE]

We have $\sum_{i=j+1}^{k+T}(i-k)=\frac{T(T+1)}{2}-\frac{(j-k+1)(j-k)}{2}$. Using the limits of the first summation we have $k\leq j\leq k+T-1$ which implies $0\leq(j-k+1)(j-k)\leq T(T-1)$. Hence $\sum_{i=j+1}^{k+T}(i-k)\leq\frac{T(T+1)}{2}$. Using this and adding the positive term $\frac{l^{2}T(T+1)}{2}\epsilon_{k+T}^{\top}P_{k+T}\epsilon_{k+T}$ to the right hand side of (32) yields (23) and completes the proof.