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

Authors: Zhouming Yang; Xin Liu; Jinyun Guo; Yaowei Xia; Xiaotao Chang
Journal: Journal of Sensors
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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 time-dependent 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 low-speed 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 second-order, 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 single-frequency 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 time-difference of the LEO satellite’s position) and STG (second-order 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.


## 2. 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 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 chi-square 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.


### 2.1. Time-Difference Observation Equation

Narrow-lane carrier phase observations can be expressed asLn=f1L1+f2L2/f1+f2, where L1 and L2 are GPS single-frequency carrier phase observations and f1 and f2 are the frequencies corresponding to carrier phase observations L1 and L2, respectively. The time-difference equation of Ln can be given as [1, 15]
(1)Lni,i+1k=Lni+1k−Lnik=ρi+1k−ρik+Cδri+1−Cδri−Cδsi+1k+Cδsik+ΔNni,i+1k+Ini,i+1k+εk,where Lni,i+1k represents the NLTD observation of satellite k and is constructed by the difference in Ln between epochs i+1 and i; Lni+1k is the narrow-lane observation of satellite k in epoch i+1; Lnik is the narrow-lane observation in epoch i; ρi+1k is the geometric distance from the receiver to satellite k in epoch i+1; ρik is the geometric distance from the receiver to satellite k in epoch i; δri+1 is the receiver clock error in epoch i+1; δri is the receiver clock error in epoch i; δsi+1k is the satellite clock error in epoch i+1; δsik is the satellite clock error in epoch i; C is the speed of light; ΔNni,i+1k is the noise difference between the ambiguities of Lni+1k and Lnik; Ini,i+1k 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 ΔNni,i+1k is 0. ρi+1k is calculated as follows:
(2)ρi+1k=Xi+ΔXi+1,i−XSi+1k2+Yi+ΔYi+1,i−YSi+1k2+Zi+ΔZi+1,i−ZSi+1k2,where ΔXi+1,i, ΔYi+1,i, and ΔZi+1,i are the differences in the corresponding coordinates of the LEO satellite between adjacent epochs; Xi, Yi, and Zi are the LEO satellite coordinates in epoch i; XSi+1k, YSi+1k, and ZSi+1k are the GPS satellite coordinates in epoch i+1; and ρik is similar to ρi+1k.


### 2.2. ATS Observation Equation

UsingL1 and Ln as examples, according to formula (1), the SFTD and NLTD observations are given as L1i,i+1bc and Lni,i+1k, respectively. The superscript bc represents the reference satellite. The ATS observation equation is given as
(3)Ln1i,i+1k,bc=Lni,i+1k−L1i,i+1bc=ρi,i+1k−ρi,i+1bc−Cδsi,i+1k+Cδsi,i+1bc+ΔNn1i,i+1k,bc+In1i,i+1k,bc+εk,bc,where Ln1i,i+1k,bc represents the ATS observation; ρi,i+1k is the difference between ρi+1k and ρik(ρi+1k−ρik) in nonreference satellite k; ρi,i+1bc is the difference between ρi+1bc and ρibc (ρi+1bc−ρibc) in reference satellite bc; δsi,i+1k is the difference between δsi+1k and δsik (δsi+1k−δsik); δsi,i+1bc is the difference between δsi+1bc and δsibc (δsi+1bc−δsibc); ΔNn1i,i+1k,bc is the difference between ΔNni,i+1k and ΔN1i,i+1bc (ΔNni,i+1k−ΔN1i,i+1bc); and In1i,i+1k,bc is the difference between Ini,i+1k and I1i,i+1bc (Ini,i+1k−I1i,i+1bc).Supposing that there areN satellites in an epoch, the satellite corresponding to L1i,i+1bc is called the reference satellite, and the satellite corresponding to Lni,i+1k0<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]:
(4)V=AdΔX−l,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ΔXi+1,idΔYi+1,idΔZi+1,iT is the parameter vector to be solved. The expression A is
(5)A=Xi+ΔXi+1,i−XSi+11ρi+11Yi+ΔYi+1,i−YSi+11ρi+11Zi+ΔZi+1,i−ZSi+11ρi+11⋮⋮⋮Xi+ΔXi+1,i−XSi+1N−1ρi+1N−1Yi+ΔYi+1,i−YSi+1N−1ρi+1N−1Zi+ΔZi+1,i−ZSi+1N−1ρi+1N−1.The expressionl is
(6)l=ρi,i+11−ρi,i+1bc−Cδsi,i+11+Cδsi,i+1bc−Ln1i,i+11,bc⋮ρi,i+1N−1−ρi,i+1bc−Cδsi,i+1N−1+Cδsi,i+1bc−Ln1i,i+1N−1,bc.SupposingP is the weight matrix of the ATS observations, then, according to the LSA criterion [16], dΔX can be expressed as
(7)dΔX=ATPA−1ATPl.The initial values ofXi, Yi, and Zi and of ΔXi+1,i, ΔYi+1,i, and ΔZi+1,i can be obtained by single-point positioning (SPP) using ionosphere-free pseudorange observations. XSi+1k, YSi+1k, and ZSi+1k are obtained by the precise ephemeris provided by the Center for Orbit Determination in Europe (CODE) through Lagrange interpolation. The satellite clock delays δsi,i+1k and δsi,i+1bc 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 ΔNn1i,i+1k,bc is 0. In1i,i+1k,bc, 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).


