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

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.

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 is $\mathcal{O}\left(\Delta TN\right)$ ; 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. 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.

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 $\mathcal{O}\left({\Delta N}^2\right)$ 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:

$$x_{t+\mathrm{1}}=f\left(x_t\right)+w_t,$$

(1)

$$y_{t+1}=g\left(x_{t+1}\right)+v_t,$$

(2)

where $x_t$ is state, $y_{t+1}$ is measurement, and $w_t,v_t~\mathcal{N}\left(\mathrm{0,1}\right)$ . Therefore, the transition probability function of state variable $x_{t+1}$ is $x_{t+1}~p_f(·\mid x_t)$ and the likelihood function of measurement $y_{t+1}$ is $y_{t+1}~p_g(·\mid x_{t+1})$ . The joint distribution of measurement information $y_{1:t}=\left(y_1,\dots ,y_t\right)$ and state sequence $x_{1:t}=\left(x_1,\dots ,x_t\right)$ is

$$\begin{array}{cc}& p\left(x_{\mathrm{1}:t},y_{\mathrm{1}:t}\right)=\pi \left(x_{\mathrm{1}}\right){\displaystyle \prod }_{k=\mathrm{1}}^t{p}_f\left(x_{k+\mathrm{1}}\mid x_k\right){\displaystyle \prod }_{k=\mathrm{1}}^t{p}_g\left(y_k\mid x_k\right),\end{array}$$

(3)

where $\pi \left(x_1\right)$ is the prior probability distribution of initial state. So the posterior distribution of $x_{1:t}$ is

$$\begin{array}{cc}& p\left(x_{\mathrm{1}:t}\mid y_{\mathrm{1}:t}\right)=\frac{p\left(x_{\mathrm{1}:t},y_{\mathrm{1}:t}\right)}{p\left(y_{\mathrm{1}:t}\right)}.\end{array}$$

(4)

3.2. Particle Filter

The aim of particle filter is to get the sequential Monte Carlo samples $X_{1:t}^i~p\left(x_{1:t}\mid y_{1:t}\right)$ , here $i=\mathrm{1,2},\dots ,N$ , so the Monte Carlo approximation of $p\left(x_{1:t}\mid y_{1:t}\right)$ is

$$\begin{array}{cc}& \widehat{p}\left(x_{\mathrm{1}:t}\mid y_{\mathrm{1}:t}\right)=\frac{\mathrm{1}}{N}{\displaystyle \sum }_{i=\mathrm{1}}^N{\delta }_{X_{\mathrm{1}:t}^i}\left(x_{\mathrm{1}:t}\right),\end{array}$$

(5)

where $\delta \left(·\right)$ 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=\mathrm{1}$ , sample ${\tilde{X}}_{\mathit{1}}^i~\pi \left(x_{\mathit{1}}\right)$ to obtain the Monte Carlo approximation of $p\left(x_{\mathit{1}}\mid y_{\mathit{1}}\right)$

             $$\tilde{p}\left(x_{\mathit{1}}\mid y_{\mathit{1}}\right)={\displaystyle \sum }_{i=\mathit{1}}^N{W}_{\mathit{1}}^i{\delta }_{{\tilde{X}}_{\mathit{1}}^i}\left(x_{\mathit{1}}\right)$$ ,

where $W_{\mathit{1}}^i\propto p_g\left(y_{\mathit{1}}\mid {\tilde{X}}_{\mathit{1}}^i\right)$ is the weight of ${\tilde{X}}_{\mathit{1}}^i$ , resample $X_{\mathit{1}}^i~\tilde{p}\left(x_{\mathit{1}}\mid y_{\mathit{1}}\right)$ to obtain

             $$\widehat{p}\left(x_{\mathit{1}}\mid y_{\mathit{1}}\right)=\frac{\mathit{1}}{N}{\displaystyle \sum }_{i=\mathit{1}}^N{\delta }_{X_{\mathit{1}}^i}\left(x_{\mathit{1}}\right)$$ .

Step  2 (Filtering). At time $t\ge \mathrm{2}$ , sample ${\tilde{X}}_t^i~p_f\left(\cdot \mid X_{t-\mathit{1}}^i\right)$ to obtain the Monte Carlo approximation of $p\left(x_t\mid y_{\mathit{1}:t}\right)$

             $$\tilde{p}\left(x_t\mid y_{\mathit{1}:t}\right)={\displaystyle \sum }_{i=\mathit{1}}^N{W}_t^i{\delta }_{{\tilde{X}}_t^i}\left(x_t\right)$$ ,

where $W_t^i\propto p_g\left(y_t\mid {\tilde{X}}_t^i\right)$ is the weight of ${\tilde{X}}_t^i$ , resample $X_t^i~\tilde{p}\left(x_t\mid y_{\mathit{1}:t}\right)$ to obtain

             $$\widehat{p}\left(x_t\mid y_{\mathit{1}:t}\right)=\frac{\mathit{1}}{N}{\displaystyle \sum }_{i=\mathit{1}}^N{\delta }_{X_t^i}\left(x_t\right)$$

The time complexity is $\mathcal{O}\left(TN\right)$ at each time step, and memory complexity is $\mathcal{O}\left(TN\right)$ .

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 time $t+\Delta$ takes the form $y_{1:t+\Delta }=\left(y_1,\dots ,y_{t+\Delta }\right)$ , the true system states are $x_{1:t}=\left(x_1,\dots ,x_t\right)$ , and $E\left(x_{1:t}\mid y_{1:t+\Delta }\right)$ is the expectation of $x_{1:t}$ based on $y_{1:t+\Delta }$ . Let ${\widehat{X}}_{t+\Delta }$ be a consistent Monte Carlo of $E\left(x_{1:t}\mid y_{1:t+\Delta }\right)$ . The mean squared difference between ${\widehat{X}}_{t+\Delta }$ and $x_{1:t}$ is $E{\left({\widehat{X}}_{t+\Delta }-x_{1:t}\right)}^2$ that follows.

