Efficient Bi-Iterative Method for Source Position and Propagation Speed Estimation Using TDOA Measurements

This paper develops an efficient bi-iterative source location and propagation speed estimation method utilizing time difference of arrival (TDOA) measurements. The source location and propagation speed estimation is a nonlinear problem due to the nonlinearity in the TDOA measurement equations. The newly developed bi-iterative method computes the source location and propagation speed alternately. The asymptotic convergence of the new bi-iterative method is theoretically analyzed. First-order perturbation analysis is applied to the newly developed solution to derive its bias and variance. The first-order analytical results show that the proposed method provides approximately unbiased source position and propagation speed estimates for low noise levels and the accuracy of these estimates approaches the Cramer-Rao lower bound (CRLB). The extension of the new bi-iterative method to the more general situation where the sensor locations are subject to random errors is also presented. Simulation studies are given to show the good performance of the proposed method.


Introduction
Passive source localization using time difference of arrival (TDOA) measurements has received considerable attention and has been widely applied in target tracking [1,2], navigation [3], sensor networks [4,5], and wireless communications [6,7]. In the past decades, a number of efficient algorithms such as those in [5][6][7][8][9] were presented for TDOA-based source localization. Nevertheless, all these works assume that the signal propagation speed is known a priori so that the obtained TDOA measurements can be converted into range differences for source positioning. In practical applications such as seismic exploration [10], tangible interface for human-computer interaction [11], and underwater acoustic [12], the propagation speed is unknown and depends strongly on the propagation medium. In this case, the unknown propagation speed needs to be estimated jointly with the source location. For this problem, Mahajan and Walworth [13] proposed an unconstrained least-squares (LS) method, in which the nonlinear TDOA measurement equations are converted into pseudolinear ones by introducing two auxiliary variables. The source location and two auxiliary variables are then jointly estimated in LS sense. Reed et al. [14] selected the LS solution as the starting value and developed a four-step method which alternately estimates the source location and propagation speed. Recently, Zheng et al. [15] proposed a three-stage approach to simultaneously compute the source location and propagation speed. Very interestingly, the accuracy of the source location and propagation speed estimates approximates the CRLB for sufficiently small noise conditions. A disadvantage of this three-stage method is that it produces two results, where only one is the true solution. Furthermore, it may generate complex solutions when finding the square root in its last stage. Annibale and Rabenstein [16] investigated the influence of a wrongly presumed propagation speed because of temperature variations on the positioning accuracy of two well-known location methods. Oyzerman and Amar [17] extended the spherical intersection (SX) technique in [5] to the joint source position and propagation speed estimation. Annibale and Rabenstein [18] derived two closed-form methods for estimation of the propagation speed using time of arrival (TOA) and TDOA measurements. However, these works [13][14][15][16][17][18] mentioned above are based on the assumption that the sensor locations are accurately known. This assumption may be sometimes invalid in modern practical scenarios; for example, in wireless sensor networks (WSNs), the sensors are frequently localized level by level [19]. Due to this level-by-level procedure, the sensor locations generally include random errors in which case the sensor location errors need to be taken into consideration [20]. In this paper, we propose an efficient bi-iterative source position and propagation speed estimation method using TDOA measurements and compare its performance with the LS solution [13], four-step method [14], three-step weighted least-squares (WLS) method [15], and the CRLB. Our contributions include the following.
(1) We develop an efficient bi-iterative method for the source position and propagation speed estimation. (2) The asymptotic convergence of the new bi-iterative estimation method is analyzed.
The proposed solution is shown to approach the CRLB for low Gaussian TDOA measurement noise levels. (4) We study the CRLB of the source position and propagation speed with sensor location errors and extend the bi-iterative method to the more general case where the sensor locations are subject to random errors.
The remainder of the paper is structured as follows. Section 2 describes the problem of computing source position and propagation speed under a Gaussian TDOA noise model. Development of the new bi-iterative estimation method is presented in Section 3. Besides, the asymptotic convergence and approximate efficiency of the bi-iterative method are analyzed under low noise level conditions. Section 4 extends the bi-iterative estimation method to the general situation where the sensor locations are subject to random errors. Simulation studies are given in Section 5 to validate the theoretical analysis and to test our bi-iterative estimation method's performance. Section 6 concludes the paper.
Throughout the paper, we utilize boldface lowercase letters to stand for vectors and boldface capitals to represent matrices. A and A denote the th to th columns and the th column of matrix A, respectively. A ⪰ B means that A − B is positive semidefinite. 0 × stands for the × zero matrix and 1 denotes the ×1 vector of ones. Moreover, diag( * ) denotes a diagonal matrix, ‖ * ‖ 2 represents Euclidean norm, and [ * ] stands for the expectation operator.

