# Flight-Test Evaluation of Kinematic Precise Point Positioning of Small UAVs

Authors: Jason N. Gross; Ryan M. Watson; Stéphane D’Urso; Yu Gu
Journal: International Journal of Aerospace Engineering
 (2016)

---

## Abstract

An experimental analysis of Global Positioning System (GPS) flight data collected onboard a Small Unmanned Aerial Vehicle (SUAV) is conducted in order to demonstrate that postprocessed kinematic Precise Point Positioning (PPP) solutions with precisions approximately 6 cm 3D Residual Sum of Squares (RSOS) can be obtained on SUAVs that have short duration flights with limited observational periods (i.e., only ~≤5 minutes of data). This is a significant result for the UAV flight testing community because an important and relevant benefit of the PPP technique over traditional Differential GPS (DGPS) techniques, such as Real-Time Kinematic (RTK), is that there is no requirement for maintaining a short baseline separation to a differential GNSS reference station. Because SUAVs are an attractive platform for applications such as aerial surveying, precision agriculture, and remote sensing, this paper offers an experimental evaluation of kinematic PPP estimation strategies using SUAV platform data. In particular, an analysis is presented in which the position solutions that are obtained from postprocessing recorded UAV flight data with various PPP software and strategies are compared to solutions that were obtained using traditional double-differenced ambiguity fixed carrier-phase Differential GPS (CP-DGPS). This offers valuable insight to assist designers of SUAV navigation systems whose applications require precise positioning.


---

## Body

## 1. Introduction

The Precise Point Positioning (PPP) technique was introduced in the late nineties [1, 2] and uses state-space GNSS satellite orbit and clock bias solutions with significantly greater accuracy than their broadcast ephemeris counterparts in order to enable the user-segment to obtain accurate positioning with undifferenced data. The use of undifferenced data means that no GNSS reference station is required to form differential data combinations. To date, PPP technology has matured to the extent that there are now multiple real-time orbit and clock solution products offered by government organizations [3], commercial entities [4, 5], and the International GNSS Service [6]. The PPP technique has revolutionized many geophysical research applications that involve static reference stations and Earth orbiting spacecraft; however, it has not been heavily exploited for applications involving low-cost small UAVs (SUAVs).Many authors have conducted studies to compare the solution accuracy of the PPP technique with double-difference CP-DGPS solutions or other ground “truths.” For example, Colombo et al. show that once a PPP filter has converged, it agrees with double-differenced GPS to within 10 cm [7]. Likewise, for a kinematic vehicular application, Honda et al. demonstrate a few decimeter-level performance with respect to CP-DGPS [8]. In another study, Zhang and Forsberg consider the use of PPP to support missions that require accuracy over very long-ranges (i.e., on the order of many hundreds of kilometers), thereby making double-differences to an individual reference station impractical [9]. In their assessment, Zhang and Forsberg use comparisons of airborne laser altimetry and satellite altimetry products to assess height solution accuracy from PPP and conclude that PPP can produce accuracy at the decimeter-level. In 2009, Bisnath and Gao [10] offered insight on the state-of-the-art of PPP and, in their assessment, demonstrate decimeter-level kinematic PPP of a static reference station after an initial convergence period. Bisnath and Gao conclude that more algorithm development and additional observables are needed to reduce PPP’s convergence period before it can be considered as an RTK alternative. Many have shown that PPP can achieve accuracy levels consistent with CP-DGPS but point out its slow convergence.More recent studies have also considered the impact of multiconstellation GNSS or other navigation systems, such as Inertial Navigation Systems (INS), to yield better accuracy with the PPP approach. For example, Cai et al. evaluated multiconstellation GNSS using the GPS, GLONASS, Galileo, and Beidou and showed improvements in kinematic PPP position solution accuracy and convergence [11]. Likewise, Yigit et al. demonstrated better positioning performance with multiconstellation GNSS, especially whenever there is a relatively short observation duration [12]. In addition, Zhang and Gao include INS within a PPP filter for a kinematic application and show faster solution convergence and accuracy comparable to conventional RTK/INS solutions [13]. More recently, Gross et al. [14] reaffirmed better overall solution accuracy as well as solution robustness to rapid changes in the tropospheric delay induced by abrupt changes in aircraft altitude by comparing both a fused PPP/INS and a PPP-only solution to a postprocessed reference solution that was the NASA Jet Propulsion Laboratory’s submission to the US National Geodetic Survey’s Kinematic GPS Processing Challenge [15]. Finally, Watson et al. [16] provided a comprehensive evaluation of the benefit of incorporating INS under various common PPP error sources by using a Monte Carlo simulation environment and showed that incorporating INS becomes more important depending on the quality of the troposphere model, multipath environment, and quality of PPP orbit and clock products.Despite the extensive literature that offers performance comparisons of different GNSS processing techniques, there is a lack of published studies that demonstrate the obtainable PPP positioning accuracy of SUAV flights that have very limited flight durations. Furthermore, while several past studies point to the slow solution convergence as a pitfall of using PPP technique, the need for real-time solutions is often irrelevant for many SUAV scientific applications (e.g., remote sensing, aerial mapping) that can simply wait for postprocessed solutions. However, because some SUAVs have very limited flight durations, which are on the order of fifteen minutes, or often less, some uncertainty remains as to whether PPP’s slow solution convergence will impact the accuracy of the postprocessed (i.e., filtered and backwards smoothed) short duration SUAV solutions.As such, to fill this knowledge gap, the contribution of this paper is to offer an experimental analysis of PPP techniques when compared to CP-DGPS with data collected onboard SUAV that has very short duration flights. This is particularly relevant to the field robotics community as SUAVs are being more regularly used for more ambitious scientific applications that have stringent requirements on platform positioning. For example, recent experimental evaluations have demonstrated the use of Light Detection and Ranging (LIDAR) for fine-scale mapping with UAVs [17–19]. Accurate georegistration requires platform positioning accuracy in the cm-scale, as these LIDAR systems have cm-scale or better resolution. Additionally, researchers at the University of Kansas have recently instrumented SUAV for radar sounding of remotely located Antarctic ice sheets [20]. While the current Kansas ice sounding UAV is operating at radar frequencies with large wavelengths, it is not a stretch to envision a similar SUAV radar system operating higher frequencies (e.g., microwave) and thus needing cm-scale positioning precision, such as a miniaturized version of NASA’s L-Band UAVSAR instrument [21]. For these SUAV applications, which are typically cited to offer the benefit of being rapidly deployable and useful for remote regions, eliminating the need for GNSS differential reference station and thus reducing overall cost and complexity of the navigation system by leveraging PPP will further open the use of SUAVs for hosting scientific payloads that have stringent positioning requirements.The rest of this paper is organized as follows. Section2 reviews the GNSS observational models for both PPP, which is being evaluated, and double-differenced ambiguity fixed CP-DGPS, which is being used to estimate reference position solutions. In addition, the details of the software implementation employed in this study are discussed. Section 3 discusses the experimental SUAV, GNSS equipment, and flight-test environment used for this study. The results of an experimental comparison study are then detailed in Section 4, and the study’s findings are summarized in Section 5.


## 2. Methodology

### 2.1. Observational Models

