Hybrid Multilateration and Triangulation

Eric Widdison, David G. Long

Department of Electrical and Computer Engineering, Brigham Young University, Provo, USA.

DOI: 10.4236/pos.2021.121001                

Abstract

Wide area multilateration algorithms suffer from stability issues related to the fact that the reference points are nearly coplanar. This paper presents a method to add elevation angle measurements to a multilateration problem and thereby reduce the error perpendicular to the plane where the measurements are taken. The resulting measurement error is significantly reduced for co- planar and nearly coplanar reference points.

Keywords

Localization, Multilateration, Linear Least Square

1. Introduction

One of the goals in aviation is to detect and track aircraft in flight. This is traditionally done using radar, which measures the distance and direction to a target from a known position ( [1], p. 3). The use of radar presents a number of limitations. Radar requires high power transmitters, large antennas, and sensitive receivers. It is dependent on receiving a reflected signal that is strong enough to detect the target.

One method to overcome some of these limitations is to use signals transmitted by the aircraft to determine its position. By using a number of receivers at known points, relative distances to the target can be determined, from which the position of the aircraft can be calculated. This approach is called multilateration.

In physical multilateration implementations with ground-based receivers spread over a wide area, the target and receivers are often nearly coplanar. This arrangement is numerically unstable in multilateration algorithms and typically produces significant error in the vertical position of the aircraft.

This paper presents a hybrid multilateration (HM) algorithm for incorporating triangulation measurements into a multilateration calculation to improve the accuracy of the calculated positions. This is done by adding measurements of the elevation angle of the target at each receiver and using those measurements to triangulate the vertical position of the aircraft while using multilateration to determine its horizontal position. Elevation angle can be determined using a linear monopulse array [2], phase difference measurements, time difference of arrival measurements [3] [4], or any other method ( [1], p. 471, 474) as long as it provides sufficient resolution for this algorithm. This paper focuses on the derivation of the HM algorithm rather than on the methods of collecting the measurements used in the algorithm. The algorithm does not depend on the method of data collection.

The development of the HM algorithm starts with the derivation of an algorithm that only uses multilateration and similar terms for using triangulation data. The multilateration and triangulation equations are combined to produce hybrid equations that use both multilateration and triangulation information to more accurately determine the position of a transmitter. The effects of errors and receiver location are analyzed.

The HM algorithm improves the accuracy of position calculations for scenarios with receivers that are coplanar or nearly coplanar. This allows for more precise measurement of the 3-dimensional position of aircraft in flight over a wide area. Most wide area multilateration systems only measure the position of objects on the ground or measure the ground position of targets in flight. This is done using ground-based sensors that are approximately coplanar and are therefore unable to distinguish between aircraft above or below the ground. Methods to measure the altitude of aircraft in flight through passive signal detection typically rely on triangulation only [5] [6].

2. Algorithm

Let $P_i=\left(x_i\mathrm{,}{y}_i\mathrm{,}{z}_i\right)\mathrm{,}i=\mathrm{0,1,}\cdots \mathrm{,}N-1$ be the position of N receivers and $P_s=\left(x_s\mathrm{,}{y}_s\mathrm{,}{z}_s\right)$ be the position of a transmitter. The transmitter emits a signal at time $t_s$ that is received at $P_i$ at time $t_i$. For each receiver, a pseudo-dis- tance $d_i=c{t}_i$ is calculated, where c is the speed of light. This pseudodistance is used instead of time in the multilateration equations. For each $P_i$, the angle of arrival of the transmitted signal relative to the x-y plane, ${\theta }_i$ is measured. Variables with two subscripts and a tilde are the difference between that parameter for the two indexed points, i.e. ${\tilde{x}}_{i\mathrm{,}j}=x_i-x_j$, ${\tilde{y}}_{i\mathrm{,}j}=y_i-y_j$, ${\tilde{z}}_{i\mathrm{,}j}=z_i-z_j$, and ${\tilde{d}}_{i\mathrm{,}j}=d_i-d_j$.

The following sections present the HM algorithm by first deriving a multi- lateration algorithm, and then incorporating the measured angle information.

2.1. Multilateration Algorithm

Given pseudodistances between fixed points and an unknown transmitter location, the possible solution space is the intersection of a number of hyperbolic surfaces. To make an algorithm that is efficient and stable the goal is to reduce the hyperbolic intersection problem into a system of linear equations.

The first step in producing the HM algorithm is to derive an algorithm for multilateration using the given measured pseudodistances. This is done by eliminating the unknown quadratic terms from the pseudodistance equations, leaving linear equations in terms of the unknown transmitter location.

The distance between points $P_i$ and $P_s$ is

$\begin{array}{c}{\tilde{d}}_{i\mathrm{,}s}=\Vert P_i-P_s\Vert =\sqrt{{\left(x_i-x_s\right)}^2+{\left(y_i-y_s\right)}^2+{\left(z_i-z_s\right)}^2}\\ =c\left(t_i-t_s\right)=d_i-d_s\mathrm{.}\end{array}$ (1)

