Precise Point Positioning on the Reliable Detection of Tropospheric Model Errors

Precise point positioning (PPP) is one of the well-known applications of Global Navigation Satellite System (GNSS) and provides precise positioning solutions using accurate satellite orbit and clock products. The tropospheric delay due to the neutral atmosphere for microwave signals is one of the main sources of measurement error in PPP. As one component of this delay, the hydrostatic delay is usually compensated by using an empirical correction model. However, how to eliminate the effects of the wet delay during a weather event is a challenge because current troposphere models are not capable of considering the complex atmosphere around the receiver during situations such as typhoons, storms, heavy rainfall, et cetera. Thus, how positioning results can be improved if the residual wet delays are taken into account needs to be investigated . In this contribution, a real-time procedure of recursive detection, identification and adaptation (DIA) is applied to detect the model errors which have the same effects on both phase and code observables; e.g., the model error caused by the tropospheric delay. Once the model errors are identified, additional parameters are added to the functional model to account for the measurement residuals. This approach is evaluated with Global Positioning System (GPS) data during two rainfall events in Darwin, Australia, proving the usefulness of compensated residual slant wet delay for positioning results. Comparisons with the standard approach show that the precision of the up component is improved significantly during the periods of the weather events; for the two case studies, 72.46% and 64.41% improvements of root mean squared error (RMS) resulted, and the precision of the horizontal component obtained by the proposed approach is also improved more than 30% compared to the standard approach. The results also show that the identified model errors are concentrated at the beginning of both heavy rainfall processes when the front causes significant spatial and temporal gradients of the integrated water vapor above the receiver.


