Validation of Multi-Year Galileo Orbits Using Satellite Laser Ranging

: Satellite laser ranging (SLR) observations provide an independent validation of the global navigation satellite system (GNSS) orbits derived using microwave measurements. SLR residuals have also proven to be an important indicator of orbit radial accuracy. In this study, SLR validation is conducted for the precise orbits of eight Galileo satellites covering four to eight years (the current longest span), provided by multiple analysis centers (ACs) participating in the multi-GNSS experiment (MGEX). The purpose of this long-term analysis (the longest such study to date), is to provide a comprehensive evaluation of orbit product quality, its inﬂuencing factors, and the effect of perturbation model updates on precise orbit determination (POD) processing. A conventional ECOM solar radiation pressure (SRP) model was used for POD. The results showed distinct periodic variations with angular arguments in the SRP model, implying certain defects in the ECOM system. Updated SRP descriptions, such as ECOM2 or the Box-Wing model, led to signiﬁcant improvements in SLR residuals for orbital products from multiple ACs. The standard deviation of these residuals decreased from 8–10 cm, before the SRP update, to about 3 cm afterward. The systematic bias of the residuals was also reduced by 2–4 cm and the apparent variability decreased signiﬁcantly. In addition, the effects of gradual SRP model updates in the POD were evident in orbit comparisons. Orbital differences between ACs in the radial direction were reduced from the initial 10 cm to better than 3 cm, which is consistent with the results of SLR residual analysis. These results suggest SLR validation to be a powerful technique for evaluating the quality of POD strategies in GNSS orbits. Furthermore, this study has demonstrated that perturbation models, such as SRP, provide a better orbit modeling for the Galileo satellites.


Introduction
Satellite laser ranging (SLR) is a space geodetic technique that offers sub-centimeter single-shot precision for ground site-to-satellite distance measurements. SLR is an independent metrology system that can provide high accuracy measurements to evaluate orbits of the global navigation satellite system (GNSS) obtained from microwave observations. The international laser ranging service (ILRS) is responsible for coordinating global SLR observation sites monitoring multiple navigation satellite systems, including the American global positioning satellite (GPS), Russian GLONASS, European Galileo, and Chinese BeiDou constellation systems [1]. SLR tracking campaigns have proven to be valuable tools for GNSS orbital accuracy assessments.
SLR validation of multi-year orbit products revealed certain characteristics that are not easily identified in short-term data. For example, the analysis of multi-year GPS data has demonstrated that both the standard deviation (STD) and systematic bias, for differences

Availability of Orbital Products
The launch of the Galileo [19]. Due to a technical problem, the first two FOC satellites did not enter the predetermined orbit, but an elliptical orbit with an eccentricity of 0.16 [20].
Four IOV satellites (IOV-1-PRN no. E11, IOV-2-PRN no. E12, IOV-3-PRN no. E19, and IOV-4-PRN no. E20) and four FOC satellites (FOC-3-PRN no. E26, FOC-4-PRN no. E22, FOC-5-PRN no. E24, and FOC-6-PRN no. E30) with the longest orbital time spans (four to eight years) were selected to assess SLR validation for long-term MGEX orbit products. These satellites exhibit near-circular orbits with a semi-major axis of 29,600 km, an inclination of 56 • , and a period of 14.08 h. The launch dates and primary orbital characteristics for these satellites are shown in Table 1. The right ascension of the Remote Sens. 2021, 13, 4634 3 of 15 ascending node, calculated according to orbital position, shows that these eight satellites are distributed in three different orbital planes (A, B, and C). Detailed naming rules for orbital planes and slots can been found in [19]. More intuitive spatial configurations for the eight satellites are shown in Figure 1, in which the two coordinates represent the right ascension of ascending node Ω and the argument of latitude u, respectively. The angle between the two coordinate axes represents the orbital inclination i. The retrograde speed of the ascending node is~10 • /year, due primarily to oblate perturbations from Earth's gravitational field. orbital characteristics for these satellites are shown in Table 1. The right ascension of the ascending node, calculated according to orbital position, shows that these eight satellites are distributed in three different orbital planes (A, B, and C). Detailed naming rules for orbital planes and slots can been found in [19]. More intuitive spatial configurations for the eight satellites are shown in Figure 1, in which the two coordinates represent the right ascension of ascending node Ω and the argument of latitude u, respectively. The angle between the two coordinate axes represents the orbital inclination i. The retrograde speed of the ascending node is ~10°/year, due primarily to oblate perturbations from Earth's gravitational field.  The international GNSS service (IGS) organizes the multi-GNSS experiment (MGEX), which provides high-quality data products for major navigation system satellites [13,14]. In May 2012, major analysis centers (ACs) around the world began providing precise orbital information for navigation satellites. The Galileo precise orbit products included in this study were acquired from the following MGEX ACs: A literature review suggested that, although POD strategies varied, orbital perturbation models at each AC improved significantly during the period ranging from the implementation of MGEX to the end of 2019 (particularly the SRP model). SLR residual time series analysis was divided into two or three stages according to the update time for POD models (see Table 2). The international GNSS service (IGS) organizes the multi-GNSS experiment (MGEX), which provides high-quality data products for major navigation system satellites [13,14]. In May 2012, major analysis centers (ACs) around the world began providing precise orbital information for navigation satellites. The Galileo precise orbit products included in this study were acquired from the following MGEX ACs:

