Robust Gaussian sum filtering with unknown noise statistics: Application to target tracking

In many real-life Bayesian estimation problems, it is appropriate to consider non-Gaussian noise distributions to model the existence of outliers, impulsive behaviors or heavy-tailed physical phenomena in the measurements. Moreover, the complete knowledge of the system dynamics uses to be limited, as well as for the process and measurement noise statistics. In this paper, we propose an adaptive recursive Gaussian sum filter that addresses the adaptive Bayesian filtering problem, tackling efficiently nonlinear behaviors while being robust to the weak knowledge of the system. The new method is based on the relationship between the measurement noise parameters and the innovations sequence, used to recursively infer the Gaussian mixture model noise parameters. Numerical results exhibit enhanced robustness against both non-Gaussian noise and unknown parameters. Simulation results are provided to show that good performance can be attained when compared to the standard known statistics case.


INTRODUCTION
The problem under study concerns the derivation of efficient and robust methods to solve the recursive Bayesian filtering problem for nonlinear systems corrupted by non-Gaussian noise with unknown statistics.The classic filtering problem involves the recursive estimation of time-varying unknown states of a system using the incoming flow of information along some prior statistical knowledge about the variations of such states.The general discrete state-space model can be expressed as where x k ∈ R nx is the hidden state of the system at time k, f k−1 (•) is a known, possibly nonlinear, function of the states; and v k is referred to as process noise; z k ∈ R nz is the measurement at time k, h k (•) is a known, possibly nonlinear, function, which relates measurements with states; and n k is referred to as measurement noise, independent of v k .The optimal Bayesian filtering solution is given by the marginal distribution p(x k |z 1:k ), which gathers all the information about the system contained in the available observations, with z 1:k = {z1, . . ., z k }.
The Kalman filter (KF) provides the closed form solution to the optimal filtering problem in linear/Gaussian systems, assumptions that not always hold, reason why suboptimal techniques have to be used.A plethora of alternatives have been proposed to solve the nonlinear estimation problem, among them, the Extended KF (EKF) and the family of Sigma-Point KFs (SPKF) [1] within the Gaussian framework, and the family of Sequential Monte Carlo (SMC) methods [2], and Gaussian Sum Filters (GSFs) [1], for arbitrary noise distributions.The latter constitutes an appealing alternative to SMC methods, which provide a powerful framework to deal with nonlinear/non-Gaussian systems at expense of a high computational load, being difficult to embed in digital light processors or in realtime applications.But the main limitation of all these methods is that they assume some a priori knowledge of the noise statistics affecting the system (i.e., not only its distribution but its parameters).In many real-life systems, Gaussian noise models do not apply and the noise statistics are unknown.In these scenarios the methods based on the standard Gaussian Kalman framework give poor performance or even diverge, and we cannot directly apply SMC methods and GSFs because we need to estimate the states together with the noise statistics.Dropping these two classical assumptions, Gaussianity and known statistics, is the starting point of this work.
In this contribution, our attention focuses on the adaptive/robust filtering problem in this context, to deal with state-space models corrupted by non-Gaussian measurement noise with unknown statistics.Several improvements to the standard GSF and related methods have been proposed in the literature [3,4,5,6], showing the promising applicability of GMMs, but all of them still considering known statistics.To the authors' knowledge, an adaptive filtering solution (e.g., adaptive filtering referring to the unknown statistics estimation and robustness concepts) within the GSF framework to deal with systems corrupted by GMM non-Gaussian noise with unknown statistics is still an open problem.
The noise statistics estimation problem has been mainly studied for linear/Gaussian systems [7], and its extension to nonlinear systems usually makes use of linearization techniques.In order to deal with more general models, the problem has also been considered into the system parameters estimation framework using SMC approaches [8,9,10].In the case of GSFs, the noise statistics estimation reduces to the estimation of the GMM parameters.To solve this problem and provide a robust filtering solution, this paper proposes the Adaptive Gaussian Sum Filter (AGSF), which adjusts the parameters of the GMM distributions adaptively with no prior statistical information.

