Impacts of inter-satellite links on the ECOM model performance for BDS-3 MEO satellites

Inter-satellite link (ISL) plays an essential role in current and future Global Navigation Satellite System (GNSS). In this study, we investigate the impact of ISL observations on precise orbit determination for BeiDou-3 Navigation Satellite System (BDS-3) Medium Earth Orbit (MEO) satellites based on different Extended CODE Orbit Models (ECOM). Thanks to the better observation geometry of the Ka-band ISL data compared to the L-band data for BDS-3 MEO satellites, the ISL solution substantially reduces Orbit Boundary Discontinuity (OBD) errors, except for C30, which suffers from unstable Ka-band hardware delay. From the external quality analysis, ISL significantly enhances the reliability of the orbit of MEO satellites manufactured by the China Academy of Space Technology (CAST). The standard deviation (STD) of the satellite laser ranging (SLR) residuals is approximately 2.5 cm, and the root mean square (RMS) is reduced by 10–23% compared to L-band solutions. Besides, the Sun-elongation angle dependent systematic error in SLR residuals nearly vanishes based on the reduced 5-parameter ECOM (ECOM1) or extended 7-parameter ECOM (ECOM2) with ISL data. This is because the ISL reduces the correlation between state parameters and solar radiation pressure (SRP) parameters as well as those among SRP parameters, leading to a more accurate estimation of both orbit and SRP perturbations, particularly those along B direction. This confirms that the deficiency of the SRP models for BDS-3 CAST satellites can be compensated by using better observation geometry from ISL data. On the other hand, for the satellite manufactured by Shanghai Engineering Center for Microsatellites (SECM), the ISL allows for a more accurate estimation of the Bc 1 parameter in the ECOM1 model. This only reduces linear systematic error, possibly because the impact generated by the satellite bus cannot be entirely absorbed by the B -direction parameters.


Introduction
The BeiDou navigation satellite system (BDS) has been fully operational, with the global constellation (BDS-3) providing Positioning, Navigation, and Timing (PNT), short message, international search and rescue, as well as other services since July 31, 2020.Currently, the BDS-3 constellation consists of 30 satellites, i.e., 3 GEO satellites, 3 IGSO satellites, and 24 MEO satellites.Additionally, the BDS-3 satellites are equipped with Inter-satellite Link (ISL) devices for communication and ranging in Ka-band.ISL is an innovative key technology, which is playing a revolutionary role in current and future GNSS (Glaser et al., 2020).Many studies explored the contribution of ISL to orbital accuracy.With ISL observations of 8 satellites, the 3D overlapping orbit differences can be reduced from approximately 1.0 m to 0.5 m based on 10 regional stations (Yang et al., 2020).Xie et al. (2019) showed a 42% or 83% improvement in orbital accuracy was achieved by combining ISL observations with 16 global or 6 regional stations, respectively.Moreover, Michalak et al. (2021) performed the full-scale simulations for the future GNSS constellation "Kepler", and the results showed that ISL can significantly improve the performance of the determined orbit, clock, and geodetic parameters.
On the other hand, precise modeling of the various perturbations acting on GNSS spacecrafts is essential, particularly the solar radiation pressure (SRP), as it is the most significant non-conservative force acting on GNSS satellites.For the SRP modeling of BDS satellites, extensive research has been conducted, with three primary approaches: empirical, analytical, and semi-analytical models.The Empirical CODE orbit model (ECOM), first introduced by Beutler et al. (1994), has been widely used for BDS satellites SRP modeling, with its simplified 5-parameter model (ECOM1).However, it has been observed that ECOM1 is not applicable for the satellites in orbit-normal (ON) attitude mode.Guo et al. (2013) identified the clear deficiencies of ECOM1 during ON mode for BDS-2 IGSO and MEO satellites.Moreover, for BDS GEO satellites, there is a strong correlation between the orbital and ECOM1 parameters due to the use of ON attitude mode and poor observation geometry.To overcome these shortcomings, Guo et al. (2016) added an empirical force parameter in the along-track direction.For BDS-3 MEO satellites, the Sun-elongation angle ( εangle, the angle between the Earth and Sun direction as observed from the satellite) dependent linear error can be clearly identified within the orbits determined with ECOM1 model (Zhao et al., 2022).Yan et al. (2019) used the extended ECOM (ECOM2) model to reduce this error.Furthermore, Duan et al. (2022) divided the 18 BDS-3 MEO satellites into four groups and established a hybrid SRP model to reduce the systematic error.
The above-mentioned research for SRP models are based on L-band data.Orbit determination relies not only on dynamic models, but also on observation geometry.Compared with L-band data, ISL has a better observation geometry.Zhao et al. (2022) demonstrated that ISL measurements can effectively reduce Sun-elongation angle dependent systematic error in the BDS-3 MEO orbits determined with ECOM1 without a prior model.This discovery served as a motivation for an in-depth investigation of the impacts of ISL on BDS-3 MEO SRP models and the reasons.We begin with an introduction of ECOM series models and the general structure of the BDS-3 MEO satellites, followed by an analysis and discussion of the influence of ISL on the BDS-3 MEO orbit determined with ECOM1 and ECOM2.Finally, the conclusion is drawn.

