Fixed Point Iteration Based Algorithm for Asynchronous TOA-Based Source Localization

Abstract

This paper investigates the problem of source localization using signal time-of-arrival (TOA) measurements in the presence of unknown start transmission time. Most state-of-art methods are based on convex relaxation technologies, which possess global solution for the relaxed optimization problem. However, computational complexity of the convex optimization–based algorithm is usually large, and need CVX toolbox to solve it. Although the two stage weighted least squares (2SWLS) algorithm has very low computational complexity, its estimate performance is susceptible to sensor geometry and threshold phenomenon. A new algorithm that is directly derived from maximum likelihood estimator (MLE) is developed. The newly proposed algorithm is named as fixed point iteration (FPI); it only involves simple calculations, such as addition, multiplication, division, and square-root. Unlike state-of-the-art methods, there is no matrix inversion operation and can avoid the unstable performance incurred by singular matrix. The FPI algorithm can be easily extended to the scenario with sensor position errors. Finally, simulation results demonstrate that the proposed algorithm reaches a good balance between computational complexity and localization accuracy.

1. Introduction

Source localization has a wide range of applications, such as navigation, monitoring, tracking, rescue and so on. There are lots of methods can be used to determine the source position, where time-of-arrival (TOA), time-difference-of-arrival (TDOA), time-delay (TD), received signal strength (RSS), angle-of-arrival (AOA) and their combinations are commonly positioning schemes [1,2,3,4,5,6,7,8].

A TOA-based localization system requires that the source and sensors be accurately time synchronized. However, it is hard to maintain the accurate synchronization when the source is non-cooperative [9]. The pseudo ranges are corrupted with an unknown distance when the start transmission time is unknown in the asynchronous networks. As a result, study on the asynchronous TOA (ATOA)-based localization problem is very important in the practice.

In the work reported in [10,11,12,13,14,15,16,17,18], the ATOA-based localization problem can be addressed by two ways: (a) convex optimization methods, including semidefinite programming (SDP) and second-order-cone programming (SOCP) [10,11,12,13,14,15,16,17]; (b) two stage weighted least squares (2SWLS) method [18].

Convex optimization methods have the merit of global solution, but their computational complexity is usually huge and need a CVX toolbox to solve the optimization problem [19]. For example, Xu et al. [10] developed two SDP-based algorithms corresponding to two step least squares (2LS) and min–max (MMA) criterion, respectively. Then in [14], Zou et al. showed that the original 2LS-based SDP algorithm cannot provide a good solution because the weighted matrix $\mathbf{G}$ is singular. In order to improve the tightness of the original SDP algorithm, second-order-cone constraints and penalty term are jointly added to the original SDP algorithm. However, the algorithm proposed in [14] needs to calculate SDP problem multiple times and results in great computational burden. Recently, in [17], Ma et al. formulated an SDP algorithm to address the constraint existing in the WLS problem, which was formulated by Huang et al. in [18].

In [18], Huang et al. proposed a 2SWLS method to jointly estimate the start transmission time and the location of source. The 2SWLS method has very low computational complexity, but its localization accuracy is highly susceptible to sensor geometry when it is uniform circular array and the source is close to the array center. Besides this, the 2SWLS method presents a threshold phenomenon.

Another challenging issue for ATOA localization is sensor position errors, which could drastically degrade the accuracy [20]. Several works have addressed this problem [10,13,14,15,16], which have shown that using convex optimization could effectively address the sensor position errors. However, their computational complexity is also large.

In this paper, we develop a fixed point iteration (FPI) algorithm for the ATOA-based localization problem. The FPI algorithm is entirely and directly developed from maximum likelihood estimator (MLE). The FPI solution satisfies the equation that gradient is equal to zero. It is different from state-of-the-art methods. First, unlike weighted least squares (WLS) based algorithms, there is no approximation in its derivation. Second, unlike semidefinite programming (SDP) based algorithms, there is also no convex relaxation in its derivation. Third, unlike Gauss–Newton or Quasi-Newton based algorithms [21,22], there is no first-order Taylor series expansion in its derivation. Besides this, the proposed FPI algorithm can be easily tailored to the scenario with sensor position errors.

