A Modified TurboEdit Cycle-Slip Detection and Correction Method for Dual-Frequency Smartphone GNSS Observation

Recently, some smartphone manufacturers have subsequently released dual-frequency GNSS smartphones. With dual-frequency observations, the positioning performance is expected to be significantly improved. Cycle-slip detection and correction play an important role in high-precision GNSS positioning, such as precise point positioning (PPP) and real-time kinematic (RTK) positioning. The TurboEdit method utilizes Melbourne–Wübbena (MW) and phase ionospheric residual (PIR) combinations to detect cycle-slips, and it is widely used in the data processing applications for geodetic GNSS receivers. The smartphone pseudorange observations are proved to be much noisier than those collected with geodetic GNSS receivers. Due to the poor pseudorange observation, the MW combination would be difficult to detect small cycle-slips. In addition, some specific cycle-slip combinations, where the ratio of cycle-slip values at different carrier frequencies is close to the frequency ratio, are also difficult to be detected by PIR combination. As a consequence, the traditional TurboEdit method may fail to detect specific small cycle-slip combinations. In this contribution, we develop a modified TurboEdit cycle-slip detection and correction method for dual-frequency smartphone GNSS observations. At first, MW and PIR combinations are adopted to detect cycle-slips by comparing these two combinations with moving-window average values. Then, the epoch-differenced wide-lane combinations are used to estimate the changes of smartphone position and clock bias, and the cycle-slip is identified by checking the largest normalized residual whether it exceeds a predefined threshold value. The process of estimation and cycle-slip identification is implemented in an iterative way until there is no over-limit residual or there is no redundant measurement. At last, the cycle-slip values at each frequency are estimated with the epoch-differenced wide-lane and ionosphere-free combinations, and the least-square ambiguity decorrelation adjustment (LAMBDA) method is adopted to further obtain an integer solution. The proposed method has been verified with 1 Hz dual-frequency smartphone GNSS data. The results show that the modified TurboEdit method can effectively detect and correct even for specific small cycle-slip combinations, e.g., (4, 3), which is difficult to be detected with the traditional TurboEdit method.


