Bayesian Covariance Estimation for Kalman Filter based Digital Carrier Synchronization in GNSS-SDR

A project for Google Summer of Code 2018 in coordination with GNSS-SDR

Introduction to GNSS and Kalman Filtering

Modern mass-market GNSS software receivers rely on tracking architectures developed from analogue signal processing techniques developed in the early 20th century. Although these methods have proven themselves to be reliable and sufficient for most GNSS applications, newer applications which aim to push the boundaries of what can be done with GNSS could benefit from more versatile, modern, tracking techniques. The basis of many of these cutting edge techniques is the Kalman filter: an iterative signal processing technique which aims to estimate the internal parameters of a system based on noisy observations of its external measurements.

Within the context of GNSS, Kalman filtering has been shown experimentally to be a powerful tool for performing carrier synchronization, an important step in the GNSS datapath in which the receiver identifies the timing properties of the received signal (i.e. carrier phase and Doppler shift) so that the signal can be decoded to determine the transmitted message. In fact, the principle of carrier synchronization is present with any digital communication, but in the case of GNSS, carrier synchronization preforms two roles: not only are the timing parameters needed in order to decode the navigation data transmitted by the satellite, carrier phase and Doppler shift are also two of the observables required to get a position fix.

"Comparison of PLL-DLL tracking to tracking to 2-state Kalman filter tracking of GPS L1 C/A signal recorded at CTTC in 2013 under nearly ideal conditions. Note the highly improved signal-to-noise ratio in the Carrier Frequency and Code Frequency observables."

Kalman filtering and more advanced techniques derived from the Kalman filter, may provide immediate benefits to performance under adverse conditions. Moreover, a transition from traditional architectures to Kalman architectures will open the gateway to potential reimaginings of the GNSS signal processing pathway which may allow GNSS receivers to operate under extremely harsh conditions, and may greatly expand upon the places in which GNSS technology may be used.

Challenges involved in Kalman Filtering

Although Kalman filtering may be a powerful tool with the right planning, it is a technique which involves working through several caveats and critical challenges in order to be properly integrated into a given application. One of the most important of these practical challenges is the assumption of a well defined model of the process and measurement noise covariances, but this information is not directly available a priori and may change as channel conditions change during system operation. Traditionally, these values are estimated a priori through time series analysis, and are implemented within the Kalman filter as constant values Q and R. In general, however, there is no requirement within the Kalman filter that they must be constant, and may instead take time-varying values Qt-1 and Rt. Several methods have been proposed to estimate and update these noise characteristics iteratively as the filter operates. These methods, however, rely on frequentist methods which do not guarantee validity of the noise covariance matrix. Recently, however a purely Bayesian technique has been developed which shows promise to outperform past techniques while guaranteeing that the estimated characteristic is a valid covariance matrix. The primary intention of this research is to test this Bayesian technique within a modern GNSS receiver using real-world data.

To this end, we are collaborating with the GNSS-SDR open source project. The GNSS-SDR open-source project is a GNSS SDR receiver intended to be a testbed for GNSS receiver designers and researchers, supporting them in the complex task of developing and testing new algorithms for the demanding multiband, multisystem GNSS receivers. Currently, the open source GNSS-SDR features an experimental branch which contains an operational tracking block built around a Kalman filter for tracking with GPS L1-C/A satellite signals. Our goal with this collaboration was to expand upon this experimental tracking block to implement an improved Kalman filter which features a second order model as well as Bayesian covariance estimation.

With this goal in mind, we proposed the following three contributions to the GNSS-SDR project as part of this year's Google Summer of Code:

  • Extension of the experimental Kalman filter tracking block previously implemented in GNSS-SDR from a two-state Kalman filter to a three-state Kalman filter (i.e. first order to second order system model)
  • Preliminary assessment of the proposed Bayesian method for estimating noise statistics from a white, stationary discrete time stochastic process performed using MATLAB
  • Implementation of a support library for incorporating the proposed Bayesian method into the experimental Kalman filter tracking block within the GNSS-SDR framework

