Data Fusion With Inverse Covariance Intersection for Prior Covariance Estimation of the Particle Flow Filter

The prior covariance estimation method based on inverse covariance intersection (ICI) is proposed to apply the particle flow filter. The proposed method has better estimate performance and guarantees consistent estimation results compared with previous works. ICI is the recently developed method of ellipsoidal intersection and is used to get the combined estimate of prior covariance. This method integrates the sample covariance estimate, which is unbiased but usually with high variance, together with a more structured but typically a biased target covariance through fusion gains. For verifying the performance of the proposed algorithm, analysis and simulations are performed. Through the simulations, the results are given to illustrate the consistency and accuracy of the proposed algorithm’s estimation and target tracking performance.


I. INTRODUCTION
Recent researches in particle flow filters (PFFs) have provided a solution to avoid degeneracy problems [1]- [3]. The solution is based on solving partial differential equations for migrating particles drawn from the prior distribution to the posterior distribution in the state space. In general, the bootstrap particle filter (classified as sequential importance resampling particle filters) draws particles from the prior distribution. It updates the weight of each particle using the likelihood of the latest measurement [4], [5]. As the Monte Carlo approximation of the posterior distribution is represented by a few particles, the weight degeneracy issue causes a poor representation of the posterior distribution [5].
There are various approaches to solve degeneracy problems effectively. However, one of the solutions involves separating the state space through factorization or partitioning [6]- [9]. These techniques are promising but rely on The associate editor coordinating the review of this manuscript and approving it for publication was Cesar Vargas-Rosales .
identifying a suitable factorization of the conditional posterior, so their applicability is restricted. A more general solution involves the incorporation of Markov Chain Monte Carlo (MCMC) methods within the particle filters [10]- [14]. However, MCMC methods are almost always computationally expensive, and their inclusion is difficult to apply as realtime filtering [5].
An alternative set of methods, called the PFF, can offer similar performance without the same computational requirements, but this comes at the cost of a more limited theoretical understanding. A framework for performing a progressive Bayesian update was introduced in [15]. In a series of papers [16]- [22], authors link the log prior and the log posterior distribution using derived partial differential equations (PDE) from guiding particles to flow from the prior distribution to the posterior distribution. Implementations of various PFFs have been reported in several publications. While conceptually being quite intuitive, PFFs can limit the estimated performance in practice due to several assumptions, made both in the theory and the implementation. In the previous VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ work [23], the key factors affecting the PFF performance are classified as pseudo-time discretization, ordinary differential equation (ODE) numerical solution, prior covariance estimation, and re-generating the particles set. Among the factors which affect the performance of the PFF, the availability of the prior covariance estimate is a key factor related to the performance of the PFF. The prior covariance estimate is used to calculate the particles' flow in the PFF. Thus, several types of research for estimating the prior covariance have been studied, and shrinkage estimationbased methods were proposed to obtain the covariance for dealing with the nonlinear model, non-Gaussian noise, and ill-condition of the covariance matrix [24]- [26]. Similar to the previous methods, we propose the prior covariance estimation method based on inverse covariance intersection (ICI) [27] which is the recently developed method of ellipsoidal intersection [27]. In this paper, ICI is used to obtain a combined estimate of prior covariance. This method integrates the sample covariance estimate, which is unbiased but usually with high variance, together with a more structured but typically a biased target covariance through fusion gains. The proposed method has better estimate performance and guarantees consistent estimation results compared with previous works of prior covariance estimation [28], [29]. Besides, extended parameter conditions of fusion gains are verified through performance analysis. The verification under the condition is satisfactory to explain that ICI-based covariance estimation is always better than the covariance shrinkage estimation.
An outline of the paper is as follows. In Section II, the whole process of the PFF is explained, and the proposed method with ICI based prior covariance estimation is introduced. Subsequently, in this section, ICI and its role in the proposed algorithm to improve the performance of the PFF are discussed. The analysis of the filter convergence is performed in Section III. In particular, the effect of the characteristics that changed the prior covariance written in Section II on the convergence of the filter in the prediction process of the proposed algorithm is analyzed. In Section IV, one simulation is performed using a representative nonlinear model. In addition, another simulation considers a multi-target tracking case with measurements of received signal strength obtained from several sensor nodes to confirm the actual filter performance compared with conventional methods. The conclusions are summarized in Section V.

