Chinese Journal of Aeronautics, (2014),27(2): 365-374

JOURNAL OF

AERONAUTICS

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics

cja@buaa.edu.cn www.sciencedirect.com

A closed-form solution for moving source localization using LBI changing rate of phase difference only

Zhang Min, Guo Fucheng *, Zhou Yiyu

College of Electronic Science and Engineering, National University of Defense Technology, Changsha 410073, China

Received 12 July 2013; revised 6 September 2013; accepted 28 November 2013 Available online 28 February 2014

KEYWORDS

Changing rate of phase difference;

Cramer-Rao lower bound; Estimation bias; Instrumental variable; Least square;

Long baseline interferometer;

Radiated source localization

Abstract Due to the deficiencies in the conventional multiple-receiver localization systems based on direction of arrival (DOA) such as system complexity of interferometer or array and amplitude/phase unbalance between multiple receiving channels and constraint on antenna configuration, a new radiated source localization method using the changing rate of phase difference (CRPD) measured by a long baseline interferometer (LBI) only is studied. To solve the strictly nonlinear problem, a two-stage closed-form solution is proposed. In the first stage, the DOA and its changing rate are estimated from the CRPD of each observer by the pseudolinear least square (PLS) method, and then in the second stage, the source position and velocity are found by another PLS minimization. The bias of the algorithm caused by the correlation between the measurement matrix and the noise in the second stage is analyzed. To reduce this bias, an instrumental variable (IV) method is derived. A weighted IV estimator is given in order to reduce the estimation variance. The proposed method does not need any initial guess and the computation is small. The Cramer-Rao lower bound (CRLB) and mean square error (MSE) are also analyzed. Simulation results show that the proposed method can be close to the CRLB with moderate Gaussian measurement noise.

© 2014 Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA.

1. Introduction

Determining the position and velocity of a moving radiated source at a single or multiple observers, which is also referred as radiated source localization or passive localization, has been

widely used in many areas such as radar, sonar, reconnaissance, and wireless network.1-4 A widely applied technique is the triangulation using the direction of arrival (DOA)5-7 measured at multiple observers. It does not need accurate synchronization and high-speed data links between observers.8 However, the emitter velocity cannot be directly deduced jointly with the source position form DOA measurements in one measurement instant. The target dynamic state should be jointly estimated using the DOA and its changing rate9 measurement obtained in one instant at multiple observers or by DOA-only measurements cumulated in multiple instants at a single or multiple observers.10,11

In addition, to measure the DOA precisely, spatial spectrum-based direction finding techniques such as MUSIC (multiple signal classification)12 and ESPRIT (estimation of

* Corresponding author. Tel.: +86 731 84573490. E-mail addresses: zhangmin1984@126.com (M. Zhang), gfcly@ 21cn.com (F. Guo), zhouyiyu@sohu.com (Y. Zhou). Peer review under responsibility of Editorial Committee of CJA.

1000-9361 © 2014 Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. http://dx.doi.Org/10.1016/j.cja.2014.02.013

signal parameter via rotation invariant subspace)13 methods can get better precision and resolving performance. However, these methods have large computation and system complexity, due to e.g., signal covariance matrix estimation and a large number of antennas for achieving good performance.14 Moreover, signal and noise should be modeled accurately, which leads to the sensitivity to the model imperfections.

An alternative direction finding technique uses a phase interferometer.14-16 It needs multiple baselines to eliminate the 2p ambiguity16 of the phase difference (PD) measurements when the baseline length of the interferometer is longer than half of the received signal wavelength, which induces hardware complexity and rigorous constraint on the geometric configuration of the antenna array. Furthermore, the array or interferometer performance is susceptibly affected by the amplitude/ phase unbalance between receiving channels. Array calibration may effectively mitigate the impacts of the model imperfections and channel unbalance.17-19 However, it also causes an increase of system complexity and a higher cost.

To utilize a longer baseline for performance improving and meanwhile reduce the receiving number of antennas and channels, a new multiple-observer localization system using changing rate of phase difference (CRPD) only is studied in this paper. The position and velocity of the moving emitter are jointly determined from the CRPDs measured by a long baseline interferometer (LBI) array at multiple observers. Since the CRPD can be measured by estimating the slope of the PD sequence, the PD bias20 between receiving channels caused by amplitude/phase unbalance can be removed without extra calibration. Moreover, the CRPD estimation merely requires that the change between two consecutively measured PDs is less than p. It means that there is no need to solve the absolute ambiguity of each PD, and only the relative ambiguity between the measured PDs is adequate for the measuring. Consequently, a single LBI is sufficient and no more compound baselines are needed21 for ambiguity resolving.

