Distributed Box Particle Filtering for Target Tracking in Sensor Networks

Distributed target tracking is a significant technique and is widely used in many applications. Combined with the interval analysis, box particle filtering (BPF) has been proposed to solve the problem of Bayesian filtering when the uncertainties in the measurements are intervals; that is, the measurements are interval-based vectors. This paper is targeted for extending the existing BPF based on a single sensor to a distributed sensor network. We propose a distributed BPF (d-BPF) that each sensor communicates with its direct neighbors to collaboratively estimate the states of the target. The feasibility of the proposed distributed BPF is justified, and some numerical simulations are presented to show its effectiveness in target tracking.


Introduction
Target tracking is a prevalent task in many scenarios, such as navigation, military, and transportation. In the literature, various approaches have been proposed for target tracking [1,2]. One of the most popular approaches is the Bayesian inference [3]. Within the Bayesian framework, the posterior probability density function (pdf) of the states of the target is calculated, conditioned on the available measurements obtained by sensors. Then, based on the obtained posterior pdf, the states of the target can be estimated, correspondingly.
Reviewing the related works of target tracking based on Bayesian inference, the Kalman filter (KF) [4][5][6] is a classical approach. However, in the presence of the nonlinearities and non-Gaussian uncertainties, the KF becomes inefficient. To tackle this problem, extended trackers such as the extended Kalman filter (EKF) [7] and the unscented Kalman filter (UKF) [8] have been developed. Although these extended trackers can handle some nonlinear cases, they become ineffective when the high-variance measurements and state noise combined with nonlinearities exist [9]. The particle filter (PF), which represents the posterior pdf of the states by sampling, shows outstanding performance, especially for nonlinear and/or non-Gaussian cases [3]. However, the efficiency and accuracy of the PF largely depend on the number of the particles. In the case where a large number of particles are used, the computational complexity is quite high.
It is worth pointing out that all of the above-mentioned approaches can only deal with the point measurements represented by fixed vectors. However, for a large amount of target tracking problems, the uncertainties may lie in a certain region, and thus the measurements are intervals. To deal with such interval-based uncertainties, the box particle filtering (BPF) [10][11][12] based on the interval framework is proposed, which propagates a set of weighted boxes in a sequential way. Instead of using a large amount of point samples, much less weighted boxes are used to approximate the desired moments of the posterior pdfs.
On the other hand, distributed in-network target tracking, where a set of distributed sensors collaboratively estimate the states of target, has received much attention recently. Distributed tracking can achieve low energy consumption as well as high robustness and scalability. Owing to these merits, various distributed target tracking algorithms have been proposed, such as the distributed KF [13][14][15] and the distributed PF [16][17][18][19].
However, the above-distributed algorithms are only used for the cases where the system measurements are point-based vectors. Decentralized BPF that can handle interval-based uncertainties has not aroused much attention. Although 2 International Journal of Distributed Sensor Networks a decentralized BPF algorithm using belief propagation (BP) like message-passing algorithm is proposed in [20], it is only focused on the context of target localization but not applied to target tracking. Considering this, a distributed BPF (d-BPF) algorithm is proposed in this paper, where only the onehop data exchange and communication between sensors and their neighbors are used to collaboratively track the states of the target. Theoretical studies and numerical simulations are then performed to show the effectiveness of the proposed algorithm.
The rest of this paper is organized as follows. In Section 2, we review some basic knowledge related to BPF. In Section 3, we derive the d-BPF and analyze the feasibility of the d-BPF. The selection of the combination weights for cooperation is also introduced in the same section. Section 4 presents some simulation examples. Finally, some conclusion is drawn in Section 5.

Preliminaries
The traditional Bayesian estimation deals with the point measurements, while the BPF can cope with interval-based measurements. The key idea of the BPF is to replace a particle by a multidimensional interval in the state space [11]. Therefore, the derivation of BPF depends on the theory of interval framework. So, before introducing the BPF, we overview some basic knowledge of interval analysis, including interval arithmetic and inclusion function, as well as constraint satisfaction problems and contraction.

Interval arithmetic. A real interval [ ] is a nonempty set of real numbers:
where and denote the lower and upper bounds of [ ], respectively. From this definition, it is natural to define settheoretic operations, such as intersection and union, and the results of these must be intervals also. For instance, the interval union of two intervals [ ] ⊔ [ ] is defined as where the symbol "∪" denotes the union on sets and the symbol "[⋅]" denotes the interval hull operator returning the smallest interval enclosing , for any set in real set R [25].

