Localization for mobile sensor network based on unscented Kalman filter with clipping and uncorrelated conversion

Abstract

This article studies localization for mobile sensor network with incomplete range measurements, that is, the ranges between some pairs of sensors cannot be measured. Different from existing works that localize sensors one by one, the localization problem in this article is solved for the whole mobile sensor network. According to whether the ranges can be measured or not, all the sensors in the network are grouped to construct basic localization units. For the sensors in basic localization units, a constrained nonlinear model is first established to formulate their relative motion, where the motion states are chosen as ranges and cosine values of angles between ranges. Then, based on the established model, a constrained unscented Kalman filter is adopted to provide motion state estimation. In the constrained unscented Kalman filter, the clipping technique is introduced to handle the model constraints, and the uncorrelated conversion technique is introduced to make full use of measurements. Hence, the estimation accuracy can be improved. Finally, the distributed multidimensional scaling-map method is used to localize the whole sensor network using the estimated ranges, and a localization algorithm is presented. The effectiveness and advantages of the proposed algorithm are demonstrated through several simulation examples.

Introduction

Wireless sensor network (WSN) is a research hotspot in recent years and has been widely used in battlefield environment awareness, environment monitoring, robots cooperative control, unmanned combat and other fields. ¹ As a prerequisite and basis procedure in the WSN applications, the localization of WSN has attracted more and more researchers whose attention is concentrated on improving localization accuracy and effectiveness.

According to the manoeuvrability of sensors, localization of WSN can be classified into two categories, that is, static sensor network localization and mobile sensor network localization. In the former category, the stationary sensors are considered, and lots of effective algorithms have been put forward, see, for example, time of arrival (TOA) localization algorithm, ² time difference of arrival (TDOA) localization algorithm, ³ angle of arrival (AOA) localization algorithm ⁴ and received signal strength indicator (RSSI) localization algorithm.^5–7 To improve the localization accuracy, some researches attend to introduce filters to design localization algorithms. The nonlinear filters (extended Kalman filter (EKF) and unscented Kalman filter (UKF)) and the multidimensional scaling-map (MDS-MAP) method are combined to design localization algorithms (see MDS-EKF and MDS-UKF algorithms), ⁸ and the accuracy comparison of different algorithms is provided. However, neither the MDS-UKF nor MDS-EKF algorithms can effectively deal with the sensors with manoeuvrability and cannot be applied for mobile sensor network localization.

For the mobile sensor network, the manoeuvrability of sensors leads to time-varying structure of WSN, hence brings more challenges to the localization with high accuracy. In the localization of mobile sensor network, the Global Positioning System (GPS)/inertial navigation system (INS) information are generally used, and many algorithms have been proposed, see for example, localization algorithm based on joint distribution state-information filter ⁹ and localization algorithm via an inertial navigation system–assisted single roadside unit. ¹⁰ However, the GPS/INS devices are not only expensive but also unreliable in some practical conditions.

Without using the GPS/INS information, the mobility of sensors is considered and the localization algorithm can be designed by introducing particle filter; see the Monte Carlo localization (MCL) algorithm ¹¹ and other improved MCL-based algorithms.^11–18 In these MCL-based algorithms, the manoeuvrability of the mobile sensors to be localized is handled by generating particles with random directions and velocities, and the detection radii of the anchor sensors (i.e. sensors with known positions) are used to filter out the impossible particles. Then, the approximate positions of mobile sensors can be calculated by fusing the positions of possible particles. However, because some detailed measurements (e.g. ranges, angles) are not used by MCL-based algorithms, the accurate positions of mobile sensors cannot be obtained.

Using the ranges and angles between the sensors to be localized and the anchor sensors, some other filters can be adopted in the localization algorithm design, see, for example, the square root-UKF (SR-UKF) ¹⁹ and the SR-cubature Kalman filter (SR-CKF). ²⁰ However, in some general practical conditions, the angles between sensors cannot be measured, even the ranges between some pairs of sensors are unavailable (i.e. the range measurements are incomplete). For these conditions, the above localization algorithms are not suitable. Motivated by this situation, this article proposes to design a localization algorithm with incomplete range measurements. Considering that some sensors to be localized may not be detected by the anchor sensors, the localization of these sensors cannot be achieved independently. Hence, the clustering and merging methods, which are illustrated to be effective for static sensor network, ²¹ are introduced into mobile sensor network in this article. This contributes to achieving overall localization of mobile sensor network. By first grouping the whole sensor network into some basic localization units, the ranges between sensors can be compensated and refined with filtering. According to the compensated and refined ranges, the local sensor networks are obtained by clustering the whole sensor network, and then be localized and merged to achieve overall localization. Following the manners of the proposed algorithm, better performance of localization is obtained.

The rest of this article is organized as follows. The localization problem for mobile sensor network is first stated in section ‘Problem statement’. Section ‘Localization for mobile sensor network’ details the design procedure of localization algorithm, with establishing the motion model in section ‘Relative motion model of basic localization unit’, proposing the filter design in ‘UKF with clipping and uncorrelated conversion’, clustering and merging for the whole sensor network in section ‘Clustering and merging’, and the complete localization algorithm in section ‘Mobile sensor network localization algorithm with incomplete measurement’. Finally, section ‘Simulation’ provides a demonstration of the advantages of the designed algorithm through several simulation examples, and section ‘Conclusion’ summarizes this article with a conclusion.

Notation. For a WSN of N sensors, we denote $\mathbb{N}=\{1,2,\dots ,N\}$ as its index set. ${\mathbb{R}}^n$ and ${\mathbb{R}}^{m\times n}$ represent the n-dimensional Euclidean space and $m\times$ n-dimensional real matrix set, respectively. ${\mathbb{I}}_{a:b}=\{a,a+1,\dots ,b-1,b\}$ denotes the integer set from a to b. ${\mathit{I}}_n$ denotes the n-dimensional identity matrix, and $0_{m\times n}$ the $m\times n$ -dimensional zero matrix. For a variable x, $\mathrm{sat}(x)$ denotes its saturation with magnitude 1, that is, $\mathrm{sat}(x)=\mathrm{sgn}(x)min\{1,\left|x\right|\}$ with $\mathrm{sgn}(·)$ being a sign function. $N(0,\mathit{R})$ denotes the normal distribution with zero mean and covariance matrix R . $\mathrm{diag}\{a_1,a_2,\dots \}$ denotes a diagonal matrix with diagonal entries $a_1,a_2,\dots$ .

Problem statement

Consider a mobile sensor network with N sensors shown in Figure 1, where the sensors move in three-dimensional (3D) space with unknown motion. At time t, the actual ranges among these sensors can be measured as

$$r_t^{\mathrm{ij}}=d_t^{\mathrm{ij}}+e_t^{\mathrm{ij}}\ge 2\rho ,\forall i,j\in \mathbb{N}$$

(1)

where $r_t^{\mathrm{ij}}\in \mathbb{R}$ , $d_t^{\mathrm{ij}}\in \mathbb{R}$ and $e_t^{\mathrm{ij}}~N(0,\mathit{R})$ denote the measured range, actual range and measurement error between sensors i and j, respectively, and $\rho$ denotes the radius of each sensor. Due to some practical conditions such as range limitation, electronic interference and obstruction, the ranges among some sensors cannot be measured, that is, the range measurements are incomplete.

Figure 1.

Mobile sensor network of N sensors moving with unknown motion.

Based on the incomplete range measurements, the localization problem to be solved in this article is to calculate the absolute positions of all sensors $i\in \mathbb{N}$ (denoted as ${\mathit{x}}_i\in {\mathbb{R}}^3$ ). To guarantee the practicability of localization without loss of generality, we assume that there are at least four anchor sensors indexed by $i_a,j_a,\dots ,\in {\mathbb{N}}_a$ , where ${\mathbb{N}}_a\subset \mathbb{N}$ denotes the index set of anchor sensors. Accordingly, the index set of sensors to be localized is denoted as ${\mathbb{N}}_l=\mathbb{N}\backslash {\mathbb{N}}_a$ .

