An Enhanced Method for Detecting and Repairing the Cycle Slips of Dual-Frequency Onboard GPS Receivers of LEO Satellites

Cycle slip detection and repair play important roles in the processing of data from dual-frequency GPS receivers onboard low-Earth orbit (LEO) satellites. To detect and repair cycle slips more comprehensively, an enhanced error method (EEM) is proposed. EEM combines single-frequency and narrow-lane carrier phase observations to construct special observations and observation equation groups. These special observations differ across time and satellite (ATS). ATS observations are constructed by three steps. The first step is differencing single-frequency and narrow-lane observations through a time difference (TD). The second step is to select a satellite as a reference satellite and other satellites as nonreference satellites. The third step is to difference the single-frequency TD observations from the reference satellite and the narrow-lane TD observations from the nonreference satellites by a satellite difference. If cycle slips occur at the reference satellite, the correction values for these ATS observations can be significantly enlarged. To process all satellites, the EEM selects each satellite as a reference satellite and builds the corresponding equation group. The EEM solves these observation equation groups according to the weighted least-squares adjustment (LSA) criterion and obtains the correction values; these correction values are then used to construct the χ2 values corresponding to different equation groups, and the EEM subsequently carries out a chi-square distribution test for these χ2. The satellite corresponding to the maximum χ2 will be marked. Then, the EEM iteratively processes the other satellites. Cycle slips can be estimated by rounding the float solutions of changes in the ambiguities of cycle slip satellites to the nearest integer. The simulation test results show that the EEM can be used to detect special cycle slip pairs such as (1, 1) and (9, 7). The EEM needs only observation data in two adjacent epochs and is still applicable to observation epochs with continuous cycle slips.


