Detecting and Repairing Inter-system Bias Jumps with Satellite Clock Preprocessing

: The key to performing successful multi-GNSS (Global Navigation Satellite System) precise point positioning is calibrating ISB (inter-system bias) in di ﬀ erent systems. We can use the method of modeling to eliminate the ISB error. However, the ISB series are commonly discontinuous, as they contain jumps and outliers caused by day boundaries, gaps, or outliers in the precise clock products, which will break the ability of modeling. Thus, before modeling the ISB, we must remove outliers and repair jumps to improve the ISB continuity and achieve a continuous and smooth ISB series. Preprocessing on precise clock products is focused on in this study for the detection of ISB jumps and their repair. From the results, a positive correlation is revealed within the residuals of satellite clock o ﬀ set and ISB di ﬀ erences between adjacent days. This ﬁnding implies ISB continuity can be improved through the preprocessing of precise clock products. It is also found that the exact reason for the occurrence of ISB jumps is the presence of extrema (i.e., maximum or minimum points) in the frequency domain. From the clock data in the frequency domain, larger extrema are identiﬁed directly. Meanwhile, for the detection of smaller extrema, a robust estimation method based on the median ﬁlter was applied. Then, all smaller extrema were classiﬁed into three types. Di ﬀ erent preprocessing methods were applied for every type. After that, a new preprocessed precise clock product was obtained. With this preprocessed satellite clock product, the ISB continuity was substantially improved, and the improvement in the ISB continuity can reach 85.1%, on the average. These results indicate that for detecting and repairing ISB jumps, the proposed preprocessing method on satellite clock products is very e ﬀ ective.


Introduction
GNSS (Global Navigation Satellite System) is developing quickly along with the modernizing of American GPS (global positioning system) and Russian GLONASS (global navigation satellite system) and the launching of Europe's Galileo, China's BDS (BeiDou navigation satellite system) with new frequency signals being transmitted. Thus, the increasing number of satellite and frequency data can be obtained and used. Together with the applications of the precise products from the IGS MGEX (International GNSS Service multi-GNSS experiment), numerous studies on the multi-system and multi-frequency combined positioning have been made [1]. A general multi-frequency modeling strategy for precise point positioning (PPP) ambiguity resolution was presented by Gu et al. [2]. Odolinski et al. [3] made a discussion on the RTK (real-time kinematic) research with four CDMA (code division multiple access) satellite systems. From the results, they found the combination of multi-system enhanced the RTK performance. Compared with the single GPS PPP, Chen et al. [4]

Estimation Model and Processing Strategies
To attain the ISB series between the GPS and BDS, we carried out integrated PPP solutions using the data from the nine stations mentioned above. The ISB parameter was estimated together with other unknowns.

