Fast 3D Node Localization in Multipath for UWB Wireless Sensor Networks Using Modified Propagator Method

Abstract

A fast three-dimensional (3D) node localization algorithm in multipath for ultra-wideband (UWB) wireless sensor networks is developed. The algorithm employs a modified propagator method (MPM) for time delay estimation and then uses a marriage algorithm of 3D Chan and Taylor for range-based multilateral localization and node position computation. In the proposed algorithm, the traditional propagator method (PM) for direction-of-arrival (DOA) estimation is extended to frequency-domain time-of-arrival (TOA) estimation in multipath, which can effectively measure the distance between an unknown sensor node and an anchor node. MPM algorithm requires neither spectral searching nor covariance matrix estimation and its eigenvalue decomposition, which reduces the computational complexity. The marriage location algorithm enhances the robustness and accuracy of node localization. The simulations validate the effectiveness of the proposed algorithm in locating multiple unknown nodes of UWB wireless sensor networks in 3D space.

1. Introduction

Ranging and localization of unknown sensor nodes in wireless sensor networks (WSNs) [1] have drawn considerable attention in many aspects such as environmental monitoring, health tracking, smart home, machine-to-machine (M2M), and body area networks (BAN). In actual environment, the sensor nodes are commonly placed in three-dimensional (3D) terrains, such as workshops, forests, and oceans. Recently, several range-based and range-free 3D node localization algorithms in WSNs have been proposed [2–6]. However, some problems should be investigated in node design and implementation, such as large computational amount, slow executing, and inaccurate localization in multipath and noise. Since the energy and power carried by node equipment are limited, effective methods with low complexity, low power consumption, and good robustness for 3D node positioning are of great significance.

This paper investigates fast 3D node localization problem based on time delay measurement in multipath in ultra-wideband (UWB) wireless sensor networks. As we know, time delay estimation problem has been studied with a variety of super-resolution subspace techniques [7–11], such as multiple signal classification (MUSIC), total least square-estimation of signal parameters via rotational invariance technique (TLS-ESPRIT), and matrix pencil (MP). Compared with the correlator-based methods, subspace-based algorithms can increase time resolution even if the time delay is smaller than a pulse width. Unfortunately, these techniques increase the complexity of WSN implementation. Specifically, in MUSIC algorithm, spectral searching is needed through all the space, which increases its computational time. ESPRIT-like algorithms require the covariance matrix estimation of the observed data using a large number of snapshots and then performing eigenvalue decomposition (EVD) or singular value decomposition (SVD) of it. In [12], we developed a unitary matrix pencil (UMP) based delay estimation algorithm for 3D localization in UWB wireless sensor networks, which can reduce the complexity load. However, its estimation accuracy is lower than many subspace-based methods.

The propagator method (PM), developed in [13], is a traditional subspace method previously used in the direction-of-arrival (DOA) estimation. In [14], we extend it to the time-of-arrival (TOA) estimation. PM algorithm only requires linear operations, and thus it avoids covariance matrix estimation and its EVD or SVD which are the main computational burden in subspace methods. However, similar to MUSIC algorithm, a spectral peak searching through all the space is needed in traditional PM algorithm.

In this paper, we develop a fast range-based 3D node localization algorithm in multipath for UWB wireless sensor networks. The algorithm employs a modified propagator method (MPM) for frequency-domain time delay estimation and then uses a marriage algorithm of 3D Chan and Taylor for range-based multilateral localization and node position computation. The MPM algorithm requires neither spectral searching nor covariance matrix estimation and its eigenvalue decomposition, which can reduce the computational load. Furthermore, the marriage localization algorithm can enhance the robustness and accuracy of node localization in comparison with the previous algorithm in [15]. The simulations validate the effectiveness of the proposed algorithm in locating multiple unknown nodes of UWB wireless sensor networks in 3D space.

This paper is organized as follows. In Section 2, the UWB received signal model for multipath time delay estimation is given. MPM-based multipath delay estimation algorithm is proposed in Section 3. Based on it, the multilateral localization and node position computation using a marriage of 3D Chan and Taylor algorithms are presented in Section 4. In Section 5, the simulation results are given to verify the performance of the proposed method. Finally, a conclusion is drawn in Section 6.

2. Signal Model

Assume that an ultra-wideband pulse is transmitted from an unknown node to an anchor node through L paths. During the qth snapshot, $q=\mathrm{1},\dots ,Q$ , the received signal can be expressed as

$$\begin{array}{cc}& y^{\left(q\right)}\left(t\right)=\sum _{l=\mathrm{1}}^L{\beta }_l^{\left(q\right)}p\left(t-{\tau }_{\text{TOA}}-\Delta {\tau }_l\right)+w^{\left(q\right)}\left(t\right),\end{array}$$

(1)

where $p(t)$ denotes a UWB signal having basic Gaussian pulse waveform $p(t)=A_p{\text{e}}^{-\mathrm{2}\pi t^{\mathrm{2}}/T_p^{\mathrm{2}}}$ , with $A_p$ and $T_p$ being the amplitude and pulse width of the UWB signal, respectively. $w^{(q)}(t)$ denotes the complex additive Gaussian white noise. ${\tau }_{\text{TOA}}$ denotes the TOA of the unknown node, and $\Delta {\tau }_l$ and ${\beta }_l^{(q)}$ represent the relative delay and time-varying complex fading amplitude of the lth path, respectively: $\Delta {\tau }_{\mathrm{1}}=\mathrm{0}$ . Denote ${\tau }_l={\tau }_{\text{TOA}}+\Delta {\tau }_l$ as the total propagation time delay of the lth path. Then, the received signal in (1) can be rewritten as

