EFFECT OF SPATIAL CORRELATION ON THE PERFORMANCES OF MODERNIZED GPS AND GALILEO IN RELATIVE POSITIONING

In the context of processing GNSS (Global Navigation Satellite System) data, it is known that the estimation of 
 the ionospheric delays decreases the strength of the observation model and makes significant the time required to fix the ambiguities 
 namely in case of long baselines. However, considering the double-differenced (DD) ionospheric delays as stochastic quantities, the 
 redundancy in this case increases and leads to the reduction of time of fixing the ambiguities. The approach developed in the present 
 paper makes two considerations: 1) the DD ionospheric delays are assumed as stochastic quantities and, 2) the spatial correlation of 
 errors is accounted for based on a simple model of correlation. A simulation is made and aims to study the effect of these two mentioned 
 considerations on the performances of the three multifrequency GNSSs; modernized GPS, Galileo and BDS which are not yet in full capability. 
 For each GNSS, dual-frequency combinations of frequencies as well as triple-frequency combination are investigated in the simulation. 
 The performances studied include: the time to fix the ambiguities with high success rate and the precision of coordinates in static 
 relative positioning with varying baseline length. A method is developed to derive what we call the spatial correlation model which 
 approximately gives the covariance between the individual errors belonging to two stations. Furthermore, the stochastic models that 
 follow from accounting and neglecting the spatial correlation are developed. The LAMBDA (Least-squares Ambiguity Decorrelation Adjustment) 
 method is implemented for ambiguity decorrelation. The results show that the time to fix the ambiguities caused by accounting the spatial 
 correlation is less than the time of fix without the spatial correlation. Also, a slight superiority of Galileo in terms of performances 
 is seen compared to the other GNSS. For all the dualfrequency combinations investigated, when processing a baseline length of 500 km with 
 accounted spatial correlation, the time needed to successfully fix the ambiguities lies between 5 and 9 min, whereas it becomes only 
 between 2.5 and 3 min for all the triple-frequency combinations, this is with a sampling time of 5 s. In addition, for all different 
 combinations, the coordinates precision is less than 8 mm even for 500 km. We think that these high performances result from: 1) the precise 
 codes of future GNSS signals, 2) the high redundancy in the observations equation and, 3) taking into account the spatial correlation 
 in the definition of the stochastic model.


Introduction
In the ongoing modernization of GPS, the new signal L5 available on the block IIF satellites which have been launched since 2010 is introduced and added to the current L1 and L2 signals. Besides this, in the upcoming few years, Galileo and BDS will be fully operational. Galileo was designed to transmit for Open Service, signals centered at E1, E5b and E5a such that the code precision on L5a is likely to be better than E1. Similarly, BDS will transmit to GNSS users three signals: B1, B2 and B3. Two frequencies are shared between GPS and Galileo: L1 (E1) and L5 (E5a), and one frequency is shared between Galileo and BDS: E5b (B2). Unlike the three GNSS cited strong and mitigates consequently these errors to a great extent. However, for medium and long baselines, the spatial correlation becomes weak and decreases as the length of baseline increases. In this case the spatially correlated errors remain and need to be either estimated or considered as residuals (i.e. resulting from double differencing) of stochastic nature. The estimation of ionospheric delays decreases the redundancy and results in slow convergence of the estimated coordinates to cm-mm precision. In other part, considering the residuals as stochastic quantities with known statistics, seems to be beneficial in the sense that it augments the redundancy. Consequently, the ambiguities are expected to be fixed with high probability and within short time span especially for long baselines.
The present study aims to examine through a simulation, the effect of spatial correlation on the time required to achieve a successful fixing of the ambiguities as well as on the position quality with respect to baseline lengths.
The subsequent subsections will give in detail the steps followed to reach our purposes. This simulation does not need real or simulated observations since the performances that we wish to evaluate by least squares method need only the design matrix and the variance covariance matrix of observations. This work is organized as follows: first, we start by deriving the spatial correlation model. Second, we present the observations equation and the stochastic model along with the assumptions made in the operation of derivation. In the next step, the parameters used in the simulation as well as the different combinations of frequencies for each GNSS are provided. Also, a brief review of how the ambiguities are decorrelated by LAMBDA method is given. The numerical results of the simulation are presented and analyzed and in the last section we give conclusions.

