Abstract

Cycle slip detection and repair play important roles in the processing of data from dual-frequency GPS receivers onboard low-Earth orbit (LEO) satellites. To detect and repair cycle slips more comprehensively, an enhanced error method (EEM) is proposed. EEM combines single-frequency and narrow-lane carrier phase observations to construct special observations and observation equation groups. These special observations differ across time and satellite (ATS). ATS observations are constructed by three steps. The first step is differencing single-frequency and narrow-lane observations through a time difference (TD). The second step is to select a satellite as a reference satellite and other satellites as nonreference satellites. The third step is to difference the single-frequency TD observations from the reference satellite and the narrow-lane TD observations from the nonreference satellites by a satellite difference. If cycle slips occur at the reference satellite, the correction values for these ATS observations can be significantly enlarged. To process all satellites, the EEM selects each satellite as a reference satellite and builds the corresponding equation group. The EEM solves these observation equation groups according to the weighted least-squares adjustment (LSA) criterion and obtains the correction values; these correction values are then used to construct the values corresponding to different equation groups, and the EEM subsequently carries out a chi-square distribution test for these . The satellite corresponding to the maximum will be marked. Then, the EEM iteratively processes the other satellites. Cycle slips can be estimated by rounding the float solutions of changes in the ambiguities of cycle slip satellites to the nearest integer. The simulation test results show that the EEM can be used to detect special cycle slip pairs such as (1, 1) and (9, 7). The EEM needs only observation data in two adjacent epochs and is still applicable to observation epochs with continuous cycle slips.

1. Introduction

Since low-Earth orbit (LEO) satellites are equipped with GPS receivers, they have become the main means for supporting LEO satellite missions. In the precise orbit determination of LEO satellites, the quality of the carrier phase observation data plays an important role. However, LEO satellites are in a state of high-speed movement, and both the altitude angle and the satellite attitude can change at any time; hence, the satellite signal can easily suffer from a loss of lock, which results in a discontinuous whole cycle count of carrier phase observations, a phenomenon known as cycle slip. The cycle slip of GPS carrier phase observations is an important factor that affects data quality. Carrier phase observation data with cycle slip issues will seriously affect the orbit determination accuracy. Therefore, to improve the orbit determination accuracy, it is necessary to detect and repair cycle slips. At present, there are many methods for cycle slip detection. Existing cycle slip detection methods are divided into two categories.

The first kind of cycle slip detection method, such as the polynomial fitting method, high-order difference method, Doppler integration method [1], TurboEdit method [2], and ionospheric residual method [3, 4], does not involve observation equations. However, due to the high-speed movement of LEO satellites, these methods cannot be fully applied to cycle slip detection. In the polynomial fitting method, carrier phase observations are regarded as time-dependent variables, and the cycle slip is determined by comparing the polynomial-fitted values with the actual observations. This method is more suitable for static or low-speed GPS receivers and is insensitive to small and continuous cycle slips. Similar to the polynomial fitting method, the high-order difference method is not suitable for detecting small cycle slips. De Lacy et al. [5] combined Bayes’ theorem with the polynomial fitting method to realize cycle slip detection, but this method is not very appropriate for LEO satellites in a high-speed moving state. The Doppler integration method requires Doppler observations, so it is not suitable for receivers with only GPS observations; furthermore, affected by the accuracy of Doppler observations, the Doppler integration method is also insensitive to small cycle slips. Blewitt [2] proposed the TurboEdit method to detect and repair cycle slips together with Melbourne–Wübbena (M-W) and geometry-free (GF) combinations. The TurboEdit method needs only data from one observation epoch, so it is widely used in GPS cycle slip detection and repair. However, due to the influence of pseudorange measurement errors, it is difficult to detect single-frequency cycle slips below 2 cycles and special cycle slip pairs such as (5, 5). Joining the variation characteristics of the ionosphere with the M-W wide-lane combination, Liu [3] proposed a cycle slip detection method based on ionospheric residuals. In the detection of small cycle slips, this ionospheric residual method improves the cycle slip detection sensitivity by comparing the difference between the residual due to ionospheric delay (IOD) and the prior tolerance. The ionospheric residual method can detect a single cycle slip; thus, when the ionosphere changes slowly, this method is more suitable than other existing approaches [3]. More recently, a method was proposed that integrates the forward and backward moving window averaging (FBMWA) algorithm with the second-order, time-difference phase ionospheric residual (STPIR) algorithm [4]. The combined FBMWA and STPIR method can detect cycle slips very well; however, when some cycle slips occur in the near field, pseudorange noise will still affect the cycle slip repair accuracy. More recently, a single-frequency cycle slip detection and repair technique based on the generalized likelihood ratio (GLR) test was presented [6]; however, this method may also be influenced by continuous cycle slips.

