TOA-Based Source Localization: A Linearization Approach Adopting Coordinate System Translation

Abstract

This paper addresses the localization of a timing signal source based on the time of arrival (TOA) measurements that are collected from nearby sensors that are position known and synchronized to each other. Generally speaking, for such TOA-based source localization, the corresponding observation equations contain nonlinear relationship between measurements and unknown parameters, which normally results in the nonexistence of any efficient unbiased estimator that attains the Cramer-Rao lower bound (CRLB). In this paper, we devise a new approach that utilizes linearization and adopts suitable coordinate system translation to eliminate nonlinearity from the converted observation equations. The performance analysis and simulation study conducted show that our proposed algorithm can achieve the CRLB when the zero-mean Gaussian and independent measurement errors are sufficiently small.

1. Introduction

Wireless sensor networks have the potential to play a very important role in various location-aware applications [1], for example event/target monitoring [2, 3]. One of such applications is the monitoring of a timing signal source, where the task is to localize the source based on the time of arrival (TOA) measurements that are collected from nearby sensors that are position known and synchronized to each other [4, 5].

Generally speaking, for such TOA-based source localization the corresponding observation equations contain nonlinear relationship between measurements and unknown parameters (i.e., the position and transmit time of the source). And such nonlinearity normally results in the nonexistence of any efficient unbiased estimator that attains the Cramer-Rao lower bound (CRLB) [6].

In [4, 5], the authors proposed several relaxation-for-convexity-based algorithms for the above-considered TOA-based source localization. However, these algorithms are biased in general and suffer from considerable localization errors especially when the source is outside the convex hull of these position-known sensors [5].

Linearization-based algorithm is another option, similar to that adopted in [7, 8] for range-based localization. The normally adopted linearization is to square the distance that is represented by a norm of corresponding vector difference such that in the converted observation equations the relationship between measurements and unknowns (including a thus produced auxiliary variable which is the nonlinear function of the basic unknown parameters) is linear. However, such linearization only produces quasilinear observation equations for the unknowns because of the introduced auxiliary variable. And the source localization algorithms based on such linearization cannot achieve the CRLB straightforward, where multistage estimation is generally needed for localization performance improvement.

In this paper, we consider the above TOA-based source localization due to its extreme importance, and our contribution in this paper can be listed as below.

• (1)

We propose a novel linearization approach that adopts suitable coordinate system translation to eliminate nonlinearity from the converted observation equations.

(2)

The performance analysis and simulation study show that our proposed algorithm can achieve the CRLB in the presence of sufficiently small zero-mean Gaussian and independent measurement errors and is not sensitive to any large or trivial absolute coordinates of the source and sensors.

The rest of this paper is arranged as follows: In Section 2, a detailed signal processing model for the considered TOA-based source localization is given. In Section 3, the usually adopted linearization for similar nonlinear observation equations is applied and the corresponding drawbacks are highlighted. In Section 4, the proposed algorithm for the considered TOA-based source localization is presented. In Section 5, our proposed algorithm is evaluated by performance analysis and simulation study. And finally, Section 6 concludes the paper.

2. Signal Processing Model

As the extension to higher dimension cases is straightforward, in this paper we only consider 2D cases, where a location/position can be represented by a column vector composed of two coordinates.

We consider a 2D scenario that can be described as follows: a timing signal source, located at an unknown position denoted as ${\mathbf{p}}_0$ , transmits a timing signal at an unknown time denoted as $t_0$ . M nearby sensors, located at known positions denoted as ${\mathbf{p}}_m$ , $m=1,\dots ,M$ and synchronized to each other, recorded their time of arrivals at these sensors, that is, TOA measurements, denoted as $t_m$ , $m=1,\dots ,M$ . All the observation equations can be represented as follows:

$$\begin{array}{cc}& t_m=t_0+\frac{1}{c}\parallel {\mathbf{p}}_0-{\mathbf{p}}_m\parallel +n_m,\hspace{1em}m=1,\dots ,M,\end{array}$$

(1)

