Improving the Orbits of the BDS-2 IGSO and MEO Satellites with Compensating Thermal Radiation Pressure Parameters

Abstract

The orbit accuracy of the navigation satellites relies on the accurate knowledge of the forces on the spacecraft, in particular the non-conservative perturbations. This study focuses on the Inclined Geosynchronous Orbit (IGSO) and Medium Earth Orbit (MEO) satellites of the regional Chinese BeiDou Navigation Satellite System (BDS-2), for which apparent deficiencies of non-conservative models are identified and evidenced in the Satellite Laser Ranging (SLR) residuals. The orbit errors derived from the empirical 5-parameter Extended CODE Orbit Model (ECOM) as well as a semi-analytical adjustable box-wing model show prominent dependency on the Sun elongation angle, even in the yaw-steering attitude mode. Hence, a periodic acceleration in the normal direction of the +X surface, presumably generated by the mismodeled thermal radiation pressure, is introduced. The SLR validations reveal that the Sun elongation angle-dependent systematic errors were significantly reduced, and the orbit accuracy was improved by 10–30% to approximately 4.5 cm and 3.0 cm for the BDS-2 IGSO and MEO satellites, respectively.

1. Introduction

With the launch of the latest Geostationary Earth Orbit (GEO) satellite of the China BeiDou Navigation Satellite System (BDS), BDS started to provide global positioning, navigation, and timing (PNT) services with 57 active satellites at the beginning of July, 2020 [1]. Among them, there are five GEO, seven Inclined Geosynchronous Orbit (IGSO), and three Medium Earth Orbit (MEO) satellites that form the regional BDS constellation (BDS-2).

For the precise orbit determination (POD) of BDS-2 satellites, much research was carried out to improve the GEO orbits with enhanced solar radiation pressure (SRP) models [2,3], observation correction models [4,5,6] or onboard observations from low-earth orbits [7] and MEO satellites [8]. With the adoption of the empirically derived SRP model for BDS GEO satellites, the radial orbit accuracy reached 10 cm as validated by Satellite Laser Ranging (SLR) [3]. For IGSO and MEO satellites with yaw-steering (YS) and orbital normal (ON) attitude mode, much more attention was paid on the yaw attitude modeling [9,10] as well as dynamic force modeling for satellites in ON mode [11,12,13].

Although significant improvement was obtained for satellites in ON mode with the improved SRP model, the orbit accuracy is still worse by a factor of two compared with that in the YS mode [14]. Moreover, the SLR residuals show the Sun elongation angle (i.e., the $ϵ$ angle between the Sun and Earth as seen from the satellite)-dependent errors related to the purely empirical 5-parameter Extended CODE Orbit model (referred to as ECOM, [15]) for BDS-2 IGSO satellites [11]. The systematic orbit errors were also identified in the International GNSS Service (IGS, [16])-combined Multi-GNSS orbit solutions [17]. The radial orbit accuracy relying on the dynamic model is essential for GNSS solutions and scientific applications [18,19,20]. In the absence of a physics-based SRP for the BDS-2 IGSO and MEOs, the studies have typically used an empirical SRP model. Fortunately, usage of analytical or semi-analytical models for BDS is possible since the satellite metadata were released by the China Satellite Navigation Office (CSNO) [21], which comprise the mass, size, and optical properties of the main satellite surfaces and can be used to improve the dynamic orbit modeling of BDS-2 satellites.

The objective of this research is to diagnose the Sun elongation angle-dependent errors for BDS-2 IGSO and MEO satellites in YS mode and to seek a possible solution for their mitigation. Following the introduction of the theory for solar radiation and thermal radiation pressure (TRP) in Section 2, the IGS Multi-GNSS (MGEX) orbits for BDS-2 IGSO and MEO from different analysis centers (ACs) are investigated to demonstrate the noticeable errors. Moreover, the standard adjustable box-wing (ABW) model is used to investigate whether they can be reduced. Unfortunately, this mismodelling deficiency cannot be compensated by the ABW model as shown in Section 3. Through analysis of the series of ECOM parameters as well as the distribution of SLR residuals on the orbital angle and the Sun elevation angle, the thermal re-radiation of the satellite bus is considered the potential cause. A periodic acceleration along the +X surface of the satellite body frame obeying the IGS convention [22] is proposed and calibrated with the real observations in Section 4. The derived parameters for each BDS-2 IGSO and MEO satellite are then used as the a prior model to augment ECOM. Furthermore, the performance of the proposed model is assessed with the metrics of SLR residuals and orbit prediction performance in Section 5. Finally, the study is summarized.

2. Orbit Dynamic Perturbations

For navigation satellites at an altitude over tens of thousands of kilometers, the orbit accuracy relies on the accurate knowledge of the forces acting on the spacecraft surfaces. Among them, non-gravitational forces require precise modeling, as they are more easily affected by factors such as incident radiation fluxes, the satellite attitude, the optical as well as thermal properties of the satellite bus, and solar panels. These non-gravitational forces mainly include SRP, TRP, earth radiation pressure, and satellite antenna thrust [22].

