## Abstract

Global navigation satellite system (GNSS) precise point positioning (PPP) has potential as an alternative or replacement for real-time kinematic (RTK) processing. In this work, we reached for RTK levels of performance without the need for local information through PPP (i.e., centimeter-level positioning that was reached near-instantaneously). This work makes use of information currently available from processing signals from global positioning system (GPS), Galileo, BeiDou-2/3, and GLONASS by fixing ambiguities for the first three constellations on all available frequencies. This processing was done using a four-frequency, four-constellation uncombined decoupled clock model (DCM) that has been expanded as part of this work. The results were tested on 1448 global datasets and showed that instantaneous convergence on average to 2.5 cm error can be achieved for 81% of the stations. These findings were reinforced by the results of epoch-by-epoch processing, as an average of 80% of all single epochs converged below 2.5 cm error at 1σ, as opposed to less than the 0.5% typically observed for classic PPP.

## 1 INTRODUCTION

Global navigation satellite systems (GNSSs) have been the go-to technology for applications requiring outdoor positioning largely due to the availability of the well-known global positioning system (GPS) and, more recently, other constellations including the Russian GLONASS, European Galileo, and Chinese BeiDou. Although GPS was initially designed to allow multi-meter positioning, specifically, tens of meters once selective availability was enabled, and a few meters after it was disabled in 2000 (GPS.gov: Selective Availability, 2021), its potential to provide precise decimeter or even centimeter-level positions was quickly unveiled via the use of measurement differencing. The first iteration came in the form of differential GPS (DGPS), which relies on the presence of a base station near the rover with known coordinates and can achieve decimeter-level positioning that relies on pseudorange measurements only. The incorporation of carrier-phase measurements into DGPS led to real-time kinematics (RTK), which became capable of centimeter-level positioning based on the possibility of resolving carrier-phase ambiguities to integers (Gao et al., 1997; van Diggelen, 1997; Zhodzishsky et al., 1999). However, the performance of RTK relies on the distance to the base station, as atmospheric delays degrade performance the more distant the rover is from the base. In an attempt to resolve this issue, network RTK (NRTK) was developed, in which a network of base stations is used to send corrections to a rover that remains within its boundaries (Geisler, 2006; Jensen & Cannon, 2000; Townsend et al., 1999). Given the number of corrections that needs to be sent to the rover, it became necessary to transition from observation space representation (OSR) to state space representation (SSR) to limit the bandwidth required to transmit this information (Wübbena et al., 2005; Wübbena & Willgalis, 2001). Instead of transmitting range corrections to the rover for each signal and from each base station, NRTK estimates the different error components as states using the network of reference stations and transmits the states to the rover. The error components include satellite orbits, clocks, biases, and ionospheric and tropospheric delays.

Both RTK and NRTK assume the existence of nearby receivers that can be used as base stations. Unfortunately, these networks are not widely available and are concentrated only in major cities and other densely populated areas. To provide precise positioning for users outside of the NRTK networks, precise point positioning (PPP) was invented (Kouba & Héroux, 2001; Muellerschoen et al., 2001; Zumberge et al., 1997). The principle underlying PPP is similar to NRTK, in the sense that SSR satellite products are computed from a network of stations and are then broadcast to the user. The main differences in the performance of PPP and NRTK lie in presence/absence of precise atmospheric corrections due to the size of the networks, as PPP with regional atmospheric corrections can provide an NRTK-like level of performance (Banville et al., 2014; Teunissen et al., 2010; Wübbena et al., 2005). PPP provides global coverage via the use of global reference stations spread around the world; it does not focus on specific areas. The existence of a worldwide network of stations permits satellite products to be estimated for all GNSS satellites, regardless of their location. This permits users at all locations to access the corrections for the satellites that are visible in their immediate areas. In addition, the reference station network is relatively sparse (as few as tens of stations (Zumberge et al., 1997)) compared to NRTK. This facilitates better decoupling and de-correlation between the various GNSS error components. However, one consequence of these sparse networks is that the atmospheric delays (i.e., those from the troposphere and ionosphere) cannot be estimated precisely. Therefore, the atmospheric delays on the PPP user side need to be estimated without specific a priori information. This significantly increases the number of states that need to be estimated. These factors lead to PPP being heavily dependent on the quality and number of measurements. Indeed, the high number of states, together with the absence of precise atmospheric corrections, leads to the fact that PPP requires minutes to tens of minutes to achieve centimeter-level positioning (Bisnath & Gao, 2009).

In an attempt to reduce the PPP convergence time and achieve NRTK levels of performance, (i.e.,centimeter-level positioning near-instantaneously), research has been performed to augment the global corrections with regional atmospheric corrections by generating precise atmospheric information from a dense network and maintaining the rover as present within the network (Nadarajah et al., 2018; Psychas & Verhagen, 2020; Zhang et al., 2011). However, while this approach allows for RTK-like performance with PPP, it still relies on the use of local reference stations, and thus does not address the need for a global precise positioning technique with instantaneous centimeter-level positioning. Other research has been done to improve PPP through the use of ambiguity resolution (AR), multiple frequencies, and multiple constellations. AR has been identified as crucial to reducing the PPP convergence time to facilitate positioning at a few centimeters within a few minutes (Bertiger et al., 2010; Collins et al., 2008; Ge et al., 2008b; Geng et al., 2012; Laurichesse et al., 2009; Mervart et al., 2008; Naciri & Bisnath, 2020 2021b; Teunissen et al., 2010). However, AR alone is not enough to reach RTK levels as the process requires several minutes to achieve this level of performance. Recent advances in GNSS constellations have led to improvements in the performance of PPP, which is a measurement-dependent technique. Additional GNSS constellations that complement GPS have emerged in recent years; all of these constellations can broadcast signals on as many as four different frequency bands (Naciri et al., 2021). This plethora of constellations and signals means that users have access to many more measurements than were available only a few years ago, thereby allowing for improvements in the performance of PPP (Geng et al., 2020; Katsigianni et al., 2019; Naciri & Bisnath, 2021a; Psychas et al., 2021; Xiao et al., 2019; Xin et al., 2020). For instance, Geng et al. (2020) showed that 6.1 minutes were required to achieve PPP-AR when using triple-frequency data from GPS, Galileo, BeiDou, and QZSS, as opposed to the 9.2 minutes required when using dual-frequency data. Similarly, Naciri & Bisnath (2021a) showed that only two minutes were required on average at 1*σ* to reach and settle below a 10 cm horizontal error when fixing ambiguities on GPS, Galileo, and BeiDou using as many as three frequencies. Hexagon, which is a commercial company, has released its latest correction service that permits centimeter-level positioning in fewer than two minutes at 2*σ* without regional corrections; details on how such performance is achieved have not been provided (RTK From The Sky | Hexagon Autonomy & Positioning, 2021). Psychas et al. (2021) showed the importance of increasing frequencies for accelerating convergence, as well as the critical role of the separation between frequencies in increasing the fixing success rate. The study was based on Galileo E1/E5a/E5b/E5/E6 and GPS L1/L2 signals. Xin et al. (2020) also highlighted the importance of a third frequency choice, as they demonstrated that the E1/E5a/E6 combination of Galileo frequencies led to the best PPP-AR performance in a three-frequency solution. These publications demonstrate the benefits of leveraging additional frequencies and constellations to improve PPP performance.