Squaring both sides gives

${\tilde{d}}_{i\mathrm{,}s}^2=d_i^2-2d_i{d}_s+d_s^2=x_i^2-2x_i{x}_s+x_s^2+y_i^2-2y_i{y}_s+y_s^2+z_i^2-2z_i{z}_s+z_s^2$ (2)

which can be rearranged to yield

$d_i^2=x_i^2-2x_i{x}_s+x_s^2+y_i^2-2y_i{y}_s+y_s^2+z_i^2-2z_i{z}_s+z_s^2+2d_i{d}_s-d_s^2\mathrm{.}$ (3)

Taking the difference between the squared pseudodistance for two points $P_i$ and $P_j$ gives

$\begin{array}{c}{d}_i^2-d_j^2=x_i^2-x_j^2-2x_i{x}_s+2x_j{x}_s+y_i^2-y_j^2-2y_i{y}_s+2y_j{y}_s\\ +z_i^2-z_j^2-2z_i{z}_s+2z_j{z}_s+2d_i{d}_s-2d_j{d}_s\mathrm{.}\end{array}$ (4)

This can be rearranged to produce a linear equation in terms of four unknowns, $x_s$, $y_s$, $z_s$, and $d_s$.

$\begin{array}{l}{x}_i^2-x_j^2+y_i^2-y_j^2+z_i^2-z_j^2-d_i^2+d_j^2\\ =2\left(x_i-x_j\right)x_s+2\left(y_i-y_j\right)y_s+2\left(z_i-z_j\right)z_s-2\left(d_i-d_j\right)d_s\mathrm{.}\end{array}$ (5)

Observing that $x_i^2-x_j^2+y_i^2-y_j^2+z_i^2-z_j^2={\Vert P_i\Vert }^2-{\Vert P_j\Vert }^2$ and using the difference forms of the variables, Equation (5) can be rewritten as

${\Vert P_i\Vert }^2-{\Vert P_j\Vert }^2-d_i^2+d_j^2=2{\tilde{x}}_{i\mathrm{,}j}{x}_s+2{\tilde{y}}_{i\mathrm{,}j}{y}_s+2{\tilde{z}}_{i\mathrm{,}j}{z}_s-2{\tilde{d}}_{i\mathrm{,}j}{d}_s\mathrm{.}$ (6)

Normally, one point is chosen as the reference point and used as $P_j$. For this derivation let $P_0$ be the reference point. Applying Equation (6) across at least 4 independent pairs of sensors gives a multilateration matrix equation

$A_M{x}_s=b_M$ (7)

where

$\begin{array}{l}{A}_M=\left[\begin{array}{cccc}2{\tilde{x}}_{\mathrm{1,0}}& 2{\tilde{y}}_{\mathrm{1,0}}& 2{\tilde{z}}_{\mathrm{1,0}}& -2{\tilde{d}}_{\mathrm{1,0}}\\ 2{\tilde{x}}_{\mathrm{2,0}}& 2{\tilde{y}}_{\mathrm{2,0}}& 2{\tilde{z}}_{\mathrm{2,0}}& -2{\tilde{d}}_{\mathrm{2,0}}\\ ⋮& ⋮& ⋮& ⋮\\ 2{\tilde{x}}_{N-\mathrm{1,0}}& 2{\tilde{y}}_{N-\mathrm{1,0}}& 2{\tilde{z}}_{N-\mathrm{1,0}}& -2{\tilde{d}}_{N-\mathrm{1,0}}\end{array}\right]\\ x_s={\left[\begin{array}{cccc}{x}_s& y_s& z_s& d_s\end{array}\right]}^{\text{T}}\\ b_M=\left[\begin{array}{c}{\Vert P_1\Vert }^2-{\Vert P_0\Vert }^2-d_1^2+d_0^2\\ {\Vert P_2\Vert }^2-{\Vert P_0\Vert }^2-d_2^2+d_0^2\\ ⋮\\ {\Vert P_{N-1}\Vert }^2-{\Vert P_0\Vert }^2-d_{N-1}^2+d_0^2\end{array}\right]\mathrm{.}\end{array}$ (8)

Using the method of least squares, the estimate of the position of the transmitter $x_s$ is

$x_s={\left(A_M^{\text{T}}{A}_M\right)}^{-1}{A}_M^{\text{T}}{b}_M\mathrm{.}$ (9)

2.2. Triangulation

The triangulation algorithm uses the measured values of the angle of the transmitter relative to vertical for each receiver, as shown in Figure 1. Adding this one-dimensional triangulation to the multilateration equations can help to significantly improve the solution in the z direction.