Introduction
In May 2016, Google released an application programming interface (API), which makes the developers access global navigation satellite system (GNSS) raw measurements, including pseudorange, phase, and doppler measurements, in Android 7.0 and up [1]. In the past years, single-frequency GNSS chipsets are dominant in the smartphone market. With a large noise level of smartphone single-frequency GNSS measurements, the positioning accuracy and reliability are thus unsatisfactory, especially under complex conditions [2,3].
Innovatively, in May 2018, the world's first dual-frequency GNSS smartphone produced by Xiaomi was launched. It is equipped with a Broadcom BCM47755 chipset, which is a dual-frequency (L1/L5+E1/E5a) GNSS chipset [4][5][6]. Subsequently, some smartphone manufacturers have recently released dual-frequency GNSS smartphones, e.g., P30/P30 Pro from Huawei, Galaxy Note 10/10+ from Samsung, One Plus 7 from OnePlus, etc., (https://developer.android.com). A series of data quality assessment and positioning algorithms based on the smartphone dual-frequency GNSS observations were carried out [5][6][7][8][9][10]. According to these published researches, the positioning performance with smartphone dual-frequency GNSS observations has been improved effectively thanks to the newly added L5/E5a signals. In GNSS high-precision positioning, such as real-time kinematic (RTK) positioning and precise point positioning (PPP), the GNSS phase measurements are necessarily applied, and they are desired to be continuously tracked. However, the discontinuities of GNSS phase observations, widely called cycle-slips, usually happen due to signal obstacles, large environment noise, strong electromagnetic interference, etc., [11]. Therefore, cycle-slip detection and correction are very crucial for phased-based GNSS high-precision positioning.
As an important and challenging issue in GNSS data preprocessing, lots of algorithms have been proposed for cycle-slip detection and correction in the past decades. These methods include phase ionospheric residual method [12], Kalman filtering based method [13], polynomial fitting method [14,15], high-order difference method [16], and wavelet transform method [17,18]. Each of these methods has its own limitations [19]. The TurboEdit algorithm, initially developed by Blewitt (1990), is an easy and effective cycle-slip detection method [20,21]. In addition, it has been widely applied in the high-precision GNSS data processing software, such as GIPSY, PANDA, and Bernese 5.x [22][23][24]. The TurboEdit algorithm utilizes Melbourne-Wübbena (MW) combination and phase ionospheric residual (PIR) combination to detect cycle-slips. The PIR combination is insensitive to special cycle-slips where the ratio of cycle-slip values at each carrier frequency is close to the frequency ratio, but the MW combination can benefit by detecting these special cycle-slips. However, the efficiency of detecting cycle-slips with the MW combination will decrease when applying it to the smartphone GNSS measurements due to its poor pseudorange. Miao et al. (2011) proposed a modified TurboEdit method for satellite-borne GPS observations by introducing a satellite elevation weighted factor to MW combination [19]. However, this method does not change the situation that the noise of MW combination is too large to detect small cycle-slips for smartphone GNSS measurements. In order to reduce the effect of the noise of MW combination, Cai et al. (2013) applied a forward and backward moving window averaging algorithm. Obviously, this method is only suitable for GNSS data post-processing [25]. In a word, the TurboEdit method needs to be advanced for cycle-slip detection of smartphone GNSS observations.
We propose a modified TurboEdit cycle-slip detection and correction method for dual-frequency smartphone GNSS observations. At first, the MW and PIR combinations are adopted to detect cycle-slips by comparing these two combinations with moving-window average values. Most of the cycle-slips are expected to be detected in this step. Then, the epoch-differenced wide-lane combinations are used to estimate the changes of smartphone position and clock bias by applying least-square adjustment. The largest normalized residual is sorted out and checked whether it exceeds the threshold value. If the largest normalized residual is bigger than the predefined threshold, both dual-frequency phase observations of the corresponding satellite are marked as cycle-slip. The process of estimation and cycle-slip identification is implemented in an iterative way until there is no over-limit residual or there is no redundant measurement. At last, the float cycle-slip values at two frequencies can be estimated by using the epoch-differenced wide-lane and ionosphere-free combinations, and then the least-square ambiguity decorrelation adjustment (LAMBDA) method is applied to obtain an integer solution of cycle-slips.
The rest of this paper is organized as follows. In Section 2, the modified TurboEdit cycle-slip detection and correction method is discussed in detail. In Section 3, the performance is evaluated with Sensors 2020, 20, 5756 3 of 14 the collected dual-frequency smartphone GNSS measurements. Finally, conclusions are summarized in the last section.

Traditional TurboEdit Algorithm
The GNSS raw code and phase measurements read as [26]: where P i , L i denote code and phase measurements; the subscript i represents the carrier frequency; ρ is the geometric distance between receiver and satellite; dt r , dt s are the clock offsets at the receiver and satellite; trop, ion represent the tropospheric and ionospheric delays; c, λ i are the speed of light and the wavelength, N i is the ambiguity; ε P i , ε L i denote the code and phase measurement noises including multipath errors.
In the traditional TurboEdit method, the cycle-slip detection is based on MW and PIR combinations. With dual-frequency smartphone GNSS observations, these two combinations can be formed as follows [27]: where L MW , L PIR represent the MW and PIR combinations, respectively; f 1 and f 5 are the carrier frequencies; λ WL = c/( f 1 − f 5 ) and N WL = N 1 − N 5 denote the wide-lane wavelength and wide-lane ambiguity.
To detect cycle-slips with MW and PIR combinations, a moving-window averaging filter is utilized to calculate the average and standard deviation (STD) values: where N WL , L PIR are the average value of N WL and L PIR in the moving-window with a window size of m; k and k − 1 denote the present and previous epochs; σ N WL , σ L PIR represent the STD of N WL and L PIR , respectively. The satellite will be marked as cycle-slip when any one of the following conditions is satisfied: The noise level of smartphone pseudorange observations is much larger than that of pseudorange observations from geodetic GNSS receivers. This will lead to a larger noise level of MW combination. σ N WL can be approximatively calculated in the following equation: At present, the noise of smartphone pseudorange observations is generally higher than 1.5 m [8], 4·σ N WL ≈ 5.70 cycles. Therefore, it is difficult to detect wide-lane cycle-slips of 1-5 cycles with MW Sensors 2020, 20, 5756 4 of 14 combination. According to Equation (2), when the cycle-slip combination (∆N 1 , ∆N 5 ) nearly meets the relationship of ∆N 1 /∆N 5 = f 1 / f 5 , it cannot be detected by PIR combination. As a result, the traditional TurboEdit method may fail to detect specific small cycle-slip combinations, e.g., (4·k, 3·k), k < 6.