2.
Center for Orbit Determination in Europe (CODE, Switzerland, ID COM).
A literature review suggested that, although POD strategies varied, orbital perturbation models at each AC improved significantly during the period ranging from the implementation of MGEX to the end of 2019 (particularly the SRP model). SLR residual time series analysis was divided into two or three stages according to the update time for POD models (see Table 2).
The radio frequency, differencing, observation sampling, elevation cutoff, ground network, arc length, SRP model, and estimated parameters varied between each AC providing the orbit products [23,25]. As such, multiple products were used for a comparative analysis of accuracy-related factors.  Figure 2 shows the time coverage for eight satellite orbits from five ACs, representing the majority of the satellite lifespan, with some small gaps included. [ 25,26] The radio frequency, differencing, observation sampling, elevation cutoff, ground network, arc length, SRP model, and estimated parameters varied between each AC providing the orbit products [23,25]. As such, multiple products were used for a comparative analysis of accuracy-related factors. Figure 2 shows the time coverage for eight satellite orbits from five ACs, representing the majority of the satellite lifespan, with some small gaps included.

SLR Observation Status
Global SLR sites typically start tracking satellites after they enter predetermined orbits. The ILRS has previously organized three global intensive monitoring campaigns for navigation systems, in order to increase the number of SLR observations available for research. This included the laser ranging to GNSS spacecraft experiment (LARGE), with activities conducted in 2014-2015, 2017, and 2018, each of which spanned several months. As of December 2019, a total of 43 SLR sites have provided worldwide SLR observations from four IOV and four FOC satellites, with time spans ranging from four to eight years.

SLR Observation Status
Global SLR sites typically start tracking satellites after they enter predetermined orbits. The ILRS has previously organized three global intensive monitoring campaigns for navigation systems, in order to increase the number of SLR observations available for research. This included the laser ranging to GNSS spacecraft experiment (LARGE), with activities conducted in 2014-2015, 2017, and 2018, each of which spanned several months. As of December 2019, a total of 43 SLR sites have provided worldwide SLR observations from four IOV and four FOC satellites, with time spans ranging from four to eight years. Fourteen of these sites were selected for this study, based on ILRS system performance data and both the quantity and quality of observations performed at each site. Basic information for the selected sites is shown in Table 3. The geographic distribution of these sites was relatively uniform, as shown in Figure 3. SLR observation data were provided by two data centers, the American Crustal Dynamics Data Information System (CDDIS) and the EUROLAS Data Center (EDC). Monthly quantities, global shares, and the proportion of daytime observations at selected sites are shown in Figure 4. Fourteen of these sites were selected for this study, based on ILRS system performance data and both the quantity and quality of observations performed at each site. Basic information for the selected sites is shown in Table 3. The geographic distribution of these sites was relatively uniform, as shown in Figure 3. SLR observation data were provided by two data centers, the American Crustal Dynamics Data Information System (CDDIS) and the EUROLAS Data Center (EDC). Monthly quantities, global shares, and the proportion of daytime observations at selected sites are shown in Figure 4.  As seen in the Figure 4, the monthly average for SLR observations at 14 selected sites mostly exceeded 200 normal points. This accounted for more than 80% of observations at each of the 43 sites, as shown in Table 4. The proportion and quality of observations at selected locations were high and the geographical distributions were relatively uniform. More than half of the selected sites offered strong daytime observation capabilities. Daytime measurements accounted for ~25-40% of total tracking data, of which IOV satellites constituted a higher percentage than FOC satellites. SLR data from the selected sites were used for subsequent analysis.    As seen in the Figure 4, the monthly average for SLR observations at 14 selected sites mostly exceeded 200 normal points. This accounted for more than 80% of observations at each of the 43 sites, as shown in Table 4. The proportion and quality of observations at selected locations were high and the geographical distributions were relatively uniform. More than half of the selected sites offered strong daytime observation capabilities. Daytime measurements accounted for~25-40% of total tracking data, of which IOV satellites constituted a higher percentage than FOC satellites. SLR data from the selected sites were used for subsequent analysis.

