USING RANDOM SAMPLING CONSENSUS (RANSAC) TO DETECT ERRORS IN
GLOBAL NAVIGATION SATELLITE SYSTEMS (GNSS) SIGNALS AND DATA
by
Nihar Shah, B.S., M.S.
December, 2021
Director of Thesis: Zhen Zhu, Ph.D.
Major Department: Mechanical Engineering
Abstract
A positioning, navigation, and timing (PNT) signal can be used to estimate a user’s po-
sition at an identified time. A global navigation satellite system (GNSS) uses the PNT
signal to provide satellite-based navigation. Advanced receivers can track multiple GNSS
constellations simultaneously. In order to have a robust and accurate solution, a user needs
to detect any faulty measurements and data, and identify which satellite provided them so
that faulty satellite can be excluded from a GNSS solution. Differencing techniques, such
as time-differenced carrier phase (TDCP), provide for error reduction. The random sample
consensus (RANSAC) method allows for the smoothing of data, even when there are a lot of
gross errors present in the data set. The residuals from RANSAC and TDCP were studied
to determine if they can be used to detect and identify error sources. A downsampling and
thresholding method was able to identify first-order biases with slopes on the order of 10?6
within minutes, while biases with slopes on the order of 10?7 were identified on the order
of one hour. The residuals from RANSAC and TDCP were ultimately able to detect and
identify error sources.

