Determination of global geodetic parameters using satellite laser ranging to Galileo, GLONASS, and BeiDou satellites

The new Global Navigation Satellite System (GNSS) satellites, including GLONASS, Galileo, and BeiDou system, are equipped with Laser Retroreflector Arrays (LRA) to support Satellite Laser Ranging (SLR) tracking, which contributes to the estimation of global geodetic parameters. In this study, we estimate the global geodetic parameters using the SLR observations to GNSS satellites and also investigate the effects of different data processing strategies on the estimated Earth Rotation Parameters (ERP), geocenter motion, and terrestrial scale. The results indicate that setting range bias parameters for each satellite-station pair can effectively account for the satellite-specific biases induced by LRAs, leading to smaller Root Mean Square Errors (RMSE) of the post-fit SLR residuals. Furthermore, estimating the range biases for each satellite-station pair improves the accuracy of the estimated station coordinates and ERP. We also examine the impact of different arc lengths on the estimates of ERP, geocenter motion, and terrestrial scale. The results show that extending arc length can significantly reduce the formal error of ERP. The 7-day strategy produces the smallest RMSEs of 473 microarcseconds and 495 microarcseconds for the estimated X - and Y -component of pole coordinates, and 52 microseconds for length-of-day, respectively. However, the estimated geocenter motion is less affected by the arc length, even the shortest 1-day arc strategy can capture the seasonal variations of geocenter motion in Z component. For scale estimation, extending the arc length notably improves the accuracy of the estimated station coordinates and scale, but this advantage becomes less noticeable in longer arcs. The 7-day solution also obtains the closet scale results compared to ITRF2014, with the RMSE of 2.10 × 10 –9 .


