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

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.


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][2][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][5][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 treelike 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 is O(Δ ); here Δ is the delayed time length, is the number of tracking steps, and is the number of samples in each step.
The remainder of this paper is organized as follows. Section 2 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.

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(Δ 2 ) in [13], where Δ is the time interval and is the number of samples in each step.

Particle
Filter. The aim of particle filter is to get the sequential Monte Carlo samples 1: ∼ ( 1: | 1: ), here = 1, 2, . . . , , so the Monte Carlo approximation of ( 1: | where (⋅) is Delta function. The particle filter generally consists of the steps shown in Algorithm 1 (take bootstrap filter algorithm as an example). The time complexity is O( ) at each time step, and memory complexity is O( ).

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 time + Δ takes the form 1: +Δ = ( 1 , . . . , +Δ ), the true system states are 1: = ( 1 , . . . , ), and ( 1: | 1: +Δ ) is the expectation of 1: based on 1: +Δ . Let̂+ Δ be a consistent Monte Carlo of ( 1: | 1: +Δ ). The mean squared difference between̂+ Δ and 1: is (̂+ Δ − 1: ) 2 that follows. 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 Proposition 1 [10]. The first method is the delayed-sample method. It designed an importance distribution, which utilizes measurement information after time to regenerate weighted samples of state at time , for example, auxiliary particle filter. The second method is the delayed-weight method, which utilized information after time to reweight samples { } =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]: 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).

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. Figure 1 illustrates the evolvement of particles in the standard particle filter method at time .
Therefore, the delayed weight̂of particle satisfieŝ 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 + 1 4 International Journal of Distributed Sensor Networks · · · · · · · · · · · · · · · · · · y 1 y 2 y t+Δ−1 y t+Δ to calculate the delayed weight of each particle at time ; for example, the delayed weight̂of iŝ Since the number of particles is at time , the time complexity of this algorithm is (Δ 2 ); when 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 (Δ ) but also greatly improves the tracking accuracy.

An Auxiliary Delayed-Weight Method
Our method consists of two key technologies: designing a backward smoothing procedure and adding auxiliary particles.
Observing the evolvement of particles in Figure 1, 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 particle and its children/grandchildren is called a generation tree, denoted by tree( ). is the root of tree( ).
As illustrated in Figure 3, the particles denoted by green circles construct the generation tree tree( ) and is the root of tree( ).
The main advantage of doing so is that we only need to make use of children/grandchildren of to calculate the delayed weight̂of , shown in Figure 4. Therefore, (8) can be reduced tô where + is the number of children/grandchildren of at time + , = 1, . . . , Δ.
Definition 5. One calls the joint probability function representing the parent-child relation between leaf node leaf and root node the backward smoothing term, denoted by ( leaf → ); for example, (X i+2 t+1 → ) = ( , X i+2 t+1 , : +1 ). Thereforêin (11) can be rewritten aŝ More importantly, each backward smoothing term can be decomposed as the product of transition probability functions and likelihood functions; for example, the backward International Journal of Distributed Sensor Networks y t+1 y t+2 y t+3 Figure 3: The particles denoted by green circles form tree( ); is the root of tree( ). smoothing term (X i+1 t+2 → ) can be decomposed as follows: 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.
Obviously, the time complexity of Procedure 1 is (Δ). Figure 5 shows the backward smoothing procedure.

Adding Auxiliary Particles.
Experiment results demonstrate that it is not enough to use particles in tree( ) to calculate the delayed weight̂of . We should add some auxiliary particles to tree( ) as illustrated in Figure 6. Definition 6. One calls the modified generation tree to which auxiliary particles have been added the auxiliary generation tree, denoted by aux tree( ), where is the root of aux tree( ).
In Figure 6, we also call these newly added auxiliary particles, like 1 +2 , 1 +3 , and 2 +3 , grandchildren of . The most important feature of aux tree( ) 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 particles 1 +2 and 1 +3 in Figure 6 are defined as follows: 6 International Journal of Distributed Sensor Networks Input: a leaf node Output: the backward smoothing term ( → root −Δ ), where root −Δ is the grandparent of at time − Δ. (1) Initialization: ( → root −Δ ) = 1; (2) temp = ; (3) for = , . . . , − Δ + 1 // the data structure of includes a pointer to his parent (4) parent = parent( temp ); //a backward smoothing //operation Procedure 1: BS. Furthermore, according to the definition of transition probability function (1), we have where since ∼ N(0, 1), ( ) = 0. According to aux tree( ) in Figure 6, the delayed weight of can be calculated as follows: where {u 1 t+3 , X i+1 t+3 , X i+2 t+3 , u 2 t+3 , X i+3 t+3 } are leaf nodes of aux tree( ). It is a little different from (11). And, according to the definition of backward smoothing term,̂can also be rewritten aŝ 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.
Procedure 2 gives the process of constructing an auxiliary generation tree (CAGF).

Lemma 8. The number of leaf nodes in tree( ) is equal to the number of leaf nodes in aux tree( ).
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. ) . (18) The proof is obvious.

Algorithm Framework.
In this section, we first give the framework of our algorithm, and then we analyze its time complexity.
tree(X i t ) aux_tree(X i t ) Figure 6: Add auxiliary particles to tree( ) to get a new tree aux tree( ). Auxiliary particles are denoted by red circles. All leaf nodes of aux tree( ) have the same depth, which equals Δ = 3.

Simulation and Results
To verify the performance of our algorithm ADWPF, we consider a wireless sensor network with 70 anchor sensors { = ( , ) } 70 =1 . 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 is = (loc , , V , , loc , , V , ) , where (loc , , loc , ) is location of the mobile senor and (V , , V , ) is its velocity. In our simulation, we set (V , , V , ) = (1, 1) . So we just need to utilize our algorithm to infer the state loc = (loc , , loc , ) of the mobile sensor at time . The communication radius of sensor is set as 5. At time , the distance measurement from anchor sensor is where the distance measurement noise ∈ follows ∼ N(0, 1).

Linear Motion
Model. The state of the mobile sensor evolves according to [16,17]: where ] .

(26)
The initial state of the mobile sensor is 0 = (0, 1, 0, 1) , the sampling interval is = 1, the number of tracking steps is = 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: where loc is the true location of the mobile sensor at time in the th simulation and est loc is the estimation of the location. Table 1 compares the average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model.

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 V and direction of the mobile target evolves as follows: where V and are white noise, is the tuning parameter used to control the randomness, and V and are constants Assume is the number of particles, is current time index, count is the number of simulation steps, Δ is the delayed time interval, is the filtered weight and̂is the modified weight. (1) for = 1, . . . , (2) Generate samples { , ,̂= 0} =1 through a forward particle filter (like Algorithm 1).
Let is the number of particles, is current time index, Δ is the delayed time length, is the forward filtered weight and̂is the delayed weight.     representing the mean value of speed and direction. The location loc = (loc , , loc , ) of the mobile target evolves as follows: loc , = loc , −1 + V cos , loc , = loc , −1 + V sin .
As illustrated in Tables 1 and 2, with 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 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 ( = 50) and Gaussian-Markov motion model ( = 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 delayweight 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.

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 time to reweight Monte Carlo samples { } =1 generated at time 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].