Modified TurboEdit Method
Except for specific small cycle-slip combinations, most of the cycle-slips can be detected with the traditional TurboEdit method. Therefore, the traditional TurboEdit method is first applied in our modified TurboEdit method. In order to further detect specific small cycle-slip combinations, the epoch-differenced wide-lane combinations are used to estimate the changes of smartphone position and clock bias, and the cycle-slip is identified by checking the largest normalized residual whether it exceeds the predefined threshold value. The wide-lane combination can be written as: where where ∆ stands for the epoch-differenced operator. The variation of wide-lane ambiguity ∆N WL is expected to zero for continuous phase observations. The satellite position and clock offset can be calculated with broadcast ephemeris. Hence, the unknown parameters remain the changes of smartphone position and clock bias. All epoch-differenced wide-lane combinations, of which satellites are not marked as cycle-slip with the traditional TurboEdit method, are adopted to estimate the unknown parameters. The measurement equation can be expressed as: where y WL is the observed-minus-computed epoch-differenced wide-lane observation residual vector; A is the design matrix of the unknown parameters; x = [∆dt r ] is the unknown parameter vector for the smartphone at a static scenario; and x = [∆x r ∆y r ∆z r ∆dt r ] for the smartphone at a kinematic scenario. ε WL is the measurement error vector and its covariance matrix is Q. Assuming that the phase observations at two frequencies are independent at the same noise level with a zenith STD of σ L , the covariance matrix Q can be derived as: Here, ; el represents the satellite elevation with respect to the smartphone. The least-square adjustment can be used to estimate the changes of smartphone position and clock bias:x where P = Q −1 denotes the weight matrix. The corrected residuals can be computed as: Sensors 2020, 20, 5756

