Improved Ultra-Rapid UT1-UTC Determination and Its Preliminary Impact on GNSS Satellite Ultra-Rapid Orbit Determination

Errors in ultra-rapid UT1-UTC primarily affect the overall rotation of spatial datum expressed by GNSS (Global Navigation Satellite System) satellite ultra-rapid orbit. In terms of existing errors of traditional strategy, e.g., piecewise linear functions, for ultra-rapid UT1-UTC determination, and the requirement to improve the accuracy and consistency of ultra-rapid UT1-UTC, the potential to improve the performance of ultra-rapid UT1-UTC determination based on an LS (Least Square) + AR (Autoregressive) combination model is explored. In this contribution, based on the LS+AR combination model and by making joint post-processing/rapid UT1-UTC observation data, we propose a new strategy for ultra-rapid UT1-UTC determination. The performance of the new strategy is subsequently evaluated using data provided by IGS (International GNSS Services), iGMAS (international GNSS Monitoring and Assessment System), and IERS (International Earth Rotation and Reference Systems Service). Compared to the traditional strategy, the numerical results over more than 1 month show that the new strategy improved ultra-rapid UT1-UTC determination by 29–43%. The new strategy can provide a reference for GNSS data processing to improve the performance of ultra-rapid products.


Introduction
In recent years, the requirement of GNSS (Global Navigation Satellite System) to have real-time precise positioning is becoming more and more urgent for real-time detection research in the field of Earth sciences and for the increasing number of real-time high-precision positioning applications [1][2][3]. GNSS satellite ultra-rapid orbit is the key in determining real-time precise orbits and achieving real-time precision point positioning (PPP) [4,5]. The accuracy and stability of GNSS satellite ultra-rapid orbit determination directly or indirectly affects and limits the performance and promotion of real-time PPP and related applications [6][7][8]. UT1 is the angle of the Earth's rotation about the CIP (Celestial Intermediate Pole) axis defined by its conventional linear relation to the Earth Rotation Angle (ERA), and UT1-UTC is the difference between the UT1 parameter derived from observations and UTC (Coordinated Universal Time) [9]. Ultra-rapid UT1-UTC is currently one of the main factors affecting the performance of GNSS satellite ultra-rapid orbit determination [10]. Compared with the present performance, the existing error in ultra-rapid UT1-UTC has a significant impact on ultra-rapid orbit determination [6,[11][12][13]. Therefore, further improvement in the processing strategy of ultra-rapid

Basic Principles of LS + AR Combination Model
The fitted equation of an LS model for periodic and linear terms of UT1R-TAI can be written as follows [29]: [c i cos(2πt/PE i ) + d i sin(2πt/PE i )] (1) Remote Sens. 2020, 12,3584 3 of 15 where f t is the time sequence of UT1R-TAI; k refers to the number of periodic terms in these series; PE i is the corresponding oscillation period obtained by spectrum analysis; and a, b, c i , and d i are the unknown parameters to be estimated. The equation of the AR(p) model for the fitted residual terms of UT1R-TAI can be written as follows [29,30]: where ϕ i is the unknown parameter to be estimated, Z t is the fitted residual terms, a t refers to white noise, p is the order of the AR model, and the search range of p is 1 ∼ √ N; the rule for estimating the AR model order is to take the minimum values of FPE (Final Prediction Error) [31], AIC (Akaike Information Criterion) [32], and BIC (Bayes Information Criterion) [33], where N represents the number of residual terms used for modeling.
In evaluating the quality of earth rotation parameters, the mean absolute error (MAS) [27] adopted by EOP PCC (prediction comparison campaign) can be written as follows: where l and N are the prediction span and the number of prediction periods, respectively. O l i and P l i refer to the reference value and predicted value of earth rotation parameters. The value of O l i − P l i is the absolute residual. Figure 1 shows the new strategy for ultra-rapid UT1-UTC determination based on the LS + AR combination model and the traditional strategy. The other details for new strategy are displayed in the following.