The rest of this paper is organized as follows. Section 2 develops an iterative algorithm when sensors positions are accurate. The algorithm is then extended to the scenario with non-accurate sensors positions in Section 3. Simulation results are presented in Section 4 to compare the performances of the proposed algorithm, state-of-art algorithms and Cramér-Rao lower bound (CRLB).

The following notations are used throughout the paper. Bold lowercase and uppercase letters denote vectors and matrices, respectively; $\mathbf{A}(:,i)$ denotes the ith column of matrix $\mathbf{A}$, and $\mathbf{A}(i,:)$ denotes the ith row of matrix $\mathbf{A}$; $\mathrm{tr}\left(\mathbf{A}\right)$ is the trace of $\mathbf{A}$; ${\mathbf{I}}_M$ is the $M\times M$ identity matrix, ${\mathbf{1}}_M$ is the column vector of M ones, and ${\mathbf{0}}_M$ is the column vector of M zeros; $∥·∥$ is the $l_2$ norm; ⊗ is the Kroneker product of two matrices; and $\mathbf{E}[·]$ is the expectation.

2. Localization with Accurate Sensor Positions

Consider a network with M sensors and one source whose location $\mathbf{u}\in {\mathbb{R}}^{m\times 1}$ is unknown and to be estimated ($m=2$ or 3). The ith sensor position ${\mathbf{s}}_i$ is accurate known.

Under line-of-sight condition, the range-of-arrival (ROA) measurements between sensor i and the source are expressed as

$$r_i=t_{0}c+∥\mathbf{u}-{\mathbf{s}}_i∥+n_i,i=1,\cdots ,M$$

(1)

where c is the known signal propagation speed, and $t_0$ is unknown start transmission time. In the next, ${\mathbf{t}}_{0}c$ will be replaced by $d_0$, and Equation (1) can be written as

$$r_i=d_0+∥\mathbf{u}-{\mathbf{s}}_i∥+n_i$$

(2)

where $n_i,i=1,\cdots ,M,$ are the ROA measurement noise, which are modeled as zero-mean Gaussian random variables with covariance matrix $\mathbf{Q}=\mathrm{diag}\left([{\sigma }_1^2,\dots ,{\sigma }_M^2]\right)$. The maximum likelihood estimator (MLE) from Equation (2) is formulated as

$$\begin{array}{cc}\hfill & \underset{\begin{array}{c}\mathbf{u},d_0\end{array}}{min}\sum _{i=1}^M\frac{{\left(r_i-d_0-∥\mathbf{u}-{\mathbf{s}}_i∥\right)}^2}{{\sigma }_i^2}.\hfill \end{array}$$

(3)

The above formulation can also be written as

$$\begin{array}{cc}\hfill \underset{\begin{array}{c}\mathbf{u},d_0,\mathbf{d}\end{array}}{min}& {(\mathbf{r}-d_0{\mathbf{1}}_M-\mathbf{d})}^T{\mathbf{Q}}^{-1}(\mathbf{r}-d_0{\mathbf{1}}_M-\mathbf{d})\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& d_i=∥\mathbf{u}-{\mathbf{s}}_i∥.\hfill \end{array}$$

(4b)

where $\mathbf{d}={[d_1,\dots ,d_M]}^T$. Let gradient of objective function in (4a) with respect to $d_0$ to zero

$$\begin{array}{c}\hfill -2{\mathbf{1}}_M^T{\mathbf{Q}}^{-1}(\mathbf{r}-d_0{\mathbf{1}}_M-\mathbf{d})=\mathbf{0}\end{array}$$

(5)

and we can obtain

