Uncertainty Model Estimation in an Augmented Data Space for Robust State Estimation

TL;DR

This paper enhances robust state estimation for autonomous systems by augmenting the data space with metadata, improving the accuracy of measurement uncertainty modeling and resulting in better localization performance.

Contribution

It introduces BCE-AD, an extension of the BCE framework that incorporates metadata into the uncertainty estimation process for more accurate state estimation.

Findings

BCE-AD significantly improves localization accuracy.

Augmenting data with metadata enhances measurement uncertainty modeling.

Experimental results validate the superiority of BCE-AD over existing methods.

Abstract

The requirement to generate robust robotic platforms is a critical enabling step to allow such platforms to permeate safety-critical applications (i.e., the localization of autonomous platforms in urban environments). One of the primary components of such a robotic platform is the state estimation engine, which enables the platform to reason about itself and the environment based upon sensor readings. When such sensor readings are degraded traditional state estimation approaches are known to breakdown. To overcome this issue, several robust state estimation frameworks have been proposed. One such method is the batch covariance estimation (BCE) framework. The BCE approach enables robust state estimation by iteratively updating the measurement error uncertainty model through the fitting of a Gaussian mixture model (GMM) to the measurement residuals. This paper extends upon the BCE approach by arguing that the uncertainty estimation process should be augmented to include metadata (e.g., the signal strength of the associated GNSS observation). The modification of the uncertainty estimation process to an augmented data space is significant because it increases the likelihood of a unique partitioning in the measurement residual domain and thus provides the ability to more accurately characterize the measurement uncertainty model. The proposed batch covariance estimation over an augmented data-space (BCE-AD) is experimentally validated on collected data where it is shown that a significant increase in state estimation accuracy can be granted compared to previously proposed robust estimation techniques.

Tables

Table 1. TABLE I : Horizontal RSOS error statistics with green and red entries corresponding to the minimum or maximum statistic, respectively.

(m.)	$L_2$	DCS	MM	BCE	BCE-AD
median	9.84	10.82	9.13	6.70	2.14
variance	0.33	0.21	0.37	0.09	0.16
max	14.84	16.10	14.11	11.75	7.79

Table 2. (a) Horizontal RSOS positioning error results for data collect 1.

(m.)	$L_2$	DCS	MM	BCE	BCE-AD
median	9.84	10.82	9.13	6.70	2.14
variance	0.33	0.21	0.37	0.09	0.16
max	14.84	16.10	14.11	11.75	7.79

Table 3. (b) Horizontal RSOS positioning error results for data collect 2.

(m.)	$L_2$	DCS	MM	BCE	BCE-AD
median	5.09	5.02	4.11	3.58	3.58
variance	613.13	342.92	673.50	393.32	393.32
max	127.29	98.05	132.49	103.64	103.64

Table 4. (c) Horizontal RSOS positioning error results for data collect 3.

(m.)	$L_2$	DCS	MM	BCE	BCE-AD
median	10.98	9.88	8.02	7.31	3.25
variance	3.06	1.49	3.36	2.55	3.25
max	17.26	17.30	15.06	15.35	12.64

Table 5. TABLE II : Utilized features for each iteration of optimization with the BCE-AD approach. The green and red cell entries correspond to the utilized and non-utilized features, respectively. Within the table, ρ 𝜌 \rho is the pseudorange residual, Φ Φ \Phi is the carrier-phase residual, El. is the elevation angle, Az. is the azimuth angle, and SS. is the signal strength.

	Iter 0.	Iter 1.	Iter 2.	Iter 3.	Iter 4.	Iter 5.	Iter 6.	Iter 7.	Iter 8.
$\rho$
$\Phi$
El.
Az.
SS.

Table 6. (a) Utilized features per iteration of optimization for data collect 1.

	Iter 0.	Iter 1.	Iter 2.	Iter 3.	Iter 4.	Iter 5.	Iter 6.	Iter 7.	Iter 8.
$\rho$
$\Phi$
El.
Az.
SS.

Table 7. (b) Utilized features per iteration of optimization for data collect 2.

	Iter 0.	Iter 1.	Iter 2.	Iter 3.	Iter 4.	Iter 5.	Iter 6.
$\rho$
$\Phi$
El.
Az.
SS.

