Real-Time Tropospheric Delay Retrieval from Multi-GNSS PPP Ambiguity Resolution: Validation with Final Troposphere Products and a Numerical Weather Model

Abstract

The multiple global navigation satellite systems (multi-GNSS) bring great opportunity for the real-time retrieval of high-quality zenith tropospheric delay (ZTD), which is a critical quality for atmospheric science and geodetic applications. In this contribution, a multi-GNSS precise point positioning (PPP) ambiguity resolution (AR) analysis approach is developed for real-time tropospheric delay retrieval. To validate the proposed multi-GNSS ZTD estimates, we collected and processed data from 30 Multi-GNSS Experiment (MGEX) stations; the resulting real-time tropospheric products are evaluated by using standard post-processed troposphere products and European Centre for Medium-Range Weather Forecasts analysis (ECMWF) data. An accuracy of 4.5 mm and 7.1 mm relative to the Center for Orbit Determination in Europe (CODE) and U.S. Naval Observatory (USNO) products is achievable for real-time tropospheric delays from multi-GNSS PPP ambiguity resolution after an initialization process of approximately 5 min. Compared to Global Positioning System (GPS) results, the accuracy of retrieved zenith tropospheric delay from multi-GNSS PPP-AR is improved by 16.7% and 31.7% with respect to USNO and CODE final products. The GNSS-derived ZTD time-series exhibits a great agreement with the ECMWF data for a long period of 30 days. The average root mean square (RMS) of the real-time zenith tropospheric delay retrieved from multi-GNSS PPP-AR is 12.5 mm with respect to ECMWF data while the accuracy of GPS-only results is 13.3 mm. Significant improvement is also achieved in terms of the initialization time of the multi-GNSS tropospheric delays, with an improvement of 50.7% compared to GPS-only fixed solutions. All these improvements demonstrate the promising prospects of the multi-GNSS PPP-AR method for time-critical meteorological applications.

1. Introduction

Global Positioning System (GPS) meteorology, which is an important approach used for the remote sensing of atmospheric water vapor, was first proposed by Bevis et al. [1]. Significant development has been achieved in the past few years [2,3,4,5,6,7,8]. In comparison with established atmospheric sounding techniques, GPS has extensive application prospects for sounding water vapor since it has the unique advantages of all-weather availability, high accuracy, long-term stability, and high time resolution [9,10]. The zenith tropospheric delay (ZTD), which is the basic observation of GPS meteorology, is commonly utilized to quantify precipitable water vapor (PWV) [11,12]. The GPS-based tropospheric delay and precipitable water vapor have been assimilated into a numerical weather model (NWM) for the improvement of nowcasting of severe rainfall [13,14]. The widespread application of GPS-derived tropospheric delay at several weather agencies has demonstrated that GPS meteorology has the potential to provide accurate monitoring quantities of water vapor and validate the contribution of tropospheric delay for time-critical meteorological studies [15,16].

To meet different requirements of meteorological applications, GPS water vapor derived from both post-processing and near-real-time modes is assimilated into regional and global numerical weather forecast models [13,17]. However, several important and time-critical meteorological applications, for example, the sounding of short-time weather variation, could benefit from the provision of the atmospheric state with a more rapid update rate, which can be from real-time processing [18]. For example, tropospheric delays with latency of a few minutes and an accuracy of 5–30 mm have been applied for nowcasting models already [19].

In terms of the real-time ZTD estimation method, precise point positioning (PPP) using precise real-time clock and orbit products has significant advantages in processing efficiency and flexibility [20]. The accuracy of real-time PWV from GPS PPP solutions better than 3 mm is qualified for weather nowcasting, which also indicates that the real-time PPP technique can monitor atmospheric water vapor effectively. The PPP ambiguity resolution (AR) technique [21], which was proposed to shorten initialization time and improve positioning accuracy [22,23], can further improve the real-time tropospheric products. The integrated water vapor obtained by GPS PPP AR with an accuracy of 1–2 mm with respect to the water vapor radiometer measurements is valuable for meteorological studies [9]. The PPP tropospheric delays estimated by three software products, including BKG Ntrip Client (BNC) 2.7, G-Nut/Tefnut, and PPP-Wizard, have revealed that centimeter-level accuracy can be achieved for PPP float solutions and mm-level improvement can be observed on the tropospheric delays by ambiguity resolution [24].