2.1. Solar Radiation Pressure

This subsection provides some classical and widely used SRP models for POD within the IGS community, including empirical and semi-analytical models.

2.1.1. Extended CODE Orbit Model (ECOM)

The ECOM model, widely used within the IGS community for modeling SRP acting on GNSS satellites, decomposes the perturbations in a Sun-orientation reference frame with the three orthogonal directions, i.e., ${\stackrel{\rightharpoonup }{e}}_D$, ${\stackrel{\rightharpoonup }{e}}_Y$, and ${\stackrel{\rightharpoonup }{e}}_B$, expressed as

$$\begin{array}{c}{\stackrel{\rightharpoonup }{e}}_D=\frac{{\overrightarrow{r}}_{\odot }-\overrightarrow{r}}{\Vert {\overrightarrow{r}}_{\odot }-\overrightarrow{r}\Vert }\\ {\stackrel{\rightharpoonup }{e}}_Y=-\frac{\overrightarrow{r}\times {\overrightarrow{e}}_D}{\Vert \overrightarrow{r}\times {\overrightarrow{e}}_D\Vert }\\ {\stackrel{\rightharpoonup }{e}}_B={\stackrel{\rightharpoonup }{e}}_D\times {\stackrel{\rightharpoonup }{e}}_Y\end{array}$$

(1)

where ${\overrightarrow{r}}_{\odot }$ and $\overrightarrow{r}$ are the geocentric vectors of the Sun and the satellite, respectively. The vector ${\overrightarrow{e}}_D$ is the unit vector in the satellite–Sun direction; ${\overrightarrow{e}}_Y$ points along the satellite’s solar panel axis and coincides with the Y axis of the satellite reference frame in the YS mode. ${\overrightarrow{e}}_B$ completes the right-hand orthogonal system. With this frame, the SRP accelerations are expressed as Fourier series using the satellite’s argument of latitude $u$as an angular argument,

$$\begin{array}{c}D\left(u\right)=D_0+{{\displaystyle \sum }}_{i=1}^n\left[D_{i,c}\cos \left(iu\right)+D_{i,s}\sin \left(iu\right)\right]\\ Y\left(u\right)=Y_0+{{\displaystyle \sum }}_{i=1}^n[Y_{i,c}\cos \left(iu\right)+Y_{i,s}\sin \left(iu\right)]\\ B\left(u\right)=B_0+{{\displaystyle \sum }}_{i=1}^n[B_{i,c}\cos \left(iu\right)+B_{i,s}\sin \left(iu\right)]\end{array}$$

(2)

where $i$ is an index, n is the terminated order and $D_{i,c}$, $D_{i,s}$, $Y_{i,c}$, $Y_{i,s}$, $B_{i,c}$, and $B_{i,s}$ are the amplitudes of the periodic signals. For the original ECOM model, n was selected as 1, i.e., three constant and six periodic parameters had to be estimated [15]. Usually, the reduced ECOM model [23], i.e., the classical ECOM, is used, and in the D and Y directions, only the constant parameters are estimated compared to the original ECOM as follows:

$$\begin{array}{c}D\left(u\right)=D_0\\ Y\left(u\right)=Y_0\\ B\left(u\right)=B_0+B_{1,c}\cos \left(u\right)+B_{1,s}\sin \left(u\right)\end{array}$$

(3)

Moreover, an a priori SRP model can be incorporated with the ECOM model, e.g., the ROCK model [24].

2.1.2. The Adjustable Box-Wing (ABW) Model

The ECOM models have deficiencies for SRP modeling in ON mode and cause noticeable errors in the GPS-derived geodetic parameters, i.e., station coordinates, earth rotation parameters [25,26]. In order to overcome this deficiency, the adjustable box-wing (ABW) model was proposed based on the physical interaction between solar photons and the satellite surface by approximating the satellite surface structure as a box (bus) with two wings (solar panels) in YS mode [27]. As we have modified this ABW model with compensated parameters discussed in Section 4, we refer to it as a standard ABW model.

The SRP force acting on a surface with area $\mathrm{A}$ depends on the relative alignment of the Sun’s direction (${\overrightarrow{e}}_{\odot }$) and the surface normal ${\overrightarrow{e}}_N$ along with the incident angle ($\theta$), the fraction of absorbed photons ($\mathsf{\alpha }$), of diffusely reflected photons ($\mathsf{\delta }$), and of specularly reflected photons ($\rho$). For an illuminated surface, the acceleration ${\overrightarrow{a }}_{SRP}$ at the distance of 1 astronomical unit (AU) can be obtained as follows:

$${\overrightarrow{a}}_{SRP}=-\frac{A}{M}\frac{S_0}{c}\cos \theta \left[\left(\mathsf{\alpha }+\mathsf{\delta }\right){\overrightarrow{e}}_{\odot }+2\left(\frac{\delta }{3}+\rho \cos \theta \right){\overrightarrow{e}}_N\right],$$

(4)