Table 8. (c) Utilized features per iteration of optimization for data collect 3.

	Iter 0.	Iter 1.	Iter 2.	Iter 3.	Iter 4.
$\rho$
$\Phi$
El.
Az.
SS.

Equations

$$\hat{X}=\underset{X}{\mathrm{argmax}}p(X|Y)$$

$$p(X|Y)\propto \prod _{n=1}^N{\psi }_n(A_n,B_n),$$

$$\hat{X}=\underset{X}{\mathrm{argmin}}\sum _{n=1}^N||r_n(X)|{|}_{{\Lambda }_n}\text{s.t.}{r}_n(X)\triangleq y_n-h_n(X),$$

$$\hat{X}=\underset{X}{\mathrm{argmin}}\sum _{n=1}^N{w}_n(e_n)||r_n(X)|{|}_{{\Lambda }_n}\text{s.t.}{e}_n\triangleq ||r_{n-1}(X)|{|}_{{\Lambda }_n}$$

$$D\sim \sum _{m=1}^M{w}_m\mathcal{N}(D|{\theta }_m)\text{s.t.}{\theta }_m\triangleq \{{\mu }_m,{\Lambda }_m\},$$

$$\log p(D)=\log \int p(D,\theta ,Z)d\mathbf{Z}d\mathbf{\theta },$$

$$q(Z{)}_{t+1}=C_z\int q(\theta {)}_t\log p(D,Z|\theta )d\theta$$

$$q(\theta {)}_{t+1}=C_{\theta }p(\theta )\int q(Z{)}_{t+1}\log p(D,Z|\theta )dZ$$

Full text

I Introduction

The expected operation domain of autonomous robotic platforms is ever-increasing. Currently, these applications span the space from structured industrial applications to autonomous exploration of novel environments [1]. This expectation has lead to astonishingly complex robotic systems with stringent safety requirements.

One of the key components of such robotic systems is the ability to accurately infer the essential states of the system when provided with a set of information (i.e., state estimation [2]). When the provided information adheres to the *a priori *models, this problem is addressed through numerous optimization frameworks (i.e., Kalman filtering [3], particle filtering [4], or graphical methods [5]). However, many scenarios emit characteristics that do not adhere to the specified assumption (e.g., the accuracy of the *a priori *models cannot be guaranteed in an novel operating environment). In such situations, the solution provided by any of the mentioned estimation frameworks can be arbitrarily biased by any single unmodeled observations111This can be shown through the calculation of the asymptotic breakdown point of the $l^{2}$-based estimators (i.e., any estimation framework that utilizes the $l^{2}\text{-norm}$ exclusively to construct the objective function) [6]..

To combat the breakdown of the $l^{2}\text{-based}$ estimators, research has been conducted in robust state estimation (i.e., state estimation methodologies that enable an increased breakdown point [6] when compared to traditional $l^{2}\text{-based}$ estimation). These robust estimation frameworks can be categorized into two paradigms: data weighting techniques, and data exclusion techniques. The data weighting techniques (e.g., robust maximum likelihood estimator [7], dynamic covariance scaling (DCS) [8], and max-mixtures (MM) [9]) enable robustness by reducing the influence (i.e., increase the corresponding uncertainty) of data which does not adhere to the *a priori *models. The data exclusion techniques (e.g., random sample consensus (RANSAC) [10], receiver autonomous integrity monitoring (RAIM) [11], and $l^{1}\text{-relaxation}$ [12]) enable robustness through selecting a subset of the data (i.e., the subset of data that adheres to the *a priori *models) and conducting state estimation with only the trusted subset.

All of the robust state estimation implementations described above share one key deficiency. Specifically, they all assume that the *a priori *uncertainty model accurately characterize the underlying system. This assumption can be problematic, for example, in adversarial environments that emit data degrading characteristics (i.e., environments where it is not feasible to accurately know *a priori *the measurement uncertainty model). To overcome this issue, the batch covariance estimation (BCE) approach was proposed within [13, 14]. The BCE approach works by iteratively updating the measurement uncertainty model between each iteration of optimization by fitting a Gaussian mixture model (GMM) to the measurement residuals. This approach has the benefit of learning the underlying system uncertainty model instead of assuming it from the outset.

