A Novel Estimator for TDOA and FDOA Positioning of Multiple Disjoint Sources in the Presence of Calibration Emitters

Multiple-target localization is extensively applied in wireless connected networks. However, sensor location uncertainty is known to degrade significantly the target localization accuracy. Fortunately, calibration emitters such as unmanned aerial vehicles (UAV) with known location can be used to reduce the loss in localization accuracy due to sensor location errors. This paper is devoted to the use of UAV calibration emitters for time differences of arrival (TDOA) and frequency differences of arrival (FDOA) positioning of multiple targets. The study starts with deriving the Cramér–Rao bound (CRB) for TDOA/FDOA-based target location estimate when several UAV calibration signals are available. Subsequently, the paper presents an iterative constrained weighted least squares (ICWLS) estimator for multiple-target joint localization using TDOA/FDOA measurements from both target sources and UAV calibration emitters. The newly proposed method consists of two stages. In the first phase, the sensor locations are refined based on the calibration measurements as well as the prior knowledge of sensor locations. The second step provides the estimate of multiple-target locations by combining the measurements of target signals as well as the estimated values in the first phase. An efficient ICWLS algorithm is presented at each stage. Both the two algorithms are implemented by using matrix singular value decomposition (SVD), which is able to provide a closed-form solution and update the weighting matrix at every iteration. Finally, the convergence behavior and estimation mean-square-error (MSE) of the new estimator are deduced. Both theoretical analysis and simulation results show that the developed method can improve the TDOA/FDOA localization accuracy obviously with the help of UAV calibration emitters.


