Abstract
The traditional triple-frequency geometry-free pseudorange minus phase (GFPMP) method is very susceptible to the influence of pseudorange observation noise, and it is difficult to detect insensitive small cycle slips. The dual-frequency phase ionospheric residual (PIR) method has the ability to detect sensitive small cycle slips, but its detection results have multivalue problems. In view of this, a GFPMP-PIR combination model method with wavelet denoising (GPW) is proposed here. The idea of this method is as follows: firstly, wavelet transform is used to denoise the pseudorange observations; then, by optimizing and selecting the combination coefficients with the minimum condition, an optimized GFPMP-PIR combination model is constructed. Hence, the source data quality is assured through wavelet denoising. And the accuracy of cycle slip detection and repair is improved by using the dual-frequency PIR combination. Moreover, the multivalue problem of the dual-frequency PIR method is solved by triple-frequency GFPMP combination; finally, the BeiDou navigation satellite system (BDS) triple-frequency measurement data is used for experimental verification. The experiment results show that the GFPMP-PIR optimization combination model with wavelet denoising can detect and repair various cycle slips, especially for insensitive small cycle slips.
1. Introduction
With the development of global navigation satellite system (GNSS), the research on the integrity, continuity, and accuracy of high-precision positioning becomes more and more important. In GNSS precision positioning, when tracking the carrier phase of a GNSS signal, if phase lock is lost during a specific observation interval, the resulting observation may show a discontinuity compared with observations made before and after. This discontinuity is referred to as small cycle slips [1], which seriously affects the accuracy and reliability of navigation and positioning [2]. It is a valid way to detect and repair the cycle slips for the realization of high-precision positioning of the system [3–5].Up to now, many methods with respect to cycle slip detection and repair had been studied and developed. The polynomial fitting method is more suitable for cycle slip detection of static GNSS data. But it is susceptible to the randomness of the observation noise, so the detection of small cycle slips is \not ideal, and when the cycle slips occur multiple times in multiple epochs, the detection effect will also deteriorate [6–8]. Although the wavelet transform method is widely used in various fields, the wavelet coefficients decomposed when detecting cycle slips cannot accurately represent the size of cycle slips, and the efficiency of cycle slip repair is low [9, 10]. The PIR combination method can detect small cycle slips, but it has the problem of multivalue in detection results, which makes it undefinable to accurately detect and repair cycle slips [11, 12]. The Melbourne–Wübeena (MW) combination method has real-time functionality, but it cannot detect the cycle slip of the combinations with equivalent coefficients [13, 14]. In order to overcome the disadvantage of the MW combination method, the TurboEdit method is proposed [15–17]. It combines MW combination method with geometry-free (GF) combination method together and is currently widely used, especially for cycle slip detection of nondifferenced data. However, due to the observation noise of pseudorange, it cannot detect the small cycle slips.
The above methods are mainly applicable to single- and dual-frequency data. With the increase of the frequency number of carrier observation, it has become a trend to detect and repair cycle slip by combining three or even multiple-frequency data with the reasonable combination coefficients. The geometry-free pseudorange minus phase (GFPMP) [18, 19] combinations and the geometry-free (GF) [20, 21] are the typical triple-frequency methods to detect and repair cycle slip. The GFPMP combination method can quickly detect cycle slips, but its three combination coefficients are complicated to determine. Moreover, it is difficult to detect the insensitive small cycle slips for the pseudorange observation noise influence. Combined with GF method, the combination coefficients of GFPMP become to easily select with minimum condition. But the condition number of the combination coefficient matrix is large, which may cause the ill-conditioned equation problem to occur. The extrawide lane, wide lane, and narrow lane combination methods can determine the cycle slips sequentially in three cascaded steps [22, 23]. However, these methods are performed well, and the ionosphere influence is reduced well. Moreover, they are also severely affected by pseudorange observation noise.
In view of the shortcomings of the above methods, in this paper, firstly, we were using wavelet denoising method to reduce the influence of pseudorange observation noise with respect to the triple-frequency GFPMP combination, then combining triple-frequency GFPMP combination with dual-frequency PIR combination to construct a GFPMP-PIR combination model, then selecting the optimal combination coefficients to solve the problem of the insensitive cycle slip detection of the GFPMP combination, and overcoming the multivalue of the detection results of the PIR combination; finally, experimental verification is carried out through BDS triple-frequency measurement data.
2. GFPMP-PIR Combination Model with Wavelet Denoising to Detect and Repair Cycle Slips
2.1. GFPMP-PIR Combination Model with Wavelet Denoising to Detect Cycle Slip
The BDS carrier phase and pseudorange observation equations can be derived as follows:
In Equations (1) and (2), the subscript denotes multiple the frequencies of BDS; the symbols , , , and are the carrier phase measurements, integer ambiguities, carrier phase multipath errors, and carrier observation noise errors, respectively; the symbols , , , and denote the carrier phase wavelengths, pseudorange measurements, pseudorange multipath errors, and pseudorange observation noise errors, respectively; the symbol denotes the geometric distance between the receiver antenna and the satellite; the symbols and denote the receiver clock error and the satellite clock error, respectively; the symbol denotes the ionospheric delay with respect to the frequency B1I of BDS; the symbol denotes the tropospheric delay error.
When using the GFPMP combination method to detect cycle slip, it is mainly affected by pseudorange observation noise, carrier phase observation noise, multipath error, ionospheric delay, and tropospheric delay [19]. In the case of a high sampling rate (1 Hz), the multipath error, ionospheric delay, and tropospheric delay of the adjacent epoch have a strong temporal and spatial correlation. Thus, the pseudorange phase combination can exactly eliminate the influence of these errors through the difference between epochs [24–26]. And due to the carrier phase, observation noise error is much smaller than the pseudorange observation noise error; therefore, after the epoch difference, the GFPMP combination is mainly affected by the pseudorange observation noise. Hence, before constructing the GFPMP-PIR combination model, the pseudorange observation errors must to be denoised.
Wavelet transform has the advantage in multiscale analysis of signal local characteristics. The db4 wavelet transform of the Daubechies wavelet system can not only eliminate noise but also keep the position of the signal’s step or mutation point unchanged. Moreover, it can ensure that most of the original important information and characters of the signal are not lost. Therefore, the db4 wavelet transform can be used to denoise the GNSS pseudorange observation errors to improve the performance of cycle slip detection.
We assume that the pseudorange observation is , let , and use db4 wavelet transform to decompose at different resolutions. The decomposition can be expressed as follows: where the subscript denotes the number of decomposition (, and the symbol denotes epoch number of the observation). The subscript denotes decomposition scale. The symbol denotes the scale coefficient. The symbol denotes the detail coefficient. The symbol denotes the filter length. The symbols and denote the high- and low-frequency filters, respectively. There is no fixed selection formula for the scale of wavelet. Generally, different decomposition scales are selected according to different processed signals. Considering the characteristics and computational efficiency of pseudorange observation signal, one-layer wavelet decomposition scale is adopted in this paper, that is, .
When the original signal with noise is processed by using wavelet transform, the maximum norm of the wavelet coefficients with respect to the useful signal will be increase with the scale is increasing, while the maximum norm of the wavelet coefficients with respect to the noise signal will be smaller with the scale is increasing. Thus, the wavelet coefficients of the useful and the noise signals are clearly distinguished by wavelet transformation. Based on this principle, the thresholded of the high-frequency coefficients can be set reasonably.
The selection and processing strategy about the threshold are the critical steps to denoise by wavelet transform. Visushrink, SureShrink, Minimax, and BayesShrink threshold selections and processing strategies are the common methods [27]. This paper adopts the Visushrink threshold, which can denoise thoroughly, and it can be described as follows: where the symbol denotes the middle number of the absolute of the high-frequency coefficients. The symbol denotes signal length.
Hard threshold function and soft threshold function are commonly used to denoise in data processing. The soft threshold function has good continuity and smoothness characters, but all the coefficients are compressed, which will cause some effective information of the signal to be lost. However, the hard threshold function does not compress the step or mutation points of the signal. So, the hard threshold processing strategy is used widely to denoise. The hard threshold function can be expressed as follows:
After the threshold is determined by the Visushrink threshold method, the high-frequency coefficient can be got by Equation (5). Further, the denoised pseudorange observation can be described as follows: where denotes the denoised pseudorange observations of B1I, B2I, and B3I of BDS; and are the conjugate of and , respectively.
When the pseudorange observation noise errors are denoised, according to the BDS dual-frequency and triple-frequency combination theories [19, 28], the triple-frequency pseudorange GFPMP combination and dual-frequency PIR combination can be derived as follows:
In Equations (7) and (8), the symbols , , and denote the carrier phase combination coefficients of B1I, B2I, and B3I of BDS, respectively; the symbols , , and denote the pseudorange combination coefficients with respect to B1I, B2I, and B3I, and their sum is 1; the symbols and denote the ionospheric delay coefficients of GFPMP combination and PIR combination, respectively; the symbol denotes the wavelength of the combined observation; the subscript symbols and denote the B1I, B2I, or B3I of BDS.
The higher the sample rate is, such as 1 Hz or smaller, the smaller of the ionospheric delay variation of the adjacent epochs is. Under this condition, the ionospheric delay influence of the adjacent epochs can be ignored [29]. According to Equations (7) and (8), we can obtain the difference equations of the adjacent epochs as follows:
In Equations (9) and (10), the symbol denotes the difference between adjacent epochs; the symbols and denote the difference GFPMP and PIR combinations of triple-frequency BDS observations, respectively, which can be defined as the cycle slip detection observations. Assuming that the carrier phase observations and pseudorange observations with respect to the three frequencies of BDS are independent each other, according to literature [28], we can set the noise error of the carrier phase cycle, and set the noise error of pseudorange observation, which has been denoised, . Based on the law of error propagation, we can get the errors of and as follows:
When the absolute value of the cycle slip detection observations is correspondingly greater than 4 times of their mean square errors (4 times of the mean square error is considered the detection threshold [30]), the confidence level is 99.9%, and cycle slip can be considered to have occurred. According to Equation (11), to minimize the influence of pseudorange observation noise, the coefficients of pseudorange combination are set to 1/3, that is, .
2.2. GFPMP-PIR Combination Model with Wavelet Denoising to Repair Cycle Slip
In order to better detect and repair insensitive cycle slip, referring to literature [31], three linearly independent combinations, including two triple-frequency GFPMP combinations and one dual-frequency PIR combination, are selected. Based on reference [30], the carrier phase coefficients of the two GFPMP combinations are selected as (0, 1, -1) and (1, 4, -5) with respect to the B1I, B2I, or B3I of BDS.
There are three ways to combine the dual-frequency PIR combination with triple-frequency GFPMP combinations. The combination coefficient matrix can be expressed as follows: where , , and denote the coefficients of PIR combination. When B1I and B2I are selected, the values of , , and are 1, , and 0, respectively; when B1I and B3I are selected, the values of , , and are 1, 0, and , respectively; when B2I and B3I are selected, , , and are 0, 1, and , respectively. In order to repair cycle slip validly, it is necessary to consider the condition number of the coefficient matrix . If the number is too large, the equation will be ill conditioned and even cause the instability of the solution result. The condition number of the coefficient matrix formed by different combinations is shown in Table 1.
It can be seen from Table 1 that the noise errors of the three PIR combinations are all similar, and the condition number of the combination coefficient matrix composed of B1I and B3I frequencies is the smallest one. Therefore, we choose , , and , as the coefficients of the PIR combination.
When the optimal combination coefficients are selected, the coefficient matrix of the GFPMP-PIR combination is , is the matrix formed of carrier phase observations with cycle slips, is the matrix formed of the cycle slip detection observation with wavelet denoising, and the relationship of them can be expressed as follows: where , , and .
The condition number of the matrix is 227.41. It is smaller than 414.24, which used in reference [30]. Hence, the matrix obtained by the combination model in this paper is more stable. Moreover, the rank of the matrix is 3, so Equation (14) has a unique solution, and the multivalue problem of the detection results is also solved.
Theoretically, according to Equation (14), the cycle slip estimate can be obtained directly from the float solution . Due to the influence of observation noise, sometimes this will cause cycle slip repair to fail. In order to ensure the success rate of cycle slip repair, on the basis of the float solution to determine the initial value of the cycle slip repair value, the space search and the 1-norm minimum principle are used to search and confirm it further. Base on reference [32], the confirmation formula is given as follows:
3. Results and Discussion
This paper is verified by using BDS triple-frequency measurement data. The data acquisition time is May 23, 2021. The acquisition location is the school of surveying and mapping of Henan Polytechnic University. The data sampling rate is 1 Hz, a total of 7200 epochs. The observation data of Beidou C1 satellite are selected for experimental analysis. In this paper, two experimental schemes are designed to verify the correctness, effectiveness, and applicability of GPW.
Scheme I. Scheme I is designed to verify by the carrier phase measurement and pseudorange data. Comparing with the GFPMP combination method in reference [30], which is referred to as GPGF hereinafter, the cycle slip detection and repair result of the basic performance with respect to GPW are verified. For getting the carrier phase measurement in certain period, the type character of the cycle slip is single. In order to further verify the detection applicability of the GPW for various cycle slips, according to scheme I, the carrier phase data without cycle slip is obtained and used as the experimental data in scheme II.
Scheme II. As mentioned above, the carrier phase data processed in scheme I is used as the phase detection sequence. And the original pseudorange observation data is used as the other experimental data here. Add artificial cycle slips to the carrier phase sequence without cycle slip. The specific method is as follows: starting from the fifth epoch, the cycle slip combinations (1, 0, 0), (2, 0, 0), ..., (9, 9, 9) are added every 5 epochs. In order to better reflect the diversity of cycle slips, special artificial cycle slip combinations (0, 10, 0), (10, 10,10), (6, 11, 10), (12, 0, 10), (13, 0, 10) are added every 400 epochs from the 5500th epoch, and 1004 artificial cycle slip combinations are added totally. GPW is used to detect and repair the carrier phase data with artificial cycle slips. And the detection and repair results are compared with those of GPGF. The purpose is to fully verify the detection and repair correctness, effectiveness, and applicability of GPW for different cycle slip combinations. Moreover, the effect of wavelet denoising is also verified.
3.1. Scheme I Results and Discussion
Figure 1 shows the cycle slip detection results using three optimized coefficient combinations of GPW. It can be seen from Figure 1 that no cycle slip has been detected by the combination coefficient (0, -1, 1), the cycle slip at the 1518th epoch has been detected by the combination coefficient (1, 4, -5), and three cycle slips at the 1514th, 1517th, and 1518th epochs have been detected by the combination coefficient (1, 0, -1.2306). Therefore, the GPW can detect cycle slips in the phase detection sequence through the optimized coefficient combinations.

