Variational Bayesian Based IMM Robust GPS Navigation Filter

Dah-Jing Jwoand Wei-Yeh Chang

1Department of Communications, Navigation and Control Engineering, National Taiwan Ocean University Keelung,202301, Taiwan

2TDK Taiwan Corp., Yangmei, Taoyuan, 326021, Taiwan

Abstract: This paper investigates the navigational performance of Global Positioning System (GPS) using the variational Bayesian (VB) based robust filter with interacting multiple model (IMM) adaptation as the navigation processor.The performance of the state estimation for GPS navigation processing using the family of Kalman filter (KF) may be degraded due to the fact that in practical situations the statistics of measurement noise might change.In the proposed algorithm, the adaptivity is achieved by estimating the timevarying noise covariance matrices based onVB learning using the probabilistic approach, where in each update step, both the system state and time-varying measurement noise were recognized as random variables to be estimated.The estimation is iterated recursively at each time to approximate the real joint posterior distribution of state using the VB learning.One of the two major classical adaptive Kalman filter (AKF) approaches that have been proposed for tuning the noise covariance matrices is the multiple model adaptive estimate (MMAE).The IMM algorithm uses two or more filters to process in parallel, where each filter corresponds to a different dynamic or measurement model.The robust Huber’s M-estimation-based extended Kalman filter (HEKF) algorithm integrates both merits of the Huber M-estimation methodology and EKF.The robustness is enhanced by modifying the filter update based on Huber’s M-estimation method in the filtering framework.The proposed algorithm, referred to as the interactive multi-model based variational Bayesian HEKF (IMM-VBHEKF), provides an effective way for effectively handling the errors with time-varying and outlying property of non-Gaussian interference errors, such as the multipath effect.Illustrative examples are given to demonstrate the navigation performance enhancement in terms of adaptivity and robustness at the expense of acceptable additional execution time.

Keywords: GPS; variational bayesian; Huber’s M-estimation; interacting multiple model; adaptive; outlier; multipath

1 Introduction

Non-Gaussian noise is often encountered in many practical environments and the estimation performance deteriorates dramatically.Known to be one of the dominant error sources in high accuracy global navigation satellite systems (GNSS) [1] positioning systems, such as the Global Positioning System (GPS) [2,3], multipath effects [2,4] occur when GPS signals arrive at a receiver site via multiple paths due to reflections from nearby objects, such as the ground and water surfaces,buildings, vehicles, hills, trees, etc.Many estimation algorithms have been studied to eliminate the positioning error caused by multipath.Since multipath errors are among uncorrelated errors that are not cancelled out during observation differencing, the performance of high precision GPS receivers are mostly limited by the multipath induced errors.One of the most important issues in GPS system performance improvement is the interference suppression techniques.

The well-known Kalman filter (KF) [5] provides optimal estimate, namely minimum mean square error (MSE), of the system state vector and has been recognized as one of the most powerful state estimation techniques.As a recursive minimum norm estimation technique, the KF exhibits sensitivity to deviations from the assumed underlying probability distributions.Due to its simple structure, stable performance and low computational complexity, the conventional adaptive filtering algorithm where the least MSE is involved has been widely used in a variety of applications in the fields of adaptive signal processing and machine learning.However, the MSE criterion is limited to the assumption of linearity and Gaussianity while most of the noise in real word is non-Gaussian.The performance deteriorates significantly in the non-Gaussian noise environment.The traditional Kalman-type filter provides the best filter estimate when the noise is Gaussian, but most noise in real life is unknown,uncertain and non-Gaussian.The extended Kalman filter (EKF) is a nonlinear version of the KF and has been widely employed as the GPS navigation processor.The fact that EKF highly depends on a predefined dynamics model forms a major drawback.Furthermore, it is not optimal when the system is disturbed by heavy-tailed (or impulsive) non-Gaussian noises.To solve the performance degradation problem with non-Gaussian errors or heavy-tailed non-Gaussian noises, some robust Kalman filters have been developed by using non-minimum MSE criterion as the optimality criterion.

