Refinement of TOA Localization with Sensor Position Uncertainty in Closed-Form

Abstract

The subject of localization has received great deal attention in the past decades. Although it is perhaps a well-studied problem, there is still room for improvement. Traditional localization methods usually assume the number of sensors is sufficient for providing desired performance. However, this assumption is not always satisfied in practice. This paper studies the time of arrival (TOA)-based source positioning in the presence of sensor position errors. An error refined solution is developed for reducing the mean-squared-error (MSE) and bias in small sensor network (the number of sensors is fewer) when the noise or error level is relatively large. The MSE performance is analyzed theoretically and validated by simulations. Analytical and numerical results show the proposed method attains the Cramér-Rao lower bound (CRLB). It outperforms the existing closed-form methods with slightly raising computation complexity, especially in the larger noise/error case.

1. Introduction

Location determination is one of the classical research fields in signal processing. With the rise of the fifth-generation communication system (5G), Internet of things, automatic pilot and unmanned aerial vehicle, localization continues to receive great attention [1,2,3,4]. Based on different type of measurements, the common localization methods estimate the source position using time of arrival (TOA), time difference of arrival (TDOA), angle of arrival (AOA), Doppler frequency, Fingerprint and received signal strength (RSS) [5,6,7,8,9,10,11,12,13]. To improving the accuracy further, a series of hybrid method is developed by combining several kinds of measurements [14,15,16]. Since TOA localization has the precision advantage in position estimate, particularly in indoor environments [17], it attracts much interest in the netted mono-static radar, multi-static radar and distributed multi-input multi-output (MIMO) system.

TOA localization is categorized into two classes: circular-TOA method and elliptic-TOA method [18]. The circular-TOA localization [19] is usually applied to estimate the emitter location whose signal is intercepted by passive sensors. The TOAs (or ranges if the propagation speed is known) determine circles, taking the sensor as the center, that trace out the possible source location. The elliptic-TOA localization [20,21] creates ellipse loci through the signals traveling times from transmitters to the target and reflected back to received sensors. The ellipses from different sensor pairs intersect, yielding an estimate of the target position. In this paper, we focus our study on the circular-TOA case.

Time synchronization error is one of the factors that degrade the TOA localization performance. An opportunistic positioning method is proposed, exploiting radio transmitters and GPS-equipped nodes jointly [22]. G. Wang [23] considers the asynchronous sensors that the start transmission time is unknown. The source location is solved through second-order cone relaxation (SOCR). Considering the synchronization and localization problem simultaneously, the cooperation improves localization performance significantly without the requirement of high anchor node densities [24]. Y. Zou [25] presents an improved semidefinite programming (SDP) method when both asynchronism and sensor uncertainty occurs. Transforming the TOA measurements into TDOA measurements may be an effective way to eradicate the clock offset of the target’s clock [26]. Using the a priori of the existing unknown clock skew, the source position is solved by a fractional programming problem, resulting in the superior performance of the state-of-the-art methods. Y. Kang [27] proposes a technique to eliminate the time synchronization error iteratively, achieving a position estimator almost without time asynchronization. To reduce the computational cost, by applying robust squared-range (R-SR) and weighted least squares (WLS) criteria, [28] converts the originally nonconvex problem into a generalized trust region subproblem (GTRS) framework. It realizes comparable performance to the current methods with significantly higher effectiveness. To deal with multiple types of errors, the artificial neural network (ANN) and radial basis function (RBF) neural network are introduced to alleviate the accuracy loss due to the measurement error, non-line of sight (NLOS) error, and synchronization error [29].

Another factor that reduces the position estimation accuracy is sensor position uncertainty. Many works of literature study the source location estimate problem without considering the sensor position errors [30]. The research in [31] discusses the effect of sensor position errors. It indicates that using statistical knowledge of the sensor position errors is able to achieve the estimator with optimum localization performance. Commonly, there are two ways to take advantage of this statistical knowledge. On the one hand, the errors are estimated in terms of the coarse source position obtained under the assumption that sensor positions are precise. After calibrating the sensors, a more accurate source position estimate is achieved. This method is more efficient for locating multiple disjoint sources [32]. On the other hand, the sensor position error can be equivalently transformed into the measurement noise, where a new weighting matrix is constructed to cancel out the effect of sensor errors [33,34,35].

Usually, the Maximum Likelihood (ML) approaches which perform exhaustive grid search or iteration, and the SDP-based solutions are computation-cumbersome. Due to the low computation cost and easy implementation in engineering, the closed-form approaches are more attractive. The two-stage method, which is perhaps the most classical closed-form solution, resorts to weighted least squares (WLS) only [34]. The projection method for TOA localization is firstly proposed by [36]. The work in [37] improves the previous study and derives a closed-form projection method in reaching the CRLB performance. The improved projection method (IPM) doesn’t introduce a redundant variable, so it performs better than the two-step WLS (2WLS) [34]. However, the new projection matrix is not appropriate for some geometries. The research in [38] extends the multidimensional scaling (MDS) method [30] to the case with sensor position uncertainties. The modified MDS method implements better performance than the 2WLS when the sensor position errors are significant. However, the accuracy when using a small wireless sensor network, in which the number of sensors is limited, is poor. The recent research in [39] improves the TDOA localization performance in this situation. To the best of our knowledge, there is no such research on TOA localization in the presence of sensor position errors. The work in [40] illustrates that reducing the bias can improve the performance of the 2WLS method significantly. Therefore, subtracting the estimation bias from the solution in the second stage resulting in a better solution when the noise is relatively large. Since the bias is from the first stage which introduces the redundant variable, eliminating the redundant variable before estimating the bias is more effective, especially when the sensor network consists of fewer sensors.

