Time Difference of Arrival Source Localization:
Exact Linear Solutions for the General 3D Problem

The time difference of arrival (TDOA) problem admits exact, purely algebraic solutions for the situation in which there are 4 and 5 sensors and a single source whose position is to be determined in 3 dimensions. The solutions are exact in the sense that there is no least squares operation (i.e., projection) involved in the solution. The solutions involve no linearization or iteration, and are algebraically transparent via vector algebra in Cartesian coordinates. The solution with 5 sensors requires no resolution of sign ambiguities; the solution with 4 sensors requires resolution of one sign ambiguity. Solutions are effected using only TDOA and not, e.g., frequency difference of arrival (FDOA) or angle of arrival (AOA) [1].

We note previous work towards achieving closed-form solutions to the TDOA problem in 2 dimensions ([2, 3]) and 3 dimensions ([4]). Our work differs from the former in that it is fully general in 3 dimensions. It differs from the latter in that we demonstrate a purely linear algebraic solution in Cartesian coordinates without auxiliary variables and show clearly how the requirement for ambiguity resolution arises for the 4-sensor case but does not appear in the 5-sensor case.

We first present the 5-sensor solution (Section 2) and then follow with the 4-sensor scenario (Section 3). Numerical experiments (Section 4) are presented showing the performance of the calculations in the case of no noise, before closing with conclusions (Section 5). Performance of the calculations below is exact within numerical error, and in the small fraction of cases in which source localization does not occur, it is driven by misidentification in resolution of sign ambiguity without priors.

Our results are suitable for self-localization applications such as simultaneous localization and mapping (SLAM) and precision navigation and timing (PNT, through, e.g., GPS), as well electromagnetic and acoustic signal source geolocation problems. For their speed, exactness, and linearity, we believe the calculations below have substantial practical utility.

First we consider the 5-sensor scenario. Let the sensors have positions given by $\textbf{r}^{\prime}_{k}$, $k=1,...,5$ and the source position, to be solved for, be given by $\textbf{r}^{\prime}_{S}$. If $\textbf{r}^{\prime}_{1}$ is set to be the origin of the coordinate system [5], then $\textbf{r}_{k}\equiv\textbf{r}^{\prime}_{k}-\textbf{r}^{\prime}_{1}$, $k=2,...,5$, and $\textbf{r}_{1}=0$.

The ranges from the source to the sensors, $\rho_{k}=|\textbf{r}_{k}-\textbf{r}_{S}|$, are unknowns, while the quantities $\rho_{k}-\rho_{j}\equiv\delta_{kj}$ are measured/inferred quantities based on time arrival differences of some waveform. That is, for measurement times $t_{k}$ and $t_{j}$ at the $k\textsuperscript{th}$ and $j\textsuperscript{th}$ sensor respectively, $\rho_{k}-\rho_{j}=\mathrm{c}(t_{k}-t_{j})$ where $\mathrm{c}$ is the speed of light and the $t$s are corrected times of arrival at each sensor. Since $\textbf{r}_{1}=0$, $\rho_{1}=\sqrt{\textbf{r}_{S}^{T}\textbf{r}_{S}}=|\textbf{r}_{S}|$.

We have for $\delta_{k1}$, $k>2$,

$\displaystyle\delta_{k1}=\rho_{k}-\rho_{1}$

$\displaystyle=\sqrt{\textbf{r}_{k}^{T}\textbf{r}_{k}-2\textbf{r}_{k}^{T}\textbf{r}_{S}+\textbf{r}_{S}^{T}\textbf{r}_{S}}-\sqrt{\textbf{r}_{S}^{T}\textbf{r}_{S}}.$

(1)

Upon rearranging and squaring, we get an expression for $\rho_{1}$ that is linear in the source position $\textbf{r}_{S}$:

$\displaystyle\rho_{1}=\frac{\textbf{r}_{k}^{T}\textbf{r}_{k}-2\textbf{r}_{k}^{T}\textbf{r}_{S}-\delta_{k1}^{2}}{2\delta_{k1}}.$

(2)

The expression for $\rho_{1}$ can be substituted into any expression for $\delta_{j1}$ ($j\neq 1,j\neq k$) to yield an expression for $\textbf{r}_{S}$:

