Recover the abnormal positioning, velocity and timing services caused by BDS satellite orbital maneuvers

The BeiDou Navigation Satellite System (BDS) provides global Positioning, Velocity, And Timing (PVT) services that are widely used in various areas. The BDS satellites frequently need the orbit maneuvers due to various perturbations to keep satellites in their designed positions. During these maneuvers, PVT services may be abnormal if the data from a maneuvering satellite is used. In this paper we developed an approach to recover the abnormal PVT services. By using BDS observations from multiple tracking stations, the orbital errors of a maneuvering satellite can be in real time obtained and corrected, thereby avoiding any influence on the performance of PVT services. The tests show that the average precision of position, velocity and timing services are improved by 0.8 m, 0.1 mm/s and 0.16 ns, respectively, using the developed orbital maneuver recovery approach. In addition, the approach can also be used for the orbital maneuver detection and monitoring.


Introduction
A navigation satellite system is an important utility, which can provide precise time and position information, and has a wide spectrum of applications (Hofmann-Wellenhof et al., 2008;Yang et al., 2011). For both global and local navigation satellite systems, the basic observable is the radio range measurement. The satellite position can be computed at any time using the broadcasts orbit ephemerides or precise orbit product. Given the satellite positions and ranging measurements with the satellite signals, the user position can be computed when four or more satellites are observed (Camacho-Lara, 2013). The positioning service is therefore based on the assumption that the satellites, signals, and receivers are in normal operation. The error of the satellite's position will bias the user's position. Thus, an abnormality in the broadcast ephemeris will have a serious influence on the user's Positioning, Velocity, And Timing (PVT) performance (Pervan & Chan, 2003).
Generally, abnormalities are caused by satellite orbital maneuvers, for which satellites use propulsion systems to change their orbits (Chen et al., 2009). Due to the effects of the non-spherical gravity field of the Earth and other perceptual factors, the orbital elements of a satellite will undergo long-term perturbations and its position will have large offset from the design position. In order to keep the satellites moving in their designed orbits, orbital maneuvering is needed, especially for Geostationary Earth Orbit (GEO) and Inclined Geosynchronous Orbits (IGSO) satellites (Cao et al., 2014;Shi et al., 2013;Wang et al., 2016). Usually, the offset of the satellite's position can reach to tens of kilometers from its normal position during the orbit maneuvering. This large offset of the satellite orbit will have serious impacts on the positioning, navigation, and timing service (Byun, 2003;Steigenberger et al., 2013;Du et al., 2015;Liu et al., 2017;Zhao et al., 2015;Geng et al., 2016). The usual approach to solve this problem is to broadcast the orbital maneuver information in the ephemeris. Users can then avoid using the satellite which is marked as maneuvering (Huang et al., 2018;Tu et al., 2020;Ye et al., 2017). Here, we examine whether the observation data at ground stations can be used to measure the orbital errors caused by BeiDou Navigation Satellite System (BDS) satellite maneuvers. We then explore whether these data can be used to correct the orbit and recover the normal PVT services.
After the introduction, section 'Methodology' introduces the methodology, followed by section 'Validation and results' about the data used and the detailed validations of the results. Then section 'Discussions' discusses the results and the concluding remarks are in section 'Conclusions' .

