# Joint Position and Velocity Estimation of a Moving Target in Multistatic Radar by Bistatic Range, TDOA, and Doppler Shifts

Authors: Lijuan Yang; Huotao Gao; Boya Li; Yunkun Yang; Guobao Ru
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

In this paper, an estimation algorithm for the position and velocity of a moving target in a multistatic radar system is investigated. Estimation accuracy is improved by using bistatic range (BR), time-difference-of-arrival (TDOA), and Doppler shifts. Multistatic radar system includes several independent receivers and transmitters of time synchronization. Different transmitters radiate signals of different frequencies, and receivers detect the Doppler shifts of the received signals. These estimation parameters, BR, TDOA, and Doppler shifts, are readily available. The proposed algorithm combines different estimated parameters and optimizes estimation accuracy by two-step weighted least squares minimisations (WLS). This estimation algorithm is analysed and verified by simulations, which can reach the Cramer–Rao lower bound (CRLB) performance under mild Gaussian noise when the measurement error is small. Numerical simulations also demonstrate the superior performance of this method.


---

## Body

## 1. Introduction

The classical problem of position and velocity estimation in the multistatic radar system receives great attentions in recent years due to its importance to detection, parameter estimation, tracking, surveillance, and navigation systems [1–3]. Generally, there are two modes for the multistatic radar system [4, 5], one is that a multistatic radar system consists of one transmitter and N receivers (T-RN) [6] and the other is that the multistatic radar system which is composed of M identical transmitters and N identical receivers, radiating in a common coverage area (TM-RN) [3]. The multistatic radar system structure of TM-RN has the advantages in antistealth and antidestruction. In this structure, transmitters radiate signals of different frequencies, and those signals are received by receivers to detect the Doppler shifts. More transmitters and receivers would obtain more estimation parameters. Therefore, the multistatic radar system offers more flexibility and better system robustness as well as improves estimation accuracy of position and velocity.Several localization algorithms about the multistatic system, including radar and sonar, have been developed such as time-of-arrival (TOA) [7], time-difference-of-arrival (TDOA) [8, 9], Doppler shift [10, 11], bistatic range (BR) [2], angle-of-arrival (AOA) [8, 12], and received signal strength (RSS) [13]. In the multistatic system, AOA needs a complex estimation algorithm. BR is equivalent to ellipse localization, which is relied on the propagation time of the transmitted signals travelling from transmitters to the target and reflected back to receivers. The transmitters and receivers are located at different sites [14]. Every single BR can form an ellipse where the target lies on, and the corresponding transmitter and receiver serve as its foci. The intersection of multiple ellipses can determine the approximate position of the target. TDOA, which is the time difference of the signal reflected from the target to different receivers, induces a hyperbola locus for the target to be located, with the associated different receivers as its foci. The intersection of the hyperbola from multiple transmitter-receiver pairs gives the target position. Ellipse and hyperbola can be obtained based on time synchronization of transmitters and receivers. Doppler shifts can be detected through receivers. A moving target creates Doppler shifts which depend on the location and speed of a target. Thus, it can be exploited to improve the localization accuracy in the multistatic radar system [12], and then the accuracy of velocity estimation can be improved.Single parameter can be used to estimate the position or velocity in some cases, but the accuracy needs to be improved. Several localization algorithms have been researched using two types of measurements, such as BR and Doppler shifts [15], BR and bearing measurements [16], and ellipse and hyperbola [17]. Also, there is an estimation algorithm using time delay, Doppler shifts, and AOA measurements [18]. Using two or three kinds of estimated parameters improves the position or velocity accuracy than only one estimated parameter.In this paper, we present a joint position and velocity estimation algorithm, which uses Doppler shifts, BR, and TDOA, to estimate the position and velocity of a moving target. This algorithm is optimized based on two-step weighted least squares minimisations (WLS) [9, 18]. In the first step, a set of pseudolinear BR, TDOA, and Doppler shift equations are established by introducing nuisance parameters, rewriting pseudolinear equations to a linear equation that can utilize the calculation formula of WLS, and then the initial estimation of position and velocity is obtained. In the second step, the error term of the first-stage solution is estimated by using the relationship between initial estimated terms and nuisance variables. Similarly, error terms are calculated by WLS, and then the final estimation of position and velocity is obtained through subtracting the error terms from initial values. Simulation results show that it has a better localization and velocity accuracy not only than the method of ellipse and hyperbola but also than BR and Doppler shifts.We shall use the common convention that bold upper and lower case letters denote matrix and column vector.Im, 0m, and 0m,n represent an m×m identity matrix, a zero vector with the length m, and m×n zero matrix, respectively. The notations AT and A−1 denote transpose and inverse of matrix A. The l2 norm is denoted by ∗. The symbol ⊗ stands for Kronecker products, and a0 denotes the true value of the variable α.