Using Random Sampling Consensus
(RANSAC) to Detect Errors in Global
Navigation Satellite Systems (GNSS) Signals
and Data
A Thesis
Presented to the Faculty of the Department of Engineering
East Carolina University®
In Partial Fulfillment of the Requirements for the Degree
Masters of Science in Mechanical Engineering
by
Nihar Shah, B.S., M.S.
December, 2021
©2021, Nihar Shah
Using Random Sampling Consensus (RANSAC) to Detect Errors in Global Navigation
Satellite Systems (GNSS) Signals and Data
by
Nihar Shah, B.S., M.S.
APPROVED BY:
Director of Thesis
Zhen Zhu, Ph.D.
Committee Member
Thad Wasklewicz, Ph.D.
Committee Member
Teresa Ryan, Ph.D.
Chair of the Department of Engineering
Barbara Muller-Borer, Ph.D.
Dean of the Graduate School
Paul J. Gemperline, Ph.D.
Contents
List of Tables v
List of Figures vi
List of Acronyms viii
1 Introduction 1
2 Literature Review 5
2.1 The Measurement Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Sources of Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3 Error Reduction and Time-Differenced Carrier Phase . . . . . . . . . . . . . 8
2.4 Trajectory Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.5 Random Sampling Consensus . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.6 Time-Variant Stochastic Statistics . . . . . . . . . . . . . . . . . . . . . . . . 11
3 Methodology 14
3.1 Corrected TDCP Measurements . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Correction Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.3 Impact of Satellite Clock and Orbit Errors . . . . . . . . . . . . . . . . . . . 16
3.4 Using RANSAC with Instantaneous Measurements . . . . . . . . . . . . . . 17
3.5 Simulated Bias in Satellite Orbit . . . . . . . . . . . . . . . . . . . . . . . . 18
3.6 Downsampling and Thresholding . . . . . . . . . . . . . . . . . . . . . . . . 19
4 Results and Analysis 21
4.1 SP3 Corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.2 Autocorrelation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.3 Allan Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.4 Elevation Angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.5 Autocorrelation of Truncated Data . . . . . . . . . . . . . . . . . . . . . . . 32
4.6 Allan Variance Analysis of Truncated Data . . . . . . . . . . . . . . . . . . . 34
4.7 Downsampling and Thresholding . . . . . . . . . . . . . . . . . . . . . . . . 37
4.8 Simulated Bias in Satellite Orbit . . . . . . . . . . . . . . . . . . . . . . . . 45
5 Conclusion 53
5.1 Detect and Identify Error Sources . . . . . . . . . . . . . . . . . . . . . . . . 53
5.2 Recommendations for Future Work . . . . . . . . . . . . . . . . . . . . . . . 54
References 56
List of Tables
2.1 Definition of terms in Equation 2.1 . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Expected log-log slopes of Allan variance for common noise types . . . . . . 12
3.1 Definition of terms in Equation 3.1 . . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Definition of terms in Equation 2.1 to account for TDCP . . . . . . . . . . . 15
4.1 Summary statistics of differences between original and Standard Product #3
(SP3)-corrected residuals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 Time (in s) to identify simulated bias of 1× 10?4 . . . . . . . . . . . . . . . 48
4.3 Time (in s) to identify simulated bias of 5× 10?5 . . . . . . . . . . . . . . . 48
4.4 Time (in s) to identify simulated bias of 1× 10?5 . . . . . . . . . . . . . . . 48
4.5 Time (in s) to identify simulated bias of 5× 10?6 . . . . . . . . . . . . . . . 49
4.6 Time (in s) to identify simulated bias of 1× 10?6 . . . . . . . . . . . . . . . 49
4.7 Time (in s) to identify simulated bias of 5× 10?7 . . . . . . . . . . . . . . . 49
4.8 Time (in s) to identify simulated bias of 1× 10?7 . . . . . . . . . . . . . . . 49
List of Figures
1.1 Finding two positions O in two-dimensional space using range from two satellites 3
2.1 Faulty orbital model leading to orbital error . . . . . . . . . . . . . . . . . . 7
2.2 Antenna receiving a direct path signal and a reflected path signal, leading to
multipath error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.1 Overview of satellite and receiver motion . . . . . . . . . . . . . . . . . . . . 14
4.1 Comparison of original and SP3-corrected residuals for satellite 14 . . . . . . 22
4.2 Comparison of original and SP3-corrected residuals for satellite 31 . . . . . . 22
4.3 Normalized autocorrelation for satellite 14 . . . . . . . . . . . . . . . . . . . 23
4.4 Normalized autocorrelation for satellite 14 with primary peak highlighted . . 24
4.5 Normalized autocorrelation for satellite 31 . . . . . . . . . . . . . . . . . . . 25
4.6 Normalized autocorrelation for satellite 31 with primary (red) and secondary
(blue) peaks highlighted . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.7 Allan variance analysis for satellite 14 . . . . . . . . . . . . . . . . . . . . . . 27
4.8 Allan variance analysis for satellite 31 . . . . . . . . . . . . . . . . . . . . . . 28
4.9 Elevation angle of satellite 14 . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.10 Corrected residuals from satellite 14, channel 1 with a 250-point centered
moving average . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.11 Corrected residuals from satellite 14, channel 1 with a 250-point centered
moving average and approximation of the moving average . . . . . . . . . . . 29
4.12 Elevation angle of satellite 31 . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.13 Corrected residuals from satellite 31, channel 1 with a 250-point centered
moving average . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.14 Corrected residuals from satellite 31, channel 1 with a 250-point centered
moving average and approximation of the moving average . . . . . . . . . . . 31
4.15 Normalized autocorrelation for truncated satellite 14 residuals . . . . . . . . 32
4.16 Comparison of original residuals and truncated residuals for satellite 14 . . . 33
4.17 Normalized autocorrelation for truncated satellite 31 residuals . . . . . . . . 33
4.18 Comparison of original residuals and truncated residuals for satellite 31 . . . 34
4.19 Allan variance analysis for truncated residuals from satellite 14 . . . . . . . . 35
4.20 Allan variance analysis for truncated residuals from satellite 31 . . . . . . . . 36
4.21 Variation of initialization and window sizes for satellite 14, Part A . . . . . . 38
4.22 Variation of initialization and window sizes for satellite 14, Part B . . . . . . 39
4.23 Percentage of points identified as outliers versus initialization period . . . . . 40
4.24 Percentage of points identified as outliers versus window size . . . . . . . . . 40
4.25 Variation of initialization and window sizes for satellite 31, Part A . . . . . . 42
4.26 Variation of initialization and window sizes for satellite 31, Part B . . . . . . 43
4.27 Percentage of points identified as outliers versus initialization period . . . . . 44
4.28 Percentage of points identified as outliers versus window size . . . . . . . . . 44
4.29 Variation of added bias, initialization period, and window size for residuals
for satellite 14, Part A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.30 Variation of added bias, initialization period, and window size for residuals
for satellite 14, Part B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.31 Variation of added bias, initialization period, and window size for residuals
for satellite 31, Part A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.32 Variation of added bias, initialization period, and window size for residuals
for satellite 31, Part B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
List of Acronyms
ADR accumulated Doppler range
AOA angle of arrival
BDS BeiDou Navigation Satellite System
CA coarse/acquisition
CDDIS Crustal Dynamics Data Information System
CEI clock, ephemeris, integrity
CL civil-long
CM civil-moderate
GLONASS Global’naya Navigatsionnaya Sputnikovaya Sistema
GNSS global navigation satellite system
GPS Global Positioning System
IGS International Global Navigation Satellite System Service
INU inertial navigation unit
IRNSS Indian Regional Navigation Satellite System
k-NN k-nearest-neighbor
NASA National Aeronautics and Space Administration
NavIC Navigation with Indian Constellation
P precision
PNT positioning, navigation, and timing
PRN pseudorandom noise
QZSS Quasi-Zenith Satellite System
RANSAC random sample consensus
RMS root-mean-square
RSS relative signal strength
RTK real-time kinematic
SP3 Standard Product #3
SPP single point positioning
SWaP size, weight, and power
TDCP time-differenced carrier phase
TDOA time difference of arrival
TOA time of arrival
TOF time of flight
Chapter 1
Introduction
A positioning, navigation, and timing (PNT) radio signal can be used to estimate a user’s
navigation solution. A PNT receiver processes the transmitter’s signal to extract ranging
measurements and other parameters of interest, which are used to compute a user’s position
and/or clock solution.29
The general term for satellite-based navigation systems is global navigation satellite sys-
tem (GNSS). In the United States, our primary GNSS is Global Positioning System (GPS).
It was designed for PNT of aircraft and other types of vehicles. Additionally, GPS is critical
in less-obvious uses, such as banking, telecommunications, weather, seismology, power syn-
chronization, and atmospheric research.23 The need for accurate civilian GPS service was
addressed by President Clinton as early as 2000.10 The problems with GPS accuracy, avail-
ability, integrity, and reliability were identified at the Presidential level as early as 2004.7
The impacts of degraded GPS continue to affect all levels of national security, as identi-
fied in multiple national strategy documents that are applicable at the time this thesis was
authored.43,45,46 Other GNSS systems around the world include:24
? BeiDou Navigation Satellite System (BDS), owned by China
? Galileo, owned by the European Union
? Global’naya Navigatsionnaya Sputnikovaya Sistema (GLONASS), owned by the Rus-
sian Federation
The Indian Regional Navigation Satellite System (IRNSS)/Navigation with Indian Constella-
tion (NavIC) (owned by the Government of India) and Quasi-Zenith Satellite System (QZSS)
(owned by the Government of Japan) operate regionally over India and Japan, respectively,
but not globally.24
As of 2019, there were over 120 GNSS satellites worldwide.51 The United States’ GPS
transmits signals at two frequencies (L1 = 1575.42 MHz and L2 = 1227.6 MHz, or 154 and
120 times the internal clock frequency of 10.23 MHz, respectively).32,33,44,52 Each signal is
encoded by mod-2 addition (can also be considered to be a logical “exclusive or” or XOR)
of a pseudorandom noise (PRN) ranging code and navigation data.33 The L1 signal has two
codes in quadrature: a coarse/acquisition (CA) and a precision (P) code.32 The P -code also
has an anti-spoofing mode that is reserved for military use.32 The L2 signal provides a more
accurate timing signal using the L2 civil-moderate (CM) and L2 civil-long (CL) codes.32
The L5 signal is the newest signal that has a faster acquisition time and provides increased
accuracy by implementing a new ranging signal, synchronization, and civil navigation data.33
To acquire signals from the satellites, receivers use the predicted position calculated with
navigation data.29 Receivers use a locally-generated copy of the code to determine the time
of reception. Based on the fact that the signal travels at a known speed (the speed of light),
receivers are able to calculate the distance between the transmitter and receiver.29 Two
satellites (S1 and S2) can be used to calculate user position (U) in the two-dimensional case,
as illustrated in Figure 1.1.57
There are two possible points that result from the calculation with S1 and S2; the closer
of the two points is defaulted to being the correct point.29 A third satellite is needed to
calculate user position in the three-dimensional case.29 All of the calculations assume that
the clocks on the transmitter and receiver are perfectly synchronized; a fourth satellite is
2
Figure 1.1: Finding two positions O in two-dimensional space using range from two
satellites
needed to resolve the difference in satellite and user timing in the three-dimensional case.29,57
At any given time, a user in the United States can expect to have six to eight GPS satellites
within line of sight.16
Advanced receivers can track multiple GNSS constellations simultaneously.66 All of these
signals can be integrated to increase the accuracy and reliability of PNT with a minimal
increase in size, weight, and power (SWaP); and cost.66 Since more GNSS constellations can
be observed, there is greater redundancy in satellite signals. Integrating and comparing data
from different GNSS allows for the identification of errors.11
In order to have a robust and accurate solution, a user needs to detect and identify any
faulty measurements and data, so that faulty satellite(s) can be excluded from a GNSS solu-
tion.66 For example, multipath is considered an error source for GNSS ranging measurement.
However, if a reflected signal is received instead of the direct line of sight signal, it is con-
sidered a faulty measurement. Another example is error in satellite clock and orbit models:
if the broadcast satellite model has an error that drifts beyond acceptable tolerance, the
satellite must be identified and detected. Satellite faults can be observed in measurement
residuals.11,60,66
Some error sources can be reduced via differencing. For example, time-differenced carrier
phase (TDCP) can reduce some of the biases in the solution calculations by differencing
3
successive carrier phase measurements.66 Additionally, methods such as random sample con-
sensus (RANSAC) can be used with multi-GNSS receivers to detect and remove bias in
TDCP measurements.66 The two research questions for this thesis are:
1. Can residuals from RANSAC and TDCP be used to detect and identify
error sources?
2. If so, what types of errors can it detect?
Chapter 2 will review the literature and concepts used for this research. Chapter 3 will
summarize the analysis methodology. Chapter 4 will summarize the analysis results. Finally,
Chapter 5 will summarize the conclusions and recommendations.
4
Chapter 2
Literature Review
2.1 The Measurement Equation
The general measurement equation of a GNSS navigation system is given by Equa-
tion 2.1.67 The terms in Equation 2.1 are defined in Table 2.1.67
z = Hx+ ?+ f (2.1)
Term Definition
z Observed position measurements
H Geometric matrix
x True state
? Nominal errors
f Faults
Table 2.1: Definition of terms in Equation 2.1
The nominal errors ? are errors that can be modeled. The faults f include biases in z
and/or H, that are not captured by ?.
2.2 Sources of Error
Previous work has identified six nominal sources of GPS carrier phase error ? in calcu-
lating user position:50,56,67
1. Ionospheric delay
2. Tropospheric delay
3. Noise
4. Satellite orbit prediction
5. Inter-system bias
6. Signal multipath
Significant variations in electron density in the ionosphere result in a contribution to
ranging errors.61 Delays in ranging code and change of carrier phase lead to ranging er-
rors, which can aggregate as trends.29 The ionospheric delay is frequency-dependent, so a
dual-frequency receiver can be used to cancel out this error; however, the process magnifies
any other measurement errors present in the solution.29 If using a single-frequency measure-
ment, the ionospheric delay may contribute anywhere from 9.8 m to 19.6 m to the overall
error.44 Similarly, tropospheric delay can increase ranging errors up to 2.4 m if they are not
corrected with a model.29,38,58 Ionospheric models are able to correct about 50% to 70% of
the ionospheric error, while tropospheric models can correct over 90% of the tropospheric
error.29 Both types of atmospheric effects can be modeled as biases or first-order deviations
in TDCP.34
Different sources of noise caused by other gradients can further contribute to error.66
Thermal noise is the dominant source of noise, and is evaluated with the carrier-to-noise ratio
(C/N 290). The C/N0 is affected by the elevation angle between the satellite and the user, since
the signal interacts with more of the atmosphere at low elevation angles of 15? or lower.29
The low elevation angle can lead to error contributions on the order of 10 m.29 Physical
6
environmental factors (e.g., thermal noise) can be modeled as a random walk process, while
non-stationary noise can be modeled as Gauss-Markov processes.62,68
Errors in satellite orbit prediction (caused by inconsistent data, spoofing, drifts, etc.) can
result in large biases in the user positioning solution, as illustrated in Figure 2.1. Between
2008 and 2014, the average nominal error in orbit prediction for GPS satellites was 52 cm.29
Faults in satellite cock or orbit can be significantly greater. For example, on April 1, 2014,
GLONASS satellites transmitted incorrect broadcast messages that resulted in errors on the
order of 200 km.3 The erroneous positioning data was likely caused by a “bad” software
update.3 A sudden change in the prediction model would appear as a bias or first-order
deviation in TDCP.34
Figure 2.1: Faulty orbital model leading to orbital error
Inter-system bias is caused by timing differences between different GNSS constellations.56
Differencing methods can be used to reduce this error, and are further described in Sec-
tion 2.3.
The broadcast signal from the GPS satellites can be reflected or diffracted just like any
other type of electromagnetic wave, leading to a GNSS ranging error source known as multi-
path.31 Multipath is illustrated in Figure 2.2.54 Since the reflected signal (red dashed line in
Figure 2.2) takes longer to reach the antenna, the range calculation produces an erroneous
value (blue dot at the end of the dotted line in Figure 2.2).31 Statistically, multipath errors
can be modeled as a Gauss-Markov stochastic process.30 Further, the multipath error model
used by an autonomous vehicle in a city will be different than that of a commercial airliner
flying over an ocean. Since multipath error is so highly dependent on the environment, there
7
are not unified models to mitigate it like there are for tropospheric delay.31
Figure 2.2: Antenna receiving a direct path signal and a reflected path signal, leading to
multipath error
2.3 Error Reduction and Time-Differenced Carrier Phase
Differencing techniques provide for error reduction. Single differencing the carrier phase
signal compares the received signal from one satellite to two antennas.20 Single differencing
cancels out the common effects shared by signals from a satellite taking different paths, such
as atmospheric effects like ionospheric delay and tropospheric delay.20 In double differencing,
two single-differenced signals are subtracted from each other.1,59 Double-differencing reduces
errors in satellite orbit prediction, clock drift, and receiver clock bias.20
Using the PRN code to determine user position is a single point positioning (SPP) so-
lution, which provides positioning accuracy on the order of a few meters.5,29 Incorporating
the carrier-phase observations can allow for a position accuracy on the order of centimeters
or less; however, the carrier-phase measurement is ambiguous because the number of inte-
ger multiples of a complete cycle is unknown.5 One of the ambiguity resolution techniques
8
is real-time kinematic (RTK), which uses some type of base station from which the mea-
surements are differenced against those from the user.29 The RTK method does not require
hours for the correct solution to converge, but it does require precise estimates of the initial
carrier-phase ambiguity to work effectively.5
Alternatively, instead of differencing a signal received by two antennas, successive carrier
phase measurements are differenced in the time-differenced carrier phase (TDCP) method.59,66
With TDCP, range change can be measured by differencing accumulated Doppler range
(ADR) over time. Using pseudorange instead of ADR leads to a more noisy result by a
factor of about 10.64 With the TDCP method, inter-constellation timing bias and several
other error sources are canceled out, allowing one measurement equation to be used for all
satellites and constellations.66 With more available satellites, the solution would offer high
reliability in near real-time, but would have a higher computational cost.66
Using TDCP brings about two new problems. First, there are more data sources to filter
through to find faulty signals. Second, since several biases are differenced out, biases that
remain close to constant over time may be obfuscated.
2.4 Trajectory Correction
The accuracy of the SPP or TDCP solution is directly related to the accuracy of the
GNSS orbit and clock information used to calculate the solution.40 Models of a satellite’s
motion can be used to compute a reference trajectory which includes satellite ephemeris and
clock estimates.29 A least-squares-fit routine converts the predicted positions (based on the
reference trajectory and the position estimates) into positioning solutions.32
Each of the GPS satellites provides a unique clock, ephemeris, integrity (CEI) data that
is needed to provide a positioning solution.32,33 The CEI data includes information such as
satellite-specific clock correction polynomial parameters, ephemeris parameters, and health
9
flags.32,33 These calculations are typically accurate from 5 m to 10 m.28
The International Global Navigation Satellite System Service (IGS) calculates high-
accuracy satellite orbit solutions and provides them to National Aeronautics and Space
Administration (NASA) Crustal Dynamics Data Information System (CDDIS) for archival
in a format called Standard Product #3 (SP3).40 The SP3 format provides satellite orbit
data that is either measured from the broadcast signal, based on an interpolation, or extrap-
olated (or predicted) from an orbit model.27 The IGS orbits are calculated using a weighted
linear combination of solutions from up to eight separate analysis centers.2,26 First, differ-
ences in reference frames are removed from each center’s solution.2 Then, a seven-parameter
Helmert transformation is estimated for each solution to further reduce errors caused by
reference frame misalignment.2
The IGS has three forms of the data: ultra-rapid, rapid, and final.40 The ultra-rapid
data are comprised of fits to observational data for the first half and predicted orbits for
the second half.26 The rapid data are daily solutions that are released approximately 17 hrs
after the end of each day.40 The final data are the most accurate solutions and are weekly
solutions released approximately 13 days after the end of each week.40
A user can use a simple GNSS receiver to collect pseudorange and carrier phase and
post-process the received measurements with corrected satellite position and clock data to
produce a more accurate solution.48 The SP3-corrected data has positional accuracies on the
order of 20 mm,26 which is 500 times more accurate than the accuracy of the CEI data alone.
Therefore, the SP3-corrected data can be used as truth data for comparative analyses.
2.5 Random Sampling Consensus
Traditional methods for fitting models to data (such as the method of least squares)
assume that there are enough “good” values to smooth out any large deviations.17 However,
10
this is not always true for some applications, include GNSS. The RANSAC method allows
for finding an error-free subset, even when there are multiple gross errors present in the data
set.17 The general RANSAC method helps identify the inliers ? from a set of N observations
to solve for state vector x, even if there are multiple simultaneous outliers.67
The process was initially developed for image analysis and automated cartography, but
was then applied to generalized model fitting to estimate unknown parameters in data with
numerous outliers to improve fault detection.17,66 The application of RANSAC methods
has been demonstrated in single GNSS constellations,8,49 but due to the limited number of
satellites, there is minimal benefit from applying RANSAC.66 In a multi-GNSS problem,
the use of RANSAC would have a more significant impact.66 RANSAC can also be used to
identify faulty measurements in GNSS.
2.6 Time-Variant Stochastic Statistics
Each GNSS satellite suffers from errors such as multipath, receiver noise, interference,
troposphere error, and ionosphere gradients.14 While TDCP can minimize these errors, there
will be some residuals from each signal.56
The root-mean-square (RMS) residuals between each individual satellite’s TDCP and
calculated truth reference value can be analyzed to identify different types of trends.56,66
There are five main random processes commonly related to GNSS data: Quantization noise,
White noise, Flicker noise, Random walk, and Gauss-Markov process.15,41 Quantization noise
comes from analog to digital conversions.25 Random white noise will have very low autocor-
relation, and can often be modeled using Gaussian distributions with a long-term average of
zero.35,47,63 Flicker noise is a measure of the bias offset.35,47 The random walk errors are “a
random process of uncertain origin, possibly known as a very low-frequency noise term”47
possibly caused by tropospheric delays.35 Gauss-Markov processes are processes that are
11
Table 2.2: Expected log-log slopes of Allan variance for common noise types
Noise type Curve Slope
Quantization noise ?1
White noise ?1/2
Flicker noise 0
Random walk +1/2
Gauss-Markov process ?1/2then+ 1/2
exponentially autocorrelated with time.47
A commonly accepted method for studying the residuals is by using Allan variance.41,47
The Allan variance method was initially used to study frequency stability in clocks and
electronics, but was later applied to navigation systems to identify processes in random drift
errors over time, to characterize inertial errors, and to identify stochastic processes.41,42,65
One advantage of the Allan variance is that it converges for most types of clock noise based
on the resulting slope, as summarized in Table 2.2.19,36,41,50 Additionally, it identifies the
impact of stochastic processes when using different time periods to average the data.19 Time
periods where the Allan variance is low identify averaging periods where the impacts of of
stochastic processes are minimized.19
Kalman filters are often applied to GNSS data to combine different data sources; however,
the Kalman filter requires the inputs to be white noise.41 Instead of assuming that the data
is white noise (or approximates white noise), the Allan variance analysis can be used to
verify the validity of the white noise assumption.41 If there are other stochastic processes
present besides white noise, the Allan variance analysis would be able to identify them and
allow the user to compensate the model accordingly.41
The work in this thesis attempts to identify zeroth- and first-order biases (points and
linear biases), not traditional GNSS anomalies and faults. There are several methods that
can be used to detect these biases. There are three general categories of detection meth-
ods: supervised, semi-supervised, and unsupervised.22 Supervised methods require a labeled
12
set of training data of “good” and “bad” data.9 Semi-supervised methods have labels for
only “good” data.9 Unsupervised methods do not require training data, and are the most
applicable to this research.9.
A simple global bias detection method identifies the distance a point is from the k nearest
neighboring data points, and is appropriately named k-nearest-neighbor (k-NN) method.22
The value of k used directly affects the effectiveness of the k-NN method. If k is too low,
there are not enough data points to produce a reliable calculation; if k is too large, the
estimate may be too coarse.22 Typically, k is in the range 10 < k < 50.22 For GNSS bias
detection, the balance of low and high k values is a balance of effectiveness and timeliness.
13
Chapter 3
Methodology
3.1 Corrected TDCP Measurements
As illustrated in Figure 3.1, the jth satellite in the kth constellation moves over a time
period from position S(t0) to S(t1). The i
th ground antenna/receiver moves from position
U(t0) to U (t1). The unit vector e
j
i identifies the direction from the ground receiver to the
satellite, while Rji gives the distance.
Figure 3.1: Overview of satellite and receiver motion
The corrected TDCP measurements from the ith antenna/receiver to the lth channel
on the jth satellite in the kth constellation is described in Equation 3.1.50,66 The terms in
Equation 3.1 are defined in Table 3.1.50,66
ADRj,k,li (t1)?ADR
j,k,l
i (t0)? corrSV ? corr
j
geo = ?ei (t1)?U+?bi +?f
j,k,l
i + ?
j,k,l
i (t) (3.1)
Term Definition
ADR ADR
corrSV Satellite Doppler correction term
corrgeo Geometric correction term
eji Unit vector from the i
th antenna/receiver U to the jth satellite S
?U Displacement of the receiver from time t0 to t1
?bi Receiver clock bias
?f j,k,li Contribution of faults on the ADR measurement
?j,k,li Noise and all remaining residual errors after time differencing
Table 3.1: Definition of terms in Equation 3.1
The residual errors ?j,k,li are comprised of time-variant sources such as troposphere resid-
uals ?T ji , ionosphere residuals ?I
j,l
i , inter-system bias ??
j,k,l, multipath errors ?M j,k,li , and
?
noise ?? j,k,l 50i . These individual components to the overall residual errors can be assumed
to be rapidly changing, but with a combined magnitude of error on the order of millimeters
per second with a one-second differencing.50
The measurement equation described in Equation 2.1 can be reused to account for the
TDCP methodology. Table 2.1 would be redefined to Table 3.2.
Term Definition
z ADR
H Geometric matrix
x TDCP measurements
? Nominal errors
f Faults
Table 3.2: Definition of terms in Equation 2.1 to account for TDCP
15
3.2 Correction Terms
The TDCP reflects changes in distances in line of sight. It needs to be corrected with
satellite motion and geometric change between the satellite and user motion directly. The
satellite correction term and the geometric correction term are given by Equation 3.2 and
Equation 3.3, respectively.59
corr jSV = S (t1) · eji (t1)? Sj(t
j
0) · ei (t0) (3.2)
corrgeo = Ui(t1) · eji (t1)?Ui(t ) · e
j
0 i (t0) (3.3)
The generalized unit vector eji (t) is defined in Equation 3.4.
59
??Sjj (t)?Ui(t)e (t) = ?i (3.4)Sj(t)?Ui(t)?
3.3 Impact of Satellite Clock and Orbit Errors
The distance between the satellite and the antenna/receiver is typically ranges from
6400 km to 67,000 km,32 so the impact of a small satellite position error on eji (t) will be
negligible. Therefore, Equation 3.2 can be rewritten as Equation 3.5. In Equation 3.5,
?Sjb(t) · e
j
i (t) represents the satellite position error projected onto the unit vector at the
identified time.
?corrSV = ?S
j
b(t1) · e
j
i (t1)??S
j
b(t0) · e
j
i (t0) (3.5)
If the satellite clock or orbit model are erroneous, there will also be an additional time-
differenced satellite clock error of ??j,k,lb . This would yield the satellite fault model in
Equation 3.6. ( )
?f j,k,l = ??j,k,li,SV b + ?S
j
b(t1) · e
j
i (t1)??S
j
b(t0) · e
j
i (t0) (3.6)
16
The term ?f j,k,li,SV is a component of ?f
j,k,l
i from Equation 3.1. The error calculated by
?j,k,li (t) is the error in satellite based on truth data. If real-time data is aggregated, there is
no direct estimate of truth ?U . Using RANSAC, the instantaneous solutions of ?U? and
?b?i can be used to calculate ?
j,k,l 50
i (t) in Equation 3.7. Because the instantaneous solutions
include error contributions from other satellites, the ?j,k,li (t) error includes contributions from
other satellites too. The statistic characteristics of ?j,k,li (t) can be used to determine how
well ?f j,k,li,SV can be observed.
ADRj,k,li (t1)? ADR
j,k,l
i (t0)? corrSV ? corr
j j,k,l j,k,l
geo + ei (t1)?U???b?i = ?fi + ?i (t) (3.7)
3.4 Using RANSAC with Instantaneous Measurements
An estimate of ?U can be used to observe ?f j,k,l j,k,li,SV and ?i (t) in this work. The term
?U can be reliably estimated from RANSAC. RANSAC can be used to identify outliers
in the solution of ?U . This allows for identification of data channels or satellites that are
providing the outlier data. These data sources can be excluded from the positioning solution
to improve accuracy. To do so, the following algorithm can be applied over ? iterations:67
1. For iteration ?, M measurements are randomly selected as an initial core set s??. Since
there are four unknowns in the TDCP solution, M = 4 in this case. The initial core
set s?? is used to solve for the initial core solution x??. If the solution is erroneous, skip
the iteration.
2. Compute the error function Err (zn,Hn, x??). The measurements with an error function
smaller than some threshold ? would be considered inliers.
3. Record the size of the set of indices of all inliers in iteration ? as criterion C?.
4. If C? is greater than the current maximum Cmax, solve for x?? using all the inliers.
17
Then, record the following parameters: ?? = ?; Cmax = C?; x??? = x??; x??? = x??; s?? = s?;
and s? ??? = s?.
The error function Err (zn,Hn, x??) from step 2 will follow a 1-degree of freedom Chi
distribution.67 If the error function is less than the selected ? (? = 2 by default), a fault-free
measurement has a 95% chance of falling into a domain D? that is based on x??, assuming
that D 67? is unbiased. The measurements that are not included in D? are part of the inlier
set. If a large portion of the measurements are fault-free, they will agree with each other
when compared to an unbiased D?. With RANSAC, a majority of the measurements are
not required to be fault-free; it only needs at least one fault-free core solution through the
? iterations to form a valid solution.67 The likelihood of RANSAC failure for different levels
of faulty measurements and the overall effectiveness of the algorithm is further explored in
Zhu et al. 66 .
3.5 Simulated Bias in Satellite Orbit
As described in Section 2.2, the different sources of error would manifest as biases of
different orders in TDCP. As previously mentioned, the atmospheric effects would manifest
as a zeroth-order bias in TDCP.29 The difference between the real-time residuals and the
SP3-corrected residuals may aggregate over time and appear as a first-order bias in TDCP.29
Under normal operations for GPS, the drift is expected be less than 0.006 m/s.44
For a more controlled test, first-order biases of different levels were added to the the
existing data for satellite 14 to simulate a greater drift in satellite orbit. The bias started
20% of the way into the whole data set and increased linearly with a defined slope. The
downsampling and thresholding analysis was re-accomplished on the residuals with added
simulated biased. The following slopes were used: 1 × 10?4, 5 × 10?5, 1 × 10?5, 5 × 10?6,
1× 10?6, 5× 10?7, and 1× 10?7. The smaller slopes represented nominal errors, while the
18
larger slopes represented phenomena like faults, broadcast errors, spoofing, etc.
3.6 Downsampling and Thresholding
The SP3-corrected data was divided into “windows” containing the following number of
points in each window: 10, 25, 50, 60, 75, 100, 150, 250. Note that there is one point every
second, so the window size can also be considered temporally (i.e., 10 s, 25 s, 50 s, 60 s,
75 s, 100 s, 150 s, 250 s). The mean of each window was used to downsample the data.
The residuals were not assumed to be independent over time and identically distributed,
nor were they assumed to be Gaussian, so standard deviation would not serve as a valid
metric. One way to evaluate temporal correlation is autocorrelation. Autocorrelation evalu-
ates how closely points are statistically related.55 The autocorrelation of white noise (which
is analogous to an independent and identically distributed continuous functions) would be a
single spike.53 If there was a peak with a finite width in the resulting data, there was some
sort of statistical dependence between the points, thereby validating the assumption that the
residuals were not independent and identically distributed. The peak halfwidth identifies the
time lag in the autocorrelation.6 The time lag may be helpful in determining ideal window
sizes for analysis.
Another method of analyzing time-variant data to identify stochastic processes is the
Allan variance described in Section 2.6. The Allan variance was calculated for satellite 14
and satellite 31 following the process described in Shah et al. 50 .
The initial windows were used to calculate thresholds that was the mean of those down-
sampled window means and three times their RMS.29The number of initial windows selected
were: 3, 5, 10, 15, 20, and 25. Any downsampled window that had an average value that
was more than three times the RMS away from the mean of the initialization period was
flagged as an outlier. Since each combination of initialization period and window size had
19
a different number of points, comparing the number of outliers was not a valid comparison.
Instead, the percentage of points identified as outliers was used as the metric for analysis.
20
Chapter 4
Results and Analysis
4.1 SP3 Corrections
Two satellites were chosen for initial analysis, since they were in view for the majority
of the flight. The original residuals for satellite 14 and 31 are illustrated in Figure 4.1
and Figure 4.2. Also illustrated on those figures are the residuals with the SP3 corrections
applied and the difference between the two measurements for each satellite. The summary
statistics for the difference between the original and SP3-corrected residuals are summarized
in Table 4.1
Figure 4.1: Comparison of original and SP3-corrected residuals for satellite 14
Figure 4.2: Comparison of original and SP3-corrected residuals for satellite 31
22
Table 4.1: Summary statistics of differences between original and SP3-corrected residuals
Satellite Channel Mean Standard Deviation
1 2.44× 10?4 6.54× 10?4
14
2 2.40× 10?4 7.15× 10?4
1 ?5.05× 10?4 4.08× 10?4
31
2 ?5.03× 10?4 4.41× 10?4
The two channels (L1, L2) in each satellite had approximately equivalent performance,
as summarized in Table 4.1. This makes logical sense, since the signals were following almost
the exact same path from satellite to receiver and were thus subject to the same conditions
along the way. As seen in Figure 4.2, both channels in satellite 31 had a growing bias
between the original measurements and the SP3 corrections. The increasing bias in satellite
31 resulted in a mean value (?5.05 × 10?4 for channel 1 and ?5.03 × 10?4 for channel 2)
that has a greater magnitude than that of satellite 14.
4.2 Autocorrelation
Figure 4.3 illustrates the autocorrelation for the residuals from satellite 14. There was a
Figure 4.3: Normalized autocorrelation for satellite 14
peak that did not resemble an impulse function (i.e., there was not a single spike), indicating
23
that the residuals were not independent over time. The primary peak in Figure 4.3 is
highlighted in Figure 4.4.
Figure 4.4: Normalized autocorrelation for satellite 14 with primary peak highlighted
The base of the peak was approximately 1800 s indicated there was a statistical de-
pendence in points with a time constant on the scale of several hundred seconds . A time
constant on that scale suggested that the dependence could be related to motion, position,
or orientation (e.g., aircraft flying in a loop that had a 30-min orbit). We hypothesize that
a signal-related autocorrelation would have had a time constant the order of seconds or even
sub-second.
Figure 4.5 illustrates the autocorrelation for the residuals from satellite 31. The primary
and secondary peaks in Figure 4.5 are highlighted in Figure 4.6.
24
Figure 4.5: Normalized autocorrelation for satellite 31
Figure 4.6: Normalized autocorrelation for satellite 31 with primary (red) and secondary
(blue) peaks highlighted
25
The primary base of the peak of approximately 500 s indicated there was a statistical
dependence in points with a time constant on the scale of several hundred seconds. A time
constant on that scale suggested that the dependence could be related to motion, position,
or orientation (e.g., aircraft flying in a loop that had an 8-min orbit). We hypothesize that
a signal-related autocorrelation would have had a time constant on the order of seconds or
even sub-seconds. The secondary base of the peak was approximately 2000 s (or 33 min),
which is similar to the time constant from satellite 14 of 1800 s (or 30 min). It is likely that
both time constants were identifying the same phenomenon.
The overall shape of the autocorrelation of both satellites was approximately triangu-
lar, which suggested the distribution approximated a Gauss-Markov process. They may be
related to residual errors in the estimated solution or satellite geometry.
4.3 Allan Variance
The relationship between Allan variance (?AV ) and the cluster length (?) for satellite
14 is illustrated in Figure 4.7. The slope of the best fit line describes the trend, based on
Table 2.2. The first segment had a slope of -0.207 and the second segment had a slope of
0.603. The first slope indicated the presence of a correlated noise process, while the second
slope suggested either white noise or a time-differenced random walk trend. However, based
on the autocorrelation results, there was not pure white noise present. Additionally, since
the initial residuals were time-differenced, the time-differenced random walk was more likely
to be process that was present.
The Allan variance had a minimum value around a time constant of ? = 200 s. There-
fore, to minimize the impact of stochastic processes, the residuals should be averaged at
approximately 200 s intervals.19
There was a distinct pattern in the Allan variance results around a time constant of
26
Figure 4.7: Allan variance analysis for satellite 14
? = 103 s = 1000 s, which was on the same order of magnitude as the 1800 s halfwidth of
the autocorrelation peak of satellite 14. Therefore, the Allan variance analysis was likely
identifying the same effects from the autocorrelation analysis illustrated in Figure 4.4.
The relationship between Allan variance (?AV ) and the cluster length (?) for satellite 31
is illustrated in Figure 4.8.
The first segment had a slope of -0.218 and the second segment had a slope of 1.057.
The first slope indicates the presence of a correlated noise process, similar to the results
of satellite 14 in Figure 4.7. The second slope also suggests either white noise or a time-
differenced random walk trend.
The Allan variance had a minimum value around a time constant of ? = 400 s. To min-
imize the impact of stochastic processes, the residuals should be averaged at approximately
400 s intervals.19Overall, the Allan variance analysis for satellite 31 agreed with the results
of satellite 14.
27
Figure 4.8: Allan variance analysis for satellite 31
4.4 Elevation Angle
The carrier-to-noise ratio C/N0 is affected by the elevation angle between the satellite
and the user, since the signal interacts with more of the atmosphere at low elevation angles of
15? or lower.29 The elevation angle for the duration of the flight for satellite 14 is illustrated
in Figure 4.9.
Figure 4.9: Elevation angle of satellite 14
28
The elevation angle fell below 15? about 5000 s (83 min) into the flight. Since both
channels in satellite 14 had similar results, channel 1 aws used for further analysis of the effect
of elevation angle. The corrected residuals from channel 1 from Figure 4.1 are reproduced
in Figure 4.10 with the addition of a 250-point centered moving average. Figure 4.10 with
an approximation of the moving average line is illustrated in Figure 4.11.
Figure 4.10: Corrected residuals from satellite 14, channel 1 with a 250-point centered
moving average
Figure 4.11: Corrected residuals from satellite 14, channel 1 with a 250-point centered
moving average and approximation of the moving average
There was a noticeable change in the slope around the 5300 s mark, resulting in a
29
first-order bias with slope of approximately ?2 × 10?6 m/s . Prior to the 5300 s mark, the
s
moving average appeared to identify a zeroth-order bias (i.e., approximately constant value)
of ?6.5 × 10?4, though there was some observable sinusoidal noise. The simultaneity of
the low elevation angle and change in slope of the residuals suggests there may have been
a causal relationship between the two. In an operational case, users would likely disregard
signals from satellites with elevation angles below 15?.
The elevation angle for the duration of the flight for satellite 31 is illustrated in Fig-
ure 4.12. The elevation angle falls below 15? about 7500 s (125 min) into the flight. The
Figure 4.12: Elevation angle of satellite 31
7500 s point also marks the end of the flight, so the low elevation angle did not appear to
significantly affect the overall results for satellite 31.
Since both channels in satellite 31 had similar results, channel 1 was used for further anal-
ysis of the effect of elevation angle. The corrected residuals from channel 1 from Figure 4.2
were reproduced in Figure 4.13 with the addition of a 250-point centered moving average.
Figure 4.13 with an approximation of the moving average line is illustrated in Figure 4.14.
30
Figure 4.13: Corrected residuals from satellite 31, channel 1 with a 250-point centered
moving average
Figure 4.14: Corrected residuals from satellite 31, channel 1 with a 250-point centered
moving average and approximation of the moving average
31
For the first 2000 s, there is a noticeable drop dip the residuals. The decrease and increase
are both roughly first-order biases with a slope of approximately ??2.5×10?4, respectively.
After to the 2000 s mark, the moving average appears to identify a zeroth-order bias (i.e.,
approximately constant value) of ?9.2 × 10?4, though there is some observable sinusoidal
noise. Unlike the case with satellite 14, there was sinusoidal noise even when the elevation
angle was above 15?.
4.5 Autocorrelation of Truncated Data
For a more representative analysis, the residuals from satellite 14 were truncated to 5300 s
(i.e., data from the satellite when it was below an elevation angle of 15? were discarded). Fig-
ure 4.15 illustrates the recalculated autocorrelation for the truncated residuals from satellite
14. For comparison, both the original and truncated autocorrelation analyses are illustrated
Figure 4.15: Normalized autocorrelation for truncated satellite 14 residuals
in Figure 4.16.
As seen in Figure 4.16, there was not a distinct wide peak in the truncated autocorrelation
like there was in the original data. The overall triangle distribution of the autocorrelation
32
Figure 4.16: Comparison of original residuals and truncated residuals for satellite 14
in Figure 4.15 suggested the truncated residuals better fits a Gaussian distribution, but it
still includes a correlation competent.
For the sake of thoroughness, the residuals for satellite 31 were truncated to 6250 s, which
was when the satellite fell below 20?. Figure 4.17 illustrates the recalculated autocorrelation
for the truncated residuals from satellite 31. For comparison, both the original and truncated
Figure 4.17: Normalized autocorrelation for truncated satellite 31 residuals
autocorrelation analyses are illustrated in Figure 4.18.
As seen in Figure 4.18, there was not a significant difference between the overall shapes
of the autocorrelation analysis for the original data and truncated data. The similarities in
33
Figure 4.18: Comparison of original residuals and truncated residuals for satellite 31
shape verifies that there were no significant impacts on the autocorrelation caused by low
elevation angle for satellite 31.
4.6 Allan Variance Analysis of Truncated Data
Similar to the autocorrelation analysis, the Allan variance analysis was re-accomplished
for satellite 14 truncated to 5300 s. The relationship between Allan variance (?AV ) and the
cluster length (?) for the truncated residuals from satellite 14 is illustrated in Figure 4.19.
The first segment had a slope of -0.285 and the second segment had a slope of 1.896.
The first slope indicated the presence of a Gauss-Markov process. The second slope did not
identify any stochastic processes. The Allan variance had a minimum value around a time
constant of ? = 300 s. To minimize the impact of stochastic processes, the residuals should
be averaged at approximately 300 s intervals.19
Notably, the distinct pattern that was present in Figure 4.7 around ? = 103 s was missing
in Figure 4.19. Therefore, the low elevation angle caused the distinct pattern. Likewise, the
autocorrelation results in Figure 4.15 were also missing the identified spike from Figure 4.3,
validating the previous conclusion that the Allan variance analysis was identifying the same
34
Figure 4.19: Allan variance analysis for truncated residuals from satellite 14
phenomenon as as the autocorrelation analysis.
The autocorrelation results in Figure 4.15 suggested the truncated residuals included ran-
dom Gaussian noise. The Allan variance analysis suggested the presence of a correlated noise
process. The truncated residuals may have a mix of Gaussian noise and a Gauss-Markov
process. However, since the time constant was large, the stochastic process could also be
modeled as a random walk process. Time differencing each successive data point would lead
to a Gaussian distribution. If the differencing occurred at a different interval, the result
would be a Gauss-Markov process. Either way, both distributions are Gaussian in nature,
so the assumption that the residuals were independent over time was valid.
Since the measurements had been time-differenced, we did not expect to see time-
correlated noise of measurements. Additionally, errors caused by tropospheric effects were
minimized in the truncated data set. It was more likely that the majority of the errors were
related to aircraft motion, position, orientation, and solution. Errors in ?U? (from Equa-
tion 3.7) could be correlated over time, which could be responsible for the correlated noise
35
components.
The relationship between Allan variance (?AV ) and the cluster length (?) for the truncated
residuals from satellite 31 is illustrated in Figure 4.20. The first segment had a slope of -0.218
Figure 4.20: Allan variance analysis for truncated residuals from satellite 31
and the second segment had a slope of 0.559. The first slope matched the slope from the
Allan variance analysis for the original data and also suggested the presence of a correlated
noise process. The second slope suggested the presence of a time-differenced random walk
process, similar to the original analysis. The Allan variance had a minimum value around
a time constant of ? = 500 s, versus ? = 400 s for the original analysis. To minimize
the impact of stochastic processes, the residuals should be averaged at approximately 500 s
intervals.19
Overall, the analysis for the truncated data for satellite 31 matched the original analysis.
Additionally, the identified processes in satellite 31 agreed with the results from satellite 14.
36
4.7 Downsampling and Thresholding
Residual errors were downsampled to reduce teh impact of error correlation. The SP3-
corrected residuals for satellite 14 were downsampled at various rates (“window sizes”). The
mean and RMS deviation of the first “initialization period” windows were used to establish
thresholds to flag subsequent windows as potential outliers. Note that these outliers are not
necessarily erroneous; rather, they are outlier points that may merit further investigation.
Points were flagged as outliers if they were more than three times the RMS from the mean.
The different combinations of window size and number of initialization periods are illustrated
in Figure 4.21 and Figure 4.22, with the outliers flagged with red circles.
37
Figure 4.21: Variation of initialization and window sizes for satellite 14, Part A
38
Figure 4.22: Variation of initialization and window sizes for satellite 14, Part B
39
Since each combinations of initialization period and window size had a unique number of
data points, counting just the number of flagged outliers was not a representative measure-
ment. Instead, the percentage of total points flagged as outliers was compared to measure
effectiveness, as illustrated in Figure 4.23 and Figure 4.24.
Figure 4.23: Percentage of points identified as outliers versus initialization period
Figure 4.24: Percentage of points identified as outliers versus window size
As seen in Figure 4.23, the percentage of identified outliers increased with a larger initial-
ization period. From the Allan variance analysis, one of the key window sizes identified was
250 s. There were no outliers identified with a window size of 250 s with an initialization
40
period of 25 windows, since almost all of the points fell within the initialization period. As
seen in Figure 4.24, there was not significant variation over the different window sizes for
each of the different selected initialization periods. The combination of a window size of
25 s and an initialization period of 25 windows appeared to identify the most outliers. This
translated to a total of 625 s or 10.4 min of total initialization time.
The different combinations of window size and number of initialization periods for satellite
31 are illustrated in Figure 4.25 and Figure 4.26, with the outliers flagged with red circles.
41
Figure 4.25: Variation of initialization and window sizes for satellite 31, Part A
42
Figure 4.26: Variation of initialization and window sizes for satellite 31, Part B
43
SThe percentage of total points flagged as outliers was compared to measure effectiveness,
as illustrated in Figure 4.27 and Figure 4.28.
Figure 4.27: Percentage of points identified as outliers versus initialization period
Figure 4.28: Percentage of points identified as outliers versus window size
As seen in Figure 4.27, the percentage of identified outliers increased with a larger ini-
tialization period. Like with satellite 14, there were no outliers identified when the case of
a window size of 250 s with an initialization period of 25 windows since almost all of the
points fell within the initialization period. As seen in Figure 4.28, there was not significant
44
variation over the different window sizes for each of the different selected initialization peri-
ods. The combination of a window size of 25 s and an initialization period of 25 windows
appeared to identify the most outliers. This translated to a total of 625 s or 10.4 min of
total initialization time.
The results in Figure 4.27 and Figure 4.28 agree with the results for satellite 14 from
Figure 4.23 and Figure 4.24.
4.8 Simulated Bias in Satellite Orbit
The difference between the original and SP3-corrected residuals is illustrated in the bot-
tom two plots of both Figure 4.1 and Figure 4.2. There was a slight bias evident in Figure 4.2;
however, for a more controlled test, the simulated first-order biases of different levels were
added to the the existing data for satellite 14 to simulate drift. The following slopes were
used: 1× 10?4, 5× 10?5, 1× 10?5, 5× 10?6, 1× 10?6, 5× 10?7, and 1× 10?7.
A subset of the combinations of bias slope, window size, and initialization period are
illustrated in Figure 4.29 and Figure 4.30. Note that the size of the window is equal to the
time of each cycle (e.g., a window size of 10 means each averaging cycle lasts 10 s). The time
to identify each simulated bias for each combination of initialization period and window size
is summarized in Table 4.2 through Table 4.8.
45
Figure 4.29: Variation of added bias, initialization period, and window size for residuals for
satellite 14, Part A
46
Figure 4.30: Variation of added bias, initialization period, and window size for residuals for
satellite 14, Part B
47
Table 4.2: Time (in s) to identify simulated bias of 1× 10?4
Initialization Window Size (s)
Period (s) 10 25 50 60 75 100 150 250
3 50 75 50 60 75 100 150 250
5 50 50 50 60 75 100 150 250
10 30 50 50 60 75 100 150 250
15 30 50 50 60 75 100 150 250
20 30 50 50 60 75 100 150 250
25 30 25 50 60 75 100 150 250
Table 4.3: Time (in s) to identify simulated bias of 5× 10?5
Initialization Window Size (s)
Period (s) 10 25 50 60 75 100 150 250
3 90 100 100 120 150 100 150 250
5 90 75 100 120 150 100 150 250
10 60 50 100 120 75 100 150 250
15 30 50 100 120 75 100 150 250
20 30 50 100 60 75 100 150 250
25 30 50 50 60 75 100 150 250
Table 4.4: Time (in s) to identify simulated bias of 1× 10?5
Initialization Window Size (s)
Period (s) 10 25 50 60 75 100 150 250
3 270 300 250 240 300 200 300 500
5 270 200 150 180 225 200 300 250
10 140 125 150 180 150 200 300 250
15 110 125 150 180 150 200 150 250
20 110 125 100 120 150 100 150 250
25 100 50 50 60 75 100 150 250
48
Table 4.5: Time (in s) to identify simulated bias of 5× 10?6
Initialization Window Size (s)
Period (s) 10 25 50 60 75 100 150 250
3 480 650 300 300 300 300 450 500
5 480 325 250 240 300 300 300 500
10 480 250 250 180 300 300 300 250
15 290 200 250 180 225 200 150 250
20 230 125 100 180 150 100 150 250
25 200 50 50 60 75 100 150 250
Table 4.6: Time (in s) to identify simulated bias of 1× 10?6
Initialization Window Size (s)
Period (s) 10 25 50 60 75 100 150 250
3 1760 2925 1750 1440 1800 1700 1350 1500
5 1670 1200 1250 1260 1275 1200 1200 1500
10 1460 875 1150 1200 1200 1100 450 500
15 1460 875 1100 1140 225 200 450 250
20 1080 875 950 780 225 200 150 250
25 1070 850 750 600 150 200 150 250
Table 4.7: Time (in s) to identify simulated bias of 5× 10?7
Initialization Window Size (s)
Period (s) 10 25 50 60 75 100 150 250
3 4380 4650 3700 3660 3525 2400 4500 4000
5 4160 4050 2600 1440 1500 1400 1800 1500
10 4100 3525 1750 1380 1350 1400 1350 500
15 3900 3375 1350 1380 1350 1300 600 500
20 3860 3250 1250 1260 1200 800 300 500
25 3260 3200 300 1200 825 400 300 250
Table 4.8: Time (in s) to identify simulated bias of 1× 10?7
Initialization Window Size (s)
Period (s) 10 25 50 60 75 100 150 250
3 Not identified Not identified 4650 4200 4575 5400 5100 5250
5 Not identified 4700 4650 4020 4050 5000 4800 4500
10 4760 4600 4150 3960 3975 4500 4050 3250
15 4690 4125 4080 3960 3975 4000 3300 2000
20 4600 4125 3900 3900 3975 3200 2550 750
25 4550 4125 3900 3660 3600 2700 1800 250
49
The thresholds were able to flag biased data points with varying levels of effectiveness.
For the larger bias slopes (1×10?4, 5×10?5, and 1×10?5), the biases were identified quickly
(less than 500 s, or 8.3 min). For the larger bias slopes, the shorter initialization periods
and smaller window sizes identified the bias faster than the larger initialization periods and
larger window sizes. Typically, the larger initialization periods and window sizes were able to
identify the bias in one or two cycles; however, each of those cycles was longer than the few
cycles needed for the shorter initialization periods and smaller window sizes. Conversely, for
the smaller bias slopes (1×10?6, 5×10?7, 1×10?7), the combination of longer initialization
periods and larger window sizes was faster at identifying the bias. For the bias slope of
5× 10?6 all combinations of initialization period and window size were similarly effective at
identifying the simulated bias.
The moving average of the residuals from Figure 4.11 identified a first-order bias with
a slope of ?2 × 10?6 m/s . In most of the cases illustrated in Figure 4.21 and Figure 4.22,
s
the deviation was identified. For a bias with a slope of ?2 × 10?6 m/s , most combinations
s
of initialization period and window size would effectively identify the bias. The longer ini-
tialization periods and larger window sizes may slightly outperform the shorter initialization
periods and smaller window sizes.
A subset of the combinations of bias slope, window size, and initialization period for
satellite 31 are illustrated in Figure 4.31 and Figure 4.32. The downsampling and threshold-
ing results for satellite 31 match those of satellite 14. Like with satellite 14, the larger bias
slopes were more quickly identified by the shorter initialization periods and smaller window
sizes. The smaller bias slopes were more quickly identified by the longer initialization periods
and larger window sizes.
50
Figure 4.31: Variation of added bias, initialization period, and window size for residuals for
satellite 31, Part A
51
Figure 4.32: Variation of added bias, initialization period, and window size for residuals for
satellite 31, Part B
52
Chapter 5
Conclusion
5.1 Detect and Identify Error Sources
The results and analysis for satellite 14 and satellite 31 were presented in Section 4.
The results and analysis for satellite 31 matched those from satellite 14, demonstrating
repeatability of the methodology. The residuals from RANSAC and TDCP were able to
detect and identify a variety of error sources. The downsampling and thresholding method
was able to identify zeroth-order deviations. Clusters of deviations were able to identify
first-order deviations with slopes on the orders of magnitude of 10?4 down to 10?7 after
seconds (for larger bias slopes) or hours (for smaller bias slopes). The autocorrelation and
Allan variance analysis were not good indicators of initialization period and window size.
Instead, a user should select combinations of initialization periods and window sizes based
on the magnitude of deviations they are looking to discover. Additionally, detection should
not be limited to just one solution methodology; instead, this method should be used in
conjunction with other detection methods.
The multipath errors manifested as random noise. The correlated noise likely resulted
from a combination of the ionospheric delays, tropospheric delays, and errors in the esti-
mation of the velocity solution. Additionally, inter-system timing biases could have lead to
further deviations. The zeroth-order and larger deviations could have been caused by errors
in satellite orbit prediction and multipath. The first-order deviations and smaller deviations
could have been caused by ionospheric and tropospheric delays.
Autocorrelation analysis and Allan variance analysis identified other underlying proper-
ties of the data, such as an autocorrelated feature with a time constant around 1000 s to
1800 s. However, after residuals from low satellite angles were removed, the autocorrelation
suggested the residuals were randomly distributed. The Allan variance analysis of the trun-
cated data suggested the presence of a correlated noise process. Therefore, the assumption
that the residuals were independent with time was valid.
5.2 Recommendations for Future Work
Once multi-GNSS navigation is feasible, the next step would be to compare it to other
sources of PNT data. Different uses of PNT information have different requirements for
precision, so PNT information does not necessarily have to come from a high-accuracy timing
signal. Instead of a GPS signal, waves from other parts of the spectrum (e.g., acoustic,
ultrasonic, visible light, radio signals, magnetic waves, cell phone signals) can be used to
determine location using time of arrival (TOA), angle of arrival (AOA), time difference
of arrival (TDOA), or time of flight (TOF) calculations.4,12,13 Pedestrian dead reckoning
uses an inertial navigation unit (INU) to keep track of changes in position from a known
starting point.4,13 The signals from from multiple known Wi-Fi routers can use relative signal
strength (RSS) to determine location.4,13,39
If the area is well mapped, then these signals of opportunity can be correlated to the
map to further increase PNT confidence anywhere from 45% to 75%.37 All of these indicators
can be used individually or in parallel to increase the confidence in all-source navigation.18
54
Instead of comparing different GNSS constellations to each other, different PNT technologies
can be compared to each other. See Tariq, et al.4 for a more detailed breakdown of each of
these technologies.
The operational effectiveness of implementing the multi-GNSS navigation corrections
should be studied in realistic settings. While ground and flight tests in operationally rel-
evant environments would be preferred, the use of modeling and simulation tools may be
a sufficient alternative. The Air Force Research Lab’s Sensors Directorate (AFRL/RY),
Aerospace Systems Directorate (AFRL/RQ), or Space Vehicles Directorate (AFRL/RV) may
be well-positioned to undertake such an analysis.
Lastly, the methods described here can be applied to other areas for error detection. Most
time-series signals can be monitored for consistency. Other approaches rely on machine-
learning algorithms to identify errors in data sets.21
55
References
[1] BA Kumar Babu. Method and system for resolving double difference GPS carrier phase
integer ambiguity utilizing decentralized Kalman filters, Sep 1995. US Patent 5,451,964.
[2] Gerhard Beutler, Jan Kouba, and Tim Springer. Combining the orbits of the IGS
analysis centers. Bulletin Ge?ode?sique, 69(4):200–222, December 1995. doi: 10.1007/
bf00806733. URL https://doi.org/10.1007/bf00806733.
[3] Gerhard Beutler, Rolf Dach, Oliver Montenbruck, Urs Hugentobler, Georg Weber, and
Elmar Brockmann. GLONASS and multi-GNSS in the IGS: Lessons learned from
GLONASS service disruptions. In 13th Meeting of the National Space-Based Posi-
tioning, Navigation, and Timing (PNT) Advisory Board. International Association of
Geodesy, 2014. doi: 10.7892/BORIS.53907. URL https://boris.unibe.ch/53907/.
[4] Zain Bin Tariq, Dost Muhammad Cheema, Muhammad Zahir Kamran, and Ijaz
Haider Naqvi. Non-GPS positioning systems. ACM Computing Surveys, 50(4):1–34,
Aug 2017. doi: 10.1145/3098207. URL https://doi.org/10.1145/3098207.
[5] Andreas Brack. Partial Carrier-Phase Integer Ambiguity Resolution for High Accuracy
GNSS Positioning. Dissertation, Technische Universita?t Mu?nchen, Mu?nchen, 2019.
[6] Piet M.T. Broersen. Automatic Autocorrelation and Spectral Analysis. Springer London,
2006. ISBN 9781846283291.
[7] George W Bush. U.S. space-based position, navigation, and timing policy. National
Security Policy Directive, 39, December 2004. URL https://fas.org/irp/offdocs/
nspd/nspd-39.htm.
[8] Gaetano Castaldo, Antonio Angrisano, Salvatore Gaglione, and Salvatore Troisi. P-
ransac: an integrity monitoring approach for GNSS signal degraded scenario. Interna-
tional Journal of Navigation and Observation, 2014, 2014.
[9] Varun Chandola, Arindam Banerjee, and Vipin Kumar. Anomaly detection. ACM
Computing Surveys, 41(3):1–58, July 2009. doi: 10.1145/1541880.1541882. URL https:
//doi.org/10.1145/1541880.1541882.
[10] William J Clinton. Statement by the President regarding the United States’ decision to
stop degrading global positioning system accuracy. In Office of Science and Technology
Policy, Washington, DC, 2000. Office of the Press Secretary.
[11] John C Cohenour. Global positioning system clock and orbit statistics and precise point
positioning. PhD thesis, Ohio University, 2009.
[12] Azusa Danjo, Yuki Watase, and Shinsuke Hara. A theoretical error analysis on in-
door TOA localization scheme using unmanned aerial vehicles. In IEEE 81st Vehicular
Technology Conference, pages 1–5. IEEE, 2015.
[13] Davide Dardari, Pau Closas, and Petar M Djuric?. Indoor tracking: Theory, methods,
and technologies. IEEE Transactions on Vehicular Technology, 64(4):1263–1278, 2015.
[14] Thomas Dautermann, Christoph Mayer, Felix Antreich, Andriy Konovaltsev, Boubeker
Belabbas, and Ulrich Kalberer. Non-gaussian error modeling for gbas integrity assess-
ment. IEEE transactions on aerospace and electronic systems, 48(1):693–706, 2012.
[15] Naser El-Sheimy, Haiying Hou, and Xiaoji Niu. Analysis and modeling of inertial sensors
using allan variance. IEEE Transactions on Instrumentation and Measurement, 57(1):
140–149, January 2008. doi: 10.1109/tim.2007.908635. URL https://doi.org/10.
1109/tim.2007.908635.
[16] Federal Aviation Administration. GNSS frequently asked questions, Dec
2016. URL https://www.faa.gov/about/office_org/headquarters_offices/ato/
service_units/techops/navservices/gnss/faq/gps/.
[17] Martin A Fischler and Robert C Bolles. Random sample consensus: a paradigm for
model fitting with applications to image analysis and automated cartography. Commu-
nications of the ACM, 24(6):381–395, 1981.
[18] Kenneth A Fisher and John F Raquet. Precision position, navigation, and timing
without the global positioning system. Air & Space Power Journal, pages 24–33, 2011.
[19] Thomas Friederichs. Analysis of geodetic time series using Allan variances. PhD thesis,
Universita?t Stuttgart, 2010. URL http://elib.uni-stuttgart.de/handle/11682/
3866.
[20] JG Garc??a, PI Mercader, and Carlos H Muravchik. Use of GPS carrier phase double
differences. Latin American applied research, 35(2):115–120, 2005.
[21] Alexander Geiger, Dongyu Liu, Sarah Alnegheimish, Alfredo Cuesta-Infante, and
Kalyan Veeramachaneni. TadGAN: Time series anomaly detection using generative
adversarial networks. arXiv preprint arXiv:2009.07769, 2020.
57
[22] Markus Goldstein and Seiichi Uchida. A comparative evaluation of unsupervised
anomaly detection algorithms for multivariate data. PLOS ONE, 11(4), April 2016.
doi: 10.1371/journal.pone.0152173. URL https://doi.org/10.1371/journal.pone.
0152173.
[23] GPS.gov. Timing applications, Nov 2019. URL https://www.gps.gov/applications/
timing/.
[24] GPS.gov. Other global navigation satellite systems (GNSS), Jan 2020. URL https:
//www.gps.gov/systems/gnss/.
[25] R.M. Gray and D.L. Neuhoff. Quantization. IEEE Transactions on Information Theory,
44(6):2325–2383, 1998. doi: 10.1109/18.720541. URL https://doi.org/10.1109/18.
720541.
[26] Jake Griffiths and Jim R. Ray. On the precision and accuracy of IGS orbits. Journal
of Geodesy, 83(3-4):277–287, February 2009. doi: 10.1007/s00190-008-0237-6. URL
https://doi.org/10.1007/s00190-008-0237-6.
[27] Steve Hilla. The extended standard product 3 orbit format (SP3-c). Technical report,
National Geodetic Survey, National Oceanic and Atmospheric Administration, 2010.
URL https://files.igs.org/pub/data/format/sp3c.txt.
[28] David C Jefferson and Yoaz E Bar-Sever. Accuracy and consistency of broadcast GPS
ephemeris data. In Proceedings of the 13th International Technical Meeting of the Satel-
lite Division of The Institute of Navigation (ION GPS 2000), pages 391–396. Institute
of Navigation, 2000.
[29] Elliott D Kaplan and Christopher J Hegarty. Understanding GPS/GNSS: Principles
and Applications. Artech House, third edition, 2017.
[30] Samer Khanafseh, Birendra Kujur, Mathieu Joerger, Todd Walter, Sam Pullen, Juan
Blanch, Kevin Doherty, Laura Norman, Lance de Groot, and Boris Pervan. GNSS
multipath error modeling for automotive applications. In Proceedings of the 31st Inter-
national Technical Meeting of The Satellite Division of the Institute of Navigation (ION
GNSS 2018). Institute of Navigation, October 2018. doi: 10.33012/2018.16107. URL
https://doi.org/10.33012/2018.16107.
[31] Tomislav Kos, Ivan Markezic, and Josip Pokrajcic. Effects of multipath reception on
GPS positioning performance. In Proceedings ELMAR-2010, pages 399–402. IEEE,
2010.
[32] Philip Kwan. NAVSTAR GPS Space Segment/Navigation User Segment Interfaces, k
edition, May 2019. IS-GPS-200K.
[33] Philip Kwan. NAVSTAR GPS Space Segment/Navigation User Segment L5 Interfaces,
f edition, May 2019. IS-GPS-705F.
58
[34] Bofeng Li, Zhetao Zhang, Yunzhong Shen, and Ling Yang. A procedure for the signif-
icance testing of unmodeled errors in GNSS observations. Journal of Geodesy, 92(10):
1171–1186, January 2018. doi: 10.1007/s00190-018-1111-9. URL https://doi.org/
10.1007/s00190-018-1111-9.
[35] Jinyu Li, Li Deng, Reinhold Haeb-Umbach, and Yifan Gong. Fundamentals of speech
recognition. Robust Automatic Speech Recognition, pages 9–40, 2016.
[36] Yiran Luo, Jian Li, Chunyang Yu, Bing Xu, You Li, Li-Ta Hsu, and Naser El-Sheimy.
Research on time-correlated errors using allan variance in a kalman filter applicable to
vector-tracking-based GNSS software-defined receiver for autonomous ground vehicle
navigation. Remote Sensing, 11(9):1026, April 2019. doi: 10.3390/rs11091026. URL
https://doi.org/10.3390/rs11091026.
[37] Mahdi Maaref and Zaher M Kassas. Ground vehicle navigation in GNSS-challenged
environments using signals of opportunity and a closed-loop map-matching approach.
Transactions on Intelligent Transportation Systems, 2019.
[38] VB Mendes and RB Langley. Tropospheric zenith delay prediction accuracy for airborne
gps high-precision positioning. In Proceedings of the Institute of Navigation 54th Annual
Meeting. Institute of Navigation, June 1998. URL https://doi.org/10.33012/2020.
17163.
[39] Ana Moragrega, Pau Closas, and Christian Ibars. Potential game for energy-efficient
RSS-based positioning in wireless sensor networks. IEEE Journal on Selected Areas in
Communications, 33(7):1394–1406, 2015.
[40] National Aeronautics and Space Administration. GNSS orbit products, Mar
2020. URL https://cddis.nasa.gov/Data_and_Derived_Products/GNSS/orbit_
products.html.
[41] Xiaoji Niu, Qijin Chen, Quan Zhang, Hongping Zhang, Jieming Niu, Kejie Chen,
Chuang Shi, and Jingnan Liu. Using Allan variance to analyze the error char-
acteristics of GNSS positioning. GPS solutions, 18(2):231–242, April 2014. doi:
10.1007/s10291-013-0324-x. URL https://doi.org/10.1007/s10291-013-0324-x.
[42] XJ Niu, Q Zhang, and HP Zhang. Use Allan variance to analyze the results of GNSS/INS
navigation systems. In Proceedings of the CPGPS 2010 Technical Forum on Satellite
Navigation and Positioning, Shanghai, China, pages 565–572, 2010.
[43] United States Department of Defense. National Defense Strategy of the United States
of America. United States Department of Defense, 2018.
[44] United States Department of Defense. Global Positioning System Standard Positioning
Service Performance Standard, fifth edition, April 2020.
59
[45] United States Chairman of the Joint Chiefs of Staff. National Military Strategy. United
States Chairman of the Joint Chiefs of Staff, 2018.
[46] United States White House Office. National Security Strategy of the United States of
America. United States White House Office, 2017.
[47] Ahmed Radi, Sameh Nassar, Maan Khedr, Naser El-Sheimy, Roberto Molinari, and
Ste?phane Guerrier. Improved stochastic modelling of low-cost GNSS receivers position-
ing errors. In 2018 IEEE/ION Position, Location and Navigation Symposium (PLANS),
pages 108–117. IEEE, 2018.
[48] Benjamin W Remondi. Extending the national geodetic survey standard GPS orbit
formulas. Technical Report NOS 133 NGS 46, National Oceanic and Atmospheric Ad-
ministration, 1989.
[49] Georg Schroth, Markus Rippl, Alexandru Ene, Juan Blanch, Boubeker Belabbas, Todd
Walter, Per Enge, and Michael Meurer. Enhancements of the range consensus algorithm
(RANCO). In Proceedings of the ION GNSS Conference, pages 93–103, 2008.
[50] Nihar N. Shah, Zhen Zhu, Eric Vinande, and Jason Pontious. The feasibility of ob-
serving user-segment faults in the residuals of a consensus-based multi-GNSS solu-
tion. In Proceedings of the 2020 International Technical Meeting of The Institute of
Navigation. Institute of Navigation, February 2020. doi: 10.33012/2020.17163. URL
https://doi.org/10.33012/2020.17163.
[51] Global Navigation Satellite Systems – What’s up?, June 2019. Space-Based Position,
Navigation & Timing National Advisory Board.
[52] J. J. Spilker. GPS signal structure and performance characteristics. Navigation, 25(2):
121–146, June 1978. doi: 10.1002/j.2161-4296.1978.tb01325.x. URL https://doi.org/
10.1002/j.2161-4296.1978.tb01325.x.
[53] Michael L. Stein. Interpolation of Spatial Data. Springer, 12 2012. ISBN 9781461214946.
[54] GDAL Team. Satellite of GDAL, 2010. Licensed under CC BY-SA.
[55] Charles Therrien and Murali Tummala. Probability and Random Processes for Electrical
and Computer Engineers. CRC Press LLC, 09 2018. ISBN 9781439897492.
[56] Bryan Townsend and Patrick Fenton. A practical approach to the reduction of pseu-
dorange multipath errors in a l1 gps receiver. In Proceedings of the 7th International
Technical Meeting of the Satellite Division of the Institute of Navigation, Salt Lake City,
UT, USA. Institute of Navigation, 1994.
[57] James B Tsui. Fundamentals of Global Positioning System Receivers: A Software Ap-
proach. John Wiley & Sons, 2 edition, 2005.
60
[58] A. Tuka and A. El-Mowafy. Performance evaluation of different troposphere delay
models and mapping functions. Measurement, 46(2):928–937, February 2013. doi: 10.
1016/j.measurement.2012.10.015. URL https://doi.org/10.1016/j.measurement.
2012.10.015.
[59] Frank van Graas and Andrey Soloviev. Precise velocity estimation using a stand-alone
GPS receiver. Navigation, 51(4):283–292, 2004.
[60] Marc Weiss, Pradipta Shome, and Ron Beard. On-board GPS clock monitoring for
signal integrity. In Proceedings of the 42nd Annual Precise Time and Time Interval
Systems and Applications Meeting, pages 465–480. National Institute of Standards and
Technology, 2010.
[61] Asmamaw Yehun, Tsegaye Kassa, Martin Vermeer, and Addisu Hunegnaw. Higher
order ionospheric delay and derivation of regional total electron content over ethiopian
global positioning system stations. Advances in Space Research, May 2020. doi: 10.
1016/j.asr.2020.04.035. URL https://doi.org/10.1016/j.asr.2020.04.035.
[62] W Yu, G Lachapelle, and S Skone. Pll performance for signals in the presence of thermal
noise, phase noise, and ionospheric scintillation. In Proceedings of the 19th International
Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS
2006), pages 1341–1357, 2006.
[63] Eva Z?a?c?ekova?, Samuel Pr??vara, Zdene?k Va?n?a, and Ji?r?? Cigler. Building modeling: on
selection of the model complexity for predictive control. In Computer Aided Chemical
Engineering, volume 32, pages 205–210. Elsevier, 2013.
[64] Chunxi Zhang, Xianmu Li, Shuang Gao, Tie Lin, and Lu Wang. Performance analysis
of global navigation satellite system signal acquisition aided by different grade inertial
navigation system under highly dynamic conditions. Sensors, 17(5):980, April 2017.
doi: 10.3390/s17050980. URL https://doi.org/10.3390/s17050980.
[65] Xin Zhang, Yong Li, Peter Mumford, and Chris Rizos. Allan variance analysis on
error characters of MEMS inertial sensors for an FPGA-based GPS/INS system. In
Proceedings of the International Symposium on GPS/GNSS, pages 127–133, 2008.
[66] Zhen Zhu, Eric Vinande, and Maarten Uijt de Haag. Multi-constellation time-
differenced carrier phase solution with protection from multiple failures. In 2018
IEEE/ION Position, Location and Navigation Symposium (PLANS), pages 336–348.
IEEE, 2018.
[67] Zhen Zhu, Eric Vinande, Jason Pontious, and Nihar Shah. A robust multi-constellation
time-differenced carrier phase solution. GPS Solutions, 2021. Unpublished.
[68] C. Zucca and P. Tavella. The clock model and its relationship with the allan and related
variances. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 52
61
(2):289–296, February 2005. doi: 10.1109/tuffc.2005.1406554. URL https://doi.org/
10.1109/tuffc.2005.1406554.