The elevation angle to the target at point $P_i$ is ${\theta }_i$. This angle can be measured using a large, directional antenna such as a monopulse radar antenna, or by using time difference of arrival or received phase difference for two vertically aligned antennas. Recall that this paper focuses on the algorithm to use the measured angle to estimate height rather than the method of obtaining the measurements.

For a given point and angle, the vertical position of the target can be given as

$z_s=z_i+\left(d_i-d_s\right)\cos {\theta }_i$ (10)

which can be rearranged to give a linear equation in two unknowns,

$z_s+d_s\cos {\theta }_i=z_i+d_i\cos {\theta }_i\mathrm{.}$ (11)

This gives a triangulation matrix equation $A_T{x}_s=b_T$ with terms

$\begin{array}{l}{A}_T=\left[\begin{array}{cccc}0& 0& 1& \cos {\theta }_0\\ 0& 0& 1& \cos {\theta }_1\\ ⋮& ⋮& ⋮& ⋮\\ 0& 0& 1& \cos {\theta }_{N-1}\end{array}\right]\\ b_T=\left[\begin{array}{c}{z}_0+d_0\cos {\theta }_0\\ z_1+d_1\cos {\theta }_1\\ ⋮\\ z_{N-1}+d_{N-1}\cos {\theta }_{N-1}\end{array}\right]\mathrm{.}\end{array}$ (12)

Figure 1. Angle measurement for triangulation algorithm. The measured angle describes a cone containing the transmitter.

With two columns that are all zero, the equation $A_T{x}_s=b_T$ is clearly deficient and cannot be used by itself to solve for $x_s$. It can be combined with $A_M$ and $b_M$, which are not deficient, to produce a linear equation $A_{MT}{x}_s=b_{MT}$ using augmented matrices

$A_{MT}=\left[\begin{array}{c}{A}_M\\ w{A}_T\end{array}\right]b_{MT}=\left[\begin{array}{c}{b}_M\\ w{b}_T\end{array}\right]$ (13)

where w is a weighting factor. The weighting factor is included to control the relative contribution of the triangulation equations to the final solution. Using $w=1$ would result in the multilateration terms in Equation (8) dominating the solution so that the triangulation terms do not contribute to the result. If w is too large then the triangulation terms will dominate the solution, skewing the x and y values to correspond to the z value determined by the triangulation terms. For the best result the weight w should be chosen so that the magnitude of the terms in $A_M$ and $w{A}_T$ are approximately the same size. The coefficients of $A_T$ are between 0 and 1 so reasonable values for w would be the same size as the coefficients in $A_M$. The coefficients have the same magnitude as the average distance between the receivers $P_i$ or the average distance from the receivers to the transmitter1.

2.3. Hybrid Multilateration/Triangulation

Equation (6) can be rewritten so that the x and y dimensions are solved using multilateration and the z dimension is solved solely by triangulation. This is done by substituting Equation (10) into Equation (6),

$\begin{array}{l}{\Vert P_i\Vert }^2-{\Vert P_0\Vert }^2-d_i^2+d_0^2\\ =2{\tilde{x}}_{i\mathrm{,0}}{x}_s+2{\tilde{y}}_{i\mathrm{,0}}{y}_s+2{\tilde{z}}_{i\mathrm{,0}}\left(z_i+\left(d_i-d_s\right)\cos {\theta }_i\right)-2{\tilde{d}}_{i\mathrm{,}j}{d}_s\mathrm{,}\end{array}$ (14)

which can be expressed in linear terms of the unknown variables as

$\begin{array}{l}{\Vert P_i\Vert }^2-{\Vert P_0\Vert }^2-d_i^2+d_0^2-2{\tilde{z}}_{i\mathrm{,0}}\left(z_i+d_i\cos {\theta }_i\right)\\ =2{\tilde{x}}_{i\mathrm{,0}}{x}_s+2{\tilde{y}}_{i\mathrm{,0}}{y}_s-2\left({\tilde{z}}_{i\mathrm{,0}}\cos {\theta }_i+{\tilde{d}}_{i\mathrm{,0}}{d}_s\right)\mathrm{.}\end{array}$ (15)

Equation (15) is a linear equation of the unknown terms, which can be combined with Equation (11) and expressed as

$A_H\mathrm{<mo>\; =\; </mo>}\left[\begin{array}{cccc}2{\tilde{x}}_{\mathrm{1,0}}& 2{\tilde{y}}_{\mathrm{1,0}}& 0& -2\left({\tilde{z}}_{\mathrm{1,0}}\cos {\theta }_1+{\tilde{d}}_{\mathrm{1,0}}\right)\\ 2{\tilde{x}}_{\mathrm{2,0}}& 2{\tilde{y}}_{\mathrm{2,0}}& 0& -2\left({\tilde{z}}_{\mathrm{2,0}}\cos {\theta }_2+{\tilde{d}}_{\mathrm{2,0}}\right)\\ ⋮& ⋮& ⋮& ⋮\\ 2{\tilde{x}}_{N-\mathrm{1,0}}& 2{\tilde{y}}_{N-\mathrm{1,0}}& 0& -2\left({\tilde{z}}_{N-\mathrm{1,0}}\cos {\theta }_{N-1}+{\tilde{d}}_{N-\mathrm{1,0}}\right)\\ 0& 0& 1& \cos {\theta }_1\\ 0& 0& 1& \cos {\theta }_2\\ ⋮& ⋮& ⋮& ⋮\\ 0& 0& 1& \cos {\theta }_{N-1}\end{array}\right]$