where $M$ is the mass of the satellite, $S_0$ is the solar flux at 1 AU; a revised value for $S_0$ based on re-analysis and re-calibration of satellite data of 1361 W/m² was used [28], $\mathrm{c}$ is the velocity of light. For materials that have zero thermal capacity and completely prevent heat transfer toward the satellite interior such as multilayer insulation (MLI), the energy absorbed by the satellite surfaces is instantaneously reradiated back to space, and the corresponding thermal re-radiation (TRR) acceleration ${\overrightarrow{a}}_{TRR}$ can be obtained approximately according to Lambert’s law as follows [29]:

$${\overrightarrow{a}}_{TRR}=-\frac{A}{M}\frac{S_0}{c}\cos \theta \frac{2}{3}\alpha {\overrightarrow{e}}_N$$

(5)

For surfaces of the satellite bus covered by MLI, the perturbations generated by SRP and TRR are the sum of Equation (4) and Equation (5):

$$\overrightarrow{a}=-\frac{A}{M}\frac{S_0}{c}\cos \theta \left[\left(\mathsf{\alpha }+\mathsf{\delta }\right)\left({\overrightarrow{e}}_{\odot }+\frac{2}{3}{\overrightarrow{e}}_N\right)+2\rho \cos \theta {\overrightarrow{e}}_N\right]$$

(6)

Note that using this equation for all the surfaces of the satellite bus means assuming the satellite box is completely covered by MLI. For solar panels, no additional TRR is considered, and Equation (4) is used to model the perturbation of SRP.

As the accuracy of the analytical model is limited by the accuracy of the satellite attitude, geometry, optical, and thermal properties of the satellite bus as well as solar panels, some parameters were fitted with the tracking observation, including the solar panel lag angle (SB), Y bias (Y0), a scaling factor of solar panels(SP), the absorption plus diffusion (AD) as well as the reflection (R) coefficients of the illuminated +X/+Z/−Z bus surface (+XAD, +ZAD and −ZAD; +XR, +ZR and −ZR), to improve the POD accuracy. The partial derivatives of the acceleration with respect to these parameters can be found in [27].

2.2. Thermal Radiation Pressure

The thermal radiation forces of a satellite mainly arise from the interaction of the spacecraft with the environment. Usually, the Sun is the primary source of radiation, and the internal subsystems are the other sources of heat within the satellite bus. In this study, the thermal force caused by the re-radiation of the absorbed energy from the Sun in the form of heat is indicated as TRR, whereas the force due to emitted heat from the thermal radiator is called thermal radiation (TR).

Since TRR and SRP show similar characteristics that originated from the Sun, the former is often incorporated with SRP in a combined form as done in Equation (6). Ref. [24] presented the SRP acceleration approximately represented as a Fourier series and stated that the TRR force was about 5% of the total solar radiation. Moreover, a detailed thermal force model was derived at the University College London (UCL) based on the analytical approach [29,30]. For solar panels, the acceleration ${\overrightarrow{a}}_{TRR SP}$ due to the thermal emission of radiation can be computed as follows:

$${\overrightarrow{a}}_{TRR\_SP}=-\frac{2A\sigma \left({\epsilon }_f{T}_f^4-{\epsilon }_b{T}_b^4\right)}{3Mc}{\overrightarrow{e}}_N$$

(7)

where $\sigma$ is the Stephan–Boltzmann constant (5.6699 × 10⁻⁹ Wm⁻² K⁻⁴), ${\epsilon }_f$ and ${\epsilon }_b$ are the emissivity of the front $f$ and back $b$ panels, and $T_f$ and $T_b$ are the absolute temperature in Kelvin of the surfaces. Usually, the same emissivity as well as temperature are assumed for the front and back surfaces of solar panels, resulting in no thermal perturbations for solar panels. On the other hand, for the satellite bus covered by MLIs, the TRR acceleration ${\overrightarrow{a}}_{TRR\_BUS}$ can be computed as follows:

$${\overrightarrow{a}}_{TRR\_BUS}=-\frac{2A\epsilon \sigma T_{MLI}^4}{3Mc}{\overrightarrow{e}}_N$$

(8)

where $T_{MLI}$ is the absolute temperature in Kelvin of the MLIs, and $T_{MLI}^4$ can be represented as follows:

$$T_{MLI}^4=\frac{\alpha S_0\cos \theta +{\epsilon }_{eff}\sigma T_{sat}^4}{\sigma \left({\epsilon }_{MLI}+{\epsilon }_{eff}\right)}$$

(9)

where ${\mathsf{\epsilon }}_{eff}$ is the effective emissivity, ${\mathsf{\epsilon }}_{MLI}$ is the emissivity of the MLI, $\alpha$ is the absorptivity of the MLI, and $T_{sat}$ is the temperature of the satellite below the MLI. According to model parameter sensitivity analysis for GPS IIR and Jason-1 satellite in Adhya (2005), the TRR acceleration is not highly dependent on ${\mathsf{\epsilon }}_{MLI}$, as a higher value of this parameter acts to decrease the temperature and reduce the force but also to increase the amount of the energy emitted and hence increase the force.

