Sequential Multidimensional Scaling with Kalman Filtering for Location Tracking

Abstract

Localization always plays a critical role in wireless sensor networks for a wide range of applications including military, healthcare, and robotics. Although the classical multidimensional scaling (MDS) is a conventionally effective model for positioning, the accuracy of this method is affected by noises from the environment. In this paper, we propose a solution to attenuate noise effects to MDS by combining MDS with a Kalman filter. A model is built to predict the noise distribution with regard to additive noises to the distance measurements following the Gaussian distribution. From that, a linear tracking system is developed. The characteristics of the algorithm are examined through simulated experiments and the results reveal the advantages of our method over conventional works in dealing with the above challenges. Besides, the method is simplified with a linear filter; therefore it suits small and embedded sensors equipped with limited power, memory, and computational capacities well.

1. Introduction

The sensor, a device that is able to measure and respond to physical input from the environment, has important roles that emerged in our daily life. In addition, a set of sensors embedded with communication capability can form a wireless sensor network to perform complicated tasks. One of the most important applications involves military areas where sensor networks are used for spying, map exploration, and topographic reconstruction. To do so, localization is critical to indicate the position of each sensor with others in a global context. Usually a wireless sensor network is a set of up to hundreds of nodes, in which each node is able to sense the surrounding environment, to perform computations and to communicate with other nodes to accomplish a task. Node localization is needed to report the origin of events in the network, to assist group querying of sensors and geographical routing, and to address network coverage. Therefore, knowing node location is crucial.

For decades, localization has mainly been based on Global Positioning System (GPS). However, GPS devices are too big to be fitted into tiny sensors. In addition, GPS cannot be implemented in indoor localization and consumes a lot of energy to be operated. A great deal of efforts over the past few years has focused on the indoor localization problem and on energy saving for small and embedded devices of sensors. A more innovative technique is angle-of-arrival (AOA) that estimates the direction of the arrival signal to point out the location. However, these techniques need antenna systems, resulting in highly complex hardware requirements. Additional drawbacks are connected to errors caused by multipath and non-line-of-sight phenomena.

The practical approach suitable for sensor nodes deals with locating nodes by measuring distances between them. The survey of existing approaches is introduced in [1]. There are many solutions available for finding the distance between two sensors. For wireless devices, the general technique involves the received signal strength indicator (RSSI) method that measures the power of the signal at receiver to locate mobile station. Other approaches, including time-of-arrival (TOA) and time-difference-of-arrival (TDOA), analyze the differences between the emitted and returned light from a source to a destination. In recent works [2, 3], non-line-of-sight errors are significantly mitigated with TOA measurement, making this approach promising for practical applications. Basically, the position of an object is determined by converting the obtained measurements to draw circles centered at known-base stations and estimating the intersection of the circles. The intersection of three circles is sufficient to locate a point corresponding to the sensor node in two dimensions (2D). It is noted that the more the distances are discovered from a sensor to its neighbors, the higher the chances will be that the position of that sensor is located exactly. In general, all pairwise distances of all sensors should be used. An effective method to recover the localization maps of sensor nodes using all pairwire distances in a network copes with the multidimensional scaling (MDS). MDS refers to a family of approaches employed for exploratory analysis of multidimensional data. The MDS method was originally proposed by Torgerson [4] where the algorithm was based on metric analyses and the law of cosine to place objects in the Euclidean coordinate system that reflects measured distances or dissimilarity between them. The technique was also applied to represent multidimensional data, especially images in a low-dimensional space for visualization. Till then many researchers attempted to introduce MDS in numerous applications like environmental monitoring, surveillance, transportation, and healthcare. Another important role we want to focus on in this paper is to address localization problems for sensor network [5] and mobile devices [6, 7] with MDS. Agrawal and Patel [5] took advantage of the MDS algorithm for localization in wireless sensor networks that consist of a large number of nodes. The approach was analyzed with numerous simulations to address shortcomings caused by anisotropic network topology and complex terrain. Yet the analysis of the communication costs, messaging complexity, and power consumption was lacking. With regard to the sensor power limits to exploit the distance of nodes in a large-scale network, in [8] the MDS-MAP algorithm has been introduced to approximate the Euclidean distance by the shortest path between nodes. Based on the known local connectivity, a 2D map is constructed and the shortest path is found with Dijkstra's algorithm. In addition, MDS is also applicable for 3D positioning through distance measurement [9–11].