Angular Argument Definitions
As discussed previously, understanding the relationship between SLR residuals and angular arguments is of great importance for the analysis of orbital accuracy in the SRP model. As such, it is necessary to clarify the definitions of relevant angular arguments, four of which are discussed below.

1.
Solar elevation over the satellite orbital plane (β)-This term represents the elevation angle for the Sun-Earth connection with respect to the satellite orbital plane. In the Galileo system, β is a slow variable that completes approximately one cycle each year, but will reach its zero value twice from the positive and negative directions. When |β| is less than the critical value β c = asin(R E /r), where R E is the radius of the Earth and r is the satellite orbital radius, the satellite will enter an eclipse season. During this period, part of the satellite orbital motion will experience a shadow from the Earth. This angle definition has been used in prior SRP models [24] and exhibits a maximum allowable range of [−23.5 • , 23.5 • ].

2.
The satellite-Sun elongation angle (γ)-This term describes the elongation of positions of the Earth and the Sun relative to the satellite. In this study, γ is a fast variable that completes approximately one cycle in each satellite orbital period. It is only approximate because the position of the Sun will change slightly during one satellite orbital period. This angle definition has been used in previous SRP models, such as the ROCK [26] and Box-Wing (BW) models [13] and it exhibits a maximum allowable The argument of latitude (u)-This term denotes the angle between the satellite position and the ascending node of the satellite orbit. Here, u is a fast variable that completes exactly one 360 • cycle in each satellite orbital period. This definition has been used in prior SRP models, such as ECOM [5].

4.
The argument of satellite latitude relative to the Sun (∆u)-This argument is defined as the difference between the argument of latitude (u) of the satellite and the solar argument of latitude in the satellite orbital plane (u 0 ). The angle ∆u = u − 0 is a fast variable that completes approximately one 360 • cycle in each satellite orbital period. This definition has been used in previous SRP models, such as ECOM2 [10]. The geometric meaning of each angle is represented in Figure 5.
used in prior SRP models, such as ECOM [5]. 4. The argument of satellite latitude relative to the Sun (Δu)-This argument is defined as the difference between the argument of latitude (u) of the satellite and the solar argument of latitude in the satellite orbital plane (u0). The angle Δu = u − 0 is a fast variable that completes approximately one 360° cycle in each satellite orbital period. This definition has been used in previous SRP models, such as ECOM2 [10]. The geometric meaning of each angle is represented in Figure 5.

SLR Residuals for COM Orbits
Precise orbits and SLR normal point data from the longest spans (as of December 2019) were selected for Galileo orbit validation. The orbital SLR residual time series for two IOV and two FOC satellites are shown for COM orbits in Figure 6, the evolution of which can be divided into three stages. Residuals in the first stage (prior to 4 January 2015) showed a non-zero mean value with a systematic bias of −6 cm, an amplitude of ~20 cm, and a standard deviation (STD) of ~8 cm. These residual terms improved significantly in the second stage (between 4 January 2015 and 6 August 2017), exhibiting a bias of roughly −4 cm, an STD of only ~4 cm, and an amplitude of less than 10 cm. In the third stage (after 6 August 2017), the bias and STD were further reduced to approximately −1 cm and 3 cm, respectively.