This paper focuses on improving the TOA localization performance when the sensor network consists of fewer sensors. Converting the sensor position errors and measurement noise together as equivalent measurement noise, the TOA-based source location estimation is formulated as a weighted optimization problem with constraint. The proposed closed-form solution deals with this problem through three stages. The first stage introduces a weighting matrix in terms of the sensor position statistical knowledge to obtain a coarse solution, which is similar to the spherical-interpolation [14], where the constraint is neglected. The covariance of the coarse solution is derived. Then, the solution is improved through the relationship between the source position and redundant variable, using the covariance of stage-1. Stage three recalls the constraint to evaluate the estimate bias in the previous stage. The source position is refined by subtracting the estimated value, resulting in the final estimator. The analytical result shows the proposed method has CRLB performance when the noise and errors are mild. The simulation part verifies the theoretical analysis. It demonstrates that the new method outperforms the compared counterparts when the noise/error level is relatively higher in a small sensor network. Moreover, the mean-squared-error (MSE) and bias are comparable to the existing methods when using more sensors.

The main contributions of this paper are:

• The weighted spherical-interpolation is derived, which solves the source position and redundant variable successively;

• An approximate expression of the theoretical covariance analysis is presented for the weighted spherical-interpolation;

• Eliminating the redundant variable, a refinement for the solution is proposed to improve the MSE and bias further;

• The simulation shows the proposed method performs better than the state-of-the-art methods when using fewer sensors.

The novelties of this paper are:

• Introducing a weighting matrix for spherical-interpolation resulting in the weighted spherical-interpolation;

• Analyzing the covariance in the small noise region;

• Refining the solution to improve the MSE and lower the bias further when using only 4 sensors.

This paper is organized as follows: Section 2 formulates the problem and clarifies the task. The proposed estimator is developed in Section 3. Section 4 presents the CRLB and theoretical performance analysis. Section 5 demonstrates the numerical results to show the advantages of the proposed when using fewer sensors. The conclusion is drawn in Section 6.

The unified notations in this paper is as following: Bold upper and lower case letter/symbol denote the matrix and vector. The vectors in this paper are column vectors. ${\mathbf{x}}^o$ is the true value of $\mathbf{x}$. $\parallel \mathbf{x}\parallel$ represents the Euclidean norm of $\mathbf{x}$. $diag\left\{\mathbf{x}\right\}$ is a diagonal matrix consisted of the elements of $\mathbf{x}$. ${\mathbf{X}}^T$ and ${\mathbf{X}}^{-1}$ are the transpose and inverse of $\mathbf{X}$. $\mathbf{x}(i:j)$ denotes a subvector constructed by the i-th to j-th elements. ⊙ is the Hadamard product. sgn is the the signum operation.

2. Problem Formulation

Consider a scenario with M sensors deployed in a N-dimension space ($N=2$ or 3), whose true positions are ${\mathbf{s}}_i^o$, $i=1,2,\cdots ,M$. Each sensor has the ability to measure the TOA from the unknown source. Denote the true value of the source position is ${\mathbf{u}}^o$. The TOA between the unknown source and sensor i is

(1)${\tau }_i^o=\parallel {\mathbf{u}}^o-{\mathbf{s}}_i^o\parallel /c,$

where c is the propagation speed. Multiplying c on both sides of (1) gives

(2)$r_i^o=\parallel {\mathbf{u}}^o-{\mathbf{s}}_i^o\parallel ,$

where $r_i^o=c{\tau }_i^o$ is called range of arrival (ROA) measurement. In practice, the measured ROA is expressed as the true value with noise,

(3)$r_i=r_i^o+n_i$

where $n_i$ is the measurement noise. Collecting all ROAs together yields the measurement vector

(4)$\mathbf{r}={\mathbf{r}}^o+\mathbf{n},$

where ${\mathbf{r}}^o={[r_1^o,r_2^o,\cdots ,r_M^o]}^T$, $\mathbf{n}$ is modeled zero-mean with Gaussian distribution, and the covariance matrix of $\mathbf{n}$ is ${\mathbf{Q}}_r$.

Noting the true sensor position is unavailable, only the erroneous counterpart is known,

(5)${\mathbf{s}}_i={\mathbf{s}}_i^o+\mathsf{\Delta }{\mathbf{s}}_i.$

The error vector $\mathsf{\Delta }\mathbf{s}={[\mathsf{\Delta }{\mathbf{s}}_1^T,\mathsf{\Delta }{\mathbf{s}}_2^T,\cdots ,\mathsf{\Delta }{\mathbf{s}}_M^T]}^T$ is Gaussian distribution with zero mean and covariance ${\mathbf{Q}}_s$. We further assume that $\mathsf{\Delta }\mathbf{s}$ is independent of the measurement noise $\mathbf{n}$, and $\mathsf{\Delta }\mathbf{s}$ is irrelevant across axis and sensors.

The localization scenario is shown in Figure 1. Our mission is locating the unknown source using the noisy TOA measurements and erroneous sensor positions, where the solution is closed-form. The performance with fewer sensors shall be improved compared with the existing closed-form method.

3. Refined Estimator

The derivation begins with the ROA expression (2). Squaring both sides of (2) and substituting (3) and (5) into the result, we arrive

(6)$\begin{array}{c}-r_i^2+{\mathbf{s}}_i^T{\mathbf{s}}_i-2{\mathbf{s}}_i^T{\mathbf{u}}^o+v^o\hfill \\ =2r_i^o{n}_i+n_i^2+2{({\mathbf{u}}^o-{\mathbf{s}}_i^o)}^T\mathsf{\Delta }{\mathbf{s}}_i-\mathsf{\Delta }{\mathbf{s}}_i^T\mathsf{\Delta }{\mathbf{s}}_i,\hfill \end{array}$

where $v^o={\mathbf{u}}^{oT}{\mathbf{u}}^o$. Stacking (6) across $i=1,2,\cdots ,M$ yields matrix form equation

(7)${\mathbf{h}}_1-\mathbf{A}{\mathit{\psi }}^o=\mathit{ϵ},$