We start by considering the generic observation models for GNSS pseudorange and carrier-phase measurements as shown in (1) and (2), which are found in many textbooks [22, 23]. Consider(1)
ρ
i
j
=
R
i
j
+
c
δ
t
i
+
T
i
j
+
I
i
j
+
c
δ
t
j
+
ϵ
ρ
j
,
(2)
Φ
i
j
=
R
i
j
+
c
δ
t
i
+
T
i
j
-
I
i
j
+
c
δ
t
j
+
λ
L
k
N
L
k
j
+
ϵ
Φ
j
.In (1) and (2), the subscript i is used to denote the user and superscript j denotes the GNSS satellite index. Common to both pseudorange and carrier-phase are several error sources, where t
i is the unknown GNSS receiver clock bias in units of seconds, t
j is the GNSS transmitter’s clock bias in units of seconds, c is the speed of light in a vacuum in units of meters per second, T
i
j is the GNSS signal delay due to the refraction in the neutral atmosphere in units of meters, and I
i
j is the phase-advance/pseudorange delay caused by signal refraction through traversing the Earth’s ionosphere (note the sign change between (1) and (2)) in units of meters. The remaining unmodeled error sources are denoted by ϵ and are in units of meters. Within the carrier-phase model (2), there is also unknown integer phase tracking ambiguity that is denoted by N
i
j and is taken from units of carrier-phase cycles to meters through a multiplication with the GNSS carrier wavelength λ
L
k, where k
=
1,2 and L
1
=
19.0 cm and L
2
≈
24.4 cm for the case of GPS. Furthermore, in (1) and (2), the geometric range between the user’s receiver antenna phase center and the satellite’s transmitter antenna phase center is denoted as R
i
j and given by(3)
R
i
j
=
x
j
-
x
i
2
+
y
j
-
y
i
2
+
z
j
-
z
i
2
,where both the user’s position and satellite’s position are modeled in the same Cartesian reference frame, typically either an Earth-Centered-Earth-Fixed frame or an Earth-Centered-Inertial frame. For the PPP technique, accurate solutions to each GNSS transmitter’s clock bias, t
j, and orbital location (
x
j
,
y
j
,
z
j
) are determined using a global network of tracking reference stations and are usedin lieu of the GNSS broadcast ephemeris. However, typically, the satellite locations are provided with respect to the center of gravity of each GNSS spacecraft. Therefore, to be consistent with (3), the attitude of each satellite must be modeled such that the lever arm offset between each satellite’s antenna phase center and its center of gravity is properly considered. To do this, GNSS satellite type dependent known lever arm offsets and attitude models [24, 25] must be included in the processing strategy.In this study, we assume access to dual-frequency GNSS data. As such, we use the dispersive nature of the ionosphere’s delay and form a linear combination of the dual-frequency signals in order to cancel the effect of the ionosphere signal refraction to the first order. This linear combination is denoted as ionospheric-free (IF) combination [22] and given by (4) and (5), for pseudorange and carrier-phase observations, respectively:(4)
ρ
I
F
j
=
ρ
L
1
j
f
1
2
f
1
2
-
f
2
2
-
ρ
L
2
j
f
2
2
f
1
2
-
f
2
2
,
(5)
ϕ
I
F
j
=
ϕ
L
1
j
f
1
2
f
1
2
-
f
2
2
-
ϕ
L
2
j
f
2
2
f
1
2
-
f
2
2
.Note that the coefficients in (4) and (5) (i.e., 2.546 and −1.546 for GPS) sum to 1.0, and thus the modeled common mode error sources, including clock biases and the troposphere delay, remain unchanged in the IF observations. However, when using the IF combination, the unmodeled random errors denoted by ϵ become amplified. For GNSS, the unmodeled error sources consist of thermal noise and receiver front end noise and multipath. For pseudorange observables ϵ
ρ is on the order of 1 meter and for carrier-phase ϵ
ϕ is on the order of 1 millimeter. As such, when using the IF combination, the measurement noise is ~3 meters for pseudorange and 3 millimeters for carrier-phase (i.e., 3
≈
2.54
6
2
+
1.54
6
2).In (1) and (2), the tropospheric delay T
i
j is typically modeled by scaling the zenith direction delay using an elevation angle dependent mapping function to reduce the number of model parameters. Furthermore, this is composed of both wet and dry components of which the dry component is ~90% of the total delay and can be well approximated with a model. As such, the zenith dry delay is typically modeled and the zenith wet delay is estimated as an unknown parameter such that the delay is modeled as shown in(6)
T
i
j
=
T
z
D
r
y
M
D
r
y
e
l
i
j
+
T
z
W
e
t
M
W
e
t
e
l
i
j
,where e
l
i
j is the elevation angle from the nominal user location to satellite j. Within this study, different mapping functions, M
D
r
y
,
M
W
e
t, are utilized depending on the solution strategy and software used as will be detailed in the next section.The remaining modeled error sources are the carrier-phase ambiguities. Traditionally, with kinematic PPP techniques, despite the fact that these are known to be an integer number, these are modeled as floating point parameters in the PPP filter/smoother. However, PPP ambiguity resolution techniques have been developed [26, 27] and this remains an active research area.An important unmodeled error source within theϵ terms is present in (1) and (2) which are the errors induced by multipath reflections.  These errors are attributed to reflections of the GNSS signals being present in the received signal in addition to the line of sight signal. Unfortunately, multipath leads to errors that are non-Gaussian and non-White in nature, such that Kalman or least-squares type estimators are not best suited to mitigate them. Instead, a total least-squares processing technique could be employed, such as the one mentioned in Juang [28]. However, the multipath induced errors for carrier-phase are multiple orders of magnitude smaller than those for pseudorange. Further, for many airborne applications, multipath induced errors, due to signal reflections from man-made objects, are less of a potential problem. In this study, to mitigate multipath induced errors, our PPP estimators will heavily rely on carrier-phase data relative to pseudorange data.


### 2.2. Kinematic PPP with JPL’s GIPSY-OASIS

