Geolocation and Tracking by TDOA Measurements Based on Space–Air–Ground Integrated Network

Abstract

Due to the development of manufacturing and launch technologies for satellites, there are now more and more satellite networks. Hence, cooperative reconnaissance is possible to implement among satellite networks, aerial vehicles and ground stations. In this paper, we study the method of geolocation and tracking by time difference of arrival (TDOA) measurements based on space–air–ground integrated (SAGI) network. We first analyze the Cramer Rao lower bound (CRLB) for the source localization accuracy in different coordinate systems. Then, we compare the effects of different system errors, such as clock synchronization error, position bias of the observers, elevation bias of the target and non-horizontal velocity of the target. Further, we also develop a maximum likelihood (ML) estimator for target position and velocity. Finally, the theoretical performance of the proposed estimator is validated via computer simulations.

1. Introduction

Passive source localization has been an important study aspect in the area of electronic reconnaissance. The localization method can be classified as single station [1,2,3] or multi-station [4,5] by the number of observers. Localization approaches can also be divided by the type of observers: ground-based [6], air-based [7] and space-based localization [8]. Commonly used signal measurements include time of arrival (TOA) [9], time difference of arrival [10,11], frequency difference of arrival (FDOA) [12] and angle of arrival (AOA) [13,14]. The attitude of observer should be measured accurately for the localization approach with AOA. Clock synchronization between observer stations should be precise to ensure position accuracy when using TDOA or FDOA. In practice, clock synchronization is commonly precise to the level of a nanosecond. Hence, many popular passive localization approaches utilize TDOA measurements. Over the years, many approaches have been proposed, such as the iterative Taylor-series method [15], two-step least squares [10], Hough transform [16], particle swarm optimization [17], semi-definite relaxation [18] and so on.

However, most current research focuses on the single-station type. Localization accuracy and detection range could not accommodate usability requirements while using aircraft or ground stations only. There are few studies on source localization using different observation platforms. Due to the development of manufacturing and launch technology for satellites, there are now more and more satellite networks. Hence, cooperative reconnaissance is possible to implement among satellite networks, aerial vehicles and ground stations. Therefore, it is highly expected to enhance passive localization performance by using SAGI networks.

With the development of communication technology, SAGI communication is envisioned to achieve significant performance gain in coverage, reliability and flexibility, which has attracted great interests over the past few years [19]. Many companies have also made extensive efforts in the area of SAGI communication [20,21], such as SpaceX and Google. Therefore, it is viable to realize source localization based on SAGI networks. Li et al. [22] developed an active cooperative localization method for SAGI networks based on TOA measurements. The authors also analyzed the effects of clock offset among ground agents. However, there is no literature about the problem of passive localization based on SAGI networks.

In this paper, we study the method of geolocation and tracking by TDOA measurements based on SAGI networks. The conference article [23] presented some preliminary work about performance analysis of geolocation based on SAGI networks. The contributions of this paper mainly include the following.

• (1). We derive the CRLB of the source localization in different coordinate systems. Most research about the CRLB of geolocation is derived in the standard Earth-centered Earth-fixed (ECEF) coordinate system [24,25]. Under the constraint conditions of source position and velocity on the Earth, the Fisher information matrix (FIM) of the source location in the ECEF coordinate system may be occasionally nonreversible. For the purpose of more accurate estimation, we also derive the CRLB for the target modeled in the geodetic coordinate system.

• (2). Our second contribution is to analyze the effects of different system errors, such as clock synchronization error, position bias of the observers, elevation bias of the target and non-horizontal velocity of the target. Meanwhile, most existing research focuses on the effect of random errors [26,27,28]. However, it should be remarked that system biases are common in the scenario of passive geolocation based on SAGI networks, where accurate synchronization among observers is an impractical assumption.

• (3). Note that most previous research concentrates on estimation of source position with TDOA measurements. Due to the spherical constraint conditions of source position and velocity, previous algorithms are not suitable for jointly estimating the source position and velocity [29,30] with TDOA measurements only. Therefore, we propose an iterative maximum likelihood estimator for both position and velocity in the scenario of geolocation based on SAGI networks.

The remainder of this paper is organized as follows. In Section 2, we develop the CRLB for the source localization in different coordinate systems and analyze the effects of different system errors which contain clock synchronization error, position bias of the observers, elevation bias of the target and non-horizontal velocity of the target. Section 3 presents the maximum likelihood estimator of the source position and velocity. Numerical investigations for three different geolocation scenarios are reported in Section 4. Conclusions are presented in Section 5.

2. Theoretical Performance

2.1. Problem Statement and Mathematical Formulation

Suppose the target is on the surface of the Earth with known constant height $h$. The original position and velocity of the target are $u_o$ and ${\dot{\mathbf{u}}}_o$, which are defined in the ECEF coordinate system. The velocity is assumed to be horizontal and fixed during the whole observation time. The constrained conditions of $u_o$ and ${\dot{u}}_o$ are shown as follows.

(1)$\left\{\begin{array}{l}{u}_o^T{u}_o={\left(N+h\right)}^2\\ u_o^T{\dot{u}}_o=0\end{array}\right.$

where $N$ is the radius of curvature in prime vertical, $N=\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}$. $a$ is the semi-major axis of the earth and $e$ is the first eccentricity. $B$ is the latitude of the target; ${\dot{u}}_o$ and $h$ will be zeros when the emitter is stationary.

Assume there are $M$ moving observers in the SAGI network. Every observer intercepts $K$ chips of signal. The true position of the $j\mathrm{th}$ observer at the time of intercepting the $k\mathrm{th}$ chip of signal is $s_{j,k}^o$. The noisy value $s_{j,k}$ is available and the random error is expressed as $\Delta s_{j,k}$.

Then, the TDOA measurement between the $j\mathrm{th}$ observer and the $0\mathrm{th}$ observer for the $k\mathrm{th}$ chip of signal is

(2)$t_{j,0,k}=\frac{‖u_k-s_{j,k}^o‖}{c}-\frac{‖u_k-s_{0,k}^o‖}{c}+\Delta t_{j,0,k},j=0,\dots ,M-1;k=1,\dots ,K.$

where $\Delta t_{j,0,k}$ is the TDOA measured noise, and $j=0,\dots ,M-1$, $k=1,\dots ,K$. Due to the assumption of uniform motion, the position and velocity of the emitter at other times are expressed as $u_k=u_0+\left(k-1\right)T{\dot{u}}_0$ and ${\dot{u}}_k={\dot{u}}_0$. $T$ is the transmitted interval of signal chips.

Collecting the TDOA measurements and observer locations yields

(3)$\begin{array}{l}z=z^o+n\\ z={\left[z_t^T,z_{\mathrm{nav}}^T\right]}^T\\ z_t={\left[z_{t,1}^T,z_{t,2}^T,\dots ,z_{t,K}^T\right]}^T\\ z_{\mathrm{t},\mathrm{k}}={\left[t_{1,0,k},t_{2,0,k},\dots ,t_{M-1,0,k}\right]}^T\\ z_{\mathrm{nav}}={\left[z_{nav,1}^T,z_{nav,2}^T,\dots ,z_{nav,K}^T\right]}^T\\ z_{\mathrm{nav},\mathrm{k}}={\left[s_{0,k}^T,s_{1,k}^T,\dots ,s_{M-1,k}^T\right]}^T\end{array}$

where $z^o$ represents the true values of TDOA measurements and observer locations and $n$ is the measurement noise vector of $z^o$.

(4)$\begin{array}{l}{z}^o={\left[z_t^{oT},z_{nav}^{oT}\right]}^T\\ z_t^o={\left[z_{t,1}^{oT},z_{t,2}^{oT},\dots ,z_{t,K}^{oT}\right]}^T\\ z_{t,k}^o={\left[t_{1,0,k}^o,t_{2,0,k}^o,\dots ,t_{M-1,0,k}^o\right]}^T\\ z_{nav}^o={\left[z_{nav,1}^{oT},z_{nav,2}^{oT},\dots ,z_{nav,K}^{oT}\right]}^T\\ z_{nav,k}^o={\left[s_{0,k}^{oT},s_{1,k}^{oT},\dots ,s_{M-1,k}^{oT}\right]}^T\end{array}$

