3D TDOA Emitter Localization Using Conic Approximation

Abstract

This paper develops a new time difference of arrival (TDOA) emitter localization algorithm in the 3D space, employing conic approximations of hyperboloids associated with TDOA measurements. TDOA measurements are first converted to 1D angle of arrival (1D-AOA) measurements that define TDOA cones centred about axes connecting the corresponding TDOA sensor pairs. Then, the emitter location is calculated from the triangulation of 1D-AOAs, which is formulated as a system of nonlinear equations and solved by a low-complexity two-stage estimation algorithm composed of an iterative weighted least squares (IWLS) estimator and a Taylor series estimator aimed at refining the IWLS estimate. Important conclusions are reached about the optimality of sensor–emitter and sensor array geometries. The approximate efficiency of the IWLS estimator is also established under mild conditions. The new two-stage estimator is shown to be capable of outperforming the maximum likelihood estimator while performing very close to the Cramer Rao lower bound in poor sensor–emitter geometries and large noise by way of numerical simulations.

1. Introduction

Time difference of arrival (TDOA) localization is a passive emitter localization method that is capable of locating non-cooperative emitters when no time-of-transmission information is available. In the 3D space, TDOA measurements taken at pairs of sensors define hyperboloids as possible emitter locations with two foci placed at the corresponding sensors. The emitter location is determined from the intersection of multiple (at least three) hyperboloids, which implies that at least four sensors are necessary. However, because of measurement noise, TDOA hyperboloids do not intersect uniquely and the intersection point must be estimated from noisy TDOA measurements. Broadly speaking, the existing estimators for TDOA localization can be grouped into (i) maximum likelihood estimator (MLE), which takes the form of a nonlinear least squares estimator for Gaussian noise and is considered to provide the benchmark performance (see e.g., [1]), (ii) constrained “linearized” solutions based on two-stage estimators [2] and generalized trust region solution (GTRS) [3], (iii) semidefinite relaxation methods [4], and (iv) angle-of-arrival (AOA) solutions based on asymptotic approximation of TDOA hyperbolae in the 2D plane [5,6]. In this paper, we extend the 2D-AOA solutions to the 3D space using conic approximation of TDOA hyperboloids, which culminates in a new 3D TDOA localization method.

While the MLE is known to be asymptotically unbiased and efficient, it does not have a closed-form solution and requires a computationally expensive numerical search algorithm. Furthermore, the MLE cost function for TDOA localization is nonconvex [5], implying that, unless an appropriate initial guess close to the final estimate is chosen, the numerical search can diverge. Being a nonlinear estimator, the MLE also exhibits the threshold effect [7], which causes the MLE performance to degrade suddenly as the noise is increased, thereby producing unreliable estimates at large noise levels. These observations have motivated the development of several alternative algorithms for TDOA localization over the last four decades, starting with the seminal work on hyperbolic location in [2]. More recently, GTRS was exploited to solve the TDOA localization problem using a quadratic cost function with quadratic constraints [3]. While the GTRS solution achieves a localization performance on par with the MLE, it has a large computational complexity [8]. In [9], a least squares estimator was presented for TDOA localization in the 3D space, accounting for the constant bias in TDOA measurements. The work in [10] presents a weighted least squares (WLS) estimator with the cone tangent plane constraint for 2D TDOA localization. The Lagrange programming neural network was applied to the TDOA localization problem in [11] using an analogue neural network model. A hybrid firefly algorithm was proposed in [12], which combines the WLS estimator with the firefly algorithm to restrict the search region, achieving reduced computational complexity and improved accuracy. In [13] closed-form 2D and 3D TDOA localization algorithms were developed in modified polar representation by minimizing a quadratic cost function with a quadratic constraint. Successive unconstrained minimization and GTRS were employed to solve the constrained optimization problem. An algorithm that solves the 3D TDOA problem uniquely in a sensor network with four sensors rather than a minimum of five sensors was proposed in [14], exploiting confidence regions. Semidefinite programming based TDOA localization algorithms were developed to solve the maximum likelihood estimation problem for emitter location, as well as joint emitter location and propagation speed estimation in the presence of sensor position errors (see, e.g., [4,15]). In [16], source localization from a set of squared noisy range difference measurements was considered. The localization problem was solved in the least squares sense by expressing the source location in polar/spherical coordinates, which leads to a quotient of two quadratic forms, whose constrained maximization yields an easy solution.

In this paper, we propose a novel 3D TDOA localization algorithm based on an approximation of TDOA hyperboloids by cones. This is similar to the approximation of TDOA hyperbolae in the 2D plane by asymptotes, which results in a bearings-only localization problem [6]. However, the localization problem in the 3D space is significantly more challenging than in the 2D plane, as unique azimuth and elevation angles for the emitter are not available from TDOA cones. To overcome this, TDOA measurements are first converted to 1D-AOA measurements [17] that define TDOA cones centred about axes connecting the corresponding TDOA sensor pairs. Then, the emitter location is calculated from the triangulation of 1D-AOAs, which is formulated as a nonlinear matrix equation and solved using a low-complexity two-stage algorithm composed of an iterative weighted least squares (IWLS) estimator that improves emitter range estimates, and a Taylor series estimator aimed at refining the final IWLS estimate. The optimality of sensor–emitter and sensor array geometries is considered, and important conclusions are reached as to what determines good geometry in terms of range differences of arrival, orientation of intersensor vectors and emitter ranges from the sensors. The approximate efficiency of the IWLS estimator is also established under mild conditions. The performance improvement of the proposed 3D TDOA localization algorithm over the MLE, which provides the benchmark performance, is demonstrated by way of numerical simulations.

The paper is organized as follows. Section 2 describes the 3D TDOA localization problem and presents the MLE and the Cramer Rao lower bound. In Section 3, conic approximations for TDOA hyperboloids are obtained in terms of 1D-AOAs and the covariance matrix of the 1D-AOA noise is derived. Section 4 converts the 3D TDOA localization problem into the triangulation of 1D-AOAs and formulates a nonlinear matrix equation in the unknown emitter location. Section 5 presents the new two-stage 3D TDOA localization algorithm to solve the system of nonlinear equations for the emitter location and analyzes the efficiency of the IWLS estimator in its first stage. Comparative simulation examples are presented to demonstrate the superior performance of the new algorithm against the MLE in Section 6. The concluding remarks are presented in Section 7.

2. TDOA Localization in 3D Space

