Distributed Joint Cooperative Self-Localization and Target Tracking Algorithm for Mobile Networks

Abstract

Location information is a key issue for applications of the Internet of Things. In this paper, we focus on mobile wireless networks with moving agents and targets. The positioning process is divided into two phases based on the factor graph, i.e., a prediction phase and a joint self-location and tracking phase. In the prediction phase, we develop an adaptive prediction model by exploiting the correlation of trajectories within a short period to formulate the prediction message. In the joint positioning phase, agents calculate the cooperative messages according to variational message passing and locate themselves. Simultaneously, the average consensus algorithm is employed to realize distributed target tracking. The simulation results show that the proposed prediction model is adaptive to the random movement of nodes. The performance of the proposed joint self-location and tracking algorithm is better than the separate cooperative self-localization and tracking algorithms.

1. Introduction

With the development of the Internet of Things (IoT) and related new information technology, wireless sensor networks (WSNs) are expected to gradually penetrate various industries and applications [1,2,3,4], influencing all aspects of people’s lives. WSNs connect the material world with the human world. Embedded in various “things”, sensors sense the states of the “things” and transmit collected information to processing terminals via the Internet to achieve monitoring and management in real-time [5]. The information collected by sensors is meaningful only when combined with location information, which helps the administrators know what, when, and where “things” happened. Therefore, cooperative self-localization and target tracking are two key issues in location-based applications of WSNs.

In WSNs, there are a small number of anchors and thousands of agents. Anchors are equipped with Global Position System (GPS) receivers or other devices with known coordinates, while the agents use a cooperative self-localization algorithm to estimate their own positions. In addition, when a target enters the WSNs, anchors and agents work cooperatively to track and locate it. In static WSNs, anchors and agents are fixed and so their coordinates do not change with time, meaning agents only need to locate themselves once. In mobile WSNs, anchors and agents move with time. Therefore, agents need to locate themselves in a sequential manner. In contrast with static WSNs, mobile WSNs alleviate several issues, such as coverage optimization and target tracking. A plethora of localization systems are proposed in the literature. In [6], the authors provided a general overview of localization in WSNs and showed the importance of different localization approaches in modern IoT applications. In [7], the authors focused on mobile WSNs and analyzed localization algorithms in mobile WSNs.

Cooperative self-localization and target tracking are closely related because target tracking can only be carried out when the locations of agents are clear. In recent years, some researches combined self-localization with target tracking and proposed a series of joint location algorithms. In [8], a simultaneous localization and tracking algorithm was introduced to solve the problem of tracking a non-cooperative target while simultaneously localizing and calibrating the static nodes in the network. In [9], a Bayesian method based on belief propagation (BP) was proposed for the distributed sequential localization of mobile networks composed of both cooperative agents and non-cooperative targets. This provided a consistent combination of cooperative self-localization and distributed tracking. However, the motion parameters (such as velocity) and the state-transition probability density function (pdf) of each node were assumed to be available. In [10], a distributed variational filtering was presented to simultaneously localize the detecting sensors and track the target by exploiting a series of measurements generated in the sensors when the target moved through the network field.

Target tracking is divided into two categories, i.e., centralized and distributed. In centralized target tracking, the nodes that detect the target transmit their measurements related to the target to a central processing node in a hop-by-hop communication mode. Then, the central node determines the positioning of the target by uniformly processing all of the measurements. However, the communication overhead of multi-hop data transmission is very large and energy consumption is too high. In distributed target tracking, each node locates the target according to all measurements related to the target. Therefore, the key problem is how to make each node obtain all of the measurements. In [11], a scheme was proposed in which agents broadcasted their own measurements, identification numbers, and additional information, such as the ID of the target. They also received these from other nodes and updated the information. All of the agents repeated the process until each agent obtained all of the measurement information related to the target. In this way, each agent was a fusion center and could perform tracing tasks. However, much redundant information existed in the network, and the required storage space was also very large. Therefore, this method is suitable for a fully connected network with a large communication radius. Consensus algorithm [12,13,14] provides another way for each agent to communicate observations with the remaining agents engaged in the tracking. Consensus is a distributed iterative algorithm relying on communication links between neighboring nodes and can achieve entire network polymerization (such as sum, average, or maximum). This does not need a routing algorithm and is robust to changes in network topology and unreliable network environments. Therefore, consensus algorithm is widely used in distributed target tracking [15,16]. In [17], each node calculated the mean and variance of the local posterior by weighted sampling. Then, the average consensus algorithm was employed to calculate the approximation of the global posterior to complete distributed target tracking. In [18], the global likelihood function was approximated to a Gaussian distribution, and each node obtained the approximation through consensus iteration. Furthermore, a consensus algorithm, which is applicable to any exponential likelihood function, was proposed in [19].

As an implementation of approximate Bayesian inference, messaging passing algorithm is particularly well-suited to distributed calculation and has been applied to cooperative self-localization and tracking of WSNs [9,20,21]. In previous work [22], we proposed a distributed cooperative self-localization algorithm by employing variational message passing (VMP) on factor graphs for static WSNs. With regard to the non-Gaussian messages caused by the nonlinear ranging model, we approximated them to Gaussian messages by exploiting second-order Taylor expansion, which significantly reduced the computational complexity. Then, we further combined the VMP approach in [22] with BP and proposed a cooperative self-localization algorithm for mobile WSNs in [23]. In this paper, we focus on mobile networks with anchors, moving agents, and targets deployed in a two-dimensional plane. In particular, the anchors have perfect location information at all times, while the cooperative agents and the non-cooperative targets are at unknown positions and can move independently. We divide the positioning process into two phases, i.e., a prediction phase and a joint self-location and tracking phase. Firstly, an adaptive prediction model is proposed to calculate the prediction message from the previous time slot. Then, in the joint self-location and tracking phase, the agents locate themselves using the measurements from neighboring nodes (agents and anchors) and the detectable targets, and the average consensus algorithm is employed to achieve distributed positioning of the targets. The simulation results show that the proposed prediction model is more adaptive and accurate than instant prediction, linear prediction, or square prediction, and the proposed joint self-location and tracking algorithm is better than the separate cooperative self-localization and target tracking algorithms.

This paper is organized as follows. In Section 2, we first describe the system model of cooperative self-localization, and tracking. Then, the corresponding probabilistic model and factor graph are analyzed. In Section 3, the proposed distributed joint self-location and tracking algorithm are described in detail, including a prediction phase and a joint self-location and tracking phase. The performance of the proposed prediction model and joint localization algorithm are evaluated in Section 4. Finally, conclusions are drawn in Section 5.

Notations: Boldface lowercase and uppercase letters denote vectors and matrices, respectively. Superscript ( )^T and symbol $||·||$ stand for transposition and Euclidian norm, respectively. The probability density function (pdf) of a 1D Gaussian distribution with mean $\mu$ and variance ${\sigma }^2$ is represented by $\mathcal{N}(x;\mu ,{\sigma }^2)$ while the pdf of a 2D Gaussian distribution with mean vector $\mathit{\mu }$ and covariance matrix $\mathit{V}$ is represented by $\mathcal{N}(\mathit{x};\mathit{\mu },\mathit{V})$.