ECOM models
The ECOM model is widely employed in the GNSS community.To account for SRP perturbations, Beutler et al. (1994) and Springer et al. (1999) proposed a Sun-oriented reference frame, which comprises three unit vectors: where u and u s are the latitude of satellite and Sun in the satellite's orbital plane, respectively.The u can be replaced by the orbital angle µ (Arnold et al., 2015;Mon- tenbruck et al., 2015), which is adopted in this study.

BDS-3 MEO satellites and ISL
The BDS-3 satellites, manufactured by China Academy of Space Technology (CAST) and Shanghai Engineering Center for Microsatellites (SECM), are distributed on three different orbital planes.Satellites C27-C30, C34, C35, C43 and C44 are in slot-A; satellites C19-C22, C32, C33, C41 and C42 are in slot-B; and satellites C23-C26, C36, C37, C45 and C46 are in slot-C.The metadata for the BDS-3 satellite was released by China Satellite Navigation Office (CSNO) at the end of 2019.Figure 1 depicts an artist's impression of the BDS-3 MEO satellites from CSNO and SECM.From Fig. 1, we can observe that CAST MEO satellites are based on the dedicated T-shaped platform with a panel for providing search and rescue (SAR).SECM MEO satellites are manufactured based on the same bus (SECM-A) with the exception of C43 and C44 (SECM-B), which exhibit slightly different shapes.Moreover, we can see from Fig. 1 that the satellite bus is non-cubic for both CAST and SECM missions.These two types of satellites are stretched in different (1) directions.And CAST MEO satellites with a mass of approximately 941-1061 kg are lighter than SECM MEO satellites with a mass of approximately 1010-1079 kg.
Additionally, BDS-3 adopts the Concurrent Spatial Time Division (CSTD) scheme to achieve inter-satellite communication and ranging.Under this design, a dual one-way measurement is completed by a pair of satellites within 3 s based on predefined timeslot schedule.Due to the Earth obstruction and the maximum nadir angle limitation (60°), each satellite is unable to establish a connection with its adjacent satellite and the farthest satellite in the same orbital plane.Instead, it can form continuous links with the second and third nearest satellites within the same orbital plane, resulting in a total of four continuous links in each plane.In terms of out-of-plane links, there are twelve discontinuous links and four continuous links.Furthermore, each IGSO and GEO satellite has the capability to establish connections with the MEO satellites.And the IGSO satellites can establish connections with one another.Zhao et al. (2022) provide a detailed description of the connection scheme.Figure 2 shows the observation geometry on day of year (DOY) 336, 2020, taking C20 as example.It can be clearly observed that the observation geometry of ISL is significantly better than that of L-band, for which the range of the nadir angle is approximately 0°~14°.The nadir angle for in-plane links is approximately 23° and 45°, whereas the nadir angle of out-plane links varies from 16° to 60°.

