Indoor Robot Localization Based on Multidimensional Scaling

Abstract

In pertinence to the application of multidimensional scaling (MDS) methods in ranging-based positioning systems, an analysis is firstly conducted by the classical MDS algorithm. Modified MDS algorithm and subspace method are presented in localization application. We also depicted the unified framework and general solutions of MDS methods. However, the least square solutions under this framework are not optimal. Their performance is still related to selection of coordinate reference points. To address this problem, a minimum residual MDS algorithm based on particle swarm optimization (PSO) is proposed to derive a new solution for indoor robot localization under the unified framework. The result of analysis indicates that the performance of minimum residual MDS method is immune to selection of reference points. Furthermore, the localization accuracy for indoor robot has been enhanced by 41% as compared with the classical MDS algorithm.

1. Introduction

Location information is very significant in such scenarios as military affairs, industry, civil use, disaster relief, and so forth [1]. However, the widely used GPS system is hardly applicable to indoor environments. Wireless sensor network (WSN) has provided a new solution for accurate localization and tracking of indoor targets, but environmental interference and influence are prone to cause larger error of localization and tracking. Therefore, high-efficient algorithms are necessary to enhance localization accuracy of WSN.

Existing algorithms in wireless sensor network localization can be broadly divided into two classes: range-based and range-free. The class of range-based algorithms is required to provide the accurate distance estimation between the target node and anchor node. Moreover, range-based algorithms require each sensor node to be equipped with more powerful CPU. There are several techniques to get the 4 distance estimation: received signal strength indicator (RSSI) [2], time of arrival (TOA), and time different of arrival (TDOA) [3]. The class of range-free algorithms can reduce the energy consumption and the demands for special hardware, but their accuracy is lower than the former one, such as APIT [4] and DV-Hop algorithms. In order to improve localization accuracy, the swarm intelligence algorithms [5] are applied by the researcher to node localization. However, these algorithms need a large number of anchor nodes. With regard to this, the literature [6] proposed the distributed localization algorithm. However, how to reduce localization errors has not been considered by the distributed localization algorithm. The distributed weighted combining scheme for cooperative spectrum sensing in cognitive radio networks was developed [7, 8]. This technique can incorporate specific weights according to different channel conditions and exhibits clear advantages with respect to channel fading, low PU transmission power, low false alarm rate, and the network size variation.

In recent years, a lot of novel ideas and solutions have emerged for WSN localization. The localization technique based on multidimensional scaling (MDS), first proposed by Shang et al. [9], has offered a new idea for node localization. MDS method was widely applied in data analysis and data visualization in fields of physics, biology, geology, brain science, and social and psychological phenomenon research, as well as in pattern recognition, machine learning, and data mining, which are concerned with dimension reduction and feature extraction, and so forth [10–12]. Compared with previous localization algorithms, the MDS-based localization algorithms can simultaneously locate multiple nodes by utilizing the associated information among all nodes within the network. It can even get a diagram of relative positions among nodes while anchor nodes are unavailable. When absolute location is in need, fewer anchor nodes are required for this algorithm. Anchor nodes are also deployed without restriction. Recently, a variety of localization methods derived from classical MDS method have been applied in sensor network [13, 14] as well as cellular network [15]. The classical MDS method takes the mass centers of all nodes (including all anchor nodes and target nodes) as coordinate reference points. Thus, Cheung and So think it is inapplicable to wireless localization without prior information about the location of target nodes. Therefore, they propose a modified MDS method with mass centers of all anchor nodes as coordinate reference points [16]. At the same time, Wan et al. provide a subspace solution used for three-station positioning systems, with the locations of target nodes as coordinate reference points [17]. On this basis, So and Chan extend the subspace method to multistation positioning systems [18]. Chen et al. provide a unified framework of MDS method and its general solution for range-based localization, among which the above methods are several particular cases. In [19], a weighted MDS method under the unified framework has been obtained. However, according to the Gaussian Markov Theorem, all methods under this framework are not optimal. Their performance still depends on selection of coordinate reference points. Furthermore, although the weighted MDS method does not depend on selection of reference points, known environmental parameters, which are very difficult to obtain, are needed in the solving process.

Against the flaws and defects of the above algorithms, this paper proposes a minimum residual MDS (MinRE MDS) algorithm under the unified framework. We also provide a unified solution with minimum residual to MDS algorithms, including the classical, the modified, and the subspace-based ones. With this framework, the result of analysis indicates that the performance of the minimum residual MDS method does not depend on selection of reference points. Therefore, it does not require environmental data and can enhance the localization accuracy greatly as compared with the classical MDS algorithm.

2. Application of MDS Algorithm in Localization Systems

Given the dimension of the target space, relative positions of all objects in the space can be determined by the classical MDS method. In order to determine the absolute position of target nodes, it is required to perform rotation and translation over the matrix of relative coordinates. In this paper, we conduct theoretical analysis of three important MDS-based localization methods, namely, the classical MDS algorithm, the modified MDS algorithm, and the subspace method.

2.1. Classical MDS

For range-based localization system, the target nodes (TN) whose positions are to be determined are located at $x_{}\in R^n$ . N anchor nodes (AN) with known positions are located at $x_i\in R^n,\hspace{0.17em}\hspace{0.17em}i=\mathrm{1},\dots ,N\hspace{0.17em}\hspace{0.17em}(N\ge \mathrm{3})$ . Thus, a centralized coordinate matrix composed of the coordinates of TN and AN can be constructed:

$$\begin{array}{cc}& X_c={\left[\begin{array}{cccc}x-x_c& x_{\mathrm{1}}-x_c& \cdots & x_N-x_c\end{array}\right]}^T,\end{array}$$

(1)

where

$$\begin{array}{cc}& x_c=\frac{\mathrm{1}}{N+\mathrm{1}}\left(x+{\displaystyle \sum }_{i=\mathrm{1}}^N{x}_i\right)\end{array}$$