2. System Model and Factor Graph

We consider a network with anchors, moving agents, and targets deployed in a two-dimensional plane, which can be defined by a vertices set $\mathcal{V}$ and an edges set $\mathcal{E}$, i.e., $\mathcal{G}=(\mathcal{V},\mathcal{E})$. Here, each vertex $i\in \mathcal{V}$ represents a node in the plane and each edge $(i,j)\in \mathcal{E}$ indicates a neighborhood structure between node i and node j. We divide set $\mathcal{V}$ into three parts: $\mathcal{V}=\mathcal{A}\cup \mathcal{M}\cup \mathcal{S}$, where $\mathcal{A}$ designates all anchors with perfect location information, $\mathcal{M}$ includes all agents at unknown positions, and $\mathcal{S}$ represents all non-cooperative targets. Moreover, the number of anchors, agents, and targets are denoted by $N_A$, $N_M$, and $N_O$, respectively.

Time is slotted and we assume that all agents and targets move independently. At the kth time slot, the position variable of node $i\in \mathcal{V}$ is denoted by ${\mathit{x}}_i^k\triangleq {\left[x_i^k,y_i^k\right]}^{\mathrm{T}}$. For simplicity of illustration, the communication link between cooperative neighbors is assumed to be symmetric, which means $(i,j)\in \mathcal{E}$ and $(j,i)\in \mathcal{E}$ if $|{\mathit{x}}_i^k-{\mathit{x}}_j^k|\le R$, where R is the communication radius. At the kth time slot, if $(i,j)\in \mathcal{E}$, nodes i can obtain the current position of node j and a ranging measurement $d_{i\leftarrow j}^k$ with Gaussian noise ${\omega }_{ij}^k\sim \mathcal{N}({\omega }_{ij}^k;0,{\left({\sigma }_{ij}^k\right)}^2)$. Consequently, the ranging measurement $d_{i\leftarrow j}^k$ can be formulated as

$$d_{i\leftarrow j}^k=||x_i^k-x_j^k||+{\omega }_{ij}^k,$$

and the likelihood function of the position variables of $x_i^k$ and $x_j^k$ is a Gaussian pdf, i.e.,

$$p(d_{i\leftarrow j}^k\left|{\mathit{x}}_i^k,{\mathit{x}}_j^k)\right.\sim \mathcal{N}(d_{i\leftarrow j}^k;\left|\right|{\mathit{x}}_i^k-{\mathit{x}}_j^k\left|\right|,{\sigma }_{ij}^k).$$

For convenience of description, we define ${\mathcal{V}}_o^k\triangleq {\mathcal{A}}_o^k\cup {\mathcal{M}}_o^k$ (${\mathcal{A}}_o^k$ for anchors and ${\mathcal{M}}_o^k$ for agents) to represent all nodes who can detect target at the kth time slot. Similarly, let ${\mathcal{V}}_m^k\triangleq {\mathcal{A}}_m^k\cup {\mathcal{M}}_m^k\cup {\mathcal{S}}_m^k$ denote all nodes from whom node $m\in \mathcal{M}$ can implement ranging measurement at the kth time slot. Typically, notation ${\mathcal{S}}_m^k$ represents all targets that node m can detect at the kth time slot. Moreover, we denote the position variables of all nodes by ${\mathit{X}}^k\triangleq \{{\mathit{x}}_i^k,\forall i\in \mathcal{A}\cup \mathcal{M}\cup \mathcal{S}\}$, the ranging measurements of node $m\in \mathcal{M}$ by ${\mathit{Z}}_m^k\triangleq \{d_{m\leftarrow j}^k,\forall j\in {\mathcal{A}}_m^k\cup {\mathcal{M}}_m^k\cup {\mathcal{S}}_m^k\}$, and the ranging measurements of target $O\in \mathcal{S}$ by ${\mathit{Z}}_o^k\triangleq \{d_{j\leftarrow \mathrm{o}}^k,\forall j\in {\mathcal{A}}_o^k\cup {\mathcal{M}}_o^k\}$, respectively. Then, let ${\mathit{X}}^k\triangleq \{{\mathit{x}}_i^k,\forall i\in \mathcal{V}\}$ and ${\mathit{Z}}^k\triangleq \{{\mathit{Z}}_i^k,\forall i\in \mathcal{M}\cup \mathcal{S}\}$ represent the position variables of all nodes and ranging measurements at the kth time slot, respectively. Furthermore, let ${\mathit{X}}^{0:K}\triangleq \{{\mathit{X}}^k,k=0:K\}$ denote the position variables of all nodes from the 0th time slot to the kth time slot and ${\mathit{Z}}^{1:K}\triangleq \{{\mathit{Z}}^k,k=1:K\}$ denote the ranging measurements from the 1st time slot to the Kth time slot, respectively.

According to Bayesian rules and the assumption that all nodes move independently, the joint a posteriori distribution of ${\mathit{X}}^{0:K}$ with given observations ${\mathit{Z}}^{0:K}$ can be formulated as

$$\begin{array}{cc}\hfill p({\mathit{X}}^{0:K}\left|{\mathit{Z}}^{1:K}\right.)& \propto p({\mathit{X}}^0)\prod _{k=1}^{K}p({\mathit{Z}}^k\left|{\mathit{X}}^k\right.)p({\mathit{X}}^k\left|{\mathit{X}}^{k-1}\right.)\hfill \\ & \propto \prod _{i\in \mathcal{V}}p({\mathit{x}}_i^0)\prod _{k=1}^K\prod _{m\in \mathcal{M}}\prod _{j\in {\mathcal{V}}_m^k}p(d_{m\leftarrow j}^k\left|{\mathit{x}}_m^k,{\mathit{x}}_j^k)\right.\prod _{o\in \mathcal{S}}\prod _{j\in {\mathcal{V}}_o^k}p(d_{j\leftarrow o}^k\left|{\mathit{x}}_o^k,{\mathit{x}}_j^k)\right.\prod _{i\in \mathcal{V}}p({\mathit{x}}_i^k|{\mathit{x}}_i^{k-1})\hfill \end{array},$$

(1)

where $p\left({\mathit{x}}_i^0\right)$ is the a prior distribution of node I, which is assumed to be a Gaussian pdf with mean ${\mu }_i^0$ and covariance matrix ${\mathit{V}}_i^0\triangleq {\sigma }_{i0}^2{\mathit{I}}_{2\times 2}$, $p(d_{i\leftarrow j}^k\left|{\mathit{x}}_i^k,{\mathit{x}}_j^k)\right.$ is the likelihood function of node i and node j, and $p\left({\mathit{x}}_i^k\right|{\mathit{x}}_i^{k-1})$ is the probabilistic state-transition function.

