Incrementally Learned Mixture Models for GNSS Localization

TL;DR

This paper presents an online variational Bayesian sensor fusion method that learns and adapts Gaussian mixture models for improved GNSS localization accuracy in the presence of non-Gaussian errors.

Contribution

It introduces an incremental variational mixture algorithm that automatically adapts the number of components to model complex error distributions in GNSS data.

Findings

Outperforms state-of-the-art methods on real-world datasets

Automatically adapts the number of mixture components

Achieves superior localization accuracy

Abstract

GNSS localization is an important part of today's autonomous systems, although it suffers from non-Gaussian errors caused by non-line-of-sight effects. Recent methods are able to mitigate these effects by including the corresponding distributions in the sensor fusion algorithm. However, these approaches require prior knowledge about the sensor's distribution, which is often not available. We introduce a novel sensor fusion algorithm based on variational Bayesian inference, that is able to approximate the true distribution with a Gaussian mixture model and to learn its parametrization online. The proposed Incremental Variational Mixture algorithm automatically adapts the number of mixture components to the complexity of the measurement's error distribution. We compare the proposed algorithm against current state-of-the-art approaches using a collection of open access real world datasets and demonstrate its superior localization accuracy.

Figures

Tables

Table 1. TABLE I : Robust Factor Graph Algorithms categorized regarding the applied Error Model

Error Model	Static	Dynamic
Gaussian	Standard Factor Graph	SC [11] DCE [12]
Non-Gaussian	DCS [13] cDCE [12] Max-Mix. [4] Sum-Mix. [5]	Self-tuning M-Est. [14] Self-tuning Mix. [6] Adaptive Mix. [7] Batch Cov. Est. [8]

Table 2. TABLE II : Legend of Fig. 3

Symbol	Name
$𝓘$	Information Matrix of the GMM Components
$𝝁$	Mean of the GMM Components
$w$	Weight of the GMM Components
$𝐞$	Measurement Error
$𝒔$	Hidden Variable
$\nu$	Wishart Degree of Freedom
$𝑽$	Wishart Scale Matrix
$\beta$	Information Matrix Scaling
$N$	Number of Measurements

Table 3. TABLE III : Accuracy and Runtime of the Final Evaluation.

	Chemnitz		Berlin PP		Berlin GM		Frankfurt MT		Frankfurt WT
Algorithm	ATE [m]	Time [s]	ATE [m]	Time [s]	ATE [m]	Time [s]	ATE [m]	Time [s]	ATE [m]	Time [s]
Gaussian	30.0	58.6	29.2	9.5	13.38	49.7	30.97	58.1	23.54	31.9
DCS [13]	4.403	54.9	25.04	14.7	19.11	64.5	13.25	69.5	11.39	35.2
cDCE [12]	4.326	54.5	17.91	14.3	14.59	65.3	14.93	73.6	11.12	34.9
SM+EM [7]	2.378	102.0	12.45	39.8	13.88	163.0	10.72	136.0	6.42	76.5
SM+VBI	3.106	119.0	12.4	33.3	13.55	186.0	11.23	131.0	4.117	84.4
SM+EM+CL	6.723	172.0	12.23	45.0	13.37	312.0	10.41	214.0	9.706	152.0
IVM (SM+VBI+CL)	2.48	122.0	11.56	59.4	10.99	352.0	8.586	248.0	4.626	146.0

Table 4. TABLE IV : Parameters of the Incremental Variational Mixture Algorithm

Parameter Name	Symbol	Value
Wishart Prior Degree of Freedom	${\nu }_0$	2
Wishart Prior Scale Matrix	${𝐕}_0$	$\mathrm{var}\left(\widehat{𝐞}\right)/\nu$
Normal Prior Information Matrix	${\beta }_0$	${10}^{-6}{m}^{-2}$
Pruning Threshold	$w_{min}$	$1/N$
Max. Iterations	$i_{max}$	$1000$
Min. $\Delta$Likelihood	$\Delta L_{min}$	${10}^{-6}$
Max. Number of Components	$K_{max}$	$8$

Equations

$${\mathbf{X}}^{\ast }=\underset{\mathbf{X}}{\mathrm{argmax}}\mathbf{P}(\mathbf{X}|\mathbf{Z})$$