The second kind of cycle slip detection involves observation equations. Bastos and Landau [7], Gao and Li [8], Colombo et al. [9], and Lee et al. [10] investigated cycle slip detection algorithms based on observation equations. However, cycle slip detection methods based on observation equations have limitations. Among the studies mentioned above, Colombo et al. [9] combined GPS and INS (internal navigation system) parameters to realize cycle slip detection, and Lee et al. [10] further studied a GPS/INS cycle slip detection approach based on the cumulative sum (CUSUM) method. However, the GPS/INS cycle slip detection method must use an INS system and thus is not suitable for cycle slip detection with LEO satellites. A method combining geometry-free (GF) and ionospheric-free (IF) double-difference observations can realize cycle slip detection better than the above mentioned techniques [11], but is suitable only for static CORS stations. To realize cycle slip detection for dynamic GPS receivers, a method combining the LAMBDA carrier phase ambiguity resolution approach with carrier phase observations was used for dynamic PPP (precise point positioning) cycle slip detection and repair [12]. More recently, a new dynamic PPP cycle slip detection method was proposed by combining the ambiguity characteristics of the single-difference observation equation with the least square method [13]. However, these dynamic PPP cycle slip detection methods are aimed at mainly ground observation stations. For cycle slip in onboard GPS receivers, a cycle slip detection method based on STP (second-order time-difference of the LEO satellite’s position) and STG (second-order time-difference geometry-free) was proposed [14]. This method does not need pseudorange observations and can reach a high detection success rate when the elevation angle of the satellite is low (less than 2.1°). However, for this method to be successful, the prior acceleration must be known, so this technique can be used as a supplement.

This research makes several contributions. First, a special time-difference and satellite-difference observation value is established. If cycle slips occur in the reference satellite, the reference satellite cycle slip error will be enhanced for all of the special observations. Second, to process the cycle slips of GPS receivers onboard LEO satellites, the enhanced error method (EEM) is proposed, which is not affected by pseudorange observations. The EEM is applicable for detecting and repairing continuous cycle slips. Third, the ability of the EEM at detecting cycle slip is verified by evaluating simulated cycle slips. We will introduce in detail the EEM cycle slip detection algorithm in the next section. Then, the method proposed will be tested in the following section and summarized in the last section.

2. Observational Model and Methodology

The EEM establishes a special observation value based on single-frequency (SF) and narrow-lane carrier phase observations. Narrow-lane carrier phase observations are built by dual-frequency carrier phase observations. This special observation differs across time and satellite (ATS) data. The EEM carries out independent time differences for single-frequency and narrow-lane carrier phase observations and establishes single-frequency time-difference (SFTD) and narrow-lane time-difference (NLTD) observation values. We select a satellite as the reference satellite and other satellites as nonreference satellites. Then, ATS observations are obtained by satellite differences, namely, the differences between the SFTD observations of the reference satellite and the NLTD observations of each nonreference satellite. If a cycle slip or an outlier occurs in the reference satellite, the ATS observations corresponding to the reference satellite will be corrupted, and thus, the correction values of the ATS observations will be significantly enlarged. To detect cycle slips in all satellites, the EEM selects each satellite as the reference satellite and builds the corresponding equation group. The EEM solves these observation equation groups according to the weighted least-squares adjustment (LSA) criterion and obtains the correction values of the ATS observations. The EEM uses these correction values to construct corresponding to different equation groups and carries out a chi-square distribution test for these values of . If cycle slips occur in several satellites, the values that do not obey the chi-square distribution exceed one. Therefore, the satellite corresponding to the maximum will be marked provisionally as a cycle slip satellite. Then, the EEM iteratively processes the other satellites. Based on the LSA criterion, cycle slips and outliers can be estimated by the float solution of changes in the ambiguities of cycle slip satellites.

In this section, NLTD and ATS observations will be introduced first, and then, the principle of the EEM, namely, the enhancement of cycle slip errors and the cycle slip detection strategy, will be explained.