For simplicity, we define $f_i^{k|k-1}\triangleq p\left({\mathit{x}}_i^k\right|{\mathit{x}}_i^{k-1})$ and $f_{ij}^k\triangleq p(d_{i\leftarrow j}^k\left|{\mathit{x}}_i^k,{\mathit{x}}_j^k)\right.$. Based on the factorization in (1), the joint a posteriori distribution of $p\left({\mathit{X}}^{0:K}\left|{\mathit{Z}}^{0:K}\right.\right)$ can be represented by a factor graph, as shown in Figure 1. In the factor graph, each variable node represents the position variable ${\mathit{x}}_i^k$ of node i and is depicted by a circle, while each factor node represents a local function $f_i^{k|k-1}$ or $f_{ij}^k$ and is drawn by a square. Moreover, if a variable is an argument of a local function, it is connected to the factor node through an edge.

Generally speaking, factor graphs are undirected digraphs, which means the edge between a factor node and a variable node is undirected. However, in Figure 1, messages only flow forward in time regarding the spatiotemporal constraints of the network. We do not calculate the messages from the present to the past because network connectivity may change and the states of the nodes may be outdated. Furthermore, the variable node of an anchor’s position ignores the messages from its neighbors as anchors’ positions are known.

3. Distributed Joint Self-Location and Tracking Algorithm

As shown in Figure 1, for the cooperative agent $m\in \mathcal{M}$ and the non-cooperative target $O\in \mathcal{S}$, the beliefs of ${\mathit{x}}_m^k$ and ${\mathit{x}}_O^k$, denoted by $b\left({\mathit{x}}_m^k\right)$ and $b\left({\mathit{x}}_o^k\right)$, consist of two kinds of message: Prediction messages from the (k−1)th time slot and cooperative messages from the neighbors. Based on message passing rules, for node $i\in \mathcal{M}\cup \mathcal{S}$, the belief $b\left({\mathit{x}}_i^k\right)$ can be described as

$$b\left({\mathit{x}}_i^k\right)=\frac{1}{Z}{m}_{f_i^{k|k-1}\to {\mathit{x}}_i^k}\left({\mathit{x}}_i^k\right)\prod _{j\in {\mathcal{V}}_i^k}{m}_{f_{ij}^k\to {\mathit{x}}_i^k}\left({\mathit{x}}_i^k\right),$$

(2)

where Z is the normalization constant and $m_{f_i^{k|k-1}\to {\mathit{x}}_i^k}\left({\mathit{x}}_i^k\right)$ and $m_{f_{ij}^k\to {\mathit{x}}_i^k}\left({\mathit{x}}_i^k\right)$ represent messages from factor node and factor nodes $f_{ij}^k$ ($\forall j\in {\mathcal{V}}_i^k$) to variable node ${\mathit{x}}_i^k$, respectively. Therefore, agent $m\in \mathcal{M}$ determines $b\left({\mathit{x}}_m^k\right)$ and $b\left({\mathit{x}}_o^k\right)$($m\in {\mathcal{V}}_o^k$) by two phases, i.e., a prediction phase and a joint self-location and tracking phase.

3.1. Adaptive Prediction Model and Prediction Message Calculating

In (1), the probabilistic state-transition function $p(x_i^k|x_i^{k-1})$ of node i is related to the position prediction model. In this paper, according to the inertia of the motion, we propose an adaptive prediction model by exploiting the correlation of trajectories of node i within a short period.

We denote the trajectory of node i from the (k−N)th time slot to the (k−1)th time slot as $T_i(k-1,N)\triangleq {\left({({\widehat{x}}_i^{k-1})}^{\mathrm{T}},{({\widehat{x}}_i^{k-2})}^{\mathrm{T}},\cdots ,{({\widehat{x}}_i^{k-N})}^{\mathrm{T}}\right)}^{\mathrm{T}}$, where ${\widehat{x}}_i^{k-N}\triangleq {[{\widehat{x}}_i^{k-N},{\widehat{y}}_i^{k-N}]}^{\mathrm{T}}$, …,${\widehat{x}}_i^{k-1}\triangleq {[{\widehat{x}}_i^{k-1},{\widehat{y}}_i^{k-1}]}^{\mathrm{T}}$ are the estimated positions of node i at the (k−N)th, …, (k−1)th time slot. Similarly, the trajectory of node i from the (k−N + 1)th time slot to the kth time slot can be described as $T_i(k,N)\triangleq {\left({({\widehat{x}}_i^k)}^{\mathrm{T}},{({\widehat{x}}_i^{k-1})}^{\mathrm{T}},\cdots ,{({\widehat{x}}_i^{k-N+1})}^{\mathrm{T}}\right)}^{\mathrm{T}}$. Then, we formulate the prediction model of node i as

$$T_i(k,N)=P^i·T_i(k-1,N) ,$$

(3)

where $P^i$ is the prediction matrix of node i and its dimension is 2N × 2N. Obviously, there are only two different elements between $T_i(k-1,N)$ and $T_i(k,N)$, i.e., ${\widehat{x}}_i^k$ and ${\widehat{x}}_i^{k-N}$. Therefore, the elements from the 3rd to the 2Nth rows are

