Algorithm for Real-Time Cycle Slip Detection and Repair for Low Elevation GPS Undifferenced Data in Different Environments

: Real-time cycle slip detection and repair is one of the key issues in global positioning system (GPS) high precision data processing and application. In particular, when GPS stations are in special environments, such as strong ionospheric disturbance, sea, and high-voltage transmission line interference, cycle slip detection and repair in low elevation GPS observation data are more complicated than those in normal environments. For low elevation GPS undifferenced carrier phase data in different environments, a combined cycle slip detection algorithm is proposed. This method uses the ﬁrst-order Gauss–Markov stochastic process to model the pseudorange multipath in the wide-lane phase minus narrow-lane pseudorange observation equation, and establishes the state equation of the wide-lane ambiguity with the pseudorange multipath as a parameter, and it uses the Kalman ﬁlter for real-time estimation and detects cycle slips based on statistical hypothesis testing with a predicted residual sequence. Meanwhile, considering there are certain correlations among low elevation, observation epoch interval, and ionospheric delay error, a second-order difference geometry-free combination cycle slip test is constructed that takes into account the elevation. By combining the two methods, real-time cycle slip detection for GPS low elevation satellite undifferenced data is achieved. A cycle slip repair method based on spatial search and objective function minimization criterion is further proposed to determine the correct solution of the cycle slips after they are detected. The whole algorithm is experimentally veriﬁed using the static and kinematic measured data of low elevation satellites under four different environments: normal condition, high-voltage transmission lines, dynamic condition in the sea, and ionospheric disturbances. The experimental results show that the algorithm can detect and repair cycle slips accurately for low elevation GPS undifferenced data, the difference between the ﬂoat solution and the true value for the cycle slip does not exceed 0.5 cycle, and the differences obey the normal distribution overall. At the same time, the wide-lane ambiguity and second-order difference GF combination sequence calculated by the algorithm is smoother, which give further evidence that the algorithm for cycle slip detection and repair is feasible and effective, and has the advantage of being immune to the special observation environments.


Introduction
High-quality global positioning system (GPS) observation data play a key role in obtaining high precision GPS positioning. However, when a satellite signal is blocked by an obstacle and fails to reach a receiver, or the satellite signal experiences temporary loss of lock due to external interference and harsh environment in which the receiver is located, an integer number of cycles in the phase observable jumps suddenly. This

Materials and Methods
In this section, the pseudorange multipath error and cycle slip detection performance of wide-lane phase minus narrow-lane pseudorange combination are firstly analyzed. Then, the influence of ionospheric residuals on the cycle slip detection method of firstorder difference geometry-free combination is discussed. On this basis, the combination algorithm of cycle slip real-time detection and repair is proposed, and the principle and steps of the whole algorithm are described in detail.

Wide-Lane Phase Minus Narrow-Lane Pseudorange (WL-NL) Combination with Pseudorange Multipath
For single epoch GPS satellites, the undifferenced observation equation of dualfrequency code pseudorange and carrier phase can be derived as follows [23][24][25]: In Equation (9), N w = N 1 − N 2 is the wide-lane ambiguity and λ w = c is the widelane carrier phase wavelength. Because the observation noise of the carrier phase is 2 mm, the maximum multipath error does not exceed 5 cm, so the item , the relationship between the combined pseudorange multipath error MP and wide-lane ambiguity N w can be obtained as It can be seen from Equation (10) that, when using it for cycle slip detection, the accuracy of cycle slip detection mainly comes from the pseudorange multipath error MP and measurement noise η; in an ideal situation, the mean value of the multipath error and measurement noise of pseudorange measurement should be close to 0, and the pseudorange measurement error M error introduced by MP and η can be approximately calculated by In Equations (11) and (12), M LW is the measurement error of the pseudorange combination including the wide-lane ambiguity, and N MW = λ w N w , N MW should be a constant when there is no cycle slip. Mean(M LW ) is the mean of M LW within the observation arc.
To intuitively analyze the magnitude of the pseudorange measurement error (M error ), the measured GPS data from the BJFS reference station are used as an example. The data were observed on 1 January 2014 and the sampling rate was 30 s. The code observation types are C/A and P2. The dual frequency observations corresponding to satellites G26 and G08 are selected and calculated according to Equation (12). The relationship between the pseudorange measurement error sequence of G26 and G08 and their elevation is shown in Figure 1. can be obtained as It can be seen from Equation (10) that, when using it for cycle slip detection, the accuracy of cycle slip detection mainly comes from the pseudorange multipath error and measurement noise ; in an ideal situation, the mean value of the multipath error and measurement noise of pseudorange measurement should be close to 0, and the pseudorange measurement error introduced by and can be approximately calculated by In Equations (11) and (12), is the measurement error of the pseudorange combination including the wide-lane ambiguity, and = , should be a constant when there is no cycle slip.
( ) is the mean of within the observation arc.
To intuitively analyze the magnitude of the pseudorange measurement error ( ), the measured GPS data from the BJFS reference station are used as an example. The data were observed on 1 January 2014 and the sampling rate was 30 s. The code observation types are C/A and P2. The dual frequency observations corresponding to satellites G26 and G08 are selected and calculated according to Equation (12). The relationship between the pseudorange measurement error sequence of G26 and G08 and their elevation is shown in Figure 1.  It can be seen from Figure 1 that, even when the satellite is at a high elevation (above 30 • ), its pseudorange measurement error can reach ±0.4 m; when the satellite is at a low elevation (5-20 • ), the effect of its pseudorange measurement error can reach ±1.0 m; therefore, when the wide-lane phase minus the narrow-lane pseudorange combination is used to calculate the ambiguity N MW , the pseudorange measurement error must be dealt with, otherwise the accuracy of ambiguity will be reduced.