### 2.3. Theory of Enhanced Cycle Slip Error

For convenience of discussion, it is assumed that the components of the observation correction vectorV 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 EV=0. If only the L1 observations of the reference satellites have cycle slips or outliers, V will obey a multivariate skewed distribution, and the expectation of V is
(8)EV=ΔNn1i,i+11,bc⋯ΔNn1i,i+1N−1,bcT.If cycle slips or outliers occur in nonreference satellitek, the expectation of V is
(9)EV=0⋯0Nn1i,i+1k,bc0⋯0T.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:
(10)χk2n−3=VTPVσ02.The original hypothesis and the alternative hypothesis are(11)H0:χk2n−3≤χα2,H1:χk2n−3>χα2,where σ0 is the prior standard deviation (STD) of V; α is the significance level; and H0 and H1 are the original hypothesis and the alternative hypothesis, respectively. When χk2n−3 meets the condition of H0, no cycle slip is found in this epoch. Otherwise, the alternative hypothesis H1 is accepted; i.e., cycle slip occurs in this epoch.Formula (10) shows that the selection of the prior σ02 value is the key to establishing χ2. We find that σ0 is affected by the sampling interval [22–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 to epoch i−1 is σ0i−1¯, and the variance of σ0i−1¯ is σσ0i−1¯2. The formula for choosing reasonable values of σ0i¯ and σσ0i¯2 for epoch i+1 is as follows:
(12)σ0i¯=σ0i−1¯+1iσ0i−σ0i−1¯,(13)σσ0i¯2=σσ0i−1¯2+1iσ0i−σ0i−1¯2−σσ0i−1¯2.Fromσ0i−1¯ and σσ0i−1¯2 of epoch i−1, combined with the PauTa criterion, the σ0i¯ of epoch i can be detected:
(14)σ0i¯−σ0i−1¯<3σσ0i−1¯2.If (13) is true, epoch i+1 can use σ0i¯ as the prior STD; otherwise, epoch i+1 still uses σ0i−1¯ as the prior STD.


### 2.4. Novel Method for Cycle Slip Detection and Repair

To determine cycle slips for satellites in epochi+1, the specific flowchart of the novel cycle slip detection method is given in Figure 1. The corresponding procedures are as follows:
(1)
Supposing that there areN satellites in adjacent epochs, equation (4) is used to construct N observation equation groups. Each group has N−1 observation equations and corresponds to a reference satellite sk0<k≤N(2)
According to equation (7), N observation equation groups are solved to obtain N observation correction vectors. According to equation (10), these N vectors are used to construct N chi-square values χk2k=1,⋯,N. The subscript k represents the chi-square value corresponding to reference satellite sk(a)
If only a singleχk2 value does not obey a chi-square distribution, provisionally mark sk as a cycle slip satellite. Move to step (4)(b)
If more than oneχk2 value fails to obey a chi-square distribution, provisionally mark the satellite that corresponds to the maximum χk2 value as the cycle slip satellite and eliminate the satellite. Set N=N−1. Move to step (1)(c)
IfN=4, move to step (3)(3)
If all cycle slip satellites have not been exactly marked, the normal satellites will be eliminated until 4 satellites remain. Considering this extreme condition, a more robust detection strategy is given(a)
Initializek=1; then, satellite Sk0<k≤N can be temporarily eliminated. According to equations (4) and (10), other satellites can be used to establish new values of χk2. Renew kk=2,⋯,N, and then, obtain the corresponding χk2(b)
Select the satellite corresponding to the minimumχ2 as a cycle slip satellite. SetN=N−1. Move to step (a) until all cycle slips have been marked(4)
These cycle slips are marked as unknown parameters. These parameters can be calculated according to the float solution of changes in the ambiguities of cycle slip satellites [25, 26]. Supposing that cycle slip occurs only in satellite k, the new equation groups can be expressed as(15)V=AdΔX+CkΔXΔN−l,where Ck=0⋯1⋯0T is the N−1×1 coefficient matrix and ΔXΔN is an unknown parameter of ΔNi,i+1k,bc. ΔNi,i+1k,bc can be directly used to calculate ΔN1 or ΔN2. All of the unknown parameters can be given as
(16)dXΔXΔN=ATPAATPCCTPACTPC−1ATPlCTPl.Figure 1
Cycle slip detection algorithm.ΔXΔN can be solved as
(17)ΔXΔN=QΔNΔNu,(18)QΔNΔN−1=CTPC−CTPAATPA−1ATPC,(19)u=CTPl−CTPAATPA−1ATPl.The cofactor matrix of unknown parameters can be given as(20)QXX=ATPAATPCCTPACTPC−1.QXX 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 distinguish outliers and cycle slips. The detection formula is written as
(21)ΔXΔN−roundΔXΔN<3σ0QΔNΔN,where QΔNΔN is the cofactor of ΔXΔN in QXX. 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 matrixP 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:
(a)
P is a unit matrix in the initial epoch P=I(b)
The prior cofactor matrix can be given asQ=P−1(c)
According to the formulaQvv=Q−AATPA−1AT, calculate the cofactor matrix of V(d)
According to the formulaP=Qvv−1, calculate the posterior weight matrix and use it as the prior weight matrix of the next epoch(e)
Use equations (12), (13), and (14) to renew the prior STDIt 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.


## 2.1. Time-Difference Observation Equation

Narrow-lane carrier phase observations can be expressed asLn=f1L1+f2L2/f1+f2, where L1 and L2 are GPS single-frequency carrier phase observations and f1 and f2 are the frequencies corresponding to carrier phase observations L1 and L2, respectively. The time-difference equation of Ln can be given as [1, 15]
(1)Lni,i+1k=Lni+1k−Lnik=ρi+1k−ρik+Cδri+1−Cδri−Cδsi+1k+Cδsik+ΔNni,i+1k+Ini,i+1k+εk,where Lni,i+1k represents the NLTD observation of satellite k and is constructed by the difference in Ln between epochs i+1 and i; Lni+1k is the narrow-lane observation of satellite k in epoch i+1; Lnik is the narrow-lane observation in epoch i; ρi+1k is the geometric distance from the receiver to satellite k in epoch i+1; ρik is the geometric distance from the receiver to satellite k in epoch i; δri+1 is the receiver clock error in epoch i+1; δri is the receiver clock error in epoch i; δsi+1k is the satellite clock error in epoch i+1; δsik is the satellite clock error in epoch i; C is the speed of light; ΔNni,i+1k is the noise difference between the ambiguities of Lni+1k and Lnik; Ini,i+1k 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 ΔNni,i+1k is 0. ρi+1k is calculated as follows:
(2)ρi+1k=Xi+ΔXi+1,i−XSi+1k2+Yi+ΔYi+1,i−YSi+1k2+Zi+ΔZi+1,i−ZSi+1k2,where ΔXi+1,i, ΔYi+1,i, and ΔZi+1,i are the differences in the corresponding coordinates of the LEO satellite between adjacent epochs; Xi, Yi, and Zi are the LEO satellite coordinates in epoch i; XSi+1k, YSi+1k, and ZSi+1k are the GPS satellite coordinates in epoch i+1; and ρik is similar to ρi+1k.


## 2.2. ATS Observation Equation

UsingL1 and Ln as examples, according to formula (1), the SFTD and NLTD observations are given as L1i,i+1bc and Lni,i+1k, respectively. The superscript bc represents the reference satellite. The ATS observation equation is given as
(3)Ln1i,i+1k,bc=Lni,i+1k−L1i,i+1bc=ρi,i+1k−ρi,i+1bc−Cδsi,i+1k+Cδsi,i+1bc+ΔNn1i,i+1k,bc+In1i,i+1k,bc+εk,bc,where Ln1i,i+1k,bc represents the ATS observation; ρi,i+1k is the difference between ρi+1k and ρik(ρi+1k−ρik) in nonreference satellite k; ρi,i+1bc is the difference between ρi+1bc and ρibc (ρi+1bc−ρibc) in reference satellite bc; δsi,i+1k is the difference between δsi+1k and δsik (δsi+1k−δsik); δsi,i+1bc is the difference between δsi+1bc and δsibc (δsi+1bc−δsibc); ΔNn1i,i+1k,bc is the difference between ΔNni,i+1k and ΔN1i,i+1bc (ΔNni,i+1k−ΔN1i,i+1bc); and In1i,i+1k,bc is the difference between Ini,i+1k and I1i,i+1bc (Ini,i+1k−I1i,i+1bc).Supposing that there areN satellites in an epoch, the satellite corresponding to L1i,i+1bc is called the reference satellite, and the satellite corresponding to Lni,i+1k0<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]:
(4)V=AdΔX−l,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ΔXi+1,idΔYi+1,idΔZi+1,iT is the parameter vector to be solved. The expression A is
(5)A=Xi+ΔXi+1,i−XSi+11ρi+11Yi+ΔYi+1,i−YSi+11ρi+11Zi+ΔZi+1,i−ZSi+11ρi+11⋮⋮⋮Xi+ΔXi+1,i−XSi+1N−1ρi+1N−1Yi+ΔYi+1,i−YSi+1N−1ρi+1N−1Zi+ΔZi+1,i−ZSi+1N−1ρi+1N−1.The expressionl is
(6)l=ρi,i+11−ρi,i+1bc−Cδsi,i+11+Cδsi,i+1bc−Ln1i,i+11,bc⋮ρi,i+1N−1−ρi,i+1bc−Cδsi,i+1N−1+Cδsi,i+1bc−Ln1i,i+1N−1,bc.SupposingP is the weight matrix of the ATS observations, then, according to the LSA criterion [16], dΔX can be expressed as
(7)dΔX=ATPA−1ATPl.The initial values ofXi, Yi, and Zi and of ΔXi+1,i, ΔYi+1,i, and ΔZi+1,i can be obtained by single-point positioning (SPP) using ionosphere-free pseudorange observations. XSi+1k, YSi+1k, and ZSi+1k are obtained by the precise ephemeris provided by the Center for Orbit Determination in Europe (CODE) through Lagrange interpolation. The satellite clock delays δsi,i+1k and δsi,i+1bc 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 ΔNn1i,i+1k,bc is 0. In1i,i+1k,bc, 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).