The first PPP approach considered in this study uses Caltech JPL’s GNSS-Inferred Positioning System and Orbit Analysis Simulation Software (GIPSY-OASIS) 6.2 [29]. GIPSY has been the primary geodetic and positioning software for NASA’s TOPEX/Poseidon [30] and JASON [31] and GRACE [32] low Earth orbiting spacecraft and is operationally used to generate JPL’s precise GPS orbits and clock products for the IGS [33]. GIPSY is licensed for free by Caltech to institutions for use in academic research.When using GIPSY for kinematic PPP, our strategy in this study is to iteratively process the position solution while varying GIPSY configuration parameters in order to converge it to an optimal solution that is free of data outliers. A block diagram of the processing strategy is shown in Figure1, which requires defining some GIPSY terms.Figure 1
Kinematic PPP strategy with GIPSY-OASIS. A wrapper software is used to interactively process the position solution.(i) GNSS Data to Positioning (GD2P). It is GIPSY’s main user interface script for PPP.(ii) Pseudorange Data to Positioning (PR2P). It is GIPSY’s script for pseudorange-only point positioning.(iii) Time Dependent Parameter (TDP). It is GIPSY’s output format for positioning solutions and other solved for parameters (e.g., clock biases, troposphere, and phase biases).(iv) QM File. It is GIPSY’s native binary GNSS measurement format.As indicated in Figure1, for the first iteration, a position solution is estimated using only pseudorange measurement with GIPSY’s PR2P. In addition, the data is translated from the Receiver Independent Exchange Format (RINEX) to a GIPSY binary QM file. During this process, a GNSS data editor is used to flag carrier-phase breaks and remove gross data outliers [34]. For data editing, the GIPSY defaults were used with the exception of the editor requiring a minimum data arc length for a given transmitter. This is due to the short overall observation window of the flights.For the remaining iterations, a subset of GIPSY processing options are varied while accepting the previous position solution (TDP file) as thea priori position solution. Within Figure 1, the configuration options that are varied for each run are as follows.(i) Data Weights. They are measurement noise ratio between pseudorange and carrier-phase measurements. This starts as one-to-one and varies to one-hundred-to-one, where carrier-phase data are modeled as 100 times more precise than pseudorange data.(ii) Postfit Residual Window. Within each GIPSY processing run, multiple passes of a Kalman filter and smoother are conducted. Between each pass, postfit data residuals are evaluated and data are marked as outlier based on defined thresholds and excluded from the next pass. At each pass, the residuals of all data, inlier and outlier, are evaluated and either added back in or excluded from the run. This process is repeated until all data meet the postfit window criteria or a maximum number of iterations are exceeded.(iii) Stochastic Models. The position and wet troposphere delay estimates can be modeled using either white noise about the nominal solution or random walk process noise. Additionally, thea priori
σ magnitude and rate of process noise updates can be set. In particular, the wet zenith delay was estimated as a random walk process, given that the SUAV evaluated is not flying over great distances. Position is initially estimated using random walk, and after a few iterations it is estimated using white noise about ana priori nominal solution.(iv) Minimum Slip. After each filter iteration, jumps in the postfit phase residuals are used to identify the possibility of a carrier-phase break that was missed by the data editor. New breaks are flagged for the next iteration.In addition to the configuration parameters listed above that are varied for each iteration, several other GIPSY options were selected and held fixed in this study. In particular, the troposphere mapping function used is the Niell mapping function [35] and the nominal troposphere delays were set using the atmospheric model populated with a nominal height for scaling and pressure and temperature from the Global Pressure and Temperature model [36]. For the several remaining available GD2P options (e.g., elevation cut-off, tide models), the defaults provided by JPL were used.The specific strategy used for this study with respect to Figure1 is listed in Table 1, for the 10 iterations conducted.Table 1
GIPSY wrapper options of kinematic PPP of SUAV data. Using the process shown in Figure1.
Iteration #
1
2
3
4
5
6
7
8
9
10
Data weight ratio phase/range
1
1
10
100
100
100
100
100
100
100
Range residual window (m)
1.0
e
3
1.0
e
2
10
5.0
5.0
2.5
2.5
2.5
2.5
2.5
Phase residual window (m)
1.0
e
3
1.0
e
2
1.0
1.0
0.5
0.25
0.25
0.25
2.5
e
-
2
2.5
e
-
2
Per-epoch  white noise (m)
1
e
3
1
e
3
1
e
2
1
e
2
1
e
2
1
e
2
50
50
5
5
Postbreak min slip (m)
1
e
12
1
e
12
1
e
12
1
e
12
1
e
12
1.0
1.0
1.0
1.0
1.0While the outlier deletion windows are selectedad hoc based upon experience, in order to assess the effectiveness of postfit data residual analysis and outlier deletion, the GIPSY log files report an approximate χ
2 statistic of the residuals for each filter/smoother and data editing iteration. The value reported is only approximate, because it is reported as if the entire data set were processed as least-squares batch as opposed to a sequential filter. More formally, the Residual Standard Error (RSE) of the estimator for each outlier editing iteration is reported as follows:(7)
σ
R
S
E
=
R
S
O
S
N
-
D
,where R
S
O
S is the Residual Sum of Squares normalized to the data weight of each type of observation, N is the total number of GNSS observations, and D is the total number of parameters estimated. In addition, while it may seem problematic to delete data based upon user selected thresholds, in practice only gross outliers are removed from the filter run. This is because the residual analysis is done in a sequential narrowing window manner, such that such that, first, only gross outliers are deleted and then the residuals are reevaluated, with the possibility of adding data that was flagged as outlier on a previous run as being marked as inlier again (e.g., to handle the case in which one large outlier pollutes all of the residuals for a given epoch). Figure 2 shows an example for the RSE as well as the total percentage of data deleted during a GIPSY run, where the RSE is shown to converge and less than 1.5% of data is shown to have been deleted.Figure 2
Example of GIPSY’s iterative outlier deletion strategy.


### 2.3. Kinematic PPP with RTKLIB

The second software package used for PPP in this study is the open-source RTKLIB [37, 38]. When using RTKLIB, several processing options are available for the user. The following list describes the adopted processing strategy of this study.(i) Kalman Filter Set-Up. The Kalman estimator was configured for both forward filtering and then backwards smoothing.(ii) Elevation Cut-Off. A 5° elevation angle cut-off was used.(iii) Troposphere Model. The troposphere was modeled using the Saastamoinen model and a residual wet zenith delay was estimated.(iv) Ionosphere. The ionosphere free linear combination was used.(v) Misc. Solid Earth tides and ocean loading tides were modeled. Phase windup was considered.(vi) Process Noise. The SUAV dynamics were modeled using a double integrator with process noise driving the acceleration states of 10 m/s2.These options were selected using RTKLIB’s graphical user interface.


### 2.4. Carrier-Phase Double-Differenced Processing for Reference Solutions

