A Robust Method to Detect BeiDou Navigation Satellite System Orbit Maneuvering/Anomalies and Its Applications to Precise Orbit Determination

The failure to detect anomalies and maneuvering of the orbits of navigation satellite sensors will deteriorate the performance of positioning and orbit determination. Motivated by the influence of the frequent maneuvering of BDS GEO and IGSO satellites, this paper analyzes the limitations of existing methods, where BDS orbit maneuvering and anomalies can be detected, and develops a method to solve this problem based on the RMS model of orbit mutual differences proposed in this paper. The performance of this method was assessed by comparison with the health flag of broadcast ephemeris, precise orbit products of GFZ, the O-C values of a GNSS station and a conventional method. The results show that the performance of the method developed in this paper is better than that of the conventional method when the periodicity and trend items are obvious. Meanwhile, three additional verification results show that the method developed in this paper can find error information in the merged broadcast ephemeris provided by iGMAS. Furthermore, from the testing results, it can be seen that the detection of anomaly and maneuvering items do not affect each other based on the robust thresholds constructed in this paper. In addition, the precise orbit of the maneuvering satellites can be determined under the circumstances that the maneuver information detected in this paper is used, and the root mean square (RMS) of orbit overlap comparison for GEO-03/IGSO-03 in Radial, Along, Cross, 1D-RMS are 0.7614/0.4460 m, 1.8901/0.3687 m, 0.3392/0.2069 m, 2.0657/0.6145 m, respectively.


Introduction
On 31 March 2015, the first New-Generation navigation satellite of BeiDou navigation satellite system (BDS) successfully launched into orbit, which marked BDS taking a key step from the service of the Asia-Pacific area to the global network. By 2020, there will be 35 BeiDou satellites in outer space, and the system will provide worldwide users with high accuracy navigation, positioning and timing services [1][2][3][4][5]. As an important part of the system, the Geostationary Orbit (GEO) and Inclined Geosynchronous Orbit (IGSO) satellites need to be frequently maneuvered to maintain geosynchronous characteristics. Due to the impacts of various perturbations while running the orbit, satellites may be extremely disturbed, which makes the abnormal condition of satellite position occur [6][7][8]. To adjust the strategies of positioning and orbit determination in a correct and timely manner, the abnormal condition and maneuvers must be found as soon as possible after they occur. Otherwise, it will cause serious impacts on the service performance of the positioning and orbit determination [9][10][11][12][13].
predicted from the broadcast ephemeris at epoch t i , and the RMS values of orbit mutual difference can be calculated from Equation (1): where, PRN is the ID of a satellite, and n is the number of epochs contained in the broadcast ephemeris which is needed to be analyzed. Considering the symmetry of RMS PRN , for simplicity, it can be discussed later, then we define the variable rms PRN (t) at the epoch t according to Equation (2): where, making rms PRN as the processing array, t re f is the selected reference epoch, t = t 1 , t 2 , t 3 , . . . , t n , and this paper is based on the variable rms PRN (t) to start the discussion and research.
The difference operator ∇ [20] is defined as follows: In order to build a model that conforms to the characteristics of the variable rms PRN (t), a large number of experiments were analyzed, but only some of them are displayed here(the data is 27 June 2016-11 July 2016, each analysis period includes a three-day broadcast ephemeris, that is, n = 72), and the others are very similar. Figures 1 and 2 shows the results for rms 2 (t)/rms 6 (t) and ∇rms 2 (t)/∇rms 6 (t) at 27 June 2016-29 June 2016; Figures 3 and 4 shows the results for rms 4 (t)/rms 5 (t) and ∇rms 4 (t)/∇rms 5 (t) at 30 June 2016-2 July 2016; Figures 5 and 6 shows the results for rms 7 (t)/rms 10 (t) and ∇rms 7 (t)/∇rms 10 (t) at 3 July 2016-5 July 2016; Figure 7 shows the results for rms 9 (t) and ∇rms 9 (t) at 9 July 2016-11 July 2016 (Section 3 will test PRN 1/3/8, so they are not displayed here).
The upper part of each figure represents rms PRN (t), and the lower part of each figure represents ∇rms PRN (t) (that is the First-order difference of rms PRN (t)).
where, making PRN rms as the processing array, ref t is the selected reference epoch, 1 2 3 , , ,..., n t t t t t  , and this paper is based on the variable ( ) PRN rms t to start the discussion and research.
The difference operator  [20] is defined as follows: In order to build a model that conforms to the characteristics of the variable                                It can be seen from the above figures that the change in rms PRN (t) over time shows a periodicity and trend characteristics obviously. Figures 1 and 5 showed that the orbital anomaly only affects rms PRN (t) at the current epoch, but the orbital maneuver has an influence on rms PRN (t) at the current epoch and the subsequent epoch, respectively. After the first-order difference, the effects of periodicity and trend can be removed, and when there is no orbital anomaly and maneuver, ∇rms PRN (t) will be displayed as the ergodic Gaussian random process.

