A novel distributed fusion algorithm for multi-sensor nonlinear tracking

The covariance intersection (CI), especially with feedback structure, can be easily combined with nonlinear filters to solve the distributed fusion problem of multi-sensor nonlinear tracking. However, this paper proves that the CI algorithm is sub-optimal, thus degrading the fusion accuracy. To avoid such an issue, a novel distributed fusion algorithm, namely Monte Carlo Bayesian (MCB) algorithm, is proposed. First, it builds a distributed fusion architecture based on the Bayesian tracking framework. Then, the Monte Carlo sampling is incorporated into this architecture to form a feasible solution to nonlinear tracking. Finally, the simulation results verify that our MCB algorithm advances the state-of-the-art distributed fusion of nonlinear tracking.

Distributed fusion [1] refers to combining the information of decentralized sensors [2, 3], in which the observation information of each sensor is processed independently. In comparison with centralized fusion, the distributed fusion can significantly save time and storage resources in the fusion center. As for distributed fusion, the conventional solution is Bar-Shalom-Campo algorithm [4]. Later, Chang et al. [5] proved that this algorithm is not optimal in terms of the root-mean-square-error (RMSE), and they further presented an optimal fusion algorithm. Nevertheless, this algorithm can only be used in the scenario of two sensors fusion. Hence, Li et al. [6] proposed a unified fusion architecture, which can be used for two and more sensors. However, all aforementioned algorithms are designed with linear state space models (SSM), which results in the inferior fusion performance of nonlinear tracking.

Recently, some intelligent techniques, such as the adaptive fuzzy backstepping control technique [7] and the adaptive neural network technique [8] with nonlinear model predictive control, are proposed to advance the application of nonlinear models in multi-sensors or sensor network control [9]. Hence, the optimal distributed fusion with nonlinear model becomes one of the most important direction of signal processing. To our best knowledge, the most effective algorithms for fusion are based on covariance intersection (CI) [10, 11]. These CI-based algorithms can be easily combined with nonlinear filers to form the state-of-the-art solutions to nonlinear problems, such as UKF-SCI [12] and DPF-ICI [13]. The former is more accurate because the unscented Kalman filter (UKF) performs better in normally nonlinear tracking, while the latter has an advantage in non-Gaussian scenarios benefitting from the particle filter (PF). However, these CI-based algorithms require the computation on additional fractional powers to calculate the fusion results [14], which causes an increment of estimation error. Although many applications utilize a feedback structure¹ to improve the fusion accuracy [15], this paper proves that the sub-optimality caused by the increment of estimation error still exists in CI. To our best knowledge, few algorithms are proposed to overcome this sub-optimality. Hence, the present study is motivated to achieve a fusion algorithm which can outperform the CI algorithm. In summary, the novelty of this paper is that our algorithm overcomes the sub-optimality problem in covariance intersection (CI) fusion, which has not been addressed in the previous work.

In this paper, a Monte Carlo Bayesian (MCB) algorithm is proposed to overcome the disadvantage of the sub-optimality in CI algorithm, such that the fusion performance can be improved. Specifically, a distributed fusion architecture is designed based on the Bayesian tracking framework (BTF). This algorithm utilizes the law of total probability to form the distributed architecture, thus avoiding the increment of estimation error. Furthermore, the Monte Carlo sampling [16] is incorporated into our distributed architecture. This sampling method offers a direct approximate inference of the expectation of the target function with respect to a probability distribution [17]. Benefiting from this sampling method, the intractable problem of nonlinear estimation in BTF is circumvented by means of the approximation algorithm with sampling particles. Finally, based on the approximation, the fusion results are achieved in the form of the mean and variance of estimation at each step. The meanings of notation in the paper are listed in Table 1. The contributions of this work are listed as follows,

• A novel distributed fusion architecture is proposed, which makes full use of the information of different sensors to produce robust fusion estimation.

  Full size table

• A MCB algorithm is developed by means of the Monte Carlo sampling, which solves the nonlinear fusion problem based on the numerical approximation.