## 2.3. Theory of Enhanced Cycle Slip Error

For convenience of discussion, it is assumed that the components of the observation correction vectorV 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 EV=0. If only the L1 observations of the reference satellites have cycle slips or outliers, V will obey a multivariate skewed distribution, and the expectation of V is
(8)EV=ΔNn1i,i+11,bc⋯ΔNn1i,i+1N−1,bcT.If cycle slips or outliers occur in nonreference satellitek, the expectation of V is
(9)EV=0⋯0Nn1i,i+1k,bc0⋯0T.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:
(10)χk2n−3=VTPVσ02.The original hypothesis and the alternative hypothesis are(11)H0:χk2n−3≤χα2,H1:χk2n−3>χα2,where σ0 is the prior standard deviation (STD) of V; α is the significance level; and H0 and H1 are the original hypothesis and the alternative hypothesis, respectively. When χk2n−3 meets the condition of H0, no cycle slip is found in this epoch. Otherwise, the alternative hypothesis H1 is accepted; i.e., cycle slip occurs in this epoch.Formula (10) shows that the selection of the prior σ02 value is the key to establishing χ2. We find that σ0 is affected by the sampling interval [22–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 to epoch i−1 is σ0i−1¯, and the variance of σ0i−1¯ is σσ0i−1¯2. The formula for choosing reasonable values of σ0i¯ and σσ0i¯2 for epoch i+1 is as follows:
(12)σ0i¯=σ0i−1¯+1iσ0i−σ0i−1¯,(13)σσ0i¯2=σσ0i−1¯2+1iσ0i−σ0i−1¯2−σσ0i−1¯2.Fromσ0i−1¯ and σσ0i−1¯2 of epoch i−1, combined with the PauTa criterion, the σ0i¯ of epoch i can be detected:
(14)σ0i¯−σ0i−1¯<3σσ0i−1¯2.If (13) is true, epoch i+1 can use σ0i¯ as the prior STD; otherwise, epoch i+1 still uses σ0i−1¯ as the prior STD.


## 2.4. Novel Method for Cycle Slip Detection and Repair

To determine cycle slips for satellites in epochi+1, the specific flowchart of the novel cycle slip detection method is given in Figure 1. The corresponding procedures are as follows:
(1)
Supposing that there areN satellites in adjacent epochs, equation (4) is used to construct N observation equation groups. Each group has N−1 observation equations and corresponds to a reference satellite sk0<k≤N(2)
According to equation (7), N observation equation groups are solved to obtain N observation correction vectors. According to equation (10), these N vectors are used to construct N chi-square values χk2k=1,⋯,N. The subscript k represents the chi-square value corresponding to reference satellite sk(a)
If only a singleχk2 value does not obey a chi-square distribution, provisionally mark sk as a cycle slip satellite. Move to step (4)(b)
If more than oneχk2 value fails to obey a chi-square distribution, provisionally mark the satellite that corresponds to the maximum χk2 value as the cycle slip satellite and eliminate the satellite. Set N=N−1. Move to step (1)(c)
IfN=4, move to step (3)(3)
If all cycle slip satellites have not been exactly marked, the normal satellites will be eliminated until 4 satellites remain. Considering this extreme condition, a more robust detection strategy is given(a)
Initializek=1; then, satellite Sk0<k≤N can be temporarily eliminated. According to equations (4) and (10), other satellites can be used to establish new values of χk2. Renew kk=2,⋯,N, and then, obtain the corresponding χk2(b)
Select the satellite corresponding to the minimumχ2 as a cycle slip satellite. SetN=N−1. Move to step (a) until all cycle slips have been marked(4)
These cycle slips are marked as unknown parameters. These parameters can be calculated according to the float solution of changes in the ambiguities of cycle slip satellites [25, 26]. Supposing that cycle slip occurs only in satellite k, the new equation groups can be expressed as(15)V=AdΔX+CkΔXΔN−l,where Ck=0⋯1⋯0T is the N−1×1 coefficient matrix and ΔXΔN is an unknown parameter of ΔNi,i+1k,bc. ΔNi,i+1k,bc can be directly used to calculate ΔN1 or ΔN2. All of the unknown parameters can be given as
(16)dXΔXΔN=ATPAATPCCTPACTPC−1ATPlCTPl.Figure 1
Cycle slip detection algorithm.ΔXΔN can be solved as
(17)ΔXΔN=QΔNΔNu,(18)QΔNΔN−1=CTPC−CTPAATPA−1ATPC,(19)u=CTPl−CTPAATPA−1ATPl.The cofactor matrix of unknown parameters can be given as(20)QXX=ATPAATPCCTPACTPC−1.QXX 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 distinguish outliers and cycle slips. The detection formula is written as
(21)ΔXΔN−roundΔXΔN<3σ0QΔNΔN,where QΔNΔN is the cofactor of ΔXΔN in QXX. 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 matrixP 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:
(a)
P is a unit matrix in the initial epoch P=I(b)
The prior cofactor matrix can be given asQ=P−1(c)
According to the formulaQvv=Q−AATPA−1AT, calculate the cofactor matrix of V(d)
According to the formulaP=Qvv−1, calculate the posterior weight matrix and use it as the prior weight matrix of the next epoch(e)
Use equations (12), (13), and (14) to renew the prior STDIt 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.


## 3. 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. Table1 lists the detailed web addresses of the data sources.Table 1
Data source details.
Type of dataReleasing agencySite of downloadData descriptionSwarm-A satellite-borne GPS observationsESAhttps://swarm-diss.eo.esa.intRINEX 3.0, sampling interval (1 s)GPS satellite precise ephemerisCODEhttp://ftp.aiub.unibe.ch/CODE/2018/Sampling interval (15 min)GPS satellite precision clock errorCODEhttp://ftp.aiub.unibe.ch/CODE/2018/Sampling interval (30 s)The certain times at which the observation data of a satellite are acquired are regarded as observation satellite epochs. A complete observation satellite epoch includes carrier phase observations (L1 and L2) and pseudorange observations (P1 and P2). Incomplete observation satellite epochs are eliminated. On DOY 145 of 2018, Swarm-A had 684087 observation satellite epochs. Complete observation satellite epochs account for 97.6% of the total. This means that the observation results of Swarm-A on DOY 145 are relatively complete and can be used as test data.


### 3.1. Special Cycle Slip Pair Detection and Repair

To test the validity of the new method, from epoch 60 (GPS time 2018:05:25:00:01:00) to epoch 670 (GPS time 2018:05:25:00:11:10), some special simulated cycle slip pairs are inserted into satellite G01. Table2 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.Table 2
Results of cycle slip detection and repair by two methods for G01.
EpochSlip (L1, L2)FBMWA and TurboEditEEMCycle slip detectableCycle slip repairableTurboEditEEMTurboEditEEM60(9, 7)(8.991, 6.991)(8.986, 6.991)YesYesYesYes120(6, 6)(5.999, 5.999)(6.006, 6.014)YesYesYesYes180(5, 5)(5.048, 5.048)(5.004, 5.014)YesYesYesYes240(3, 4)(3.027, 4.027)(3.018, 4.030)YesYesYesYes300(4, 3)(3.988, 2.988)(3.997, 2.996)YesYesYesYes360(1, 0)(1.010, 0.010)(0.998, 0.007)YesYesYesYes410(1, 1)(-3.548, 2.548)(1.010, 1.021)YesYesNoYes420(1, 1)(0.988, 0.988)(1.006, 0.987)YesYesYesYes480(77, 60)(77.042, 60.042)(76.986, 59.984)YesYesYesYes520(4, 3)(8.580, 6.580)(4.019, 3.033)YesYesNoYes660(3, 4)(2.971, 3.971)(3.014, 4.006)YesYesYesYes670(5, 5)(9.454, 8.454)(5.035, 5.015)YesYesNoYesTo 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. Table3 lists the detection and repair results. The results show that the EEM algorithm can detect and repair cycle slips that occur in continuous epochs.Table 3
Results of detecting and repairing continuous cycle slips in G01.
EpochSimulated cycle slip pairs (L1, L2)New method300(1, 1)(1.103, 1.096)301(1, 1)(1.096, 1.019)302(1, 1)(0.901, 1.000)303(9, 7)(8.996, 6.892)304(9, 7)(8.808, 6.993)


### 3.2. 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–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 ΔN1 and ΔN2. Figures 2(a)–2(d) show the variations in the IOD (sampling interval: 1 s). Figures 3(a)–3(h) show the variations in ΔN1 and ΔN2. 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 ΔN1 and ΔN2 are significantly small because the estimates of ΔN1 and ΔN2 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 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 time-difference carrier phase observations in equations (3) and (17) can be ignored for the proposed cycle slip detection and repair method.Figure 2
Elevation angles and variations in the IOD for (a) G01, (b) G08, (c) G15, and (d) G17 on DOY 145 of 2018.
(a)(b)(c)(d)Figure 3
Elevation angles,ΔN1, and ΔN2 for (a, b) G01, (c, d) G08, (e, f) G15, and (g, h) G17 on DOY 145 of 2018.
(a)(b)(c)(d)(e)(f)(g)(h)To investigate the performance of the proposed method at low elevation angles, we select G01 and G17 as examples (sampling interval: 1 s). From epoch 360 (GPS time 2018:05:25:00:01:00) to epoch 930 (GPS time 2018:05:25:00:15:30), 20 pairs of the simulated cycle slip pair (9, 7) are added to G01. The elevation angles in the period from epoch 600 to epoch 930 are all lower than 15°. From epoch 360 to epoch 570 (GPS time 2018:05:25:00:9:30), 8 pairs of the simulated cycle slip pair (1, 1) are added to G17. The elevation angles in the period from epoch 360 to epoch 570 are all lower than 15°. These cycle slip pairs are described in Figures4(a) and 4(b).Figure 4
Distribution of simulated cycle slip pairs at an interval of every 30 epochs.
(a)(b)The cycle slip repair results for G01 and G17 are shown in Figures5(a)–5(f). In each plot, the values of ΔN and the variation in the elevation angle are given. To visualize the results better, values of ΔN larger than 1.2 cycles were set as 1.2 cycles. These cycle slip epochs are indicated by vertical lines. 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 ΔN1 and ΔN2 are close to the simulated cycle slip pair (1, 1). This result verifies that these ΔN can be easily rounded to their nearest integers.Figure 5
Elevation angle,ΔN1, and ΔN2 for G01 and G17.


### 3.3. Standard Deviation of the Observation Correction Vector

Constructingχ2 requires the prior STD of the observation correction vector (V). To study the distribution of the STD, the STD values corresponding to each epoch and reference satellite are calculated for the Swarm-A observation data on DOY 145 according to sampling intervals of 1 s, 5 s, 10 s, 15 s, 30 s, and 60 s. According to the formula σ0=VTPV/n−3, the mean posterior STD values corresponding to different reference satellites in each epoch are calculated [1]. The statistical results related to the posterior STD are listed in Tables 4 and 5, which show that the posterior STD of the observation correction vector under different reference satellites is closely related to the data sampling interval, and the value of the posterior STD increases with a decrease in the sampling interval. At the same sampling 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]. 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.Table 4
STD values ofV under different basic satellites (unit: m).
GPS PRNSampling interval1 s5 s10 s15 s30 s60 sG010.00290.01070.02010.02900.05090.0794G020.00270.01070.02020.02960.05340.0819G030.00270.00970.01820.02640.04780.0788G050.00280.01080.02060.03010.05560.0858G060.00260.01010.01910.02770.04930.0801G070.00280.01060.02000.02920.05300.0796G080.00290.01100.02060.02930.05250.0792G090.00280.01050.01980.02870.04970.0788G100.00320.01120.01960.02800.05170.0839G110.00290.01110.02080.03020.05320.0760G120.00270.01050.01940.02810.05090.0807G130.00280.01080.02030.02950.05310.0794G140.00260.00940.01750.02510.04470.0765G150.00320.01180.02070.02970.05350.0815G160.00270.00980.01830.02660.04940.0826G170.00280.01090.02050.02950.05320.0860G180.00300.01080.02020.02920.05150.0854G190.00280.01100.02080.03030.05460.0880G200.00360.01230.02140.03020.05420.0862G210.00310.01110.02040.02960.05410.0869G220.00290.01030.01920.02790.04780.0809G230.00270.01010.01910.02780.04980.0803G240.00310.01150.02010.02850.05130.0817G250.00270.01010.01870.02710.04880.0795G260.00270.00990.01880.02730.04840.0816G270.00300.01050.01870.02670.04760.0760G280.00300.01150.02150.03090.05720.0863G290.00300.01080.01950.02810.05010.0793G300.00280.01020.01920.02790.05040.0770G310.00260.00930.01740.02520.04550.0750G320.00300.01080.01920.02690.04880.0808Table 5
Statistic results of the STD (unit: m).
Sampling intervalσ0MaxMinMean1 s0.00360.00260.00295 s0.01230.00930.010710 s0.02150.01740.019715 s0.03100.02520.028530 s0.05730.04480.051160 s0.08800.07500.0812