The work presented within this paper provides an extension to the BCE approach. Specifically, this work removes the assumption that the measurement uncertainty model can be accurately characterized exclusively utilizing information contained within the measurement residuals. Instead, this work argues that the uncertainty estimation problem can effectively be augmented to incorporate metadata, leading to a more accurate characterization of the observation uncertainty model. As shown in Section IV the increase in positioning accuracy achieved by including metadata in the residual classification procedure can be significant.

The remainder of this paper proceeds in the following manner. In section II, a succinct overview of nonlinear least squares (NLLS) is provided, with a specific emphasis placed on robust estimation. In Section III, the batch covariance estimation over an augmented data-space (BCE-AD) robust estimation methodology is discussed. In Section IV the proposed approach is validated with multiple kinematic global navigation satellite system (GNSS) data-sets. Finally, the paper concludes with final remarks and proposed future research efforts.

II State Estimation

In this section, we briefly review of state estimation and its robust variants. For a more thorough discussion, the reader is referred to Section II of [14]. To begin the discussion, the state estimation problem can generically be defined as the process of calculating a set of states $X$ that – in some sense – are in best agreement with the provided information $Y$. Within this work, the metric utilized to quantify agreement is the maximization of the posterior distribution (i.e., the maximum a posteriori (MAP) state estimate $\hat{X}$), as presented in Eq. 1.

[TABLE]

The implementation a MAP estimator can be achieved through the utilization of the factor graph [5] formulation. The factor graph is a probabilistic graphical model framework which enables the factorization of the posterior distribution into a product of functions that operate on a reduced domain, as shown in Eq. 2

[TABLE]

where, $\psi_{n}(A_{n},B_{n})$ is an application specific domain reduced function (i.e., a factor in the factor graph model), which operates on $A_{n}\subseteq\{X_{1},X_{2}\ldots,X_{n}\}$, and $B_{n}\subseteq\{Y_{1},Y_{2}\ldots Y_{m}\}$.

To facilitate a computationally efficient implementation, it is commonly assumed that each factor within the models adheres to a Gaussian noise assumption. When this assumption is enforced, the estimation problem presented in Eq. 1 is reduced to finding the set of states which minimizes the squared sum of weighted residuals [15], as presented in Eq. 3

[TABLE]

where $r_{n}(X)$ is an observation residual, $h_{n}$ is a function that maps the state estimate to the observation domain, $\Lambda_{n}$ is a covariance (i.e., residual weighting) matrix, and $||\ *\ ||$ is defined as the $l^{2}\text{-norm}$.

II-A Robust State Estimation

Based upon the discussion provided within [14], we can view robust state estimation schemes as variants of iteratively re-weighted least squares (IRLS) [16]. The implementation of the IRLS formulation is presented in Eq. 4, where $w_{n}$ is a real valued function that attempts to appropriately weight observations based upon the previous iteration’s residuals.

[TABLE]

The specific implementation of weighting functions has been extensively explored. One such implementation, commonly utilized within the robotics community, is the DCS approach [8] which implements a redescending222A weighting function that approaches [math] as the magnitude of the error approaches $\infty$ (i.e., $\lim_{e_{n}\to\infty}w(e_{n})=0$)

m-estimators type weighting function. Another commonly utilized framework is the MM approach [9], which enable the utilization of a GMM error uncertainty model through the approximation of the summation operation with the maximum operation333When a GMM is utilized to characterize the error uncertainty model, the objective function of the corresponding estimation problem does not reduce to a NLLS form, which increases computation complexity. To reduce this additional complexity, the maximum operation can be utilized in place of the summation operation in the GMM (i.e., $\operatorname{p}(y_{i}\ |\ X)=\sum_{m=1}^{M}w_{m}\mathscr{N}(\mu_{m},\Lambda_{m})\approx\max_{i}w_{i}\mathscr{N}(\mu_{i},\Lambda_{i}))$.. Finally, the MM approach was extended within [13, 14, 17] to allow for non-parametrically learned uncertainty models (i.e., the uncertainty model is not assumed to be static) based upon the measurement residuals in the BCE approach.

III Proposed Approach

The estimation framework proposed in this paper is an extension of the BCE methodology as detailed within [14]. In this section we have two primary objectives. Firstly, to expound the differentiating factors between the two implementations. Secondly, to provide an overview of the proposed framework entitled BCE-AD.

