Research Article  Open Access
Instantaneous TripleFrequency GPS CycleSlip Detection and Repair
Abstract
A realtime algorithm to detect, determine, and validate the cycleslips for triplefrequency GPS is proposed. The cycleslip detection is implemented by simultaneously applying two geometryfree phase combinations in order to detect more insensitive cycleslips, and it is applicable for high data rate applications. The cycleslip determination adaptively uses the predicted phase data and the code data. LAMBDA technique is applied to search for the cycleslip candidates. The cycleslip validation provides strict test criteria to identify the cycleslip candidates under low phase noise. The reliability of the proposed algorithms is tested in different simulated scenarios.
1. Introduction
Failure of the GPS receiver, signal interruption, low signaltonoise ratio, and high receiver dynamics can cause a sudden jump of a specific integer number of cycles in the carrier phase measurement, which is referred to as cycleslip. Cycleslips may occur independently on each carrier frequency per GPS satellite and remain in the phase data.
The handling of cycleslips is conventionally composed of four sequential stages: ( 1) cycleslip detection, which checks the occurrence of cycleslips, ( 2) cycleslip determination, which quantifies the sizes of cycleslips, ( 3) cycleslip validation, which tests whether the cycleslips are correctly resolved, and ( 4) cycleslip removal, which removes the cycleslips from the phase measurement.
The occurrence of cycleslips is a random event, and therefore the cycleslip detection should be applied epochbyepoch. For this reason, it should be a rapid algorithm with small computational burden. Cycleslip determination and validation will be performed for the phase data contaminated by the cycleslips. After the cycleslip values are fixed and pass the validation, the final stage, namely, cycleslip removal, can be simply realized by mathematical subtraction.
The processing of cycleslips can be either applied to standalone positioning with only one receiver used or to differential positioning (DGPS) where two or more receivers are involved. In each case, the cycleslip detection and determination can be further categorized according to the number of available signals. Listed below are the commonly used cycleslip detection and determination methods, including: (a) polynomial fitting [1], (b) a high order betweenepoch phase differences [2], (c) Kalman filter prediction [3], (d) phase combinations [4, 5], (e) phase/code combination [2], and (f) quality control [6, 7]. In Table 1, these methods are categorized according to the application scenarios.