where the unknown is ${\mathit{\psi }}^o={[{\mathbf{u}}^{oT},v^o]}^T$,

(8)${\mathbf{h}}_1=\left[\begin{array}{c}-r_1^2+{\mathbf{s}}_1^T{\mathbf{s}}_1\\ -r_2^2+{\mathbf{s}}_2^T{\mathbf{s}}_2\\ ⋮\\ -r_M^2+{\mathbf{s}}_M^T{\mathbf{s}}_M\end{array}\right],$

(9)$\mathbf{A}={\left[\begin{array}{cccc}2{\mathbf{s}}_1& 2{\mathbf{s}}_2& \cdots & 2{\mathbf{s}}_M\\ 1& 1& \cdots & 1\end{array}\right]}^T,$

(10)$\mathit{ϵ}=\mathbf{B}\mathbf{n}+\mathbf{C}\mathsf{\Delta }\mathbf{s}+{\mathit{ϵ}}^{\prime },$

and

(11)$\mathbf{B}=2diag\left\{{\mathbf{r}}^o\right\},$

(12)$\mathbf{C}=2\left[\begin{array}{ccc}{({\mathbf{u}}^o-{\mathbf{s}}_1^o)}^T& & \\ & \ddots & \\ & & {({\mathbf{u}}^o-{\mathbf{s}}_M^o)}^T\end{array}\right],$

(13)${\mathit{ϵ}}^{\prime }=\mathbf{n}\odot \mathbf{n}-{\left[\begin{array}{c}\mathsf{\Delta }{\mathbf{s}}_1^T\mathsf{\Delta }{\mathbf{s}}_1\\ \mathsf{\Delta }{\mathbf{s}}_2^T\mathsf{\Delta }{\mathbf{s}}_2\\ ⋮\\ \mathsf{\Delta }{\mathbf{s}}_M^T\mathsf{\Delta }{\mathbf{s}}_M\end{array}\right]}^T$

is the second order noise term. Based on the least square criterion, finding ${\mathbf{u}}^o$ is solving the optimization problem

(14)$\underset{\mathit{\psi }}{min}{({\mathbf{h}}_1-\mathbf{A}\mathit{\psi })}^T\mathbf{W}({\mathbf{h}}_1-\mathbf{A}\mathit{\psi })$

(15)$\mathrm{s}.\mathrm{t}.\psi (N+1)=\mathit{\psi }{(1:N)}^T\mathit{\psi }(1:N).$

$\mathbf{W}$ is the weighting matrix approximating by

(16)$\mathbf{W}=E{\left[\mathit{ϵ}{\mathit{ϵ}}^T\right]}^{-1}\simeq {(\mathbf{B}{\mathbf{Q}}_r\mathbf{B}+\mathbf{C}{\mathbf{Q}}_s{\mathbf{C}}^T)}^{-1}.$

Traditional method [34] first minimizes (14) by ignoring the constraint. Then the constraint is utilized in the second stage to improve the accuracy. In this paper, we shall follow a similar step to obtain the initial solution, but the processing is different from the 2-stage WLS (2WLS) [34]. A further stage is needed to refine the bias resulting in the final solution.

3.1. Coarse Solution

Rewriting (7) as

(17)${\mathbf{h}}_1-{\mathbf{G}}_1{\mathbf{u}}^o+\mathbf{z}{v}^o\simeq \mathbf{B}\mathbf{n}+\mathbf{C}\mathsf{\Delta }\mathbf{s}$

where ${\mathbf{G}}_1=\mathbf{A}(:,1:N)$ and $\mathbf{z}={\mathbf{1}}_M$. The second-order noise term ${\mathit{ϵ}}^{\prime }$ is neglected. The existing sphere interpolation (SI) method [14] is suboptimal since it doesn’t take the different error in each row of (7) into account. We introduce a weighting matrix W that is set equal to $E{\left[\mathit{ϵ}{\mathit{ϵ}}^T\right]}^{-1}$ to balance the error in each element of e. Using such a W leads to the weighted sphere interpolation, reducing the covariance in the estimated $\mathit{\psi }$ [41]. The optimization problem in (14) and (15) becomes

(18)$\underset{\mathbf{u},v}{min}{({\mathbf{h}}_1-{\mathbf{G}}_1\mathbf{u}+\mathbf{z}v)}^T\mathbf{W}({\mathbf{h}}_1-{\mathbf{G}}_1\mathbf{u}+\mathbf{z}v)$

(19)$\mathrm{s}.\mathrm{t}.v={\mathbf{u}}^T\mathbf{u}.$

We now turn to solving the constrained optimization problem (18) and (19). If the constraint between $\mathbf{u}$ and v is ignored, the WLS solution of ${\mathbf{u}}^o$ in terms of v is

(20)$\mathbf{u}=\mathbf{H}({\mathbf{h}}_1+\mathbf{z}v),$

where

(21)$\mathbf{H}={\left({\mathbf{G}}_1^T\mathbf{W}{\mathbf{G}}_1\right)}^{-1}{\mathbf{G}}_1^T\mathbf{W}.$

Putting (20) into the right side of (17) yields

(22)${\mathbf{h}}_1-{\mathbf{G}}_1\mathbf{u}+\mathbf{z}v=\mathbf{R}({\mathbf{h}}_1+\mathbf{z}{v}^o),$

where $\mathbf{R}={\mathbf{I}}_M+{\mathbf{G}}_1\mathbf{H}$. Using (22) to the original cost function (18), the optimization problem given by (18) and (19) is equivalent to

(23)$\underset{v}{min}{({\mathbf{h}}_1+\mathbf{z}v)}^T{\mathbf{R}}^T\mathbf{W}\mathbf{R}({\mathbf{h}}_1+\mathbf{z}v),$

subject to (19). Ignoring the constraint (19), finding the suboptimal v is solving the unconstrained optimization problem (23). Taking the derivative of (23) with respect to v and equaling the result to zero, the estimation of $v^o$ is given by the WLS technology

