A hierarchical combination algorithm for real-time cycle slip detection and repair in low satellite elevation and high ionospheric activity conditions

To enhance the accuracy and robustness of cycle slip detection and repair for triple-frequency data while minimizing the adverse effects of low satellite elevation and high ionospheric activity, a hierarchical combination algorithm for real-time cycle slip detection and repair is proposed. This algorithm begins by prioritizing the reduction of noise and ionospheric delay coefficients. It determines the optimal coefficients for the combination of observations from the BeiDou Navigation Satellite System’s (BDS) Extra-Wide Lane (EWL), Wide Lane (WL), and Narrow Lane (NL). Leveraging the longer wavelength characteristics of the EWL combination, it simultaneously conducts cycle slip detection on the EWL combination alongside the pseudorange combination. Following this, based on the detection outcomes from the EWL combination, cycle slip detection is carried out on the WL combination. Finally, using the detection findings from the WL combination, cycle slip detection is executed on the NL combination. Given the NL combination’s shorter wavelength and higher susceptibility to ionospheric delay, a dynamic ionospheric prediction model is applied to the NL combination to further mitigate the impact of ionospheric disturbances. After completing the cycle slip detection process, the results from the EWL, WL, and NL combinations are integrated and solved. Experimental results clearly demonstrate that, even in scenarios characterized by low satellite elevation and active ionospheric conditions, this algorithm consistently delivers outstanding detection performance for cycle slip instances, particularly for small cycle slip (less then two cycles). Remarkably, this performance is achieved without the need for intricate searches during cycle slip repair.

The BeiDou Navigation Satellite System (BDS) has injected fresh vigor into the advancement of Global Navigation Satellite Systems (GNSS).BDS boasts widespread applications in diverse domains such as aviation, maritime navigation, vehicle navigation, surveying, and geographic information systems [1][2][3][4] .Nevertheless, the practical use of GNSS receivers is often fraught with challenges due to the multitude of error sources during signal propagation.One such significant challenge is the cycle slip issue.A cycle slip denotes a sudden and discontinuous shift in the phase of satellite signals received by a GNSS receiver, stemming from various causes 5,6 .This phase disruption can result in inaccuracies in positioning, subsequently compromising navigation precision and dependability.To confront the cycle slip predicament, numerous researchers have put forth various cycle slip detection and repair methodologies.Common cycle slip detection approaches include those grounded in combination models, statistical features, and filtering techniques.
The cycle slip detection method based on combination models identifies cycle slips by skillfully amalgamating diverse models or multiple data sources to establish a comprehensive observation equation or parameter inversion model.This encompasses techniques such as the TurboEdit algorithm [7][8][9] , the pseudorange phase combination method [10][11][12] , and the ionospheric residual method [13][14][15] , among others.The TurboEdit algorithm merges the Melbourne-Wübbena (MW) combination with the Geometry-Free (GF) combination, effectively mitigating

