Improving Real-Time Position Estimation Using Correlated Noise Models †

We provide algorithms for inferring GPS (Global Positioning System) location and for quantifying the uncertainty of this estimate in real time. The algorithms are tested on GPS data from locations in the Southern Hemisphere at four significantly different latitudes. In order to rank the algorithms, we use the so-called log-score rule. The best algorithm uses an Ornstein–Uhlenbeck (OU) noise model and is built on an enhanced Kalman Filter (KF). The noise model is capable of capturing the observed autocorrelated process noise in the altitude, latitude and longitude recordings. This model outperforms a KF that assumes a Gaussian noise model, which under-reports the position uncertainties. We also found that the dilution-of-precision parameters, automatically reported by the GPS receiver at no additional cost, do not help significantly in the uncertainty quantification of the GPS positioning. A non-learning method using the actual position measurements and employing a constant uncertainty does not even converge to the correct position. Inference with the enhanced noise model is suitable for embedded computing and capable of achieving real-time position inference, can quantify uncertainty and be extended to incorporate complementary sensor recordings, e.g., from an accelerometer or from a magnetometer, in order to improve accuracy. The algorithm corresponding to the augmented-state unscented KF method suggests a computational cost of O(dx2dt), where dx is the dimension of the augmented state-vector and dt is an adjustable, design-dependent parameter corresponding to the length of “past values” one wishes to keep for re-evaluation of the model from time to time. The provided algorithm assumes dt=1. Hence, the algorithm is likely to be suitable for sensor fusion applications.


Introduction
Global Positioning System receivers output a time-series of position measurements, but this signal suffers from errors due to many reasons, such as atmospheric and multipath effects [1][2][3] or insufficient coverage of satellite constellations. Rather than providing a full quantification of uncertainty, GPS devices generally provide Dilution-Of-Precision (DOP) variables, dimensionless multiplicative factors of the position error, designed to flag bad satellite coverage. In this paper, we develop we formulate the scoring rules to evaluate algorithm performance. In Section 5, we present the results of the inferences and uncertainty quantification. We conclude in Section 6.

Data Analysis
The data on which we construct and test our methods were recorded by a single-frequency GlobalSat BU-353S4 unit operating on the L1 frequency, 1575.42 MHz [17]. We used this GPS receiver, fixed in position at one location in Dunedin, New Zealand, to record seven time-series sampled every second over time scales between approximately one and 25 days between the 6 August and the 24 November 2014. The same receiver recorded four additional time-series while fixed in position in three locations in Australia: North Stradbroke Island, Port Douglas and Thursday Island (see Figure 1). These additional time-series were also sampled every second and had time scales between approximately one and three days. The variety of locations and recording dates allowed us to check for possible geographical and temporal variation in the GPS signal, e.g., due to the varying position of GPS constellations. Additional data were taken at the same times and locations using tag-style GPS units, which also used the L1 frequency. This additional data had the same noise characteristics as the original data [8], showing that the noise properties, which we analyse below, are not dependent on the particulars of the GPS receiver, but rather attributed to environmental effects. The results in this paper are therefore expected to be applicable to any single-frequency GPS device operating on the L1 frequency. Map of GPS receiver locations for the datasets considered. A single GPS receiver was located at different times at the labelled locations: Dunedin (New Zealand) and North Stradbroke Island, Port Douglas and Thursday Island (Australia). The map demonstrates the spread of the four testing sites, indicating that the observed correlation patterns cannot be attributed solely to an unfortunate "bad" geographic location. The partial autocorrelation functions,ρ p (h), for latitude, longitude and altitude are shown on the left for the four geographical locations as a function of the lag, h, measured in seconds.

Noise Correlations
The datasets were analysed in detail in Reference [8], Here, we briefly present the findings relevant for this study. When positioned in a new location, the GPS device gives readings offset from the correct location over a distance scale of approximately 100 m. The readings tend to the correct location within a time scale of seconds to minutes. In this paper, all further analysis and inference is performed on datasets with the first 300 s truncated in order to ignore this initial transient behaviour. Subsequently, the readings take continual excursions with typical length scales of metres around what we assume to be the "true" location (see Figures 2 and 3). These excursions have typical time scales of minutes to hours, causing the data to be autocorrelated [8]. Note that the GPS data are recorded in degrees to eight significant figures; these data are converted to meters, which, at the locations analysed in this paper, leads to an apparent quantisation of approximately 0.1 m. This quantised length scale is larger than the typical length scales for the processes, and thus, the data apparently jump between these quantised locations in tens to hundreds of seconds. None of the further treatment in this paper considers this apparent quantisation, which may be an avenue for further work. Figure 2. (a-c) First 600 s of GPS time-series residuals taken at Port Douglas, illustrating the rapid initial convergence towards the "true" value; (d-f) corresponding later GPS time-series residuals illustrating the time-correlated deviations from the "true" value (zero). The apparent position quantization due to rounding is also visible.

