BeiDou satellite radiation force models for precise orbit determination and geodetic applications

—China’s BeiDou satellite navigation system (BDS) has completed its full constellation in orbit since June 2020. Services have been evolved from regional (BDS-2) to global (BDS- 3). This contribution evaluates the impact of solar radiation pressure (SRP) modeling on satellite orbits and geodetic parameters. To that end, we process 2 years of BDS observations (2019-2021), collected by a network of 100 ground stations. A physical a priori box-wing (bw) model based on the estimated optical properties is introduced. Various physical effects, such as yaw bias, self-shadowing, radiator emission and thermal radiation of solar panels are considered. The ECOM (Empirical CODE orbit Model, 5 parameters), ECOM+along-track and ECOM2 (both 7 and 9 parameters) models are employed on top of the a priori box-wing model in the experiment. We show that without the use of the a priori box-wing model, the ECOM+along-track model shows clear better orbit solutions during eclipse seasons for BDS-3 satellites. This is proven to be mainly due to the thermal radiation of solar panels. However, the along-track acceleration is highly correlated with LOD (length of day) and ECOM parameters. LOD estimates in this case are contaminated. The STD (standard deviation) of daily LOD estimates with respect to IERS-C04-14 series increases from 40 µ s (ECOM) to 85 µ s (ECOM+along-track). After the consideration of the a priori box- wing model, satellite orbital errors are greatly reduced for all the ECOM models. For instance, orbit misclosures of BDS-3 CAST (China Academy of Space Technology) satellites improve by a factor of two for the ECOM model during eclipse seasons; dependencies of SLR (satellite laser ranging) residuals on the sun elongation angle almost vanish for BDS-3 satellites. Furthermore, the use of the a priori box-wing model mitigates a great majority of the spurious signals in the geodetic parameters. In particular, the total amplitude of the 1, 3, 5, 7 cpy signals for the geocenter Z component has been reduced by a factor of 4.5 for the ECOM model. In general, the combination of the introduced physical a priori box-wing model and the ECOM model is preferred for BDS satellites.


