A new large-area hierarchical PPP-RTK service strategy

Introducing precise atmosphere information into precise point positioning enables rapid ambiguity resolution and introduces a significant accuracy improvement. However, it can only be implemented in regions with dense networks and stable communication links. For larger areas, e.g., an intercontinental level, there is a conflict between the accuracy of corrections and the amount of atmosphere information to be disseminated. We develop a hierarchical augmentation mode to combine the advantages of the fitting model and region interpolation model to relieve the communication burden. Relying on the fitting model with fewer coefficients applied over large areas as the essential information, the unmodeled errors are calculated at each reference station, and further correction information is optional compensation depending on the magnitude of the unmodeled residuals. We perform the proposed models on 103 EUREF Permanent Network stations with 200-km station spacing and 84 stations as the external validation. The ionosphere and troposphere fitting models have an average accuracy of about 4.2 and 1.3 cm, respectively, under meteorologically calm conditions. The unmodeled error transmission determined by the magnitude of residuals can be reduced by 61% and 96% for the ionospheric and tropospheric delays, respectively, with respect to the legacy interpolation mode. Further compensation implemented, i.e., unmodeled residuals, can achieve instantaneous convergence for 83.6% of all solutions, and the overall initialization time is within 1.0 min. Thus, the proposed hierarchical positioning mode satisfies real-time positioning convergence requirements and significantly reduces massive corrections in communication.


