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

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.


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: 1.
The weighted spherical-interpolation is derived, which solves the source position and redundant variable successively; 2.
An approximate expression of the theoretical covariance analysis is presented for the weighted spherical-interpolation; 3.
Eliminating the redundant variable, a refinement for the solution is proposed to improve the MSE and bias further; 4.
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: 1.
Introducing a weighting matrix for spherical-interpolation resulting in the weighted spherical-interpolation; 2.
Analyzing the covariance in the small noise region; 3.
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. x o is the true value of x.
x represents the Euclidean norm of x. diag{x} is a diagonal matrix consisted of the elements of x. X T and X −1 are the transpose and inverse of X. 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.

Problem Formulation
Consider a scenario with M sensors deployed in a N-dimension space (N = 2 or 3), whose true positions are s o i , i = 1, 2, · · · , M. Each sensor has the ability to measure the TOA from the unknown source. Denote the true value of the source position is u o . The TOA between the unknown source and sensor i is where c is the propagation speed. Multiplying c on both sides of (1) gives where r o i = cτ o i is called range of arrival (ROA) measurement. In practice, the measured ROA is expressed as the true value with noise, where n i is the measurement noise. Collecting all ROAs together yields the measurement vector where n is modeled zero-mean with Gaussian distribution, and the covariance matrix of n is Q r .
Noting the true sensor position is unavailable, only the erroneous counterpart is known, The error vector ∆s = [∆s T 1 , ∆s T 2 , · · · , ∆s T M ] T is Gaussian distribution with zero mean and covariance Q s . We further assume that ∆s is independent of the measurement noise n, and ∆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.

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 where v o = u oT u o . Stacking (6) across i = 1, 2, · · · , M yields matrix form equation where the unknown is and is the second order noise term. Based on the least square criterion, finding u o is solving the optimization problem W is the weighting matrix approximating by 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.

Coarse Solution
Rewriting (7) as where G 1 = A(:, 1 : N) and z = 1 M . The second-order noise term 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[ T ] −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 ψ [41]. The optimization problem in (14) and (15) becomes We now turn to solving the constrained optimization problem (18) and (19). If the constraint between u and v is ignored, the WLS solution of u o in terms of v is where Putting (20) into the right side of (17) yields where R = I M + G 1 H. Using (22) to the original cost function (18), the optimization problem given by (18) and (19) is equivalent to 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 To obtain the source location estimate, the value of v calculated from (24) is substituted to (20).

Error Reduction
We next explore the constraint to reduce the error of the coarse u and v given by (20) and (24). Before formulating the problem in this stage, the covariance expression of 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 u o yields Using the measurement Equation (17), G 1 u o can be expressed by Substituting (26) to (25) and ∆v is obtained by subtracting both sides of v (24) by the true value where Inserting (28) where Denoting ∆ψ = [∆u T , ∆v] T , the covariance is 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 u o and the residual variable v o , Commonly, it is reasonable to consider the errors in the estimated unknown parameters, e.g., ∆u and ∆v, as additional [33,34,40,42]. Thus, we can express the redundant variable v as Substituting (33) to (34) results Squaring both sides of where the second-order term of ∆u is dropped. Together with the equation obtained in (35), we have where Thus, the estimate of ψ o is the one that satisfies which is given by where the weighting matrix W is approximated as where B o is replaced by B using u instead. W is related to W by (32). The estimation of u o is the square root of ψ in term of (39). However, there is sign ambiguity since using only ψ can't determine the sign of 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 where u is the result given in the first stage (20), and in terms of (33).

Refinement
The last stage will refine the solution above. Express the u as where ∆u is the estimation error. So the expression of v o in terms of u and ∆u is Substituting (48) and (49) to (17) and dropping the second order noise term yield where The estimation of error ∆u is The final solution isû = u − ∆û .
Remark 1. The weighting matrix W is determined by the true value of the source position. A remedy is setting B as identity and C = I M ⊗ 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 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.

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 where n = n + B −1 C∆s .
In (56), the sensor errors are linearly transformed as an additional term of measurement noise. Thus, 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 Then, the CRLB is given by [41]

Covariance
The estimator in the first and the second stage can be expressed as Substituting (60) to (54) yieldsû Therefore, the final estimation error is −δu, After neglecting the noise/error terms higher than second order, the covariance ofû is approximately

Comparison
Under the small noise conditions Substituting (16) and (57) to (63), the inverse of the theoretical covariance ofû is For simplicity, we denote Substituting (11) and (52) into (66), after some algebraic manipulations, and comparing the result with (59) yield Therefore, we conclude that (63) equals to the CRLB if the noise is mild.

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 [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, where u k is the estimates at ensemble run k.

Fewer Sensors
The motivation of this research is reducing the MSE and bias when using fewer sensors. Figures 2-5 illustrate the MSE and bias as the TOA measurement noise increases, where the sensor network consist of 4 sensors, s 1 , s 2 , s 3 and s 4 . The sensor position error level is fixed at σ 2 s = 10 −5 m 2 . The scaling vector is ρ = [10, 2, 10, 40] T . For the near source, the 2WLS, IPM and the proposed algorithm have comparable MSE and bias if σ 2 r 10 −2 m 2 , 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 2 , 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 σ 2 r > 1 m 2 . The result is similar for the far source, the proposed method deviates the CRLB when σ 2 r > 10 −1.5 m 2 , while the performance of the 2WLS, IPM and MDS become worse at σ 2 r = 10 −3 m 2 . The proposed method outperforms the 2WLS and IPM, even after deviating the CRLB. And the proposed algorithm lowers the bias, especially when the σ 2 s > 10 −2 m 2 for near target and σ 2 r > 10 −3 m 2 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 Figures 6-9. The sensors' deployment is the same as the experiments above. The TOA measurement noise power is set as σ 2 r = 10 −3 m 2 . If the source is near, Figure 6 indicates that the 2WLS, IPM and MDS fail to reach the CRLB after σ 2 s > 10 −4 m 2 , 10 −4 m 2 and 10 −4.5 m 2 respectively, while the proposed method delays the divergence to σ 2 s > 10 −2 m 2 . 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 σ 2 s > 10 −2 m 2 . It lowers the bias over 10 dB than the compared three closed-form methods at most, which is comparable to the MLE. Moreover, Figures 7 and 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 σ 2 s = 10 −3 m 2 and have bias less than that of the three existing closed-form methods since σ 2 s > 10 −5 m 2 . 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 Tables 2 and 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.

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 ρ = [10, 2, 10, 40, 20, 3] T . To limit the number of illustrations, we only present the results for single near source in Figures 10-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 σ 2 r > 10 2 m 2 or σ 2 s > 10 m 2 . 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.

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.

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 Figures 3, 5, 7 and 9. That is the reason why the proposed outperforms the 2WLS and MDS if the noise increases from 10 −2 m 2 to 1 m 2 for a near source or from 10 −3 m 2 to 0.1 m 2 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.

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.