Introduction
The troposphere is the lowest portion of the Earth's atmosphere, and tropospheric delay due to the neutral atmosphere is one of the main error sources of the Global Navigation Satellite System. This delay can cause up to 2.5 m at zenith direction of the Global Navigation Satellite System (GNSS) signal transmission and over 20 m when satellites are at low elevation angles; e.g., below 10 degree [1,2]. The tropospheric delay is commonly expressed with the following model [3] T(e, α) = M h (e) · Z h + M w (e) · Z w + M g (e) · cot(e) · (G N · cos(α) + G E · sin(α)) where e and α are respectively the elevation and the azimuth angle of a specific satellite. The total tropospheric slant delay T between receiver and satellite at an elevation angle e is the sum of three portions: a hydrostatic portion, a wet portion and a gradients portion. Z h and Z w are the zenith hydrostatic delay and zenith wet delay, respectively. M h and M w are the mapping functions for the zenith hydrostatic and wet delay, respectively. G N and G E are the gradients which account for the azimuthally inhomogeneous troposphere in north-south and east-west directions with the corresponding gradient mapping function M g . The hydrostatic delay due to the refractivity of the dry gases in the troposphere can be corrected by the conventional models such as Saastamoinen [1] and Hopfield [4], which can model the hydrostatic delay at the millimeter level in the zenith direction [5]. Collins and Langley [6] proposed a neutral atmosphere model designed for Wide Area Augmentation System (WAAS) users, which is the so-called UNB model series (UNB1 through UNB4) and has been assessed for the use in North America [7,8], Europe [9] and Japan [10]. Li et al. [11,12] developed a multi-dimensional grid model, IGGtrop, to provide tropospheric delay corrections for the users of the BeiDou Navigation Satellite System (BDS) and the area augmentation system based on BDS in China. Although the models mentioned above can correct the wet delay to some extent, the accuracy varies from centimeter to decimeter level, which is still insufficient for high precision positioning and navigation. In addition, using the empirical atmospheric information obtained from the profile of global pressure and temperature may reduce the accuracy of the troposphere models due to the high spatial and temporal variability of water vapor [13,14]. Therefore, the zenith wet delay is usually estimated as an unknown parameter at each epoch or within a certain time span.
When the zenith tropospheric delays are estimated or provided, the slant delays to each visible satellite are obtained by assuming a specific relation between the zenith and slant direction in which the troposphere is assumed to be symmetrical about the vertical direction of the receiver. The relation between zenith and slant delay can be modeled by a so-called mapping function, as already shown in Equation (1), for which a wide range of mapping functions have been developed in the past. The Niell mapping function (NMF) [15] and the global mapping function (GMF) [16] consist of easy-to-handle formulae which only need the input parameters of approximate latitude, height and day of year. On the other hand, the isobaric mapping function (IMF) [17] and the Vienna mapping function 1 (VMF1) [18] provide support for mapping functions derived from numerical weather models (NWM) by applying the ray-tracing technique and/or climatological data. The crucial variable in mapping functions is the elevation angle. Most mapping functions are azimuth-independent, which reveals the underlying assumption that the troposphere is azimuthally homogeneous. A successful application of an azimuthally inhomogeneous tropospheric delay modeling in GPS geodesy and very long baseline interferometry (VLBI) was proposed by MacMillan [19] and Chen and Herring [20], in which the so-called horizontal gradients are considered in addition to a mapping of the zenith to slant delays. In this way, a linear asymmetry of the troposphere is accounted for by introducing a tilted direction instead of the zenith direction. For an extensive review of the troposphere model and mapping function, see Teunissen and Montenbruck [21].
Positioning in severe weather conditions has received more attentions in recent years. Yasyukevich et al. [22] investigated the influence of solar flares on the GNSS and high-frequency propagation. Luo et al. [23] analyzed the performance of double and single-frequency base PPP during three typical geomagnetic storms. As for the tropospheric delay, the standard troposphere model is capable of estimating the tropospheric delay with centimeter accuracy in normal weather conditions [24,25]; however, it should be investigated how positioning results can be improved if the residuals of the tropospheric delay caused by weather events are taken into account. The issue is that satellites at the same elevation angle would be compensated by almost the same tropospheric delay correction based on the standard mapping function approach. However, the symmetrical troposphere about the zenith direction of the receiver is not realistic when it suffers from the complex weather situation. The performance of the horizontal gradients is also limited, because they can only consider a linear asymmetry of the troposphere around the geodetic site [26]. Kleijer [27] analyzed that significant biases can be introduced in the estimated ZWD when the atmosphere is not symmetrical. However, the suggestion of using an accurate wet mapping function is still limited by the assumption of the homogeneous troposphere. Li et al. [28] assessed the impacts of the tropospheric biases on the integer ambiguity resolution and gave the recommendations of under which conditions the tropospheric biases can be ignored. However, only zenith tropospheric biases are taken into account, without considering the biases caused by the inhomogeneous troposphere. Hobiger et al [29] proposed a method to combine the mesoscale and fine-mesh numerical weather model to provide the ray-traced tropospheric slant delay during a typhoon passage. The result shows that the height repeatability is improved up to 30% compared to standard data processing. However, this could still be insufficient for high-precision positioning, and it is not possible to provide the fine-mesh numerical weather model to worldwide users in (near) real-time.
The detection, identification and adaptation (DIA) procedure was first demonstrated by Baarda [30] and Teunissen [31,32]. Teunissen [33] introduced this method into GNSS to detect, identify and adapt the mismodeled errors, and then it was applied in a wide range of GNSS applications; for example, kinematic GNSS surveying [34], permanent station resolution [35] and observation quality control [36,37]. In this contribution, a real-time recursive DIA procedure is implemented to detect the model errors which have the same effects on both phase and code observables, and once the errors are identified, additional parameters will be added to the functional model to account for the measurement residuals. One of the applications of this approach is to detect model errors caused by the tropospheric delay; therefore it was evaluated with GPS data during two different rainfall events in Darwin, Australia, proving the usefulness of compensated residual slant tropospheric delay for positioning results. Comparisons with the standard approach show that the precision of the up component is improved significantly during the period of the weather events, and the precision of the horizontal component is also improved.
This article is organized as follows. Section 2 reviews the standard functional model for PPP data processing and the theory of DIA and the construction of the improved functional model, which takes into account the model errors. Section 3 analyzes the performance of the proposed procedure via two case studies during a weather event. Section 4 contains the summary and conclusions.