III-A The Data Model

To initiate a discussion of the proposed estimation framework, we will start with the assumed data model. In our previous work [13, 14], it is assumed that a given set of residuals $R=\{r_{1},r_{2},\ldots,r_{n}\}$ with $r_{n}=y_{n}-h_{n}(X)\in\mathbb{R}^{d}$ can be accurately partitioned into groupings of similar instances. To depict this visually, Fig. 1(a) provides the GNSS measurement residuals (i.e., the pseudorange and carrier-phase residuals) for a typical localization application. As illustrated, there is no obvious partitioning of the set.

To increase the likelihood of a unique partitioning, we can project the residual dataset into a higher-dimensional space [18]. The specific projection utilized within this study is the augmentation of the original dataset with a set of metadata (i.e., for each calculated residual, there is an additional set of features $F=\{f_{1},f_{2},\ldots,f_{n}\}$ with $f_{n}\in\mathbb{R}^{f}$). Where a feature considered for inclusion in the metadata-set is an observed quantity that is known to correlate to the quality of the collected sensor observation444For GNSS applications, signal strength and elevation angle of a GNSS observation are commonly utilized. Other sensor have similarly useful features (e.g., the mean illumination of an image for computer vision applications)

. To depict the benefit granted by the incorporation of metadata, Fig. 1(b) provides the carrier-phase residuals augmented by two additional features (i.e., elevation angle and signal strength of the collected observation). From this figure, it can be seen that the inclusion of metadata provides a more obvious partitioning.

Given this set of augmented data $D=\{d_{1},d_{2},\ldots,d_{n}\}$ with $d_{n}=\{r_{n},f_{n}\}\in\mathbb{R}^{d+f}$. We will assume that the augmented set can be partitioned into similar groupings (i.e., $\bigcup_{m=1}^{M}C_{m}=D$), where each group, $C_{m}$, can be characterized by a Gaussian distribution (i.e., $C_{m}\sim\mathscr{N}(\mu_{m},\Lambda_{m})$). With this assumed model, the augmented dataset is fully characterized as a GMM, as depicted in Eq. 5

[TABLE]

where, $m$ is the number of components in the mixture model555The number of mixture components $M$ in the estimated model is not assumed to be known *a priori *. Instead, a truncation level $M^{*}$ is set [19] (i.e., a maximum number of components, $M^{*}\geq M$) and the number of utilized components is autonomously selected based upon the variational free energy of the estimated model [20].

, $w$ is the set of mixture weights with the constraint that $\sum_{m}w_{m}=1$, and $\theta_{m}$ is the mixture components sufficient statistics.

III-B Variational Clustering

As discussed within [14], to fit the a GMM to a provided dataset there are two broad classes of algorithms. There are the sampling based frameworks (i.e., Monte Carlo Based approaches [21]), and the optimization based frameworks (i.e., the variational based approaches [22]). Within this study, the variational clustering framework [22, 20] is utilized to reduce the computational complexity of the fitting process.

Specifically, the variational clustering approach is utilized to estimate the model parameters that maximize the $\operatorname{log\ marginal\ likelihood}$ of the data, as provided in Eq. 6

[TABLE]

where $Z=\{z_{1},z_{2},\ldots z_{m}\}$ with $z_{m}\in\mathbb{R}^{M}$ is a set of assignment variables (i.e., Z provides the explicit assignment of each data instance to a component within the GMM). To enable the tractable computation of the GMM, the mean-field assumption [23] (i.e., $\operatorname{p}(D,\theta,Z)\approx\operatorname{q}(\theta)\operatorname{q}(Z)$), is utilized to construct a lower-bound on the true $\operatorname{log\ marginal\ likelihood}$.

With the mean-field assumption in place, the GMM parameters can be estimated with the variational framework in a iterative fashion, where the sequence of iteration is presented in Eqs. 7 and 8. The assignment parameters $Z$ are updated by optimizing Eq. 7, where $C_{z}$ is the normalizing constant to $q(z)$. Then model parameters $\theta$ are updated by holding the assignment parameters fixed and optimizing Eq. 8, where $C_{\theta}$ is the normalizing constant for $q(\theta)$. This iterative process is continued until the Kullback-Leibler divergence between the approximating and true distributions is minimized.

