Introduction

Precise point positioning (PPP) (Malys and Jensen 1990; Zumberge et al. 1997) is a commonly used technique for real-time Global Navigation Satellite Systems (GNSS) applications that employ precise orbits and clocks derived from the global network to achieve high-precision positioning (Fotopoulos and Cannon 2001). Currently, thanks to the availability of uncalibrated phase delay (UPD) products (Ge et al. 2008), the real-time PPP ambiguity resolution (AR) with a convergence time of about 20–30 min can be achieved all over the world (Dow et al. 2009; Cui et al. 2021). However, the convergence time is still a limiting factor to further promote the real-time PPP applications.

One of the main errors limiting the PPP convergence time is the atmospheric delay error (Teunissen and Khodabandeh 2015). The satellite signals travel through the atmosphere and suffer considerable delays in the ionosphere and the troposphere. The ionospheric delay is usually mitigated by the ionosphere-free (IF) combination for dual-frequency users or estimated as an epoch-wise parameter for the uncombined solution (Zhao et al. 2019). The zenith tropospheric delay (ZTD) is usually divided into zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD) (Davis et al. 1985). ZHD can be precisely modeled given the surface pressure, temperature, and humidity observations, whereas ZWD is related to the water vapor content in the atmosphere, which can hardly be precisely modeled because of its rapid variations in time and space (Bevis et al. 1992; Hadas et al. 2020). Due to the strong correlation between ZTD, station coordinates, and receiver clocks, precise modeling and estimation of tropospheric delay in PPP is critical to both the convergence time and positioning precision and is still a big challenge (Yao et al. 2016; Wang and Liu 2019).

To decorrelate the ZTD and station coordinate parameters during the GNSS processing, precise external tropospheric delay products are used as a priori constraint, e.g., the tropospheric delay from numerical weather model (NWM) (Hobiger et al. 2008; Andrei and Chen 2009; Lu et al. 2017; Wilgan et al. 2017), empirical tropospheric delay models such as global pressure and temperature 2 (GPT2w) and GZTD2 (Böhm et al. 2015; Yao et al. 2016; Chen et al. 2020), in situ instrument measurements such as water vapor radiometer (Ware et al. 1993; Alber et al. 1997), and Raman lidar (Bock et al. 2001; Bosser et al. 2009). On the other hand, using the tropospheric delay estimates of a regional GNSS reference network, a correction model can be generated and broadcasted to the users within the region, that is, the tropospheric delay regional augmentation (Takeichi et al. 2009; Fund et al. 2011; Kalinnikov et al. 2012). The regional tropospheric delay model is usually determined by either a regional grid (Rózsa et al. 2020; Li et al. 2021) or the polynomial function (Zou et al. 2018; Zhou et al. 2020). Depending on the service region and the resolution, the distances between grid points can be large and cause communicating burden, whereas the polynomial function is usually determined with up to ten coefficients (second-order) and can be easily broadcasted to users with one-way communication, which is more feasible for small-area.

Thanks to the high precision tropospheric delay (around 6.5–13 mm) retrieved from real-time PPP solutions (Hadas et al. 2013, 2020), the regional tropospheric delay augmentation can provide precise products with a rapid update rate of 5-min, and can be widely used given the well-distributed reference stations. The contribution of regional tropospheric delay augmentation to the float GPS and BDS PPP solutions was investigated in the mainland China regions (Zheng et al. 2018; Zhang et al. 2018), where the convergence speed is improved by 29%, and the positioning accuracy is also improved by 18%. Other studies investigate enhancing the PPP using regional tropospheric delay augmentation with improvements up to 1% (east), 20% (north), and 5% (up) (Oliveira et al. 2017), the modeling precision under different weather conditions (Shi et al. 2014), the instantaneous GPS + BDS positioning enhanced with dynamic atmospheric constraints (Zou et al. 2018), and the impact of troposphere augmentation on single-frequency PPP (Zhou et al. 2020).

One important issue in tropospheric delay regional modeling is that the tropospheric delay vertical decrease has to be precisely modeled due to the correlation between atmospheric pressure and altitude. The vertical decrease of pressure can be approximated by linear or quadratic functions, for instance, in the optimal fitting coefficients (OFC) model, which is based on the polynomial functions (Shi et al. 2014), or modeled using the exponential function which corresponds to the exponential decrease of water vapor (Troller 2004; Dousa and Elias 2014; Yu et al. 2017). However, the linear or quadratic approximation can cause large modeling errors for a region with large altitude variations, for instance, up to 39 mm in the France region (Oliveira et al. 2017). Note that in previous studies, the performance of the OFC model combined with the exponential vertical decrease has not been investigated yet.