BDS/GPS PPP Model for ISB Estimation
Usually, to eliminate the first-order error of ionospheric delay, IF (ionosphere-free) linear combined phase and pseudorange observations are implemented in PPP. The observation functions in a single system within a receiver and a satellite are illustrated by [17]: IF combined phase and pseudorange observations are L and P. The range between the satellite and receiver is expressed by ρ. The light speed is shown by c. Receiver clock error is expressed by dt r and the satellite clock error is expressed by dt s . The mapping function is M and zenith tropospheric delay is ZTD. The subscripts r and s denote the parameters related to the receivers and satellites, respectively. B r , B s expresse phase hardware delay biases, b r and b s are pseudorange hardware delay biases. Ambiguity is denoted by N and ς, ε is the residual unmodeled errors, including the measurement noise.
The pseudorange hardware delay biases b r and b s can be absorbed by the clock errors c · (dt r − dt s ). In data processing, the carrier phase hardware delay biases B s and B r are not considered because they can be assimilated into ambiguity, as they are the satellite-dependency and stability over time [18,19]. Equation (1) can be expressed after implementing precise satellite orbit and clock products as: In this case, d t r is the newly formed receiver clock error and N is the re-formed ambiguity.
While utilizing the precise products in Equation (2), satellite hardware delays can be neglected and are assimilated by precise satellite clock [20]. Equation (3) shows that ambiguity is not an integer feature anymore because of including bias term B r − b r . It is defined as the uncalibrated phase delay error [21].
Comparing to the single system observation equations, the estimation of an extra ISB parameter is required in the BDS/GPS combined PPP model. The corresponding combined observation functions are illustrated as: Here, d t r is the common receiver clock error in GPS and BDS. GPS/BDS inter-system bias parameter is expressed by ISB, which is also treated as an unknown parameter, and simultaneously estimated with other parameters, like receiver clock errors as well as coordinates. The superscript C means BDS, and G is GPS system. Thus, the vector of estimated parameters is: In Equation (5), receiver position increments are separately expressed in three components x, y, z. The German Research Centre for Geosciences (GFZ) MGEX precise orbit and clock products, named GBM products, were applied in the processing. For the estimation of unknown parameters, the method of extended Kalman filter was employed. Along with other unknowns, GPS/BDS ISB was estimated as a piece-wise constant parameter every 5 min. For data processing, the code observation types C1C and C2W, and corresponding phase observation types L1C and L2W were applied in the GPS system, which are abbreviated as L1 and L2. For the BDS system, the code observation types C1I and C7I, together with phase types L1I and L7I were utilized with the abbreviation of B1 and B2. IF linear combination was implemented so that the first-order effect of ionospheric delay could be removed. The cutoff elevation angle was set to 7 • . Estimation of zenith wet delays was treated as a random walk process. The constraint on processing noise was empirically set to 1 cm 2 /h. Elevation-dependent weighting for the observations was implemented while the elevations were below 30 • , otherwise, the weight was set to 1. In BDS and GPS, the initial standard deviation of carrier phase observations was set to 0.003 m, as well as the pseudorange observations were set to 0.3 m. In BDS/GPS combined PPP, following the empirical model, the weights between GPS, MEO in BDS, IGSO in BDS, and GEO in BDS were set to 10:10:5:1, as orbit accuracy in GEO (geostationary orbit) satellites is at a lower level than the IGSO (inclined geostationary earth orbit) and the MEO (medium earth orbit) satellites [22]. In the case of each continuous satellite arc, phase ambiguity was estimated as a float constant. BDS/GPS observation models and data processing strategies are summarized in Table 2.