$$\begin{array}{c}\hfill d_0=\frac{{\mathbf{1}}_M^T{\mathbf{Q}}^{-1}(\mathbf{r}-\mathbf{d})}{{\mathbf{1}}_M^T{\mathbf{Q}}^{-1}{\mathbf{1}}_M}.\end{array}$$

(6)

Putting above $d_0$ back to Equation (4), we obtain

$$\begin{array}{cc}\hfill \underset{\begin{array}{c}\mathbf{u},\mathbf{d}\end{array}}{min}& {(\mathbf{r}-\mathbf{d})}^T\mathbf{G}(\mathbf{r}-\mathbf{d})\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& d_i=∥\mathbf{u}-{\mathbf{s}}_i∥\hfill \end{array}$$

(7b)

where

$$\begin{array}{c}\hfill \mathbf{G}={\mathbf{I}}_M-\frac{{\mathbf{1}}_M{\mathbf{1}}_M^T{\mathbf{Q}}^{-1}}{{\mathbf{1}}_M^T{\mathbf{Q}}^{-1}{\mathbf{1}}_M}.\end{array}$$

(8)

Equation (7) is a non-convex optimization problem, and it can be resorted to convex relaxation method [14]. However, the convex optimization–based algorithm usually has very high computational complexity, and need CVX toolkit to solve it. Next, we develop an alternative solution for the above problem. First, Equation (7) can be expressed as the form of unconstrained optimization problem

$$\begin{array}{c}\hfill \underset{\begin{array}{c}\mathbf{u}\end{array}}{min}\sum _{i=1}^M\sum _{j=1}^{M}G(i,j)(r_i-∥\mathbf{u}-{\mathbf{s}}_i∥)(r_j-∥\mathbf{u}-{\mathbf{s}}_j∥).\end{array}$$

(9)

The gradient of objective function in (9) with respect to $\mathbf{u}$ is

$$\begin{array}{c}\hfill \mathbf{g}=-2\sum _{i=1}^M\sum _{j=1}^{M}G(i,j)(r_i-∥\mathbf{u}-{\mathbf{s}}_i∥)\frac{\mathbf{u}-{\mathbf{s}}_j}{∥\mathbf{u}-{\mathbf{s}}_j∥}.\end{array}$$

(10)

The necessary condition for the optimization problem in (9) is

$$\begin{array}{c}\hfill \mathbf{g}=0\end{array}$$

(11)

i.e.,

$$\begin{array}{cc}\hfill & \sum _{i=1}^{M}G(i,i)\left(r_i\frac{\mathbf{u}-{\mathbf{s}}_i}{∥\mathbf{u}-{\mathbf{s}}_i∥}-(\mathbf{u}-{\mathbf{s}}_i)\right)+\sum _{i=1}^M\sum _{j\ne i}^{M}G(i,j)(r_i-∥\mathbf{u}-{\mathbf{s}}_i∥)\frac{\mathbf{u}-{\mathbf{s}}_j}{∥\mathbf{u}-{\mathbf{s}}_j∥}=0.\hfill \end{array}$$

(12)

Next, move the linear function of $\mathbf{u}$ to one side of the equation, and the constant term and nonlinear function of $\mathbf{u}$ to another side of the equation

$$\begin{array}{cc}\hfill & \sum _{i=1}^{M}G(i,i)\mathbf{u}=\sum _{i=1}^{M}G(i,i)\left(r_i\frac{\mathbf{u}-{\mathbf{s}}_i}{∥\mathbf{u}-{\mathbf{s}}_i∥}+{\mathbf{s}}_i\right)+\sum _{i=1}^M\sum _{j\ne i}^{M}G(i,j)(r_i-∥\mathbf{u}-{\mathbf{s}}_i∥)\frac{\mathbf{u}-{\mathbf{s}}_j}{∥\mathbf{u}-{\mathbf{s}}_j∥}.\hfill \end{array}$$