In 3D TDOA localization, the objective is to estimate the location of an emitter at $\mathit{s}={[x,y,z]}^T$ (where ${}^T$ denotes the matrix transpose) using TDOA measurements obtained from N sensors ($N\ge 4$) positioned at ${\mathit{r}}_i={[x_i,y_i,z_i]}^T$, $i=1,\dots ,N$ (see Figure 1). The TDOA measurements at sensors i and j are given by

$${\tau }_{ij}={\tau }_j-{\tau }_i,i,j\in \{1,\dots ,N\},i\ne j$$

(1)

where ${\tau }_i$ is the time it takes for the signal transmitted from the emitter to arrive at sensor i:

$${\tau }_i=\frac{\parallel {\mathit{d}}_i\parallel }{c}.$$

(2)

Here, $\parallel ·\parallel$ denotes the Euclidean norm, ${\mathit{d}}_i$ is the emitter range vector from the sensor at ${\mathit{r}}_i$:

$${\mathit{d}}_i=\mathit{s}-{\mathit{r}}_i$$

(3)

and c is the speed of propagation (speed of light in free space).

Using (1) and (2), the range difference of arrival (RDOA), $g_{ij}$, becomes

$$\begin{array}{cc}\hfill g_{ij}& =\parallel {\mathit{d}}_j\parallel -\parallel {\mathit{d}}_i\parallel \hfill \end{array}$$

(4)

$$\begin{array}{c}=c{\tau }_{ij}.\hfill \end{array}$$

(5)

Noting that TDOA and RDOA only differ by a scaling factor c, we will use TDOA and RDOA interchangeably.

Each RDOA defines a hyperboloid of possible emitter locations, as depicted in Figure 1. It is common practice to nominate one of the sensors as the reference receiver and take all TDOA measurements with respect to it [2]. We assume that the sensor at ${\mathit{r}}_1$ is the reference sensor. Given the RDOAs with respect to the reference sensor, $g_{1i}$, $i=2,\dots ,N$, the emitter location $\mathit{s}$ is obtained from the intersection of $N-1$ hyperboloids:

$$\begin{array}{cc}\hfill \parallel \mathit{s}-{\mathit{r}}_2\parallel -\parallel \mathit{s}-{\mathit{r}}_1\parallel & =g_{12}\hfill \\ \hfill \parallel \mathit{s}-{\mathit{r}}_3\parallel -\parallel \mathit{s}-{\mathit{r}}_1\parallel & =g_{13}\hfill \\ \hfill ⋮\\ \hfill \parallel \mathit{s}-{\mathit{r}}_N\parallel -\parallel \mathit{s}-{\mathit{r}}_1\parallel & =g_{1N}.\hfill \end{array}$$

(6)

To solve the above set of nonlinear equations for $\mathit{s}$, a minimum of three equations are required (i.e., $N\ge 4$) since there are three unknowns, viz., the x, y and z coordinates of the emitter location $\mathit{s}$. However, in practice, more than four sensors may be necessary and desirable to ensure a unique solution and to improve the accuracy of the solution.

True RDOAs are unknown and only their estimates obtained from noisy received signals are available. RDOAs can be estimated using the method of generalized cross-correlation [18]. Noisy RDOA measurements are modelled as

$${\tilde{g}}_{1i}=g_{1i}+n_{1i}$$

(7)

where the RDOA noise $n_{1i}$ is assumed to be zero-mean Gaussian. For i.i.d. additive Gaussian noise at each sensor, the covariance of RDOA noise $n_{1i}$ is

$$\begin{array}{cc}& \mathsf{\sum }=E\left\{\left[\begin{array}{c}{n}_{12}\\ n_{13}\\ ⋮\\ n_{1N}\end{array}\right]\left[\begin{array}{ccc}{n}_{12}& n_{13}\cdots & n_{1N}\end{array}\right]\right\}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & =\frac{{\sigma }_n^2}{2}({\mathit{I}}_{N-1}+{\mathbf{1}}_{N-1})\hfill \end{array}$$

(9)

where ${\sigma }_n^2$ is the RDOA noise variance, ${\mathit{I}}_N$ is the $N\times N$ identity matrix and ${\mathbf{1}}_N$ is the $N\times N$ matrix of ones. Note that $\mathsf{\sum }$ is not a diagonal matrix, which means that the $n_{1i}$ are correlated.

The maximum likelihood estimator (MLE) for the emitter location is constructed by maximizing the joint probability density function of the noisy RDOA measurements conditioned on the emitter location. The MLE takes the form of a nonlinear weighted least squares estimator for Gaussian noise:

$${\widehat{\mathit{s}}}_{\mathrm{MLE}}=\underset{\mathit{s}}{\arg \min }{\mathit{e}}^T\left(\mathit{s}\right){\mathsf{\sum }}^{-1}\mathit{e}\left(\mathit{s}\right)$$

(10)

where

$$\mathit{e}\left(\mathit{s}\right)=\left[\begin{array}{c}{\tilde{g}}_{12}\\ {\tilde{g}}_{13}\\ ⋮\\ {\tilde{g}}_{1N}\end{array}\right]-\left[\begin{array}{c}\parallel \mathit{s}-{\mathit{r}}_2\parallel -\parallel \mathit{s}-{\mathit{r}}_1\parallel \\ \parallel \mathit{s}-{\mathit{r}}_3\parallel -\parallel \mathit{s}-{\mathit{r}}_1\parallel \\ ⋮\\ \parallel \mathit{s}-{\mathit{r}}_N\parallel -\parallel \mathit{s}-{\mathit{r}}_1\parallel \end{array}\right].$$

(11)

The MLE has the desirable properties of being asymptotically unbiased and efficient (i.e., its covariance becomes identical to the Cramer–Rao lower bound as N tends to infinity). However, for a finite number of sensors, it is only approximately unbiased and efficient. Being a nonlinear estimator, it is subject to the threshold effect (sudden degradation in estimation performance) as the measurement noise variance is increased.

As the MLE does not have a closed-form solution, it requires a numerical search algorithm. For this purpose, several methods have been used, such as the Gauss–Newton (GN) algorithm, the Levenberg–Marquardt algorithm and the Nelder–Mead simplex method (see, e.g., [19,20,21,22]).

For Gaussian measurement noise, the Cramer–Rao lower bound (CRLB) for TDOA localization is given by

