A Two-Stage Aerial Target Localization Method Using Time-Difference-of-Arrival Measurements with the Minimum Number of Radars

Abstract

Distributed radar systems promise to significantly enhance target localization by virtue of the superiority of multi-view observations from widely separated radars, compared to their monostatic counterparts. Nevertheless, when the radar number is limited, performing target localization bears the brunt of the parameter identifiability requirement that the parameter number must be no less than the number of independent measurements. In this way, the canonical two-stage target localization method, as well as its developments, is no longer appropriate for direct application. Hence, in this paper, we propose a novel target localization method using time-difference-of-arrival (TDOA) measurements with the minimum number of radars under platform position uncertainties. The referred distributed system is a bistatic multi-receiver system, where the primary signal is transmitted by a geostationary Earth orbit (GEO) satellite while receivers are equipped on several unmanned aerial vehicles (UAVs). In the first stage, the reference range from the reference radar to the target is estimated by a quadratic function, and then the weighted least squares (WLS) solution of the target location is updated by substituting the range estimate back into it. In the second stage, we invoke the Taylor series approximation to further refine the target localization obtained by the first stage. It can be foreseen that the developed method is beneficial for scenarios with a limited number of radars, including engineering projects such as fire control, surveillance, and guidance, to support high-accuracy target localization. The simulation results show the superiority of the localization performance of the proposed method over other existing methods.

1. Introduction

Target localization with spatially distributed radars has drawn more and more attention to provide a location for the effective implementation of information-snooping electronic countermeasures (ECM) and even-final precision strikes in the background of modern warfare [1,2,3]. Generally, target localization with multiple radars can be categorized into direct methods and indirect methods. For the former type, target parameters are directly estimated from received echo signals [4,5,6]. For the latter type, the target measurements are extracted from the returns followed by solving linearized equations to acquire those parameters [7,8,9,10]. Although the direct approach seems to have a more intuitive formulation, solving the direct localization problem is more difficult. Moreover, the direct approach is vulnerable to parametric model mismatch and bears a heavy computational burden. Therefore, we resort to the indirect approach for our work in this paper.

Similar to the application of satellite navigation and positioning [11,12], the key to both problems lies in determining the location of the target of interest using received target echoes or obtained target measurements. However, satellite navigation and positioning are different from target localization in that the relationship between the radar system and the object is cooperative. In this case, the ground observation station sends the precise satellite ephemerides to the user, and then the user calculates its current location itself by combining the ephemerides and the received echoes/measurements, or by saying that the user obtains its 3D coordinates in an active way. For the target localization, this relationship is noncooperative. The target is detected first and then localized using the received data at the processing center in a passive manner.

Before performing target localization, whether or not target parameters are identifiable is a required prerequisite. If time-of-arrival (TOA) measurements are employed to localize an aerial target, there are at least three radars required excluding pathological configurations. On the other hand, if time-difference-of-arrival (TDOA) measurements are used, the number should be no less than four [2]. Furthermore, it is worthwhile to note that the parameter identifiability condition does not necessarily satisfy the localization accuracy. If one stipulates a localization accuracy, the radar network should be optimized involving more radars [13,14,15].

The TDOA target localization problem has been intensively investigated in the fields of communication, radar, and sonar. Although the measurements obtained are distinct, their localization processes are analogous. In early times, Fang et al. developed a closed-form solution [16] for when the number of TDOA measurements is equal to the target parameter number, whereby an additional piece of information is invoked. The more general case was researched by Smith [7], Friedlander [17], and Schau [18]. However, their developments were not optimum. With an initial guess, the Taylor series method [19,20] is a good candidate, but its estimates suffer from local minima once the initial guess deviates from the true values to some extent. The two-stage hyperbolic location estimator [21] proposed by Chan and Ho is one of the most efficient methods for target localization using TDOA measurements. In the first stage, the target location, as well as the reference range, i.e., the range from the reference radar to the target, is collected as a parameter vector to be estimated, and the arranged parameter vector is resolved by transforming the group of nonlinear equations into linear equations. In the second stage, the target location can be refined by a combination of the estimated reference range and the obtained statistics of target location estimates in the first stage. Although the case of three sensors was considered in this study, the sensor positions were exactly known. To take the platform position uncertainty into account, at the 2004 International Symposium on Circuits and Systems (ISCAS) conference, Ho and Parith further improved the two-stage method to dispose of the problem related to erroneous positions [22]. Accordingly, a two-stage target localization method using both TDOA and frequency-difference-of-arrival (FDOA) measurements [23], where positions are imprecisely known, was developed by Ho and Lu to improve target location accuracy to a certain extent, in addition to the extra ability to acquire the target velocity estimation using FDOA measurements. However, the introduction of the stated nonlinear relationships led to a notable degradation of the localization accuracy with the increase in platform position errors. Moreover, the authors of [24,25] applied the first-order Taylor series approximation process, instead of the refinement in the second stage, to avoid a nonlinear operation for the target location. As a result, the localization accuracy was ameliorated even with large position errors. Another improved method was proposed in [26], dedicated to improving target localization accuracy by directly exploiting the connection between the reference range and the target location in the first stage. However, it was designed for the case of absent position errors, while also lacking the target location refinement in the second stage. Furthermore, TDOA localization for locating multiple disjoint sources was studied in [27], where the fact that TDOA measurements from different targets produce the same position displacements was exploited to improve performance. Furthermore, Sun dealt with this problem using both TDOA and FDOA measurements in [28], while Lu improved the technique in [29] via new second-stage processing. To handle robust target localization, Cheung and Ma established a network of sensors within the framework of maximum likelihood estimation (MLE) [30]. By converting the original nonconvex problems into convex ones, the SDP solver was employed to yield a satisfying solution. The authors of [31,32] further tackled this problem from the perspective of convex optimization, and the estimation results matched the optimal performance. Additionally, in a similar way, TDOA-based target localization was extensively investigated in distributed multiple-input multiple-output (MIMO) systems by Amiri and Behnia [33,34,35,36], and they also extended their research to the case of clock synchronization errors [37].