The goal of this research study is to reach for RTK levels of performance (i.e., within few centimeters and near-instantaneously) using only global corrections without local augmentation. The methodology involves achieving the best possible PPP solution using as many available measurements as possible and comparing the performance to RTK thresholds and requirements. Additionally, it is not only necessary to use as many possible constellations and frequencies, it is also important to fix ambiguities associated with these additional signals because the fixing reduces the PPP convergence time. Given the varying number of frequencies that are broadcast by and within each constellation, it is critical to have a flexible PPP-AR model that is capable of adjusting to the different available signals and fixing ambiguities to align with their correct integer values. This research assesses how close PPP is to becoming a single-receiver technique capable of instantaneously providing centimeter-level positions using only global corrections. The novelty of this research lies in:

The expansion of the classic dual-frequency combined decoupled clock model (DCM) into an uncombined four-frequency model;

The flexible processing of multiple frequency combinations, depending on their availability, even within the same constellation;

The processing of all four GNSS (GPS, GLONASS, Galileo, and BeiDou-2/3), and the capacity to perform ambiguity resolution on GPS, Galileo, and BeiDou-2/3;

The use of PPP performance thresholds that are stricter than those typically used in the literature to facilitate a better comparison with RTK-quality performance. This analysis involves stricter convergence and root mean square (RMS) error thresholds (e.g., convergence to 2.5 cm rather than 10 cm), and epoch-by-epoch processing; and

The analysis of the equivalence in redundancy/degrees of freedom between the DCM and classic PPP, as well as an analysis of the redundancy in the multi-frequency, multi-constellation context.

This manuscript begins with a description of the newly-developed PPP-AR model that facilitates flexible multi-frequency multi-GNSS PPP-AR. The mathematical model section continues with a redundancy analysis in which equivalence between the DCM and the classic model is established and the impact of new signals on the redundancy of the PPP model is evaluated. A description of the data and processing strategies is included, followed by results from the processing using a typical PPP analysis. The penultimate section contains comparisons with RTK-quality performance using the strict thresholds imposed as part of the analysis. This section also contains an analysis of the possibility of robust centimeter-level positioning in an epoch-by-epoch PPP solution. The manuscript ends with conclusions and a discussion of future work.

## 2 MODEL DESCRIPTION AND REDUNDANCY ANALYSIS

Ambiguity resolution plays a critical role in this work, as a flexible model that can process multiple constellations and frequencies is needed. This section describes a mathematical model that was developed for this work by expanding the classic combined dual-frequency DCM model to create an uncombined four-frequency model. The newly-developed model requires the choice of one reference satellite per constellation, as well as the estimation of numerous receiver biases; these are added to the total number of states to be estimated. Despite these additional state terms, the system is equivalent to a classic PPP model; this will be demonstrated in a redundancy/degrees of freedom analysis in which both the equivalence and the impact of the additional constellations and frequencies are demonstrated.

### 2.1 Mathematical Model

Some of the most popular models used for ambiguity resolution-enabled PPP include the fractional cycle bias (FCB) model (Ge et al., 2008a), the integer recovery clock (IRC) model (Laurichesse et al., 2009), and the decoupled clock model (DCM) (Collins et al., 2008). The DCM is used in this research. This model is based on the assumption that carrier-phase and pseudorange measurements are not synchronized to the same level of precision. In order to solve this difference in the respective levels of precision, two sets of clocks - pseudorange clocks, and carrier-phase clocks - are estimated. Hence, the clocks are decoupled. Decoupling the clocks leads to a loss of datum from the carrier-phase measurements, which was previously set by the pseudorange clocks. The phase datum is recovered by selecting one satellite as a reference and setting its ambiguities to arbitrary values against which other satellites’ ambiguities are implicitly differenced.

The DCM was initially introduced by Collins (2008) for GPS processing. It was based on the formation of the ionosphere-free combination of pseudorange and carrier-phase measurements, as well as the Melbourne-Wübbena (MW) combination. These linear combinations allow for the fixing of the narrowlane (NL) and widelane (WL) ambiguities. Given that the model is based on the formation of linear combinations, the ionospheric delay is eliminated through the ionophere-free combination. To allow for the estimation of the state term (for ionospheric-constraining), the model was later extended to an extended DCM (EDCM), in which the MW combination is split into two separate combinations, i.e., the NL pseudorange and WL carrier-phase (Collins et al., 2012).

In the current GNSS context, a variety of signals are broadcast by multiple constellations with each satellite broadcasting on as many as four or five frequencies. This has led to a significant increase in the number of possible linear combinations between measurements if one were to use more than two frequencies. Therefore, uncombined processing, in which raw measurements are processed without forming linear combinations, becomes the preferable approach. A version of the DCM in which dual-frequency uncombined measurements were used instead of combined measurements has been described by Naciri & Bisnath (2021b); Seepersad (2018). As part of this research, the DCM has been expanded to include uncombined triple-frequency processing, with the mathematical derivation detailed by Naciri & Bisnath (2021b) for the Galileo constellation. Additionally, the flexible multi-constellation, dual- and triple-frequency model was introduced by Naciri & Bisnath (2021a), in which dual-frequency and triple-frequency measurements are processed and fixed simultaneously for GPS, Galileo, and BeiDou. In the current manuscript, simultaneous processing of dual-, triple-, and quadruple-frequency measurements is performed via an upgrade of the model introduced by Naciri & Bisnath (2021a ,2021b). To avoid repeating the same mathematical derivations, only the final set of equations is shown, as the upgrade from triple-frequency to quadruple-frequency DCM is similar to the derivation shown in Section 4 of Naciri & Bisnath (2021b). The complete uncombined quadruple-frequency DCM derived as part of this work is summarized in Equation (1):

1

In Equation (1), for satellite *s*, receiver *r*, and frequency *i* ∈ {1,2,3,4}, and are the pseudorange and carrier-phase measurements, respectively, is the geometric range between the receiver and the satellite, *c* is the speed of light, *γ _{i}* is the ratio of frequencies between the first frequency and frequency