Determining the source position and velocity from the CRPD measurements obtained in a single time instant is not a trivial task. This is because the source location is nonlinearly related to the CRPD measurements. To find a solution from the CRPDs, the most straightforward method is the exhaustive search22 in the solution space. This is robust but computationally expensive, which is hard for real-time applications. A possible alternative is to recast the localization problem as a nonlinear least square (NLS) problem by using the maximum likelihood (ML) solution. In general, iterative gradient search techniques, such as the Taylor series expansion23 and the Gauss-Newton (GN) algorithm6, are employed to find the numerical ML solution. The iteration method, however, requires a proper initial position and velocity guess close to the true solution, which may not be easy to satisfy in practice. Furthermore, under some circumstances when the localization geometry is poor, this method may fail to converge.

For DOA-only also called as bearing-only (BO) localization problems, the pseudolinear least square (PLS) methods have been an active research area for several decades.24-27 The pioneering work is the Stansfield5 estimator. Most of the current BO emitter localization algorithms have been developed from this algorithm, such as the instrumental variable (IV) method,24 the weighted instrumental variable (WIV) method,25 the total least square (TLS) method,26 the constrained total least square (CTLS) method,27 and so on.

However, these methods cannot be directly applied to the CRPD localization problem because of the high nonlinearity between the unknowns and the measurements.

To make use of the geometrical characteristic intrinsically in the CRPD measurements, this paper proposes a computationally attractive two-stage closed-form solution to determine the source position and velocity using the CRPD measurements. The proposed algorithm does not require an initial solution guess and it involves PLS minimization only. In the first stage, the DOA and its changing rate are estimated from the CRPD of each observer by PLS minimization. In the second stage of the algorithm, the source position and velocity are found by another PLS minimization.

Because the measurement matrix of PLS is generally correlated with the noise vector, the pseudolinear estimator is not consistent and therefore biased.24-28 To reduce the estimation bias caused by the correlation, an IV method using CRPD only is proposed. A weighted version is also developed to reduce the estimation variance of the IV estimator. The proposed method attains the Cramer-Rao lower bound (CRLB) at a moderate noise level when the distribution of the CRPD measurement noise is Gaussian.

The paper is organized as follows. The next section presents the multiple-observer localization model using CRPD only. Section 3 reviews the ML estimator and develops the proposed two-stage solution. Section 4 derives the bias of the two-stage algorithm. Section 5 gives the IV estimator and its weighted version. Section 6 derives the CRLB of the unknowns in the two stages. Section 7 gives simulations to support the theoretical developments of the proposed estimator and compares it with different localization methods. Finally, conclusions are given in Section 8.

2. Problem formulation

2.1. Localization model

Without loss of generality, a two-dimensional (2D) scenario in which the emitter is far from the observers is considered in this paper. In particular, N geographically separated observers on the ground with known positions xn — [ xn yn ]T and velocities vn — [ vxn vyn ]T, n = 1,2, ■■■, N, are used to determine the unknown position xt — [xt yt ]T and velocity vt — [ vxt vyt ]T of a moving source. This is illustrated in Fig. 1. The carrier frequency of the emitted signal is f, corresponding to the wavelength k — cjf of the emitter signal, where c is the speed

Fig. 1 Localization sketch using CRPD only.

of light. An LBI array is mounted at each observer. It is postulated that M long baselines with different baseline pointing vectors are formed by each LBI array and a CRPD can be measured from each long baseline. The location problem requires at least N = 2 observers and at least M = 2 long baselines for each LBI array.

Suppose that the LBI arrays are horizontally deployed. The unambiguous PD of baseline m at observer n is

0nm — jdnm cOs(ßn - 0m

where k — 2p/k, dnm is the baseline length, 0nm is the baseline azimuth, and jn is the DOA with respect to observer n defined

ßn — atan( X-— )

V^t- yn)

Let h — ^xT vT]T, and the corresponding CRPD20 of baseline m at observer n is

^nm fnm(h) jdnm sin(jn hnm)jjn ^ dnm (3)

where dnm is the zero mean white Gaussian measurement noise with covariance a^m. To simplify the analysis, the noise is assumed to be independent identity distribution (IID) with covariance r\m — r2, jjn is the DOA rate of changing with respect to observer n, which is the first-order derivative of jjn with respect to time and equal to

(yt - yn)(vxt - Vxn)-(xt - Xn)(Vyt - Vyn)

where rn = ||xt — xn\\2 is the distance between the emitter and observer n, in which ||-||2 denotes the Euclidean norm.

Stacking the measured CRPDs of all the N observers yields

z = f[0)+e (5)