of 14
Assuming that the number of observed satellites is n and the number of unknown parameters is m, the formal STD can be calculated as:σ The cycle-slips will be identified by checking the largest normalized residual whether it exceeds a predefined threshold [28,29]: wherev k represents the k th corrected residual; σv k =σ 0 · Q[k][k] and Q[k][k] denotes k th diagonal element of Q; η is the predefined threshold value. If that is the case, the satellite with the largest normalized residual will be then isolated from the least-square adjustment. An iterative way is implemented until the condition of Equation (14) is not satisfied or there is no redundancy. Once all possible cycle-slips are detected, the cycle-slip values at each frequency will be estimated with all epoch-differenced wide-lane and ionosphere-free combinations including both cycle-slip and non-cycle-slip observations, and an integer estimation of cycle-slip values can be resolved with the LAMBDA method.
In case that the cycle-slip has been detected for a specific satellite, the estimation parameters for cycle-slip values should be introduced to the corresponding measurement equation. As a result, the unknown parameters are the changes of smartphone position, clock bias, and cycle-slip parameters. Hence, the measurement equation formed by the epoch-differenced wide-lane observations can be written as: where x represents the vector composed by the change of smartphone position and clock bias; ∆N is the cycle-slip vector; A WL and B WL are the corresponding design matrices; ε WL is the measurement error vector. The ionosphere-free combination formed by the carrier phase measurement at two frequencies can be expressed as [30,31]: where stands for the ionosphere-free ambiguity in meters; ε L IF = α·ε L 1 + (1 − α)·ε L 5 denotes measurement noise of the ionosphere-free phase combination. Once the epoch-differencing is employed to high-rate smartphone GNSS ionosphere-free observations, the tropospheric delay would be nearly eliminated. Consequently, the epoch-differenced ionosphere-free observation can be written as: Broadcast ephemeris can be adopted to calculate the satellite position and clock offset. Similarly, all epoch-differenced ionosphere-free observations at the current epoch are employed to estimate the unknown parameters. The measurement equation can be expressed as: where A IF and B IF are the design matrices with respect to x and ∆N; ε IF is the measurement error vector of ionosphere-free phase combination. Combining Equations (15) and (18), we yield: Sensors 2020, 20, 5756 6 of 14 The float estimation of cycle-slips can be calculated by applying the least-square adjustment to Equation (19). With float estimation of cycle-slips, the LAMBDA technique can be employed to obtain a integer solution of cycle-slips [32]. Furthermore, the integer estimation of cycle-slips is validated by applying the ratio test [33,34].
The procedure of the modified TurboEdit method is shown as Figure 1. At first, the MW and PIR combinations are used to detect the cycle-slip satellite by a satellite. Most of the cycle-slips will be detected in this step except for some specific small cycle-slip combinations. Then, the epoch-differenced wide-lane combination is applied to estimate the changes of position and clock bias. The cycle-slip is identified by checking the largest normalized residual whether it exceeds the predefined threshold value. If yes, the satellite with the largest normalized residual will be isolated from the least-square adjustment. The process of estimation and cycle-slip identification is implemented in an iterative way until there is no over-limit residual or there is no redundant measurement. Finally, the cycle-slip values at each frequency are estimated with the epoch-differenced wide-lane and ionosphere-free combinations, and the LAMBDA method is adopted to further obtain an integer cycle-slip solution.
Sensors 2020, 20, x FOR PEER REVIEW 6 of 14 threshold value. If yes, the satellite with the largest normalized residual will be isolated from the least-square adjustment. The process of estimation and cycle-slip identification is implemented in an iterative way until there is no over-limit residual or there is no redundant measurement. Finally, the cycle-slip values at each frequency are estimated with the epoch-differenced wide-lane and ionosphere-free combinations, and the LAMBDA method is adopted to further obtain an integer cycle-slip solution.

Experimental Results and Analysis
In this section, the performance of the modified TurboEdit method will be presented by using the real dual-frequency smartphone GNSS data. The raw GNSS dataset, including GPS L1/L5 and Galileo E1/E5a code and phase measurements, was collected with a Xiaomi Mi-8 smartphone under static condition from 06:40:00 to 07:30:00, on August 20, 2019, in GPS time. The GNSS observations were recorded at a sample interval of 1 s. Two data subsets, from 06:42:00 to 06:46:59 and from 06:50:00 to 06:54:59, are extracted from the collected dataset. One is used to analyze the performance of detecting and correcting small cycle-slips, and the other is used to verify the effectiveness of detecting and correcting specific small cycle-slip combinations. These two data subsets are assumed to be collected in kinematic scenes during the following process of cycle-slip detection and correction. The cutoff elevation is set to 7 degrees. Figure 2 shows the number of visible GPS and Galileo satellites

Experimental Results and Analysis
In this section, the performance of the modified TurboEdit method will be presented by using the real dual-frequency smartphone GNSS data. The raw GNSS dataset, including GPS L1/L5 and Galileo E1/E5a code and phase measurements, was collected with a Xiaomi Mi-8 smartphone under static condition from 06:40:00 to 07:30:00, on August 20, 2019, in GPS time. The GNSS observations were recorded at a sample interval of 1 s. Two data subsets, from 06:42:00 to 06:46:59 and from 06:50:00 to Sensors 2020, 20, 5756 7 of 14 06:54:59, are extracted from the collected dataset. One is used to analyze the performance of detecting and correcting small cycle-slips, and the other is used to verify the effectiveness of detecting and correcting specific small cycle-slip combinations. These two data subsets are assumed to be collected in kinematic scenes during the following process of cycle-slip detection and correction. The cutoff elevation is set to 7 degrees. Figure 2 shows the number of visible GPS and Galileo satellites with dual-frequency observations.

