GPS/BDS triple-frequency cycle slip detection and repair based on moving window global search method

Cycle slip detection and repair are crucial steps in achieving high accuracy in Global Navigation Satellite System (GNSS) data processing. The use of Global Positioning System (GPS) and BeiDou Navigation Satellite System (BDS) triple frequency observations allows for more accurate detection and repair of cycle slips compared to single or dual frequency. This study presents a moving window global search method by selecting three sets of combined coefficients to construct geometry-free (GF) models to minimize the influence of the ionosphere, using a moving window to update the standard deviation of cycle slip estimation, applying the "3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upsigma$$\end{document}σ" criterion to constrain the range, and utilizing a global search method to detect and repair triple-frequency cycle slips. Through five sets of 1 Hz GNSS data experiments, the results demonstrate the effectiveness of this method in determining the position and size of triple-frequency cycle slips while avoiding multi-value problems. The detection success rate for GPS ranges from 98.0 to 100.0%, while BDS ranges from 92.0 to 100.0%. On average, GPS achieves a detection rate of 99.2%, and BDS reaches 96.7%, which is 0.8% and 1.8% higher than the direct rounding method, respectively. Compared to existing methods, it is also effective for the vast majority of small cycle slips within 2 cycles.

precision positioning (PPP).Currently, the TurboEdit method [19][20][21][22][23] is a widely used dual-frequency cycle slip detection and repair method that uses Melbourne-Wübbena (MW) and Geometry-Free (GF) models.Although this method is widely used in the practical world of engineering, it is not without its own problems.For instance, the dual-frequency combination wavelength is still relatively short, and at times, the pseudo-range error can significantly affect the detection outcome.At the same time, it is insensitive to cycle slips of the same size or special combination cycle slips on the dual-frequency, and it is prone to missed or false detections.Recently, scholars 24 have made new attempts to detect cycle slips through graphical structural constraints.With the development of navigation systems, many studies have been conducted on the combination of GPS and BDS multi-frequency 25 .Currently, the combination coefficient method is widely used to detect triple frequency cycle slips.The selection of combination coefficients must ensure the integer characteristics of cycle slips, such as long wavelengths, low noise impact, and minimal influence from the ionosphere.Dai et al. 26 propose a dual GF phase model for detecting and repairing triple-frequency cycle slips.The model utilises the least squares method to remove ambiguity correlation and determine the optimal value of cycle slips.While the model can detect most cycle slips, there are still a few insensitive ones that cannot be detected or repaired.Some scholars [27][28][29][30][31][32] use a combination of pseudorange-phase and GF models to detect triple-frequency cycle slips in BDS.The optimal combination coefficient is selected and possible cycle slips are searched for.Although most cycle slips can be detected, small cycle slips may be difficult to detect due to observation noise or special combination cycle slips.
This article presents a novel moving window global search algorithm for detecting and repairing GPS/BDS triple-frequency cycle slips.The algorithm utilizes three sets of frequency coefficient combinations, implements a moving window approach, applies integer constraints to cycle slip estimation, updates its range, and employs the global search method to simultaneously detect and repair triple-frequency cycle slips.Compared to other methods, this approach offers a simpler model and superior detection capabilities.One key advantage of this method is its ability to eliminate the need to consider the integer characteristics of combined cycle slips.It also effectively separates small cycle slips from observation noise, accurately constrains the real-time standard deviation of cycle slip sequences, and determines the position and size of each frequency's cycle slip without external conditions.Furthermore, it avoids the issue of multiple values for cycle slips.

Method Cycle slip valuation
Given the triple-frequency f i (i = 1, 2, 3) of GNSS, the equations for pseudo-range and phase observations can be expressed as follows: where, P i represents pseudo-range observation value, i represents carrier wavelength, ϕ i represents the carrier phase observation, ρ represents the geometric range, ε P i and ε ϕ i represent the observation noise of pseudo-range and phase respectively, N i is the integer ambiguity, η i = f 2 1 /f 2 i , I 1 is the first-order ionospheric term parameters of f 1 .
By combining Eqs.(1) and ( 2) and ignoring the influence of the ionosphere, we can calculate the difference between adjacent epochs t 1 and t 2 .This will allow us to obtain the estimated cycle slip at epoch t 2 : where, If the receiver receives pseudo-range P i on f i (i = 1, 2, 3) , it can take the average of the received multi-fre- quency pseudo-range observations instead of their respective P i .If the receive can receive the GNSS P-code, replace the pseudo-range with the P-code.The triple-frequency cycle slip estimation represented by Eq. ( 3) can be expressed as:

Combined GF model
Linear combinations are employed in multi-frequency GNSS observations to mitigate the impact of certain parameters.The combination can be represented based on the triple-frequency phase observations: where α represents the combination coefficient.Expand the above Eq.( 5) as: (1) If α 1 + α 2 + α 3 = 0 , the influence of geometric distance ρ can be eliminated, and the difference between adjacent epochs can be obtained as: If ( 1 η 1 + 2 η 2 + 3 η 3 )�I 1 ≈ 0 , the influence of ionospheric residuals can be ignored.Combining the integer feature of cycle slips, the cycle slip estimation in Eq. ( 4) can be used to calculate the integer solution that satisfies Eq. (7), it is the cycle slip value.

α coefficient optimization
In order to eliminate the influence of geometric distance and weaken the influence of ionospheric residuals, multiple sets of coefficients α can be found within a certain range, meeting the following conditions: To prevent error amplification, select α ∈ [−1010] as the range for combination, and through comparison, obtain the relevant information of the combination coefficients of GPS and BDS triple-frequency phase observations, as shown in Tables 1 and 2. The table shows that the combination coefficient can eliminate the distance term and reduce the influence of the ionosphere by up to 80% compared to its original effect I 1 .Based on the selected coefficients, cycle slip detection Eqs. ( 9) and ( 10) can be established, where j in �∅ j α represents the combination coefficient number.

Moving window global search method
The cycle slip will cause an outlier in �∅ α , and the location of the outlier is the epoch of the suspected cycle slip.By constraining the value range of N i , the integer solution satisfying Eqs. ( 9) and (10) can be calculated to obtain the cycle slip value of GPS and BDS.Therefore, it is important to determine the location of the cycle slip in the sequence �∅ α and the value range of N i .
Figure 1 illustrates the process of identifying outliers in a sequence X(n) of length n.The first step is to set a window of width w, starting from the first element of the sequence.This window consists of the first w elements of the sequence and is used to determine if the w + 1 element is an outlier.The window is then moved one element back, and the next w elements (2 ~ w + 1) are used to determine if the w + 2 element is an outlier, and so on.This method involves a global search, which involves locating outliers within the moving window, setting constraints on parameter values, and detecting and repairing triple-frequency cycle slips.( 7) www.nature.com/scientificreports/

Cycle slip detection and repair process
The algorithm flowchart is shown in Fig. 2, and the specific steps are as follows: (1) Based on the selected coefficients α , construct �∅  11) and ( 13) below is the cycle slip.Detect and calculate t epoch N i according to Eqs. ( 11) and ( 13); Simultaneously repair �N i (t) and �∅ j α (t) in the original sequence according to Eqs. ( 12) and ( 14); recalculate the standard deviations

Data source and experimental description
To validate the GPS/BDS triple-frequency cycle slip detection and repair method described in this article, we utilized data from two Hong Kong CORS stations HKKT and HKMW, which were observed on the 1st and 54th days of 2023, respectively.The HKKT station was observed for 5 h from 0:00 to 5:00, while the HKMW station was observed for 8 h from 5:00 to 13:00 with a sampling rate of 1 Hz.Table 3 displays the essential information of the satellites used in the experiment, and the data can be downloaded from ftp:// ftp.geode tic.gov.hk/.The moving window width was set to w = 30 epochs, which is equivalent to 30 s.The experiment was divided into five groups, as shown in Table 4.During data processing, P-code was used to replace pseudo-range observations in GPS, and the average of triple-frequency pseudo-range was used to replace pseudo-range observations in BDS.

