# Distributed Particle Flow Filter for Target Tracking in Wireless Sensor Networks

Authors: Junjie Wang; Lingling Zhao; Xiaohong Su
Journal: Journal of Sensors
 (2018)

---

## Abstract

We propose, in this paper, a fully distributed tracking algorithm based on particle flow filter over sensor networks based on the max-consensus. The presented distributed particle flow filter is particularly suitable for the sensor network with limited sensing range and consists of two phases: the estimation phase and consensus phase. The local estimation results are obtained via particle flow filter in the estimation phase; then the sensor nodes agree on the best estimation based on max-consensus protocol in the consensus phase. Numerical simulations and comparisons with other distributed target tracking algorithms are carried out to show the effectiveness and feasibility of our approach.


---

## Body

## 1. Introduction

Distributed target tracking focuses on using a group of sensors to collect and process information about environment status. Compared with the central target tracking, the distributed target tracking has the following characteristics: scalability, flexibility, robustness, and fault tolerance. Due to these characteristics, the distributed target tracking has played an import role in many applications such as pedestrian tracking [1], biology [2], and environmental monitoring [3].The distributed target tracking algorithms can be classified into three types: fusion center (FC) based, leader agent (LA) based, and consensus based [4]. In the FC-based approaches, each sensor node uses the local measurement to estimate the local states by filtering algorithms and then transmits the local estimation to a single FC, where a global estimate is calculated based on all the local estimates. In a LA-based filter, only a subset of sensors are activated in a special manner and the information about the target is accumulated along a path formed by selected sensors. While in a consensus-based filter, all sensor nodes are simultaneously active and process the local data to get local posterior; then each one communicates with its neighboring agents [4] to agree with global posterior using consensus algorithms.The Kalman filter is an optimal target tracking algorithm in the linear Gaussian situation. Recently, the distributed Kalman filter (DKF) [5] for the track-to-track fusion has been proposed as an optimal solution. However, the track-to-track DKF needs to handle the multiple information paths. Another well-known distributed Kalman filter is the Kalman consensus filter (KCF) [6] based on the consensus fusion algorithm which can avoid addressing the multiple information paths. The KCF performs well when all the nodes can get the measurement of a target. In the realistic scenario with limited sensing capability of sensors, some nodes become naive about the target state at some time instants [7]; the performance of the KCF will deteriorate as each node weighs its neighbors’ estimates in an equal manner. To overcome this issue, the generalized Kalman consensus filter (GKCF) [7] was proposed utilizing the weighted averaging consensus. A node selection [8, 9] strategy was proposed to select the best estimate and propagate it to all the nodes rather than fusing all the sensor-estimated results. The distributed Kalman filter with node selection [9] selects the most accurate estimate to propagate through the network rather than fusing all the sensor-estimated results. The distributed Kalman filters have the desirable property of computational simplicity in linear systems, but it is still needed to develop distributed target tracking algorithms for nonlinear non-Gaussian systems.The distributed particle filters (DPFs) [4, 10–12] have been proposed to track a target in a wireless sensor network for the nonlinear non-Gaussian system. In [11], a distributed particle filter computes the product of likelihood function over the network using iterative average consensus. A kind of information-weighted consensus-based distributed particle filter [12] can avoid the divergence of the consensus error introduced by the naive nodes, but it performs in a low convergence rate. In the distributed particle filter proposed in [13, 14], the local posterior probability density function (PDF) is approximated by a Gaussian distribution. Then the local PDF parameters are fused into the global posterior PDF’s parameters by average consensus. The work in [15] approximates the local likelihood functions by a Gaussian function and builds the global likelihood through exchange of information with neighboring nodes.The DPFs based on averaging consensus mentioned above have two defects. The first drawback is that they require a large number of particles or samples for a given level of accuracy. To address this drawback of DPFs, we adapt another nonlinear filter, namely, particle flow filter proposed by Daum and Huang [16–19] to estimate the target state. The particle flow filter can achieve a good performance with fewer particles compared with the particle filter, especially in the high-dimensional state space. The principle behind the particle flow filter is to sample a set of particles from the prior distribution and use a stochastic method to move them such that they are then distributed according to the posterior. In other words, particles are migrated smoothly using a particle flow derived from a log homotopy relating the prior and the posterior [20]. Compared with the particle filter, the particle flow filter can yield a significant reduction of the number of particles especially in the high-dimensional case. Another issue of DPFs based on averaged consensus is that they are not suitable for the scenario which there exist some naive nodes in the wireless sensor network (WSN). In this scenario, it may happen that only a minority of sensors have measurements. Therefore, the WSN involves many unreliable sensors in the fusion step which may cause a divergent error. This is because of the fact that the average consensus algorithm gives all the nodes equal weights; even the naive nodes get less information about the target.In this paper, inspired by the particle flow filter framework in [17] and the fusion rule in [21], we propose a distributed particle flow filter (DPFF) algorithm for WSN. We approximate local posterior as a Gaussian distribution and fuse the local posterior via a max-consensus protocol. To address the challenge in limited energy and sensor range of sensor nodes, the particle flow filter [17] is utilized to approximate local posterior. Also, the proposed DPFF seeks consensus on the best local posterior, rather than on the average of local posteriori. To the best of our knowledge, particle flow filter has not been yet investigated in WSNs.The rest of the paper is organized as follows. Section2 introduces the background of the sensor network model and consensus theory. The details of DPFF are presented in Section 3. Section 4 evaluates the DPFF performance via two numerical examples. Finally, Section 5 gives the conclusion of this paper and the possible future work.


## 2. Background

### 2.1. Network Model

We consider a single-target tracking in the WSN which consists ofNc sensors with restricted monitoring area and communication ranges. Assume that the WSN can be modeled as an undirected connected graph G=C,ℰ, in which each vertex c∈C represents a sensor node and each edge i,j∈ℰ denotes the link between different nodes. If an edge i,j∈ℰ denotes that the node j can get information from the node i and vice versa, the set of neighbors of node i is denoted as Ni=j∈C:j,i∈ℰ.Figure1 shows a network with four nodes. The distributed WSN has no central unit, and, thus, the sensors locally process their measurements. We restrict the monitoring region of each sensor node defined as an area within a dotted circle of radius ρ, in the sense that a sensor may not detect the target over the tracking period. Also, we assume that each sensor can only directly communicate with its neighbors in a certain communication range.Figure 1
Abstract representation of the distributed tracking scenario. Four nodessa,sb,sc,andsd are randomly placed in the tracking environment. Each node has a local sensing area (dotted circles around the node). The target trajectory is represented by the dashed line. The communication topology is denoted by the solid black line between nodes.


### 2.2. Target Model

The state of a target is represented by a vectorx. For a maneuvering target, the state vector x contains the information about the position and velocity. The dynamic transition of a target is modeled as
(1)xk=fxk−1+wk,where f is the state transition function and wk is the processing noise.At timek, each sensor can only obtain measurement when a target appears in its sensing area. The measurement function of sensor i is
(2)zki=hkixk+vki,where hki denotes the observation function for the ith sensor node and vki stands for the measurement noise. The measurements of nodes are assumed independent over the network.


### 2.3. Max-Consensus

The max-consensus [21] is a well-known distributed algorithm which makes all the sensor nodes agree on the maximum of the value of their initial state through finite iterations. In the max-consensus algorithm, each sensor node initializes its state value as ui0=ui and iteratively communicates with its neighbors Ni based on the update rule as follows to update its state:
(3)uit+1=maxj∈Ni∪iujt,where ui is the state value of node i and t is the iterative step index.According to [21], we define the max-consensus as follows:Definition 1 (max-consensus).
Consider a WSN withC nodes, connected over an undirected graph G=C,ℰ. Each node has a state variable ui0,i∈C. The discrete time max-consensus protocol is defined as
(4)uix=ujk=maxu10,…,un0∀k≥K,∀i,j∈C.
If (4) holds for all u0∈ℝn, strong max-consensus is achieved. If (4) only holds for a subset of initial states, weak max-consensus is achieved [22].


### 2.4. Average Consensus-Based Distributed Particle Filter in WSNs

The average consensus-based distributed particle filter consists two steps: local particle filter and average consensus filter. The local particle filter uses the local observation to get the local estimation. The output of the local particle filter is the local posteriori approximated as a Gaussian distributionN·∣mki,local,Pki,local. For the sake of using outputs of local filter among networks more effectively, each node maintains an average consensus filter. The aim of the average consensus filter is to fuse these local posteriori between neighbor nodes. According to [13], node i at time k, based on the local measurement zki, runs a local particle filter to obtain the parameters of local posterior mki,local,Pki,local as follows:
(5)mki,local=∑j=1Mwkijxkij,Pki,local=∑j=1Mwkijxkij−mki,localxkij−mki,localT,where wkij and xkij denote the weights and particles of the local particle filter, respectively. Then the average consensus algorithm is run between linked neighbor nodes as the following equation to obtain the global estimation results:
(6)yi,l+1=yi,l+ϵ∑j∈Niyj,l−ji,l+ui,l−yi,l,where ϵ is the updating rate. yi,l can be represented as a parameter of local posterior such as mki,local or Pki,local.


## 2.1. Network Model

We consider a single-target tracking in the WSN which consists ofNc sensors with restricted monitoring area and communication ranges. Assume that the WSN can be modeled as an undirected connected graph G=C,ℰ, in which each vertex c∈C represents a sensor node and each edge i,j∈ℰ denotes the link between different nodes. If an edge i,j∈ℰ denotes that the node j can get information from the node i and vice versa, the set of neighbors of node i is denoted as Ni=j∈C:j,i∈ℰ.Figure1 shows a network with four nodes. The distributed WSN has no central unit, and, thus, the sensors locally process their measurements. We restrict the monitoring region of each sensor node defined as an area within a dotted circle of radius ρ, in the sense that a sensor may not detect the target over the tracking period. Also, we assume that each sensor can only directly communicate with its neighbors in a certain communication range.Figure 1
Abstract representation of the distributed tracking scenario. Four nodessa,sb,sc,andsd are randomly placed in the tracking environment. Each node has a local sensing area (dotted circles around the node). The target trajectory is represented by the dashed line. The communication topology is denoted by the solid black line between nodes.


## 2.2. Target Model

The state of a target is represented by a vectorx. For a maneuvering target, the state vector x contains the information about the position and velocity. The dynamic transition of a target is modeled as
(1)xk=fxk−1+wk,where f is the state transition function and wk is the processing noise.At timek, each sensor can only obtain measurement when a target appears in its sensing area. The measurement function of sensor i is
(2)zki=hkixk+vki,where hki denotes the observation function for the ith sensor node and vki stands for the measurement noise. The measurements of nodes are assumed independent over the network.


## 2.3. Max-Consensus

The max-consensus [21] is a well-known distributed algorithm which makes all the sensor nodes agree on the maximum of the value of their initial state through finite iterations. In the max-consensus algorithm, each sensor node initializes its state value as ui0=ui and iteratively communicates with its neighbors Ni based on the update rule as follows to update its state:
(3)uit+1=maxj∈Ni∪iujt,where ui is the state value of node i and t is the iterative step index.According to [21], we define the max-consensus as follows:Definition 1 (max-consensus).
Consider a WSN withC nodes, connected over an undirected graph G=C,ℰ. Each node has a state variable ui0,i∈C. The discrete time max-consensus protocol is defined as
(4)uix=ujk=maxu10,…,un0∀k≥K,∀i,j∈C.
If (4) holds for all u0∈ℝn, strong max-consensus is achieved. If (4) only holds for a subset of initial states, weak max-consensus is achieved [22].


## 2.4. Average Consensus-Based Distributed Particle Filter in WSNs