Strategy
Since 2020, there have been over 100 monitoring stations with BDS-3 MEO satellites tracking capability.In this study, the data of approximately 1 year (DOY 007, 2020 ~ DOY 337, 2020) are selected for processing.In our analysis, we select approximately 15 International GNSS Monitoring and Assessment System (iGMAS) stations and 100 IGS stations for precise orbit determination (POD) based on the Position and Navigation Data Analyst (PANDA) software (Liu & Ge, 2003).Figure 3 shows the distribution of these stations.
A two-step approach is applied for POD.We first estimate station coordinates, receiver clock offsets, and tropospheric zenith delays by utilizing the Precision Point Positioning (PPP) method based on GPS L1 and L2 observations and IGS precise orbits as well as 30 s clocks.Then, these parameters are fixed, and BDS-3 B1I as well as B3I observations with a sampling rate of 300 s are used to estimate orbits, satellite clocks, SRP parameters, and ambiguity parameters in a 24-h batch mode.The POD strategy is nearly identical to the one used by Guo et al. (2023), with the exception of SRP models.In the ISL data processing, the original dual one-way observation is first transformed into a single observable for orbit determination.The transformation algorithm is presented by Tang et al. (2018).The reference time t 0 used for the reduction of dual one-way ISL data in this study is the middle epoch  where P BA (t 0 ) and P AB (t 0 ) are the transformed forward and backward observations, respectively.R A (t 0 ) and R B (t 0 ) are the position vectors of satellites A and B at t 0 , respectively.c is the speed of light in vacuum.δ A and (2) δ B are the hardware delays of satellites A and B, respec- tively.In our processing, the hardware delay is estimated as a constant parameter for each arc.ε is the noise of the ISL measurement.The detailed strategy can be found in Table 1.Since all ISLs between satellites are used, there is no need to set a cutoff angle.Additionally, the troposphere delay is not considered, as the ISLs between satellites do not pass through the troposphere.In the GNSS community, it is widely acknowledged that the precisions of phase and code are in the millimeter level (about 3 mm) and decimeter level (about 30 cm for precise code).Considering the accuracy of the ISL is approximately 3 cm, the weight ratio was set as 10,000:100:1 for the phase, ISL, and code measurements in this study (Yang et al., 2020).
To assess the impact of ISL on the orbital characteristics of BDS-3 MEOs, we conducted experiments that differ in terms of (1) the usage of data and (2) the SRP models (as outlined in Table 2).Three primary strategies were employed in this study.Strategy "E1" refers to the ECOM1 model.Strategy "E2" represents the ECOM2 model with 7 parameters.And strategy "EA" incorporates an additional along-track constant perturbation based on the ECOM1 model.Moreover, the last letter "I" in the solution denotes the integration of ISL data and L-band data for joint POD.

Results and analysis
This section analyzes the results of the BDS-3 POD based on L-band and ISL data.We evaluate the accuracy of the orbit solution and explore the impact of the introduction of ISL observations for POD.

Orbit boundary discontinuity (OBD)
The OBD describes the difference between the satellite orbits from two adjacent orbit arcs at the same daily boundary epoch, and it indicates the internal quality of orbit solutions.Figure 4 illustrates the OBD of BDS-3 MEO satellites for the six solutions.Here, we selected C42, C44, and C46 to represent the orbit performance for the satellites with PRN beyond 40.Significantly different performances can be seen among the different satellites.For those solutions with L-band data only (E1, E2 and EA solutions), the same phenomenon can be observed.Specifically, satellites C19-C25 and C32-C37 demonstrate the highest level of consistency, followed by C26-C30.However, orbits of C42-C46 show the worst performance, primarily due to the limited tracking data available in the selected study period.Once the ISL observations are used, the median decreases by approximately 5-30 cm.Among the solutions, E2I shows the best performance, and the median of the OBD is 6.6 cm with the inter-quartile range (IQR) at the level of 4.7 cm.Meanwhile, the Q3 + 1.5•IQR (top whisker) decreases to below 25 cm, except for C30 and C44.The improvement can be due to the better observation geometry.
To seek the reason for the degradation of orbit quality by using ISL for C30, Fig. 5 further illustrates the OBD errors in the along-track, cross-track and radial direction.For L-band solutions, i.e., EA, E1 and E2, they show almost similar results in cross-track and radial direction, while E2 has a better performance in the along-track direction.Furthermore, there is no significant Sun-elevation angle ( β-angle, the angle between the Sun direc- tion and the satellite orbital plane) dependent systematic error and the degradation of orbital accuracy during the eclipse periods.Once the ISL data are used, the OBDs are significantly reduced during DOY 147-340, 2020.However, before DOY 147, 2020, degenerated performance can be observed.To illustrate the reason, we redetermined the orbit by modeling the hardware delay as random parameters instead of the constants, and the estimated values are plotted in Fig. 6.The top two subplots in Fig. 6 present the estimated constant hardware delay for C20 and C30 during the study period, while the bottom two display the epoch-wise estimates for C20 and C30 during DOY 146 to 147, 2020.In general, the hardware delay of C20 is quite stable with an STD of about 0.13 m.However, unstable changes before DOY 147, 2020 and a   5, the stable hardware delay shows a strong correlation with the better OBD.This confirms that the stability of the Ka-band hardware delay has a significant impact on the performance of the ISL solutions.However, the cause of changes in hardware delay before DOY 147, 2020, requires further investigation.

SLR validation
SLR analysis is an effective and independent validation method for radial orbit.Every BDS-3 spacecraft is equipped with a laser-reflected array, and the International Laser Ranging Service (ILRS) has tracked all BDS-3 MEO satellites since 2023.However, during the study period, only C20 and C21 manufactured by CAST as well as C29 and C30 produced by SECM were tracked by ILRS.The LRA offsets with respect to the Center of Mass (CoM) in the body-fixed frame disclosed by CSNO (2019) are used in this study for validation.Table 4 presents the statistical results of SLR residuals for the six solutions.The residuals larger than 0.3 m were identified as outliers and removed.For C30, the results from the periods with unstable hardware delays (before DOY 147) are excluded here.As a result, for C20, C21, C29, and C30, there are a total of 3457, 3604, 2651, and 1087 normal points out of 3458, 3611, 2657, and 1089, respectively.
In general, the four satellites exhibit different performances.However, the satellites from the same manufacture show similar performance.Specifically, CAST satellites have positive mean value up to around 6.5-7.4 cm, whereas it is about − 1.2 ~ − 3.3 cm for SECM satellites.For L-band solutions, E1 and EA demonstrate similar performance, while E2 has the best performance.
Once the ISL data are used, the mean of SLR residuals is reduced from around 7.0 cm to 5.9 cm for C20 and C21, and slightly improved from − 2.3 cm to − 2.5 cm for C29 and C30.Meanwhile, the STD is reduced from around 3.6 cm to 2.5 cm for C20 and C21, and 4.1 cm to 3.5 cm for C29 and C30.It clearly confirms that ISL improves the orbit quality.
More importantly, SLR can be used for detecting systematic errors in mismodeling orbit perturbations.The SLR residuals of the four MEO satellites are significantly larger in the range of |β|< 20°.Once the ISL data are added, the systematic ε-angle dependence of C20 and C21 almost disappeared, i.e., the slopes of the SLR residuals of C20 and C21 are decreased from − 0.073 cm/degree and − 0.078 cm/degree for the solution E1 to 0.006 cm/degree and − 0.010 cm/degree for the solution E1I, respectively.While a slight reduction is obtained for C29 and C30, a noticeable error can still be observed in the SLR residuals for SECM satellites.The reasons will be explained below.
Additionally, the ECOM2 was reported to effectively reduce the Sun-elongation angle dependent systematic error (Li et al., 2020;Yan et al., 2019).Hence, Fig. 8 depicts the SLR residuals with respect to ε-angles for the four BDS-3 MEO satellites determined with ECOM2 model.It is clearly observed that the linear error is reduced compared to the solution E1 with slopes of − 0.013 cm/degree, − 0.017 cm/degree, 0.034 cm/degree and 0.014 cm/degree for C20, C21, C29 and C30, respectively.By using ISL data, the slopes for C20, C29 and C30 are further reduced to − 0.005 cm/degree, − 0.015 cm/ degree and − 0.009 cm/degree, respectively, while almost no change is observed for C21.