Introduction
Although global navigation satellite systems (GNSSs) precise point positioning (PPP) technology has the advantage of providing high-precision station coordinates in a global reference frame without the need for local reference stations (Malys and Jensen 1990;Zumberge et al. 1997), the long convergence time for achieving a centimeter-level accuracy makes it vulnerable. The convergence time can be shortened by the ambiguity resolution (AR) technique for a single station (Ge et al. 2007;Li et al. 2010), whereas the rapid or instantaneous AR is achievable if PPP is augmented with precise atmosphere corrections (Wübbena et al. 2005).
Generally, the augmentation information can be generated from a reference network in the form of a fitting function model (de Oliveira et al. 2016;Boisits et al. 2020;Cui et al. 2022), a grid model (Li et al. 2021b), or the correction values of the reference stations (Gao et al. 2021;Li et al. 2022). Although these models can be applied for both large and regional areas, their accuracy depends on the number of model parameters and the area of its coverage for a certain reference network (Zha et al. 2021;Banville et al. 2022).
In order to achieve a rapid, even instantaneous PPP-AR, the high-precision correction usually relies on dense station networks in regional areas (Psychas et al. 2019;Li et al. 2021b). Currently, the atmosphere augmentation mode is prevalent in the multi-frequency and multi-GNSS (European GNSS Agency 2019). For example, Centimeter-Level Augmentation System (CLAS) for Quasi-Zenith Satellite System is implemented by dense stations via 30-km grid points (Cabinet Office 2021;Hirokawa et al. 2021). Moreover, the atmospheric delays from nearby stations also can be directly used to provide the corrections (Tao and Jan 2015). However, high representation accuracy requires high-resolution grids and high-order models for grid and fitting methods. It requires many parameters, especially for dense networks. Thus, the high-precision augmentation mode is currently only available for areas with very high communication capacity (Teunissen and Khodabandeh 2014;Psychas and Verhagen 2020;Banville et al. 2022). In addition, ensuring nearby stations or dense grid point information promptly and continuously to provide the required high-frequent corrections in larger areas is still challenging due to several reasons, e.g., communication interruptions, geological disasters, etc. The users need to change participating reference stations based on their real-time availability (Geng 2010;Li et al. 2010). Therefore, high demands are placed on real-time computation power for server and communication capability for larger areas.
In contrast, the large-area augmentation system, implemented by fitting models and broadcast via satellite with fewer communication volumes, provides more extensive services to achieve rapid positioning in most areas without a stable communication network. The ionospheric and tropospheric delays need to be considered and modeled separately. The tropospheric wet delay can be modeled well using a fitting model due to the significant altitudinal correlation in the troposphere (de Oliveira et al. 2016;Cui et al. 2022). However, the ionospheric delay modeling is more sophisticated than that of the troposphere because of its active variations (Wang et al. 2015(Wang et al. , 2020a. Currently, the large-area modeling usually adopts the vertical total electron content (VTEC) mode, while the imperfect mapping function could jeopardize the accuracy of high-precision slant TEC (STEC) (Boisits et al. 2020). The global ionosphere model (GIM) with a 5° × 2.5° sparse grid generated by the fitting function has a poor precision of about 2.7 TECU (Liu et al. 2021). As a result, precise modeling provided in large areas using sparse station networks is still challenging. Although the fitting model with lower dependence on communication could continuously and steadily provide augmentation information in large areas, ensuring the atmosphere information available anytime, their modeling precision is insufficient to enable instantaneous AR. In other words, only employing the fitting model is difficult to meet the demand for large-area coverage and high-accuracy service simultaneously (Li et al. 2019;Zhao et al. 2021).
Aiming to provide high-precision atmosphere corrections with low data communication volume and achieve high precision in large areas, we propose a hierarchical augmentation mode to combine fitting models and residual unmodeled compensation information from the region to achieve rapid/ instantaneous convergence in large areas with reduced correction volume. The combined mode takes the fitting model as the essential information in large areas, e.g., intercontinental, and the regional interpolation as optional corrections depending on the error magnitude and network communication capability. Additionally, the areas without regional corrections can also perform the augmentation with the fitting model by the satellite broadcast.
Thanks to the stable and high-precision single-station VTEC modeling method (Wang et al. 2020b;Li et al. 2021a), the satellite-plus-receiver (SPR) biases can be precisely and stably estimated, and thus, a more flexible and precise large-area satellite-wise STEC model can be generated. Given that the fitting model already corrects the majority of atmospheric errors, only the unmodeled residuals need to be compensated for achieving high-accuracy corrections. Therefore, further compensation for the large-area fitting model can be determined automatically according to the magnitude of unmodeled residuals to reduce the burden of data communication. We investigate the large-area atmosphere modeling performance generated by a sparse station-spacing network and analyze the data broadcast volume. In addition, considering the strong correction between ionospheric delay modeling and SPRs, we also evaluate the ionosphere modeling precision under mixed-receiver-type and fixed-receiver-type to analyze the impact of SPRs on ionospheric delay modeling. Moreover, the multi-GNSS positioning performance is also investigated in Europe using the proposed hierarchical augmentation mode.
In the methodology part, we first formulate the observation equations of the undifferenced and uncombined (UDUC)-PPP-AR method ) and then describe the estimation and modeling methods of the tropospheric and ionospheric delays, including the fitting and interpolation models. Then, the hierarchical augmentation mode is described in detail. In the experimental validation part, we briefly present the data processing strategy and introduce the data used in this study, which is used for both real-time uncalibrated phase delay (UPD), atmospheric product estimation, and PPP-AR validation. Then, the realtime orbit, clock, and UPD are analyzed, and the ionospheric delay modeling performance is investigated. Finally, we investigate the performance of multi-level positioning modes with multi-GNSS.

Methodology
We start with the GNSS observation equations to elaborate on the positioning and atmospheric delay calculation methods. Once the atmospheric delays are derived by PPP with ambiguity resolution, the ionospheric slant delay and tropospheric zenith wet delay (ZWD) models can be generated and applied for unmodeled error calculation on all reference stations. Then, the atmosphere correction method applied to the user is presented.

Raw observation equation
The observation equations of UDUC-PPP carrier phase L and pseudo-range P can be expressed as follows, where indices s , r , and f refer to the satellite, receiver, and frequency, respectively; s r is the range from satellite to receiver; dt s and dt r are satellite and receiver clock offsets, respectively; T s r is the tropospheric delay; I s r,1 is the ionospheric delay of the signal path at frequency 1 and f = f 2 1 ∕f 2 f ; f is the wavelength; N f is the integer ambiguity; b r,f and b s f are the receiver-and satellite-dependent UPDs, respectively; d r,f and d s f are the code biases of the receiver and satellite, respectively; P,f and L,f are measurement noise of the pseudo-range and carrier phase observations, respectively. It should be noted that the phase center offsets and variations (Rebischung and Schmid 2016), station displacements (Petit and Luzum 2010), phase wind-up (Wu et al. 1993), and relativistic delays (Ashby 2003) also should be corrected, although they are not included in the equations.
The slant tropospheric delay consists of the dry and wet components, both of which could be expressed by zenith delays with mapping functions (Boehm et al. 2006(Boehm et al. , 2014. The a priori zenith hydrostatic delay (ZHD) can be calculated using meteorological data and the Saastamoinen model (Saastamoinen 1972), and the ZWD is estimated as an unknown parameter.
For multi-GNSS observations, the additional inter-system bias (ISB) parameter between systems should be considered. Additionally, the inter-frequency bias (IFB) parameter should be added to all GLONASS satellites because the frequency division multiple access (FDMA) technology is used. Relying on the received real-time satellite orbit and clock as well as model correction items, the linearized (1) and (2) can be simplified as follows, where l s r,j and p s r,j denote observed-minus-computed (OMC) carrier phase and pseudo-range observations, respectively; u s r is the unit vector from receiver to satellite; Δr denotes the vector of the receiver position increment; m s,Sys r is the mapping function of ZWD Z r ; B s r,f = N s r,f + b r,f − b s f is the undifferenced ambiguity; Sys represents the satellite system, including GPS, GLONASS, Galileo, and BDS; ISB Sys represents the ISB parameters of Galileo and BDS relative to GPS; IFB s,R is the IFB parameter of each GLONASS satellite relative to GPS.
For each epoch, the estimated parameters are, The Kalman filter is employed for the parameter estimation. Z r and I s r,1 are estimated as random walk process noise. The receiver clock t r is estimated as epoch-wise white noise and the carrier phase ambiguities B s r , ISB Sys and IFB s parameters are estimated as constant over time.
Fixing the carrier phase ambiguity is essential to achieve stable and high-precision atmospheric delay from all reference stations and the key to realizing rapid positioning for the user (Ge et al. 2007;Cui et al. 2021). For users, the estimated UPD is used to implement integer AR, where wide-lane (WL) ambiguity is fixed by rounding to the nearest integer, and narrow-lane (NL) ambiguity is solved using the least-squares ambiguity decorrelation adjustment (LAMBDA) method (Teunissen 1995).

Atmospheric delay model
By applying the real-time UPD, orbit, and clock products, the undifferenced ambiguity can be fixed at all reference stations, and the slant ionospheric delay and tropospheric ZWD can be derived and used for modeling. (3)

Tropospheric delay modeling
For the tropospheric delay modeling, we use the modified optimal fitting coefficients (MOFC) model to model the estimated tropospheric ZWD : where the a 0 , a 1 , a 2 , a 3 , a 4 , a 5 are the fitting coefficients of the model, the e dh i H refers to the altitude-related correction, and H is the scale height, d i , d i , dh i are differences in latitude, longitude, and altitude between users and reference points.

Ionospheric delay modeling
The ionospheric delay estimation in (3) introduces the satellite and receiver differential code biases. These biases can be canceled in the interpolation method by the satellite common view. In the large-area modeling, these biases should be removed in advance to model the "clean" ionospheric delay. The estimated ionospheric value I s r,1 can be expressed as, where M(e) is the ionospheric delay mapping function; Î s r,1 is the "clean" slant ionospheric delay. The local ionospheric VTEC model on each station is introduced to model all satellite STEC and to separate the SPR biases from the ionospheric delays. The local ionosphere model is illustrated in (7) (Wang et al. 2020b), where d and d are the latitude and longitude differences between the ionospheric pierce point (IPP) and ground station; t loc is the local time; p max and q max are the maximum degree of the polynomials; k max is the maximum degree of the finite Fourier series; E pq , C k , and S k are the model coefficients to be estimated; H ion is the altitude of the ionospheric single-layer shell; and R E is the mean radius of the earth. In this study, the maximum degrees of the polynomials and the finite Fourier series are set as 2 and 4, respectively, and the height of the assumed single layer is set as 350 km. After estimating and eliminating the SPRs, a "clean" ionospheric slant delay can be used for modeling. We build the large-area satellite-wise STEC fitting model using all observed reference stations. A second-order polynomial model is described in (8).
are the fitting coefficients of the model; d i and d i are differences in latitude and longitude between the reference point and station, respectively.

Regional unmodeled correction
Although the fitting model can correct the majority of the tropospheric and ionospheric delays, the unmodeled part remains on each reference station. Hence, the unmodeled residuals at these stations can further be disseminated to users.
where Ĩ s r,1 and Z r are ionospheric and tropospheric modeled values using the received coefficients, respectively; ls r and lt r are ionospheric and tropospheric unmodeled residuals, which can be further interpolated to users. Generally, the inverse distance weighting (IDW) algorithm is used to interpolate the unmodeled values from the selected reference stations (usually using three nearby stations) to users.
where ̂l pos is the interpolated atmospheric residual part on the user side; l s i (ref) and w s i are the atmospheric residual delay and weight of reference station i , respectively; d s i is the geometric distance from the user to the reference station i. Figure 1 presents the structure of the hierarchical augmentation mode, including UPD estimation, SPRs estimation, atmosphere modeling, unmodeled correction generation, and positioning. After receiving the real-time orbits and clocks, the static float PPP is first performed on all global (8)

Hierarchical positioning system
Once the WL and NL UPDs are successfully estimated using global reference stations, the fixed solution can be implemented at all reference stations. Then, the atmospheric delays can be precisely derived and used for modeling. We use a sparse network in the large area to build tropospheric ZWD and satellite-wise ionospheric STEC delay models. Before the ionosphere modeling, the SPRs should be estimated and removed from the PPP-derived ionospheric delays station by station. We estimate the SPRs using daily slant ionospheric delays of all satellites on each station. Additionally, we adopt a seven-day sliding window to smooth the SRP estimates, which will further improve the stability of SPRs. The fitting model coefficients can be broadcast via satellite communication links for the all-time all-region service to provide augmentation information for users. Given that the fitting performance in some areas is affected by topography variations, meteorology changes, etc., part of unmodeled errors still exists at reference stations. After subtracting the ionospheric and tropospheric modeled values, the unmodeled values are complementarily provided from the nearby stations to the user by the IDW interpolation model. It should be mentioned that the AR at the reference station is very critical, as only corrections with fixed ambiguity can be used (Li et al. 2010).
As for the positioning, with different augmentation information listed in Table 1, a three-level augmentation mode is performed. In the first level, i.e., PPP and PPP-AR, orbit, clock, and UPD products are provided for global users as essential corrections. The second level, i.e., PPP-AR with wide-area augmentation (PPP-WA), includes first-level and troposphere ZWD and satellite-wise slant ionospheric delay fitting models. The third level, i.e., PPP-AR with both widearea and region augmentation (PPP-WRA), is based on the first-and second-level and combines the unmodeled corrections from three nearby reference stations.
At the server, the atmosphere parameters are estimated using the random walk process. In contrast, at the user, atmospheric delays from the fitting model or interpolation model can be used as virtual observations with an a priori constraint. Since the fitting model aims to achieve the overall optimal solution in large areas, the model precision hardly reflects actual precision in all areas. Therefore, we set three times the modeling precision as the constraint value for the fitting model, while the constraints can be directly determined by inter-satellite cross-verification after implementing unmodeled correction ). The user should impose a constraint on the corresponding parameters with the following equations, where x est and x ref are estimated and external correction values, respectively; 2 ref is the a priori variance factor from the fitting model or error function, which serves as the constraint.
Compared to the legacy interpolation mode, the hierarchical augmentation mode can significantly reduce the burden of data communication volume. Because the fitting model has already corrected the majority of tropospheric and ionospheric delays in large areas, not all unmodeled errors of all epochs need to be broadcast. Therefore, the unmodeled errors lower than 2.5 cm, which is lower than a quarter of one cycle NL wavelength, can be considered insignificant for AR, this is, only residual unmodeled errors larger than this threshold need to be provided to users (Psychas et al. 2019). Combining the region interpolation optional compensation and fitting model is crucial to reduce the communication burden, accelerate convergence, and achieve better positioning performance in large areas.

Experimental validations
This section presents the datasets and data processing strategy at the server and user of the experiments. The global and large-area reference stations applied for UPD estimation and atmospheric delay derivation are presented. In addition, more details at the server, including UPD estimation, atmosphere modeling, unmodeled errors calculation, and user positioning methods, are described.

Data set
Three months (DOY 120-210, 2022) of multi-GNSS dualfrequency data from 83 globally distributed MGEX permanent stations (Fig. 2) and 187 EUREF Permanent GNSS Network (EPN) stations (Fig. 3) are processed. The global stations are used for UPD estimation. Moreover, 103 EPN stations (all with GPS and Galileo observations) are used as servers performing static UDUC-PPP-AR to calculate the atmospheric delays, and the rest 84 stations (all with multi-GNSS observations) serve as the positioning validation stations.
As shown in Fig. 3, the 103 service stations (blue dots) with 200 km station spacing enable uniform coverage of the Europe region. Moreover, to evaluate the performance of the proposed hierarchical method in the large area, we limit the distance of PPP-WRA between reference stations and users larger than 100 km.

Data processing
We use the Positioning And Navigation Data Analyst (PANDA) data processing platform (Liu and Ge 2003) for UPD estimation. The orbit and clock products are from GFZ real-time product stream, which are also processed by PANDA. Orbit products are updated every three-hour using the latest 24-h observations, and the clock products are updated every five-second. The data processing implements multi-GNSS, while the ionospheric delay augmentation and PPP-AR are only performed with GPS and Galileo systems because only part of the stations can observe BDS satellites. In addition, GLONASS FDMA mode will introduce more complex AR solutions. WL and NL UPDs, fitting models, and unmodeled correction interpolation are updated every 30 s. Since ionospheric variations are strongly influenced by latitude, we divide the large area (Europe region) into four sub-regions according to the latitude 50°N and longitude 15°E lines for ionospheric delay modeling to achieve the consistent accuracy across regions. The average station coordinates from PPP daily static solutions are used as the reference for positioning verification. The details of server and user data processing settings are listed in Table 2.

Real-time products precision
This section evaluates the precision of real-time orbit, clock, and UPD products. As we mainly analyze the positioning performance in the European region, only the results of the orbit and clock of the medium Earth orbit (MEO) satellites are presented. Figure 4 shows the signal-in-space ranging error (SISRE) value of the real-time products with respect to the GBM products for all satellites during 90 days (Montenbruck et al. 2014). Due to the datum differences in clock precision calculation, we set the average satellite clock offsets as the reference. The SISRE values are 22.7, 29.2, 31.0, and 37.1 mm for GPS, GLONASS, Galileo, and BDS satellites, respectively. Among them, the GPS satellites have the best precision, followed by the Galileo and GLONASS satellites. Figure 5 provides the WL and NL UPDs time series for all observed GPS and Galileo satellites on DOY 150, 2022. Ionosphere-free combination for UPD estimation UDUC method for the atmospheric delay, unmodeled correction calculation, and positioning Ionospheric delay Server-side Estimated with random walk process to extract the values Positioning The additional augmentation information is used to provide the a priori ionospheric delay and the corresponding constraints Tropospheric delay Server-side A priori ZHD from Saastamoinen equation with GPT2w meteorological data input, ZWD estimated with random walk process, GPT2w&VMF1 for the mapping function Positioning The priori delay and mapping function are the same as at the server, but the additional augmentation information is used to provide the a priori ZWD value and the corresponding constraints   In this section, the precision of large-area troposphere and ionosphere modeling, as well as unmodeled correction, are evaluated. First, the tropospheric ZWD and ionospheric slant delay fitting model are assessed by the differences between PPP-derived and model-fitted values. Then, once the reference stations receive the fitting model coefficients, the unmodeled errors are obtained, and the volume of data transmission is analyzed.

Troposphere model precision
For the tropospheric delays modeling precision in the large area, we calculate the differences between modeled and estimated ZWD values on all stations over 90-day and present the RMS of the differences in Fig. 6. On all server and user stations, fitting residuals are less than 2.8 cm for more than 90 percent of the stations, and the average precision is 1.3 cm. A few stations located in the boundary of the modeling region or with large height differences have residuals up to 4.6 cm. Overall, the results show that the ZWD model achieves a nearly homogeneous accuracy for the stations located in different altitudes and areas.

Ionosphere model precision
Unlike the tropospheric delay modeling, the ionospheric delay is modeled satellite-wise, i.e., one set of parameters for each satellite. Before ionospheric delay modeling, the SPRs should be separated from the PPP-derived ionospheric delay on each station. We first give an example of the estimated SPRs values of GPS and Galileo at station 0LOV in Fig. 7.
We can see that the SPRs are very stable over the 30 days for both GPS and Galileo satellites; fluctuations are usually within 5 cm. Different types of receivers have different biases, which could introduce the biases and jeopardize the accuracy of ionospheric delay modeling if the SPRs are not well considered. Considering the similar performance of the same type of receiver, we also compare the ionosphere modeling performance of five receiver types. As an example, the slant ionospheric delay differences between modeled and estimated values of fixed-receiver-type with (Fig. 8a) and without removing SPRs (Fig. 8b) at each IPP on DOY 150 14:00-15:00 (each point denotes the average value of ten epoch residuals), 2022 are given.
From Fig. 8, the average modeling precision of fixedreceiver-type with and without removing SPRs are 3.7 and 9.3 cm, respectively. From the top panel, after removing SPR values, the "clean" ionospheric delays are used to perform the satellite-wise STEC modeling. More than 95% of  the IPPs have the model residuals of less than 4.5 cm, and only a few located in the boundary areas or at the beginning of the observation periods have large residuals of up to 17.4 cm. However, the bottom panel shows relatively larger biases than the scheme of considering SPR correction. It demonstrates that even if the same type of receiver is used, different hardware versions and/or other impacts, e.g., local temperature, also could impact the SPRs (Cui et al. 2021).
We also calculate the RMS of daily ionospheric delay differences of mixed-receiver-type with SPRs removed and fixed-receiver-type with and without SPRs removed for all stations over 90 days and present the statistical results by boxplot in Fig. 9. We further compare the selected five types and show them in Table 3. The average modeling precision of GPS and Galileo satellites are 4.6 and 3.9 cm, respectively, when considering the SPRs correction using mixed-type. In contrast, the fixed-receiver-type without removing SPRs presents a slightly larger bias than the scheme of mixed-receiver-type by removing SPRs.
Among them, the LEICA GR50 type has the best precision, followed by the LEICA GR30, LEICA GR25, and the Trimble NetR9 receiver types. However, with removing the SPRs, the fixed-receiver-type schemes show the best performance, indicating that removing the SPRs can significantly mitigate the impact of receiver-related biases. In addition, some receiver-type number is lower than others, but they are concentrated in a smaller area. Overall, the Galileo satellites have better precision and consistent modeling performance among different satellites than that of GPS satellites, which could be caused

Fig. 10
Selected stations for the demonstration of different positioning modes (red star) and the corresponding server stations used for regional unmodeled error calculation (blue dot). The distances between user stations and reference stations are given

Residual unmodeled error analysis
For a comprehensive evaluation of the residual unmodeled values, we first present the results of 0LOV, BUTE, EUSK, and ZIM3 stations as an example, as shown in Fig. 10. The distances with respect to corresponding server stations vary between 70 and 350 km, and the average distances are After receiving the fitting model coefficients, the obtained modeled values are compared with those estimated at user stations derived in the same way as for the reference stations to evaluate the accuracy of the fitting models. We present the differences between modeled and PPP-derived values of tropospheric ZWD in Fig. 11 and that of ionospheric delay in Fig. 12.
From Figs. 11 and 12, the unmodeled values of the tropospheric and ionospheric delay models in these four stations are presented. The RMS of ZWD differences between modeled and PPP-derived are 0.82, 0.43, 0.73, and 0.89 cm for stations 0LOV, BUTE, ZIM3, and EUSK, respectively. On the other hand, the variation of ionospheric delays is more significant than that of tropospheric delays. Therefore, the model precision of ionospheric delay is lower than that of tropospheric delay. The RMS ionospheric delay differences between modeled and derived values for these four stations are 3.48, 2.73, 3.97, and 2.01 cm for Galileo and 3.86, 2.84, 4.39, and 2.77 cm for GPS. The performance of ionospheric delay modeling for Galileo satellites is better than that of GPS satellites, which can be attributed to the more stable SPRs estimation and can be related to hardware performance of the Galileo satellites.

Unmodeled error volume analysis
From the fitting model results, most atmospheric delays are corrected. Users need only the residual unmodeled tropospheric and ionospheric delays for precise region correction. We further count the unmodeled values for all epochs of all stations and present the distribution in Fig. 13.
From Fig. 13, the majority of unmodeled values are within 3.0 and 20.0 cm for the tropospheric and ionospheric delays. When the unmodeled values are less than the set threshold (2.5 cm), there is no need to provide further compensation. This value is mainly set for the unmodeled ionospheric delay errors to ensure that the possibility and correctness of ambiguity fixing would not be significantly impacted. For tropospheric delay, most of the unmodeled residuals (about 82%) are within 1.5 cm and about 15% of residuals are between 1.5 and 2.5 cm. For these larger residuals, most of the stations are located in the higher altitude areas, with larger differences with respect to the mean level. In addition, iterative data processing is also introduced to detect outliers and continuously de-weight them to reduce their impact on the solution. Therefore, the impact of the same threshold applied for the tropospheric and ionospheric delay on the positioning could be mitigated. Compared with the legacy interpolation mode, corrections of tropospheric and ionospheric can achieve 96.8% and 61.1% bandwidth savings. Although the large-area fitting model can achieve the cm-level precision, it still has around 5.1% remaining ionospheric unmodeled errors greater than 10 cm. Taking the fitting model as an essential correction for large areas and the unmodeled errors as compensation, the proposed hierarchical model achieves a lower data  communication volume than the OMC-based and atmosphere-based interpolation-only modes. In order to compare the data transmission volume of different augmentation modes, the data communication volume of different modes is presented in Table 4. It should be noted that the percentile number applied in the hierarchical mode is taken from Fig. 13. Among them, the fitting mode includes only tropospheric and ionospheric delay fitting model coefficients, i.e., 7 for troposphere and 6 × S for ionosphere, which takes the minimum data in augmentation. At the same time, it can only achieve a faster convergence due to the poor correction precision in the large-area service. In contrast, higher accuracy interpolation usually requires massive correction when applied in dense networks over large areas, e.g., N and S × N for troposphere and ionosphere in the atmosphere mode and N × S × F in OMC mode.

Positioning performance verification
For a comprehensive evaluation of the different levels of positioning performance, Fig. 14 presents the positioning results of PPP, PPP-AR, PPP-WA, and PPP-WRA modes using four stations (shown in Fig. 10) as an example. We divide daily observation into 24 sub-sessions, i.e., one hour per sub-session, resulting in a total of 181,440 sub-sessions in 90 days. More details about data processing strategies are presented in Table 2. The float PPP scheme has the slowest convergence of about 13 min and the worst positioning precision of about 1.8 cm and 3.5 cm for horizontal and vertical components, respectively. The PPP-AR solution, which enables AR technology, can accelerate the convergence to 8 min and improve the precision. The PPP-WA can further reduce the convergence time and time to first fix (TTFF) to around 3.0 min, and the positioning precision is improved to 0.8 cm and 1.3 cm in the horizontal and vertical components, respectively. TTFF indicates the time when the ambiguity is first successfully and continuously fixed (ten epochs). Convergence is defined as the time needed to reach a positioning error of less than ten centimeters and last for at least ten epochs. When the unmodeled information is available from nearby reference stations, in the case of PPP-WRA mode, the ambiguity resolution can be achieved in the first epoch at 0LOV and EUSK stations. For ZIM3 and BUTE, ambiguity resolution is achieved at the second epoch, mainly due Moreover, we inspect the number of satellites that are tracked and those with ambiguity fixed in different positioning modes and present the results in Fig. 15. In general, the ambiguity of half of the satellites can be fixed, corresponding to the number of GPS and Galileo satellites, as we do not perform AR on the GLONASS and BDS satellites. The major impact caused by using large areas or region augmentation information is during the convergence period, where more satellites are fixed with this aid.
In addition, from Fig. 12, unmodeled values of some epochs are significantly larger than other periods, especially between 12 and 13 o'clock on station 0LOV. Therefore, we also present the corresponding positioning results, as shown in Fig. 16. It takes about 10 min to achieve successful convergences for horizontal and vertical components when only relying on the PPP-WA. In contrast, this convergence time drops drastically at the first epoch after PPP-WRA is achieved, and the vertical component shows the most significant improvement with respect to other modes. This result verifies that the precision of the hierarchical augmentation mode, though deteriorated by considerable ionospheric delays, still has improved convergence and precision.
Finally, we summarize the result of all stations and present the 68th (good) and 90th (poor) percentile positioning results, convergence time, and the TTFF in Fig. 17 and Table 5. Positioning precision of the 68th percentile improves by 14% and 21% for PPP-WA with respect to PPP-AR and PPP-WRA with respect to PPP-WA, and that of the 90th percentile is 6% and 7%, respectively. Moreover, convergence times of the 68th percentile improve by 80% and 68% for PPP-WA and PPP-WRA, and that of the 90th percentile is 51% and 61%, respectively. The 68th percentile TTFF of PPP-WA and PPP-WRA are 2.0 min and 0.5 min, and that of the 90th percentile are 8.0 min and 4.0 min, respectively. Compared with the PPP-AR and PPP-WA, the performance of PPP-WRA is greatly improved. The instantaneous AR can be achieved in the station with average reference station distance of less than 180 km. Moreover, for all solutions, 84% of solutions can achieve instantaneous AR. It should be noted that the range of nearby stations applied in the PPP-WRA interpolation distances is 23-187 km, and the average distance is about 113 km. Due to the strong correlation between the atmospheric delay variation and the meteorological conditions, the modeling and interpolation corrections and the unmodeled volume are also affected by the different meteorological conditions. In this experimental verification analysis, we generate the ionosphere and troposphere models in a relatively calm variation without extreme meteorological conditions. It should be noted that a larger volume of data could appear under the extreme meteorological condition.

Conclusions and remarks
We present a large-area hierarchical augmentation positioning strategy combining the fitting model and supplementary interpolation mode and analyze their performance. The large-area sparse reference stations generate the satellitewise ionospheric slant delay and tropospheric ZWD fitting models as the essential augmentation. The unmodeled errors can be optionally provided according to the magnitude of the residuals, which significantly relieves the communication pressure. This study can be taken as a demonstration of the real-time PPP services efficiency in the perspective of future complete global applications.
In total, 90 days of data in 2022 from 83 stations are used for UPD estimation and 187 stations were selected to investigate atmospheric delay modeling and positioning performance. We found that WL and NL UPDs are quite stable over time with STD of less than 0.02 and 0.03 cycles, respectively. The troposphere fitting model achieves an average precision of 1.3 cm. Moreover, SPRs daily biases are quite stable over time, with STD of less than 0.04 m and 0.03 m for GPS and Galileo satellites, respectively. This favorable temporal stable property facilitates their precise predictions for real-time ionospheric delay modeling. The large-area ionospheric delay models by removing the SPR biases, have an average precision of about 4.2 cm. In contrast, the differences across receivers could still jeopardize the model accuracy even if using a consistent type of receiver for modeling.
Overall, the PPP-WA positioning accuracy at the 68th percentile after receiving the modeling coefficients can reach 0.8 cm and 1.5 cm for the horizontal and vertical components, respectively, while those of PPP-AR solutions are 1.0 cm and 1.7 cm. As a result, 7.5, 6.0, 1.0, and 0.5 min and 10.5, 9.5, 2.0, and 0.5 min for PPP, PPP-AR, PPP-WA, and PPP-WRA are correspondingly required to achieve horizontal and vertical positioning errors less than 10 cm. Thanks to the improved precision of atmospheric delays obtained from the optional residual compensation, PPP-WRA can achieve instantaneous AR for 84% of all solutions, and the overall mean initialization time is within 1.0 min. In contrast, PPP-WA has already reached rapid  Finally, from the previous analysis, it is clearly shown that the fusion of multi-level augmentation can significantly improve the accuracy, continuity, and reliability of positioning with less communication burden. This is crucial since such capability will dramatically increase the GNSS applicability in various geoscience and commercial applications.