To provide reference position solutions for the kinematic PPP solutions, a carrier-phase double-differenced integer ambiguity fixed processing strategy was used. In this case, the GNSS error sources present in (2), in common with a GNSS base station, which is located at a well-known location within a short baseline separation, are canceled via data differencing [22]. In this scenario, the relative position of the UAV with respect to the base station, r
U
,
B, is estimated, and the known location of the base station is added to recover an absolute position estimate of the UAV. For this technique, first, two carrier-phase measurements from the same satellite, j, are differenced between the two user receivers, denoted here as U, for UAV, and B, for base station, to form single-differenced carrier-phase measurements:(8)
Δ
ϕ
U
,
B
j
=
λ
-
1
r
U
,
B
j
+
c
λ
δ
t
i
U
,
B
+
N
U
,
B
j
+
ϵ
ϕ
,
U
,
B
j
,where the satellite clock bias errors, the troposphere delay, and ionosphere delay are eliminated, and the ambiguity N
U
,
B
j remains an integer. Next, to further eliminate the error attributed to the combined users’ receiver clock biases, δ
t
U
,
B, two of the single-differenced carrier-phase measurements from satellites j and k are subtracted to form double-differenced carrier-phase measurements, ∇
Δ
ϕ
U
,
B
j
,
k. The double-differencing process requires the selection of a reference satellite which is typically selected to be a high-elevation satellite (indicated herein with index k) and subtracting its single-differenced measurement from all other available single-differenced measurements. With double-differenced measurements, the only errors that remain are the multipath errors, which are small for carrier-phase observables and airborne applications, and the ambiguity, which is known to be an integer number of wavelengths of the carrier frequency.To resolve the integer ambiguity, a Kalman filter first estimates the carrier-phase ambiguities as floating point parameters, and, then, a separate technique is used to to take advantage of the fact that these states are an integer number of wavelengths. In particular, thede facto technique that is used in this study is known as the Least-squares AMBiguity Decorrelation and Adjustment (LAMBDA) method [39, 40]. The objective of the LAMBDA method is to find an Integer Least Squared Solution (ILS) with respect to the estimated float ambiguities, N
^
ϕ, and a corresponding variance-covariance estimate of the phase ambiguities, P
N
^
ϕ
,
N
^
ϕ [39]. This is done by finding an orthogonal transformation that preservers integer values and decorrelates P
N
^
ϕ
,
N
^
ϕ to a diagonal covariance matrix, such that simple rounding can be employed to estimate the integer states. Following the rounding process, the transformation is reversed to get back into the state-space domain. The cost function that the LAMBDA method optimizes is given by [41](9)
F
N
ϕ
=
N
^
ϕ
-
N
ϕ
T
P
N
^
ϕ
,
N
^
ϕ
-
1
N
^
ϕ
-
N
ϕ
≤
χ
2
,
N
ϕ
∈
Z
n
,where the integer grid search space is defined by χ
2. Once the integer fixed biases have been determined, the relative navigation states are adjusted by assuming that the integer-fixing process is deterministic, such that the nonambiguity states are corrected as(10)
x
p
o
s
i
t
i
o
n
f
i
x
=
x
p
o
s
i
t
i
o
n
f
l
o
a
t
-
P
p
o
s
i
t
i
o
n
,
N
^
ϕ
P
N
^
ϕ
,
N
^
ϕ
N
^
ϕ
-
N
ϕ
,where P refers to the variance-covariance matrix for the floating point ambiguities that were estimated by the Kalman filter, of which particular sections are identified by the subscripts position, which refers to the position states, and N
^
ϕ, which refers to phase ambiguity states.It is important to point out that integer bias fixing will not always be a success in the presence of errors. As such, there must be a validation process. The specific acceptance-test employed in the RTKLIB implementation used in this study was the ratio-test [42]. The ratio-test evaluates how close the float ambiguity estimates are to the best integer ambiguity estimates when compared to the next best integer ambiguity candidate. The best candidate, N
ϕ
1
s
t, and the second best candidate, N
ϕ
2
n
d, are defined as those that minimize the cost function F
(
N
ϕ
) of (9):(11)
Accept
N
ϕ
1
s
t
iff
F
N
ϕ
1
s
t
F
N
ϕ
2
n
d
≤
1
C
,where C is the critical value, which can be derived on the fly to allow a fixed failure rate or set to a constant [42]. For this study, C was set to 3 and held constant. Three is often used in practice, though this is only empirically justified [42]. Throughout this study, only the epochs that were successfully fixed and passed the ratio-test criterion, which ranged from 65 to 85% of all epochs of each flight, were used as the reference solution for conducting the error analysis presented below.


## 2.1. Observational Models

We start by considering the generic observation models for GNSS pseudorange and carrier-phase measurements as shown in (1) and (2), which are found in many textbooks [22, 23]. Consider(1)
ρ
i
j
=
R
i
j
+
c
δ
t
i
+
T
i
j
+
I
i
j
+
c
δ
t
j
+
ϵ
ρ
j
,
(2)
Φ
i
j
=
R
i
j
+
c
δ
t
i
+
T
i
j
-
I
i
j
+
c
δ
t
j
+
λ
L
k
N
L
k
j
+
ϵ
Φ
j
.In (1) and (2), the subscript i is used to denote the user and superscript j denotes the GNSS satellite index. Common to both pseudorange and carrier-phase are several error sources, where t
i is the unknown GNSS receiver clock bias in units of seconds, t
j is the GNSS transmitter’s clock bias in units of seconds, c is the speed of light in a vacuum in units of meters per second, T
i
j is the GNSS signal delay due to the refraction in the neutral atmosphere in units of meters, and I
i
j is the phase-advance/pseudorange delay caused by signal refraction through traversing the Earth’s ionosphere (note the sign change between (1) and (2)) in units of meters. The remaining unmodeled error sources are denoted by ϵ and are in units of meters. Within the carrier-phase model (2), there is also unknown integer phase tracking ambiguity that is denoted by N
i
j and is taken from units of carrier-phase cycles to meters through a multiplication with the GNSS carrier wavelength λ
L
k, where k
=
1,2 and L
1
=
19.0 cm and L
2
≈
24.4 cm for the case of GPS. Furthermore, in (1) and (2), the geometric range between the user’s receiver antenna phase center and the satellite’s transmitter antenna phase center is denoted as R
i
j and given by(3)
R
i
j
=
x
j
-
x
i
2
+
y
j
-
y
i
2
+
z
j
-
z
i
2
,where both the user’s position and satellite’s position are modeled in the same Cartesian reference frame, typically either an Earth-Centered-Earth-Fixed frame or an Earth-Centered-Inertial frame. For the PPP technique, accurate solutions to each GNSS transmitter’s clock bias, t
j, and orbital location (
x
j
,
y
j
,
z
j
) are determined using a global network of tracking reference stations and are usedin lieu of the GNSS broadcast ephemeris. However, typically, the satellite locations are provided with respect to the center of gravity of each GNSS spacecraft. Therefore, to be consistent with (3), the attitude of each satellite must be modeled such that the lever arm offset between each satellite’s antenna phase center and its center of gravity is properly considered. To do this, GNSS satellite type dependent known lever arm offsets and attitude models [24, 25] must be included in the processing strategy.In this study, we assume access to dual-frequency GNSS data. As such, we use the dispersive nature of the ionosphere’s delay and form a linear combination of the dual-frequency signals in order to cancel the effect of the ionosphere signal refraction to the first order. This linear combination is denoted as ionospheric-free (IF) combination [22] and given by (4) and (5), for pseudorange and carrier-phase observations, respectively:(4)
ρ
I
F
j
=
ρ
L
1
j
f
1
2
f
1
2
-
f
2
2
-
ρ
L
2
j
f
2
2
f
1
2
-
f
2
2
,
(5)
ϕ
I
F
j
=
ϕ
L
1
j
f
1
2
f
1
2
-
f
2
2
-
ϕ
L
2
j
f
2
2
f
1
2
-
f
2
2
.Note that the coefficients in (4) and (5) (i.e., 2.546 and −1.546 for GPS) sum to 1.0, and thus the modeled common mode error sources, including clock biases and the troposphere delay, remain unchanged in the IF observations. However, when using the IF combination, the unmodeled random errors denoted by ϵ become amplified. For GNSS, the unmodeled error sources consist of thermal noise and receiver front end noise and multipath. For pseudorange observables ϵ
ρ is on the order of 1 meter and for carrier-phase ϵ
ϕ is on the order of 1 millimeter. As such, when using the IF combination, the measurement noise is ~3 meters for pseudorange and 3 millimeters for carrier-phase (i.e., 3
≈
2.54
6
2
+
1.54
6
2).In (1) and (2), the tropospheric delay T
i
j is typically modeled by scaling the zenith direction delay using an elevation angle dependent mapping function to reduce the number of model parameters. Furthermore, this is composed of both wet and dry components of which the dry component is ~90% of the total delay and can be well approximated with a model. As such, the zenith dry delay is typically modeled and the zenith wet delay is estimated as an unknown parameter such that the delay is modeled as shown in(6)
T
i
j
=
T
z
D
r
y
M
D
r
y
e
l
i
j
+
T
z
W
e
t
M
W
e
t
e
l
i
j
,where e
l
i
j is the elevation angle from the nominal user location to satellite j. Within this study, different mapping functions, M
D
r
y
,
M
W
e
t, are utilized depending on the solution strategy and software used as will be detailed in the next section.The remaining modeled error sources are the carrier-phase ambiguities. Traditionally, with kinematic PPP techniques, despite the fact that these are known to be an integer number, these are modeled as floating point parameters in the PPP filter/smoother. However, PPP ambiguity resolution techniques have been developed [26, 27] and this remains an active research area.An important unmodeled error source within theϵ terms is present in (1) and (2) which are the errors induced by multipath reflections.  These errors are attributed to reflections of the GNSS signals being present in the received signal in addition to the line of sight signal. Unfortunately, multipath leads to errors that are non-Gaussian and non-White in nature, such that Kalman or least-squares type estimators are not best suited to mitigate them. Instead, a total least-squares processing technique could be employed, such as the one mentioned in Juang [28]. However, the multipath induced errors for carrier-phase are multiple orders of magnitude smaller than those for pseudorange. Further, for many airborne applications, multipath induced errors, due to signal reflections from man-made objects, are less of a potential problem. In this study, to mitigate multipath induced errors, our PPP estimators will heavily rely on carrier-phase data relative to pseudorange data.