(2)

is mass-center coordinate. Real distance between the target node and the ith anchor node can be expressed as

$$\begin{array}{cc}& d_i={‖x-x_i‖}_{\mathrm{2}}\hspace{1em}i=\mathrm{1},\dots N,\end{array}$$

(3)

while the distance between anchor nodes can be expressed as

$$\begin{array}{cc}& d_{ij}={‖x_i-x_j‖}_{\mathrm{2}}\hspace{1em}i,j=\mathrm{1},\dots N.\end{array}$$

(4)

Due to the existence of measurement noise, the distance between the target node and the ith anchor node is usually expressed as

$$\begin{array}{cc}& r_i=d_i+{\eta }_i,\end{array}$$

(5)

where ${\eta }_i$ is the measurement noise as is usually supposed as Gaussian white noise under line-of-sight (LOS) settings.

When the classical MDS method is applied for localization in this system, detailed procedures are as follows.

Since the space encompassing the target is known, the dimension n can be determined. $n=2$ corresponds to 2D localization, whereas $n=3$ corresponds to 3D localization. In this paper, we only focus on the application of MDS methods in 2D localization. When the information of distance measure is obtained, the matrix of distances can be constructed as follows:

$$\begin{array}{cc}& \tilde{D}=\left[\begin{array}{ccccc}\mathrm{0}& r_{\mathrm{1}}^{\mathrm{2}}& r_{\mathrm{2}}^{\mathrm{2}}& \cdots & r_N^{\mathrm{2}}\\ r_{\mathrm{1}}^{\mathrm{2}}& 0& d_{\mathrm{12}}^{\mathrm{2}}& \cdots & d_{\mathrm{1}N}^{\mathrm{2}}\\ r_{\mathrm{2}}^{\mathrm{2}}& d_{\mathrm{21}}^{\mathrm{2}}& 0& \cdots & d_{\mathrm{1}N}^{\mathrm{2}}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ r_N^{\mathrm{2}}& d_{N\mathrm{1}}^{\mathrm{2}}& d_{N\mathrm{2}}^{\mathrm{2}}& \cdots & \mathrm{0}\end{array}\right].\end{array}$$

(6)

Apply double central transformation to get

$$\begin{array}{cc}& {\tilde{B}}_c=-\frac{\mathrm{1}}{\mathrm{2}}J\tilde{D}{J}^T,\end{array}$$

(7)

where

$$\begin{array}{cc}& J=I_{N+\mathrm{1}}-\frac{\mathrm{1}}{N+\mathrm{1}}{\mathrm{1}}_{N+\mathrm{1}}{\mathrm{1}}_{N+\mathrm{1}}^T,\end{array}$$

(8)

and $I_{N+\mathrm{1}}$ is a $\left(N+\mathrm{1}\right)\times \left(N+\mathrm{1}\right)$ identity matrix. Then, eigen decomposition is applied to ${\tilde{B}}_c$ . Since $n=\mathrm{2}$ , let

$$\begin{array}{cc}& {\tilde{X}}_c={\tilde{V}}_{c\mathrm{1}}{\tilde{\Lambda }}_{c\mathrm{1}}^{\mathrm{1}/\mathrm{2}}={\left[\begin{array}{cccc}\tilde{x}& {\tilde{x}}_{\mathrm{1}}& \cdots & {\tilde{x}}_N\end{array}\right]}^T,\end{array}$$

(9)

where

$$\begin{array}{cc}& {\tilde{\Lambda }}_{c\mathrm{1}}^{\mathrm{1}/\mathrm{2}}=\mathrm{diag}\left\{\sqrt{{\tilde{\lambda }}_{\mathrm{1}}},\sqrt{{\tilde{\lambda }}_{\mathrm{2}}}\right\},\\ {\tilde{V}}_{c\mathrm{1}}=\left[\begin{array}{ccc}{\tilde{v}}_{\mathrm{1}}& \cdots & {\tilde{v}}_{\mathrm{2}}\end{array}\right].\end{array}$$

(10)

${\tilde{\lambda }}_{\mathrm{1}}\ge {\tilde{\lambda }}_{\mathrm{2}}>\mathrm{0}$ are two greatest positive eigenvalues of ${\tilde{B}}_c$ ; ${\tilde{v}}_{\mathrm{1}},{\tilde{v}}_{\mathrm{2}}$ are orthogonal eigenvectors in correspondence to ${\tilde{\lambda }}_{\mathrm{1}},{\tilde{\lambda }}_{\mathrm{2}}$ . ${\tilde{V}}_{c\mathrm{1}}$ is actually a signal subspace of ${\tilde{B}}_c$ .

The coordinates obtained from Formula (9) are merely relative position of the nodes within the space. In order to obtain the absolute position of target node, transformation is required by making further use of the position information of anchor nodes. Without measurement error, there is such a relation between ${\tilde{X}}_c$ and the centralized coordinate matrix $X_c$ [20]:

$$\begin{array}{cc}& X_c={\tilde{X}}_c\Omega ,\end{array}$$

(11)

where

$$\begin{array}{cc}& \Omega ={\left({\tilde{X}}_c^T{X}_c{X}_c^T{\tilde{X}}_c\right)}^{\mathrm{1}/\mathrm{2}}{\left(X_c^T{\tilde{X}}_c\right)}^{-\mathrm{1}}\end{array}$$

(12)

is an orthogonal matrix, which has undergone rotation transformation with $X_c$ as the center. Through the rotation transformation and translation, the position estimation of target node can be obtained as

$$\begin{array}{cc}& \widehat{x}={\Omega }^T\tilde{x}+x_c.\end{array}$$

(13)