Modeling and Filtering
The undifferenced, ionosphere-free (IF), linear combinations of phase and code are used as the basic observables [38] ∆φ s r, where ∆φ s r,IF and ∆p s r,IF represent the so-called observed-minus-computed IF combinations for phase and code observable in meters, respectively. Notice that the a priori hydrostatic delay has been corrected for these observations; u s r denotes a unit line of sight vector from satellite s to receiver r; ∆r r contains the increments of geodetic coordinate; dt r refers to the receiver clock offset. The satellite clock offset has been corrected by a priori precise products. The wet tropospheric delay Z w is the main interest in this study. The notations of two horizontal gradients are M g,N := M g (e) · cot(e) · cos(α), M g,E := M g (e) · cot(e) · sin(α), and the definitions of e, α, G N and G E have been illustrated in Equation (1). λ IF denotes the wavelength of the IF combination and N s r,IF the IF ambiguity. Note that both receiver and satellite hardware delays have been ignored because they are not the main parameters of interest. ε s r and e s r are phase and code measurement errors, respectively. It is worth noting that this research only applies the traditional IF combination, and for a more rigorous model, one needs to consider a third observable; that is, the difference between the wide lane phase and the narrow lane pseudorange [39].
After collecting the observed-minus-computed observables ∆φ s r,j and ∆p s r,j at epoch k, the IF combination vector of the phase and code y k can be formed from dual-frequency observations. One can symbolize a linear model of the compact formula of Equation (2) as where A k is the so-called design matrix and x k the n−dimensional state vector containing the unknown estimable parameters; e k refers to the measurement noise vector with e k ∼ N(0, R k ). The linear dynamic model describing the time evolution of the unknown parameters is given as where x k and x k−1 refer to the state vectors of the system at epochs k and k − 1, respectively; Φ k,k−1 represents the transition matrix between two epochs. This matrix is regarded as the identity matrix, because the dynamic system is described by the differential equations of a first-order linearized system, and the identity matrix is obtained by solving the first-order vectorial differential equations. d k represents the system noise at epoch k with d k ∼ N(0, Q k ) and is assumed to be uncorrelated in time.
The initial state of the system and its variance matrix can be given aŝ The time update state vector and its variance matrix are given aŝ wherex k|k−1 is the predication of the unknown parameters at epoch k, and P k|k−1 its corresponding predicted variance matrix. The predicted residual vector and its variance matrix can be given as Using the predicted residual vector, the updated state and its variance matrix are given aŝ With the gain matrix

Detection
The objective of the detection step is to detect the mismodeling errors of the mathematical model. The functional model of Equation (2) will be tested at each epoch k to detect the presence of model errors; e.g., unmodeled outliers in one or more observations. The distributional property of the predicted error of the unbiased functional model can be expressed as Equation (10) [40,41], which is the so-called null hypothesis model.
Otherwise, if the functional model is biased, the distributional property then turns to where C v k is an m × p matrix with p additional unknown parameters and p−vector ∇ is assumed to be unknown. Equation (11) is the so-called alternative hypothesis model. The test statistic for detecting model errors reads as where r k is the redundancy at epoch k. Equation (12) is also referred to as the local overall model (LOM) test, and model errors are considered to be present at epoch k if where F α (r k , ∞, 0) is the critical value based on the central F−distribution with the level of significance α and two degrees of freedom r k , ∞.