Small Cycle-slip Combination Detection and Correction
Using the first data subset, six cases are simulated to analyze the performance of detecting small cycle-slips. In each case, the simulated cycle-slip is added to raw GNSS phase observations manually. These six simulated small cycle-slip cases include 1, 3, 5 cycles on the G26 L1 carrier phase signal at the 151th epoch, and the same cycle-slips on the E01 E1 carrier phase signal at the same epoch. The initial STD values of and are set to 3.0 cycles and 0.028 m, respectively, which are approximately calculated according to the noise level of smartphone measurements at present [8]. The time series of ( ) − ( − 1) in six different cases are shown in Figures 3 and 4. The blue points represent the variation of ( ) − ( − 1) and the red curves denote the cycle-slip boundaries formed by ±4 . It is clear that the variation of ( ) − ( − 1) is well bounded by the red curves. Hence, no cycle-slip is identified by MW combination according to the cycle-slip judgment condition in Equation (5). The failure to identify the small cycle-slips is because of the poor smartphone pseudorange. Shown in Figures 5 and 6 are the cycle-slip detection results for the G26 and E01 by using the PIR combination. The blue points present the variation of ( ) − ( − 1) and red curves are the corresponding boundaries yielded by ±4 . Obviously, the simulated cycle-slips are uniquely identified by PIR combination. Therefore, all simulated small cycle-slips can be correctly detected with the traditional TurboEdit method, and this also indicates that most of the cycle-slips can be identified by applying the traditional TurboEdit method.

Small Cycle-Slip Combination Detection and Correction
Using the first data subset, six cases are simulated to analyze the performance of detecting small cycle-slips. In each case, the simulated cycle-slip is added to raw GNSS phase observations manually. These six simulated small cycle-slip cases include 1, 3, 5 cycles on the G26 L1 carrier phase signal at the 151th epoch, and the same cycle-slips on the E01 E1 carrier phase signal at the same epoch. The initial STD values of N WL and L PIR are set to 3.0 cycles and 0.028 m, respectively, which are approximately calculated according to the noise level of smartphone measurements at present [8]. The time series of N WL (k) − N WL (k − 1) in six different cases are shown in Figures 3 and 4. The blue points represent the variation of N WL (k) − N WL (k − 1) and the red curves denote the cycle-slip boundaries formed by ±4σ N WL . It is clear that the variation of N WL (k) − N WL (k − 1) is well bounded by the red curves. Hence, no cycle-slip is identified by MW combination according to the cycle-slip judgment condition in Equation (5). The failure to identify the small cycle-slips is because of the poor smartphone pseudorange. Shown in Figures 5 and 6 are the cycle-slip detection results for the G26 and E01 by using the PIR combination. The blue points present the variation of L PIR (k) − L PIR (k − 1) and red curves are the corresponding boundaries yielded by ±4σ N PIR . Obviously, the simulated cycle-slips are uniquely identified by PIR combination. Therefore, all simulated small cycle-slips can be correctly detected with the traditional TurboEdit method, and this also indicates that most of the cycle-slips can be identified by applying the traditional TurboEdit method.