ECOM parameters
Figure 9 illustrates the estimated ECOM1 parameters with respect to β-angles for C20 (CAST) and C29 (SECM) of the solutions E1 and E1I.As to the D 0 parameter of C20, its absolute value exhibits a positive correlation with the |β| beyond 25°, while it varies slightly for |β|< 25°.For C29, the absolute value of D 0 shows a negative correlation with |β| during the noneclipse season, while it has a positive correlation during the eclipse season.After including the ISL data, the estimates of D 0 become smoother for both satellites, particularly in high |β| regimes.Regarding the Bc 1 , the two satellites show distinctly different variations.The direction of Bc 1 variations for C20 is reversed when |β|> 60°, while it remains unchanged for C29, mainly due to the maximum |β| of 34°.Moreover, a noticeable bias can be observed for both.For C20, it is up to approximately 6.7 nm/s 2 , whereas this deviation for C29 is relatively smaller, ranging from only about 0 to 4.6 nm/s 2 .For the B 0 and Bs 1 , the inclusion of ISL data makes the estimates more stable.Figure 10 presents the estimated ECOM2 parameters of the solutions E2 and E2I for C20 and C29.The variations of D 0 are similar to those in Fig. 9, however, the estimates of D 0 become different as two parameters are added, i.e., Dc 2 and Ds 2 .With ISL data, the estimates become more precise and concentrated, particularly for C29 in the eclipse season.For the parameters in the B direction, the variations of the three parameters, i.e., B 0 , Bc 1 , and Bs 1 , for C20 are entirely different from those of C29.The patterns as well as signs of those parameters show clear opposite variations for both satellites, and this is possibly caused by the different orientation of the stretched surfaces of satellite buses.Once the ISL data are used, the estimates of these parameters become more precise, as shown in Fig. 10.However, different from the ECOM1 solution, no noticeable bias can be observed in Bc 1 and Bs 1 estimates.

Correlations
To explain the effect of ISL on parameter estimation, we investigated the correlation between the state and estimated ECOM1 parameters of the solution E1 and E1I for C20 at different β-angles.As shown in Fig. 11, ISL data markedly reduces the correlation between ECOM1 parameters and state parameters.This enables a more precise calculation of the satellite orbit and the perturbations.Moreover, for low to medium β-angles, the correlation between SRP parameters is reduced, resulting in a more accurate estimation of the parameters in the DYB direction.This has a significant impact on the estimates of Bc 1 parameter, resulting in the improvement of its accuracy and the noticeable bias.Similarly, Fig. 12 shows the correlation among the state and ECOM2 parameters of the solution E2 and E2I for C29.It is clearly observed that ISL results in a reduction of the correlation among most of the parameters at β = 0°.This confirms that the inconsist- ency between the estimated parameters of E2I and E2 is due to more precise estimation of SRP parameters from ISL during the eclipse period.Furthermore, the strong correlation between Bc 1 and Dc 2 indicates that the inclusion of Dc 2 parameter leads to a more accurate estimation of the Bc 1 parameter, not generating the noticeable biases in the Bc 1 estimates for both E2I and E2 solutions.