[TABLE]

[TABLE]

III-C Feature Selection

The inclusion of additional features in the assumed data model have the benefit of increasing the likelihood of a unique partitioning in the measurement domain. However, clustering over an augmented dataset (i.e., a dataset with increased dimensionality) could have the detrimental side effect of increased computation complexity. Thus, to remain computational tractable, the set of augmenting data must be intelligently selected.

The enabling frameworks for intelligent feature selection can be classified as either offline, or online. The offline approach can either be implemented algorithmically of through the utilization of an area expert (i.e., someone with extensive experience with the utilized sensing modality). A primary drawback of the offline feature selection implementation is a static feature space model (i.e., the utilized features are provided *a priori *and remain fixed over the duration of the estimation process).

Due to the undesirable properties of an offline feature selection approach, an online framework is utilized within this study. Specifically the framework developed within [24] is utilized to autonomously select the most relevant features between each iteration of optimization. In brief, the utilized feature selection algorithm can be described through three primarily steps: 1) construct a nearest neighbor (NN) graph from the provided augmented dataset, 2) conduct eigenvalue decomposition on the NN graph to measure the importance of each feature for partitioning the dataset, 3) perform least angle regression (LARS) [25] on the calculated eigenvectors, to find the most important features. For a more thorough discussion on the feature selection algorithm, the reader is referred to [24].

III-D Algorithm Overview

With the previously detailed topics, the discussion can now proceed to an overview of the proposed estimation framework. From Fig. 2, it is shown that the proposed algorithm is comprised of two primary segments. The first component is the initialization of the estimator. This process begins by constructing the factor graph representation of the NLLS optimization problem – as thoroughly detailed within [5] – from the a priori state and uncertainty information and the provided observations. Then, an initial iteration of optimization is conducted with a NLLS estimation algorithm (e.g., for this study, the Levenberg-Marquardt [26] implementation was utilized) to update the *a priori *state estimate.

The second component of the proposed framework commences with the calculation of the measurement residuals given the previously estimated set of states. The set of calculated residuals are augmented with the provided set of metadata, as discussed in section III-A. Using the augmented dataset for the current iteration of optimization, the most relevant features are selected using the framework discussed in section III-C. Using the relevant set of features, the variational cluster framework – as disused in section III-B – is utilized to assign each instance in the augmented dataset to a component within an estimated GMM.

Utilizing the partitioning of the augmented data space, the measurement error uncertainty model can be calculated. This is achieved by first partitioning the measurement residuals with the assignment vector $Z$ estimated in the augmented data space. Then, the sufficient statistics for each grouping in the measurement residual domain are calculated. With this $m\text{-dimensional}$ GMM in the measurement residual domain, the measurement uncertainty model of the factor graph is updated (i.e., each measurement’s uncertainty model is updated to the sufficient statistics of the assigned GMM component). With the updated measurement uncertainty model, a new iteration of optimization is conducted. This process is iterated until a measure of convergence (e.g., the error decrease between consecutive iterations is less than a user defined threshold) – or a limit on the number of iterations – has been reached.

IV Experimental Validation

IV-A Data Collection

To enable the validation of the proposed estimation framework, three kinematic GNSS data collects666These three GNSS data collects were made publicly available within [14]. For a detailed description of the data collection procedure, the reader is referred to Section IV of [14] are utilized. The ground-trace for the three driving data collects is provided in Fig. 3. For each of the depicted data collects, the GNSS binary in-phase and quadrature (IQ) observations were recorded with a LabSat-3 GPS record and playback device [27]. The collected IQ data were then played back into two GNSS receivers: a geodetic-grade (Novatel OEM-638) GNSS receiver, and an open-source global positioning system (GPS) software defined radio (SDR).

The observations generated by the geodetic-grade receiver, in conjunction with the observations generated by an additional static GNSS receiver, were utilized to calculate a reference positioning solution for each data collect. The reference solution was calculated through a differential (i.e., real time kinematic (RTK)) filter-smoother framework implemented within the open-source package RTKLIB [28].

The observations generated by the GPS SDR are utilized by the estimation framework to validate the proposed methodology. These observations are utilized for validation because they were intentionally degraded by altering the GNSS receiver’s tracking parameters. For a more detailed discussion on the specific GNSS tracking parameter settings, the reader is referred to Section IV of [14].

