GNSS/RNSS Integrated PPP Time Transfer: Performance with Almost Fully Deployed Multiple Constellations and a Priori ISB Constraints Considering Satellite Clock Datums

: Currently, the space segment of all the ﬁve satellite systems capable of providing precise time transfer services, namely BDS (including BDS-3 and BDS-2), GPS, GLONASS, Galileo, and Quasi-Zenith Satellite System (QZSS), has almost been fully deployed, which will deﬁnitely beneﬁt the precise time transfer with satellite-based precise point positioning (PPP) technology. This study focuses on the latest performance of the BDS/GPS/GLONASS/Galileo/QZSS ﬁve-system combined PPP time transfer. The time transfer accuracy of the ﬁve-system integrated PPP was 0.061 ns, and the frequency stability was 1.24 × 10 − 13 , 2.28 × 10 − 14 , and 8.74 × 10 − 15 at an average time of 10 2 , 10 3 , and 10 4 s, respectively, which signiﬁcantly outperforms the single-system cases. We also veriﬁed the outstanding time transfer performance of the ﬁve-system integrated PPP at locations with limited sky view. In addition, a method is proposed to mitigate the day-boundary jumps of inter-system bias (ISB) estimates by considering the difference in the satellite clock datums between two adjacent days. After applying a priori ISB constraints, the time transfer accuracy of the ﬁve-system integrated PPP can be improved by 37.9–51.6%, and the frequency stability can be improved by 14.8–21.6%, 5.3–7.6% and 20.0–29.6% at the three average times, respectively.


Introduction
In addition to the positioning and navigation services, the booming Global Navigation Satellite System (GNSS), which has the advantages of high accuracy, low cost, all-weather capability, and flexibility, can also provide precise time-related services (timing and remote time transfer) to users around the world.With the modernization of GNSS, the combination of observation data from multi-GNSS constellations can adequately utilize the on-orbit satellites, so that the accuracy and frequency stability of time transfer can be greatly improved.When the observation conditions are poor, the improvement brought by the multi-GNSS combination is considerable.
The GNSS time transfer methods can be divided into common view (CV) time transfer, all-in view (AV) time transfer, and precise point positioning (PPP) time transfer.The CV and AV methods are based on the pseudo-range observations [1], hence the time transfer accuracy is limited.The CV method strictly requires the inter-visibility of the two places for the same satellites, and the corresponding time transfer accuracy is only 1-10 ns.Following Zhang et al. [2], the accuracy of optimized CV time transfer could be merely 2.5 ns at a distance of 3000 km between two stations.The AV method does not require the intervisibility, and its time transfer accuracy can reach the nanosecond level.The PPP [3] time transfer has attracted increasing attention from many scholars, as it is independent of intervisibility, and employs high-accuracy carrier phase observations.Orgiazz et al. [4] assessed the performance of the GPS PPP time transfer, and the results showed that the PPP time transfer could be as accurate as the two-way satellite time and frequency transfer (TWSTFT).Kuang et al. [5] also assessed the performance of PPP time transfer based on the GPS data, and the result demonstrated that the frequency stability could reach 2 × 10 −14 within a day.Petit and Jiang [6] adopted the GPS PPP to realize the inter-station International Atomic Time (TAI) calculation between station UNSO and station PTB, and it was proved that the PPP would be fully capable of replacing the traditional time transfer technique.With the development and maturity of the GLONASS, Galileo, and BDS, many scholars focused on the application of other GNSSs in precise time transfer.Guang et al. [7] successfully realized the PPP time transfer based on BDS-2 dual-frequency observations, and their results showed that the frequency stability could reach 5.59 × 10 −14 when the average time was 15,360 s.Ge et al. [8] rigorously investigated the influence of Inter-Frequency Code Bias (IFCB) on GLONASS-only PPP, and established a feasible GLONASS-only PPP time transfer model.Zhang et al. [9] evaluated the enhanced performance of Galileo-only PPP time transfer with a priori constraints of receiver coordinates and zenith tropospheric delays.
With the gradual progress of the Multi-GNSS Experiment (MGEX), more and more scholars began to focus on the PPP time transfer performance of multi-GNSS combinations.Defraigne and Baire [10] adopted the PPP and CV methods for time transfer based on GPS/GLONASS dual-system observations.The results showed that the dual-system combination could improve the short-term frequency stability of both two methods as well as the time transfer accuracy of PPP method, but the accuracy of CV time transfer was not significantly affected.Wang et al. [11] evaluated the performance of real-time PPP time transfer based on the GPS, GLONASS, and Galileo triple-system combination.The results showed that the GPS-only case performed better than the other single-system cases, and the triple-system case had the best time transfer performance.Ge et al. [12] further assessed the performance improvement of the PPP time transfer based on the GPS/GLONASS/BDS-2/Galileo quad-system combination compared with the single-system cases, and the gain of PPP time transfer brought by a modified model considering the correlation of receiver clock offsets among epochs.The results showed that, compared with GPS-only, Galileo-only, and BDS-2-only cases, the frequency stability of the quad-system PPP time transfer could be improved by 20.3%, 45.4%, and 84.0%, respectively, and the application of the modified time transfer model could further improve the frequency stability.Zhang et al. [13] evaluated the impact of different post-processed precise products from three analysis centers, namely the Center for Orbit Determination in Europe (CODE), Wuhan University (WHU), and GeoForschungsZentrum Potsdam (GFZ), on the PPP time transfer performance of multi-GNSS combinations.The results showed that the PPP with precise products provided by CODE and GFZ had better time transfer performance than that with WHU.Tu et al. [14] also investigated the PPP time transfer performance of the GPS/GLONASS/BDS-2/Galileo quad-system combination and verified that constraining PPP with a priori inter-system bias (ISB) could improve the time transfer accuracy.
As of 31 July 2020, the construction of the full BDS-3 constellation was completed [15].The nominal BDS-3 constellation includes 3 geostationary orbit (GEO) satellites, 3 inclined geosynchronous orbit (IGSO) satellites, and 24 medium earth orbit (MEO) satellites [16].Some researchers evaluated the PPP time transfer performance of BDS-3 [17,18].It is worth noting that most of the atomic clocks carried by the BDS-3 satellites are hydrogen masers.Compared with the satellite clocks of the BDS-2 satellites, those of the BDS-3 satellites are more stable.In addition, the quality of BDS-3 observations is better.Therefore, it is generally believed that there is also an ISB between BDS-2 and BDS-3.Jiao et al. [19] adopted four strategies (i.e., estimating ISB as constants, random walk process, and white noise process, as well as ignoring ISB) to deal with the ISB between BDS-2 and BDS-3, and compared the influence of the four strategies on PPP time transfer performance.The result showed that under the four strategies, the frequency stability with the legacy signals B1I and B3I was improved by 23.89%, 23.96%, 11.46%, and 20.18%, respectively, compared with the time transfer results using only BDS-2 observations.For the frequency stability with B1C and B2a, which are newly broadcast by BDS-3, the corresponding enhancements of the four strategies were 20.22%, 17.53%, −3.69%, and −50.77%, respectively.This provides a foundation for the selection of ISB estimation strategies in BDS-2/BDS-3 PPP time transfer.
With the modernization of GPS, the gradual revival of GLONASS, the completion of the full constellation construction of BDS-3, the rapid progress of Galileo, and the vigorous development of Quasi-Zenith Satellite System (QZSS), more and more in-orbit GNSS and Regional Navigation Satellite System (RNSS) satellites that are capable of providing precise services are available and the observation quality is continuously improved.It is foreseeable that the performance of multi-system integrated PPP time transfer will be improved accordingly.The previous studies on PPP time transfer based on multi-constellation combinations could not fully utilize all available satellites, as the satellites were continuously launched with the modernization of GNSS and RNSS.The derived conclusions might be obsolete.In order to fully exploit the advantages of multi-constellation combinations in PPP time transfer under the current constellation status, this contribution focuses on the BDS/GPS/GLONASS/Galileo/QZSS five-system integrated PPP time transfer, in which the ISB parameter between BDS-2 and BDS-3 is carefully considered, and a rigorous multisystem PPP time transfer mathematical model is established.Although some previous studies tried to enhance the multi-system PPP time transfer performance by introducing the a priori ISBs into the observation model [14], the day-boundary jumps of ISB estimates due to the effects of satellite clock datums were ignored.Therefore, the five-system integrated PPP time transfer with a priori ISB constraints considering satellite clock datums is also investigated in this study.This paper starts with the materials and methods, including the five-system integrated PPP time transfer model and the proposed PPP time transfer model with a priori ISB constraints considering satellite clock datums.Subsequently, the results are presented and analyzed.Then, the discussion is conducted.Finally, the main conclusions are summarized.

