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

Authors: Bin Yang; Zeyu Yang; Ding Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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 Section2, 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.Table 1
Mathematical notation.
NotationExplanationON×nN×n matrix with all-zero entriesONN×N matrix with all-zero entriesINN dimensional identity matrix11×N or 1N×11×N or N×1 vector with all-one entriesblkdiag⋅Block diagonal matrix consisting of the matrix or vector⋅2The euclidean norm of a vectoriNnN×N matrix of zeros except for the ith diagonal element is equal to unity⋅†The pseudo-inverse of a matrix that satisfies the Moore-Penrose conditions⊗Kronecker product


## 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, whereM moving sensors are used to determine the N moving sources of interest with unknown positions un=xnynznT and velocities u˙n=x˙ny˙nz˙nT. The three-dimensional scenario can be straightly transformed into two-dimensional scenario generally. The location of the nth source can be expressed as u¯n=unTu˙nTT and the composite source location vector can be written as u¯=u¯1Tu¯2T…u¯NTT. The sensor positions sm=xmymzmT and velocities s˙m=x˙my˙mz˙mT are assumed known perfectly, and we can collect the sensor positions and velocities as the vector s¯m=smTs˙mTT. The composite sensor location vector can be written as s¯=s¯1Ts¯2T…s¯MTT. 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, sensor1 can be set as the reference sensor. Let rm,np be the true distance between the source n and sensor m:(1)rm,np=un−sm2.Then, the TDOA measurements between sensor pairm and 1 is(2)r^m1,np=rm1,np+εm1,np=rm,np−r1,np+εm1,np=fmpu¯n,s¯+εm1,np2≤m≤M1≤n≤N,where fmpu¯n,s¯ is a function with respect to un and s¯, εm1,np denotes the TDOA noise. Although velocity u˙n is useless here, u¯n is still used rather than un 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:(3)r^np=r^21,npr^31,np…r^M1,npT=rnp+εnp=fpu¯n,s¯+εnp1≤n≤N,where(4)rnp=r21,npr31,np…rM1,npT,εnp=ε21,npε31,np…εM1,npT,fpu¯n,s¯=f2pu¯n,s¯f3pu¯n,s¯…fMpu¯n,s¯T.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(5)r^m1,nv=rm1,nv+εm1,nv=rm,nv−r1,nv+εm1,nv=un−smTu˙n−s˙mun−sm2−un−s1Tu˙n−s˙1un−s12+εm1,nv=fmvu¯n,s¯+εm1,nv2≤m≤M1≤n≤N,where εm1,nv is the FDOA noise; fmvu¯n,s¯ is a function with respect to un and s¯. Like the expression of TDOA, the composite vector of the whole sensor pairs from source n can be written as(6)r^nv=r^21,nvr^31,nv…r^M1,nvT=rnv+εnv=fvu¯n,s¯+εnv1≤n≤N,where(7)rnv=r21,nvr31,nv…rM1,nvT,εnv=ε21,nvε31,nv…εM1,nvT,fvu¯n,s¯=f2vu¯n,s¯f3vu¯n,s¯…fMvu¯n,s¯T.Combining TDOA and FDOA of sourcen together, the measurements vector is given:(8)r^n=r^npTr^nvTT=fu¯n,s¯+εn,where(9)rn=rnpTrnvTT,εn=εnpTεnvTT,rn=fu¯n,s¯=fpTu¯n,s¯fvTu¯n,s¯T,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:(10)Qn=EεnεnT=QnpQnvpQnpvQnv,where Qnp=EεnpεnpT, Qnv=EεnvεnvT, and the cross-correlation between the TDOA and FDOA noise can be characterized:(11)Qnpv=EεnpεnvT,Qnvp=EεnvεnpT,Qnpv=QnpvT.Then putting the all measurements together, TDOA and FDOA measurements fromN sources yields(12)r¯^=r^1Tr^2T…r^NTT=r¯+ε¯=f¯u¯,s¯+ε¯,where(13)r¯=r1Tr2T…rNTT,ε¯=ε1Tε2T…εNTT,r¯=f¯u¯,s¯=fu¯1,s¯Tfu¯2,s¯T…fu¯N,s¯T.The composite vectorε¯ follows zero-mean Gaussian distribution with covariance matrix:(14)Q=Eε¯ε¯T=blkdiagQ1Q2…QN.