where z — [ zj

/(h)— [/T f T

_ /n1 /n2

Zn2(0)

■ /nu\ ,

■ fnM(h)]T,

fn — [fm(0)

e — [eT eT ■■■ eTN]T, and en — [ 8nl dnl ■■■ dnM ]T.

We shall focus on estimation of the source position and velocity using the CRPD measurements only.

2.2. Measuring the PD

The CRPD is usually obtained by fitting the slope of the PD sequence. The PD can be measured using a two-element radio interferometer as shown in Fig. 2.

As shown in Fig. 2, antennas A1 and A2 receive the plane wave signals S1 and S2 and send the signals to the corresponding receiving channels. A receiving channel includes a low-noise amplifier (LNA), a down converter (D/C), an analog-to-digital converter (A/DC), and a fast Fourier

Channel2

Fig. 2 Phase difference observation model.

transformer (FFT), as well as a common reference oscillator (RO) for equipment-delay variations and a local oscillator (LO) for phase synchronizing.15

The signals in channels 1 and 2 are Fourier transformed, respectively to X1(x) and X2(x) over a known bandwidth of x. Their complex product is calculated as

Y(x) — Xl(a)X*2(a) (6)

For a narrow band signal, the PD is obtained as

/ — arg(F(wmax)) (7)

where arg is the argument of a complex number and wmax is the frequency at the maximum of the power spectrum | Y(ra)|.

For a wide band signal, the PD can be obtained by estimating the slope of arg(F(w)) around the band-center frequency.15 These techniques are referred to as the frequency-domain correlation known as the FX correlation processing.29

2.3. Measuring the CRPD

However, when the baseline length is longer than half of the signal wavelength, the estimated PD /m(t) at baseline m of observer n may appear 2p ambiguity. Mathematically, we have

CmW — /nm(t) - 2pk(t) + bnm + fnm 2 (0, 2p] (8)

where k(t) is an unknown ambiguous integer of each measured PD at t, bnm is the system bias caused by the channel unbalance, and fnm is random measurement noise.

It is necessary to eliminate the impact of k(t) before measuring the CRPD. Because the target speed is a limited value (for example less than 600 m/s), the change caused by the source movement between consecutively measured PDs would be less than p within a short interval (such as the radar pulse recurrence interval). This can be derived from Eqs. (2)-(4) as follows

I _ I ^ j I _ I jdnm vn i • : 0 \ i , jdnm vn \<Pnm\ 6 Kdnm\ßn\ — -| Sin(ßn - <*n)\ 6 -

where an — atan

vxt - vx

is the course of the relative veloc-

ity between the emitter and observer n, and vn = ||vt — vn||2 is the corresponding speed.

The change of the PD during interval At is

i _ i a . jdnm vn » . \ /nm \ Dt 6 -Dt

Taking f = 10 GHz, dnm = 50 m, vn = 300 m/s, rn

km, and At = 1/100 s as an example, p

lA/nml 6 p

As shown in Fig. 3, once an abrupt change exceeding p appears between two consecutively measured PDs, a 2p jump occurs15 and a modified PD sequence /m(t) can be obtained by the detected jumping as

/Um(t) — €m(t) + 2/(t)P — /mn(t)- 2koP

where l(t) is the relatively ambiguous integer of the measured PD with respect to the first PD and k0 is the unknown ambiguous integer of the first PD.

Although the phase ambiguity may be still present in sequence /m(t), it is a smooth function of time, in contrast to the original discontinuous sequence /m(t). As a result, the

where i = 0,1,..., ÔML(0) is an initial guess that must be sufficiently close the ML solution to ensure convergence, P(i)=Z - f(ÔмL(i)), and J(i)= /(fl)^^^ in which the expression of J(h) is given in Eq. (72).

The GN algorithm is stopped at the maximum iteration

(JT(is)K-1J(4))- JT(4)K-1p(4)\\ < e

Fig. 3 Phase difference unwrapping.

CRPD of sequence /mn(t) can be determined equivalently by estimating the slope of sequence /m(t)

/UmW = /mn (t)

The CRPD can be estimated using the least square or Kal-man filter method.21 It is worth to note that the PD bias can be removed by estimating the slope of the PD sequence.

3. Localization algorithm

3.1. Maximum likelihood estimator

The likelihood function for the CRPD measurements is given by the joint probability density function conditioned on the emitter location

p(z|fl)=-

(2p)NM/2 det (K)

f(e))TK-1(z - /(h))

where K = diag(K1, K2, ■■■,KN) is the NM • NM diagonal covariance matrix of the CRPD noise of N observers, in which Kn = diag(rn1; ■■■ ; rnM) is the M • M diagonal covariance [-Kdnm cos 6nm Kdnm sin 0nm ]T.

where e is a threshold.

The ML estimator is asymptotically unbiased and efficient. However, it requires a proper initial position and velocity guess close to the true solution, which may not be easy to satisfy in practice. Furthermore, under some circumstances when the localization geometry is poor, the iteration method may fail to converge. In addition, the computational complexity of the ML estimator is large than that of the linear LS estimators because of the iterative nature of the ML estimation process.

3.2. Two-stage pseudolinear least square estimator

A two-stage pseudolinear estimator can be derived from the CRPD measurements. The DOA and its changing rate from the CRPD of each observer are estimated during the first stage. The source position and velocity are found from the DOAs and their changing rates during the second stage.

3.2.1. Stage 1—DOA and its changing rate estimation The CRPD equation of baseline m of observer n can be expressed as

nm = jdnm ßn cos hnm sin ßn

+ Kdnmßn sin Onm cos ßn + &nm (19)

Let yn = [ßn sin ßn ßn cos ßj T and the pseudolinear relationship between the CRPD and yn becomes clear. We shall estimate cn via the PLS method.

Stacking the M CRPDs of observer n gives

Zn = AnCn + en (20)

where An =[ a„i an2 ■■■ anM ]T, in which anm =

matrix of the CRPD noise of observer n.

The maximum likelihood (ML) estimator of the emitter location hML is obtained from the maximization of the log-likelihood function lnp(z\6) with respect to h, which can be equiv-alently written as

hML = argminJML(h) (15)

where JML(h) is the ML cost function given by

Jml(A)=1 (z - f(0))TK-l(z - AO))

The minimization of JML(O) is an NLS problem without a closed-form solution. In general, iterative gradient search techniques, such as the GN algorithm, are employed to find the numerical ML solution. The GN algorithm is given by the recursion

The PLS solution of yn is

"in — (A„An) AnZn

The source DOA and its changing rate can be obtained by using yn via

ßn — atanf ^ Vrn(2,1)

(16) b n — ±^y2(1,1) + cn(2,1)

0ml(î + 1) — ôml(Î) + (JT(i)K-1J(i))- JT(i)K-1q(i)

where the sign ambiguity can be removed by yn(1,1) / sinpn or

y„(2,1)/cos bn.

3.2.2. Stage 2—Position and velocity estimation From Eqs. (2) and (4), the DOAs and their changing rates are nonlinearly dependent on the unknown position and velocity. However, a pseudolinear equation with respect to the

time imax or at iteration is if

unknowns can still be derived, which enables the application of the closed-form LS technique.

Specifically, Eq. (2) can be rewritten as, after some manipulations,

Xt cos ßn - yt sin ßn — bßn +

where bpn — xn cos bn — yn sin bn and e^n is the equation error due to the DOA estimation error.

Similarly, Eq. (4) can be transformed into

nXt sin ßn + ßnyt cos ßn - Vxt cos ßn + Vyt sin ß

— bh + eh

In + ßnyn cos ßn - Vxn cos ßn + Vyn sin ß

where b^ — b„x„ sin^ and e^ is the equation error due to the noise in the source DOA and its changing rate estimations.

Stacking Eqs. (24) and (25) for the available N observers yields

b — H6 + e

where H — [ HT Ht2 ■■■ HTN ]T, Hn — [hb„ ^ ]" hl — [ cos bn — sin bn 0 0] ,

h'jl„ — [ bn sin bn bn cos bn — cos bn sin b n ] ,

— [bßn bßn

]T, and Sen — [eßn eßJT.

b = [bTi bT2 ■■■ bl] , bhn e — [ eT1 eT2

The source position and velocity can be found via

ÖPLS — (HTH)—Vb (27)

The proposed method does not need any initial guess and the computation load is not heavy. However, as shown in the next section, the pseudolinear estimator is biased due to the correlation between H and e.

4. Bias analysis of the pseudolinear estimator

4.1. Bias of stage 1

The bias of estimator for yn in the stage 1 is given as

bias(y„) — E(y„) — yn — —E[(ATnAn)—1ATnen] (28)

Since the geometry of the LBI array at observer n is known, the matrix An is determinate. The bias of yn is therefore

bias (in) — (ATAn)- ATE(en) — 0

Under the assumption that the estimating noise of CRPDs is small, expand ßn at the true value yn according to the multivariable Taylor series up to the first order term

ßn « atanf+ f

Vrn(2; 1)J @in

••xtT-« ßn _ r I__\ /_> !

where —[(cosßn)/'ßn (- sinßn)/'ßn] is the gradient of ß„ with respect to yn and dy — I'd.

y — L°yn(1.1) °yn(2,1).

is the estimate

error of yn.

Similarly, for the unambiguous estimation b - n, we can get

ßn «^yn(1; 1)+cn(2; 1) + mdyn

where — [ sin ßn cos ßn ].

@ yn n n

Let gn = [bn j}„] . Under the above assumption as well as the conclusion in Eq. (29), the bias of gn is therefore given as

BT E(\

bias(gn) — E(gn

is the partial

(cos ßn)/ßn (- sin ßn)/ßn sin ßn cos ßn

derivative of gn with respect to yn.

It can be seen from the above analysis that the estimation of gn is approximatively unbiased under small CRPD noise.

4.2. Bias of stage 2

Generally, when the measurement matrix H is correlated with the noise vector s, the pseudolinear estimator is not consistent, that is, ePLS does not convergence to the true location vector e and is therefore biased.

The bias of the pseudolinear estimator of stage 2 is

bias(Öpls) — E(hpLs)- 6 — -E[(HTH) Hts]

In order to analyze the asymptotically biased property, enough large N is considered in the following. The pseudolinear estimation bias can be approximated by30

bias(6P

The bias is the correlation between the entries of the measurement matrix H and the least squares noise vector e. Because of the high nonlinearity of the problem, the bias analysis is more difficult than the BO localization or tracking estimator as in Ref. 25 or 28. In what follows, we provide a bias analysis for the proposed pseudolinear estimator.

Applying the Taylor series expansion up to the first order term at the true value bn transforms Eq. (24) into

ebn ~ rndbn (35)

where Spn is the DOA measurement error.

Similarly, for Eq. (25), after the first-order approximation at the true value gn, we get

Cndßn - fndßn

in + Vyt cos ßn - Vxn sin ßn - Vyn cos ß

where cn — —( vxt siny

For the measurement matrix, expand h^ at the true value gn according to the multivariable Taylor series up to the first order term

hß « +in dßn

dltßn _ r

hßn « hßn

- cos ßn 0 0 ]T.

Expand hb using the similar transformation at the true va-

lue gn

hßn « hßn + ^ dßn + d ß

where — cos bn — bn sin bn sin bn cos bjT and

dh■ °bn

-h — [ sin bn cos bn 0 0 ]T.

Substitute Eqs. (35)-(38) into Eq. (34) and get the terms of the noiseless measurements and the measurement noise

E(Ht) =-£e(hT)

dh= dh

«1*2 + u2r2

ßnßn 1 "3rÀ

dhßr dh

where «1 ß + Cn"^T' «2 ß C^-rfi - rn^, «3 ß-rn--^,

dßn dßn dßn dßn dßn

4 ß E(d2n )' rß n ß E(d2n )' rßnßn ß E(dßndßn )' E(dßn )-0, and

E(dß n )~ 0.

The autocorrelation matrix E(HTH/N) can be written as

HTH\ 1 N

(V1 rßn + ^ n + V3 r

@hbn @hbn T @hbß @hbß T

where V—SIWJ + @fXwJ, V—2,

and V3 — ^ (dJ%X T

dp A d_n

Substituting Eqs. (39) and (40) into Eq. (34) gives the estimation bias for the pseudolinear estimator, in terms of the geometric parameters and measurement noise variance.

It can be seen from the above analysis that the bias is dependent not only on the correlation vector E(HTe/N), but also on the verse of the autocorrelation matrix E(HTH/N). It can be observed that the bias can be alleviated by reducing the CRPD noise or by using well-conditioned observer geometries for good observability.

5. Weighted instrumental variable estimator

5.1. Instrumental variable estimator

As shown in the above analysis, the second-stage pseudolinear estimator is biased. The method of instrumental variables (IV) can be used to reduce the bias.24'25

According to the instrumental variable method, the normal equations are given as

■dh'ßnf dhß,

GTH0W = GTb

where G is the instrumental variable matrix. The corresponding solution is

hIV — (GTH)—1GTb

Similarly, for enough large N, the estimation bias can be approximated by

bias(hi i

GTH\ 1 (GTe

Matrix G should be chosen in such a way that E(GTH/N) is nonsingular and E(GTe/N) = 0, and then the IV estimate hIV will be asymptotically unbiased.

The selection of G is problem-dependent,31 and it is not easy to formulate a proper matrix G. In the following, a suboptimal closed-form IV estimator using the pseudolinear estimate is proposed.

The optimal instrumental variable matrix is given by the noise-free version of the H matrix. However, it is not available in practice. In the proposed IV estimator, the DOA and its

changing rate for each observer are calculated using the pseudolinear estimate hPLS = [ xt yt vxt vyt ]T as follows

~n = atan

i Xt — Xn

V yt— y„

(yt — yn)(i)xt — Vxn) — (Xt — Xn)(Vyt — Vyn)

where rn — ||it — x„||2.

The instrumental variable matrix can be constructed as

G — G1T G2T where Gn —

and S_„ — [~n sin cos — cos sin ~n] .

Since the DOA and its changing rate are obtained from rPLS, matrix G is asymptotically uncorrelated with e, i.e.,

?Nf (46)

T [ ~ ~ T

, in which gßn = ^cos ~n — sin ~,00]

For the observable localization scenarios, E(HTH/N) is nonsingular and hPLS can be calculated. It means that E(GTH/N) will be nonsingular. From the above analysis, it is seen that hIV is asymptotically unbiased for the instrumental variable given in Eq. (46).

5.2. WIV estimator

To reduce the estimation variance of the IV estimator, a weighting matrix is introduced into the proposed IV estimator.

In the first stage, a weighting matrix Wn can be chosen to minimize the error covariance of yn, i.e., E[e„eJ] 1.32 Assuming that the CRPD noise variance r^m is known, the weighting matrix Wn can be obtained from

Wn — E[e„eT]—1 — diag(r-2, r—2, ■■■, r—M) (48)

The WIV solution of yn is

yn — (AJnWnAn)—1 AJnWnZn (49)

In the second stage, the weighted matrix W can also be chosen to minimize the error covariance of h, i.e., W = E[eeT|h]—1. For this purpose, we evaluate W as follows. From Eqs. (35) and (36), we have

' Dnegn

and e„ — \6ß

egn L °bn 0'Hn\ .

Collecting the results for all the N observers yields e — Pe„

where eg = [ eT1

P = diag(Dx,

The weighting matrix is therefore W = E[ggT]—1 = (PQPT)—1 = P—TQ—1P—1

(51) and

where Q = diag(Q1, Q2, ■■■, QN), in which the expression of Qn is given in Eq. (58).

The source position and velocity can be found via

hWIV = (GTWH) 1GTWb

The evaluation of the weighting matrix requires the knowledge on the true emitter position and velocity, which are unknown. To bypass this difficulty, W is first obtained from 0pls and then the computation of W is repeated using 0wiv to improve the estimation accuracy.

5.3. MSE of the WIV estimator

The corresponding covariance matrix of hwiv using the weighting matrix can be computed as

MSE(hwiv) - CRLB(h)

MSE(0wiv) - ( HT WH) 1 = ( HTP TQ-1P-1H)-

where H « G for the moderate noise level.

Note that MSE( 0WIV) has the same form as the CRLB in Eq. (71). Their relationship is derived as follows.

In the first stage of the proposed algorithm, the error covariance matrix of yn is given as

MSE (y,,) — r2 (ATA,)"1 (55)

The estimation error covariance matrix of gn therefore is MSE (gn) — BnMSE (y,)BTn (56) The inverse matrix of Bn can be simply evaluated as

fin cos ^ sin L n sin Pn cos Pn Substituting Eq. (57) into Eq. (56) yields

Qn — MSE (,n) = r2 (B;TATAn B"1)"1 (58)

It is easy to get the following result according to the matrix multiplication

AW-1- f

AnBn — ...

din @f

where -1- is defined in Eq. (69).

The inverse matrix of Q in Eq. (52) is therefore Q"1 — diag (Q"1, Q"1,Q"1) —r"2 f f

In the second stage of the proposed algorithm, the inverse matrix of D„ is

.cn/rl -1/r

P-1 = diag ( 0nxn, 0nxn, ■■■, D-1,..., On

It is easy to get the following result according to the matrix multiplication

P- H —

After some straightforward algebraic manipulations, the result is

where CRLB(h) is given in Eq. (73).

As a result, it shows that the estimation accuracy of the proposed method attains the CRLB under a Gaussian noise assumption when using the proposed weighted matrix.

6. Cramer-Rao lower bound analysis

6.1. CRLB of the DOA and its changing rate

The CRLB gives the lowest possible variance that an unbiased estimator can achieve. The CRLB for the DOA and its changing rate in the first stage is firstly analyzed. The likelihood function of the CRPD measurements given by the joint probability density function conditioned on the emitter location for observer n is

P ( Zn\gn)--

( 2p)M/n det (Kn)

x exp 2 ( Zn " fn( 1n))TK"1 ( z " fn( in))) ( 66)

The CRLB is equal to the inverse of the Fisher matrix32 defined as

F(gn) - E

idln p ( Zn\gn)\Y dln P ( Zn\gn

V dgn / V dgn

The CRLB of gn is therefore

CRLB (gn) - F-1 ( gn) - rn( JT(gn)J(gn))-1 where the Jacobi matrix J(g„) is equal to

n \ dfn

J( gn)-f

d4 dJk

dgn dgn

df.M dgn

dfnm dfnm dfnm

, dgn Wn

(68) (69)

dfnm ■ dfnm

— jdnm cos(Pn " 0nm)Pn, and —P — jdnm sin(Pn " ®nm). @Pn d_n

6.2. CRLB of the position and velocity

The CRLB of the source position and velocity is defined as

m — Ejp^p^j (70)

The CRLB of h is therefore

CRLB (h) — F"1 ( h) — r2 (J ( h)J( h))"1 ( 71)

where J(h) is the Jacobi matrix, which, according to the chain rule, is given as

From Eqs. (59) and (64), a conclusion can be drawn as

J( dfjdf dg ( ) dh dg dh

The CRLB can be computed as

CRLB(h)-n( f If f g)"

<r Moving emitter ir -

A Stationary observers Target

Al ■

■ A3 A2 ,

-40 -20 0 20 40 60 80 X (km)

Fig. 4 Simulated emitter localization geometry.

The partial derivates — and — are given as follows. In particular, we have g

df—diagff ■■■ f

dg V@g1 ,dg2 ; ,dgN

The partial derivative of g with respect to h is

dgT dgT

where ^ — dh

dj. dj, dj\

dxt dyt dvxt

dK djjn d^

dxt dxt d vxt

dk dvyt d jn

d jn = yt — yn

dvxt rl ,

.dh dv„

dyt dvyt

~ _ xt — xn —2 jn

dfjn — djn — 0 din —

dvvt '

, d jjn

, and —— —

vxt — v

dxt yt — yn

7. Simulation

The emitter localization geometry employed in the simulation is shown in Fig. 4. A scenario with three stationary observers on the ground is considered. The positions of the observers are xi — [0 22.5 ]T km, x2 — [22.5 0 ]T km and x3 — [—22.5 0]T km, respectively. Each observer has an LBI array including three antennas fixed on the vertex of an equilateral triangle with a side length of 45 m. The position and velocity of the moving emitter is xt — [ 70.7 70.7 ]T km and vt — [ 100 330 ]T m/s. The carrier frequency of the emitter is 10 GHz and the pulse repetition frequency (PRF) is 1000 Hz. The observation time for measurement is 0.5 s.

7.1. Performance of the DOA and its changing rate estimation

Firstly, the estimation performance of the DOA and its changing rate in stage 1 is simulated. The bias of the estimation is defined by the distance between the mean estimate and the true value, i.e.,

L —") (L —u

where ul is an estimate of u in simulation l and L is the total simulation number.

The root mean square error (RMSE) of the estimator is defined by

RMSE =

L XX (»■

- u)T(Ul-

The simulation results for the bias of the DOA and its changing rate at observer 1 as a function of the CRPD noise standard deviation is shown in Fig. 5 and the RMSE in Fig. 6. The results were obtained using 50,000 Monte Carlo simulations.

The simulated results in Fig. 5 show that the DOA bias is less than0.01 degree, only 0.18& of the true value and the DOA rate bias is less than 1 x 10—4 (0)/s, only 0.7& of the true value. It is compatible to the analysis in Section 4.1 that the DOA and its changing rate are approximatively unbiased. The small bias in the results is caused by the nonlinear transformation in Eqs. (22) and (23) as well as the finite baseline number.

Fig. 5 DOA and its changing rate bias vs. CRPD noise standard deviation.

Fig. 6 DOA and its changing rate RMSE vs. CRPD noise standard deviation.

5 10 15 20 25 30 35 40 45 50 (b) Velocity bias as function CRPD noise standard deviation a ((°)/s)

Fig. 7 Localization bias vs. CRPD noise standard deviation.

Fig. 8 Localization RMSE vs. CRPD noise standard deviation.

As shown in Fig. 6, the DOA and its changing rate RMSE obtained from stage 1 can closely approach to the CRLB when the CRPD noise standard deviation is between 5 (0)/s and 45 (o)/s.

7.2. Performance of localization estimation

In this simulation, we compare the performances of the ML estimator, the PLS estimator and the WIV estimator. The ML estimator was implemented using the GN algorithm. It was initialized to the true value and iterated three times.

The simulation results for the bias of the localization algorithms as a function of the CRPD noise standard deviation is shown in Fig. 7 and the corresponding RMSE is shown in Fig. 8. The bias and RMSE estimates were obtained using 50,000 Monte Carlo simulations.

The ML estimator has the smallest bias. The PLS estimator exhibits significantly larger bias than those of the other two estimators because of the correlation between the measurement matrix and the measurement noise. While the bias of the WIV estimator is not as small as that of the ML estimator, it is nonetheless significantly smaller than that of the PLS estimator.

It can be seen from Fig. 8 that the RMSE of the LS estimator is significantly larger than those of the WIV and ML estimators. The RMSE performance of the WIV estimator closely matches that of the optimal ML estimator and reaches the CRLB when the measured CRPD has an error below 35 (0)/s. It suffers from the thresholding effect at the high noise level due to ignoring the second and high-order error terms in the algorithm development. We also note that the complexity of the WIV estimator is comparable to that of the pseudolinear estimator while the ML estimator has larger computation than the WIV estimator.

8. Conclusions

This paper proposes a two-stage closed-form solution to determine the position and velocity of a moving source using CRPD measurements only. The proposed technique employs least squares only in the two stages and is computationally attractive. It does not have convergence and initialization problems. The proposed two-stage closed-form PLS estimator however is biased because of the correlation between the measurement matrix and the measurement noise. The proposed WIV estimator can reduce the bias as well as the estimation variance. It can attain the CRLB at moderate Gaussian noise.

Acknowledgments

The authors are grateful to Assoc. Prof. Yang Le of Jiangnan University, Wuxi, PR China for discussions and spelling review. They would also like to thank the anonymous reviewers for their critical and constructive review of the manuscript. This study was co-supported by the Foundation of National Defense Key Laboratory of China (No. 9140C860304) and the National High Technology Research and Development Program of China (No. 2011AA7072048).

References

1. Weinstein E. Optimal source localization and tracking from passive array measurements. IEEE Trans Acoust Speech Signal Process 1982;30(1):69-76.

2. Carter GC. Time delay estimation for passive sonar signal processing. IEEE Trans Acoust Speech Signal Process 1981;29(3): 462-70.

3. Yang R, Foo PH, Ng BP, Ng GW. RF emitter geolocation using amplitude comparison with auto-calibrated relative antenna gains. IEEE Trans Aerosp Electron Syst 2011;47(3):2098-110.

4. Rappaport TS, Reed JH, Woerner BD. Position location using wireless communications on highways of the future. IEEE Commun Mag 1996;34(3):33-41.

5. Stansfield RG. Statistical theory of DF fixing. J Inst Electr Eng Part IIIA Radiocommun 1947;94(15):762-70.

6. Torrieri DJ. Statistical theory of passive location systems. IEEE Trans Aerosp Electron Syst 1984;20(2):183-98.

7. Alba PZ, Vidal J, Brooks DH. Closed-form solution for positioning based on angle of arrival of arrival measurements, The 13th IEEE international symposium on personal, indoor and mobile radio communications, Sept 15-18, 2002, Lisbon, Portugal. IEEE Press; 2002.

8. Fowler ML, Chen M, Johnson JA, Zhou Z. Data compression using SVD and Fisher information for radar emitter location. Signal Process 2010;90(7):2190-202.

9. Sathyan T, Sinha A. A two-stage assignment-based algorithm for asynchronous multisensor bearings-only tracking. IEEE Trans Aerosp Electron Syst 2011;47(3):2153-67.

10. Lindgren AG, Gong KF. Position and velocity estimation via bearing observations. IEEE Trans Aerosp Electron Syst 1978;14(4):564-77.

11. Zhang YJ, Xu GZ. Bearings-only target motion analysis via instrumental variable estimation. IEEE Trans Signal Process 2010;58(11):5523-33.

12. Schmidt R. Multiple emitter location and signal parameter estimation. IEEE Trans Antennas Propag 1986;34(3):276-80.

13. Roy R, Kailath T. ESPRIT - estimation of signal parameters via rotational invariance techniques. IEEE Trans Acoust Speech Signal Process 1989;37(7):984-95.

14. Shieh CS, Lin CT. Direction of arrival estimation based on phase differences using neural fuzzy network. IEEE Trans Antennas Propag 2000;48(7):1115-24.

15. Kawase S. Radio interferometer for geosynchronous satellite direction finding. IEEE Trans Aerosp Electron Syst 2007;43(2):443-9.

16. Jacobs E, Ralston EW. Ambiguity resolution in interferometry. IEEE Trans Aerosp Electron Syst 1981;17(6):766-80.

17. See CMS. Sensor array calibration in the presence of mutual coupling and unknown sensor gains and phases. Electron Lett 1994;30(5):373-4.

18. Ng BC, See CMS. Sensor array calibration using a maximum-likelihood approach. IEEE Trans Antennas Propag 1996;44(6): 827-35.

19. Liao B, Zhang ZG, Chan SC. DOA estimation and tracking of ULAs with mutual coupling. IEEE Trans Aerosp Electron Syst 2012;48(1):891-904.

20. Sakamoto Y, Nishio M. Orbit determination using radio interferometer of small-diameter antennas for LEO satellites. IEEE Trans Aerosp Electron Syst 2011;47(3):2111-8.

21. Wan F, Ding JJ, Yu CL. New measurement method for phase difference change rate of radar pulse signals. Syst. Eng. Electron. 2011;33(6):1257-304 [in Chinese].

22. Laneuville D, Vignal H. Grid based target motion analysis, IEEE aerospace conference, March 3-10, 2007, Big Sky, MT. IEEE Press; 2007. p. 1-7.

23. Foy WH. Position location solution by Taylor series estimation. IEEE Trans Aerosp Electron Syst 1976;12(2):187-94.

24. Nardone SC, Lindgren AG, Gong KF. Fundamental properties and performance of conventional bearings-only target motion analysis. IEEE Trans Autom Control 1984;29(9):775-87.

25. Dogancay K. Passive emitter localization using weighted instrumental variables. Signal Process 2004;84(3):487-97.

26. Rao KD, Reddy DC. A new method for finding electromagnetic emitter location. IEEE Trans Aerosp Electron Syst 1994;30(4): 1081-5.

27. Dogancay K. Bearings-only target localization using total least squares. Signal Process 2005;85(9):1695-710.

28. Dogancay K. Bias compensation for the bearings-only pseudolinear target track estimator. IEEE Trans Signal Process 2006;54(1): 59-68.

29. Woods A. Accelerating software radio astronomy FX correlation with GPU and FPGA co-processors [dissertation]. Cape Town (WC): University of Cape Town; 2010.

30. Mendel JM. Lessons in estimation theory for signal processing, communications, and control. Englewood Cliffs (NJ): Prentice-Hall; 1995.

31. Ljung L. System identification: theory for the user. 2nd ed. Englewood Cliffs (NJ): Prentice-Hall; 1999.

32. Kay SM. Fundamentals of statistical signal processing, estimation theory. Englewood Cliffs (NJ): Prentice Hall; 1993.

Zhang Min is a Ph.D. student in the College of Electronic Science and Engineering at National University of Defense Technology (NUDT), Changsha, PR China. He received his B.S. degree in information engineering from Beijing Institute of Technology (BIT) in 2007 and M.S. degree in information and communication engineering from NUDT in 2009, respectively. His main research interests are statistical signal processing, passive localization, target tracking, etc.

Guo Fucheng is a professor and M.S. advisor in the College of Electronic Science and Engineering at NUDT, Changsha, PR China, where he received his Ph.D. degree in 2002. His area of research includes passive localization, target tracking, etc.

Zhou Yiyu is a professor and Ph.D. advisor in the College of Electronic Science and Engineering at NUDT, Changsha, PR China. He received his Ph.D. degree from NUDT in 1992. His area of research includes passive localization, target tracking, etc.