### 3.4. Performance of the EEM to Enhance Cycle Slip Errors

To verify the ability of the EEM to amplify cycle slip errors, a simulated outlier of 0.5 cycles and a cycle slip of 1 cycle are inserted into theL1 carrier phase of G11 in the observation epoch corresponding to GPS time 2018:05:25:00:06:30. Then, we calculate the χ2 values corresponding to the different 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.Table 6
Values ofχ2 for different basic satellites.
GPS PRNValue ofχ2No cycle slipOutlier of 0.5 cycles1 cycle slipG010.2740.23464.206G080.2739.71962.558G110.27278.374443.99G130.38133.1252.491G150.73945.65762.921G170.27332.95452.071G280.433.42553.079G301.14958.97779.322Figure 6
Values ofχ2 for different reference satellites under an outlier of 0.5 cycles.


## 3.1. Special Cycle Slip Pair Detection and Repair

To test the validity of the new method, from epoch 60 (GPS time 2018:05:25:00:01:00) to epoch 670 (GPS time 2018:05:25:00:11:10), some special simulated cycle slip pairs are inserted into satellite G01. Table2 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.Table 2
Results of cycle slip detection and repair by two methods for G01.
EpochSlip (L1, L2)FBMWA and TurboEditEEMCycle slip detectableCycle slip repairableTurboEditEEMTurboEditEEM60(9, 7)(8.991, 6.991)(8.986, 6.991)YesYesYesYes120(6, 6)(5.999, 5.999)(6.006, 6.014)YesYesYesYes180(5, 5)(5.048, 5.048)(5.004, 5.014)YesYesYesYes240(3, 4)(3.027, 4.027)(3.018, 4.030)YesYesYesYes300(4, 3)(3.988, 2.988)(3.997, 2.996)YesYesYesYes360(1, 0)(1.010, 0.010)(0.998, 0.007)YesYesYesYes410(1, 1)(-3.548, 2.548)(1.010, 1.021)YesYesNoYes420(1, 1)(0.988, 0.988)(1.006, 0.987)YesYesYesYes480(77, 60)(77.042, 60.042)(76.986, 59.984)YesYesYesYes520(4, 3)(8.580, 6.580)(4.019, 3.033)YesYesNoYes660(3, 4)(2.971, 3.971)(3.014, 4.006)YesYesYesYes670(5, 5)(9.454, 8.454)(5.035, 5.015)YesYesNoYesTo 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. Table3 lists the detection and repair results. The results show that the EEM algorithm can detect and repair cycle slips that occur in continuous epochs.Table 3
Results of detecting and repairing continuous cycle slips in G01.
EpochSimulated cycle slip pairs (L1, L2)New method300(1, 1)(1.103, 1.096)301(1, 1)(1.096, 1.019)302(1, 1)(0.901, 1.000)303(9, 7)(8.996, 6.892)304(9, 7)(8.808, 6.993)