$$p_{r,s}^i=\left\{\begin{array}{ll}0,& r\ge 3 \mathrm{and} r\ne s+2\\ 1,& r\ge 3 \mathrm{and} r=s+2\end{array}\right..$$

(4)

According to the inertia of the motion of node i, we assume that the prediction matrix $P^i$ is unchanged in a short period. Based on the prediction model in (4), we have

$$T_i(k-1,N)=P^i·T_i(k-2,N) .$$

(5)

Consequently, it can be formulated as

$$T_i(k,N)={(P^i)}^n·T_i(k-n,N).$$

(6)

Based on Equations (3) and (5), we have

$$T_i(k,N)-T_i(k-1,N)=P^i·[T_i(k-1,N)-T_i(k-2,N)] .$$

(7)

Obviously, translational motion will not change the prediction matrix $P^i$. Therefore, it can be expressed as

$$T_i(k,N)-H=P^i·[T_i(k-1,N)-H] ,$$

(8)

and $H$ satisfies the constraint as follows

$$H={(P^i)}^n·H,$$

(9)

where $H\triangleq {[a,b,\dots ,a,b]}_{2N\times 1}^{\mathrm{T}}$ and its two base vectors are $h_1\triangleq {[1 0 1 0\dots 1 0]}^{\mathrm{T}}$ and $h_2\triangleq {[0 1 0 1\dots 0 1]}^{\mathrm{T}}$. Then we have $h_1=P^i·h_1$ and $h_2=P^i·h_2$.

Assume that each node stores M ($M>N$) estimated positions before the kth time slot, which constitute $M-N+1$ trajectories, i.e., $T_i(k-1,N)$, $T_i(k-2,N)$,…, $T_i(k-M+N-1,N)$. Based on the prediction model in Equation (3), we obtain an equation set which can be expressed as

$$\left\{\begin{array}{c}{T}_i(k-1,N)=P^i\times T_i(k-2,N)\\ T_i(k-2,N)=P^i\times T_i(k-3,N)\\ ⋮\\ T_i(k-M+N,N)=P^i\times T_i(k-M+N-1,N)\end{array}\right..$$

(10)

According to Equations (4) and (10), we have

$$\{\begin{array}{l}\left(\begin{array}{cc}1& 0\end{array}\right)=h_1^{\mathrm{T}}\times \left(\begin{array}{cc}{(p_{1,\ast }^i)}^{\mathrm{T}}& {(p_{2,\ast }^i)}^{\mathrm{T}}\end{array}\right)\\ \left(\begin{array}{cc}0& 1\end{array}\right)=h_2^{\mathrm{T}}\times \left(\begin{array}{cc}{(p_{1,\ast }^i)}^{\mathrm{T}}& {(p_{2,\ast }^i)}^{\mathrm{T}}\end{array}\right)\end{array},$$

(11)

where $p_{1,\ast }^i$ and $p_{2,\ast }^i$ are the first and second rows of $P^i$. For convenience of description, we define

$$B_{i,M-N}\triangleq \left[\begin{array}{cc}\begin{array}{c}{\widehat{x}}_i^{k-1}\\ \begin{array}{c}{\widehat{x}}_i^{k-2}\\ ⋮\end{array}\\ {\widehat{x}}_i^{k-M+N}\\ \begin{array}{c}1\\ 0\end{array}\end{array}& \begin{array}{c}{\widehat{y}}_i^{k-1}\\ \begin{array}{c}{\widehat{y}}_i^{k-2}\\ ⋮\end{array}\\ {\widehat{y}}_i^{k-M+N}\\ \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\right], T_{i,M-N}\triangleq \left[\begin{array}{c}{\left(T_i(k-2,N)\right)}^{\mathrm{T}}\\ \begin{array}{c}{\left(T_i(k-3,N)\right)}^{\mathrm{T}}\\ ⋮\end{array}\\ {\left(T_i(k-M+N-1,N)\right)}^{\mathrm{T}}\\ \begin{array}{c}{h}_1^{\mathrm{T}}\\ h_2^{\mathrm{T}}\end{array}\end{array}\right].$$

Then, according to Equations (4), (10), and (11), we have

$$B_{i,M-N}=T_{i,M-N}·\left(\begin{array}{cc}{(p_{1,\ast }^{\mathrm{i}})}^{\mathrm{T}}& {(p_{2,\ast }^{\mathrm{i}})}^{\mathrm{T}}\end{array}\right).$$

(12)

Consequently, $p_{1,\ast }^i$ and $p_{2,\ast }^i$ can be determined based on the least square method, i.e.,

$$\left(\begin{array}{l}{p}_{1,\ast }^i\\ p_{2,\ast }^i\end{array}\right)=B_{i,M-N}^{\mathrm{T}}·T_{i,M-N}·{\left[{(T_{i,M-N}^{\mathrm{T}}·T_{i,M-N})}^{-1}\right]}^{\mathrm{T}}.$$

(13)

Finally, the prediction matrix Pⁱ is given by Equations (4) and (13). Let ${\tilde{\mathit{x}}}_i^k\triangleq {[{\tilde{x}}_i^k,{\tilde{y}}_i^k]}^{\mathrm{T}}$ denote the predicted position of node i. According to the prediction model in Equation (3), ${\tilde{x}}_i^k$ can be given by

$${\tilde{x}}_i^k=\left(\begin{array}{l}{p}_{1,\ast }^{\mathrm{i}}\\ p_{2,\ast }^{\mathrm{i}}\end{array}\right)·T_i(k-1,N)$$

(14)

Then, we use the predicted position ${\tilde{x}}_i^k$ to formulate the prediction message $m_{f_i^{k|k-1}\to x_i^k}(x_i^k)$ as

$$m_{f_i^{k|k-1}\to {\mathit{x}}_i^k}\left({\mathit{x}}_i^k\right)=\mathcal{N}({\mathit{x}}_i^k;{\tilde{\mathit{x}}}_i^k,{\tilde{\mathit{V}}}_i^k),$$

(15)

where the covariance matrixes are ${\tilde{V}}_i^k={\tilde{\sigma }}_{i,k}^2{I}_{2\times 2}$ and ${\tilde{\sigma }}_{i,k}^2=||{\widehat{x}}_i^{k-1}-{\tilde{x}}_i^{k-1}|{|}^2$.

3.2. Joint Self-Location and Tracking Phase

After calculating the prediction message based on the adaptive prediction model, agent $m\in \mathcal{M}$ locates itself using ranging measurements, not only from its neighbor anchors and agents $j\in {\mathcal{A}}_i^k\cup {\mathcal{M}}_i^k$, but also from targets $O\in {\mathcal{S}}_i^k$ that it can detect. Meanwhile, agents and anchors that can detect target O share information with each other based on consensus iteration and determine the position of the target O. Compared with the separate cooperative self-localization and algorithms, agents use the ranging measurements from targets to improve the self-location accuracy, which promotes tracking performance in turn.

3.2.1. Self-Location of Agents Based on Variational Message Passing

At the kth time slot, agent $m\in \mathcal{M}$ obtains its prediction message $m_{f_m^{k|k-1}\to {\mathit{x}}_m^k}\left({\mathit{x}}_m^k\right)=\mathcal{N}({\mathit{x}}_m^k;{\tilde{\mathit{x}}}_m^k,{\tilde{\mathit{V}}}_m^k)$ based on the adaptive prediction model described in Section 3. Then, according to VMP update rules, cooperative messages from its neighbors are expressed as follows:

$$\begin{array}{cc}\hfill m_{f_{ma}^k\to {\mathit{x}}_m^k}\left({\mathit{x}}_m^k\right)& \triangleq exp\left\{\int m_{{\mathit{x}}_a^k\to f_{ma}^k}\left({\mathit{x}}_a^k\right)ln{f}_{ma}^k({\mathit{x}}_a^k,{\mathit{x}}_m^k)\mathrm{d}{\mathit{x}}_a^k\right\}\hfill \\ & =\mathcal{N}\left(d_{m\leftarrow a}^k;∥{\mathit{x}}_m^k-{\mu }_a^k∥,{\left({\sigma }_{ma}^k\right)}^2\right),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill m_{f_{mj}^k\to {\mathit{x}}_m^k}\left({\mathit{x}}_m^k\right)& =exp\left\{\int m_{{\mathit{x}}_j^k\to f_{mj}^k}lnp\left(d_{m\leftarrow j}^k\right|{\mathit{x}}_m^k,{\mathit{x}}_j^k)\mathrm{d}{\mathit{x}}_j^k\right\}\hfill \\ & =exp\left\{\int b\left({\mathit{x}}_j^k\right)ln\mathcal{N}\left(d_{m\leftarrow j}^k;∥{\mathit{x}}_m^k-{\mathit{x}}_j^k∥,{\left({\sigma }_{mj}^k\right)}^2\right)\mathrm{d}{\mathit{x}}_j^k\right\}.\hfill \end{array}$$

(17)

In (16), ${\mu }_a^k$ is the position of anchor $a\in \mathcal{A}$ at the kth time slot. Consequently, the belief $b(x_m^k)$ in Equation (2) can be rewritten as

$$\begin{array}{cc}\hfill b\left({\mathit{x}}_m^k\right)& =\frac{1}{Z}\mathcal{N}\left({\mathit{x}}_m^k;{\tilde{\mathit{x}}}_m^k,{\tilde{\mathit{V}}}_m^k\right)\times \prod _{a\in {\mathcal{A}}_m^k}\mathcal{N}\left(d_{m\leftarrow a}^k;∥{\mathit{x}}_m^k-{\mu }_a^k∥,{\left({\sigma }_{ma}^k\right)}^2\right)\hfill \\ & \times \prod _{j\in {\mathcal{M}}_m^k\cup {\mathcal{S}}_m^k}exp\left\{\int b\left({\mathit{x}}_j^k\right)ln\mathcal{N}\left(d_{m\leftarrow j}^k;∥{\mathit{x}}_m^k-{\mathit{x}}_j^k∥,{\left({\sigma }_{mj}^k\right)}^2\right)\mathrm{d}{\mathit{x}}_j^k\right\}.\hfill \end{array}$$

(18)

Considering the communication overhead, the messages passed on in Figure 1 are restricted to Gaussian distribution. However, $b(x_m^k)$, given by Equation (18), is non-Gaussian because of the nonlinear ranging model. In [20], we propose a low-complexity approximation method based on second-order Taylor expansion. In this way, the belief $b(x_m^k)$ is approximated into a Gaussian function denoted by $\widehat{b}(x_m^k)\sim \mathcal{N}\left(x_m^k;{\widehat{\mu }}_m^k,{\widehat{V}}_m^k\right)$, where ${\mu }_i^k$ and $V_i^k$ are as follows:

$$V_m^k\triangleq {\left\{{({\tilde{V}}_m^k)}^{-1}+{\displaystyle \sum _{a\in {\mathcal{A}}_m^k}\frac{I_{2\times 2}-d_{m\leftarrow a}^k{\nabla }_{F_{ma}}^2}{{({\sigma }_{ma}^k)}^2}}+{\displaystyle \sum _{j\in {\mathcal{M}}_i^k}\frac{1}{{({\sigma }_{mj}^k)}^2}[I_{2\times 2}-d_{m\leftarrow j}^k\frac{{\partial }^2{F}_{mj}}{\partial {(x_m^k)}^2}]}\right\}}^{-1},$$

(19)

$${\mu }_m^k\triangleq {\widehat{V}}_m^k\left\{{({\tilde{V}}_m^k)}^{-1}{\tilde{x}}_m^k\right. +{\displaystyle \sum _{a\in {\mathcal{A}}_m^k}\frac{{\mu }_a^k-d_{i\leftarrow a}^k({\nabla }_{F_{ma}}-{\nabla }_{F_{ma}}^2{\mu }_m^{k\ast })}{{({\sigma }_{ma}^k)}^2}} +\left.{\displaystyle \sum _{j\in {\mathcal{M}}_m^k\cup {\mathcal{S}}_m^k}\frac{1}{{({\sigma }_{mj}^k)}^2}[{\mu }_j^{k\ast }+d_{m\leftarrow j}^k(\frac{\partial F_{mj}}{\partial x_m^k}-\frac{{\partial }^2{F}_{mj}}{\partial {(x_m^k)}^2}{\mu }_m^{k\ast })]}\right\}.$$

(20)

In Equations (19) and (20), the notations ${\mu }_m^{k\ast }$ and ${\mu }_j^{k\ast }$ represent the estimated positions of nodes m and j in the latest iteration, respectively, notations ${\nabla }_{F_{ma}}$ and ${\nabla }_{F_{ma}}^2$ are the first-order gradient and the Hessian matrix of $F_{ma}^k\triangleq \left|\right|{\mathit{x}}_m^k-{\mu }_a^k\left|\right|$ at ${\mu }_m^{k\ast }$, respectively, and the notations $\frac{\partial F_{mj}}{\partial x_m^k}$ and $\frac{{\partial }^2{F}_{mj}}{\partial {(x_m^k)}^2}$ are the first-order partial derivatives of $F_{mj}^k\triangleq ||x_m^k-x_j^k||$ and $\frac{1}{\partial x_m^k}(\frac{\partial F_{mj}}{\partial x_m^k})$ at (${\mu }_m^{k\ast }$,${\mu }_j^{k\ast }$), respectively (please see [22] for detailed derivation). Consequently, agent $i\in \mathcal{M}$ obtains its location based on maximum a posteriori (MAP), i.e., ${\widehat{x}}_i^k={\mu }_i^k$.

3.2.2. Target Tracking Based on Average Consensus

At kth time slot, node $i\in {\mathcal{S}}_o^k$, which can detect target $O\in \mathcal{S}$, computes the prediction message $m_{f_o^{k|k-1}\to {\mathit{x}}_o^k}\left({\mathit{x}}_o^k\right)=\mathcal{N}({\mathit{x}}_o^k;{\tilde{\mathit{x}}}_o^k,{\tilde{\mathit{V}}}_o^k)$ and the local message $m_{f_{oi}^k\to {\mathit{x}}_o^k}\left({\mathit{x}}_o^k\right)$. Then, $m_{f_{oi}^k\to {\mathit{x}}_o^k}\left({\mathit{x}}_o^k\right)$ is approximated into Gaussian function $\mathcal{N}({\mathit{x}}_o^k;{\mu }_{oi}^k,{\mathit{V}}_{oi}^k)$ using the approximation method in [22]. Therefore, the Gaussian approximation of $b\left({\mathit{x}}_o^k\right)$ can be expressed as

$$\widehat{b}\left({\mathit{x}}_o^k\right)=\frac{1}{Z}\mathcal{N}({\mathit{x}}_o^k;{\tilde{\mathit{x}}}_o^k,{\tilde{\mathit{V}}}_o^k)\times \prod _{i\in {\mathcal{V}}_o^k}\mathcal{N}({\mathit{x}}_o^k;{\mu }_{oi}^k,{\mathit{V}}_{oi}^k).$$

(21)

The mean vector ${\mu }_o^k$ and the covariance matrix $V_o^k$ are given by ${(V_o^k)}^{-1}={({\tilde{V}}_o^k)}^{-1}+{\displaystyle \sum _{i\in {\mathcal{V}}_o^k}{(V_{oi}^k)}^{-1}}$ and ${\mu }_o^k={\mathit{V}}_o^k\left({\left({\tilde{\mathit{V}}}_i^k\right)}^{-1}{\mu }_i^k+\sum _{i\in {\mathcal{V}}_o^k}{\left({\mathit{V}}_{oi}^k\right)}^{-1}{\mu }_{oi}^k\right)$, respectively. Each node $i\in {\mathcal{S}}_o^k$ can calculate ${\tilde{V}}_o^k$ and $V_{oi}^k$ locally. However, node $i\in {\mathcal{S}}_o^k$ cannot determine $\widehat{b}(x_o^k)$ unless all other neighbors of target $O$ share their local information with it.

For convenience of discussion, let $W_o^k\triangleq {(V_o^k)}^{-1}$, ${\tilde{W}}_o^k\triangleq {({\tilde{V}}_o^k)}^{-1}$ and $W_{oi}^k\triangleq {(V_{oi}^k)}^{-1}$; therefore, we have

$$W_o^k={\tilde{W}}_o^k+{\displaystyle \sum _{i\in {\mathcal{V}}_o^k}{W}_{oi}^k},$$

(22)

$${\mathit{W}}_o^k{\mu }_o^k={\tilde{\mathit{W}}}_o^k{\tilde{\mathit{x}}}_o^k+\sum _{i\in {\mathcal{V}}_o^k}{\mathit{W}}_{oi}^k{\mu }_{oi}^k.$$

(23)

In order to locate the target in a distributed manner, we use the average consensus based on Metropolis weights proposed in [13] to diffuse information across the network. During the process, each node updates its data with a weighted average of its neighbors’ data, which converges to the global average. Let t denote the iteration index and $N_i^k$ represent the number of neighbors of node i, respectively. The consensus iteration process is as follows.

Step 1: Each node $i\in {\mathcal{V}}_o^k$ uses local information to initialize the global average ${\mathit{W}}_i^{k,\left(0\right)}$ and ${\mathit{W}}_i^{k,\left(0\right)}{\mu }_i^{k,\left(0\right)}$ as

$${\mathit{W}}_i^{k,\left(0\right)}={\mathit{W}}_{oi}^k,{\mathit{W}}_i^{k,\left(0\right)}{\mu }_i^{k,\left(0\right)}={\mathit{W}}_{oi}^k{\mu }_{oi}^k,$$

while each node $i\notin {\mathcal{V}}_o^k$ initializes ${\mathit{W}}_i^{k,\left(0\right)}$ and ${\mathit{W}}_i^{k,\left(0\right)}{\mu }_i^{k,\left(0\right)}$ as

$${\mathit{W}}_i^{k,\left(0\right)}=0,{\mathit{W}}_i^{k,\left(0\right)}{\mu }_i^{k,\left(0\right)}=0.$$

Step 2: Each node computes the Metropolis weight matrix, which is defined as

$${\xi }_{i,j}^{\left(t\right)}=\left\{\begin{array}{cc}\frac{1}{1+max\{N_i^k,N_j^k\}}\hfill & \mathrm{if} j\in {\mathcal{V}}_i^k,\hfill \\ 1-\sum _{l\in {\mathcal{V}}_i^k}{\xi }_{i,l}^{\left(t\right)}\hfill & \mathrm{if} i=j,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(24)

Step 3: Each node broadcasts its global average computed in the latest iteration and collects the information broadcasted by its neighbors. Then, each node updates the global average $W_i^{k,(t+1)}$ and $W_i^{k,(t+1)}{\mu }_i^{k,(t+1)}$ by

$${\mathit{W}}_i^{k,(t+1)}={\xi }_{i,i}^{\left(t\right)}{\mathit{W}}_i^{k,\left(t\right)}+\sum _{j\in {\mathcal{V}}_i^k}{\xi }_{i,j}^{\left(t\right)}{\mathit{W}}_j^{k,\left(t\right)},$$

(25)

$${\mathit{W}}_i^{k,(t+1)}{\mu }_i^{k,(t+1)}={\xi }_{i,i}^{\left(t\right)}{\mathit{W}}_i^{k,\left(t\right)}{\mu }_i^{k,\left(t\right)}+\sum _{j\in {\mathcal{V}}_i^k}{\xi }_{i,j}^{\left(t\right)}{\mathit{W}}_j^{k,\left(t\right)}{\mu }_j^{k,\left(t\right)}.$$

(26)

Step 4: If $t=t_{\max }$ (where $t_{\max }$ is the maximum number of iteration), the consensus process ends. Otherwise, go to Step 2.

After sufficient consensus iterations, all nodes obtain the global average of ${\displaystyle \sum _{i\in {\mathcal{V}}_o^k}{W}_{oi}^k}$ and ${\displaystyle \sum _{i\in {\mathcal{V}}_o^k}}{\mathit{W}}_{oi}^k{\mu }_{oi}^k$. Therefore, each node $i\in {\mathcal{V}}_o^k$ can compute ${\displaystyle \sum _{i\in {\mathcal{V}}_o^k}}{\mathit{W}}_{oi}^k$ and ${\displaystyle \sum _{i\in {\mathcal{V}}_o^k}{W}_{oi}^k {\mu }_{oi}^k}$ by

$$\sum _{i\in {\mathcal{V}}_o^k}{\mathit{W}}_{oi}^k=N_o^k\times {\mathit{W}}_i^{k,\left(t_{max}\right)},$$

(27)

$$\sum _{i\in {\mathcal{V}}_o^k}{\mathit{W}}_{oi}^k{\mu }_{oi}^k=N_o^k\times {\mathit{W}}_i^{k,\left(t_{max}\right)}{\mu }_i^{k,\left(t_{max}\right)}.$$

(28)

Submitting Equations (27) and (28) to Equations (22) and (23), each node $i\in {\mathcal{V}}_o^k$ has

$$W_o^k={\tilde{W}}_o^k+N_o^k\times W_i^{k,(t_{\max })},$$

(29)

$$W_o^k{\mu }_o^k={\tilde{W}}_o^k{\tilde{x}}_o^k+N_o^k\times W_i^{k,(t_{\max })}{\mu }_i^{k,(t_{\max })}.$$

(30)

Therefore, each node $i\in {\mathcal{V}}_o^k$ can determine the position of target $O$ based on MAP, i.e., ${\widehat{x}}_o^k={\mu }_o^k$.

3.3. Implementation Scheme of the Joint Self-Location and Tracking Algorithm

The proposed distributed joint self-location and tracking algorithm based on VMP and average consensus is stated as follows.

Algorithm 1: Distributed joint self-location and tracking algorithm
1	for $k=1:k_{\max }$(k is the time slot index)
2	Each agent $m\in \mathcal{M}$ updates M historical position information of itself and the targets before the kth time slot.
3	Each agent $m\in \mathcal{M}$ calculates its prediction matrix $P^m$ using Equation (4) and Equation (13), then predicts its position ${\tilde{x}}_m^k$ using Equations (14) and obtains its prediction message $m_{f_m^{k|k-1}\to x_m^k}(x_m^k)$ using Equation (15).
4	Each node $i\in {\mathcal{S}}_o^k$ predicts the position of target $O$ using Equation (4), Equation (13), and Equation (14), i.e., ${\tilde{x}}_o^k$. Then, the prediction message $m_{f_o^{k|k-1}\to x_o^k}(x_o^k)$ is calculated using Equation (15).
5	Initialization: Each agent $m\in \mathcal{M}$ initializes its position as $x_m^k={\tilde{x}}_m^k$ and initializes the position of the target that it detects as $x_o^k={\tilde{x}}_o^k$.
6	for $q=1:q_{\max }$(where q is the iteration index)
7	(1) Each agent $m\in \mathcal{M}$ collects ranging measurements from its neighbors and calculates cooperative messages $m_{f_{ma}^k\to {\mathit{x}}_m^k}^{\left(q\right)}\left({\mathit{x}}_m^k\right)$($\forall a\in {\mathcal{A}}_m^k$) from neighbor anchors using Equation (16) and $m_{f_{mj}^k\to {\mathit{x}}_m^k}^{\left(q\right)}\left({\mathit{x}}_m^k\right)$($\forall j\in {\mathcal{M}}_m^k\cup {\mathcal{S}}_m^k$) from neighbor agents and targets using Equation (17), respectively. Then, each agent $m\in \mathcal{M}$ calculates ${\mathit{\mu }}_m^{k,\left(q\right)}$ and ${\mathit{V}}_m^{k,\left(q\right)}$ of ${\widehat{b}}^{\left(p\right)}\left({\mathit{x}}_m^k\right)$ using Equations (19) and (20), respectively.
8	(2) Each node $i\in {\mathcal{S}}_o^k$ receives the local message $m_{f_{oi}^k\to x_o^k}^{(q)}(x_o^k)$ and approximates it into Gaussian function $\mathcal{N}(x_o^k;{\mu }_{oi}^{k,(q)},V_{oi}^{k,(q)})$ based on the Taylor expansion approximation method in [20]. Then the global average is calculated based on average consensus iteration from Step 1 to Step 4, which are depicted in Section 3.2.2. Then, each node $i\in {\mathcal{S}}_o^k$ computes $W_o^{k,(q)}$ and $W_o^{k,(q)}{\mu }_o^{k,(q)}$ using Equations (27)–(30).
9	(3) Each agent $m\in \mathcal{M}$ broadcasts the mean ${\mu }_m^{k,(q)}$ and the covariance matrix $V_m^{k,(q)}$ of ${\widehat{b}}^{(p)}(x_m^k)$. Each node $i\in {\mathcal{S}}_o^k$ broadcasts $W_o^{k,(q)}{\mu }_o^{k,(q)}$ and $W_o^{k,(q)}$. At the same time, they collect the information broadcasted by its neighbors.
10	end for q
11	Each agent $m\in \mathcal{M}$ determines its position based on MAP and broadcasts ${\mu }_m^k$ and $V_m^k$. Each node $i\in {\mathcal{S}}_o^k$ determines the position of target O based on based on MAP and broadcast $W_o^k$ and $W_o^k{\mu }_o^k$.
12	end for k

4. Stimulation Results and Analysis

In this section, we evaluate the performance of the proposed adaptive prediction model and the distributed joint self-location and tracking algorithm (denoted by JSLT).

4.1. Stimulation Results and Performance Analysis of the Proposed Adaptive Prediction Model

In order to assess the performance of the proposed adaptive prediction model, we compare it with linear prediction and square prediction. We consider a 50 m × 50 m plane with 10 moving nodes. The prediction error is defined as the distance between the true position and the prediction position. The average prediction errors in Figure 2, Figure 3 and Figure 4 are obtained by averaging over 10,000 realizations, with each realization continuing 20 time slots.

We employ the Gauss–Markov (GM) mobility model [24] to simulate the motion of nodes. In the GM model, the velocity and direction of nodes are time-dependent Gauss–Markov processes, which means the velocity and direction at the kth time slot is related to that at the (k−1)th time slot. Let $d_i^k$ and ${\theta }_i^k$ denote the velocity and direction of node i from the kth time slot to the (k−1)th time slot, respectively, then we have

$$\left\{\begin{array}{l}{d}_i^k=\alpha d_i^{k-1}+(1-\alpha ){\overline{d}}_i+\sqrt{1-{\alpha }^2}{\omega }_d\\ {\theta }_i^k=\alpha {\theta }_i^{k-1}+(1-\alpha ){\overline{\theta }}_i+\sqrt{1-{\alpha }^2}{\omega }_{\theta }\end{array}\right.,$$

(31)

where constants ${\overline{d}}_i$ and ${\overline{\theta }}_i$ are the velocity and direction of node i when $k\to \infty$, ${\omega }_d$ and ${\omega }_{\theta }$ are Gaussian random variables independent of ${\overline{d}}_i$ and ${\overline{\theta }}_i$, and $\alpha \in [0,1]$ is an adjustment parameter to control the randomness of motion. When $\alpha$ increases, the impact of $d_i^{k-1}$ and ${\theta }_i^{k-1}$ rises, while the impact of ${\omega }_d$ and ${\omega }_{\theta }$ declines. As discussed above, the Gauss–Markov mobility model not only reflects the relationship between the current moment and a historical moment, but also reflects the randomness of motion. Therefore, it is very suitable for the simulation of high random motion, such as the travel of birds and the diffusion of gas.

The performance of the proposed adaptive prediction model with different numbers of historical estimated positions and $\alpha$ is illustrated in Figure 2. Simulation configurations are as follows. For the GM mobility model shown in Equation (31), we set ${\overline{d}}_i=5 \mathrm{m}/\mathrm{s}$ and ${\overline{\theta }}_i={60}^{°}$. For the prediction using Equation (13), the parameters are $N=1$ and $M=3/4/5/6$. As we know, the prediction algorithm is based on the historical estimated positions. In initialization, we use the true positions along with Gaussian noise (mean = 0 and variance = 4 m²) as the estimated positions. As shown in Figure 2, there is little difference in the average prediction error using different historical positions with different $\alpha$. Comparatively speaking, when the value of $\alpha$ is smaller, more historical positions are used and the more accurate the obtained prediction is; when the value of $\alpha$ is bigger, more historical positions are used and a more accurate prediction is achieved. This is because a smaller $\alpha$ value indicates a greater randomness of motion, therefore the impact of the (k−1)th time slot on the kth time slot is less. More historical positions can reflect the trend of the movement. However, a bigger $\alpha$ value leads to a stronger correlation between the kth time slot and the (k−1)th time slot. In this situation, the contribution of the earlier time slots to the current moment is small, and using too much historical information even influences the accuracy of prediction. Moreover, greater the number of historical positions used, the more information the nodes store and the greater the energy consumption. In subsequent simulations, we use three historical positions, meaning each node stores three historical positions.

To further evaluate the performance of the proposed adaptive prediction model, we compare the proposed model with linear prediction, square prediction, and instant prediction, which use the movements in x axis and y axis at the (k−1)th time slot to predict the position at the kth time slot. Simulation configurations are $\alpha =0.7$, ${\overline{d}}_i=3 \mathrm{m}/\mathrm{s}$, ${\overline{\theta }}_i={60}^{°}$, $N=1$, and $M=3$. The performance comparison of the four prediction models with different standard deviations ${\sigma }_{{\omega }_d}={\sigma }_{{\omega }_{\theta }}=\sigma$ of the random variables ${\omega }_d$ and ${\omega }_{\theta }$ is shown in Figure 3. It is observed that the average prediction error of the four prediction models increases with an increase in $\sigma$. Specifically, the performance of instant prediction is the worst and the performance of the proposed adaptive prediction is the best. Moreover, the advantage of the proposed adaptive prediction model is more obvious when $\sigma$ is bigger. Figure 4 illustrates the performance of the four prediction models with different adjustment parameter $\alpha$ values. Simulation configurations are $\alpha =0.7$, ${\overline{d}}_i=3 \mathrm{m}/\mathrm{s}$, ${\overline{\theta }}_i={60}^{°}$, $\sigma =1 \mathrm{m}$ and $M=3$. It can be observed that the average prediction error of the proposed adaptive prediction model is smaller than the other three prediction models and the difference is bigger with increasing $\alpha$. In fact, from the Gauss–Markov mobility model in Equation (31), we see that the randomness of motion rises with the decline of $\alpha$ and the increase in $\sigma$. Therefore, the proposed adaptive prediction model better adapts to different motion patterns than the other three prediction models.

4.2. Stimulation Results and Performance Analysis of the Proposed Distributed Joint Self-Location and Tracking Algorithm

The proposed distributed joint self-location and tracking (JSLT) algorithm is evaluated by comparison with the distributed cooperative self-localization algorithm and the centralized target tracking algorithm in this section. We consider a 50 m × 50 m plane with $N_A=5$ anchors, $N_M$ agents, and $N_O$ targets. The communication radius and ranging radius of anchors and agents are $R_{\mathrm{c}}=25 \mathrm{m}$, the detection radius is $R_{\mathrm{d}}=25 \mathrm{m}$, and the standard variance of ranging measurement noise is 1 m. At the 0th time slot, the variance matrix of the prior of agents and targets are set to be $25I_{2\times 2}$. The iterations of the average consensus algorithm and the joint location algorithm are 20 and 10, respectively. The results in Figure 5, Figure 6 and Figure 7 are obtained by performing 10,000 realizations over 20 time slots.

Self-localization performance of the proposed JSLT algorithm and the cooperative self-localization (CSL) algorithm in [23] with different numbers of agents and targets are shown in Figure 5 and Figure 6. In Figure 5, $N_M=20$, $N_A=1$ or 10; in Figure 6, $N_M=10$, $N_A=1$ or 5. It can be seen that the localization performance of both algorithms improves as time increases. This is because the agents’ and their neighbors’ positions become more accurate with the increase in time. However, a higher k value does not always correspond to improved positioning accuracy. The estimated accuracy becomes stable as time goes on. Furthermore, in the two scenarios, the proposed JSLT algorithm has a higher positioning accuracy than the CSL algorithm. This is because, in the JSLT algorithm, agents not only use the distance observation from the neighboring anchors and agents as in the CSL algorithm, but they also make use of the ranging measurements from the targets that they can detect, which improves the self-localization accuracy. Moreover, the performance of the two algorithms is quite different with fewer agents and more targets. When the density of agents is high, agents have enough neighbors. Therefore, they can localize themselves well by cooperating with neighboring anchors and agents. However, when the density of agents is low, many agents have few neighboring anchors and agents to localize themselves. In this case, the ranging measurements from targets play an important role in the self-localization of agents. The computational complexity of both algorithms depends on the number of neighbors. For agent $m\in \mathcal{M}$, we define $N_m^A$, $N_m^M$, and $N_m^O$ to denote the number of neighboring agents, anchors, and targets that can detected at a certain time. The computational complexity of the CSL algorithm and the proposed JSLT algorithm are $\mathcal{O}(N_i^A+N_i^M)$ and $\mathcal{O}(N_i^A+N_i^M+N_i^O)$, respectively. With regard to communication overhead, as all messages are approximated to Gaussian function, each agent broadcasts the mean and the covariance to its neighbors and collects the information broadcasted by its neighbors. The information collected by each agent is $\mathcal{O}(N_i^A+N_i^M)$ and $\mathcal{O}(N_i^A+N_i^M+N_i^O)$ for the CSL algorithm and the proposed JSLT algorithm, respectively. Since $N_i^O\ll N_i^A+N_i^M$, the complexity and the communication overhead of the two algorithms are similar.

Figure 7 shows the tracking performance of the proposed JSLT algorithm and the centralized tracking (CT) algorithm is illustrated with $N_M=20$ and $N_A=1$. In the CT algorithm, agents first localize themselves using the CSL algorithm in [23], then anchors and agents track the targets by sending information to a control center. As shown in Figure 7, the tracking performance of the proposed distributed JSLT algorithm is very close to that of CT algorithm, but the communication overhead of JSLT algorithm is much less than that of CT algorithm. In the CT algorithm, the information collected and broadcasted by each node (anchor or agent) are both proportional to the number of neighbors. In the proposed JLST algorithm, each node broadcasts its global average to its neighbors and collects global averages broadcasted by its neighbors, which is proportional to the number of neighbors. Therefore, the communication overhead sharply reduces by half. Moreover, in the CT algorithm, anchors and agents transmit information to the control center in a multi-hop manner, which leads to a fair amount of redundant information in the networks and information delay.

In order to visualize the tracking performance, a single trial performance of target tracking is shown in Figure 8. We can see that the trajectory recognized by the proposed algorithm is close to the actual trajectory of the target, except for at the initial time and at the moment of a sharp turn.

5. Conclusions

In this paper, we built an adaptive prediction model by exploiting the correlation of trajectories and then proposed a joint self-location and tracking algorithm based on VMP and average consensus. Simulation results and analysis illustrated that, compared with instant prediction, linear prediction, and square prediction, the proposed prediction model was more adaptive to the movement of the nodes. The self-location performance of the proposed algorithm is better than the separated cooperative self-localization, and the tracking performance of the proposed algorithm is close to the centralized tracking algorithm with less communication overhead.

The proposed algorithm is suitable for tracking bird migration patterns, monitoring gas diffusion, etc. In these applications, the motion of the targets is random. The following research aimed to apply the proposed approach to practices. In practical networks, delays in data transmission and communication errors are two important issues which affect positioning accuracy and real-time tracking. Future research will focus on the issue of how to deal with the impact of these two issues.