Two-epoch centimeter-level PPP-RTK without external atmospheric corrections using best integer-equivariant estimation

Ambiguity resolution enabled precise point positioning (PPP-AR or PPP-RTK) without atmospheric corrections requires the user to estimate tropospheric and ionospheric delay parameters. The presence of the unconstrained ionosphere parameters impedes fast and reliable ambiguity resolution, so a time-to-first-fix of around 30 min for GPS-only solutions is generally reported, which can, to some extent, be reduced when combining multiple GNSS. In this contribution, we investigate the capabilities of almost instantaneous PPP-RTK, using only a few observation epochs at a sampling interval of 30 s, with the ionosphere-float model. The considered key elements are (a) the MSE-optimal best integer-equivariant estimator, (b) a combination of dual-frequency GPS, Galileo, BDS, and QZSS, (c) an area with good visibility of BDS and QZSS, and (d) a proper weighting of the PPP-RTK corrections. We provide a formal and simulation-based analysis of kinematic and static PPP-RTK with perfect, i.e., deterministic, clock and bias corrections as well as corrections computed from only a single reference station. The results indicate that, on average, one can expect centimeter-level positioning results with just slightly more than two epochs already with single-station corrections. This is confirmed with real four-system GNSS data, for which the availability of two-epoch centimeter-level horizontal positioning results is 99.7% during an exemplary day.