In summary, the stated localization techniques involve a variety of system perturbations in practical scenarios. However, the fact that the number of deployed radars is limited was ignored during their design. Actually, the radar number is extraordinarily valuable, especially for some military applications. Therefore, it is desired to achieve the required localization accuracy by reducing the radar number as much as possible. In this regard, we find out that the state-of-the-art two-stage methods add the reference range into the parameter vector of interest; thus, the required minimum number of deployed radars decreases to five, which brings about an additional burden of required radars. Therefore, in this paper, a novel two-stage aerial target localization method using TDOA measurements with a minimum number of receiving radars is proposed. By exploiting the relationship between the reference range and target location, as well as a priori statistics of platform position errors, we update the reference range and the target location in an iterative way to avoid the estimation of the reference range. As a result, the required minimum number of radars is reduced. Lastly, the effectiveness of our proposed methods is verified by simulation results.

The remaining sections are structured as follows: in Section 2, the multi-radar localization model and the proposed two-stage method using TDOA measurements are discussed in detail; numerical results are presented in Section 3 to verify the effectiveness of the proposed method, and the corresponding advantages are summarized in Section 4; lastly, we draw a brief conclusion in Section 5.

Notation: In this paper, vectors are denoted by boldface lowercase letters and matrices are represented by boldface uppercase letters. Superscripts T and H on a matrix or a vector denote the transpose and Hermitian transpose operations, respectively. $E\{\cdot \}$ is used to denote the statistical expectation operator, and ${()}^{-1}$ the inverse operator. $I_m$ is an $m\times m$ identity matrix, while $1_{m\times 1}$ is a vector of all ones. Lastly, $\left|\right|\cdot \left|\right|$ represents the $l_2$ Euclidean norm.

2. Materials and Methods

2.1. Basic Signal Model

Consider a space–air integrated distributed radar system composed of a geostationary Earth orbit (GEO) spaceborne radar, served as a transmitter and K unmanned aerial vehicle (UAV) airborne radars as receivers, as shown in Figure 1. When the measurement error and the satellite orbit error are considered, we can obtain the GEO satellite navigation position $\overline{s}={[{\overline{s}}_x,{\overline{s}}_y,{\overline{s}}_z]}^{\mathrm{T}}$ by invoking the imprecise information from the Global Navigation Satellite System (GNSS). However, the actual position $s={[s_x,s_y,s_z]}^{\mathrm{T}}$ is unknown and is expressed by $s=\overline{s}+\mathsf{\Delta }s$, where $\mathsf{\Delta }s={[\mathsf{\Delta }{s}_x,\mathsf{\Delta }{s}_y,\mathsf{\Delta }{s}_z]}^{\mathrm{T}}$ represents the three-dimensional (3D) GEO satellite position error. Without loss of generality, $\mathsf{\Delta }s$ obeys Gaussian distribution with zero mean and the following covariance matrix:

$$Q_{\mathrm{s}}=E\left\{\mathsf{\Delta }s\mathsf{\Delta }{s}^{\mathrm{T}}\right\}.$$

(1)

Likewise, the actual position of the k-th UAV position $u_k={[u_{xk},u_{yk},u_{zk}]}^{\mathrm{T}}$ is perturbed from the navigation position ${\overline{u}}_k={[{\overline{u}}_{xk},{\overline{u}}_{yk},{\overline{u}}_{zk}]}^{\mathrm{T}}$ by its 3D platform position error $\mathsf{\Delta }{u}_k={[\mathsf{\Delta }{u}_{xk},\mathsf{\Delta }{u}_{yk},\mathsf{\Delta }{u}_{zk}]}^{\mathrm{T}}$ due to the measurement uncertainty and nonstationary flight motion. Here, position errors among different UAVs are assumed uncorrelated with each other, and the k-th UAV position error is Gaussian distributed with zero mean and the following covariance matrix:

$$Q_{\mathrm{u}}=E\left\{\mathsf{\Delta }u\mathsf{\Delta }{u}^{\mathrm{T}}\right\},$$

(2)

where $\mathsf{\Delta }u={[\mathsf{\Delta }{u}_1^{\mathrm{T}},\mathsf{\Delta }{u}_2^{\mathrm{T}},\cdots ,\mathsf{\Delta }{u}_K^{\mathrm{T}}]}^{\mathrm{T}}$.

