Document [original]





Citation: Wang, C.; Guo, J.; Zhao, Q.;

Ge, M. Improving the Orbits of the

BDS-2 IGSO and MEO Satellites with

Compensating Thermal Radiation

Pressure Parameters. Remote Sens.

2022,14, 641. https://doi.org/

10.3390/rs14030641

Academic Editor: Yunbin Yuan

Received: 1 December 2021

Accepted: 25 January 2022

Published: 28 January 2022

Publisher’s Note: MDPI stays neutral

with regard to jurisdictional claims in

published maps and institutional affil-

iations.

Licensee MDPI, Basel, Switzerland.

This article is an open access article

distributed under the terms and

conditions of the Creative Commons

Attribution (CC BY) license (https://

creativecommons.org/licenses/by/

4.0/).

remote sensing

Article

Improving the Orbits of the BDS-2 IGSO and MEO Satellites

with Compensating Thermal Radiation Pressure Parameters

Chen Wang 1,2,3, Jing Guo 2,*, Qile Zhao 2,4 and Maorong Ge 3,5

1School of Geology and Geomatics, Chang’an University, Yanta Road, Xi’an 710000, China;

[email protected]

2GNSS Research Center, Wuhan University, Luoyu Road, Wuhan 430079, China; [email protected]

3Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, Strasse des 17. Juni 135,

10623 Berlin, Germany; [email protected]

4Collaborative Innovation Center of Geospatial Technology, Luoyu Road, Wuhan 430079, China

5Deutsches GeoForschungsZentrum GFZ, Telegrafenberg, 14473 Potsdam, Germany

*Correspondence: [email protected]

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.

Keywords: BDS-2; precise orbit determination; thermal radiation pressure; solar radiation pressure

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 [

]. 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) mod-

els [

], observation correction models [

–

] or onboard observations from low-earth

orbits [

] and MEO satellites [

]. With the adoption of the empirically derived SRP model

for BDS GEO satellites, the radial orbit accuracy reached 10 cm as validated by Satel-

lite Laser Ranging (SLR) [

]. 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–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 [

]. 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

Remote Sens. 2022,14, 641. https://doi.org/10.3390/rs14030641 https://www.mdpi.com/journal/remotesensing

Remote Sens. 2022,14, 641 2 of 19

ECOM, [

]) for BDS-2 IGSO satellites [

]. The systematic orbit errors were also identified

in the International GNSS Service (IGS, [

])-combined Multi-GNSS orbit solutions [

The radial orbit accuracy relying on the dynamic model is essential for GNSS solutions and

scientific applications [

–

]. 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) [

], 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 [

] 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.,

eD,

eY, and

eB, expressed as

eD=

→

r−→

k→

r−→

eY=−

→

r×→

k→

r×→

eDk

eB=

eD×

(1)

where

→

r

and

→

are the geocentric vectors of the Sun and the satellite, respectively. The

vector

→

is the unit vector in the satellite–Sun direction;

→

points along the satellite’s

solar panel axis and coincides with the Yaxis of the satellite reference frame in the YS mode.

→

completes the right-hand orthogonal system. With this frame, the SRP accelerations

Remote Sens. 2022,14, 641 3 of 19

are expressed as Fourier series using the satellite’s argument of latitude

as an angular

argument,

D(u)=D0+∑n

i=1[Di,ccos(iu)+Di,ssin(iu)]

Y(u)=Y0+∑n

i=1[Yi,ccos(iu)+Yi,ssin(iu)]

B(u)=B0+∑n

i=1[Bi,ccos(iu)+Bi,ssin(iu)]

(2)

where

is an index, nis the terminated order and

Di,c

Di,s

Yi,c

Yi,s

Bi,c

, and

Bi,s

are the

amplitudes of the periodic signals. For the original ECOM model, nwas selected as 1, i.e.,

three constant and six periodic parameters had to be estimated [

]. Usually, the reduced

ECOM model [

], 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:

D(u)=D0

Y(u)=Y0

B(u)=B0+B1,ccos(u)+B1,ssin(u)

