Real-time Triple-Frequency Cycle Slip Detection and Repair Algorithm Based on the Optimal Fixing Probability

Global Navigation Satellite System (GNSS) is now being speedily expanded to our daily life, but the positioning precision still can hardly meet the demands of many applications, such as approaching landing system on airports. Due to the development of GNSS, triple-frequency signals are now available which can contribute to positioning precision. Positioning precision cannot be improved by triple-frequency carrier phases until cycle slips are detected and repaired. Traditional cycle slip detection and repair algorithms choose detection combinations with long wavelength, weak ionospheric delay and small combination noise separately. However, these three conditions cannot be satisfied simultaneously. In this paper, these three conditions are not considered separately. On the contrary, the eventual fixing probability of cycle slip is set as the optimal goal to determine the three detection combinations. The combined ionospheric delay and noise in cycles can be regarded as bias and variance respectively. The proposed algorithm has been tested on observations with simulated and real cycle slips. The results show that the proposed algorithm can detect and repair even single cycle slips in real time effectively.


Introduction
The Global Navigation Satellite System (GNSS) is being widely used in many fields for precise positioning. Because of pseudorange noise, carrier phases are the main observation for highly accurate positioning. However, because of a low signal-to-noise ratio or jamming, receivers may lose track of navigation signals so that cycle slips occur. Without careful consideration of such cycle slips, the precision of positioning is seriously affected [1][2][3][4][5][6].
Over the past decades, many cycle slip detection and repair methods have been developed. One method to detect and repair cycle slips using double differenced observations can be found in Chen et al. (2016) [1]. This method needs two receivers at least. The Hatch-Melbourne-Wuebbena (HMW) (Hatch, 1982) linear combination has been widely applied to the cycle slip detection in the dual-frequency receiver [2]. In addition, the ionospheric residual combination, which uses the difference of carrier phase observations on two frequencies, can also be used. However, it has several insensitive combinations (Cai et al. 2013) [3]. Together with the ionospheric combination, the HMW linear combination was used in the TurboEdit algorithm to detect and repair cycle slips (Blewitt, 1990) [4]. Liu (2011) used the ionospheric total electron contents rate (TECR) and the HMW linear combination to uniquely determine the cycle slip on both L1 and L2 frequencies [5]. Banville et al. (2013) proposed a geometry-based approach with the rigorous handling of the ionosphere [6]. Cai et al. (2013) used a forward and backward moving window averaging algorithm and a second-order, time-difference phase ionospheric residual algorithm to detect and repair cycle slips [3].
In the case of triple-frequencies, the additional frequency results in more linear combinations with longer wavelengths, weaker ionospheric delays, and smaller noise [7]. Hatch et al. (2000) introduced the benefits of the third frequency signal on cycle slip correction [8]. Dai et al. (2009) used three linear combinations and the LAMBDA algorithm to detect cycle slip on triple-frequency data [9]. De Lacy et al. (2012) defined five geometry-free linear combinations used in three steps for real-time cycle slip detection and repair on triple-frequency GNSS data [10]. Huang et al. (2015) used two geometry-free phase combinations and one geometry-free pseudorange minus phase linear combination to detect and correct cycle slip in real time [11]. Zhao et al. (2015) used three combinations, namely extra-wide lane, wide lane, and narrow lane, to determined cycle slips sequentially in three cascaded steps and results showed that this method performed well under high ionospheric activity [12]. Zeng et al. (2018) used triple-frequency combinations to detect and repair cycle slip under ionospheric disturbance with BDS data [13].
To enlarge the fixing probability of cycle slip detection, researches above choose detection combinations with long wavelength, weak ionospheric delay and smaller combination noise separately. However, these three conditions cannot be satisfied simultaneously. There is a compromise between these three conditions, but there are almost no researches about which three linear combinations are the optimal yet. In this research, a real-time triple-frequency cycle slip detection and repair algorithm based on optimal fixing probability is proposed. We do not consider these three conditions separately when we determine the three detection combinations. On the contrary, we take the fixing probability as the optimal goal to determine the three detection combinations. The combined ionospheric delay and noise in cycles can be regarded as bias and variance respectively when we calculate the fixing probability. When those three detection combinations are determined, cycle slips on the original carrier phase observations can be uniquely identified.
The following section describes the choices of those three combinations in detail and introduces the method of recovering cycle slip on each frequency based on the chosen three combinations. Then, the BDS and GPS triplefrequency data with simulated and real cycle slips are used to test the performance of the proposed algorithm. Finally, conclusions are drawn.