The methods (a) and (b) are based on the phase data of the previous epochs, and normally 4 or more observation epochs are needed. Therefore, they are applicable in static or low dynamic applications, where the phase data show a smoothed curve when plotted with time. But if the antenna is undergoing a complex motion, the phase is hard to be predicted, and therefore these methods might fail. In comparison with the singlefrequency approach, the dualfrequency methods are superior as they do not rely on the motion of the antenna and are very sensitive to detect the small cycleslips. There are, however, some specific cycleslip pairs which cannot be readily detected by using the dualfrequency phase combination. The detection of these cycleslip pairs is still an issue for dualfrequency GPS.
Due to the introduction of the new GPS signal, the traditional approaches dealing with cycleslip problems should be expanded to triplefrequency case. Literature regarding triplefrequency cycleslip detection and correction is still scarce, and this motivates us to make an investigation in this field. As the algorithms suitable for standalone positioning can also be applied to differentialpositioning, our approach is then oriented to the single GPS receiver in order to have a wide range of application scenarios. For cycleslip detection, we mainly refine the approach based on the geometryfree phase combination [2] and expand it to the triplefrequency case. Major modifications are implemented based on the fact that there are more geometryfree combinations in case of triplefrequency signals. The number of insensitive cycleslips is decreased with respect to the dualfrequency case by properly choosing two geometryfree combinations. In comparison with the singlefrequency method, the proposed triplefrequency approach also has the advantage that it better fits dynamic applications. Considering that the cycleslips and carrier phase ambiguities have a similar integer nature, we apply the LAMBDA technique to search for the cycleslip candidates in the cycleslip determination. Besides that, an efficient cycleslip validation approach is also introduced.
In Sections 2, 3, and 4, the cycleslip detection, validation, and determination will be presented, respectively. In Section 5, simulations are carried out to test the performance of the algorithms in various scenarios.
2. CycleSlip Detection
2.1. CycleSlip Detection Model
We first define the cycleslips occurred on the triplefrequency signal of a GPS satellite as cycleslip groups , where denotes the cycleslip value on the signal.
The carrier phase observation equation for each signal at epoch can be expressed as
where (subscript) indicates the corresponding signal; is the carrier phase observable in cycles; is the geometric distance from the GPS receiver's antenna phase center at the epoch of signal reception to the GPS satellite's antenna phase center at the epoch of signal transmission; is the integer ambiguity in units of cycles; is the wavelength; is the ionospheric delay; is the tropospheric delay; is the satellite orbit bias; is the satellite clock bias; is the receiver clock bias; is the thermal noise contained in the carrier phase data in units of cycles; is the multipath error. All terms except for , , and are given in units of length.
Suppose that a cycleslip arises on the signal at the next epoch . We have
By differencing (1) and (2) the betweenepoch observation equation is obtained as
where the operator represents the differencing between epochs and .
For triplefrequency GPS signals, a linear combination of betweenepoch observation equations can be expressed as
where , , and are scalars. The core idea of cycleslip detection is to derive the relation between the phase measurement ( terms) and the cycleslip values ( terms). For this purpose, the nondispersive error term and the dispersive error term need to be cancelled from (4). The term can be eliminated by assuring that the sum of the scalars , , and equals zero, that is, by forming a geometryfree phase combination. The term is difficult to be eliminated but can be possibly minimized by a proper choice of the scalars. The selection of the scalars will be discussed later, and we will first introduce the cycleslip detection principle.
Assuming that the term is eliminated, and neglecting the betweenepoch ionospheric error and multipath error in , (4) can be rewritten as
The carrier phase noise is assumed to be white Gaussian noise. It is further assumed that the signals have the same resolution in units of cycles, that is, [8]. Applying the variance propagation law yields the noise of the lefthand side of (5), expressed by :
where reflects the betweenepoch differencing.
By choosing a proper confidence level, we have a critical value to test the occurrence of the cycleslips, where the scalar denotes the multiplication factor of and is always chosen as (% confidence level) or (% confidence level) in GPS applications. Based on the prior discussions, once the following inequality holds true, we can conclude that a cycleslip group arises on this satellite:
2.2. The Construction of Optimized GeometryFree Combinations
For triplefrequency GPS, there are theoretically an infinite number of scalar groups which can be used to form the geometryfree combinations, and therefore we need to give some constraints for the choice of the scalars. Firstly, the scalars can be either integer values or float values. However, the float values can be converted to their integer counterparts through the following relation:
As the term will exist on both sides of (4), it will not affect the cycleslip detection. In this sense, the float scalars are equivalent to the integer scalars, we therefore only consider the integer scalars. The second case to be avoided arises when there is an integer common divisor contained in each scalar group , because the common divisor can be canceled simultaneously from the both sides of (4). Summarizing the discussions above, we should only consider the integer scalars having no common divisor. For example, as the effect of the scalar group is equivalent with and , we only take into account.
In order to clearly interpret the minimization of the term , a rigorous expression of (7) should be formulated first by taking the term into account:
Formula (9) reveals that the term to be minimized is actually not itself but divided by the term . After substituting the expression of given in (4) into formula (9), we can start discussing the minimization of the ionosphere term, the thermal noise term, and the multipath term, respectively. For simplicity, we use to denote in the following.
The reduction of thermal noise can be formulated as
The magnitude of thermal noise depends on the carrier signals and varies with time. As we have assumed the thermal noise to be white Gaussian noise, the term can be replaced by the corresponding standard deviations. According to the prior assumptions that each signal has the same phase resolution, the relation
should hold true, implying that should be assigned with possibly small values.
The ionospheric delay can be minimized through the following relation:
The betweenepoch ionospheric error on a single carrier signal is usually at millimeter level between two closespaced epochs with high data rate [9], and furthermore, the selected phase combination can guarantee that the combined ionospheric error is ignorable.
The multipath error reduction is more difficult in comparison with the thermal noise and ionospheric error reduction, since the multipath error depends on the environment nearby and generally cannot be modeled as a Gaussian distribution. For these reasons, we will not explore the effect of multipath on the cycleslip detection mathematically. The performance of cycleslip detection under high multipath environment will be shown later using simulations.
According to the constraints given before, we can choose the scalars to construct the geometryfree combinations. The thermal noise reduction requires that should be small, and therefore we fix the search range of each from to cycles. Within this range, the scalars presented in Table 2 yield the geometryfree combinations having relatively small ionospheric residuals, where the “ionospheric residuals” represents the result of the lefthand side of (12).