The interval intersection of two intervals [ ] ⊓ [ ] is defined as
where the symbol "∩" denotes the intersection on sets.
Interval arithmetic operations are also defined on IR such that the interval result encloses all possible real results, where IR denotes a set of all real-valued intervals. Given [ ] = [ , ] and [ ] = [ , ], the binary operations {+, −, ×, /} can be defined by where ⊗ denotes the binary operations.
Some basic interval arithmetic operations are summarized in Table 1 [25]. An interval vector [x] is defined as a vector with interval components and the space of all dimension interval vectors is denoted by Similarly, interval arithmetic to boxes can be extended straightforwardly [25].

Inclusion Function.
As mentioned above, in the BPF, each particle is replaced by a multidimensional interval in the state space. In order to keep consistency, it is necessary to guarantee that the result of interval operation is still an interval or an interval vector. Fortunately, an inclusion function can meet our need. To make it simple, assume f is a function from R to R . In the BPF, when propagating box particles via the continuous nonlinear function f, the mapping of the box particle is not necessary a box. But it can be approximated by an inclusion function. An interval function [f] from IR to IR is said to be an inclusion function for f if An inclusion function may be very pessimistic. But the minimal inclusion function for f denoted by [f] * is unique [11].

Constraint Satisfaction Problems and Contraction
Define the vector x as x = ( 1 , . . . , ) , where denotes transposition, and let f is a vector function whose components are , = 1, . . . , . Then, the constraints can be rewritten as This corresponds to a constraint satisfaction problem (CSP), which is denoted as H; that is, and the solution set of H is defined as Note that the set S is not necessary a box. Within the interval framework, the term contracting refers to solving a CSP, that is, finding the smallest box [x] ⊂ [x] that constitutes an outer approximation of S such that S ⊆ [x] ⊆ [x]. In the BPF, a general and well-known contraction method, the constraints propagation (CP) technique is applied [12]. Illustrative examples of CP method can be found in [11].

Target Tracking with Distributed Box Particle Filtering
Assume that there are static sensors and a moving target. These sensors are randomly distributed over a geometrical field with known positions and can periodically obtain measurements from the target. Our goal is to track the target in a distributed way so as to improve the tracking accuracy. We firstly derive the d-BPF algorithm and then give the Bayesian justification for it. Finally, some choices for the combination weight are given.

Distributed Box Particle
Filtering. Consider the following state-space model: where f : R → R is a nonlinear transition function from time to + 1 with being the dimension of the states. k +1 denotes an independent identically distributed (i.i.d.) process noise. Function g : R → R defines the relationship between the states and the measurements in node at time + 1, where is the dimension of the measurements. w , +1 is an i.i.d. measurement noise of node . Here, we consider interval-based uncertainties and thus interval measurements. That is, each sensor provides interval measurements [y , +1 ], which can be expressed as [y , +1 ] = [y , +1 + Δ 1 , y , +1 + Δ 2 ], where Δ 1 and Δ 2 denote the lower bound and the higher bound uncertainties, respectively. It is assumed that each node knows the locations of its neighbors at each time .
In [10,11], target tracking with a single sensor using BPF algorithm has been studied. Different from the PF, in the BPF, each particle is a box, which takes into account the imprecision caused by the model errors [11]. In this manner, the contraction step is similar to the variance matrix measurement update step of the Kalman filtering. Therefore, we follow a similar guideline of distributed Kalman filtering [26] to derive a distributed BPF algorithm for target tracking in a sensor network.
The d-BPF algorithm can be divided into five steps: initialization, time update and measurements broadcast, measurement update, estimation fusion, and resampling, where the second to the fourth steps are shown schematically in Figure 1.
(i) Initialization. At the beginning, the box particles in each node have to be initialized. Taking node as an example, it can be initialized by equally weighted and mutually disjoint boxes represented by ,0 ] denotes the box particles and ( ) ,0 denotes the corresponding weights. The BPF is superior to the PF in that it can cover larger prior uncertainty region with much fewer box particles.
(ii) Time Update and Measurements Broadcast. Referring to Figure 1(a), each node obtains a set of predicted box particles according to the transition function locally and broadcasts measurements to its close neighbors simultaneously. Knowing a set of weighted box particles {[x ( ) , ], ( ) , } =1 at time step and assuming that the system noise is known to be enclosed in the interval [k +1 ], the predicted boxes at step +1 can be calculated by It is worth mentioning that the inclusion function [f] needs to be used for ensuring that the predicted box particle is still an interval when propagating the box particles.
(iii) Measurement Update. As shown in Figure 1(b), after performing time update locally and broadcasting measurements to the neighbors, the weights of the predicted box particles in each node are then updated using the new available measurements from its neighbors (including itself) at time step + 1. This is an important step in the d-BPF, which is different from the BPF. For all box particles ( = 1 ⋅ ⋅ ⋅ ) in node , we have the predicted box measurements expressed as where [g ] is an inclusion function for g and N represents the neighbors of node . Likelihood factors are calculated using innovation quantities, which reflects the proximity between the measured box and the predicted box measurement [11,20]. The innovation is represented by the intersection between these two boxes, expressed as The likelihood factor is calculated based on the following idea: if predicted measurement of a box particle is included in the new measurement box, its updated weight should be larger. In contrast, if its predicted measurement does not have any intersection with the new measurement box, its updated weight should be equal to zero [11]. Furthermore, the contraction step is used to remove the intervals in the box particles that are not consistent with the measurements.