$$\mathrm{CRLB}={\left({\mathit{J}}_{\mathit{o}}^T{\mathsf{\sum }}^{-1}{\mathit{J}}_{\mathit{o}}\right)}^{-1}$$

(12)

where ${\mathit{J}}_{\mathit{o}}$ is the Jacobian of $\mathit{e}\left(\mathit{s}\right)$ computed at the true emitter location $\mathit{s}$:

$${\mathit{J}}_{\mathit{o}}=\left[\begin{array}{c}{\mathit{v}}_1^T\left(\mathit{s}\right)-{\mathit{v}}_2^T\left(\mathit{s}\right)\\ {\mathit{v}}_1^T\left(\mathit{s}\right)-{\mathit{v}}_3^T\left(\mathit{s}\right)\\ ⋮\\ {\mathit{v}}_1^T\left(\mathit{s}\right)-{\mathit{v}}_N^T\left(\mathit{s}\right)\end{array}\right],{\mathit{v}}_i^T\left(\mathit{s}\right)=\frac{\mathit{s}-{\mathit{r}}_i}{\parallel \mathit{s}-{\mathit{r}}_i\parallel }.$$

(13)

3. Conic Approximation of TDOA Hyperboloids

In many applications of TDOA localization, the TDOA measurements are processed to obtain directional information about the signal source. However, in the 3D space, TDOAs correspond to 1D-AOAs [17] that define cones centred about the axes connecting the TDOA sensor pairs. Therefore, they do not contain precise three-dimensional direction information. The cones approximate TDOA hyperboloids with increased accuracy in the far field (as the emitter range increases). For a sufficiently large emitter range from the sensors (e.g., the emitter range is an order of magnitude larger than the maximum separation between the sensors), the TDOA hyperboloids can be substituted by TDOA cones with negligible error, which are defined by

• The mid-point between TDOA sensor pairs

  $${\mathit{m}}_{1i}=\frac{1}{2}({\mathit{r}}_1+{\mathit{r}}_i),i=2,3,\dots ,N$$

  (14)
• The unit vector for the TDOA sensor pair direction

  $${\mathit{u}}_{1i}=\frac{{\mathit{r}}_{1i}}{\parallel {\mathit{r}}_{1i}\parallel },{\mathit{r}}_{1i}={\mathit{r}}_i-{\mathit{r}}_1$$

  (15)
• The 1D-AOA

  $${\theta }_{1i}={cos}^{-1}\left(-\frac{g_{1i}}{\parallel {\mathit{r}}_{1i}\parallel }\right),0\le {\theta }_{1i}\le \pi .$$

  (16)

Geometrically, the 1D-AOA, ${\theta }_{1i}$, is the angle that the TDOA cone makes with its axis passing through the sensor pair ${\mathit{r}}_1$ and ${\mathit{r}}_i$, as depicted in Figure 2.

In 2D TDOA localization, the 1D-AOA corresponds to a bearing angle and has sign ambiguity, as ${\theta }_{1i}$ and $-{\theta }_{1i}$ both give the same TDOA in the far field. This ambiguity needs to be resolved, e.g., by resorting to the clustering of hyperbolic asymptotes as described in [6]. However, in 3D TDOA localization, no such ambiguity exists, as the localization problem boils down to the intersection of cones rather than hyperbolic asymptotes, and the cones are not influenced by the sign ambiguity in 1D-AOAs.

Substituting the noisy RDOA measurements for the true RDOAs in (16) and using a truncated Taylor series expansion to retain the linear terms, we obtain

$${\tilde{\theta }}_{1i}={\theta }_{1i}+{ϵ}_{1i},i=2,3,\dots ,N$$

(17)

where ${\tilde{\theta }}_{1i}$ denotes the noisy 1D-AOA measurement

$${\tilde{\theta }}_{1i}={cos}^{-1}\left(-\frac{{\tilde{g}}_{1i}}{\parallel {\mathit{r}}_{1i}\parallel }\right),0\le {\tilde{\theta }}_{1i}\le \pi$$

(18)

and ${ϵ}_{1i}$ is the 1D-AOA noise approximated by

$$\begin{array}{cc}\hfill {ϵ}_{1i}& \approx \frac{n_{1i}}{\sqrt{\parallel {\mathit{r}}_{1i}{\parallel }^2-g_{1i}^2}}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & \approx \frac{n_{1i}}{\sqrt{\parallel {\mathit{r}}_{1i}{\parallel }^2-{\parallel {\mathit{r}}_{1i}\parallel }^2{cos}^2{\theta }_{1i}}}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & \approx \frac{n_{1i}}{\parallel {\mathit{r}}_{1i}\parallel sin{\theta }_{1i}}\hfill \end{array}$$

(21)

which is zero-mean Gaussian. We have used (16) in (20) to arrive at the final expression (21).

We remark that, based on (19), the 1D-AOA noise is minimized if $g_{1i}=0$ with ${\theta }_{1i}=\pi /2$ rad; i.e., the emitter is equidistant from the sensor pair ${\mathit{r}}_1$ and ${\mathit{r}}_i$. This suggests that a necessary condition for optimal sensor–emitter geometry would be to have all TDOA sensors equidistant from the emitter, assuming of course that prior knowledge of the emitter location is available. Referring to (21) we see that ${ϵ}_{1i}$ tends to infinity if ${\theta }_{1i}=0$ or ${\theta }_{1i}=\pi$, which happens if the emitter and TDOA sensor pair ${\mathit{r}}_1$ and ${\mathit{r}}_i$ are collinear. This represents the worst TDOA geometry and should be avoided. If $|{\tilde{g}}_{1i}|/\parallel {\mathit{r}}_{1i}\parallel >1$ due to large measurement noise and/or the TDOA sensor pair being almost collinear with the emitter, ${\tilde{\theta }}_{1i}$ cannot be computed from (18) as the arccosine function will have an argument with an absolute value exceeding one. In such cases, dropping those TDOA measurements with $|{\tilde{g}}_{1i}|/\parallel {\mathit{r}}_{1i}\parallel >1$ from the estimation process (if $N>4$) or using a different reference sensor may resolve the problem. To sum up, an optimal sensor–emitter geometry would have $0\le |{\tilde{g}}_{1i}|/\parallel {\mathit{r}}_{1i}\parallel \ll 1$, $i=2,3,\dots ,N$.

The covariance of the 1D-AOA noise ${ϵ}_{1i}$ is