In this section, our distributed fusion architecture is developed based on BTF. Let us consider the I sensors fusion scenario with feedback structure. The state and the ith observation sequences are denoted as \(\left \{\boldsymbol {x}_{k};k\in \mathbb {N},\boldsymbol {x}_{k}\in \mathbb {R}^{d_{x}}\right \}\) and \(\left \{\boldsymbol {z}^{i}_{k};k\in \mathbb {N},\boldsymbol {z}^{i}_{k}\in \mathbb {R}^{d_{z^{i}}}\right \}\), where k is the time step, i={1,2,…I} is the number of sensors, d _x and \(\phantom {\dot {i}\!}d_{z^{i}}\) are the dimensions of the state and the i-th observation vectors, respectively. Their relationship can be represented by the SSM as follows,

where p(·) is the transition distribution, and q(·) is the likelihood of x _k in the i-th sensor. According to the BTF, the fusion process is to obtain the posterior distribution \(p\left (\boldsymbol {x}_{k}|\boldsymbol {z}^{f}_{1:k}\right)\) based on the SSM, where \(\boldsymbol {z}^{f}_{1:k} = \left \{\boldsymbol {z}^{i}_{1:k}\right \}^{I}_{i=1}\) is the integrated observation set combining all the sensors information from step 1 to k. Notice that, in the feedback structure, we normally have \(\boldsymbol {z}^{i}_{1:k}=\left \{\boldsymbol {z}^{i}_{k},\boldsymbol {z}^{f}_{1:k-1}\right \}\) because the prior observation information has been shared between sensors. Then, the fusion is calculated iteratively with two steps: the prediction and update steps.

The prediction step:

The update step:

where

is the marginal likelihood.

In normal cases, the current observation from each sensor is independent of others. Thus, the posterior distribution of (4) can be recalculated according to the law of total probability as follows,

Then, our distributed architecture is obtained by utilizing (6) to decompose the computation of (4) to each sensor, which is summarized in Fig. 1.

The distributed fusion architecture. The figure describes the fusion architecture in which sensors and fusion center are linked with information flow

Full size image

Figure 1 shows that, in the fusion process of time step k, the fusion center makes a prediction with (3). Then, the information of prediction is sent to each distributed sensor, in which the likelihood is calculated. Afterwards, each sensor sends the likelihood back to the fusion center to yield the final result with (6). Finally, the fusion result feeds back to the fusion center as the prior information of the next step.

This section analyzes the sub-optimality of the existing CI algorithm with feedback structure, by comparing it with the aforementioned distributed fusion in Section 2. This paper only concentrates on the most common scenarios that the distribution in SSM is Gaussian [18]. Without loss of generality, the fusion with I sensors is considered, where I≥2. According to [19], for any w∈(0,1), the posterior distribution of two sensors fusion in CI algorithm holds:

Lemma 1.

Let \(\boldsymbol {x}_{k}\in \mathbb {R}^{d_{x}}\), \(\boldsymbol {z}^{i}_{k}\in \mathbb {R}^{d_{z^{i}}}\), where \(k\in \mathbb {N}\); \(i=\{1,2,\ldots I\},I\in \mathbb {N},I\geq 2\); d _x and \({d_{z^{i}}}\phantom {\dot {i}\!}\) are the dimensions of vectors x _k and \({\boldsymbol {z}^{i}_{k}}\). Let \(\boldsymbol {z}^{f}_{1:k} = \left \{\boldsymbol {z}^{i}_{1:k}\right \}^{I}_{i=1}\), where \(\boldsymbol {z}^{i}_{1:k} = \left \{\boldsymbol {z}^{i}_{k},\boldsymbol {z}^{f}_{1:k-1}\right \}\). Given (7), the posterior distribution \(p\left (\boldsymbol {x}_{k}|\boldsymbol {z}^{f}_{1:k}\right)\) can be derived as follows,

where w _i ∈(0,1) and \(\sum _{i=1}^{I}w_{i}=1\).

Proof.

We decompose \(p\left (\boldsymbol {x}_{k}|\boldsymbol {z}^{f}_{1:k}\right)\) into a product form according to (7) as follows

where ω0′≡0, and \(w^{\prime }_{i}\in (0,1),~~i=1,2\ldots I\). Then, (9) can be rewritten as