*i*, and

*λ*is the wavelength for frequency

_{i}*i*. The satellite corrections from the network side are, in addition to the satellite orbits in , the satellite clocks

*dt*, the satellite code biases , and the satellite phase biases . The estimated atmospheric delays are the zenith wet tropospheric delay

^{s}*T*to which the function is applied to map it to the slant direction, as well as the ionospheric delay on the first frequency , which is mapped to the delay on frequency

_{r}*i*through the factor

*γ*. A more detailed look at the receiver biases and clocks is necessary to understand the newly-developed DCM.

_{i}and are the decoupled pseudorange and carrier-phase receiver clocks, respectively. The estimation of two receiver clocks is due to the decoupling of clocks. As a consequence, the carrier-phase measurements lose their datum. The datum in those measurements is recovered through selecting a reference satellite, not estimating its ambiguities, and setting them to arbitrary integer values. As a consequence, the datum associated with each frequency’s carrier-phase measurements is set by the reference satellite’s ambiguity at that frequency. Therefore, is the integer single-differenced ambiguity for satellite

*s*relative to the reference satellite. The differencing is done implicitly when choosing one satellite as reference, as there is no need to manually difference measurements or ambiguities;*δt*_{12},*δt*_{13}, and*δt*_{14}are referred to as the second, third and fourth frequency’s receiver phase biases. The three state terms are all receiver-related and they contain combinations of receiver biases as well as the reference satellite’s ambiguities. These three state terms are direct consequences of the DCM and the implicit differencing, as they show up when arranging the reference satellite’s measurement equations. Each of the three state terms is specific to the carrier-phase measurements at a specific frequency;*IFB*_{r,3}and*IFB*_{r,4}are inter-frequency pseudorange biases. These biases are common to all multi-frequency models, and they arise from the fact that other state terms cannot absorb all of the code biases at all of the frequencies. The first and second frequency’s receiver code biases are absorbed by the pseudorange receiver clock and the ionospheric delays, but the third and fourth frequency’s receiver code biases have to be estimated as additional state terms, as they would otherwise show up in the post-fit residuals.

Given that the DCM does not make assumptions on the nature of the biases or their behavior, all of the receiver bias state terms are estimated as white noise processes, the same way that clock terms are estimated. In addition to the pseudorange and carrier-phase clocks, these include the inter-frequency code biases *IFB*_{r,3} and *IFB*_{r,4}, and the receiver phase biases *δt*_{12}, *δt*_{13}, and *δt*_{14}.

The model above has been described for a quadruple-frequency satellite. In case of a dual-frequency or triple-frequency satellite, the model remains the same, with the exception of the third- and fourth-frequency and fourth frequency measurements being omitted for dual-frequency, and triple-frequency satellites, respectively. As a consequence, some of the state terms do not need to be estimated when processing fewer frequencies for a specific satellite. For example, *δt*_{14} and *IFB*_{r,4} will not be estimated for a triple-frequency satellite. Another consequence of the capacity for flexible processing of multiple frequencies is that careful consideration must be given to the reference satellite. Given that the datum on each frequency’s phase measurements is set by the reference satellite’s ambiguity on that frequency, it is important to ensure that the reference satellite has a number of frequencies that matches the maximum number of frequencies being processed. Other strategies can be adopted, such as selecting one reference satellite per frequency, rather than having one reference satellite for all frequencies. A more detailed explanation of the reference satellite choice based on frequency is provided by Naciri & Bisnath (2021a).

### 2.2 Redundancy Analysis

In PPP, one typically estimates constellation-dependent receiver-related states; in this case, one set of clocks is estimated per constellation. This is because constellations are typically tracked by different processing units in the receivers, and because of the different time scales between constellations and the presence of inter-system time biases (El-Mowafy et al., 2016; Hakansson et al., 2017). The presence of these constellation-related biases implies that, in addition to receiver clocks, receiver biases are also constellation-dependent. Therefore, in the model described in Equation (1), there is one set of receiver clocks and biases per constellation, consisting of . As described in the previous section, the receiver phase bias state terms are correlated to the reference satellite’s integer ambiguities. Therefore, it is also required to have one reference satellite per constellation, as the receiver biases are constellation-dependent. Furthermore, both BeiDou-2 and BeiDou-3 experience clock and time delay biases between generations; this means that the biases can either be estimated or that the generations need to be treated as two different constellations (Jiao et al., 2020). The second option was the one selected for this work; therefore, one reference satellite was chosen for each constellation for which AR is performed. These include GPS, Galileo, BeiDou-2, and BeiDou-3, with each reference satellite having measurements on the maximum number of frequencies being tracked by all satellites in its constellation.

Fortunately, the presence of numerous reference satellites has no impact on the PPP model and its degrees of freedom. As explained above, the reference satellite’s ambiguities are not estimated but are set to arbitrary integers. Therefore, the additional receiver-related state terms which are specific for the DCM (, *IFB*_{r,3} *IFB*_{r,4}, *δt*_{12}, *δt*_{13}, *δt*_{14}) are compensated by the reference satellite’s ambiguities that do not need to be estimated. In fact, the DCM is equivalent to a classic PPP model in terms of redundancy and degrees of freedom, as the DCM is a reformulation of the classic PPP model in order to access the integer nature of the ambiguities. Next, we examine the aforementioned equivalence. To achieve this goal, let *N ^{c}* be the number of constellations in use,

*N*the total number of satellites regardless of the constellation, and

^{s}*N*the number of frequencies processed. For this demonstration, we consider a subset of satellites that are broadcasting on the same number of frequencies. This demonstration will remain valid and the derivations will remain the same for the more realistic scenarios in which satellites are broadcasting signals on a varying number of frequencies. The number of measurements and of states to be estimated depending on the number of constellations, satellites and frequencies is:

^{f}Receiver coordinates

*X, Y, Z*(**3**), and zenith wet tropospheric delay*T*(_{r}**1**) are estimated regardless of*N*,^{c}*N*, or^{s}*N*.^{f}Receiver code and phase clocks are estimated as long as

*N*≥ 2. One set of states is estimated per constellation that includes^{f}**N**sets.^{c}The third frequency

*IFB*_{r,3}and fourth frequency*IFB*_{r,4}receiver code biases are only estimated if*N*≥ 3 or^{f}*N*= 4, respectively. If these conditions are fulfilled, one state term is estimated per constellation. Therefore, the number of receiver code biases to be estimated is (^{f}**N**– 2)^{f}**N**.^{c}The second frequency

