# A Consistent Track-to-Track Fusion Method Based on Copula Theory

Authors: Kelin Lu; K. C. Chang; Rui Zhou
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

This paper addresses the problem of distributed fusion when the conditional independence assumptions on sensor measurements or local estimates are not met. A new data fusion algorithm called Copula fusion is presented. The proposed method is grounded on Copula statistical modeling and Bayesian analysis. The primary advantage of the Copula-based methodology is that it could reveal the unknown correlation that allows one to build joint probability distributions with potentially arbitrary underlying marginals and a desired intermodal dependence. The proposed fusion algorithm requires no a priori knowledge of communications patterns or network connectivity. The simulation results show that the Copula fusion brings a consistent estimate for a wide range of process noises.


---

## Body

## 1. Introduction

The rapid growth of sensing and computational capabilities has advanced the development of distributed estimation [1] and information fusion [2] technologies for wireless sensor networks. Sensor networks provide greater coverage than single sensors by improving detection and estimation accuracy and extending geographical reach. The distributed fusion architecture [3] has many advantages over a centralized one such as lower bandwidth and increased reliability. In a typical distributed sensor network, estimates are often computed locally and transmitted to a fusion site through a communication network for fusion. However, one of the main challenges for distributed fusion in dynamic wireless sensing networks is to model the unknown correlation between the local estimates [4]. In optimal Bayesian fusion, the key is to identify the common information that has to be removed to avoid double counting [5]. The common information is usually from common target process noise or the information shared during the past communication.In a distributed estimation scenario where the local nodes process their local information and exchange their results (local tracks) to each other for fusion, the common process noise from the target of interest causes the dependency between the local estimation errors [6]. One solution proposed to solve this track-to-track correlation and fusion problem is to calculate the cross-covariances between the local estimation errors by communicating the local historical Kalman gains [7]. In addition to the common process noise, the local estimates may hold common information due to information exchange patterns in a distributed network. The Bayesian fusion equation [5] can be used to derive optimal fusion formulae for the state of interest as long as the measurements are conditionally independent given the state and where the common information can be extracted based on the processing architecture [8]. While the removal of duplicate information is straightforward in the theoretical formulation, identification of duplicate information for distributed estimation systems is difficult in dynamic networks. A recognized approach is to use graphical models [9] to represent dependent information with spatial and temporal relationships. A fundamental issue associated with this method is the amount of data needed to express correlations between the local estimates. In particular, long pedigree information might need to be communicated in order to achieve optimality. Since the optimal result is difficult to achieve in practice, several scalable fusion rules such as channel filter, the covariance intersection method [10], Chernoff fusion, and Bhattacharyya fusion [11] are proposed to reduce the communication loads. Furthermore, it is also difficult to combine data from multiple heterogeneous sources having different probability distributions, some of which may be non-Gaussian. For example, the dependency structure for Gaussian mixture estimates [12] has been demonstrated similarly to the linear Gaussian case [6]. Several approximate methods for the fusion of Gaussian mixtures have been presented based on Chernoff fusion [4, 13].Most traditional fusion approaches rely on assumptions of conditionally independent measurements. When such condition is not met, the Bayesian fusion equation is not optimal. In that case, a complete knowledge of the joint probability distribution of the observations is required in order to produce the optimal fusion results. However, the derivation of the joint distribution is intractable in general. We propose using Copula theory to model the correlations between local estimates in a consistent manner.The primary advantage of the Copula-based methodology is that it can characterize modal dependencies regardless of the respective marginal distributions [14, 15]. This property allows one to build joint probability distributions with potentially arbitrary underlying marginals and a desired intermodal dependence [16–18]. This is particularly suitable for distributed fusion with only knowledge of the marginal distributions given local sensor observations but not sufficient knowledge on their joint distribution.This paper proposes a new fusion method to deal with unknown correlation in a distributed sensing network. The method will function in a mathematically consistent manner while limiting data exchange and processing requirements. The paper is organized as follows. Section2 formulates the problem and Section 3 presents the Copula fusion method. The algorithm implementation details are presented in Section 4. Section 5 details the simulation results followed by the concluding remarks given in Section 6.


## 2. Problem Formulation

### 2.1. Kinematics Model