where \(\omega _{i} \triangleq \omega '_{i}\prod _{j=0}^{i-1}(1 - \omega '_{j}),~~i=1,2\ldots I-1\), and \(\omega _{I} \triangleq \prod _{l=1}^{I}(1 - \omega '_{l})\). It can be easily known that w _i ∈(0,1) and \(\sum _{i=1}^{I}w_{i}=1\). Based on (10), we finally have

This completes the proof of Lemma 1.

In the feedback structure, Lemma 1 provides a generalized form of posterior distribution \(p\left (\boldsymbol {x}_{k}|\boldsymbol {z}^{f}_{1:k}\right)\) calculated in CI fusion. As seen from (8) and (6), the CI algorithm adds a fractional power ω _i to each likelihood.

Lemma 2.

Assuming that likelihoods obey Gaussian distribution, i.e., \(q\left (\boldsymbol {z}^{i}_{k}|\boldsymbol {x}_{k}\right) = \mathcal {N}\left (\boldsymbol {z}^{i}_{k}; h^{i}(\boldsymbol {x}_{k}),{\sigma ^{2}_{i}}\right)\) for each i={1,2,…I}, where h ⁱ(·) is the mapping from \(\mathbb {R}^{d_{x}}\) to \(\mathbb {R}^{d_{z^{i}}}\), and \({\sigma ^{2}_{i}}\) is the corresponding variance. Then, the variance of likelihood in each sensor becomes \({\sigma ^{2}_{i}}/\omega _{i}\) in CI algorithm.

Proof.

According to Lemma 1, the posterior distribution of CI is calculated with (8) as

Hence, the variance of each likelihood in CI fusion has turned into \({\sigma ^{2}_{i}}/\omega _{i}\) in the process of calculating the posterior distribution. This completes the proof of Lemma 2.

Since there exists

ω _i increases the variance of likelihood in CI fusion. Hence, compared with (6) in the distributed fusion of Section 2, the estimation (8) calculated in CI contains more uncertainty under the same observations. Such uncertainty leads to the increment of estimation error. As a result, the sub-optimality still exists in CI fusion with feedback structure.

To avoid the increment of estimation error, this paper utilizes the distributed fusion architecture in Section 2 to develop the fusion algorithm. In the most common applications, the observation is nonlinear [20]. Therefore, the posterior distributed in (6) cannot be simply calculated by joint Gaussian method with the procedure of squares completion [21]. To solve this problem, our MCB algorithm is proposed. This algorithm incorporates the Monte Carlo sampling into the distributed fusion architecture, by which the posterior distribution can be directly approximated with particles.

4.1 Monte Carlo sampling in MCB algorithm

At time k, N independent random particles are drawn from prediction distribution \(p\left (\boldsymbol {x}_{k}|\boldsymbol {z}^{f}_{1:k-1}\right)\) in (3):

where \(\boldsymbol {X}^{n}_{k}\) is the n-th particle drawn from \(p\left (\boldsymbol {x}_{k}|\boldsymbol {z}^{f}_{1:k-1}\right)\). According to the approximation inference method of numerical sampling [17], the prediction distribution \(p\left (\boldsymbol {x}_{k}|\boldsymbol {z}^{f}_{1:k-1}\right)\) can be approximated as follows,

In (14), \(\delta _{\boldsymbol {X}^{n}_{k}}(\cdot)\) denotes the Dirac delta function,

According to (14), the posterior distribution in (6) can be approximated as,

where C is the normalization constant.

4.2 Fusion with MCB algorithm

According to the distributed fusion architecture in section 2, the estimation of state x _k is calculated iteratively by (3) and (6).

First, in the fusion center, at time step k, given the previous state estimation with the mean \(\boldsymbol {\bar {x}}_{k-1}\) and variance P _k−1, the unscented sigma point set χ _k−1 is calculated as follows,

In (17), χ _i,k−1∈χ _k−1 denotes the ith point in the sigma point set, and λ is the scaling parameter [22]. Then, the prediction (3) is calculated in the form of mean \(\boldsymbol {\bar {x}}_{k|k-1}\) and variance P _k|k−1 [22] as

where f(·) is the state transition function, Q is the variance of corresponding noise, \({W^{m}_{i}}\in \boldsymbol {W}^{m}\) and \({W^{c}_{i}}\in \boldsymbol {W}^{c}\) are unscented transformation parameters calculated before estimation [22]. Then, N particles are drawn from the Gaussian approximation of prediction:

Second, the particles are sent to each sensor node. After achieving the observation \(\boldsymbol {z}^{i}_{k}\) in the ith sensor, the likelihood is approximated by the set: \(\left \{q\left (\boldsymbol {z}^{i}_{k}|\boldsymbol {X}^{n}_{k}\right)\right \}^{N}_{n=1}\). Then, each sensor sends its own approximated likelihood back to the fusion center.

Finally, the fusion center receives all the approximated likelihoods, and calculates the fusion results in the form of the mean and variance with (16) as follows,

In this section, the simulation results are provided to validate the performance of our MCB algorithm. In our simulation, a classic two-dimensional (2-D) fusion scenario with the nonlinear observation of one active and two passive radars is considered. In such a scenario, the fusion process is simulated by tracking a single target under several cases. For each case, different noise and kinematic models of transition equation are applied. Finally, 100 Monte Carlo simulations have been run for each case, and the fusion performance of nonlinear tracking is compared between our MCB algorithm and the UKF-SCI algorithm [12] with feedback structure. Note that, in this paper, we concentrate on Gaussian tracking scenarios, in which the UKF-SCI outperforms the DPF-ICI [13].

5.1 Simulation setup

In the common radar system, the transition Eq. (1) is considered as

where x _k =[d _x,k , d _y,k , v _x,k , v _y,k ]^T is the column vector, indicating the 2-D distance and velocity of a single target in the x−y plane at time step k. F is the state transition matrix with sampling interval s _T . In this paper, three state transition matrices are taken into account for three kinematic models: constant-velocity (CV) and constant-turn (CT) with known turn rates ω=2.5°/s and ω=5°/s [23]:

for CV model and

for CT model [24]. The trajectories of three kinematic models in our scenario are shown in Fig. 2.

The trajectories of target for three dynamic models. There are three sub-figures a, b, and c in this figure. a is the trajectory drew according CV model. b is the trajectory drew according CT model with turn rate ω=2.5°/s. c is the trajectory drew according CT model with turn rate ω=5°/s

Full size image

In addition, n _k =[n _d , n _d , n _v , n _v ]^T represents the state transition noise in distance and velocity. For fair comparison, this paper utilizes the same noise for different models which obeys the Gaussian distribution, i.e., \(n_{d} \sim \mathcal {N}(n_{d}; 0, {\sigma ^{2}_{d}})\) and \(n_{v} \sim \mathcal {N}(n_{v}; 0, {\sigma ^{2}_{v}})\). Here, σ _d and σ _v =2σ _d /s _T are standard deviations of transition noise of distance and velocity, respectively.

In our simulation, the observation equation in (2) is replaced by active and passive radar observation equations defined as,

where \(n_{\theta } \sim \mathcal {N}(n_{\theta }; 0,\sigma ^{2}_{\theta })\), \(n_{r} \sim \mathcal {N}(n_{r}; 0,{\sigma ^{2}_{r}})\); σ _θ and σ _r are standard deviations for tracking azimuth and distance. Then, we evaluate the performance of fusion with three radars, i.e., one active and two passive radars. Following [25], the aforementioned parameters are set as s _T =0.1(s), σ _θ =0.0001(r a d), σ _r =0.5(m), and σ _d ={0.25,0.3,0.35,0.4}(m).

The simulation setting may be applied to a composite guidance scenario [26], in which there are one active and two passive radars [27, 28]. Although we simplify the observation equation in a 2-D scenario, the simulation scenario is still suitable for the practical application of tracking and surveillance in network centric warfare (NCW) [29].

5.2 Evaluation

In our simulation, the single target tracking is performed with the aforementioned kinematic models. For each model, both MCB and UKF-SCI algorithms are utilized to fuse the information of three radars (one active and two passive) together, and obtain the final tracking result. Furthermore, in the tracking process of each model, the fusion performance is evaluated with four standard deviations of transition noise (i.e., {0.25,0.3,0.35,0.4}(m)). Note that 100 Monte Carlo runs are applied in our simulation for each scenario, and the particle number N in our MCB algorithm is set to be 200.

The simulation results of CV model are shown in Figs. 3 and 4. In Fig. 3, two scenarios with σ _d ={0.3,0.4} are selected to show the fusion root-mean-squared-error (RMSE) averaged over all 100 Monte Carlo simulations, along with the tracking time. Obviously, we can see from Fig. 3 that our MCB algorithm reduces the RMSE of both azimuth and distance, along with the tracking time. Such RMSE reduction becomes larger, when the transition noise increases. Furthermore, in Fig. 4, we show the average RMSE (from tracking time 3 to 50 s) of MCB and UKF-SCI algorithms with all scenarios of four transition noises. As seen from this figure, our MCB algorithm outperforms the UKF-SCI algorithm in all cases except the azimuth case of σ _d =0.25. In addition, the improvement of our MCB algorithm increases, when the transition noise becomes larger.

The fusion RMSE of azimuth and distance for both MCB and UKF-SCI algorithms, with CV model. In the upper figures, σ _d =0.3(m), and in the bottom figures, σ _d =0.4(m). There are four sub-figures a, b, c, and d in this figure, in which, under the CV model, the changes of fusion RMSE of azimuth and distance are described along with time steps. Meanwhile, a and b show the comparison between MCB and UKF-SCI algorithms under the case of standard deviations of distance σ _d =0.3(m); b and d show the comparison between MCB and UKF-SCI algorithms under the case of standard deviations of distance σ _d =0.4(m)

Full size image

The comparison of azimuth and distance RMSE at different standard deviations of distance, with CV model. There are two sub-figures a and b in this figure, in which, under the CV model, the changes of fusion RMSE of azimuth and distance are described along with different standard deviations of distance. The fusion RMSE in a and b is the mean of all results covering all fusion time steps

Full size image

The simulation results of CT model are shown in Figs. 5, 6, 7, and 8. Figures 5 and 6 show the fusion RMSE comparison of MCB and UKF-SCI algorithms with the known turn rate ω=2.5°/s. Figures 7 and 8 show the same things with the known turn rate ω=5°/s. All the results in these figures are similar to the ones in Figs. 3 and 4. Hence, our MCB algorithm also outperforms the UKF-SCI algorithm with CT models. In addition, our MCB algorithm can perform better when the turn rate increases. The details are shown in Tables 2 and 3.

The fusion RMSE of azimuth and distance for both MCB and UKF-SCI algorithms, with CT model of turn rate 2.5°/s. In the upper figures, σ _d =0.3(m), and in the bottom figures, σ _d =0.4(m). There are four sub-figures a, b, c, and d in this figure, in which, under the CT model of turn rate ω=2.5°/s, the changes of fusion RMSE of azimuth and distance are described along with time steps. Meanwhile, a and b show the comparison between MCB and UKF-SCI algorithms under the case of standard deviations of distance σ _d =0.3(m); b and d show the comparison between MCB and UKF-SCI algorithms under the case of standard deviations of distance σ _d =0.4(m)

Full size image

The comparison of azimuth and distance RMSE at different standard deviations of distance, with CT model of turn rate 2.5°/s. There are two sub-figures (a) and (b) in this figure, in which, under the CT model of turn rate 2.5°/s., the changes of fusion RMSE of azimuth and distance are described along with different standard deviations of distance. The fusion RMSE in (a) and (b) is the mean of all results covering all fusion time steps

Full size image

The fusion RMSE of azimuth and distance for both MCB and UKF-SCI algorithms, with CT model of turn rate 5°/s. In the upper figures, σ _d =0.3(m), and in the bottom figures, σ _d =0.4(m). There are four sub-figures a, b, c, and d in this figure, in which, under the CT model of turn rate ω=5°/s, the changes of fusion RMSE of azimuth and distance are described along with time steps. Meanwhile, a and b show the comparison between MCB and UKF-SCI algorithms under the case of standard deviations of distance σ _d =0.3(m); b and d show the comparison between MCB and UKF-SCI algorithms under the case of standard deviations of distance σ _d =0.4(m)

Full size image

The comparison of azimuth and distance RMSE at different standard deviations of distance, with CT model of turn rate 5°/s. There are two sub-figures (a) and (b) in this figure, in which, under the CT model of turn rate 5°/s., the changes of fusion RMSE of azimuth and distance are described along with different standard deviations of distance. The fusion RMSE in (a) and (b) is the mean of all results covering all fusion time steps

Full size image

Tables 2 and 3 depict the proportions of reduced RMSE for different kinematic models and noises. As seen in these two tables, the reduction of RMSE becomes larger when the transition noise increases. This is consistent with the results shown in aforementioned figures. Moreover, the reduction of RMSE with CT model of turn rate 5°/s is largest among all three kinematic models. That means, our MCB algorithm has more advantage for high maneuvering targets.

In summary, based on the BTF, our MCB algorithm makes full use of the information of all observations, and fuses it to obtain more accurate estimation on target tracking. Hence, compared with the CI algorithm, there are two advantages in our MCB algorithm: (1) When the transition noise is large, the fusion RMSE of azimuth and distance is still small; (2) When the turn rate is large, the fusion RMSE of azimuth and distance is small as well. In other words, in the high maneuvering cases such as large transition noise and turn rate, our MCB outperforms the state-of-the-art CI algorithm, in terms of the fusion RMSE.

5.3 Computational complexity

In this section, the computational complexity of our MCB algorithm is analysed. According to Section 4.2, in the fusion center, the algorithm contains two steps for each iteration: the prediction and update steps. In the prediction step, according to (18) and (19), the mean and variance are calculated with a summed form, in which the computational complexity is proportional to the dimension of tracking state. Moreover, N particles are drawn by sampling precess whose complexity is O(N). In the update step, firstly, we need to multiple all I sensors information of each particle to form a fusion likelihood. The computational complexity of this process is O(I). Then, the likelihood is used to compute the fusion results with the complexity being O(N), according to (20) and (21). In summary, the computational complexity in the update step is O(I·N).

To further evaluate the computational complexity of our MCB algorithm, we have recorded the computational time of the prediction and update steps, respectively, for one iteration in the simulation. Specially, the computer used for the test is with Intel Core i7-3770 CPU at 3.4 GHz and 4 GB RAM. In the aforementioned tracking case in Sections 5.1 and 5.2, the dimension of tracking state is 4, the particle number is 200 and the sensor number is 3. Through the simulation, we found out that the prediction and update steps take around 0.563 and 2.3 ms for one iteration. In other word, we only need 2.863 ms to compute the tracking results in the fusion center at each iteration. Hence, this algorithm is fast enough to utilize in the fields of radar tracking and fusion.

In this paper, we have proposed a novel MCB algorithm to achieve the distributed fusion estimation of nonlinear tracking. First, the distributed fusion architecture is set up based on BTF. Second, the sub-optimality in CI algorithms is proved. Then, to solve the estimation problem of nonlinear tracking, the Monte Carlo sampling method is incorporated into the distributed architecture. Benefiting from this sampling method, the approximation of fusion results is obtained through random particles. Simulation results verify that our MCB algorithm outperforms the state-of-the-art CI algorithm.

In summary, there are three directions of the future work in our paper. (1) Our MCB algorithm only offers a distributed calculation on the update step. Hence, a total distributed fusion structure is needed to further reduce the computation and communication overhead. (2) The time-discrete SSM used in our paper is actually a special case of continuous-time SSM. Therefore, we can extend our method to the exponential tracking scenario, in which the filtering can be processed with partially unknown and uncertain transition probabilities [30–33] with Markovian jump system. (3) Unknown inputs which represent the faults can be added into SSM and the residuals are calculated [34]. Hence, the fault detection algorithms [35] and the fuzzy model [36] can also be incorporated into our fusion system to strengthen the reliability of fusion process in the future work.

¹ The feedback structure means a structure which feeds the tracking result back to each sensor as the prior knowledge for the next step.

This work was supported by the National Nature Science Foundation of China projects under Grants 61471022 and 61573037.

Authors and Affiliations

1. School of Electrical and Information Engineering, Beihang University, XueYuan Road No.37, Haidian District, Beijing, 100191, China

   Jingxian Liu, Zulin Wang & Mai Xu

2. Collaborative Innovation Center of Geospatial Technology, 129 Luoyou Road, Wuhan, 430079, China

   Zulin Wang

Corresponding author

Correspondence to Mai Xu.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

This work is a product of inputs and work of all the authors to the research concept and experiment.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions