# An Advanced Auxiliary Delay-Weight Particle Filter with Linear Computation Cost

Authors: Chen Li; Lin Sun; Zengwei Zheng; Dan Chen
Journal: International Journal of Distributed Sensor Networks
 (2016)

---

## Abstract

We investigate the problem of tracking mobile targets in wireless sensor networks. We propose an advanced auxiliary delayed-weight particle filter algorithm (ADWPF). We make a deep study on the evolvement of particles and formally define the tree-like structure relationship among particles based on observations. Most importantly, we add some auxiliary particles to these structures formed by sampled particles in order to obtain more efficient structures. Based on the newly tree-like structures formed by auxiliary particles and sampled particles, we design a well efficient delayed-weight algorithm with linear computation cost. Experiment results demonstrate that our algorithm can greatly improve the tracking accuracy of a mobile target, compared with bootstrap filter, auxiliary particle filter, and another delayed-weight particle filter.


---

## Body

## 1. Introduction

In recent years, tracking mobile targets in wireless sensor networks has attracted more and more attention of the researchers in academia and industry [1–3]. This technology has fairly good application prospects in indoor environments, such as shopping mall and underground parking.Obviously, particle filter is one of the most important technologies for tracking mobile targets in wireless senor networks. Many researchers in [4–6] have utilized particle filter to track the mobile targets in wireless sensor networks. Obviously, the performance of particle filter algorithm has a great impact on the tracking accuracy. So it is very meaningful and important to design an efficient particle filter algorithm. This is also the motivation of our research.In our paper, we propose an auxiliary delayed-weight particle filter algorithm (ADWPF). Our main contributions are as follows:(i)
We make a deep study on the evolvement of particles and give the formal definition of the tree-like structure relationship among particles based on observations.(ii)
We add some auxiliary particles to the well-defined tree-like structures formed by sampled particles in order to get more efficient structures. These auxiliary particles can greatly help to improve tracking accuracy.(iii)
Based on the newly tree-like structures formed by auxiliary particles and sampled particles, we design a well efficient delayed-weight algorithm with linear computation cost, whose time complexity isOΔTN; here Δ is the delayed time length, T is the number of tracking steps, and N is the number of samples in each step.The remainder of this paper is organized as follows. Section2 reviews the related work in the area of particle filter and delayed-state filtering. In Section 3, we briefly overview the particle filter and describe the key technologies applied in our algorithm. In Section 4, we describe the deduction and the framework of our algorithm. Section 5 presents simulation results for a self-navigation problem in wireless sensor networks. Section 6 makes the conclusion.


## 2. Related Work

The particle filter algorithm has been widely used to infer the latent states of stochastic dynamic systems in [7]. It uses the sequential nature of stochastic dynamic systems to generate sequentially weighted Monte Carlo samples, and these samples are used to infer the latent states. For example, in the problem of self-navigation for a mobile target in wireless sensor networks, the position of the mobile target at each time step is the latent state to be inferred.The bootstrap filter in [8] was the first successful application of the particle filter to the field of nonlinear filtering. It utilizes transition density function as the importance distribution to generate Monte Carlo samples of latent states. Then it assigns the weight of each sample based on the likelihood of the measurement. However, outliers that occurred in the measurement data will decrease the performance of the bootstrap filter as pointed in [8].In order to increase the robustness of the sampling procedure and obtain approximate samples from the optimal importance distribution, which could minimize the variance of importance weights, [9] introduced the auxiliary particle filter. This algorithm can efficiently reduce the variance of samples’ weights.Reference [10] pointed out that the “future” measurement information can be utilized to improve the estimation accuracy of current system state and mitigate degeneracy of the importance weights over time. So the approximations of the “earlier” marginal distributions of latent states will not collapse as time increases. Reference [10] introduced some delayed-sample methods, such as the Exact Lookahead Sampling method in [11] and the Block Sampling method in [12]. The key technology of these methods is to design an efficient importance distribution, which utilizes “future” information to generate weighted samples of the current states.In fact, the delayed-sample algorithms which utilize “future” information to generate samples of current states are very similar to the smoothing algorithms which utilize current received measurement information to smooth the marginal filtering distributions of historical hidden states. Reference [10] also pointed that these two methods have very close relationship. Traditional smoothing algorithms mainly include the forward filtering-backward smoothing algorithm and the two-filter smoothing algorithm in [13]. The former is relying on the particle implementation of the forward filtering-backward smoothing formula to reweight the historical filtered samples, and the latter is based on the particle implementation of the two-filter smoothing formula to regenerate weighted samples for the historical hidden states. These two methods have time complexity OΔN2 in [13], where Δ is the time interval and N is the number of samples in each step.


## 3. Particle Filter Algorithm

### 3.1. System Description

Suppose that the dynamic system [14] is given by the following state and measurement equations:(1)xt+1=fxt+wt,(2)yt+1=gxt+1+vt,where xt is state, yt+1 is measurement, and wt,vt~N0,1. Therefore, the transition probability function of state variable xt+1 is xt+1~pf(·∣xt) and the likelihood function of measurement yt+1 is yt+1~pg(·∣xt+1). The joint distribution of measurement information y1:t=y1,…,yt and state sequence x1:t=x1,…,xt is (3)px1:t,y1:t=πx1∏k=1tpfxk+1∣xk∏k=1tpgyk∣xk,where πx1 is the prior probability distribution of initial state. So the posterior distribution of x1:t is(4)px1:t∣y1:t=px1:t,y1:tpy1:t.


### 3.2. Particle Filter

The aim of particle filter is to get the sequential Monte Carlo samplesX1:ti~px1:t∣y1:t, here i=1,2,…,N, so the Monte Carlo approximation of px1:t∣y1:t is(5)p^x1:t∣y1:t=1N∑i=1NδX1:tix1:t,where δ· is Delta function. The particle filter generally consists of the steps shown in Algorithm 1 (take bootstrap filter algorithm as an example).Algorithm 1:Bootstrap filter algorithm.
Step  1 (Initialization). At time t=1, sample X~1i~πx1 to obtain the Monte Carlo approximation of px1∣y1
p~x1∣y1=∑i=1NW1iδX~1ix1,
whereW1i∝pgy1∣X~1i is the weight of X~1i, resample X1i~p~x1∣y1 to obtain
p^x1∣y1=1N∑i=1NδX1ix1.
Step  2 (Filtering). At time t≥2, sample X~ti~pf⋅∣Xt-1i to obtain the Monte Carlo approximation of pxt∣y1:t
p~xt∣y1:t=∑i=1NWtiδX~tixt,
whereWti∝pgyt∣X~ti is the weight of X~ti, resample Xti~p~xt∣y1:t to obtain
p^xt∣y1:t=1N∑i=1NδXtixtThe time complexity isOTN at each time step, and memory complexity is OTN.


### 3.3. Delay Strategy

The key idea of this strategy is that the dynamic system offers more and more information of the latent states as time increases. In order to formalize this fact, assume that the measurement information available at timet+Δ takes the form y1:t+Δ=y1,…,yt+Δ, the true system states are x1:t=x1,…,xt, and Ex1:t∣y1:t+Δ is the expectation of x1:t based on y1:t+Δ. Let X^t+Δ be a consistent Monte Carlo of Ex1:t∣y1:t+Δ. The mean squared difference between X^t+Δ and x1:t is EX^t+Δ-x1:t2 that follows.Proposition 1.
When the Monte Carlo sample sizeN is sufficiently large, EX^t+Δ-x1:t2 decreases as Δ increases.For the proof one can refer to [10].Based on the statement of this proposition, we can design an algorithm that utilizes more and more information to improve the inference about the latent state in a specific time. This lays down the fundament of our work.There are mainly two methods to implement the idea of Proposition1 [10]. The first method is the delayed-sample method. It designed an importance distribution, which utilizes measurement information after time t to regenerate weighted samples of state at time t, for example, auxiliary particle filter. The second method is the delayed-weight method, which utilized information after time t to reweight samples Xtii=1N.


### 3.4. Relation between Delayed-Sample Method and Delayed-Weight Method

#### 3.4.1. The Problem of Delayed-Sample Method

We know that the most important technology in particle filter is to choose an efficient importance distribution. The optimal importance distribution of the delayed-sample method has the following form [10]:(6)X^ti~qtxt∣Xt-1i,yt:t+Δ∝∫pXt-1i,xt,xt+1,…,xt+Δ,yt:t+Δdxt+1⋯dxt+Δ=pXt-1i,xt,yt:t+Δ.This optimal importance distribution can generate more effective samples than standard particle filter method, which is based on the one-at-a-time sampling strategy. For continuous state space, it is very hard to adopt this method, since it is not easy to get the explicit form of (6).


#### 3.4.2. Delayed-Weight Method Is a Better Choice

In order to solve the problem of the delayed-sample method, we can discretize the continuous state space using sampled particles, which is the core idea of the delayed-weight method. Figure1 illustrates the evolvement of particles in the standard particle filter method at time t.Figure 1
The evolvement of particles in standard particle filter method at timet.Assume that we observe measurementsy1,…,yt+Δ and have got weighted particles X1i,w1ii=1N,…,Xt+Δi,wt+Δii=1N using a forward particle filter method until time t+Δ; then we can use these particles to represent the continuous state space, as illustrated in Figure 2.Figure 2
Continuous state space can be discretized by weighted particles.The main advantage of doing so is that the optimal importance distributionpXt-1i,xt,yt:t+Δ in (6) can be approximated by(7)pXt-1i,xt,yt:t+Δ≈∑it+1=1N⋯∑it+Δ=1NpXt-1i,xt,Xt+1it+1,…,Xt+Δit+Δ,yt:t+Δ.Therefore, the delayed weightw^ti of particle Xti satisfies (8)w^ti∝pXt-1i,Xti,yt:t+Δ=∑it+1=1N⋯∑it+Δ=1NpXt-1i,Xti,Xt+1it+1,…,Xt+Δit+Δ,yt:t+Δ.Obviously, this delayed weight can be calculated by the forward filtering-backward smoothing method. We name this method DWFB in [11, 15].The main disadvantage of this method in [15] is that, in each iteration, we must make use of all particles at time k+1 to calculate the delayed weight of each particle at time k; for example, the delayed weight w^ki of Xki is(9)w^ki∝∑j=1Nw^k+1jpfXk+1j∣Xki,i=1,…,N.Since the number of particles isN at time k, the time complexity of this algorithm is OΔTN2; when N is large, the computation cost is very high. In order to reduce the time complexity, we propose an advanced delayed-weight particle filter weight method, which not only reduces the time complexity to OΔTN but also greatly improves the tracking accuracy.


## 3.1. System Description

Suppose that the dynamic system [14] is given by the following state and measurement equations:(1)xt+1=fxt+wt,(2)yt+1=gxt+1+vt,where xt is state, yt+1 is measurement, and wt,vt~N0,1. Therefore, the transition probability function of state variable xt+1 is xt+1~pf(·∣xt) and the likelihood function of measurement yt+1 is yt+1~pg(·∣xt+1). The joint distribution of measurement information y1:t=y1,…,yt and state sequence x1:t=x1,…,xt is (3)px1:t,y1:t=πx1∏k=1tpfxk+1∣xk∏k=1tpgyk∣xk,where πx1 is the prior probability distribution of initial state. So the posterior distribution of x1:t is(4)px1:t∣y1:t=px1:t,y1:tpy1:t.


## 3.2. Particle Filter

The aim of particle filter is to get the sequential Monte Carlo samplesX1:ti~px1:t∣y1:t, here i=1,2,…,N, so the Monte Carlo approximation of px1:t∣y1:t is(5)p^x1:t∣y1:t=1N∑i=1NδX1:tix1:t,where δ· is Delta function. The particle filter generally consists of the steps shown in Algorithm 1 (take bootstrap filter algorithm as an example).Algorithm 1:Bootstrap filter algorithm.
Step  1 (Initialization). At time t=1, sample X~1i~πx1 to obtain the Monte Carlo approximation of px1∣y1
p~x1∣y1=∑i=1NW1iδX~1ix1,
whereW1i∝pgy1∣X~1i is the weight of X~1i, resample X1i~p~x1∣y1 to obtain
p^x1∣y1=1N∑i=1NδX1ix1.
Step  2 (Filtering). At time t≥2, sample X~ti~pf⋅∣Xt-1i to obtain the Monte Carlo approximation of pxt∣y1:t
p~xt∣y1:t=∑i=1NWtiδX~tixt,
whereWti∝pgyt∣X~ti is the weight of X~ti, resample Xti~p~xt∣y1:t to obtain
p^xt∣y1:t=1N∑i=1NδXtixtThe time complexity isOTN at each time step, and memory complexity is OTN.


## 3.3. Delay Strategy

The key idea of this strategy is that the dynamic system offers more and more information of the latent states as time increases. In order to formalize this fact, assume that the measurement information available at timet+Δ takes the form y1:t+Δ=y1,…,yt+Δ, the true system states are x1:t=x1,…,xt, and Ex1:t∣y1:t+Δ is the expectation of x1:t based on y1:t+Δ. Let X^t+Δ be a consistent Monte Carlo of Ex1:t∣y1:t+Δ. The mean squared difference between X^t+Δ and x1:t is EX^t+Δ-x1:t2 that follows.Proposition 1.
When the Monte Carlo sample sizeN is sufficiently large, EX^t+Δ-x1:t2 decreases as Δ increases.For the proof one can refer to [10].Based on the statement of this proposition, we can design an algorithm that utilizes more and more information to improve the inference about the latent state in a specific time. This lays down the fundament of our work.There are mainly two methods to implement the idea of Proposition1 [10]. The first method is the delayed-sample method. It designed an importance distribution, which utilizes measurement information after time t to regenerate weighted samples of state at time t, for example, auxiliary particle filter. The second method is the delayed-weight method, which utilized information after time t to reweight samples Xtii=1N.


## 3.4. Relation between Delayed-Sample Method and Delayed-Weight Method

### 3.4.1. The Problem of Delayed-Sample Method

We know that the most important technology in particle filter is to choose an efficient importance distribution. The optimal importance distribution of the delayed-sample method has the following form [10]:(6)X^ti~qtxt∣Xt-1i,yt:t+Δ∝∫pXt-1i,xt,xt+1,…,xt+Δ,yt:t+Δdxt+1⋯dxt+Δ=pXt-1i,xt,yt:t+Δ.This optimal importance distribution can generate more effective samples than standard particle filter method, which is based on the one-at-a-time sampling strategy. For continuous state space, it is very hard to adopt this method, since it is not easy to get the explicit form of (6).


### 3.4.2. Delayed-Weight Method Is a Better Choice

In order to solve the problem of the delayed-sample method, we can discretize the continuous state space using sampled particles, which is the core idea of the delayed-weight method. Figure1 illustrates the evolvement of particles in the standard particle filter method at time t.Figure 1
The evolvement of particles in standard particle filter method at timet.Assume that we observe measurementsy1,…,yt+Δ and have got weighted particles X1i,w1ii=1N,…,Xt+Δi,wt+Δii=1N using a forward particle filter method until time t+Δ; then we can use these particles to represent the continuous state space, as illustrated in Figure 2.Figure 2
Continuous state space can be discretized by weighted particles.The main advantage of doing so is that the optimal importance distributionpXt-1i,xt,yt:t+Δ in (6) can be approximated by(7)pXt-1i,xt,yt:t+Δ≈∑it+1=1N⋯∑it+Δ=1NpXt-1i,xt,Xt+1it+1,…,Xt+Δit+Δ,yt:t+Δ.Therefore, the delayed weightw^ti of particle Xti satisfies (8)w^ti∝pXt-1i,Xti,yt:t+Δ=∑it+1=1N⋯∑it+Δ=1NpXt-1i,Xti,Xt+1it+1,…,Xt+Δit+Δ,yt:t+Δ.Obviously, this delayed weight can be calculated by the forward filtering-backward smoothing method. We name this method DWFB in [11, 15].The main disadvantage of this method in [15] is that, in each iteration, we must make use of all particles at time k+1 to calculate the delayed weight of each particle at time k; for example, the delayed weight w^ki of Xki is(9)w^ki∝∑j=1Nw^k+1jpfXk+1j∣Xki,i=1,…,N.Since the number of particles isN at time k, the time complexity of this algorithm is OΔTN2; when N is large, the computation cost is very high. In order to reduce the time complexity, we propose an advanced delayed-weight particle filter weight method, which not only reduces the time complexity to OΔTN but also greatly improves the tracking accuracy.


## 3.4.1. The Problem of Delayed-Sample Method

We know that the most important technology in particle filter is to choose an efficient importance distribution. The optimal importance distribution of the delayed-sample method has the following form [10]:(6)X^ti~qtxt∣Xt-1i,yt:t+Δ∝∫pXt-1i,xt,xt+1,…,xt+Δ,yt:t+Δdxt+1⋯dxt+Δ=pXt-1i,xt,yt:t+Δ.This optimal importance distribution can generate more effective samples than standard particle filter method, which is based on the one-at-a-time sampling strategy. For continuous state space, it is very hard to adopt this method, since it is not easy to get the explicit form of (6).


## 3.4.2. Delayed-Weight Method Is a Better Choice

In order to solve the problem of the delayed-sample method, we can discretize the continuous state space using sampled particles, which is the core idea of the delayed-weight method. Figure1 illustrates the evolvement of particles in the standard particle filter method at time t.Figure 1
The evolvement of particles in standard particle filter method at timet.Assume that we observe measurementsy1,…,yt+Δ and have got weighted particles X1i,w1ii=1N,…,Xt+Δi,wt+Δii=1N using a forward particle filter method until time t+Δ; then we can use these particles to represent the continuous state space, as illustrated in Figure 2.Figure 2
Continuous state space can be discretized by weighted particles.The main advantage of doing so is that the optimal importance distributionpXt-1i,xt,yt:t+Δ in (6) can be approximated by(7)pXt-1i,xt,yt:t+Δ≈∑it+1=1N⋯∑it+Δ=1NpXt-1i,xt,Xt+1it+1,…,Xt+Δit+Δ,yt:t+Δ.Therefore, the delayed weightw^ti of particle Xti satisfies (8)w^ti∝pXt-1i,Xti,yt:t+Δ=∑it+1=1N⋯∑it+Δ=1NpXt-1i,Xti,Xt+1it+1,…,Xt+Δit+Δ,yt:t+Δ.Obviously, this delayed weight can be calculated by the forward filtering-backward smoothing method. We name this method DWFB in [11, 15].The main disadvantage of this method in [15] is that, in each iteration, we must make use of all particles at time k+1 to calculate the delayed weight of each particle at time k; for example, the delayed weight w^ki of Xki is(9)w^ki∝∑j=1Nw^k+1jpfXk+1j∣Xki,i=1,…,N.Since the number of particles isN at time k, the time complexity of this algorithm is OΔTN2; when N is large, the computation cost is very high. In order to reduce the time complexity, we propose an advanced delayed-weight particle filter weight method, which not only reduces the time complexity to OΔTN but also greatly improves the tracking accuracy.


## 4. An Auxiliary Delayed-Weight Method

Our method consists of two key technologies: designing a backward smoothing procedure and adding auxiliary particles.


### 4.1. A Backward Smoothing Procedure

Definition 2.
If one particleXt at time t is sampled from the importance distribution qxt∣Xt-1, where Xt-1 is a particle at time t-1 and t=1,…,T, one says that Xt is the child of Xt-1, denoted by Xt=childXt-1 or, conversely, that Xt-1 is the parent of Xt, denoted by Xt-1=parent(Xt). Furthermore, if Xt has a child Xt+1, then Xt+1 is the grandchild of Xt-1, denoted by Xt+1=grandchildXt-1, and Xt-1 is the grandparent of Xt+1, denoted by Xt-1=grandparentXt+1.
Observing the evolvement of particles in Figure1, some particles have more than one child, while some have no child. Particularly, we observe that a particle and its child/grandchild form a tree-like structure, which is straightforward in Figure 2.Definition 3.
The tree that consists of particleX and its children/grandchildren is called a generation tree, denoted by tree(X). X is the root of tree(X).
As illustrated in Figure3, the particles denoted by green circles construct the generation tree tree(Xti) and Xti is the root of tree(Xti).Figure 3
The particles denoted by green circles formtree(Xti); Xti is the root of tree(Xti).Assumption 4.
If the particleXt+Δ is not the child/grandchild of Xt, then pfXt+k∣Xt=0, where pf· is the transition probability function, t=1,…,T, and k>0.The main advantage of doing so is that we only need to make use of children/grandchildren ofXti to calculate the delayed weight w^ki of Xki, shown in Figure 4. Therefore, (8) can be reduced to(10)w^ti∝pXt-1i,Xti,yt:t+Δ≈∑it+1=1Nt+1⋯∑it+Δ=1Nt+ΔpXt-1i,Xti,Xt+1it+1,…,Xt+Δit+Δ,yt:t+Δ,where Nt+j is the number of children/grandchildren of Xti at time t+j, j=1,…,Δ.Figure 4
We only need to use children/grandchildren ofXti, which are denoted by green circles, to calculate the delayed weight w^ki of Xki.Observing the structure oftree(Xti) in Figure 4(b), we design an advanced delayed-weight algorithm whose time complexity is OΔTN to calculate w^ki. Take tree(Xti) in Figure 4(b) as an example; here Δ=3. Since there are 5 leaf nodes in tree(Xti), the delayed weight w^ki of Xki can be calculated as follows:(11)w^ti∝pXti,Xt+1i+2,yt:t+1+pXti,Xt+1i+1,Xt+2i+1,yt:t+2+pXti,Xt+1i+1,Xt+2i+2,Xt+3i+1,yt:t+3+pXti,Xt+1i+1,Xt+2i+2,Xt+3i+2,yt:t+3+pXti,Xt+1i+3,Xt+2i+3,Xt+3i+3,yt:t+3,where Xt+2i+1,Xt+3i+1,Xt+3i+2,Xt+1i+2,Xt+3i+3 are leaf nodes of tree(Xti). Importantly, we find that every term in the right side of (11) includes one leaf node.Definition 5.
One calls the joint probability function representing the parent-child relation between leaf nodeXleaf and root node X the backward smoothing term, denoted by PXleaf→X; for example, PXt+1i+2→Xti=pXti,Xt+1i+2,yt:t+1.
Thereforew^ti in (11) can be rewritten as (12)w^ti∝PXt+1i+2⟶Xti+PXt+2i+1⟶Xti+PXt+3i+1⟶Xti+PXt+3i+2⟶Xti+PXt+3i+3⟶Xti.
More importantly, each backward smoothing term can be decomposed as the product of transition probability functions and likelihood functions; for example, the backward smoothing termPXt+2i+1→Xti can be decomposed as follows:(13)PXt+2i+1⟶Xti=pXti,Xt+1i+1,Xt+2i+1,yt:t+2=pXtipfXt+1i+1∣XtipfXt+2i+1∣Xt+1i+1∗pgyt∣Xtipgyt+1∣Xt+1i+1pgyt+2∣Xt+2i+1.
Other backward smoothing terms can be decomposed in the same way. Moreover, from (13), we can see that it is very straightforward to use an iterate method to calculate each backward smoothing term, which is the core idea of our backward smoothing algorithm. Procedure 1 gives the backward smoothing (BS) process from one leaf node to the root node.Procedure 1:BS.
Input: a leaf node Xti
Output: the backward smoothing term PXti→Xt-Δroot, where Xt-Δroot is the grandparent of Xti at time t-Δ.
(1) Initialization: PXti→Xt-Δroot=1;
(2) Xtemp=Xti;
(3) for k=t,…,t-Δ+1
// the data structure of Xkj includes a pointer to his parent
(4)   Xparent=parent(Xtemp); //a backward smoothing
//operation
(5)   PXti→Xt-Δroot=PXti→Xt-Δroot∗pgyk∣Xtemp∗pfXtemp∣Xparent;
(6)   Xtemp=Xparent;
end
(7) return PXti→Xt-Δroot;Obviously, the time complexity of Procedure1 is OΔ. Figure 5 shows the backward smoothing procedure.Figure 5
The backward smoothing procedure of Procedure1. The calculation of each backward smoothing term is iterated from one leaf node to the root node.


### 4.2. Adding Auxiliary Particles

Experiment results demonstrate that it is not enough to use particles intree(Xti) to calculate the delayed weight w^ki of Xki. We should add some auxiliary particles to tree(Xti) as illustrated in Figure 6.Figure 6
Add auxiliary particles totree(Xti) to get a new tree aux_tree(Xti). Auxiliary particles are denoted by red circles. All leaf nodes of aux_tree(Xti) have the same depth, which equals Δ=3.Definition 6.
One calls the modified generation tree to which auxiliary particles have been added the auxiliary generation tree, denoted byaux_tree(X), where X is the root of aux_tree(X).
In Figure6, we also call these newly added auxiliary particles, like ut+21, ut+31, and ut+32, grandchildren of Xti. The most important feature of aux_tree(Xti) is that all of its leaf nodes have the same depth, which equals the delayed time length Δ.Definition 7.
The state of auxiliary particle is defined as the expectation of transition probability function; that is, the auxiliary particlesut+21 and ut+31 in Figure 6 are defined as follows:(14)ut+21=Epf·∣Xt+1i+2,ut+31=Epf·∣ut+21.
Furthermore, according to the definition of transition probability function (1), we have (15)ut+21=fXt+1i+2+Ewt=fXt+1i+2,ut+31=fut+21+Ewt=fut+21,where since wt~N0,1, Ewt=0.
According toaux_tree(Xti) in Figure 6, the delayed weight w^ki of Xki can be calculated as follows:(16)w^ti∝pXti,Xt+1i+1,Xt+2i+1,ut+31,yt:t+3+pXti,Xt+1i+1,Xt+2i+2,Xt+3i+1,yt:t+3+pXti,Xt+1i+1,Xt+2i+2,Xt+3i+2,yt:t+3+pXti,Xt+1i+2,ut+31,ut+32,yt:t+1+pXti,Xt+1i+3,Xt+2i+3,Xt+3i+3,yt:t+3,where ut+31,Xt+3i+1,Xt+3i+2,ut+32,Xt+3i+3 are leaf nodes of aux_tree(Xti). It is a little different from (11). And, according to the definition of backward smoothing term, w^ti can also be rewritten as(17)w^ti∝Put+31⟶Xti+PXt+3i+1⟶Xti+PXt+3i+2⟶Xti+Put+32⟶Xti+PXt+3i+3⟶Xti.
Absolutely, each backward smoothing term in (17) can be calculated by Procedure 1.The main advantage of adding auxiliary particles to the generation tree is that they can greatly improve the filtering accuracy, which is verified in our experiments.Procedure2 gives the process of constructing an auxiliary generation tree (CAGF).Procedure 2:CAGF.
LetΔ is the height of treeXti, Nk is the number of child/grandchildren of Xti at time k, where t+1≤k≤t+Δ.
Input: treeXti
Output: an auxiliary generation tree aux_treeXti
(1)Initialization: index=0; //index auxiliary particles
(2)if Xti has no child
(3)X=Xti;
(4)for j=k+1,…,t+Δ
(5) ujindex=Epf⋅∣X;
(6) X=ujindex;
(7) index=index+1;
end
end
(8)else for k=t+1,…,t+Δ-1
(9)for n=1,…,Nk
(10)  if Xkn is a leaf node
(11)   X=Xkn;
//add auxiliary particles
(12)   for j=k+1,…,t+Δ
(13)    ujindex=Epf⋅∣X;
(14)    X=ujindex;
(15)    index=index+1;
end
end
end
end
(16) return aux_treeXtiLemma 8.
The number of leaf nodes intreeXti is equal to the number of leaf nodes in aux_treeXti.Proof.
According to step (13) in Procedure 2, only one child is added to every leaf node. So the number of leaf nodes will not increase.Lemma 9.
LetNk denote the number of children/grandchildren of Xti in treeXti at time k, where t+1≤k≤t+Δ. The time complexity TCAGFXti of Procedure 2 satisfies (18)O∑k=t+1t+ΔNk≤TCAGFXti≤OΔ+1∑k=t+1t+ΔNk.The proof is obvious.


### 4.3. Algorithm Framework

In this section, we first give the framework of our algorithm, and then we analyze its time complexity.Theorem 10.
The time complexityTADWPF of Algorithm 3 is OCN, C≪N.Proof.
Obviously, the time cost of Algorithm3 mainly consists of steps (4) and (6). Let TCAGFXt-Δi, 1≤i≤N, denote the time cost of step (4) in ADWPF, and Nki denotes the number of children/grandchildren of Xt-Δi in treeXt-Δi at time k, where t-Δ+1≤k≤t. According to Lemma 9(19)∑i=1NTCAGFXt-Δi≤∑i=1NOΔ+1∑k=t-Δ+1tNki=OΔ+1∑i=1N∑k=t-Δ+1tNki=OΔ+1N,where the second equation is based on the fact that we only sampled N particles at each time step according to step (1) in ADWPF. And according to Lemma 8(20)∑i=1NTCAGFXt-Δi≥∑i=1NO∑k=t-Δ+1tNki=O∑i=1N∑k=t-Δ+1tNki=ON.
LetTBSXt-Δi, 1≤i≤N, denote the time cost of step (6) in ADWPF. The sum number of leaf nodes of all generation trees treeXt-Δ1,…,treeXt-ΔN is less than ΔN and greater than N, so the sum number of leaf nodes of all auxiliary generation trees aux_treeXt-Δ1,…,aux_treeXt-ΔN is less than ΔN and greater than N according to Lemma 8. Since the time cost of Procedure 1 is OΔ for each leaf node, we have (21)N∗OΔ≤∑i=1NTBSXt-Δi≤ΔN∗OΔ=OΔ2N.
So for the time complexityTADWPF of Algorithm 2 we have(22)TADWPF=T∗∑i=1NTCAGFXt-Δi+∑i=1NTBSXt-Δi≤T∗OΔ+1N+OΔ2N=OΔ2+Δ+1TN.
Meanwhile,(23)TADWPF=T∗∑i=1NTCAGFXt-Δi+∑i=1NTBSXt-Δi≥T∗ON+OΔN=OΔ+1TN,where Δ,T≪N. So the time complexity of ADWPF is OCN, C≪N.Algorithm 2:DWFB.
AssumeN is the number of particles, t is current time index, count is the number of simulation steps, Δ is the delayed time interval,
w
t
i is the filtered weight and w^ti is the modified weight.
(1)  for t=1,…,T
(2)  Generate samples Xti,wti,w^ti=0i=1N through a forward particle filter (like Algorithm 1).
(3)  if t>Δ
(4)  for i=1:N
(5)   w^ti=wti
(6)   for k=T-1:T-∆
(7)   for i=1:N
w^ki=∑j=1Nw^k+1jwkipfXk+1j∣Xki∑lNwklpfXk+1j∣Xkl
end
end
end
end
The inference result of system statexT-Δ is
x^t-Δ=∑i=1Nw^t-Δi⋅Xt-ΔiAlgorithm 3:ADWPF.
LetN is the number of particles, t is current time index, Δ is the delayed time length,wti is the forward filtered weight
andw^ti is the delayed weight.
(
1) for t=1,…,T
Generate samples Xti,wti,w^ti=0i=1N through a forward particle filter (like Algorithm 1).
end
(2) if t>Δ
(3)  for i=1,…,N
//Procedure 2.
(4)   aux_treeXt-Δi=CAGFtreeXt-Δi,Δ
(5)   for each leaf node Xtj  of aux_treeXt-Δi
//Procedure 1.
(6)    PXtj→Xt-Δi=BSXtj;
(7)    w^t-Δi=w^t-Δi+PXtj→Xt-Δi;
(8)    X=Xkn;
end
end
end


## 4.1. A Backward Smoothing Procedure

Definition 2.
If one particleXt at time t is sampled from the importance distribution qxt∣Xt-1, where Xt-1 is a particle at time t-1 and t=1,…,T, one says that Xt is the child of Xt-1, denoted by Xt=childXt-1 or, conversely, that Xt-1 is the parent of Xt, denoted by Xt-1=parent(Xt). Furthermore, if Xt has a child Xt+1, then Xt+1 is the grandchild of Xt-1, denoted by Xt+1=grandchildXt-1, and Xt-1 is the grandparent of Xt+1, denoted by Xt-1=grandparentXt+1.
Observing the evolvement of particles in Figure1, some particles have more than one child, while some have no child. Particularly, we observe that a particle and its child/grandchild form a tree-like structure, which is straightforward in Figure 2.Definition 3.
The tree that consists of particleX and its children/grandchildren is called a generation tree, denoted by tree(X). X is the root of tree(X).
As illustrated in Figure3, the particles denoted by green circles construct the generation tree tree(Xti) and Xti is the root of tree(Xti).Figure 3
The particles denoted by green circles formtree(Xti); Xti is the root of tree(Xti).Assumption 4.
If the particleXt+Δ is not the child/grandchild of Xt, then pfXt+k∣Xt=0, where pf· is the transition probability function, t=1,…,T, and k>0.The main advantage of doing so is that we only need to make use of children/grandchildren ofXti to calculate the delayed weight w^ki of Xki, shown in Figure 4. Therefore, (8) can be reduced to(10)w^ti∝pXt-1i,Xti,yt:t+Δ≈∑it+1=1Nt+1⋯∑it+Δ=1Nt+ΔpXt-1i,Xti,Xt+1it+1,…,Xt+Δit+Δ,yt:t+Δ,where Nt+j is the number of children/grandchildren of Xti at time t+j, j=1,…,Δ.Figure 4
We only need to use children/grandchildren ofXti, which are denoted by green circles, to calculate the delayed weight w^ki of Xki.Observing the structure oftree(Xti) in Figure 4(b), we design an advanced delayed-weight algorithm whose time complexity is OΔTN to calculate w^ki. Take tree(Xti) in Figure 4(b) as an example; here Δ=3. Since there are 5 leaf nodes in tree(Xti), the delayed weight w^ki of Xki can be calculated as follows:(11)w^ti∝pXti,Xt+1i+2,yt:t+1+pXti,Xt+1i+1,Xt+2i+1,yt:t+2+pXti,Xt+1i+1,Xt+2i+2,Xt+3i+1,yt:t+3+pXti,Xt+1i+1,Xt+2i+2,Xt+3i+2,yt:t+3+pXti,Xt+1i+3,Xt+2i+3,Xt+3i+3,yt:t+3,where Xt+2i+1,Xt+3i+1,Xt+3i+2,Xt+1i+2,Xt+3i+3 are leaf nodes of tree(Xti). Importantly, we find that every term in the right side of (11) includes one leaf node.Definition 5.
One calls the joint probability function representing the parent-child relation between leaf nodeXleaf and root node X the backward smoothing term, denoted by PXleaf→X; for example, PXt+1i+2→Xti=pXti,Xt+1i+2,yt:t+1.
Thereforew^ti in (11) can be rewritten as (12)w^ti∝PXt+1i+2⟶Xti+PXt+2i+1⟶Xti+PXt+3i+1⟶Xti+PXt+3i+2⟶Xti+PXt+3i+3⟶Xti.
More importantly, each backward smoothing term can be decomposed as the product of transition probability functions and likelihood functions; for example, the backward smoothing termPXt+2i+1→Xti can be decomposed as follows:(13)PXt+2i+1⟶Xti=pXti,Xt+1i+1,Xt+2i+1,yt:t+2=pXtipfXt+1i+1∣XtipfXt+2i+1∣Xt+1i+1∗pgyt∣Xtipgyt+1∣Xt+1i+1pgyt+2∣Xt+2i+1.
Other backward smoothing terms can be decomposed in the same way. Moreover, from (13), we can see that it is very straightforward to use an iterate method to calculate each backward smoothing term, which is the core idea of our backward smoothing algorithm. Procedure 1 gives the backward smoothing (BS) process from one leaf node to the root node.Procedure 1:BS.
Input: a leaf node Xti
Output: the backward smoothing term PXti→Xt-Δroot, where Xt-Δroot is the grandparent of Xti at time t-Δ.
(1) Initialization: PXti→Xt-Δroot=1;
(2) Xtemp=Xti;
(3) for k=t,…,t-Δ+1
// the data structure of Xkj includes a pointer to his parent
(4)   Xparent=parent(Xtemp); //a backward smoothing
//operation
(5)   PXti→Xt-Δroot=PXti→Xt-Δroot∗pgyk∣Xtemp∗pfXtemp∣Xparent;
(6)   Xtemp=Xparent;
end
(7) return PXti→Xt-Δroot;Obviously, the time complexity of Procedure1 is OΔ. Figure 5 shows the backward smoothing procedure.Figure 5
The backward smoothing procedure of Procedure1. The calculation of each backward smoothing term is iterated from one leaf node to the root node.


## 4.2. Adding Auxiliary Particles

Experiment results demonstrate that it is not enough to use particles intree(Xti) to calculate the delayed weight w^ki of Xki. We should add some auxiliary particles to tree(Xti) as illustrated in Figure 6.Figure 6
Add auxiliary particles totree(Xti) to get a new tree aux_tree(Xti). Auxiliary particles are denoted by red circles. All leaf nodes of aux_tree(Xti) have the same depth, which equals Δ=3.Definition 6.
One calls the modified generation tree to which auxiliary particles have been added the auxiliary generation tree, denoted byaux_tree(X), where X is the root of aux_tree(X).
In Figure6, we also call these newly added auxiliary particles, like ut+21, ut+31, and ut+32, grandchildren of Xti. The most important feature of aux_tree(Xti) is that all of its leaf nodes have the same depth, which equals the delayed time length Δ.Definition 7.
The state of auxiliary particle is defined as the expectation of transition probability function; that is, the auxiliary particlesut+21 and ut+31 in Figure 6 are defined as follows:(14)ut+21=Epf·∣Xt+1i+2,ut+31=Epf·∣ut+21.
Furthermore, according to the definition of transition probability function (1), we have (15)ut+21=fXt+1i+2+Ewt=fXt+1i+2,ut+31=fut+21+Ewt=fut+21,where since wt~N0,1, Ewt=0.
According toaux_tree(Xti) in Figure 6, the delayed weight w^ki of Xki can be calculated as follows:(16)w^ti∝pXti,Xt+1i+1,Xt+2i+1,ut+31,yt:t+3+pXti,Xt+1i+1,Xt+2i+2,Xt+3i+1,yt:t+3+pXti,Xt+1i+1,Xt+2i+2,Xt+3i+2,yt:t+3+pXti,Xt+1i+2,ut+31,ut+32,yt:t+1+pXti,Xt+1i+3,Xt+2i+3,Xt+3i+3,yt:t+3,where ut+31,Xt+3i+1,Xt+3i+2,ut+32,Xt+3i+3 are leaf nodes of aux_tree(Xti). It is a little different from (11). And, according to the definition of backward smoothing term, w^ti can also be rewritten as(17)w^ti∝Put+31⟶Xti+PXt+3i+1⟶Xti+PXt+3i+2⟶Xti+Put+32⟶Xti+PXt+3i+3⟶Xti.
Absolutely, each backward smoothing term in (17) can be calculated by Procedure 1.The main advantage of adding auxiliary particles to the generation tree is that they can greatly improve the filtering accuracy, which is verified in our experiments.Procedure2 gives the process of constructing an auxiliary generation tree (CAGF).Procedure 2:CAGF.
LetΔ is the height of treeXti, Nk is the number of child/grandchildren of Xti at time k, where t+1≤k≤t+Δ.
Input: treeXti
Output: an auxiliary generation tree aux_treeXti
(1)Initialization: index=0; //index auxiliary particles
(2)if Xti has no child
(3)X=Xti;
(4)for j=k+1,…,t+Δ
(5) ujindex=Epf⋅∣X;
(6) X=ujindex;
(7) index=index+1;
end
end
(8)else for k=t+1,…,t+Δ-1
(9)for n=1,…,Nk
(10)  if Xkn is a leaf node
(11)   X=Xkn;
//add auxiliary particles
(12)   for j=k+1,…,t+Δ
(13)    ujindex=Epf⋅∣X;
(14)    X=ujindex;
(15)    index=index+1;
end
end
end
end
(16) return aux_treeXtiLemma 8.
The number of leaf nodes intreeXti is equal to the number of leaf nodes in aux_treeXti.Proof.
According to step (13) in Procedure 2, only one child is added to every leaf node. So the number of leaf nodes will not increase.Lemma 9.
LetNk denote the number of children/grandchildren of Xti in treeXti at time k, where t+1≤k≤t+Δ. The time complexity TCAGFXti of Procedure 2 satisfies (18)O∑k=t+1t+ΔNk≤TCAGFXti≤OΔ+1∑k=t+1t+ΔNk.The proof is obvious.


## 4.3. Algorithm Framework

In this section, we first give the framework of our algorithm, and then we analyze its time complexity.Theorem 10.
The time complexityTADWPF of Algorithm 3 is OCN, C≪N.Proof.
Obviously, the time cost of Algorithm3 mainly consists of steps (4) and (6). Let TCAGFXt-Δi, 1≤i≤N, denote the time cost of step (4) in ADWPF, and Nki denotes the number of children/grandchildren of Xt-Δi in treeXt-Δi at time k, where t-Δ+1≤k≤t. According to Lemma 9(19)∑i=1NTCAGFXt-Δi≤∑i=1NOΔ+1∑k=t-Δ+1tNki=OΔ+1∑i=1N∑k=t-Δ+1tNki=OΔ+1N,where the second equation is based on the fact that we only sampled N particles at each time step according to step (1) in ADWPF. And according to Lemma 8(20)∑i=1NTCAGFXt-Δi≥∑i=1NO∑k=t-Δ+1tNki=O∑i=1N∑k=t-Δ+1tNki=ON.
LetTBSXt-Δi, 1≤i≤N, denote the time cost of step (6) in ADWPF. The sum number of leaf nodes of all generation trees treeXt-Δ1,…,treeXt-ΔN is less than ΔN and greater than N, so the sum number of leaf nodes of all auxiliary generation trees aux_treeXt-Δ1,…,aux_treeXt-ΔN is less than ΔN and greater than N according to Lemma 8. Since the time cost of Procedure 1 is OΔ for each leaf node, we have (21)N∗OΔ≤∑i=1NTBSXt-Δi≤ΔN∗OΔ=OΔ2N.
So for the time complexityTADWPF of Algorithm 2 we have(22)TADWPF=T∗∑i=1NTCAGFXt-Δi+∑i=1NTBSXt-Δi≤T∗OΔ+1N+OΔ2N=OΔ2+Δ+1TN.
Meanwhile,(23)TADWPF=T∗∑i=1NTCAGFXt-Δi+∑i=1NTBSXt-Δi≥T∗ON+OΔN=OΔ+1TN,where Δ,T≪N. So the time complexity of ADWPF is OCN, C≪N.Algorithm 2:DWFB.
AssumeN is the number of particles, t is current time index, count is the number of simulation steps, Δ is the delayed time interval,
w
t
i is the filtered weight and w^ti is the modified weight.
(1)  for t=1,…,T
(2)  Generate samples Xti,wti,w^ti=0i=1N through a forward particle filter (like Algorithm 1).
(3)  if t>Δ
(4)  for i=1:N
(5)   w^ti=wti
(6)   for k=T-1:T-∆
(7)   for i=1:N
w^ki=∑j=1Nw^k+1jwkipfXk+1j∣Xki∑lNwklpfXk+1j∣Xkl
end
end
end
end
The inference result of system statexT-Δ is
x^t-Δ=∑i=1Nw^t-Δi⋅Xt-ΔiAlgorithm 3:ADWPF.
LetN is the number of particles, t is current time index, Δ is the delayed time length,wti is the forward filtered weight
andw^ti is the delayed weight.
(
1) for t=1,…,T
Generate samples Xti,wti,w^ti=0i=1N through a forward particle filter (like Algorithm 1).
end
(2) if t>Δ
(3)  for i=1,…,N
//Procedure 2.
(4)   aux_treeXt-Δi=CAGFtreeXt-Δi,Δ
(5)   for each leaf node Xtj  of aux_treeXt-Δi
//Procedure 1.
(6)    PXtj→Xt-Δi=BSXtj;
(7)    w^t-Δi=w^t-Δi+PXtj→Xt-Δi;
(8)    X=Xkn;
end
end
end


## 5. Simulation and Results

To verify the performance of our algorithm ADWPF, we consider a wireless sensor network with 70 anchor sensorsai=axi,ayiTi=170. These sensors are randomly distributed in a 50×50 deployment area. Our task is to track a mobile sensor that moves within the area. The state of the mobile sensor at time t is xt=locx,t,vx,t,locy,t,vy,tT, where locx,t,locy,tT is location of the mobile senor and vx,t,vy,tT is its velocity. In our simulation, we set vx,t,vy,tT=1,1T. So we just need to utilize our algorithm to infer the state loct=locx,t,locy,tT of the mobile sensor at time t. The communication radius of sensor is set as 5. At time t, the distance measurement yt from anchor sensor ai is(24)yt=axi-locx,t2+ayi-locy,t2+n,where the distance measurement noise n∈R follows n~N0,1.


### 5.1. Linear Motion Model

The state of the mobile sensor evolves according to [16, 17]:(25)xt=Φxt-1+Γwt,where wt∈R2, wt~N00,0.0001000.0001, and Φ and Γ have the following form:(26)Φ=1100010000110001,Γ=0.501000.501.The initial state of the mobile sensor isx0=0,1,0,1T, the sampling interval is γ=1, the number of tracking steps is T=50, the number of simulations is num=50, and we use the mean squared error (MSE) as the performance measurement. It is defined as follows:(27)MSE=1T∗num∑j=1num∑t=1Tloctj-est_loctj2,where loctj is the true location of the mobile sensor at time t in the jth simulation and est_loctj is the estimation of the location.Table1 compares the average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model.Table 1
Average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model.
MSE
N
=
50
N
=
100
N
=
200
N
=
500
N
=
800
Bootstrap particle filter
0.129
0.081
0.055
0.042
0.039
Auxiliary particle filter
0.112
0.071
0.053
0.040
0.040
DWFBΔ=1
0.123
0.070
0.055
0.042
0.039
DWFBΔ=3
0.107
0.065
0.044
0.030
0.029
DWFBΔ=5
0.106
0.058
0.034
0.025
0.024
ADWPFΔ=1
0.117
0.063
0.047
0.035
0.032
ADWPFΔ=3
0.097
0.053
0.034
0.019
0.017
ADWPFΔ=5
0.098
0.049
0.024
0.015
0.012


### 5.2. Gaussian-Markov Motion Model

In this test, the target node is moving according to the Gaussian-Markov motion model [18]. In this model, we can use one tuning parameter to vary the degree of randomness of the movement, and the target node will not move out of the deployment area [19]. The speed vt and direction dt of the mobile target evolves as follows:(28)vt=αvt-1+1-αv-+1-α2nv,dt=αdt-1+1-αd-+1-α2nd,where nv and nd are white noise, α is the tuning parameter used to control the randomness, and v- and d- are constants representing the mean value of speed and direction. The location loct=locx,t,locy,tT of the mobile target evolves as follows:(29)locx,t=locx,t-1+vtcos⁡dt,locy,t=locy,t-1+vtsin⁡dt.We setα=0.75, v-=1, d-=90°, locx,0=25, locy,0=25, and v0=1. Table 2 compares the average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in Gaussian-Markov motion model.Table 2
Average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in Gaussian-Markov motion model.
MSE
N
=
50
N
=
100
N
=
200
N
=
500
N
=
800
Bootstrap particle filter
2.189
2.038
1.908
1.870
1.859
Auxiliary particle filter
2.194
2.063
1.993
2.021
2.028
DWFBΔ=1
2.577
2.254
2.159
2.103
2.101
DWFBΔ=3
2.598
2.399
2.399
2.312
2.344
DWFBΔ=5
2.029
1.911
1.781
1.763
1.766
ADWPFΔ=1
1.880
1.546
1.441
1.370
1.366
ADWPFΔ=3
1.204
0.999
0.880
0.795
0.781
ADWPFΔ=5
1.178
0.893
0.748
0.606
0.603As illustrated in Tables1 and 2, with N fixed, the MSE values of our method ADWPF with different Δ are all smaller than those of bootstrap particle filter and auxiliary particle filter. Furthermore, with N and Δ fixed, the MSE value of ADWPF is smaller than that of DWFB. Tables 3 and 4 compare the total time cost of bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model (T=50) and Gaussian-Markov motion model (T=50). Since our method ADWPF takes the delay strategy, its time cost is a little higher than bootstrap particle filter and auxiliary particle filter, but our method efficiently utilizes the extra time cost to improve tracking accuracy, which is verified in Tables 1 and 2. Comparing to the delay-weight method DWFB in [15], our method ADWPF takes full advantage of the structure relationship among particles, so the time cost of ADWPF is far less than that of DWFB. These experiment results demonstrate the excellent performance of our algorithm.Table 3
Total time (sec.) cost of bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model (T=50).
MSE
N
=
50
N
=
100
N
=
200
N
=
500
N
=
800
Bootstrap particle filter
0.017
0.028
0.056
0.141
0.230
Auxiliary particle filter
0.031
0.062
0.121
0.301
0.497
DWFBΔ=1
0.246
0.927
3.713
22.322
56.934
DWFBΔ=3
0.663
2.592
10.494
63.828
164.757
DWFBΔ=5
1.040
4.096
16.320
101.648
259.860
ADWPFΔ=1
0.042
0.066
0.124
0.312
0.512
ADWPFΔ=3
0.084
0.153
0.317
0.803
1.300
ADWPFΔ=5
0.154
0.308
0.610
1.520
2.458Table 4
Total time (sec.) cost of bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in Gaussian-Markov motion model (T=50).
MSE
N
=
50
N
=
100
N
=
200
N
=
500
N
=
800
Bootstrap particle filter
0.011
0.021
0.037
0.101
0.161
Auxiliary particle filter
0.021
0.038
0.069
0.184
0.301
DWFBΔ=1
0.271
0.925
3.654
22.719
58.071
DWFBΔ=3
0.679
2.632
10.448
65.003
166.134
DWFBΔ=5
1.051
4.167
16.664
103.563
263.259
ADWPFΔ=1
0.044
0.055
0.109
0.266
0.444
ADWPFΔ=3
0.087
0.147
0.287
0.709
1.135
ADWPFΔ=5
0.142
0.276
0.547
1.764
2.179


## 5.1. Linear Motion Model

The state of the mobile sensor evolves according to [16, 17]:(25)xt=Φxt-1+Γwt,where wt∈R2, wt~N00,0.0001000.0001, and Φ and Γ have the following form:(26)Φ=1100010000110001,Γ=0.501000.501.The initial state of the mobile sensor isx0=0,1,0,1T, the sampling interval is γ=1, the number of tracking steps is T=50, the number of simulations is num=50, and we use the mean squared error (MSE) as the performance measurement. It is defined as follows:(27)MSE=1T∗num∑j=1num∑t=1Tloctj-est_loctj2,where loctj is the true location of the mobile sensor at time t in the jth simulation and est_loctj is the estimation of the location.Table1 compares the average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model.Table 1
Average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model.
MSE
N
=
50
N
=
100
N
=
200
N
=
500
N
=
800
Bootstrap particle filter
0.129
0.081
0.055
0.042
0.039
Auxiliary particle filter
0.112
0.071
0.053
0.040
0.040
DWFBΔ=1
0.123
0.070
0.055
0.042
0.039
DWFBΔ=3
0.107
0.065
0.044
0.030
0.029
DWFBΔ=5
0.106
0.058
0.034
0.025
0.024
ADWPFΔ=1
0.117
0.063
0.047
0.035
0.032
ADWPFΔ=3
0.097
0.053
0.034
0.019
0.017
ADWPFΔ=5
0.098
0.049
0.024
0.015
0.012


## 5.2. Gaussian-Markov Motion Model

In this test, the target node is moving according to the Gaussian-Markov motion model [18]. In this model, we can use one tuning parameter to vary the degree of randomness of the movement, and the target node will not move out of the deployment area [19]. The speed vt and direction dt of the mobile target evolves as follows:(28)vt=αvt-1+1-αv-+1-α2nv,dt=αdt-1+1-αd-+1-α2nd,where nv and nd are white noise, α is the tuning parameter used to control the randomness, and v- and d- are constants representing the mean value of speed and direction. The location loct=locx,t,locy,tT of the mobile target evolves as follows:(29)locx,t=locx,t-1+vtcos⁡dt,locy,t=locy,t-1+vtsin⁡dt.We setα=0.75, v-=1, d-=90°, locx,0=25, locy,0=25, and v0=1. Table 2 compares the average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in Gaussian-Markov motion model.Table 2
Average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in Gaussian-Markov motion model.
MSE
N
=
50
N
=
100
N
=
200
N
=
500
N
=
800
Bootstrap particle filter
2.189
2.038
1.908
1.870
1.859
Auxiliary particle filter
2.194
2.063
1.993
2.021
2.028
DWFBΔ=1
2.577
2.254
2.159
2.103
2.101
DWFBΔ=3
2.598
2.399
2.399
2.312
2.344
DWFBΔ=5
2.029
1.911
1.781
1.763
1.766
ADWPFΔ=1
1.880
1.546
1.441
1.370
1.366
ADWPFΔ=3
1.204
0.999
0.880
0.795
0.781
ADWPFΔ=5
1.178
0.893
0.748
0.606
0.603As illustrated in Tables1 and 2, with N fixed, the MSE values of our method ADWPF with different Δ are all smaller than those of bootstrap particle filter and auxiliary particle filter. Furthermore, with N and Δ fixed, the MSE value of ADWPF is smaller than that of DWFB. Tables 3 and 4 compare the total time cost of bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model (T=50) and Gaussian-Markov motion model (T=50). Since our method ADWPF takes the delay strategy, its time cost is a little higher than bootstrap particle filter and auxiliary particle filter, but our method efficiently utilizes the extra time cost to improve tracking accuracy, which is verified in Tables 1 and 2. Comparing to the delay-weight method DWFB in [15], our method ADWPF takes full advantage of the structure relationship among particles, so the time cost of ADWPF is far less than that of DWFB. These experiment results demonstrate the excellent performance of our algorithm.Table 3
Total time (sec.) cost of bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model (T=50).
MSE
N
=
50
N
=
100
N
=
200
N
=
500
N
=
800
Bootstrap particle filter
0.017
0.028
0.056
0.141
0.230
Auxiliary particle filter
0.031
0.062
0.121
0.301
0.497
DWFBΔ=1
0.246
0.927
3.713
22.322
56.934
DWFBΔ=3
0.663
2.592
10.494
63.828
164.757
DWFBΔ=5
1.040
4.096
16.320
101.648
259.860
ADWPFΔ=1
0.042
0.066
0.124
0.312
0.512
ADWPFΔ=3
0.084
0.153
0.317
0.803
1.300
ADWPFΔ=5
0.154
0.308
0.610
1.520
2.458Table 4
Total time (sec.) cost of bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in Gaussian-Markov motion model (T=50).
MSE
N
=
50
N
=
100
N
=
200
N
=
500
N
=
800
Bootstrap particle filter
0.011
0.021
0.037
0.101
0.161
Auxiliary particle filter
0.021
0.038
0.069
0.184
0.301
DWFBΔ=1
0.271
0.925
3.654
22.719
58.071
DWFBΔ=3
0.679
2.632
10.448
65.003
166.134
DWFBΔ=5
1.051
4.167
16.664
103.563
263.259
ADWPFΔ=1
0.044
0.055
0.109
0.266
0.444
ADWPFΔ=3
0.087
0.147
0.287
0.709
1.135
ADWPFΔ=5
0.142
0.276
0.547
1.764
2.179


## 6. Conclusion

Particle filter is a very important technology for tracking mobile targets in wireless sensor networks. In order to improve the tracking accuracy of the mobile node, we propose an advanced auxiliary delay-weight particle filter algorithm. By adding some auxiliary particles to the generation trees formed by sampled particles, our algorithm ADWPF can efficiently utilize measurement information after timet to reweight Monte Carlo samples Xtii=1N generated at time t without creating a substantially high computation burden. Experiment results demonstrate that our algorithm ADWPF can greatly improve the tracking accuracy of a mobile node, compared with bootstrap filter, auxiliary particle filter, and another delayed-weight particle filter in [15].



---

*Source: 4535963-2016-05-11.xml*