Interval analysis-based Bi-iterative algorithm for robust TDOA-FDOA moving source localisation

In this study, an interval extension method of a bi-iterative is proposed to determine a moving source. This method is developed by utilising the time difference of arrival and frequency difference of arrival measurements of a signals received from several receivers. Unlike the standard Gaussian noise model, the time difference of arrival - frequency difference of arrival measurements are obtained by interval enclosing, which avoids convergence and initialisation problems in the conventional Taylor-series method. Using the bi-iterative strategy, the algorithm can alternately calculate the position and velocity of the moving source in interval vector form. Simulation results indicate that the proposed scheme significantly outperforms other methods, and approaches the Cramer-Rao lower bound at a sufficiently high noise level before the threshold effect occurs. Moreover, the interval widths of the results provide the confidence degree of the estimate.


Introduction
Passive source localisation based on time difference of arrival (TDOA) and frequency difference of arrival (FDOA) measurements has become a topic of interest owing to its wide application in radar, 1 navigation, 2 interference localisation, 3 and wireless networks. 4 For the high nonlinearity implied in the measurement equations, many available methods for solving TDOA-and FDOA-based moving source localisation have been proposed. The traditional method in Foy, 5 and Lu and Ho 6 linearises the equations via Taylor-series expansion. Two drawbacks must also be considered. First, an initial solution is required to obtain the source estimate. Second, the convergence is not guaranteed. To avoid these issues, researchers have investigated algebraic methods, which allow the independence of the initial estimate and have high computationally efficient. 7 Among these methods, the two-step weight least-square (TSWLS) method, 8,9 is popular, and can provide an algebraic solution without an initial guess. Interestingly, the closed-form result of the TSWLS reaches the Cramer-Rao lower bound (CRLB) at the low-noise level. However, it shows poor localisation accuracy when the target approaches the axis of the reference sensor. The method discussed in, Zou and Liu 10 and Zhou, et al. 11 is proposed to transform the maximum likelihood estimator problem into a convex optimisation problem, but its results are more robust at high computational complexity. The CTLS algorithm 12 overcomes the problem that the total least squares algorithm cannot achieve the optimal solution when the noise components are statistically correlated, and can achieve CRLB when the measurement noise is moderate. However, due to the influence of the target and the location distribution of the external radiation source, the coefficient matrix may have ill-conditioned problems, resulting in huge fluctuations in the results of the equation solution due to small observation errors, and degrading the performance of the CTLS algorithm. The bi-iterative technique recently proposed in, Zhu and Feng 13 and Zhu et al. 14 can reduce the computational cost by alternately calculating the source location and velocity. The mentioned methods assume that the TDOA and FDOA noises both accord with the Gaussian distribution. They only provide the point estimates of the source position and velocity, and not a confidence interval.
From the Taylor-series expansion, this work derives the interval extension of the bi-iterative method for a moving source localisation based on the bounded error framework. The proposed algorithm combines the interval analysis technique in, Abdallah et al., 15 and Jaulin and Walter 16 with a bi-iterative strategy 13,14 to estimate the target position and velocity. Compared to the traditional iterative method, the novel scheme obtains estimations of the source in interval vector form, which has a 95% probability 17 of containing the true vales. It allows global convergence and avoids the initialisation problem. Furthermore, the interval widths of the results provide the estimated confidence.
The rest of this paper is organised as follows: Section 2 defines the symbols and notations. Section 3 discusses the measurement model. The details of the interval analysis-based bi-iterative algorithm are shown in Section 4. Section 5 deduces the CRLB and the mean square error (MSE) of the source location and velocity estimates. Section 6 provides the simulation results, and valuable conclusions are presented in Section 7.