International Journal of Distributed Sensor Networks
Each node carries out time update step locally and exchanges observations with neighbors.
Each node carries out measurement update step locally using all the available observations and finally calculates the local state estimation.
Each node combines all the available local estimations together to get a new estimation.
Exchange local state estimation with neighbors.  Based on these concepts, for each node , the likelihood of a box particle at time + 1 can be expressed in the form where Note that the term [x ( ) , +1 ( )] in (16) represents the th dimension of updated th box particle after contraction using all the new available measurements that node receives. | ⋅ | denotes the length of the interval. The implementation of contraction step is given in Algorithm 1.
Algorithm 1 (contraction in measurement update). Consider the following.
International Journal of Distributed Sensor Networks 5 Based on the likelihood in (15), the update weights for each node at time + 1 can be calculated by and then the weights are normalized by the following: Hence, a set of updated weighted box particles , +1 } =1 are obtained. Meanwhile, local estimation in each node can be approximated by the weighted particles as where c ( ) , +1 is the center of the th box particlex ( ) , +1 at time step + 1 in node .
It is remarked that the box particles in each node are different after contraction, since the new available measurements are different in each node. So, it is improper to transmit weights to the neighbors as that performed in the distributed PF. Instead, here we let nodes only broadcast measurements to their closed neighbors in a distributed way and the box particles at each node are obtained based on the information in its neighborhood.
(iv) Estimation Fusion. In Figures 1(c) and 1(d), a combination step is applied and it attempts to achieve the new estimatê x , +1 that approximates the global performance via local node interactions; that is, where , is the ( , ) element of the combination matrix ∈ R × denoting the weight between and , which satisfies where 1 is an × 1 column vector with unity entries. There are several ways to select , , and this will be described elaborately later.
(v) Resampling. The resampling is performed to introduce variety while preventing the degeneration of the box particles. Similar to the BPF algorithm with a single sensor [10,11], the multinomial resampling combined with a subdivision step is applied, in which a box particle with high weight is divided by the corresponding number of smaller boxes, instead of obtaining the duplicated particles after the resampling step. The d-BPF algorithm implements these steps in each time slot and a series of weighted box particles and the state estimation can be calculated sequentially in each node. To summarize, its implementation procedure is given in Algorithm 2.
Algorithm 2 (distributed box particle filtering). Consider the following.
where ( ), Weight Update. Consider Weight Normalization. Consider Local State Estimate. Consider (4) Estimation Fusion. Consider (5) Resampling. Combine subdivision method for resampling to create new box particles with the same weights.
(6) Set = + 1, and go to step 2 until a predefined number end .

6
International Journal of Distributed Sensor Networks