With these three datasets, the proposed approach (i.e., BCE-AD) is validated against four additional estimators. The first comparison algorithm is a traditional (i.e., non-robust) estimation framework that utilizes a $l^{2}\text{-norm}$ cost function. Additionally, the DCS robust estimation framework [8] is utilized. The third comparison algorithm is the MM estimation framework [9] with a static measurement uncertainty model. The final comparison algorithm is the BCE approach [14].

To enable the utilization of the BCE-AD approach, the set of features $F$ need to be specified. For GNSS applications, several useful metadata can be calculated for each observation. For this evaluation, it was elected to use three such features: the signal strength (SS), the elevation angle (EL), and the azimuth angle (AZ). These features were selected for two primary reasons: 1) these features can be easily extracted for each GNSS observation, 2) these features are known to correlate to the quality of the recorded GNSS signal.

IV-B Results

To begin the evaluation of the proposed methodology, the horizontal residual-sum-of-squares (RSOS) positioning error for the three data collects will be examined. This evaluation is presented visually in Fig. 4, in the form of a box-plot. From Fig. 4 it is shown that the proposed BCE-AD significantly reduces the median horizontal RSOS positioning error for data collects 1 and 3, and performs comparably well to the other robust estimators on data collect 2. The specific statistics representing the analysis presented in Fig. 4 are provided in Table I(c).

To examine the reason why the BCE-AD approach provided no additional positioning performance benefit when compared to the BCE approach for data collect 2, the measurement domain partitioning for the two approaches can be evaluated. This evaluation is provided in Fig. 5 for all three data collects. From the provided figure, it can be hypothesized that no additional positioning accuracy is granted when the measurement residuals can be easily partitioned utilizing exclusively the information in the measurement residual domain. This is depicted visually by evaluating Fig. 5(c), where it is shown that the clustering utilizing only the information in the measurement domain can easily partition the provided measurement residuals. However, when evaluating the measurement residuals for data collects 1 and 3 (i.e., the datasets where the greatest positioning performance benefit is granted by the BCE-AD approach), it is apparent that clustering over the two different domains provides significantly different partitionings.

Finally, the utilized features with the BCE-AD approach for each dataset are evaluated. This evaluation is depicted visually in Table II(c). From this visual, it can be noted that the most relevant features not only vary from dataset to dataset, but also from one optimization iteration to the next within a given dataset. Additionally, the hypothesis presented in the previous paragraph (i.e., that the BCE-AD only grants increased positioning performance when sufficient information is not present in the measurement residual domain to partition the residuals) is further verified in Table II(b) where it is shown that only features utilized for data collect 2 are the measurement residuals.

V Conclusion

This paper presents an extension of the previously proposed batch covariance estimation (BCE) technique to enable robust state estimation. The BCE approach enables robust state estimation through the iterative estimation of a measurement error uncertainty model based upon the previous iterations measurement residuals. Where, the estimated measurement error uncertainty model is characterized by a Gaussian mixture model (GMM) which is fit to the measurement residuals through variational clustering. After fitting the GMM to the current optimization iterations residuals, the uncertainty model of each observation is updated to the sufficient statistics of the assigned cluster within the GMM.

The approach proposed within this work extends the BCE approach on one front. Specifically, it removes the assumption that the measurement error uncertainty model can be accurately characterized utilizing information exclusively contained within the measurement residual domain. Instead, this paper argues that the uncertainty estimation process should be augmented to include additional metadata. The modification of the uncertainty estimation process to an augmented data space increases the likelihood of a unique partitioning in the measurement residual domain and thus provides the ability to more accurately characterize the measurement uncertainty model.

To verify the proposed batch covariance estimation over and augmented data space (BCE-AD) approach, three GNSS data sets were utilized. The utilized data sets provide varying levels of degradation to quantify the robustness of the proposed algorithm against other state-of-the-art robust estimators. Utilizing these data sets, it is shown that the proposed approach provides comparable or improved state estimation accuracy when compared to other robust estimation techniques. In conclusion, it should be noted that, while the approach was validated on GNSS observations for this study, the proposed approach is generically extensible to other estimation domains.

Acknowledgment

The authors would like to thank MacAuly-Brown Inc. for funding the work presented within this article.