On the potential contribution of BeiDou-3 to the realization of the terrestrial reference frame scale

Since the release of the phase center calibrations for both the receivers and the satellites, the BeiDou Navigation Satellite System (BDS) became a new potential contributor to the realization of the terrestrial reference frame (TRF) scale of future International Terrestrial Reference Frame releases. This study focuses on the evaluation of the potential usage of the BDS-3 Medium Earth Orbit (MEO) constellation to the definition of the TRF scale. To that aim, we used ground calibrated BDS-3 satellite PCOs provided by the China Satellite Navigation Office and multi-GNSS robot calibrations for the ground antennas conducted by Geo++. Two ionosphere-free linear combinations of signals, namely B1I/B3I and B1C/B2a, have been investigated to find out whether using different frequencies may lead to different TRF scale definitions. Differences between the z components of the satellite phase offsets as given by manufacturer calibrations and those estimated based on IGS14 scale amount to 6.55 ± 12.56 cm and − 0.32 ± 10.99 cm for B1I/B3I and B1C/B2a frequency pairs, respectively. On the one hand, the substantial deviation from the mean reflects the disparities in the quality of calibrations for the individual spacecraft, especially those manufactured by the Shanghai Engineering Center for Microsatellites (SECM). On the other hand, the difference between the two frequency pairs arises to a great extent from the doubtful quality of the SECM PCO calibrations, which certainly do not reflect the frequency dependence of the PCOs. Eventually, the mean scale bias with respect to IGS14 equals + 0.546 ± 0.085 ppb and + 0.026 ± 0.085 ppb for B1I/B3I and B1C/B2a solutions, respectively, when using all 24 BDS-3 MEO satellites.


Introduction
Thanks to the release of the Galileo phase center calibrations for both, the ground (receivers) and space (satellite antennas) segment, Global Navigation Satellite Systems (GNSSs) became a new potential contributor to the realization of the terrestrial reference frame (TRF) scale of future International Terrestrial Reference Frame releases (ITRF, Altamimi et al. 2016;Villiger et al. 2020). The phase center offset (PCO) is a vector between the antenna phase center and a well-defined physical reference point. The PCO is defined by two horizontal components, i.e., PCO-x and PCO-y, and one vertical component, i.e., PCO-z, following the spacecraft axis definitions (Montenbruck et al. 2015b) used within the International GNSS Service (IGS; Johnston et al. 2017). In principle, thze information about the PCO in the z direction (PCO-) is essential, as this component, ideally pointing toward the center of the earth, is in a straight-line relationship with the TRF scale (Zhu et al. 2003;Ge et al. 2005). Villiger et al. (2020) reported that the Galileo-based scale difference with respect to ITRF2014 is 1.4 parts per billion (ppb) at the epoch of January 1, 2018. For comparison, Altamimi et al. (2016) stated that the individual solutions from Satellite Laser Ranging (SLR) and Very Long Baseline Interferometry (VLBI) had a corresponding offset of 0.4 and − 0.4 ppb, respectively. Meanwhile, thanks to the improved procedure of the SLR range bias handling in International Laser Ranging Service (ILRS) solutions, the ITRF2020 results show a substantial scale consistency improvement, as the level of agreement between SLR and VLBI intrinsic scales is now at the level of 0.15 ± 0.05 ppb (Altamimi 2022). Thus, despite the availability of the GNSS-based scale solution, the scale of the ITRF2020 has been established, excluding the GNSS contribution, as the Galileobased scale is now inconsistent with VLBI and SLR to the level of approximately 0.5-0.6 ppb.
The robot calibrations for the ground antennas provided by Geo++ in 2019 covered a wide range of multi-GNSS signal frequencies, including all the GPS, GLONASS, Galileo, BeiDou, and QZSS frequencies (Wübbena et al. 2019). Henceforth, any GNSS can be incorporated into the scale estimation provided that the responsible agency would deliver ground calibrated satellite antenna patterns.
Since GPS week 2056, the BeiDou-3 (BDS-3) frequency and block-specific PCOs informally provided by the Test and Assessment Research Center of the China Satellite Navigation Office (CSNO-TARC; IGS Mail 7782) were available for BDS-3 orbit determination within the IGS. Only half a year later, CSNO officially published metadata information on the BeiDou-2 (BDS-2) and BDS-3, including satelliteand frequency-specific PCOs from ground calibrations. The information about the PCOs for the BDS satellites opened up a space for the second GNSS able to provide an independent realization of the TRF scale. Xia et al. (2020) provided the first validation of the PCO values released by CSNO. For BDS-2, it has been recognized that the CSNO-TARC manufacturer calibrations are in poor agreement with the empirical estimates and the type mean values presently used in the current IGS antenna model (igs14.atx; Rebischung and Schmid 2016) do a much better job. This finding has also been confirmed by Villiger et al. (2021), who found discrepancies of 1-2 m between manufacturer values and their own BDS-2 PCO estimates. Along with the regional nature of the BDS-2 system, the disappointingly low quality of the BDS-2 manufacturer calibrations makes it unreasonable to use them for the TRF scale realization purposes.
For BDS-3 Medium Earth Orbit (MEO) satellites, the results of Xia et al. (2020) and Villiger et al. (2021) confirmed a fairly decent agreement between the PCO estimates from the GNSS network processing and the manufacturer values. Unfortunately, the empirical PCO estimates are so far only available for B1I/B3I but not B1C/B2a and are limited to PRNs of up to C37. The preliminary estimation of the BDS-based TRF scale provided by Qu et al. (2021) demonstrates the inconsistency between BDS and Galileo at the level of about + 1.854 ± 0.191 ppb. Hence, more efforts should be made to bridge the gap and obtain a consistent scale determined by different GNSS constellations.
This study aims at evaluating the potential contribution of the BDS-3 MEO constellation in the definition of the TRF scale based on the ground calibrated BDS-3 satellite PCOs (CSNO 2019a) and the absolute multi-frequency GNSS receiver antenna calibrations delivered by Geo++ (Wübbena et al. 2019). Given the fairly poor orbit quality and regional nature of the BDS-3 Inclined Geosynchronous Orbit (IGSO) and Geostationary (GEO) spacecraft (Zhao et al. 2022), we limited our analysis to BDS-3 MEO satellites to provide the best possible quality of the results.
First, we assessed the quality of the manufacturerreleased PCOs within the BDS-3 MEO constellation, which includes 24 satellites. Second, we evaluated the differences between the TRF scale derived from the BDS PCOs released by CSNO and the scale of IGS14 (Rebischung and Schmid 2016). Two linear combinations of signals, B1I/B3I and B1C/B2a, have been investigated to determine whether using different frequencies may lead to different TRF scale definitions. Since proper modeling of the direct solar radiation pressure (SRP) is a prerequisite for the accurate determination of PCOs, special attention was given to selecting the best processing strategy regarding the latest BDS-3 orbit modeling advances (Duan et al. 2022).