The Wiener process is widely used to model unknown inputs (maneuvers) in state estimation/tracking problem. Consider a 1D Wiener random process (or a discrete time Brownian motion) as follows:(1)xk+1=xk+vk,vk~N0,Q,where xk represents target position at time step k and vk is a white Gaussian process noise with variance Q.


### 2.2. Sensor Configuration and Measurement Model

The sensor configuration involves two range measuring sensors, where each sensori collects a set of measurements Zik=zi,1,…,zi,k about the target state x. The measurement equation is given by(2)zi,k=xk+wi,kwi,k~N0,Rifori=1,2,where zi,k represents the measurement at time step k for sensor i. wi,k is a white Gaussian measurement noise with variance Ri for sensor i.


### 2.3. Fusion Architecture

Typically, each sensor computes an estimate of the state of interestpx∣Zik based on its local observations Zik and sends the local estimation to a fusion node. The objective of the fusion node is to estimate the fused state given all the information available through communication, namely, px∣Z1k,Z2k. In this paper, we focus on hierarchical fusion where two sensors communicate their local estimates to the fusion site to be combined periodically. Denote Z1k and Z2k as the two information sets consisting of measurements. Then fusing state estimates are straightforward when these pieces of information are mutually exclusive and conditionally independent; that is,(3)pZ1k,Z2k∣x=pZ1k∣xpZ2k∣x.Then(4)px∣Z1k,Z2k=pZ1k∣xpZ2k∣xpxpZ1k,Z2k=1cpx∣Z1kpx∣Z2kpx,where the normalizing constant is given by(5)c=pZ1k,Z2kpZ1kpZ2k.When the conditional independence assumption is not met, the fusion equation (4) is only approximate. For example, in hierarchical fusion architecture, when two sensors communicate at less than full rate, the measurements are correlated given the target at the previous observation time due to the common target process noise [11]. In that case, a complete knowledge of the joint probability distribution of the observation conditioned on the target state is required in order to obtain the optimal results as shown in the following equation [5]:(6)px∣Z1k,Z2k=pZ1k,Z2k∣xpxpZ1k,Z2k.However, the derivation of the joint distributionpZ1k,Z2k∣x may be intractable in general. In the next section, we propose a method based on Copula theory to model the unknown dependence in order to derive the joint probability density.


## 2.1. Kinematics Model

The Wiener process is widely used to model unknown inputs (maneuvers) in state estimation/tracking problem. Consider a 1D Wiener random process (or a discrete time Brownian motion) as follows:(1)xk+1=xk+vk,vk~N0,Q,where xk represents target position at time step k and vk is a white Gaussian process noise with variance Q.


## 2.2. Sensor Configuration and Measurement Model

The sensor configuration involves two range measuring sensors, where each sensori collects a set of measurements Zik=zi,1,…,zi,k about the target state x. The measurement equation is given by(2)zi,k=xk+wi,kwi,k~N0,Rifori=1,2,where zi,k represents the measurement at time step k for sensor i. wi,k is a white Gaussian measurement noise with variance Ri for sensor i.


## 2.3. Fusion Architecture

Typically, each sensor computes an estimate of the state of interestpx∣Zik based on its local observations Zik and sends the local estimation to a fusion node. The objective of the fusion node is to estimate the fused state given all the information available through communication, namely, px∣Z1k,Z2k. In this paper, we focus on hierarchical fusion where two sensors communicate their local estimates to the fusion site to be combined periodically. Denote Z1k and Z2k as the two information sets consisting of measurements. Then fusing state estimates are straightforward when these pieces of information are mutually exclusive and conditionally independent; that is,(3)pZ1k,Z2k∣x=pZ1k∣xpZ2k∣x.Then(4)px∣Z1k,Z2k=pZ1k∣xpZ2k∣xpxpZ1k,Z2k=1cpx∣Z1kpx∣Z2kpx,where the normalizing constant is given by(5)c=pZ1k,Z2kpZ1kpZ2k.When the conditional independence assumption is not met, the fusion equation (4) is only approximate. For example, in hierarchical fusion architecture, when two sensors communicate at less than full rate, the measurements are correlated given the target at the previous observation time due to the common target process noise [11]. In that case, a complete knowledge of the joint probability distribution of the observation conditioned on the target state is required in order to obtain the optimal results as shown in the following equation [5]:(6)px∣Z1k,Z2k=pZ1k,Z2k∣xpxpZ1k,Z2k.However, the derivation of the joint distributionpZ1k,Z2k∣x may be intractable in general. In the next section, we propose a method based on Copula theory to model the unknown dependence in order to derive the joint probability density.


