Analysing the Zenith Tropospheric Delay Estimates in On-line Precise Point Positioning (PPP) Services and PPP Software Packages

As Global Navigation Satellite System (GNSS) signals travel through the troposphere, a tropospheric delay occurs due to a change in the refractive index of the medium. The Precise Point Positioning (PPP) technique can achieve centimeter/millimeter positioning accuracy with only one GNSS receiver. The Zenith Tropospheric Delay (ZTD) is estimated alongside with the position unknowns in PPP. Estimated ZTD can be very useful for meteorological applications, an example is the estimation of water vapor content in the atmosphere from the estimated ZTD. PPP is implemented with different algorithms and models in online services and software packages. In this study, a performance assessment with analysis of ZTD estimates from three PPP online services and three software packages is presented. The main contribution of this paper is to show the accuracy of ZTD estimation achievable in PPP. The analysis also provides the GNSS users and researchers the insight of the processing algorithm dependence and impact on PPP ZTD estimation. Observation data of eight whole days from a total of nine International GNSS Service (IGS) tracking stations spread in the northern hemisphere, the equatorial region and the southern hemisphere is used in this analysis. The PPP ZTD estimates are compared with the ZTD obtained from the IGS tropospheric product of the same days. The estimates of two of the three online PPP services show good agreement (<1 cm) with the IGS ZTD values at the northern and southern hemisphere stations. The results also show that the online PPP services perform better than the selected PPP software packages at all stations.


Introduction
GNSS data is widely used for positioning and navigation in mass-market and engineering applications [1] and for altitude determination [2], moreover, it can also be used for monitoring the atmosphere. The electron content in the ionosphere and the air density in the electrically neutral atmosphere (troposphere) affect GNSS signals propagating through the atmosphere. The influence of the troposphere is described by the total refractivity N, which depends on pressure, temperature and water vapor partial pressure [3]. An example of the use of GNSS data for applications other than positioning and navigation is the integration of the GNSS-derived Path Delay with microwave radiometer measurements to find a precise wet tropospheric correction for altimetric products [4]. Another use of GNSS data is the remote sensing of the atmosphere, where GNSS signals can be used packages. The days chosen for this assessment are days 027,118,208 and 300 of year 2016 and days 027,117,207 and 299 of year 2017, covering different meteorological conditions during different seasons. This assessment is relevant for GNSS meteorology applications such as water vapor calculation using ZWD and Numeric Weather Prediction models where a precise tropospheric delay is needed.

Tropospheric Delay
GNSS signals propagating through the atmosphere are delayed due to the free electron content in the ionosphere and by the air density in the electrically neutral atmosphere also called troposphere. The refractive index n or the total refractivity N of the troposphere is described by: The total refractivity of the troposphere can be separated into two main components, the hydrostatic or dry (N dry ) and the wet component (N wet ) caused by dry gases and the water vapor respectively [19] and it can be expressed as a function of meteorological parameters such as air pressure p, temperature T and water vapor partial pressure e [3]: where k 1 = 77.689 K·h·Pa −1 , k 2 = 71.295 K·h·Pa −1 and k 3 = 375,463 K 2 ·h·Pa −1 are empirically determined coefficients [20]. The troposphere causes a delay to the signal ∆ PD which can be expressed as an integral of the total refractivity N along the propagation path s from receiver r to the satellite w: The tropospheric delay can also be separated in the hydrostatic and the wet component. Therefore, Equation (3) can be written as: ∆ PD = 10 −6 w r N dry ds + 10 −6 w r N wet ds (4) The total tropospheric delay in slant path delay can be mapped to the zenith direction, yielding the Zenith Tropospheric Delay (ZTD) using a mapping function depending on the elevation angle of the satellite. The ZTD is defined as the addition of the Zenith Hydrostatic Delay (ZHD) and the Zenith Wet Delay (ZHD) which are the left and the right side of Equation (4): where m h (E) and m w (E) are the hydrostatic and wet mapping functions depending on the elevation angle. The ZTD can be determined as an integral of N in the zenith direction [20]: Equation (7) indicates that the Zenith Tropospheric Delay and the refractivity of the troposphere are related. Since the refractivity depends on meteorological conditions along the signal path, the ZTD can be related to these conditions. The ZTD is a parameter estimated in the Precise Point Positioning technique.

Precise Point Positioning
Undifferenced pseudorange and carrier phase observations from a single GPS receiver are processed in the PPP technique [21]. The technique is named "precise" because precise a priori information such as satellites orbits and clock errors from different sources like the International GNSS Service (IGS) [22], the Jet Propulsion Lab (JPL) or Natural Resources Canada (NRCan) [23] are used in the data processing. The advantages of this technique are the accuracies obtained (cm-level with only one receiver [23]) and that it eliminates the need to acquire simultaneous tracking data from a reference station or from a network of stations [22]. Moreover, the use of a single receiver reduces equipment cost and makes the processing less labor and resources intensive.