$$\begin{array}{cc}& y^{\left(q\right)}\left(t\right)=\sum _{l=\mathrm{1}}^L{\beta }_l^{\left(q\right)}p\left(t-{\tau }_l\right)+w^{\left(q\right)}\left(t\right).\end{array}$$

(2)

The channel impulse response of UWB signal in multipath is a sum of L components shifted from the corresponding time delays. Upon discrete Fourier transform and matched filtering based on (2), the frequency-domain representation of the identified channel signal can be written as

$$\begin{array}{l}{H}^{\left(q\right)}\left(k\right)=\sum _{l=\mathrm{1}}^L{\beta }_l^{\left(q\right)}{e}^{-j(\mathrm{2}\pi k{\tau }_l/KT)}+V^{\left(q\right)}\left(k\right)\\ =\sum _{l=\mathrm{1}}^L{\beta }_l^{\left(q\right)}{z}_l^k+V^{\left(q\right)}\left(k\right),\end{array}$$

(3)

$k=\mathrm{0,1},\dots ,K-\mathrm{1}$ , with K denoting the number of discrete Fourier transforms and T being the sampling period. $V^{(q)}(k)$ denotes the frequency-domain noise, ${\tau }_l$ is the estimated delay parameter, and $z_l=e^{-j(\mathrm{2}\pi {\tau }_l/KT)}$ .

Collecting the data of all the Q snapshots from (3) yields an observed data matrix:

$$\begin{array}{cc}& \mathbf{H}=\mathbf{Z}\left(\tau \right)\mathbf{B}+\mathbf{V}\mathbf{,}\end{array}$$

(4)

where $\mathbf{H}\in {ℂ}^{K\hspace{0.17em}\times \hspace{0.17em}Q}$ and $\mathbf{V}\in {ℂ}^{K\hspace{0.17em}\times \hspace{0.17em}Q}$ denote the observed matrix and noise matrix, respectively. $\mathbf{Z}(\tau )$ denotes the time delay matrix:

$$\begin{array}{cc}& \mathbf{Z}\left(\tau \right)={\left[\begin{array}{ccc}\mathrm{1}& \cdots & 1\\ z_{\mathrm{1}}& \cdots & z_L\\ ⋮& \ddots & ⋮\\ z_{\mathrm{1}}^{K-\mathrm{1}}& \cdots & z_L^{K-\mathrm{1}}\end{array}\right]}_{K\hspace{0.17em}\times \hspace{0.17em}L},\end{array}$$

(5)

and $\mathbf{B}\in {ℂ}^{L\hspace{0.17em}\times \hspace{0.17em}Q}$ denotes the fading amplitude matrix:

$$\begin{array}{cc}& \mathbf{B}\hspace{0.17em}\mathbf{=}\left[\begin{array}{ccc}{\beta }_{\mathrm{1}}^{\left(\mathrm{1}\right)}& \cdots & {\beta }_{\mathrm{1}}^{\left(Q\right)}\\ ⋮& \cdots & ⋮\\ {\beta }_L^{\left(\mathrm{1}\right)}& \cdots & {\beta }_L^{\left(Q\right)}\end{array}\right].\end{array}$$

(6)

The problem interest is to estimate the parameter ${\tau }_l$ based on the observed matrix in (4), $l=\mathrm{1},\dots ,L$ .

3. Multipath Delay Estimation Using Modified Propagator Method

3.1. Review of Traditional Propagator Method

Let us first review the traditional propagator method. The propagator method is previously presented in DOA estimation [13]. We introduce its application in TOA estimation [14].

From (5), it is noted that $\mathbf{Z}(\tau )$ , as a Vandermonde matrix, is of full column rank with L linearly independent rows. The other $K-L$ rows can be expressed as a linear combination of these L rows. Therefore, $\mathbf{Z}(\tau )$ can be partitioned into

$$\begin{array}{cc}& \mathbf{Z}\left(\tau \right)=\left[\begin{array}{c}{\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)\\ {\mathbf{Z}}_{\mathrm{2}}\left(\tau \right)\end{array}\right]\begin{array}{c}\}L\\ \}K-L,\end{array}\end{array}$$

(7)

where ${\mathbf{Z}}_{\mathrm{1}}(\tau )\in {ℂ}^{L\hspace{0.17em}\times \hspace{0.17em}L}$ and ${\mathbf{Z}}_{\mathrm{2}}(\tau )\in {ℂ}^{(K-L)\hspace{0.17em}\times \hspace{0.17em}L}$ . The propagator is the unique linear operator P of ${ℂ}^{K-L}$ into ${ℂ}^L$ such that

$$\begin{array}{cc}& {\mathbf{P}}^H{\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)={\mathbf{Z}}_{\mathrm{2}}\left(\tau \right)\end{array}$$

(8)

or

$$\begin{array}{cc}& \left[{\mathbf{P}}^H,-{\mathbf{I}}_{K-L}\right]\mathbf{Z}\left(\tau \right)={\mathbf{Q}}^H\mathbf{Z}\left(\tau \right)=\mathbf{0}\mathbf{,}\end{array}$$