New Strategy of Ultra-Rapid UT1-UTC Determination for Ultra-Rapid Orbit
Ultra-rapid UT1-UTC includes observation and prediction parts, and a description of the UT1-UTC observation data used by our new strategy for different periods is summarized in Figure 2. Among them, the green line and red line represent the observation and prediction parts of the ultra-rapid ERP, respectively. In fact, UT1-UTC can only be determined from the VLBI (Very Long Baseline Interferometry) technique and derived from Lunar Laser Ranging (LLR) in a limited scope [34]. Different from other parameters in ultra-rapid ERP, Figure 2 also shows that UT1-UTC cannot be solved directly by GNSS [35] and that the observation part in ultra-rapid UT1-UTC should also be predicted. In Figure 2, 00, 06, 12, and 18 represent time periods of four ultra-rapid products per day. In addition, Figures 1 and 2 show the general idea of the new strategy for ultra-rapid UT1-UTC determination, which includes joint UT1-UTC data in Figure 1 that comes from EOP 14 C04 and IERS gpsrapid.daily.
Based on the joint UT1-UTC data sequence and time delay of the observed UT1-UTC data, in order to determine ultra-rapid UT1-UTC, at least 2 days of UT1-UTC need to be predicted ( Figure 2); thus, two processing methods are designed to determinate ultra-rapid UT1-UTC: Method 1 (A)-The UT1-UTC data sequence is predicted for 2 days, and then the second-order Lagrange interpolation algorithm is used to obtain the observation and prediction parts of the ultra-rapid UT1-UTC parameter, which is updated once a day.

Method 2 (B)-
The UT1-UTC data sequence is predicted for one day, and then the predicted value of this step is added to the original sequence again for a 1-day prediction. Finally, the second-order Lagrange interpolation algorithm is used to obtain the observation and prediction parts of the ultra-rapid UT1-UTC parameter, which is updated once a day. Lagrange interpolation algorithm is used to obtain the observation and prediction parts of the ultrarapid UT1-UTC parameter, which is updated once a day.   Lagrange interpolation algorithm is used to obtain the observation and prediction parts of the ultrarapid UT1-UTC parameter, which is updated once a day.   Method 2 (B) is smaller, but the final observed UT1-UTC data in the second iteration are simulated, whereas the UT1-UTC data in Method 1 (A) are all observed. Therefore, the performance of the two methods should be evaluated by experimentation and analysis, and the best processing method will be considered in the new strategy.

Results
The quality of ultra-rapid UT1-UTC based on different cases was firstly analyzed in detail. In this experiment, data from EOP 14 C04, IERS gpsrapid.daily, and IGS (International GNSS Services) and iGMAS (international GNSS Monitoring and Assessment System) ultra-rapid products were used. Here, MJD (Modified Julian Date) 58356-58400 was the test period, and when calculating MAS, the reference values were from EOP 14 C04 and determined (predicted) values, including observation and prediction parts from different cases. The initial MJDs were 54,466 (EOP 14 C04), 58,326 (IERS gpsrapid.daily), and 58,326 (ultra-rapid products).
To analyze and compare the performance of new and traditional strategies, four cases [Case 1 (A), Case 2 (B), Case 3 (iGMAS), and Case 4 (IGU)] were conducted. In this part, ultra-rapid UT1-UTC for different cases was first obtained, and then the absolute residuals and MAS of different cases were subsequently calculated by comparing these ultra-rapid UT1-UTCs with EOP 14 C04.
Cases 1 (A) and 2 (B) are designed to analyze the performance of Methods 1 (A) and 2 (B) in the new strategy. Cases 3 (iGMAS) and 4 (IGU) are designed to analyze performance related to the traditional strategy. Here, Case 1 (A) is where the ultra-rapid UT1-UTC parameter is generated by Method 1 (A), and Case 2 (B) is where the ultra-rapid UT1-UTC parameter is generated by Method 2 (B). Both Cases 1 (A) and 2 (B) generate ultra-rapid UT1-UTC based on more than 10-year joint products and the LS+AR combination model. Additionally, Case 3 (iGMAS) is where the ultra-rapid UT1-UTC parameter is provided by iGMAS ultra-rapid products, and Case 4 (IGU) is where the ultra-rapid UT1-UTC parameter is provided by IGS ultra-rapid products. Both iGMAS ultra-rapid and IGS ultra-rapid products generate ultra-rapid UT1-UTC based on piecewise linear functions and consecutive rapid UT1-UTC observation data for several days. Figure 3 shows the results determined for ultra-rapid UT1-UTC for different cases, in which the top panel of Figure 3 shows the absolute residuals of the observed part under different processed cases, and the bottom panel of Figure 3 is the absolute residuals of the predicted part. The legends in Figure 3 are the mean absolute errors of UT1-UTC parameters for observed and predicted parts in the ultra-rapid product. The residual distributions for Cases 2 (B) and 4 (IGU) are also displayed in Appendix A. Table 1 presents the summary of Figure 3 as well as the relative gain (RG) of ultra-rapid UT1-UTC for case 2 compared to IGS and iGMAS ultra-rapid products.  Compared to other cases, either observed UT1-UTC or predicted UT1-UTC, Figure 3 and Appendix A show that case 2 (B) had the best performance. This indicates that the prediction span seriously impacts the prediction performance of ultra-rapid UT1-UTC. After shortening the prediction span, the performance of iterative predictions was better. In addition, from Table 1, compared to Case 4 (IGU), we can see that case 2 (B) improved by 42.70% and 29.17% for the observed and predicted parts, respectively. At the same time, compared to Case 3 (iGMAS), case 2 (B) improved by 42.05% and 37.61%. From the above results, it was concluded that the new strategy for ultra-rapid UT1-UTC determination should follow Method 2 (B).