$$\hat{\mathbf{X}}=\underset{\mathbf{X}}{\mathrm{argmax}}\mathbf{P}(\mathbf{Z}|\mathbf{X})\mathbf{P}(\mathbf{X})$$

$$\begin{array}{rl}\mathbf{P}(\mathbf{Z}|\mathbf{X})\propto & \prod _t\left[\prod _i\mathbf{P}({\mathbf{z}}_{t,i}^{pr}|{\mathbf{x}}_t^{x,y,z,\delta })\right]\\ \cdot & \prod _t\mathbf{P}({\mathbf{z}}_t^{odo}|{\mathbf{x}}_{t+1}^{x,y,z,\varphi },{\mathbf{x}}_t^{x,y,z,\varphi })\\ \cdot & \prod _t\mathbf{P}({\mathbf{x}}_{t+1}^{\delta ,\dot{\delta }}|{\mathbf{x}}_t^{\delta ,\dot{\delta }})\end{array}$$

$$\hat{\mathbf{X}}=\underset{\mathbf{X}}{\mathrm{argmin}}\sum _n-\ln (\mathbf{P}({\mathbf{z}}_n|{\mathbf{x}}_n))$$

$$\hat{\mathbf{X}}=\underset{\mathbf{X}}{\mathrm{argmin}}\sum _n\frac{1}{2}{\Vert {\mathcal{I}}_n^{\frac{1}{2}}\left({\mathbf{e}}_n-{\mathbf{\mu }}_n\right)\Vert }^2$$

$$\hat{\mathbf{X}}=\underset{\mathbf{X}}{\mathrm{argmin}}\sum _n{\Vert \sqrt{-\ln \left(\frac{\mathbf{P}({\mathbf{z}}_n|{\mathbf{x}}_n)}{\gamma }\right)}\Vert }^2$$

$${\mathbf{e}}_{t,i}^{pr}\sim \sum _{k=1}^K{w}_k\cdot \mathcal{N}({\mathbf{\mu }}_k,{\mathcal{I}}_k)\text{ with }\sum _{k=1}^K{w}_k=1$$

$$\begin{array}{rl}\mathbf{P}({\mathbf{z}}_n|{\mathbf{x}}_n)=& \sum _k{c}_k\cdot \exp \left(-\frac{1}{2}{\Vert {\mathcal{I}}_k^{\frac{1}{2}}\left(\mathbf{e}-{\mathbf{\mu }}_k\right)\Vert }^2\right)\\ \text{with }{c}_k=& w_k\cdot \det \left({\mathcal{I}}_k^{\frac{1}{2}}\right)\end{array}$$

$${\mathbf{\theta }}^{\ast }=\underset{\mathbf{\theta }}{\mathrm{argmax}}\mathbf{P}(\mathbf{\theta }|\mathbf{e})$$

$$\begin{array}{rl}& {\mathbf{\theta }}^{\ast }=\underset{\mathbf{\theta }}{\mathrm{argmax}}\mathbf{P}(\mathbf{e}|\mathbf{\theta })\\ \text{with }& \mathbf{P}(\mathbf{e}|\mathbf{\theta })=\int \mathbf{P}(\mathbf{e},\mathbf{s}|\mathbf{\theta })\mathrm{d}\mathbf{s}\end{array}$$

$${\mathbf{\theta }}^{\ast }=\underset{\mathbf{\theta }}{\mathrm{argmax}}\frac{\mathbf{P}(\mathbf{e}|\mathbf{\theta })\mathbf{P}(\mathbf{\theta })}{\int \mathbf{P}(\mathbf{e}|\mathbf{\theta })\mathbf{P}(\mathbf{\theta })\mathrm{d}\mathbf{\theta }}$$

$$\mathbf{\mu }\sim \mathcal{N}(0,{\beta }_{0}I)$$

$$\mathcal{I}\sim \mathcal{W}({\mathbf{V}}_0,{\nu }_0)$$

$$\mathbf{P}(\mathbf{\theta }|\mathbf{e})\approx Q(\mathbf{\theta })$$

$$Q(\mathbf{\theta })=\prod _j{Q}_j({\mathbf{\theta }}_j)=Q_{\mathcal{I}}(\mathcal{I})Q_{\mathbf{\mu }}(\mathbf{\mu })Q_{\mathbf{s}}(\mathbf{s})$$