2.1. Time-Difference Observation Equation

Narrow-lane carrier phase observations can be expressed as , where and are GPS single-frequency carrier phase observations and and are the frequencies corresponding to carrier phase observations and , respectively. The time-difference equation of can be given as [1, 15] where represents the NLTD observation of satellite and is constructed by the difference in between epochs and ; is the narrow-lane observation of satellite in epoch ; is the narrow-lane observation in epoch ; is the geometric distance from the receiver to satellite in epoch ; is the geometric distance from the receiver to satellite in epoch ; is the receiver clock error in epoch ; is the receiver clock error in epoch ; is the satellite clock error in epoch ; is the satellite clock error in epoch ; is the speed of light; is the noise difference between the ambiguities of and ; is the variation in the IOD between epochs and ; and represents the noise in the observation data. The time-difference equation can eliminate the ambiguity parameter. If there is no cycle slip or outlier, the value of is 0. is calculated as follows: where , , and are the differences in the corresponding coordinates of the LEO satellite between adjacent epochs; , , and are the LEO satellite coordinates in epoch ; , , and are the GPS satellite coordinates in epoch ; and is similar to .

2.2. ATS Observation Equation

Using and as examples, according to formula (1), the SFTD and NLTD observations are given as and , respectively. The superscript represents the reference satellite. The ATS observation equation is given as where represents the ATS observation; is the difference between and () in nonreference satellite ; is the difference between and () in reference satellite ; is the difference between and (); is the difference between and (); is the difference between and (); and is the difference between and ().

Supposing that there are satellites in an epoch, the satellite corresponding to is called the reference satellite, and the satellite corresponding to is called the nonreference satellite. ATS observations can be constructed by combining the reference and nonreference satellites. The observation equation group for ATS observations is established as follows [1]: where is the coefficient matrix; is the constant matrix; is the correction vector of ATS observations to be solved; and is the parameter vector to be solved. The expression is

The expression is

Supposing is the weight matrix of the ATS observations, then, according to the LSA criterion [16], can be expressed as

The initial values of , , and and of , , and can be obtained by single-point positioning (SPP) using ionosphere-free pseudorange observations. , , and are obtained by the precise ephemeris provided by the Center for Orbit Determination in Europe (CODE) through Lagrange interpolation. The satellite clock delays and are obtained by the precise clock error provided by the CODE through Lagrange interpolation [17, 18]. If there is no cycle slip, the theoretical value of is 0. , which is the variation in the IOD, is very smooth in LEO satellites when the elevation angle is larger than zero [19, 20]; however, when the sample interval is 1 s or higher, the variation in the IOD is distributed mainly in the range of cycles at low elevation angles. Hence, the variation in the IOD can be ignored (details are presented in Section 3.3).

2.3. Theory of Enhanced Cycle Slip Error

For convenience of discussion, it is assumed that the components of the observation correction vector are independent from each other. If no cycle slip occurs, conforms to the multivariate normal distribution [1, 21], and the expectation of is . If only the observations of the reference satellites have cycle slips or outliers, will obey a multivariate skewed distribution, and the expectation of is

If cycle slips or outliers occur in nonreference satellite , the expectation of is

According to formula (8), the EEM expands the cycle slips and outlier errors of the reference satellites to affect all the observations so that no ATS observations obey a normal distribution. In this case, the occurrence of cycle slips in the epoch can be effectively detected by using to construct and conducting a chi-square distribution test. According to formula (9), if the selected reference satellite is not a cycle slip satellite, the corresponding may not obey a chi-square distribution, but its value is less than the value corresponding to the cycle slip satellite as the reference satellite. The statistics are constructed as follows:

The original hypothesis and the alternative hypothesis are where is the prior standard deviation (STD) of ; is the significance level; and and are the original hypothesis and the alternative hypothesis, respectively. When meets the condition of , no cycle slip is found in this epoch. Otherwise, the alternative hypothesis is accepted; i.e., cycle slip occurs in this epoch.

Formula (10) shows that the selection of the prior value is the key to establishing . We find that is affected by the sampling interval [2224]. The value of is closely related to the data sampling interval and fluctuates greatly, so it should not be set as a fixed prior value. Therefore, an automatic method to select is given here: supposing the current epoch is , the mean prior STD from the initial epoch to epoch is , and the variance of is . The formula for choosing reasonable values of and for epoch is as follows:

From and of epoch , combined with the PauTa criterion, the of epoch can be detected:

If (13) is true, epoch can use as the prior STD; otherwise, epoch still uses as the prior STD.

2.4. Novel Method for Cycle Slip Detection and Repair

To determine cycle slips for satellites in epoch , the specific flowchart of the novel cycle slip detection method is given in Figure 1. The corresponding procedures are as follows: (1)Supposing that there are satellites in adjacent epochs, equation (4) is used to construct observation equation groups. Each group has observation equations and corresponds to a reference satellite (2)According to equation (7), observation equation groups are solved to obtain observation correction vectors. According to equation (10), these vectors are used to construct chi-square values . The subscript represents the chi-square value corresponding to reference satellite (a)If only a single value does not obey a chi-square distribution, provisionally mark as a cycle slip satellite. Move to step (4)(b)If more than one value fails to obey a chi-square distribution, provisionally mark the satellite that corresponds to the maximum value as the cycle slip satellite and eliminate the satellite. Set . Move to step (1)(c)If , move to step (3)(3)If all cycle slip satellites have not been exactly marked, the normal satellites will be eliminated until 4 satellites remain. Considering this extreme condition, a more robust detection strategy is given (a)Initialize ; then, satellite can be temporarily eliminated. According to equations (4) and (10), other satellites can be used to establish new values of . Renew , and then, obtain the corresponding (b)Select the satellite corresponding to the minimum as a cycle slip satellite. Set. Move to step (a) until all cycle slips have been marked(4)These cycle slips are marked as unknown parameters. These parameters can be calculated according to the float solution of changes in the ambiguities of cycle slip satellites [25, 26]. Supposing that cycle slip occurs only in satellite , the new equation groups can be expressed aswhere is the coefficient matrix and is an unknown parameter of . can be directly used to calculate or . All of the unknown parameters can be given as

can be solved as

The cofactor matrix of unknown parameters can be given as

is used to divide cycle slips and outliers.

While the float solution of a cycle slip is estimated, fixing the cycle slip as an integer is important. The LAMBDA method is generally adopted to round cycle slips to the nearest integer [12, 27]. However, as shown in (17), the accuracy of is not affected by pseudorange observations. Therefore, the cycle slip can be accurately rounded to the nearest integer. We use the STD values of the float solutions to distinguish outliers and cycle slips. The detection formula is written as where is the cofactor of in . On the one hand, while (21) is true, round () can be viewed as a cycle slip. On the other hand, if (21) is false, the satellite corresponding to will be marked as an outlier satellite. (5)There is a correlation between the ATS observation values, and the corresponding prior weight matrix is no longer a diagonal matrix. Therefore, cannot be simply defined as a unit matrix [28]. The strategy for processing is given below: (a) is a unit matrix in the initial epoch (b)The prior cofactor matrix can be given as (c)According to the formula , calculate the cofactor matrix of (d)According to the formula , calculate the posterior weight matrix and use it as the prior weight matrix of the next epoch(e)Use equations (12), (13), and (14) to renew the prior STD

It should be mentioned that the validity of the EEM depends on the quality of observations in the present epoch. If the number of satellites without cycle slips is less than 5, the EEM is not valid. Therefore, to detect and repair cycle slips more robustly, the EEM can be combined with other methods to address cycle slips.

3. Simulated Case Study

The case data are taken from the onboard GPS receiver of the Swarm-A satellite provided by the European Space Agency (ESA). The precise GPS ephemeris with a 15 min sampling interval and the precise satellite clock error file with a 30 s sampling interval are provided by the CODE. The date on which the Swarm-A satellite data were acquired is DOY 145 of 2018 (GPS time 2018:05:25), and the sampling interval is 1 s. Table 1 lists the detailed web addresses of the data sources.

The certain times at which the observation data of a satellite are acquired are regarded as observation satellite epochs. A complete observation satellite epoch includes carrier phase observations ( and ) and pseudorange observations ( and ). Incomplete observation satellite epochs are eliminated. On DOY 145 of 2018, Swarm-A had 684087 observation satellite epochs. Complete observation satellite epochs account for 97.6% of the total. This means that the observation results of Swarm-A on DOY 145 are relatively complete and can be used as test data.

3.1. Special Cycle Slip Pair Detection and Repair