The average consensus-based distributed particle filter consists two steps: local particle filter and average consensus filter. The local particle filter uses the local observation to get the local estimation. The output of the local particle filter is the local posteriori approximated as a Gaussian distributionN·∣mki,local,Pki,local. For the sake of using outputs of local filter among networks more effectively, each node maintains an average consensus filter. The aim of the average consensus filter is to fuse these local posteriori between neighbor nodes. According to [13], node i at time k, based on the local measurement zki, runs a local particle filter to obtain the parameters of local posterior mki,local,Pki,local as follows:
(5)mki,local=∑j=1Mwkijxkij,Pki,local=∑j=1Mwkijxkij−mki,localxkij−mki,localT,where wkij and xkij denote the weights and particles of the local particle filter, respectively. Then the average consensus algorithm is run between linked neighbor nodes as the following equation to obtain the global estimation results:
(6)yi,l+1=yi,l+ϵ∑j∈Niyj,l−ji,l+ui,l−yi,l,where ϵ is the updating rate. yi,l can be represented as a parameter of local posterior such as mki,local or Pki,local.


## 3. Distributed Particle Flow Filter

The average consensus-based distributed particle filter is robust to time-varying network topologies [4]. In the average-based distributed particle filter, each node runs a local particle filter to estimate the target state. This method inherits the drawbacks of the particle filter, namely, the curse of dimensionality, and particle degeneracy in the highly informative scenario. The other drawback of the particle filter is its requirement to maintain a large number of particles to attain good performance, leading to the waste of the sensors’ energy.Motivated by some desirable properties of the particle flow filter (see for example, [23, 24]) such as sufficient accuracy and low computational complexity, uniqueness of the solutions, we propose a consensus-based distributed particle flow filter (DPFF) algorithm, which consists of two main phases: estimation and consensus. According to the sensing range of sensors, the estimation phase can be divided into two branches. If a target is measured by the node i at time k, its estimation is carried out via a particle flow filter. On the contrary, the node can not detect the target at time k; only the prediction part of the particle flow filter is run. Then each sensor computes its perception confidence value γik (refer to (15)) based on the estimated posterior covariance matrix or prior covariance matrix. In the consensus phase, a max-consensus algorithm is utilized to make all the sensors agree on the best estimated sensor.Note that our algorithm requires the synchronization of clocks over the sensor networks. We now give the details of the two phases of DPFF at timek and assumed that each node obtains the best estimation results μ¯k−1,P¯k−1 over the WSN at time k−1.


### 3.1. Estimation Phase

In the DPFF, the posterior distribution is approximated by a particle setxk−1i,wk−1ii=1Np which is sampled from Gaussian distribution N·,μ¯k−1,P¯k−1 with same weights wk−1i=1/Np. Then all particles are transmitted to x^kii=1Np at time step k by the dynamic model. Therefore, the prior distribution can be represented by particles x^kii=1Np.The overall process of the estimation phase for each node is outlined in Figure2. As a target may move in or out of the sensing area of the node i, if the sensor i has detected the target in its sensing range, then its estimation is carried out by the particle flow filter (Algorithm 1). The local output of each node is the local posterior approximated as a Gaussian distribution. On the contrary, in the case that the measurement is not available to the node i, then only the prediction of the particle flow filter will be executed (lines 1–3 of Algorithm 1).Algorithm 1: Particle flow filter for node i at time k.
Input:μ¯ik−1,P¯ik−1 Output:μ^ik,P^ik
1: redraw particlesx^ijkj=1Np from the Gaussian
distribution whichμ¯ik−1,P¯ik−1 are the mean and
covariance, respectively.
2:forj=1,…,Npdo
3:   propagate particlesxijk=fkx^ijk+vk
4: computePi−k=covxk
5:forj=1,…,Nλdo
6:   Calculate the mean particles valuex¯k
7:   setλ=jΔλ
8:   linearize the measurement function atx¯k to get
Hxk
9:   CalculateA and b from (12) and (13), respectively, using
Pi−k,x¯k,Hxk;
10:forj=1,…,Npdo
11:   migrate particlesxijk=xijk+ΔλAxijk+b
12: Compute the mean value and error matrixμ^ik,P^ik of the particles by (14)Figure 2
The block diagram of the estimation phase for each node.In the following, the particle flow filter is presented in details. To avoid confusion, we will omit particle state indices. Particle flow filter is used to guide particles by the current measurement so that they can more accurately approximate the posterior distribution. Lethxk=pyk∣xk and gxk=pxk∣y1:k−1 denote the likelihood and prior functions, respectively. The log homotopy is defined as
(7)logpx,λ=loggxk+λloghx−logKλ,where λ is a real number that varies from 0 to 1. Kλ is the normalization constant. When λ=0, we obtain the predicted density function gxk, and when λ=1, we get the posterior distribution. Suppose that the flow of particles is guided by the Bayes’ rule according to
(8)dx=fx,λ+dw.Combined by (7), fx,λ can be computed by the Fokker-Planck equation
(9)∂px,λ∂λ=−∇·fx,λ,px,λ+12trQ∂2px,λ∂x2,where Q is the covariance of the process noise. For simplicity, it is assumed that Q=0. According to (7) and (9), the following equation can be derived:
(10)loghx+∇logpx,λT·fx,λ=−∇fx,λ.Assuming thatgx,hx are Gaussian PDFs, then a closed-form solution termed as exact flow filter [25] can be derived as
(11)dxdλ=Aλx+bλ,where
(12)A=−12PHTλHPHT+R−1H,(13)b=I+2λAI+λAPHTR−1z+Ax¯.Letx¯ and P represent the predicted mean vector and the prior covariance matrix, respectively. H denotes the measurement function matrix and R is the covariance of the measurement noise. For nonlinear models, the measurement function matrix H can be obtained by linearization of the measurement model. For more details on the implementation and analysis of the exact flow filter, please refer to [20, 25]. We summarize the exact flow filter in Algorithm 1.It is important to note that there exist several different realizations for the particle flow filter, such as nonzero diffusion particle flow filter [26], incompressible flow filter [16], and exact Daum and Huang (EDH) filter [18]. In some cases, the Gromov’s method was explained in [27–29], which would improve the particle flow filter performance. The performance of the particle flow filter also is influenced by the discretization of pseudotime λ. As the particle flow filter is described by an ordinary differential equation, a suitable discretization is essential to capture the flow dynamics [30]. In this paper, we set a sequence of discrete steps with uniform step size Δλ=0.1.In summary, when the sensor node received the measurement, the node can obtain the estimation resultμ^ik,P^ik according to the particle flow filter (Algorithm 1). On the other hand, if the target is not detected by the ith sensor node, the target state is estimated according to the particles x^ki,jj=1Np from prediction. Then the estimated state and error covariance matrix can be computed as
(14)μ^ik=1Np∑j=1Npx^ki,j,P^ik=1Np∑j=1Npx^ki,j−μ^ikx^ki,j−μ^ikT.In order to measure the quality of the estimated target state of each node, the perception confidence valueγik needs to be calculated at the end of the estimation phase. If the target is detected by the sensor i, γik is calculated based on the posterior error covariance matrix P^ik; otherwise, γik is calculated based on the prior error covariance matrix P^ik. The specific calculation formula of γik is as follows
(15)γik=1traceP^ik,where trace· is the matrix trace operator. It is clear that γik grows with the reliability of the estimation performed by sensor node i at time k [22].At the end of the estimation phase, each sensor node will obtain the value of perception confidence valueγik, the local estimate μ^ik, and the local error covariance matrix P^ik. These values can be employed to reach consensus in the consensus phase.


