Accurate and Rapid Broadcast Ephemerides for Beidou-Maneuvered Satellites

Abstract

The geostationary earth orbit (GEO) and inclined geosynchronous orbit (IGSO) satellites of the Beidou navigation satellite system are maneuvered frequently. The broadcast ephemeris can be interrupted for several hours after the maneuver. The orbit-only signal-in-space ranging errors (SISREs) of broadcast ephemerides available after the interruption are over two times larger than the errors during normal periods. To shorten the interruption period and improve the ephemeris accuracy, we propose a two-step orbit recovery strategy based on a piecewise linear thrust model. The turning points of the thrust model are firstly determined by comparison of the kinematic orbit with an integrated orbit free from maneuver; afterward, precise orbit determination (POD) is conducted for the maneuvered satellite by estimating satellite orbital and thrust parameters simultaneously. The observations from the IGS Multi-Global Navigation Satellite System (GNSS) Experiment (MGEX) network and ultra-rapid products of the German Research Center for Geosciences (GFZ) are used for orbit determination of maneuvered satellites from Sep to Nov 2017. The results show that for the rapidly recovered ephemerides, the average orbit-only SISREs are 1.15 and 1.0 m 1 h after maneuvering for GEO and IGSO respectively, which is comparable to the accuracy of Beidou broadcast ephemerides in normal cases.

1. Introduction

During nominal operations, Beidou broadcast ephemerides including satellite orbit and clock information are generated and uploaded regularly by the ground segment of the Beidou navigation satellite system (BDS), with an update rate of every hour [1]. Each satellite broadcasts its ephemeris to the users to enable calculation of satellite position, velocity and clock error for arbitrary instants. However, in the case of a satellite maneuver, the satellite orbit cannot be well determined and predicted in the absence of exact knowledge of the thrust force.

Due to the equatorial ellipticity of the Earth, the 1:1 orbital resonance effect caused by the commensurable angular motion of geostationary earth orbit/ inclined geosynchronous orbit (GEO/IGSO) satellites and rotational velocity of the Earth can lead to periodic variation and secular evolution of the satellite Keplerian elements. Beidou GEO and IGSO satellites are more prone to drift away from their assigned orbit slots [2,3,4]. Moreover, the lunisolar resonance and solar radiation pressure (SRP) also have a perturbation effect for orbits at the geosynchronous distance [5,6,7]. Frequent maneuvers are executed for Beidou GEO and IGSO satellites, generally, about 9 to 15 times per year for each GEO and twice per year for each IGSO, whereas no maneuver was found for Beidou medium earth orbit (MEO) satellites in a one-year period studied [8]. In the vicinity of a maneuver, the Beidou satellite broadcast ephemerides are marked as unhealthy for several hours, and the accuracy of the healthy ephemerides available immediately after maneuvering is lower than that in normal cases.

In routine operation, the Beidou maneuvered satellite is removed from the system precise orbit determination (POD) process, and its orbit is determined separately with only observations after maneuvering [9]. Without using observations before and during the maneuvering process, it usually takes as long as 24 h for the orbit accuracy of the maneuvered satellite to arrive at a stable level [10]. Rapid orbit recovery for Beidou maneuvered satellites is crucial to improving the satellite usability and system performance.

The China Beijing Satellite Navigation Center (BSNC) has developed a fast short-arc orbit recovery method with prior information, i.e., SRP parameters and system bias, from the long-arc POD before maneuvering, and the clock errors estimated beforehand. Validation by satellite laser ranging (SLR) showed an orbit accuracy of 0.71 m in the radial direction within 4 h after the maneuver, and 1.17 m for 1 h prediction [10]. The maneuvered orbit can be effectively recovered with this method; however, a gap of at least 4 h still exists in the broadcast ephemeris. Rapid orbit recovery strategies with short-arc tracking data using kinematic orbit determination have also been proposed. Orbit determination with pseudo-ranges of 10 min has obtained a 3-dimensional (3D) accuracy of 3.27, 8.19 and 5.6 m for the GEO, IGSO and MEO satellites, respectively. By orbit fitting of the kinematic orbits based on a simplified Keplerian ephemeris model and then orbit prediction, the predicted orbits could maintain a user range error (URE) of better than 3 m for 10 min and 3–6 min for GEO and IGSO/MEO satellites, respectively [11]. Apart from its low accuracy, with this short-arc orbit determination method, the ephemerides of GEO and IGSO/MEO have to be updated every 10 min and 3–6 min, more frequent than the standard ephemeris update rate of every hour.