BDS triple-frequency basic combination
The pseudorange and carrier observation equations for BDS are expressed as follows: where:n = 1, 2, 3 denotes the different frequencies of BDS.For IGSO/MEO satellites, they represent the L1X, L5X, and L6I respectively, and for MEO satellites, they represent the L2I, L7D, and L6I respectively;ρ denotes the geometric distance between the receiver antenna and the satellite;c denotes the speed of light;δT r and δT denote the receiver clock and the satellite clock, respectively;P and ϕ denote pseudorange observations and carrier phase observations, respectively;γ denotes the ionospheric delay coefficient;δI and trop denote the ionospheric and the tropospheric delays, respectively; denotes carrier wavelength; N(t 0 ) denotes ambiguity of whole cycles ; ε denotes observation noise.Combining three frequencies, the triple-frequency pseudorange and carrier combination can be represented respectively as: (1) where: f denotes the carrier frequency; (a, b, c) denote coefficients of the pseudorange combination; (i, j, k) denote coefficients of the carrier combination; f (i,j,k) = if 1 + jf 2 + kf 3 denotes the frequency of the carrier combination.
To reduce errors, such as the receiver and satellite clock biases, differences are computed using observations between adjacent epochs 17 .The differences in pseudorange observations between adjacent epochs of Eq. ( 3) are expressed as: where: denotes difference between adjacent epochs; γ P denotes the ionospheric delay coefficients of pseudor- ange combination; σ P denotes the noise of pseudorange combination.
The differences in carrier observations between adjacent epochs of Eq. ( 4) are expressed as where: γ l denotes the ionospheric delay coefficients of carrier combination; (i,j,k) denotes the wavelength of the carrier combination; N (i,j,k) (t 0 ) denotes the ambiguity of carrier combination.σ l denotes the noise of carrier combination.
The specific expressions for each of these quantities are as shown in Table 1: In order to find the optimal linear combination of carriers, it is essential to establish effective selection criteria.As per Eq. ( 6), the residual ionospheric delay and observation noise play a crucial role in determining phase ambiguities.Taking into consideration the different combination wavelengths, the following two criteria are adopted: (1)Minimize the ionospheric delay coefficients to reduce the impact of ionospheric variability; (2)Minimize the interference of observation noise on data as much as possible.These criteria can be expressed as follows: Based on the geometric analysis in reference 30 , simultaneously mitigating errors caused by observation noise and ionospheric delay presents a set of conflicting factors.Therefore, in practice, there is a need to strike a balance between these two factors by selecting suitable combination coefficients through a compromise.To limit the error magnitude, integer coefficients should be determined within the range of − 5 to 5. Assuming that the carrier observation noise for the three frequencies is independently and identically distributed with the same standard deviation 31,32 , and given that the carrier phase noise for different BDS frequency bands is the same(i.e.,ε ϕ 1 = ε ϕ 2 = ε ϕ 3 = 0.002m),the coefficients l (0,−1,1) , , l (1,0,−1) and l (1,0,0) are selected as the best combina- tion based on the selection criteria.The BDS carrier frequencies are L1X, L2I, L5X, L6I, and L7D, with IGSO/ MEO satellites using L1X, L5X, and L6I frequencies, and GEO satellites using L2I, L7D, and L6I frequencies.The Extra-Wide Lane (EWL), Wide Lane (WL), and Narrow Lane (NL) combinations parameters for the threefrequency BDS signals are as shown in Table 2.
(3) www.nature.com/scientificreports/According to Eq. ( 5), pseudorange combination is primarily affected by pseudorange observation noise.Similar to carrier observation noise 33 , assume that the pseudorange noise for BDS is the same (i.e.,ε P 1 = ε P 2 = ε P 3 = ε P = 0.3m ).Referring to Table 1, it is apparent that when a = b = c = 1/3 , the combi- nation noise is minimized.Therefore, selecting P (1/3,1/3,1/3) as the pseudorange combination coefficient.After determining the coefficients for carrier and pseudorange combination, a hierarchical combination model is established to facilitate cycle slip detection and repair.

EWL combination cycle slip detection model
The EWL combination retains the integer cycle slips of the combination ambiguity, and since the wavelength is large, it is not affected much by the residual errors.First, cycle slips of the EWL combination are detected and repaired.Combined with the pseudorange combination, the pseudorange and carrier phase combination is constructed according to Eqs. ( 5) and ( 6), which is expressed as: The ionospheric coefficients η 1 for GEO satellites and IGSO/MEO satellites in the EWL combination are 0.040 and 0.066, respectively.Correspondingly, the observation noise values σ 1 are 0.037 and 0.054, respectively.When the EWL combination observation exceeds 0.5, indicating that �N (0,−1,1) > 0.5 , it is considered a cycle slip.In such instances, the cycle slip value for the EWL combination is simply rounded to the nearest integer, denoted as � N(0,−1,1) = round[�N (0,−1,1) ] .Given the negligible magnitude of η 1 , it can effectively be disregarded when using a 30 s sampling interval 31 .Consequently, the success rate of cycle slip detection for the EWL combination is as follows:

WL combination cycle slip detection model
After the cycle slip of the EWL combination is determined, the repaired EWL combination is utilized to detect the observations of the WL combination: The ionospheric coefficients η 2 for GEO satellites and IGSO/MEO satellites in the WL combination are -0.352 and -0.431, respectively.Correspondingly, the observation noise values σ 2 are 0.068 and 0.040 respectively.When the WL combination observation exceeds 0.5, indicating that �N (1,0,−1) > 0.5 , it is considered a cycle slip.In such instances, the cycle slip value for the WL combination is simply rounded to the nearest integer, denoted as � N(1,0,−1) = round[�N (1,0,−1) ].Similarly to the EWL combination, since the η 2 value remains small, it can be ignored with a 30 s sampling interval for WL combination.Therefore, the success rate of cycle slip detection for the WL combination is:

NL combination cycle slip detection model
After the cycle slip of the WL combination is determined, the repaired WL combination is utilized to detect the observations of the NL: Vol.:(0123456789)The ionospheric coefficients η 3 for GEO satellites and IGSO/MEO satellites in the NL combination are − 11.941 and − 11.781, respectively.The observation noise values σ 3 are 0.058 and 0.070, respectively.It is evident that when the value of η 3 experiences a significant increase, direct integer rounding of NL's cycle slip values may not be suit- able during periods of heightened ionospheric activity.Polynomial regression models are commonly employed to capture temporal variations.In this context, a polynomial regression function incorporating time variables can be introduced to model the time series of NL combination values.This facilitates the prediction of epochto-epoch ionospheric delay within the NL combination, thereby mitigating the influence of ionospheric delay.The general form of a polynomial regression fit with a window size of n and a polynomial order of p is as follows: where: i = 1, 2, 3, • • • , m;x and y denote polynomial regression coefficients and NL combination observations, respectively.t denotes the time of the epoch.The choice of the polynomial fit order p in polynomial regression should be based on the variations of each NL combination value with respect to time within the data cycle.According to prior research: During periods of active ionospheric conditions, p is set to 2. And during stable ionospheric conditions, p is set to 1 34 .The decreasing correlation between observations over time has been confirmed through multiple tests.It has been found that a time interval of 5 to 15 min is sufficient to meet the requirements of the prediction window under normal conditions.However, when satellite elevation is relatively low, and noise effects are more pronounced, longer time intervals are needed for prediction.Therefore, in the calculation of ionospheric delay, a dynamic window of length m is employed, with the maximum window size set to m max = 30.The specific implementation strategy is as follows: (1) When the receiver initially receives satellite signals, especially when the satellite elevation is low, and effective data support is needed for subsequent NL combination observations detection, the window size should be as large as possible.Therefore, in this case, the window size is set to the maximum window count, which is m = m max ; (2) When the satellite elevation EL is greater than 30 • , the variations in multipath error and noise error tend to be relatively smooth.Therefore, in this situation, it is advisable to use a fixed window size for data processing, which is typically expressed as m = 0.5m max .(3) When EL is less than 30 • but greater than 15 • , where noise error has a significant impact, it is advisable to dynamically adjust the window size to better accommodate the changing elevation.In this case, the window size can be set as m = m max (1 − sin el).(4) When EL is less than 15 • , noise error significantly increases, and in such conditions, it is advisable to use the maximum window size.Additionally, when EL is less than 10 • , it becomes challenging to capture satellite signals, and the signal quality is very poor.Therefore, for EL less than 15 • greater than 10 • the window size can be set as m = m max .
In summary, the dynamic window size can be expressed as follows: After determining the order of polynomial fitting and the window, the following equation can be derived based on Eq. ( 16): Equation ( 18) can be simplified into matrix form as A NL X NL = Y NL .Subsequently, the error equation matrix can be expressed as V NL = AX NL − Y NL .Employing the least squares principle and the function-free extrema algorithm, we can eliminate V NL to obtain the polynomial fitting coefficients XNL , denoted as: After obtaining the polynomial fitting coefficients, the predicted ionospheric layer delay of the ( m + 1)th epoch can be expressed as: (15) At the same time, in order to ensure the accuracy of ionospheric prediction results, it should be ensured that there is no cycle slip in the data within the initial window.The judgment algorithm is as follows: (1) According to the Eq. ( 19), obtain the polynomial fitting coefficients T within the initial window, then the predicted ionospheric delay I NL within the initial window can be expressed as

