1. INTRODUCTION
The Precise Point Positioning (PPP) technique has been widely and successfully applied in precise orbiting, Global Positioning System (GPS) meteorology, precise timing, etc. (Zumberge et al., Reference Zumberge, Heflin, Jefferson, Watkins and Webb1997; Kouba, Reference Kouba2005). Conventional algorithms for PPP calculate the parameters of tropospheric zenith wet delay (ZWD), ambiguity, and receiver clock error together with the station coordinate (Keshin and Marel, Reference Keshin and van der Marel2008). Part of the remaining errors, including zenith hydrostatic delay (ZHD), antenna phase windup, tidal gravity and antenna phase centre offset are corrected using specific models, while the others such as satellite clock error are corrected with prior given values (Wu et al., Reference Wu, Wu, Hajj, Bertiger and Lichten1992; Kouba and Heroux, Reference Kouba and Héroux2001). Ionosphere delays are eliminated by ionosphere-free combination. It is crucial for the development of PPP to shorten the convergence time and refine its parameter estimation models (Bisnath and Gao, Reference Bisnath and Gao2008).
Many scholars have completed research to reduce the PPP convergence period. For conventional PPP using ionosphere-free measurements, the convergence takes about thirty minutes to get ten-centimetre-level accuracy (Bisnath and Gao Reference Bisnath and Gao2008; Zhang et al., Reference Zhang, Li and Guo2011). In recent years, some ambiguity-fixing methods have been proposed to improve PPP positioning accuracy (Ge et al., Reference Ge, Gendt, Rothacher, Shi and Liu2008; Collins et al., Reference Collins, Bisnath, Lahaye and Héroux2010; Loyer et al., Reference Loyer, Perosanz, Mercier, Capdeville and Marty2012; Li et al., Reference Li, Zhang and Ge2011), but they did not decrease the convergence time. Combining Near Real Time Kinematics (NRTK) and PPP, Teunissen et al. (Reference Teunissen, Odijk and Zhang2010) presented the PPP-RTK method based on the Continuous Operational Reference System (CORS) network, while Ge et al. (Reference Ge, Zou, Dick, Jiang, Wickert and Liu2010) and Zou et al. (Reference Zou, Ge, Tang, Shi and Liu2013) put forward an undifferenced network RTK method, which sped the convergence as well as assuring the accuracy. These two methods, which depend on the corrections broadcasted by the control centre, need a reference stations network. Cai and Gao (Reference Cai and Gao2007; Reference Cai and Gao2012) proposed the GPS/GLONASS PPP model distinct from the PPP model only using GPS observations. The GLONASS observations cut down the convergent period significantly, however, the positioning accuracy is not improved much when the GPS satellites are relatively abundant and the geometrical strength is favourable. Geng and Bock (Reference Geng and Bock2013) suggested a method using a triple-frequency GPS signal that can accelerate convergence and fix ambiguity. All of those studies provide valuable improvement for speeding up PPP, whereas their practical applications are still limited because of their reliance on the CORS network and of multi-GNSS and multi-frequency measurements.
This study decreases the convergent time through improving the algorithm model. Tropospheric delay is relatively significant among all errors that are corrected by models. The total tropospheric delay in the path of the GPS signal can be expressed to be the sum of the products of ZWD, ZHD and their respectively mapping function (Davis et al., Reference Davis, Herring, Shapiro, Rogers and Elgered1985) that is:
$$d_{trop} = M_d (elv_i^k ) \cdot ZDD + M_w (elv_i^k ) \cdot ZWD$$where M w(elv ik) and M d(elv ik) are the mapping function of ZWD and ZHD with the elevation angle of elv ik, respectively. The subscript i and superscript k, respectively, are the receiver and satellite; numerical values of both mapping functions closely approach each other. d trop is the total tropospheric delay along the slant path.
The conventional PPP algorithm deals with tropospheric delay using the strategy which corrects ZHD with the empirical Hopfield or Saastamoinen model and estimates ZWD as a parameter with random walk or piecewise linear process. However, the empirical models cannot obtain an accurate ZHD, and they need measured meteorological parameters above the station or approximations represented by standard meteorological parameters. Bar-sever et al. (Reference Bar-Sever, Kroger and Borjesson1998) improved the estimation accuracy of tropospheric delay and positioning by calculating tropospheric gradients, while two more parameters of gradients were added in this way, thus requiring at least seven satellites to obtain the full rank coefficient matrix.
Research on estimating zenith tropospheric delay (ZTD) based on the ZTD estimation model without meteorological elements has been well developed. Collins and Langley (Reference Collins and Langley1997) built the UNB model (a neutral atmosphere model developed by the University of New Brunswick) for the US Wide Area Augmentation System (WAAS), which has been developed into a series of models such as UNB1–4, UNB3 m, and UNB.na. The European Geostationary Navigation Overlay Service (EGNOS) ZTD estimation model has an accuracy generally comparable to the Hopfield and Saastamoinen model with measured meteorological parameters (Qu et al., Reference Qu, Zhu, Song and Ping2008). Li et al. (Reference Li, Yuan, Ou, Li and Li2012) constructed the latest global empirical model IGGtrop with analysed data from the National Centers for Environmental Prediction (NCEP) thereby improving the accuracy. While by using spherical harmonic functions and global tropospheric delay grids data offered by the Global Geodetic Observing System (GGOS), Yao et al. (Reference Yao, He, Zhang and Xu2013) established a global ZTD model without meteorological parameters. This global ZTD model, though having equal accuracy with IGGtrop, is more concise in model expression. All of these models get ZTD estimation functions by interpolating the global tropospheric delay data collected in a certain period and then estimate the parameters with a particular mathematical model. This paper will probe into their application in PPP.
Owing to the limited accuracy and spatial resolution of ZTD estimation models, the information derived from these models contains not only ZTD, but also error information, leading to a lower ZTD estimation model application efficiency. In order to overcome this drawback, we propose a new way to utilise the ZTD estimation model, namely the virtual observation method which treats ZTD calculated by the ZTD estimation model as a kind of observation instead of correcting it directly, combining it with phase and pseudo-range observations to build equations on the basis of conventional PPP. Since the virtual observations provide highly precise a priori initial values of the tropospheric delay, the separation of ambiguity and zenith coordinate parameters will also speed up. By better applying the model error correction to PPP, it is also effective in dealing with other error terms which are corrected by models, such as phase wind up, the tidal and even satellite clock bias.
To verify the validity of this new method, some experiments are carried out by using precise ephemeris, precise clock and observation files provided by the International GNSS Service (IGS) as well as comparing it with the conventional method, then verifying its accuracy and convergence time on a global scale.
2. POSITIONING MODEL
This section mainly introduces the two PPP models that are studied in this paper, namely the conventional model and the virtual observation model, presents their observation equations and analyses their advantages and disadvantages respectively.
GPS observations contain pseudo-range and phase, the basic observation equations refer to formulas given by Teunissen and Kleusberg (Reference Teunissen and Kleusberg1996). Assuming that the satellite and receiver antenna phase centre offset, the tidal relativistic effect, antenna phase wind up and zenith hydrostatic delay (calculated by the Saastamoinen model) have already been corrected in advance, the observation equations can be expressed as:
$$P_i^k = \rho _i^k + c(dt_i - dt^k ) + M_w (elv_i^k ) \cdot T_i^k + \varepsilon _i^k $$
$$L_i^k = \rho _i^k + c(dt_i - dt^k ) + M_w (elv_i^k ) \cdot T_i^k + \lambda b_i^k + e_i^k $$where c represents the speed of light; The subscript i and superscript k, respectively, are the receiver and satellite; L and P are the ionosphere-free combination of phase and pseudo-range measurements; ε and e denote the corresponding observation noise; ρ ik is the geometric distance between the satellite and receiver; T ik is ZWD, M w (elv ik) is the wet projection function when the elevation angle is elv ik. Receiver and satellite clock errors are dt i and dtk, respectively. b ik is phase ambiguity; λ is the wavelength of the corresponding non-ionosphere combination; Uncalibrated phase delay of the phase observations is absorbed by the ambiguity, while other errors, such as multi-path effect and pseudo-range hardware delay, are absorbed by the receiver clocks, ionospheric delays, and ambiguities.
2.1. Conventional model
The conventional model corrects the ZHD in advance by empirical model (here we use the Saastamoinen model) and estimates the ZWD parameter together with other parameters. Linearizing the basic Equations (2) and (3), the observation equation of conventional model can be expressed as:
$$\left[{\matrix{ {{\bi P}_{\bi i}^{\bi 1}} \cr \vdots \cr {{\bi P}_{\bi i}^{\bi m}} \cr {{\bi L}_{\bi i}^{\bi 1}} \cr \vdots \cr {{\bi L}_{\bi i}^{\bi m}} \cr}} \right] = \left[ {\matrix{ {{\bi - \mu} _{\bi i}} & {\bi I} & {\bi D} & {{\bi M}_{\bi w}} \cr}} \right]\left[ {\matrix{ {{\bi X}_{\bi i}} \cr {cdt_i} \cr {{\bi b}_{\bi i}} \cr {zwd_{\bi i}} \cr}} \right] + \left[ {\matrix{ {{\bi \varepsilon} _{\bi i}} \cr {{\bi e}_{\bi i}} \cr}} \right], {\bi P}_{\bi L} $$where μ i is a matrix composed of the unit vectors between the receiver and the satellite, corresponding to the station coordinates vector x i, I is a matrix of 2*m rows and one column, of which each element is one, corresponding to the clock parameter cdt i; Matrix D corresponds to the ambiguity parameters; M w is a matrix composed of wet projection function and corresponds to parameter zwd i; P L is the weight matrix of pseudo-range and phase observations. Satellite clock errors and satellite coordinates vector have been corrected by precise ephemeris and clock error products provided by the IGS.
This paper employs a Kalman filter to realize the algorithm. The ambiguity parameter remains constant when there is no cycle slip and the coordinate parameters of the state are determined according to the station movement. The IGS stations are relatively stable and their movement has no effect on our experiment. Receiver clock error is considered as white noise while the ZWD parameter is estimated by random walk process.
2.2. Virtual Observation model
We propose a new Virtual Observation Model to utilise the ZTD estimation model in PPP. This treats the ZTD as a virtual observation, combining it with phase and pseudo-range to build the observation equation. In this case, the observation equation can be expressed as:
$$\left[{\matrix{ {{\bi P}_{\bi i}^{\bi 1}} \cr \vdots \cr {{\bi P}_{\bi i}^{\bi m}} \cr {{\bi L}_{\bi i}^{\bi 1}} \cr \vdots \cr {{\bi L}_{\bi i}^{\bi m}} \cr {\Delta ZTD_i} \cr}} \right] = \left[ {\matrix{ {{\bi - \mu} _{\bi i}} & {\bi I} & {\bi D} & {{\bi M}_{\bi w}} \cr 0 & 0 & 0 & 1 \cr}} \right]\left[ {\matrix{ {{\bi X}_{\bi i}} \cr {cdt_i} \cr {{\bi b}_{\bi i}} \cr {zwd_i} \cr}} \right] + \left[ {\matrix{ {{\bi \varepsilon} _{\bi i}} \cr {{\bi e}_{\bi i}} \cr {\omega _i} \cr}} \right],{\rm} {\bi P}_{\bi L}, {\rm} P_{ZTD} $$
$$\Delta ZTD_i = ZTD_i - ZDD_i $$where ZTD i is calculated by ZTD estimation model; ZDD i is ZHD calculated by the Saastamoinen model in conventional PPP; M w is a matrix composed of the wet projection function and corresponds to parameter zwd i; ωi is the noise of the ZTD estimation model; the weight of ZTD i is denoted by P ZTD, the accuracy of the ZTD estimation model can reach several centimetres and its weight is about ten times higher than the pseudo-range measurements. The other symbols are the same as in Equation (4).
The determination of the virtual observations' weight matrix has a significant impact on the effectiveness of the algorithm. The result will be identical with that of the conventional method if the weights are too small, thus the role of the ZTD estimation model cannot be fully performed. The accuracy of the ZTD estimation model that this paper adopts is centimetre-level, which is ten times higher than the accuracy of pseudo-range observation, so the P ZTD should be ten times larger than the pseudo-range's weight.
As can be seen from Equation (5), the parameters in the observation equation of the new method have not been reduced compared with the conventional PPP algorithm, however, because of the existence of virtual observation, it can still complete the positioning when the satellite number m=4 so redundant observations are relatively increased. In addition, it does not reduce the number of parameters to be estimated and will not weaken the stability of the solution. By introducing the ZTD estimation model as virtual observations, one can utilise it more comprehensively.
3. COMPUTING STRATEGY
In order to verify the effectiveness of the algorithm, data from 196 IGS stations on 2 June 2009 (Figure 1) and of station ALBH during 2009–2011 were selected to conduct the analysis in space and time. The data that contains gross errors is eliminated through TEQC translation, editing and quality checking (TEQC). Then the Passion software, which is developed by our lab, was employed to compute PPP by using the new method, and conventional method respectively. The adopted post-processed sp3 file and precise clock file are offered by IGS with time intervals of 15 minutes and 30 seconds respectively, while the observation file sampling interval is 30 seconds. The coordinate true value of the IGS station uses the values under the International Terrestrial Reference Frame (ITRF) framework post-processed by IGS. The global ZTD model, a global zenith tropospheric delay model without meteorological parameters constructed by Yao et al. (Reference Yao, He, Zhang and Xu2013) by using GGOS Atmosphere, which has a global mean bias of 0·2 cm and an accuracy of 3·7 cm, is utilised in this paper as the ZTD estimation model. Actually, the choice of ZTD estimation model (from global ZTD model, IGGtrop, etc.) is not important because we can set a suitable weight for their estimations. The mapping function, identical for the three methods, is the Niell Mapping Function (Niell, Reference Niell1996).
Figure 1. Global distribution of the selected 196 stations.
We used daily observations of the selected IGS stations that are 24 hours of 30 sec samples. Our PPP algorithm computes a solution for every epoch. In order to obtain a better solution for each day, we conduct a smoothing program defined as:
$$Bias = \displaystyle{1 \over n}\sum\limits_{i = 1}^n {\left( {X_i - X} \right)} $$where X i is the i-th epoch positioning result; X is the true value provided by the IGS Analysis Center; not all epochs of a day are involved in this smoothing because the accuracy of the first few epochs are as large as several decimetres, or even several metres. So symbol n in the formula is the epoch number involved in the smoothing progress.
For the convergence time, we set a search condition: from the beginning epoch to an epoch whose horizontal component bias is less than 10 cm and zenith component bias is less than 15 cm, if the average deviation of its next 20 consecutive epochs also satisfies this requirement, then the observation time from the first epoch to this epoch is the convergence time.
4. VALIDATION AND ANALYSIS
This section analyses the computed results, including both accuracy and convergence time, obtained by the two methods mentioned above.
4.1. Accuracy analysis
Figure 2 shows the single day differential bias of all stations calculated by the two methods described above respectively. As can be seen, the bias of the new method and conventional method is roughly identical with differential less than 1 mm, especially for the horizontal component. Most stations have a quite small differential between the two methods. Large differences mainly appear in high altitudes because of poor quality of the ZTD estimate model in these areas.
Figure 2. The Bias difference of two methods for all stations.
In the following part we mainly analyse the accuracy and convergence time of the up component. Taking full advantage of the ZTD provided by the ZTD estimation model and its accuracy information and estimating ZWD as parameters, the new method is comparable with the conventional method in accuracy and stability.
4.2. Effects of the number of visible satellites on the positioning results
The statistical information of bias indicates the accuracy of the new method has rarely been improved compared to the conventional method. As mentioned in the Introduction, the positioning can be accomplished with only four satellites, so theoretically it would excel the conventional method in positioning accuracy. The reason why the improvement is not so obvious lies in the excellent observing conditions and a large number of visible satellites, which narrows the difference in positioning results caused by geometric strength. In order to highlight the advantages of positioning with four satellites by using new method, we used TEQC to edit observation data of the ARLT station on 2 June 2009. By removing some satellites, in the first three observing hours, only four satellites are available for a long period of time. Figure 3 shows the results of the three-hour positioning. The above figure indicates the bias curve of North, East, and up component as the observation time increases. The horizontal axis represents the positioning time and the vertical axis represents the positioning bias of each epoch. Figure 3 illustrates the number of visible satellites in each corresponding epoch. As can be seen, in the situation where there are only four satellites as well as the tropospheric delay, the new method can also be carried out and smoothly converge to centimetre-level accuracy within the ZWD parameter being estimated. On the other hand, the accuracy of the up component is much lower than that of the horizontal component. But as the number of satellites increases, up component accuracy is improved, while the horizontal component accuracy has not been improved much, which suggests the decrease in the number of satellites affects the zenith component notably, yet the horizontal component is much less influenced. Compared with the conventional method demanding at least five satellites to locate, the new method shows a substantial improvement.
Figure 3. At the positioning initial period for the ARLT station with only four satellites.
4.3. Convergence time analysis
In the conventional method, the first epoch often has a large bias, mostly in the sub-metre or even metre level, especially in the zenith component. Compared with conventional method, due to the use of the ZTD estimation model, the new method's first epoch bias is small, generally several decimetres or even centimetres. At the initial stage of positioning, the introduction of tropospheric delay as virtual observations contributes to the separation of ambiguity and tropospheric delay, accelerating the convergence time of the up component. In short, the convergence time of new method has been improved.
Figure 4 shows the improvement in the convergence time of four IGS stations' up component coordinate parameters. The horizontal axis represents the positioning time, and the vertical axis represents the positioning bias. As can be seen, the initial positioning bias of the new method is smaller than that of the conventional method particularly in the first epoch. The new method obtains convergence faster than the conventional method. In addition, before the convergence, its bias changes more smoothly without any dramatic fluctuations. As observation time grows, especially after convergence, the positioning results of the new method and the conventional method tend to be consistent. Among the four stations, the BAKO station has been improved quite remarkably. At the beginning of the positioning, the new method can quickly converge to the centimetre level with no large bias fluctuations, but a large fluctuation of 40 cm appears in the conventional method, which slows down the convergence time. With a relatively small improvement, the NYA1 station has a large fluctuation. Despite that, the fluctuation reaches approximately 40 cm in the conventional method, while reduced by nearly a half, it is only 20 cm in the new method.
Figure 4. Comparison of four stations' zenith coordinates convergence time.
Figure 5 gives a contrast of ambiguity convergence time between the PRN 14 satellite and the PRN 18 satellite observed by BAKO station. The horizontal axis represents the positioning time, and the vertical axis represents the ambiguity of each epoch. The ambiguity parameter here is the floating solution, and there is a slight difference between the calculated ambiguities of the two methods. During the converging period of the new method, the ambiguity parameters change flatly and converge fast. Like the coordinate parameters of the zenith component, the new method becomes consistent with the conventional method as the observation time increases. This fully shows that the new method can separate ambiguity parameters and other parameters faster and converge quickly, which is significant for fast ambiguity fixing.
Figure 5. Comparison of portion satellites ambiguity parameters convergence time comparison at BAKO station.
Table 1 and Figure 6 show distributions of the 196 global stations' convergence time. The stations that converge in less than 20 minutes by the new method clearly outnumber those converging by the conventional method, while the number of those requiring one hour to converge is smaller than that of the conventional method. From the statistical data of Table 1, the average convergence time of the new method is four minutes faster than the conventional method, improved by 15·4%. For some stations whose effects are notable, the convergence time can be shortened by about ten minutes, which indicates the new method can effectively accelerate the convergence time of PPP resolution. Moreover, faster convergence time is a great help for fixing ambiguity and improving the accuracy of single-day solutions.
Figure 6. The convergence time statistics of 196 global stations.
Table 1. Convergence time statistics.
5. DISCUSSION AND SUMMARY
The zenith component coordinates parameter, ambiguity parameter, ZTD and receiver clock bias parameter are highly correlated in the conventional PPP algorithm. The key factor that limits the convergence time is the isolation and fixing of the above three parameters. At the beginning of the positioning, because of the virtual observations, the new method can separate ZTD from the other three parameters earlier and fix it quickly, thus fastening the convergence of parameters in PPP, especially the zenith component coordinates and ambiguity parameters. Also, the method can achieve a relatively high positioning accuracy (generally in the decimetre-level) in the first epoch, being one order of magnitude higher than that of the conventional method. On the other hand, with a suitable weight matrix introduced to the virtual observation, it can analyse the model residuals, rendering PPP less dependent on the ZTD estimation model and avoiding extremely adverse impacts on the positioning results when a large deviation occurs in the ZTD estimation model. As the observation time grows, especially when the parameters to be estimated have been converged, the subsequent epoch calculation results will be consistent with the conventional method, unless there is a larger cycle slip or data interrupt requiring re-convergence, as the virtual observation's effect on positioning is much weaker after the tropospheric parameters have been essentially fixed.
This paper mainly addresses how to efficiently apply the ZTD estimation model to the PPP algorithm so as to accelerate the convergence time. First we analysed some problems existing in the conventional method, then we proposed a new method to improve the convergence of PPP. After comparing the pros and cons of the two methods in the accuracy and convergence time and verifying them by experiment, we reached the following conclusions:
• On accuracy, the new method is roughly comparable to the conventional method.
• If the satellites are few, especially when there are only four available, the new method excels the conventional method in its positioning accuracy. The new method has greater robustness for PPP positioning under poor observing conditions.
• Compared to the conventional method, the new method speeds the convergence time by 15·4% worldwide in the up component and converges more stations in less than 20 minutes with higher stability. From the perspective of combined accuracy and convergence time, the new method does use the ZTD estimation model more efficiently, therefore raising the convergence of the PPP resolution.
In the conventional PPP model, various errors are corrected by the direct use of the model, such as antenna phase winding error and the tidal correction, while others are corrected directly by the information that is already known, such as satellite orbit and clock error. These errors can be seen being corrected by directly using the model. Undoubtedly, the solution is not perfect because model errors are inevitable. So the new method proposed in this paper can be used to improve the solution, while the relevant effects need to be further studied.
ACKNOWLEDGEMENTS
The authors would like to thank the IGS Data Center for related data. This research was supported by the National Natural Science Foundation of China (41174012 and 41274022) and The National High Technology Research and Development Program of China (2013AA122502).