$b_H=\left[\begin{array}{c}{\Vert P_1\Vert }^2-{\Vert P_0\Vert }^2-d_1^2+d_0^2-2z_{\mathrm{1,0}}\left(z_1+d_1\cos {\theta }_1\right)\\ {\Vert P_2\Vert }^2-{\Vert P_0\Vert }^2-d_2^2+d_0^2-2z_{\mathrm{2,0}}\left(z_2+d_2\cos {\theta }_2\right)\\ ⋮\\ {\Vert P_{N-1}\Vert }^2-{\Vert P_0\Vert }^2-d_{N-1}^2+d_0^2-2z_{N-\mathrm{1,0}}\left(z_{N-1}+d_{N-1}\cos {\theta }_{N-1}\right)\\ z_{\mathrm{1,0}}+d_{\mathrm{1,0}}\cos {\theta }_1\\ z_{\mathrm{2,0}}+d_{\mathrm{2,0}}\cos {\theta }_2\\ ⋮\\ z_{N-\mathrm{1,0}}+d_{N-\mathrm{1,0}}\cos {\theta }_{N-1}\end{array}\right]$ (16)

and

$A_H{x}_s=b_H$ (17)

the solution of which is

$x_s={\left(A_H^{\text{T}}{A}_H\right)}^{-1}{A}_H^{\text{T}}{b}_H\mathrm{.}$ (18)

3. Stability

Multilateration algorithms are sensitive to the positions where the measurements are taken. If the reference point locations are not sufficiently distributed, then the algorithm can fail. This is a problem with terrestrial multilateration systems that are tracking aircraft over a large area. In a system, the reference points are likely to be approximately coplanar, which makes it difficult for a multilateration algorithm to distinguish between points that are at orthogonally opposite positions relative to that plane.

In analyzing the algorithm, it is assumed that the reference points are coplanar or nearly coplanar. For most scenarios the transmitter is not in the same plane. The stability of the approach depends on the stability of the matrix ${\left(A^{\text{T}}A\right)}^{-1}$ given by

$\begin{array}{l}{A}^{\text{T}}A\\ =4{\displaystyle \underset{i=1}{\overset{N}{\sum }}}\left[\begin{array}{cccc}{\tilde{x}}_{i\mathrm{,0}}^2& {\tilde{x}}_{i\mathrm{,0}}{\tilde{y}}_{i\mathrm{,0}}& 0& -{\tilde{x}}_{i\mathrm{,0}}\left({\tilde{d}}_{i\mathrm{,0}}+{\tilde{z}}_{i\mathrm{,0}}{C}_i\right)\\ {\tilde{x}}_{i\mathrm{,0}}{\tilde{y}}_{i\mathrm{,0}}& {\tilde{y}}_{i\mathrm{,0}}^2& 0& -{\tilde{y}}_{i\mathrm{,0}}\left({\tilde{d}}_{i\mathrm{,0}}+{\tilde{z}}_{i\mathrm{,0}}{C}_i\right)\\ 0& 0& 0& 0\\ -{\tilde{x}}_{i\mathrm{,0}}\left({\tilde{d}}_{i\mathrm{,0}}+{\tilde{z}}_{i\mathrm{,0}}{C}_i\right)& -{\tilde{y}}_{i\mathrm{,0}}\left({\tilde{d}}_{i\mathrm{,0}}+{\tilde{z}}_{i\mathrm{,0}}{C}_i\right)& 0& {\left({\tilde{d}}_{i\mathrm{,0}}+{\tilde{z}}_{i\mathrm{,0}}{C}_i\right)}^2\end{array}\right]\\ +{\displaystyle \underset{i}{\sum }}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 1& \cos {\theta }_i\\ 0& 0& \cos {\theta }_i& {\cos }^2{\theta }_i\end{array}\right]\end{array}$ (19)

for the ${\tilde{x}}_{i\mathrm{,0}}=x_i-x_0$, ${\tilde{y}}_{i\mathrm{,0}}=y_i-y_0$, ${\tilde{z}}_{i\mathrm{,0}}=z_i-z_0$, ${\tilde{d}}_{i\mathrm{,0}}=d_i-d_0$, and $C_i=\cos {\theta }_i$ values used.

The following sections analytically demonstrate the instability of the algorithm for a few specific special cases of receiver layouts. It also considers several other receiver configurations to gain insight into the stability of the algorithm.