## 3.2. 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–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 ΔN1 and ΔN2. Figures 2(a)–2(d) show the variations in the IOD (sampling interval: 1 s). Figures 3(a)–3(h) show the variations in ΔN1 and ΔN2. 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 ΔN1 and ΔN2 are significantly small because the estimates of ΔN1 and ΔN2 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 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 time-difference carrier phase observations in equations (3) and (17) can be ignored for the proposed cycle slip detection and repair method.Figure 2
Elevation angles and variations in the IOD for (a) G01, (b) G08, (c) G15, and (d) G17 on DOY 145 of 2018.
(a)(b)(c)(d)Figure 3
Elevation angles,ΔN1, and ΔN2 for (a, b) G01, (c, d) G08, (e, f) G15, and (g, h) G17 on DOY 145 of 2018.
(a)(b)(c)(d)(e)(f)(g)(h)To investigate the performance of the proposed method at low elevation angles, we select G01 and G17 as examples (sampling interval: 1 s). From epoch 360 (GPS time 2018:05:25:00:01:00) to epoch 930 (GPS time 2018:05:25:00:15:30), 20 pairs of the simulated cycle slip pair (9, 7) are added to G01. The elevation angles in the period from epoch 600 to epoch 930 are all lower than 15°. From epoch 360 to epoch 570 (GPS time 2018:05:25:00:9:30), 8 pairs of the simulated cycle slip pair (1, 1) are added to G17. The elevation angles in the period from epoch 360 to epoch 570 are all lower than 15°. These cycle slip pairs are described in Figures4(a) and 4(b).Figure 4
Distribution of simulated cycle slip pairs at an interval of every 30 epochs.
(a)(b)The cycle slip repair results for G01 and G17 are shown in Figures5(a)–5(f). In each plot, the values of ΔN and the variation in the elevation angle are given. To visualize the results better, values of ΔN larger than 1.2 cycles were set as 1.2 cycles. These cycle slip epochs are indicated by vertical lines. 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 ΔN1 and ΔN2 are close to the simulated cycle slip pair (1, 1). This result verifies that these ΔN can be easily rounded to their nearest integers.Figure 5
Elevation angle,ΔN1, and ΔN2 for G01 and G17.