Despite the previous studies which mainly focus on some specific areas using the OFC model, e.g., the French (Oliveira et al. 2017), China Hubei (Shi et al. 2014), or the large region using the grid model, e.g., in mainland China (Zheng et al. 2018; Zhang et al. 2018) and Australia (Li et al. 2021), the performance of the polynomial model in wide-area has not been fully investigated. Zus et al. (2019) and Selbesoglu (2019) applied their interpolation model in the German and Turkey region using a relatively dense network of stations, but the corresponding impact on the positioning was not given. Moreover, the previous studies mainly focus on the GPS-only PPP or RTK service and the GPS + BDS solutions, whereas the performance of Galileo and GPS + Galileo PPP-AR with regional tropospheric augmentation has not been reported yet. The Galileo system with 24 satellites already available (as of July 15, 2021) can already provide full positioning capacity. Glaner and Weber (2021) evaluated the GPS + Galileo PPP-AR performance using Multi-GNSS Experiment (MGEX) stations with products from different International GNSS Service (IGS) real-time analysis centers. On average, the horizontal coordinates can converge to the centimeter level within seven minutes. However, the performance of regional troposphere augmentation in Galileo-only PPP-AR solution has not been investigated. Moreover, the previous studies cover Mainland China (Zheng et al. 2018; Zhang et al. 2018; Li et al. 2020), Australia (Li et al. 2021), Japan (Takeichi et al. 2009), whereas other regions, such as Europe (up to 2.9 km) is still under investigation.

The purpose of this study is to investigate the wide-area tropospheric delay augmentation over the Europe region in different station-spacing networks with a focus on the vertical decrease modeling and the impact of troposphere augmentation in the GPS and Galileo PPP-AR, including the GPS-only, Galileo-only, and GPS + Galileo solutions. We deliberate the altitude-related variations to combine polynomial with exponential functions to improve the troposphere model accuracy in wide-area using sparse station-spacing. Moreover, the performance of PPP-AR augmented with the proposed model in the Europe region in terms of positioning accuracy and convergence is evaluated and analyzed for GPS, Galileo, and GPS + Galileo constellations.

In the method part, we first formulate the observation equations for the undifferenced and uncombined (UDUC) PPP model and then describe the tropospheric delay modeling, including the OFC and MOFC models. In the experimental validations part, we briefly present the data processing strategy, and introduce the EUREF Permanent GNSS Network (EPN) and EPN densification network observations, which are used for both tropospheric delay modeling and PPP-AR validation. In the performance of tropospheric delay augmented PPP-AR part, we investigate the performance of both OFC and MOFC regional tropospheric augmentation in GPS and Galileo PPP-AR using about 200 stations.

Method

We start with the UDUC PPP model to elaborate on the methods of extracting the precise tropospheric delays, which are used in the regional modeling at the server, and then present the OFC and MOFC models.

GNSS observation equations

The linearized equations for PPP with raw observations can be expressed as,