With the development and improvement of multiple GNSSs (multi-GNSS), 88 satellites (32 GPS, 24 Global Orbiting Navigational Satellite System (GLONASS), 18 Galileo, and 14 BeiDou Navigation Satellite System (BDS)) are available currently and the number of navigation satellites will exceed 120 when the GNSS is fully operational. To make full use of the upcoming new satellite systems and signals, the Multi-GNSS Experiment (MGEX) program has been implemented since 2012. With the increased number of navigation satellites, multi-GNSS can provide great potential for geosciences applications. The European Co-operation in the field of Scientific and Technical Research (COST) Action ES1206 “Advanced Global Navigation Satellite Systems tropospheric products for monitoring severe weather events and climate (GNSS4SWEC)” is about advanced real-time tropospheric products retrieved from multi-GNSS observations [25,26]. The ZTDs retrieved from multi-GNSS observations with several millimeters accuracy would complement existing atmospheric observation systems and be widely applied in meteorological applications [27,28].

Currently, multi-GNSS-based ZTD is mainly retrieved from float PPP solutions. In this contribution, for the first time, the multi-GNSS PPP-AR approach is developed for retrieving high-quality real-time ZTD. Observation data collected from 30 MGEX stations are processed to evaluate our new zenith tropospheric delay products. The post-processed final tropospheric delay products of the Center for Orbit Determination in Europe (CODE) and the U.S. Naval Observatory (USNO) and products from the European Centre for Medium-Range Weather Forecasts analysis (ECMWF) are utilized to validate the real-time tropospheric delay products retrieved from multi-GNSS PPP-AR. The contribution of the multi-GNSS PPP-AR method to improvements of real-time tropospheric delay retrieval is demonstrated.

2. Real-Time Sensing of Tropospheric Delay from Multi-GNSS PPP-AR

The ionosphere-free (IF) combination is usually applied to the PPP model to correct the first-order ionospheric delays. The IF PPP model for multi-GNSS processing can be formulated as