The tracking block which we will be referencing to throughout this document, and which will be the subject of most of the changes referenced as part of this Google Summer of Code contirbution may be found within the GNSS-SDR repository here.

EDIT 2018-08-21: The working branch 'kf' including all of the changes for this project was merged into the development branch 'next' through commit a833abb.

Project development

The following sections contain a brief discription of the work done and findings for each of these contributions. A PDF document containing a more comprehensive description of the theorretical aspects of this work, including the theory of the Bayesian covariance estimator can be found here.

Extending the 2-State Kalman Filter to 3-States

One of the advantages in using a system architecture based on the Kalman filter is that it is relatively trivial to update the architecture to reflect changes in the signal model. In particular, it may be of interest in certain applications to extend a lower order state space model to a higher order in order to perform a better estimate at each time instance, or to make inferences about behavior. In the case of GNSS tracking, it is sufficient to perform 2-state estimation (code delay and doppler frequency) for the purposes of obtaining a position fix, however, a more robust estimate of these parameters may be obtained by extending to a second order model. Furthermore, by extending to a third order model, it may be possible to make inferences about atmospheric scintillation effects which may be used to mitigate interference as part of future work.

"Block diagram of signal tracking loop with direct-state KF for overall signal tracking loop."

J-H. Won, T. Pany, and B. Eissfeller. Characteristics of Kalman filters for GNSS tracking loops. IEEE Trans. on Aerospace and Electronic Systems, 48(4):3671–3681, Oct. 2012

Implemetation of the 3-state standard Kalman filter within the GNSS-SDR framework was made easier by the existance of a previously implemented 2-state standard Kalman filter added during a past contribution. As a result, the following steps were taken to implement a configurable Kalman filter:

  • Add a configuration parameter allowing the user to select the order of the Kalman filter
  • Add a synchronization parameter between the acquisition stage and tracking stage which communicates the Doppler bin width used for a particular acquisition
  • Update the tracking block initialization to properly configure the Kalman filter parameters for 2-state or 3-state operation
  • Update the tracking block output to record the third state (Doppler rate) in the dump file. Note that this value is always recorded as a floating point number, which takes values of zero during 2-state operation

With the exception of the final step, all of the steps to produce a 3-state configuration for the GPS L1 C/A Kalman filter tracking block were accomplished in commits 0dd99e3 and d565d65 in the GNSS-SDR repository. Further changes to the 3-state Kalman filter parameters required to obtain a position fix, as well as minor bug fixes were made in 032e73e.

3-State Kalman Filter Test

Testing of the 3-state Kalman filter was performed by way of a recorded data test. Data recorded using a hardware receiver at CTTC in 2013 was played back through the GNSS-SDR datapath using the 3-state Kalman filter configuration settings. A position fix was obtained, and the Kalman filter outputs were recorded to a datafile and plotted using MATLAB for visual analysis. The Kalman filter outputs for the 2-stage configuration were also recorded and plotted for comparison. The following plots show those outputs for two satellites, one with a high carrier-to-noise ratio and one with a low carrier-to-noise ratio.

"Comparison of Kalman filter outputs for 2-state Kalman filter tracking to 3-state Kalman filter tracking for GPS L1 C/A signal with PRN 11"

"Comparison of Kalman filter outputs for 2-state Kalman filter tracking to 3-state Kalman filter tracking for GPS L1 C/A signal with PRN 27"

In the data used for this comparison, PRN 11 was received with a high signal to noise ratio. Note that the Carrier Phase and Doppler Frequency estimates between the two configurations are similar under these conditions, but that the 3-state Kalman filter provides a relatively stable estimate of Doppler frequency rate, which may be a useful measurement for certain GNSS applications.

On the other hand, in the data used for this comparison, PRN 27 was received with a low signal-to-noise ratio, leading to signal several signal dropouts in the 2-state estimation. By contrast, for the same signal from the same data set, the 3-state Kalman filter was able to mitigate some of the effects of high noise (low carrier-to-noise ratio) conditions by relying on the internal model of Doppler frequency propagation.

Experimental Bayesian Covariance Estimation