## 2. Problem Formulation for the Multistatic Radar System

We consider the localization of a moving target in a three-dimensional (3D) scenario using the multistatic radar system network with M transmitters and N receivers, whose positions are denoted byti=xit,yit,zitT, i=1,2,…,M and sj=xjs,yjs,zjsT, j=1,2,…,N. The target position u0=x,y,zT, and the target velocity v0=vx,vy,vzT. The geometry structure of the multistatic radar system is shown in Figure 1.Figure 1
Geometry for multistatic radar.Each transmitter radiates a signal of different frequency, and all receivers observe the signal through two paths, one is the direct propagation from the transmitter and the other is the indirect reflection of a moving target. Signals from the two paths for each transmitter-receiver pair provide differential delay and Doppler shift measurements, where the latter comes from the target movement. The BRsri+rj are obtained based on the time synchronization of transmitters and receivers. TDOA is the time difference of signals reflected from the target to spatially separated receivers; it can be obtained due to the time synchronization of receivers. We will take advantage of these readily available parameters, BR, TDOA, and Doppler shift, to estimate the position and velocity of a moving target.The range between the direct and the indirect signal from transmitteri to receiver j is given by(1)rij0=cτ=ri+rj−dij,where c is the speed of signal transmission, τ is the differential delay between the direct and the indirect signal transmission, and rj=u0−sj, ri=u0−ti, and dij=ti−sj; therefore, the BRs of the multistatic radar system are rij0+dij=cτ+dij.Collecting the range measurements gives(2)r=r1T,r2T,…,rMTT=r0+Δr,where ri=ri1,ri2,…,riNT=ri0+Δri, Δri=Δri1,Δri2,…,ΔriNT.It is assumed that receivers1 is the reference receiver. After multiplying TDOA by signal propagation speed, the range difference between the receivers and s1 can be expressed as rj10=rj−r1, j=2,3,…,N and the measured version of rj10 is rj1=rj10+Δrj1, where Δrj1 is the additive noise. We can get the N − 1 hyperbolic measurements in vector form as(3)r˜=r21,r31,…,rN1T=r˜0+Δr˜.Doppler shift is defined as the rate of the change of total path length of the signal to the time because transmitters and receivers are stationary, and the true Doppler shift is(4)fij0=1λdridt+drjdt=fcicρu0tiTv0+ρu0sjTv0,where ρu0ti=u0−ti/u0−ti denotes a unit vector from ti to u0 and fci is the carrier frequency of transmitter ti; the measured fij0 is fij=fij0+Δfij, where Δfij is the additive noise. Collecting the Doppler measurements gets(5)f=f1T,f2T,…,fMTT=f0+Δf,where fi=fi1,fi2,…,fiNT=fi0+Δfi, Δfi=Δfi1,Δfi2,…,ΔfiNT.Definingm=rT,r˜T,fTT, m is a 2MN+N−1×1 measurement vector, the error vector Δm=ΔrT,Δr˜T,ΔfTT is assumed to be zero-mean Gaussian with covariance matrix Qm=diagQr,Qr˜,Qf, where Qr, Qr˜, and Qf are covariance matrices of r, r˜, and f. We will use BR, TDOA, and Doppler shifts to estimate the position and velocity of a moving target.The proposed algorithm is optimized by two steps. First, get the initial estimates through establishing a series of pseudolinear equations which introduce nuisance parameters. Second, get the final estimates by exploring the estimation error and estimating the amount of correction to refine the solution.