Methodology
We prepared the GPS + BDS-3 global network processing covering the whole year of 2021 based on the observations collected by up to 145 globally distributed ground stations from the IGS multi-GNSS network (Montenbruck et al. 2017). The capabilities of the selected network in tracking the BDS-3 MEO constellation are discussed in the next subsection. The processing was performed using the NAPEOS software (Springer 2009). Table 1 lists the key characteristics of the processing scheme. Our analysis consists of solutions that differ in some individual processing features. The solutions will be introduced in the next subsections.

Estimation of satellite antenna PCO-z and TRF scale
Satellite antenna PCOs can be estimated from GNSS precise orbit determination network solutions. When forming a linear combination of frequencies, the estimated PCO values refer to the combination without a straightforward way to separate the offsets between the individual frequencies. We prepared two kinds of solutions aiming to verify whether the scale realization is consistent for different pairs of frequencies forming a linear combination. The first group of solutions uses B1I/B3I signals (denoted as B1B3), which are heritage signals that are interoperable with BDS-2 satellites. The second group incorporates B1C/B2a signals (denoted as B1B2), which are interoperable with GPS L1/L5 and Galileo E1/E5a. When estimating PCO-z in global network solutions, one should remember about strong correlations between TRF scale, station heights, and the PCO-z; thus, none of these can be estimated simultaneously without additional constraints (Ge et al. 2005). Based on the analysis of GPS network solutions, Zhu et al. (2003) showed that a scale change of + 7.8 ppb, which corresponds to a + 5 cm increase in station heights, can be compensated by a GPS satellite PCO change of − 1 m. Therefore, the PCO estimation requires a priori information about the TRF scale.
The scale of the TRF can be introduced to the GNSS solution twofold. First, the TRF scale can be constrained by a no-net-scale condition with respect to the a priori coordinates of the subset of core stations in the network ). In the second option, the scale of the reference frame is introduced by fixing a given set of satellite PCO values in the solution. For example, introducing the PCOs of one GNSS defines the PCOs of all other GNSSs in the solution via the scale change of the ground network. A totally independent method to estimate GNSS PCOs can be achieved by using GNSS observations on low earth orbiters with well-calibrated phase center locations and the gravitational constraints from the orbital dynamics (Haines et al. 2015;Huang et al. 2020).
In this study, the PCOs of BDS-3 MEO satellites are estimated with the scale introduced from the GPS PCO (igs14_2178.atx), which should strictly realize the IGS14/ ITRF2014 scale. The PCO values are estimated together with satellite orbits, earth rotation parameters, and station coordinates. The station coordinates are estimated with 10 mm constraints imposed on the X, Y, and Z components.

