An Optimized Method to Detect BDS Satellites’ Orbit Maneuvering and Anomalies in Real-Time

The orbital maneuvers of Global Navigation Satellite System (GNSS) Constellations will decrease the performance and accuracy of positioning, navigation, and timing (PNT). Because satellites in the Chinese BeiDou Navigation Satellite System (BDS) are in Geostationary Orbit (GEO) and Inclined Geosynchronous Orbit (IGSO), maneuvers occur more frequently. Also, the precise start moment of the BDS satellites’ orbit maneuvering cannot be obtained by common users. This paper presented an improved real-time detecting method for BDS satellites’ orbit maneuvering and anomalies with higher timeliness and higher accuracy. The main contributions to this improvement are as follows: (1) instead of the previous two-steps method, a new one-step method with higher accuracy is proposed to determine the start moment and the pseudo random noise code (PRN) of the satellite orbit maneuvering in that time; (2) BDS Medium Earth Orbit (MEO) orbital maneuvers are firstly detected according to the proposed selection strategy for the stations; and (3) the classified non-maneuvering anomalies are detected by a new median robust method using the weak anomaly detection factor and the strong anomaly detection factor. The data from the Multi-GNSS Experiment (MGEX) in 2017 was used for experimental analysis. The experimental results and analysis showed that the start moment of orbital maneuvers and the period of non-maneuver anomalies can be determined more accurately in real-time. When orbital maneuvers and anomalies occur, the proposed method improved the data utilization for 91 and 95 min in 2017.


Introduction
BeiDou Navigation Satellite System (BDS) has been providing continuous Positioning, Navigation, and Timing (PNT) services across the whole Asia-Pacific region since 27 December 2012. It aims to serve global users upon its completion in 2020 [1][2][3][4][5][6][7][8]. The BDS is influenced by Earth's non-spherical gravity and other perceptual factors, which lead to long-term perturbations of the offset of the satellite location and orbital elements. In order to keep the satellite in the normal range of the design orbit, orbit maneuvering is necessary to adjust the orbit of the satellite using the propulsion systems. In maintaining geosynchronous characteristics, the frequency of orbit maneuvering for Geostationary Orbit (GEO) and Inclined Geosynchronous Orbit (IGSO) satellites is higher than in Medium Earth Orbit (MEO) satellites. The position of satellite will vary by tens of kilometers because of orbit maneuvering, causing serious impact on the orbit determination and the service performance of PNT [9][10][11]. In addition, the abnormal condition of the satellite position would occur because of the impact of various perturbations while running the orbit. The earlier the maneuver is detected, the sooner the strategies of positioning and orbit determination can be adjusted [12][13][14][15][16][17][18]. The information on satellite maneuvers and anomalies is not available publicly, and are broadcasted on a time-interval basis

The Detection Theory of Orbital Maneuvers/Anomaly Using a One-Step Method
When the satellite is maneuvered or has an anomaly, the real-time satellite position calculated by the broadcast ephemeris is not correct given the increasing errors. The strategy of the BeiDou system is to mark the satellite health status as unhealthy about one hour before orbit maneuvering. This means that the real-time positions of satellite calculated by the broadcast ephemeris are not correct before about one hour, which reduces the utilization rate of available satellite data. On the other hand, the health marks for BDS satellites from the broadcast ephemeris are misidentified or sometimes missing, which leads to the information being received by common users to be unreliable. Considering that the real-time pseudo-range observations are not influenced by maneuvers or anomalies, the absolute residuals between pseudo-range observations and the pseudo-range calculated by the broadcast ephemeris will increase during the period of orbital maneuver or anomaly. Thus, a robust method using the pseudo-range residuals is proposed to detect the orbital maneuver and anomaly. This is an effective way that could be useful to common users on every static station.
The pseudo-range residual will contain increasing errors because of the satellite position being changed by orbit maneuvering. When the satellite orbit has an anomaly, the pseudo-range residual will change suddenly in the corresponding time period. These characteristics of different situations can be used for detection. Two factors-namely the orbital maneuver detection factor (L M ) and the anomaly detection factor (S M )-were defined and used to detect the start time of orbit maneuvering and the period of non-maneuvering anomaly.