I. INTRODUCTION
Determining the target locations from source signal measurements collected by a wireless connected network (or an array) The associate editor coordinating the review of this manuscript and approving it for publication was Bo Zhang . of sensors at a time instance has become a fundamental issue during the past few decades [1]- [15]. It has inspired numerous applications in wireless communications, unmanned system, wireless connected robot swarms, object tracking and many other areas. Unfortunately, target localization is a nontrivial problem, mainly because the target positions and the VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ signal measurements are generally nonlinearly related. The localization task is usually accomplished using a two-step approach. In the first step, certain signal parameters are extracted from the radiated signal at the output of an array of spatially distributed sensors. In the second step, the target location is estimated by solving a set of non-linear equations defined by the signal parameters obtained in the first step.
In this work, we restrict our attention to the second step. Typical positioning parameters include time of arrival (TOA) [16], time difference of arrival (TDOA) [17], angle of arrival (AOA) [18], received signal strength (RSS) [19], and a combination of them [20]. When there is relative motion between the target and the sensors, frequency differences of arrival (FDOA) can be exploited not only to further improve the estimation accuracy of the target position, but also to provide an estimate of the target velocity. Besides, combinations of TDOA and FDOA measurements can facilitate localization of the target source in a wide range of bandwidths. As a consequence, FDOA is often combined with TDOA for moving target localization. However, finding the target position and velocity from TDOA and FDOA measurements obtained at a single time instant is not easy to complete as the equations are nonlinear. Another undesirable factor is that there are often random errors in the positions and velocities of the mobile sensors, which will result in obvious losses in estimation accuracy. To limit the scope of this paper, we concentrate on the TDOA/FDOA-based target localization problem.
Over the years, a number of algorithms for TDOA and FDOA positioning have been presented in the literature. Among the existing methods, the two-stage weighted least squares (TSWLS) algorithm [21] is very popular because it has closed-form solution and asymptotic efficiency. However, this classical TSWLS algorithm has poor localization accuracy when the target is located near any coordinate axis of the reference sensor. To remove this shortcoming, an improved version of the popular TSWLS algorithm is proposed in [22]. An alternative approach that can yield explicit solutions for target position and velocity is based on multidimensional scaling (MDS) analysis [23], which has been developed for data analysis in fields such as physics, geography and biology. This method is designed to minimize the cost function dependent on the scalar product matrix in the MDS framework. Theoretical analysis and simulation results demonstrate that this method is more robust against high measurement noise than the classical TSWLS method. Certainly, the closedform methods have low computational complexity and do not suffer from local minima and divergence problems. However, the accuracy of these estimators can attain Cramér-Rao bound (CRB) only when the noise level is sufficiently low. To improve the estimation performance in the high noise region, iterative localization approaches can be exploited. The most common iterative technique is perhaps the Taylorseries linearization algorithm [24], which can be applied to almost any localization scenario. A point to note is that this method requires a proper initial position and velocity guess close to the true solution, and such a good guess may not be easy to obtain because the pseudo-linear equations are not available in this algorithm. Additionally, the semidefinite relaxation technique has been employed to relax the maximum likelihood (ML) estimator to semidefinite programming (SDP) [25]. The advantage of the SDP-based method lies in the fact that the cost function has no local minima or saddle points and thus the iterative sequence always convergence to the global minimum. However, this kind of method generally has a complex mathematical model and requires high computational complexity for reasonable localization accuracy. An efficient constrained total least squares (CTLS) algorithm for estimating the position and velocity of a moving target is developed in [26], and a constrained weighted least squares (CWLS) algorithm for TDOA and FDOA positioning is presented in [27]. These two estimators show superior performance and can reach CRB accuracy under moderate noise level conditions. However, both methods are implemented based on Newton's iteration. Hence, global convergence cannot be ensured unless the initial guess is close enough to the optimal solution. In [28], a polynomial root-findingbased CWLS algorithm is developed for TDOA/FDOA-based emitter locations. It makes use of the altitude constraint to reduce the localization errors. Unfortunately, the altitude of the emitter may not be accurately known or measured in practice. Recently, an iterative constrained weighted least squares (ICWLS) algorithm has been presented in [29]. It iteratively performs a linearization procedure on the quadratic equality constraints to get approximate programming with linear constraints. As a consequence, a closed-form solution is available at every iteration. Both theoretical analysis and simulation studies reveal that this novel method has better convergence properties and is able to delay the thresholding effect where the performance drops suddenly.
In addition to the estimation errors in TDOA/FDOA measurements, the uncertainties in sensor position and velocity can also cause the location performance to deteriorate dramatically, regardless of the approaches used to determine the solution. Indeed, a slight error in sensor location may lead to a large error in the target location estimate. It must be pointed out that sensor position and velocity errors often occur in real scenarios [30]- [37]. Recently, an emitter localization problem that accounts for uncertainties in sensor location has attracted considerable attention. In [30], [31], the increase in the mean-square error (MSE) in the target location estimate is derived when the sensor positions and velocities are assumed correct, but in fact they are not accurate due to random errors. To obtain the asymptotic optimum performance for this localization scenario, TSWLS and the CWLS estimators are proposed in [30] and [33], respectively. Both methods can achieve CRB accuracy over a small noise region. On the other hand, the problem of locating multiple disjoint sources using TDOAs and FDOAs is also studied in the literature, especially for intelligent robot swarms. Note that the positioning parameters (e.g., TDOA and FDOA) of the disjoint sources can be estimated separately because of the disjointness in time, frequency or both. An important observation about multiple disjoint sources is that they are subject to the same sensor location displacements. Thus it is possible to reduce the performance loss resulting from sensor location uncertainties by joint localization. This procedure can be interpreted as multiple-target cooperative positioning. In [34], a closedform solution is proposed to locate multiple disjoint sources in the presence of random sensor location errors. This solution has been proven analytically to reach CRB performance when the amount of noise is small. However, this estimator cannot provide optimum accuracy for sensor positions and velocities. Following the work of [34], an improved closedform solution for simultaneously locating multiple disjoint sources and refining erroneous sensor positions and velocities is proposed in [35], [36], where the source and the receiver location estimates can asymptotically attain CRB. In addition, an ICWLS estimator for TDOA/FDOA positioning of multiple disjoint sources is developed in [37]. This method can also yield solutions for both target and sensor locations, and moreover, it can tolerate higher noise levels before the thresholding effect starts to happen.
To further optimize the localization performance in the presence of erroneous sensor locations, we can exploit a set of calibration emitters whose positions and/or velocities are known to an estimator in advance. Extensive studies have shown that the use of calibration emitters can significantly reduce the loss in target localization accuracy when the available sensor locations have random errors. In [38], the use of a single calibration emitter with known position is investigated to mitigate the TDOA-based localization errors caused by sensor position uncertainties. The gain in localization accuracy resulting from a calibration signal is evaluated based on the CRB analysis results. Additionally, a closedform solution is presented for target localization using TDOA measurements from both the target source and the calibration emitter. Reference [39] extends the work in [38] to a more practical scenario where the position of a calibration emitter is not known accurately. An asymptotically efficient solution is developed by incorporating statistical knowledge of the calibration position errors. It can be easily observed that the performance gain from a calibration emitter depends on its position. Consequently, the calibration position is optimized in [40] by improving the Fisher information matrix of the target location estimate. The use and geometric effects of multiple calibration emitters are studied in [41], which presents the analytical conditions under which sensor position errors can be completely eliminated in target localization through the use of multiple calibration emitters. A Taylor-series TDOA localization method with calibration emitters in the presence of synchronization clock bias and sensor location errors is proposed in [42]. Note that the above-mentioned approaches only consider a stationary source location using TDOA measurements. For a moving target localization, some efficient methods with calibration signal are also developed in the literature. Specifically, reference [43] develops an explicit solution for moving target location using TDOA and FDOA from both the target source and the calibration emitter. The work in [44] focuses on the use of calibration sensors for locating a moving target based on TDOA/FDOA measurements. The calibration sensors can broadcast calibration signals to other sensors and can also receive signals from the target source and other calibration sensors. It presents a closed-form estimator with low complexity. It must be emphasized, however, that the solutions in [43], [44] are designed for the localization of a single source only. Moreover, their noise thresholds could be further improved because they directly follow the idea behind the classical TSWLS estimator. This paper focuses on developing a novel estimator for the TDOA/FDOA multiple-target positioning for wireless connected network in the presence of calibration emitters. In the localization scenario studied here, the calibration emitters are stalled on the cooperative unmanned aerial vehicles (UAVs), whose locations can be known accurately. In addition, the target sources may be placed on any mobile vehicles, such as aircrafts or mobile robots. The study starts by deriving the CRB for the TDOA/FDOA-based target location estimate when some UAV calibration emitters with known locations are available. The insight gained from the CRB illustrates that UAV calibration emitters are very useful in reducing the effects of sensor position and velocity errors. Moreover, the cooperative gain is considerable when the multiple moving targets are located simultaneously in the presence of UAV calibration emitters, even if the target sources are disjoint. To achieve the optimum estimation accuracy, this work proceeds to develop an ICWLS localization approach using TDOA/FDOA measurements from both target sources and calibration emitters. Similar to [29], the estimator presented here also uses a set of linear equality constraints instead of the quadratic constraints to yield an explicit solution in each iteration. The difference lies in the fact that our method has two stages and can utilizes the calibration measurements to locate multiple targets in the presence of sensor position and velocity errors. Specifically, the first stage uses the calibration TDOA/FDOA measurements as well as the statistical characteristics of the noisy sensor locations to refine the sensor locations. The second stage provides estimates of the multiple-target locations by combining the measurements of target signals and the values estimated in the first step. An efficient ICWLS algorithm is designed at each stage. Both algorithms are implemented using matrix singular value decomposition (SVD), which can be achieved through some simple and efficient numerical techniques. Moreover, in the proposed iterative algorithms, a closed-form solution is available at every iteration and the weighting matrix can be updated recursively. The theoretical performance of the proposed ICWLS estimator is also studied. The performance analysis covers convergence property and MSE of the estimation. The results indicate that the proposed method, if it converges, can yield the optimal solution to the formulated non-convex CWLS problem. Besides, its estimation MSE can achieve CRB with UAV calibration emitters at moderate noise levels. Numerical experiments are conducted to confirm the effectiveness of the theoretical analysis and demonstrate the advantages of the proposed ICWLS method.
The remainder of this paper is organized as follows. In Section II, the localization scenario is described and the measurement model is formulated. The CRB expressions for parameter estimation are derived in Section III. Section IV deduces the pseudo-linear equations. The quadratic equation constraints are presented in Section V. Section VI develops a novel ICWLS localization method by the use of UAV calibration emitters. The theoretical analysis of the proposed estimator is derived in Section VII. Simulation results are reported in Section VIII. The conclusions are drawn in Section IX, and the proofs of the main results are given in the Appendices.
In this paper, lowercase and uppercase boldface letters are used to denote vectors and matrices, respectively. In addition, the notational conventions and matrix identities that are used throughout this paper are listed in Tables 1 and 2, respectively.