Proposition 1.

When the Monte Carlo sample size N is sufficiently large, $E{\left({\widehat{X}}_{t+\Delta }-x_{1:t}\right)}^2$ 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 Proposition 1 [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 ${\left\{X_t^i\right\}}_{i=1}^N$ .

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]:

$$\begin{array}{cc}& {\widehat{X}}_t^i~q_t\left(x_t\mid X_{t-\mathrm{1}}^i,y_{t:t+\Delta }\right)\propto \int p\left(X_{t-1}^i,x_t,x_{t+1},\dots ,x_{t+\Delta },y_{t:t+\Delta }\right)d{x}_{t+1}\cdots d{x}_{t+\Delta }=p\left(X_{t-1}^i,x_t,y_{t:t+\Delta }\right).\end{array}$$

(6)

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. Figure 1 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 time t.

Assume that we observe measurements $\left(y_1,\dots ,y_{t+\Delta }\right)$ and have got weighted particles ${\left\{X_1^i,w_1^i\right\}}_{i=1}^N,\dots ,{\left\{X_{t+\Delta }^i,w_{t+\Delta }^i\right\}}_{i=1}^N$ using a forward particle filter method until time $t+\Delta$ ; 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 distribution $p\left(X_{t-1}^i,x_t,y_{t:t+\Delta }\right)$ in (6) can be approximated by

$$\begin{array}{cc}& p\left(X_{t-\mathrm{1}}^i,x_t,y_{t:t+\Delta }\right)\approx {\displaystyle \sum }_{i_{t+1}=1}^N\cdots {\displaystyle \sum }_{i_{t+\Delta }=1}^{N}p\left(X_{t-1}^i,x_t,X_{t+1}^{i_{t+1}},\dots ,X_{t+\Delta }^{i_{t+\Delta }},y_{t:t+\Delta }\right).\end{array}$$

(7)

Therefore, the delayed weight ${\widehat{w}}_t^i$ of particle $X_t^i$ satisfies

$$\begin{array}{cc}& {\widehat{w}}_t^i\propto p\left(X_{t-\mathrm{1}}^i,X_t^i,y_{t:t+\Delta }\right)={\displaystyle \sum }_{i_{t+1}=1}^N\cdots {\displaystyle \sum }_{i_{t+\Delta }=1}^{N}p\left(X_{t-1}^i,X_t^i,X_{t+1}^{i_{t+1}},\dots ,X_{t+\Delta }^{i_{t+\Delta }},y_{t:t+\Delta }\right).\end{array}$$

(8)

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+\mathrm{1}$ to calculate the delayed weight of each particle at time k; for example, the delayed weight ${\widehat{w}}_k^i$ of $X_k^i$ is

$$\begin{array}{cc}& {\widehat{w}}_k^i\propto {\displaystyle \sum }_{j=\mathrm{1}}^N{\widehat{w}}_{k+\mathrm{1}}^j{p}_f\left(X_{k+\mathrm{1}}^j\mid X_k^i\right),\hspace{1em}i=\mathrm{1},\dots ,N.\end{array}$$

(9)

Since the number of particles is N at time k, the time complexity of this algorithm is $O\left(\Delta T{N}^2\right)$ ; 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\left(\Delta TN\right)$ 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 particle $X_t$ at time t is sampled from the importance distribution $q\left(x_t\mid X_{t-1}\right)$ , where $X_{t-1}$ is a particle at time $t-1$ and $t=1,\dots ,T$ , one says that $X_t$ is the child of $X_{t-1}$ , denoted by $X_t=\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{l}\mathrm{d}\left(X_{t-1}\right)$ or, conversely, that $X_{t-1}$ is the parent of $X_t$ , denoted by $X_{t-1}=\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}(X_t)$ . Furthermore, if $X_t$ has a child $X_{t+1}$ , then $X_{t+1}$ is the grandchild of $X_{t-1}$ , denoted by $X_{t+1}=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{l}\mathrm{d}\left(X_{t-1}\right)$ , and $X_{t-1}$ is the grandparent of $X_{t+1}$ , denoted by $X_{t-1}=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}\left(X_{t+1}\right)$ .

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 X and its children/grandchildren is called a generation tree, denoted by $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X)$ . X is the root of $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X)$ .

As illustrated in Figure 3, the particles denoted by green circles construct the generation tree $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ and $X_t^i$ is the root of $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ .

Figure 3

The particles denoted by green circles form $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ ; $X_t^i$ is the root of $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ .

Assumption 4.

If the particle $X_{t+\Delta }$ is not the child/grandchild of $X_t$ , then $p_f\left(X_{t+k}\mid X_t\right)=\mathrm{0}$ , where $p_f\left(·\right)$ is the transition probability function, $t=1,\dots ,T$ , and $k>0$ .

The main advantage of doing so is that we only need to make use of children/grandchildren of $X_t^i$ to calculate the delayed weight ${\widehat{w}}_k^i$ of $X_k^i$ , shown in Figure 4. Therefore, (8) can be reduced to