It should be stressed that the scalars should be nonzero values. A zerovalued scalar implies that the cycleslip detection on the corresponding signal is excluded.
2.3. The Selection of Two Optimized GeometryFree Combinations
The geometryfree combinations formed by the scalars given in Table 2 can be used to detect the cycleslips. However, there are some special cycleslip groups which cannot be readily detected using a single geometryfree phase combination. We define them as insensitive cycleslip groups, which usually fulfill the following inequality:
It can be understood that the insensitive cycleslip groups are related to the scalars , and hence there exist different insensitive cycleslip groups depending on the scalars used. Nevertheless, the following cycleslip groups are always undetectable as they are proportional to their individual frequencies:
We define the cycleslip groups in (14) as the most insensitive cycleslip groups.
The insensitive cycleslip groups ranging from to cycles belonging to the combination are given in Table 3. Some statistical results of the insensitive cycleslip groups within different ranges are presented in Table 4, where the occurrence probability means the ratio of the insensitive cycleslip groups to all cycleslip groups.


In case of dualfrequency GPS, these insensitive cycleslips are usually ignored due to their low probability. In triplefrequency GPS, the number of insensitive cycleslips can be reduced due to the fact that the insensitive cycleslip groups corresponding to a specific geometryfree combination could be detected by applying other geometryfree combinations. Nevertheless, it is not computationally efficient to apply many geometryfree combinations simultaneously for cycleslip detection. Our goal is then to detect maximal number of cycleslip groups using minimal number of geometryfree combinations. From the study we found out that all the cycleslip groups except for the most insensitive ones can be detected by properly choosing two geometryfree combinations from Table 2. These dualcombinations are presented in Table 5, where the last column shows the sum of the ionospheric residuals.

By comparison, we will use the geometryfree combinations constructed by the scalars and simultaneously for cycleslip detection, because they contribute the smallest ionospheric residuals. We define them as the first optimal phase combination and the second optimal phase combination, respectively.
2.4. TCAR for CycleSlip Detection
The Three Carrier Ambiguity Resolution (TCAR), as described in [10], sequentially applies the dualfrequency ambiguity resolution technique based on geometryfree phase combinations. Due to the use of the geometryfree combinations, TCAR is also applicable for cycleslip detection in triplefrequency GPS. The basic idea is then to sequentially detect the cycleslips on each combined dualfrequency signal. The procedure canalso be expressed using the scalar groups as shown in Table 6.
 