Orbital Maneuver Detection Factor
The orbital maneuver detection factor is calculated by subtracting the absolute residual value of the pseudo-range from the empirical threshold of the residual value of the pseudo-range. The modified pseudo-range (D j i (k)) is corrected for the receiver clock offsets, the ionospheric delay, and the tropospheric delay from the raw observations ( D j i (k)).D j i (k) is the credible variable to reflect the distance between the satellite and the station during the orbital maneuver period. The distance between the station and satellite is calculated using the broadcast ephemeris, which contains the gross error because of the unreliable orbital parameters provided by the broadcast ephemeris. L j i (k) is the absolute residual value of observations which could be used to detect the orbit maneuvering and is computed by where | | is the function of the absolute value.D  (2) where the L j Max is the empirical threshold of the pseudo-range observations belonging to satellite j. The empirical threshold L j Max of satellite j is given in advance. When L j M (k) is bigger than 0 and keeps a continuous upward trend during a period (i.e., 10 min chosen in this paper), the corresponding time of L j M (k) is considered to be the start moment of the orbit maneuver. For the empirical threshold of the observations for satellite j, L j Max is the key to determine the start moment of orbit maneuvering. Considering that the residual of the observations obeys normal distribution, these residuals in the normal condition of the BDS satellite are all in the confidence interval of the 99.74% (3σ) confidence coefficient. In other words, if the factor is out of the corresponding confidence interval, it would be considered abnormal. Thus, the L j Max for the BDS satellite is provided by the upper confidence interval of the residuals of the observations in the 99.74% confidence coefficient.
Because there is no uniform correction model of the receiver clock offset, the least-squares method can be used to estimate the receiver clock offset. It needs to be emphasized that the least-squares method would be evidently affected by the blunders. Once the orbital parameters in the broadcast ephemeris are unusable and have big biases, the receiver clock offset cannot be calculated correctly. Then, the pseudo-range observations corrected by this receiver clock offset will have big residuals compared with the calculated pseudo-range, even if the corresponding satellites are healthy and usable. This leads to a difficulty to confirm the PRN of the maneuvering satellite. Therefore, this paper proposed an optimized method of a robust weight matrix for the observations, in order to avoid disturbing the orbit deviation due to the pseudo-range residuals of normal satellites. This robust least-squares method is used to estimate the receiver clock offset as follows.
Equation (3) is the observation equation for estimating the receiver clock offset, where B is the coefficient matrix of the estimated parameters δt i , l is a constant term, V is the corrections for observations, and P is the robust equivalence weight matrix. Considering the rapid changes of L j i (k) after the orbital maneuver or anomaly, the two-stage function of the robust equivalent weight matrix is structured as follows [27,28].
where the p j (k) is the diagonal element of the observation weight matrix for the j satellite at the k epoch, which is in accordance to the weight rule of the satellite elevation angle. p j (k) is the diagonal element of robust equivalence weight matrix P and r is the critical index, which is assigned as 1.0 in this work.
If one satellite j is suspected to be the maneuvering orbit or non-maneuvering anomaly at the epoch k, p j (k) would be assigned as 0 (i.e., 0 for doubtful maneuvering or anomaly, 1 for normal). Then the receiver clock offset at the next epoch k + 1 would be estimated with p j (k) valued at 0. Thus, the L j i (k + 1) of other normal satellites assigned as 1 would not be influenced. When L j i (k) is bigger than 0 and shows a continuous growth trend within a time period, the first epoch with L j i (k) > 0 is considered to be the start time of the maneuver.

Orbital Anomaly Detection Factor
When the satellite orbit has a non-maneuvering anomaly, the pseudo-range residual will change suddenly in the corresponding time period. If this change is too great, this leads to the increase of L j i (k) for many satellites, and the p j (k) of these satellites would be marked as 0. When L j i (k) is more than 0, the satellite is considered a suspicious object for anomaly. The number of suspect objects is defined as Q in the same epoch. According to the magnitude of Q, the non-maneuvering anomaly is divided into two categories: weak anomaly and strong anomaly. When the value of Q is less than the redundant observation number, it is considered a weak anomaly. Otherwise, it is considered a strong anomaly. The redundant observation number for estimating the receiver clock offset is defined as R. R is equal to N − 1, and N is the number of satellites in one epoch. The critical index r in Equation (5) is assigned as 3.0 in this work.