Experiment 1
The GPS and BDS data utilize different combination coefficients.Specifically, GPS uses [1 − 6 5]  and N i before and after cycle slips were added to the G10, C04, C12, and C16.The cycle slip detection method is compared with the direct rounding method in Table 5.
Figures 3, 4, 5 and 6 demonstrate that cycle slips have a significant impact on the amplitude of �∅ j α and N i .The experiment achieved a detection success rate of over 94.0%, but there were some cases of false and non-detections.The analysis shows that the BDS system has more false detections than GPS due to the use of pseudo-range for cycle slip estimation calculation, while GPS uses P-code.Furthermore, it is important to note (13) that small cycle slips may be obscured by observation noise.The experimental results show that the average success rate is 1.3% higher than the direct rounding method used.

Experiment 2
The experiment used the same original observation data and combination coefficients as Experiment 1 for each satellite, but with added cycle slip sizes of [−10 −4] and [4 10] cycles.Figures 7, 8, 9 and 10 compared the G10, C04, C12 and C16 before and after adding cycle slips.Table 6 shows the detection effect of this experiment.From Figs. 7, 8, 9 and 10, it can be seen that before and after adding cycle slips, the amplitudes of �∅ j α and N i significantly increase.However, GPS has a higher success rate than BDS in terms of detection effectiveness,  with BDS experiencing false detections, non-detections, and multiple detections.Despite this, the overall detection success rate is over 97.0%.The experimental results demonstrate that this method has an average success rate 1.6% higher than the direct rounding method.

Experiment 3
The experiment used the same original observations and combination coefficients for each satellite as in Experiment 1, but with the addition of cycle slip sizes of [−4 −2] and [2 4] cycles.Figures 11, 12, 13 and 14 compare the G10, C04, C12 and C16 satellites before and after adding the cycle slips.The detection results are presented in Table 7.  α with a relatively large amplitude, N i shows a slight increase in amplitude.Undetected situations were observed for multiple system satellites in terms of detection effectiveness.It is analysed that the added cycle slip size is submerged by noise, resulting in false and multiple detections by BDS.The accuracy of the pseudo-range is significantly lower than that of the P-code.However, the overall success rate exceeds 96.0%.The experimental results demonstrate that this method has an average success rate 0.7% higher than that of the direct rounding method.

Experiment 4
To test the effectiveness of the method described in this article in detecting and repairing small cycle slips, cycle slip sizes of [−2 2] cycles were added in the experimental design.The observation data and combination coefficients used for each satellite in this experiment are identical to those used in Experiment 1. Figures 15, 16, 17  and 18 compare the G10, C04, C12, and C16 satellites before and after adding the cycle slips.The detection effect of this experiment is shown in Table 8.
From Figs. 15, 16, 17 and 18, it can be seen that before and after adding cycle slips, �∅ j α changes significantly and the curve has sharp parts, the change in N i is gentle, with almost the same shape before and after.From the perspective of detection effectiveness, the method of using P-code to constrain GPS cycle estimation for the added small cycle slips is significantly higher than that of BDS using pseudo-range average constraint.In terms of the number of undetected cycles, BDS satellites are significantly higher than GPS satellites.The overall success rate is over 92.0%.From the experimental results, it can be seen that the average success rate of this method is 2.5% higher than that of the direct rounding method.9. Figures 19, 20 and 21 show that the same size cycle slip method with fixed epochs did not result in any false or multiple detections.However, there were still some cases of undetected detections.Experiment 5 demonstrates that the detection success rate of this method is consistent across different types of BDS satellites, with a success rate of over 95.0%.The experimental results demonstrate that this method has an average success rate 2.0% higher than the direct rounding method.

