Abstract
This 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. The 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.
1. 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–5], microseismic sources [6–8], and underwater acoustics [9–11]. Typical measurements including the angle of arrival (AOA) [12, 13], time of arrival (TOA) [14, 15], time difference of arrival (TDOA) [15–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. Thus, 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. The 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]. This 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–31]. The 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–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. The 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. The 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.
The 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. The 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. The details of complicated matrices can be seen in Appendix. The main notations used in this paper are listed in Table 1 for convenience.
2. Localization Scenario and Problem Formulation
2.1. Signal Model for Localization Based on TDOA and FDOA Measurements
In this paper, we consider a three-dimensional localization scenario, where moving sensors are used to determine the moving sources of interest with unknown positions and velocities . The three-dimensional scenario can be straightly transformed into two-dimensional scenario generally. The location of the th source can be expressed as and the composite source location vector can be written as . The sensor positions and velocities are assumed known perfectly, and we can collect the sensor positions and velocities as the vector . The composite sensor location vector can be written as . Only overdetermined scenario can be used to estimate the source location precisely, and the underdetermined scenario is not the orientation of this paper. Usually, at least four sensors are used for localization using TDOA and FDOA measurements.
Note that range difference of arrival (RDOA) measurements can be converted into TDOA measurements using signal propagation speed. Without loss of generality, sensor can be set as the reference sensor. Let be the true distance between the source and sensor :
Then, the TDOA measurements between sensor pair and iswhere is a function with respect to and , denotes the TDOA noise. Although velocity is useless here, is still used rather than in the function for unification. In this way, the composite function of the whole sensor pairs can be written as a form of source :where
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 and can be expressed aswhere is the FDOA noise; is a function with respect to and . Like the expression of TDOA, the composite vector of the whole sensor pairs from source can be written aswhere
Combining TDOA and FDOA of source together, the measurements vector is given:wherein which is the TDOA and FDOA noise vector of source . Measurement noise is assumed to be a zero-mean Gaussian random vector with covariance matrix:where , , and the cross-correlation between the TDOA and FDOA noise can be characterized:
Then putting the all measurements together, TDOA and FDOA measurements from sources yieldswhere
The composite vector follows zero-mean Gaussian distribution with covariance matrix:
2.2. The Inequality Constraints on Source Location
This paper considers the scenario of organized sources, such as aircraft formation and tactical squads. These 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–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 sources in a formation, there are at least distance constraints among them. Using prior information such as the safe distance, the two sources and can be limited to a range that is greater than or equal to the true distance. A normal distance constraint between source and can be described as
To represent the position and velocity of source with vector , the matrix is defined to select the parameters. In addition, the position of source can be written as . The distance constraint can be rewritten as
The 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. Therefore, 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
The velocity of source can be written as where the matrix .
In contrast to the distance constraints, the true velocity correlation should be larger than 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.
(a)
(b)
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 , 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. Y- and Z-direction component can also be set as the limiting factors. The formulas are given as follows:
The constraint for the X-component can be straightforwardly extended to give constraints in the Y- and Z-directions. If all three directions are considered at the same time, another two limiting factors will be incorporated in the estimator. The constraints in the X-direction are sufficient to determine the velocity correlation.
3. 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. The 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
The details of are shown in Appendix A. Each element in observation vector is governed by the Gaussian probability distribution, and vector lies in the unconstrained parameter space . If the parameter is restricted by the constraints, the unconstrained parameter space is reduced to a subset of called .
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 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. The 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 . Let parameter space 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 and can be rewritten as functionswhere represents the three kind of constraints respectively.
There are functions about distance constraint overall. For example, the distance constraint of sensor is given in the vector
Then, we have the function vector putting all the distance constraints together:
It is also true for the other constraints, which can be written aswhere
Finally, combining (24) and (25), the composite vector is given
As is proved in Theorem 1 [40], the gradient matrix , whose rank is smaller than the number of constraints , error covariance matrix satisfieswhere
Matrix is defined aswhere is the gradient matrix of (27). The details of are given in Appendix B. In (29), matrix contains the prior information of multiple sources location. Then, substituting (30) into (29), the expression of is obtained
Matrix clearly can reduce the original CRLB, because the second term of (31) is a semidefinite matrix. The constrained CRLB is called here.
Obviously, when the constraints are active, which means that is both in feasible sets and satisfies the constraints, the solution of 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. Theorem 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. The 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.
(a)
(b)
4. The Constrained Estimator Based on ML Estimator and Lagrangian Multiplier Technique with Newton’s Method
4.1. Constrained Estimator for Source Localization
The ML estimator of the sources’ position and velocity is to minimize . With (22), we can establish the following constrained estimator:
For all sources in a group, there will be 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.
4.2. 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 (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. The 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 everywhere. We substitute the constraint functions into an exponential function:
So exponential penalty functions are producedwherein which
Then, the exponential Lagrangian function is givenwherein which
Function is the new cost function. Overall, the constrained estimator is transformed into an unconstrained estimator, and we just need to compute as
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 as a multiplier and as a positive penalty parameter with the constraint function . The 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 , as , penalty term tends to for all infeasible and to zero for all feasible ; otherwise, for a fixed , as the corresponding penalty term goes to zero whether is feasible or not. The exponential Lagrangian cost function contains the original cost function (32) and the penalty terms.
Remark 3. The 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. The 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. Thus, 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.
The 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 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 super-linear 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.
4.3. 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 with respect to , gradient vector is obtained:where
Hessian matrix of is given by
The details of the gradient and Hessian matrix can be seen in Appendixes A and B. The 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. The simulations in Section 5 prove that this initial guess can provide good results.
5. 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 and , respectively. The noise power and can be obtained from the results of the TDOA and FDOA estimation before the location estimation [45, 46]. The covariance matrices of TDOA and FDOA were and , which is one of the typical values of the TDOA and FDOA estimation in practice. The matrix is defined as a matrix, the diagonal elements of which are equal to unity and 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 near- and far-field scenarios, we perform a better analysis and obtain more believable conclusions.
5.1. Experiment 1: Near-Field Scenario
In this subsection, the proposed method is compared with the performance in near-field scenario. The positions and velocities of the sensors are listed in Table 3. The 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 m, m/s and m, m/s, respectively. The sensor network is employed to locate the two sources. In this case, we set for the covariance matrices of TDOA and FDOA. The 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 . Both will lead to a clearly incorrect estimated result. The 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. Therefore, we set m, , and m/s for the three constraints. The 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 . The Lagrangian multiplier vector is set to and penalty parameter vector 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.
(a)
(b)
(a)
(b)
5.2. Experiment 2: Far-Field Scenario
In this case, a far-field scenario is considered, in which the three sources’ true positions are m, m, and m. Their velocities are m/s, m/s, and m/s, respectively. We use the same sensor location in [27]. The positions and velocities of the sensors are listed in Table 4. The locations of the sensors and unknown sources can be seen in Figure 9.
Here, the covariance matrix of TDOA is assumed to be ten times of FDOA. The distance among the three far-field sources is less than 20 m. The velocity correlation and the X-direction difference in velocity are the same as in the previous simulation. Therefore, we set m, , and m/s for the three constraints. The Lagrangian multiplier is set to and penalty parameter vector which are the same as those in near-field scenarios. Figures 10–12 depict the results.
(a)
(b)
(a)
(b)
(a)
(b)
The results of the three-source localization are similar to that of two sources. The 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. The 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. The 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. The 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. Therefore, 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.
5.3. Experiment 3: Changing Positions and Velocities
5.3.1. 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 m with an azimuth of rad and are almost static relative to one another while moving along an arc from an elevation angle of rad to rad with an angular velocity of rad/s. The distance between the two sensors is 10 m from beginning to end, and we set 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 for the velocity correlation, where the true value is one. The TDOA and FDOA noise powers are and , respectively. The stationary sensor positions are listed in Table 5.
We conducted 5000 Monte Carlo experiments for every point in its trajectory. The Lagrangian multiplier vector was set to and penalty parameter vector was at the beginning. The 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 increases, the improvement of the proposed method becomes more apparent. The TS-WLS method and Gauss–Newton method can achieve the CRLB in all directions, but the ICWLS method cannot estimate the source location effectively.
(a)
(b)
The iterative processes of the Lagrangian multipliers and the Euclidean norm of the gradient vector at the starting point of Figure 14, 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. The trend of the lines indicates that the proposed method can converge within 10 iterations, which shows the convergence is also good and fast.
(a)
(b)
(c)
(d)
5.3.2. 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. The location of geometry is shown in Figure 16. The true starting coordinates of the formation sources are m, m, and m. The sources were moving from a near-field area to a far-field area with the same velocity. Their initial velocities of them are all m/s and their acceleration is m/s2. During the movement, they remain relatively static. The TDOA noise variance was and FDOA noise was in this case. The Lagrangian multiplier is and penalty parameter is . The conditions of the constraints are set as m, , and m/s. The simulation results are obtained via 5000 Monte Carlo experiments.
Figures 17–19 illustrate the RMSEs of the proposed method and the existing method corresponding to the three sources, respectively. The accuracy of position and velocity estimate of the proposed method in terms of RMSE is shown as the time goes. The estimation accuracy of the proposed method reaches the constrained CRLB from the beginning to the end. The 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, in this case. This is reflected in Figures 17–19, and the proposed method is still significantly superior to the other three methods. The formation constraints can effectively eliminate ambiguous estimation.
(a)
(b)
(a)
(b)
(a)
(b)
5.4. 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. The simulations were all conducted using MATLAB on a computer with an Intel Core i7-9700K CPU @ 3.6 GHz and 16 GB RAM. The conditions of these experiments are the same as those in Sections 5.1 and 5.2, respectively. The 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. The 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.
6. Conclusion
This 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. The constrained CRLB was derived in this paper. This 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. Three similar techniques (TS-WLS, ICWLS, and the Gauss–Newton method) were compared with the proposed method. The results show that the proposed method can achieve the constrained CRLB and has a higher threshold than the other approaches.
Future work will focus on developing more robust algorithms with erroneous sensor positions and physical obstacles.
Appendix
A
Let , then, in (21) can be expressed aswhere
Sub-matrices in (A.2) is given abovewhere
B
Function vector contains all the constraints defined in Section 3. Differentiating with respect to yieldsand the submatrix can be expressed asin which
The gradient of distance constraint function is given by
Hessian matrix for is given bywherein which
Then, the gradient of velocity correlation constraint function for is given bywhere
The Hessian matrix for is given bywhere
The details of and are given by
Finally, the gradient of the X-axis velocity constraint function is given by
In Hessian matrix , the element which is on , is equal to unity, and , are equal to negative unity, otherwise zero.
Data Availability
The data can be obtained from the corresponding author of this paper.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge support from National Natural Science Foundation of China (Grant nos. 61772548, 61401513, and 61201381), China Postdoctoral Science Foundation (Grant no. 2016M592989), the Self-Topic Foundation of Information Engineering University (Grant no. 2016600701), and the Outstanding Youth Foundation of Information Engineering University (Grant no. 2016603201).