## 2.2. Kinematic PPP with JPL’s GIPSY-OASIS

The first PPP approach considered in this study uses Caltech JPL’s GNSS-Inferred Positioning System and Orbit Analysis Simulation Software (GIPSY-OASIS) 6.2 [29]. GIPSY has been the primary geodetic and positioning software for NASA’s TOPEX/Poseidon [30] and JASON [31] and GRACE [32] low Earth orbiting spacecraft and is operationally used to generate JPL’s precise GPS orbits and clock products for the IGS [33]. GIPSY is licensed for free by Caltech to institutions for use in academic research.When using GIPSY for kinematic PPP, our strategy in this study is to iteratively process the position solution while varying GIPSY configuration parameters in order to converge it to an optimal solution that is free of data outliers. A block diagram of the processing strategy is shown in Figure1, which requires defining some GIPSY terms.Figure 1
Kinematic PPP strategy with GIPSY-OASIS. A wrapper software is used to interactively process the position solution.(i) GNSS Data to Positioning (GD2P). It is GIPSY’s main user interface script for PPP.(ii) Pseudorange Data to Positioning (PR2P). It is GIPSY’s script for pseudorange-only point positioning.(iii) Time Dependent Parameter (TDP). It is GIPSY’s output format for positioning solutions and other solved for parameters (e.g., clock biases, troposphere, and phase biases).(iv) QM File. It is GIPSY’s native binary GNSS measurement format.As indicated in Figure1, for the first iteration, a position solution is estimated using only pseudorange measurement with GIPSY’s PR2P. In addition, the data is translated from the Receiver Independent Exchange Format (RINEX) to a GIPSY binary QM file. During this process, a GNSS data editor is used to flag carrier-phase breaks and remove gross data outliers [34]. For data editing, the GIPSY defaults were used with the exception of the editor requiring a minimum data arc length for a given transmitter. This is due to the short overall observation window of the flights.For the remaining iterations, a subset of GIPSY processing options are varied while accepting the previous position solution (TDP file) as thea priori position solution. Within Figure 1, the configuration options that are varied for each run are as follows.(i) Data Weights. They are measurement noise ratio between pseudorange and carrier-phase measurements. This starts as one-to-one and varies to one-hundred-to-one, where carrier-phase data are modeled as 100 times more precise than pseudorange data.(ii) Postfit Residual Window. Within each GIPSY processing run, multiple passes of a Kalman filter and smoother are conducted. Between each pass, postfit data residuals are evaluated and data are marked as outlier based on defined thresholds and excluded from the next pass. At each pass, the residuals of all data, inlier and outlier, are evaluated and either added back in or excluded from the run. This process is repeated until all data meet the postfit window criteria or a maximum number of iterations are exceeded.(iii) Stochastic Models. The position and wet troposphere delay estimates can be modeled using either white noise about the nominal solution or random walk process noise. Additionally, thea priori
σ magnitude and rate of process noise updates can be set. In particular, the wet zenith delay was estimated as a random walk process, given that the SUAV evaluated is not flying over great distances. Position is initially estimated using random walk, and after a few iterations it is estimated using white noise about ana priori nominal solution.(iv) Minimum Slip. After each filter iteration, jumps in the postfit phase residuals are used to identify the possibility of a carrier-phase break that was missed by the data editor. New breaks are flagged for the next iteration.In addition to the configuration parameters listed above that are varied for each iteration, several other GIPSY options were selected and held fixed in this study. In particular, the troposphere mapping function used is the Niell mapping function [35] and the nominal troposphere delays were set using the atmospheric model populated with a nominal height for scaling and pressure and temperature from the Global Pressure and Temperature model [36]. For the several remaining available GD2P options (e.g., elevation cut-off, tide models), the defaults provided by JPL were used.The specific strategy used for this study with respect to Figure1 is listed in Table 1, for the 10 iterations conducted.Table 1
GIPSY wrapper options of kinematic PPP of SUAV data. Using the process shown in Figure1.
Iteration #
1
2
3
4
5
6
7
8
9
10
Data weight ratio phase/range
1
1
10
100
100
100
100
100
100
100
Range residual window (m)
1.0
e
3
1.0
e
2
10
5.0
5.0
2.5
2.5
2.5
2.5
2.5
Phase residual window (m)
1.0
e
3
1.0
e
2
1.0
1.0
0.5
0.25
0.25
0.25
2.5
e
-
2
2.5
e
-
2
Per-epoch  white noise (m)
1
e
3
1
e
3
1
e
2
1
e
2
1
e
2
1
e
2
50
50
5
5
Postbreak min slip (m)
1
e
12
1
e
12
1
e
12
1
e
12
1
e
12
1.0
1.0
1.0
1.0
1.0While the outlier deletion windows are selectedad hoc based upon experience, in order to assess the effectiveness of postfit data residual analysis and outlier deletion, the GIPSY log files report an approximate χ
2 statistic of the residuals for each filter/smoother and data editing iteration. The value reported is only approximate, because it is reported as if the entire data set were processed as least-squares batch as opposed to a sequential filter. More formally, the Residual Standard Error (RSE) of the estimator for each outlier editing iteration is reported as follows:(7)
σ
R
S
E
=
R
S
O
S
N
-
D
,where R
S
O
S is the Residual Sum of Squares normalized to the data weight of each type of observation, N is the total number of GNSS observations, and D is the total number of parameters estimated. In addition, while it may seem problematic to delete data based upon user selected thresholds, in practice only gross outliers are removed from the filter run. This is because the residual analysis is done in a sequential narrowing window manner, such that such that, first, only gross outliers are deleted and then the residuals are reevaluated, with the possibility of adding data that was flagged as outlier on a previous run as being marked as inlier again (e.g., to handle the case in which one large outlier pollutes all of the residuals for a given epoch). Figure 2 shows an example for the RSE as well as the total percentage of data deleted during a GIPSY run, where the RSE is shown to converge and less than 1.5% of data is shown to have been deleted.Figure 2
Example of GIPSY’s iterative outlier deletion strategy.


## 2.3. Kinematic PPP with RTKLIB

The second software package used for PPP in this study is the open-source RTKLIB [37, 38]. When using RTKLIB, several processing options are available for the user. The following list describes the adopted processing strategy of this study.(i) Kalman Filter Set-Up. The Kalman estimator was configured for both forward filtering and then backwards smoothing.(ii) Elevation Cut-Off. A 5° elevation angle cut-off was used.(iii) Troposphere Model. The troposphere was modeled using the Saastamoinen model and a residual wet zenith delay was estimated.(iv) Ionosphere. The ionosphere free linear combination was used.(v) Misc. Solid Earth tides and ocean loading tides were modeled. Phase windup was considered.(vi) Process Noise. The SUAV dynamics were modeled using a double integrator with process noise driving the acceleration states of 10 m/s2.These options were selected using RTKLIB’s graphical user interface.