To cater to a more practical localization scenario, the range measurements from the GEO satellite to an aerial target of interest, to UAVs are gathered to perform the target localization, with the assistance of the prior statistics of position errors, in our designed framework. In view of the relationship between the TDOA and the range difference of arrival (RDOA), $\tau =r/c$, where $c$ is the speed of light. It is worthwhile to note that we interchange the terms TDOA and RDOA hereafter for their equivalence. The aerial target to be located is at the coordinate $p={[x_p,y_p,z_p]}^{\mathrm{T}}$; then, the bistatic range between the navigation position from the GEO satellite and the k-th UAV navigation position is expressed as

$$r_{Tk}^o=r_T^o+r_k^o={‖\overline{s}-p‖}_2+{‖{\overline{u}}_k-p‖}_2.$$

(3)

Furthermore, the corresponding bistatic pseudo-range for the transmit/receive pair $(T,k)$, where T denotes the transmitter, is given by

$$r_{Tk}={\overline{r}}_T+r_k,$$

(4)

where ${\overline{r}}_T=\left|\right|s-p{\left|\right|}_2$ represents the noiseless range from the GEO satellite to the target, and $r_k={\overline{r}}_k+n_k$ indicates the noisy range measurement from the target to the k-th receiving radar; ${\overline{r}}_k=\left|\right|u_k-p{\left|\right|}_2$ is the noiseless range, and $n_k$ is the measurement noise. By differencing the range measurement $r_{Tk}$ to the reference one, i.e., $r_{T1}$, we can obtain the k-th TDOA measurement as follows:

$$r_{k1}=r_{Tk}-r_{T1}={‖u_k-p‖}_2-{‖u_1-p‖}_2+n_{k1},$$

(5)

where $n_{k1}=n_k-n_1$ is the TDOA noise measurement. Collecting the TDOA measurements from all receiving radars into a vector, Equation (5) can be written as

$$r_1={\left[r_{21},\cdots ,r_{K1}\right]}^{\mathrm{T}}={\overline{r}}_1+n_1,$$

(6)

where ${\overline{r}}_1={[{\overline{r}}_{21},\cdots ,{\overline{r}}_{K1}]}^{\mathrm{T}}$ represents the TDOA measurement vector for the navigation position with ${\overline{r}}_{k1}=\left|\right|u_k-p{\left|\right|}_2-\left|\right|u_1-p{\left|\right|}_2$, and $n_1={[n_{21},\cdots ,n_{K1}]}^{\mathrm{T}}$ represents the TDOA noise measurement vector obeying the zero-mean Gaussian distribution with the following covariance matrix:

$$Q_{\mathrm{n}1}=E\left\{n_1{n}_1^{\mathrm{T}}\right\}.$$

(7)

The TDOA noise measurement vector $n_1$ and the position error vector $\mathsf{\Delta }u$ are mutually uncorrelated.

The objective of this paper is to perform 3D target localization with the number of receiving radars as low as possible in the presence of the GEO satellite and UAV position uncertainties. In the next section, we propose a two-stage WLS solution to tackle this problem. It should be emphasized that our localization method is not difficult to extend to the case of multiple transmitters, without needing many modifications.

2.2. Stage 1: Estimation of the Reference Range and the Target Location in an Iterative Way

The noiseless TDOA measurement ${\overline{r}}_{k1}={\overline{r}}_{Tk}-{\overline{r}}_{T1}$ and the bistatic range ${\overline{r}}_{Tk}={\overline{r}}_T+{\overline{r}}_k$ can be combined, yielding

$${\overline{r}}_{k1}={\overline{r}}_k-{\overline{r}}_1.$$

(8)

From Equation (8), we show that the influence of the GEO satellite position error on target localization can be removed by employing TDOAs. Rearranging Equation (8) as ${\overline{r}}_{k1}+{\overline{r}}_1={\overline{r}}_k$ and squaring both sides, we can get the following [21]:

$${\overline{r}}_{k1}^2+2{\overline{r}}_{k1}{\overline{r}}_1=S_k-S_1-2{\left(u_k-u_1\right)}^{\mathrm{T}}p,$$

(9)

where $S_k=u_k^{\mathrm{T}}{u}_k$. Moreover, substituting $u_k={\overline{u}}_k+\mathsf{\Delta }{u}_k$ and ${\overline{r}}_{k1}=r_{k1}-n_{k1}$ into Equation (9) and ignoring the second-order error terms [22,23], the equation error of Equation (9) is approximately expressed as

$${\varsigma }_{k1}\triangleq 2n_{k1}{r}_1^o-2{\left({\overline{u}}_k-p\right)}^{\mathrm{T}}\mathsf{\Delta }{u}_k+2{\left({\overline{u}}_1-p\right)}^{\mathrm{T}}\mathsf{\Delta }{u}_1\approx r_{k1}^2-{\overline{S}}_k+{\overline{S}}_1+2{\left({\overline{u}}_k-{\overline{u}}_1\right)}^{\mathrm{T}}p+2r_{k1}{r}_1^o,$$

(10)

where ${\overline{S}}_k={\overline{u}}_k^{\mathrm{T}}{\overline{u}}_k$. Stacking Equation (10) for all transmit/receiver pairs in a matrix form yields

$${\mathsf{\varsigma }}_1=h_1-G_1{\mathsf{\phi }}_1,$$