Small Cycle-slip Combination Detection and Correction
Using the first data subset, six cases are simulated to analyze the performance of detecting small cycle-slips. In each case, the simulated cycle-slip is added to raw GNSS phase observations manually. These six simulated small cycle-slip cases include 1, 3, 5 cycles on the G26 L1 carrier phase signal at the 151th epoch, and the same cycle-slips on the E01 E1 carrier phase signal at the same epoch. The initial STD values of and are set to 3.0 cycles and 0.028 m, respectively, which are approximately calculated according to the noise level of smartphone measurements at present [8]. The time series of ( ) − ( − 1) in six different cases are shown in Figures 3 and 4. The blue points represent the variation of ( ) − ( − 1) and the red curves denote the cycle-slip boundaries formed by ±4 . It is clear that the variation of ( ) − ( − 1) is well bounded by the red curves. Hence, no cycle-slip is identified by MW combination according to the cycle-slip judgment condition in Equation (5). The failure to identify the small cycle-slips is because of the poor smartphone pseudorange. Shown in Figures 5 and 6 are the cycle-slip detection results for the G26 and E01 by using the PIR combination. The blue points present the variation of ( ) − ( − 1) and red curves are the corresponding boundaries yielded by ±4 . Obviously, the simulated cycle-slips are uniquely identified by PIR combination. Therefore, all simulated small cycle-slips can be correctly detected with the traditional TurboEdit method, and this also indicates that most of the cycle-slips can be identified by applying the traditional TurboEdit method.                Once cycle-slips have been detected, all epoch-differenced wide-lane and ionosphere-free observations at the same epoch are adopted to estimate the changes of smartphone clock bias and cycle-slip parameters. The differences between the cycle-slip float estimations and the real values are 0.021, 0.018, 0.023 cycles for G26 and 0.023, 0.026, 0.023 cycles for E01. The LAMBDA method is further employed to obtain an integer solution of cycle-slips, which is presented in Table 1. The integer solution is validated by applying the ratio test with a threshold of 3.0. The ratio values are 438.12, 453.80, 408.80 for G26 and 382.60, 368.79, 395.97 for E01, respectively. Therefore, the integer solution of cycle-slip can be accepted. It is clear that by checking the largest residual of exceeding the predefined threshold can also uniquely identify the small cycle-slips (range from 1 to 5 cycles) for dual-frequency smartphone GNSS observations.   Once cycle-slips have been detected, all epoch-differenced wide-lane and ionosphere-free observations at the same epoch are adopted to estimate the changes of smartphone clock bias and cycle-slip parameters. The differences between the cycle-slip float estimations and the real values are 0.021, 0.018, 0.023 cycles for G26 and 0.023, 0.026, 0.023 cycles for E01. The LAMBDA method is further employed to obtain an integer solution of cycle-slips, which is presented in Table 1. The integer solution is validated by applying the ratio test with a threshold of 3.0. The ratio values are 438.12, 453.80, 408.80 for G26 and 382.60, 368.79, 395.97 for E01, respectively. Therefore, the integer solution of cycle-slip can be accepted. It is clear that by checking the largest residual of exceeding the predefined threshold can also uniquely identify the small cycle-slips (range from 1 to 5 cycles) for dual-frequency smartphone GNSS observations.