(24)$v=-{\left({\mathbf{z}}^T{\mathbf{R}}^T\mathbf{W}\mathbf{R}\mathbf{z}\right)}^{-1}{\mathbf{z}}^T{\mathbf{R}}^T\mathbf{W}\mathbf{R}{\mathbf{h}}_1.$

To obtain the source location estimate, the value of v calculated from (24) is substituted to (20).

3.2. Error Reduction

We next explore the constraint to reduce the error of the coarse $\mathbf{u}$ and v given by (20) and (24). Before formulating the problem in this stage, the covariance expression of $\mathbf{u}$ and v shall be detailed. To the best knowledge of the authors, there is no such analysis for weighted spherical-interpolation.

Subtracting both sides of (20) by the true value ${\mathbf{u}}^o$ yields

(25)$\mathsf{\Delta }\mathbf{u}=\mathbf{H}({\mathbf{h}}_1+\mathbf{z}v-{\mathbf{G}}_1{\mathbf{u}}^o).$

Using the measurement Equation (17), ${\mathbf{G}}_1{\mathbf{u}}^o$ can be expressed by

(26)${\mathbf{G}}_1{\mathbf{u}}^o={\mathbf{h}}_1+\mathbf{z}{v}^o-\mathbf{B}\mathbf{n}-\mathbf{C}\mathsf{\Delta }\mathbf{s}.$

Substituting (26) to (25) and using $\mathsf{\Delta }v=v-v^o$ yield

(27)$\mathsf{\Delta }\mathbf{u}=\mathbf{H}(\mathbf{z}\mathsf{\Delta }v+\mathbf{B}\mathbf{n}+\mathbf{C}\mathsf{\Delta }\mathbf{s}).$

$\mathsf{\Delta }v$ is obtained by subtracting both sides of v (24) by the true value

(28)$\begin{array}{cc}\mathsf{\Delta }v& =-{\left({\mathbf{z}}^T{\mathbf{R}}^T\mathbf{W}\mathbf{R}\mathbf{z}\right)}^{-1}{\mathbf{z}}^T{\mathbf{R}}^T\mathbf{W}\mathbf{R}({\mathbf{h}}_1+\mathbf{z}{v}^o)\hfill \\ & ={\mathit{\rho }}^T(\mathbf{B}\mathbf{n}+\mathbf{C}\mathsf{\Delta }\mathbf{s}),\hfill \end{array}$

where

(29)$\mathit{\rho }=-{\left({\mathbf{z}}^T{\mathbf{R}}^T\mathbf{W}\mathbf{R}\mathbf{z}\right)}^{-1}{\mathbf{R}}^T\mathbf{W}\mathbf{R}\mathbf{z}.$

Inserting (28) into (25) yields the expression of $\mathsf{\Delta }\mathbf{u}$,

(30)$\mathsf{\Delta }\mathbf{u}=\mathbf{P}(\mathbf{B}\mathbf{n}+\mathbf{C}\mathsf{\Delta }\mathbf{s})$

where

(31)$\mathbf{P}=\mathbf{H}({\mathbf{I}}_M+\mathbf{z}{\mathit{\rho }}^T).$

Denoting $\mathsf{\Delta }\mathit{\psi }={[\mathsf{\Delta }{\mathbf{u}}^T,\mathsf{\Delta }v]}^T$, the covariance is

(32)$\begin{array}{cc}cov\left(\mathit{\psi }\right)& =E\left[\mathsf{\Delta }\mathit{\psi }\mathsf{\Delta }{\mathit{\psi }}^T\right]\hfill \\ & =E\left\{\left[\begin{array}{cc}\mathsf{\Delta }\mathbf{u}\mathsf{\Delta }{\mathbf{u}}^T& \mathsf{\Delta }\mathbf{u}\mathsf{\Delta }v\\ \mathsf{\Delta }v\mathsf{\Delta }{\mathbf{u}}^T& \mathsf{\Delta }{v}^2\end{array}\right]\right\}\hfill \\ & =\left[\begin{array}{cc}\mathbf{P}{\mathbf{W}}^{-1}{\mathbf{P}}^T& \mathbf{P}{\mathbf{W}}^{-1}\mathit{\rho }\\ {\mathit{\rho }}^T{\mathbf{W}}^{-1}{\mathbf{P}}^T& {\mathit{\rho }}^T{\mathbf{W}}^{-1}\mathit{\rho }\end{array}\right].\hfill \end{array}$

To reduce the errors in the coarse solution, a new formula shall be established in terms of the constraint (19). Recall the relationship between the source position ${\mathbf{u}}^o$ and the residual variable $v^o$,

(33)$v^o={\mathbf{u}}^{oT}{\mathbf{u}}^o.$

Commonly, it is reasonable to consider the errors in the estimated unknown parameters, e.g., $\mathsf{\Delta }\mathbf{u}$ and $\mathsf{\Delta }v$, as additional [33,34,40,42]. Thus, we can express the redundant variable v as

(34)$v=v^o+\mathsf{\Delta }v.$

Substituting (33) to (34) results

(35)$v-{\mathbf{u}}^{oT}{\mathbf{u}}^o=\mathsf{\Delta }v.$

Squaring both sides of

(36)$\mathbf{u}={\mathbf{u}}^o+\mathsf{\Delta }\mathbf{u}$

elementwise results in

(37)$\mathbf{u}\odot \mathbf{u}-{\mathbf{u}}^o\odot {\mathbf{u}}^o=2{\mathbf{u}}^o\odot \mathsf{\Delta }\mathbf{u},$

where the second-order term of $\mathsf{\Delta }\mathbf{u}$ is dropped. Together with the equation obtained in (35), we have

(38)${\mathbf{h}}_1^{\prime }-{\mathbf{G}}_1^{\prime }{\mathit{\psi }}^{\prime o}={\mathbf{B}}^{\prime o}\mathsf{\Delta }\mathit{\psi },$