## 3. Distributed Fusion with Copula

Theorem 1 (Sklar’s Theorem for continuous conditional distributions [14, 15]).
LetH be a conditional bivariate distribution function with continuous marginals F and G, and let I be some conditioning set. Then there exists a unique conditional Copula C:0,1×0,1→0,1 such that(7)Hx,y∣I=CFx∣I,Gy∣I∣I,hx,y∣I=∂2Hx,y∣I∂x∂y=cu,v∣I·fx∣I·gy∣I,where u=Fx∣I and v=Gy∣I.Theorem1 indicates that the unknown joint density term in (6) can be constructed as(8)pZ1k,Z2k∣x=cu,v∣xpZ1k∣xpZ2k∣xu=FZ1k∣xv=FZ2k∣x,where u,v are marginal distributions and cu,v∣x is the Copula function. Substituting (8) into (6) we have(9)px∣Z1k,Z2k=cFZ1k∣x,FZ2k∣x∣xpZ1k∣xpZ2k∣xpxpZ1k,Z2k=cFZ1k∣x,FZ2k∣x∣x·px∣Z1kpZ1k/px·px∣Z2kpZ2k/px·pxpZ1k,Z2k=pZ1kpZ2kpZ1k,Z2k·cFZ1k∣x,FZ2k∣x∣x·px∣Z1kpx∣Z2kpx.Sklar’s Theorem states that a multivariate joint distribution can be written in terms of univariate marginal distribution functions and a Copula function which describes the dependence structure between the variables. Besides, when the marginal distribution functions are continuous, then the Copula function is unique.While there exists a unique Copula, its exact form is not available to us when we want to construct a joint distribution function with only marginal distribution functions. In this regard, we are supposed to choose a suitable Copula function to “approximate” the unknown true one. There are families of Copula functions, like Gaussian Copula, Archimedean Copula, and Student-t Copula. In this work, we choose Gaussian Copula since the local estimates are also Gaussian. The bivariate Gaussian Copula function [16, 17] is defined as(10)cu,v∣x≜11-ρ2exp⁡-ρ2Φ-1u2+Φ-1v2-2ρ·Φ-1uΦ-1v21-ρ2,where u=FZ1k∣x and v=FZ2k∣x. ρ∈-1,1 represents the correlation coefficient, Φ is the cumulative distribution function of N0,1, and Φ-1 is its functional inverse [17].This result implies that the Copula fusion equation (9) is similar to the original Bayesian fusion equation given in (4), the difference is the intermodal dependence term represented by the Copula function cFZ1k∣x,FZ2k∣x∣x, which is determined by the information of local sensor observation marginal densities.


## 4. Approximate Implementation of Copula Fusion

### 4.1. Pseudo Measurements

The exact implementation of fusion rule represented by (9) requires transmitting a full pedigree of the measurements from each sensor between communications, which is not practical in reality. In order to implement the Copula fusion rule in an autonomous manner for Gaussian Copula, we first approximate the information contained in the measurement history Zik=zi1,…,zik with a pseudo measurement zi,kpseudo; namely(11)FZik∣x≈Fzi,kpseudo∣x=121+erfzi,kpseudo-x2Rifori=1,2with(12)zi,kpseudo=x^i,k∣k-n+Pi,k∣k-n+RiPi,k∣k-n-1x^i,k∣k-x^i,k∣k-n,where n is the number of observations between communications. The derivation of (12) is given in the Appendix. With the pseudo measurement, we could approximate the Copula function (10) and derive the fused results (9) given the two local estimates.


### 4.2. Correlation Coefficientρ