## 2.4. Carrier-Phase Double-Differenced Processing for Reference Solutions

To provide reference position solutions for the kinematic PPP solutions, a carrier-phase double-differenced integer ambiguity fixed processing strategy was used. In this case, the GNSS error sources present in (2), in common with a GNSS base station, which is located at a well-known location within a short baseline separation, are canceled via data differencing [22]. In this scenario, the relative position of the UAV with respect to the base station, r
U
,
B, is estimated, and the known location of the base station is added to recover an absolute position estimate of the UAV. For this technique, first, two carrier-phase measurements from the same satellite, j, are differenced between the two user receivers, denoted here as U, for UAV, and B, for base station, to form single-differenced carrier-phase measurements:(8)
Δ
ϕ
U
,
B
j
=
λ
-
1
r
U
,
B
j
+
c
λ
δ
t
i
U
,
B
+
N
U
,
B
j
+
ϵ
ϕ
,
U
,
B
j
,where the satellite clock bias errors, the troposphere delay, and ionosphere delay are eliminated, and the ambiguity N
U
,
B
j remains an integer. Next, to further eliminate the error attributed to the combined users’ receiver clock biases, δ
t
U
,
B, two of the single-differenced carrier-phase measurements from satellites j and k are subtracted to form double-differenced carrier-phase measurements, ∇
Δ
ϕ
U
,
B
j
,
k. The double-differencing process requires the selection of a reference satellite which is typically selected to be a high-elevation satellite (indicated herein with index k) and subtracting its single-differenced measurement from all other available single-differenced measurements. With double-differenced measurements, the only errors that remain are the multipath errors, which are small for carrier-phase observables and airborne applications, and the ambiguity, which is known to be an integer number of wavelengths of the carrier frequency.To resolve the integer ambiguity, a Kalman filter first estimates the carrier-phase ambiguities as floating point parameters, and, then, a separate technique is used to to take advantage of the fact that these states are an integer number of wavelengths. In particular, thede facto technique that is used in this study is known as the Least-squares AMBiguity Decorrelation and Adjustment (LAMBDA) method [39, 40]. The objective of the LAMBDA method is to find an Integer Least Squared Solution (ILS) with respect to the estimated float ambiguities, N
^
ϕ, and a corresponding variance-covariance estimate of the phase ambiguities, P
N
^
ϕ
,
N
^
ϕ [39]. This is done by finding an orthogonal transformation that preservers integer values and decorrelates P
N
^
ϕ
,
N
^
ϕ to a diagonal covariance matrix, such that simple rounding can be employed to estimate the integer states. Following the rounding process, the transformation is reversed to get back into the state-space domain. The cost function that the LAMBDA method optimizes is given by [41](9)
F
N
ϕ
=
N
^
ϕ
-
N
ϕ
T
P
N
^
ϕ
,
N
^
ϕ
-
1
N
^
ϕ
-
N
ϕ
≤
χ
2
,
N
ϕ
∈
Z
n
,where the integer grid search space is defined by χ
2. Once the integer fixed biases have been determined, the relative navigation states are adjusted by assuming that the integer-fixing process is deterministic, such that the nonambiguity states are corrected as(10)
x
p
o
s
i
t
i
o
n
f
i
x
=
x
p
o
s
i
t
i
o
n
f
l
o
a
t
-
P
p
o
s
i
t
i
o
n
,
N
^
ϕ
P
N
^
ϕ
,
N
^
ϕ
N
^
ϕ
-
N
ϕ
,where P refers to the variance-covariance matrix for the floating point ambiguities that were estimated by the Kalman filter, of which particular sections are identified by the subscripts position, which refers to the position states, and N
^
ϕ, which refers to phase ambiguity states.It is important to point out that integer bias fixing will not always be a success in the presence of errors. As such, there must be a validation process. The specific acceptance-test employed in the RTKLIB implementation used in this study was the ratio-test [42]. The ratio-test evaluates how close the float ambiguity estimates are to the best integer ambiguity estimates when compared to the next best integer ambiguity candidate. The best candidate, N
ϕ
1
s
t, and the second best candidate, N
ϕ
2
n
d, are defined as those that minimize the cost function F
(
N
ϕ
) of (9):(11)
Accept
N
ϕ
1
s
t
iff
F
N
ϕ
1
s
t
F
N
ϕ
2
n
d
≤
1
C
,where C is the critical value, which can be derived on the fly to allow a fixed failure rate or set to a constant [42]. For this study, C was set to 3 and held constant. Three is often used in practice, though this is only empirically justified [42]. Throughout this study, only the epochs that were successfully fixed and passed the ratio-test criterion, which ranged from 65 to 85% of all epochs of each flight, were used as the reference solution for conducting the error analysis presented below.


## 3. Experimental Set-Up

The Phastball SUAV airframe [43] was developed as a modular research platform and has been used for multiple sensor-fusion studies [44–46]. The Phastball Zero SUAV is shown in Figure 3.Figure 3
Phastball SUAV research platform. The two GNSS antennas are separated by an 85.3 cm baseline.For this study, its payload features two Novatel OEM-615 dual-frequency GNSS receivers with antennas separated 85.3 cm along the airframe’s longitudinal axis. A block diagram of Phastball Zero’s data acquisition structure is shown in Figure4.Figure 4
Phastball SUAV data acquisition set-up. Two GNSS receivers are logged in addition to IMU, mechanical vertical gyroscope, and ground pilot input data.In this set-up, one Novatel Receiver is GPS and GLONASS capable and the other is GPS only. Raw GNSS dual-frequency pseudorange and carrier-phase observables are recorded at a rate of 10 Hz using OpenLog serial microSD data-loggers. The GPS/GLONASS Novatel OEM-615 receiver is wirelessly connected through a 900 MHz modem to a base station that sends up RTK differential correctors. Additionally, a mechanical vertical gyroscope that directly measures the UAV’s pitch and roll is included onboard, as well as the pilot input and four Analog Devices ADIS-16405 Microelectromechanical Systems (MEMS) IMUs, which are interfaced through two Netburner MOD-5213 microcontrollers logging data with OpenLog serial stream data-loggers. For postprocessing, the recorded GNSS data and data streams that are interfaced through the Netburner MOD-5213 microcontroller are synchronized by recording the state of the Pulse-Per-Second (PPS) signal from one of the GNSS receivers.In addition to the Phastball SUAV, the experimental set-up consisted of another Novatel OEM-615 dual-frequency receiver serving as a static reference receiver with the antenna mounted on a tripod. Experimental flight-tests for this study were conducted at WVU’s Jackson’s Mill airfield. A typical flight pattern is shown in Figure5.Figure 5
Birdseye view of typical Phastball SUAV flight at the WVU Jackson’s Mill airfield. Image was generated using Google Earth.In total, six data sets of short duration SUAV flights were collected for this study. These consisted of three flight-tests with two data sets per flight. Table2 lists the total flight duration from take-off to landing of the three flights.Table 2
SUAV flight test durations.
Flight #
Duration (s)
1
156.2
2
342.4
3
264.0


## 4. Results and Discussion

### 4.1. Validation of Reference Position Solutions