$$\begin{array}{cc}& {\widehat{w}}_t^i\propto p\left(X_{t-\mathrm{1}}^i,X_t^i,y_{t:t+\Delta }\right)\approx {\displaystyle \sum }_{i_{t+1}=1}^{N_{t+1}}\cdots {\displaystyle \sum }_{i_{t+\Delta }=1}^{N_{t+\Delta }}p\left(X_{t-1}^i,X_t^i,X_{t+1}^{i_{t+1}},\dots ,X_{t+\Delta }^{i_{t+\Delta }},y_{t:t+\Delta }\right),\end{array}$$

(10)

where $N_{t+j}$ is the number of children/grandchildren of $X_t^i$ at time $t+j$ , $j=1,\dots ,\Delta$ .

Figure 4

We only need to use children/grandchildren of $X_t^i$ , which are denoted by green circles, to calculate the delayed weight ${\widehat{w}}_k^i$ of $X_k^i$ .

Observing the structure of $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ in Figure 4(b), we design an advanced delayed-weight algorithm whose time complexity is $O\left(\Delta TN\right)$ to calculate ${\widehat{w}}_k^i$ . Take $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ in Figure 4(b) as an example; here $\Delta =3$ . Since there are 5 leaf nodes in $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ , the delayed weight ${\widehat{w}}_k^i$ of $X_k^i$ can be calculated as follows:

$$\begin{array}{cc}& {\widehat{w}}_t^i\propto p\left(X_t^i,{\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{1}}^{\mathit{i}\mathbf{+}\mathbf{2}},y_{t:t+\mathrm{1}}\right)+p\left(X_t^i,X_{t+1}^{i+1},{\mathit{X}}_{\mathit{t}+\mathit{2}}^{\mathit{i}+\mathit{1}},y_{t:t+2}\right)+p\left(X_t^i,X_{t+1}^{i+1},X_{t+2}^{i+2},{\mathit{X}}_{\mathit{t}+\mathit{3}}^{\mathit{i}+\mathit{1}},y_{t:t+3}\right)+p\left(X_t^i,X_{t+1}^{i+1},X_{t+2}^{i+2},{\mathit{X}}_{\mathit{t}+\mathit{3}}^{\mathit{i}+\mathit{2}},y_{t:t+3}\right)+p\left(X_t^i,X_{t+1}^{i+3},X_{t+2}^{i+3},{\mathit{X}}_{\mathit{t}+\mathit{3}}^{\mathit{i}+\mathit{3}},y_{t:t+3}\right),\end{array}$$

(11)

where $\left\{X_{t+2}^{i+1},X_{t+3}^{i+1},X_{t+3}^{i+2},X_{t+1}^{i+2},X_{t+3}^{i+3}\right\}$ are leaf nodes of $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ . 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 node $X_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}$ and root node X the backward smoothing term, denoted by $P\left(X_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}\to X\right)$ ; for example, $P\left({\mathit{X}}_{\mathit{t}+\mathit{1}}^{\mathit{i}+\mathit{2}}\to X_t^i\right)=p\left(X_t^i,{\mathit{X}}_{\mathit{t}+\mathit{1}}^{\mathit{i}+\mathit{2}},y_{t:t+1}\right)$ .

Therefore ${\widehat{w}}_t^i$ in (11) can be rewritten as

$$\begin{array}{cc}& {\widehat{w}}_t^i\propto P\left({\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{1}}^{\mathit{i}\mathbf{+}\mathbf{2}}⟶X_t^i\right)+P\left({\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{2}}^{\mathit{i}\mathbf{+}\mathbf{1}}⟶X_t^i\right)+P\left({\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathit{i}\mathbf{+}\mathbf{1}}⟶X_t^i\right)+P\left({\mathit{X}}_{\mathit{t}+\mathit{3}}^{\mathit{i}+\mathit{2}}⟶X_t^i\right)+P\left({\mathit{X}}_{\mathit{t}+\mathit{3}}^{\mathit{i}+\mathit{3}}⟶X_t^i\right).\end{array}$$

(12)

More importantly, each backward smoothing term can be decomposed as the product of transition probability functions and likelihood functions; for example, the backward smoothing term $P\left({\mathit{X}}_{\mathit{t}+\mathit{2}}^{\mathit{i}+\mathit{1}}\to X_t^i\right)$ can be decomposed as follows:

$$\begin{array}{cc}& P\left({\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{2}}^{\mathit{i}\mathbf{+}\mathbf{1}}⟶X_t^i\right)=p\left(X_t^i,X_{t+\mathrm{1}}^{i+\mathrm{1}},{\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{2}}^{\mathit{i}\mathbf{+}\mathbf{1}},y_{t:t+\mathrm{2}}\right)=p\left(X_t^i\right)p_f\left(X_{t+1}^{i+1}\mid X_t^i\right)p_f\left(X_{t+2}^{i+1}\mid X_{t+1}^{i+1}\right)\ast p_g\left(y_t\mid X_t^i\right)p_g\left(y_{t+1}\mid X_{t+1}^{i+1}\right)p_g\left(y_{t+2}\mid X_{t+2}^{i+1}\right).\end{array}$$

(13)

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 $X_t^i$

Output: the backward smoothing term $P\left(X_t^i\to X_{t-\mathrm{\Delta }}^{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t}}\right)$ , where $X_{t-\mathrm{\Delta }}^{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t}}$ is the grandparent of $X_t^i$ at time $t-\mathrm{\Delta }$ .

( $\mathit{1}$ )  Initialization: $P\left(X_t^i\to X_{t-\mathrm{\Delta }}^{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t}}\right)=\mathit{1}$ ;

( $\mathit{2}$ )   $X_{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}}=X_t^i$ ;

( $\mathit{3}$ )  for $k=t,\dots ,t-\mathrm{\Delta }+\mathit{1}$

    // the data structure of $X_k^j$ includes a pointer to his parent