Noise Models
We studied and evaluated the performance of various noise models for GPS data [8]; specifically, we focused on uncorrelated noise, OU, autoregressive processes (x t = − ∑ p k=1 h k x t−k + E t ), moving-average processes (x t = µ − ∑ q =0 g E t− ) and their combination, an ARMA(p = 1, q = 1) Here, x t represents the signal at time-step t, while E t denotes the noise at the same time-step. The coefficients h 1 , g 0 and g 1 are constants specific to the system under investigation. Using the standard Akaike Information Criterion (AIC) [18], we compared these model.
One of the conclusions of the comparison was that all noise models proved to be better than the the uncorrelated noise model. In this paper, we continue our investigations and compare techniques including or omitting correlation in noise models. An OU process is also included in these analyses, because it has the following beneficial properties: mean-reverting, continuous and having only a few parameters; hence, it lends itself to fast embedded computation and sensor fusion applications.
The OU process is defined by the stochastic differential equation: where W t represents a Wiener process, while θ, µ and σ are parameters to be fitted. The associated Fokker-Planck equation has an analytical solution [11]: Here, N µ n , σ 2 n is the Gaussian distribution with mean µ n and variance σ 2 n . This probability distribution for x t is initially a delta function located at x = x 0 , with the limit, N (µ, σ 2 /(2θ)), as t → ∞. The process is mean-reverting for θ > 0; i.e., its limit is a stationary solution to the Fokker-Planck equation with mean µ, as t → ∞. The time-discretised version of (1) for time-step ∆t is given by: where t , the noise term, follows a Gaussian distribution with zero mean and variance σ 2 2θ 1 − e −2θ∆t . During data processing, we assumed that the long-term average represents the "true" value for a random variable; hence, after subtracting this long-term average, the residual data must have zero mean, µ = 0. Hence, from each residual, we obtained maximum likelihood estimates σ 2 and θ. For each coordinate, the logarithms θ and σ 2 were independently modelled, where the coordinate type (latitude, longitude and altitude) was assumed to have a fixed effect, while the geographic location and time have random effects. The geographic location seemed to have no significant effect on θ, but it weakly influenced σ 2 ; though, we believe, that this effect has no practical significance. The type of coordinate, on the other hand, seems to affect both θ and σ 2 . Parameter θ was statistically significantly smaller for the altitude than for the other two coordinates; simultaneously, σ 2 was significantly larger for altitude, than for latitude or longitude. These results essentially confirm the conclusions of a visual inspection of Figure 4. Other models: higher order Autoregressive (AR), Moving-Average (MA) and, their combination, ARMA processes [8] are outside the scope of the current paper, but are promising candidates for the extension of this study (see Section 6). Qualitatively, this can be understood by the larger variability of the altitude residuals (see Figure 3).

Position and Uncertainty Estimation Methods
Before introducing Bayesian inference methods, we create a basic model that uses the raw position as the centre of a prediction interval, while the OU process determines the interval width. Specifically, our prediction for latitude, longitude and altitude are Gaussian distributions with variances σ 2 /θ, where σ and θ 2 are the means of the fitted OU parameters for all three coordinates observed in Dunedin (see Section 2.2). We consider this model basic, because previous data do not inform later position estimates, hence cannot take into account the noise correlations. By using the parameters derived from the Dunedin time-series for the uncertainty quantification on the Australian time-series, we tested that the method works for datasets independent of those on which it was developed.
Receivers also often report DOP figures, horizontal and vertical, for the latitude, longitude and altitude coordinates. These numbers, being related to the diagonal elements of the covariance matrix, carry some limited information about time correlation. We incorporated this complementary data-stream into the statistical modelling. In Section 5, we present the effect of the HDOPand VDOPnumbers of the parameter estimates. If these estimates improve, DOP parameters are worthwhile to include in the statistical modelling of the noise.