II. MEASUREMENT MODEL AND PROBLEM FORMULATION A. MEASUREMENT MODEL FOR TARGET SOURCE
This paper considers a three-dimension (3D) localization scenario based on wireless connected network as shown in Fig.1. It consists of a total of D moving targets whose positions and velocities are to be determined. The target sources may be stalled on any mobile vehicles, such as aircrafts or mobile robots. The calibration emitters can be placed on the cooperative UAVs, whose locations are obtained accurately. In addition, the sensors are also moving, and their locations can not be known accurately. Let u d andu d be the true position and velocity of the dth target, respectively. Then, the location vector for the dth target source is denoted asū In order to implement multiple-target cooperative positioning, we need to define a high-dimensional location vector asū = [ū T 1ū T 2 · · ·ū T D ] T , which is found through the use of M moving sensors. The actual position and velocity of the mth sensor are denoted by s m andṡ m , respectively. So, the location vector of the mth sensor is represented as T m ] T , and the composite sensor location vector can be written ass = [s T 1s T 2 · · ·s T M ] T . The TDOAs and FDOAs are exacted from the received signals to locate the multiple targets simultaneously. Notice that this work focuses on the case where the target sources are disjoint so that for It is worthy to point out that TDOA and FDOA can be replaced by range difference of arrival (RDOA) and range difference rate of arrival (RDROA), respectively, in a constantvelocity propagation medium. In the sequel, we use RDOA and RDROA to describe the proposed estimator and conduct the performance analysis for easy of presentation. Without loss of generality, let the first sensor be the reference. Then, the RDOA measurement between sensor pair m and 1 from the dth target is given bŷ where ε p,dm is the RDOA error; r dm is the true RDOA; f p,m (ū d ,s) = ||u d − s m || 2 − ||u d − s 1 || 2 is a function with respect toū d ands. The collection of RDOA measurements from the dth target yieldŝ where Taking time derivative of (1), the RDROA measurement between sensor pair m and 1 from the dth target can be written aŝ where ε v,dm is the RDROA noise;ṙ dm is the true RDROA; is a function of bothū d and. Combining all the RDROA measurements from the dth target produceŝ where Putting (2) and (5) together, we have the total measurement vector for target d as follows: T is the composite noise vector that follows zero-mean Gaussian distribution with covariance matrix It must be emphasized that E pv,d characterizes the crosscorrelation between the RDOA and RDROA measurements.
Collecting the RDOAs and RDROAs from all the target sources produceŝ where · · · (f(ū D ,s)) T ] T (10) It is easy to see that the high-dimensional errorε obeys the Gaussian distribution with zero mean and covariance matrix where ξ m is the error inŝ m . For convenience purpose, let us collect {ŝ m } 1≤m≤M to form the 6M × 1 sensor location vector as follows:ŝ whereξ = [ξ T 1 ξ T 2 · · · ξ T M ] T is the sensor location error vector, which is zero-mean Gaussian distributed with covariance matrixĒ B . VOLUME 8, 2020

B. MEASUREMENT MODEL FOR CALIBRATION EMITTER
It is well known that sensor location uncertainties can degrade the target localization accuracy considerably. To reduce the effects of sensor location errors, it is desirable to use some UAV calibration emitters, whose locations are known to an estimator in advance. The number of calibration emitters is set to N and the position and velocity of the nth calibration emitter are assumed to be w d andẇ d , respectively. Moreover, the calibration emitters are also uncorrelated with each other. Similar to (2), the noisy RDOA measurement vector for the nth calibration signal can be written aŝ where ε (c) p,n is the error vector; r Following (5), the RDROA measurement vector for the nth calibration signal is given bŷ Combining (13) and (15) produceŝ wherer (c) ,n ] T , which is modeled as a zero-mean Gaussian random vector with covariance matrix Putting all the N equations in (17) together yieldŝ where Clearly,ε (c) is a zero-mean Gaussian random vector with covariance matrixĒ

C. PROBLEM FORMULATION
The multiple-target cooperative localization problem in this work can be stated as follows: Given the available erroneous sensor locations and the RDOA/RDROA measurements from both target sources and UAV calibration emitters, identify the composite location vector of the multiple targets as accurately as possible.

III. CRB DERIVATION AND ANALYSIS
This section is devoted to the derivation of the CRB on the estimation of the parameters of interest. The obtained results can provide some valuable insights into the performance gain for target localization through the use of UAV calibration emitters. In addition, we also compare the obtained CRB with the one when the target locations are identified separately. The comparison of these two CRBs can indicate the improvement in localization accuracy resulted from multiple-target cooperative positioning.

A. CRB DERIVATION AND ANALYSIS BASED ON ALL THE TDOA/FDOA MEASUREMENTS
In this subsection, the CRB on the covariance matrix of parameter estimation are deduced based on the TDOA/FDOA measurements from both target sources and calibration emitters. In this situation, the observations consist ofr,r (c) and s, and the unknowns includeū ands. Applying the result in [38], the CRB matrix for joint estimation ofū ands is given in (21), as shown at the bottom of the next page, wherē F 1 (ū,s) = ∂f(ū,s) ∂ū T ,F 2 (ū,s) = ∂f(ū,s) ∂s T andF (c) (s) = ∂f (c) (s) ∂s T . Then, using matrix identities (I)-(III) in Table 2 yields Before proceeding, three remarks are made. Remark 1: According to the result in [34], the CRB ofū without calibration emitters can be written as where the subscript ''n'' is used to indicate that this CRB corresponds to the case of no calibration signals. Since A ) −1F (c) (s)) −1 ≤Ē B , it can be seen from (22) and (24) that CRB(ū) ≤ CRB n (ū), which means that improving the best localization accuracy is possible by the use of calibration emitters. In fact, the improvement could be significant at typical measurement noise level, as is illustrated in simulation section.

Remark 2:
Applying the result of [35], the CRB for the estimate ofs in the absence of calibration emitters is given by (25) immediately results in CRB(s) ≤ CRB n (s). Therefore, the best achievable estimation accuracy for sensor locations can also be improved with the help of calibration emitters.
Remark 3: It can be readily seen from (22) that the CRB for the estimate of the dth target location is given in (26), as shown at the bottom of this page. It is noteworthy that this CRB corresponds to the scenario of multiple-target cooperative localization. If each target source is located separately, the CRB of u d becomes where The subscript ''i'' is used to indicate that these CRB matrices are obtained for the case of non-cooperative localization.
. This result reveals that, compared with the noncooperative positioning, the joint localization of multiple targets can achieve higher estimation accuracy. The subsection presents the CRB on the estimation of parameters based on the TDOA/FDOA measurements from the calibration emitters only. In this case, the observations are composed ofr (c) andŝ, and the unknown iss. The CRB matrix fors is denoted as CRB o (s). Then, the CRB matrix is given by (23) and (28) that CRB o (s) ≥ CRB(s). Hence, the TDOA/FDOA measurements from target sources can be exploited to improve the optimum accuracy of sensor location estimates.

IV. PSEUDO-LINEAR EQUATIONS FOR TDOA/FDOA MEASUREMENTS
The aim of this section is to derive the pseudo-linear equations for TDOA/FDOA measurements from both target sources and UAV calibration emitters.

