Robust TDOA Localization Algorithm for Asynchronous Wireless Sensor Networks

Abstract

This paper presents a tag localization algorithm based on the time-difference-of-arrival (TDOA) of mobile tag signal for asynchronous wireless sensor network (WSN) with N anchors (nodes with known locations) and a large number of mobile tags. To obtain time synchronization, all anchors broadcast signals periodically; relative clock offsets and skews of anchor pairs are estimated by the least-square (LS) method using the times-of-arrival (TOAs) of broadcast signals at anchors. When a tag transmits signal, the TOA of tag signal at each anchor is stamped and errors in original TDOAs of tag signal due to relative clock offsets and skews of anchor pairs are eliminated. Based on Gaussian noise model, maximum likelihood estimation (MLE) for the tag position is obtained. Performance issues are addressed by evaluating the Cramér-Rao lower bound of synchronization and localization algorithms. Since the tag can be located via a single transmission, least power consumption of tag is required, and large number of tags can be served in WSN. The proposed algorithm is simple and effective, with performance close to that of synchronous TDOA algorithm.

1. Introduction

With the advances of wireless communications and microelectronics technologies, wireless sensor networks (WSNs) have received more attention over the recent decade due to their potential application to a wide variety of diverse areas [1], such as environmental monitoring, space exploration, military applications, target tracking, and healthcare.

One of the most important issues for WSN is that the location of each sensor node must be determined. However, it is expensive and impractical to equip all sensors in the network with a global positioning system (GPS) receiver [2]. To obtain accurate location estimates of sensor nodes, range-based localization algorithms are more favorable than range-free ones [3–5]. In general, the range-based algorithms always follow two steps [6, 7]: they first measure some metrics bearing location information, the so-called ranging or bearing, and second estimate the positions based on those metrics, the so-called location information fusion. There are mainly four metrics: time-of-arrival (TOA) or time-of-flight (TOF) [8], time-difference-of-arrival (TDOA) [5, 6], angle-of-arrival (AOA) [9], and received signal strength (RSS) [10]. The ranging methods using RSS can be implemented by energy detectors, but they can only achieve a coarse resolution. Antenna arrays are required for AOA-based methods, which encumbers their popularity. On the other hand, high accuracy and potentially low cost implementation make TOA or TDOA based on ultrawideband system (UWB) and broadband spread spectrum system a promising ranging method [11, 12].

Since TOA and TDOA measurements are time-based, clock synchronization is essential for high accuracy localization and low cost implementation. Many protocols have been proposed to resolve the synchronization problem of WSN, such as the reference broadcast synchronization (RBS) protocol [13], the timing-sync protocol for sensor networks (TPSN) [14], and flooding time synchronization protocol (FTSP) [15]. On the other hand, clock synchronization can also be handled by estimating the relative clock offsets and skews among sensor nodes using maximum likelihood estimator (MLE) [3–5, 7]; the theoretical performance limits are evaluated by Cramér-Rao lower bound (CRLB).

Due to the close relationship between synchronization and localization problems, joint estimation of the relative clock skew, the relative clock offset, and the positions of source nodes are proposed in [4, 5] for network with asynchronous anchors. The broadcasting property is exploited to jointly estimate clock parameters and source position using asymmetrical time-stamping and a passive listening protocol in [5, 16]. Although these algorithms achieve good performance, complex protocols, such as two-way ranging protocol, are necessary. Moreover, frequent transmissions of source node may cause more energy dissipation and limit the number of tags in the service area of WSN.

In this paper, we propose a robust TDOA localization algorithm for asynchronous wireless sensor networks. In WSN, sensor nodes are divided into two categories, anchor and tag. Tags are mobile sensor nodes; anchor nodes are used as reference points to localize tag nodes and assumed to have known positions and good power supplies. There are many limitations for tag nodes, such as size, cost, and energy dissipation which is concerned intensively for practical implementation. In order to synchronize anchor clocks, all anchors broadcast signal periodically; TOAs of broadcasting signal at other anchors are stamped according to their own local clocks, respectively. A localization server on the backbone network collects all these timestamps to estimate the relative clock offsets and skews of anchor pairs. When anchors receive a broadcasting signal from a tag, the TOAs of tag signal at anchors are also stamped and delivered to the localization server. The error in original TDOA of tag signal due to relative clock offsets and skews among anchors can be easily canceled by a compensation operation. Based on a Gaussian measurement noise model, maximum likelihood estimation (MLE) for tag position is obtained. Performance issues are addressed by evaluating the Cramér-Rao lower bound (CRLB) of synchronization and localization algorithm. Only anchor-anchor synchronization is required in this paper and tag node is localized via a signal transmission, so power consumption of tag can be minimized and the number of tags to be served can be maximized.

This paper is organized as follows. The sensor network model for localization problem considered in the paper is given in Section 2. Estimation of relative clock offsets and skews of anchor pairs and performance analysis are presented in Section 3. The localization of tags is studied in Section 4, and the statistical performance analysis of the localization algorithm is also presented in this section. Section 5 presents some numerical results and Section 6 concludes the paper.

2. Sensor Network Model for Localization

Wireless sensor network localization is the process by which sensor nodes determine their location. In simple terms, localization is a mechanism for discovering spatial relationships between objects. The various approaches taken in literature to solve this localization problem differ in the assumptions they make about their respective network and sensor capabilities. A detailed, but not exhaustive, list of assumptions made include assumptions about network infrastructure, device hardware, signal propagation models, timing and energy requirements, operational environment, time synchronization, communication costs, error requirements, and node mobility.

In this paper, sensor nodes are divided into two categories: anchor node and tag node. Anchor nodes are used as reference points to locate tag nodes; they are always connected to network infrastructure and well powered. To reduce network complexity and costs, the quantity of anchor is limited and positions of anchor nodes are planned in advance to guarantee the coverage of service area. The tags are sensor nodes to be located and must be designed for low cost and low power consumption and operation for a long lifetime. To further reduce the complexity and power consumption of tag, the processing capability is restricted for practical implementation. For example, a WSN based patient management consists of limited anchors which serve for large quantity of tags carried by patients.

The radio sensor network model considered in the paper is depicted in Figure 1; all anchors receive broadcasting signals from other anchors or tags, stamp TOAs according to their own clocks, and route the TOA information to the localization server where synchronization and localization algorithms are implemented. The only work of tag is to transmit signals with certain period which is determined by working condition and requirement, measuring TOA and extra processing function are not needed in tag node. For tags with low mobility, a large value of period can be chosen to prolong the working time.