Number “1” implies that only the most insensitive cycleslip group is undetectable. 
Note that the zerovalued scalar in each scalar group makes the triplefrequency phase combination be equivalent to the dualfrequency combination. The first three rows contain only the dualfrequency phase combinations, whereas each row of the last three rows is composed of a dualfrequency combination and a triplefrequency combination. From the last column we know that it is still possible to detect all cycleslip groups except for the most insensitive ones using TCAR. Nevertheless, TCAR contributes larger ionospheric delay in comparison with the optimal combinations and . For this reason, we still use the aforementioned two optimal combinations for cycleslip detection.
3. CycleSlip Validation
Traditionally, the cycleslip validation is the next step following the cycleslip determination. But in this study, the cycleslip validation is embedded into the cycleslip determination to test the cycleslip candidates. Thus we will first introduce the cycleslip validation in this section and then the cycleslip determination in the next section.
The aforementioned cycleslip detection approach can serve as the cycleslip validation approach with a slight modification, as follows: In comparison with formula (9), the difference of (15) lies in the superscript “repair”, which implies that the carrier phase data are already corrected by subtracting the calculated cycleslip values from the original phase data. Once the repaired phase data fail in the test, it means that this cycleslip candidate under test cannot be the true value.
4. CycleSlip Determination
Once the cycleslips on a satellite have been detected, the next step is then to quantify the integer cycles of the slips. A general model for cycleslip determination can be formulated by rewriting the betweenepoch phase observation equations expressed in (4):
where the nondispersive errors, including tropospheric delay, satellite, and receiver clock bias, are put into the term since they contribute the same amount to all phase observables; column vector includes the measuring noise, multipath error, and the remaining ionospheric error after being differenced. By comparing (16) with (4) it can be seen that the term is shifted to the lefthand side of the equation as a known value in order to avoid an underdetermined model. However, term is still unknown and should be fixed a priori using additional measurement immune to cycleslips. Then the estimated float cycleslips and the corresponding covariance matrix will be obtained, so that a search space for the integer cycleslip candidates can be defined. The true cycleslip value can be identified by testing all the candidates using the cycleslip validation procedure.
For a single receiver, the following data can be employed to estimate : (a) predicted phase data based on a polynomial fitting, (b) code/phase combination, and (c) Doppler frequency data.
The Doppler frequency data, however, are not available for some receivers. As the Doppler data reflect the instantaneous phase rate between two adjacent epochs, they can also be treated as a first order polynomial fitting. For these reasons, we mainly focus on the first two types of measurement in the following text.
4.1. Using Predicted Phase Data for CycleSlip Determination
The phase data at epoch are assumed to fit an order polynomial expressed by
The coefficients can be estimated using leastsquares principle based on the previous phase data of at least epochs. The residuals of the leastsquares estimation reflect the quality of the polynomial fitting and can be used as the variance of . Using the calculated coefficients, the term is calculated from: where represents the current epoch. The phase prediction based on polynomial fitting is the most commonly used cycleslip detection and determination method for singlefrequency receivers. As stated before, it provides highquality phase estimation in static or low dynamic case and is less affected by multipath error, so that the estimated float value is very close to the true value. But once the antenna is undergoing a complex motion, this method might provide unexpected results. Similar as the predicted phase data, it is also hard to acquire the Doppler frequency data with sufficient accuracy for moving antennas [11].
4.2. Using Combined Code/Phase Data for CycleSlip Determination
Using the code data, the term can be calculated from the following relation:
where the symbol represents the betweenepoch code data on the signal; stands for ; contains the thermal noises of the code data.
The model (19) employs the code data on three signals. Actually, only the code data on a single signal enough to calculate . The can be obtained according to the leastsquares principle once the covariance matrix of code data is assigned a priori.
The resulted from the code data is not affected by the motion status of the antenna, but using code data will cause the problems in the following two aspects. Firstly, the high noise of the code data will result in a less accurate float estimate of the cycleslips and a large search space, somehow degrading the efficiency. Secondly, in the richmultipath environment, the multipath error contained in the code data will bias the estimated float cycleslip values to a large extent. Consequently, the origin of the search space will be significantly deviated from the true value, and the true value might not exist inside the search space.
4.3. The Search of the CycleSlip Candidates
The search of cycleslip candidates can be implemented by applying the LAMBDA technique [12]. Although the LAMBDA technique is related to the ambiguity resolution, it is still reasonable for the cycleslip determination as the cycleslips and the phase ambiguities are both integers. A slight difference lies in the validation of the candidates. The ambiguity candidates are ranked according to their leastsquares residuals, and the best candidate should offer significantly lower leastsquares residuals than the second best one. But in this study, the cycleslip validation can identify the cycleslip candidates rapidly and sensitively, and hence each candidate can be tested by directly invoking the cycleslip validation. In other words, we only use the decorrelation and searching functions of LAMBDA, but not the selection criteria.
Two parameters are needed in order to invoke the LAMBDA technique: the covariance matrix of the cycleslips and the estimated float cycleslip values. In (16), the covariance matrix of known values (denoted as ) and that of cycleslips (denoted as ) read
where is the variance of and depends on the measurement used; is the covariance matrix of the carrier phase data on the three signals; is a threedimensional identity matrix; the factor 2 in the expression of reflects the betweenepoch differencing; the shorthand notations , and are introduced in model (16).
The estimated float cycleslips can be obtained according to the leastsquares principle:
In LAMBDA, the threedimensional cycleslip search space is defined as
where is the float cycleslip estimate on signal, and is the corresponding covariance matrix. The decorrelation of the search space is achieved by iteratively applying the integer approximations of the conditional leastsquares transformations [12]. The change of the cofactor matrix and the motion of the ellipsoid center after a specific number of transformations (denoted as ) are expressed as
where the operator reflects the motion of search space center due to the matrix transformation. The construction of the transformation matrix is discussed in [7]. The fixed true cycleslip value needs to be retransformed by multiplying .
4.4. Summary of CycleSlip Determination
The proposed cycleslip determination is composed of three steps. The first step is to calculate the term using either the predicted phase data or phase/code combination. The second step is then to search for the cycleslip candidates using LAMBDA technique. The third step is to test each the cycleslip candidate using the cycleslip validation.
The predicted phase value and the code/phase combination can be adaptively used to calculate . The first choice is the predicted carrier phase data. If only one candidate passes the validation, we can say that the cycleslip value has been successfully determined. If this is not the case, the code data can be used instead to calculate . If the true cycleslip value is still not resolvable, it can be concluded that the cycleslip determination is not capable of offering reliable cycleslip estimation, and hence the phase data of current epoch cannot be further used for positioning. In this case, the handling of cycleslips should be shifted to the next epoch.
TCAR technique can also be applied for the cycleslip determination. However, LAMBDA will always perform better or at least as good as TCAR, in a sense that it offers the highest probability of success [13]. Thus, we only consider the use of LAMBDA technique in this study.
5. Simulation
Some simulations are carried out to test the algorithms. The code and carrier phase data are generated using the commercial GNSS software simulator SatNav Toolbox for MATLAB by GPSoft. In this software simulator, the errorfree phase and code data are generated first according to the simulated satellite coordinates and the antenna location. Then the errors are added to the phase and code data, where the thermal noise on the phase and code data is produced as white Gaussian noise. The other errors including the ionospheric error, tropospheric error, and multipath error are generated according to the corresponding parameters.
We assume that the antenna is undergoing a movement along a trajectory illustrated in Figure 1. Given in the axes are the coordinates in and direction in an EarthCenteredEarthFixed frame. We assume that the antenna is always moving at the ellipsoidal height of 360 m; that is, value is a constant value. The numbers marked close to the red points represent the epochs. The motion is described in Table 7.