( $\mathit{4}$ )     $X_{\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}}=\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}(X_{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}})$ ; //a backward smoothing

              //operation

( $\mathit{5}$ )     $P\left(X_t^i\to X_{t-\mathrm{\Delta }}^{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t}}\right)=P\left(X_t^i\to X_{t-\mathrm{\Delta }}^{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t}}\right)\ast p_g\left(y_k\mid X_{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}}\right)\ast p_f\left(X_{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}}\mid X_{\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}}\right)$ ;

( $\mathit{6}$ )     $X_{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}}=X_{\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}}$ ;

    end

( $\mathit{7}$ )  return $P\left(X_t^i\to X_{t-\mathrm{\Delta }}^{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t}}\right)$ ;

Obviously, the time complexity of Procedure 1 is $O\left(\Delta \right)$ . Figure 5 shows the backward smoothing procedure.

Figure 5

The backward smoothing procedure of Procedure 1. 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 in $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ to calculate the delayed weight ${\widehat{w}}_k^i$ of $X_k^i$ . We should add some auxiliary particles to $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ as illustrated in Figure 6.

Figure 6

Add auxiliary particles to $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ to get a new tree $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ . Auxiliary particles are denoted by red circles. All leaf nodes of $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ have the same depth, which equals $\Delta =3$ .

Definition 6.

One calls the modified generation tree to which auxiliary particles have been added the auxiliary generation tree, denoted by $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X)$ , where X is the root of $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X)$ .

In Figure 6, we also call these newly added auxiliary particles, like $u_{t+2}^1$ , $u_{t+3}^1$ , and $u_{t+3}^2$ , grandchildren of $X_t^i$ . The most important feature of $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ 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 $u_{t+2}^1$ and $u_{t+3}^1$ in Figure 6 are defined as follows:

$$\begin{array}{cc}& u_{t+\mathrm{2}}^{\mathrm{1}}=E_{p_f}\left(·\mid X_{t+\mathrm{1}}^{i+\mathrm{2}}\right),\\ u_{t+3}^1=E_{p_f}\left(·\mid u_{t+2}^1\right).\end{array}$$

(14)

Furthermore, according to the definition of transition probability function (1), we have

$$\begin{array}{cc}& u_{t+\mathrm{2}}^{\mathrm{1}}=f\left(X_{t+\mathrm{1}}^{i+\mathrm{2}}\right)+E\left(w_t\right)=f\left(X_{t+\mathrm{1}}^{i+\mathrm{2}}\right),\\ u_{t+3}^1=f\left(u_{t+2}^1\right)+E\left(w_t\right)=f\left(u_{t+2}^1\right),\end{array}$$

(15)

where since $w_t~\mathcal{N}\left(\mathrm{0,1}\right)$ , $E\left(w_t\right)=0$ .

According to $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ in Figure 6, the delayed weight ${\widehat{w}}_k^i$ of $X_k^i$ can be calculated as follows:

$$\begin{array}{cc}& {\widehat{w}}_t^i\propto p\left(X_t^i,X_{t+\mathrm{1}}^{i+\mathrm{1}},X_{t+\mathrm{2}}^{i+\mathrm{1}},{\mathit{u}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathbf{1}},y_{t:t+\mathrm{3}}\right)+p\left(X_t^i,X_{t+\mathrm{1}}^{i+\mathrm{1}},X_{t+\mathrm{2}}^{i+\mathrm{2}},{\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathit{i}\mathbf{+}\mathbf{1}},y_{t:t+\mathrm{3}}\right)+p\left(X_t^i,X_{t+\mathrm{1}}^{i+\mathrm{1}},X_{t+\mathrm{2}}^{i+\mathrm{2}},{\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathit{i}\mathbf{+}\mathbf{2}},y_{t:t+\mathrm{3}}\right)+p\left(X_t^i,X_{t+\mathrm{1}}^{i+\mathrm{2}},u_{t+\mathrm{3}}^{\mathrm{1}},{\mathit{u}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathbf{2}},y_{t:t+\mathrm{1}}\right)+p\left(X_t^i,X_{t+\mathrm{1}}^{i+\mathrm{3}},X_{t+\mathrm{2}}^{i+\mathrm{3}},{\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathit{i}\mathbf{+}\mathbf{3}},y_{t:t+\mathrm{3}}\right),\end{array}$$

(16)

where $\left\{{\mathit{u}}_{\mathit{t}+\mathit{3}}^{\mathit{1}},{\mathit{X}}_{\mathit{t}+\mathit{3}}^{\mathit{i}+\mathit{1}},{\mathit{X}}_{\mathit{t}+\mathit{3}}^{\mathit{i}+\mathit{2}},{\mathit{u}}_{\mathit{t}+\mathit{3}}^{\mathit{2}},{\mathit{X}}_{\mathit{t}+\mathit{3}}^{\mathit{i}+\mathit{3}}\right\}$ are leaf nodes of $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}(X_t^i)$ . It is a little different from (11). And, according to the definition of backward smoothing term, ${\widehat{w}}_t^i$ can also be rewritten as

$$\begin{array}{cc}& {\widehat{w}}_t^i\propto P\left({\mathit{u}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathbf{1}}⟶X_t^i\right)+P\left({\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathit{i}\mathbf{+}\mathbf{1}}⟶X_t^i\right)+P\left({\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathit{i}\mathbf{+}\mathbf{2}}⟶X_t^i\right)+P\left({\mathit{u}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathbf{2}}⟶X_t^i\right)+P\left({\mathit{X}}_{\mathit{t}\mathbf{+}\mathbf{3}}^{\mathit{i}\mathbf{+}\mathbf{3}}⟶X_t^i\right).\end{array}$$

