1. INTRODUCTION
Today, Global Navigation Satellite System (GNSS) Precise Point Positioning (PPP) uses available GNSS orbit and clock offset correction products to perform point positioning by a single GNSS receiver. The PPP technique requires a number of accurate correction models to obtain a high-precision solution with undifferenced observations (Jin et al., Reference Jin, Wang, Zhang and Zhu2004). Tropospheric delay is one of the most significant error sources in GNSS positioning, navigation and geodetic applications (Jin et al., Reference Jin, Luo and Gleason2009; Reference Jin, Luo and Ren2010; Reference Jin, Han and Cho2011; Jin and Park, Reference Jin and Park2006). GNSS signals are affected by the Earth's atmosphere as they pass through the troposphere, causing a propagation delay of over two metres in the zenith direction (Jin et al., Reference Jin, Li and Cho2008). In precise positioning technology, the tropospheric delay is usually treated as Zenith Hydrostatic Delay (ZHD) and Zenith Wet Delay (ZWD) (Kouba, Reference Kouba2009); the latter is estimated directly as an unknown parameter, but the former is corrected by a non-meteorological model or a meteorological model in data processing. Atmospheric parameters such as pressure, temperature and water vapour pressure are needed for the meteorological model to calculate the Zenith Tropospheric Delay (ZTD) as in the Hopfield (Reference Hopfield1969) or Saastamoinen (Reference Saastamoinen1972) models, etc. Atmospheric parameters can be acquired from meteorological or some partial non-meteorological models at intervals of one day. Many GNSS stations are beginning to provide real observed meteorological data with a sampling rate of 30 s, which can explicitly show the variation of meteorological parameters and help study the impact of real measured atmospheric data on PPP.
Many scholars have assessed the performances of different tropospheric models. Leandro et al. (Reference Leandro, Santos and Langley2006) concluded that the prediction errors of the University of New Brunswick (UNB3m) model have a mean bias of −0·5 cm and standard deviation of 4·9 cm based on ray-tracing analyses of 703,711 profiles from 223 stations in North America and surrounding territories from 1990 to 1996. Zhang et al. (Reference Zhang, Yuan, Li, Li and Chai2016) conducted a comprehensive study of the performances of a global three-dimensional grid-based ZTD model called IGGtrop, and European Geo-stationary Navigation Overlay System (EGNOS) and UNB3m models in China where the IGGtrop model performed the best. Hadas et al. (Reference Hadas, Kaplon, Bosy, Sierny and Wilgan2013) applied various ZTD models to PPP using Bernese Global Positioning System (GPS) software and assessed the quality of the models. Ding et al. (Reference Ding, Teferle, Kazmierski, Laurichesse and Yuan2017) estimated the tropospheric delay and analysed the initialisation period and accuracy of the GNSS troposphere estimates. Yao et al. (Reference Yao, He, Zhang and Xu2013) analysed the temporal-spatial variations of Global Zenith Troposphere Delay (GZTD) using the time series of Four-Dimensional (4D)-grid ZTD from 2002 to 2009 provided by the Global Geodetic Observing System (GGOS) compared to a series of University of New Brunswick (UNB) models.
There are numerous models that can calculate meteorological data including the standard atmospheric model, UNB3m model, Global Pressure and Temperature (GPT) model and Global Pressure and Temperature-2 (GPT2) model (Leandro et al., Reference Leandro, Santos and Langley2006; Reference Leandro, Langley and Santos2008; Böhm et al., Reference Böhm, Heinkelmann and Schuh2007; Kouba, Reference Kouba2009; Lagler et al., Reference Lagler, Schindelegger, Böhm, Krásná and Nilsson2013). Among them, the accuracy of the standard atmospheric model is the worst. Accordingly, this is also reflected in the PPP accuracy. The results of other meteorological models make little difference on PPP. Steigenberger et al. (Reference Steigenberger, Boehm and Tesmer2009) compared station coordinates with a global distribution obtained by different troposphere models, which showed that the site height differences computed with the GPT model and the European Centre for Medium-Range Weather Forecasts (ECMWF) numerical weather model data are below 1 mm in general, and the horizontal differences are even smaller. Real meteorological data can be obtained after daily observation collected by meteorological sensors and it can be applied for post-processing PPP. Precise hydrostatic delay calculated from meteorological data with high reliability and accuracy can improve the performance of PPP in satisfying the need for high precision navigation and positioning.
This paper addresses issues related to meteorological parameters and its influence on PPP. In Section 2, the different meteorological models including the standard atmospheric model, UNB3m, GPT, and GPT2 are introduced. In Section 3, the multi-GNSS PPP processing strategy is proposed. In Section 4, experiments of different GNSS PPP with real meteorological data and the meteorological model to correct the ZHD are implemented to analyse the improvement of multi-GNSS PPP with real observed meteorological data. Finally, some conclusions and perspectives are given in Section 5.
2. METEOROLOGICAL MODELS
2.1. Standard atmospheric model
The height-dependent values of pressure (P), temperature (T), and relative humidity (R h) derived from a standard atmospheric model may be obtained by the equations:
where P 0, T 0 and R h0 are standard pressure, temperature and humidity at the reference height H 0 (H 0=0 m, P 0=1013·25 mbar, T 0=18°C, R h0=50%) and H represents the height of the station.
The partial pressure of water vapor e (in mb) can be obtained by the equation:
The meteorological parameters derived from the standard atmospheric model have low accuracy, which cannot reflect the temporal-spatial variations of atmospheric parameters (Kouba, Reference Kouba2009).
2.2. UNB3m model
Trenberth (Reference Trenberth1981) derived the UNB2 model based on the analysis of the Saastamoinen (Reference Saastamoinen1972) model considering the variation of the water vapour profile and latitude. With the latitude of the site according to the model, it is possible to calculate the average value of the meteorological parameters, including pressure, temperature, relative humidity, temperature lapse rate and water vapour pressure, but their temporal variation cannot be reflected (Swanson and Trenberth, Reference Swanson and Trenberth1981). Collins and Langley (Reference Collins and Langley1997) derived the UNB3 model using the standard meteorological datum to acquire the average value and yearly amplitude of five parameters. The user can obtain the value of meteorological parameters according to the latitude and epoch of the site to acquire the tropospheric delay. However, the part of the relative humidity converted from the water vapor pressure in the UNB3 model will exceed 100% (Collins and Langley, Reference Collins and Langley1999). Hence, Leandro et al. (Reference Leandro, Langley and Santos2008) refined the UNB3m model, which can calculate the relative humidity between 75% and 85%. In the UNB3m model, five output parameters including pressure, temperature, relative humidity, temperature lapse rate and water vapour pressure height factor can be acquired from the position of the site and the day of year. The meteorological parameter values for a particular latitude and day of the year can be obtained using a look-up table (Leandro et al., Reference Leandro, Langley and Santos2008). The annual average and amplitude of given parameter can be computed as:
where ϕ is the latitude of the station of interest, Avg ϕ is the computed average, Amp ϕ stands for the computed amplitude, i is the index of the nearest lower tabled latitude and Lat represents latitude from the look-up table.
After the average and amplitude are computed at the given latitude, the parameter values P 0, T 0, RH 0, β0 and λ0 can be estimated in the desired day of year according to:
where X ϕ,doy represents each computed parameter value for latitude ϕ and day of year doy. The conversion between relative humidity and water vapour pressure e 0 can be carried out as follows:
where the saturation vapour pressure e s can be computed as:
The enhancement factor f w can be determined as follows:
Hence, the meteorological parameters can be represented as:
where the power value of water vapour pressure e p=g/R/β 0, g is the surface acceleration of gravity (9·80665 ms−2) (Tenzer et al., Reference Tenzer, Chen, Tsoulis, Bagherbandi, Sjoberg and Novak2015), R is the gas constant for dry air (287·054 J kg−1 K−1), λ′=1+λ and H is the orthometric height in m.
2.3. GPT model
In order to express the variation of the meteorological information in different places and times, Böhm et al. (Reference Böhm, Heinkelmann and Schuh2007) proposed the GPT model based on three years (September 1999 to August 2002) of 15°×15° global grids of monthly mean profiles for pressure and temperature from ECMWF reanalysis data. The corresponding grid of orthometric heights of the Earth surface with respect to mean sea level is also available from the ECMWF. The parameters computed by the GPT model according to the station coordinate and day of the year can be expressed as:
where P 0 and T 0 are the mean value of pressure and temperature at mean sea level and doy is the day of year. The mean values a 0 and the annual amplitudes A expressed by degree and order nine are available:
where P nm is the Legendre polynomials, φ and λ are latitude and longitude and A nm and B nm are the coefficients for degree n and order m. Hence, the meteorological parameter can be acquired by Equation (14).
2.4. GPT2 model
The GPT2 model is based on the global monthly mean profile of pressure, temperature, specific humidity and geopotential from ECMWF reanalysis of data from 2001 to 2010 (Lagler et al., Reference Lagler, Schindelegger, Böhm, Krásná and Nilsson2013). It can provide the pressure P 0, temperature T 0, lapse rate dT and water vapour pressure e 0 at 5°×5° grid intervals. With regard to the changing meteorological parameters, the mean values A 0, annual values (A 1, B 1) and semi-annual values (A 2, B 2) various calculations expressed by parameter r(t) can be carried out as follows:
The inputs of the model are the coordinates of the site and epoch. The correction of the pressure P, temperature T and water vapour pressure e can be carried out as follows:
where Q is the specific humidity, gm is the gravity (gm=9·80665 m/s2), dM r is the molar mass of dry air (dM r=0·0028965 kg/mol) and R g is the universal gas constant (R g=8·3143 J/K/mol). H′ is the orthometric height. The parameter c can be represented as:
3. PPP PROCESSING STRATEGY
The traditional model of PPP can be represented as (Zumberge et al., Reference Zumberge, Heflin, Jefferson, Watkins and Webb1997):
where P r,ifs is the Ionosphere-Free (IF) combination of pseudorange P 1 and P 2 observations at two distinct signal frequencies f 1 and f 2, $\Phi_{r,IF}^{s}$ is the IF combination of the corresponding carrier-phases Φ1 and Φ2 and ρrs is the geometrical range from the satellite position at the signal emission epoch to the receiver position at its arrival epoch. δ t r is the receiver clock offset, δ t s is the satellite clock offset from the GNSS system time and c is the vacuum speed of light. T rs is the troposphere delay, A IF is the non-integer ambiguity of the IF carrier phase combination, λIF is the IF combination of the carrier phase wavelengths and e IF and εIF are the unmodelled multipath error and the relevant measurement noise components.
The observations are weighted with the satellite elevation angle. The stochastic model can be represented by a sinusoidal function (Witchayangkoon, Reference Witchayangkoon2000):
where E is the satellite elevation angle and σ0 is the prior variance of the observations. When the observations of the pseusorange and carrier phase are used at the same time, the variance-covariance expression is:
The experiment is carried out by using the multi-GNSS software developed independently. When the ionosphere effect is mitigated through the IF combination, the unknown parameters of the traditional PPP model are receiver position coordinates, receiver clock, ZHD and IF carrier phase ambiguities (non-integer). The GPS and GLONASS satellite Phase Centre Offset (PCO) and Phase Centre Variations (PCV) are corrected by the ‘igs08_1861.atx’ file generated by IGS. The BDS and Galileo antenna offset and variation recommended by Multi-GNSS Experiments (MGEX) are used to correct the PCO and PCV of BDS and Galileo satellites. The receiver PCO is derived from the components provided by the antenna file (Rizos et al., Reference Rizos, Montenbruck, Weber, Weber, Neilan and Hugentobler2013). When integrated multi-GNSS PPP is carried out, additional system time difference is needed for each newly added GNSS system. The cut-off angle is set to 5°. A Kalman filter is used in the process of parameter estimation. The GPS and GLONASS code observation precision is set to 0·3 m and the phase observation precision is set to 0·003 m. The BDS and Galileo code observation precision is set to 0·6 m and the phase observation precision is set to 0·004 m (Cai et al., Reference Cai, Gao, Pan and Zhu2015). The values of parameters are determined by the precision of observations. The parameter estimation strategy in PPP is shown in Table 1. In data processing, the meteorological models listed in Section 2 and the real observed surface atmospheric parameters provided by the International GNSS Service (IGS) are used to provide the atmospheric parameters required for the improved Hopfield tropospheric model (Hopfield, Reference Hopfield1969), to compute the ZHD, and to evaluate the impact on PPP. The ZWD is represented by a random-walk process and the receiver clock offset is estimated epoch-by-epoch. The power density of the tropospheric and receiver clock noise is set to 5 mm hr−1/2 and 100 m s−1/2, respectively (Hadas et al. Reference Hadas, Teferle, Kazmierski, Hordyniec and Bosy2017). The modified Hopfield model for calculating ZHD can be summarised as:
where z is the zenith angle of the satellite, T is the temperature at the station (in units of Kelvin [K]), P is the atmospheric pressure (in units of millibars, mb), e is the partial pressure of water vapour (in mb), R E is the Earth's radius, ZTD is the zenith tropospheric delay, ZD dry is the zenith tropospheric dry delay, ZD wet is the zenith tropospheric wet delay, N is atmospheric refractivity for dry or wet component, r is the distance from the centre of the Earth to the station, h is the height and i represents dry or wet part of the delay. The convergence criterion is defined when the component of positioning coordinates is less than 0·1 m and keeping within 0·1 m in the subsequent 20 epochs. The positioning errors are calculated from the convergence epoch to the end epoch of a day.
4. RESULTS AND ANALYSIS
4.1. Data description
The datasets were collected at eight IGS stations for the period of 1 to 15 May 2015 (day of year from 121 to 135). Figure 1 shows the distribution of the selected stations. All of the stations can receive quad-system observations from GPS, BDS, GLONASS and Galileo constellations. The processing strategy to obtain a PPP solution is listed in Table 1.
The observation data of XMIS station on 1 May 1 2015 was analysed as an example in the following subchapters. The XMIS station is located on Christmas Island, Australia. The type of antenna was JAVRINGANT_DM+ NONE. The data sampling rate of was 30 s. Figure 2 shows the partially observed satellite trajectory of BDS, GPS, GLONASS and Galileo on 1 May 1 2015. It shows that the satellite orbits of 13 BDS satellites (C01~C14, except C13), 32 GPS satellites and 24 GLONASS satellites can be tracked. All BDS satellites, GPS satellites (G01~G07), GLONASS satellites (R01~R07) and Galileo satellites (E11, E12 and E19) are shown in the figure. The integration of different GNSS can significantly increase the number of visible satellites and obtain better zenith coverage. The BDS Geosynchronous Earth Orbit (GEO) satellites had the longest tracking periods, for example, C01~ C04 can be observed all day. The BDS Inclined Geosynchronous Orbit (IGSO) satellites C06~ C10 were tracked for a shorter period than GEO satellites, the Medium Earth Orbit (MEO) satellites were tracked for the shortest period. The number of visible satellites and the Position Dilution Of Precision (PDOP) of each combination (BDS, GPS, GLONASS, BDS+GPS, and four-system) at XMIS on 1 May 2015, are shown in Figure 3. At each epoch, The PDOP values for BDS, GPS and GLONASS are 0·5~3, 1·5~4 and 1~ 6 respectively. The improvement of the spatial geometry is obvious under integrated multi-GNSS. It can be seen that the PDOP values in four-system PPP are below 1·0 in this period.
4.2. Performance of different PPP solutions
Figure 4 shows the results of different GNSS PPP solutions with real meteorological data-calculated ZHD at the XMIS station. At XMIS, the coordinate accuracy of BDS-only PPP was 5·7 mm, 12·9 mm and −14·6 mm in the north, east and up components while for GPS-only PPP the values were 2·9 mm, 4·3 mm and −5·6 mm, which is similar to BDS+GPS PPP, the position errors in three directions were 2·8 mm, 4·2 mm and −5·4 mm. For GLONASS-only PPP, the final positioning accuracy values were 4·3 mm, 5·5 mm and 9·1 mm. The positioning accuracy for different GNSS PPP can reach the level of several millimetres to centimetres. The statistical convergence time of BDS-only, GPS-only, GLONASS-only and BDS+GPS was 53·0 min, 32·0 min, 36·5 min and 25·5 min, respectively.
In order to further evaluate the influence of the ZHD correction calculated by surface observed atmospheric parameters on BDS-only PPP, GPS-only PPP, GLONASS-only PPP and BDS+GPS+GLONASS+Galileo PPP, datasets collected at eight stations for 15 days were used. Figure 5 shows the average Root Mean Square (RMS) of position errors in the north, east and up components at the eight selected MGEX stations. Figure 6 gives the average convergence time of the corresponding stations.
According to Figures 5 and 6, we can draw the conclusion that the positioning accuracy of BDS-only PPP using surface observed meteorological data for ZHD corrections is the worst, which is mainly limited by the current distribution and number of BDS satellites and the orbit determination accuracy of the BDS constellation. Except for MRO1 and NNOR stations, the RMS of position errors for BDS-only PPP in three directions is below 2·5 cm and the GPS-only PPP accuracy is slightly better than GLONASS-only PPP. The positioning accuracy of GPS-only PPP, BDS+GPS PPP, and BDS+GPS+GLONASS+Galileo PPP are approximately the same with observed meteorological data and the coordinate accuracy is better than 5 mm, 5 mm and 6 mm in the north, east and up components. Hence, integrated multi-GNSS PPP accuracy is not significantly improved in comparison with GPS-only PPP.
We can also conclude that the convergence time of BDS-only PPP is also the longest of the single systems and GPS-only is the best. Combination systems are obviously better than a single system. The average convergence time for different types of PPP is 55·89 min, 25·88 min, 33·30 min, 20·50 min and 15·71 min for BDS-only PPP, GPS-only PPP, GLONASS-only PPP and four systems PPP, respectively. Integrated quad-system PPP can significantly accelerate convergence, which is one of the toughest problems in current PPP solutions.
4.3. Performance of PPP with different atmospheric parameters
In order to evaluate the performance of improved Hopfield ZHD correction with surface observed meteorological data in GNSS PPP, the results of BDS+GPS+GLONASS+Galileo quad-system PPP are compared with the atmospheric parameters provided by the standard atmospheric model, UNB3m, GPT and GPT2 for ZHD correction in PPP. The meteorological parameters computed from four models at the XMIS station on 1 May 2015 are listed in Table 2, and Figure 7 gives the time series of surface observed meteorological parameters provided by IGS, whose sampling rate is 30 s.
From Table 2, it can be seen that the atmospheric parameters derived from meteorological models are the daily average. Hence, the model outputs cannot reflect the sub-diurnal variation characteristics of atmospheric parameters. There are also significant deviations between the standard atmospheric model and the other models. Taking the surface measured data as the reference, the accuracy of water vapour pressure provided by GPT2 is slightly better than that from UNB3m and GPT. The performances of integrated multi-PPP are examined and the meteorological parameters needed in the improved Hopfield model for ZHD correction provided by observed data and meteorological models are evaluated in a PPP solution. Other model corrections for PPP are listed in Table 1. Figure 8 shows the position error series of BDS+GPS+GLONASS+Galileo PPP in the north, east and up components. Where StandardModel stands for the results using the standard atmospheric model, UNB3m represents that the UNB3m model is adopted, GPT stands for the PPP strategy with GPT model outputs for ZHD correction and GPT2 and MetData stand for the results using the GPT2 model and observed meteorological data.
The positioning results of the BDS+GPS+GLONASS+Galileo PPP with ZHD correction using the atmospheric parameters provided by real meteorological data and the output values from the standard atmospheric model, UNB3m, GPT and GPT2 are analysed. Table 3 and Figure 9 show the RMS values of position errors in the north, east and up components at eight stations. Figure 10 shows the average convergence time in three directions at corresponding stations. Figure 11 gives the statistical results of position accuracy and convergence time with different atmospheric parameter sources for ZHD correction.
By comparing Table 3, Figure 9 and Figure 11, it can be seen that the RMS of positioning error in the up component is obviously bigger than that in the north and east components. When the atmospheric parameters provided by the standard atmospheric model for computing ZHD corrections in muti-GNSS PPP is used, the RMS of positioning error in the up component is ten times bigger than the horizontal components. The result of PPP using real atmospheric data is the best and the RMS of position error for the other three models is almost the same. The statistical coordinate differences in the north and east components for all of the atmospheric parameter sources are lower than 1·4 mm. Hence, the ZHD model with accurate atmospheric parameters benefits the position accuracy in the up component. Compared with the observed meteorological parameters, the RMS of position error with the model parameters from the standard atmospheric model, UNB3m, GPT, and GPT2 can be decreased by 90·6%, 33·0%, 22·2% and 19·8%, respectively. The convergence time of integrated multi-GNSS PPP in the horizontal directions are almost the same, as shown in Figures 10 and 11. The convergence time for the four meteorological models in the up component is 28·48 min, 20·68 min, 20·59 min and 20·53 min, respectively. It only takes 13·90 min for surface observed atmospheric parameters to converge and the convergence efficiency is improved by 51·2%, 32·8%, 32·5% and 32·3%, respectively.
5. CONCLUSION
The aim of this paper is to study the performance of PPP using real meteorological data for ZHD correction. Firstly, this contribution summarises the standard atmospheric model, UNB3m, GPT and GPT2 meteorological models, which are usually adopted for providing atmospheric parameters for calculating the ZHD correction. Then the multi-GNSS PPP experiments are carried out to evaluate the impact of ZHD correction with different sources of atmospheric parameters on the positioning solution.
The accuracy of different GNSS PPP solutions using real atmospheric parameters for ZHD correction can reach several millimetres. The statistical results show that the BDS-only PPP accuracy is the worst for a single GNSS system, the average RMS of positioning errors is 13·6 mm, 13·5 mm and 20·5 mm in the north, east and up components while the ones for GPS-only PPP are the best. BDS+GPS PPP and BDS+GPS+GLONASS+Galileo PPP can slightly improve the positioning accuracy, and the RMS of positioning errors in three directions is lower than 5 mm, 5 mm and 6 mm. According to our convergence criterion, the average convergence time for BDS-only PPP, GPS-only PPP, GLONASS-only PPP, BDS+GPS PPP and BDS+GPS+GLONASS+Galileo PPP is 55·89 min, 25·88 min, 33·30 min, 20·50 min and 15·71 min, respectively; integrated multi-GNSS PPP significantly contributes to improving the convergence time, which is of high value in improving the utility of PPP.
Compared with atmospheric parameters provided by meteorological models for ZHD correction computing in PPP, the improvement of positioning accuracy with real atmospheric data is not obvious in horizontal components. As we know that there is a strong link between the coordinate solution in the up component and tropospheric delay, our experimental results confirm that using externally observed meteorological parameters for troposphere hydrostatic correction in PPP can greatly benefit the accuracy of coordinates in the up component, compared to the parameters provided by the standard atmospheric model, UNB3m, GPT and GPT2 computing ZHD for a PPP solution. The vertical RMS of positioning errors can be improved by 90·6%, 33·0%, 22·2% and 19·8% respectively. The convergence time in the up component is decreased by 51·2%, 32·8%, 32·5% and 32·3%, respectively.
ACKNOWLEDGEMENTS
This research is supported by the National Natural Science Foundation of China (NSFC) Project (Grant No. 11573052). The authors gratefully acknowledged the GFZ for providing products and IGS for the data from MGEX.