(a)

(b)

(c)
Figure 2 shows the cycle slip detection results using the three optimized coefficient combinations of GPGF. As shown in Figure 2, the detection result of GPGF is consistent with those of the GPW. Cycle slips are detected in the 1514th, 1517th, and 1518th epochs. Therefore, the cycle slip detection capability of the two methods is equivalent.

(a)

(b)

(c)
According to Equations (14) and (15), the float solution and the estimate about the cycle slip with respect to GPW and GPGF can be obtained, as shown in Table 2. The closer the float solution is to the cycle slip estimate, the higher the accuracy of the cycle slip repair result is. It can be seen from Table 2 that the two methods can correct cycle slips, but the float solution of GPW is closer to the cycle slip estimate value. Therefore, comparing with the GPGF, the GPW has a higher accuracy of the cycle slip repair. After the measurement data is repaired by the two methods, the change of the cycle slip detection observations is within the detection threshold. Due to limited space, the restored cycle slip detection diagram will no longer be displayed here.
3.2. Scheme II Results and Discussion
Through scheme I, we got the phase data without cycle slip. In order to further verify the ability of GPW to detect various cycle slips, we add artificial cycle slips according to scheme II plan above and got the experimental data for scheme II.
Figure 3 shows the artificial cycle slip detection results using three optimized coefficient combinations of GPW. According to the statistics, GPW with the combination coefficients of (0, -1, 1), (1, 4, -5) cannot detect the cycle slip combinations similar to (1, 1, 1), (2, 2, 2),..., (10, 10, 10) and can detect other types of cycle slip combinations; GPW with the combination coefficients of (1, 0, -1.2306) cannot detect the cycle slip combinations similar to (0, 1, 0), (0, 2, 0),..., (0, 10, 0) and can detect other types of cycle slip combinations. Therefore, through the combination of three optimized different coefficients, GPW can detect all kinds of cycle slips.