3.1. Linear Deployment

If the receivers are in a line where $x_i=x_j\forall i\mathrm{,}j$ then the first row and column of $A^{\text{T}}A$ are uniformly zero, and $A^{\text{T}}A$ has no inverse. The same holds for when $y_i=y_j\forall i\mathrm{,}j$, which leads to all zeros in the second row and column.

For other linear arrangements, consider a line of receivers located at $\left(x_i\mathrm{,}{y}_i\mathrm{,}{z}_i\right)=\left(a{k}_i+x_0\mathrm{,}b{k}_i+x_0\mathrm{,}{z}_i\right)$ for some $a\mathrm{,}b\mathrm{,}{k}_1\mathrm{,}\cdots \mathrm{,}{k}_n$ where $k_i$ is a scalar representing the location of the receiver on the line. Then ${\tilde{x}}_{i\mathrm{,0}}=a{k}_i$ and ${\tilde{y}}_{i\mathrm{,0}}=b{k}_i$. Let $K={\displaystyle \sum }{k}_i^2$. Then

$A^{\text{T}}A=\left[\begin{array}{cccc}4a^{2}K& 4abK& 0& -{\displaystyle {\sum }_{i}a{k}_i{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}}\\ 4abK& 4b^{2}K& 0& -{\displaystyle {\sum }_{i}b{k}_i{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}}\\ 0& 0& N& {\displaystyle {\sum }_i\cos {\theta }_i}\\ -{\displaystyle {\sum }_{i}a{k}_i{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}}& -{\displaystyle {\sum }_{i}b{k}_i{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}}& {\displaystyle {\sum }_i\cos {\theta }_i}& {\displaystyle {\sum }_i\left({\cos }^2{\theta }_i+{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}\right)}\end{array}\right]$ (20)

where ${{\tilde{d}}^{\prime }}_{i\mathrm{,0}}={\tilde{d}}_{i\mathrm{,0}}+{\tilde{z}}_{i\mathrm{,0}}\cos {\theta }_i$. The first two columns of that matrix are $ax$ and $bx$ where

$x={\left[\begin{array}{cccc}4aK& 4bK& 0& -{\displaystyle {\sum }_i{k}_i{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}}\end{array}\right]}^{\text{T}}\mathrm{.}$ (21)

This means that those two columns are not linearly independent and that $A^{\text{T}}A$ is deficient. Therefore any linear arrangement of sensors produces a deficient $A^{\text{T}}A$ matrix.

3.2. Circular Deployment

Let the receivers be located on a circle, such that the location of each sensor is $\left(\hat{x}+r\cos {\varphi }_i\mathrm{,}\hat{y}+r\sin {\varphi }_i\mathrm{,}\hat{z}\right)$, and the transmitter located at $\left(\hat{x}\mathrm{,}\hat{y}\mathrm{,}{z}_s\right)$. The values of ${\tilde{d}}_{i\mathrm{,}s}=\sqrt{r^2+{\left(z-\hat{z}\right)}^2}+d_s$ are the same for all of the sensors, so that ${\tilde{d}}_{i\mathrm{,0}}=0$. Similarly, ${\tilde{z}}_{i\mathrm{,0}}=0$ for all sensors. Since the horizontal and vertical distance from the transmitter to each sensor is the same,

${\theta }_i={\tan }^{-1}\frac{z_s-\hat{z}}{r}$

is the same for all sensors. This results in

$\begin{array}{c}{A}^{\text{T}}A={\displaystyle \underset{i=1}{\overset{N}{\sum }}}\left[\begin{array}{cccc}4{\tilde{x}}_{i\mathrm{,0}}^2& 4{\tilde{x}}_{i\mathrm{,0}}{\tilde{y}}_{i\mathrm{,0}}& 0& -{\tilde{x}}_{i\mathrm{,0}}{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}\\ 4{\tilde{x}}_{i\mathrm{,0}}{\tilde{y}}_{i\mathrm{,0}}& 4{\tilde{y}}_{i\mathrm{,0}}^2& 0& -{\tilde{y}}_{i\mathrm{,0}}{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}\\ 0& 0& 1& \cos {\theta }_i\\ -{\tilde{x}}_{i\mathrm{,0}}{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}& -{\tilde{y}}_{i\mathrm{,0}}{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}& \cos {\theta }_i& {\cos }^2{\theta }_i+{{\tilde{d}}^{\prime }}_{i\mathrm{,0}}^2\end{array}\right]\\ ={\displaystyle \underset{i=1}{\overset{N}{\sum }}}\left[\begin{array}{cccc}4{\tilde{x}}_{i\mathrm{,0}}^2& 4{\tilde{x}}_{i\mathrm{,0}}{\tilde{y}}_{i\mathrm{,0}}& 0& 0\\ 4{\tilde{x}}_{i\mathrm{,0}}{\tilde{y}}_{i\mathrm{,0}}& 4{\tilde{y}}_{i\mathrm{,0}}^2& 0& 0\\ 0& 0& 1& \cos {\theta }_i\\ 0& 0& \cos {\theta }_i& {\cos }^2{\theta }_i\end{array}\right]\\ =\left[\begin{array}{cccc}{\displaystyle {\sum }_{i=1}^{N}4{\tilde{x}}_{i\mathrm{,0}}^2}& {\displaystyle {\sum }_{i=1}^{N}4{\tilde{x}}_{i\mathrm{,0}}{y}_{i\mathrm{,0}}}& 0& 0\\ {\displaystyle {\sum }_{i=1}^{N}4{\tilde{x}}_{i\mathrm{,0}}{y}_{i\mathrm{,0}}}& {\displaystyle {\sum }_{i=1}^{N}4{\tilde{y}}_{i\mathrm{,0}}^2}& 0& 0\\ 0& 0& N& N\cos {\theta }_0\\ 0& 0& N\cos {\theta }_0& N{\cos }^2{\theta }_0\end{array}\right]\end{array}$ (22)