where c is the known propagation speed. We consider the estimation of ${\mathbf{p}}_0$ only, although $t_0$ can be estimated at the same time.

If the estimated position of the source is ${\hspace{0.17em}\hspace{0.17em}\widehat{\mathbf{p}}}_0$ , a suitable performance measure is the mean square (m.s.) localization error:

$$\begin{array}{cc}& \text{m.s.}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{locazation error}=E\left\{{\parallel {\widehat{\mathbf{p}}}_0-{\mathbf{p}}_0\parallel }^2\right\}.\end{array}$$

(2)

In the presence of sufficiently small zero-mean Gaussian and independent measurement errors, the above m.s. localization error is lower bounded by the CRLB, which can be derived based on the probability density function (pdf) of all TOA measurements given all unknowns and all sensors' positions [5]. For convenience and comparison later, we just simply list it here.

The corresponding pdf can be represented as follows:

$$\begin{array}{l}f\left(t_1,\dots ,t_m\mid t_0,{\mathbf{p}}_0,{\mathbf{p}}_1,\dots ,{\mathbf{p}}_M\right)\\ \hspace{1em}\hspace{1em}=\frac{1}{\sqrt{{\left(2\pi {\sigma }_n^2\right)}^N}}{e}^{{-\left(1/2{\sigma }_n^2\right){\sum }_{m=1}^M\left(t_0-t_m+\left(1/c\right)\parallel {\mathbf{p}}_0-{\mathbf{p}}_m\parallel \right)}_2}.\end{array}$$

(3)

The Fisher information matrix is as follows:

$$\begin{array}{cc}& I=-E\left\{\left(\begin{array}{cc}\frac{{\partial }^2}{\partial {\mathit{p}}_0{\partial }^T{\mathbf{p}}_0}\ln f& \frac{{\partial }^2}{\partial t_0{\partial }^T{\mathbf{p}}_0}\ln f\\ \frac{{\partial }^2}{\partial {\mathbf{p}}_0\partial t_0}\ln f& \frac{{\partial }^2}{{\partial }^2{t}_0}\ln f\end{array}\right)\right\},\end{array}$$

(4)

where

$$\begin{array}{cc}& -E\left\{\frac{{\partial }^2}{\partial {\mathbf{p}}_0{\partial }^T{\mathbf{p}}_0}\ln f\right\}=\frac{1}{c^2{\sigma }_n^2}\sum _{m=1}^M\frac{\left({\mathbf{p}}_0-{\mathbf{p}}_n\right){\left({\mathbf{p}}_0-{\mathbf{p}}_n\right)}^T}{\parallel {\mathbf{p}}_0-{\mathbf{p}}_n\parallel },\\ -E\left\{\frac{{\partial }^2}{{\partial }^2{t}_0}\ln f\right\}=\frac{M}{{\sigma }_n^2},\\ -E\left\{\frac{{\partial }^2}{\partial {\mathbf{p}}_0\partial t_0}\ln f\right\}=\frac{1}{c^2{\sigma }_n^2}\sum _{m=1}^M\frac{\left({\mathbf{p}}_0-{\mathbf{p}}_n\right)}{\parallel {\mathbf{p}}_0-{\mathbf{p}}_n\parallel }.\end{array}$$

(5)

The CRLB is then as follows:

$$\begin{array}{cc}& \text{CRLB}=\text{TR}\left\{{\mathbf{I}}^{-1}\right\},\end{array}$$

(6)

where $\text{TR\hspace{0.17em}}\{·\}$ is the operator to calculate the trace of a matrix.

3. Normally Adopted Linearization

Normally, the observation equations in (1) are linearized as follows. Firstly, these equations are manipulated and squared on their two sides as follows:

$$\begin{array}{cc}& {\parallel {\mathbf{p}}_0-{\mathbf{p}}_m\parallel }^2=c^2{\left(t_0-t_m+n_m\right)}^2,\hspace{1em}m=1,\dots ,M\end{array}$$

(7)

Then, they are further manipulated as follows:

$$\begin{array}{l}{\parallel {\mathbf{p}}_m\parallel }^2-c^2{t}_m^2\\ \hspace{1em}\hspace{1em}=2{\mathbf{p}}_m^T{\mathbf{p}}_0-2c^2{t}_m{t}_0-\left({\parallel {\mathbf{p}}_0\parallel }^2-c^2{t}_0^2\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}+2c^2\left(t_0-t_m\right)n_m+c^2{n}_m^2,\hspace{1em}m=1,\dots ,M,\end{array}$$

(8)

or in the following vector-matrix form:

$$\begin{array}{l}\left(\begin{array}{l}{\parallel {\mathbf{p}}_1\parallel }^2-c^2{t}_1^2\\ ⋮\\ {\parallel {\mathbf{p}}_M\parallel }^2-c^2{t}_M^2\end{array}\right)\\ \hspace{1em}\hspace{1em}=\left(\begin{array}{lll}2{\mathbf{p}}_1^T& -2c^2{t}_1& 1\\ ⋮& ⋮& ⋮\\ 2{\mathbf{p}}_M^T& -2c^2{t}_M& 1\end{array}\right)\left(\begin{array}{l}{\mathbf{p}}_0\\ t_0\\ {\parallel {\mathbf{p}}_0\parallel }^2-c^2{t}_0^2\end{array}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}+\left(\begin{array}{l}2c^2(t_0-t_1)n_1+c^2{n}_1^2\\ ⋮\\ 2c^2(t_0-t_M)n_M+c^2{n}_M^2\end{array}\right).\end{array}$$

(9)

We can see that if ${\parallel {\mathbf{p}}_0\parallel }^2-c^2{t}_0^2$ is regarded as an auxiliary unknown variable (although dependent on the basic unknowns ${\mathbf{p}}_0\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}{t}_0$ ), these equations would become quasi-linear observation equations. So a two-stage estimation would be applicable to estimate the basic unknowns. In the first stage, ${\mathbf{p}}_0$ , $t_0,$ and ${\parallel {\mathbf{p}}_0\parallel }^2-c^2{t}_0^2$ all are estimated coarsely (the estimation is coarse because the relationship between the auxiliary unknown and the basic unknowns has not been utilized); in the second stage, more accurate ${\mathbf{p}}_0$ and $t_0$ can be obtained by utilizing the relationship and corresponding estimation errors.

We also observe that these converted observation equations have errors (with different variances); that is, $2c^2(t_0-t_m)n_m+c^2{n}_m^2$ for $m=1,\dots ,M$ , which contain unknown $t_0$ . This also needs a two-stage estimation: in the first stage, an LS estimation with equal weights can be performed. In the second stage, as the unknowns have been coarsely estimated, a WLS estimation with approximately known weights can be performed.

4. Proposed Algorithm for Source Localization

We firstly present the algorithm framework and then list the steps of the algorithm at the end of this section.

Before the first stage of estimation of our proposed algorithm, a suitable coordinate system translation is performed to translate the origin of the coordinate system to the position and the received time of the sensor closest (in terms of TOA measurements) to the source. The coordinate system translation can be performed as follows if the suitable sensor is sensor k:

$$\begin{array}{cc}& {\mathbf{p}}_i^{\prime }={\mathbf{p}}_i-{\mathbf{p}}_k,\hspace{1em}i=\mathrm{0,1},\dots ,M\end{array}$$

(10)

$$\begin{array}{cc}& t_i^{\prime }=t_i-t_k,\hspace{1em}i=\mathrm{0,1},\dots ,M,\end{array}$$

(11)

where

$$\begin{array}{cc}& k=\arg \hspace{0.17em}{\min }_m{t}_m.\end{array}$$

(12)

After the coordinate system translation, the observation equations in (8) become as follows:

$$\begin{array}{l}{\parallel {\mathbf{p}}_m^{\prime }\parallel }^2-c^2{t}_m^{\prime 2}\\ \hspace{1em}\hspace{1em}=2{\mathbf{p}}_m^{\prime T}{\mathbf{p}}_0^{\prime }-2c^2{t}_m^{\prime }{t}_0^{\prime }-\left({\parallel {\mathbf{p}}_0^{\prime }\parallel }^2-c^2{t}_0^{\prime 2}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}+2c^2\left(t_0^{\prime }-t_m^{\prime }\right)n_m+c^2{n}_m^2,\hspace{1em}m=1,\dots ,M.\end{array}$$

(13)

Considering that after the above coordinate system translation, the observation equation in (1) for $m=k$ becomes as follows:

$$\begin{array}{cc}& 0=t_0^{\prime }+\frac{1}{c}\parallel {\mathbf{p}}_0^{\prime }\parallel +n_k,\end{array}$$

(14)

so we have the following:

$$\begin{array}{cc}& {\parallel {\mathbf{p}}_0^{\prime }\parallel }^2-c^2{t}_0^{\prime 2}=2c^2{t}_0^{\prime }{n}_k+c^2{n}_k^2.\end{array}$$

(15)

We can see that its right-hand side only contains the items of measurement errors. The converted observation equations can be expressed in matrix-vector form as follows:

$$\begin{array}{l}\left(\begin{array}{l}{\parallel {\mathbf{p}}_1^{\prime }\parallel }^2-c^2{t}_1^{\prime 2}\\ ⋮\\ {\parallel {\mathbf{p}}_M^{\prime }\parallel }^2-c^2{t}_M^{\prime 2}\end{array}\right)\\ \hspace{1em}\hspace{1em}=\left(\begin{array}{ll}2{\mathbf{p}}_1^{\prime T}& -2c^2{t}_1^{\prime }\\ ⋮& ⋮\\ 2{\mathbf{p}}_M^{\prime T}& -2c^2{t}_M^{\prime }\end{array}\right)\left(\begin{array}{l}{\mathbf{p}}_0^{\prime }\\ \hspace{0.17em}\\ t_0^{\prime }\end{array}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}+\left(\begin{array}{l}-2c^2{t}_0^{\prime }{n}_k-c^2{n}_k^2+2c^2\left(t_0^{\prime }-t_1^{\prime }\right)n_1+c^2{n}_1^2\\ ⋮\\ -2c^2{t}_0^{\prime }{n}_k-c^2{n}_k^2+2c^2\left(t_0^{\prime }-t_M^{\prime }\right)n_M+c^2{n}_M^2\end{array}\right),\end{array}$$

(16)

or simply as follows:

$$\begin{array}{cc}& {\mathbf{b}}^{\prime }={\mathbf{A}}^{\prime }{\mathbf{q}}^{\prime }+{\mathbf{w}}^{\prime }.\end{array}$$

(17)

We can see that the nonlinearity has been eliminated. As the variance of the errors in the above equations are not known in the first stage of the estimation, so the LS solution with equal weights is adopted here:

$$\begin{array}{cc}& {\widehat{\mathbf{q}}}^{\prime }={\left({\mathbf{A}}^{\prime T}{\mathbf{A}}^{\prime }\right)}^{-1}{\mathbf{A}}^{\prime T}{\mathbf{b}}^{\prime }=\left(\begin{array}{c}{\mathbf{p}}_0^{\prime }\\ \hspace{0.17em}\\ t_0^{\prime }\end{array}\right).\end{array}$$

(18)

We can see that when there is no measurement error, the above-estimated unknowns, that is, ${\widehat{\mathbf{p}}}_0^{\prime }$ and  ${\widehat{t}}_0^{\hspace{0.17em}\prime }$ , must be the correct ones.

Before the second stage of estimation, the suitable coordinate system translation is performed again to translate the origin of the coordinate system to the estimated position and the estimated transmitted time of the source in the first stage of estimation:

$$\begin{array}{cc}& {\mathbf{p}}_i^{\prime \prime }={\mathbf{p}}_i^{\prime }-{\widehat{\mathbf{p}}}_0^{\prime },\hspace{1em}i=\mathrm{0,1},\dots ,M,\\ t_i^{\prime \prime }=t_i^{\prime }-{\widehat{t}}_0^{\hspace{0.17em}\prime },\hspace{1em}i=\mathrm{0,1},\dots ,M.\end{array}$$