A. PSEUDO-LINEAR EQUATIONS ASSOCIATED WITH TARGET SOURCES
From (1), the RDOA equation of the dth target source can be transformed into The matrix form of (29) can be described as where Here, β 1 (ū d ,s) and β 2 (ū d ,s) should be considered as auxiliary variables, which plays a very crucial role in ensuing the equations are linear. Subsequently, to get the pseudo-linear equations with respect to the RDROA measurements, we take the time derivative of (29), leading to Similarly, the matrix form of (33) can be expressed as where Stacking (30) and (34) yields the following composite matrix equation for the dth target from both RDOA and RDROA measurements To achieve multiple-target cooperative localization, we need to collect all the D equations in (36) together to form the high-dimensional pseudo-linear equation: where

B. PSEUDO-LINEAR EQUATIONS ASSOCIATED WITH CALIBRATION EMITTERS
Using (14), we can transform the RDOA equation of the nth calibration emitter into Eq. (40) can be restated in a matrix form as (41) [see (42), as shown at the bottom of the next page].
Taking the time derivative of (40) leads to the following equations related to the RDROA measurementṡ whose matrix form is given in (45) and (46), as shown at the bottom of the next page.
Combining (41) and (45), we get the following composite matrix equation for the nth calibration emitter where Stacking (47) for n = 1, 2, · · · , N yields the entire matrix equation with respect to the composite unknown vector as follows:

V. QUADRATIC EQUATION CONSTRAINTS FOR TDOA/FDOA LOCALIZATION
This section is devoted to presenting the quadratic equation constraints on vectorst and η. They play important role in the development of the proposed estimator.

A. QUADRATIC EQUATION CONSTRAINTS ON VECTORt
In this subsection, a set of quadratic equation constraints on vectort is given. To this end, we first need to deduce the quadratic equations that vectors satisfy. Since vector t d includes two intermediate variables, namely, β 1 (ū d ,s) and β 2 (ū d ,s), there exist two constraints related to t d . It follows from the functional forms of β 1 (ū d ,s) and Obviously, the element in 1 and 2 is either one or zero. Both matrices are symmetric. From (50), we arrive at It is noteworthy that (51) describes the final quadratic equation constraints ont. The number of constraints is equal to 2D, which is also the number of nuisance variables int.
On the other hand, differentiating both sides of the first equality in (51) with respect toū yields wherē Following an analogous derivation, we get form the second equality in (51) the relation Notice that (52) and (54) are both crucial for the performance analysis in Section VII.

B. QUADRATIC EQUATION CONSTRAINTS ON VECTOR η
This subsection shows the quadratic equation constraints on vector η. It can be observed from the definition of η that there exist four kinds of instrumental variables in η and the number of these variables equals to 2(M +N ). So, we can form 2(M + N ) quadratic equation constraints with respect to η.
Differentiating both sides of (55)-(58) with respect tos leads to . We stress that (59)-(62) are useful for the performance analysis in Subsection VI.B.

VI. PROPOSED ICWLS ESTIMATOR BASED ON TDOA/FDOA MEASUREMENTS
The objective of this section is to present a novel TDOA/FDOA localization approach for multiple targets when a set of UAV calibration emitters with known locations are available. The proposed estimator is based on the pseudolinear equations derived in Section IV as well as the quadratic equation constraints given in Section V. The new method has two stages. In the first stage, the sensor locations are refined based on the UAV calibration measurements as well as the prior knowledge of sensor locations. In the second stage, the multiple-target cooperative localization is performed by combining the measurements form the target signals with the estimated values in the first phase. It is noteworthy that in each stage a novel ICWLS algorithm is developed. The algorithms are implemented by using matrix SVD, which has robust numerical performance. Moreover, the algorithm can provide closed-form solutions and update the weighting matrices accurately in every iteration.

A. STAGE-1 OF THE PROPOSED METHOD
In the first stage, the measurement vectorsr (c) andŝ are combined to estimates. For this purpose, a novel ICWLS estimator is formulated by using the weighted least squares (WLS) criterion, which is asymptotically efficient under Gaussian noise assumption.

1) OPTIMIZATION MODEL FOR ESTIMATION OF SENSOR LOCATIONS
Notice that the functional forms ofĀ (c) (·) andb (c) (·) in (49) are known, but vectorr (c) is not exactly known and only its noisy value (i.e.,r (c) ) is available. To construct the cost function, we introduce an error vector as Applying a first-order Taylor series expansion ofb (c) (r (c) ) and wherē It follows from (67) thatδ (c) is approximately Gaussian distributed with zero mean and covariance matrix For the purpose of exploiting the noisy measurements of sensor locations, we should introduce an augmented error vector asδ it can be seen thatδ (c) approximately obeys a zero-mean Gaussian distribution with covariance matrix Putting (55)-(58), (70) and (73) together, we can formulate the following constrained minimization problem in (74), as shown at the bottom of this page.
Problem (74) is certainly a quadratic programming with 2M + 2N quadratic indefinite equality constraints, which are nonconvex. Hence, (74) can be viewed as a CWLS problem. It does not seem possible to solve (74) in closed form because of its nonlinearity. Therefore, we have to develop an iterative algorithm to obtain the solution of (74). Besides, it is important to point out that the weighting matrix˜ (c) is dependent on the unknown vector η so that, strictly speaking, it should be denoted as˜ (c) (η). However, since the weighting matrix is updated recursively in the proposed iterative algorithm, we denote it as˜ (c) for the sake of brevity.