Inaccurate Detection of Cycle Slip by Conventional MW Combination
The analysis in Section 2.1 indicates that the pseudorange multipath in the wide-lane phase minus narrow-lane pseudorange combination (i.e., MW combination) is the main factor affecting the accuracy of the wide-lane ambiguity calculation; when the conventional MW combination is used to detect cycle slips, the mean recursive algorithm is used to obtain the mean of the wide-lane ambiguity N i MW and the root mean square σ i for each epoch, and the cycle slip testing is performed by the following equation [12,18,26]: Because the accuracy of the P code pseudorange is approximately 0.3 m, and the MW combination is also affected by the pseudorange measurement error, this method is not sensitive to small cycle slips (two cycles or below), and thus inaccurate cycle slip detection may occur. The above mentioned G08 satellite is used as an example again; at the 50th and 250th epochs of the L1 carrier phase observation, cycle slips of injecting one and two cycles are simulated. The conventional MW combination cycle slip detection results are shown in Figure 2. It can be seen from the results in Figure 2 that, for low elevation satellites, the cycle slips identified by conventional MW combination for two cycles or less are invalid. It can be seen from Figure 1 that, even when the satellite is at a high elevation (above 30°), its pseudorange measurement error can reach ±0.4 m; when the satellite is at a low elevation (5°-20°), the effect of its pseudorange measurement error can reach ±1.0 m; therefore, when the wide-lane phase minus the narrow-lane pseudorange combination is used to calculate the ambiguity , the pseudorange measurement error must be dealt with, otherwise the accuracy of ambiguity will be reduced.

Inaccurate Detection of Cycle Slip by Conventional MW Combination
The analysis in Section 2.1 indicates that the pseudorange multipath in the wide-lane phase minus narrow-lane pseudorange combination (i.e., MW combination) is the main factor affecting the accuracy of the wide-lane ambiguity calculation; when the conventional MW combination is used to detect cycle slips, the mean recursive algorithm is used to obtain the mean of the wide-lane ambiguity and the root mean square for each epoch, and the cycle slip testing is performed by the following equation [12,18,26]: Because the accuracy of the P code pseudorange is approximately 0.3 m, and the MW combination is also affected by the pseudorange measurement error, this method is not sensitive to small cycle slips (two cycles or below), and thus inaccurate cycle slip detection may occur. The above mentioned G08 satellite is used as an example again; at the 50th and 250th epochs of the L1 carrier phase observation, cycle slips of injecting one and two cycles are simulated. The conventional MW combination cycle slip detection results are shown in Figure 2. It can be seen from the results in Figure 2 that, for low elevation satellites, the cycle slips identified by conventional MW combination for two cycles or less are invalid.