$$\begin{array}{c}\mathit{K}=E\left\{\left[\begin{array}{c}{ϵ}_{12}\\ ⋮\\ {ϵ}_{1N}\end{array}\right]\left[\begin{array}{ccc}{ϵ}_{12}& \cdots & {ϵ}_{1N}\end{array}\right]\right\}\hfill \end{array}$$

(22)

$$\begin{array}{c}=\mathit{D}\mathsf{\sum }\mathit{D}\hfill \end{array}$$

(23)

where $\mathsf{\sum }$ is the $(N-1)\times (N-1)$ covariance matrix of TDOA noise defined in (9) and $\mathit{D}$ is a diagonal matrix of size $N-1$ given by

$$\mathit{D}=\mathrm{diag}\left(\frac{1}{\parallel {\mathit{r}}_{12}\parallel sin{\theta }_{12}},\frac{1}{\parallel {\mathit{r}}_{13}\parallel sin{\theta }_{13}},\cdots ,\frac{1}{\parallel {\mathit{r}}_{1N}\parallel sin{\theta }_{1N}}\right).$$

(24)

4. Triangulation of 1D-AOAs

The 3D TDOA localization problem is converted to an equivalent AOA localization problem using the 1D-AOAs ${\theta }_{1i}$ obtained from the conic approximation discussed in the previous section, as depicted in Figure 3. The $N-1$ TDOA sensor pair midpoints ${\mathit{m}}_{1i}$ act as virtual sensors. The 1D-AOA based localization problem considered here is different from the conventional 3D AOA localization problem. While in 3D localization, azimuth and elevation angles define directional vectors emanating from the AOA sensors, the equivalent AOA localization problem in Figure 3 essentially employs cones defined by the 1D-AOAs. This difference rules out a straightforward application of 3D AOA localization algorithms to the 3D TDOA problem. In 2D localization, however, the 1D-AOAs reduce to bearing angles in the 2D plane, which allows the application of existing 2D AOA localization algorithms after simple modification (see, e.g., [6]).

Referring to Figure 3, we note that the emitter location $\mathit{s}$ is related to the 1D-AOAs via

$${\mathit{u}}_{1i}^T\mathit{s}=d_{1i}\left(\mathit{s}\right)cos{\theta }_{1i}+{\mathit{u}}_{1i}^T{\mathit{m}}_{1i},i=2,3,\dots ,N$$

(25)

where $d_{1i}\left(\mathit{s}\right)=\parallel {\mathit{d}}_{1i}\left(\mathit{s}\right)\parallel$ is the emitter range from ${\mathit{m}}_{1i}$ with ${\mathit{d}}_{1i}\left(\mathit{s}\right)=\mathit{s}-{\mathit{m}}_{1i}$ denoting the emitter range vector originating from ${\mathit{m}}_{1i}$, and ${\mathit{u}}_{1i}$ is the unit vector pointing from the reference sensor ${\mathit{r}}_1$ to ${\mathit{r}}_i$ along the TDOA cone axis:

$${\mathit{u}}_{1i}=\frac{{\mathit{r}}_i-{\mathit{r}}_1}{\parallel {\mathit{r}}_i-{\mathit{r}}_1\parallel }.$$

(26)

Substituting ${\theta }_{1i}={\tilde{\theta }}_{1i}-{ϵ}_{1i}$ into (25), we have

$$\begin{array}{cc}{\mathit{u}}_{1i}^T\mathit{s}\hfill & =d_{1i}\left(\mathit{s}\right)\left(cos{\tilde{\theta }}_{1i}cos{ϵ}_{1i}+sin{\tilde{\theta }}_{1i}sin{ϵ}_{1i}\right)+{\mathit{u}}_{1i}^T{\mathit{m}}_{1i}\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & \approx d_{1i}\left(\mathit{s}\right)\left(cos{\tilde{\theta }}_{1i}+{ϵ}_{1i}sin{\tilde{\theta }}_{1i}\right)+{\mathit{u}}_{1i}^T{\mathit{m}}_{1i}\hfill \end{array}$$

(28)

where we assume ${ϵ}_{1i}\approx 0$ so that $cos{ϵ}_{1i}\approx 1$ and $sin{ϵ}_{1i}\approx {ϵ}_{1i}$. Thus, using the noisy 1D-AOAs, (25) is replaced by

$${\mathit{u}}_{1i}^T\mathit{s}=d_{1i}\left(\mathit{s}\right)cos{\tilde{\theta }}_{1i}+{\mathit{u}}_{1i}^T{\mathit{m}}_{1i}+{\eta }_{1i}\left(\mathit{s}\right)$$

(29)

where the noise term ${\eta }_{1i}\left(\mathit{s}\right)$ is given by

$${\eta }_{1i}\left(\mathit{s}\right)\approx {ϵ}_{1i}{d}_{1i}\left(\mathit{s}\right)sin{\theta }_{1i}.$$

(30)

Substituting (21) into the above equation yields

$${\eta }_{1i}\left(\mathit{s}\right)\approx \frac{d_{1i}\left(\mathit{s}\right)}{\parallel {\mathit{r}}_{1i}\parallel }{n}_{1i}.$$

(31)

Stacking (29) for $i=2,3,\dots ,N$, we obtain the following matrix equation which is nonlinear in $\mathit{s}$:

$$\underset{\mathit{A}}{\underbrace{\left[\begin{array}{c}{\mathit{u}}_{12}^T\\ {\mathit{u}}_{13}^T\\ ⋮\\ {\mathit{u}}_{1N}^T\end{array}\right]}}=\underset{\mathit{f}\left(\mathit{s}\right)}{\underbrace{\left[\begin{array}{c}{d}_{12}\left(\mathit{s}\right)cos{\tilde{\theta }}_{12}\\ d_{13}\left(\mathit{s}\right)cos{\tilde{\theta }}_{13}\\ ⋮\\ d_{1N}\left(\mathit{s}\right)cos{\tilde{\theta }}_{1N}\end{array}\right]}}+\underset{\mathit{b}}{\underbrace{\left[\begin{array}{c}{\mathit{u}}_{12}^T{\mathit{m}}_{12}\\ {\mathit{u}}_{13}^T{\mathit{m}}_{13}\\ ⋮\\ {\mathit{u}}_{1N}^T{\mathit{m}}_{1N}\end{array}\right]}}+\underset{\mathit{\eta }\left(\mathit{s}\right)}{\underbrace{\left[\begin{array}{c}{\eta }_{12}\left(\mathit{s}\right)\\ {\eta }_{13}\left(\mathit{s}\right)\\ ⋮\\ {\eta }_{1N}\left(\mathit{s}\right)\end{array}\right]}}.$$