II. PFF WITH ICI
This section provides a brief review of how deterministic particle flow can be used to address the nonlinear filtering problem with the particle filter. Besides, the prior covariance estimate using ICI, which is one of the important techniques for improving the performance of the PFF, is explained. The proposed estimation of prior covariance is used to calculate the particles' flow in the PFF, and its performance analysis is also performed.

A. FLOW TYPES OF PFF
After propagating particles using the dynamic model, a set of , and it represents the predictive posterior distribution at time k [5]. Particle flow, which is modeled as a stochastic background process in a pseudo time interval, λ ∈ [0, 1] is then used to migrate the particles so that they approximate the posterior distribution at time k. The stochastic process of i-th realization is defined as η i λ to design the particle flow (η i 0 = x i k , i = 1, 2, . . . N ). The zero diffusion PFFs(deterministic particle flow) [16]- [22] involve no random displacements of particles, and its flows are deterministic.
The trajectory of η i λ for i-th realization follows the ODE is written as where ς : R d → R d is governed by the Fokker-Planck equation and additional flow constraints [17]. The Fokker-Planck equation with zero diffusion is defined as where p η i λ , λ is the probability density of η i λ at the pseudo time λ of the particle flow. div ( ) refers to the differential operator (divergence). By imposing different constraints on the flow, (2) can lead to a variety of particle flow filters [5].
Among the various deterministic flows, the flow trajectory in the resultant exact Daum and Huang (EDH) [16] and the localized exact Daum and Huang filter (LEDH) [30] are the representative flow trajectory used in the PFF. To deal with nonlinear systems, a linearization of the filtering model is performed at the η i λ in the process of LEDH. For the i-th particle, the flow trajectory of LEDH becomes: where In the above equations, P is the predictive covariance (also called prior covariance) and R is the measurement covariance. H i (λ) refers to the measurement matrix and is expressed as h (η, 0) refers to the measurement model at the pseudo time λ = 0 and z is the measurement. Also, e i (λ) is the linearization error:

B. IMPORTANT FACTORS IN PFF
The equations (4) and (5) related to the flow require the availability of the prior covariance estimate (P). There are various methods to obtain the prior covariance, and one of the simplest ways is to estimate the covariance matrix using the prior particles, which is referred to as the sample covariance estimate (P). The sample covariance estimate is the result of the maximum likelihood estimate if the process data is Gaussian distributed and is an unbiased estimator of the true prior covariance (P). However, in the case of the nonlinear models or non-Gaussian noises, the Gaussian assumption may not remain valid andP could progressively get ill-conditioned [23]. An alternative method suggested by authors in [16] is to run an EKF/UKF in parallel to the PFF, and those filters generate the prior covariance matrix. While this approach is better than using the sample covariance estimate with the raw data, it ties the PFF estimation accuracy to that of the EKF/UKF [23] and this method may have limitations when applied to a system with high nonlinearity or non-Gaussian case because of the Gaussian assumption used in EKF/UKF structures. In addition, EKF causes an underestimate or overestimate situation due to error covariance update with Kalman filter iteration and therefore risks becoming inconsistent in the statistical sense. In the case of UKF, it has the limitation that it does not apply to general non-Gaussian distributions. Besides, prior covariance results obtained from EKF or UKF could also exhibit a wide spread of the eigenvalues [23] (this feature does not guarantee tight boundaries of error covariance).
To deal with those problems, there is another approach used in the multivariate statistics literature for the estimation of the covariance matrices, known as the shrinkage estimation. The use of such methods dates back to the work of Stein [31]. The main idea is to merge the raw estimate of the sample covariance (P), which is close to unbiased but typically with high variance, together with a more structured but typically a biased target (P − ) through a scale factor (ρ, also called fusion gain), to get the combined estimate covariance. There are several shrinkage estimators mentioned in the literature [23], with different target covariance matrices. Shrinkage estimators are defined through a convex combination of the matrices P − andP. The objective becomes to find an optimal shrinkage intensity that minimizes the cost function min ρ E P * −P 2 (6) where P * = ρP − + (1 − ρ)P The covariance shrinkage estimation can be categorized as a kind of covariance intersection (CI) method in terms of multiple data fusion. CI Algorithm is proposed to combine two quantities in the presence of unknown correlation, and to provide an appropriate estimate of the resulting covariance matrix. In addition, the CI algorithm requires ρ to be optimized at every step by minimizing the cost function (for example, trace or the determinant of error covariance). Thus, since the configuration of the CI algorithm is similar to the shrinkage estimation method described above, it is shown in this paper as a CI algorithm for convenience of performance comparison.
CI is introduced in the previous work [32] and is a representative fusion rule. There are various researches on CI, and it has been implemented in applications [27]. According to the analysis results in the previous work [33], it is shown that CI tightly bounds the entirety of possible error covariance matrices. It means that if the correlations between two estimation values to be fused are entirely unknown, CI encompasses the optimal fusion method in terms of a minimum mean squared error and also other optimality criteria [27]. However, CI often provides too conservative fusion results as typical estimation tasks, in general, prevent extremal correlation terms from occurring [27].
To obtain a less conservative result compared to CI, ellipsoidal intersection (EI) is proposed [29], and it is applied to many applications [34], [35]. However, EI does not guarantee consistent fusion results in the presence of unknown common information [27]. Thus, in recent years, the novel approach achieves a conservative fusion result by computing the bound on the intersection of inverse covariance ellipsoids, which gives reason to name it inverse covariance intersection [27].
Thus, we propose the prior covariance estimation method based on inverse covariance intersection. Still, the proposed method has better estimate performance and guarantees consistent estimation results compared with previous works based on the covariance shrinkage estimation method.

