A Lagrangian Multiplier Method for TDOA and FDOA Positioning of Multiple Disjoint Sources with Distance and Velocity Correlation Constraints

-is paper considers the source localization problem using time differences of arrival (TDOA) and frequency differences of arrival (FDOA) for multiple disjoint sources moving together with constraints on their distances and velocity correlation. To make full use of the synergistic improvement of multiple source localization, the constraints on all sources are combined together to obtain the optimal result. Unlike the existing methods that can achieve the normal Cramér-Rao lower bound (CRLB), our object is to further improve the accuracy of the estimation with constraints. On the basis of maximum likelihood criteria, a Lagrangian estimator is developed to solve the constrained optimization problem by iterative algorithm. Specifically, by transforming the inequality constraints into exponential functions, Lagrangian multipliers can be used to determine the source locations via Newton’s method. In addition, the constrained CRLB for source localization with distance and velocity correlation constraints is also derived.-e estimated accuracy of the source positions and velocities is shown to achieve the constrained CRLB. Simulations are included to confirm the advantages of the proposed method over the existing methods.


Introduction
Passive source localization is a classic problem that has been investigated over recent decades. It has received significant attention in various areas of signal processing research including unmanned aerial vehicle (UAV) [1,2], passive radar [3][4][5], microseismic sources [6][7][8], and underwater acoustics [9][10][11]. Typical measurements including the angle of arrival (AOA) [12,13], time of arrival (TOA) [14,15], time difference of arrival (TDOA) [15][16][17], frequency of arrival (FOA) [18,19], frequency difference of arrival (FDOA) [20,21], received signal strength (RSS) [22,23], and gain ratios of arrival (GROA) [24,25] (or some combination of them) can be used to estimate the source location. To the best of our knowledge, TDOA positioning is one of the most widely used schemes. Additionally, when there is a relative motion between the source and the sensors, FDOA can be exploited to further improve the estimation accuracy of the source position and provide an estimate of the source velocity. us, FDOA is often used in addition to TDOA for moving source localization.
In this paper, we discuss a constrained estimator for multiple moving source localization based on TDOA and FDOA measurements. e combination of TDOA and FDOA is an appropriate solution for locating the moving sources considered in this study. To improve the accuracy of the satellite systems only using TDOA measurements, Ho and Chan [26] first proposed a closed-form solution with TDOA and FDOA measurements. By introducing nuisance parameters, the well-known two-step weighted least-squares (TS-WLS) method was derived for positioning moving source [27]. is explicit algebraic solution can be used to estimate the source position and velocity with small computational burden. Although the estimate results of TS-WLS method are shown to attain the Cramér-Rao lower bound (CRLB) at low noise levels, the classical TS-WLS method suffers from the "threshold effect", which means the estimate accuracy decreases dramatically at moderate or higher noise level. To eliminate the threshold effects, many methods were proposed to overcome the drawbacks of the original TS-WLS method. Inspired by the multidimensional scaling (MDS) technique, a closed-form solution was proposed based on the optimization of the cost function to the scalar product matrix [28]. Utilizing Newton's method, a constrained WLS method was developed, which exploits the known relation between the intermediate variable and the source location coordinates explicitly. Subsequently, a convex relaxation techniques was utilized in TDOA and FDOA localization [29][30][31]. e method in [29] formulated the semidefinite programming (SDP) using the maximum likelihood (ML) criteria, and it can approach the CRLB better than WLS estimator. Besides, the semidefinite relaxation (SDR) was also applied to reformulated WLS problem to obtain an SDP [30], which made it applicable in real-time applications even under large measurement noise conditions. Different from the formulation methods in [29,30], a SDR-based method was developed based on the robust LS criterion, and it can determine the source location without approximation and initial estimate [31].
Besides the abovementioned closed-form methods using TDOA and FDOA measurements, some iterative algorithms can also be used for solving this type of problem. In [32], an improved algorithm was developed based on the original Taylor-series method [33]. By using Newton's method, a constrained weighted least-squares (CWLS) estimator derived from the WLS function is proposed in [34]. To reduce the computational cost in [34], a bi-iterative algorithm alternately calculates the position and velocity, which requires much fewer floating-point operations than the CWLS method [35,36]. Recursively relaxing the quadratic constraints into linear constraints, an iterative algorithm was developed to solve the CWLS optimization [37], namely, ICWLS method. In addition, the ICWLS method has a significant improvement in both accuracy and robustness. Meanwhile, Aziz [38] introduced the cuckoo search technique, which can avoid the convergence problem and outperform the previous methods [27]. In order to avoid expensive computation of grid searching, an importance sampling algorithm [39] was developed based on theorem of Pincus, and it can guarantee the global convergence. Based on A * search algorithm, a set of practical localization methods were developed to estimate the location of the microseismic source/acoustic emission (MS/AE) in complex structure [6][7][8].
To the best of my knowledge, most of the methods are shown to achieve the CRLB at low noise levels, and the negative influence of the threshold effects appear sooner or later. However, when the multiple disjoint sources are moving together, such as a group of UAVs or aircraft formation, certain distance and velocity constraints must exist among sources in the group, which may be possible to further enhance the accuracy of the source location. e constraints of distance and velocity correlation are obtained from prior information among the formation sources. Based on the information theory, the accuracy of estimation can be expected to improve when some prior information is available to the estimator [40].
In this study, we consider the localization problem of moving source with three types of constraints on their positions and velocities. Different from the previous methods, the ML technique is used to derive an unbiased estimator without approximation. e comparison between the CRLB and constrained CRLB indicates that the accuracy of the estimation can be improved by introducing the inequality constraints among the sources into the estimator. On the basis of Lagrangian multiplier technique, an efficient augmented cost function is established based on the exponential-based Lagrangian function, and the inequalities specifying the distance and velocity constraints are transformed into a series of exponential penalty terms in the cost function. An iterative technique based on Newton's method is employed to determine the location of multiple sources, because it is quite difficult to develop a closed-form solution for the nonconvex and highly nonlinear optimization problem without linear approximation. Combining all the whole constraints and estimating sources together takes the advantage of this synergistic improvement. Furthermore, extensive simulations illustrated the good performance of the proposed method and verified its efficiency in attaining the constrained CRLB. e rest of this paper is organized as follows. In Section 2, we formulate the localization problem for multiple sources that are moving together and introduce the constraints on source positions and velocities. e constrained CRLB is derived in Section 3. Section 4 provides the expandable estimator and the iteration algorithm. In Section 5, the simulations, which contain near-field and far-field source cases, are conducted. Finally, the conclusions are given in Section 6. e details of complicated matrices can be seen in Appendix. e main notations used in this paper are listed in Table 1 for convenience.