(9)

where H denotes the Hermitian transpose. ${\mathbf{I}}_{K-L}\in {ℝ}^{(K-L)\hspace{0.17em}\times \hspace{0.17em}(K-L)}$ and $\mathbf{0}\in {ℝ}^{(K-L)\hspace{0.17em}\times \hspace{0.17em}L}$ denote the identity matrix and null matrix, respectively, and the matrix Q is given by

$$\begin{array}{cc}& \mathbf{Q}\triangleq \left[\begin{array}{c}\mathbf{P}\\ -{\mathbf{I}}_{K-L}\end{array}\right].\end{array}$$

(10)

Equation (9) reveals that $\mathbf{Z}(\tau )$ is orthogonal to the columns of Q. Therefore, the following TOA estimator can find the peaks of the spectrum:

$$\begin{array}{cc}& F_{\text{PM}}\left(\tau \right)=\frac{\mathrm{1}}{{\mathbf{Z}}^H\left(\tau \right)\mathbf{Q}{\mathbf{Q}}^H\mathbf{Z}\left(\tau \right)}.\end{array}$$

(11)

In (11), Q can be effectively estimated from the linear operations of the observed data. We similarly partition the observed matrix H in (4) into

$$\begin{array}{cc}& \begin{array}{cc}\mathbf{H}=\left[\begin{array}{c}{\mathbf{H}}_{\mathrm{1}}\\ {\mathbf{H}}_{\mathrm{2}}\end{array}\right]& \begin{array}{c}\}L\\ \}K-L,\end{array}\end{array}\end{array}$$

(12)

where ${\mathbf{H}}_{\mathrm{1}}\in {ℂ}^{L\hspace{0.17em}\times \hspace{0.17em}Q}$ and ${\mathbf{H}}_{\mathrm{2}}\in {ℂ}^{(K-L)\hspace{0.17em}\times \hspace{0.17em}Q}$ . In the noise-free case, from (4), (8), and (12), it yields

$$\begin{array}{cc}& {\mathbf{H}}_{\mathrm{2}}={\mathbf{P}}^H{\mathbf{H}}_{\mathrm{1}}.\end{array}$$

(13)

In the presence of noise, a least squares (LS) estimate of the propagator operator P is

$$\begin{array}{cc}& \widehat{\mathbf{P}}={\left({\mathbf{H}}_{\mathrm{2}}{{\mathbf{H}}_{\mathrm{1}}}^{+}\right)}^H,\end{array}$$

(14)

where ${{\mathbf{{\rm H}}}_{\mathrm{1}}}^{+}$ is the Moore-Penrosepseudo inverse of ${\mathbf{{\rm H}}}_{\mathrm{1}}$ , ${{\mathbf{{\rm H}}}_{\mathrm{1}}}^{+}\hspace{0.17em}\mathbf{=}\hspace{0.17em}\hspace{0.17em}\mathbf{(}{{\mathbf{{\rm H}}}_{\mathrm{1}}}^H{\mathbf{{\rm H}}}_{\mathrm{1}}{\mathbf{)}}^{\mathbf{-}\mathrm{1}}{{\mathbf{{\rm H}}}_{\mathrm{1}}}^H$ . Thus the estimate of Q can be calculated from (10).

3.2. The Modified Propagator Method

The traditional propagator method can reduce computational load compared with the common subspace techniques such as MUSIC and ESPRIT since it uses linear operations on the observed data, avoiding the EVD or SVD of covariance matrix. However, it reveals from (11) that long time is needed for spectral peak searching through all the scope of time delays. To avoid it, we combine the idea of shift invariance with propagator method in the paper, presenting a modified propagator method (MPM) for multipath delay estimation.

Let ${\mathbf{H}}_{\mathrm{1}}$ and ${\mathbf{H}}_{\mathrm{2}}$ be two $(K-\mathrm{1})\hspace{0.17em}\times \hspace{0.17em}Q$ submatrices of the matrix H consisting of the first and last $K-\mathrm{1}$ rows of H; then,

$$\begin{array}{cc}& {\mathbf{H}}_{\mathrm{1}}={\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)\mathbf{B}+{\mathbf{V}}_{\mathrm{1}},\hspace{1em}\hspace{1em}{\mathbf{H}}_{\mathrm{2}}={\mathbf{Z}}_{\mathrm{2}}\left(\tau \right)\mathbf{B}+{\mathbf{V}}_{\mathrm{2}},\end{array}$$

(15)

where ${\mathbf{Z}}_{\mathrm{1}}(\tau )$ and ${\mathbf{Z}}_{\mathrm{2}}(\tau )$ are the first and last $K-\mathrm{1}$ rows of the matrix $\mathbf{Z}(\tau )$ , respectively. Thus, we have

$$\begin{array}{cc}& {\mathbf{Z}}_{\mathrm{2}}\left(\tau \right)={\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)\mathbf{\Phi },\end{array}$$

(16)

where $\mathbf{\Phi }=\mathrm{diag}\{e^{-j(\mathrm{2}\pi {\tau }_{\mathrm{1}}/KT)},\dots ,e^{-j(\mathrm{2}\pi {\tau }_L/KT)}\}$ is a diagonal matrix, which contains the estimated time delay parameters ${\tau }_{\mathrm{1}},\dots ,{\tau }_L$ .