*δt*_{12}, third frequency*δt*_{13}, and fourth frequency*δt*_{14}receiver phase biases are only estimated if*N*≥ 2,^{f}*N*≥ 3, or^{f}*N*= 4, respectively. If these conditions are fulfilled, one state term is estimated per constellation. Therefore, the number of receiver phase biases to be estimated is (^{f}**N**– 1)^{f}**N**.^{c}One slant ionospheric delay , is estimated per satellite meaning

**N**ionospheric state terms. There are also as many ambiguities per satellite as there are frequencies (^{s}*N*ambiguity state terms per satellite). Since ambiguities are not estimated for the^{f}*N*reference satellites, the total number of ambiguity state terms is^{c}**N**(^{f}**N**–^{s}**N**).^{c}The measurements include one code and one phase measurement per frequency per satellite. Therefore, the total number of measurements is 2

**N**.^{s}N^{f}

Therefore, considering that *N ^{meas}* measurements are available and

*N*states are to be estimated, the number of degrees of freedom

^{states}*M*in the uncombined multi-frequency DCM is summarized in Equation (2):

2

In a classic PPP model, many of the state terms listed above are not estimated. The states estimated in such a model include **3** receiver coordinates, **1** wet zenith tropospheric delay, **N ^{c}** receiver clocks, (

**N**– 2)

^{f}**N**receiver code biases,

^{c}**N**ionospheric delays, and

^{s}**N**ambiguities. Therefore, the redundancy in a classic PPP model can be calculated as shown in Equation (3):

^{f}N^{s}3

As shown in Equations (2) and (3), both models have the same degrees of freedom, despite the fact that DCM requires the choice of one reference satellite per constellation. The equivalence is due to the additional state terms in the DCM that appear as a consequence of not estimating the reference satellites’ ambiguities, thereby balancing the redundancy equation.

The degree of freedom equation can be used to understand the effect of the additional frequencies and constellations on PPP models, and subsequently their performance. In light of the current GNSS context in which many constellations and frequencies are available to users, it is important to understand how these additional measurements affect PPP processing. Figure 1 shows the global distribution of the number of visible satellites on average over 24 hours. As shown, an average of nearly 50 satellites can be observed on the Asian continent over the course of one day, with a minimum of over 25 satellites detected in the rest of the world. This high number of visible satellites ensures better satellite geometry for users around the globe, and thus better position dilution of precision (PDOP). This increase in the number of satellites is accompanied by an increase in the number of broadcast frequencies, as shown in Table 1. The table shows that all constellations have satellites that are broadcasting on three different frequencies, with GLONASS mainly broadcasting on two frequencies; GLONASS is currently being upgraded with newer generation multi-frequency satellites. Some constellations (e.g., Galileo and BeiDou-3) broadcast on four frequencies, highlighting the significance of the work presented in this manuscript.

To gain a full understanding of the impact of all the additional satellites, constellations, and frequencies, the degree of freedom function was plotted as a function of the three parameters (number of satellites, constellations and frequencies) as shown in Figure 2.

A clear trend can be noticed with respect to the number of satellites, as redundancy increases as the number of satellites increases, regardless of the number of constellations or frequencies. This trend is expected, as additional satellites introduce additional measurements that help with the estimation of the state terms, some of which are common to all satellites (e.g., receiver coordinates, clocks, and biases). An increase in the number of frequencies has a similar effect on the redundancy, as the latter also increases with the number of frequencies processed. Furthermore, the effect of the number of frequencies on the redundancy is more noticeable with a higher number of satellites. Adding more frequencies increases the slope of the redundancy plots as a function of the number of satellites. Therefore, having a higher number of frequencies associated with a higher number of satellites leads to a greater increase in the redundancy compared to what would be observed with a lower number of satellites. This observation makes sense, as each frequency adds two measurements (pseudorange and carrier-phase) and one additional state to be estimated (ambiguity). Therefore, this increases the redundancy by one, considering that receiver-related biases have already been estimated via the other satellites. As a consequence, adding one frequency when 10 satellites are available leads to an increase in the redundancy of 10. This can be contrasted to the addition of one frequency when 40 satellites are available, which leads to an increase in the redundancy of 40.

The final parameter to be analyzed is the impact of the number of constellations on the redundancy. When evaluating a single number of frequencies, it appears that adding more constellations has an adverse effect on the redundancy, as each additional constellation shifts the redundancy graph toward the lower values. This shift is a consequence of the existence of additional receiver state terms that need to be estimated whenever a new constellation is added, including the receiver clocks and code and phase biases. These additional state terms need to be estimated and cannot be modeled easily, as their behavior is not stable and is receiver-dependent (Wang et al., 2020). However, one must keep in mind that tracking a higher number of satellites requires the use of multiple constellations, as only a limited number of satellites can be tracked with a single constellation due to the relatively limited number of satellites in orbit per constellation. Therefore, an ideal solution would be to process satellites from all constellations on as many frequencies as available. This will permit to achieve the highest redundancy in the PPP model. Including as many measurements in the processing allows for the robustness of the solution and introduces potential improvements in performance, although ever-increasing the redundancy is not a goal in itself due to the law of diminishing returns.

## 3 PROCESSING AND ANALYSIS STRATEGY

The model and algorithms described in the previous section were implemented in the in-house multi-frequency, multi-constellation York-PPP engine at the GNSS lab of York University. To perform the analysis, a global set of stations was selected and observations from these stations were processed with the York-PPP engine. A map of stations used in the processing is shown in Figure 3. Twenty-seven IGS stations were selected by ensuring that each supported three GPS frequencies, four Galileo frequencies, as well as the recent BeiDou-3. Observation data from these stations were processed over one week between day 305 and day 311 of the year 2021. The data were processed at a rate of 30 seconds, and were split into three hour-long datasets that were processed individually in the kinematic mode without assuming that they were static. A total of 1448 three-hour datasets were processed. This number was considered high enough to permit us to reach conclusions on the performance of each processing strategy. The parameters and strategies used to estimate the various states are summarized in Table 2. The frequencies that were selected for processing are as follows:

GPS: L1 (1575.42 MHz), L2 (1227.60 MHz), and L5 (1176.45 MHz);

GLONASS: G1 , and G2 ;

Galileo: E1 (1575.42 MHz), E5b (1207.140 MHz), E5a (1176.45 MHz), and E6 (1278.75 MHz);

BeiDou-2: B1-2 (1561.098 MHz), B2b (1207.140 MHz), and B3 (1268.52 MHz); and

BeiDou-3: B1 (1575.42 MHz), B2a (1176.45 MHz), and B2b (1207.140 MHz).