Usually, the heat generated by satellite internal subsystems is vented to space via radiators or louvers on the surface of the bus. According to [31], this emission will also produce a recoil acceleration as follows:

$${\overrightarrow{a}}_{TR}=-\frac{2q_{dissipated}}{3Mc}{\overrightarrow{e}}_N$$

(10)

where $q_{dissipated}$ is the power that is generated by the payload but dissipated to space. Hence, the thermal radiator details for a certain satellite are required for TR modeling. However, this kind of information is not disclosed by spacecraft manufacturers. In general, the radiators are potentially mounted on the unilluminated surfaces, i.e., +Y, −Y, and −X according to the IGS attitude conventions [22]. Recently, TR force due to emitted heat from the thermal radiator caused by the onboard atomic clock generating the orbit errors in eclipse seasons for Galileo satellites was identified, and an empirical model for accelerations caused by the spacecraft’s thermal radiation was developed in [31]. By incorporating the model in the CODE (Center for Orbit Determination in Europe) MGEX solution, an improvement of Galileo orbits by about 14% during eclipse seasons was achieved [26]. Moreover, based on the analysis of the constant Y parameter (Y0) of the ECOM model, the potential radiators on the −X surface were identified for GLONASS satellites [32]. With this parameter as well as the adjusted box-wing parameters, the orbit difference between two consecutive arcs for both GLONASS-M and GLONASS-K satellites was reduced. In particular, the improvement for GLONASS-M satellites was greater than 30% for the ECOM model during eclipse seasons.

3. Status of BDS-2 IGSO and MEO Orbits

In this section, the orbits of the BDS-2 IGSO and MEO satellites from three IGS MGEX ACs [22,26,33,34], i.e., CODE, GeoForschungsZentrum, and Wuhan University (WHU), are evaluated with SLR to show the Sun elongation angle-dependent systematic errors. Afterwards, the POD strategies are presented including solutions with the four SRP models.

3.1. BDS-2 Orbits from MGEX ACs

The BDS-2 IGSO and MEO orbit solutions from three MGEX ACs, i.e., CODE, GFZ, and WHU, over the full year period of 2018 were collected and validated with SLR. Figure 1 shows the SLR residuals of the C13 satellite with respect to the Sun elongation angle ($ϵ$) for each AC. The empirical ECOM SRP model without a prior model was employed for BDS-2 IGSO and MEO data analysis by all the three ACs. Although all of the BDS-2 satellites are equipped with Laser Ranging Array (LRA), only the C08, C10, C11, and C13 satellite are tracked by the International Laser Ranging Service (ILRS, [35]). Note that SLR residuals exceeding an absolute value of 50 cm were excluded in this study. As can been seen from Figure 1, although the scatter distribution of the SLR residuals for the three ACs’ solutions differ slightly, they all show pronounced systematic orbit errors dependent on the $ϵ$ angle, which indicates the modeling deficiencies of non-conservative perturbations. We checked the three other BDS-2 IGSO and MEO satellites tracked by ILRS, and a similar pattern was identified. One might argue that these $ϵ$-angle-dependent errors are caused by using wrong geometric corrections of LRA coordinates. However, the simulated results indicate that a horizontal LRA coordinate error of up to 30 cm will only lead to a deviation of 1 cm along the radial direction, and a vertical LRA coordinate error will not lead to discrepancies in the function of the $ϵ$ angle but a bias in the SLR residuals. For further explanation, we list the statistical values for each satellite in

Table 1. Slope of the SLR residuals with respect to the Sun elongation angle ($ϵ$ ) for BDS-2 IGSO and MEO satellites of three MGEX ACs. (unit: cm·(100 deg)⁻¹).

Solutions	C08	C10	C13	C11
CODE	−7.0	−13.6	−16.8	−2.3
GFZ	−9.1	−9.6	−15.2	−5.4
WHU	−6.3	−6.1	−13.9	−3.2

, where the slope of the SLR residuals with respect to the Sun elongation is also listed.