(19)

So ${\mathbf{p}}_i^{\mathrm{\prime \prime }}$ and $t_i^{\prime \prime }$ only contain the items of measurement errors.

After the coordinate system translation, the observation equations in (8) or (13) become as follows:

$$\begin{array}{l}{\parallel {\mathbf{p}}_m^{\prime \prime }\parallel }^2-c^2{t}_m^{\prime \prime 2}\\ \hspace{1em}=2{\mathbf{p}}_m^{\prime \prime T}{\mathbf{p}}_0^{\prime \prime }-2c^2{t}_m^{\prime \prime }{t}_0^{\prime \prime }-\left({\parallel {\mathbf{p}}_0^{\prime \prime }\parallel }^2-c^2{t}_0^{\prime \prime 2}\right)\\ \hspace{1em}\hspace{1em}+2c^2\left(t_0^{\prime \prime }-t_m^{\prime \prime }\right)n_m+c^2{n}_m^2,\hspace{1em}m=1,\dots ,M.\end{array}$$

(20)

Considering that after the above coordinate system translation. ${\mathbf{p}}_0^{\mathrm{\prime \prime }}$ and $t_0^{\prime \prime }$ only contain the items of measurement errors, So, when the measurement errors are sufficiently small, item ${\parallel {\mathbf{p}}_0^{\mathrm{\prime \prime }}\parallel }^2$ only contains the items of the second-order measurement errors and can be ignored as compared to item $2{\mathbf{p}}_m^{\mathrm{\prime \prime }T}{\mathbf{p}}_0^{\mathrm{\prime \prime }}$ that coexists in (20). Similarly, item $c^2{t}_0^{\prime \prime 2}$ only contains the items of the second-order measurement errors and can be ignored as compared to item $2c^2{t}_m^{\prime \prime }{t}_0^{\prime \prime }$ that coexists in (20). Thus, ${\parallel {\mathbf{p}}_0^{\mathrm{\prime \prime }}\parallel }^2-c^2{t}_0^{\prime \prime 2}$ in (20) is an item that can be ignored when the measurement errors are sufficiently small. The approximation to remove ${\parallel {\mathbf{P}}_0\parallel }^2$ and $c^2{t}_0^{\prime \prime 2}$ can be justified as follows. Take $c^2{t}_0^{\prime \prime 2}$ as an example. Such a term is insignificant compared to $2c^2{t}_m^{\prime \prime }{t}_0^{\prime \prime }$ , which can be demonstrated as follows. Consider a scenario with a fixed c and $t_m^{\prime \prime }$ , and the estimation error $t_0^{\prime \prime }$ is complex Gaussian distributed with zero mean and variance ${\sigma }^2$ , that is, $\text{CN}(0,{\sigma }^2)$ . Consider the following probability:

$$\begin{array}{cc}& P\left(\frac{\left|c^2{t}_0^{\prime \prime 2}\right|}{2\left|c^2{t}_m^{\prime \prime }{t}_0^{\prime \prime }\right|}<\epsilon \right),\end{array}$$

(21)

where ε is a fixe positive number, that is, a threshold to ignore the denominator. For example, a threshold of $\epsilon =0.01$ will be a sufficient condition to ignore $c^2{t}_0^{\prime \prime 2}$ for many practical scenarios. And the statement is equivalent to the following equation:

$$\begin{array}{cc}& P\left(\frac{\left|c^2{t}_0^{\prime \prime 2}\right|}{2\left|c^2{t}_m^{\prime \prime }{t}_0^{\prime \prime }\right|}<\epsilon \right)⟶1,\end{array}$$

(22)

when the variance of the estimation error goes to zero ${\sigma }^2\to 0$ . This probability can be rewritten as follows:

$$\begin{array}{cc}& P\left(\left|t_0^{\prime \prime }\right|<2\left|t_m^{\prime \prime }\right|\epsilon \right)=1-e^{-\left(2{\left|t_m^{\prime \prime }\right|}^2{\epsilon }^2/{\sigma }^2\right)},\end{array}$$