Satellite products are required to perform PPP-AR on these data. In this work, rapid products were used, with satellite orbits and clocks provided by GFZ (GeoForschungsZentrum) (Deng et al., 2017) and the compatible satellite code and phase biases provided by CNES (Centre Nationale d’Etudes Spatiales) (Laurichesse & Blot, 2016). The products are provided in post-processing RINEX format. Code and phase biases are available in an observable specific bias (OSB) format. This format (as opposed to differential biases) simplifies the use of satellite biases on the user side. Satellite and receiver antenna corrections were applied using the IGS14 ANTEX file (Schmid et al., 2016). For the signals with antenna corrections available in the ANTEX file, the corrections were applied. For satellites with missing antenna corrections, corrections from nearby available signals were used. For example, because GPS L5 antenna corrections are missing from the ANTEX file, GPS L2 antenna corrections were used on the assumption that they are likely to be similar to the L5 corrections. A similar issue exists with respect to BeiDou-3, as ANTEX corrections are not available for its B1 and B2a signals. In this case, B1-2 and B2b corrections were applied, respectively, as the frequencies are relatively close. Furthermore, all pseudorange and carrier-phase measurements from the same satellite were given the same weight regardless of frequency. Tuning the measurement weight based on measurement noise will be addressed in future work. Additionally, BeiDou measurements were assigned half the weight of other constellations. To avoid any unwanted behavior because of the quality of the products or its signals due to BeiDou-3 being a relatively new constellation, the decision to assign less weight to this constellation ensures that it does not have a negative impact the solution while still permitting to benefit from the additional measurements. The choice of a 50% weight was empirical and was based on processing experience. Future work will include additional investigation into the fine-tuning of these weights. In addition to this stochastic modeling, due to the existence of inter-system biases between BeiDou-2 and BeiDou-3, the two generations were treated as two different constellations in the PPP engine. This means that different receiver states and different reference satellites were chosen for each generation. Furthermore, the BeiDou-2 group delay variations were corrected using the frequency- and elevation-dependent corrections provided by Wanninger & Beer (2015). Other corrections, such as the phase wind-up, relativistic effect, and Earth rotation were applied following the IERS conventions (Kouba & Mireault, 1998).

Ambiguity resolution was performed using the modified Least-squares AMBiguity Decorrelation Adjustment (mLAMBDA) method (Chang et al., 2005). The uncombined float ambiguities on each frequency were fed into the mLAMBDA function together with their covariances. Since the estimated ambiguities were already implicitly single-differenced, they do not contain any non-integer biases; thus, additional differencing is not required. The resulting fixed ambiguities from mLAMBDA were used to update the state terms and a standard ratio test was used to validate the fixed ambiguities. Ambiguity resolution was performed on all the ambiguities without selecting a subset to highlight the impact of the additional frequencies on the success rate of the fixing process. Partial ambiguity fixing will be considered as future work.

An epoch-by-epoch solution was also generated as part of this analysis. This solution is similar to the other solutions presented in this manuscript with respect to processing strategy, except that the filter was reset for each epoch. This way, information from previous epochs is not carried forward so that the filter behaves as though it is a single point positioning solution (SPP) with additional measurement processing enhancements.

In terms of ionospheric activity during the period of processing, Figure 4 shows the hourly average Kp and f10.7 indices, which were retrieved from https://omni-web.gsfc.nasa.gov/form/dx1.html. The Kp index characterizes the magnitude of the geomagnetic storms and is used to describe the disturbances in the Earth’s magnetic field. Its value ranges from 0 to 9. The figure shows that, for most days, the index was below 4, indicating quiet to unsettled conditions, except for a one-day period during which the Kp index reached 7, indicating a strong geomagnetic storm. By contrast, the f10.7 index is an indicator of solar activity with values ranging from 50 to greater than 400. Given that the period of interest was during a period of low solar activity, the f10.7 index was relatively low.

## 4 RESULTS

The datasets described in Section 3 were processed using the algorithms described in Section 2 following the strategies discussed in Section 3. The results of the processing are summarized in the subsections to follow. The first subsection focuses on methods used to analyze the overall performance of the proposed algorithms in a typical PPP analysis. This analysis is done by examining the RMS and convergence time metrics and comparing them in different frequency and processing modes, as well as to values reported in the literature. The second subsection focuses on drawing comparisons between the performance of the proposed algorithms and RTK. This comparison is done using stricter comparison thresholds that are equivalent to those that could be used in RTK, as well as by analyzing the performance of an epoch-by-epoch PPP solution with all the measurement enhancements and determining its performance compared to that expected for RTK.

### 4.1 Solution Performance

In this subsection, the results of an analysis of the average performance of the flexible multi-constellation multi-frequency PPP-AR solution using two, two/three, and two/three/four-frequency measurements are presented. Both float and fixed solutions are analyzed. The datasets described in Section 3 were first processed in their entirety without any pre-filtering or removal of outliers. Given that the goal from this work is to assess the average performance from the different processing modes, the statistics in this section were generated by analyzing the average solutions based on all datasets. Because no outliers were removed in the processing procedure, rather than assessing only the averages across all datasets (which would contain outlier results), the average solutions based on three different percentiles (100^{th}, 95^{th} and 67^{th}) were evaluated. For example, the average 100^{th} solution is generated by computing the epoch-by-epoch average of all datasets, while the average 95^{th} solution is generated by computing the epoch-by-epoch average of the 95% lowest position errors. The results from averaging out positioning errors from all the datasets as described in Section 3 are summarized in Figure 5. Given our goal of achieving centimeter-level positioning near-instantaneously, the initial processing period is of primary interest with a focus on sub-decimeter level accuracies. Thus, the figure includes only the first 20 minutes of processing and up to 20 cm of error.

Figure 5 shows the average horizontal and vertical error time series from all 1448 datasets at the 100^{th}, 95^{th} (2*σ*), and 67^{th} (1*σ*) percentiles. Both float and fixed solutions (with and without ambiguity resolution, respectively) are shown for each frequency combination, although only one float solution is shown regardless of the frequency combination. The reason for only showing one float solution is that the additional frequencies have no impact on the float solutions, as previously described by Geng et al. (2020); Naciri & Bisnath (2021a ,2021b); Psychas et al. (2021) who showed that additional frequencies do not lead to improvements in e.g., the satellite geometry. By contrast, the additional frequencies have a noticeable impact on the fixed solutions at all percentiles because adding frequencies leads to correct fixing of ambiguities earlier in the processing. Indeed, adding more frequencies not only reduces convergence time, but it also reduces initial errors. The initial errors of the dual, triple, and quadruple-frequency solutions (referred to as ’all freq.’ in the figure) reach values as low as 12 cm at the 100^{th} percentile, 8 cm at the 95^{th} percentile, and 1 cm at the 67^{th} percentile for the horizontal position. These initial values are significantly lower compared to the other frequency combinations, as well as to the float solution. This reduction in initial errors was previously reported in the literature, for example by Li et al. (2020), who reported improvements up to 42% with a five-frequency solution compared to three- and four-frequency solutions when fixing widelane and extra-widelane ambiguities for Galileo. As a consequence of the decrease in position errors with the addition of frequencies, both convergence time and RMS errors were also reduced with the addition of more frequencies, as shown in Figure 6.