(32)

Finding a closed-form estimate for $\mathit{s}$ in $\mathit{A}\mathit{s}\approx \mathit{f}\left(\mathit{s}\right)+\mathit{b}$, akin to least squares estimation, is not possible because of the nonlinear dependence of the vector $\mathit{f}\left(\mathit{s}\right)$ on the unknown $\mathit{s}$. A nonlinear least squares solution may be attempted using an iterative numerical search algorithm. However, this would not be attractive from the computational complexity point of view, especially when compared with the MLE.

The covariance of the noise vector $\mathit{\eta }\left(\mathit{s}\right)$ is

$$\begin{array}{cc}& \mathit{C}=E\left\{\mathit{\eta }\left(\mathit{s}\right){\mathit{\eta }}^T\left(\mathit{s}\right)\right\}\hfill \end{array}$$

(33)

$$\begin{array}{c}=\mathit{L}\mathsf{\sum }\mathit{L}\hfill \end{array}$$

(34)

where $\mathit{L}$ is the diagonal matrix:

$$\mathit{L}=\mathrm{diag}\left(\frac{d_{12}\left(\mathit{s}\right)}{\parallel {\mathit{r}}_{12}\parallel },\frac{d_{13}\left(\mathit{s}\right)}{\parallel {\mathit{r}}_{13}\parallel },\cdots ,\frac{d_{1N}\left(\mathit{s}\right)}{\parallel {\mathit{r}}_{1N}\parallel }\right).$$

(35)

Note that in (34), $\mathsf{\sum }$ can be replaced by $(\mathit{I}+\mathbf{1})$ by dropping the proportionality factor ${\sigma }_n^2/2$ with no performance penalty. In other words, prior knowledge of RDOA noise variance ${\sigma }_n^2$ is not necessary.

5. New Algorithm for 3D TDOA Localization Using 1D-AOAs

5.1. Algorithm Derivation

We propose a two-stage algorithm to solve (32) for $\mathit{s}$. In the first stage, (32) is approximated as a linear matrix equation with two unknowns $\mathit{s}$ and d by assuming that $d=d_{12}\left(\mathit{s}\right)=d_{13}\left(\mathit{s}\right)=\cdots =d_{1N}\left(\mathit{s}\right)$ (i.e., equal emitter range from all sensor pair midpoints ${\mathit{m}}_{1i}$), which yields

$$\begin{array}{c}\mathit{As}=d\underset{\mathit{c}}{\underbrace{\left[\begin{array}{c}cos{\tilde{\theta }}_{12}\\ cos{\tilde{\theta }}_{13}\\ ⋮\\ cos{\tilde{\theta }}_{1N}\end{array}\right]}}+\mathit{b}+\mathit{\eta }\left(\mathit{s}\right)\end{array}$$

(36)

$$[\mathit{A}-\mathit{c}]\left[\begin{array}{c}\mathit{s}\\ d\end{array}\right]\approx \mathit{b}.$$

(37)

Since the unknown vector now includes a nuisance parameter d, (37) becomes overdetermined with a unique solution if there are five or more sensors; therefore, we assume $N\ge 5$. In addition, the unit vectors ${\mathit{u}}_{1i}$ that form the matrix $\mathit{A}$ must be linearly independent to avoid rank deficiency in $\mathit{A}$. This necessarily rules out all sensors forming a linear array or those that are on the same 2D horizontal plane. The estimation accuracy will improve if the condition number of $\mathit{A}$ (i.e., the ratio of its largest singular value to the smallest) is small and as close to one as possible. This is also achieved by selecting sensor locations that produce linearly independent ${\mathit{u}}_{1i}$. In summary, an optimal sensor array geometry would minimize the condition number of $\mathit{A}$.

Equation (37) is easily solved in the weighted least squares sense using

$$\left[\begin{array}{c}{\widehat{\mathit{s}}}_0\\ \widehat{d}\end{array}\right]={\left({[\mathit{A}-\mathit{c}]}^T{\mathit{W}}_0[\mathit{A}-\mathit{c}]\right)}^{-1}{[\mathit{A}-\mathit{c}]}^T{\mathit{W}}_0\mathit{b}$$

(38)

where the weighting matrix is

$$\begin{array}{c}{\mathit{W}}_0={\mathit{C}}_0^{-1},{\mathit{C}}_0={\mathit{L}}_0\mathsf{\sum }{\mathit{L}}_0\end{array}$$

(39)

$$\begin{array}{c}{\mathit{L}}_0=\mathrm{diag}\left(\frac{1}{\parallel {\mathit{r}}_{12}\parallel },\frac{1}{\parallel {\mathit{r}}_{13}\parallel },\cdots ,\frac{1}{\parallel {\mathit{r}}_{1N}\parallel }\right)\end{array}$$

(40)

which is obtained from the covariance matrix $\mathit{C}$ in (34) by dropping the constant proportionality factor d. The assumption of $d_{12}\left(\mathit{s}\right)=d_{13}\left(\mathit{s}\right)=\cdots =d_{1N}\left(\mathit{s}\right)$ is a particularly good approximation for large emitter range to sensor baseline ratio situations characterized by a tight clustering of the sensors away from the emitter.

Next, we refine the weighted least squares estimate of the emitter location in (38) by re-estimating its weighting matrix. Using the initial estimate ${\widehat{\mathit{s}}}_0$ in (38), the new weighting matrix becomes

$${\mathit{W}}_1={\mathit{C}}_1^{-1},{\mathit{C}}_1={\mathit{L}}_1\mathsf{\sum }{\mathit{L}}_1$$

(41)

where

$${\mathit{L}}_1=\left(\frac{d_{12}\left({\widehat{\mathit{s}}}_0\right)}{\parallel {\mathit{r}}_{12}\parallel },\frac{d_{13}\left({\widehat{\mathit{s}}}_0\right)}{\parallel {\mathit{r}}_{13}\parallel },\cdots ,\frac{d_{1N}\left({\widehat{\mathit{s}}}_0\right)}{\parallel {\mathit{r}}_{1N}\parallel }\right)$$

(42)

with

$$d_{1i}\left({\widehat{\mathit{s}}}_0\right)=\parallel {\widehat{\mathit{s}}}_0-{\mathit{m}}_{1i}\parallel .$$