Introduction
Since low-Earth orbit (LEO) satellites are equipped with GPS receivers, they have become the main means for supporting LEO satellite missions. In the precise orbit determination of LEO satellites, the quality of the carrier phase observation data plays an important role. However, LEO satellites are in a state of high-speed movement, and both the altitude angle and the satellite attitude can change at any time; hence, the satellite signal can easily suffer from a loss of lock, which results in a discontinuous whole cycle count of carrier phase observations, a phenomenon known as cycle slip. The cycle slip of GPS carrier phase observations is an important factor that affects data quality. Carrier phase observation data with cycle slip issues will seriously affect the orbit determination accuracy. Therefore, to improve the orbit determination accuracy, it is necessary to detect and repair cycle slips. At present, there are many methods for cycle slip detection. Existing cycle slip detection methods are divided into two categories.
The first kind of cycle slip detection method, such as the polynomial fitting method, high-order difference method, Doppler integration method [1], TurboEdit method [2], and ionospheric residual method [3,4], does not involve observation equations. However, due to the high-speed movement of LEO satellites, these methods cannot be fully applied to cycle slip detection. In the polynomial fitting method, carrier phase observations are regarded as timedependent variables, and the cycle slip is determined by comparing the polynomial-fitted values with the actual observations. This method is more suitable for static or lowspeed GPS receivers and is insensitive to small and continuous cycle slips. Similar to the polynomial fitting method, the high-order difference method is not suitable for detecting small cycle slips. De Lacy et al. [5] combined Bayes' theorem with the polynomial fitting method to realize cycle slip detection, but this method is not very appropriate for LEO satellites in a high-speed moving state. The Doppler integration method requires Doppler observations, so it is not suitable for receivers with only GPS observations; furthermore, affected by the accuracy of Doppler observations, the Doppler integration method is also insensitive to small cycle slips. Blewitt [2] proposed the TurboEdit method to detect and repair cycle slips together with Melbourne-Wübbena (M-W) and geometry-free (GF) combinations. The TurboEdit method needs only data from one observation epoch, so it is widely used in GPS cycle slip detection and repair. However, due to the influence of pseudorange measurement errors, it is difficult to detect single-frequency cycle slips below 2 cycles and special cycle slip pairs such as (5,5). Joining the variation characteristics of the ionosphere with the M-W wide-lane combination, Liu [3] proposed a cycle slip detection method based on ionospheric residuals. In the detection of small cycle slips, this ionospheric residual method improves the cycle slip detection sensitivity by comparing the difference between the residual due to ionospheric delay (IOD) and the prior tolerance. The ionospheric residual method can detect a single cycle slip; thus, when the ionosphere changes slowly, this method is more suitable than other existing approaches [3]. More recently, a method was proposed that integrates the forward and backward moving window averaging (FBMWA) algorithm with the secondorder, time-difference phase ionospheric residual (STPIR) algorithm [4]. The combined FBMWA and STPIR method can detect cycle slips very well; however, when some cycle slips occur in the near field, pseudorange noise will still affect the cycle slip repair accuracy. More recently, a singlefrequency cycle slip detection and repair technique based on the generalized likelihood ratio (GLR) test was presented [6]; however, this method may also be influenced by continuous cycle slips.
The second kind of cycle slip detection involves observation equations. Bastos and Landau [7], Gao and Li [8], Colombo et al. [9], and Lee et al. [10] investigated cycle slip detection algorithms based on observation equations. However, cycle slip detection methods based on observation equations have limitations. Among the studies mentioned above, Colombo et al. [9] combined GPS and INS (internal navigation system) parameters to realize cycle slip detection, and Lee et al. [10] further studied a GPS/INS cycle slip detection approach based on the cumulative sum (CUSUM) method. However, the GPS/INS cycle slip detection method must use an INS system and thus is not suitable for cycle slip detection with LEO satellites. A method combining geometry-free (GF) and ionospheric-free (IF) double-difference observations can realize cycle slip detection better than the above mentioned techniques [11], but is suitable only for static CORS stations. To realize cycle slip detection for dynamic GPS receivers, a method combining the LAMBDA carrier phase ambiguity resolution approach with carrier phase observations was used for dynamic PPP (precise point positioning) cycle slip detection and repair [12]. More recently, a new dynamic PPP cycle slip detection method was proposed by combining the ambiguity characteristics of the single-difference observation equation with the least square method [13]. However, these dynamic PPP cycle slip detection methods are aimed at mainly ground observation stations. For cycle slip in onboard GPS receivers, a cycle slip detection method based on STP (second-order timedifference of the LEO satellite's position) and STG (secondorder time-difference geometry-free) was proposed [14]. This method does not need pseudorange observations and can reach a high detection success rate when the elevation angle of the satellite is low (less than 2.1°). However, for this method to be successful, the prior acceleration must be known, so this technique can be used as a supplement.
This research makes several contributions. First, a special time-difference and satellite-difference observation value is established. If cycle slips occur in the reference satellite, the reference satellite cycle slip error will be enhanced for all of the special observations. Second, to process the cycle slips of GPS receivers onboard LEO satellites, the enhanced error method (EEM) is proposed, which is not affected by pseudorange observations. The EEM is applicable for detecting and repairing continuous cycle slips. Third, the ability of the EEM at detecting cycle slip is verified by evaluating simulated cycle slips. We will introduce in detail the EEM cycle slip detection algorithm in the next section. Then, the method proposed will be tested in the following section and summarized in the last section.

Observational Model and Methodology
The EEM establishes a special observation value based on single-frequency (SF) and narrow-lane carrier phase observations. Narrow-lane carrier phase observations are built by dual-frequency carrier phase observations. This special observation differs across time and satellite (ATS) data. The EEM carries out independent time differences for single-frequency and narrow-lane carrier phase observations and establishes single-frequency time-difference (SFTD) and narrow-lane time-difference (NLTD) observation values. We select a satellite as the reference satellite and other satellites as nonreference satellites. Then, ATS observations are obtained by satellite differences, namely, the differences between the SFTD observations of the reference satellite and the NLTD observations of each nonreference satellite. If a cycle slip or an outlier occurs in the reference satellite, the ATS observations corresponding to the reference satellite will be corrupted, and thus, the correction values of the ATS observations will be significantly enlarged. To detect cycle slips in all satellites, the EEM selects each satellite as the reference satellite and builds the corresponding equation group. The EEM solves these observation equation groups according to the weighted least-squares adjustment (LSA) criterion and obtains the correction values of the ATS observations. The EEM uses these correction values to construct χ 2 corresponding to 2 Journal of Sensors different equation groups and carries out a chi-square distribution test for these values of χ 2 . If cycle slips occur in several satellites, the χ 2 values that do not obey the chisquare distribution exceed one. Therefore, the satellite corresponding to the maximum χ 2 will be marked provisionally as a cycle slip satellite. Then, the EEM iteratively processes the other satellites. Based on the LSA criterion, cycle slips and outliers can be estimated by the float solution of changes in the ambiguities of cycle slip satellites. In this section, NLTD and ATS observations will be introduced first, and then, the principle of the EEM, namely, the enhancement of cycle slip errors and the cycle slip detection strategy, will be explained.

Time-Difference Observation Equation.
Narrow-lane carrier phase observations can be expressed as L n = ð f 1 L 1 + f 2 L 2 Þ/ðf 1 + f 2 Þ, where L 1 and L 2 are GPS single-frequency carrier phase observations and f 1 and f 2 are the frequencies corresponding to carrier phase observations L 1 and L 2 , respectively. The time-difference equation of L n can be given as [1,15] L k where L k nði,i+1Þ represents the NLTD observation of satellite k and is constructed by the difference in L n between epochs i + 1 and i; L k nði+1Þ is the narrow-lane observation of satellite k in epoch i + 1; L k nðiÞ is the narrow-lane observation in epoch i; ρ k i+1 is the geometric distance from the receiver to satellite k in epoch i + 1; ρ k i is the geometric distance from the receiver to satellite k in epoch i; δ rði+1Þ is the receiver clock error in epoch i + 1; δ rðiÞ is the receiver clock error in epoch i; δ k sði+1Þ is the satellite clock error in epoch i + 1; δ k sðiÞ is the satellite clock error in epoch i; C is the speed of light; ΔN k nði,i+1Þ is the noise difference between the ambiguities of L k nði+1Þ and L k nðiÞ ; I k nði,i+1Þ is the variation in the IOD between epochs i + 1 and i; and ε k represents the noise in the observation data. The time-difference equation can eliminate the ambiguity parameter. If there is no cycle slip or outlier, the value of ΔN k nði,i+1Þ is 0. ρ k i+1 is calculated as follows: where ΔX ði+1,iÞ , ΔY ði+1,iÞ , and ΔZ ði+1,iÞ are the differences in the corresponding coordinates of the LEO satellite between adjacent epochs; X i , Y i , and Z i are the LEO satellite coordinates in epoch i; X k Sði+1Þ , Y k Sði+1Þ , and Z k Sði+1Þ are the GPS satellite coordinates in epoch i + 1; and ρ k i is similar to ρ k i+1 .

ATS Observation Equation.
Using L 1 and L n as examples, according to formula (1), the SFTD and NLTD observations are given as L bc 1ði,i+1Þ and L k nði,i+1Þ , respectively. The superscript bc represents the reference satellite. The ATS observation equation is given as where L ðk,bcÞ n1ði,i+1Þ represents the ATS observation; ρ k ði,i+1Þ is the difference between ρ k i+1 and ρ k i (ρ k i+1 − ρ k i ) in nonreference satellite k; ρ bc ði,i+1Þ is the difference between ρ bc i+1 and ρ bc ). Supposing that there are N satellites in an epoch, the satellite corresponding to L bc 1ði,i+1Þ is called the reference satellite, and the satellite corresponding to L k nði,i+1Þ ð0 < k ≤ N, k ≠ bcÞ is called the nonreference satellite. N − 1 ATS observations can be constructed by combining the reference and nonreference satellites. The observation equation group for N − 1 ATS observations is established as follows [1]: where A is the ðN − 1Þ × 3 coefficient matrix; l is the ðN − 1 Þ × 1 constant matrix; V is the correction vector of ATS observations to be solved; and dΔX = ½dΔX ði+1,iÞ dΔY ði+1,iÞ dΔZ ði+1,iÞ T is the parameter vector to be solved. The expression A is 3