$$\begin{array}{rl}& \hat{\mathbf{X}}=\underset{\mathbf{X}}{\mathrm{argmax}}\mathbf{P}(\mathbf{X}|\mathbf{Z})\\ \text{with }& \mathbf{P}(\mathbf{X}|\mathbf{Z})=\int \mathbf{P}(\mathbf{X},\mathbf{\theta }|\mathbf{Z})\mathrm{d}\mathbf{\theta }\end{array}$$

$$\hat{\mathbf{\theta }}=\underset{\mathbf{\theta }}{\mathrm{argmax}}\mathbf{P}(\mathbf{\theta }|{\hat{\mathbf{e}}}_{t-1})$$

$${\hat{\mathbf{X}}}_t=\underset{\mathbf{X}}{\mathrm{argmax}}\mathbf{P}(\mathbf{X}|{\mathbf{Z}}_t,E({\hat{\mathbf{\theta }}}_t))$$

$${\hat{\mathbf{e}}}_t=\mathbf{e}({\hat{\mathbf{X}}}_t,{\mathbf{Z}}_t)$$

$$E(\mathcal{I})={\nu }_0\cdot {\mathbf{V}}_0^{-1}$$

$${\mathbf{V}}_0=\frac{\mathrm{var}(\hat{\mathbf{e}})}{{\nu }_0}$$

Full text

We introduce a novel sensor fusion algorithm based on variational Bayesian inference, that is able to approximate the true distribution with a Gaussian mixture model and to learn its parametrization online. The proposed Incremental Variational Mixture algorithm automatically adapts the number of mixture components to the complexity of the measurement’s error distribution. We compare the proposed algorithm against current state-of-the-art approaches using a collection of open access real world datasets and demonstrate its superior localization accuracy.

I Introduction

A reliable localization is crucial for autonomous systems like self-driving cars or mobile robots. While global navigation satellite systems (GNSS) can provide this information for open-field scenarios, their accuracy and consistency is decreased in urban areas. The reason are signal reflections on tall buildings, which cause non-line-of-sight (NLOS) measurements. Recent publications like [1] demonstrated the non-Gaussian characteristic of the resulting error distributions. Therefore, a lot of effort has already been invested to make GNSS more robust against NLOS effects [2, 3].

One approach to reduce the influence of these errors, is to incorporate knowledge about their true distribution into the state estimation process. With robust factor graph optimization, there already is a class of algorithms that allows to include such distributions either in a predefined [4, 5] or self-tuning [6, 7, 8] way. To describe the non-Gaussian errors, Gaussian mixture models (GMM) are the established state-of-the-art. However, existing approaches are not able to fully adapt the GMM during the online estimation process. Therefore, they do not address applications where the error characteristic changes significantly over time. For GNSS this can easily happen, since the NLOS error depends on the receiver’s surrounding.

Based on our recently published Adaptive Mixture algorithm [7], we want to introduce a novel approach that combines variational Bayesian inference [9] with non-linear graph optimization to achieve an even better robustness and accuracy. Due to an incremental construction of the GMM, we preserve the real-time capability of the original approach and eliminate existing drawbacks.

Since the proposed algorithm is the result of an incremental research process, we want to give an overview over prior work in Sec. II. How GNSS localization and Gaussian mixtures are connected to graph based optimization is explained in Sec. III, while Sec. IV gives an overview about the algorithms that are used to estimate Gaussian mixtures from data. Then we combine these techniques to the proposed Incremental Variational Mixture approach in Sec. V before we evaluate it in Sec. VI using several real world GNSS datasets and finally conclude.

II Prior Work

Several recently published algorithms try to solve the problem of non-Gaussian sensor fusion in different ways. From our point of view, they can be separated by two properties of their error models: If they are using static or dynamic ones and whether the models are Gaussian or not. We summarize existing work in Tab. I, categorized by the error model’s properties. For example, the standard least squares approach [10] is a static Gaussian method since it uses a Gaussian error model that is predefined.

With Switchable Constraints (SC) [11] and it’s modification Dynamic Covariance Estimation (DCE) [12] exist two algorithms that implement Gaussian but dynamic error models. While SC introduces tunable weights for each Gaussian, DCE estimates the covariance of the distribution. Both are able to handle Gaussian distributions with outliers but tend to fail for heavily multimodal distributions that can occur when GNSS is used in urban canyons.

