Wide-area Multilateration Airspace Surveillance with Unsynchronized Low-Cost ADS-B Receivers Using TDOA Observations

Clemens Sonnleitner; Thomas Hobiger

doi:10.33012/navi.704

Abstract

This paper proposes and evaluates a novel approach for wide area multilateration (WAM) airspace surveillance based on time difference of arrival (TDOA) navigation. Unlike commercial Automatic Dependent Surveillance-Broadcast (ADS-B)-based WAM solutions, which require high-grade clock synchronization, the framework proposed here achieves airspace surveillance without the need for highly stable clocks or time synchronization between ground stations. In the proposed approach, aircraft 3D positions and velocities and the relative clock offsets of the receivers are estimated consistently using an extended Kalman filter (EKF). The accuracy of the 3D aircraft position estimates was tested using simulated ADS-B messages across a variety of different ground station network configurations.

Keywords

1 INTRODUCTION

Automatic Dependent Surveillance-Broadcast (ADS-B), one of the key components of the Next Generation Air Transportation System (NextGen) in the U.S. and Single European Sky ATM Research Programme (SESAR) in Europe, has become a widespread technology for airspace management and control. According to Eurocontrol (2024), the share of flights performed by aircraft equipped with an ADS-B v2-compliant transmitter reached 96.3 % in February 2024.

Compared to classic radar systems, ADS-B has higher accuracy, higher update rates, and lower costs. However, the ADS-B standard is quite old and relies on technologies from the early- to mid-20th century. It was developed without any security concerns in mind, making it vulnerable to different kinds of attacks that could have severe impacts on air traffic management safety (Schäfer et al., 2013, 2014).

Many of these disadvantages can be mitigated by wide area multilateration (WAM), which uses ADS-B transmissions and other aircraft signals for hyperbolic navigation. The WAM approach obtains precise and robust aircraft positions using time difference of arrival (TDOA) measurements between receivers that are spread over a large area (Neven et al., 2005). However, a major challenge of such systems is that local clocks, which drive the A/D converters in the ADS-B ground-stations, must be synchronized across large distances in order to calculate the transmitters’ positions.

Solutions to the mentioned problems have been suggested by e.g. Wu (2012). This paper presents a novel approach that eliminates the stringent requirement for precise clock synchronization. It shows that the relative clock offsets between the different receivers, together with the positions and velocities of the aircraft, can be estimated in a single, fully self-consistent extended Kalman filter (EKF). The paper also reports the results of simulations performed to demonstrate the feasibility of this approach, and it evaluates the accuracy of the estimated aircraft positions and velocities.

The remainder of this paper is structured as follows: In Section 2 classic WAM is introduced. Furthermore the feasibility of WAM without clock synchronization and the theory behind the novel approach of this paper, called unsynchronized wide area multilateration (uWAM), is presented. In Section 3 the foundations of the simulations that were performed to proof the applicability of uWAM are described. In Section 4 the simulation results are presented and discussed in detail. Conclusions and an outlook are presented in the final Section 5.

2 WIDE AREA MULTILATERATION

2.1 Classical WAM

WAM is a surveillance technique for air traffic control over large areas. It serves as a complementary surveillance tool to ADS-B or as a replacement for secondary radar. The basic operating principle of WAM is depicted in Figure 1. In brief, a network of ground stations receives, decodes, and time-stamps signals, such as ADS-B and transponder replies (e.g., Mode-S Orlando, 1989), that are transmitted from aircraft within a given airspace volume (EUROCAE, 2010). Multilateration based on TDOA (optionally combined with barometric altitude information from the aircraft) enables accurate position estimation (FAA, 2023; Neven et al., 2005). For simplicity, this paper focuses on ADS-B messages, but most of the methods also apply to other signals such as transponder replies.

FIGURE 1

The WAM measurement principle is based on signals transmitted by an aircraft and received at time-synchronized ground stations g_i. From these time-stamped time of arrival (TOA) measurements, the TDOA observables can be calculated. If the signal has been received by four or more geometrically well-distributed ground stations, the aircraft position can be determined using hyperbolic navigation.

Using TDOA cancels out the actual transmission time and any internal delays in the ADS-B transmitter by differencing the signal arrival times at two or more ground stations. However, these ground stations must be accurately synchronized so that the TDOA measurements can be interpreted as differences in geometric distance. A major advantage of WAM is that it uses existing aircraft transmission standards without the need for additional active or passive radar signals. The range difference Δd between the TDOA reference station g_tr and another ground station g_i for a message sent by aircraft a_k can be calculated from the difference in arrival times according to O’Keefe (2017):

$\begin{matrix} Δ d & = & c \cdot Δ t \\ = & c \cdot (Δ t_{a_{k} g_{i}} - Δ t_{a_{k} g_{t r}}) \\ = & {‖ r_{a_{k}} - r_{g_{i}} ‖}_{2} - {‖ r_{a_{k}} - r_{g_{t r}} ‖}_{2} \end{matrix}$ 1

where $r_{a_{k}}$ is the 3D transmitter position, $r_{g_{i}}$ and $r_{g_{tr}}$ are the 3D receiver positions, and $Δ t_{a_{k} g_{i}}$ and $Δ t_{a_{k} g_{t r}}$ are the signal propagation delays from the aircraft to the ground stations for signals traveling at velocity c.

One of the major shortcomings of current WAM implementations is that clock synchronization between ground stations is complicated and costly. Moreover, many commercial solutions achieve time synchronization using global navigation satellite system (GNSS), which makes the system vulnerable to jamming and spoofing (Neven et al., 2005; Kujur et al., 2020).

2.2 Feasibility of WAM Without Clock Synchronization

This paper discusses a novel approach to overcome this requirement for precise clock synchronization.

In general, relative precise clock synchronization between the ground stations is required, but individual clocks do not need to align to an absolute time reference like Coordinated Universal Time (UTC) or Global Positioning System (GPS) time. However, if the relative clock offsets of the individual stations together with all aircraft positions could be estimated in a single adjustment process, external synchronization of the station clocks would no longer be necessary.

The following sections investigate whether these unknown parameters can be estimated simultaneously, given a sufficient number of aircraft and ground stations.

2.2.1 Single Epoch Observability

To assess the observability of the system, the proposed method begins with a simple dimension check. In a scenario with N_g ground stations for which the relative clock offsets need to be estimated and N_a tracked aircraft positions, the system has

$N = 3 \cdot N_{a} + N_{g} - 1$ 2

unknown variables. Assuming that the TDOA observations are formed from one ADS-B message per aircraft per epoch, the total number of observations per epoch would be

$M = (N_{g} - 1) \cdot N {}_{a}$ 3

assuming no signal losses. The difference Δ between the number of unknowns and the number of observations equals:

$\begin{matrix} Δ & = & M - N \\ = & N_{g} N_{a} - 4 N_{a} - N_{g} + 1 \end{matrix}$ 4

If Δ<0 (i.e., more unknowns than observations), the system is underdetermined. The impact of different numbers of aircraft and ground stations on system observability, based on this equation, is visualized in Figure 2.

FIGURE 2

Visualization of Equation (4), which describes the difference Δ between the number of observations and the number of unknowns in the system. For scenarios where Δ < 0 (red, orange, and yellow dots), the system is underdetermined. The black line indicates where Δ = 0. The blue and green dots mark setups for which Δ ≥ 0.

However, an overdetermined system is not necessarily solvable in a least-square sense if the system suffers from rank deficits. The second check therefore involves determining how the number of randomly distributed receivers and aircraft affect the matrix rank. If there is a rank deficit, the observation equation cannot be inverted in a least-square estimator, so not all unknown parameters can be estimated. For example, Figure 3 shows the same setup as Figure 2 but with the rank deficit N – rank(A^TA), where A is the linearized observation operator that maps from the state space to the observation space. Both analyses indicate that, given a sufficient number of aircraft and ground stations, the aircraft positions and relative clock offsets can be estimated simultaneously. In other words, for airspace volumes that are covered by a sufficient number of ground stations and in which a sufficient number of ADS-B-equipped aircraft are operating, the requirement for external clock synchronization can be avoided. Instead, aircraft positions and relative clock offsets can be estimated in a single adjustment process.

FIGURE 3

Difference between the number of unknowns and the rank of the A-matrix for a simulation with randomly placed ground stations and aircraft. Blue dots indicate system setups that do not suffer from rank deficit.

However, this first observability analysis is based on a simplified setup. In a realistic scenario, especially in WAM, not all ground stations will receive the signal (refer to Section 3.1). and the tree spanned by the relative clock offsets of the ground stations might become fragmented and fall into multiple separated networks. This fragmentation can be monitored using graph theory to check the number of spectral clusters (von Luxburg, 2007).

Additionally, realistic scenarios must be dynamic because aircraft positions change over time. For airspace surveillance applications, the velocity vector of the aircraft is particularly important for continuously monitoring and predicting the airspace over time. Instead of solving for aircraft positions and relative clock offsets independently for each epoch using a least-squares estimator, the method proposed above can be extended by using an EKF to estimate all unknown parameters, including the aircraft 3D velocity, in a dynamic system. The EKF accounts for dynamic changes in the aircraft position and the relative clock offsets and can even provide state estimates for single epochs in which the system is underdetermined.

2.2.2 Extended Kalman Filter

The following paragraphs briefly introduce the EKF, a well-established variant of the Kalman filter (KF) for nonlinear systems. Most of the derivations and notation are based on Groves (2013), where ${\hat{x}}_{n}$ is the estimated KF state vector at epoch n and P_n is the corresponding error covariance matrix. The state vector and covariance matrix are denoted as ${\hat{x}}_{n}^{-}$ and $P_{n}^{-}$ after the prediction step and ${\hat{x}}_{n}^{+}$ and $P_{n}^{+}$ after a measurement update.

PREDICTION

Assuming linear state transition models, the state vector $\hat{x}$ can be predicted from the previous epoch n –1 to the current one n using the transition matrix $Φ_{n - 1}$ . This matrix is a function of the time interval τ_s, such that the time-discrete prediction is given as Equation (5). This equation assumes that the expected system noise contributions are all zero.

In addition to the new state prediction, the error covariance matrix P must be propagated to the new epoch using Equation (6), where Q_n–1 is the process noise covariance matrix. This matrix increases the state uncertainties over time and accounts for any errors from neglected dynamics and other random sources in the dynamic system.

The equations for the predicted state vector and covariance matrix are:

${\hat{x}}_{n}^{-} = Φ_{n - 1} {\hat{x}}_{n - 1}^{+}$ .5

$P_{n}^{-} = Φ_{n - 1} P_{n - 1}^{+} Φ_{n - 1}^{T} + Q_{n - 1} .$ 6

UPDATE

The difference between the observed measurements z_n and the expected measurements $h ({\hat{x}}_{k}^{-}, t_{n})$ , the latter of which are derived from the predicted state estimate, is called the measurement innovation:

$δ z_{n}^{-} = z_{n} - h ({\hat{x}}_{k}^{-}, t_{n}) \cdot$ 7

The measurement innovation includes both state estimation errors and measurement errors. The measurement errors are assumed to have a zero-mean Gaussian distribution with a corresponding measurement noise covariance matrix R. To map the non-linearity between the state and observation space using a linear model, this relationship is approximated using the measurement model matrix:

$H_{n} = {(\frac{\partial h (x, t_{n})}{\partial x})}_{x = {\hat{x}}_{k}}$ 8

For the update step, the Kalman gain matrix K_n can thus be calculated from the measurement model matrix using:

$K_{n} = P_{n}^{-} H_{n}^{T} \underset{S_{n}}{\underset{︸}{(H_{n} P_{n}^{-} H_{n}^{T} + R_{n})}} - 1$ 9

where S_n is the innovation covariance matrix. The Kalman gain K_n maps the innovation space to the state space corrections by balancing the effects of process and observation noise. The updated state estimate ${\hat{x}}_{n}^{+}$ and associated error covariance matrix $P_{n}^{+}$ can then be calculated using K_n and H_n in the following equation, shown below in the so-called “Joseph form”:

${\hat{x}}_{n}^{+} = {\hat{x}}_{n}^{-} + K_{n} δ z_{n}^{-}$ 10

$P_{n}^{+} = (𝟙 - K_{n} H_{n}) P_{n}^{-} {(𝟙 - K_{n} H_{n})}^{T} + K_{n} R_{n} K_{n} .$ 11

In practice, the covariance matrix P is often initialized with very large values to reflect large a priori uncertainties. However, in the case of EKF, these high initial values can also contribute to filter divergence (Groves, 2013).

2.3 Unsynchronized Wide-area Multilateration (uWAM)

The ADS-B standard is well-suited for semi-cooperative, TOA-based navigation in which the signal structure is known but the transmission times are not. Another advantage of ADS-B is that every ADS-B message encodes the International Civil Aviation Organization (ICAO) address of the aircraft, making it easy to associate the arrival time of the same message at different ground stations and form the associated TDOA observations.

Having demonstrated above that unsynchronized ADS-B ground stations can be used for WAM by estimating relative clock offsets as additional parameters (i.e., as states in a Kalman filter), this paper now explores the concept of unsynchronized wide-area multilateration (uWAM). This approach builds on the findings from Section 2.2 by using a fully consistent dynamic state estimator, in this case an EKF, to provide aircraft positions and relative clock states at any given epoch.

The KF state vector consists of the state entries for the N_a currently tracked aircraft a_k, where $k \in (1, 2, ..., N_{a})$ , and for the N_g – 1 ground stations g_i that are currently receiving ADS-B messages, with $i \in (1, 2, ..., N_{g} - 1)$ . The state vector entry for each aircraft a_k includes a 3D position vector r_k and a 3D velocity vector ${\dot{r}}_{k}$ , both in Earth-centered Earth-fixed (ECEF) coordinates. For the ground stations, only the relative clock offsets $Δ t_{c r, i}$ to an arbitrarily chosen clock reference ground station g_cr are estimated, so the reference station g_cr does not appear in the state vector.

The full state vector contains a total of $6 \cdot N_{a} + N_{g} - 1$ entries and can be expressed as

$\hat{x} = {(\overset{g_{1}}{\overset{︷}{Δ t_{c r, g_{1}}}} \overset{g_{2}}{\overset{︷}{Δ t_{c r, g_{2}}}} \dots \overset{a_{1}}{\overset{︷}{\underset{r_{1}}{\underset{︸}{\begin{matrix} x_{1} & y_{1} & z_{1} \end{matrix}}} \underset{{\dot{r}}_{1}}{\underset{︸}{\begin{matrix} {\dot{x}}_{1} & {\dot{y}}_{1} & {\dot{z}}_{1} \end{matrix}}}}} \overset{a_{1}}{\overset{︷}{\underset{r_{2}}{\underset{︸}{\begin{matrix} x_{2} & y_{2} & z_{2} \end{matrix}}} \underset{{\dot{r}}_{2}}{\underset{︸}{\begin{matrix} {\dot{x}}_{a_{2}} & {\dot{y}}_{2} & {\dot{z}}_{2} \end{matrix}}}}} \dots)}^{T} .$ 12

Up to the authors’ knowledge, a similar approach has not been proposed so far. Similar to the idea of using TDOA data for estimating relative clock offsets in WAM scenarios is the work of Leonardi (2019). But in their work the multilateration data is explicitly not used for estimating the transmitter positions, only for anomaly detection and removal of erroneous ADS-B messages.

2.3.1 Batch Processing

Typical WAM scenarios cover a large airspace with numerous aircraft each emitting multiple messages per second. It is therefore not computationally feasible to predict and update all states for each TDOA observation or each batch of TDOA observations from a single ADS-B message.

This paper proposes an alternative approach in which observations are collected over a time span τ_s and then a single prediction and update is performed at $t_{n} = t {}_{n - 1}+ τ_{s}$ using all collected observations. However, in an example where messages are collected over 1 s, the position errors for an aircraft flying at 900 km h^–1 would be up to 250 m if no correction is applied for the unsynchronized transmission times. To mitigate this error, the observations $z_{n, i}^{k}$ are predicted to time t_n and denoted as virtual observations $z_{n, i}^{k, *}$ (see Figure 5).

FIGURE 4

Kalman filter flow diagram showing the standard filter steps (in blue) alongside additional steps required for batched processing (in green).

FIGURE 5

Schematic showing how virtual observations are calculated. Observations z_n,i collected over a batch interval $(t_{n - 1}, t_{n})$ (depicted here for only a single aircraft k) are predicted to so-called virtual observations $z_{n, i}^{*}$ at epoch time t_n. Colors indicate messages with different ADS-B type code (TC).

For the simulations performed here, the batch time τ_s was set to 1 s as a compromise between near-real-time performance and the amount of data that can be processed by the filter.

All types of ADS-B messages can be used to form TDOA observations in the update step of the KF. The specific message type can be distinguished by its type code (TC), which is part of the ADS-B message. However, the transmitted information is only relevant if the digitized message content (e.g., the transmitted GNSS position (TC 9 ... 18 and 20 ... 22) or velocity (TC 19)) is used to support the state estimation. In this case, the transmitted aircraft position can be compared to the position determined from TDOA measurements, thus providing a method to monitor the integrity of the transmitted ADS-B messages. This idea is briefly discussed in Section 5.

The following paragraphs describe the derivation of the different virtual observations.

TDOA

The batch of virtual TDOA observations $z_{T, n}^{k, *}$ of aircraft a_k at epoch t_n, which was originally registered at time t_o, can be formulated as:

$z_{T, n}^{k, *} = z_{T, n}^{k} + \underset{Geometric}{\underset{︸}{f_{1} (t_{o}, t_{n})}} + \underset{Clock error (\approx 0)}{\underset{︸}{f_{2} (t_{o}, t_{n})}} .$ 13

Function f₁ accounts for the movement of aircraft a_k over the time span $t_{n} - t_{o}$ , and f₂ corrects for clock errors over the same time interval. The latter term is assumed to be zero, as the relative clocks are modeled by a random walk process (see Section 3.2). The contribution of f₁ can be calculated using the predicted KF state ${\hat{x}}_{n}^{-}$ :

$\Rightarrow f_{1} (t_{o}, t_{n}) = h ({\hat{x}}_{n}^{-}) - h ({\hat{x}}_{o}^{-}) .$ 14

Substituting this equation into the original expression for the virtual TDOA observations gives:

$z_{T, n}^{k, *} = z_{T, n}^{k} + h ({\hat{x}}_{n}^{-}) - h ({\hat{x}}_{o}^{-}) .$ 15

The TDOA measurement innovation δz^– can therefore be written as:

$\begin{matrix} δ z_{T, n}^{k, -} & = & z_{T, n}^{k, *} - h ({\hat{x}}_{n}^{-}) \\ (13), (14) \\ = & z_{T, n}^{k} + h ({\hat{x}}_{n}^{-}) - h ({\hat{x}}_{o}^{-}) - h ({\hat{x}}_{n}^{-}) \\ = & z_{T, n}^{k} - h ({\hat{x}}_{o}^{-}) . \end{matrix}$ 16

This equation indicates that measurement innovation can be determined by evaluating h using the predicted estimated state vector ${\hat{x}}_{o}^{-}$ at the observation time t_o rather than at the epoch time t_n.

Position

An airborne aircraft transmits its position every 0.5 s in the form of ADS-B messages, which can be used as additional measurements. In this case, the virtual positional observation of aircraft k can be calculated as:

$z_{P, n}^{k, *} = z_{P, t_{o}}^{k} + {\hat{\dot{r}}}_{k}^{-} (t_{n} - t_{o})$ 17

where $\overset{\cdot}{{\hat{r}}_{k}^{-}}$ is the predicted velocity from the state vector.

Velocity

This paper does not account for ADS-B velocity messages, but this information could also be used. Because the state vector contains only aircraft positions and velocities and the dynamic model does not consider higher-order positional derivatives (see Section 2.3.2), aircraft velocities are assumed to be constant during a given batch time τ_s. The corresponding virtual observation can therefore be written as:

$z_{V, n}^{k, *} = z_{V, t_{o}}^{k} .$ 18

This simplification is justified as long as the batch time τ_s is short compared to the dynamics of the surveilled aircraft. For example, this assumption holds for τ_s= 1 s in the public transportation sector. However, if the aircraft’s state space considers accelerations or higher-order terms, their predicted values should be accounted for when mapping the actual observations to virtual observations at epoch t_n.

2.3.2 Prediction

In the KF prediction step, the positions of multiple aircraft are propagated to the next epoch. Assuming that the aircraft dynamics are expressed by an integrated random-walk model and that the relative clocks can be modeled as a simple random walk process, the velocities and relative clock offsets remain constant during the prediction step, whereas the position changes linearly with velocity over time. Under these assumptions, the prediction for a single aircraft a_k would be:

$\underset{a_{k, n}^{-}}{\underset{︸}{(\begin{matrix} r_{k}^{-} \\ {\dot{r}}_{k}^{-} \end{matrix})}} = \underset{Φ_{a_{k}}}{\underset{︸}{(\begin{matrix} 𝟙_{3 \times 3} & τ_{s} \\ 0 & 𝟙_{3 \times 3} \end{matrix})}} \underset{a_{k, n - 1}^{+}}{\underset{︸}{(\begin{matrix} r_{k}^{+} \\ {\dot{r}}_{k}^{+} \end{matrix})}}$ 19

where $τ_{s} = τ_{s} \cdot 𝟙_{3 \times 3}$ , with τ_s representing the time interval over which the prediction is made. The global transition matrix Φ_n–1 is assembled by placing the matrices $Φ_{a_{k}}$ for each aircraft at their corresponding positions along the diagonal and adding entries of 1 for the relative clock offset propagation:

$Φ_{n - 1} = (\begin{matrix} 1 & 0 & \dots & 0 & 0 & \dots \\ 0 & 1 & \dots & 0 & 0 & \dots \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ \\ 0 & 0 & \dots & Φ_{a_{1}} & 0 & \dots \\ 0 & 0 & \dots & 0 & Φ_{a_{2}} & \dots \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ \end{matrix}) .$ 20

The process noise covariance matrix for the aircraft state entries $Q_{n - 1, a_{k}}$ can be calculated following Groves (2013). Assuming an integrated random walk as the underlying noise process for each aircraft, the covariance matrix can be expressed as:

$Q_{n - 1, a_{k}} = (\begin{matrix} σ_{a_{k}, \dot{x}}^{2} \frac{τ_{s}^{3}}{3} & 0 & 0 & σ_{a_{k}, \dot{x}}^{2} \frac{τ_{s}^{2}}{2} & 0 & 0 \\ 0 & σ_{a_{k}, \dot{y}}^{2} \frac{τ_{s}^{3}}{3} & 0 & 0 & σ_{a_{k}, \dot{y}}^{2} \frac{τ_{s}^{2}}{2} & 0 \\ 0 & 0 & σ_{a_{k}, \dot{z}}^{2} \frac{τ_{s}^{3}}{3} & 0 & 0 & σ_{a_{k}, \dot{z}}^{2} \frac{τ_{s}^{2}}{2} \\ σ_{a_{k}, \dot{x}}^{2} \frac{τ_{s}^{2}}{2} & 0 & 0 & σ_{a_{k}, \dot{x}}^{2} τ_{s} & 0 & 0 \\ 0 & σ_{a_{k}, \dot{y}}^{2} \frac{τ_{s}^{2}}{2} & 0 & 0 & σ_{a_{k}, \dot{y}}^{2} τ_{s} & 0 \\ 0 & 0 & σ_{a_{k}, \dot{z}}^{2} \frac{τ_{s}^{2}}{2} & 0 & 0 & σ_{a_{k}, \dot{z}}^{2} τ_{s} \end{matrix})$ 21

where $σ_{a_{k}, \dot{x}}^{2}$ , $σ_{a_{k}, \dot{y}}^{2}$ , and $σ_{a_{k}, \dot{z}}^{2}$ are the power spectral densities (PSDs) of the velocities of aircraft k in ECEF coordinates. Because the values are usually given in a local coordinate system, they must be converted to the ECEF frame. For the relative clock offsets of the ground stations, a random walk process is assumed, and the corresponding covariance matrix can be calculated as:

$Q_{n - 1, g_{i}} = σ_{g_{c r}}^{2} τ_{s} + σ_{g_{i}}^{2} τ_{s}$ 22

where $σ_{g_{c r}}^{2}$ is the power spectral density (PSD) random walk parameter of the ground station clock to which all other relative clock offsets are referenced, and $σ_{g_{i}}^{2}$ is the random walk parameter of ground station g_i.

The overall process noise covariance matrix Q_n–1 is built by placing the individual Q-matrices for the ground stations and airplanes on the diagonal:

$Q_{n - 1} = (\begin{matrix} Q_{n - 1, g_{1}} & 0 & \dots & 0 & 0 & \dots \\ 0 & Q_{n - 1, g_{2}} & \dots & 0 & 0 & \dots \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ \\ 0 & 0 & \dots & Q_{n - 1, a_{1}} & 0 & \dots \\ 0 & 0 & \dots & 0 & Q_{n - 1, a_{2}} & \dots \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ \end{matrix}) .$ 23

Using the Φ-matrix (Equation (20)) and the Q-matrix (Equation (23)), the Kalman filter prediction of the state and the P-matrix can be performed following Equations (5) and (6). The result is the predicted state estimate ${\hat{x}}_{k}^{-}$ and predicted covariance matrix estimate $P_{k}^{-}$ .

2.3.3 Update

For the calculation of measurement innovation $δ z^{-,^{*}} = z_{n}^{*} - h ({\hat{x}}_{n}^{-}, t_{n})$ , the entries h_j need to be evaluated for all observations. This paper only presents the derivations for the TDOA and position observations (the latter of which can be extracted from the ADS-B messages). Other information encoded in ADS-B messages, such as aircraft velocity, can be used in the update step; however, that is beyond the scope of this paper.

TDOA

For an ADS-B message transmitted by aircraft a_k and received by N_a receivers, a set $z_{T}^{k,^{*}}$ of N_a – 1 TDOA observations $z_{j}^{*}$ can be computed relative to an arbitrarily chosen reference receiver g_tr. A single TDOA measurement of ground station g_i is denoted as $z_{T, t r, i}^{k,^{*}}$ . For TDOA observations, h is a nonlinear function, and its entries can be calculated as:

$\begin{matrix} h_{T, t r, i}^{k} & = & (Δ t_{a_{k} g_{i}} + Δ t_{g_{c r} g_{i}}) - (Δ t_{a_{k} g_{t r}} - Δ t_{g_{c r} g_{t r}}) \\ = & Δ t_{g_{c r} g_{i}} - Δ t_{g_{c r} g_{t r}} + \frac{1}{c} ({‖ {\hat{r}}_{a_{k}} - r_{g_{i}} ‖}_{2} - {‖ {\hat{r}}_{a_{k}} - r_{g_{t r}} ‖}_{2}) \end{matrix}$ 24

$\begin{matrix} = & Δ t_{g_{c r} g_{i}} - Δ t_{g_{c r} g_{t r}} \\ + \frac{1}{c} (\sqrt{{(x_{a_{k}} - x_{g_{i}})}^{2} + {(y_{a_{k}} - y_{g_{i}})}^{2} + {(z_{a_{k}} - z_{g_{i}})}^{2}} \\ - \sqrt{{(x_{a_{k}} - x_{g_{t r}})}^{2} + {(y_{a_{k}} - y_{g_{t r}})}^{2} + {(z_{a_{k}} - z_{g_{t r}})}^{2}}) \end{matrix}$ 25

where $r_{g_{i}}$ and $r_{g_{t r}}$ are the known ADS-B receiver positions and ${\hat{r}}_{a_{k}}$ is the estimated transmitter position from the KF state.

The measurement innovation for a single TDOA entry can then be computed as:

${(δ z_{T, t r, i}^{k})}^{-} = z_{T, t r, i}^{k} - h_{T, t r, i}^{k}$ 26

in accordance with Equation (7). The observation model matrix H for the TDOA measurements has entries defined as:

$H_{T} = (\begin{matrix} \frac{\partial h_{1}}{\partial (Δ t_{g_{c r} g_{1}})} & \frac{\partial h_{1}}{\partial (Δ t_{g_{c r} g_{2}})} & \dots & \frac{\partial h_{1}}{\partial (Δ t_{g_{c r} g_{t r}})} & \dots & \frac{\partial h_{1}}{\partial x_{a_{k}}} & \frac{\partial h_{1}}{\partial y_{a_{k}}} & \frac{\partial h_{1}}{\partial z_{a_{k}}} & \frac{\partial h_{1}}{\partial {\dot{x}}_{a_{k}}} & \frac{\partial h_{1}}{\partial {\dot{y}}_{a_{k}}} & \frac{\partial h_{1}}{\partial {\dot{z}}_{a_{k}}} & \dots \\ \frac{\partial h_{2}}{\partial (Δ t_{g_{c r} g_{1}})} & \frac{\partial h_{2}}{\partial (Δ t_{g_{c r} g_{2}})} & \dots & \frac{\partial h_{2}}{\partial (Δ t_{g_{c r} g_{t r}})} & \dots & \frac{\partial h_{2}}{\partial x_{a_{k}}} & \frac{\partial h_{2}}{\partial y_{a_{k}}} & \frac{\partial h_{2}}{\partial z_{a_{k}}} & \frac{\partial h_{2}}{\partial {\dot{x}}_{a_{k}}} & \frac{\partial h_{2}}{\partial {\dot{y}}_{a_{k}}} & \frac{\partial h_{2}}{\partial {\dot{z}}_{a_{k}}} & \dots \\ \frac{\partial h_{3}}{\partial (Δ t_{g_{c r} g_{1}})} & \frac{\partial h_{3}}{\partial (Δ t_{g_{c r} g_{2}})} & \dots & \frac{\partial h_{3}}{\partial (Δ t_{g_{c r} g_{t r}})} & \dots & \frac{\partial h_{3}}{\partial x_{a_{k}}} & \frac{\partial h_{3}}{\partial y_{a_{k}}} & \frac{\partial h_{3}}{\partial z_{a_{k}}} & \frac{\partial h_{3}}{\partial {\dot{x}}_{a_{k}}} & \frac{\partial h_{3}}{\partial {\dot{y}}_{a_{k}}} & \frac{\partial h_{3}}{\partial {\dot{z}}_{a_{k}}} & \dots \\ \frac{\partial h_{4}}{\partial (Δ t_{g_{c r} g_{1}})} & \frac{\partial h 4}{\partial (Δ t_{g_{c r} g_{2}})} & \dots & \frac{\partial h_{4}}{\partial (Δ t_{g_{c r} g_{t r}})} & \dots & \frac{\partial h_{4}}{\partial x_{a_{k}}} & \frac{\partial h_{4}}{\partial y_{a_{k}}} & \frac{\partial h_{4}}{\partial z_{a_{k}}} & \frac{\partial h_{4}}{\partial {\dot{x}}_{a_{k}}} & \frac{\partial h_{4}}{\partial {\dot{y}}_{a_{k}}} & \frac{\partial h_{4}}{\partial {\dot{z}}_{a_{k}}} & \dots \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \end{matrix})$ 27