### 3.2. Consensus Phase

The aim of the consensus phase is to select the best estimation over the sensor networks and propagate the selected estimation with correlation information (error covariance matrix). It is noted that, in order to reduce the communication cost, we exchange the error covariance matrix and state estimate rather than the whole particle set. Therefore, at the next time step, each sensor node needs to redraw particles based on the best state estimate and error covariance matrix. The max-consensus for nodei at time k is reported in Algorithm 2.Algorithm 2: Max-consensus for node i at time k.
Input:γik,μ^ik,P^ik Output:μ¯ik,P¯ik
1:Γi0=γik
2:Ui0=μ^ik
3:∏i0=P^ik
4:fort=1,…,Ddo
5:   sendΓit−1,Uit−1,∏it−1
6:   receive the information with connected sensors
setsNi to obtain Γjt−1,Ujt−1,∏jt−1, ∀j∈
Ni∪i
7:Γit=maxj∈Ni∪iΓjt−1
8:α=argmaxj∈Ni∪iΓjt−1
9:Uit=Uαt−1
10:∏it=∏αt−1
11:μ¯ik=UiD
12:P¯ik=∏iDAlgorithm2 works as follows: node i obtains the values from the estimation phase, γik,μ^ik,P^ik. And then sensor node i initializes its variables Γi0,Ui0,∏i0 with γik,μ^ik,P^ik, respectively. After initializing all variables, node i exchanges the variables Γit−1,Uit−1,∏it−1 with its neighbors (lines 5-6 of Algorithm 2). Then node i selects the max perception confidence value which corresponds to the best estimate from its neighbors (lines 7-8 of Algorithm 2). The node i will replace its estimation results to the corresponding estimate and covariance estimate matrices Uit and ∏it (lines 9 and 10, resp.). At the end of the phase, each sensor will agree with the best state estimate of the target μ¯ik with related error covariance matrix P¯ik(lines 12-13 of Algorithm 2). These two variables will be used in the next time k+1 of the estimation phase, in order to let the particle flow filter of each node start from the best estimation results and therefore to improve the algorithm’s performance.It can be proved [21] that the node converges during D steps, where D is the diameter of the sensor graph.


### 3.3. Convergence Analysis

In the DPFF algorithm, each node will converge within finite discrete steps through the max-consensus algorithm. The convergence of the DPFF is proved as follows:Assuming that each node runs the DPFF algorithm, after consensus phases at timek, node i and node j hold the equation
(16)μ¯ik=μ¯jk=μ¯maxk,∀i,j∈C,where C is the set of nodes and μ¯maxk is estimated mean value corresponding to the max perception confidence value γk.In each iteration of Algorithm2, line 7 is to select the max γk from the neighbor nodes which is the update rule for the max-consensus algorithm. The max-consensus algorithm is guaranteed to converge in a finite number of iterations. Therefore, in D steps, Algorithm 2 will be guaranteed that
(17)Γin=Γjn=maxl∈ℐγlk,∀i,j∈C.In each iteration, the variablesUi and ∏i store the corresponding mean value and error covariance matrix, respectively. Thus, at the end of iterations, (16) is workable. The convergence issues are discussed in [9] in detail.Remark 1.
In this paper, we approximate the posterior distribution with a Gaussian distribution. In the non-Gaussian system, the posterior distribution can be represented by the mixture of multiple Gaussian distributions. So, the DPFF can not be applied directly. A new particle flow filter called Gaussian mixture particle flow [31] can be employed to cater for the non-Gaussian situation. When the posterior distribution is represented by the Gaussian mixture model (GMM) ∑j=1NN·∣mkj,Pkj, the corresponding γk is calculated by the following:
(18)γk=1∑j=1NtracePkj.


## 3.1. Estimation Phase