Compose the two matrices ${\mathbf{H}}_{\mathrm{1}}$ and ${\mathbf{H}}_{\mathrm{2}}$ to form a $\mathrm{2}(K-\mathrm{1})\hspace{0.17em}\times \hspace{0.17em}Q$ matrix:

$$\begin{array}{cc}& \mathbf{X}=\left[\begin{array}{c}{\mathbf{H}}_{\mathrm{1}}\\ {\mathbf{H}}_{\mathrm{2}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)\\ {\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)\mathbf{\Phi }\end{array}\right]\mathbf{B}+\left[\begin{array}{c}{\mathbf{V}}_{\mathrm{1}}\\ {\mathbf{V}}_{\mathrm{2}}\end{array}\right]=\mathbf{D}\left(\tau \right)\mathbf{B}+\mathbf{N},\end{array}$$

(17)

where $\mathbf{D}(\tau )$ is a $\mathrm{2}(K-\mathrm{1})\hspace{0.17em}\times \hspace{0.17em}L$ delay matrix and N is the noise matrix upon composition.

From the idea of propagator method, if $\mathbf{D}(\tau )$ is partitioned into

$$\begin{array}{cc}& \begin{array}{cc}\mathbf{D}\left(\tau \right)=\left[\begin{array}{c}{\mathbf{D}}_{\mathrm{1}}\left(\tau \right)\\ {\mathbf{D}}_{\mathrm{2}}\left(\tau \right)\end{array}\right]& \begin{array}{c}\}L\\ \}2\left(K-L\right)-L,\end{array}\end{array}\end{array}$$

(18)

then there exists a unique propagator operator P satisfying

$$\begin{array}{cc}& {\mathbf{P}}^H{\mathbf{D}}_{\mathrm{1}}\left(\tau \right)={\mathbf{D}}_{\mathrm{2}}\left(\tau \right).\end{array}$$

(19)

Similarly, the matrix X can be partitioned into two matrices ${\mathbf{X}}_{\mathrm{1}}$ and ${\mathbf{X}}_{\mathrm{2}}$ in the same way as $\mathbf{D}(\tau )$ . If the noise matrix N is ignored, then

$$\begin{array}{cc}& {\mathbf{X}}_{\mathrm{1}}={\mathbf{D}}_{\mathrm{1}}\left(\tau \right)\mathbf{B},\hspace{1em}\hspace{1em}{\mathbf{X}}_{\mathrm{2}}={\mathbf{D}}_{\mathrm{2}}\left(\tau \right)\mathbf{B}.\end{array}$$

(20)

Premultiply ${\mathbf{P}}^H$ in ${\mathbf{X}}_{\mathrm{1}}$ yielding

$$\begin{array}{cc}& {\mathbf{P}}^H{\mathbf{X}}_{\mathrm{1}}={\mathbf{P}}^H{\mathbf{D}}_{\mathrm{1}}\left(\tau \right)\mathbf{B}={\mathbf{D}}_{\mathrm{2}}\left(\tau \right)\mathbf{B}={\mathbf{X}}_{\mathrm{2}}.\end{array}$$

(21)

In the noisy environment, P can be achieved by the LS solution; that is,

$$\begin{array}{cc}& \mathbf{P}={\left({\mathbf{X}}_{\mathrm{2}}{\mathbf{X}}_{\mathrm{1}}^H{\left({\mathbf{X}}_{\mathrm{1}}{\mathbf{X}}_{\mathrm{1}}^H\right)}^{-\mathrm{1}}\right)}^H.\end{array}$$

(22)

Let $\tilde{\mathbf{P}}=\left[\begin{array}{c}{\mathit{I}}_L\\ {\mathbf{P}}^H\end{array}\right]$ , where ${\mathit{I}}_L$ is a $L\times L$ identity matrix; then

$$\begin{array}{cc}& \tilde{\mathbf{P}}{\mathbf{D}}_{\mathrm{1}}\left(\tau \right)=\left[\begin{array}{c}{\mathit{I}}_L\\ {\mathbf{P}}^H\end{array}\right]{\mathbf{D}}_{\mathrm{1}}\left(\tau \right)=\left[\begin{array}{c}{\mathbf{D}}_{\mathrm{1}}\left(\tau \right)\\ {\mathbf{D}}_{\mathrm{2}}\left(\tau \right)\end{array}\right]=\mathbf{D}\left(\tau \right).\end{array}$$

(23)

From (23), we observe that $\tilde{\mathbf{P}}$ spans signal subspace. Evenly dividing $\tilde{\mathbf{P}}$ into two $(K-\mathrm{1})\times L$ matrices ${\tilde{\mathbf{P}}}_{\mathrm{1}}$ and ${\tilde{\mathbf{P}}}_{\mathrm{2}}$ and based on (17), we have

$$\begin{array}{cc}& \left[\begin{array}{c}{\tilde{\mathbf{P}}}_{\mathrm{1}}\\ {\tilde{\mathbf{P}}}_{\mathrm{2}}\end{array}\right]{\mathbf{D}}_{\mathrm{1}}\left(\tau \right)=\left[\begin{array}{c}{\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)\\ {\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)\mathbf{\Phi }\end{array}\right].\end{array}$$

(24)

Therefore,