Based on MDS and its variant MDS-MAP, improved MDSs have been proposed. Wu et al. [12] developed a dynamic mobility-assisted MDS localization technique for a mobile sensor network. This method highly depends on the degree of a network since virtual nodes are added to acquire precise locations, causing more data to be transferred over a network and increasing computation cost. In [13], the classical MDS was combined with a neural network to improve location accuracy. Hierarchical MDS (HMDS) [14], ordinal MDS [15], universal MDS [16], improved MDS [17, 18], and deterministic annealing MDS [19] focus on the optimization of MDS. Cluster MDS [20, 21] divided a large network into subgroups so that the estimation on the big matrices is mitigated. Distributed MDS [22] helped share the computational workloads with different sensors in the network. In another aspect, there has been increasing interest on extending the classical MDS algorithm to deal with the environment artifacts on distance measurements.

Previously, filtering algorithms have been developed to attenuate measurement errors caused by tracking equipment such as radars and biosensors. However, there is a lack of research focusing on tracking nodes when the distance measurements are available. The main obstacle involves the nonlinear equation expressed to estimate location with respect to the distance between nodes. One of the methods of avoiding nonlinear calculation in tracking involves nonparametric method of the particle filter (PF) [23] in which the node distribution is described by a set of particles. A probabilistic model was established to constrain the distance between nodes. To perform the inference with both unitary and pairwise elements in the probabilistic model, the mean field was used [24]. Since mean field is variational inference, it requires a loop procedure to approximate and to replace the pairwise elements with the unitary ones. The computational time needed for particle filter is highly dependent on the convergence of the mean field approximation. Another kind of nonparametric inference is based sampling process of Markov chain Monte Carlo [25, 26]. Yet, sampling data from a highly multivariate distribution of sensor nodes is complicated. Concerning the parametric based methods, in [27], the nonlinear equation presenting the relationship between the locations of nodes and their distance measurement is linearized by partial derivatives. Consequently extended Kalman (EKL) filter is applied for localization. However, with regard to a large number of nodes, EKL easily leads to a suboptimal solution. Alternatively, the Kalman (KL) filter [23] constructs a linear filter for the tracking problem. Since MDS presents a good model for positioning nodes with distances, the linear KL filter combined with MDS tracking is able to extend tracking with independent nodes to a model where the distances of nodes are taken into account. Besides, the effective implementation of a linear filter makes it suitable for broader practical applications.

In this paper, we propose a linear KL filter integrated with MDS for tracking the location of nodes. We show an effective way to find the changing location of a node when the corresponding pairwise distance is measured; from that the KL filter is designed. The paper is organized as follows. Section 2 presents the essentials of classical MDS, KL filter, EKL filter, and our proposed tracking filter. Section 3 describes experimental results. Finally, we conclude our work and present further discussion in Section 4.

2. Methods

2.1. Classical MDS for Localization

So far, MDS [5–7, 28] has been proposed to establish a linear transformation to locate a set of nodes in multidimensional space to fit the similarity between nodes. Suppose that there are N-sensor nodes in a network noted by $x_1,x_2,\dots$ and let $d_{ij}$ be the dissimilarity measure between ith and jth nodes denoted by the $N\times N$ pairwise dissimilarity matrix of $d_{ij}$ . For the sensor localization problem, the dissimilarity is presented by the Euclidean distance between two points $x_i={[x_{i1},x_{i2},\dots ,x_{im}]}^T$ and $x_j={[x_{j1},x_{j2},\dots ,x_{jm}]}^T$ :

$$\begin{array}{cc}& d_{ij}=‖x_i-x_j‖=\sqrt{{\displaystyle \sum }_{k=\mathrm{1}}^m{\left(x_{ik}-x_{jk}\right)}^{\mathrm{2}}},\end{array}$$

(1)

where m is the number of dimensions. The target of classical MDS is to find an assignment of nodes $[x_1,x_2,\dots ,x_N]$ , $x_i={[x_{i1},x_{i2},\dots ,x_{im}]}^T$ , in m-dimensional space that minimizes a target function $F(\mathit{Y})$ expressed as

$$\begin{array}{cc}& F\left(\mathit{Y}\right)=\frac{1}{2}{\displaystyle \sum }_i^N{\displaystyle \sum }_j^N{\left(d_{ij}-‖y_i-y_j‖\right)}^2\\ \mathit{X}=\underset{\mathit{Y}}{\mathrm{argmin}}\hspace{0.17em}F\left(\mathit{Y}\right),\end{array}$$

(2)

where $y_i={[y_{i1},y_{i2},\dots ,y_{im}]}^T$ is a vector in m-dimensional space. In MDS, the elements of X are calculated in two stages including estimating the double centered squared distance matrix and applying eigenvalue decomposition on this double centered squared distance matrix. First, the double centered squared distance matrix is denoted by a matrix $\mathit{B}=-\left(1/2\right){\mathit{C}}_N\mathit{D}{\mathit{C}}_N$ , where ${\mathit{C}}_N$ is a centering operator and is computed as ${\mathit{C}}_N={\mathit{I}}_N-\left(1/N\right)\mathit{1}{\mathit{1}}^T$ , $\mathit{1}$ is the column-vector of N ones, and ${\left[\mathit{D}\right]}_{ij}=d_{ij}^2$ . The elements of matrix B are defined as