Cycle Slip Detection with Fused Kalman Filter and Hypothesis Testing
As GPS receivers cannot usually obtain code observation, in general, C/A code is used instead of code to perform the MW combination in Equation (10). The observation accuracy of C/A code is lower than that of code, which is greatly affected by pseudorange multipath, so the accuracy of the calculated wide-lane ambiguity value is reduced. Therefore, the pseudorange multipath error in Equation (11) must be separated from the wide-lane ambiguity . As the pseudorange multipath error is not constant and correlated with a certain time interval, the first-order Gauss-Markov process can be used to estimate the multipath error [27,28]:

Cycle Slip Detection with Fused Kalman Filter and Hypothesis Testing
As GPS receivers cannot usually obtain P 1 code observation, in general, C/A code is used instead of P 1 code to perform the MW combination in Equation (10). The observation accuracy of C/A code is lower than that of P 1 code, which is greatly affected by pseudorange multipath, so the accuracy of the calculated wide-lane ambiguity value is reduced. Therefore, the pseudorange multipath error MP in Equation (11) must be separated from the wide-lane ambiguity N MW . As the pseudorange multipath error is not constant and correlated with a certain time interval, the first-order Gauss-Markov process can be used to estimate the multipath error [27,28]: Remote Sens. 2021, 13, 2078 6 of 25 whereMP(t) is the change rate of multipath error and τ is the correlation time. Meanwhile, when the carrier phase has no cycle slip, the ambiguity is equal, and a state equation with the wide-lane ambiguity and pseudorange multipath error as the state parameters can be established: Using Equation (11) as the observation equation of the discrete system, combined with Equation (14), the following Kalman filter model can be established [29][30][31]: where X k and X k−1 are the state parameter vectors, φ k,k−1 is the state transition matrix , and A k is the coefficient matrix and A k = 1 1 . ω k and Q ω k are the state noise and its variance matrix, respectively, Q k = and v k and R v k are the observation noise and its variance matrix, respectively [32,33]. The Kalman filter recursive equation is used to achieve the separation of wide-lane ambiguity and pseudorange multipath error by epoch, and the best estimation of wide lane ambiguity is obtained. Meanwhile, assuming that there is no cycle slip in the carrier phase observation of the first k − 1 epochs and cycle slip occurs at the kth epoch, the Kalman filter model in Equation (17) becomes The prediction residual of the kth epoch is then It can be seen from Equation (18) that, when a cycle slip occurs, the prediction residual sequence V k will shift and no longer obey the zero-mean normal distribution [34,35]. Thus, cycle slip detection is performed by constructing a statistical hypothesis testing containing the predicted residual sequence, with the null hypothesis being that the cycle slip does not occur and the alternative hypothesis being that the cycle slip occurs, namely, The testing is performed based on the following criteria: Here, n is a weighting coefficient, which can be three (99.7% confidence level) or four (99.9% confidence level). Cycle slip detection is performed up to the last epoch using this criterion.

Detecting Cycle Slips with Second-Order Difference Geometry-Free Combination Considering Elevation
Based on the MW combination and using the fused method of the Kalman filter and hypothesis testing, cycle slip detection can be performed effectively. However, when the Remote Sens. 2021, 13, 2078 7 of 25 equal cycle slip occurs at frequencies L1 and L2, detection blind spots may appear. Therefore, a difference geometry-free (GF) combination taking into account the elevation factor is proposed here to correct the blind spots of the cycle slip detection method described in Section 2.3. By subtracting Equation (4) from Equation (3), the GF combination observation equation can be obtained to eliminate the clock error and tropospheric delay error [36,37]: Here, the multipath (mp 1 − mp 1 ) and observation noise (ε 1 − ε 2 ) on frequencies L1 and L2 are approximately equal, so the (mp 1 − mp 1 + ε 1 − ε 2 ) term on the right-hand side of the equation can be ignored, and the corresponding GF combination only contains the effects of the ionosphere and the ambiguity. When a cycle slip occurs in the carrier phase observation of the ith epoch, differencing the above equation between adjacent epochs yields the first-order difference GF combination [18,25] containing the ionospheric residual I e−res = (I e (i) − I e (i − 1)) and cycle slip information (λ 2 ∆N 2 − λ 1 ∆N 1 ): Under normal circumstances, the ∆L GF sequence appears as a rapidly changing curve due to cycle slip. Meanwhile, it can be seen from Equation (21) that the cycle slip information is affected by the ionospheric residual I e−res , and especially when GPS satellites are in an active period of the ionosphere, the ionospheric residual term will inevitably cause disturbances in the ∆L GF sequence of Equation (21), resulting in inaccurate cycle slip detection [20,21]. Therefore, to minimize the effect of the ionospheric residual term, considering the correlation between ∆L GF and the elevation (especially the low elevation), and the epoch time interval, the second-order difference GF combination cycle slip detection quantity is constructed that takes into account the elevation, as expressed by the following equation: The root mean square of the first i epochs is recursively calculated by Eqution (23): Then, the cycle slip detection threshold condition is expressed by the following equation: The cycle slip detection of low elevation GPS carrier phases in different environments is performed based on the combination of the above two methods.