(5)$\begin{array}{l}n={\left[n_t^T,n_{nav}^T\right]}^T\\ n_t={\left[n_{t,1}^T,n_{t,2}^T,\dots ,n_{t,K}^T\right]}^T\\ n_{t,k}={\left[\Delta t_{1,0,k},\Delta t_{2,0,k},\dots ,\Delta t_{M-1,0,k}\right]}^T\\ n_{nav}={\left[n_{nav,1}^T,n_{nav,2}^T,\dots ,n_{nav,K}^T\right]}^T\\ n_{nav,k}={\left[\Delta s_{0,k}^T,\Delta s_{1,k}^T,\dots ,\Delta s_{M-1,k}^T\right]}^T\end{array}$

Suppose the observer location error $n_{nav}$ and TDOA measurement error $n_t$ are zero-mean Gaussian random vectors with covariance matrices $C_{nav}=diag\left\{C_{nav,1},\dots ,C_{nav,K}\right\}$ and $C_t=diag\left\{C_{t,1},\dots ,C_{t,K}\right\}$. Hence, the covariance matrix of $n$ is

(6)$C=\left[\begin{array}{cc}{C}_t& 0\\ 0& C_{nav}\end{array}\right]$

The task of this study is to identify the target’s position and velocity based on the SAGI network, where the measurement vector is (3) and the constrained conditions are (1).

2.2. CRLB

CRLB presents the best accuracy than an unbiased estimator could achieve. Due to the constrained conditions of the target’s position and velocity on the earth, the derivations of CRLB can be classified into two types. In this section, we evaluate the two types of CRLB for passive geolocation and tracking based on SAGI network as follows. The first type of method is developed for the target denoted in the ECEF coordinate system. We first derive the FIM without using the constrained conditions. Then, the constrained CRLB is generated by the inversion of FIM and gradient matrix of the constraints [31]. However, the FIM may be occasionally nonreversible. For the purpose of more accurate estimation, the second type of the CRLB is derived for the target modeled in the geodetic coordinate system.

2.2.1. CRLB in ECEF Coordinate System

Define the unknown parameter as ${\mathsf{\theta }}_1={\left[{\mathsf{\theta }}_u^T,{\mathsf{\theta }}_{nav}^T\right]}^T$, where the target’s location vector is ${\mathsf{\theta }}_u={\left[u_o^T,{\dot{u}}_o^T\right]}^T$, and ${\mathsf{\theta }}_{nav}=z_{nav}^o$ is the location vector of the observers. The FIM of the unknown parameter without using the constrained conditions in (1) can be derived as

(7)$I_1\left({\mathsf{\theta }}_1\right)=J_1^T{C}^{-1}{J}_1$

where

(8)$J_1=\frac{\partial z}{\partial {\mathsf{\theta }}_1}=\left[\begin{array}{cc}{J}_u& J_{t,nav}\\ 0& E_{nav}\end{array}\right]$

(9)$J_u=\frac{\partial z_t}{\partial {\mathsf{\theta }}_u},J_{t,nav}=\frac{\partial z_t}{\partial {\mathsf{\theta }}_{nav}},E_{nav}=\frac{\partial z_{nav}}{\partial {\mathsf{\theta }}_{nav}}$

Then, the CRLB of ${\mathsf{\theta }}_1$ without using the constraint conditions is