The Incentive of Precise Satellite Clock Preprocessing
The estimated BDS/GPS ISB series of all nine stations are expressed in Figure 2. Three types of series can be noticed here. The first one is equipped with Trimble NetR9 receivers, which include three different stations such as JFNG, SIN1, and CUT0. The second one consists of three stations HKSL, HKWS, and HKOH, which are installed with Leica GRX1200+GNSS receivers. The last one includes the rest of the stations NNOR, MAJU, and NAUR, which are equipped with SEPT PolaRx4/4TR receivers. Figure 2 illustrates that the ISB series for the same category receiver are similar to each other. That is because between two satellite systems, the hardware delay is the major component of ISB, and for the same category receiver, the hardware delay is theoretically similar. ISBs in epochs 4800 to 6400 are included in the red box, where two jumps are observed. In the enlarged view, they are more evident.
Both the gaps appear in the junction of the day-pair 17/18 and 19/20 August. The irregular variations in the satellite clock are causing these gaps. There are mainly three kinds of these irregular variations: one is satellite clock PM (phase modulation), which is considered as the deviation between the initial and true phase, the other one is FM (frequency modulation), and the last one is the switch of the reference clock and the outliers. Usually, the satellite clock error is formed as a linear function with a first-order item and a constant term. The value of the constant term is considered as the phase term of satellite clock error, and the coefficient of first-order term, i.e., slope, is regarded as the frequency term. Therefore, PM means variation of constant term, and FM is the change of the coefficient of first-order term. Both of them can lead the irregular variation of satellite clock. The satellite clock is substituted as a known term into the processing equation for PPP processing. Estimated ISB varies correspondingly when satellite clock experiences jump due to linear dependence between them.
considered as the deviation between the initial and true phase, the other one is FM (frequency modulation), and the last one is the switch of the reference clock and the outliers. Usually, the satellite clock error is formed as a linear function with a first-order item and a constant term. The value of the constant term is considered as the phase term of satellite clock error, and the coefficient of first-order term, i.e., slope, is regarded as the frequency term. Therefore, PM means variation of constant term, and FM is the change of the coefficient of first-order term. Both of them can lead the irregular variation of satellite clock. The satellite clock is substituted as a known term into the processing equation for PPP processing. Estimated ISB varies correspondingly when satellite clock experiences jump due to linear dependence between them.
Analysis of correlation coefficients between clock offsets and ISB was carried out to verify the linear dependence between them. From the MGEX GFZ analysis center (GBM), BDS precise clock products during August 2015 are implemented here. The values of BDS satellite clock offsets from GBM products for C01-C14 are represented by Figure 3, except C13. We can see from Figure 3 that large data jumps and DIs (data interruptions) continually appear in the clock offsets of BDS satellites. The large clock jumps are defined as PMs. They happen more frequently in GEO satellites (C01-C05). Probable cause is that the GEO satellites perform more orbital maneuvers. As the magnitudes of PM and DI are huge in the satellite clock, they can be easily detected and marked in Figure 3. Table 3 summarizes these corresponding detected times of PM and DI.  Analysis of correlation coefficients between clock offsets and ISB was carried out to verify the linear dependence between them. From the MGEX GFZ analysis center (GBM), BDS precise clock products during August 2015 are implemented here. The values of BDS satellite clock offsets from GBM products for C01-C14 are represented by Figure 3, except C13. We can see from Figure 3 that large data jumps and DIs (data interruptions) continually appear in the clock offsets of BDS satellites. The large clock jumps are defined as PMs. They happen more frequently in GEO satellites (C01-C05). Probable cause is that the GEO satellites perform more orbital maneuvers. As the magnitudes of PM and DI are huge in the satellite clock, they can be easily detected and marked in Figure 3. Table 3 summarizes these corresponding detected times of PM and DI. Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 19 Analyzing Figure 2, Figure 3 and Table 3, it can be noticed that, in Figure 2, the time when ISB jumps occur is not simultaneous with the time of PM happening from Figure 3 and Table 3. Based on the findings, it can be said that the true reason for ISB discontinuity is not the clock PM.  Analyzing Figure 2, Figure 3 and Table 3, it can be noticed that, in Figure 2, the time when ISB jumps occur is not simultaneous with the time of PM happening from Figure 3 and Table 3. Based on the findings, it can be said that the true reason for ISB discontinuity is not the clock PM.
For further revelation of the linear dependency between satellite clock offsets and ISB, numerical verification was carried out. The test period for numerical verification was from 17 August to 21 August, which means that the experiments were carried out on the four day-pairs: August 17/18, 18/19, 19/20, and 20/21, 2015. Residual of the satellite clock offsets ∆t in a day-pair was calculated by making the difference between Pt (predicted satellite clock from the polynomial fitting on basis of the satellite clocks from the last ten epochs on 1st day) and Rt, which is the real satellite clock in the 1st epoch of 2nd day. The whole process is summarized as: Remote Sens. 2020, 12, 850 9 of 18 where F(Rt 1 , Rt 2 , . . . , Rt 10 ) denotes polynomial fitting using the satellite clock offsets from the last 10 epochs in the 1st day. By the extrapolation of F(Rt 1 , Rt 2 , . . . , Rt 10 ), Pt is derived. f e is the 1st epoch on the 2nd day. ISB difference between adjacent days is derived by: The daily mean ISB value is indicated as ISBmean. ISB value is expressed by ISB with epoch number n ranging from 1 to 288. The subscript of a day-pair is expressed by i and i − 1. ISB difference between the 1st and 2nd day is denoted as ∆ISB. Figure 4 shows the ISB differences along with the satellite clock offset residuals in the four day-pairs. It is shown that ISB differences in these nine stations express different magnitudes in different day-pairs. In day-pairs 19/20 and 20/21, the most noticeable changes appear. It indicates that a large jump of satellite clock offsets occurred on 20 August, 2015. We checked the description comments on the reference clock in the header of the GBM product file, and concluded that the large jump is ascribed to the reference clock switch in the GBM precise clock product. The level of ISB differences in the nine stations is close in every day-pair. In addition, like the ISB difference, satellite clock offsets express the same phenomenon, which is for all satellites, the level of clock offsets residual is also close in each day-pair.
The daily mean ISB value is indicated as ISBmean . ISB value is expressed by ISB with epoch number n ranging from 1 to 288 . The subscript of a day-pair is expressed by i and 1 i  . ISB difference between the 1st and 2nd day is denoted as ISB  . Figure 4 shows the ISB differences along with the satellite clock offset residuals in the four daypairs. It is shown that ISB differences in these nine stations express different magnitudes in different day-pairs. In day-pairs 19/20 and 20/21, the most noticeable changes appear. It indicates that a large jump of satellite clock offsets occurred on 20 August, 2015. We checked the description comments on the reference clock in the header of the GBM product file, and concluded that the large jump is ascribed to the reference clock switch in the GBM precise clock product. The level of ISB differences in the nine stations is close in every day-pair. In addition, like the ISB difference, satellite clock offsets express the same phenomenon, which is for all satellites, the level of clock offsets residual is also close in each day-pair.  Table 4. At the same time, Table 5 includes satellite clock offset residuals and their statistics. Table 5 also includes some blank entries, which are due to the DI of the satellite clock in the adjacent days. Comparing Table 4 and Table 5, it is noticed that among all the four day-pairs, ISB differences at every station are correlated with the residuals of the satellite clock in all satellites. We can also see that if the satellite clock offset jumps take place, at the same time ISB jumps will appear. In the next part, we will make an analysis of the corresponding correlation coefficients.
Correlation coefficients were implemented as indicators of the correlation degree between two parameters. The correlation coefficients between the residuals of satellite clock offset in every satellite and ISB differences in every station were calculated from Table 4 and Table 5, and the results of correlation coefficients are shown in Table 6. Correlation analysis was only carried out using the satellites with residual data. It can be found from Table 6 that there is a strong positive correlation between every pair of the residual of satellite clock offset and ISB difference, since all the correlation coefficients are close to 1. In accordance with these findings, the incentive of satellite clock  Table 4. At the same time, Table 5 includes satellite clock offset residuals and their statistics. Table 5 also includes some blank entries, which are due to the DI of the satellite clock in the adjacent days. Comparing Tables 4 and 5, it is noticed that among all the four day-pairs, ISB differences at every station are correlated with the residuals of the satellite clock in all satellites. We can also see that if the satellite clock offset jumps take place, at the same time ISB jumps will appear. In the next part, we will make an analysis of the corresponding correlation coefficients. Correlation coefficients were implemented as indicators of the correlation degree between two parameters. The correlation coefficients between the residuals of satellite clock offset in every satellite and ISB differences in every station were calculated from Tables 4 and 5, and the results of correlation coefficients are shown in Table 6. Correlation analysis was only carried out using the satellites with residual data. It can be found from Table 6 that there is a strong positive correlation between every pair of the residual of satellite clock offset and ISB difference, since all the correlation coefficients are close to 1. In accordance with these findings, the incentive of satellite clock preprocessing becomes patently clear, which is to improve ISB continuity. In the next sections, we will present a detailed preprocessing process on the satellite clock to enhance ISB continuity. Table 6. The given value is 1-ρ, and ρ means the correlation coefficient between the residuals of satellite clock offset and ISB differences. So the smaller the given value, the closer the coefficient is to 1; unit: 10 −5 .