Formula (12) reveals that the transformation matrix Ω is dependent on the centralized coordinate matrix $X_c$ , which contains the unknown target node. Meanwhile, the center-of-mass coordinate $x_c$ in (13) also contains the coordinate matrix of unknown target node. Therefore, it is bound to cause a larger error when applying the classical MDS method to estimate the target node's position. For this problem, a modified MDS method applicable to mobile station localization was presented in literature [16].

2.2. Modified MDS

The modified MDS algorithm selects the mass centers of all anchor nodes as coordinate reference points, with respect to which the coordinate matrix composed of the target node and anchor nodes can be expressed as

$$\begin{array}{cc}& X_b={\left[\begin{array}{cccc}x-x_b& x_{\mathrm{1}}-x_b& \cdots & x_N-x_b\end{array}\right]}^T,\end{array}$$

(14)

where

$$\begin{array}{cc}& x_b=\frac{\mathrm{1}}{N}{\displaystyle \sum }_{i=\mathrm{1}}^N{x}_i.\end{array}$$

(15)

By separating the expression of the component in Formula (14) which is related to anchor nodes, we can get

$$\begin{array}{cc}& X_{lb}={\left[\begin{array}{ccc}{x}_{\mathrm{1}}-x_b& \cdots & x_N-x_b\end{array}\right]}^T.\end{array}$$

(16)

Divide ${\tilde{B}}_c$ in Formula (7) into

$$\begin{array}{cc}& {\tilde{B}}_c=\left[\begin{array}{cc}{B}_{\mathrm{11}}& b^T\\ b& {\tilde{B}}_{\mathrm{1}}\end{array}\right],\end{array}$$

(17)

where $B_{\mathrm{11}}$ is the first diagonal element of matrix ${\tilde{B}}_c$ , b is an $N\times \mathrm{1}$ vector, and ${\tilde{B}}_{\mathrm{1}}$ is an $N\times N$ matrix. Then, the following matrix can be constructed by Formulas (16) and (17):

$$\begin{array}{cc}& {\tilde{B}}_b=\left[\begin{array}{cc}{\left(\frac{N+\mathrm{1}}{N}\right)}^{\mathrm{2}}{B}_{\mathrm{11}}& \left(\frac{N+\mathrm{1}}{N}\right)b^T+\left(\frac{N+\mathrm{1}}{N^{\mathrm{2}}}\right)B_{\mathrm{11}}{\mathrm{1}}_N^T\\ \left(\frac{N+\mathrm{1}}{N}\right)b+\left(\frac{N+\mathrm{1}}{N^{\mathrm{2}}}\right)B_{\mathrm{11}}{\mathrm{1}}_N& X_{lb}{X}_{lb}^T\end{array}\right].\end{array}$$

(18)

It has been shown in [16] that without measuring error

$$\begin{array}{cc}& {\tilde{B}}_b=X_b{X}_b^T.\end{array}$$

(19)

In 2D localization, eigen decomposition is firstly applied to matrix ${\tilde{B}}_b$ , which is constructed by Formula (18). Then, the coordinate matrix estimation ${\tilde{X}}_b$ can be obtained according to Formula (9) and divided into

$$\begin{array}{cc}& {\tilde{X}}_b={\left[\begin{array}{cc}\tilde{x}& {\tilde{X}}_{lb}\end{array}\right]}^T.\end{array}$$

(20)

Thus, the transformation matrix with $x_b$ as the center can be expressed as

$$\begin{array}{cc}& {\Omega }_b={\left({\tilde{X}}_{lb}^T{X}_{lb}{X}_{lb}^T{\tilde{X}}_{lb}\right)}^{\mathrm{1}/\mathrm{2}}{\left(X_{lb}^T{\tilde{X}}_{lb}\right)}^{-\mathrm{1}}.\end{array}$$

(21)

Therefore, the estimated coordinate of the target node can be obtained as

$$\begin{array}{cc}& \widehat{x}={\Omega }_b^T\tilde{x}+x_b.\end{array}$$

(22)

2.3. Subspace-Based MDS

Unlike the modified MDS method in which the mass centers of all anchor nodes are used as the coordinate reference points, the subspace-based MDS localization algorithm takes the coordinate of target node as the coordinate reference point. The corresponding coordinate matrix can be expressed as

$$\begin{array}{cc}& X_s={\left[x_{\mathrm{1}}-x\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{x}_{\mathrm{2}}-x\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\cdots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{x}_N-x\right]}^T.\end{array}$$

(23)

Construct matrix ${\tilde{B}}_s$ , in which the element at the ith row and jth column is

$$\begin{array}{cc}& {\left[{\tilde{B}}_s\right]}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left(r_i^{\mathrm{2}}+r_j^{\mathrm{2}}-d_{ij}^{\mathrm{2}}\right).\end{array}$$

(24)

Apply eigen decomposition to ${\tilde{B}}_s$ ; we can get

$$\begin{array}{cc}& {\tilde{X}}_s={\tilde{V}}_{s\mathrm{1}}{\tilde{\Lambda }}_{s\mathrm{1}}^{\mathrm{1}/\mathrm{2}}={\left[{\tilde{x}}_{\mathrm{1}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\tilde{x}}_{\mathrm{2}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\cdots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\tilde{x}}_N\right]}^T\end{array}$$

(25)

in 2D localization in the same way, where ${\tilde{\Lambda }}_{s\mathrm{1}}$ is the diagonal matrix composed of the two greatest positive eigenvalues of ${\tilde{B}}_s$ and ${\tilde{V}}_{s\mathrm{1}}$ is the matrix composed of the eigenvectors in correspondence to both eigenvalues. Equally, a rotation transformation with the target node as the center exists between $X_s$ and ${\tilde{X}}_s$ ; namely,

$$\begin{array}{cc}& X_s\approx {\tilde{X}}_s{\Omega }_s.\end{array}$$

(26)

${\Omega }_s$ is an unknown variable. Thus, the least square solution of (26) is