The Weak Anomaly Detection Factor
The method to calculate the weak anomaly detection factor is the same as the method for the orbit maneuvering. The weak anomaly detection factor can be calculated by where the S j M − (k) is the weak anomaly detection factor for j satellite at k epoch. When S j M − (k) is bigger than 0 and shows abnormal jumps in a time period, the corresponding period of S j M − (k) > 0 is considered to be the period of non-maneuvering weak anomaly.

Strong Anomaly Detection Factor
When Q is greater than or equal to the redundant observation number (R), it is considered a strong anomaly. This will cause the receiver clock offset in the next epoch to have a low reliability or to not be solved. Thus, we used the median robust method to solve this problem as follows.
where V is the correction for observations of the satellite i. l is a constant term. N is the total number of observations. D j i (k) is corrected using δt i in Equation (8), then the new L j i (k) and p j (k) are updated using Equations (1)- (5). Then, the δt i will be calculated with the new p j (k) using Equations (3) and (4). The following steps are the same with the steps in Section 2.1. The strong anomaly detection factor can be calculated by S j M where S j M + (k) is the factor of non-maneuvering strong anomaly for the j satellite at the k epoch.
When Specifically, when Q marked by the new p j (k) is equal to N, it is considered the anomaly for the user terminal.

Steps for Orbital Maneuver and Anomaly Detection
The following steps are used for orbital maneuver detection in this work.  The following steps are used for orbital anomaly detection in this work.
(2) If the S j M (k) is greater than 0, p j (k) and Q are calculated by Equation (5).
(3) Compare the Q with the R. If Q is smaller than R, go to step (4); otherwise, go to step (5).   When the Q marked by the new p j (k) is equal to N, it is considered as the anomaly for user terminal.

Data Description
The data of the Multi-GNSS Experiment (MGEX) station SIN1 was selected to analyze the detection results of BeiDou GEO and IGSO satellites' orbital maneuvers and anomalies. The station SIN1 is located in Singapore (1 • 20 N, 103 • 40 E, see Figure 1) with an LEIAR25.R3 GNSS receiver. The sample interval is 30 s. The coordinates of SIN1 can be achieved from the IGS Solution Independent Exchange (SINEX) product or using the Precise Point Positioning technique. In this work, the SINEX coordinates of SIN1 are chosen to verify the detection performance of GEO and IGSO. The distributions of trajectories on the station SIN1 are shown in Figure 1.  From Figure 1, the locations of the SIN1 station are marked as red circles. The satellites moving with trajectories in "8" shapes are IGSO satellites, and the satellites with trajectories in a point shape are GEO satellites. The satellites with distributions of trajectories with in the yellow circle can be observed by SIN1.
The thresholds of the residuals of the pseudo-range for different satellites of the SIN1 station on doy (day of year) 029 (29 January), doy 140 (20 May), doy 229 (17 August), and doy 296 (23 October) of 2017 are shown in Table 1. The blue marked PRNs are the GEO satellites, and the unmarked PRNs are the IGSO satellites. Consider that the thresholds of the residuals of the pseudo-range keep steady only for a short period of time, which is not over one week [27]. Table 1 shows that the values of are different in different days. Therefore, should be updated frequently, and we suggest updating it every three days.
In addition, the selected stations for maneuver detection of BDS MEO satellites are show in Table 2. The distributions of trajectories on the ground and the location of the selected stations are shown in Figure 2. From Figure 1, the locations of the SIN1 station are marked as red circles. The satellites moving with trajectories in "8" shapes are IGSO satellites, and the satellites with trajectories in a point shape are GEO satellites. The satellites with distributions of trajectories with in the yellow circle can be observed by SIN1.
The thresholds of the residuals of the pseudo-range for different satellites of the SIN1 station on doy (day of year) 029 (29 January), doy 140 (20 May), doy 229 (17 August), and doy 296 (23 October) of 2017 are shown in Table 1. The blue marked PRNs are the GEO satellites, and the unmarked PRNs are the IGSO satellites. Consider that the thresholds of the residuals of the pseudo-range keep steady only for a short period of time, which is not over one week [27]. Table 1 shows that the values of L Max are different in different days. Therefore, L Max should be updated frequently, and we suggest updating it every three days.
In addition, the selected stations for maneuver detection of BDS MEO satellites are show in Table 2. The distributions of trajectories on the ground and the location of the selected stations are shown in Figure 2.   In Figure 2, the locations of the selected stations are marked as red circles. The red trajectory is for a MEO satellite on the ground during a regression period. The area of green circle is the operating range of MEO satellites monitored by the corresponding station, which is analyzed by the Systems Tool Kit (STK). The area above the purple line is the range monitored by the UCAL station (The station UCAL is located in Canada with the TRIMBLE NETR9 GNSS receiver), and the area under the white line is the range monitored by GAMB station (The station GAMB is located in French Polynesia with the TRIMBLE NETR9 GNSS receiver). The seven selected stations using the proposed method ensure that the operating range of MEO satellites can be monitored for orbital maneuver and anomaly detection all the time.