Results and Discussion
In this section, a preprocessing method on the precise satellite products is proposed. First, the satellite clock product is converted from time domain into frequency domain, since it is more sensitive to the detection of ISB jumps. Then, extrema in the frequency domain clock data are detected and classified. After that, clock preprocessing is made to get a new preprocessed satellite clock product.
At last, a test on the ISB continuity improvement is carried out after getting re-estimated ISB using processed satellite clock products.

Conversion from Time Domain into Frequency Domain
Larger clock jumps caused by considerable PMs can be easily detected in the time domain of the satellite clock in Figure 3. But as said before, the ISB continuity will not be destroyed by this type of clock jump. Thus, the true reason causing the ISB discontinuity needs to be found out. In addition, some smaller clock jumps can take place due to reference clock substitutions, weak PMs, clock day boundaries, which can be drowned in the large-scale satellite clock offset series. In order to overcome the above problems, the satellite clock is converted into the frequency domain from the time domain, since clock jumps, specifically small clock jumps, can be detected more sensitively by the satellite clock in the frequency domain.
Firstly, clock offset time series, excluding DIs as well as numerous PM epochs, are converted into CFD (clock frequency data). Fractional frequency data can be derived by taking the first differences x i+1 − x i of the phase data and dividing by the sampling interval τ, so it can be expressed as [23]: Clock offset in the time domain in two adjacent epochs is expressed by x i and x i+1 . Interval τ is set to five minutes. In the frequency domain, clock data is denoted as y i .
CFD series in August 2015 are shown in Figure 5. From Figure 5, we can see that the CFD series of every satellite has two extrema in epochs 5760 and 5472, the corresponding epoch date is the start time on 20 and 21 August, 2015. They are synchronous with the occurrence times of the most remarkable ISB gaps. Therefore, the true reason for those ISB jumps is the extrema in the CFD. It can be also noted that ISB continuity can be destroyed by other smaller extrema in the CFD as well. In another scope, it can be said that the most significant residuals of clock offset are indicated by the CFD extrema, and ISB differences between two adjacent epochs are considered as those ISB jumps. They have a great correlation. This finding deeply provides evidence of a strong correlation of the residuals of satellite clock offset with ISB differences.