For each sensor $i\in {\mathbb{N}}_l$ to be localized, most of the existing range-based approaches can calculate its absolute position ${\mathit{x}}_i$ based on the assumption that the information of at least one anchor sensor is available for sensor i. However, incomplete range measurements may violate this assumption and then disable these localization approaches. Motivated by this condition, this article contributes to present a localization approach for the mobile sensor network with incomplete range measurements.

Localization for mobile sensor network

Using the incomplete measured ranges and the positions of anchor sensors ${\mathbb{N}}_a$ , the whole mobile sensor network can be overall localized by introducing clustering and merging. Whereas these measured ranges are not precise due to ranging errors and mobility of sensors. By directly using these measured ranges, the localization with clustering and merging will result in low accuracy. Aiming at alleviation of these errors in measured ranges, a proper filter is designed to calculate more accurate estimated ranges. Out of requirement of filtering, a motion model should be first established to formulate the relative motion of sensor network. However, the motion model with high dimensions will lead to tremendous computational complexity of filtering. Therefore, the whole mobile sensor network $\mathbb{N}$ is preferred to be grouped into M basic localization units ${\mathbb{N}}_p^{\mathrm{unit}}$ ’s, $p\in {\mathbb{I}}_{1:M}$ . For each ${\mathbb{N}}_p^{\mathrm{unit}}$ , the relative motion model is established as follows.

Relative motion model of basic localization unit

Based on graph rigidity, at least four common sensors between clusters are required for clustering and merging in 3D space. Thus, five sensors are selected to construct a basic localization unit ${\mathbb{N}}_p^{\mathrm{unit}}$ as the most simplified scheme.

For each basic localization unit at time t, the relative motion model can be established in two different forms:

• Model (1) relative range model, which describes the relative motion between the ranges among sensors.

• Model (2) relative coordinate model, which describes the relative motion between the position coordinates of sensors.

Compared with model (1), the nonlinearity of model (2) is more significant, which introduces more complexity in the filtering process and leads to less precision of the localization results. Therefore, for each basic localization unit ${\mathbb{N}}_p^{\mathrm{unit}}$ , model (1) is adopted to construct a relative motion model.

As shown in Figure 2, a basic localization unit composed of five sensors (i.e. A, B, C, D and E) is denoted as ${\mathbb{N}}_p^{\mathrm{unit}}=\{A,B,C,D,E\}$ ; the angles between sensors are denoted as ${\theta }_t^1=\angle \mathrm{ADC}$ , ${\theta }_t^2=\angle \mathrm{ADB}$ , ${\theta }_t^3=\angle \mathrm{BDC}$ , ${\theta }_t^4=\angle \mathrm{BAC}$ , ${\theta }_t^5=\angle \mathrm{AEB}$ , ${\theta }_t^6=\angle \mathrm{AEC}$ , ${\theta }_t^7=\angle \mathrm{BEC}$ and ${\theta }_t^8=\angle \mathrm{AED}$ at time t. As the sensors move from positions $A,B,C,D,E$ at time t to positions $A\prime ,B\prime ,C\prime ,D\prime ,E\prime$ at time $t\prime =t+\Delta t$ , the relative motion between them can be modelled in continuous-time form as (more details of the model establishment are provided in Appendix 1)

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{\stackrel{~}{\mathit{X}}}}_t=[\begin{array}{ccccccccc}{\stackrel{·}{\xi }}_t^1\hfill & \cdots \hfill & {\stackrel{·}{\xi }}_t^{10}\hfill & {\stackrel{·}{\xi }}_t^{11}\hfill & \cdots \hfill & {\stackrel{·}{\xi }}_t^{20}\hfill & {\stackrel{·}{\xi }}_t^{21}\hfill & \cdots \hfill & {\stackrel{·}{\xi }}_t^{28}\hfill \end{array}{]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{=\left[\begin{array}{c}{\xi }_t^{11}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\xi }_t^{20}}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{\frac{{\xi }_t^3{\xi }_t^{13}+{\xi }_t^8{\xi }_t^{18}-{\xi }_t^2{\xi }_t^{12}-({\xi }_t^3{\xi }_t^{18}+{\xi }_t^8{\xi }_t^{13}){\xi }_t^{21}}{{\xi }_t^3{\xi }_t^8}}\\ \multicolumn{1}{c}{\frac{{\xi }_t^3{\xi }_t^{13}+{\xi }_t^6{\xi }_t^{16}-{\xi }_t^1{\xi }_t^{11}-({\xi }_t^3{\xi }_t^{16}+{\xi }_t^6{\xi }_t^{13}){\xi }_t^{22}}{{\xi }_t^3{\xi }_t^6}}\\ \multicolumn{1}{c}{\frac{{\xi }_t^6{\xi }_t^{16}+{\xi }_t^8{\xi }_t^{18}-{\xi }_t^5{\xi }_t^{15}-({\xi }_t^6{\xi }_t^{18}+{\xi }_t^8{\xi }_t^{16}){\xi }_t^{23}}{{\xi }_t^6{\xi }_t^8}}\\ \multicolumn{1}{c}{\frac{{\xi }_t^1{\xi }_t^{11}+{\xi }_t^2{\xi }_t^{12}-{\xi }_t^5{\xi }_t^{15}-({\xi }_t^1{\xi }_t^{12}+{\xi }_t^2{\xi }_t^{11}){\xi }_t^{24}}{{\xi }_t^1{\xi }_t^2}}\\ \multicolumn{1}{c}{\frac{{\xi }_t^4{\xi }_t^{14}+{\xi }_t^7{\xi }_t^{17}-{\xi }_t^1{\xi }_t^{11}-({\xi }_t^4{\xi }_t^{17}+{\xi }_t^7{\xi }_t^{14}){\xi }_t^{25}}{{\xi }_t^4{\xi }_t^7}}\\ \multicolumn{1}{c}{\frac{{\xi }_t^4{\xi }_t^{14}+{\xi }_t^9{\xi }_t^{19}-{\xi }_t^2{\xi }_t^{12}-({\xi }_t^4{\xi }_t^{19}+{\xi }_t^9{\xi }_t^{14}){\xi }_t^{26}}{{\xi }_t^4{\xi }_t^9}}\\ \multicolumn{1}{c}{\frac{{\xi }_t^7{\xi }_t^{17}+{\xi }_t^9{\xi }_t^{19}-{\xi }_t^5{\xi }_t^{15}-({\xi }_t^7{\xi }_t^{19}+{\xi }_t^9{\xi }_t^{17}){\xi }_t^{27}}{{\xi }_t^7{\xi }_t^9}}\\ \multicolumn{1}{c}{\frac{{\xi }_t^4{\xi }_t^{14}+{\xi }_t^{10}{\xi }_t^{20}-{\xi }_t^3{\xi }_t^{13}-({\xi }_t^4{\xi }_t^{20}+{\xi }_t^{10}{\xi }_t^{14}){\xi }_t^{28}}{{\xi }_t^4{\xi }_t^{10}}}\end{array}\right]+\mathit{G}\left[\begin{array}{c}{\tilde{w}}_t^1\hfill \\ \multicolumn{1}{c}{{\tilde{w}}_t^2}\\ \multicolumn{1}{c}{{\tilde{w}}_t^3}\\ \multicolumn{1}{c}{{\tilde{w}}_t^4}\\ \multicolumn{1}{c}{{\tilde{w}}_t^5}\\ \multicolumn{1}{c}{{\tilde{w}}_t^6}\\ \multicolumn{1}{c}{{\tilde{w}}_t^7}\\ \multicolumn{1}{c}{{\tilde{w}}_t^8}\\ \multicolumn{1}{c}{{\tilde{w}}_t^9}\\ \multicolumn{1}{c}{{\tilde{w}}_t^{10}}\end{array}\right]}\\ \multicolumn{1}{c}{=f({\tilde{\mathit{X}}}_t)+\mathit{G}{\tilde{\mathit{W}}}_t}\end{array}$$