The Orbital Maneuver Detection for BDS
Using the one-step method, the orbital maneuver was detected for C03 on 23 October 2017. The detection results of orbit maneuvering for SIN1 are firstly given, with the orbital maneuver detection factor series shown in Figure 3.
In Figure 3, the detection factor series shows a trend continually on the rise of about 4.5 h. The value of the orbital maneuver detection factor is greater than 0 from 9:38:30. The start time of the orbit maneuvering of C03, determined by the orbital maneuver detection factor on 23 October 2017, is 9:38:30.
We used the health marks for satellites (i.e., 0 for healthy status and 1 for unhealthy status) from the broadcast ephemeris provided by MGEX and the precise orbit products from international GNSS Monitoring & Assessment System (iGMAS) to verify the correction of the orbit maneuvering detection results.
The health marks for satellites in the broadcast ephemeris of C03 and the header information of the precise orbit products are shown in  In Figure 2, the locations of the selected stations are marked as red circles. The red trajectory is for a MEO satellite on the ground during a regression period. The area of green circle is the operating range of MEO satellites monitored by the corresponding station, which is analyzed by the Systems Tool Kit (STK). The area above the purple line is the range monitored by the UCAL station (The station UCAL is located in Canada with the TRIMBLE NETR9 GNSS receiver), and the area under the white line is the range monitored by GAMB station (The station GAMB is located in French Polynesia with the TRIMBLE NETR9 GNSS receiver). The seven selected stations using the proposed method ensure that the operating range of MEO satellites can be monitored for orbital maneuver and anomaly detection all the time.

The Orbital Maneuver Detection for BDS
Using the one-step method, the orbital maneuver was detected for C03 on 23 October 2017. The detection results of orbit maneuvering for SIN1 are firstly given, with the orbital maneuver detection factor series shown in Figure 3.
In Figure 3, the detection factor series shows a trend continually on the rise of about 4.5 h. The value of the orbital maneuver detection factor is greater than 0 from 9:38:30. The start time of the orbit maneuvering of C03, determined by the orbital maneuver detection factor on 23 October 2017, is 9:38:30. We used the health marks for satellites (i.e., 0 for healthy status and 1 for unhealthy status) from the broadcast ephemeris provided by MGEX and the precise orbit products from international GNSS Monitoring & Assessment System (iGMAS) to verify the correction of the orbit maneuvering detection results.
The health marks for satellites in the broadcast ephemeris of C03 and the header information of the precise orbit products are shown in Figures 4 and 5.         In Figure 4, the health marks for the satellites from the broadcast ephemeris could be marked unhealthy ahead of orbit maneuvering. Specifically, C03 is marked unhealthy from 8:00:00, while it is detected as orbit maneuvering at 9:38:30. Satellite orbit cannot be precisely determined by the conventional model of orbital mechanics due to the intervention of the orbiting maneuvering force. In Figure 4, the health marks for the satellites from the broadcast ephemeris could be marked unhealthy ahead of orbit maneuvering. Specifically, C03 is marked unhealthy from 8:00:00, while it is detected as orbit maneuvering at 9:38:30. Satellite orbit cannot be precisely determined by the conventional model of orbital mechanics due to the intervention of the orbiting maneuvering force. The C03 satellite misses based on the header of the precise orbit products for the corresponding date, indicating the correction of orbit maneuvering detection results.
Also, the SPP results of another site named XMIS (located in Christmas Island, AU 10 • 26 S, 105 • 41 E) were used to further validate the correction of the detection results. The SPP deviations of XMIS on doy 296 of 2017 are shown in Figure 6. In Figure 6, the detected start time of orbit maneuvering and the start time of the unhealthy marks are pointed out. The deviations marked blue and green are in same level, while the deviations marked red are obviously higher and show a rapid upward trend. Thus, the satellite starts to maneuver at 1-1.5 h after it is marked unhealthy in the broadcast ephemeris. This strategy results in the loss of observations for 1-1.5 h and reduces the utilization of valid observations. On the other hand, after the detected start time of orbit maneuvering, the SPP deviation gradually increases to a tendency where the positioning result is not reliable. The above positioning results can further prove the correction of the detected start time of orbit maneuvering. In order to justify that the start time of orbital maneuver detected by this method is more accurate than previous research (Huang et al., 2017), the start time of orbit maneuvering for C03 on 23 October 2017 (which calculated by the previous method) is 10:01:30. There is a delay of about 23 min from the start time detected by the two-steps method.
The applicability of the optimized method for detecting the orbit maneuvering of IGSO and MEO was verified by the detected results, which are not shown in this paper because of the reasonable length limitation.