Analysis
(1) The results of Experiments 1 to 5 demonstrate that the moving window global search method is effective in detecting and repairing cycle slips within 10 cycles, 4-10 cycles, 2-4 cycles, and within 2 cycles, with average detection success rates of 97.2%, 98.8%, 97.5%, 95.0% and 97.0%, respectively.No issues with multiple values were observed.(2) The average detection success rate across all five experiments was 97.1%, which is 1.5% higher than the success rate of the direct rounding method at 95.6%.(3) The success rate of GPS cycle slip detection using P-code constraints is higher than that of BDS cycle slip detection using pseudo range constraints, with an average increase of about 1.8%.( 4) Experiments 1 to 4 showed that using P-code instead of pseudo-range for GPS cycle slip detection resulted in an average success rate of 99.2%, which is 0.7% higher than the direct rounding method of 98.5%.The detection effect is equivalent.Additionally, using the average of triple-frequency pseudo-range instead of BDS for cycle slip detection resulted in an average success rate of 96.6%, which is 1.8%  higher than the direct rounding method of 94.8%.(5) In general, the moving window global search method has a higher success rate than the direct rounding method.

Conclusion
Based on the moving window global search method, the following conclusions were drawn from the experimental analysis: (1) The GPS and BDS systems both utilize different combination coefficients to form the GF model can effectively reduce the impact of ionospheric interference.By implementing a moving window and constraining the search range using the repaired sequence standard deviation, cycle slips can be accurately detected and repaired.The success rate range from 92.0% to 100.0%, demonstrating its feasibility and effectiveness.(2) Applying the "3σ" criterion to constrain the search range in a moving window can effectively avoid the multi value problem in cycle slip detection process.(3) The moving window global search method is able to detect and repair both small cycle slips within 2 cycles and larger cycle slips exceeding 2 cycles.
The detection and repair of cycle slips in BDS with more than triple-frequency, dynamic modes, and different sampling conditions will be our future focus of research.

Figure 3 .
Figure 3.Comparison of �∅ j α and N i before and after adding cycle slips to G10 (the first column represents �∅ j α and the second column represents N i ; The first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 4 .
Figure 4. Comparison of �∅ j α and N i before and after adding cycle slips to C04 (the first column represents �∅ j α and the second column represents N i ; The first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 5 .
Figure 5.Comparison of �∅ j α and N i before and after adding cycle slips to C12 (the first column represents �∅ j α and the second column represents N i ; The first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 6 . 5 Figures 11 ,
Figure 6.Comparison of �∅ j α and N i before and after adding cycle slips to C16 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Table 5 .Figure 7 . 5 5
Figure 7.Comparison of �∅ j α and N i before and after adding cycle slips to G10 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 8 .
Figure 8.Comparison of �∅ j α and N i before and after adding cycle slips to C04 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 9 .
Figure 9.Comparison of �∅ j α and N i before and after adding cycle slips to C10 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 10 .
Figure 10.Comparison of �∅ j α and N i before and after adding cycle slips to C16 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 11 .
Figure 11.Comparison of �∅ j α and N i before and after adding cycle slips to G10 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 12 .
Figure 12.Comparison of �∅ j α and N i before and after adding cycle slips to C04 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 13 .
Figure 13.Comparison of �∅ j α and N i before and after adding cycle slips to C12 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 14 .
Figure 14.Comparison of �∅ j α and N i before and after adding cycle slips to C16 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 15 .
Figure 15.Comparison of �∅ j α and N i before and after adding cycle slips to G10 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 16 .Figure 17 .
Figure 16.Comparison of �∅ j α and N i before and after adding cycle slips to C04 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 18 .
Figure 18.Comparison of �∅ j α and N i before and after adding cycle slips to C16 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 19 .
Figure 19.Comparison of �∅ j α and N i before and after adding cycle slips to C04 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 20 .
Figure 20.Comparison of �∅ j α and N i before and after adding cycle slips to C12 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Figure 21 .
Figure 21.Comparison of �∅ j α and N i before and after adding cycle slips to C13 (the first column represents �∅ j α and the second column represents N i ; the first three rows and the last three rows represent before and after adding the cycle slips respectively).

Table 9 .
Detection effect and comparison of experiment 5.

Table 3 .
Basic information of experimental satellites.

Table 4 .
Basic information of experimental design.

Table 6 .
Detection effect and comparison of experiment 2.

Table 7 .
Detection effect and comparison of experiment 3.
PRN Total

Table 8 .
Detection effect and comparison in experiment 4.