Ionospheric Delay Handling for Relative Navigation by Carrier-Phase Differential GPS

The paper investigates different solutions for ionospheric delay handling in high accuracy long baseline relative positioning by Carrier-Phase Differential GPS (CDGPS). Standard literature approaches are reviewed and the relevant limitations are discussed. Hence, a completely ionosphere-free approach is proposed, in which the differential ionospheric delays are cancelled out by combination of dual frequencyGPSmeasurements.The performance of this approach is quantified over real-world spaceborneGPS data made available by the Gravity Recovery and Climate Experiment (GRACE) mission and compared to the standard solution.


Introduction
Carrier-Phase Differential GPS (CDGPS) is a proven technology in several fields of application.CDGPS has been already employed for relative positioning of Low Earth Orbit (LEO) satellites flying in formation [1][2][3][4], of aircraft with respect to runways [5] and for cooperative self-separation of general aviation aircraft [6].The capability to achieve high accuracy by CDGPS is based on the possibility to exploit the integer nature of Double Difference (DD) carrierphase ambiguities [7].However, as the separation among the satellites increases, the correlation of ionospheric delays among the receivers decreases [8].As a result, DD GPS observables are affected by significant errors that complicate the integer resolution task.This paper investigates the effects of different strategies for ionospheric delay compensation on the accuracy in the relative positioning of GPS receivers in LEO over long baselines.Its results can be extended, at least in principle, to other formation flying applications [9], such as those involving Very Light Jets [10] and/or Unmanned Aerial Systems [11][12][13].
Different approaches exist in the literature for dealing with ionospheric delays.In high accuracy, postprocessing applications with dual frequency data [14,15], the DD ionospheric delays are estimated within a dynamic filter [16], for example, the Extended Kalman Filter (EKF), and are modelled by very simple stochastic models, typically using random walk processes in the filter's state vector.As an alternative, delays are modelled by Lear's model [1,17] which allows relating the  slant ionospheric delays to the Vertical Total Electron Content (VTEC) above the receivers.Even though modelling the ionospheric delays helps to increase their observability, and thus to aid in the ambiguity resolution task, Lear's model is known to be structurally capable of reproducing actual ionospheric delays only to a limited extent [8].
It is well known that linear combinations of GPS dual frequency measurements can be used to delete the ionospheric delays to first order.When those combinations are used as an input to the EFK, it can be expected that the magnitude of the ionospheric delays does not affect the achievable relative positioning accuracy.However, the use of an ionosphere-free approach is known to complicate significantly the integer ambiguities resolution task if compared to approaches that attempt to model the ionosphere [18,19].Hence, the choice between model-based methods (e.g., based on Lear's model) and ionosphere-free approaches is in general not trivial [20], depending on the relative distance between the receivers, on relative dynamics, and on the status of ionospheric activity.In this sense, the ionospheric activity plays a major role in determining the set of GPS measurements and combinations to process to improve the relative positioning accuracy.

International Journal of Aerospace Engineering
In this paper, a completely ionosphere-free approach is pursued, in which the ionospheric delays are cancelled out by combination of dual frequency GPS measurements.Several alternative combinations are investigated based on the ionosphere-free combination of pseudorange and carrierphase observables, but also on GRoup And PHase Ionospheric Correction (GRAPHIC) and Melbourne-Wubbena combinations and thus the best combinations are selected.Based on a relative positioning scheme previously developed by the authors [17], the performance of each approach is quantified over real-world spaceborne GPS data made available by the Gravity Recovery and Climate Experiment (GRACE) mission.
The paper is organized as follows.First the conventional approach is presented in which the ionospheric delays are estimated as part of the state vector through Lear's model.Then ionosphere-free observables are derived from the DD observation models, and two different combinations are selected for being integrated in an ionosphere-free formulation of the EKF.The developed filter is finally tested on GRACE data.

Ionospheric Delays Estimated by Lear's Model
The most common approach for handling the differential ionospheric delays in high accuracy long baseline applications is to estimate the DD ionospheric delays within the EKF modelling them using random walk processes in the state vector [14,15].This scheme usually integrates an EKF with an Integer Least Squares estimator based on the LAMDBA method.This approach can be not sufficiently accurate in real-time on-board implementations when accurate stochastic and dynamic models cannot be used.In such conditions, a reliable way to proceed is to model the ionospheric delays through the VTEC above the receivers by Lear's mapping function [1,8].With specific reference to a formation of two satellites, this approach leads to the following state and measurement vectors: where b and ḃ are the baseline and the relative velocity vectors expressed in the Earth-Centred Earth-Fixed (ECEF) reference frame and N 1 stack together the cycle ambiguities on the L1 frequency for the  DD couples: where the pivot satellite, that is, the reference GPS satellite selected for calculating DD observables, is indicated with 0 for simplicity.N WL , instead, stacks together wide-lane cycle ambiguities.In ( where  1 is L1 carrier frequency in Hz, VTEC  is the (scalar) vertical total electron content for modeling the ionospheric delay of the receiver  and is expressed in number of electrons per square meter, and    is the elevation of the GPS satellite  with respect to the receiver .
The main advantage of Lear's model is the capability to predict zero difference ionospheric delays, relevant to different tracked GPS satellites, as a function of a single unknown parameter (i.e., VTEC A ), which is a desirable property for navigation filters [1].
This reference model of ( 1)-( 3) has shown satisfactory observability features [17], since it is capable of delivering good estimates of the Integer Ambiguities (IA) even in case of only 3 DD observations.However, this VTEC-based EKF has some inherent accuracy limitations, due to its inability of rejecting deviations of the true ionosphere from Lear's model, which appear as additional error terms in the baseline estimate [8].In what follows, the possibility is discussed of overcoming the limitations of the VTEC-based EKF by deleting the ionospheric delays through measurement combinations.

Ionosphere-Free Observables
Dual frequency DD carrier-phase and pseudorange observables can be combined in different ways to generate ionosphere-free measurements.The complete DD observation model is thus presented before deriving the relevant ionosphere-free observables.

DD Observation Model.
The DD measurements have the following expressions [21]: where (i)  =  2 / 1 is ratio between the L2 and L1 frequencies, and the  and  superscripts refer to the GPS satellites radiating the navigation signal; (ii)  1 ,  2 are the wavelengths of L1 and L2 signals, respectively; (iii) is the integer ambiguity on the L1 frequency (on L2 which is analogous);

𝑗
is the  1 noise term on the L1 frequency (on L2 which is analogous).
Each of the four observables in ( 4) is assumed to be independent from the other ones.However, the  DD measurements of the same kind at a certain time epoch are mutually correlated due to the presence of the pivot satellite in all the measurements.More precisely, denoting generically by  the observation type, that is,  =  1 , . . .,  2 , and by   its standard deviation, we have ) .(5)