(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 notice-

able errors in the GPS-derived geodetic parameters, i.e., station coordinates, earth rotation

parameters [

]. In order to overcome this deficiency, the adjustable box-wing (ABW)

model was proposed based on the physical interaction between solar photons and the satel-

lite surface by approximating the satellite surface structure as a box (bus) with two wings

(solar panels) in YS mode [

]. 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

depends on the relative alignment of

the Sun’s direction (

→

e

) and the surface normal

→

along with the incident angle (

the fraction of absorbed photons (

), of diffusely reflected photons (

), and of specularly

reflected photons (

). For an illuminated surface, the acceleration

→

aSRP

at the distance of

1 astronomical unit (AU) can be obtained as follows:

→

aSRP =−A

ccos θ(α+δ)→

e+2δ

3+ρcos θ→

eN, (4)

where

is the mass of the satellite,

is the solar flux at 1 AU; a revised value for

based on re-analysis and re-calibration of satellite data of 1361 W/m

was used [

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

→

aTRR

can be obtained

approximately according to Lambert’s law as follows [29]:

→

aTRR =−A

ccos θ2

3α→

eN(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):

→

a=−A

ccos θ(α+δ)→

e+2

→

eN+2ρcos θ→

eN(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.

Remote Sens. 2022,14, 641 4 of 19

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 [

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. [

]

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 [

]. For solar panels, the acceleration

→

aTRR SP

due to the thermal emission of

radiation can be computed as follows:

→

aTRR_SP =−2AσεfT4

f−εbT4

b

3Mc

→

eN(7)

where

is the Stephan–Boltzmann constant (5.6699

−9

−2

−4

εf

and

εb

are the

emissivity of the front

and back

panels, and

and

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

→

aTRR_BUS can be computed as follows:

→

aTRR_BUS =−2AεσT4

MLI

3Mc

→

eN(8)

where

TMLI

is the absolute temperature in Kelvin of the MLIs, and

MLI

can be represented

as follows:

MLI =αS0cos θ+εe f f σT4

sat

σεMLI +εe f f (9)

where

εe f f

is the effective emissivity,

εMLI

is the emissivity of the MLI,

is the absorptivity

of the MLI, and Tsat 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

ε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 [

], this emission will also

produce a recoil acceleration as follows:

→

aTR =−2qdissipated

3Mc

→

eN(10)

Remote Sens. 2022,14, 641 5 of 19

where

qdissipated

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 [

]. 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 [

]. 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 [

]. 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 [

]. 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 [

], 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, [

]). 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, where the slope of the SLR residuals with

respect to the Sun elongation is also listed.

Remote Sens. 2022,14, 641 6 of 19

Remote Sens. 2022, 14, x FOR PEER REVIEW 6 of 18

Remote Sens. 2022, 14, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing

Figure 1. SLR residuals of C13 with distribution on the Sun elongation angle (ϵ) for CODE (top),

GFZ (middle), and WHU (bottom).

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/per-

formance/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 As-

sessment 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 main-

land. Almost the same POD strategy was used as for the Wuhan University MGEX prod-

ucts [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 ap-

proach. 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–11]. For the modeling aspect of antenna

thrust for BeiDou IGSO and MEO satellites, measured values from [38] were adopted.

Table 2 summarizes the background geophysical models used, and Table 3 lists the four

orbit solutions along with the estimated parameters.

Figure 1.

SLR residuals of C13 with distribution on the Sun elongation angle (

) for CODE (

top

), GFZ

(middle), and WHU (bottom).

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)−1).

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

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 [

] and the International GNSS Monitoring and As-

sessment System (iGMAS, [

]) 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 2shows 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 [

Remote Sens. 2022,14, 641 7 of 19

except for the SRP model. All data were processed by a modified version of Position And

Navigation Data Analyst (PANDA) software [

] using a two-step approach. For obser-

vation 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 con-

tinuous YS model was used [

–

]. For the modeling aspect of antenna thrust for BeiDou

IGSO and MEO satellites, measured values from [

] were adopted. Table 2summarizes

the background geophysical models used, and Table 3lists the four orbit solutions along

with the estimated parameters.

Remote Sens. 2022, 14, x FOR PEER REVIEW 7 of 18

Remote Sens. 2022, 14, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing

Figure 2. Distributions of the approximately 90 ground stations used in this study. Red circles indi-

cate IGS MGEX stations. Blue triangles represent iGMAS stations.

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, Ura-

nus, 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]

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

Solution SRP Model Estimated Parameters

ECOM The 5-parameter ECOM

model 𝐷, 𝑌, 𝐵, 𝐵 and 𝐵.

ABW The standard ABW

𝑆𝑃, 𝑆𝐵, 𝑌,+𝑋𝐴𝐷,−𝑋𝐴𝐷,+𝑋𝑅,−𝑋𝑅,+𝑍𝐴𝐷,

−𝑍𝐴𝐷,+𝑍𝑅 and −𝑍𝑅.

ABW + TRR The modified ABW with TRR

parameter 𝑆𝑃, 𝑆𝐵, 𝑌,+𝑋𝐴𝐷,−𝑋𝐴𝐷,+𝑋𝑅,−𝑋𝑅,+𝑍𝐴𝐷,

−𝑍𝐴𝐷,+𝑍𝑅,−𝑍𝑅 and 𝑘.

ECOM + TRR ECOM with additional TRR

parameter as the a prior 𝐷, 𝑌, 𝐵, 𝐵 and 𝐵.

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

3.3. The Orbits Determined with the Standard ABW Model

The standard ABW model can accurately compensate the perturbing acceleration act-

ing 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

Figure 2.

Distributions of the approximately 90 ground stations used in this study. Red circles

indicate IGS MGEX stations. Blue triangles represent iGMAS stations.

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]

Remote Sens. 2022,14, 641 8 of 19

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

Solution SRP Model Estimated Parameters

ECOM The 5-parameter ECOM model D0,Y0,B0,BcandBs.

ABW The standard ABW SP,SB,Y0,+XAD,−XAD,+XR,

−XR,+ZAD,−ZAD,+ZR and −ZR.

ABW + TRR The modified ABW with TRR parameter SP,SB,Y0,+XAD,−XAD,+XR,−XR,

+ZAD,−ZAD,+ZR,−ZR andk.

ECOM + TRR

ECOM with additional TRR parameter as the a prior

D0,Y0,B0,Bcand Bs.

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

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 [

]. 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 [

]. While this requires proper knowledge of the optical, physical, thermal,

and geometrical properties of the satellite surfaces, these values were released by [

] 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 [

]

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 [

], for the AD and R

parameters of the illuminated surfaces, respectively.

Figure 3shows 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.

Remote Sens. 2022,14, 641 9 of 19

Remote Sens. 2022, 14, x FOR PEER REVIEW 8 of 18

Remote Sens. 2022, 14, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing

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, ther-

mal, 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-depend-

ent 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 ad-

dition to mismodeled SRP perturbation, as they are in the same orbital plane, and the

geometry as well as perturbation should be similar.

Figure 3. Distribution of the SLR residual with respect to the Sun elongation angle (ϵ) for C08 (top)

and C13 (bottom) determined with the ABW model.

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-grav-

itational perturbations with similar characteristics as SRP. TRR may be the candidate as it

Figure 3.

Distribution of the SLR residual with respect to the Sun elongation angle (

) for C08 (

top

)

and C13 (bottom) determined with the ABW model.

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.

Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 18

Remote Sens. 2022, 14, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing

mainly originated from the Sun and was not fully captured in the ABW or empirical SRP

model.

Figure 4. SLR residuals distributed with respect to β and μ for the C13 ABW orbit solution.

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 parame-

ters (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 bal-

anced or are very small.

Figure 4. SLR residuals distributed with respect to βand µfor the C13 ABW orbit solution.

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.

Remote Sens. 2022,14, 641 10 of 19

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 4was caused

by TR, the orbits were re-determined with the same strategy as [

] based on the ECOM

SRP model. Figure 5illustrates the variations in the estimates of five SRP parameters (see

Equation (3)), i.e., three constant parameters along the D,Y, and Baxes (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.

Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 18

Remote Sens. 2022, 14, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing

Figure 5. Estimated parameters of the ECOM SRP model for C13 orbit solutions (The areas shaded

in grey indicate eclipse seasons).

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 intro-

duced 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.

Figure 6. Distribution of the SLR residual with respect to the Sun elongation angle (ϵ) for C13 orbits

determined with the ABW model plus a constant acceleration along the X direction.

Figure 5.

Estimated parameters of the ECOM SRP model for C13 orbit solutions (The areas shaded in

grey indicate eclipse seasons).

Furthermore, the possibility of TR along the Xaxis was also analyzed. For GLONASS

and Galileo, one additional constant radiator parameter along the -X direction was intro-

duced to account for the potential TR force generated by the radiator along the Xaxis [

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 Xaxis cannot explain the orbit deficiency. Based on the results, we conclude that TR

perturbation is not the root of the orbit errors.

Remote Sens. 2022,14, 641 11 of 19

Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 18

Remote Sens. 2022, 14, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing

Figure 5. Estimated parameters of the ECOM SRP model for C13 orbit solutions (The areas shaded

in grey indicate eclipse seasons).

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 intro-

duced 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.

Figure 6. Distribution of the SLR residual with respect to the Sun elongation angle (ϵ) for C13 orbits

determined with the ABW model plus a constant acceleration along the X direction.

Figure 6.

Distribution of the SLR residual with respect to the Sun elongation angle (

) for C13 orbits

determined with the ABW model plus a constant acceleration along the X direction.

4.2. Compensated TRR Model

Usually, thermal perturbations for solar panels can be absorbed with the parameter

along the Yaxis 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 Zaxis correspond to the radial direction, the

mismodeling perturbation along the Zaxis 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 7shows 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 4show a similar pattern as the

TRR-induced accelerations in Figure 7.

Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 18

Remote Sens. 2022, 14, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing

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 mismodel-

ing 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 im-

portantly, 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.

Figure 7. Simulated TRR-induced accelerations (unit: nm/s

) along the D direction with distribution

along the Sun elevation angle β and orbital angle μ for C13.

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

two parts:

𝑎

 =− 2

𝐴

𝜀

3𝑀𝑐𝜀 +𝜀(𝜀𝜎𝑇

+𝛼𝑆cos𝜃)𝑒 (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:

𝑎

 =− 2

𝐴

𝜀

3𝑀𝑐𝜀 +𝜀𝛼𝑆cos𝜃𝑒=−𝑘cos𝜃𝑒 (12)

with 𝑘=











. By adding Equations (6) and (12), the SRP and TRR perturbation

acting on the +X surface can be obtained:

𝑎=−

𝐴

𝑀𝑆

𝑐𝑐𝑜𝑠𝜃(α+δ)𝑒⊙+2

3𝑒+2𝜌𝑐𝑜𝑠𝜃𝑒−𝑘cos𝜃𝑒 (13)

Figure 7.

Simulated TRR-induced accelerations (unit: nm/s

) along the D direction with distribution

along the Sun elevation angle βand orbital angle µfor C13.

Remote Sens. 2022,14, 641 12 of 19

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

two parts:

→

aTRR_BUS =−2Aε

3McεMLI +εe f f εσT4

sat +αS0cos θ→

eN(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 [

]. Hence, it can be omitted without significant

loss of accuracy, leaving the remaining part as follows:

→

aTRR_BUS =−2Aε

3McεMLI +εe f f αS0cos θ→

eN=−kcos θ→

eN(12)

with

2AεαS0

3Mc(εMLI +εe f f )

. By adding Equations (6) and (12), the SRP and TRR perturbation

acting on the +X surface can be obtained:

→

a=−A

ccos θ(α+δ)→

e+2

→

eN+2ρcos θ→

eN−kcos θ→

eN(13)

Therefore, the compensated TRR scale factor

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

, the correlations with the orbital as well as ABW

parameters were investigated. Figure 8shows the correlations among those parameters for

C13 on Day of Year (DOY) 200, 2018. Generally, the correlations of the thermal scale factor

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.

Remote Sens. 2022, 14, x FOR PEER REVIEW 12 of 18

Remote Sens. 2022, 14, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing

Therefore, the compensated TRR scale factor 𝑘 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 𝑘, 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 𝑘 with the ABW parameters were less than 0.5 except for those with the SP param-

eter, 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 𝑘.

Figure 8. Correlations between the satellite state, ABW as well as the thermal scale factor 𝑘 param-

eters of C13 on DOY 200, 2018 (β angle of approximately 20 deg).

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 sat-

ellites showed a similar behavior. The averaged estimates derived from the daily estimates

outside the eclipse regime for each satellite are listed in Table 4. 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.

Figure 9. Daily estimates of the thermal scale factors for C13 with distribution on the β angle. (The

grey represents the eclipse season).

Figure 8.

Correlations between the satellite state, ABW as well as the thermal scale factor

parameters

of C13 on DOY 200, 2018 (βangle of approximately 20 deg).

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 9demonstrates 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. 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

Remote Sens. 2022,14, 641 13 of 19

the largest estimate of approximately 2.6 mm/s

, and the smallest was about 0.9 nm/s

for

C14.