where ${{\tilde{d}}^{\prime }}_{i\mathrm{,0}}={\tilde{d}}_{i\mathrm{,0}}+{\tilde{z}}_{i\mathrm{,0}}\cos {\theta }_i=0$ for all i. The third and fourth columns of that matrix are equal to $N{\left[\begin{array}{cccc}0& 0& 1& \cos {\theta }_0\end{array}\right]}^{\text{T}}$ and $N\cos {\theta }_0{\left[\begin{array}{cccc}0& 0& 1& \cos {\theta }_0\end{array}\right]}^{\text{T}}$, so the columns of $A^{\text{T}}A$ are linearly dependent and $A^{\text{T}}A$ is deficient.

Therefore, an arrangement of receivers on a circle cannot produce a solution when the transmitter is located exactly equidistant from the receivers at the center of the circle. This holds whether the receivers are evenly or unevenly spaced, and whether the receivers are around the whole circle or located exclusively along an arc of the circle. This pole of instability is only present when the transmitter is at the center of the circle. As the transmitter moves away from the center of the circle the instability disappears.

The behavior of the HM algorithm with a circular configuration of receivers can be seen in Figure 2. The HM algorithm is compared to a similar algorithm that only uses multilateration (ML) [7]. Both plots show the deviation as a function of transmitter position relative to the receivers at the positions shown. Gaussian error with a standard deviation of 1 m was added to the true pseudodistances and a Gaussian error with a standard deviation of 1˚ was added to the true elevation values2. Note that the scales are different for the two images.

This shows a problem with using a multilateration algorithm in that scenario that needs to be corrected to give realistic comparable results between the two algorithms. A comparison of the scales in Figure 2 shows that HM algorithm out-performs the ML algorithm by at least 60 dB meters, which is a factor of 1,000,000:1. This is due to the fact that the receivers were coplanar. In any algorithm that relies solely on multilateration, having all of the reference points in a single plane creates an ambiguity. The algorithm cannot distinguish between points above the plane and points below the plane. To resolve the simulations were performed using receiver locations with ${\sigma }_z=20\text{m}$ of Gaussian random variation added to their height. This is the near-coplanar configuration used in all of the remaining simulations in this paper. The results are shown in Figure 3. The HM algorithm still produces significantly lower error except at its point of instability.

The choice of ${\sigma }_z=20\text{m}$ is based on the underlying scenario that the ground level varies by less than 20 m. The performance of both algorithms is affected by this vertical perturbation of receiver positions, but for small variation of height only the ML algorithm is significantly affected. A plot of the position error as a function of ${\sigma }_z$ is shown in Figure 4. Accuracy of the ML algorithm increases as the vertical separation between receivers increases. The HM algorithm performs equally well across a range of local topological variation. This means that the HM algorithm performs better than the ML algorithm when the receivers are more planar and the ML performs better when the receivers have significant vertical variation. In the scenario analyzed in Figure 4 the vertical relief must have a standard deviation of about ${\sigma }_z=500\text{meters}$ (27 dBmeters) before the two algorithms perform equally well. The near-coplanar simulations in this paper use vertical standard deviation of ${\sigma }_z=20\text{meters}$ (13 dBmeters).

Figure 2. Standard deviation of a multilateration (ML) algorithm and the hybrid multi- lateration (HM) algorithm error as a function of transmitter location with coplanar receivers. The ML algorithm has sizeable error with all receivers on the same plane. It cannot distinguish between transmitters with a positive z coordinate and those with a negative z coordinate. Note that the two figures use different scales for the standard deviation of the error.

Figure 3. Standard deviation of ML and HM algorithm error as a function of transmitter location with near-coplanar receivers. The ML error is about 6 to 10 times larger at all points apart from the center where the HM algorithm is unstable.