Estimation of horizontal PCO components
The estimation of the horizontal components of the PCO depends on the elevation of the sun above the orbital plane (β angle; Steigenberger et al. 2016). When the sun is perpendicular to the orbital plane (β = ± 90°), the spacecraft's y-axis lines up with the along-track direction, thus causing a full correlation of the PCO-y and the orbit parameters. Considering the same geometry state, the satellite x-axis points toward the sun, which introduces the correlations between the PCO-x and the direct SRP term D 0 when using the Empirical CODE Orbit Model (ECOM, Beutler et al. 1994;Arnold et al. 2015) for SRP modeling. Despite the negligible impact of the horizontal PCO on the scale realization, we discuss the results of the horizontal PCO estimation as an inherent part of the quality assessment of the adopted methodology and the completion of the validation of the BDS-3 PCO model provided by CSNO.

Observations
In general, the tracking capability of the MGEX network for tracking BDS satellites has been extended continuously over See "Solar radiation pressure modeling" section the last decade due to receiver firmware updates and receiver replacements. A rapid change was observed in the middle of 2018 due to the deployment of Trimble Alloy and the update of Trimble NetR9, Septentrio PolaRx5/PolaRx5TR receivers. However, the two main issues remain nowadays in tracking BDS-3 satellites. First, not all deployed BDS-3 satellites can be tracked by an equal number of receivers, mainly due to the limitation of tracking channels or outdated firmware. Second, despite a wide range of frequencies and signals transmitted by BDS satellites, not all of them are supported by the receivers. Figure 1 shows the time series of the number of stations used in daily solutions for processing GPS L1/L2 and BDS-3 B1C/B2a and B1I/B3I observations. We selected a network of 145 stations, all of which track GPS satellites. In 2021, the network included on average 109 and 77 stations tracking BDS-3 B1I/B3I and B1C/B2a signal pairs, respectively. As the availability of the receiver PCO calibration is essential for the estimation of the TRF scale, at least 88% of the selected stations tracking BDS-3 make use of antennas with the multi-GNSS calibration provided by Geo++. The remaining part of the stations employs GPS L1/L2 calibration; however, because of the marginal participation in the network, these should not affect the conclusions of this study. One should note that until February 6, 2021, the number of stations supporting B1C/B2a tracking in the network was too poor to deliver solutions of the comparable quality to the rest of the year; thus, we excluded this set from the next analyses. Figure 2 presents the network of ground sites selected for the analysis with BDS-3 tracking capabilities for the particular frequency pairs. The ground receivers that can track B1I/B1C/ B2a/B3I signals from all the 24 available BDS-3 MEO satellites are Javad TRE_3, Leica GR50, Septentrio AsteRx4 and PolaRx5, and Trimble Alloy. However, some of the PolaRx5 receivers in the network do not provide observations from PRNs greater than C37, most possibly because of outdated firmware or other issues. Trimble NetR9 receivers track only BDS-3 B1I/B3I signals from the satellite channels up to C30. The B1C/B2a signals are not tracked also by individual Trimble Alloys and Septentrio PolaRx5s in the network; thus, one may conclude that even with the technical possibility of tracking specific system frequencies, the receiver has to be properly configured by the responsible receiver operator. Due to the differences in the tracking capabilities of the particular groups of the BDS PRNs, the number of observations from the individual satellites also differs significantly. Figure 3 shows the percentage of observations available for the individual BDS-3 satellites compared to the average number of observations per GPS satellite throughout the analysis period. For B1B3, the mean number of observations with respect to GPS varies from 80% (PRNs up to C30) through 70% (PRNs C32-C37) to 58% (for PRNs larger than C37). In the case of B1B2, the number of observations provided for different groups of satellites is 58% and 52% for PRNs up to C37 and higher, respectively.