Figure 1

Radio WSN model for tag localization (anchor nodes are represented by red rectangles and tags are represented by green circles).

In Figure 1, all sensor nodes broadcast signals periodically with respect to their own local clocks. The message of broadcasting signal does not include a timestamp generated by the sender, nor is it important exactly when it is sent. As a matter of fact, the broadcast of tag does not even need to be a dedicated packet for localization. Almost any extant broadcast can be used to tag localization, for instance, ARP packets, or the broadcast control traffic in wireless networks (e.g., RTS/CTS exchanges or route discovery packets).

Consider a radio WSN with N anchors at known, fixed positions ${\mathbf{x}}_i$ for $i=\mathrm{1},\dots ,N$ and a tag at an unknown position ${\mathbf{x}}_{\mathrm{0}}$ (since tags are located using TOAs arriving at anchors independently, the proximity of tags cannot be taken into consideration; for simplicity, only one tag is considered in the paper). The position vectors, ${\mathbf{x}}_i$ , $i=\mathrm{0,1},\dots ,N$ , are all 2D or 3D. Throughout the paper, the subscripts $i=\mathrm{1},\dots ,N$ refer to anchor i and subscript 0 refers to the tag. Each anchor has a quartz clock that operates at the same nominal frequency $f_{\mathrm{0}}$ with unknown clock skew ${\epsilon }_i$ . The frequency $f_i$ of anchor i is [6]

$$\begin{array}{cc}& f_i=f_{\mathrm{0}}\left(\mathrm{1}+{\epsilon }_i\right).\end{array}$$

(1)

Considering random clock drift due to frequency fluctuation, the frequency $f_i$ may change slightly over times in practice. In this paper, the following assumptions are made for wireless sensor network localization problem.

• (1)

The coverage range of sensor node is within 300 m, which is a typical value for short range radio transmitter working in ISM (industry, science, and medicine) band, such as IEEE 802.11 b/n/g and IEEE 802.15.4a.

(2)

All sensor nodes are equipped with low cost quartz oscillators; the typical clock skew is in the range of $[\mathrm{1}{\mathrm{0}}^{-\mathrm{6}},\mathrm{1}{\mathrm{0}}^{-\mathrm{4}}]$ .

(3)

Since the anchors are connected to wired network and well powered, the power consumption of anchor is not concerned in this paper. However, the power consumption of tag is critical and there is no limitation on signal transmitting period and moving velocity of tag. Localization server where synchronization and localization algorithms are implemented is also assumed to have enough processing capability (if the number of tag exceeds a threshold, more localization servers running on desktop computers can be adopted to balance the computation load; there is no need to revise the design of anchor and tag).

3. Estimation of Relative Clock Offsets and Skews

3.1. TOA and TDOA for Anchor Blink

Suppose anchors broadcast signals by turns with a period of T with respect to their own local clocks. The broadcasting packet or pulse of the anchor is called “anchor blink.” Two consecutive blinks and TOA measurements are plotted in Figure 2.

Figure 2

Anchor blinks and local TOA measurements.

Assuming anchor q broadcasts kth blink, the TOA at anchor j is

$$\begin{array}{cc}& l_{j,q}\left(k\right)={\tilde{l}}_{j,q}\left(k\right)+n_j\left(k\right),\end{array}$$

(2)

where ${\tilde{l}}_{j,q}(k)$ is the true value of TOA and $n_j(k)$ is measurement noise that can be modeled as additive white Gaussian noise (AWGN) with mean of zero and variance of ${\sigma }_n^{\mathrm{2}}$ [6, 17]. In this paper, we focus on the line-of-sight (LOS) propagation between sensor nodes. Define ${\phi }_{ij}(k)$ and $y_{ij}(k)$ to be the true and measured relative clock offset between anchor i and anchor j when anchor q transmits blink k. Consider

$$\begin{array}{cc}& {\phi }_{ij}\left(k\right)={\tilde{l}}_{i,q}\left(k\right)-{\tilde{l}}_{j,q}\left(k\right)-\left(d_{iq}-d_{jq}\right),\end{array}$$

(3)

$$\begin{array}{cc}& y_{ij}\left(k\right)=l_{i,q}\left(k\right)-l_{j,q}\left(k\right)-\left(d_{iq}-d_{jq}\right),\end{array}$$

(4)

where $d_{jp}=(\mathrm{1}/c)‖{\mathbf{x}}_j-{\mathbf{x}}_p‖$ is propagation delay between anchor j and anchor p, and c is the velocity of light. Substituting (2) and (3) into (4), we have

$$\begin{array}{cc}& y_{ij}\left(k\right)={\phi }_{ij}\left(k\right)+n_i\left(k\right)-n_j\left(k\right).\end{array}$$

(5)

The true value $\tilde{T}(k+\mathrm{1},k)$ of the interval of blink k can be obtained from the TOAs of two consecutive blinks at each anchor

$$\begin{array}{l}\tilde{T}\left(k+\mathrm{1},k\right)=\frac{{\tilde{l}}_{i,p}\left(k+\mathrm{1}\right)-d_{ip}-\left[{\tilde{l}}_{i,q}\left(k\right)-d_{iq}\right]}{\mathrm{1}+{\epsilon }_i}\\ =\frac{{\tilde{l}}_{j,p}\left(k+\mathrm{1}\right)-d_{jp}-\left[{\tilde{l}}_{j,q}\left(k\right)-d_{jq}\right]}{\mathrm{1}+{\epsilon }_j}.\end{array}$$

(6)

Since the clock skew of a typical low cost quartz oscillator is in the range of $[\mathrm{1}{\mathrm{0}}^{-\mathrm{6}},\mathrm{1}{\mathrm{0}}^{-\mathrm{4}}]$ , the first-order approximation of (6) is

$$\begin{array}{cc}& \frac{{\tilde{l}}_{i,p}\left(k+\mathrm{1}\right)-{\tilde{l}}_{i,q}\left(k\right)-\left[d_{ip}-d_{iq}\right]}{{\tilde{l}}_{j,p}\left(k+\mathrm{1}\right)-{\tilde{l}}_{j,q}\left(k\right)-\left[d_{jp}-d_{jq}\right]}\approx \mathrm{1}+{\epsilon }_i-{\epsilon }_j.\end{array}$$