To test the validity of the new method, from epoch 60 (GPS time 2018:05:25:00:01:00) to epoch 670 (GPS time 2018:05:25:00:11:10), some special simulated cycle slip pairs are inserted into satellite G01. Table 2 shows the distribution of the simulated cycle slips and the detection and repair results. The detection results indicate that these special cycle slip pairs can be detected by the EEM and by the FBMWA and TurboEdit methods. On the one hand, for these cycle slip pairs inserted every 30 epochs, the noise in the pseudorange observations is effectively smoothed by the FBMWA algorithm. However, several repair errors occur in epochs 410, 520, and 670. The size of the moving window is 25 seconds [4]. It may be that the short size of the moving window is not enough to reduce the noise level in the pseudorange observations. On the other hand, the EEM is sensitive to special cycle slip pairs such as (1, 1), (3, 4), and (9, 7). According to equation (17), these simulated cycle slip pairs can be repaired well. These estimated cycle slip values are quite close to the simulated cycle slips. This is because equation (17) does not use pseudorange observations to repair cycle slips. As a result, these cycle slips can be accurately rounded to the nearest integer. Therefore, combining the EEM with the FBMWA and TurboEdit algorithms is more robust at processing cycle slips than using any of these methods alone.

To verify the ability of the EEM to detect and repair continuous cycle slips from epoch 300 (GPS time 2018:05:25:00:05:00) to epoch 304, continuous cycle slips are inserted into G01. Table 3 lists the detection and repair results. The results show that the EEM algorithm can detect and repair cycle slips that occur in continuous epochs.

3.2. Variation in the Ionospheric Delay

The ionospheric delay (IOD) commonly varies quickly in LEO satellites, as shown in previous studies [4, 14]. However, rapid variation in the IOD is not always observed when the elevation angle is larger than zero [19, 20]. Generally, onboard GPS receivers can provide a high sampling interval, such as 1 s [2931]. Therefore, a significant question is raised as to whether the variation in the IOD can be neglected in the processing of cycle slips for LEO satellites with a high sampling interval (1 s) and a low elevation angle. To solve this question, we use a geometry-free linear combination to assess the variation in the IOD [20] and use equation (17) to estimate and . Figures 2(a)2(d) show the variations in the IOD (sampling interval: 1 s). Figures 3(a)3(h) show the variations in and . As these figures show, we can reach several conclusions. First, the ambiguity and IOD vary slowly when the elevation angle is larger than 10°. These variations are distributed mainly in the range of cycles. Second, the maximum variation is lower than 0.1 cycles when the elevation angle approaches zero. This finding confirms that the IOD variation is stable. Finally, the variations in and are significantly small because the estimates of and are not affected by pseudorange noise, IOD variations, or other errors, even at low elevation angles. These phenomena suggest that (1) the variations in the ambiguity and IOD obey a normal distribution and that (2) the variation in the IOD does not obviously corrupt time-difference carrier phase observations. Hence, the variation in the IOD of time-difference carrier phase observations in equations (3) and (17) can be ignored for the proposed cycle slip detection and repair method.

To investigate the performance of the proposed method at low elevation angles, we select G01 and G17 as examples (sampling interval: 1 s). From epoch 360 (GPS time 2018:05:25:00:01:00) to epoch 930 (GPS time 2018:05:25:00:15:30), 20 pairs of the simulated cycle slip pair (9, 7) are added to G01. The elevation angles in the period from epoch 600 to epoch 930 are all lower than 15°. From epoch 360 to epoch 570 (GPS time 2018:05:25:00:9:30), 8 pairs of the simulated cycle slip pair (1, 1) are added to G17. The elevation angles in the period from epoch 360 to epoch 570 are all lower than 15°. These cycle slip pairs are described in Figures 4(a) and 4(b).

The cycle slip repair results for G01 and G17 are shown in Figures 5(a)–5(f). In each plot, the values of and the variation in the elevation angle are given. To visualize the results better, values of larger than 1.2 cycles were set as 1.2 cycles. These cycle slip epochs are indicated by vertical lines. All simulated cycle slip pairs (9, 7) and (1, 1) in G01 and G17 were detected and repaired. As shown in Figures 5(e) and 5(f), in the simulated cycle slip epochs, the values of and are close to the simulated cycle slip pair (1, 1). This result verifies that these can be easily rounded to their nearest integers.

3.3. Standard Deviation of the Observation Correction Vector