(23)

where we have used the cumulative distribution function of Rayleigh distribution. As can be observed from (23), for any fixed ε, by decreasing ${\sigma }^2$ to zero, we can have the following:

$$\begin{array}{cc}& P\left(\frac{\left|c^2{t}_0^{\prime \prime 2}\right|}{2\left|c^2{t}_m^{\prime \prime }{t}_0^{\prime \prime }\right|}\right)⟶1,\end{array}$$

(24)

which means that the term $|c^2{t}_0^{\prime \prime 2}|$ is insignificant compared to $2|c^2{t}_m^{\prime \prime }{t}_0^{\prime \prime }|$ .

Therefore these equations can be approximated as follows when the measurement errors are sufficiently small:

$$\begin{array}{l}{\parallel {\mathbf{p}}_m^{\prime \prime }\parallel }^2-c^2{t}_m^{\prime \prime 2}\approx 2{\mathbf{p}}_m^{\prime \prime T}{\mathbf{p}}_0^{\prime \prime }-2c^2{t}_m^{\prime \prime }{t}_0^{\prime \prime }\\ \phantom{{\parallel {\mathbf{p}}_m^{\prime \prime }\parallel }^2-c^2{t}_m^{\prime \prime 2}\approx }-2c^2{t}_m^{\prime \prime }{n}_m,\hspace{1em}m=1,\dots ,M\end{array}$$

(25)

or in vector-matrix form as follows:

$$\begin{array}{l}\left(\begin{array}{l}{\parallel {\mathbf{p}}_1^{\prime \prime }\parallel }^2-c^2{t}_1^{\prime \prime 2}\\ ⋮\\ {\parallel {\mathbf{p}}_M^{\prime \prime }\parallel }^2-c^2{t}_M^{\prime \prime 2}\end{array}\right)\\ \hspace{1em}\hspace{1em}\approx \left(\begin{array}{ll}2{\mathbf{p}}_0^{\prime \prime T}& -2c^2{t}_1^{\prime \prime }\\ ⋮& ⋮\\ 2{\mathbf{p}}_M^{\prime \prime T}& -2c^2{t}_M^{\prime \prime }\end{array}\right)\left(\begin{array}{l}{\mathbf{p}}_0^{\prime \prime }\\ \hspace{0.17em}\\ t_0^{\prime \prime }\end{array}\right)+\left(\begin{array}{l}-2c^2{t}_1^{\prime \prime }{n}_1\\ ⋮\\ -2c^2{t}_M^{\prime \prime }{n}_M\end{array}\right),\end{array}$$

(26)

or simply

$$\begin{array}{cc}& {\mathbf{b}}^{\mathrm{\prime \prime }}=\hspace{0.17em}{\mathbf{A}}^{\mathrm{\prime \prime }}{\mathbf{q}}^{\mathrm{\prime \prime }}\hspace{0.17em}+\hspace{0.17em}{\mathbf{w}}^{\mathrm{\prime \prime }}.\end{array}$$

(27)

With approximately known weights, the WLS solution can be written as follows:

$$\begin{array}{cc}& {\widehat{\mathbf{q}}}^{\mathrm{\prime \prime }}={\left({\mathbf{A}}^{\mathrm{\prime \prime }T}{{\mathbf{W}}^{\mathrm{\prime \prime }}}^{-1}{\mathbf{A}}^{\mathrm{\prime \prime }}\right)}^{-1}{\mathbf{A}}^{\prime T}{{\mathbf{W}}^{\mathrm{\prime \prime }}}^{-1}{\mathbf{b}}^{\prime }=\left(\begin{array}{c}{\mathbf{p}}_0^{\prime \prime }\\ \hspace{0.17em}\\ t_0^{\prime \prime }\end{array}\right),\end{array}$$

(28)

where