Journal of Sensors
The expression l is Supposing P is the weight matrix of the ATS observations, then, according to the LSA criterion [16], dΔX can be expressed as The initial values of X i , Y i , and Z i and of ΔX ði+1,iÞ , Δ Y ði+1,iÞ , and ΔZ ði+1,iÞ can be obtained by single-point positioning (SPP) using ionosphere-free pseudorange observations. X k Sði+1Þ , Y k Sði+1Þ , and Z k Sði+1Þ are obtained by the precise ephemeris provided by the Center for Orbit Determination in Europe (CODE) through Lagrange interpolation. The satellite clock delays δ k sði,i+1Þ and δ bc sði,i+1Þ are obtained by the precise clock error provided by the CODE through Lagrange interpolation [17,18]. If there is no cycle slip, the theoretical value of ΔN ðk,bcÞ n1ði,i+1Þ is 0. I ðk,bcÞ n1ði,i+1Þ , which is the variation in the IOD, is very smooth in LEO satellites when the elevation angle is larger than zero [19,20]; however, when the sample interval is 1 s or higher, the variation in the IOD is distributed mainly in the range of ±0:05 cycles at low elevation angles. Hence, the variation in the IOD can be ignored (details are presented in Section 3.3).

Theory of Enhanced Cycle Slip
Error. For convenience of discussion, it is assumed that the components of the observation correction vector V are independent from each other. If no cycle slip occurs, V conforms to the multivariate normal distribution [1,21], and the expectation of V is EðVÞ = 0. If only the L 1 observations of the reference satellites have cycle slips or outliers, V will obey a multivariate skewed distribution, and the expectation of V is If cycle slips or outliers occur in nonreference satellite k, the expectation of V is According to formula (8), the EEM expands the cycle slips and outlier errors of the reference satellites to affect all the observations so that no ATS observations obey a normal distribution. In this case, the occurrence of cycle slips in the epoch can be effectively detected by using V to construct χ 2 and conducting a chi-square distribution test. According to formula (9), if the selected reference satellite is not a cycle slip satellite, the corresponding χ 2 may not obey a chi-square distribution, but its value is less than the χ 2 value corresponding to the cycle slip satellite as the reference satellite. The χ 2 statistics are constructed as follows: The original hypothesis and the alternative hypothesis are where σ 0 is the prior standard deviation (STD) of V; α is the significance level; and H 0 and H 1 are the original hypothesis and the alternative hypothesis, respectively. When χ 2 k ðn − 3Þ meets the condition of H 0 , no cycle slip is found in this epoch. Otherwise, the alternative hypothesis H 1 is accepted; i.e., cycle slip occurs in this epoch.
Formula (10) shows that the selection of the prior σ 2 0 value is the key to establishing χ 2 . We find that σ 0 is affected by the sampling interval [22][23][24]. The value of σ 0 is closely related to the data sampling interval and fluctuates greatly, so it should not be set as a fixed prior value. Therefore, an automatic method to select σ 0 is given here: supposing the current epoch is i, the mean prior STD from the initial epoch If (13) is true, epoch i + 1 can use σ i 0 as the prior STD; otherwise, epoch i + 1 still uses σ i−1 0 as the prior STD.