Spatial correlation model
In this subsection, we give a derivation of a simple mathematical model of spatial correlation that approximates the covariance between individual errors belonging to two stations. This model is derived from the known variances of undifferenced (UD) and DD errors. In addition, the effect of accounting or neglecting the spatial correlation on the expression of the variance covariance matrix of DD errors will be demonstrated.
The fact that certain errors are mitigated in double-differencing, this proves the presence of non-null covariance between them. We try in this section to model approximately this covariance, following a matrix fashion in the derivation and taking the orbit error as typical example, as it will be described below.
First, we apply the covariance propagation law on the vector of UD orbital errors is a matrix of order m -1 whose the diagonal elements are 2 and ones elsewhere. We know that each diagonal element in the matrix DD Q is the variance of DD error, i.e.
Therefore, from Equation (7), the covariance between UD orbital errors 1,2 ( ) σ belonging to station 1 and 2 for any satellite, can be evaluated approximately by The covariance 1,2 σ is here positive since .
DD σ > σ The simple model (8) is what we call the spatial correlation model, it means that if we know the variance of the UD and DD error, we can derive approximately the value of the covariance between individual errors. This covariance is usually required when attempting to develop a stochastic model for relative positioning.
Using the spatial correlation model, the expression of DD Q of equation (6) becomes Now, when the spatial correlation is neglected, the matrix DD Q will be diagonal and takes the form In this study, the same form of the variance covariance matrix (9) (with spatial correlation) and (10) (without spatial correlation) that have been designed for orbital error is also adopted for ionospheric delay. Note that the variance 2 DD σ is assumed constant regardless of the satellite elevation angle.

Observations equation
In relative positioning, double-differenced code and carrier phase observables are usually formed because many systematic errors existing in GPS measurements cancel out or decrease in size. This produces on one part, a simplified mathematical model and ensure on the other part an integer nature of ambiguities. When the integer ambiguities are parameterized within a strong model and then fixed by an appropriate ambiguity resolution technique, this allows to achieve a fast and precise positioning.
The linearized observations equation at a single epoch k for a multi-frequency GNSS, is written as where ⊗ denotes the Kronecker product (see Schäcke, 2013;Zhang & Ding, 2013) for more information about the properties of Kronecker products). , k k R Φ denotes respectively the vector of observed minus computed DD code and carrier phase measurements. f is the number of frequencies used.
[1, 1, ..., 1] T f e = a column vector of length f. The matrix k F is a concatenation of the Jacobian matrix that follows from the derivation of DD geometric distances with respect to the unknown station coordinates, and the column vector of Niell Mapping Functions (NMFs) of the tropospheric wet component. an identity matrix of order m -1 such that m is the number of satellites. The vector (4 1) g × comprises the unknown geometric terms: station coordinates and relative zenith tropospheric wet component. a is the vector of DD integer ambiguities. The vectors , are the sum of unmodeled DD errors for code and phase respectively, such that where the DD error vectors (1) , i δρ and (2) i denote respectively: the orbital error, the first and second order ionospheric delay. Any other systematic error in equation (11) that can exist is considered properly corrected. At each epoch k, the matrix k F in equation (11) needs to be updated as the satellite positions change with time. The ionospheric delay in equation (12) has been approximated to a second order according to Datta-Barua et al. (2008) and Alizadeh et al. (2013).

Stochastic model
The stochastic model describes the precision of measurements and the eventual covariances that exist between them. In this section we give one stochastic model in general form that can reflect the two cases to be investigated: when the spatial correlation is accounted (using Equation (9)) and when it is neglected (using Equation (10)). By applying the covariance propagation law on the error vector where the blocks with: nal matrix containing the variances of the DD orbital error, first and second order ionospheric delay. The matrix is as defined in section 2 and 1 m C I − = if it is neglected. The diagonal matrices R D and D Φ are of order f that contain respectively the variances of UD code noise and carrier phase noise at zenith (their magnitudes are provided in Table 1 below). Ω is the double-difference operator matrix. k M denotes the diagonal matrix of noise mapping functions at epoch k.
It is worth to notice that during the development of the stochastic model, the following considerations and simplifications are made: -We consider no correlation between different types of errors since they come from independent physical processes. For example, there is no need to consider that the ionosphere and orbital effect are correlated; -The temporal correlation of noise is neglected as at least a sample time of 5 s is chosen in the simulation, according to (Nardo et al., 2015); -The temporal correlations of ionospheric residuals are neglected in the stochastic model to avoid inverting a large and fully populated variance covariance matrix of observations; -The uncorrected DD ionospheric delays are normally distributed and centered at zero even for long baselines (up to 407 km) as stated (Liu, 2001, p. 100); -No cross-correlation between signals exists, since we assume that the signals in either modernized GPS, Galileo or BDS will be likely non-correlated. -The mapping function for code noise is assumed identical to carrier phase noise, it is formulated as where ε is the elevation angle in degree. Note that we computed the parameters of this mapping function by fitting a curve of exponential function 0 -/ ( ) a be ε ε σ ε = + through the observed RMS of noise. These observed values of RMS result from the study conducted by (Cai et al., 2016), for a Trimble NetR9 receiver are exploited for accurate assessment.