Problem Formulation
Consider a two-dimensional plane where a sensor array is composed of sensors at exactly known locations s = [ , ] , = 1, 2, . . . , . It is straightforward to generalize the methodology to three-dimensional case. We postulate that all these sensors are not lying on a straight line, which guarantees that the matrix A in the proposed method has full column rank. Let u = [ , ] be the position of the source which emits signal in free-field conditions to be determined.
Denote by the true distance from the source to the th sensor: Without losing generality, we set the first sensor to be the reference sensor. The TDOA measurement obtained from sensor pair s and s 1 would be [20], after taking into account the measurement noise, where 1 = 1 / , 1 = − 1 , and is the unknown propagation speed. 1 is the additive measurement noise. For simplicity, we presume that all the TDOA measurements used for source localization come from line-of-sight propagation.
Equation (2) can be arranged into matrix form where t = [ 21 , 31 , . . . , 1 ] and n = [ 21 , 31 , . . . , 1 ] . n is assumed to be Gaussian distributed with zero-mean and known covariance matrix Q . This assumption has been frequently adopted in literatures such as [14,15,18] for source location and propagation speed estimation. Notice that TDOA measurement equations (2) are not linearly related to the source position u and propagation speed . The problem of interest is to estimate the source position u and propagation speed as accurately as possible via exploiting measurement equations (3).