$$\begin{array}{cc}& b_{ij}=-\frac{\mathrm{1}}{\mathrm{2}}\left(d_{ij}^{\mathrm{2}}-\frac{\mathrm{1}}{N}{\displaystyle \sum }_{r=\mathrm{1}}^N{d}_{rj}^{\mathrm{2}}-\frac{\mathrm{1}}{N}{\displaystyle \sum }_{r=\mathrm{1}}^N{d}_{ir}^{\mathrm{2}}+\frac{\mathrm{1}}{N^{\mathrm{2}}}{\displaystyle \sum }_{r=\mathrm{1}}^N\hspace{0.17em}{\displaystyle \sum }_{s=\mathrm{1}}^N{d}_{rs}^{\mathrm{2}}\right).\end{array}$$

(3)

It has been shown that B is a matrix of the dot products of all vectors $x_i$ , where $d_{ij}$ is calculated using the Euclidean distance [29]. The coefficient $1/N$ appearing in this equation is used to locate the origin of the space where all vectors $x_i$ belong. Since B is positive semidefinite matrix, its eigenvectors always exist.

In the second stage, the matrix X is obtained by calculating matrices U and V which result from performing singular decomposition (SVD) to the matrix B :

$$\begin{array}{cc}& \mathit{B}=\mathit{U}\mathit{V}{\mathit{U}}^T,\\ \mathit{X}={\mathit{V}}^{\mathrm{1}/\mathrm{2}}{\mathit{U}}^{\mathit{T}},\end{array}$$

(4)

where U and V contain all eigenvectors and eigenvalues of the double centered squared distance matrix B . Assume that the eigenvectors are ordered by decreasing eigenvalues and m eigenvectors are retained corresponding to the dimension of the space. The jth coordinate of the vector i is given by

$$\begin{array}{cc}& x_{ij}=\sqrt{{\lambda }_j}{U}_{ij},\end{array}$$

(5)

where ${\lambda }_j$ is the jth eigenvalue of B and $U_{ij}$ is the element of U .

Performing the SVD of a matrix B with size $N\times N$ is time consuming [30, 31] since it costs $O(N^3)$ computations with QR decomposition to estimate all eigenvectors and eigenvalues of B . Fortunately, the location problems commonly deal with two- or three-dimensional spaces. Therefore, the complete SVD is unnecessary. In this work, we utilize the fast-fixed point algorithm to estimate the largest eigenvectors of B . In particular, we express the eigenvector estimation as the optimization problem as

$$\begin{array}{cc}& \underset{w}{\mathrm{minimize}}\hspace{0.17em}\hspace{0.17em}{w}^T\mathit{B}w\end{array}$$

(6)

subject to $‖w‖=1$ . This optimization problem is effectively solved by the fast-fixed point updates as follows:

$$\begin{array}{cc}& w⟵\mathit{B}w.\end{array}$$

(7)

The normalization is applied at each iteration, $w\leftarrow w/‖w‖$ . Consequently, the pth eigenvector is found in a similar manner, while the Gram-Schmidt process is used to make the pth eigenvector uncorrelated with the other eigenvectors:

$$\begin{array}{cc}& w_p⟵\mathit{B}{w}_p,\\ w_p⟵w_p-{\displaystyle \sum }_{j=\mathrm{1}}^{p-\mathrm{1}}\left(w_p^T{w}_j\right)w_j.\end{array}$$

(8)

The eigenvalue corresponding to each eigenvector is simply calculated as

$$\begin{array}{cc}& {\lambda }_p=‖\mathit{B}{w}_p‖.\end{array}$$

(9)

The incremental updates presented in (8) and (9) need several iterations to converge.

2.2. Kalman Filtering

The KL filter [23] is the estimate of hidden states of a system given its measurements (or observations). This filter deals with the appearance of Gaussian noise in both state transition and measurement process. The KL filter is modeled by

$$\begin{array}{cc}& {\theta }_k={\mathit{F}}_k{\theta }_{k-\mathrm{1}}+w_{k-\mathrm{1}},\\ Z_k={\mathit{H}}_k{\theta }_k+{\nu }_k,\end{array}$$

(10)