(7)

Define $e_{ij}={\epsilon }_i-{\epsilon }_j$ to be the relative clock skew between anchor i and anchor j. Substituting (3) into (7), we have

$$\begin{array}{cc}& {\phi }_{ij}\left(k+\mathrm{1}\right)={\phi }_{ij}\left(k\right)+e_{ij}\left[{\tilde{l}}_{j,p}\left(k+\mathrm{1}\right)-{\tilde{l}}_{j,q}\left(k\right)-\left(d_{ip}-d_{iq}\right)\right].\end{array}$$

(8)

Notice that the coverage range of sensor node is within 300 m, $(d_{ip}-d_{iq})$ is smaller than 1 μs, and $(d_{ip}-d_{iq})e_{ij}$ is much smaller than $\mathrm{1}{\mathrm{0}}^{-\mathrm{10}}$ and may be ignored. If the noise of TOA measurement is also much smaller than 1 μs (TOA measurement standard deviation is inversely proportional to the effective bandwidth of signal [18]; for ultrawideband system (UWB) complying with IEEE 802.15.4a standards, the measurement noise may be only several nanoseconds [12]), we have

$$\begin{array}{cc}& {\phi }_{ij}\left(k+\mathrm{1}\right)\approx {\phi }_{ij}\left(k\right)+e_{ij}\left[l_{j,p}\left(k+\mathrm{1}\right)-l_{j,q}\left(k\right)\right].\end{array}$$

(9)

Comparing (8) and (9), it can be seen that the true value of blink interval can be replaced by local TOA measurement. Define $T_j(k+\mathrm{1},k)=l_{j,p}(k+\mathrm{1})-l_{j,q}(k)$ to be the time interval of two consecutive TOAs at anchor j with respect to local clock. The relative clock offset becomes

$$\begin{array}{cc}& {\phi }_{ij}\left(k+\mathrm{1}\right)={\phi }_{ij}\left(k\right)+e_{ij}{T}_j\left(k+\mathrm{1},k\right).\end{array}$$

(10)

3.2. Relative Clock Offset and Skew Estimation

Suppose anchor i transmits blink $k+i-\mathrm{1}$ , and we can obtain the $N\mathrm{\ast }(N-\mathrm{1})$ noisy measurement of relative clock offset from N blinks:

$$\begin{array}{l}{y}_{i\mathrm{1}}\left(k+m\right)={\phi }_{i\mathrm{1}}\left(k+m\right)+n_i\left(k+m\right)-n_{\mathrm{1}}\left(k+m\right)\\ i=\mathrm{2},\dots ,N,\hspace{1em}m=\mathrm{0},\dots ,N-\mathrm{1}.\end{array}$$

(11)

From (10), we have

$$\begin{array}{cc}& {\phi }_{i\mathrm{1}}\left(k+m\right)={\phi }_{i\mathrm{1}}\left(k+N-\mathrm{1}\right)-e_{i\mathrm{1}}{T}_{\mathrm{1}}\left(k+N-\mathrm{1},k+m\right).\end{array}$$

(12)

Since clock frequency changes slightly over times, the relative clock skews of anchors are assumed to be constant during N blinks if value of blink interval is appropriately selected. Substituting (12) into (11), we can get $N-\mathrm{1}$ measurement equations for $i=\mathrm{2},\dots ,N$

$$\begin{array}{l}{y}_{i\mathrm{1}}\left(k+m\right)={\phi }_{i\mathrm{1}}\left(k+N-\mathrm{1}\right)-e_{i\mathrm{1}}{T}_{\mathrm{1}}\left(k+N-\mathrm{1},k+m\right)\\ +n_i\left(k+m\right)-n_{\mathrm{1}}\left(k+m\right).\end{array}$$

(13)

Stacking $N-1$ measurement equation into matrix form,

$$\begin{array}{cc}& \mathbf{y}\left(k+m\right)=\mathbf{h}\left(k+m\right)\mathbf{s}\left(k\right)+{\mathbf{C}}_{\mathrm{0}}\mathbf{n}\left(k+m\right),\end{array}$$

(14)

where

$$\begin{array}{cc}& \mathbf{y}\left(k+m\right)={\left[y_{\mathrm{21}}\left(k+m\right),\dots ,y_{N\mathrm{1}}\left(k+m\right)\right]}^T,\\ \mathbf{s}\left(k\right)={\left[\begin{array}{cccc}{\mathbf{s}}_{\mathrm{21}}^T\left(k\right)& {\mathbf{s}}_{\mathrm{31}}^T\left(k\right)& \cdots & {\mathbf{s}}_{N\mathrm{1}}^T\left(k\right)\end{array}\right]}^T,\\ {\mathbf{s}}_{i\mathrm{1}}\left(k\right)={\left[{\phi }_{i\mathrm{1}}\left(k+N-\mathrm{1}\right),e_{i\mathrm{1}}\right]}^T,\\ \mathbf{h}\left(k+m\right)=I_{N-\mathrm{1}}\otimes \left[\mathrm{1},-T\left(k+N-\mathrm{1},k+m\right)\right],\\ {\mathbf{C}}_{\mathrm{0}}={\left[\begin{array}{ccccc}-\mathrm{1}& \mathrm{1}& & & \\ -\mathrm{1}& & \mathrm{1}& & \\ ⋮& & & \ddots & \\ -\mathrm{1}& & & & \mathrm{1}\end{array}\right]}_{\left(N-\mathrm{1}\right)\times N},\\ \mathbf{n}\left(k+m\right)={\left[\begin{array}{cccc}{n}_{\mathrm{1}}\left(k+m\right)& n_{\mathrm{2}}\left(k+m\right)& \cdots & n_N\left(k+m\right)\end{array}\right]}^T,\\ n_i\left(k+m\right)~N\left(\mathrm{0},{\sigma }_n^{\mathrm{2}}\right).\end{array}$$

(15)

It should be addressed that ${\mathbf{I}}_N$ represents $N\times N$ matrix of identity and the operator “⊗” denotes Kronecker product, respectively.

After each anchor blinks once, we can get $N\mathrm{\ast }(N-\mathrm{1})$ measurement equations for consecutive N blinks; combining all these equations into matrix form, we have

