Algorithm for Geodetic Positioning

Based on Angle-of-Arrival of
Automatic Dependent Surveillance-Broadcasts
Richard A. Gross Nicholas A. Baine, Ph.D., P.E.
School of Engineering School of Engineering
Grand Valley State University Grand Valley State University
Grand Rapids MI, USA Grand Rapids MI, USA
grossri@[Link] bainen@[Link]

Abstract—This paper develops a non-precision geolocation I. INTRODUCTION

algorithm for airborne vehicles that can serve as a redundant
navigation system for use during locally limited Global Navigation
Since their deployment, Global Navigation Satellite
Satellite Systems (GNSS) availability, that can be used to validate Systems (GNSS) have set the standard for geodetic positioning.
on-board satellite navigation systems to detect local spoofing The highly engineered nature of these systems can produce
attempts, and that can be used to validate Automatic Dependent geodetic position estimates with errors on the order of 10
Surveillance – Broadcast (ADS-B) position reports. The algorithm meters. In addition, various augmentation methods have been
utilizes the proliferation of ADS-B equipped aircraft as airborne developed to further improve the accuracy of GNSS including
navigation aids in a radio frequency angle-of-arrival (AOA) based but not limited to: Wide Area Augmentation System (WAAS),
geodetic positioning algorithm. Differential Global Positioning System (DGPS), and Real Time
Kinematic (RTK) GNSS. Utilizing these augmentation
The navigation algorithm is loosely based on Simultaneous
methods can improve the accuracy of the GNSS position
Localization and Mapping (SLAM) in that it tracks ADS-B
estimate to the centimeter range. Despite the unparalleled
capable aircraft to refine their ADS-B reported position estimates
while simultaneously determining the geodetic position and
accuracy of GNSS systems, the low power of the satellite-based
velocity of the host vehicle. Unlike SLAM, where the absolute Radio Frequency (RF) signals required to formulate the
location – latitude/longitude – of the landmarks is unknown and geodetic position estimates leaves GNSS susceptible to a lack
must be estimated as the vehicle encounters them, the absolute of availability and spoofing. The lack of GNSS availability may
position of the airborne navigation aids is reasonably well-known be due to natural phenomena, obstructions in the line of sight to
and periodically reported in the ADS-B data set. Because the the satellite constellation, or malicious intent. Spoofing, on the
absolute position of the navigation aids is known, the resulting host other hand, is the intentional introduction of a higher power
vehicle position will also be an absolute, rather than a relative ‘look-a-like’ GNSS signal that causes the GNSS receiver to
position. Secondarily, the continuous tracking of the airborne report an incorrect position estimate. It has been widely
navigation aids allows reported ADS-B positions to be validated theorized that spoofing can be used to ‘take control’ of a GNSS
against the estimated navigation aid position; thereby, guided vehicle; therefore, significant research has been
concurrently accomplishing ADS-B validation and host vehicle performed to detect spoofing attempts.
geolocation.
This paper develops a non-precision geolocation algorithm
It has been demonstrated through simulation that the for airborne vehicles that can serve as a redundant navigation
algorithm is capable of generating valid position estimates and a system for use during locally limited GNSS availability, can be
reliable estimate of its accuracy across a variety of anticipated used to validate on-board satellite navigation systems to detect
input conditions. With multiple GNSS quality navigation aids local spoofing attempts, and can be used to validate Automatic
available, mean position errors well below 500 meters were
Dependent Surveillance – Broadcast (ADS-B) position reports.
observed. As the quality of the navigation aids decreased, so too
did the accuracy of the algorithm. Utilizing navigation aids with
The algorithm utilizes the proliferation of ADS-B equipped
an accuracy of 4 nautical miles (95% containment) resulted in aircraft as airborne navigation aids in an RF angle-of-arrival
mean position errors on the order of 1.25 nautical miles. These (AOA) based geodetic positioning algorithm.
results demonstrate that the method is feasible, and even under The navigation algorithm presented herein is loosely based
worst case conditions, the accuracy of the position estimate on Simultaneous Localization and Mapping (SLAM) in that it
generated by the algorithm was sufficient to allow an aircraft to
navigate to its destination.
tracks ADS-B capable aircraft to refine their ADS-B reported
position estimates while simultaneously determining the
Keywords—geodetic position, geo-location, position fix, ADS-B, geodetic position and velocity of the host vehicle. Unlike
SLAM, angle of arrival SLAM, where the absolute location – latitude/longitude – of the
landmarks is unknown and must be estimated as the vehicle
encounters them, the absolute position of the airborne navigation

978-1-5386-1647-5/18/$31.00 ©2018 IEEE 990

