Functional model modification of precise point positioning considering the time-varying code biases of a receiver

Precise Point Positioning (PPP), initially developed for the analysis of the Global Positing System (GPS) data from a large geodetic network, gradually becomes an effective tool for positioning, timing, remote sensing of atmospheric water vapor, and monitoring of Earth’s ionospheric Total Electron Content (TEC). The previous studies implicitly assumed that the receiver code biases stay constant over time in formulating the functional model of PPP. In this contribution, it is shown this assumption is not always valid and can lead to the degradation of PPP performance, especially for Slant TEC (STEC) retrieval and timing. For this reason, the PPP functional model is modified by taking into account the time-varying receiver code biases of the two frequencies. It is different from the Modified Carrier-to-Code Leveling (MCCL) method which can only obtain the variations of Receiver Differential Code Biases (RDCBs), i.e., the difference between the two frequencies’ code biases. In the Modified PPP (MPPP) model, the temporal variations of the receiver code biases become estimable and their adverse impacts on PPP parameters, such as ambiguity parameters, receiver clock offsets, and ionospheric delays, are mitigated. This is confirmed by undertaking numerical tests based on the real dual-frequency GPS data from a set of global continuously operating reference stations. The results imply that the variations of receiver code biases exhibit a correlation with the ambient temperature. With the modified functional model, an improvement by 42% to 96% is achieved in the Differences of STEC (DSTEC) compared to the original PPP model with regard to the reference values of those derived from the Geometry-Free (GF) carrier phase observations. The medium and long term (1 × 104 to 1.5 × 104 s) frequency stability of receiver clocks are also significantly improved.