Basic Observations
The pseudorange and carrier observation equations can be expressed as: L 1 satisfies Gaussian distribution, mean is N i,j,k and variance is 2σ 2 L1 . It can be expressed as follows: where N(μ, σ 2 ) expresses Gaussian distribution. The cycle slips are fixed to an integer by rounding the float estimation. Thus, the threshold is set as 0.5 cycles and the probability of successfully determining cycle slip, which is defined as fixing probability (F.P.) in this research, is expressed as follows: The coefficient of pseudorange observations (l 1 , l 2 , l 3 ) can be determined when the coefficient of carrier phase observations (i,j,k) is given: Then the optimal solution is as follows:    To guarantee that the noise is not too large, we vary the range of (i,j,k) from -5 to +5. Table 1 lists the five best pseudorange/carrier phase combinations according to the fixing probability (F.P.) when , σ 2 ρ3 and σ 2 Φ 3 are set as 2, 0.3 m and 0.003 m, respectively. Figure 1 shows the results of fixing probability as a function of pseudorange standard deviation on f 3 and ratio , when the carrier phase standard deviation is set as 0.003 m. As the figure shows, when the value of σ ρ3 is growing, the fixing probabilities of all combinations, expect for (0, -1, 1), decrease no matter the value of the ratio  for both BDS and GPS. As to (0, -1, 1), the fixing probability is almost 100% even when the value of σ ρ3 is 0.5 m. So, we choose (0, -1, 1) as the first carrier phase combination for both BDS and GPS. It should be mentioned that we do not need to repair cycle slips on each frequency when using the first carrier phase combination to provide a geometrical reference for the second carrier phase combination which will be introduced in the following part. We just need to repair the combined cycle slips on the first carrier phase combination. When we repair the combined cycle slips, cycle slips on each frequency will not cause discontinuity on the geometry of the first carrier phase combination. For example, when cycle slips exist in the form of (0, N, N) (N can be any integers), this type of cycle slip pair will not cause discontinuity on the geometry of the first carrier phase combination so that we do not need to repair it when we use the first carrier phase combination to provide geometry to the second carrier phase combination.

Second Detection Combination
The second detection combination in cycles consisting of the first fixed carrier phase combination and the second carrier phase combination can be expressed as follows where (i,j,k) is the coefficient of the first carrier phase combination and (m,n,t) is the coefficient of the second carrier phase combination. The second detection combination eliminates the geometry but the ionospheric delay still exists and the value of ionospheric delay in cycles is as follows , , In accordance with the law of variance-covariance propagation, the variance of L 2 in cycles can be expressed as follows 2   3   2   2  2  2  , , , The difference L 2 between continuous epochs, which is expressed as L 2 , can be calculated as follows L 2 satisfies Gaussian distribution with average N m,n,t + I L2 and variance is 2σ 2 L2 . It can be expressed as follows The fixing probability can be expressed as follows To enlarge the fixing probability of cycle slip detection, I L 2  and σ 2 L 2 should be as small as possible. But both I L 2  and σ 2 L 2 are determined by (m,n,t) so the two cannot achieve minimum value simultaneously when (m,n,t) is an integer combination. So, we do not consider I L 2  and σ 2 L 2 separately and we take the eventual fixing probability as our optimizing goal. Just as the first detection combination, the range of (m,n,t) varies from -5 to +5 in this research. As to TECR, Liu et al. (2009) showed that in the equatorial region of Hong Kong the ionospheric slant total electron content rate (TECR) was about 0.01 TECU/s during quiet ionosphere periods and to 0.03 TECU/s during disturbed periods [15]. We set the TECR as 0.03 TECU/s to assess the validity of different carrier phase combinations under high ionospheric activity. Table 2 lists the ten best carrier phase combinations according to the fixing probability for BDS and GPS; the carrier phase standard deviation is set as 0.003 m. Figure 2 shows the fixing probability varying as a function of TECR. As the figure shows, and in agreement with Tab. 2, the fixing probabilities of the ten carrier phase combinations are larger than 99.95% for BDS and 99.1% for GPS even when TECR is 0.03 TECU/s. So, these combinations can determine cycle slips efficiently under high ionospheric activity. Because they have the same fixing probability, we cannot determine the best combination during the second stage. We set all of them as alternative combinations and provide a geometrical reference for the third carrier phase combination. Just as with the first carrier phase combination, the second carrier phase combination also does not need to repair cycle slips on each frequency  when it is used to provide a geometrical reference for the third carrier phase combination. We just need to guarantee that the combined cycle slips have been repaired.

