Evaluation of strategies for the ultra-rapid orbit prediction of BDS GEO satellites

ABSTRACT

The quality of BeiDou Navigation Satellite System (BDS) Geostationary Earth Orbit (GEO) ultra-rapid products is unsatisfactory because GEO satellites are nearly stationary relative to ground stations. To optimize the quality of these ultra-rapid orbit products, we investigated the effects of the fitting arc length, an a priori Solar-Radiation Pressure (SRP) model, and the along-track empirical acceleration on the prediction of BDS GEO satellite orbits. The predicted orbit arcs of 24-h were evaluated through comparisons with the corresponding observed orbit arc and Satellite Laser Ranging (SLR) observations. In both eclipse and non-eclipse seasons, accuracy of the orbit predictions obtained using a 48-h fitting arc length were better than those obtained using 24-h and 72-h fitting arc lengths. Although the overlapping precision of predicted orbits exhibited no obvious improvement when an a priori SRP model was employed, the systematic bias in the SLR residuals was significantly reduced. Specifically, the mean value of SLR residuals decreased from −0.248 m to −0.024 m during non-eclipse seasons and from −0.333 m to −0.041 m during eclipse seasons, respectively. In addition, when an empirical acceleration in the along-track direction was introduced, the three-Dimensional Root-Mean-Square (3D RMS) of overlapping orbits during eclipse seasons decreased from 2.964 to 1.080 m, which is comparable to that during non-eclipse seasons. Furthermore, the Standard Deviation (STD) of SLR residuals decreased from 0.419 to 0.221 m during eclipse seasons. The analysis of SRP estimates shows that the stability of SRP parameters was significantly enhanced after the introduction of along-track empirical acceleration in eclipse seasons. The optimal BDS GEO ultra-rapid orbit prediction products were yielded by using a 48-h fitting arc length, an a priori SRP model and an along-track empirical acceleration.

KEYWORDS:

• BDS GEO satellites
• orbit prediction
• solar radiation pressure model

1. Introduction