Figure 6 summarizes the statistics for the results shown in Figure 5. The convergence time shown here is defined as the time required to reach and remain below 10 cm horizontal error until the last epoch of dataset processing. Similar performance is noted for all float solutions, whether in terms of convergence time or RMS error. This behavior is expected as was previously discussed, as the additional frequencies have only limited effect on the float solutions; and most of their effect is in improving ambiguity fixing. It can be noted that fixing ambiguities leads to a 73% improvement in convergence time for the dual-frequency solution at the 100^{th} percentile. This significant decrease in convergence time demonstrates the importance of ambiguity resolution in improving PPP performance and its capacity to perform comparably to RTK. The significant decrease is emphasized when adding more frequencies to the processing procedure. The addition of one frequency decreases the convergence time by one minute to reach two minutes; the addition of two frequencies decreases the convergence time further to reach one minute. Keep in mind that the measurement rate in this case is 30 seconds, therefore the statistics are presented in 30 second increments. These improvements in convergence time are reached by determining the average of all datasets. At both the 95^{th} and 67^{th} percentiles, the all-frequency solution converges instantaneously, with convergence time values of zero (i.e., convergence that occurs with the very first epoch of measurements). The improvements from the additional frequencies are attributed to the fact that more ambiguities per satellite were included in the fixing process. Given that the ambiguities from the same satellite correlate with one another, the inclusion of additional ambiguities provides more information to mLAMBDA. This facilitates improved ambiguity fixing and leads to improvements in PPP-AR performance (Geng & Bock, 2013; Naciri & Bisnath, 2021a 2021b). These results highlight the possibility of instantaneous centimeter-level performance using PPP without additional local augmentations. This possibility has also been shown by Laurichesse & Banville (2018), who reported that dual-frequency GPS and four-frequency Galileo measurements could be used to demonstrate instantaneous centimeter-level positioning using only global products. The authors demonstrated that these results could be obtained using a limited sample size using one receiver in Australia over one hour while making use of global ionospheric maps (GIM) as a priori constraints for the ionosphere. In the current manuscript, the study is broadened to include more signals and constellations and demonstrate the performance of a global distribution of receivers and its consistency across multiple stations.

Similar conclusions can be reached from an assessment of the RMS error statistics. As can be seen from the figure, the RMS error is comparable between all float solutions; ambiguity resolution improves the RMS error while adding more frequencies decreases the RMS error. It is interesting to note that the all-frequency solution has an overall horizontal RMS error that is approximately half that of the dual-frequency solution. This significant decrease can be attributed to the errors at the initial stages of processing which are decreased with additional frequencies; this leads to a smaller overall RMS error and a reduced convergence time. As expected, ambiguity resolution has less of an effect on the vertical component than it does on the horizontal component. The vertical component is harder to estimate in general as a consequence of the satellite distribution around the antenna.

To obtain a better understanding of the performance shown in Figure 5 and Figure 6, Figure 7 highlights the 95* ^{th}* percentile performance of each station. Each bar is based on all seven days of processing from each station and is computed for each station as the convergence time of the solution that is generated by taking the epoch-by-epoch average of the 95% lowest horizontal errors. The performance metric of interest is the convergence time to a 10 cm horizontal error.

Figure 7 corroborates the observations presented by the overall average results. Some variation in performance between stations can be observed with respect to the float solution, as stations converge within 6.5 to 15 minutes. This finding is in accordance with previous results, as the average float solution regardless of station converges within 9.5 minutes and thus falls between the maximum and minimum performance of stations. Fixing ambiguities on two frequencies leads to a noticeable decrease in convergence time, as most stations converge within less than five minutes, and many converge within two minutes or less, which is the average convergence time based on all datasets shown in Figure 6. Finally, fixing ambiguities on as many as four frequencies reveals the possibility of instantaneous convergence for numerous stations. As shown, 63% of stations converge instantaneously, and 89% do so within one minute or less. This means that only three of the stations required more than one minute (three epochs) to converge. This very promising result is consistent with the average results shown in Figure 6 in which the average solution converged instantaneously. The dual- and triple-frequency results were omitted from the figure to maintain clarity. As anticipated, the results from the mix of dual and triple-frequency data fall in between the yellow and red lines shown in the figure. Collectively, these results demonstrate that modern PPP-AR can be used for applications requiring sub-decimeter-level horizontal (95%), and near-instantaneous positioning/navigation.

### 4.2 Benchmark Against Real-Time Kinematics (RTK)

The current subsection performs a more rigorous comparison against RTK-like performance through: 1) the use of a more strict convergence threshold and 2) an analysis of an epoch-by-epoch solution and the comparison of its performance with RTK-like requirements. Figure 8 summarizes the performance for each station by analyzing the convergence time to 2.5 cm rather than 10 cm, as was shown in Figure 7. This approach was selected to analyze the possibility of achieving the (N) RTK-like requirement of near-instantaneous 2.5 cm horizontal error at 1*σ* (Bisnath et al., 2013).

The findings shown in Figure 8 reveal that the requirement for instantaneous convergence to 2.5 cm is achieved by many of the stations. While all float solutions require at least 16.5 minutes to reach 2.5 cm at each station, both the dual-frequency, and the dual, triple, and quadruple-frequency (referred to as “all freq.”) fixed solutions result in substantial improvements. As shown, 81% of the stations converge instantaneously using all frequencies, compared to only 7% of the fixed dual-frequency solutions. Furthermore, 89% of the all-frequency stations converge within one minute, compared to 37% of the stations with fixed dual-frequency solutions. These improvements highlight the importance of the additional frequencies to reach the few centimeter level of error and to maintain that level of accuracy. While fixing ambiguities on two frequencies helps tremendously in the effort to achieve rapid centimeter positioning, the inclusion of more frequencies allows for more reliable fixing in the initial processing stages, and thus leads to instantaneous centimeter-level positioning.