For comparing the estimation performances of various PPP and single receiver point positioning algorithms, Kalman filter/smoother CP-DGPS solutions were generated for each flight-test with respect to a local GPS reference station set-up at the airfield. When processing the CP-DGPS data, in order to obtain the absolute coordinates of the base station’s tripod, ~6 hours of static data was processing using GIPSY-OASIS PPP in static mode. The software that was used for generating the CP-DGPS reference solutions was RTKLIB 2.4.2.While it is well known that double-difference integer fixed CP-DGPS solutions can obtain cm-level accuracy and precision, the Phastball flight-test configuration provides the opportunity to validate this level of accuracy and precision for the reference position solutions. In particular, the Phastball SUAV is outfitted with two GPS antennas separated by a known baseline distance of 83.5 cm. Note that, with this set-up, it is possible to further leverage the known baseline constraint between the multiple antennas to improve estimation performance of position and attitude [47], which can further be extended to improvements with an array of multiple antennas (i.e., more than 2) [48, 49]. While the array-aided PPP is an interesting research direction that the authors intend to pursue with this SUAV platform, for the present study, the goal is to assess the accuracy of single antenna PPP for SUAVs. Therefore, in order to validate the reference solutions used for each separate antenna, the two ambiguity fixed CP-DGPS solutions for each receiver on the SUAV were differenced, and the magnitude of the difference was used to evaluate how well the known antenna separation distance is estimated. Figure 5 shows an example of the result of this analysis for our second SUAV flight-test.As shown in Figure6, the separation between the two antennas was estimated with both cm-level accuracy and precision. Table 3 shows the antenna separation estimation error for all three flights, where sub-cm accuracy and cm-level precision were obtained.Table 3
85.3 cm antenna baseline estimation performance using the RTK reference solutions for the three SUAV flight tests considered.
Flight #
μ (cm)
σ (cm)
1
0.41
1.4
2
0.68
1.09
3
0.22
0.45Figure 6
Phastball Zero SUAV GNSS antenna baseline separation as estimated by differencing the two RTK GPS solutions that are being used as reference solutions for PPP comparisons and by measuring with a tape measure.


### 4.2. Performance Comparison Study

In order to provide a context, in addition to the PPP solution methods described in Section2, position solutions obtained by processing the pseudorange data only with the GPS broadcast ephemeris using both the GIPSY software and the RTKLIB software are presented. These solutions represent the expected level of accuracy one would obtain from a GPS receiver position solution without any DGPS or PPP products. Short name identifiers for the processing strategies that were compared are shown in Table 4.Table 4
Short names and definitions for the PPP solution strategies evaluated.
Short name
Solution strategy description
PR2P-PP
Dual-frequency pseudorange measurement-only  point positioning using broadcast ephemeris
RTKLIB-PP
Single receiver point positioning using RTKLIB with broadcast ephemeris
GD2P-PPP
Kinematic PPP using JPL final GPS orbit and clock products and the processing method described in Section2.2
RTKLIB-PPP
Kinematic PPP using International GNSS Service combined final GPS orbit and clock products and method described in Section2.3Table5 shows the mean and standard deviation error statistics for the four processing methods considered over the six data sets, and Table 6 shows the overall average mean error and standard deviation error over the six flights, where all solutions are compared to the ambiguity fixed CP-DGPS reference solutions.Table 5
Mean and standard deviation positioning error statistics by positioning approach for each flight.
Receiver
Flight 1
Flight 2
Flight 3
Front
Back
Front
Back
Front
Back
Algorithm
μ (m)
σ (m)
μ (m)
σ (m)
μ (m)
σ (m)
μ (m)
σ (m)
μ (m)
σ (m)
μ (m)
σ (m)
PR2P-PP
1.98
1.32
3.11
2.80
3.20
3.82
2.57
0.85
1.85
0.84
2.27
0.96
RTKLIB-PP
1.77
1.17
3.68
2.23
3.57
3.75
2.99
1.80
1.69
0.64
2.05
0.83
GD2P-PPP
0.61
0.038
0.91
0.20
0.63
0.03
0.99
0.082
1.391
0.029
1.33
0.026
RTKLIB-PPP
1.99
0.048
2.23
0.18
0.59
0.013
0.65
0.066
0.97
0.018
1.013
0.025Table 6
Mean and standard deviation error overall average of the 6 data sets.
Algorithm
⁢Overall average
μ (m)
σ (m)
PR2P-PP
2.50
1.76
RTKLIB-PP
2.62
1.74
GD2P-PPP
0.98
0.067
RTKLIB-PPP
1.23
0.058As shown in Tables5 and 6, when using pseudorange based point positioning, 2-3 meters-level precision and accuracy can be expected. However, when adopting kinematic PPP, decimeter-level accuracy and centimeter-level precision were obtained, despite the fact that these are flight data sets with only a few minutes of observations. Further, to obtain this level of positioning performance, all the user needs to do is collect the raw GNSS measurements for postprocessing.While the level of precision on short duration SUAV flights is centimeter-level, decimeter-level biases remain in the solution. It is expected that these are likely due to the short duration data set and the known slow convergence of PPP. However, an important insight provided by this study is that when considering short duration postprocessing applications of PPP, the slow convergence yields a constant bias, that is, when considering the error sources that must be modeled or mitigated by PPP in (1) and (2) (e.g., troposphere, orbit, clock, and phase ambiguities). Despite the fact that there is no enough data to converge upon the absolute values of these error sources, each of these error sources remains quite constant over the few minutes’ time-scale. Therefore, these decimeter-level biases can easily be mitigated by starting the SUAV at a known location (e.g., by placing a receiver at a take-off location for static PPP) or simply allowing the SUAV to sit static on the runway for an initialization period.To further substantiate why kinematic PPP with short observational periods is still very precise but biased, an additional five-minute SUAV flight that was simulated using a PPP simulation tool developed for a previous study [16] was analyzed. With this tool, all of the GNSS error sources discussed in Section 2.1 are modeled, but because they are simulated, perfect knowledge about all error sources is available for analysis. As such, the performance with respect to a perfectly known truth was assessed. Figure 7 shows an example kinematic PPP forward filter’s position solution errors and phase bias estimation errors, as well as the backwards Rauch-Tung-Striebel [50] Kalman smoothed position and phase bias errors.Figure 7
Simulated position and phase bias errors of a 5-minute SUAV PPP filter (a) and the smoothed position and phase bias errors (b).
(a)
(b)As shown in Figure7, a similar level of positioning performance as presented in the flight data analysis is shown for the simulated flight, that is, decimeter-level accuracy with centimeter-level precision. Further, the phase biases, which are estimated as random constants, do converge but maintain approximately 10 to 15 cm biases after the 5 minutes. Then, during the backwards smoother, the phase bias estimates are nearly constant, leading to considerably smoother positioning errors, as shown, but with decimeter-level biases remaining present. In particular, the mean errors between the forward filter and smoother are only reduced by 1% from 44 cm to 45 cm; however, the standard deviation of the position errors was reduced by more than 25% with a 3D σ reduction from 14.7 cm to 11 cm.


## 4.1. Validation of Reference Position Solutions

