Cycle Slip Detection during High Ionospheric Activities Based on Combined Triple-Frequency GNSS Signals

Abstract

The current cycle slip detection methods of Global Navigation Satellite System (GNSS) were mostly proposed on the basis of assuming the ionospheric delay varying smoothly over time. However, these methods can be invalid during active ionospheric periods, e.g., high Kp index value and scintillations, due to the significant increase of the ionospheric delay. In order to detect cycle slips during high ionospheric activities successfully, this paper proposes a method based on two modified Hatch–Melbourne–Wübbena combinations. The measurement noise in the Hatch–Melbourne–Wübbena combination is minimized by employing the optimally selected combined signals, while the ionospheric delay is detrended using a smoothing technique. The difference between the time-differenced ambiguity of the combined signal and this estimated ionospheric trend is adopted as the detection value, which can be free from ionospheric effect and hold the high precision of the combined signal. Five threshold determination methods are proposed and compared to decide the cycle slip from the magnitude aspect. This proposed method is tested with triple-frequency Global Navigation Satellite System observations collected under high ionospheric activities. Results show that the proposed method can correctly detect and fix cycle slips under disturbed ionosphere.

1. Introduction

High precision Global Navigation Satellite System (GNSS) positioning requires the availability of high quality carrier phase observations. However, carrier phase signals are easily affected by high ionospheric activities, e.g., ionospheric scintillation, which occurs because of plasma density irregularities and causes rapid phase and amplitude variations of GNSS and other radio signals [1,2,3]. This phenomenon leads to the degradation of GNSS signals and on occasion can also lead to complete loss of locks. The existing studies already show a strong correlation of scintillation and cycle slip occurrence [4,5,6,7,8]. These cycle slips need to be taken into careful consideration in order to maintain the precision and reliability of position.

When high ionospheric activity occurs, the sudden changes of the noise level and ionospheric delay can be problematic to the methods developed based on the legacy signals, e.g., high-order phase differencing [9], Doppler integration [10], and some of the latest proposed methods, such as Huang et al. [11], de Lacy [12] and Zhao et al. [13], as these methods are raised on the basis of the assumption that the ionospheric delay varies smoothly over time. More literature on cycle slip detection is available in Section 2. Although the method proposed by Zhao et al. [13] cannot provide a robust detection performance during high ionospheric activities, one of the proposed combined signals (4, 0, −5) in the modified Hatch–Melbourne–Wübbena (HMW) combinations can maintain a high precise detection value even when it is biased by the ionosphere, revealing the potential to expand the application of the modified HMW combination in the scenario of high ionospheric activities. Inspired by this, this paper proposes newly constrained conditions to select the optimally combined signals of the modified HMW combinations to detect cycle slips in the observations collected by geodetic receivers during high ionospheric activities, as shown in Section 3. The measurement noise is minimized in these combined signals, but the ionospheric effect is amplified. In double-differenced mode, the ionospheric bias can be eliminated using the method proposed by Sieradzki and Paziewski [14,15], while the cycle slip is usually checked in a time-differenced mode. This paper employs a locally weighted scatter-plot smoothing technique to extract the ionospheric trend in the detection value, which will work with the Generalized AutoRegressive Conditional Heteroscedasticity (GARCH) with a partly fixed threshold determination method to detect all the cycle slips, as expressed in Section 4. The performance of this proposed method is extensively evaluated and tested with the data collected with different ionospheric conditions on account of its ability to detect existing and artificially added cycle slips in Section 5 and Section 6, respectively. Finally, in the last section, conclusions will be drawn.

2. Literature Review

Some efforts have been conducted to account for the effect of high ionospheric activities in cycle slip detection and correction using single-/dual-frequency signals. The impact of traveling ionosphere disturbances on cycle slip detection was analyzed by Tang et al. [16] using dual-frequency geometry-free combinations. Ji et al. [17] tried to expand the high-order differencing method to the scenario of ionospheric scintillation by introducing a non-geometry-free but ionosphere-free dual-frequency combination, but their method may not be very reliable according to their detection results. Liu [18] detected the cycle slips by calculating the difference between ionospheric total electron contents rate of the current epoch and that predicted by previous epochs. Although this method was tested to be able to detect any cycle slip, even under high ionospheric activities, GPS data with a sampling rate of 1 Hz or higher should be applied in order to keep the sensitivity to small cycle slips. Cai et al. [19] expanded the method of Liu by using a forward and backward moving window averaging algorithm, but this method may not successfully detect cycle slips in conditions with significant measurement noise. Zhang and Li [20], and Banville and Langley [21] tried to repair cycle slips to reduce the re-initializations in precise point positioning (PPP) using dual-frequency signals, with the assistance of external information, e.g., precise orbit and clock corrections.

With the development of GNSS, the availability of triple-frequency signal provides new opportunities for cycle slip detection under high ionospheric activities by introducing additional combined signals with longer wavelengths, weaker ionospheric delays and smaller measurement noise [22,23]. Zhao et al. [24] tried to detect cycle slips under high ionospheric activities by predicting the ionospheric delay based on the assumption that ionospheric delay remains unchanged within a period. This assumption cannot be valid under high ionospheric activities, e.g., strong ionospheric scintillation, or for the observations with low sampling rate, e.g., 30 s. Besides that, the method to eliminate the ionospheric bias may lead to the measured value being inaccurate and thus losing sensitivity to cycle slips with one cycle magnitude. A similar ionospheric prediction method was also applied by Yao et al. [25], in which several previous carrier-phase observations without cycle slips are required to correct the ionospheric bias. This prediction method can only effectively predict the time-differenced ionospheric delay when the data epoch interval is no larger than 5 s [26]. When the data epoch interval is larger than 5 s, second-order ionospheric variations estimated from the previous 10 epochs should be applied by a linear fitting model [26]. Liu et al. [27] also applied the second-order ionospheric variation to predict the ionospheric delay but using less consecutive epochs. The ionospheric bias mainly affects the determination of the threshold; however, all of the previous studies did not pay attention to the determination of the threshold. Therefore, this paper will discuss the method to decide the optimal threshold for the detection value.

3. Selection of Optimally Combined Signals