Ultra-rapid products are essential for real-time and near-real-time applications of Global Navigation Satellite System (GNSS), such as rapid orbit determination for satellites in low-Earth orbit, clock monitoring, atmospheric monitoring, and precise point positioning (Montenbruck, Gill, and Kroes 2005; Heo, Cho, and Heo 2011; Springer and Hugentobler 2001; Li et al. 2019). The International GNSS Service (IGS) has been officially providing ultra-rapid products since Global Positioning System (GPS) week 1087 (starting 5 November 2000; https://kb.igs.org/hc/en-us/articles/203866356). To integrate new constellations, such as the BeiDou Navigation Satellite System (BDS) (Lu, Guo, and Su 2020), Galileo, and the Quasi-Zenith Satellite System (QZSS), into the IGS (Beutler et al. 1999) processing routine, the multi-GNSS experiment (MGEX) was initiated in 2013. Ultra-rapid products from GPS, Globalnaya Navigazionnaya Sputnikovaya Sistema (GLONASS), Galileo, BDS, and QZSS satellites are currently provided by several IGS Analysis Centers (ACs), such as the Center for Orbit Determination in Europe (CODE), GeoForschungsZentrum (GFZ), the European Space Agency (ESA), and Wuhan University (WHU) (Montenbruck et al. 2017a). WHU also provides BDS hourly ultra-rapid products including Geostationary Earth Orbit (GEO) satellites (http://mgex.igs.org/IGS_MGEX_Products.php). Each ultra-rapid product consists of observation data of the previous 24 h of orbits and predictions for the next 24 h.

The BDS GEO satellites, which are distributed over the Indian and Pacific Oceans, play an important role in geometry strength enhancement (Hao et al. 2018). Moreover, they support Two-Way Satellite Time and Frequency Transfer (TWSTFT) (Gao et al. 2011). These are exploited in multi-GEO TWSTFT links between two stations, which can improve time and frequency transfer (Yang et al. 2014). However, the application of BDS GEO satellites is limited by the poor accuracy of their orbit products. BDS GEO satellites are characterized by a near-static observation geometry; thus, to ensure synchronization with Earth rotation they maneuver frequently (Tu et al. 2021), which impedes the determination and prediction of BDS GEO satellite orbits (Cao et al. 2014). To improve Precise Orbit Determination (POD) for BDS GEO satellites, several studies have been conducted to enhance observation geometry and perturbation modeling. Onboard observations of Low Earth Orbit (LEO) or Medium Earth Orbit (MEO) satellites have proven to be an effective approach for improving observation geometry (Zhao et al. 2017; Ge et al. 2022). Ge et al. (2017) showed that accuracy of BDS GEO satellites could be considerably improved from meters to decimeters mainly in the along-track direction. In addition, the Orbit-Normal (ON) attitude mode used by BDS GEO satellites hinders the accurate modeling of Solar Radiation Pressure (SRP). IGS ACs widely use the extended Empirical CODE Orbit Model (ECOM, ECOM2) for SRP modeling. However, the ECOM was initially developed for GPS satellites with a Yaw-Steering (YS) attitude mode, and thus an optimal POD for BDS GEO satellites with the ON mode cannot be achieved using ECOM. Lou et al. (2014) identified that the ECOM Y-axis component in the ON attitude cannot fully span the SRP force component perpendicular to the solar panel. Thus, for POD of BDS GEO satellites, these researchers adopted an ECOM based on redefined forms of the three orthogonal axes, which considerably improved the orbit determination accuracy. The proposed redefined three ECOM orthogonal axes for ON satellites was then adopted in Prange et al. (2020) to construct an empirical SRP model for improved orbit determination of QZSS satellites. Duan, Hugentobler, and Selmke (2018) applied an a priori box-wing model to enhance the ECOM-ON model for the POD of BDS GEO satellites. This afforded better orbit quality than that obtained with the pure ECOM. However, BDS GEO satellites still retained a large Satellite Laser Ranging (SLR) residual offset of up to 20 cm. BDS GEO satellites carry a large Communication Antenna (CA) on their +X panel, and the CA extends over their +Z and −Z surfaces and generates perturbation along the Z-axis. As the CA is positioned on the +X side and has a large area, self-shadowing effects on the +X surface have to be accounted for. Wang et al. (2019) proposed an a priori adjustable box-wing SRP model for BDS GEO satellites, which considered the perturbation generated by the CA, as well as its shadowing effect. This afforded a Root-Mean-Square (RMS) of the BDS GEO satellite SLR residuals that was four to five times smaller than that of the ECOM, with a near-zero SLR residual bias. Zhao et al. (2013) introduced an along-track empirical acceleration parameter in the POD processing of BDS GEO satellites to compensate for the inaccuracies of SRP modeling of BDS GEO satellites. This led to a RMS of 68.5 cm for the SLR residuals from 25 normal points. However, these studies have mainly focused on post POD, and have not investigated orbit prediction performance for BDS GEO satellites.

In the process of ultra-rapid orbit prediction, the position and velocity of a satellite at any time can be obtained through integration if the dynamic model and initial values are given. Thus, the observed orbit arc and the dynamic model are crucial factors that determine the quality of ultra-rapid orbit products. Choi et al. (2013) studied the effects of the arc length of fitted observed orbits and SRP parameterization on GPS satellites orbit prediction. Their results showed that the most stable and accurate predictions were obtained using observed arc lengths of 40–45 h and a nine-parameter ECOM model. Geng et al. (2018) evaluated the impact of the fitting arc length of observed orbits and SRP model on the orbit prediction performance for multi-GNSS, including GPS, GLONASS, Galileo, and BDS satellites; however, the BDS GEO satellites were not considered. Currently, the orbit accuracy of BDS GEO ultra-rapid products provided by MGEX ACs is relatively low and insufficiently stable (Dai et al. 2019). Therefore, to improve the performance of ultra-rapid orbit products, there is a need for an efficient strategy for BDS GEO orbit prediction.

In this study, we investigated the effect of three factors – fitting orbit length, the a priori SRP model and the along-track empirical acceleration – on BDS GEO ultra-rapid orbit prediction, and used our findings to develop an improved orbit prediction method. We first present the SRP models and empirical constant-acceleration models, followed by our analysis of the determination of an a priori constraint for the SRP and empirical acceleration parameters. Three POD strategies were designed to analyze which factors affect orbit solution, and processing strategies used in BDS GEO ultra-rapid orbit prediction are presented. Then, we assessed the ultra-rapid orbit-prediction performance of various strategies in terms of orbit overlapping analysis and SLR validation. Finally, we developed a strategy that can improve the accuracy of ultra-rapid orbit prediction, which we present in the last section.

2. Methodology

In this section, the SRP models and empirical constant-acceleration model used in this study are presented first. Then, values of the constraints for the empirical acceleration parameter are discussed. Finally, the experiment and data processing strategy used for ultra-rapid POD of BDS GEO satellites are given.

2.1 The SRP model

Beutler et al. (1994) and Springer, Beutler, and Rothacher (1999) introduced the DYB frame (Figure 1) to describe SRP accelerations for GPS satellites operating in the YS mode. The DYB frame comprises three unit vectors:

Figure 1. Illustration of DYB (Sun-fixed) and XYZ (body-fixed) orthogonal frames (Rodriguez-Solano, Hugentobler, and Steigenberger 2012).

$e_D=e_{\odot }$

(1) $e_Y=\frac{e_{\odot }\times r}{\left|e_{\odot }\times r\right|}$(1)

$e_B=e_D\times e_Y$

where $e_{\odot }$ is the unit vector that points from the satellite to the sun, $e_Y$ is the unit vector along the solar-panel axis of a spacecraft, and $r$ is the geocentric vector from Earth to the satellite. $e_B$, $e_D$, and $e_Y$ constitute a right-hand coordinate system.

The reduced five-parameter ECOM (ECOM5) (Springer, Beutler, and Rothacher 1999) used in this study decomposes the SRP accelerations in the above three orthogonal directions for the YS attitude mode (subsequently designated as ECOM-YS), as follows:

$a_D=D_0$

(2) $a_Y=Y_0$(2)

$a_B=B_0+B_C\cdot cosu+B_S\cdot sinu$

where $D_0$, $Y_0$, and $B_0$ are three constant-acceleration parameters; $B_C$ and $B_S$ denote one-cycle-per-revolution (1-cpr) acceleration, and $u$ is the satellite–Earth ascending node angle.

The ECOM was initially developed for GPS satellites operating in the YS attitude mode. However, as BDS GEO satellites operate in the ON attitude mode, their performance cannot be optimized by using only the ECOM in POD. In the ON attitude mode, satellite solar panels are not perpendicular to the Sun direction. Thus, analogous to the DYB frame used for YS mode, the $\overline{\mathrm{D}\mathrm{Y}\mathrm{B}}$ frame, comprising three unit vectors ($e_{\bar{D}}$, $e_{\bar{Y}}$, and $e_{\bar{B}}$), is introduced to describe SRP accelerations for satellites in the ON mode (Montenbruck, Steigenberger, and Darugna 2017b):

$e_{\bar{D}}=e_{\bar{Y}}\times e_{\bar{B}}$

(3) $e_{\bar{Y}}=-\frac{r\times v}{\left|r\times v\right|}$(3)

$e_{\bar{B}}=\frac{e_{\odot }\times e_{\bar{Y}}}{\left|e_{\odot }\times e_{\bar{Y}}\right|}$

where, $e_{\odot }$ and $r$ have the same definitions as in Equation (1)(1) $e_Y=\frac{e_{\odot }\times r}{\left|e_{\odot }\times r\right|}$(1) ; $e_{\bar{D}}$, $e_{\bar{Y}}$, and $e_{\bar{B}}$ constitute the right-hand coordinate system; and $v$ is the satellite velocity. Likewise, the SRP acceleration parameters in each component can be expressed as follows:

$a_{\bar{D}}={\bar{D}}_0$

(4) $a_{\bar{Y}}={\bar{Y}}_0$(4)

$a_{\bar{B}}={\bar{B}}_0+{\bar{B}}_C\cdot cosu+{\bar{B}}_S\cdot sinu$

For the ECOM described in the ON attitude mode (subsequently designated as ECOM-ON), ${\bar{D}}_0$, ${\bar{Y}}_0$, and ${\bar{B}}_0$ are constant-acceleration parameters, and ${\bar{B}}_C$ and ${\bar{B}}_S$ are 1-cpr acceleration parameters. However, the orbital inclination of BDS GEO satellites is nearly 0°; consequently, the orbital plane coincides with the equatorial plane. In this case, the $e_{\bar{Y}}$ vector is aligned in the Z-direction of the celestial reference frame, and a strong linear correlation exists between the constant-acceleration parameter ${\bar{Y}}_0$ and the satellite initial position on the Z-axis of the celestial reference frame (the satellite position parameter in the Z-direction). Wang et al. (2019) found that the ECOM-ON model with a strong constraint on Y₀ parameter did not outperform the ECOM-YS model in terms of daily Orbit Boundary Discontinuities (OBDs) values and SLR residuals. Therefore, the ECOM-YS model was adopted in the current study. Wang et al. (2019) found that the introduction of the a priori model considerably improved BDS GEO satellite POD. This a priori model for the DYB frame, which they used to enhance the ECOM, is formulated as follows:

$a_{D_P}=D_{0P}+D_{1P}cos\epsilon +D_{2P}cos(2\epsilon )+D_{4P}cos(4\epsilon )$

(5) $a_{Y_P}=Y_{1P}sin\mu$(5)

$a_{B_P}=B_{0P}+B_{1P}cos\epsilon +B_{3P}cos(3\epsilon )$

where $D_{0P}$ and $B_{0P}$ are constant-acceleration parameters; $D_{1P}$, $Y_{1P}$, and $B_{1P}$ are 1-cpr acceleration parameters; $D_{2P}$, $B_{3P}$, and $D_{4P}$ are 2-cpr, 3-cpr and 4-cpr acceleration parameters, respectively; and $\epsilon$ is the Earth–satellite–Sun angle. $D_{0P}$, $Y_{1P}$, and $B_{1P}$ have $\beta$-angle-dependent variations, and are therefore fitted using a polynomial function of the $\beta$ angle. The coefficients derived for the model are listed in Table 1.

Table 1. Values of parameters in the a priori SRP model for BDS GEO satellites (Wang et al. 2019; units: nm/s²; $\beta$ angle in degrees).

Parameter	Value nm/s^2
$D_{0P}$	0.0136 ${\beta }^2-$ 112.1
$D_{1P}$	0.32
$D_{2P}$	−10.4
$D_{4P}$	−1.6
$Y_{1P}$	−0.416 $\beta$
$B_{0P}$	1.4
$B_{1P}$	−0.022 ${\beta }^2-$ 6.53
$B_{3P}$	−1.51

2.2 Constant-acceleration model

BDS GEO satellites suffer from poor orbit quality, particularly in the along-track direction (Wang et al. 2018). Thus, to compensate for the deficiencies of SRP modeling, empirical constant-acceleration is commonly introduced into the satellite orbit equation (Zhao et al. 2013; Guo 2014), which is expressed as

(6) ${\overrightarrow{a}}_{emp}=\left[\begin{array}{c}{a}_A\\ a_C\\ a_R\end{array}\right]\cdot {\overrightarrow{e}}_{ACR}$(6)

where ${\overrightarrow{a}}_{emp}$ is the empirical constant-acceleration perturbation in the inertial frame; $a_A$, $a_C$, and $a_R$ are the empirical constant-acceleration parameters in the along-track, cross-track, and radial directions, respectively; and ${\overrightarrow{e}}_{ACR}$ is the unit vector of the along-track, cross-track, and radial directions in the inertial frame. In principle, the existence of too many empirical acceleration parameters will reduce estimability of the POD system. Therefore, in the current study, only an empirical constant-acceleration parameter in the along-track direction was considered. The priori constraint of the empirical parameter strongly affects the influence of the dynamic model of Equation (6)(6) ${\overrightarrow{a}}_{emp}=\left[\begin{array}{c}{a}_A\\ a_C\\ a_R\end{array}\right]\cdot {\overrightarrow{e}}_{ACR}$(6) on the estimated orbits. Accordingly, we adjusted the constraint to better compensate for the influence of strong correlations between the along-track empirical constant-acceleration and SRP parameters, which will be discussed in more detail in Section 2.3.

2.3 POD and orbit prediction strategies

The precision of observed orbits will affect that of the predicted orbits. To maintain the consistency of the dynamic model used for orbit calculation and prediction, we collected tracking data from 97 MGEX stations (Montenbruck et al. 2014, 2017a; Villiger and Dach 2021) and 13 International GNSS Monitoring and Assessment System (iGMAS) (Jiao 2014) stations from Day of Year (DOY) 001, 2018, to DOY 365, 2019, to determine the observed orbits. Figure 2 shows the distribution of 110 stations selected from the MGEX and iGMAS tracking networks. For POD data processing, we used Positioning and Navigation Data Analyst (PANDA) software (Liu and Mao 2003; Shi et al. 2008) developed by the GNSS Research Center of Wuhan University. During POD processing, the code and phase observations on GPS L1/L2 and BDS B1I/B3I were used to form ionosphere-free combinations that were jointly processed (Lou et al. 2016). The batch estimation least-squares processing mode was adopted for the POD. Table 2 summarizes the main orbit dynamic models used in the POD strategies.

Figure 2. Distributions of MGEX (blue) and iGMAS (red) stations used in ultra-rapid orbit prediction experiment.

Table 2. Orbit dynamic models.

Models	Description
Geopotential (static)	EGM 08 up to 12 × 12 (Pavlis et al. 2012)
N-body	Sun, Moon, Jupiter, Venus, Mars Mercury, Uranus, Neptune, Saturn, Pluto JPL DE405 ephemeris (Standish et al. 1998)
Solid earth tides	International Earth rotation and Reference systems Service (IERS) conventions 2010 (Petit and Luzum 2010)
Ocean tides	IERS conventions 2010
Solid earth pole tides	IERS conventions 2010
Relativistic effects	IERS conventions 2010
Solar radiation pressure	Five-parameter ECOM with and without an a priori model (Wang et al. 2019; Springer, Beutler, and Rothacher 1999)
Earth radiation pressure	None
Orbit integration	Runge–Kutta–Fehlberg 7(8), Adams–Bashforth, and Adams–Moulton methods (Montenbruck and Gill 2012)

First, the a priori constraint of the empirical along-track acceleration parameter in orbit determination was defined. The orbit solution quality will deteriorate if inappropriate constraints are applied. An a prior constraint of 0.1 nm/s² (Guo 2014) is widely used for empirical along-track acceleration parameter; however, the correlations between the empirical along-track acceleration parameter and SRP parameters have never been analyzed. In this study, these correlations were investigated using three solutions with different arc lengths, including 24-h, 48-h, and 72-h. Figure 3 shows the correlation coefficients between the empirical along-track acceleration and SRP parameters as a function of the elevation angle of the Sun above the satellite orbital plane ($\beta$ angle) for the BDS C01 satellite from DOY 179 to 356, 2018, these days contain a complete season of eclipse and non-eclipse. Interestingly, the correlations between the empirical along-track acceleration and D₀ or B₀ parameters were nearly equal to 1 in both eclipse and non-eclipse seasons for the 24-h arc length solution. With 48-h and 72-h arc lengths, the correlations between the empirical along-track acceleration and D₀,B₀,B_C, and B_S parameters were less than 0.5. In addition, with 24-h, 48-h, and 72-h arc lengths, the correlations between the empirical along-track acceleration and Y₀ parameters exceed 0.5 in non-eclipse seasons. However, the empirical along-track acceleration and Y₀ parameters were distinguishable in eclipse seasons when the $\beta$ angle was nearly equal to zero degree. Thus, we adjusted the a prior constraint for the empirical along-track acceleration parameter to 0.05 nm/s² for the 24-h arc length solution, so that the correlation coefficients between the empirical along-track acceleration and SRP parameters were all less than 0.5. For the 48-h and 72-h arc length solutions, a loose a priori constraint of 10 nm/s² was adopted for the empirical along-track acceleration parameter, to better compensate for the orbit modeling error. Moreover, in non-eclipse seasons, the a priori constraint of 0.1 nm/s² was added to the Y₀ parameter, to distinguish the Y₀ and empirical acceleration parameters.

Figure 3. Correlations between empirical along-track acceleration and SRP parameters with 24-h, 48-h, and 72-h orbit arc lengths. Other BDS GEO satellites show similar performances.

During the orbit prediction processing, 24-h observed orbits were concatenated to build observation arcs from 24-h to 72-h to predict the forward 24-h orbits. Furthermore, 24-h observed orbits were generated by different strategies, using BDS and GPS daily observations. Although the accuracies of the observed orbits obtained from these various strategies differed, the observed orbits were still adopted to generate the predicted orbits, as this maintained the consistency of the model used in orbit determination and prediction. To determine which factors influenced the orbit solution, we designed a group of POD strategies for comparison. Strategy 1 uses only the five-parameter ECOM SRP model. In previous studies, the introduction of an a priori model markedly improved the performance of the BDS GEO satellite POD (Duan, Hugentobler, and Selmke 2018; Wang et al. 2019). Thus, strategy 2 employs ECOM5 with an a priori model. To further improve the accuracy of BDS GEO satellite orbit determination, strategy 3 introduces an along-track empirical constant-acceleration parameter.

After 24-h-observed orbit products were generated using these strategies, the predicted orbits were computed using the corresponding models. First, the 24-h-observed orbits were concatenated into 24-h, 48-h, and 72-h arcs. The observed orbits were rotated from the Earth-centered Earth-fixed (ECEF) frame to the Earth-Centered Inertial (ECI) frame using known Earth rotation parameters from IERS Bulletin A, to form pseudo-observations. The satellite precise initial orbit parameters were adjusted according to the orbit dynamic model. Then, 24-h, 48-h, and 72-h fitted orbits were integrated to generate a 24-h-predicted orbit using the precise initial parameters. Finally, the orbits were rotated back to the ECEF frame, and ultra-rapid orbit products including the 24-h-predicted arc were produced. The 24-h-predictions of orbit arcs were evaluated by comparisons with the corresponding observed orbit arc and SLR observations. The workflow is illustrated in Figure 4.

Figure 4. Workflow of orbit prediction and comparison.

3. Results and discussion

In this section, we first report our assessments of the precision for the observed orbit. Then, we describe our analyses for the accuracy of the predicted orbits, which we determined by overlapping the predicted orbits with post-processed observation arcs and comparing the predicted orbits with SLR observations. Finally, the SRP parameters estimated through different POD strategies are analyzed.

3.1 The observed orbit quality

The precision of the observed orbit was evaluated via OBDs analysis and SLR validation. The daily OBDs time series of observed orbits for strategies 1 (ECOM5), 2 (ECOM5 + a priori SRP model), and 3 (ECOM5 + a priori SRP model + empirical along-track acceleration) for BDS C01 satellite are shown in Figure 5. The OBDs accuracies of the three strategies for the along-track direction were comparable. However, in the radial direction, considerable periodic systematic errors were evident. The three strategies, especially strategy 1, showed considerably reduced orbit accuracies in eclipse seasons. Moreover, the systematic error in strategy 3 was more pronounced after the application of empirical along-track acceleration. The daily OBDs RMS and three-Dimensional (3D) RMS statistical results of all BDS GEO satellites are listed in Table 3. Compared with strategy 1, strategy 2, which includes an a priori SRP model, exhibited better OBDs precision in the radial and cross-track directions, and a slightly variation in the along-track direction. Strategy 3, which includes empirical along-track acceleration, obtained more accurate OBDs than strategy 2 in both the cross-track and along-track directions for all BDS GEO satellites, especially in the cross-track direction, but slightly less accurate OBDs in the radial direction for BDS C01, C04, and C05 satellites.

Figure 5. Daily OBDs time series of observed orbits for strategies 1, 2, and 3 in the radial (R), cross-track (C), and along-track (A) directions for BDS C01 satellite (unit: m). The gray blocks indicate the satellite eclipse seasons. The black lines indicate the elevation angle of the Sun above the satellite orbital plane. Other BDS GEO satellites showed similar performance.

Table 3. The daily OBDs RMS values of BDS GEO satellites.

Satellite	Strategy	R	C	A	3D
C01	1	0.345	0.133	0.490	0.614
	2	0.218	0.124	0.499	0.559
	3	0.258	0.053	0.490	0.563
C02	1	0.304	0.151	0.487	0.594
	2	0.228	0.143	0.486	0.556
	3	0.174	0.069	0.469	0.505
C03	1	0.415	0.147	0.664	0.797
	2	0.402	0.138	0.622	0.753
	3	0.324	0.056	0.620	0.702
C04	1	0.419	0.149	0.560	0.715
	2	0.385	0.131	0.570	0.701
	3	0.397	0.055	0.560	0.689
C05	1	0.410	0.138	0.582	0.725
	2	0.370	0.130	0.581	0.703
	3	0.379	0.068	0.563	0.682

The SLR technique is an independent validation tool for orbit products (Combrinck 2010). The SLR validation approach for GNSS orbits is well known and applied by several IGS ACs (Prange et al. 2017; Uhlemann et al. 2015). The International Laser Ranging Service (ILRS) was established in September 1998 to support geodetic and geophysical research programmes (Pearlman, Degnan, and Bosworth 2002). The current ILRS network includes more than 40 SLR stations (Wilkinson et al. 2019), which typically use three different types of detectors: Compensated Single-Photon Avalanche Diode (CSPAD), Micro-Channel Plate (MCP) detectors, and Photo-Multiplier Tube (PMT) detectors. The observation qualities of CSPAD and MCP detectors are similar, whereas PMT detectors have the worst observation quality, as verified by Strugarek et al. (2021). Although all BDS GEO satellites are equipped with Laser Retroreflector Arrays (LRAs), the coordinate parameters of the laser reflection prism relative to the satellite center of mass are only available for the BDS C01 satellite from the ILRS. Thus, the BDS C01 satellite was selected for SLR validation.

Five stations equipped with a CSPAD detector (KOML, STL3, SHA2, CHAL, and BEIL) and one MCP-equipped station (YARL), can track the BDS C01 satellite. Because the stations used for SLR validation in this study do not use PMT detectors, separate statistics and combined statistics are slightly different. The SLR residuals, which are the differences between SLR observations and the range calculated from microwave-based stations to satellite positions from the predicted orbit arc, typically showed GNSS orbit accuracy in the radial direction. From 2018 to 2019, there were 1,430 SLR normal points in non-eclipse seasons and 708 in eclipse seasons, which we used in our study. SLR residuals with absolute values exceeding 2.5 m were excluded. The time series of SLR residuals for strategies 1 to 3 for BDS C01 satellite are presented in Figure 6. Table 4 also presents the mean and Standard Deviation (STD) values of SLR residuals for strategies 1 to 3. Strategy 2, which includes an a priori SRP model, obtained less scattered SLR residuals than strategy 1. The mean and STD of the residuals decreased from −0.303 m to −0.026 m, and 0.210 m to 0.103 m, respectively. After the empirical along-track acceleration was added, the STD of SLR residuals was almost unchanged, while the mean SLR residual decreased to −0.019 m. The external accuracy was improved, although the inner accuracy was slightly decreased in the radial direction (Table 3).

Figure 6. SLR residual time series of observed orbits for strategies 1–3 for BDS C01 satellite (unit: m).

Table 4. Mean and STD of SLR residuals for observed orbits for BDS C01 satellite (unit: m), as determined by strategies 1–3.

Strategy	Mean	STD
1	−0.303	0.210
2	−0.026	0.103
3	−0.019	0.104

3.2 Predicted orbit quality

The orbit overlap difference reflects the inner consistency of orbit strategies, which is a useful measure of orbit accuracy. We generated orbit differences with an overlapping arc length of 24-h, and a comparison of the last 24-h predicted orbit with the observed orbit revealed orbit differences in the radial, cross-track, and along-track directions. Figure 7 shows the daily RMS values of overlap differences in the radial, cross-track, and along-track directions for the different fitting arc lengths of BDS C01 satellite, and the $\beta$ angle. For strategies 1–3, the RMS values of orbit overlap differences in the radial, cross-track, and along-track directions varied with the fitting arc length. In non-eclipse seasons, the RMS values of overlap differences obtained through strategies 1–3 were similar. In the radial and along-track directions, the 48-h fitting arc length contributed the best overlap precision, followed by the 72-h and 24-h fitting arc lengths. In the cross-track direction, strategy 1 achieved the best overlap precision for a 24-h fitting arc length, while in strategies 2 and 3, the 24-h fitting arc length resulted in the worst precision. In eclipse seasons, strategies 1 and 2 showed similar performance. Moreover, in contrast to the non-eclipse seasons, the overlap precision in all directions considerably deteriorated for 48-h and 72-h fitting arc lengths, but slightly changed for the 24-h fitting arc length. The largest orbit deviation occurred at the epochs when the elevation angle of the Sun is nearly zero degrees. Strategy 3 with the 48-h or 72-h fitting arc length obtained higher overlap precision in the three directions than the other two strategies, and the overlap precision achieved the same level as that in non-eclipse seasons. It is notice that the β-dependent periodic errors in the cross-track direction were obvious in all strategies which may be caused by the imperfect SRP model. Firstly, although BDS GEO satellites are in ON mode, the ECOM-YS model was adopted in this study for BDS GEO POD. According to the study by Hauschild et al. (2012), the received SRP in ON-mode is approximately proportional to cos(β). It means that the difference of received SRP between YS-mode and ON-mode reduces when the β angle gets smaller and it increases when the β angle gets larger. In addition, the correlations between SRP parameters and the PZ (satellite position parameter in the Z direction of the celestial reference frame) parameter were relatively larger than the other directions as shown in Figure 8. And the errors in the PZ direction were more projected to the cross-track direction than along-track and radial directions. This could be a possible cause of the systematic errors in the cross-track direction. Therefore, the systematic errors cannot be solved by only introduction of empirical accelerations in the three orbit directions. More investigation is needed for further optimizing of the SRP model.

Figure 7. Daily RMS of BDS C01 orbit overlap differences in the radial, cross-track, and along-track directions for the fitting arc lengths of 24-h, 48-h, and 72-h with different strategies. Other BDS GEO satellites show similar performances.

Figure 8. Correlations between SRP parameters and the PX, PY, PZ (satellite position parameter in the X, Y and Z direction of the celestial reference frame) parameters, VX, VY VZ (satellite velocity parameter in the X, Y and Z direction of the celestial reference frame) parameters of 48-h orbit arc length for BDS C01 satellite when β is equal to 21 degrees, parameter correlations of 24-h and 72-h orbit arc length show similar characteristics.

Table 5 further presents the corresponding average RMS values of all BDS GEO satellites in the three directions, and the 3D RMS for strategies 1–3 with fitting arc lengths of 24-h, 48-h, and 72-h, in eclipse and non-eclipse seasons. In the solutions of three strategies, the 48-h fitting arc length contributes lower 3D RMS values than the 24-h and 72-h fitting arc lengths in both eclipse and non-eclipse seasons. The 3D RMS of strategy 2 obtained using the 48-h fitting arc length was similar to that of strategy 1, indicating that the a priori SRP model did not improve the orbit overlap precision. Strategy 3, which includes along-track empirical acceleration, led to considerably improved overlap precision in eclipse seasons. The overlap precision in eclipse seasons and non-eclipse seasons in the radial, cross-track, and along-track directions, and the 3D RMS, were improved. In eclipse seasons, strategy 3 with the 48-h fitting arc length obtained average RMS values of 0.268 m and 1.080 m in the radial and 3D components, respectively, 52.5% and 63.6% less than the corresponding values for strategy 2. Moreover, strategy 3 achieved the same level of overlap precision for the eclipse and non-eclipse seasons, with 3D RMS values of 1.080 m and 0.916 m, respectively.

Table 5. Average RMS values of all BDS GEO satellites in the three directions and the 3D component for strategies 1–3 with fitting arc lengths of 24-h, 48-h, and 72-h for eclipse (E) and non-eclipse (NE) seasons (unit: m).

Strategy	Fitting arc length (h)	R		C		A		3D
		E	NE	E	NE	E	NE	E	NE
1	24	1.220	0.951	0.117	0.197	5.630	3.145	5.775	3.287
	48	0.510	0.192	0.857	0.248	2.689	0.765	2.946	0.851
	72	0.499	0.272	1.515	0.361	2.927	1.031	3.514	1.168
2	24	1.432	1.325	0.210	0.464	5.624	5.478	5.842	5.683
	48	0.564	0.242	0.763	0.282	2.703	0.845	2.964	0.959
	72	0.550	0.264	1.504	0.365	2.978	0.985	3.558	1.130
3	24	1.432	1.330	0.209	0.462	5.627	5.528	5.845	5.727
	48	0.268	0.227	0.185	0.215	1.007	0.836	1.080	0.916
	72	0.287	0.265	0.117	0.169	1.042	0.984	1.105	1.050

The above analysis only reflected internal accuracy of all BDS GEO satellites. The same SLR validation method as used in Section 3.1 was adopted to verify external accuracy of orbit prediction for a comprehensive inspection. SLR residuals computed from the 24-h predicted orbits of strategies 1–3 for BDS G01 satellite are shown in Figure 9. The lack of SLR observations resulted in discontinuities in the residual series. With a 24-h fitting arc length, SLR residuals exhibited the most pronounced dispersion and systematic error. During non-eclipse seasons and with 48-h and 72-h fitting arc lengths, the SLR residuals of strategy 2 were less dispersed and contained less systematic error than those of strategy 1, while the results for eclipse seasons were still worse. During eclipse seasons and with 48-h and 72-h fitting arc lengths, the SLR residuals of strategy 3, which included the empirical constant acceleration in along-track direction, were considerably less dispersed than those of strategy 2.

Figure 9. Time series of SLR residuals for predicted orbits with 24-h, 48-h, and 72-h fitting arc lengths of BDS C01 satellites (unit: m). Other BDS GEO satellites show similar performances.

The mean and STD of SLR residuals for predicted orbits for BDS G01 satellite are presented in Table 6. For strategies 1–3, the STDs with a 48-h fitting arc length were less than those with a 24-h or 72-h fitting arc length. With a 48-h fitting arc length, strategy 2 obtained smaller biases and STDs than strategy 1. Moreover, the STDs of the SLR residuals of strategy 2 in the eclipse and non-eclipse seasons (0.419 m and 0.213 m, respectively) were less than those of strategy 1 (0.491 m and 0.291 m, respectively). This means that introducing the a priori SRP model significantly reduced the large negative SLR biases. In strategy 3, the introduction of empirical along-track acceleration further reduced the STD in eclipse seasons; the STD of SLR residuals (0.221 m) was 47.0% less than that of strategy 2 (0.419 m). However, the difference was not pronounced under non-eclipse conditions.

Table 6. Means and STDs of SLR residuals of orbit prediction for the three strategies (unit: m).

Strategy	Fit length (h)		Eclipse	Non-eclipse
1	24	Mean	−0.556	−0.964
		STD	0.726	0.872
	48	Mean	−0.248	−0.333
		STD	0.491	0.291
	72	Mean	−0.288	−0.249
		STD	0.528	0.328
2	24	Mean	−0.399	−0.467
		STD	0.869	0.797
	48	Mean	−0.024	−0.041
		STD	0.419	0.213
	72	Mean	−0.039	0.002
		STD	0.497	0.257
3	24	Mean	−0.395	−0.453
		STD	0.869	0.801
	48	Mean	0.038	−0.088
		STD	0.221	0.226
	72	Mean	−0.094	0.018
		STD	0.233	0.242

3.3 SRP parameters

The SRP parameters were generated together with orbit predictions; these were estimated as constant parameters and updated once a day. The SRP parameter time series of BDS C01 satellite over 2 years and the elevation angle of the Sun above the satellite orbital plane are shown in Figure 10. SRP parameters of strategies with 48-h and 72-h fitting arc lengths showed similar performance in non-eclipse seasons. The estimated SRP parameters derived from the 24-h fitting arc length considerably differed from those estimated from the 48-h and 72-h fitting arc lengths, as only one orbital arc length data was used in the former estimations. Moreover, for solutions with the three types of fitting arc lengths, two jumps occurred in the D₀ time series over our experimental period. These jumps occurred during the north–south orbit maneuver period, in which the orbit inclination changes. In strategy 2, which included an a priori SRP model, the D₀ estimates decrease in the range of −1.2 nm/s² to −0.3 10⁻⁵ nm/s², because of the large a priori correction in D₀ direction.

Figure 10. Estimated SRP parameters (D₀,Y₀,B₀,B_C, and B_S) time series with fitting arc lengths of 24-h, 48-h, and 72-h for strategies 1–3 for BDS C01. Other BDS GEO satellites show similar performances.

The estimated SRP parameters and the elevation angle of the Sun are correlated. In eclipse seasons, especially when the elevation angle of the Sun approaches zero degrees, the values of the estimated D₀,Y₀,B₀,B_C, and B_S parameters exhibited considerable biases in the solutions of strategies 1 and 2 with 48-h and 72-h fitting arc lengths, whereas the 24-h fitting arc length provided more consistent results. This was possibly because the SRP model error is larger when the elevation angle of the Sun is close to zero; thus, a longer fitting arc length will introduce more inconsistencies to the SRP parameter estimates. In strategy 3, which includes the empirical acceleration in the along-track direction, the biases of estimated SRP parameters in eclipse seasons were considerably less than those in solution of strategy 2. Thus, the addition of the empirical constant-acceleration parameter improved the stability of the SRP parameters, especially the Y₀ parameter in eclipse seasons. In order to better explain the obvious variation of Y₀ parameter, we analyzed the correlations for estimated SRP parameters. The correlation coefficients are dependent on the β angle. In particular, Y₀ was strongly correlated with B_S when the absolute value of the β angle was close to zero. The introduction of empirical constant acceleration weakened the correlation between Y₀ and B_S when 48-h or 72-h fitting arc length was used, as displayed in Figure 11. When the β angle was nearly equal to zero degree, correlations have been significantly reduced. Consequently, strategy 3 yields improved orbit determination in eclipse seasons (Figure 7 and Figure 9).

Figure 11. Correlations coefficients between the Y₀ and B_S parameters of BDS C01. The black lines indicate the elevation angle of the Sun.

4. Discussion and conclusion

To improve the ultra-rapid orbit products of BDS GEO satellites, we investigated the influences of the fitting arc length, an a priori SRP model, and along-track empirical acceleration on orbit prediction. To determine which factors influenced the orbit determination solution, we designed three POD strategies for comparison: strategy 1 involves only ECOM5; strategy 2 combines ECOM5 with an a priori model; and strategy 3 also includes an along-track empirical constant-acceleration parameter. The 24-h predicted orbit arcs were evaluated through comparisons with the corresponding observed orbit arc and SLR observations.

Although different dynamic models were applied, the BDS GEO satellite orbit prediction with a 48-h fitting arc length was considerably better than those with 24-h and 72-h fitting arc lengths, in both eclipse and non-eclipse seasons. With a 48-h fitting arc length, strategy 2, with an a priori SRP model, yielded systematic biases for SLR residuals of −0.024 m and −0.041 m in eclipse and non-eclipse seasons, respectively, which were 90.3% and 87.7% less than those of strategy 1 (i.e. without an a priori SRP model). Similarly, the STD of the SLR residuals of strategy 2 under eclipse and non-eclipse conditions (0.419 and 0.213 m, respectively) were less than those of strategy 1 (0.491 and 0.291 m, respectively). However, the introduction of the a priori SRP model did not improve the overlap precision, because the overlap precision mainly reflects the consistency of the applied models rather than their absolute accuracy. These results also suggest that the a priori SRP model needs further improvement.

Furthermore, when 48-h or 72-h fitting arc length was used, the addition of an empirical acceleration (i.e. strategy 3) weakened the correlation between Y₀ and B_S and improved the stability of estimated SRP parameters in eclipse seasons. The overlap precision was 0.268 m and 1.080 m in the radial and 3D components in eclipse seasons, respectively, which were 52.5% and 63.6% less than the results obtained without empirical acceleration. Meanwhile, the STD of the SLR residuals of strategy 3 (0.221 m) was 47.0% less than that of strategy 2 (0.419 m). The orbit prediction performance in eclipse seasons was comparable to that in non-eclipse seasons. Overall, the use of a 48-h fitting arc length, an a priori SRP model, and along-track empirical acceleration yielded the optimal BDS GEO ultra-rapid orbit prediction products.

However, it should be mentioned that in the actual processing, the precision of ultra-orbit prediction will be affected by the unavoidable errors associated with the prediction of the Earth rotation parameters that is needed in the coordinate translation between the ECEF and ECI frames. Meanwhile, despite improvements of BDS GEO predicted orbits achieved by using strategy 3 in this study, the orbit precision is still worse than that of BDS Inclined Geosynchronous Orbit (IGSO) and MEO satellites, especially in the along-track direction. Thus, further investigation is needed. Efforts can be expected from further optimizing the SRP model, modeling the Earth radiation pressure, as well as introducing onboard observations from LEO or MEO satellites into the ultra-rapid processing of BDS GEO satellites.

Data available

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

The experiment data are collected from the MGEX, iGMAS networks. The numerical calculations in this paper have been done on the supercomputing system in the Supercomputing Center of Wuhan University.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [grant number:41904021].

Notes on contributors

Wenxi Zhao

Xiaolei Dai is currently an associate professor at Wuhan University. She received her Ph.D. degree from Wuhan University in 2016. Her current research interests are real-time precise orbit determination and clock estimation of GNSS, orbit attitude and force modeling for BDS and real-time data processing in the RTS system.

Xiaolei Dai

Wenxi Zhao is currently a postgraduate in GNSS Research Center, Wuhan University. Her area of research currently focuses on real-time GNSS orbit determination.

Yidong Lou

Yidong Lou is currently a professor at GNSS Research Center, Wuhan University. He received his Ph.D. in Geodesy and Surveying Engineering from the Wuhan University in 2008. His current research interest is in the real-time precise GNSS orbit determination and real-time GNSS PPP.

Yaquan Peng

Yaquan Peng is currently a Ph.D. candidate in GNSS Research Center, Wuhan University. His area of research currently focuses on multi-GNSS orbit determination.

Xueyong Xu

Xueyong Xu is a professorate senior engineer in North Information Control Research Academy Group Co. Ltd. China. He received his Ph.D. degree from University of Science and Technology of China. His current research focus on high-precision GNSS applications.