$$\begin{array}{cc}& {\mathbf{Z}}_{\mathrm{1}}\left(\tau \right)\mathbf{\Phi }=\mathbf{\Psi }{\mathbf{Z}}_{\mathrm{1}}\left(\tau \right),\end{array}$$

(25)

where $\mathbf{\Psi }={\tilde{\mathbf{P}}}_{\mathrm{2}}{({{\tilde{\mathbf{P}}}_{\mathrm{1}}}^H{\tilde{\mathbf{P}}}_{\mathrm{1}})}^{-\mathrm{1}}{{\tilde{\mathbf{P}}}_{\mathrm{1}}}^H$ . It is shown from (25) that Φ is the eigenvalue matrix of Ψ. Therefore, the time delays of L paths can be calculated from L main diagonal elements of Φ; that is,

$$\begin{array}{cc}& {\widehat{\tau }}_l=\frac{KT}{\mathrm{2}\pi }\text{angle}\left({\lambda }_l\right),\end{array}$$

(26)

where ${\lambda }_l$ is the lth eigenvalue of Φ, for $l=\mathrm{1},\dots ,L$ .

The procedure of MPM-based time delay estimation algorithm can be summarized as follows.

• (i)

Divide the observed matrix H into ${\mathbf{H}}_{\mathrm{1}}$ and ${\mathbf{H}}_{\mathrm{2}}$ consisting of the first and last $K-\mathrm{1}$ rows of H, respectively.

(ii)

Compose ${\mathbf{H}}_{\mathrm{1}}$ and ${\mathbf{H}}_{\mathrm{2}}$ to form a $\mathrm{2}(K-\mathrm{1})\times Q$ matrix X. Partition X into $\mathbf{X}=[{\mathbf{X}}_{\mathrm{1}}^T,{\mathbf{X}}_{\mathrm{2}}^T{]}^T$ such that ${\mathbf{X}}_{\mathrm{1}}\in {ℂ}^{L\hspace{0.17em}\times \hspace{0.17em}Q}$ and ${\mathbf{X}}_{\mathrm{2}}\in {ℂ}^{(\mathrm{2}(K-\mathrm{1})-L)\hspace{0.17em}\times \hspace{0.17em}Q}$ .

(iii)

Obtain the propagator $\mathbf{P}={({\mathbf{X}}_{\mathrm{2}}{\mathbf{X}}_{\mathrm{1}}^H{({\mathbf{X}}_{\mathrm{1}}{\mathbf{X}}_{\mathrm{1}}^H)}^{-\mathrm{1}})}^H$ . Let $\tilde{\mathbf{P}}=\left[\begin{array}{c}{\mathit{I}}_L\\ {\mathbf{P}}^H\end{array}\right]$ , and evenly divide $\tilde{\mathbf{P}}$ into two $(K-\mathrm{1})\hspace{0.17em}\times \hspace{0.17em}L$ matrices ${\tilde{\mathbf{P}}}_{\mathrm{1}}$ and ${\tilde{\mathbf{P}}}_{\mathrm{2}}$ .

(iv)

Calculate the eigenvalue matrix Φ of the matrix $\mathbf{\Psi }={\tilde{\mathbf{P}}}_{\mathrm{2}}{({{\tilde{\mathbf{P}}}_{\mathrm{1}}}^H{\tilde{\mathbf{P}}}_{\mathrm{1}})}^{-\mathrm{1}}{{\tilde{\mathbf{P}}}_{\mathrm{1}}}^H$ . Then, ${\widehat{\tau }}_l$ can be obtained from the diagonal elements of Φ.

4. Node Position Computation Algorithm

4.1. Multilateral Localization

From the results of MPM-based multipath delay estimation given in (26), the ${\tau }_{\text{TOA}}$ and distance r from an unknown node to anchor nodes can be determined. In the following procedure of 3D node position computation, we consider multilateral localization other than trilateral localization to improve the accuracy. The number of anchor nodes is no less than five in our algorithm. A geometry relationship of the cross point of spheres with five anchor nodes (A, B, C, D, E) is shown in Figure 1.

Figure 1

The geometric relationship of 3D multilateral localization with five anchor nodes.

Assume that $(x,y,z)$ denotes the coordinate of an unknown node. $(X_{\mathrm{0}},Y_{\mathrm{0}},Z_{\mathrm{0}})$ represents the coordinate of a main anchor node, and $(X_i,Y_i,Z_i)$ represents the coordinate of the other anchor nodes, $i=\mathrm{1,2},\dots ,M-\mathrm{1}$ , with M denoting the number of anchor nodes. Let t be the time of the transmitted UWB signal by the unknown node and let $t_i$ be the time of the received first-path signal by the ith anchor node. Then, the distance between the unknown node and the ith anchor node is $r_i=c(t_i-t)=c{{\tau }_{\text{TOA}}}_i$ . The equation of 3D node position computation can be written by

$$\begin{array}{l}{r}_{\mathrm{0}}^{\mathrm{2}}={\left(x-X_{\mathrm{0}}\right)}^2+{\left(y-Y_{\mathrm{0}}\right)}^2+{\left(z-Z_{\mathrm{0}}\right)}^2,\\ r_i^{\mathrm{2}}={\left(x-X_i\right)}^2+{\left(y-Y_i\right)}^2+{\left(z-Z_i\right)}^2,\\ \Delta r_i=r_i-r_{\mathrm{0}}=c{\tau }_{\text{TDOA}i},\end{array}$$