According to Table 1, all SLR residuals showed a systematic negative slope ranging from −2.3 to −16.8 cm/100 deg with respect to the Sun elongation angle. This indicates that BDS-2 IGSO and MEO have modeling deficiencies of orbit dynamic perturbations. Despite that the slope value for the same satellite from the three ACs is different, the MEO C11 satellite exhibits the smallest slope value compared to the other three IGSO satellites, i.e., C08, C10, and C13. In addition, C13 had the largest slope value among the four satellites, and the slope value of C13 was almost twice that of C08 for the same AC, although C13 and C08 are in the same orbit plane and were manufactured based on the same satellite platform (Launched in 2011 and 2016, respectively. See http://www.csno-tarc.cn/en/performance/track, accessed on 15 November 2021). This fact suggests that there might be a force that is inaccurately modeled or not considered in current force modelling.

3.2. Strategy for the POD with IGS and iGMAS Observations

In this study, tracking data from the IGS Multi-GNSS Experiment (MGEX) performed by the International GNSS Service [16] and the International GNSS Monitoring and Assessment System (iGMAS, [36]) in 2018 were used to determine orbit solutions using B1 and B2 signals with a sampling interval of 300 s. A POD arc length of 24-h and 10-degree cutoff elevation and elevation-dependent weighting for the observations under 30 degrees were applied. Figure 2 shows the distribution of the approximately 100 stations, among which about 20 are from iGMAS to further fulfill the ground coverage in China’s mainland. Almost the same POD strategy was used as for the Wuhan University MGEX products [34], except for the SRP model. All data were processed by a modified version of Position And Navigation Data Analyst (PANDA) software [37] using a two-step approach. For observation corrections that depend on the satellite attitude, such as phase center corrections of satellites, phase wind-up corrections, and LRA offsets, the YS-ON along with the continuous YS model was used [9,10,11]. For the modeling aspect of antenna thrust for BeiDou IGSO and MEO satellites, measured values from [38] were adopted.

Table 2. Summary of the models and information used in this study.

Models	Description
Geopotential (static)	EGM2008 with 12 × 12
N-body	Sun, Moon, Jupiter, Venus, Mars, Mercury, Uranus, Neptune, Saturn, Pluto, Charon as point mass JPL DE405 ephemeris used
Solid Earth tides	IERS conventions 2010 [39]
Solid Earth pole tides	IERS conventions 2010 [39]
Ocean tides	None
Ocean pole tides	None
Relativistic effects	IERS conventions 2010 [39]
Solar radiation pressure	See Table 3
Antenna thrust	Applied. [38]
Earth radiation pressure	Applied. [40]

summarizes the background geophysical models used, and

Table 3. Solutions with the SRP model used and estimated parameters.

Solution	SRP Model	Estimated Parameters
ECOM	The 5-parameter ECOM model	$D_0$$, Y_0$$, B_0$$, B_c$$\mathrm{and}{B}_s$.
ABW	The standard ABW	$SP$$, SB$$, Y_0$$,+XAD$$,-XAD$$,+XR$$,-XR$$,+ZAD$$,-ZAD$$,+ZR$$\mathrm{and}-ZR$.
ABW + TRR	The modified ABW with TRR parameter	$SP$$, SB$$, Y_0$$,+XAD$$,-XAD$$,+XR$$,-XR$$,+ZAD$$,-ZAD$$,+ZR$$,-ZR$$\mathrm{and}k$.
ECOM + TRR	ECOM with additional TRR parameter as the a prior	$D_0$$, Y_0$$, B_0$$, B_c$$\mathrm{and} B_s$.

k indicates the compensated TRR parameter in Equation (13).

lists the four orbit solutions along with the estimated parameters.

3.3. The Orbits Determined with the Standard ABW Model

The standard ABW model can accurately compensate the perturbing acceleration acting on satellites. By estimating the optical properties of the satellite’s surface along with other parameters (see Table 3), this model can provide a similar orbit performance as that generated with the ECOM for satellites in YS attitude mode [27]. ABW can be normally used for identifying the deficiencies in SRP modeling and for providing a precise orbit solution based on the physical interaction between solar photons and satellite surfaces, as done by [11] and [32]. While this requires proper knowledge of the optical, physical, thermal, and geometrical properties of the satellite surfaces, these values were released by [21] and can be used as the a priori values in the ABW model for BDS-2 satellites. However, CSNO did not release the specular and diffuse reflection coefficients; therefore, values from [11] were adopted. Although the sum of absorbed, diffusely reflected, and specularly reflected coefficients should theoretically equal one, we did not implement this constraint, which allowed absorption of the unmodeled perturbations as much as possible. A loose constraint of 0.1 and a tight constraint of 0.01 were adopted, as applied in [27], for the AD and R parameters of the illuminated surfaces, respectively.

Figure 3 shows the SLR residuals with respect to the $ϵ$ angle for C08 and C13. These two IGSO satellites were selected and compared because they are in the same orbital plane and were manufactured based on the same satellite platform. Hence, they have almost identical observation conditions from ground networks. It can be observed clearly that the orbit solution determined by the ABW model also showed pronounced $ϵ$-angle-dependent patterns, and the slope of SLR residuals with respect to $ϵ$ for C13 was greater than that of C08, which is similar to ECOM. This was not expected, as the use of the standard ABW model can usually reduce the periodic β- and μ-angle-dependent variations in SLR residuals, where β is the Sun elevation angle above the orbital plane, and the orbital angle μ is the argument of the latitude of the satellite with respect to midnight. The results of these two satellites determined by the standard ABW model obviously showed that the $ϵ$-angle-dependent radial orbit errors did not originate from ground observations, in addition to mismodeled SRP perturbation, as they are in the same orbital plane, and the geometry as well as perturbation should be similar.

To reveal the characteristics of the orbit errors, the SLR residuals of C13 with respect to the β and μ angles are shown in Figure 4. It can be observed that apparent β- and μ-angle-dependent errors exist. The maximum positive bias occurred around the midnight point (μ approached zero), while the largest negative bias was near the noon point (μ was close to 180 deg). This indicates that there are some mismodeled or unmodeled non-gravitational perturbations with similar characteristics as SRP. TRR may be the candidate as it mainly originated from the Sun and was not fully captured in the ABW or empirical SRP model.

4. The Compensated TRR Model for BDS-2 IGSOs and MEOs

In this section, a periodic acceleration along the +X surface to compensate the TRR is proposed and calibrated with the real observations based on the standard ABW model.

4.1. TR

As aforementioned, we divided the thermal forces into TRR and TR. The former is caused by the re-radiation of the absorbed energy from the Sun in the form of heat, whereas the latter is caused by the emitted heat from a thermal radiator. The radiators are used to release heat, and in most cases, are located on the −X or ±Y surfaces as these are not illuminated in YS mode. To evaluate whether the orbit deficiency in Figure 4 was caused by TR, the orbits were re-determined with the same strategy as [34] based on the ECOM SRP model. Figure 5 illustrates the variations in the estimates of five SRP parameters (see Equation (3)), i.e., three constant parameters along the D, Y, and B axes (D0, Y0, B0) and two periodic parameters along B axes (Bc and Bs) for satellite C13. It can be clearly seen that all parameters did not show noticeable changes in and out of the eclipse seasons except for Bc, for which an apparent bias up to 1–2 nm/s² can be identified in the eclipse seasons. Most important, the estimates of the Y0 parameter were almost constant with a magnitude less than 0.5 nm/s². Therefore, radiator effects on ±Y surfaces are either balanced or are very small.

Furthermore, the possibility of TR along the X axis was also analyzed. For GLONASS and Galileo, one additional constant radiator parameter along the -X direction was introduced to account for the potential TR force generated by the radiator along the X axis [31,32]. Similarly, the constant parameter was also estimated for BDS-2 IGSO and MEO satellites in this study. However, the SLR residuals of C13 still showed noticeable $ϵ$-angle-dependent systematic values as illustrated in Figure 6. Hence, the constant TR acceleration along the X axis cannot explain the orbit deficiency. Based on the results, we conclude that TR perturbation is not the root of the orbit errors.

4.2. Compensated TRR Model

Usually, thermal perturbations for solar panels can be absorbed with the parameter along the Y axis in a satellite body frame. Hence, only the TRR for the surfaces of the satellite-bus should be considered. For satellites in YS mode, the +X, +Z, and −Z surfaces are illuminated. As the perturbations along the Z axis correspond to the radial direction, the mismodeling perturbation along the Z axis can be compensated mostly by the satellite clock parameters in a zero-difference network solution. Hence, we focus on the mismodeling perturbations on the +X surface.

Figure 7 shows the distribution of the TRR acceleration derived from Equation (5) with respect to the Sun elevation angle β and orbital angle μ for C13 in the satellite–Sun direction (D). Generally, the TRR perturbation along the D direction varies periodically with both the μ and β angles. Hence, the constant term along the D direction in ECOM cannot absorb this perturbation completely once the TRR becomes significant, and the remaining unmodeled perturbation will result in deficiencies in orbit determination. The magnitude of the TRR force is dependent on the geometrical and physical properties of the satellite surface, which are usually slightly different among the individual satellites. This probably explains the orbit performance differences between C08 and C13. More importantly, we can observe that the orbit errors revealed by SLR residuals in Figure 4 show a similar pattern as the TRR-induced accelerations in Figure 7.

With Equations (8) and (9), the TRR acting on the satellite bus can be divided into two parts:

$${\overrightarrow{a}}_{TRR\_BUS}=-\frac{2A\epsilon }{3Mc\left({\epsilon }_{MLI}+{\epsilon }_{eff}\right)}\left(\epsilon \sigma T_{sat}^4+\alpha S_0\cos \theta \right){\overrightarrow{e}}_N$$

(11)

The first part is related to the 4th power of the temperature under the MLIs, whereas the other is proportional to the cosine of the incident angle. The magnitude of the first part is less than 1% of the TRR perturbation [30]. Hence, it can be omitted without significant loss of accuracy, leaving the remaining part as follows:

$${\overrightarrow{a}}_{TRR\_BUS}=-\frac{2A\epsilon }{3Mc\left({\epsilon }_{MLI}+{\epsilon }_{eff}\right)}\alpha S_0\cos \theta {\overrightarrow{e}}_N=-k\cos \theta {\overrightarrow{e}}_N$$

(12)

with $k$=$\frac{2A\epsilon \alpha S_0}{3Mc\left({\epsilon }_{MLI}+{\epsilon }_{eff}\right)}$. By adding Equations (6) and (12), the SRP and TRR perturbation acting on the +X surface can be obtained:

$$\overrightarrow{a}=-\frac{A}{M}\frac{S_0}{c}\cos \theta \left[\left(\mathsf{\alpha }+\mathsf{\delta }\right)\left({\overrightarrow{e}}_{\odot }+\frac{2}{3}{\overrightarrow{e}}_N\right)+2\rho \cos \theta {\overrightarrow{e}}_N\right]-k\cos \theta {\overrightarrow{e}}_N$$

(13)

Therefore, the compensated TRR scale factor $k$ on the +X surface can be estimated with the parameters of the ABW model by fitting the observations. To obtain a reliable estimation of the thermal scale factor $k$, the correlations with the orbital as well as ABW parameters were investigated. Figure 8 shows the correlations among those parameters for C13 on Day of Year (DOY) 200, 2018. Generally, the correlations of the thermal scale factor $k$ with the ABW parameters were less than 0.5 except for those with the SP parameter, with values up to 0.7. For the sake of reliable estimation, an a priori constraint of 5 nm/s² was applied for the thermal scale factor $k$.

Based on equation (13), the daily thermal scale factors for the BDS-2 IGSO and MEO satellites were fitted with the real observations in 2018, and Figure 9 demonstrates the estimates for C13 as an example. Noticeable positive and nearly constant values were identified for the estimates out of the eclipse season, whereas the estimated value of k dropped to −2 nm/s² when satellite crossed the eclipse regime. Other IGSO and MEO satellites showed a similar behavior. The averaged estimates derived from the daily estimates outside the eclipse regime for each satellite are listed in

Table 4. The averaged estimates of the thermal scale factor for the BDS-2 IGSO and MEO satellites. (unit: nm/s²).

IGSOs	C06	C07	C08	C09	C10	C13
Value	1.5	1.6	1.2	2.3	1.8	2.6
MEOs	C11	C12	C14
Value	1.6	1.5	0.9

. It can be observed that the estimates did not show clear dependency on the satellite type, i.e., IGSO and MEO. In general, the averaged estimates varied from about 0.9 to 2.6 nm/s². C13 had the largest estimate of approximately 2.6 mm/s², and the smallest was about 0.9 nm/s² for C14.

5. Validation

Based on the analysis above, the parameters in Table 4 were used as a prior model to augment ECOM for the BDS-2 IGSO and MEO satellites. This morel introduces a periodic acceleration along the +X axis to compensate possible model deficiencies in orbit determination. For the purpose of comparison and validation, all four solutions were determined with same data set. The SLR residuals and predicted orbits were used as the metrics for validation.

Figure 10 illustrates the one-way SLR residuals of the BDS-2 IGSO and MEO satellites from the four solutions listed in Table 3. It can be clearly observed that the orbit solution of ECOM augmented by TRR had the best performance. The SLR residuals resembled a normal distribution, and many were located around zero and were distributed equally around zero. In addition, we observed distinct improvements by introducing the TRR parameter with respect to the counterpart solutions. The solutions based on the semi-analysis method, i.e., ABW and ABW + TRR, showed a scatter distribution of the SLR residuals compared to the ECOM or ECOM + TRR solutions. Moreover, this was more evident for the IGSOs than the MEOs, as the higher orbit height of the IGSO satellites results in a smaller nadir angle, leading to a stronger correction between the estimated parameters. The statistics, i.e., the bias and Root Mean Square (RMS), are listed in

Table 5. The SLR residual statistics for the four orbit solutions of C08, C10, C11, and C13. (unit: cm).

	ECOM		ABW		ABW + TRR		ECOM + TRR
	Mean	RMS	Mean	RMS	Mean	RMS	Mean	RMS
C08	1.7	5.7	4.0	7.8	1.3	6.1	−0.1	5.2
C10	3.2	5.6	3.2	7.8	2.0	7.5	0.1	4.4
C11	0.0	3.4	0.1	3.9	0.0	3.2	0.0	2.9
C13	1.9	6.4	2.1	8.9	0.1	6.6	0.1	4.3

. The ABW + TRR again had lower performance compared to the solutions using the empirical SRP model, i.e., ECOM/ECOM + TRR. When applying the TRR parameter presented in Table 4, the radial orbits of C13 were improved by approximately 30% from about 6 cm to 4 cm. For the three other BDS-2 satellites, the improvement was not as significant as C13, but it still reached 10% to 20%. In order to further clarify this, the SLR residuals of ECOM + TRR solution are presented with respect to the Sun elongation angle in Figure 11. As can be seen, in contrast to the results shown in Figure 1, the slope of the SLR residuals as a function of the $ϵ$ angle almost vanished, which contributed to the improvements revealed by the SLR measurements.

Besides the SLR measurements, the orbit prediction was used to evaluate the model. In this study, the 24-h predicted orbits were generated for all four solutions. Then, the accuracy of the predicted orbits was obtained and compared with the orbits of the next day. It should be mentioned that the predicted orbit is affected by accuracy of the Earth rotation parameter, which was fixed in this study. Hence, the predicted orbits in this way showed somewhat better performance compared with the actual condition.

Table 6. The RMS of the 24-h predicted orbits for the four orbit solutions. A represents the along-track, C the cross-track, and R the radial direction. (unit: cm).

		ECOM			ABW			ABW + TRR			ECOM + TRR
		A	C	R	A	C	R	A	C	R	A	C	R
IGSO	C06	127.5	8.0	29.6	158.6	9.5	35.5	137.5	9.0	30.7	113.6	7.6	28.0
	C07	84.5	10.1	25.5	136.8	10.9	30.2	123.5	9.7	25.4	95.6	10.3	25.7
	C08	126.0	8.7	31.3	156.7	10.8	33.2	153.1	10.7	37.3	134.3	10.0	25.7
	C09	95.8	8.8	25.1	222.4	8.9	40.0	175.3	9.0	36.6	108.9	9.2	24.5
	C10	85.9	10.6	19.4	104.1	8.8	29.7	123.3	9.1	28.5	117.8	10.2	18.5
	C13	171.6	7.8	36.6	107.1	7.5	28.3	95.9	6.6	25.3	160.8	9.0	36.1
MEO	C11	19.7	5.7	4.4	20.0	4.3	4.3	21.5	5.2	4.5	21.7	5.6	4.3
	C12	23.5	4.9	6.0	30.4	5.4	6.1	20.6	5.8	5.7	23.0	4.9	5.7
	C14	26.8	5.0	5.3	26.9	4.8	5.1	24.4	5.1	4.8	26.9	5.0	5.2

lists the RMS values for the four solutions. Similar to the results of SLR residuals, the contribution of the TRR parameter to orbit accuracy was mainly in the radial direction, while degeneration along the track was observed for several satellites, especially for the IGSOs. This needs further investigation. Moreover, the accuracy of predicted orbits from ABW/ABW + TRR for BDS-2 MEOs was comparable to that of the ECOM/ECOM + TRR solution, while the IGSO satellites had worse performance due to poor tracking geometry compared to that of MEOs. In order to illustrate this, Figure 12 plots the average RMS values for the BDS-2 IGSO and MEO groups.

6. Discussion

With the analysis of the SLR residuals for BDS-2 IGSO and MEO orbits from three IGS MGEX ACs, i.e., CODE, GFZ, and WHU, the noticeable systematic orbit errors dependent on the Sun elongation angle were identified for each satellite, but with different magnitude. Moreover, these errors cannot be reduced by using the standard ABW model, and different patterns were identified for the orbit solution of C08 and C13, which are in the same orbital plane and are manufactured based on the same satellite platform. This fact suggests that this $ϵ$-angle-dependent radial orbit errors do not origin from tracking geometry as well as mismodeling SRP perturbation. Hence, the TRP is considered as the most possible root.

Through analysis of the estimated parameters of ECOM, no significant differences were observed for the Y0 parameter in and out of the eclipse region. Hence, the force generated by a potential radiator along the Y axis can be excluded. Similar to [31,32], a constant acceleration along the X axis was introduced into ABW model to account for the potential thermal radiation generated by the radiator along the X axis. However, the Sun elongation angle-dependent radial orbit errors still exist. This makes us think that the TR perturbation can hardly be the reason for this type of orbit error.

Finally, a periodic acceleration along the illuminated +X axis was introduced to compensate for possible TRR-induced forces that are not fully absorbed by the parameter of the SRP model. This compensating TRR parameter was then adjusted with real ground tracking data. When using the a priori acceleration to augment ECOM, a noticeable improvement validated by SLR was observed for the BDS-2 IGSO and MEO satellites. More importantly, the Sun elongation angle-dependent orbit error pattern almost vanished. In terms of the orbit prediction performance, the proposed model seemed to deteriorate the along-track for the ECOM in many cases, especially for the IGSO satellites, while for the ABW, there was an improvement in view of orbit prediction. This needs further investigation.

Despite the success of the a priori model to compensate the TRR force for BDS-2 IGSO and MEO satellites, further improvements can be expected from more detailed information of the internal thermal emission mechanism and solar activity. In addition, the accuracy of the BDS orbit solution reaching a similar level of the GPS satellites might still need long-term reprocessing and modification, and further investigations can be expected from the improvements in orbit dynamic models for the derived geophysical parameters [25,41,42].

7. Conclusions

In this contribution, we focus on investigating the Sun elongation angle-dependent orbit errors for BDS-2 IGSO and MEO satellites, which is evidenced in the SLR residuals when using the empirical 5-parameter ECOM, and seeking a possible solution for their mitigation. Unfortunately, this mismodelling deficiency cannot be compensated by the standard ABW model. Based on analysis of the series of ECOM parameters as well as the distribution of SLR residuals on the orbital angle and the Sun elevation angle, the thermal re-radiation of the satellite bus was considered the potential cause. Hence, a periodic acceleration in the normal direction of the +X surface, presumably generated by the mismodeling TRP, was accordingly proposed and calibrated with real observations. The validation results of the four solutions listed in Table 5 indicate that when the TRR parameter was applied, the Sun elongation angle-dependent systematic errors were significantly reduced. Moreover, SLR validation indicated that the orbit accuracy improved by around 10–30% to approximately 4.5 cm and 3.0 cm for the BDS-2 IGSO and MEO satellites, respectively. In view of the metric of the one-day predicted orbits, orbit accuracy improvements can be also observed in the radial direction.