Accelerations
To investigate the impacts of ISL on the estimated SRP perturbations, Fig. 13 presents the differences of SRP accelerations between E1 and E1I solutions for C20 in the along-track, cross-track, radial, and D, Y, B directions.In general, the differences in the along-track and cross-track direction are larger than those in the radial direction.Clear antisymmetric variations on µ-angles or β-angles can be observed for the acceleration differences in the along-track and cross-track, respectively.However, the radial acceleration shows symmetrical variation on both β and µ-angles.
Considering that SLR validation mainly indicates the radial accuracy, we focus on the radial acceleration in this study.In addition, it can be clearly observed from Fig. 13 that the effect of ISL on SRP is mainly concentrated in B direction.The differences in the D direction are only up to 1.5 nm/s 2 in high β regimes, while those in the B direction can reach 6 nm/s 2 at midnight or noon points.This is consistent with the estimates of ECOM1 shown in Fig. 9.An important finding is that the acceleration differences in the B direction present an almost negative pattern compared to the radial acceleration differences.This is because the angle between B and radial direction is always greater than 90°.Meanwhile, the same variation trend confirms that the radial orbit differences mainly come from the perturbations in B direction.On the other hand, the ECOM2 model can effectively reduce the systematic Sun-elongation angle dependent linear trend in the SLR residuals.Hence, we investigated the differences of SRP perturbations of the solutions E1 and E2 within the DYB frame, as shown in Fig. 14.As expected, the large differences along the D direction are observed, as two additional periodic parameters are estimated in D.Moreover, noticeable differences in the B direction can also be observed, and the variation in the B direction is very similar to that shown in Fig. 13(f ).From the above correlation analysis, it can be observed that there is a significant correlation between Bc 1 and Dc 2 , as well as between Bs 1 and Ds 2 .The inclusion of Dc 2 and Ds 2 plays a key role in absorbing the effects of Bc 1 and Bs 1 , leading to more accurate estimation.This, in turn, contributes to the differences in acceleration of the B direction.For CAST satellites, ISL compensates for the deficiencies of the ECOM1 model by directly enhancing the accuracy of the SRP estimation in the B direction.
Considering the radial orbit are highly correlated with the clock, the linear clock-fit (LCF) residuals can be used to demonstrate the radial orbit errors.In Fig. 15, the differences of the LCF residuals between the solutions E1 and E1I w.r.t.β-angles and μ-angles are shown for the selected representative satellites, i.e., C20, C21, C29 and C30.As reported by the TARC of CSNO, the Rubidium Atomic Frequency Standard (RAFS) clocks are used as the primary onboard clock for C20 and C21, whereas According to the theory of the analytical SRP model (Milani et al., 1987;Montenbruck et al., 2015;Rodriguez-Solano et al., 2012), SRP results from SPs and satellite surfaces illuminated by the Sun.Only the + X, ± Z surfaces and SPs are illuminated in the nominal yaw steering mode for BDS-3 MEO satellites.The perturbation generated by the SPs primarily contributes to the D direction, while the SRP in the B direction mainly originates from the satellite bus.As shown in Fig. 16, the projection of SRP generated by the CAST satellite bus in the D direction is approximately 10 to 25% of the SPs, whereas that of SECM ranges from 13 to 45%.The variation in the projection of SRP generated by the CAST satellite bus in the B direction is larger than that of the SECM, as also illustrated in Fig. 9.These explain why improving the B direction estimation of the ECOM1 model with ISL causes the linear ε-angle depend- ence to almost disappear for the CAST satellite, whereas it is only reduced for the SECM satellite.

