Sage Journals: Discover world-class research

Abstract

When carrier phase observations are applied to high-precision positioning, how to handle the cycle slip is an unavoidable problem. For cycle slip correction, detection combination noise and the ionospheric delay are two crucial factors. Specifically, the drastic changes in the ionosphere and the increased noise of code observations will increase the failure probability of cycle slip detection. To reduce the influence of code observation noise and ionospheric bias, a novel cycle slip detection method for BDS-3 satellites is proposed. Considering that code measurement noise is closely related to the satellite elevation angle, an elevation-based model is built to evaluate the code measurement noise. Firstly, two modified code-phase combinations are selected optimally based on 1% missed detection rate and 99% success detection rate to minimize the effects of code measurement noise. However, the second modified code-phase combination is more affected by ionospheric delay bias, so ionospheric bias of current epoch needs to be corrected. To reduce the influence of ionospheric bias, two moving windows of time-differenced ionospheric delay are introduced to correct the ionospheric bias of the second code-phase combination. Experiments with BeiDou-3 data are implemented in three different scenarios. To verify the effectiveness of the algorithm in the environment of high code observations noise, Gaussian noise is added to the code observations in the first scenario, and the results demonstrate that the success rate of cycle slip detection and repair is still greater than 95% when the standard deviation of Gaussian noise is 0.8 m. The second scenario is carried out under low ionospheric activity, and results indicate that the proposed method significantly reduces the times of failed detection and repair. Moreover, in the third scenario, BeiDou-3 data with cycle slips of different types under high ionospheric activity are tested, and all cycle slips can be correctly detected and corrected.

Keywords

BeiDou-3 cycle slip triple-frequency signals elevation-dependent model ionospheric bias compensation code observation noise

Introduction

Global navigation satellite system (GNSS) plays an important role in military confrontations,¹ transportation,^2
–4 and measurement.⁵ The Chinese BeiDou Navigation Satellite System (BDS) includes regional system (BDS-2) and global system (BDS-3). BDS-2 has already been accomplished in 2012 and triple-frequency signals can be transmitted by all satellites.^6,7 BDS-3 can provide global services and it is expected to be completely built in 2020. The release of four signals including B1I, B3I, B1C, and B2a is beneficial to the ambiguity resolution, cycle slip detection, and high-precision positioning.⁶

Millimeter precision is the advantage of carrier phase observations compared to code observations. However, the existence of cycle slips is an unavoidable problem in the application of carrier phase observations. Cycle slips often occur in poor environment, such as signal interruption and interference. The quality of satellite data preprocessing is directly related to the accuracy and reliability of navigation and positioning. If cycle slips are ignored or cannot be correctly repaired, it will bring deviations to the ultimate results.

For cycle slip detection, the basic idea is to obtain a sequence of detections that reflects the change of the cycle slip, and the position of the cycle slip can be determined from the detection sequence. For single-frequency observations, Beutler et al.⁸ proposed a polynomial fitting method to detect cycle slip. Meanwhile Hofmann et al.⁹ put forward a high-order inter-epoch phase differentiation method. However, they need several observations from previous epochs for detection. Zangeneh-Nejad et al.¹⁰ improved the single-frequency GPS cycle slip correction method, which is based on the generalized likelihood ratio test. Momoh et al.¹¹ proposed to detect cycle slip through the adaptive difference sequences of observations, which can achieve single-frequency cycle slip detection and correction with the slips of one cycle.

For GPS dual-frequency data, Hatch–Melbourne–Wübbena (HMW) combinations^12
–14 are extensively applied in cycle slip detection due to their excellent characteristics. Blewitt¹⁵ developed the Turbo Edit algorithm and used undifferenced observations to detect cycle slip. Specifically, Turbo Edit algorithm consists of ionospheric combination and HMW combination. Considering the gross error of code observations may result in a wrong detection result, a method is proposed by Li et al.,¹⁶ which is relied on satellite orbit and smoothed code observations. Ionospheric total electron contents rate (TECR) method was proposed by Liu,¹⁷ which was based on smoothly changing ionosphere. However, the TECR method does not work well under high ionospheric activity. To solve this problem, Cai et al.¹⁸ proposed improved TurboEdit method, which based on the second-order, time-difference phase ionospheric residual (STPIR) algorithm, and the proposed approach was able to determine the cycle slip of each signal under high ionospheric activities. Banville and Langley¹⁹ developed the dual-frequency model, and least-squares adjustment was adopted to reduce the impact of the ionosphere in the case of ionospheric disturbance. Hu et al.²⁰ improved Turbo Edit algorithm and a polynomial function to fit the ionospheric delays in short times.

Up to now, many satellites have supported the transmission of triple-frequency signals, which benefits ambiguity resolution^21,22 and cycle slip detection. Triple-frequency signals can present more high-performance observation combinations with longer wavelengths, lower noise, and ionospheric errors. For GPS triple-frequency observations, Dai et al.²³ constructed geometry-free (GF) phase observations as detection sequence. De Lacy et al.²⁴ proposed a real-time detection method and five first-order time-difference GF linear combinations were selected to detect different types of cycle slips. Zhao et al.²⁵ and Gu and Zhu²⁶ extended three carrier ambiguity resolution algorithm from ambiguity resolution to cycle slip correction, and extra-wide-lane, wide-lane, and narrow-lane (NL) combinations are utilized to determine cycle slip.

For BDS-2 triple-frequency cycle slip detection, three GF combinations are used by Huang et al.²⁷ to determine the position of the cycle slip, and the size of the cycle slip will be determined by least-squares ambiguity decorrelation adjustment.²⁸ Xiao et al.²⁹ modified the geometry-based detection model and a new strategy was put forward to reduce false alarms. Zeng et al.³⁰ divided the special cycle slip groups into two types and used three GF linear combinations which were collected to detect and repair. Li and Melachroinos³¹ proposed an enhanced repair algorithm, and the real-time implementation can be realized by Kalman filter. To decrease the influence of code observation noise, Zhao et al.³² improved the HMW combination and determined different coefficients for satellites of different orbital types. Xiao et al.³³ applied three GF combinations and 20° was applied as the threshold of the elevation angle to ensure a high success rate. To eliminate the bias of ionosphere, prediction functions were used by Yao et al.³⁴ and Li et al.³⁵ to predict ionospheric variation. Pu et al.³⁶ used moving window filter function to predict ionospheric variation, and the cycle slip will be re-search when the predicted value exceeds the threshold. Xiao et al.³³ used previous $Δ I$ to compensate the NL combination of current epoch when the second-order time-differenced ionospheric delay ( $Δ Δ I$ ) became abnormal. The second-order time-difference algorithm was used by Liu et al.³⁷ to identify cycle slip during high ionospheric activity; however, the approach will increase noise of the detection combinations.