A common robust method, not only in optimization, is the application of predefined non-Gaussian weight functions, so called M-estimators. These static non-Gaussian methods include not only the well-known Huber or Tukey estimators but also techniques like Dynamic Covariance Scaling (DCS) [13] or the closed form of DCE (cDCE) [12]. Usually they come with a tuning parameter to adjust the trade-off between robustness and well-behaved convergence. Beside M-estimators, there are two approaches to include different distributions in the optimization problem. Max-Mixtures [4] allows approximated inference over a Gaussian mixture by replacing the sum-operator by a max-operator. The work of Rosen [5] includes all kinds of continuous probability distributions including Gaussian mixtures. However, all methods with static error models require exact knowledge about the expected error distribution. For GNSS in changing environments, this information is often not available.

Algorithms that apply non-Gaussian error models and estimate their parameter dynamically are a relatively new class of robust solutions to non-linear optimization. The authors of [14] were the first who described an algorithm that is capable to estimate the optimal parameter of a certain group of M-estimators, using adaptive importance sampling. Later in [6], the parameters of a Gaussian mixture model (GMM) were optimized simultaneously during the estimation process. In recent work [7], Expectation-Maximization (EM) was applied to estimate the GMM in a more reliable way. Also, Wang et al. introduced a similar algorithm in the field of laser scan matching [15]. They estimate a mixture of exponential power to describe the residuals of the point-to-point registration process. The prerequisite for this type of algorithms is the selection of a suitable model respectively its complexity. For GNSS applications this can be a limitation, since the complexity (e.g. the number of modes) can change with the environment over time.

Watson et al. currently proposed a novel algorithm, based on the infinite Gaussian mixture model, which is adaptive regarding its number of components [8]. Since the applied Gibbs sampling is not real-time capable, the approach can only be applied to batch problems. For the GNSS localization of an autonomous system this is a harsh limitation, since it usually has to be solved in real time. Also, the algorithm applies one common model for the whole dataset, so it cannot address time dependent changes.

While existing algorithms tune a predefined distribution to approximate the estimation error, we want to go one step further and learn its structure and parametrization online. Therefore, we want to propose an algorithm that is able to exploit the changing complexity of the error distribution in urban scenarios. Compared to the most recent self-tuning mixtures algorithm [7], our approach is novel in the following ways:

The proposed algorithm incrementally constructs a Gaussian mixture with a (theoretically) unlimited number of components. 2. 2.

We are able to remove unused Gaussian components to preserve real-time capability. 3. 3.

To ensure a numerically stable solution, we apply variational Bayesian inference instead of Expectation-Maximization to estimate the GMM.

With this combination of techniques, we propose the first real-time capable sensor fusion algorithm, that is able to learn the complexity and parametrization of the measurement’s error distribution over time.

III GNSS Localization as Factor Graph

In the following section, we want to give a brief overview how the GNSS problem can be described as factor graph and how non-Gaussian error distributions can be represented. Factor graphs, a graphical representation of least squares problems, are wildly used in robotics. Therefore, we assume basic knowledge about their theoretical properties and refer readers to [10] for a more detailed introduction.

III-A Factor graphs and Least Squares

The GNSS localization of a moving system can be described as optimization problem (missing) 1, where $\mathbf{Z}$ is a set of pseudorange and odometry measurements and $\mathbf{X}$ is the set of estimated states including the vehicle’s position. $\mathbf{X}^{\ast}$ is the most likely set of states according to the set of measurements. By applying the Bayes theorem, it can be written as maximum a posteriori (MAP) problem (missing) 2 with $\mathbf{\hat{X}}$ as MAP estimator of the true state variables.

[TABLE]

[TABLE]

At each time step $t$, one odometry measurement $\mathbf{z}^{odo}_{t}$ and multiple pseudorange measurements $\mathbf{z}^{pr}_{t,i}$, each corresponding to the $i\text{th}$ satellite, are available. The set of state variables can be separated in two subsets that consist of the 3D pose $\mathbf{x}^{x,y,z,\phi}_{t}$ and the GNSS specific clock error state $\mathbf{x}^{\delta,\dot{\delta}}_{t}$. The Cartesian coordinates $x,y,z$ are in the earth-centered-earth-fixed (ECEF) frame, while $\phi$ denotes the rotation around the vehicles upright axis. The clock error $\delta$ and its derivation $\dot{\delta}$ are required for the constant clock error drift model (CCED). By assuming uninformative priors, $\mathbf{P}(\mathbf{X})$ can be omitted and the posterior likelihood is expanded to (missing) 3. A detailed explanation of these error terms can be found in [6]. Fig. 2 shows a small example of the corresponding factor graph.