## 3.3. Standard Deviation of the Observation Correction Vector

Constructingχ2 requires the prior STD of the observation correction vector (V). To study the distribution of the STD, the STD values corresponding to each epoch and reference satellite are calculated for the Swarm-A observation data on DOY 145 according to sampling intervals of 1 s, 5 s, 10 s, 15 s, 30 s, and 60 s. According to the formula σ0=VTPV/n−3, the mean posterior STD values corresponding to different reference satellites in each epoch are calculated [1]. The statistical results related to the posterior STD are listed in Tables 4 and 5, which show that the posterior STD of the observation correction vector under different reference satellites is closely related to the data sampling interval, and the value of the posterior STD increases with a decrease in the sampling interval. At the same sampling 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]. 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.Table 4
STD values ofV under different basic satellites (unit: m).
GPS PRNSampling interval1 s5 s10 s15 s30 s60 sG010.00290.01070.02010.02900.05090.0794G020.00270.01070.02020.02960.05340.0819G030.00270.00970.01820.02640.04780.0788G050.00280.01080.02060.03010.05560.0858G060.00260.01010.01910.02770.04930.0801G070.00280.01060.02000.02920.05300.0796G080.00290.01100.02060.02930.05250.0792G090.00280.01050.01980.02870.04970.0788G100.00320.01120.01960.02800.05170.0839G110.00290.01110.02080.03020.05320.0760G120.00270.01050.01940.02810.05090.0807G130.00280.01080.02030.02950.05310.0794G140.00260.00940.01750.02510.04470.0765G150.00320.01180.02070.02970.05350.0815G160.00270.00980.01830.02660.04940.0826G170.00280.01090.02050.02950.05320.0860G180.00300.01080.02020.02920.05150.0854G190.00280.01100.02080.03030.05460.0880G200.00360.01230.02140.03020.05420.0862G210.00310.01110.02040.02960.05410.0869G220.00290.01030.01920.02790.04780.0809G230.00270.01010.01910.02780.04980.0803G240.00310.01150.02010.02850.05130.0817G250.00270.01010.01870.02710.04880.0795G260.00270.00990.01880.02730.04840.0816G270.00300.01050.01870.02670.04760.0760G280.00300.01150.02150.03090.05720.0863G290.00300.01080.01950.02810.05010.0793G300.00280.01020.01920.02790.05040.0770G310.00260.00930.01740.02520.04550.0750G320.00300.01080.01920.02690.04880.0808Table 5
Statistic results of the STD (unit: m).
Sampling intervalσ0MaxMinMean1 s0.00360.00260.00295 s0.01230.00930.010710 s0.02150.01740.019715 s0.03100.02520.028530 s0.05730.04480.051160 s0.08800.07500.0812


## 3.4. Performance of the EEM to Enhance Cycle Slip Errors

To verify the ability of the EEM to amplify cycle slip errors, a simulated outlier of 0.5 cycles and a cycle slip of 1 cycle are inserted into theL1 carrier phase of G11 in the observation epoch corresponding to GPS time 2018:05:25:00:06:30. Then, we calculate the χ2 values corresponding to the different 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.Table 6
Values ofχ2 for different basic satellites.
GPS PRNValue ofχ2No cycle slipOutlier of 0.5 cycles1 cycle slipG010.2740.23464.206G080.2739.71962.558G110.27278.374443.99G130.38133.1252.491G150.73945.65762.921G170.27332.95452.071G280.433.42553.079G301.14958.97779.322Figure 6
Values ofχ2 for different reference satellites under an outlier of 0.5 cycles.


## 4. Conclusion

A method for detecting and repairing onboard GPS receiver cycle slips is presented. This method combines single-frequency 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.



---

*Source: 8817626-2020-11-26.xml*