The next issue to consider is how to derive the “consistent” correlation coefficientρ. When the system parameters, namely, the fusion rate, Q, and R, are given, we can perform exhaustive search over the grid in -1,1 and conduct extensive Monte Carlo simulations to find an appropriate ρ such that the resulting estimation variance corresponding to px∣Z1k,Z2k is consistent; namely, the result passes the NEES [18] test. However, this method only works in simulation, since the true state is not available in practice.In this paper, we conduct an off-line process to approximate the coefficientρ with the help of the analytical expression of steady state estimation variance for MAP fusion derived in [19]. A numerical search algorithm is applied to select the coefficient ρ such that the resulting estimation variance corresponding to px∣Z1k,Z2k is consistent with the theoretical result produced by MAP fusion.Note that (9) boils down to finding the joint density function of px∣Z1k and px∣Z2k. Since we do not keep the full pedigree of the measurements from each sensor, we use a single equivalent one to represent the history of measurements at the cost of losing information. MAP fusion is to find px∣x^1,k∣k,x^2,k∣k, since x^1,k∣k and x^2,k∣k only carry the cumulative effect from the past observations Z1k and Z2k, just like the partial information carried in the pseudo measurement, not the complete information in all measurements history, therefore we cannot reconstruct the optimal centralized fusion results with MAP fusion, and we can only get the suboptimal results. So theoretically, with an appropriate ρ, we are finding the joint density given the pseudo measurements so that the resulting estimate is consistent. On the other hand, MAP fusion is finding the best estimate given the two local estimates, which are equivalent to two pseudo measurements and the resulting fused estimate is also consistent.Although it is hard to prove the equivalence analytically, we show that the Copula fusion and MAP fusion obtain very similar results in a numerical manner. Here we conduct simulations at different fusion rate, where the coefficientρ is predetermined from Monte Carlo simulations and NEES tests.As shown in Tables1–3, the two fusion methods have similar results on both estimation variance and averaged MSE. For the reason that the steady state estimation variance for MAP fusion can be predicted in an analytic form, we could find an appropriate ρ such that the resulting estimation variance from Copula fusion is consistent with the one from MAP fusion.Table 1
Fusion rate = 1/2,Q=R=1, and ρ=0.1.
Estimation variance Averaged MSE Copula 0.3754 0.3749 MAP 0.3704 0.3723Table 2
Fusion rate = 1/4,Q=R=1, and ρ=0.25.
Estimation variance Averaged MSE Copula 0.3834 0.3822 MAP 0.3818 0.3806Table 3
Fusion rate = 1/8,Q=R=1, and ρ=0.35.
Estimation variance Averaged MSE Copula 0.3849 0.3841 MAP 0.3826 0.3829


### 4.3. Algorithm Summary

In this section, we summarize the Copula fusion algorithm. As shown in Figure1, at each fusion step, with the two local estimates, we first calculate the pseudo measurement zi,kpseudo, i=1,2, according to (12). Since the density function of the fused estimate px∣Z1k,Z2k cannot be derived in an analytical closed form, here we approximate the distribution with a Gaussian distribution Nxc,k∣k,Pc,k∣k through sampling.Figure 1
Copula fusion algorithm description.Suppose that, at time stepk, we generate a sampling set XS=x1S,x2S,…,xNS based on the prior probability distribution Nxck∣k-n,Pck∣k-n, where N is the number of samples. Then, for each element xiS in the sampling set, we can calculate the likelihood function (weight) of each sample point, namely, pxiS∣Z1k,Z2k, based on (9), where Fz1pseudo∣xiS and Fz2pseudo∣xiS can be obtained from (11) and the local posterior distributions and fused prior distribution(13)px∣Zik=Nx^ik∣k,Pik∣kpx=px∣Z1k-n,Z2k-n=Nxck∣k-n,Pck∣k-ncan be used to calculate pxiS∣Zik and pxiS∣Z1k-n,Z2k-n, respectively. Finally, we can approximate the resulting PDF as a Gaussian distribution and obtain its mean xc,k∣k and variance Pc,k∣k accordingly.


## 4.1. Pseudo Measurements

The exact implementation of fusion rule represented by (9) requires transmitting a full pedigree of the measurements from each sensor between communications, which is not practical in reality. In order to implement the Copula fusion rule in an autonomous manner for Gaussian Copula, we first approximate the information contained in the measurement history Zik=zi1,…,zik with a pseudo measurement zi,kpseudo; namely(11)FZik∣x≈Fzi,kpseudo∣x=121+erfzi,kpseudo-x2Rifori=1,2with(12)zi,kpseudo=x^i,k∣k-n+Pi,k∣k-n+RiPi,k∣k-n-1x^i,k∣k-x^i,k∣k-n,where n is the number of observations between communications. The derivation of (12) is given in the Appendix. With the pseudo measurement, we could approximate the Copula function (10) and derive the fused results (9) given the two local estimates.


