New Adaptable All-in-One Strategy for Estimating Advanced Tropospheric Parameters and Using Real-Time Orbits and Clocks

We developed a new strategy for a synchronous generation of real-time (RT) and near real-time (NRT) tropospheric products. It exploits the precise point positioning method with Kalman filtering and backward smoothing, both supported by real-time orbit and clock products. The strategy can be optimized for the latency or the accuracy of NRT production. In terms of precision, it is comparable to the traditional NRT network solution using deterministic models in the least-square adjustment. Both RT and NRT solutions provide a consistent set of tropospheric parameters such as zenith total delays, horizontal tropospheric gradients and slant delays, all with a high resolution and optimally exploiting all observations from available GNSS multi-constellations. As the new strategy exploits RT processing, we assessed publicly precise RT products and results of RT troposphere monitoring. The backward smoothing applied for NRT solution, when using an optimal latency of 30 min, reached an improvement of 20% when compared to RT products. Additionally, multi-GNSS solutions provided more accurate (by 25%) tropospheric parameters, and the impact will further increase when constellations are complete and supported with precise models and products. The new strategy is ready to replace our NRT contribution to the EUMETNET EIG GNSS Water Vapour Programme (E-GVAP) and effectively support all modern multi-GNSS tropospheric products.


Introduction
The exploitation of Global Navigation Satellite System (GNSS) observations for monitoring the troposphere in support of meteorology has been proposed by [1].The initial parameter of interest was the Zenith Total Delay (ZTD), which represents a GNSS signal path delay due to transecting the neutral atmosphere in a vertical column above the station.The ZTD depends mainly on atmospheric pressure and partial water vapor pressure [2].First computation of ground-based GNSS ZTDs in near-real time (NRT) were demonstrated in 2000 [3,4], i.e., shortly after establishing GNSS hourly data flow in 1999 and after providing precise orbits in ultra-rapid mode by the International GNSS Service (IGS, http://www.igs.org) in 2000 [5].The COST Action 716 [6] then played an important role in developing and evaluating methods of GNSS NRT troposphere monitoring and defining the standard format (COST-716) for the product dissemination.The operational production of GNSS ZTDs was organized within the COST-716 Demonstration Project [7], and it has never been closed because its coordination was handed over to the newly established EUMETNET EIG GNSS Water Vapour Programme (E-GVAP, http://egvap.dmi.dk) in 2004.Operational assimilations of ZTDs from E-GVAP into mesoscale Numerical Weather Models (NWM) for weather forecasting have been performed by Météo-France [8,9] and UK Met Office [10] since 2007.Both meteorological agencies started also an assimilation of the first E-GVAP global NRT ZTD product [11].Until present days, the standard GNSS ZTD monitoring in E-GVAP is performed in hourly update rate with a maximum product delay of 45 min with respect to the last GNSS processed observation.These requirements are a prerequisite of the ZTD assimilation performed usually with a 3-6 h time resolution.
Until now, a majority of E-GVAP analysis centers (ACs) uses a double-difference observation processing in a network solution.This strategy eliminates clock errors at GNSS receiver and satellite and was compulsory while public products were not available in NRT.The situation has changed in 2013, when the IGS introduced the Real-Time Service (RTS, http://rts.igs.org)providing GPS and GLONASS orbit and clock corrections by combining contributions from several IGS real-time analysis centers [12].The IGS RTS aims at supporting real-time (RT) analyses with the Precise Point Positioning (PPP) method [13].The PPP is based on original observations or their linear combination without differencing between receivers or satellites.Though German Research Centre for Geosciences (GFZ) has provided a NRT PPP ZTD product [14,15] since 2001, it was possible only thanks to their two-step processing approach consisting of (1) a global NRT solution for determining consistent satellite clock and orbit products and (2) a distributed PPP processing for ZTD estimated at each station individually.With the availability of global real-time data flow, software and standards specified for precise product dissemination, the PPP is becoming more popular for the troposphere monitoring.Compared to the traditional approach in E-GVAP dominated by the double-difference network processing, the PPP offers several advantages: (a) an easy production in real-time or NRT fashion, (b) flexible use of central or distributed processing scheme including a receiver built-in solution, (c) an estimation of tropospheric parameters in the absolute sense with a high spatio-temporal resolution, and (d) an optimal support of all satellite constellations and new signals including multiple frequencies; all profiting from a highly efficient and autonomous processing approach.The price for mentioned advantages is however paid by several disadvantages.Compared to the strategy using double differences, all observation models need to be carefully applied to reach the best accuracy.In addition, integer ambiguity resolution is possible only if precise observation phase biases are available, thus often non-integer-fixed ambiguities are usually estimated.
First ZTD products using the PPP method and new IGS RTS products were demonstrated in 2014 [16][17][18].Though the ZTD estimates already reached an internal precision below 10 mm when compared to GNSS final products, the studies also revealed significant station-/product-specific systematic errors attributed to inconsistencies in precise models applied in PPP and to the precise product generation.Fortunately, a large portion of the systematic errors are changing slowly over time and are thus not critical for an assimilation into NWMs which were designed to identify and remove biases on monthly basis by comparing station/product-specific GNSS ZTDs with their counterparts from the model background information [10].However, the precision of the real-time ZTDs obtained using the Kalman filter still remained worse by a factor of 1.5 compared to the E-GVAP standard NRT ZTD products.The reasons were twofold: (1) dependence on the quality of the RT products, and (2) use of the strategy for real-time data analysis.The E-GVAP solutions utilizing the batch processing and data files can be characterized with a standard deviation of 3-6 mm and 3-8 mm, for regional and global products, respectively, and a systematic error within 1-3 mm [19].The accuracy evaluated with respect to external data, such as radiosonde, corresponds to 1-2 mm in the precipitable water vapor [15].Anyway, real-time tropospheric products are ready to support assimilation in the rapid update cycle of NWM prediction or nowcasting, and target a short-term weather prediction or severe weather events monitoring [20].
Nowadays, multi-GNSS offers many satellites and signals that are expected to strengthen all estimated parameters, in particular the ZTD and horizontal tropospheric gradients.For the E-GVAP production, data from the US NASTAR Global Positioning System (GPS) have been initially used.In 2008, the IGS started a provision of ultra-rapid precise orbit products for the Russian GLONASS satellites too.However, GLONASS data could have been included in NRT analyses only after resolving a 1.5 mm systematic difference in the ZTD when estimated independently from GPS and GLONASS data [21].The bias was caused by using the IGS05 model of phase center offsets (PCO) for GLONASS satellite antennas [22], and the problem has been eliminated in 2012 [19] by adopting consistent PCO models for both constellations in new IGS08 realization [23].Besides others, general limitations for the use of GNSS data from other global systems, European Galileo and Chinese BeiDou, persist mainly in (1) incompleteness of the constellations, (2) lack of precise models and calibrations for new signals, receiver and satellite instrumentations, and (3) lack of precise orbit and clock products supporting the ultra-fast processing mode.The situation will change soon as both global systems will become operational in next years.Since 2012, the IGS Multi-GNSS Experiment (MGEX, http://mgex.igs.org) has been successfully filling the gaps in data, metadata, models, formats, standards and products for an optimal exploitation of all global satellite constellations and other regional augmentations.Although Galileo and BeiDou systems have not been completed yet, several groups reported an initial positive impact when using multi-constellation for estimating ZTDs and horizontal tropospheric gradients [24][25][26][27].
The GNSS ZTD estimated from a single hour of data is not stable and accurate enough for the E-GVAP usage due to high correlations with coordinates and initial phase ambiguities.At least 12-h data batch is usually used for the NRT processing in E-GVAP, which means a high degree of redundancy in data processing updated on an hourly basis.In past, the issue was usually solved by reducing the data processing window to 1-6 h and by combining normal equations for the tropospheric parameter estimation [11,28].It was also reported, that ZTD values at the edges of the processing interval are about 20% less accurate compared to those in the middle of the interval and that the main impact is actually observed during the first and last hours of the interval [29].The E-GVAP expects estimated ZTDs within the last (production) hour, thus a trade-off between the accuracy and the latency is important.
A piece-wise linear function for modeling the tropospheric parameters within each hour is used by a majority of E-GVAP contributors using the Bernese GNSS Software V52 [30].In addition, horizontal tropospheric gradients are usually not estimated in operational solutions because of two reasons.First, there is presently no operational assimilation of gradients into NWM.Second, the estimation of high-resolution horizontal gradients reflecting a high spatio-temporal variability of local humidity increases the number of estimated parameters in the network and, consequently, the computation time by a factor of 2-3 at least.The PPP with an epoch-wise filtering supported by the IGS RTS products is an optimal strategy for generating advanced GNSS tropospheric products for future meteorological applications [20,31] such as high-resolution ZTDs, horizontal tropospheric gradients, and slant tropospheric delays.
In this paper, we present a new processing strategy adaptable in the operational mode when prioritizing product latency (RT) or product accuracy (NRT).Additionally, all the advanced tropospheric parameters for RT and NRT products are derived within a single continuously operating RT PPP engine supported by the IGS RTS orbit and clock corrections.We believe that such unified processing system will replace soon our existing NRT contribution to E-GVAP.The strategy exploits all abovementioned advantages including benefits of full multi-constellation in near future.
The paper is organized as follows: after the introduction we assess real-time precise orbit and clock products available from the IGS RTS, evaluate an operational prototype of GNSS ZTD real-time production, both as pre-requisites for the new strategy.In Section 3, we introduce the new strategy for an adaptable PPP processing solution for estimating ZTDs, horizontal tropospheric gradients and tropospheric slant delays using precise real-time products.In Sections 4 and 5, we assess the results of estimated tropospheric parameters driven by the backward data smoothing including study of impacts of IGS real-time products and multiple constellations.In the last section, we close the paper with conclusions.

Assessment of Available RT Orbit and Clock Products and RT ZTDs
The accuracy of the real-time ZTD calculated with the PPP method strongly depends on the quality of real-time GNSS orbit and clock corrections [16,32,33].Our new strategy for the troposphere monitoring is based on global RT products and analysis, thus we first evaluate publicly available global RT products for its support.Second, we summarize our ZTD contributions to the Real-time Demonstration Campaign initiated in 2015 by the COST Action ES1206 [20]-Advanced global navigation satellite systems tropospheric products for monitoring severe weather events and climate (GNSS4SWEC, http://gnss4swec.knmi.nl),later coordinated by the real-time troposphere monitoring Working Group 4.3.7 of the International Association of Geodesy (IAG, http://www.iag-aig.org).Third, we studied the impact of IGS RTS products for the simulated RT ZTD estimates within the GNSS4SWEC Benchmark Campaign [31].

Assessment of Real-Time Orbit and Clock Corrections
We investigate the performance of four real-time products being collected and archived at Geodetic Observatory Pecny (GOP) using the BNC Software [34] since 2013.Three of them (IGS01, IGS02, IGS03) are the official IGS RTS combined products [12], while the last one (CNS91) is provided by individual CNES RT analysis center [35].Two different strategies and software are used for combining the IGS RTS [33].While IGS01 is generated on the basis of epoch-wise combination, the Kalman filer technique is exploited for producing IGS02 and IGS03.Two different agencies are responsible for these RT products: European Space Agency (ESA) provides IGS01 and Bundesamt für Kartographie und Geodäsie (BKG) IGS02 and IGS03.All mentioned streams include orbits/clock corrections to GPS satellites, and only IGS03 supports also GLONASS constellation.Navigation data from MGEX [36] were used together with the RT corrections to recover the precise satellite orbit and clocks.First, the availability and completeness of RT corrections were checked and, second, satellite orbit and clocks were compared to IGS final orbit and clocks, both during the period of 2013-2017.
Figures 1 and 2 depict monthly completeness of RT corrections for all GPS satellites from the IGS01 and IGS02 combined products.The IGS03 product is not shown here because it is almost identical to IGS02; slightly more outages are related to GLONASS satellites.Figure 3 shows then the same picture for the CNES product.From the comparison, we can classify problems into three groups: (1) temporal unavailability period of some satellites, e.g., G03, G04, (2) source-specific unavailability, e.g., G01 for CNS91, and (3) satellite-specific incompleteness.The first group is usually caused by the loss of observations due to the upgrade of a satellite, such as replacing the old Block IIA satellites with the new Block IIF satellites or a maintenance identified by satellite unhealthy status.The second and third groups of gaps are caused by data unavailability from a global network and the processing strategy including outlier detection in the product generation.The availability of the corrections is significantly lower for some months (June 2015, December 2016) compared to others, which was caused by the Internet connection failures at GOP when receiving the streams.The source-specific loss of data at IGS02 and CNS91 streams are visible in June 2015 and, these are mainly due to the inconsistent navigation message Issue of Date (IOD) available from the MGEX broadcast and those referred by RT corrections.It can be thus recommended to use consistent RT navigation data and precise correction streams optimally guaranteed by the same provider.In general, the availability of RT corrections is well over 90% for most satellites, which agrees with findings in [33].It indicates that the RT corrections were provided continuously for use in troposphere monitoring, however, problems can be expected in a kinematic positioning which is more sensitive to the product incompleteness.
Apart from the availability of the corrections, the precision is critical for the user performance.The orbits are compared in 5-min intervals for three components: radial, along-track and cross-track while the clock comparison is based on the second order difference method.The IGS08 and the IGS14 model is used to correct satellite PCOs prior and after 29 January 2017, respectively, corresponding to the adoption of the IGS2014 reference frame [37].The clock datum is estimated by calculating a mean over all satellites clocks at each epoch.The datum inconsistencies are then eliminated through single-differences between individual satellite clocks and the clock datum.The single-differences from the real-time clocks are compared to those from the IGS final product.The root-mean-square errors (RMSE) of RT orbits and clocks are calculated for each day while outliers are removed using a fixed threshold.Although there is a strong correlation between clocks and radial orbit component, we haven't corrected this dependency.Table 1 gives summary statistics for all products over all days.
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 22 root-mean-square errors (RMSE) of RT orbits and clocks are calculated for each day while outliers are removed using a fixed threshold.Although there is a strong correlation between clocks and radial orbit component, we haven't corrected this dependency.Table 1 gives summary statistics for all products over all days.root-mean-square errors (RMSE) of RT orbits and clocks are calculated for each day while outliers are removed using a fixed threshold.Although there is a strong correlation between clocks and radial orbit component, we haven't corrected this dependency.Table 1 gives summary statistics for all products over all days.root-mean-square errors (RMSE) of RT orbits and clocks are calculated for each day while outliers are removed using a fixed threshold.Although there is a strong correlation between clocks and radial orbit component, we haven't corrected this dependency.Table 1 gives summary statistics for all products over all days.The orbit difference in radial component shows the smallest RMSE for all products, whereas the along-/cross-track components reached slightly larger values.The IGS01 orbit shows the best agreement with respect to the IGS final orbits.Largest differences are observed for the orbits from IGS03, which might be attributed to a different outlier detection method applied when including GLONASS satellites.Time evolution of the orbit comparison for each product and specific component is shown in Figure 4. Coordinate differences greater than 30 cm are plotted at the top horizontal lines of each graph.Orbits from the IGS01 stream are less affected by the outliers compared to IGS02 and IGS03 products, as indicated by outliers mainly during March 2015.The switch from the IGS08 to the IGS14 PCO model (28 January 2017) can be observed in statistics of the radial component.It seems that the CNS91 product used the new IGS14 model as of 9 March 2017, while official IGS solutions are difficult to recognize due to most likely asynchronous switches by different contributing providers.Otherwise, the orbit accuracy for all products shows an overall good consistency over the period.The orbit difference in radial component shows the smallest RMSE for all products, whereas the along-/cross-track components reached slightly larger values.The IGS01 orbit shows the best agreement with respect to the IGS final orbits.Largest differences are observed for the orbits from IGS03, which might be attributed to a different outlier detection method applied when including GLONASS satellites.Time evolution of the orbit comparison for each product and specific component is shown in Figure 4. Coordinate differences greater than 30 cm are plotted at the top horizontal lines of each graph.Orbits from the IGS01 stream are less affected by the outliers compared to IGS02 and IGS03 products, as indicated by outliers mainly during March 2015.The switch from the IGS08 to the IGS14 PCO model (28 January 2017) can be observed in statistics of the radial component.It seems that the CNS91 product used the new IGS14 model as of 9 March 2017, while official IGS solutions are difficult to recognize due to most likely asynchronous switches by different contributing providers.Otherwise, the orbit accuracy for all products shows an overall good consistency over the period.Table 1 also summarizes RMSE and standard deviation (SDEV) of the real-time clock corrections.The former represents the accuracy relevant for the processing of code pseudoranges while the latter characterizes the precision important for the carrier-phase processing.It can be also interpreted from the PPP point of view combining both observation types as follows-the former have a positive impact on the PPP convergence time while the latter enable more precise positioning within already converged solution [38].Obviously, this is the case of IGS01 and CNS91 products when the first is more accurate, but the second more precise for the PPP application.The IGS02 and IGS03 products performs slightly worse in terms of both RMSE and SDEV.
Figure 5 finally shows time series of the clock SDEV and RMSE statistics.The former (top plot) indicates a comparable high quality over the period for IGS01 and CNS91, while more outliers are  The former represents the accuracy relevant for the processing of code pseudoranges while the latter characterizes the precision important for the carrier-phase processing.It can be also interpreted from the PPP point of view combining both observation types as follows-the former have a positive impact on the PPP convergence time while the latter enable more precise positioning within already converged solution [38].Obviously, this is the case of IGS01 and CNS91 products when the first is more accurate, but the second more precise for the PPP application.The IGS02 and IGS03 products performs slightly worse in terms of both RMSE and SDEV.

Impact of IGS RTS Products on ZTD Estimates
The impact of the IGS RTS products on PPP ZTD estimates was assessed by exploiting the GNSS4SWEC Benchmark campaign with 400 GNSS stations in central Europe during the period of May-June 2013.The ZTD was calculated using the G-Nut/Tefnut software in the post-processing mode when supported with two precise products: (1) IGS final orbits and clocks, and (2) IGS RTS orbit and clock corrections.Two reference solutions were provided for the benchmark using different software and processing strategies [31].GOP used the Bernese GNSS Software and the double-difference processing (DD) and GFZ used the EPOS software and the PPP method.Statistics from the comparison of both testing solutions with respect to both reference products are given in Table 2. Generally, the results indicate a good agreement, however, the impact of the IGS RTS products (IGS01) on ZTDs is clearly visible in two aspects: (a) a common systematic error of 2.4-2.8mm, and (b) a lower precision of 13-17%.Interestingly, a better agreement in terms of SDEV is reached between 10-20% when using two PPP solutions (G-Nut/Tefnut vs. EPOS software) compared to the processing strategies (DD vs. PPP).The results also showed that input products and the processing strategy might result in a similar impact on the ZTD estimates, which can reach up to 20% in terms of accuracy.Finally, it should be noted that the PPP ZTD estimation used a stochastic model and an epoch-wise filtering method in the G-Nut/Tefnut software, while a deterministic model with the least-squares batch adjustment used in the EPOS software.

Impact of IGS RTS Products on ZTD Estimates
The impact of the IGS RTS products on PPP ZTD estimates was assessed by exploiting the GNSS4SWEC Benchmark campaign with 400 GNSS stations in central Europe during the period of May-June 2013.The ZTD was calculated using the G-Nut/Tefnut software in the post-processing mode when supported with two precise products: (1) IGS final orbits and clocks, and (2) IGS RTS orbit and clock corrections.Two reference solutions were provided for the benchmark using different software and processing strategies [31].GOP used the Bernese GNSS Software and the double-difference processing (DD) and GFZ used the EPOS software and the PPP method.Statistics from the comparison of both testing solutions with respect to both reference products are given in Table 2. Generally, the results indicate a good agreement, however, the impact of the IGS RTS products (IGS01) on ZTDs is clearly visible in two aspects: (a) a common systematic error of 2.4-2.8mm, and (b) a lower precision of 13-17%.Interestingly, a better agreement in terms of SDEV is reached between 10-20% when using two PPP solutions (G-Nut/Tefnut vs. EPOS software) compared to the processing strategies (DD vs. PPP).The results also showed that input products and the processing strategy might result in a similar impact on the ZTD estimates, which can reach up to 20% in terms of accuracy.Finally, it should be noted that the PPP ZTD estimation used a stochastic model and an epoch-wise filtering method in the G-Nut/Tefnut software, while a deterministic model with the least-squares batch adjustment used in the EPOS software.The RT ZTD from the demonstration campaign is evaluated for 18 European stations during the initial year of the GNSS4SWEC Real-time Demonstration campaign.Two GOP solutions using the IGS03 product are compared with respect to the EUREF 2nd reprocessing combined tropospheric product [39]: (1) GOPR-standalone GPS solution, and (2) GOPQ-GPS + GLONASS solution.
In Table 3, we can observe a systematic error in ZTD of about 2 mm in the long-term evaluation, similar as observed in the simulated real-time processing in the benchmark campaign, see Section 2.2.Although GLONASS observations are down-weighted by a factor of 2 in our solution in order to reflect the lower quality of GLONASS precise products, a small positive impact on the ZTD is observed in terms of mean bias (10%) and mean SDEV (7%), both calculated over 18 stations.Figure 6 shows the comparison of the GOPR solution with respect to the EUREF 2nd reprocessing combined product during the first year of the RT demonstration campaign.Monthly mean ZTD biases, standard deviations and their 1-sigma scatter calculated over all 18 stations indicate a long-term stability of the operational real-time production with a small seasonal effect in SDEV due to a less accurate troposphere modeling during the summer period [19].

Long-Term Quality of Operational RT ZTD Production
The RT ZTD from the demonstration campaign is evaluated for 18 European stations during the initial year of the GNSS4SWEC Real-time Demonstration campaign.Two GOP solutions using the IGS03 product are compared with respect to the EUREF 2nd reprocessing combined tropospheric product [39]: (1) GOPR-standalone GPS solution, and (2) GOPQ-GPS + GLONASS solution.In Table 3, we can observe a systematic error in ZTD of about 2 mm in the long-term evaluation, similar as observed in the simulated real-time processing in the benchmark campaign, see Section 2.2.Although GLONASS observations are down-weighted by a factor of 2 in our solution in order to reflect the lower quality of GLONASS precise products, a small positive impact on the ZTD is observed in terms of mean bias (10%) and mean SDEV (7%), both calculated over 18 stations.

Solution Description BIAS (mm) mean ± sdev SDEV (mm) mean ± sdev RMSE (mm) mean ± sdev GOPQ-GPS+GLO
1.8 ± 2.9 6.7 ± 1.2 7.5 ± 2.5 GOPR-GPS 2.0 ± 2.8 7.2 ± 1.0 7.9 ± 1.5 Figure 6 shows the comparison of the GOPR solution with respect to the EUREF 2nd reprocessing combined product during the first year of the RT demonstration campaign.Monthly mean ZTD biases, standard deviations and their 1-sigma scatter calculated over all 18 stations indicate a long-term stability of the operational real-time production with a small seasonal effect in SDEV due to a less accurate troposphere modeling during the summer period [19].

New Adaptable Strategy for RT and NRT Troposphere Monitoring
Optimizing the accuracy and the timeliness for the tropospheric products was the first motivation for developing a new adaptable processing strategy.The second goal aimed at effectively supporting various users by simultaneously combining RT and NRT analysis modes.The third goal was to retrieve all the advanced tropospheric parameters consistently when using a unique and flexible processing scheme and exploit optimally multi-GNSS data.Additionally, strategy targets all

New Adaptable Strategy for RT and NRT Troposphere Monitoring
Optimizing the accuracy and the timeliness for the tropospheric products was the first motivation for developing a new adaptable processing strategy.The second goal aimed at effectively supporting various users by simultaneously combining RT and NRT analysis modes.The third goal was to retrieve all the advanced tropospheric parameters consistently when using a unique and flexible processing scheme and exploit optimally multi-GNSS data.Additionally, strategy targets all the advantages of PPP when using epoch-wise analysis supported by real-time precise products, e.g., by IGS RTS.The new strategy has been implemented in the G-Nut/Tefnut software designed for estimating ZTDs, horizontal tropospheric gradients and tropospheric slant path delays.Currently, the parameters are estimated using dual-frequency ionosphere-free linear combinations of pseudorange and carrier-phase observations.
3.1.Epoch-Wise Filtering vs. Batch Processing, PPP vs. Network Approach Nowadays, a common procedure of NRT analysis in E-GVAP is based on the least-squares adjustment (LSQ) analysis and a piece-wise linear function for modeling the tropospheric parameters within the processing interval.Figure 7 depicts the estimated tropospheric parameters by black dots and the corresponding piece-wise deterministic model by dashed lines connecting parameters within each hourly product update; not necessarily connected at update boundaries.According to E-GVAP conventions, the last product value is shifted by 1 min (HR:59) in order to avoid a duplicate value with the next hour product update.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 22 the advantages of PPP when using epoch-wise analysis supported by real-time precise products, e.g., by IGS RTS.The new strategy has been implemented in the G-Nut/Tefnut software designed for estimating ZTDs, horizontal tropospheric gradients and tropospheric slant path delays.Currently, the parameters are estimated using dual-frequency ionosphere-free linear combinations of pseudorange and carrier-phase observations.

Epoch-Wise Filtering vs. Batch Processing, PPP vs. Network Approach
Nowadays, a common procedure of NRT analysis in E-GVAP is based on the least-squares adjustment (LSQ) analysis and a piece-wise linear function for modeling the tropospheric parameters within the processing interval.Figure 7 depicts the estimated tropospheric parameters by black dots and the corresponding piece-wise deterministic model by dashed lines connecting parameters within each hourly product update; not necessarily connected at update boundaries.
According to E-GVAP conventions, the last product value is shifted by 1 min (HR:59) in order to avoid a duplicate value with the next hour product update.Although temporal resolution of estimated parameters can be increased in the LSQ analysis, it means a significant increase of size of normal equations, in particular for the network processing, and consequently extending the processing time due to the inversion of the large matrix of the solution.On the other hand, the Kalman filter is a powerful method for real-time applications because state vectors and covariance matrices are calculated recursively epoch by epoch.One of the consequences is that only previous observations contribute to a current estimate in the Kalman filter, and thus the solution needs a certain time to converge.An accuracy of parameters estimated during the initial convergence, or any later re-convergence, can be then improved using the backward smoothing algorithm [40].This algorithm exploits all observations from past and future epochs for estimating the state vector at any epoch of the processed window.It improves significantly parameters such as kinematic coordinates or tropospheric parameters [41].While the Kalman filter was designed for real-time signal processing, the backward smoothing filter was developed for the post-processing to achieve the accuracy of all parameters similar to those from the LSQ solution.Important advantage of data filtering is that the parameters are estimated epoch by epoch with applying individual stochastic properties and the temporal resolution can reach up to the original data sampling.
For the real-time data processing, the PPP method has been implemented using the Extended Kalman filter in our software.In the first step of the filter, parameters are predicted via adding a particular amount of noise to diagonal elements of the variance-covariance matrix belonging to dynamic parameters.The state vector itself remains the same in this step.When new observations are available the state vector with its variance-covariance matrix is updated.Since the standard Kalman filter can suffer from numerical instabilities due to the round-off error we have implemented an alternative form of the filter using a square root of variance-covariance matrix instead of the original one.This modification guarantees that the estimated variance-covariance matrix is positive semidefinite.For the post-processing purposes, solutions from the Kalman filter are additionally smoothed by the filter running in the backward direction.Storing the current results Although temporal resolution of estimated parameters can be increased in the LSQ analysis, it means a significant increase of size of normal equations, in particular for the network processing, and consequently extending the processing time due to the inversion of the large matrix of the solution.On the other hand, the Kalman filter is a powerful method for real-time applications because state vectors and covariance matrices are calculated recursively epoch by epoch.One of the consequences is that only previous observations contribute to a current estimate in the Kalman filter, and thus the solution needs a certain time to converge.An accuracy of parameters estimated during the initial convergence, or any later re-convergence, can be then improved using the backward smoothing algorithm [40].This algorithm exploits all observations from past and future epochs for estimating the state vector at any epoch of the processed window.It improves significantly parameters such as kinematic coordinates or tropospheric parameters [41].While the Kalman filter was designed for real-time signal processing, the backward smoothing filter was developed for the post-processing to achieve the accuracy of all parameters similar to those from the LSQ solution.Important advantage of data filtering is that the parameters are estimated epoch by epoch with applying individual stochastic properties and the temporal resolution can reach up to the original data sampling.
For the real-time data processing, the PPP method has been implemented using the Extended Kalman filter in our software.In the first step of the filter, parameters are predicted via adding a particular amount of noise to diagonal elements of the variance-covariance matrix belonging to dynamic parameters.The state vector itself remains the same in this step.When new observations are available the state vector with its variance-covariance matrix is updated.Since the standard Kalman filter can suffer from numerical instabilities due to the round-off error we have implemented an alternative form of the filter using a square root of variance-covariance matrix instead of the original one.This modification guarantees that the estimated variance-covariance matrix is positive semidefinite.For the post-processing purposes, solutions from the Kalman filter are additionally smoothed by the filter running in the backward direction.Storing the current results from the forward filter is necessary for the backward smoothing.An accuracy of the state vector is improved using the special recurrent algorithm running from the last epoch with estimated parameters toward the beginning of the processing period.Numerical issues are also critical as in the case of the Kalman filter; therefore an advanced technique using singular value decomposition of the variance-covariance matrix was implemented [40].
By employing the PPP method, all GNSS data can be analyzed station by station, either at a central server with a separation into processing threads or through a decentralized manner including the processing at original observation sites.The PPP method effectively supports the utilization of multi-GNSS observations if consistent precise models and products are available for all satellites and systems.In future, our approach will be extended to a highly flexible processing of a multiple frequencies and signals when optimally exploiting raw GNSS observations [42].

Combining RT and NRT Processing Supported by Observations from Files or Streams
In E-GVAP, the NRT products are updated on hourly basis and the requirement for a delivery to the E-GVAP server at UK Met Office is 45 min after the last observation used in the analysis.The NRT LSQ batch processing, initiated every hour when obtaining a majority of data files, is thus a relevant solution for this purpose when using hourly data files.Although the backward smoothing approach is the most beneficial for the post-processing solutions, it can be effectively used in NRT applications too.We have thus used it for the new adaptable strategy combining a continuously running forward filter with a regularly triggered backward smoothing filter.The former is aimed for the estimating epoch-wise tropospheric parameters in real time indicated by white points in Figure 8.In addition, the filter solutions are stored from all epochs relevant for a possible improvement in the second step.The backward smoothing is started periodically at a pre-defined time stamp as indicated at HR:00 in the figure.It uses the initial state vector from the same epoch provided by the real-time filter for recalculating the past parameters from the Kalman filter.Obviously, such recalculation is able to significantly refine older parameters, but cannot improve parameters at the initial epochs of the backward smoothing.The length of the smoothing period can be set flexibly to reflect an actual user preference for a higher accuracy or a shorter latency of the product.In such a way, the standard E-GVAP NRT tropospheric product can also be provided on hourly basis as shown by blue points in the figure .Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 22 from the forward filter is necessary for the backward smoothing.An accuracy of the state vector is improved using the special recurrent algorithm running from the last epoch with estimated parameters toward the beginning of the processing period.Numerical issues are also critical as in the case of the Kalman filter; therefore an advanced technique using singular value decomposition of the variance-covariance matrix was implemented [40].
By employing the PPP method, all GNSS data can be analyzed station by station, either at a central server with a separation into processing threads or through a decentralized manner including the processing at original observation sites.The PPP method effectively supports the utilization of multi-GNSS observations if consistent precise models and products are available for all satellites and systems.In future, our approach will be extended to a highly flexible processing of a multiple frequencies and signals when optimally exploiting raw GNSS observations [42].

Combining RT and NRT Processing Supported by Observations from Files or Streams
In E-GVAP, the NRT products are updated on hourly basis and the requirement for a delivery to the E-GVAP server at UK Met Office is 45 min after the last observation used in the analysis.The NRT LSQ batch processing, initiated every hour when obtaining a majority of data files, is thus a relevant solution for this purpose when using hourly data files.Although the backward smoothing approach is the most beneficial for the post-processing solutions, it can be effectively used in NRT applications too.We have thus used it for the new adaptable strategy combining a continuously running forward filter with a regularly triggered backward smoothing filter.The former is aimed for the estimating epoch-wise tropospheric parameters in real time indicated by white points in Figure 8.In addition, the filter solutions are stored from all epochs relevant for a possible improvement in the second step.The backward smoothing is started periodically at a pre-defined time stamp as indicated at HR:00 in the figure.It uses the initial state vector from the same epoch provided by the real-time filter for recalculating the past parameters from the Kalman filter.Obviously, such recalculation is able to significantly refine older parameters, but cannot improve parameters at the initial epochs of the backward smoothing.The length of the smoothing period can be set flexibly to reflect an actual user preference for a higher accuracy or a shorter latency of the product.In such a way, the standard E-GVAP NRT tropospheric product can also be provided on hourly basis as shown by blue points in the figure.Initially, the new adaptable strategy was designed to use the precise orbit and clock corrections disseminated through RT streams, it can however exploit observations coming from both real time (streams) and near-real time (hourly or sub-hourly files).While the former is necessarily used in a simultaneous RT and NRT product generation, the latter is applicable for NRT only.In any case, both data flows can be mixed and analyzed for each station independently.Additional advantage of the NRT analysis utilizing the backward smoothing and RT observations may profit from the 45-min requirement in E-GVAP for the product delivery in NRT.By starting the backward smoothing shortly before the delivery request, indicated in Figure 8 by the red arrow starting at 12:40, new Initially, the new adaptable strategy was designed to use the precise orbit and clock corrections disseminated through RT streams, it can however exploit observations coming from both real time (streams) and near-real time (hourly or sub-hourly files).While the former is necessarily used in a simultaneous RT and NRT product generation, the latter is applicable for NRT only.In any case, both data flows can be mixed and analyzed for each station independently.Additional advantage of the NRT analysis utilizing the backward smoothing and RT observations may profit from the 45-min requirement in E-GVAP for the product delivery in NRT.By starting the backward smoothing shortly before the delivery request, indicated in Figure 8 by the red arrow starting at 12:40, new observations can further improve the accuracy of NRT product, especially last parameters of the NRT product and still reduce systematic errors typical for data interval boundaries.

Estimating High-Resolution ZTDs and Horizontal Gradients
The extended model for the GNSS slant tropospheric delay (δT) for a single receiver and a satellite is defined as a function of Zenith Hydrostatic Delay (ZHD), Zenith Wet Delay (ZWD) and horizontal tropospheric gradients pointing to North (G N ) and East (G E ) using an azimuth (A), which direction is counted from the North [43] (1) The ZHD and the ZWD are hydrostatic and wet contributions to the ZTD which represents a vertical delay due to the neutral atmosphere in a profile above the station.Note that ZHD, ZWD and horizontal gradients are projected to any satellite direction by using m H , m W and m G mapping functions, respectively.In our software, the Global Mapping Function [44] is used for this mapping.The a priori ZTD is represented by the ZHD in (1) which is calculated from the Global Pressure and Temperature model, GPT [45] or from an augmentation tropospheric model exploiting actual NWM analysis [41].In (1), the ZWD represents then a correction of the estimated ZTD parameter.The ZTD and horizontal tropospheric gradients are estimated every epoch with applying the random-walk process.All these epoch-wise model parameters are recalculated during the backward smoothing.
Figure 9 demonstrates the behavior of the ZTD from different solutions during the fast change in the troposphere, indicated with a sudden decrease of ZTD by 6-7 cm during a 2.5-h interval at POTS station on 10 May 2013.Precise IGS orbit and clocks are used for the demonstration.Obviously, the post-processing 15-min ZTDs from the GFZ solution using the EPOS software (gray points) and the daily smoothed 30-s stochastic ZTDs from our software (black dots) are in a very good agreement.Note that ZTDs from our solution are resampled to 5 min in the figure.Due to the use of past observations only, the real-time Kalman filter (red points) shows a delay in the change of estimated parameter.The random walk for ZTD was set to 5 mm/sqrt(hour), a commonly used value for any change in the atmosphere.The figure reveals a characteristic behavior for the real-time processing when using past observations in stochastic parameter estimation.We can observe significant improvements mainly in reducing the systematic behavior for the 1-h backward smoothing (green points).The main improvement can be reached within 15-30 min, as indicated for the ZTD from the backward smoothing at 12:00, 13:00 and 14:00.Interestingly, the ZTD from the 2-h smoothing already shows a very good agreement with the post-processing solution.It shows a potential improvement discussed in the previous section and postponing the NRT smoothing by at least 30 min if RT observations are available.
observations can further improve the accuracy of NRT product, especially last parameters of the NRT product and still reduce systematic errors typical for data interval boundaries.

Estimating High-Resolution ZTDs and Horizontal Gradients
The extended model for the GNSS slant tropospheric delay (δT) for a single receiver and a satellite is defined as a function of Zenith Hydrostatic Delay (ZHD), Zenith Wet Delay (ZWD) and horizontal tropospheric gradients pointing to North (GN) and East (GE) using an azimuth (A), which direction is counted from the North [43] ( ) ( ) The ZHD and the ZWD are hydrostatic and wet contributions to the ZTD which represents a vertical delay due to the neutral atmosphere in a profile above the station.Note that ZHD, ZWD and horizontal gradients are projected to any satellite direction by using H m , W m and G m mapping functions, respectively.In our software, the Global Mapping Function [44] is used for this mapping.The a priori ZTD is represented by the ZHD in (1) which is calculated from the Global Pressure and Temperature model, GPT [45] or from an augmentation tropospheric model exploiting actual NWM analysis [41].In (1), the ZWD represents then a correction of the estimated ZTD parameter.The ZTD and horizontal tropospheric gradients are estimated every epoch with applying the random-walk process.All these epoch-wise model parameters are recalculated during the backward smoothing.Figure 9 demonstrates the behavior of the ZTD from different solutions during the fast change in the troposphere, indicated with a sudden decrease of ZTD by 6-7 cm during a 2.5-h interval at POTS station on 10 May 2013.Precise IGS orbit and clocks are used for the demonstration.Obviously, the post-processing 15-min ZTDs from the GFZ solution using the EPOS software (gray points) and the daily smoothed 30-s stochastic ZTDs from our software (black dots) are in a very good agreement.Note that ZTDs from our solution are resampled to 5 min in the figure.Due to the use of past observations only, the real-time Kalman filter (red points) shows a delay in the change of estimated parameter.The random walk for ZTD was set to 5 mm/sqrt(hour), a commonly used value for any change in the atmosphere.The figure reveals a characteristic behavior for the real-time processing when using past observations in stochastic parameter estimation.We can observe significant improvements mainly in reducing the systematic behavior for the 1-h backward smoothing (green points).The main improvement can be reached within 15-30 min, as indicated for the ZTD from the backward smoothing at 12:00, 13:00 and 14:00.Interestingly, the ZTD from the 2-h smoothing already shows a very good agreement with the post-processing solution.It shows a potential improvement discussed in the previous section and postponing the NRT smoothing by at least 30 min if RT observations are available.

Retrieving Slant Tropospheric Delays from Both RT and NRT Processing
The G-Nut/Tefnut supports SINEX_TRO format Version 2 [46] which was newly designed for the exchange of all tropospheric parameters including slant tropospheric delays (SD), meteorological parameters or other auxiliary parameters for reconstructing all components of the slant tropospheric delay model for a single epoch, station and satellite.Using the PPP method, the SDs are retrieved by adding undifferenced carrier-phase post-fit residuals (RES) and, optimally, by subtracting any systematic errors (MPT) such as carrier-phase multipath and errors in antenna phase center variations [47] SD = δT + RES − MPT The estimation of systematic errors, e.g., by generating the so-called multipath maps [48], is possible only through the stacking of carrier-phase post-fit residuals from multiple days for a single station.The post-fit residual stacking hasn't been completed yet in our software, however, we shortly discuss a procedure planned for the continuous production of all tropospheric parameters in both RT and NRT.Due to the temporal stability of the multipath maps, and a necessary presence of time series of post-fit residuals over two weeks at least, it is practical to separate the procedure of post-fit residuals processing from the actual RT and NRT data processing.
For this purpose, we are developing a dedicated tool that will be regularly triggered over archived NRT products for updating multipath maps for all available stations and GNSS constellations.The results will be stored in the SINEX_TRO file and regularly ingested to the new RT/NRT operational processing.The retrievals of the raw slant tropospheric delays consisted of the model slant delay and raw carrier-phase residuals, i.e., the first two terms from (2), could be provided immediately from the RT/NRT analysis.The systematic term in (2) will be introduced later for the stations and satellites/constellations supported from independently provided multipath maps.All parameters, corrections, mapping factors and other auxiliary information such as satellite number, elevation and azimuth angles, are provided in the SINEX_TRO format too, thus a user can decide how to use the product optimally.
In order to demonstrate the impact of the backward smoothing on the slant tropospheric delay retrievals, we exploited a dual GNSS station (POTM, POTS) from the benchmark campaign.We compared differences in slant delays from these two stations located only 2.5 m from each other when considering three retrievals strategies [47]: (1) model slant delays estimated from model parameters, δT (model SD), (2) raw slant delays, δT + RES (raw SD), and (3) clean slant delays, δT + RES − MPT (clean SD).The estimated slant delays at both stations should theoretically represent the same tropospheric effects for each particular satellite, while non-zero differences reflect mainly the uncertainty from the data processing or station-specific systematic errors.
Figure 10 shows the elevation dependence of slant delay biases, standard deviations and the normalized statistics calculated from the differences divided by actual slant delay; all provided for the period of May-June 2013.The three variants of slant delay retrievals are plotted by solid lines (model SDs), dotted lines (raw SDs), and dashed lines (clean SDs).Additionally, different colors of solutions depict the utilization of real-time analyses (F-Kalman filter) or post-processing (S-backward smoothing) and the use of IGS precise products (PP) or real-time products (RT).Generally, performances of biases for the model and the clean SDs are in a very good agreement and normalized biases do not vary with the elevation.On the other hand, a clear elevation-specific pattern is visible for the normalized biases of raw slant delays indicating systematic effects at one of the stations [47].There are negligible differences in terms of variants using the solution (F vs. S) or input products (PP vs. RT).Standard deviations show however different behavior for all three variants of slant delays.As expected, the model SDs are in the best agreement, but it should be stressed that they do not contain full information about the atmosphere, but only that approximated by the ZTD and horizontal gradients.The clean SD performs much better compared to the raw SDs; both contain additional spatio-temporal information about the water vapor content in the troposphere [47].
The most interesting is the impact of the backward smoothing which shows the largest improvement for the model SD compared to the clean SDs and the raw SDs.A similar impact is visible for the solutions using precise and real-time products.The backward smoothing has a dominating effect on improving the estimated parameters, which has a positive impact on slant delays too.Finally, Figure 11 shows similar statistics for the GFZ post-processing solution (GFZ) and our solution (GOP_S), both previously evaluated in [47].The new post-processing variants from the above figure are also included for the comparison purposes.Interestingly, we have significantly improved the original GOP_S solution mainly by resolving specific problems related to the backward smoothing combined with the SD retrievals from the post-fit carrier-phase residuals.Thanks to the troposphere stochastic modeling, the new solution now shows the best performance, even when compared to the deterministic model and the LSQ with a lower-parameter sampling rate.

Assessment of New Method Compared to the Existing E-GVAP Processing
The new strategy has been initially developed and assessed using the GNSS4SWEC Benchmark dataset [31] and firstly compared to the GOP NRT tropospheric solution contributing operationally to E-GVAP.Though the new strategy can provide RT and NRT products in high temporal resolution, we compared only the ZTD as a product of HH:00 and HH:59 time stamps in every hour representing the standard NRT E-GVAP product, see Figure 7.
Table 4 summarizes results of three strategies and six ZTD solutions using 13 EUREF stations selected from the benchmark campaign and exploiting the EUREF combined ZTD product as a reference for all comparisons [39].The table show summary statistics indicating similar improvements in terms of the standard deviation over all ZTDs estimated at HH:00 in NRT independently of applied products (IGS RTS vs. IGS final products), processing strategies and software (G-Nut/Tefnut PPP vs. Bernese V52 DD).Compared to the Kalman filter ZTD estimated at the last epoch (HH:59), the backward smoothing running on hourly basis showed the improvement of 20% and 24% for the IGS RTS and the IGS final product, respectively.The E-GVAP/GOP product demonstrates a similar improvement (24%) comparing ZTD from HH:00 against HH:59, which corresponds to our previous results [29].The new adaptable PPP solution using the IGS final orbits Finally, Figure 11 shows similar statistics for the GFZ post-processing solution (GFZ) and our solution (GOP_S), both previously evaluated in [47].The new post-processing variants from the above figure are also included for the comparison purposes.Interestingly, we have significantly improved the original GOP_S solution mainly by resolving specific problems related to the backward smoothing combined with the SD retrievals from the post-fit carrier-phase residuals.Thanks to the troposphere stochastic modeling, the new solution now shows the best performance, even when compared to the deterministic model and the LSQ with a lower-parameter sampling rate.Finally, Figure 11 shows similar statistics for the GFZ post-processing solution (GFZ) and our solution (GOP_S), both previously evaluated in [47].The new post-processing variants from the above figure are also included for the comparison purposes.Interestingly, we have significantly improved the original GOP_S solution mainly by resolving specific problems related to the backward smoothing combined with the SD retrievals from the post-fit carrier-phase residuals.Thanks to the troposphere stochastic modeling, the new solution now shows the best performance, even when compared to the deterministic model and the LSQ with a lower-parameter sampling rate.

Assessment of New Method Compared to the Existing E-GVAP Processing
The new strategy has been initially developed and assessed using the GNSS4SWEC Benchmark dataset [31] and firstly compared to the GOP NRT tropospheric solution contributing operationally to E-GVAP.Though the new strategy can provide RT and NRT products in high temporal resolution, we compared only the ZTD as a product of HH:00 and HH:59 time stamps in every hour representing the standard NRT E-GVAP product, see Figure 7.
Table 4 summarizes results of three strategies and six ZTD solutions using 13 EUREF stations selected from the benchmark campaign and exploiting the EUREF combined ZTD product as a reference for all comparisons [39].The table show summary statistics indicating similar improvements in terms of the standard deviation over all ZTDs estimated at HH:00 in NRT independently of applied products (IGS RTS vs. IGS final products), processing strategies and software (G-Nut/Tefnut PPP vs. Bernese V52 DD).Compared to the Kalman filter ZTD estimated at the last epoch (HH:59), the backward smoothing running on hourly basis showed the improvement of 20% and 24% for the IGS RTS and the IGS final product, respectively.The E-GVAP/GOP product demonstrates a similar improvement (24%) comparing ZTD from HH:00 against HH:59, which corresponds to our previous results [29].The new adaptable PPP solution using the IGS final orbits

Assessment of New Method Compared to the Existing E-GVAP Processing
The new strategy has been initially developed and assessed using the GNSS4SWEC Benchmark dataset [31] and firstly compared to the GOP NRT tropospheric solution contributing operationally to E-GVAP.Though the new strategy can provide RT and NRT products in high temporal resolution, we compared only the ZTD as a product of HH:00 and HH:59 time stamps in every hour representing the standard NRT E-GVAP product, see Figure 7.
Table 4 summarizes results of three strategies and six ZTD solutions using 13 EUREF stations selected from the benchmark campaign and exploiting the EUREF combined ZTD product as a reference for all comparisons [39].The table show summary statistics indicating similar improvements in terms of the standard deviation over all ZTDs estimated at HH:00 in NRT independently of applied products (IGS RTS vs. IGS final products), processing strategies and software (G-Nut/Tefnut PPP vs. Bernese V52 DD).Compared to the Kalman filter ZTD estimated at the last epoch (HH:59), the backward smoothing running on hourly basis showed the improvement of 20% and 24% for the IGS RTS and the IGS final product, respectively.The E-GVAP/GOP product demonstrates a similar improvement (24%) comparing ZTD from HH:00 against HH:59, which corresponds to our previous results [29].The new adaptable PPP solution using the IGS final orbits reached the same accuracy as the E-GVAP/GOP product using the IGS ultra-rapid orbits and NRT DD network solution from the Bernese GNSS Software.On the other hand, the use of IGS RTS products instead of IGS final products in PPP indicates a degradation of 18% in ZTD SDEV and a 2.5 mm bias, which agrees with the summary in Sections 2.2 and 2.3.It should be finally noted, that the E-GVAP/GOP solution and the reference EUREF solution are based on a similar processing strategy and the software, while the new strategy is significantly different.Figure 12 shows standard deviations and biases individually for all stations comparing the first ZTDs (HH:00) and the last ZTDs (HH:59).The statistics of ZTDs from the E-GVAP/GOP solution (PPP:DD-ultra) are plotted in red and pink for HH:00 and HH:59, respectively.The results from the new strategy using the PPP with IGS final orbit and clock products (PPP:IGS_final) are shown in dark and light blue and, using IGS RTS (PPP:IGS03) in black and grey.Standard deviations for all the stations show a similar improvement in the ZTD SDEV over all the strategies, software and precise products.However, there is no significant impact of the strategy on systematic errors, and we can observe only a common positive bias attributed to the use of the IGS RTS products.
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 22 reached the same accuracy as the E-GVAP/GOP product using the IGS ultra-rapid orbits and NRT DD network solution from the Bernese GNSS Software.On the other hand, the use of IGS RTS products instead of IGS final products in PPP indicates a degradation of 18% in ZTD SDEV and a 2.5 mm bias, which agrees with the summary in Sections 2.2 and 2.3.It should be finally noted, that the E-GVAP/GOP solution and the reference EUREF solution are based on a similar processing strategy and the software, while the new strategy is significantly different.Figure 12 shows standard deviations and biases individually for all stations comparing the first ZTDs (HH:00) and the last ZTDs (HH:59).The statistics of ZTDs from the E-GVAP/GOP solution (PPP:DD-ultra) are plotted in red and pink for HH:00 and HH:59, respectively.The results from the new strategy using the PPP with IGS final orbit and clock products (PPP:IGS_final) are shown in dark and light blue and, using IGS RTS (PPP:IGS03) in black and grey.Standard deviations for all the stations show a similar improvement in the ZTD SDEV over all the strategies, software and precise products.However, there is no significant impact of the strategy on systematic errors, and we can observe only a common positive bias attributed to the use of the IGS RTS products.

Impact Assessments on Estimated Parameters at Collocated Stations
The impact of various aspects of the new strategy on all the tropospheric parameters can be optimally assessed using closely collocated GNSS stations, e.g., within few meters.Although different instrumentation-specific effects, such as phase center modeling and the quality of a receiver tracking, can affect analyses at both stations, the station should principally observe the same tropospheric delays.For the purposes of our evaluation, we selected two IGS station pairs, ZIM2-ZIMJ (Zimmerwald, Switzerland) and MAT1-MATE (Matera, Italy), all collecting data from GPS, GLONASS and Galileo systems within 10 m and 2.5 m in horizontal and vertical distances,

Impact Assessments on Estimated Parameters at Collocated Stations
The impact of various aspects of the new strategy on all the tropospheric parameters can be optimally assessed using closely collocated GNSS stations, e.g., within few meters.Although different instrumentation-specific effects, such as phase center modeling and the quality of a receiver tracking, can affect analyses at both stations, the station should principally observe the same tropospheric delays.For the purposes of our evaluation, we selected two IGS station pairs, ZIM2-ZIMJ (Zimmerwald, Switzerland) and MAT1-MATE (Matera, Italy), all collecting data from GPS, GLONASS and Galileo systems within 10 m and 2.5 m in horizontal and vertical distances, respectively.In the following sections, we evaluate impacts of multi-constellation analyses and various input products for PPP during September 2017.We used the GFZ MGEX (GBM) product [49] and CODE MGEX (COM) product [50] as multi-GNSS final products and, additionally, we compared them with the utilization the IGS RTS (IGS03) product for the operation in real time.

Impact of Backward Smoothing on Adaptable RT/NRT Solutions
The new strategy was evaluated using multi-constellation GBM and COM products.Dual-station differences were calculated for the ZTD and horizontal linear gradients when considering various delays (5-60 min) applied for the backward smoothing.In other words, the dependency on delay characterizes a possible improvement due to the increased latency of the product.Figure 13 shows the statistics for both collocated stations using the standalone GPS, and multi-constellation (GPS + GLO and GPS + GLO + GAL) solutions.The positive impact is observed when comparing the backward smoothing and the Kalman filter, and it increases steadily for both the ZTD and the east horizontal gradient (the north gradient is not shown here, but similar to the east gradient) along with the delay for triggering the backward smoothing.The impact becomes significant after the 20-min delay and is larger for the collocated stations at Matera.A very similar impact is observed for multi-constellations and GPS-only solutions and resulting in smaller discrepancies of the tropospheric parameters at both collocated stations.
Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 22 various input products for PPP during September 2017.We used the GFZ MGEX (GBM) product [49] and CODE MGEX (COM) product [50] as multi-GNSS final products and, additionally, we compared them with the utilization the IGS RTS (IGS03) product for the operation in real time.

Impact of Backward Smoothing on Adaptable RT/NRT Solutions
The new strategy was evaluated using multi-constellation GBM and COM products.Dual-station differences were calculated for the ZTD and horizontal linear gradients when considering various delays (5-60 min) applied for the backward smoothing.In other words, the dependency on delay characterizes a possible improvement due to the increased latency of the product.Figure 13 shows the statistics for both collocated stations using the standalone GPS, and multi-constellation (GPS + GLO and GPS + GLO + GAL) solutions.The positive impact is observed when comparing the backward smoothing and the Kalman filter, and it increases steadily for both the ZTD and the east horizontal gradient (the north gradient is not shown here, but similar to the east gradient) along with the delay for triggering the backward smoothing.The impact becomes significant after the 20-min delay and is larger for the collocated stations at Matera.A very similar impact is observed for multi-constellations and GPS-only solutions and resulting in smaller discrepancies of the tropospheric parameters at both collocated stations.

Tropospheric Parameters from Multi-GNSS Analyses
Improved solutions using multi-constellations and the backward smoothing are demonstrated in the time series of ZTD and horizontal gradient differences obtained from the two collocated stations, ZIM2 and ZIMJ, Figure 14.Results of the single-/multi-constellations are visualized by different colors: (1) standalone GPS in red, (2) GPS + GLO in green, and (3) GPS + GLO + GAL in blue.A positive effect is visible for all parameters, and is similar using both the Kalman filter (a) and the backward smoothing (b).The scatters for multi-constellations are smaller compared to the standalone GPS solution.An even more significant effect of the smoothing is visible for the gradient parameters.Theoretically, zero differences are expected for the collocated stations with the same antenna height.However, a vertical difference between ZIM2 and ZIMJ is about 2 m which can cause about 0.5 mm difference in ZTDs when considering the pressure decreases approximately by 11.3 Pa per meter near the geoid and the 100 Pa difference in the atmospheric pressure causes a 2.27 mm difference in ZHD [51].As we observe a ZTD difference about 2-3 mm, it can be still attributed to remaining station-specific systematic errors, e.g., such as phase center offset and variation models.

Tropospheric Parameters from Multi-GNSS Analyses
Improved solutions using multi-constellations and the backward smoothing are demonstrated in the time series of ZTD and horizontal gradient differences obtained from the two collocated stations, ZIM2 and ZIMJ, Figure 14.Results of the single-/multi-constellations are visualized by different colors: (1) standalone GPS in red, (2) GPS + GLO in green, and (3) GPS + GLO + GAL in blue.A positive effect is visible for all parameters, and is similar using both the Kalman filter (a) and the backward smoothing (b).The scatters for multi-constellations are smaller compared to the standalone GPS solution.An even more significant effect of the smoothing is visible for the gradient parameters.Theoretically, zero differences are expected for the collocated stations with the same antenna height.However, a vertical difference between ZIM2 and ZIMJ is about 2 m which can cause about 0.5 mm difference in ZTDs when considering the pressure decreases approximately by 11.3 Pa per meter near the geoid and the 100 Pa difference in the atmospheric pressure causes a 2.27 mm difference in ZHD [51].As we observe a ZTD difference about 2-3 mm, it can be still attributed to remaining station-specific systematic errors, e.g., such as phase center offset and variation models.Numerical statistics (biases and standard deviations) characterize an impact of single-and multi-constellation solutions on the estimated ZTDs and horizontal tropospheric gradients when using the Kalman filter, Table 5.It should be noted that GLONASS (R) and Galileo (E) observations were down-weighted by a factor of 2 with respect GPS (G) to reflect less accurate precise products and models.A positive effect of multi-constellation is visible at all parameters, and particularly in terms of the standard deviation, while the impact of GLONASS is more significant compared to Galileo.That is expected, as the GLONASS is the operational system longer supported by precise models and products in the scientific community as well as due to a lower number of operational Galileo satellites.As already discussed, ZIM2-ZIMJ differences indicate a bias of about 2-3 mm in ZTD.The improvements in all parameters reached 15-30% in terms of RMSE at both dual-stations.Table 6 then shows the impact of the backward smoothing on all solutions using single-or multi-constellation data.All the above-mentioned characteristics are similar to the Kalman filter, and the backward smoothing then improved mainly standard deviations (by about 25%).Numerical statistics (biases and standard deviations) characterize an impact of single-and multi-constellation solutions on the estimated ZTDs and horizontal tropospheric gradients when using the Kalman filter, Table 5.It should be noted that GLONASS (R) and Galileo (E) observations were down-weighted by a factor of 2 with respect GPS (G) to reflect less accurate precise products and models.A positive effect of multi-constellation is visible at all parameters, and particularly in terms of the standard deviation, while the impact of GLONASS is more significant compared to Galileo.That is expected, as the GLONASS is the operational system longer supported by precise models and products in the scientific community as well as due to a lower number of operational Galileo satellites.As already discussed, ZIM2-ZIMJ differences indicate a bias of about 2-3 mm in ZTD.The improvements in all parameters reached 15-30% in terms of RMSE at both dual-stations.Table 6 then shows the impact of the backward smoothing on all solutions using single-or multi-constellation data.All the above-mentioned characteristics are similar to the Kalman filter, and the backward smoothing then improved mainly standard deviations (by about 25%).

Impact of Precise Products on ZTD and Gradient Estimates
Table 7 summarizes the ZTD and horizontal gradient statistics indicating the impact of different orbit and clock products and, again, the backward smoothing, all using the same processing strategy and the software.The first value in each table column represents the result of the Kalman filter processing and the second value represents the result of the backward smoothing.Comparison of the GPS-only solutions to the COM and GBM MGEX products are included too, both indicating that the usage of IGS RTS products might still degrade the accuracy up to a factor of 2. A significant difference is observed also in statistics of IGS RT products (IGS01 and IGS03).The IGS03 solution performed better than IGS01 during September 2017, which seems to be in contradiction to the results of a long-term evaluation in Section 2.1.Indeed, the precision (SDEV) of the IGS RTS clocks was 3.5 cm and 3.0 cm during this month for IGS01 and IGS03, respectively, indicating a high variability in the actual performance of the real-time products.

Carrier-Phase Post-Fit Residuals and Slant Delays
Figure 15 shows the carrier-phase post-fit residuals when using the Kalman filter PPP (a) and the backward smoothing PPP (b) for multi-GNSS solutions supported with the COM (top) and GBM (bottom) MGEX products.The carrier-phase residuals are useful indicators of an overall performance of the solution including the quality of input products and models.Showing plots for the ZIM2 station only, below discussed characteristics are common to other stations too.First, we observe a common elevation-dependent pattern of characteristics of post-fit residuals when using elevation-dependent observation weighting, 1/ sin 2 (e).Second, the backward smoothing does not change the distribution of the carrier-phase post-fit residuals significantly.The main effect of the backward smoothing is thus understood mainly as improved accuracy of the estimated parameters.The tropospheric slant delays reconstructed from the model parameters and post-fit residuals will thus benefit primarily from the improvement of the parameters.Third, the GPS residuals (black) are the smallest and compact compared to other systems indicating actual quality of precise models and products.Galileo shows the largest residuals, however, we had to substitute various precise models, in particular station antenna phase center offsets and variations by using the values from GPS models.Due to the same reason, we may notice systematic changes in the elevation-dependent redistribution of Galileo residuals (red) after applying the backward smoothing.Fourth, we notice post-fit residuals GLONASS (green) about double the size when using the GBM product compared to the COM product.As the characteristics are common to all the stations, it indicates a lower quality of GLONASS orbits and clocks from the GBM product or some inconsistent models used for the product generation and in the PPP software.
Remote Sens. 2018, 10, x FOR PEER REVIEW 18 of 22 models.Due to the same reason, we may notice systematic changes in the elevation-dependent redistribution of Galileo (red) after applying the backward smoothing.Fourth, we notice post-fit residuals GLONASS (green) about double the size when using the GBM product compared to the COM product.As the characteristics are common to all the stations, it indicates a lower quality of GLONASS orbits and clocks from the GBM product or some inconsistent models used for the product generation and in the PPP software.

Conclusions
We described a new strategy for generating advanced RT and NRT tropospheric products exploiting the PPP method and, simultaneously, the Kalman filter and the backward smoothing supported with the IGS real-time products.The strategy can be used to provide RT and NRT products synchronously, while it can be further optimized either for the latency or the accuracy of the NRT product.Although the NRT solution is generated as a side product of the real-time analysis, in terms of the precision, it is comparable to the traditional NRT ZTDs using the LSQ adjustment and batches of double-difference observations in a network mode, still commonly used within the E-GVAP.In addition, both the RT and NRT products can produce a consistent set of all parameters, namely the ZTD, horizontal tropospheric gradients and slant delays, all provided in high resolution and using optimally all observations from multi-constellations if precise products are available.
A long-term assessment of the IGS RTS in terms of the quality demonstrated that the products were available over 90% during the period 2013-2017 with 2-3 longer gaps only, and reached the 3D orbit accuracy of 4-6 cm and the clock precision of 2-4 cm.All the products were stable and were usable for the purpose of the troposphere monitoring, though some temporal variability was observed in the quality.An assessment of 1-year ZTDs generated routinely in the real-time PPP solution reached a stable precision of 6-8 mm over 18 global stations with about 10% improvements when additionally including GLONASS observations.A simulated analysis of the impact of the backward smoothing on the ZTD and horizontal tropospheric gradients in NRT were evaluated within the GNSS4SWEC Benchmark campaign.Improvements of 20% for NRT ZTDs compared to RT ZTDs were demonstrated using the new strategy, i.e., the same improvement as observed within the piece-wise linear NRT ZTD estimation using the LSQ data processing.The overall improvements of the backward smoothing algorithm reached over 20%, independently if the final or real-time precise products are used.The impact of RT precise products on ZTD estimates indicates a systematic error of 2.4-2.8mm and a degradation of

Conclusions
We described a new strategy for generating advanced RT and NRT tropospheric products exploiting the PPP method and, simultaneously, the Kalman filter and the backward smoothing supported with the IGS real-time products.The strategy can be used to provide RT and NRT products synchronously, while it can be further optimized either for the latency or the accuracy of the NRT product.Although the NRT solution is generated as a side product of the real-time analysis, in terms of the precision, it is comparable to the traditional NRT ZTDs using the LSQ adjustment and batches of double-difference observations in a network mode, still commonly used within the E-GVAP.In addition, both the RT and NRT products can produce a consistent set of all parameters, namely the ZTD, horizontal tropospheric gradients and slant delays, all provided in high resolution and using optimally all observations from multi-constellations if precise products are available.
A long-term assessment of the IGS RTS in terms of the quality demonstrated that the products were available over 90% during the period 2013-2017 with 2-3 longer gaps only, and reached the 3D orbit accuracy of 4-6 cm and the clock precision of 2-4 cm.All the products were stable and were usable for the purpose of the troposphere monitoring, though some temporal variability was observed in the quality.An assessment of 1-year ZTDs generated routinely in the real-time PPP solution reached a stable precision of 6-8 mm over 18 global stations with about 10% improvements when additionally including GLONASS observations.A simulated analysis of the impact of the backward smoothing on the ZTD and horizontal tropospheric gradients in NRT were evaluated within the GNSS4SWEC Benchmark campaign.Improvements of 20% for NRT ZTDs compared to RT ZTDs were demonstrated using the new strategy, i.e., the same improvement as observed within the piece-wise linear NRT ZTD estimation using the LSQ data processing.The overall improvements of the backward smoothing algorithm reached over 20%, independently if the final or real-time precise products are used.The impact of RT precise products on ZTD estimates indicates a systematic error of 2.4-2.8mm and a degradation of 13-17% when compared with the use of final products, however, up to a factor of 2 worse in the accuracy when comparing parameters from collocated stations.
We further evaluated differences of ZTDs and horizontal tropospheric gradients from two collocated GNSS stations.We observed a significant impact of the backward smoothing on estimated tropospheric parameters when applied with 30-min latency in near-real time.A similar effect was observed for single-and multi-constellation solutions.An improvement of about 25% was then achieved when using multi-constellation compared to the standalone GPS, though observations from GLONASS and Galileo systems were down-weighted and their precise products and models are less accurate.Models were verified by visualizing carrier-phase post-fit residuals at individual stations.It showed that the backward smoothing improves mainly adjusted parameters, but does not affect essentially the distribution of post-fit residuals, however, retrieved slant delays still benefit from improved estimated parameters.

Figure 1 .
Figure 1.Monthly statistics of availability of satellite corrections from IGS01 RT stream.

Figure 2 .
Figure 2. Monthly statistics of availability of satellite corrections from IGS02 RT stream.

Figure 3 .
Figure 3. Monthly statistics of availability of satellite corrections from CNS91 RT stream.

Figure 1 .
Figure 1.Monthly statistics of availability of satellite corrections from IGS01 RT stream.

Figure 1 .
Figure 1.Monthly statistics of availability of satellite corrections from IGS01 RT stream.

Figure 2 .
Figure 2. Monthly statistics of availability of satellite corrections from IGS02 RT stream.

Figure 3 .
Figure 3. Monthly statistics of availability of satellite corrections from CNS91 RT stream.

Figure 2 .
Figure 2. Monthly statistics of availability of satellite corrections from IGS02 RT stream.

Figure 1 .
Figure 1.Monthly statistics of availability of satellite corrections from IGS01 RT stream.

Figure 2 .
Figure 2. Monthly statistics of availability of satellite corrections from IGS02 RT stream.

Figure 3 .
Figure 3. Monthly statistics of availability of satellite corrections from CNS91 RT stream.Figure 3. Monthly statistics of availability of satellite corrections from CNS91 RT stream.

Figure 3 .
Figure 3. Monthly statistics of availability of satellite corrections from CNS91 RT stream.Figure 3. Monthly statistics of availability of satellite corrections from CNS91 RT stream.

Figure 4 .
Figure 4. Daily RMS of real-time orbits with respect to IGS final orbit.

Figure 4 .
Figure 4. Daily RMS of real-time orbits with respect to IGS final orbit.

Figure 5
Figure 5 finally shows time series of the clock SDEV and RMSE statistics.The former (top plot) indicates a comparable high quality over the period for IGS01 and CNS91, while more outliers are observed for IGS02 and IGS03 including the problematic period in 2015 identified in the orbit availability evaluation.The clock RMSE from IGS01 is the lowest and the most stable compared to the others during the period while, the RMSE of CNS91 clocks was more accurate during 2015 when compared to the other years.Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 22

Figure 6 .
Figure 6.Monthly summary biases and standard deviations of real-time ZTDs over 18 stations.

Figure 6 .
Figure 6.Monthly summary biases and standard deviations of real-time ZTDs over 18 stations.

Figure 7 .
Figure 7. Tropospheric parameters in NRT products using the deterministic model (dotted lines).

Figure 7 .
Figure 7. Tropospheric parameters in NRT products using the deterministic model (dotted lines).

Figure 8 .
Figure 8. Stochastic modeling of the troposphere represented by real-time Kalman filter (white points), hourly backward smoothing (blue points), and 45 min postponed hourly backward smoothing (red points).

Figure 8 .
Figure 8. Stochastic modeling of the troposphere represented by real-time Kalman filter (white points), hourly backward smoothing (blue points), and 45 min postponed hourly backward smoothing (red points).

Figure 9 .
Figure 9. Real time (red), near-real time (green, blue) and reference (gray, black) ZTD estimates during a fast change in the troposphere.

Figure 9 .
Figure 9. Real time (red), near-real time (green, blue) and reference (gray, black) ZTD estimates during a fast change in the troposphere.

Figure 11 .
Figure 11.Elevation-dependent absolute (a) and normalized (b) statistics for GOP slant delay differences at POTM and POTS stations using post-processing (PP) product and Kalman filter (F) or backward smoothing (S), GFZ slant delays and old GOP product (GOP_S).Line types indicates model SD (solid), clean SD (dashed), raw SD (dotted).

Figure 11 .
Figure 11.Elevation-dependent absolute (a) and normalized (b) statistics for GOP slant delay differences at POTM and POTS stations using post-processing (PP) product and Kalman filter (F) or backward smoothing (S), GFZ slant delays and old GOP product (GOP_S).Line types indicates model SD (solid), clean SD (dashed), raw SD (dotted).

Figure 11 .
Figure 11.Elevation-dependent absolute (a) and normalized (b) statistics for GOP slant delay differences at POTM and POTS stations using post-processing (PP) product and Kalman filter (F) or backward smoothing (S), GFZ slant delays and old GOP product (GOP_S).Line types indicates model SD (solid), clean SD (dashed), raw SD (dotted).

Figure 12 .
Figure 12.ZTD standard deviations and biases for 13 EUREF stations and six processing strategies compared to the reference EUREF product.

Figure 12 .
Figure 12.ZTD standard deviations and biases for 13 EUREF stations and six processing strategies compared to the reference EUREF product.

Figure 13 .
Figure 13.Standard deviations for ZTDs (a) and east gradients (b) using different delays in backward smoothing at two dual-stations-Zimmerwald (top) and Matera (bottom).

Figure 13 .
Figure 13.Standard deviations for ZTDs (a) and east gradients (b) using different delays in backward smoothing at two dual-stations-Zimmerwald (top) and Matera (bottom).

Figure 14 .
Figure 14.Time series of ZTD (top), north gradient (middle) and east gradient (bottom) differences at Zimmerwald dual-station when using the Kalman filter (a) and the backward smoothing (b).

Figure 14 .
Figure 14.Time series of ZTD (top), north gradient (middle) and east gradient (bottom) differences at Zimmerwald dual-station when using the Kalman filter (a) and the backward smoothing (b).

Figure 15 .
Figure 15.Carrier-phase post-fit residuals from the Kalman filter (a) and the backward smoothing (b) using COM (top) and GBM (bottom) products at ZIM2 station.

Figure 15 .
Figure 15.Carrier-phase post-fit residuals from the Kalman filter (a) and the backward smoothing (b) using COM (top) and GBM (bottom) products at ZIM2 station.

Table 1 .
RT clocks and orbits components compared to the IGS final products (units: cm).

Table 1 .
RT clocks and orbits components compared to the IGS final products (units: cm).

Table 1
also summarizes RMSE and standard deviation (SDEV) of the real-time clock corrections.

Table 2 .
Summary statistics from the comparison of PPP ZTD results using two inputs (IGS01 RT vs. IGS final)

Table 3 .
Summary statistics over 18 stations from routine RT product using GPS and GPS + GLO data w.r.t.EUREF reprocessing.

Table 2 .
Summary statistics from the comparison of PPP ZTD results using two inputs (IGS01 RT vs.

Table 3 .
Summary statistics over 18 stations from routine RT product using GPS and GPS + GLO data w.r.t.EUREF reprocessing.

Table 4 .
Summary statistics of three processing strategies and six ZTD solutions compared to EUREF combined tropospheric product.

Table 4 .
Summary statistics of three processing strategies and six ZTD solutions compared to EUREF combined tropospheric product.

Table 7 .
Summary statistics for dual stations using different products: real-time (IGS01, IGS03) and final MGEX (COM, GBM), the Kalman filter and, in addition, the backward smoothing.