Specific Small Cycle-Slip Combination Detection and Correction
Most of the cycle-slips can be detected by applying the traditional TurboEdit method, but it may fail to identify some specific small cycle-slip combinations. For instance, if the cycle-slip combinations satisfy the conditions that |∆N 1 − ∆N 5 | is small than six cycles and (∆N 1 /∆N 5 ) ≈ ( f 1 / f 5 ), the cycle-slip may not be identified by the traditional TurboEdit method. In this part, a group of specific small cycle-slips, which are (4, 3), (12,9), and (20, 15) are simulated. These three specific small cycle-slip combinations are chosen because they are difficult to be detected by applying the traditional TurboEdit method alone. Based on the second data subset, six cases are also simulated by adding the cycle-slip combinations of (4, 3), (12,9), and (20, 15) on the G26 L1 and L5 carrier phase signals at the 151th epoch, and adding the same cycle-slip combinations on the E01 E1 and E5a carrier phase signal at the same epoch.
The proposed method is employed to detect and correct cycle-slip combinations with the simulated six sub-datasets. Figures 9 and 10 show the time series of N WL (k) − N WL (k − 1) and the cycle-slip boundaries formed by ±4σ N WL . The blue points still represent the variation of N WL (k) − N WL (k − 1) and the red curves denote the corresponding boundaries. The values of L PIR (k) − L PIR (k − 1) and the cycle-slip boundaries formed by ±4σ N PIR at different epochs are shown in Figures 11 and 12. It is apparent that the values of N WL (k) − N WL (k − 1) and L PIR (k) − L PIR (k − 1) are well constrained by the cycle-slip boundaries. In other words, both MW combination and PIR combination have failed to identify these specific small cycle-slip combinations. Therefore, we can conclude that these specific small cycle-slip combinations cannot be detected by the traditional TurboEdit method with dual-frequency smartphone GPS/Galileo observations. satisfy the conditions that |∆ − ∆ | is small than six cycles and (∆ ∆ ⁄ ) ≈ ( ⁄ ), the cycleslip may not be identified by the traditional TurboEdit method. In this part, a group of specific small cycle-slips, which are (4, 3), (12,9), and (20,15) are simulated. These three specific small cycle-slip combinations are chosen because they are difficult to be detected by applying the traditional TurboEdit method alone. Based on the second data subset, six cases are also simulated by adding the cycle-slip combinations of (4, 3), (12,9), and (20, 15) on the G26 L1 and L5 carrier phase signals at the 151th epoch, and adding the same cycle-slip combinations on the E01 E1 and E5a carrier phase signal at the same epoch. The proposed method is employed to detect and correct cycle-slip combinations with the simulated six sub-datasets. Figures 9 and 10 show the time series of ( ) − ( − 1) and the cycle-slip boundaries formed by ±4 . The blue points still represent the variation of ( ) − ( − 1) and the red curves denote the corresponding boundaries. The values of ( ) − ( − 1) and the cycle-slip boundaries formed by ±4 at different epochs are shown in Figures 11 and 12. It is apparent that the values of ( ) − ( − 1) and ( ) − ( − 1) are well constrained by the cycle-slip boundaries. In other words, both MW combination and PIR combination have failed to identify these specific small cycle-slip combinations. Therefore, we can conclude that these specific small cycle-slip combinations cannot be detected by the traditional TurboEdit method with dual-frequency smartphone GPS/Galileo observations.   slip may not be identified by the traditional TurboEdit method. In this part, a group of specific small cycle-slips, which are (4, 3), (12,9), and (20,15) are simulated. These three specific small cycle-slip combinations are chosen because they are difficult to be detected by applying the traditional TurboEdit method alone. Based on the second data subset, six cases are also simulated by adding the cycle-slip combinations of (4, 3), (12,9), and (20, 15) on the G26 L1 and L5 carrier phase signals at the 151th epoch, and adding the same cycle-slip combinations on the E01 E1 and E5a carrier phase signal at the same epoch. The proposed method is employed to detect and correct cycle-slip combinations with the simulated six sub-datasets. Figures 9 and 10 show the time series of ( ) − ( − 1) and the cycle-slip boundaries formed by ±4 . The blue points still represent the variation of ( ) − ( − 1) and the red curves denote the corresponding boundaries. The values of ( ) − ( − 1) and the cycle-slip boundaries formed by ±4 at different epochs are shown in Figures 11 and 12. It is apparent that the values of ( ) − ( − 1) and ( ) − ( − 1) are well constrained by the cycle-slip boundaries. In other words, both MW combination and PIR combination have failed to identify these specific small cycle-slip combinations. Therefore, we can conclude that these specific small cycle-slip combinations cannot be detected by the traditional TurboEdit method with dual-frequency smartphone GPS/Galileo observations.   . Specific small cycle-slip combination detection results with PIR combination for G26. Figure 11. Specific small cycle-slip combination detection results with PIR combination for G26.
The epoch-differenced wide-lane combinations, formed by GPS L1/ L5 and Galileo E1/E5a phase observations, are utilized to estimate the changes of smartphone positions and clock bias by applying the least-square adjustment. Figures 13 and 14 show the time series of the residuals of G26 and E01.
The over-limit residuals are 0.639, 2.122, 3.535 m for G26 and 0.742, 2.187, 3.227 m for E01 at the 151th epoch. That is to say, the simulated specific small cycle-slip combinations (4,3), (12,9), (20,15) can be clearly detected by checking the largest normalized residual of exceeding the predefined threshold.   Figures 13 and 14 show the time series of the residuals of G26 and E01. The over-limit residuals are 0.639, 2.122, 3.535 m for G26 and 0.742, 2.187, 3.227 m for E01 at the 151th epoch. That is to say, the simulated specific small cycle-slip combinations (4, 3), (12,9), (20,15) can be clearly detected by checking the largest normalized residual of exceeding the predefined threshold.     Figures 13 and 14 show the time series of the residuals of G26 and E01. The over-limit residuals are 0.639, 2.122, 3.535 m for G26 and 0.742, 2.187, 3.227 m for E01 at the 151th epoch. That is to say, the simulated specific small cycle-slip combinations (4, 3), (12,9), (20,15) can be clearly detected by checking the largest normalized residual of exceeding the predefined threshold.  After all cycle-slip combinations are detected, the float cycle-slip values at each frequency can be estimated with the epoch-differenced wide-lane and ionosphere-free observations. The biases between the cycle-slip float estimations and the real values are 0.021, 0.020, 0.023 cycles for G26 and 0.020, 0.013, 0.020 cycles for E01. By further employing the LAMBDA method, an integer solution of cycle-slips can be obtained which is shown in   After all cycle-slip combinations are detected, the float cycle-slip values at each frequency can be estimated with the epoch-differenced wide-lane and ionosphere-free observations. The biases between the cycle-slip float estimations and the real values are 0.021, 0.020, 0.023 cycles for G26 and 0.020, 0.013, 0.020 cycles for E01. By further employing the LAMBDA method, an integer solution of cycle-slips can be obtained which is shown in It illustrates that the modified TurboEdit method can correctly detect and repair the simulated cycle-slip combinations. Compared to the traditional TurboEdit method, the modified TurboEdit method is more effective in detecting specific small cycle-slip combinations for dual-frequency smartphone GNSS observations. In summary, given the poor smartphone pseudorange observations, it is difficult to detect a specific small cycle-slip combination with the traditional TurboEdit method when the ratio of cycle-slips at each frequency is close to that of carrier frequencies. The experiment results show that all the simulated cycle-slips are successfully detected and corrected by the modified TurboEdit method, and the modified TurboEdit method is effective even for specific small cycle-slip combinations, which is difficult to be identified with the traditional TurboEdit method.