Introduction
For precise point positioning (PPP), a global navigation satellite systems (GNSS) user requires precise satellite orbit and clock products. Static PPP is capable of reaching a centimeter-level positioning accuracy but only with a convergence time of several hours (Zumberge et al. 1997;Kouba and Héroux 2001;Bisnath and Gao 2008). When also provided with satellite phase biases, a single-receiver user can recover the integer property of the carrier phase ambiguities and reduce the long convergence times through ambiguity resolution (AR), known as PPP-RTK or PPP-AR (Wübbena et al. 2005;Laurichesse et al. 2009;Mervart et al. 2008). Further strategies and methods for PPP-RTK have, for instance, been formulated and demonstrated in Collins (2008), Ge et al. (2008), Bertiger et al. (2010), Teunissen et al. (2010), Zhang et al. (2011), or Geng et al. (2012. An overview and comparison are provided in Teunissen and Khodabandeh (2015), and the interoperability of different PPP-RTK corrections is shown in Banville et al. (2020) for products of the International GNSS Service (IGS).
Successful AR instantly leads to centimeter-level PPP-RTK results, but fast and reliable AR still remains a challenge, mainly due to the presence of ionospheric delays. A time-to-first-fix (TTFF) of about 30 min with 1 Hz GPS data is reported in Geng et al. (2011), and a very similar value in Zhang et al. (2019) with 30-s data. Similar results of 34/22 min with 30-s GPS data for kinematic/static PPP-RTK are obtained in Li and Zhang (2014), which are reduced to 20/16 min when integrating GLONASS data. Geng and Shi (2017) reported convergence times of 25 and 6 min with GPS and GPS + GLONASS using partial AR, and Li et al. (2018) showed that the convergence time of static PPP-RTK can be further reduced to around 9 min when also including BDS and Galileo data. In Li et al. (2020), Galileo PPP-RTK 12 Page 2 of 11 with two to five frequencies is analyzed, with resulting convergence times between 22 and 15 min with 30-s data. In a simulation study, Psychas et al. (2020) show that sub-decimeter horizontal positioning errors can be expected within 25 min for GPS and 6.5 and 4.5 min for dual and triple frequency GPS + Galileo + BeiDou with partial AR and deterministic PPP-RTK corrections. Different criteria used in the above studies when defining the TTFF or the convergence time make it difficult to compare those numbers, but a general trend to faster solutions when combining GNSS constellations or using advanced AR strategies is clearly visible.
Faster centimeter-level positioning solutions can be obtained when also ionospheric corrections are provided to the user. However, as demonstrated in Psychas et al. (2018), these have to be at the level of a few centimeters to noticeably accelerate reliable AR. This requirement is not met by current global ionosphere solutions, such as the total electron content maps provided by the IGS, which have an uncertainty of several decimeters for GNSS frequencies (Brack et al. 2021b). Instantaneous PPP-RTK results using (interpolated) ionospheric corrections from one or more nearby reference stations are reported in Teunissen et al. (2010), Banville et al. (2014), and Psychas et al. (2022), however, at the cost of requiring rather dense local or regional GNSS networks.
An alternative to fixing the ambiguities to integers is to use the best integer-equivariant (BIE) ambiguity estimator (Teunissen 2003). The resulting position estimates are MSE optimal, and a closed-form expression for normally distributed data is provided in Teunissen (2003). An extension of the BIE principle for elliptically contoured distributions is presented in Teunissen (2020), and a sequential scalar approximation is proposed in Brack et al. (2014). An evaluation of the BIE estimator based on simulations is given in Verhagen and Teunissen (2005), and its performance for multi-GNSS single-baseline RTK positioning is analyzed in Odolinski and Teunissen (2020).
In Brack et al. (2021a), it is demonstrated that when combining all available GNSS with a partial AR approach, one can reach centimeter-level single-baseline RTK results within 3.3 epochs of 30-s data when ionospheric delays have to be estimated. A similar positioning performance should also be feasible with PPP-RTK without atmospheric corrections, which is the topic of this contribution. We will make use of (a) the MSE-optimal BIE estimator, (b) a combination of dual-frequency GPS, Galileo, BDS2/3, and QZSS, (c) an area with good visibility of BDS and QZSS, and (d) a proper weighting of the PPP-RTK corrections. Simulated GNSS data in the area of Perth, Australia, will be used to show that with corrections from only a single reference station just slightly more than two epochs are required on average to reach centimeter-level static or kinematic PPP-RTK results, whereas using deterministic corrections often only one epoch is sufficient. This result will be confirmed with real GNSS data, where centimeter-level horizontal position estimates are obtained after two epochs with an availability of 99.7%, thus demonstrating that almost instantaneous PPP-RTK is feasible with the current GNSS constellations, even without any atmospheric corrections.

PPP-RTK observation model
The single-system undifferenced, uncombined code and carrier phase observations p s r,f and s r,f between receiver r and satellite s on frequency f can be modeled as with E[⋅] being the expectation operator, s r the geometric range between satellite s and receiver r, dt r and dt s the receiver and satellite clock offsets, r the residual zenith tropospheric delay (ZTD) at receiver r with the mapping function m s r , i s r the first-order slant ionospheric delay on the first frequency with the coefficient f = 2 f 2 1 depending on the wavelengths f , d r,f and d s ,f the receiver and satellite code biases, r,f and s ,f the respective phase biases, and a s r,f the integer phase ambiguity.
When processing the data of a GNSS receiver or network, not all parameters in (1) can be unbiasedly estimated due to rank deficiencies in the underlying system model. Estimable combinations of the parameters resulting in a full-rank model can be determined using S-system theory (Baarda 1973;Teunissen 1985) by constraining a minimum set of parameters.
We assume that precise orbital positions are available. We further assume that for the reference stations providing the PPP-RTK corrections, the coordinates are a priori known, so that the ranges s r can be removed from the model. A set of estimable parameters together with their definitions is given in Table 1, where the first receiver and satellite are chosen as pivot (cf. Odijk et al. 2016). In the case of a local network, the tropospheric mapping functions are (almost) identical so that the ZTDs can only be estimated relative to the pivot receiver, whose tropospheric slant delays are then included in the satellite clock corrections. This also holds if only one reference receiver is used. The PPP-RTK corrections provided to the user contain the estimates of the satellite clocks dt s , the satellite phase biases ̃s ,f , and the satellite code biases d s ,f for f > 2. On the user side, the observed-minus-computed observation equations follow from (1) by replacing the ranges s r with g s,T r Δx r , where g s r is the satellite-to-receiver unit vector, Page 3 of 11 12 and Δx r is the unknown increment of the user position from the initial value. The PPP-RTK corrections are applied to the user observations, and the corresponding parameters are removed. Other than that, the user parameters are the same as defined in Table 1 for r ≠ 1.
For the multi-GNSS model, we assume that all parameters except the user position Δx r and the ZTDs r or 1r are set up per constellation. This implies that the biases are assumed constellation specific. A calibration of inter-system biases on overlapping frequencies that would enable the use of a common pivot satellite is not considered (Odijk and Teunissen 2013;Paziewski and Wielgosz 2015). With the satellite phase biases applied, the user can recover the integer estimable double-difference ambiguities with the pivot receiver and use AR techniques.

Solving the PPP-RTK user positioning model
Let the integer ambiguity parameters of the PPP-RTK user model be collected in the vector a ∈ ℤ n and all real-valued parameters, including the position increment Δx r , in b ∈ ℝ p . We consider three estimators in this contribution, which are briefly introduced in the following.
The first one is the ambiguity float estimator b , for which the integer property of the ambiguities is disregarded. For short observation time spans, its precision is driven by the code data. The high precision of the carrier phases only starts to have an impact when multiple epochs with timeconstant ambiguities are combined.
The second estimator is the ambiguity fixed estimator b , for which all ambiguity parameters are fixed to integers. The integer least-squares (ILS) ambiguity estimator is given by Q̂a with â the float ambiguity solution and Q̂a its covariance matrix. It is optimal in the sense of maximizing the probability of correct ambiguity estimates (Teunissen 1999) and is efficiently implemented in the LAMBDA method (Teunissen 1995). If the ambiguities are resolved correctly, the positioning precision is immediately driven by the carrier phases, but wrong ambiguity estimates can lead to large positioning errors, so one will usually prefer the float solution b if the ambiguity success rate is too low.
The third estimator is the BIE estimator b (Teunissen 2003). For normally distributed data, the BIE ambiguity estimates are a weighted sum of integers and b is the conditional least-squares estimator assuming the ambiguities given by a . For computational reasons, the infinite set ℤ n in (3) has to be replaced by a finite set, which is chosen as the set of integers within an ellipsoidal region around â . Its radius is defined in the metric of Q̂a and is derived from a central Chi-squared distribution with n degrees of freedom and a significance level of = 10 −5 , see Teunissen (2005). The BIE solution b is MSE-optimal and unbiased within the class of integer-equivariant estimators, which also contains the float and fixed solutions b and b (Teunissen 2003). As the variances of the BIE estimates b are equal or smaller than those of the float and any admissible fixed or partially fixed solution, they can serve as a benchmark for analyzing the best possible performance of a GNSS model. This is done in the following in order to evaluate the limits of the convergence time of the PPP-RTK model that can be achieved with the current GNSS constellations.
( Table 1 Estimable parameter combinations with singlesystem undifferenced and uncombined code and carrier phase observations; the receiver and satellite with index 1 are chosen as pivot;  (2023) 27:12 12 Page 4 of 11

Experimental setup
The multi-GNSS PPP-RTK positioning capabilities of combined GPS (G), Galileo (E), BDS 2 + 3 (C), and QZSS (J) are analyzed with the ambiguity float, ILS and BIE solutions. The station PERT (Trimble Alloy) in Perth, Australia, acts as the user receiver, and the corrections are provided by the stations CUT0, CUTA, CUTB, and CUTC (all Trimble NetR9, distance to PERT 22.4 km) forming a local receiver array or by the individual station NNOR (Septentrio PolaRx5TR, distance to PERT 88.5 km). The number of visible satellites with an elevation angle greater than 10° during April 1, 2022, is shown in Fig. 1. The observations are weighted according to exponential noise amplification factors as defined in Euler and Goad (1991). The zenith-referenced standard deviations of the considered dual-frequency code observations are estimated from a different day of double-differenced code data from the local CUT* receiver array in a least-squares sense (Teunissen and Amiri-Simkooei 2008) and given in Table 2. For the tropospheric modeling, we use the global mapping function (Boehm et al. 2006) and the a priori corrections from the blind MOPS model (MOPS, 1999). The GFZ precise multi-GNSS orbit products are applied (Deng et al. 2017).
The PPP-RTK strategy used in this contribution is implemented as follows: The above-mentioned PPP-RTK corrections are computed on an epoch-by-epoch basis from either a single reference receiver (CUT0 or NNOR) or from the local CUT* receiver array. AR techniques are not applied on the network side. In the simulations, also perfect, i.e., deterministic, corrections are considered. The user station PERT then applies these corrections together with their uncertainty as given by their covariance matrix. Neglecting the uncertainty of the corrections can strongly degrade the position performance (Psychas et al. 2022). No delay is assumed between the generation of the PPP-RTK corrections and the user positioning. The float solution is computed with a recursive least-squares implementation, and the ILS and BIE ambiguity solutions ǎ and a are computed anew in each epoch. We consider both static and kinematic positioning, where for the latter, the coordinates are assumed to be completely unlinked in time. The user ambiguities and receiver phase biases are assumed time-constant, and the ZTD is modeled as a random walk with a process noise of 2 mm∕ √ h . All other parameters including the ionospheric delays are assumed unlinked in time, so that the results do not depend on the ionospheric activity.

Formal and simulation analysis
The average ambiguity float positioning precision against the number of epochs is shown in Fig. 2 for the east component, where the estimation is started every 10 min during April 1, 2022. We can clearly see the benefit of combining multiple systems, so that, for instance, the kinematic precision with single-station corrections reaches 30 cm after around 15 epochs for GPS, 6 epochs for GPS + Galileo, and 4 epochs when also including BDS and QZSS. The static counterparts    converge slightly faster. Khodabandeh and Teunissen (2015) demonstrate that PPP-RTK with single-station corrections is equivalent to single-baseline RTK positioning. Applying the theoretically assumed perfect corrections instead of the single-station corrections translates to having no measurement noise at the reference receiver for single-baseline RTK positioning. Assuming a similar noise level at both receivers, the uncertainty of the PPP-RTK user observations with the corrections applied is then essentially reduced to half, so that the same precision values are obtained faster. When extending the single reference station to an array or network of receivers, the user positioning precision will be between these two shown extremes. In all presented cases, however, a centimeter-level precision within at most a few epochs cannot be expected.
In comparison, the ambiguity fixed precision values in the first epoch are already between 9 mm for the GPS-only case with single-station corrections and 3 mm for the foursystem case with perfect corrections. These values are the conditional standard deviations assuming the ambiguities known and do not reflect the uncertainty of the ambiguity estimators, but rather the highest precision that can possibly be reached with any ambiguity estimator. A formal measure of the strength of the GNSS model for ambiguity resolution is the ambiguity dilution of precision (ADOP, Teunissen 1997). The average ADOP values of the above PPP-RTK examples are shown in Fig. 3. Again, we see an improvement when combining systems or increasing the quality of the corrections. Odijk and Teunissen (2008) found that an ADOP ≤ 0.12 generally allows for reliable ILS ambiguity resolution with a failure rate of less than 0.1%. This seems promising for the four-system cases with perfect corrections, but when we look at the actual average TTFF with ILS and a failure rate constraint of 0.1% in Table 3 (the integer bootstrapping failure rate as a tight upper bound is used, see Verhagen et al. 2013), we see that on average still around 15 epochs are required in this case, which is far from the envisioned almost instantaneous solutions. The increase of the convergence time of the four system solution compared to GPS + Galileo is attributed to more frequently rising satellites. As the BIE estimator does not fix the ambiguities to integers, the concept of a success rate does not apply. In the following, it is investigated how this increased model strength as measured by the ADOP translates into high positioning precision with the BIE estimator.
The    Fig. 5; see also Verhagen and Teunissen (2005), Brack (2019), or Odolinski and Teunissen (2020) for similar results. While the BIE solution has the smallest RMS error, the ILS solution has a much higher probability of very small positioning errors but is also more likely to result in large errors. Interestingly, the ILS success rate of this example is 26.6%, implying that correct ambiguity estimates ǎ are not necessarily required for centimeterlevel errors, which occur at about 50%. Finally, the 3D RMS positioning errors of the above example relative to the float solution against the number of epochs are shown in Fig. 6; see also Brack (2019) or Odolinski and Teunissen (2020). The fixed solution is worse than the float solution in the first epochs due to the very low ILS success rate and eventually, with increasing success rate, converges to the dashed line, for which the ambiguities are assumed known. The BIE estimator always has the minimum RMS error. Its results are close to the float solution for very low ambiguity precision in the first epochs and close to the fixed solution for high ambiguity precision, which was proven in Teunissen (2003). The capabilities of the three estimators for kinematic and static PPP-RTK are analyzed by means of the average simulated RMS positioning error, shown for the east component in Fig. 7, and by means of the average probability of obtaining an absolute positioning error of less than 3 cm for the horizontal components and 15 cm for the vertical component in Fig. 8. Again, solid and dashed lines represent the setup with single-station and deterministic corrections, respectively. The float results confirm the formal analysis in Fig. 2, i.e., fast centimeter-level results cannot be expected. The BIE results always lead to the smallest RMS positioning errors. For strong positioning models in which centimeterlevel results are actually possible, the RMS errors with ILS are very close to the BIE results. Almost instantaneous centimeter-level results are only obtained with the four-system combination, where the RMS error of the BIE estimator with single-station corrections is at the sub-decimeter level using two epochs and the sub-centimeter level using three epochs. With deterministic corrections, about one fewer epoch is required. This result is confirmed by the complementary availability values in Fig. 8, showing a very high availability of precise positioning solutions with two to three epochs with single-station corrections and one to two epochs with perfect corrections. This figure also demonstrates again that the ILS solution can have a higher probability of very precise solutions than the BIE solution, cf. Figure 5. That is, although the RMS errors of the BIE and ILS estimators are very similar, their characteristics are still quite different.

Real-data analysis
Judging from the simulation results in the previous section, we can expect centimeter-level four system PPP-RTK results in generally not more than three epochs, quite often with one or two epochs, depending on the quality of the corrections. This is now verified with real GNSS data. We focus on the kinematic case, but a big benefit of the static case can anyway not be expected when only a few epochs are used.
Positioning solutions are initialized at every 30 s observation epoch.
With real data, it is obviously not possible to generate perfect PPP-RTK corrections that can be assumed deterministic. To demonstrate the improvements in the positioning capabilities with corrections of higher quality, we compare the results using single-station corrections from CUT0 and corrections from the four-station array CUT0/A/B/C. The same set of satellites is available in both cases. The resulting empirical RMS positioning errors of the three estimators are given in Table 4 for different numbers of epochs. The optimality of the BIE estimator is also confirmed with the real data, where for five and ten epochs, the fixed solution is essentially equally good in terms of the RMS error. Already by using a local array to compute the corrections instead of a single station, the PPP-RTK performance can be considerably improved, with an absolute improvement of around one meter for the vertical components in the first epoch and a relative improvement of the BIE RMS errors between 7.5% and 77.3%. The obtained empirical RMS positioning errors cannot keep up with the simulations, which promised subcentimeter values already for three epochs. The reason is that the NetR9 receivers at the CUT* stations do not track   BDS satellites with a PRN larger than 30, so the number of visible satellites is clearly reduced. In a second example, we therefore make use of the station NNOR tracking all BDS satellites to compute the PPP-RTK corrections. The empirical RMS positioning errors are given in Table 5. The horizontal RMS errors using the optimal BIE estimator are at the decimeter level with one observation epoch, at the centimeter level with two observation epochs, and at the sub-centimeter level with three or more epochs, confirming that two-epoch centimeter-level PPP-RTK without atmospheric corrections is possible.
The positioning errors of the east, north, and up component with one and two observation epochs are shown in Fig. 9. The corresponding empirical rates of a positioning error of less than 3 cm for the horizontal components and 15 cm for the up component are given in Table 6, where for the values in the parenthesis only the horizontal components are considered. With one epoch, the horizontal positioning error of the ILS (red) and BIE (blue) solutions is less than 3 cm in 97.6% and 87.6% of the cases. With two epochs, only six and nine of the 2,879 epochs exceed these limits, and with three epochs, this limit is always met with ILS and exceeded twice with the BIE estimator. The results in Table 6 confirm the simulation results in the previous section: While the BIE solutions are RMS optimal, see Tables 4 and 5, ILS often leads to a higher availability of very small positioning errors.

Conclusions
The user performance of PPP-RTK with the current GNSS constellations and without atmospheric corrections was analyzed. In the literature, often convergence times of several tens of minutes are reported for this model in order to reach centimeter-level positioning results. This contribution focused on the feasibility of almost instantaneous precise solutions. The main conclusions can be summarized as follows: The BIE estimator provides MSE optimal positioning results within the class of integer-equivariant estimators and is in this sense superior to the ambiguity float and any admissible ambiguity fixed solution. Its results can therefore be seen as benchmark results to evaluate the highest achievable performance of a given GNSS model. To this end, a simulation analysis of kinematic and static PPP-RTK with perfect corrections and single-station corrections was conducted. It showed that a high availability of almost instantaneous centimeter-level results with one to three observation epochs can only be expected when combining all available   systems, namely GPS, Galileo, BDS, and QZSS, but is generally possible. This result was confirmed with real data recorded at the station PERT, with corrections provided by the station CUT0, a local array of four receivers around CUT0, or the station NNOR. Although perfect PPP-RTK corrections can, of course, not be provided, the results showed how the performance of the PPP-RTK user model could already be improved when using a local array of receivers instead of a single receiver, thereby improving the empirical BIE RMS position errors by 7.5% to 77.3%. Using PPP-RTK corrections from the station NNOR, horizontal positioning errors of less than 3 cm are obtained with the BIE estimator in 87.6% of the cases with one observation epoch, and in 99.7% of the cases with two observation epochs. Two-epoch centimeter-level horizontal PPP-RTK results without atmospheric corrections are, therefore, indeed feasible. While this precision cannot be reached for the vertical component, its empirical two-epoch RMS error using the BIE estimator was reduced to 8.5 cm as compared to 14.2 cm with ILS.
It has to be stressed, however, that with the current constellations such short observation spans require a good satellite visibility. For instance, BDS satellites with a PRN larger than 30 are not tracked by the CUT* stations. As a consequence of removing these satellites, centimeter-level RMS positioning errors with two epochs were no longer possible, see Table 4. In view of the developments of the GNSS constellations in the previous years, however, the number of available GNSS satellites can be expected to increase further in future so that similar results will also be possible on a global scale. Further improvements toward instantaneous results can be expected when observations on additional frequencies are included. Psychas et al. (2021) concluded that the frequency separation is more important than the number of frequencies in terms of the ambiguity resolution performance, so that, for instance, adding E6 observations to a Galileo E1/E5a solution is more beneficial than adding both E5b and E5. As the computational burden of the BIE estimator increases with the ambiguity dimension, this result can prove very beneficial for extending the presented study to multi-GNSS solutions with three or more frequencies.
Author contributions The study was designed, and the data analysis was performed by AB. The manuscript was written by AB. All authors contributed to the discussion about the content and provided comments on the manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.