Discussion
The above results show that using the joint products and LS+AR combination model improved the accuracy and consistency of ultra-rapid UT1-UTC determination and IERS EOP 14 C04. In order to analyze and discuss the direct and extended applications of this new strategy, the preliminary impact on ultra-rapid orbit determination is firstly assessed. In addition, the preliminary prediction and quality assessment of ultra-rapid LOD is performed.

Discussion of Ultra-Rapid Orbit Determination
This section only assesses the coordinate transformation impact of ultra-rapid UT1-UTC on orbits coordinate of prediction part in the ultra-rapid product. The experimental steps designed to analyze the impact of using ultra-rapid UT1-UTC in transforming coordinates from the celestial reference frame (defined by the mean pole and equinox at J2000.0) to the terrestrial frame (IGS14) are as follows: (1) Obtain the simulated true satellite orbit state-the initial state of satellite orbit and the related state transition matrix are fitted by the final precise orbit and ERP products. (2) Obtain the simulated true satellite orbit coordinates in J2000.0-the coordinates of the satellite in ultra-rapid orbit in the celestial reference frame are calculated by integrating the satellite state parameters and related state transition matrix obtained in the first step and the final precise ERP products. Compared to other cases, either observed UT1-UTC or predicted UT1-UTC, Figure 3 and Appendix A show that case 2 (B) had the best performance. This indicates that the prediction span seriously impacts the prediction performance of ultra-rapid UT1-UTC. After shortening the prediction span, the performance of iterative predictions was better. In addition, from Table 1, compared to Case 4 (IGU), we can see that case 2 (B) improved by 42.70% and 29.17% for the observed and predicted parts, respectively. At the same time, compared to Case 3 (iGMAS), case 2 (B) improved by 42.05% and 37.61%. From the above results, it was concluded that the new strategy for ultra-rapid UT1-UTC determination should follow Method 2 (B).

Discussion
The above results show that using the joint products and LS+AR combination model improved the accuracy and consistency of ultra-rapid UT1-UTC determination and IERS EOP 14 C04. In order to analyze and discuss the direct and extended applications of this new strategy, the preliminary impact on ultra-rapid orbit determination is firstly assessed. In addition, the preliminary prediction and quality assessment of ultra-rapid LOD is performed.