Five-System Integrated PPP Time Transfer Model
In this contribution, the dual-frequency code and phase observations from GPS, GLONASS, Galileo, BDS (BDS-2 and BDS-3), and QZSS are used to form the ionosphericfree (IF) combination to eliminate the influence of first-order ionospheric delay, which can be expressed as: where s and r denote a GNSS (or RNSS) satellite and a receiver, respectively, Q is the navigation satellite system corresponding to the used satellite in the specific equation, the subscript IF denotes the dual-frequency IF combination, P and L are code and phase observations in meters, respectively, ρ is the geometric distance between satellite and receiver, cdt r and cdt s represent the receiver clock offset and satellite clock offset, respectively, T denotes the slant tropospheric delay, N is the phase ambiguity in meters, b r and b s represent the receiver-dependent and satellite-dependent code hardware delays, respectively, B r and B s represent the phase hardware delays at the receiver and satellite, respectively, and ε denotes the measurement noise including multipath error.The dual-frequency IF combined observables are usually used for precise clock estimation (PCE).As the frequency-independent physical clock errors of the satellite cannot be estimated in an unbiased manner because of rank deficiency, the satellite clock estimates can be formulated as: where dD AC is the satellite clock datum introduced by the respective analysis center (AC) during PCE.
In view that the time scales as well as the satellite clock datums of each GNSS (or RNSS) are defined differently, and the corresponding receiver-dependent code hardware delay differs from each other, it is necessary to introduce ISB parameters to compensate for the inconsistency.Besides, it is generally believed that the ISB between BDS-3 and BDS-2 should be considered due to the different quality of observations.This contribution takes the receiver clock offset of GPS as the datum in the five-system combined PPP time transfer model, and introduces ISB parameters into the observation equations for other satellite systems.Similarly, in the BDS-only observation model, the receiver clock offset of BDS-3 is used as the datum, and the ISB parameters are introduced into the observation equations of BDS-2.After applying the precise satellite orbit and clock offset products and using a priori model to correct the tropospheric dry delay, the linearized five-system PPP time transfer model can be described as follows: where G, R, E, C2, C3, and J denote the GPS, GLONASS, Galileo, BDS-2, BDS-3, and QZSS, respectively, p and l denote the observed-minus-computed (OMC) values of code and phase observations, respectively, µ denotes the unit vector in the direction from the receiver to the satellite, x represents the three-dimensional (3D) coordinates of the receiver, cdt G r is the receiver clock estimate including the receiver-dependent code hardware delay as well as the satellite clock datum of GPS, ISB G,Q r denotes the ISB between GPS and other satellite systems, N s,Q IF,r is the phase ambiguity that loses its integer character after being mixed with various hardware delays, m indicates the wet mapping function, and Z denotes the tropospheric zenith wet delay (ZWD).
To enhance the PPP time transfer, the station coordinates in the Solution Independent Exchange (SINEX) weekly solution files provided by the International GNSS Service (IGS) can be used as a priori constraints.The following pseudo-observation equation can be introduced into the observation model of PPP time transfer: where x SNX denotes the a priori receiver coordinates in three dimensions.The parameters that need to be estimated include the 3D coordinates of the receiver, receiver clock offset, ISB, tropospheric ZWD, and phase ambiguity, which can be expressed as: where O denotes the vector of estimates.
In addition to a rigorous functional model, a suitable stochastic model is also crucial for high-accuracy PPP time transfer.Assuming that the observations from different satellites or on different frequencies are independent of each other, an elevation-dependent strategy is employed to determine the weights, which can be expressed as follows: where σ 2 (ele) is the variance of observations at the elevation angle ele, k 1 and k 2 are constants, which can be set to 3 mm for phase observations and 3 dm for code observations, and a represents the weight factor.The weight factor does not remain consistent for all the GNSS and RNSS constellations due to different observation qualities.The weight factors of GPS, GLONASS, Galileo, BDS, and QZSS are set to 1.0, 1.5, 1.0, 1.0, and 2.0, respectively.The GLONASS and QZSS observations are down-weighted due to the complex influence of frequency division multiple access (FDMA) modulation scheme or the relatively lower accuracy of precise products.In addition, because of the inferior orbit accuracy of the BDS GEO and IGSO satellites, the corresponding factors are multiplied by 10 times and 2 times when determining the weights of their observations.As for the stochastic model for the pseudo-observation, it can be determined as: where σ 2 x is the variance of a priori 3D receiver coordinates, and k 3 is also a constant.Although the receiver coordinates released by the IGS usually can reach an accuracy level of a few millimeters, the numerical value of k 3 is taken as 2 cm here, so as to avoid the biased estimates.

PPP Time Transfer Model with a Priori ISB Constraints Considering Satellite Clock Datums
In view that the ISB estimates are usually very stable, it is deduced that the a priori ISB estimates, such as the estimates of the previous day, can be used as constraints to enhance the PPP time transfer performance.Following Equation (4), the ISB estimates are composed of the differences in time scales, satellite clock datums, and receiver-dependent code hardware delays among different satellite systems.As the time system is aligned to the GPS time for all satellite systems in both multi-system observation files and precise satellite clock files, the effects of time scales are marginal and can be ignored.As for the receiver-dependent code hardware delays, they are considered to be stable within a short period of time, for example, several days [20,21].In contrast, the satellite clock datums are redefined day-by-day.Thus, the different satellite clock datums between two adjacent days will result in the distinct ISB estimates for the two days, which is an obstacle for the introduction of a priori ISB constraints in the data processing of multi-system integrated PPP time transfer.If the variations of the satellite clock datums between two adjacent days can be computed in advance, this problem will be solved.
With the use of the satellite clock corrections at the GPS time 23:58:30, 23:59:00, and 23:59:30 in the precise satellite clock file of the previous day, those at the GPS time 00:00:00 on this day can be predicted based on the second-order polynomial fitting.The factors that represent the stability of satellite clocks in the time-frequency domain include phase, frequency, and frequency drift.Usually, a second-order polynomial model containing the three factors can be selected as the clock offset prediction model [22].The complicated periodic terms can be ignored with the datasets of only several minutes.Due to the effects of different satellite clock datums, the satellite clock corrections at the GPS time 00:00:00 in the precise satellite clock file of this day differ from the predicted ones.However, the difference between them is roughly consistent for all the satellites (see Equation ( 2)).It is noted that the satellite-dependent code hardware delays are also stable within a short period of time.Thus, the difference in the satellite clock datums between two adjacent days can be calculated as: where dD AC,Q t is the satellite clock datum of this day, dD AC,Q p is the satellite clock datum of the previous day, cdt s,Q IF,e is the satellite clock corrections at the GPS time 00:00:00 in the precise satellite clock file of this day, cdt s,Q IF,p denotes the predicted ones, and n is the number of satellites.
To constrain the ISB estimates with the consideration of satellite clock datums, the following pseudo-observation equation can be introduced into the observation model of multi-system integrated PPP time transfer: where ISB G,Q r,p denotes the fully converged ISB estimates of the previous day.As to the stochastic model for the pseudo-observation, the STD of a priori ISBs including the adjustment of satellite clock datums is empirically set to 1 ns.

Data Sets and Processing Strategies
In order to assess the PPP time transfer performance of the BDS/GPS/GLONASS/ Galileo/QZSS five-system combination, we selected 13 MGEX stations around the world, which are all equipped with high-performance hydrogen atomic clocks, for analysis.The station BRUX operated by the time laboratory was chosen as the central node to form 12 time links with other stations.Since the orbital repeat period is one day for GPS, QZSS, and BDS non-MEO satellites and several days for other satellites, the employed data sets span seven days (also taken in [23,24]) from day of year (DOY) 248 to 254 in 2021.The geographical distribution of the selected 13 MGEX stations is shown in Figure 1, and Table 1 lists the detailed information for these stations.Table 2 summarizes the processing strategy of five-system integrated PPP time transfer.