(13)

Then

$$\begin{array}{cc}\hfill & \mathbf{u}=\frac{1}{\mathrm{tr}\left(\mathbf{G}\right)}\left(\sum _{i=1}^{M}G(i,i)\left(r_i\frac{\mathbf{u}-{\mathbf{s}}_i}{∥\mathbf{u}-{\mathbf{s}}_i∥}+{\mathbf{s}}_i\right)+\sum _{i=1}^M\sum _{j\ne i}^{M}G(i,j)(r_i-∥\mathbf{u}-{\mathbf{s}}_i∥)\frac{\mathbf{u}-{\mathbf{s}}_j}{∥\mathbf{u}-{\mathbf{s}}_j∥}\right).\hfill \end{array}$$

(14)

Moreover,

$$\begin{array}{c}\hfill \mathrm{tr}\left(\mathbf{G}\right)=\mathrm{tr}({\mathbf{I}}_M-\frac{1_M{1}_M^T{\mathbf{Q}}^{-1}}{1_M^T{\mathbf{Q}}^{-1}{1}_M})=M-1.\end{array}$$

(15)

As a result, the FPI computation of $\mathbf{u}$ is

$$\begin{array}{cc}\hfill & {\mathbf{u}}_k=f\left({\mathbf{u}}_{k-1}\right)=\frac{1}{M-1}\left(\sum _{i=1}^{M}G(i,i)\left(r_i\frac{{\mathbf{u}}_{k-1}-{\mathbf{s}}_i}{∥{\mathbf{u}}_{k-1}-{\mathbf{s}}_i∥}+{\mathbf{s}}_i\right)+\sum _{i=1}^M\sum _{j\ne i}^{M}G(i,j)(r_i-∥{\mathbf{u}}_{k-1}-{\mathbf{s}}_i∥)\frac{{\mathbf{u}}_{k-1}-{\mathbf{s}}_j}{∥{\mathbf{u}}_{k-1}-{\mathbf{s}}_j∥}\right).\hfill \end{array}$$

(16)

It should be noted that the above derivations are based on the necessary condition. In other words, the result from above FPI may be a stationary point, i.e., a local maximum, a local minimum point or a saddle point. In order to obtain the global minimum point, we adopt the following steps. First is to use multiple initial values to calculate FPI in Equation (16). Without loss of generality, the set of initial values can be set as follows (from (16), it can be seen that the initial values should be not equal to the location of sensors. The main principle for the setting of initial values is that it includes points located inside and outside the convex hull of sensors).

$$\begin{array}{c}\hfill {\mathbf{U}}_0=[\frac{1}{M}\sum _{i=1}^M{\mathbf{s}}_i,2{\mathbf{s}}_1-\frac{1}{M}\sum _{i=1}^M{\mathbf{s}}_i,\dots ,2{\mathbf{s}}_M-\frac{1}{M}\sum _{i=1}^M{\mathbf{s}}_i].\end{array}$$

(17)

Then, we calculate the corresponding cost function for each initial value

$$\begin{array}{c}\hfill c_f\left(n\right)={(\mathbf{r}-{\mathbf{d}}_n^{*})}^T\mathbf{G}(\mathbf{r}-{\mathbf{d}}_n^{*}),n=1,\dots ,M+1\end{array}$$

(18)

where ${\mathbf{d}}_n^{*}={[∥{\mathbf{u}}_n^{*}-{\mathbf{s}}_1∥,\dots ,∥{\mathbf{u}}_n^{*}-{\mathbf{s}}_M∥]}^T$. Finally, we choose the result that minimizes the cost function.

Besides this, in order to reduce the calculation amount, we take the following measure. If $c_f\left(n\right)<\mu$ for a certain n, then we stop calculating the iterations for the rest initial values. Next, we show how to determine the threshold $\mu$. Reviewing (3), we can obtain