where

(39)${\mathit{\psi }}^{\prime o}={\mathbf{u}}^o\odot {\mathbf{u}}^o,$

(40)${\mathbf{h}}_1^{\prime }=\left[\begin{array}{c}\mathbf{u}\odot \mathbf{u}\\ v\end{array}\right],$

(41)${\mathbf{G}}_1^{\prime }=\left[\begin{array}{c}{\mathbf{I}}_N\\ {\mathbf{1}}_N^T\end{array}\right],$

(42)${\mathbf{B}}^{\prime o}=\left[\begin{array}{cc}2diag\left\{{\mathbf{u}}^{oT}\right\}& \\ & 1\end{array}\right].$

Thus, the estimate of ${\mathit{\psi }}^{\prime o}$ is the one that satisfies

(43)$\underset{{\mathit{\psi }}^{\prime }}{min}{({\mathbf{h}}_1^{\prime }-{\mathit{\psi }}^{\prime })}^T{\mathbf{W}}^{\prime }({\mathbf{h}}_1^{\prime }-{\mathbf{G}}_1^{\prime }{\mathit{\psi }}^{\prime }),$

which is given by

(44)${\mathit{\psi }}^{\prime }={\left({\mathbf{G}}_1^{\prime T}{\mathbf{W}}^{\prime }{\mathbf{G}}_1^{\prime }\right)}^{-1}{\mathbf{G}}_1^{\prime T}{\mathbf{W}}^{\prime }{\mathbf{h}}_1^{\prime },$

where the weighting matrix ${\mathbf{W}}^{\prime }$ is approximated as

(45)${\mathbf{W}}^{\prime }={({\mathbf{B}}^{\prime }cov\left(\psi \right){\mathbf{B}}^{\prime })}^{-1},$

where ${\mathbf{B}}^{\prime o}$ is replaced by ${\mathbf{B}}^{\prime }$ using $\mathbf{u}$ instead. ${\mathbf{W}}^{\prime }$ is related to $\mathbf{W}$ by (32).

The estimation of ${\mathbf{u}}^o$ is the square root of ${\mathit{\psi }}^{\prime }$ in term of (39). However, there is sign ambiguity since using only ${\mathit{\psi }}^{\prime }$ can’t determine the sign of $\mathbf{u}$ elementwise. One way to solve this sign ambiguity is to use the sign given by the first stage [10]. Thus, the final error-reduced estimate is

(46)${\mathbf{u}}^{\prime }=\mathrm{sgn}\left(\mathbf{u}\right)\odot \sqrt{{\mathit{\psi }}^{\prime }}$

where $\mathbf{u}$ is the result given in the first stage (20), and

(47)$v^{\prime }={\mathbf{u}}^{\prime T}{\mathbf{u}}^{\prime }$

in terms of (33).

3.3. Refinement

The last stage will refine the solution above. Express the ${\mathbf{u}}^{\prime }$ as

(48)${\mathbf{u}}^{\prime }={\mathbf{u}}^o+\mathsf{\Delta }{\mathbf{u}}^{\prime }$

where $\mathsf{\Delta }{\mathbf{u}}^{\prime }$ is the estimation error. So the expression of $v^o$ in terms of ${\mathbf{u}}^{\prime }$ and $\mathsf{\Delta }{\mathbf{u}}^{\prime }$ is

(49)$v^o={({\mathbf{u}}^{\prime }-\mathsf{\Delta }{\mathbf{u}}^{\prime })}^T({\mathbf{u}}^{\prime }-\mathsf{\Delta }{\mathbf{u}}^{\prime })\approx {\mathbf{u}}^{\prime T}{\mathbf{u}}^{\prime }-2{\mathbf{u}}^{\prime T}\mathsf{\Delta }{\mathbf{u}}^{\prime }.$

Substituting (48) and (49) to (17) and dropping the second order noise term yield

(50)${\mathbf{h}}_2-{\mathbf{G}}_2\mathsf{\Delta }{\mathbf{u}}^{\prime }\simeq \mathbf{B}\mathbf{n}+\mathbf{C}\mathsf{\Delta }\mathbf{s},$

where

(51)${\mathbf{h}}_2={\mathbf{h}}_1-{\mathbf{G}}_1{\mathbf{u}}^{\prime }+\mathbf{z}{v}^{\prime },$

(52)${\mathbf{G}}_2=2{\left[\begin{array}{c}{({\mathbf{u}}^{\prime }-{\mathbf{s}}_1)}^T\\ {({\mathbf{u}}^{\prime }-{\mathbf{s}}_2)}^T\\ ⋮\\ {({\mathbf{u}}^{\prime }-{\mathbf{s}}_M)}^T\end{array}\right]}^T.$

The estimation of error $\mathsf{\Delta }\mathbf{u}$ is

(53)$\mathsf{\Delta }{\widehat{\mathbf{u}}}^{\prime }={\left({\mathbf{G}}_2^T\mathbf{W}{\mathbf{G}}_2\right)}^{-1}{\mathbf{G}}_2^T\mathbf{W}{\mathbf{h}}_2.$

The final solution is

(54)$\widehat{\mathbf{u}}={\mathbf{u}}^{\prime }-\mathsf{\Delta }{\widehat{\mathbf{u}}}^{\prime }.$

Remark 1.

The weighting matrix $\mathbf{W}$ is determined by the true value of the source position. A remedy is setting $\mathbf{B}$ as identity and $\mathbf{C}={\mathbf{I}}_M\otimes {\mathbf{I}}_N/N$ for a coarse solution. Then, the weighting matrix can be updated through the coarse solution and a better initial solution result. In the refinement stage, the updated $\mathbf{W}$ is sufficient for an accurate final solution.

Remark 2.