where $w_{k-\mathrm{1}}$ and ${\nu }_k$ are the Gaussian noises coming from the transition process from state $k-1$ to k and from measurement, ${\theta }_k$ is the state of the system at the time index k, and $Z_k$ is the measurement of ${\mathit{H}}_k{\theta }_k$ . The matrix ${\mathit{F}}_k$ is used to predict the state from the previous one and the matrix ${\mathit{H}}_k$ presents the measurements given the state. Suppose that the process error $w_{k-1}$ and measurement error $v_k$ are zero mean Gaussian noises with covariance matrices of ${\mathit{Q}}_k$ and ${\mathit{R}}_k$ , respectively. The predictions thereby are performed by

$$\begin{array}{cc}& {\theta }_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}={\mathit{F}}_k{\theta }_{k-\mathrm{1}},\\ {\mathit{P}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}={\mathit{F}}_k{\mathit{P}}_{k-\mathrm{1}}{\mathit{F}}_k^T+{\mathit{Q}}_k.\end{array}$$

(11)

Meanwhile, a set of updated equations [23] are given by

$$\begin{array}{cc}& \delta Z_k=Z_k-{\mathit{H}}_k{\theta }_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}},\\ {\mathit{S}}_k={\mathit{H}}_k{\mathit{P}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\mathit{H}}_k^T+{\mathit{R}}_k,\\ {\mathit{K}}_k={\mathit{P}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\mathit{H}}_k^T{\mathit{S}}_k^{-\mathrm{1}},\\ X_k=X_{k\hspace{0.17em}|\hspace{0.17em}k-1}+{\mathit{K}}_k\delta Z_k,\\ {\mathit{P}}_k=\left(\mathit{I}-{\mathit{K}}_k{\mathit{H}}_k\right){\mathit{P}}_{k\hspace{0.17em}|\hspace{0.17em}k-1},\end{array}$$

(12)

where I is the identity matrix, $\delta Z_k$ is the difference between predicted observation and measurement, ${\mathit{P}}_k$ is the covariance matrix of state errors, and ${\mathit{K}}_k$ is the Kalman gain.

2.3. MDS with Kalman Filter for Tracking Location of Nodes

Obviously, the function to locate sensor nodes from distance realized in (5) is nonlinear and is based on eigenvalue decomposition. Therefore, it is complicated to establish a direct model to reflect the pairwise distance errors with the node locations. In this work, we track the changes of the double centered square distance matrix ${\mathit{B}}_k$ from the given measurement ${\mathit{D}}_k$ . The state of the tracking system is the vector

$$\begin{array}{cc}& {\theta }_k=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\left({\mathit{B}}_k\right).\end{array}$$

(13)

Meanwhile, the measurement vector $Z_k$ is the double centered square distance matrix of the real distance measurement ${\tilde{\mathit{D}}}_k$ :

$$\begin{array}{cc}& Z_k=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\left(-\frac{\mathrm{1}}{\mathrm{2}}\mathit{J}{\tilde{\mathit{D}}}_k\mathit{J}\right).\end{array}$$

(14)

It is obvious that the noise model to reflect the measurement error ${\mathit{Z}}_k$ no longer follows Gaussian noise and the KL filter cannot be applied. In order to find the error ${\mathit{\nu }}_k$ of estimating ${\mathit{Z}}_k$ , we start from the measurement error of the distances between two sensor nodes i and j:

$$\begin{array}{cc}& {\tilde{d}}_{ij,k}=d_{ij,k}+{\omega }_{ij,k},\end{array}$$

(15)

where ${\omega }_{ij,k}~\mathcal{N}(0,{\sigma }_{ij,k}^2)$ is independent and identically distributed (i.i.d) noise and ${\tilde{d}}_{ij,k}$ is the actual distance measured by sensors. Let $\mathcal{N}(\mathrm{0},{\sigma }_{ij,k}^2)$ be Gaussian noise with zero mean and variance ${\sigma }_{ij,k}^2$ . A model widely used in the literature assumes that artifacts increase with a longer measurement distance [22, 32]. In this case, ${\sigma }_{ij,k}^2=d_{ij,k}^2/{\gamma }^2$ , where γ is a constant value. In this work, errors to estimate the squared distance matrix ${\mathit{D}}_k$ , ${\left[\mathit{D}\right]}_{ij,k}=d_{ij,k}^2$ are expressed by a vector ${\xi }_k=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\mathbf{\left(}{\mathcal{E}}_k\mathbf{\right)}$ , where

$$\begin{array}{cc}& {\left[\mathcal{E}\right]}_{ij,k}=r_{ij,k}^{\mathrm{2}}-d_{ij,k}^{\mathrm{2}}=\mathrm{2}{d}_{ij,k}{\omega }_{ij,k}+{\omega }_{ij,k}^{\mathrm{2}}\approx \mathrm{2}{d}_{ij,k}{\omega }_{ij,k}.\end{array}$$

(16)

