Approximate MLE Algorithm for Source Localization Based

Approximate MLE Algorithm for Source Localization Basedon TDOA Measurements∗Guoxiang GuDepartment of Electrical and Computer EngineeringLouisiana State University,Baton Rouge,LA70803-5901July2006,Revised January2007ABSTRACTSource localization is investigated based on noisy measurements of TDOA(time-difference of arrival)in which the measurement noises are assumed to be Gauss distributed.The solution to the constrained WLS(weighted least-squares)is derived and applied to the source localization problem based on TDOAs.The proposed algorithm is shown to be an approximate MLE(maximum likelihood estimation)algorithm under some mild condition.The simulation results show that the proposed approximate MLE algorithm compares favorably with the existing solution methods for source localization based on TDOA measurements.Keywords:Passive source localization,TDOA,AOA,MLE,constrained LS1.INTRODUCTIONPassive source localization has been an important research topic in the signal processing society.It has received greater attention in recent years due to the renewed interest in wireless location(cf.Sayed et al.11and the references therein).This paper considers source localization based on TDOA measurements that has been studied in several papers from Abel and Smith1to Friedlander,6in Hahn and Tretter,8in papers from Ho and Xu10to Schau and Robinson,12and in Torrieri.15Despite the fact that the measurement noises are additive and Gaussian,such a source localization problem is nonlinear in nature and the exact MLE solution is difficult to compute.For this reason,the problem was initially investigated for the case when the sensors are arranged in a linear fashion1,2,8where different optimum processing techniques are proposed under various assumptions. When the sensors are distributed irregularly,the optimum solutions are much harder to obtain.A linearization approach based on Taylor series is proposed in5,15to compute the optimum solution iteratively.Because of the existence of local minima,the Taylor series method requires an initial solution close to the global minimum.A more appealing approach is the quasi-linear method employed in6,12,13that converts the nonlinear measurement equations into quasi-linear measurement equations by treating the nonlinear term as a parameter.Two different procedures are developed independently in6,13that turn out be mathematically equivalent.A further progress is made in3,10which treats the nonlinear term as an independent variable to be estimated that is solved via WLS procedure together with the location variables.A second WLS is then introduced to minimize the mismatch between the nonlinear term and the location estimate.This two-step WLS procedure is shown to be more effective and closer to the true optimum solution than other existing approaches.In this paper we consider the same passive source localization problem using the same quasi-linear measure-ment equations.Different from the existing solution methods,the nonlinear term in the quasi-linear measurement ∗This research is supported in part by US Air Force.Further author information E-mail:ggu@,Telephone:(225)578-5534Wireless Sensing and Processing II, edited by Raghuveer M. Rao, Sohail A. Dianat, Michael D. Zoltowski,Proc. of SPIE Vol. 6577, 657707, (2007) · 0277-786X/07/$18 · doi: 10.1117/12.719865equations is treated as a constraint which is in fact a quadratic constraint.A numerical procedure based on simultaneous diagonalization9is developed to compute the WLS solution under the quadratic constraint.This new procedure is shown to be an approximate MLE solution under the assumption that the ratios of the distances between all but one sensor and the target to the noise variance in TDOA measurements are suitably large.The assumption is rather mild and is satisfied so long as the relative distances among all the sensors are suitably large and the noise variance is suitably small which hold in practice for localization based on TDOA measurements. Our results are fairly complete and complement the existing ones in the literature.The simulation results il-lustrate the effectiveness of the proposed solution procedures that compare favorably with the existing results. However due to the space limite,the simulation results are not presented which are available in.72.PROBLEM FORMULATION AND PRELIMINARY ANALYSISFor simplicity we consider only two-dimensional localization although the results are applicable to the three-dimensional case as well.Denote(x k,y k)as the position of the k th sensor and(x T,y T)as the position of the stationary target.Our goal is to estimate the position(x T,y T)based on the TDOA measurement data collected by n sensors located at{(x k,y k)}n k=1.LetR k,T=(x T−x k)2+(y T−y k)2,R T=R1,T,(1)be the distance between the target and the k th sensor.Then the TDOAs are defined by∆t k,1:=(R k,T−R T)/c=T−x k2T−y k2−T−x12T−y12/c(2)for k=2,3,···,n with c the speed of light.The target location is embedded in TDOAs or∆t k,1=∆t k,1(x T,y T). As in the literature we assume that the measurements of TDOAs are given by∆ˆt k,1=∆t k,1(x T,y T)+δt k,1(3)where{δt k,1}n k=2are uncorrelated Gaussian random variables with mean zero and varianceσ2t.That is,the measurement noises have the joint probability density function(PDF)asp∆(δt;x T,y T)=1(2π)n−1σn−1texp−nk=212σ2t∆ˆt k,1−∆t k,1(x T,y T)2.(4)It follows that source localization based on TDOAs is a nonlinear estimation problem.For the joint PDF in(4),an MLE solution to the source localization problem seeks(x T,y T)such that p(δt;x T,y T)achieves its global maximum that is an unbiased estimate under some regularity condition.14The well-known Cram´e r-Rao bound dictates that the error covariance associated with any unbiased estimate is bounded below by the inverse of the corresponding FIM(Fisher information matrix).The computation of the FIM can be complex but has simpler formulas in the case of Gaussian random variables.For source localization based on TDOA measurements,the corresponding FIM can be obtained by applying the Slepian formula14that leads toFIM(∆t)=nk=21c2σ2tcos(βk)−cos(β1)sin(βk)−sin(β1)cos(βk)−cos(β1)sin(βk)−sin(β1)(5)whereβk is the bearing parameter given byβk=tan−1[(y T−y k)/(x T−x k)]and1≤k≤n.Denoteρk=x2k+y2k and rewrite(2)equivalently as(x T−x k)2+(y T−y k)2=(x T−x1)2+(y T−y1)2+c∆t k,12.(6)Then with appropriate scaling we arrive at the following equation:1√cρ2k−ρ21=2√cx k−x1y k−y1xTy T+c√c∆t2k,1+2√cR T∆t k,1(7)for2≤k≤n.The TDOA data in(3)can be substituted in to obtain a quasi-linear model as in.3,6,12,13Indeeddenotea k=1√cρ2k−ρ21−c2(∆ˆt2k,1+µσ2t),ηk=−2√c(R T+c∆ˆt k,1)δt k+c√c(δt2k−µσ2t)(8)for2≤k≤n whereµ≥0.Letθ=x T y T R Tbe the location parameter vector.Thena=Hθ+η⇐⇒⎡⎢⎢⎢⎢⎣a2a3...a n⎤⎥⎥⎥⎥⎦=2√c⎡⎢⎢⎢⎢⎣x2−x1y2−y1c∆ˆt2,1x3−x1y3−y1c∆ˆt3,1.........x n−x1y n−y1c∆ˆt n,1⎤⎥⎥⎥⎥⎦⎡⎢⎣x Ty TR T⎤⎥⎦+⎡⎢⎢⎢⎢⎣η2η3...ηn⎤⎥⎥⎥⎥⎦(9)The above is the same as the quasi-linear model in the literature with{a k}pseudo-measurements and{ηk}the corresponding noises.However it imposes a constraint to the parameter vectorθas(θ−θ0) Q(θ−θ0)=0,Q=diag(1,1,−1),(10)whereθ0=x1y10.Denote E{·}and E{·|·}as the expectation and conditional expectation.Forµ=1,E{ηk|∆ˆt k,1}=0,E{η2k|∆ˆt k,1}=4c(R T+c∆ˆt k,1)2σ2t+2c3σ4t(11) where E{δt3k}=0and E{δt4k}=3σ4t are used.The question is under what condition{ηk}are Gauss distributed, assuming that{δt k}are Gauss distributed with zero mean.The next lemma is useful.Lemma2.1.Let Z=−2αV+V2−σ2whereα>0and Z is a Gauss random variable with mean zero and varianceσ2.If the ratio ofαtoσis sufficiently large,then Z is approximately Gauss distributed in the neighborhood of−σ2with mean zero and variance4α2σ2.Proof:It is easy to verify that E{Z}=0andσ2Z=E{Z2}=4α2σ2+2σ4=2σ2(2α2+σ2).(12) By noting that Z=(V−α)2−(σ2+α2),we haveV=−α±22Z≥−(σ2+α2).(13)Since V is Gauss distributed,the PDF of Z is given byp Z(z)=1√2πσ2⎧⎨⎩e−12σ2α−√z+(σ2+α2)22z+(σ2+α2)+e−12σ2α+√z+(σ2+α2)22z+(σ2+α2)⎫⎬⎭.(14)Figure1.PDF p Z(z)versus Gauss PDF with x-axis scaled byσ−1,y-axis byσ,andα=10σLet U=22The integral of thefirst term in p Z(z)isI1=1√2∞−(σ2+α2)e−12σ2α−√z+(σ2+α2)2222dz=1√∞−ασe−v2/2dv=:Qασ.Hence it is concluded that ifασis sufficiently large,then I1≈1and thus p Z(z)is dominated by thefirst term.It is also easy to see that whenασis sufficiently large,σ2Z is dominated by4α2σ2.In the case of z+σ2= withclose to zero,we havep Z(z)≈1√2⎧⎨⎩e−12σ2α−√z+(σ2+α2)2222⎫⎬⎭≈1√22exp−(z+σ2)28σα(15)concluding that Z is approximately Gauss distributed in the neighborhood of z=−σ2.Although Lemma2.1states that Z is approximately Gauss distributed near−σ2,it does not imply that its PDF peaks at Z=−σ2.In fact it can be shown by straightforward calculation that the PDF of Z peaks at −3σ2,E{Z}≈0,and E{Z2}≈(2ασ)2provided thatασ>>1.Figure1shows the approximate PDF of Z given byp Z(z)≈1√2⎧⎨⎩e−12α−√ 2222⎫⎬⎭(16)in solid line and the PDF of Gauss random variable with mean−σ2and variance(2ασ)2in dashed line where α=10σis used.Therefore Z is approximately Gaussian not only in the neighborhood of Z=−σ2but in a much greater interval.We comment that the approximate PDF of Z and the exact of PDF of Z are indistinguishable for the caseα≥10σ.We note that Z is Gauss distributed if and only ifγZ is Gauss distributed withγa nonzero constant.In addition it is reasonable to assume that∆ˆt k,1≈∆t k,1for2≤k≤n.Hence by applying Lemma2.1to theelements of the noise vector in(9)we can conclude that{ηk}are approximately Gauss distributed provided thatασ=R T+c∆ˆt k,1cσt≈R T+c∆t k,1cσt=R k,Tcσt>>1.(17)Such a condition holds as long as the relative distances among all the sensors are suitably large and the noise variance is suitably small.It is thus concluded that by takingµ=0,{ηk}in(9)are approximately Gauss distributed and can be well approximated by Gauss variable with mean zero and variance4c(R k,T+c∆ˆt k,1)2σ2t≈4cR2k,Tσ2t provided that the condition in(17)holds.In factµ=0is what employed in.3,10 Because{ηk}are approximately Gauss distributed with mean zero and variance{4cR2k,Tσ2t}in the caseµ=0, the joint PDF in(4)can be replaced by the following PDFp Z(η;x T,y T)≈1(2π)n−1det(Σ)exp−12(a−Hθ) Σ−1(a−Hθ)(18)whereΣ≈4cσ2t diagR22,T,R23,T,···,R2n,T.Recall the constraint(10)for the parameter vector.Hence anapproximate MLE for the location parameter vectorθis the one that minimizesJ=12(a−Hθ) Σ−1(a−Hθ)(19)subject to the constraint(10)under the condition in(17).A solution will be provided in the next section to such constrained WLS problems.In the following we provide an approximate expression of the FIM for the quasi-linear model(9)when{ηk}are approximately Gaussin distributed.Lemma2.2.Suppose that the noise vector in the quasi-linear model(9)admits an approximate joint PDF in (18).Then the associated FIM with respect to the location parameter(x T,y T)is approximately the same as in(5).Proof:DenoteθT=x T y T.Then Hθ=H1θT+bR T whereH1=2√c⎡⎢⎢⎢⎢⎣x2−x1y2−y1x3−x1y3−y1......x n−x1y n−y1⎤⎥⎥⎥⎥⎦,b=2√c⎡⎢⎢⎢⎢⎣∆ˆt2,1∆ˆt3,1...∆ˆt n,1⎤⎥⎥⎥⎥⎦.(20)Thus the quadratic function in(19)that is the exponent of the joint PDF in(18)has the formJ=12(a−bR T−H1θT) Σ−1(a−bR T−H1θT).(21)Because both Hθ=bR T+H1θT andΣin the joint PDF involve the location parameters(x T,y T),the Slepian-Bangs formula14can be applied to compute the associated FIM that is given byFIM=nk=2c24R2k,Tσ2t2x kc2+b k cos(β1)2y kc2+b k sin(β1)2x kc2+b k cos(β1)2y kc2+b k sin(β1)+nk=22R2k,Tcos(β1)sin(β1)cos(β1)sin(β1).(22)Under the condition(17),the above is dominated by thefirst summation leading toFIM≈nk=214c2σ2t R2k,T2x k+c2b k cos(β1)2y k+c2b k sin(β1)2x k+c2b k cos(β1)2y k+c2b k sin(β1)=nk=21c2σ2tcos(βk)−cos(β1)sin(βk)−sin(β1)cos(βk)−cos(β1)sin(βk)−sin(β1)(23)that is identical to the expression in(5).We comment that the second summation in(22)is resulted from the Bangs’formula.Under the condition(17), it is insignificant compared with thefirst summation implying thatFIM is insensitive to(x T,y T)embedded in the covarianceΣ.In addition Lemma2.2shows thatFIM associated with the quasi-linear model is approximately the same as the FIM associated with the nonlinear measurement model.This fact justifies the formulation of the constrained WLS for computing an approximate MLE solution to source localization based on TDOA measurements.3.SOLUTION TO THE CONSTRAINED WLS PROBLEMMinimization of J in(19)subject to the constraint(10)has been tackled in6,12,13and in.3,10Several solution procedures have been developed but they all bypass the difficulty in solving the constrained WLS problem directly.In this section,we provide a general solution tomin θ Qθ=012(a−Hθ) Σ−1(a−Hθ)(24)The solution to the above constrained WLS is applicable to the case when the constraint is given by(θ−θ0) Q(θ−θ0)=0withθ0=0.See Remark3.2.In(24),the matrix Q is an indefinite nonzero matrix with size × and matrix H has rank that is its number of columns.We employ the method of Lagrange multiplier to solve(24).Letλbe real and considerJ=12(Hθ−a) Σ−1(Hθ−a)+θ Qθ.(25)Then the necessary condition for optimality yields the conditionH Σ−1[Hθ−a]+λQθ=0⇐⇒θ=[H Σ−1H+λQ]−1H Σ−1a.(26) An optimal solution needs to satisfy the constraintθ Qθ=0leading toa Σ−1H[H Σ−1H+λQ]−1Q[H Σ−1H+λQ]−1H Σ−1a=0.(27)The solution algorithm hinges to the computation of the real rootλfrom the above equation and there can be more than one such real root.The result of simultaneous diagonalization in9can be employed for this purpose. BecauseΣ=Σ >0and Q=Q ,there exists a nonsingular matrix S such that H Σ−1H=SDΣS and Q=SD Q S where DΣand D Q are both diagonal.It is noted that DΣand D Q have the same inertia asΣand Q,respectively.It follows that(27)is equivalent to(S−1H Σ−1a) (λI+DΣD−1Q )−1D−1Q(λI+DΣD−1Q)−1(S−1H Σ−1a)=0.(28)Let D −1Q =diag(q 1,q 2,···,q ).Then it has the same number of negative and positive elements as D =D ΣD −1Q =diag(d 1,d 2,···,d )by the positivity of Σand D Σ.In fact q i d i >0.The matrices S and D can be obtained by eigenvalue decomposition of A Σ−1AQ −1=SDS −1.Let v i be the i th element of S −1H Σ−1a .Then (28)is equivalently converted into the following:(S −1H Σ−1a ) (λI +D −1Q D Σ)−1D −1Q (λI +D ΣD −1Q )−1(S −1H Σ−1a )=i =1q i v 2i (λ+d i )2=0.(29)Because Q is indefinite,there is at least one strictly positive and one strictly negative element from {q i } i =1.Hence the above equation has at least one real root.However there are only finitely many real λvalues satisfying(29).In fact all the roots of the nonlinear equation in (29)are roots of the following polynomial with degree 2( −1):i =1q i v 2i k =i(λ+d k )2=0.(30)For source localization based on TDOAs,the constraint is given by (10)in which case =3implying that the computational complexity due to the roots computation is rather modest.Denote {λk }r k =1as the r real roots of (30).Now by (26),Hθ−a =[H (H Σ−1H +λk Q )−1H Σ−1−I ]a = HQ −1(λk I +H Σ−1HQ −1)−1H Σ−1−I a = (λk I +HQ −1H Σ−1)−1HQ −1H Σ−1−I a=−λk (λk I +HQ −1H Σ−1)−1a =−λk Σ(λk Σ+HQ −1H )−1aSubstituting the above into the performance index leads toJ =12λ2k a (λk Σ+HQ −1H )−1Σ(λk Σ+HQ −1H )−1a.(31)Let λopt be one of the r real roots that minimizes J over {λk }r k =1.Then in light of (26),the optimal θis obtained as θ=θopt given by θopt =[H Σ−1H +λopt Q ]−1H Σ−1a.(32)The above solution procedure is summarized into the following theorem.Theorem 3.1.Consider the constrained WLS in (24).Let H Σ−1HQ −1=SDS −1be the eigenvalue decompo-sition with {d i }eigenvalues.Then there holdH Σ−1H =SD ΣS ,Q =SD Q S ,(33)with D Σand D Q diagonal and D =D ΣD −1Q.Denote v =S −1H Σ−1a and {λk }as real roots of (30).If λopt minimizes J in (31)over {λk },then the solution to the constrained WLS is given by (32).Remark 3.2.If the quadratic constraint is δθ Qδθ=(θ−θ0) Q (θ−θ0)=0which involves θ0=0,then by Hθ=Hθ0+Hδθ,we may consider to solve alternativelymin δθ Qδθ=012[(a −Hθ0)−Hδθ] Σ−1[(a −Hθ0)−Hδθ].The above is identical to the constrained WLS in (24)with δθas unknown to which the solution procedure developed earlier is applicable.After the optimal δθopt is available,θopt =δθopt +θ0can also be obtained.It is noted that in applying the solution procedure in Theorem3.1to source localization based on TDOAs, the noise covariance matrixΣis unavailable that depends on the target location(x T,y T).As in the literature, we may setΣ=I to compute an initial location solution(ˆx(0)T,ˆy(0)T)via the constrained WLS procedure which can then be substituted intoΣto compute the solutionθopt to the corresponding constrained WLS.Although further iterations can be employed to obtain a more accurate solution,simulation results show that they are unnecessary.In fact by the proof of Lemma2.2,the second summation in the expression ofFIM in(23)is due to the Bang’s formula obtained by computing derivative ofΣwith respect to(x T,y T).This term is insignificant compared with thefirst summation implying that the minimum achievable error variance by unbiased estimators is insensitive to the location parameter(x T,y T)embedded inΣ.Therefore it is adequate to compute the solution to the constrained WLS twice.Finally we would like to comment that the constrained WLS solution is an approximate MLE solution to the source localization problem based on TDOAs by treating measured TDOAs{∆ˆt k,1}n k=2as deterministic quantities.That is,it is an approximate MLE solution conditioned on the measured TDOAs.However it is also important to observe that even if we treat{∆ˆt k,1}n k=2as random,the constrained WLS solution is still an approximate MLE solution.Recall that a WLS solution can be obtained by choosing a perturbation vectorˆ∆a that has the smallest Euclidean norm such that16rankΣ−1/2HΣ−1/2a+ˆ∆a=rankΣ−1/2H(34)and then solveΣ−1/2Hθ=Σ−1/2a+ˆ∆a.Basically the imperfect measurements in a are removed byˆ∆a.It is noted that the third column of H,or b in(20),involves imperfect measurements as well.However the variance of the measurement noise in vector b is only4cσ2t,much smaller than4cσ2t R2k,T(variances ofηk)under the condition(17).Hence it is unnecessary to introduce the perturbation to the third column of H that does not affect the approximate MLE estimate that is validated in our simulation study,although not reported here.4.CONCLUDING REMARKSThe problem of passive localization has been investigated based on TDOA measurements which is a nonlinear estimation problem.It has been shown that the nonlinear term in the quasi-linear measurement equations can be treated as a constraint which is in fact a quadratic constraint.A numerical procedure based on simultaneous diagonalization is developed to compute the weighted LS(LS)solution under the quadratic constraint.This new procedure is shown to be an approximate MLE solution under the assumption that the ratios of the distances between all but one sensor and the target to the noise variance in TDOA measurements are suitably large.Because of the space limit,the simulation results are not presented in the present conference version.As concluding remarks,we comment that the assumption on uncorrelated TDOA measurement noises can be removed.Such an assumption is in fact not used in deriving the constrained WLS solution.That is,the solution procedure summarized in Theorem3.1is applicable to the case whenΣis an arbitrary symmetric positive definite matrix. See also the results in3for the derivation of the various covariance matrices when the TDOA measurement noises are correlated.We also comment that the moving source localization studied in10involves two quadratic constraints to which a similar constrained WLS solution can be derived.However it involves roots computation for two2-variate polynomials which has high complexity.How to bypass such a high complexity is currently under study.REFERENCES1.J.S.Abel and J.O.Smith,“Source range and depth estimation from multipath range difference measure-ments,”IEEE Trans.Acoust.Speech,Signal Processing,vol.37,pp.1157-1165,Aug.1989.2.G.C.Carter,“Time delay estimation for passive sonar signal processing,”IEEE Trans.Acoust.Speech,Signal Processing,vol.29,pp.463-470,June1981.3.Y.T.Chan and K.C.Ho,“A simple and efficient estimator for hyperbolic location,”IEEE Trans.SignalProcessing,vol.42,pp.1905-1915,Aug.1994.4.B.T.Fang,“Simple solutions for hyperbolic and related positionfixes,”IEEE Trans.Aerosp.Electron.Syst.,vol.26,pp.748-753,Sept.1990.5.W.H.Foy,“Position-location solutions by Taylor-series estimation,”IEEE Trans.Aerosp.Electron.Syst.,vol.12,pp.187-194,Mar.1976.6.B.Friedlander,“A passive localization algorithm and its accuracy analysis,”IEEE J.Ocean.Eng.,vol.12,pp.234-245,Jan.1987.7.G.Gu,“Target localization based on TDOAs and FDOAs,”Quarterly report to Sensor Directorate,WPAFB,March2006.8.W.R.Hahn and S.A.Tretter,“Optimum processing for delay-vector estimation in passive signal arrays,”IEEE rm.Theory,vol.19,pp.608-614,Sept.1973.9.R.A.Horn and C.R.Johnson,Matrix Analysis,Cambridge University Press,reprinted in1999.10.K.C.Ho and W.Xu,“An accurate algebraic solution for moving source location using TDOA and FDOAmeasurements,”IEEE Trans.Signal Processing,vol.52,pp.2453-2463,Sept.2004.11.A.H.Sayed,A.Tarighat,and N.Khajehnouri,“Network-based wireless location,”IEEE Signal ProcessingMagazine,vol.22,pp.24-40,July2005.12.H.C.Schau and A.Z.Robinson,“Passive source localization employing intersecting spherical surfaces fromtime-of-arrival differences,”IEEE Trans.Acoust.Speech,Signal Processing,vol.35,pp.1223-1225,Aug.1987.13.J.O.Smith and J.S.Abel,“Closed-form least-squares source location estimation from range-difference mea-surements,”IEEE Trans.Acoust.Speech,Signal Processing,vol.35,pp.1661-1669,Dec.1987.14.P.Stoica and R.Moses,Introduction to Spectral Analysis,Prentice-Hall,Upper-Saddle River,NJ,1997.15.D.J.Torrieri,“Statistical theory of passive location systems,”IEEE Trans.Aerosp.Elect.Syst.,vol.20,pp.183-198,Mar.1984.16.S.van Huffel and J.Vandewalle,The Total Least Squares problem:Computational Aspects and Analysis,SIAM Publisher,1991.。