I. INTRODUCTION
BeiDou satellite navigation system (BDS, originally called COMPASS) is developed in three generations (BDS-1, BDS-2, BDS-3). BDS-1 consists of two operational GEO (Geostationary orbit) satellites and an additional one for backup. BDS-1 started providing positioning, timing and short message services in 2003 [1]. With range measurements from two satellites, users cannot compute their 3D positions directly. Bingbing  The measured ranges are first transmitted to the ground control center, then are mapped with the topographic map, and finally the computed 3D positions are transmitted back to the users [2]. The positioning accuracy from BDS-1 is about 20 m [3]. BeiDou-1 was decommissioned at the end of 2012. BDS-2 aggregates 5 GEO, 5 IGSO (Inclined Geosynchronous Orbits) and 4 MEO (Medium Earth Orbits) satellites. BDS-2 became operational in 2012, serving mainly the Asia-Pacific area. As presented by [4]- [8], the radial accuracy of BDS-2 satellites is better than 10 cm for the IGSO and MEO satellites while is about 50 cm for the GEO satellites. However, during the eclipse seasons, especially during the attitude-turn maneuver periods, orbit radial accuracy can be two times worse. The global BDS-3 consists of 3 GEO, 3 IGSO and 24 MEO satellites. Before the launch of the formal BDS-3 system, a demonstration (or experimental) BDS-3S system with two IGSO and three MEO satellites was established between 2015 and 2016 [3]. The purpose is to test new payloads (hydrogen maser clocks), new signals (for instance the S-band signal) and new techniques (inter-satellite link) [9]- [14]. With this demonstration system, Zhang et al. [15] confirm that the satellite-induced code biases present in BDS-2 satellites [16] almost vanish for BDS-3S satellites. Also, triple-frequency signals perform better than dual-frequency signals for ambiguity resolution [17]- [20].
With the launch of the last GEO satellite in June 2020, BDS-3 has completed its full constellation in orbit (http://www.beidou.gov.cn). Since 2019, BDS-3 ground tracking capabilities of the IGS (International GNSS Service) MGEX (Multi-GNSS Experiment) network have been significantly enhanced. As reported by Steigenberger and Montenbruck [21], four types of receivers (Javad TRE 3 3.7.6, Septentrio PolaRx 5 5.3.0, Trimble NetR9 5.42 and Trimble Alloy 5.42) track BDS-3 B1I and B3I signals for the PRNs up to C37. Only Trimble Alloy and the latest released Javad 3.7.8 (or newer version) receivers are capable of tracking BDS-3 PRNs up to C63. To know the current tracking capabilities of the IGS network for BDS-3 satellites, we download 1-day RINEX files of all the IGS multi-GNSS tracking stations in June 2020. More than 100 IGS stations track BDS-3 signals on B1I (C2I/L2I) and B3I (C6I/L6I) frequencies simultaneously for PRNs up to C37. However, tracking stations for PRNs up to C67 are not yet sufficient, with less than 10 stations observed simultaneously for one satellite (for instance PRN C41 in June 2020). Consequently, studies in this contribution focus on BDS satellites with PRNs up to C37.
In addition to the capability of ground tracking network, solar radiation pressure (SRP) is another significant factor affecting satellite orbits and geodetic parameters. This perturbation can be modeled both analytically and empirically. Milani et al. [22] formulate the physical interaction between SRP and satellite attitude, dimensions, total mass and optical properties. In December 2019, CSNO (China Satellite Navigation Office) published this metadata information for BDS satellites [23], [24]. However, the diffusion and reflection properties of satellite surfaces and solar panels are not given. In the absence of a precise analytical model, the empirical ECOM model is widely used [25]. The ECOM parameters are defined in a Sun-oriented DYB frame, with D pointing towards the Sun, Y along the solar panel axis and B completing the right-handed system. The classical 5-parameter ECOM model consists of three constants (D0, Y0, B0) and two firstorder Fourier coefficients (B1C and B1S) in the B direction [26]. This model is later extended as a new ECOM2 model considering higher order Fourier terms in the D direction [27]. It is confirmed by Prange et al. [28] that the ECOM2 model shows much better performance than the ECOM model for Galileo and QZSS (Quasi Zenith Satellite System) satellite orbits, which both have clear elongated shape for the satellite body. For GLONASS satellites, CODE suggests excluding the 4-th order terms (D4C, D4S) and using a 7-parameter ECOM2 model since solutions could be contaminated by unnecessary parameters due to correlations [29].
As proven by [30]- [37], an a priori box-wing model not only reduces orbital errors but also mitigates spurious signals in the spectra of the geodetic parameters. Therefore, the main purpose of this research is to set up a physical box-wing model for BDS satellites. We first estimate satellite surface optical properties based on BDS measurements over two years. Various physical effects, such as yaw-bias, self-shadowing, radiator emission and thermal radiation of solar panels are considered in the adjustment. All the adjusted parameters are then employed in a physical box-wing model, which is jointly used with the ECOM, ECOM+along-track and ECOM2 model, respectively. Finally, performances of all the tested SRP models are assessed by looking into the ECOM estimates, orbit misclosures, orbit predictions, SLR residuals and geodetic parameters.

II. BEIDOU SATELLITE RADIATION FORCE MODELS
The box-wing model is widely used to describe SRP and earth radiation (albedo and infrared) of GNSS satellites, (1) where A denotes the surface area, M the total mass of the satellite, S 0 the solar flux, c the vacuum velocity of light, κ the thermal reradiation factor (0 for solar panel and 1 for satellite body surfaces in this contribution), α, δ, ρ the fractions of absorbed, diffusely scattered, and specularly reflected photons. Furthermore, e D denotes the incident vector direction, e N the surface normal vector (depending on satellite attitudes), and θ the angle between both vectors. The attitude mode of BDS-2 IGSO and MEO satellites switches from yaw-steering to orbit normal when the absolute value of the sun β angle is smaller than 4 deg [38]. However, for the new BDS-2 IGSO-6 satellite, Dilssner [39] shows that it does not change to the orbit normal mode but switches to a smoothed yaw-steering mode when the absolute β angle is smaller than 2.8 deg. BDS-3 satellites are built by two manufactures, CAST (China Academy of Space Technology) and SECM (Shanghai Engineering Center for Microsatellites). BDS-3 CAST MEO satellites take the same attitude law as the new BDS-2 IGSO-6 satellite [40], whereas, for BDS-3 SECM MEO satellites, the yaw angle computation is based on a β angle fixed to 3 deg (with the same sign as the true β angle) when the absolute β angle is smaller than 3 deg [41].
Solar panels and satellite body surfaces are treated differently in equation (1). The satellite body is usually externally covered with multilayer insulation (MLI) blankets in order to isolate the satellite from the external environment. External heat on satellite surfaces in this case is assumed to be immediately radiated back into space. However, electronic devices (such as clocks) generate heat inside the satellite. This heat raises the internal temperature and could exceed the acceptable limits. Therefore, it is necessary to release this heat through openings in the blanket, for instance through optical solar reflectors (radiators). As GNSS satellites keep transmitting signals continuously, we can assume that heat generated by the equipment is constant over time. In this case, the radiator effect can be easily considered as an additional constant acceleration, where R is the radiator acceleration. Thermal radiation of solar panels is more complicated. If the satellite is continuously illuminated thermal radiation of solar panels should be constant over time. The ECOM D0 parameter can fully absorb this effect. However, if the satellite is inside eclipse seasons thermal radiation of solar panels must have to be carefully modeled, because thermal radiation performs differently when solar panels start cooling down (enter into shadows) and heating up (exit shadows). To model thermal radiation of solar panels, we first need to compute the temperature difference of the two sides of the solar panels. Then, according to the emissivity on both sides of the solar panels we compute the thermal radiation acceleration.
The temperature variation of the solar panel may be modeled by balancing the input energy flux from the solar radiation and the thermal radiation on both sides of the array where C A is the heat capacity per unit area, i.e., the heat capacity per mass multiplied by the density and thickness of the layer, accumulated for the different material layers of the panel. The first term on the right-hand side of equation (3) gives the heating input energy flux containing the sun illumination factor γ, the electric efficiency of the panel η and the absorption coefficient α. The second term describes the radiative cooling on both sides of the panel with the Stefan-Boltzmann constant σ and the emissivities ε on the front and back side of the panel. The parameter ξ models the temperature difference of the front and back side of the solar panel, which is related to the thermal conductivity through the panel.
For simplicity, we assume that the temperature difference is proportional to the temperature while in principle the heat diffusion equation should be solved, a partial differential equation of second order. The simplifications are compensated by estimating a scaling parameter s. In order to get a stationary solution, the differential equation (3) has to be solved as a boundary value problem, requesting the same temperature after one satellite revolution. With the given temperature the thermal acceleration can be computed as As described by ShangHai Institute of Space Power-Sources, BDS-3 satellites use the triple junction GaAs solar cells (http://en.811sisp.com). This cell type is a GaInP 2 /GaAs/Ge assembly on Ge substrate. The total thickness of the solar panel is about 2 mm, the absorption property is 0.92 and the efficiency is more than 30%. Thermal properties of the same materials are taken from [42], [43], as shown in Table I. We assume that the solar cells are mounted on a carbon fiber honeycomb structure to provide the mechanical stability but having a negligible contribution to the heat capacity of the solar panel. Emissivity information of BDS solar panels is unknown. We assume an emissivity of 0.75 on the front side and 0.90 on the back side. In addition we use a difference of front-and backside temperatures of 2% and a panel thickness of 2 mm. In order to compensate for errors in these assumptions and simplifications in the model we estimate an acceleration scaling s as part of the orbit determination procedure. Computing orbit solutions for different panel thicknesses from 2 mm to 1 cm we verified that the scaling parameter in fact absorbs differences in the model assumptions and that the resulting orbits show a low sensitivity for varying model input parameters. Fig. 1 shows the temperature variation of solar panels for one BDS-3 CAST satellite when entering into and exiting the shadow based on the above assumptions. The obtained constant value in the sunlight is 316 K (43 • C) for the front side and is 309 K (36 • C) for the back side. The temperature behavior is very similar to that reported for GPS Block IIA and IIR solar panels [44], [45]. Accordingly, the thermal emission difference is 26 W/m 2 , causing a thermal thrust acceleration of 1.7 nm/s 2 away from the sun, as shown in Fig. 2. The cooling down part in Fig. 2 can be obtained as an exact analytical solution of the differential equation (3) while the heating up part has no analytical solution but can be approximated by an exponential function. For the sake of a simple usage, we provide directly the function coefficients of these two parts. This allows an easy adoption of these functions in other software packages without the need for a deep understanding of the physical background.
The analytical solution of the cooling down part is  where acc shd denotes the thermal acceleration inside shadows, p and τ c the function coefficients, t the current time epoch in seconds, t 0 the shadow entering epoch and s a scaling factor. The empirical exponential function for the heating up part is with a time shift parameter assuring optimal fit to the solution of the differential equation after shadow exit where acc sun denotes the thermal acceleration outside shadows, t h and τ h the function coefficients and t sha the shadow duration in seconds. The cooling down and heating up function coefficients for different BDS satellites are given in Table II.  The table also gives the scaling parameter mentioned above that compensates for uncertain thermal properties in Table I and inaccurate assumed emissivity of the solar panels. The parameter is estimated for each type of BDS satellites based on BDS measurements as part of the orbit determination procedure. The impact of thermal radiation on satellite orbits is shown in section IV. We compute precise satellite orbits using the ECOM model as a first step. The geodetic datum is defined by constraining station coordinates tightly to values in the IGS14 frame of the same station and epoch. The satellite body-fixed frame follows the IGS convention [46]. Earth rotation parameters are fixed to the values of IERS Bulletin A in this step. Earth radiation [47] and satellite antenna thrust [48] are considered. BDS satellite L-band transmit power is taken from the IGS metadata SINEX file (IGSMAIL-8015, from Peter Steigenberger). Satellite antenna phase center offsets and variations (PCOs, PCVs) are corrected using the igs14.atx file (IGSMAIL-7782, from Arturo Villiger). Station PCO and PCV corrections for BDS B1I and B3I frequencies are assumed to be the same as those for GPS L1 and L2 frequencies. BDS-2 and BDS-3 satellite undifferenced phase ambiguities are resolved individually at single difference level [49] because there are clear receiver type specific system-like biases between BDS-2 and BDS-3 pseudorange signals [50], [51]. Detailed settings are shown in Table III. In the second step, we estimate BDS satellite box-wing parameters while introducing integer ambiguities and station related parameters determined in the first step as known. All the computed BDS satellites have been in orbit for more than half a year at the computation starting epoch. Outgassing effects are greatly reduced in this case. The total mass of each BDS satellite is fixed to the officially published value except  Fig. 4. In this case, the red color part in the +X surface could be shaded by the blue color step part if the sun elongation angle is between 0 and 90 degrees. Part of the -Z surface is shaded as well expect for the case when the elongation angle is close to 180 degrees. Moreover, Chen et al. [54] show that there could be an additional panel on the -X surface (green color in Fig. 4) flipped up when the satellite is in orbit. The additional surface will not cause shadows in the +X surface but will greatly change the cross section area of the +Z and -Z surfaces. Furthermore, self-shadow effects in the -Z surface increase as well. To obtain correct shape and dimensions, we test three cases in the optical parameter adjustment. (1) The CSNO published dimensions without self-shadow effects. (2) Dimensions in Fig. 4 considering self-shadow effects. (3) Case 2 considering additionally the flipped surface. Same as we did for other GNSS satellites [55], [56], we estimate satellite optical properties (absorption+diffusion, α + δ and reflection, ρ) in +X and ±Z surfaces of the satellite. The satellite solar panel ρ is estimated to model the constant thermal radiation in  sunlight. Solar sensor bias, solar panel rotation lag, radiator effect and scaling factor of solar panel thermal radiation are estimated for each of the satellite [57], [58]. Fig. 5 shows the adjusted optical properties of BDS-3 CAST MEO satellites for three cases. Satellite specific solutions are averaged into three groups according to the optical property estimates. For satellite group C19 to C24 and C36 to C37, α + δ + ρ in case 2 is more close to 1, the physically correct value. For satellite group C32 to C33, case 1 and case 2 show clear unphysical negative estimates of α+δ in the +X surface while case 3 shows realistic estimates. More information from CSNO demonstrates that satellites C32 and C33 for the first time carry SAR (Search and Rescue) antennas, which might be mounted on the additional surface. Thus, in our following computation, we use case 3 for satellites C32 and C33 and use case 2 for the other BDS-3 CAST MEO satellites. All the dimensions, total mass and the estimated box-wing parameters are given in Table IV, V, VI and VII. In addition to the optical properties, we estimate additional box-wing parameters such as radiator accelerations and attitude biases. Different than optical properties, these parameters are similar for the satellites of the same type. BDS-3 CAST satellites show a clear radiator acceleration of -1.4 nm/s 2 along the body-fixed -X direction, corresponding to an emission power of around 600 W. By checking the CSNO documents, we find that the -X surface of BDS-3 CAST satellites in fact consists of two parts. It seems that one part with an area of 1.11 m 2 is a glass radiator. Table VIII shows solar sensor bias (that causes a yaw bias) and solar panel rotation lag of BDS satellites. The scaling factors of solar panel thermal radiation for each type of BDS satellites are given in Table II. The precision of all the adjusted box-wing parameters is evaluated by assessing the repeatability of three solutions. The three solutions use data every 3 days over 2 years, each shifted by one day with respect to the other. All the three solutions agree with each other fairly well, for instance the mean STD of all the optical property parameters is about 0.012.

MODELS
We use all the estimated box-wing parameters in a physical a priori box-wing model, and jointly use it with ECOM, If the a priori box-wing model is perfect the estimated  Orbit misclosures at the day boundaries between consecutive arcs are traditionally used to assess the orbit internal consistency. In order to see the impact of thermal radiation of solar panels, We first test ECOM (no box-wing model), ECOM considering thermal radiation of solar panels and ECOM+along-track model (ECOM model with an additional constant acceleration parameter in the along-track direction), respectively. We compute orbit differences at the day boundaries between consecutive arcs and define the orbit misclosure of two sets of orbits as where i and i + 1 refer to days, n the number of satellites  Table X. The bw+ECOMA model shows the best orbit solutions, followed by the bw+ECOM model and then the bw+9ECOM2 model. The reason could be that thermal radiation of MLI has a small effect as well, and an along-track acceleration can absorb it fairly well.
To further assess each of the SRP models, we predict BDS The SLR technique is usually adopted as an external reference to evaluate satellite orbit products [60], [61]. Four BDS-3 satellites (two CAST C20, C21 and two SECM C29, C30) have SLR measurements during the experiment time period. Antenna information of the LRA (laser retroreflector array) follows values published by CSNO [23], [24]. Fig. 9 and Fig.  10 show the dependencies of SLR residuals on the elongation angle for the ECOM model with and without the box-wing model. It is evident that such dependencies almost vanish for all the BDS-3 satellites when using the box-wing model. Mean offsets and STD values regarding different SRP models and satellite groups are given in Table XII (SLR residuals outside ±50 cm are not used for the statistic). The box-wing model reduces mean offsets and STD values of SLR residuals for all the BDS satellites and ECOM models. The bw+ECOM model results in the best STD value of SLR residuals.

V. IMPACT OF SRP MODELS ON GEODETIC PARAMETERS
GNSS is one of the four space geodetic techniques contributing to the construction of the ITRF (International Terrestrial Reference Frame). The advantage of the GNSS technique is the continuously available measurements and the globally distributed tracking stations. Altamimi et al. [62] show that the GNSS-based (GPS and GLONASS) X and Y pole coordinates dominate the combined solution. However, the GNSS technique is not able to estimate the UT1 parameter due to the high correlations with the rotation of the orbital planes. Nonetheless, with the GNSS technique, the LOD (length of day), i.e., the negative time derivative of UT1-UTC can be estimated. Zajdel et al. [36] present the performance of Galileo satellites for the earth rotation parameter (ERP) estimation. In this section, we show the BDS satellite-type-specific geodetic parameter solutions with different SRP models. All the ERP estimates are compared to the values from IERS-C04-14 [63]. The reason is that BDS-3 consists of 18 satellites and all are MEO satellites while BDS-2 consist of 10 satellites and 7 are IGSO satellites. The BDSALL based solution performs slightly better than the BDS-3 based solution. The use of the a priori box-wing model makes all the differences more smooth over time. Table XIII shows the mean offset and standard deviation of the ERP residuals with respect to the IERS-C04-14. The BDS-3 based X and Y pole coordinates are about two times better than those of the BDS-2 based solution. It makes sense to use all the BDS satellites for ERP estimation in particular if the a priori box-wing model is not considered. Also, interestingly, LOD results from the ECOMA and bw+ECOMA models show a more than two times higher scatter than the solutions from other SRP models, despite the fact that we observe better orbit solutions in this case. This is because LOD estimation is sensitive to the uncertainty in orbit modeling. High correlations between the along-track acceleration and the LOD parameter contaminate the LOD estimates. As the computation time period is only two years, the draconitic year related signals cannot yet be extracted.
Geocenter (GCC) motion is the movement of the center of mass of the total earth system with respect to the center of figure of the solid earth surface, in response to various geophysical fluid displacements within the earth system [62], [64]. The origin of ITRF 2014 is defined solely by SLR data collected from the two LAGEOS satellites over more than 20 years. The reason is that SLR measurements are absolute and unambiguous, with an accuracy of sub-centimeter. Nongravitational forces on LAGEOS satellites may be accurately modeled at sub-centimeter level without solving for any model parameter. In this case, less parameters need to be estimated and the orbit modeling errors will not affect the geocenter coordinates. Moreover, with small force errors, long arcs (for instance 7-day-arc) can be used to further de-correlate parameters of different types [65].   GNSS technique does not have these advantages and it is well-known that geocenter coordinates are sensitive to radiation pressure modeling deficiencies, especially for the Z component. Fig. 14 and Fig. 15 show the time series and amplitudes of the BDS-2, BDS-3 and BDSALL based geocenter coordinates in the Z direction for the ECOM model with and without the box-wing model. The difference between BDS-2 and BDS-3 based solutions is very clear, the former case shows much more scattered estimates. The difference between BDS-3 and BDSALL is however minor. The use of the a priori box-wing model mitigates a great majority of the periodical signals.   Fig. 11. Time series of X pole coordinates with respect to IERS-C04-14 (mas) with (red color) and without (blue color) the box-wing model. 5, 7 cpy and the sum of these amplitudes for the BDSALL solution with all the SRP models. It is clear that the annual signal reduces from about 28 mm -60 mm to about 6 mm. The sum of all the amplitudes is reduced by more than a factor of two for all the ECOM models when using the box-wing model.

VI. CONCLUSION
BDS satellites have completed the full constellation since June 2020. The optimal SRP model and the contribution of BDS satellites to the geodetic parameters are not yet fully clear. BDS satellite metadata information is officially published at the end of 2019. However, only one of the three optical properties is given. Nevertheless, with the published   Fig. 13. Time series of LOD with respect to IERS-C04-14 (ms) with (red color) and without (blue color) the box-wing model.  metadata information, we can estimate optical parameters based on true BDS measurements. Various physical effects, such as yaw bias, self-shadowing, radiator emission and thermal radiation of solar panels are considered in the adjustment. We find that to obtain more physical optical parameters, it is necessary to consider the T-shape body of BDS-3 CAST MEO satellites and to take the self-shadowing effects into account. Also, the BDS-3 CAST MEO satellites should be divided into several sub-groups considering the optical property estimates. Satellites C32 and C33 seem to have an additional earth-facing surface, which may host the SAR antennas.
Thermal radiation of solar panels for BDS-3 satellites is more important than the radiator acceleration. Because BDS-3 satellites keep maintaining yaw-steering attitude, ECOM/ECOM2 parameters are able to absorb the constant radiator acceleration in the -X direction even inside shadows. However, thermal radiation of solar panels performs differently and needs to be well modeled. We have shown that an additional along-track acceleration absorbs part of this effect and thus results in better orbit solution than the ECOM-only solution.
All the physical effects and the estimated optical properties are considered in a physical box-wing model, and is jointly used with different empirical models in satellite orbit determination. The improvement by using the a priori boxwing model is obvious for all the empirical models. The best solution of orbit misclosures is 4.9 cm, 2.1 cm and 2.1 cm for BDS-2, BDS-3 CAST and BDS-3 SECM satellites respectively. The best 24-hour orbit predictions is 13.8 cm, 8.5 cm and 8.8 cm for BDS-2, BDS-3 CAST and BDS-3 SECM satellites respectively. Dependencies of SLR residuals of BDS-3 satellites on the elongation angle almost vanish by using the box-wing model. The smallest STD value of SLR residuals is 3.7 cm for both BDS-3 CAST and SECM satellites.
The geodetic parameters are sensitive to different constellations and different SRP models. Polar motion parameters estimated with BDS-3 satellites are two times closer to IERS-C04-14 than those determined with BDS-2 satellites. Results using all the BDS satellites show a slight improvement over the BDS-3 based results. Therefore, it makes sense to process BDS-2 and BDS-3 satellites together for the geodetic applications. The use of the box-wing model improves geodetic parameters in all cases. The results based on the ECOM model with additionally estimated along-track acceleration shows much worse LOD estimates than the other models because the along-track acceleration is highly correlated with LOD and other ECOM parameters. The smallest STD value compared to IERS-C04-14 is 72.4 µas 71.0 µas and 36.3 µs for X pole, Y pole and LOD parameters respectively. The sensitivity of geocenter coordinates on different satellites and different SPR models is evident. BDS-3 satellites show clear better performance than BDS-2 satellites. The total amplitude has been reduced by more than a factor of three on average for all the empirical models after the consideration of the physical a priori model.