The Orbital Maneuver Detection for GPS
In order to verify the applicability of the detection method for orbital maneuvers in other GNSS In Figure 6, the detected start time of orbit maneuvering and the start time of the unhealthy marks are pointed out. The deviations marked blue and green are in same level, while the deviations marked red are obviously higher and show a rapid upward trend. Thus, the satellite starts to maneuver at 1-1.5 h after it is marked unhealthy in the broadcast ephemeris. This strategy results in the loss of observations for 1-1.5 h and reduces the utilization of valid observations. On the other hand, after the detected start time of orbit maneuvering, the SPP deviation gradually increases to a tendency where the positioning result is not reliable. The above positioning results can further prove the correction of the detected start time of orbit maneuvering. In order to justify that the start time of orbital maneuver detected by this method is more accurate than previous research (Huang et al., 2017), the start time of orbit maneuvering for C03 on 23 October 2017 (which calculated by the previous method) is 10:01:30. There is a delay of about 23 min from the start time detected by the two-steps method.
The applicability of the optimized method for detecting the orbit maneuvering of IGSO and MEO was verified by the detected results, which are not shown in this paper because of the reasonable length limitation.

The Orbital Maneuver Detection for GPS
In order to verify the applicability of the detection method for orbital maneuvers in other GNSS MEO constellations using the one-step method, the orbital maneuver is detected for G03 satellite of GPS on 10 January 2017. The orbit maneuvering has been detected by the BRUN station. The results of BRUN are shown, with the orbital maneuver detection factor series shown in Figure 7.
In Figure 7, the detection factor series shows a continually rising trend for hours. The value of the orbital maneuver detection factor is greater than 0 from 17:22:00. The start time of the orbit maneuvering of G03 is determined by the orbital maneuver detection factor on 10 January 2017, and is 17:22:00. We used the health status marks for satellites (i.e., 0 for healthy status and 1 for unhealthy status) from the broadcast ephemeris provided by MGEX and the precise orbit products from iGMAS to verify the correction of the orbit maneuvering detection results. This is not shown because of the reasonable length limitation. The health status marks of satellite G03 in the broadcast ephemeris is marked as unhealthy from 15:59:44 to 21:57:52. In addition, the satellites' orbit cannot be precisely determined using the conventional model of orbital mechanics due to the intervention of the orbiting maneuvering force. The G03 satellite misses based on the header of the precise orbit products from iGMAS for the corresponding date, indicating the correction of orbit maneuvering detection results. Therefore, the new one-step method also can be used to determine the start time of orbital maneuvers of other GNSS MEO constellations.