(11)

where ${\mathsf{\phi }}_1={[p^{\mathrm{T}},r_1^o]}^{\mathrm{T}}$, ${\mathsf{\varsigma }}_1={[{\varsigma }_{21},\cdots ,{\varsigma }_{K1}]}^{\mathrm{T}}$, and

$$h_1=\left[\begin{array}{c}{r}_{21}^2-{\overline{S}}_2+{\overline{S}}_1\\ ⋮\\ r_{K1}^2-{\overline{S}}_K+{\overline{S}}_1\end{array}\right],$$

(12)

$$G_1=-2\left[\begin{array}{cc}{\left({\overline{u}}_2-{\overline{u}}_1\right)}^{\mathrm{T}}& r_{21}\\ ⋮& ⋮\\ {\left({\overline{u}}_K-{\overline{u}}_1\right)}^{\mathrm{T}}& r_{K1}\end{array}\right].$$

(13)

Note that the number of unknown parameters in ${\mathsf{\phi }}_1$ is four, while the obtained TDOA measurement number is $(K-1)$. As a result, if we apply the canonical Chan-Ho technique to solve the problem directly, then $K\ge 5$ receiving radars are deployed to meet the requirement of the parameter identifiability. Otherwise, the larger number of unknown parameters than TDOA measurements would induce a trivial localization result. Hence, we rewrite Equation (11) as

$${\mathsf{\varsigma }}_1=h_1+2r_1{r}_1^o-{\stackrel{⌣}{\mathbf{G}}}_{1}p,$$

(14)

where

$${\stackrel{⌣}{\mathbf{G}}}_1=-2\left[\begin{array}{c}{\left({\overline{u}}_2-{\overline{u}}_1\right)}^{\mathrm{T}}\\ ⋮\\ {\left({\overline{u}}_K-{\overline{u}}_1\right)}^{\mathrm{T}}\end{array}\right].$$

(15)

Invoking the Gauss–Markov theorem [38] in Equation (14), we can obtain the WLS estimate of $p$ by minimizing ${\mathsf{\varsigma }}_1^{\mathrm{T}}{W}_1{\mathsf{\varsigma }}_1$ as