Bi-Iterative Source Position and Propagation Speed
Estimation Method. Achieving high source position and propagation speed estimation accuracy is far from being straightforward due to the nonlinearity implied in the TDOA localization equations. According to (2), the TDOA measurement equation can be rearranged as̃1 + 1 = + 1 . Squaring both sides of̃1 + 1 = + 1 and considering (1), we obtain By stacking (4) for = 2, 3, . . . , together, we have where G = [s 2 −s 1 , s 3 −s 1 , . . . , s −s 1 ] , h = 0.5[‖s 2 −s 1 ‖ 2 2 , ‖s 3 − s 1 ‖ 2 2 , . . . , ‖s − s 1 ‖ 2 2 ] , and the symbol ⊙ represents the Schur product operator (element-by-element multiplication). On the right-hand side of (5), = Bn + 2 n ⊙ n is the equation error vector and B is defined as B = diag( The source position u and propagation speed can be achieved by minimizing the following WLS problem: . When the TDOA measurement noises 1 , = 2, 3, . . . , , are relatively small, the second-order noise terms can be ignored and the weighting matrix can be approximated by W ≈ (BQ B ) −1 . It should be noted that the source position u and propagation speed are coupled together by 1 . In order to eliminate this coupling effect, we herein follow the bi-iterative scheme utilized in [22] to alternately estimate the source position and propagation speed.
Firstly, let us consider the situation where the source position u is fixed. From (5) and (6), we have where a 0 = G(u − s 1 ) − h, a 1 = 1t , and a 2 = 0.5t ⊙t. Thus, the cost function (u, ) is a quartic function with respect to the propagation speed . By taking the derivative of (u, ) with respect to and equating the result to zero, we obtain where the coefficients are 0 = a 1 Wa 0 , 1 = 2a 2 Wa 0 + a 1 Wa 1 , 2 = 3a 2 Wa 1 , and 3 = 2a 2 Wa 2 . The solution of the cubic polynomial equation in (8) can be obtained efficiently utilizing Cardano's method [23]. Secondly, consider another situation where the propagation speed is fixed; we arrive at where H = [G,t], y = [(u − s 1 ) , 1 ] , and h = h − 0.5 2t ⊙t. It is noticed that (9) is nonlinear and nonconvex with respect to u since (u − s 1 ) is related to 1 through (1). Here we apply the semidefinite relaxation (SDR) technique [24] to convert (9) to an approximate but convex problem. The minimization problem can be reexpressed as where y(1 : 2) and y(3) denote the first to second elements and the third element of vector y, respectively. Expressing and writing the cost function in (10) into where trace( * ) denotes the trace operator, F =  (10) can be equivalently written in the following form: The minimization problem (12) is just as difficult to solve as problem (10) because of the nonconvexity in the second equality constraint. Thus, replacing the equality constraint Y = yy by an inequality Y ⪰ yy to meet the convex specification and establishing the equivalence between Y ⪰ yy and [ Y y y 1 ] ⪰ 0 4×4 utilizing the Schur complement [24], we have the following SDP: That is to say, the nonconvex optimization problem in (9) is transformed into a convex one. The optimal value of (13) is less than or equal to the optimal value of (12), since we minimize the same cost function over a larger set. The global optimum can be obtained by applying existing numerical methods, for instance, the interior-point methods. Based upon the above analysis, the procedures of the proposed bi-iterative method for source position and propagation speed estimation are given as follows.
(1) Obtain the initial source position and propagation speed estimates u 0 and 0 using the solution in [13].
(2) Put −1 into (9) and obtain the source position estimate u by solving the convex optimization problem (13) utilizing the convex optimization toolbox CVX [24].
(3) Substituting u into (8) produces a cubic equation in , and the roots of this cubic equation can be obtained efficiently utilizing Cardano's method. We select the one that gives the smallest (u , ) in (7).

Asymptotic Convergence.
Here, we show that the biiterative estimation method has asymptotic convergence. It is apparent that the cost function (u, ) given in (6) is continuous because it is differentiable. From the definition of (u, ), we have 4 International Journal of Distributed Sensor Networks Using the fact that A has full column rank and applying the Cauchy-Schwartz inequality [25] to the second term in the right-hand side of (14), we obtain where min is the smallest eigenvalue of A WA. From (15), it can be concluded that the set {u, | min ‖ ‖ 2 2 − 2‖ ‖ 2 ‖A Wh‖ 2 + h Wh ≤ (u, ) < } is bounded for any finite positive constant . In addition, steps (2) and (3) lead to Thus, for = 1, 2, . . ., the cost function reduces or stays the same. As a result, the cost function (u, ) is the so-called Lyapunov function [26], which indicates that the proposed bi-iterative method has asymptotic convergence.

Performance Analysis.
We evaluate the bias and covariance matrix of the source position and propagation speed estimatesû and̂through using the first-order perturbation method [27]. From (16), when the iteration procedure starting from appropriate initial source position and propagation speed estimates is terminated, the gradient of (u, ) atp = [û ,̂] satisfies where J = [G +t u,s 1 ,t ⊙t + 1t ] is the ( − 1) × 3 Jacobian matrix and a,b = (a − b)/‖(a − b)‖ 2 denotes a unit vector. A, u, and̂can be expressed as A = A + ΔA,û = u + Δu, and = + Δ , respectively, where Δu and Δ are the estimation bias inû and̂. A is obtained by replacingt with t in A, and ΔA is the deviation of A from A .
Ignoring the perturbations that are of higher order than 1, we have where the fact that ΔA 12 = 0 ( −1)×2 is utilized. Assigning Δp = [Δu , Δ ] and preserving only the linear perturbation terms, we obtain from (17) where J = [A 12 + A 4 u,s 1 , 2 A 3 + 1 A 4 ]. Then the estimation bias Δp can be approximated as We have [Δp] ≈ 0 3×1 by using the assumption that n is a zero-mean random vector, which indicates that the source position and propagation speed estimates are approximately unbiased for low noise levels. Recalling the definition of W in (6) we have Note that (21) is identical to the CRLB derived in [14], which means that the theoretical performance of the bi-iterative method closely meets the CRLB. The estimation bias and covariance matrix derived here are valid only if the TDOA measurement noises are relatively small.

Algorithm Extension
In practical environments, the sensor locations need to be measured and generally include random errors. In this section, we first examine the increase in the CRLB of the source position and propagation speed due to sensor location errors and then extend the bi-iterative method to the situation where the sensor locations are subject to random errors. Let s = s + Δs denote the erroneous location of the th sensor, where s is the true sensor location and Δs is the random location error, = 1, 2, . . . , . Arrange these sensor location errors Δs 1 , Δs 2 , . . . , Δs into matrix form as Δs = [Δs 1 , Δs 2 , . . . , Δs ] . Δs is assumed to be Gaussian distributed with zero-mean and covariance matrix Q s = [ΔsΔs ] and to be mutually independent of the measurement noise vector n.

CRLB.
We obtain the CRLB of p = [u , ] with sensor location errors by following the same approach as in [20]: where  where ( t/ p ) and ( t/ s ) are ] .
The term X −1 in (22) is the CRLB of p when there exists no sensor location error [14]. The expression (22) where the second-order error terms have been neglected. It is noteworthy that the variable 1 in (25) depends on the true sensor position s 1 . Hence, expand 1 around the erroneous sensor location s 1 and preserve only the first two terms as where 1 = ‖u − s 1 ‖ 2 . Inserting in (25) Stacking (27) for = 2, 3, . . . , gives Then we have the following cost function: where the weighting matrix is defined as W = (BQ B + DQ D ) −1 .
Once the cost function (u, ) is obtained, we can use the same bi-iterative scheme presented in Section 3.1 to compute the source position and propagation speed.

Simulation Study
In this section, a set of Monte Carlo simulations are carried out to assess the bi-iterative method's performance by comparing with the LS solution [13], the four-step method [14], the three-step WLS method [15], the CRLB, and the theoretical bias [28]. We simply follow [15] to consider a tangible acoustic interface application of interactive displays [29]. = 2000 is the total number of Monte Carlo trials.

