Robust Estimation Fusion in Wireless Senor Networks with Outliers and Correlated Noises

This paper addresses the problem of estimation fusion in a distributed wireless sensor network (WSN) under the following conditions: (i) sensor noises are contaminated by outliers or gross errors; (ii) process noise and sensor noises are correlated; (iii) cross-correlation among local estimates is unknown. First, to attack the correlation and outliers, a correlated robust Kalman filtering (coR2KF) scheme with weighted matrices on innovation sequences is introduced as local estimator. It is shown that the proposed coR2KF takes both conventional Kalman filter and robust Kalman filter as a special case. Then, a novel version of our internal ellipsoid approximation fusion (IEAF) is used in the fusion center to handle the unknown cross-correlation of local estimates. The explicit solution to both fusion estimate and corresponding covariance is given. Finally, to demonstrate robustness of the proposed coR2KF and the effectiveness of IEAF strategy, a simulation example of tracking a target moving on noisy circular trajectories is included.


Introduction
Estimation fusion, or data fusion for estimation, has widespread applications in many practical situations that data from multiple sources are involved, for example, guidance, defense, robotics, integrated navigation, target tracking, and GPS positioning [1].Combining the results of multiple sensors can provide more accurate estimation than using a single sensor [2].There are two basic fusion architectures [1,3]: centralized and distributed (referred to as measurement fusion and estimate fusion in target tracking, resp.), depending on whether raw data are sent to the fusion center or not.Both architectures have pros and cons in terms of optimality, channel requirements, requirements, reliability, information sharing, and so forth.For the distributed fusion, it has been realized for many years that local estimates (track) have correlated errors thereafter the original work of [3].How to counter this cross-correlation has been a central topic in distributed fusion.One problem with the Kalman Filtering is that it requires either that the measurements are independent or that the cross-covariance is known [4,5].As is well known, the independent assumption can be relaxed in the case of correlated data, if the cross-covariance information is available.The optimal KF-based approach that the KF maintains cross-covariance information between updates is proposed considering the correlation among the local estimates [3,[6][7][8].
A common simplification is to assume the crosscovariance to be zero; that is, the measurements are independent.However, the simplification may cause that the KF produces nonconservative covariance.This leads to an artificially high confidence value, which can lead to filter divergence [9].Although under the assumption that the cross-correlation is known, the optimal KF-based approach scales quadratically with the number of updates.This makes optimal KF-based approach impractical [10].Several ingenious techniques, such as the tracklets technique [11], covariance intersection [12], and those based on the information graph [13], have been developed.Unfortunately, if the crosscorrelation information is missing or incomplete, the Kalman filter cannot be applied.In such situations, to allow the use of the Kalman filter, the independence of the sources is 2 International Journal of Distributed Sensor Networks often assumed and the correlation is simply ignored in the fusion process, for example, the simple convex combination (SCC) method [14].This makes the filter over optimistic in its estimation, which may lead to divergence [7].Recently proposed covariance intersection (CI) filtering [8,10] is based on convex combination of information matrices, that is, inverse covariance matrices and the corresponding information states.The algorithm provides a general framework for information fusion with incomplete knowledge about the signal sources since it yields consistent estimates for any degree of cross-correlation.Since covariance intersection filtering requires optimization of a nonlinear cost function and instead of underestimation of the actual covariance matrix, the covariance intersection method overestimates it, which obviously results in a significant decrease in performance.To avoid both the inconsistency of the basic convex combination and the lack of performance of the covariance intersection method, internal ellipsoid approximation fusion (IEAF) has been proposed [15].For this approach, the largest volume ellipsoid within the intersection of two ellipsoids can be computed by internal ellipsoid approximation.
On the other hand, the distribution of noise arising in application deviates frequently from assumed Gaussian model, often being characterized by skewed (asymmetric) or heavier tails generating the outlier.A sufficiently far away located outlier can completely cause the least squares estimator or the Kalman filter to break down [16,17].This will degrade the fusion performance greatly.Therefore it is of practical interest to consider filters which are robust to perform fairly well in non-Gaussian environment especially in the presence of outliers, and some results have been obtained during the last decade.Robust statistical procedures provide formal methods to spot the outlying data points and reduce their influence.Most of the contributions in this area have been directed toward censoring data; namely, if an observation differs sufficiently from its predicted value, then it is discarded.For example, an M-estimate filter for robust adaptive filtering in impulse noise is proposed in [18]; a recursive adaptive algorithm and a robust threshold estimation method are derived employing an M-estimate cost function.Djurovic and Kovacevic established the equivalence between the Kalman filter algorithm and a particular least squares regression problem.Based on the equivalence it solved the robust estimation with unknown noise statistics with the help of M-estimate method, and the equivalence between the Kalman filter and proposed technique is established [19].Recently, in [20] a deemphasis weighting approach is used to suppress the effect of outliers in background samples during the formation of a sample covariance matrix.When outliers are present when comparing the results from processing simulated and real coherent radar data using the proposed approach with results using no outlier suppression and censored sample matrix inversion pruning methods, the deemphasis techniques are shown to produce the most robust diction performance.Very recently, by reformulating the traditional Kalman filter into a least square form, a novel version of RKF has been proposed using L 1 -regularized optimization [21].To the best of the authors' knowledge, however, there is little result discussing how to eliminate or reduce the influence of outliers on the fusion performance, which remains a challenging problem so far.
In this paper, the problem of estimation fusion in a distributed architecture under the following conditions is addressed: (i) process noise and sensor noises are correlated; (ii) sensor noises are contaminated by outliers or gross errors; (iii) cross-correlation among local estimates is unknown.First, to attack the correlation and outliers, a novel robust Kalman filtering (coR 2 KF) scheme with weighted matrices of innovation sequences is introduced as local estimator.It is shown that the proposed coR 2 KF takes both conventional Kalman filter and robust Kalman filter as a special case.Then, a novel version of our internal ellipsoid approximation fusion (IEAF) is used in the fusion center to handle the unknown cross-correlation of local estimates.

Problem Statement and Lemmas
Consider the discrete linear stochastic system with N sensors in the network: where () ∈ R  is the state vector,   () ∈ R   is the th measurement in the sampling period ; () ∈ R  is the disturbance input or system noise with zero mean and variance matrix where G  is the zero-mean Gaussian density and ΔG  is some unknown symmetric function representing the impulsive part of the noise density or outliers.
Remark 3. The problem of outliers is of practical importance in a target tracking system using multiple sensors (radar or infrared) [21,22] and communication applications where non-Gaussian (heavy-tailed) noise occurs, such as in underwater acoustics, and satellite communications through the ionosphere [23].It is commonly used to describe unmodeled measurement uncertainties that relate to the sensor failure, spikes, or jamming.
Assumption 4. The initial state (0) is independent of () and   (),  = 1, 2, 3, . . ., , and The problem is to find the estimation fusion x0 () of the state () in terms of local robust Kalman filter based on the measurement (  (), . . .,   (1)),  = 1, 2, 3, . . ., .The estimation fusion should have the desirable properties of efficiency and robustness; that is, it (i) yields the estimation fusion with a high accuracy for normal distributed observation while keeping the solution in high efficiency; (ii) reduces the bad effect of moderate errors on filtering and fusion in the way of weighting innovation; (iii) is robust in the sense that outliers and correlation do not affect the solution by setting the weighted matrix of innovation to be zeroes and further suppress the influence of outliers to the performance of fusion.
We start with some definitions and lemmas.Definition 5. Ellipsoid ( 0 , ) in   with center  0 and shape matrix  is the set where  > 0 might be standing for the covariance matrix of the estimation or fusion error.
Lemma 6.Under the Assumptions 1 and 4, for th sensor subsystem of the system (1)-( 2) without outlier (i.e.,   () = 0), the local optimal Kalman filters [24] one has: where Φ  () = () − Θ  ()  (), Θ  () = ()  () −1  , and   stands for a -dimensional unit diagonal matrix.  ( + 1 | ) and   ( | ) are the one-step prediction and filtering error variance, respectively, and   () is the filtering gain matrix.As described in Section 3, the innovation   () plays an important role in robust Kalman filter recursive process.Remark 7. If the process noise and measuring noises are not correlated, that is,   () = 0, it is obvious that Θ  () = 0, Φ  () = (), and (11) reduced to   ( + 1 | ) = ()  ( | )  () +   .This makes the correlated Kalman filter recursions in Lemma 6 taking the traditional Kalman filter in case of uncorrelated noises as a special case.Lemma 8. Checking whether two ellipsoids ( 0 1 ,  1 ) and ( 0 2 ,  2 ) have nonempty intersection can be cast as to the following Quadratic Programming (QP) problem with quadratic constraints [25,26]: where  1 and  2 are invariant with respect to affine coordinate transformation and describe the position of ellipsoids ( 0 1 ,  1 ), ( 0 2 ,  2 ) with respect to each other:  In another point of view, the conventional Kalman filter can also be formulated as a solution to a particular weighted least squares problem [21].Unfortunately, it is not robust because extreme outliers with arbitrarily large residuals can have an infinitely large influence on the resulting estimate.To handle this, the M-estimator, one of the most sophisticated approaches among the robust statistics approaches, is proposed [16].Further, the M-estimator has an advantage of less computational effort as it can be computed by a standard least squares algorithm with minor modifications [27].
M-estimators attempt to suppress the influence of outliers by replacing the square of the residuals with a less rapidly increasing loss function, which is where   (),   (), and ℎ  stand for theth row of   (),   (), and   , respectively (cf. ( 9)).(⋅) is a scalar robust International Journal of Distributed Sensor Networks convex function that has to cut off the outliers.Particularly, if one chooses (⋅) to be a quadratic function, the estimator according to (15) reduces to the least squares estimator or Kalman filter solution (7) [21,27].
Equating the first partial derivatives with respect to the state to be estimated () leads to the following relation: The score function (⋅) is usually nonnegative and symmetric, and (⋅), the derivative of (⋅), is often called the influence (score) function, since it describes the influence of the measurement errors on the solutions.Now, ( 16) can be rewritten as Letting (  ) = (  )/  , then ( 17) can be reformulated as the matrix form where   () = diag{( 1 ), ( 2 ), . . ., (  )}.
In the light of the above comparison and analysis of conventional Kalman filtering and M-estimator, the proposed coR 2 KF is given in Theorem 9 as follows.
Proof.The formula (19) can be derived from above directly, and the covariance of weighted innovation is from which we have the robust Kalman gain matrix as (20).
Substituting (19) into the filtering error equation where x ( + 1 | ) is the one-step prediction residual, and after mathematical manipulation, the robust filter covariance can be computed as This completes the proof.
Remark 10. (⋅) is a robust -estimate function for suppressing the outliers and is important for the estimation performance.Different (⋅) will result in different -estimate and fusion performance.Say, for a given density , the choice () = − log () yields the maximum likelihood estimator.Several robust cost functions have been used in the robust statistics setting, such as Huber's robust cost function, Andrews' method, Vapnik's loss function, or the biweight approach.Here, we propose a more general -estimate function generated from Huber's robust cost function: where  and  have to be chosen to provide the desired efficiency at the Gaussian model while possess robustness at the non-Gaussian model.Usually, they are chose empirically [18,22].For simplicity, we let  = 3 √  Remark 11.It can be seen that (⋅) is an even real-valued function and it is quadratic when   () is smaller than , which is just the same as the maximum likelihood (ML) cost function and keeps the efficiency of the M-estimate.For larger values of   () in the interval [, ], the function is linear and increase more slowly than ML.For residuals greater than , the function is equal to a constant.Based on (25), the weighted matrix of innovations can be formulated as The three different intervals of   (⋅) serve to deal with different kinds of residuals.In order to keep the accuracy and efficiency, when |  ()| ≤ ,   (⋅) are set to be 1; when sampling from the moderate innovations,   (⋅) are decreased with the residuals while sampling from a heavytailed distribution or outliers, the weighted matrix is set to be zeroes.
Remark 12.If () and   (),  = 1, 2, 3, . . ., , are correlated white Gaussian noises with variance matrices Q and   , respectively, we can see that   =   from Remark 11.In this case, the proposed coR 2 KF reduced to the correlated Kalman filter formulated in Lemma 6, which in turn takes the traditional Kalman filter with uncorrelated noises as a special case.

Robust Internal Ellipsoid Approximation Fusion
Once the local estimation is obtained by the subsystems, we are facing the problem of how to fuse the estimation in a right way in the higher level, that is, a cluster head or a sink node.In this paper, our internal ellipsoid approximation fusion (IEAF) method [15] is adopted to fuse local estimate.For convenience, the IEAF is reformulated in the following theorem.
Remark 14.After the parameters  1 and  2 are determined by Lemma 8, the center and shape parameters of the best internal ellipsoid ( x− 0 , P− ) can be calculated by ( 27)-( 28) and (30).
Remark 15.A graphic demonstration of the relation between the correlated source ellipsoids and the fused one by IEAF and CI are shown in Figure 1.It can be seen that the actual covariance ellipsoid always lies within the region defined by the intersection of the covariance ellipsoids of the fused sources regardless the degree of correlation between these sources.The CI, on the contrary, results in a covariance matrix that will always be greater than the actual one.This makes the CI some degree of conservativeness and a suboptimal estimate.The proposed IEAF slightly underestimates the intersection region, instead of overestimating this region like CI.In fact, the IEAF computes the largest ellipsoid contained within the intersection region and results in an increased performance.Also the consistency of the IEAF can be observed from Figure 1 directly.
Note that the center of internal approximation ellipsoid can be easily calculated from (27).Hence, the only difficulty remained is to solve (28) for the shape matrix.In this paper, we derive the explicit solution to (28) as explained in the following theorem in light of the symmetric positive-define property of  − ,  1 , and  2 .
Theorem 16.Given two local estimations  0 1 and  0 2 with the error covariance matrices  1 and  2 , respectively, according to the internal ellipsoid approximation method, the fusion estimation and its covariance are where and 0 ≤  1 ,  2 ≤ 1 can be calculated using Lagrange multiplier method from Lemma 8 and (30).
Proof.Note that all the parameters in ( 28) are given except  − .Letting and  =  − , now we are in the position to calculated  from the following equation: Noting that  and  are given symmetric positive-define matrices, we have that is,  −     = (−  ) −1  +      − .
Letting  =  −1  and using the properties of Kronecker product, we have where vec(⋅) represents vectorization function which is column stacking of the matrix.
Then  and  can be derived as follows, respectively, as the same of (35)-(37): where unvec(⋅) transfers an vector to corresponding matrix.

A Simulation Example
To evaluate the robustness of the proposed coR 2 KF and effectiveness of the IEAF approach, a simulation example of tracking a target moving on noisy circular trajectories is given in this section.The objectives of the simulation examples are twofold: (a) to verify the robustness of the proposed coR 2 KF and (b) to demonstrate the performance superiority of the IEAF method.Suppose N sensors are randomly deployed in the ROI, which is 50 m × 50 m with the coordinate from (−25, −25) to (25,25).The layout of the simulation scenarios is illustrated in Figure 2, where a five-point star stands for the location of a sensor.
Consider a target with the following dynamics: where  = [  ,   ]  denotes the state variable with   ,   stands for the target position on x-and y-direction, respectively,   = [ 0 −2; 2 0 ], and   = [1, 1]  .We use the discrete-time model of the target with parameters: The step-size is  = 0.025.Each sensor makes noisy measurements of the position of the target according to (2) with  =  2 .We add outliers into   (),  = 1, 2, all the simulation period.The outlier-corrupted noise   () is with the contaminated Gaussian density function F  = (1 − )(0,  2  ) + (0,  2  ), where   is set as 5.Moreover,   () are correlated with the process noise (), which is standard Gaussian white noise.The correlation coefficient   () is supposed to be [2,2] , [2,1], and [1,2], respectively.Our goal is to find the estimation fusion based on the local robust Kalman filters suppressing the effect of outliers on the estimation performance in case of cross-correlation among local estimates is unknown.We use the Root of Mean Square Error (RMSE), that is, RMSE = √ ∑  =1 [ x ()  x ()]/, as a performance criterion, where M is the number of Monte Carlo runs and x () represents the estimation error.
Scenario 1 (one sensor case, Robustness of coR 2 KF).To verify the robustness of proposed coR 2 KF under condition of both outlier-corrupted measuring noises and process-measuring noises correlation, we first consider that just a single sensor is used to track the target.For the coR 2 KF, we set  = 3  and  = 5  in (26).The RMSE over 500 Monte Carlo runs with different PDF function is shown in Table 1.
From Table 1, we can see in case of no outlier presented that the proposed coR 2 KF performs a little poorer than traditional KF.This is because coR 2 KF has deweighted the elements of matrix D in case of larger innovations.However, when outliers present, the adaptive weight of matrix D according to different innovation make the proposed coR 2 KF more robust than traditional KF.For example, in case of  = 0.2 and  = 1000, the degradation of coR 2 KF and traditional KF are, respectively, (2.97 − 2.27)/2.27= 30.84%and (6.04−2.25)/2.25 = 168.44%.Therefore, the performance discrepancy of coR 2 KF with different  and  is not very large, which demonstrates the robustness of proposed coR 2 KF.Scenario 2 (multiple sensors case, advantage of IEAF).In order to demonstrate the performance superiority of the IEAF method, we then consider multiple-sensor tracking of the target in a distributed way.Suppose that there are 3 sensors randomly deployed in the ROI as shown in Figure 2, while other parameters are the same as in Scenario 1.
In fact, the performance of IEAF has been shown to be less conservative than CI; both are not necessarily to have the knowledge of cross-correlation [15].Also, this can be observed from the results of this scenario, where the proposed coR 2 KF is adopted as the local estimator.The simulation is performed over 500 Monte Carlo runs each with 200 time steps.In this scenario, only the parameters  = 0.2 and  = 1000 are considered.In Table 2, the performance of the fast covariance intersection (FCI) [28], the simple convex combination (SCC), and the proposed IEAF are compared with the centralized fusion, which is optimal by directly fusing the measurements from local systems.
Obviously, using coR 2 KF as local estimators, proposed IEAF and the FCI achieve almost identical performance; both perform very closely to the optimal centralized fusion by coR 2 KF.Besides, both IEAF and FCI are much superior to the SCC approach, which ignore the cross-correlation among local estimates.To demonstrate the results visually and clearly, the RMSE of four fusion methods mentioned above versus time sampling periods is shown in Figure 3.Moreover, the real track along x-direction, measurement of each sensor and fused estimate by coR 2 KF-IEAF are shown in Figure 4.The results along -axis are similar, which are omitted for space reason.From Figure 4, we can see that the measurement of each sensor fluctuates seriously as the outliers appear.However, the outliers are suppressed by the coR 2 KF, as can be seen from Table 2 and both figures.In Table 2, we also give the results on average computational complexity over 500 Monte Carlo runs for whole 200 simulation steps.It is obvious that SCC has the smallest complexity since no cross-correlation needed.Our proposed IEAF has almost identical computational burden as FCI approach.This  can be expected since there is no extra computation except Y in Theorem 16.

Conclusion and Future Work
The problem of estimation fusion in a distributed WSN has been addressed under the conditions of correlation between process noise and sensor noises, outlier-corrupted sensor noises, and unknown cross-correlation among local estimates.To attack the correlation and outliers, a novel robust Kalman filtering (coR 2 KF) scheme with weighted matrices of innovation has been introduced as local estimator.It has been shown that the proposed coR 2 KF takes both conventional Kalman filter and robust Kalman filter as a special case.Then, an improved version of our previously proposed IEAF has been used in the fusion center to handle the unknown crosscorrelation of local estimates.To demonstrate the robustness of proposed coR 2 KF and the effectiveness of IEAF strategy, a simulation example of tracking a target moving on noisy circular trajectories has been included.Future work will be focused on posterior Cramer-Rao lower bounds (pCRLBs) analysis [29] and comparison with other outlier processing method (e.g., sensor validation technique) in the framework of estimation fusion [6].

Figure 1 :
Figure 1: Comparison among the proposed IEAF and CI, and actual covariance (AC).

Figure 2 :
Figure 2: A particular setup of the simulation scenarios.

Figure 4 :
Figure 4: The real track, each measurement, and the fused estimate along -axis.The measurements of the 3 sensors are denoted by SM1, SM2, and SM3, respectively, while the fused position estimate by the black dotted line, which tracks the real position on -axis very closely (cf. the red solid line).

Table 2 :
Comparison of fusion accuracy and average computational burden for 200 steps.