Earlier Work
Since the late 1990s, few researchers outlined alternative ways to estimate distributional parameters of time-series of GPS signals and thereby analysed the noise inherent in these signals [19][20][21][22]. Johnson and Agnew even demonstrated, although for crustal velocities, that temporal correlation in observation can decisively influence quantities estimated from the measurement [23]. In general, the temporal correlation reduces the uncertainties of the geophysical parameters inferred from the signal; thus, one can be more certain of a potentially incorrect estimate if the temporal correlation is not considered. Most approaches relied on Maximum Likelihood Estimators (MLE) either in fitting ad-hoc models to the data or in estimating the covariance matrix of the signal components [19,[24][25][26]. This MLE approach was considered to provide robust results [27] compared to linear regression. Similar to Johnson's earlier study, Santamaría-Gómez et al. [28] and Masson et al. [29] concluded that incorrectly modelled noise in synthetic time-series introduces bias in the velocity data, and even worse, increasing the length of the time-series does not guarantee diminishing bias.
While some steps have been made to depart from the use of the white noise model, e.g., replacing it with time-correlated (i.e., coloured) noise [19,30], it is still often assumed that the noise can be modelled by independent, identically distributed Gaussian random variables. Bos et al. [21] developed a technique that substantially increases the efficiency of the MLE methods, which has since been improved [25]. Other, more advanced and computationally costly methods, e.g., Markov chain Monte Carlo methods [22,31] or wavelet analysis [32], have also been proposed to obtain non-biased probability distributions for the noise within the time-series.

Kalman Filter with Higher Order Noise
We formulate sophisticated methods that use Bayesian inference of the position and its uncertainties. Kalman filters [4] are a type of sequential Bayesian inference algorithms, often useful for real-time state estimation. An improved version, the Unscented Kalman Filter (UKF), allows state estimations on nonlinear systems by approximating the prior distribution by well-chosen "sigma" points [33,34]. The UKF, like other Kalman filters, assumes that a state, x n , changes between two consecutive measurements, and this change can be described f : x n → x n+1 , called the forward map. Furthermore, the measured quantities fully determine the state x, i.e., g : x →ŷ. Since one's knowledge of the system of interest is incomplete, one often models it by a stochastic process, in which a random variation, known as process noise, is also added to the forward map. Likewise, g may also be contaminated by noise, representing random measurement errors. The Augmented-state UKF (AUKF) is a variant of the UKF, which describes the measurement and process noise distributions by "sigma" points [33][34][35][36][37]. The AUKF is described in Algorithm A1 located in Appendix A. Most Kalman filter schemes require the process and measurement noise to be additive; however, within the AUKF, the forward map and observation map are allowed in any form. The AUKF is therefore sufficiently flexible to handle quite general processes. Such processes include simultaneous estimation of OU parameter θ, which we can constrain to be positive by estimating its logarithm, and possibly moving-average processes [8,33,[38][39][40], although moving-average processes would require a different choice of sigma points from those described in Algorithm A1. The AUKF will also be able to handle nonlinear forward maps due to some complicated sensor dynamics.

Uncorrelated Gaussian Noise
For a GPS measurement with iid Gaussian noise (both in different coordinates and in time), the system state is represented by one position coordinate for each spatial direction considered (latitude, longitude and altitude). We, therefore, treat the coordinates separately: x = x; hence, the forward map, f , is the identity with vanishing process noise. The measurement noise is considered to be additive and uncorrelated: g(x) = x + , with ∼ N (0, Σ v ), and Σ v is associated with a stationary OU distribution Σ v = σ 2 /θ. Here, σ 2 and θ are the same as for the non-learning method above.
At the start of the data processing, we initialise the prior distribution with state mean µ 0 and covariance K 0 . We chose µ 0 to be 5 m from the "true" position, and the variance was taken to be K xx,0 = 20 m 2 . This choice was motivated by the fact that for a Gaussian prior distribution, the true mean value lies near the edge of the central 90% probability region, since K xx,0 ≈ 4.47 m. An inference algorithm, we expect, would lead to a narrower posterior distribution around the true value.
The above procedure is used generally if a model diverges from the "true" value even if evolved by the exact dynamics. In such a case, some additive process noise is incorporated into the inference process. For the iid Gaussian model, the usage of such complementary process noise is equivalent to a Brownian motion around the "true" position. In order to find an optimal value for the process noise, we vary its magnitude and carry out inference with these Brownian models. Although these Brownian models contain correlated noise, their performance is inferior to a well-constructed OU model.

Ornstein-Uhlenbeck Noise
Correlated noise can be included in the Kalman filter schemes by extending the state vector to include one or more parameters to represent the noise model. When an OU process generates the noise in the state variable, x is augmented by a new variable, x t , determined by the OU process. The corresponding forward and measurement maps are x t → (x t − µ) exp (−θ∆t) + , where ∼ N 0, σ 2 1 − e −2θ∆t /2θ and g(x) = x t . The noise is thus shifted to the process noise from the observation noise (cf. Section 3.2.1); hence, the former can now be treated as inconsequential. Specifically, we choose a mean position coordinate µ and variance K µµ,0 with the same values as x and K xx,0 in the prior for the uncorrelated noise model. The state also includes the OU variate x; we chose a prior for x equal to the first data point of the GPS time-series and therefore started the filtering from the second data point onwards. Consequently, we chose a small value for the prior variance in the noise parameter K xx,0 = 0.1 m 2 . We also included the reversion rate parameter θ (or, more specifically, its logarithm) in the state, with the prior mean/variance of θ equal to the mean/variance of the fitted parameters for the Dunedin GPS datasets for latitude, longitude and altitude, respectively (see Section 2.2). Therefore, the inferences for the North Stradbroke Island, Port Douglas and Thursday Island time-series (which had no role in the parameter fitting) serve as a test of the general applicability of the algorithm in independent locations.
It is not possible to estimate the OU variance parameter, σ 2 , using orthogonal sigma point selection for which the covariances of the state parameters with the process noise are set to zero (see Algorithm A1, Step 2). Hence, here, σ 2 was set for each coordinate as the mean of the fitted σ 2 for the Dunedin GPS datasets. The simultaneous estimation of σ 2 with the state parameters is left for future research.

Quantifying the Filtering Algorithm Performance
Similar to the method in Section 3.2.1, which issues time-series estimates of position in the form of Gaussian distributions, the UKF produces {x n } together with {Σ n }. We interpret these as the mean and variance parameters of Gaussian posteriors. With this interpretation, the marginal distributions {q(x)} of position are also Gaussian.
To compare and rank algorithms with different noise models according to their performance, we introduce the log-score [15,16,41] ( = 1 m): This quantity is often interpreted as a "surprise" of true value x. If q is large around the true value, the corresponding log-score, S, is low, i.e., the outcome is not surprising. An algorithm producing lower log-scores will be regarded as superior compared to those achieving higher S values.

Results
We analyse the predictions of the methods described in Section 3. Figure 5 shows the estimates for GPS data collected in Port Douglas. We note that this dataset is independent of the Dunedin datasets on which the model parameters were based. Port Douglas is a representative example, and the discussion below could be repeated for the other locations as well.
The uncertainty of the unprocessed location data is well captured by the Gaussian prediction obtained from the OU stationary distribution. Taking into account the reported DOP parameters led to small modifications only. Although occasional large spikes did appear in the uncertainties, these did not seem to correspond to large deviations in the position data. The time-series of the log scores of these predictions with and without DOP are very similar; cf. Figure 6d-f. Figure 7a,b shows the time-averaged log scores, S, for both methods to be very similar, which holds for all the locations tested.
The raw location data jump around the correct value, similarly to the marginals from the UKF inference obtained assuming a Brownian model. Due to the process noise, the Brownian model does not converge to the correct value. On the other hand, the position marginals from the UKF augmented with OU noise converge to the "true" value (see Figure 5d-f). The improvement in performance is supported by the log-scores. The posteriors produced from the OU noise score better by three orders of magnitude (see Figure 6). It is clear that the KF with the OU noise model is the superior method and the only method to produce negative scores, which continue to decrease with time. Figure 8 shows the scores for all the time-series from all locations averaged over the run-time. With respect to this average score, all KFs with OU noise outperform all other models using iid Gaussian noise models and all inferences using the non-learning method. The single exception is the altitude measurement for Thursday Island. This discrepancy is due to mischaracterised noise; σ 2 in this run was much larger than the one calculated from the Dunedin dataset. We thus have doubled σ 2 in the OU model for all the inferences. This manual intervention brought σ closer to the fitted value for the Thursday Island datasets, having negligible effect on the other inferences. We, therefore, suggest that extending the UKF to simultaneously estimate σ 2 during the inference process, in an embedded calculation, could allow the algorithm to adjust automatically to the variation of σ 2 with time and location.
As expected based on Figures 3 and 4, the inferences of altitude are worse than those of the latitude and longitude. We believe that the otherwise well-characterised noise provides close to optimal quantification of uncertainty for the altitude coordinate, especially if one takes into account the intrinsic character of the altitude time-series due to the relative position of the satellite constellation.
We include for completeness the θ marginals in Figure 9. The θ estimates from the Bayesian filtering have the same orders of magnitude as the fitted parameters (see Figure 4), with the θ for the altitude coordinate smaller than the others, as expected. However, the θ for the longitude is larger than that for the latitude, which is not expected from the fitted version-it is perhaps an artefact of the imperfect model fit.

Discussion and Future Work
We developed algorithms for the sequential inference of the location of a GPS receiver, taking into account the correlated nature of the GPS noise. We found that the DOP parameters do not help in uncertainty quantification of the GPS positioning. This noise, which may be attributable to atmospheric and multipath effects, as well as the satellite configuration, is found to be correlated in time and can be modelled as an OU process. The AUKF algorithm using an OU noise model provides an accurate series of posterior distributions for GPS position with improved uncertainty quantification compared to inference assuming iid Gaussian noise provided a large enough value is assumed for the OU noise variance. The improved posterior distributions offer many practical advantages in applications in which knowledge of the correct uncertainty is important, e.g., automated navigation, where one must be very confident of avoiding obstacles. The algorithm is also easily extendible to incorporate data from additional sensors, e.g., accelerometers and gyroscopes to inform the inference (so-called sensor fusion). We expect that using a correlated model for the GPS noise in a sensor fusion application will extend its potential to improve accuracy, as a well-quantified position uncertainty will allow position corrections from other sensors to be accepted.
There are clear avenues for further improvement of the techniques in this paper. An algorithm similar to Algorithm A1 may be used to simultaneously estimate the OU variance σ 2 during the inference process, by including σ 2 as a state parameter. However, such an algorithm requires a different choice of sigma points from that in Algorithm A1. Additional sigma point should be added to the scheme with mixed components from σ 2 and the process noise, such that the mean and co-variance of the sigma points does not change. Furthermore, since we have earlier showed that more sophisticated noise models, e.g., the higher order AR processes and ARMA processes could be superior to the OU model, we expect that upgrading the noise models will lead to further improved position inference. The AUKF algorithm presented could be straightforwardly upgraded to use the higher order AR processes, taking care to parametrise the model such that it represents a stationary process. The forward map of an ARMA process also includes multiplicative noise, so to use such a forward map, similar care must be taken with the choice of sigma points.
We mention an extension of our current work to engineering applications. The effect of time-correlated noise on the applicability of classical Kalman filters has been previously analysed to a certain degree [42]. The noise models proposed here, however, can be seamlessly incorporated into a Kalman filter or Bayesian inference algorithm providing both an estimate for the location and the associated uncertainty quantified. Although such extensions require extra numerical calculation, the computational load is acceptable and provides seamless integration for sensor fusion, treating all sensors on equal footing. Incorporating the aforementioned noise models into such algorithms is in preparation.
Finally, although for the sake of transparency and simplicity, we modelled latitude, longitude and altitude separately as one-dimensional time-series, a more natural approach could be to combine them into one vector quantity. Building an inference algorithm similar to that in Section 3.2 would lead to a compact model that calculates and takes into account the covariances between these coordinates. This could further improve the algorithm performance.

Conclusions
We collected location data with a stationary GPS sensor at four locations: Dunedin (New Zealand) and North Stradbroke Island, Port Douglas and Thursday Island (Australia). The durations of these time-series range from a day to 25 days. We found appreciable autocorrelation within the noise present in the signal; thus novel, noise-models were constructed and compared for performance. In order to rank the performance of different noise models, we used the log-score. We found that an enhanced Kalman filter augmented with the Ornstein-Uhlenbeck noise model performed the best, and its pseudocode is provided. The computational cost of the corresponding algorithm is modest, and it generalises to sensor fusion application easily; thus, it is promising for embedded computing.

Conflicts of Interest:
The authors declare no conflict of interest. The funding agencies had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results.