Localization Performance without Sensor Location Errors.
We first examine the performance of the bi-iterative estimation method when there exists no sensor location error.  small as 10 −6 . The remaining three methods closely approach the CRLB for both source position and propagation speed estimates when ≤ 4.6 × 10 −6 , which verifies the analysis in Section 3.3 on the approximate efficiency of the bi-iterative method. As further increases, the bi-iterative, three-step WLS and four-step methods deviate from the CRLB. This is the so-called threshold phenomenon, which occurs in nonlinear estimation when the noises exceed a threshold value [8]. Note that the estimation RMSEs of the bi-iterative method are below the CRLB for high noise levels. A reasonable explanation for this is that the bi-iterative method gives biased estimates for high noise levels. It is possible that when an estimate is biased, its RMSE can be smaller than the CRLB [21]. Figure 2 compares the estimation biases of the bi-iterative method with the remaining three methods and the theoretical bias [28] for source at u = [0.2, 0.1] m. As expected, the source position and propagation speed estimation biases of the bi-iterative method are all reasonably small for low noise levels, which illustrates the unbiasedness of the bi-iterative method for small measurement noises. As increases, the estimation biases of the bi-iterative method deviate from  the theoretical bias apparently. The bi-iterative method is superior to the remaining three methods. Figures 3 and 4 show the estimation RMSE results for two sources at u = [0.5, 0.1] m and u = [0.2, 0.5] m, respectively. Each of these two sources has one coordinate being identical to that of the reference sensor. More specifically, the -coordinate of [0.5, 0.1] and the -coordinate of [0.2, 0.5] are equivalent to the corresponding coordinates of s 1 . Clearly, the estimation RMSEs of the bi-iterative and fourstep methods closely meet the CRLB for ≤ 10 −4 , while those of the LS and three-step WLS methods are significantly larger than the CRLB. The most likely reason for this is that a matrix ill-conditioned problem occurs in computing the weighting matrix, which indicates that the LS and three-step WLS methods are more sensitive to the source positions. In the following, we evaluate the computational complexity of the bi-iterative method. In step (2), the complexity of calculating propagation speed is about (2 + 6 2 ), where is the number of sensors and is the dimension of the source position. In step (3), the complexity induced by calculating F is roughly (2( 2 + 2 )). Denote by and V the number of equality constraints and the optimization problem size, respectively; the worst case computational complexity of solving SDP is ( In the SDP (13), = 2 and V = + 2. Therefore, the total complexity is about (2( 2 + 2 ) + ( 3 The complexity is (2( 2 + 2 ) + (6 + 7) 2 ) for the four-step method. It is clear that the bi-iterative and four-step methods have comparable computational complexity since we have = 2 or 3.

Localization Performance with Sensor Position Errors.
Here, we investigate the CRLB of source position and propagation speed and evaluate the bi-iterative method's performance with sensor location errors. The sensor location errors at different sensors are set to be unequal to make the comparison more general, where Q = 2 diag( [4,4,2,2,8,8,20,20,10,10]). Fix at 10 −5 and Figure 7 depicts the trace of the upper-left 2 × 2 submatrix and the third diagonal element of CRLB(p) as increases for source at u = [0.2, 0.1] m. The CRLBs of source position and propagation speed without sensor location errors are also plotted for comparison purpose. It is clear that the deviations of the CRLB of source position and propagation speed with sensor location errors from the case without sensor location errors become larger and larger as increases.
We set to be 10 −2 and check the influence of TDOA measurement noise on the estimation accuracy in the presence of sensor location errors. Figure 8 plots the RMSE results of the bi-iterative and LS methods as increases for source at u = [0.2, 0.1] m. Again, the LS solution cannot achieve the CRLB even when = 10 −6 , while the biiterative method approaches the CRLB until it suffers from the thresholding effect at = 4.3 × 10 −5 . A comparison of Figures 1 and 8 reveals that a small error in the sensor locations can significantly degrade the estimation accuracy when the TDOA measurement noise is relatively small.
We next keep at 10 −5 and investigate the effect of on the source position and propagation speed estimation accuracy. Figure 9 plots the RMSE results against for the source at u = [0.2, 0.1] m. Clearly, neither the source position nor propagation speed estimation accuracy in the LS solution can achieve the CRLB performance even when is as small as 10 −3 , while the bi-iterative method closely approximates the CRLB when is less than 3.6 × 10 −2 .

Conclusions
The problem of estimating the source location and the propagation speed using TDOA measurements was considered in this paper. Using the bi-iterative scheme, an efficient method that alternately computes the source position and propagation speed is developed. The asymptotical convergence of the proposed bi-iterative estimation method is analyzed. The proposed method is shown by first-order perturbation analysis to approach the CRLB for low noise levels. In contrast to other existing estimation methods, the threshold effect of the new bi-iterative method occurs later as the measurement noise increases. Simulation studies are given to illustrate the effectiveness of the proposed bi-iterative method.