Detection of Small Extrema in CFD
The median value is better able to reveal the characteristics of most data than the mean value. The median value more effectively avoids the impact of outliers. CFD extrema detection is more improved while applying the method of median filter-based robust estimation [24]. Small ISB jumps, which are not possible to observe directly from the CFD series, are detected by implementing this method. Following the median principle, the root mean square error of unit weight, marked with σ 0 , is expressed as the normalized median for the absolute deviation between CFD and the median of CFD [25]. It is illustrated as: CFD at epoch i are indicated by y i . Median of y i is denoted by m. The deviation between y i as well as m is expressed as v i . The root mean square error of unit weight is expressed by σ 0 . n · σ 0 is considered as the threshold value. y i can be considered as an extremum of CFD if equation |v i | ≥ n · σ 0 is satisfied. If not, y i is considered as a regular CFD point. An empirical model specifies the variability of n. In this study, it is selected as 10, since we use the method of exhaustion with the value from 1 to 10 to check the result of extrema detection, and when the value is selected as 10, the performance of detection is best.
After the detections of CFD extrema, they will be classified by different causes. Depending on the types of CFD extrema, different preprocessing strategies on satellite clock products will be performed in the next section.

Detection of Small Extrema in CFD
The median value is better able to reveal the characteristics of most data than the mean value. The median value more effectively avoids the impact of outliers. CFD extrema detection is more improved while applying the method of median filter-based robust estimation [24]. Small ISB jumps,