Cycle Slip Repair Based on Spatial Search and Objective Function Minimization Criterion
The methods proposed in Sections 2.3 and 2.4 can achieve real-time detection of cycle slips. By analyzing Equations (17) and (18), it can be seen that the value of the cycle slips can be expressed by residuals. Equation (22) shows that the cycle slip value can be expressed by the difference between the epochs of L GF . However, when the two frequencies have equal cycle slips, or the cycle slips of the two frequencies are proportional to the wavelength, a detection blind zone that is not sensitive to cycle slips manifests when the two methods are used alone. Therefore, in actual application, the two methods need to be combined effectively. Together, they make up for each others' shortcomings, and the occurring frequency and value of cycle slips can be determined by Equation (25).
Cycle slips are expressed in integer numbers, while the solution given by Equation (22) is the float solution of cycle slips [38]; to avoid the effect of rounding errors caused by direct rounding, a method using spatial search and objective function minimization criteria is proposed. To obtain the integer solution of cycle slip, the specific steps are as follows: 1.
Round the float value calculated by Equation (25) as the initial value of the cycle slip. 2.
The effects of low elevation and ionospheric disturbance period are fully considered, the initial value of cycle slip is taken as the center, ±5 cycles is the search range, and one cycle is the search step, in order to form a cycle slip candidate combination. 3.
Use the cycle slip candidate combination constructed in Step 2 to repair the carrier phase observation in real time, then re-check it through the cycle slip detection methods proposed in this study, and store the cycle slip combination that the two cycle slip detection methods meet the condition of no cycle slip simultaneously.

4.
Use the cycle slip combination that meets the detection condition to repair the carrier phase observation, and substitute the repaired carrier phase into the objective function, so that the set of cycle slip combinations that makes the objective function in Equation (26) meet the minimization criterion is the correct cycle slip solution.

Results
In this section, the GPS data in four different environments are firstly introduced, then the pseudorange multipath error of low elevation satellite is analyzed. Based on the experimental schemes of small cycle slips, large cycle slips and special cycle slips were added into the simulation at different epochs of GPS data, and the cycle slip detection and repair method proposed in this study is verified and analyzed.

Measured GPS Data in Four Different Environments
To analyze and verify the feasibility of the cycle slip detection and repair methods proposed in this paper, the following four different environments and different sampling rates of GPS static and kinematic measured data were used for testing.

Normal Environment
The GPS data of the P090 reference station measured on 3 January 2010 were selected for testing. The sampling rate of the data was 15 s. The satellite G09 is a low elevation satellite. The elevation and pseudorange measurement error sequence are shown in Figure 3, from which it can be seen that, even if the IGS reference station is in a good location environment, the influence of low elevation on the pseudorange measurement error is still large.

High-Voltage Transmission Lines Environment
The measured data of the GPS receiver at station TB06, which is located directly under high-voltage transmission lines and is closest to the transmission tower (see Figure 4), were selected for testing. The data were taken on 2 June 2016 with a sampling rate of 1 s. G03 is a low elevation satellite; its elevation and pseudorange measurement error sequence are shown in Figure 5. Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 25

High-Voltage Transmission Lines Environment
The measured data of the GPS receiver at station TB06, which is located directly under high-voltage transmission lines and is closest to the transmission tower (see Figure  4), were selected for testing. The data were taken on 2 June 2016 with a sampling rate of 1 s. G03 is a low elevation satellite; its elevation and pseudorange measurement error sequence are shown in Figure 5.

High-Voltage Transmission Lines Environment
The measured data of the GPS receiver at station TB06, which is located directly under high-voltage transmission lines and is closest to the transmission tower (see Figure  4), were selected for testing. The data were taken on 2 June 2016 with a sampling rate of 1 s. G03 is a low elevation satellite; its elevation and pseudorange measurement error sequence are shown in Figure 5.