aids is reasonably well-known and periodically reported in the The airborne position message provides the transmitting
ADS-B data set. Because the absolute position of the navigation aircraft’s ICAO address, altitude, latitude, and longitude
aids is known, the resulting host vehicle position will also be an expressed in the WGS-84 ellipsoid [1]. AC 20-165B indicates
absolute, rather than a relative position. Secondarily, the that the transmitting aircraft must transmit its position within 2.0
continuous tracking of the airborne navigation aids allows seconds of measurement [2]; however, the position that is
reported ADS-B positions to be validated against the estimated transmitted must be latent by less than 0.6 seconds, and may not
navigation aid position; thereby, concurrently accomplishing be overcompensated by more than 0.2 seconds [2]. This is taken
ADS-B validation and host vehicle geolocation. Finally, unlike to imply that a position latency in the range [-0.6, 0.2] seconds
is expected. The Radio Technical Commission for Aeronautics
GNSS systems that utilize low power RF signals, ADS-B
(RTCA) recommends a latency of less than 0.4 seconds to
transmissions are relatively high power, ranging from 70 W to
support future ADS-B applications [2].
200 W [1], making them very difficult to jam or spoof.
The airborne velocity message provides the transmitting
In summary, this paper presents a non-precision means of aircraft’s ICAO address, velocity over the ground in a Cartesian
determining an airborne vehicle’s geodetic position, even in the north/east frame, the vertical rate, and a measure of horizontal
absence of GNSS. In addition, this algorithm can be used to velocity error expressed at the 95% probability level. The
validate GNSS position reports in a GNSS spoofed horizontal velocity error, referred to as Navigation Accuracy
environment, and can serve as a means to validate ADS-B Category for Velocity (NACv), has the interpretation given in
position reports from suitably equipped aircraft operating in Table I.
proximity to the host vehicle.
The aircraft operational status message provides information
II. BACKGROUND related to the current status of the aircraft [1]. Of interest to this
application are the ICAO address, the Navigation Accuracy
The algorithm described in this paper requires an Category for Position (NACp), and the Geometric Vertical
understanding of ADS-B, RF AOA, SLAM, direction finding Accuracy (GVA).
position uncertainty, and statistical conversions between
Gaussian and Rayleigh distributions. This section presents an The NACp field provides a measure of the aircraft’s
overview of these areas at a sufficient level of detail to enable Estimated Position Uncertainty (EPU) expressed as a 95% radial
the reader to understand the technical details of the method and value [1]. EPU defines the radius of a circle, centered on the
to support the assumptions relative to the algorithm. reported position, where the probability of the actual position
lying inside the circle is 95% [1]. A cross reference between the
A. Automatic Dependent Surveillance - Broadcast various NACp values and the corresponding EPU is provided in
ADS-B is a modern technology that has been designed to Table II. Similarly, the GVA field describes the one-
enhance air traffic control’s situational awareness of aircraft dimensional, 95% vertical uncertainty as described in Table III.
operations by providing a three-dimensional depiction of each
Each ADS-B transmitter operating on the network emits
ADS-B equipped aircraft’s intended flight path [2]. To
spontaneous broadcasts limited to a maximum of 6.2 messages
accomplish this, each ADS-B equipped aircraft periodically
per second as averaged over a 60 second interval [1]. The
transmits its identification, position, altitude, velocity, and other
aircraft specific information [2]. Ground stations and suitably TABLE I. NAVIGATION ACCURACY CATEGORY FOR VELOCITY [1]
equipped aircraft may receive these transmissions, enhancing air
traffic control and allowing aircraft to maintain adequate NACv Horizontal Velocity Uncertainty (95%)
separation [3]. 0 Unknown or ≥ 10 meters/second
1 < 10 meters/second
The Federal Aviation Administration (FAA) has mandated 2 < 3.0 meters/second
that all aircraft operating within certain airspace segments over 3 < 1.0 meters/second
the United States be ADS-B compliant by January 1, 2020 [3]. 4 < 0.3 meters/second
According to 2015 data released by the Bureau of Transportation
Statistics, there are over 6,800 commercial aircraft and 210,000 TABLE II. NAVIGATION ACCURACY CATEGORY FOR POSITION [1]
general aviation aircraft registered in the United States [4]. A
NACp EPU (95%) NACp EPU (95%)
vast majority of these aircraft will be subject to the ADS-B
0 EPU ≥ 10 NM 6 EPU < 0.3 NM
mandate. The large quantity of ADS-B equipped aircraft make
1 EPU < 10 NM 7 EPU < 0.1 NM
them ideal candidates to serve as airborne navigation aids in the 2 EPU < 4 NM 8 EPU < 0.05NM
proposed navigation system. 3 EPU < 2 NM 9 EPU < 30 meters
4 EPU < 1 NM 10 EPU < 10 meters
Of particular interest to this application are the ADS-B
5 EPU < 0.5 NM 11 EPU < 3 meters
airborne position, airborne velocity, and aircraft operational
status messages. These messages are summarized in the
remaining paragraphs of this section. It is worth noting that each TABLE III. GEOMTERIC VERTICAL ACCURACY [1]
ADS-B message contains the transmitting aircraft’s GVA Vertical Uncertainty (95%)
International Civil Aviation Organization (ICAO) address, 0 Unknown or > 150 meters
which allows the message to be unambiguously associated with 1 ≤ 150 meters
the transmitting aircraft. 2 ≤ 45 meters
3 Reserved