For comparing the estimation performances of various PPP and single receiver point positioning algorithms, Kalman filter/smoother CP-DGPS solutions were generated for each flight-test with respect to a local GPS reference station set-up at the airfield. When processing the CP-DGPS data, in order to obtain the absolute coordinates of the base station’s tripod, ~6 hours of static data was processing using GIPSY-OASIS PPP in static mode. The software that was used for generating the CP-DGPS reference solutions was RTKLIB 2.4.2.While it is well known that double-difference integer fixed CP-DGPS solutions can obtain cm-level accuracy and precision, the Phastball flight-test configuration provides the opportunity to validate this level of accuracy and precision for the reference position solutions. In particular, the Phastball SUAV is outfitted with two GPS antennas separated by a known baseline distance of 83.5 cm. Note that, with this set-up, it is possible to further leverage the known baseline constraint between the multiple antennas to improve estimation performance of position and attitude [47], which can further be extended to improvements with an array of multiple antennas (i.e., more than 2) [48, 49]. While the array-aided PPP is an interesting research direction that the authors intend to pursue with this SUAV platform, for the present study, the goal is to assess the accuracy of single antenna PPP for SUAVs. Therefore, in order to validate the reference solutions used for each separate antenna, the two ambiguity fixed CP-DGPS solutions for each receiver on the SUAV were differenced, and the magnitude of the difference was used to evaluate how well the known antenna separation distance is estimated. Figure 5 shows an example of the result of this analysis for our second SUAV flight-test.As shown in Figure6, the separation between the two antennas was estimated with both cm-level accuracy and precision. Table 3 shows the antenna separation estimation error for all three flights, where sub-cm accuracy and cm-level precision were obtained.Table 3
85.3 cm antenna baseline estimation performance using the RTK reference solutions for the three SUAV flight tests considered.
Flight #
μ (cm)
σ (cm)
1
0.41
1.4
2
0.68
1.09
3
0.22
0.45Figure 6
Phastball Zero SUAV GNSS antenna baseline separation as estimated by differencing the two RTK GPS solutions that are being used as reference solutions for PPP comparisons and by measuring with a tape measure.


## 4.2. Performance Comparison Study

In order to provide a context, in addition to the PPP solution methods described in Section2, position solutions obtained by processing the pseudorange data only with the GPS broadcast ephemeris using both the GIPSY software and the RTKLIB software are presented. These solutions represent the expected level of accuracy one would obtain from a GPS receiver position solution without any DGPS or PPP products. Short name identifiers for the processing strategies that were compared are shown in Table 4.Table 4
Short names and definitions for the PPP solution strategies evaluated.
Short name
Solution strategy description
PR2P-PP
Dual-frequency pseudorange measurement-only  point positioning using broadcast ephemeris
RTKLIB-PP
Single receiver point positioning using RTKLIB with broadcast ephemeris
GD2P-PPP
Kinematic PPP using JPL final GPS orbit and clock products and the processing method described in Section2.2
RTKLIB-PPP
Kinematic PPP using International GNSS Service combined final GPS orbit and clock products and method described in Section2.3Table5 shows the mean and standard deviation error statistics for the four processing methods considered over the six data sets, and Table 6 shows the overall average mean error and standard deviation error over the six flights, where all solutions are compared to the ambiguity fixed CP-DGPS reference solutions.Table 5
Mean and standard deviation positioning error statistics by positioning approach for each flight.
Receiver
Flight 1
Flight 2
Flight 3
Front
Back
Front
Back
Front
Back
Algorithm
μ (m)
σ (m)
μ (m)
σ (m)
μ (m)
σ (m)
μ (m)
σ (m)
μ (m)
σ (m)
μ (m)
σ (m)
PR2P-PP
1.98
1.32
3.11
2.80
3.20
3.82
2.57
0.85
1.85
0.84
2.27
0.96
RTKLIB-PP
1.77
1.17
3.68
2.23
3.57
3.75
2.99
1.80
1.69
0.64
2.05
0.83
GD2P-PPP
0.61
0.038
0.91
0.20
0.63
0.03
0.99
0.082
1.391
0.029
1.33
0.026
RTKLIB-PPP
1.99
0.048
2.23
0.18
0.59
0.013
0.65
0.066
0.97
0.018
1.013
0.025Table 6
Mean and standard deviation error overall average of the 6 data sets.
Algorithm
⁢Overall average
μ (m)
σ (m)
PR2P-PP
2.50
1.76
RTKLIB-PP
2.62
1.74
GD2P-PPP
0.98
0.067
RTKLIB-PPP
1.23
0.058As shown in Tables5 and 6, when using pseudorange based point positioning, 2-3 meters-level precision and accuracy can be expected. However, when adopting kinematic PPP, decimeter-level accuracy and centimeter-level precision were obtained, despite the fact that these are flight data sets with only a few minutes of observations. Further, to obtain this level of positioning performance, all the user needs to do is collect the raw GNSS measurements for postprocessing.While the level of precision on short duration SUAV flights is centimeter-level, decimeter-level biases remain in the solution. It is expected that these are likely due to the short duration data set and the known slow convergence of PPP. However, an important insight provided by this study is that when considering short duration postprocessing applications of PPP, the slow convergence yields a constant bias, that is, when considering the error sources that must be modeled or mitigated by PPP in (1) and (2) (e.g., troposphere, orbit, clock, and phase ambiguities). Despite the fact that there is no enough data to converge upon the absolute values of these error sources, each of these error sources remains quite constant over the few minutes’ time-scale. Therefore, these decimeter-level biases can easily be mitigated by starting the SUAV at a known location (e.g., by placing a receiver at a take-off location for static PPP) or simply allowing the SUAV to sit static on the runway for an initialization period.To further substantiate why kinematic PPP with short observational periods is still very precise but biased, an additional five-minute SUAV flight that was simulated using a PPP simulation tool developed for a previous study [16] was analyzed. With this tool, all of the GNSS error sources discussed in Section 2.1 are modeled, but because they are simulated, perfect knowledge about all error sources is available for analysis. As such, the performance with respect to a perfectly known truth was assessed. Figure 7 shows an example kinematic PPP forward filter’s position solution errors and phase bias estimation errors, as well as the backwards Rauch-Tung-Striebel [50] Kalman smoothed position and phase bias errors.Figure 7
Simulated position and phase bias errors of a 5-minute SUAV PPP filter (a) and the smoothed position and phase bias errors (b).
(a)
(b)As shown in Figure7, a similar level of positioning performance as presented in the flight data analysis is shown for the simulated flight, that is, decimeter-level accuracy with centimeter-level precision. Further, the phase biases, which are estimated as random constants, do converge but maintain approximately 10 to 15 cm biases after the 5 minutes. Then, during the backwards smoother, the phase bias estimates are nearly constant, leading to considerably smoother positioning errors, as shown, but with decimeter-level biases remaining present. In particular, the mean errors between the forward filter and smoother are only reduced by 1% from 44 cm to 45 cm; however, the standard deviation of the position errors was reduced by more than 25% with a 3D σ reduction from 14.7 cm to 11 cm.


## 5. Conclusions

An experimental evaluation of kinematic PPP using short duration SUAV data has been presented, as such a study was missing in the literature and the PPP approach has yet to be heavily leveraged in the field robotics and SUAV communities. Through comparison with CP-DGPS reference solutions, it has been shown that approximately 6 cm 3D positioning precision with decimeter-to-meter-level 3D accuracy is obtainable when using PPP even when the flights are only a few minutes in duration. This is of benefit because PPP does not require the user to have access to a differential reference station. This result has been demonstrated with six SUAV flight data sets and two popular GNSS processing software packages. The results of this study are of benefit to many potential SUAV science applications that require precise positioning.



---

*Source: 1259893-2016-07-13.xml*