2) PROPOSED ITERATION ALGORITHM
The aim of this subsection is to present an iterative algorithm to solve the quadratic programming in (74). In optimization problem (74), the Hessian matrix of the quadratic objective function is positive semidefinite; thus the cost function is convex. But, all the equality constraints in (74) are homogenous, indefinite and non-convex, so the key difficulty focuses on the constraints. Inspired by the work of [29], we would like to make the non-convex constraints become a set of linear equality constraints. As a result, a new convex programming can be obtained as an approximation of (74).
For each homogenous quadratic equality constraint in (74), if one of the variable η is replaced with its last iteration's value, then the non-convex quadratic constraints become linear and convex. The resulted programming is a convex approximation of (74). To be more specific, in the k + 1th iteration we can formulate a convex optimization problem in (75), as shown at the bottom of this page. whereη(k) is the estimated result in the kth iteration. An important advantage of the problem (75) lies in that its optimal solution can be analytically expressed. We are now ready to state the main result.
Proposition 1 can be proved by applying Lagrange multiplier method in a direct manner, so it is omitted due to limited space. It should be emphasized that Proposition 1 holds only when matrixÃ (c) (r (c) ) is of full column rank, which means that there are at least two calibration emitters that must be used for target localization. Next, we want to present an alternative solution of (75) by using matrix SVD technique, with no requirement on the number of calibration emitters. Moreover, it requires less computation and has better numerical stability compared with the solution in Proposition 1, due to the advantage of matrix SVD. Proposition 2: Performing matrix SVD on G (c) (η(k)) leads to (c) 1 (k + 1) is a diagonal matrix, whose diagonal elements are the singular values of G (c) (η(k)). Subsequently, performing matrix SVD on (˜ (c) ) −1/2Ã (c) (r (c) )Q  (c) 2 (k + 1) is a diagonal matrix, whose diagonal entries are the singular values of (˜ (c) ) −1/2Ã (c) (r (c) )Q (c) 12 (k + 1). Further, let us define a vector x 1 (k + 1) that satisfies a set of linear equations: and define a vector x 2 (k + 1) that satisfies the linear relationship: (c) Then, the optimal solution of (75) is given bŷ The proof of Proposition 2 is shown in Appendix B. Certainly, whenÃ (c) (r (c) ) has full column rank, the solution in Proposition 2 is consistent with that in Proposition 1. However, if this condition is not satisfied, for example, only a single calibration emitter is present, then only the solution in Proposition 2 can be applied.
Although the solution to problem (75) is an approximation of the optimal solution of the original problem (74), it incorporates 2M + 2N approximate linear equality constraints to improve the estimation performance. An intuitive update strategy is to choose the optimal solution (84) to be the k + 1th iteration's value. However, this update strategy requires minor modification to avoid possible divergence problem. We would like to adopt the update procedure introduced in [29], which determines the current iteration's value using the linear combination of the solution to problem (75) and the estimation result in the last iteration. Besides, since the weighting matrix˜ (c) is dependent on η, it should also be updated at each iteration step. The first stage of the proposed ICWLS method is formally summarized in Table 3.
The following remarks concern the presented algorithm described above.
Remark 5: If the optimization problem is non-convex, it is very crucial to choose a reasonable proper starting point to guarantee that the iterative algorithm converges to the global minimum of the optimization problem. For the proposed algorithm, we can exploit the prior measurementŝ to form the initial guessη(0) in the following way: The simulation results in Section VIII reveal that this initial estimate is precise enough for global convergence. Remark 6: It can be easily seen from step 10 that the updated value in the k +1th iteration is the linear combination of the solution to problem (75) and the estimation result in the kth iteration. Extensive simulation results indicate that as long as w 1 is not close to 1, the final estimation accuracy is not sensitive to the numerical value of w 1 . We set w 1 = w 2 = 0.5 in the simulation section.
Remark 7: It is impossible to obtain the true weighting matrix˜ (c) , although the algorithm updates this matrix at each iteration step. This is because the calculation of˜ (c) requires true measurementr (c) , which is not available in practice. So we have to compute˜ (c) usingr (c) in place ofr (c) . Fortunately, plentiful simulation results show that the estimation accuracy is relatively insensitive to errors in the weighting matrix and the performance degradation is insignificant due to the approximation of the weighting matrix. This observation is consistent with the findings in [45]- [47].

B. PERFORMANCE ANALYSIS FOR STAGE-1
The aim of this subsection is to provide the theoretical performance analysis for the proposed iterative algorithm in stage-1. First, its convergence property is studied based on the optimality conditions in optimization theory. We prove that if the iterative algorithm converges, then it is able to converge to the optimal solution of the original CWLS problem (74). Next, the MSE expression of the proposed solution is derived by exploiting the first-order perturbation analysis. Moreover, the MSE is proved to achieve the corresponding CRB given in Subsection III.B at moderate error level before the thresholding effect occurs.

1) CONVERGENCE ANALYSIS FOR STAGE-1
Here, the convergence property of the algorithm is investigated. Notice that it is impossible to provide a strict proof as to whether the iteration procedure converges or not, because the original problem (74) is not convex. Therefore, we have to simplify the analysis and prove that if the iterative algorithm converges, then it must converge to the minimizer of problem (74). The convergence analysis begins with the second-order sufficient condition for the optimal solution to (74). It is formally described as follows: Moreover, for any vector y belonging to the null space of (G (c) (η opt )) T we have Then, the vectorη opt is the strictly optimal solution to the CWLS problem (74). The proof of Lemma 1 is given in [48]. The scalars 4n } 1≤n≤N are commonly referred to as Lagrange multipliers. Subsequently, we continue describing the first-order necessary condition for the minimizer of problem (75).
Lemma 2: Assuming that the optimal solution to problem (75) is denoted byη opt (k + 1), there exist {λ The proof of Lemma 2 can also be found in [48]. Moreover, {λ    4n (k + 1)} 1≤n≤N are also called Lagrange multipliers. We are now ready to prove that if the iteration sequence {η opt (k + 1)} 0≤k≤+∞ converges, it will converge toη opt , which is the strictly optimal solution to the CWLS problem (74). For this purpose, three propositions are presented in turn.
Proof: From step 10 in Table 3 we have which completes the proof. Proposition 4: If the sequence {η(k + 1)} 0≤k≤+∞ converges, then the sequences {λ    Appendix D shows the proof of Proposition 5, which reveals that if the proposed iterative algorithm converges, then it can provide the optimal solution for problem (74). It should be noted, however, that there is no guarantee at all that the iterative algorithm always converges. Fortunately, the simulation results in Section VIII demonstrate that in most cases the convergence can be achieved after a few iterations.

2) ESTIMATION MSE FOR STAGE-1
The aim of this subsection is to derive the MSE expression of the proposed solution for stage-1. For this purpose, the firstorder perturbation analysis is applied to get the linear relation between the estimation errors and the measurement errors, which consist ofξ andε (c) .
Assume that the iterative algorithm described above converges toη f , which is also the minimizer for (74). According to step 12 in Table 3, the final estimate for sensor locations is given byŝ f = [I 6M O 6M ×(2M +2N ) ]η f . Therefore, in order to get the analytical expression for the MSE ofŝ f , it is necessary to derive the theoretical MSE ofη f first. Let the estimation error inη f be η f . From the theoretical frame of the firstorder error analysis, η f must be the optimal solution of the following constrained optimization problem [28]-[see (92), as shown at the bottom of the next page.] Assuming the matrix SVD of G (c) (η) is given by η f must lie in the range space of Q (c) 12 , which is orthogonal to that of G (c) (η). As a consequence, we can express η f as It can be observed from (94) thatη f is an approximately unbiased estimate of η and, moreover, the estimation MSE ofη f is given by whereP (c) (r (c) ) = (Ã (c) (r (c) )) T (˜ (c) ) −1Ã (c) (r (c) ).
Before proceeding, three remarks are concluded in order.
It is noteworthy that (96) is useful in the proof of Proposition 6 given below. Next, we derive the theoretical MSE of the estimatê s f . Assuming that the estimation error inŝ f is defined as s f , it is straightforward to check that Notice that the solutionŝ f is an asymptotically efficient estimate ofs because it can achieve the associated CRB accuracy under moderate noise level, as shown below. Proposition 6: Under the first-order approximation, the relation MSE(ŝ f ) = CRB o (s) holds.
The proof of Proposition 6 is provided in Appendix E, which reveals that the sensor location estimate obtained in the first stage is asymptotically optimal. The estimator developed in the second stage will exploit the estimation result from stage-1.

C. STAGE-2 OF THE PROPOSED METHOD
In the second phase, we combine the measurement vectorr with the estimateŝ f provided in the first stage to perform multiple-target cooperative positioning. Besides, the sensor location is further refined in comparison with the estimatê s f . Similar to the first stage, an alternative ICWLS algorithm is proposed, which is also asymptotically efficient under Gaussian noise assumption.

1) OPTIMIZATION MODEL FOR MULTIPLE-TARGET COOPERATIVE LOCALIZATION
It is easy to see that the functional forms ofĀ(·, ·) andb(·, ·) in (38) are known, but the values ofr ands can not be obtained accurately because only their noisy values (i.e.,r andŝ f ) are available. To formulate the optimization model for localization, it is necessary to first introduce the following error vectorδ Applying a first-order Taylor series expansion ofb(r,ŝ f ) and A(r,ŝ f ) around the true valuesr ands yields It is clear from (101) thatδ is approximately Gaussian distributed with zero mean and the covariance matrix =C 1 (t,r,s)Ē A (C 1 (t,r,s)) T +C 2 (t,r,s) · MSE(ŝ f ) · (C 2 (t,r,s)) T (104) To delay the noise threshold before performance deviates suddenly from the CRB, the joint estimation ofū ands should be performed [35,49]. For this purpose, we need to introduce an augmented parameter vector as Then, from (98) and (104), the criterion function with respect tot can be expressed as Applying the WLS criterion, the weighting matrix˜ is given in (108), as shown at the bottom of the next page. Besides, it can be verified from (51) that vectort satisfies the constraints Combining (106) and (109) yields the following CWLS problem (110), as shown at the bottom of the next page. Similar to (74), (110) is also a nonconvex problem because the quadratic equality constraints are indefinite. Undoubtedly, the closed-form solution for (110) is analytically intractable; hence we have to solve (110) in an iterative way. In the sequel, an alternative efficient ICWLS algorithm is presented. On the other hand, it is worth noting that the weighting matrix˜ is related to the unknown vectort and so it should be updated recursively in the proposed procedure.