(43)

Replacing (37) with

$$\mathit{As}\approx \mathit{f}\left({\widehat{s}}_0\right)+\mathit{b}$$

(44)

in accordance with (32), the new weighted least squares estimate for the emitter location is

$${\widehat{\mathit{s}}}_1={\left({\mathit{A}}^T{\mathit{W}}_1\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{W}}_1\left(\mathit{f}\left({\widehat{s}}_0\right)+\mathit{b}\right).$$

(45)

Starting with ${\widehat{\mathit{s}}}_0$ in (38), (45) is computed iteratively using the following equations in each iteration for k$=1,2,\dots$

$$\begin{array}{c}{\mathit{W}}_k={\mathit{C}}_k^{-1},{\mathit{C}}_k={\mathit{L}}_k\mathsf{\sum }{\mathit{L}}_k\end{array}$$

(46)

$$\begin{array}{c}{\mathit{L}}_k=\mathrm{diag}\left(\frac{d_{12}\left({\widehat{\mathit{s}}}_{k-1}\right)}{\parallel {\mathit{r}}_{12}\parallel },\frac{d_{13}\left({\widehat{\mathit{s}}}_{k-1}\right)}{\parallel {\mathit{r}}_{13}\parallel },\cdots ,\frac{d_{1N}\left({\widehat{\mathit{s}}}_{k-1}\right)}{\parallel {\mathit{r}}_{1N}\parallel }\right)\end{array}$$

(47)

$${\widehat{\mathit{s}}}_k={\left({\mathit{A}}^T{\mathit{W}}_k\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{W}}_k\left(\mathit{f}\left({\widehat{s}}_{k-1}\right)+\mathit{b}\right)$$

(48)

which produces an iterative weighted least squares (IWLS) estimator. The iterations are stopped when the difference between successive estimates is sufficiently small; i.e., $\parallel {\widehat{\mathit{s}}}_k-{\widehat{\mathit{s}}}_{k-1}\parallel <\gamma$ where $\gamma$ is a threshold, or the maximum number of iterations is reached.

The second stage uses the final IWLS estimate, denoted ${\widehat{\mathit{s}}}_{\mathrm{IWLS}}$, calculated from the iterations (46)–(48) as an initial estimate for the Taylor series estimator [1] based on the MLE:

$$\widehat{\mathit{s}}={\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\left({\mathit{J}}^T{\mathsf{\sum }}^{-1}\mathit{J}\right)}^{-1}{\mathit{J}}^T{\mathsf{\sum }}^{-1}\mathit{e}$$

(49)

where $\mathit{J}$ is the Jacobian matrix, defined as

$$\mathit{J}=\left[\begin{array}{c}{\mathit{v}}_1^T\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)-{\mathit{v}}_2^T\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)\\ {\mathit{v}}_1^T\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)-{\mathit{v}}_3^T\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)\\ ⋮\\ {\mathit{v}}_1^T\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)-{\mathit{v}}_N^T\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)\end{array}\right],{\mathit{v}}_i^T\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)=\frac{{\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\mathit{r}}_i}{\parallel {\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\mathit{r}}_i\parallel }$$

(50)

and

$$\mathit{e}=\left[\begin{array}{c}{\tilde{g}}_{12}\\ {\tilde{g}}_{13}\\ ⋮\\ {\tilde{g}}_{1N}\end{array}\right]-\left[\begin{array}{c}\parallel {\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\mathit{r}}_2\parallel -\parallel {\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\mathit{r}}_1\parallel \\ \parallel {\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\mathit{r}}_3\parallel -\parallel {\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\mathit{r}}_1\parallel \\ ⋮\\ \parallel {\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\mathit{r}}_N\parallel -\parallel {\widehat{\mathit{s}}}_{\mathrm{IWLS}}-{\mathit{r}}_1\parallel \end{array}\right].$$

(51)

5.2. Performance Analysis of the IWLS Estimator

In this subsection, we establish that the IWLS estimator is approximately efficient under the following assumptions:

Assumption 1. 

The TDOA noise is sufficiently small.

Assumption 2. 

All sensors are at approximately the same distance from the emitter.

Assumption 3. 

Intersensor distances are small compared with the emitter range.

To begin with, consider the final IWLS estimate

$${\widehat{\mathit{s}}}_{\mathrm{IWLS}}={\left({\mathit{A}}^T{\mathit{W}}_{\mathrm{IWLS}}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{W}}_{\mathrm{IWLS}}\left(\mathit{f}\left({\widehat{s}}_{\mathrm{IWLS}}\right)+\mathit{b}\right)$$

(52)

where, under Assumption 1, the weighting matrix is

$$\begin{array}{cc}\hfill {\mathit{W}}_{\mathrm{IWLS}}& =\mathrm{diag}\left(\frac{\parallel {\mathit{r}}_{12}\parallel }{d_{12}\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)},\frac{\parallel {\mathit{r}}_{13}\parallel }{d_{13}\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)},\cdots ,\frac{\parallel {\mathit{r}}_{1N}\parallel }{d_{1N}\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)}\right){\mathsf{\sum }}^{-1}\hfill \\ \hfill & \times \mathrm{diag}\left(\frac{\parallel {\mathit{r}}_{12}\parallel }{d_{12}\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)},\frac{\parallel {\mathit{r}}_{13}\parallel }{d_{13}\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)},\cdots ,\frac{\parallel {\mathit{r}}_{1N}\parallel }{d_{1N}\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}\right)}\right)\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill & \approx \mathrm{diag}\left(\frac{\parallel {\mathit{r}}_{12}\parallel }{d_{12}\left(\mathit{s}\right)},\frac{\parallel {\mathit{r}}_{13}\parallel }{d_{13}\left(\mathit{s}\right)},\cdots ,\frac{\parallel {\mathit{r}}_{1N}\parallel }{d_{1\mathit{s}N}\left(\mathit{s}\right)}\right){\mathsf{\sum }}^{-1}\mathrm{diag}\left(\frac{\parallel {\mathit{r}}_{12}\parallel }{d_{12}\left(\mathit{s}\right)},\frac{\parallel {\mathit{r}}_{13}\parallel }{d_{13}\left(\mathit{s}\right)},\cdots ,\frac{\parallel {\mathit{r}}_{1N}\parallel }{d_{1N}\left(\mathit{s}\right)}\right)\hfill \end{array}$$

(54)