Discussion of Ultra-Rapid Orbit Determination
This section only assesses the coordinate transformation impact of ultra-rapid UT1-UTC on orbits coordinate of prediction part in the ultra-rapid product. The experimental steps designed to analyze the impact of using ultra-rapid UT1-UTC in transforming coordinates from the celestial reference frame (defined by the mean pole and equinox at J2000.0) to the terrestrial frame (IGS14) are as follows: (1) Obtain the simulated true satellite orbit state-the initial state of satellite orbit and the related state transition matrix are fitted by the final precise orbit and ERP products. (2) Obtain the simulated true satellite orbit coordinates in J2000.0-the coordinates of the satellite in ultra-rapid orbit in the celestial reference frame are calculated by integrating the satellite state parameters and related state transition matrix obtained in the first step and the final precise ERP products. (3) Obtain the satellite orbit coordinates in IGS14 with different ultra-rapid ERPs [36]-the coordinates of the satellite in ultra-rapid orbit obtained in the second step and the ultra-rapid ERPs (IGG ERP and IGU ERP, where IGG ERP represents ultra-rapid UT1-UTC and is generated by the above new strategy; other parameters in ultra-rapid ERP come from the IGS ultra-rapid product, and IGU ERP represents ultra-rapid ERP and is provided by the IGS ultra-rapid product) obtained in the above part are used to transform the coordinates of the satellite in ultra-rapid orbit in the terrestrial frame. (4) Comparative analysis-the 6-hour prediction for each arc of satellite in ultra-rapid orbit obtained in the third step is compared with the final precise orbit products, and the residual of each satellite and its maximum (Max), average (Ave), and RMS (Root Mean Square) are calculated.
The experimental period was MJD 58356-58400, and the final precise orbit and ERP products are provided by IGS. Figure 4 represents the coordinate residuals between the 6-hour prediction of ultra-rapid orbits and the IGS final precise orbits in the X, Y and Z directions when IGG ERP is used. Similarly, Figure 5 represents these coordinate residuals when IGU ERP is used. (3) Obtain the satellite orbit coordinates in IGS14 with different ultra-rapid ERPs [36]-the coordinates of the satellite in ultra-rapid orbit obtained in the second step and the ultra-rapid ERPs (IGG ERP and IGU ERP, where IGG ERP represents ultra-rapid UT1-UTC and is generated by the above new strategy; other parameters in ultra-rapid ERP come from the IGS ultra-rapid product, and IGU ERP represents ultra-rapid ERP and is provided by the IGS ultra-rapid product) obtained in the above part are used to transform the coordinates of the satellite in ultrarapid orbit in the terrestrial frame. (4) Comparative analysis-the 6-hour prediction for each arc of satellite in ultra-rapid orbit obtained in the third step is compared with the final precise orbit products, and the residual of each satellite and its maximum (Max), average (Ave), and RMS (Root Mean Square) are calculated.
The experimental period was MJD 58356-58400, and the final precise orbit and ERP products are provided by IGS. Figure 4 represents the coordinate residuals between the 6-hour prediction of ultrarapid orbits and the IGS final precise orbits in the X, Y and Z directions when IGG ERP is used. Similarly, Figure 5 represents these coordinate residuals when IGU ERP is used.   Figure 4 shows that the maximum residuals for all satellites are 2.0 and 1.9 cm in the X and Y directions when IGG ERP is used. Figure 5 shows the maximum residuals are 2.4 and 2.5 cm when IGU ERP is used. Therefore, compared to the maximum residuals for all satellites when IGU ERP is used, the maximum residuals on coordinate transformation in ultra-rapid orbit determination are improved by 17% and 24% in the X and Y directions, respectively, when IGG ERP is used. Both Figures 4 and 5 show that the maximum residuals for all satellites are 1.2 cm. The error in the Z direction is caused by polar motion; polar motion is not the subject of our discussion here, so only the X and Y directions are displayed when the maximum, average, and RMS of residuals of each satellite are calculated in next analysis.
The top panel in Figure 6 represents the residual statistics of the 6-hour prediction for each arc for all satellites in the X and Y directions, respectively. The bottom panel in Figure 6 represents the improved percentages of residual statistics for all satellites in the X and Y directions, respectively. Figure 6 shows that the residuals when IGG ERP is used are better than the residuals when IGU ERP is used, in terms of maximum and average residuals. However, the expected improvement in the RMS of residuals was not seen, and this phenomenon may be due to the better consistency of UT1-UTC parameters between the final and ultra-rapid ERP products provided by IGS. Therefore, the UT1-UTC residuals and statistics between the final ERP products provided by IGS and EOP 14 C04 provided by IERS are also calculated here, as shown in Figure 7. From Figure 7, it can be seen that the UT1-UTC  Figure 5. Coordinate residuals between the 6-hour predicted part of ultra-rapid orbits and the IGS final precise orbits in the X, Y, and Z directions when IGU ERP is used. MAX_X, MAX_Y, and MAX_Z are the maximum residuals for all satellites in X, Y, and Z directions, respectively. Figure 4 shows that the maximum residuals for all satellites are 2.0 and 1.9 cm in the X and Y directions when IGG ERP is used. Figure 5 shows the maximum residuals are 2.4 and 2.5 cm when IGU ERP is used. Therefore, compared to the maximum residuals for all satellites when IGU ERP is used, the maximum residuals on coordinate transformation in ultra-rapid orbit determination are improved by 17% and 24% in the X and Y directions, respectively, when IGG ERP is used. Both Figures 4 and 5 show that the maximum residuals for all satellites are 1.2 cm. The error in the Z direction is caused by polar motion; polar motion is not the subject of our discussion here, so only the X and Y directions are displayed when the maximum, average, and RMS of residuals of each satellite are calculated in next analysis.
The top panel in Figure 6 represents the residual statistics of the 6-hour prediction for each arc for all satellites in the X and Y directions, respectively. The bottom panel in Figure 6 represents the improved percentages of residual statistics for all satellites in the X and Y directions, respectively. Figure 6 shows that the residuals when IGG ERP is used are better than the residuals when IGU ERP is used, in terms of maximum and average residuals. However, the expected improvement in the RMS of residuals was not seen, and this phenomenon may be due to the better consistency of UT1-UTC parameters between the final and ultra-rapid ERP products provided by IGS. Therefore, the UT1-UTC residuals and statistics between the final ERP products provided by IGS and EOP 14 C04 provided by IERS are also calculated here, as shown in Figure 7. From Figure 7, it can be seen that the UT1-UTC parameter sequence provided by IGS in the final ERP products had a deviation of about 63.15 µs from the UT1-UTC parameter sequence published by IERS in EOP 14 C04. On the other hand, in terms of the absolute RMS and the accuracy of ultra-rapid orbits, the RMS of residuals can be considered equivalent when IGG ERP and IGU ERP are used.