Essentially based on modification of the quadratic cost function in the filter framework, the Huber M-estimation methodology [6-9] provides between smooth norm properties for small residuals and robust norm properties for large residuals.The technique that relies on Huber’s generalized maximum likelihood (ML) methodology exhibits robustness against deviations from the commonly assumed Gaussian probability density functions and can solve the non-Gaussian distribution problem efficiently and has been successfully employed for robust state estimation, inertial navigation system and visual tracking applications.The robust Huber’s M-estimation-based EKF (HEKF) algorithm integrates both merits of the Huber M-estimation methodology and EKF.For the signals contaminated with non-Gaussian noises or outliers, a robust scheme combining the Huber M-estimation methodology and the EKF framework is beneficial where the Huber M-estimation methodology is used to reformulate the measurement information of the EKF to handle outliers and provide robustness against deviations from Gaussianity.The measurement update can be modified to enhance robustness using the Huber M-estimation methodology.

The variational Bayesian (VB) approach [10-15] is an inference method that utilizes a simple distribution to approximate the true posterior distribution of hidden variables, usually assuming that the hidden variables are independent of each other.The VB has been developed for a wide range of models to perform approximate posterior inference at low computational cost in comparison with the sampling methods.These methods assume a simpler, analytically tractable form for the posterior.The purpose of the variation is to find a variational distribution of the posterior probability density function (pdf) to minimize the Kullback-Leibler divergence (KLD) between the true posterior pdf.Two main approaches are either to derive a factored free form distribution or to assume a fixed form posterior distribution (e.g., a multivariate Gaussian, with possibly a suitable parametrization of the model).

The classical way of solving the problem of uncertain parameters is to use adaptive filters [16]where the model parameters or the noise statistics are estimated together with the dynamic states.The classical noise adaptive filtering approaches can be divided into Bayesian, maximum likelihood,correlation and covariance matching methods.The Bayesian approach is the most general of these and the other approaches can often be interpreted as approximations to the Bayesian approach.Examples of Bayesian approaches to noise adaptive filtering include state augmentation based methods,multiple model methods such as the interacting multiple model (IMM) algorithm [17-21].The IMM algorithm has the configuration that runs in parallel several model-matched state estimation filters, which exchange information (namely interact) at each sampling time.The IMM algorithm has been originally applied to target tracking, and recently extended to navigation application.A model probability evaluator calculates the current probability of the vehicle being in each of the possible modes.A global estimate of the vehicle’s state is computed using the latest mode probabilities.The algorithm carries out a soft-switching between the various modes by adjusting the probabilities of each mode, which are used as weightings in the combined global state estimate.The covariance matrix associated with this combined estimate takes into account the covariances of the mode-conditioned estimates as well as the differences between these estimates.

The algorithm proposed in this paper intends to provide an effective way for effectively handling the errors with time-varying and outlying property of non-Gaussian interference errors.The method utilizes the VB learning to approximate the noise strength for time-varying noise covariances, the Huber M-estimation methodology to enhance robustness especially for overcoming the problem of contamination distribution or outliers, and the IMM algorithm to furtherly tune the noise covariance matrices.The remaining of this paper is organized as follows.In Section 2, preliminary background on the Huber’s M-estimation-based extended Kalman filter is reviewed.The variational Bayesian algorithm is discussed in Section3.In Section4,the interacting multiple model algorithm is presented.In Section 5, simulation experiments are carried out to evaluate the performance and effectiveness.Conclusions are given in Section 6.

2 Huber’s M-Estimation-Based Extended Kalman Filter

For the nonlinear single model equations in discrete-time form, we have

where the state vector xk∈Rn, process noise vector wk∈Rm, measurement vector zk∈Rm,measurement noise vector vk∈Rm,Qkis the process noise covariance matrix and Rkis the measurement noise covariance matrix.The vectors wkand vkin Eqs.(1) and (2) are zero mean Gaussian white sequences has zero cross-correlation with each other:

whereE[·]represents expectation, and superscript “ T”denotes matrix transpose.The symbol δikstands for the Dirac function:

Based on robust strategy, the Huber M-estimation methodology possesses the ability to solve pollution distribution or outliers by improving filter update.The essence of Huber’s robust estimation method is to solve the filter estimation value when the minimum cost function is obtained.For further applications,based on the generalized maximum likelihood estimation, amore general object function is given by Huber