Zhao et al. [13] proposed the following modified HMW combination for cycle slip detection,

$$\Delta N_{(i,j,k)}=\frac{1}{{\lambda }_{(i,j,k)}}(a\Delta P_1+b\Delta P_2+c\Delta P_3+d\Delta {\check{\Phi }}_{(0,1,-1)}-\Delta {\Phi }_{(i,j,k)}),$$

(1)

where Δ denotes the time differenced operation between two consecutive epochs; N is the ambiguity; the subscripts indicate the frequencies; $(i,j,k)$ indicates the combined signal with the coefficients i, j and k for the frequencies [23]; $\lambda$ is the wavelength; P and $\Phi$ represent the code measurement and the phase measurement in units of length respectively; ${\check{\Phi }}_{(0,1,-1)}$ is the combined phase measurement $(0,1,-1)$ where cycle slips have been corrected; a, b, c and d are the coefficients. Compared to the traditional HMW combination, this modified HMW combination can have a better manipulation of the ionospheric bias and measurement noise [13]. The difficulty of this modified HMW combination is the determination of the optimally combined signals, which should be able to provide the most precise time-differenced ambiguities.

In order to select the optimally combined signals for detecting cycle slips during high ionospheric activities, we propose the following three constrained conditions. The first is the geometry-free condition, which means the sum of the coefficients a, b, c and d should be equal to 1. The second is called the linearly independent condition, which means the selected signals should be linearly independent from each other. The third constrained condition is to minimize the effect of the measurement noise on $\Delta N_{(i,j,k)}$. Based on the empirical assumptions for the phase and code measurement noise given in Zhao et al. [13], the effect of the measurement noise on $\Delta N_{(i,j,k)}$ can be calculated as follows:

$${\sigma }_{\Delta N_{(i,j,k)}}=\frac{1}{{\lambda }_{(i,j,k)}}\sqrt{a^2{\sigma }_{\Delta {\epsilon }_1}^2+b^2{\sigma }_{\Delta {\epsilon }_2}^2+c^2{\sigma }_{\Delta {\epsilon }_3}^2+2(d^2{\mu }_{(0,1,-1)}^2+{\mu }_{(i,j,k)}^2){\sigma }_{\Phi }^2},$$

(2)

where $\sigma$ denotes the standard deviation of the measurement noise; $\epsilon$ denotes the code measurement noise. ${\mu }^2$ is the phase noise factor. Determining the optimally combined signals in Equation (1) with the previous three constraints is an equality constrained optimization problem, which can be solved using the method proposed by Zhao et al. [13].

After considering all the combined signals with their i, j and k taking the values from −50 to 50, two combined signals least affected by measurement noise are selected for each satellite system, i.e., the combined signals (−3, 1, 3) and (4, −5, 0) for GPS and Galileo, the combined signals (−3, 4, 0) and (−4, 0, 5) for BeiDou Navigation Satellite System (BDS).

Table 1. Minimal noise level of the detection theoretically and ISF. The unit of the noise is in cycle.

Combined Signals	${\mathit{\sigma }}_{\Delta {\mathit{N}}_{(\mathit{i},\mathit{j},\mathit{k})}}$			ISF
	GPS	Galileo	BDS	GPS	Galileo	BDS
(−3, 1, 3) or (−3, 4, 0)	0.057	0.057	0.072	119.81	359.32	24.91
(4, −5, 0) or (−4, 0, 5)	0.096	0.130	0.087	−21.55	−13.30	35.78

shows the minimal effect of measurement noise on each combined signal and the ionospheric scale factor (ISF), while the coefficients are listed in

Table 2. Coefficients of the selected signals. The unit of the interval is in seconds.

SAT Type	Interval	a	b	c	d
GPS	1	0.0304	0.0475	0.0621	0.8600
	2	0.0271	0.0469	0.0540	0.8720
	3	0.0243	0.0461	0.0474	0.8822
	4	0.0219	0.0454	0.0419	0.8908
	5	0.0198	0.0445	0.0373	0.8983
	6	0.0180	0.0437	0.0335	0.9048
	7	0.0164	0.0429	0.0302	0.9106
	8	0.0151	0.0420	0.0273	0.9156
	9	0.0138	0.0412	0.0248	0.9201
	≥10	0.0128	0.0403	0.0227	0.9242
Galileo	1	0.1536	0.2945	0.3456	0.2062
	2	0.1554	0.2635	0.3207	0.2604
	3	0.1534	0.2367	0.2959	0.3139
	4	0.1490	0.2133	0.2722	0.3655
	5	0.1430	0.1927	0.2500	0.4143
	6	0.1361	0.1745	0.2294	0.4600
	7	0.1289	0.1584	0.2104	0.5023
	8	0.1216	0.1440	0.1932	0.5412
	9	0.1143	0.1313	0.1775	0.5768
	≥10	0.1074	0.1200	0.1634	0.6092
BDS GEO	1∼30	0.0344	0.0611	0.2208	0.6837
BDS IGSO	1	0.0461	0.0720	0.2453	0.6367
	2	0.0432	0.0698	0.2249	0.6620
	3	0.0406	0.0676	0.2068	0.6850
	4	0.0381	0.0653	0.1907	0.7059
	5	0.0358	0.0630	0.1762	0.7250
	6	0.0337	0.0607	0.1632	0.7424
	7	0.0317	0.0585	0.1515	0.7583
	8	0.0299	0.0561	0.1409	0.7731
	9	0.0282	0.0542	0.1313	0.7864
	10	0.0266	0.0523	0.1224	0.7988
	11	0.0251	0.0503	0.1146	0.8100
	12	0.0238	0.0483	0.1077	0.8202
	13	0.0226	0.0465	0.1012	0.8298
	14	0.0214	0.0447	0.0952	0.8387
	15	0.0203	0.0431	0.0897	0.8469
	16	0.0193	0.0415	0.0847	0.8545
	17	0.0184	0.0400	0.0800	0.8616
	18	0.0175	0.0385	0.0758	0.8682
	19	0.0167	0.0371	0.0717	0.8745
	≥20	0.0159	0.0358	0.0680	0.8802