Since the value of $d_{ij,k}$ is unknown, the estimation of   ${\xi }_k$ is nontrivial. We use ${\tilde{d}}_{ij,k}$ instead to further approximate ${\mathbf{\left[}\mathcal{E}\mathbf{\right]}}_{ij,k}$ as ${\mathbf{\left[}\mathcal{E}\mathbf{\right]}}_{ij,k}\approx 2{\tilde{d}}_{ij,k}{\omega }_{ij,k}$ . Therefore, the covariance matrix of this measurement error is defined as

$$\begin{array}{cc}& {\mathbb{R}}_k=\mathbb{E}\left\{{\xi }_k{\xi }_k^T\right\}=\mathrm{4}{\mathit{\Psi }}_k{\mathrm{\Xi }}_k{\mathit{\Psi }}_k,\\ {\mathit{\Psi }}_k=\mathrm{diag}\left({\tilde{d}}_{\mathrm{11},k},{\tilde{d}}_{\mathrm{12},k},\dots ,{\tilde{d}}_{NN,k}\right),\\ {\mathrm{\Xi }}_k=\mathrm{diag}\left({\sigma }_{\mathrm{11},k}^{\mathrm{2}},{\sigma }_{\mathrm{12},k}^{\mathrm{2}},\dots ,{\sigma }_{NN,k}^{\mathrm{2}}\right).\end{array}$$

(17)

Based on the definition of ${\mathit{B}}_k=-\left(1/2\right){\mathit{C}}_N{\tilde{\mathit{D}}}_k{\mathit{C}}_N$ , we have

$$\begin{array}{cc}& Z_k=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\left({\mathit{B}}_k\right)=-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{C}}_N\otimes {\mathit{C}}_N\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\left({\tilde{\mathit{D}}}_k\right),\end{array}$$

(18)

where ⊗ is the Kronecker product. Let $\mathit{\Lambda }={\mathit{C}}_N\otimes {\mathit{C}}_N$ ; the covariance matrix of estimating $Z_k$ is given by

$$\begin{array}{cc}& \mathbb{E}\left\{\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\left({\mathit{B}}_k\right)\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}{\left({\mathit{B}}_k\right)}^T\right\}={\mathit{\Lambda }\mathit{\Psi }}_k{\mathrm{\Xi }}_k{\mathit{\Psi }}_k\mathit{\Lambda }.\end{array}$$

(19)

Now, we establish the state equation and the measurement error for tracking sensor position based on MDS. Let the definition of the state of the system be a vector with size $N^2\times 1$ , ${\theta }_k=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}({\mathit{B}}_k)$ . The corresponding state equation can be presented as

$$\begin{array}{cc}& {\theta }_k=\mathit{F}{\theta }_{k-\mathrm{1}}+w_k,\end{array}$$

(20)

where $\mathit{F}={\mathit{I}}_{N^2}$ . A variety of the model can be concerned to describe the transition process from ${\theta }_{k-1}$ to ${\theta }_k$ . Nevertheless, we employ a low motion estimate in which a moving step is modeled by a random walk process with variance ${\mathit{Q}}_k={{\sigma }_w^2\mathit{I}}_{N^2}$ . Meanwhile, $Z_k$ holds the measurements of the pairwise distance matrix after being doubled centered:

$$\begin{array}{cc}& Z_k=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\left({\tilde{\mathit{B}}}_k\right)=\mathit{H}{\theta }_k+{\nu }_k,\end{array}$$

(21)

where H is simply defined as ${\mathit{I}}_{N^2}$ . Finally, the covariance matrix of ${\nu }_k$ is computed by

$$\begin{array}{cc}& {\mathbf{R}}_k=\mathbb{E}\left\{{\nu }_k{\nu }_k^T\right\}={\mathit{\Lambda }\mathit{\Psi }}_k{\mathrm{\Xi }}_k{\mathit{\Psi }}_k\mathit{\Lambda }.\end{array}$$

(22)

After all matrices ${\mathit{F}}_k$ , ${\mathit{H}}_k$ , ${\mathit{Q}}_k$ , and ${\mathit{R}}_k$ are well defined, prediction and update are performed normally as described in (11) and (12), respectively. Since the state of the system described by ${\mathit{B}}_k$ is tracked on each frame, the corresponding eigenvector and eigenvalue are computed and the node location is found. Note that the nodes from ${\mathit{B}}_k$ are centered around the origin since ${\mathit{B}}_k$ is doubly centered. Therefore, some beacons with known locations are used to perform registration of rigid transformations including translation, rotation, and reflection to move a point cloud to the right position.

Let ${\mathit{X}}_A$ and ${\mathit{X}}_B$ be the coordinates of beacons of two consecutive snapshots. The registration of rotation/reflection matrix R is computed by minimizing