Symbols and notations
We denote the punctual column vectors and matrices in bold lower and upper-case letters, respectively. For n 3 1 vector x, x i , i = 1, 2, :::, n, represents the ith component. For n 3 m matrix A, A i, j is the entry in the ith row and jth column, where i = 1, 2, :::, n and j = 1, 2, :::, m.
In the interval number system, ½x denotes the closed set ½x, x = fx 2 R, x\x\ xg and is referred to as an interval scalar ½Ã.
x and x are called the lower and upper endpoints of ½x, respectively. If x = x, then ½x degenerates to a scalar (i.e. in this case, ½x = x). The midpoint of ½x is equal to mid(½x) = 1 2 (x + x), and its width is defined as w(½x) = x À x.
To operate on the interval, we extend the basic calculation of the floating-point numbers and the set operations, such as + , À, Ã, 4, \ and [, into the interval analysis. 18 The interval correspondence of functions is usually impossible to calculate. Thus, the concept of inclusion function is proposed. The inclusion function ½f (Á) for a function f (Á) indicates that the image of interval ½x by ½f (Á) is an interval containing f (½x). 19

Measurements model
We consider a scenario of M (M.4) moving or stationary sensors that collaborate to determine a moving target with an unknown position u = ½x, y, z T and velocity _ u = ½_ x, _ y, _ z T using TDOA and FDOA measurements in three dimensional space. Sensor position s i = ½x i , y i , z i T and velocity _ s i = ½_ x i , _ y i , _ z i T , i = 1, 2, :::, M, are accurately derived by an estimator. Without loss of generality, the first sensor is set as the reference.
The accurate TDOA measurement between sensor pair i (i = 1, 2, :::, M) and 1 multiplied by the speed of propagation, known as the range difference of arrival, is is the distance from the target to the ith sensor, and c is the signal propagation speed. Similarly, the FDOA measurement between sensor i (i = 1, 2, :::, M) and the reference sensor multiplied by the speed of propagation is is the time derivative of (2), and f 0 is the carrier frequency.

Localisation algorithm
Let u = ½u T , _ u T T be the unknown parameter of the true location and velocity of the moving source, using TDOA-FDOA measurements to estimate u that can be expressed as where Dm = ½Dd T , D _ d T T is an approximately zeromean Gaussian noise and its covariance matrix is where Q d is the covariance matrix for TDOA measurement noise and Q _ d is the covariance matrix for FDOA measurement noise. Noise vectors Dd and D _ d are assumed to be independent with one another. Original interval vectors ½u and ½ _ u can be obtained from prior knowledge, such as the sensors' coverage area and the maximum source speed. Based on the relationship of the TDOA-FDOA measurements in equations (1) and is assumed as the initial position solution, equation (8) can be solved by the Taylor-series expansion of m (u) around u m as shown in In addition, Jacobian matrix J 1 is shown in and the derivation of partial derivative where the (i-1)th rows of (∂d=∂u) and (∂ _ d=∂u) are given by, i = 2, :: In accordance with equation (10), the least square (LS) estimator of u is shown in where Q r is the covariance matrix for TDOA measurement noise. Given that m o 2 ½m o andm 2 ½m, ½m o and ½m are used instead of m o and ½m in (15), we obtain The inverse of (½J 1 T Q À1 r ½J 1 ) is required to solve equation (16), which increases the computational cost. As for the nonlinearity implied in interval analysis, the idea of a midpoint test is proposed in Jaulin and Walter. 18 Essentially, reducing the computational complexity is effective. Therefore, we simplify ½J 1 with J 1, m = mid(½J 1 ) in the case of lower location errors. The source position interval vector ½u can be solved using the following equation Second, we consider another situation where source position u 2m = mid(½u) is fixed, ½u is updated by equation (17), and m (u) is expanded around _ u m = mid(½ _ u), as shown in m andm are provided with the same functional form, Jacobian matrix J 2 is shown in and partial derivative ( According to equation (18), we can obtain the LS estimator of the source velocity, as shown in where Q _ r is the covariance matrix for FDOA measurement noise.
Similar to the solution process of ½u, m o is first replaced with ½m o in equation (23) to obtain the following formula Then, J 2, m = mid(½J 2 ) is used to approximate ½J 2 in equation (24), and ½ _ u can be obtained, as shown in On the basis of the above description, ½u and ½ _ u can be alternately substituted in equations (17) and (25). Repeating the two steps can guarantee the convergence of the proposed algorithm because of the convergence of the Taylor-series algorithm near the solution. Moreover, instead of estimating the target position and velocity simultaneously, the developed method can alternately calculates the target position and velocity, which can simplify the derivative where the Jacobian matrix J is shown in and the derivation of the partial derivative ( Thus, the position estimation accuracy of the proposed method is nearly identical to the CRLB under low noise conditions.

Performance analysis
In this subsection, we deduce the MSE of the source location and velocity estimates through the Taylorseries method at low TDOA and FDOA measurement noise levels. If the iteration procedure in the localisation algorithm is terminated, then we can obtain the one-order Taylor-series expansion 21 at approximately u = ½u T , _ u T T by combing equations (17) and (25). The equation can be obtained as shown in where u m = ½u T m , _ u T m T is the final midpoint of position interval vector ½u and velocity interval vector ½ _ u, and J m = ½J 1, m , J 2, m is the approximate values of ½J 1 and ½J 2 under low noise conditions. Let Du = uÀu m denote the estimation errors of the source position and velocity, which is a zero-mean random vector. 22

Equation (29) can be transformed into
The covariance matrix of Du can be approximated by Thus, the position and velocity estimation accuracy of the proposed method is nearly identical to the CRLB under low noise conditions. When the noise increases, the developed algorithm cannot guarantee the high accuracy of position and velocity estimation results, but the interval results can still estimate the range of the target position and velocity. Therefore, the proposed method provides anther idea for target localisation in high noise conditions.

Simulation
In this section, a set of Monte Carlo simulations evaluate the performance of the proposed algorithm by comparing it with TSWLS 8 and the maximum likelihood (ML) estimator 23 and CTLS. 12 The localisation accuracy is evaluated in terms of the root mean square error (RMSE) and bias of the source location and velocity. The estimation bias in terms of the norm of estimation bias is defined as  Table 1. The TDOA and FDOA measurements are generated by adding zero-mean Gaussian noises with the covariance matrix Q = diag(d 2 Q r , 0:1d 2 Q r ), where Q r = 0:5(I + I T I), and I is a unit array. In this simulation, we consider nearand far-field scenarios in the moving source localisation. The estimation bias and accuracy of the source are investigated with the increase in TDOA and FDOA measurement noise.

Near-field scenario
In the near-field case, the true position and velocity of the target are u = ½285, 320, 275 T m and _ u = ½À20, 15, 40 T m Á s À1 , respectively. 13 The noise level d 2 varies from -30 in log scale to 40 in log scale. We set the original interval vectors to ½u = ½½À500, 500, ½À500, 500, ½À500, 500 T m and ½ _ u = ½½À50, 50, ½À50, 50, ½À50, 50 T m Á s À1 . Figure 1 shows that the estimation bias norm of the proposed method is significantly smaller than those of TSWLS, ML and CTLS, especially for high measurement noise levels. Figure 2 plots the accuracy of the position and velocity estimates of TSWLS, ML, CTLS and the proposed method, in terms of the RMSE with increasing d 2 . The three methods can approach the CRLB at a low noise level. For the source position estimation, the threshold effect of the proposed method occurs at 28 dB, which is later than the others. In the velocity estimation, TSWLS, ML and CTLS also provide inaccurate estimates that are apparently earlier than that of the  1  300  100  150  30  220  20  2  400  150  100  230  10  20  3  300  500  200  10  22  1 0  4  350  200  100  10  20  30  5 2100 2100 2100 220 10 10 proposed algorithm. For ML, it is very complicated to obtain initial estimation of position and velocity when noise increases. Excessive noise affects the construction of the weight matrix in TSWLS and CTLS, resulting in a large deviation of the estimated value.

Far-field scenario
In the far-field case, we set the true target position u = ½1000, 1500, 2000 T m and velocity _ u = ½À20, 15, 40 T m Á s À1 , 10 whereas the noise level d 2 varies from -30 in log scale to 20 in log scale. The original interval vectors of the moving source position and velocity are ½u= ½½À3000, 3000, ½À3000, 3000, ½À3000, 3000 T m and ½ _ u=½½À50, 50, ½À50, 50, ½À50, 50 T m Á s À1 , respectively. Figure 3 depicts that the estimation bias norms of the four methods in the far-field case are more unstable than those in the near-field case, but the location and velocity estimation biases of the proposed algorithm    are relatively small even at the high measurement noise level.
As shown in Figure 4, the accuracy of the position and velocity estimates of TSWLS, ML, CTLS and the proposed method also decrease in terms of RMSE with increasing d 2 in the far-field case. The three methods can reach the CRLB at a low noise level, but their threshold effects occur earlier than those in the nearfield case because the location can be uniquely evaluated to a single coordinate point when the source is near the sensors where the wavefront is curved. 24,25 However, if the target is far from the sensors, the methods ignore the curvature of the wavefront. In the source position estimation, the threshold effect of the proposed method occurs at 25 dB, which is also later than the others. For the velocity estimation, the performance of the proposed algorithm is similar to that of the location estimation.

Confidence probability
For source location and velocity estimations, the proposed method can not only provide a point estimate, but also the confidence interval. When the simulation scenario and parameter settings are consistent with the near-field case, the algorithm provides the interval vector of the position and velocity of the source, that is., ½u and ½ _ u, where they include their true values. Confidence probability measures the probability that the interval vector contains the true solution and is defined as n N , where N is the total number of independent trials, and n is the number of trials whose interval vectors include the true position and velocity of the source. As shown in Figure 5, the proposed algorithm can guarantee that the true source position and velocity are contained in the interval vector when d 2 is lower than 20 in log scale.  The specific boundaries of the location and velocity estimates of the proposed method in the X-axis are also shown in Figure 6. When d 2 is lower than 30 in log scale, the boundaries of the interval vector also remain the same. This result indicates that the proposed algorithm can provide a reliable range of position and velocity estimates at low and moderate noise levels.

Computational complexity
In this paper, the multiplicative times of four algorithms are used as a measure of computational complexity. As shown in Table 2, TSWLS and CTLS have the same computational complexity O (N 3 ) . In contrast, the computational complexity of the proposed algorithm and ML algorithm are O ((N iter + 1)N 2 ) , where N iter is the number of iterations. After N = 5000 independent simulation experiments, the average running time of the proposed method, ML, TSWLS and CTLS are 3:9 3 10 À3 s, 4:9 3 10 À3 s, 2 3 10 À3 s and 4:5 3 10 À2 s, respectively. The results show that the proposed method can reduce computation complexity in comparison with the traditional ML estimator and CTLS. Furthermore, the proposed method provides a more complex solution compared with TSWLS, but it is superior in accuracy performance.

Conclusion
In this study, we consider the problem of estimating a moving source using the TDOA-FDOA measurements obtained from multiple sensors based on the bounded error framework. By combining interval analysis with a bi-iterative strategy, we develop an efficient method that alternately calculates the source position and velocity interval vectors that enclose the true values. Simulation results show that the algorithm has superior performance over other methods and approaches the CRLB at a sufficiently high noise level. ML: maximum likelihood; TSWLS: two-step weight least-square; CTLS: constrained total least-square.