2) PROPOSED ITERATION ALGORITHM
The numerical technique adopted in Subsection VI.A is also applied here. First, we need to transform problem (110) into a convex programming. To this end, a set of linear equality constraints are employed instead of the homogenous and nonconvex constraints. Specifically, for each quadratic equality constraint in (110), if we replace one of the variablet with its last iteration's estimate, then the non-convex constraints become linear and convex. The resulted programming is a convex approximation of the CWLS problem (110). Mathematically, the approximate problem in the k + 1th iteration can be formulated [see (111)], as shown at the bottom of the next page. wheret(k) is the estimation result in the kth iteration. Similar to problem (75), the closed-form solution for (111) also exists and it can be obtained by means of matrix SVD. We are now ready to show the corresponding result.
Proposition 7: Lett(k) denote the vector composed of the first 8D elements oft(k). Defining the matrix and performing matrix SVD on G(t(k)) lead to where [Q 1 (k) Q 2 (k)] is an orthogonal matrix; R 1 (k) is an orthogonal matrix; 1 (k) is a diagonal matrix, whose diagonal entries are the singular values of G(t(k)). Besides, let us defineQ 2 (k) = blkdiag[Q 2 (k) I 6M ×6M ]. Then, the optimal solution for (111) is given bŷ The proof of Proposition 7 is similar to that of Proposition 2, so it is omitted due to limited space. Although the solution to problem (111) is an approximation of the optimal solution of the CWLS problem (110), it incorporates 2D approximate linear equality constraints to improve the localization accuracy. Intuitively, the optimal solution (114) can be directly used as the k + 1th iteration's value. But, to avoid the resulted solution sequence switching between two points, we adopt the update strategy in [29] once again. Specifically, the current iteration's value is determined using the linear combination of the solution to problem (111) and the estimation result in the last iteration. Furthermore, the weighting matrix˜ should also be updated at each iteration step because it is related tot.
The second stage of the proposed ICWLS method is formally outlined in Table 4.
At this point, we make three important remarks about the proposed algorithm.
Remark 11: Here, it is also very important to find a proper initial guess for the iterative algorithm. Notice that the unknown vectort is formed byt ands; hence we must choose good starting points for both of them. First, the estimation result in the first stage (i.e.,ŝ f ) can be directly considered as the initial value ofs. Second, the common least squares (LS) criterion can be applied to determine the starting point oft, which can written as t(0) = ((Ā(r,ŝ f )) TĀ (r,ŝ f )) −1 (Ā(r,ŝ f )) Tb (r,ŝ f ) (115) Then, the initial estimate oft is given byt(0) = [(t(0)) TŝT f ] T . Extensive Monte Carlo simulation tests show that this initial value can yield satisfactory localization accuracy.
Remark 12: As shown in Step 7, the updated value in the k + 1th iteration is the linear combination of the optimal solution (114) and the estimation result in the kth iteration. A large number of simulation results show that the localization performance is not sensitive to the numerical value of w 1 . In our simulation, w 1 and w 2 are both set to 0.5.
Remark 13: In optimization problem (111), the weighting matrix˜ is not accurately known because it relies ont,s and MSE(ŝ f ). It can be seen from Proposition 6 that matrix MSE(ŝ f ) also depends ons. Therefore, it is necessary to update˜ at each iteration step using the current iteration's value. Besides, matrix MSE(ŝ f ) can be obtained from (97). Theoretical analysis demonstrates that such an approximation of the weighting matrix does not affect the asymptotic properties of the estimator.

D. DISCUSSION ON THE PROPOSED ICWLS ESTIMATOR
As stated above, the proposed ICWLS method comprises two stages: a first stage for refining the sensor locations, and a second stage for multiple-target cooperative localization. At each stage, an efficient iterative algorithm is proposed. Moreover, it can be seen that the two iterative algorithms are developed based on a unified theoretical framework.
It should be pointed out that although the proposed estimator requires iteration, the convergence rate is fast, and moreover, the closed-form solution is available at every iteration. From our simulation results, it can be found that ten iterations are enough to achieve the convergence criterion. In addition, since the weighting matrices are updated recursively, more accurate weighting matrices are obtained, leading to a higher noise threshold than some existing methods before performance breaks away from the CRB.