Availability of Five-System Integration
Based on the number of theoretically visible satellites, we assessed the availability of the BDS/GPS/GLONASS/Galileo/QZSS five-system combination.Following Yang et al. [25], the whole world can be evenly divided into 72 × 72 grids; each of them extends for 2.5 • and 5 • in terms of latitude and longitude, respectively, and it is assumed that there is a virtual receiver at an altitude of 25 m in the center of each grid.The receiver sampling interval is set to 30 s, and the elevation cut-off angle is set to 10 • .With the use of the satellite coordinates provided by GFZ precise products, the theoretical availability of five-system combination can be assessed by analyzing the relative geometry of each satellite and the virtual receiver.
The theoretical distribution of the average number of visible satellites (over a day) under the five-system combination case on DOY 254, 2021, is shown in Figure 2. Obviously, a diamond-shaped area exists in the figure.It centers around the point at 120 • E longitude at the equator, and extends from 60 • E to 180 • in the east-west direction and from 60 • N to 60 • S in north-south direction.The number of visible satellites in the diamond-shaped regions far exceeds other regions, with an average of 48.7.In contrast, there is also a typical region extending from 165 • W to 45 • E in the east-west direction, and from 60 • N to 60 • S in the north-south direction.The number of visible satellites is significantly reduced in this region, with an average number of 33.9.It is clear that the number of visible satellites decreases with the increasing distance from the diamond-shaped area, which is caused by the GEO and the IGSO satellites.Under ideal conditions, the average number of visible satellites on a global scale can be 38.9, with a minimum statistic of 29.7, implying that the reliable PPP time transfer can be ensured even at locations with limited sky view.
The theoretical distribution of the average number of visible satellites (over under the five-system combination case on DOY 254, 2021, is shown in Figure 2 ously, a diamond-shaped area exists in the figure.It centers around the point a longitude at the equator, and extends from 60°E to 180° in the east-west directi from 60°N to 60°S in north-south direction.The number of visible satellites in t mond-shaped regions far exceeds other regions, with an average of 48.7.In contras is also a typical region extending from 165°W to 45°E in the east-west direction, an 60°N to 60°S in the north-south direction.The number of visible satellites is signif reduced in this region, with an average number of 33.9.It is clear that the number of satellites decreases with the increasing distance from the diamond-shaped area, w caused by the GEO and the IGSO satellites.Under ideal conditions, the average n of visible satellites on a global scale can be 38.9, with a minimum statistic of 29.7, im that the reliable PPP time transfer can be ensured even at locations with limited sk For further analysis, the actual number of visible satellites and time dilution cision (TDOP) values at stations BRUX, ROAG, and SPT0 on DOY 254, 2021, are pre in Figure 3.The smaller TDOP values are beneficial for the high-accuracy time tran addition to the five-system integration (GRECJ), the GPS-only (GPS), GLONAS (GLO), Galileo-only (GAL), and BDS-only (BDS) cases are also provided in Figu terms of the actual satellite number, it is generally less than 15 for each GNSS For further analysis, the actual number of visible satellites and time dilution of precision (TDOP) values at stations BRUX, ROAG, and SPT0 on DOY 254, 2021, are presented in Figure 3.The smaller TDOP values are beneficial for the high-accuracy time transfer.In addition to the five-system integration (GRECJ), the GPS-only (GPS), GLONASS-only (GLO), Galileo-only (GAL), and BDS-only (BDS) cases are also provided in Figure 3.In terms of the actual satellite number, it is generally less than 15 for each GNSS single-system case, while the average satellite number at stations BRUX, ROAG, and SPT0 under the five-system combination is 32.0, 29.7, and 36.4,respectively.For the TDOP, the average value at three stations is 0.87, 0.96, and 0.86 for the GPS-only case, 1.53, 1.83, and 1.82 for the GLONASS-only case, 1.53, 1.35, and 1.25 for the Galileo-only case, and 0.91, 0.99, and 0.71 for the BDS-only case, respectively.Notably, the TDOP value of the GLONASS-only and Galileo-only cases is extremely unstable and has a wide range of variation due to the relatively smaller number of available satellites.In contrast, the average TDOP value at three stations under the five-system combination case is 0.42, 0.44, and 0.40, respectively.Thus, compared with the GNSS single-system cases, the five-system combination can significantly increase the number of visible satellites, reduce the TDOP values, and improve the stability of the time series of epoch-wise TDOP values, which will definitely benefit the PPP time transfer.In addition, the average number of visible satellites for the five-system combination at station TID1 on DOY 254 of 2021 can be 42.3.In conjunction with the results shown in Figure 2 as well as the information of stations provided in Table 1, it is indicated that the actual satellite numbers are roughly close to the theoretical ones.This verifies the correctness of our theoretical analysis about the availability of five-system integration.the stability of the time series of epoch-wise TDOP values, which will definitely benefit the PPP time transfer.In addition, the average number of visible satellites for the fivesystem combination at station TID1 on DOY 254 of 2021 can be 42.3.In conjunction with the results shown in Figure 2 as well as the information of stations provided in Table 1, it is indicated that the actual satellite numbers are roughly close to the theoretical ones.This verifies the correctness of our theoretical analysis about the availability of five-system integration.

Performance of the Five-System PPP Time Transfer under Open Sky
To assess the remote PPP time transfer performance of the five-system combination under open sky, the elevation cut-off angle is set to 10°.For comparison, the results of BDS-only, GPS-only, GLONASS-only, and Galileo-only PPP time transfer are also provided.The IGS final products, which are generated by combining the GPS-only or GPS/GLONASS network solutions of all analysis centers, are used as the reference.Figure 4 shows the difference in the receiver clock offset between PPP estimates and IGS final products for the five-system combination case and four GNSS single-system cases at stations BRUX and SPT0 on DOY 252, 2021.The effects of different time datums are eliminated by subtracting the average value.The results in the first two hours are excluded to ensure that only the fully converged results can be involved in the analysis.It is seen that

Performance of the Five-System PPP Time Transfer under Open Sky
To assess the remote PPP time transfer performance of the five-system combination under open sky, the elevation cut-off angle is set to 10 • .For comparison, the results of BDS-only, GPS-only, GLONASS-only, and Galileo-only PPP time transfer are also provided.The IGS final products, which are generated by combining the GPS-only or GPS/GLONASS network solutions of all analysis centers, are used as the reference.Figure 4 shows the difference in the receiver clock offset between PPP estimates and IGS final products for the five-system combination case and four GNSS single-system cases at stations BRUX and SPT0 on DOY 252, 2021.The effects of different time datums are eliminated by subtracting the average value.The results in the first two hours are excluded to ensure that only the fully converged results can be involved in the analysis.It is seen that the receiver clock difference usually varies within a range of 1 ns.In addition, the receiver clock difference exhibits some consistent trends at the two stations, which may be attributed to the commonly employed precise satellite products.Figure 5 further displays the corresponding difference of the time link BRUX-SPT0.The five-system combination case has the smallest time link difference, followed by the GPS-only and Galileo-only cases, while the varying range of the time link difference can be 0.25 ns for the BDS-only case and 0.5 ns for the GLONASS-only case.Benefiting from the strong satellite geometry and the considerable number of visible satellites under an open-sky environment, as well as the high-accuracy satellite orbit and clock products, the observation information is sufficient even for GPS-only and Galileo-only cases.Thus, the improvement on the time link difference for the five-system combination case over the GPS-only and Galileo-only cases is not significant.
Usually, the standard deviation (STD) statistics of the time link differences between PPP estimates and IGS final products can be taken as an important index to evaluate the PPP time transfer accuracy.We first calculated the STD statistics of the daily time link differences, and then obtained the average STD statistics over a week.Table 3 illustrates the derived results of four different time links for the five-system combination case and the four GNSS single-system cases.The time transfer accuracy of the five-system integrated PPP is better than that of the GNSS single-system PPP for most of the cases.According to the average results over the four time links, the time transfer accuracy is 0.074 ns for the GPS-only PPP, 0.359 ns for the GLONASS-only PPP, 0.068 ns for the Galileo-only PPP, and 0.346 ns for the BDS-only PPP.The five-system integrated PPP can improve the time transfer accuracy by 18% over the GPS-only case, 83% over the GLONASS-only case, 10% over the Galileo-only case, and 82% over the BDS-only case, to 0.061 ns.
the receiver clock difference usually varies within a range of 1 ns.In addition, the receiver clock difference exhibits some consistent trends at the two stations, which may be attributed to the commonly employed precise satellite products.Figure 5 further displays the corresponding difference of the time link BRUX-SPT0.The five-system combination case has the smallest time link difference, followed by the GPS-only and Galileo-only cases, while the varying range of the time link difference can be 0.25 ns for the BDS-only case and 0.5 ns for the GLONASS-only case.Benefiting from the strong satellite geometry and the considerable number of visible satellites under an open-sky environment, as well as the high-accuracy satellite orbit and clock products, the observation information is sufficient even for GPS-only and Galileo-only cases.Thus, the improvement on the time link difference for the five-system combination case over the GPS-only and Galileo-only cases is not significant.Usually, the standard deviation (STD) statistics of the time link differences between PPP estimates and IGS final products can be taken as an important index to evaluate the PPP time transfer accuracy.We first calculated the STD statistics of the daily time link differences, and then obtained the average STD statistics over a week.Table 3 illustrates the receiver clock difference usually varies within a range of 1 ns.In addition, the receiver clock difference exhibits some consistent trends at the two stations, which may be attributed to the commonly employed precise satellite products.Figure 5 further displays the corresponding difference of the time link BRUX-SPT0.The five-system combination case has the smallest time link difference, followed by the GPS-only and Galileo-only cases, while the varying range of the time link difference can be 0.25 ns for the BDS-only case and 0.5 ns for the GLONASS-only case.Benefiting from the strong satellite geometry and the considerable number of visible satellites under an open-sky environment, as well as the high-accuracy satellite orbit and clock products, the observation information is sufficient even for GPS-only and Galileo-only cases.Thus, the improvement on the time link difference for the five-system combination case over the GPS-only and Galileo-only cases is not significant.Usually, the standard deviation (STD) statistics of the time link differences between PPP estimates and IGS final products can be taken as an important index to evaluate the PPP time transfer accuracy.We first calculated the STD statistics of the daily time link differences, and then obtained the average STD statistics over a week.Table 3 illustrates Table 3.Average STD statistics of the time link differences over a week for the five-system combination case and the four GNSS single-system cases (unit: ns).All the time links involved in this contribution take the station BRUX as the center node, but the two stations that constitute the time link have different time and frequency references.As such, it is difficult to directly compare the PPP time transfer performance on stability.Therefore, this contribution adopts the Modified Allan Deviation (MDEV) to evaluate the frequency stability of PPP time transfer.Figure 6 illustrates the MDEV values of three time links, BRUX-MGUE, BRUX-PTBB, and BRUX-SPT0, using a cut-off elevation angle of 10 • on DOY 252, 2021.The frequency stability of the GLONASS-only and BDS-only PPP time transfer is obviously inferior to other time transfer cases, and the MDEV statistics of the five-system combination case are superior to the four single-GNSS cases.When the average time is larger than 2 × 10 4 s, the MDEV values are not reliable, which can be attributed to the small number of data samples, and thus they are not presented in the figure .BRUX All the time links involved in this contribution take the station BRUX as the center node, but the two stations that constitute the time link have different time and frequency references.As such, it is difficult to directly compare the PPP time transfer performance on stability.Therefore, this contribution adopts the Modified Allan Deviation (MDEV) to evaluate the frequency stability of PPP time transfer.Figure 6 illustrates the MDEV values of three time links, BRUX-MGUE, BRUX-PTBB, and BRUX-SPT0, using a cut-off elevation angle of 10° on DOY 252, 2021.The frequency stability of the GLONASS-only and BDSonly PPP time transfer is obviously inferior to other time transfer cases, and the MDEV statistics of the five-system combination case are superior to the four single-GNSS cases.When the average time is larger than 2 × 10 4 s, the MDEV values are not reliable, which can be attributed to the small number of data samples, and thus they are not presented in the figure.To investigate the effects of the selected datum of the receiver clock offset on the five-system combined PPP time transfer performance, we compared the MDEV values of the time link BRUX-MGUE taking the receiver clock offset of GPS and BDS-3 as the datum on DOY 252, 2021.The frequency stability for the five-system combined PPP time transfer with GPS datum can be 1.14 × 10 −13 , 1.70 × 10 −14 , and 9.88 × 10 −15 at an average time of 10 2 , 10 3 , and 10 4 s, respectively, while the corresponding frequency stability for that with BDS-3 datum can be 1.25 × 10 −13 , 1.48 × 10 −14 , and 9.99 × 10 −15 , respectively.The difference in the MDEV values between the two datums is not significant.Therefore, the receiver clock offset of GPS is taken as the datum in all other five-system combined PPP time transfer results.
Figures 7, 8 and 9 illustrate the MDEV values of three time links, BRUX-MGUE, BRUX-PTBB, and BRUX-SPT0, on DOY 248 to 254, 2021, when the average time is 10 2 , 10 3 , and 10 4 s, respectively.The MDEV values at the average time 10 2 , 10 3 , and 10 4 s are taken as typical results to reflect the short-term, medium-term, and long-term frequency stability, respectively [26].It is seen that the GLONASS-only PPP time transfer shows the largest day-to-day scattering of the MDEV values, followed by the Galileo-only and BDS-only cases, which may be attributed to the relatively lower accuracy of precise satellite products or the insufficient number of available satellites.In contrast, the day-to-day scattering of the MDEV values is not obvious for GPS-only and five-system integrated PPP time transfer, especially for the five-system combination case.For further analysis, Table 4 provides the average value of day-to-day variation ratios of MDEV statistics over seven days (DOY 248 to 254, 2021) and three specific average times (10 2 , 10 3 , and 10 4 s) for the three employed time links.According to the average results over the three time links, the variation ratios of frequency stability between two adjacent days for the GLONASS-only PPP time transfer can be up to 30.5%, while the corresponding variation ratios under the GPS, Galileo, and BDS single-GNSS cases are 7.7%, 10.3%, and 10.0%, respectively.As for the five-system combination case, the average variation ratios of frequency stability are only 4.8%, implying its ability to provide a steady time transfer service.
as typical results to reflect the short-term, medium-term, and long-term frequency stability, respectively [26].It is seen that the GLONASS-only PPP time transfer shows the largest day-to-day scattering of the MDEV values, followed by the Galileo-only and BDS-only cases, which may be attributed to the relatively lower accuracy of precise satellite products or the insufficient number of available satellites.In contrast, the day-to-day scattering of the MDEV values is not obvious for GPS-only and five-system integrated PPP time transfer, especially for the five-system combination case.For further analysis, Table 4 provides the average value of day-to-day variation ratios of MDEV statistics over seven days (DOY 248 to 254, 2021) and three specific average times (10 2 , 10 3 , and 10 4 s) for the three employed time links.According to the average results over the three time links, the variation ratios of frequency stability between two adjacent days for the GLONASS-only PPP time transfer can be up to 30.5%, while the corresponding variation ratios under the GPS, Galileo, and BDS single-GNSS cases are 7.7%, 10.3%, and 10.0%, respectively.As for the five-system combination case, the average variation ratios of frequency stability are only 4.8%, implying its ability to provide a steady time transfer service.Figure 10 displays the average MDEV values over a week for 12 time links at an average time of 10 2 , 10 3 , and 10 4 s.It is seen that the difference in the MDEV statistics among different time links is most obvious for the GLONASS-only case, followed by the Galileo-only and BDS-only cases, which may be attributed to the distinct number of visible satellites for the stations located in different latitude and longitude regions.In contrast, the MDEV values of the five-system combination case at an average time of 10 2 and 10 3 s seem to be at the same level for different time links.Except for GLONASS-only PPP time transfer, the average MDEV values for most of the time links are on the order of 10 −13 when the average time is 10 2 s; the corresponding values are on the order of 10 −14 when the average time increases to 10 3 s; and the corresponding values are on the order of 10 −14 to 10 −15 when the average time is as long as 10 4 s.As for the GLONASS-only case, the average MDEV values are usually on the order of 10 −12 to 10 −13 , 10 −13 to 10 −14 , and 10 −13 to 10 −14 when the average time is 10 2 , 10 3 , and 10 4 s, respectively.to be at the same level for different time links.Except for GLONASS-only PPP time transfer, the average MDEV values for most of the time links are on the order of 10 −13 when the average time is 10 2 s; the corresponding values are on the order of 10 −14 when the average time increases to 10 3 s; and the corresponding values are on the order of 10 −14 to 10 −15 when the average time is as long as 10 4 s.As for the GLONASS-only case, the average MDEV values are usually on the order of 10 −12 to 10 −13 , 10 −13 to 10 −14 , and 10 −13 to 10 −14 when the average time is 10 2 , 10 3 , and 10 4 s, respectively.Figure 11 further shows the average MDEV results over a week and all 12 time links at an average time of 10 2 , 10 3 , and 10 4 s.The frequency stability of the five-system integrated PPP time transfer can be 1.24 × 10 −13 , 2.28 × 10 −14 , and 8.74 × 10 −15 at an average time of 10 2 , 10 3 , and 10 4 s, respectively.Compared with the GPS-only case, GLONASS-only case, Galileo-only case, and BDS-only case, the average improvement rates of the fivesystem combination case on the frequency stability are 18.4%, 82.2%, 34.7%, and 54.1%   Figure 11 further shows the average MDEV results over a week and all 12 time links at an average time of 10 2 , 10 3 , and 10 4 s.The frequency stability of the five-system integrated PPP time transfer can be 1.24 × 10 −13 , 2.28 × 10 −14 , and 8.74 × 10 −15 at an average time of 10 2 , 10 3 , and 10 4 s, respectively.Compared with the GPS-only case, GLONASS-only case, Galileo-only case, and BDS-only case, the average improvement rates of the five-system combination case on the frequency stability are 18.4%, 82.2%, 34.7%, and 54.1% when the average time is 10 2 s; the corresponding improvement rates are 19.4%,94.0%, 62.6%, and 71.6% when the average time is 10 3 s; and the average improvement rates are 21.3%, 85.9%, 23.3%, and 59.7% when the average time is 10 4 s.It is concluded that even under an open-sky environment, the five-system combination case can significantly improve the time transfer performance (in terms of the frequency stability as well as the time transfer accuracy) compared with the single-GNSS cases.

Performance of Five-System PPP Time Transfer at Locations with Limited Sky View
In the previous section, we assess the PPP time transfer performance under good observation conditions, but in the practical applications, the users also have the requirement of high-accuracy time transfer at locations with limited sky view.Therefore, this section focuses on the PPP time transfer performance of the BDS/GPS/GLONASS/Galileo/QZSS five-system combination under the environments with limited satellite visibility.To simply simulate the environments with limited sky view, the cut-off elevation angles are set from 20° to 40° in steps of 10°.For comparative analysis, the results of the four single-GNSS cases are also given.
Table 5 shows the service rates of each PPP time transfer case under observation conditions with limited sky view.The service rate for a time link is roughly defined as the ratio of the number of epochs in which the receiver clock differences between two stations can be successfully solved over the total number of epochs during the whole analysis period.The average results over all the 12 time links are computed and provided in Table 5.Under an elevation cut-off angle of 20°, the seven-day service rate of the GPS-only case, GLONASS-only case, Galileo-only case, and BDS-only case is 100.0%,70.6%, 79.3%, and 91.3%, respectively.After further setting the elevation cut-off angle to 30°, the corresponding service rate of the GLONASS-only case, Galileo-only case, and BDS-only case significantly decreases.When the elevation cut-off angle is further set to 40°, the service rate is

Performance of Five-System PPP Time Transfer at Locations with Limited Sky View
In the previous section, we assess the PPP time transfer performance under good observation conditions, but in the practical applications, the users also have the requirement of high-accuracy time transfer at locations with limited sky view.Therefore, this section focuses on the PPP time transfer performance of the BDS/GPS/GLONASS/Galileo/QZSS five-system combination under the environments with limited satellite visibility.To simply simulate the environments with limited sky view, the cut-off elevation angles are set from 20 • to 40 • in steps of 10 • .For comparative analysis, the results of the four single-GNSS cases are also given.
Table 5 shows the service rates of each PPP time transfer case under observation conditions with limited sky view.The service rate for a time link is roughly defined as the ratio of the number of epochs in which the receiver clock differences between two stations can be successfully solved over the total number of epochs during the whole analysis period.The average results over all the 12 time links are computed and provided in Table 5.Under an elevation cut-off angle of 20 • , the seven-day service rate of the GPS-only case, GLONASS-only case, Galileo-only case, and BDS-only case is 100.0%,70.6%, 79.3%, and 91.3%, respectively.After further setting the elevation cut-off angle to 30 • , the corresponding service rate of the GLONASS-only case, Galileo-only case, and BDS-only case significantly decreases.When the elevation cut-off angle is further set to 40 • , the service rate is reduced to less than 36% for the GPS-only case, and less than 6% for the other three single-GNSS cases.In contrast, the service rate of the five-system combination case can still maintain 100.0%.For further analysis, the service rates of single-GNSS PPP time transfer cases for each time link are detailed in Figure 12.It is seen that the service rates dramatically vary among the 12 time links.Actually, the low service rates of single-GNSS PPP time transfer cases are related to the number of visible satellites.
Table 5. Service rates of each PPP time transfer case using the cut-off elevation angle of 20  Table 6 further provides the success rates of each PPP time transfer case under observation conditions with limited sky view.The success rate is defined as the ratio of the number of epochs with successfully solved receiver clock differences over the number of epochs with sufficient satellites (at least four satellites) using the datasets of all the stations and days.Similar to the service rates (see Table 5), the success rates always maintain 100.0%for the five-system combination case, and decrease with the increasing cut-off elevation angles for the four single-GNSS cases.When using the cut-off elevation angle of 40°, the success rates are still larger than 90% for the GPS-only, Galileo-only, and BDSonly cases, while those of the GLONASS-only case are reduced to only 37.6%, which may be attributed to the worse quality of observations and the lower accuracy of precise satellite products in comparison with other satellite systems.
Table 6.Success rates of each PPP time transfer case using the cut-off elevation angle of 20°, 30°, and 40° (unit: %).Table 6 further provides the success rates of each PPP time transfer case under observation conditions with limited sky view.The success rate is defined as the ratio of the number of epochs with successfully solved receiver clock differences over the number of epochs with sufficient satellites (at least four satellites) using the datasets of all the stations and days.Similar to the service rates (see Table 5), the success rates always maintain 100.0%for the five-system combination case, and decrease with the increasing cut-off elevation angles for the four single-GNSS cases.When using the cut-off elevation angle of 40 • , the success rates are still larger than 90% for the GPS-only, Galileo-only, and BDS-only cases, while those of the GLONASS-only case are reduced to only 37.6%, which may be attributed to the worse quality of observations and the lower accuracy of precise satellite products in comparison with other satellite systems.

Elevation
Table 6.Success rates of each PPP time transfer case using the cut-off elevation angle of 20 Figure 13 shows the difference in the receiver clock offset between PPP estimates and IGS final products for five PPP time transfer cases at stations BRUX and SPT0 under an elevation cut-off angle of 20 • on DOY 252, 2021, and Figure 14 displays the corresponding difference of the time link BRUX-SPT0.Obviously, when the elevation cut-off angle is increased to 20 • , the performance of the receiver clock offset estimates and time link estimates of each single-GNSS case is worse than that under an elevation cut-off angle of 10 • .The fluctuation range gets larger for the errors of receiver clock offset estimates and time link estimates, and the results at some epochs are absent.By contrast, the receiver clock offset difference and time link difference of the five-system combination case are similar to the results at a mask elevation of 10 to the results at a mask elevation of 10°, and show good performance in terms of data integrity.Tables 7 and 8 provide the average STD statistics of time link differences over a week under an elevation cut-off angle of 20° and 30°, respectively.According to the average Tables 7 and 8 provide the average STD statistics of time link differences over a week under an elevation cut-off angle of 20 • and 30 • , respectively.According to the average results over the four time links, the time transfer accuracy is 0.234, 0.849, 0.254, and 0.571 ns for the GPS-only case, GLONASS-only case, Galileo-only case, and BDS-only case when the elevation cut-off angle is set to 20 • .The corresponding accuracy of the five-system combination case is 0.227 ns.In the case of the elevation cut-off angle of 30 • , the time transfer accuracy is 0.368 ns for the GPS-only case, 1.029 ns for the GLONASS-only case, 0.318 ns for the Galileo-only case, and 0.725 ns for the BDS-only case.As for the fivesystem integrated PPP time transfer, the corresponding accuracy is 0.264 ns.With the increment of the elevation cut-off angle, the PPP time transfer accuracy of each single-GNSS case significantly degrades.Although the PPP time transfer accuracy of the five-system combination case also declines with the increasing elevation masks, it is much more reliable than that of the single-GNSS cases.When the elevation mask is further increased to 40 • , a time transfer accuracy of 0.376 ns can still be achieved for the five-system integrated PPP.Tables 7 and 8 provide the average STD statistics of time link differences over a week under an elevation cut-off angle of 20° and 30°, respectively.According to the average results over the four time links, the time transfer accuracy is 0.234, 0.849, 0.254, and 0.571 ns for the GPS-only case, GLONASS-only case, Galileo-only case, and BDS-only case when the elevation cut-off angle is set to 20°.The corresponding accuracy of the five-system combination case is 0.227 ns.In the case of the elevation cut-off angle of 30°, the time transfer accuracy is 0.368 ns for the GPS-only case, 1.029 ns for the GLONASS-only case, 0.318 ns for the Galileo-only case, and 0.725 ns for the BDS-only case.As for the five-system integrated PPP time transfer, the corresponding accuracy is 0.264 ns.With the increment of the elevation cut-off angle, the PPP time transfer accuracy of each single-GNSS case significantly degrades.Although the PPP time transfer accuracy of the five-system combination case also declines with the increasing elevation masks, it is much more reliable than that of the single-GNSS cases.When the elevation mask is further increased to  To analyze the frequency stability of PPP time transfer under an elevation cut-off angle of 20 • , the MDEV statistics of three time links, BRUX-MGUE, BRUX-PTBB, and BRUX-SPT0, on DOY 252, 2021, are displayed in Figure 15.The MDEV value of the GLONASS-only case has a significant increase after the cut-off elevation is set to 20 • , and thus the GLONASSonly PPP cannot provide the reliable time transfer service for the three time links.The other three single-GNSS PPP time transfer cases also show an increase in the MDEV values compared with those at an elevation cut-off angle of 10 • (see Figure 6).For the five-system combined PPP time transfer, the increased cut-off elevation of 20 • has less effect on the frequency stability, and its time transfer results are considered to be more reliable than the single-GNSS cases.When the elevation cut-off angle is further increased to 30 • , the service rate of single-GNSS cases is reduced to a low level (see Table 5), and the time transfer results of single-GNSS cases are incomplete.Thus, they cannot provide a reliable time transfer service.Based on the MDEV statistics derived from the incomplete time links, it is difficult to analyze the frequency stability of single-GNSS PPP time transfer.By contrast, the service rate of the five-system combined PPP time transfer can still maintain 100.0%, and the MDEV values can still keep within a reasonable range.In view of the superior performance of the five-system integrated PPP time transfe under high cut-off elevations, more related results under the elevation cut-off angle of 20° 30°, and 40° are provided.Figure 16 shows the epoch-wise average number of visible sat ellites and TDOP value over all the employed days and stations under the elevation cut off angle of 20°, 30°.and 40°.The visible satellite numbers at the common epoch of al stations (13 stations) from all the days (7 days) are used to compute the epoch-wise aver age number of visible satellites, and it is also the case for the epoch-wise average TDOP value.According to the average results of epoch-wise statistics over 24 h, the visible sat ellite number is 22.3 for the cut-off angle 20°, 18.7 for the cut-off angle 30°, and 11.6 for th cut-off angle 40°.Regarding the TDOP, the numerical value is 0.82 for the cut-off angl 20°, 2.44 for the cut-off angle 30°, and 3.23 for the cut-off angle 40°.It is clear that both th visible satellite number and TDOP value get worse with the increasing elevation mask angles.Despite this, the five-system combination case can still guarantee the sufficien number of visible satellites and acceptable TDOP values even under a cut-off elevation angle of 40°. Figure 17 shows the MDEV value for three time links (BRUX-MGUE, BRUX-PTBB, and BRUX-SPT0) derived from the five-system integrated PPP when the elevation cut-off angle is set to 20  Figure 18 illustrates the average MDEV values over a week and all 12 time links at an average time of 10 2 , 10 3 , and 10 4 s under environments with limited sky view (at an elevation cut-off angle of 20 • , 30 • , and 40 • ).Based on the average MDEV values, the frequency stability of five-system integrated PPP time transfer can be 1.17 × 10 −13 , 9.64 × 10 −13 , and 8.42 × 10 −12 under the elevation cut-off angle of 20 • , 30 • , and 40 • when the average time is 10 2 s; the corresponding frequency stability can be 2.19 × 10 −14 , 1.29 × 10 −13 , and 4.80 × 10 −13 when the average time is 10 3 s; the corresponding frequency stability can be 8.85 × 10 −15 , 2.40 × 10 −14 , and 3.46 × 10 −14 when the average time is 10 4 s.According to these results, it can be concluded that the five-system combination can still provide a stable PPP time transfer service even in the case of observation conditions with limited sky view.

Enhanced Performance with a Priori ISB Constraints Considering Satellite Clock Datums
As both BDS-2 and QZSS constellations mainly provide the services over Asia-Pacific areas, only the ISBs of the GLONASS, Galileo, and BDS-3 against GPS are employed to investigate the performance enhancement owing to a priori ISB constraints for time transfer.The elevation mask angle is set to 10 • .Figure 19 illustrates the time series of original ISB estimates at station BRUX on DOY 248 to 254, 2021.It is seen that the day-boundary jumps of original ISB estimates are dramatic.For comparison, the adjusted ISBs with the estimated differences of satellite clock datums between two adjacent days are also provided in Figure 19.After considering the effects of satellite clock datums, the consistency of ISBs is significantly improved.This creates a precondition for the multi-system integrated PPP time transfer with a priori ISB constraints.
Remote Sens. 2023, 15, x FOR PEER REVIEW 22 of 27 Figure 18 illustrates the average MDEV values over a week and all 12 time links at an average time of 10 2 , 10 3 , and 10 4 s under environments with limited sky view (at an elevation cut-off angle of 20°, 30°, and 40°).Based on the average MDEV values, the frequency stability of five-system integrated PPP time transfer can be 1.17 × 10 −13 , 9.64 × 10 −13 , and 8.42 × 10 −12 under the elevation cut-off angle of 20°, 30°, and 40° when the average time is 10 2 s; the corresponding frequency stability can be 2.19 × 10 −14 , 1.29 × 10 −13 , and 4.80 × 10 −13 when the average time is 10 3 s; the corresponding frequency stability can be 8.85 × 10 −15 , 2.40 × 10 −14 , and 3.46 × 10 −14 when the average time is 10 4 s.According to these results, it can be concluded that the five-system combination can still provide a stable PPP time transfer service even in the case of observation conditions with limited sky view.

Enhanced Performance with a Priori ISB Constraints Considering Satellite Clock Datums
As both BDS-2 and QZSS constellations mainly provide the services over Asia-Pacific areas, only the ISBs of the GLONASS, Galileo, and BDS-3 against GPS are employed to investigate the performance enhancement owing to a priori ISB constraints for time transfer.The elevation mask angle is set to 10°. Figure 19 illustrates the time series of original ISB estimates at station BRUX on DOY 248 to 254, 2021.It is seen that the day-boundary jumps of original ISB estimates are dramatic.For comparison, the adjusted ISBs with the estimated differences of satellite clock datums between two adjacent days are also provided in Figure 19.After considering the effects of satellite clock datums, the consistency of ISBs is significantly improved.This creates a precondition for the multi-system integrated PPP time transfer with a priori ISB constraints.The a priori ISB estimates of the previous day are introduced into the data proces of the five-system integrated PPP time transfer of the next day.The time transfer res for three time links, BRUX-PTBB, BRUX-SPT0, and BRUX-USN7, on DOY 249 to 254, 2 are obtained.Table 9 lists the average STD statistics of daily time link differences betw PPP estimates and IGS final products over six days for the five-system combination with and without ISB constraints for three time links.Table 10 presents the average MD values over six days for three time links derived from the five-system integrated PPP and without ISB constraints at an average time of 10 2 , 10 3 , and 10 4 s.The time tran accuracy can be improved by 51.6%, 37.9%, and 45.3% for the three time links, res tively, after considering ISB constraints.As for the frequency stability, it can be impro The a priori ISB estimates of the previous day are introduced into the data processing of the five-system integrated PPP time transfer of the next day.The time transfer results for three time links, BRUX-PTBB, BRUX-SPT0, and BRUX-USN7, on DOY 249 to 254, 2021, are obtained.Table 9 lists the average STD statistics of daily time link differences between PPP estimates and IGS final products over six days for the five-system combination case with and without ISB constraints for three time links.Table 10 presents the average MDEV values over six days for three time links derived from the five-system integrated PPP with and without ISB constraints at an average time of 10 2 , 10 3 , and 10 4 s.The time transfer accuracy can be improved by 51.6%, 37.9%, and 45.3% for the three time links, respectively, after considering ISB constraints.As for the frequency stability, it can be improved by 14.8%, 21.6%, and 15.0% at an average time of 10 2 s, 6.1%, 5.3%, and 7.6% at an average time of 10 3 s, and 29.6%, 20.0%, and 29.5% at an average time of 10 4 s, respectively, for the three time links for the five-system integrated PPP with ISB constraints over that without ISB constraints.

Discussion
PPP is one of the satellite-based positioning technologies.Due to the strong satellite geometry and the increased number of available satellites, both the accuracy and frequency stability of PPP time transfer can be improved by including more satellite systems, especially for the frequency stability.Since the users also have the requirement of high-accuracy time transfer at locations with limited sky view in the practical applications, the considerable service rates should be ensured.The service rates of single-GNSS PPP time transfer cases with high cut-off elevation angles are low, while those of the five-system integrated PPP time transfer remain at 100.0%.The satisfactory time transfer performance of the fivesystem integrated PPP can facilitate the applications of PPP-based time transfer.
There is a strong correlation among the estimated parameters in the PPP.The dayboundary jumps of ISB estimates can be mitigated by considering the difference in the satellite clock datums between two adjacent days.The adjusted ISBs can be taken as the a priori constraints to reduce the correlation among the estimated parameters, so that the time transfer performance can be enhanced.
With the rise of large-scale low earth orbit (LEO) constellations, the LEO-based performance enhancement for GNSS applications has become a research focus.The PPP time transfer technology can also benefit from the integrated data processing with the satellites at high, medium, and low orbits.On the one hand, the inclusion of tens of LEO satellites can further increase the number of available satellites and improve the satellite geometry, and thus the accuracy, frequency stability, and service rates of PPP time transfer can be enhanced.On the other hand, the separation of various estimated parameters will be easier, as the LEO satellites travel faster over stations and show rapid changes in spatial

7 Figure 1 .
Figure 1.Geographical distribution of 13 MGEX stations employed for time transfer.The red refers to the center node of all the time links.

Figure 1 .
Figure 1.Geographical distribution of 13 MGEX stations employed for time transfer.The red point refers to the center node of all the time links.

Figure 2 .
Figure 2. Average number of visible satellites of the five-system combination on DOY 254,

Figure 2 .
Figure 2. Average number of visible satellites of the five-system combination on DOY 254, 2021.

Figure 4 .
Figure 4. Difference in the receiver clock offset between PPP estimates and IGS final products for the five-system combination case and the four GNSS single-system cases at stations BRUX and SPT0 on DOY 252, 2021.

Figure 5 .
Figure 5. Difference in the time link BRUX-SPT0 between PPP estimates and IGS final products for the five-system combination case and the four GNSS single-system cases on DOY 252, 2021.

Figure 4 .
Figure 4. Difference in the receiver clock offset between PPP estimates and IGS final products for the five-system combination case and the four GNSS single-system cases at stations BRUX and SPT0 on DOY 252, 2021.

Figure 4 .
Figure 4. Difference in the receiver clock offset between PPP estimates and IGS final products for the five-system combination case and the four GNSS single-system cases at stations BRUX and SPT0 on DOY 252, 2021.

Figure 5 .
Figure 5. Difference in the time link BRUX-SPT0 between PPP estimates and IGS final products for the five-system combination case and the four GNSS single-system cases on DOY 252, 2021.

Figure 5 .
Figure 5. Difference in the time link BRUX-SPT0 between PPP estimates and IGS final products for the five-system combination case and the four GNSS single-system cases on DOY 252, 2021.

Figure 6 .
Figure 6.MDEV values of three time links under a cut-off elevation angle of 10° on DOY 252, 2021.Figure 6. MDEV values of three time links under a cut-off elevation angle of 10 • on DOY 252, 2021.

Figure 6 .
Figure 6.MDEV values of three time links under a cut-off elevation angle of 10° on DOY 252, 2021.Figure 6. MDEV values of three time links under a cut-off elevation angle of 10 • on DOY 252, 2021.

Figure 10 .
Figure 10.Average MDEV values over a week for 12 time links at an average time of 10 2 , 10 3 , and 10 4 s.The MDEV results for the time links with CEBR, CRO1, KOKB, MATE, MGUE and TID1 stations are provided in (a), and those for the time links with PTBB, ROAG, SPT0, USN7, WTZR and ZIM3 stations are provided in (b).

Figure 10 .
Figure 10.Average MDEV values over a week for 12 time links at an average time of 10 2 , 10 3 , and 10 4 s.The MDEV results for the time links with CEBR, CRO1, KOKB, MATE, MGUE and TID1 stations are provided in (a), and those for the time links with PTBB, ROAG, SPT0, USN7, WTZR and ZIM3 stations are provided in (b).

Figure 11 .
Figure 11.Average MDEV values over a week and all 12 time links at an average time of 10 2 , 10 3 , and 10 4 s.

Figure 11 .
Figure 11.Average MDEV values over a week and all 12 time links at an average time of 10 2 , 10 3 , and 10 4 s.

Figure 12 .
Figure 12.Service rates of single-GNSS PPP time transfer cases for each time link using the cut-off elevation angle of 20°, 30°, and 40°.

igure 12 .
Service rates of single-GNSS PPP time transfer cases for each time link using the cut-off elevation angle of 20 • , 30 • , and 40 • .

Figure 13 .
Figure 13.Difference of receiver clock offset between PPP estimates and IGS final products for five PPP time transfer cases at stations BRUX and SPT0 under an elevation cut-off angle of 20° on DOY 252, 2021.

Figure 14 .
Figure 14.Difference of the time link BRUX-SPT0 between PPP estimates and IGS final products for five PPP time transfer cases under an elevation cut-off angle of 20° on DOY 252, 2021.

Figure 13 .
Figure 13.Difference of receiver clock offset between PPP estimates and IGS final products for five PPP time transfer cases at stations BRUX and SPT0 under an elevation cut-off angle of 20 • on DOY 252, 2021.

Figure 13 .
Figure 13.Difference of receiver clock offset between PPP estimates and IGS final products for five PPP time transfer cases at stations BRUX and SPT0 under an elevation cut-off angle of 20° on DOY 252, 2021.

Figure 14 .
Figure 14.Difference of the time link BRUX-SPT0 between PPP estimates and IGS final products for five PPP time transfer cases under an elevation cut-off angle of 20° on DOY 252, 2021.

Figure 14 .
Figure 14.Difference of the time link BRUX-SPT0 between PPP estimates and IGS final products for five PPP time transfer cases under an elevation cut-off angle of 20 • on DOY 252, 2021.

Figure 15 .
Figure 15.MDEV values of three time links under an elevation cut-off angle of 20° on DOY 252 2021.

Figure 15 .
Figure 15.MDEV values of three time links under an elevation cut-off angle of 20 • on DOY 252, 2021.In view of the superior performance of the five-system integrated PPP time transfer under high cut-off elevations, more related results under the elevation cut-off angle of 20 • , 30 • , and 40 • are provided.Figure 16 shows the epoch-wise average number of visible satellites and TDOP value over all the employed days and stations under the elevation cut-off angle of 20 • , 30 • .and 40 • .The visible satellite numbers at the common epoch of all stations (13 stations) from all the days (7 days) are used to compute the epoch-wise average number of visible satellites, and it is also the case for the epoch-wise average TDOP value.According to the average results of epoch-wise statistics over 24 h, the visible satellite number is 22.3 for the cut-off angle 20 • , 18.7 for the cut-off angle 30 • , and 11.6 for the cut-off angle 40 • .Regarding the TDOP, the numerical value is 0.82 for the cut-off angle 20 • , 2.44 for the cut-off angle 30 • , and 3.23 for the cut-off angle 40 • .It is clear that both the visible satellite number and TDOP value get worse with the increasing elevation mask angles.Despite this, the five-system combination case can still guarantee the sufficient number of visible satellites and acceptable TDOP values even under a cut-off elevation angle of 40 • .Figure17shows the MDEV value for three time links (BRUX-MGUE, BRUX-PTBB, and BRUX-SPT0) derived from the five-system integrated PPP when the elevation cut-off angle is set to 20 • , 30 • , and 40 • on DOY 252, 2021.It is found that the MDEV value increases with the increment of the cut-off elevation.The corresponding MDEV value of the three time links under a cut-off angle of 40 • can reach 8.30 × 10 −15 , 1.41 × 10 −14 , and 9.86 × 10 −15 at an average time of 10 4 s.Thus, the MDEV value of the three time links at an elevation mask of 40 • still remains within a reasonable range, and the corresponding PPP time transfer results can still have satisfactory frequency stability.Figure18illustrates the average MDEV values over a week and all 12 time links at an average time of 10 2 , 10 3 , and 10 4 s under environments with limited sky view (at an elevation cut-off angle of 20 • , 30 • , and 40 • ).Based on the average MDEV values, the frequency stability of five-system integrated PPP time transfer can be 1.17 × 10 −13 , 9.64 × 10 −13 , and 8.42 × 10 −12 under the elevation cut-off angle of 20 • , 30 • , and 40 • when the average time is 10 2 s; the corresponding frequency stability can be 2.19 × 10 −14 , , 30 • , and 40 • on DOY 252, 2021.It is found that the MDEV value increases with the increment of the cut-off elevation.The corresponding MDEV value of the three time links under a cut-off angle of 40 • can reach 8.30 × 10 −15 , 1.41 × 10 −14 , and 9.86 × 10 −15 at an average time of 10 4 s.Thus, the MDEV value of the three time links at an elevation mask of 40 • still remains within a reasonable range, and the corresponding PPP time transfer results can still have satisfactory frequency stability.

27 Figure 16 .
Figure 16.Epoch-wise average number of visible satellites and TDOP value over all the employed days and stations under the elevation cut-off angle of 20°, 30°, and 40°.

Figure 17
Figure 17 shows the MDEV value for three time links (BRUX-MGUE, BRUX-PTBB, and BRUX-SPT0) derived from the five-system integrated PPP when the elevation cut-off angle is set to 20°, 30°, and 40° on DOY 252, 2021.It is found that the MDEV value increases with the increment of the cut-off elevation.The corresponding MDEV value of the three time links under a cut-off angle of 40° can reach 8.30 × 10 −15 , 1.41 × 10 −14 , and 9.86 × 10 −15 at an average time of 10 4 s.Thus, the MDEV value of the three time links at an elevation mask of 40° still remains within a reasonable range, and the corresponding PPP time transfer results can still have satisfactory frequency stability.

Figure 17 .
Figure 17.MDEV value of three time links derived from five-system integrated PPP under the elevation cut-off angle of 20°, 30°, and 40° on DOY 252, 2021.

Figure 16 . 27 Figure 16 .
Figure 16.Epoch-wise average number of visible satellites and TDOP value over all the employed days and stations under the elevation cut-off angle of 20 • , 30 • , and 40 • .

Figure 17
Figure 17 shows the MDEV value for three time links (BRUX-MGUE, BRUX-PTBB, and BRUX-SPT0) derived from the five-system integrated PPP when the elevation cut-off angle is set to 20°, 30°, and 40° on DOY 252, 2021.It is found that the MDEV value increases with the increment of the cut-off elevation.The corresponding MDEV value of the three time links under a cut-off angle of 40° can reach 8.30 × 10 −15 , 1.41 × 10 −14 , and 9.86 × 10 −15 at an average time of 10 4 s.Thus, the MDEV value of the three time links at an elevation mask of 40° still remains within a reasonable range, and the corresponding PPP time transfer results can still have satisfactory frequency stability.

Figure 17 .
Figure 17.MDEV value of three time links derived from five-system integrated PPP under the elevation cut-off angle of 20°, 30°, and 40° on DOY 252, 2021.

Figure 17 .
Figure 17.MDEV value of three time links derived from five-system integrated PPP under the elevation cut-off angle of 20 • , 30 • , and 40 • on DOY 252, 2021.

Figure 18 .
Figure 18.Average MDEV values over a week and all the 12 time links at an average time of 10 2 , 10 3 , and 10 4 s under environments with limited sky view.

Figure 18 . 27 Figure 19 .
Figure 18.Average MDEV values over a week and all the 12 time links at an average time of 10 2 , 10 3 , and 10 4 s under environments with limited sky view.Remote Sens. 2023, 15, x FOR PEER REVIEW 23 of 27

Figure 20
Figure 20 illustrates the difference of the time link BRUX-SPT0 between PPP estimates and IGS final products for the five-system combination case with and without ISB constraints on DOY 249, 2021.The MDEV values of the time link BRUX-SPT0 derived from five-system integrated PPP with and without ISB constraints on DOY 249, 2021, are provided in Figure 21.After taking the a priori ISB constraints into account, both the time

Figure 19 .
Figure 19.Time series of original ISB estimates and those considering satellite clock datums at station BRUX on DOY 248 to 254, 2021.

Figure 20
Figure 20 illustrates the difference of the time link BRUX-SPT0 between PPP estimates and IGS final products for the five-system combination case with and without ISB constraints on DOY 249, 2021.The MDEV values of the time link BRUX-SPT0 derived from five-system integrated PPP with and without ISB constraints on DOY 249, 2021, are provided in Figure 21.After taking the a priori ISB constraints into account, both the time transfer accuracy and frequency stability can be significantly improved.

Figure 19 .
Figure 19.Time series of original ISB estimates and those considering satellite clock datums at tion BRUX on DOY 248 to 254, 2021.

Figure 20
Figure 20  illustrates the difference of the time link BRUX-SPT0 between PPP e mates and IGS final products for the five-system combination case with and without constraints on DOY 249, 2021.The MDEV values of the time link BRUX-SPT0 deri from five-system integrated PPP with and without ISB constraints on DOY 249, 2021, provided in Figure21.After taking the a priori ISB constraints into account, both the t transfer accuracy and frequency stability can be significantly improved.

Figure 20 .
Figure 20.Difference of the time link BRUX-SPT0 between PPP estimates and IGS final products the five-system combination case with and without ISB constraints on DOY 249, 2021.

Figure 20 .
Figure 20.Difference of the time link BRUX-SPT0 between PPP estimates and IGS final products for the five-system combination case with and without ISB constraints on DOY 249, 2021.Remote Sens. 2023, 15, x FOR PEER REVIEW 24

Figure 21 .
Figure 21.MDEV values of the time link BRUX-SPT0 derived from the five-system integrated with and without ISB constraints on DOY 249, 2021.

Figure 21 .
Figure 21.MDEV values of the time link BRUX-SPT0 derived from the five-system integrated PPP with and without ISB constraints on DOY 249, 2021.

Table 1 .
Detailed information of 13 MGEX stations.

Table 1 .
Detailed information of 13 MGEX stations.

Table 2 .
Processing strategy of five-system integrated PPP time transfer.

Table 4 .
Average value of day-to-day variation ratios of MDEV statistics over seven days (DOY 248 to 254, 2021) and three specific average times (10 2 , 10 3 , and 10 4 s) for three time links.
• , and show good performance in terms of data integrity.

Table 7 .
Average STD statistics of time link differences over a week under an elevation cut-off angle of 20 • (unit: ns).

Table 8 .
Average STD statistics of time link differences over a week under an elevation cut-off angle of 30 • (unit: ns).

Table 9 .
Average STD statistics of daily time link differences between PPP estimates and IGS final products over six days for the five-system combination case with and without ISB constraints for three time links (unit: ns).