Discussion
There are a few issues to be further addressed.For previous analysis, the weight ratio was set as 10,000:100:1 for the phase, ISL, and code measurements.However, we also observed that the different weighting ratio settings for the phase, code, and ISL data can affect the results.For instance, when the weight ratio of the ISL is set to 10 or 1000 instead of 100, the sun-elongation angle dependent systematic errors are not eliminated for all the four satellites of these two solutions, as shown in Fig. 17.For the two SECM satellites, i.e., C29 and C30, the slope of linear variation is increased from around 0.058 cm/degree to 0.094 cm/degree.While the noticeable decrease of linear slope can be observed compared to the solution E1 for both CAST satellites, i.e., C20 and C21.The sign of the slope even changes.This demonstrates the ISL has more pronounced impacts on the orbit determined based on ECOM1 for CAST satellites than that of SECM satellites.The reason for this discrepancy may be because the SECM satellite body structure contributes strongly to the acceleration in the D direction and is less susceptible to receiving changes in the B direction.This also clearly indicates the significant impact of the stochastic model on the orbit solution.Therefore, it is essential to use the appropriate stochastic model.According to the analysis by Furthermore, as to the potential effects from space weather conditions, the ionospheric delay and tropospheric delay do not have influences on the ISL measurements between the satellites used in this study for data processing.Regarding the effects of other space weather, such as solar wind, geomagnetic storm, etc., on the orbit are absorbed and reflected by the SRP model.However, the deficiencies in the ECOM1 for BDS-3 are not mainly from the mismodelling of the space weather' effect, as shown by the results, and ISL data or ECOM2 can be used to reduce the deficiencies of ECOM1.For the timedependent variations in the SLR residuals of the orbit determined with improved SRP models, these are possibly related to the space weather.From the perspective of orbit accuracy, it appears that ECOM2 demonstrates good performance in POD for BDS-3 MEO satellites without the a priori model.Hence, the ECOM2 model is recommended for BDS-3 MEO POD if no a prior SRP model is used.

Conclusion
The introduction of ISL data for GNSS satellites offers a new perspective for studying the characteristics of satellites orbits.In this study, we investigated the impact of the ECOM1 and ECOM2 models on the BDS-3 MEO POD.We explored the significant effect of ISL data on the orbital characteristics of the BDS-3 MEO satellites and provided a physical explanation for this phenomenon.
The internal quality of solutions based on L-band data (E1, EA and E2) is at a similar level.By using the ISL observations, the median of 3D OBD is decreased by approximately 5-30 cm compared to the L-band solutions.However, the solutions E1I and EAI are characterized with significantly lower internal stability to the C30 than the solutions E1 and EA.We have confirmed that this issue is attributed to the unstable Kaband hardware delay of the C30.This finding emphasizes the substantial impact of Ka-band hardware delay stability on the performance of joint ISL solutions.Additionally, it is worth noting that the unstable hardware delay for C30 occurs before DOY 147, 2020, when the BDS-3 constellation was under the in-orbit validation phase.
Once the system was commissioned on July 31, 2020, it remains stable as shown by the results.
The external quality analysis shows that SLR residuals of the CAST satellites are characterized by a positive bias whereas the SECM satellites are characterized by a negative bias.This is consistent with the results presented by the GFZ analysis center with the ECOM1 model (Steigenberger et al., 2023).ISL clearly makes the CAST satellites orbit results more reliable, i.e., the STD of the SLR residuals is about 2.5 cm and the RMS is reduced by 10-23% compared to L-band solutions.This result is approximately 29% better than the POD solutions of Duan et al. (2022).The benefits of the ISL data are visible when considering the SLR residuals in relation to the ε-angles.For the solution E1I, the systematic ε-angle dependence of C20 and C21 almost disappeared and is decreased for C29 and C30.From our analysis of the SRP accelerations in both the L-band and ISL solutions, we observe that the acceleration differences in the B direction exhibit a negative pattern compared to those in the radial direction.Additionally, the variation introduced in the B direction by ECOM2 is similar to that introduced by the ISL.This is due to the strong correlation between Bc 1 and Dc 2 , as well as between Bs 1 and Ds 2 .The ECOM2 model increases the accuracy of the estimates in the B direction by adding two periodic parameters in the D direction.Furthermore, the ISL reduces the correlation between state parameters and SRP parameters as well as among SRP parameters, leading to a more accurate estimation of both orbit and SRP perturbations.Nonetheless, for the SECM satellite, characterized by small SPs and larger main illuminated surfaces, the ISL can only reduce the linear systematic error.This is because the satellite body structure contributes significantly to the acceleration in the D direction, and this contribution cannot be fully absorbed by the parameters in the B direction.Generally, the deficiency of the SRP model can be compensated by using better observation geometry from ISL data.
This study provides a new perspective to view the SRP perturbation.Considering that future GNSS constellations will consist of LEO satellites and optical ISLs will be used, more accurate dynamic models and contribution from better geometry can be expected.The estimated accelerations and orbits are available upon request.The ISL observation data are available from the corresponding author upon reasonable request and with the permission of the China Satellite Navigation Office (CSNO).
to form the right-hand reference system.To reduce the draconitic signals of the geodetic parameters derived from GNSS, Arnold et al. (2015) summarized and extended the ECOM model.In this model, the accelerations are expressed as: The ECOM1 model can be derived by taking n D = 0 and n B = 1.With n D = 1 and n B = 1, we get ECOM2 model.The constant perturbation along the ⇀ e D direction is repre- sented by D 0 , while Y 0 and B 0 denote the biases in the ⇀ e Y and ⇀ e B direction, respectively.D 2i,c and D 2i,s denote the even-order perturbations along the ⇀ e D direction.B 2i−1,c and B 2i−1,s represent the odd-order perturbations along the ⇀ e B direction.The angular parameter u ≈ u − u s ,