A further comparison of the multi-constellation multi-frequency PPP-AR performance to RTK-like performance was performed by testing the proposed model using an epoch-by-epoch filter. PPP has traditionally used sequential estimation filters, in which information from previous epochs is used in epochs that follow. The use of sequential filters introduces the constraint of the ambiguities as constants and permits the use of previous information and satellite movement to estimate parameters more accurately. However, given the results shown in Figure 8, which document that instantaneous centimeter-level positioning is achievable, it is interesting to consider the possibility of an epoch-by-epoch solution in which each epoch is treated independently and previous information is not used. This type of solution would rely on additional constellations and frequencies to fix ambiguities instantaneously to achieve the required accuracy. Similar work was done by Geng & Guo (2020) in which extra-widelane (EWL) and widelane (WL) ambiguities were fixed for Galileo and BeiDou-3 to achieve single-epoch sub-decimeter positioning. The novelty of this work is the fact that these findings go further by presenting the possibility of centimeter-level positioning by fixing all ambiguities (up to four) for GPS, Galileo, and BeiDou-2/3. The availability of single-epoch precise solutions would mean more robust PPP, as earlier faulty epochs would not be used for current epochs in case of, for example, undetected cycle-slips, or high-level multi-path. The distribution of the horizontal errors resulting from the epoch-by-epoch processing is shown in Figure 9 at both the 100^{th} and 67^{th} percentiles.

The results from Figure 9(a) show that single-epoch, decimeter-level positioning is possible when fixing ambiguities on as many frequencies and constellations as are available. The mean horizontal error for the float dual-frequency results is 46 cm, compared to 44 cm for the fixed dual-frequency results and 12 cm for the fixed all-frequency results. These mean values, along with the shape of the graphs, show that single-epoch fixing in the dual-frequency mode does not lead to significant improvements in performance. In fact, fixing in the dual-frequency mode appears to degrade many of the float solutions, as the means are comparable but the standard deviation of the fixed results is 40 cm, compared to 26 cm for the float result. Therefore, single-epoch processing in dual-frequency mode is not an option to achieve higher precision positioning, particularly because fixing ambiguities in a single epoch is a difficult task to perform. However, as described in the previous subsection, adding more frequencies helps to achieve more reliable fixing, as shown in Figure 9. The mean of the all-frequency fixed distribution is approximately four times lower than the float and dual-frequency fixed solutions, with a smaller standard deviation (18 cm). This is highlighted by the fact that although 96% of all epochs are below half a meter, 66% and 61% of all epochs are below a decimeter and 2.5 cm, respectively. This performance is noticeably better than the other solutions and is attributed to the benefits provided by the additional frequencies toward more reliable fixing.

The results shown in panel (a) were based on all epochs of processing from all stations. Results focused on the 67^{th} percentile (1*σ*) are highlighted in panel (b), in which the benefits of the all-frequency solution are clear. Conclusions from panel (a) still hold, as fixing ambiguities on two frequencies leads to worse performance over many epochs, and the all-frequency solution has lower mean and standard deviation. In this 1*σ* case, 99.9% of epochs are below half a meter, 86% are below a decimeter, and 81% are below 2.5 cm. These results show that PPP is capable of RTK-like performance (i.e., centimeter-level epoch-by-epoch positioning) in the context of an epoch-by-epoch solution. These results are in agreement with and expand on the results shown by Laurichesse & Banville (2018) in which centimeter-level solutions were achieved with dual-frequency GPS and quadruple-frequency Galileo based on an epoch-by-epoch solution generated using a sample dataset. Results presented in this earlier publication also predicted that adding more constellations and frequencies will help make this level of performance a reality for more datasets. This prediction is realized by the findings presented here, as performance that exceeds 2.5 cm can be achieved for 81% of the epochs (1*σ*) using single-epoch processing with as many as four frequencies and constellations is shown. This level of performance is believed to be comparable to that typically achieved using RTK (i.e., epoch-by-epoch centimeter-positioning) and is a tremendous improvement over typical dual-frequency PPP results.

As shown in Figure 9(a) the all-frequency solution has a spike at the left-most side of the x-axis; this is followed by a Gaussian curve with a smaller magnitude. In order to understand this behavior, the station-specific percentage of epochs below 2.5 cm is plotted as shown in Figure 10. Both the 100^{th} and 67^{th} percentiles are included on different panels.

Figure 10(a) documents the percentage of epochs from each station that have a horizontal error less than 2.5 cm. Each bar is based on one week of station data. This metric (i.e., horizontal errors below 2.5 cm in each epoch) is equivalent to typical performance expected from RTK, in which a high percentage of epochs is expected to reach this degree of accuracy instantaneously. For example, Bock et al. (2000) showed that 16.8% of epochs fall outside of the 5.7 cm, 4.3 cm, and 38.4 cm outlier bounds in north, east and up directions, respectively given a 37 km baseline. This means that 83.2% of the epochs are within the defined intervals. Figure 10 shows that the float solutions are unlikely to reach instantaneous centimeter-level accuracy, as only 0.4% of epochs at most have an error below 2.5 cm for all stations. Ambiguity resolution significantly increases the percentage of stations with error below 2.5 cm in a dual-frequency mode, as these stations have 29% of of their epochs below 2.5 cm on average, including a minimum of 9.5% of the epochs for station CHPG and a maximum of 86.2% of the epochs for station ARHT. These results are improved further when considering all frequencies, as the mean number of epochs increases to 61%; the minimum percentage of converged epochs is reached by station PIE1 at 33.8% and a maximum of 96.9% is reached by station ARHT.

The average number of epochs below 2.5 cm error for both dual-frequency and dual, triple, and quadruple-frequency results corroborate the results shown in Figure 9(a). As shown, approximately 30% of the epochs are below 2.5 cm in the dual-frequency mode as compared to the approximately 60% observed with all frequencies. Therefore, the results included in Figure 10(a) can be used to understand the distribution shown in Figure 9(a). One noticeable point is the very strong performance exhibited by station ARHT in which 97% of the epochs are below 2.5 cm. Many of the epochs from this station drive the initial peak shown in Figure 9. However, this station’s epochs are not the only source of this peak; many other stations have more than 60% of their epochs below this strict threshold.