SLR Residuals for COM Orbits
Precise orbits and SLR normal point data from the longest spans (as of December 2019) were selected for Galileo orbit validation. The orbital SLR residual time series for two IOV and two FOC satellites are shown for COM orbits in Figure 6, the evolution of which can be divided into three stages. Residuals in the first stage (prior to 4 January 2015) showed a nonzero mean value with a systematic bias of −6 cm, an amplitude of~20 cm, and a standard deviation (STD) of~8 cm. These residual terms improved significantly in the second stage (between 4 January 2015 and 6 August 2017), exhibiting a bias of roughly −4 cm, an STD of only~4 cm, and an amplitude of less than 10 cm. In the third stage (after 6 August 2017), the bias and STD were further reduced to approximately −1 cm and 3 cm, respectively. Modifications to the Galileo satellite POD strategy caused significant changes in SLR validation. Prior to 4 January 2015, COM orbit products used the ECOM SRP model. The updated ECOM2 SRP model was then adopted. This update had obvious effects on improvements to orbital accuracy, which reduced residuals by 50%. After 6 August 2017, updated Earth albedo radiation pressure and satellite antenna thrust models were used in the orbit determination process [17], leading to evident improvements in the systematic bias of SLR residuals (~3 cm). Table 5 summarizes the validation results for three-stage Modifications to the Galileo satellite POD strategy caused significant changes in SLR validation. Prior to 4 January 2015, COM orbit products used the ECOM SRP model. The updated ECOM2 SRP model was then adopted. This update had obvious effects on improvements to orbital accuracy, which reduced residuals by 50%. After 6 August 2017, updated Earth albedo radiation pressure and satellite antenna thrust models were used in the orbit determination process [17], leading to evident improvements in the systematic bias of SLR residuals (~3 cm). Table 5 summarizes the validation results for three-stage COM orbits from all eight Galileo satellites. The red curve in Figure 6 represents changes in the Sun elevation angle (β) over time, with the grey area indicating the eclipse season. It is evident from the plot that changes in the amplitude of SLR residuals follow a synchronization pattern for satellites in the same orbital plane (e.g., E12 and E26 in plane B), since variations in the angle β are temporally synchronized. Orbital accuracy is also low, while SLR residuals are high during the eclipse season. SLR residuals could exhibit a local maximum when the angle β is close to 0 • (i.e., when the Sun passes through the satellite orbital plane). This is possibly due to sharp increases in the yaw rate of the Galileo satellites (when β = 0 • ), inducing attitude control system mechanisms and corrective maneuvers. This is reflected in both the SRP model and the estimation of D0 (a fitting coefficient in the satellite-Sun vector direction). It also produced local abnormalities [13], which in turn caused satellite orbital errors to increase in the radial direction and eventually led to local increases in SLR residual amplitudes, which are typically small when the angle β is large. In contrast, these residuals were only 25-50% of their value during the eclipse season (when β was at its maximum).
The SLR residuals of satellites in different orbital planes showed similar changes. In addition, the eclipse season and maximum value of residuals manifested in different temporal regions, due to varying trends in the angle β with time. SLR residuals also became less dependent on β from the first to the third stage of the orbit determination strategy. However, in the third stage, the residuals from some orbit products still exhibited a weak dependence on β, suggesting imperfections in current orbital perturbation models and a need for further refinement.