Localization Scenario and Problem Formulation
en, the TDOA measurements between sensor pair m and 1 is where f (p) m (u n , s) is a function with respect to u n and s, ε (p) m1,n denotes the TDOA noise. Although velocity _ u n is useless here, u n is still used rather than u n in the function for unification. In this way, the composite function of the whole sensor pairs can be written as a form of source n: where r (p) When there is a relative motion between the source and the sensors, FDOA can be exploited not only to further improve the estimation accuracy of the source position, but also to provide an estimate of the source velocity. As a result, FDOA is often used in addition to TDOA for moving source localization. After taking the time derivative of (2), the FDOA measurements between sensor pair m and 1 can be expressed as where ε (v) m1,n is the FDOA noise; f (v) m (u n , s) is a function with respect to u n and s. Like the expression of TDOA, the composite vector of the whole sensor pairs from source n can be written as where Combining TDOA and FDOA of source n together, the measurements vector is given: where in which ε n is the TDOA and FDOA noise vector of source n. Measurement noise is assumed to be a zero-mean Gaussian random vector with covariance matrix: n , and the cross-correlation between the TDOA and FDOA noise can be characterized: en putting the all measurements together, TDOA and FDOA measurements from N sources yields where e composite vector ε follows zero-mean Gaussian distribution with covariance matrix: 2.2. e Inequality Constraints on Source Location. is paper considers the scenario of organized sources, such as aircraft formation and tactical squads. ese sources can be located in a certain space and move at approximately equal velocities. To stay information, each source in the group has to transmit and receive signals from the others unless they keep radio silence. Hence, the probability of detection increases. If the prior information is used, information theory tells us the accuracy of estimation must be improved [40][41][42]. In this case, we consider a series of distance constraints and velocity constraints for the multiple sources.
To improve the accuracy of estimation, the following three kinds of constraints are given. Figure 1 shows one of the source location scenarios considered in this paper. Identify the formation of a set of aircraft is a typical multiple source localization problem. To establish a constrained estimator, we first transform the mathematical relationship into inequality constraints.
In practice, prior range constraints must exist for the multiple sources that accompany each other. As is shown in Figure 2, when we have N sources in a formation, there are at least N(N − 1)/2 distance constraints among them. Using prior information such as the safe distance, the two sources n 1 and n 2 can be limited to a range d n 1 ,n 2 that is greater than or equal to the true distance. A normal distance constraint between source n 1 and n 2 can be described as To represent the position and velocity of source n with vector u, the matrix J (p) is defined to select the parameters. In addition, the position of source n can be written as J (p) n u. e distance constraint can be rewritten as e FDOA measurements depend on the relative movement between sensors and sources. To maintain communications or ensure safety, certain rules must be specified to control the velocity for the sources moving in formation. erefore, it is much more reasonable to define the velocity correlation using the velocity constraints when describing the relationship between the two sources shown in Figure 3, which is given by  e velocity of source n can be written as . In contrast to the distance constraints, the true velocity correlation should be larger than ρ n 1 ,n 2 in the formulas. However, the velocity correlation cannot distinguish velocities that are completely opposite or doubled and redoubled. For instance, in the example shown in Figure 4, all the vectors have the highest possible correlation, but the directions or absolute values can be very different from those of the others. Hence, correlation is not sufficient for describing the similarity of the velocities.
To avoid the fuzzy estimation caused by opposite directions and describe the velocity accurately, we consider that the difference between two velocities in one dimension should be less than some threshold, so the constraint is given Using the matrix i (1) 3 , the auxiliary constraint can be rewritten as We choose the X-direction component to be the limiting factor. Obviously, the choice will not make any difference. Yand Z-direction component can also be set as the limiting factors. e formulas are given as follows: e constraint for the X-component can be straightforwardly extended to give constraints in the Y-and Zdirections. If all three directions are considered at the same time, another two limiting factors will be incorporated in the estimator. e constraints in the X-direction are sufficient to determine the velocity correlation.