Observation Equations
The ionosphere causes a delay to the GNSS signal propagating through the atmosphere which can be reduced greatly with dual-frequency data. The ionosphere-free combinations of dual-frequency GPS pseudoranges (P) and carrier-phase observations (ϕ) are related to user position, clock, troposphere and ambiguity parameters according to the following simplified observation equations [12]: where l p is the ionosphere-free combination of L1 and L2 pseudoranges, l ϕ is the ionosphere-free combination of L1 and L2 carrier-phases, dt is the receiver's clock offset from GPS time, dT is the satellite clock offset from GPS time, C is the speed of light in vacuum, Tr is the signal path delay due to the neutral-atmosphere (primarily the troposphere), λ is the carrier wavelength, B is the ambiguity of the carrier-phase ionosphere-free combination, and ε p and ε ϕ are the measurement noise components.
ρ is the geometrical range computed as a function of satellite (Xs,Ys,Zs) and station coordinates (x,y,z) defined as: the position, the tropospheric delay, the ambiguity and the clock offsets need to be estimated using estimators such as the Extended Kalman Filter (EKF) or the Least Squares technique.

Materials and Methods
The PPP technique has been implemented by different software packages which follow different estimation strategies or use precise ephemeris from different sources. There are online services accessible through the Internet and software packages that have to be installed and run locally in a computer. For this study, six PPP post-processing software packages were used in total, three of them are online services, APPS, CSRS-PPP and MagicGNSS and three stand-alone post-processing software packages gLAB, POINT and RTKLIB. A summary of the characteristics of the online-services and the software used for this study can be found in Table 1.
The Automatic Precise Positioning Service (APPS) is an online service by the Jet Propulsion Lab which can estimate position coordinates as a single set in Static Mode or a time series in Kinematic Mode. In APPS the receiver clock states are estimated as white noise with updates every measurement epoch and the Zenith Wet delay (ZWD) is estimated as a random walk with variance of 3 mm 2 per hour. Moreover, the wet delay gradient is estimated as a random walk with variance of 0.3 mm 2 per hour and the phase ambiguities are estimated as real numbers [24].
APPS can take only dual-frequency GPS observations. The service allows the user to decide whether to use Final, Rapid and Ultra-rapid type products from the Jet Propulsion Lab (JPL) for corrections of obits and clocks of satellites. The ZTD is estimated applying the Global Mapping Function troposphere mapping function with an a priori hydrostatic delay of 1.013 * 2.27 * e −0.000116 * h meters where h is the station height above the ellipsoid in meters and a priori wet delay of 0.1 m. The wet delay is estimated together with positioning unknowns.
The Canadian Spatial Reference System (CSRS-PPP service) by Natural Resources Canada uses a dynamic filter to estimate the station position in static or kinematic mode, the station-clock states, the local tropospheric zenith delays and the carrier-phase ambiguities. The approach used in CSRS-PPP for ZTD estimation is to smooth the estimates by a backward substitution with the final converged satellite ambiguity parameters held fixed for all epochs. This approach is implemented to obtain optimal station Zenith Path Delay time series based on all observations within the observation session [12]. Precise corrections to orbits and clocks of satellites used are made available by the IGS. The mapping function used in CSRS-PPP is the Global Mapping function. As an input, single or dual-frequency GNSS data can be used. The user may choose NAD83 or ITRF2008 frame of reference to determine coordinates [12].
The service of MagicGNSS operated by the company GMV Aerospace and Defense is made available through their website where the user can process data in static and kinematic mode at two frequencies. The user can choose to use final and rapid products for corrections of orbits and clocks of satellites made accessible by the IGS or GMV, the user can choose either rapid or final products. The current version can process data from constellations, GPS, GLONASS, Galileo, BeiDou. Coordinates of the calculated position can be determined in two frames of reference ITRF2008 or ETRS89. MagicGNSS does not take into account in calculation parameters of the phase center antenna.
The following three software packages had to be installed in the computer and had to be run locally in the processing computer. For all of them, it is necessary to load the observation file, the navigation file and for PPP, the precise ephemeris are also needed. Additionally, other corrections such as Ocean Tide Loading or the parameters of the phase center antenna can be included. Since the user has to load all the files, there is flexibility to use corrections from different sources. For this study, all the final products from the IGS were used.
gLAB is an advanced interactive educational multipurpose package to process and analyze GNSS data [25] developed by Catalonia Technical University and the European Space Agency. It can process either single or dual-frequency GPS only data. The tropospheric delay is defined in terms of the elevation angle (El) of the satellite as: T r (El) = T rz,dry M dry (El) + T rz,wet M wet (El) (11) where T rz,dry and T rz,wet are the dry and wet slant tropospheric delay which can be estimated with a simple model: T rz,dry =∝ e βH (12) T rz,wet = T rz0,wet + ∆T rz,wet (13) where ∝= 2.3m, β = 0.116·10 −3 H is the height above sea level, in meters. T rz0,wet = 0.1 m and ∆T rz,wet is estimated as a random walk process in the navigation-Kalman-filter together with the coordinates and other parameters [25][26][27]. Or with the UNB-3 model. M dry and M wet are the dry and wet mapping function (Neill mapping function) which does not require any meteorological data. The multiplication of the mapping function and the slant delay yield the Zenith Troposheric Delay.
POINT is a software package developed by the University of Nottingham. It is capable of processing L1 and L2 GPS data. It implements an Extended Kalman Filter for positioning employing double differences observables [28]. The hydrostatic component of the ZTD is calculated using a model such as Saastamoinen, Hopfield or IFADIS and the Neill mapping function. The wet component is an unknown in the Extended Kalman Filter, the total zenith tropospheric delay is calculated as the sum of wet and dry components: RTKLIB is an open source positioning software developed by Takatsu. It can implement different positioning techniques, among them, PPP which can be computed in static or kinematic mode. All the corrections are input to the software via its Graphic User Interface. The effect of the troposphere is modelled using a mapping function and zenith tropospheric delays. The mapping function in terms of the elevation angle (El) and the azimuth angle (Az) between the satellite and the receiver is calculated as: and the tropospheric delay is calculated as: where Z T is the tropospheric zenith total delay in meters. This parameter is estimated from the Extended Kalman Filter together with the north component of tropospheric gradient (G N, ) and the east component of tropospheric gradient (G E, ). Z H is the tropospheric zenith hydro-static delay in meters which is calculated using a tropospheric model, either Saastamoinen, Hopfield or modified Hopfield model with the zenith angle z = 0 and relative humidity hrel = 0. M h (El) and M w (El) are the hydro-static and wet mapping function respectively. The Niell Mapping Function (NMF) is used in both cases. A summary of the capabilities of the software used for this study is shown in Table 1.
The three PPP software used for this experiment were: gLAB, POINT and RTKLIB which require precise ephemeris and other corrections to be input manually. In all cases, the corresponding IGS final clock and ephemeris files were used which contain data every 15 min. The three PPP software required the ANTEX file for antenna phase correction. No decimation was chosen for neither of the software, therefore, a solution was found for every epoch available in the observation file and the ephemeris file. The elevation mask was set to 10 degrees as indicated in Table 1 and GPS-only data was used in all software. The UNB-3 Nominal model was used together with the Neil Mapping function to model the effect of the troposphere in POINT and gLAB.
POINT requires more files to process data using the PPP technique. The differential code bias product (DCB) computed by the Institute of Geodesy and Geophysics of the Chinese Academy of Sciences was used. The Ocean Tide Loading (OTL) file obtained with the GOT00.2 model was also included and the same ANTEX file as for gLAB was used as well. The Saastamoinen model is used to compute the hydrostatic component of the tropospheric delay and the wet component is estimated as an unknown in the Extended Kalman Filter. RTKLIB had as input the same ephemeris file (Final from the IGS) and the same corrections as POINT, the same DCB, OTL and ANTEX files were used, it was run in PPP kinematic mode, the elevation mask was set to 10 degrees. The tropospheric effect was calculated as an unknown in the Extended Kalman Filter.
The data used in this experiment was collected from the IGS which operates over 400 GNSS stations across the world, it provides daily and hourly observation and navigation files for each station. Furthermore, the IGS provides other products, such as satellite ephemeris, earth rotation parameters and tropospheric delay with different latencies. The tropospheric delay product is generated from ground-based GNSS data with the Bernese GPS Software version 5.0, a cut-off angle of 7, IGS final satellite, orbit and EOP products are used for the computation [29]. Therefore, the IGS ZTD product is available approximately after three weeks after the observation date once the final products are available. The product contains the estimation of clock, position of the receiver which is presented as a constant, zenith delay in millimeters, which is estimated as a random walk with variance of 3 cm/h. Also included are the atmospheric gradients estimated as a random walk with variance of 0.3 cm/h. The temporal resolution of zenith day estimates is 5 min and the mapping function is the Global Mapping Function (GMF). For this study, observation data from 9 IGS stations listed in Table 2 were used as well as the IGS tropospheric product for 8 days. The days chosen for the study were the 27th calendar day of the month of January, April, July and October of the year 2016 and 2017 because these dates cover weather conditions during the four seasons in the different hemispheres.