(17)

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).

Procedure 2: CAGF.

Let $\mathrm{\Delta }$ is the height of $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_t^i\right)$ , $N_k$ is the number of child/grandchildren of $X_t^i$ at time k, where $t+\mathrm{1}\le k\le t+\mathrm{\Delta }$ .

Input: $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_t^i\right)$

Output: an auxiliary generation tree $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_t^i\right)$

( $\mathit{1}$ ) Initialization: $\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}=\mathit{0}$ ; //index auxiliary particles

( $\mathit{2}$ ) if $X_t^i$ has no child

( $\mathit{3}$ )   $X=X_t^i$ ;

( $\mathit{4}$ )  for $j=k+\mathit{1},\dots ,t+\mathrm{\Delta }$

( $\mathit{5}$ )    $u_j^{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}}=E_{p_f}\left(\cdot \mid X\right)$ ;

( $\mathit{6}$ )    $X=u_j^{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}}$ ;

( $\mathit{7}$ )    $\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}=\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}+\mathit{1}$ ;

    end

   end

( $\mathit{8}$ ) else for $k=t+\mathit{1},\dots ,t+\mathrm{\Delta }-\mathit{1}$

( $\mathit{9}$ )  for $n=\mathit{1},\dots ,N_k$

( $\mathit{10}$ )   if $X_k^n$ is a leaf node

( $\mathit{11}$ )    $X=X_k^n$ ;

      //add auxiliary particles

( $\mathit{12}$ )   for $j=k+\mathit{1},\dots ,t+\mathrm{\Delta }$

( $\mathit{13}$ )     $u_j^{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}}=E_{p_f}\left(\cdot \mid X\right)$ ;

( $\mathit{14}$ )     $X=u_j^{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}}$ ;

( $\mathit{15}$ )     $\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}=\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}+\mathit{1}$ ;

     end

     end

    end

   end

( $\mathit{16}$ ) return $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_t^i\right)$

Lemma 8.

The number of leaf nodes in $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_t^i\right)$ is equal to the number of leaf nodes in $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_t^i\right)$ .

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.

Let $N_k$ denote the number of children/grandchildren of $X_t^i$ in $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_t^i\right)$ at time k, where $t+\mathrm{1}\le k\le t+\Delta$ . The time complexity $T_{CAGF}\left(X_t^i\right)$ of Procedure 2 satisfies

$$\begin{array}{cc}& \mathcal{O}\left({\displaystyle \sum }_{k=t+\mathrm{1}}^{t+\Delta }{N}_k\right)\le T_{CAGF}\left(X_t^i\right)\le \mathcal{O}\left(\left(\Delta +\mathrm{1}\right){\displaystyle \sum }_{k=t+\mathrm{1}}^{t+\Delta }{N}_k\right).\end{array}$$

(18)

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 complexity $T_{ADWPF}$ of Algorithm 3 is $O\left(CN\right)$ , $C\ll N$ .

Proof.

Obviously, the time cost of Algorithm 3 mainly consists of steps (4) and (6). Let $T_{\mathrm{C}\mathrm{A}\mathrm{G}\mathrm{F}}\left(X_{t-\Delta }^i\right)$ , $\mathrm{1}\le i\le N$ , denote the time cost of step (4) in ADWPF, and $N_k^i$ denotes the number of children/grandchildren of $X_{t-\Delta }^i$ in $\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_{t-\Delta }^i\right)$ at time k, where $t-\Delta +\mathrm{1}\le k\le t$ . According to Lemma 9

$$\begin{array}{cc}& {\displaystyle \sum }_{i=\mathrm{1}}^N{T}_{\mathrm{C}\mathrm{A}\mathrm{G}\mathrm{F}}\left(X_{t-\Delta }^i\right)\le {\displaystyle \sum }_{i=\mathrm{1}}^N\mathcal{O}\left(\left(\Delta +\mathrm{1}\right){\displaystyle \sum }_{k=t-\Delta +\mathrm{1}}^t{N}_k^i\right)=\mathcal{O}\left(\left(\Delta +\mathrm{1}\right){\displaystyle \sum }_{i=\mathrm{1}}^N\hspace{0.17em}{\displaystyle \sum }_{k=t-\Delta +\mathrm{1}}^t{N}_k^i\right)=\mathcal{O}\left(\left(\Delta +\mathrm{1}\right)N\right),\end{array}$$

(19)

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

$$\begin{array}{cc}& {\displaystyle \sum }_{i=\mathrm{1}}^N{T}_{\mathrm{C}\mathrm{A}\mathrm{G}\mathrm{F}}\left(X_{t-\Delta }^i\right)\ge {\displaystyle \sum }_{i=\mathrm{1}}^N\mathcal{O}\left({\displaystyle \sum }_{k=t-\Delta +\mathrm{1}}^t{N}_k^i\right)=\mathcal{O}\left({\displaystyle \sum }_{i=1}^N\hspace{0.17em}{\displaystyle \sum }_{k=t-\Delta +1}^t{N}_k^i\right)=\mathcal{O}\left(N\right).\end{array}$$

(20)

Let $T_{\mathrm{B}\mathrm{S}}\left(X_{t-\Delta }^i\right)$ , $\mathrm{1}\le i\le N$ , denote the time cost of step (6) in ADWPF. The sum number of leaf nodes of all generation trees $\left\{\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_{t-\Delta }^1\right),\dots ,\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_{t-\Delta }^N\right)\right\}$ is less than $\Delta N$ and greater than N, so the sum number of leaf nodes of all auxiliary generation trees $\left\{\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_{t-\Delta }^1\right),\dots ,\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_{t-\Delta }^N\right)\right\}$ is less than $\Delta N$ and greater than N according to Lemma 8. Since the time cost of Procedure 1 is $O\left(\Delta \right)$ for each leaf node, we have