Derivation for the Constrained CRLB
As is well-known, the multiple-parameter CRLB is widely used for describing estimator performance [43], especially in multidimensional parameter estimation. e CRLB is the lowest possible variance that an unbiased estimator can achieve. When there is no prior constraint, the CRLB of the error covariance can be described as the inverse of the Fisher information matrix, which is defined as e details of zr/zu T are shown in Appendix A. Each element in observation vector r is governed by the Gaussian probability distribution, and vector u lies in the unconstrained parameter space Θ � R 6N×1 . If the parameter is restricted by the constraints, the unconstrained parameter space is reduced to a subset of Θ called Θ C .
Note that in general constrained optimization, there are both equality constraints and inequality constraints. As in definition in [40], the constraints that appear in the optimization decompose the constraint set Θ C into two sets: feasible points, which are all interior points or on the boundary of the inequality constraints, and infeasible points. In the method proposed in this paper, we only use pure inequality constraints, which are continuously differentiable, as prior information to further improve the estimated accuracy. In addition, an equality constraint can be transformed into two inequality constraints which have sign reversal. e estimator can be straightly expended to the mix of inequality and equality constraints. Hence, the solutions of the constrained positions and velocities are the feasible solutions in the closure of the interior of Θ C . Let parameter space Θ C be defined by all the constraints mentioned in Section 2.2. To develop a unified form to establish the vector function, the three inequality constraints for source pair n 1 and n 2 can be rewritten as functions where c i (i � 1, 2, 3) represents the three kind of constraints respectively. ere are N(N − 1)/2 functions about distance constraint overall. For example, the distance constraint of sensor n is given in the vector en, we have the function vector putting all the distance constraints together: Figure 3: Diagram for velocity constraints.