Constructing requires the prior STD of the observation correction vector (). To study the distribution of the STD, the STD values corresponding to each epoch and reference satellite are calculated for the Swarm-A observation data on DOY 145 according to sampling intervals of 1 s, 5 s, 10 s, 15 s, 30 s, and 60 s. According to the formula , the mean posterior STD values corresponding to different reference satellites in each epoch are calculated [1]. The statistical results related to the posterior STD are listed in Tables 4 and 5, which show that the posterior STD of the observation correction vector under different reference satellites is closely related to the data sampling interval, and the value of the posterior STD increases with a decrease in the sampling interval. At the same sampling interval, the mean difference in the posterior STD of the ATS observations constructed by different reference satellites is small. The mean difference in the posterior STD between the sampling intervals of 1 s and 60 s reaches an order of magnitude. The accuracy of is better than 5 cm when the sampling interval of the observation data is higher than or equal to 15 s, whereas the accuracy of is less than 5 cm when the sampling interval of the observation data is less than or equal to 30 s. This may be because with a decrease in the sampling interval, the variation in the IOD between epochs is expanded and cannot be greatly weakened by the difference between epochs at a low sampling interval [4]. The above data analysis shows that the sensitivity of the EEM algorithm increases with the increase in the sampling interval of the observation data, so the algorithm is more suitable for data acquired with a high sampling frequency. At the same time, this analysis confirms that formula (12) should be used to estimate the prior STD.

3.4. Performance of the EEM to Enhance Cycle Slip Errors

To verify the ability of the EEM to amplify cycle slip errors, a simulated outlier of 0.5 cycles and a cycle slip of 1 cycle are inserted into the carrier phase of G11 in the observation epoch corresponding to GPS time 2018:05:25:00:06:30. Then, we calculate the values corresponding to the different inserted simulated cycle slips and outliers. Table 6 lists the values. The results show that when G11 is selected as the reference satellite, its value is much larger than the values of the other observation satellites selected as the reference satellite. To obtain experimental results that are more statistically representative, an outlier of 0.5 cycles is inserted into satellite G11 every 10 s during the continuous observation epoch from GPS time 2018:05:25:00:00:20 to 2018:05:25:00:16:00 with a sampling interval of 1 s. Figure 6 shows the values for the satellites with inserted and no inserted simulated values. The results show that when satellite G11 is selected as the reference satellite with an inserted outlier of 0.5 cycles, the corresponding value is greater than the values of the other reference satellites without inserted simulated outliers. These results verify that the EEM algorithm can effectively enlarge the errors caused by small cycle slips and outliers.

4. Conclusion

A method for detecting and repairing onboard GPS receiver cycle slips is presented. This method combines single-frequency and narrow-lane observations to construct ATS observations, which can significantly enlarge the small errors caused by cycle slips and outliers. By amplifying the influences caused by cycle slip and outliers, the EEM algorithm can detect cycle slips more stably than other existing algorithms.

The sensitivity of the EEM algorithm is evaluated by inserting simulated cycle slips into different reference satellites. The results show that the EEM, FBMWA, and TurboEdit methods have similar repair accuracies for cycle slips. However, compared with the FBMWA and TurboEdit methods, the EEM is not affected by pseudorange observations, which makes up for the shortcomings of the FBMWA and TurboEdit algorithms. The proposed algorithm is applicable even to observation satellite epochs with continuous cycle slips and outliers.

The proposed EEM cycle slip detection method can effectively detect and repair cycle slips and outliers and is still applicable to observation epochs with continuous cycle slips. However, for low-frequency observation data (sampling interval less than 30 s), the EEM algorithm cannot eliminate the impact of variation in the IOD. In addition, the EEM algorithm depends on the quality of the observation data from other satellites at the same time. Therefore, to thoroughly detect onboard GPS receiver cycle slips, the EEM, FBMWA, and TurboEdit algorithms should be combined to process cycle slips simultaneously.

Data Availability

The space-borne GPS observation data used in our work can be freely accessed at https://swarm-diss.eo.esa.int/. The GPS satellite ephemeris and clock products are obtained from ftp://cddis.nasa.gov/gnss/products.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

We thank ESA for providing space-borne GPS data of Swarm-A and CODE for providing precise GPS ephemeris and satellite clock error file. This research is supported by the National Natural Science Foundation of China (grant numbers 41774001, 41774021, and 41874091) and the SDUST Research Fund (grant number 2014TDJH101).