VII. PERFORMANCE ANALYSIS OF THE PROPOSED ICWLS ESTIMATOR
This section is devoted to the performance analysis of the proposed ICWLS estimator. The theoretical derivation is mainly conducted for the iterative algorithm in stage-2, which provides the final estimation results. Similar to Subsection VI.B, the work also consists of two parts: (1) analysis on convergence property; (2) derivation for estimation MSE. Since the convergence analysis presented here follows the same theoretical framework as in Subsection VI.B, we just describe the main results due to limited space. Besides, the estimation MSE is deduced using the perturbation approach. More importantly, the MSE is proved to asymptotically attain the associated CRB given in Subsection III.A.

A. CONVERGENCE ANALYSIS FOR THE PROPOSED ICWLS ESTIMATOR
Here, the convergence analysis is similar to that in Subsection VI.B. Therefore, we only state the main results without detailed proofs. The objective is to illustrate that as long as the iterative algorithm converges, it must converge to the optimal solution for the CWLS problem (110).
According to the optimality of constrained optimization theory, we first get the following two Lemmas.
The proofs of Propositions 8-10 are similar to those of Propositions 3-5, respectively. Proposition 10 shows that if the proposed ICWLS algorithm converges, then it will provide the optimal solution for (110). There is no general approach, however, to ensure that the iterative algorithm must converge. Fortunately, the simulation results in Section VIII reveal that the convergence criterion can be reached within ten iterations in most cases.