991
transmission of each individual message from the ADS-B TABLE IV. UNIFORMLY DISTRIBUTED ADS-B TRANSMIT INTERVALS [1]
message set is staggered and broadcast at a uniformly distributed Message Minimum Interval Maximum Interval
random interval [1]. The sparseness of transmissions, along with Airborne Position 0.4 seconds 0.6 seconds
the random/staggered interval, reduces the number of data Airborne Velocity 0.4 seconds 0.6 seconds
collisions experienced at the receiver and prevents multiple Operational Status 2.4 seconds 2.6 seconds
aircraft from continuously masking one another due to
transmission synchronization. The uniform distribution
TABLE V. NOTATION RELATIVE TO SLAM
parameters for the messages of interest to this application are
given in Table IV [1]. Control input at time k
Additive zero mean white noise with covariance R
B. Angle of Arrival/Direction Finding Additive zero mean white noise with covariance Q
Initial vehicle state vector at time 0
Angle of Arrival is a method of determining the direction Vehicle state vector at time k
from which an RF wave originated. AOA estimation can be Observation at time k
considered as a far-field plane wave impinging on an array of , ,…, Complete set of all landmarks
antennas [5]. Based on the impinging angle, the signals received U , ,…, Complete history of control inputs
at adjacent antennas will have a phase offset [5]. This phase Z , ,…, Complete set of all observations
relationship can then be used to estimate the AOA of the (⋅) Vehicle motion model function
(⋅) Measurement geometry function
received signal [5].
State transition matrix (Jacobian of )
Christoph Reck and his team of researchers have thoroughly H Measurement matrix (Jacobian of )
studied AOA estimation from ADS-B transmissions. Their
research primarily consisted of comparing the AOA of the ( | , ) (3)
received ADS-B transmission with the expected AOA; where It is common to express the SLAM problem in terms of a
the expected AOA was computed from the receiving antenna state-space model with additive Gaussian noise and to solve the
location and the ADS-B position report [5]. One embodiment of problem using an Extended Kalman Filter (EKF) [6]. In EKF
their research utilized a dedicated target aircraft to validate the based solutions, the motion model of the vehicle, shown in (2),
accuracy of the AOA estimates [5]. The target aircraft contained can be expressed as given in (4); and the observation model,
a highly accurate navigation sensor suite capable of estimating shown in (3), can be expressed as given in (5) [6].
its position to the sub-meter accuracy level [5]. This ensured that
the observed AOA errors were a function of the AOA ( | , )⇔ = ( , )+ (4)
measurement sensor and not induced by errors present in the
transmitting aircraft’s position report. This experiment ( | , )⇔ = , + (5)
empirically demonstrated that an AOA RMSE on the order of
The resulting EKF time update equations are then as
0.66° could be achieved for aircraft operating at a variety of
presented in (6) and (7); while the EKF observation update
ranges [5].
equations are as given in (8) through (10).
C. Simultaneous Localization and Mapping = ( , ) (6)
SLAM is an algorithmic process that allows a mobile vehicle
= + (7)
to generate a map of the environment while simultaneously
using the map to determine the location of the vehicle within the
= + − , (8)
mapped environment [6]. In a SLAM implementation, both the
motion of the vehicle and the location of the landmarks are
=( − ) (9)
estimated without previous knowledge of their locations [6]. An
introductory tutorial of SLAM is available in [6], a summary of = + (10)
which is provided here. Preliminary notion relative to SLAM is
provided in Table V. D. Uncertainty in Airborne Direction Finding
Expressed in probabilistic terms, SLAM attempts to Ancker published a method to determine the uncertainty in a
determine the probability density function (PDF) given in (1) lateral position that was determined from the intersection of
[6]. This density function describes the joint distribution of the multiple lines of bearing measured with respect to airborne
vehicle state and landmark locations given the initial vehicle direction-finding (DF) stations [7]. Ancker used the variance of
state, all control inputs, and all observations [6]. The motion the AOA measurement along with knowledge of each DF
model for the vehicle is assumed to be a Markov process; this station’s position and variance to construct a PDF that represents
implies that the next vehicle state can be determined using only the statistical distribution of the position errors [7]. Ancker
the previous state and control input [6]. The motion model given concluded that the variance of the resulting position fix along an
in terms of a probability is expressed in (2). The observation arbitrary X and Y axis could be expressed as given in (11)
model, given in (3), expresses the probability of the observation, through (14) with parameters defined as:
zk, when the system state and landmark locations are known [6].
n = Number of DF stations used.
( , | , , ) (1) = True bearing angle to station j.
( | , ) (2) = Variance of bearing angle error to station j.

992
 = Distance to station j. utilizing fixed landmarks typical in a SLAM implementation,
= Variance of navigation error of station j. ADS-B capable aircraft will act as airborne navigation aids
(SLAM landmarks). The system determines its best estimate of
= = (11) the position of each airborne navigation aid through fusion of
the reported ADS-B information and the AOA of the received
=∑ (12) signal. Tracking each navigation aid independently allows the
system to detect erroneous ADS-B position reports and account
for them while utilizing the data to determine the position of the
=∑ (13)
host vehicle. Detection of erroneous ADS-B position reports
could then be forwarded to other aircraft and ADS-B ground
=∑ (14) stations, although this is considered beyond the scope of this
paper.
E. Statistical conversions A. Assumptions Relevant to the Method
The accuracy and uncertainty values provided in the ADS-B For purposes of presenting this method, it is assumed that the
data set are expressed at a 95% containment level. This does not ADS-B AOA information will be available referenced to a
satisfy the basic assumption of the EKF-SLAM filter that the locally level, true north frame. That is, the ADS-B AOA receiver
process and measurement noise are uncorrelated, zero mean, and is assumed to be responsible for (capable of) rotating the
normally distributed with covariance Q and R respectively [6]. measurements into a locally level, true north oriented frame
Therefore, these containment values must be adjusted to reflect based on the vehicle orientation. This assumption allows the
a 1-sigma level for use in the filter. EKF-SLAM algorithm to be presented in a more concise manner
The horizontal and vertical uncertainty values indicate the as it does not need to concern itself with the additional burden
uncertainty in a single dimension. As such, these values are of rotating the AOA measurements into a true north frame.
taken to be Gaussian random variables. The Gaussian It is further assumed that utilizing a spherical Earth model in
distribution is well studied and it is understood that 1-sigma place of the more complex and realistic oblate spheroid model
encloses roughly 68% of all samples, while 2-sigma encloses will provide sufficient accuracy for the purposes of this analysis.
roughly 95% of all samples. For purposes of this paper, the The level of relative error due to the Earth’s flattening over the
proportional relationship given in (15) will be assumed. The somewhat small distance to the radio horizon is expected to
letter N in the superscript denotes the normal distribution and the result in an absolute error that can be treated as an uncertainty in
number denotes the containment level percentage. the EKF-SLAM implementation.
= 1.96 (15) B. Filter Formulation
The navigation aid EPU provided in the ADS-B data defines 1) State vector: The state vector for the host vehicle,
a radius of containment. Therefore, it is assumed to be a presented in (20), contains estimates for the host vehicle’s
Rayleigh random variable formulated as the sum of two zero position, velocity, and acceleration in the Earth Centered, Earth
mean, independent Gaussian random variables [8]. The PDF and Fixed (ECEF) frame.
cumulative distribution function (CDF) of a Rayleigh random
variable formed in this manner are given in (16) and (17) It is assumed that the dominant errors present in the ADS-B
respectively, where is the variance of the underlying position and velocity data are due to a combination of latency
Gaussian distribution, and is an integration variable [8]. and navigational errors inherent in the transmitting aircraft’s
navigation sensor suite. This causes these errors to appear as
( )= (16)
biased measurements. As such, the mathematical model of these
errors does not satisfy the fundamental Kalman Filter
assumption that error sources are normally distributed, zero
( )= ( ) =1 − (17) mean Gaussian random variables. Therefore, these measurement
The CDF can then be used to determine the probability of a biases are tracked as filter states. The resulting state vector for
random variable being within n standard deviations of the mean, each airborne navigation aid then contains estimates for the
as shown in (18). Setting the probability to 95% and solving for navigation aid’s position, velocity, position measurement bias,
n indicates that a 95% probability level represents 2.4477-sigma. and velocity measurement bias, all expressed in the ECEF
This establishes the proportional relationship given in (19). frame. The state vector for an arbitrary airborne navigation aid
is given in (21).
( )=1− =1− (18) Finally, the host vehicle state and the state of each observed
airborne navigation aid are expressed as the single combined
= 2.4477 (19)
state vector shown in (22). It is evident that the size of the state
vector is a function of the number of navigation aids being
III. METHOD tracked. Accordingly, the state vector must expand and contract
The method described in this paper implements an to accommodate the changing number of navigation aids within
EKF-SLAM-like algorithm that will track the position and the radio horizon. According to research performed by
velocity of airborne navigation aids while simultaneously Dissanayake, Durrant-Whyte, and Bailey, it is possible to
calculating the geodetic position of the host vehicle. Instead of remove a landmark from a SLAM map “without affecting the