Identification
This step is to identify the most likely model error to account for the unexpected effects. For simplification, the case of one single model error at an epoch is considered for each recursion of the DIA process, so the matrix C v k of Equation (11) reduces to a vector, and the test statistic reads as Commonly, this test is applied to test for outliers in a single observation only. In order to test for an outlier in the ith observation, the c i -vector should be a 1 as its ith element, and zero otherwise. In this study, the goal is to also consider a model error associated with a residual tropospheric delay for one of the satellites due to asymmetry of the troposphere, which implies that the model error has the same influence on both phase and code observables. Thereby, the vector c T i turns to c T i = [0, ..., 1, ..., 0, ..., 1, ..., 0] which denotes all the elements in this vector are 0 except for the two 1s corresponding to the ionosphere free combined phase and code observation for one particular satellite. As for the uncombined dual-frequency observable, the c vector can extend to four 1s to account for the same model error on two phase observations and two code observations. v k is the predicted residual vector and Q v k v k its corresponding variance matrix; see Equation (7). After computing all the test statistics of the alternative functional models, i.e., for each of the visible satellites, the likelihood of the most likely model error can be determined by comparing the t i with the critical value N α/2 (0, 1). For each DIA recursion, among all the satellites with |t i | ≥ N α/2 (0, 1), the one with the maximum absolute value |t i | is then considered to be the most likely satellite affected by the model errors. After the adaptation step, the DIA process is repeated to test whether additional satellites are affected by model errors.

Adaptation
In case of an outlier, adaption implies disregarding the affected observation. However, in case the model error is affecting both the phase and code observation, it is better to adapt for the error in the functional model as where b T j = [0, ..., 1, ..., 0, ..., 1, ..., 0] T ; all the elements in this vector are 0 except for the elements of phase and code observation of jth satellite is 1, which implies a single additional parameter τ j is added to account for the unexpected model error of jth satellite. In terms of redundancy, this is better than disregarding both observations. However, as explained above, several outliers may occur or several satellites may be affected by a weather event at the same time, and in this case, more than one model error might be identified in a recursive DIA process. Then, the vector b j extends to a matrix B and the unknown parameter τ j extends to an unknown parameter vector τ. The adapted model at k epoch then reads as 0 · · · 1 · · · 0 · · · 0 · · · 0 · · · 1 · · · 0 · · · 0 · · · 0 0 · · · 0 · · · 1 · · · 0 · · · 0 · · · 0 · · · 1 · · · 0 · · · 0 . . . 0 · · · 0 · · · 0 · · · 1 · · · 0 · · · 0 · · · 0 · · · 1 · · · 0 where B is an m × q matrix, for which m is the dimension of the observation vector and q is the number of additional parameters. In each column of B, all the elements are 0. except for two 1s, which correspond to the phase and code observation of one satellite. as identified in the recursive DIA procedure.
The main procedure of this detection, identification and adaptation can be seen in Figure 1