Ionospheric Free Combination.
The most common combination for eliminating the ionospheric delay, referred to as Ionosphere-Free (IF) [21], is concerned with combining observations of the same type on the two carrier frequencies, exploiting the frequency dependence of the first-order ionospheric delay effect.More precisely, the IF combinations are obtained as Thus, two IF observables per each DD couple  out of the four measurements in (4) are The IF observation model can be obtained combining ( 4), ( 6), and (7) which yields Since the observation types are assumed to be independent, assuming white Gaussian measurements noises, the IF combinations are affected by Gaussian white noise with variance Hence, the noise is increased compared to the original uncombined observations.

GRAPHIC Combinations. GRAPHIC combinations
exploit the asymmetry of the ionospheric effect on group and phase propagation [22].In practice, they combine pseudorange and carrier-phase measurements on each frequency, as follows: From the above equations, the GRAPHIC observation model and variance read where  1 ,  2 indicate the noise terms of  1 and  2 combinations, respectively.

Melbourne-Wubbena Combinations.
Melbourne-Wubbena (MW) combinations combine all four observable types for cancelling out the ionospheric delay [23,24].They build upon the definition of the wide and narrow lane (NL) wavelengths: where  is the speed of light in vacuum.The MW combinations are obtained as and have the following observation model and variance: where  represents the noise term of the combination. 2