$$\begin{array}{cc}& \mathbf{R}=\underset{R}{\mathrm{minimize}}{‖{\mathit{X}}_A-\mathit{R}{\mathit{X}}_B‖}^{\mathrm{2}}.\end{array}$$

(23)

The solution for this problem is acquired by estimating the SVD of matrix ${\mathit{X}}_A{\mathit{X}}_B^T={\mathit{U}}_{AB}{\mathit{\Lambda }}_{AB}{\mathit{V}}_{AB}^T$ and

$$\begin{array}{cc}& \mathit{R}={\mathit{U}}_{AB}{\mathit{V}}_{AB}^T.\end{array}$$

(24)

Let ${\mu }_A$ and ${\mu }_B$ be the centers of the two sets A and B, respectively. With regard to the above transformation, the whole set of nodes ${\mathit{Y}}_B$ is registered to the whole set ${\mathit{Y}}_A$ by

$$\begin{array}{cc}& {\mathit{Y}}_B=\mathit{R}\left({\mathit{Y}}_B-{\mu }_A\right)+{\mu }_B.\end{array}$$

(25)

2.4. Tracking with Extended Kalman Filtering

In EKL [23], the matrices ${\mathit{F}}_k$ and ${\mathit{H}}_k$ in (10) are described by nonlinear functions:

$$\begin{array}{cc}& {\theta }_k=f\left({\theta }_{k-\mathrm{1}}\right)+w_{k-\mathrm{1}},\\ Z_k=h\left({\theta }_k\right)+{\nu }_k.\end{array}$$

(26)

The prediction and update steps are performed normally in which the Jacobian matrices of f and g are used to present ${\mathit{F}}_k$ and ${\mathit{H}}_k$ , respectively:

$$\begin{array}{cc}& {\mathit{F}}_k={\left.\frac{\partial f}{\partial \theta }\right|}_{{\theta }_{k-\mathrm{1}}},\\ {\mathit{H}}_k={\left.\frac{\partial h}{\partial \theta }\right|}_{{\theta }_k}.\end{array}$$

(27)

Besides, there are some modifications on (11) and (12) given as follows:

$$\begin{array}{cc}& {\theta }_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}=f\left({\theta }_{k-\mathrm{1}}\right)\\ \delta Z_k=Z_k-h\left({\theta }_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right).\end{array}$$

(28)

In the work presented in [27], the EKL filter is applied for tracking nodes with the pairwise measurements. Similarly, we establish the testing environment where the distance measurements are available. Therefore, the measurement vector $Z_k$ , a vector with a size $N^2\times 1$ , is expressed by

$$\begin{array}{cc}& Z_k=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\left(\left[\begin{array}{ccc}{\tilde{d}}_{\mathrm{11},k}& \cdots & {\tilde{d}}_{\mathrm{1}N,k}\\ ⋮& \ddots & ⋮\\ {\tilde{d}}_{N\mathrm{1},k}& \cdots & {\tilde{d}}_{NN,k}\end{array}\right]\right)={\left[\begin{array}{cccc}{d}_{\mathrm{11},k}& d_{\mathrm{12},k}& \cdots & d_{NN,k}\end{array}\right]}^T+{\nu }_k.\end{array}$$

(29)

As aforementioned, ${\tilde{d}}_{ij,k}$ is the actual range measured by sensors. The same additive noises as presented by (15) are used to model the artifacts on the measurements. Thus the noise vector ${\nu }_k$ is formed by concatenating the i.d.d noise sources ${\omega }_{ij,k}~\mathcal{N}(0,{\sigma }_{ij,k}^2)$ , ${\sigma }_{ij,k}^2\approx {\tilde{d}}_{ij,k}^2/{\gamma }^2$ , and

$$\begin{array}{cc}& v_k=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\left(\left[\begin{array}{ccc}{\omega }_{\mathrm{11},k}& \cdots & {\omega }_{\mathrm{1}N,k}\\ ⋮& \ddots & ⋮\\ {\omega }_{N\mathrm{1},k}& \cdots & {\omega }_{NN,k}\end{array}\right]\right)={\left[\begin{array}{cccc}{\omega }_{\mathrm{11},k}& {\omega }_{\mathrm{12},k}& \cdots & {\omega }_{NN,k}\end{array}\right]}^T.\end{array}$$

(30)

The corresponding covariance matrix ${\mathit{R}}_k$ of ${\nu }_k$ is modeled by

$$\begin{array}{cc}& {\mathit{R}}_k={\left[\begin{array}{ccc}\frac{{\tilde{d}}_{\mathrm{11},k}^{\mathrm{2}}}{{\gamma }^{\mathrm{2}}}& \cdots & \mathrm{0}\\ ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & \frac{{\tilde{d}}_{NN,k}^{\mathrm{2}}}{{\gamma }^{\mathrm{2}}}\end{array}\right]}_{{N^{\mathrm{2}}\times N}^{\mathrm{2}}}.\end{array}$$