### 2.1. First Step

We calculate the initial position and velocity by introducing nuisance parameters and establishing a set of pseudolinear equations about BR, TDOA, and Doppler shift.Rearrangingrij0=ri+rj−dij, substituting rij=rij0+Δrij, and squaring both sides yield the elliptic measurement equations. It can be obtained that(6)2rjΔrij=rij2+2rijdij+2tiTti−sj−2ti−sjTu0−2rij+dijri.Note that BR and ellipse localization will be used interchangeably in this paper.Similarly, rearrangingrj10=rj−r1, substituting rj1=rj10+Δrj1, and squaring both sides yield the hyperbolic measurement equations. It can be expressed as(7)2rj1Δrj1=rj12−sjTsj+s1Ts1+2sj−s1Tu0+2rj1u0−s1.Note that TDOA and hyperbola will be used interchangeably in this paper.Then, rearranging (4) and fij=fij0+Δfij, it can be expressed as(8)2cΔfijrj+2fciρu0sjTv0Δrij=2cfijrij+2cfijdij−2fciti−sjTv0−2cfijri−2fijrij+dijρu0tiTv0.In (6)–(8), the second-order noise terms have been ignored, and the unknown vector in the first stage can be φ0=u0T,ri0,u0−s1,v0T,βi0, where βi0=ρu0tiTv0 which contains 2M+7 nuisance variables. Equations (6)–(8) can be rewritten in the matrix form as follows:(9)B1Δm=h1−G1φ0.This is a typical equation that can be solved using WLS, whereh1=h1rT,h1r˜T,h1fTT, h1r=rij2+2rijdij+2tiTti−sj, h1r˜=rj12−sjTsj+s1Ts1, and h1f=2cfijrij+2cfijdij;(10)G1=G1rT,G1r˜T,G1fTT,G1r=2ti−sjT,0i−1T,rij+dij,02M−i+4T,G1r˜=2s1−sjT,0MT,−s1−sj,0M+3T,G1f=20i+2T,cfij,0M−i+1T,fciti−sjT,0i−1T,fcirij+dij,0M−iT,B1=blkdiagB1r,B1r˜,B1f,B1r=2IM⊗diagrj,B1r˜=2diagrj1,B1f=diag2fciρu0sjTv0,0MN,N−1,2cIMdiagrj.The first-step WLS solution is [19, 20](11)φ1=G1TW1G1−1G1TW1h1,where W1=B1−TQ−1B1−1 is the weighting matrix. Since the variable B1 depends on the estimated position and velocity, however, the initial estimations are unknown in the beginning; hence, we can first set W1=Qm−1 to produce an initial solution, then use the initial solution to recalculate B1, and finally expected W1 can be obtained. When the noise is small, the covariance matrix of φ1 is approximately equal to covφ1≈G1TW1G1−1.


### 2.2. Second Step

The error terms are estimated by using the relationship between the solution of the first stage and nuisance variables.We denote the first-stage estimation asφ1=u˜T,r˜i,u˜−s1,v˜T,β˜i=φ10+Δφ1, where Δφ1=Δu˜T,Δr˜i,Δu˜−s1,Δv˜T,Δβ˜i, ri=r˜i−ΔriRearranging the parameters ofφ1 and Δφ1, we get the equations about the nuisance variables and the estimation values from the first step:(12)2r˜iΔr˜i=r˜i2−u˜Tu˜−tiTti+2tiTu˜+2u˜−tiTΔu,(13)2u˜−s1Δu˜−s1=u˜−s12−s1Ts1−u˜Tu˜+2s1Tu˜+2u˜−s1TΔu,(14)β˜iΔr˜i+r˜iΔβ˜i=r˜iβ˜i−u˜Tv˜+tiTv˜+v˜TΔu+u˜−tiTΔv.Stacking (12)–(14) and putting them together yield the matrix form equation (15) which can be used to obtain estimation error terms by WLS:(15)B2Δφ1=h2−G2Δϑ,where(16)Δϑ=ΔuT,ΔvTT,B2=B21T,B22TT,B21=blkdiagI3,2ri,2u˜−s1,0M+4,M+3,B22=0M,3,diagβi,0M,4,diagri,h2=01,3T,h2rT,h2r˜T,01,3T,h2fTT,h2r=r˜i2−u˜Tu˜−tiTti+2tiTu˜,h2r˜=u˜−s12−s1Ts1−u˜Tu˜+2s1Tu˜,h2f=r˜iβ˜i−u˜Tv˜+tiTv˜,G2=−I3,03,3T,G2rT,G2r˜T,03,3,−I3T,G2fT,G2r=−2u˜−tiT,0M,3,G2r˜=−2u˜−s1T,01,3,G2f=−v˜T,−u˜−tiT.The second-step WLS solution is(17)φ2=Δϑ=G2TW2G2−1G2TW2h2,where W2=B2−Tcovφ1B2−1. Δu=φ21:3 and Δv=φ24:6. Under small noise conditions, we can get(18)covφ2≈G2TW2G2−1,and u˜=φ11:3, v˜=φ14+M:6+M can be obtained from (11).Finally, the estimations of position and velocity are(19)u=u˜−Δu,v=v˜−Δv.