$$\begin{array}{cc}& \mathbf{Y}\left(k\right)=\mathbf{H}\left(k\right)\mathbf{s}\left(k\right)+\mathbf{C}·\mathbf{u}\left(k\right),\end{array}$$

(16)

where

$$\begin{array}{cc}& \mathbf{Y}\left(k\right)={\left[{\mathbf{y}}^T\left(k\right),{\mathbf{y}}^T\left(k+\mathrm{1}\right),\dots ,{\mathbf{y}}^T\left(k+N-\mathrm{1}\right)\right]}^T,\\ \mathbf{H}\left(k\right)={\left[\begin{array}{cccc}{\mathbf{h}}^T\left(k\right)& {\mathbf{h}}^T\left(k+\mathrm{1}\right)& \cdots & {\mathbf{h}}^T\left(k+N-\mathrm{1}\right)\end{array}\right]}^T,\\ \mathbf{u}\left(k\right)={\left[\begin{array}{cccc}{\mathbf{n}}_{}^T\left(k\right)& {\mathbf{n}}_{}^T\left(k+\mathrm{1}\right)& \cdots & {\mathbf{n}}_{}^T\left(k+N-\mathrm{1}\right)\end{array}\right]}^T,\\ \mathbf{C}={\mathbf{I}}_N\otimes {\mathbf{C}}_{\mathrm{0}}.\end{array}$$

(17)

Since the noise of TOA measurement is independent,

$$\begin{array}{cc}& {\mathbf{Q}}_u={\sigma }_n^{\mathrm{2}}·{\mathbf{I}}_{N^{\mathrm{2}}}.\end{array}$$

(18)

The LS estimation of $\widehat{\mathbf{s}}(k)$ is [19]

$$\begin{array}{cc}& \widehat{\mathbf{s}}\left(k\right)={\left[{\mathbf{H}}^T\left(k\right){\mathbf{Q}}_c^{-\mathrm{1}}\mathbf{H}\left(k\right)\right]}^{-\mathrm{1}}{\mathbf{H}}^T\left(k\right){\mathbf{Q}}_c^{-\mathrm{1}}\mathbf{Y}\left(k\right),\end{array}$$

(19)

where ${\mathbf{Q}}_c^{}={\mathbf{C}}^T{\mathbf{Q}}_u\mathbf{C}$ . The covariance of $\widehat{\mathbf{s}}(k)$ is

$$\begin{array}{cc}& {\mathbf{Q}}_{s\left(k\right)}={\left[{\mathbf{H}}^T\left(k\right){\mathbf{Q}}_c^{-\mathrm{1}}\mathbf{H}\left(k\right)\right]}^{-\mathrm{1}}.\end{array}$$

(20)

For the Gaussian linear model, the CRLB of element n of $\widehat{\mathbf{s}}(k)$ is diagonal element n of ${\mathbf{Q}}_{s(k)}$ [19].

The procedure of the relative clock offset and skew estimation is described as follows.

• (1)

The localization server collects the TOAs of N continuous anchor blinks (each anchor blinks once).

(2)

Evaluate the value of blink interval using $T_j(k+\mathrm{1},k)=l_{j,p}(k+\mathrm{1})-l_{j,q}(k)$ , and form measurement matrix $\mathbf{H}(k)$ .

(3)

Calculate measurement values of anchor relative clock offset using (4), and combine these values into measurement vector $\mathbf{Y}(k)$ .

(4)

Calculate covariance matrix of measurement noise of relative clock offsets using ${\mathbf{Q}}_c^{}={\mathbf{C}}^T{\mathbf{Q}}_u\mathbf{C}$ .

(5)

Get the LS estimation and CRLB of relative clock offset and skew based on (19) and (20), respectively.

4. Localization of Tag

4.1. Calibration of TDOA Measurement

Since the tag transmits signal based on its own clock, some kind of media access method may be adopted to avoid collision of node transmissions. The TOAs of tag signals at anchors are depicted in Figure 3.

Figure 3

TOAs of tag signals at anchors.

Define ${\psi }_{ij}(kT+{\tau }_{\mathrm{0}})$ as the original TDOA of the tag signal between anchor i and anchor j

$$\begin{array}{l}{\psi }_{ij}\left(kT+{\tau }_{\mathrm{0}}\right)=l_{i,\mathrm{0}}\left(kT+{\tau }_{\mathrm{0}}\right)-l_{j,\mathrm{0}}\left(kT+{\tau }_{\mathrm{0}}\right),\end{array}$$

(21)

$$\begin{array}{l}{\psi }_{ij}\left(kT+{\tau }_{\mathrm{0}}\right)={\phi }_{ij}\left(kT+{\tau }_{\mathrm{0}}\right)+d_{i\mathrm{0}}-d_{j\mathrm{0}}\\ \hspace{1em}+n_{i,\mathrm{0}}\left(kT+{\tau }_{\mathrm{0}}\right)-n_{j,\mathrm{0}}\left(kT+{\tau }_{\mathrm{0}}\right),\end{array}$$

(22)

where ${\tau }_{\mathrm{0}}$ is the elapsed time from the last anchor blink. From (22), we can see that the original TDOA is corrupted by relative clock offset which should be removed to get an accurate TDOA. The estimated relative clock offset at time $kT+{\tau }_{\mathrm{0}}$ is

$$\begin{array}{cc}& {\widehat{\phi }}_{i\mathrm{1}}\left(kT+{\tau }_{\mathrm{0}}\right)={\widehat{\phi }}_{i\mathrm{1}}\left(k\right)+{\tau }_{\mathrm{0}}{\widehat{e}}_{i\mathrm{1}}\left(k\right).\end{array}$$

(23)

Comparing (23) and (9), the elapsed time from the last anchor blink can be replaced by local TOA measurements of anchor #1; we have

$$\begin{array}{cc}& {\tau }_{\mathrm{0}}\approx l_{\mathrm{1,0}}\left(kT+{\tau }_{\mathrm{0}}\right)-l_{\mathrm{1},i}\left(k\right).\end{array}$$

(24)

The calibrated TDOA $z_{i\mathrm{1}}(kT+{\tau }_{\mathrm{0}})$ is obtained by subtracting the estimated relative clock offset ${\widehat{\phi }}_{i\mathrm{1}}(kT+{\tau }_{\mathrm{0}})$ from ${\psi }_{i\mathrm{1}}(kT+{\tau }_{\mathrm{0}})$