(31)

The state of a system ${\theta }_k$ is a vector of a size $mN\times 1$ created by $\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}({\mathit{\varphi }}_k)$ , where

$$\begin{array}{cc}& {\mathit{\varphi }}_k=\left[x_{\mathrm{1},k},x_{\mathrm{2},k},\dots ,x_{N,k}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{11},k}& x_{\mathrm{21},k}& \cdots & x_{N\mathrm{1},k}\\ x_{\mathrm{12},k}& x_{\mathrm{22},k}& \cdots & x_{N\mathrm{2},k}\\ ⋮& ⋮& ⋮& ⋮\\ x_{\mathrm{1}m,k}& x_{\mathrm{2}m,k}& \cdots & x_{Nm,k}\end{array}\right].\end{array}$$

(32)

It is noted that $x_{i,k}={[x_{i1,k},x_{i2,k},\dots ,x_{im,k}]}^T$ is the coordinate of the point i in the snapshot k. The transition from ${\theta }_{k-\mathrm{1}}$ to ${\theta }_k$ is simply presented by the random walk process with the variance ${\mathit{Q}}_k={{\sigma }_w^{\mathrm{2}}\mathit{I}}_{N^{\mathrm{2}}}$ , ${\mathit{F}}_k={\mathit{I}}_{m\times N}$ , and $f\left({\theta }_{k-1}\right)={\mathit{F}}_k{\theta }_{k-1}$ . However, the function $h({\theta }_k)$ is nonlinear:

$$\begin{array}{cc}& h\left({\theta }_k\right)={\left[\begin{array}{cccc}{d}_{\mathrm{11},k}& d_{\mathrm{12},k}& \cdots & d_{NN,k}\end{array}\right]}^T\end{array}$$

(33)

in which $d_{ij,k}=‖x_{i,k}-x_{j,k}‖=\sqrt{{\sum }_{l=\mathrm{1}}^m{(x_{il,k}-x_{jl,k})}^{\mathrm{2}}}$ . Thus, the Jacobian matrix ${\mathit{H}}_k={(\left.\partial h/\partial \theta )\right|}_{{\theta }_k}$ is a matrix of partial derivatives $\partial d_{ij,k}/\partial x_{il,k}$ and $\partial d_{ij,k}/\partial x_{jl,k}$ , $j=\mathrm{1,2},\dots ,m$ .

From the defined ${\mathit{F}}_k$ , ${\mathit{H}}_k$ , ${\mathit{Q}}_k$ , and ${\mathit{R}}_k$ , the predict and update steps of EKL are performed normally as KL filter. However, as ${\mathit{H}}_k$ is the Jacobian matrix of the nonlinear function in high dimensional space, the update of EKL may approach the suboptimal solution. Thus, the number of nonlinear constraints is usually relaxed by applying EKL on separate groups of nodes and the tracking results are combined for the whole network.

3. Experimental Results

This section shows the performance of the proposed tracking algorithm when compared with conventional approaches of localization without tracking. First of all, a fully connected network of $N=20$ sensor nodes in 2D is used for evaluation. The nodes are initially distributed in an area of $100\hspace{0.17em}\mathrm{m}\times 100\hspace{0.17em}\mathrm{m}$ and four beacons at the coordinates (20, 20) m, (20, 80) m, (80, 20) m, and (80, 80) m are fixed for the rigid registration of nodes. At each time snapshot, the nodes move away from the beacons in both vertical and horizontal directions. The moving step is uniformly assigned in the range (0, 2) m. With regard to noise increases proportional to the distance between nodes, the accuracy of location estimation decays overtime. To compare the performance of alternative approaches, the mean error distance (MED) at the k-snapshot is calculated as

$$\begin{array}{cc}& {\mathrm{M}\mathrm{E}\mathrm{D}}_k=\frac{{\sum }_{i=\mathrm{1}}^N{d}_{i,k}}{N},\end{array}$$

(34)

where $d_{i,k}$ is distance difference between the real ( ${\tilde{x}}_{i,k},{\tilde{y}}_{i,k}$ ) and the estimated location $(x_{i,k},y_{i,k})$ of the nodes $d_{i,k}=\sqrt{{(x_{i,k}-{\tilde{x}}_{i,k})}^{\mathrm{2}}+{(y_{i,k}-{\tilde{y}}_{i,k})}^{\mathrm{2}}}$ .