(27)

where c denotes the velocity of wave, $\Delta r_i$ and ${\tau }_{\text{TDOA}i}$ denote the distance difference and time-difference-of-arrival (TDOA) from the unknown node to the ith anchor node and to the main anchor node, and ${\tau }_{\text{TDOA}i}={{\tau }_{\text{TOA}}}_i-{{\tau }_{\text{TOA}}}_{\mathrm{0}}$ , $i=\mathrm{1,2},\dots ,M-\mathrm{1}$ .

Chan algorithm [16] and Taylor algorithm [17] are two traditional methods to solve nonlinear positioning equations for radio location in a 2D space. Chan algorithm has low computational complexity and high accuracy in high signal-noise-ratio (SNR) and Gaussian noise environment. However, with the increase of TOA measurement error, the performance of Chan algorithm in low SNR environment will degrade. Taylor algorithm has higher accuracy and good robustness in noise. However, the performance of Taylor algorithm is highly dependent on the initial estimate of iterative computation. An improper initial estimate will lead to no convergence.

In the paper, we present a marriage algorithm of 3D Chan and Taylor location to calculate the 3D position of unknown node. Figure 2 shows the procedure of locating an unknown node with $M=5$ anchor nodes.

Figure 2

Procedure of locating an unknown node with five anchor nodes.

4.2. $\text{3D}$ Chan Algorithm

The traditional Chan algorithm can be extended from 2D space to 3D space. From the positioning equation in (27), we have

$$\begin{array}{cc}& \left(X_{\mathrm{0}}-X_i\right)x+\left(Y_{\mathrm{0}}-Y_i\right)y+\left(Z_{\mathrm{0}}-Z_i\right)z=k_i+r_{\mathrm{0}}\Delta r_i,\end{array}$$

(28)

where $k_i=(\mathrm{1}/\mathrm{2})[\Delta r_i^{\mathrm{2}}-(X_i^{\mathrm{2}}+Y_i^{\mathrm{2}}+Z_i^{\mathrm{2}})+(X_{\mathrm{0}}^{\mathrm{2}}+Y_{\mathrm{0}}^{\mathrm{2}}+Z_{\mathrm{0}}^{\mathrm{2}})]\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}$ for $i=\mathrm{1,2},\dots ,M-\mathrm{1}$ . The coordinate $(x,y,z)$ of unknown node can be obtained using the weighted LS (WLS) solution [18] for in (28), and thus,

$$\begin{array}{l}\left[\begin{array}{l}x\\ y\\ z\end{array}\right]=-{\left[\begin{array}{lll}{X}_{\mathrm{1,0}}& Y_{\mathrm{1,0}}& Z_{\mathrm{1,0}}\\ X_{\mathrm{2,0}}& Y_{\mathrm{2,0}}& Z_{\mathrm{2,0}}\\ ⋮& ⋮& ⋮\\ X_{M-\mathrm{1,0}}& Y_{M-\mathrm{1,0}}\mathrm{\hspace{0.17em}\hspace{0.17em}}& Z_{M-\mathrm{1,0}}\end{array}\right]}^{+}\\ \times \left\{\frac{\mathrm{1}}{\mathrm{2}}\left[\begin{array}{l}\Delta r_{\mathrm{1}}^{\mathrm{2}}-K_{\mathrm{1}}+K_{\mathrm{0}}\\ \Delta r_{\mathrm{2}}^{\mathrm{2}}-K_{\mathrm{2}}+K_{\mathrm{0}}\\ ⋮\\ \Delta r_{M-\mathrm{1}}^{\mathrm{2}}-K_{M-\mathrm{1}}+K_{\mathrm{0}}\end{array}\right]+\left[\begin{array}{l}\Delta r_{\mathrm{1}}\\ \Delta r_{\mathrm{2}}\\ ⋮\\ \Delta r_{M-\mathrm{1}}\end{array}\right]r_{\mathrm{0}}\right\},\end{array}$$

(29)

where $X_{i,\mathrm{0}}=X_i-X_{\mathrm{0}}$ , $Y_{i,\mathrm{0}}=Y_i-Y_{\mathrm{0}}$ , $Z_{i,\mathrm{0}}=Z_i-Z_{\mathrm{0}}$ , $K_i=X_i^{\mathrm{2}}\hspace{0.17em}+\hspace{0.17em}{Y}_i^{\mathrm{2}}\hspace{0.17em}+\hspace{0.17em}{Z}_i^{\mathrm{2}}$ , $K_{\mathrm{0}}=X_{\mathrm{0}}^{\mathrm{2}}\hspace{0.17em}+\hspace{0.17em}{Y}_{\mathrm{0}}^{\mathrm{2}}\hspace{0.17em}+\hspace{0.17em}{Z}_{\mathrm{0}}^{\mathrm{2}}$ , and “+” denotes pseudoinverse of a matrix.

In (29), since $r_{\mathrm{0}}$ contains unknown coordinate $(x,y,z)$ , substituting (29) into (27) yields a second-order equation with respect to $r_{\mathrm{0}}$ . Obtaining the WLS solution of $r_{\mathrm{0}}$ and substituting its positive root into (29), we can attain an accurate position estimate $(\widehat{x},\widehat{y},\widehat{z})$ of unknown node.