The relationship (33) has been used in Section 3.2, but it is insufficient to obtain the best estimate. Normally, the constrained weighted least squares (CWLS) is better than the two-stage WLS [43,44,45,46], which means dealing with the quadratic cost function and the constraint subsequently is suboptimal. Although the 2WLS approaches the CRLB in small noise region, its performance degradation is distinct as the noise level increases. The work in [40] illustrates that the large bias is the reason for performance degradation. Reducing the bias of the first stage can improve performance significantly. Therefore, subtracting the estimation bias from the solution in the second stage can improve the performance when the noise is relatively large. Since the bias is from the first stage which introduces the redundant variable, eliminating the redundant variable before estimating the bias is more effective when the sensor network consists of fewer sensors.

4. Analysis

4.1. CRLB

The CRLB is the most widely applied technique for bounding the ability to estimate the interest parameters from the given data. Although the CRLB of TOA localization when the sensor positions are uncertain has been developed in [34], we will present the CRLB from the other perspective for the theoretical performance comparison.

Equation (10) can be rewritten as

(55)$\mathit{ϵ}=\mathbf{B}\tilde{\mathbf{n}},$

where

(56)$\tilde{\mathbf{n}}=\mathbf{n}+{\mathbf{B}}^{-1}\mathbf{C}\mathsf{\Delta }\mathbf{s}.$

In (56), the sensor errors are linearly transformed as an additional term of measurement noise. Thus, $\tilde{\mathbf{n}}$ is called the equivalent noise. This transformation combines the noise and sensor position errors, which simplifies the algorithm derivation and analysis expression.

In terms of the equivalent noise, the corresponding covariance matrix is

(57)$\mathbf{Q}=E\left[\tilde{\mathbf{n}}{\tilde{\mathbf{n}}}^T\right]={\mathbf{Q}}_r+{\mathbf{B}}^{-1}\mathbf{C}{\mathbf{Q}}_s{\mathbf{C}}^T{\mathbf{B}}^{-1}.$

Then, the CRLB is given by [41]

(58)$\mathbf{CRLB}\left({\mathbf{u}}^o\right)={\left(\frac{\partial {\mathbf{r}}^{oT}}{\partial {\mathbf{u}}^o}\mathbf{Q}\frac{\partial {\mathbf{r}}^o}{\partial {\mathbf{u}}^{oT}}\right)}^{-1}$

where

(59)$\frac{\partial {\mathbf{r}}^o}{\partial {\mathbf{u}}^{oT}}=\left[\begin{array}{c}{({\mathbf{u}}^o-{\mathbf{s}}_1^o)}^T/r_1^o\\ {({\mathbf{u}}^o-{\mathbf{s}}_1^o)}^T/r_1^o\\ ⋮\\ {({\mathbf{u}}^o-{\mathbf{s}}_M^o)}^T/r_M^o\end{array}\right].$

4.2. Covariance

The estimator in the first and the second stage can be expressed as

(60)$\begin{array}{cc}\hfill {\mathbf{u}}^{\prime }& ={\mathbf{u}}^o+\mathsf{\Delta }{\mathbf{u}}^{\prime },\hfill \\ \hfill \mathsf{\Delta }{\widehat{\mathbf{u}}}^{\prime }& =\mathsf{\Delta }{\mathbf{u}}^{\prime }+\delta {\mathbf{u}}^{\prime }.\hfill \end{array}$

Substituting (60) to (54) yields

(61)$\widehat{\mathbf{u}}={\mathbf{u}}^o-\delta {\mathbf{u}}^{\prime }.$

Therefore, the final estimation error is $-\delta \mathbf{u}$,

(62)$\begin{array}{cc}-\delta {\mathbf{u}}^{\prime }& =\widehat{\mathbf{u}}-{\mathbf{u}}^o\hfill \\ & =-{\left({\mathbf{G}}_2^T\mathbf{W}{\mathbf{G}}_2\right)}^{-1}{\mathbf{G}}_2^T\mathbf{W}(\mathbf{B}\mathbf{n}+\mathbf{C}\mathsf{\Delta }\mathbf{s}).\hfill \end{array}$

After neglecting the noise/error terms higher than second order, the covariance of $\widehat{\mathbf{u}}$ is approximately

(63)$cov\left({\widehat{\mathbf{u}}}^{\prime }\right)=E{\left[\delta {\mathbf{u}}^{\prime }\delta {\mathbf{u}}^{\prime T}\right]}^{-1}\simeq {\left({\mathbf{G}}_2^T\mathbf{W}{\mathbf{G}}_2\right)}^{-1}.$

4.3. Comparison

Under the small noise conditions

(64)$\begin{array}{c}1)n_i/r_i^o\approx 0,\hfill \\ 2)\parallel \mathsf{\Delta }{\mathbf{s}}_i\parallel /r_i^o\approx 0,\hfill \end{array}$

Substituting (16) and (57) to (63), the inverse of the theoretical covariance of $\widehat{\mathbf{u}}$ is

(65)$cov{\left(\widehat{\mathbf{u}}\right)}^{-1}\approx {\mathbf{G}}_2^{oT}{\mathbf{B}}^{-1}\mathbf{Q}{\mathbf{B}}^{-1}{\mathbf{G}}_2^o.$

For simplicity, we denote

(66)${\mathbf{G}}_3^o={\mathbf{B}}^{-1}{\mathbf{G}}_2^o.$

Substituting (11) and (52) into (66), after some algebraic manipulations, and comparing the result with (59) yield

(67)${\mathbf{G}}_3^o=\frac{\partial {\mathbf{r}}^o}{\partial {\mathbf{u}}^{oT}}.$

Therefore, we conclude that (63) equals to the CRLB if the noise is mild.

5. Simulations