## 2.1. First Step

We calculate the initial position and velocity by introducing nuisance parameters and establishing a set of pseudolinear equations about BR, TDOA, and Doppler shift.Rearrangingrij0=ri+rj−dij, substituting rij=rij0+Δrij, and squaring both sides yield the elliptic measurement equations. It can be obtained that(6)2rjΔrij=rij2+2rijdij+2tiTti−sj−2ti−sjTu0−2rij+dijri.Note that BR and ellipse localization will be used interchangeably in this paper.Similarly, rearrangingrj10=rj−r1, substituting rj1=rj10+Δrj1, and squaring both sides yield the hyperbolic measurement equations. It can be expressed as(7)2rj1Δrj1=rj12−sjTsj+s1Ts1+2sj−s1Tu0+2rj1u0−s1.Note that TDOA and hyperbola will be used interchangeably in this paper.Then, rearranging (4) and fij=fij0+Δfij, it can be expressed as(8)2cΔfijrj+2fciρu0sjTv0Δrij=2cfijrij+2cfijdij−2fciti−sjTv0−2cfijri−2fijrij+dijρu0tiTv0.In (6)–(8), the second-order noise terms have been ignored, and the unknown vector in the first stage can be φ0=u0T,ri0,u0−s1,v0T,βi0, where βi0=ρu0tiTv0 which contains 2M+7 nuisance variables. Equations (6)–(8) can be rewritten in the matrix form as follows:(9)B1Δm=h1−G1φ0.This is a typical equation that can be solved using WLS, whereh1=h1rT,h1r˜T,h1fTT, h1r=rij2+2rijdij+2tiTti−sj, h1r˜=rj12−sjTsj+s1Ts1, and h1f=2cfijrij+2cfijdij;(10)G1=G1rT,G1r˜T,G1fTT,G1r=2ti−sjT,0i−1T,rij+dij,02M−i+4T,G1r˜=2s1−sjT,0MT,−s1−sj,0M+3T,G1f=20i+2T,cfij,0M−i+1T,fciti−sjT,0i−1T,fcirij+dij,0M−iT,B1=blkdiagB1r,B1r˜,B1f,B1r=2IM⊗diagrj,B1r˜=2diagrj1,B1f=diag2fciρu0sjTv0,0MN,N−1,2cIMdiagrj.The first-step WLS solution is [19, 20](11)φ1=G1TW1G1−1G1TW1h1,where W1=B1−TQ−1B1−1 is the weighting matrix. Since the variable B1 depends on the estimated position and velocity, however, the initial estimations are unknown in the beginning; hence, we can first set W1=Qm−1 to produce an initial solution, then use the initial solution to recalculate B1, and finally expected W1 can be obtained. When the noise is small, the covariance matrix of φ1 is approximately equal to covφ1≈G1TW1G1−1.


## 2.2. Second Step