Mathematical Problems in Engineering
It is also true for the other constraints, which can be written as Finally, combining (24) and (25), the composite vector is given As is proved in eorem 1 [40], the (3N(N − 1)/2) × 6N gradient matrix ∇G u , whose rank is smaller than the number of constraints 3N(N − 1)/2, error covariance matrix Σ u satisfies where Matrix Q u is defined as where ∇G u is the gradient matrix of (27). e details of ∇G u are given in Appendix B. In (29), matrix Q u contains the prior information of multiple sources location. en, substituting (30) into (29), the expression of CRLB C (u) is obtained Matrix Q u clearly can reduce the original CRLB, because the second term of (31) is a semidefinite matrix. e constrained CRLB is called CRLB C (u) here.
Obviously, when the constraints are active, which means that u is both in feasible sets and satisfies the constraints, the solution of u is not the optimal result that we want. If the initial value is an infeasible point, the algorithm will keep searching for the feasible point until convergence. eorem 1 in [40] indicates that the bound can be further reduced with the constraints. In contrast, the theorem also proves that the constraints can only reduce the CRLB on the component error variances. Hence, the constrained CRLB is the lowest bound that the method can achieve with the corresponding constraints. e associated CRLB and constrained CRLB will be used to evaluate the performance of the proposed method and existing methods in Section 5.
In Figure 5, the CRLB and constrained CRLB for different numbers of sources are shown to validate the accuracy improvement. A network composed of several random moving sensors is used to locate multiple sources with constraints. As the number of sources increases, more constraints are incorporated. It is evident that the accuracy of the position and velocity estimates for the same source improves as the number of sources increases.

Constrained Estimator for Source Localization.
e ML estimator of the sources' position and velocity is to minimize ε T Q − 1 ε. With (22), we can establish the following constrained estimator: For all N sources in a group, there will be 3N(N − 1)/2 inequality constraints in total. Overall, the ML estimator is transformed into a constrained estimator. On the other hand, an equality constraint can be transformed into two inequality constraints with sign reversal.

Lagrangian Multiplier Technique for Constrained
Estimator. To solve a constrained optimization problem that contains equality constraints and inequality constraints, the usual approach is to transform the constrained optimization model into an unconstrained optimization model by applying a Lagrangian multiplier. Because of the nonlinearity of the constraint function for the source location in y . u 1 .  Mathematical Problems in Engineering (32), it is difficult to derive closed-form or near closed-form solutions using the results of previous studies. Here, the constraints are converted into penalty items and an iterative method using Lagrangian multipliers is proposed. e Lagrangian multiplier technique provides several different ways to establish an augmented Lagrangian function, including quadratic and exponential types. In contrast with the quadratic augmented Lagrangian function for inequality constraints, the exponential augmented Lagrangian function used in the proposed method is twice differentiable. Hence, the exponential function is adopted to transform the inequality constraints into exponential penalty terms in the cost function. Note that the cost function and constraint functions in (32) are twice differentiable with respect to u everywhere. We substitute the constraint functions into an exponential function: So exponential penalty functions are produced in which en, the exponential Lagrangian function is given where in which

Mathematical Problems in Engineering
Function L(u, λ, τ) is the new cost function. Overall, the constrained estimator is transformed into an unconstrained estimator, and we just need to compute u as u � arg min u∈R 6N×1 L(u, λ, τ) � arg min u∈R 6N×1 L(u).
Remark 1. Note that (33) has following properties: As is shown in (ii), when the parameter is on the boundary of constrained space, the penalty term is equal to zero, as expected. In addition, it is easier for the exponential function to take derivative.
We associate λ (c i ) n 1 ,n 2 > 0 as a multiplier and τ (c i ) n 1 ,n 2 as a positive penalty parameter with the constraint function g (c i ) n 1 ,n 2 (u). e unconstrained estimator can then be solved with an iterative method, and the multipliers are updated at the end of each minimization using the first equation in (42). After taking the derivative of the exponential function and substituting the constraint functions into the exponential function (33), we can obtain the second equation in (42): Remark 2. Note that for a fixed λ (c i ) n 1 ,n 2 > 0, as τ n 1 ,n 2 tends to ∞ for all infeasible u and to zero for all feasible u; otherwise, for a fixed τ (c i ) n 1 ,n 2 , as λ (c i ) n 1 ,n 2 ⟶ 0 the corresponding penalty term goes to zero whether u is feasible or not. e exponential Lagrangian cost function contains the original cost function (32) and the penalty terms.