Simulations run with $M=6$ sensors. It is worth to note that the proposed method is applicable to both two-dimensional and three-dimensional cases. Only the results for the more complicated case of 3-D are illustrated in this Section. The certain sensor positions are listed in Table 1. We use a near source and a far source for the experiments: ${\mathbf{u}}_{near}^o={[400,350,550]}^T$ m and ${\mathbf{u}}_{far}^o={[2000,1750,2250]}^T$ m. The covariance matrix of TOAs and sensor positions are ${\mathbf{Q}}_r={\sigma }_r^2{\mathbf{I}}_M$ and ${\mathbf{Q}}_s={\sigma }_s^{2}diag\left\{\mathit{\rho }\right\}\otimes {\mathbf{I}}_N$, where $\mathit{\rho }$ is a scaling vector for sensor position errors. The proposed method is compared with three closed-form solutions, 2WLS [34], improved projection method (IPM) [37] and multidimensional scaling (MDS) based method [30], and the maximum likelihood estimator (MLE) implemented through Gauss-Newton iteration. In addition, the CRLB is introduced as a benchmark. The number of ensemble runs is $K=1000$, unless stated otherwise. The performance is evaluated by MSE and bias,

(68)$MSE(\mathbf{u})=\frac{1}{K}\sum _{k=1}^K{({\mathbf{u}}_k-{\mathbf{u}}^o)}^T({\mathbf{u}}_k-{\mathbf{u}}^o),$

(69)$bias(\mathbf{u})=\frac{1}{K}\sum _{k=1}^K\parallel {\mathbf{u}}_k-{\mathbf{u}}^o\parallel ,$

where ${\mathbf{u}}_k$ is the estimates at ensemble run k.

5.1. Fewer Sensors

The motivation of this research is reducing the MSE and bias when using fewer sensors. Figure 2, Figure 3, Figure 4 and Figure 5 illustrate the MSE and bias as the TOA measurement noise increases, where the sensor network consist of 4 sensors, ${\mathbf{s}}_1$, ${\mathbf{s}}_2$, ${\mathbf{s}}_3$ and ${\mathbf{s}}_4$. The sensor position error level is fixed at ${\sigma }_s^2={10}^{-5}$ m². The scaling vector is $\mathit{\rho }={[10,2,10,40]}^T$. For the near source, the 2WLS, IPM and the proposed algorithm have comparable MSE and bias if ${\sigma }_r^2⩽{10}^{-2}$ m², while the MDS is more robust than 2WLS, IPM and even the proposed in relatively larger error region, although the sensor position errors are not taken into account. When the noise increases, the 2WLS, IPM and MDS deviate the CRLB after the noise power higher than ${10}^{-2}$ m², while the proposed method diverges over 15 dB later than the 2WLS and IPM. If the noise increases further, the proposed method keeps the MLE-level MSE and bias even when the noise is significantly large. The MLE deviates the MLE after ${\sigma }_r^2>1$ m². The result is similar for the far source, the proposed method deviates the CRLB when ${\sigma }_r^2>{10}^{-1.5}$ m², while the performance of the 2WLS, IPM and MDS become worse at ${\sigma }_r^2={10}^{-3}$ m². The proposed method outperforms the 2WLS and IPM, even after deviating the CRLB. And the proposed algorithm lowers the bias, especially when the ${\sigma }_s^2>{10}^{-2}$ m² for near target and ${\sigma }_r^2>{10}^{-3}$ m² for far source. The MLE is the best but delays the threshold only 5 dB.

The performance behavior as the position error level of sensors varies is demonstrated in Figure 6, Figure 7, Figure 8 and Figure 9. The sensors’ deployment is the same as the experiments above. The TOA measurement noise power is set as ${\sigma }_r^2={10}^{-3}$ m². If the source is near, Figure 6 indicates that the 2WLS, IPM and MDS fail to reach the CRLB after ${\sigma }_s^2>{10}^{-4}$ m², ${10}^{-4}$ m² and ${10}^{-4.5}$ m² respectively, while the proposed method delays the divergence to ${\sigma }_s^2>{10}^{-2}$ m². And the MLE delays the threshold 5 dB further. In Figure 7, five methods asymptotically have the same level bias in small error region. The proposed method can yield a much smaller bias about 15 dB around the error level at ${\sigma }_s^2>{10}^{-2}$ m². It lowers the bias over 10 dB than the compared three closed-form methods at most, which is comparable to the MLE. Moreover, Figure 7 and Figure 9 imply the reason that 2WLS, IPM and MDS deviate the CRLB earlier. The raising bias deteriorates the MSE performance. The results when locating a far source are similar to before. The proposed achieves the CRLB accuracy until ${\sigma }_s^2={10}^{-3}$ m² and have bias less than that of the three existing closed-form methods since ${\sigma }_s^2>{10}^{-5}$ m².

To show the effectiveness of the proposed work more clearly, the performance loss is presented with the percentage with respect to the CRLB. The results for near source and far source are listed in Table 2 and Table 3. The proposed method provides comparable MSE performance to the MLE before the divergence. The MSE of the proposed and MLE is very close to the CRLB where the difference is less than 2%. Noting the ensemble runs is limited, and the bias has been subtracted from the generated noise, the MSE may be slightly lower than the CRLB. 2WLS and IPM have larger performance loss with respect to the CRLB. The MDS is more robust than 2WLS, IPM and even the proposed in relatively larger error region, although the sensor position errors are not taken into account.

5.2. More Sensors

The expectation of the proposed solution is improving the performance when using fewer sensors, while the accuracy under a normal sensor network with sufficient sensors should not be worse than the existing methods. The simulations above illustrate the advantages of the proposed closed-form solution on both MSE and bias. In this subsection, we’d like to show the new approach has comparable MSE and bias to the compared closed-form solutions. The scaling vector is $\mathit{\rho }={[10,2,10,40,20,3]}^T$. To limit the number of illustrations, we only present the results for single near source in Figure 10, Figure 11, Figure 12 and Figure 13. They display similar performance of all methods as the measurement noise power or sensor error level increases, albeit the proposed is slightly better than the 2WLS, IPM and MDS when ${\sigma }_r^2>{10}^2$ m² or ${\sigma }_s^2>10$ m². The simulation results verify that the proposed method improves the performance when using 4 sensors but does as good as the existing methods when using more sensors, rather than designed for 4 sensors-case only.