Determination on the RMS Model of Orbit Mutual Difference
Based on above analysis, the RMS value of the BDS orbit mutual difference consists of trend items (TR PRN ), periodicity items (P PRN ), orbit anomaly items (δrms PRN ), orbit maneuver items (Crms PRN ) and gauss white noise (ε). Then, we believe that the RMS value of orbit mutual difference meets the model at any epoch t ( t = t 1 , t 2 , t 3 , . . . , t n ): where, E(ε) = 0 (E stands for mathematical expectation), but also δrms PRN , Crms PRN and ε are independent to each other, that is: From Section 2.1.1 and Equation (3), we can see that the effects of periodicity and trend can be removed after the first-order difference, that is: ≈ ∇δrms PRN (t) + ∇Crms PRN (t) + ∇ε(t) E(∇rms PRN (t)) = E(∇δrms PRN (t)) + E(∇Crms PRN (t)) + E(∇ε(t)) According to the invariance property of linear transformations of Gaussian distributions, if maneuvering and anomalies do not occur in the satellite orbit, ∇rms PRN (t) is a Gaussian random variable with zero expectation value as well: E(∇rms PRN (t)) = E(∇δrms PRN (t)) + E(∇Crms PRN (t)) + E(∇ε(t)) = E(∇ε(t)) = 0 In contrast, ∇rms PRN (t) will destruct Gaussian distributions, and based on this, a new method to detect maneuvering and anomalies is proposed.

A Robust Method of Detection of Orbit Anomalies and Maneuvering for BDS Satellite
Supposing the orbital arc (size of the sliding ephemeris window) of a satellite contains data of n epochs, and the epoch t m (t m > 1) is detected, then the steps for the detection of orbit anomalies and maneuvering for a BDS satellite are as follows: (1) Select n as the size of the sliding ephemeris window.
(2) Use the corresponding data of the broadcast ephemeris to calculate the orbital coordinate X ij , Y ij , Z ij at the epoch t j which are predicted from the broadcast ephemeris at epoch t i (where i, j = 1, 2, 3, . . . , n). When satellite maneuvering occurs at the epoch t m , the rms PRN (t) value is calculated as follows: where: after the first-order difference operation, ∇rms PRN (t) can obtain: When satellite anomalies occur at the epoch t m , the rms PRN (t) value is calculated as follows: After a one order difference operation, ∇rms PRN (t) can obtain: where: 3 When the satellite is within normal status at the epoch t m , the rms PRN (t) value is calculated as follows: After a one order difference operation, ∇rms PRN (t) can obtain: Because the median is provided with high robust and high breakdown contamination rates [21][22][23][24][25], this paper uses this characteristic to suppress the influence of anomalies when maneuvering is detected. Considering the characteristics of the processing array, the absolute value of the array ∇rms PRN (t) is first obtained by step (5), and then, according to the relationship between median and mean square error, the array of absolute value was used to calculate robust variance factors: (7) According to past experience, the paper chooses four times (it can be appropriate to enlarge if needed) the variance factors as thresholds of detection: (8) For the detection of maneuvering and anomalies at epoch t m , when satellite maneuvering or anomalies occur, ∇rms PRN (t) will destruct Gaussian distributions; according to this principle, the detection criteria were set as follows: 1 When satellite maneuvering occurs at the epoch t m , the following criteria can be obtained by Equations (8)- (11): 2 When satellite anomalies occur at the epoch t m , the following criteria can be obtained by Equations (12)-(15): 3 When the satellite stays within normal status at the epoch t m , the following criteria can be obtained by Equations (16)- (18): The anomalies and maneuvering of a BDS satellite orbit can be detected by the above steps (1-9). The following sections will use the BDS broadcast ephemeris to validate and analyze this proposed method.