$$\begin{array}{cc}& z_{i\mathrm{1}}\left(kT+{\tau }_{\mathrm{0}}\right)={\psi }_{i\mathrm{1}}\left(kT+{\tau }_{\mathrm{0}}\right)-\left[{\widehat{\phi }}_{i\mathrm{1}}\left(k\right)+{\tau }_{\mathrm{0}}{\widehat{e}}_{i\mathrm{1}}\left(k\right)\right].\end{array}$$

(25)

Combine the $N-1$ calibrated TDOA into matrix form

$$\begin{array}{cc}& \mathbf{z}=f\left({\mathbf{x}}_{\mathrm{0}}\right)+\mathbf{w}\mathbf{,}\end{array}$$

(26)

where $\mathbf{z}={\left[z_{\mathrm{21}}\left(kT+{\tau }_{\mathrm{0}}\right),\dots ,z_{N\mathrm{1}}(kT+{\tau }_{\mathrm{0}})\right]}^T$ , and

$$\begin{array}{cc}& f\left({\mathbf{x}}_{\mathrm{0}}\right)={\left[d_{\mathrm{2,0}}-d_{\mathrm{1,0}},d_{\mathrm{3,0}}-d_{\mathrm{1,0}},\dots ,d_{N,\mathrm{0}}-d_{\mathrm{1,0}}\right]}^T,\end{array}$$

(27)

$$\begin{array}{cc}& \mathbf{w}=\mathbf{D}{\mathbf{s}}_e\left(k\right)+{\mathbf{C}}_{\mathrm{0}}{\mathbf{u}}_{\mathrm{0}}\left(kT+{\tau }_{\mathrm{0}}\right),\end{array}$$

(28)

where

$$\begin{array}{l}{\mathbf{s}}_e\left(k\right)=\mathbf{s}\left(k\right)-\widehat{\mathbf{s}}\left(k\right),\\ {\mathbf{u}}_{\mathrm{0}}\left(kT+{\tau }_{\mathrm{0}}\right)\\ \hspace{1em}={\left[n_{\mathrm{1,0}}\left(kT+{\tau }_{\mathrm{0}}\right),n_{\mathrm{2,0}}\left(kT+{\tau }_{\mathrm{0}}\right),\dots ,n_{N,\mathrm{0}}\left(kT+{\tau }_{\mathrm{0}}\right)\right]}^T,\\ {\mathbf{u}}_{\mathrm{0}}\left(kT+{\tau }_{\mathrm{0}}\right)~N\left(\mathrm{0},{\sigma }_n^{\mathrm{2}}{\mathbf{I}}_N\right),\\ \mathbf{D}={\left[\begin{array}{llllllll}\mathrm{1}& {\tau }_{\mathrm{0}}& \mathrm{0}& & \cdots & \cdots & & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& {\tau }_{\mathrm{0}}& \mathrm{0}& \cdots & & ⋮\\ \cdots & & & \cdots & \cdots & & & \mathrm{0}\\ \mathrm{0}& & & \cdots & & \mathrm{0}& \mathrm{1}& {\tau }_{\mathrm{0}}\end{array}\right]}_{\left(N-\mathrm{1}\right)\times \left(\mathrm{2}N-\mathrm{2}\right)}.\end{array}$$

(29)

Since ${\mathbf{s}}_e(k)$ and ${\mathbf{u}}_{\mathrm{0}}(kT+{\tau }_{\mathrm{0}})$ are independent, the covariance matrix of w is

$$\begin{array}{cc}& {\mathbf{Q}}_w=E\left(\mathbf{w}{\mathbf{w}}^T\right)=\mathbf{C}{\mathbf{Q}}_u{\mathbf{C}}^T+\mathbf{D}{\mathbf{Q}}_{s\left(k\right)}{\mathbf{D}}^T.\end{array}$$

(30)

The likelihood function is

$$\begin{array}{l}p\left(\mathbf{z},{\mathbf{x}}_{\mathrm{0}}\right)=\frac{\mathrm{1}}{{\left(\mathrm{2}\pi \right)}^{\left(N-\mathrm{1}\right)/\mathrm{2}}\det {\left({\mathbf{Q}}_w\right)}^{\mathrm{1}/\mathrm{2}}}\\ ·\exp \left[-\frac{\mathrm{1}}{\mathrm{2}}{\left(\mathbf{z}-f\left({\mathbf{x}}_{\mathrm{0}}\right)\right)}^T{\mathbf{Q}}_w^{-\mathrm{1}}\left(\mathbf{z}-f\left({\mathbf{x}}_{\mathrm{0}}\right)\right)\right].\end{array}$$

(31)

Notice that the clock skews of anchors are assumed to be constant in (19) and (25). Considering random clock drift due to frequency fluctuation, blink interval T must be selected carefully to guarantee the accuracy of localization. The choice of blink interval will be further discussed in Section 5.

4.2. Maximum Likelihood Estimation

The MLE of ${\mathbf{x}}_{\mathrm{0}}$ is [19]

$$\begin{array}{cc}& {\widehat{\mathbf{x}}}_{\mathrm{0}}=\arg \hspace{0.17em}\min \left[{\left(\mathbf{z}-f\left({\mathbf{x}}_{\mathrm{0}}\right)\right)}^T{\mathbf{Q}}_w^{-\mathrm{1}}\left(\mathbf{z}-f\left({\mathbf{x}}_{\mathrm{0}}\right)\right)\right].\end{array}$$

(32)

Under the Gaussian noise assumption, the MLE has a least-squares interpretation. But a closed-form solution of (32) does not exist in general due to the nonlinear function $f({\mathbf{x}}_{\mathrm{0}})$ . Numerical minimization is thus needed; a successive linearization procedure [6, 17, 19, 20] is summarized as follows.

(1) Let the estimation at the mth iteration be ${\widehat{\mathbf{x}}}_{\mathrm{0}}(m)$ and ${\mathbf{x}}_{\mathrm{0}}={\widehat{\mathbf{x}}}_{\mathrm{0}}(m)+\Delta (m)$ . Linearizing $f({\mathbf{x}}_{\mathrm{0}})$ around ${\widehat{\mathbf{x}}}_{\mathrm{0}}(m)$ yields