Third Detection Combination
The third detection combination in cycles consisting of the second fixed carrier phase combination and the third carrier phase combination can be expressed as follows: where (m, n, t) is the coefficient of the second alternative carrier phase combination listed in Tab. 2 and (u, v, w) is the coefficient of the third carrier phase combination. Just as for the second detection combination, geometry has been eliminated in the third detection combination, but the ionospheric delay still exists. In accordance with the law of variance-covariance propagation, the variance of L 3 in cycles can be expressed as follows: (28) The third detection combination is sensitive to the ionospheric delay and we use the second-order time-difference method to reduce the effect of ionospheric delay as has been introduced by (Cai et al. 2013) [3]. In order to recover cycle slips on each frequency, the matrix composed by (i, j, k), (m, n, t) and (u, v, w) should be reversible, which means that (i, j, k), (m, n, t) and (u, v, w) should be linearly independent. Under this condition, we take, again, the eventual fixing probability as our optimizing goal as we did with the second detection combination. (m, n, t) is chosen from Tab. 2 and (u, v, w) varies from -5 to +5. Assuming that cycle slip occurs at epoch k, and there are no cycle slips or the cycle slips have been repaired at epoch k -1 and k -2, the second-order time difference of L 3 between three continuous epochs, which expressed as L 3 , can be calculated as: The fixing probability can be expressed as follows:   Table 3 lists the best carrier phase combinations according to the fixing probability when the carrier phase standard deviation is set as 0.003 m. As the table shows, the fixing probability can be up to 99.958% for BDS when the second carrier phase combination is chosen as (1 -5 4) and the third carrier phase combination is chosen as (-3 0 4) or (4 -5 0), and 99.996% for GPS when the second carrier phase combination is chosen as (1 4 -5) and the third carrier phase combination is chosen as (-3 2 2). The combined noise is enlarged by two times compared to the first-order time-difference combination, so it is necessary to assess fixing the probability as it varies with the carrier phase standard deviation. Figure 3 shows the fixing probability of the best carrier phase combinations varying with carrier phase standard deviation. The fixing probability decreases when the carrier phase standard deviation increases. But the fixing probability is larger than 99.1% for BDS and 99.7% for GPS even when the carrier phase standard deviation is 0.004 m. The results show that the second-order time-difference method is efficient when applied to fix cycle slip using the carrier phase combinations listed in Tab. 3.

Recover Cycle Slip on Each Frequency
For BDS, there are two optimal combinations listed in Tab. 3. The first and the second carrier phase combinations are selected as (0 -1 1) and (1 -5 4) and the third carrier phase combination can be selected as (-3 0 4) or (4 -5 0) which has the same fixing probability. We select (-3 0 4) as the third carrier phase combination in this research. Actually, as Tables 1, 2 and 3 show, when we determinate each optimal combination, all factors affecting the fixing probability, including ionospheric delay and combined noise, are in units of cycles. Although the wavelength and the combined noise of (-3 0 4) and (4 -5 0) are different in meters, their combined noise in cycles is the same so that their fixing probabilities are the same as those shown in Fig. 3. As a result, we can select (4 -5 0) also as the third carrier phase combination in this research. To be brief, we choose (-3 0 4) as the third combination. For GPS, there is one optimal combination listed in Tab. 3 and the first, the second, and the third carrier phase combination are selected as (0 -1 1), (1 4 -5) and (-3 2 2), respectively. Then the detection equation for BDS can be expressed as follows: and the detection equation for GPS can be expressed as follows: int(x) means that x is rounded to the nearest integer.