Distributed Box Particle Filtering Analysis within Bayesian
Framework. In this section, similar to the proof of BPF [10,11], we present the Bayesian justification for d-BPF algorithm.
To simplify the representation, we denote the set of all the states and measurements obtained until time by X = {x , = 1, . . . , } and Y = {y , = 1, . . . , }, respectively. Our objective is to estimate the predicted distribution (x +1 | Y ) and the posterior distribution (x +1 | Y +1 ), recursively. Once the posterior pdf (x +1 | Y +1 ) is estimated, the system states can be calculated, correspondingly.
Within the framework of Bayesian inference, the posterior pdf of the state vector (x +1 | Y +1 ) is given by is a normalization factor. The recursion can be updated with an initial prior pdf (x 0 ), for example, a uniform distribution over some region of the state space. Equations (31) and (32) are known as the measurement update and the time update, respectively. Let a set of box particles be represented by a mixture of uniform pdfs. Comparing with the Gaussian mixture pdf, the uniform pdf is not only simple, but also owns a remarkable adaptation ability to strong nonlinear cases. Based on the interval analysis and by choosing box supports, d-BPF can propagate pdfs through linear, nonlinear, and even differential functions. More specifically, since the particles in the d-BPF algorithm are boxes and with the help of the interval analysis, the box particles can propagate sequentially with the uniform pdf.
Referring to [11], the uniform pdf sum representation of x can be given as where [x ( ) ] denotes the uniform pdf of the box particle [x ( ) ] and denotes the number of mixture components. According to (33), the relative pdfs in the Bayesian inference can be represented by a set of weighted box particles. Now we can deduce the Bayesian expression of the d-BPF algorithm. Assume that a set of weighted box particles , ] (x , ), ( ) , } =1 of node is obtained at time . Consequently, the pdf of the state at time can be represented by Then, the time update step at each node is carried out locally. Taking node as an example, according to (32), the pdf of the states can be predicted by Referring to [11], we can derive Based on (36), (35) can be rewritten as where the second equality holds according to (12). So far, we derive that the predicted pdf after time update in each node can still be represented by a set of weighted box particles, where the box particles change according to the transition function, while the weights keep invariant. In the measurement update step, each node spreads its observations to its neighbors. When nodes receive information from their neighbors, the measurement update step can be implemented. Consequently, according to (31), the posterior pdf of the states can be computed by where +1 = ∫ R (x , +1| | Y ∈N , ) (y ∈N , +1 | x , +1| ) x , +1| is a normalized factor. Referring to Algorithm 2, the measurement update step consists of two main steps, that is, contraction and weight update. The former is used to contract a set of new box particles and the latter is used to calculate updated weights.
International Journal of Distributed Sensor Networks 7 Once a set of new weighted box particles is obtained, the estimation can be calculated naturally. We first focus on the contraction step. Since the measurements are intervals, the measurement likelihood function (y , +1 | x , +1| ) is taken to be one component with uniform distribution, which is expressed as Consequently, (38) can be rewritten as We use [y ] to denote [y , +1 ] (g (x ( ) , +1| )) for short. Then, (40) can be rewritten as where ∈ N , = 1 ⋅ ⋅ ⋅ with being the degree of node . The term (41) is a function with a support being the following set: ∈ [y , +1 ]} . , +1| ] (x , +1| )}}} is also a support or a box particle. In Table 2: Several choices of combination matrix . In all cases, , = 0 if ∉ N , and , is chosen such that ∑ =1 , = 1 for all .
After contraction, the corresponding weights can be calculated according to (16). Finally, the posterior pdf at time + 1 can be expressed as which has the same form as (x , | Y ∈N , ) at time . That is to say, the posterior distribution (x , +1 | Y ∈N , +1 ) can be estimated recursively following the framework of Bayesian inference, and thus the feasibility of the d-BPF algorithm is justified.

Remark 3.
It is also remarked that the d-BPF is superior to the noncooperation BPF (Nc-BPF). On the one hand, it can be seen from (42) that the contraction uses the measurements to remove the inconsistent part of each box particle. As a result, combined with measurements of neighbors, the box particles after contraction in each node would have less or at least no more inconsistent part than those without cooperation; that is, the uncertainty of box particle is decreased. On the other hand, according to (16), in the d-BPF algorithm, we can deduce that the box particles with large consistent part have larger normalized weights while those with small consistent part acquire smaller normalized weights, comparing with those in a noncooperation algorithm. As a conclusion, a better set of normalized weights of the box particles can be obtained than that without cooperation and thus the estimate outperforms the noncooperation one.