Novel Method for Cycle Slip Detection and Repair.
To determine cycle slips for satellites in epoch i + 1, the specific flowchart of the novel cycle slip detection method is given in Figure 1. The corresponding procedures are as follows:  [25,26]. Supposing that cycle slip occurs only in satellite k , the new equation groups can be expressed as where ΔX ΔN can be solved as

Journal of Sensors
The cofactor matrix of unknown parameters can be given as Q XX is used to divide cycle slips and outliers.
While the float solution of a cycle slip is estimated, fixing the cycle slip as an integer is important. The LAMBDA method is generally adopted to round cycle slips to the nearest integer [12,27]. However, as shown in (17), the accuracy of ΔX ΔN is not affected by pseudorange observations. Therefore, the cycle slip can be accurately rounded to the nearest integer. We use the STD values of the float solutions to dis-tinguish outliers and cycle slips. The detection formula is written as where Q ΔNΔN is the cofactor of ΔX ΔN in Q XX . On the one hand, while (21) is true, round (ΔX ΔN ) can be viewed as a cycle slip. On the other hand, if (21) is false, the satellite corresponding to ΔX ΔN will be marked as an outlier satellite.
(5) There is a correlation between the ATS observation values, and the corresponding prior weight matrix P is no longer a diagonal matrix. Therefore, P cannot be simply defined as a unit matrix [28]. The strategy for processing P is given below:      (12), (13), and (14) to renew the prior STD It should be mentioned that the validity of the EEM depends on the quality of observations in the present epoch. If the number of satellites without cycle slips is less than 5, the EEM is not valid. Therefore, to detect and repair cycle slips more robustly, the EEM can be combined with other methods to address cycle slips.