Validation and Analysis
To test the ability of the proposed method, which is used to detect maneuvering and anomalies, this paper uses the data of the merged broadcast ephemeris provided by iGMAS and utilizes the following programs to detect maneuvering and anomalies of GEO and IGSO. Considering the fact that it is not currently available to get real information on BDS orbit anomalies and maneuvering. So the test results of the proposed method (method one) are mainly compared with the health flag of broadcast ephemeris, precise orbit products of GFZ, and the O-C values of station BJF1/WUH1, which is obtained by that the station-satellite distance minus the observation distance, the station coordinate is provided by the observation file, and the satellite coordinate is fitted using BDS broadcast ephemeris. It is possible to verify whether or not the RMS model of orbit mutual difference is reliable. Moreover, it can also be shown whether the method is effective or not. Furthermore, in order to analyze the performance of the RMS model and the detection method proposed in this paper, the proposed method will also be compared with the conventional method (method two) using rms PRN (t) obtained by step (4) to directly detect maneuvering, in which the threshold of detection is the empirical value, which is equal to 5000, as proposed by [19]. In addition, in order to analyze the influence of orbital maneuver on precise orbit determination of BeiDou satellites, the validity of the maneuver information detected in this section will be analyzed in the last part of this section from the perspective of precise orbit determination.
Since the solutions of precise orbit determination often use single-day, three-day and seven-day broadcast ephemeris products as the initial orbit, the size of the windows selected for use in the test were single-day, three-day and seven-day.
Program one: The testing for the detecting ability of BDS orbit maneuvering of the robust method proposed in this paper. For GEO-03, the time is from 5 January 2015 00:00 to 7 January 2015 23:00; for IGSO-03, the time is from 9 January 2015 05:00 to 10 January 2015 05:00. The data comes from the BDS merged broadcast ephemeris provided by iGMAS.