Non-Maneuvering Anomaly Detection for a Weak Anomaly
The non-maneuvering anomaly was detected as a weak anomaly on 29 January 2017. The performance of the non-maneuvering weak anomaly detection for SIN1 is firstly given, with an anomaly detection factor series shown in Figure 8. We used the health status marks for satellites (i.e., 0 for healthy status and 1 for unhealthy status) from the broadcast ephemeris provided by MGEX and the precise orbit products from iGMAS to verify the correction of the orbit maneuvering detection results. This is not shown because of the reasonable length limitation. The health status marks of satellite G03 in the broadcast ephemeris is marked as unhealthy from 15:59:44 to 21:57:52. In addition, the satellites' orbit cannot be precisely determined using the conventional model of orbital mechanics due to the intervention of the orbiting maneuvering force. The G03 satellite misses based on the header of the precise orbit products from iGMAS for the corresponding date, indicating the correction of orbit maneuvering detection results. Therefore, the new one-step method also can be used to determine the start time of orbital maneuvers of other GNSS MEO constellations.

Non-Maneuvering Anomaly Detection for a Weak Anomaly
The non-maneuvering anomaly was detected as a weak anomaly on 29 January 2017. The performance of the non-maneuvering weak anomaly detection for SIN1 is firstly given, with an anomaly detection factor series shown in Figure 8.

Non-Maneuvering Anomaly Detection for a Weak Anomaly
The non-maneuvering anomaly was detected as a weak anomaly on 29 January 2017. The performance of the non-maneuvering weak anomaly detection for SIN1 is firstly given, with an anomaly detection factor series shown in Figure 8.  In Figure 8, the detection factor series shows that the period of non-maneuvering anomaly for C04 detected by the weak anomaly detection factor is 8:20:00 to 8:46:00 on 29 January 2017.
In order to justify the correction of the detection results for SIN1, the precise orbit products provided by iGMAS and the pseudo-range SPP bias of the other station, XMIS, are used as references.
We used the information of precise orbit products from iGMAS to verify the correction of the orbital anomaly detection results. The pseudo-range SPP bias of the other station, XMIS, was used as a reference.
The information of the header file from the precise orbit products provided by iGMAS contained the C04 satellite in Figure 9. It is proved that there was no orbital maneuver for the C04 satellite. In Figure 8, the detection factor series shows that the period of non-maneuvering anomaly for C04 detected by the weak anomaly detection factor is 8:20:00 to 8:46:00 on 29 January 2017.
In order to justify the correction of the detection results for SIN1, the precise orbit products provided by iGMAS and the pseudo-range SPP bias of the other station, XMIS, are used as references.
We used the information of precise orbit products from iGMAS to verify the correction of the orbital anomaly detection results. The pseudo-range SPP bias of the other station, XMIS, was used as a reference.
The information of the header file from the precise orbit products provided by iGMAS contained the C04 satellite in Figure 9. It is proved that there was no orbital maneuver for the C04 satellite.   In Figure 8, the detection factor series shows that the period of non-maneuvering anomaly for C04 detected by the weak anomaly detection factor is 8:20:00 to 8:46:00 on 29 January 2017.
In order to justify the correction of the detection results for SIN1, the precise orbit products provided by iGMAS and the pseudo-range SPP bias of the other station, XMIS, are used as references.
We used the information of precise orbit products from iGMAS to verify the correction of the orbital anomaly detection results. The pseudo-range SPP bias of the other station, XMIS, was used as a reference.
The information of the header file from the precise orbit products provided by iGMAS contained the C04 satellite in Figure 9. It is proved that there was no orbital maneuver for the C04 satellite.  The red square marks the unhealthy status and the corresponding time. C04 is unhealthy and is marked as 1 between 8:00:00 and 9:12:00, which is not shown because of the reasonable length limitation. Figure 10. The C04 satellite is marked unhealthy in the broadcast ephemeris on 29 January 2017. The red square marks the unhealthy status and the corresponding time. C04 is unhealthy and is marked as 1 between 8:00:00 and 9:12:00, which is not shown because of the reasonable length limitation.  In Figure 11, the detected periods of non-maneuvering anomaly and the periods of unhealthy marks are pointed out. The deviations marked blue and green are in same level, while the deviations marked red are obviously higher and show a jumping trend. The deviations of SPP results between the marked period and the detected period are usable and not influenced by the non-maneuvering anomaly. Figure 11 can verify that the usable data is lost, with the unhealthy status marked in the broadcast ephemeris; the detection method can improve the utilization rate of the observed data. In addition, in Figure 11, the deviations after the detected period jumped sharply due to the non-maneuvering anomaly, rendering the SPP results unreliable. Also, this period (from 8:20:00 to 8:46:00) is definitely close to the real period of the non-maneuvering anomaly.