$$\left\{ \begin{gathered} P_{{r,f}}^{s} = \vec{e}_{r}^{s} \cdot \vec{x} + c\left( {dt_{r} - dt^{s} } \right) + T_{r}^{s} + \gamma _{f} I_{{r,1}}^{s} + c\left( {d_{{r,f}} + d_{f}^{s} } \right) + \varepsilon _{{P,f}} \hfill \\ L_{{r,f}}^{s} = \vec{e}_{r}^{s} \cdot \vec{x} + c\left( {dt_{r} - dt^{s} } \right) + T_{r}^{s} - \gamma _{f} I_{{r,1}}^{s} + \lambda _{f} N_{f} + c\left( {b_{{r,f}} + b_{f}^{s} } \right) + \varepsilon _{{L,f}} , \hfill \\ \end{gathered} \right.$$
(1)

where superscript \(s\), subscript \(r\), and \(f\) denote a specific satellite, receiver, and frequency band, respectively; \(P_{r,f}^{s}\) and \(L_{r,f}^{s}\) are code and carrier phase observation values; \(\vec{e}_{r}^{s}\) is the unit vector from satellite to receiver; \(\vec{x} = \left[ {x y z} \right]\) is the receiver coordinates vector; \(c\) is the light speed; \(dt_{r}\) and \(dt^{s}\) are receiver and satellite clock offset; \(T_{r}^{s}\) is the slant tropospheric delay; \(I_{r,1}^{s}\) is ionospheric delay along with the signal from the satellite to a receiver at the first frequency and \(\gamma_{f} = f_{1}^{2} /f_{f}^{2}\), i.e., the factor of frequency \(L_{1}\) and other frequency; \(\lambda_{f}\) is wavelength; \(N_{f}\) is phase ambiguity; \(d_{r,f}\) and \(d_{f}^{s}\) are the frequency-dependent code biases for the receiver and satellite, and \(b_{r,f}\) and \(b_{f}^{s}\) are the biases in phase; \(\varepsilon_{P,f}\) and \(\varepsilon_{L,f}\) include the multipath effect and measurement noise for the code and phase, respectively.

To fix the integer ambiguity for PPP, the UPD has to be removed from the carrier phase ambiguities. Although in UDUC PPP, we usually estimate the ambiguities on L1 and L2 frequencies, the estimates are not stable due to their high correlation with the ionospheric parameters. Therefore, they are converted to IF ambiguities for ambiguity resolution and UPD estimation as in the IF-PPP (Gu et al. 2015). The wide-lane (WL) UPDs are calculated from the Melbourne–Wübbena (MW) combination, and the narrow-lane (NL) UPDs are estimated based on the NL ambiguities derived from IF ambiguities with fixed WL integers, and their usability in UDUC-PPP has been demonstrated by Zhang et al. (2019) and Du et al. (2020). Generally, \(\overline{N}_{r,IF}^{s}\) is decomposed into the following combination of integer WL \(N_{r,WL}^{s}\) and float NL \(\overline{N}_{r,NL}^{s}\) ambiguities for ambiguity fixing (Cui et al. 2021). The relationship among the IF, L1, L2, WL, and NL ambiguities can be expressed as follows to understand the principle for UPD estimation at the server and using UPDs for ambiguity fixing at both the server and user.

$$\overline{N}_{r,IF}^{s} = \frac{{f_{1}^{2} }}{{f_{1}^{2} - f_{2}^{2} }}\overline{N}_{r,1}^{s} - \frac{{f_{1} f_{2} }}{{f_{1}^{2} - f_{2}^{2} }}\overline{N}_{r,2}^{s} = \frac{{f_{1} }}{{f_{1} + f_{2} }}\overline{N}_{r,NL}^{s} - \frac{{f_{1} f_{2} }}{{f_{1}^{2} - f_{2}^{2} }}N_{r,WL}^{s},$$
(2)

with

$$\begin{gathered} \overline{N}_{r,WL}^{s} = N_{r,WL}^{s} + d_{r,WL} - d_{WL}^{s} \hfill \\ \overline{N}_{r,NL}^{s} = N_{r,NL}^{s} + d_{r,NL} - d_{NL}^{s} , \hfill \\ \end{gathered}$$
(3)

where \(N_{r,WL}^{s}\), \(\overline{N}_{r,WL}^{s}\), \(N_{r,NL}^{s}\), and \(\overline{N}_{r,NL}^{s}\) are WL integer and float, and NL integer and float ambiguities, respectively; \(d_{r,WL}\), \(d_{WL}^{s}\), \(d_{r,NL}\), and \(d_{NL}^{s}\) are WL and NL UPDs at the receiver and satellite, respectively. From (2), either the estimated IF ambiguity or L1 and L2 ambiguities can be decomposed into the corresponding WL and NL ambiguities. If the float estimates of all the ambiguities are obtained, the UPD can be estimated according to (3) for WL and NL, respectively (Ge et al. 2008). With the estimated WL and NL UPDs, the ambiguities can be fixed by the user. It should be mentioned that more accurate tropospheric delays can be obtained by fixing the ambiguity at the server for better tropospheric delay modeling (Lu et al. 2018).

Specifically, the total slant tropospheric delay at an elevation angle \(el\) and azimuth angle \(az\) can be written as:

$$\begin{aligned} T_{r} \left( {el,az} \right) & = mf_{h} \left( {el} \right) \cdot {\text{ZHD}} + mf_{w} \left( {el} \right) \cdot {\text{ZWD}} \\ & + mf_{g} \left( {el} \right)\left( {\cos \left( {az} \right) \cdot G_{n} + \sin \left( {az} \right) \cdot G_{e} } \right), \\ \end{aligned}$$
(4)

where \(mf_{h} \left( {el} \right)\), \(mf_{w} \left( {el} \right)\), and \(mf_{g} \left( {el} \right)\) are mapping functions of \({\text{ZHD}}\), \(\mathrm{ZWD}\), and horizontal gradients, respectively; \({G}_{n}\) and \({G}_{e}\) are the horizontal gradients in the north and east directions, respectively. Generally, the ZHD can be corrected using the GPT2w model (Böhm et al. 2015) and the Saastamoinen model (Saastamoinen 1972), whereas the ZWD and gradients have to be estimated together with the user position during the data processing.

As for the ZHD, it can be precisely calculated given the accurate atmospheric pressure. On the other hand, if the empirical tropospheric delay model such as GPT is used to provide the a priori ZHD, the precision is usually at the level of less than 3 cm, and the residual ZHD goes into the ZWD estimates. In this case, we address that the consistent a priori ZHD model must be used at the server and user ends. The a priori ZHD modeling information will also be broadcasted to the users as keywords.

Optimal fitting coefficients model

In PPP augmentation, it is more feasible to broadcast the fitted coefficients to users than to broadcast the gridded corrections or nearby station corrections since the first one is suitable in real-time communication. Usually, in real-time PPP server processing, a second-order polynomial model is usually adapted (Shi et al. 2014).

$$\begin{aligned} {\text{ZWD}}\left( {dB,dL,h} \right) & = a_{0} + a_{1} \cdot dB + a_{2} \cdot dL + a_{3} \cdot h + a_{4} \cdot dB \cdot dL + a_{5} \cdot dB \cdot h \\ & + a_{6} \cdot dL \cdot h + a_{7} \cdot dB^{2} + a_{8} \cdot dL^{2} + a_{9} \cdot h^{2} , \\ \end{aligned}$$
(5)

where the terms (\(a_{0} ,a_{1} , \ldots ,a_{9}\)) represent the fitting coefficients; \(dB\) and \(dL\) are the latitude and longitude differences of the site with respect to the reference point, respectively, and \(h\) is the ellipsoid height. The internal quality metric for the model is the root mean square (RMS) of ZWD fitting residuals on all reference stations. The criterion for obtaining the optimal fitting coefficient is to minimize the modeling RMS.

It must be noted that the height-related tropospheric delay correction in the OFC model is a mathematical approximation, and the physical variation of tropospheric delay has not been taken into consideration (Dousa and Elias 2014). In fact, the average pressure structure of the atmosphere in the vertical direction can be modeled by an exponential decay function as a consequence of hydrostatic balance. Therefore, the ZWD variations are likely to change exponentially along with the height. On the other hand, it might be difficult to achieve an effective fit for areas with large altitude variations with a polynomial model.

Modified OFC model

To better approximate the altitude-related variations of ZWD, the OFC model is modified by introducing the exponential function. The modified OFC model, referred to as MOFC hereafter, is proposed as:

$$ZWD\left( {dB,dL} \right) = a_{0} + a_{1} \cdot dB + a_{2} \cdot dL + a_{3} \cdot dB \cdot dL + a_{4} \cdot dB^{2} + a_{5} \cdot dL^{2} ,$$
(6)
$$ZWD\left( {dB,dL,h} \right) = ZWD\left( {dB,dL} \right) \cdot e^{{\left( \frac{h}{H} \right)}} ,$$
(7)

where the second-order polynomial models the horizontal variation, the function \({\text{e}}^{{\left( \frac{h}{H} \right)}}\) describes the change in height, and \(H\) is the ZWD scale height which describes the altitude dependence of the wet delay and calculated in the least-squares estimation in modeling (Dousa and Elias 2014). The number of the coefficients of this model is reduced from ten to seven compared to the OFC model with the second-order because the exponential term replaces the height-related terms in the OFC model.

During the parameter estimation, gross errors should be effectively detected and eliminated to achieve the high-accuracy in modeling. Before estimation, the quality control in preprocessing is divided into two steps. The first step is to select the station that has a high ambiguity fixing rate as the reference station for modeling. In this study, we chose the station which has an ambiguity fixing rate of higher than 95%. The second step is based on the median detection method to verify extracted ZWDs on each station. Since the model coefficients are broadcasted every 5 min, we choose the data from the previous 5 min of the fitted time point for modeling. Therefore, to obtain stable data for modeling, the ZWDs of multiple epochs will be tested and outliers will be rejected. To ensure the accuracy and numerical stability of the estimation, the optimal parameters could be obtained by iteration. The model criterion is based on minimizing the RMS of the model fitting residuals.

In each iteration, the model coefficients in (6) and (7) are estimated in the least-squares adjustment, where the ZWD estimates at each station serve as the observation. If the residual of one observation is larger than three times the threshold, it is eliminated as an outlier. The threshold is different at different station spacing. Once there are no more outliers and the RMS of the fitting residuals is within the criteria, the iteration stops. Note that we define a maximum iteration number to be ten. If the convergence cannot be achieved within ten iterations, the modeling at this epoch is considered as failed and we broadcast the modeling coefficients of the last successful epoch.

In our study, convergence can always be achieved within six iterations. For the RMS threshold, the value is adjusted by the station-spacing and the number of stations involved in the modeling. From our experience, we set the threshold as 4 mm, 5 mm, and 6 mm for 100 km, 150 km, and 200 km station-spacing, respectively.

The delays from all reference stations are used at the server to estimate and generate the wide-area ZWD MOFC model, as expressed by (7). The model coefficients are broadcasted to the user where the ZWD will be calculated and added as a virtual observation to enhance the PPP solution,

$$v_{{{\text{zwd}}}} = {\text{ZWD}}_{W} - {\text{ZWD}}_{{\text{w,fit}}} , \sigma_{{{\text{model}}}}$$
(8)

where \({\text{ZWD}}_{{\text{w,fit}}}\) is the ZWD value calculated with the received coefficients and \(v_{{{\text{zwd}}}}\) are the residuals, \(\sigma_{{{\text{model}}}}\) is the a priori variance factor which serves as the constraint.

Figure 1 shows the data processing of this paper. The UPD WL and NL products are calculated by orbit, clock products, and observation data with globally distributed stations. Once the UPD products are obtained, stable and accurate atmospheric delays can be extracted by PPP with a fixed solution in a wide area of stations. Eventually, the extracted atmospheric delays are modeled and broadcasted to the user as coefficients for fast and high accuracy positioning.

Fig. 1
figure 1

Flowchart of tropospheric delay modeling and PPP augmentation processing. Globally stations estimated UPD products are used to achieve an ambiguity fixed solution in regional stations to extract atmospheric delays. User-side can achieve fast resolution ambiguity after obtaining the augmentation information

Full size image

Experimental validations

In this section, the data collection and processing strategies for the experiments are first presented, and the performance of the troposphere model is then assessed.

Data set

To evaluate the precision of the proposed method and its dependency on the station density of the reference network, the EPN network (Bruyninx et al. 2019), including its densification with about 460 stations in total is used. The station distribution is shown in Fig. 2. Three reference networks that will be used to generate the tropospheric correction models are defined by selecting stations with an averaged inter-station distance of 100 km, 150 km, and 200 km, marked by blue, green, and red solid cycles, respectively. The number of stations for the three reference networks are 260, 165, and 100, respectively. The rest 200 stations marked by purple stars are used as user stations for assessing the positioning performance.

Fig. 2
figure 2

EPN network, including its densification with about 460 stations. The blue, green, and red circles denote the reference network with the station-spacing of about 100 km, 150 km, and 200 km, respectively, with a corresponding station number of 260, 165, and 100. For PPP validation 200 stations in the purple star are used

Full size image

The three networks represent a dense, medium, and sparse network, respectively, and the number of stations of the latter two are reduced by 37% and 62% with respect to the dense network. As shown in Fig. 2, all three networks show an approximate uniform distribution and can cover the region of Europe well.

In total, data over 7 weeks in 2020 from the EPN network are used, including four weeks in winter, doy 001 to 028, and one week in each of the other seasons, i.e., for spring from doy 090 to 096, summer from doy 180 to 186, and autumn from doy 270 to 276. We use observations across different seasons to fully consider the different water vapor content and validate our method.

The major part is to perform the PPP with ambiguity resolution on each reference station to extract the tropospheric delay and estimate the OFC and MOFC models coefficients. After receiving the coefficients on the user, the PPP-AR with troposphere augmented can be performed. The tropospheric delay parameters are updated every 5 min, which is enough to model the temporal variations. At the server, the parameter is set as a random walk process to estimate, while at the user, tropospheric delay obtained by the model can be used as an a priori constraint. The OFC and MOFC models are obtained based on the estimated ZWDs at the reference stations. The fitting RMS and the PPP-AR accuracy and convergence time are calculated and statistically investigated. The data processing strategy is shown in Table 1.

Table 1 PPP data processing strategies
Full size table

Analysis of the tropospheric delay modeling precision

It is well known that the correlation between tropospheric delay and altitude is significant. In this study, among all stations, the highest one is PIMI (2923.43 m, located in Bareges La Mongie, France) and the lowest one is VERG (29.99 m, located in Vergi, Estonia), indicating a significant altitude difference of 2893.44 m. The corresponding ZWD differences between PPP extracted values and those fitted from OFC and MOFC models are presented in Fig. 3.

Fig. 3
figure 3

Station-specific RMS of the ZWD differences between PPP estimated and that calculated from OFC and MOFC modeling using the dense (top left), medium (top right), and sparse (bottom left) station-spacing in Europe on doy 001, 2020. Bottom right shows the number of stations in each altitude interval

Full size image

According to Fig. 3, most stations are below 500 m. Only 51 stations are between 500 and 1000 m, 22 within 1000–2000 m, and 3 above 2000 m. Worth mentioning that in the regional modeling, the selected reference stations cover those above 2000 m (except for station JANU) to improve the model precision at high altitude. The overall RMS of ZWD from the OFC model is about 15.4 mm, 18.3 mm, and 23.7 mm for networks with 100 km, 150 km, and 200 km station-spacing, respectively. The RMS of the OFC model shows a significant linear increase with the station-spacing, especially for the larger station-spacing network. On the other hand, in the MOFC model, the precisions are improved to about 13.9 mm, 15.8 mm, and 16.5 mm for the 100 km, 150 km, and 200 km networks, respectively. Nevertheless, with the MOFC model, the linear trend of the RMS with respect to station height disappears and the modeling precision for the latter two reference networks are almost the same.

We divided all stations into four groups according to their altitudes to further investigate the correlation between fitting accuracy and station altitude. The mean RMS values of the two models for seven weeks in different heights ranging from 0 to 3000 m are shown in Fig. 4.

Fig. 4
figure 4

Tropospheric delay modeling precision of OFC and MOFC at different altitude intervals using 460 stations during seven weeks cross four seasons

Full size image

From Fig. 4, the RMS of the OFC and MOFC are similar for stations under 500 m altitude and they increase along with the altitude, but with a very different rate: gently for MOFC and rapidly for OFC. The RMS of OFC reaches values larger than 60 mm for stations with an altitude over 2000 m, while that of MOFC is always below 18 mm for all the three reference networks, indicating an improvement of 70%–80%. These results reflect that we cannot achieve high-precision tropospheric modeling in areas with considerable altitude changes with the OFC model. Besides, ZWD modeling from MOFC using the dense, medium, and sparse network configurations has a similar accuracy. Because the MOFC model has better modeling of the altitude dependence, it is possible to reconstruct a ZWD correction model with homogeneous accuracy for stations at different altitudes. Figure 5 further illustrates the improvements of the OFC and MOFC models with respect to the PPP estimated values during all test periods.

Fig. 5
figure 5

Improvements of three types of station-spacing in different altitude stations between OFC and MOFC model. The different colors represent different station altitudes

Full size image

From Fig. 5, both models show similar accuracy for the stations below 500 m of three networks. As the altitude increases, the RMS of OFC model increases dramatically faster than that of the MOFC model. It is evident that the significant improvement is in higher altitude areas. The mean RMS for all assessed periods of the dense, medium, and sparse networks are 16.2 mm, 17.1 mm, and 19.1 mm for the OFC model, and 13.5 mm, 14.4 mm, and 15.0 mm for the MOFC model, respectively. In general, compared with the result of the traditional OFC model, which has been widely used in stand-alone network-RTK and PPP-RTK mode, the modeling precision of MOFC improves by 15.6–21.4% in different station-spacing networks. In addition, the differences in temperature and humidity by inter-seasonal variability can cause ZWD fluctuation. However, compared with modeling and positioning precision, with the high-frequency model coefficients calculated and broadcast, the inter-epoch tropospheric delay variation is very small. Hence, with the rapid update of augmentation information from the server, the residual ZWD at the user does not present significant variation in different seasons. The ZWD delays and inter-epoch difference values on 213, 2020 at three different altitude stations are given in Fig. 6.

Fig. 6
figure 6

ZWD values (in black) and the between-epoch differences (in blue) for three stations with different altitudes on doy 213, 2020. Top ARGI (62.0º N, − 06.8º E, 110.20 m), middle: GRAS (43.8º N, 06.92º E, 1318.43 m), bottom ESCO (42.5º N, 00.98º E, 2508.43 m)

Full size image

Performance of tropospheric delay augmented PPP-AR

To quantify the impacts of the improved troposphere model, we processed the 200 selected stations (purple stars in Fig. 2) in GPS-only, Galileo-only, and combined GPS + Galileo kinematic PPP-AR solutions augmented by the OFC and MOFC models determined from the three reference networks. In addition, the float and fixed solutions without troposphere augmentation are also calculated and compared. The station coordinates from average daily GPS + Galileo static PPP solutions are used as reference. Forming a single inter-satellite difference by selecting a reference satellite can eliminate receiver-side biases to resolution PPP ambiguity on a single station. In the real-time PPP processing with ambiguity fixing, since the solution is performed unidirectionally epoch-by-epoch, the first successful fixing is considered when the ambiguity parameter can be continuously and successfully fixed (Feng and Wang 2008; Teunissen 2001). In addition, we consider the successfully fixed when the ratio value larger than the threshold we set and the number of satellites fixed in a single epoch exceeds 75%. Hence, the convergence time is defined as the time to first fix (TTFF), which is achieved when the ambiguity can be fixed continuously. The a priori variance of augmented tropospheric ZWD is set according to the RMS of modeling.

The estimated ZWD models with three types of station-spacing are used for the kinematic PPP test on EMM1 (32.53 m) and JANU (2583.89 m) stations, which have been marked in Fig. 2. It should be noted that these two stations are not used in the modeling and thus can serve as an independent reference. Figure 7 (Galileo-only), Fig. 8 (GPS-only), and Fig. 9 (GPS + Galileo) present the results of the first two hours (240 epochs) position differences in the north (top panel), east (middle panel), and up (bottom panel) components, respectively. In Figs. 7, 8, and 9, positioning observation data are collected from 12:00 to 14:00 on doy 180, 2020.

Fig. 7
figure 7

Position differences of kinematic PPP using Galileo-only constellations in the north (top), east (middle), and up (bottom) components for station EMM1 (left panels) and JANU (right panels). The solutions are float PPP, standard PPP-AR, and troposphere augmented PPP-AR solutions

Full size image
Fig. 8
figure 8

Position differences of kinematic PPP using GPS-only constellations in the north (top), east (middle), and up (bottom) components for station EMM1 (left panels) and JANU (right panels) with float PPP, standard PPP-AR, and troposphere augmented PPP-AR solutions

Full size image
Fig. 9
figure 9

Position differences of kinematic PPP using combined GPS + Galileo constellations in the north (top), east (middle), and up (bottom) components for station EMM1 (left panels) and JANU (right panels) with float PPP, standard PPP-AR, and troposphere augmented PPP-AR solutions

Full size image

For the GPS-only and Galileo-only solutions shown in Figs. 7 and 8, the MOFC-100 km solution has the best performance in TTFF and in precision. Both OFC-200 km and float solutions have relatively poor performance with a longer convergence time. In station EMM1 (left), the fixed solutions take about 35 and 30 min, while it can be improved to 25 and 25 min by MOFC model using Galileo and GPS constellations, respectively. The effectiveness of MOFC model is especially remarkable in JANU (right) station on the up component, while OFC solutions are degraded due to the large modeling errors. Compared to the Galileo-only solution (Fig. 7), the GPS-only solution (Fig. 8) performs better in all results. This is due to the fact that GPS currently has more satellites than Galileo, resulting in a more robust observation geometry.

As for the GPS + Galileo combined PPP-AR solution shown in Fig. 9, the north and east components of station EMM1 with troposphere augmentation converge slightly faster than those without augmentation. Compared to the solution without augmentation, which takes up to 30 min for the first fix, the TTFF of OFC and MOFC is 25 min and 20 min, respectively. On the other hand, the TTFF and accuracy on the JANU station north and east components show worse performance in the OFC-150 and 200 km solutions. In the up component, the 200 km solutions converged after 40 min, while it takes only 20 min for the MOFC models using different networks.

In addition, Fig. 10 shows the residuals of ZWD estimates in PPP-AR solutions with OFC and MOFC augmentations at EMM1 (left) and JANU (right) stations on doy 1, 2020. As we can see, the three types of OFC and MOFC models have comparable residuals at station EMM1, but for the station JANU with a higher altitude, the OFC models have extremely large residuals up to 15 cm. Hence, the results indicate that the MOFC model can well fit the altitude-dependent tropospheric ZWD variations to reduce the convergence time up to 15 min, especially for the high-altitude station.

Fig. 10
figure 10

The residuals between ZWD modeling values obtained from OFC and MOFC models and PPP estimates on three types of station-spacing. Top is EMM1 station, and bottom is JANU station

Full size image

To further evaluate the MOFC model in different station-spacing networks, the kinematic daily positioning accuracy and TTFF of the 200 stations over seven weeks are further calculated and shown in Fig. 11.

Fig. 11
figure 11

Kinematic position precision and average TTFF statistics with float PPP, standard PPP-AR, and troposphere augmented PPP-AR solutions in seven weeks across four seasons. The top and bottom lines (red line) in error bars denote the max and min variability of results after outlier detection and removal. The horizontal axis shows the different PPP solutions, “Float” and “Fixed” for PPP without and with ambiguity resolution, and the rest for the PPP-AR using tropospheric augmentation with different tropospheric delay modeling method, including OFC and MOFC models with 100 km, 150 km, and 200 km inter-station distances

Full size image

For the GPS + Galileo combined fixed solution, the TTFF is about 31 min and positioning accuracies are 2.3 cm, 2.2 cm, and 3.6 cm in the north, east, and up components, respectively. When applying ZWDs from OFC modeling using dense and medium network configurations, the corresponding TTFF improvements are 10 min (33.0%) and 2 min (6.4%) for combined GPS + Galileo solutions. However, using the OFC sparse network with worse ZWD model accuracy, the TTFF is slightly increased. The positioning accuracy also worsens with the RMS values of 3.1 cm, 2.9 cm, and 4.6 cm on the north, east, and up components. Among them, the combined performance is significantly better than the single system solution.

Nevertheless, using the a priori ZWD from the MOFC model can improve the PPP-AR performance no matter which network configuration is used. Compared to the OFC model, when MOFC models with dense, medium, and sparse networks are used, the corresponding TTFF improvements are about 1 min (5.6%), 8.4 min (28.3%) and 13.6 min (38.8%), respectively. The gain in positioning accuracy on the north and east components with respect to using OFC are slightly improved by about 0.1 cm and 0.2 cm when derived from the dense network, 0.3 cm and 0.5 cm from the medium network, and 0.9 cm and 1.2 cm from the sparse network, respectively. It indicates that the high accuracy regional tropospheric corrections using MOFC model can be applied as reliable constraints to enhance PPP performance.

From these results, the MOFC shows significantly better performance than the OFC models. The OFC model with relatively larger modeling errors in regions with large altitude variations degrades the PPP-AR solutions, whereas using the MOFC model for regional tropospheric delay augmentation, the PPP-AR solutions are always enhanced in both the precision and the convergence time in different station-spacing networks.

Moreover, Fig. 12 presents the improvement of the combined GPS and Galileo positioning compared to processing with a single constellation. For Galileo, the solutions without additional tropospheric information are improved by 5.9% after being combined with GPS observations. For the MOFC model, the improvement is significantly greater than that of the OFC model. Besides, the effect of the combination is lower for GPS improvement. This is because GPS has more robust observation geometry with more satellites than Galileo, and as a result, the relative improvement if more significant for Galileo than GPS.

Fig. 12
figure 12

Relative precision improvement of GPS + Galileo combined solution with respect to the GPS-only and Galileo-only solutions using the kinematic PPP approach

Full size image

Conclusions

In this study, an improved method for generating a tropospheric delay correction model to augment PPP-AR over wide areas and high altitude variations is developed, where the exponential function is adopted to model the altitude dependence of tropospheric delay. We evaluate the troposphere modeling accuracy in the Europe region with large altitude variations and further analyze the performance in Galileo and GPS PPP-AR.

We use seven weeks of GPS + Galileo observations for validation, collected from about 460 stations of the EPN and EPN densification networks and covering the four seasons. We investigate the modeling precision of different network-spacing, including dense (100 km), medium (150 km), and sparse (200 km) networks. The proposed MOFC model demonstrates improved precision than the traditional OFC model, especially for the high-altitude stations. Compared to OFC, the proposed MOFC improves modeling precision by 16.1%, 15.5%, and 21.4% for dense, medium, and sparse networks. For different altitudes, the improvements are 5%, 20%, 55%, and 70% for 0–500 m, 500–1000 m, 1000–2000 m, and 2000–3000 m, respectively. Moreover, the kinematic positioning accuracy using MOFC improves up to 45%, and the TTFF reduces up to 12 min.

The proposed model is capable of achieving high accuracy tropospheric modeling and positioning in a wide area with large altitude variations using fewer stations and coefficients, which satisfy the current real-time precise positioning service.