$$\begin{array}{cc}& N\ast O\left(\Delta \right)\le {\displaystyle \sum }_{i=\mathrm{1}}^N{T}_{\mathrm{B}\mathrm{S}}\left(X_{t-\Delta }^i\right)\le \Delta N\ast O\left(\Delta \right)=O\left({\Delta }^{\mathrm{2}}N\right).\end{array}$$

(21)

So for the time complexity $T_{\mathrm{A}\mathrm{D}\mathrm{W}\mathrm{P}\mathrm{F}}$ of Algorithm 2 we have

$$\begin{array}{cc}& T_{\mathrm{A}\mathrm{D}\mathrm{W}\mathrm{P}\mathrm{F}}=T\ast \left({\displaystyle \sum }_{i=\mathrm{1}}^N{T}_{\mathrm{C}\mathrm{A}\mathrm{G}\mathrm{F}}\left(X_{t-\Delta }^i\right)+{\displaystyle \sum }_{i=\mathrm{1}}^N{T}_{\mathrm{B}\mathrm{S}}\left(X_{t-\Delta }^i\right)\right)\le T\ast \left(\mathcal{O}\left(\left(\Delta +\mathrm{1}\right)N\right)+O\left({\Delta }^{\mathrm{2}}N\right)\right)=\mathcal{O}\left(\left({\Delta }^{\mathrm{2}}+\Delta +\mathrm{1}\right)TN\right).\end{array}$$

(22)

Meanwhile,

$$\begin{array}{cc}& T_{\mathrm{A}\mathrm{D}\mathrm{W}\mathrm{P}\mathrm{F}}=T\ast \left({\displaystyle \sum }_{i=\mathrm{1}}^N{T}_{\mathrm{C}\mathrm{A}\mathrm{G}\mathrm{F}}\left(X_{t-\Delta }^i\right)+{\displaystyle \sum }_{i=\mathrm{1}}^N{T}_{\mathrm{B}\mathrm{S}}\left(X_{t-\Delta }^i\right)\right)\ge T\ast \left(\mathcal{O}\left(N\right)+O\left(\Delta N\right)\right)=\mathcal{O}\left(\left(\Delta +\mathrm{1}\right)TN\right),\end{array}$$

(23)

where $\Delta ,T\ll N$ . So the time complexity of ADWPF is $O\left(CN\right)$ , $C\ll N$ .

Algorithm 2: DWFB.

Assume N is the number of particles, t is current time index, count is the number of simulation steps, $\mathrm{\Delta }$ is the delayed time interval,

$w_t^i$ is the filtered weight and ${\widehat{w}}_t^i$ is the modified weight.

( $\mathit{1}$ )  for $t=\mathit{1},\dots ,T$

( $\mathit{2}$ )   Generate samples ${\left\{X_t^i,w_t^i,{\widehat{w}}_t^i=\mathit{0}\right\}}_{i=\mathit{1}}^N$ through a forward particle filter (like Algorithm 1).

( $\mathit{3}$ )  if $t>\mathrm{\Delta }$

( $\mathit{4}$ )   for $i=\mathit{1}:N$

( $\mathit{5}$ )     ${\widehat{w}}_t^i=w_t^i$

( $\mathit{6}$ )    for $k=T-\mathit{1}:T-∆$

( $\mathit{7}$ )     for $i=\mathit{1}:N$

             $${\widehat{w}}_k^i={\displaystyle \sum }_{j=\mathit{1}}^N{\widehat{w}}_{k+\mathrm{1}}^j\frac{w_k^i{p}_f\left(X_{k+\mathit{1}}^j\mid X_k^i\right)}{\left[{\sum }_l^N{w}_k^l{p}_f\left(X_{k+\mathit{1}}^j\mid X_k^l\right)\right]}$$

      end

     end

    end

   end

 The inference result of system state $x_{T-\mathrm{\Delta }}$ is

             $${\widehat{x}}_{t-\mathrm{\Delta }}={\displaystyle \sum }_{i=\mathit{1}}^N{\widehat{w}}_{t-\mathrm{\Delta }}^i\cdot X_{t-\mathrm{\Delta }}^i$$

Algorithm 3: ADWPF.

Let N is the number of particles, t is current time index, $\mathrm{\Delta }$ is the delayed time length, $\mathrm{\hspace{0.17em}\hspace{0.17em}}{w}_t^i$ is the forward filtered weight

and ${\widehat{w}}_t^i$ is the delayed weight.

$(\mathit{1}$ ) for $t=\mathit{1},\dots ,T$

   Generate samples ${\left\{X_t^i,w_t^i,{\widehat{w}}_t^i=\mathit{0}\right\}}_{i=\mathit{1}}^N$ through a forward particle filter (like Algorithm 1).

   end

( $\mathit{2}$ ) if $t>\mathrm{\Delta }$

( $\mathit{3}$ )  for $i=\mathit{1},\dots ,N$

    //Procedure 2.

( $\mathit{4}$ )    $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_{t-\mathrm{\Delta }}^i\right)=\mathrm{C}\mathrm{A}\mathrm{G}\mathrm{F}\left(\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_{t-\mathrm{\Delta }}^i\right),\mathrm{\Delta }\right)$