Figure 4. Algorithm comparison as a function of vertical variance using a central receiver in the same configuration shown in Figure 6 with the transmitter directly above the central receiver.

Moving one of the sensors outward from the circle by a small distance creates two poles of instability along the axis that contains the point that was moved. Moving one sensor inward by a small distance creates two poles of instability that are located on a line perpendicular to the axis that contains the point that was moved. Plots of the convergence are shown in Figure 5. This shows that small changes from a strictly circular configuration aren’t enough to remove the instability.

The circular non-convergence problem can be resolved by moving one of the sensors to the center of the circle and rearranging the other sensors evenly around a circle. This eliminates the pole of instability at the center of the circle, as shown in Figure 6. The HM algorithm performs better than the ML algorithm in this configuration, producing errors less than 5% of the ML algorithm errors at most points.

3.3. Other Receiver Configurations

Circular or near-circular receiver arrangements are undesirable because they create a pole of instability. This instability isn’t present in other receiver configurations. This section considers the stability of the HM algorithm with other receiver configurations.

Figure 5. Standard deviation of the HM algorithm error. The right-most point is displaced 400 m in the x-direction from the circle containing the other points. This shows that the HM algorithm still has points of instability when small changes are made to the circular receiver configuration.

Figure 6. Comparison of the standard deviation of the ML and HM algorithm error for near-coplanar receivers with a central receiver. The HM algorithm no longer has a point of instability. The ML error is 10 times larger than the HM error at every point.

A delta, or V, configuration places the receivers along two lines that meet at a point. This works relatively well with the ML algorithm, which performs roughly the same as it did in the circular configurations. With the HM algorithm this configuration induces a line of relative instability in the direction that the V is pointing, as shown in Figure 7. The HM algorithm still out-performs the ML algorithm, producing about 10% as much error even in the places where it was not performing as well.

Another configuration is a wye or Y configuration. This consists of a central point and the receivers placed on three lines radiating outward from the center. The Y configuration is very similar to the circular configuration with a receiver in the center, as seen in Figure 8. The HM algorithm produces significantly less error at every point.

The performance of the scenarios in this section and Figure 9 is summarized in Table 1. Note that these ranges include the values from outside the sensor perimeter.

4. Implementation Factors

A number of factors should be taken into consideration when implementing the HM algorithm. They include the desired coverage area, the number and location of receivers, and the utilization of those receivers in the calculations. Each of these are addressed in the following section.

Table 1. Range of error in the scenarios presented in this section. All values are standard deviation of error and are given in dB meters.

Figure 7. Comparison of the standard deviation of the ML and HM algorithm error for near-coplanar receivers in a delta configuration. The HM algorithm performs well over most of the area, but is relatively unstable as the transmitter moves away from the receivers in one direction. It still outperforms the ML algorithm at every point.

Figure 8. Comparison of the standard deviation of the ML and HM algorithm error for near-coplanar receivers in a Y configuration. The HM algorithm performs better at every point.

Figure 9. Comparison of the standard deviation of the ML and HM algorithm error for near-coplanar receivers in a triangular lattice configuration. The patterns in the ML figure are a result of the vertical deviation that was added to the receiver locations. Both algorithms show an increase in error as the transmitter moves outside the region containing the receivers.

The HM algorithm works best when the transmitter is located somewhere between the receivers. As such, the optimal arrangement is to have receivers surrounding the area of interest.

4.1. Coverage Area

The coverage area places some requirements or constraints on the employment of this algorithm. The algorithm is built around the assumption that all of the angles are measured relative to a common vertical direction. This can be a problem if the coverage area gets wide enough for the curvature of the Earth to start introducing error into the angle measurements. The two solutions to this are to limit the coverage area or to orient the receivers to a shared vertical that may be different from the local vertical.

The alternative is to orient all the sensors so that they measure angle with respect to parallel vertical axes. This presents a significant challenge installing and maintaining the receiver antennas as it is relatively easy to measure whether an object is aligned to vertical and difficult to align it to a specific deviation from vertical. In a wide area deployment the receivers need to be aligned to an angle that is slightly off from vertical or bias is introduced to the resulting measured angles which produces more position error from the algorithm.

The size of a limited coverage area affects the accuracy of the calculated position. An inherent problem with triangulation is that triangulation errors are proportional to the error in measuring the angle and to the distance from the vertex of the measured angles. This triangulation error must be considered when determining the size of the coverage area and the spacing of sensors within it. A denser deployment of sensors allows for more angular error because the distance from transmitter to receiver does not contribute as much to the overall triangulation error. Receivers deployed less than 44 km apart have less than 1˚ difference in their local vertical direction which may be sufficient for a given application, depending on the tolerance for angular error.

Coverage area can also be expanded by using a wider array of receivers and then only including the closest measurement to the transmitter, determined by the time of arrival of the signal, in the calculations. This keeps the triangulation error in any given calculation small while allowing for a larger coverage area at the expense of requiring more receivers.