In the DPFF, the posterior distribution is approximated by a particle setxk−1i,wk−1ii=1Np which is sampled from Gaussian distribution N·,μ¯k−1,P¯k−1 with same weights wk−1i=1/Np. Then all particles are transmitted to x^kii=1Np at time step k by the dynamic model. Therefore, the prior distribution can be represented by particles x^kii=1Np.The overall process of the estimation phase for each node is outlined in Figure2. As a target may move in or out of the sensing area of the node i, if the sensor i has detected the target in its sensing range, then its estimation is carried out by the particle flow filter (Algorithm 1). The local output of each node is the local posterior approximated as a Gaussian distribution. On the contrary, in the case that the measurement is not available to the node i, then only the prediction of the particle flow filter will be executed (lines 1–3 of Algorithm 1).Algorithm 1: Particle flow filter for node i at time k.
Input:μ¯ik−1,P¯ik−1 Output:μ^ik,P^ik
1: redraw particlesx^ijkj=1Np from the Gaussian
distribution whichμ¯ik−1,P¯ik−1 are the mean and
covariance, respectively.
2:forj=1,…,Npdo
3:   propagate particlesxijk=fkx^ijk+vk
4: computePi−k=covxk
5:forj=1,…,Nλdo
6:   Calculate the mean particles valuex¯k
7:   setλ=jΔλ
8:   linearize the measurement function atx¯k to get
Hxk
9:   CalculateA and b from (12) and (13), respectively, using
Pi−k,x¯k,Hxk;
10:forj=1,…,Npdo
11:   migrate particlesxijk=xijk+ΔλAxijk+b
12: Compute the mean value and error matrixμ^ik,P^ik of the particles by (14)Figure 2
The block diagram of the estimation phase for each node.In the following, the particle flow filter is presented in details. To avoid confusion, we will omit particle state indices. Particle flow filter is used to guide particles by the current measurement so that they can more accurately approximate the posterior distribution. Lethxk=pyk∣xk and gxk=pxk∣y1:k−1 denote the likelihood and prior functions, respectively. The log homotopy is defined as
(7)logpx,λ=loggxk+λloghx−logKλ,where λ is a real number that varies from 0 to 1. Kλ is the normalization constant. When λ=0, we obtain the predicted density function gxk, and when λ=1, we get the posterior distribution. Suppose that the flow of particles is guided by the Bayes’ rule according to
(8)dx=fx,λ+dw.Combined by (7), fx,λ can be computed by the Fokker-Planck equation
(9)∂px,λ∂λ=−∇·fx,λ,px,λ+12trQ∂2px,λ∂x2,where Q is the covariance of the process noise. For simplicity, it is assumed that Q=0. According to (7) and (9), the following equation can be derived:
(10)loghx+∇logpx,λT·fx,λ=−∇fx,λ.Assuming thatgx,hx are Gaussian PDFs, then a closed-form solution termed as exact flow filter [25] can be derived as
(11)dxdλ=Aλx+bλ,where
(12)A=−12PHTλHPHT+R−1H,(13)b=I+2λAI+λAPHTR−1z+Ax¯.Letx¯ and P represent the predicted mean vector and the prior covariance matrix, respectively. H denotes the measurement function matrix and R is the covariance of the measurement noise. For nonlinear models, the measurement function matrix H can be obtained by linearization of the measurement model. For more details on the implementation and analysis of the exact flow filter, please refer to [20, 25]. We summarize the exact flow filter in Algorithm 1.It is important to note that there exist several different realizations for the particle flow filter, such as nonzero diffusion particle flow filter [26], incompressible flow filter [16], and exact Daum and Huang (EDH) filter [18]. In some cases, the Gromov’s method was explained in [27–29], which would improve the particle flow filter performance. The performance of the particle flow filter also is influenced by the discretization of pseudotime λ. As the particle flow filter is described by an ordinary differential equation, a suitable discretization is essential to capture the flow dynamics [30]. In this paper, we set a sequence of discrete steps with uniform step size Δλ=0.1.In summary, when the sensor node received the measurement, the node can obtain the estimation resultμ^ik,P^ik according to the particle flow filter (Algorithm 1). On the other hand, if the target is not detected by the ith sensor node, the target state is estimated according to the particles x^ki,jj=1Np from prediction. Then the estimated state and error covariance matrix can be computed as
(14)μ^ik=1Np∑j=1Npx^ki,j,P^ik=1Np∑j=1Npx^ki,j−μ^ikx^ki,j−μ^ikT.In order to measure the quality of the estimated target state of each node, the perception confidence valueγik needs to be calculated at the end of the estimation phase. If the target is detected by the sensor i, γik is calculated based on the posterior error covariance matrix P^ik; otherwise, γik is calculated based on the prior error covariance matrix P^ik. The specific calculation formula of γik is as follows
(15)γik=1traceP^ik,where trace· is the matrix trace operator. It is clear that γik grows with the reliability of the estimation performed by sensor node i at time k [22].At the end of the estimation phase, each sensor node will obtain the value of perception confidence valueγik, the local estimate μ^ik, and the local error covariance matrix P^ik. These values can be employed to reach consensus in the consensus phase.


## 3.2. Consensus Phase

The aim of the consensus phase is to select the best estimation over the sensor networks and propagate the selected estimation with correlation information (error covariance matrix). It is noted that, in order to reduce the communication cost, we exchange the error covariance matrix and state estimate rather than the whole particle set. Therefore, at the next time step, each sensor node needs to redraw particles based on the best state estimate and error covariance matrix. The max-consensus for nodei at time k is reported in Algorithm 2.Algorithm 2: Max-consensus for node i at time k.
Input:γik,μ^ik,P^ik Output:μ¯ik,P¯ik
1:Γi0=γik
2:Ui0=μ^ik
3:∏i0=P^ik
4:fort=1,…,Ddo
5:   sendΓit−1,Uit−1,∏it−1
6:   receive the information with connected sensors
setsNi to obtain Γjt−1,Ujt−1,∏jt−1, ∀j∈
Ni∪i
7:Γit=maxj∈Ni∪iΓjt−1
8:α=argmaxj∈Ni∪iΓjt−1
9:Uit=Uαt−1
10:∏it=∏αt−1
11:μ¯ik=UiD
12:P¯ik=∏iDAlgorithm2 works as follows: node i obtains the values from the estimation phase, γik,μ^ik,P^ik. And then sensor node i initializes its variables Γi0,Ui0,∏i0 with γik,μ^ik,P^ik, respectively. After initializing all variables, node i exchanges the variables Γit−1,Uit−1,∏it−1 with its neighbors (lines 5-6 of Algorithm 2). Then node i selects the max perception confidence value which corresponds to the best estimate from its neighbors (lines 7-8 of Algorithm 2). The node i will replace its estimation results to the corresponding estimate and covariance estimate matrices Uit and ∏it (lines 9 and 10, resp.). At the end of the phase, each sensor will agree with the best state estimate of the target μ¯ik with related error covariance matrix P¯ik(lines 12-13 of Algorithm 2). These two variables will be used in the next time k+1 of the estimation phase, in order to let the particle flow filter of each node start from the best estimation results and therefore to improve the algorithm’s performance.It can be proved [21] that the node converges during D steps, where D is the diameter of the sensor graph.