As previously described, one of the biggest challenges in working with Kalman filtering techniques is in modelling the noise characteristsics of the sytem. Namely, optimal filter preformance is contingent on an accurate knowledge of the time-varying statespace covariance matrices Qt-1 and Rt. To this end, we have developed a Bayesian technique based on the concept of conjugate priors, which allows us to make iterative real-time estimates of Rt within the framework of the standard linear Kalman filter. For brevity, the details of this technique will not be described here. For an introduction to the Bayesian covariance estimation technique, please refer to the complete report on our Google Summer of Code findings here, and for a more complete explanation please check out our presentation at the ION GNSS+ 2018 conference in September 2018.

Preliminary Analysis of Bayesian Covariance Estimation using Modelled GPS L1 C/A signals in Matlab

Preliminary assessment of the Bayesian Covariance Estimation method was performed in MATLAB using simulated GPS L1 signals rather than within the GNSS-SDR C++ framework. It was determined that prototyping and experimentation should be performed in MATLAB with data exhibiting completely known behavior before moving onto implementation within the GNSS-SDR framework and testing using real-world data.

Implemetation of Bayesian covariance estimation within MATLAB was performed by creating a function which performs a series of computations given a sample taken from a locally stationary white noise stochastic process, and parameters characterizing the prior knowledge about the distribution of the covariance of the process in the form of a normal-inverse-Wishart distribution. In this implementation, the four parameters characterizing the prior and posterior distributions are maintained in the calling function, a workflow which is well suited to MATLAB's multiple-input-multiple-output function definition scheme.

Testing of this method was performed according to the following steps:

  • Compute initial signal parameter states according to assumed distributions for the possible initial values of the GPS L1 C/A signal properties (i.e. carrier phase, Doppler shift and Doppler rate)
  • Perform monte carlo iterations using noise generated according to chosen values of Qt-1 and Rt to produce states and phase error observations for a simulated GPS L1 signal
  • Iteratively perform linear Kalman filtering according to the standard Kalman filtering equations over each element of the simulated measurement sequence
  • Simultaneously, iteratively perform Bayesian covariance estimation using the innovations computed during the update step of the Kalman filter, keeping track of priori and posteriori normal-inverse-Wishart parameters
  • Once the Kalman filter has converged, update the innovation covariance term St using the mean of the posteriori distribution as described in place of the standard innovation covariance estimate within the Kalman filter

The code used to perform this prototoyping and preliminary assessment were not included in the GNSS-SDR repository, and may instead be found in a personal repository which may be found here. Of particular interest is the function BayesianCovEst.m, which is the prototype Bayesian covariance estimation method, and the script KF_carrier_phase.m, which performs experimental analysis of the method, including: Monte Carlo simulation of GPS L1 C/A signal propagation, Kalman filtering, maintenance of the prior and posterior distributions for the Bayesian covariance estimation function, and logging and plotting of the Kalman and Bayesian filter outputs. This script also uses the Myers_cov.m function, which implements a popular method for iterative covariance estimation based on sample covariance averaging to provide a benchmark for the new Bayesian method.

Bayesian Covariance Estimation MATLAB Test

We tested the proposed methodology in an experiment in which the proposed Kalman filter with Bayesian covariance estimation was applied to the same simulated linear multivariate GNSS signal model previously described under two different conditions: during static conditions where the signal-to-noise ratio maintained a constant value of 30 dB-Hz throughout the entire 10 second duration of the test, and during a simulated fading event in which the signal to noise ratio drops to 19 dB-Hz at the 6th second of the test.

"Variance estimation results for the linear/Gaussian multivariate system (N=1). Comparison of Myers approach and the proposed methodology."

"Root mean squared error of variance estimation for the linear/Gaussian multivariate system for stationary process noise (N=1000). Comparison of Myers approach and the proposed methodology."

The first of the above figures shows an example of the covariance estimates computed using the proposed method, as well as the results using a popular frequentist method for performing covariance estimation known as the Myers method, for one of the trials of the stationary experiment. The second figure shows the root mean squared (RMS) covariance estimation error over one-thousand trials. Upon converging, the Bayesian estimate varies significantly less than the Myers method, and on average maintains a more accurate estimate of the covariance.