993
statistical consistency of the underlying estimation process” [9]. 0 0 0
Therefore, this expansion and contraction of the state is 0 0 0
considered a viable method to limit the computational = (27)
0 0 ⋱ 0
complexity of the filter. 0 0 0
b) Process noise matrix: The process noise matrix is
= (20) constructed using the piecewise white noise model described in
Where: = Host vehicle position estimate in ECEF. [10]. In this model, the highest order state transition term is
assumed constant for the duration of each time step [10]. The
= Host vehicle velocity estimate in ECEF.
resulting single ECEF dimension process noise for the host
= Host vehicle acceleration estimate in ECEF.
vehicle is presented in (28); while the one-dimensional process
noise for an airborne navigation aid is given in (30). The
= (21) empirically derived terms and represent the variance of
Where: = Navigation aid position estimate in ECEF. acceleration for the host vehicle and navigation aids
= Navigation aid velocity estimate in ECEF. respectively. These values, given in (29) and (31), are functions
= Navigation aid position measurement bias in ECEF. of the navigation aid EPU, adjusted to the 1-sigma level using
= Navigation aid velocity measurement bias in ECEF. the relationship defined in (19). Finally, the host vehicle
process noise matrix and the process noise matrix for each
= ⋯ (22) observed airborne navigation aid are expressed as the single
combined process noise matrix shown in (32).
2) Time propagation: A 10 Hz propagation rate has been
selected for this system. This rate allows a linear state transition
to be assumed because the non-linearity due to the curvature of = (28)
the Earth over the distance traveled in 100 miliseconds at , ,

subsonic speeds is negligible. Given the assumption of local 1

linearity, the a priori estimate of the system state and
covariance are determined using the linear Kalman Filter time
update equations given in (23) and (24), where A is the state = (29)
transition matrix, P is the state covariance matrix, and Q is the
Where: = Mean navigation aid position uncertainty (68%).
process noise matrix.
545 = Empirically derived constant.
= (23)

= + (24) 0 0

a) State transition matrix: The state transition matrix ( ) 0 0

models the linear kinematic equations of motion. For purposes , ,
= (30)
of this model, constant acceleration is assumed for the host 0 0 0
vehicle, while constant velocity is assumed for the airborne
0 0 0
navigation aids. An expression for the state transition matrix for
a single ECEF dimension of host vehicle motion is given in
(25). Similarly, (26) gives the expression for a single ECEF = (31)
dimension of an airborne navigation aid. In both cases, the term
dt represents the time propagation period. The navigation aid Where: = Navigation aid position uncertainty, (68%).
bias states are assumed constant for purposes of time 170 = Empirically derived constant.
propagation. These single dimension expressions are then
extended to three dimensions to populate the corresponding 0 0 0
transition matrix. The host vehicle state transition matrix and 0 0 0
= (32)
the state transition matrix for each observed navigation aid are 0 0 ⋱ 0
then expressed as the single combined state transition matrix 0 0 0
shown in (27). 3) Observation update: Observation updates are also
1 0.5 applied using a 10 Hz frame rate. A measurement update is
, , = 0 1 (25) applied for each navigation aid that has provided fresh ADS-B
0 0 1 data in the preceding 100 miliseconds. ADS-B position and
velocity messages are used to update the corresponding
1 0 0 navigation aid state, while AOA measurements are
0 1 0 0 simultaneously applied to both the host vehicle and the source
, , = (26)
0 0 1 0 navigation aid. Observations are applied using the following
0 0 0 1
order of precedence: position messages, velocity messages,
then operational status messages. That is, if both a velocity
message and a position message have been received from a

994
given navigation aid, then only the position message is used to measurement is given in (39), where the magnitude of the host
update the filter. vehicle position vector is also estimated as part of the pseudo
measurement to hold the host vehicle’s ECEF position on the
a) Measurement vector: The measurement vector is a surface of the Earth.
column vector of the measured data that is to be applied to the
filter. Since it contains measured data, its contents vary as a = = − (37)
function of the type of ADS-B packet being processed. The
AOA measurements are applied independently of the position = = − (38)
or velocity data because position and velocity measurements
are applied only to the airborne navigation aid being tracked; = = (39)
while AOA measurements are simultaneously applied to both
the navigation aid and the host vehicle. Prior to constructing the The i by j measurement matrix given as (40) is then
measurement vector, the received position or velocity data must formulated as the Jacobian of h, where i is the number of rows
be expressed in the ECEF frame to correspond with the in h and j is the number of elements in the state vector. The
coordinate system used by the system state vector. (33) partial derivatives for the non-zero elements of H when an
provides the expression for the measurement vector for an ADS-B position measurement is being applied are given in (41),
ADS-B position message, while (34) gives the measurement (42), and (43). Similarly, the non-zero elements for an ADS-B
vector for an ADS-B velocity message. Finally, (35) shows the velocity message are provided in (44) through (46). The partial
measurement vector for AOA measurements. Note that (35) derivatives for the AOA measurement, provided by Leick et al.,
includes a term for the radius of the Earth. This term applys a are expressed in (47), (48), and (49) [11]. The term s is taken to
pseudo measurement to the host vehicle to hold the ECEF be the magnitude of the vector from the host vehicle to the
position on the surface of the Earth. navigation aid [11]. Finally, the partial derivatives for the
pseudo altitude measurement are given in (50) through (52).
= ̃ ̃ ̃ (33)
Where: ̃ , ,
= Navigation aid ADS-B reported position in ECEF. =∇ = (40)
,