Remark 3.
e exponential Lagrangian function can only converge to the minimum value in the feasible sets if the Lagrangian multipliers converge to a limit. Otherwise, the penalty terms will block the convergence process. e constrained estimator is indeed an NP-hard problem. We can only guarantee that the proposed method converges to the optimal solution in the feasible sets. As expected, the numerical results in Section 5 indicate that the proposed method converges.
As is commonly known, the quadratic augmented Lagrangian function for inequality constraints has a theoretical advantage for constrained minimization problems. To convert a constrained problem into an unconstrained problem, the inequality constraint functions must be transformed into associated equalities by introducing auxiliary variables. Hence, the iterative method for searching for the optimal result needs to iterate both on the multipliers and auxiliary variables. Simultaneously, the quadratic terms are allowed to employ the augmented Lagrangian function only if the iteration is not within the constrained parameter space.
us, the augmented Lagrangian function is no longer differentiable, and an extra judgment on the falling direction must be performed for each iteration. For multiple parameters, the complexity of the program and huge amount of computation required are impractical. e principal motivation for using the exponential function is that, in contrast to the usual quadratic augmented Lagrangian function for inequality constraints, the cost function in (37) is twice differentiable if ε T Q − 1 ε and the constraint functions are also twice differentiable. In this way, the corresponding unconstrained minimization can be solved by using Newton's method with guaranteed superlinear convergence without determining if the falling direction is right for each iteration. Once the exponential Lagrangian function is established, the constrained optimization can be converted to an unconstrained estimator and Newton-type methods can be used to iterate the optimization directly.

Newton's Iteration Algorithm of the Lagrangian Estimator.
For the above optimization, the closed-form solution is difficult to derive. Hence, we introduce Newton's iteration, which is widely used in optimization. Differentiating L(u) with respect to u, gradient vector is obtained: where Hessian matrix of L(u) is given by 8 Mathematical Problems in Engineering e details of the gradient and Hessian matrix can be seen in Appendixes A and B. e process of proposed algorithm can be seen in Table 2.
Remark 4. we set the rough estimation given by grid search method as the initial position and velocity. Although the estimation is not asymptotically optimal, it is a sufficiently precise initial guess for our algorithm to converge. e simulations in Section 5 prove that this initial guess can provide good results.

Simulation Result and Analysis of Proposed Estimator
To examine the performance of the proposed method, we apply it to various scenarios and compare the results with those given by three similar methods, including TS-WLS method in [27], Gauss-Newton method in [33], and ICWLS method in [44]. As described in Section 2, zero-mean Gaussian noise is assumed to disturb the TDOA and FDOA measurements. TDOA and FDOA estimates were generated by adding to the values zero-mean Gaussian noise. For one source, the covariance matrices of the TDOA and FDOA noise are σ 2 d R and σ 2 f R, respectively. e noise power σ 2 d and σ 2 f can be obtained from the results of the TDOA and FDOA estimation before the location estimation [45,46]. e covariance matrices of TDOA and FDOA were σ 2 d R and σ 2 f R, which is one of the typical values of the TDOA and FDOA estimation in practice.
e matrix R is defined as a (M − 1) × (M − 1) matrix, the diagonal elements of which are equal to unity and 0.5 otherwise.
We conduct 5000 Monte Carlo experiments for all the following simulations. To evaluate the localization accuracy of the proposed method and existing methods, the root mean-square error (RMSE) which is equal to the variance of the estimate results of Monte Carlo experiments is defined: Usually, localization scenarios are classified as either near-field or far-field scenarios. In near-field scenarios, the distance between the source and the sensors is greater than the distance between pairs of sensor. In far-field scenarios, the distance between the source and the sensors is less than the distance between sensor pairs. By examining both nearand far-field scenarios, we perform a better analysis and obtain more believable conclusions.