Simulation parameters
The positions of GPS satellites are computed from almanac data using YUMA files for May 3, 2018, from 19:39 to 20:32. The full GPS constellation is assumed to contain 31 satellites. Because the Galileo and the MEO (Medium Earth Orbit) of BDS satellites are not yet fully deployed, the positions of satellites are simulated based only on their nominal constellation configurations described as follows. For Galileo: 3 orbital planes, 9 satellites equally distributed and one spare satellite in each orbital plane, nearly circular orbit, inclination: 56° and semimajor axis of 29,600 km. For BDS, 3 orbital planes, 8 satellites equally distributed and one spare satellite in each orbital plane, nearly circular orbit, inclination: 55° and semi-major axis of 27,900 km. The full Galileo and BDS constellations are each assumed to contain 27 satellites. The period of time is chosen such that the same set of 10 satellites still in view in each GNSS. The reference station is located at latitude 34° 40′ N and longitude 3° 14′ E (located at Djelfa, Algeria). The second station is located in the west of the reference station and has the same latitude. The elevation mask angle is set at 10°.
The variation of standard deviation of DD orbital errors ( ) DDδρ σ with respect to the length of baseline b can be formulated as where q takes the values of 4.9, 4.6 and 4.3 for GPS, BDS and Galileo respectively. δρ σ the standard deviation of the UD ephemerides. See in the appendix how the formula (18) is derived. For GPS, δρ σ is set to a value of 1 m as it corresponds to the IGS (International GNSS Service) broadcast ephemerides precision. Similarly, we assume the same magnitude of δρ σ for the Galileo and BDS. Concerning the ionospheric errors, the standard deviations used in the simulation are taken from (Feng 2008) and provided in Table 1. In this table, we see that the first order DD ionospheric delay standard deviation grow with distance by about 2 ppm. However, in the study conducted by (Liu & Lachappelle, 2002), it is only 1.5ppm for an active day and without any ionospheric corrections. Therefore, the values of standard deviation of Table 1 have been likely obtained without correcting the ionospheric delays. For a given baseline length that does not correspond to the baselines mentioned in Table 1, we performed a linear interpolation to predict the DD Ionospheric standard deviation.  Deprez and Warnant (2016) (Tables 1, 2 and 3) provide the average standard deviations of noise corresponding to a Trimble NetR9 receiver for GPS, Galileo and BDS. Since we need these values at zenith, we divide the average standard deviations by a factor of 1.5 to obtain approximately the standard deviations of noise at zenith. The following Table results from this operation.
In the previous two subsections we presented a general form of the observations equation and the associated stochastic model designed for a multi-frequency combination of signals. For this study, dual-frequency combinations as well as triple-frequency combinations are planned for comparison purpose. These combinations are summarized in the following Table 3.