B. ESTIMATION MSE FOR THE PROPOSED ICWLS ESTIMATOR
This subsection is devoted to the derivation of the MSE expression for the proposed ICWLS estimator. Moreover, the obtained theoretical MSE is analytically shown to equal the CRB given by (21). The first-order perturbation approach is adopted for this purpose.
Suppose that the proposed ICWLS algorithm in the second stage converges tot s , which is also the optimal solution for (110) as shown in Proposition 10. Since the final location estimates for target sources and sensors are obtained by û ŝ s s = t s , we should first deduce the theoretical MSE oft s . Let the estimation error int s be denoted as t s . In the theoretical frame of the first-order perturbation analysis, t s can be approximately expressed as a linear function with respect to the measurement noiseε as well as the estimation error s f in the first stage. Indeed, by taking the first-order approximation, t s is the optimal solution for the following constrained minimization problem [28]: The matrix SVD ofG(t) can be described as Due to the fact that t s is orthogonal to the range space of G(t), it is easy to see that t s belongs to the range space of ]. Therefore, we can express t s as t s =Q 2 (Q  Since the mean value of t s is equal to zero,t s is an approximately unbiased estimate oft. From (108) and (122), the estimation MSE oft s is given by whereP(r,s) = (Ã(r,s)) T˜ −1Ã (r,s). Three remarks are drawn in the sequel. Note that the number of columns ofH 1 (ū,s) is equal to that ofQ 2 , so it can be verified that range{H 1 (ū,s)} = range{Q 2 }.
It is noteworthy that (124) plays an important role in the proof of Proposition 11 described below.
To proceed further, the theoretical MSE of the joint estimate û ŝ s s is presented. Let its estimation error be ū s s s . Then, we can get which, combined with (123), yields The proof of Proposition 11 is described in Appendix F. The result in Proposition 11 demonstrates that the novel estimator is asymptotically efficient. In Section VIII, some simulations are conducted to support the theoretical development as well as the advantages of the newly proposed method.

VIII. NUMERICAL EXPERIMENT AND RESULT ANALYSIS
In this section, we present a set of Monte Carlo simulations in order to verify the analytical results and evaluate the accuracy of the proposed localization method. The estimation performance is assessed in terms of root-mean-square-error (RMSE) and radius of circular error probability (CEP) [50]. The RMSEs are obtained from the average of 5000 independent runs.
In the first set of experiments, the radius of CEP of the proposed method is compared with that of the TSWLS method in [34], which can locate multiple targets simultaneously but does not utilize the UAV calibration emitters. From this comparison, the performance improvement resulted from the use of the calibration emitters can be clearly observed. Consider a 3-D localization scenario consisting of 8 moving sensors whose true positions and velocities are tabulated in Table 5 Figs.2(a)-(d) depict the scatter plots of the target location estimates obtained from the proposed method as well as the TSWLS method in [34] through 2000 independent experiments. The radiuses of CEP for the two localization methods are also provided in these figures. Additionally, the position iterative sequence of the proposed method is plotted in Figs.3(a)-(d), where we perform 10 independent Monte-Carlo runs, each of which has different noise values. Meanwhile, the globally optimal solutions for target positions are also displayed. Note that the optimal solutions can be obtained by a direct grid search and they are generally not equal to the true values due to the presence of measurement errors.
From Figs.2 and 3 we can draw the following conclusions: (1) From the starting points described in Remarks 5 and 11, the new method is able to find the globally optimum solutions. Moreover, the convergence criterion can be reached within 10 iterations in most cases. Therefore, the validity of the convergence analysis presented in Subsections VI.B and VII.A is verified. (2) The CEP radius of the proposed method is much smaller than that of the TSWLS method in [34]. This finding reveals that the location accuracy could be significantly improved if the measurements from the calibration emitters are fully utilized.
In the second set of simulations, the newly developed method is compared with the TSWLS methods in [34], [35] and the Taylor-series iterative method extended from [51] in terms of estimation RMSE. Since these compared methods do not exploit the calibration emitters, we can accurately assess the performance gain due to the use of the calibration emitters from this comparison. On the other hand, as mentioned in [38], the differential calibration (DC) technique is commonly used in global positioning systems (GPS) to mitigate the effect of uncertainties in satellite position and velocity. This method can be easily extended to the localization scenario studied here. So, it is reasonable to compare our method with the DC technique. The simulation settings are the same as those in the first set of experiments, except that the covariance matrices (1) The proposed method is able to reach the corresponding CRB accuracy under a moderate noise level. This is because the new estimator is established based on the WLS criterion, which yields an asymptotically efficient solution under Gaussian model. Note that we prove analytically using first order analysis that the estimation performances in both stage-1 and stage-2 asymptotically attain the relevant CRBs. Hence, the simulation results confirm the theoretical development in this work. (2) The advantages of the new method over the TSWLS methods in [34], [35] and the Taylor-series iterative method extended from [51] are remarkable. The reason is that none of the compared methods exploits the measurements from the calibration signals. In other words, this significant performance gap mainly comes from the use of the calibration emitters. In addition, both the TSWLS methods and the Taylor-series iterative method can achieve CRB without calibration signal, before the thresholding effect starts to occur.
(3) Our method is superior to the DC technique and, moreover, the RMSE improvement increases as σ 1 increases. Besides, the RMSE of the DC method is larger than the associated CRB because its cost function is not asymptotically efficient. This observation is consistent with the analytical result in [38]. More importantly, the DC approach can not further refine the sensor position and velocity, while the proposed method can.
The third set of simulations is carried out to demonstrate that cooperation positioning for multiple targets can produce considerable performance gain compared to decoupled localization. For this purpose, the localization method presented in [43] is implemented for comparison. Since this method is developed for single-target localization scenario where the TDOA/FDOA measurements from both target source and calibration emitters are used, its performance can be chosen as benchmark to show performance improvement due to joint localization. The simulation scenario has an array of   (1) The performance improvement in positioning accuracy resulted from cooperative localization is noticeable. (2) The RMSE gap is more pronounced as σ 2 increases. This is because when σ 2 is large, the uncertainties in sensor position and velocity dominate the performance and then the advantage of cooperative positioning over decoupled localization becomes significant. (3) The estimation accuracy of the method in [43] is able to attain CRB for non-cooperative localization at moderate noise and error level. This observation confirms the analytical result in [43].
The fourth set of experiments studies the effect of the number of UAV calibration emitters. The simulation parameters are the same as those used to produce Fig.6. We examine the estimation accuracy of the new method for three cases.
In the first case, single calibration emitter is used for target localization and it is located at w 1  (Ā(r,s)) T (C 1 (t,r,s)) −TĒ −1 A (C 1 (t,r,s)) −1Ā (r,s) −(Ā(r,s)) T (C 1 (t,r,s)) −TĒ −1 A (C 1 (t,r,s)) −1C 2 (t,r,s) −(C 2 (t,r,s)) T (C 1 (t,r,s)) −TĒ −1 A (C 1 (t,r,s)) −1Ā (r,s) (MSE(ŝ f )) −1 + (C 2 (t,r,s)) T (C 1 (t,r,s)) −T ×Ē −1 A (C 1 (t,r,s)) −1C 2 (t,r,s) estimation and sensor velocity estimation as a function of σ 2 , respectively. From Figs.7 we can draw the following conclusions: (1) The proposed estimator is shown to yield the solution reaching the CRB accuracy under moderate noise level. This finding can further support the theoretical development in Section VII. (2) As expected, the localization accuracy will improve when more calibration signals are available for target localization. Furthermore, the higher the noise level, the greater the performance gain in estimation accuracy resulted from the increase in the number of calibration emitters. However, no matter how many calibration emitters exist, the positioning accuracy will not be higher than the bound for the case of no sensor location errors.

IX. CONCLUSION
This work focuses on the use of UAV calibration emitters for TDOA/FDOA multiple-target positioning in a wireless connected network. The study starts with deriving the CRB for TDOA/FDOA-based target location estimate when some UAV calibration emitters with known locations are available. The insight gained from the CRB reveals that the calibration signals are rather useful in reducing the effects of the uncertainties in sensor position and velocity. Moreover, the cooperative gain is significant when the multiple moving targets are located simultaneously in the presence of calibration emitters, even if the target sources are disjoint. The disjointness may be in time, frequency or both. The paper then proceeds to present a novel ICWLS estimator for multiple-target joint localization based on TDOA/FDOA measurements from both target sources and calibration emitters. The newly proposed method is composed of two stages. Specifically, the first stage refines the sensor locations by using the calibration TDOA/FDOA measurements as well as the statistical characteristic of the noisy sensor locations. The second stage provides the estimate of multiple-target locations by combining the measurements of target signals and the estimated values in the first stage. An efficient ICWLS algorithm is designed at each stage. Both the two algorithms are implemented based on matrix SVD, which is able to find a closed-form solution and update the weighting matrix at every iteration. Performance analysis of the proposed ICWLS estimator is also conducted, which consists of two parts: (1) analysis on convergence property; (2) derivation for estimation MSE. The results demonstrate that the proposed method, if converges, can produce the optimal solution of the formulated nonconvex CWLS problem. Besides, its estimation MSE can achieve CRB with calibration emitters at moderate noise. Simulation results verify the validity of the analytical results in this paper and confirm the superiority of the proposed method over some existing localization methods. It is worthy to point out that the proposed estimator can be extended to reduce the effects of the clock offsets and frequency deviations on target localization accuracy. In future work, we will focus on this topic.
First, it is readily checked that For convenience, we define the matrix and partition it into blocks (submatrices) as follows: Then, using the second equality in (A.1) yields Further, combining matrix identity (I) in Table 2, the first equality in (A.1) and (22) produces (A.5), as shown at the bottom of the next page. where It can be verified from (A.6) that YX 2 X −1 3 = X −1 1 X 2 Z. Then, from (A.5) and matrix identity (I) in Table 2, we get whereũ r = [ū T 2ū T 3 · · ·ū T D ] T , which does not comprise the first location vector u 1 .
On the other hand, it readily follows from matrix identity (II) in Table 2 that According to the definition of orthogonal projection matrix, (A.10) can be rewritten as Substituting (A.11) into the first equality in (A.9) and using (A.4) lead to which implies CRB(ū 1 ) ≤ CRB i (ū 1 ). The derivation described above can be easily extend to prove CRB(ū d ) ≤ CRB i (ū d ) for 2 ≤ d ≤ D, as long as we exchange the order ofū 1 andū d in composite location vectorū. At this point, the proof is completed.

B. PROOF OF PROPOSITION 2
It can be easily checked from (80) that where the second equality uses the relation (Q 11 (k +1)x 1 (k +1) satisfies the equality constraint in (75). As a consequence, we can write the optimal solution of (75) aŝ whereη o (k +1) fulfils the linear equation (G (c) (η(k))) Tη o (k + 1) = O and, moreover, it should minimize the cost function in (75). We are now ready to prove that the vectorη o (k + 1) is equal to Q (c) 12 (k + 1)x 2 (k + 1). Sinceη o (k + 1) belongs to the null space of (G (c) (η(k))) T , it can be expressed aŝ Substituting (83) into (B.5) leads to (B.6), as shown at the bottom of the next page, which implies that the objective function reaches its minimum value at z 0 = x 2 (k + 1). Then, combining (B.3) and (B.4) we havê At this point, the proof is ended.

C. PROOF OF PROPOSITION 4
It can be seen from (76), (78) and (90) that (C.2), as shown at the bottom of the next page.