Validation and Analysis
To test the ability of the proposed method, which is used to detect maneuvering and anomalies, this paper uses the data of the merged broadcast ephemeris provided by iGMAS and utilizes the following programs to detect maneuvering and anomalies of GEO and IGSO. Considering the fact that it is not currently available to get real information on BDS orbit anomalies and maneuvering. So the test results of the proposed method (method one) are mainly compared with the health flag of broadcast ephemeris, precise orbit products of GFZ, and the O-C values of station BJF1/WUH1, which is obtained by that the station-satellite distance minus the observation distance, the station coordinate is provided by the observation file, and the satellite coordinate is fitted using BDS broadcast ephemeris. It is possible to verify whether or not the RMS model of orbit mutual difference is reliable. Moreover, it can also be shown whether the method is effective or not. Furthermore, in order to analyze the performance of the RMS model and the detection method proposed in this paper, the proposed method will also be compared with the conventional method (method two) using (t) PRN rms obtained by step (4) to directly detect maneuvering, in which the threshold of detection is the empirical value, which is equal to 5000, as proposed by [19]. In addition, in order to analyze the influence of orbital maneuver on precise orbit determination of BeiDou satellites, the validity of the maneuver information detected in this section will be analyzed in the last part of this section from the perspective of precise orbit determination.
Since the solutions of precise orbit determination often use single-day, three-day and seven-day broadcast ephemeris products as the initial orbit, the size of the windows selected for use in the test were single-day, three-day and seven-day.
Program one: The testing for the detecting ability of BDS orbit maneuvering of the robust method proposed in this paper. For GEO-03, the time is from 5 January 2015 00:00 to 7 January 2015 23:00; for IGSO-03, the time is from 9 January 2015 05:00 to 10 January 2015 05:00. The data comes from the BDS merged broadcast ephemeris provided by iGMAS. (1) Testing and analysis in GEO-03(PRN3). The testing results of methods one and two are shown in Figure 8 (X-axis represents time, unit is hours; Y-axis represents value calculated by Equation (3), unit is meters); from 5 January 2015 9:00 to 5 January 2015 14:00, the health flag of the broadcast ephemeris showed that GEO-03 is in an unhealthy state; on 5 January 2015, among the products of precise orbits, GFZ did not calculate the coordinates of GEO-03; at station WUH1, the O-C values of GEO-03 are shown in Figure 9 (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance of L3, unit is meters); From Figure 8 and the other three verification methods for GEO-03, the testing results from 5 January 2015, at 14:00-15:00, are shown in Table 1.  (1) Testing and analysis in GEO-03(PRN3). The testing results of methods one and two are shown in Figure 8 (X-axis represents time, unit is hours; Y-axis represents value calculated by Equation (3), unit is meters); from 5 January 2015 9:00 to 5 January 2015 14:00, the health flag of the broadcast ephemeris showed that GEO-03 is in an unhealthy state; on 5 January 2015, among the products of precise orbits, GFZ did not calculate the coordinates of GEO-03; at station WUH1, the O-C values of GEO-03 are shown in Figure 9 (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance of L3, unit is meters); (1) Testing and analysis in GEO-03(PRN3). The testing results of methods one and two are shown in Figure 8 (X-axis represents time, unit is hours; Y-axis represents value calculated by Equation (3), unit is meters); from 5 January 2015 9:00 to 5 January 2015 14:00, the health flag of the broadcast ephemeris showed that GEO-03 is in an unhealthy state; on 5 January 2015, among the products of precise orbits, GFZ did not calculate the coordinates of GEO-03; at station WUH1, the O-C values of GEO-03 are shown in Figure 9 (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance of L3, unit is meters); From Figure 8 and the other three verification methods for GEO-03, the testing results from 5 January 2015, at 14:00-15:00, are shown in Table 1.  From Figure 8 and the other three verification methods for GEO-03, the testing results from 5 January 2015, at 14:00-15:00, are shown in Table 1. The results in Table 1 show that on 5 January 2015, at 14:00-15:00, both methods detect that GEO-03 maneuvering may occur: the health flag of broadcast ephemeris shows that GEO-03 is in a unhealthy state, but in the previous five hours, the health flag shows that GEO-03 is also in an unhealthy state; GFZ might think that GEO-03 is in an unhealthy state; and, there is a jump in the O-C value at WUH1, which may represent that GEO-03 is in an unhealthy state. Furthermore, it can be concluded from Figure 8 that the detection result of the conventional method is affected by the selected reference epoch, and the result of the epoch with far from the reference epoch is obviously affected by the periodicity and trend term, so that the detection result is not completely correct. But the new method is not affected by these factors: (2) Testing and analysis in IGSO-03(PRN8). The testing results of methods one and two are shown in Figure 10. From 9 January 2015 14:00 to 9 January 2015 20:00, the health flag of broadcast ephemeris showed that IGSO-03 is in an unhealthy state, while at 11:00, the health flag is in a healthy state. On 9 January 2015, among the products of precise orbits, GFZ did not calculate the coordinate of IGSO-03. At station BJF1, the O-C values of IGSO-03 are shown in Figure 11 (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance of B1, unit is meters); The results in Table 1 show that on 5 January 2015, at 14:00-15:00, both methods detect that GEO-03 maneuvering may occur: the health flag of broadcast ephemeris shows that GEO-03 is in a unhealthy state, but in the previous five hours, the health flag shows that GEO-03 is also in an unhealthy state; GFZ might think that GEO-03 is in an unhealthy state; and, there is a jump in the O-C value at WUH1, which may represent that GEO-03 is in an unhealthy state. Furthermore, it can be concluded from Figure 8 that the detection result of the conventional method is affected by the selected reference epoch, and the result of the epoch with far from the reference epoch is obviously affected by the periodicity and trend term, so that the detection result is not completely correct. But the new method is not affected by these factors: (2) Testing and analysis in IGSO-03(PRN8). The testing results of methods one and two are shown in Figure 10. From 9 January 2015 14:00 to 9 January 2015 20:00, the health flag of broadcast ephemeris showed that IGSO-03 is in an unhealthy state, while at 11:00, the health flag is in a healthy state. On 9 January 2015, among the products of precise orbits, GFZ did not calculate the coordinate of IGSO-03. At station BJF1, the O-C values of IGSO-03 are shown in Figure 11 (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance of B1, unit is meters); From Figure 10 and the other three verification methods for IGSO-03, from 9 January 2015, at 10:00-11:00 and 19:00-20:00, the testing results are shown in Tables 2 and 3, respectively. The results in Table 1 show that on 5 January 2015, at 14:00-15:00, both methods detect that GEO-03 maneuvering may occur: the health flag of broadcast ephemeris shows that GEO-03 is in a unhealthy state, but in the previous five hours, the health flag shows that GEO-03 is also in an unhealthy state; GFZ might think that GEO-03 is in an unhealthy state; and, there is a jump in the O-C value at WUH1, which may represent that GEO-03 is in an unhealthy state. Furthermore, it can be concluded from Figure 8 that the detection result of the conventional method is affected by the selected reference epoch, and the result of the epoch with far from the reference epoch is obviously affected by the periodicity and trend term, so that the detection result is not completely correct. But the new method is not affected by these factors: (2) Testing and analysis in IGSO-03(PRN8). The testing results of methods one and two are shown in Figure 10. From 9 January 2015 14:00 to 9 January 2015 20:00, the health flag of broadcast ephemeris showed that IGSO-03 is in an unhealthy state, while at 11:00, the health flag is in a healthy state. On 9 January 2015, among the products of precise orbits, GFZ did not calculate the coordinate of IGSO-03. At station BJF1, the O-C values of IGSO-03 are shown in Figure 11 (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance of B1, unit is meters); From Figure 10 and the other three verification methods for IGSO-03, from 9 January 2015, at 10:00-11:00 and 19:00-20:00, the testing results are shown in Tables 2 and 3, respectively. From Figure 10 and the other three verification methods for IGSO-03, from 9 January 2015, at 10:00-11:00 and 19:00-20:00, the testing results are shown in Tables 2 and 3, respectively.  Figure 10, it can be seen that anomalies do not affect the detection of maneuvering. From the results in Table 3, it can be seen that on 9th January 2015, at 19:00-20:00, both methods detect that IGSO-03 maneuvering may occur: the health flag of broadcast ephemeris shows that IGSO-03 is in an unhealthy state, while in the previous five hours and the later one hour, the health flag shows that IGSO-03 is also in an unhealthy state; GFZ may think that IGSO-03 is in an unhealthy state; and, there is a jump in the O-C value at BJF1, which may represent that IGSO-03 is in an unhealthy state.
Program two: Testing for the detecting ability of BDS orbit anomalies by the robust method proposed in this paper. For GEO-01, the time was from 4th January 2015 00:00 to 10th January 2015 23:00. The data comes from the BDS merged broadcast ephemeris provided by iGMAS.
Testing and analysis in GEO-01(PRN1). The testing results of method one and two are shown in Figure 12. On 8th January 2015 19:00, the health flag of broadcast ephemeris showed that GEO-01 is in a healthy state. Among the products of precise orbit, GFZ calculated the coordinates of GEO-01. At station BJF1, the O-C values of GEO-01 are shown in Figure 13 (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance of B1, unit is meters).    Figure 10, it can be seen that anomalies do not affect the detection of maneuvering. From the results in Table 3, it can be seen that on 9th January 2015, at 19:00-20:00, both methods detect that IGSO-03 maneuvering may occur: the health flag of broadcast ephemeris shows that IGSO-03 is in an unhealthy state, while in the previous five hours and the later one hour, the health flag shows that IGSO-03 is also in an unhealthy state; GFZ may think that IGSO-03 is in an unhealthy state; and, there is a jump in the O-C value at BJF1, which may represent that IGSO-03 is in an unhealthy state.
Program two: Testing for the detecting ability of BDS orbit anomalies by the robust method proposed in this paper. For GEO-01, the time was from 4th January 2015 00:00 to 10th January 2015 23:00. The data comes from the BDS merged broadcast ephemeris provided by iGMAS.
Testing and analysis in GEO-01(PRN1). The testing results of method one and two are shown in Figure 12. On 8th January 2015 19:00, the health flag of broadcast ephemeris showed that GEO-01 is in a healthy state. Among the products of precise orbit, GFZ calculated the coordinates of GEO-01. At station BJF1, the O-C values of GEO-01 are shown in Figure 13 (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance of B1, unit is meters).  From Figure 12 and other three verification ways for GEO-01, on 8 January 2015, at 18:00-19:00, the testing results are shown in Table 4. From the results in Table 4, it can be seen that on 8 January 2015, at 18:00-19:00, both methods detected that GEO-01 anomalies may occur that are caused by decoding errors or other reasons: the health flag of broadcast ephemeris show that GEO-01 is in a healthy state; GFZ may think that GEO-01 is in a healthy state; and, there is a jump in the O-C value at BJF1, which may represent that GEO-01 is in an unhealthy state. Moreover, it can be seen from Figure 12 that the results at both ends of the detection window due to the impact of the periodicity and trend term is clearly wrong.
From Figure 10, it can be seen that maneuvering does not affect detection of anomalies. From the results in Table 2, it can be seen that on 9 January 2015, at 10:00-11:00, both methods detect that IGSO-03 may occur anomalies that are caused by decoding errors or other reasons: the health flag of broadcast ephemeris show that IGSO-03 is in a healthy state; GFZ may think that IGSO-03 is in an unhealthy state; and, there is a jump in the O-C value at BJF1, which may represent that IGSO-03 is in an unhealthy state.
Program three: In order to further test the effect of orbital maneuvering information on the precise orbit determination, the orbit maneuvering information is used for precise orbit determination. As neither GFZ nor WHU provides the precise orbits of the two satellites when these two satellites maneuvered their orbits, orbit overlap comparison is used to analyze the precision of the orbit determination.
(1) GEO-03 precise orbit maneuvering test. In this test, the BDS measured data of 20 stations from MGEX and 10 stations from iGMAS from 3 to 6 January 2015 were used to obtain the final precision orbit using three-day long arc POD. The comparison strategy is shown in Figure 14, and the statistical results of orbit overlap comparison in radial direction, along direction, cross direction, 1D-RMS are shown in Table 5;  Table 4. From the results in Table 4, it can be seen that on 8 January 2015, at 18:00-19:00, both methods detected that GEO-01 anomalies may occur that are caused by decoding errors or other reasons: the health flag of broadcast ephemeris show that GEO-01 is in a healthy state; GFZ may think that GEO-01 is in a healthy state; and, there is a jump in the O-C value at BJF1, which may represent that GEO-01 is in an unhealthy state. Moreover, it can be seen from Figure 12 that the results at both ends of the detection window due to the impact of the periodicity and trend term is clearly wrong.
From Figure 10, it can be seen that maneuvering does not affect detection of anomalies. From the results in Table 2, it can be seen that on 9 January 2015, at 10:00-11:00, both methods detect that IGSO-03 may occur anomalies that are caused by decoding errors or other reasons: the health flag of broadcast ephemeris show that IGSO-03 is in a healthy state; GFZ may think that IGSO-03 is in an unhealthy state; and, there is a jump in the O-C value at BJF1, which may represent that IGSO-03 is in an unhealthy state.
Program three: In order to further test the effect of orbital maneuvering information on the precise orbit determination, the orbit maneuvering information is used for precise orbit determination. As neither GFZ nor WHU provides the precise orbits of the two satellites when these two satellites maneuvered their orbits, orbit overlap comparison is used to analyze the precision of the orbit determination.
(1) GEO-03 precise orbit maneuvering test. In this test, the BDS measured data of 20 stations from MGEX and 10 stations from iGMAS from 3 to 6 January 2015 were used to obtain the final precision orbit using three-day long arc POD. The comparison strategy is shown in Figure 14, and the statistical results of orbit overlap comparison in radial direction, along direction, cross direction, 1D-RMS are shown in Table 5;  (2) IGSO-03 precise orbit maneuvering test. In this test, the BDS measured data from 7 January to 10 2015 were used to obtain the final precision orbit. The comparison strategy is shown in Figure 15, and the statistical results of orbit mutual difference of overlapping arc in radial direction, along direction, cross direction, 1D-RMS are shown in Table 6;  From the results of GEO-03 and IGSO-03 precise orbit maneuvering tests, we can see that it is difficult to provide precise orbit determination results for GEO-03 and IGSO-03 when there is no satellite orbit maneuvering information. However, when using the maneuver information detected in this paper, the precise orbit determination results of these two satellites not only can be  (2) IGSO-03 precise orbit maneuvering test. In this test, the BDS measured data from 7 to 10 January 2015 were used to obtain the final precision orbit. The comparison strategy is shown in Figure 15, and the statistical results of orbit mutual difference of overlapping arc in radial direction, along direction, cross direction, 1D-RMS are shown in Table 6;  (2) IGSO-03 precise orbit maneuvering test. In this test, the BDS measured data from 7 January to 10 2015 were used to obtain the final precision orbit. The comparison strategy is shown in Figure 15, and the statistical results of orbit mutual difference of overlapping arc in radial direction, along direction, cross direction, 1D-RMS are shown in Table 6;  From the results of GEO-03 and IGSO-03 precise orbit maneuvering tests, we can see that it is difficult to provide precise orbit determination results for GEO-03 and IGSO-03 when there is no satellite orbit maneuvering information. However, when using the maneuver information detected in this paper, the precise orbit determination results of these two satellites not only can be  From the results of GEO-03 and IGSO-03 precise orbit maneuvering tests, we can see that it is difficult to provide precise orbit determination results for GEO-03 and IGSO-03 when there is no satellite orbit maneuvering information. However, when using the maneuver information detected in this paper, the precise orbit determination results of these two satellites not only can be determined, but also the accuracy of the orbit determination is equivalent to that of the same type of orbit determination (The accuracy of GEO/IGSO are meters and decimeters, respectively). What's more, the GEO-03 is more accurate than the GEO-02 in terms of orbit determination accuracy, and the differences in 1D-RMS accuracy between IGSO-03 and IGSO-02, IGSO-04 1D-RMS are in the order of decimeters.
The above programs can verify that the ability of the robust method proposed in this paper, which can detect BDS orbit maneuvering and anomalies, is effective and that they can also verify the reliability of the proposed RMS model. In the process of satellite orbit determination, if the maneuvering moments detected in this paper is segmented, this satellite can perform normal precise orbit determination. The results of precision orbit determination when the maneuver information detected in this paper was used are shown that the 1D-RMS is 2.0657 m before the maneuver of GEO-03, and the 1D-RMS is 0.6145 m before the maneuver of IGSO-03. It should be said that due to the limited resolution of the maneuver detection, the maneuver time points aren't determined precisely, so, only the orbit arcs before the maneuver time points are tested.