Ionosphere-Free Relative Positioning
This section is concerned with establishing which of the measurement combinations presented in Section 3 is suitable for computing the relative position of the two receivers in long baseline applications.The measurement combinations available with no ionospheric effects are the   IF combinations, the   IF combinations, the  GRAPHIC combinations on L1 and L2, and the  MW combinations, for a total of 5n measurements.However, the 5 measurements are not linearly independent.In particular, each group of 5 observables per each of the  DD couples can be seen as a linear transformation of the 4 uncombined measurements of (4).More precisely, the cancellation of the ionosphere from (4) can be seen as a linear projection of the 4dimensional measurement vector [ 1 ,  2 ,  1 ,  2 ]  onto a 3dimensional hyperplane.This implies that no more than three linear combinations of the four dimensional vector exist being linearly independent, and thus suitable for use as measurement vector in an EKF.In addition, whatever set of three linearly independent vectors lying on the hyperplane can be used as a basis for describing all the vectors belonging to the hyperplane.Even though any of such bases would yield theoretically equivalent positioning solutions, significant differences exist in practical applications due to unmodeled systematic errors.In order to find such "best" basis, its performance in accurately estimating both the geometric term and the cycle ambiguity is analyzed.A natural indicator to evaluate the accuracy achieved by a specific measurement combination is the standard deviation (STD) of the noise affecting the measurements.
Hence, a specific analysis has been performed.Starting from different candidate noise levels for the uncombined measurements, that is,  1 ,  2 ,  1 , and  2 (see the first three columns in Table 1), reflecting typical performance of spaceborne GPS receivers, the expected STD for the estimation of the geometric term   is derived as listed in Table 1.The rationale of this analysis is that different combinations of GPS measurements may modify, that is either amply or reduce, the uncertainty of the geometric term.Combinations able to reduce this uncertainty should be preferred for implementation in an ionosphere-free EKF.The same analysis is then repeated for the integer ambiguity on L1.Again, the smaller the uncertainty on this ambiguity, the larger the probability to correctly fix it.The results are shown in Table 2.In this case, a scale factor is introduced to represent the STD of combined measurements as a fraction of L1 cycles.In addition to this theoretical analysis, the dispersion of actual GRACE data is computed and the relevant STD is calculated (see the last row in Tables 1 and 2) for gaining further insight into the true-world accuracy of the various combinations.
According to the results presented in Tables 1 and 2, the set of measurements including  IF ,  1 , and MW combinations is the one allowing a higher potential accuracy on both the DD geometric term and the cycle ambiguities.This set of measurements is called set "3" in what follows.A potential limitation of set 3 is that cycle ambiguities are present in all the three observation models.Hence, when the number of fixed ambiguities is smaller than three, the observation model provides no useful information to the EKF, which is thus forced to rely only on its propagation model.In this particular situation, the EKF is expected to perform poorly.To cope with this problem an additional set is considered, called set "4" and made of  IF ,  1 , and MW combinations, plus the uncombined  1 measurement.The addition of  1 observations allows the filter to rely on relative positioning information even when less than three cycle ambiguities are fixed.
Because the various combinations have different observation models selecting different sets of measurements affects the observability of the EKF state vector differently.The observability of the system at hand is extremely involved to be analyzed.Indeed, the dynamical system is nonlinear and evolves on a time-varying trajectory.Its linearization is thus time-varying and depends on the state estimate.The observation model is nonlinear, due to the DD geometry terms.Its linearization depends on the state trajectory, and thus it is time-varying as well.Both the dynamical and the observation models depend on an exogenous vector parameter (e.g., the absolute position of one of the receivers) [17].Finally, both the EKF state vector and the measurements depend on the number of DD couples , which is timevarying as well.In such a condition, instead of developing a complicated observability analysis a different approach is pursued, which is based on numerically evaluating the EKF performances using both set 3 and set 4 as shown in the next section.

No-Iono EKF Dynamic Model.
The no-iono EKF state vector is defined as Concerning the process model, the baseline dynamics is described by the nonlinear relative Keplerian model with J2 effects included.This model is equal to the one used in the VTEC-based EKF [17] and is denoted at the time epoch   as follows: where  is the nonlinear baseline propagation function, r  denotes the receiver  position vector in ECEF, Φ is the time-discrete baseline propagation matrix, w  is the baseline process noise vector, and E indicates the mathematical expectation.
Cycle ambiguities are instead assumed to be constant.However, a small process noise w  might be injected, letting them be modeled as a random constant plus random walk.The model is equal to the one used in the VTEC-based EKF [17] and is denoted by (1  is the -by- identity matrix)

No-Iono EKF Observation
Model.The two different observation models share some common structure, which is discussed first.Then, details on each of the two models will be given.The no-iono EKF measurement vector  is made of a series of measurement combinations, depending on the chosen set.Each measurement combination stacks together the combinations for each of the  DD couples: The covariance matrix of  has a nonsparse structure, due to the presence of the various uncombined measurements in more combinations, and because of the full covariance matrix of each of the uncombined measurements X.In order to have the desired covariance matrix, cov(), it is possible to introduce the vector  that stacks the uncombined DD observations  1 , . . .,  2 : International Journal of Aerospace Engineering The chosen combined measurements set can be expressed as a linear function of the uncombined measurement vector , by means of a matrix  depending on the chosen set.The matrix cov() can be obtained by cov() applying the linear transformation: The observation model relating the EKF state to the observables is nonlinear due to the geometric term.The geometric term is independent from the chosen measurement set and thus is discussed herein.The  DD geometric terms are related to the baseline by the following equations (  denotes the th GPS satellite position vector in ECEF): The nonlinear observation model and the linear delta observation model, needed for the computation of the Kalman gain in the EKF, are obtained by linearization of the propagated state vector, and it can be expressed generically as Thus, at each time epoch, from the observation equations, the measurement vector  is computed combining uncombined measurements.The choice of the measurement set to consider thus changes only the following: (1)  matrix for obtaining the EKF combined measurement vector  from the uncombined one , (2) ℎ() nonlinear observation function, and (3)  matrix of the linearized observation model.These items are discussed in the following subsections for each measurement set.