Classification of CFD Extrema and Clock Preprocessing
After the analysis, it is clear that clock jumps are caused by many reasons. For each of the reasons, the application of a different clock preprocessing strategy is necessary. These reasons are categorized after the detection of these CFD extrema. In the time domain, small PMs are the reasons for one kind of clock jump. They are not detected in the satellite clock offsets series. As PMs do not affect ISB continuity, preprocessing is not implemented for these types of extrema. The day boundaries of the clock caused by the daily solution strategy for the precise clock products result in the other kind of clock jump. The best solution to these types of extrema is to extend the estimation time period of satellite clock, e.g., one week or even one month. For these types of extrema, satellite clock offsets are not re-estimated in a long period in this discussion. Instead, two different strategies are implemented. Polynomial fitting using the previous data is utilized when the clock jump is over threshold value. We employ the clock data from the previous day as the essential data to fit a polynomial, then use this polynomial to predict the clock data for the next day. After that, the original satellite clock will be replaced with this predicted clock data to get a new clock product. If not, the original clock product is kept, as well as no processing implemented. The switch of the reference clock causes the last kind of clock jump. For these types of extrema, it is also solved by implementing the method of polynomial fitting. Figure 6 illustrates the categorized types of extrema within CFD. We can see from Figure 6 that three different reasons mainly cause the extrema in CFD, which are: switches of the reference clock, clock day boundaries, and small PMs. When the reference clock is converted, extrema (jumps) will appear in each satellite. For these types of extrema, the magnitude is remarkable, so a method to generate a new satellite clock product by applying polynomial fitting is employed. In the current scenario of clock day boundaries, clock continuity will not be destroyed if the error is less than the threshold value, we apply no processing. Otherwise, the same polynomial fitting preprocessing procedure that is used for reference clock switches is utilized. The source of this threshold value is also from Equation (9), but in this part n is selected as 3. In the case of the CFD extrema detection, the median filter is applied to detect the hidden smaller PMs in the time domain of satellite clock offsets. However, continuity of the estimated ISB is not affected by those small clock PMs, so we also make no processing.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 19 A smooth and new clock product is gained once the preprocessing of satellite clock data is completed. After that, the ISB is re-estimated using the preprocessed satellite clock product, Figure 7 shows the new estimated ISB. Compared to the original ISB in Figure 2, the re-estimated ISB from this newly precise clock product has almost no jump, which is much more continuous and smoother.  A smooth and new clock product is gained once the preprocessing of satellite clock data is completed. After that, the ISB is re-estimated using the preprocessed satellite clock product, Figure 7 shows the new estimated ISB. Compared to the original ISB in Figure 2, the re-estimated ISB from this newly precise clock product has almost no jump, which is much more continuous and smoother.