where ρ(·) is an appropriate positive symmetric function,υk,iis the residual.The cost function given by Huber is

where γ is the default parameter.Amounts of experiences show that the performance is best when γ= 1.345.For Huber’s filtering robustness framework, reweighted average is used by

where υk,iis the residual,φ(υk,i) =ρ′(υk,i) is the influence function.

The Huber M-estimation methodology makes use of the weighting matrix to recast the measurement information.The weighting matrix can be interpreted in terms of reweighting the residual error covariance and constructing the measurement process.The modified covariance has the form

The HEKF algorithm is summarized as follow:

Stage 1: Prediction steps/time update equations

(1) Predict the state vector

(2) Predict the state error covariance matrix

Stage 2: Correction steps/measurement update equations

(3) Compute the residual

(4) Compute the reweighting covariance matrix for measurement noise

(5) Compute Kalman gain matrix

(6) Update state vector

(7) Update state error covariance matrix

3 The Variational Bayesian Approach

The variational Bayesian approach is an inference method that utilizes a simple distribution to approximate the true posterior distribution of hidden variables to find a variational distribution of the posterior probability density function (pdf) to minimize the KLD between the true posterior pdf.The variable decibel parameter is estimated by iterative calculations in which the gradient descent algorithm is involved.The difference between the lower bounds of two adjacent repeated calculations is used as the basis for judging convergence.

3.1 Variational Bayesian Extended Kalman Filter

Based on the Bayesian rule, the posterior distribution of the system state and the measurement noise covariance can be represented as

which can be simplified top(zk|xk,Rk).For the first-order Markov process

the Chapman-Kolmogorov equation can be obtained through the integration.

Given the next measurement zk, the predictive distribution is updated to a posterior distribution using the Bayes rule

The conditional probability distribution can be approximated as a product of Gaussian and Inverse-Gamma distributions

whereN(ˆk|k-1,Pk|k-1) is the Gaussian distribution for prediction xkwith parametersˆk|k-1and Pk|k-1.αk|k-1,i,βk|k-1,iare two parameters of the Inverse-Gamma distribution for, thei-th diagonal element of Rk.

The variational Bayesian extended Kalman filter (VBEKF) employes a heuristic approach in the calculation process and predicts the distribution parameter by the first-order approximation

The parameters of the new distribution is updated through

and the measurement noise variance can be estimated

The VBEKF algorithm is summarized as follows:

Stage 1: Prediction steps/time update equations

(1) Predict the state vector:

(2) Predict the state error covariance matrix

(3) Predict the distribution parameter

Stage 2: Correction steps/measurement update equations

(4) Compute measurement noise variance

(5) Compute Kalman gain matrix

(6) Update state vector

(7) Update state error covariance matrix

(8) Update the measurement distribution parameter

3.2 Variational Bayesian Huber’s M-Estimation Based Extended Kalman Filter

The positioning accuracy is degraded due to the presence of outliers and time-varying noise strength.Incorporation of the VB and Huber M-estimation methodologies into the EKF yielding the variational Bayesian Huber’s M-estimation based extended Kalman filter (VBHEKF) algorithm can effectively overcome the outliers and time-varying variance in the measurement noise and improve the positioning accuracy.The VBHEKF algorithm is summarized as follows:

Stage 1: Prediction steps/time update equations

(1) Predict the state vector

(2) Predict the state error covariance matrix

(3) Predict the distribution parameters

Stage 2: Correction steps/measurement update equations

(4) Compute measurement noise variance:

(5) Compute the residual

(6) Compute the reweighting covariance matrix for measurement noise

(7) Compute Kalman gain matrix

(8) Update state vector

(9) Update state error covariance matrix

(10) Update the measurement distribution parameter

4 Interacting Multiple Model Extended Kalman Filter

Based on pseudo-Bayesian theory, the IMM algorithm employs two or more filters to process in parallel where each filter corresponds to a different state space or measurement model.It is an adaptive mechanism which dynamically identifies uncertainties or modeling errors can be adopted and is mainly an algorithm involving Markov chain switching coefficients between different models.In each model an EKF is running, and the IMM algorithm makes uses of model probabilities to weight the inputs and output of a bank of parallel filters at each time instant.