GAUSSIAN SUM FILTERING
In many applications, the Gaussian assumption does not hold [11,12] and non-Gaussian distributions have to be used in order to capture the statistical behavior of the system dynamics and/or the prior knowledge.On the basis of the Wiener approximation theory, any probability density function (pdf) can be expressed, or approximated with a given level of accuracy, using a finite sum of Gaussian densities, what is known as a Gaussian mixture model (GMM): px(x) = L i=1 αiN (x; mi, Bi), where L i=1 αi = 1, αi ≥ 0 for i = 1, . . ., L, and N (x; mi, Bi) represents a Gaussian distribution with mean mi and covariance matrix Bi (i.e., for the sake of simplicity we use the notation p(x) ∼ {αi, mi, Bi} L i=1 ).Therefore, any non-Gaussian distribution of interest can be expressed with a GMM, and considering GMM distributed additive noises, the predictive and posterior densities can also be expressed using Gaussian sums by using the following equalities.Assume that the posterior distribution, p(x k |z 1:k ), can be represented with a GMM, ), then the mean and covariance matrix can be easily computed as, To obtain the recursive Bayesian filter in the GMM framework, the initial state (prior knowledge), and both process and measurement noise distributions, are assumed to follow a GMM distribution: Let us also assume that the posterior distribution at time k − 1 is given by the following GMM, The predictive density p(x k |z 1:k−1 ) is obtained using the transition density and the Chapman-Kolmogorov equation, and approximated within the KF framework using a bank of n k = s k l k−1 time update steps of the EKF, where for each Gaussian term the state prediction and the corresponding error covariance matrix are with i=1 and the posterior distribution, p(x k |z 1:k ), as done before with the predictive distribution, can be approximated using a GMM and a Kalman approach with a bank of r k = m k s k l k−1 measurement update steps of the EKF as For each term the state estimation and the corresponding error co- . Refer to Algorithm 2 for the computation of the Kalman gain K ijl,k , the predicted measurement ẑijl,k|k−1 , and the innovations covariance matrix Σ zz,ijl,k|k−1 .The weighting factor can be written as with i=1 α ijl,k|k xijl,k|k .Notice that within one time step k, the number of Gaussian terms increased from l k−1 to m k s k l k−1 .In practice, Gaussian mixture reduction techniques [13,3] have to be used to reduce the number of Gaussian terms after each iteration (in this work the implementation reported in [14] was used).

ADAPTIVE RECURSIVE BAYESIAN FILTER
This section proposes a new AGSF to solve the robust adaptive Bayesian filtering problem for nonlinear discrete state-space models corrupted by non-Gaussian measurement noise.This method estimates the statistical parameters of the GMM measurement noise using the innovations provided by the GSF.

Innovations and Measurement Noise
The innovations are denoted as zk and the predicted state estimates are obtained as xk|k−1 = E(x k |z 1:k−1 ), which are available from the Time Update step of the GSF.The innovations sequence is de- . If the measurement function h k is linear, the above approximation will turn into equality.Therefore, the pdf of innovations is equal to the convolution of the one-step prediction error's distribution and the measurement noise distribution, where ⊗ stands for the convolution operator.It is easy to see that both p(H k xk|k−1 |z 1:k−1 ) and p(n k ) are Gaussian mixture distributed.Therefore, the essence of estimating p(n k ) from zk is a deconvolution problem.Using (4), the pdf of xk|k−1 follows a GMM and is given by, where To simplify, we assume that the prediction error, xk|k−1 is Gaussian distributed, which means that we are only interested in the first and second moments of the distribution: Then, we conclude that Using ( 10) and (8), Therefore, the GMM innovations distribution has the same number of Gaussian terms than the measurement noise distribution.