C. PERFORMANCE ANALYSIS OF ICI COMPARED WITH FILTER PREDICTED VALUE
A consistent combination of the estimates x − k , P − k obtained in the filter prediction process and x − k ,P − k obtained from the sample mean and covariance of particles is provided by x − k ,P − k with ICI as follows [27]: for any ρ ∈ [0, 1]. The gains in (31) are set as In (7) and (8), the estimates of state variables and its covariance x − k , P − k are obtained by Kalman filter structure as follows: where f k () is the nonlinear system model of the filter and x + k−1 refers to the posterior estimation results of the state variables at time k − 1.
is the Jacobian matrix of the system model, P + k−1 is the posterior error covariance at time k−1, and Q k is the noise covariance of the system model. The posterior estimation results x + k−1 , P + k−1 are calculated by after weight update and resampling process of the particles in the PFF as follows: where x i k refers to i−th particle of the state variables indicates the weights of the i−th particle at time k − 1. The update process of weights is similar to that of the PF and is written at line 20 of Table 1. After passing through the resampling process, each weight value becomes 1/N shown at line 28 of Table 1.
In the case of the sample mean and covariance of particles are calculated aŝ Finally, the result of prior covarianceP − k written in (8) is substituted into the P in (4) and (5).
According to lemma 9 in the previous work [27], it is verified that the ICI approach provides more accurate fusion results than CI under the specific condition: ρ ICI = 1 − ρ CI . ρ ICI , ρ CI refer to the scale parameters in ICI and CI, respectively. In this particular condition, ICI based covariance estimation method can have better estimation performance than the covariance shrinkage estimation (CI-based covariance estimation method). The verification under the condition is unsatisfactory to explain that ICI-based covariance estimation is always better than the covariance shrinkage estimation. Therefore, in this paper, we have found extended parameter conditions related to ρ ICI , ρ CI satisfying the verification result through additional performance analysis.
A simple analysis is performed to verify that using ICI based covariance estimation method is better than using the covariance shrinkage estimation.
If there is a parameter ρ ICI ∈ [0, 1] for each ρ CI ∈ [0, 1], comparison of the two covariances is explained as Equation (17) is converted into a diagonal matrix form for ease of analysis. Multiply the transformation matrix T k on both sides of components in (17) to convert into a diagonal matrix as follows: are the diagonal matrices. T k is a transformation matrix which can be computed with the aid of an eigenvalue decomposition as in [27].
If the equation (19) is rearranged by the last term , always positive), (19) can be written as follows: By using inequality a 2 + b 2 ≥ 2 |ab|, it follows: The expression in braces in (20) can be expressed as Finally, the diagonal entries of the two covariances' difference are written as If the system we're dealing with is nonlinear and non-Gaussian, then the inequality d − k j ≥ d − k j is likely to be valid. Thus, if 0 ≤ ρ ICI + ρ CI ≤ 1, the diagonal entries of the two covariances' differences are always positive ( D k jj ≥ 0). This means that ICI-based covariance estimation is always more accurate than CI-based covariance estimation. Table 1 shows the overall structure of the PFF proposed, including ICI-based covariance estimation. In this paper, the particle resampling process is performed by the systematic method [36] as an option as in Algorithm 1. The resampling process is written in Table 1 is performed when the effective number of particles (N eff ) [37] obtained by the equation of line 26 in Table 1 is less than the maximum number of particles (N max = N /2). The user may use the resampling technique that has been studied to improve computational efficiency. An adaptive resampling method called KLD-resampling [38], [39] is proposed that determines the number of particles to resample based on the Kullback-Leibler measure of the fit of the posterior distribution represented by weighted particles for increasing efficiency of the algorithm.
In this paper, the performance analysis of the proposed covariance estimation technique is given priority, and the performance improvement by the resampling technique is not analyzed. However, in order to efficiently apply the proposed algorithm to real systems, it is necessary to study the optimal resampling technique of the PFF in the future.