$$\begin{array}{c}\hfill \mathbf{E}[\sum _{i=1}^M\frac{{\left(r_i-d_0-∥\mathbf{u}-{\mathbf{s}}_i∥\right)}^2}{{\sigma }_i^2}]=\mathbf{E}[\sum _{i=1}^M\frac{n_i^2}{{\sigma }_i^2}]=M\end{array}$$

(19)

and

$$\begin{array}{c}\hfill \mathbf{E}\left[{(\mathbf{r}-\mathbf{d})}^T\mathbf{G}(\mathbf{r}-\mathbf{d})\right]=M\end{array}$$

(20)

for true source location $\mathbf{u}$. As a result, the threshold is set to $\mu =M$.

The above proposed algorithm procedures are summarized in Algorithm 1.

Algorithm 1: Localization with scenario 1: accurate sensor position.

3. Localization with Non-Accurate Sensor Positions

In this section, we will extend the FPI localization algorithm to the scenario with non-accurate sensor positions. First, the ith true but unknown sensor location is ${\mathbf{s}}_i^0$. The available sensor location can be denoted as

$${\mathbf{s}}_i={\mathbf{s}}_i^0+{\beta }_i,i=1,\cdots ,M$$

(21)

where ${\beta }_i$ is a zero-mean white Gaussian vector with covariance matrix ${\delta }_i^2{\mathbf{I}}_m$ [23]. The MLE is formulated as

$$\begin{array}{cc}\hfill & \underset{\begin{array}{c}\mathbf{u},d_0,{\mathbf{s}}_i^0\end{array}}{min}\sum _{i=1}^M\frac{{\left(r_i-d_0-∥\mathbf{u}-{\mathbf{s}}_i^0∥\right)}^2}{{\sigma }_i^2}+\sum _{i=1}^M\frac{∥{\mathbf{s}}_i-{\mathbf{s}}_i^0∥}{{\delta }_i^2}.\hfill \end{array}$$

(22)

Putting (6) into (22), this results in

$$\begin{array}{cc}\hfill \underset{\begin{array}{c}\mathbf{u},\mathbf{d},{\mathbf{s}}_i^0\end{array}}{min}& {(\mathbf{r}-\mathbf{d})}^T\mathbf{G}(\mathbf{r}-\mathbf{d})+\sum _{i=1}^M\frac{{∥{\mathbf{s}}_i-{\mathbf{s}}_i^0∥}^2}{{\delta }_i^2}\hfill \end{array}$$

(23a)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& d_i=∥\mathbf{u}-{\mathbf{s}}_i^0∥.\hfill \end{array}$$

(23b)

Let $\mathbf{x}=[\mathbf{u};{\mathbf{s}}_1^0;\dots ;{\mathbf{s}}_M^0]$, $\mathbf{s}=[{\mathbf{s}}_1;\dots ;{\mathbf{s}}_M]$, ${\mathbf{W}}_s={\mathbf{Q}}_s\otimes {\mathbf{I}}_m$, ${\mathbf{Q}}_s=\mathrm{diag}\left([\frac{1}{{\delta }_1^2},\dots ,\frac{1}{{\delta }_M^2}]\right)$, ${\mathbf{P}}_i=[{\mathbf{I}}_m,{\mathbf{0}}_{m,m(i-1)},-{\mathbf{I}}_m,{\mathbf{0}}_{m,m(M-i)}]$, $\mathbf{B}=[{\mathbf{0}}_{mM,m},{\mathbf{I}}_{mM}]$. Equation (23) can be written as

$$\begin{array}{cc}\hfill \underset{\begin{array}{c}\mathbf{x}\end{array}}{min}& \sum _{i=1}^M\sum _{j=1}^{M}G(i,j)(r_i-∥{\mathbf{P}}_i\mathbf{x}∥)(r_j-∥{\mathbf{P}}_j\mathbf{x}∥)+\sum _{i=1}^{mM}{W}_s(i,i){\left(s\left(i\right)-\mathbf{B}(i,:)\mathbf{x}\right)}^2.\hfill \end{array}$$