Introduction
Satellite Laser Ranging (SLR) stands out as a groundbased optical space geodetic technique which is not affected by the ionosphere (Pearlman et al., 2019).Leveraging its remarkable observation accuracy, SLR plays a pivotal role in various domains, such as the determination of the orbits and validations for Global Navigation Satellite System (GNSS) and Low Earth Orbit (LEO) satellites (Hackel et al., 2015;Prange et al., 2017;Sośnica et al., 2018a;Strugarek et al., 2022;Urschl et al., 2007) and the global geodetic parameters including Earth Rotation Parameters (ERP), geocenter motion, and scale (Bruni et al., 2018;Bury et al., 2019a;Li et al., 2021Li et al., , 2023a)).Due to its sensitivity in estimating geocenter motion and terrestrial scale, SLR is also used for the realization of the origin and scale of the International Terrestrial Reference Frame (ITRF) (Altamimi et al., 2016).
The global SLR tracking activities are managed and distributed by the International Laser Ranging Service (ILRS) (Pearlman et al., 2019).In the early stages, the ILRS network primarily tracked geodetic satellites such as LAGEOS-1 and LAGEOS-2 for the determination of the global geodetic parameters (Appleby, 1998;Pearlman et al., 2002;Sośnica et al., 2019).With more LEO and GNSS satellites equipped with SLR prism reflectors and Laser Retroreflector Arrays (LRA), ILRS has expanded its tracking capabilities, encompassing more LEO satellites like Etalon and Swarm, and GNSS satellites such as GPS, GLONASS, Galileo, and BeiDou Navigation Satellite System (BDS).For GNSS satellites, ILRS has organized campaigns to increase the SLR tracking missions (Bury et al., 2019b;Pearlman et al., 2015).Currently, more than 60 active GNSS satellites have been tracked by ILRS stations (https:// ilrs.gsfc.nasa.gov/ missi ons/ satel lite_ missi ons/), including the new Galileo satellites and BDS-3 satellites.The total number of SLR observations to GNSS satellites has exceeded that of LAGEOS.Over the period from 2019 to 2021, the number of normal points to LAGEOS-1/2 approached 400,000, whereas the total normal points to GLONASS, Galileo, and BDS exceeded 600,000.Additionally, some SLR stations such as Altay (1879, ALTL) and Baikonur (1887, BAIL) have also shown a predominant focus on tracking GNSS satellites (Bury et al., 2021), which will enhance the effectiveness of SLR observations to GNSS in certain scenarios.
SLR observations to GNSS satellites are also used in the estimation of global geodetic parameters.Thaller et al. (2011) indicated that for the stations with reasonable number of observations, the combined GNSS and SLR solution exhibits the agreement at the level of 2 cm for SLR coordinates with the a priori TRF, and 1-2 cm for the GNSS coordinates corrected by local ties.Sośnica et al. (2018b) exhibited that the combination of SLR observations to GNSS and LAGEOS obtained an improvement of 79.5 microseconds in Length of Day (LoD) Root Mean Square Errors (RMSEs) compared to the LAGEOS only solution.The improvement is more significant than those in combination of the observations to Starlette, Stella, and Ajisai satellites (Sośnica et al. 2014a) or those gained by building a new station (Otsubo et al., 2016).However, the SLR observations to GNSS satellites used in most previous studies were limited.On average, the number of SLR observations to GLONASS, Galileo, and BDS systems in 2020 witnessed an increase of over 700 per week compared to that in 2018.Furthermore, the number of active satellites before 2017 was comparatively smaller.The SLR observations to GNSS in current period hold notable potential for the determination of global geodetic parameters.Appropriate utilization of these observations can improve the accuracy of the estimated global geodetic parameters.On the other hand, the range biases, as one of the major error sources in SLR processing (Drożdżewski & Sośnica, 2021;Thaller et al., 2015), can be corrected with prior values and estimated as station-, system-and satellite-specific biases (Appleby et al., 2016;Guo et al., 2018;Sośnica et al., 2014a;Sośnica et al., 2019).Proper handling of range biases can effectively improve the estimates of station coordinates and geodetic parameters.Sośnica (2014b) exhibited that estimating range bias can entirely absorb the biases of about − 0.5 mm in height component caused by the Earth radiation pressure on the station coordinate estimation and compensate for the deficiencies in center-of-mass corrections for LEO satellites.Luceri et al. (2019) proved that estimating range biases can obtain the scale realized by SLR closer to the scale of ITRF by ~ 1 × 10 -9 .Sośnica et al. (2019) exhibited the estimated ERPs by employing the SLR data to GNSS.Fixed range biases by the annual value obtained the X-and Y-pole RMSE with respect to IERS EOP 14C04 (Bizouard et al., 2019), on average 35 microarcseconds smaller than those estimated in systemstation pair, and 124 microarcseconds smaller than those estimated in satellite-station pair.Therefore, it is highly significant to explore the determination of consequent global geodetic parameters using the SLR observations to GNSS satellites in the current scenario.
In this study, we determined ERP, geocenter motion, and scale using the SLR observations to 22 Galileo, 20 GLONASS, 4 BDS-2, and 4 BDS-3 satellites between 2019 and 2021.The model of SLR observations to GNSS satellites is first demonstrated.Secondly, we introduce the data and process strategies.Then, we present the results of geodetic parameters estimation results using SLR observations, followed by the conclusions.

Method
The observations of SLR can be expressed as follows: where c is the velocity of light, τ is the two-way travel time of the laser, d s r is the one-way distance between the station receiver and the center of satellite LRA, d trop is the troposphere delay while d rel is the relativity effect, d LRA represents the range correction from the center of the LRA to the satellite center of mass, d RB is range bias which is one of the major errors of SLR process, and ε SLR represents the other unmodeled errors and random errors in the observation.
The models of troposphere delay (Mendes & Pavlis, 2004), relativity effect (Petit and Luzum, 2010), LRA offsets correction (https:// ilrs.gsfc.nasa.gov/ missi ons/ satel lite_ missi ons/ index.html), and post-seismic (1) (Altamimi et al., 2016) are used, while the range bias is estimated.By denoting the oneway range observation cτ 2 with the bias corrections as RSLR , and expanding the one-way distance d s r as the function of SLR station and satellite positions, the observation model can be expressed as: where ψ s r is the direction cosines, �(t, t 0 ) s is the state- transition matrix from the reference epoch t 0 to the measuring epoch t which is obtained from orbit integration, o 0 s represents the initial status of satellite at the time t 0 , R represents the rotation matrix from terrestrial reference frame to celestial reference frame, r r is the coordinates of SLR stations, G represents the distance vector of the instantaneous Earth center of mass relative to the center of the reference frame.

Dataset
The ILRS presently encompasses a global SLR station network to monitor the satellite orbits and sample global geodetic parameters properly (Pearlman et al., 2019).We employ 29 stations which track GLONASS, Galileo, or BDS satellites, including 4 receiver types: 5 Micro-Channel Plate (MCP) stations, 13 photo-multiplier tube (PMT) stations, 4 Single Photon Avalanche Diode (SPAD) stations, and 7 Compensated Single Photon Avalanche (2) Diode (C-SPAD) stations which are the improved version of the SPAD. Figure 1 shows the distribution of these stations, which are mainly located in the northern hemisphere.Table 1 provides the information about their monuments, codes, types, and locations.The further information can be obtained from the ILRS website (https:// ilrs.gsfc.nasa.gov/ netwo rk/ stati ons/ active/ index.html).
We employ the ILRS normal points data to 22 Galileo, 20 GLONASS, 4 BDS-2, and 4 BDS-3 satellites during 2019-2021 for the geodetic parameter estimation.Considering that not all observations have high ranging accuracy, we selected the normal points data of these satellites.We obtained the positions of satellites using the results from a prior microwave-based precise orbit determination and calculated the corresponding satellite-station ranges at the observation epochs.After employing the delays and corrections for SLR observations, we screen away data whose difference between observation and satellite-station range is greater than 0.15 m.Overall, 1.05% of the data were screened away from for the period of three years.The daily count of selected observations for each satellite can be observed in Fig. 2. We used Pseudo-Random Noise (PRN) to represent satellites in the axes, where PRN ranges from E01 to E36 for Galileo, R01 to R24 for GLONASS, and C08 to C30 for BDS.noise level, leading to a reduction in the number of available observations which can be observed in Fig. 2 (Thaller et al., 2012).

Data processing strategies
As one of the primary error sources in observations, SLR range biases encompass the photoelectric hardware delay and other unmodeled errors that occur in station or satellite such as the orbits errors and the signature effect (Otsubo et al., 2001;Sośnica et al., 2018a).The strategies to estimate range biases will directly affect the estimates.In the previous studies, range bias is usually estimated for each system or fixed to an annual averaged value due to the limited number of SLR observations to GNSS satellites (Bury et al., 2019a;Glaser et al., 2015;Thaller et al., 2011).However, these are the compromise solutions to address the issue of limited number of observations.In these cases, satellite-specific characteristics of range biases are neglected, which may eventually lead to unexpected errors in SLR solutions.In recent years, the SLR tracking to GNSS satellites has been rapidly growing, making it possible to estimate the range bias for each satellite-station pair.In this study, we used five strategies to investigate the optimal range bias process strategies as listed in Table 2. Range biases have been individually addressed for each station.In the SYS solution, we set a common range bias parameter to be estimated for each system-station pair.The CON solution, as an extension of the SYS solution, introduces additional initial values and constraints for the range bias estimation.These initial values and constraints are derived from the annual average values and Standard Deviation (STD) of the estimates obtained in the SYS solution.In the FIX_SYS solution, we utilize the annual average values obtained from the SYS solution to fix the range biases for each system, thereby reducing the number of estimated parameters.In the SAT solution, we estimate range biases for each satellitestation pair, which can better consider the satellite-specific characteristics, but resulting in the most estimated parameters.In the FIX_SAT solution, the range biases for each satellite are fixed using the annual average values derived from the SAT solution.The effects of different range bias process strategies are given in "Effects of different range bias processing strategies" section.
On the other hand, the arc length is another important factor in SLR processing, can multiply the number of observations in the process and directly affect the precision of the estimation results.And the determination of global geodetic parameters with SLR is typically conducted using long arcs, such as the 7-day arc (Sośnica et al., 2018b;Thaller et al., 2014), to ensure adequate observations in the process.However, with an increase in the number of observations, the short arc strategy may also yield the results comparable to the long arc strategy.We design four schemes with the arc length of 1-day, 3-day, 5-day, and 7-day to explore the effect of arc length scheme.The results are given in "Determination of global geodetic parameters" section.
The SLR analysis and parameter estimation in this study are completed with the GNSS + Research, Application and Teaching (GREAT) software (Li et al., 2021).Table 3 shows the processing strategies for the global geodetic parameter estimation with GNSS SLR observations.
The GNSS satellite orbits are fixed to the prior microwave-based solution which are also obtained by the GREAT software.Table 4 shows the processing strategies for the prior GNSS satellite orbits.These orbits are expressed in the celestial reference frame and have not been transformed into the terrestrial reference frame.
Table 5 presents the average RMSE of orbit differences in along-, cross-track, radial directions and 3D directions between the GNSS orbits and the products of the Multi-GNSS Experiment (MGEX) products from the Center for Orbit Determination in Europe (CODE) and Wuhan University.These results show that the discrepancies between the orbits calculated by GREAT software and the products provided by MGEX institutions are smaller than 6 cm, 9 cm, 10 cm and 15 cm for Galileo, GLONASS, BDS MEO and BDS IGSO in 3D direction, which provide sufficient supports for the SLR process in this study.

Effects of different range bias processing strategies
Figure 3 displays the range bias estimates derived from SAT and SYS solutions for six SLR stations.Stations YARL, CHAL, and GRZL have the most observations during the period from 2019 to 2021.Specifically, YARL recorded over 60,000 observations over the three-year period, CHAL recorded approximately 50,000 observations, and GRZL recorded nearly 30,000 observations.In contrast, stations STL3, BEIL, and MONL exhibit comparatively lower observation counts over the same period, amounting to approximately 14,000, 10,000, and 3,000, respectively.Overall, the stations with more      However, it is noteworthy that the satellites with limited observations, such as R01 and C08, averaging only four observations per day, may exhibit relatively higher STD values in the estimated range biases when employing the SAT solution.
Notably, for station YARL, most Galileo satellites exhibit the mean range bias around 0 mm in SAT solution, whereas the mean range bias in SYS solution is 17.6 mm.This difference could be attributed to E12 and E19, which have more observations and smaller average values, introduce a systematic bias affecting other satellites range biases and influence the outcome of SYS solution.The same phenomenon can also be observed in BDS.The range biases obtained in SYS solution have an average of approximately 20 mm, but in SAT solution, the C21 and C20 satellites manufactured by China Academy of Space Technology (CAST) exhibit an average range bias of about 40 mm, while C29 and C30 satellites manufactured by Shanghai Engineering Center for Microsatellites (SECM) show an opposite value of about -40 mm.Although these four satellites are all BDS-3 MEO satellites, the differences in LRA sizes, prisms, and optical coatings could lead to different performances (Li et al., 2023b).In addition, other satellite-specific errors such as orbital errors can also be absorbed into the range bias.Therefore, estimating range biases as a common parameter for the satellites within one system will neglect the errors caused by satellite-specific differences and lead to significant systematic errors.
Figure 4 displays the post-fit SLR residuals of Galileo satellites at station YARL in the three years with the five range bias strategies.The nadir angle on the X axis represents the laser incidence angle with respect to the normal vector of the LRA.Slopes of about -0.4 mm/° can be observed in the solutions except the SAT.These residual variations with nadir angle are caused by the signature effect, which is a representation of the range deviation that occurs due to the reflection of photons at different locations of the LRAs.When the laser is not perpendicular to the LRA flat, the photons reflected at the edge of the plate will return earlier and result in a shorter ranging time, leading to a ranging error.This effect is theoretically a function of the satellite nadir angle (Otsubo et al., 2001), which also differs from LRAs.As the FIX_SAT solution fixes range biases using the annual average values derived from SAT estimates, it provides a degree of absorption of the biases arising from satellite characteristics.Nevertheless, given the varying incidence angles of the laser on the LRA during different arcs, the time-varying component of the signature effect remains inadequately corrected.Consequently, in FIX_SAT, the residuals also decrease as the nadir angle increases.A similar situation appears in the FIX_SYS solution, and due to its inadequate correction for satellite-specific biases, FIX_SYS yields a larger RMSE.In stations SYS and CON, estimating range biases in each system-station pair and each arc can absorb the errors that are common within the system in each processing arc, which leads to smaller RMSE, but they still fail to address the satellite-specific errors.Only estimated range bias in each satellite-station pair and in each arc can better account for the satellite-specific biases, as well as the time-varying component of signature effect among different arcs.Therefore, the SAT solution can most effectively resolve the influence of the signature effect and obtain the residuals with smaller RMSE.
Figure 5 and Table 6 show the station coordinates with respect to the SLRF2014 which is the ILRS realization of ITRF2014, employing five range bias strategies.Both for the stations with sufficient observations and for the stations with relatively limited numbers of observations, the coordinate RMSEs obtained in the SAT and two FIX solutions exhibit greater proximity and alignment with those in SLRF2014, in contrast to the other two solutions.This shows that a reduction in estimated parameters or estimating range bias for each satellite can both enhance the precision of the station coordinate results.Simultaneously, within the three strategies, the RMSEs, whiskers, and interquartile range (IQR) values derived from the SAT and FIX_SAT strategies exhibit greater similarity.In comparison to these two approaches, the RMSE values obtained by FIX_SYS are on average 4 mm larger in the The box plots of up component errors for 6 stations in 2019-2021, with different range bias process strategies and the 1-day/3-day arc length strategies three components.Moreover, the median values in FIX_ SYS exhibit slight biases that align more closely with the SYS strategy.From the perspective of station coordinate estimation, unaddressed satellite-specific biases are more likely to impact the accuracy of the station coordinates compared to the estimation of range biases.
For the SYS and CON solutions, the compromise between accounting specific errors and reducing parameter number actually leads to more scattered coordinate results.Furthermore, the CON solution obtains more dispersed results with the initial values and constraints.Considering that even satellites within the same system may be equipped with different LRAs, the STD derived from the results in the SYS solution cannot effectively serve as a common constraint in the solutions.
Conversely, for the stations with fewer observations, the advantages of the SAT and FIX_SAT solutions become less obvious relatively, but they still outperform other strategies.For instance, for station MONL with the fewest observations these two strategies maintain an advantage of about 2-10 mm compared to other solutions.
Table 7 presents the RMSE values of ERP and geocenter motion with different range bias process strategies.In terms of ERP, the estimation of range bias in SAT solution obtains RMSE values of microarcseconds and 605 microarcseconds for pole coordinates and 151 microseconds for LoD in the 1-day arc.The FIX_SAT solution obtains a smaller Y-pole RMSE, and the similar X-pole and LoD RMSE compared to the SAT solution.The FIX_SYS solution exhibits the smallest RMSE values for Y-pole coordinates and LoD, which are 560 microarcseconds and 125 microseconds, respectively.With the arc length extended to 3-day, the SAT strategy achieves the smallest RMSE values for ERP, which benefits from the threefold increase in observation number and a less significant increase in the errors of parameters.Regarding geocenter motion, the various range bias process strategies have less impact on the results with the differences in RMSE values generally falling within 1 mm.But it is worth noting that the RMSE disparity between our determined geocenter motion and the ILRS official backup combined product (ILRS_b) (Pearlman et al., 2019) series is approximately 7-8 mm in the X and Y components, and 10-11 mm in the Z component.Furthermore, these differences with respect to ILRS_b tend to increase with the extension of the processing arc.This is attributed to the radiation pressure-related characteristics inherent in GNSS orbits, which is presented and explained in Fig. 7.
In addition to the aforementioned methods, range biases can also be grouped based on satellite types and LRA types.We also explore this strategy which is named as GRP, categorizing the three systems into 9 groups (Galileo IOV and FOC; GLONASS M, M + , and K; BDS-2 MEO/IGSO, BDS-3 MEO CAST/SECM, respectively).And the range biases are estimated for each group-station pair within each arc.The post-fit residual RMSE and the slope over nadir angle are both larger than those in the SAT strategy, but smaller than those in other strategies.Comparatively, for the six stations mentioned in Fig. 5 and Table 6, the RMSE of up-component coordinate differences w.r.t.SLRF2014 is 10 mm larger than that in SAT strategy but 4 mm smaller than that in SYS strategy with the 1-day arc length.Regarding ERP parameters, the GRP strategy yielded RMSE differences w.r.t.IERS EOP 14C04 of 559 microarcseconds, 569 microarcseconds, and 179 microseconds for X-pole, Y-pole, and LoD in the 1-day arc.When the arc extended to three days, the ERP differences are 576 microarcseconds, 530 microarcseconds, and 110 microseconds, respectively.In general, the results obtained by the GRP strategy show an improvement compared to the SYS strategy, but the extents are less than the SAT strategy.
Overall, the appropriate correction or estimation of range bias for each satellite in the GNSS SLR process proves effective in addressing satellite-specific biases, leading to improved post-fit residuals and global geodetic parameter results.Furthermore, the estimation of range bias in each satellite-station pair does not adversely impact the results, despite the increase in estimated parameters, given that both the SAT and FIX_SAT strategies yield favorable outcomes.And the SAT strategy accounts better for the time-varying aspects of the signature effect across different arc lengths, resulting in post-fit residuals with smaller RMSE and slope of nadir angle.We employ this strategy in the subsequent section to explore the impact of arc length on the estimation of global geodetic parameters.

Determination of global geodetic parameters
Based on the above analysis, it is evident that estimating range bias in each satellite-station pair yields centralized SLR residuals and optimal results.With the SAT range bias process strategy, the performance for the 1-day, 3-day, 5-day, and 7-day arc length strategies are evaluated from the aspects of ERP, geocenter motion, station coordinates, and scale.

Earth rotation parameters
Figure 6 illustrates the formal errors of X-pole, Y-pole, and LoD parameters under different arc lengths.For the pole coordinates, the 1-day solution yields the average formal errors of 13.49 milliarcseconds and 13.65 milliarcseconds for the X and Y components.As the arc length increases, the formal errors of the pole coordinates for the 3-day, 5-day, and 7-day solutions are decreased by 16.2%, 25.2%, and 31.3%, respectively.The 7-day solution can achieve an average formal error of better than 10 milliarcseconds for the X and Y components.For the LoD, there is also a significant reduction in the formal errors as the arc length extends from 1 to 7-day.Specifically, the formal error of the 7-day solution is 42.9% smaller than that of the 3-day solution, with an average of 0.016 ms.Meanwhile, the strong correlation between the formal error and the number of observations can be observed.Notably, larger formal errors tend to occur on the days with fewer observations.Moreover, extending the arc length to increase the number of observations can yield a considerable reduction in formal errors.Table 8 presents the results of estimated ERP parameters with respect to the IERS EOP 14C04 series from the IERS combination of operational EOP series derived from GNSS, SLR, VLBI and DORIS.For pole coordinates, extending the arc length evidently improves the quality of estimation.In the 1-day arc solution, the RMSEs of X and Y components are 523 microarcseconds and 605 microarcseconds, respectively.The extension of arc length has a significant impact on improving the results.The most accurate pole coordinates are obtained  in the 7-day solution with RMSEs of 473 microarcseconds and 495 microarcseconds for X and Y components, which are improved by 9.6% and 18.2% compared to the 1-day solution.The biases represented by the estimated means are another noteworthy aspect.In X-pole, the minimum bias occurs in the 5-day and 7-day solutions with values of 48 microarcseconds.In the Y-pole direction, the smaller biases are obtained in the 1-day and 3-day solutions with values of − 3 microarcseconds and − 11 microarcseconds, respectively.However, considering the larger RMSEs in the 1-day and 3-day solutions, the 7-day solution with a closer bias of − 15 microarcseconds is recommended.For LoD results, the RMSE for the 1-day solution is 151 microseconds, while an evident improvement by over 50% can be observed in the 3-day solution.Along with the increase of arc length, the accuracy improvement of LoD gradually becomes smaller.
The 5-day solution can almost achieve a similar accuracy as the 7-day solution.The best accuracy of LoD is found in the 7-day solution with an average of 2 microseconds and a RMSE of 52 microseconds.

Geocenter motion
Figure 7 displays geocenter motion series after the Vondrak filter with the cut-off period of 1 day (Vondrak, 1969), using four arc length strategies.The geocenter motion published by ILRS is also presented with the yellow line (Pearlman et al., 2019).It is visible that different arc lengths do not significantly affect the results.
Even for the shortest 1-day solution, we can still observe an evident seasonal variation of geocenter motion.On the other hand, the 1 cpy signal of our solution has a smaller amplitude than that in the ILRS product series, especially in X and Y components.The amplitudes of the 7-day solution in X and Y components are 1.78 mm and 1.71 mm, while the ILRS series has the amplitudes of 3.73 mm and 3.18 mm.In the Z component, the 1 cpy signals in our solutions are similar to that in ILRS product.The amplitude of 1 cpy signal in our 7-day solution is 4.35 mm, while 4.84 mm in the ILRS product.
On the other hand, a significant 3 cpy signal can be observed in our solutions, which has an amplitude of 9.8 mm in the 7-day solution.However, the 3 cpy signal cannot be observed in the ILRS product.This is related to the periodic variation of Sun elevation angles above and below the orbital planes of GNSS satellites, which affects the solar radiation pressure on satellites (Meindl et al., 2013).Figure 8 shows the Z component of geocenter motion from the 7-day solution, ILRS product, and Center for Orbit Determination in Europe (CODE) series Meanwhile, the series can be affected by solar radiation pressure due to the variation of Sun elevation angles.

Station coordinates and scale
Figure 9 shows the differences of the up components of the station coordinates w.r.t.SLRF2014 under different arc length strategies.The PMT stations exhibit larger whiskers and IQRs compared to C-SPAD and MCP stations, primarily due to their relatively smaller number of observations.Different from the geocenter motion, the extension of arc length leads to a significant improvement in station coordinates.In the 3-day solution, the whiskers are about 10 mm smaller than those in the 1-day solution, and the outliers are also significantly decreased.The 5-day solution shows an improvement with respect to the 3-day solution, but it is less noticeable.In the 7-day solution, the station coordinates are almost the same as in the 5-day solution with only slight reductions of whiskers and IQR.
The scale results in this study are obtained by Helmert transformation.Figure 10 shows the scale results with respect to SLRF2014 and the number of stations used in Helmert transformation in different arc lengths.Extending the arc length significantly increases the number of stations involved in the Helmert transformation.The station numbers of the 3-day, 5-day, and 7-day solutions are increased by 33.4%, 44.6%, and 46.8% compared to the 1-day solution.The 7-day solution exhibits the minimum RMSE of scale at 2.10 × 10 -9 , but the improvement over the 5-day solution is only 0.09 × 10 -9 .Additionally, biases of approximately -0.6 × 10 -9 can be observed in the four series.This is attributed to the fact that the scale of ITRF is determined by both SLR and VLBI techniques.Therefore, the scale determined by SLR-only will show biases with respect to that provided by ITRF.Taking a comprehensive view, the 5-day solution is deemed adequate to fulfill the scale estimation requirement.

Conclusions
The SLR technique plays an irreplaceable role in the realization of the ITRF.However, only the SLR observations to LAGEOS and Etalon satellites are used in the current realization of ITRF.With the increasing of GNSS SLR observing activities, the total number of SLR observations to GNSS satellites exceeds that to LAGEOS.These data hold the potential to contribute to the estimation of global geodetic parameters.This study shows a detailed investigation on range biases and arc length used in GNSS SLR data processing strategy.This paper aims to identify a processing strategy better suited for the current abundance of GNSS SLR observations, which is of significance in incorporating GNSS SLR observations into ITRF realization and enhancing the accuracy of estimated global geodetic parameter.
By investigating five strategies for range bias, we demonstrate that the estimated range bias in SAT solution (estimated in each satellite-station pair) can yield more residuals with smaller RMSE and the slope closer to zero than other solutions, which is attributed to the better absorption of signature effects.Additionally, the estimated range bias STD in SAT solution are smaller than those in SYS solution (estimated in each system -station pair) for 72.2% GLONASS, 97.0%Galileo, and 74.4% BDS satellites.Then we prove that the SAT and FIX_SAT solutions obtain the coordinates that are more stable and closer to SLRF2014.From the perspective of ERP, the RMSEs of pole coordinates reach 495 microarcseconds and 522 microarcseconds for X and Y components, while the LoD reaches 73 microseconds in the 3-day arc and SAT solution, which are smaller than those in other solutions.For geocenter motion, different range bias processing strategies exhibit minor impact.The RMSE differences among the five strategies are generally within 1 mm.Overall, our results suggest that the SAT solution is a more effective method for estimating the range bias in GNSS SLR processing.
Then we analyze the effect of arc length on the determination of earth rotation parameters, geocenter motion, and scale using the three-year SLR observations to Galileo, GLONASS, and BDS satellites.The results of the estimated ERP parameters can be greatly improved if the arc length is extended.The results with the 7-day solution are improved by 9.6% and 18.2% compared to the 1-day solution with RMSEs of 473 microarcseconds The optimal LoD result is also obtained with the 7-day solution, with an average of 2 microseconds and a RMSE of 52 microseconds.However, the geocenter motion results are not significantly affected by the arc length, and the seasonal variation can also be observed even in the shortest 1-day arc solution.On the other hand, a 3 cpy signal in the Z component can be observed in GNSS solutions, including our results and CODE series, which is caused by the periodic change of solar radiation pressure due to the variation of the Sun's elevation angle.Different from the geocenter motion, the extension of arc length leads to a significant increase in the number of stations used in the process and the reduction of station coordinate errors w.r.t.SLRF2014 in short arcs, but these improvements become less noticeable with the arc length extended.The scale obtained by Helmert transformation shows a similar trend.The minimum RMSE of scale are exhibited in the 7-day solution, which is 2.10 × 10 -9 with respect to ITRF2014, but the improvement compared to the 5-day solution is only 9 × 10 -11 .
With the increasing number of GNSS satellites tracked by ILRS, the determination of global geodetic parameters will benefit significantly.For instance, ILRS has tracked more BDS-3 MEO satellites since February 2023, bringing the current number of tracked BDS satellites to 28 (https:// ilrs.gsfc.nasa.gov/ missi ons/ satel lite_ missi ons/ curre nt_ missi ons/ index.html).These new observations hold huge potential for enhancing the determination of global geodetic parameters.Additionally, the new SLRF2020 frame employed the Station Systematic Error Monitoring (SSEM) approach (Luceri et al., 2019), which can evidently mitigate the systematic errors.These progressive changes are anticipated to drive further advancements in GNSS SLR processing solutions, contributing to the accurate determination of global geodetic parameters.

Fig. 2
Fig. 2 Daily SLR observation number of each GNSS satellite between 2019 and 2021

Fig. 3
Fig. 3 Average and STD values of range bias estimation results with 1-day solution in 6 SLR stations.The blue and red lines are the average and STD value reference lines in the SYS solution

Fig. 4
Fig. 4 SLR residuals of Galileo satellites and YARL between 2019 and 2021 with the five different range bias strategies Fig. 5The box plots of up component errors for 6 stations in 2019-2021, with different range bias process strategies and the 1-day/3-day arc length strategies

Fig. 6
Fig. 6The formal error of ERP estimated in different arc length (up), and the observation number of process in each week (bottom)

Fig. 7
Fig. 7 Geocenter motion in three components with different arc length strategies in time domain (left) and frequency domain (right), the yellow lines represent the result of ILRS product

Fig. 8
Fig.8Geocenter motion from the 7-day solution, CODE series, and ILRS series, the lines represent the result smoothed with Vondrak filter

Fig. 9
Fig. 9 The box plots of station residuals w.r.t.SLRF2014 in up component with different arc length strategies.The three stations in each type are the stations with the most observations

Fig. 10
Fig. 10 Scale differences series w.r.t.ITRF2014 in different arc length (up), and the station number used in Helmert transformation each day with different arc length strategies (bottom)

Table 1
Monument, code, type, and location of 29 ILRS stations with SLR to GNSS used in this study of SLR observations compared to the other two systems.On average, each GLONASS satellite has more than 12 observations per day, while each Galileo and BDS satellite have had an average of 7 and 6 observations per day, respectively.Moreover, it is worth noting that the number of observations tends to be relatively smaller from June to August each year, which is mainly due to the station distribution.During the summer in the northern hemisphere, SLR stations experience the influence of longer daytime periods and increased

Table 2
The range bias process strategies of SLR to GNSS SYSEstimated for system-station pair in each arc CON Estimated with prior value and constraint in each arc FIX_SYS Fixed as an annual bias for each system SAT Estimated for satellite-station pair in each arc FIX_SAT Fixed as an annual bias for each satellite-station pair observations exhibited larger STD values o the estimated range biases.Notably, across all six stations, the STD values of the range biases obtained by the SAT solution are mostly smaller than those derived from the SYS solution.This trend persists for the stations with a relatively lower observation count, such as station MONL.

Table 4
The process strategies for the prior GNSS satellite orbits

Table 5
Average RMSE of orbit differences (2019-2021)in along-, cross-track, and radial directions and 3D between the orbits and the products ofCOD and WUM (cm) Satellite type Differences w.r.t.COD products Differences w.r.t.WUM products

Table 6
The RMSE of ENU component errors (mm) w.r.t.SLRF2014 for 6 stations in 2019-2021, with different range bias process strategies and the 1-day/3-day arc length strategies

Table 7
The RMSE of ERP w.r.t.IERS EOP 14C04 in different range bias process strategies (left), and the RMSE of geocenter motion w.r.t.ILRS_b in different range bias process strategies (right)

Table 8
The ERP result w.r.t.IERS EOP 14C04 in different arc length strategies