In general, there are many methods for cycle slip detection; code measurement noise and ionospheric bias are two crucial factors. Many methods assumed that code measurement noise is constant; however, code measurement noise will significantly increase when the satellite elevation angle becomes lower.³⁸ Based on the assumption that the ionosphere changes smoothly, the failure probability of cycle slip detection will increase. Therefore, a new method for cycle slip detection and correction is presented in this article. Firstly, in order to evaluate the code measurement noise, an elevation-based model is built. Secondly, HMW combination is used as the first detection combination due to its extremely high success rate. To minimize the combined observation noise, two optimally code-phase combinations are selected. Then, in order to decrease the impact of ionospheric bias, two moving windows of time-differenced ionospheric delay are introduced to correct the ionospheric bias.

The structure of this article is as follows. Triple-frequency linear combination model is introduced in the second section. The third section analyzes the error affecting the detection of cycle slip, and elevation-based model is built. In fourth and fifth sections, three cycle slip detection combinations and corresponding coefficients are determined. In the sixth section, two moving windows are proposed to correct the ionospheric bias of the third detection combination. In the seventh section, the approach is tested in three different scenarios. Finally, conclusions are summarized in the last section. Appendix 1 shows all symbols and their definitions used in this article.

Triple-frequency linear combination

Code and phase observations at frequency f_i $(i = 1, 2, 3)$ can be expressed as³⁹

P_{i} = ρ + c (d t_{r} - d t^{s}) + T + k_{1 i} I + ε_{P_{i}}

Φ_{i} = ρ + c (d t_{r} - d t^{s}) + T - k_{1 i} I - N_{i} λ_{i} + ε_{Φ_{i}}

where P and Φ are code and phase observations, respectively. ρ denotes the geometric distance from satellite to receiver, and it is a frequency-independent item. c is the speed of light, $d t_{r}$ is the receiver clock error, and $d t^{s}$ is the satellite clock error. T denotes the error of the troposphere. N is the integer ambiguity and the unit is cycle. I is the ionospheric bias of frequency f ₁, and $k_{1 i} = f_{1}^{2} / f_{i}^{2}$ denotes ionospheric scale factor (ISF) for different frequencies. λ is the wavelength; $ε_{P}$ and $ε_{Φ}$ denote the measurement noise of code and phase observations, respectively.

The triple-frequency carrier phase combination can be expressed as²¹

Φ_{(i, j, k)} = \frac{i f_{1} Φ_{1} + j f_{2} Φ_{2} + k f_{3} Φ_{3}}{i f_{1} + j f_{2} + k f_{3}} = ρ - β_{(i, j, k)} I - λ_{(i, j, k)} N_{(i, j, k)} + ε_{Φ (i, j, k)}

where $i, j, k$ are the coefficients of the combination and they are integers, $β_{(i, j, k)}$ is the ISF of the triple-frequency phase combination, and it can be derived as

β_{(i, j, k)} = \frac{f_{1}^{2} (i / f_{1} + j / f_{2} + k / f_{3})}{i f_{1} + j f_{2} + k f_{3}}

The wavelength and integer ambiguity of the combination are

λ_{(i, j, k)} = \frac{c}{i f_{1} + j f_{2} + k f_{3}}

N_{(i, j, k)} = i N_{1} + j N_{2} + k N_{3}

Error analysis

In the time-differenced model, the bias of troposphere, the clock error, and the hardware delay will be eliminated. However, the code measurement noise and the ionospheric bias are difficult to be eliminated by the time-differenced model, which are the main obstacle to the detection of cycle slip. To evaluate the code measurement noise and ionospheric bias, we analyzed the data of the six stations from International GNSS Service (IGS). The information of the six stations is listed in Table 1.

Table 1.

The information of the six stations from IGS.

Station name	Location (city, country)	Observation time	Sampling interval (s)
ULAB00MNG (ULAB)	Ulaanbataar, Mongolia	June 19, 2019 00:00:00–23:59:30	30
SGOC00LKA (SGOC)	Narahenpita, Colombo, Sri Lanka	June 19, 2019 00:00:00–23:59:30	30
POTS00DEU (POTS)	Potsdam, Germany	June 3, 2019 00:00:00–23:59:30	30
MIZU00JPN (MIZU)	Mizusawa, Japan	June 3, 2019 00:00:00–23:59:30	30
URUM00CHN (URUM)	Urumqi, China	June 19, 2019 00:00:00–23:59:30	30
SUTM00ZAF (SUTM)	Sutherland, South Africa	June 19, 2019 00:00:00–23:59:30	30

IGS: International global navigation satellite system service.

Satellite elevation-based measurement noise model of code observations

The measurement noise is correlated with observation environments, signal–noise ratio, and satellite elevation angle. Multipath effect is the main cause of the measurement noise. In this article, code multipath errors are used to evaluate code measurement noise. Multipath combination related to code range can be expressed as⁴⁰

M P_{i} = p_{i} - \frac{f_{i}^{2} + f_{j}^{2}}{f_{i}^{2} - f_{j}^{2}} Φ_{i} + \frac{2 f_{j}^{2}}{f_{i}^{2} - f_{j}^{2}} Φ_{i} - B_{i}

where f is the frequency. $i, j (i, j = 1, 2, 3)$ represent different frequency and $i \neq j_{.} B_{i}$ mainly contains the linear combination of the phase ambiguities and the constant part of hardware delays. If no cycle slip occurs, B_i can be considered a constant.

The standard deviation of multipath error corresponding to each elevation angle is selected to create multipath-elevation sequence for modeling.

By analyzing and processing real data of C20–C36 at six IGS stations, we obtain the elevation-based measurement noise model of code observations

σ_{p_{i}} = \frac{c_{0}}{{sin}^{2} θ + c_{1}} (5^{\circ} \leq θ \leq 90^{\circ})