Conclusions
Cycle-slip detection and correction play an important role in high-precision GNSS positioning. It is a great challenge to detect and correct specific small cycle-slip combinations for smartphone GNSS observations. In this contribution, a modified TurboEdit method has been proposed for dual-frequency smartphone GNSS observations. In the modified TurboEdit method, MW and PIR combinations are first utilized to detect the cycle-slip. Most of the cycle-slips are expected to be detected in this step. Then, the epoch-differenced wide-lane combinations are used to estimate the changes of smartphone position and clock bias, and the cycle-slip is identified by checking the largest normalized residual whether it exceeds the predefined threshold value. The process of estimation and cycle-slip identification is implemented in an iterative way until there is no over-limit residual or there is no redundant measurement. Once all cycle-slips are detected, the float values of cycle-slips at each frequency are estimated with the epoch-differenced wide-lane and ionosphere-free combinations, and LAMBDA is further employed to obtain an integer solution of cycle-slips.
The performance of the modified TurboEdit method has been evaluated with dual-frequency smartphone GNSS observations collected by Xiaomi Mi-8. The results indicate that the modified TurboEdit method can effectively detect and correct the simulated cycle-slip combinations. Compared with the traditional TurboEdit method, the modified TurboEdit method is more effective in detecting specific small cycle-slip combinations.