Introduction
In exploring the potential of Global Positioning System (GPS) for a variety of applications, Precise Point Positioning (PPP) has been developed as a tool for processing code and phase observations from a stand-alone GPS receiver at the undifferenced level (Zumberge et al. 1997a) along with the use of precise satellite orbit and clock products provided by the International GNSS (Global Navigation Satellite System) Service (IGS) (Kouba and Héroux 2001). The PPP can deliver various types of parameters, including station positions, receiver clock offsets, Zenith Troposphere Delays (ZTDs), and slant ionosphere delays, which are of great importance for practical uses. The positioning accuracy is at the level of a few centimeters in static mode and below one decimeter in kinematic mode (Bisnath and Gao 2009). This suggests PPP can be used for crustal deformation monitoring (Wright et al. 2012;Xu et al. 2013), marine surveying (Alkan and Öcalan 2013;Geng et al. 2010), and land-vehicle navigation (Rabbou and El-Rabbany 2015;Wielgosz et al. 2005). In recent years, the efficient and flexible PPP technique has also played a crucial role in several non-positioning applications, especially in atmospheric (ionosphere and troposphere) sounding (Rovira-Garcia et al. 2015;Shi et al. 2014;Yuan et al. 2014) as well as time and frequency transfer Ge et al. 2019;Orgiazzi et al. 2005;Tu et al. 2019).
In the implementation of PPP, one needs to formulate the functional model (i.e., observation equations), relating the GPS observations to the parameters to be estimated. One assumption underlying this formulation is that the receiver code biases do not change significantly over time (Banville and Langley 2011b;Håkansson et al. 2017). A vast amount of work has cast considerable doubt on the validity of this assumption (Bruyninx et al. 1999;Coster et al. 2013;Wanninger et al. 2017) and found that the phenomenon of receiver code bias variation, which is closely related to the ambient temperature, is widespread. For example, the Geometry-Free (GF) receiver code biases, known as receiver Differential Code Biases (DCBs) (Montenbruck et al. 2014;Sardon et al. 1994), has been found to exhibit an apparent intraday variability of 4 ns to 9 ns (Ciraolo et al. 2007), far exceeding the code noise level (depending on the receiver type, but generally smaller than 1 ns). If the assumption of time-constant receiver code biases is taken, PPP solutions will be biased and hence the performances in PPP applications will be degraded. For the ionospheric Slant Total Electron Content (STEC) retrieval, the Modified Carrier-to-Code Leveling (MCCL) method , as well as the integer-levelling procedure (Banville and Langley 2011a;Banville et al. 2012) can both effectively eliminate the effect of receiver code bias variations. As far as the GPSbased timing application is concerned, there are still a limited amount of research focusing on how to cope with the time-varying receiver code biases.
In this study, a modified version of PPP is proposed, whose functional model contains time-varying receiver code biases, and is formulated based on raw observations instead of the Ionosphere-Free (IF) combinations. The rank deficiencies are removed, which are caused by the functional model and lead to the non-estimability of all the parameters, by fixing a minimum set of parameters to their prior values, resulting in a full-rank system of observation equations (Odijk et al. 2016;Teunissen 1985). In this system of observation equations, the temporal variations of the receiver code biases on two frequencies become estimable. It will be shown that the ambiguity parameters, ionospheric delay, and receiver clock offset directly absorb the combination of timevarying receiver code biases, which are normally considered as constant in subsequent ionospheric modeling (Liu et al. 2017(Liu et al. , 2018 and time transfer (Ge et al. 2019;Tu et al. 2019). This latter assumption causes adverse effects in these applications will be shown. In this regard, the Modified PPP (MPPP) functional model in ionospheric STEC retrieval and timing will be the focus of the research. The precise satellite orbit and clock products from an external provider such as the IGS are applied as deterministic corrections when conducting PPP, hence the MPPP proposed in this work is to be considered as a geometry-based method. The model can simultaneously obtain the variations of receiver code biases at two different frequencies using raw observations. In contrast, the MCCL method that is actually a kind of geometry-free method  can only detect the betweenepoch fluctuations of Receiver Differential Code Biases (RDCBs) and cannot be applicable to time or frequency transfer because the receiver clock parameter is eliminated in the geometry-free combinations of code and phase observations.
The organization of this work proceeds as follows. In "Methods" section, the full-rank functional model for the original as well as the modified PPP are constructed, followed by the interpretation of the estimable parameters. "Results" secrtion provides numerical insights into the ability of the MPPP to exclude the effects that the shortterm temporal variability of the receiver code biases has on ambiguity estimation, ionospheric STEC retrieval and timing. "Conclusions" section concludes the study with a short discussion.

Methods
In this section, the rank-deficient system of GPS observation equations is taken as a starting point, and then the two PPP models are presented (i.e., original and modified), focusing mainly on developing the functional model and interpreting the estimable parameters.

GPS observation equations
Let us consider a scenario where m satellites, transmitting signals on f frequencies, are tracked by a single receiver over t epochs. Under this scenario, the system of observation equations has the following form (Leick et al. 2015), with r , s = 1, . . . , m , j = 1, . . . , f and i = 1, . . . , t being the receiver, satellite, frequency and epoch indices, respectively, and where p s r,j (i) and φ s r,j (i) denote, respectively, the code and phase observables, l s r (i) the receiversatellite range, τ r (i) the zenith troposphere delay and m s r the corresponding troposphere mapping function, dt r (i) and dt s (i) the receiver and satellite clock offsets, d r,j and d s ,j the frequency-dependent receiver and satellite code biases, ι s r (i) the (first-order) slant ionosphere delay on the first frequency ( µ j = 2 j 2 1 ) and a s r,j the (non-integer) ambiguity with j the wavelength of frequency j . Note that all quantities are in unit of range, and the time-constant parameters do not have an epoch index i.
For positioning purposes the primary parameter of interest is the three-dimensional receiver position vector. To employ the least squares adjustment one needs to linearize the above observation equations. The linearized form of Eq. (1) is as follows (Teunissen and Kleusberg 2012): where �p s r,j (i) and �φ s r,j (i) are the code and phase observables that are corrected for the approximate receiver-satellite ranges and the satellite clocks. Note that the satellite positions and clocks computed using the IGS final products are not part of the parameters, and the receiver coordinates are fixed to the values in the IGS Solution Independent Exchange (SINEX) product and do not need to be estimated. Also the derived model shall not change if one further considers the receiver positions as unknown parameters. The ionosphere-free satellite code bias, d s ,IF , appears in both code and phase observations, because of its presence in the introduced satellite clocks (Zumberge et al. 1997b).
In the following two sections, only the dual-frequency case is considered (that is, f = 2 ) to keep the presentation of the functional model as simple as possible.

The original version of PPP
Equation (2) represents a rank-deficient system, implying that some of the parameters are not unbiased estimable, but only combinations of them. This becomes clearer if we write , with d r,IF the ionosphere-free receiver code bias, d r,GF and d s ,GF the geometry-free receiver and satellite code biases. It then follows, d r,IF is not separable from dt r (i) , and the ι s r (i) , d r,GF and d s ,GF are not separable from one another. The idea is therefore to reduce the number of parameters by lumping some of them together. For instance, with the lumped parameters, the code observation equation is of full-rank and reads, with dt r (i) the biased receiver clock and ι s r (i) the biased slant ionosphere delay in (4).
Inserting dt r (i) and ι s r (i) into the phase observation equation and lumping a s r,j with d s ,IF , we eliminate the rank deficiency and then obtain with being the biased ambiguity.
Equations (5) and (6) represent the full-rank functional model for the original PPP, in which all estimable parameters are interpreted by Eqs. (4) and (7). Recall that it is a common practice to use the (recursive) leastsquares estimator to solve for the parameters. Bearing this in mind, a few remarks on the effects of time-varying receiver code biases are thus in order. Firstly, the timeconstant of a s r,j , an implicit assumption made to exploit the very precise phase observable, becomes invalid. This, in turn, can introduce larger errors in comparison with the formal errors of all parameters. Secondly, the use of dt r (i) for timing and time transfer can be susceptible to batch boundary discontinuities (Collins et al. 2010;, because time averaging of d r,IF induces clock datum changes between batches. The variability of d r,IF , which does not average out to zero, limits the use of dt r (i) for frequency transfer, inducing worse frequency stability than one would expect theoretically (Bruyninx et al. 1999). Thirdly, when using the thin-layer ionosphere model for determining vertical TEC from ι s r (i) , a temporal variation of d r,GF is partially responsible for the occurrence of the model error effects (Brunini and Azpilicueta 2009).

The modified version of PPP
Let us consider again Eq.
(2), but replace d r,j by d r,j (i) , implying that the receiver code bias is allowed to vary freely over time. Thus Note that if the d r,j (i) is ignored, the (rank-deficient) design matrix of Eq. (8) turns out to be the same as that of Eq. (2). In the following d r,j (i) is assumed as a timevarying parameter. But there will be a rank deficiency between dt r (i) , d r,j (i) and a s r,j parameters. The receiver code biases at the first epoch are chosen as datum to eliminate the rank deficiency, thus the estimated d r,j (i) are the variations of receiver code bias with respect to the first epoch, see Eq. (9). As for other parameters, the procedures described in the preceding subsection are closely followed to overcome the rank deficiency problem, thereby yielding with where a tilde marks the biased, but estimable parameters. Note the difference to the estimable parameters in Eqs. (5) and (6).
Equation (10) is a full-rank system, representing the full-rank functional model constructed for the modified PPP. Note that, the estimable receiver code biases d r,j (i) only appear in the observation equations at the second epoch and beyond, i.e., i ≥ 2 , for those at the first epoch get lumped with the parameters given in Eq. (11). The estimability of d r,j (i) implies that any temporal variations in receiver code biases shall fully enter into d r,j (i) , and thus have no negative impacts on the estimation of remaining parameters.

Results
The present section starts with a description of the experimental data and processing strategies, then proceeds with the results and analysis, and ends with some concluding remarks.

Experiment setup
The data for this analysis were collected at four stations with dual-frequency observations and a 30-s sampling interval. The detailed information is given in Table 1. Note that stations ALIC and MTDN have thermometers to gauge the air temperatures, and the receivers at stations NOT1 and MEDI were both connected to a high performance external frequency standard (H-maser).
Our data processing strategies are as follows. Firstly, the original as well as the modified PPP were conducted, solving for the parameters by means of Kalman filter. The estimated parameters included the ZTDs (a random walk process with a variance rate of 1 × 10 −7 m 2 /s), the biased receiver clocks (time-varying), the biased slant ionosphere delay (time-varying), the biased ambiguities (time-constant), and the biased receiver code biases (time-varying for the modified PPP). The elevation cutoff angle was set as 10°, and the P1-C1 satellite DCB corrections provided by the IGS were applied to the C1 observations. The IGS final orbit and clock products were employed to correct for the satellite orbital and clock errors. The other critical corrections to raw observations were also considered, including the solid Earth tide, the phase wind-up effects, and the satellite and receiver phase center offsets and variations.
Second, The variations of the code observations were estimated on each frequency epoch by epoch. The results would allow us to verify the ability of the modified PPP to directly measure the temporal variability of the receiver code biases for each observable type.
Third, the time-varying receiver code biases are not considered in Global Ionosphere Maps (GIMs) provided by the IGS, and the GF carrier phase observables are almost unaffected by measurements noise and multipath. The Difference of Slant Total Electron Content (DSTEC) obtained from the GF carrier phase observations is thus chosen as the reference to evaluate the performance of Fourthly, in order to validate the feasibility and effectiveness of the MPPP model in timing, the Allan DEViation (ADEV) was used to evaluate the frequency stability of the receiver clock solutions. The Pearson Correlation Coefficients (PCC) was used to measure the correlation between the time-varying receiver code biases and the intraday temperature (Zha et al. 2019). The corresponding PCC are all above 0.95. The ionosphere-free (gold lines) and GF (purple lines) combinations of P1 and P2 code bias variations also have a strong correlation with the intraday temperature (red lines), which is consistent with the previous research (Zha et al. 2019). Note that the variation values of receiver code biases and their GF/IF combinations in subplots (c), (d), (e) and (f ) in Fig. 1 were reversed to make them consistent with the variations of temperature. Therefore, their variation series look consistent with the temperature, but their PCC values are negative. Figure 2 shows the estimates of P1 (red lines) and P2 (green lines) receiver code biases by the MPPP proposed in this work and their ionosphere-free (gold line) and GF (blue line) combinations on an epoch-by-epoch basis at stations NOT1 and MEDI in two consecutive days.

Intraday variations of receiver code biases
It can be inferred from Fig. 2 that the principal factor responsible for this change cannot be the multipath effects, which should have a sidereal repeat period (Agnew and Larson 2007;Axelrad et al. 2005), but should be the difference of code bias variations for the two consecutive days.
According to the figures (Figs. 1, 2), the variations of P1 and P2 code biases range from a few nanoseconds to tens of nanoseconds. In this test the most significant change occurs at station MTDN on day 005 in 2018, with a peakto-peak range of about 30 ns for the variation of P1 and P2 code biases as well as their ionosphere-free combinations. This is bound to bring disastrous effects on the estimated parameters, especially for ambiguity parameters, ionospheric STEC retrieval and receiver clock offsets, refer to Eq. (11).
It can be seen from Eq. (7) that the GF and ionospherefree combinations of receiver code biases both enter ambiguity parameters when conducting the original PPP. The time-varying receiver code biases will therefore impair the ambiguity estimation performance and make them not constant because of the incorrect functional model. Let us take station MTDN associated with the largest variation of receiver code biases as an example. Figures 3 and 4 show the estimated ambiguities of L1 frequency based on the original PPP and MPPP models, respectively.
With the original PPP model, ambiguity estimates are subject to variations, not constant (see Fig. 3). Furthermore, the residuals obtained do not follow the normal distribution (see Fig. 5) owing to the effect of receiver code biases.
This will inevitably induce leveling errors in the STEC estimation. In contrast, with the MPPP model (see Fig. 4)  where the variations of receiver code biases become estimable, the ambiguity estimates converge to constant reasonably and the residuals of code observations follow the normal distribution (see Fig. 6).
At this stage it is proved that MPPP model can get correct ambiguity estimates and reasonable code observation residuals. In the following sections the improvements in STEC retrieval and timing based on the MPPP model will be investigated.

Retrieval of STEC parameter
As mentioned earlier, also see Ciraolo et al. (2007), if the short-term variability of receiver DCBs is not properly managed, it can have significant impact on the STEC results. In order to evaluate the performance of the MPPP in retrieving STEC, we compared the DSTECs obtained by the original and modified PPP models was compared with those derived by using the GF carrier phase observations to minimize any impact from multipath or code biases. The corresponding results are shown in Fig. 7, where different colors represent different satellites.
It can be clearly seen that the STECs retrieved by the MPPP (right column) are more accurate and stable than those by the original PPP model (left column). Table 2 summarizes the Root-Mean-Square (RMS) values for the original PPP-DSTEC as well as MPPP-DSTEC. Note that the most obvious improvement of 96% occurs at station MTDN, which also has the significant variations of code biases on P1 and P2 as shown in Fig. 1. The other three stations (ALIC, NOTI and MEDI) also show considerable improvements. Considering the fact that the MPPP results are free of the effects due to variability of code biases, we can thus interpret the differences as the errors to which the original PPP results are subject.

Timing
Generally, the colored signature of the code noise is an interfering factor for the medium and long-term (from some hours to some days) stability of receiver clock solutions (Martínez-Belda et al. 2011;Martinez-belda et al. 2012). The MPPP (see Eq. 10) can theoretically eliminate the effect caused by the time-varying receiver code biases. To validate its effectiveness, the data at station NOT1 on day 001 in 2018 and station MEDI on day 009 in 2018 were used. Figure 8 shows the Allan deviation of receiver clock offsets at station NOT1 and MEDI, respectively, which are both connected to the external Hydrogen atomic clock.
The blue lines and red lines represent the Allan deviation of the PPP-derived and MPPP-derived receiver clock offsets, respectively. It can be clearly seen that the frequency stability of the MPPP-derived solutions is improved significantly relative to the original   For stations NOT1 and MEDI, the medium and long term stabilities (1 × 10 4 to 1.5 × 10 4 s) of receiver clock are improved from 2.18 × 10 −14 to 4.75 × 10 −15 and 4.08 × 10 −14 to 1.04 × 10 −14 , respectively. Therefore, the MPPP model can greatly improve the medium and long term (1 × 10 4 to 1.5 × 10 4 s) frequency stability of PPP-based timing.

Conclusions
As well known, receiver code biases may exhibit dramatic variations in hours or less. Therefore, the original PPP model cannot give optimal results, because it implicitly assumes receiver code biases are constant over the time  course of interest. To account for this, the PPP functional model was modified by assuming that the receiver code biases can vary freely in time. By means of re-parameterization the rank deficiencies were overcome, resulting in a full-rank functional model in which the variations of receiver code biases on both frequencies between epochs are directly estimated, and thus have no impact on the performance. A series of experiments were carried using the observations at four stations equipped with dual-frequency receivers. By comparing the results using the original and modified PPP, the following conclusions can be drawn. First, the modified PPP is capable of determining the short-term temporal variability of the bias associated with a single receiver and a code observable. The intraday variations of receiver code biases are indeed quite significant in some cases (see MTDN in Fig. 1), varying in a magnitude of a few nanoseconds to tens of nanoseconds and showing a strong correlation with the intraday temperature. Secondly, for the estimated STEC parameters the accuracy of MPPP-based STECs was improved by 46% to 96% compared with the PPP-based STECs. This was evaluated by comparing with a set of reference DSTEC values obtained using the GF carrier phase observations. Thirdly, for the timing, after eliminating the influence of the time-varying receiver code biases with the modified PPP model, the medium and long term (1 × 10 4 to 1.5 × 10 4 s) frequency stability of the MPPPderived receiver clock was improved remarkably relative to the PPP-derived ones.