( $\mathit{5}$ )   for each leaf node $X_t^j$ of $\mathrm{a}\mathrm{u}\mathrm{x}\_\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\left(X_{t-\mathrm{\Delta }}^i\right)$

     //Procedure 1.

( $\mathit{6}$ )     $P\left(X_t^j\to X_{t-\mathrm{\Delta }}^i\right)=\mathrm{B}\mathrm{S}\left(X_t^j\right)$ ;

( $\mathit{7}$ )     ${\widehat{w}}_{t-\mathrm{\Delta }}^i={\widehat{w}}_{t-\mathrm{\Delta }}^i+P\left(X_t^j\to X_{t-\mathrm{\Delta }}^i\right)$ ;

( $\mathit{8}$ )     $X=X_k^n$ ;

    end

   end

   end

5. Simulation and Results

To verify the performance of our algorithm ADWPF, we consider a wireless sensor network with 70 anchor sensors ${\left\{a_i={\left(a_x^i,a_y^i\right)}^T\right\}}_{i=1}^{70}$ . These sensors are randomly distributed in a $\mathrm{50}\times \mathrm{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 $x_t={\left({\mathrm{l}\mathrm{o}\mathrm{c}}_{x,t},v_{x,t},{\mathrm{l}\mathrm{o}\mathrm{c}}_{y,t},v_{y,t}\right)}^T$ , where ${\left({\mathrm{l}\mathrm{o}\mathrm{c}}_{x,t},{\mathrm{l}\mathrm{o}\mathrm{c}}_{y,t}\right)}^T$ is location of the mobile senor and ${\left(v_{x,t},v_{y,t}\right)}^T$ is its velocity. In our simulation, we set ${\left(v_{x,t},v_{y,t}\right)}^T={\left(1,1\right)}^T$ . So we just need to utilize our algorithm to infer the state ${\mathrm{l}\mathrm{o}\mathrm{c}}_t={\left({\mathrm{l}\mathrm{o}\mathrm{c}}_{x,t},{\mathrm{l}\mathrm{o}\mathrm{c}}_{y,t}\right)}^T$ of the mobile sensor at time t. The communication radius of sensor is set as 5. At time t, the distance measurement $y_t$ from anchor sensor $a_i$ is

$$\begin{array}{cc}& y_t=\sqrt{{\left(a_x^i-{\mathrm{l}\mathrm{o}\mathrm{c}}_{x,t}\right)}^{\mathrm{2}}+{\left(a_y^i-{\mathrm{l}\mathrm{o}\mathrm{c}}_{y,t}\right)}^{\mathrm{2}}}+n,\end{array}$$

(24)

where the distance measurement noise $n\in R$ follows $n~\mathcal{N}\left(\mathrm{0,1}\right)$ .

5.1. Linear Motion Model

The state of the mobile sensor evolves according to [16, 17]:

$$\begin{array}{cc}& x_t=\Phi x_{t-\mathrm{1}}+\Gamma w_t,\end{array}$$

(25)

where $w_t\in R^{\mathrm{2}}$ , $w_t~\mathcal{N}\left(\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\end{array}\right],\left[\begin{array}{cc}\mathrm{0.0001}& \mathrm{0}\\ \mathrm{0}& \mathrm{0.0001}\end{array}\right]\right)$ , and Φ and Γ have the following form:

$$\begin{array}{cc}& \Phi =\left[\begin{array}{cccc}1& 1& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 1\\ 0& 0& 0& 1\end{array}\right],\\ \Gamma =\left[\begin{array}{cc}0.5& 0\\ 1& 0\\ 0& 0.5\\ 0& 1\end{array}\right].\end{array}$$

(26)

The initial state of the mobile sensor is $x_0={\left(\mathrm{0,1},0,1\right)}^T$ , the sampling interval is $\gamma =\mathrm{1}$ , the number of tracking steps is $T=\mathrm{50}$ , the number of simulations is $\mathrm{n}\mathrm{u}\mathrm{m}=\mathrm{50}$ , and we use the mean squared error (MSE) as the performance measurement. It is defined as follows:

$$\begin{array}{cc}& \mathrm{M}\mathrm{S}\mathrm{E}=\frac{\mathrm{1}}{T\ast \mathrm{n}\mathrm{u}\mathrm{m}}{\displaystyle \sum }_{j=\mathrm{1}}^{\mathrm{n}\mathrm{u}\mathrm{m}}{\displaystyle \sum }_{t=\mathrm{1}}^T{\left({\mathrm{l}\mathrm{o}\mathrm{c}}_t^j-{\mathrm{e}\mathrm{s}\mathrm{t}\_\mathrm{l}\mathrm{o}\mathrm{c}}_t^j\right)}^{\mathrm{2}},\end{array}$$

(27)

where ${\mathrm{l}\mathrm{o}\mathrm{c}}_t^j$ is the true location of the mobile sensor at time t in the jth simulation and ${\mathrm{e}\mathrm{s}\mathrm{t}\_\mathrm{l}\mathrm{o}\mathrm{c}}_t^j$ 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.

Table 1

Average MSE for bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in linear motion model.