4.3. $\text{3D}$ Taylor Algorithm

Similarly, Taylor algorithm can also be extended to 3D space. Starting with an initial estimate $(x^{(\mathrm{0})},y^{(\mathrm{0})},z^{(\mathrm{0})})$ , the 3D Taylor algorithm iteratively calculates the position error Δ using WLS:

$$\begin{array}{cc}& \mathbf{\Delta }=\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]={\left({\mathbf{G}}^T{\mathbf{Q}}^{-\mathrm{1}}\mathbf{G}\right)}^{-1}{\mathbf{G}}^T{\mathbf{Q}}^{-\mathrm{1}}\mathbf{h}\mathbf{,}\end{array}$$

(30)

where Q is the covariance matrix of TDOA measurement. Consider the following:

$$\begin{array}{cc}& \mathbf{G}=\left[\begin{array}{ccc}\frac{\left(X_{\mathrm{0}}-x\right)}{r_{\mathrm{0}}}-\frac{\left(X_{\mathrm{1}}-x\right)}{r_{\mathrm{1}}}& \frac{\left(Y_{\mathrm{0}}-y\right)}{r_{\mathrm{0}}}-\frac{\left(Y_{\mathrm{1}}-y\right)}{r_{\mathrm{1}}}& \frac{\left(Z_{\mathrm{0}}-z\right)}{r_{\mathrm{0}}}-\frac{\left(Z_{\mathrm{1}}-z\right)}{r_{\mathrm{1}}}\\ \frac{\left(X_{\mathrm{0}}-x\right)}{r_{\mathrm{0}}}-\frac{\left(X_{\mathrm{2}}-x\right)}{r_{\mathrm{2}}}& \frac{\left(Y_{\mathrm{0}}-y\right)}{r_{\mathrm{0}}}-\frac{\left(Y_{\mathrm{2}}-y\right)}{r_{\mathrm{2}}}& \frac{\left(Z_{\mathrm{0}}-z\right)}{r_{\mathrm{0}}}-\frac{\left(Z_{\mathrm{2}}-z\right)}{r_{\mathrm{2}}}\\ ⋮& ⋮& ⋮\\ \frac{\left(Z_{\mathrm{0}}-z\right)}{r_{\mathrm{0}}}-\frac{\left(Z_{\mathrm{2}}-z\right)}{r_{\mathrm{2}}}& \frac{\left(Y_{\mathrm{0}}-y\right)}{r_{\mathrm{0}}}-\frac{\left(Y_{M-\mathrm{1}}-y\right)}{r_{M-\mathrm{1}}}& \frac{\left(Z_{\mathrm{0}}-z\right)}{r_{\mathrm{0}}}-\frac{\left(Z_{M-\mathrm{1}}-z\right)}{r_{M-\mathrm{1}}}\end{array}\right],\\ \mathbf{h}=\left[\begin{array}{c}\Delta r_{\mathrm{1}}-(r_{\mathrm{1}}-r_{\mathrm{0}})\\ \Delta r_{\mathrm{2}}-(r_{\mathrm{2}}-r_{\mathrm{0}})\\ ⋮\\ \Delta r_{M-\mathrm{1}}-(r_{M-\mathrm{1}}-r_{\mathrm{0}})\end{array}\right].\end{array}$$

(31)

4.4. Marriage of 3D Chan/Taylor Location

Based on the two algorithms mentioned above, a marriage algorithm of 3D Chan and Taylor location is presented to calculate the 3D position of the unknown node. The procedure is as follows. In the first step, 3D Chan algorithm is employed to calculate the estimated coordinate of an unknown node based on (29), denoted by $(\tilde{x},\tilde{y},\tilde{z})$ . In the following step, taking the result of 3D Chan algorithm as an initial value of Taylor series expansion, 3D Taylor algorithm is used to iteratively calculate the location error Δ in (30). The iterative procedure has ended until Δ is smaller than a threshold. Then, an accurate estimate of the node position can be calculated by

$$\begin{array}{cc}& \left[\begin{array}{c}\widehat{x}\\ \widehat{y}\\ \widehat{z}\end{array}\right]={\left[\begin{array}{c}\tilde{x}\\ \tilde{y}\\ \tilde{z}\end{array}\right]}_{\text{Chan}}+{\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]}_{\text{Taylor}}.\end{array}$$

(32)

Therefore, the 3D location method inherits the advantages of Chan and Taylor algorithms, having small computational complexity, high accuracy, and good robustness in noise.

5. Simulations

In the simulations, the performance improvement achieved by the MPM is presented. Then, node localization accuracy using MPM-based multipath delay estimation and the marriage of 3D Chan/Taylor location is verified.

5.1. MPM-Based Multipath Delay Estimation Results

The transmitted UWB signal ( $T_p=\mathrm{0.1}\hspace{0.17em}\text{ns}$ , $A_p=\mathrm{1}$ ) and the received five-path superposition signal are shown in Figure 3.

Figure 3

The transmitted UWB signal and the received multipath superposition signal.

Table 1 illustrates the results of time delay estimation using the proposed MPM algorithm. $\text{SNR}=0\hspace{0.17em}\text{dB}$ . $L=5$ . $Q=100$ . $K=100$ . $T=0.04$  ns. The simulations confirm high-resolution performance of the proposed MPM-based delay estimation algorithm in multipath environment, even at lower SNRs.