Ambiguity decorrelation
By least-squares adjustment, we can obtain the float solution g a         and its associated variance covariance matrix After obtaining the float solution, the ambiguities need to be fixed on the correct integers: this step is called: ambiguity resolution. An advantage of this technique is its ability to improve the precision of other estimates of interest as for example the baseline components. In this work, the ambiguity resolution method LAMBDA (Least-squares AMbiguity Decorrelation Adjustment) is implemented. This method was first introduced in 1993 by P. J. G Teunissen in his paper (Teunissen, 1993) and discussed in detail in Teunissen (1995), it is known among other methods by its high success probability of resolution, see the advantages of LAMBDA method over other methods in Teunissen et al. (2002) and the study conducted by Jazaeri et al. (2014) that compares different decorrelation techniques. The LAMBDA method is performed into two steps: the first step achieves a decorrelation of ambiguities i.e. reduction of correlation of ambiguities by using Z-transformation and the second step consists of finding the most probable integer candidate inside the so-called: ambiguity search space. The role of the first step is to speed up the second step and to make the search process more efficient (Chang et al., 2005). The Z-transformation consists to find the matrix Z, such that Matrix Z should obey three conditions: 1), it should minimize the product of diagonal elements of the transformed covariance matrix Z Q  , 2): the entries of Z are integers and 3): det( ) 1, Z = ± see Teunissen (1994) for more details. L is a unit lower triangular matrix. D diagonal matrix containing the so-called conditional variances.
Once the ambiguities are fixed with high success probability (i.e. when a probability very close to 1 is obtained), the variance covariance matrix of the fixed solution g  reads: -1 , , -.
To measure the success of ambiguity resolution, usually the so-called success probability or success rate is computed. The success rate that corresponds to the search technique ILS (Integer Least Squares) is difficult to be computed but the computation of success rate related to the search technique called bootstrapping is possible. The bootstrapped success rate denoted s P is regarded as the lower bound of ILS success rate and is computed by where n -number of ambiguities; i σ -the standard deviation of conditional ambiguity i that appear in matrix D, for further details see Teunissen (1998), Chang et al. (2005), Verhagen and Li (2012). Arora (2012) states that the ambiguities were considered successfully resolved only after a minimum of 0.999 (or 99.9%) of probability is obtained from integer bootstrap. This threshold probability is chosen in the study conducted by Odijk et al. (2014) and is adopted also here. In this work, the time to fix the ambiguities for a given baseline length is determined by accumulating the data until a success rate greater than or equals the threshold 99.9% is attained.

Results and discussion
The present section provides the numerical results relative to several combinations of frequencies for each GNSS and gives an insight of their performances corresponding to a certain period of time. The GNSS performances that we focus on in the simulation are: the time required to fix the ambiguities and the precision of station coordinates with respect to baseline lengths. The effect of non-considering the spatial correlation on the performances is highlighted only in the Galileo plots. Figure 1 below shows the plots of the time required to fix the ambiguities against the baseline length with a 5 s sampling time. In Galileo plot, it is clear that when the spatial correlation is considered (solid lines), this causes a shorter time of fix compared to when it is not considered (dashed lines). All the Galileo frequency combinations (solid lines) do better than their counterparts in GPS and BDS, this is due to the high precisions of the codes of Galileo compared to those in the other GNSSs as shows Table 2. BDS also do better than GPS in terms of time of fix of ambiguities.
For all GNSS, we can see that the triple-frequency combination in general outperforms the dual-frequency combination in terms of time to fix the ambiguities. The time of fix that results from all dual-frequency combinations, lies between 5 and 9 min for 500 km baseline, whereas, for triple-frequency combinations, this time lies only between 2.5 and about 3 min with a superiority of Galileo (2.5 min). Also, we can remark in the GPS plot that the combination (L1+L2) provides a longer time of fix compared to (L1+L5), though, the code precision of L2 and L5 are almost the same (20 cm). Based on this remark, we can deduce a rule for the dual-frequency combinations. This rule is stated as follows: as the two frequencies get closer, the time of fixing the ambiguities get longer and vice-versa. In Galileo plot for example, where the combination (E1+E5b) gives a longer time than (E1+E5a), this effect does not result from only the fact that the code E5b is less precise than E5a, but also because the frequencies (E1, E5b) are closer to each other than (E1, E5b). Note that, the combination of the second and third frequency, which are the closest frequencies in each GNSS, has not been visualized in the plots because of their longer time of fix. For instance, the combination (E5b+E5a) yields for 500 km, a time of 19 min even though E5b and E5a have the highest code precision. Such kind of combination may be seen as a reinforcement of the correctness of the rule deduced above.
The precisions of fixed coordinates, arising from the different combinations are depicted in Figure 2. Note that the precisions in this figure are computed in analogy to the ADOP (Ambiguity Dilution Of Precision) concept. For more information about ADOP and its properties see (Teunissen and Odijk 1997). Following this concept, the precision of coordinates is represented uniquely by one measure Q  a 3×3 variance covariance matrix of baseline components. The advantage of computing the precision by such formula is its ability to use the complete information inside the matrix b Q  i.e. considering the off-diagonal elements instead of considering only the individual diagonal elements of b Q  as precision indicators. In Figure 2, we see for Galileo case, that there is no significant difference in positioning precision between accounting and no accounting, the spatial correlation.
The precisions do not exceed 8 mm for GPS, 5.5 mm for Galileo and 5 mm for BDS. The dual-frequency combinations in each GNSS give relatively the same magnitude of precision regardless of the variation of the baseline length. Unlike the GPS case, the dual-frequency combinations in Galileo or BDS do better than the triple-frequency combination. We can also remark in Galileo and BDS at about 200 km, a progressive decrease in precision with baseline length. The progressive decrease in precision that starts at about 200 km can be explained by the fact that the decrease of precision of ionospheric delays with increasing baseline length, starts to outweigh at 200 km, the gain of precision of estimates caused by the increasing number of observations. Additionally, and despite of the ADOP-based precision results, the horizontal precision is evaluated in general to less than few millimeters and the height precision to less than few centimeters for all combinations studied.
We think that the resulting short time to fix the ambiguities and the high positioning precisions are essentially due to: 1) the high redundancy caused by the non-parameterization of the ionospheric delays, 2) the precise code used and third, taking into account the spatial correlation in the stochastic model. This study is based on a simple model of spatial correlation, thought, it needs more refinements and improvements. Despite of whether the spatial correlation model established is realistic or not, we believe that accounting for the spatial correlation in GNSS data processing will have as consequences: a shorter time span to fix the ambiguities and an improved positioning precision.