[TABLE]

Instead of maximizing the likelihood (missing) 2 directly, $\mathbf{\hat{X}}$ is usually estimated by minimizing the negative log likelihood:

[TABLE]

One common index $n$ is used for all factors here. By defining $\mathbf{P}(\mathbf{Z}|\mathbf{X})$ as normal distribution $\mathcal{N}({\boldsymbol{\mu}}_{n},\boldsymbol{\mathcal{I}}_{n})$ over the non-linear error function $\mathbf{e}_{n}=f(\mathbf{x}_{n},\mathbf{z}_{n})$ with mean ${\boldsymbol{\mu}}_{n}$ and information matrix $\boldsymbol{\mathcal{I}}_{n}$, the problem can be formulated as non-linear least squares estimation:

[TABLE]

As already mentioned, the Gaussian assumption does not apply to all GNSS measurements. Therefore, a reliable state estimation requires more advanced distributions.

III-B Factor Graphs and Gaussian Mixtures

Since empirical sensor distribution can be asymmetric or even multimodal, a class of distributions is required that is able to represent these properties. For this work we choose GMMs, because they fulfill this condition and have several positive properties regarding their application in autonomous systems [16]. To apply them to a least squares problem, we use a solution proposed in [5] that allows almost arbitrary non-Gaussian distributions $\mathbf{P}(\mathbf{z}_{n}|\mathbf{x}_{n})$, using the following:

[TABLE]

The normalization constant $\gamma$ ensures a numerically stable solution. We reference the approach as Sum-Mixture (SM) in difference to Max-Mixture of [4], which only approximates a GMM. In this work, a $K\text{-component}$ Gaussian mixture defined by (missing) 7 is used to describe the pseudorange error $\mathbf{e}^{pr}_{t,i}=f(\mathbf{x}^{x,y,z,\delta}_{t},\mathbf{z}^{pr}_{t,i})$. Each Gaussian component $k$ is scaled with a normalized weight $w_{k}$.

[TABLE]

Since our approach is not limited to the pseudorange error, we use the generic error $\mathbf{e}$ instead of $\mathbf{e}^{pr}_{t,i}$ for the following equations.

The likelihood of the GMM is defined by:

[TABLE]

Therefor, the normalization constant can be set to $\gamma=\sum_{k}c_{k}$. For a more detailed derivation of these equations we want tor refer the reader to the original work [5] and to [7] for its application to GMMs.

IV Mixture Model Estimation

As shown in prior work [14, 8, 7], the residual $\mathbf{e}(\mathbf{x},\mathbf{z})$ of the estimation problem can be used to get an approximation of the measurement’s true distribution. To describe the measurement error with a $K$-component Gaussian mixture model, the GMM’s parameters $\boldsymbol{\theta}=\left\{w,{\boldsymbol{\mu}},\boldsymbol{\mathcal{I}}\right\}$ have to be estimated with (missing) 9 from $N$ samples of $\mathbf{e}$.

[TABLE]

The set of parameters $\boldsymbol{\theta}$ includes a weight $w_{k}$, mean ${\boldsymbol{\mu}}_{k}$ and information matrix $\boldsymbol{\mathcal{I}}_{k}$ for each component $k$ of $K$. Since the probability that measurement $n$ belongs to component $k$ is not known, the corresponding hidden parameter $\boldsymbol{s}=\left\{s_{kn}\right\}$ have to be included in the estimation problem. The directed graphs of two possible interpretations are shown in Fig. 3 and explained in the following subsections.

IV-A Expectation Maximization

One interpretation, shown in 3(a), is to treat $\boldsymbol{\theta}$ as parameters and describe the estimation as maximum-likelihood problem:

[TABLE]