Data Analysis
The active state of the ionosphere can be judged by Kp index [16]. We have analyzed the Kp index during recent years retrieved from http://omniweb.gstc.nasa.gov/form/dx1.html. We found that the Kp indexes were high on March 17, 2013 and September 8, 2017 compared to other days. Figure 4 shows the geomagnetic Kp index on those two days.
In order to assess the validity of the proposed method, triple-frequency BDS and GPS data from stations DYNG, FINA, JFNG, CUT0 and HARB on September 8, 2017 are analyzed. Table 4 lists the antenna types and receiver types    Figure 6 shows the first-order and the second-order time-differenced TEC variations of C04 at FINA, C09 at JFNG, C11 at DYNG, G01 at HARB, and G03 at CUT0 on September 8, 2017, and G01 observed at JFNG on March 17, 2013. The first-order ionospheric delay variations for all involved satellites were fairly significant. However, the trend component is not obvious in the second-order timedifference TEC variations for all satellites involved. The second-order time-difference TEC variations are clustered around zero over the entire observation time. This indicates that the impact of ionospheric delay can be significantly reduced using the second-order time-difference method.  Figure 7 shows the L 1 , L 2 and L 3 of the satellites on each station calculated by observations from September 8, 2017. Satellites (C04, C09, C11, G01, G03) observed from (FINA, JFNG, DYNG, HARB, CUT0), respectively, are analyzed here. The figure shows the noise of those three combinations varies with elevation, but the detected results are lower than 0.5 cycles even when elevation is 10 degrees. This means that our algorithm performs well during the entire observation time. Figure 8 shows the comparison of our algorithm and that of Liu et al. (2018).  G01 observed from JFNG on March 17, 2013 are analyzed. They used three optimal linearly independent geometryfree pseudorange minus phase combinations to detect and repair cycle slips. Their algorithm performs well when the elevation is high. But as the right column shows, when the elevation is low, there are some results detected exceeding 0.5 cycles due to the large pseudorange noise. As to our algorithm, detected results are between in 0.5 cycles for the whole observation time even when elevation is low. The results show that our algorithm is more reliable when elevation is low.

Simulated Cycle Slip Test
In order to assess the performance of the proposed method, we simulated several cycle slip pairs on the original carrier phase observations. Large cycle slips can be detected and repaired easily so that we do not discuss them in this research. The small cycle slips range from 1 to 2 and they can occur on any frequencies. Several random cycle slip pairs are simulated on the original carrier phase observations to test the algorithm.   Figure 9 shows the random cycle slip detection results on the three combinations for each satellite. Cycle slip pairs (0, 1, 1) and (1, 1, 1) are simulated to observations. In addition, several other random cycle slip pairs are also simulated. As Figure 9 shows, the cycle slips, despite the small cycle size of 1, cause significant jump on the three combinations. The third detection combination, which is insensitive to ionospheric delay, can detect cycle slips efficiently during the entire observation time when the ionospheric activity is high. The accurate determination of the cycle slips on the three combinations easily resolves the original observations. Table 5 summarizes all the detection results in the test. As listed, the proposed method correctly detected and repaired all the random cycle slips.

Real Cycle Slip Test
To test the performance of the proposed algorithm, we use it to detect and repair real cycle slips using observations of BDS C10 at JFNG on March 17, 2013. The method by Cai et al. (2013) to detect dual-frequency cycle slips performs well when the ionosphere is active [3]. We use this method to detect B1, B2 and B1, B3 dual-frequency cycle slip, respectively. The results show that there occurs a (-2, -2, -2) cycle slip pair at epoch 2075. Then we use our algorithm to process the same observation. Figure 10 shows the results. The detection result of the third combination is -1.76 and it can be rounded to -2 at epoch 2075. The detection results of the first and the second combination do not exceed 0.5 cycles so that they can be rounded to 0 at epoch 2075. Then the cycle slip pair can be calculated according to (34) and the detected result is (-2, -2, -2) which is exactly the same as the result detected by the method of Cai et al. (2013). Consequently, the reliability of our algorithm is further tested.

Conclusions
A real-time triple-frequency cycle slip detection and repair algorithm under high ionospheric activity based on optimal fixing probability has been proposed. The effect of ionospheric delay and combined noise are not considered separately and we take the eventual fixing probability of cycle slips as optimizing goal. We analyze the fixing probability of cycle slips quantitatively in detail and choose three detection combinations in accordance with fixing probability. In addition, the proposed algorithm forms the three detection combinations in three cascaded steps so that we can guarantee that each combination is optimal in accordance with the fixing probability. The proposed algorithm has been tested on real 30-second triple-frequency static observations of BDS and GPS on September 8, 2017 and March 17, 2013 when the ionospheric activity was high. Simulated cycle slips and real cycle slips are tested. The results show that the proposed algorithm can detect and repair cycle even single slips in real time effectively.