Conclusions and Outlook
In this study, an RMS model of orbit mutual differences was established according to characteristics of orbital prediction accuracy, and a robust method to detect BDS orbit anomalies and maneuvering was proposed based on the RMS model.
BDS orbit anomalies and maneuvering were detected using the proposed method in this study, and the testing results were then verified by a health flag of broadcast ephemeris, precise orbit products of GFZ, the O-C values of station BJF1/WUH1 and the conventional method. By comparison with the results of the conventional method, it can be seen that when the periodicity and trend items are obvious, the performance of the method proposed in this paper is more stable than the conventional method. Moreover, the empirical threshold used by the conventional method does not always meet the actual requirements, while the robust thresholds constructed in this paper have a higher reliability. From the other three verification results, it can be seen that the RMS model of orbit mutual difference proposed in this paper is reliable; furthermore, the detecting method of maneuvering and anomalies developed in this paper is effective. In addition, this method can also detect inaccurate information in the merged broadcast ephemeris provided by iGMAS, such as the error health flag of broadcast ephemeris. The testing results also show that, based on the robust thresholds constructed in this paper, anomalies do not affect the detection of maneuvering; meanwhile, maneuvering does not affect the detection of anomalies, and this method has greater consistency compared with the other verification results. During this experiment, it was shown that the sliding window and the reference epoch can be adjusted according to the actual demand; this also reflects that the method proposed in this paper has greater flexibility.
In addition, the maneuver information detected in this paper will have a significant impact on the precise orbit determination of satellites, which can not only realize the precise orbit determination of maneuvered satellites, but also the accuracy of the orbit determination is equivalent to that of the same type of satellite orbit determination. The experimental results show that the precise orbit determination results of maneuvered satellites can't be obtained before using the satellite maneuver information. However, after using the maneuver information detected in this paper, the results of orbit overlap comparison of GEO-03 in radial direction, along direction, cross direction, 1D-RMS are 0.7614 m, 1.8901 m, 0.3392 m, 2.0657 m, respectively, and the results of orbit overlap comparison of IGSO-03 in radial direction, along direction, cross direction, 1D-RMS are 0.4460 m, 0.3687 m, 0.2069 m, 0.6145 m, respectively.
Although the proposed robust method to detect BDS orbit anomalies and maneuvering was successful, the results of this paper depend on BDS merged broadcast ephemeris, which is updated every hour by iGMAS. As such, the time resolution is limited, so a method to detect maneuvering and anomalies in higher time resolutions will be addressed in future work.