Experiment 1: Near-Field Scenario.
In this subsection, the proposed method is compared with the performance in near-field scenario. e positions and velocities of the sensors are listed in Table 3. e locations of the sensors and unknown sources can be seen in Figure 6 in which the arrows represent their velocity and direction.
We assume that there are two moving sources whose positions and velocities are u 1 � 654 803 708 T m, _ u 1 � 20. 23 15.18 40.34 T m/s and u 2 � 652 808 710 T m, _ u 2 � 20 15 38 T m/s, respectively. e sensor network is employed to locate the two sources. In this case, we set σ 2 f � 0.25σ 2 d for the covariance matrices of TDOA and FDOA. e distance between two sources in a group can be small. As for the velocity constraints, the velocity correlation of the two sources will have extremely high similarity. Each one will have nearly same the velocity as the others. If the initial value of Lagrangian multiplier in λ is too large, the method will fail to converge to the true solution. Otherwise, the penalty functions are invalid and the solution is not in Θ C . Both will lead to a clearly incorrect estimated result. e true distance between the two sources is about 6 m, the velocity correlation is 0.9998, and the difference in the sources velocities in the X-direction is 0.23 m/s. erefore, we set d 1,2 � 10 m, ρ 1,2 � 0.999, and Δ _ x 1,2 � 0.5 m/s for the three constraints. e first and third are approximately two times their true values. When the two sources have the same direction of motion and maintain a similar relative distance, the velocity correlation approaches 1 asymptotically, so we set ρ 1,2 � 0.999. e Lagrangian multiplier vector λ is set to [2 1 1] T and penalty parameter vector τ � 10 10 10 T initially. Figures 7 and 8 show the accuracy of the two parameter estimates of the proposed method in terms of RMSE as the noise power increases. Comparing the proposed method to the existing methods, both with respect to position and velocity, all the methods can achieve the CRLB and the constrained CRLB, respectively.  [27]. e positions and velocities of the sensors are listed in Table 4. e locations of the sensors and unknown sources can be seen in Figure 9. τ , positive penalty parameter which are fixed in iteration. η, the maximum number of iterations.(η � 20 in this paper) μ (0) , initial step for iterations, σ attenuation factor. ξ, threshold, a sufficiently small and positive value (10 −6 in this paper). τ) , initialize the value for comparison and save gradient vector ∇L(u (0) , λ (0) , τ) and hessian matrix ∇ 2 L(u (0) , λ (0) , τ).

End
Output: u (final) as the optimal solution    Here, the covariance matrix σ 2 d R of TDOA is assumed to be ten times σ 2 f R of FDOA. e distance among the three far-field sources is less than 20 m. e velocity correlation and the X-direction difference in velocity are the same as in the previous simulation. erefore, we set d n 1 ,n 2 � 25 m, ρ n 1 ,n 2 � 0.999, and Δ _ x 1,2 � 1 m/s for the three constraints. e Lagrangian multiplier λ is set to 21 1×3 1 1×3 1 1×3 T and penalty parameter vector τ � 101 1×3 101 1×3 101 1×3 T which are the same as those in near-field scenarios. Figures 10-12 depict the results. e results of the three-source localization are similar to that of two sources. e estimate result of the Gauss-Newton method decreases dramatically when the noise power is larger than a certain value, usually 0 dB. Because the second-order error terms in the solution are ignored, the threshold effect of the Gauss-Newton method appears at a lower noise power level than it does for the other methods. It is distressing that the ICWLS method had several patchy estimates. Both of the existing methods can usually achieve the associated CRLB at high SNR.
e proposed method has a higher threshold of noise than those existing methods and can always the constrained CRLB as shown in the figures.
As shown in the statistical results, the proposed method can fail to achieve the constrained CRLB in low noise levels, because the initial guess could easily meet the constraints when the disturbance is weak. In this way, the constraints are inactive. Hence, a lower noise power means that the constraints have less influence in the convergence processing. As the noise power increases, more initial guesses are feasible points. e penalty function forces the iterative parameters to move into the constrained parameter space and helps to convergence to the optimal location.
Overall, the constrained CRLB is the lowest bound that the method can ultimately achieve. e derivation of the      constrained CRLB (31) shows that the lower bound is not related to the value of the constraints. After taking the derivation in the both side of the constraint formulations, the constants of the prior information for distance and velocity constraints cannot be found in the constrained CRLB anymore. erefore, the prior distance and velocity correlation do not change the constrained CRLB, which just concerns the mathematical expression of the inequalities. Whether the method can achieve the lowest bound is dependent on the relationship between the true location and the prior information in the constraints. When the distance constraints and velocity correlation are strictly satisfied by the true location, the method is valid and can achieve the constrained CRLB.