### 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. Figure1 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.Figure 1
Location scenario.In practice, prior range constraints must exist for the multiple sources that accompany each other. As is shown in Figure2, when we have N sources in a formation, there are at least NN−1/2 distance constraints among them. Using prior information such as the safe distance, the two sources n1 and n2 can be limited to a range dn1,n2 that is greater than or equal to the true distance. A normal distance constraint between source n1 and n2 can be described as(15)un1−un22≤dn1,n21≤n1<n2≤N.Figure 2
Diagram for distance constraints.To represent the position and velocity of sourcen with vector u¯, the matrix Jnp=O3×6n−1I3O3×3O3×6N−n is defined to select the parameters. In addition, the position of source n can be written as Jnpu¯. The distance constraint can be rewritten as(16)Jn1pu¯−Jn2pu¯2≤dn1,n21≤n1<n2≤N.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 Figure3, which is given by(17)u˙n1Tu˙n2u˙n12u˙n22=Jn1vu¯TJn2vu¯Jn1vu¯2Jn2vu¯2≥ρn1,n21≤n1<n2≤N.Figure 3
Diagram for velocity constraints.The velocity of sourcen can be written as Jnvu¯ where the matrix Jnv=O3×6n−1O3×3I3O3×6N−n.In contrast to the distance constraints, the true velocity correlation should be larger thanρn1,n2 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.Figure 4
Diagram for Velocity Correlation. (a) Same direction. (b) Opposite direction.
(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(18)x˙n1−x˙n2≤Δx˙n1,n21≤n1<n2≤N.Using the matrixi31, the auxiliary constraint can be rewritten as(19)i31Jn1vu¯−Jn2vu¯2≤Δx˙n1,n21≤n1<n2≤N.We choose theX-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:(20)i32Jn1vu¯−Jn2vu¯2≤Δy˙n1,n2i33Jn1vu¯−Jn2vu¯2≤Δz˙n1,n2,1≤n1<n2≤N.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 theX-direction are sufficient to determine the velocity correlation.


## 2.1. Signal Model for Localization Based on TDOA and FDOA Measurements

In this paper, we consider a three-dimensional localization scenario, whereM moving sensors are used to determine the N moving sources of interest with unknown positions un=xnynznT and velocities u˙n=x˙ny˙nz˙nT. The three-dimensional scenario can be straightly transformed into two-dimensional scenario generally. The location of the nth source can be expressed as u¯n=unTu˙nTT and the composite source location vector can be written as u¯=u¯1Tu¯2T…u¯NTT. The sensor positions sm=xmymzmT and velocities s˙m=x˙my˙mz˙mT are assumed known perfectly, and we can collect the sensor positions and velocities as the vector s¯m=smTs˙mTT. The composite sensor location vector can be written as s¯=s¯1Ts¯2T…s¯MTT. 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, sensor1 can be set as the reference sensor. Let rm,np be the true distance between the source n and sensor m:(1)rm,np=un−sm2.Then, the TDOA measurements between sensor pairm and 1 is(2)r^m1,np=rm1,np+εm1,np=rm,np−r1,np+εm1,np=fmpu¯n,s¯+εm1,np2≤m≤M1≤n≤N,where fmpu¯n,s¯ is a function with respect to un and s¯, εm1,np denotes the TDOA noise. Although velocity u˙n is useless here, u¯n is still used rather than un 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:(3)r^np=r^21,npr^31,np…r^M1,npT=rnp+εnp=fpu¯n,s¯+εnp1≤n≤N,where(4)rnp=r21,npr31,np…rM1,npT,εnp=ε21,npε31,np…εM1,npT,fpu¯n,s¯=f2pu¯n,s¯f3pu¯n,s¯…fMpu¯n,s¯T.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(5)r^m1,nv=rm1,nv+εm1,nv=rm,nv−r1,nv+εm1,nv=un−smTu˙n−s˙mun−sm2−un−s1Tu˙n−s˙1un−s12+εm1,nv=fmvu¯n,s¯+εm1,nv2≤m≤M1≤n≤N,where εm1,nv is the FDOA noise; fmvu¯n,s¯ is a function with respect to un and s¯. Like the expression of TDOA, the composite vector of the whole sensor pairs from source n can be written as(6)r^nv=r^21,nvr^31,nv…r^M1,nvT=rnv+εnv=fvu¯n,s¯+εnv1≤n≤N,where(7)rnv=r21,nvr31,nv…rM1,nvT,εnv=ε21,nvε31,nv…εM1,nvT,fvu¯n,s¯=f2vu¯n,s¯f3vu¯n,s¯…fMvu¯n,s¯T.Combining TDOA and FDOA of sourcen together, the measurements vector is given:(8)r^n=r^npTr^nvTT=fu¯n,s¯+εn,where(9)rn=rnpTrnvTT,εn=εnpTεnvTT,rn=fu¯n,s¯=fpTu¯n,s¯fvTu¯n,s¯T,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:(10)Qn=EεnεnT=QnpQnvpQnpvQnv,where Qnp=EεnpεnpT, Qnv=EεnvεnvT, and the cross-correlation between the TDOA and FDOA noise can be characterized:(11)Qnpv=EεnpεnvT,Qnvp=EεnvεnpT,Qnpv=QnpvT.Then putting the all measurements together, TDOA and FDOA measurements fromN sources yields(12)r¯^=r^1Tr^2T…r^NTT=r¯+ε¯=f¯u¯,s¯+ε¯,where(13)r¯=r1Tr2T…rNTT,ε¯=ε1Tε2T…εNTT,r¯=f¯u¯,s¯=fu¯1,s¯Tfu¯2,s¯T…fu¯N,s¯T.The composite vectorε¯ follows zero-mean Gaussian distribution with covariance matrix:(14)Q=Eε¯ε¯T=blkdiagQ1Q2…QN.


## 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. Figure1 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.Figure 1
Location scenario.In practice, prior range constraints must exist for the multiple sources that accompany each other. As is shown in Figure2, when we have N sources in a formation, there are at least NN−1/2 distance constraints among them. Using prior information such as the safe distance, the two sources n1 and n2 can be limited to a range dn1,n2 that is greater than or equal to the true distance. A normal distance constraint between source n1 and n2 can be described as(15)un1−un22≤dn1,n21≤n1<n2≤N.Figure 2
Diagram for distance constraints.To represent the position and velocity of sourcen with vector u¯, the matrix Jnp=O3×6n−1I3O3×3O3×6N−n is defined to select the parameters. In addition, the position of source n can be written as Jnpu¯. The distance constraint can be rewritten as(16)Jn1pu¯−Jn2pu¯2≤dn1,n21≤n1<n2≤N.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 Figure3, which is given by(17)u˙n1Tu˙n2u˙n12u˙n22=Jn1vu¯TJn2vu¯Jn1vu¯2Jn2vu¯2≥ρn1,n21≤n1<n2≤N.Figure 3
Diagram for velocity constraints.The velocity of sourcen can be written as Jnvu¯ where the matrix Jnv=O3×6n−1O3×3I3O3×6N−n.In contrast to the distance constraints, the true velocity correlation should be larger thanρn1,n2 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.Figure 4
Diagram for Velocity Correlation. (a) Same direction. (b) Opposite direction.
(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(18)x˙n1−x˙n2≤Δx˙n1,n21≤n1<n2≤N.Using the matrixi31, the auxiliary constraint can be rewritten as(19)i31Jn1vu¯−Jn2vu¯2≤Δx˙n1,n21≤n1<n2≤N.We choose theX-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:(20)i32Jn1vu¯−Jn2vu¯2≤Δy˙n1,n2i33Jn1vu¯−Jn2vu¯2≤Δz˙n1,n2,1≤n1<n2≤N.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 theX-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(21)CRLBu¯=FIM−1=∂r¯∂u¯TTQ−1∂r¯∂u¯T−1.The details of∂r¯/∂u¯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 Θ=ℝ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. 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 Θ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 n1 and n2 can be rewritten as functions(22)gn1,n2c1u¯=Jn1pu¯−Jn2pu¯2−dn1,n2≤0,gn1,n2c2u¯=−Jn1pu¯TJn2pu¯Jn1pu¯Jn2pu¯2+ρn1,n2≤0,gn1,n2c3u¯=i31Jn1vu¯−Jn2vu¯2−Δx˙n1,n2≤0,1≤n1<n2≤Nwhere cii=1,2,3 represents the three kind of constraints respectively.There areNN−1/2 functions about distance constraint overall. For example, the distance constraint of sensor n is given in the vector(23)gnc1u¯=gn,n+1c1u¯gn,n+2c1u¯…gn,Nc1u¯T1≤n≤N−1.Then, we have the function vector putting all the distance constraints together:(24)g¯c1u¯=g1c1Tu¯g2c1Tu¯…gN−1c1Tu¯T.It is also true for the other constraints, which can be written as(25)g¯ciu¯=g1ciTu¯g2ciTu¯…gN−1ciTu¯T,where(26)gnciu¯=gn,n+1ciu¯gn,n+2ciu¯…gn,Nciu¯T1≤n≤N−1.Finally, combining (24) and (25), the composite vector is given(27)Gu¯=g¯c1Tu¯g¯c2Tu¯g¯c3Tu¯T.As is proved in Theorem 1 [40], the 3NN−1/2×6N gradient matrix ∇Gu¯, whose rank is smaller than the number of constraints 3NN−1/2, error covariance matrix Σu¯ satisfies(28)∑u¯≥CRLBCu¯,where(29)CRLBCu¯=Qu¯CRLBu¯.MatrixQu¯ is defined as(30)Qu¯=I−CRLBu¯∇Gu¯T∇Gu¯CRLBu¯∇Gu¯T†∇Gu¯,where ∇Gu¯ is the gradient matrix of (27). The details of ∇Gu¯ are given in Appendix B. In (29), matrix Qu¯ contains the prior information of multiple sources location. Then, substituting (30) into (29), the expression of CRLBCu¯ is obtained(31)CRLBCu¯=CRLBu¯−CRLBu¯∇Gu¯T∇Gu¯CRLBu¯∇Gu¯T†∇Gu¯CRLBu¯.MatrixQu¯ clearly can reduce the original CRLB, because the second term of (31) is a semidefinite matrix. The constrained CRLB is called CRLBCu¯ here.Obviously, when the constraints are active, which means thatu¯ 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. 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 Figure5, 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.Figure 5
CRLB and constrained CRLBs as the number of sources varies.
(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ε¯TQ−1ε¯. With (22), we can establish the following constrained estimator:(32)minu¯∈ℝ6N×1ε¯TQ−1ε¯s.t.gn1,n2ciu¯≤01≤n1<n2≤N,1≤i≤3.For allN sources in a group, there will be 3NN−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.


### 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 u¯ everywhere. We substitute the constraint functions into an exponential function:(33)φt=et−1.So exponential penalty functions are produced(34)ℓu¯,λci,τci=∑1≤n1<n2≤Nλn1,n2ciτn1,n2ciφτn1,n2cign1,n2ciu¯,where(35)λci=λ1ciTλ2ciT…λN−1ciTT,τci=τ1ciTτ2ciT…τN−1ciTT,in which(36)λnci=λn,n+1ciλn,n+2ci…λn,NciTτnci=τn,n+1ciτn,n+2ci…τn,NciT1≤n≤N−1.Then, the exponential Lagrangian function is given(37)Lu¯,λ¯,τ¯=ε¯TQ−1ε¯+ℓu¯,λ¯,τ¯,where(38)ℓu¯,λ¯,τ¯=∑i=13ℓu¯,λci,τci,in which(39)λ¯=λc1Tλc2Tλc3TT,τ¯=τc1Tτc2Tτc3TT.FunctionLu¯,λ¯,τ¯ is the new cost function. Overall, the constrained estimator is transformed into an unconstrained estimator, and we just need to compute u¯^ as(40)u¯^=argminu¯∈ℝ6N×1Lu¯,λ¯,τ¯=argminu¯∈ℝ6N×1Lu¯.Remark 1.
Note that (33) has following properties:(41)(i)For all thet∈ℝ,∇2φt>0;iiφ0=0,∇φ0>0;iiiWhent<0,φt<0,whent>0,φt>0.
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λn1,n2ci>0 as a multiplier and τn1,n2ci as a positive penalty parameter with the constraint function gn1,n2ciu¯. 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):(42)λn1,n2cik+1=λn1,n2cik∇φτn1,n2cign1,n2ciu¯k+1=λn1,n2cikexpτn1,n2cign1,n2ciu¯k+1.Remark 2.
Note that for a fixedλn1,n2ci>0, as τn1,n2ci⟶∞, penalty term λn1,n2ciφτn1,n2cign1,n2ciu¯/τn1,n2ci tends to ∞ for all infeasible u¯ and to zero for all feasible u¯; otherwise, for a fixed τn1,n2ci, as λn1,n2ci⟶0 the corresponding penalty term goes to zero whether u¯ 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 Section5 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 ε¯TQ−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 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. DifferentiatingLu¯ with respect to u¯, gradient vector is obtained:(43)∇Lu¯=∂Lu¯∂u¯=2−∂r¯∂u¯TTQ−1r¯^−r¯+∂ℓu¯,λ¯,τ¯∂u¯,where(44)∂ℓu¯,λ¯,τ¯∂u¯=∂ℓu¯∂u¯=∑i=13∑1≤n1<n2≤Nλn1,n2ciexpτn1,n2cign1,n2ciu¯∂gn1,n2ciu¯∂u¯.Hessian matrix ofLu¯ is given by(45)∇2Lu¯=∂2Lu¯∂u¯∂u¯T=2∂r¯∂u¯TTQ−1∂r¯∂u¯T+2r¯^−r¯TQ−1⊗I6N∂∂u¯Tvec−∂r¯∂u¯TT+∂2ℓu¯,λ¯,τ¯∂u¯∂u¯T,(46)∂2ℓu¯,λ¯,τ¯∂u¯∂u¯T=∂2ℓu¯∂u¯∂u¯T=∑i=13∑1≤n1<n2≤Nλn1,n2ciexpτn1,n2cign1,n2ciu¯τn1,n2ci∂n1,n2ciu¯∂u¯T∂gn1,n2ciu¯∂u¯+∂2gn1,n2ciu¯∂u¯∂u¯T.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 Table2.Table 2
Process of proposed algorithm.
Initialization and Input:u¯0, initial value for Newton’s iteration.r¯^, TDOA and FDOA measurements.λ¯0, positive multiplier vector.τ¯ , 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).Preparation:k=0Lmin=Lu¯0,λ¯0,τ¯ , initialize the value for comparison and save gradient vector ∇Lu¯0,λ¯0,τ¯ and hessian matrix ∇2Lu¯0,λ¯0,τ¯.Iteration process:Fork+1≤η:Step 1: Update the estimate byu¯k+1=u¯k−μk∇2Lu¯k,λ¯k,τ¯−1∇Lu¯k,λ¯k,τ¯ and turn to step 2.Step 2: if Lu¯k+1<Lmin, renew Lmin:=Lu¯k+1, and update the multipliers by using (42), then k:=k+1, record the updated ∇Lu¯k,λ¯k,τ¯ and ∇2Lu¯k,λ¯k,τ¯, turn to step 4.Step 3: if Lu¯k+1,λ¯k+1,τ¯≥Lmin, μk+1:=σμk, thenk:=k+1, turn to step 4.Step 4: when ∇Lu¯k+1,λ¯k+1,τ¯2≤ξ, output u¯k+1, otherwise, turn to step 1.EndOutput:u¯final as the optimal solutionRemark 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 Section5 prove that this initial guess can provide good results.


## 4.1. Constrained Estimator for Source Localization

The ML estimator of the sources’ position and velocity is to minimizeε¯TQ−1ε¯. With (22), we can establish the following constrained estimator:(32)minu¯∈ℝ6N×1ε¯TQ−1ε¯s.t.gn1,n2ciu¯≤01≤n1<n2≤N,1≤i≤3.For allN sources in a group, there will be 3NN−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.


## 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 u¯ everywhere. We substitute the constraint functions into an exponential function:(33)φt=et−1.So exponential penalty functions are produced(34)ℓu¯,λci,τci=∑1≤n1<n2≤Nλn1,n2ciτn1,n2ciφτn1,n2cign1,n2ciu¯,where(35)λci=λ1ciTλ2ciT…λN−1ciTT,τci=τ1ciTτ2ciT…τN−1ciTT,in which(36)λnci=λn,n+1ciλn,n+2ci…λn,NciTτnci=τn,n+1ciτn,n+2ci…τn,NciT1≤n≤N−1.Then, the exponential Lagrangian function is given(37)Lu¯,λ¯,τ¯=ε¯TQ−1ε¯+ℓu¯,λ¯,τ¯,where(38)ℓu¯,λ¯,τ¯=∑i=13ℓu¯,λci,τci,in which(39)λ¯=λc1Tλc2Tλc3TT,τ¯=τc1Tτc2Tτc3TT.FunctionLu¯,λ¯,τ¯ is the new cost function. Overall, the constrained estimator is transformed into an unconstrained estimator, and we just need to compute u¯^ as(40)u¯^=argminu¯∈ℝ6N×1Lu¯,λ¯,τ¯=argminu¯∈ℝ6N×1Lu¯.Remark 1.
Note that (33) has following properties:(41)(i)For all thet∈ℝ,∇2φt>0;iiφ0=0,∇φ0>0;iiiWhent<0,φt<0,whent>0,φt>0.
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λn1,n2ci>0 as a multiplier and τn1,n2ci as a positive penalty parameter with the constraint function gn1,n2ciu¯. 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):(42)λn1,n2cik+1=λn1,n2cik∇φτn1,n2cign1,n2ciu¯k+1=λn1,n2cikexpτn1,n2cign1,n2ciu¯k+1.Remark 2.
Note that for a fixedλn1,n2ci>0, as τn1,n2ci⟶∞, penalty term λn1,n2ciφτn1,n2cign1,n2ciu¯/τn1,n2ci tends to ∞ for all infeasible u¯ and to zero for all feasible u¯; otherwise, for a fixed τn1,n2ci, as λn1,n2ci⟶0 the corresponding penalty term goes to zero whether u¯ 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 Section5 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 ε¯TQ−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 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. DifferentiatingLu¯ with respect to u¯, gradient vector is obtained:(43)∇Lu¯=∂Lu¯∂u¯=2−∂r¯∂u¯TTQ−1r¯^−r¯+∂ℓu¯,λ¯,τ¯∂u¯,where(44)∂ℓu¯,λ¯,τ¯∂u¯=∂ℓu¯∂u¯=∑i=13∑1≤n1<n2≤Nλn1,n2ciexpτn1,n2cign1,n2ciu¯∂gn1,n2ciu¯∂u¯.Hessian matrix ofLu¯ is given by(45)∇2Lu¯=∂2Lu¯∂u¯∂u¯T=2∂r¯∂u¯TTQ−1∂r¯∂u¯T+2r¯^−r¯TQ−1⊗I6N∂∂u¯Tvec−∂r¯∂u¯TT+∂2ℓu¯,λ¯,τ¯∂u¯∂u¯T,(46)∂2ℓu¯,λ¯,τ¯∂u¯∂u¯T=∂2ℓu¯∂u¯∂u¯T=∑i=13∑1≤n1<n2≤Nλn1,n2ciexpτn1,n2cign1,n2ciu¯τn1,n2ci∂n1,n2ciu¯∂u¯T∂gn1,n2ciu¯∂u¯+∂2gn1,n2ciu¯∂u¯∂u¯T.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 Table2.Table 2
Process of proposed algorithm.
Initialization and Input:u¯0, initial value for Newton’s iteration.r¯^, TDOA and FDOA measurements.λ¯0, positive multiplier vector.τ¯ , 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).Preparation:k=0Lmin=Lu¯0,λ¯0,τ¯ , initialize the value for comparison and save gradient vector ∇Lu¯0,λ¯0,τ¯ and hessian matrix ∇2Lu¯0,λ¯0,τ¯.Iteration process:Fork+1≤η:Step 1: Update the estimate byu¯k+1=u¯k−μk∇2Lu¯k,λ¯k,τ¯−1∇Lu¯k,λ¯k,τ¯ and turn to step 2.Step 2: if Lu¯k+1<Lmin, renew Lmin:=Lu¯k+1, and update the multipliers by using (42), then k:=k+1, record the updated ∇Lu¯k,λ¯k,τ¯ and ∇2Lu¯k,λ¯k,τ¯, turn to step 4.Step 3: if Lu¯k+1,λ¯k+1,τ¯≥Lmin, μk+1:=σμk, thenk:=k+1, turn to step 4.Step 4: when ∇Lu¯k+1,λ¯k+1,τ¯2≤ξ, output u¯k+1, otherwise, turn to step 1.EndOutput:u¯final as the optimal solutionRemark 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 Section5 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 σd2R and σf2R, respectively. The noise power σd2 and σf2 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 σd2R and σf2R, which is one of the typical values of the TDOA and FDOA estimation in practice. The matrix R is defined as a M−1×M−1 matrix, the diagonal elements of which are equal to unity and 0.5 otherwise.(47)Qn=σd2ROM−1OM−1σf2R.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:(48)RMSEun=15000∑i=15000Jnpu¯^−Jnpu¯2,RMSEu˙n=15000∑i=15000Jnvu¯^−Jnvu¯2.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 Table3. The locations of the sensors and unknown sources can be seen in Figure 6 in which the arrows represent their velocity and direction.Table 3
Locations of the sensors.
Sensor no.mSensor network123456x150−100200−150−100−300y3001500200400250z–200200100100350100x˙102025-101520y˙20−30101520−15z˙−155−202530−10Figure 6
Geometry of sensor network and sources location.We assume that there are two moving sources whose positions and velocities areu1=654803708T m, u˙1=20.2315.1840.34T m/s and u2=652808710T m, u˙2=201538T m/s, respectively. The sensor network is employed to locate the two sources. In this case, we set σf2=0.25σd2 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 ΘC. 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 d1,2=10 m, ρ1,2=0.999, and Δx˙1,2=0.5 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 ρ1,2=0.999. The Lagrangian multiplier vector λ¯ is set to 211T and penalty parameter vector τ¯=101010T initially.Figures7 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.Figure 7
Performance of source 1 in near field.
(a)(b)Figure 8
Performance of source 2 in near field.
(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 areu1=126014401193T m, u2=125614501195T m, and u3=125114561188T m. Their velocities are u˙1=20.2315.1840.34T m/s, u˙2=201538T m/s, and u˙3=19.7414.7941.39T 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.Table 4
Locations of the sensors.
Sensor no.mSensor network123456x300400300350−100200y100150500200−100−300z150100200100−100−200x˙30−301010−20020y˙−2010−20−2010010z˙2020103010010Figure 9
Geometry of sensor network and sources location.Here, the covariance matrixσd2R of TDOA is assumed to be ten times σf2R 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 dn1,n2=25 m, ρn1,n2=0.999, and Δx˙1,2=1 m/s for the three constraints. The Lagrangian multiplier λ¯ is set to 211×311×311×3T and penalty parameter vector τ¯=1011×31011×31011×3T which are the same as those in near-field scenarios. Figures 10–12 depict the results.Figure 10
Performance of source 1 in far field.
(a)(b)Figure 11
Performance of source 2 in far field.
(a)(b)Figure 12
Performance of source 3 in far field.
(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 ofR=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 γ of −2π/9 rad to 2π/9 rad with an angular velocity of 0.02 rad/s. The distance between the two sensors is 10 m from beginning to end, and we set d1,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. The TDOA and FDOA noise powers are σd2=0.025/c2 and σf2=0.0025/c2, respectively. The stationary sensor positions are listed in Table 5.Table 5
Location of the sensors.
Sensor no.m123456x0400−400000y00040000z0000400−400We conducted 5000 Monte Carlo experiments for every point in its trajectory. The Lagrangian multiplier vectorλ¯ was set to 211T and penalty parameter vector was τ¯=101010T at the beginning. The proposed method is compared with the existing methods.As shown in Figure13, 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.Figure 13
Standard deviation for the methods to localize two far-field sources. (a) Performance of source 1. (b) Performance of source 2.
(a)(b)The iterative processes of the Lagrangian multipliers and the Euclidean norm of the gradient vector at the starting point of Figure14, γ=−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. The trend of the lines indicates that the proposed method can converge within 10 iterations, which shows the convergence is also good and fast.Figure 14
Geometry of sensors and sources location.Figure 15
Lagrangian multiplier and the Euclidean norm of gradient vector. (a) The Lagrangian multiplierλn1,n2c1 for distance constraint. (b) The Lagrangian multiplier λn1,n2c2 for X-axis component of velocity. (c) The Lagrangian multiplier λn1,n2c3 for velocity constraint. (d) The Euclidean norm of gradient vector ∇Lu¯k+1,λ¯k+1,τ¯2.
(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 Table6. The location of geometry is shown in Figure 16. The true starting coordinates of the formation sources are u1=000T m, u2=7−103T m, and u3=2101T 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 203/3203/3203/3T m/s and their acceleration is 23/323/3−23/3T m/s2. During the movement, they remain relatively static. The TDOA noise variance was σd2=0.25/c2 and FDOA noise was σf2=0.025/c2 in this case. The Lagrangian multiplier is λ¯=211T and penalty parameter is τ¯=101010T. The conditions of the constraints are set as dn1,n2=25 m, ρn1,n2=0.999, and Δx˙1,2=1 m/s. The simulation results are obtained via 5000 Monte Carlo experiments.Table 6
Location of the sensors.
Sensor no.m123456x200−2000000y00200−20000z0000200−200Figure 16
Geometry of sensors and sources location changes over time.Figures17–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, 0 0 0T 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.Figure 17
Performance of source 1 changes over time.
(a)(b)Figure 18
Performance of source 2 changes over time.
(a)(b)Figure 19
Performance of source 3 changes over time.
(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 Sections5.1 and 5.2, respectively. The results of average run time of one Monte Carlo experiment are listed in Table 7.Table 7
Comparison of run time.
Localization methodAverage CPU computational time (ms)2 sources (near-field)3 sources (far-field)TW-WLS method1.672.38ICWLS method1.492.17Gauss–Newton method2.893.85Proposed method3.823.57With 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, Table7 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.


## 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 Table3. The locations of the sensors and unknown sources can be seen in Figure 6 in which the arrows represent their velocity and direction.Table 3
Locations of the sensors.
Sensor no.mSensor network123456x150−100200−150−100−300y3001500200400250z–200200100100350100x˙102025-101520y˙20−30101520−15z˙−155−202530−10Figure 6
Geometry of sensor network and sources location.We assume that there are two moving sources whose positions and velocities areu1=654803708T m, u˙1=20.2315.1840.34T m/s and u2=652808710T m, u˙2=201538T m/s, respectively. The sensor network is employed to locate the two sources. In this case, we set σf2=0.25σd2 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 ΘC. 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 d1,2=10 m, ρ1,2=0.999, and Δx˙1,2=0.5 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 ρ1,2=0.999. The Lagrangian multiplier vector λ¯ is set to 211T and penalty parameter vector τ¯=101010T initially.Figures7 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.Figure 7
Performance of source 1 in near field.
(a)(b)Figure 8
Performance of source 2 in near field.
(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 areu1=126014401193T m, u2=125614501195T m, and u3=125114561188T m. Their velocities are u˙1=20.2315.1840.34T m/s, u˙2=201538T m/s, and u˙3=19.7414.7941.39T 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.Table 4
Locations of the sensors.
Sensor no.mSensor network123456x300400300350−100200y100150500200−100−300z150100200100−100−200x˙30−301010−20020y˙−2010−20−2010010z˙2020103010010Figure 9
Geometry of sensor network and sources location.Here, the covariance matrixσd2R of TDOA is assumed to be ten times σf2R 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 dn1,n2=25 m, ρn1,n2=0.999, and Δx˙1,2=1 m/s for the three constraints. The Lagrangian multiplier λ¯ is set to 211×311×311×3T and penalty parameter vector τ¯=1011×31011×31011×3T which are the same as those in near-field scenarios. Figures 10–12 depict the results.Figure 10
Performance of source 1 in far field.
(a)(b)Figure 11
Performance of source 2 in far field.
(a)(b)Figure 12
Performance of source 3 in far field.
(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 ofR=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 γ of −2π/9 rad to 2π/9 rad with an angular velocity of 0.02 rad/s. The distance between the two sensors is 10 m from beginning to end, and we set d1,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. The TDOA and FDOA noise powers are σd2=0.025/c2 and σf2=0.0025/c2, respectively. The stationary sensor positions are listed in Table 5.Table 5
Location of the sensors.
Sensor no.m123456x0400−400000y00040000z0000400−400We conducted 5000 Monte Carlo experiments for every point in its trajectory. The Lagrangian multiplier vectorλ¯ was set to 211T and penalty parameter vector was τ¯=101010T at the beginning. The proposed method is compared with the existing methods.As shown in Figure13, 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.Figure 13
Standard deviation for the methods to localize two far-field sources. (a) Performance of source 1. (b) Performance of source 2.
(a)(b)The iterative processes of the Lagrangian multipliers and the Euclidean norm of the gradient vector at the starting point of Figure14, γ=−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. The trend of the lines indicates that the proposed method can converge within 10 iterations, which shows the convergence is also good and fast.Figure 14
Geometry of sensors and sources location.Figure 15
Lagrangian multiplier and the Euclidean norm of gradient vector. (a) The Lagrangian multiplierλn1,n2c1 for distance constraint. (b) The Lagrangian multiplier λn1,n2c2 for X-axis component of velocity. (c) The Lagrangian multiplier λn1,n2c3 for velocity constraint. (d) The Euclidean norm of gradient vector ∇Lu¯k+1,λ¯k+1,τ¯2.
(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 Table6. The location of geometry is shown in Figure 16. The true starting coordinates of the formation sources are u1=000T m, u2=7−103T m, and u3=2101T 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 203/3203/3203/3T m/s and their acceleration is 23/323/3−23/3T m/s2. During the movement, they remain relatively static. The TDOA noise variance was σd2=0.25/c2 and FDOA noise was σf2=0.025/c2 in this case. The Lagrangian multiplier is λ¯=211T and penalty parameter is τ¯=101010T. The conditions of the constraints are set as dn1,n2=25 m, ρn1,n2=0.999, and Δx˙1,2=1 m/s. The simulation results are obtained via 5000 Monte Carlo experiments.Table 6
Location of the sensors.
Sensor no.m123456x200−2000000y00200−20000z0000200−200Figure 16
Geometry of sensors and sources location changes over time.Figures17–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, 0 0 0T 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.Figure 17
Performance of source 1 changes over time.
(a)(b)Figure 18
Performance of source 2 changes over time.
(a)(b)Figure 19
Performance of source 3 changes over time.
(a)(b)


## 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 ofR=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 γ of −2π/9 rad to 2π/9 rad with an angular velocity of 0.02 rad/s. The distance between the two sensors is 10 m from beginning to end, and we set d1,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. The TDOA and FDOA noise powers are σd2=0.025/c2 and σf2=0.0025/c2, respectively. The stationary sensor positions are listed in Table 5.Table 5
Location of the sensors.
Sensor no.m123456x0400−400000y00040000z0000400−400We conducted 5000 Monte Carlo experiments for every point in its trajectory. The Lagrangian multiplier vectorλ¯ was set to 211T and penalty parameter vector was τ¯=101010T at the beginning. The proposed method is compared with the existing methods.As shown in Figure13, 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.Figure 13
Standard deviation for the methods to localize two far-field sources. (a) Performance of source 1. (b) Performance of source 2.
(a)(b)The iterative processes of the Lagrangian multipliers and the Euclidean norm of the gradient vector at the starting point of Figure14, γ=−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. The trend of the lines indicates that the proposed method can converge within 10 iterations, which shows the convergence is also good and fast.Figure 14
Geometry of sensors and sources location.Figure 15
Lagrangian multiplier and the Euclidean norm of gradient vector. (a) The Lagrangian multiplierλn1,n2c1 for distance constraint. (b) The Lagrangian multiplier λn1,n2c2 for X-axis component of velocity. (c) The Lagrangian multiplier λn1,n2c3 for velocity constraint. (d) The Euclidean norm of gradient vector ∇Lu¯k+1,λ¯k+1,τ¯2.
(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 Table6. The location of geometry is shown in Figure 16. The true starting coordinates of the formation sources are u1=000T m, u2=7−103T m, and u3=2101T 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 203/3203/3203/3T m/s and their acceleration is 23/323/3−23/3T m/s2. During the movement, they remain relatively static. The TDOA noise variance was σd2=0.25/c2 and FDOA noise was σf2=0.025/c2 in this case. The Lagrangian multiplier is λ¯=211T and penalty parameter is τ¯=101010T. The conditions of the constraints are set as dn1,n2=25 m, ρn1,n2=0.999, and Δx˙1,2=1 m/s. The simulation results are obtained via 5000 Monte Carlo experiments.Table 6
Location of the sensors.
Sensor no.m123456x200−2000000y00200−20000z0000200−200Figure 16
Geometry of sensors and sources location changes over time.Figures17–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, 0 0 0T 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.Figure 17
Performance of source 1 changes over time.
(a)(b)Figure 18
Performance of source 2 changes over time.
(a)(b)Figure 19
Performance of source 3 changes over time.
(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 Sections5.1 and 5.2, respectively. The results of average run time of one Monte Carlo experiment are listed in Table 7.Table 7
Comparison of run time.
Localization methodAverage CPU computational time (ms)2 sources (near-field)3 sources (far-field)TW-WLS method1.672.38ICWLS method1.492.17Gauss–Newton method2.893.85Proposed method3.823.57With 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, Table7 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.



---

*Source: 1382321-2020-06-16.xml*