The dots indicate that the entries must be placed at the correct positions in the matrix. The partial derivatives for the ground station entries for a single TDOA observation can be written as:

$\frac{\partial h_{T, j}^{k} (g_{o}, g_{t r, o})}{\partial (Δ t_{g_{c r} g_{i}})} = {\begin{matrix} 1, & if g_{i} = g_{o} \land g_{i} \neq g_{c r} \\ - 1, & if g_{i} = g_{t r, o} \land g_{i} \neq g_{c r} \\ 0, & otherwise \end{matrix}$ 28

where g_o and g_tr,o are the ground stations involved in entry h_T,j. For the position derivations of the aircraft, the distances between the aircraft a_k and each of the ground stations g_i and g_tr can be given as:

$s_{g_{i}} = \sqrt{{(x_{a_{k}} - x_{g_{i}})}^{2} + {(y_{a_{k}} - y_{g_{i}})}^{2} + {(z_{a_{k}} - z_{g_{i}})}^{2}}$ 29

$s_{g_{t r}} = \sqrt{{(x_{a_{k}} - x_{g_{t r}})}^{2} + {(y_{a_{k}} - y_{g_{t r}})}^{2} + {(z_{a_{k}} - z_{g_{t r}})}^{2}} \cdot$ 30

Using these distances, the partial derivatives for the position components of aircraft a_k can be written as:

$\frac{\partial h_{T, j}}{\partial x_{a_{k}}} = \frac{1}{c} (\frac{x_{a_{k}} - x_{g_{i}}}{s_{g_{i}}} - \frac{x_{a_{k}} - x_{g_{t r}}}{s_{g_{t r}}})$ 31

$\frac{\partial h_{T, j}}{\partial y_{a_{k}}} = \frac{1}{c} (\frac{y_{a_{k}} - y_{g_{i}}}{s_{g_{i}}} - \frac{y_{a_{k}} - y_{g_{t r}}}{s_{g_{t r}}})$ 32

$\frac{\partial h_{T, j}}{\partial z_{a_{k}}} = \frac{1}{c} (\frac{z_{a_{k}} - z_{g_{i}}}{s_{g_{i}}} - \frac{z_{a_{k}} - z_{g_{t r}}}{s_{g_{t r}}})$ 33

$\frac{\partial h_{T, j}}{\partial {\dot{x}}_{a_{k}}}, \frac{\partial h_{T, j}}{\partial {\dot{y}}_{a_{k}}}, \frac{\partial h_{T, j}}{\partial {\dot{z}}_{a_{k}}} = 0.$ 34

The measurement noise covariance matrix $R_{T, s_{i}}$ for a set of TDOA observations s_i from aircraft a_k can be derived from the covariance matrix $\tilde{R}$ of the corresponding TOA observations. $\tilde{R}$ is a diagonal matrix, with the TOA measurement variances $σ_{TOA, g_{i}}^{2}$ of each ground station along its diagonal. Based on the propagation of uncertainty principle:

$R_{T, s_{i}} = J \tilde{R} J^{T}$ 35

where the entries of the Jacobian matrix J are given by $J_{j i} = \frac{\partial h_{T, j}}{\partial (Δ t_{a_{k} g_{j}})}$ .

Using the TDOA observation equation in the time range domain (24), the measurement noise covariance matrix can be derived. For the simplified case in which all ground stations have the same TOA measurement uncertainty $σ_{T O A, g_{1}}^{2} = σ_{T O A, g_{2}}^{2} = \dots = σ_{TOA}^{2}$ , the measurement noise covariance matrix for a set of TDOA observations is given by:

$R_{T} = σ_{TOA}^{2} (\begin{matrix} 2 & 1 & 1 & \dots \\ 1 & 2 & 1 & \dots \\ 1 & 1 & 2 & \dots \\ ⋮ & ⋮ & ⋮ & ⋱ \end{matrix})$ 36

where the size of $R_{T, s_{i}}$ is equal to the number of TDOA observations in s_i. An example derivation is given in Section A.

The overall measurement noise covariance matrix R_T for the TDOA observations is obtained by “concatenating” the individual $R_{T, s_{i}}$ matrices:

$R_{T, s_{i}} = (\begin{matrix} R_{T, s_{1}} & 0 & 0 & \dots \\ 0 & R_{T, s_{2}} & 0 & \dots \\ 0 & 0 & R_{T, s_{2}} & \dots \\ ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}) .$ 37

ADS-B Position

In addition to the TDOA observations, the GNSS positions of aircraft (which are encoded in the transmitted ADS-B messages) can also be used as observations. For this observation type, the observation model h is the position ${\hat{r}}_{k}$ of the aircraft and can be taken directly from the predicted state vector. The measurement innovation is therefore given as:

${(δ z_{P}^{k, *})}^{-} = z_{P}^{k, *} - {\hat{r}}_{k} = z_{P}^{k, *} - (\begin{matrix} x_{k} \\ y_{k} \\ z_{k} \end{matrix}) \cdot$ 38

The observation model matrix is the 3×3 identity matrix 𝟙_3×3. One approach to populating the measurement noise error matrix is to derive its entries from the decoded ADS-B message position uncertainties. In the simulations for this paper a fixed value for the accuracy of the encoded GNSS position was assumed. According to Tesi and Pleninger (2017) and summarized in Table 1, 97 % of the ADS-B messages (ICAO version 1) report a navigation accuracy category-position (NACp) level of 8 and above. Defined in the Radio Technical Commission for Aeronautics (RTCA) DO-260 standard, this relates to a 95 % accuracy bound of 93 m which corresponds to $σ_{N, E, D} \approx \frac{93 m}{2} = 46.5 m$ if one assumes the errors are Gaussian distributed. For ICAO version 2 over 97 % of the observed ADS-B messages report a NACp level of 9 and above. This relates to a 95 % horizontal accuracy bound of 30 m (in ICAO version 2 the vertical accuracy is not reflected anymore by the NACp level). For the scope of this paper, the position error of the ADS-B messages to the true aircraft position was simulated Gaussian distributed with $σ_{N, E} = 23.6 m$ and $σ_{D} = 33.4 m$ based on the insights mentioned beforehand. This leads to a simulated standard deviation of the 3D position error of $σ_{P} \approx 47.2 m$ .

View this table:

TABLE 1.

Cumulative number of ADS-B messages per ICAO version and navigation accuracy category-position (NACp) level. (Tesi & Pleninger, 2017)

2.3.4 Initialization

For all analyses in this paper, the initial values for the aircraft position were simulated with an uncertainty of $σ_{r, 0} \approx 236 m$ , which approximately corresponds to an ADS-B NACp level 5 (95 % confidence interval < 0.5 NM). The horizontal and vertical initial velocity errors were simulated with $σ_{{\dot{r}}_{h}, 0} \approx 5.2 m s^{- 1}$ and $σ_{{\dot{r}}_{v}, 0} \approx 7.8 m s^{- 1}$ , which corresponds to about twice the accuracy of navigation accuracy category-velocity (NACv) level 1. These values are conservatively based on the insights of Tesi and Pleninger (2017), who derived ADS-B message statistics for data in 2015 and 2016. Although this paper uses simulated observations, real-world implementations could take initial values directly from the decoded ADS-B messages.

The EKF relies on an accurate state vector initialization to ensure filter convergence and construct a realistic Jacobian matrix for mapping between the state and observation spaces. This is especially important for the relative clock offsets, for which no precise a priori information is available given the assumption that each station clock is not monitored or compared against the other clocks in the network. This paper therefore derives a simple scheme for obtaining initial values for the clock offsets between the ground stations.

For initialization, the reported aircraft positions are treated as true positions. The relative clock offsets can then be calculated by comparing the observed TDOA values to the values calculated from the aircraft positions and ground receiver locations.

However, it is necessary to approximate the error introduced by the uncertainty in the reported aircraft position and the fact that not all ADS-B messages are transmitted at the same epoch but rather are collected for batch processing (see Section 2.3.1).

For example, assuming a batch time of 1 s in which the ADS-B position messages from all aircraft are analyzed, and assuming a maximum aircraft velocity of 1000 km h^–1, the resulting position error due to aircraft movement within this batch time would be less than 278 m. This error, combined with the intrinsic uncertainty in the reported aircraft position (see above) due to error propagation, results in a total position uncertainty of approximately 364 m (which corresponds to approximately 1.2 μs in the time domain). In cases where some reported aircraft positions might be spoofed or unreliable, the initial relative clock station offset uncertainty could be set at a more conservative 2 μs. For the analysis described in the following sections, the relative clock states in the filter were initialized using this more conservative approach, with the uncertainty in the covariance matrix set to 2 μs.

2.3.5 Dynamic Switching of Relative Clock Offset Reference Ground Station

As aircraft enter or exit the surveilled airspace, they are dynamically added or removed from the KF state vector to maintain the observability of the KF (see Section 2.2.1). The same principle also applies to ground stations that have not recorded any observations for an extended period. In the simulations conducted for this paper, tracked aircraft and ground station entries were removed from the KF if they had not been observed for at least 10 s.

If this requirement removed the relative clock offset reference ground station g_cr , a different ground station was selected as the new relative clock offset reference ground station. The new state ${\hat{x}}_{new}$ with respect to the new ground station can be calculated from the old state by defining a transformation matrix G such that:

${\hat{x}}_{new} = G {\hat{x}}_{old} .$ 39

The following shows an example transformation matrix G for the case where the reference ground station g_cr was initially g₀ but should now be switched to g₂:

$\underset{{\hat{x}}^{*}}{\underset{︸}{(\begin{matrix} Δ t_{g_{2}, g_{1}} \\ Δ t_{g_{2}, g_{0}} \\ Δ t_{g_{2}, g_{3}} \\ Δ t_{g_{2}, g_{4}} \\ Δ t_{g_{2}, g_{5}} \\ ⋮ \\ a_{1} \\ ⋮ \end{matrix})}} = (\begin{matrix} Δ t_{g_{2}, g_{1}} \\ - Δ t_{g_{0}, g_{2}} \\ Δ t_{g_{2}, g_{3}} \\ Δ t_{g_{2}, g_{4}} \\ Δ t_{g_{2}, g_{5}} \\ ⋮ \\ a_{1} \\ ⋮ \end{matrix}) = (\begin{matrix} Δ t_{g_{0}, g_{1}} - Δ t_{g_{0}, g_{2}} \\ 0 - Δ t_{g_{0}, g_{2}} \\ Δ t_{g_{0}, g_{3}} - Δ t_{g_{0}, g_{2}} \\ Δ t_{g_{0}, g_{4}} - Δ t_{g_{0}, g_{2}} \\ Δ t_{g_{0}, g_{5}} - Δ t_{g_{0}, g_{2}} \\ ⋮ \\ a_{1} \\ ⋮ \end{matrix}) = \underset{G}{\underset{︸}{(\begin{matrix} 1 & - 1 & 0 & 0 & 0 & \dots \\ 0 & - 1 & 0 & 0 & 0 & \dots \\ 0 & - 1 & 1 & 0 & 0 & \dots \\ 0 & - 1 & 0 & 1 & 0 & \dots \\ 0 & - 1 & 0 & 0 & 1 & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & 0 & \dots & 𝟙_{6 \times 6} & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \end{matrix})}} \cdot (\begin{matrix} Δ t_{g_{0}, g_{1}} \\ Δ t_{g_{0}, g_{2}} \\ Δ t_{g_{0}, g_{3}} \\ Δ t_{g_{0}, g_{4}} \\ Δ t_{g_{0}, g_{5}} \\ ⋮ \\ a_{1} \\ ⋮ \end{matrix})$ 40

As shown above, G is a unit matrix with entries of –1 in the column that corresponded to the new g_cr in the old state vector but only in the rows corresponding to ground stations in the new state vector. This same transformation matrix G can be used to update the estimation covariance matrix. The new P_new matrix is obtained through error propagation as:

$P_{new} = G P_{old} G^{T} .$

3 SETUP

To achieve realistic simulations, the aircraft trajectories and ground receiver networks in the simulation must represent real-world conditions and account for the physical boundary conditions of signal propagation. The simulations in this paper account for the capability to receive signals over larger distances and the need for line-of-sight (LOS) visibility, but they disregard any message loss due to garbling. The following sections describe key parameters about the simulation methods implemented here.

3.1 Visibility

Receiver visibility is a critical consideration for a valid simulation. The simulations in this paper account for two major effects related to receiver visibility: signal loss due to insufficient reception power or a low signal-to-noise ratio (SNR), and signal blockage when there is no direct LOS connection between the transceiver and receiver.

3.1.1 Signal Strength

The maximum range from which aircraft ADS-B signals can be observed by receivers was set to ≈ 220 km. This value is based on the maximum distance at which aircraft are detected by the low-cost ADS-B receiver installed on roof of the institute where this research was conducted. Measurements of receivers located further than ≈ 220 km from an aircraft were neglected in the simulations. Compared to the ranges of up to 450 km reported by Kaune et al. (2012) and Petchenik (2015), the range chosen here is realistic if not conservative.

3.1.2 Signal Loss due to Missing LOS

The simulations reported here account for LOS-based signal blockage due to the horizon. However, the effects of local terrain are neglected because signal blockage due to terrain features likely plays a minor role in the ADS-B WAM scenario. Airplanes are mostly at high altitudes, ADS-B equipment is mounted thoughtfully, and the landing phase occurs near airports, which usually have good ADS-B receiver coverage.

Determining which signals remain above the horizon line of different receivers requires algorithms similar to those developed for ray casting in computer graphic rendering. The algorithm used here is derived from Ring (2013), who refer to this problem as horizon culling.

Figure 6 depicts an overview of the general approach to horizon culling; Ring (2013) provide a more in-depth explanation. In brief, receiver antennas (depicted by circles) in the gray area cannot detect the signal because their LOS to the aircraft is blocked by the horizon. The specific horizon culling algorithm described by Ring (2013) proceeds in three steps.

FIGURE 6

Overview of the horizon culling algorithm implemented here. Receiver antennas (depicted by circles) in the green areas labeled and (i.e., ground stations g₂, g₃, and g₄) will detect transmitted signals, whereas for receiver antennas in the gray area (i.e., ground station g₁), the signals are blocked by the curvature of the Earth.

Step 1: “Coordinate rescaling”

The simulations conducted here assume an elliptical earth. In the first step of horizon culling, the elliptical coordinates are transformed into a scaled unit-sphere coordinate system using:

$r^{*} = (\begin{matrix} \frac{1}{a} & 0 & 0 \\ 0 & \frac{1}{a} & 0 \\ 0 & 0 & \frac{1}{b} \end{matrix}) . r$ 41

where a = 6378137 m and b = 6356752 m are the semi-major and semi-minor axes, respectively, of the WGS84 reference ellipsoid (Ring, 2013).

Step 2: “Plane Test”

The second step of horizon culling involves checking whether each ground station receiver g_i is located in area . This area is confined by plane Π₁ which is defined by the tangent points of the cone from the aircraft to the sphere (dashed lines). In 3D space, those points form a circle. This area can be described by the inequality (Ring, 2013):

$- r_{g_{i}}^{*} r_{a}^{*} + r_{a}^{*} r_{a}^{*} \leq {(r_{a}^{*})}^{2} - 1.$ 42

For each receiver where this inequality is fulfilled, the receiver itself is on the same side of the plane as the transceiver, so there is no blockage due to curvature of the Earth. No further analyses of these receivers is necessary because it can be concluded that the aircraft signals will be received by the ground station.

Step 3: “Cone Test”

However, for receivers where Equation (42) is not fulfilled, a third step must be performed to determine whether the receiver antenna (depicted by circles) is above or below the horizon (i.e., above or below the dashed lines representing the cone in 3D). If the antenna is in area (i.e., above the horizon), it has a direct LOS to the aircraft, and the ground station will detect signals from the aircraft (e.g., g₂). However, if the antenna is in area (i.e., below the horizon), it lacks a direct LOS, and the ground station will not detect signals from the aircraft. Expressed mathematically, this step involves checking whether α ≤ β, following (Ring, 2013):

$\frac{{(- r_{g_{i}}^{*} r_{a}^{*} + r_{a}^{*} r_{a}^{*})}^{2}}{{(r_{g_{i}}^{*} - r_{a}^{*})}^{2}} \leq {(r_{a}^{*})}^{2} - 1.$ 43

If the inequality holds, the receiver is above the horizon and will receive the signals emitted by the aircraft transceiver.

3.2 Clocks

This paper uses a simple random walk model to simulate the ground station clocks. The clock parameters and behavior are based on a Skyworks SI532 low-jitter crystal oscillator priced at about 10 €, which reflects the goal of implementing a ground station network with low-cost receivers. Figure 7 shows an example of the various simulated clock biases over time for simulation T_G42-S, and Figure 8 depicts an Allan deviation (ADEV) plot for a representative clock with a slope of −½, typical of a random walk process.

FIGURE 7

An example of the simulated clock offsets over time for a clock simulated by a random walk process. After ≈30 minutes, the biases are distributed over up to ≈4 μs.

FIGURE 8

An ADEV for a single simulated clock showing the general clock behavior. The slope of −½ is typical of a random walk process.

4 SIMULATIONS AND RESULTS

Several simulation studies were designed and implemented to evaluate the effectiveness of the uWAM approach proposed above. These simulations tested two different methods for ground station distribution. In the first approach, an artificial regular grid of ground stations was used to guarantee stable coverage over the entire area of interest. In the second approach, receivers were placed in large cities based on the assumption that larger cities are most likely to have ADS-B-equipped airports. The distance between neighboring ground stations placed in a grid respectively the number of ground stations placed in cities have been chosen based on multiple practical considerations. One aspect was the computational complexity which rises with the number of ground stations as more relative ground station clock offsets have to be estimated in the KF and more TDOA observations have to be processed. Another consideration that was made relates to the visibility and the fact that the ADS-B messages have to be received by a sufficiently large number of ground stations to apply multilateration (MLAT). As a third aspect one needs to consider that implementing such a system gets more costly the more ground stations have to be deployed in the field. The chosen distances respectively the number of ground stations placed in cities in the presented simulations are a trade-off between those different aspects.

ADS-B data from the OpenSky Network (Schäfer et al., 2014) was used to ensure realistic aircraft trajectories. Trajectories were based on the first 30 minutes of the dataset “states_2022-01-31-16” (→ 2022-01-31 from 16:00 to 16:30), which can be downloaded from https://opensky-network.org/datasets/#states/2022-01-31/16/. The data was filtered to remove obvious outliers and altitudes below 600 m. The ADS-B messages and their transmit times were then simulated based on these trajectories. For simplicity, only ADS-B aircraft identification and airborne position messages were simulated. The corresponding message rates were 0.2 Hz and 2 Hz, respectively, following the ADS-B version 2 standard (Sun, 2021). Both message types were used to derive TDOA observations, and the airborne position messages were used in some simulations to provide additional direct position observations.

View this table:

TABLE 2

This paper investigates several simulation scenarios defined by the parameters listed here. Parameters common across all scenarios are described in Section 3 and Section 4. The simulations discussed in the main text are highlighted, while the others are available online as supplemental material. Simulations are assigned a unique identifier following the naming convention shown, in which the individual parts of the identifier describe the most relevant parameters that differ between the simulations.

An overview of the tested simulations is given in Table 2. Each simulation was assigned a unique identifier containing the relevant simulation parameters following the naming convention next to Table 2.

The ground station placements and simulated flight trajectories for the simulations are depicted in Figure 9. Plots of the dilution of precision (DOP) values for the simulated placing types, derived from the observation model matrix H (Bard & Ham, 1999), are provided in the supplemental material (Auxiliary plots.pdf) which is available online and additionally on https://q-wertz.github.io/2024-04_navigation_auxiliary. One can see, that for the grid placed ground stations most of the area of interest has a horizontal DOP below five. The optimization of ground station locations is a non-trivial task and is out of the scope of this paper. In general one could address the problem of ground station distribution by simply spreading them equally over the area. This rule of thumb is justified by the fact that ground stations that are located close together receive messages almost with the same travel time, which in turn do not carry new information in terms of hyperbolic navigation and thus do not improve the geometry. Naghavi and Norouzi (2017) and Yang and Scheuing (2005) provide more elaborated solutions of such an optimization problem. From the DOP plots one can conclude, that the accuracy in the vertical direction is worse than in the horizontal direction. This originates from the fact that ground stations are only spread over the horizontal domain and their distribution in the vertical direction is limited.