An assessment of the 67^{th} percentile performance shown in Figure 10(b) helps to gain further insight into the results, as lower percentiles eliminate some of the most poorly-performing epochs. The figure documents the impressive performance that can be achieved by fixing ambiguities on as many frequencies from as many constellations as available. These results show that many stations (18 of 27) can have more than 80% of their epochs within a 2.5 cm error using only a single epoch of data. On average, 80.4% of the epochs are below 2.5 cm error, which is consistent with the results shown in Figure 9(b). These results are impressive, as they highlight the level of single-epoch performance that can be achieved by adding more frequencies while fixing their ambiguities. These results show that it is possible to achieve instantaneous, single-epoch precise positioning. However, certain stations clearly perform better than others. To provide additional insight into the differences in performance between these stations, the findings shown in Figure 11 document the time series of ambiguity-fixed horizontal errors for three stations based on their performance as shown in Figure 10. Shown here is station ARHT, which is an example of a highly-performing station, station GOPE, which is an example of an average-performing station, and station BOGT, which is representative of a poorly-performing station. It should be noted that the percentages shown in Figure 11 are based on all epochs and correspond to percentages at the 100^{th} percentile. Note that all three solutions correspond to the “all-frequency” processing discussed in previous figures. This figure shows that more than 90% of the epochs in stations ARHT and GOPE have their ambiguities fixed correctly with the solution converging in every epoch, despite the fact that the epochs were processed independently. Additionally, the same percentages can be seen for both the 10 and 2.5 cm error thresholds for station ARHT due to the fact that fixing the ambiguities correctly leads to a centimeter-level position solution. However, station BOGT has many incorrect fixes. Only 49.7% of the epochs are below one decimeter, with the remaining epochs at a few decimeters. These decimeter-level epochs are similar to a single point positioning solution, as information from previous epochs is not used.

More research is needed to understand why some stations, (e.g., ARHT), perform better than others (e.g., BOGT). This information is critical in the attempt to bring most results up to a quality that is similar to station ARHT. Doing so would lead the way to more robust PPP. Attempts were made in this research to understand the different behaviors of various stations. For example, different types of receivers and antennas were assessed to verify that good or poor performance was not due to the antenna or receiver types; no correlations were found. Additionally, the locations of the receivers were evaluated to eliminate the possibility that a better or worse performance might be due to a different satellite view. However, many adjacent stations were found to have different levels of performance, thus ruling out this possibility. An analysis based on the number of satellites in view and of the PDOP was also performed; these findings revealed no direct correlations with the more poorly-performing stations. Other correlations were studied but none were found; further analysis will be required. For instance, in this manuscript, the same level of pseudorange and carrier-phase noise was assumed across all stations. This may not be a realistic assumption, as stations use different receiver and antenna combinations and are located in different environments, which are both factors that might affect the contributions of noise to their measurements. Additionally, the same weight is given to measurements on different frequencies, which may be improved with signal-specific values. Therefore, efforts toward fine-tuning the stochastic model and adapting it to individual stations might be expected to reduce the difference in performance between stations even further and thus lead to improvements in the consistency of the results.

## 5 CONCLUSION AND FUTURE WORK

For many years, RTK has been the standard for precise positioning because of its ability to achieve instantaneous centimeter-level positioning for short baselines or within GNSS networks. However, RTK relies on information from local reference stations with known coordinates. These local reference stations are not always available, as they are costly to maintain. PPP has superseded RTK in areas without RTK coverage, as it does not require additional infrastructure other than a sparse worldwide network of receivers. Since PPP relies only on the use of global corrections which can be transmitted either through Internet, satellites, or from the constellations themselves, PPP can be used anywhere in the world where GNSS signals can be tracked. The only caveat is that PPP may require tens of minutes to reach acceptable sub-decimeter accuracy. The work presented here uses the current multi-constellation multi-frequency context to provide as many measurements as available to PPP, which is a measurement-dependent technique. These additional measurements, together with ambiguity resolution, permit PPP to achieve RTK levels of performance without relying on local information. To perform such processing, a novel flexible uncombined quadruple-frequency decoupled clock model is derived and used to process data from receivers that are distributed around the world.

The results presented in this research show that instantaneous centimeter-level positioning using PPP is possible, as the average solution based on 1448 datasets converged instantaneously and remained below 10 cm horizontal error at both 2 *σ* and 1*σ*. These solutions were computed by taking the epoch-by-epoch average across all datasets after eliminating the 5% and 23% of the datasets with the highest error rates. This permitted the assessment of the average performance after eliminating outlier results. This impressive performance was confirmed by individual stations. On average, 63% of the stations converged instantaneously, while 89% converged within one minute when using all frequencies from all constellations. These results show that PPP can be used to achieve instantaneous sub-decimeter accuracy provided that ambiguities are fixed on as many frequencies and constellations as are available. These results are even more impressive when the thresholds were tightened, as it was found that 81% of the stations converged instantaneously below 2.5 cm error instantaneously at 1*σ*, and that 89% converge below 2.5 cm error within one minute at 1*σ*. However, one of the strengths of RTK is its ability to provide centimeter-level epoch-by-epoch positions for short to medium baselines. These epoch-by-epoch solutions are generated with the proposed multi-frequency multi-constellation PPP solution. Using this approach, 66% and 61% of all epochs were below one decimeter and 2.5 cm, respectively. This was achieved by processing data from every epoch individually. This is all the more impressive, as only 0.4% of epochs would ordinarily converge below a decimeter for a typical all-constellation dual-frequency float solution. Furthermore, it was showed that, on average, 60% of the epochs from individual stations were below 2.5 cm error in single-epoch processing. Remarkably, 96% of the epochs from station AHRT instantaneously converged to 2.5 cm. These results show that the 1*σ* epoch-by-epoch solutions for 18 of 27 stations have at least 80% of their epochs within 2.5 cm. Collectively, these results demonstrate the impressive RTK-like performance that can be achieved when leveraging many of the constellations and signals currently available to GNSS users.

These results represent drastic improvements in performance and highlight the benefit of leveraging all the constellations in orbit and all the frequencies broadcast, as they show that PPP can achieve instantaneous centimeter-level accuracy without needing any local information or infrastructure. These results demonstrate that it is possible to achieve centimeter-level positioning anywhere in the world using only a single receiver together with global satellite corrections. Additional work will be needed to generate a larger understanding of the differences in performance between stations for epoch-by-epoch processing with the goal of having all stations perform as well as station ARHT. Furthermore, only three frequencies were used for BeiDou-3. It would be very interesting to include the additional two frequencies currently broadcast by this constellation. Future work will also include an analysis of system performance using corrections transmitted by Galileo’s high accuracy service, as this introduces the possibility of having corrections transmitted by constellations themselves. This will allow users to obtain high-precision position information via signals broadcast by constellations alone.

## HOW TO CITE THIS ARTICLE

Naciri, N., & Bisnath, S. (2023) RTK-quality positioning with global precise point positioning corrections. *NAVIGATION, 70*(3). https://doi.org/10.33012/navi.575

## CONFLICT OF INTEREST

The authors declare no potential conflict of interests.

## ACKNOWLEDGEMENTS

The authors would like to thank the International GNSS Service (IGS) for the observational data, GeoForschungsZentrum (GFZ) for the satellite orbit and clock products, the Centre National d’Etudes Spatiales (CNES) for the satellite code and phase biases, and the Natural science and Engineering Research Council of Canada (NSERC) and York University for research funding.

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