III. CONVERGENCE OF THE PFF
In the previous works [40], [41], it is possible to find bounds for the mean square error of the SMC-PHD filter at each state of the algorithm. These depend on some considerations related to the designed parameters of the target system's model. Besides, the analysis of the filter convergence is based on the fact that the sum of two sequences converges weakly to the sum of the limits of those sequences, which follows from a basic result of real analysis on the convergence of sequences of real numbers [40]- [43]. In this paper, those analysis results are extended and applied to the proposed algorithm to verify VOLUME 8, 2020 the convergence of the proposed algorithm. In particular, the effect of the characteristic that changed the prior covariance written in Section II on the convergence of the filter in the prediction process of the proposed algorithm is analyzed. Before proving the convergence of the proposed filter, some explanation and consideration are given as follows.
If θ N is a sequence of measures that depend on the number of particles, N , then θ N converges to θ when ∀ϕ ∈ B R d [41].
where B R d is the set of bounded Borel measurable function on R d . d is the dimension of the space. When the measure in the inner product . , . is discrete, it defines the summation inner product, In (24), D N k|k is the density of the PFF and refers to To prove convergence of the proposed algorithm, D k|k−1 is defined as the density propagated from the previous time step (the prediction process of the filter), and the boundary of D N k|k , ϕ − D k|k−1 , ϕ is analyzed. By the triangle inequality, where η k refers to the stochastic process of realization (propagated particles' density) at the pseudo time, λ = 1. Let ς k−1 be the σ -algebra generated by the particles x i . Then, the first term of the right side of (26) is expected as The above equation is expressed in sum using the independence of each particle as follows: According to (28), the boundary of the equation is expressed as [40], [41] Using Minkowski's inequality, (26) is rearranged as If there are no newly created targets at time k, the boundary of the (30) is finally expressed as (31) In (31), c k−1|k−1 refers to the boundary of the mean square error at the previous time step k − 1 and is obtained from the assumption that [42] In addition, some assumptions are held for verifying the boundedness of η k written in (29). There exist real constants which are related to the system model, measurement model, and the following bounds on matrices of filter models are satisfied for every time index, k as follows: The boundary of η k is expressed as When η 1 ≤ η max , the boundaries of (I + λA (1)) η 1 and λb (1) are written as (36) and (37), as shown at the bottom of the next page, Thus, the boundary of η k is related to the bounds on matrices of filter models. In the case of the proposed algorithm, it is confirmed that the convergence of the mean square error depends on N , the boundary of η k and c k−1k−1 .

IV. SIMULATIONS
Two simulations are performed to analyze the performance of the proposed algorithm. In the first simulation, a representative nonlinear model is used, and in the second simulation, multiple nonlinear measurements are considered in multitarget tracking.

A. UNIVARIATE NONSTATIONARY GROWTH MODEL
A nonlinear model called univariate nonstationary growth model (UNGM) is mainly used as a benchmark in the previous works [44], [45]. The reason why this model is widely used is highly nonlinear and bimodal properties, so it is challenging for conventional filtering methods.
Thus, to illustrate some of the advantages of the proposed algorithm, this model is used in the first simulation for verifying the performance of the proposed algorithm.
The dynamic system and measurement model for UNGM can be written as where v k is system noise and its distribution is zero mean and one variance of Gaussian distribution N (0, 1). Also, w k refers to the measurement model with Gaussian distribution N (0, 1). Some previous works are used for comparing and evaluating the performance of the proposed algorithm through this simulation case as follows: the PF with the different number of particles (1000, 10K), the EDH PFF [5], [16] and LEDH PFF [5], [30] with the sample covariance estimation, LEDH PFF with EKF for estimating the predictive covariance, LEDH PFF with CI which combines the sample covariance using PF and the predicted covariance using EKF, and LEDH PFF with ICI which combines the sample covariance obtained by PF and the predicted covariance obtained by EKF (proposed algorithm). The number of particles in all Monte Carlo based algorithms is set to 1000 (N = 1000). In Fig. 1, the absolute errors and 3σ confidence intervals of each filtering methods result in the case of the first simulation. In this case, it seems to be very crucial, as the model is highly nonlinear and multi-modal. It can be shown that conventional methods are overoptimistic in many cases. While the proposed algorithm has better estimation performance when its estimation results are unreliable (when 3σ values increase). In addition, the lower estimation error of the proposed algorithm can be seen from Fig. 1. Lastly, the root mean square errors (RMSE) of each simulated method of 1000 Monte Carlo runs are calculated in Table 2. Likewise, the proposed algorithm is superior over all other simulated methods because of multiple covariance value fusion based on ICI even though estimation results are unreliable.

B. MULTI-TARGET TRACKING
A multi-target tracking case is considered with a relatively large state space and highly informative measurement based on the previous paper [5] for evaluating the performance of the proposed algorithm and comparison with the previous works. The simulation setup is the same as the previous study [5], and only the measurement model is changed to the received signal strength. In the simulations, the number of targets is set to 4 (c = 4), and targets independently move in the simulation region [5]. State variables of the tracking system are set x c k y c kẋ c kẏ c k (2D position and velocity of the c-th target) and a system model of the multi-target tracking case is modeled as where the system model is the constant velocity model and the matrix of the system model can be expressed as In equation (38), v c k is the c-th target's system white noise of which the probability density function (pdf) is assumed to be the Gaussian pdf with zero means and pre-designed covariance matrices Q. At each time step, all targets receive the signal sent from the sensors and theirs transmit signal power P c 0 . Attenuated signal's powers are measured by the c-th target's system. Thus, the measurement model for the c-th target located at x c k , y c k with the s-th sensor located at (x s , y s ) is designed as [46] Equation (42) is based on the free space model between transmitter and receiver in the radio frequency communication system [46], [47]. This model is generally used to explain free space loss with some conditions: the isotropic transmit antenna of sensors, which radiates signal equally in all directions. The isotropic receive antenna captures power equal to the density times the area of the antenna when the ideal area of the antenna is set to λ c 4π (λ 2 c is the wavelength of the signal's carrier ) [47]. In assumption, the power density of the effective area is set to [47]. Where is the Euclidean norm, and the number of sensors (N s ) is set to 25 located at grid intersection within the simulation region, as shown in Fig. 2. Finally, the received power of the signal came from the sensor is obtained by multiplying the ideal area of the antenna, and the power density of the effective area, as shown in (42). In addition, the measurements are perturbed by the white noise of which the probability density function (pdf) is assumed to be the Gaussian pdf with z s (x k ) and variance σ 2 s . The filters used in this simulation for performance comparison are the same as those used in the first simulation.
The optimal mass transfer (OMAT) [48] is selected as the comparison metric in this simulation case with a fixed number of target tracking situation. The OMAT metric d p X ,X between two arbitrary sets X = {x 1 , x 2 , . . . , x c } andX = x 1 ,x 2 , . . . ,x c is defined as where the scalar p is set to 1, so the OMAT metric assigns targets using the permutation that minimizes the Euclidean distance d to the true target positions. is the set of possible permutations of {1, 2, · · · , C}, and d x c ,x c is the Euclidean distance between x c andx c .       algorithm. The algorithms used in the comparative analysis are as follows: EKF, PF (particle filter) [37], EDH PFF [5], [16], LEDH PFF [5], [17], LEDH PFF with CI, and LEDH PFF with ICI (proposed).
According to Fig. 3, the EDH PFF is less accurate than the LEDH PFF because the proposal distribution constructed using the EDH flow does not provide an excellent match to the posterior distribution of the nonlinear model. In addition, the LEDH PFF has the operational advantage brought by the importance sampling step in the filter. There is a performance difference of the PFF algorithm depending on the presence of a prior covariance estimation algorithm in the PFF. Especially, using the ICI based algorithm shows a more accurate estimation performance than the case of using the CI based algorithm. Thus, the proposed algorithm (LEDH PFF with ICI) has the smallest tracking error and reduced the average OMAT below 1 meter with a few time steps.
However, even though the ICI based covariance estimation scheme is applied to the LEDH PFF for improving the estimation performance, the efficiency of the proposed algorithm is lowered because of the increased amount of computation by the additional algorithm. Table 3 shows the average OMAT and the computation time of each algorithm. In order to apply the proposed algorithm efficiently to the nonlinear system model, it is necessary to study a technique that can reduce the amount of computation. Although the tracking performance of multi-targets is compared using OMAT parameters, the target tracking errors of each simulated filter are shown in Figs. 4, 5, 6, and 7. By comparing the estimation error results of the position for each target trajectory, it is shown that the proposed algorithm tracks dynamic objects well.

V. CONCLUSION
In this paper, the data fusion method between sample covariance and estimated covariance obtained from the EKF using ICI is proposed to apply the particle flow filter. The proposed method has better estimate performance and guarantees consistent estimation results compared with conventional algorithms such as the covariance intersection based method and EKF. To verify the performance of the proposed algorithm, analysis and simulations are performed. According to analysis in Section III, the analysis of the filter convergence is based on the fact that the sum of two sequences converges weakly to the sum of the limits of those sequences. It is confirmed that the convergence of the mean square error depends on the number of particles and the boundary of particle flows. Besides, through the simulations in the case of multiple target tracking using received signal strength from the multiple sensors, the simulation results are given to illustrate the accuracy of the proposed algorithm's estimate performance compared to other conventional algorithms. However, even though the ICI based covariance estimation scheme is applied to the PFF for improving the estimation performance, the efficiency of the proposed algorithm is lowered because of the increased amount of computation by the additional algorithm. To apply the proposed algorithm efficiently to the nonlinear system model, it is necessary to study a technique that can reduce the amount of computation.