Simulated Case Study
The case data are taken from the onboard GPS receiver of the Swarm-A satellite provided by the European Space Agency (ESA). The precise GPS ephemeris with a 15 min sampling interval and the precise satellite clock error file with a 30 s sampling interval are provided by the CODE. The date on which the Swarm-A satellite data were acquired is DOY 145 of 2018 (GPS time 2018:05:25), and the sampling interval is 1 s. Table 1 lists the detailed web addresses of the data sources.  Table 2 shows the distribution of the simulated cycle slips and the detection and repair results. The detection results indicate that these special cycle slip pairs can be detected by the EEM and by the FBMWA and TurboEdit methods. On the one hand, for these cycle slip pairs inserted every 30 epochs, the noise in the pseudorange observations is effectively smoothed by the FBMWA algorithm. However, several repair errors occur in epochs 410, 520, and 670. The size of the moving window is 25 seconds [4]. It may be that the short size of the moving window is not enough to reduce the noise level in the pseudorange observations. On the other hand, the EEM is sensitive to special cycle slip pairs such as (1, 1), (3,4), and (9, 7). According to equation (17), these simulated cycle slip pairs can be repaired well. These estimated cycle slip values are quite close to the simulated cycle slips. This is because equation (17) does not use pseudorange observations to repair cycle slips. As a result, these cycle slips can be accurately rounded to the nearest integer. Therefore, combining the EEM with the FBMWA and TurboEdit algorithms is more robust at processing cycle slips than using any of these methods alone.
To verify the ability of the EEM to detect and repair continuous cycle slips from epoch 300 (GPS time 2018:05:25:00:05:00) to epoch 304, continuous cycle slips are inserted into G01. Table 3 lists the detection and repair results. The results show that the EEM algorithm can detect and repair cycle slips that occur in continuous epochs.

Variation in the Ionospheric
Delay. The ionospheric delay (IOD) commonly varies quickly in LEO satellites, as shown in previous studies [4,14]. However, rapid variation in the IOD is not always observed when the elevation angle is larger than zero [19,20]. Generally, onboard GPS receivers can provide a high sampling interval, such as 1 s [29][30][31]. Therefore, a significant question is raised as to whether the variation in the IOD can be neglected in the processing of cycle slips for LEO satellites with a high sampling interval (1 s) and a low elevation angle. To solve this question, we use a geometry-free linear combination to assess the variation in the IOD [20] and use equation (17) to estimate ΔN 1 and ΔN 2 . Figures 2(a)-2(d) show the variations in the IOD (sampling interval: 1 s). Figures 3(a)-3(h) show the variations in ΔN 1 and ΔN 2 . As these figures show, we can reach several conclusions. First, the ambiguity and IOD vary slowly when the elevation angle is larger than 10°. These variations are distributed mainly in the range of ±0:05 cycles. Second, the maximum variation is lower than 0.1 cycles when the elevation angle approaches zero. This finding confirms that the IOD variation is stable. Finally, the variations in ΔN 1 and Δ N 2 are significantly small because the estimates of ΔN 1 and ΔN 2 are not affected by pseudorange noise, IOD variations, or other errors, even at low elevation angles. These phenomena suggest that (1) the variations in the ambiguity and IOD