FIGURE 9

Simulations included different ground station placements following two different placement strategies. The first strategy simulated uniform coverage, such as the grid with 120 km spacing shown here. The second strategy placed ground stations in large cities; the 40 and 60 ground stations placed in the largest cities are shown here. Ground station placements are shown in Mercator projection. Simulated aircraft trajectories were consistent across all scenarios and are therefore only shown once. Plots of the DOP values are provided in the online supplemental material.

Most simulation runs assumed a standard deviation of $100 ns \hat{\approx} 30 m$ for the uncertainty of the TOA of an ADS-B message at any given ground station (see Table 2). This number is a conservative estimate based on the findings from Steffes et al. (2011) and is supported by a zero-baseline test conducted by Chen et al. (2013) on their TDOA equipment, which yielded a standard deviation of approximately $8.8 ns \hat{\approx} 2.6 m$ . However, a test with a larger measurement noise of 1 μs (simulation T_G42-S-1u) was also performed to demonstrate the effect of this stochastic property on parameter estimation accuracy.

Figure 10 provides an overview of the 3D position error statistics across the different simulations. Several key observations emerge from this analysis. First, for all simulations except T_G42-S-1u, for which the TOA measurement noise was 1 μs, the average position error was well below 100 m, with a median error less than 50 m and a third quartile (Q₃) below 100 m. These low errors demonstrate the feasibility of the uWAM approach while also indicating that performance scales with TOA measurement noise. Furthermore, as expected, a higher number of more well-distributed ground stations reduces average errors and limits outliers. Third, comparing T_G42-S to T_G42-P and T_C60-S to T_C60-P shows that having perfect clocks does not significantly improve overall performance. This finding further suggests that the uWAM approach proposed here can achieve a level of performance comparable to commercial solutions that rely on accurately synchronized clocks. However, the fact that the mean error exceeds the median in every simulation indicates the presence of a small number of large outliers that degrade the overall performance. These outliers are discussed further in the following section comparing T_C40-S to TP_C40-S, the latter of which also includes ADS-B position observations (Figures 15 and 16). Finally, accounting for position observations (i.e., comparing simulation T_G42-S to TP_G42-S and T_C40-S to TP_C40-S) reduces the magnitude of the errors and especially the number of outliers, an effect that is discussed in Section 4.2 below for simulations T_C40-S and TP_C40-S.

FIGURE 10

Box-and-whisker plot showing the 3D position errors in the different simulation scenarios, with an additional dashed line indicating the mean error. Numbers list the median error, and whiskers indicate the last value within 1.5 times the interquartile range (IQR) (where IQR = Q₃ – Q₁), measured from the first (Q₁) and third (Q₃) quartile. Outliers are shown as dots further beyond the whiskers. In some cases, these outliers appear as thick lines due to grouping and the large amount of data.

The following sections discuss specific simulations to demonstrate the performance and applicability of the proposed algorithm. The weighted root-mean-square (WRMSE) errors are also depicted for these simulations, as the KF provides the variances of the state estimates and therefore provide an indication of their reliability. In particular, the DOP is contained within the KF update, so using WRMSE instead of unweighted errors reduces the effect of ground station placement. More figures, especially for the simulations that are not discussed here, are available in the online supplemental material (Auxiliary plots.pdf) and at https://q-wertz.github.io/2024-04_navigation_auxiliary.

4.1 Grid-Placed Ground Stations: T_G42-S

In this simulation, the receivers are placed in a grid with an edge length of approximately 120 km to provide uniform coverage. The aircraft trajectories and ground station positions are shown in Figure 9, and the WRMSE errors are shown in Figure 11. Outside of the KF initialization phase, the WRMSE of the horizontal position estimates is below 30 m for most parts of the simulated period, and the WRMSE for altitude is below 50 m. Figure 12 shows that > 90% of the errors for the horizontal components are below 20 m, a promising result that further supports the uWAM approach. In addition, the symmetry around zero indicates no systematic biases in the solutions. From the figures which depict the position errors, the sensitivity of MLAT to transmitter and receiver placement is evident. All setups in which the placement is not uniform lack accuracy in the respective direction (in our case the vertical).

The WRMSE for most of the relative clock offsets, also depicted in Figure 11, is between 0.2…0.7×10^–7 s, which matches the simulated noise of the clocks.

FIGURE 11

The WRMSEs of the estimated parameters in simulation T_G42-S, calculated based on the variances estimated by the EKF. The low errors demonstrate the applicability of the proposed uWAM algorithm: most horizontal position errors are below 30 m, and the vertical position errors are below 50 m.

FIGURE 12

Histogram of the position errors of simulation T_G42-S. In this simulation, ground stations are placed on a regular grid, and only TDOA observations are used.

FIGURE 13

Histogram of the horizontal and vertical position errors in simulation T_C40-S (which used only TDOA observations). Large numbers of outliers can be seen, especially in the altitude direction.

FIGURE 14

Histogram of the horizontal and vertical position errors in simulation TP_C40- S (which used both TDOA and position observations). Compared to simulation T_C40-S, the inclusion of ADS-B position observations reduces both the magnitude of the position errors and the number of outliers.

4.2 Irregular Placed Ground Stations: T_C40-S and TP_C40-S

Several simulations investigated irregularly placed ground stations. Two such sim-ulations (T_C40-S and TP_C40-S), in which receivers were located in the 40 largest cities in the surveilled area, are presented here as a comparison to simulations that contained approximately the same number of receivers distributed in an even grid. Figure 9 shows that some areas in the western part of Germany have increased receiver density while other areas have almost no coverage. In simulations that used only TDOA observations (e.g., T_C40-S), aircraft could not be tracked in some regions because their signals were not received by enough ground stations to perform WAM.

Comparing T_C40-S to T_G42-S shows that the WRMSE in the horizontal direc-tions is more or less the same regardless of the grid placement strategy, but the standard deviation of the WRMSE for altitude is about twice as large when receiv-ers are distributed irregularly in large cities.

As soon as the ADS-B position observations are taken into account (simulation TP_C40-S) the errors shrink significantly. This can be seen especially when comparing the absolute position errors (Figures 15 and 16) and the WRMSE errors of all estimated parameters (Figures 17 and 18). In particular the position accuracy in the vertical direction benefits from the absolute position observations, as MLAT in general suffers from the fact that ground stations are usually only located at very low altitudes, which is not in favor of an equally sampling in the vertical domain. From the comparison of the position errors one can also see that the number of extreme outliers is reduced dramatically when considering absolute position observations.

FIGURE 15

Plots of the individual position errors over time for simulation T_C40-S (which only uses TDOA observations) show that the worst position results consistently originate from a select few aircraft. The altitude error is particularly affected by bad DOP values due to the non-optimal ground station distribution.

FIGURE 16

Plots of the individual position errors over time for simulation TP_C40-S (which uses both TDOA and position observations) show the improved stability caused by the additional use of ADS-B positions. In particular, the altitude outliers are significantly reduced compared to Figure 15.

FIGURE 17

The WRMSE over time for simulation T_C40-S, in which ground stations were not optimally distributed, shows that the position errors are generally noisier and, on average, worse compared to simulations with uniformly distributed ground stations.

FIGURE 18

The WRMSE over time for simulation TP_C40-S shows the improved accuracy that results from including the ADS-B position observations.

Despite this variation in estimation accuracy, the relative clock offset errors are almost the same across all three simulations. This consistency is likely due to the simulated error of the position observations (approximately 58 m; Section 2.3.3) being in the same range as the TOA measurement errors $(100 n s \tilde{\approx} 30 m)$ .

4.3 On the Correspondence Between Formal Errors and Position Errors

The previous subsections focused on position errors (i.e., the difference between the observed states and the simulated states). However, because the KF provides not only state estimates but also the associated covariance estimates, it is possible to compare this stochastic information to the position errors. Using the estimated covariance matrix of each aircraft’s position estimate, the 3D formal position error $σ_{r_{k}}$ of an aircraft’s position estimate can be calculated as:

$σ_{r_{k}} = \sqrt{tr (P_{n}^{r_{k}})} .$ 44

FIGURE 19

Visualizing the 3D error against the formal errors estimated by the KF in simulation T_C40-S shows that the 3D errors correspond well with the estimated formal errors.

FIGURE 20

The plot of the 3D error against the estimated formal errors for simulation TP_C40-S show that the formal errors underestimate the actual errors. The errors associated with the initialization of the aircraft are visible as vertical distributed points.

This calculation allows the 3D position errors to be plotted against their corresponding formal errors, as depicted in Figure 19 and Figure 20. Three main conclusions can be drawn from these figures.

First, as discussed in the previous section, including position observations (i.e., the aircraft positions reported in the ADS-B messages) in addition to TDOA observations leads to better position estimates. This improvement is evidenced by the smaller 3D position errors in Figure 20 compared to Figure 19.

Second, when only TDOA observations are used (Figure 19), the formal 3D errors closely match the position errors. This agreement between the errors suggests that this approach could potentially be used to detect GNSS spoofing attacks, which would be visible as differences between the reported aircraft position and the position estimated by uWAM.

Third, when additional positional observations are considered (e.g., Figure 20), the formal errors become overly optimistic, underestimating the actual error by at least a factor of two. Because the KF was tuned with the same stochastic information that was used to generate random errors for the observations, the formal errors and actual errors would be expected to correspond more closely. This mismatch can instead be explained by the batch processing mode, which assumes that any observation made within a defined interval can be linearly propagated the next update epoch. However, this assumption only holds for aircraft traveling in straight lines between two consecutive update epochs. This assumption is often violated, especially by aircraft flying at lower altitudes where more maneuvers are typically performed. Formal errors smaller than 10 m therefore might not be meaningful given the assumptions made for the batch processing approach.

5 CONCLUSIONS AND FUTURE WORK