BDS MEO	1	0.0284	0.0726	0.2018	0.6972
	2	0.0250	0.0614	0.1560	0.7576
	3	0.0219	0.0522	0.1236	0.8024
	4	0.0191	0.0445	0.1000	0.8364
	5	0.0168	0.0382	0.0824	0.8626
	6	0.0148	0.0331	0.0687	0.8833
	7	0.0131	0.0289	0.0583	0.8996
	8	0.0117	0.0254	0.0500	0.9129
	9	0.0105	0.0226	0.0433	0.9237
	≥10	0.0095	0.0201	0.0378	0.9326

. For a given observation interval, the same coefficient can be applied to both of the two combined signals for each satellite system. With the selected combined signals, the standard deviation (STD) of the detection value for nearly all the satellite systems can achieve about 0.1 cycle, which can guarantee a 1% missed detection rate at most and a 99.994% success rate (SR) when the ionospheric bias are totally eliminated based on the discussion of Zhao et al. [13]. However, from Table 1, all the selected combined signals magnify the ionospheric bias by tens or even hundreds of times compared to the effect on L1, E1 or B1 of GPS, Galileo and BDS, respectively. Therefore, the SR of detecting cycle slips in real observations mainly depends on the precision of the estimated ionospheric trend in the combined signals.

GEostationary Orbit (GEO); Inclined GeoSynchronous Orbit (IGSO); Medium Earth Orbit (MEO).

4. Using the Proposed Combined Signals to Detect Cycle Slips

This section introduces the process of cycle slip detection using the combined signals proposed in Section 3 to detect cycle slips. The first detection process employs the extra-wide-lane (EWL) combined signal, i.e., (0, 1, −1) for GPS and Galileo, (0, −1, 1) for BDS, to detect cycle slips and provide a cycle slip free combined signal, which will work with the proposed combined signal to form the modified HMW combinations. Besides eliminating the ionospheric bias, the detection value provided by this EWL combined signal can still be precise and valid to detect cycle slips. The second and the third detection processes employ the time-differenced ambiguities of the proposed combined signals which can be calculated based on Equation (1). Different from Zhao et al. [13], these time-differenced ambiguities cannot be used as the detection value directly, as they are highly biased by the disturbed ionosphere. The robust locally weight scatter plot smoothing (RLOWESS) technique is adopted to estimate the smoothed ionospheric trend in the time-differenced ambiguities. The difference between the time-differenced ambiguities and the smoothed value can nearly eliminate the effect of ionospheric biases and hold the high precision character of the proposed combined signals, thus this difference is used as the detection value. Once the magnitude of the cycle slips detected by the combined signals are known, those on the original signals can be determined by solving the equation set, as shown in [13].

4.1. Elimination of Outliers

Outliers and large cycle slips in the time-differenced ambiguities ($\Delta N$) should be removed, as they can affect the estimation of the ionospheric trend. The detection value in this process is the first-order time-differenced ambiguities ($\Delta \dot{N}$), which can be calculated as follows:

$$\Delta \dot{N}\left(k\right)=\Delta N\left(k\right)-\Delta N(k-1),$$

(3)

where k denotes epoch k, while the subscript (i, j, k) indicting the combined signals is omitted for the sake of simplicity. The reason for adopting $\Delta \dot{N}$ as the detection value is that the ionospheric residuals amplified by the combined signals can affect the detection of the outliers using the time-differenced ambiguities directly [27]. Although $\Delta \dot{N}$ can amplify the measurement noise, it is still precise enough to detect the outliers.

The threshold of how large an $\Delta \dot{N}$ can be recognized as an outlier is selected as four times of the STD [18]. The STD of $\Delta \dot{N}$ at epoch (k) can be calculated as below:

$${\sigma }_{\Delta \dot{N}}\left(k\right)=\sqrt{E[{(\Delta \dot{N}\left(k\right))}^2]-{(E[\Delta \dot{N}\left(k\right)])}^2},$$

(4)

where $\sigma$ denotes the STD operation here; $E[\Delta \dot{N}\left(k\right)]$ and $E[{(\Delta \dot{N}\left(k\right))}^2]$ are the mean value and the mean squared value of $\Delta \dot{N}\left(k\right)$ respectively, and can be recursively calculated as follows:

$$E[\Delta \dot{N}\left(k\right)]=E[\Delta \dot{N}(k-1)]+\frac{1}{k}\{\Delta \dot{N}\left(k\right)-E[\Delta \dot{N}(k-1)]\},$$

(5)

$$E\left[\begin{array}{c}{(\Delta \dot{N}\left(k\right))}^2\end{array}\right]=E\left[\begin{array}{c}{(\Delta \dot{N}(k-1))}^2\end{array}\right]+\frac{1}{k}\{{(\Delta \dot{N}\left(k\right))}^2-E[\begin{array}{c}{(\Delta \dot{N}(k-1))}^2\end{array}]\}.$$

(6)

The calculation of the STD in Equation (4) is the same as the method adopted by Blewitt [28] and Liu [18]. However, their methods adopt the observations of all the previous epochs to calculate STD, while the latest 80 epochs are just adopted to maintain the sensibility to the sudden fluctuation of $\Delta \dot{N}$, which is caused by the fast varying ionosphere. This proposed method to determine the threshold is called a local backward method in this paper.

Figure 1 shows the process of eliminating the outliers in the time-differenced ambiguities epoch by epoch. If the absolution of $\Delta \dot{N}\left(k\right)$ is larger than the threshold, it is regarded as an outlier. Then, $\Delta \dot{N}\left(k\right)$ is set to zero, which is the expectation of $\Delta \dot{N}$, in order to eliminate the effect of the outliers on the estimation of $\Delta N$ in the following epochs. The outliers are removed by setting $\Delta N\left(k\right)$ as its estimation, which can be calculated based on the measurements of the previous epochs, as all the epochs before the current epoch (k) are free from outliers. If any, they have been repaired by our epoch-by-epoch based method. Therefore, the estimation of $\Delta N\left(k\right)$ can be estimated as follows:

$$\Delta \tilde{N}\left(k\right)=\Delta N(k-1)+\Delta \dot{N}(k-1),$$

(7)