Assuming that a target with r states, corresponding to r models, set the system state equation and system observation equation represented by thej-th model are as follows:

where the state vector xk∈Rn, process noise vector wk∈Rm, measurement vector zk∈Rm,measurement noise vector vk∈Rm,Qkis the process noise covariance matrix and Rkis the measurement noise covariance matrix.The vectors wkand vkin Eqs.(35) and (36) are zero mean Gaussian white sequences has zero cross-correlation with each other:

whereE[·]represents expectation, and superscript “ T”denotes matrix transpose.The symbol δikstands for the Dirac function:

The IMM-EKF algorithm is summarized as follows:

Stage 1: Model interaction/mixing

Compute mixed state and corresponding covariances of modelj

Stage 2: Model filtering

Stage 3: Model probabilities update

Compute the likelihood function

wherevj= zk-and.The model probability is updated according to the model likelihood and model transition probability

Stage 4: Combination of model estimates and covariances

The algorithm performs real-time detection by setting model filters corresponding to the number of possible models, setting weight coefficients and model update probabilities for each filter, and finally performing weighted fusion to calculate the current optimal system state estimation to achieve the purpose of model adaptation.

5 Illustrative Examples

To validate the effectiveness of the proposed approaches, simulation experiments have been carried out to evaluate the performance of the proposed approach in comparison with the other conventional methods for GPS navigation processing.The computer codes were developed by the authors using the Matlab®software.The commercial software Satellite Navigation (SatNav) Toolbox by GPSoft LLC [22] was employed for generation of the GPS satellite orbits/positions and thereafter, the satellite pseudoranges, carrier phase measurement, and constellation, required for simulation.

The simulated pseudorange error sources corrupting GPS measurements include ionospheric delay, tropospheric delay, receiver noise and multipath.It is assumed that the differential GPS(DGPS)mode is available and therefore most of the receiver-independent common errors can be corrected,while the multipath and receiver thermal noise cannot be eliminated.The multipath interferences are added into the GPS pseudorange observation data during the vehicle moving.Two scenarios dealing with two types of interferences in pseudorange observables are carried out, which covers pseudorange errors involving (1) outlier type of multipath interferences during the vehicle moving, and (2) timevarying variance in the measurement noise with additional outliers.

Fig.1 shows the test trajectory for the simulated vehicle and the skyplot during the simulation.In the simulation, there are 8 GPS satellites available.Since the research focus on the mitigation of multipath errors, the influence of measurement noise is relatively critical.A vehicle is designed to perform the uniform accelerated motion to reduce the impact caused by unmodeling system dynamic errors.First of all, performance comparison for the four estimator/filter including EKF,VBEKF, HEKF and VBHEKF is shown.Secondly, performance enhancement based on the IMM configuration is presented.

Figure 1: Test trajectory for the simulated vehicle (left) and the skyplot (right) during the simulation

5.1 Scenario 1: Pseudorange Observable Involving Outlier Type of Multipath Errors

In Scenario 1, mitigation of the errors for pseudorange observable involving outlier type of multipath interferences is investigated.There are totally seven time durations where additional randomly generated errors are intentionally injected into the GPS pseudorange observation data during the vehicle moving.Tab.1 provides the information of the standard deviation for the 5 outliers.

Table 1: Information of the standard deviation for the five outliers

5.1.1 Performance Comparison for EKF and its Variants

Comparison of GPS navigation accuracy for the four schemes: EKF, VBEKF, HEKF and VBHEKF is shown in Fig.2.The results show that either the VB or Huber’s algorithms can assist EKF to effectively deal with the outliers in the pesudorange observables individually and combination of the two algorithms can furtherly enhance the performance.

Figure 2: Comparison of positioning accuracy for extended Kalman filter (EKF) and its variants-Scenario 1

5.1.2 Performance Enhancement Based on the Interacting Multiple Model Configuration