Fig. 4
Fig. 4 Boxplots of 3D OBD errors for the E1 (red), E1I (purple), EA (green), EAI (blue), E2 (yellow), and E2I (cyan).The box is divided into three parts: the upper boundary (Q3 quantile), the middle groove (median) and the lower boundary (Q1 quantile).The length of the entire box is inter-quartile range (IQR).The upper and lower short lines of the box extension represent Q3 + 1.5•IQR and Q1 − 1.5•IQR, respectively Figure  7illustrates the SLR residuals with respect to the Sun elongation angle for the solutions E1 and E1I.As expected, SLR residuals clearly show a linear systematic error with respect to the ε-angles, and the fitting slopes are similar (around 0.074 cm/degree) with negative values for CAST and positive values for SECM satellites.

Fig. 5
Fig. 5 Time series of OBD in the along-track (red), cross-track (green) and radial (blue) direction for C30 of the six solutions, i.e., E1, E1I, EA, EAI, E2, and E2I.The shaded area represents the eclipse period

Fig. 6
Fig. 6 Variations of Ka-band hardware delays of C30 (left) and C20 (right) for E1I solution.The top two subplots present the estimated constant hardware delay for C20 and C30 during the study period, while the bottom two display the epoch-wise estimates for C20 and C30 during DOY 146 to 147, 2020

Fig. 7
Fig. 7 SLR residuals with respect to ε-angles for the BDS-3 CAST C20 and C21 as well as SECM C29 and C30 satellites of E1 (a-d) solution and E1I (e-h) solution.The results from the periods with unstable hardware delays (before DOY 147) are excluded for C30

Fig. 8 Fig. 9
Fig. 8 SLR residuals with respect to ε-angles for the BDS-3 CAST C20 and C21 as well as SECM C29 and C30 satellites of E2 (a-d) solution and E2I (e-h) solution.The results from periods with unstable hardware delays (before DOY 147) are excluded for C30

Fig. 13
Fig. 13 SRP acceleration differences between the solution E1 and E1I for C20 in the along-track (a), cross-track (b), radial (c) and D (d), Y (e), B (f) direction Fig. 14 SRP acceleration differences between the solution E1 and E2 for C20 in the D (a), Y (b), B (c) direction

Fig. 17
Fig. 17 SLR residuals with respect to ε-angles for the ISL measurements weight ratio of 10 (a-d) or 1000 (e-h)

Table 1
Strategy for L-band and ISL data analysis as well as POD of BDS-3 MEO satellites Solid Earth pole tide: IERS Conventions 2010 (Petit and Luzum, 2010) Third-body gravitation Sun, Moon, Mercury, Venus, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto JPL DE405 is used SRP model See the Table 2 Antenna trust Considered Earth albedo Considered

Table 2
Information on the solutions as well as the corresponding SRP and data used

Table 4
Statistical results of SLR residuals for six solutions (unit: cm)The results from the periods with unstable hardware delays (before DOY 147) are excluded for C30