MSE	$$N=\mathrm{50}$$	$$N=\mathrm{100}$$	$$N=\mathrm{200}$$	$$N=\mathrm{500}$$	$$N=\mathrm{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 $\left(\mathrm{\Delta }=\mathrm{1}\right)$	0.123	0.070	0.055	0.042	0.039
DWFB $\left(\mathrm{\Delta }=\mathrm{3}\right)$	0.107	0.065	0.044	0.030	0.029
DWFB $\left(\mathrm{\Delta }=\mathrm{5}\right)$	0.106	0.058	0.034	0.025	0.024
ADWPF   $\left(\mathrm{\Delta }=\mathrm{1}\right)$	0.117	0.063	0.047	0.035	0.032
ADWPF   $\left(\mathrm{\Delta }=\mathrm{3}\right)$	0.097	0.053	0.034	0.019	0.017
ADWPF   $\left(\mathrm{\Delta }=\mathrm{5}\right)$	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 $v_t$ and direction $d_t$ of the mobile target evolves as follows:

$$\begin{array}{cc}& v_t=\alpha v_{t-\mathrm{1}}+\left(\mathrm{1}-\alpha \right)\stackrel{-}{v}+\sqrt{\left(\mathrm{1}-{\alpha }^{\mathrm{2}}\right)}{n}_v,\\ d_t=\alpha d_{t-1}+\left(1-\alpha \right)\stackrel{-}{d}+\sqrt{\left(1-{\alpha }^2\right)}{n}_d,\end{array}$$

(28)

where $n_v$ and $n_d$ are white noise, α is the tuning parameter used to control the randomness, and $\stackrel{-}{v}$ and $\stackrel{-}{d}$ are constants representing the mean value of speed and direction. The location ${\mathrm{l}\mathrm{o}\mathrm{c}}_t={\left({\mathrm{l}\mathrm{o}\mathrm{c}}_{x,t},{\mathrm{l}\mathrm{o}\mathrm{c}}_{y,t}\right)}^T$ of the mobile target evolves as follows:

$$\begin{array}{cc}& {\mathrm{l}\mathrm{o}\mathrm{c}}_{x,t}={\mathrm{l}\mathrm{o}\mathrm{c}}_{x,t-\mathrm{1}}+v_t\cos d_t,\\ {\mathrm{l}\mathrm{o}\mathrm{c}}_{y,t}={\mathrm{l}\mathrm{o}\mathrm{c}}_{y,t-1}+v_t\sin d_t.\end{array}$$

(29)

We set $\alpha =0.75$ , $\stackrel{-}{v}=\mathrm{1}$ , $\stackrel{-}{d}={90}^{°}$ , ${\mathrm{l}\mathrm{o}\mathrm{c}}_{x,0}=25$ , ${\mathrm{l}\mathrm{o}\mathrm{c}}_{y,0}=25$ , and $v_0=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=\mathrm{50}$$	$$N=\mathrm{100}$$	$$N=\mathrm{200}$$	$$N=\mathrm{500}$$	$$N=\mathrm{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 $\left(\mathrm{\Delta }=\mathrm{1}\right)$	2.577	2.254	2.159	2.103	2.101
DWFB $\left(\mathrm{\Delta }=\mathrm{3}\right)$	2.598	2.399	2.399	2.312	2.344
DWFB $\left(\mathrm{\Delta }=\mathrm{5}\right)$	2.029	1.911	1.781	1.763	1.766
ADWPF   $\left(\mathrm{\Delta }=\mathrm{1}\right)$	1.880	1.546	1.441	1.370	1.366
ADWPF   $\left(\mathrm{\Delta }=\mathrm{3}\right)$	1.204	0.999	0.880	0.795	0.781
ADWPF   $\left(\mathrm{\Delta }=\mathrm{5}\right)$	1.178	0.893	0.748	0.606	0.603

As illustrated in Tables 1 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=\mathrm{50}$ ) and Gaussian-Markov motion model ( $T=\mathrm{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=\mathrm{50}$ ).

MSE	$$N=\mathrm{50}$$	$$N=\mathrm{100}$$	$$N=\mathrm{200}$$	$$N=\mathrm{500}$$	$$N=\mathrm{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 $\left(\mathrm{\Delta }=\mathrm{1}\right)$	0.246	0.927	3.713	22.322	56.934
DWFB $\left(\mathrm{\Delta }=\mathrm{3}\right)$	0.663	2.592	10.494	63.828	164.757
DWFB $\left(\mathrm{\Delta }=\mathrm{5}\right)$	1.040	4.096	16.320	101.648	259.860
ADWPF   $\left(\mathrm{\Delta }=\mathrm{1}\right)$	0.042	0.066	0.124	0.312	0.512
ADWPF   $\left(\mathrm{\Delta }=\mathrm{3}\right)$	0.084	0.153	0.317	0.803	1.300
ADWPF   $\left(\mathrm{\Delta }=\mathrm{5}\right)$	0.154	0.308	0.610	1.520	2.458

Table 4

Total time (sec.) cost of bootstrap particle filter, auxiliary particle filter, DWFB, and our method ADWPF in Gaussian-Markov motion model ( $T=\mathrm{50}$ ).

MSE	$$N=\mathrm{50}$$	$$N=\mathrm{100}$$	$$N=\mathrm{200}$$	$$N=\mathrm{500}$$	$$N=\mathrm{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 $\left(\mathrm{\Delta }=\mathrm{1}\right)$	0.271	0.925	3.654	22.719	58.071
DWFB $\left(\mathrm{\Delta }=\mathrm{3}\right)$	0.679	2.632	10.448	65.003	166.134
DWFB $\left(\mathrm{\Delta }=\mathrm{5}\right)$	1.051	4.167	16.664	103.563	263.259
ADWPF   $\left(\mathrm{\Delta }=\mathrm{1}\right)$	0.044	0.055	0.109	0.266	0.444
ADWPF   $\left(\mathrm{\Delta }=\mathrm{3}\right)$	0.087	0.147	0.287	0.709	1.135
ADWPF   $\left(\mathrm{\Delta }=\mathrm{5}\right)$	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 time t to reweight Monte Carlo samples ${\left\{X_t^i\right\}}_{i=1}^N$ 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].