Comparison for Changing Positions.
It is necessary to analyze the performance of the methods as the source positions and velocities change. When the geometric dilution of precision (GDOP) is large, it is difficult to attain the corresponding CRLB, especially in far-field scenarios. To examine the effect of poor location geometry, the multiple sources are located at a range of R � 2400 m with an azimuth of π/4 rad and are almost static relative to one another while moving along an arc from an elevation angle c of −2π/9 rad to 2π/9 rad with an angular velocity of 0.02 rad/s. e distance between the two sensors is 10 m from beginning to end, and we set d 1,2 � 20 m as the distance constraint. If the two sources are moving in different directions, they will move away from each other over time. Hence, we set ρ 1,2 � 0.999 for the velocity correlation, where the true value is one. e TDOA and FDOA noise powers are σ 2 d � 0.025/c 2 and σ 2 f � 0.0025/c 2 , respectively. e stationary sensor positions are listed in Table 5.
We conducted 5000 Monte Carlo experiments for every point in its trajectory. e Lagrangian multiplier vector λ was set to [2 1 1] T and penalty parameter vector was τ � [10 10 10] T at the beginning. e proposed method is compared with the existing methods.
As shown in Figure 13, the proposed method can generally achieve the constrained CRLB. As the elevation angle c increases, the improvement of the proposed method becomes more apparent. e TS-WLS method and Gauss-Newton method can achieve the CRLB in all directions, but the ICWLS method cannot estimate the source location effectively. e iterative processes of the Lagrangian multipliers and the Euclidean norm of the gradient vector at the starting point of Figure 14, c � −2π/9 are shown in Figure 15. It can be seen that each line represents a Monte Carlo experiment and various colors are used to differentiate the experiment results, and the colors do not have obvious physical meaning. e trend of the lines indicates that the proposed method can converge within 10 iterations, which shows the convergence is also good and fast.

Comparison for Changing
Velocities. Now, we consider a scenario that the velocity of the formation changes over time. An array of six fixed receivers was chosen, and their positions are listed in Table 6. e location of geometry is shown in Figure 16.
e true starting coordinates of the formation sources are u 1 � 0 0 0 T m, u 2 � 7 −10 3 T m, and u 3 � 2 10 1 T m. e sources were moving from a near-field area to a far-field area with the same velocity. eir initial velocities of them are all 20 T m/s and their acceleration is  Figures 17-19 illustrate the RMSEs of the proposed method and the existing method corresponding to the three sources, respectively. e accuracy of position and velocity estimate of the proposed method in terms of RMSE is shown as the time goes. e estimation accuracy of the proposed method reaches the constrained CRLB from the beginning to the end. e results from the TS-WLS method, Gauss-Newton method, and ICWLS method are also given. Simulation results indicate that the three existing methods fail to locate the sources at the beginning, because the TDOA and FDOA measurement are equal to zero when the sources are around the center of the network, [0 0 0] T in this case.
is is reflected in Figures 17-19, and the proposed method is still significantly superior to the other three methods. e formation constraints can effectively eliminate ambiguous estimation.

Experiment 4: Comparison of the Average CPU Run Time.
We compare the run time of the proposed method with that of the other three methods mentioned above. e simulations were all conducted using MATLAB on a computer with an Intel Core i7-9700K CPU @ 3.6 GHz and 16 GB RAM. e conditions of these experiments are the same as those in Sections 5.1 and 5.2, respectively. e results of average run time of one Monte Carlo experiment are listed in Table 7.
With a good initial guess, the iterative methods need more time to converge to the optimal result (or even fail) as the SNR increases. As a result, the iterative methods take longer than the closed-form or near closed-form methods. Compared with TS-WLS and ICWLS, the iterative process in the proposed method takes a relatively long time. e run time of the Gauss-Newton method is less than that of the proposed method because the          proposed method takes more steps to converge in the same scenarios. However, Table 7 indicates that the average CPU time is at the millisecond level on the MATLAB platform. With advances in processor technology and the development of programmability, the estimation operation will become less time-consuming.

Conclusion
is paper focused on multiple disjoint sources localization with constraints on distance and velocity correlation. To improve the accuracy using the constraints, an expandable estimator was developed and augmented cost function was established to solve the problem. In contrast to the usual use of Lagrangian multiplier, the constraints were transformed into an exponential function as a penalty term. In this way, we ensure that the cost function is still twice differentiable so that Newton's method may be used. e constrained CRLB was derived in this paper. is derivation shows that the bound has nothing to do with the value of the constraint range, so the constrained CRLB is the lowest bound that the method can achieve. Finally, comprehensive simulations were conducted in different scenarios. ree similar techniques (TS-WLS, ICWLS, and the Gauss-Newton method) were compared with the proposed method. e results show that the proposed method can achieve the constrained CRLB and has a higher threshold than the other approaches.
Data Availability e data can be obtained from the corresponding author of this paper.

Conflicts of Interest
e authors declare that they have no conflicts of interest.