where $σ_{p}$ is standard deviation of code observations measurement noise, θ is the satellite elevation angle, and c ₀, c ₁ are coefficients of the model. $5^{\circ}$ is set as the cutoff elevation angle, because the data are easily interrupted when the elevation angle is less than $5^{\circ}$ . Table 2 lists the values of c ₀ and c ₁ for BDS-3 satellites. The unit of c ₀ and c ₁is meter.

Table 2.

Coefficients of satellite elevation-based model for BDS-3 satellites.

B1C		B1I		B3I		B2a
c ₀	c ₁	c ₀	c ₁	c ₀	c ₁	c ₀	c ₁
0.17	0.16	0.446	0.364	0.23	0.36	0.14	0.226

In this article, all of the carrier phase measurement noises $σ_{Φ_{1}}, σ_{Φ_{2}}$ , and $σ_{Φ_{3}}$ are assumed to be equal ( $σ_{Φ_{1}} = σ_{Φ_{2}} = σ_{Φ_{3}} = 0.003 m$ ) for each satellite of BDS. B1C, B3I, and B2a are selected as three frequencies of cycle slip detection. $f_{1}, f_{2}$ , and f ₃ represent the frequency of B1C, B3I, and B2a, respectively.

Time-differenced ionospheric delay

The ionospheric variation is a non-negligible factor for cycle slip detection. The ionospheric variation of two consecutive between-epoch on the B1C signal can be calculated as follows²⁵

Δ I = \frac{Δ Φ_{1} - Δ Φ_{3}}{k_{13} - 1}

where Δ denotes the first-order time-difference. $Δ I$ is a steady change in high ionospheric activity.

The second-order time-differenced ionospheric delay is

Δ Δ I = Δ I (t) - Δ I (t - 1)

where $Δ Δ$ represents the second-order time-difference and t is the current epoch.

De Lacy et al.²⁴ considered $Δ I$ as normal distribution with mean zero, and standard deviation of $Δ I$ is used to evaluate the magnitude of the ionospheric variation. Standard deviation of $Δ I$ in 24 h for BDS-3 C20–C36 satellites is computed, and the result is shown in Figure 1.

Figure 1.

Standard deviation of $Δ I$ for BDS-3 satellites.

Figure 1 reflects the overall level of ionospheric variation within a day, and it shows that standard deviation of $Δ I$ is less than 0.02 m for BDS-3 C20–C36 satellites.

Cycle slip detection combinations

First detection combination

The first detection combination is HMW combination. HMW combination is an optimal choice due to longer wavelength and higher reliability. In HMW combination, geometric distance, the bias of troposphere, and ionosphere will be eliminated. For HMW combination, the cycle slip can be accurately detected. HMW combination for BDS-3 satellites can be expressed as

N_{(0, 1, - 1)} = N_{2} - N_{3} = (\frac{f_{2} P_{2} + f_{3} P_{3}}{f_{2} + f_{3}} - Φ_{(0, 1, - 1)}) / λ_{(0, 1, - 1)}

where $Φ_{(0, 1, - 1)} = (f_{2} Φ_{2} - f_{3} Φ_{3}) / (f_{2} - f_{3})$ .

By differencing equation (11) between adjacent epochs, the magnitude of the cycle slip can be derived as