Since the integral is intractable, the iterative Expectation-Maximization (EM) algorithm [17] splits the problem into two subproblems: The E-step estimates the hidden variable $\boldsymbol{s}$ based on an initial guess of $\boldsymbol{\theta}$. The M-step estimates the maximum likely set of parameters $\boldsymbol{\theta}$ based on the hidden variables that were estimated before. As demonstrated in [7], the EM can be applied in an alternating sequence with the least squares optimization of the original state estimation problem (missing) 1. On the downside, the choice of the parameter $K$ is not trivial because it depends on the complexity of the true distribution of $\mathbf{e}$. While a too small $K$ results in less robustness, a too large number of components leads to numerical instabilities. These instabilities can occur, if only a few samples are responsible for a specific component. There exists a variety of EM derivations that aim to overcome this limitation like the greedy EM [18] or the split and merge EM [19], but their parametrization can be even more difficult. Even a simple merging criterion is not easy to chose, as the author of [20] demonstrates.

IV-B Bayesian Methods

A fully Bayesian interpretation of the GMM optimization (missing) 9, treats the distribution parameters $\boldsymbol{\theta}$ as random variables with corresponding distributions. After adding a set of priors $\mathbf{P}(\boldsymbol{\theta})$, the estimation problem can be written as:

[TABLE]

There exist different suggestions for the choice of this priors. In difference to the more common model used in [21, 16, 8] or [22, p.475], we follow the suggestion of the authors of [9] and omit the Dirichlet prior of the mixture’s weight $w$. 3(b) shows the corresponding graph with the remaining Wishart prior (missing) 13 for the information matrix and the normal prior (missing) 12 for the mean of the GMM. The matrix $\operatorname{\boldsymbol{I}}$ is the identity matrix here.

[TABLE]

[TABLE]

To omit the prior of $w$ comes with two advantages, one obvious is the reduced number of hyperparameters which makes the algorithm easier to apply. The other one is the possibility to reduce weights to zero, which allows us to remove unused components from the model. Since the integral in (missing) 11 can not be calculated directly, approximate solutions like Markov Chain Monte Carlo (MCMC) methods or variational Bayesian inference (VBI) have to be applied. MCMC algorithms, such as Gibbs sampling [23] used in [16] and [8], approximate the integral over the set of parameters $\boldsymbol{\theta}$ by generating many samples of them from their estimated distribution. This process is computationally expensive and therefor rather unsuitable for real-time applications, which leads us to variational inference.

The key idea of VBI is to approximate the posterior $\mathbf{P}(\boldsymbol{\theta}|\mathbf{e})$ with a distribution $\operatorname{\mathbf{Q}}(\boldsymbol{\theta})$.

[TABLE]

To find the distributions over all hidden variables $\boldsymbol{\theta}$, it is assumed that they can be partitioned in $j$ independent groups $\left\{\boldsymbol{\theta}_{j}\right\}$. According to these groups, $\operatorname{\mathbf{Q}}(\boldsymbol{\theta})$ can be factorized, as shown in (missing) 15, and each factor $\operatorname{\mathbf{Q}}_{j}(\boldsymbol{\theta}_{j})$ can be estimated separately. For the model proposed in [9], the groups are mean ${\boldsymbol{\mu}}$, information matrix $\boldsymbol{\mathcal{I}}$ and the correspondence variable $\boldsymbol{s}$.

[TABLE]

Since we use exactly the same model that is proposed in [9], we refer interested readers to the original work. For a more gentle introduction we can also recommend the tutorial paper [24].

Nevertheless, there are still some drawbacks of Bayesian inference in general and of the concrete VBI solution. Bayesian estimators require a set of predefined priors to represent the knowledge that is available in advance. If there is little to no prior knowledge about the estimated parameters, it can be hard to specify meaningful priors. This work can not provide a general solution for this problem, but we try to find a good trade-off between uninformative priors and prior knowledge that exists about NLOS errors. Variational inference require a predefined number of Gaussian components similar to the maximum-likelihood approach. However, our approach is robust against an exaggerated number of components, since the weight of unused ones drops to zero without causing singularities.

As a last statement of this section, we want to emphasize that to our knowledge no “best” algorithm exists to estimate a Gaussian mixture from empirical data. Nevertheless, VBI seems to be well suited for the problem we want to address.

V Learning Complexity Online

In this section, we describe the complete algorithm which is composed of robust factor graph optimization and variational Bayesian inference. The factor graph applies the GMM to robustly estimate the system’s state and the VBI estimates the mixtures based on the empirical distribution of the sensor. An overview of the general approach is given by Fig. 1 on the first page and an algorithmic description is provided with Alg. 1 and Alg. 2.