Table 1

MPM-based delay estimation results when SNR = 0 dB.

Path delays	True (ns)	Estimate (ns)
Path 1	1.60	1.594
Path 2	1.65	1.643
Path 3	1.70	1.695
Path 4	1.75	1.752
Path 5	1.80	1.797

Figure 4 presents the comparison of root mean square error (RMSE) and run time between the MPM algorithm and TLS-ESPRIT algorithm [10]. There are $L=3$ paths with time delays 1.6 ns, 1.7 ns, and 1.8 ns. $K=40$ . $\hspace{0.17em}\hspace{0.17em}T=0.1$  ns, with 100 Monte Carlo trials. In Figure 4(b), we evaluate the computational complexity using TIC and TOC instructions in MATLAB. Figure 4 demonstrates that MPM algorithm and TLS-ESPRIT algorithm have similar RMSE performance under medium and high SNR. However, MPM algorithm has less computational complexity compared with TLS-ESPRIT algorithm. The computation load of TLS-ESPRIT mainly relies on estimating the covariance matrix ${\mathbf{R}}_{\mathbf{X}\mathbf{X}}\in {ℂ}^{\mathrm{2}(K-\mathrm{1})\hspace{0.17em}\times \hspace{0.17em}\mathrm{2}(K-\mathrm{1})}$ and performing eigenvalue decomposition for three times, while MPM algorithm avoids the covariance matrix estimation, whose computation load mainly relies on the calculation of propagator P and an eigenvalue decomposition.

Figure 4

Performance comparison between MPM and TLS-ESPRIT.

5.2. $\text{3D}$ Position Computation Results

In Figure 5, we assume that $M=7$ anchor nodes are arranged in a $\mathrm{2}\hspace{0.17em}\text{m}\times \mathrm{2}\hspace{0.17em}\text{m}\times \mathrm{2}\hspace{0.17em}\text{m}$ space. Among them, a main anchor node is in the center of the cube, the coordinate of which is (0, 0, 0). The coordinates of the other six anchor nodes are (−1, 0, 0), (0, 1, 0), (0, 0, 1), (−1, 0, 0), (0, −1, 0), and (0, 0, −1), respectively. Multiple unknown nodes are randomly generated within the space. Assume that UWB signal is sent from an unknown node to an anchor node through $L=2$ paths. SNR = 0 dB. At first, we can obtain ${\tau }_{\text{TOA0}}~{\tau }_{\text{TOA6}}$ using MPM-based estimation algorithm, in which $Q=100$ . $K=20$ . $T=1$  ns. Then, the positions of multiple unknown nodes can be calculated using the marriage algorithm of 3D Chan/Taylor location.

Figure 5

Node position computation result.

The results of 3D position computation of 10 or 100 unknown nodes are demonstrated in Figures 5(a) and 5(b), respectively. The results show that the coordinates of multiple unknown nodes can be effectively estimated with high accuracy, even with low density of anchor nodes.

5.3. Accuracy of the Proposed Algorithm

We investigate the localization accuracy and robustness of the proposed algorithm by simulation. The node position error can be calculated by

$$\begin{array}{l}\text{Position}\hspace{0.17em}\hspace{0.17em}\text{Error}\\ \hspace{0.5em}=\sqrt{\frac{\mathrm{1}}{\mathrm{3}N}\sum _{n=\mathrm{1}}^N\left({\left({\widehat{x}}_n-x\right)}^2+{\left({\widehat{y}}_n-y\right)}^2+{\left({\widehat{z}}_n-z\right)}^2\right)},\end{array}$$

(33)

where N is the number of unknown nodes and ( ${\widehat{x}}_n,{\widehat{y}}_n,{\widehat{z}}_n$ ) is the coordinate estimate of ( $x_n,y_n,z_n$ ) for the nth unknown node. 10 unknown nodes are randomly generated in the space. When SNR = 0 dB, the performance of node position error versus the number of anchor nodes is shown in Figure 6, in which the number of anchor nodes varies from 5 to 8. The simulation result from Figure 6 indicates that the accuracy is promoted with the increase of the number of anchor nodes.

Figure 6

Position error versus the number of anchor nodes (SNR = 0 dB).

Furthermore, the performance of the proposed algorithm is compared with the matrix pencil based localization algorithm in [12] and the previously proposed localization algorithm in [15], as shown in Figure 7. The simulation results reveal that the performance of the proposed algorithm outperforms the other algorithms. By employing the MPM-based time delay estimation algorithm, the TOA measurement accuracy of our method can be greatly improved compared with matrix pencil based TOA measurement algorithm in [12]. In addition, by the use of the marriage of 3D Chan/Taylor location, the proposed method has better robustness in low SNR compared with the previous algorithm in [15].

Figure 7

The comparison of localization accuracy (7 anchor nodes).

6. Conclusion

In this paper, we propose a fast 3D node localization algorithm in multipath for UWB wireless sensor networks. It utilizes MPM-based multipath time delay estimation and the marriage of 3D Chan/Taylor location. Compared with the traditional range-based methods, our method has greatly reduced the computational complexity and enhanced the robustness and localization accuracy in multipath and noise environment. The proposed algorithm has key superiorities in fast executing, high accuracy in multipath, and low SNR, which is applicable for node localization in 3D space.