SLR Residuals for GBM Orbits
The residual time series for GBM orbit products from four satellites is shown in Figure 7. This sequence was divided into two stages using 25 October 2016 as the boundary. The ECOM SRP model was applied in the first stage [8] and an a priori Box-Wing (BW) model was adopted (in addition to ECOM) in the second stage [27]. The orbital SLR residual STD in Stage 2 (after 25 October) was significantly smaller than that of Stage 1. In addition, the dependence of Stage 2 residuals on the angle β was nearly eliminated, indicating the SRP model update was successful. Systematic bias in GBM residuals was also reduced significantly after adopting the a priori BW model, indicating a partial correction Remote Sens. 2021, 13, 4634 9 of 15 of model errors caused by albedo radiation pressure from the Earth and satellite antenna thrust. Table 6 summarizes these validation results for GBM orbits.
The residual time series for GBM orbit products from four satellites is shown in Fig-ure 7. This sequence was divided into two stages using 25 October 2016 as the boundary. The ECOM SRP model was applied in the first stage [8] and an a priori Box-Wing (BW) model was adopted (in addition to ECOM) in the second stage [27]. The orbital SLR residual STD in Stage 2 (after 25 October) was significantly smaller than that of Stage 1. In addition, the dependence of Stage 2 residuals on the angle β was nearly eliminated, indicating the SRP model update was successful. Systematic bias in GBM residuals was also reduced significantly after adopting the a priori BW model, indicating a partial correction of model errors caused by albedo radiation pressure from the Earth and satellite antenna thrust. Table 6 summarizes these validation results for GBM orbits.

Radial Differences in the COM and GBM Orbits
In the Galileo system observation geometry, the line of sight from the ground to satellite is always at an angle of less than 14 • from the radial orbit direction. In addition, SLR residuals were highly consistent with orbital clock corrections, which map orbital radial errors [13,22]. This suggests that SLR residuals can generally be regarded as a measure of radial errors for Galileo satellites, which were investigated in this study using orbital comparisons. Figure 8 shows radial differences in the COM and GBM orbits of the four Galileo satellites. This comparison also demonstrates an obvious dependence on β and indicates that orbital errors increased significantly during eclipse season.
ellite is always at an angle of less than 14° from the radial orbit direction. In addition, SLR residuals were highly consistent with orbital clock corrections, which map orbital radial errors [13,22]. This suggests that SLR residuals can generally be regarded as a measure of radial errors for Galileo satellites, which were investigated in this study using orbital comparisons. Figure 8 shows radial differences in the COM and GBM orbits of the four Galileo satellites. This comparison also demonstrates an obvious dependence on β and indicates that orbital errors increased significantly during eclipse season. This orbit comparison also demonstrates the existence of three distinct time segments, separated by boundaries occurring on 25 October 2016 (when the GBM SRP model was updated) and 6 August 2017 (the Earth albedo radiation pressure and satellite antenna thrust model updates for COM orbit products). The STD of differences in radial components were ~10 cm in Segment 1, generally reduced below 5 cm in Segment 2, and further improved to less than 3 cm in Segment 3. The absolute value of systematic biases along orbital radial directions was also reduced from 4 cm to nearly less than 1 cm, representing a significant improvement. Table 7 summarizes the statistical results of orbit This orbit comparison also demonstrates the existence of three distinct time segments, separated by boundaries occurring on 25 October 2016 (when the GBM SRP model was updated) and 6 August 2017 (the Earth albedo radiation pressure and satellite antenna thrust model updates for COM orbit products). The STD of differences in radial components were~10 cm in Segment 1, generally reduced below 5 cm in Segment 2, and further improved to less than 3 cm in Segment 3. The absolute value of systematic biases along orbital radial directions was also reduced from 4 cm to nearly less than 1 cm, representing a significant improvement. Table 7 summarizes the statistical results of orbit comparison for each of these three segments. Comparing Table 7 with Tables 5 and 6 suggests that SLR residuals are in good agreement with radial orbit errors.