= (34)
Where: = Navigation aid ADS-B reported velocity in ECEF.
, , =− =1 (41)
= (35)
Where: = AOA to navigation aid. =− =1 (42)
= Mean radius of the Earth.
b) Measurement matrix: The measurement matrix relates =− =1 (43)
the state of the system to the measured data. In less formal
terms, the measurement matrix is used by the EKF-SLAM
algorithm to generate a prediction of the expected measurement =− =1 (44)
values based on the current system state. These predicted
(estimated) measurements are then compared with the actual
measurements to update the system state as previously shown =− =1 (45)
in (8).
Position and velocity information for each airborne =− =1 (46)
navigation aid is included in the state vector, allowing their
values to be estimated directly from the filter state. Developing
an AOA estimate from the state vector requires the =− =
( ) ( ) ( ) ( ) ( )
(47)
implementation of the non-linear function given in (36) [11].
( ) ( ) ( ) ( ) ( )
= tan
( ) ( )
(36) =− = (48)
( ) ( ) ( ) ( ) ( )

Where: = Latitude of the surface frame origin (host vehicle). ( ) ( )

=− = (49)
= Longitude of the surface frame origin (host vehicle).
∆ ̂ − ̂
∆ = − = ̂ − ̂ ̂
∆ ̂ − ̂ = (50)

Formulation of the non-linear vector function h varies based ̂

on the type of ADS-B data being processed. The measurement = (51)
vector function for a position measurement is given in (37) and
the measurement vector function for a velocity measurement is ̂
= (52)
provided in (38). Finally, the measurement function for an AOA

995
 c) Measurement noise: The measurment noise matrix all been received from the navigation aid, and the NACp, NACv,
expresses the variance of the measurements being applied to the and GVA values are all non-zero. When a new valid navigation
filter. Again, the contents of the matrix vary based on the type aid is detected, the state vector and covariance matrix are
of measurment being processed. The navigation aid lateral expanded to accommodate the additional navigation aid states.
position measurment uncertainty is established from the EPU,
expressed at the 1-sigma level. The vertical measurement The starting state vector position and velocity values for the
uncertainty is established from the GVA, also expressed at the navigation aid are the ADS-B position and velocity
1-sigma level. As an intermediate step, the lateral and vertical re-referenced to the ECEF frame. The measurement bias values
variance values are expressed in a North-East-Down (NED) are initialized to zero.
surface frame, then rotated in to the ECEF frame for inclusion The initial navigation aid lateral position uncertainty is
in the measurement noise matrix. The final ECEF expression is established from the EPU associated with the received NACp
given in (53). value. The EPU is adjusted from a 95% containment level to a
The navigation aid horizontal velocity measurement 1-sigma normally distributed value using the conversion
uncertainty is a function of the received ADS-B NACv value constant defined in (19). The navigation aid vertical uncertainty
and its related velocity uncertainty, adjusted to the 1-sigma is established from the GVA, where the corresponding vertical
uncertainty is adjusted to a 1-sigma level using the conversion
level. The vertical velocity uncertainty is assumed to be 1.5
times the horizontal uncertainty. The horizontal and vertical constant defined in (15). An additional lateral position
velocity variance values in the NED frame are rotated in to the uncertainty is included during initialization to account for the
ECEF frame for inclusion in the measurement noise matrix maximum ADS-B latency of 0.6 seconds and to account for the
given as (54). time delta between the current system time and the time that the
ADS-B position message was received. Following initialization,
As discussed in the background section, Reck et al. these latencies are accounted for by the position measurement
determined that an AOA RMSE of 0.66° was achievable when bias states. The initial lateral north and east position variance is
the AOA measurement was made with respect to a highly as shown in (56). The lateral and vertical variance values in the
accurate aircraft operating at a variety of ranges [5]. Based on NED frame are then expressed in ECEF for inclusion in the
this empirical result, 0.7° is used as the assumed AOA covariance matrix.
uncertainty as shown in (55). 2
= + (0.6 + ∆ ) (56)
Where: = Navigation aid position uncertainty (68%).
= (53) 0.6 = Maximum ADS-B latency.
Where: = ECEF navigation aid position uncertainty (68%).
∆ = Time since ADS-B position message received.
, ,
= Navigation aid ground speed.

= (54) The preliminary navigation aid velocity variance is a

function of the received ADS-B NACv value and its related
Where: , ,
= ECEF navigation aid velocity uncertainty (68%). velocity uncertainty, adjusted to represent a 1-sigma value. This
value is then assumed for the lateral dimensions, while the
= 0.7 1 (55) vertical dimension is assumed to be 1.5 times the lateral
uncertainty. The NED referenced velocity variance is then
rotated into ECEF for inclusion in the covariance matrix.
C. Filter Initialization
The filter is initialized when two or more navigation aids The initial position and velocity bias variance values were
become available and the angle created by the intersection of the empirically established and are set to 5% and 50% of the
measured lines of bearing from any pair of navigation aids to the corresponding position and velocity variance respectively.
host vehicle is in the range [10°, 170°]. This allowable angle-of-
intersection (AOI) range is somewhat arbitrary and exists to E. Removing Navigation Aids from the Filter
ensure that a unique intersection can be computed from the lines As an airborne navigation aid moves beyond the radio
of bearing. horizon of the host vehicle, the reception of ADS-B messages
from that navigation aid will cease. When an ADS-B message
The initial host vehicle position is calculated as the
hasn’t been received from a navigation aid in the previous 20
intersection of the two lines of bearing from the navigation aid
seconds, the navigation aid is removed from the filter. Removal
pair who’s AOI with the host vehicle is within the allowable
of a navigation aid from the filter requires that all corresponding
range and nearest to 90°. The initial host vehicle lateral
state vector elements be removed and the state vector collapsed.
uncertainty is computed as in (11). The initial position and
Likewise, all corresponding rows and columns of the covariance
uncertainty are then expressed in the ECEF frame for inclusion
matrix are removed and the matrix collapsed.
in the filter state and covariance, where a zero altitude value is
assumed to complete the conversion to ECEF. F. Estimated Position Uncertainty
The instantaneous EPU for the host vehicle, expressed as a
D. Adding Airborne Navigation Aids to the Filter Rayleigh random variable at the 95% containment level, is
A navigation aid is considered valid for use by the filter shown in (57). and are north and east components of
when a position, velocity, and operational status message have