Discussion on Ultra-Rapid LOD Determination
Ultra-rapid LOD includes observation and prediction parts, and a description of LOD observation data used by the new strategy for ultra-rapid LOD prediction in different periods is summarized in Figure 8.

Discussion on Ultra-Rapid LOD Determination
Ultra-rapid LOD includes observation and prediction parts, and a description of LOD observation data used by the new strategy for ultra-rapid LOD prediction in different periods is summarized in Figure 8. Different from ultra-rapid UT1-UTC determination, Figure 8 shows that the observation part in ultra-rapid LOD can be solved directly by GNSS, and only the prediction part should be predicted. Therefore, two processing methods are designed to obtain the prediction part in ultra-rapid LOD as follows: Method 1 (A)-LOD data sequences are predicted for one day, which is updated every 6 h.  Different from ultra-rapid UT1-UTC determination, Figure 8 shows that the observation part in ultra-rapid LOD can be solved directly by GNSS, and only the prediction part should be predicted. Therefore, two processing methods are designed to obtain the prediction part in ultra-rapid LOD as follows: Method 1 (A)-LOD data sequences are predicted for one day, which is updated every 6 h. The data from EOP 14 C04 [MJD: 54,466-58,400 and MJD: 58,356-38,400 are the test period as the reference (true) value compared to the determined value, also namely prediction part value, of different cases], and IGS and iGMAS ultra-rapid products (MJD: 58,326-58,400) are used, and four cases [Case 1 (A), Case 2 (B), Case 3 (iGMAS), and Case 4 (IGU)] are conducted similarly. Figure 9 shows the absolute residuals of the LOD parameter for predicted parts under different processed cases, and the legend in Figure 9 shows mean absolute errors of the LOD parameter for predicted parts in the ultra-rapid product. Table 2 presents a summary of Figure 9 and the relative gain of the LOD parameter for the predicted part in case 1 compared to IGS and iGMAS products. From Table 2, compared to Case 4 (IGU), we can see that the improvement in case 1 (A) was 20.83%. At the same time, compared to Case 3 (iGMAS), the improvement in case 1 (A) was 64.81%.
Remote Sens. 2020, 12, 3584 11 of 16 Figure 9. Results of the prediction part of ultra-rapid LOD. Therefore, based on the LS + AR combination model, the above assessment shows that the new strategy further improves the performance of ultra-rapid UT1-UTC determination and the prediction performance of ultra-rapid LOD when using joint observation data.
On the one hand, ultra-rapid UT1-UTC determined by the new strategy can directly apply to GNSS satellite ultra-rapid orbit determination and can directly constrain the known value. The prior value in ultra-rapid LOD can also be derived directly from ultra-rapid UT1-UTC. In addition, ultrarapid LOD prediction directly by the new strategy should depend on the observation part in ultrarapid LOD, so iterations are essential, and the time consumed needs to be balanced. Therefore, the impact of the new strategy for ultra-rapid LOD on ultra-rapid orbit determination needs to be further analyzed and researched. In addition, typically a combination of GNSS [37] and Satellite Laser Ranging (SLR) observations in the LAser GEOdynamics Satellite (LAGEOS) and Etalon satellites are used to determine LOD [38]. In a limited scope, Doppler Orbitography and Radio-positioning Integrated by Satellite (DORIS) observations can be employed to derive LOD values as well [39].  Therefore, based on the LS + AR combination model, the above assessment shows that the new strategy further improves the performance of ultra-rapid UT1-UTC determination and the prediction performance of ultra-rapid LOD when using joint observation data.
On the one hand, ultra-rapid UT1-UTC determined by the new strategy can directly apply to GNSS satellite ultra-rapid orbit determination and can directly constrain the known value. The prior value in ultra-rapid LOD can also be derived directly from ultra-rapid UT1-UTC. In addition, ultra-rapid LOD prediction directly by the new strategy should depend on the observation part in ultra-rapid LOD, so iterations are essential, and the time consumed needs to be balanced. Therefore, the impact of the new strategy for ultra-rapid LOD on ultra-rapid orbit determination needs to be further analyzed and researched. In addition, typically a combination of GNSS [37] and Satellite Laser Ranging (SLR) observations in the LAser GEOdynamics Satellite (LAGEOS) and Etalon satellites are used to determine LOD [38]. In a limited scope, Doppler Orbitography and Radio-positioning Integrated by Satellite (DORIS) observations can be employed to derive LOD values as well [39]. These observations will be used to enrich the strategy in the future.
On the other hand, it is worth noting that EOP 14 C04 and IERS gpsrapid.daily sometimes have increased latency in practice, which will result in a slower performance of the new strategy. However, as long as the two products are released on time, the performance will be ensured Compared to the LS+AR combination model, the other improved prediction models and methods (e.g., [16,17]) can also be used in this new strategy after appropriate modifications, and the performance should be better. However, this was not assessed and verified here, and future works are required.