SLR Residuals Results for WUM, GRM, and TUM Orbits
In addition to the residual validation results for COM and GBM orbits listed above, Tables 8-10 summarize results for WUM, GRM, and TUM products. Improvements in the orbit product quality, due to SRP model updates, were reflected in the SLR validation process as the STD of residuals generally decreased from about 10 cm to 3-4 cm. The absolute value of systematic biases also improved, decreasing by 2-4 cm from the initial to the latter stages. In addition, the bias in some orbit products was close to zero and the residual STD using the ECOM2 model was significantly lower than that of the ECOM model. Orbits processed using only ECOM2 exhibited small but obvious bias, while orbits processed using both BW and ECOM models showed significantly lower bias.    Figure 9 shows the relationship between SLR residuals for COM and GBM orbits and the satellite-Sun elongation angle γ at each stage. These residuals showed a pronounced correlation with γ before the SRP model update and were typically positive when the satellite was near the Sun (compared to the Earth), corresponding to γ > 90 • . Conversely, residuals were typically negative when the satellite was in the far-Sun position (γ < 90 • ). This correlation between SLR residuals and the angle γ was reduced significantly after the SRP model update. The combined BW-ECOM model used for GBM was more effective in reducing the correlation with orbital residuals than that of the ECOM2 model for COM. The dispersion and systematic residual bias were further reduced when the COM orbit was supplemented with albedo radiation pressure and antenna thrust model updates. Figure 9 shows the relationship between SLR residuals for COM and GBM orbits and the satellite-Sun elongation angle γ at each stage. These residuals showed a pronounced correlation with γ before the SRP model update and were typically positive when the satellite was near the Sun (compared to the Earth), corresponding to γ > 90°. Conversely, residuals were typically negative when the satellite was in the far-Sun position (γ < 90°). This correlation between SLR residuals and the angle γ was reduced significantly after the SRP model update. The combined BW-ECOM model used for GBM was more effective in reducing the correlation with orbital residuals than that of the ECOM2 model for COM. The dispersion and systematic residual bias were further reduced when the COM orbit was supplemented with albedo radiation pressure and antenna thrust model updates.  Figure 10 shows the relationship between SLR residuals and the satellite argument of latitude u. Prior to the SRP model update, residuals from both COM and GBM orbit products varied with the angle u (using the ECOM model). This period of fluctuation was Figure 9. SLR residuals as a function of the Earth-satellite-Sun elongation γ.   Figure 11. As a result, residuals were mostly positive, with a maximum value of ~15-20 cm, when Δu was close to 0° or 360° (i.e., the satellite and the Sun were aligned in approximately the same direction). In contrast, residuals were mostly negative, with a minimum value of roughly −20 cm to −25 cm, when Δu was close to 180° (i.e., the satellite and the Sun were aligned in approximately opposite directions relative to the Earth). Some studies attributed the relationship between the sign of the SLR residuals (positive or negative) and the angle Δu to surface area characteristics on the X-and Z-side of the satellite body [18,28,29]. This correlation decreased significantly after the SRP model update. Figure 10. SLR residuals as a function of the argument of latitude u. Prior to the update, SLR residuals also exhibited a clear dependence on the argument of satellite latitude with respect to the Sun (∆u), as shown in Figure 11. As a result, residuals were mostly positive, with a maximum value of~15-20 cm, when ∆u was close to 0 • or 360 • (i.e., the satellite and the Sun were aligned in approximately the same direction). In contrast, residuals were mostly negative, with a minimum value of roughly −20 cm to −25 cm, when ∆u was close to 180 • (i.e., the satellite and the Sun were aligned in approximately opposite directions relative to the Earth). Some studies attributed the relationship between the sign of the SLR residuals (positive or negative) and the angle ∆u to surface area characteristics on the X-and Z-side of the satellite body [18,28,29]. This correlation decreased significantly after the SRP model update.

Relationship between Orbital SLR Residuals on the Argument of Satellite Latitude
Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 17 Figure 11. SLR residuals as a function of relative Sun-satellite latitude Δu.

Dependence of Orbital SLR Residuals on the Local Time of Observation
The dependence of orbital SLR residuals on the local time of observation is shown in Figure 12. The satellite and Sun were both above the horizon of the tracking site when conducting daytime observations (i.e., on the same side of the Earth). This was similar to the relative positions of the satellite and the Sun when |Δu−180°| > 90° and γ > 90°. In this configuration, most of the orbital residuals are positive. During nighttime observations, the position of the Sun is on the other side of the tracking site, similar to the relative position of the satellite and the Sun when |Δu−180°| < 90° and γ < 90°. In this configuration, most of the residuals are negative. Residual variations with a site's local time were reduced significantly after the updates. This is similar to the relationship between orbital SLR residuals and the angles Δu or γ. Daytime observations accounted for ~25-40% of the total data, as mentioned previously. In other words, the share of day and night observations was roughly equal and such residual distributions could not be attributed to observation time. These results suggest that diurnal variations in the SLR residuals from Galileo satellite orbits are strongly correlated with the quality of SRP modeling. In fact, the included model appears to be the primary factor affecting the dependence of orbital residuals on local time, as opposed to differences in the accuracy of SLR observations during both daytime and nighttime measurements. Figure 11. SLR residuals as a function of relative Sun-satellite latitude ∆u.