## 4.2. Correlation Coefficientρ

The next issue to consider is how to derive the “consistent” correlation coefficientρ. When the system parameters, namely, the fusion rate, Q, and R, are given, we can perform exhaustive search over the grid in -1,1 and conduct extensive Monte Carlo simulations to find an appropriate ρ such that the resulting estimation variance corresponding to px∣Z1k,Z2k is consistent; namely, the result passes the NEES [18] test. However, this method only works in simulation, since the true state is not available in practice.In this paper, we conduct an off-line process to approximate the coefficientρ with the help of the analytical expression of steady state estimation variance for MAP fusion derived in [19]. A numerical search algorithm is applied to select the coefficient ρ such that the resulting estimation variance corresponding to px∣Z1k,Z2k is consistent with the theoretical result produced by MAP fusion.Note that (9) boils down to finding the joint density function of px∣Z1k and px∣Z2k. Since we do not keep the full pedigree of the measurements from each sensor, we use a single equivalent one to represent the history of measurements at the cost of losing information. MAP fusion is to find px∣x^1,k∣k,x^2,k∣k, since x^1,k∣k and x^2,k∣k only carry the cumulative effect from the past observations Z1k and Z2k, just like the partial information carried in the pseudo measurement, not the complete information in all measurements history, therefore we cannot reconstruct the optimal centralized fusion results with MAP fusion, and we can only get the suboptimal results. So theoretically, with an appropriate ρ, we are finding the joint density given the pseudo measurements so that the resulting estimate is consistent. On the other hand, MAP fusion is finding the best estimate given the two local estimates, which are equivalent to two pseudo measurements and the resulting fused estimate is also consistent.Although it is hard to prove the equivalence analytically, we show that the Copula fusion and MAP fusion obtain very similar results in a numerical manner. Here we conduct simulations at different fusion rate, where the coefficientρ is predetermined from Monte Carlo simulations and NEES tests.As shown in Tables1–3, the two fusion methods have similar results on both estimation variance and averaged MSE. For the reason that the steady state estimation variance for MAP fusion can be predicted in an analytic form, we could find an appropriate ρ such that the resulting estimation variance from Copula fusion is consistent with the one from MAP fusion.Table 1
Fusion rate = 1/2,Q=R=1, and ρ=0.1.
Estimation variance Averaged MSE Copula 0.3754 0.3749 MAP 0.3704 0.3723Table 2
Fusion rate = 1/4,Q=R=1, and ρ=0.25.
Estimation variance Averaged MSE Copula 0.3834 0.3822 MAP 0.3818 0.3806Table 3
Fusion rate = 1/8,Q=R=1, and ρ=0.35.
Estimation variance Averaged MSE Copula 0.3849 0.3841 MAP 0.3826 0.3829


## 4.3. Algorithm Summary

In this section, we summarize the Copula fusion algorithm. As shown in Figure1, at each fusion step, with the two local estimates, we first calculate the pseudo measurement zi,kpseudo, i=1,2, according to (12). Since the density function of the fused estimate px∣Z1k,Z2k cannot be derived in an analytical closed form, here we approximate the distribution with a Gaussian distribution Nxc,k∣k,Pc,k∣k through sampling.Figure 1
Copula fusion algorithm description.Suppose that, at time stepk, we generate a sampling set XS=x1S,x2S,…,xNS based on the prior probability distribution Nxck∣k-n,Pck∣k-n, where N is the number of samples. Then, for each element xiS in the sampling set, we can calculate the likelihood function (weight) of each sample point, namely, pxiS∣Z1k,Z2k, based on (9), where Fz1pseudo∣xiS and Fz2pseudo∣xiS can be obtained from (11) and the local posterior distributions and fused prior distribution(13)px∣Zik=Nx^ik∣k,Pik∣kpx=px∣Z1k-n,Z2k-n=Nxck∣k-n,Pck∣k-ncan be used to calculate pxiS∣Zik and pxiS∣Z1k-n,Z2k-n, respectively. Finally, we can approximate the resulting PDF as a Gaussian distribution and obtain its mean xc,k∣k and variance Pc,k∣k accordingly.