$\displaystyle 2\left(\textbf{r}_{k}^{T}-\frac{\delta_{k1}}{\delta_{j1}}\textbf{r}_{j}^{T}\right)\textbf{r}_{S}=-\left(\delta_{k1}^{2}-\frac{\delta_{k1}}{\delta_{j1}}\delta_{j1}^{2}\right)+\left[\textbf{r}_{k}^{T}\textbf{r}_{k}-\frac{\delta_{k1}}{\delta_{j1}}\left(\textbf{r}_{j}^{T}\textbf{r}_{j}\right)\right].$

(3)

The vector transpose multiplying $\textbf{r}_{S}$ on the left hand side of Eq. (3) is built up into a matrix B, while the scalars on the right hand side are built up into a vector x so that $\textbf{B}\textbf{r}_{S}=\textbf{x}$. Since we are solving for a $3\times 1$ array, we need 3 pairs of indices $(k,j)$ that will maintain the full rank of B. We choose $(k,j)=\{(3,2),(4,3),(5,4)\}$, so that explicitly, we have

$\displaystyle\textbf{B}=\left[\begin{array}[]{c}\rule[2.15277pt]{10.76385pt}{0.5pt}\leavevmode\nobreak\ 2\left(\textbf{r}_{3}^{T}-\displaystyle\frac{\delta_{31}}{\delta_{21}}\textbf{r}_{2}^{T}\right)\rule[2.15277pt]{10.76385pt}{0.5pt}\\ \rule[2.15277pt]{10.76385pt}{0.5pt}\leavevmode\nobreak\ 2\left(\textbf{r}_{4}^{T}-\displaystyle\frac{\delta_{41}}{\delta_{31}}\textbf{r}_{3}^{T}\right)\rule[2.15277pt]{10.76385pt}{0.5pt}\\ \rule[2.15277pt]{10.76385pt}{0.5pt}\leavevmode\nobreak\ 2\left(\textbf{r}_{5}^{T}-\displaystyle\frac{\delta_{51}}{\delta_{41}}\textbf{r}_{4}^{T}\right)\rule[2.15277pt]{10.76385pt}{0.5pt}\end{array}\right]$

(7)

and

$\displaystyle\textbf{x}=\left[\begin{array}[]{c}-\left(\delta_{31}^{2}-\displaystyle\frac{\delta_{31}}{\delta_{21}}\delta_{21}^{2}\right)+\left[\textbf{r}_{3}^{T}\textbf{r}_{3}-\displaystyle\frac{\delta_{31}}{\delta_{21}}\left(\textbf{r}_{2}^{T}\textbf{r}_{2}\right)\right]\\ -\left(\delta_{41}^{2}-\displaystyle\frac{\delta_{41}}{\delta_{31}}\delta_{31}^{2}\right)+\left[\textbf{r}_{4}^{T}\textbf{r}_{4}-\displaystyle\frac{\delta_{41}}{\delta_{31}}\left(\textbf{r}_{3}^{T}\textbf{r}_{3}\right)\right]\\ -\left(\delta_{51}^{2}-\displaystyle\frac{\delta_{51}}{\delta_{41}}\delta_{41}^{2}\right)+\left[\textbf{r}_{5}^{T}\textbf{r}_{5}-\displaystyle\frac{\delta_{51}}{\delta_{41}}\left(\textbf{r}_{4}^{T}\textbf{r}_{4}\right)\right]\end{array}\right].$

(11)

Thus,

$\displaystyle\textbf{r}_{S}=\textbf{B}^{-1}\textbf{x}.$

(12)

The solution is exact (so long as sensor geometry ensures B has rank 3) and requires no resolution of sign ambiguities. Note that since the coordinates have been referenced to $\textbf{r}^{\prime}_{1}$, we have the actual source location $\textbf{r}^{\prime}_{S}$ as

$\displaystyle\textbf{r}^{\prime}_{S}=\textbf{r}_{S}+\textbf{r}^{\prime}_{1}.$

(13)