$$\begin{array}{cc}& {\mathbf{W}}^{\mathrm{\prime \prime }}=\mathrm{diag}\left\{\left[4c^4{t}_1^{\prime \prime 2},\dots ,4c^4{t}_M^{\prime \prime 2}\right]\right\}.\end{array}$$

(29)

After the aforementioned two-stage estimation, the corresponding coordinate system translation will be applied again to restore the original coordinate system:

$$\begin{array}{cc}& {\widehat{\mathbf{p}}}_{\mathbf{0}}={\widehat{\mathbf{p}}}_{\mathbf{0}}^{\mathrm{\prime \prime }}+{\widehat{\mathbf{p}}}_{\mathbf{0}}^{\prime }+{\mathbf{p}}_k,\end{array}$$

(30)

$$\begin{array}{cc}& {\widehat{t}}_0={\widehat{t}}_0^{\hspace{0.17em}\prime \prime }+{\widehat{t}}_0^{\hspace{0.17em}\prime }+t_k.\end{array}$$

(31)

We note that our adopted coordinate system translation can avoid any large or trivial absolute coordinates encountered (e.g., large values in positions, TOA measurements, and transmit time) that may result in the instability of the algorithm implementation.

In summary, we can list our proposed algorithm for the considered TOA-based source localization as below.

• (1)

Perform the first round of coordinate system translation according to (10)–(12) for $i=1,\dots ,M$ .

(2)

Perform the LS estimation with equal weights according to (16)–(18) in the first stage of estimation.

(3)

Perform the second round of coordinate system translation according to (19).

(4)

Perform the WLS estimation according to (25)–(28) in the second stage of estimation.

(5)

Perform the coordinate system restoration according to (30).

As for the computational complexity, we count it in terms of multiplications and additions as follows:

• (i)

Coordinate system translation for M sensors: $3M$ additions;

(ii)

LS estimation: around $15M$ multiplications and $14M$ additions;

(iii)

WLS estimation (diagonal weight matrix): around $18M$ multiplications and $16M$ additions.

Therefore, the total computational complexity should be around $33M$ multiplications and $36M$ additions.

Below, we propose two optional subsequent operations, for an improved performance or for the comparison later in this paper:

• (1)

Subsequent third stage of estimation: after the second stage of estimation, that is, the WLS estimation, perform the coordinate system translation similar to (19) and the WLS estimation similar to (25)–(28) and then perform the coordinate system restoration similar to (30). This is a three-stage algorithm.

(2)

Subsequent iterative LS minimization: an efficient version of such iterative LS minimization can be proposed as below:

$$\begin{array}{cc}& {\widehat{\mathbf{p}}}_0\left(k+1\right)=\frac{1}{M}\sum _{m=1}^M{\mathbf{p}}_{0m}\left(k⟶k+1\right),\end{array}$$

(32)

where the coordination can be updated as follows:

$$\begin{array}{cc}& {\mathbf{p}}_{0m}\left(k⟶k+1\right)={\mathbf{p}}_m+\frac{{\widehat{\mathbf{p}}}_0\left(k\right)-{\mathbf{p}}_m}{\parallel {\widehat{\mathbf{p}}}_0\left(k\right)-{\mathbf{p}}_m\parallel }\left[t_m-{\widehat{t}}_0\left(k\right)\right]c\end{array}$$

(33)

and the time update is completed by using the following:

$$\begin{array}{cc}& {\widehat{t}}_0\left(k+1\right)=\frac{1}{M}\sum _{m=1}^M{t}_{0m}\left(k⟶k+1\right),\end{array}$$

(34)

where

$$\begin{array}{cc}& t_{0m}\left(k⟶k+1\right)=t_m-\frac{1}{c}\parallel {\widehat{\mathit{p}}}_0\left(k\right)-{\mathit{p}}_m\parallel .\end{array}$$

(35)

5. Evaluation of Proposed Algorithm

5.1. Performance Analysis

In the previous section, we have explained the reasonable approximations, when measurement errors are sufficiently small, that result in the elimination of the nonlinearity from the converted observation equations, that is, in (25). If the nonlinearity has been eliminated, the corresponding relationship between measurements and basic unknowns in the converted observation equations would be linear. And the following theorem shows the relevant conditions for the existence of efficient unbiased estimator that can attain the CRLB in the presence of Gaussian measurement errors.