where $\Delta \tilde{N}\left(k\right)$ is the estimated value of $\Delta N\left(k\right)$; $\Delta \dot{N}(k-1)$ is estimated using measurements of previous 30 epochs in this study to smooth the noise by averaging.

4.2. Estimation of the Smoothed Ionospheric Bias

The aim of estimating the smoothed ionospheric bias is to provide the ionosphere-free cycle slip detection value based on $\Delta N$, where the outliers have been eliminated. Usually, $\Delta N$ is assumed to be a zero mean normal distribution; however, this is not the case under a fast varying ionosphere, which makes $\Delta N$ normally distribute along a trend fluctuating around zero. A non-parametric smoothing strategy called LOWESS [29] is applied to estimate the ionospheric trend of $\Delta N$. Non-parametric means the type of distribution needs not to be assumed in advance. A non-parametric strategy provides a flexible approach to representing data without the danger of leading to fitting a smooth curve that misrepresents the data. Thus, LOWESS can be used to create a smooth line to fit to $\Delta N$ in order to foresee the ionospheric variation trend.

Using LOWESS to obtain the ionospheric trend is introduced as follows: LOWESS smooths data by replacing each measurement with the weighted regression value determined by neighboring data points within the span. The regression weight ${\omega }_i$ is defined by the tri-cube function, shown as follows:

$${\omega }_i={(1-|\frac{x-x_i}{d\left(x\right)}{|}^3)}^3,$$

(8)

where x is the measurement to be smoothed; $x_i$ is the neighboring data of x within the span, and $d\left(x\right)$ is the distance along the abscissa from x to $x_i$. The measurement to be smoothed has the largest weight, so it affects the most on the fit. The weight is decreased with the increase of the distance. The smoothed value is given by performing a weighted linear least-squares regression with a first degree polynomial.

Although continuous and big outliers have been removed from $\Delta N$ in Section 4.1, the remaining outliers can still distort the smoothed values and affect the behavior of the neighboring measurements. In order to smooth the data without the influence of small outliers, a robust procedure is conducted as follows based on the previous LOWESS. The differences between the response measurement and the fit is also called the residuals ($r_i$). The median absolute deviation of the residuals ($MAD$) is calculated as follows:

$$MAD=\mathrm{median}\left(\right|r\left|\right).$$

(9)

$MAD$ measures the statistical dispersion of the residuals. According to the relation between the residual of each measurement and $MAD$, the robust regression weight (${\omega }_{ri}$) can be given by the following bi-square function:

$${\omega }_{ri}=\left\{\begin{array}{cc}{(1-{(r_i/6MAD)}^2)}^2,\hfill & |r_i|<6MAD,\hfill \\ 0,\hfill & |r_i|\ge 6MAD,\hfill \end{array}\right.$$

(10)

where the subscript i means the ith data point. The outliers are excluded from the smoothing process by setting the weights to zero for those measurements with the residuals greater than $6MAD$, while the smoothed values only reflect the behavior of those measurements with small residuals compared to $6MAD$, as the weights of these measurements are close to 1. Then, $\Delta N$ is smoothed again using the combined weight, which is the local regression weight (${\omega }_i$) multiplied with the robust regression weight (${\omega }_{ri}$). This robust procedure needs to be repeated for several times to get the optimal fit. The residuals between $\Delta N$ and the fit is adopted as the cycle slip detection value.

4.3. Determination of Detection Threshold

The last process of cycle slip detection is to set a proper threshold for the detection value. In cycle slip detection, the detection threshold is the magnitude that must be exceeded for the detection value to be manifested as a cycle slip. An optimal threshold should always be sensitive to small cycle slips. As the observations are post-processed in this paper, the threshold can be decided by either previous epochs or following epochs. In Section 4.1, the outliers are detected by a proposed threshold determination method, namely the local backward method, which is also valid to be used to detect cycle slips. Besides this method, the other four methods are summarized and proposed as follows.

The first method is called the sample variation method [30], which is to determine the threshold by four times the standard deviation of all the detection values during the observation period. This threshold determination method is simple for calculation and can guarantee a significance level of 0.006%. However, in the case of highly changing ionosphere, the changes of the code noise and the ionospheric effect can abruptly increase the standard deviation of the detection value. In such a case, the threshold determined by the standard deviation of the detection value may lead to the insensitiveness to small cycle slips or a great number of false alarms in the period where ionospheric scintillation occurs.

Another method to decide the threshold is GARCH process, which is widely used in modeling financial time series. Due to its ability to exhibit time-varying volatility clustering, GARCH can be used to estimate the volatility of the detection value, which is assumed to be conditional distributed and satisfies the following condition:

$$\Delta N_t|{\psi }_{t-1}\sim N(0,{\sigma }_t^2),$$

(11)

where ${\psi }_{t-1}$ denotes the conditions, including both the past conditional variances ${\sigma }_1^2$, ${\sigma }_2^2$, $···$, ${\sigma }_{t-1}^2$ and the past innovations $\Delta N_1$, $\Delta N_2$, $···$, $\Delta N_{t-1}$. The conditional variance is defined as follows in the GARCH model:

$${\sigma }_t^2=\sqrt{\varpi +\sum _{i=1}^q{\alpha }_i\Delta N_{t-i}^2+\sum _{i=1}^p{\beta }_i{\sigma }_{t-i}^2},$$

(12)

where p is the order of the GARCH terms ${\sigma }^2$ and q is the order of the AutoRegressive Conditional Heteroskedasticity (ARCH) terms $\Delta N^2$ in the GARCH(p, q) model. $\varpi$, $\alpha$ and $\beta$ are unknown parameters to be estimated. Equation (12) shows how the threshold estimated by GARCH can present the volatility of the detection value. If the former detection value of ARCH legs q has an unusually large absolute value, then ${\sigma }_t^2$ is larger than usual. The former conditional variances of GARCH legs p makes the estimated volatility more persistent.

The process to apply GARCH to estimate the threshold is listed as follows:

• Create a GARCH model. This is to determine the length of the GARCH legs and ARCH legs. The detailed process on how to establish the optimal length can refer to [31,32]. In most cases, the most widely used model GARCH(1, 1) can satisfy the requirement of conditional variance modeling for the detection series [30] because they have a nonzero conditional mean offset and exhibit volatility clustering.
• Use the maximum likelihood method to estimate all the unknown parameters in the proposed GARCH model in step 1.
• Simulate conditional variance from the GARCH model, specified in step 2.

The threshold decided by GARCH can detect all the cycle slips in low ionospheric cases. However, it is easily affected by the sudden increase of ionosphere and the epochs with cycle slips, which will make the threshold lose the capability of detecting small cycle slips in and around these epochs. This problem also occurs on the threshold given by the local backward method. Thus, these two threshold determination methods should be applied by setting the detection value with magnitude larger than 1 to 0 or after the detection value are checked with the threshold given by sample variance method for several times. However, both ways have flaws. The threshold determined by GARCH or the local backward method can still fail to detect small cycle slip even after applying the sample variance method for several times, while that given by setting the detection value to zero can lead to incorrect detections during the period with a number of continuous cycle slips. In order to determine a suitable threshold, another two approaches are proposed here, i.e., GARCH or local backward with a partly fixed threshold method. In these two methods, the threshold is estimated by GARCH or the local backward method, while the threshold will be fixed to one cycle when the value provided by GARCH or the local backward method is larger than one cycle. Thus, this threshold can take all sizes of cycle slips into account.

5. Empirical Test Using the Dataset with Strong Scintillation

This section is to validate the proposed methods in Section 4 using real observations affected by strong ionospheric scintillation. Ionospheric scintillation is a short-period, localized and strong ionosphere activity. The middle and bottom panels of Figure 2 show that applying the two combined signals from the proposed method in Zhao et al. [13] to this dataset leads to numerous missed detections. The top panel shows the elevation angle of satellite PRN 25 and the corresponding S4 value, which is a quantified index for amplitude scintillation. In order to reveal the validation of the proposed methods clearly, only part of the dataset affected by strong scintillation is displayed in this section. If the methods proposed in Section 4 are effective during scintillation events, it can also be applied to the data collected under weak ionospheric activities.

5.1. Validation of the Proposed Outliers Elimination Method

The performance of the proposed outliers elimination method is shown in Figure 3. After applying the proposed method to $\Delta N$, most of the outliers are removed or reduced, especially those appearing continuously, which can affect the estimation of the smoothed ionospheric bias. Although several outliers are still there, they will not affect the estimation of the ionospheric trend. Figure 3 also shows that the proposed method cannot eliminate the outliers or cycle slips with a magnitude between one to two cycles, thus a more precise method needs to be proposed to detecting small cycle slips.

5.2. Validation of the Proposed Smoothing Technique

The optimal smoothing technique should be able to reflect the ionospheric changing trend, which means the residuals should be unbiased and precise. The top panel of Figure 4 shows the ionospheric trend given by applying LOWESS and RLOWESS with a span of 20 to the outliers removed $\Delta N$. Both of these methods can generally reflect the deviation caused by ionosphere changes to $\Delta N$; however, the smoothed value of RLOWESS is much smaller than that of LOWESS at the epochs containing a cycle slip, which indicates that RLOWESS has a much stronger ability to dispense from small cycle slips. Differences between LOWESS and RLOWESS can be neglected for observations with low ionospheric effect. The bottom panel displays the RLOWESS fit residual, which is the difference between the outliers removed $\Delta N$ and the estimated ionospheric trend. Compared to the outliers removed $\Delta N$ shown in the top panel, the tilting trend is removed in the residual. In order to analyze the bias level and precision, the existing cycle slips should be removed, so the residuals with a magnitude of larger than 1 cycle are deleted. The mean of the residuals is zero, indicating that the ionospheric bias has been eliminated. The precision of the detection value can be reflected by the STD of the residuals, which is 0.23 cycle for this data sample. According to Zhao et al. [13], the detection value with this bias level can achieve a 100% SR, while the precision can guarantee a 32.34% and 0.1‰ missed detection rate for detecting cycle slips with the magnitude of one and two cycles, respectively. Overall, using the residual as the detection value can achieve a robust cycle slip detection performance for nearly all types of cycle slips, except one cycle, even under extreme ionospheric condition.

5.3. Validation of the Proposed Threshold

Figure 5 shows the detection value of the combined signal (−3, 1, 3) from the original observations and the corresponding thresholds given by the five proposed threshold determination methods. As shown in Panel A, the threshold determined by the sample variation method does not change with the variation of the ionosphere, so it loses sensitivity to the cycle slips with the magnitude of one cycle. The threshold determined by sample variation method is highly affected by the existing cycle slips. Although the sample variation method has been applied once, cycle slips can still be detected sometimes due to the improvement of the data quality after cycle slips are partly fixed. Thus, the sample variance method needs to be applied several times in order to achieve a reliable detection result. Panel B and D reveal that both GARCH and the local backward method can determine a threshold that varies with the ionospheric bias and can detect any cycle slips under the condition of weak ionospheric scintillation. With regard to the the top panel of Figure 2, the thresholds given by these two methods can exceed one cycle under the condition of the ionospheric scintillation index S4 being larger than 0.5, losing the ability to detect small cycle slips. When the partly fixed method is applied to both methods, all the cycle slips can be detected, as revealed in Panels C and E. By comparing Panels B, D and C, E, respectively, it can be found that the GARCH method performs better than the local backward method, as the threshold that GARCH is involved in determining responds more rapidly with the changes of ionospheric scintillation. This fast-changing threshold has two benefits. When the ionospheric scintillation varies from strong to weak, the threshold can still maintain the sensitivity to small cycle slips, while it can reduce the occurrence of false detection when the ionospheric scintillation gets stronger suddenly, although both benefits are not revealed during the test of this data.

6. Extensive Empirical Test Using Observations with Different Kp Indices

The performance of the proposed cycle slip detection method is evaluated by checking the detection rate of artificially added cycle slips. Data samples with large S4 value (greater than 0.5), such as the observations adopted at Section 5, are not suitable for analyzing the performance by checking the detection rate as the artificial cycle slips have a high likelihood to be added to existing cycle slips, making the detection rate be lower than it should be. Besides that, severe scintillations are only frequently observed around the equatorial and high-latitude regions, while they are rarely experienced in the mid-latitude [33]. The proposed cycle slip detection method is expected to work with all the observations globally, thus the global geomagnetic storm index (Kp index) is adopted to reveal the level of the daily ionospheric activity [18].

6.1. Data Description

The effectiveness of this proposed cycle slip detection method was extensively tested using twelve sets of 24-h GNSS observations, which were collected with different receivers and antennas as shown in

Table 3. Description of the datasets. The Kp index have been daily averaged.

Dataset	Site	City	Lat. (Deg.)	Long. (Deg.)	Receiver	Antenna	Collection Date	Kp Index
1	CUT0	Perth, Australia	−32.0039	115.8948	TRIMBLE NETR9	TRM59800.00	9 September 2017	0.625
2	CUT0	Perth, Australia	−32.0039	115.8948	TRIMBLE NETR9	TRM59800.00	11 September 2017	2.625
3	CUT0	Perth, Australia	−32.0039	115.8948	TRIMBLE NETR9	TRM59800.00	8 September 2017	6
4	PNGM	Lombrum, Papua New Guinea	−2.0432	147.3660	TRIMBLE NETR9	TRM59800.00	9 September 2017	0.625
5	PNGM	Lombrum, Papua New Guinea	−2.0432	147.3660	TRIMBLE NETR9	TRM59800.00	11 September 2017	2.625
6	PNGM	Lombrum, Papua New Guinea	−2.0432	147.3660	TRIMBLE NETR9	TRM59800.00	8 September 2017	6
7	METG	Metsahovi, Finland	60.2419	24.3841	TRIMBLE NETR9	TRM59800.00	9 September 2017	0.625
8	METG	Metsahovi, Finland	60.2419	24.3841	TRIMBLE NETR9	TRM59800.00	11 September 2017	2.625
9	METG	Metsahovi, Finland	60.2419	24.3841	TRIMBLE NETR9	TRM59800.00	8 September 2017	6
10	UNNC	Ningbo, China	29.8030	121.5568	JAVAD TR_VS	LEIAR20	9 September 2017	0.625
11	UNNC	Ningbo, China	29.8030	121.5568	JAVAD TR_VS	LEIAR20	11 September 2017	2.625
12	UNNC	Ningbo, China	29.8030	121.5568	JAVAD TR_VS	LEIAR20	8 September 2017	6

. In order to test the applicability of the proposed method to the observations globally, four stations located in different countries, including Australia, Papua New Guinea, Finland and China, were used. Except the data from Station UNNC, the rest are from the Multi-GNSS Experiment and Pilot Project (MGEX) [34]. Three of the four sites are located in Asia-Pacific region in order to have a better observation of the Beidou satellites. Stations CUT0, PNGM and METG support triple-frequency GPS, Galileo and BDS, while only triple-frequency GPS, Galileo and dual-frequency BDS are available at UNNC. Two different types of receivers and antennas were adopted, making it suitable to study the impact of the receiver and antenna on the proposed cycle slip detection method. The effect of the data epoch interval on the performance of the cycle slip detection using the modified HMW combinations has been tested in Zhao et al. [13], resulting in the cycle slips in the observations with low data rate being more difficult to detect. Thus, all the datasets adopted in this section were recorded at an interval of 30 s for the sake of simplicity.

The selected data sets have very good ionospheric diversities, varying from quiet (Kp = 0 to 2), unsettled and active ionospheric activity (Kp = 2 to 4), minor geomagnetic storm (Kp = 5), moderate geomagnetic storm (Kp = 6), strong geomagnetic storm (Kp = 7) to severe geomagnetic storm (Kp = 8). The Kp index value is released by German Research Centre for Geosciences (GFZ). The largest daily average Kp index during the year of 2016 and 2017 is only 6, which occurs on 8 September 2017; however, the Kp index during a certain period of that day, as shown in Figure 6 can reach the level of a severe storm. A severe storm like this can only be observed for approximately 60 days during every solar cycle [35]. As the Kp index is a global index, it cannot reveal the local ionospheric changes accurately. However, on account of the good geographical distribution of the selected observation stations over the world, the data sets can be affected by this severe effect, especially on the two stations PNGM and METG, which are located in the equatorial region and arctic circle, respectively. The datasets in Table 3 are listed in the order of the increase of Kp index for the observations from the same station, in order to clearly reveal the impact of the ionospheric activity on the proposed cycle slip detection method.

6.2. Ability to Detect Cycle Slips of Different Magnitudes

The adopted datasets mainly have four factors biasing the successful detection of the artificially added cycle slips. The first is the missing observation due to the loss of locks, which is because the strength of the GNSS signal drops below the tracking threshold. During the period with severe ionospheric activities, e.g., scintillation, motion of the GNSS signal pierce point along the ionosphere and motion of the receiver can lead the power fades of the signals to be under 30 dB-Hz, a condition where most of geodetic GNSS receivers stop tracking signals. The receiver needs several seconds to resume tracking, although the fade of the signal lasts less than one second. In most cases, only one or two satellites are experiencing these power fades at the same time; therefore, these satellites in our test are excluded to reduce the danger of adding the artificial cycle slips to the missing observations. Associated with these large power fades, another two factors are the existing cycle slips and half-cycle phase jumps [36,37]. The existing cycle slips in the datasets can affect the detection by fixing the cycle slip to a value that is not equal to the expected one. Thus, these cycle slips need to be corrected by the proposed cycle slip detection method before the data can be used in the following analysis. As our proposed cycle slip detection method cannot provide a reliable detection and correction of these half-cycle phase jumps, the influence of the half-cycle phase jump is still there even though the datasets are preprocessed by our proposed method. One possible way to improve the detection ability to these half-cycle phase jumps is to change the threshold from four times STD to three times STD, in which one drawback is leading to undetected cycle slips with regard to the noise of the detection value. The undetected cycle slips here include both the cases in which the algorithm fails to alarm the cycle slip and those fixing the cycle slips to the wrong magnitudes.

In order to analyze the performance of the proposed method in detecting different magnitude cycle slips, the datasets are configured as follows. Firstly, the observations from GPS satellites supporting only dual-frequency signals and those with loss of lock are excluded. Before preprocessing with the proposed cycle slip detection method to detect and correct any existing cycle slips, the remaining observations are filtered with the elevation cutoff angle, which is set as 10${}^{\circ }$ in this test, to reduce the effect of multipath. Then, the observations are added the cycle slip sets with different magnitudes, which are simulated as follows. Firstly, a matrix with seven cycle slip sets with the magnitude of one cycle on each frequency is simulated, as listed in

Table 4. Simulated cycle slip sets with the magnitude of one cycle. The unit is in cycles.

$\mathbf{Mag}.\mathbf{on}{\mathit{f}}_1$	$\mathbf{Mag}.\mathbf{on}{\mathit{f}}_2$	$\mathbf{Mag}.\mathbf{on}{\mathit{f}}_3$
0	0	1
0	1	0
0	1	1
1	0	0
1	0	1
1	1	0
1	1	1

. Then, this matrix is multiplied by the magnitude ranging from 2 to 10 to obtain other cycle slip sets with different magnitudes. The simulated cycle slips of all these ten matrices are manually added to the carrier phase observations of each satellite, respectively. For the proposed cycle slip detection method, the difficulty of successful detection mainly depends on the magnitude of the cycle slip. Smaller magnitude is harder to detect. Therefore, any other cycle slip sets with magnitudes above 10 cycles will not be simulated, e.g., (77, 66, 0), which is one of the ionosphere-free cycle slip sets for GPS observations.

The number of simulated cycle slips and success rates (SRs) for GPS, Galileo and BDS are summarized in

Table 5. Results of cycle slip detection for GPS, Galileo and BDS. The SRs are obtained with the magnitude of one cycle.

Dataset	Kp Index	GPS		Galileo		BDS
		No. of Simulations	SR (%)	No. of Simulations	SR (%)	No. of Simulations	SR (%)
1	0.625	84	100	112	100	98	100
2	2.625	84	100	105	100	98	100
3	6	84	100	105	100	98	100
4	0.625	70	100	105	99.05	77	100
5	2.625	77	100	105	99.05	84	98.81
6	6	84	98.81	98	100	84	98.81
7	0.625	84	100	112	100	70	98.57
8	2.625	84	100	105	100	70	98.57
9	6	77	94.81	84	97.62	56	100
10	0.625	77	100	91	100	N/A	N/A
11	2.625	77	100	77	100	N/A	N/A
12	6	84	100	84	100	N/A	N/A

. For the cycle slips with different magnitudes, the epoch where the cycle slip is added is the same for a certain satellite. The number of cycle slips added to each satellite arc is 7. The reason for the difference in the number of added cycle slips each day is the exclusion of satellites with the loss of lock caused by the local ionospheric activity, such as scintillation. The number of observed triple-frequency satellites on PNGM, located in the equatorial region, and METG, located in the arctic region, are 12 and 16 for GPS and Galileo, respectively. Comparing the number of added cycle slip sets to the dataset 4 to 9, the number of omitted satellites for PNGM decreases, while that for METG increases with the growth of Kp index. This phenomenon can be clearly observed in GPS and BDS satellites, but may be not in Galileo satellites observed at PNGM. However, taking the detection results into account, it can be found that one Galileo satellite is still affected by local ionospheric activity in datasets 4 and 5. The reason for this phenomenon is due to the correlation of the Kp index and the scintillation parameters. The correlations between the scintillation occurrence and the global geomagnetic storm are different to latitudes. In high latitude areas, such as Alaska and Northern Europe, strong scintillation events have a strong linear correlation with geomagnetic storms [38,39], while the geomagnetic disturbance does not trigger the occurrence of scintillation in low latitude areas, such as Sanya [40].

The detection results in Table 5 clearly illustrate that the proposed cycle slip detection method is effective in detecting cycle slips in the three GNSS systems despite a few undetected cycle slips. This table displays the SRs obtained with the magnitude of one cycle. It should be noted that the SRs with the magnitudes of 2–10 cycles are the same as that given with one cycle. The undetected cycle slips always occur in the same epoch where there is a half-cycle phase jump. Overall, except for a few undetected cycle slips, the proposed method can detect and repair the cycle slips robustly even during high ionospheric activities, especially in the middle latitude regions.

6.3. Application in PPP

As the proposed cycle slip detection method requires observations from one receiver only and checks the cycle slips on a satellite-by-satellite basis, it is applicable in the PPP process. In order to take full advantage of all the available triple-frequency satellites from GPS, Galileo and BDS, dataset 3 located in the Asia-Pacific region is selected from Table 3 to serve as an application situation of PPP. Dataset 3 is a 24-h GNSS observations with a 30 s data epoch interval. As the proposed method supports only triple-frequency signals, all the dual-frequency GPS satellites are excluded. As listed in

Table 6. Simulated cycle slip sets and the corresponding start time in dataset 3. The unit of the cycle slip is in cycle. The time is GPS time.

Simulated Cycle Slip Set			Start Time
$\mathbf{Mag}.\mathbf{on}{\mathit{f}}_{\mathbf{1}}$	$\mathbf{Mag}.\mathbf{on}{\mathit{f}}_{\mathbf{2}}$	$\mathbf{Mag}.\mathbf{on}{\mathit{f}}_{\mathbf{3}}$
1	1	0	02:45
1	0	0	05:45
0	2	0	08:45
77	60	0	11:45
9	7	0	14:45
3	4	0	20:45

, six cycle slip sets are simulated to assess the effect of the cycle slip on PPP. These six pairs include three small cycle slips, two large cycle slips and one ionosphere-free cycle slip set for GPS, which was not considered in Section 6.2. In order to reveal the improvement in the precision of PPP after the cycle slips are repaired by the proposed method, each cycle slip set is added to all the signals at the same epoch of all the satellites and iterated 12 times at an interval of every four epochs. The start times of adding each cycle slip set are shown in Table 6. Finally, dataset 3 is configured to the observations with intensive cycle slips over several half-hour time spans.

RTKLIB is applied to perform PPP in this section. RTKLIB is an open source program package which can conduct precise positioning with GNSS dual-frequency signals on $f_1$ and $f_2$ [41], which is the reason why no cycle slips are added to $f_3$ signal in Table 6. In order to perform PPP in kinematic mode using RTKLIB, the precise satellite ephemeris and clock files are required besides the previous configured observation file [42]. The positioning options are configured as follows. The cutoff elevation angle is set to 15 degrees, in order to eliminate the effect of multipath on the accuracy of PPP. Solid earth tides correction is applied. Ionosphere-free linear combination with dual-frequency measurements is used for ionospheric correction, while the tropospheric parameters (ZTD) are estimated as the states of the extended Kalman filter. The antenna phase center offset (PCO) and phase center variation (PCV) of the satellite antenna and the receiver antenna are corrected using the ANTEX antenna parameters file provided by the International GNSS Service (IGS) [43]. The errors caused by earth rotation are corrected by GNSS final combined Earth Rotation Parameter (ERP) Product available from the Crustal Dynamics Data Information System (CDDIS) [44]. The default values provided by RTKLIB are adopted for any other options.

Our proposed method is integrated into RTKLIB by preprocessing cycle slip detection and correction before conducting PPP. Figure 7 shows the coordinate errors in longitude (Lon), latitude (Lat) and altitude (Alt) by PPP using the data with and without being preprocessed by our proposed method, respectively. The coordinate errors are obtained by subtracting the single PPP static position value, which is given by RTKLIB using all the available satellites, to the epoch-by-epoch kinematic PPP solutions. Small cycle slips of 1 or 2 cycles. e.g., the ones introduced at 2:45 a.m., 5:45 a.m. and 8:45 a.m. lead to coordinate errors of up to 20 m in horizontal position and up to 50 m in the elevation, so RTKLIB can only perform reliable single point positioning in the absence of cycle slips and RTKLIB does not have the ability to repair cycle slips. Because the adopted version of RTKLIB is 2.4.2, which does not support PPP ambiguity resolution (AR), it needs at least half an hour to converge to the positioning precision of PPP in most cases. After the cycle slips are detected and corrected, the coordinate errors in position would be in the sub-meter level. The precision of PPP, which is the STD of the coordinate errors, can achieve 0.2 m in the horizontal direction and 0.5 m in the vertical direction. To conduct triple-frequency PPP, the number of available satellites should be six to seven in our experiment. At around 1:00 p.m., the coordinate errors increase to several tens of meters in all the three directions, which is because RTKLIB excluded most of the satellites. The exclusion might be due to the large residuals of these satellites. The detailed reason for this deviant exclusion needs more research, but it does not affect us drawing the final conclusion. Overall, RTKLIB detects but does not repair cycle slips, while our proposed method cycle slip detection method is effective in fixing cycle slips and then decreasing coordinate errors of PPP caused by cycle slips.

7. Conclusions

Since the cycle slip detection method proposed in Zhao et al. [13] was invalidated under the condition of strong scintillation, this paper expanded the application of the method to the scenario of high ionospheric activities by reselecting the combined signals and modifying the detection process. The combined signals which were least affected by the measurement noise were selected as the optimal, without considering the ionospheric bias. Strong ionospheric scintillations can produce many cycle slips with large magnitude in succession. These large cycle slips had to be removed in order to reduce their biases on the estimated ionospheric trend in the HMW combination by RLOWESS. The residual of subtracting the estimated ionospheric trend from the HMW combination was adopted as the detection value to determine the magnitude of cycle slips. Under the condition of strong ionospheric activities, the detection value will fluctuate greatly with the change of ionospheric, which affects the threshold on the judgment of cycle slips. Thus, five threshold determination methods were proposed and examined, including the sample variation method, GARCH, the local backward method, GARCH with the partly fixed method and the local backward with the partly fixed method. All the proposed methods in the detection process, including the outliers elimination method, applying RLOWESS to estimate the ionospheric trend and the threshold determination methods, were validated using the observations collected under strong ionospheric scintillation with the following conclusions.

• The proposed outliers elimination method was effective in removing the outliers, but failed in detecting the cycle slips with a magnitude of one or two cycles.
• RLOWESS could estimate the ionospheric trend in the HMW combinations without the influence by the dispersion value in the detection value.
• The sample variation method was not suitable to be adopted under strong ionospheric scintillation, as cycle slips could still be detected after one detection process. Both GARCH and the local backward method can correctly repair all the detected slips, but they might lose the sensitivity to the cycle slips with the magnitude of one or two cycles under strong ionospheric scintillation. All the cycle slips can be detected and correctly repaired using either the GARCH with the partly fixed method or the local backward with the partly fixed method. Compared to the threshold determined by the local backward method, that given by GARCH has a quicker response to the changes of the ionosphere. Thus, the GARCH with the partly fixed threshold determination method is optimal.

The effectiveness of the whole proposed cycle slip detection method was extensively further tested using twelve GNSS datasets, which were collected with a data rate of 30 s from the stations located in different countries globally, on the dates with different levels of Kp index. A high success rate was achieved despite a few undetected cycle slips, due to the half-cycle phase jump brought by scintillation. As the proposed method detects and corrects cycle slips on a satellite-wise basis, it can be used to improve the precision and continuity of PPP.

The proposed cycle slip detection and correction method mainly have three problems, which are left for future research. The first is to determine the optimal threshold for the observations under strong ionosphere scintillation. The two relatively tight threshold determination methods, i.e., the GARCH with the partly fixed threshold method and the local backward threshold method, can edit the GNSS data as cleanly as possible, but they might be overly conservative in deciding whether a cycle slip has occurred during the process of PPP AR, as frequently fixing cycle slips would lead to unnecessary ambiguity resets, which would weaken the positioning geometry and reliability [45]. The other problem is the handling of the half-cycle phase jump. Correctly detecting and repairing all the half-cycle phase jumps will greatly improve the success rate of this cycle slip detection and correction method. The last drawback is that the proposed method cannot detect and repair cycle slips instantaneously because the adopted smoothing technique needs both backward and forward observations to estimate the optimal ionospheric trend in the detection value. Thus, a smoothing technique requiring only backward data needs to be proposed in order to detect cycle slips in real time under high ionospheric activities.