Conclusions
In order to explore the potential to improve the accuracy and consistency of ultra-rapid UT1-UTC determination to better serve ultra-rapid orbit determination of GNSS satellites, this study has proposed a new strategy related to ultra-rapid UT1-UTC. Furthermore, to assess the performance of the new strategy, we used UT1-UTC data of EOP 14 C04 provided by IERS, daily rapid products provided by IERS, and ultra-rapid products provided by IGS and iGMAS.
The ultra-rapid UT1-UTC test results of the new strategy showed an overall better performance than the ultra-rapid UT1-UTC test results with the traditional strategy. Compared to the UT1-UTC parameter of IGS ultra-rapid products, the improvements to observation and prediction part accuracies were 42.70% and 29.17%, respectively, while compared to the UT1-UTC parameter of iGMAS ultra-rapid products, the improvements were 42.05% and 37.61%, respectively. Results show that the new strategy can improve the maximum and average of residuals for ultra-rapid orbit determination, but the RMS of residuals are almost equivalent. Additionally, the extended application results show that the new strategy can improve the prior value in ultra-rapid LOD.
Different from the traditional strategy, if the complexity of the new strategy is tolerated during ultra-rapid orbit determination, this strategy can provide an available reference to improve the performance of the ultra-rapid UT1-UTC determination and reduce the impact on coordinate transformation between the celestial (inertial) reference frame and the terrestrial frame in ultra-rapid orbit determination.

Appendix A
Systems Service (IERS), Multi-GNSS Experiment (MGEX), the Crustal Dynamics Data Information System (CDDIS) of the International GNSS Services (IGS), and Chinese iGMAS networks for providing the data and products.

Conflicts of Interest:
The authors declare no conflicts of interest.