Non-Maneuvering Anomaly Detection for a Strong Anomaly
The non-maneuvering anomaly was detected as a strong anomaly on 20 May 2017. The performance of the non-maneuvering strong anomaly detection for SIN1 is given, with the non-maneuvering strong anomaly detection factor series shown in Figure 12. In Figure 11, the detected periods of non-maneuvering anomaly and the periods of unhealthy marks are pointed out. The deviations marked blue and green are in same level, while the deviations marked red are obviously higher and show a jumping trend. The deviations of SPP results between the marked period and the detected period are usable and not influenced by the non-maneuvering anomaly. Figure 11 can verify that the usable data is lost, with the unhealthy status marked in the broadcast ephemeris; the detection method can improve the utilization rate of the observed data. In addition, in Figure 11, the deviations after the detected period jumped sharply due to the nonmaneuvering anomaly, rendering the SPP results unreliable. Also, this period (from 8:20:00 to 8:46:00) is definitely close to the real period of the non-maneuvering anomaly.

Non-Maneuvering Anomaly Detection for a Strong Anomaly
The non-maneuvering anomaly was detected as a strong anomaly on 20 May 2017. The performance of the non-maneuvering strong anomaly detection for SIN1 is given, with the nonmaneuvering strong anomaly detection factor series shown in Figure 12. In Figure 12, the detection factor series shows that the period of non-maneuvering anomaly for C02 detected by the strong anomaly detection factor is 2:00:00 to 3:00:00 on 20 May 2017.
The information of precise orbit products from iGMAS is as shown in Figure 13. The SPP deviations of the other station, XMIS, are used as references. In Figure 12, the detection factor series shows that the period of non-maneuvering anomaly for C02 detected by the strong anomaly detection factor is 2:00:00 to 3:00:00 on 20 May 2017.
The information of precise orbit products from iGMAS is as shown in Figure 13. The SPP deviations of the other station, XMIS, are used as references. In Figure 12, the detection factor series shows that the period of non-maneuvering anomaly for C02 detected by the strong anomaly detection factor is 2:00:00 to 3:00:00 on 20 May 2017.
The information of precise orbit products from iGMAS is as shown in Figure 13. The SPP deviations of the other station, XMIS, are used as references. Figure 13. The information of header file from the precise orbit products on 20 May 2017. The red squares mark the corresponding date and highlight the precise orbit of C02 has been determined.
The information of the header file from the precise orbit products provided by iGMAS contained the C02 satellite in Figure 13. It is proven that there is no orbit maneuvering for the C02 satellite. The information of the header file from the precise orbit products provided by iGMAS contained the C02 satellite in Figure 13. It is proven that there is no orbit maneuvering for the C02 satellite.
Also, the health marks of C02 and the SPP results of another site, XMIS, is used to further validate the correction of the detection results. The SPP deviations of XMIS on doy 140 of 2017 is shown in Figures 14 and 15.  In Figure 15, the detected periods of non-maneuvering anomalies are pointed out. The SPP deviations after the detected period jumped sharply because of the non-maneuvering anomaly, leading the SPP results to be unreliable. Also, this period (from 2:00:00 to 3:00:00) is definitely close to the real period of the non-maneuvering anomaly. From Figure 14, the healthy status of C02 is marked as 0 between 2:00:00 and 3:00:00 in the broadcast ephemeris. Therefore, the real-time robust detection method not only improves the utilization of observations for common users with effective data, but also makes an effective supplement for the health status of the BeiDou satellites in the broadcast ephemeris.  In Figure 15, the detected periods of non-maneuvering anomalies are pointed out. The SPP deviations after the detected period jumped sharply because of the non-maneuvering anomaly, leading the SPP results to be unreliable. Also, this period (from 2:00:00 to 3:00:00) is definitely close to the real period of the non-maneuvering anomaly. From Figure 14, the healthy status of C02 is marked as 0 between 2:00:00 and 3:00:00 in the broadcast ephemeris. Therefore, the real-time robust detection method not only improves the utilization of observations for common users with effective data, but also makes an effective supplement for the health status of the BeiDou satellites in the In Figure 15, the detected periods of non-maneuvering anomalies are pointed out. The SPP deviations after the detected period jumped sharply because of the non-maneuvering anomaly, leading the SPP results to be unreliable. Also, this period (from 2:00:00 to 3:00:00) is definitely close to the real period of the non-maneuvering anomaly. From Figure 14, the healthy status of C02 is marked as 0 between 2:00:00 and 3:00:00 in the broadcast ephemeris. Therefore, the real-time robust detection method not only improves the utilization of observations for common users with effective data, but also makes an effective supplement for the health status of the BeiDou satellites in the broadcast ephemeris.