The 4-sensor solution has its starting point in Eq. (2). Taking $k=2,3,4$, we have

	$\displaystyle 2\delta_{21}\rho_{1}$	$\displaystyle=$	$\displaystyle\textbf{r}_{2}^{T}\textbf{r}_{2}-2\textbf{r}_{2}^{T}\textbf{r}_{S}-\delta_{21}^{2}$		(14)
	$\displaystyle 2\delta_{31}\rho_{1}$	$\displaystyle=$	$\displaystyle\textbf{r}_{3}^{T}\textbf{r}_{3}-2\textbf{r}_{3}^{T}\textbf{r}_{S}-\delta_{31}^{2}$		(15)
	$\displaystyle 2\delta_{41}\rho_{1}$	$\displaystyle=$	$\displaystyle\textbf{r}_{4}^{T}\textbf{r}_{4}-2\textbf{r}_{4}^{T}\textbf{r}_{S}-\delta_{41}^{2},$		(16)

or, defining

$\displaystyle\textbf{z}=\left[\begin{array}[]{c}2\delta_{21}\\ 2\delta_{31}\\ 2\delta_{41}\end{array}\right],\leavevmode\nobreak\ \leavevmode\nobreak\ \textbf{y}=\left[\begin{array}[]{c}\textbf{r}_{2}^{T}\textbf{r}_{2}-\delta_{21}^{2}\\ \textbf{r}_{3}^{T}\textbf{r}_{3}-\delta_{31}^{2}\\ \textbf{r}_{4}^{T}\textbf{r}_{4}-\delta_{41}^{2}\end{array}\right],\leavevmode\nobreak\ \leavevmode\nobreak\ \textbf{C}=\left[\begin{array}[]{c}\rule[2.15277pt]{10.76385pt}{0.5pt}-2\textbf{r}_{2}^{T}\rule[2.15277pt]{10.76385pt}{0.5pt}\\ \rule[2.15277pt]{10.76385pt}{0.5pt}-2\textbf{r}_{3}^{T}\rule[2.15277pt]{10.76385pt}{0.5pt}\\ \rule[2.15277pt]{10.76385pt}{0.5pt}-2\textbf{r}_{4}^{T}\rule[2.15277pt]{10.76385pt}{0.5pt}\end{array}\right],$

(26)

we have equivalently

$\displaystyle\rho_{1}\textbf{z}=\textbf{y}+\textbf{C}\textbf{r}_{S}.$

(27)

If the sensor geometry ensures C has rank 3, then we have

$\displaystyle\textbf{r}_{S}=\rho_{1}\textbf{C}^{-1}\textbf{z}-\textbf{C}^{-1}\textbf{y}.$

(28)

Define for convenience $\boldsymbol{\xi}\equiv\textbf{C}^{-1}\textbf{z}$ and $\boldsymbol{\eta}\equiv\textbf{C}^{-1}\textbf{y}$. Then using $\textbf{r}_{S}^{T}\textbf{r}_{S}=\rho_{1}^{2}$ gives a quadratic equation for $\rho_{1}$

$\displaystyle\rho_{1}^{2}\left(\boldsymbol{\xi}^{T}\boldsymbol{\xi}-1\right)-2\boldsymbol{\xi}^{T}\boldsymbol{\eta}\rho_{1}+\boldsymbol{\eta}^{T}\boldsymbol{\eta}=0,$

(29)

with a solution

$\displaystyle\rho_{1}=\frac{\boldsymbol{\xi}^{T}\boldsymbol{\eta}\pm\sqrt{(\boldsymbol{\xi}^{T}\boldsymbol{\eta})^{2}-\left(\boldsymbol{\xi}^{T}\boldsymbol{\xi}-1\right)\boldsymbol{\eta}^{T}\boldsymbol{\eta}}}{\left(\boldsymbol{\xi}^{T}\boldsymbol{\xi}-1\right)}.$

(30)

Once $\rho_{1}$ is solved for, the solution for $\textbf{r}_{S}$ is effected through Eq. (28), and the actual source position via Eq.(13).

The sign ambiguity in Eq. (30) can usually be resolved by taking the two source positions corresponding to the two $\rho_{1}$ solutions ($\textbf{r}_{S}(\rho_{1}^{(1)})$ and $\textbf{r}_{S}(\rho_{1}^{(2)})$, say), recalculating the corresponding $\rho_{k}-\rho_{1}$, and determining the one of two solutions that minimizes $\sum_{k=2}^{4}\left[(\rho_{k}-\rho_{1})-\delta_{k1}\right]^{2}$. The presence of priors (positions constrained, e.g., by fixed beam widths) can also be used to constrain the solution.