Results
In order to assess the quality of the ZTD estimates of the 6 PPP software previously described, observation data from the 9 IGS stations were processed with each of the software packages and the ZTD was estimated in kinematic mode. The Root Mean Square Error (RMSE) is used as the quality indicator in this performance assessment, which is computed as: where n is the total number of ZTD estimates available. The RMSE was computed for each station and each software per day using Equation (17). All estimates available were used. The results are shown in Figures 1-8. The first three stations are in the northern hemisphere, the next three are near the equator and the last three are in the southern hemisphere, the data was grouped according to the station location and the RMSE of the differences for each group was calculated and the results are shown in Tables 3-6. Finally, Table 7 shows the RMSE of all differences estimated with each software and each online PPP service. Figures 1 and 2 represent the 27th of January of 2016 and 2017 respectively. This is the winter in the northern hemisphere and the summer in the southern hemisphere. In both figures, it can be seen that the estimates obtained with APPS and POINT for the stations MAL2 and NAUR are in both cases very far away from the value from the IGS tropospheric product. Also, station MAL2 and RIOP produce very high values of RMSE with POINT and APPS. However, in both days, the RMSE obtained with CSRS-PPP and MagicGNSS are lower than 5 cm for all stations. Furthermore, high RMSE values for the stations PARC, MAW1 and MAC1 are obtained with RTKLIB and POINT.
April 27th 2016 and 2017 are days 118 2016 and 117 2017 respectively. This is a day in spring in the northern hemisphere or autumn in the southern hemisphere. Therefore, mild changes of temperature are expected. The RMSE values of the differences between the estimated ZTD and the IGS tropospheric product's value are depicted in Figures 3 and 4 from where it can be seen that APPS, POINT and RTKLIB obtain a high RMSE value for the station MAL2 and NAUR. In contrast, the RMSE value is very low in all stations (<5 cm) for the estimates obtained with CSRS-PPP and MagicGNSS.                   Figures 1-6 show a trend that in most cases the RMSE obtained with data from stations near the equator is higher than for the other stations with most of the software used for implementing the PPP technique. In order to further study the quality of ZTD estimation at different latitudes, the stations were grouped in three groups according on their latitude and the GNSS data from all the stations in the same latitude was compared to the respective IGS Tropospheric Product by calculating the difference between both values. The RMSE of all the differences from the stations in the group was   Figures 1-6 show a trend that in most cases the RMSE obtained with data from stations near the equator is higher than for the other stations with most of the software used for implementing the PPP technique. In order to further study the quality of ZTD estimation at different latitudes, the stations were grouped in three groups according on their latitude and the GNSS data from all the stations in the same latitude was compared to the respective IGS Tropospheric Product by calculating the difference between both values. The RMSE of all the differences from the stations in the group was  Figures 1-6 show a trend that in most cases the RMSE obtained with data from stations near the equator is higher than for the other stations with most of the software used for implementing the PPP technique. In order to further study the quality of ZTD estimation at different latitudes, the stations were grouped in three groups according on their latitude and the GNSS data from all the stations in the same latitude was compared to the respective IGS Tropospheric Product by calculating the difference between both values. The RMSE of all the differences from the stations in the group was calculated. These results are presented in Tables 3-6. The regions are defined as: North: ALGO, REYK and TIXI, Center: MAL2, RIOP and NAUR and South: PARC, MAW1 and MAC1. According to the results shown in Tables 3-6 CSRS-PPP is the online software that performs the best for all stations because the RMSE is always lower than RMSE from other software. According to the same data, most of the software used for this study had the highest RMSE value with data from the equatorial stations. APPS, POINT and RTKLIB had their highest value of RMSE for the equatorial stations. MAGIC had most of its highest RMSE values in stations in the equatorial region except with data from July 2017, October 2016 and 2017. Similarly, CSRS had the highest RMSE for stations near the equator for five days, three days the highest RMSE was found for stations in the southern hemisphere. In contrast, GLAB did not have a clear pattern, three days the highest RMSE was found for the stations in the North, two days for the equatorial region and three days for the stations in the southern hemisphere.
In order to evaluate the quality of ZTD estimates of each software, the RMSE of all the differences (ZTD estimated from all stations with the same software minus IGS tropospheric product) was calculated. Its results are shown in Table 7 with all values in centimeters.  Table 7 shows that the two online services CSRS-PPP and MAGIC estimate the ZTD to a closer value to the IGS tropospheric product than the three other software packages for most cases the RMSE is equal or less than 1 cm (not the case with data from Day 27 of year 2017). From the three PPP software run locally GLAB is the one that had the lowest RMSE.