Dependence of Orbital SLR Residuals on the Local Time of Observation
The dependence of orbital SLR residuals on the local time of observation is shown in Figure 12. The satellite and Sun were both above the horizon of the tracking site when conducting daytime observations (i.e., on the same side of the Earth). This was similar to the relative positions of the satellite and the Sun when |∆u−180 • | > 90 • and γ > 90 • . In this configuration, most of the orbital residuals are positive. During nighttime observations, the position of the Sun is on the other side of the tracking site, similar to the relative position of the satellite and the Sun when |∆u−180 • | < 90 • and γ < 90 • . In this configuration, most of the residuals are negative. Residual variations with a site's local time were reduced significantly after the updates. This is similar to the relationship between orbital SLR residuals and the angles ∆u or γ. Daytime observations accounted for~25-40% of the total data, as mentioned previously. In other words, the share of day and night observations was roughly equal and such residual distributions could not be attributed to observation time. These results suggest that diurnal variations in the SLR residuals from Galileo satellite orbits are strongly correlated with the quality of SRP modeling. In fact, the included model appears to be the primary factor affecting the dependence of orbital residuals on local time, as opposed to differences in the accuracy of SLR observations during both daytime and nighttime measurements.

Conclusions
In this study, SLR observations were used to validate the accuracy of orbit products for Galileo satellites with the longest temporal spans currently provided by MGEX analysis centers. The characteristics of these residual time series were analyzed, with orbital comparisons used as auxiliary support for SLR validation. The relationship between residuals and solar elevation, satellite-Sun elongation angle, the argument of satellite latitude, the argument of satellite latitude with respect to the Sun, and local observation time

Conclusions
In this study, SLR observations were used to validate the accuracy of orbit products for Galileo satellites with the longest temporal spans currently provided by MGEX analysis centers. The characteristics of these residual time series were analyzed, with orbital comparisons used as auxiliary support for SLR validation. The relationship between residuals and solar elevation, satellite-Sun elongation angle, the argument of satellite latitude, the argument of satellite latitude with respect to the Sun, and local observation time were discussed in detail. The results showed that multi-year SLR validation is of great benefit for accuracy assessment in Galileo satellite orbits. The presented methodology is also critical for evaluating the quality of orbit determination models such as SRP and other POD perturbation modeling strategies. As demonstrated, the SRP update significantly improved the quality of Galileo satellite orbit products. The STD of current SLR residuals, produced using updated perturbations, was reduced from 8-10 cm in the early data products to 3-4 cm in the recent products. The systematic bias was also reduced significantly (~2-4 cm). However, the Galileo SRP model and other perturbation techniques should be improved and evaluated further.
Both the quantity and quality of SLR tracking data have an obvious effect on the orbit evaluation. As such, this study only selected site data exhibiting sufficiently high observation quality and uniform spatial distributions for orbital SLR validation. Even so, at SLR sites using multi-photon detectors, beams along the direction of inclination for the on-board retro-reflector array can cause asymmetric broadening of the echo pulse. This produced a functional relationship between the mean value of SLR residuals and the satellite nadir angle, known as the SLR signature effect [30], which resulted in a difference of~1 cm between the bias at single-photon and multi-photon sites [18]. This small discrepancy does not affect the overall conclusion of the study, but it is likely worthy of further investigation. It should also be given more attention in the process of orbital SLR validation, due to recent improvements in GNSS orbital accuracy.
Author Contributions: E.T. performed experiments, analyzed the data, and prepared the manuscript. X.Z. and N.G. provided crucial guidance and made important suggestions on data processing. K.X. and B.W. revised the papers. All authors have read and agreed to the published version of the manuscript.