(24)

The gradient of objective function in Equation (24) with respect to $\mathbf{x}$ is

$$\begin{array}{cc}\hfill \mathbf{g}=& -2\sum _{i=1}^M\sum _{j=1}^{M}G(i,j)(r_i-∥{\mathbf{P}}_i\mathbf{x}∥)\frac{{\mathbf{P}}_j^T{\mathbf{P}}_j\mathbf{x}}{∥{\mathbf{P}}_j\mathbf{x}∥}-2\sum _{i=1}^{mM}{W}_s(i,i)\left(s\left(i\right)-\mathbf{B}(i,:)\mathbf{x}\right)\mathbf{B}{(i,:)}^T.\hfill \end{array}$$

(25)

The necessary condition for the optimization problem in (24) is

$$\begin{array}{c}\hfill \mathbf{g}=0.\end{array}$$

(26)

Equivalent to

$$\begin{array}{cc}\hfill & \sum _{i=1}^{M}G(i,i)r_i\frac{{\mathbf{P}}_i^T{\mathbf{P}}_i\mathbf{x}}{∥{\mathbf{P}}_i\mathbf{x}∥}-\sum _{i=1}^{M}G(i,i){\mathbf{P}}_i^T{\mathbf{P}}_i\mathbf{x}+\sum _{i=1}^M\sum _{j\ne i}^{M}G(i,j)(r_i-∥{\mathbf{P}}_i\mathbf{x}∥)\frac{{\mathbf{P}}_j^T{\mathbf{P}}_j\mathbf{x}}{∥{\mathbf{P}}_j\mathbf{x}∥}+\hfill \\ & \sum _{i=1}^{mM}{W}_s(i,i)s\left(i\right)\mathbf{B}{(i,:)}^T-\sum _{i=1}^{mM}{W}_s(i,i)\mathbf{B}{(i,:)}^T\mathbf{B}(i,:)\mathbf{x}=0\hfill \end{array}$$

(27)

i.e.,

$$\begin{array}{cc}\hfill & \left(\sum _{i=1}^{M}G(i,i){\mathbf{P}}_i^T{\mathbf{P}}_i+\sum _{i=1}^{mM}{W}_s(i,i)\mathbf{B}{(i,:)}^T\mathbf{B}(i,:)\right)\mathbf{x}=\sum _{i=1}^{M}G(i,i)r_i\frac{{\mathbf{P}}_i^T{\mathbf{P}}_i\mathbf{x}}{∥{\mathbf{P}}_i\mathbf{x}∥}+\sum _{i=1}^M\sum _{j\ne i}^{M}G(i,j)(r_i-∥{\mathbf{P}}_i\mathbf{x}∥)\frac{{\mathbf{P}}_j^T{\mathbf{P}}_j\mathbf{x}}{∥{\mathbf{P}}_j\mathbf{x}∥}\hfill \\ \hfill & +\sum _{i=1}^{mM}{W}_s(i,i)s\left(i\right)\mathbf{B}{(i,:)}^T.\hfill \end{array}$$

(28)

Finally,

$$\begin{array}{cc}\hfill & \mathbf{x}={\left(\sum _{i=1}^{M}G(i,i){\mathbf{P}}_i^T{\mathbf{P}}_i+\sum _{i=1}^{mM}{W}_s(i,i)\mathbf{B}{(i,:)}^T\mathbf{B}(i,:)\right)}^{-1}·\hfill \\ \hfill & [\sum _{i=1}^{M}G(i,i)r_i\frac{{\mathbf{P}}_i^T{\mathbf{P}}_i\mathbf{x}}{∥{\mathbf{P}}_i\mathbf{x}∥}+\sum _{i=1}^M\sum _{j\ne i}^{M}G(i,j)(r_i-∥{\mathbf{P}}_i\mathbf{x}∥)\frac{{\mathbf{P}}_j^T{\mathbf{P}}_j\mathbf{x}}{∥{\mathbf{P}}_j\mathbf{x}∥}+\sum _{i=1}^{mM}{W}_s(i,i)s\left(i\right)\mathbf{B}{(i,:)}^T].\hfill \end{array}$$