$$\begin{array}{cc}& {\Omega }_s={\left({\tilde{X}}_s^T{\tilde{X}}_s\right)}^{-\mathrm{1}}{\tilde{X}}_s^T{X}_s.\end{array}$$

(27)

Substitute Formula (25) into the above expression; we can get

$$\begin{array}{cc}& {\Omega }_s={\tilde{\Lambda }}_{s\mathrm{1}}^{\mathrm{1}/\mathrm{2}}{\tilde{V}}_{s\mathrm{1}}^T{X}_s.\end{array}$$

(28)

With Formulas (26) and (28),

$$\begin{array}{cc}& X_s\approx {\tilde{V}}_{s\mathrm{1}}{\tilde{V}}_{s\mathrm{1}}^T{X}_s\end{array}$$

(29)

since

$$\begin{array}{cc}& I_N-{\tilde{V}}_{s\mathrm{1}}{\tilde{V}}_{s\mathrm{1}}^T={\tilde{V}}_{s\mathrm{2}}{\tilde{V}}_{s\mathrm{2}}^T,\end{array}$$

(30)

where

$$\begin{array}{cc}& {\tilde{V}}_{s\mathrm{2}}=\left[\begin{array}{ccc}{\tilde{v}}_{\mathrm{3}}& \cdots & {\tilde{v}}_N\end{array}\right]\end{array}$$

(31)

is the noise subspace of ${\tilde{B}}_s$ and $I_N$ is an $N\times N$ identity matrix. Substitute Formula (23) into (29) to get

$$\begin{array}{cc}& {\tilde{V}}_{s\mathrm{2}}{\tilde{V}}_{s\mathrm{2}}^T{\mathrm{1}}_N{x}^T\approx {\tilde{V}}_{s\mathrm{2}}{\tilde{V}}_{s\mathrm{2}}^T{\left[\begin{array}{ccc}{x}_{\mathrm{1}}& \cdots & x_N\end{array}\right]}^T.\end{array}$$

(32)

Thus, the least square solution to this equation can pinpoint the estimated position of the target node as

$$\begin{array}{cc}& \widehat{x}=\left[\begin{array}{ccc}{x}_{\mathrm{1}}& \cdots & x_N\end{array}\right]\frac{{\tilde{V}}_{s\mathrm{2}}{\tilde{V}}_{s\mathrm{2}}^T{\mathrm{1}}_N}{{\mathrm{1}}_N^T{\tilde{V}}_{s\mathrm{2}}{\tilde{V}}_{s\mathrm{2}}^T{\mathrm{1}}_N}.\end{array}$$

(33)

3. MDS under Unified Framework

From the above analysis, it may be discovered that the difference among the classical MDS method, the modified MDS method, and the subspace method lies in selecting coordinate reference points when constructing the inner product matrix. In fact, according to the definition of inner product matrix, any point within the space encompassing the target nodes and all anchor nodes can be taken as the coordinate reference point. Namely, the coordinate reference point can be expressed as

$$\begin{array}{cc}& x_r=\left[\begin{array}{cccc}x& x_{\mathrm{1}}& \cdots & x_N\end{array}\right]W,\end{array}$$

(34)

where

$$\begin{array}{cc}& W={\left[\begin{array}{cccc}w& w_{\mathrm{1}}& \cdots & w_N\end{array}\right]}^T\end{array}$$

(35)

is a $\left(N+\mathrm{1}\right)\times \mathrm{1}$ weight vector, and $W^T{\mathrm{1}}_{N+\mathrm{1}}=\mathrm{1}$ , whereby the determined coordinate reference point is called generalized mass center.

With respect to a given generalized mass-center reference point, the coordinate matrix composed of target node and anchor nodes can be expressed as

$$\begin{array}{cc}& X_r={\left[\begin{array}{cccc}x-x_r& x_{\mathrm{1}}-x_r& \cdots & x_N-x_r\end{array}\right]}^T.\end{array}$$

(36)

Namely,

$$\begin{array}{ll}{X}_r& ={\left[\begin{array}{llll}x& x_{\mathrm{1}}& \cdots & x_N\end{array}\right]}^T-{\mathrm{1}}_{N+\mathrm{1}}{W}^T{\left[\begin{array}{llll}x& x_{\mathrm{1}}& \cdots & x_N\end{array}\right]}^T\\ & =P_w{\left[\begin{array}{llll}x& x_{\mathrm{1}}& \cdots & x_N\end{array}\right]}^T,\end{array}$$

(37)

where

$$\begin{array}{cc}& P_w=I_{N+\mathrm{1}}-{\mathrm{1}}_{N+\mathrm{1}}{W}^T.\end{array}$$

(38)

Without measurement error, the inner product matrix under the generalized mass center is

$$\begin{array}{ll}{B}_r& =X_r{X}_r^T\\ & =P_w{\left[\begin{array}{llll}x& x_{\mathrm{1}}& \cdots & x_N\end{array}\right]}^T\left[\begin{array}{llll}x& x_{\mathrm{1}}& \cdots & x_N\end{array}\right]P_w^T.\end{array}$$

(39)

The method of getting the inner product matrix $B_r$ out of the Euclidean distance matrix is presented in the following. Note

$$\begin{array}{cc}& P_w{\mathrm{1}}_{N+\mathrm{1}}={\mathrm{0}}_{N+\mathrm{1}}.\end{array}$$

(40)

Then, for $X_c$ in Formula (1), there is

$$\begin{array}{cc}& P_w{X}_c{X}_c^T{P}_w^T=P_w{\left[\begin{array}{cccc}x& x_{\mathrm{1}}& \cdots & x_N\end{array}\right]}^T\left[\begin{array}{cccc}x& x_{\mathrm{1}}& \cdots & x_N\end{array}\right]P_w^T.\end{array}$$

(41)

Compare Formulas (39) with (41) to get