High-Voltage Transmission Lines Environment
The measured data of the GPS receiver at station TB06, which is located directly under high-voltage transmission lines and is closest to the transmission tower (see Figure  4), were selected for testing. The data were taken on 2 June 2016 with a sampling rate of 1 s. G03 is a low elevation satellite; its elevation and pseudorange measurement error sequence are shown in Figure 5.

Dynamic Environment in the Sea Area
The measured kinematic data of GPS receiver installed on a buoy in a sea region of Indonesia were selected for testing. The GPS trajectory on the buoy is shown in Figure 6. The measured kinematic data of GPS receiver installed on a buoy in a sea region of Indonesia were selected for testing. The GPS trajectory on the buoy is shown in Figure 6. The data were taken on 3 March 2010 with a sampling rate of 5 s. G10 is a low elevation satellite; the corresponding elevation and pseudorange measurement error sequence are shown in Figure 7.

Ionospheric Disturbance Environment
Cai et al. [18] show that the solar flare intensity on 15 April 2001 was 14.4, which was the highest flare intensity in the last maximum solar activity cycle from the year 2000 to 2002. The ionosphere was thus in a period of strong disturbance. The actual measurement data of IGS reference station AMC2 in this environment were used for testing, and its sampling rate was 30 s. G26 is a low elevation satellite; the corresponding elevation and pseudorange measurement error sequence are shown in Figure 8.

Dynamic Environment in the Sea Area
The measured kinematic data of GPS receiver installed on a buoy in a sea region of Indonesia were selected for testing. The GPS trajectory on the buoy is shown in Figure 6. The data were taken on 3 March 2010 with a sampling rate of 5 s. G10 is a low elevation satellite; the corresponding elevation and pseudorange measurement error sequence are shown in Figure 7.

Ionospheric Disturbance Environment
Cai et al. [18] show that the solar flare intensity on 15 April 2001 was 14.4, which was the highest flare intensity in the last maximum solar activity cycle from the year 2000 to 2002. The ionosphere was thus in a period of strong disturbance. The actual measurement data of IGS reference station AMC2 in this environment were used for testing, and its sampling rate was 30 s. G26 is a low elevation satellite; the corresponding elevation and pseudorange measurement error sequence are shown in Figure 8.

Ionospheric Disturbance Environment
Cai et al. [18] show that the solar flare intensity on 15 April 2001 was 14.4, which was the highest flare intensity in the last maximum solar activity cycle from the year 2000 to 2002. The ionosphere was thus in a period of strong disturbance. The actual measurement data of IGS reference station AMC2 in this environment were used for testing, and its sampling rate was 30 s. G26 is a low elevation satellite; the corresponding elevation and pseudorange measurement error sequence are shown in Figure 8.
It can be seen from Figures 3, 5, 7 and 8 that the pseudorange measurement error of each low elevation satellite varies greatly in the different special environmental conditions, which will inevitably affect the accuracy of conventional MW and GF combination cycle slip detection. Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 25 It can be seen from Figures 3, 5, 7, and 8 that the pseudorange measurement error of each low elevation satellite varies greatly in the different special environmental conditions, which will inevitably affect the accuracy of conventional MW and GF combination cycle slip detection.

Validation of Cycle Slip Detection and Repair Results
To test the feasibility of the cycle slip detection and repair algorithm developed in this study, large cycle slips, small cycle slips, and special cycle slip pairs were added into the simulation at different epochs of the low elevation measured data sequence under the above four different environmental conditions. The cycle slip pairs are shown in Table 1.