## 3.3. Convergence Analysis

In the DPFF algorithm, each node will converge within finite discrete steps through the max-consensus algorithm. The convergence of the DPFF is proved as follows:Assuming that each node runs the DPFF algorithm, after consensus phases at timek, node i and node j hold the equation
(16)μ¯ik=μ¯jk=μ¯maxk,∀i,j∈C,where C is the set of nodes and μ¯maxk is estimated mean value corresponding to the max perception confidence value γk.In each iteration of Algorithm2, line 7 is to select the max γk from the neighbor nodes which is the update rule for the max-consensus algorithm. The max-consensus algorithm is guaranteed to converge in a finite number of iterations. Therefore, in D steps, Algorithm 2 will be guaranteed that
(17)Γin=Γjn=maxl∈ℐγlk,∀i,j∈C.In each iteration, the variablesUi and ∏i store the corresponding mean value and error covariance matrix, respectively. Thus, at the end of iterations, (16) is workable. The convergence issues are discussed in [9] in detail.Remark 1.
In this paper, we approximate the posterior distribution with a Gaussian distribution. In the non-Gaussian system, the posterior distribution can be represented by the mixture of multiple Gaussian distributions. So, the DPFF can not be applied directly. A new particle flow filter called Gaussian mixture particle flow [31] can be employed to cater for the non-Gaussian situation. When the posterior distribution is represented by the Gaussian mixture model (GMM) ∑j=1NN·∣mkj,Pkj, the corresponding γk is calculated by the following:
(18)γk=1∑j=1NtracePkj.


## 4. Experiments

In this section, we evaluate the performance of the proposed DPFF algorithm in the simulated environment and compare it with other approaches: the centralized particle filter (CPF) where we use the performance of the CPF as the base performance, the distributed particle filter based on average consensus [15] (DPF-AV), and the information weight average consensus-based distributed particle filter (DPF-WAV) [12].


### 4.1. Example 1: Grid Network

Consider such a simulation that a target moves in a 300 m × 500 m area with 15 sensors. The sensors have overlapped monitored regions. The monitored region of each sensor node is assumed to be a circle region of 75 m radius whose center is at the sensor’s location. A sensor has a measurement of the target only if the target appears in the sensor’s sensing area. Figure3 illustrates the sensors and network connectivity. Connections between the sensors are shown as grey dashed lines. Each sensor can only communicate directly with its neighboring sensors whose distance to it is less than the communication range. The state vector of the target is represented as xk=xk,yk,vx,k,vy,k,wk, where xk,yk and vx,k,vy,k represent the target position and velocity, respectively. wk denotes the turn rate of the target. In this case, a nearly coordinated turn model with the known constant turn rate wk and the unknown velocity vk is considered. This model is able to account for the motion of complicated maneuverable targets. The nonlinear scenario is used in [32], and the motion of target is modeled according to
(19)xk+1=Fωkxk+Gwk,ωk+1=ωk+Δuk,where
(20)Fω=1sinωΔω0−1−cosωΔω0cosωΔ0−sinωΔ01−cosωΔω1sinωΔω0sinωΔ0cosωΔ,Gω=Δ220Δ00Δ220Δ.wk~N·;0,σw2I,uk~N·;0,σu2I,Δ=1 s, σw=15 m/s2, and σu=π/180 rad/s. The observation is a noisy bearing and range vector given by
(21)zk=arctanpx,kpy,kpx,k2+py,k2+εk,where εk~N·;0,Rk, with Rk=diagσθ2,σr2T and σθ=2π/180,σr=10. px,py are the distance from the target to the sensor.Figure 3
Object trajectory, sensor node, and connectivity considered in the experiments. The connectivity among the nodes is shown using the grey dashed lines between nodes.The number of particles for each node is to set 500 in the DPF-WAV and DPF-AV algorithms, while the DPFF algorithm only uses 50 particles in each node, while the CPF algorithm collects the entire available measurements from all sensor nodes and uses 500 particles to estimate.The error metric that we have computed is the root mean square error (RMSE) between the true and estimated target positions at each time instantk from all sensors. Let pkx,pky and p^kx,p^ky denote the true and estimated target positions, respectively, at time k. The RMSE value at time k is calculated over a number of Monte Carlo NMC runs according to
(22)RMSE=1NMC∑i=1NMCpkx−p^kx2+pky−p^ky2,where NMC=500 is the number of Monte Carlo runs.A sample of the estimated target track of DPFF is shown in Figure4, which also shows the estimated track from the CPF, DPF-WAV, and DPF-AV algorithms. It is obvious that the DPF-AV lost the track while other algorithms work well when tracking the target.Figure 4
One example of estimated trajectory using DPFF (black curve), CPF (green curve), DPF-AV (cyan curve), and DPF-WAV (blue curve). The red curve indicates the ground truth.Figure5 shows the temporal evolution of RMSE. It is observed that the DPF-AV algorithm performs worst in these methods, while results of CPF, DPFF, and DPF-WAV are fairly close to each other. The performance of DPFF and DPF-WAV is almost as good as the CPF and better than DPF-AV. Particularly, there is serious divergence with DPF-AV. For clear comparison of the proposed algorithm with DPF-WAV, we plot only DPF-WAV, CPF, and DPFF in Figure 6. Figure 6 illustrates the relationship of RMSE with timestamp of the DPF-WAV, CPF, and DPFF algorithms. As the DPF-WAV has a sharp fluctuation, it is obvious that the DPFF is better than the DPF-WAV and closer to the CPF. Especially, the DPF-WAV needs more particles than the DPFF and usually requires a significant amount of computing resources.Figure 5
Root mean square error (RMSE) versus time.Figure 6
Root mean square error (RMSE) versus time.Moreover, we also computed the averaged RMSE (ARMSE) for all methods. The ARMSE is given by(23)ARMSE=1TNMC∑i=1NMC∑j=1Tpkx−p^kx2+pky−p^ky2.Table1 lists the ARMSE, standard deviation of RMSE values, and the average execution time per node over 500 simulation trials. The first column in Table 1 is the number of particles for different algorithms. It can be seen that the proposed DPFF only with 50 particles gets better estimate results than the DPF-WAV and DPF-AV, almost equal to the CPF with one-tenth the number of particles. For the average execution time per node, the DPFF is less than a quarter of the DPF-WAV and DPF-AV. In summary, the DPFF has higher accuracy and less execution time.Table 1
Average results over 500 runs with different methods.
Number
Methods
Average RMSE
Std RMSE
Average runtime (sec)
500
CPF
6.33
4.43
0.3919
DPF-WAV
9.16
4.35
0.2851
DPF-AV
126.62
27.51
0.2705
1000
CPF
5.98
4.03
0.7815
DPF-WAV
8.97
4.12
0.3334
DPF-AV
124.12
21.48
0.3102
50
DPFF
6.32
4.76
0.0672