As explained in [7] the estimation problem (missing) 16 can be described as EM algorithm with the set of measurements $\mathbf{Z}$ as observed variable and the parameters of the mixture model $\boldsymbol{\theta}$ as hidden ones.

[TABLE]

Again, the intractable integral is solved by alternately estimating the expected value of the hidden parameters $\boldsymbol{\theta}$ (E-step) and the state variables (M-step). Since the problem is formulated as sliding window estimation, one iteration per time step is enough to achieve good convergence.

In the E-step, the posterior probability of $\boldsymbol{\theta}$ is calculated with VBI based on a previously determined MAP estimate of the error $\mathbf{\hat{e}}$:

[TABLE]

The desired expected value $\operatorname{\mathbb{E}}(\boldsymbol{\hat{\theta}}_{t})$ is computed during the variational inference. Again, further details are available in [9].

The M-step uses the estimated distribution parameter $\operatorname{\mathbb{E}}(\boldsymbol{\hat{\theta}}_{t})$ and applies it to the least squares optimization according to the Sum-Mixture formulation (missing) 6. After the optimization, the non-linear error function $\mathbf{e}$ can be used to get a MAP estimate of the error $\mathbf{e}$ as shown in:

[TABLE]

[TABLE]

The complete procedure is described with Alg. 2.

To adapt the mixture’s number of components $K$ to the complexity of the empirical distribution, we add one new component in each time step, starting with $K=2$ at $t=0$. This leads to a fast rising number of components. However, through the variational inference, the weight of unused components drops to zero. Therefore, components with a weight below a predefined threshold $w_{min}$ can be removed. We set this pruning threshold to $\nicefrac{{1}}{{N}}$ with $N$ as the number of measurements. Each added component is initialized according to its prior distribution with zero mean $\mu_{0}=0$ and the information matrix $\boldsymbol{\mathcal{I}}_{0}=\operatorname{\nu_{0}}\operatorname{\boldsymbol{V}_{0}^{-1}}$. With this simple approach, we are able to choose the right number of components just by “saturating” the mixture with components. Since the saturation point depends on the parametrization of the Gauss and Wishart prior, a careful parametrization is essential for a good performance.

To ensure real time capability even for critical applications, we specify an optional upper limit $K_{max}$ for the number of components. If this limit is exceeded, the component with the lowest weight is removed before a new one is added. The online complexity learning is described by Alg. 1.

VI Evaluation

In this section, we demonstrate the performance of the proposed Incremental Variational Mixture (IVM) algorithm on several real world GNSS datasets. The resulting trajectories are shown in Fig. 4. We do not only compare against the robust state-of-the-art approaches DCS [13] and cDEC [12] but also show the impact of VBI and the proposed Complexity Learning (CL) mechanism. Beside a comparison based on position error metrics, we want to demonstrate the sensitivity of our approach to its parametrization, which is important for practical applications. To proof the performance under real time conditions, we perform the estimation under online conditions. This means, results are calculated and stored without information from future measurements.

VI-A The Datasets

We use five different datasets, that consist of raw pseudorange measurements from a mass market receiver and wheel odometry from measurement vehicle. These datasets are collected in different urban areas in Germany and have an accumulated length of $18.6\text{\,}\mathrm{km}$ respectively $92\text{\,}\mathrm{min}$. A precise ground truth is provided by a combination of differential GNSS and a tactical grade inertial measurement unit. Technical details are published in [25] and the datasets themselves are available online 111http://mytuc.org/GNSS.

VI-B Parametrization

Since the choice of parameters has a significant impact on the performance of the proposed algorithms, we want to explain them clearly. The noise parameters of the factor graph are identical to our previous work, therefore we want to refer to [7] for detailed information.

To apply the variational inference, several parameters have to be specified that are summarized in Tab. IV. Not all of them are equally important, especially the convergence criteria and the pruning threshold $w_{min}$ have just a small impact on the overall performance. The optional maximum of Gaussian components $K_{max}$ is also easy to parametrize, it has to be high enough to reflect the worst case complexity of the empirical distribution. A value of $K_{max}=8\dots 12$ should be enough for all practical applications and if the algorithm is applied to a fully unknown problem, the limit can be omitted. The information matrix of the normal prior $\beta_{0}$ defines the initial uncertainty about each component’s mean. Since NLOS errors over $1\text{\,}\mathrm{km}$ are very unlikely, a value of \beta_{0}=\nicefrac{{1}}{{(1\text{,}\mathrm{km})^{2}}} seems to be sufficient.