(29)

The FPI computation of $\mathbf{x}$ is

$$\begin{array}{cc}\hfill & {\mathbf{x}}_k={\left(\sum _{i=1}^{M}G(i,i){\mathbf{P}}_i^T{\mathbf{P}}_i+\sum _{i=1}^{mM}{W}_s(i,i)\mathbf{B}{(i,:)}^T\mathbf{B}(i,:)\right)}^{-1}·\hfill \\ \hfill & [\sum _{i=1}^{M}G(i,i)r_i\frac{{\mathbf{P}}_i^T{\mathbf{P}}_i{\mathbf{x}}_{k-1}}{∥{\mathbf{P}}_i{\mathbf{x}}_{k-1}∥}+\sum _{i=1}^{mM}{W}_s(i,i)s\left(i\right)\mathbf{B}{(i,:)}^T+\sum _{i=1}^M\sum _{j\ne i}^{M}G(i,j)(r_i-∥{\mathbf{P}}_i{\mathbf{x}}_{k-1}∥)\frac{{\mathbf{P}}_j^T{\mathbf{P}}_j{\mathbf{x}}_{k-1}}{∥{\mathbf{P}}_j{\mathbf{x}}_{k-1}∥}].\hfill \end{array}$$

(30)

The proposed algorithm for scenario 2 is shown in Algorithm 2. It includes two stages. The first stage is to carry out Algorithm 1 to obtain an initial guess of source position by using the non-accurate sensor positions. The second stage is to update $\mathbf{x}$ using Equation (30) and taking the result from stage one as initial value.

Algorithm 2: Localization with scenario 2: non-accurate sensor positions.

4. Simulation Results

In this section, two scenarios with accurate and non-accurate sensor positions are respectively considered to demonstrate the performance of two proposed algorithms.

In the first scenario, we compare several algorithms: the proposed Algorithm 1 (with one middle point as initial value labeled as ‘Proposed-Mid’, with multiple initial values labeled as ‘Proposed-Multi’), the Gauss–Newton method (in the simulation, the result of Gauss–Newton method may be NaN due to the matrix is singular in the calculation of matrix inversion. In order to avoid this case, it is set that, if the distance between the updated estimate and the initial value is large than 1000, then the iteration will be stopped and output the initial value as its result) with one middle point as initial value labeled as ‘Gauss–Newton-Mid’, with multiple initial values labeled as ‘Gauss–Newton-Multi’; the derivation for the Gauss–Newton method is shown in Appendix A), 2SWLS [18], SDP-Zou [14], SDP–Ma [17] and CRLB.

In the second scenario, we compare the proposed Algorithm 2 (labeled as ‘Proposed’) with the Gauss–Newton method (labeled as ‘Gauss–Newton’), SDP–Zou [14] and CRLB when sensor positions are non-accurate. The derivation of the Gauss–Newton method with position errors for the problem of (24) is given in Appendix B.