Selection of Combination Weight.
The entities of combination matrix represent the weights used to combine the estimates of the neighbors. Some popular choices in the context of graph theory, average consensus, and diffusion adaptive networks are summarized in Table 2, where denotes the degree of node . In this paper, we employ the metropolis rule as it determines the weights only according to the degrees of the neighbors. According to [23], weights chosen by the metropolis rule can guarantee convergence provided that the graph is not bipartite.
Moreover, considering that the sensors are usually adopted with mobility, a dynamic network topology is also 8 International Journal of Distributed Sensor Networks considered. To model this, we assume that the links and nodes are random entities for simplicity. That is, at every time instant , the undirected random link weight , ( ) = , ( ), which connects node to node , is set to either a nominal value , = , with probability , = , or zero with probability , = 1 − , ; that is, the combination matrix A is selected by where , is selected according to the metropolis rule mentioned above.
Remark 4. Note that when the sensor network is fully connected, the combination matrix selected by the Metroplis rule becomes A = (1/ )I , where I is an identity × matrix. In this case, the estimate of each node is the same since every node collects the information of all the sensors. That is to say, each sensor can be considered as the fusion center in the centralized processing, and consequently the d-BPF reduces to the centralized BPF (c-BPF).

Simulation
In this section, we perform two examples of target tracking using the d-BPF algorithm proposed in Section 3. In these two examples, we consider the target tracking in a 2D plane. Assume there are = 50 sensors randomly deployed in the area of 300 m × 500 m. There is one target in the area and its state information at time is denoted by x = [ ,̇, ,̇] , where and represent the position of the target andȧ nḋdenote the velocity of the target, respectively. In the following, we will discuss these two examples, respectively.

Example 1.
In the first example, we consider a linear and Gaussian system that can be expressed as The target is moving according to the nearly constant velocity model with a transition density function (x +1 | x ) = N(x +1 , Fx , Q): where Q is the covariance matrix of the system noise, I 2 is 2 × 2 identity matrix, ⊗ is the Kronecker product, denotes the sampling interval (here set to 1 second), and denotes the process noise intensity (here set to 0.005). Initially, at = 0, the target is located at position (550, 300) m and is moving with the velocity (−5, −8.5) m/s. The target is tracked for = 50 s. A nonlinear and Gaussian measurement function g is also used. Each node uses the same measurement model expressed as where +1 ,̇+ 1 , +1 denote range, range rate, and azimuth, respectively. The measurement noise is ordinary additive zero-mean white Gaussian denoted by w , +1 , with a covariance matrix Σ = diag ( 2 , 2̇, 2 ), where = 1 m,̇= 0.01 m/s, and = 0.025 rad. The sensors provide interval measurements with interval length Δ = [Δ , Δ, Δ ] , where Δ = 20 m, Δ̇= 0.2 m/s, and Δ = 0.4 rad are the length of intervals in range, range rate, and azimuth, respectively. That is, the interval-based measurement of node at time + 1 lies within where Δ 1 = −(1/2)Δ and Δ 2 = (1/2)Δ denote the lower bound and the higher bound uncertainties, respectively.
In the following, both the simulations on static network topology and dynamic network topology are represented, respectively.