$$\{\begin{array}{l}{P}_{r,IF}^G={{\rho }_r^G}_g+c(t_r-t^G)+c(b_{r,IF}^G-b_{IF}^{s,G})+T_r^G+e_{r,IF}^G\\ P_{r,IF}^C={{\rho }_r^C}_g+c(t_r-t^C)+c(b_{r,IF}^C-b_{IF}^{s,C})+T_r^C+e_{r,IF}^C\\ P_{r,IF}^{R_k}={{\rho }_r^{R_k}}_g+c(t_r-t^{R_k})+c(b_{r,IF}^{R_k}-b_{IF}^{s,R_k})+T_r^{R_k}+e_{r,IF}^{R_k}\\ P_{r,IF}^E={{\rho }_r^E}_g+c(t_r-t^E)+c(b_{r,IF}^E-b_{IF}^{s,E})+T_r^E+e_{r,IF}^E\end{array}$$

(1)

$$\{\begin{array}{l}{L}_{r,IF}^G={{\rho }_r^G}_g+c(t_r-t^G)+{\lambda }_{IF,G}(N_{r,IF}^G+B_{r,IF}^G-B_{IF}^{s,G})+T_r^G+{\epsilon }_{r,IF}^G\\ L_{r,IF}^C={{\rho }_r^C}_g+c(t_r-t^C)+{\lambda }_{IF,C}(N_{r,IF}^C+B_{r,IF}^C-B_{IF}^{s,C})+T_r^C+{\epsilon }_{r,IF}^C\\ L_{r,IF}^{R_k}={{\rho }_r^{R_k}}_g+c(t_r-t^{R_k})+{\lambda }_{IF,R_k}(N_{r,IF}^{R_k}+B_{r,IF}^{R_k}-B_{IF}^{s,R_k})+T_r^{R_k}+{\epsilon }_{r,IF}^{R_k}\\ L_{r,IF}^E={{\rho }_r^E}_g+c(t_r-t^E)+{\lambda }_{IF,E}(N_{r,IF}^E+B_{r,IF}^E-B_{IF}^{s,E})+T_r^E+{\epsilon }_{r,IF}^E\end{array}$$

(2)

where the superscripts $G$, $C$, $R$, and $E$ refer to the GPS, BDS, GLONASS, and Galileo satellites, respectively; $c$ refers to the speed of light in vacuum; $k$ denotes the frequency factor and $r$ refers to the receiver; $P_{r,IF}$ and $L_{r,IF}$ represent the IF observations for the code and phase, respectively, in meters; ${\lambda }_{IF}$ denotes the wavelength of the IF combination in meters; $t^s$ is the satellite clock offset while $t_r$ refers to the receiver clock error in seconds; $b_{r,IF}$ and $b_{IF}^s$ are code hardware delays in seconds while $B_{r,IF}$ and $B_{IF}^s$ are the phase delays of the IF combination; ${\rho }_{r,g}$ denotes the geometric distance; and $N_{r,IF}$ is the IF ambiguity. The code hardware delays ($b_{r,IF}^G$, $b_{r,IF}^{R_k}$, $b_{r,IF}^E$, $b_{r,IF}^C$) and phase delays ($B_{r,IF}^G$, $B_{r,IF}^{R_k}$, $B_{r,IF}^E$, $B_{r,IF}^C$) of different satellite systems are different due to the different frequencies. The GPS code hardware delay and phase delay are set to zero and the inter-system biases (ISBs) are introduced for the BDS and Galileo systems, respectively, while the inter-frequency bias (IFB) is estimated for each frequency of GLONASS in multi-GNSS PPP processing. $T_r^s$, the slant tropospheric delay, which is of main interest in this study, can be expressed as

$$T_r^s=M{h}_r^s\cdot ZH{D}_r+M{w}_r^s(ZW{D}_r+\cot (e)(G_N\cdot \cos (a)+G_E\cdot \sin (a)))$$

(3)

in which $M{h}_r^s$ denotes the hydrostatic coefficient of the mapping function and $ZH{D}_r$ denotes the zenith hydrostatic delay (ZHD) in meters; $M{w}_r^s$ denotes the wet coefficients of the mapping function and $ZW{D}_r$ denotes the zenith wet delay (ZWD) in meters; $e$ and $a$ denote the elevation and azimuth angle, respectively; and $G_E$ and $G_N$ refer to the east and north gradient, respectively. The Saastamoinen model [29] is used to model the ZHD, while the ZWD and gradient components can be estimated in PPP processing.

The uncalibrated phase delay (UPD), which destroys the integer nature of ambiguity, should be precisely estimated to enable ambiguity fixing and then provided to users. For multi-GNSS PPP-AR processing, the UPD products of the four systems should be estimated first [30,31]. The IF combination ambiguity in dual-frequency PPP is commonly described as the combination of wide-lane (WL) and narrow-lane (NL) ambiguity:

$$N_{r,IF}=\frac{f_1}{f_1+f_2}{N}_{r,nl}+\frac{f_1{f}_2}{f_1^2-f_2^2}{N}_{r,wl}$$

(4)

where $N_{r,nl}$ denotes the NL ambiguity and $N_{r,wl}$ represents the WL ambiguity. We obtain the WL ambiguity from Hatch–Melbourne–Wübbena (HMW) combination observations [32,33,34]. With multipath, noise, and delays included, the WL ambiguities do not have an integer nature here. After the correction of WL UPDs, the WL ambiguity can be fixed to an integer by the round strategy [35]. Then, the fixed WL ambiguity and IF combination ambiguity from the PPP solutions are used for the narrow-lane ambiguity. When the WL ambiguity and the NL ambiguity are both fixed, the fixed solution of the IF ambiguity for PPP can be obtained.

A four-system UPD estimation model is proposed in this paper to obtain the real-time wide-lane and narrow-lane UPDs. If there are $n$ stations and $m$ satellites that can be observed at each station, the four-system UPD estimation model proposed here can be formulated as the following [30]:

$$\left[\begin{array}{c}{\mathbf{D}}_1\\ {\mathbf{D}}_2\\ ⋮\\ {\mathbf{D}}_n\end{array}\right]=\left[\begin{array}{ccccc}{\mathbf{R}}_{1G}& {\mathbf{R}}_{1C}& {\mathbf{R}}_{1R}& {\mathbf{R}}_{1E}& {\mathbf{s}}_1\\ {\mathbf{R}}_{2G}& {\mathbf{R}}_{2C}& {\mathbf{R}}_{2R}& {\mathbf{R}}_{2E}& {\mathbf{s}}_2\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathbf{R}}_{nG}& {\mathbf{R}}_{nC}& {\mathbf{R}}_{nR}& {\mathbf{R}}_{nE}& {\mathbf{s}}_n\end{array}\right]•\left[\begin{array}{c}{\mathbf{d}}_G\\ {\mathbf{d}}_C\\ {\mathbf{d}}_R\\ {\mathbf{d}}_E\\ {\mathbf{d}}^s\end{array}\right]$$

(5)

where the matrix $\mathbf{D}$ with $m$ rows and 1 column denotes the ambiguity fractional part for each satellite-station link; ${\mathbf{d}}_G$, ${\mathbf{d}}_R$, ${\mathbf{d}}_E$, and ${\mathbf{d}}_C$ represent matrices of UPDs at the receiver site for the GPS, BDS, GLONASS, and Galileo systems; and ${\mathbf{d}}^s$ denotes the matrices of satellite UPDs. ${\mathbf{R}}_i$ ($m$ rows and $n$ columns) and $S_i$ ($m$ rows and $m$ columns) refer to the coefficient matrices of UPDs at the receiver and satellite sides. All of the ambiguities of the global stations at the current epoch are processed together to estimate satellite and station UPDs based on the least squares method. To eliminate the rank-deficiency for a multi-GNSS UPD estimation equation, one station or one satellite is selected as the datum for each individual satellite system with the corresponding UPD set as zero.

The satellite UPDs can be estimated epoch-by-epoch based on an IF PPP model and then applied back as corrections to achieve real-time ambiguity resolution. Several specific issues of the multi-GNSS UPD estimation for each system should be mentioned. Owing to the frequency division multiple access (FDMA) strategy of GLONASS, the signal delays of different satellites vary and the IFB problems will affect the estimation of WL UPDs. To solve the problems caused by the IFB, observations from the homogeneous receivers are used here to estimate the GLONASS UPDs. For BDS, the satellite-induced code bias of Inclined Geosynchronous Orbit (IGSO) and Medium Earth Orbit (MEO) satellites should be eliminated by the elevation-dependent model proposed by Wanninger and Beer [36]. Meanwhile, for GEO satellites, code bias is corrected by a sidereal wavelet filter.

To obtain real-time ZTD retrievals from multi-GNSS PPP-AR, multi-GNSS real-time orbits and clock and UPD corrections are determined first using observation streams from the global network. With the satellite UPD corrected, the fractional parts of all corrected WL ambiguities are averaged to obtain the WL UPD of a receiver using the formulation by Gabor and Nerem [37,38]. The WL UPDs at the satellite side and the receiver side can be further applied to the corrected WL ambiguities at the current epoch. The receiver NL UPD can be separated from the NL ambiguity with the selected reference ambiguity set as its nearest ambiguity. With the WL and NL UPDs separated, the integer property of the WL and NL ambiguities can be recovered and then the ambiguities can be fixed to integers. The coordinates of stations are commonly tightly constrained to known values to obtain the real-time ZTD retrievals. After applying the real-time multi-GNSS orbit, clock, and UPDs corrections, the parameters $\mathbf{X}$ to be estimated will be

$$\mathbf{X}=(t_r,b_{r,IF}^C,b_{r,IF}^{R_k},b_{r,IF}^E,ZW{D}_r,G_N,G_E,{\mathbf{N}}_r^s).$$

(6)

The multi-GNSS processing is performed by a sequential least squares estimator to obtain the estimated parameters. The receiver-dependent clock offset $t_r$ is estimated as white noise. We set the code bias of GPS as zero, i.e., $b_{r,IF}^C$, $b_{r,IF}^E$, and $b_{r,IF}^{R_k}$ are relative to the GPS bias $b_{r,IF}^G$. These ISBs and IFBs are estimated as constants; ${\mathbf{N}}_r^s$ denotes the integer ambiguities, which are fixed to integers using the LAMBDA method [39]. The ratio test was used to validate the ambiguity validation with a threshold of 2 [40]. In real-time data processing, both the ZWD and horizontal gradient are estimated as a random walk process. For ZWD, it can be expressed as

$$ZW{D}_r(t_{i+1})=ZW{D}_r(t_i)+\epsilon ,\epsilon ~N(0,{\delta }_0^2\cdot \Delta t)$$

(7)

where $t_i$ and $t_{i+1}$ are two adjacent epochs and $\epsilon$ is the temporal variation of ZWD between $t_i$ and $t_{i+1}$; ${\delta }_0$ refers to the noise intensity of the zenith wet delay, which is set as $3\mathrm{mm}/\sqrt{\mathrm{hour}}$ in this study; and $\Delta t$ is the time difference between two adjacent epochs. The multi-GNSS data processing strategy used for real-time ZWD estimation is summarized in

Table 1. Multi-GNSS processing strategy for real-time ZWD estimation.

Item	Processing Strategies
Estimator	Sequential least squares estimator
Observations	Observation from GPS/GLONASS/BDS/Galileo
Signals	GPS and GLONASS: L1/L2; BDS: B1/B2; Galileo: E1/E5a
Sampling rate	30 s
Elevation cutoff	7°
Weight for observations	Elevation-dependent weighting strategy.
Satellite orbit	Fixed
Satellite clock	Fixed
Zenith Tropospheric delay	Initial model + random walk model
Tropospheric gradients	Random walk model
Mapping function	Global Mapping Function (GMF) [41]
Phase-windup effect	Corrected
Receiver clock	Estimated, white noise
ISB and IFB	Estimated as constant, GPS as reference
Station displacement	Solid Earth tide, pole tide, ocean tide loading International Earth Rotation and Reference Systems Service (IERS) Convention 2010 [42]
Satellite antenna phase center	Corrected
Receiver antenna phase center	Corrected
Coordinates of stations	Fixed
Ambiguities	PPP ambiguity resolution is applied

.

With an a priori zenith hydrostatic delay and wet delay parameter, the ZTD can be reconstructed according to

$$ZTD=ZHD+ZWD.$$

(8)

To retrieve the actual ZWD from the PPP-derived ZTD, a new zenith hydrostatic delay is calculated by using pressure parameters from the ECMWF analysis. Based on the accurate ZTD and ZHD, the ZWD, which will be converted into precipitable water vapor, can be calculated precisely. The multi-GNSS observation processing procedure for real-time tropospheric delay retrieval is shown in the Figure 1. The stochastic model is based on the precision of observations with the elevation-dependent weighting strategy $Q=1/{\sin }^2(ele)$, in which $ele$ is the satellite elevation. The elevation cut-off angle was set as 7°. The precision for phase raw observables and code raw observables for both the GPS and GLONASS systems is 3 mm and 0.3 m, respectively, while it is 6 mm and 0.6 m for the BDS and Galileo systems, respectively.

3. Data Collection

3.1. Multi-GNSS Data

To capture new satellite signals and achieve a good geographic coverage for all constellations, International GNSS Service (IGS) [43] has initiated the Multi-GNSS Experiment since 2012 [44]. With the development of MGEX, multi-GNSS stations have been deployed around the world, working with existing IGS stations to provide abundant observation data for users. Since August 2017, 212 globally distributed MGEX stations equipped with all kinds of GNSS receivers, such as Trimble NetR9, Septentrio, and Leica receivers [45], have been providing high-quality multi-GNSS observation data to users incessantly (http://www.igs.org/network). As for the estimation of multi-GNSS tropospheric delay, we collected observation data from 30 MGEX stations where GPS, GLONASS, Galileo, and BDS observations were simultaneously available. The distribution of stations in our study is shown in Figure 2.

3.2. Final Troposphere Products

IGS operates a global network of GNSS ground stations, data centers, and data analysis centers (ACs) to provide data and generate data products, such as GNSS satellite orbit and clock products, global ionosphere maps (GIM), and ZTD products. More than 200 organizations and institutes participate in MGEX, providing and sharing GNSS data, products, and services to support high-precision GNSS applications. To date (August 2017), there are two ACs which generate two types of post-processed tropospheric delay products to provide a reference of tropospheric estimates, and they are the U.S. Naval Observatory and the Center for Orbit Determination in Europe.

The CODE ZTD products sample every 2 h. The data obtained from about 250 stations are processed in a post-processing network solution by using Bernese GNSS Software Version 5.3. The elevation cutoff angle is set as 3°. The global pressure and temperature (GPT) empirical model and the global mapping function (GMF) [41] are used in CODE products for modeling the tropospheric delays. The final ZTD products using the Bernese GPS Software Version 5.0 based on GNSS observational data of 350 stations are provided by USNO. The time resolution for USNO ZTD products is 5 min. The GMF model is also adopted by USNO [46]. The accuracy of the two mentioned ZTD products can be up to 4 mm with respect to the tropospheric products generated from different geodetic measurements, for instance, very long baseline interferometry (VLBI), Delft object-oriented radar interferometric software (DORIS), and radiosondes [19,45]. The Standard Deviations (STD) of the CODE ZTDs are about 0.3–1.0 mm while the corresponding statistics for the USNO ZTDs are about 1–2 mm [47]. Since both products are calculated from GNSS observations, which is the same as our approach, it would be a good opportunity to evaluate our data processing strategy and results.

3.3. ECMWF Data

Numerical weather models (NWM) with various meteorological observations are able to provide neutral atmosphere information and compute weather forecasts. The meteorological parameters for any location at any time can be obtained from the NWM, and tropospheric delay is calculated for each station-satellite link with the trace algorithms [48,49]. As important NWM-based tropospheric products, parameters from ECMWF (http://www.ecmwf.int/) were used for the provision of ZTD products in this study. Using the final ZTD product as reference, ZTD parameters obtained from the ECMWF were also evaluated and a good consistency was discovered between the two ZTD products [50]. Therefore, it is reasonable to use the ECMWF analysis data to evaluate our multi-GNSS PPP-AR ZTD estimates.

4. Results and Validations

A 30-day experiment from 1–30 January 2017 was conducted for the performance evaluation of real-time ZTD, which is estimated using our proposed multi-GNSS PPP-AR method. Thirty globally distributed MGEX stations were selected with simultaneous availability of GPS, GLONASS, Galileo, and BDS observations. Different data-processing modes, as listed in

Table 2. List of Data Processing Models.

Models	Details
G	Float PPP solution based on GPS-only
GR	Float PPP solution based on GPS/GLONASS
GE	Float PPP solution based on GPS/Galileo
GC	Float PPP solution based on GPS/BDS
GREC	Float PPP solution based on GPS/GLONASS/Galileo BDS
G-AR	Fixed PPP solution based on GPS-only
GR-AR	Fixed PPP solution based on GPS/GLONASS
GE-AR	Fixed PPP solution based on GPS/Galileo
GC-AR	Fixed PPP solution based on GPS/BDS
GREC-AR	Fixed PPP solution based on GPS/GLONASS/Galileo/BDS

, were adopted to assess the benefits of multi-GNSS and ambiguity resolution. The real-time tropospheric products were estimated every 30 s. Then, the estimated real-time ZTDs were validated with the final tropospheric products and ECMWF data. The performance in terms of initialization time and accuracy of estimated ZTD products from all processing modes was analyzed and evaluated.

4.1. Initialization Analysis

To evaluate the performance of the retrieved real-time ZTD based on multi-GNSS PPP-AR, we investigated the initialization time of the ZTD series first. Here, initialization time denotes the time it takes for the accuracy of ZTD estimates to be better than 20 mm and keep within 20 mm with respect to the USNO ZTD products. All collected observation data are processed every 2 h separately. About 360 2-h sessions for each station for a period of 30 days are available to obtain the average initialization time. Here, six processing modes, including GPS-only float solution (G), GPS-only fixing solution (G-AR), GPS+GLONASS fixing solution (GR-AR), GPS+Galileo fixing solution (GE-AR), GPS+BDS fixing solution (GC-AR), and multi-GNSS fixing solution (GREC-AR), were applied to the initialization time evaluation of real-time ZTD products retrieved from multi-GNSS PPP-AR.

The real-time ZTDs of stations DUND and ZIM2 in all data-processing modes are presented in Figure 3. An initialization period is clearly visible, and the initialization time was found to be within 10 min for all processing modes. Compared to GPS-only, the dual-system fixed solutions require less initialization time, especially for the GR fixed-solution. The initialization time is effectively shortened for the GREC fixing solution with the minimum being less than 5 min. The multi-GNSS PPP-AR solutions ensure more stable and reliable results and can significantly mitigate sudden fluctuations, which are commonly observed in GPS-only float solutions.

Figure 4 presents real-time PPP-inferred troposphere results at the ZIM2 station from 1:00 a.m. to 7:00 a.m. on Day of Year (DOY) 002 of 2017. The estimator is restarted every hour to make the initialization process clearer. As shown in the figure, the troposphere results of single-, dual-, and four-system PPP solutions (G, G-AR, GR-AR, GE-AR, and GREC-AR) can all converge to a stable value after an initialization period of 1 h. Under the dual- and four-system environments, the initialization process becomes shorter and the four-system PPP AR results present the fastest initialization process. It can be noticed that there are still two sessions in which the initialization time exceeds 30 min for GPS-only solutions, while the GREC-AR solutions can converge to a stable value in a very short period in all sessions.

The initialization time of each station in different processing modes is plotted as boxplots in Figure 5. The box represents the 25%, 50%, and 75% quantiles. The whiskers are set to the 5% and 95% quantiles. The rhombus corresponds to the mean value. It is visible that the four-system PPP AR solutions show the smallest mean value and 50% quantiles with fewer outliers, which means that the fastest initialization process can be achieved by the multi-GNSS PPP-AR. The average initialization time of the single-, dual-, and four-system retrieved ZTD is summarized as

Table 3. Average root mean square (RMS) of the initialization time of real-time ZTD for single-, dual-, and four-system solutions with respect to the U.S. Naval Observatory (USNO) product (unit: min).

Solutions	Average Initialization Time
G-F	10.1
G-AR	9.8
GR-AR	5.1
GE-AR	9.1
GC-AR	8.9
GREC-AR	4.8

. For the GPS-only float results of ZTD estimates, the average initialization time is 10.1 min. Through the use of ambiguity resolution, the initialization time becomes shorter for most stations and the average initialization time is 9.8 min. In comparison with the ZTD series from GPS PPP AR, the initialization process is accelerated by multi-GNSS PPP ambiguity resolution. The average initialization time for all 2 h sessions of GR-, GE-, and GC-derived ZTD is 5.1, 9.1, and 8.9 min, with improvements of 47.9%, 7.1%, and 9.1%, respectively. Moreover, real-time ZTDs from four-system PPP AR require the shortest time for initialization, for which initialization can be achieved on average after 4.8 min. Compared to the GPS-only fixed solution, the initialization speed is improved by 50.7% by multi-GNSS PPP-AR processing. It can be found that although ambiguity resolution can accelerate the initialization process, improvement on observing geometry is more effective than ambiguity resolution to accelerate the initialization process of troposphere estimates.

4.2. Accuracy Validation with Final Troposphere Products

To analyze the accuracy of ZTDs derived from different solutions, a continuous period of 30 days was processed and compared with two different types of reference data; one comes from final troposphere products and the other is from ECMWF data. Here, the accuracy is defined as root mean square (RMS) values of difference between the estimated ZTD and reference data. Since the sample interval is different for different reference data, only the data in the same epochs were compared to eliminate the impact of interpolation on the evaluation results. In this section, the accuracy of real-time tropospheric delays with respect to two final troposphere products (CODE and USNO) are compared and analyzed first.

Figure 6 shows the real-time ZTDs derived from GPS PPP AR and multi-GNSS PPP-AR during the 30-day period from 1 January 2017 to 30 January 2017 at two MGEX stations BOR1 and CHTI. As a reference, the corresponding ZTD time-series of USNO and CODE are also shown in Figure 6. As expected, the real-time ZTDs from GPS PPP AR and multi-GNSS PPP-AR all reveal a great agreement with the two standard post-processed products from CODE and USNO. However, the ZTD series from GPS PPP AR exhibits larger noise and more outliers while ZTDs derived from multi-GNSS PPP-AR show more stable results with fewer outliers. The real-time ZTD differences at the KIRI station with respect to the final troposphere products from USNO are presented in Figure 7. The G-AR, GR-AR, GC-AR, GE-AR, and GREC-AR solutions are shown as red, green, purple, yellow, and blue dots. It can be found that the average ZTD difference of all solutions is close to zero without obvious deviation, which indicates that the real-time ZTDs from single-, dual-, and multi-GNSS ambiguity-fixed solutions all agree well with the USNO products. Additionally, the real-time ZTDs from the four-system PPP AR show the most stable result with the ZTD difference less than 10 mm during the period.

The root mean square (RMS) statistics of the deviations between the real-time ZTDs and the IGS final tropospheric products from CODE and USNO are computed for each station during the period from 1 January 2017 to 30 January 2017. The average RMS values of all 30 stations are calculated and results from DOY 001–010, 2017 with respect to CODE and USNO are presented in Figure 8 and Figure 9, respectively. The RMS statistics of ZTDs of different days agree well with each other. Results from all ten days reveal that the RMS values of real-time ZTDs were decreased by multi-GNSS fusion and PPP ambiguity resolution can further promote the results of real-time ZTDs. The RMS value of each station during the 30 days is also calculated and presented in Figure 10. The RMS values of real-time ZTD with respect to the USNO products for 30 stations range from 3 mm to 15 mm. The dual- and the four-system solutions exhibit a smaller RMS than the GPS-only solutions due to the larger number of satellites and the improvement in the observation geometry. After ambiguity resolution, the accuracy of the real-time ZTD can be further improved and the best performance can be achieved by the multi-GNSS PPP-AR.

The average RMS values of 30 days and 30 stations for all of the processing modes are calculated and shown in Figure 11 and

Table 4. Average RMS of ZTDs calculated from observations at 30 stations for 30 days (unit: mm).

System	PPP_CODE	PPP_AR_CODE	PPP_USNO	PPP_AR_USNO
G	7.1	6.6	8.9	8.5
GR	5.3	5.2	7.9	7.6
GE	6.4	5.9	8.7	8.3
GC	5.6	5.2	8.5	8.1
GREC	4.8	4.5	8.0	7.1

. From Figure 11, the RMS values of all solutions appear to be smaller than 8 and 10 mm with respect to CODE and USNO, respectively. After applying ambiguity resolution, the accuracy of the GPS-only and multi-GNSS solutions are all obviously improved. Taking the four-system solution as an example, the average RMS values of the GREC float solutions are 4.8 and 8.0 mm relative to the CODE and USNO products, respectively. With the ambiguity fixed, the GREC fixed solutions are improved up to 4.5 and 7.1 mm, showing an improvement of 6.6% and 11.3%, respectively. For both float and fixed results, ZTD estimates derived from dual-system solutions obtain higher accuracy than those from GPS-only solutions. Furthermore, the four-system PPP-derived ZTD achieves the best accuracy thanks to the most visible satellites and the best observable geometry. Taking the CODE product as a reference, the RMS value of the ZTD estimates for GPS-only fixed solutions is 6.6 mm, whereas for dual-system fixed results the values decrease to 5.2, 5.2, 5.9, and 4.5 mm for GC, GR, GE, and GREC, respectively. A significant improvement can be noticed for the accuracy of ZTDs derived from GREC PPP AR with respect to the single- and dual-system results. The ZTDs derived from GC PPP AR show similar accuracy to the GR solutions, while the GE results are a little worse. Compared to the GPS-only fixed results, the RMS values are reduced by 16.7% and 31.7% with respect to USNO and CODE by applying the multi-GNSS PPP-AR method.

4.3. Accuracy Validation with ECMWF Data

In this analysis, real-time tropospheric delays, retrieved from multi-GNSS PPP-AR processing, were validated by the ZTDs generated by the ECMWF parameters. Since ECMWF data are available every 6 h, we only compare the results at times 00:00, 06:00, 12:00, and 18:00 UTC. The ZTD values retrieved from multi-GNSS PPP-AR and ECMWF at station AUCK and PARK during DOY 001–030, 2017 are shown in Figure 12. It can be found that GNSS-derived ZTDs exhibit great agreement with the ZTD series from ECMWF. Because of fast changes in the humidity around the station, the ZTD time-series changes rapidly, which can be tracked accurately by both multi-GNSS and ECMWF results. We also calculated the RMS of the ZTDs differences for single-, dual-, and four-system solutions with respect to the ECMWF data. The average RMSs of all data-processing modes within our experiment span are summarized in

Table 5. Average RMS of the ZTDs in all data-processing modes with respect to ECMWF data (unit: mm).

System	PPP	PPP_AR
G	14.3	13.3
GR	13.6	12.6
GE	14.0	12.9
GC	13.7	12.7
GREC	13.2	12.5

. It can be obviously noticed that the ambiguity resolution contributes to GNSS-based ZTD estimates for all processing modes. The accuracy of the four-system fixed solution is 12.5 mm with respect to ECMWF data while the average RMS is 14.3 mm for GPS-only float results.

5. Conclusions

A new approach for real-time ZTD estimates based on multi-GNSS PPP ambiguity resolution was proposed in this paper. Observations from 30 globally distributed stations during a period of 30 days were employed in this paper to validate the contribution of the proposed method to real-time ZTD retrieval. Considering the post-processed final CODE/USNO ZTD products and ECMWF data as the references, we evaluated the initialization time and accuracy of real-time ZTD estimates for various PPP float or fixed solutions. Some specific validation results were given in this study.

The real-time ZTDs derived from multi-GNSS fixed solutions show a high consistency with CODE and USNO products. Accuracies of 4.5 mm and 7.1 mm with respect to CODE and USNO products were achieved after an initialization process of less than 5 min. Compared to the GPS-only fixed solutions, the real-time ZTDs with multi-GNSS PPP ambiguity resolution processing are significantly more stable with an obviously shorter initialization time. Utilizing the multi-GNSS PPP-AR processing, an improvement of 50.7% was achieved for initialization speed. Moreover, the RMS values for the ZTD deviation with respect to USNO and CODE final products decrease by 16.7% and 31.7%, respectively. Furthermore, we find from the comparisons with the ECMWF data that the accuracy of the GPS-only float solution is 14.3 mm while the average RMS is 12.5 mm for four-system fixed results, which indicates a great consistency between our GNSS-based ZTDs and ECMWF products. The new real-time tropospheric products with improved accuracy and reliability based on multi-GNSS PPP-AR will undoubtedly contribute to time-critical meteorological and geodetic applications.

Due to the availability of the four-system GNSS data, the employed stations in this study were mainly for Europe and the Asia-Pacific region. With the development of the global BDS3 system and the MGEX, more and more stations can provide high-quality four-system observations and contribute to GNSS meteorology. The observations of 30 days are sufficient for the validation of the effectiveness of the proposed approach, but further analysis (e.g., the seasonal variations in tropospheric delay) with a long period of several years will be carried out in the near future.