This paper demonstrates the general applicability of a modified WAM scheme that is based on TOA measurements and does not require external clock synchronization. Especially when taking into account the data that is encoded in the ADS-B messages (simulated with errors that have been reasoned in Section 2.3.3 paragraph “ADS-B position”) the position errors are well below 100 m. If a suitable ground station network (in terms of density and geometry) is available, errors are smaller than 50 m. This indicates the huge potential of the presented approach. Optimization of the placement of ground stations can further improve the accuracy as discussed by e.g. Naghavi and Norouzi (2017) and Yang and Scheuing (2005). Moreover, the usage of position information improves the filter convergance and accuracy of estimated aircraft positions significantly. For example in the simulations with 42 ground stations placed on a grid the median error was reduced from 43 m to 24 m (see Figure 10) when the encoded ADS-B positions were simulated with an accuracy of σ_P ≈ 42.7 m. Assuming sufficient coverage of ground stations equipped with ADS-B receivers, aircraft positions can be estimated using TDOA measurements and compared to the GNSS positions reported in each ADS-B message. Because each estimated position includes an associated uncertainty derived from the state covariance matrix, the reported and estimated positions can be statistically compared, thereby providing an integrity monitoring system that can flag or report jammed or spoofed aircraft positions. Given the increasing number of reported jamming or spoofing events (Kujur et al., 2020), the approach proposed here supports the development of a low-cost network capable of identifying and removing unreliable aircraft position reports.

Our future work will accordingly focus on implementing an algorithm for detecting jamming and spoofing events as well as identifying and removing outliers more generally. First results of an implemented outlier detection hypothesis test, based on the normalized innovation squared (NIS) metric (Bar-Shalom et al., 2002) were presented on the European Navigation Conference 2024 (ENC2024) (Sonnleitner & Hobiger, 2024). This can be achieved by checking whether the NIS $ϵ_{z_{n}} = (δ z_{n}^{-}) S^{- 1} {(δ z_{n}^{-})}^{T}$ is central Chi-square distributed. Outliers, such as spoofed ADS-B messages contradict the assumption that the innovations have a zero mean and that their magnitudes are consistent with the innovation covariance. First test reveal, that spoofed ADS-B messages can be reliably detected and removed by this approach.

Future work must also address the algorithm’s stability and performance, particularly for epochs with only a small number of observable aircraft in the surveilled airspace and under conditions of message loss due to signal garbling. For example, Schäfer et al. (2014) found that message loss due to signal garbling can reach up to 50 % at receivers in highly congested airspaces with over 60 aircraft. Furthermore, to make the proposed method competitive with existing alternatives, more research will be needed to address proper filter initialization and the general problem of estimating state vectors with several thousands of states. Studies could also investigate whether including estimates of the relative clock drift enhances the results. Finally, because relative clock offsets are estimated as a by-product of the proposed algorithm, future studies could explore the use of this method for time synchronization or time transfer.

A well-known problem of ADS-B in crowded airspaces is the massive spectrum load of the unsynchronized transmitters, which results in overlapping signals and higher loss rates. (Strohmeier et al., 2014; Verbraak et al., 2017) While not been taken into account in our simulations so far, based on some simple assumptions such as the arrival time and signal power some insights regarding this topic could be gathered from improved simulation runs. For the proposed uWAM algorithm it is also beneficial, that it is sufficient to reconstruct only part of the ADS-B message, as the message content correctness is of inferior importance.

All results presented here were obtained from self-consistent simulations so that the estimated parameters can be compared to the simulated truth. No real-world data has yet been incorporated, so a necessary next step will be testing this algorithm in a real-world scenario. This testing will ideally require real-time access to data streams from a ground station network in which ADS-B receivers are not synchronized using GNSS receivers. Although the assumptions in the simulations, like the maximum signal range of 220 km (Section 3.1.1) or the measurement error (Section 4), were based on the most up-to-date knowledge of the individual components, testing the algorithm on real-world data would provide better insights into potential pitfalls not captured in the simulations.

HOW TO CITE THIS ARTICLE

Sonnleitner, C., & Hobiger, T. (2025). Wide-Area multilateration airspace surveillance with unsynchronized low-cost ADS-B receivers using time difference of arrival observations. NAVIGATION, 72(3). https://doi.org/10.33012/navi.704

CONFLICT OF INTEREST

The authors declare no potential conflicts of interest.

ACKNOWLEDGMENTS

The authors want to thank the OpenSky Network for providing open access to their air traffic dataset. We further acknowledge support from German Federal Ministry for Economic Affairs and Climate Action (BMWK) project number 20V1904B (CNS-ALPHA).

Appendix

EXAMPLE DERIVATION OF THE J AND R-MATRIX

This appendix provides an example of how the measurement noise covariance matrix of the TDOA observations R_T is derived from the measurement noise of the TOA observations. The example is based on a network of seven ground stations (g₀,⋯,g₆) observing two aircraft (a₀,a₁). Aircraft a₀ has been observed by $(g_{0,} g_{1,} g_{2,} g_{3,} g_{5})$ , and the TDOA reference station g_tr for a₀ is g₅. Aircraft a₁ has been observed by $(g_{0,} g_{1,} g_{4,} g_{2,} g_{6,} g_{3})$ , and the TDOA reference station g_tr for a₁ is g₃. The relative clock offset reference station g_cr is g₆.

For clarity, the observations of the first aircraft are sorted numerically by ground station (excluding stations g₄ and g₆, which did not observe this aircraft). The observations of the second aircraft have been deliberately shuffled. For this example, the clock reference station also makes an observation to demonstrate that the J-matrix must be constructed carefully based on the structure of the H-matrix:

$h = (\begin{matrix} h_{a_{0}} \\ h_{a_{1}} \end{matrix}) = (\begin{matrix} h_{0} \\ h_{1} \\ h_{2} \\ h_{3} \\ h_{4} \\ h_{5} \\ h_{6} \\ h_{7} \\ h_{8} \end{matrix}) = (\begin{matrix} Δ t_{a_{0} g_{0}} - Δ t_{a_{0} g_{t r}} + Δ t_{g_{c r} g_{0}} - Δ t_{g_{c r} g_{t r}} \\ Δ t_{a_{0} g 1} - Δ t_{a_{0} g_{t r}} + Δ t_{g_{c r} g_{1}} - Δ t_{g_{c r} g_{t r}} \\ Δ t_{a_{0} g_{2}} - Δ t_{a_{0} g_{t r}} + Δ t_{g_{c r} g_{2}} - Δ t_{g_{c r} g_{t r}} \\ Δ t_{a_{0} g_{3}} - Δ t_{a_{0} g_{t r}} + Δ t_{g_{c r} g_{3}} - Δ t_{g_{c r} g_{t r}} \\ Δ t_{a_{1} g_{0}} - Δ t_{a_{1} g_{t r}} + Δ t_{g_{c r} g_{0}} - Δ t_{g_{c r} g_{t r}} \\ Δ t_{a_{1} g_{1}} - Δ t_{a_{1} g_{t r}} + Δ t_{g_{c r} g_{1}} - Δ t_{g_{c r} g_{t r}} \\ Δ t_{a_{1} g_{4}} - Δ t_{a_{1} g_{t r}} + Δ t_{g_{c r} g_{4}} - Δ t_{g_{c r} g_{t r}} \\ Δ t_{a_{1} g_{2}} - Δ t_{a_{1} g_{t r}} + Δ t_{g_{c r} g_{2}} - Δ t_{g_{c r} g_{t r}} \\ Δ t_{a_{1} g_{6}} - Δ t_{a_{1} g_{t r}} + Δ t_{g_{c r} g_{6}} - Δ t_{g_{c r} g_{t r}} \end{matrix}) \underline{\underline{\land}} (\begin{matrix} Δ t_{a_{0} g_{0}} - Δ t_{a_{0} g_{5}} + Δ t_{g_{6} g_{0}} - Δ t_{g_{6} g_{5}} \\ Δ t_{a_{0} g_{1}} - Δ t_{a_{0} g_{5}} + Δ t_{g_{6} g_{1}} - Δ t_{g_{6} g_{5}} \\ Δ t_{a_{0} g_{2}} - Δ t_{a_{0} g_{5}} + Δ t_{g_{6} g_{2}} - Δ t_{g_{6} g_{5}} \\ Δ t_{a_{0} g_{3}} - Δ t_{a_{0} g_{5}} + Δ t_{g_{6} g_{3}} - Δ t_{g_{6} g_{5}} \\ Δ t_{a_{1} g_{0}} - Δ t_{a_{1} g_{3}} + Δ t_{g_{6} g_{0}} - Δ t_{g_{6} g_{3}} \\ Δ t_{a_{1} g_{1}} - Δ t_{a_{1} g_{3}} + Δ t_{g_{6} g_{1}} - Δ t_{g_{6} g_{3}} \\ Δ t_{a_{1} g_{4}} - Δ t_{a_{1} g_{3}} + Δ t_{g_{6} g_{4}} - Δ t_{g_{6} g_{3}} \\ Δ t_{a_{1} g_{2}} - Δ t_{a_{1} g_{3}} + Δ t_{g_{6} g_{2}} - Δ t_{g_{6} g_{3}} \\ Δ t_{a_{1} g_{6}} - Δ t_{a_{1} g_{3}} + 0 - Δ t_{g_{6} g_{3}} \end{matrix})$

The matrix J, defined as $J_{j i} = \frac{\partial h_{T, j}}{\partial (Δ t_{a_{k} g_{j}})}$ , can be divided into individual block matrices per aircraft, as all other entries are 0:

$J = (\begin{matrix} J_{a_{0}} & 0 \\ 0 & J_{a_{1}} \end{matrix}) \to J_{a_{0}} = (\begin{matrix} 1 & 0 & 0 & 0 & 0 & - 1 & 0 \\ 0 & 1 & 0 & 0 & 0 & - 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & - 1 & 0 \end{matrix}) a n d J_{a_{1}} = (\begin{matrix} 1 & 0 & 0 & - 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & - 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & - 1 & 0 & 0 & 1 \end{matrix}) .$

For the simplified case in which all ground stations have the same observation precision $(\to \tilde{R} = σ_{TOA}^{2} \cdot 𝟙)$ , the measurement noise covariance matrix can be further simplified using $R_{T, a_{i}} = J_{a_{i}} \cdot \tilde{R} \cdot J_{a_{i}}^{T} = σ_{TOA}^{2} \cdot J_{a_{i}} \cdot J_{a_{i}}^{T}$ :

$\begin{matrix} \to R_{T, a_{0}} = σ_{TOA}^{2} \cdot J_{a_{0}} \cdot J_{a_{0}}^{T} = σ_{TOA}^{2} \cdot (\begin{matrix} 2 & 1 & 1 & 1 \\ 1 & 2 & 1 & 1 \\ 1 & 1 & 2 & 1 \\ 1 & 1 & 1 & 2 \end{matrix}) and \\ R_{T, a_{1}} = σ_{TOA}^{2} \cdot J_{a_{1}} \cdot J_{a_{1}}^{T} = σ_{TOA}^{2} \cdot (\begin{matrix} 2 & 1 & 1 & 1 & 1 \\ 1 & 2 & 1 & 1 & 1 \\ 1 & 1 & 2 & 1 & 1 \\ 1 & 1 & 1 & 2 & 1 \\ 1 & 1 & 1 & 1 & 2 \end{matrix}) \end{matrix}$