$$\begin{array}{c}\approx {\mathit{C}}^{-1}.\end{array}$$

(55)

The error covariance of the IWLS estimator is given by

$${\mathit{C}}_{\mathrm{IWLS}}=E\left\{\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}-\mathit{s}\right){\left({\widehat{\mathit{s}}}_{\mathrm{IWLS}}-\mathit{s}\right)}^T\right\}$$

(56)

where, using ${\mathit{W}}_{\mathrm{IWLS}}\approx {\mathit{C}}^{-1}$ in (55), we have

$$\begin{array}{cc}\hfill {\widehat{\mathit{s}}}_{\mathrm{IWLS}}-\mathit{s}& ={\left({\mathit{A}}^T{\mathit{W}}_{\mathrm{IWLS}}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{W}}_{\mathrm{IWLS}}\left(\mathit{f}\left({\widehat{s}}_{\mathrm{IWLS}}\right)+\mathit{b}\right)-{\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{C}}^{-1}\hfill \\ \hfill & \times \left(\mathit{f}\left(\mathit{s}\right)+\mathit{b}+\mathit{\eta }\left(\mathit{s}\right)\right)\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill & \approx {\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{C}}^{-1}\left(\mathit{f}\left(\mathit{s}\right)+\mathit{b}\right)-{\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{C}}^{-1}\hfill \\ \hfill & \times \left(\mathit{f}\left(\mathit{s}\right)+\mathit{b}+\mathit{\eta }\left(\mathit{s}\right)\right)\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill & \approx -{\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{C}}^{-1}\mathit{\eta }\left(\mathit{s}\right).\hfill \end{array}$$

(59)

Substituting (59) into (56), we obtain

$$\begin{array}{cc}\hfill {\mathit{C}}_{\mathrm{IWLS}}& \approx E\left\{{\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{C}}^{-1}\mathit{\eta }\left(\mathit{s}\right){\mathit{\eta }}^T\left(\mathit{s}\right){\mathit{C}}^{-1}{\mathit{A}}^T{\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}\right\}\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill & \approx {\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathit{CC}}^{-1}{\mathit{A}}^T{\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill & \approx {\left({\mathit{A}}^T{\mathit{C}}^{-1}\mathit{A}\right)}^{-1}.\hfill \end{array}$$

(62)

Using (34), the above equation can be rewritten as

$${\mathit{C}}_{\mathrm{IWLS}}\approx {\left({\left({\mathit{L}}^{-1}\mathit{A}\right)}^T{\mathsf{\sum }}^{-1}{\mathit{L}}^{-1}\mathit{A}\right)}^{-1}$$

(63)

where

$$\begin{array}{cc}\hfill {\mathit{L}}^{-1}\mathit{A}& =\mathrm{diag}\left(\frac{\parallel {\mathit{r}}_{12}\parallel }{d_{12}\left(\mathit{s}\right)},\frac{\parallel {\mathit{r}}_{13}\parallel }{d_{13}\left(\mathit{s}\right)},\cdots ,\frac{\parallel {\mathit{r}}_{1N}\parallel }{d_{1N}\left(\mathit{s}\right)}\right)\left[\begin{array}{c}{\mathit{u}}_{12}^T\\ {\mathit{u}}_{13}^T\\ ⋮\\ {\mathit{u}}_{1N}^T\end{array}\right]\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill & ={\left[\begin{array}{cccc}\frac{{\mathit{r}}_2-{\mathit{r}}_1}{d_{12}\left(\mathit{s}\right)},& \frac{{\mathit{r}}_3-{\mathit{r}}_1}{d_{13}\left(\mathit{s}\right)},& \cdots & \frac{{\mathit{r}}_N-{\mathit{r}}_1}{d_{1N}\left(\mathit{s}\right)}\end{array}\right]}^T\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill & ={\left[\begin{array}{cccc}\frac{{\mathit{r}}_2-\mathit{s}}{d_{12}\left(\mathit{s}\right)}-\frac{{\mathit{r}}_1-\mathit{s}}{d_{12}\left(\mathit{s}\right)},& \frac{{\mathit{r}}_3-\mathit{s}}{d_{13}\left(\mathit{s}\right)}-\frac{{\mathit{r}}_1-\mathit{s}}{d_{13}\left(\mathit{s}\right)},& \cdots & \frac{{\mathit{r}}_N-\mathit{s}}{d_{1N}\left(\mathit{s}\right)}-\frac{{\mathit{r}}_1-\mathit{s}}{d_{1N}\left(\mathit{s}\right)}\end{array}\right]}^T.\hfill \end{array}$$

(66)

Under Assumptions 2 and 3, i.e.,

$$d_{12}\left(\mathit{s}\right)\approx d_{13}\left(\mathit{s}\right)\approx \cdots \approx d_{1N}\left(\mathit{s}\right)\approx \parallel \mathit{s}-{\mathit{r}}_1\parallel \approx \parallel \mathit{s}-{\mathit{r}}_2\parallel \approx \parallel \mathit{s}-{\mathit{r}}_3\parallel \approx \cdots \approx \parallel \mathit{s}-{\mathit{r}}_N\parallel$$

(67)

Equation (66) can be approximately written as

$$\begin{array}{cc}\hfill {\mathit{L}}^{-1}\mathit{A}& \approx -{\left[\begin{array}{cccc}\frac{{\mathit{r}}_1-\mathit{s}}{\parallel \mathit{s}-{\mathit{r}}_1\parallel }-\frac{{\mathit{r}}_2-\mathit{s}}{\parallel \mathit{s}-{\mathit{r}}_2\parallel },& \frac{{\mathit{r}}_1-\mathit{s}}{\parallel \mathit{s}-{\mathit{r}}_1\parallel }-\frac{{\mathit{r}}_3-\mathit{s}}{\parallel \mathit{s}-{\mathit{r}}_3\parallel },& \cdots & \frac{{\mathit{r}}_1-\mathit{s}}{\parallel \mathit{s}-{\mathit{r}}_1\parallel }-\frac{{\mathit{r}}_N-\mathit{s}}{\parallel \mathit{s}-{\mathit{r}}_N\parallel }\end{array}\right]}^T\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill & \approx -{\mathit{J}}_{\mathit{o}}\hfill \end{array}$$

(69)

where ${\mathit{J}}_{\mathit{o}}$ is the Jacobian matrix in (13). Plugging (69) into (63) finally gives

$${\mathit{C}}_{\mathrm{IWLS}}\approx {\left({\mathit{J}}_{\mathit{o}}^T{\mathsf{\sum }}_{-1}{\mathit{J}}_{\mathit{o}}\right)}^{-1}$$

(70)

which is identical to the CRLB in (12). Thus, we conclude that the IWLS estimator is approximately efficient under the assumptions of small TDOA noise and tightly clustered sensors with roughly the same distance from the emitter (see (67)).

6. Simulation Studies

Monte Carlo (MC) simulations have been carried out to evaluate the performance of the new estimator developed in Section 5 in comparison with the MLE computed using the GN algorithm. For performance comparison, the bias and RMSE of the estimators are considered:

$$\begin{array}{c}\mathrm{Bias}\mathrm{norm}=∥\frac{1}{M}\sum _{k=1}^M({\widehat{\mathit{s}}}_k-\mathit{s})∥\end{array}$$

(71)

$$\begin{array}{c}\mathrm{RMSE}=\sqrt{\mathrm{tr}\frac{1}{M}\sum _{k=1}^M({\widehat{\mathit{s}}}_k-\mathit{s}){({\widehat{\mathit{s}}}_k-\mathit{s})}^T}\end{array}$$

(72)

where M is the number of MC simulation runs, ${\widehat{\mathit{s}}}_k$ is the emitter location estimate at the kth run and tr denotes the matrix trace (sum of diagonal entries). A total of 10,000 MC runs were used in each simulation.

Figure 4 shows the simulated TDOA localization geometry, which resembles a low earth orbit (LEO) satellite localization scenario. The emitter (LEO satellite) is positioned at $\mathit{s}={[400,500,300]}^T$ km. A collective of seven sensors ($N=7$) are used for localizing the emitter and are placed at the locations ${\mathit{r}}_1={[-100,100,10]}^T$ km, ${\mathit{r}}_2=[-80,80,20]$ km, ${\mathit{r}}_3={[-60,60,0]}^T$ km, ${\mathit{r}}_4={[-40,40,10]}^T$ km, ${\mathit{r}}_5={[-20,20,20]}^T$ km, ${\mathit{r}}_6={[0,30,0]}^T$ km and ${\mathit{r}}_7={[20,40,10]}^T$ km. This is a relatively poor localization geometry as a result of a large emitter range to baseline ratio. However, it is somewhat improved by small $|g_{1i}|/\parallel {\mathit{r}}_{1i}\parallel$, as alluded to in Section 3. For the simulated sensor array geometry, the condition number of $\mathit{A}$ is 9.6, which indicates that the unit vectors ${\mathit{u}}_{1i}$ are linearly independent to a satisfactory extent. The ratios $|g_{1i}|/\parallel {\mathit{r}}_{1i}\parallel$ and the emitter ranges from the sensors are listed in

Table 1. TDOA localization geometry parameters for emitter location $\mathit{s}={[400,500,300]}^T$.

Sensor Index i	$|{\mathit{g}}_{1\mathit{i}}|/\parallel {\mathit{r}}_{1\mathit{i}}\parallel$	$\parallel \mathit{s}-{\mathit{r}}_{\mathit{i}}\parallel$ (km)
1	–	702.9225
2	0.2120	696.5630
3	0.0136	703.7045
4	0.0403	699.4998
5	0.0560	696.5630
6	0.1364	686.2215
7	0.2946	663.4003

. We observe that $|g_{1i}|/\parallel {\mathit{r}}_{1i}\parallel$ are mostly close to zero, which is desirable. As the sensors are closely spaced, the emitter ranges $\parallel \mathit{s}-{\mathit{r}}_i\parallel$ do not vary greatly, thereby satisfying Assumptions 2 and 3, in Section 5.2

We have simulated the bias and RMSE of the new two-stage algorithm and IWLS, developed in Section 5, and the MLE computed using the GN algorithm for three different emitter heights: $\mathit{s}={[400,500,300]}^T$ km, $\mathit{s}={[400,500,500]}^T$ km and $\mathit{s}={[400,500,700]}^T$ km. The simulation results and CRLB, calculated by taking the square root of the trace of the CRLB matrix, are shown in Figure 5, Figure 6 and Figure 7 for RDOA noise standard deviation in the range $0.1\le {\sigma }_n\le 1.1$ km. Considering an emitter at $\mathit{s}={[400,500,500]}^T$ km with carrier frequency 2 GHz (L band), bandwidth 1 MHz and time interval of 1 ms for TDOA measurements, the simulated RDOA noise range approximately corresponds to an effective transmit power of 500 mW at ${\sigma }_n=0.1$ km and 43.5 mW at ${\sigma }_n=1.1$ km, which is practical for LEO satellite localization. To arrive at these emitter transmit powers, we have used a link budget analysis based on the TDOA variance bound derived in [23]. In the simulations, the IWLS is iterated four times. The GN algorithm uses 10 iterations and is initialized to the true target location.

Referring to Figure 5, Figure 6 and Figure 7, it is clear that the new algorithm enjoys a stable performance and relatively small bias which is much lower than the MLE at large noise. Being a nonlinear estimator and therefore subject to the threshold effect, the MLE suddenly loses its stability as the RDOA noise is increased. The IWLS estimator exhibits a large bias. In terms of the RMSE performance, the new algorithm performs similarly to the MLE at small RDOA noise and stays close to the CRLB at large noise. The IWLS estimator and the new algorithm have a similar RMSE performance.

7. Conclusions

A new two-stage TDOA localization algorithm was proposed based on finding the intersection of RDOA cones in the far field, defined by 1D-AOAs readily available from TDOA measurements. The proposed algorithm has low complexity and exhibits good performance compared to the MLE, which is considered to be the benchmark given its asymptotic unbiasedness and efficiency. It is composed of an IWLS estimator incorporating approximate linearization and a Taylor series estimator aimed at refining the IWLS estimate. The IWLS estimator was shown to be approximately efficient, achieving the CRLB, under the condition of small noise and closely spaced sensors with the emitter in the far field. In general, the effectiveness of TDOA localization algorithms strongly depends on the localization geometry, and more specifically the relative sensor–emitter geometry. The new algorithm was shown to perform well, outperforming the MLE, in a large emitter range to baseline ratio scenario, which is known to represent a poor geometry. In general, a TDOA sensor array geometry with relatively small RDOAs compared with intersensor distances and linearly independent intersensor vectors will ensure a good localization performance.