996
the host vehicle position variance estimates rotated from the Finally, the AOA value reported by the simulation defines
ECEF frame to the NED frame. the true north referenced azimuth from the host vehicle to the
navigation aid from which ADS-B data is being simulated. The
= 2.4477 max , (57) AOA value is disturbed by zero mean, uncorrelated, Gaussian
noise with a standard deviation of 0.7°; where the choice of 0.7°
is based on Reck’s empirical results described in the
The final EPU output is produced by smoothing the Background section.
instantaneous EPU using a first order filter employing a 20
second time constant. A comparison between the simulated true position and
ADS-B reported position for an arbitrary airborne navigation aid
IV. SIMULATION is presented as Fig. 1. A NACp value of 2 was selected for this
The method detailed in this paper represents a theoretical navigation aid to emphasize the magnitude of ADS-B reported
algorithmic approach to determine a host vehicle’s geodetic position error; two position bias changes are evident in the
position based on AOA measurements of ADS-B data. Although reported position plot.
Reck’s research has demonstrated the ability to determine AOA
from ADS-B data, all analysis performed to validate this method V. TEST METHOD AND EVALUATION CRITERIA
has been conducted using simulations rather than real world In an effort to enable black box testing of the ADS-B
data. To that end, a software library was developed in C# to geodetic positioning algorithm, a test harness was written in
simulate ADS-B capable aircraft. The simulation allows for the MATLAB that serves as the interface between the simulation
host vehicle and multiple airborne navigation aids to be library and the algorithm. The physical boundary between the
programmed, each with a user defined flight profile and two applications ensures that the algorithm under test only
corresponding navigation accuracy. The simulation publishes operates on the ADS-B data received from the test harness. A
the true position and true velocity information for all defined top level overview of the test harness and its interfaces to the
aircraft, as well as the AOA and pseudo ADS-B data for those simulation library and the algorithm is provided in Fig. 2.
navigation aids that are considered to be in range of the host
The algorithmic method described in this paper presents an
vehicle.
untested approach to geodetic positioning that is highly
The preceding paragraph refers to the information published dependent on navigation aid availability and geometry relative
by the simulation as pseudo ADS-B data. This is due to the fact to the host vehicle. Because of these limitations, there are no
the data content available from the simulation does not strictly clear or predefined expectations for the performance or accuracy
follow the bit-level ADS-B protocol; however, the data is of the algorithm. Therefore, the accuracy of the geodetic
consistent in content, precision, and transmit rate. This method algorithm was evaluated by analyzing the measured radial
has been selected to emulate the data that would likely be position error (RPE) and the EPU output by the filter.
provided by an ADS-B AOA receiver should one be
The terms EPU, MPU, and RPE are used liberally
constructed.
throughout the remainder of this document; they are described
The ADS-B data provided by the simulation is subject to here for convenience:
various perturbations including: latency, bias, and Gaussian
noise. Upon initialization of each simulated aircraft, a random
latency value is selected from a uniform distribution over the
range [-0.6, 0.2] seconds. This latency value remains fixed for
the duration of the simulation run and is applied to all ADS-B
position reports from that aircraft.
The initial aircraft position and velocity bias is also drawn
for each aircraft during initialization. The magnitude of the bias
is selected from a zero mean, two-degree of freedom, chi
distribution with standard deviation equal to the defined 1-sigma
NACp uncertainty (NACv for velocity bias). The direction of
the bias is selected from a uniform distribution over the range
(0, 360] degrees. To emulate an inflight change in GNSS
accuracy due to a constellation change or other phenomena, a
new bias is selected for each aircraft at a 900 second interval. In
order to stagger the bias reselection between aircraft, the first
bias reselection time is randomly selected for each aircraft from
a uniform distribution on the range [0, 900] seconds.
In addition to bias, the reported position and velocity values
are subject to additive, zero mean, uncorrelated, disturbances
generated from the same distributions as the bias values.
However, the standard deviation of the noise is limited to 10%
of the defined 1-sigma NACp or NACv value. Fig. 1. Simulated ADS-B position error with NACp = 2 (EPU < 4 nautical
miles) to emphasize the change in bias.

997
 Each Monte Carlo test consisted of 50-runs, with each run