Theorem 1 (see [6]).

Suppose that $\mathrm{\Phi }$ is an open subset of $R^2$ and that $u\in C(\varphi )$ . Let

$$\begin{array}{cc}& \mathbf{r}\hspace{0.17em}\mathbf{=}\hspace{0.17em}\mathbf{u}\left(\mathbf{s}\right)\mathbf{+}\hspace{0.17em}\mathbf{w}\mathbf{,}\end{array}$$

(36)

and suppose that w is a Gaussian vector. Then an efficient unbiased estimator of s exists if and only if

$$\begin{array}{cc}& \mathbf{u}\left(\mathbf{s}\right)\mathbf{=}\hspace{0.17em}\mathbf{H}\mathbf{s}\hspace{0.17em}\mathbf{+}\hspace{0.17em}\mathbf{v}\end{array}$$

(37)

or some full-rank matrix H and some vector v.

See proof in [6].

According to this theorem, we know that our proposed algorithm for the considered TOA-based source localization can achieve the CRLB in the presence of sufficiently small zero-mean Gaussian independent measurement errors, as the linearization adopting suitable coordinate system translation has eliminated the nonlinearity from the converted observation equations.

5.2. Simulation Studies

We evaluate, by simulations, our proposed algorithm for the considered TOA-based source localization. We assume that the source and all five sensors are randomly deployed in a square area cornered by $(-300,-300)$ , $(-300,300)$ , $(300,300)$ ,  and $(300,-300)$ (in meters) and that the source transmit time is randomly set in a time period $(-1000,1000)$ (in nanoseconds). The TOA measurement errors are assumed as zero-mean Gaussian errors and independent of each other.

We compare three versions of our proposed algorithm, that is, the one-stage algorithm (only containing the first stage of coordinate system translation and estimation as well as the subsequent coordinate system restoration), the two-stage algorithm (containing two stages), and the three-stage algorithm (containing three stages, i.e., repeating the second stage), the two-stage algorithm with subsequent iterative LS minimization (converged performance), with corresponding CRLB.

Figure 1 shows the results when the source-sensor deployment and the source transmit time are randomly generated and then fixed, whereas Figure 2 shows the results for randomly generated source-sensor deployment and transmit time. The curves are obtained by averaging over 1000 independent runs.

Figure 1

The simulated localization error against the measurement error: the source-sensor deployment and source transmit time are randomly generated and then fixed. The curves are obtained by averaging over 1000 independent runs, in each of which all measurement errors are generated independently.

Figure 2

The simulated localization error against the measurement error: the source-sensor deployment and source transmit time are randomly generated. The curves are obtained by averaging over 1000 independent runs, in each of which the measurement errors, source-sensor deployment and source transmit time are generated independently.

From the figures we observe that the two-stage version of our proposed algorithm has better localization performance than the one-stage algorithm and can achieve the CRLB when the measurement errors are sufficiently small (up to a considerable large range). The three-stage algorithm achieves the localization performance that is only slightly better than the two-stage algorithm when the measurement errors are relatively larger. The two-stage algorithm with subsequent iterative LS minimization has a converged performance also slightly better than the two-stage algorithm. From the figures, we also can see that, in the upper range of the measure errors, the localization error of our proposed two-stage algorithm are away from the CRLB. This is mainly due to two reasons: the condition for the approximation in (25) does not hold as the measurement errors are comparable to the source-sensor distances, the derived CRLB based on pdf only holds in the vicinity of the correct values where the measurement errors should be sufficiently small.

6. Conclusion

In this paper, we have proposed a novel algorithm for the considered TOA-based source localization. The proposed algorithm utilizes linearization and adopts suitable coordinate system translation that can eliminate nonlinearity from the converted observation equations, thus, can achieve the CRLB in the presence of sufficiently small zero-mean Gaussian independent measurement errors. This is shown by the performance analysis and simulation study.