PPP-RTK considering the ionosphere uncertainty with cross-validation

With the high-precision products of satellite orbit and clock, uncalibrated phase delay, and the atmosphere delay corrections, Precise Point Positioning (PPP) based on a Real-Time Kinematic (RTK) network is possible to rapidly achieve centimeter-level positioning accuracy. In the ionosphere-weighted PPP–RTK model, not only the a priori value of ionosphere but also its precision affect the convergence and accuracy of positioning. This study proposes a method to determine the precision of the interpolated slant ionospheric delay by cross-validation. The new method takes the high temporal and spatial variation into consideration. A distance-dependent function is built to represent the stochastic model of the slant ionospheric delay derived from each reference station, and an error model is built for each reference station on a five-minute piecewise basis. The user can interpolate ionospheric delay correction and the corresponding precision with an error function related to the distance and time of each reference station. With the European Reference Frame (EUREF) Permanent GNSS (Global Navigation Satellite Systems) network (EPN), and SONEL (Système d'Observation du Niveau des Eaux Littorales) GNSS stations covering most of Europe, the effectiveness of our wide-area ionosphere constraint method for PPP-RTK is validated, compared with the method with a fixed ionosphere precision threshold. It is shown that although the Root Mean Square (RMS) of the interpolated ionosphere error is within 5 cm in most of the areas, it exceeds 10 cm for some areas with sparse reference stations during some periods of time. The convergence time of the 90th percentile is 4.0 and 20.5 min for horizontal and vertical directions using Global Positioning System (GPS) kinematic solution, respectively, with the proposed method. This convergence is faster than those with the fixed ionosphere precision values of 1, 8, and 30 cm. The improvement with respect to the latter three solutions ranges from 10 to 60%. After integrating the Galileo navigation satellite system (Galileo), the convergence time of the 90th percentile for combined kinematic solutions is 2.0 and 9.0 min, with an improvement of 50.0% and 56.1% for horizontal and vertical directions, respectively, compared with the GPS-only solution. The average convergence time of GPS PPP-RTK for horizontal and vertical directions are 2.0 and 5.0 min, and those of GPS + Galileo PPP-RTK are 1.4 and 3.0 min, respectively.


Introduction
The Precise Point Positioning (PPP) technique has been widely used in various fields, such as timing, meteorology, and low-earth-orbit satellite orbit determination (Zumberge et al. 1997). The biggest limitation of the traditional Global Positioning System (GPS) dual-frequency ambiguity-float PPP is the long convergence time of 30 min or more to achieve a position solution at the centimeter-level. To overcome the limitation, the PPP based on a Real-Time Kinematic (RTK) network with ambiguity-fixing and optionally external high-precision atmosphere corrections from regional reference stations is proposed. The PPP-RTK technique uses the state-space representation (Wübbena et al. 2005) to provide the users with the individual Global Navigation Satellite System (GNSS) related errors instead of the raw and/or corrected observations. It takes the full advantages of PPP and RTK and enables the integer Ambiguity Resolution (AR) of the PPP. Since then, the PPP-RTK technique has attracted a lot of attention in the GNSS community.
Several equivalent methods with different variants were proposed by Collins et al. (2008), Ge et al. (2008), and Laurichesse et al. (2009), and the satellite phase biases are provided as the corrections to the users to construct double-differenced ambiguities by using the S-system theory (Khodabandeh and Teunissen 2019). PPP-AR can significantly improve the positioning precision anywhere around the world with an extent of 30-60% but the contribution to shortening convergence time is marginal (Ge et al. 2008;Hu et al. 2020;Li et al. 2016). Furthermore, researchers use a regional troposphere model to enhance the PPP solutions. For example, Shi et al. (2014) introduced a method to determine the Optimal Fitting Coefficients (OFCs) of local troposphere models to augment the GPS PPP-AR. With the OFCs model, De Oliveira et al. (2017) tested the local troposphere model in a larger area over France with a maximum height difference of 1 651 m. The modeled Zenith tropospheric Wet Delays (ZWDs) present a precision of around 1.3 cm with respect to the final Zenith Tropospheric Delay (ZTD) products. In addition, Zheng et al. (2017) established a real-time tropospheric grid model based on the Global Pressure and Temperature 2 wet (GPT2w) to provide the ZTD with a precision of about 1.2 cm for the users in China to accelerate the BeiDou Navigation Satellite System (BDS) and GPS PPP.
In the studies mentioned above, the ionosphere is eliminated by forming the ionosphere-free combination or estimated in the Un-Differenced and Un-Combined (UDUC) model without any external constraint, which is the so-called ionosphere-float PPP-RTK model. On the one hand, the number of ionosphere parameters is much larger than that of the troposphere parameters. On the other hand, the ionosphere is highly variable with the time and area and hard to be precisely modeled in a larger area. The precision of the global ionosphere map is a few decimeters. Therefore, tens of minutes are needed to obtain a convergent solution for the ionosphere-float PPP-RTK model, as the unmodelled ionosphere is a major barrier to rapid ambiguity resolution.
Moreover, researchers studied the ionosphere-fixed PPP-RTK in which the ionospheric delays of a user are directly corrected by the interpolated values from its nearby reference stations. This model is mainly used in a small or medium-scale region with a inter-station distance of such as 50 km (Li et al. 2011(Li et al. , 2014, about 100 km (Psychas et al. 2019), and up to 200 km under a calm ionosphere condition . Although the fast, even single-epoch instantaneous PPP AR has been demonstrated achievable, it is not applicable in many areas where no dense CORS stations exist or the ionosphere conditions are severe.
The ionosphere-weighted model is a general model, from which the aforementioned two models can be produced. It becomes the ionosphere-float model if there is no prior ionosphere information, and the ionospherefixed model if the ionospheric corrections are precise enough. To make the PPP-RTK applicable for wider areas, the high-accuracy ionospheric delays must be interpolated at the user stations with the precise information. Teunissen and Odijk (2010) used the ionosphere standard deviation (STD) of 5 mm and 1 cm in two small-scale networks with inter-station distances of around 20 km and 60 km, respectively, in different locations and achieved the corresponding success rate of 99.8% and 91.8% for single-epoch ambiguity resolution. Using two networks with inter-station distances ranging about 50 km in different locations, and with the fixed ionosphere STD of 1 dm and 1.4 dm, Zhang et al. (2011) showed that the time-to-first-fix was about 5 min for both single-and dual-frequency PPP-RTK. Wang et al. (2020) comprehensively assessed six low-order surface models and one distance-based linear interpolation method and used them to generate interpolated corrections to conduct PPP-RTK applications and to develop proper interpolation methods for different network scales, terrain, and numbers of reference stations. For each ionosphere interpolation model, a fixed STD value is calculated. These fixed ionosphere STDs are mainly applicable in the case studies for the networks with a similar size.  empirically determined the variance of the ionosphere corrections by comparing the predicted ionospheric delay and the ambiguity fixed ionospheric delay at a user location, and achieved PPP ambiguity resolution within 3 min with the atmosphere corrections from the network whose inter-station distance is about 100 km. In addition, Nadarajah et al. (2018) used the distance lineardependent ionosphere STD for high-end geodetic receivers and low-cost single-frequency mass-market receivers to provide numerical insights into the role taken by the multi-GNSS integration in delivering fast and high-precision PPP-RTK solutions. In the study of Zha et al. (2021), an empirical distance exponential-dependent ionosphere stochastic model is used, and the ionosphere-weighted UDUC PPP-RTK is assessed with the network with an inter-stations distance of about 100 km. Fixed coefficients for the model are employed in their studies. These ionosphere STDs will change with the network size but are time-irrelevant and may be imprecise in the active ionosphere period. Furthermore, Li et al. (2020) evaluated the BeiDou-3 Navigation Satellite System (BDS-3), BeiDou-2 Navigation Satellite System (BDS-2), GPS stand-alone, and combined GNSS PPP-RTK performance based on three days of observations and showed that BDS can provide high-accuracy positioning services independently for the users in Europe. They calculated the interpolated ionosphere STD by comparing the ionosphere at each reference station with the weighted value. However, this STD may be vulnerable to gross error because it is calculated with few ionosphere observations. Wang et al. (2022) analytically studied the form of the variancecovariance matrix of ionosphere interpolation errors for both accuracy and integrity purposes, which considers the processing noise, the ionosphere activities, and the network scale. How to properly determine the a priori precision for the interpolated ionospheric delay still needs further investigation.
The goal of this study is to build a proper empirical stochastic model for the interpolated ionosphere in wide-area PPP-RTK. By the cross-validation of the slant ionospheric delays between reference stations from near to far, a distance-dependent and 5 min piece-wise linear function is proposed to model the interpolated ionosphere error RMS. Furthermore, the ionosphere precision map is generated to reflect the precision of the ionosphere correction for the region covered by the PPP-RTK service. It is expected to improve the wide-area PPP-RTK performance with over 100 km inter-station distance in which case the ionosphere interpolated error can no longer be ignored. This paper is organized as follows. First, the relevant PPP-RTK studies on ionosphere interpolation are reviewed. The next section formulates the mathematical model for PPP-RTK data processing at both server and user ends, addresses the methods applied to build a stochastic model which is a function of distance and time for ionosphere correction. Then the GNSS data and the experiment scheme are described, the troposphere and ionosphere precision map are shown, the PPP-RTK performance with our automatic stochastic model is compared with the fixed precision strategies, and the contribution of multi-GNSS on PPP-RTK is analyzed as well. Finally, the conclusions and remarks are presented.

Method
The UDUC PPP model can directly estimate and constrain the tropospheric and ionospheric delays, and is flexible to process single-, dual-, or multi-frequency observations. It is commonly used in PPP-RTK.

GNSS observation model
For a satellite s observed by a receiver r, the UDUC PPP equations of pseudorange P and carrier phase L are as follows (Zumberge et al. 1997): where i represents the carrier frequency ( i = 1, 2 ); ρ is the geometric distance from the satellite antenna phase center to the receiver phase center; c is the speed of light in vacuum; dt r and dt s are the receiver clock offset and satellite clock offset, respectively; m is the tropospheric mapping function; T r is the tropospheric zenith wet delay; is the ionospheric coefficient on frequency i; f is the frequency; I s r,1 is the slant ionospheric delay from the receiver to the satellite on frequency 1; i is the wavelength of frequency i; N s r,i is the integer carrier phase ambiguity; b r,i and b s i are the receiver and satellite pseudorange hardware delays on frequency i, respectively; B r,i and B s i are the receiver and satellite phase hardware delays on frequency i, respectively; e s r,i and ε s r,i are the observation noise of pseudorange and carrier phase, respectively. Other system error items, such as phase windup, Phase Center Offset (PCO) and Phase Center Variation (PCV), tide loading, and relativistic effect, are not presented in Eqs. (1) and (2), and assumed to be corrected by the corresponding model (Petit et al. 2010).
The precise products from International GNSS Service (IGS) are usually used in the PPP data processing to correct satellite orbit and satellite clock offset errors. The precise clock offset includes the satellite pseudorange hardware delay deviation. After corrected with the precise products, due to the correlation among receiver clock offset, satellite clock offset, ionospheric delay, and hardware delay parameters of carrier phase and pseudorange, Eqs. (1) and (2) are of rank deficiency and cannot be directly solved. The method of estimating the reparameterized parameters instead of the original ones is usually used to eliminate the rank deficit (Odijk et al. 2016(Odijk et al. , 2017. Therefore, the reparameterized UDUC PPP equation for PPP-RTK is (Zhang et al. 2019) (Table 1): where ρ s r is the geometric distance from the satellite antenna phase center to the receiver phase center corrected by the precise satellite clock and orbit as well as other error models; x, y, z are the receiver coordinates; c · dt r,IF 12 is the receiver clock offset after reparameterization; I s r,1 and N s r,i are the slant ionospheric delay and (1) phase ambiguity after reparameterization, respectively; α and β are the coefficients of ionosphere-free combination; D r,12 DCB and D s,12 DCB are the differenced code hardware delays between frequency 1 and 2 for receiver and satellite, respectively.
Furthermore, on both the server and user sides, correct ambiguity fixing is the key to besting the accuracy of estimated parameters. The satellite Uncalibrated Phase Delay (UPD) bias should be corrected to recover the integer property of the PPP ambiguities (Hu et al. 2020). In this study, the GBM integer clock provided by German Research Centre for Geosciences (GFZ) is used for PPP and ambiguity resolution (Deng et al. 2016). The widelane UPD is given in the header of the GBM precise clock file while the narrow-lane UPD is grouped into the satellite clock correction.
The server side adopts Eqs.
(3)-(4) to obtain the UDUC ambiguity-fixed PPP at each reference station to extract the precise atmospheric delay . In turn, the information on the accurate tropospheric and ionospheric delays can be obtained from the reference station. The troposphere is precisely modelled on the server side, then the coefficient and precision of the troposphere model are broadcasted to users. The slant ionospheric delays of the reference stations and their uncertainty are provided to interpolate the ionospheric corrections for the users. The detailed procedure for the troposphere and ionosphere modelling will be elaborated in the following two sub-sections. The atmospheric information can be used to enhance the positioning results on the user side, as the following equations: where the T u and I s u,i represent the tropospheric and ionospheric delays at user ends, respectively.
The high precision atmospheric delays generated on the server side based on the UDUC ambiguity-fixed PPP solution of Eqs. (3) and (4) are used as the input to build a high-precision correction model. The user side can utilize the high-precision atmospheric information to achieve high-precision positioning solutions based on Eqs. (3), (4), and (5). After the user obtains the accurate atmospheric information broadcast by the ground-based reference network, the troposphere and ionosphere corrections are added to the UDUC PPP normal equation as virtual observation equations (Psychas and Verhagen 2020; Psychas et al. 2018;Zha et al. 2021). The enhancement effect of atmospheric information on positioning can be adjusted according to the size of the virtual observation errors. The higher the matching degree between the virtual observation and its error, the better the positioning enhancement effect will be. Because the receiver Differential Code Bias (DCB) which is generally unknown for the user station is absorbed into the undifferenced ionosphere parameter, it is impossible to interpolate the unbiased and undifferenced ionosphere to constrain the PPP solution for a user. Hence, the single-difference-between-satellite ionosphere argumentation is used instead.

Regional troposphere modeling
To consider both the horizontal and vertical variations of the tropospheric ZWDs, a Modified OFC (MOFC) model is used. The MOFC model effectively matches the tropospheric characteristics in altitude by adding an exponential term of altitude dependence. The model is written as Cui et al. (2022): (6) T w (B, L, h) = a 0 + a 1 · dB + a 2 · dL + a 3 · dB · dL + a 4 · dB 2 + a 5 · dL 2 · e dH i H where a 0 , a 1 , a 2 , a 3 , a 4 and a 5 refer to the horizontalpart polynomial coefficients, dH i is the height difference with respect to the reference height, and H is the scale height. dB , dL are the longitude and latitude differences with respect to the reference coordinates. The number of the coefficients in this model is reduced compared to the traditional second-order OFC model. For the horizontalpart, a second-order polynomial model is chosen, while the altitude-part is replaced as an independent exponential component.
On the server side, the delays extracted from each reference station are used to calculate and broadcast the wide-area ZWD model coefficients. The model coefficients ( a 0 , a 1 , a 2 , a 3 , a 4 , a 5 and H) as well as the modeling precision (residuals RMS) are broadcasted to the users. The user can calculate the tropospheric ZWDs according to its locations. At the same time, modeling precision will be added as a constraint to improve the PPP solution.

Temporal and spatial modeling of the ionosphere interpolation error
In order to get the accurate ionospheric delay constraint on the user side, the slant ionospheric delays of three reference stations close to the user are used to interpolate the ionospheric information with the Inverse Distance Weighting (IDW) algorithm (Gao 1997) as follows: where I(u) is the ionospheric delay on the user side; I u, s i ref and w u, s i ref are the ionospheric delay and weight of reference station i around the user, respectively; d u, s i ref is the geometric distance from the user to the reference station i. Based on the above IDW algorithm, the ionospheric correction on user side can be obtained.
Unlike ZWD, the ionosphere is highly irregular, and its modeling accuracy is limited for a large area (Pi et al. 1997). For example, the accuracy of the current publicly available ionosphere global modeling products is only about 2.7 TECU (Total Electron Content Unit; Liu et al. 2021). In addition, the ionosphere is time-and spacedependent, therefore the interpolation error varies with different regions, different baseline lengths, and different time periods.
In order to calculate the RMS of the interpolated ionospheric delay, the precision of the ionospheric delay from each reference station and its attenuation or variation along with the baseline length should be known. Short  (Hu et al. 2003). Therefore, the key is to determine the relationship between the ionosphere correction error and the local station-space distance.
This study proposes to determine the ionosphere interpolating error by the cross-validation of the high-precision single-differenced-between-satellite ionospheric delay between the reference stations. The dependence of the ionospheric delay error on distance for each reference station is investigated. For reference station i, different reference networks from small to large scale formed by the nearby stations are selected to interpolate ionospheric delays. It looks like the different virtual reference stations derived from different reference networks. Two assumptions are made in our study: (1) the average distance between i and a reference network is the baseline distance between i and that virtual reference station; (2) the interpolated ionospheric delay is treated as the 'true' value (Psychas and Verhagen 2020). Under these two assumptions, the baseline length D i and the RMS of the interpolated error R i for each network can be calculated and input for building a mathematical function R = f (D) between the two variables for reference station i. With this function, the ionospheric delay error can be calculated using this reference station to interpolate the value. This cross-validation is performed at all reference stations to derive their RMS errors. Furthermore, the ionospheric delay error anywhere within the PPP-RTK service region can be calculated with the RMS function for the selected reference stations. As a result, the ionosphere delay precision map can be obtained, as shown in later section.
It is worth noting that the precision map reflects the situation of the current station configuration. If the stations are denser, the precision will be higher. One should note that the ionosphere activity rapidly changes with time, thus it is reasonable to calculate the error RMS function for each station and generate the precision map piece-wisely. In this study, the precision map is generated on a five-minute basis.

Algorithm and strategy summary
The flowchart of the overall system structure is shown in Fig. 1. The data processing on the server side Fig. 1a and user side Fig. 1b can be summarized as follows: For the server side, 1. Preparing the GNSS observation data, the precise orbit and clock products, DCB and ERP, and Antenna Exchange format (ANTEX) correction files; 2. Fixing the station coordinates and performing UDUC PPP-AR for each reference station to generate the single-station ZWD and ionosphere estimation; 3. The ZWD values of all stations and their coordinates are used as inputs to model the regional ZWD, and then generate the model correction coefficients, and obtain the residual RMS index; 4. Cross-validating the slant ionospheric delay for each reference station, and then fitting a linear error function with respect to the distance; 5. Optionally, generating the ionosphere error map by interpolating the error at each grid point with the linear error function using the data of the nearest three reference stations.
For the user side, 1. Filtering out the effective reference stations which can observe required target constellations and are within the required distance to the rover station. Constructing the reference network whose center is close to the rover. Choosing the network with the smallest average distance between the rover and reference stations; 2. Forming the GNSS observation functions with the measurements and the precise orbit and clock files, DCB and Earth Rotation Parameters (ERP) and ANTEX correction files; 3. Adding ZWD constraint according to the ZWD model; 4. Interpolating ionosphere value and its error, and then adding ionosphere constraints; 5. Executing the ambiguity-float PPP with robust estimation and then fixing the PPP ambiguity. Fig. 2a, to generate the high-precision atmosphere corrections. The EPN network is operated under the International Association of Geodesy (IAG) regional reference frame sub-commission for Europe to support the research on tectonic deformation, climate monitoring, weather prediction, Sea-level variations, etc.. But its coastline coverage is limited, so the SONEL network is involved to densify the overall server network in this paper. The stations with low availability of measurements or the solutions with a lower fixing rate will be filtered out. There are about 400 stations in the experiment. 40 stations are selected as the user stations, whose distribution is illustrated in Fig. 2b. These user stations are not included in the atmospheric modeling. The remaining server stations are densely distributed with an average spacing of 136 km. While in north Europe, the average spacing is about 295 km. The GNSS observation data with a sampling rate of 30 s in the above-mentioned network from the Day of Year (DOY) 213 to 243 in 2021 are used to perform the assessment. On the user side, the one-day observations are divided into 24 one-hour sections to get enough samples for generating reliable convergence curves, while on the server side, the daily observations are processed continuously to generate reliable atmosphere corrections. Excluding the missing observation files, the total number of hourly user samples is over 29,000. Precise orbit and integer-recovery clock products for the same days from GFZ, as well as the DCB products and ANTEX file from IGS, are applied in the processing.  In the processing procedure, some common strategies are adopted on both user and server sides: the elevation-angle-based weighting scheme is applied, cycle slips are detected by the Turbo-Edit method without repairing, the ionosphere delay and receiver clock bias are estimated as white noise, the ambiguity is estimated as a constant parameter if no cycle slip occurs, and the troposphere ZWD is estimated as random-walk process with a process noise of 1 cm 2 /h. To ensure the precision of the ambiguity parameters, the ambiguity with a lower elevation angle of 15 degrees or with a STD value larger than one cycle will be rejected in the ambiguity resolution. In AR processing, the WL ambiguity is first fixed by rounding, and then the NL ambiguity by LAMBDA (Least-square AMBiguity Decorrelation Adjustment), with the ratio value threshold of 3 for ambiguity validation (Teunissen and Verhagen 2004;Wang et al. 2022).

A combination of networks EPN (EUREF permanent GNSS network; Bruyninx et al. 2019) and SONEL (Système d'Observation du Niveau des Eaux Littorales) is used, as shown in
The position is estimated as a constant parameter on server side while as white noise on the user side. In addition, on the server side, the corrections of troposphere and ionosphere are delivered with an interval of 5 min. On the user side, three stations with a minimum average distance from user position are selected as reference stations in PPP-RTK processing. In order to investigate the performance of our PPP-RTK method for a wide-area network, we limit the reference station spacing to over 100 km. The inter-station distance of 40 groups of the network ranges from 110 to 350 km with the median of 171 km.

Results
The precision of the regional troposphere modelling and the interpolated ionosphere value are evaluated first. Then the GPS PPP-RTK positioning performance analysis is conducted. Finally, the contribution of multi-GNSS on PPP-RTK is investigated.

Evaluation of troposphere modeling precision
It is well known that the height positioning precision is much related to the tropospheric delay. In this part, the troposphere modelling precision is evaluated with the residuals of all service stations. With the coefficients from  Fig. 3. Figure 3 tells that more than 90 percent of the stations have the fitting residuals within 2.5 cm, and overall RMS is 2.6 cm. A few stations located in the boundary and large height difference areas have the residuals up to 8 cm. Overall, the results show that the regional ZWD model can achieve a homogeneous accuracy for the stations located in different altitudes and areas with the MOFC model.

Evaluation of ionosphere error map
The ionosphere error maps at 0 and 12 o'clock on DOY 215, 2021 are shown in Fig. 4. The color bar range is between 0 and 0.2 m, and no color is shown in the place where the precision is worse than 0.2 m. This is the reason why the coverage of the two precision maps is not the same. It is obvious that the ionosphere interpolation error at hour 0 is much higher than that at hour 12. In the Fig. 4(a), the ionosphere error is 20 cm or more in northern Europe. In such a case, it is difficult to augment the PPP-RTK with the imprecise interpolated ionospheric delay. At the same time, the ionosphere error at the lberian Peninsula is about 10 cm. However, in the Fig. 4(b), the errors in these two areas are much smaller, about 2-5 cm. The ionosphere error in Germany and its neighboring countries is stable and small. One reason is that there are lots of reference stations in this area and thus the station density is higher than that in other areas. For any user site in this area, it is easy to choose three nearby reference stations to provide the high-precision ionosphere correction. Figure 5 displays the time series of the difference between the interpolated ionospheric delays and the estimated ambiguity-fixed ionospheric delays at station MALA in Spain. The difference reaches up to 1 dm between 0 and 5 o'clock and is generally less than 5 cm between 6 and 18 o'clock, and slightly increases after hour 18. These results are consistent with the two error maps in Fig. 4. It indicates that a dynamic and adaptive ionosphere precision should be assigned for the PPP-RTK solution.

Positioning performance analysis
GPS PPP-RTK kinematic solutions are generated with four ionosphere weighting strategies, namely, the fixed RMS value of 1, 8 and 30 cm, and the automatically , and 30 cm serve as the extremely tight, reasonable, and conservative example, respectively, for the comparison purpose. They are abbreviated as "G-1 cm", "G-8 cm", "G-30 cm", and "G-auto". In these four solutions, the same troposphere constraints based on Eq. (5) are utilized. In addition, two other groups of solutions are also performed for comparison. They are the GPS PPP-AR without troposphere or ionosphere constraints, named as "G-AR", and the GPS PPP-AR with only ionosphere constraint using GIM product, named as "G-GIM". The "G-GIM" solution uses the final GIM product from Center for Orbit Determination in Europe (CODE) analysis center because the final CODE GIM provides quite reasonable RMS maps, and the distribution of the actual error is properly bounded by the normal distribution derived from the RMS map. The high ionospheric integrity can be ensured with the CODE GIM and is helpful to enhance the positioning performance of PPP . The time series of the 90th percentile of the horizontal and vertical errors for all test samples are shown in Fig. 6. The GPS PPP-AR solution takes the longest convergence time of over 60 min to reach 1 dm precision in two directions. It also has the largest positioning errors in the two directions, no matter how many epochs are processed. The "G-GIM" solution can systematically accelerate the convergence in two directions, especially in the horizontal direction and during the first half an hour. The "G-1 cm" PPP-RTK solution can significantly improve the performance compared with the "G-AR" and "G-GIM" solutions. The troposphere and ionosphere constraints can accelerate the ZWD and ionosphere parameter estimation on user side to reach a high-precision level and reduce their correlations with coordinate parameters. Even though the STD value of 1 cm for ionosphere error may be over-optimized as indicated in the ionosphere precision map shown in Fig. 4, its variance can somehow be adjusted by the iteratively extended Kalman filter. However, the effect of robust estimation would be impacted if there are many over-optimized observations. The conservative ionosphere constraint is always beneficial to UDUC PPP and the more conservative the constraint, the less significant the effect will be. The advantage of "G-30 cm" solution over the "G-AR", "G-GIM", and "G-1 cm" solutions indicates that the ionosphere corrections with a precision of 30 cm is still helpful for UDUC PPP. The "G-8 cm" PPP-RTK solution is better than the "G-1 cm" and "G-30 cm" solutions and can further speed up the convergence. The ionosphere precision map shows that the 8 cm ionosphere STD is reasonable or sometimes a little conservative for many user stations.
Among all six solutions, our "G-auto" PPP-RTK solution adaptively determining the ionosphere STD achieves the best performance. It slightly reduces the convergence time needed to reach the precision better than 1 dm. But after convergence, it can further improve the positioning precision obviously. As given in Table 2, compared with "G-8 cm" solution, it slightly reduces the convergence time from 4.5 and 23.5 to 4.0 and 20.5 min in horizontal and vertical directions, respectively. The average convergence time is also a commonly used indicator in many studies. We can see that the average convergence time of PPP-AR is 22.1 and 18.4 min while those of 'G-auto' PPP-RTK solution are 2.0 and 5.0 min, in horizontal and vertical directions, respectively.

GPS PPP-RTK versus multi-GNSS PPP-RTK
Multi-GNSS PPP positioning outperforms the singlesystem one in both convergence and precision. This has been well demonstrated for ambiguity-float and -fixed PPP (Hu et al. 2020;Kiliszek & Kroszczynski, 2020). We also use our PPP-RTK method to process the GPS + Galileo dual-system observations to quantitatively assess how much the performance can be improved. Aside from GPS and Galileo, most of the user stations can track GLObal NAvigation Satellite System (GLONASS) signals. The GLONASS observations are involved in the data processing on user side without ionosphere constraint or ambiguity resolution since the inter-frequency bias caused by FDMA (Frequency Division Multiple Access) signal structure impedes the GLONASS PPP ambiguity fixing. These two solutions are named as "GE-auto", "GE-auto + R".
The 90th percentile of the horizontal and vertical coordinate errors for all test samples are shown in Fig. 7. The convergence time of the 90th percentile for three solutions is given in Table 3, which indicates PPP-RTK, adding the Galileo system can significantly accelerate the convergence in the beginning period of 60 min. The "GE-auto" solution takes only 2.0 and 9.0 min to achieve the convergence in horizontal and vertical directions, respectively. After initialization of half an hour, the positioning error is stable in both directions which reaches 2.1 cm in the horizontal and 4.4 cm in the vertical direction. Compared with the "G-auto" solution, the "GE-auto" solution improves the convergence time by 50.0% in horizontal and 56.1% in vertical. The further integration with GLONASS only shows slight improvement when the observation length is around 20 min. In the beginning and after 30 min period, "GE-auto + R" nearly overlaps with the "GEauto" time series. It indicates that for GPS and Galileo combined PPP-RTK with atmosphere constraints, the effect of adding GLONASS in the ambiguity-float PPP on positioning accuracy is negligible. As for the average convergence time, the "GE-auto" solution takes only 1.4 and 3.0 min to converge, with an improvement of 30.0% and 40.0% compared with the "G-auto" solution in horizontal and vertical directions, respectively.
It is well known that integer ambiguity resolution is the key to fast and high-precision GNSS positioning and navigation. The Galileo observations are added together with the high-precision ionosphere constraints in this study. The Galileo PPP ambiguity resolution can be realized within a few minutes and contributes much to the user position estimation. While for GLONASS, the PPP ambiguity is kept float, and no high-precision ionosphere is derived on the server side. Therefore, the contribution of GLONASS is much smaller than Galileo.
In this contribution, another global satellite system BDS is not involved because of only 40% of the stations tracking BDS signals. In the future, we can expect a faster PPP-RTK solution with the multi-system of GPS, Galileo, BDS, and GLONASS.

Conclusions and remarks
For the ionosphere-weighted PPP-RTK model, both the a priori value of ionosphere and its precision are important to calculate the user position. The ionosphere is highly irregular and difficult to be modeled accurately. This study proposed a practical method to model the error of the interpolated slant ionosphere by cross-validation. A five-minute piecewise and distance-dependent function are built to represent the stochastic model of the slant ionosphere derived from each reference station. This function can be further used to generate the dynamic ionosphere precision map for the service area.
The EPN + SONEL networks are used in this study to provide the PPP-RTK service for the Europe mainland, and to assess the proposed ionosphere cross-validation method. The server side extracts the tropospheric and ionospheric delays from each station by the UDUC PPP AR solution. The regional MOFC troposphere model is built with a residual RMS of 2.6 cm to provide the troposphere augmentation for the users. Furthermore, Table 2 Average and 90th percentile convergence time for "AR", "G-GIM", "G-1 cm", "G-8 cm", "G-30 cm" and "G-auto" solutions  Table 3 Average and 90-percentile convergence time for "G-auto", "GE-auto", and "GE-auto + R" solutions the ionosphere correction precision map is generated to reflect the precision of the interpolated ionosphere value within the coverage of the PPP-RTK service. The RMS of interpolated slant ionospheric delay is generally less than 5 cm in most areas and periods, while it can reach 10 cm or more for some periods in the areas with sparse reference stations. Therefore, it is more reasonable to dynamically determine the ionosphere correction precision in wide-area PPP-RTK. The performance of our ionosphere precision map generated by cross-validation is also evaluated with the PPP-RTK solutions using 31 days of GNSS observations at 40 test stations. The hourly PPP-RTK was processed using three fixed ionosphere correction STD (i.e., 1 cm, 8 cm, and 30 cm) and the dynamic STD derived with our method, respectively. The comparison shows that our PPP-RTK solution outperforms all other solutions in terms of convergence time and positioning accuracy. The 90th percentile convergence time of our GPS kinematic solution is 4.0 and 20.5 min for horizontal and vertical directions, respectively. After integrating the Galileo into PPP-RTK, the 90th percentile convergence time of our GPS and Galileo combined kinematic solution is 2.0 and 9.0 min for horizontal and vertical directions, respectively. The average convergence time for horizontal and vertical directions, of our GPS kinematic PPP-RTK solution is 2.0 and 5.0 min and that of GPS + Galileo kinematic PPP-RTK solution is 1.4 and 3.0 min, respectively.
Currently, there is an increasing number of available GNSS frequencies for the newly launched GNSS satellites. Integrating the observations of multiple frequencies provides a stronger positioning model, potentially resulting in the faster PPP-RTK ambiguity resolution. The studies on investigating the PPP AR with multi-GNSS and multi-frequency observations have demonstrated that the horizontal convergence of the kinematic Galileo + GPS multi-frequency PPP AR needs 3.0 min on average and 5.0 min in the 90% of cases (Psychas et al. 2021), and the convergence can be achieved within 3 min on average if there are more than 10 Galileo + BeiDou satellites involved in PPP AR without any ionospheric corrections (Guo et al. 2021). These results imply that the rapid convergence of a few minutes is achievable by multi-GNSS and multi-frequency PPP AR even without the a priori ionospheric constraints. It is worth further investigating and comparing the contribution of multi-GNSS, multi-frequency, and precise atmosphere constraints on near-instantaneous cm-level PPP. Our future work will be also focused on utilizing the multi-frequency observations to improve the atmosphere estimation on the server side.
Availability of data and materials GNSS observations in the experiment were downloaded from the EUREF Permanent GNSS Network (ftp:// gnss. bev. gv. at/ pub/ obs) and SONEL network