$$\begin{array}{cc}& f\left({\mathbf{x}}_{\mathrm{0}}\right)\approx f\left({\widehat{\mathbf{x}}}_{\mathrm{0}}\left(m\right)\right)+G\left({\widehat{\mathbf{x}}}_{\mathrm{0}}\left(m\right)\right)\Delta \left(m\right),\end{array}$$

(33)

where $G({\widehat{\mathbf{x}}}_{\mathrm{0}}(m))$ is the Jacobian matrix

$$\begin{array}{cc}& \mathbf{G}\left({\mathbf{x}}_{\mathrm{0}}\right)=\frac{\partial f\left({\mathbf{x}}_{\mathrm{0}}\right)}{\partial {\mathbf{x}}_{\mathrm{0}}}.\end{array}$$

(34)

Substituting (33) into (32) and solving the linearized minimization problem for $\Delta (m)$ yield

$$\begin{array}{l}\widehat{\Delta }\left(m\right)={\left[{\mathbf{G}}^T\left({\widehat{\mathbf{x}}}_{\mathrm{0}}\left(m\right)\right){\mathbf{Q}}_w^{-\mathrm{1}}\mathbf{G}\left({\widehat{\mathbf{x}}}_{\mathrm{0}}\left(m\right)\right)\right]}^{-\mathrm{1}}\\ \hspace{1em}\hspace{1em}·{\mathbf{G}}^T\left({\widehat{\mathbf{x}}}_{\mathrm{0}}\left(m\right)\right){\mathbf{Q}}_w^{-\mathrm{1}}\left[z-f\left({\widehat{\mathbf{x}}}_{\mathrm{0}}\left(m\right)\right)\right],\end{array}$$

(35)

where

$$\begin{array}{cc}& \mathbf{G}\left({\mathbf{x}}_{\mathrm{0}}\right)=\frac{\mathrm{1}}{c}{\left[{\mathbf{r}}_{\mathrm{2}}^T\left({\mathbf{x}}_{\mathrm{0}}\right)-{\mathbf{r}}_{\mathrm{1}}^T\left({\mathbf{x}}_{\mathrm{0}}\right),\dots ,{\mathbf{r}}_N^T\left({\mathbf{x}}_{\mathrm{0}}\right)-{\mathbf{r}}_{\mathrm{1}}^T\left({\mathbf{x}}_{\mathrm{0}}\right)\right]}^T.\end{array}$$

(36)

${\mathbf{r}}_i({\mathbf{x}}_{\mathrm{0}})$ defines a unit-norm direction vectors:

$$\begin{array}{cc}& {\mathbf{r}}_i\left({\mathbf{x}}_{\mathrm{0}}\right)=\frac{{\mathbf{x}}_{\mathrm{0}}-{\mathbf{x}}_i}{‖{\mathbf{x}}_{\mathrm{0}}-{\mathbf{x}}_i‖},\hspace{1em}i=\mathrm{1},\dots ,N.\end{array}$$

(37)

(2) The estimation at the $(m+\mathrm{1})$ th iteration is

$$\begin{array}{cc}& {\widehat{\mathbf{x}}}_{\mathrm{0}}\left(m+\mathrm{1}\right)={\widehat{\mathbf{x}}}_{\mathrm{0}}\left(m\right)+\widehat{\Delta }\left(m\right).\end{array}$$

(38)

The iteration starts with an initial guess ${\widehat{\mathbf{x}}}_{\mathrm{0}}(\mathrm{0})$ and terminates at convergence. When a tag transmits signals periodically, previous estimate can serve as good initial guess of the current value of ${\widehat{\mathbf{x}}}_{\mathrm{0}}$ . However, the iteration may stop at a local minimum and may not converge when $[{\mathbf{G}}^T({\widehat{\mathbf{x}}}_{\mathrm{0}}(m)){\mathbf{Q}}_w^{-\mathrm{1}}\mathbf{G}({\widehat{\mathbf{x}}}_{\mathrm{0}}(m)){]}^{-\mathrm{1}}$ is large. A two-step process, starting with a coarse grid search and continuing with an iterative procedure, can be adopted to search for the global minimum.

The procedure of the tag localization algorithm is described as follows.

• (1)

The localization server collects the tag signal's TOAs at anchors and calculates original $N-1$ TDOA value using (21).

(2)

Calibrate TDOAs by eliminating relative clock offset errors using (22).

(3)

Calculate covariance matrix ${\mathbf{Q}}_w$ , set $m=1$ , and give an initial value ${\widehat{\mathbf{x}}}_{\mathrm{0}}(\mathrm{1})$ of tag position randomly.

(4)

Calculate $f({\widehat{\mathbf{x}}}_{\mathrm{0}}(m))$ and step vector $\widehat{\Delta }(m)$ using (27) and (35), respectively.

(5)

Update tag position with (38); if position change is larger than 0.1 mm, set $m=m+1$ and go to step (4) to continue iteration. Otherwise, the iteration stops and estimate of tag position is achieved.

4.3. CRLB for Tag Localization Algorithm

The localization performance of the proposed algorithm is analyzed using the CRLB. The CRLB is the inverse of the Fisher information matrix. For Gaussian noise, the Fisher information matrix has the form [6, 19]

$$\begin{array}{cc}& J\left({\mathbf{x}}_{\mathrm{0}}\right)={\left(\frac{\partial f\left({\mathbf{x}}_{\mathrm{0}}\right)}{\partial {\mathbf{x}}_{\mathrm{0}}}\right)}^T{\mathbf{Q}}_w^{-\mathrm{1}}\left(\frac{\partial f\left({\mathbf{x}}_{\mathrm{0}}\right)}{\partial {\mathbf{x}}_{\mathrm{0}}}\right)=G^T\left({\mathbf{x}}_{\mathrm{0}}\right){\mathbf{Q}}_w^{-\mathrm{1}}G\left({\mathbf{x}}_{\mathrm{0}}\right).\end{array}$$

(39)

Assuming that $J^{-\mathrm{1}}({\mathbf{x}}_{\mathrm{0}})$ exists and the nth diagonal element of $J^{-\mathrm{1}}({\mathbf{x}}_{\mathrm{0}})$ is denoted by $[J^{-\mathrm{1}}({\mathbf{x}}_{\mathrm{0}}){]}_{nn}$ for $n=\mathrm{1,2},\mathrm{3}$ , the variance of any element $x_{\mathrm{0},n}$ of ${\mathbf{x}}_{\mathrm{0}}$ is bounded below by $[J^{-\mathrm{1}}({\mathbf{x}}_{\mathrm{0}}){]}_{nn}$ . The CRLB of the position estimation is

