On the Problem of Double-Filtering in PPP-RTK

To obtain single-receiver Global Navigation Satellite System (GNSS) parameter solutions, the PPP-RTK user-filter combines measurements with time-correlated corrections that are separately computed by the filter of an external provider. The consequence of exercising such double-filtering is that the Kalman filter’s standard assumption of having uncorrelated measurements in time becomes violated. This leads the user-filter to lose its ‘minimum variance’ property, thereby delivering imprecise parameter solutions. The solutions’ precision-loss becomes more pronounced when one experiences an increase in the correction latency, i.e., the delay in time after the corrections are estimated and the time they are applied to the user measurements. In this contribution, we propose a new multi-epoch formulation for the PPP-RTK user-filter upon which both the uncertainty and the temporal correlation of the corrections are incorporated. By a proper augmentation of the user-filter state-vector, the corrections are jointly measurement-updated with the user parameter solutions. Supported by numerical results, the proposed formulation is shown to outperform its commonly used counterpart in the minimum-variance sense.


Introduction
Integer carrier-phase ambiguity resolution-enabled precise point positioning, PPP-RTK, has the potential to enormously benefit from the state-space packaging of the positioning corrections to reduce their transmission rate, i.e., the frequency with which the corrections are to be provided to single-receiver global navigation satellite system (GNSS) users, see, e.g., [1][2][3][4][5]. However, a reduction in the transmission rate comes at the cost of delivering time-delayed corrections. Consequently, the user is required to time-predict the corrections so as to bridge the gap between the corrections' generation time and the user positioning time.
Therefore, next to the intrinsic uncertainty brought by the randomness of GNSS measurements, PPP-RTK corrections also inherit extra uncertainty that is associated with their time-prediction [6]. Should the characteristics of such correctional uncertainty, e.g., in terms of corrections' (co)variance matrices, be made available, the PPP-RTK user would then be in a position to employ rigorous estimation methods so as to achieve the mostprecise (minimum-variance) parameter solutions [7]. In practice, however, a proper quality description of the corrections is often not provided to the user [8]. This is because the statespace representation of the corrections is aimed at minimizing the amount of information required to be made available. In the absence of such information, the approach commonly taken is to assume that the corrections are precise enough so that they can be treated as non-random [9]. The consequence of this practice is that the user estimation method becomes suboptimal in the sense that it loses its minimum-variance property and fails to provide the correct quality description of parameter solutions [10].
In our earlier contributions [7][8][9], we showed how the PPP-RTK user can limit the sub-optimality level of their parameter solutions when dealing with a 'single epoch' of data. Such single-epoch solutions do benefit from either closed-form expressions of the corrections' variance matrix or their approximate versions. The topic of the present contribution concerns the scenario where the user aims to estimate their parameters using multiple epochs of data in a near real-time manner. The well-known candidate for the real-time computation of such multi-epoch solutions is the Kalman filter [11,12]. In this contribution, we show why existing formulations of the user's Kalman filter are misspecified, and propose an alternative formulation that can largely avoid the precision loss of the user filtered solutions.
The structure of the paper is organized as follows. The problem of double-filtering in PPP-RTK is briefly reviewed in Section 2. The dependence of the user-filter on an external provider-filter is characterized, identifying the factors that lead to a misspecified stochastic model of the user-filter. In Section 3, we discuss potential choices with which the user can weight their time-correlated, 'corrected' observation vectors. Simulated examples are given to provide numerical insight into the consequences of such choices. Section 4 is devoted to the new formulation of the user-filter. The recursive structure of the filter is presented and the approximation on which the filter's optimality is based is highlighted. To numerically demonstrate the superiority of the filter over existing formulations, a single-station PPP-RTK setup [7,13] is employed in Section 5. It is thereby shown how the filter responds to rather high correction latency, i.e., the delay in time after the corrections are generated and the time they are applied to the user data. Finally, summary and concluding remarks are provided in Section 6.
The following notation will be used in this paper. The underscore symbol indicates the 'randomness' of quantities. Thus, x is random, while x is not. The hat· and check· symbols indicate the solutions of unknown parameters. Thus,x (orx) is a solution of x. The subscript t|τ ofx t|τ indicates thatx t|τ is a solution of x t which is obtained based on all the observations collected up to and including the time-instant (epoch) τ. The covariance operator is denoted by Cov(·, ·), while the capital Q is reserved for (co)variance matrices.