### 4.2. Example 2: Random Network

In the second scenario, we consider a large-scale sparse WSN as shown in Figure7. That is a good benchmark since most of the WSN is spread over a large scale in the real world. Figure 7 shows the large scale sparse WSN in our simulation. There exists 100 sensor nodes and positions of sensors are randomly placed. The communication and sensing ranges of each sensor node are 20 m and 10 m, respectively. In Figure 7, the red line denotes true trajectory of the target. The state dynamic function and measurement function are similar to (19) and (21) whereas the parameter of process noise σw=1 m/s2 and the parameters of measurement noise σθ=pi/180,σr=2. Then we compared the CPF, DPF-WAV, and DPF-AV with the DPFF as well. In this simulation, the DPFF maintains 50 particles per node for each time step while CPF, DPF-WAV, and DPF-AV use 1000 particles.Figure 7
Large-scale sparse WSN.Figure8 illustrates the estimated tracks of all algorithms. It can be seen that the estimated trajectory of DPF-AV (cyan curve) is far away from the target trajectory (red curve). Compared with the DPF-WAV (blue curve), the estimated result (black curve) of the DPFF is more close to the target trajectory. The RMSE values of four algorithms versus time are shown in Figure 9. A remarkable fact is that the DPF-AV clearly diverge after 60 time steps. We note from Figure 9 that the RMSE of DPF-WAV is fluctuating with time depending on the number of nodes which detected the target. In addition, to provide an overall indication of the comparative performance of the different methods, Table 2 gives the average and standard deviation RMSE and the average runtime over four methods. In Table 2, it can be seen that the average runtime per node of DPFF is 0.1295 seconds, which accounts for approximately one-tenth of the DPF-WAV (0.9928 s). The most time-consuming algorithm is the CPF, which runs a particle filter in one central node based on all the measurements. In the end, the proposed algorithm significantly outperforms the alternative algorithms in this experiment.Figure 8
One example of estimated trajectory using DPFF (black curve), CPF (green curve), DPF-AV (cyan curve), and DPF-WAV (blue curve). The red curve indicates the ground truth in the large-scale sparse WSN.Figure 9
RMSE versus time.Table 2
Average results over 500 runs with different methods for random network.
Number
Methods
Average RMSE
Std RMSE
Average runtime (sec)
1000
CPF
0.19
0.13
5.7027
DPF-WAV
1.29
1.73
0.9928
DPF-AV
14.84
9.29
2.6869
50
DPFF
0.97
0.43
0.1295As a whole, the DPFF algorithm achieves a steady tracking with high accuracy with a few particles.


## 4.1. Example 1: Grid Network