The calculations above have been numerically tested with generally excellent, exact performance. Suppose the source location is varied in fully 3-dimensional space ($\pm x$, $\pm y$, and $\pm z$) over 1,000 Monte Carlo instances, and the locations of up to 5 sensors are also varied in 3-dimensional space. The sensor positions $\textbf{r}_{k}^{\prime}$ ($k=1,...,5$) and source position $\textbf{r}_{S}^{\prime}$ are drawn from uniform distributions $\mathcal{U}(0,1)$. The source positions are multiplied by a scaling factor $\mathtt{SOURCE\_SCALE}$ relative to the sensor locations:

	$\displaystyle\textbf{r}_{k}^{\prime}$	$\displaystyle\sim$	$\displaystyle\left\{\left[\mathcal{U}(0,1),\mathcal{U}(0,1),\mathcal{U}(0,1)\right]^{T}-0.5\right\},k=1,...,5$		(31)
	$\displaystyle\textbf{r}_{S}^{\prime}$	$\displaystyle\sim$	$\displaystyle\mathtt{SOURCE\_SCALE}\times\left\{\left[\mathcal{U}(0,1),\mathcal{U}(0,1),\mathcal{U}(0,1)\right]^{T}-0.5\right\}.$		(32)

In each experiment, the inferred source position $\tilde{\textbf{r}}_{S}^{\prime}$ is compared to the truth source position $\textbf{r}_{S}^{\prime}$ with no measurement noise. For each Monte Carlo instance, if the relative error $|\tilde{\textbf{r}}_{S}^{\prime}-\textbf{r}_{S}|/|\textbf{r}_{S}|$ is less than some threshold $\mathcal{T}$, the calculation is considered a success. For a given $\mathtt{SOURCE\_SCALE}$, the fraction of successes over all 1,000 Monte Carlo runs is calculated. In Fig. 1 below, we show results from numerical experiments for cases with $\mathcal{T}=10^{-6},10^{-3}$ and $\mathtt{SOURCE\_SCALE}\in[10^{-6},1]$ .

The very small fraction of errors that does exist is driven by numerical matrix inversion errors for cases in which sensor positions have poor geometric diversity and when the source positions are highly constrained compared to the sensors. For the 4-sensor solution, errors that exist for $\mathtt{SOURCE\_SCALE}$ have been verified to be driven entirely by the incorrect choice in $\rho_{1}$ when solving Eq. (30). Hence, in all cases (4- or 5-sensor), the inferred solution for no noise is within numerical error of the exact value, and even in incorrect inferences in the limiting 4-sensor case, the correct solution is still one of two tested values that may potentially be easily selected based on priors. Depending on the application, the existence of this error may therefore have negligible practical impact.

We have presented here fully algebraic solutions to the TDOA problem. For the case of 5 sensors, the solution method requires only inverting a set of 3 linear equations. For the case of 4 sensors, the solution method requires inverting a set of 3 linear equations and solving one quadratic equation. Numerical experiments have been carried out using the techniques derived for varying sensor and source positions, each over 1,000 Monte Carlo instances with excellent, exact performance in the noiseless case. To the best of our knowledge, this is the first presentation of exact, linear solutions to the TDOA problem for 4- or 5-sensor scenarios in 3 dimensions.

That 4 sensors require ambiguity resolution while 5 sensors do not is perhaps evident from geometrical considerations, although the presence of a purely linear algebraic solution to the 5-sensor case might be somewhat unexpected. For 2 sensors in 3 dimensions, the TDOA problem possess a degeneracy about the surface of a hyperboloid. Introducing a third sensor reduces this degeneracy to a line (the intersection between two hyperboloids; Fig. 2), while a fourth sensor reduces the source location to up to two points, depending on the shape of the line and the new hyperboloid. The introduction of a fifth sensor ensures that if there are two possible solutions from the 4-sensor case, only the viable one is selected.

The fact that our solutions are linear make them amenable to rapid, on-platform calculation. The lack of need for iteration obviates concerns about convergence. We believe there is substantial practical utility in implementing these solutions for a range of localization and tracking problems.