Optimal Provider-Filter vs. Misspecified User-Filter
We commence with the (linearized) observation equations of a single-receiver PPP-RTK user at epoch i where the user observation vector u i , together with the zero-mean random noise n i , is linked to the user's unknown parameter vector b i and the unknown correction vector c i through the full-rank design matrices B i and C i . The observation vector u i may contain the GNSS carrier-phase and pseudorange (code) measurements, with b i containing the position coordinates, carrier-phase ambiguities, receiver clock parameters, and instrumental biases. On the other hand, the correction vector c i may contain estimable forms of satellite orbit and clock parameters, atmospheric parameters, and phase/code biases [14][15][16].
With the sole use of their measurements, the single-receiver user cannot unbiasedly determine both the unknown vectors b i and c i . In other words, the augmented design matrix [B i , C i ] is rank-defect, meaning that part of b i (or of c i ) has to be held fixed as S-basis so that (1) becomes solvable for biased versions of b i and c i [15]. This of course does not suit user positioning. To obtain b i unbiasedly, the PPP-RTK user needs to receive an unbiased solution of the correction vector c i from an external provider, e.g., a network of permanent GNSS stations [1].

Non-Random Corrections as a Basic Principle
For the sake of argument, let us first assume that the correction vector c i is determined by the provider with no uncertainty and is made available to the user. Accordingly, the user would correct their observation Equation (1) as follows Since matrix B i is of full-column rank, 'single-epoch' solutions of the user parameters b i can then be unbiasedly computed. Instead of sticking to one single-epoch of data, the user may improve the precision of their parameter solutions by incorporating the temporal behaviour of some involved parameters. For instance, the user phase ambiguities remain constant in time unless a cycle-slip occurs. The rather stable instrumental biases can be linked in time by a random-walk process [17], while a polynomial dynamic model may be employed to capture the motion of the user position [7]. Such constraints between the user parameter vectors b h (h = i, i + 1 . . .) can be expressed by the following dynamic models [18] o where the randomness of the zero-sampled pseudo-observation o b h is characterized by the process noises w b h . The transition matrix Φ b links the user parameters between two successive epochs.
The constraints (3) serve as extra observations to increase the user model's redundancy, which is, the number of user observations minus the number of the estimable parameters involved. An increase in the redundancy strengthens the user model, thus improving the user parameter solutions [19]. In the context of PPP-RTK, however, the user parameter solutions are to be computed in a near real-time manner, demanding recursive estimation methods. The Kalman-filter [11,12] is known to be an optimal estimation method to handle such recursive computation in a minimum-variance sense. This optimality property relies on a key assumption, however, namely that the corrected observation vectors u h − C h c h (h = i, i + 1, . . .) must be uncorrelated in time, and uncorrelated with the process noises . Likewise, the process noises w b h must be uncorrelated in time. Provided that the user collects their measurements independently over time, the observations u h (and therefore their corrected versions u h − C h c h ) fulfill such an assumption. This serves as a basic principle of PPP-RTK that is commonly exercised in practice [8]. As the will be shown below, however, such a key assumption is violated if the non-random correction vector c h in u h − C h c h is replaced by a 'random' correction solution that has been computed by the external provider filter. As a consequence, the underlying model of the PPP-RTK user-filter becomes misspecified, losing its minimum-variance optimality property [10].

Optimal Provider-Filter
The role of a PPP-RTK provider is to estimate the unknown correction vectors c h (h = i, i + 1, . . .) and make the corresponding solutions available to the user. As with the user, the provider formulates their own Kalman filter with the observation equations (at epoch t) and the dynamic models where the provider observation vector y t , together with the zero-mean measurement noise e t , are linked to the unknown corrections vector c t through the full-rank design matrix A t . The transition matrix Φ c links the correction vectors between two successive epochs with the zero-sampled pseudo-observation o c t and time-uncorrelated process noises w c t (t = 2, 3 . . .). If the observation vectors y t (t = 1, 2, . . .) are also time-uncorrelated, and uncorrelated with the process noises w c t , the provider would then be in a position to run their optimal minimum-variance filter.
The three-step structure of the provider-filter is presented in the left-panel of Figure 1. The structure follows from an application of the least-squares principle to (4) and (5), see, e.g., [18]. At the initialization step (first epoch t = 1), the provider initializes their filter by the least-squares solutionĉ 1|1 where Q y t denotes the variance matrix of the observation vector y t (t = 1, 2, . . .). The application of the (co)variance propagation law also gives the solution's error-variance matrix as At the time-update (TU) step, the correction vector of the upcoming epoch can be time-predicted asĉ t|t−1 = Φ cĉ t−1|t−1 (t = 2, 3, . . .), with the error-variance matrix where Q w c t is the variance matrix of w c t . Finally, at the measurement-update (MU) step, the filter makes use of the upcoming observation vector y t to recursively update the correction solution asĉ t|t =ĉ t|t−1 + K c t (y t −A tĉt|t−1 ) with the error-variance matrix Qĉ t|t = (I − K c t A t )Qĉ t|t−1 , in which the Kalman gain matrix is given by where Q y t denotes the variance matrix of the observation vector y t (t = 1, 2, . . .). The application of the (co)variance propagation law also gives the solution's error-variance matrix as At the time-update (TU) step, the correction vector of the upcoming epoch can be time-predicted asĉ t|t−1 = Φ cĉ t−1|t−1 (t = 2, 3, . . .), with the error-variance matrix where Q w c t is the variance matrix of w c t . Finally, at the measurement-update (MU) step, the filter makes use of the upcoming observation vector y t to recursively update the correction solution asĉ t|t =ĉ t|t−1 + K c t (y t −A tĉt|t−1 ) with the error-variance matrix Qĉ t|t = (I − K c t A t )Qĉ t|t−1 , in which the Kalman gain matrix is given by User-filter