Case 1: Static Network
Topology. Firstly, we consider a fixed network topology as shown in Figure 2. The positions of these sensors are perfectly known and each sensor has the knowledge of the positions of its neighbors. The experiment is implemented using MATLAB and the so-called toolbox INTLAB [27]. The results are averaged over 50 experiments. To evaluate the performance of the proposed d-BPF, the rootmean-square error (RMSE) is used, which is defined as where Γ denotes the total number of the experiments and the termx , | represents the estimation of node at time in the th experiment.
In this example, we use 32 box particles for the d-BPF algorithm. Figure 3 shows tracking trajectory of the d-BPF of node 1 for a single run based on measurements generated from (47). Circle marks represent the true target trajectory, while the "plus" represents the estimated trajectory. From Figure 3, it is observed that the d-BPF is able to track the trajectory of target accurately. Note that the results for different sensors are similar, and thus they are omitted here. Figure 4 shows the averaged performance of the entire network over 50 Monte Carlo simulations. For performance comparison, the results of the Nc-BPF, the c-BPF, and distributed PF (d-PF) with 300 particles and 5000 particles are also presented. For clarity, the results are presented by taking the Napierian logarithm. It is observed that the d-BPF has significantly reduced the RMSE comparing to the Nc-BPF. Although it is a little inferior to the c-BPF from the viewpoint of RMSE, less communication cost is taken since the number of links is greatly smaller than that in the c-BPF.
The performance of the d-PF is largely dependent on the number of the particles. As shown in Figure 4, the d-BPF with only 32 box particles outperforms the d-PF with 300 particles. Although the accuracy of d-BPF is worse than the d-PF with 5000 particles, it takes much less computation time than the d-PF by using fewer box particles. It is also worth mentioning that, in order to make a comparison between d-BPF and d-PF algorithms, the noise we used in the simulation is Gaussian. It is a special case that we can use cumulative distribution function (cdf) to calculate weights so that the d-PF can deal with the interval measurements. This does not contradict with the aforementioned statement that the d-PF is in general unsuitable for the interval measurements, since the prior pdf of noise may be non-Gaussian or even unknown. However, the d-BPF can handle these cases very well. Therefore, it is where the d-BPF is expert.

Case 2: Dynamic Network Topology.
In the second case, we simulate the d-BPF over a dynamic sensor network. In the simulation, the propagation model is also the constant velocity model defined in (45). The measurement model defined in (47) is used. We also use 32 box particles for d-BPF algorithm. Assume that all the link probabilities between two nodes are equivalent and denoted as = , and keeps invariant in each simulation. Note that when is close to 0, the network almost reduces to a noncooperation one, that is, Nc-BPF. When = 1, the network is fully connected, and the result corresponds to the c-BPF. Figures 5(a) and 5(b) show the average position RMSE and average velocity RMSE, respectively, using different link probabilities over 50 Monte Carlo simulations, while Table 3 shows the corresponding average RMSE for ∈ [20, 50] s. From Figure 5, we can see that the performance of d-BPF over dynamic network depends on the network topology in each time slot, which is dependent on the value of the link probability . With the increasing of , the performance of the d-BPF is increased. Moreover, initially, when is relatively small, for example, < 0.2, the average RMSE is reduced greatly compared to the Nc-BPF. After that, with the increasing of , the performance of d-BPF improves slightly. When is larger than about 0.6, further increasing the value of , the performance of d-BPF does not improve much and converges to the c-BPF. It is noted that although the networks with larger link probabilities show better estimation performance, their communication cost is also increased. Therefore, in practice, we should make a tradeoff between the estimation performance and the cost.

Example 2.
In the second example, we consider the target tracking of a variable velocity model, expressed as So the target is moving with a transition density function (x +1 | x ) = N(x +1 , Fx + Gu, Q): where F and Q are defined in (46). The measurement function g is selected the same as (47). In the simulation, the same network as given in Figure 2 is used. The initialization and other parameters are set the same as those in Example 1. We also use 32 box particles for the d-BPF algorithm in this example. Figure 6 shows the tracking trajectory of the d-BPF of node 1 for a single run. From Figure 6, it is observed that the d-BPF is also successful in tracking the trajectory of the target even if the velocity of target is time-varying. Figure 7 shows the averaged performance of the entire network over 50 Monte Carlo simulations. The results of the Nc-BPF, the c-BPF, and distributed PF (d-PF) are also presented. Note that the results are also presented by taking the Napierian logarithm. Similar to Example 1, it is observed that the d-BPF still significantly reduces the RMSE comparing with the Nc-BPF in this case.

Conclusion
In this paper, based on the traditional box particle filtering and distributed processing, we have proposed a distributed box particle filter for target tracking, namely, d-BPF, to handle the cases where the measurements are intervals-based vectors. The feasibility of the proposed d-BPF algorithm is justified. The simulations for tracking both constant and variable velocity models are carried out. From the simulation results, we find that the tracking performance of our proposed d-BPF is significantly superior to that of the Nc-BPF. Moreover, from the simulation results with a dynamic network topology, we find that the tracking performance is improved with the increasing of the link probability. It is close to the c-BPF when the link probability is much larger, while the communication cost is increased simultaneously. Therefore, in practice, we should make a balance between the tracking performance and its cost.