Results and discussion
The data comes from the IGS WUH2 station TRIMBLE ALLOY GNSS station receiver.The collection time is September 04, 2022, and the sampling rate is 30 s.A total of 2880 epochs of BeiDou satellite observation data were selected for experimental analysis.The geomagnetic Kp index on the day of the experimental data was obtained from the Space Environment Prediction Center (http:// www.sepc.ac.cn) as shown in Fig. 2: According to the calculation, the average Kp index on September 4, 2022 is 4.75, exceeding 5 most of the time, which means a major geomagnetic storm occurred that day and the ionosphere was in a violent activity state.
The experiment is divided into two schemes, as follows: Scheme I: Verify with the original station data.The algorithm of hierarchical combination cycle slip detection and repair (hereinafter referred to as Algorithm 1) is directly used to detect and repair cycle slips in the original data, and compared with the three-frequency Geometry-Free and Ionosphere-Free(GFIF) combination algorithm (hereinafter referred to as Algorithm 2) detection and repair results in Reference 16 to verify the basic performance of Algorithm 1. Algorithm 2 incorporates MW combination, GIGF combination, and PIR combination, building upon the TurboEdit algorithm, which is currently one of the most widely used cycle slip detection algorithms.In order to further ensure the accuracy of the Algorithm 2 for cycle slip repair, the least-squares ambiguity decorrelation adjustment (LAMBDA) 17 is used to search for cycle slip candidates of the Algorithm 2. Because for a carrier data segment, the type of cycle slips occurred is singular.To further verify the applicability of Algorithm 1 to various types of cycle slips, obtain the carrier data without cycle slips through Scheme 1 for the next experiment.
Scheme II: Taking the carrier data without cycle slips from Scheme I as the basis, first add different types of cycle slip combinations from the ( m max + 1)th epoch of the carrier data at intervals of 5 epochs.The first 8 cycle slip combinations are insensitive small cycle slips of EWL, WL or NL combinations, which are (1,0,0), (0,1,0), (0,0,1), (0, 2, 2), (3, 0, 3), (5, 4, 0), (5, 5, 5) and (13, 10, 0).Other cycle slip combinations are random cycle slip combinations of 0-9 cycles randomly generated for each frequency (not all zeros for the 3 frequencies).Algorithm 1 and Algorithm 2 are used to detect and repair cycle slips in the experimental data, with the purpose of fully verifying the correctness, effectiveness and applicability of Algorithm 1 in detecting and repairing different cycle slip combinations (especially insensitive cycle slips).

Results and analysis of scheme 1
Figure 3 shows the detection results of Algorithm 1 for satellites C19, C38, and C59.It can be seen from Fig. 3 that the EWL combination is not affected by the satellite elevation and ionosphere, and the combination observations are within the detection threshold range without significant fluctuations.The WL combination is not affected by the ionosphere.Due to the influence of satellite elevation, the fluctuation amplitude of the combination observation increases at low elevation, but the fluctuation range is still within the detection threshold.The NL combination is affected by both elevation and ionosphere, so the fluctuation range of the combination observations changes dramatically, with a considerable portion exceeding the detection threshold.However, the NL combination improved by the dynamic ionospheric prediction model is not affected by elevation and ionosphere.Its combination observations are within the detection threshold range without significant fluctuations.
Figure 4 shows the detection results of Algorithm 2 for satellites C19, C38, and C59.It can be seen from Fig. 4 that the GFIF combination is not affected by satellite elevation and ionosphere.The combination observations are within the detection threshold range without significant fluctuations.The MW combination is not affected by the ionosphere.Due to the severe influence of satellite elevation, the fluctuation range of the combined values increases significantly at low elevation, with some exceeding the threshold.The PIR combination has an obviously larger fluctuation range in the combination observations due to the dual influence of elevation and ionosphere, and some of the combination observations exceed the detection threshold.
The cycle slip detection results of Algorithm 1 are shown in Table 3. Algorithm 1 did not detect cycle slips on satellites C19 and C59, and detected one group of cycle slips on satellite C38 which can be corrected successfully, without false detection.The unimproved algorithm (without using ionospheric prediction model), a total of 183 sets of cycle slips were detected, and none of them were successfully repaired.Therefore, the ionospheric Therefore, compared with Algorithm 2, Algorithm 1 has higher accuracy and better applicability in cycle slip detection and repair.After the measured data was repaired by Algorithm 1, the changes of all combination observations were within the detection threshold.The cycle slip detection figures after repair are not shown due to limited space.
Figure 6 shows the detection results of Algorithm 2 for simulated cycle slips added to satellites C19, C38 and C59.It can be seen from Fig. 6 that in Algorithm 2, GFIF combination is insensitive to the combined cycle slips (1,0,0), (0,2,2) and (5,5,5), but can detect other combined cycle slips.The MW combination is insensitive to the combined cycle slips (0,1,0), (3,0,3) and (5,5,5).At the same time, due to the influence of satellite elevation, it still cannot accurately detect the remaining cycle slip combinations.The PIR combination is insensitive to the combined cycle slip (13,10,0), but still cannot detect all the remaining cycle slip combinations due to the influence of ionospheric activity.Therefore, by combining GFIF, MW and PIR, Algorithm 2 can also detect the added simulated cycle slips, but there will be false detections.
The repair results of simulated cycle slips by Algorithm 1 are shown in Table 5. Statistically, Algorithm 1 can correctly repair all the 1123 groups of simulated cycle slips added to the 3 satellites, especially for insensitive small cycle slip combinations, with a 100% correction accuracy.
The repair results of simulated cycle slips by Algorithm2 are shown in Table 6.Among the 1123 groups of simulated cycle slips added to the 3 satellites, Algorithm 2 correctly repaired 1113 groups and incorrectly repaired

Figure 1 .
Figure 1.Flowchart of the cycle slip detection and repair algorithm.

Figure 5 .
Figure 5. Simulated cycle slip detection results of Algorithm 1.

Figure 6 .
Figure 6.Simulated cycle slip detection results of Algorithm 2.

Table 5 .
Simulated cycle slip repair results of Algorithm 1.

Table 6 .
Simulated cycle slip repair results of Algorithm 2.