Improvement of the ISB Continuity
For validating ISB continuity by applying the preprocessed satellite clock product, the statistical hypothesis technique is applied. The Equation (10) shows the statistics of the test (Zeng et al. 2017): where T N represents the continuity factor. Variance on days a and b is denoted as σ 2 a and σ 2 b . The weighted ISB mean values on days a and b are x a and x b , which can be calculated from: σ i 2 is the variance of ISB, and x i is the estimated ISB. The daily weighted mean value for ISB is expressed by x. Usually, N α (0, 1) means the critical value T N follows the normal distribution with the significance level α. If the equation T N > N α (0, 1) is satisfied, the ISB is discontinuous in two adjacent days, which means it varies greatly, or else ISB is continuous with no remarkable difference for two days. We can note that if the value of ISB continuity index T N is lower, the smoothness of ISB in two adjacent days is higher.
18 August and 20 August precise clock products were preprocessed previously in the above section. Thus, the test on 17-21 August was chosen here to analyze ISB continuity and its improvement.  Figure 9 shows the new ISB series after applying with the preprocessed precise clock products, along with the daily mean values. It is clear that the ISB discontinuity on 18 and 20 August was removed by the method of satellite clock preprocessing. adjacent days, which means it varies greatly, or else ISB is continuous with no remarkable difference for two days. We can note that if the value of ISB continuity index N T is lower, the smoothness of ISB in two adjacent days is higher.
18 August and 20 August precise clock products were preprocessed previously in the above section. Thus, the test on 17-21 August was chosen here to analyze ISB continuity and its improvement. ISB series for HKSL, HKWS, HKOH, JFNG, SIN1, CUT0, NNOR, MAJU, and NAUR stations are illustrated in the top panel of Figure 8. Taking MAJU, CUT0 and HKOH stations as an example, their mean values and standard deviations are illustrated by the panel at the bottom. We can find from Figure 8 that on 18 and 20 August, the top panel shows the discontinuity of the ISB series. The mean of daily ISB values from the bottom panel further confirms this finding. Figure 9 shows the new ISB series after applying with the preprocessed precise clock products, along with the daily mean values. It is clear that the ISB discontinuity on 18 and 20 August was removed by the method of satellite clock preprocessing.   We designed two schemes to reveal the improvement of ISB continuity by the satellite clock preprocessing. In scheme a, the solution is using the original clock product. In scheme b, the results are applying with preprocessed clock products. The continuity factors for both schemes a and b are determined by Equation (10). The orange line in Figure 10 shows the calculated value of significant level 0.01 in normal distribution, which is 2.327. We designed two schemes to reveal the improvement of ISB continuity by the satellite clock preprocessing. In scheme a, the solution is using the original clock product. In scheme b, the results are applying with preprocessed clock products. The continuity factors for both schemes a and b are determined by Equation (10). The orange line in Figure 10 shows the calculated value of significant level 0.01 in normal distribution, which is 2.327. We designed two schemes to reveal the improvement of ISB continuity by the satellite clock preprocessing. In scheme a, the solution is using the original clock product. In scheme b, the results are applying with preprocessed clock products. The continuity factors for both schemes a and b are determined by Equation (10). The orange line in Figure 10 shows the calculated value of significant level 0.01 in normal distribution, which is 2.327. Figure 10. The values of the continuity factor for the two designed schemes (a, b). Three different signs are used to express MAJU, CUT0, and HKOH stations. The critical value for normal distribution with the significant value of 0.01 is expressed by the orange line.
ISB is regarded as continuous and stable in the adjacent days if the continuity factor is less than the critical value (2.327, the orange line). In scheme a, ISB in each adjacent day during the period of ISB is regarded as continuous and stable in the adjacent days if the continuity factor is less than the critical value (2.327, the orange line). In scheme a, ISB in each adjacent day during the period of August 17-21 is discontinuous, as the continuity factor is all above the line of critical value. On the contrary, for scheme b, all the ISB continuity factors between adjacent days of August 17-21 are lower than the critical value, which means after satellite clock products are preprocessed, ISB continuity can be markedly improved. Table 7 includes the detailed statistics of the continuity factor. It is seen that, the same as in Figure 10, the ISB continuity can be significantly increased using scheme b. For day-pairs 17/18, 18/19, 19/20 and 20/21 August, the improvement rates are 62.4%, 79.5%, 99.6% and 99.0%. It can be concluded that it is possible to remarkably improve ISB continuity by applying preprocessing on precise satellite clock products.

Conclusions
In the time domain, large PMs as well as DIs in satellite clock products can be detected easily. However, it was found that PMs in satellite clock product were not the true cause for ISB jumps. Again, in the time domain, the linear correlation in satellite clock offsets and ISB was discussed. From the numerical results, a strong positive correlation was found between the residuals of satellite clock offset and ISB differences. Because of this, the method of making preprocessing on precise satellite clock products was proposed, which was used to remove the ISB discontinuity. Satellite clock data were converted from the time domain into the frequency domain. It made the detection of small jumps drowning in large-scale clock data much more sensitive. For detecting the clock jumps in the frequency domain data, a robust estimation method was applied, which was based on the median filter. After applying this method, we found that extrema in the data of the frequency domain were the real reason for ISB jumps. It was also found that clock day boundaries, reference clock switches, and small PMs were three individual reasons, which are categorized by the detection of smaller clock jumps. For all these three separate reasons, different strategies of clock preprocessing were implemented. After that, the re-estimation of ISB was performed using the preprocessed satellite clock product. From the ISB continuity improvement test, it was known that it was possible to improve the ISB continuity by 85.1% with the proposed preprocessing method on clock products. It showed the appropriateness of precise clock product preprocessing in terms of detecting and repairing the ISB jumps.
When we model ISB, ISB series need to be isolated into a number of components if the ISB series are discontinuous. But this problem can be overcome since we can apply the method of preprocessing on satellite clock products presented here. With this method, the smoother and continuous ISB series can be achieved. In that case, it is possible to obtain more sample points for modeling, as a result, the accuracy of modeling can also be enhanced. A future study on ISB modeling will be made using this method on the preprocessing of clock products.