Misspecified User-Filter
In contrast to the provider-filter, the user-filter cannot stand on its own and requires the output of the provider-filter, that is, the correction solutionsĉ h|h (h = i, i + 1, . . .). Due to the intrinsic randomness of the provider observations,ĉ h|h is accompanied with an amount of uncertainty characterized by the error-variance matrix Qĉ h|h . With reference to (2), this implies that the random filtered solutionsĉ h|h replace their non-random versions c h to form the user corrected observation vectors u h − C hĉh|h . As the observation vectors u h are time-uncorrelated, it is the correction solutionsĉ h|h that dictate whether the Kalman filter's key assumption holds. This is followed by applying the covariance propagation law to corrected observation vectors of any two distinct epochs j and i, that is in which the covariance matrix Qĉ j|j ,ĉ i|i between the correction solutionsĉ j|j andĉ i|i is shown to read (Appendix A).
The nonzero covariance matrix above indicates that the user-corrected observation vectors u h − C hĉh|h (h = i, i + 1, . . .) are indeed time-correlated, making the user-filter misspecified.
Although the PPP-RTK user-filter is not minimum-variance, and thus suboptimal, previous studies have demonstrated that the filter can still deliver successful ambiguity-

Misspecified User-Filter
In contrast to the provider-filter, the user-filter cannot stand on its own and requires the output of the provider-filter, that is, the correction solutionsĉ h|h (h = i, i + 1, . . .). Due to the intrinsic randomness of the provider observations,ĉ h|h is accompanied with an amount of uncertainty characterized by the error-variance matrix Qĉ h|h . With reference to (2), this implies that the random filtered solutionsĉ h|h replace their non-random versions c h to form the user corrected observation vectors u h − C hĉh|h . As the observation vectors u h are time-uncorrelated, it is the correction solutionsĉ h|h that dictate whether the Kalman filter's key assumption holds. This is followed by applying the covariance propagation law to corrected observation vectors of any two distinct epochs j and i, that is in which the covariance matrix Qĉ j|j ,ĉ i|i between the correction solutionsĉ j|j andĉ i|i is shown to read (Appendix A).
The nonzero covariance matrix above indicates that the user-corrected observation vectors u h − C hĉh|h (h = i, i + 1, . . .) are indeed time-correlated, making the user-filter misspecified.
Although the PPP-RTK user-filter is not minimum-variance, and thus suboptimal, previous studies have demonstrated that the filter can still deliver successful ambiguityresolved positioning solutions when the duration of the provider-filter initialization, i.e., the time-difference between the epoch i and the initial epoch t = 1, becomes sufficiently large (e.g., ∼1 h), as can be seen in, e.g., [6,7,9]. In other words, the filtered solutionĉ i|i can become sufficiently precise so as to neglect its uncertainty relative to that of the user observations, i.e., Qĉ i|i ≈ 0. By making such an approximation, the covariance matrices in (6) and (7) become zero, meaning that the user-filter is expected to deliver minimumvariance solutions. While the duration of the provider-filter initialization can be ensured to be sufficiently long to make that approximation, the corrections cannot be instantaneously transferred to the user due to the limited data-transmission bandwidth [1]. This is all the more so as PPP-RTK is meant to take advantage of the efficient 'state-space' packaging of the corrections to reduce their transmission rate, i.e., the frequency with which the corrections are to be provided to the user. Consequently, PPP-RTK corrections ranging from orbits and clocks to phase biases are stored in, e.g., an Internet server, each having its own sampling period τ. In contrast to RTK for which the correction latency is less than 4 s [20] (i.e., τ ≤ 4), PPP and PPP-RTK state-space corrections are to be provided with higher time-delays. For instance, current PPP real-time service of IGS (https://igs.org/rts/) (accessed on 18 December 2022) disseminates state-space corrections with a typical latency of 5-10 s, see, e.g., [21,22]. In the following, it is shown how such state-space packaging brings additional correctional uncertainty.
The red-box in Figure 1 indicates the 'correction packs'ĉ kτ|kτ (k = 1, 2, . . .) that are generated and stored by the provider every τ s. As a result, the user gains access to the correction packs with a time-delay or 'latency', i.e., kτ ≤ i. The correction latency i − kτ ranges from 0 to τ − 1 s. As shown in the right-panel of Figure 1, the user picks the 'most recent' correction pack and time-predicts the correction asĉ i|kτ = Φ c (i−kτ)ĉ kτ|kτ to feed their model at epoch i. Thus, even if the approximation Qĉ kτ|kτ ≈ 0 would be plausible, the uncertainty associated with the time-predicted correctionĉ i|kτ may not be negligible. This is indeed the case when the time-difference i − kτ is considered to be large. To see this, consider the error-variance matrix ofĉ i|kτ as (Appendix A) While the first term in (8) may be considered negligible for a long duration of the provider-filter initialization, the second term increases as the latency or the time-difference i − kτ increases, leading to a misspecified user-filter. Note, for the sake of presentation, that we did not distinguish between each individual correction type (e.g., satellite orbits versus clocks) in Figure 1. We instead only show one common sampling period τ for all correction types. In practice however, each individual correction can, of course, have its own sampling period τ.
The three-step structure of the misspecified user-filter follows, in a way analogous to that of the provider-filter, be it that the role of the correction parameters c t is replaced by the user parameters b i , the design matrix A t by B i , and the observation vectors y t by the user corrected observation vectors u i − C iĉi|kτ (Figure 1, right-panel). In contrast to the provider who weights the observation vectors y t using the inverse of their variance matrix Q y t , the user may not have access to the variance matrix of the correctionĉ i|kτ to properly weight the corrected observation vectors u i − C iĉi|kτ . Instead, let us assume that the user takes a given inverse-weight matrixQ u i to replace its provider-counterpart Q y t . As shown in [10], the misspecified user-filter would then report incorrect error-variance matrices as where the Kalman gain matrix is evaluated as Here and in the following, the·-symbol over the Capital Q is used to distinguish incorrect variance matrices. Thus, not only does the user-filter deliver suboptimal solutions, but it also fails to provide the correct quality description of the user parameter solution b i|i . In the following, we discuss a choice of the inverse-weight matrixQ u i which is often adopted in practice.

On the Choice of the Inverse-Weight MatrixQ u i
In the previous section, it was shown that the main rinput of the user-filter, i.e., the corrected observation vectors u h − C hĉh|kτ (h = i, i + 1, . . .), are time-correlated cf. (6), prohibiting the recursive computation of the most precise parameter solutions. By adopting the inverse-weight matrixQ u h for u h − C hĉh|kτ , however, the user can still recursively compute suboptimal parameter solutions using their misspecified filter (cf. Figure 1). Relying on the assumption that the filtered correctionsĉ kτ|kτ are sufficiently precise to make the approximation Qĉ kτ|kτ ≈ 0, one may chooseQ u h in a way to limit the sub-optimality level of the parameter solutions. For instance, if the uncertainty involved in the timeprediction of the time-delayed correctionsĉ h|kτ can be ignored (i.e., the second term in (8) is neglected), the inverse-weight matrixQ u h can then be set to the variance matrix of the user observations u h , that is For this case, the user only takes the variance matrix of their own data, i.e., Q u h , for the weighting of the corrected data u h − C hĉh|kτ . In other words, the external correctionŝ c h|kτ are considered sufficiently precise to be treated as non-random, the scenario that is commonly exercised in practice [7,8]. This is because the corrections' error variance matrix Qĉ kτ|kτ is often not provided to the user. After all, the purpose of using state-space PPP-RTK corrections is to minimize the amount of information to be transmitted to the user [1]. The following example shows the consequence of this choice, i.e., Case 1 (10).

Example 1.
To give primary numerical insight into the consequence of choosing (10), consider a single-receiver user with a known location who is tracking the L1/L2 dual-frequency code data of a pair of satellites to determine the corresponding single-differenced (SD) slant ionospheric delay over 100 epochs with 1 Hz measurement sampling-rate. The user is given a satellite clock offset-and rate-corrections every τ = 10 s. Thus, the correction latency ranges from 0 to 9 s. In this simulation example, the filtered correctionsĉ kτ|kτ are assumed and simulated to be non-random. Thus, Qĉ kτ|kτ = 0. However, the user still has to time-predict the satellite clock correctionsĉ i|kτ at every epoch i to compute their filtered ionospheric solutionb i|i . The timebehaviour of the undifferenced ionospheric delays is modelled by a random-walk process with a standard-deviation of 1 mm/ √ s, whereas the undifferenced satellite clocks follow a constantvelocity dynamic model with a standard-deviation of 1 cm/ √ s 3 , as can be seen in, e.g., [6].
The standard-deviation of the undifferenced code data is set to 20 cm. To measure how the misspecified user filter under Case 1 performs, 1000 normally distributed samples of both the correction and user data are simulated over the fixed 100 epochs. The corresponding samples of the ionospheric estimation-error b h −b h|h (h = 1, . . . , 100), i.e., the difference between the true simulated ionospheric parameter and its filtered solution, are shown in the left-panel of Figure 2 (grey lines). The black solid lines indicate the associated 99.9% confidence-intervals, whereas the dashed lines represent the corresponding incorrect confidence-intervals which are reported by the user filter. As the number of user epochs increases, the dispersion of the estimation-error becomes smaller, showing that the precision of the user filtered solution improves. However, more than 100 epochs of data are required to have the absolute value of the 99.9% confidence-intervals smaller than 1 dm. On the contrary, the filter reports rather optimistic confidence-intervals. The considerable gap between the correct and incorrect confidence-intervals of the user parameter solutions is due to the choice made in (10), i.e., both the variance and time-correlation of the time-predicted corrections are ignored. As the random correctionsĉ i|kτ make the observation vectors u i − C iĉi|kτ time-correlated, one may employ the technique of first-order Markov state-augmentation to handle time-correlated measurements of the Kalman filter, as can be seen in, e.g., [23] (p. 180) or [24]. In that technique, the time-correlation of the corrections is 'approximated' by an exponential auto-correlation function [18], that is (compare with (7))Qĉ where the matrixQ c is to capture the uncertainty of the correction at every epoch. Next to the inverseweight matrixQ u h in (10), the state-augmentation technique incorporates the corrections' variance and time-correlation via the matricesQĉ j|j ,ĉ i|i to weight both the uncorrected user observations u h and the time-predicted correctionsĉ h|kτ . The coefficient α governs the magnitude of the correlation between the observation vectors. The larger the coefficient α, the larger the time-correlation is assumed between the observation vectors. To run their filter in recursive form, the user would need to make the following modifications [18] for epochs h > i, where the notation 'A → B' means 'replace A by B'. Thus, the user state-vector b h is augmented with the parameter vector a h whose process noises' time-correlation exponentially decays over time. For the user initial epoch h = i, the filter is initialized by the augmented state- Thus, the initial user solutionb i|i is identical to that of Case 1 (10). For the upcoming epochs h > i, however, the filter aims to capture the randomness of the corrections by updating the solutionsâ h|h over time.
To see the extent to which the state-augmentation technique can alleviate the effect of the time-correlated corrections, we evaluate the ionospheric estimation errors using the coefficients α = 50 andQ c = 0.02 m 2 . These coefficients are empirically chosen so as to approximate the time-correlation of the corrections. The results, together with their correct and incorrect 99.9% confidence intervals, are depicted in the right-panel of Figure 2. As shown, the precision of the user filtered solution slightly improves, that is, the absolute value of its 99.9% confidence-intervals reaches 1 dm after 75 epochs. The gap between the correct and incorrect confidence-intervals also becomes smaller.
As shown in Example 1, treating the time-delayed PPP-RTK corrections as non-random leads to an incorrect and optimistic quality description of the user parameter solutions. Alternatively, the state-augmentation technique may be used to 'approximate' both the variance and time-correlation of the corrections via the exponential auto-correlation function (11). Such a technique does, however, not take advantage of the information contained in the corrections' dynamic model (5), thereby discarding the structure of the error-variance matrix (8). We now consider another case that incorporates (8) into the user-filter. As stated before, the approximation Qĉ kτ|kτ ≈ 0 may hold by letting the duration of the provider-filter initialization be sufficiently long. The substitution of Qĉ kτ|kτ = 0 into (8) gives an inverseweight matrixQĉ i|kτ for weighting the time-predicted correctionsĉ i|kτ . Accordingly, the user inverse-weight matrixQ u h in (10) is modified as follows Thus, the inverse-weight matrixQĉ i|kτ only considers the second term in (8), i.e., the uncertainty due to the time-prediction of the corrections. Although the randomness of PPP-RTK corrections is taken into account under Case 2 (13), their nonzero time-correlation (7) is still dismissed. As with Case 1 (10), such a time-correlation dismissal is required to run the user-filter in its recursive form. In the next section, we show how Case 2 can be generalized so as to account for the corrections' time-correlation (7), yet allowing the recursive estimation of user parameters.

User-Filter with Correctional Update
To date, the focus has been restricted to the choice of the user inverse-weight matrix Q u h that can be made in Cases 1 and 2. Depending on how the user weights their corrected observation vectors u h −ĉ h|kτ (h = i, i + 1, . . .) throughQ u h , their misspecified filter can be recursively run in accordance with the right-panel of Figure 1. In both Cases 1 and 2, however, the time-predicted correctionsĉ h|kτ do not benefit from the information contained in the user observations u h . In other words, the corrections are merely derived from the correction packsĉ kτ|kτ , that isĉ h|kτ = Φ c (h−kτ)ĉ kτ|kτ over the epochs h = i, . . . , (k + 1)τ − 1. As with the first-order Markov state-augmentation (12), the idea is to augment the user state-vector, with a difference, so that the correction parameter vector c h now replaces the parameter vector a h . Thus, instead of approximating the corrections' time-correlation by an exponential auto-correlation function (11), the goal is now to directly incorporate the corrections' dynamic model (5) into the user-filter. To this end, we again assume that the approximation Qĉ kτ|kτ ≈ 0 holds. The initialization of the user-filter is then executed through the following modifications where the inverse-weight matrixQĉ i|kτ follows from that of (13). In contrast to Cases 1 and 2 in which the corrections are merely time-updated fromĉ kτ|kτ , we now let the user-filter, next to the time-updating, also measurement-update the correctionsĉ i|kτ . We therefore employ the discriminating notation· (instead of·) to denote the user-augmented parameter solutions During the period in which no newer correction pack is available, the user-filter recursively performs the time-and measurement-updates via the following settings At epoch i = (k + 1)τ, when the new correction packĉ (k+1)τ|(k+1)τ is made available by the provider, the new time-predicted correctionsĉ i|(k+1)τ are to replace the existing user-filtered correctionsč i|i . This follows from the approximation Qĉ (k+1)τ|(k+1)τ ≈ 0 (k = 1, 2, . . .), namely that the new correctionsĉ i|(k+1)τ contain all the information contained in the user-filtered correctionsč i|i . The structure of our proposed user-filter is presented in Figure 3. As shown, the user-filtered correctionsč i|i have to be initialized and replaced byĉ i|kτ every time a newer correction packĉ kτ|kτ is made available. Note that such a filter is still misspecified, thus delivering suboptimal solutions. However, its sub-optimality level is only dictated by the extent to which the error-variance matrix Qĉ kτ|kτ is different from zero. This is in marked contrast to Cases 1 and 2, whose solutions' loss of precision is also driven by the correction latency i − kτ. To compare the performance of the proposed filter with those of Cases 1 and 2, let us again consider the problem in Example 1. At epoch i = (k + 1)τ, when the new correction packĉ (k+1)τ|(k+1)τ is made available by the provider, the new time-predicted correctionsĉ i|(k+1)τ are to replace the existing user-filtered correctionsč i|i . This follows from the approximation Qĉ (k+1)τ|(k+1)τ ≈ 0 (k = 1, 2, . . .), namely that the new correctionsĉ i|(k+1)τ contain all the information contained in the user-filtered correctionsč i|i . The structure of our proposed user-filter is presented in Figure 3. As shown, the user-filtered correctionsč i|i have to be initialized and replaced byĉ i|kτ every time a newer correction packĉ kτ|kτ is made available. Note that such a filter is still misspecified, thus delivering suboptimal solutions. However, its sub-optimality level is only dictated by the extent to which the error-variance matrix Qĉ kτ|kτ is different from zero. This is in marked contrast to Cases 1 and 2, whose solutions' loss of precision is also driven by the correction latency i − kτ. To compare the performance of the proposed filter with those of Cases 1 and 2, let us again consider the problem in Example 1.
The filter does not only deliver the user parameter solutionb i|i , but also recursively updates the corrections asč i|i (i := i + 1). At every MU step, if a newer correction packĉ kτ|kτ is made available by the provider (i.e., if k > i−1 τ ), then the user-filtered correctionsč i|i are initialized byĉ i|kτ . Example 2 (Continuation of Example 1). Given the same 1000 normally distributed samples of both the correction and user data in Example 1, we evaluate the ionospheric estimation-errors b h −b h|h (h = 1, . . . , 100) by the user-filter using the settings om Case 2 (13) and our proposed settings as outlined in Figure 3. The corresponding results are presented in Figure 4. In comparison to the results in Figure 2, Case 2 outperforms Case 1, showing a similar performance to that of the state-augmentation technique (12), i.e., it requires at least 70 epochs to have solutions more precise than 1 dm (with 99.9% confidence). Similar to Case 1 and the state-augmentation technique, there is also a gap between the correct and incorrect confidence intervals. The proposed filter does, however, deliver the most precise results in the sense that the absolute value of their 99.9% confidence-intervals becomes smaller than 1 dm after 50 epochs. At the same time, there is no gap between the correct and reported confidence intervals. This is because the correction packsĉ kτ|kτ are simulated in a way to be non-random, Qĉ kτ|kτ = 0. In the next section, the performance of the proposed user-filter is studied for a real-world GNSS data-set for which the equality Qĉ kτ|kτ = 0 does not hold. Instead, the filter has to rely on the approximate counterpart Qĉ kτ|kτ ≈ 0.
The filter does not only deliver the user parameter solutionb i|i , but also recursively updates the corrections asč i|i (i := i + 1). At every MU step, if a newer correction packĉ kτ|kτ is made available by the provider (i.e., if k > i−1 τ ), then the user-filtered correctionsč i|i are initialized byĉ i|kτ . Example 2 (Continuation of Example 1). Given the same 1000 normally distributed samples of both the correction and user data in Example 1, we evaluate the ionospheric estimation-errors b h −b h|h (h = 1, . . . , 100) by the user-filter using the settings om Case 2 (13) and our proposed settings as outlined in Figure 3. The corresponding results are presented in Figure 4. In comparison to the results in Figure 2, Case 2 outperforms Case 1, showing a similar performance to that of the state-augmentation technique (12), i.e., it requires at least 70 epochs to have solutions more precise than 1 dm (with 99.9% confidence). Similar to Case 1 and the state-augmentation technique, there is also a gap between the correct and incorrect confidence intervals. The proposed filter does, however, deliver the most precise results in the sense that the absolute value of their 99.9% confidence-intervals becomes smaller than 1 dm after 50 epochs. At the same time, there is no gap between the correct and reported confidence intervals. This is because the correction packsĉ kτ|kτ are simulated in a way to be non-random, Qĉ kτ|kτ = 0. In the next section, the performance of the proposed user-filter is studied for a real-world GNSS data-set for which the equality Qĉ kτ|kτ = 0 does not hold. Instead, the filter has to rely on the approximate counterpart Qĉ kτ|kτ ≈ 0.

A Single-Station PPP-RTK Example
In this section, we make use of a Galileo dual-frequency (E1/E5a) data-set to study the positioning performance of the misspecified user-filter for three variants: Case 1 (10), Case 2 (13), and Case 3 ( Figure 3). The data-set was collected with a 1 Hz sampling-rate on 21 January 2022 by two permanent GNSS stations: CUT0 and UWA0, both located in Western Australia. The precise orbital corrections are a priori applied to the data.
State-space corrections (i.e., clock, ionospheric and phase-bias corrections) are generated via a single-station PPP-RTK setup [7,13]. The setup is visualized in Figure 5. Station CUT0 serves as correction-provider, whereas the station UWA0 serves as a user that is approximately 8 km away from the provider. To emphasize the performance of our filter formulation (Case 3) in handling time-delayed corrections, we consider rather high correction latencies. Accordingly, the clock correction packs are assumed to be made available to the user every 10 s, ionospheric correction packs every 30 s, and phase-bias correction packs every 10 min. In order to make the approximation Qĉ kτ|kτ ≈ 0 plausible, the provider sends the corrections only after 1 h from the time of their filter-initialization. The user then time-predicts each correction and runs their filter over 20 min (1200 epochs). To form their observation variance matrix, the user employs the sinusoidal satellite elevationweighting strategy with the zenith-referenced standard-deviations of 20 cm and 2 mm for their undifferenced code and phase measurements, respectively.  Over the fixed 1200 epochs, we consider both the ambiguity-float and -fixed positioning performance of the user. To show the role played by the correction latency, we also compare the user positioning results with those obtained with zero latency. Let us first focus on the ambiguity-float results, i.e., before resolving the user float ambiguities to their integers. The results are presented in Figure 6. The left-panel corresponds to the case where no correction latency is experienced by the user. For the zero-latency case, the results of Cases 1 and 2 coincide with those of Case 3 (green lines). This can be understood as follows. Recall that Case 1 (red lines) treats the corrections as non-random, and Case 2 (blue lines) models the variance of the corrections, while Case 3 models both the variance and time-correlation of the corrections. As both Cases 2 and 3 rely on the approximation Qĉ kτ|kτ ≈ 0, they deliver results identical to those of Case 1 in the case of zero latency (i.e., when i = kτ). Thus, the underlying difference of these three cases only lies in the way of correctional uncertainty due to the fact that the latency is handled.
The effect of nonzero correction latency is shown in the right-panel of Figure 6. Since i ≥ kτ, the user has to time-predict the corrections, thereby requiring that it take the associated correctional uncertainty into account. As Case 1 does not consider such uncertainty, it delivers positioning results with the largest root-mean-squared (RMS) errors in all the three East-North-Up directions. For the horizontal components, Case 2 delivers similar results to those of Case 3 in an RMS sense. However, Case 3 experiences less discontinuities in its positioning results than Case 2. Such discontinuities are owed to the periodic increase in the latency that it varies from 0 to τ − 1 seconds for every correction pack. Let us now turn our attention to the user ambiguity-fixed results, i.e., after mapping the float ambiguities to their integers. At every epoch, full ambiguity resolution using the integer least-squares estimation is conducted without an extra ambiguity validation procedure [25]. The results are presented in Figure 7. As before, both the zero-latency (left-panel) and nonzero-latency (right-panel) cases are considered. While the three cases 1, 2, and 3 perform the same when the user receives corrections with no time-delay, their performances become distinct when treating time-delayed corrections. Case 1 (red lines) exhibits several large ambiguity-resolved positioning errors (≥1 dm) in all three East-North-Up directions, leading to large positioning RMS errors. In other words, the corresponding ambiguity success-rate is not sufficiently high to deliver successful ambiguity-resolved positioning solutions. On the contrary, Cases 2 and 3 outperform Case 1 in the RMS sense.
Similar to their ambiguity-float counterparts, the RMS errors of Cases 2 and 3 are the same for the horizontal components. Observing the results in Figures 6 and 7, one may be inclined to conclude that the performance of the proposed user-filter (Case 3) is the same as that of Case 2. From a computational point of view, Case 2 seems to be more attractive as it does not need the augmentation of the user state-vector with the corrections. This is in contrast to the formulation of Case 3 which needs to recursively update both the user and correction parameter solutions. One should, however, remark that those results only represent one realization of the user-filter's parameter solutions. In order to infer the overall performance of the user-filter under the formulations offered by Cases 1, 2, and 3, one needs to instead consider the distribution of the user parameter solutions. To that end, we generate 300 different realizations of the results by shifting the user-filter starting epoch i every 15 s. The time-series of the medians (i.e., 50% percentiles) of these realizations are presented in Figure 8 for both the user ambiguity-float (left) and -fixed (right) cases. The medians of the positioning errors corresponding to Cases 2 and 3 are shown to be considerably smaller than those of Case 1. The results also indicate that Case 3 outperforms Case 2 as it, on average, delivers smaller medians of the positioning errors. Note also the presence of periodic jumps of the medians for all the three cases. This behaviour is due to the periodic nature of the correction latencies that vary from zero to τ − 1 s for each correction pack. It is, however, observed that Case 3 exhibits smaller periodic jumps. This is because Case 3 accounts for the uncertainty of the time-delayed corrections by modelling their variance and time-correlation through their dynamic model (5).  Figure 8. Medians (50% percentiles) of the absolute positioning errors corresponding to 300 user-filter realizations before (left) and after (right) ambiguity-fixing in the East-North-Up directions when the correction packs are provided with a time-delay (cf. Figure 5). The results of Cases 1, 2, and 3 are indicated in red, blue, and green, respectively.

Conclusions and Outlook
In this contribution, a new multi-epoch formulation for the PPP-RTK user-filter was proposed ( Figure 3). Under the proposed formulation, the user-filter state-vector b h is augmented with the correction parameter vector c h as [b T h , c T h ] T , thereby allowing the corrections to be jointly measurement-updated with the user parameter solutions.
It was shown why the user-filter always remains misspecified, and therefore, suboptimal in the minimum-variance sense, no matter which formulation is adopted. Since the PPP-RTK provider has the freedom to disseminate their filtered correctionsĉ kτ|kτ only after a sufficiently long period from the time of their filter-initialization, the user can, however, benefit from the approximation Qĉ kτ|kτ ≈ 0 to limit the sub-optimality level of their filter. Instead of current formulations in which the time-correlated corrected observation vectors u h − C hĉh|kτ (h = i, i + 1, . . .) serve as input of the user-filter, our proposed formulation treats the correctionsĉ h|kτ as additional observation vectors and works with the augmented observation vector [u T h ,ĉ T h|kτ ] T . Relying on the approximation Qĉ kτ|kτ ≈ 0, the observation vectors [u T h ,ĉ T h|kτ ] T (k = 1, 2, . . .) would then become time-uncorrelated, leading the user-filter to deliver close-to-minimum-variance solutions. Under the proposed formulation, however, the user needs to know the dynamic model that the provider uses for the corrections. For further standardization of the State Space Representation (SSR) corrections (https://www.igs.org/formats-and-standards (accessed on 18 December 2022)), Radio Technical Commission for Maritime Services (RTCM) committees may therefore request PPP-RTK providers to disclose the dynamic models underlying their SSR corrections.
A single-station PPP-RTK setup was employed to numerically demonstrate the superiority of the proposed filter over existing formulations. To emphasize the performance of our formulation, rather high correction latencies were considered, e.g., ionospheric correction packs were provided every 30 s. It was observed that the proposed filter delivers smaller estimation errors (in an RMS sense) when handling time-delayed corrections.
In the present contribution, attention was focused on the formulation of the user-filter only. Addressing open research questions such as 'how to develop measures for assessing