As aforementioned, the distance measurements are affected by additive noise due to the power loss of the signal during transmission. The amplitude or variance of the noise is proportional to the distance between nodes by the equation ${\sigma }_{ij}^2=d_{ij}^2/{\gamma }^2$ . For the first evaluation, we compare MED over the snapshot for different value of γ, in particular, $\mathrm{10}\hspace{0.17em}{\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{10}\hspace{0.17em}\mathrm{d}\mathrm{B}$ , $1\mathrm{0}\hspace{0.17em}{\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{20}\hspace{0.17em}\mathrm{d}\mathrm{B}$ , and $\mathrm{10}\hspace{0.17em}{\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{30}\hspace{0.17em}\mathrm{d}\mathrm{B}$ , respectively, and show the compared results in Figure 1. Once the nodes move far away from the origin, the noises raise more effects on the range measurement. Thus the localization outcomes acquired by MDS are less accurate. Meanwhile, the errors specified by MED of our MDS integrated with KL filter are attenuated and are lower than that of conventional MDS. Besides, MED returned by conventional MDS may be as large as several hundreds of meters. The MED estimated by our approach is approximately less than 5 m with $\mathrm{10}\hspace{0.17em}{\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{30}\hspace{0.17em}\mathrm{d}\mathrm{B}$ , less than 20 m with $\mathrm{10}\hspace{0.17em}{\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{20}\hspace{0.17em}\mathrm{d}\mathrm{B}$ , and less than 50 m with $\mathrm{10}\hspace{0.17em}{\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{10}\hspace{0.17em}\mathrm{d}\mathrm{B}$ . This shows the efficiency of our proposed approach to enhance the performance of location tracking.

Figure 1

MED versus snapshot for (a) $10{\hspace{0.17em}\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{10}\hspace{0.17em}\mathrm{d}\mathrm{B}$ , (b) $10\hspace{0.17em}{\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{20}\hspace{0.17em}\mathrm{d}\mathrm{B}$ , and (c) $10{\hspace{0.17em}\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{30}\hspace{0.17em}\mathrm{d}\mathrm{B}$ .

The qualitative evaluation of the MDS combined with KL filter is conducted by showing the positions between the original nodes and their estimations. An example of distribution of sensor nodes at the snapshot 200 and $\gamma =\mathrm{50}$ is depicted in Figure 2. Due to noises, the estimated positions are shifted far from the true coordinates of nodes. From another aspect, an example of trajectories of a node at the snapshot 300 found by MDS and our method are compared with the ground truth as described in Figure 3 ( $\mathrm{10}\hspace{0.17em}{\mathrm{l}\mathrm{o}\mathrm{g}}_{\mathrm{10}}\left(\gamma \right)=\mathrm{20}\hspace{0.17em}\mathrm{d}\mathrm{B}$ ). Obviously, the moving trajectory of a node found by MDS presents rapid changes even between two consecutive snapshots. The node trajectory found by MDS filtered with KL is smoother and more fitted to the real movement of nodes.

Figure 2

True and estimated position of sensor nodes.

Figure 3

Estimated trajectories of sensor nodes.

Finally, we evaluate the efficiency of the proposed approach against conventional MDS [5–7] described in Section 2.1 and EKL [27] in Section 2.3 for tracking nodes based on the range measurements. Two fully connected networks of $N=10$ and of $N=20$ sensors in 2D are examined. To estimate the overall accuracy of locating the position of nodes over 500 snapshots, we compute the average MED. The overall comparison is shown in Figure 4 where different values of γ are validated. Although EKL mitigates the artifacts of range measurements, the proposed MDS with tracking is always superior to both conventional MDS and EKL to address localization problems.

Figure 4

Averaged MED over 500 snapshots versus γ.

4. Conclusions and Future Work

In this paper, an effective tracking algorithm has been proposed for tracking the position of sensors by measuring the distances between them. Basically using three beacons with known locations, a sensor is approximately located. In the case a pairwise distance matrix of all sensor nodes is provided, MDS is utilized. However, due to the energy loss of transmitting signals from a source to a detector, the distance measurement is affected with noise. In particular, the noise level increases when a sensor node is far away from the others. This correspondingly reduces the accuracy of the MDS approach. Therefore, we have combined conventional MDS with the KL filter for location tracking. Our proposed approach is validated with simulation experiments and we showed that our algorithm is superior to the conventional MDS and EKL tracking on the localization accuracy. The algorithm is effective and significantly reduces the noise effects from a real environment. Therefore, the proposed method has potential to address mobile sensor localization with the practical requirement of high accuracy localization.

In future work, we target applying our algorithm to partially connected networks where the distance measurement is feasible for only close-by nodes. Besides, the KL filter can be extended into the Bernoulli filter [33] to address the localization problem of switching the states (on/off) of sensors.