{\begin{cases} Δ N_{(0, 1, - 1)} = Δ N_{2} - Δ N_{3} = (\frac{f_{2} Δ P_{2} + f_{3} Δ P_{3}}{f_{2} + f_{3}} - Δ Φ_{(0, 1, - 1)}) / λ_{(0, 1, - 1)} \\ Δ {\hat{N}}_{(0, 1, - 1)} = round [Δ N_{(0, 1, - 1)}] \end{cases}

where $Δ {\hat{N}}_{(0, 1, - 1)}$ is the integer cycle slip of the first combination. When $| Δ N_{(0, 1, - 1)} | > 0.5$ cycle, we recognize that a cycle slip occurs.

The value of cycle slip can be easily determined, and HMW combination can be repaired immediately. The combined signal $Δ {\hat{Φ}}_{(0, 1, - 1)}$ without cycle slip can be expressed as

Δ {\hat{Φ}}_{(0, 1, - 1)} = Δ Φ_{(0, 1, - 1)} + λ_{(0, 1, - 1)} Δ {\hat{N}}_{(0, 1, - 1)}

Assuming $Δ N_{(0, - 1, 1)}$ is normally distributed with zero mean and standard deviation at $σ_{Δ N_{(0, 1, - 1)}}$ .²⁴ The success detection rate of HMW combination is

P_{(0, 1, - 1)} = P (| Δ N_{(0, 1, - 1)} - Δ {\hat{N}}_{(0, 1, - 1)} | < 0.5)

From Figure 2, we can know that $P_{(0, 1, - 1)}$ can achieve 100% when the elevation degree exceeds $30^{\circ}$ . When the elevation degree is $5^{\circ}$ , $P_{(0, 1, - 1)}$ is 99.14%. It can be seen that HMW combination is still reliable at low elevation angles.

Figure 2.

The success detection rate of the first detection combination.

Second and third detection combinations

The traditional code-phase cycle slip detection combination is defined as

Δ N_{(i, j, k)} = (l Δ P_{1} + m Δ P_{2} + n Δ P_{3} - Δ Φ_{(i, j, k)}) / λ_{(i, j, k)}

where $l, m, n$ are the coefficients of the differential code observations. The disadvantage of this detection combination is that it will be seriously affected by the measurement noise of code observations. Two modified code-phase combinations are applied as the second and third detection combinations. Two virtual carrier phase observations are used. The cycle slip of the two combinations can be shown as follows

{\begin{cases} Δ N_{(i, j, k)} = (a Δ \bar{P} + b Δ P_{3} + c Δ {\hat{Φ}}_{(0, 1, - 1)} - Δ Φ_{(i, j, k)}) / λ_{(i, j, k)} \\ Δ {\hat{N}}_{(i, j, k)} = round [Δ N_{(i, j, k)}] \end{cases}

where $Δ \bar{P} = (Δ P_{1} + Δ P_{2} + Δ P_{3}) / 3$ , $a, b, c$ are the coefficients, and $Δ {\hat{Φ}}_{(0, 1, - 1)}$ is the phase combination whose cycle slip has been repaired in the first detection combination. By rounding $Δ N_{(i, j, k)}$ , an integer cycle slip $Δ {\hat{N}}_{(i, j, k)}$ can be obtained.

Assuming $σ_{Φ_{1}} = σ_{Φ_{2}} = σ_{Φ_{3}} = 0.003$ m, $σ_{Δ P_{1}}, σ_{Δ P_{2}}, σ_{Δ P_{3}}$ and $σ_{Δ {\hat{Φ}}_{(0, 1, - 1)}}$ are shown in Figure 3.

Figure 3.

Standard deviations of measurement noises $Δ P_{1}$ , $Δ P_{2}$ , $Δ P_{3}$ , $Δ \bar{P}$ , and $Δ {\hat{Φ}}_{(0, 1, - 1)}$ .

From Figure 3, we can know that standard deviations of measurement noises $Δ \bar{P}$ , $Δ P_{3}$ , and $Δ {\hat{Φ}}_{(0, 1, - 1)}$ are smaller than those of $Δ P_{1}$ , $Δ P_{2}$ , and $Δ P_{3}$ , especially at low elevation angles, which effectively reduces the influence of code observation noise on detection combinations. Since cycle slip of the first combination can always be successfully fixed, $σ_{Δ {\hat{Φ}}_{(0, 1, - 1)}}$ can be guaranteed at a normal level. The second and third detection combinations are GF combinations; hence, the coefficients $a, b, c$ satisfy the following equation

a + b + c = 1

Determination of the second and third detection combinations

In this section, the second and third detection combinations and the corresponding coefficients will be determined. The level of the time-differenced ionospheric delay and combination noise must be considered when choosing the optimal detection combinations.

According to equations (4) and (16), the ionospheric delay bias of the second and third combinations can be derived as

{\begin{cases} β_{Δ N_{(i, j, k)}} = \frac{a}{3} (β_{(1, 0, 0)} + β_{(0, 1, 0)} + β_{(0, 0, 1)}) + b β_{(0, 0, 1)} - c β_{(0, 1, - 1)} + β_{(i, j, k)} \\ δ_{Δ N_{(i, j, k)}} = β_{Δ N_{(i, j, k)}} \frac{Δ I}{λ_{(i, j, k)}} \end{cases}

where $β_{Δ N_{(i, j, k)}}$ is the ISF of the second and third detection combinations.

We can obtain the noise of the second and third detection combinations as

σ_{Δ N_{(i, j, k)}} = \frac{\sqrt{a^{2} σ_{Δ \bar{P}}^{2} + b^{2} σ_{Δ P_{3}}^{2} + c^{2} σ_{Δ {\hat{Φ}}_{(0, 1, - 1)}}^{2} + σ_{Δ Φ_{(i, j, k)}}^{2}}}{λ_{(i, j, k)}}

Success detection rate and missed detection rate

Assuming $Δ N_{(i, j, k)}$ without cycle slip obeys a normal distribution with mean $δ_{Δ N_{(i, j, k)}}$ and standard deviation $σ_{Δ N_{(i, j, k)}}$ . $4 σ_{Δ N_{(i, j, k)}}$ is chosen as the threshold. When $| Δ N_{(i, j, k)} - δ_{Δ N_{(i, j, k)}} | > 4 σ_{Δ N_{(i, j, k)}}$ , it is considered that a cycle slip occurs. The success detection rate P_S can be expressed as³²

P_{S} = \int_{- 4 σ_{Δ N_{(i, j, k)}}}^{4 σ_{Δ N_{(i, j, k)}}} \frac{1}{\sqrt{2 π} σ_{Δ N_{(i, j, k)}}} exp (- \frac{{(x - δ_{Δ N_{(i, j, k)}})}^{2}}{2 σ_{Δ N_{(i, j, k)}}^{2}}) d x

The missed detection rate P_M represents the probability that the detection combination does not detect the existence of a cycle slip when cycle slip occurs, it can be indicated as³²

P_{M} = P (Z > (1 + δ_{Δ N_{(i, j, k)}} - 4 σ_{Δ N_{(i, j, k)}}))

We analyzed the relationship among $δ_{Δ N_{(i, j, k)}}$ , $σ_{Δ N_{(i, j, k)}}$ , and P_S or P_M , and the results can be seen in Figure 4.

Figure 4.

The three-dimensional relationship diagram. (a) The three-dimensional relationship diagram of $δ_{Δ N_{(i, j, k)}}$ , $σ_{Δ N_{(i, j, k)}}$ , and P_S . (b) The three-dimensional relationship diagram of $δ_{Δ N_{(i, j, k)}}$ , $σ_{Δ N_{(i, j, k)}}$ , and P_M .

Figure 4(a) shows the relationship among $δ_{Δ N_{(i, j, k)}}$ , $σ_{Δ N_{(i, j, k)}}$ , and P_S . The yellow surface represents equation (20); two red lines are boundaries between the surface and the plane $P_{M} = 99 %$ . The absolute value of slopes of the two red lines $| k | = | k' | = 1.6753$ . If 99% is adopted as the acceptable degree of success rate, the following condition must be satisfied:

| δ_{Δ N_{(i, j, k)}} | \leq 1.6753 σ_{Δ N_{(i, j, k)}}

Figure 4(b) shows the relationship among $δ_{Δ N_{(i, j, k)}}$ , $σ_{Δ N_{(i, j, k)}}$ , andP_M . It is obvious that the missed detection rate is only related to the noise of detection combination. If the noise of detection noise is too large, the cycle slip is likely to be overwhelmed and cannot be detected. The blue vertical line represents that the missed detection rate achieves 1% when $σ_{Δ N_{(i, j, k)}} = 0.1581$ . To ensure $P_{M} \leq 1 %$ , the following condition must be satisfied:

σ_{Δ N_{(i, j, k)}} \leq 0.1581

Selection of optimal detection combinations

According to the above analysis, 99% success detection rate and 1% missed detection rate are adopted as the selection criterion of the second and third detection combinations. Hence, the second and third detection combinations must satisfy equations (17), (22), and (23). In addition to these three conditions, the three detection combinations must be linearly independent. To select the optimal virtual signals $Φ_{(i, j, k)}$ and the corresponding coefficients $a, b, c$ , the values $i, j, k$ of $Φ_{(i, j, k)}$ are tested within the range of [−20, 20], as shown in Figure 5. Here, $Δ I$ is configured as 0.02 m.

Figure 5.

Flow diagram of selecting the optimally detection combinations and corresponding coefficients.

As shown in Figure 5, each group of $i, j, k$ will be tested for satellite elevation angle $θ = 5^{\circ}$ . If $i, j, k$ and $a, b, c$ can satisfy the lowest elevation angle, other elevation angles can also be satisfied. The minimum of $σ_{Δ N_{(i, j, k)}}$ will be calculated. Only when the two conditions (22) and (23) are satisfied at the same time, the values of $i, j, k$ and $a, b, c$ can be output. Otherwise, the next group of $i, j, k$ will be tested.

When the above process is completed, all possible combinations will be output. Notice that the three detection combinations must be linearly independent, so the signals (1, −4, 3) and (−4, 6, −1) are selected as the second and third detection combinations, respectively, and the specific results are shown in Table 3.

Table 3.

Coefficients, wavelengths, and ISF of the second and third detection combinations.

Detection combinations	a	b	c	$λ / m$	$β_{Δ N}$
(1, −4, 3)	0.0157	0.0086	0.9757	9.768	4.210
(−4, 6, −1)	0.0157	0.0086	0.9757	2.254	26.686

ISF: ionospheric scale factor.

Figure 6 shows P_S and P_M of the second and third detection combinations.

Figure 6.

P_S and P_M of the second and third detection combinations. (a) The success detection rate. (b) The missed detection rate.

According to Figure 6, it is known that the success detection rate of the second detection combination (1, −4, 3) is close to 100%, and the missed detection rate is close to 0%. The second detection combination benefits from longer wavelength and smaller ISF, which are minimally affected by the code measurement noise and ionosphere. Due to shorter wavelength and larger ISF, the combined signal (−4, 6, −1) is more vulnerable to the ionospheric bias. When $Δ I = 0.03 m$ , P_S of (−4, 6, −1) will drop below 94% for all possible elevation angles. Since the missed detection rate is only related to $σ_{Δ N_{(i, j, k)}}$ , it can maintain a low level regardless of the ionospheric bias. To eliminate the ionospheric bias of (−4, 6, −1) and improve the success detection rate, the ionospheric bias must be compensated. The detection value of the detection combination (−4, 6, −1) corrected for ionospheric bias can be expressed as

Δ {\tilde{N}}_{(- 4, 6, - 1)} = Δ N_{(- 4, 6, - 1)} - β_{Δ N_{(- 4, 6, - 1)}} Δ I / λ_{(- 4, 6, - 1)}

Cycle slips of three phase observations can be derived as

[\begin{matrix} Δ {\hat{N}}_{1} \\ Δ {\hat{N}}_{2} \\ Δ {\hat{N}}_{3} \end{matrix}] = {[\begin{matrix} 0 & 1 & - 1 \\ 1 & - 4 & 3 \\ - 4 & 6 & - 1 \end{matrix}]}^{- 1} {[\begin{matrix} Δ {\hat{N}}_{(0, 1, - 1)} & Δ {\hat{N}}_{(1, - 4, 3)} & Δ {\hat{N}}_{(- 4, 6, - 1)} \end{matrix}]}^{T}

Time-differenced ionospheric delay compensation algorithm based on moving windows

The precise value of $Δ {\tilde{N}}_{(- 4, 6, - 1)}$ is expected to obtain, and how to get the accurate value of $Δ I$ in current epoch is the key issue. If there is no cycle slip on $Φ_{1}$ and $Φ_{3}$ , $Δ I$ of current epoch can be used to update equation (24), and the cycle slip of the third combination will be determined correctly. When cycle slips occur on $Φ_{1}$ or $Φ_{3}$ , $Δ I$ will be biased by the cycle slips. If $Δ I$ of the current epoch is used to update equation (24), it is likely to cause misjudgment of cycle slip. However, the first and second detection combinations are insensitive to cycle slip combination $Δ N_{1} = Δ N_{2} = Δ N_{3}$ ; if identical cycle slips occur on three carrier phases, the first and second combinations cannot detect cycle slips and $Δ I$ may not be estimated correctly. To determine if there are identical cycle slips occur on the three carrier phases, the second-order time-differenced ionospheric delay $Δ Δ I$ is introduced. When there is no cycle slip, $Δ Δ I$ is very small. If identical cycle slips $Δ N_{1} = Δ N_{2} = Δ N_{3}$ occur on the three carrier phase observations, the value of $Δ Δ I$ will become abnormal.

First moving window

The first moving window of $Δ Δ I$ is used to check whether $Δ Δ I$ has a significant increase. When there is no cycle slip, we assume that $Δ Δ I$ is normally distributed with zero mean. $W F σ_{Δ Δ I}$ is the standard deviation of $Δ Δ I$ in the moving window. Cycle slips including $Δ N_{1} = Δ N_{2} = Δ N_{3}$ may occur when $Δ Δ I_{c u r}$ satisfies the following condition:

| Δ Δ I_{c u r} | > 4 \cdot W F σ_{Δ Δ I}

where $Δ Δ I_{c u r}$ is $Δ Δ I$ of current epoch, and $4 \cdot W F σ_{Δ Δ I}$ is adopted as the threshold. The number of moving windows is set to 20. If $Δ Δ I_{c u r}$ satisfies equation (26) or $Δ {\hat{N}}_{(0, 1, - 1)} \neq 0$ or $Δ {\hat{N}}_{(1, - 4, 3)} \neq 0$ , equation (24) cannot be updated by $Δ I$ of current epoch, and the second moving window is adopted to evaluate $Δ I$ of current epoch.

Second moving window

To evaluate $Δ I$ of the current epoch, the second moving window of $Δ I$ is adopted. The average value of $Δ I$ in the second moving window is used to evaluate $Δ I$ of current epoch. The mean of $Δ I$ in the moving window can be expressed as

W S_{Δ \bar{I}} = \frac{1}{n} \sum_{k = 1}^{n} Δ I_{k}

where $W S_{Δ \bar{I}}$ is the average value of $Δ I$ in the moving window, and the number of moving window is $n = 10$ . When $| Δ Δ I_{c u r} | > 4 \cdot W F σ_{Δ Δ I}$ or $Δ {\hat{N}}_{(0, 1, - 1)} \neq 0$ or $Δ {\hat{N}}_{(1, - 4, 3)} \neq 0$ , $W S_{Δ \bar{I}}$ is used as $Δ I$ of current epoch to update equation (24).

Process of compensating for ionospheric bias

Figure 7 shows that the process of the third combination compensates for the ionospheric bias. When $Δ {\hat{N}}_{(0, 1, - 1)} \neq 0$ or $Δ {\hat{N}}_{(1, 4, - 3)} \neq 0$ , we can be sure that there must be a cycle slip due to success detection rates of the first and second detection combinations which are close to 100%. However, the first and second detection combinations are insensitive to cycle slip combination $Δ N_{1} = Δ N_{2} = Δ N_{3}$ ; hence, the first moving window is introduced to detect identical cycle slips. When $| Δ Δ I_{c u r} | > 4 \cdot W F σ_{Δ Δ I}$ , cycle slips including identical cycle slips may occur, and equation (24) will be updated by $W S_{Δ \bar{I}}$ . If all three conditions in Figure 7 are not satisfied, equation (24) will be updated by $Δ I_{c u r}$ .

Figure 7.

The third combination compensates for the ionospheric bias.

Figure 8 shows $Δ I$ and $Δ Δ I$ of C29 from MIZU station. From Figure 8(a), we can see that the threshold $Δ Δ I + 4 W F σ_{Δ Δ I}$ and $Δ Δ I - 4 W F σ_{Δ Δ I}$ of the first moving window is constantly changing with $Δ Δ I$ , and a cycle slip may occur when $Δ Δ I_{c u r}$ exceeds the threshold. Figure 8(b) shows the evaluation value of the second moving window for $Δ I$ of the current epoch. It can be seen that $Δ I$ and $W S_{Δ \bar{I}}$ have the same trend. Except few epochs, the difference between $Δ I$ and $W S_{Δ \bar{I}}$ is small. From Figure 8(a) and (b), we can know that as the elevation angle of satellite decreases, $| Δ I |$ and $| Δ Δ I |$ will increase, and the corresponding threshold of the first moving window will also increase.

Figure 8.

$Δ I$ and $Δ Δ I$ of C29 from MIZU station. (a) $Δ Δ I$ , $Δ Δ I + 4 W F σ_{Δ Δ I}$ , $Δ Δ I - 4 W F σ_{Δ Δ I}$ , and elevation angle of C29. (b) $Δ I$ , $W S_{Δ \bar{I}}$ , and elevation angle of C29.

Numerical tests and analysis

To verify the effectiveness of our proposed approach, we test it in three different scenarios.

Test with simulated noise of code observations

In a poor observation environment, the noise of code observations will increase significantly. Due to the existence of code observation in the detection combination, the noise of code observations is bound to affect the detection of cycle slip. To test the ability of the proposed method to detect cycle slip in the case of high code observation noise, the noise is added to the original code observations. Gaussian noise with zero mean is added to the three code observations in the second experiment. We simulated two scenarios, in which the standard deviation of the added noise is 0.5 and 0.8 m, respectively. The data of BDS-3 satellites from SGOC, ULAB, and STUM stations on June 19, 2019 are used and the data are clean without cycle slip. The sampling interval is 30 s. At every 10 epochs, cycle slip groups from (0, 0, 0) to (2, 2, 2) will be added to the data in turn. The results of SGOC station are shown in Figure 9. Table 4 shows the total number of simulated cycle slips and the number of incorrect detection and repair.

Figure 9.

The times of failed detection and correction for SGOC station. (a) Test using noise with standard deviation of 0.5 m. (b) Test using noise with standard deviation of 0.8 m. (c) Test under normal situation.

Table 4.

Cycle slips detection results of SGOC, SUTM, and ULAB stations.

Name of station	Total number of cycle slips	Success rate
Name of station	Total number of cycle slips	Scenario #1 (%)	Scenario #2 (%)	Scenario #3 (%)
SGOC	1595	99.56	98.75	95.11
SUTM	1607	100	99.50	95.15
ULAB	1410	99.36	98.30	96.45

In Table 4, scenario #1 is normal and no noise is added. Simulated noise is added to scenario #2 with standard deviation of 0.5 m. In scenario #3, simulated noise is added using noise with standard deviation of 0.8 m. Scenario #3 has higher noise level than scenario #2.

From Figure 9 and Table 4, we can know that the number of failures increases as the noise level increases, but the success rate is still greater than 95%. From Table 3, it is known that the code observations have a very small weight in the detection combination; hence, the proposed approach is less affected by code observations.

Test under normal circumstances

The second experiment was carried out under normal circumstances, the noise of code observations is at normal level, and the ionosphere changes slowly. We compare the proposed approach with the approach proposed by Zhao et al.³² They used three code-phase combinations to detect and correct cycle slip, and in order to eliminate the ionospheric bias in the third combination, they propose

Δ {\tilde{N}}_{(i, j, k)} = Δ N_{(i, j, k)} - κ \cdot β_{Δ N_{(i, j, k)}} Δ I / λ_{(i, j, k)}

where $κ (0 \leq κ \leq 1)$ is scale factor to balance the corrected percentage of the ionospheric bias and the amplification of the measurement noise. When $κ = 1$ , the ionospheric bias is totally corrected. When $κ = 0$ , the ionospheric bias is not corrected. When $0 \leq κ \leq 1$ , the ionospheric bias is partly corrected.

The real BDS-3 data from SGOC station on June 19, 2019 are used to test the effectiveness of these two methods. The original observations have been repaired and have no cycle slips beforehand. Cycle slip groups from (0, 0, 0) to (2, 2, 2) are added to the data in turn, at every 10 epochs. The sampling interval is 30 s. The approach of Zhao et al.³² is tested when $κ = 0, 0.5, 1$ , respectively. The results of C20–C36 satellites can be seen in Figure 10, in which C31 satellite of SGOC station has no data.

Figure 10.

The times of incorrect cycle slip detection and repair. (a) Test using the method proposed by Zhao et al.³² when $κ =0$ . (b) Test using the method proposed by Zhao et al.³² when $κ = 0.5$ . (c) Test using the method proposed by Zhao et al.³² when $κ =1$ .

In Figure 10, the graph consisting of two red triangles and a green square represents failure of detection and repair. The number in the right coordinate area represents the times of failed detection and correction. As shown in Figure 10, the times of failed detection and correction in Figure 9(a) have a significant reduction compared with Figure 10(a) to (c). Therefore, the method proposed by this article has stronger ability to eliminate the effect of ionospheric bias. From Figure 10, it can be known that all of the failed detection and repair of cycle slip occurs at low elevation angle. At low elevation angles, the ionospheric bias and code measurement noise may increase. If the ionospheric bias cannot be handled correctly, it will lead to erroneous results.

Test under high ionospheric activity

In theory, as long as the prediction of $Δ I$ under high ionospheric activity can reach a required precision, the cycle slip will be detected. Many methods use STPIR algorithm to perform the second-order time-difference on detection combination, which will increase the noise of detection combination. To verify the validity of the approach, the data of BDS-3 during the period of active ionosphere will be used. The observations of BDS-3 are collected from WUH2 station on March 21, 2019, and the ionosphere is active for a period of time. The sampling interval is 30 s. Figure 11 shows $Δ I$ and $Δ Δ I$ of C28. In Figure 11(a), the ionospheric delay variation of C28 can reach 0.18 m. The prediction values $W S_{Δ \bar{I}}$ for ionospheric delay variation are shown in Figure 11(b). Different types of cycle slips are added to original observations, and the results are shown in Figure 12 and Table 5.

Figure 11.

(a) $Δ I$ and $Δ Δ I$ of C28 from WUH2 station. (b) $Δ I$ and $W S_{Δ \bar{I}}$ of C28 from WUH2 station.

Figure 12.

Cycle slips detection of C28 with simulated cycle slips.

Table 5.

Added cycle slips epoch, type, and detection results.

Epoch	Simulated cycle slip	Cycle slip type	$(Δ N_{(0, 1, - 1)}, Δ N_{(1, - 4, 3)}, Δ {\tilde{N}}_{(- 4, 6, - 1)})$	Fixed cycle slip
20	(1, 0, 1)	S	(−0.883, 3.81, −4.74)	(1, 0, 1)
35	(11, 13, 20)	L	(−6.85, 19, 14.1)	(11, 13, 20)
80	(1, 0, 0)	S&P	(−0.020, 0.981, −3.76)	(1, 0, 0)
90	(1, 1, 1)	S&P	(0.089, 0.080, 1.16)	(1, 1, 1)
100	(0, 0, 1)	S	(−1, 3.08, −1.14)	(0, 0, 1)
110	(1, 1, 0)	S	(1.01, −2.95, 1.84)	(1, 1, 0)
120	(0, 1, 0)	S	(0.95, −3.99, 5.86)	(0, 1, 0)
130	(0, 1, 1)	S&P	(0.015, −0.997, 4.94)	(0, 1, 1)
140	(5, 8, 10)	L	(−1.94, 2.98, 17.9)	(5, 8, 10)

Small, particular, and large cycle slips are simulated in the third experiment. In Figure 12, $β_{Δ N_{(- 4, 6, - 1)}} Δ I / λ_{(- 4, 6, - 1)}$ denotes the deviation of the cycle slip detection caused by the ionospheric delay variation, and the unit is cycle. From Figure 12, we can see that $| β_{Δ N_{(- 4, 6, - 1)}} Δ I / λ_{(- 4, 6, - 1)} |$ is more than 0.5 cycles in many epochs. The bias caused by ionospheric delay can reach two cycles in 80th and 90th epochs. Although the ionosphere has a significant impact on the third detection combination, all simulated cycle slips can be correctly corrected.

From Figure 12 and Table 5, we can know that cycle slips (1, 0, 0), (1, 1, 1), and (0, 1, 1) are insensitive for the first detection combination (0, 1, −1). The second detection combination cannot detect cycle slip (1, 1, 1). However, all cycle slips can be detected with the combined three combinations.

Conclusions

A new method of cycle slip detection and repair for BDS-3 triple-frequency observations is proposed. Firstly, two modified code-phase combinations are used to minimize the effect of noise of code measurement. Then, 1% missed detection rate and 99% success detection rate are used as the criterion to optimally select the second and third detection combinations by combining BDS elevation-based model and HMW combination. Since the third detection combination is easily biased by ionosphere, two moving windows are used to correct the ionospheric bias of the third combination. The method has been tested in three different scenarios. In the first scenario, simulated Gaussian noise is added to the three code observations and the proposed method is still effective in the environment of high code observation noise. Compared with the method proposed by Zhao et al.,³² the proposed method has obvious advantages in reducing the times of failed detection and correction. From the first and second experiments, we can know that most of the detection failures occur at low elevation angles. We can increase the cutoff elevation angle of satellites to reduce the times of failed detection and repair. Further, BDS-3 data under high ionospheric activity is used to verify the validity of the algorithm, cycle slips of different types are added to the original carrier phase observations, and all cycle slips can be correctly detected and fixed.

Footnotes

Acknowledgement

The authors would like to thank the National Natural Science Foundation of China.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (no. 41574010).

ORCID iD

Xiangxiang Fan

Appendix 1

All symbols and their definitions used in this article are shown in the following table.

References

Bian

Cao

, et al. GNSS spoofing detection based on new signal quality assessment model. GPS Solut 2017; 22(1): 28.

Liu

Cai

Wang

. Cooperative localization of connected vehicles: integrating GNSS with DSRC using a robust cubature Kalman filter. IEEE. T Intell Transp 2017; 1–15.

Gao

Zhuang

, et al. Robust Kalman filter aided GEO/IGSO/GPS raw-PPP/INS tight integration. Sensors 2019; 19(2): 417.

Duan

Sun

Ouyang

, et al. Reducing the effect of positioning errors on kinematic raw doppler (RD) velocity estimation using BDS-2 precise point positioning. Sensors 2019; 19(13): 3029.

Zheng

Li , et al. Improving the GNSS-R specular reflection point positioning accuracy using the gravity field normal projection reflection reference surface combination correction method. Remote Sensing 2018; 11(1): 33.

Jin

. Triple-frequency carrier phase precise time and frequency transfer models for BDS-3. GPS Solut 2019; 23(3): 1–12.

Montenbruck

Hauschild

Steigenberger

, et al. Initial assessment of the COMPASS/BeiDou-2 regional navigation satellite system. GPS Solut 2012; 17(2): 211–222.

Beutler

Davidson

Langley

, et al. Some theoretical and practical aspects of geodetic positioning using carrier phase difference observations of GPS satellites. In: Surveying Engineering, University of New Brunswick: Fredericton, NB, Canada, 1985.

Hofmann

Lichtenegger

Collins

. GPS Theory and Practice. Springer Wien: New York, NY, USA, 2001.

10.

Zangeneh-Nejad

Amiri-Simkooei

Sharifi

, et al. Cycle slip detection and repair of undifferenced single-frequency GPS carrier phase observations. GPS Solut 2017; 21: 1593–1603.

11.

Momoh

Bhattarai

Ziebart

. Receiver clock jump and cycle slip correction algorithm for single-frequency GNSS receivers. GPS Solut 2019; 23(2): 1–12.

12.

Hatch

. The synergism of GPS code and carrier measurements. In: Proceedings of the International Geodetic Symposium on Satellite Doppler Positioning, Las Cruces, NM, USA, 8–12 February 1982; pp. 1213–1231.

13.

Melbourne

. The case for ranging in GPS-based geodetic systems. In: Proceedings of the First International Symposium on Precise Positioning with the Global Positioning System, Rockville, MD, USA, 15–19 April 1985.

14.

Wübbena

. Software developments for geodetic positioning with GPS using TI-4100 code and carrier measurements. In: Proceedings of the First International Symposium on Precise Positioning with the Global Positioning System, Rockville, MD, USA, 15–19 April 1985, vol. 19.

15.

Blewitt

. An automatic editing algorithm for GPS data. Geophys Res Lett 1990; 17: 199–202.

16.

Gao

, et al. step cycle slip detection and repair method based on double-constraint of ephemeris and smoothed pseudorange. ACTA Geody ET Geom 2019; 16(4): 337–348.

17.

Liu

. A new automated cycle slip detection and repair method for a single dual-frequency GPS receiver. J Geodesy 2010; 85: 171–183.

18.

Cai

Liu

Xia

, et al. Cycle slip detection and repair for undifferenced GPS observations under high ionospheric activity. GPS Solut 2012; 17: 247–260.

19.

Banville

Langley

. Mitigating the impact of ionospheric cycle slips in GNSS observations. J Geodesy 2013; 87: 179–193.

20.

Wang

. New-generation BeiDou (BDS-3) experimental satellite precise orbit determination with an improved cycle-slip detection and repair algorithm. Sensors 2018; 18(5): 1402.

21.

Feng

. GNSS three carrier ambiguity resolution using ionosphere-reduced virtual signals. J Geodesy 2008; 82: 847–862.

22.

Zhang

. Benefits of the third frequency signal on cycle slip correction. GPS Solut 2016; 20: 451–460.

23.

Dai

Knedlik

Loffeld

. Real-time cycle-slip detection and determination for multiple frequency GNSS. In: Proceedings of the 2008 5th Workshop on Positioning, Navigation and Communication (WPNC 2008), Hannover, Germany, 27 March 2008, pp. 37–43.

24.

De Lacy

Reguzzoni

Sansò

. Real-time cycle slip detection in triple-frequency GNSS. GPS Solut 2011; 16: 353–362.

25.

Zhao

Sun

Dai

, et al. Real-time detection and repair of cycle slips in triple-frequency GNSS measurements. GPS Solut 2014; 19: 381–391.

26.

Zhu

. Detection and correction of cycle slip in triple-frequency GNSS positioning. IEEE Acc 2017; 5: 12584–12595.

27.

Huang

Zhai

, et al. A new triple-frequency cycle slip detecting algorithm validated with BDS data. GPS Solut 2015; 20: 761–769.

28.

Teunissen

. The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation. J Geod 1995; 70: 65–82.

29.

Xiao

Mayer

Heck

, et al. Improved time-differenced cycle slip detect and repair for GNSS undifferenced observations. GPS Solut 2018; 22(1): 1–13.

30.

Zeng

Sui

, et al. Real-time triple-frequency cycle slip detection and repair method under ionospheric disturbance validated with BDS data. GPS Solut 2018: 22(3): 1–13.

31.

Melachroinos

. An enhanced cycle slip repair algorithm for real-time multi-GNSS, multi-frequency data processing. GPS Solut 2018; 23(1): 1–11.

32.

Zhao

Roberts

Hancock

, et al. A triple-frequency cycle slip detection and correction method based on modified HMW combinations applied on GPS and BDS. GPS Solut 2019; 23(1): 1–13.

33.

Xiao

Zhiqiang

Ying

, et al. Real-time cycle slip correction for a single triple-frequency BDS receiver based on ionosphere-reduced virtual signals. Adv Space Res 2018; 62: 2381–2392.

34.

Yao

Gao

Wang

, et al. Real-time cycle-slip detection and repair for BeiDou triple-frequency undifferenced observations. Surv Rev 2016; 48: 367–375.

35.

Qin

, et al. Undifferenced cycle slip estimation of triple-frequency BeiDou signals with ionosphere prediction. Mar Geodesy 2016; 39: 348–365.

36.

Xiong

. An improved algorithm based on combination observations for real time cycle slip processing in triple frequency BDS measurements. Adv Space Res 2019; 63: 2796–2808.

37.

Liu

Jin

, et al. A new real-time cycle slip detection and repair method under high ionospheric activity for a triple-frequency GPS/BDS receiver. Sensors 2018; 18(2): 427.

38.

. Stochastic modeling of triple-frequency BeiDou signals: estimation, assessment and impact analysis. J Geodesy 2016; 90: 593–610.

39.

Eueler

Goad

. On optimal filtering of GPS dual frequency observations without using orbit information. Bull Geodesique 1991; 65: 130–143.

40.

Wang

De Jong

Zhao

, et al. Multipath analysis of code measurements for BeiDou geostationary satellites. GPS Solut 2015; 19(1): 129–139.

Cycle slip detection and repair for BeiDou-3 triple-frequency signals

Abstract

Keywords

Introduction

Triple-frequency linear combination

Error analysis

Satellite elevation-based measurement noise model of code observations

Time-differenced ionospheric delay

Cycle slip detection combinations

First detection combination

Second and third detection combinations

Determination of the second and third detection combinations

Success detection rate and missed detection rate

Selection of optimal detection combinations

Time-differenced ionospheric delay compensation algorithm based on moving windows

First moving window

Second moving window

Process of compensating for ionospheric bias

Numerical tests and analysis

Test with simulated noise of code observations

Test under normal circumstances

Test under high ionospheric activity

Conclusions

Footnotes

Acknowledgement

Declaration of conflicting interests

Funding

ORCID iD

Appendix 1

References