$$\begin{array}{cc}& B_r=P_w{X}_c{X}_c^T{P}_w^T.\end{array}$$

(42)

Refer to Formula (7) and $B_c=X_c{X}_c^T$ ; we can get

$$\begin{array}{cc}& B_r=P_w{B}_c{P}_w^T=-\frac{\mathrm{1}}{\mathrm{2}}\left(P_{w}J\right)D{\left(P_{w}J\right)}^T.\end{array}$$

(43)

It is easy to prove

$$\begin{array}{cc}& P_{w}J=P_w.\end{array}$$

(44)

Therefore, for reference points determined by any weight vector that satisfies $W^T{\mathrm{1}}_{N+\mathrm{1}}=\mathrm{1}$ , the inner product matrix can be obtained directly out of the Euclidean distance matrix; namely,

$$\begin{array}{cc}& B_r=-\frac{\mathrm{1}}{\mathrm{2}}{P}_{w}D{P}_w^T.\end{array}$$

(45)

Evidently, the double central transformation is a special instance of the above transformation.

Refer to the noise subspace $V_{r\mathrm{2}}$ of $B_r$ to get

$$\begin{array}{cc}& V_{r\mathrm{2}}{V}_{r\mathrm{2}}^T{X}_r={\mathrm{0}}_{(N+\mathrm{1})\times \mathrm{2}}.\end{array}$$

(46)

Then, arrange the above expression according to Formula (36) to get

$$\begin{array}{cc}& V_{r\mathrm{2}}{V}_{r\mathrm{2}}^T{Q}_r{x}^T=V_{r\mathrm{2}}{V}_{r\mathrm{2}}^T{P}_r,\end{array}$$

(47)

where

$$\begin{array}{cc}& Q_r={\left[\begin{array}{c}w-\mathrm{1}\\ w\\ ⋮\\ w\end{array}\right]}_{\left(N+\mathrm{1}\right)\times \mathrm{1}},\\ P_r={\left[\begin{array}{c}-{\displaystyle \sum }_{i=\mathrm{1}}^N{w}_i{x}_i^T\\ x_{\mathrm{1}}^T-{\displaystyle \sum }_{i=\mathrm{1}}^N{w}_i{x}_i^T\\ ⋮\\ x_N^T-{\displaystyle \sum }_{i=\mathrm{1}}^N{w}_i{x}_i^T\end{array}\right]}_{(N+\mathrm{1})\times \mathrm{2}}.\end{array}$$

(48)

When there exists measurement noise, by Formula (47), the least square estimation of the target node's position can be determined as

$$\begin{array}{cc}& \widehat{x}=\frac{P_r^T{\tilde{V}}_{r\mathrm{2}}{\tilde{V}}_{r\mathrm{2}}^T{Q}_r}{Q_r^T{\tilde{V}}_{r\mathrm{2}}{\tilde{V}}_{r\mathrm{2}}^T{Q}_r},\end{array}$$

(49)

where ${\tilde{V}}_{r\mathrm{2}}$ is the estimator of $V_{r\mathrm{2}}$ . The solution derived from the above formula is also a general solution of MDS method under the unified framework.

Under the above unified framework, the proposed method may change by selecting a different weight vector W. For instance, let

$$\begin{array}{cc}& W=\frac{\mathrm{1}}{N+\mathrm{1}}{\mathrm{1}}_{N+\mathrm{1}}\end{array}$$

(50)

to get the classical MDS method with all nodes’ mass centers as coordinate reference points. The solution of classical MDS method can be acquired by Formula (49). Consider

$$\begin{array}{cc}& w=\mathrm{0},\hspace{1em}\hspace{1em}{w}_{\mathrm{1}}=w_{\mathrm{2}}=\cdots =w_N=\frac{\mathrm{1}}{N}.\end{array}$$

(51)

The inner product matrix constructed by Formula (51) corresponds to the one by the modified MDS method (18), whereby the solution of the modified MDS method can be obtained. If the coordinate of the target node is selected as the coordinate reference point, namely,

$$\begin{array}{cc}& w=\mathrm{1},\hspace{1em}\hspace{1em}{w}_{\mathrm{1}}=w_{\mathrm{2}}=\cdots =w_N=\mathrm{0},\end{array}$$

(52)

the derived solution agrees with the solution derived from the subspace method. Apart from the above methods, a larger variety of MDS methods can be derived from the definition of generalized mass center.

4. Minimum Residual MDS Algorithm under Unified Framework

From the above analysis, it can be known that the solving process of MDS algorithm under the unified framework is actually to find the least square solution. However, the least square algorithm inevitably introduces unavoidable errors during the linearization process. To optimize the performance of MDS method, this paper proposes a minimum residual based MDS algorithm, extends MDS method to the circumstance of generalized mass center, and gets the minimum residual solution under the unified framework.

4.1. Theory of Minimum Residual Localization

The idea of minimum residual localization is to find a coordinate point within the region to be located so as to minimize the residual sum of squares of the distances between it and anchor nodes AN and the corresponding measured distances, namely, to maximize the matched degree between this coordinate point and the measured distances, and finally to take the coordinate point as the estimated position of the target.

According to what is assumed in the above experimental system, for N anchor nodes, the distance between the target node and the ith anchor node containing measurement noise errors is expressed as $r_i$ as in Formula (5), throughout the localization process, where $d_i$ denotes the real distance between TN and $\text{A}{\text{N}}_i$ and ${\eta }_i$ denotes the measurement noise. In practical application, this measurement noise is dependent on complicated factors such as hardware device, ranging technique, and environment. Its distribution parameters are difficult to obtain.

According to the idea of minimum residual localization, a target function can be built:

$$\begin{array}{cc}& f\left(x\right)={\displaystyle \sum }_i^N{\left\{‖x-\text{A}{\text{N}}_i‖-r_i\right\}}^{\mathrm{2}},\end{array}$$

(53)