Orbit Maneuvering and Anomaly Detection for the BeiDou Satellites in 2017
We used the results from 001 to 300 doy of 2017 to verify the applicability of the optimized method for long-term orbit maneuvering and anomaly detection. The detected results are compared with the health marks from the broadcast ephemeris. The detected orbit maneuverings are all marked unhealthy in the period. The results are shown in Table 3. For the non-maneuvering anomaly, some health marks in the broadcast ephemeris are missed for BDS satellites, as shown in Table 4. Table 3. Results of orbit maneuvering detection in 2017 for BDS (measure in hour: minute: second). The detected satellites, the times of orbital maneuvers, and the average time differences between the detected start moment and the marked moment are shown.  Table 3 shows that there were 46 orbital maneuvers of BDS satellites in 2017 (001 to 300 doy). IGSO satellites had 13 orbital maneuvers. MEO satellites had only one orbital maneuver. This is because the GEO satellites need more maneuvers to maintain the static character with Earth because of weak observation geometry and the lack of orbit dynamics [29]. The frequency of orbit maneuvering for GEO satellites is higher than that for IGSO and MEO satellites. The average time difference between the detected start moment using the method and the marked started moment from the broadcast ephemeris for BeiDou satellites is 91 min, which is the time period the data is still available in. Therefore, the optimized method could detect orbit maneuvering more accurately in real-time and improve the utilization of observations. Table 4. The detection results for non-maneuvering anomalies in 2017 for BeiDou satellites (measured in hour: minute: second). The detected satellites, the times of anomalies, and the average duration differences between the marked period and the detected period are shown.  Table 4 shows that the BeiDou GEO satellites had 102 non-maneuvering anomalies in year 2017. IGSO satellites had 13 non-maneuvering anomalies. MEO satellites had 12 non-maneuvering anomalies. The frequency of non-maneuvering anomalies for GEO satellites is higher than that for IGSO and MEO satellites. The average duration in difference between the marked period and the detected period is 95 min, which is the time period when the observational data is available. There were 17 unmarked health marks for BDS in the broadcast ephemeris. Therefore, the proposed detection method for an orbital anomaly not only improves the utilization of observations, but also supplements the information of health status for satellites in the broadcast ephemeris.
Orbit maneuvering and non-maneuvering anomalies appeared in the same day (C01 and C04 on 9 January 2017). This can be detected by the orbital maneuver detection factor and the anomaly detection factor in real-time. Thus, the proposed method performs well even when the orbital maneuvers and non-maneuvering anomalies appear on the same day.

Discussion and Conclusions
In this paper, an optimized robust method is proposed to detect orbit maneuvering and non-maneuvering anomalies for BeiDou GEO, IGSO, and MEO satellites using only the pseudo-range observations, the broadcast ephemeris, and the known station coordinates. The data of MGEX are analyzed, and the results show that the satellite orbital maneuvers and the non-maneuvering anomalies can both be detected accurately. In addition, the average time of orbit maneuvering differences between the marked started time in the broadcast ephemeris and the detected start time was 91 min in 2017. It indicates that the presented detection method for orbit maneuvering can enlarge the available observational data by about 1.52 h when the satellite maneuvers. The average duration of the non-maneuvering anomaly difference between the marked period of the broadcast ephemeris and the detected period was 94.95 min in 2017, which indicates that the presented detection method for orbital anomaly could enlarge the available observational data by about 1.58 h when the satellite has an anomaly. The proposed method of orbit maneuvering detection and anomaly detection could be used well together in real-time.
Also, the presented method is more effective and could be implemented more easily in any BeiDou static station. The proposed method could improve the PNT service during orbit maneuvers and anomalies; in addition, it also can be programmed as a data preprocessing tool to enhance the quality and reliability of GNSS real-time services in the applications.