having a duration of 3600 seconds. Each run utilized a different
seed for the random number generator to ensure that the inputs
to the algorithm varied between runs. The lateral flight profile
for the host vehicle and all navigation aids used for the Monte
Carlo tests is given in Fig. 3.
During each Monte Carlo run, at fixed and predetermined
times, the availability of pairs of navigation aids was altered to
simulate the detection of new navigation aids and the loss of
existing navigation aids. Table VII provides a listing of the
available navigation aids at each inflection point during the test.
Every 15 seconds during each run of the test, the following
host vehicle data was stored: time, EPU, RPE, and number of
Fig. 2. Top level ADS-B simulation test harness block diagram. navigation aids in use. Upon completion of the final run in the
Monte Carlo test, the host vehicle mean EPU, mean RPE, RPE
• Estimated Position Uncertainty (EPU) – The radius of a standard deviation, and MPU were computed for each saved
circle, centered on the reported position, where the time step. To allow evaluation of the filter’s performance, these
probability of the actual position lying inside the circle is TABLE VI. REPORTED ADS-B ACCURACY STATISTICS
95% [1]. This is computed by the algorithm to indicate
the accuracy of the host vehicle position estimate. Parameter Mean Median Maximum Minimum
NACp 9.66 10 11 6
• Radial Position Error (RPE) – The magnitude of the NACv 1.48 1 4 1
lateral position error between the true position and the GVA 1.98 2 2 1
reported position.
• Measured Position Uncertainty (MPU) – Position TABLE VII. SUMMARY OF NAVIGATION AID AVIALABILITY
uncertainty based on the standard deviation of the Time [sec] Navigation Aid Availability
measured RPE, adjusted to the 95% containment level 0 None available
using the conversion constant defined in (19). 100 1 & 2 available
600 1, 2, 3, 4 available
A nearly infinite number of test scenarios could be derived 1100 1, 2, 3, 4, 5, 6 available
to evaluate the accuracy of the navigation algorithm; however, 1600 1, 2, 3, 4, 5, 6, 7, 8 available
as an initial assessment of a somewhat unique navigation 2100 3, 4, 5, 6, 7, 8 available
solution, the test cases selected for this analysis were limited to 2600 5, 6, 7, 8 available
those that represent nominal conditions. Nominal conditions are 3100 7 & 8 available
those cases where airborne navigation aids are available,
providing positions within the bounds of their reported
uncertainty, and located to provide a favorable geometry.
To gauge the nominal reported accuracy of ADS-B
transmissions for use in these test scenarios, a survey the NACp,
NACv, and GVA values from 2417 unique ADS-B capable
aircraft was conducted. The hardware consisted of a NooElec
USB RTLSDR ADS-B receiver and the software was a locally
modified version of the open source Dump1090 Mode S
decoder. The software was configured to record the ICAO,
NACp, NACv, and GVA values for aircraft that were
transmitting non-zero accuracy values. Once recorded, the data
was post processed to determine the mean, median, maximum,
and minimum NACp, NACv, and GVA values. The survey
results are presented as Table VI.
The accuracy of the ADS-B geodetic positioning algorithm
was evaluated against a varying number of available navigation
aids by performing three Monte Carlo tests. The first test
configured all simulated airborne navigation aids to use the
median NACp, NACv, and GVA values from Table VI. The
second test utilized the minimum values given in Table VI. The
final test simulated the case of total GNSS loss by utilizing a
NACp of 2 and a NACv and GVA of 1. Fig. 3. Simulated lateral flight profile. Host vehicle path is in red, airborne
navigation aid paths are numbered and in cyan. Asterisks indicate the starting
position, with the arrow denoting the initial direction of travel.

998
values were then plotted on a single graph along with an the second interval, even though both intervals utilized the same
annotation of the number of navigation aids in use. number of navigation aids.
VI. RESULTS By definition, the EPU is expected to represent a 95%
The results of the median navigation aid accuracy (NACp = probability of containing the actual position. The results clearly
10) test are presented as Fig. 4. Note that there were no indicate that the mean EPU completely contained the mean RPE
navigation aids available to the algorithm for the first 100 for this sample of 50 runs. A maximum mean EPU of 4676
seconds of the flight, upon which time navigation aids 1 and 2 meters and a minimum of 548 meters were observed.
became available. Once these navigation aids were available, the Additionally, the mean EPU completely contained the MPU for
filter self-initialized its position and uncertainty using the all time steps. This is significant because the MPU is taken to be
methods described for filter initialization. the best estimate of the true position uncertainty because it is
calculated from measured data. Ideally, the EPU and MPU
It is clear from Fig. 4 that the accuracy of the algorithm – would be of similar magnitude; these results indicate that the
reflected in both the EPU and RPE – is a function of the number filter’s estimate of its position uncertainty is conservative.
of navigation aids that are available to the filter. Distinct step
changes in the EPU are evident when the number of available Fig. 5 shows the results for the minimum observed
navigation aids changes, with the minimum EPU of 546 meters navigation aid accuracy Monte Carlo test (NACp = 6). Because
occurring when eight navigation aids were in use. Likewise, the the same flight profile is used, the results follow a similar pattern
minimum observed mean RPE of 173 meters occurred when the to the median accuracy case given in Fig. 4. However, due to the
full complement of navigation aids were available; while the reduced navigation aid EPU (556 meters versus 10 meters), the
maximum mean RPE of 2211 meters occurred during resulting host vehicle position errors increased. In this case, a
initialization. maximum mean RPE of 2510 meters and a minimum mean RPE
of 256 meters were observed. This constitutes a 13.5% increase
The more subtle effect of navigation aid geometry relative to in the maximum mean RPE and a 47.9% increase in the
the host vehicle is also evident in Fig. 4. Two navigation aids are minimum mean RPE.
available in the intervals (100, 600] and (3100, 3600]; however,
a mean RPE on the order of 1000 meters was observed during Although the algorithm clearly produced a less accurate
the first interval, while a mean RPE on the order of 350 meters result in this case, it correctly reported the increase in position
was observed during the second. This is largely due to the uncertainty by increasing the EPU. The mean EPU again
positioning of the navigation aids relative to the host vehicle. completely contained the mean RPE, with a maximum mean
EPU of 4727 meters and a minimum of 797 meters. The mean
During the first interval, the navigation aids are on the order EPU also continued to contain the MPU for all time steps.
of 270 kilometers from the host vehicle and their AOA
measurements over the interval result in an average AOI with Fig. 6 provides the results for the total GNSS loss Monte
the host vehicle of approximately 34°. This constitutes a usable, Carlo test (NACp = 2). Again, the same flight profile was used,
but less than ideal geometry. During the second interval, the so the results follow a similar error profile. In this case, the
navigation aids are roughly 130 kilometers away, and have an navigation aid EPU is 4 nautical miles, resulting in an increase
average AOI of 93°. Because the navigation aids are closer to in the measured mean RPE for the host vehicle. The maximum
the host vehicle, and the AOI is nearly at the ideal AOI of 90°, mean RPE was 8578 meters, and the minimum was 1686 meters.
the host vehicle position uncertainty is significantly lower over For the portion of the flight beyond 1100 seconds, the RPE was