The individual $R_{T, a_{i}}$ -matrices are thus filled with values of 1, with entries of 2 on the diagonal. These individual matrices are then assembled into the TDOA measurement noise covariance matrix R_T:

$\Rightarrow R_{T} = σ_{TOA}^{2} \cdot (\begin{matrix} 2 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 2 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 2 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 & 2 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 & 1 & 2 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 & 2 & 1 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 2 \end{matrix}) \cdot$

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

↵
Bard, J., & Ham, F. (1999). Time difference of arrival dilution of precision and applications. IEEE Transactions on Signal Processing, 47(2), 521–523. https://doi.org/10.1109/78.740135
Google Scholar
↵
Bar-Shalom, Y., Li, X.-R., & Kirubarajan, T. (2002). Estimation with applications to tracking and navigation: Theory, algorithms and software (1st ed.). Wiley. https://doi.org/10.1002/0471221279
Google Scholar
↵
Chen, Y.-H., Lo, S., Akos, D. M., Wong, G., & Enge, P. (2013). A testbed for studying Automatic Dependent Surveillance–Broadcast (ADS-B) based range and positioning performance to support alternative position navigation and timing (APNT). Proc. of the 26th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2013), Nashville, TN, 263–273. https://www.ion.org/publications/abstract.cfm?articleID=11258
Google Scholar
↵
EUROCAE. (2010). Technical specification for wide area multilateration (WAM) systems. https://standards.globalspec.com/std/1288603/eurocae-ed-142
Google Scholar
↵
Eurocontrol. (2024). Automatic Dependent Surveillance -Broadcast airborne equipage monitoring. https://www.eurocontrol.int/service/adsb-equipage
Google Scholar
↵
FAA. (2023). ADS-B wide area multilateration (WAM). Federal Aviation Administration. https://www.faa.gov/air_traffic/technology/adsb/atc/wam
Google Scholar
↵
Groves, P. D. (2013). Principles of GNSS, inertial, and multisensor integrated navigation systems (2nd ed). Artech House. https://us.artechhouse.com/Principles-of-GNSS-Inertial-and-Multisensor-Integrated-Navigation-Systems-Second-Edition-P2046.aspx
Google Scholar
↵
Kaune, R., Steffes, C., Rau, S., Konle, W., & Pagel, J. (2012). Wide area multilateration using ADS-B transponder signals. Proc. of the 2012 15th International Conference on Information Fusion, Singapore, 727–734. https://ieeexplore.ieee.org/document/6289874
Google Scholar
↵
Kujur, B., Khanafseh, S., & Pervan, B. (2020). Detecting GNSS spoofing of ADS-B equipped aircraft using INS. Proc. of the 2020 IEEE/ION Position, Location and Navigation Symposium (PLANS), Portland, OR, 548–554. https://doi.org/10.1109/PLANS46316.2020.9109966
Google Scholar
↵
Leonardi, M. (2019). ADS-B anomalies and intrusions detection by sensor clocks tracking. IEEE Transactions on Aerospace and Electronic Systems, 55(5), 2370–2381. https://doi.org/10.1109/TAES.2018.2886616
Google Scholar
↵
Naghavi, S. H., & Norouzi, Y. (2017). TDOA sensor array optimization using digital elevation model. Proc. of the 2017 Iranian Conference on Electrical Engineering (ICEE), Tehran, Iran, 1757–1762. https://doi.org/10.1109/IranianCEE.2017.7985335
Google Scholar
↵
Neven, W. H. L., Quilter, T. J., Weedon, R., & Hogendoorn, R. A. (2005). Wide area multilateration report. https://www.eurocontrol.int/sites/default/files/2019-05/surveilllance-report-wide-area-multilateration-200508.pdf
Google Scholar
↵
O’Keefe, B. (2017). Finding location with time of arrival and time difference of arrival techniques (Tech Notes). https://sites.tufts.edu/eeseniordesignhandbook/files/2017/05/FireBrick_OKeefe_F1.pdf
Google Scholar
↵
Orlando, V. A. (1989). The Mode S beacon radar system. The Lincoln Laboratory Journal, 2(3), 345–362. https://www.ll.mit.edu/sites/default/files/publication/doc/mode-s-beacon-radar-system-orlando-ja-6373.pdf
Google Scholar
↵
Petchenik, I. (2015). How we track flights with ADS-B. Flightradar24 Blog. https://www.flightradar24.com/blog/how-we-track-flights-with-ads-b/
Google Scholar
↵
Ring, K. (2013). Horizon culling. Cesium. https://cesium.com/blog/2013/04/25/horizon-culling/
Google Scholar
↵
Schäfer, M., Lenders, V., & Martinovic, I. (2013). Experimental analysis of attacks on next generation air traffic communication. In M. Jacobson, M. Locasto, P. Mohassel, & R. Safavi- Naini (Eds.), Applied Cryptography and Network Security, 253–271. Springer. https://doi.org/10.1007/978-3-642-38980-1_16
Google Scholar
↵
Schäfer, M., Strohmeier, M., Lenders, V., Martinovic, I., & Wilhelm, M. (2014). Bringing up OpenSky: A large-scale ADS-B sensor network for research. Proc. of the IPSN-14 13th International Symposium on Information Processing in Sensor Networks, Berlin, Germany, 83–94. https://doi.org/10.1109/IPSN.2014.6846743
Google Scholar
↵
Sonnleitner, C., & Hobiger, T. (2024). Resilient ADS-B based airspace surveillance by means of TDOA. [European Navigation Conference 2024]. ESA ESTEC, Noorwijk, Netherlands. https://enc-series.org/2024/wp-content/uploads/2025/02/ProgramSCHEDULE_ENC2024_A5_13MEI_2.pdf
Google Scholar
↵
Steffes, C., Kaune, R., & Rau, S. (2011). Determining times of arrival of transponder signals in a sensor network using GPS time synchronization. In INFORMATIK 2011 - Informatik schafft Communities, 481–481. Gesellschaft für Informatik e.V. https://www.user.tu-berlin.de/komm/CD/paper/100123.pdf
Google Scholar
↵
Strohmeier, M., Schäfer, M., Lenders, V., & Martinovic, I. (2014). Realities and challenges of NextGen air traffic management: The case of ADS-B. IEEE Communications Magazine, 52(5), 111–118. https://doi.org/10.1109/MCOM.2014.6815901
Google Scholar
↵
Sun, J. (2021). The 1090 megahertz riddle. TU Delft OPEN Books. https://doi.org/10.34641/mg.11
Google Scholar
↵
Tesi, S., & Pleninger, S. (2017). Analysis of quality indicators in ADS-B messages. MAD - Magazine of Aviation Development, 5(3), 6–12. https://doi.org/10.14311/MAD.2017.03.01
Google Scholar
↵
Verbraak, T. L., Ellerbroek, J., Sun, J., & Hoekstra, J. M. (2017). Large-scale ADS-B data and signal quality analysis. Proc. of the 12th Seminar Papers: 12th USA/Europe Air Traffic Management Research and Development Seminar, Seattle, WA. https://drive.google.com/file/d/1G0c2TdGNAOjL0isrYL2KojU3ixQ_7FvU/view
Google Scholar
↵
von Luxburg, U. (2007). A tutorial on spectral clustering. Statistics and Computing, 17(4), 395–416. https://doi.org/10.1007/s11222-007-9033-z
Google Scholar CrossRef
↵
Wu, R. H. (2012). TCAS-aided multilateration for terminal surveillance with improved accuracy. Proc. of the 2012 IEEE/AIAA 31st Digital Avionics Systems Conference (DASC), Williamsburg, VA, 5B6-1–5B6-12. https://doi.org/10.1109/DASC.2012.6382364
Google Scholar
↵
Yang, B., & Scheuing, J. (2005). Cramer-Rao bound and optimum sensor array for source localization from time differences of arrival. Proc. of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP ’05), Philadelphia, PA, iv/961–iv/964. https://doi.org/10.1109/ICASSP.2005.1416170
Google Scholar

PDF

Previous Next

Wide-area Multilateration Airspace Surveillance with Unsynchronized Low-Cost ADS-B Receivers Using TDOA Observations

Abstract

1 INTRODUCTION

2 WIDE AREA MULTILATERATION

2.1 Classical WAM

2.2 Feasibility of WAM Without Clock Synchronization

2.2.1 Single Epoch Observability

2.2.2 Extended Kalman Filter

PREDICTION

UPDATE

2.3 Unsynchronized Wide-area Multilateration (uWAM)

2.3.1 Batch Processing

TDOA

Position

Velocity

2.3.2 Prediction

2.3.3 Update

TDOA

ADS-B Position

2.3.4 Initialization

2.3.5 Dynamic Switching of Relative Clock Offset Reference Ground Station

3 SETUP

3.1 Visibility

3.1.1 Signal Strength

3.1.2 Signal Loss due to Missing LOS

Step 1: “Coordinate rescaling”

Step 2: “Plane Test”

Step 3: “Cone Test”

3.2 Clocks

4 SIMULATIONS AND RESULTS

4.1 Grid-Placed Ground Stations: T_G42-S

4.2 Irregular Placed Ground Stations: T_C40-S and TP_C40-S

4.3 On the Correspondence Between Formal Errors and Position Errors

5 CONCLUSIONS AND FUTURE WORK

HOW TO CITE THIS ARTICLE

CONFLICT OF INTEREST

ACKNOWLEDGMENTS

Appendix

EXAMPLE DERIVATION OF THE J AND R-MATRIX

REFERENCES

In this Issue

NAVIGATION: Journal of the Institute of Navigation

Alerts for this Article

Email This Article

Citation Tools

Wide-area Multilateration Airspace Surveillance with Unsynchronized Low-Cost ADS-B Receivers Using TDOA Observations

Citation Manager Formats

Share This Article

Wide-area Multilateration Airspace Surveillance with Unsynchronized Low-Cost ADS-B Receivers Using TDOA Observations

SAVE ITEM IN SELECTED FOLDER.

Jump to section

Related Articles

Cited By...

More in this TOC Section

Similar Articles

Keywords