Validation of Cycle Slip Detection and Repair Results
To test the feasibility of the cycle slip detection and repair algorithm developed in this study, large cycle slips, small cycle slips, and special cycle slip pairs were added into the simulation at different epochs of the low elevation measured data sequence under the above four different environmental conditions. The cycle slip pairs are shown in Table 1. The corresponding data of each low elevation in the above four different environmental conditions are used to test the cycle slip detection and repair method proposed in this paper. The cycle slip detection results are shown in Figures 9-12 The corresponding data of each low elevation in the above four different environmental conditions are used to test the cycle slip detection and repair method proposed in this paper. The cycle slip detection results are shown in Figures 9-12. (a) epoch of 50th            It can be seen from Figures 9-12 that the cycle slip combination detection algorithm proposed in this paper accurately detects the epoch position in which the cycle slips occur, and is not affected by the sampling rate of the observation data. Thus, this algorithm is feasible for cycle slip detection of low elevation dual-frequency undifferenced GPS phase data in different special environments, and it is easy to implement.
At the same time, to evaluate the cycle slip repair method based on the spatial search and objective function minimization criterion proposed in this paper, Table 2 shows the cycle slip repair results of each low elevation under the aforementioned four different It can be seen from Figures 9-12 that the cycle slip combination detection algorithm proposed in this paper accurately detects the epoch position in which the cycle slips occur, and is not affected by the sampling rate of the observation data. Thus, this algorithm is feasible for cycle slip detection of low elevation dual-frequency undifferenced GPS phase data in different special environments, and it is easy to implement.
At the same time, to evaluate the cycle slip repair method based on the spatial search and objective function minimization criterion proposed in this paper, Table 2 shows the cycle slip repair results of each low elevation under the aforementioned four different environments. The statistics of the difference between the float solution and its respective true value are shown in Figure 13. environments. The statistics of the difference between the float solution and its respecti true value are shown in Figure 13.

Environment and Satellite Number Epoch Float Solution Objective Function Value Cycle Slip Search Resul
Normal, G09  It can be seen from Table 2 and Figure 13 that, when the cycle slip repair meth proposed in this paper is used to obtain the initial float solution, the difference betwe It can be seen from Table 2 and Figure 13 that, when the cycle slip repair method proposed in this paper is used to obtain the initial float solution, the difference between the float solution and the true value for the cycle slip does not exceed 0.5 cycles, and the differences obey the normal distribution overall. This indicates that the float solution has high accuracy, and the cycle slip results obtained through spatial search and objective function minimization criterion are accurate and reliable, which further demonstrates the effectiveness of the cycle slip repair method.

Discussion
According to the experimental results in Section 3.2, for different sampling rates of low elevation GPS phase data, it is necessary to effectively combine the two cycle slip detection methods proposed in this paper to achieve accurate cycle slip detection and repair in four different environments. Only at the 500th epoch in Figure 11, the cycle slip detection method with the fused Kalman filter and hypothesis testing based on MW combination cannot detect the small cycle slip of one cycle. The reason is that the acquisition environment of the dynamic data is located in the sea, which is one of the main error sources of multipath error, and as there is no P1 code observation in the data, the observation noise is large when the C/A code is used to replace the P code for MW combination. However, the cycle slip is accurately detected by the second-order differential GF method considering the elevation, and the correct cycle slip value is finally found. Therefore, for small cycle slips such as (1,0), if the observation data contains P code, the performance of the combined cycle slip detection algorithm proposed in this paper will be perfect.
To further evaluate and analyze the performance of the cycle slip detection algorithm, satellite G10 was used under the dynamic environment in the sea area, and satellite G26 in the ionospheric disturbance environment as examples to give the ambiguity sequence after separating the pseudorange multipath based on Kalman filter, and the second-order difference GF combination sequence taking into account the satellite elevation; the corresponding results are shown in Figures 14 and 15 the float solution and the true value for the cycle slip does not exceed 0.5 cycles, and the differences obey the normal distribution overall. This indicates that the float solution has high accuracy, and the cycle slip results obtained through spatial search and objective function minimization criterion are accurate and reliable, which further demonstrates the effectiveness of the cycle slip repair method.