(2)

with the state vector

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\mathit{X}}}_t={[\begin{array}{ccccc}{\xi }_t^1\hfill & {\xi }_t^2\hfill & \cdots \hfill & {\xi }_t^{27}\hfill & {\xi }_t^{28}\hfill \end{array}]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{=[\begin{array}{ccccccc}{d}_t^{\mathrm{AB}}\hfill & d_t^{\mathrm{AC}}\hfill & d_t^{\mathrm{AD}}\hfill & d_t^{\mathrm{AE}}\hfill & d_t^{\mathrm{BC}}\hfill & d_t^{\mathrm{BD}}\hfill & d_t^{\mathrm{BE}}\hfill \end{array}}\\ \multicolumn{1}{c}{\begin{array}{ccccccc}{d}_t^{\mathrm{CD}}\hfill & d_t^{\mathrm{CE}}\hfill & d_t^{\mathrm{DE}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{AB}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{AC}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{AD}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{AE}}\hfill \end{array}}\\ \multicolumn{1}{c}{\begin{array}{cccccc}{\stackrel{·}{d}}_t^{\mathrm{BC}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{BD}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{BE}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{CD}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{CE}}\hfill & {\stackrel{·}{d}}_t^{\mathrm{DE}}\hfill \end{array}}\\ \multicolumn{1}{c}{\begin{array}{cccc}\cos ({\theta }_t^1)\hfill & \cos ({\theta }_t^2)\hfill & \cos ({\theta }_t^3)\hfill & \cos ({\theta }_t^4)\hfill \end{array}}\\ \multicolumn{1}{c}{\begin{array}{cccc}\cos ({\theta }_t^5)\hfill & \cos ({\theta }_t^6)\hfill & \cos ({\theta }_t^7)\hfill & \cos ({\theta }_t^8)\hfill \end{array}{]}^{\mathrm{T}}}\end{array}$$

(3)

the disturbance (seemed as the unknown relative acceleration)

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\mathit{W}}}_t={[\begin{array}{ccc}{\tilde{w}}_t^1\hfill & \cdots \hfill & {\tilde{w}}_t^{10}\hfill \end{array}]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{=[\begin{array}{ccccc}{\stackrel{··}{d}}_t^{\mathrm{AB}}\hfill & {\stackrel{··}{d}}_t^{\mathrm{AC}}\hfill & {\stackrel{··}{d}}_t^{\mathrm{AD}}\hfill & {\stackrel{··}{d}}_t^{\mathrm{AE}}\hfill & {\stackrel{··}{d}}_t^{\mathrm{BC}}\hfill \end{array}}\\ \multicolumn{1}{c}{\begin{array}{ccccc}{\stackrel{··}{d}}_t^{\mathrm{BD}}\hfill & {\stackrel{··}{d}}_t^{\mathrm{BE}}\hfill & {\stackrel{··}{d}}_t^{\mathrm{CD}}\hfill & {\stackrel{··}{d}}_t^{\mathrm{CE}}\hfill & {\stackrel{··}{d}}_t^{\mathrm{DE}}\hfill \end{array}{]}^{\mathrm{T}}}\end{array}$$

(4)

and the disturbance gain matrix

$$\mathit{G}=\left[\begin{array}{c}{0}_{10\times 10}\hfill \\ \multicolumn{1}{c}{1_{10\times 10}}\\ \multicolumn{1}{c}{0_{8\times 10}}\end{array}\right]$$

(5)

Figure 2.

The relative location of sensors in basic localization unit.

The measurement model can be formulated as

$${\tilde{\mathit{Z}}}_t=\mathit{H}{\tilde{\mathit{X}}}_t+{\tilde{\mathit{V}}}_t$$

(6)

with the measurement vector ${\tilde{\mathit{Z}}}_t=[\begin{array}{ccccccccc}{r}_t^{\mathrm{AB}}\hfill & r_t^{\mathrm{AC}}\hfill & r_t^{\mathrm{AD}}{r}_t^{\mathrm{AE}}\hfill & r_t^{\mathrm{BC}}\hfill & r_t^{\mathrm{BD}}\hfill & r_t^{\mathrm{BE}}\hfill & r_t^{\mathrm{CD}}\hfill & r_t^{\mathrm{CE}}\hfill & r_t^{\mathrm{DE}}\hfill \end{array}{]}^{\mathrm{T}}$ , the measurement matrix $\mathit{H}=[\begin{array}{cc}{\mathit{I}}_{10}\hfill & 0_{10\times 18}\hfill \end{array}]$ and measurement error vector ${\tilde{\mathit{V}}}_t=[\begin{array}{ccc}{e}_t^{\mathrm{AB}}\hfill & \dots \hfill & e_t^{\mathrm{DE}}\hfill \end{array}]~\mathrm{N}(0,\mathit{R})$ .

Remark 1

According to rigidity of graph, all the 10 ranges among unit ${\mathbb{N}}_p^{\mathrm{unit}}$ are required to uniquely determine its topology. In this article, these ranges can be completed by the combination of model (2) and filter if nine of them are measured. This advantage is introduced by considering the relationship between the ranges and cosine values of angles in model (2) (see equation (46) in Appendix 1). The more the angles’ cosine values are incorporated in state, the less the ranges are required to be measurable. However, the dimension increase in state will lead to heavy computation burden and less measurements of range will lead to reduction in localization accuracy.

Because the range between each pair of sensors is greater than $2\rho$ , and the cosine values of angles are between $[-1,1]$ , the state vector ${\tilde{\mathit{X}}}_t$ should satisfy the constraints

$$\{\begin{array}{cc}{\tilde{\mathit{X}}}_t(i)\ge 2\rho ,\hfill & i\in {\mathbb{I}}_{1:10}\hfill \\ \multicolumn{1}{c}{-1\le {\tilde{\mathit{X}}}_t(i)\le 1,}& i\in {\mathbb{I}}_{21:28}\hfill \end{array}$$

(7)

For the convenience of filtering, we prefer to discretize models (2) and (6) by sampling it at instant $t_k$ , as

$$\{\begin{array}{c}{\mathit{X}}_{k+1}=F({\mathit{X}}_k)+{\mathit{G}}_k{\mathit{W}}_k\hfill \\ \multicolumn{1}{c}{{\mathit{Z}}_k={\mathit{H}}_k{\mathit{X}}_k+{\mathit{V}}_k}\end{array}$$

(8)

with the discrete-time transition function $F(·)$ , the matrices ${\mathit{G}}_k=\delta \mathit{G}$ and ${\mathit{H}}_k=\mathit{H}$ , sampling period $\delta =t_{k+1}-t_k$ , ${\mathit{W}}_k={\tilde{\mathit{W}}}_{t_k}~\mathrm{N}(0,\mathit{Q})$ and ${\mathit{V}}_k={\tilde{\mathit{V}}}_{t_k}~\mathrm{N}(0,\mathit{R})$ . Here, the Runge–Kutta method is adopted to calculate the transition function as

$$F({\mathit{X}}_k)={\mathit{X}}_k+\frac{\delta (F_k^1+2F_k^2+2F_k^3+F_k^4)}{6}$$

(9)

with $F_k^1=f({\mathit{X}}_k)$ , $F_k^2=f({\mathit{X}}_k)+(\delta F_k^1/2)$ , $F_k^3=f({\mathit{X}}_k)+(\delta F_k^2/2)$ and $F_k^4=f({\mathit{X}}_k)+\delta F_k^3$ .

UKF with clipping and uncorrelated conversion

Due to the nonlinearity of model (2), a nonlinear filter should be used to provide estimation for ranges. However, for a nonlinear model, the original measurement cannot be fully used by filters based on linear minimum mean square error estimation, ²² for example, Kalman-based filters. In this article, the uncorrelated conversion (UC), which is proved to be an effective approach to extract additional measurement information, ²² is introduced to further improve the estimation performance. To additionally guarantee the physical constraints (7) for estimation, the constraint handling techniques (see the clipping technique ²³ and quadratic programming (QP) technique ²⁴ ) are required. Here, the clipping technique is used because the nonlinear converted measurement model (induced by UC) prevents the QP technique to be used. In this section, the UKF with appropriate computational complexity and acceptable accuracy is chosen, and the filter process with incorporating clipping and UC is detailed as follows:

• 1. One-step prediction with clipping

Based on unscented transformation (UT), the sigma points of ${\hat{\mathit{X}}}_{k-1}$ are generated as

$${\mathit{\chi }}_{k-1}=[{\mathit{\chi }}_{k-1,j}]=[\begin{array}{cc}{\hat{\mathit{X}}}_{k-1}\hfill & {\hat{\mathit{X}}}_{k-1}\pm \gamma (n)\sqrt{{\mathit{P}}_{k-1}^X}\hfill \end{array}]$$

(10)

where ${\hat{\mathit{X}}}_{k-1}$ and ${\mathit{P}}_{k-1}^X$ are the estimated state vector and corresponding covariance matrix at time instant $t_{k-1}$ , respectively; ${\mathit{\chi }}_{k-1,j}$ ’s, $j\in {\mathbb{I}}_{0:2n}$ are the generated sigma points, with $n=28$ being the dimension of the state vector; and $\gamma (n)=\sqrt{n+\lambda }$ is the suggested scaling factor,^24,25 with $\lambda =(n+\kappa ){\alpha }^2-n$ , $\kappa =3-n$ and $0\le \alpha \le 1$ .

To guarantee the practical constraints (7) in filter, all the sigma points before and after propagation, respectively, are clipped according to the bounds of constraints, that is, for all $j\in {\mathbb{I}}_{0:2n}$

$${\mathit{\chi }}_{k-1,j}^{+}(i)=max\{2\rho ,{\chi }_{k-1,j}(i)\},\forall i\in {\mathbb{I}}_{1:10}$$

(11)

$${\mathit{\chi }}_{k-1,j}^{+}(i)=\mathrm{sat}({\chi }_{k-1,j}(i)),\forall i\in {\mathbb{I}}_{21:28}$$

(12)

and

$${\mathit{\chi }}_{k|k-1,j}^{+}(i)=max\{2\rho ,{\mathbf{\chi }}_{k|k-1,j}(i)\},\forall i\in {\mathbb{I}}_{1:10}$$

(13)

$${\mathit{\chi }}_{k|k-1,j}^{+}(i)=\mathrm{sat}({\mathbf{\chi }}_{k|k-1,j}(i)),\forall i\in {\mathbb{I}}_{21:28}$$

(14)

where ${\mathit{\chi }}_{k|k-1,j}=F({\mathit{\chi }}_{k-1,j}^{+})$ denotes the propagation of the jth sigma point.

Remark 2

Note in physical constraint (7) that each element of state ${\tilde{\mathit{X}}}_t$ is independent of each other. Hence, the clipping technique adopted here can be seemed as a special case of estimate projection, ²⁶ that is, equations (11) and (12) (equations (13) and (14)) are analytical solution of

$${\mathit{\chi }}_{(·)}^{+}=\underset{{\mathit{\chi }}_{(·)}^{+}}{\arg min}\{{({\mathit{\chi }}_{(·)}-{\mathit{\chi }}_{(·)}^{+})}^{\mathrm{T}}\mathbf{\Pi }({\mathit{\chi }}_{(\mathbf{·})}-{\mathit{\chi }}_{(\mathbf{·})}^{+})\}$$

(15)

$$\mathrm{s}.\mathrm{t}.\{\begin{array}{cc}{\mathit{\chi }}_{(·)}^{+}(i)\ge 2\rho \hfill & i\in {\mathbb{I}}_{1:10}\hfill \\ \multicolumn{1}{c}{-1\le {\mathit{\chi }}_{(·)}^{+}(i)\le 1}& i\in {\mathbb{I}}_{21:28}\hfill \end{array}$$

(16)

with any positive definite weighting matrix $\mathbf{\Pi }$ .

After clipping, the one-step prediction of state can be calculated as

$${\hat{\mathit{X}}}_{k|k-1}=\sum _{j=0}^{2n}{\omega }_j^m{\mathit{\chi }}_{k|k-1,j}^{+}$$

(17)

with the weighting scalars

$${\omega }_j^m=\{\begin{array}{cc}\frac{\lambda }{n+\lambda },\hfill & j=0\hfill \\ \multicolumn{1}{c}{\frac{1}{2(n+\lambda )},}& j\in {\mathbb{I}}_{1:2n}\hfill \end{array}$$

(18)

The corresponding covariance matrices can be calculated as

$${\mathit{P}}_{k|k-1}^X=\sum _{j=0}^{2n}{\omega }_j^{\mathrm{c}}({\mathit{\chi }}_{k|k-1,j}^{+}-{\hat{\mathit{X}}}_{k|k-1})({\mathit{\chi }}_{k|k-1,j}^{+}-{\mathit{X}}_{k|k-1}{)}^{\mathrm{T}}$$

(19)

And the one-step prediction of measurement is

$${\hat{\mathit{Z}}}_{k|k-1}={\mathit{H}}_k{\hat{\mathit{X}}}_{k|k-1}$$

(20)

Accordingly, the corresponding covariance matrices can be calculated as

$${\mathit{P}}_{k|k-1}^Z=\sum _{j=0}^{2n}{\omega }_j^{\mathrm{c}}({\hat{\mathit{Z}}}_{k|k-1,j}^{+}-{\hat{\mathit{Z}}}_{k|k-1})({\hat{\mathit{Z}}}_{k|k-1,j}^{\chi }-{\hat{\mathit{Z}}}_{k|k-1}{)}^{\mathrm{T}}$$

(21)

$${\mathit{P}}_{k|k-1}^{\mathrm{XZ}}=\sum _{j=0}^{2n}{\omega }_j^{\mathrm{c}}({\mathit{\chi }}_{k|k-1,j}^{+}-{\hat{\mathit{X}}}_{k|k-1})({\hat{\mathit{Z}}}_{k|k-1,j}^{+}-{\hat{\mathit{Z}}}_{k|k-1}{)}^{\mathrm{T}}$$

(22)

with ${\hat{\mathit{Z}}}_{k|k-1,j}^{+}={\mathit{H}}_k{\mathit{\chi }}_{k|k-1,j}^{+}$ and

$${\omega }_j^{\mathrm{c}}=\{\begin{array}{cc}\frac{\lambda }{n+\lambda }+(1-{\alpha }^2+\beta ),\hfill & j=0\hfill \\ \multicolumn{1}{c}{\frac{1}{2(n+\lambda )},}& j\in {\mathbb{I}}_{1:2n}\hfill \end{array}$$

(23)

• 2. Measurement augmentation with UC

To deeply explore the information of the measurement, the original measurement is converted uncorrelatively by introducing an effective approach and then augmented by the converted measurement. Because the converted measurement is uncorrelated with the original one, the augmentation of the original measurement can provide more useful information and then improve the accuracy of range estimation. The augmented measurement can be generated as follows.

First, the reference distribution-based UC ²² is introduced to uncorrelatively convert the original measurement. Denote ${\mathit{Z}}_k^r$ as a reference random variable with the first two moments ${\overline{\mathit{Z}}}_k^r$ and ${\mathit{P}}_k^{Z^r}$ . To incorporate more information of ${\mathit{Z}}_k$ into ${\mathit{Z}}_k^r$ , the first two moments of ${\mathit{Z}}_r$ are chosen as

$${\overline{\mathit{Z}}}_k^r=E({\mathit{Z}}_k)={\hat{\mathit{Z}}}_{k|k-1}$$

(24)

$${\mathit{P}}_k^{Z^r}=\mathrm{cov}({\mathit{Z}}_k,{\mathit{Z}}_k)={\mathit{P}}_{k|k-1}^Z$$

(25)

Then, the UC of the original measurement can be defined as ²²

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{Y}}_k\stackrel{\Delta }{=}{E}_{{\mathit{Z}}_k^r}[g({\mathit{Z}}_k-{\mathit{Z}}_k^r)]=\int p({\mathit{Z}}_k^r)g({\mathit{Z}}_k-{\mathit{Z}}_k^r)\mathrm{d}{\mathit{Z}}_k^r}\\ \multicolumn{1}{c}{\approx \sum _{j=0}^{2m}{\eta }_j^{m}g({\mathit{Z}}_k-{\mathit{Z}}_{k,j}^r)}\end{array}$$

(26)

where $g(\mathit{z})=e^{-{({\mathit{P}}_{k|k-1}^Z)}^{-1/2}\mathit{z}/b}$ is a nonlinear function with a selected positive scalar b, and $p({\mathit{Z}}_k^r)$ denotes the specific probability density function of ${\mathit{Z}}_k^r$

$$[{\mathit{Z}}_{k,j}^r]=[\begin{array}{cc}{\overline{\mathit{Z}}}_k^r\hfill & {\overline{\mathit{Z}}}_k^r\pm \gamma (m)\sqrt{{\mathit{P}}_k^{Z^r}}\hfill \end{array}]$$

(27)

are the sigma points of ${\mathit{Z}}_k^r$ generated by UT, and

$${\eta }_j^m=\{\begin{array}{cc}\frac{\lambda }{m+\lambda },\hfill & j=0\hfill \\ \multicolumn{1}{c}{\frac{1}{2(m+\lambda )},}& j\in {\mathbb{I}}_{1:2m}\hfill \end{array}$$

(28)

are the associated weighting scalars.

Remark 3

In equation (26), the converted measurement ${\mathit{Y}}_k$ is defined as the expectation of function $g({\mathit{Z}}_k-{\mathit{Z}}_k^r)$ . When the probability density function $p({\mathit{Z}}_k^r)$ is unknown, the expectation can be calculated approximately by different numerical methods, for example, UT, Gaussion hypothesis quadrature rules or cubature rules. Here, the UT approximation is adopted.

According to the definition in equation (26), the mean of the UC ${\mathit{Y}}_k$ can be determined as

$$\begin{array}{c}\multicolumn{1}{c}{{\overline{\mathit{Y}}}_k\stackrel{\Delta }{=}{E}_{{\mathit{Z}}_k}[{\mathit{Y}}_k]=\int p({\mathit{Z}}_k)\int p({\mathit{Z}}_k^r)g({\mathit{Z}}_k-{\mathit{Z}}_k^r)\mathrm{d}{\mathit{Z}}_k^r\mathrm{d}{\mathit{Z}}_k}\\ \multicolumn{1}{c}{\approx \sum _{i=0}^{2m}{\eta }_i^{\mathrm{m}}\underset{{\mathit{Y}}_{k,i}^Z}{\underbrace{\sum _{j=0}^{2m}{\eta }_j^{\mathrm{m}}g({\mathit{Z}}_{k,i}^y-{\mathit{Z}}_{k,j}^r)}}}\end{array}$$

(29)

where $p({\mathit{Z}}_k)$ denotes the specific probability density function of ${\mathit{Z}}_k$

$$[{\mathit{Z}}_{k,i}^y]=[\begin{array}{cc}{\hat{\mathit{Z}}}_{k|k-1}\hfill & {\hat{\mathit{Z}}}_{k|k-1}\pm \gamma (m)\sqrt{{\mathit{P}}_{k|k-1}^Z}\hfill \end{array}]$$

(30)

are the sigma points of ${\mathit{Z}}_k$ generated by UT, ${\eta }_i^{\mathrm{m}}$ ’s are the associated weighting scalars as defined in equation (28), and ${\mathit{Y}}_{k,i}^Z=\sum _{j=0}^{2m}{\eta }_j^{\mathrm{m}}g({\mathit{Z}}_{k,i}^y-{\mathit{Z}}_{k,j}^r)$ denotes the UC of the sigma point ${\mathit{Z}}_{k,i}^y$ . Notice that the equivalence of the first moments between ${\mathit{Z}}_k$ and ${\mathit{Z}}_r$ (i.e. equations (24) and (25)) implies ${\mathit{Z}}_{k,i}^r={\mathit{Z}}_{k,i}^y$ for all $i\in {\mathbb{I}}_{0:2m}$ .

Now, the augmented measurement can be constructed as the combination of the original and uncorrelatively converted measurements, that is

$${\mathit{Z}}_k^a=[\begin{array}{cc}{\mathit{Z}}_k^{\mathrm{T}}\hfill & {\mathit{Y}}_k^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$$

(31)

The mean of the augmented measurement ${\mathit{Z}}_k^a$ can be easily determined as ${\overline{\mathit{Z}}}_k^a=[\begin{array}{cc}{\hat{\mathit{Z}}}_{k|k-1}^{\mathrm{T}}\hfill & {\overline{\mathit{Y}}}_k^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ . The cross-covariance between ${\mathit{X}}_k$ and ${\mathit{Z}}_k^a$ and the covariance of ${\mathit{Z}}_k^a$ can be calculated by UT as

$${\mathit{P}}_{k|k-1}^{X{Z}^a}=[\begin{array}{cc}{\mathit{P}}_{k|k-1}^{\mathrm{XZ}}\hfill & {\mathit{P}}_{k|k-1}^{\mathrm{XY}}\hfill \end{array}]$$

(32)

$${\mathit{P}}_{k|k-1}^{Z^a}=\left[\begin{array}{cc}{\mathit{P}}_{k|k-1}^Z\hfill & {\mathit{P}}_{k|k-1}^{\mathrm{ZY}}\hfill \\ \multicolumn{1}{c}{{({\mathit{P}}_{k|k-1}^{\mathrm{ZY}})}^{\mathrm{T}}}& {\mathit{P}}_{k|k-1}^Y\hfill \end{array}\right]$$

(33)

with

$${\mathit{P}}_{k|k-1}^{\mathrm{XY}}=\sum _{i=0}^{2n}{\omega }_i^{\mathrm{c}}({\mathit{\chi }}_{k|k-1,i}^{+}-{\hat{\mathit{X}}}_{k|k-1})({\mathit{Y}}_{k|k-1,i}^{+}-{\hat{\mathit{Y}}}_{k|k-1}{)}^{\mathrm{T}}$$

$${\mathit{P}}_{k|k-1}^{\mathrm{ZY}}=\sum _{i=0}^{2m}{\eta }_i^c({\mathit{Z}}_{k,i}^y-{\hat{\mathit{Z}}}_{k|k-1})({\mathit{Y}}_{k,i}^Z-{\overline{\mathit{Y}}}_k{)}^{\mathrm{T}}$$

$${\mathit{P}}_{k|k-1}^Y=\sum _{i=0}^{2m}{\eta }_i^{\mathrm{c}}({\mathit{Y}}_{k,i}^Z-{\overline{\mathit{Y}}}_k)({\mathit{Y}}_{k,i}^Z-{\overline{\mathit{Y}}}_k{)}^{\mathrm{T}}$$

where ${\mathit{Y}}_{k|k-1,i}^{+}=\sum _{j=0}^{2m}{\eta }_j^{\mathrm{m}}g({\hat{\mathit{Z}}}_{k|k-1,i}^{+}-{\mathit{Z}}_{k,j}^r)$ denotes the UC of ${\hat{\mathit{Z}}}_{k|k-1,i}^{+}={\mathit{H}}_k{\mathit{\chi }}_{k|k-1,i}^{+}$ , ${\hat{\mathit{Y}}}_{k|k-1}=\sum _{i=0}^{2n}{\omega }_i^{\mathrm{m}}{\mathit{Y}}_{k|k-1,i}^{+}$ , and

$${\eta }_i^c=\{\begin{array}{cc}\frac{\lambda }{m+\lambda }+(1-{\alpha }^2+\beta )\hfill & i=0\hfill \\ \multicolumn{1}{c}{\frac{1}{2(m+\lambda )}}& i\in {\mathbb{I}}_{1:2m}\hfill \end{array}$$

(34)

• 3. Estimation update

With the augmented measurement ${\mathit{Z}}_k^a$ , the one-step prediction of state ${\hat{\mathit{X}}}_{k|k-1}$ can be updated by linear minimum mean square error estimator as

$${\hat{\mathit{X}}}_k={\hat{\mathit{X}}}_{k|k-1}+{\mathit{P}}_{k|k-1}^{X{Z}^a}({\mathit{P}}_{k|k-1}^{Z^a}{)}^{-1}({\mathit{Z}}_k^a-{\overline{\mathit{Z}}}_k^a)$$

(35)

and the associated covariance matrix is calculated as

$${\mathit{P}}_k^X={\mathit{P}}_{k|k-1}^X-{\mathit{P}}_{k|k-1}^{X{Z}^a}({\mathit{P}}_k^{Z^a}{)}^{-1}({\mathit{P}}_{k|k-1}^{X{Z}^a}{)}^{\mathrm{T}}$$

(36)

The whole process of UKF with clipping and UC is summarized as Algorithm 1.

Algorithm 1. UKF with clipping and UC.
Input: the updated measurement range ${\mathit{Z}}_k$ at time $t_k$ ;
the estimated state vector ${\hat{\mathit{X}}}_{k-1}$ at time $t_{k-1}$ ;
the covariance matrix ${\mathit{P}}_{k-1}^X$ at time $t_{k-1}$ .
Output: the estimated state vector ${\hat{\mathit{X}}}_k$ at time $t_k$ ;
the covariance matrix ${\mathit{P}}_k^X$ at time $t_k$ .
1: Generate the sigma points ${\mathit{\chi }}_{k-1,j}$ of state vector ${\hat{\mathit{X}}}_{k-1}$ as equation (10).
2: Clip the sigma points ${\mathit{\chi }}_{k-1,j}$ to the bounds of constraints as equations (11) and (12).
3: Calculate propogated sigma points ${\mathit{\chi }}_{k|k-1,j}$ .
4: Clip the propogated sigma points ${\mathit{\chi }}_{k|k-1,j}$ as equations (13) and (14).
5: Using ${\mathit{\chi }}_{k|k-1,j}^{+}$ to calculate the one-step prediction of state vector ${\hat{\mathit{X}}}_{k|k-1}$ in equation (17) and measurement vector ${\hat{\mathit{Z}}}_{k|k-1}$ as equation (20). And the corresponding covariance matrix ${\mathit{P}}_{k|k-1}^X$ in equation (19), ${\mathit{P}}_{k|k-1}^Z$ in equation (21), and ${\mathit{P}}_{k|k-1}^{\mathrm{XZ}}$ in equation (22) can also be calculated.
6: Generate uncorrelated conversion of the original measurement ${\mathit{Y}}_k$ as equation (26).
7: Calculate the mean of the uncorrelated conversion ${\overline{\mathit{Y}}}_k$ as equation (29).
8: Calculate the corresponding covariance matrix ${\mathit{P}}_{k|k-1}^{X{Z}^a}$ as equation (32) and ${\mathit{P}}_{k|k-1}^{Z^a}$ as equation (33).
9: Calculate the estimation of ${\hat{\mathit{X}}}_k$ as equation (35) and ${\mathit{P}}_k^X$ as equation (36).

Clustering and merging

According to estimated ranges in the basic localization unit ${\mathbb{N}}_p^{\mathrm{unit}}$ , the distance matrix ${\mathit{D}}_{p,k}^{\mathrm{unit}}$ of ${\mathbb{N}}_p^{\mathrm{unit}}$ at time instant $t_k$ is constructed as

$${\mathit{D}}_{p,k}^{\mathrm{unit}}=\left[\begin{array}{ccccc}0\hfill & {\hat{\mathit{X}}}_k(1)\hfill & {\hat{\mathit{X}}}_k(2)\hfill & {\hat{\mathit{X}}}_k(3)\hfill & {\hat{\mathit{X}}}_k(4)\hfill \\ \multicolumn{1}{c}{{\hat{\mathit{X}}}_k(1)}& 0\hfill & {\hat{\mathit{X}}}_k(5)\hfill & {\hat{\mathit{X}}}_k(6)\hfill & {\hat{\mathit{X}}}_k(7)\hfill \\ \multicolumn{1}{c}{{\hat{\mathit{X}}}_k(2)}& {\hat{\mathit{X}}}_k(5)\hfill & 0\hfill & {\hat{\mathit{X}}}_k(8)\hfill & {\hat{\mathit{X}}}_k(9)\hfill \\ \multicolumn{1}{c}{{\hat{\mathit{X}}}_k(3)}& {\hat{\mathit{X}}}_k(6)\hfill & {\hat{\mathit{X}}}_k(8)\hfill & 0\hfill & {\hat{\mathit{X}}}_k(10)\hfill \\ \multicolumn{1}{c}{{\hat{\mathit{X}}}_k(4)}& {\hat{\mathit{X}}}_k(7)\hfill & {\hat{\mathit{X}}}_k(9)\hfill & {\hat{\mathit{X}}}_k(10)\hfill & 0\hfill \end{array}\right]$$

(37)

where ${\hat{\mathit{X}}}_k(i),i\in \{1,2,\dots ,28\}$ is the state variable at time instant $t_k$ . According to ${\mathit{D}}_{p,k}^{\mathrm{unit}}$ , the distance matrix ${\mathit{D}}_k$ of the global sensor network $\mathbb{N}$ can be constructed with entrances

$${\mathit{D}}_k(i,j)=\{\begin{array}{cc}\frac{1}{m_{\mathrm{ij}}}\sum _{\begin{array}{c}p\in {\mathbb{I}}_{\mathrm{I}:M},\forall i,j\in {\mathbb{N}}_p^{\mathrm{unit}}\hfill \end{array}}{\hat{d}}_{p,k}^{\mathrm{ij}},\hfill & i\ne j\mathrm{and}{m}_{\mathrm{ij}}\ne 0\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{otherwise}\hfill \end{array}$$

(38)

where ${\hat{d}}_{p,k}^{\mathrm{ij}}$ is the estimated range between sensors i and j in ${\mathit{D}}_{p,k}^{\mathrm{unit}}$ , and $m_{\mathrm{ij}}$ is the number of basic localization units that include i and j.

With the estimated ranges in $\mathbb{N}$ , a distributed MDS-MAP localization algorithm is adopted in this article to localize the mobile sensor network. MDS-MAP method is first proposed by Shang et al. ²⁷ of Columbia University to localize sensors in networks. Directly using the MDS-MAP method in the localization of the whole sensor network $\mathbb{N}$ with incomplete measurement will lead to errors. An effective approach is first dividing $\mathbb{N}$ into O local sensor networks with complete measurement, then using the MDS-MAP method to calculate their accurate relative positions and finally merging these local sensor networks to obtain the relative positions of $\mathbb{N}$ . The absolute positions can be further acquired using the anchor sensors ${\mathbb{N}}_a$ . The specific procedure of localization stage is as follows:

• 1. Clustering and local-localization

According to the distance matrix ${\mathit{D}}_k$ , the entrances of the connection matrix ${\mathit{L}}_k$ can be calculated as

$${\mathit{L}}_k(i,j)=\{\begin{array}{cc}1,\hfill & {\mathit{D}}_k(i,j)\ne 0\hfill \\ \multicolumn{1}{c}{1,}& j=i\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{otherwise}\hfill \end{array}$$

(39)

According to the connection matrix ${\mathit{L}}_k$ , the global network can be clustered into clusters. Considering requirement in local-localization with the MDS-MAP method, all-pair estimated ranges between sensors in each cluster should be obtained. To this end, the process of clustering and local-localization is executed as Algorithm 2.

Algorithm 2. Clustering and local-localization.
Input: the distance matrix ${\mathit{D}}_k$ ;
the connection matrix ${\mathit{L}}_k$ .
Output: the local sensor networks ${\mathbb{C}}_s,s\in {\mathbb{I}}_{0:O-1}$ ;
the relative positions ${\mathit{F}}_s$ of sensor networks ${\mathbb{C}}_s$ .
1: for all node $i\in \mathbb{N}$ do
2:  for all node $j\in \mathbb{N}$ do
3:   if $L_k(i,j)=1$ then
4:    Add node j to the set ${\mathbb{S}}_i=\{j\in \mathbb{N}|L_k(i,j)=1\}$
5:   end if
6:  end for
7: end for
8: for all node $i\in \mathbb{N}$ do
9:  Update ${\mathbb{S}}_i$ as ${\mathbb{S}}_i={\mathbb{S}}_i\cap {\mathbb{S}}_j$ , $\forall j\in {\mathbb{S}}_i$ .
10: end for
11: Search for the set ${\mathbb{S}}_i,i\in \mathbb{N}$ of largest cardinality, denoted as ${\mathbb{S}}_i^{*}$ . All nodes in ${\mathbb{S}}_i^{*}$ construct a reference sensor network ${\mathbb{C}}_0$ .
12: Calculate the relative positions of ${\mathbb{C}}_0$ using the MDS-MAP method.
13: Set $s=0$ , ${\mathbb{N}}_0=\mathbb{N}-{\mathbb{S}}_i^{*}$ .
14: while ${\mathbb{N}}_0\ne \varnothing$ do
15:  Update $s=s+1$ and ${\mathbb{N}}_0={\mathbb{N}}_0\backslash {\mathbb{S}}_i^{*}$ .
16:  Search for the set ${\mathbb{S}}_i,i\in \mathbb{N}$ of largest cardinality, which can be denote as ${\mathbb{S}}_i^{*}$ . All nodes in ${\mathbb{S}}_i^{*}$ construct a local sensor network ${\mathbb{C}}_s$ .
17:  Calculate the relative positions of ${\mathbb{C}}_s$ using the MDS-MAP method.
18: end while

For each local sensor network ${\mathbb{C}}_s$ , the MDS-MAP method in Algorithm 2 is briefly stated as follows. According to the distance matrix ${\mathit{D}}_k$ , construct the inner product matrix ${\mathit{B}}_s$ with entrances

$${\mathit{B}}_s(s_i,s_j)=\frac{1}{2}({\mathit{D}}_k^2(i,1)+{\mathit{D}}_k^2(j,1)-{\mathit{D}}_k^2(i,j))$$

(40)

where $s_i$ and $s_j$ are the indices of sensors i and j in ${\mathbb{C}}_s$ , respectively. By decomposing the inner product matrix as ${\mathit{B}}_s=\mathit{UE}{\mathit{U}}^{\mathrm{T}}$ , the relative position matrix ${\hat{\mathit{F}}}_k^s$ of the local sensor network ${\mathbb{C}}_s$ can be determined as

$$\begin{array}{c}{\hat{\mathit{F}}}_k^s={[{\hat{\mathit{x}}}_i^s,{\hat{\mathit{x}}}_j^s,\dots ]}^{\mathrm{T}}\hfill \end{array}$$

(41)

$$=\mathrm{diag}\{\sqrt{{\lambda }_1},\sqrt{{\lambda }_2},\sqrt{{\lambda }_3}\}[{\mathit{u}}_1,{\mathit{u}}_2,{\mathit{u}}_3]$$

(42)

where ${\hat{\mathit{x}}}_i^s=[x_i^s,y_i^s,z_i^s{]}^{\mathrm{T}}$ is the relative position of sensor i in ${\mathbb{C}}_s$ ; ${\lambda }_1$ , ${\lambda }_2$ and ${\lambda }_3$ are the largest three eigenvalues of E ; ${\mathit{u}}_1$ , ${\mathit{u}}_2$ and ${\mathit{u}}_3$ are the corresponding eigenvectors in U .

• 2. Merging of local sensor networks

According to the relative positions of the common sensors in different local sensor networks, the global map can be constructed using the rotation and translation method. Namely, merge other local sensor networks into the reference sensor network (i.e. the relative positions of the reference network ${\mathbb{C}}_0$ ) one by one.

Assume that common nodes of the local reference network ${\mathbb{C}}_0$ and another local sensor network ${\mathbb{C}}_s$ are denoted by $\{a,b,c,d,\dots \}$ . According to equation (43), each sensor $i\in {\mathbb{C}}_s\backslash ({\mathbb{C}}_s\cap {\mathbb{C}}_0)$ can be merged into the reference network ${\mathbb{C}}_0$ through rotation and translation method as

$${\hat{\mathit{x}}}_i^0={\mathbf{\Omega }}_{0s}{\hat{\mathit{x}}}_i^s+{\mathit{t}}_{0s}$$

(43)

with the rotation matrix

$${\mathbf{\Omega }}_{0s}=\frac{({\hat{\mathit{F}}}_k^0-1_{c\times 3}{({\overline{\mathit{F}}}_k^0)}^{\mathrm{T}}){({\hat{\mathit{F}}}_k^s-1_{c\times 3}{({\overline{\mathit{F}}}_k^s)}^{\mathrm{T}})}^{\mathrm{T}}}{({\hat{\mathit{F}}}_k^s-1_{c\times 3}{({\overline{\mathit{F}}}_k^s)}^{\mathrm{T}}){({\hat{\mathit{F}}}_k^s-1_{c\times 3}{({\overline{\mathit{F}}}_k^s)}^{\mathrm{T}})}^{\mathrm{T}}}$$

and the translation vector ${\mathit{t}}_{0s}=({\overline{\mathit{F}}}_k^0{)}^{\mathrm{T}}-{\mathbf{\Omega }}_{0s}({\overline{\mathit{F}}}_k^s{)}^{\mathrm{T}}$ , where ${\overline{\mathit{F}}}_k^s=1_{1\times c}{\hat{\mathit{F}}}_k^s/c$ and ${\overline{\mathit{F}}}_k^0=1_{1\times c}{\hat{\mathit{F}}}_k^0/c$ are the central positions of ${\mathbb{C}}_s$ and ${\mathbb{C}}_0$ , respectively, and c is the number of common sensors.

• 3. Global absolute positions generated

Using the absolute position ${\mathit{F}}_k^a$ of anchor sensors in ${\mathbb{N}}_a$ and equation (43), the relative positions of the global sensor network ${\widehat{\mathit{F}}}_k^0=[{\widehat{\mathit{x}}}_1^0,{\widehat{\mathit{x}}}_2^0,\dots ,{\widehat{\mathit{x}}}_N^0]$ can be converted to the absolute positions ${\hat{\mathit{F}}}_k=[{\hat{\mathit{x}}}_1,{\hat{\mathit{x}}}_2,\dots ,{\hat{\mathit{x}}}_N]$ of the global sensor network.

Mobile sensor network localization algorithm with incomplete measurement

For a mobile sensor network $\mathbb{N}$ of N sensors, the localization algorithm for the whole mobile sensor network is given in Algorithm 3.

Algorithm 3. Mobile sensor network localization.
Input: the incomplete measurement range ${\mathit{Z}}_k,k=1,2,\dots$ ;
the absolute positions ${\mathit{F}}_k^a,k=1,2,\dots$ of anchor nodes ${\mathbb{N}}_a$ .
Output: the positions of the sensor network ${\hat{\mathit{F}}}_k,k=1,2,\dots$ .
1: Group $\mathbb{N}$ into M basic localization units ${\mathbb{N}}_p^{\mathrm{unit}},p\in {\mathbb{I}}_{1:M}$ according to the measurement information.
2: Set $k=1$ .
3: for all ${\mathbb{N}}_p^{\mathrm{unit}},p\in {\mathbb{I}}_{1:M}$ do
4:  Calculated ${\hat{\mathit{X}}}_k$ and ${\mathit{P}}_k^X$ as Algorithm 1.
5:  Calculated ${\mathit{D}}_{p,k}^{\mathrm{unit}}$ as equation (37).
6: end for
7: Calculated ${\mathit{D}}_k$ and ${\mathit{L}}_k$ as equations (38) and (39).
8: Clustering into O clusters ${\mathbb{C}}_s,s\in \mathbb{O}={\mathbb{I}}_{0:O-1}$ and calculate the relative positions ${\mathit{F}}_k^s,s\in \mathbb{O}$ as Algorithm 2.
9: Set $\mathrm{mergetime}=1$
10: while $\mathrm{mergetime}\le O-1$ do
11:  for all ${\mathbb{C}}_s,s\in \mathbb{O}$ do
12:   if the number of common sensors of ${\mathbb{C}}_s$ and ${\mathbb{C}}_0$ is no less than 4 then
13:    Merging the ${\mathbb{C}}_s$ into ${\mathbb{C}}_0$ by equation (43).
14:    Update ${\mathbb{C}}_0={\mathbb{C}}_0\cup {\mathbb{C}}_s$ .
15:    Update $\mathbb{O}=\mathbb{O}\backslash \{s\}$ .
16:   end if
17:  end for
18: end while
19: Using ${\mathit{F}}_k^a,k=1,2,\dots$ to transform the relative positions ${\hat{\mathit{F}}}_k^0$ into absolute ones ${\hat{\mathit{F}}}_k$ by equation (43).

Simulation

To verify the proposed localization algorithm for mobile sensor network, several algorithms are compared in this section:

1. Sequential Monte Carlo localization (SMCL) algorithm. ¹³

2. MDS-MAP algorithm with EKF (MDS-EKF). ⁸

3. MDS-MAP algorithm with UKF (MDS-UKF). ⁸

4. The proposed localization algorithm with QP-UKF. ²⁴

5. The proposed localization algorithm without UC.

In the comparison, the ranging error ratio is set separately at 1% and 10%, and we take the root mean square error (RMSE)

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}\left[\sum _{i=1}^N{({\mathit{x}}_i-{\hat{\mathit{x}}}_i)}^2\right]}$$

(44)

to describe the localization error, and stress

$$\mathrm{Stress}=\sqrt{\frac{\sum _{i=1}^N\sum _{i=1}^N{({d_k^{\mathrm{ij}}}^2-{{\hat{d}}_k^{\mathrm{ij}}}^2)}^2}{\sum _{i=1}^N\sum _{j=1}^N{d_k^{\mathrm{ij}}}^4}}$$

(45)

to describe the topology error of mobile sensor network. The smaller the RMSE and stress are calculated, the higher the localization accuracy is achieved.

Considering the sensor network composed of 10 mobile nodes in 3D space, and each node is free to move at an average speed of 3 m/s. Assuming that ranges among nodes in sensor network can be completely measured. The localization results of different localization algorithms under different range error ratios (denoted as e) are shown in Figures 3–12.

Figure 3.

RMSE of localization algorithm with 1% range error.

Figure 4.

Stress of localization algorithm with 1% range error.

Figure 5.

RMSE of localization algorithm with 10% range error.

Figure 6.

Stress of localization algorithm with 10% range error.

Figure 7.

RMSE of the localization algorithm with different types of UKF.

Figure 8.

Stress of the localization algorithm with different types of UKF.

Figure 9.

The boxplot for RMSE of different algorithms with e = 1%.

Figure 10.

The boxplot for stress of different algorithms with e = 1%.

Figure 11.

The boxplot for RMSE of different algorithms with e = 10%.

Figure 12.

The boxplot for stress of different algorithms with e = 10%.

In the simulation, the ranges are gradually increasing due to free movement of nodes in sensor network. The pre-set relative range error ratio leads to the increase in localization error (i.e. RMSE), as shown in Figures 3, 5 and 7. In Figures 3–6, both the average RMSE (0.37 m under 1% range error and 0.78 m under 10% range error) and the average stress (0.004 under 1% range error and 0.031 under 10% range error) of the proposed algorithm are much smaller than those of MDS-EKF, MDS-UKF and SMCL algorithms. This verifies the higher localization accuracy of the proposed localization algorithm.

In the framework of the proposed localization algorithm, the UKF with different techniques (i.e. UKF with QP, ²⁴ UKF with clipping and UKF with both clipping and UC) is compared to demonstrate the effectiveness of introducing clipping and UC. The performance comparison is displayed in Figures 7 and 8. Simulation results show that the average of RMSE is 0.61 and 1.35 m (4.17 and 5.98 m) for UKF with QP and clipping, respectively, under $e=1\%$ ( $e=10\%$ ); and the average of stress is 0.014 and 0.313 (0.077 and 0.218) for UKF with QP and clipping, respectively, under $e=1\%$ ( $e=10\%$ ). All of them are larger than UKF with both clipping and UC, which verifies the better performance of introducing clipping and UC.

In the form of boxplot, the average values, 25th and 75th percentiles, and the maximum and minimum of RMSE and stress in different algorithms are illustrated in Figures 9–12. In these boxplots, the central mark denotes the average value, the edges of the box denote the 25th and 75th percentiles, the upper and lower lines denote the maximum and minimum, and the separate points are the abnormal data. The shape of the box clearly shows the distribution of the data. It is easy to notice that the average of RMSE and stress of the proposed localization algorithm are smaller than those of other localization algorithms and so are the distributions. This demonstrates the accuracy and stability of the proposed algorithm.

The comparison for time cost of these algorithms is provided in Table 1. It is easily seen that the proposed algorithm possesses the highest localization accuracy with an acceptable increase in time cost. Compared with the algorithm without UC, the proposed algorithm is more accurate, with a slight increase in time cost, which verifies the effectiveness of introducing UC.

Table 1.

The comparison for average time cost of different algorithms.

Algorithm	Average RMSE with 1% range error (m)	Average time cost before convergence (s)	Average time cost after convergence (s)	Average time cost (s)
SMCL	16.7226	0.5635	0.0393	0.1441
MDS-EKF	2.5344	0.0166	0.0033	0.0034
MDS-UKF	1.5132	0.0304	0.0061	0.0062
QP-UKF	0.9872	0.5458	0.2622	0.2913
Algorithm without UC	1.3511	0.0157	0.0136	0.0136
The proposed localization algorithm	0.3741	0.0466	0.0324	0.0316

RMSE: root mean square error; SMCL: Sequential Monte Carlo localization; MDS-EKF: multidimensional scaling extended Kalman filter; MDS-UKF: multidimensional scaling unscented Kalman filter; QP-UKF: quadratic programming; UC: uncorrelated conversion.

To verify the proposed node localization performance with incomplete measurement, we take the connectivity parameter $C=l/l_0$ to describe the degree of connectivity in the whole wireless sensor network, with l denoting the number of the measurable ranges and $l_0$ denoting the number of the ranges in the whole WSN. The localization results of proposed algorithm with different connectivities (ranging from $53.3\%$ to $100\%$ ) are shown in Figures 13 and 14.

Figure 13.

The average RMSE of proposed algorithm with different connectivities.

Figure 14.

The average stress of proposed algorithm with different connectivities.

It can be easily seen from Figures 13 and 14 that the lowest connectivity of the proposed algorithm can reach $53.3\%$ . Note that when the sensor network connectivity is less than $66.7\%$ , the average RMSE is about $20-30\mathrm{m}$ and the average stress is about $0.4-0.5$ under different range errors. This is because the numbers of ranges in some basic localization units are less than 10, and the effectiveness of the proposed algorithm is guaranteed by clustering and merging methods and the relative motion model (see Remark 1). With the increase in connectivity, the accuracy of localization is gradually enhanced.

Conclusion

In this article, the localization problem of mobile sensor network with incomplete measurement is investigated, and a localization algorithm is proposed. The effectiveness of the proposed localization algorithm is achieved by the relative motion model with physical constraints, the UKF with clipping and UC, and the clustering and merging methods. Comparison with other existing algorithms demonstrates that the proposed algorithm can achieve highest localization accuracy. For mobile sensor network with different connectivities, the proposed localization algorithm is still valid, and its localization accuracy is gradually improved with the increase in connectivity.