Long-arc orbit determination taking advantage of the dynamic information is the standard way of obtaining high-precision satellite orbit, especially in terms of orbit prediction. Efforts have been made towards continuous orbit determination of the maneuvered satellite with observations before, during and after maneuvering [12,13,14]. Corresponding thrust models have to be used in this case based on satellite control information which contains or can deduce the maneuver start and end time, satellite mass, pulse width and period, force efficiency, and so forth. The control information is available from the telemetry and telecontrol data or Xi’an satellite measurement and control center [12,13]. Orbit recovery for Beidou maneuvered satellites has been investigated with various thrust models, such as the widely adopted constant thrust model [12,13,14], and the piecewise constant thrust model complemented with a first-order Markov stochastic process [15]. Su has conducted GEO maneuvered satellite orbit determination in a simulation study and proposed a linear force model combined with a second-order polynomial [16].

Recent orbit recovery studies using constant thrust model by the institutions mainly in charge of Beidou satellite control or monitoring have shown good results. The C-band ranging data within the Chinese Area Positioning System (CAPS) have high accuracy and do not contain satellite and station clock errors. The maximum of the thrust force in the radial along- and cross-track (RAC) directions can be obtained from the telemetry data. By supposing a sinusoid change of the thrust force, an averaged constant vector average = max × 2/π defined in the maneuvering directions may be used to represent the thrust. However, this force vector contains errors and cannot enable high-precision satellite orbit determination. Therefore, with the maneuver start and end time from the telemetry data and C-band observations, constant force parameters in the RAC directions during the maneuvering period were estimated by BSNC. The obtained 3D and radial orbit accuracies were around 30 m and 3 m 2 h after the maneuver, 20 m and 2.5 m 4 h after [12]. Using C-band observations one day before and 2 h after maneuvering as well as the pulse information, continuous orbit determination of the GEO satellites achieved a predicted URE of less than 3 m 3–4 h after maneuvering [13]. Though the orbit recovery performance has been dramatically improved, the UREs of predicted orbits in both studies [12,13] exceed the standard positioning requirement of 2 m [17] even several hours after maneuvering.

Evaluation of the normal Beidou broadcast ephemerides from February 2014 to October 2016 indicated an accuracy of 0.85, 9.23 and 3.57 m in the RAC directions for GEO, and 0.66, 2.71 and 1.75 m for IGSO [18]. The orbit-only URE of the whole constellation achieved an accuracy of 0.8 m in the time frame of 2013–2016 [19]. The recent analysis of normal broadcast ephemerides in 2017 showed that the signal-in-space ranging errors (SISREs), including orbit, clock and bias errors, are around 1 m for all the satellites, and 0.7 m while excluding GEO satellites [20]. Compared to the accuracy of normal broadcast ephemerides, the predicted orbits after maneuvers in the studies discussed above are not accurate enough and should be further improved.

Previously, we proposed a kinematic orbit determination combined with orbit integration method for maneuver detection and thrust force estimation. By comparing the kinematic orbit determined using ground observations and the predicted orbit considering all the possible dynamic forces but the thrust force, we could obtain the thrust-induced velocity changes as well as the accelerations. We found that, typically, the thrust forces increase steadily at the beginning of the maneuver, then stay nearly constant for around 20 min, and come back to zero abruptly at the end of the maneuver [21]. The maneuver behavior has been explicitly investigated in the previous study; however, a mature thrust model has not been set up, and the orbit has not been recovered by an accurate dynamical orbit determination method.

Here, according to this thrust-changing tendency revealed in the previous study [21], a piecewise linear function is proposed to model the thrust forces in each of the RAC directions, so as to realize rapid orbit recovery after maneuvering. We begin with an introduction of the thrust model and the detailed orbit determination strategy. Subsequently, the accuracy evaluation methods for our recovered ephemeris and the current healthy broadcast ephemeris after maneuvering are described, followed by the evaluation results of both ephemerides. A short conclusion is given in the end.

2. Methodology

To provide broadcast ephemeris for the maneuvered satellite, high-precision dynamic orbit determination uniting observations before, during, and after maneuver should be conducted. A general piecewise linear model has been proposed to model satellite thrust forces.

2.1. Thrust Force Model

The thrust forces in the RAC directions are estimated in three segments, as shown in Figure 1: the initial stage of the maneuver, the main period of the maneuver and the final phase of the maneuver, which are defined by four turning points t₀, t₁, t₂, and t₃ in the thrust model. The method to determine these points is described in the next subsection. To compensate for the discrepancies between the assumed and real force models and thrusting epochs, we suppose that the forces at t₀ and t₃ are not zero. Therefore, four constant force vectors $F_0,F_1,F_2$ and $F_3$ at these turning points are estimated. Although the thrust forces are mainly in the along-track direction for an in-plane maneuver and along- and cross-track directions for an out-of-plane one [8], it is not easy to judge whether there are any forces in other directions without knowledge of the telemetry and telecontrol data or the location of thrusters onboard the satellite. The telemetry and telecontrol data have indicated that the thrust forces can be in all of the three directions [12,22]. Thus, we estimate the components in the RAC directions of each constant force vector, and 12 thrust parameters are estimated by the proposed model.