where x is a certain coordinate point $(a,b)$ within the space and $\left|\right|x-\text{A}{\text{N}}_i\left|\right|$ is the Euclidean distance between the target x and anchor node $\text{A}{\text{N}}_i(a_i,b_i)$ . Consider

$$\begin{array}{cc}& ‖x-\text{A}{\text{N}}_i‖=\sqrt{{\left(a-a_i\right)}^{\mathrm{2}}+{\left(b-b_i\right)}^{\mathrm{2}}}.\end{array}$$

(54)

Then, function $f(x)$ represents the residual sum of squares of the distances between a certain coordinate point within the space and anchor nodes and the measured distances. If no noise exists along the measured distance,

$$\begin{array}{cc}& \underset{w_i\to \mathrm{0}}{\lim }f\left(x\right)={\displaystyle \sum }_i^N{\left\{‖x-\text{A}{\text{N}}_i‖-r_i\right\}}^{\mathrm{2}}={\displaystyle \sum }_i^N{\left\{‖x-\text{A}{\text{N}}_i‖-r_i\right\}}^{\mathrm{2}}.\end{array}$$

(55)

When ${\lim }_{w_i\to \mathrm{0}}f(x)=\mathrm{0}$ , the x coordinate is the real position of target node. However, due to the existence of measuring error, in general cases $f(x)\ne \mathrm{0}$ . Therefore, x which minimizes $f(x)$ may be taken as the estimated position of the target; namely,

$$\begin{array}{cc}& \tilde{x}=\arg \hspace{0.17em}\min f\left(x\right).\end{array}$$

(56)

4.2. Application of Particle Swarm Optimization (PSO) Algorithm in Minimum Residual Localization

According to the above analysis, the minimum residual localization algorithm is to find the function in an unconstrained optimization problem. Consider

$$\begin{array}{ll}\tilde{x}& =\arg \hspace{0.17em}\min f\left(x\right)\\ & =\arg \hspace{0.17em}\min {\displaystyle \sum }_i^N{\left\{\sqrt{{\left(a-a_i\right)}^{\mathrm{2}}+{\left(b-b_i\right)}^{\mathrm{2}}}-r_i\right\}}^{\mathrm{2}}.\end{array}$$

(57)

Given the coordinates of all anchor nodes AN within the positioning system, the measured distance from the target node TN to each AN is to be obtained in each locating cycle.

The idea of PSO algorithm is to adopt multiple particles to simulate the target's location. The updating direction of the particles will also be determined by comparing their individual optimum with group optimum. Through several iterations, the group optimum of the particle swarm is selected as the optimal solution. Assume that the number of particles is M and the maximum count of iteration is K. The PSO-based minimum residual localization algorithm works throughout the solving procedure as follows.

Step 1 (initialization).

Randomly generate particles and their initial speeds within the positioning space:

$$\begin{array}{cc}& x_m^{\mathrm{1}}=\left[a_m^{\mathrm{1}},b_m^{\mathrm{1}}\right],\\ V_m^{\mathrm{1}}=\left[v{a}_m^{\mathrm{1}},v{b}_m^{\mathrm{1}}\right],\end{array}$$

(58)

Compute the adaptive value $f(x_m^{\mathrm{1}})$ of each particle, which is taken as the initial historical optimal solution of this particle; namely,

$$\begin{array}{cc}& B_m^{\mathrm{1}}=x_m^{\mathrm{1}}.\end{array}$$

(59)

Meanwhile, the historical optimal solution of the particle swarm is the particle which has achieved the minimum adaptive value:

$$\begin{array}{cc}& B_g^{\mathrm{1}}=\underset{m=\mathrm{1}:M}{\arg \hspace{0.17em}\min }f\left(x_m^{\mathrm{1}}\right).\end{array}$$

(60)

Step 2 (the kth iterative update of particles ( $\mathrm{2}\le k\le K$ )).

Update the particles according to their individual historical optimal solution and group historical optimal solution:

$$\begin{array}{cc}& V_m^k=\omega V_m^{k-\mathrm{1}}+c_{\mathrm{1}}\xi \left(B_m^{k-\mathrm{1}}-x_m^{k-\mathrm{1}}\right)+c_{\mathrm{2}}\mu \left(B_g^{k-\mathrm{1}}-x_m^{k-\mathrm{1}}\right),\\ x_m^k=x_m^{k-\mathrm{1}}+V_m^{k-\mathrm{1}},\end{array}$$

(61)

where ω is called inertia weight, $c_{\mathrm{1}}$ and $c_{\mathrm{2}}$ are called learning factors, and $\xi ,\mu \in U[\mathrm{0,1}]$ are pseudorandom numbers following uniform distribution within the interval $[0,1]$ .

Step 3.

Update the individual historical optimal solution and group historical optimal solution of the particles:

$$\begin{array}{cc}& B_m^k=\left\{\begin{array}{cc}{B}_m^{k-\mathrm{1}},\hfill & f\left(x_m^k\right)>f\left(B_m^{k-\mathrm{1}}\right)\hfill \\ x_m^k,\hfill & f\left(x_m^k\right)\le f\left(B_m^{k-\mathrm{1}}\right)\hfill \end{array}\right.\\ B_g^k=\left\{\begin{array}{cc}{B}_g^{k-\mathrm{1}},\hfill & \underset{m=\mathrm{1}:M}{\min }f\left(x_m^k\right)>f\left(B_g^{k-\mathrm{1}}\right)\hfill \\ \underset{m=\mathrm{1}:M}{\arg \hspace{0.17em}\min }f\left(x_m^k\right),\hfill & \underset{m=\mathrm{1}:M}{\min }f\left(x_m^k\right)\le f\left(B_g^{k-\mathrm{1}}\right).\hfill \end{array}\right.\end{array}$$

(62)

Step 4.

Decide whether the maximum count of iteration K is achieved.

If no, let $k=k+1$ and skip to Step 2.

If yes, the initially estimated position of the target at the moment is

$$\begin{array}{cc}& \tilde{x}\left(t\right)=B_g^K.\end{array}$$

(63)

The dimensionality is n. In the process of initialization, the time complexity is $O\left(Mn\right)$ . The time complexity is $O\left(Mn\right)$ when the fitness function of M particles is calculated. The time complexity is $O\left(M\right)$ when the individual value of M particles is updated. The time complexity is $O\left(M\right)$ when the best value is chosen from M individual extremums as the global extremum. The time complexity is $O\left(Mn\right)$ when the velocity and position of each particle are updated. In the above steps, the worst time complexity is $O\left(Mn\right)$ . The time complexity after circulation is $O\left(\text{ite}{\text{r}}_{\max }Mn\right)$ when the iteration time of algorithm is $\text{ite}{\text{r}}_{\max }$ . Therefore, the time complexity of PSO algorithm is $O\left(\text{ite}{\text{r}}_{\max }Mn\right)$ .

It has been proven that $x_m$ of PSO algorithm converges to the center weighted of $B_g$ and $B_m$ when $\sqrt{2\left(1+\omega -\phi \right)-4\omega }<2\hspace{0.17em}\hspace{0.17em}(\phi ={\phi }_1+{\phi }_2,\hspace{0.17em}\hspace{0.17em}{\phi }_1=c_1\xi ,\hspace{0.17em}\hspace{0.17em}{\phi }_2=c_2\mu )$ and $B_g,B_m$ remain unchanged in the process of evolution; namely, $x_m\to \left({\phi }_1{B}_m+{\phi }_2{B}_g\right)/\phi$ [21]. In fact, $B_g$ and $B_m$ are constantly updated according to their adaptive values in the process of evolution. It also has been proven that the convergence of the PSO algorithm can be completely guaranteed if only the parameter chosen by us meet the above-mentioned requirements [22].

4.3. MDS-Based Minimum Residual Localization Algorithm (MinRE MDS)

In fact, Formula (47) can be rewritten as

$$\begin{array}{cc}& \left[\begin{array}{c}1\\ x\end{array}\right]Q_r^T{V}_{r\mathrm{2}}{V}_{r\mathrm{2}}^T=\left[\begin{array}{c}{Q}_r^T\\ P_r^T\end{array}\right]V_{r\mathrm{2}}{V}_{r\mathrm{2}}^T.\end{array}$$

(64)

Postmultiply both sides of the above equation by vector $Q_r$ , and divide it by $Q_r^T{V}_{r\mathrm{2}}{V}_{r\mathrm{2}}^T{Q}_r$ to produce

$$\begin{array}{cc}& \left[\begin{array}{c}1\\ x\end{array}\right]=\left[\begin{array}{c}{Q}_r^T\\ P_r^T\end{array}\right]\frac{V_{r\mathrm{2}}{V}_{r\mathrm{2}}^T{Q}_r}{Q_r^T{V}_{r\mathrm{2}}{V}_{r\mathrm{2}}^T{Q}_r}.\end{array}$$

(65)

Namely,

$$\begin{array}{cc}& {\left[\begin{array}{c}{Q}_r^T\\ P_r^T\end{array}\right]}^{†}\left[\begin{array}{c}1\\ x\end{array}\right]=\frac{V_{r\mathrm{2}}{V}_{r\mathrm{2}}^T{Q}_r}{Q_r^T{V}_{r\mathrm{2}}{V}_{r\mathrm{2}}^T{Q}_r}.\end{array}$$

(66)

Here, the superscript † denotes the Moore-Penrose inverse. Assuming that the signal subspace of inner product matrix $B_r$ is $V_{r\mathrm{1}}$ , we can get $B_r=V_{r\mathrm{1}}{\Lambda }_{r\mathrm{1}}{V}_{r\mathrm{1}}^T$ . Given the orthogonality of $V_{r\mathrm{1}}$ and $V_{r\mathrm{2}}$ , premultiply both sides of Formula (66) by $B_r$ to produce

$$\begin{array}{cc}& B_r{\left[\begin{array}{c}{Q}_r^T\\ P_r^T\end{array}\right]}^{†}\left[\begin{array}{c}\mathrm{1}\\ x\end{array}\right]={\mathrm{0}}_{N+\mathrm{1}}.\end{array}$$

(67)

When measurement noise exists, there is

$$\begin{array}{cc}& \tilde{D}=D+\Delta D,\\ {\tilde{B}}_r=B_r+\Delta B_r.\end{array}$$

(68)

Ignore the high-order terms to get

$$\begin{array}{cc}& \Delta D\approx \left[\begin{array}{cccc}\mathrm{0}& \mathrm{2}{d}_{\mathrm{1}}{\eta }_{\mathrm{1}}& \cdots & \mathrm{2}{d}_N{\eta }_N\\ \mathrm{2}{d}_{\mathrm{1}}{\eta }_{\mathrm{1}}& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{2}{d}_N{\eta }_N& 0& \cdots & 0\end{array}\right],\\ \Delta B_r=-\frac{\mathrm{1}}{\mathrm{2}}{P}_w\Delta D{P}_w^T.\end{array}$$

(69)

Thus, the residual vector ε of Formula (67) can be expressed as

$$\begin{array}{ll}\epsilon & =\Delta B_r{\left[\begin{array}{l}{Q}_r^T\\ P_r^T\end{array}\right]}^{†}\left[\begin{array}{l}\mathrm{1}\\ x\end{array}\right]\\ & =-\frac{\mathrm{1}}{\mathrm{2}}{P}_w\Delta D{P}_w^T{\left[\begin{array}{l}{Q}_r^T\\ P_r^T\end{array}\right]}^{†}\left[\begin{array}{l}\mathrm{1}\\ x\end{array}\right].\end{array}$$

(70)

Now, the sum of residual vectors ε makes the same physical sense as Formula (53).

Let

$$\begin{array}{cc}& A_r={\left[\begin{array}{c}{Q}_r^T\\ P_r^T\end{array}\right]}^{†}.\end{array}$$

(71)

And divide $A_r$ into

$$\begin{array}{cc}& A_r=\left[\begin{array}{ccc}{A}_{r\mathrm{1}}& \cdots & A_{r\mathrm{2}}\end{array}\right]\end{array}$$

(72)

by column, where $A_{r\mathrm{1}}$ contains only the first column of $A_r$ , while $A_{r\mathrm{2}}$ contains the other columns of $A_r$ . Then, the least square solution of target node's position is

$$\begin{array}{cc}& \widehat{x}=-\left(A_{r\mathrm{2}}^T{\tilde{B}}_r^T{\tilde{B}}_r{A}_{r\mathrm{2}}\right)A_{r\mathrm{2}}^T{\tilde{B}}_r^T{\tilde{B}}_r{A}_{r\mathrm{1}}.\end{array}$$

(73)

In order to solve the optimal position of the target, substitute $\widehat{x}$ obtained from Formula (73) into (57) to perform optimal iterations over the particle swarm. When the minimum residual solution is solved, the ultimate coordinate estimation of mobile nodes can be obtained. The classical MDS method, the modified MDS method, and the subspace-based method under the unified framework are derived by selecting a different weight vector W.

The above procedure can avoid eigen decomposition in the classical MDS algorithm, overcome errors introduced in least square algorithm as a result of linearization, and work out MDS-based minimum solution of the target coordinate.

5. Experimental Verification and Effect Analysis

5.1. Experimental Platform Construction

To verify the performance of the minimum residual MDS localization algorithm under the unified framework proposed in this paper, confirmatory experiment is implemented at Lab 104 in Science Building of Northeastern University (NEU). The whole experimental environment is arranged as in Figure 1. The environment extends a range of 10 m ∗ 6 m. The positions of the anchor nodes are AN1(8, 1), AN2(1, 1), AN3(1, 5), AN4(8, 5), and AN5(4.5, 3), respectively. The motion trail of the target, which moves 6.5 m ahead, turns right and moves 3 m ahead before continuing to turn right and move 6.5 m ahead, and starts at the point (8,1.5). The total length of the trail covered by the target is 16 m; the moving speed is 1 m/s; the sampling frequency of RSSI values is twice per second (0.5 Hz). The total sampling number is 33.

Figure 1

Environment deployment.

In the experiment, ZigBee nodes are adopted as anchor nodes and target node in WSN. The nodes are arranged above a tripod with the height of 1.6 m. The target node is carried above the mobile robot iRobot, as shown in Figure 2.

Figure 2

ZigBee nodes deployment.

5.2. Localization Result and Analysis

Iterative recursive weighted average filter [23] is adopted to process the initial measured values of RSSI, which are converted into estimated values of distance via RSSI signal attenuation model. After obtaining the estimated distances between the anchor nodes and target node, estimation of node positions can be performed.

A comparison is made between the minimum residual classical MDS method, minimum residual modified MDS method, and minimum residual subspace method under the unified framework and the classical MDS method and modified MDS method as well as subspace method under the unified framework. For the experimental environment and data at this time, the obtained result of localization is shown in Figures 3–8. According to the localization effect diagram, it can be seen that there exists a larger error within the least square MDS localization in which there exists a great deflection between the estimation result and real value, whereas the PSO-based minimum residual MDS algorithm gives significantly smaller bias in position estimation.

Figure 3

Classic MDS localization result.

Figure 4

MinRE classic MDS localization result.

Figure 5

Modified MDS localization result.

Figure 6

MinRE modified MDS localization result.

Figure 7

Subspace localization result.

Figure 8

MinRE subspace localization result.

Figures 9 and 10 are analysis diagrams of error statistics. It can be easily seen that the error of the PSO-based minimum residual MDS localization algorithm is significantly smaller than that of the least square MDS localization algorithm. In other words, it can achieve a more perfect effect of localization.

Figure 9

Localization error.

Figure 10

Cumulative distribution function.

The concrete parameter analysis of localization error is shown in Table 1, which reveals that the minimum residual MDS localization algorithm has enhanced accuracy dramatically as compared with the least square MDS. Experimental results show that the classical MDS method, modified MDS method, and subspace method under the unified framework have demonstrated different performances in accordance with different selections of coordinate reference points. However, there is not a selection scheme of coordinate reference points which outperforms the minimum residual MDS algorithm under the same reference points. The minimum residual MDS method significantly outperforms the least square MDS method, with better localization accuracy up to 41%. Furthermore, the former performance is immune to selection of reference nodes.

Table 1

Localization error statistic.

Error	Mean	Variance	Maximum
Classic MDS	1.4527	0.4819	2.6186
Min. classic MDS	0.8458	0.4039	1.6942
Modified MDS	1.4825	0.5168	2.7258
Min. modified MDS	0.8661	0.4132	1.7080
Subspace	2.7742	1.7439	9.2058
Min. subspace	0.8498	0.4129	1.7132

6. Conclusion

In this paper, three available MDS algorithms, including the classical MDS algorithm, modified MDS algorithm, and subspace method, are analyzed. Their difference mainly lies in the selection of reference nodes. Thereby, the unified framework is used to solve MDS algorithms. On the basis of existing methods, MDS localization algorithm based on minimum residual is put forward. Its general solution is provided under the above unified framework. In order to verify the performance of the proposed algorithm, a small-scale WSN network with ZigBee nodes is adopted to implement practical experiment. Experiment results show that the minimum residual MDS algorithm proposed in this paper can effectively enhance the localization accuracy of moving target. What is more, the performance of the minimum residual MDS algorithm is immune to the selection of reference nodes.