Conclusions
In the present work we consider the effect of the spatial correlation of errors in studying the performances of the modernized GPS, Galileo and BDS. The performances investigated here are: the time to fix the ambiguities and the positioning precision with varying baseline length. From the given variances of DD errors that depend on baseline length, we proved in matrix fashion the presence of a positive covariance between the individual errors and then we derived what we call the spatial correlation model. Also, when neglecting the spatial correlation, the resulting variance covariance matrix of DD errors is provided. The stochastic model is developed and written in general form to include the two cases: accounting or neglecting the spatial correlation. Furthermore, for each GNSS, several combinations of frequencies are planned to show the response of these combinations to the spatial correlation. The ambiguities are decorrelated using LAMBDA method and they are considered successfully resolved once a 99.9% success rate is attained. The results that is based on simulated data, show that accounting the spatial correlation yields a shorter time to fix the ambiguities compared to when the spatial correlation is neglected in the stochastic model. Galileo system is slightly superior in terms of performances, compared to the two other GNSSs. For all the dual-frequency combinations studied, when a long baseline of 500 km is processed with a 5 s interval, the time required to successfully fix the ambiguities lies in general between 5 and 9 min, whereas, it lies only between 2.5 and 3 min for all the triple-frequency combinations. The ADOP-based precision of positioning is high and does not exceed 8 mm for all combinations even for 500 km baseline. We believe that these considerable performances result from: 1) the high redundancy, 2) the precise codes of future GNSS signals and 3) accounting for the spatial correlation in the stochastic model.
We provide in this appendix an approximation of the differenced orbital error or rather providing the maximum value of the differenced orbital error attained with respect to the baseline length. In the Figure 3 below; S and 0 S designate the true and the computed (broadcast) position of a satellite.
(3 1) δρ × is the vector of orbit error, note that S and 0 S are not in the same plane. θ is the angle between the two receivers. 1 2 , ρ ρ are respectively the distances between receivers and 0 S . b denote the baseline. The unit vectors 1 2 , u u point from the receivers towards 0 S . An orbital error in the direction receiver-satellite is the projection of the vector δρ onto this direction. Following this rule, the error vectors in their respective directions receiver-satellite are expressed as: var( ) δρ is the (3×3) covariance matrix of orbit error. If var( ) δρ is approximated to 2 I σ where σ denotes the average standard deviation of UD orbital error and I an identity matrix, then 2 ∆δρ σ will be 2 2 2 2 1 2 1 2 1 2 2 2 1 2 2 2 1 2 ( -)( -) 2 (1-) -2 (1-cos ) 2 1-2 by using the trigonometric theorem of Al-Kashi, i.e. In the case of GPS, for example, the distance receiversatellite ρ varies between 20,200 km (satellite at zenith) and 25,800 km (satellite on the horizon). Now, by adopting a minimal value of ρ , the maximum of the SD orbital error max( ) ∆δρ σ reaches a great amount. Based on the considerations stated above, the standard deviations for will be at worst as note that max( ) ∇δρ σ and σ are in unit of meter and the distance b is in km.