Adaptive Estimation of the Noise Statistics
The adaptive filtering solution essentially refers to the estimation of the GMM parameters of the measurement noise based on the innovations sequence.The proposed method uses an Expectation Maximization (EM) solution [15], which is applied to process the 3: Estimate the m k means: ηl,k = ηl,k−1 + G l,k (z k − ηl,k−1 ).4: Estimate the m k covariance matrices: innovations, which obey a Gaussian mixture distribution, and finally get the GMM parameters estimates.Assume that K innovations are collected from time instant 1 to K, and stacked in ZK = (z1, z2, . . ., zK ).The Gaussian mixture parameters for the l th Gaussian term at time instant k are θ l,k = (γ l,k , η l,k , S l,k ), and the full set is ).In this case, the parameters of each Gaussian term and the number of terms are assumed to be constant within the innovations sequence, so the time index is dropped to simplify.The parameters estimated in the previous iteration are denoted as θ g l = {γ g l , η g l , S g l } m l=1 .From the likelihood expression, an EM method is applied to find the ML estimates of the GMM measurement noise parameters [15].Two steps are performed recursively until convergence: l ) for l = 1, . . .m, and k = 1, . . ., K.
The final estimates are obtained as γ l,k = γl , η l,k = ηl , and Σ n,l,k = Ŝl − H k Σ k|k−1 H T k , for k = 1, . . ., K.This is a good solution if data can be stored and processed off-line, but in standard filtering applications the measurement noise parameters need to be estimated online.The sequential counterpart is presented hereafter.First, the auxiliar weights are obtained from the previous estimates and the current innovation, , and then the new estimates are obtained as: One iteration (at time instant k) of the noise statistics estimation method is sketched in Algorithm 1 (i.e., the previous estimates and the current innovation zk are available), and marked in red in the complete AGSF sketched in Algorithm 2.

COMPUTER SIMULATIONS
In order to provide illustrative numerical results, the performance of the proposed method was evaluated in a 2-D radar target tracking application.The measurements were range and azimuth, z k = [r k ψ k ] T , and the states to be tracked were position and velocity of the target, respectively gathered in vector 1: Initialization: xi,0 ∼ N (µi,0, Σx,i,0) for i = 1, . . ., l0. 2: for k = 1 to ∞ do

4:
for j = 1 to s k do 5: Estimate the predicted states:

7:
Estimate the error covariance matrices:

8:
end for

14:
for l = 1 to m k do 15: for j = 1 to s k do 16: for i = 1 to l k do 17: Estimate the predicted measurements: ẑijl,k|k−1 = h k (x ij,k|k−1 ) + η l,k .

21:
Estimate the updated states:

22:
Estimate the corresponding error covariances:

23:
end for

24:
end for

25:
end for
28: end for [p x,k , v x,k , p y,k , v y,k ] T .A two-terms GMM measurement noise was considered to model possible outliers in the measurement and the clutter effect in radar tracking systems.Notice that the standard GSF with known statistics was used as the benchmark, being the ultimate performance achievable for the proposed method.The dynamics of the target were the sampling period and w k ∼ N (0, 0.01 • I2) a Gaussian process noise.The covariance of the resulting process noise can be expressed as The measurements were modeled as where n k followed a GMM distribution,  initialized as γi,0 = γi + 0.05δγi, where δγi ∼ U [0, 1], ηi,0 = ηi + δηi • [0.1 0.01] T , where δηi ∼ N (0, 1), Si,0 = 2Si, for i = 1, 2.
Figure 1 (top) shows the filtering results of the AGSF when m = 1, 2, 3, and the GSF with known statistics.In this case, when considering 2 or 3 terms for the measurement noise approximation, good performances were obtained.The best results were obtained when the number of terms was exactly the correct one (2 in this case), but with a higher number of terms the loss in performance was not significant.In these cases the AGSF estimates the states of the system and the measurement noise statistics accurately, with limited performance degradation on the state estimation performance when compared to the benchmark, i.e., the GSF with known statistics.
It seems obvious that the filtering results considering only one term were the worst, showing the improvement and therefore the interest of using the proposed method rather than an Adaptive KF when dealing with non-Gaussian measurement noise.Figure 1 (bottom) depicts the approximation of the true measurement noise pdf for three different approximations, m = 1, 2, 3 terms, verifying the previous performance results because the approximation of the noises' pdf with 2 and 3 terms are almost the same, and much better than the approximation using only 1 Gaussian term.

CONCLUSIONS
This paper presented a solution to the robust Bayesian filtering problem.The new Adaptive Gaussian Sum Filter, based on an Expectation-Maximization type solution, is able to deal with non-Gaussian measurement noise distributions with unknown statistics.This adaptive scheme can get rid of the limitations on the knowledge of the measurement noise, which means that the filter estimates the parameters of the Gaussian mixture model measurement noise exploiting information from the innovations.The method was validated by computer simulations, showing that the proposed algorithm can efficiently cope with the state estimation problem while correctly dealing with the unknown statistics of the non-Gaussian measurement noise.Future work goes towards formally proving the convergence of the algorithm.