## 5. Simulation Results

To validate the proposed Copula fusion algorithm, a simulation scenario is developed with target and sensors models described in (1) and (2), respectively. The simulation is set up with a hierarchical fusion with feedback architecture. As shown in Figure 2, each sensor obtains the measurement and updates its local tracks with a sampling interval T=1 s. The fusion center fuses the local tracks with a period of nT. When the fusion process is completed, the sensors will receive the feedback from the fusion center.Figure 2
A hierarchical fusion scenario with a fusion period of2T.To verify the performance, we compare the results with the information matrix (IM) fusion algorithm. It is well known that, except in the full-rate fusion case (fusion after each measurement update), the IM fusion is only suboptimal and could be inconsistent [19]. 1,000 Monte Carlo simulation trials were conducted for each variance Q of the process noise varying from 10-4 to 103.Figures3–6 compare the MSE results between Copula fusion and IM fusion at full, 1/2, 1/4, and 1/8 rate, respectively. It can be seen that the estimation accuracies between the above two algorithms are very close.Figure 3
Positional MSE for Copula fusion and IMF at full rate (Q=R1=R2=1).Figure 4
Positional MSE for Copula fusion and IMF at 1/2 rate (Q=R1=R2=1).Figure 5
Positional MSE for Copula fusion and IMF at 1/4 rate (Q=R1=R2=1).Figure 6
Positional MSE for Copula fusion and IMF at 1/8 rate (Q=R1=R2=1).Figures7–9 compare the distributions of the fused result at steady state between Copula fusion and the IM fusion at 1/2, 1/4, and 1/8 rate communication, respectively. It can be seen that the two posterior distributions are not the same and the corresponding variance for IMF is smaller than that of Copula fusion.Figure 7
Probability density of fused result at steady state for Copula fusion and IMF at 1/2 rate (Q=R1=R2=1).Figure 8
Probability density of fused result at steady state for Copula fusion and IMF at 1/4 rate (Q=R1=R2=1).Figure 9
Probability density of fused result at steady state for Copula fusion and IMF at 1/8 rate (Q=R1=R2=1).Figures10–12 show the ratios of the perceived variance over the true variance. It can be seen that, except in the full-rate case, the fusion algorithm based on IM fusion is inconsistent when the variance of the process noise (target dynamic) is close to the variance of the measurement noise. Unfortunately, this is a critical operating range. As shown in the figures, the perceived variance is much smaller compared to the true ones for IM fusion.Figure 10
The ratio of perceived variance over the true variance (1/2 rate for Copula fusion).Figure 11
The ratio of perceived variance over the true variance (1/4 rate for Copula fusion).Figure 12
The ratio of perceived variance over the true variance (1/8 rate for Copula fusion).With the Copula fusion algorithm, the results shown in Figures10–12 demonstrate that the estimation results are near optimal and they are “consistent” for a wide range operating region with different process noises. This is because Copula function offers a flexible and reliable representation for modeling the unknown dependency between the local estimates. While the target dynamic and sensor observation models are assumed to be Gaussian in this paper, the proposed general method could be extended to model unknown dependency between potentially arbitrary marginal distributions to fuse data in a dynamic sensing network.


## 6. Conclusion

This paper presents a novel sensor fusion methodology based on Copula and Bayesian probabilistic theory for track-to-track fusion with unknown dependency. In the method, a mathematical characterization of the dependence structure of the local estimates is constructed using Copula statistical modeling. The recursive version of distributed Copula fusion is implemented by approximating the full pedigree of the local measurements with a pseudo measurement. The simulation demonstrates that, unlike other traditional approaches, the resulting fused estimates are consistent in the sense that the perceived uncertainty characterized by the estimation error variance is close to the true uncertainty. This is particularly important because overly optimistic uncertainty assessment could mislead a critical decision. Compared with the existing work, the proposed Copula fusion requires no detailed knowledge of communications patterns and is potentially applicable to fusion process with networked disparate sensors. A natural future research direction is to generalize the methodology coupled with the proven scalable “channel filter” algorithm to the fusion problems of heterogeneous and correlated measurements in ad hoc sensing networks.



---

*Source: 3751959-2016-12-19.xml*