Fig. 4. Comparison of estimated position uncertainty, measured position Fig. 5. Comparison of estimated position uncertainty, measured position
uncertainty, and radial position error for the median navigation aid accuracy uncertainty, and radial position error for the minimum observed navigation aid
(NACp=10) Monte Carlo test case. The number of airborne navigation aids in accuracy (NACp=6) Monte Carlo test case. The number of airborne navigation
use is indicated along the top. aids in use is indicated along the top.

999
 computing a valid position estimate with a faithful estimate of
its uncertainty across a variety of expected input conditions.

VII. CONCLUSION
This work has demonstrated that a geodetic positioning
algorithm based on AOA measurements from ADS-B capable
aircraft is technically feasible given a reliable method from
which to compute the AOA. These results also demonstrate that
the EPU computed by the algorithm faithfully represents a
conservative estimate of the actual 95% radial containment
across a variety of expected input conditions.
When six or more navigation aids with an EPU of 4 nautical
miles were used, the algorithm generated a navigation solution
with a mean radial position error on the order of 1.25 nautical
miles and reported an EPU of 1.75 nautical miles. Using six or

Fig. 6. Comparison of estimated position uncertainty, measured position

uncertainty, and radial position error for the GNSS unavailable (NACp=2)
Monte Carlo test case. The number of airborne navigation aids in use is
indicated along the top.

between 1 and 1.25 nautical miles. This level of accuracy is

consistent with an inertial navigation system that has been
navigating for a couple of hours, implying that it is minimally
sufficient for airborne navigation.
The EPU again increased to correctly denote the reduced
navigation aid accuracy, and the EPU again completely
contained the mean RPE. In this case, the magnitudes of the EPU
and MPU are similar, indicating that the EPU became less
conservative as the navigation aid accuracy diminished.
These results demonstrate that the algorithm is capable of
generating a geodetic position estimate and faithfully reporting
its position uncertainty as the number of navigation aids change, Fig. 7. Estimated position uncertainty versus radial position error for a single
the geometry of the navigation aids change, and as the accuracy run of the median navigation aid accuracy level (NACp = 10) test scenario. The
of the navigation aids diminishes. Even in a worst case scenario number of airborne navigation aids in use is indicated along the top.
of total GNSS loss, the accuracy of the position estimate
generated by the algorithm was sufficient to allow an aircraft to
navigate to its destination.
During execution of the nominal accuracy Monte Carlo run,
a random sample was drawn from the set to evaluate the
instantaneous (rather than aggregate) performance of the filter.
The position uncertainty plot for this random single run is given
as Fig. 7. A maximum RPE of 2559 meters and a minimum RPE
of 0 meters was observed, and the EPU again reflected a
minimum 95% containment of the instantaneous RPE.
Although there are no quantifiable performance or accuracy
expectations for this algorithm, it is generally expected that the
reported position will follow the true path of the vehicle without
significant discontinuities, even when the host vehicle
undergoes a change in state. To this end, Fig. 8 provides a
comparison between the true host vehicle position and the
reported host vehicle position as the vehicle encounters a turn.
Fig. 8 demonstrates that the reported position tracks the true
position without significant discontinuities in the path.
Even though the position generated by this algorithm lacks Fig. 8. Comparison of the estimated host vehicle position with the true host
the fidelity of the highly engineered modern GNSS systems, vehicle position with four median accuracy navigation aids (NACp = 10)
these results demonstrate that the algorithm is capable of available.

1000
more GNSS quality navigation aids drastically improved the estimation,” International Journal of Microwave and
performance of the algorithm, leading to a mean radial position Wireless Technologies, vol. 4, no. 2, pp. 181–186, Feb.
error well below 500 meters and a reported EPU below 1000 2012.
meters. [6] H. Durrant-Whyte and T. Bailey, “Simultaneous
Taken together, these results demonstrate that the method is localization and mapping: part I,” IEEE Robotics &
capable of producing a non-precision navigation solution that is Automation Magazine, vol. 13, no. 2, pp. 99–110, Jun.
capable of guiding an aircraft to its destination, even under the 2006.
worst case input conditions of a simulated loss of the GNSS [7] C. Ancker, “Airborne direction finding - the theory of
network. navigation errors,” IRE Transactions on Aeronautical and
Navigational Electronics, vol. ANE-5, no. 4, pp. 199–210,
REFERENCES Dec. 1958.
[8] J. Shynk, Probability, random variables, and random
[1] RTCA, Inc., “Minimum operational performance processes : Theory and signal processing applications.
standards for 1090 MHz extended squitter automatic Hoboken, New Jersey: Wiley-Interscience, 2012.
dependent surveillance - broadcast and traffic information [9] G. Dissanayake, H. Durrant-Whyte, and T. Bailey, “A
services - broadcast,” Washington, DC, DO-260B, Dec. computationally efficient solution to the simultaneous
2011. localisation and map building (SLAM) problem,” in
[2] Federal Aviation Administration, Department of Robotics and Automation, 2000. Proceedings. ICRA ’00.
Transportation, “Airworthiness approval of automatic IEEE International Conference on, San Francisco, CA,
dependent surveillance - broadcast out systems,” Advisory 2000, vol. 2, pp. 1009–1014.
Circular AC 20-165B, Dec. 2015. [10] R. Labbe, “Kalman and Baysian filters in Python,” 15-Jul-
[3] Federal Aviation Administration, Department of 2017. [Online]. Available:
Transportation, “Automatic dependent surveillance - [Link]
broadcast operations,” Advisory Circular AC 90-114A, Hc0SjduNzg/view. [Accessed: 14-Aug-2017].
Oct. 2014. [11] A. Leick, L. Rapoport, and D. Tatarnikov, GPS satellite
[4] U.S. Department of Transportation - Bureau of surveying, 4th ed. Hoboken, NJ: Wiley, 2015.
Transportation Statistics, “National transportation
statistics,” Jul. 2017.
[5] C. Reck, M. S. Reuther, A. Jasch, and L.-P. Schmidt,
“Verification of ADS-B positioning by direction of arrival

1001