$$\begin{array}{cc}& {\stackrel{-}{\sigma }}_x^{\mathrm{2}}={{\displaystyle \sum }}_{n=\mathrm{1}}^{\mathrm{3}}{\left[J^{-\mathrm{1}}\left({\mathbf{x}}_{\mathrm{0}}\right)\right]}_{nn}.\end{array}$$

(40)

Another performance measure commonly used in source localization systems is the geometric dilution of precision (GDOP) [21]. The GDOP is the magnification in localization error due to the geometric relationship between the anchors and tags. Let $\text{Cov}({\widehat{\mathbf{x}}}_{\mathrm{0}})$ be the covariance matrix of an unbiased position estimate of ${\widehat{\mathbf{x}}}_{\mathrm{0}}$ . The GDOP is defined and related to the CRLB as

$$\begin{array}{cc}& \text{GDOP}=\frac{\sqrt{\mathrm{tr}\left(\text{Cov}\left({\widehat{\mathbf{x}}}_{\mathrm{0}}\right)\right)}}{c{\sigma }_n}=\frac{{\stackrel{-}{\sigma }}_x}{c{\sigma }_n},\end{array}$$

(41)

where $\mathrm{tr}(·)$ denotes matrix trace.

5. Numerical Results

In this section, numerical results for both estimations of anchors' clock parameters and tag's localization are presented. A two-dimensional localization scene is considered in this paper; four anchors are placed evenly at coordinates of (0, 0), (100 m, 0), (0, 100 m), and (100 m, 100 m), respectively. The anchors stamp TOAs of signal emitted by other anchors or tags and deliver the TOA information to localization server where estimation of clock parameters and tag location are implemented.

Firstly, Monte Carlo simulations were carried out to evaluate the statistical performance of the proposed algorithm without taking the clock drift into consideration. The simulation result unfolds some theoretical features of the proposed algorithm. Secondly, simulations in the presence of clock drift are carried out and the appropriate value of blink interval is analyzed. At last, the GDOP performances for different layouts of anchors are evaluated.

5.1. Performance without Clock Drift

5.1.1. Results of Relative Clock Offset and Skew Estimation

Without considering clock drift, the performance of synchronization was evaluated as a function of TOA measurement. The actual root mean square error (RMSE) and root CRLB (RCRLB) of relative clock offset and skew estimates with different values of blink interval T are plotted in Figures 4 and 5. The RMSE results were averaged over 50000 independent noise measurements.

Figure 4

RCRLB and RMSE for relative clock offset estimation versus TOA measurement standard deviation (clock drift not concerned).

Figure 5

RCRLB and RMSE for clock skews estimation versus TOA measurement standard deviation (clock drift not concerned).

From Figure 4, it can be seen that the value of blink interval has no effect on the performance of relative clock offset estimation. We can see that the RMSE and CRLB of relative clock skew estimation are inversely proportional to the value of blink interval in Figure 5. It is inferred in (20) that estimate of relative clock skew is a weighted average of TOAs over time span $(N-\mathrm{1})T$ , so larger value of blink interval means less error of relative clock skew estimate. Although we can get more accurate estimate with larger value of blink interval theoretically, it is not practical because it is implied in (20) that clock skews of anchors are assumed to be constant and clock drifts are not taken into account.

The performance of tag localization is given in Figure 6. It seems that blink interval has no effect on the performance of tag localization. It has been concluded that the RMSE and RCRLB of relative clock skew estimation are inversely proportional to the value of blink interval; it can be seen in (23) that ${\tau }_{\mathrm{0}}{\widehat{e}}_{i\mathrm{1}}(k)$ is proportional to the value of elapsed time from the last blink. Since the average value of ${\tau }_{\mathrm{0}}$ is $T/\mathrm{2}$ , the standard deviation of residual clock offset in calibrated TDOA $z_{i\mathrm{1}}(kT+{\tau }_{\mathrm{0}})$ remains constant whatever the value of blink interval is. Notice that the performance curves in Figure 5 are achieved under the assumptions that clock skews of anchors are constant and clock drifts are not taken into account. It is unreasonable to calibrate the original TDOA with the timing estimates of an hour ago. The clock drift must be taken into consideration for practical implementation.

Figure 6

RCRLB and RMSE for tag position estimation versus TOA measurement for different blink intervals (clock drift not concerned).

5.2. Performance with Clock Drift

5.2.1. Results of Clock Difference Estimation

Clock drift is inevitable for all kinds of clock source; it must be taken into account for real implementation. In this paper, clock frequency fluctuation of quartz oscillator is modeled as AWGN [22]

$$\begin{array}{cc}& {\dot{l}}_i\left(t\right)=f_i\left(t\right)+n_f,\end{array}$$

(42)

where $n_f$ is AWGN with variance of ${\sigma }_f^{\mathrm{2}}$ , which is unknown and may depend on the characteristics of the quartz oscillator, the ambient temperature, the working voltage, and so on.

The performance of synchronization in the presence of clock drift was also evaluated as a function of TOA measurement. Since blink interval T greatly influences the accuracy of localization, the actual root mean square error (RMSE) and root CRLB (RCRLB) of relative clock offset and skew estimates with different values of blink interval T are plotted in Figures 7 and 8. The RMSE results were averaged over 50000 independent noise measurements.

Figure 7

RCRLB and RMSE for relative clock offset estimation versus TOA measurement standard deviation, ${\sigma }_f=\mathrm{1}{\mathrm{0}}^{-\mathrm{9}}$ .

Figure 8

RCRLB and RMSE for clock skews estimation versus TOA measurement standard deviation, ${\sigma }_f=\mathrm{1}{\mathrm{0}}^{-\mathrm{9}}$ .

There is an obvious gap between the RMSE and corresponding RCRLB in Figures 7 and 8, because clock skew is treated as constant in the computation procedure of RCRLB. However, frequency fluctuations are inevitable in practice; the impact of frequency fluctuations is more significant when TOA measurement standard deviation ${\sigma }_n$ is close to ${\sigma }_f$ in (42), and ${\sigma }_f$ forms low bound of clock skew estimate. The RCRLB curves of relative clock offset in Figure 7 for different values of T converge when noise level increases, because ${\sigma }_f$ is much smaller than ${\sigma }_n$ .