The multipath error generated by the software GPS simulator is depicted in Figures 2 and 3, respectively. Figure 2 shows the multipath error contained in the carrier phase data and the code data in low multipath environment, whereas Figure 3 shows the multipath error in high multipath environment. Both have the same pattern but different magnitude.
(a)
(b)
(a)
(b)
Some parameters need to be specified before establishing the simulation. In formula (7), we set as , meaning that a sigma standard deviation is adopted here. In the cycleslip determination based on the predicted phase data, we use order polynomial, that is, in formula (17). In (19), we assume that code is available. The standard deviation of the code error is related to that of carrier phase by for low code noise and for high code noise, where is the fundamental GPS frequency of MHz, and is the carrier frequency on signal. For example, in case of the low code noise, equals on signal. For C/A code, this value should be enlarged depending on the technology used. In (17), is so set that the search space contains two cycleslip candidates. The sampling rate of the observation is second. Note that in the following figures, the notation “Detection value” represents the result of the lefthand side of formula (7).
We designed different scenarios to test the proposed approaches. A comparison of these scenarios is given in Table 8.

The original phase data are confirmed as cycleslip free, and hence some cycleslip groups should be added to the original carrier phase data to test the algorithm. After the cycleslips are fixed, they will be removed directly from the raw phase data and hence will not affect the following epochs any more. Therefore, even when the cycleslips occur epochbyepoch, the algorithm still works properly.
5.1. CycleSlip Detection in Low Multipath Environment
Small cycleslips ranging from to cycles, that is, from to , have been added to the phase data starting from the epoch with the interval of epochs. Figure 4 shows the detection values of these cycleslip groups using only the first optimal phase combination. Since and cycles as stated before, the threshold in (7) is then cycles.
Except for the cycleslip group at the epoch, the detection values of the other small cycleslip groups exceed thethreshold; therefore this detection algorithm can be regarded as a sensitive algorithm for small cycleslips.
Problems with detection are shown by using the cycleslip groups from Table 3. According to the previous analysis, these cycleslips cannot be detected by using the first optimal phase combination. This also explains why the cycleslip group lies within the thresholds in Figure 4. We add these insensitive cycleslips into the original phase data at the epochs marked in Figure 5. Figure 6 demonstrates the different detection values when detecting cycleslips using the aforementioned two optimal phase combinations.
(a)
(b)
From the Figure 6(b) it can be seen that the detection values from the second optimal phase combination apparently exceed the threshold of . This implies that these cycleslip groups insensitive to the first optimal phase combination have been detected by using the second optimal phase combination.
5.2. Test of CycleSlip Determination in Low Multipath Environment
We use the cycleslip groups in Table 3 again to test the cycleslip determination and add them to the phase data at the epochs given in Figure 5. The results obtained from the phase prediction and from the code/phase combination are given in Table 9, respectively. Multipath error is illustrated in Figure 2.