The error terms are estimated by using the relationship between the solution of the first stage and nuisance variables.We denote the first-stage estimation asφ1=u˜T,r˜i,u˜−s1,v˜T,β˜i=φ10+Δφ1, where Δφ1=Δu˜T,Δr˜i,Δu˜−s1,Δv˜T,Δβ˜i, ri=r˜i−ΔriRearranging the parameters ofφ1 and Δφ1, we get the equations about the nuisance variables and the estimation values from the first step:(12)2r˜iΔr˜i=r˜i2−u˜Tu˜−tiTti+2tiTu˜+2u˜−tiTΔu,(13)2u˜−s1Δu˜−s1=u˜−s12−s1Ts1−u˜Tu˜+2s1Tu˜+2u˜−s1TΔu,(14)β˜iΔr˜i+r˜iΔβ˜i=r˜iβ˜i−u˜Tv˜+tiTv˜+v˜TΔu+u˜−tiTΔv.Stacking (12)–(14) and putting them together yield the matrix form equation (15) which can be used to obtain estimation error terms by WLS:(15)B2Δφ1=h2−G2Δϑ,where(16)Δϑ=ΔuT,ΔvTT,B2=B21T,B22TT,B21=blkdiagI3,2ri,2u˜−s1,0M+4,M+3,B22=0M,3,diagβi,0M,4,diagri,h2=01,3T,h2rT,h2r˜T,01,3T,h2fTT,h2r=r˜i2−u˜Tu˜−tiTti+2tiTu˜,h2r˜=u˜−s12−s1Ts1−u˜Tu˜+2s1Tu˜,h2f=r˜iβ˜i−u˜Tv˜+tiTv˜,G2=−I3,03,3T,G2rT,G2r˜T,03,3,−I3T,G2fT,G2r=−2u˜−tiT,0M,3,G2r˜=−2u˜−s1T,01,3,G2f=−v˜T,−u˜−tiT.The second-step WLS solution is(17)φ2=Δϑ=G2TW2G2−1G2TW2h2,where W2=B2−Tcovφ1B2−1. Δu=φ21:3 and Δv=φ24:6. Under small noise conditions, we can get(18)covφ2≈G2TW2G2−1,and u˜=φ11:3, v˜=φ14+M:6+M can be obtained from (11).Finally, the estimations of position and velocity are(19)u=u˜−Δu,v=v˜−Δv.


## 3. Performance Analysis

The CRLB for estimating the unknownϑ=uT,vTT is equal to the inverse of the Fisher information matrix (FIM):(20)CRLBϑ=∇TQm−1∇−1,where ∇=∂ϑ/∂u,∂ϑ/∂v.The covariance matrix ofϑ is obtained from (21) as [1]:(21)covϑ=covφ2=G3TQm−1G3−1,where G3=B1−1G1B2−1G2.It can be obtained thatG3≈∇ when Δrij≪rij, Δrj1≪rj1, Δfij≪fij, Δrij≪u0−sj, Δrj1≪u0−sj, and cΔfij≪u0−sj. The calculation of this follows similar procedures as in [2, 3, 15, 20–23].Consequently, we conclude from (20) and (21) that(22)covϑ≈CRLBϑ,when those small noise conditions are satisfied.


## 4. Simulation Results