(a)

(b)

(c)
Figure 4 shows the artificial cycle slip detection results using three optimized coefficient combinations of GPGF. Statistical results show that the detection results of the GPGF are similar to those of the GPW. Hence, both methods can detect all artificial cycle slips. In other words, the cycle slip detection capability of the two methods is equivalent.

(a)

(b)

(c)
According to Equation (14), the float solution of the cycle slip can be obtained. To an extent, the difference between the float solution and the artificially cycle slip reflects the influence of observation noise. The greater noise has, the greater the difference is. For BDS B1I frequency of C1 satellite, 1004 artificial cycle slip combinations, which are designed above, the differences between the float solution and the artificial cycle slip by using the two methods are shown in Figure 5. It can be seen from Figure 5 that, comparing with GPGF, the float solution of GPW is correspondingly closer to the artificial cycle slip in general. The standard deviations of the difference of the two methods are shown in Table 3. It can be seen from Table 3 that the standard deviations of GPW are low. The main reason is that the wavelet threshold denoising effectively reduces the influence of pseudorange observation noise. Therefore, the cycle slip repair performance of GPW is better than that of GPGF.

4. Conclusions
Aiming at the problem of BDS triple-frequency cycle slip detection and repair, a GFPMP-PIR combination model with wavelet denoising (GPW) is constructed, two verification schemes are designed, and the measurement data of BDS is used to verify the performance of GPW. The main conclusions are summarized as follows: (1)The wavelet denoising method can reduce the influence of the pseudorange observation noise of the GFPMP combination, and the GPW can detect various cycle slips, including insensitive small cycle slips(2)Comparing with GPGF, the cycle slip repair performance of the GPW is better
In further application research, the timeliness and robustness of the GPW will be the key issues to be considered.
Data Availability
The datasets analyzed in this study are managed by the School of Surveying and Land Information Engineering, Henan Polytechnic University, and can be available on request from the corresponding author.
Conflicts of Interest
The authors declare no conflict of interest.
Authors’ Contributions
K.L. designed the experiments and wrote the main manuscript. H.B. and Y.J. reviewed the paper. S.L. edited the paper. All components of this research were carried out under K.L. All authors have read and agreed to the published version of the manuscript.
Acknowledgments
This work is funded by the National Natural Science Foundation of China (No. 41774039) and the State Key Lab Project of China (No. 6142210200104). The authors would like to thank C.K, who is a professor at The Ohio State University, for providing the valuable advice.