5.3. Computation Time

The processing time of the proposed method and the compared counterparts is collected by Matlab 2018b on a typical PC with AMD R5 3500X CPU. The running time of 20000 independent tests is listed in Table 4. The proposed method is less computation-complexity than MLE, IPM and MDS and is comparable with the 2WLS method.

5.4. Summary of Simulations

We shall summarize the simulation results. The proposed method performs better than the existing methods because it gives a more accurate coarse location than IPM and refines the bias that 2WLS and MDS do not. The coarse solution of IPM can’t reach the CRLB in some cases, e.g., the configuration used in this paper. After refinement, the performance is still worse than the proposed. However, the coarse solution of the proposed is able to attain the CRLB when the noise is quite small, as well as the final solution of 2WLS and MDS. The refining step reduces bias further, as shown in Figure 3, Figure 5, Figure 7 and Figure 9. That is the reason why the proposed outperforms the 2WLS and MDS if the noise increases from 10⁻² m² to 1 m² for a near source or from 10⁻³ m² to 0.1 m² for a far source.

The proposed method performs worse than MLE and MDS when the noise power is relatively high. MLE is implemented by the Gauss-Newton iteration, which is the best for parameter estimation. The MDS has been verified to be robust for the source localization. It has better MSE than the proposed method in larger noise region, but it is inferior to the proposed both MSE and bias in other cases.

Moreover, the proposed method is less computation-consumption but the performance is better than the existing in relatively large noise and error conditions.

6. Conclusions

A refined TOA localization solution is proposed in this paper, with better MSE and bias than the existing methods when the sensors composing the network are few. This new solution is closed-form, whose three stages take the WLS only. The MSE is analyzed theoretically, which is proven mathematically attaining the CRLB if the measurement noise and sensor errors are small. The MSE and bias are examined in the simulations, verifying the analytical result and illustrating the superiority of the proposed method.

Locating a source in a N-dimensional (N-D) space using TOA requires N sensors at least, which means the minimum sensor numbers for a N-D space target positioning is N. However, the proposed method, as well as the 2WLS and IPM, are not applicable when there are only N sensors since the redundant parameter is introduced for pseudo-linearization. The TOA localization using the minimal sensor network may be an interesting subject for further study.

Author Contributions

Conceptualization, Methodology and Project administration, Y.G. Formal analysis and Investigation, X.C.; Validation and Writing—original draft, Y.S.; Writing—review & editing, X.C. and Y.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The authors would like to thank the anonymous reviewers for providing valuable comments and suggestions that have helped improving the quality of the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Localization scenario. Open circles are the unknown true sensor positions and open triangles represent the available imprecise sensor positions.

Figure 2. MSE performance as noise power increasing (near source, 4 sensors).

Figure 3. Bias performance as noise power increasing (near source, 4 sensors).

Figure 4. MSE performance as noise power increasing (far source, 4 sensors).

Figure 5. Bias performance as noise power increasing (far source, 4 sensors).

Figure 6. MSE performance as error level increasing (near source, 4 sensors).

Figure 7. Bias performance as error level increasing (near source, 4 sensors).

Figure 8. MSE performance as error level increasing (far source, 4 sensors).

Figure 9. Bias performance as error level increasing (far source, 4 sensors).

Figure 10. MSE performance as noise power increasing (near source, 6 sensors).

Figure 11. Bias performance as noise power increasing (near source, 6 sensors).

Figure 12. MSE performance as error level increasing (near source, 6 sensors).

Figure 13. Bias performance as error level increasing (near source, 6 sensors).

Table 1

Certain sensor positions.

Sensor	${\mathit{s}}_1$	${\mathit{s}}_2$	${\mathit{s}}_3$	${\mathit{s}}_4$	${\mathit{s}}_5$	${\mathit{s}}_6$
x (m)	−100	200	400	350	300	300
y (m)	100	−300	150	200	500	100
z (m)	−100	−200	100	100	200	150

Table 2

MSE increasing percentage (near source).

Noise Power (dB)	−60	−50	−40	−30	−20	−10	0	10	20
Proposed (%)	0.42	−1.42	−0.57	−1.62	0.34	1.10	3.04 × 103	4.55 × 104	9.61 × 104
2WLS (%)	3.54	1.74	1.38	3.37	34.15	367.12	2.96 × 104	7.55 × 104	1.52 × 104
IPM (%)	5.59	4.36	1.31	5.48	50.56	538.84	5.78 × 104	2.23 × 105	2.62 × 104
MLE (%)	0.40	−1.22	−0.94	−1.38	0.34	1.05	1.27	3.60 × 104	8.39 × 103
MDS (%)	54.86	44.10	12.83	5.44	49.05	508.94	8.24 × 103	2.89 × 104	6.79 × 103

Table 3

MSE increasing percentage (far source).

Noise Power (dB)	−60	−50	−40	−30	−20	−10	0
Proposed (%)	0.92	−1.84	−1.27	−1.06	0.93	1.64 × 103	3.12 × 103
2WLS (%)	14.34	10.58	6.43	19.37	178.63	1.30 × 104	6.73 × 104
IPM (%)	18.34	14.94	8.85	27.36	237.50	1.31 × 104	7.90 × 104
MLE (%)	0.86	−0.88	−1.32	−1.68	0.88	1.98	2.09 × 104
MDS (%)	62.80	53.94	25.53	30.20	225.49	3.20 × 103	1.59 × 104

Table 4

Computation time.

Method	Proposed	2WLS	MLE	IPM	MDS
Time (s)	29.90	21.69	142.07	72.80	112.15
Rel. Time	1	0.73	4.75	2.43	3.75