Consider such a simulation that a target moves in a 300 m × 500 m area with 15 sensors. The sensors have overlapped monitored regions. The monitored region of each sensor node is assumed to be a circle region of 75 m radius whose center is at the sensor’s location. A sensor has a measurement of the target only if the target appears in the sensor’s sensing area. Figure3 illustrates the sensors and network connectivity. Connections between the sensors are shown as grey dashed lines. Each sensor can only communicate directly with its neighboring sensors whose distance to it is less than the communication range. The state vector of the target is represented as xk=xk,yk,vx,k,vy,k,wk, where xk,yk and vx,k,vy,k represent the target position and velocity, respectively. wk denotes the turn rate of the target. In this case, a nearly coordinated turn model with the known constant turn rate wk and the unknown velocity vk is considered. This model is able to account for the motion of complicated maneuverable targets. The nonlinear scenario is used in [32], and the motion of target is modeled according to
(19)xk+1=Fωkxk+Gwk,ωk+1=ωk+Δuk,where
(20)Fω=1sinωΔω0−1−cosωΔω0cosωΔ0−sinωΔ01−cosωΔω1sinωΔω0sinωΔ0cosωΔ,Gω=Δ220Δ00Δ220Δ.wk~N·;0,σw2I,uk~N·;0,σu2I,Δ=1 s, σw=15 m/s2, and σu=π/180 rad/s. The observation is a noisy bearing and range vector given by
(21)zk=arctanpx,kpy,kpx,k2+py,k2+εk,where εk~N·;0,Rk, with Rk=diagσθ2,σr2T and σθ=2π/180,σr=10. px,py are the distance from the target to the sensor.Figure 3
Object trajectory, sensor node, and connectivity considered in the experiments. The connectivity among the nodes is shown using the grey dashed lines between nodes.The number of particles for each node is to set 500 in the DPF-WAV and DPF-AV algorithms, while the DPFF algorithm only uses 50 particles in each node, while the CPF algorithm collects the entire available measurements from all sensor nodes and uses 500 particles to estimate.The error metric that we have computed is the root mean square error (RMSE) between the true and estimated target positions at each time instantk from all sensors. Let pkx,pky and p^kx,p^ky denote the true and estimated target positions, respectively, at time k. The RMSE value at time k is calculated over a number of Monte Carlo NMC runs according to
(22)RMSE=1NMC∑i=1NMCpkx−p^kx2+pky−p^ky2,where NMC=500 is the number of Monte Carlo runs.A sample of the estimated target track of DPFF is shown in Figure4, which also shows the estimated track from the CPF, DPF-WAV, and DPF-AV algorithms. It is obvious that the DPF-AV lost the track while other algorithms work well when tracking the target.Figure 4
One example of estimated trajectory using DPFF (black curve), CPF (green curve), DPF-AV (cyan curve), and DPF-WAV (blue curve). The red curve indicates the ground truth.Figure5 shows the temporal evolution of RMSE. It is observed that the DPF-AV algorithm performs worst in these methods, while results of CPF, DPFF, and DPF-WAV are fairly close to each other. The performance of DPFF and DPF-WAV is almost as good as the CPF and better than DPF-AV. Particularly, there is serious divergence with DPF-AV. For clear comparison of the proposed algorithm with DPF-WAV, we plot only DPF-WAV, CPF, and DPFF in Figure 6. Figure 6 illustrates the relationship of RMSE with timestamp of the DPF-WAV, CPF, and DPFF algorithms. As the DPF-WAV has a sharp fluctuation, it is obvious that the DPFF is better than the DPF-WAV and closer to the CPF. Especially, the DPF-WAV needs more particles than the DPFF and usually requires a significant amount of computing resources.Figure 5
Root mean square error (RMSE) versus time.Figure 6
Root mean square error (RMSE) versus time.Moreover, we also computed the averaged RMSE (ARMSE) for all methods. The ARMSE is given by(23)ARMSE=1TNMC∑i=1NMC∑j=1Tpkx−p^kx2+pky−p^ky2.Table1 lists the ARMSE, standard deviation of RMSE values, and the average execution time per node over 500 simulation trials. The first column in Table 1 is the number of particles for different algorithms. It can be seen that the proposed DPFF only with 50 particles gets better estimate results than the DPF-WAV and DPF-AV, almost equal to the CPF with one-tenth the number of particles. For the average execution time per node, the DPFF is less than a quarter of the DPF-WAV and DPF-AV. In summary, the DPFF has higher accuracy and less execution time.Table 1
Average results over 500 runs with different methods.
Number
Methods
Average RMSE
Std RMSE
Average runtime (sec)
500
CPF
6.33
4.43
0.3919
DPF-WAV
9.16
4.35
0.2851
DPF-AV
126.62
27.51
0.2705
1000
CPF
5.98
4.03
0.7815
DPF-WAV
8.97
4.12
0.3334
DPF-AV
124.12
21.48
0.3102
50
DPFF
6.32
4.76
0.0672


## 4.2. Example 2: Random Network

In the second scenario, we consider a large-scale sparse WSN as shown in Figure7. That is a good benchmark since most of the WSN is spread over a large scale in the real world. Figure 7 shows the large scale sparse WSN in our simulation. There exists 100 sensor nodes and positions of sensors are randomly placed. The communication and sensing ranges of each sensor node are 20 m and 10 m, respectively. In Figure 7, the red line denotes true trajectory of the target. The state dynamic function and measurement function are similar to (19) and (21) whereas the parameter of process noise σw=1 m/s2 and the parameters of measurement noise σθ=pi/180,σr=2. Then we compared the CPF, DPF-WAV, and DPF-AV with the DPFF as well. In this simulation, the DPFF maintains 50 particles per node for each time step while CPF, DPF-WAV, and DPF-AV use 1000 particles.Figure 7
Large-scale sparse WSN.Figure8 illustrates the estimated tracks of all algorithms. It can be seen that the estimated trajectory of DPF-AV (cyan curve) is far away from the target trajectory (red curve). Compared with the DPF-WAV (blue curve), the estimated result (black curve) of the DPFF is more close to the target trajectory. The RMSE values of four algorithms versus time are shown in Figure 9. A remarkable fact is that the DPF-AV clearly diverge after 60 time steps. We note from Figure 9 that the RMSE of DPF-WAV is fluctuating with time depending on the number of nodes which detected the target. In addition, to provide an overall indication of the comparative performance of the different methods, Table 2 gives the average and standard deviation RMSE and the average runtime over four methods. In Table 2, it can be seen that the average runtime per node of DPFF is 0.1295 seconds, which accounts for approximately one-tenth of the DPF-WAV (0.9928 s). The most time-consuming algorithm is the CPF, which runs a particle filter in one central node based on all the measurements. In the end, the proposed algorithm significantly outperforms the alternative algorithms in this experiment.Figure 8
One example of estimated trajectory using DPFF (black curve), CPF (green curve), DPF-AV (cyan curve), and DPF-WAV (blue curve). The red curve indicates the ground truth in the large-scale sparse WSN.Figure 9
RMSE versus time.Table 2
Average results over 500 runs with different methods for random network.
Number
Methods
Average RMSE
Std RMSE
Average runtime (sec)
1000
CPF
0.19
0.13
5.7027
DPF-WAV
1.29
1.73
0.9928
DPF-AV
14.84
9.29
2.6869
50
DPFF
0.97
0.43
0.1295As a whole, the DPFF algorithm achieves a steady tracking with high accuracy with a few particles.


## 5. Conclusions

We presented a distributed particle flow filter algorithm for wireless sensor networks. At each sensor, a local particle flow filter computes a local state estimate that only depends on the local measurement. Then a perception confidence value is calculated from the particle flow filter. A max-consensus is used to make all the nodes agree on the best estimate of the target position. In the proposed distributed particle flow filter, each node just communicates with its neighboring sensor nodes and does not require any routing protocols. We applied the proposed distributed particle flow filter in two target tracking scenarios and demonstrated experimentally that its performance is better than the distributed averaged-based particle filter; moreover, it needs less computation time and samples. An extension of the distributed particle flow filter to multiple target tracking in the WSN remains a potential topic for the future.



---

*Source: 5387142-2018-04-11.xml*