Case Studies and Results
Two case studies will be presented in which we know there was severe weather during the observation period, so as to evaluate the capability of the proposed approach to account for associated model errors due to the asymmetrical behavior of the tropospheric delays during such events. The data of the first event is from the Australian Continuous Operational Reference Station 00NA in Darwin on 14 November 2017. In the sequence, it is referred to as Event 1. The second data set is from 24 March 2018 of an IGS permanent station DARW which is also located in the same region; in the sequence it is referred to as Event 2. In both cases there was heavy rainfall with thunderstorms on the specific days. Temperature and humidity for both days were obtained from https://www.wunderground.com/ and are shown in Figures 2 and 3, respectively. The GPS data and IGS products are used to ensure highly precise orbit and clock corrections. Since this study focuses on the model error detection, the final orbit and clock products are applied in the data processing to eliminate associated errors as much as possible. Configuration of the data processing for the real-time PPP can be seen in Table 1, in which the significance level defines the critical region where the value for test statistic lies in the null hypothesis is rejected.   As for Event 1 of Figure 2, the sun rises at UTC 21:30 (6:00 local time) and then temperature increases while humidity decreases; rainfall appears from 5:00 to 8:00 (UTC) with the temperature dropping 10 degrees within 1 h. A similar effect can also be seen in the change of humidity of that day. For Event 2 of Figure 3, the weather event appears from 2:00 to 5:00 (UTC). The area between the red lines shows the period of the weather event during which a significant influence on the positioning is present. The shadow highlights the period of the dramatic change of temperature and humidity. The hydrostatic delay depends only on the total density of the air, and the change of temperature and humidity would somehow affect the density; thus, with a high probability, the rapid shift in temperature and humidity will impact the hydrostatic delay. However, the inaccuracy of the zenith hydrostatic delay would not be a problem for the proposed model because the residuals of the hydrostatic delay will be lumped into the wet delay. In this case, DIA is to identify the model errors caused by the combined wet delay and residuals of the hydrostatic delay.   Figure 4 shows the statistics of the LOM test exceeding the threshold are mostly concentrated at the beginning of the event when the front is passing through, which causes significant spatial and temporal gradients in the integrated water vapor above the receiver. The number of subsequently identified model errors is shown as well, and here, at most two model errors are identified at one epoch, which means only one or two satellites are affected by the event at the same time. Besides, it can be seen that the statistics of the LOM test are below the threshold after the DIA procedure, indicating that there is no indication for remaining undetected model errors.
Similar behavior of the LOM test and identified model errors can also be seen in Figure 5 for Event 2; the rejected LOM test and identified model errors are mostly concentrated at the beginning of the weather event. However, outside the period of this event, there is one satellite detected to be biased at 12 : 00 UTC. Although it is difficult to prove that these model errors are caused by the tropospheric delay, the results of the up component around 12 : 00 UTC are also significantly improved, which means the proposed approach is suitable to adapt for model errors which have the same effects on both phase and code measurement.     With the proposed adapted model of Equation (16) during the weather event, the additional parameters which represent the residuals of the slant wet delay are considered to account for the model errors due to the tropospheric delay. Tables 2 and 3 are the mean and root mean squared error (RMS) of the phase and code residuals of satellites being identified with the model error during the weather events of Event 1 and Event 2, respectively. As can be seen in these two tables, residuals of the phase observations of the affected satellites are mostly reduced because the adjusted functional model is more reliable with the additional parameters accounting for the model errors. As for the residuals, the improvement of the phase observation is more significant than that of the code observation, because the value of the additional parameter mainly depends on the phase observation due to its much higher weight compared to the code observation.  On the contrary, the vertical positioning errors with the proposed method are reduced because the residual slant wet delays have been compensated by the additional parameters. Although there is still a systematic error in the east direction, the performance of the proposed method in horizontal displacement is better than that with the standard approach during the weather event. As can be seen in the skyplot of Figure 7, most of the affected satellites are located in the east part of the skyplot, which leads to a partially biased horizontal component after adjusting.   Table 4 summarizes the improvement of the results obtained from the proposed method compared to the standard approach. Significant improvement can be seen in the up component, since it is known that the tropospheric delay is one of the main error sources in the vertical direction due to high correlation. The horizontal precision of the proposed method is also improved by about 30%. For Event 2, most of the satellites affected by the weather event are located in the east part of the skyplot (Figure 9).
The distribution of the influenced satellites is shown in the skyplot of Figure 7. Similarly, the positioning errors of the up component of Event 2 in Figure 8 are also reduced, since the effects of the weather event have been removed. This approach also works well for the aforementioned model errors identified outside the period of the weather event at around 12 : 00 UTC, indicating that the model errors can be compensated if they have the same influence on the phase and code observables. Table 5 shows a significant improvement of the up component, which is the same as Event 1. Meanwhile, a system error still exists in the east component, though the precision of the horizontal component is also improved. From the skyplot of Figure 9, one can see that most of the influenced satellites are concentrated in the west part of the site, which partly causes the east-west bias of the horizontal component.   Figure  7. When PRN2 and PRN12 approached this area, they also identified with model errors, which means signals are affected by an extra tropospheric delay in this direction compared to any other azimuth angles. Then, the front moved from 150 degrees to 60 degrees, and thus satellite PRN5 and PRN20 were affected, after which it dissipated. It is worth noting that wrong identification might be present; e.g., PRN6 is affected by the weather event for a long period of time and among which PRN2 at almost the same azimuth angle as PRN6 is identified with the model error for several epochs.   Figure 9. Occurrences of the model errors (indicated in red) during Event 2 as function of (a) time and azimuth and (b) elevation and azimuth (skyplot). The blue lines represent the trajectories of the satellites, and the red points indicate that model errors are identified at those epochs. Figure 10 shows values of the estimated additional parameters which are due to the unmodeled slant wet delays caused by the weather events for this specific case. For Event 1 on the left side, the model error of PRN19 reaches up to more than 40 cm, as this satellite goes down to a low elevation angle. It seems reasonable, since the wet delay may lead to a delay of up to several meters at a low elevation angle. The mean and RMS phase residuals of PRN19 shown in Table 2 reduce, respectively, to −0.01 mm and 0.04 mm when compensated by the estimated additional parameters. The negative values of the additional parameters are because of the mismodeled hydrostatic delay. The mean and RMS phase residual of PRN6 also drop to 0.31 mm and 7.3 mm, respectively.
As for Event 2, values of PRN15 change rapidly from +20 cm to -10 cm, and they are stable for a long time span; even these additional parameters are considered to be epoch independent. This implies a further implementation of the global overall model test which takes into account the test statistics over a period of time rather than a single epoch. As can be seen in Table 3, the mean and RMS phase residuals of PRN15 reduce to -2.24 mm and 14.73 mm, respectively, which still show significant improvements compared to the standard PPP without the DIA procedure.