Discussion
Every software used for the analysis presented in this study uses a similar strategy to estimate the Zenith Tropospheric Delay, which is using a model to estimate the hydrostatic slant delay, use a mapping function to estimate the delay in the zenith direction and estimate the wet delay as an unknown in the parameter estimation process typically done with an Extended Kalman Filter. The online PPP software use the Global Mapping Function based on numerical weather model data while the locally run PPP processing software implement the Niell Mapping Function which depends only on the site coordinates and day of the year. Because the GMF involves the use of weather model data, it models better the delay caused by the troposphere which as seen in Equation (2) is affected by the meteorological variables near the receiver such as temperature, pressure, and partial water vapor pressure. The use of the GMF is one of the reasons why online services obtain an estimation closer to the IGS tropospheric product.
The effect of the ionosphere is another reason why the estimated value and the value of ZTD from the IGS tropospheric product are different. The model used to correct the ionosphere effect used by the online PPP services is not stated, for the other three locally-run software no ionosphere model is used, only the carrier phase and code combinations is done to obtain ionosphere-free pseudoranges, this combination only eliminates the first order ionosphere effect but residual effects are not eliminated and they can cause an effect on to the signal.
A third reason for the discrepancies between estimated ZTD and the values in the IGS tropospheric product is the different cut-off angle set up in the software. The IGS tropospheric product has a cut-off angle of 7. All the software packages and online PPP services used for this study allow to set the cut-off angle. However, the RTKLIB version used for this study allows to use only angles multiple of 5, therefore, 10 was used as the closest option. Also, in the other software the cut-off angle was set to 10 degrees. It is possible that this 3 degrees difference has an effect on the ZTD estimation because some satellites might be discarded for the solution. Furthermore, multipath affects the signal as well.
According to the results presented in this study, POINT and RTKLIB had very high RMSE values for stations near the equator which means that the model currently used does not clearly represent the tropospheric effect at these latitudes, possible reasons is that the thickness of the troposphere on the equatorial region is higher than in the polar region and different weather conditions near the equator. Also, APPS and MAGIC obtained the highest RMSE values for stations in the center for days 27 2016 and 2017, 118 2016, 117 2017 and 209 2016, however for the other days the highest values were found for the southern and northern hemisphere respectively which confirms that the equatorial region has specific atmospheric conditions that are not properly accounted for with models and the parameters used for the estimation. CSRS-PPP and GLAB obtained high RMSE values with data from different latitudes at different days. However, the results with CSRS-PPP were always less than 1 cm (except for day 27 2017) and in the case of GLAB the results were in the range of 2 and 6 cm always which means that both software obtained estimates very close to the IGS tropospheric value with all data. The solution with GLAB takes one epoch to converge, so the first epoch is not considered in the RMSE analysis.
This study only included eight days of data, however the same days in different years were chosen, also, the stations chosen are distributed throughout the globe. It is expected to find similar weather conditions the same days of two different years, therefore a trend can be found of how close the estimations are to the IGS tropospheric product. The stations are located in different latitudes which allows to study how the different models used for the tropospheric model are influenced by the latitude of the station.

Conclusions
In this paper, a comparison analysis of the estimated Zenith Tropospheric Delay (ZTD) obtained with 6 Precise Point Positioning (PPP) post-processing software and the International GNSS Service (IGS) tropospheric product is presented. The estimated ZTD obtained with APPS, CSRS-PPP, MagicGNSS, POINT, RTKLIB and gLAB were compared with the ZTD provided by IGS. The Root Mean Square Error (RMSE) was used as the indicator of accuracy of the estimation because it indicates how different the estimated value is from the ground truth.
Three trends were found in this study, first, it was found that CSRS-PPP obtains ZTD estimates very close to the value from the IGS Tropospheric product. Second, it was found that the tropospheric models currently implemented in RTKLIB and POINT do not account properly to the weather and atmospheric conditions in the equatorial region. The corrections used by CSRS-PPP and MAGIC are very precise so estimates closer to the truth value were found. The third trend found was that GLAB also estimates the ZTD to a value very close to the IGS Troposheric product.
The season change did not have a big impact on the ZTD estimation by PPP software. With the selected data sets. If precise ZTD estimates are needed for GNSS meteorology or numeric weather models, CSRS-PPP can provide very accurate estimates followed by MagicGNSS and gLAB.