13
Journal of Sensors obey a normal distribution and that (2) the variation in the IOD does not obviously corrupt time-difference carrier phase observations. Hence, the variation in the IOD of timedifference carrier phase observations in equations (3) and (17) can be ignored for the proposed cycle slip detection and repair method.
To investigate the performance of the proposed method at low elevation angles, we select G01 and G17 as examples All simulated cycle slip pairs (9,7) and (1,1) in G01 and G17 were detected and repaired. As shown in Figures 5(e) and 5(f), in the simulated cycle slip epochs, the values of ΔN 1 and ΔN 2 are close to the simulated cycle slip pair (1,1). This result verifies that these ΔN can be easily rounded to their nearest integers.  14 Journal of Sensors interval, the mean difference in the posterior STD of the ATS observations constructed by different reference satellites is small. The mean difference in the posterior STD between the sampling intervals of 1 s and 60 s reaches an order of magnitude. The accuracy of V is better than 5 cm when the sampling interval of the observation data is higher than or equal to 15 s, whereas the accuracy of V is less than 5 cm when the sampling interval of the observation data is less than or equal to 30 s. This may be because with a decrease in the sampling interval, the variation in the IOD between epochs is expanded and cannot be greatly weakened by the difference between epochs at a low sampling interval [4].

Standard Deviation of the Observation Correction Vector.
The above data analysis shows that the sensitivity of the EEM algorithm increases with the increase in the sampling interval of the observation data, so the algorithm is more suitable for data acquired with a high sampling frequency. At the same time, this analysis confirms that formula (12) should be used to estimate the prior STD.     15 Journal of Sensors inserted simulated cycle slips and outliers. Table 6 lists the χ 2 values. The results show that when G11 is selected as the reference satellite, its χ 2 value is much larger than the values of the other observation satellites selected as the reference satellite. To obtain experimental results that are more statistically representative, an outlier of 0.5 cycles is inserted into satellite G11 every 10 s during the continuous observation epoch from GPS time 2018:05:25:00:00:20 to 2018:05:25:00:16:00 with a sampling interval of 1 s. Figure 6 shows the χ 2 values for the satellites with inserted and no inserted simulated values. The results show that when satellite G11 is selected as the reference satellite with an inserted outlier of 0.5 cycles, the corresponding χ 2 value is greater than the values of the other reference satellites without inserted simulated outliers. These results verify that the EEM algorithm can effectively enlarge the errors caused by small cycle slips and outliers.

Conclusion
A method for detecting and repairing onboard GPS receiver cycle slips is presented. This method combines singlefrequency and narrow-lane observations to construct ATS observations, which can significantly enlarge the small errors caused by cycle slips and outliers. By amplifying the influences caused by cycle slip and outliers, the EEM algorithm can detect cycle slips more stably than other existing algorithms.
The sensitivity of the EEM algorithm is evaluated by inserting simulated cycle slips into different reference satellites. The results show that the EEM, FBMWA, and TurboEdit methods have similar repair accuracies for cycle slips. However, compared with the FBMWA and TurboEdit methods, the EEM is not affected by pseudorange observations, which makes up for the shortcomings of the FBMWA and TurboEdit algorithms. The proposed algorithm is applicable even to observation satellite epochs with continuous cycle slips and outliers.
The proposed EEM cycle slip detection method can effectively detect and repair cycle slips and outliers and is still applicable to observation epochs with continuous cycle slips. However, for low-frequency observation data (sampling interval less than 30 s), the EEM algorithm cannot eliminate the impact of variation in the IOD. In addition, the EEM algorithm depends on the quality of the observation data from other satellites at the same time. Therefore, to thoroughly detect onboard GPS receiver cycle slips, the EEM, FBMWA, and TurboEdit algorithms should be combined to process cycle slips simultaneously.

Data Availability
The space-borne GPS observation data used in our work can be freely accessed at https://swarm-diss.eo.esa.int/. The GPS satellite ephemeris and clock products are obtained from ftp://cddis.nasa.gov/gnss/products.

Conflicts of Interest
The authors declare that they have no conflicts of interest.