"Signal parameters and state estimation results for the linear/Gaussian multivariate system for time varying process noise with strong fading at six seconds (N=1). Comparison of standard Kalman filtering approach and the proposed methodology."

"Variance estimation results for the linear/Gaussian multivariate system or time varying process noise with strong fading at six seconds (N=1). Comparison of standard Kalman filtering approach to the true value of the process noise variance and the static value of the standard Kalman filter."

The first of these two figures shows the behavior of the state evolution, as well as Kalman state estimation with and without Bayesian covariance estimation for a single trial of the fading experiment. Filter behavior both with and without covariance estimation is expectedly identical to the behavior from the first experiment prior to the fading event. Once the signal becomes degraded however, the Kalman filter with Bayesian estimation is able to adapt to the change in noise characteristics and maintain a lower overall estimation error for the remainder of the experiment. This improvement in performance is explained in the second figure, which shows that while the Kalman filter with Bayesian estimation gradually adapts towards the new covariance, the standard Kalman filter continues to operate with an underestimate of the true covariance leading to a loss of optimality.

Implementation of Bayesian Covariance Estimation within GNSS-SDR

Implementation within the GNSS framework presented two important challenges: operation of the Bayesian estimation needed to be optional and configurable, and results on the performance needed to be demonstrable. Unlike the MATLAB implementation described earlier, which took the form of a standalone function which operated on externally maintained variables containing the distribution parameters, Bayesian estimation within the GNSS framework was implemented as an external library containing a class which maintains its own internal state. Performing Bayesian estimation with this class mimics other filter implementations: the class is initialized with a priori values when tracking begins, and updates its internal parameters and estimates based on data input through an update function. This library can be found within the GNSS-SDR repository here.

The first iteration of this library was implemented in 486ac19 with fixes and integration work made in e42467a and MATLAB plotting functionality added in 0fd98b0. Unfortunately, as of August 14, 2018 an unresolved runtime error (segmentation fault) triggered by integration of the bayesian_estimation class within the gnuradio block prevents use (including testing) of this library for processing recorded GPS L1 C/A data. Instead, 05aea34 adds a GTest unit test to the GNSS-SDR unit test library which checks that the Bayesian_estimator class produces a valid output for a large number of iterations with random input.

"Debug results for Bayesian Estimation Positivity Test"

We anticipate that the segmentation fault preventing testing of the Bayesian estimation library with real world data may be resolved within the coming weeks by taking a look at the multithreading protections present in the GNSS-SDR block architecture. We are confident based on our preliminary results that the Bayesian_estimator module will perform as expected.

Conclusions and Future Work

Over the course of Google Summer of Code 2018 we have made several contributions to GNSS-SDR open-source GNSS receiver project. As the report shows, however, the ultimate goals of this project, testing the Bayesian covariance estimation method with real-world GNSS data, was not fully accomplished. Nevertheless, significant progress has been made toward that goal through experimental design, modelling and implementation, and it is our firm belief that with a bit more effort, completion of that final goal is well within reach.

In summary, the contributions we have made this summer are expansion of the experimental Kalman filter tracking block to include a configurable 3-state Kalman filter which operates according to a second order system model. Following this, we performed preliminary assessment of a proposed Bayesian method for estimating noise statistics from a white, stationary discrete time stochastic process, which promises to increase the versatility of the Kalman filtering approach during suboptimal noise conditions. Finally, we implemented and integrated a support library for performing the proposed Bayesian method into the GNSS-SDR framework.

The immediate future work is self evident: issues involved with integrating the Bayesian estimation library into the GPS L1 C/A tracking gnuradio block must be resolved. Following this however, the primary work should be to expand upon the Kalman filter framework to a joint code/carrier framework which replaces both DLL in addition to the PLL from the traditional locked-loop framework. In addition to this, research should be done into using nonlinear extensions of the Kalman filter to perform Kalman filter using pilot signals rather than the arctangent discriminator. Finally, these concepts may be expanded to include other signals beyond GPS L1 C/A, such as those associated with the European Galileo GNSS system, the Russian GLONASS system and the Chinese BeiDou system.