The simulation scenario is as follows. A moving target is atu0=1000,1500,3000Tm and its velocity is v0=120,150,100Tm/s. The multistatic radar system has T = 7 transmitters and R = 8 receivers. The transmitters are at t1=3000,1000,0Tm, t2=1000,1500,1000Tm, t3=3000,5000,2000Tm, t4=3500,2000,1000Tm, t5=3500,2000,0Tm, t6=0,0,0Tm, and t7=−1000,3000,0Tm. The receivers’ position are at s1=0,5000,0Tm, s2=5000,0,1000Tm, s3=−5000,0,1500Tm, s4=0,−5000,1000Tm, s5=0,0,1000Tm, s6=0,1000,1000Tm, s7=0,−2000,0Tm, and s8=−5000,2000,0Tm. The signal propagation speed c is 300000000 m/s. The carrier frequencies of the transmitted signals are fci=20,19,12,26,28,16,8MHz. The covariance matrix Qm of BR, TDOA, and Doppler shift measurements is Qm=σ2blkdiagIM,IN−1,105IM, σ represents the level of noise. The localization and velocity accuracy are assessed via root-mean-square error (RMSE). The results in Figures 2–8 are obtained from Monte Carlo simulations of 5,000 ensemble runs. We use CRLB as a benchmark for performance evaluation.Figure 2
Comparison of the position RMSE of the proposed method with AE + Doppler, AE + H, and CRLB whenT = 1, R = 8. The accuracy is lower as the noise power increases.Figure 3
Comparison of the position RMSE of the proposed method with AE + Doppler, AE + H, and CRLB whenT = 3, R = 8. The CRLB is close to zero, but it is going up as the noise power increases.Figure 4
Comparison of the position RMSE of the proposed method with AE + Doppler, AE + H, and CRLB whenT = 7, R = 8. The CRLB is close to zero, but it is going up as the noise power increases.Figure 5
Comparison of the position RMSE of the proposed method with AE + Doppler, AE + H, and CRLB whenT = 3, R = 7. The CRLB is close to zero, but it is going up as the noise power increases.Figure 6
Comparison of the velocity RMSE of the proposed method with AE + Doppler and CRLB whenT = 1, R = 8. The accuracy is lower as the noise power increases.Figure 7
Comparison of the velocity RMSE of the proposed method with AE + Doppler and CRLB whenT = 3, R = 8. The CRLB is close to zero, but it is going up as the noise power increases.Figure 8
Comparison of the velocity RMSE of the proposed method with AE + Doppler and CRLB whenT = 3, R = 7. The CRLB is close to zero, but it is going up as the noise power increases.We compared performance of the proposed method with the CRLB and two kinds of algorithms, one is the hybrid of ellipse and hyperbola algorithm and the other is the algorithm of joint BR and Doppler shifts. The comparison algorithms were shown to outperform the same kind of algorithms. Simulation results are summarised in Figures2–8, in which the curve marked with circle corresponds to the hybrid of ellipse, hyperbola, and Doppler shifts (AE + H + Doppler), and curves marked with triangle and plus correspond to the joint ellipse and Doppler shifts (AE + Doppler) and the joint ellipse and hyperbola (AE + H), respectively. The dotted line represents the CRLB under the corresponding setting.In this paper, two kinds of modes are simulated, one is one transmitter and multiple receivers and the other is multiple transmitters and multiple receivers. The results show that the proposed method is able to reach the CRLB when the measurement error is small. As is shown in Figures2–5, under the condition of changing the number of transmitters or receivers, the proposed method has a better performance in the positioning compared to the hybrid of elliptic and hyperbolic method and the joint method of BR and Doppler shifts. Also, with the increase of σ, the position RMSE of the proposed method is obviously lower than that of the hybrid of elliptic and hyperbolic method and the joint method of BR and Doppler shifts. Then, in Figures 6–8, the velocity RMSE is also analysed, and with the increase of σ, its gap between the proposed method and the hybrid of BR and Doppler shifts method becomes larger. The advantage of estimation performance is more obvious with the increase of measurement error. Moreover, when the number of receivers or transmitters is changed, the proposed algorithm has much better velocity accuracy than the joint method of BR and Doppler shifts. In Figures 3–5, 7, and 8, CRLB is close to zero when the number of transmitters and receivers is large. It is illustrated that the proposed algorithm has much better performance not only in the positioning but also in the velocity. Moreover, the estimation performance is affected by the number of transmitters and receivers, and estimation accuracy is improved when the number of transmitters and receivers is increased.


## 5. Conclusion

In this paper, we investigated the location and velocity of a moving target in the multistatic radar system by combining BR, TDOA, and Doppler shifts. This nonlinear estimation problem is solved by introducing nuisance variables and using two-step WLS to refine the estimate. It is found that the proposed method can obtain the CRLB accuracy over a small error region in the case of Gaussian noise. Simulation results have also shown its effectiveness. The proposed algorithm in this paper enriches the positioning methods of moving targets in the multistatic radar system.



---

*Source: 4943872-2019-11-23.xml*