5.2.2. Tag Localization Result

The RMSE and RCRLB for tag localization with different blink interval T are plotted in Figure 9. For simplicity, ${\tau }_{\mathrm{0}}$ is set to be $T/\mathrm{2}$ , which is the average elapsed time from the last anchor blink. There is also an obvious gap between RMSE and corresponding RCRLB because relative clock skew is treated as constant in (19). The RCRLBs for different values of T converge when noise level increases, because residual clock drift in calibrated TDOA is much smaller than noise level of TOA measurement.

Figure 9

RCRLB and RMSE for tag position estimation versus TOA measurement for different anchor blink intervals. The simulation results for algorithms in [4, 7] are also given: ${\sigma }_f=\mathrm{1}{\mathrm{0}}^{-\mathrm{9}}$ .

The performance of the synchronous TDOA algorithm in [7] is presented for comparison in Figure 9. The dotted red lines with markers “▿” and “+” represent RCRLB and RMSE of [7], respectively. We can conclude that the RMSE and RCRLB values of our algorithm decrease and approach those of the synchronous TDOA algorithm if a smaller value of T is selected. This is due to the fact that the impact of clock drift on calibrated TDOA becomes less significant as the time interval of anchor blinks is reduced. Although more accurate estimation of the positions of tags can be achieved with a smaller blink interval T, the number of tags to be served will decrease due to the large number of anchor blink signals. We can conclude that a value in the range $[0.1,0.4]$ is adequate for practical implementation.

Furthermore, the algorithm in this paper is compared with the algorithms in [4] where a joint synchronization and localization algorithm for asynchronous WSN was proposed. The RCRLB curve of this algorithm is plotted with blue marker “◊” in Figure 9. The joint algorithm in [4] performs better than our algorithm due to two-way message exchange scheme and the averaging of noisy TOAs of multiple signals from the same tag. We can also see that the difference between pairs of algorithms in localization precision is within 3 dB, which has little effect in practice for typical standard deviation in the range of $[\mathrm{1}{\mathrm{0}}^{-\mathrm{9}},\mathrm{1}{\mathrm{0}}^{-\mathrm{8}}]$ [12]. Since the tag performs two-way message exchange with each anchor, TOA measurement capability and complex protocol are needed for tags, and energy dissipation of tag may increase with the number of anchors. However, the algorithm in this paper is simple and effective, only periodical broadcasting of anchor is required, and the tag can be located via single transmission which means least power consumption is required for mobile tag node.

Although the algorithm proposed in this paper cannot reach the performance of algorithms in [4, 7], there are less assumptions required in this paper. Clock skews in [4, 7] are treated as constant. Furthermore, it is assumed that all anchors are synchronized and their positions are known in [7]; it is also assumed that there are $L\hspace{0.17em}\hspace{0.17em}(L\ge \mathrm{3})$ anchors with known timing and locations in [4]. The locations of anchors can be achieved easily in practice, but it is difficult and impractical to get accurate timing parameters of anchors equipped with low cost quartz oscillators because quartz clock skew is influenced by working voltage, ambient temperature, and initial working condition. On the contrary, the algorithm in this paper deals with this problem and estimates the relative clock offsets and skews of anchor pairs. There is no need for external sophisticated and expensive clock sources such as GPS receivers and atom clocks, and suitable for low cost implementation. Compared with the synchronized TDOA algorithm of [4, 7], when an appropriate blink interval is selected, performance loss is within 3 dB which is acceptable for many applications.

5.3. Performance of GDOP

The GDOP of the proposed asynchronous localization method is plotted in Figure 10 as a function of the tag position and it is compared with that of the synchronous algorithm in [7]. As expected, the asynchronous method had a worse geometric condition, but we notice that the tag can be located with GDOP performance (less than 2) close to synchronous algorithm when the tag is in the rectangle formed by four anchors.

Figure 10

Constant contour of GDOPs as a function of the tag location for both the asynchronous algorithm in this paper r (solid line) and synchronous algorithm in [7] (red dotted line).

The GDOP with different layout of anchors is plotted in Figure 11. Three anchors locate on the vertex of an equilateral triangle whose edge is 100 m; the fourth anchor locates on the center of gravity. The GDOP of the synchronous algorithm in [7] is plotted for comparison. We can also conclude that GDOP performance (less than 2) is close to synchronous algorithm when the tag is in the region formed by the outside three anchors.

Figure 11

Constant contour of GDOPs (anchors locate on (0, 0), (100, 0), (50, 86.6), and (50, 28.9)) for both the asynchronous algorithm in this paper (solid line) and synchronous algorithm in [7] (red dotted line).

The influence of the number of anchors is also studied in the simulation. For 2D localization, at least three anchors are needed for tag localization, and four anchors are needed for 3D localization. The GDOP of triangular layout of three anchors is plotted in Figure 12; three anchors locate on the vertex of an equilateral triangle whose edge is 100 m. The GDOP of four anchors (the fourth anchor locates on the center of gravity) is also plotted for comparison. It can be seen that the GDOP performance of the three anchors is worse than that of the four anchors. Statistically, more TOA measurements are helpful to reduce the locating error of tag. Furthermore, redundant TOA measurements can also be used to cope with the problem of non-LOS by selecting the TOA measurements with stronger receiving signal or shorter propagation delay.

Figure 12

Comparison of GDOP with three anchors (blue solid line; anchor locates on (0, 0), (100, 0), and (50, 86.6)) and GDOP with four anchors (red dotted line; anchor locates on (0, 0), (100, 0), (50, 86.6), and (50, 28.9)).

6. Conclusion

In this paper, a novel TDOA tag localization algorithm for WSNs is presented. To synchronize clocks of anchors, each anchor broadcasts blink signals periodically; relative clock offsets and skews of anchor pairs are estimated by the LS method using the TOAs of broadcast signals at anchors. When a tag transmits a signal, the TDOA error due to the relative clock offset of the anchor pair can be eliminated using a compensation operation. Moreover, a linearized MLE is adopted to estimate the position of the tag. Compared with previous methods, the algorithm proposed in this paper is simple, energy-efficient, and particularly suitable for low cost and fully asynchronous WSNs.