The thrust force out of the maneuvering period is zero, and the force within [t₀, t₃] can be expressed as:

$$\mathbf{F}=\{\begin{array}{l}{\mathbf{F}}_0+({\mathbf{F}}_1-{\mathbf{F}}_0)\frac{t-t_0}{t_1-t_0},t_0\le t<t_1\\ {\mathbf{F}}_1+({\mathbf{F}}_2-{\mathbf{F}}_1)\frac{t-t_1}{t_2-t_1},t_1\le t<t_2\\ {\mathbf{F}}_2+({\mathbf{F}}_3-{\mathbf{F}}_2)\frac{t-t_2}{t_3-t_2},t_2\le t\le t_3\end{array}$$

(1)

These force parameters $F_i=(f_{ri},f_{ai},f_{ci}),(i=0,1,2,3)$ are estimated together with satellite initial orbital state elements and other dynamic parameters, such as the SRP. Similarly, the partial derivatives of the observations to the force parameters should also be obtained to set up the observation equation. The orbit related part in the linearized observation equation can be represented by:

$$\frac{\partial \rho }{\partial (\mathbf{r},\dot{\mathbf{r}})}\frac{\partial (\mathbf{r},\dot{\mathbf{r}})}{\partial {\mathbf{X}}_0}\Delta \mathbf{X}$$

(2)

where ρ is the range or range rate observation, $(r,\dot{r})$ denotes satellite position and velocity vector; $X_0=(r_0,{\dot{r}}_0,p_1,p_2,\dots ,p_n,f_{r0},f_{a0},f_{c0},\dots ,f_{r3},f_{a3},f_{c3})$ is the vector containing these to be estimated orbital parameters, in which $(f_{ri},f_{ai},f_{ci}),(i=0,1,2,3)$ are the thrust parameters and $p_1,p_2,\dots ,p_n$ represent other dynamic parameters; $\Delta X$ is the satellite orbital correction vector. The satellite Keplerian elements $\left(a,e,i,\Omega ,\omega ,u\right)$ are used to replace $\left(r_0,{\dot{r}}_0\right)$ within $X_0$ in some orbit determination software packages, such as Bernese global navigation satellite system (GNSS) software [23]. The partial derivatives of the observation $\rho$ with respect to satellite position and velocity $(r,\dot{r})$ can be deduced from the observation geometry, and $(r,\dot{r})$ can be known from integration of the satellite motion equation; the partial derivatives of $(r,\dot{r})$ to $X_0$ should be derived by resolving the initial value problem of the so-called variation equations. The detailed resolution process for maneuver-free cases can be found in the literature [24,25].

Assuming a unit satellite mass of m = 1 kg, the partial derivatives of the force vector $f=\ddot{r}$ to the explicit maneuver parameters at different maneuver stages can be expressed. Below, the non-zero partial derivatives of $f$ with respect to thrust parameters in the along-track direction during each stage of maneuver are listed as an example.

$$\begin{array}{l}\frac{\partial \mathbf{f}}{\partial f_{a0}}=(1-\frac{t-t_0}{t_1-t_0}){\mathbf{e}}_a(t),t_0\le t<t_1\\ \\ \frac{\partial \mathbf{f}}{\partial f_{a1}}=\{\begin{array}{l}\frac{t-t_0}{t_1-t_0}{\mathbf{e}}_a(t),t_0\le t<t_1\\ (1-\frac{t-t_1}{t_2-t_1}){\mathbf{e}}_a(t),t_1\le t<t_2\end{array}\\ \\ \frac{\partial \mathbf{f}}{\partial f_{a2}}=\{\begin{array}{l}\frac{t-t_1}{t_2-t_1}{\mathbf{e}}_a(t),t_1\le t<t_2\\ (1-\frac{t-t_2}{t_3-t_2}){\mathbf{e}}_a(t),t_2\le t\le t_3\end{array}\\ \\ \frac{\partial \mathbf{f}}{\partial f_{a3}}=\frac{t-t_2}{t_3-t_2}{\mathbf{e}}_a(t),t_2\le t\le t_3\end{array}$$

(3)

where $e_a(t)$ is the unit vector in the along-track direction of the inertial system. For the partial derivatives to thrust parameters in the other two directions, we shall replace $e_a(t)$ with $e_r(t)$ or $e_c(t)$. It can be deduced that the partial derivative matrix of $(r,\dot{r})$ with respect to the thrust parameters is a zero-matrix before maneuvering, i.e., $t\le t_0$, and non-zero during and after the maneuver, i.e., $t>t_0$.

2.2. Two-Step Orbit Recovery Strategy

To carry out the POD of the maneuvered satellite with the proposed thrust model, the turning points of different maneuver stages have to be estimated. These points can be obtained by comparing the kinematic orbit with an integrated one. The kinematic orbit is determined with ground tracking observations and represents the actual trajectory of the satellite; whereas the orbit integration is conducted considering all the possible perturbations but the thrust force. The difference between these two orbits can show the effect of satellite maneuver to some extent. Therefore, we adopt the two-step method of firstly determining the turning points and then conducting dynamic orbit determination with the thrust model.

2.2.1. Determination of the Thrust Turning Points

Kinematic orbit determination should be carried out following a method similar to that used for maneuver detection; however, the initial orbit of the maneuvered satellite prior to maneuvering period is more precise than a previous study [21]. The orbit determined in an arc containing one day preceding the maneuver and the period before maneuvering on the maneuver day is adopted as the initial orbit, rather than the previously used broadcast ephemerides. With the accurate initial orbit of the maneuvered satellite at the epochs before maneuvering, ultra-rapid orbits of other satellites, known station coordinates, and different measurement corrections applied, the ambiguities of the double difference ionosphere-free observations can be fixed to the average of the ones calculated in the maneuver-free period in the kinematic orbit determination if no cycle slip happens. The kinematic orbit in the period during and after maneuvering can achieve a cm-dm level accuracy in the along- and cross-track directions. Due to the poor geometry and observation noise, the orbit error in the radial direction can jump to several meters [21].

By extrapolation of the orbit determined before maneuvering without considering thrust forces, the reference orbit during the maneuver period can be obtained. Satellite velocity and acceleration changes due to maneuver can be extracted by comparing the reference orbit with the kinematic one. Due to the observation noise, the extracted velocity and acceleration changes are very noisy with large variations or jumps (as shown in Figure 2) and, therefore, the turning points t₀, t₁, t₂, and t₃ are difficult to be directly identified from the extracted accelerations. In this study, the following method has been adopted for turning points determination. During the initial stage of maneuver [t₀, t₁], the acceleration changes are quite small; the extracted values still vary around zero and, therefore, cannot be used to identify t₀. On the other hand, the velocity changes can start to show the effect of the thrust by increasing or decreasing in a particular direction. Therefore, t₀ is estimated using a 2-min sliding window in the velocity changing domain. Considering the observation noise, t₀ is determined as the epoch 2 min before the first sliding window in which the velocity changes are of the same sign. During the main maneuver period [t₁, t₂], the accelerations have explicit bias from zero and t₁ is determined as the first epoch of the first 2-min window in which all the accelerations are negative or positive, and t₂ is the last epoch of the last window satisfying the same sign condition. The last points, t₃, is set as the epoch 3 min after t₂, which is set empirically based on our data-processing results.

2.2.2. Precise Orbit Determination (POD) of the Maneuvered Satellite

As analyzed in the last subsection, the partial derivatives of ($r$, $\dot{r}$) with respect to the thrust parameters are not zero since the start of maneuver, and the measurements observed since then will contribute to the estimation of thrust parameters [26]. Meanwhile, to better estimate other unknown parameters, the orbit should be determined in a long arc and measurements before maneuvering should also be taken into account. It has been shown that comparatively better orbit recovery performance can be achieved when observations from 24 h before maneuvering to 2 h after were used [13]. In this study, we determine the orbit of the maneuvered satellite in a 24-h sliding window, including data before, during and after the maneuver, and then carry out orbit prediction for further 1 h to obtain the broadcast ephemeris.

As is shown in Figure 3, assuming that the maneuver happens during 9:35–9:55 on Day 2, as indicated by the red shadow, the first POD process can be conducted using observations from 11:00 of Day 1 to 11:00 of Day 2. The initial orbit for the maneuvered satellite used in this first sliding window is a concatenation of the precise orbit before maneuvering and kinematic orbit during and 1 h after the maneuver. With orbit determined in the sliding window, the orbit for the hour following the window can be predicted and uploaded at epoch 11:00 of Day 2. Allowing for the time needed for the POD, orbit prediction, and fitting, we conduct the POD in an earlier time, using data from 10:50 of Day 1 to 10:50 of Day 2. With incoming observations and the orbit obtained from the previous window as initial orbit, the POD can be continued in the sliding window that moves forward by a step of one hour, as shown by the second row of Figure 3. Orbit determined in the second window will be used to generate ephemeris at epoch 12:00 of Day 2. This process will be repeated with the sliding window stepping forward hour by hour.

The orbit integration algorithm should be selected carefully in the case of maneuvered satellite orbit determination. Due to the large magnitudes of the thrust, mismatch of the turning points and the integral nodes will severely deteriorate the orbit accuracy. A modified Adams–Cowell integrator that combines the multi-step and single-step method around the maneuver has been proposed and applied for the LEO satellite orbit determination [26]. Although the Beidou thrust is not continuous and usually characterized as pulses [13], it can last 20 min, and an averaged constant force is generally assumed. In our study, with observation sampling intervals of 30 s, the maneuver timing uncertainty is 30 s. Therefore, we shorten the integration step-size, with an interval of 30 s for the motion equation and 1 min for the variation equations.

The Bernese GNSS software is modified for the maneuvered satellite POD in this research. The reference frames and dynamic models are listed in

Table 1. Reference frames and dynamic orbit models in POD.

Reference Frames
Inertial frame	International celestial reference frame (ICRF) at J2000.0
Terrestrial frame	International terrestrial reference frame (ITRF) 2014 [32]
Precession model	International astronomical union (IAU) 2000 Precession theory
Nutation model	IAU 2000R06 Nutation theory [33]
Orbit models
Geo-potential (static)	Earth gravitational model 2008 up to degree and order 12 (+C21+S21) [34]
Tidal variations in geo-potential	Solid earth tides: International earth rotation and reference systems service (IERS) 2010 [33]
	Ocean tides: Finite element solutions (FES) 2004 model [35]
	Solid earth and oceanic pole tide: IERS 2010
Third-body	Sun, Moon, Jupiter, Venus, Mars as point masses
Ephemeris used	Jet Propulsion Laboratory (JPL) DE421 [36]
SRP model	Reduced ECOM and a constant along-track acceleration [37]
Relativistic effects	IERS 2010
Thrust forces	Proposed piecewise thrust model
Attitude model	Orbit normal (ON) for GEO, yaw steering (YS) and ON for IGSO and MEO [38]

[23]. The double differenced ionosphere-free linear combination of B1 and B2 phase and code observations are used in the POD process to eliminate the satellite and station clock errors. Only data files which contain observations of the maneuvered satellite are used, and the orbits of other satellites are not estimated but fixed to the ultra-rapid orbits of GFZ. It improves the computation efficiency and furthermore avoids the possible deterioration of other orbits by the maneuvered satellite. Considering the comparatively low accuracy of GEOs, the double difference is formed between the maneuvered satellite and IGSO/MEO satellites. The earth orientation parameters (EOPs) are also from GFZ ultra-rapid products. The station coordinates and tropospheric delays are estimated beforehand using the Globe Mapping Function (GMF) [27] and the gradient model proposed by Chen and Herring (CHENHER) [28] in a precise point positioning (PPP) process where the ultra-rapid orbits are also used. Phase center offsets (PCOs) of some Beidou satellites are obtained from the estimated values [29], and the nominal value of (0.6, 0.0, 1.1 m) is adopted for others; satellite phase center variation (PCV) corrections are not applied. The ground antenna PCO and PCV of Beidou satellites are assumed the same as those of GPS from igs14.atx [30]. The estimated orbital parameters are the 6 Keplerian elements $a,e,i,\Omega ,\omega ,u_0$, 5 SRP parameters of the reduced Empirical CODE Orbit Model (ECOM), a constant acceleration in the along-track direction [31], and 12 thrust parameters of the proposed model.

2.3. Accuracy Evaluation Methods

The broadcast ephemerides are generated following the process of orbit determination, orbit prediction, and orbit fitting [39]. The final step, i.e., orbit fitting to estimate the elements of broadcast ephemerides, has been widely studied for normal satellites and the orbit fitting accuracy is on an mm–cm level [40,41,42]. Orbit fitting of 2 h, 1 h of measured orbit and 1 h of predicted orbit, showed better than 0.04 m position accuracy for GEO and IGSO satellites, and 0.02 m for MEO using the traditional 16-element ephemeris model [40]. With the proposed 16- and 18-element non-singular ephemeris models, the root mean squares (RMSs) of the orbit fitting UREs are less than 0.05 m and 1 mm for GEOs in a 2-h arc, respectively [41,42]. In theory, the equivalent fitting accuracy can be obtained for maneuvered satellite as long as it is conducted in an arc free from thrust, such as the first 2 h after maneuvering with 1 h of measured orbit and 1 h of predicted orbit. Considering the high orbit-fitting accuracy, we focus on the first two steps of broadcast ephemeris generation, and the 1-h predicted orbit after the orbit determination window represents the broadcast ephemeris. The accuracy of the ephemeris is evaluated by comparing with the precise orbit determined in an arc after the maneuver.

Before evaluating our fast-recovered ephemeris, we firstly analyze current healthy broadcast ephemerides available several hours after the maneuver from BDS, which gives a baseline for comparison. All Beidou satellite maneuvers of one year are studied, and the accuracy is evaluated by comparing the ephemerides to the reference orbits, which are obtained by orbit-fitting of the available precise ephemerides on the two successive days after maneuver and then backward orbit extrapolation to the time just after the maneuver. The backward predicted orbit in the period after maneuvering on the maneuver day has an accuracy of 0.6, 3.6 and 0.7 m in the RAC directions for GEO, and 0.3, 0.6 and 0.5 m for the IGSO, which is obtained by comparing the predicted orbits with accurate orbits determined after maneuvering for all the maneuvering cases in October 2017.

3. Results and Discussion

We first detected all the Beidou satellite maneuvers in 2017.

Table 2. Beidou geostationary earth orbit (GEO) satellite maneuver dates and broadcast ephemeris unhealthy periods in 2017.

	C01			C02			C03			C04			C05
January	09	05–10	6 h	25	09–14	6 h	12	08–13	6 h				05	00–05	6 h
February	13	10–17	8 h				06	08–13	6 h				10	00–05	6 h
	14	03–14	10 h
March	21	00–16	11 h	01	09–14	6 h	04	08–13	6 h	07	08–13	6 h	17	01–05	5 h
							30	08–13	6 h
April	21	05–10	6 h	10	09–14	6 h	25	08–13	6 h				18	00–05	6 h
May	22	05–09	5 h	16	09–15	7 h	19	08–14	7 h	31	02–23	8 h	25	00–04	5 h
June	23	09–14	6 h	27	09–14	6 h	14	08–13	6 h	01	06–11	6 h
July	07	09–14	6 h				11	08–13	6 h				01	00–05	6 h
August	28	09–14	6 h	04	10–02	17 h	07	09–00	13 h				10	00–05	6 h
September	30	09–14	6 h	11	09–14	6 h	01	09–14	6 h	04	10–14	5 h	18	00–05	6 h
							27	09–14	6 h
October	31	09–14	6 h	19	09–14	6 h	23	09–14	6 h				09	22–04	7 h
													10	23–04	6 h
November	30	09–14	6 h	20	19–02	8 h	17	10–14	5 h				10	00–05	6 h
				21	10–16	7 h
December	27	09–14	6 h	08	10–14	5 h	11	09–14	6 h	19	07–12	6 h	15	00–05	6 h

and

Table 3. Beidou inclined geosynchronous orbit (IGSO) satellite maneuver dates and ephemeris availability in 2017.

	First Maneuver			Second Maneuver
C06	April 13	22–03	6 h	September 24	10–14	5 h
C07	April 27	23–03	5 h	October 16	10–15	6 h
C08	January 22	13–18	6 h	August 01	01–05	5 h
C09	March 27	01–05	5 h	September 21	10–16	7 h
C10	June 09	00–05	6 h	November 27	10–16	7 h
C13	February 20	13–21	9 h	August 17	01–07	7 h

list the maneuver dates and broadcast ephemeris abnormal periods for Beidou GEO and IGSO satellites separately.

In Table 2, each of the three sub-columns for each GEO satellite represents the day of the month, the period of unhealthy broadcast ephemerides, and the total number of hours during this period, respectively. On average, the GEO satellites are maneuvered about once per month and the broadcast ephemeris unhealthy periods usually range around 5–8 h and can be as long as 17 h for some individual cases. Two maneuvers are detected for each of the IGSO satellite and are shown explicitly in Table 3.

3.1. Healthy Broadcast Ephemerides

The accuracy of the healthy broadcast ephemerides after maneuvering was evaluated with all the maneuvering cases in 2017.

Table 4. Root mean square (RMS) errors of the first hour of healthy broadcast ephemerides after maneuvering in 2017 for each Beidou satellite.

		Radial (m)	Along (m)	Cross (m)	3d (m)
	C01	4.51	15.08	11.83	19.69
	C02	1.54	8.01	9.54	12.55
GEO	C03	2.66	5.64	10.96	12.61
	C05	3.38	8.50	25.15	26.76
	Mean	3.02	9.31	14.37	17.90
	C06	1.42	8.63	3.50	9.42
	C07	1.50	11.09	8.19	13.87
IGSO	C08	2.67	7.65	11.48	14.05
	C09	0.91	2.70	1.78	3.36
	C10	0.78	1.52	5.58	5.83
	C13	0.77	7.71	3.55	8.53
	Mean	1.34	6.55	5.68	9.18

lists the RMS errors of the healthy broadcast ephemerides in their first available hour after maneuvering of all the maneuvering cases for each satellite, except for C04 due to the absence of precise orbit products needed for accuracy evaluation. It can be seen from the table that, apart from the unavailability of 5–8 h after the maneuver, the accuracy of the healthy ephemerides available later is still poorer than the normal cases of 0.85, 9.23 and 3.57 m in the RAC directions for GEO, and 0.66, 2.71 and 1.75 m for IGSO [18].

Figure 4 shows the average orbit-only SISREs of all the GEO and IGSO broadcast ephemerides in each hour after maneuvering in 2017. Ephemerides with gross errors are not taken into account. It can be observed from the figure that the SISREs of GEO are larger than 2 m from the 6th to 10th h after maneuvering and a decreasing tendency is displayed over time in the first 24 h. The dependence on time is not apparent for SISREs of IGSO, and they vary in a range of 1.3 to 2.5 m in the 30 h.

3.2. Orbit Recovery Results

In this paper, to test our orbit recovery method, we implemented orbit recovery for the IGSO maneuvers happened from September to November 2017 and GEO maneuvers in October and November 2017, using daily observations from the MGEX network and ultra-rapid products of GFZ. Approximately 50 stations near the longitude of each GEO satellite were used for its POD, and tracking stations of the IGSO satellite have a superior worldwide location, with a total number of more than 80 (Figure 5).

3.3.1. Maneuver Characteristics

The velocity changes due to maneuvers of the studied cases are listed in

Table 5. Velocity changes of Beidou maneuvered satellites in September–November 2017.

Satellite-Day	dVR (m/s)	dVA (m/s)	dVC (m/s)
C05- October 09	0.2120	0.3458	−28.9157
C05- October 10	−0.0129	−0.2205	−0.0015
C02- October 19	−0.0129	−0.0730	−0.0065
C03- October 23	−0.0089	−0.1466	−0.0025
C01- October 31	−0.0138	−0.1215	−0.0113
C05- November 10	0.0035	0.1282	0.0003
C03- November 17	0.0003	0.1202	0.0030
C01- November 30	0.0136	0.1023	0.0095
C09- September 21	0.0197	0.3929	0.1254
C06- September 24	0.0217	0.5311	0.0844
C07- October 16	−0.0305	−0.5638	0.4951
C10- November 27	0.0155	0.3763	0.2390

. The IGSO maneuvers were executed mainly in the along- and cross-track directions, with velocity increase or decrease of several dm/s. GEO maneuvers are in-plane, and velocities change for about 0.1–0.2 m/s in the along-track direction, except the huge out-of-plane one on 09 October. The GEO out-of-plane maneuver is more complicated and happens at a rather low frequency, 4 out of the nearly 70 maneuvers detected in 2017, and we did not perform orbit recovery for it in this study.

The turning points of the thrust model have been determined according to the proposed method, and they are given in

Table 6. Turning points of the studied maneuvering cases in September–November 2017.

Satellite-Day	t0	t1	t2	t3
C05- October 10	23:21:30	23:23:30	23:57:30	00:00:30
C02- October 19	09:16:30	09:39:30	09:51:00	09:54:00
C03- October 23	09:09:00	09:29:30	09:50:30	09:53:30
C01- October 31	09:10:00	09:36:00	09:54:30	09:57:30
C05- November 10	00:13:00	00:37:30	00:52:30	00:55:30
C03- November 17	09:28:30	09:30:30	09:47:00	09:50:00
C01- November 30	09:35:30	09:37:30	09:53:00	09:56:00
C09- September 21	09:32:30	09:35:00	09:45:00	09:48:00
C06- September 24	09:22:30	09:23:30	09:37:30	09:40:30
C07- October 16	10:19:30	10:21:00	10:39:30	10:42:30
C10- November 27	09:22:30	09:24:30	09:36:00	09:39:00

. With these turning points, the dynamic orbit determination was carried out within a 24-h sliding window using data before, during, and 1 h, 2 h, …, up to 8 h after the maneuver, respectively.

3.3.2. Orbit Errors in Radial Along- and Cross-Track (RAC) Directions

The 1-h predicted orbits after the sliding window were compared with the precise orbits. We concatenate the orbit differences of each hour and plot them in Figure 6 and Figure 7 for some satellites. The accuracy of healthy broadcast ephemerides available in the recovery periods was also displayed in the same figures for comparison.

Figure 6 shows the accuracies of orbit recovery results and healthy broadcast ephemerides for the GEO satellites in October 2017. The recovered GEO errors are generally below 2 m in the radial direction, except for the comparatively sizeable radial orbit error of 2–4 m for C01 and C02 during 4–6 h. The orbit errors in the along- and cross-track directions have been significantly reduced, less than 4 and 1 m, respectively. The improvement of orbit accuracy in the along- and cross-track directions can be benefited from the accurate estimation of the thrust and fixing of other unknown parameters in the POD process. The orbit recovery results are much better than the healthy broadcast ephemerides during 6–9 h. The orbit of satellite C05 is not successfully determined in the first sliding window, which may due to the effect of the previous out-of-plane maneuver happened at the end of 09 October.

Figure 7 shows the accuracy analysis of the IGSO satellites. We have noticed that although the broadcast ephemerides are healthy at epoch 9:00 for C06 on 24 September, C09 on 21 September and C10 on 27 October, the maneuvers actually happen between 9:00 and 10:00. The unavailable periods of the ephemerides are actually 1 h more than those shown in Table 3. It can be seen that despite the large radial errors of 3–4 m for C09 in the first several hours, the recovered orbits generally have an accuracy superior to 2 m in the radial direction. The recovered ephemerides are more accurate than the broadcast ephemerides 6 h or 8 h after maneuvering as well. The RMSs of the recovered orbit errors in the RAC components are below 1.55, 1.62 and 0.35 m, respectively.

3.3.3. Orbit-only Signal-in-Space Ranging Errors (SISREs)

The average orbit-only SISREs of these recovered and broadcast orbits in each hour after maneuvering are demonstrated in Figure 8. It shows that the orbit-only SISREs of both GEO and IGSO are less than 2.2 m throughout the recovery period; the average SISREs are 1.15 and 1.0 m for GEO and IGSO respectively. Furthermore, as indicated by the last three groups of bars, in the hours that the healthy broadcast ephemerides are available, the average SISREs of the recovered GEO and IGSO are 0.84 and 0.45 m, respectively. The recovered orbit accuracy is about three times better than the broadcast ephemerides in the common analysis period.

To analyze the impact of the update interval on the recovered orbits, further experiments have been done with update intervals of 0.5 and 2 hours, respectively.

Table 7. Average orbit-only SISREs of the recovered GEO and IGSO with different update intervals.

	GEO SISRE (m)			IGSO SISRE (m)
Time	0.5-hour	1-hour	2-hour	0.5-hour	1-hour	2-hour
2 h	1.72	2.11	2.11	0.93	1.05	1.05
3 h	0.98	1.08	1.69	1.24	1.35	1.98
4 h	1.08	1.19	1.25	1.87	2.03	2.03
5 h	1.13	1.20	2.12	1.18	1.27	2.51
6 h	0.83	0.92	0.88	0.46	0.47	0.47
7 h	0.82	0.89	1.31	0.52	0.54	0.61
8 h	0.63	0.71	0.91	0.31	0.35	0.35

lists the average orbit-only SISREs of the recovered orbits in each hour after maneuvering in different updating cases. It shows that the overall orbit accuracy improves when the update interval is shortened, and vice versa, especially during the 2–5 h after maneuvering. Shorter update intervals, i.e., 0.5 or 1 hour, should be preferred if we would like to recover the orbit soon after maneuvering.

Although the data were post-processed in this study, the POD algorithms can be adopted for real-time processing if real-time GNSS data are available. With the increasing number of real-time data streams available at the MGEX NTRIP (Networked Transport of RTCM via Internet Protocol) caster of BKG (Bundesamt für Kartographie und Geodäsie) [43], our method could be used for real-time ephemeris generation and the ephemeris could be accessed through the internet.

As the data we used has good global coverage, it is unfair to compare the recovered orbit accuracy with Beidou broadcast ephemeris directly, as they are generated mainly from tracking stations in China. Further studies are needed to evaluate the accuracy of the POD using our method with tracking stations similar to BDS tracking network.

4. Conclusions

Beidou GEO/IGSO satellites are frequently maneuvered to keep their positions in the designed orbits. Due to the unknown thrust forces, the Beidou healthy broadcast ephemerides are generally not available in the 5–8 h after maneuvers. Also, the accuracy of the ephemeris available after the gap is much worse than that in normal cases.

In this paper, we proposed a new Beidou POD algorithm which includes the maneuver periods, and the predicted orbit can be obtained one hour after the maneuver with much higher accuracy than the current Beidou broadcast ephemeris 5–8 hours later. A piece-wise linear model is proposed for modeling the thrust forces in the radial, along- and cross-track directions. The turning points of the models are estimated through the comparison of the integrated orbit without including the thrust forces and the kinematic orbit. The parameters of the thrust models are estimated together with other orbital parameters, such as the initial position and velocity, and the SRP parameters.

Using MGEX network GNSS data during September–November 2017, we have demonstrated that accurate and rapid orbits can be obtained after one hour of the maneuvers. The average orbit-only SISREs of the recovered ephemerides are 1.15 and 1.0 m for GEO and IGSO respectively, which is comparable to the accuracy of Beidou orbits in normal operation. Both the availability and accuracy of the broadcast ephemeris have been improved with our proposed orbit recovery method.

In this study, we have mainly focused on GEO in-plane maneuvers and IGSO maneuvers. For GEO out-of-plane maneuvers, the thrust forces are more complicated and the maneuver periods are much longer [21]. Further studies will be carried out to model the thrust forces during the out-of-plane maneuvers.