Discussion
According to the experimental results in Section 3.2, for different sampling rates of low elevation GPS phase data, it is necessary to effectively combine the two cycle slip detection methods proposed in this paper to achieve accurate cycle slip detection and repair in four different environments. Only at the 500th epoch in Figure 11, the cycle slip detection method with the fused Kalman filter and hypothesis testing based on MW combination cannot detect the small cycle slip of one cycle. The reason is that the acquisition environment of the dynamic data is located in the sea, which is one of the main error sources of multipath error, and as there is no P1 code observation in the data, the observation noise is large when the C/A code is used to replace the P code for MW combination. However, the cycle slip is accurately detected by the second-order differential GF method considering the elevation, and the correct cycle slip value is finally found. Therefore, for small cycle slips such as (1,0), if the observation data contains P code, the performance of the combined cycle slip detection algorithm proposed in this paper will be perfect.
To further evaluate and analyze the performance of the cycle slip detection algorithm, satellite G10 was used under the dynamic environment in the sea area, and satellite G26 in the ionospheric disturbance environment as examples to give the ambiguity sequence after separating the pseudorange multipath based on Kalman filter, and the second-order difference GF combination sequence taking into account the satellite elevation; the corresponding results are shown in Figures 14 and 15, respectively.   It can be seen from Figures 14 and 15 that the wide-lane ambiguity in the MW combination estimated by the Kalman filter is smoother than the wide-lane ambiguity sequence determined by the conventional MW combination, which indicated that the Kalman filter method used in this paper can effectively separate the pseudorange multipath error, improve the reliability of the wide-lane ambiguity estimation, and improve the cycle slip detection capability of the MW combination. At the same time, the second-order difference GF combination method that considers the satellite elevation is more stable than the conventional first-order difference GF combination. In particular, in the case of ionospheric disturbance, the second-order difference GF combination that takes into account the low satellite elevation is better in reducing the effect of ionospheric error on the cycle slip testing quantities, and further improves the GF combination cycle slip detection capability under special environments.
Analysis of experimental results further shows that, on the basis of MW combination, the wide-lane ambiguity can be separated from the pseudorange multipath through the Kalman filter, thereby improving the accuracy of float solutions of wide-lane ambiguity and providing a reference data processing method for subsequent ambiguity fixes. At the same time, the second-order difference GF combination taking into account the low elevation can effectively reduce the ionospheric residual error, and has the advantage of being immune to the strong ionospheric disturbance observation environment.

Conclusions
In this study, for the cycle slip detection problem in GPS undifferenced carrier phase data at low elevation in different environments, on the basis of MW combination, cycle slip detection is carried out by constructing cycle slip test quantities, fused Kalman filter, hypothesis testing, and a second-order difference GF combination is proposed that takes into account the elevation factor. The experimental test results of data with different sampling rates at various low elevation under four special environments show that the cycle slip combination detection algorithm proposed in this paper can detect the epoch positions of small and large cycle slips accurately in real time, and it is feasible for cycle slip detection and easy to implement.
After the cycle slips are accurately detected, the experimental results of cycle slip repair show that the proposed method based on spatial search and objective function minimization criterion can achieve accurate results for the cycle slip integer solution. This method is effective for repairing cycle slips of GPS low elevation satellite data. It can be seen from Figures 14 and 15 that the wide-lane ambiguity in the MW combination estimated by the Kalman filter is smoother than the wide-lane ambiguity sequence determined by the conventional MW combination, which indicated that the Kalman filter method used in this paper can effectively separate the pseudorange multipath error, improve the reliability of the wide-lane ambiguity estimation, and improve the cycle slip detection capability of the MW combination. At the same time, the second-order difference GF combination method that considers the satellite elevation is more stable than the conventional first-order difference GF combination. In particular, in the case of ionospheric disturbance, the second-order difference GF combination that takes into account the low satellite elevation is better in reducing the effect of ionospheric error on the cycle slip testing quantities, and further improves the GF combination cycle slip detection capability under special environments.
Analysis of experimental results further shows that, on the basis of MW combination, the wide-lane ambiguity can be separated from the pseudorange multipath through the Kalman filter, thereby improving the accuracy of float solutions of wide-lane ambiguity and providing a reference data processing method for subsequent ambiguity fixes. At the same time, the second-order difference GF combination taking into account the low elevation can effectively reduce the ionospheric residual error, and has the advantage of being immune to the strong ionospheric disturbance observation environment.

Conclusions
In this study, for the cycle slip detection problem in GPS undifferenced carrier phase data at low elevation in different environments, on the basis of MW combination, cycle slip detection is carried out by constructing cycle slip test quantities, fused Kalman filter, hypothesis testing, and a second-order difference GF combination is proposed that takes into account the elevation factor. The experimental test results of data with different sampling rates at various low elevation under four special environments show that the cycle slip combination detection algorithm proposed in this paper can detect the epoch positions of small and large cycle slips accurately in real time, and it is feasible for cycle slip detection and easy to implement.
After the cycle slips are accurately detected, the experimental results of cycle slip repair show that the proposed method based on spatial search and objective function minimization criterion can achieve accurate results for the cycle slip integer solution. This method is effective for repairing cycle slips of GPS low elevation satellite data.

Data Availability Statement:
The datasets analyzed in this study are managed by the College of Geology Engineering and Geomatics, Chang'an University and can be available on request from the corresponding author.