In Table 9, the notation “None” implies that the cycleslip value cannot been correctly determined, or in other words, none of the results outputted from LAMBDA has passed the cycleslip validation procedure. It can be seen that the method based on the phase prediction fails at the and epochs. At epoch, the antenna is changed from a sinuous motion to a straight line motion. At epoch, the antenna has just made a sharp turn between two straight line motions. Since we use a order polynomial fitting, once the antenna has undergone a significant change of the motion direction in the past epochs, the phase prediction will probably provide a wrong result. However, the incorrectly estimated cycleslip values can be still filtered out by the cycleslip validation.
As discussed before, the cycleslip determination based on the phase/code combinations is independent of the motion status of the antenna. Under the low multipath environment, the multipath error on the code data will not severely bias the float cycleslip estimate, so that the LAMBDA technique will provide the correct integer cycleslip values. The results listed in the last column of Table 9 reveal that these cycleslips are correctly identified since the estimated integer values equal the corresponding true values.
5.3. Test of the CycleSlip Detection in High Multipath Environment
The multipath error on the carrier phase data is several orders of magnitude lower than that contained on the code data. The detection part is mainly carrier phase related and hence less affected by the multipath.
We first use insensitive cycleslips given in Table 3 to check the different detection values in different multipath environments. The detection results using the first optimal phase combination are depicted in Figure 7. Although different multipath errors yield different detection values, these detection values are still limited within the thresholds. Figure 8 shows the detection values when adding small cycleslips ranging from to into the original phase data with the interval of epochs. Comparing it with the results given in Figure 4 we can see similar detection values in low and high multipath environments, and also, all the cycleslips, except for the insensitive cycleslip at the epoch, are detectable. The results obtained from Figures 7 and 8 agree with our analysis that the multipath error does not affect the cycleslip detection significantly.
5.4. Test of the CycleSlip Determination in High Multipath Environment
The method based on the predicted phase data is less affected by the multipath, as only the phase data are involved. Once the phase/code combinations are employed to determine the cycleslips, the term will be severely affected by the large multipath error on the code data. A significant change of multipath error on coda data between closespaced epochs will bias the center of the search space to a large extent, so that the initial search scope might not contain the true cycleslip values.
We add the same insensitive cycleslips used in lowmultipath environment into the phase data. The multipath error is depicted in Figure 3. The determination results are given in Table 10. By comparing Tables 9 and 10 we can see that the cycleslip determination based on the phase prediction shows a similar performance in low and high multipath environments. When using the phase/code combinations for cycleslip determination, the cycleslips at the , and epoch cannot be fixed in high multipath environment, whereas these cycleslips can be correctly identified in the low multipath environment as shown in Table 9.