Solar radiation pressure modeling
Proper modeling of the SRP is a prerequisite for an accurate determination of PCO values (Steigenberger et al. 2016). In the GNSS literature, three concepts are distinguished in terms of handling direct SRP: (1) analytical (Li et al. 2018;Bhattarai et al. 2022); (2) empirical (Beutler et al. 1994;Bar-Sever and Kuang 2004;Arnold et al. 2015); and (3) hybrid or semi-empirical approach standing somewhere in between to combine a physical understanding of SRP with empirical adjustment of parameters (Rodriguez-Solano et al. 2012a; Montenbruck et al. 2015a;Duan et al. 2019Duan et al. , 2022. Three variants of solutions with different SRP models have been prepared in this study to establish the optimum processing strategy. The first solution is purely empirical and employs a 7-parameter ECOM model (see Table 2). The other two SRP variants utilize the information about the satellites to build an a priori box-wing model and adjust the empirical 5-parameter ECOM model on top of them (see Table 2; Bury et al. 2019Bury et al. , 2020. In this study, the ECOM model is always consistent in definition with the Extended ECOM model as defined by Arnold et al. (2015) with the angular argument Δu, i.e., the argument of latitude of the satellite with respect to the sun's argument of latitude. The geometry and the mass of the satellites are taken from metadata in either case, while the two solutions differ in terms of the adopted optical surface properties. The E1 + BOX META solution relies on the CSNO released metadata. Lack of the reflection and diffusion properties limits the accuracy of SRP modeling; however, the assignment of the remaining coefficient can be concluded based on the research on the properties of individual spacecraft elements Duan et al. 2022;Zhao et al. 2022). Another simplification made in the BDS metadata file suggests that the satellite body of BDS-3 CAST MEO satellites is a standard cuboid, while Zhang et al. (2020) implied that a dedicated T-shaped platform had been developed. In response to the imperfection of the official metadata information, Duan et al. (2022) provided the empirically fitted optical coefficients accounting for self-shadowing caused by a T-shaped bus geometry (denoted as TUM in this study). Although Duan et al. (2022) provided different optical properties for C32/C33 and C36/C37, we decided to use the same a priori model for all the satellites and compensate for the model deficiencies with the ECOM parameters estimated on top. At this point, it is worth noting that along with the box-wing properties, Duan et al. (2022) estimated solar sensor bias, solar panel rotation lag, radiator effect, and a scaling factor of solar panel thermal radiation for each satellite separately. The ideal reuse of the TUM model would require the user to introduce all these values in the processing chain; however, we took the optical properties only to get the impression of the model's usefulness without any major software modification. In general, we should expect better performance obtained for those models with adjusted optical coefficients. However, this is not so straightforward as long as not all the parameters of the model are taken into account. In addition, in the E1 + BOX TUM solution, we considered a constant acceleration of − 1.4 nm/s 2 (600 W) in the − X direction emerging from the glass radiator operation. Table 2 provides a summary of the naming and description of key differences between the solutions. The comparison of the optical coefficients adopted in both box-wing-based models is given in Table 3. One of the main differences between the models concerns the + X surface of the CAST satellites. Based on the discussion given by Zhao et al. (2022), the + X surface is fully covered by multilayer insulation blankets with an absorbed coefficient (α) of 0.35 and a specularly reflected coefficient (ρ) of 0.0. On the other hand, the empirical adjustment of the optical coefficients demonstrated by Duan et al. (2022) points out that most of the photons are specularly reflected with ρ = 0.604. One should also note that the sum of the optical coefficients for the individual surfaces in the TUM model is not equal to 1 as opposed to the META model (Table 3).
Aware of the deficiencies of the ECOM2 model in terms of the PCO determination (Steigenberger et al. 2016), we prepared the ECOM2 solution only in the B1B3 test case. The two remaining E1 + BOX solutions are processed in both B1B2 and B1B3 test cases.

Orbit model validation
The evaluation of the SRP modeling quality was carried out with (1) an independent validation of the orbit using SLR observations and (2) the analysis of the estimated ECOM parameters.

Independent orbit validation using SLR observations
The CSNO BeiDou satellite information file (CSNO 2019b) indicates that all the BDS-3 satellites are equipped with laser retroreflector arrays (LRAs); however, only two BDS SECM-A (C29 and C30) and two BDS CAST (C20 and C22) are actively tracked by the International Laser Ranging Service (ILRS; Pearlman et al. 2019) due to the lack of orbit predictions and no official requirements for providing such observations so far. We generated a dataset of SLR residuals including all available normal points over the analysis period. The results of the analysis are given in Fig. 4. First of all, there are only minor differences between the corresponding B1B3 and B1B2 solutions (in favor of B1B3), which may arise from the differences in the number and geometry of observations used for the orbit determination. Concerning the SRP modeling, using the E1 + BOX TUM solution model results in the smallest offset and standard deviation of SLR residuals for both BDS CAST and SECM-A satellites. Despite the aforementioned simplifications in our use of the Duan et al. (2022) model, the standard deviation of SLR residuals at the level of 24-28 mm is approximately 20-25% better than reported originally by those authors, which may be due to the differences in the outlier screening method. The offset of SLR residuals is positive for the BDS CAST satellites and negative for BDS SECM-A satellites.
However, it is difficult to judge what causes this bias. The most important finding is that the box-wing parameters provided by Duan et al. (2022) can be successfully applied in a straightforward way leading to the best possible orbit quality, surpassing even the usage of the official metadata.

Analysis of the ECOM parameters
The a priori box-wing model should, ideally, account for all non-gravitational perturbing forces acting on a satellite. In both E1 + BOX solutions, the ECOM coefficients are estimated on top of an a priori box-wing model. Therefore, any deviations from zero in estimated ECOM parameters reflect the box-wing model deficiencies. Figure 5 presents the ECOM parameters as a function of the elevation of the sun above the orbital plane (β) for different groups of satellites. The mean value of the estimates has been removed for each satellite separately. The analysis of individual BDS-3 MEO satellites reveals that we may distinguish up to ten different groups of satellites, which are placed on a given orbital plane and are characterized by similar patterns in the estimated ECOM parameters. Official metadata specify only one set of parameters for all the CAST satellites, and apart from that, we know of four CAST satellites equipped with SAR antennas, i.e., C32/C33, C34/C35 (Li et al. 2021). One may expect some distinct SRP signatures because of the SAR antennas, which extend beyond the satellite body. Within the SECM group of satellites, metadata differentiate between SECM-A and SECM-B satellites, which have different dimensions according to the CSNO metadata file. However, the difference is primarily in the Y cross section, while the X and Z areas are almost identical. Moreover, BDS-3 SECM MEO satellites are stretched in a different direction than BDS-3 CAST MEO satellites. Despite using the box-wing model, the Y 0 anomalies for BDS-3 satellites during the eclipse seasons are not well handled as opposed to the findings given by Duan et al. (2022). This indicates that it is not controlled by the box-wing itself but by one of the additional parameters estimated by Duan et al. (2022), e.g., solar panel rotation lag. Using an a priori box-wing model flattens the pattern visible in D 0 , B 0 , and B c ; however, none of the two models diminish the estimated values completely to zero. A key exception is the satellite pair C45 and C46, for which using the box-wing model increases the spread of estimated ECOM parameters, especially for the D 0 term. Besides, the estimated B 0 term for this specific pair shows different variations compared with other CAST satellites, which only confirms its different structure. Figure 6 complements the previous analysis and shows the mean values of the estimated ECOM parameters for each satellite. Approximately 88-97% and 82-93% of a total SRP in the sun direction are covered by the a priori TUM and META box-wing models, respectively, for the individual satellite, while the remaining fraction of 3-18% needs to be compensated by the empirical D 0 ECOM term. Interestingly, there are some substantial differences between the SRP in sun direction estimates for the individual CAST satellites ranging from 126 to 143 nm/s 2 . ECOM parameters as a function of the elevation of the sun above the orbital plane (β) for different groups of satellites based on the satellite manufacturer, orbital plane, and common patterns in the estimated parameters. The mean value has been removed from the data for each satellite separately. All values are given in nm/s 2 Therefore, the estimation of additional empirical parameters on top of the box-wing model is essential for a proper SRP modeling, and the CSNO metadata need more details concerning the disparities within the constellation, aiming for the best results of the orbit modeling and determination of global geodetic parameters (Zajdel et al. 2020b, a).

Results
We estimated the PCO values in different processing scenarios and compared them with the reference values as released in the metadata file. The PCOs are provided by CSNO for each frequency individually; thus, the a priori values for an ionosphere-free linear combination of the frequencies are calculated based on the equation: where f 1 and f 2 are the wavelengths of the frequencies taken for the linear combination, and PCO f 1 , PCO f 2 , and PCO LC are the PCOs for individual frequencies and the linear (1) combination, respectively. The a priori/reference PCOs are given in Table 4. The PCO-x and PCO-y values agree for the two linear combinations at the few-millimeter level. Concerning the PCO-z, there are only a millimeter-level difference between the B1I/B3I and B1C/B2a linear combinations for the BDS-3 SECM satellites, and generally a difference of a few centimeters for the CAST satellites. It is striking to note that for the BDS-3 SECM satellites, the factory-calibrated values for B2 and B3 frequency bands are virtually the same as those of B1. As a consequence, the ionospheric-free combination of the nominal PCO values is almost identical for B1I/B3I and B1C/B2a. Table 5 provides the formal estimation error of PCOs in different solutions. Table 5 is supported by Fig. 7, which presents the time series of PCO estimation errors in E1 + BOX TUM and both B1B2 and B1B3 test cases. The formal errors mainly reflect the determinability of the PCOs due to the sun-earth-satellite geometry. First, we see the difference in the estimation errors between the corresponding B1B2 and B1B3 test cases, which may be attributed to the change in Fig. 6 Mean values of the estimated ECOM for different groups of satellites. The grouping is based on satellite manufacturer, orbital plane, and common patterns in the estimated ECOM parameters according to Fig. 5. All values in nm/s 2 the number of observations and the geometry of the network; however, it affects the estimation error only at the level of a few millimeters. In the case of high PRN numbers such as C45/C46, which are tracked by almost the same number of stations in both B1B2 and B1B3 solutions, the formal errors are almost the same (Fig. 7). There is also no difference in the formal error between the CAST and SECM satellites. The estimation error for the horizontal PCO components is growing in parallel with the elevation of the sun above the orbital plane, especially for the PCO-y. The BDS-3 MEO plane A is the most privileged in terms of the estimation of horizontal PCO, with the β angles always below 30° (Fig. 7). Using ECOM2 causes an increase in the errors for the PCO-x approximately by a factor of 4 compared to the E1 + BOX solutions. The formal errors of the PCO-y and PCO-z do not depend on the SRP modeling approach. The secular decrease in the PCO-z formal error reflects the gradual increase in the number of tracking stations in the network (Fig. 2). It amounts to roughly 5 mm or 10% of the error near the beginning of the analysis period. Finally, by the representative example of the B1B3 E1 + BOX TUM solution, we may state that the formal estimation error is at the level of 0.6-1.3, 0.9-3.3, and 3.2-3.6 cm for PCO-x, PCO-y, and PCO-z, respectively. Figure 8 illustrates PCOs estimated in the E1 + BOX TUM solution with respect to reference values from CSNO metadata. Only periods corresponding to the |β| ≤ 45° provide stable horizontal PCO estimates, whereas variations of up to ± 20 cm appear for periods outside this range. The only exception is the group of the satellites C45/C46, for which all the estimated PCO components vary in time depending on the β angle. That confirms the different structures of these satellites and the increased impact of the estimated ECOM parameters, which have to compensate for the box-wing model deficiencies. Figure 9 presents the PCO estimation results in box-whisker plots for all the satellites and all the considered solutions. We subdivided the satellites into groups consistent with the previous analysis of the patterns in the ECOM estimates. Additionally, for comparison, we added next to our results the PCO values obtained by Qu et al. (2021) for the C19-C37 satellites based on the global network processing using the linear combination of B1I/B3I signals. For PCO-x, the estimated values agree with the ground calibrations within 1 and 2 cm, for the CAST and SECM satellites, respectively. This excludes the C45/C46 pair, for which an 8 cm offset is visible. The mean offset in the PCO-x estimates for this specific pair of satellites disappears when using ECOM2; however, this SRP approach is not suited for the determination of PCO-x, as visible in the spread of the estimated values for all the remaining satellites. For PCO-y, the estimated values are on average consistent to the level of 1 cm with the ground calibrations, but with the standard deviation of estimates reaching 8 cm for the satellites on the orbital planes B and C, with wide β angle ranges.

Analysis of the estimated PCOs
More attention should be devoted to the PCO-z estimates, which directly translate into the potential contribution of the BDS-3 MEO constellation to the TRF scale definition. In general, the PCO-z estimates have a larger scatter than the horizontal PCOs, but they do not exhibit pronounced systematic effects. The standard deviations of the time series are on the 4-6 and 5-7 cm levels for B1B2 and B1B3 solutions, respectively. A comparison of these with formal errors reveals that those were too optimistic by a factor of two. Even though the estimation error is lower for B1B3 linear combination, the less scattered estimates for the B1B2 solution may arise from the fact that using B1C and B2a signals in the processing causes better signal strength and higher accuracy for both code and carrier phase observations (Zhang et al. 2021).
Besides the previous division of the BDS satellites based on the ECOM patterns, we may distinguish three main subgroups of CAST satellites based on notably different values of PCO-z from the ground calibration. While the first eight CAST satellites (C19-C24, C32-C33) have a B1 factory PCO-z value close to 1.50 m, it is around 1.35 m for the next four satellites and 1.05 m for the last two (group C41/C42). One should also note that the difference between B1B2 and B1B3 is very similar for all CAST satellites; thus, it appears reasonable to assume that all CAST satellites share the same type of antenna array.
Concerning CAST satellites, the major clash between official calibrations and estimated PCOs refers to the pair of C41/C42 satellites reaching about 30-35 cm. Interestingly, the LRA offsets in the Z direction amount to 1.25 m, 1.10-1.17 cm, and 1.10 cm for the three aforementioned groups. This corresponds to the PCO difference between the first two subgroups, but the third group is more close to the second one rather than forming a separate group. Therefore, we may assume that it is far more likely that the C41/C42 satellites suffer from a 30 cm change in the antenna reference point used for antenna calibration rather than a true change in the PCO. Clarification of this apparent inconsistency has been requested from CSNO-TARC, but was not received at the time of writing.
Another interesting observation is that we see an orbital plane dependency in the mean offset observed in the PCO-z estimates for the CAST satellites. The CAST satellites orbiting plane B have a bias in the estimated values at the level of − 10 to − 8 cm, while the bias for the CAST satellites on plane C is close to zero. The differences between the B1I/ B3I and B1C/B2a corrections with respect to the nominal values are within 0-4 cm for all the CAST satellites. At present, it is unclear whether the orbital plane dependence of the estimated PCOs is merely coincidental or hints at a hidden deficiency in the SRP modeling, such as the neglect of shading or unknown surface elements. Within the SECM satellites, we already differentiate between SECM-A and SECM-B satellites. Nevertheless, the estimated PCO-z of the SECM satellites show a substantial scatter for individual satellites with no orbital plane or satellite subtype dependence. Furthermore, the estimated PCO-z values show a 10-20 cm difference between B1I/B3I and B1C/B2a. For SECM-B satellites, the difference is slightly larger than for the SECM-A satellites. Ultimately, all that we can conclude is that the observed PCO-z for the BDS SECM satellites show notable scatter and systematic deviations from the nominal values. In contrast to the ground calibrated nominal values, the observed SECM PCO-z exposes an obvious frequency dependence. We might speculate that the SECM satellites suffer from large near-field effects causing systematic PCO differences between factory calibrations of the standalone antenna array and the antenna array as integrated with the satellite.
Comparing the estimated values with the results presented by Qu et al. (2021), we may notice that our PCO estimates for the CAST satellites are less scattered, especially for the group of C32-C37, for which Qu et al. (2021) found the discrepancies with respect to the nominal values reaching 20 cm. For the SECM satellites, on the other hand, there is an almost perfect match between our estimates and those provided by Qu et al. (2021), which only confirms the questionable quality of the ground calibrations provided by SECM.
Taking the constellation as the whole, the mean PCO-z offset with respect to nominal values is 6.55 ± 12.56 and − 0.32 ± 10.99 cm for B1I/B3I and B1C/B2a frequency pairs, respectively. The "misbehavior" of the SECM calibrations should then be indicated as the leading source of the inconsistency in mean PCO-z shift and the scale definition depending on the usage of either B1I/B3I or B1C/B2a frequency pairs. The constellation mean PCO-z offset considering only CAST satellites except for C41/C42 changes to Finally, we considered two methods to derive a single set of PCO values per satellite, i.e., (1) weighted average based on formal estimation errors, (2) unweighted average with the exclusion of periods with systematic effects (|β| > 45°). Both methods agree within a few millimeters; thus, we chose the first approach to deliver the final values, as shown in Table 6. The formal errors of the PCO estimates are on the 1 cm level for the horizontal PCOs and the 3-4 cm level for the PCO-z. Therefore, the PCOs presented in Table 6 were rounded to 1 cm. Figure 10 shows the time series of the TRF scale change derived from the Helmert transformation between the two global network solutions with the scales defined by the GPS PCOs (igs14_2178.atx, IGS14/ITRF2014 scale) and the ground calibrated BDS-3 PCOs released by CSNO. In both cases, the station coordinates are estimated with 10 mm constraints imposed on the X, Y, and Z components. The mean scale bias equals 0.546 ± 0.085 ppb and 0.026 ± 0.085 ppb, for B1B3 and B1B2 solutions, respectively. We may assume that the scale discrepancy between the B1B3 and B1B2 solutions arises to a great extent from the uncertain quality of the SECM PCO calibrations, which certainly do not reflect the frequency dependence of the PCOs. Moreover, the bias of 0.546 ppb delivers a decent agreement with Qu et al. (2021), who reported the 0.446 ± 0.153 ppb offset with respect to the ITRF2014 based on 30 days of processing B1I/B3I observations from both BDS-2 and BDS-3 satellites.

Realization of the terrestrial reference frame scale based on the CSNO BDS-3 PCO model
Furthermore, we analyzed how the scale change affects the height component of the station coordinates. Figure 11 presents the time series of differences (2)-(1) in the height component of station coordinates between the two global network solutions: (1) with GPS PCOs fixed (igs14_2178. atx, IGS14/ITRF2014 scale) and (2) the ground calibrated BDS-3 PCOs released by CSNO fixed (consistent with Fig. 10). Different selections of stations have been tested to find out, to which extent GPS + BDS-3 and GPS-only stations contribute to the scale transfer within the network. The mean difference observed in the height component equals approximately 3.4 ± 0.6 and 0.2 ± 0.4 mm for B1B3 and B1B2 solutions, respectively. One should note here that the mean offset is almost the same whether all the stations are selected, only those which track GPS + BDS signals, or even those which do not support BDS tracking at all. It confirms that the scale defined by BDS-3 PCO-z is transferred to the whole network, including also GPS-only stations. When BDS-3 PCOs are tightly constrained in the solution, the estimated GPS PCO-z are adjusted to the new TRF scale standard based on the common network of ground stations. Then, the scale transfer can also proceed to the stations that track only GPS signals.
The variations of the station heights and the scale of the TRF are approximately proportional to the mean offset in PCO-z. Therefore, based on the mean PCO-z offset, the corresponding scale offsets, and the mean offset in the height component of the station coordinates, we may formulate the relationship between them similarly to Zhu et al. (2003). For the B1B3 solution, the 0.546 ppb scale change corresponds to the scale factor of 8.3 ppb/m concerning the mean PCO-z shift of 0.0655 m. The factor of 8.3 is slightly higher than the 7.8 ppb/m reported by Zhu et al. (2003). Going further, the ratio of station height change and BDS PCO-z offset is − 0.052, i.e., − 5.2%. Contrary, Zhu et al. (2003) stated a 5% shift in station height relative to the GPS PCO-z shift, which is also slightly lower than reported here. If all satellites have identical PCO-z shifts, we would expect mainly the scalelike systematics (Zhu et al. 2003). Unfortunately, because of the substantial differences in mean PCO-z shift values of various satellites within the BDS-3 constellation (especially SECM satellites), the formulation of conclusive and precise numerical dependencies may not be so straightforward for BDS and would require further investigations.
To sum up, Fig. 12 presents the calculated PCO values versus the initial values given by CSNO. Each point represents one of the 24 MEO satellites with blue indicating CAST satellites and red indicating SECM satellites. Points above the diagonal mean that the factory calibration gives a larger PCO than observed with the scale defined by the Fig. 9 Differences between the PCO values estimated in different test cases with respect to the initial values provided by CSNO (2019a). The individual satellites are described by the manufacturer, orbital plane, and PRN number IGS14-based GPS PCO. Thus, switching to the BDS-3 factory values would imply a reduction in the station height and cause a scale reduction. Vice versa, a positive scale change would be obtained when using factory calibrations for satellites below the diagonal. The gray dashed lines refer to the ± 1 ppb scale change. For the ratio between PCO shift and TRF scale/station height change, we adopted a value of − 5.2%, i.e., a 1 m PCO-z increase implies a 5.2 cm station height reduction or −8.3 ppb scale change, as derived empirically from the analysis of the B1B3 solution.
Comparing the estimated PCOs, we obtain the consistency within ± 1 ppb of the CSNO values with the estimates for all the CAST satellites except for the pair of C41 and C42 satellites. Despite mostly identical spacecraft, the offsets of observed and ground calibrated PCOs for SECM exhibit a substantial scatter and are notably different for the two frequency pairs.

Conclusions and discussion
This study focuses on the evaluation of the potential contribution of the BDS-3 MEO constellation in the definition of the TRF scale. To that aim, we used ground calibrated BDS-3 satellite PCOs delivered in 2019 by CSNO and multi-GNSS robot calibrations for the ground antennas provided First, we assessed the quality of the manufacturerreleased PCO values considering the full constellation of 24 BDS-3 MEO satellites. The PCOs of BDS-3 MEO satellites were estimated with the scale introduced from the GPS PCOs (igs14_2178.atx), which strictly realizes the IGS14/ITRF2014 scale. The mean of the estimated PCO-z offsets is 6.55 ± 12.56 cm and − 0.32 ± 10.99 cm for B1I/ B3I and B1C/B2a frequency pairs, respectively. The substantial dispersion from the mean reflects the disparities in the quality of calibrations for the individual spacecraft. Besides the nominal division of the BDS MEO satellites by manufacturers (CAST and SECM), we may distinguish three main subgroups of CAST satellites based on notably different values of PCO-z from the ground calibration. The most notable discrepancy between official calibrations and estimated PCOs occurs for the pair of C41/C42 satellites, reaching about 30-35 cm. Based on the comparison with LRA coordinates, we suspect a 30 cm change in the antenna reference point used for antenna calibration rather than a true change of the PCO. The CAST satellites orbiting plane B have a consistent bias in the estimated PCO-z at the level of − 10 to − 8 cm, while the bias for the CAST satellites on plane C is close to zero. The differences between the B1I/ B3I and B1C/B2a corrections with respect to the nominal values are within 0-4 cm for all the CAST satellites. On the other hand, the estimated PCO-z estimates for the SECM satellites show a large scatter for individual satellites with neither orbital plane nor satellite subtype dependence. Furthermore, the estimated PCO-z values reveal a 10-20 cm difference between B1I/B3I and B1C/B2a.
Finally, we assessed the scale difference as realized by BDS-3 MEOs when using the PCO model provided by CSNO. The mean scale bias with respect to ITRF2014 equals + 0.546 ± 0.085 ppb and + 0.026 ± 0.085 ppb, for B1I/B3I and B1C/B2a solutions, respectively. The difference between the two frequency pairs results from the doubtful quality of the SECM PCO calibrations, which certainly do not reflect the actual frequency dependence of the PCOs. On the one hand, the presented results are in poor agreement with Galileo, which realizes the TRF scale with a bias of approximately 1.4 ppb with respect to ITRF2014 when making use of factory-calibrated transmit antenna PCOs. On the other hand, the mean scale bias of 0.546 ppb is seemingly close to the scale difference between ITRF2020 vs. ITRF2014, which accounts for 0.42 ppb (https:// itrf. ign. fr/ en/ solut ions/ trans forma tions). However, it should be noted that using, for example, only the most consistent group of BDS-3 CAST satellites, for which the mean PCO-z offset is at the level of − 3 cm, would lead to another scale definition with a bias of roughly − 0.3 ppb. Moreover, the scale discrepancy between the B1I/B3I and B1C/B2a frequency pairs puts into question the quality of satellite PCO calibrations published by CSNO, especially concerning the SECM