More difficult is the parametrization of the Wishart prior $\mathcal{W}(\operatorname{\boldsymbol{V}_{0}},\operatorname{\nu_{0}})$. The scale matrix $\operatorname{\boldsymbol{V}_{0}}$ in combination with the degree of freedom $\operatorname{\nu_{0}}$ determines the expected information matrix of the estimated mixture components:

[TABLE]

A coarse guess of the information matrix can be estimated directly from the error $\mathbf{\hat{e}}$ by calculating its variance. So, a meaningful prior can be defined with:

[TABLE]

For the positive integer $\operatorname{\nu_{0}}$, a straightforward approximation is not possible. The authors of the variational model [9] recommend an informative Wishart prior, which means value of $\operatorname{\nu_{0}}=1$ or at least close to 1. We evaluated the proposed algorithm with different values and visualized the localization error in Fig. 5. Therefore, we chose a prior with $\operatorname{\nu_{0}}=2$ to achieve the best performance over all datasets. The reason for the strong impact of $\operatorname{\nu_{0}}$ to the overall result is the convergence of the VBI algorithm. An increased $\operatorname{\nu_{0}}$ improves the convergence speed and increases the number of removed mixture components. With $\operatorname{\nu_{0}}=1$, the number of components over-fits the empirical distribution which leads to more local minima in the state estimation problem and therefor reduced accuracy. Although, the performance of different values of $\operatorname{\nu_{0}}$ is not fully consistent over all datasets, the choice of $\operatorname{\nu_{0}}=2$ generalizes at least for the tested ones. The comparison to the EM based algorithm in Fig. 5 shows that its fixed number of components $K$ also have to be chosen carefully. We set $K=3$ for a fair comparison since it is the best choice over all datasets.

VI-C Implementation Details

We implemented the algorithm as part of our robust sensor fusion library libRSF222http://mytuc.org/libRSF. The VBI algorithm will be published as open source in a future release. The non-linear least squares optimization is based on the Ceres solver [26] and the GMM estimation is implemented in C++ using Eigen [27] without multi-threading. All tests were performed on a Intel i7-7700 system.

VI-D Results

Metric of our comparison is the absolute trajectory error (ATE). We define the ATE as euclidean distance between estimated trajectory and ground truth in the local XY-plane. Tab. III summarizes the mean ATE as well as the runtime of the evaluated algorithms. Along with the proposed IVM, we compare different combinations of self-tuning algorithms. Common is the Sum-Mixture (SM) method to represent a GMM, but we differentiate between EM and VBI estimation and whether the proposed Complexity Learning (CL) is applied or not. Please notice that SM+EM equals the Adaptive Mixture approach from [7] and SM+VBI+CL is the proposed IVM algorithm.

Both VBI and CL seem to have a positive influence regarding the localization accuracy. Although, only in combination, they are able to consistently reduce the ATE over all datasets. Therefore, we would prefer the proposed IVM for applications where a maximum precision is required. Drawback of the CL is increased runtime of the algorithm, nevertheless it is 5 times faster than the recording time of the datasets.

If runtime is more important, we would tend to prefer the Adaptive Mixture (SM+EM) approach, since it offers a solid performance without the trouble of parameterizing the variational priors.

The results of the M-estimators DCS and cDCE show improvements compared to the non-robust factor graph, but both fall behind the GMM based methods.

VII Conclusion

The proposed algorithmic approach, allows to apply least squares optimization to sensor fusion problems with non-Gaussian and time dependent error distributions like GNSS localization in urban environments. We not only adapt the parameters of the Gaussian mixture representation, we also learn the right number of components to represent its complexity. The applied variational Bayesian inference allows a numerically stable and real time capable solution.

The comparison against several state-of-the-art algorithms were performed on a set of open access GNSS datasets that reflect the localization of an autonomous vehicle in an urban environment. Our approach demonstrated a superior estimation quality over the majority of datasets in combination with an increased but acceptable computation time.

The described connection between state estimation and adaptive error models opens a wide field of possible future improvements. Since the parametrization of the variational priors is the most critical drawback currently, we will address it in future work. Variational Bayesian inference also offers further methods to adaptively split and merge Gaussian components which could be explored in this context. Finally, a broader evaluation on different domains like SLAM would be useful to validate the approach and the choice of parameters.