Conclusion
In this contribution, a DIA procedure was implemented to identify the model errors which have the same impact on both the phase and code observables; one of the applications is to account for model errors caused by tropospheric delays. An improved functional model was proposed with the additional parameters to account for the model errors. Although precise orbit and clock products were applied in the data testing to avoid any other model errors, this troposphere identification model can be easily implemented in real-time PPP, and the DIA procedure can be processed in real-time. This procedure was evaluated by two case studies of weather events, during which the tropospheric delay might be azimuthal asymmetric around the receiver, and thus the model errors due to the inhomogeneous troposphere can be detected by this procedure. The phase residuals of the satellites identified with model errors are compared to the standard approach during the weather events, since the unmodeled wet delay can at least be partly absorbed in the additional parameters. The positioning results are also improved during the events, and the improvement is most significant for the up component (72.46% and 64.41% improvement of RMS for two weather events) since the tropospheric delay is one of the main error sources in the vertical direction. The positioning performance of the horizontal component obtained from the proposed method is also improved (more than 30% improvement of RMS) compared to the standard PPP. The values of the additional parameters indicate the model errors due to the tropospheric delay can reach 40 cm when the satellite is at a low elevation angle. At most two model errors are identified at one epoch in the two case studies, which indicates that not too many satellites are affected by the asymmetrical troposphere, even during the weather events. Despite the complexity of extreme weather, the identified model errors are concentrated at the beginning of both heavy rainfall processes when the front causes significant spatial and temporal gradients of the integrated water vapor above the receiver. Besides, the satellites affected by the events are concentrated within a certain range of azimuth angle, which is related to the path of the front. This proposed procedure can also be used in monitoring severe weather. If the outliers detected by this method increase dramatically, it may indicate the front line of weather event is passing through. More testing shows that the proposed procedure may not always bring a very large improvement, but as least it does not deteriorate the positioning solutions; meanwhile, it does prevent severe impacts in some cases.