$$\hat{p}={({\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{\stackrel{⌣}{\mathbf{G}}}_1)}^{-1}{\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1\left(h_1+2r_1{r}_1^o\right),$$

(16)

where $W_1=E{\left\{{\mathsf{\varsigma }}_1{\mathsf{\varsigma }}_1^{\mathrm{T}}\right\}}^{-1}$ is the weighting matrix. Similar to the techniques in [21,26], we substitute $\hat{p}$ in Equation (16) back into the reference range $r_1^o$, yielding

$${‖{({\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{\stackrel{⌣}{\mathbf{G}}}_1)}^{-1}{\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1\left(h_1+2r_1{r}_1^o\right)-{\overline{u}}_1‖}_2\approx r_1^o.$$

(17)

After some mathematical manipulations, we obtain

$$A{r}_1^o{}^2+B{r}_1^o+C=0,$$

(18)

where

$$A=4{‖{({\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{\stackrel{⌣}{\mathbf{G}}}_1)}^{-1}{\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{r}_1‖}_2^2-1,$$

(19)

$$B=4{\left[{({\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{\stackrel{⌣}{\mathbf{G}}}_1)}^{-1}{\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{h}_1-{\overline{u}}_1\right]}^{\mathrm{T}}{({\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{\stackrel{⌣}{\mathbf{G}}}_1)}^{-1}{\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{r}_1,$$

(20)

$$C={‖{({\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{\stackrel{⌣}{\mathbf{G}}}_1)}^{-1}{\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{h}_1-{\overline{u}}_1‖}_2.$$

(21)

The weighted matrix $W_1$ can be derived as a function of the statistics of position error in the a priori information. From Equation (10), rearranging ${\mathsf{\varsigma }}_1$ yields

$${\mathsf{\varsigma }}_1=A_1{n}_1+B_1\mathsf{\Delta }u=\left[\begin{array}{cc}{A}_1& B_1\end{array}\right]\left[\begin{array}{c}{n}_1\\ \mathsf{\Delta }u\end{array}\right],$$

(22)

where

$$A_1=2\left[\begin{array}{ccc}{r}_2^o& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & r_K^o\end{array}\right],$$

(23)

$$B_1=2\left[\begin{array}{cccc}{\left({\overline{u}}_1-p\right)}^{\mathrm{T}}& -{\left({\overline{u}}_k-p\right)}^{\mathrm{T}}& \cdots & 0_{1\times 3}\\ {\left({\overline{u}}_1-p\right)}^{\mathrm{T}}& ⋮& \ddots & ⋮\\ {\left({\overline{u}}_1-p\right)}^{\mathrm{T}}& 0_{1\times 3}& \cdots & -{\left({\overline{u}}_k-p\right)}^{\mathrm{T}}\end{array}\right],$$

(24)

In this way, the weighted matrix $W_1$ is

$$W_1={\left(A_1{Q}_{\mathrm{n}1}{A}_1^{\mathrm{T}}+B_1{Q}_{\mathrm{u}}{B}_1^{\mathrm{T}}\right)}^{-1}.$$

(25)

Nevertheless, the nominal ranges and directions from the target to the receiving radars in $A_1$ and $B_1$ are unknown for the WLS solution in Equation (16). Here, we would introduce an iterative way to update $W_1$ and $\hat{p}$ alternatively. For the first iteration, we assume that the target ranges and directions are approximately equal; then, $A_1$ and $B_1$ can be simplified to

$$A_1\approx I_{K-1},$$

(26)

$$B_1\approx \left[\begin{array}{cccc}{1}_{1\times 3}& -1_{1\times 3}& \cdots & 0_{1\times 3}\\ 1_{1\times 3}& ⋮& \ddots & ⋮\\ 1_{1\times 3}& 0_{1\times 3}& \cdots & -1_{1\times 3}\end{array}\right].$$

(27)

Solving the quadratic Equation (18) to produce the estimate of reference range ${\hat{r}}_1^o$ and substituting it back into Equation (16), we obtain the first-stage target localization estimate:

$${\hat{p}}_I={({\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1{\stackrel{⌣}{\mathbf{G}}}_1)}^{-1}{\stackrel{⌣}{\mathbf{G}}}_1^{\mathrm{T}}{W}_1\left(h_1+2r_1{\hat{r}}_1^o\right).$$

(28)

Most importantly, we emphasize here that there must be only four radars participating in the first stage even though excess radars are available.

2.3. Stage 2: Refinement of the Target Location Using Taylor Series Approximations

Since the aerial target localization is implemented by a lower number of receiving radars, the target localization accuracy can be further improved by applying the first-order Taylor series expansion [24,25,39] to $r_{k1}$ at the first-stage estimate of the target location, i.e.,

$$r_{k1}\approx {‖{\overline{u}}_k-{\hat{p}}_I‖}_2-{‖{\overline{u}}_1-{\hat{p}}_I‖}_2+{\left({\mathsf{\alpha }}_k-{\mathsf{\alpha }}_1\right)}^{\mathrm{T}}\mathsf{\Delta }{p}_I+{\mathsf{\alpha }}_k^{\mathrm{T}}\mathsf{\Delta }{u}_k-{\mathsf{\alpha }}_1^{\mathrm{T}}\mathsf{\Delta }{u}_1+n_{k1},$$

(29)

where $\mathsf{\Delta }{p}_I=p-{\hat{p}}_I$ is the estimation bias of the target location in the first stage, and ${\mathsf{\alpha }}_k=({\overline{u}}_k-{\hat{p}}_I)/\left|\right|{\overline{u}}_k-{\hat{p}}_I{\left|\right|}_2$ is the unit steering vector from the target to the k-th receiving radar. Similar to Equation (10), the equation error of Equation (29) is

$${\varsigma }_{k2}=n_{k1}+{\mathsf{\alpha }}_k^{\mathrm{T}}\mathsf{\Delta }{u}_k-{\mathsf{\alpha }}_1^{\mathrm{T}}\mathsf{\Delta }{u}_1=r_{k1}-{‖{\overline{u}}_k-{\hat{p}}_I‖}_2+{‖{\overline{u}}_1-{\hat{p}}_I‖}_2-{\left({\mathsf{\alpha }}_k-{\mathsf{\alpha }}_1\right)}^{\mathrm{T}}\mathsf{\Delta }{p}_I.$$

(30)

Once again, stacking Equation (30) leads to

$${\mathsf{\varsigma }}_2=h_2-G_2\mathsf{\Delta }p,$$

(31)

where ${\mathsf{\varsigma }}_1={[{\varsigma }_{21},\cdots ,{\varsigma }_{K1}]}^{\mathrm{T}}$,

$$h_2=\left[\begin{array}{c}{r}_{21}-{‖{\overline{u}}_2-{\hat{p}}^I‖}_2+{‖{\overline{u}}_1-{\hat{p}}^I‖}_2\\ ⋮\\ r_{K1}-{‖{\overline{u}}_K-{\hat{p}}^I‖}_2+{‖{\overline{u}}_1-{\hat{p}}^I‖}_2\end{array}\right],$$

(32)

$$G_2=\left[\begin{array}{c}{\left({\mathsf{\alpha }}_2-{\mathsf{\alpha }}_1\right)}^{\mathrm{T}}\\ ⋮\\ {\left({\mathsf{\alpha }}_K-{\mathsf{\alpha }}_1\right)}^{\mathrm{T}}\end{array}\right].$$

(33)

Then, the WLS solution of $\mathsf{\Delta }{p}_I$ remaining in the first stage is

$$\mathsf{\Delta }{\hat{p}}_I={\left(G_2^{\mathrm{T}}{W}_2{G}_2\right)}^{-1}{G}_2^{\mathrm{T}}{W}_2{h}_2,$$

(34)

where $W_2=E{\left\{{\mathsf{\varsigma }}_2{\mathsf{\varsigma }}_2^{\mathrm{T}}\right\}}^{-1}$. Now, let us calculate the weighted matrix $W_2$. By rearranging ${\mathsf{\varsigma }}_2$, we get

$${\mathsf{\varsigma }}_2=A_2{n}_1+B_2\mathsf{\Delta }u=\left[\begin{array}{cc}{A}_2& B_2\end{array}\right]\left[\begin{array}{c}{n}_1\\ \mathsf{\Delta }u\end{array}\right].$$

(35)

Similarly, the iteration technique is used again, and the initial guesses of $A_2$ and $B_2$ are given by

$$A_2=A_1,$$

(36)

$$B_2=\left[\begin{array}{cccc}-{\mathsf{\alpha }}_1^{\mathrm{T}}& {\mathsf{\alpha }}_2^{\mathrm{T}}& \cdots & 0_{1\times 3}\\ -{\mathsf{\alpha }}_1^{\mathrm{T}}& ⋮& \ddots & ⋮\\ -{\mathsf{\alpha }}_1^{\mathrm{T}}& 0_{1\times 3}& \cdots & {\mathsf{\alpha }}_K^{\mathrm{T}}\end{array}\right].$$

(37)

Hence, the weighted matrix $W_2$ is

$$W_2={\left(A_2{Q}_{\mathrm{n}1}{A}_2^{\mathrm{T}}+B_2{Q}_{\mathrm{u}}{B}_2^{\mathrm{T}}\right)}^{-1}.$$

(38)

Lastly, we combine Equations (28) and (34), and the second-stage target localization becomes

$${\hat{p}}_{\mathit{II}}={\hat{p}}_I+\mathsf{\Delta }{\hat{p}}_I.$$

(39)

Furthermore, the flowchart of the proposed TDOA-based target localization method is shown in Figure 2.

2.4. Cramer–Rao Lower Bound (CRLB)

The CRLB is commonly used as a benchmark to place a lower bound of any unbiased estimator for the parameter of interest [38]. In this section, the CRLB for target localization is derived in order to evaluate the performances of different methods.

The parameter vector of interest is defined as follows:

$$\mathsf{\theta }\triangleq {\left[p^{\mathrm{T}},u^{\mathrm{T}}\right]}^{\mathrm{T}},$$

(40)

which is collected from the actual position of receiving radars. The logarithmic version of the probability density function (PDF) is constructed by utilizing TDOA measurements and a priori statistics of position errors:

$$u={\left[u_1^{\mathrm{T}},u_2^{\mathrm{T}},\cdots ,u_K^{\mathrm{T}}\right]}^{\mathrm{T}},$$

(41)

$$\ln f\left(\mathsf{\theta }\right)=-\frac{1}{2}{\left(r_1-{\overline{r}}_1\right)}^{\mathrm{H}}{Q}_{\mathrm{n}1}^{-1}\left(r_1-{\overline{r}}_1\right)-\frac{1}{2}{\left(u-\overline{u}\right)}^{\mathrm{H}}{Q}_{\mathrm{u}}^{-1}\left(u-\overline{u}\right),$$

(42)

where all UAV navigation positions are arranged into a vector as follows:

$$\overline{u}={[{\overline{u}}_1^{\mathrm{T}},{\overline{u}}_2^{\mathrm{T}},\cdots ,{\overline{u}}_K^{\mathrm{T}}]}^{\mathrm{T}}.$$

(43)

Then, the CRLB of parameter vector $\mathsf{\theta }$, calculated by the inversion of the Fisher information matrix (FIM), is given by

$$CRLB\left(\mathsf{\theta }\right)=J^{-1}\left(\mathsf{\theta }\right)=-E{\left[\frac{{\partial }^2\ln f\left(\mathsf{\theta }\right)}{\partial \mathsf{\theta }\partial {\mathsf{\theta }}^{\mathrm{T}}}\right]}^{-1}={\left[\begin{array}{cc}X& Y\\ Y^{\mathrm{T}}& Z\end{array}\right]}^{-1},$$

(44)

where

$$X=-E{\left[\frac{{\partial }^2\ln f\left(\mathsf{\theta }\right)}{\partial p\partial p^{\mathrm{T}}}\right]}^{-1}={⟮\frac{\partial {\overline{r}}_1}{\partial p}⟯}^{\mathrm{T}}{Q}_{\mathrm{n}1}^{-1}⟮\frac{\partial {\overline{r}}_1}{\partial p}⟯,$$

(45)

$$Y=-E{\left[\frac{{\partial }^2\ln f\left(\mathsf{\theta }\right)}{\partial p\partial u^{\mathrm{T}}}\right]}^{-1}={⟮\frac{\partial {\overline{r}}_1}{\partial p}⟯}^{\mathrm{T}}{Q}_{\mathrm{n}1}^{-1}⟮\frac{\partial {\overline{r}}_1}{\partial u}⟯,$$

(46)

and

$$Z=-E{\left[\frac{{\partial }^2\ln f\left(\mathsf{\theta }\right)}{\partial u\partial u^{\mathrm{T}}}\right]}^{-1}={⟮\frac{\partial {\overline{r}}_1}{\partial u}⟯}^{\mathrm{T}}{Q}_{\mathrm{n}1}^{-1}⟮\frac{\partial {\overline{r}}_1}{\partial u}⟯+Q_{\mathrm{u}}^{-1},$$

(47)

with the Jacobian transformation matrices

$$⟮\frac{\partial {\overline{r}}_1}{\partial p}⟯=\left[\frac{{\overline{u}}_2-p}{\left|\right|u_2-p{\left|\right|}_2}-\frac{{\overline{u}}_1-p}{\left|\right|u_1-p{\left|\right|}_2}\hspace{0.33em},\hspace{0.33em}\dots \hspace{0.33em},\hspace{0.33em}\hspace{0.33em}\frac{{\overline{u}}_K-p}{\left|\right|u_K-p{\left|\right|}_2}-\frac{{\overline{u}}_1-p}{\left|\right|u_1-p{\left|\right|}_2}\right],$$

(48)

and

$$⟮\frac{\partial {\overline{r}}_1}{\partial u}⟯=\left[\begin{array}{ccc}\frac{{\overline{u}}_1-p}{\left|\right|u_1-p{\left|\right|}_2}& \dots & \frac{{\overline{u}}_1-p}{\left|\right|u_1-p{\left|\right|}_2}\\ -\frac{{\overline{u}}_2-p}{\left|\right|u_2-p{\left|\right|}_2}& & \mathbf{0}\\ & \ddots & \\ \mathbf{0}& & -\frac{{\overline{u}}_K-p}{\left|\right|u_K-p{\left|\right|}_2}\end{array}\right].$$

(49)

Furthermore, we apply the matrix inverse lemma [40] to Equation (44). The CRLB for target localization is modified as follows [16]:

$$CRLB\left(p\right)=X^{-1}+X^{-1}Y{\left(Z-Y^{\mathrm{T}}{X}^{-1}Y\right)}^{-1}{Y}^{\mathrm{T}}{X}^{-1}.$$

(50)

Lastly, the minimum variance of the target location error is calculated by the matrix trace in Equation (50).

3. Results

This section provides two cases of experiments, where the minimum number of receiving radars for TDOA aerial target localization, i.e., $K=4$ is evaluated first, followed by the results for the cases of more receiving radars. The WLS-Ho method [22,23] and the WLS-Mao method [25] are introduced as the compared methods. It is assumed that the radar transmitter is equipped on a GEO satellite, with an altitude of 36,000 km. Unless otherwise specified, the signal carrier frequency is 1.25 GHz. The unequal target SNRs are measured with respect to different bistatic ranges. Here, the 3D position errors of the satellite transmitter and receiving radars are both Gaussian-distributed. For each dimension, both types of position errors are zero mean, while their variances are fixed as 100 m and 0.25 m, respectively. All results are averaged over 2000 independent Monte Carlo trials.

3.1. TDOA Target Localization with the Minimum Number of Radars

First, we consider the target localization in the absence of position errors. The UAV navigation positions are at the coordinates [100 × 10³, 100 × 10³, 23 × 10³]^T m, [200 × 10³, 300 × 10³, 22 × 10³]^T m, [−400 × 10³, 100 × 10³, 23 × 10³]^T m, and [200 × 10³, 300 × 10³, 25 × 10³]^T m, respectively. Under the minimum number of four receiving radars, the RMSEs of the target location versus the relative SNR are depicted in Figure 3. The relative SNR is defined as the SNR calculated relative to the reference receiving radar. The RMSE is computed as $RMSE=\sqrt{{\displaystyle {\sum }_{u=1}^U{‖p_u-p‖}_2^2}/U}$, where $p_u$ is the u-th target location estimation, and $U$ is the trail number. The black solid lines denote the CRLB, and the red dashed lines represent the proposed method, where the cases of the target inside and outside the radar network are marked by the plus sign and square, respectively. In Figure 3, we show that, as the relative SNR increases, the RMSE of the target location is decreased, or the localization accuracy is improved. In high-SNR regions, the proposed method can approach the theoretical CRLB. Compared to other WLS methods, the proposed method additionally considers the relationship between the target location and the reference range. Hence, by substituting the location parameters into the reference range, the estimation progresses of the target location and the reference range can be separated, which yields a lower required number of radars for the localization. After the target location is estimated in Stage 1, we can further refine the estimation results using the Taylor series approximation method in Stage 2, regardless of whether the position of the reference radar is uncertain or not. In this way, the influence of the position uncertainty weakens. Meanwhile, it is not difficult to demonstrate that a more distributed radar configuration yields better localization accuracy. Owing to the parameter unidentifiability of either the WLS-Ho method or the WLS-Mao method under four radars, their results are trivial and are not displayed here.

Then, we consider the case of the radar network with position errors. Still under four radars, Figure 4 exhibits the RMSE of the localization versus the relative SNR. We employ the plus sign and circle to represent the localization in the absence and in the presence of the position errors, respectively. Again, in both cases, our proposed method could attain the CRLB criteria as the relative SNR increases, implying that this method is robust enough against the position uncertainty. Moreover, we show that the localization performance of the former case in high-SNR regions deviates from that of the latter case, which is explained by the result of limited position accuracy. Once the position accuracy is improved, one can expect a more relatively accurate localization.

Next, we study the case of the RMSE of the localization versus UAV position error variance in Figure 5. The upward-pointing triangle, point, and downward-pointing triangle represent the cases with error variance underestimation, exactly known, and overestimation, where the three cases correspond to 10% of the error variance, unchanged error variance, and tenfold the real error variance, respectively. As the position error variance increases, we can show that the localization accuracies gradually shrink. Nevertheless, we find that the proposed method can still approximate the performance bound to a certain extent. In Figure 5, the RMSEs of these three cases perform nearly the same localization; thus, we conclude that the erroneous error variance may not play a major role in localization performance, which paves a way for the practical engineering applications.

Furthermore, to evaluate the localization accuracy for a given location of the target under a specified radar network, the GDOP metric [14,39] is invoked to address this issue. The GDOP metric is especially applied in the GPS system which is used to measure how the localization accuracy of an object in a fixed position can be attained. The GDOP is defined as

$$GDOP=\frac{1}{c}\sqrt{\frac{{\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2}{{\sigma }_{\epsilon }^2}},$$

(51)

where ${\sigma }_x^2$, ${\sigma }_y^2$, and ${\sigma }_z^2$ are the x-, y-, and z-dimensional estimation variances, respectively, and ${\sigma }_{\epsilon }^2$ is the TDOA noise measurement variance. From Equation (51), it is easy to conclude that the GDOP is purely a geometric localization evaluation tool, where the SNR factor is precluded. Under our case of the minimum number of receiving radars, the CRLB-GDOP result and WLS-proposed-GDOP result are shown in Figure 6, where the red pentagrams are the positions of the receiving radars, and the contours represent the localization accuracy that a given target can obtain. A relative SNR of 20 dB is assumed in the simulation. From Figure 6, we can see that, if the aerial target is located inside the virtual polygon by the receiving radars in a relative way, then the proposed method enables a satisfied TDOA localization accuracy. On the contrary, if the object is outside that polygon, the estimation result is badly deteriorated, especially along the line formed by any two radars. To explain this, we notice that the observation angles between the target and the radars are quite different for the former case, whereas the angles among them appear somewhat approximately equal for the latter. With more spatial diversity of the system, the method can attain better localization accuracy. These different observation angles are the foundation for the localization.

3.2. TDOA Target Localization with more Receiving Radars

Next, we consider the number of five receiving radars. The UAVs’ nominal positions are unchanged except that a new one is added at the coordinates [0 × 10³, 500 × 10³, 25 × 10³] ^T m. The RMSE of the target location versus the relative SNR is demonstrated in Figure 7. The CRLB, WLS-Ho method, WLS-Mao method, and our WLS-proposed method are represented by the black line, blue dashed line marked with a cross, green dashed line marked with a right-pointing triangle, and red dash-dotted line marked with a pentagram, respectively. We observe that the estimation accuracy of the WLS-proposed method is slightly better than that of the WLS-Ho method and WLS-Mao method. Resulting from the approximate error shown in Equation (17), the localization accuracy may be disappointing with low SNRs. Nevertheless, from the view of practical engineering, the SNRs of receiving radars would be designed to exceed 13 dB in most cases. Thus, we can conclude that the influence of this approximation is very slight.

In Figure 8, we show that, as the position error variance increases, the RMSEs of the target location are larger. For the WLS-Ho method, the reference range is the nonlinear function of the target location; hence, the deviations occur as the variance increases. The localization accuracies of the WLS-Mao method and WLS-proposed method could nearly attain that of the CRLB, whereas our method performed slightly better.

Lastly, we present the CRLB-GDOP, WLS-Ho-GDOP, WLS-Mao-GDOP, and WLS- proposed-GDOP results to close the discussion. As shown in Figure 9, the proposed method achieves better localization accuracy compared to other state-of-the-art methods. Again, more angle divergence means better localization accuracy. As for the influence of the receiving radar selection for optimal aerial target localization, it deserves further investigation.

4. Discussions

In this paper, a two-stage target localization method using TDOA measurements was developed with a minimum number of radars. Under position uncertainties, the novel TDOA-based localization method is capable of locating targets under only four radars with comparative localization accuracy in comparison to the CRLB. The localization process can be mainly divided into two stages: (1) the 3D WLS target estimates are substituted into the reference range equation to obtain a quadratic function of the reference range, and then the reference range and the target coordinates are updated in an iterative manner; (2) the target coordinates are further refined by applying the· first-order Taylor series expansion at the estimated target location in the first stage.

It should be emphasized that only four radars must be used to localize the target in the first stage even though the total radar number exceeds that number. Otherwise, the target localization would badly deteriorate. All radars involved can be considered in the second stage to improve the localization.

Compared to the WLS-Ho method, our proposal does not introduce the nonlinear approximate error between the reference range and the target location, thus alleviating the sensitivity to position errors. In contrast to the WLS-Mao method, we have more accurate target estimates in the first stage, which is vital for the second-stage refinement.

In the case of multiple transmitting radars, the proposed method can also be applied after slight modifications by taking both the measurement errors and position errors into consideration.

5. Conclusions

For most cases of two-stage target TDOA localization methods, there are five or more receiving radars to support the necessary parameter identifiability. However, the additional deployment of radars is extremely expensive or impossible in some practical scenarios. Hence, in this paper, we propose a novel two-stage WLS target localization method to alleviate exceptional employment. We first formulate the iteration between the target location and the reference range, and then the Taylor series approximation is made to refine the target location. Simulation results show that the proposed method outperformed the existing WLS techniques under four radars and was competitive with them when using more radars. The proposed method is also applicable to localization using near-space platforms and ground-based platforms. Our future work will concentrate on localization involving clock synchronization errors and phase errors to generalize our development.