4.2. Receiver Positions

The stability analysis in Section 1 demonstrated that the algorithm works best when the transmitter is within a polygon surrounding the receivers. Some coverage outside that boundary may be possible, any positions calculated outside the perimeter contain additional error due to the geometric problems with multilateration and the fact that a fixed angular measurement produces more error the farther it is from the point of measurement.

Spacing of sensors affects the accuracy of the algorithm. The inclusion of triangulation in the algorithm means that more distance between the transmitter and receiver induces a larger position error for the same angle error. The optimal spacing depends on the accuracy of angular measurements and the required positional accuracy.

4.3. Number of Receivers

The best performance is achieved when the transmitter is inside the perimeter created by the receivers. The performance is also affected by the spacing between the receivers. These two constraints will dictate the number of receivers required.

The most efficient spacing may be an equilateral triangular lattice pattern, which puts receivers close to any transmitter within its coverage area. A plot of this type of scenario and the resulting stability is shown in Figure 9. The simulations show that this lattice arrangement gives consistent coverage over the area bounded by the sensors, with less stable responses outside of that area.

4.4. Selection of Sensors

The HM algorithm requires at least 5 receivers to accurately calculate the position of a transmitter. If a system has more than 5 receivers then there are three options for dealing with the extra sensors.

1) The extra sensors can be omitted from the calculations. In this case, the sensors with the highest signal-to-noise ratio are typically the best ones to use in the calculation.

2) The extra sensors can be included in the calculations. The algorithm can easily accommodate additional measurements. If there is error in the extra measurements it could degrade the quality of the calculations.

3) The calculations can be done using weighted measurements. Weighting the measurements according to a factor such as the signal to noise ratio can be a way to include the extra information from additional sensors without harming the calculated result. This is done by weighting each row of $A_H$ and $b_H$ according to the quality of the data in that row.

5. Conclusions

Incorporating triangulation measurements into a multilateration system produces significant improvement in the accuracy of the ML algorithm. Its ability to accurately calculate height while using sensors on or near the ground makes it useful in real-world scenarios where objects are being tracked in 3 dimensions. It can make multilateration useful for tracking aircraft in flight over a wide area. This could include tracking aircraft in flight near airports or drone swarms.

Except for a few particular cases the algorithm performs well in a variety of sensor configurations. When the receivers are nearly coplanar it consistently out-performs a multilateration-only algorithm. Adapting the algorithm for receivers that do not share a common vertical direction is an area for further research.

NOTES

¹For example, in the scenarios in Section 1 the transmitters are spaced 5000 m from the center of the figure so $w\approx 5000$ produces a good balance between the multilateration and triangulation terms.

²Additionally, a Gaussian error with a standard deviation of 100 m was added to the pseudodistances used by the HM algorithm. This term did not affect the results because it was constant across all of the pseudodistances and therefore cancelled when the ${\tilde{d}}_{i\mathrm{,}j}$ terms were calculated.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Stimson, G.W. (1998) Introduction to Airborne Radar. 2nd Edition. SciTech, Mendham, 3, 471, 474. https://doi.org/10.1049/SBRA101E
[2]	Skolnik, M.I. (1962) Introduction to Radar Systems. McGraw-Hill, New York, 175-176.
[3]	Liu, Y., Jiao, Y. and Ma, H. (2020) Joint 2-D Angle Estimation Using TDOA in Distributed Multi-Antenna System. 2020 IEEE International Conference on Power, Intelligent Computing and Systems (ICPICS), Shenyang, 22 September 2020, 642-650. https://doi.org/10.1109/ICPICS50287.2020.9202370
[4]	Soliman, M.S., Morimoto, T. and Kawasaki, Z.-I. (2006) Three-Dimensional Localization System for Impulse Noise Sources Using Ultra-Wideband Digital Interferometer Technique. Journal of Electromagnetic Waves and Applications, 20, 515-530. https://doi.org/10.1163/156939306776117027
[5]	Szüllö, á. and Seller, R. (2016) Maneuvering Target Tracking in Wide Area Multilateration Radar System. 2016 17th International Radar Symposium (IRS), Krakow, 23 June 2016, 1-4. https://doi.org/10.1109/IRS.2016.7497345
[6]	Darabseh, A., Bitsikas, E., Tedongmo, B. and Pöpper, C. (2020) On ADS-B Sensor Placement for Secure Wide-Area Multilateration. Proceedings of 8th Open Sky Symposium, 59, 3. https://doi.org/10.3390/proceedings2020059003
[7]	Norrdine, A. (2015) An Algebraic Solution to the Multilateration Problem. 2012 International Conference on Indoor Positioning and Indoor Navigation, Sydney, 13-15 November 2012, 1-4. https://doi.org/10.13140/RG.2.1.1681.3602