Set 3𝑛.
The  matrix, the ℎ() function, and the  matrix take the following form: where    vector stacks together  DD geometric terms, ∇ b is the gradient operator with respect to the baseline vector, and 0  is the -by- null matrix.

Set 4𝑛.
The  matrix, the ℎ() function, and the  matrix take the following form: where the subscript "3" is used to indicate matrix and functions relevant to set 3.It is worth noting that, since the ionospheric delay is cancelled and not included in the state vector, the implemented observation model for  1 measurement is Comparing ( 4) to (26), it is evident that set 4 implicitly assumes the DD ionospheric delay modelled as an additional Gaussian white noise affecting the  1 measurements.In other words, considering that the DD ionospheric delay is only a small fraction of the total ionospheric delays, the DD  1 observable noise is modelled as if the ionospheric delay was not time correlated but included within the white noise error.This approach resembles that one used in [7,14] where ionospheric delays are modelled as random walk processes characterized by large process noises and correlation times in the order of the filter time step.

Results
The reference EKF and the no-iono solutions have been run on GRACE data.Comparisons among the solutions are performed looking at the baseline magnitude estimation error and at the capability to obtain correct estimates of the integer ambiguities.Both of them are estimated thanks to the high accuracy reference solutions, which can be computed processing GRACE data [8].As expected, application to mild and moderate ionospheric conditions, not shown herein for brevity, suggest that the filter with Lear's model is capable to achieve a baseline estimation error which is the half of the no-iono approach (both set 3 and set 4).This behaviour changes under strong ionospheric activity.In order to show this, DOY286 (October 13th) 2011 has been selected for the analysis, in which the two GRACE spacecrafts are separated of more than 200 km, and ionospheric activity is extremely intense.Figures 1-3 report the baseline magnitude estimation error and the IA estimation performance for the reference solution and the no-iono solutions, respectively.Specifically, only L1 estimation performance is shown because wide lane fail rates are significantly lower than those on L1. Figure 1 shows that, under strong ionospheric activity, Lear's model is not able to support the integer estimation, so a significant number of ambiguities is wrong and baseline magnitude estimation errors larger than 1 m occur.IA estimation performance increases notably when set 3 is considered (see Figure 2), and the baseline estimation error reduces to decimetre scale in this case.Finally, with reference to set 4, it is evident that the addition of the  1 measurements is not a successful strategy when the ionospheric activity is intense.Indeed, DD ionospheric delays dominate over  1 measurement noise and so they cannot be modelled as additional Gaussian white noise.

Concluding Remarks
This paper has compared various strategies for compensating ionospheric delays in long baseline relative positioning applications for LEO spacecraft using dual-frequency CDGPS techniques.Based on a relative positioning technique previously developed by the authors, problem settings that do not require ionospheric delays estimation have been investigated.The ionosphere-free alternatives have been also  compared with a reference solution making use of a standard ionospheric delay model for spaceborne applications, originally developed by Lear.Results suggest that, among the possible alternatives and under strong ionospheric activity, best performance is achieved using the observation model International Journal of Aerospace Engineering made of ionosphere-free combinations of carrier-phase measurements, GRAPHIC combinations on the L1 frequency, and Melbourne-Wubbena combinations.Future work will deal with the extension of these results to a more complete data set, including different levels of the ionospheric activity.The ultimate goal is the introduction of quantitative tools to guide the selection between model-based and ionospherefree approaches depending on the ionospheric activity and the available ionospheric models.
=    −    −    +    is the DD geometric term between the two receivers; International Journal of Aerospace Engineering 3 (iv)    =    −    −    +    is the DD ionospheric delay on the L1 frequency, which is denoted simply by ionospheric delay in the following;

Figure 1 :
Figure 1: Baseline magnitude estimation error and L1 ambiguities estimation performance for the reference solution.

Figure 2 :
Figure 2: Baseline magnitude estimation error and L1 ambiguities estimation performance for set 3.

Figure 3 :
Figure 3: Baseline magnitude estimation error and L1 ambiguities estimation performance for set 4.

Table 1 :
Expected and estimated standard deviations of the various combinations for   estimation (NA: not available).

Table 2 :
Expected and estimated standard deviations of the various combinations for  1  estimation (NA: not available).