Estimation of the epoch variation of the receiver clock error by carrier phase observations
The epoch-differenced velocity estimation model is used for the estimation of the epoch variation of the receiver clock error, which can be expressed as follows (Colosimo et al., 2011;Tu et al., 2016Tu et al., , 2020: where represents the epoch difference; is the wavelength of the carrier phase observation; φ represents the carrier phase measurement; the superscripts s and r refer to satellite and receiver respectively; A is the designed coefficient matrix; v is the receiver's three Dimension (3D) velocity; t s and t r represent the satellite clock error and receiver clock error, respectively; ρ is the geometric distance from the satellite to the receiver; M represents the modeled error, which contains the Earth's rotation, relativity effect, the solid Earth tide, ocean tide loading, the ionosphere error, tropospheric error, and phase wind-up; U represents the un-modeled errors such as the atmospheric residual, ephemeris residual, and multipath error; and ε represents the measurement noise (Tu et al., 2016).
The ionosphere error is processed by using the dualfrequency ionosphere-free observation, and the other modeled errors are corrected using the proper models (Kouba & Pierre, 2001). The un-modeled errors are then absorbed by the estimated parameters. (1) The stochastic model is determined by the measurement error and satellite elevation, which is expressed as follows (Xu et al., 2013): where P is the weight matrix; e represents the satellite elevation angle (unit: rad), a represents the measurement error, whose value is set to be 0.002-0.003 m for the carrier phase measurements.
By combining Eqs. (1) and (2), the least squares method can be used to estimate the parameters while more than four satellites are tracked. The unknown parameters are the 3D velocities and the epoch-differenced receiver clock error (Tu et al., 2016). Moreover, the 3D velocities are constrained to zero values as the station is static and the observation of the monitored satellite m should be removed for the parameter estimation. Then, the epoch variation of the receiver clock error for each monitoring station can be obtained, and will be used in the subsection 'Estimation of the dynamic variation of the orbital maneuver with multi-station data' .
In addition, while the epoch-differenced observation model is used for the velocity estimation, the cycle slip must be considered and the relevant data should be deleted. Using the classic cycle slip detection approach, most of the cycle slips can be detected except for some small ones. As multi-stations data are used and compared, thus the small cycle slip of one satellite will be discovered by the observation residuals and will not affect the results.

Estimation of the dynamic variation of the orbital maneuver with multi-station data
The 3D dynamic variation of the orbital maneuver can be estimated with the data from more than three stations tracking the same monitored satellite. The model can be expressed as follows (Tu et al., 2021): (2) sin(e) e < π/6 0.5 e ≥ π/6 where t r is the estimated epoch variation of the receiver clock error, which is obtained in the subsection 'Estimation of the epoch variation of the receiver clock error by carrier phase observations' , and m represents the monitored satellite. The error corrections and the stochastic model determination are the same as in subsection 'Estimation of the epoch variation of the receiver clock error by carrier phase observations' . The least squares method is used to estimate the 3D dynamic variation of the orbital maneuver at each epoch. According to Eq. (3), the broadcast ephemeris is of low precision, but the orbit error is very small as the epochdifferenced observation is used. As shown in the former studies, if the epoch-differenced observation is used for the velocity estimation, both the broadcast ephemeris and precise ephemeris can give nearly the same accuracy (Colosimo et al., 2011;Tu et al., 2016). Referring to the former study, the precision of the estimated satellite velocity is about 1-2 cm per 30 s during the no maneuver period (Tu et al., 2020), which is accurate enough for the broadcast orbit correction.

Calculation of the orbit correction for the orbit maneuver satellite
After estimating the 3D dynamic variation of the orbital maneuver, the orbit corrections can be calculated by the integration: where b represents the calculated 3D orbit correction, T k represents the time of the current epoch, T 0 is the time of the initial epoch, and int represents the data sampling interval, 30 s in this study.
Generally, the error character of the estimated satellite velocity is white noise (Tu et al., 2020). During the integration of the velocity, the errors of the velocity will not increase exponentially, thus the orbit corrections have no obvious drift error. But for the long-time integration, the atmosphere residuals may lead to a small drift error which can be corrected by a simple linear model (Tu et al., 2016).

Correction of the orbital maneuver error for users
During the orbit maneuver period, the satellite's position will have a large offset from its normal designed orbit. Usually, the broadcast ephemeris only marks the healthy information of the satellite, and does not provide the orbit correction information, thus the users cannot use the data of this satellite. In this study, the true position of the satellite can be recovered by adding the estimated orbit correction for the orbital maneuver satellite, which can be expressed as follows.
where, o represents the true orbit position, and r represents the raw orbit position which is calculated by the broadcast ephemeris. After applying the orbit correction, the orbital maneuver satellite can be normally used for the PVT services.

Datasets and process
Eight BDS sites were used in this study: CUT0, JFNG, PNGM, DARW, HOB2, MAJU, MRO1, and ULAB, which are the stations of the International GNSS (Global Navigation Satellite System) Service (IGS) network, their distribution are shown in Fig. 1. Three reference stations (marked as triangles) were used to estimate the orbital maneuver, and another five user stations (marked as circles) were used for the positioning, velocity, and timing validations. The receiver positions and descriptions for the sites are given in Table 1. All data are publicly available via the internet (ftp:// cddis. gsfc. nasa. gov/). The carrier phase and pseudorange data were analyzed with the software developed by our group, which can conduct Single Point Positioning (SPP), Precise Point Positioning (PPP), Velocity Estimation (VE), and Orbit Maneuver Error Solution (OMES). For SPP the broadcast ephemeris and the pseudorange data are used, and the estimated parameters are the receiver clock and the receiver position of each epoch. PPP uses the fixed orbits and clocks provided by the IGS (http:// igscb. jpl. nasa. gov) and estimates the receiver clock, the constant zenith troposphere delay, the receiver's static positions, and the carrier phase ambiguities. Using one week data, we estimated the average coordinates as the references for all the stations using PPP. For VE the broadcast ephemeris and the carrier phase data are used. It also uses the epochdifferenced observations of the two adjacent epochs, and the estimated parameters are the velocity and epoch variation of the receiver clock error. OMES also uses the broadcast ephemeris and the carrier phase data. First, we constrain the station's velocity (0.001 m per 30 s via the a priori standard deviation). We use the epoch-differenced observations for the velocity estimation (not for the maneuvering satellite) to determine the epoch variation of the receiver clock error at each reference station. We then use the epoch-differenced observations (only for the maneuvering satellite) from three or more reference stations to estimate the orbit maneuver variation, after the epoch variation of the receiver clock error is corrected. Then, the estimated maneuver variations are integrated into orbital error corrections. Users can then use the orbital error corrections to recover the normal positioning, velocity, and timing services. The ephemeris, position, velocity, and orbital error corrections are referenced to the WGS84 Earth-Centered Earth-Fixed (ECEF) coordinate system in the International Terrestrial Reference Frame (ITRF) 2000 (Altamimi et al., 2002). The BDS carrier phase observation requires ambiguity estimation. For PPP, the ambiguity is estimated as a constant in a continuous period but needs to be reinitialized when a cycle slip occurs. For VE and OMES, the epoch differenced observations are used, thus no ambiguity parameter is in the data solution. The stochastic model is determined by the measurement error and satellite elevations (Xu et al., 2013). The dual-frequency ionosphere-free observations are used to correct the ionosphere errors and the other errors such as earth tide, ocean tide, phase center offset and variation, earth rotation effects, relativistic effects, and phase wind-up can be corrected using existing models (Genrich & Bock, 1992).

Estimation of the corrections of orbit errors for the maneuver satellite
The BDS C01 satellite undertook an orbital maneuver on January 09, 2020, namely the Day of Year (DOY) 009 in 2020, from Coordinated Universal Time (UTC) 10:00 to 16:00. This was marked in the broadcast ephemeris. Thus, the corresponding datasets of these eight stations during this period were used for the data solution and validation. In general, the actual maneuver period was about 30 min (Qiao et al., 2019), but the abnormal of broadcast ephemeris may last several hours. The total time span when the satellite has abnormal ephemeris is called the orbital maneuver period in this study. Figure 2a shows a time series of the estimated epoch variation for the C01 satellite. In general, the estimated epoch variation refers to the positional difference between the true orbit and the broadcast orbit during one sample interval. Thus, if there is no orbital maneuver, the estimated relative variation of the satellite should be zero as in this case the broadcast orbit is the true Fig. 1 Locations of the three reference stations (marked as triangles) and five user stations (marked as circles). All stations can track the BDS C01 satellites. In total, eight to ten BDS satellites were observed during the data collection period orbit. Figure 2a shows that before the orbital maneuver occurred or the maneuver was finished, the epoch variation was close to zero. Figure 2b shows the details of the velocity variation during the actual maneuver period. It can be seen that the adjustment is in the radial components as the cross and along track components are nearly zero values. In addition, Table 2 gives the total velocity changes in the radial, along-track, and cross-track components of the five satellite maneuver events (i.e., C01, DOY 009 in 2020; C02, DOY 003 in 2020; C03, DOY 020 in 2020; C04, DOY 308 in 2019; C05, DOY 023 in 2020). It can be concluded that the GEO satellites are mainly maneuvered in the cross-track and radial components.
The velocity variation was integrated to produce the orbital error correction at the corresponding time for this satellite. Figure 3 shows the time series of the correction of the orbital error for BDS C01. For this maneuver event, the orbital error correction reached up to 3,000 m in the 3D components. The accuracy of the estimated orbit correction can be evaluated during the no maneuver period. In this experiment, the Root Mean Square (RMS) of the correction's error is 4.2, 5.2, 5.3 cm in the X, Y, and Z components. When the orbital maneuver started, the estimated correction began to deviate from the zero value in the 3D components. While the orbital maneuver finished, the estimated variation went back to zero value, as the true orbit again overlapped with the new broadcast orbit.

Recovery of the abnormal PVT service
We analyzed the single point positioning performance with the data from five stations (Fig. 4). The station coordinates obtained with the weekly precise point positioning solution were used as references whose precision is 1-2 cm (Li et al., 2015;Tu et al., 2013;Zumberge et al., 1997). The results show that the positioning biases in all three components exist and gradually increase to thousands of kilometers during the maneuvering period. These values then returned to the normal after the end of the maneuver. In comparison, the positioning bias seemed normal when the computed orbit error corrections were used. Furthermore, whether a user is inside or outside the reference network has no effect on the recovered position results. In general, height is more poorly determined by BDS than the horizontal position due to the strong correlation between atmospheric delay and elevation.  The BDS provides not only a positioning service, but also velocity and timing services. The satellite ephemeris is needed for computing velocity and timing, thus these variables are also influenced by satellite maneuvers. Our results also show that these abnormalities in velocity and timing services can be eliminated by using the orbital error corrections (Figs. 5 and 6). In general, the velocity is affected by the position variation between the two adjacent epochs, and the magnitude of the orbital error influence for this variable is relatively small (several cm/s to dm/s). The magnitude for the timing is on the scale of thousands of nanoseconds. After repairing the orbital maneuver errors, both velocity and timing estimations are recovered to the normal.
In addition, another four orbital maneuvers were also tested. They were satellites C02 (DOY 003 in 2020), C03 (DOY 020 in 2020), C04 (DOY 308 in 2019) and C05 (DOY 023 in 2020). Table 3 shows the average Standard Deviation (STD) of the positioning, velocity, and timing biases with different approaches used to process the maneuvering data. It can be concluded that the PVT results have large biases while the satellite with maneuvers. Deleting or repairing the maneuver satellite data can recover the normal position service, and the repairing has a little better performance than that of the deleting. In all, the average precision of position, velocity, and timing services are improved by 0.8 m, 0.1 mm/s and 0.16 ns, respectively, while the orbital maneuver recovery approach is used. In some cases where the number of observed satellites is limited, the repairing the maneuver satellite is a better practice.

Determination and monitoring of orbital maneuvering
Our results demonstrate the method presented here provides another important tool for the determination and monitoring of orbital maneuvering. Figure 7 shows an example of the orbit maneuver detection for satellite C01 (DOY 009 in2020) at the station of CUT0. Specifically, it shows that the Observed Minus Computed (OMC) values of the epoch-differenced observation are stable and have a small STD when no orbital maneuver occurs, but these values increase significantly during a maneuver. Thus, we used the OMC values and the empirical STD values of the OMC time series for the orbital maneuvers detecting. If the absolute value of the OMC is larger than three times of the empirical STD, and the variation lasts for at least five minutes, we define this as the start of the orbital maneuver. While an orbital maneuver has already started, the absolute value of the OMC smaller than three times of the empirical STD and also lasts for at least five minutes, we recorded it as the end of the orbital maneuver. Note that the defined end time is not the actual end time of the orbital maneuver, but the end time of the abnormal ephemeris. In addition, the empirical value of STD is calculated by the OMC time series during the non-orbital maneuvering period. As Fig. 7 shown, the solid line represents the detected time, which is swifter than the broadcast ephemeris records. Table 4 gives a comparison between the detected time and the recorded time by the broadcast ephemeris. For these test results, the recorded start time is about 24 min later than the detected start time, and the recorded end time is about 29 min later than the detected end time.

Discussions
Orbital maneuver occurs frequently and will influence not only the user's PVT performance, but also the satellite system's health. Thus, except for the detection of the orbital maneuvers, monitoring the dynamic variations of the satellites is also important. As this approach can obtain the variation of the orbital maneuver in real-time, it is possible to evaluate whether the orbital adjustment has reached the target from the dynamic monitoring results. In addition, it also provides a new way for the dynamic monitoring of a satellite's orbital status, which is useful for healthy operation of a navigation system.
Another possible application is for precise orbit determination. Usually, the precise orbit product cannot be provided for the whole day if the satellite undergoes a maneuver, because the thrust force cannot be precisely modeled. As the velocity variation of the orbit error is directly caused by the thrust force, the acceleration of the thrust force can be obtained by differentiating the velocity. Then, the obtained thrust force's acceleration can be used as prior information for the orbit integration, and precise orbit of the maneuvering satellite can be obtained. In this case, the maneuvering satellite can be used normally, and satellite usability will be increased for users' applications (Qiao & Chen, 2018). There are at least three techniques for improving BDS orbit maneuver error estimates that should be investigated in the future. Station specific errors can be removed by using BDS estimates for the same station from the previous days that were shifted in sidereal time (Dach et al., 2007). An alternate approach involves solving for the common errors among multiple stations, i.e., a spatial filter (Wdowinski et al., 1997). And the combination of multi-systems and multi-frequencies observations to improve the precision and accuracy.
In addition, the method of using BDS data to recover the normal positioning service has the advantage of low costs. The cost of a geodetic BDS receiver is approximately several thousand dollars, and all the continuously operating BDS receivers can be used (Di Federico, 2014). It is also easy to use this method in real-time applications, as the data can be obtained anytime from the broadcast ephemeris and three or more ground stations. The real-time data stream from the station is also easy to realize, which can then be employed for the long term and continuous dynamic operation. Finally, the resolution of high-frequency BDS and more global stations has not yet been fully examined.

Conclusions
BDS satellite orbital maneuver usually occurs, which will influence the normal PVT service. The conventional method to solve the orbital maneuver is to mark and delete this satellite, at the expense of decreasing the amount of observation data. In this study, we proposed a method to recover the abnormal PVT services caused by satellite orbital maneuvers.
Using the epoch-differenced velocity estimation model and the high-precision carrier phase observation, the correction of the satellite orbit can be calculated, which can be used to correct the orbit for the orbital maneuver satellite.  Table 3 The average STDs of the positioning, velocity, and timing biases for the different approaches with five orbital maneuvers ("Contain", "Repair", "Delete" means containing, repairing and deleting maneuver satellite, respectively) The real data from three BDS stations are used for calculating the orbit correction, and another five BDS stations are used for the PVT services test. The results show that the estimated orbit correction can be used to obtain the normal PVT services during the satellite orbital maneuvers period. Compared with the traditional method, our approach can improve the average precision of position, velocity, and timing services by 0.8 m, 0.1 mm/s and 0.16 ns, respectively. In addition, the approach can also be used for the orbital maneuver detection and monitoring.