The setting of parameters are the same as in [14]: a network with four sensors, and one source is simulated. The positions of the sensors are ${[-10,-10]}^{T}m$, ${[-10,10]}^{T}m$, ${[10,-10]}^{T}m$, ${[10,10]}^{T}m$, and the location of source is randomly chosen from a square $[-15,15]m\times [-15,15]m$. The start transmission time $t_0$ is drawn from uniform distribution $\mathrm{U}[10,40]ns$. Sensors’ position covariance matrix is ${\delta }^2{\mathbf{I}}_2$. The range measurement covariance matrix is $\mathbf{Q}={\sigma }^2{\mathbf{I}}_M$. The initial set for Algorithm 1 is ${\mathbf{U}}_0=[0,-20,-20,20,20;0,-20,20,$$-20,20]$. Five penalty factors ${\eta }_1,{\eta }_2,{\eta }_3,{\eta }_4,{\eta }_5$ equal to ${10}^{-4}$, ${10}^{-3}$, ${10}^{-2}$, ${10}^{-1}$, ${10}^0$ are used for the SDP–Zou algorithm. The SDP–Zou and SDP–Ma algorithms are implemented by CVX toolbox using SeDuMi as a solver and with best precision [24]. The Gauss–Newton methods and the proposed FPI algorithm have same initials. Root mean square errors (RMSEs) are drawn from 4000 Monte Carlo realizations in the following simulations.

From Figure 1, we have following observations: (1) only the proposed algorithm with multiple initial values and the SDP–Zou algorithm can reach CRLB; (2) the ‘Proposed-Multi’ algorithm is superior to the ‘Proposed-Mid’ algorithm, which validates the existence of other stationary point except for the global minimum point in the problem of MLE; (3) the Gauss–Newton methods have inferior performance due to its first-order approximation in the derivations; (4) the 2SWLS algorithm has the worst performance even with small noise condition, which is due to the sensor geometry is uniform circular array and the source is close to the array center occasionally [25]; (5) the SDP–Ma algorithm has poor performance, which is due to the number of sensors, namely, 4, which results in rank-1 solution being hard to obtain [17].

Besides this, the average running times of different algorithms are given in

Table 1. The average running time [ms] of the considered algorithms with the two scenarios, in which a star denotes that the algorithm is not fitting for a certain scenario. CPU: i7-8700 3.2 GHz.

	Scenarios	1	2
Algorithms
2SWLS		0.19	*
Proposed-Multi		8.6	*
Gauss–Newton-Multi		0.94	*
Proposed-Mid		5.3	*
Gauss–Newton-Mid		0.67	*
SDP–Ma		412	*
SDP–Zou		1773	1861
Proposed		*	28.6
Gauss–Newton		*	2.6

. It can be seen that: (1) the 2SWLS algorithm has the lowest computation time, because it only involves three times computation of WLS solution; (2) both SDP–Zou and SDP–Ma algorithms are time-consuming to obtain the solutions by CVX toolbox; moreover, the SDP–Zou algorithm needs to solve the SDP problem five times; (3) the average running time of ‘Proposed-Multi’ is large than but less than twice of the ‘Proposed-Mid’, which validates that the stop criterion in (20) can help to reduce the calculation amount. In the simulation, we found that if $\mathsf{\mu }=2M$ the calculation amount can be further reduced; (4) the Gauss–Newton algorithms are faster than the proposed FPI algorithms.

Figure 2 and Figure 3, respectively, evaluate the source and sensors estimate performances of different algorithms versus $\delta$ when the locations of sensors are non-accurate. From these two figures, it is observed that the performance of the proposed algorithm is superior to the other two algorithms. Besides this, in Table 1, it can be seen that the average running time of the proposed algorithm is much less than the SDP–Zou algorithm, but greater than the Gauss–Newton algorithm.

Finally, from the above results and analyses, we can obtain a conclusion: The proposed FPI algorithms reach a good balance between estimation accuracy and calculation amount.

5. Conclusions

Two FPI algorithms for ATOA-based localization problem are respectively developed corresponding to the scenarios with accurate and non-accurate sensor positions. First, the problem of MLE for the unknown start transmission time and source position is formulated. Next, the gradient of objective function in MLE is derived, then we simplify the equation with gradient equal to zero. Finally, we extract the linear variable to one side of the equation and obtain the FPI algorithm. The simulation results validate the performances of the two proposed algorithms in terms of accuracy and computation.