(10)$\mathrm{CRLB}\left({\mathsf{\theta }}_1\right)={\left\{{\left[\begin{array}{cc}{J}_u& J_{t,nav}\\ 0& E_{nav}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{cc}{C}_t^{-1}& 0\\ 0& C_{nav}^{-1}\end{array}\right]\left[\begin{array}{cc}{J}_u& J_{t,nav}\\ 0& E_{nav}\end{array}\right]\right\}}^{-1}={\left[\begin{array}{cc}{A}_1& U_1\\ V_1& D\end{array}\right]}^{-1}$

Hence, the CRLB of ${\mathsf{\theta }}_u$ without constrained conditions is

(11)$\mathrm{CRLB}\left({\mathsf{\theta }}_u\right)={\mathit{A}}_1^{-1}+{\mathit{A}}_1^{-1}{\mathit{U}}_1{\left(D-V_1{\mathit{A}}_1^{-1}{U}_1\right)}^{-1}{V}_1{\mathit{A}}_1^{-1}$

${\mathit{A}}_1^{-1}$ is the CRLB of ${\mathsf{\theta }}_u$ if the observer locations are known accurately. ${\mathit{A}}_1^{-1}{U}_1{\left(D-V_1{A}_1^{-1}{U}_1\right)}^{-1}{V}_1{\mathit{A}}_1^{-1}$ represents the performance degradation due to the random errors of observer locations. $A_1$, $U_1$, $V_1$ and $D$ are defined as follows.

(12)$A_1=J_u^{\mathrm{T}}{C}_t^{-1}{J}_u={\displaystyle \sum _{k=1}^K{J}_{u,k}^{\mathrm{T}}{C}_{t,k}^{-1}{J}_{u,k}}$

(13)$U_1=J_u^{\mathrm{T}}{C}_{\mathrm{t}}^{-1}{J}_{t,nav}=diag\left\{J_{u,1}^{\mathrm{T}}{C}_{t,1}^{-1}{J}_{t,nav,1},\dots ,J_{u,K}^{\mathrm{T}}{C}_{t,K}^{-1}{J}_{t,nav,K}\right\}$

(14)$V_1=J_{t,nav}^{\mathrm{T}}{C}_t^{-1}{J}_u=diag\left\{J_{t,nav,1}^{\mathrm{T}}{C}_{t,1}^{-1}{J}_{u,1},\dots ,J_{t,nav,K}^{\mathrm{T}}{C}_{t,K}^{-1}{J}_{u,K}\right\}$

(15)$\begin{array}{l}D=J_{t,nav}^{\mathrm{T}}{C}_t^{-1}{J}_{t,nav}+C_{nav}^{-1}\\ =diag\left\{J_{t,nav,1}^{\mathrm{T}}{C}_{t,1}^{-1}{J}_{t,nav,1}+C_{nav,1}^{-1},\dots ,J_{t,nav,K}^{\mathrm{T}}{C}_{t,K}^{-1}{J}_{t,nav,K}+C_{nav,K}^{-1}\right\}\end{array}$

where the definitions of $J_{\mathrm{u},\mathrm{k}}$ and $J_{t,nav,k}$ can be expressed as

(16)$J_{u,k}=\left[\frac{\partial t_{1,0,k}}{\partial {\mathsf{\theta }}_u};\frac{\partial t_{2,0,k}}{\partial {\mathsf{\theta }}_u};\dots ;\frac{\partial t_{M-1,0,k}}{\partial {\mathsf{\theta }}_u}\right]$

(17)$\frac{\partial t_{j,0,k}}{\partial {\mathsf{\theta }}_u}=\left[\frac{\partial t_{j,0,k}}{\partial u_k}\times \frac{\partial u_k}{\partial u_o},\frac{\partial t_{j,0,k}}{\partial u_k}\times \frac{\partial u_k}{\partial {\dot{u}}_o}\right]$

(18)$\frac{\partial t_{j,0,k}}{\partial u_k}=\frac{u_k^T-s_{j,k}^T}{c‖u_k-s_{j,k}‖}-\frac{u_k^T-s_{0,k}^T}{c‖u_k-s_{0,k}‖}$

(19)$\frac{\partial u_k}{\partial u_o}=I_{3\times 3},\frac{\partial u_k}{\partial {\dot{u}}_o}=\left(k-1\right)T\times I_{3\times 3}$

(20)$J_{t,nav,k}=\frac{\partial z_{t,k}}{\partial z_{nav,k}}=\left[\frac{\partial t_{1,0,k}}{\partial z_{nav,k}};\frac{\partial t_{2,0,k}}{\partial z_{nav,k}};\dots ;\frac{\partial t_{M-1,0,k}}{\partial z_{nav,k}}\right]$

(21)$\frac{\partial t_{j,0,k}}{\partial z_{nav,k}}=\left[\frac{u_k^T-s_{0,k}^T}{c‖u_k-s_{0,k}‖},0_{1\times 3(\mathrm{j}-1)},-\frac{u_k^T-s_{j,k}^T}{c‖u_k-s_{j,k}‖},0_{1\times 3(\mathrm{M}-1-\mathrm{j})}\right]$

The constrained conditions in (1) can be rewritten as

(22)$\mathit{G}\left({\mathsf{\theta }}_1\right)=\left[\begin{array}{c}{u}_0^T{u}_0-{\left(N+h\right)}^2\\ u_0^T{\dot{u}}_0\end{array}\right]$

The gradient matrix of (22) with respect to the unknown parameter ${\mathsf{\theta }}_1$ is

(23)$F={\left[\frac{\partial G\left({\mathsf{\theta }}_1\right)}{\partial {\mathsf{\theta }}_1^T}\right]}^{\mathrm{T}}=\left[\begin{array}{c}{F}_u\\ O\end{array}\right]$

where

(24)$F_u={\left[\frac{\partial G\left({\mathsf{\theta }}_1\right)}{\partial {\mathsf{\theta }}_u^T}\right]}^{\mathrm{T}}={\left[\begin{array}{cc}2u_0^T& 0\\ {\dot{u}}_0^T& u_0^T\end{array}\right]}^{\mathrm{T}}$

The constrained CRLB of ${\mathsf{\theta }}_1$ is derived as [31]

(25)$\mathrm{CRLB}\left({\mathsf{\theta }}_1\right)=I_1^{-1}-I_1^{-1}F{\left(F^{\mathrm{T}}{I}_1^{-1}F\right)}^{-1}{F}^{\mathrm{T}}{I}_1^{-1}$

2.2.2. CRLB in Geodetic Coordinate System

The position and velocity of the target in the ECEF coordinate system can be transformed to the geodetic coordinate system. The dimension of the initial location vector of the emitter can be reduced from 6 to 4 by the use of the constrained conditions in (1). The initial position and velocity of the emitter in the geodetic coordinate system are ${\left[L_o,B_o,h\right]}^T$ and ${\left[V_L,V_B,0\right]}^T$, where $h$ is the known constant. Hence, the unknown position and velocity vectors are ${\alpha }_o={\left[L_o,B_o\right]}^T$ and ${\dot{\alpha }}_o={\left[V_L,V_B\right]}^T$, respectively. Due to the assumption of uniform motion, the position and velocity of the target at other times are expressed as ${\alpha }_k={\alpha }_0+\left(k-1\right)T{\dot{\alpha }}_0$ and ${\dot{\alpha }}_k={\dot{\alpha }}_0$. $T$ is the transmitted interval of signal chips. The transformation relation from geodetic to ECEF can be expressed as

(26)$\left\{\begin{array}{l}{x}_o=\left(N+h\right)\cos B_o\cos L_o\\ y_o=\left(N+h\right)\cos B_o\sin L_o\\ z_o=\left(N\left(1-e^2\right)+h\right)\sin B_o\end{array}\right.$

Define the unknown parameter ${\mathsf{\theta }}_2={\left[{\mathsf{\theta }}_{\alpha }^T,{\mathsf{\theta }}_{nav}^T\right]}^T$, where the vector of emitter location is ${\mathsf{\theta }}_{\alpha }={\left[{\alpha }_o^T,{\dot{\alpha }}_0^T\right]}^T$, and the vector of observer location is ${\mathsf{\theta }}_{nav}=z_{nav}^o$. Hence, the FIM can be defined as

(27)$I_2\left({\mathsf{\theta }}_2\right)=J_2^T{C}^{-1}{J}_2$

where

(28)$J_2=\frac{\partial z}{\partial {\mathsf{\theta }}_2}=\left[\begin{array}{cc}{J}_{\alpha }& J_{t,nav}\\ 0& E_{nav}\end{array}\right]$

(29)$J_{\alpha }=\frac{\partial z_{\mathrm{t}}}{\partial {\mathsf{\theta }}_{\alpha }},J_{t,nav}=\frac{\partial z_{\mathrm{t}}}{\partial {\mathsf{\theta }}_{nav}},E_{nav}=\frac{\partial z_{nav}}{\partial {\mathsf{\theta }}_{nav}}$

Then, the CRLB of ${\mathsf{\theta }}_2$ will be

(30)$\mathrm{CRLB}\left({\mathsf{\theta }}_2\right)={\left\{{\left[\begin{array}{cc}{J}_{\alpha }& J_{t,nav}\\ 0& E_{nav}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{cc}{C}_t^{-1}& 0\\ 0& C_{nav}^{-1}\end{array}\right]\left[\begin{array}{cc}{J}_{\alpha }& J_{t,nav}\\ 0& E_{nav}\end{array}\right]\right\}}^{-1}={\left[\begin{array}{cc}{A}_2& U_2\\ V_2& D\end{array}\right]}^{-1}$

Hence, the CRLB of the emitter location ${\mathsf{\theta }}_{\alpha }$ can be expressed as

(31)$\mathrm{CRLB}\left({\mathsf{\theta }}_{\alpha }\right)=A_2^{-1}+A_2^{-1}{U}_2{\left(D-V_2{A}_2^{-1}{U}_2\right)}^{-1}{V}_2{A}_2^{-1}$

$A_2^{-1}$ is the CRLB of ${\mathsf{\theta }}_{\alpha }$ if the observer locations are known accurately. $A_2^{-1}{U}_2{\left(D-V_2{A}_2^{-1}{U}_2\right)}^{-1}{V}_2{A}_2^{-1}$ represents the performance degradation due to random errors of observer locations. $D$ is defined in (15), and $A_2$, $U_2$ and $V_2$ are defined as follows.

(32)$A_2=J_{\alpha }^{\mathrm{T}}{C}_t^{-1}{J}_{\alpha }={\displaystyle \sum _{k=1}^K{J}_{\alpha ,k}^{\mathrm{T}}{C}_{t,k}^{-1}{J}_{\alpha ,k}}$

(33)$U_2=J_{\alpha }^{\mathrm{T}}{C}_t^{-1}{J}_{t,nav}=diag\left\{J_{\alpha ,1}^{\mathrm{T}}{C}_{t,1}^{-1}{J}_{t,nav,1},\dots ,J_{\alpha ,K}^{\mathrm{T}}{C}_{t,K}^{-1}{J}_{t,nav,K}\right\}$

(34)$V_2=J_{t,nav}^{\mathrm{T}}{C}_t^{-1}{J}_{\alpha }=diag\left\{J_{t,nav,1}^{\mathrm{T}}{C}_{t,1}^{-1}{J}_{\alpha ,1},\dots ,J_{t,nav,K}^{\mathrm{T}}{C}_{t,K}^{-1}{J}_{\alpha ,K}\right\}$

The definition of $J_{t,nav,k}$ is given in (20), and $J_{\alpha ,\mathrm{k}}$ is defined as

(35)$J_{\alpha ,k}=\left[\frac{\partial t_{1,0,k}}{\partial {\mathsf{\theta }}_{\alpha }};\frac{\partial t_{2,0,k}}{\partial {\mathsf{\theta }}_{\alpha }};\dots ;\frac{\partial t_{M-1,0,k}}{\partial {\mathsf{\theta }}_{\alpha }}\right]$

(36)$\frac{\partial t_{j,0,k}}{\partial {\mathsf{\theta }}_{\alpha }}=\left[\frac{\partial t_{j,0,k}}{\partial u_k}\times \frac{\partial u_k}{\partial {\alpha }_k}\times \frac{\partial {\alpha }_k}{\partial {\alpha }_o},\frac{\partial t_{j,0,k}}{\partial u_k}\times \frac{\partial u_k}{\partial {\alpha }_k}\times \frac{\partial {\alpha }_k}{\partial {\dot{\alpha }}_o}\right]$

(37)$\frac{\partial t_{j,0,k}}{\partial u_k}=\frac{u_k^T-s_{j,k}^T}{c‖u_k-s_{j,k}‖}-\frac{u_k^T-s_{0,k}^T}{c‖u_k-s_{0,k}‖}$

(38)$\frac{\partial u_k}{\partial {\alpha }_k}=H\left({\alpha }_k\right)=\left[\begin{array}{cc}-\left(N+h\right)\cos B_k\sin L_k& \frac{\partial x_k}{\partial B_k}\\ -\left(N+h\right)\cos B_k\cos L_k& \frac{\partial y_k}{\partial B_k}\\ 0& \frac{\partial z_k}{\partial B_k}\end{array}\right]$

(39)$\frac{\partial x_k}{\partial B_k}=-\left(N+h\right)\sin B_k\cos L_k+\frac{a{e}^2\sin B_k{\cos }^2{B}_k\cos L_k}{{\left(\sqrt{1-e^2{\sin }^2{B}_k}\right)}^3}$

(40)$\frac{\partial y_k}{\partial B_k}=-\left(N+h\right)\sin B_k\sin L_k+\frac{a{e}^2\sin B_k{\cos }^2{B}_k\sin L_k}{{\left(\sqrt{1-e^2{\sin }^2{B}_k}\right)}^3}$

(41)$\frac{\partial z_k}{\partial B_k}=-\left(N\left(1-e^2\right)+h\right)\cos B_k+\frac{a{e}^2{\sin }^2{B}_k\cos B_k}{{\left(\sqrt{1-e^2{\sin }^2{B}_k}\right)}^3}$

(42)$\frac{\partial {\alpha }_k}{\partial {\alpha }_o}=I_{2\times 2},\frac{\partial {\alpha }_k}{\partial {\dot{\alpha }}_o}=\left(k-1\right)T\times I_{2\times 2}$

2.3. The Effect of System Errors

System errors are familiar in the scenario of passive geolocation based on SAGI networks, where accurate synchronization among observers is an impractical assumption. The system errors may include clock synchronization error, position bias of the observers, elevation bias of the target and non-horizontal velocity of the target. In this section, we analyze the effects of different system errors for the emitter modeled in the ECEF coordinate system.

Introducing Lagrange multiplier $\lambda$ to analyze the effect, the cost function is

(43)$\xi ={\left(𝔃-{𝔃}^o\left({\mathsf{\theta }}_1\right)\right)}^{\mathrm{T}}{C}^{-1}\left(𝔃-{𝔃}^o\left({\mathsf{\theta }}_1\right)\right)+{\lambda }^{\mathrm{T}}G\left({\mathsf{\theta }}_1\right)$

By using the perturbation method, substituting ${\mathsf{\theta }}_1$ with ${\mathsf{\theta }}_1+\Delta {\mathsf{\theta }}_1$ in (43) yields

(44)$\xi ={\left(n-J_1\Delta {\mathsf{\theta }}_1\right)}^{\mathrm{T}}{C}^{-1}\left(n-J_1\Delta {\mathsf{\theta }}_1\right)+{\lambda }^{\mathrm{T}}\left(G\left({\mathsf{\theta }}_1\right)+F^{\mathrm{T}}\Delta {\mathsf{\theta }}_1\right)$

Applying a partial derivative to (44) with respect to $\Delta {\mathsf{\theta }}_1$ and $\lambda$, we can obtain

(45)$\left[\begin{array}{cc}2J_1^T{C}^{-1}{J}_1& F\\ F^{\mathrm{T}}& 0\end{array}\right]\left[\begin{array}{c}\Delta {\mathsf{\theta }}_1\\ \lambda \end{array}\right]=\left[\begin{array}{c}2J_1^T{C}^{-1}n\\ -G\left({\mathsf{\theta }}_1\right)\end{array}\right]$

Solving (45) yields

(46)$\left[\begin{array}{c}\Delta {\mathsf{\theta }}_1\\ \mathsf{\theta }\end{array}\right]=\left[\begin{array}{c}2A_{11}{J}_1^T{C}^{-1}n-A_{12}G\left(\mathsf{\theta }\right)\\ 2A_{21}{J}_1^T{C}^{-1}n-A_{22}G\left(\mathsf{\theta }\right)\end{array}\right]$

where $A_{11}$ and $A_{12}$ are defined as

(47)$\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]={\left[\begin{array}{cc}2J_1^T{C}^{-1}{J}_1& F\\ F^{\mathrm{T}}& 0\end{array}\right]}^{-1}$

Hence, the corrected estimate of $\Delta \mathsf{\theta }$ is

(48)$\Delta {\mathsf{\theta }}_1=2A_{11}{J}_1^T{C}^{-1}n-A_{12}G\left({\mathsf{\theta }}_1\right)$

If there are no system errors, $E\left(n\right)=0$, $G\left({\mathsf{\theta }}_1\right)=0$. The localization accuracy in (48) is

(49)$\mathrm{MSE}\left({\mathsf{\theta }}_1\right)=4A_{11}{J}_1^T{C}^{-1}{J}_1{A}_{11}^{\mathrm{T}}$

It can be shown that (49) is equivalent to (25) when $I_1=J_1^T{C}^{-1}{J}_1$ is reversible. If system errors are non-negligible, $E\left(n\right)\ne 0$, $G\left({\mathsf{\theta }}_1\right)\ne 0$, the mean of $\Delta {\mathsf{\theta }}_1$ is biased. The localization accuracy will be degraded by the system errors.

2.3.1. System Error of Clock Synchronization

If there is only clock synchronization error, define $E\left({\mathsf{n}}_t\right)={\mathsf{\delta }}_t$; hence, $E\left(n\right)={\left[{\gamma }_t^T,0\right]}^T$, $G\left({\mathsf{\theta }}_1\right)=0$. Then, the localization accuracy is obtained from (48) as

(50)$\begin{array}{l}\mathrm{MSE}\left({\mathsf{\theta }}_1\right)=4A_{11}{J}_1^T{C}^{-1}{J}_1{A}_{11}^{\mathrm{T}}+4A_{11}{J}_1^T{C}^{-1}\left[\begin{array}{cc}{\gamma }_t{\gamma }_t^T& 0\\ 0& 0\end{array}\right]C^{-1}{J}_1{A}_{11}^{\mathrm{T}}\\ =4A_{11}{J}_1^T{C}^{-1}{J}_1{A}_{11}^{\mathrm{T}}+4A_{11}{B}_1{A}_{11}^{\mathrm{T}}\end{array}$

where $B_1$ is defined as

(51)$B_1=\left[\begin{array}{cc}{J}_u^T{C}_t^{-1}{\gamma }_t{\gamma }_t^T{C}_t^{-1}{J}_u& J_u^T{C}_t^{-1}{\gamma }_t{\gamma }_t^T{C}_t^{-1}{J}_{t,nav}\\ J_{t,nav}^T{C}_t^{-1}{\gamma }_t{\gamma }_t^T{C}_t^{-1}{J}_u& J_{t,nav}^T{C}_t^{-1}{\gamma }_t{\gamma }_t^T{C}_t^{-1}{J}_{t,nav}\end{array}\right]$

$4A_{11}{B}_1{A}_{11}^{\mathrm{T}}$ is the amount of degradation for localization accuracy introduced by clock synchronization system error.

2.3.2. System Error of Observer Location

When the system error of observer location is present, the mean of $n_{nav}$ is $E\left(n_{nav}\right)={\mathsf{\delta }}_{nav}$. Hence, the mean of $n$ is $E\left(n\right)={\left[0,{\gamma }_{nav}^T\right]}^T$ and $G\left({\mathsf{\theta }}_1\right)=0$. The localization accuracy is obtained from (48) as

(52)$\begin{array}{l}\mathrm{MSE}\left({\mathsf{\theta }}_1\right)=4A_{11}{J}_1^T{C}^{-1}{J}_1{A}_{11}^{\mathrm{T}}+4A_{11}{J}_1^T{C}^{-1}\left[\begin{array}{cc}0& 0\\ 0& {\gamma }_{nav}{\gamma }_{nav}^T\end{array}\right]C^{-1}{J}_1{A}_{11}^{\mathrm{T}}\\ =4A_{11}{J}_1^T{C}^{-1}{J}_1{A}_{11}^{\mathrm{T}}+4A_{11}{B}_2{A}_{11}^{\mathrm{T}}\end{array}$

where $B_2$ is defined as

(53)$B_2=\left[\begin{array}{cc}0& 0\\ 0& J_{nav}^T{C}_{nav}^{-1}{\gamma }_{nav}{\gamma }_{nav}^T{C}_{nav}^{-1}{J}_{nav}\end{array}\right]$

$4A_{11}{B}_2{A}_{11}^{\mathrm{T}}$ is the amount of degradation for localization accuracy introduced by the system error of observer location.

2.3.3. System Error of Target’s Elevation Bias

When the target’s elevation bias is present, the constrained conditions in (1) become $u_o^T{u}_o\ne {\left(R+h\right)}^2$. Define ${\delta }_h$ as the elevation bias; hence, the mean of $n$ is $E\left(n\right)=0$ and $G\left({\mathsf{\theta }}_1\right)={\left[{\delta }_h^2+2\left(N+h\right){\delta }_h,0\right]}^T$. The localization accuracy is obtained from (48) as

(54)$\mathrm{MSE}\left({\mathsf{\theta }}_1\right)=4A_{11}{J}_1^T{C}^{-1}{J}_1{A}_{11}^{\mathrm{T}}+A_{12}\left[\begin{array}{cc}{\left({\delta }_h^2+2\left(N+h\right){\delta }_h\right)}^2& 0\\ 0& 0\end{array}\right]A_{12}^{\mathrm{T}}$

$A_{12}\left[\begin{array}{cc}{\left({\delta }_h^2+2\left(N+h\right){\delta }_h\right)}^2& 0\\ 0& 0\end{array}\right]A_{12}^{\mathrm{T}}$ is the amount of degradation for localization accuracy introduced by the target’s elevation bias.

2.3.4. System Error of Non-Horizontal Velocity

When the motion of the target is non-horizontal, the second constrained condition in (1) is non-zero. Define ${\delta }_v$ as the bias of non-horizontal velocity; the second condition in (1) becomes $u_o^T{\dot{u}}_o={\delta }_v‖u_o‖={\delta }_v\left(N+h\right)$. The localization accuracy is obtained from (48) as

(55)$\mathrm{MSE}\left({\mathsf{\theta }}_1\right)=4A_{11}{J}_1^T{C}^{-1}{J}_1{A}_{11}^{\mathrm{T}}+A_{12}\left[\begin{array}{cc}0& 0\\ 0& {\delta }_v^2{\left(N+h\right)}^2\end{array}\right]A_{12}^{\mathrm{T}}$

$A_{12}\left[\begin{array}{cc}0& 0\\ 0& {\delta }_v^2{\left(N+h\right)}^2\end{array}\right]A_{12}^{\mathrm{T}}$ is the amount of degradation for localization accuracy introduced by the system error of non-horizontal velocity.

2.3.5. Total System Errors

If the above system errors exist in total, the mean of $n$ is $E\left(n\right)=\gamma$ and $G\left({\mathsf{\theta }}_1\right)=\delta$. The localization accuracy is obtained from (48) as

(56)$\begin{array}{c}\mathrm{MSE}\left({\mathsf{\theta }}_1\right)=4A_{11}{J}_1^T{C}^{-1}{J}_1{A}_{11}^{\mathrm{T}}+4A_{11}{J}_1^T{C}^{-1}\left(\gamma {\gamma }^T\right)C^{-1}{J}_1{A}_{11}^{\mathrm{T}}\\ +A_{12}\delta {\delta }^{\mathrm{T}}{A}_{12}^T-4A_{11}{J}_1^T{C}^{-1}\gamma {\delta }^{\mathrm{T}}{A}_{12}^T\end{array}$

3. Maximum Likelihood Estimator

Most previous research concentrates on the estimation of the target position with TDOA measurements. Due to the constrained conditions of target, previous algorithms are not suitable for jointly estimating the source position and velocity [29,30] with TDOA measurements only. In this section, we propose an iterative maximum likelihood estimator for both position and velocity in the scenario of geolocation based on SAGI networks.

Define the unknown parameter ${\mathsf{\theta }}_2={\left[{\mathsf{\theta }}_{\alpha }^T,{\mathsf{\theta }}_{nav}^T\right]}^T$ which contains the locations of target ${\mathsf{\theta }}_{\alpha }={\left[{\alpha }_o^T,{\dot{\alpha }}_0^T\right]}^T$ and observers ${\mathsf{\theta }}_{nav}=z_{nav}^o$. Suppose ${\mathsf{\theta }}_{2,g}={\left[{\mathsf{\theta }}_{\alpha ,g}^T,{\mathsf{\theta }}_{nav,g}^T\right]}^T$ is the initial solution. The observation equation can be approximated as

(57)$\begin{array}{l}\Delta z=z-z^o\\ =z-\left(z{|}_{{\mathsf{\theta }}_2={\mathsf{\theta }}_{2,g}}\right)-\left(\frac{\partial z}{\partial {\mathsf{\theta }}_2}{|}_{{\mathsf{\theta }}_2={\mathsf{\theta }}_{2,g}}\right)\left({\mathsf{\theta }}_2-{\mathsf{\theta }}_{2,g}\right)\end{array}$

The least squares solution of $\delta {\widehat{\mathsf{\theta }}}_2$ is given by

(58)$\delta {\widehat{\mathsf{\theta }}}_2={\left[{\left(\frac{\partial 𝔃}{\partial {\mathsf{\theta }}_2}{|}_{{\mathsf{\theta }}_2={\mathsf{\theta }}_{2,g}}\right)}^T{C}^{-1}\left(\frac{\partial 𝔃}{\partial {\mathsf{\theta }}_2}{|}_{{\mathsf{\theta }}_2={\mathsf{\theta }}_{2,g}}\right)\right]}^{-1}{\left(\frac{\partial 𝔃}{\partial {\mathsf{\theta }}_2}{|}_{{\mathsf{\theta }}_2={\mathsf{\theta }}_{2,g}}\right)}^T{C}^{-1}\left({\mathsf{\theta }}_2-{\mathsf{\theta }}_{2,g}\right)$

where $C$ is defined in (6), and $\frac{\partial z}{\partial {\mathsf{\theta }}_2}$ is given by (28). The corrected estimation of ${\mathsf{\theta }}_2$ is

(59)${\mathsf{\theta }}_2^{(1)}={\mathsf{\theta }}_{2,g}+\delta {\widehat{\mathsf{\theta }}}_2$

Substituting the corrected result into (57) by multiple steps of iteration, the estimation will converge and the amended variable $‖\delta {\widehat{\mathsf{\theta }}}_2‖$ will be smaller than the threshold. Then ${\mathsf{\theta }}_2^{(\mathrm{n})}$ is the final estimation.

The best localization accuracy of ${\widehat{\mathsf{\theta }}}_2$ is

(60)$MSE\left({\widehat{\mathsf{\theta }}}_2\right)={\left[{\left(\frac{\partial 𝔃}{\partial {\mathsf{\theta }}_2}\right)}^T{C}^{-1}\left(\frac{\partial 𝔃}{\partial {\mathsf{\theta }}_2}\right)\right]}^{-1}$

It can be easily shown that (60) is equivalent to (30). Hence, the accuracy of the proposed estimator can achieve the CRLB. The estimation of the target’s position and velocity in (59) is modeled in the geodetic coordinate system. The transformation matrix from geodetic to ECEF is defined in (26).

4. Discussion

We consider the SAGI passive localization geometry used in [23]. The observer network contains GEO satellite, LEO satellite, airplane and ground station. The initial positions are (105°, 0°, 36,000 km), (101°, 0°, 500 km), (103°, 3°, 10 km) and (105°, 5°, 0 km), respectively. The initial velocity and course angle of the airplane are 2000 km/h and 60°. The initial locations of these observers under the ECEF coordinate system are shown in Table 1.

In GDOP (geometric dilution of precision) figures, the X-axis is the longitude and the Y-axis is the latitude. The curve in the GDOP figure means the localization precision of targets, whose location is on the curve. The red hexagram denotes the substar of the GEO, the red five-pointed star denotes the substar of the LEO, the red rhombus denotes the sub-point of the airplane, and the red circle denotes the location of the ground station.

Scenario 1: The sensor network contains the GEO satellite and the LEO satellite. Simulation results are shown in Figure 1, Figure 2 and Figure 3. Figure 1 places emphasis on the localization performance of total sub-satellite area, while Figure 2 and Figure 3 places emphasis on the effects of different system errors for the target fixed on a point. The random errors of satellite positions are (100, 100, 100) m which represents the covariance matrix $Q_{s_{j,k}}=diag\left\{{\left(100 \mathrm{m}\right)}^2,{\left(100 \mathrm{m}\right)}^2,{\left(100 \mathrm{m}\right)}^2\right\}$. The emitter is on the surface of Earth with zero height. The detection number $K$ is 10. The Monte Carlo number of the ML estimation is 1000.

The random error of TDOA measurement is 500 ns in Figure 1, which represents the TDOA covariance matrix $Q_{t_{j,k}}=diag\left\{{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2,\dots ,{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2\right\}$. Figure 1a presents the GDOP without any system error. Figure 1b shows the GDOP with 1 us bias of clock synchronization. Figure 1c gives the GDOP with satellite position bias of (2000, 2000, 2000) m. Figure 1d shows the GDOP with 2 km elevation bias of the target. From the simulation results in Figure 1, we can make the following conclusions: (1) The localization performance degrades seriously in the area around the moving direction of the LEO satellite. (2) The effects of different system errors are also diverse. Comparing the four GDOPs in Figure 1, we can observe that clock synchronization bias is the most influential factor for this scenario. The elevation bias of the target is the second most influential factor. The last one is position bias of the satellites.

In Figure 2, the random errors of TDOA measurements are variable as ${\sigma }_t$ ns which means the TDOA covariance matrix is $Q_{t_{j,k}}=diag\left\{{\left(\frac{{\sigma }_t}{{10}^9}\right)}^2,\dots ,{\left(\frac{{\sigma }_t}{{10}^9}\right)}^2\right\}$. The emitter location is fixed on $\left(95°,-5°,0 \mathrm{km}\right)$. Figure 2 compares theoretic performances of geolocation with different system errors. “Theoretical accuracy without system error” is generated by (25), “Theoretical accuracy with clock synchronization error” is generated by (50), “Theoretical accuracy with observer position bias” is generated by (52), “Theoretical accuracy with target elevation bias” is generated by (54) and “ML localization accuracy without system error” is generated by (58).

The effects of different system errors are clearly shown in Figure 2. We can observe that clock synchronization bias is the most influential factor. After that are target’s elevation bias and satellites’ position bias. The simulation results are similar to those in Figure 1. The localization accuracy of the proposed ML estimator is also provided in Figure 2 for the case of no system error. It can be seen from Figure 2 that localization accuracy of the proposed ML estimator agrees with the CRLB very well when TDOA noise ${\sigma }_t$ is less than 0.7 us.

Figure 3 compares localization accuracy using the GEO and LEO satellites by varying clock synchronization error. The random error of TDOA measurement is 500 ns, which represents the TDOA covariance matrix $Q_{t_{j,k}}=diag\left\{{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2,\dots ,{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2\right\}$. Other errors are assumed nonexistent. The emitter location is fixed on $\left(95°,-5°,0 \mathrm{km}\right)$. The effect of clock synchronization error is clearly shown in Figure 3. The localization performance of the proposed ML estimator agrees with the theoretical accuracy very well.

Scenario 2: The sensor network contains the LEO satellite and the airplane. Simulation results are shown in Figure 4, Figure 5 and Figure 6. The random errors of observers’ positions are (100, 100, 100) m which represents the covariance matrix $Q_{s_{j,k}}=diag\left\{{\left(100 \mathrm{m}\right)}^2,{\left(100 \mathrm{m}\right)}^2,{\left(100 \mathrm{m}\right)}^2\right\}$. The target is on the surface of Earth with zero height. The detection number $K$ is 10. The Monte Carlo number of the ML estimation is 1000.

The random error of TDOA measurement is 500 ns in Figure 4, which represents the TDOA covariance matrix $Q_{t_{j,k}}=diag\left\{{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2,\dots ,{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2\right\}$. Figure 4a presents the GDOP without any system error. Figure 4b shows the GDOP with 1 us bias of clock synchronization. Figure 4c gives the GDOP with satellite position bias of (2000, 2000, 2000) m. Figure 4d shows the GDOP with 2 km elevation bias of the target. From the simulation results in Figure 3, we can observe the following conclusions: (1) The localization accuracy in the area around the sub-satellite point is more precise than the area around the sub-airplane point. (2) There are two ambiguous location areas around the sub-airplane point. (3) The effects of different system errors are also diverse.

In Figure 5, the random errors of TDOA measurements are variable as ${\sigma }_t$ ns which means the TDOA covariance matrix is $Q_{t_{j,k}}=diag\left\{{\left(\frac{{\sigma }_t}{{10}^9}\right)}^2,\dots ,{\left(\frac{{\sigma }_t}{{10}^9}\right)}^2\right\}$. The target’s location is fixed on $\left(101°,3°,0 \mathrm{km}\right)$. “Theoretical accuracy without system error” is generated by (25), “Theoretical accuracy with clock synchronization error” is generated by (50), “Theoretical accuracy with observer position bias” is generated by (52), “Theoretical accuracy with target elevation bias” is generated by (54) and “ML localization accuracy without system error” is generated by (58).

We can observe from Figure 5 that clock synchronization bias is the most influential factor. After that are elevation bias and position bias of the satellites. The simulation results are similar to those in Figure 2. The localization accuracy of the proposed ML estimator is also provided in Figure 5 for the case of no system error. It can be seen from Figure 5 that localization accuracy of the proposed ML estimator agrees with the CRLB very well when TDOA noise ${\sigma }_t$ is less than 0.5 us.

Figure 6 compares localization accuracy using the LEO and airplane by varying observer position bias. The random errors of TDOA measurements are 500 ns, which represents the TDOA covariance matrix $Q_{t_{j,k}}=diag\left\{{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2,\dots ,{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2\right\}$. Other errors are assumed nonexistent. The emitter location is fixed on $\left(101°,3°,0 \mathrm{km}\right)$. The effect of observer position bias is clearly shown in Figure 6. The localization performance of the proposed ML estimator agree with the theoretical accuracy very well.

Scenario 3: The sensor network contains the LEO satellite, airplane and ground station. Simulation results are shown in Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12. The random errors of satellite positions are (100, 100, 100) m which represents the covariance matrix $Q_{s_{j,k}}=diag\left\{{\left(100 \mathrm{m}\right)}^2,{\left(100 \mathrm{m}\right)}^2,{\left(100 \mathrm{m}\right)}^2\right\}$. The target is moving on the surface of the Earth with a constant speed of 300 m/s. The course angle of the target is 30°. The height of the target is 10 km. The detection number $K$ is 50. The Monte Carlo number of the ML estimation is 1000.

The random errors of TDOA measurement are 500 ns in Figure 7 and Figure 8, which represent the TDOA covariance matrix $Q_{t_{j,k}}=diag\left\{{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2,\dots ,{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2\right\}$. Figure 7 and Figure 8 show the GDOPs of position estimation and velocity estimation, respectively.

Figure 7a presents the GDOP of position estimation with 1 us bias of clock synchronization. Figure 7b shows the GDOP of position estimation with observers’ position bias of (2000, 2000, 2000) m. Figure 7c gives the GDOP of position estimation with 2000 m elevation bias of the target. Figure 7d shows the GDOP of position estimation when the non-horizontal velocity of the target is 20 m/s. From the simulation results in Figure 7, we can make the following conclusions: (1) The ambiguous location area is around the line between the sub-LEO point and the sub-airplane point. (2) The effects of different system errors are also diverse. We can observe that elevation bias of the target is the most influential factor for this scenario. The position bias of the observers is the second influential factor. The effects of clock synchronization bias and non-horizontal velocity of the target are not serious.

Figure 8a presents the GDOP of velocity estimation with 1 us bias of clock synchronization. Figure 8b shows the GDOP of velocity estimation with observers’ position bias of (2000, 2000, 2000) m. Figure 8c gives the GDOP of velocity estimation with 2000 m elevation bias of the target. Figure 8d shows the GDOP of velocity estimation when the non-horizontal velocity of the target is 20 m/s. From the simulation results in Figure 8, we can make the following conclusions: (1) The ambiguous area of velocity estimation is around the line between the sub-LEO point and the sub-airplane point. (2) The effects of different system errors are also diverse. We can observe that non-horizontal velocity of the target is the most influential factor for this scenario. This is different from the conclusion in Figure 7 where non-horizontal velocity is the lowest influential factor. The differences in clock synchronization bias and position bias are difficult to distinguish from Figure 8a,b.

In Figure 9 and Figure 10, the random errors of TDOA measurements are variable as ${\sigma }_t$ ns which means the TDOA covariance matrix is $Q_{t_{j,k}}=diag\left\{{\left(\frac{{\sigma }_t}{{10}^9}\right)}^2,\dots ,{\left(\frac{{\sigma }_t}{{10}^9}\right)}^2\right\}$. The emitter initial position is $\left(101°,3°,10 \mathrm{km}\right)$. “Theoretical accuracy without system error” is generated by (25), “Theoretical accuracy with clock synchronization error” is generated by (50), “Theoretical accuracy with observer position bias” is generated by (52), “Theoretical accuracy with target elevation bias” is generated by (54), “Theoretical accuracy with target non-horizontal velocity bias” is generated by (55) and “ML localization accuracy without system error” is generated by (58).

Figure 9 presents the theoretical performances of geolocation with different system errors. We can observe from Figure 9 that elevation bias of the target is the most influential factor. After that are position bias of observers and clock synchronization bias. Due to the negligible effect, the target’s non-horizontal velocity can be ignored for geolocation. The localization accuracy of the proposed ML estimator is also provided in Figure 9 for the case of no system error. It can be seen from Figure 9 that localization accuracy of the proposed ML estimator agrees with the CRLB very well when TDOA noise ${\sigma }_t$ is less than 0.9 us.

Figure 10 shows the theoretic performances of velocity estimation with different system errors. We can observe from Figure 10 that the target’s non-horizontal velocity is the most influential factor, which is different from the conclusion in Figure 9. After that are clock synchronization bias, observers’ position bias and target’s elevation bias, respectively. The velocity estimation accuracy of the proposed ML estimator is also provided in Figure 10 for the case of no system error. It can be seen from Figure 10 that velocity estimation accuracy of the proposed ML estimator agrees with the CRLB very well when TDOA noise ${\sigma }_t$ is less than 0.7 us.

In Figure 11 and Figure 12, the random errors of TDOA measurements are 500 ns, which represent the TDOA covariance matrix $Q_{t_{j,k}}=diag\left\{{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2,\dots ,{\left(\frac{500 \mathrm{ns}}{{10}^9}\right)}^2\right\}$. The emitter initial position is $\left(101°,3°,10 \mathrm{km}\right)$.

Figure 11 compares localization accuracy by varying the target elevation bias. Only TDOA measurement error and target elevation bias are present. The effect of target elevation bias is clearly shown in Figure 11. The localization performance of the proposed ML estimator agrees with the theoretical accuracy very well.

Figure 12 compares localization accuracy by varying the non-horizontal velocity bias of the target. Only TDOA measurement error and non-horizontal velocity bias are present. The effect of non-horizontal velocity bias is clearly shown in Figure 12. The localization performance of the proposed ML estimator agrees with the theoretical accuracy very well.

5. Conclusions

This paper investigated the method of geolocation and tracking by TDOA measurements based on SAGI networks. The CRLB for target localization and tracking accuracy was derived in different coordinate systems. Then, the effects of different system errors were analyzed, such as clock synchronization error, position bias of observers, elevation bias of target and non-horizontal velocity of target. Further, a ML estimator for the target’s position and velocity was proposed to compare with theoretical performance. Finally, extensive numerical and experimental results were presented to validate the theoretical analyses.

Author Contributions

Conceptualization, J.L.; methodology, S.L. (Shouye Lv); software, C.W.; validation, Y.L. and S.L. (Shuai Liao); writing—original draft preparation, Y.J.; writing—review and editing, J.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The locations of these observers are generated by computer simulation. It is easy to generate simulation with the method in the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. GDOPs of position estimation using the GEO and LEO satellites.

Figure 2. Localization accuracy comparison using the GEO and LEO satellites by varying TDOA measurement noise.

Figure 3. Localization accuracy comparison using the GEO and LEO satellites by varying clock synchronization error.

Figure 4. GDOPs of position estimation using the LEO and airplane.

Figure 5. Localization accuracy comparison using the LEO and airplane by varying TDOA measurement noise.

Figure 6. Localization accuracy comparison using the LEO and airplane by varying observer position bias.

Figure 7. GDOPs of position estimation for the scenario of using LEO satellite, airplane and ground station.

Figure 8. GDOPs of velocity estimation for the scenario of using LEO satellite, airplane and ground station.

Figure 9. Localization accuracy comparison using LEO satellite, airplane and ground station by varying TDOA measurement noise.

Figure 10. Comparison of velocity estimation using LEO satellite, airplane and ground station by varying TDOA measurement noise.

Figure 11. Localization accuracy comparison using LEO satellite, airplane and ground station by varying the target elevation bias.

Figure 12. Comparison of velocity estimation using LEO satellite, airplane and ground station by varying non-horizontal velocity bias.

Table 1

The initial locations of the observers in ECEF.

Title 1	X (m)	Y (m)	Z (m)	Vx (m/s)	Vy (m/s)	Vz (m/s)
GEO satellite	−10,912,884.14	40,727,438.07	−44,151.31	0	0	0
LEO satellite	−1,306,917.55	6,752,827.56	−7416.31	−4793.46	−933.62	−5380.76
Airplane	−1,435,061.91	6,215,936.03	332,097.67	−264.62	−86.94	480.47
Ground station	−1,644,543.41	6,137,519.57	552,183.96	0	0	0

References

1. Tzoreff, E.; Bobrovsky, B.Z.; Weiss, A.J. Single receiver emitter geolocation based on signal periodicity with oscillator instability. IEEE Trans. Signal Process.; 2014; 62, pp. 1377-1385. [DOI: https://dx.doi.org/10.1109/TSP.2014.2298831]

2. Guoqing, Q.; Yinya, L.; Andong, S.; Xinying, G. An investigation on observer minimum circumnavigation velocity for bearings-only target tracking. IEEE Access; 2020; 8, pp. 126805-126813. [DOI: https://dx.doi.org/10.1109/ACCESS.2020.3006034]

3. Xiang, F.; Wang, J.; Yuan, X. Research on passive detection and location by fixed single observer. Proceedings of the 2020 International Conference on Information Science, Parallel and Distributed Systems (ISPDS); Xi’an, China, 14–16 August 2020; pp. 35-39.

4. Geeraert, J.L.; McMahon, J.W. Dual-satellite geolocation with ephemeris correction and uncertainty mapping. IEEE Trans. Aerosp. Electron. Syst.; 2020; 56, pp. 723-735. [DOI: https://dx.doi.org/10.1109/TAES.2019.2920045]

5. Fan, L.; Zhou, C. 3D Lightning Location Method Based on Range Difference Space Projection. Remote Sens.; 2022; 14, 1003. [DOI: https://dx.doi.org/10.3390/rs14041003]

6. Torney, D.C. Localization and Observability of Aircraft via Doppler Shifts. IEEE Trans. Aerosp. Electron. Syst.; 2007; 43, pp. 1163-1168. [DOI: https://dx.doi.org/10.1109/TAES.2007.4383606]

7. Li, J.; He, Y.; Zhang, X.; Wu, Q. Simultaneous Localization of Multiple Unknown Emitters Based on UAV Monitoring Big Data. IEEE Trans. Ind. Inform.; 2021; 17, pp. 6303-6313. [DOI: https://dx.doi.org/10.1109/TII.2020.3048987]

8. Pattison, T.; Chou, S.I. Sensitivity analysis of dual-satellite geolocation. IEEE Trans. Aerosp. Electron. Syst.; 2000; 36, pp. 56-71. [DOI: https://dx.doi.org/10.1109/7.826312]

9. Wang, Y.; Zhou, T.; Yi, W.; Kong, L. A GDOP-Based Performance Description of TOA Localization with Uncertain Measurements. Remote Sens.; 2022; 14, 910. [DOI: https://dx.doi.org/10.3390/rs14040910]

10. Chan, Y.T.; Ho, K.C. A simple and efficient estimator for hyperbolic location. IEEE Trans. Signal Process.; 1994; 42, pp. 1905-1915. [DOI: https://dx.doi.org/10.1109/78.301830]

11. Vidal-Valladares, M.G.; Díaz, M.A. A Femto-Satellite Localization Method Based on TDOA and AOA Using Two CubeSats. Remote Sens.; 2022; 14, 1101. [DOI: https://dx.doi.org/10.3390/rs14051101]

12. Wang, Y.; Sun, G.-C.; Wang, Y.; Yang, J.; Zhang, Z.; Xing, M. A Multi-Pulse Cross Ambiguity Function for the Wideband TDOA and FDOA to Locate an Emitter Passively. Remote Sens.; 2022; 14, 3545. [DOI: https://dx.doi.org/10.3390/rs14153545]

13. Wang, G.; Ho, K.C.; Chen, X. Bias Reduced Semidefinite Relaxation Method for 3-D Rigid Body Localization Using AOA. IEEE Trans. On Signal Process.; 2021; 69, pp. 3415-3430. [DOI: https://dx.doi.org/10.1109/TSP.2021.3086360]

14. Xu, S.; Dogancy, K. Optimal Sensor Placement for 3-D Angle-of-Arrival Target Localization. IEEE Trans. Aerosp. Electron. Syst.; 2017; 53, pp. 1196-1211. [DOI: https://dx.doi.org/10.1109/TAES.2017.2667999]

15. Foy, W.H. Position location solutions by Taylor-series estimation. IEEE Trans. Aerosp. Electron. Syst.; 1976; 12, pp. 187-194. [DOI: https://dx.doi.org/10.1109/TAES.1976.308294]

16. Mikhalev, A.; Ormondroyd, R.F. Comparison of hough transform and particle filter methods of emitter geolocation using fusion of TDOA data. Proceedings of the 2007 4th Workshop on Positioning, Navigation and Communication; Hannover, Germany, 22 March 2007; pp. 121-127.

17. Bin, Y.Z.; Yan, Q.; Nan, L.A. PSO based passive satellite localization using TDOA and FDOA measurements. Proceedings of the 2011 10th IEEE/ACIS International Conference on Computer and Information Science; Sanya, China, 18 May 2011; pp. 251-254.

18. Wang, G.; Li, Y.; Ansari, N. A semidefinite relaxation method for source localization using TDOA and FDOA measurements. IEEE Trans. Veh. Technol.; 2013; 62, pp. 853-862. [DOI: https://dx.doi.org/10.1109/TVT.2012.2225074]

19. Liu, J.; Shi, Y.; Fadlullah, Z.M.; Kato, N. Space-air-ground integrated network: A survey. IEEE Commun. Surv. Tutor.; 2018; 20, pp. 2714-2741. [DOI: https://dx.doi.org/10.1109/COMST.2018.2841996]

20. Wang, Y.; Xu, Y. Hybrid satellite-aerial-terrestrial networks in emergency scenarios: A survey. China Commun.; 2017; 14, pp. 1-13. [DOI: https://dx.doi.org/10.1109/CC.2017.8010971]

21. Zhou, S.; Wang, G.; Zhang, S.; Niu, Z.; Shen, X.S. Bidirectional mission offloading for agile space-air-ground integrated networks. IEEE Trans. Wirel. Commun.; 2019; 26, pp. 38-45. [DOI: https://dx.doi.org/10.1109/MWC.2019.1800290]

22. Li, W.; Liu, Y.; Li, X.; Shen, Y. Three-dimensional cooperative localization via space-air-ground integrated networks. China Commun.; 2022; 19, pp. 253-263. [DOI: https://dx.doi.org/10.23919/JCC.2022.01.018]

23. Li, J.; Lv, S.; Liu, Y.; Zhang, K.; Yang, L.; He, Y. Performance analysis for space-air-ground integrated passive localization using TDOA measurements. Proceedings of the 2022 IEEE International Conference on Unmanned Systems; Guangzhou, China, 28–30 October 2022; pp. 199-203.

24. Ho, K.C.; Chan, Y.T. Solution and performance analysis of geolocation by TDOA. IEEE Trans. Aerosp. Electron. Syst.; 1993; 29, pp. 1311-1322. [DOI: https://dx.doi.org/10.1109/7.259534]

25. Hmam, H. Optimal sensor velocity configuration for TDOA-FDOA geolocation. IEEE Trans. Signal Process.; 2017; 65, pp. 628-637. [DOI: https://dx.doi.org/10.1109/TSP.2016.2621724]

26. Ho, K.C.; Lu, X.; Kovavisaruch, L. Source localization usingTDOA and FDOA measurements in the presence of receiver locationerrors: Analysis and solution. IEEE Trans. Signal Process.; 2007; 55, pp. 684-696. [DOI: https://dx.doi.org/10.1109/TSP.2006.885744]

27. Ho, K.C.; Yang, L. On the use of a calibration emitter for source localization in the presence of sensor position uncertainty. IEEE Trans. Signal Process.; 2008; 56, pp. 5758-5772. [DOI: https://dx.doi.org/10.1109/TSP.2008.929870]

28. Yang, L.; Ho, K.C. Alleviating Sensor Position Error in Source Localization Using Calibration Emitters at Inaccurate Locations. IEEE Trans. Signal Process.; 2010; 58, pp. 67-83. [DOI: https://dx.doi.org/10.1109/TSP.2009.2028947]

29. Ho, K.C.; Chan, Y.T. Geolocation of a known altitude object from TDOA and FDOA measurements. IEEE Trans. Aerosp. Electron. Syst.; 1997; 33, pp. 770-783. [DOI: https://dx.doi.org/10.1109/7.599239]

30. Huang, S.; Pun, Y.M.; So, A.M.C.; Yang, K. A provably convergent projected gradient-type algorithm for TDOA-based geolocation under the quasi-parabolic lonosphere model. IEEE Signal Process. Lett.; 2020; 27, pp. 1335-1339. [DOI: https://dx.doi.org/10.1109/LSP.2020.3010676]

31. Marzetta, T.L. Simple derivation of the constrained multiple parameter Cramer-Rao bound. IEEE Trans. Acoust. Speech Signal Process.; 1993; 41, pp. 2247-2249. [DOI: https://dx.doi.org/10.1109/78.218151]