Incorporation of the IMM into the various filters is utilized for further performance enhancement.Fig.3 illustrates the positioning accuracy when the IMM configuration is involved for various algorithms: IMM-EKF, IMM-VBEKF, IMM-HEKF and IMM-VBHEKF.For the case that the VB does not possess sufficient capability to resolve the outlier type of interference, the IMM-VBHEKF demonstrates substantial performance improvement in navigation accuracy with acceptable extra computational expense.Tab.2 summarizes the estimation performance and execution time for various algorithms.

Figure 3: Comparison of positioning accuracy when the interacting multiple model (IMM) configuration is involved for Scenario 1

Table 2: Performance comparison for various algorithms-Scenario 1

5.2 Scenario 2: Errors Involving Time-Varying Variance in Measurement Noise with Additional Outliers

Scenario 2 is designed for investigating the performance when dealing with the time-varying Gaussian measurement noise with additional outliers.Description of time varying measurement variances in the five time intervals is shown in Tab.3.The variances of time-varying measurement noiserρiin the five time intervals for this scenario are assumed to follows the variation:rρi= 1→10→1→17→1 (in unit of m2), where the‘arrows (→)’is employed for indicating the time-varying trajectory of variances.The set of unknown time-varying statistical parameters of noise needs to be estimated with the system states and the error covariance.

Table 3: Description of the variation of noise strength and number of additional outliers in the five time intervals for Scenario 2

5.2.1 Performance Comparison for EKF and its Variants

Comparison of GPS navigation accuracy for the EKF, VBEKF, HEKF and VBHEKF is shown in Fig.4 where comparison of positioning accuracies without outlier and with outlier, respectively is provided.The results show that both VB and Huber based EKF can effectively improve the positioning performance.As can be seen, both the VB and Huber algorithms can be adopted into the EKF to improve GPS navigation accuracy in the environment of time-varying Gaussian noise.When additional outliers are involved in addition to the time-varying interferences, neither VB nor Huber only algorithm possesses sufficient capability to resolve the problem for such kind of challenge.The VBHEKF demonstrates substantial performance improvement in navigation accuracy at the expense of acceptable extra computational cost.

Figure 4: Comparison of positioning accuracy for EKF and its variants-Scenario 2 (a) without outlier(b) with outlier

5.2.2 Performance Enhancement Based on the Interacting Multiple Model Configuration

Comparison of positioning accuracies for the four algorithms: IMM-EKF, IMM-VBEKF, IMMHEKF and IMM-VBHEKF is shown in Fig.5.Positioning accuracies for the cases without and with outliers, respectively, are presented.The Huber’s approach does not catch up the variation of noise strength very well.With the assistance of VB, the VBHEKF can further improve the performance.From the other view point, the adaptation capability of noise variance for the VBHEKF has been improved with the assistance of the IMM mechanism.Tabs.4 and 5 summarize the performance for various algorithms without and with outliers, respectively.Among the various approaches, the IMMVBHEKF provides the best positioning accuracy with only limited amount of extra execution time.

Table 4: Performance comparison for Scenario 2-without outliers

Table 5: Performance comparison for Scenario 2-with outliers

Figure 5: Comparison of positioning accuracy for Scenario 2 when the IMM configuration is involved when IMM algorithm is incorporated (a) without outlier (b) with outlier

6 Conclusions

In this paper, an interacting multiple model based variational Bayesian robust extended Kalman filter is adopted for processing GPS navigation states.The method utilizes the variational Bayesian learning to approximate the noise strength and to provide a strong tracking capability for time-varying noise covariances; the Huber M-estimation methodology based on the robust strategy is employed to enhance robustness especially for overcoming the problem of contamination distribution or outliers;the IMM algorithm is introduced for further tuning of the noise covariance matrices.Combination the three merits leads to the estimator referred to as the IMM-VBHEKF.Illustrative examples have been provided for validation of the method.Two scenarios are presented representing two possible peudorange error patterns, including (1) observable involving outlier type of multipath errors;(2) observable involving time-varying variance in measurement noise with additional outliers.The proposed method has demonstrated performance enhancement in terms of adaptivity and robustness at the expense of acceptable additional execution time.

Funding Statement:This work has been partially supported by the Ministry of Science and Technology,Taiwan [Grant Numbers MOST 108-2221-E-019-013 and MOST 109-2221-E-019-010].

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.