sensors Article

A Robust Method to Detect BeiDou Navigation Satellite System Orbit Maneuvering/Anomalies and Its Applications to Precise Orbit Determination Fei Ye 1,2, *, Yunbin Yuan 1 , Bingfeng Tan 1 and Jikun Ou 1 1

2

*

State Key Laboratory of Geodesy and Earth’s Dynamics, Institute of Geodesy and Geophysics, 340 Xudong Rd., Wuhan 430077, China; [email protected] (Y.Y.); [email protected] (B.T.); [email protected] (J.O.) University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing 100049, China Correspondence: [email protected]

Received: 29 March 2017; Accepted: 12 May 2017; Published: 16 May 2017

Abstract: 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. Keywords: BDS; detection of anomalies and maneuvering; robust threshold; broadcast ephemeris; precise orbit determination (POD)

1. 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–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–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–13].

Sensors 2017, 17, 1129; doi:10.3390/s17051129

www.mdpi.com/journal/sensors

Sensors 2017, 17, 1129

2 of 17

The study of monitoring the orbits of space objects has long been conducted [14–16]; however, few studies focus on detecting the maneuvering and anomalies of navigational satellites, including GEO and IGSO satellites. At present, it is mainly based on using tracking data of orbit determination to detect the maneuvering and anomalies of GEO [17]. In general, previous research can be divided into four categories. First, radar data is used to analyze the tracking of space targets. It uses spaceborne laser/ground radar to measure the state of motion of space targets and then track these spacecrafts. An advantage is the high accuracy of measurement, but currently, the data from spaceborne laser/ground radar that can be used for navigation satellites is limited. Therefore, it is not conducive to the daily detection of the anomalies and maneuvers of navigation satellites. Second, it uses the wavelet analysis method to identify anomalies [18]; during such studies, a specific frequency signal is extracted by wavelet decomposition, after which the orbital anomalies can be distinguished using the detection theory of signal singularity. However, it is difficult to determine the wavelet basis and decomposition scale adaptively; thus, it is not easy to adaptively distinguish orbital anomalies. Third, the tracking data of orbit determination is used to monitor the maneuvering and anomalies of GEO [17]. It is based on highly accurate satellite two-way time comparison that can effectively detect orbital maneuvering, but it is too expensive to get tracking data of orbit determination for the BDS navigation satellite. It has proven difficult to popularize the application, so a method is needed that can cut costs to detect maneuvering and anomalies. Fourth, the satellite orbit calculated by the ephemeris is used to detect maneuvering and anomalies according to the jump in satellite position [19]. Raw data is easily accessible, but such research is still in its infancy, and few documents propose the orbit mutual difference RMS model of BDS broadcast ephemeris. The thresholds that can detect BDS maneuvering and anomalies do not have the robust characteristics as well. Therefore, the current body of research has its limitations in the detection of BDS anomalies and maneuvers. This paper analyzes in depth the fourth type of research described above and attempts to utilize the BDS merged broadcast ephemeris provided by iGMAS to build an RMS model of orbit mutual difference. Furthermore, this paper proposes a robust method to detect BDS orbit anomalies and maneuvering that is sufficiently convenient to popularize the application. To verify the reliability of the established model and the validity of the proposed method according to the characteristics of BDS maneuvering and anomalies, solutions will be used to test the ability of the method in detecting maneuvering and anomalies. In addition, a third set of program will be used to analyze the impact of the maneuvering information detected in this paper on the precise orbit determination of BeiDou satellites. 2. A Robust Method and Its Implementation Scheme of Detection of BDS Orbit Anomalies and Maneuvering BDS broadcast ephemeris can reflect information about satellite orbit at a fast update frequency (1 h). In this paper, BDS broadcast ephemeris was selected as the original experimental data and a method was proposed to detect BDS maneuvering and anomalies. It can be summarized as follows: use satellite orbit coordinates predicted from different epochs to the same reference epoch to calculate the RMS value of orbit mutual difference; get values of epoch-differences from the selected processing array; select robust thresholds and construct criteria according to the assumed characteristics of the RMS model of orbit mutual difference; and finally, detect orbit maneuvering and anomalies. 2.1. Introduction to the RMS Model of Orbit Mutual Difference 2.1.1. Analysis of the RMS Values of Orbit Mutual Difference Before introducing the RMS model of orbit mutual difference, it is necessary to analyze the RMS values of orbit mutual difference in detail, and then summarize the model according to their
characteristics. Assuming that the orbital coordinates are Xij , Yij , Zij at the epoch t j which are

Sensors 2017, 17, 1129

3 of 17

Sensors 2017, 17, 1129

3 of 17

predicted from the broadcast ephemeris at epoch t i , and the RMS values of orbit mutual difference can be calculated from Equation (1): predicted from the broadcast ephemeris at epoch ti , and the RMS values of orbit mutual difference can 2 2 2 k n k n k n be calculated from Equation (1): X ik  X jk    k 1 Yik  Yjk    k 1  Zik  Z jk  )   k 1 (1) RMSPRN  ti v ,t j   3n u  2 2  2 u k=n k=n k=n
t ∑k=1 Xik −nX jk + ∑k=1 Yik − Yjk + ∑k=1 Zik − Zjk where, PRN is the ID of a satellite, and is the number of epochs contained in the broadcast RMSPRN ti , t j = (1) 3n ephemeris which is needed to be analyzed. , for simplicity, it cancontained be discussed then ephemeris we define Considering PRNthe where, PRN is the the ID ofsymmetry a satellite,ofandRMS n is number of epochs in thelater, broadcast

the variable at the epoch t according to Equation (2): rmsPRN t ) analyzed. which is needed to (be Considering the symmetry of RMSPRN , for simplicity, it can be discussed later, then we define   R M S PRN ( t ref , t1 ) RM S PRN ( t ref , t 2 ) ... R M S PRN ( t ref , t n )   rm s the variable rms PRN (PRN t) at the epoch t according to Equation (2):  (2) rm s PRN ( t )  RM S PRN ( t ref , t )  i ) h rms PRN = RMSPRN (tre f , t1 ) RMSPRN (tre f , t2 ) . . . RMSPRN (tre f , tn ) (2) where, making rmsPRN as the processing array, tref is the selected reference epoch, rms PRN (t) = RMSPRN (tre f , t) t  t1, t2 , t3 ,..., tn , and this paper is based on the variable rmsPRN (t ) to start the discussion and where, making rms PRN as the processing array, tre f is the selected reference epoch, t = t1 , t2 , t3 , . . . , tn , research. 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: The difference operator ∇ [20] is defined as follows: rms PRN ( tm )  rms PRN ( tm )  rms PRN ( tm 1 )  ) (3) k( t m ) = rms PRN ( t m k ) 1 − rms PRN ( t  ) ∇rms PRN m − 1  rms PRN ( tm )   (  rms PRN ( tm ))  (3) k k − 1 ∇ rms PRN (tm ) = ∇(∇ rms PRN (tm )) In order to build a model that conforms to the characteristics of the variable rmsPRN (t ) , a large In order to build a model conforms thesome characteristics of the variablehere(the rms PRN (data t), a large number of experiments were that analyzed, but to only of them are displayed is 27 number of experiments were analyzed, but only some of them are displayed here(the data is June 2016–11 July 2016, each analysis period includes a three-day broadcast ephemeris, that 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 rms2 (t ) / rms6 (t ) and nrms = 72), and the others are very similar. Figures 1 and 2 shows the results for rms2 (t)/rms6 (t) 2 (t ) / rms6 (t ) at 27 June 2016–29 June 2016; Figures 3 and 4 shows the results for rms4 (t ) / rms5 (t ) and ∇rms2 (t)/∇rms6 (t) at 27 June 2016–29 June 2016; Figures 3 and 4 shows the results for and rms4 (t ) /  rms5 (t ) at 30 June 2016–2 July 2016; Figures 5 and 6 shows the results for rms7 (t ) / rms4 (t)/rms 5 ( t ) and ∇rms4 ( t )/∇rms5 ( t ) at 30 June 2016–2 July 2016; Figures 5 and 6 shows the rms10 (t) and rms7 (t ) / rms10 (t ) at 3 July 2016–5 July 2016; Figure 7 shows the results for rms9 (t ) results for rms7 (t)/rms10 (t) and ∇rms7 (t)/∇rms10 (t) at 3 July 2016–5 July 2016; Figure 7 shows the and for 9 July 2016–11 July 2016 (Section 3 will test PRN 1/3/8, so they are not displayed rms9rms (t ) at results 9 ( t ) and ∇rms9 ( t ) at 9 July 2016–11 July 2016 (Section 3 will test PRN 1/3/8, so they are here). not displayed here). The upper upper part part of of each eachfigure figurerepresents representsrms represents rmsPRN t )t),,and The and the the lower lower part part of of each each figure represents PRN (( ∇ rms ( t ) (that is the First-order difference of rms ( t ) ). rmsPRN (t ) (that is the First-order difference of rmsPRN (t ) ). PRN PRN 3

x 10

6

rm s 2 /m

2

1

F irs t-o rd e r d iffe re n c e rm s2 /m

0 0

1

x 10

10

20

30

40

50

60

70

80

10

20

30

40 number/t

50

60

70

80

6

0

-1 -2 -3 0

Figure Theresults resultsfor forrmsrms (t ) and rms2 (2t()t.). Figure 1.1.The 2 ( t2) and and ∇rms

Sensors 2017, 17,17, 1129 Sensors 2017, 1129 Sensors 2017, 17, 1129

of 17 4 of417 4 of 17 4000 4000

s /m rm srm 6 / m6

3000 3000 2000 2000 1000 1000 0 0 0 0

10 10

20 20

30 30

40 40

50 50

60 60

70 70

80 80

10 10

20 20

30 30

40 number/t 40 number/t

50 50

60 60

70 70

80 80

t -oe rrddeiffe r d iffe c e srm/ms6 /m F ir sFt ir-osrd re n creenrm 6

2000 2000 1000 1000 0 0

-1000 -1000 -2000 0 -2000 0

Figure 2. The results for rms6 (t ) and rms6 (t ) . Figure for rms Figure2.2. The The results for rms66 ((tt)) and ∇ rms rms66((tt))..

s /m rm srm 4 /m4

5000 5000 4000 4000 3000 3000 2000 2000 1000 1000 0 0 0 0

10 10

20 20

30 30

40 40

50 50

60 60

70 70

80 80

10 10

20 20

30 30

40 number/t 40 number/t

50 50

60 60

70 70

80 80

irsrd t-oerrdderifferen d ifferen ce srm/ms4 /m F irsFt-o ce rm 4

2000 2000 1000 1000 0 0

-1000 -1000 -2000 0 -2000 0

Figure 3. The Theresults resultsfor for rms rms4 4(t()t)and and ∇ rms44((t )t).. Figure 3.  rms Figure 3. The results for rms4 (t ) and  rms4 ( t ) . 8000 8000

s /m rm srm 5 /m5

6000 6000 4000 4000 2000 2000 0 0 0 0

10 10

20 20

30 30

40 40

50 50

60 60

70 70

80 80

10 10

20 20

30 30

40 number/t 40 number/t

50 50

60 60

70 70

80 80

irsrd t-oe rrddeiffe r d iffe c e srm/ms5 /m F irsFt-o re n cree nrm 5

2000 2000 1000 1000 0 0

-1000 -1000 -2000 0 -2000 0

Figure 4. 4. The rms5 ((t )t) and rms Figure The results resultsfor for rms and ∇ rms5 ((t )t). . Figure 4. The results for rms55 (t ) and rms55( t ) .

Sensors 2017, 17, 1129 Sensors 2017, 17, 1129

5 of 17 5 of 17

Sensors 2017, 17, 1129 Sensors 2017, 17, 1129

5 of 17

7

5 of 17

x 10 4 x 107 4

rmrm s 7 s/m 7 /m

3 3 x 107 4 2 23

rm s 7 /m

1 12 0 00 1 0

10 10

30 30

40 40

20

30

40

20 20

30 30

20

30

40 40 number/t number/t 40 number/t

10

50 50

60 60

50

70 70

60

80 80

70

80

2 x 107 24

F irs t-o rd e r d iffe re n c e rm s7 /m

F irs d iffe e rm s7 /m F irs t-ot-o rd rd e r edriffe re nrecne crm s7 /m

7

x010 4 x 0107 4

20 20

0 2 0

-2 0 -2 -4 -2 -40 0

10 10

-4 0

10

50 50

60 60

50

70 70

60

80 80

70

80

Figure 5. The results for rms7 (t ) and rms7 (t ) . Figure 5. The results for rms7 (t ) and rms7 (t ) . Figure 5. The Theresults resultsfor for rms rms7 (7t()t)and and  ∇rms rms7 (7t()t.). Figure 5. 4000 4000

10 10 10

20 20 20

30 30 30

3000 3000 3000 2000 2000 2000 1000 1000 1000 0 00 -1000 -1000 -1000 -2000 0 -2000 -2000 00

10 1010

20 20 20

30 30 30

F irs t-o rd e r d iffe re n c e rm s1 0 /m

F irs d iffe e rm F irs t-ot-o rd rd e r edr iffe re nrecne crm s1 0s1/m 0 /m

10 rm s 1 0 /m

rmrm s s 1/m 0 /m

4000 3000 3000 3000 2000 2000 2000 1000 1000 1000 0 00 00 0

40 40

40

40

4040 number/t number/t number/t

50 50

50

60

50

60 60

70

60

5050

70 70

80

70

60 60

80 80

80

70 70

80 80

Figure for rms ∇rms rms rms10 ((tt)) and  rms10(10t()(t.)t).. Figure6. 6. The The results results for Figure 6. The results for rms and  rms101010(t()t) and rms 10 10 (t ) . 4000 4000 4000

rms9 /m

rms/m /m rms 9 9

3000 3000 3000 2000 2000 2000 1000 1000 1000 0 0 00 0 0

20 20 20

30 30

40 40

30

40 number/t 40

40

5050

6060

70 70

80 80

50

60

70

80

50

60

70

80

2000 2000 2000

First-order difference rms9 /m

First-order difference rms/m /m First-order difference rms 9 9

1010 10

1000 1000 1000

0 0

0

-1000

-1000 -1000

-2000

-2000 0 -20000 0

10

10 10

20

20 20

30 30

40 number/t number/t

50 50

60 60

Figure 7. The results for rms9 (t ) and rms9 (t ) . Figure rms ((t).. Figure7.7.The Theresults resultsfor for rms rms ((t t)) and ∇rms Figure 7. The results for rms999(t ) and rms999 (t ) .

70 70

80 80

Sensors 2017, 17, 1129

6 of 17

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. 2.1.2. 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 (PPRN ), 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 = t1 , t2 , t3 , . . . , tn ): rms PRN (t) = TR PRN (t) + PPRN (t) + δrms PRN (t) + Crms PRN (t) + ε(t)

(4)

where, E(ε) = 0 (E stands for mathematical expectation), but also δrms PRN , Crms PRN and ε are independent to each other, that is:  E(δrms PRN (t), Crms PRN (t)) = 0   E(δrms PRN (t), ε(t)) = 0   E(Crms PRN (t), ε(t)) = 0

(5)

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) = ∇ TR PRN (t) + ∇ PPRN (t) + ∇δrms PRN (t) + ∇Crms PRN (t) + ∇ε(t)   (6) ≈ ∇δ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

(7)

In contrast, ∇rms PRN (t) will destruct Gaussian distributions, and based on this, a new method to detect maneuvering and anomalies is proposed. 2.2. 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 tm (tm > 1) is detected, then the steps for the detection of orbit anomalies and maneuvering for a BDS satellite are as follows: (1) (2)

(3) (4) (5)

Select n as the size of the sliding ephemeris window. Use the corresponding data of the broadcast ephemeris to calculate the orbital coordinate
Xij , Yij , Zij at the epoch t j which are predicted from the broadcast ephemeris at epoch ti (where i, j = 1, 2, 3, . . . , n). Calculate the RMSPRN (i, j) value of orbit mutual difference according to Equation (1). Select the reference epoch tre f and get the values of the variable rmsPRN (t) according to Equation (2). The first-order difference operation is used for rmsPRN (t) according to Equation (3), and making ∇rmsPRN (t1 ) = 0:

Sensors 2017, 17, 1129

1

7 of 17

When satellite maneuvering occurs at the epoch tm , the rmsPRN (t) value is calculated as follows:   rmsPRN (ta ) = ε(ta )  (8) rmsPRN (tm ) = CrmsPRN (tm ) + ε(tm )   rmsPRN (tb ) = CrmsPRN (tm ) + ε(tb ) where: t1

t2

|

... {z

tm

tm−1 }

ta

tm+1 |

tm+2 {z

tn

...

(9)

}

tb

after the first-order difference operation, ∇rmsPRN (t) can obtain:   ∇rmsPRN (ta ) = ∇ε(ta )  ∇rmsPRN (tm ) = CrmsPRN (tm ) + ∇ε(tm )   ∇rmsPRN (tb ) = ∇ε(tb )

(10)

E(∇rmsPRN (ta )) = E(∇ε(ta )) = 0 E(∇rmsPRN (tm )) = E(CrmsPRN (tm )) + E(∇ε(tm )) = E(CrmsPRN (|e − r|)) E(∇rmsPRN (tb )) = E(∇ε(tb )) = 0 2

  

(11)

 

When satellite anomalies occur at the epoch tm , the rmsPRN (t) value is calculated as follows:   rmsPRN (ta ) = ε(ta )  (12) rmsPRN (tm ) = δrmsPRN (tm ) + ε(tm )   rmsPRN (tb ) = ε(tb ) After a one order difference operation, ∇rmsPRN (t) can obtain:

∇rmsPRN (ta ) = ∇ε(ta ) ∇rmsPRN (tm ) = δrmsPRN (tm ) + ∇ε(tm ) ∇rmsPRN (tm+1 ) = −δrmsPRN (tm ) + ∇ε(tm+1 ) ∇rmsPRN (tc ) = ∇ε(tc )

    

(13)

   

E(∇rmsPRN (ta )) = E(∇ε(ta )) = 0 E(∇rmsPRN (tm )) = E(δrmsPRN (tm )) + E(∇ε(tm )) = E(δrmsPRN (|e − r|)) E(∇rmsPRN (tm+1 )) = E(−δrmsPRN (tm )) + E(∇ε(tb )) = −E(δrmsPRN (|e − r|)) E(∇rmsPRN (tc )) = E(∇ε(tc )) = 0

    

(14)

   

where: t1 | 3

t2

... {z ta

tm

tm−1 }

tm+1

tm+2 |

tm+3 {z tc

...

tn }

(15)

When the satellite is within normal status at the epoch tm , the rmsPRN (t) value is calculated as follows:  rmsPRN (ta ) = ε(ta )   (16) rmsPRN (tm ) = ε(tm )  rmsPRN (tb ) = ε(tb )  After a one order difference operation, ∇rmsPRN (t) can obtain:

∇rmsPRN (ta ) = ∇ε(ta ) ∇rmsPRN (tm ) = ∇ε(tm ) ∇rmsPRN (tb ) = ∇ε(tb )

    

(17)

Sensors 2017, 17, 1129

8 of 17

  E(∇rmsPRN (ta )) = E(∇ε(ta )) = 0  E(∇rmsPRN (tm )) = E(∇ε(tm )) = 0  E(∇rmsPRN (tb )) = E(∇ε(tb )) = 0  (6)

(18)

Because the median is provided with high robust and high breakdown contamination rates [21–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 ∇rmsPRN (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: σ0 = med(abs(∇rms PRN (t)))/0.6745

(19)

(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: ) T1 = 4 × σ0 (20) T2 = −4 × σ0

(8)

For the detection of maneuvering and anomalies at epoch tm , 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 tm , the following criteria can be obtained by Equations (8)–(11):

∇rms PRN (tm ) > T1 ∇rms PRN (tm+1 ) ≥ T2 2

) (21)

) or

∇rms PRN (tm ) < T2 ∇rms PRN (tm+1 ) > T1

) (22)

When the satellite stays within normal status at the epoch tm , the following criteria can be obtained by Equations (16)–(18):

∇rms PRN (tm ) ≤ T1 ∇rms PRN (tm+1 ) ≥ T2 (9)

or

∇rms PRN (tm ) ≤ T1 ∇rms PRN (tm+1 ) < T2

When satellite anomalies occur at the epoch tm , the following criteria can be obtained by Equations (12)–(15):

∇rms PRN (tm ) > T1 ∇rms PRN (tm+1 ) < T2 3

)

) or

∇rms PRN (tm ) ≥ T2 ∇rms PRN (tm+1 ) ≤ T1

) (23)

The whole window moves backward at an epoch, repeating steps (2–8), until all epochs are finished.

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. 3. 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

3. Validation and Analysis

F i r s t - o r d e r d i f f e r e n c e r m s3 / m

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 Sensors 17,currently 1129 9 ofSo 17 that it 2017, is not available to get real information on BDS orbit anomalies and maneuvering. 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, 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 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 coordinate is provided by the observation file, and the satellite coordinate is fitted using BDS broadcast broadcast ephemeris. It is possible to verify whether or not the RMS model of orbit mutual ephemeris. It is possible to verify whether or not the RMS model of orbit mutual difference is reliable. difference is reliable. Moreover, it can also be shown whether the method is effective or not. Moreover, it can also be shown whether the method is effective or not. Furthermore, in order to analyze Furthermore, in order to analyze the performance of the RMS model and the detection method the performance of the RMS model and the detection method proposed in this paper, the proposed proposed in this paper, the proposed method will also be compared with the conventional method method will also be compared with the conventional method (method two) using rms PRN (t) obtained (method two) using rmsPRN (t) obtained by step (4) to directly detect maneuvering, in which the by step (4) to directly detect maneuvering, in which the threshold of detection is the empirical value, threshold of detection the empiricalbyvalue, is equal 5000, proposed [19]. In of addition, which is equal to 5000,isas proposed [19]. which In addition, in to order toas analyze the by influence orbital in order to analyze the influence of orbital maneuver on precise orbit determination of BeiDou maneuver on precise orbit determination of BeiDou satellites, the validity of the maneuver information satellites, validity thebemaneuver detected this section will be analyzed in the detected inthe this sectionofwill analyzed information in the last part of this in section from the perspective of precise last part of this section from the perspective of precise orbit determination. orbit determination. Since Since the the solutions solutions of of precise precise orbit orbit determination determination often often use use single-day, single-day, three-day three-day and and seven-day seven-day broadcast ephemeris products as the initial orbit, the size of the windows selected for use broadcast ephemeris products as the initial orbit, the size of the windows selected for use in in the the test test were single-day, three-day and seven-day. were single-day, three-day and seven-day. Program one: one:The Thetesting testing detecting ability of orbit BDS maneuvering orbit maneuvering of themethod robust Program forfor thethe detecting ability of BDS of the robust method proposed in thisFor paper. For GEO-03, thefrom time5 is from 52015 January to 72015 January 2015 proposed in this paper. GEO-03, the time is January 00:002015 to 700:00 January 23:00; for 23:00; for IGSO-03, the time is from 9 January 2015 05:00 to 10 January 2015 05:00. The data comes IGSO-03, the time is from 9 January 2015 05:00 to 10 January 2015 05:00. The data comes from the BDS from thebroadcast BDS merged broadcast ephemeris provided by iGMAS. merged ephemeris provided by iGMAS.

2

x 10

method one(tref=36)

4

0

-2

First-order difference rms3 T

-4 -6 -8 0

8

x 10

10

20

30

40

50

60

70

80

method two(tref=36)

4

rm s3 /m

6 rms3

4

T

2 0 0

10

20

30

40 Time/h

50

(a) The results when tref  36 Figure 8. Cont.

60

70

80

Sensors 2017, 17, 1129 Sensors 2017, 17, 1129

F i r s t - o r d e r d if f e r e n c e r m s3 / m

10 of 17 10 of 17

6

x 10

method one(tref=10)

4

First-order difference rms3 T

4

2 1129 Sensors 2017, 17,

10 of 17

F i r s t - o r d e r d if f e r e n c e r m s3 / m

0

-2 0

8

x 10

4

6

x 10

method one(tref=10)

4

First-order difference rms3

10

4

20

r m s3 / m

50

60

T

70

8

method two(tref=10) rms3 T

10

x 10

20

30

40

50

60

70

80

method two(tref=10)

4

rms3 T

r m s3 / m

6 4

0 0

80

0

4 2

40

2

-2 0

6

30

10

20

30

40 Time/h

2 0 0

50

(b) The results when tref  10

10

20

30

40 Time/h

50

60

60

70

70

80

80

Figure for method methodone one((∇rms )/method two(rms ( rms3 (t ) ) when tref  36 and tref  10 , rms3((tt))/method Figure 8. 8. The The results results for 3 results when two tref  10 3 (t )) when tre f = 36 and tre f = 10, (b) The respectively. respectively. Figure 8. The results for method one ( rms3 ( t ) )/method two ( rms3 (t ) ) when tref  36 and tref  10 , respectively.

(1) Testing and analysis in GEO-03(PRN3). The testing results of methods one and two are shown (1) Testing and analysis in GEO-03(PRN3). The testing results of methods one and two are shown in in Figure 8 (X-axis represents time, unit hours; Y-axis represents value calculated by (1) Testing and analysis in GEO-03(PRN3). The is testing results of methods one and two are shown Figure 8 (X-axis represents time, unittime, is hours; represents value value calculated by Equation (3), 8 meters); (X-axis represents unit2015 isY-axis hours; represents calculated by Equation in (3),Figure unit is from 5 January 9:00 Y-axis to 5 January 2015 14:00, the health flag of unit is meters); from 5 January 2015 to 5 January 2015 14:00,2015 the14:00, health theofbroadcast Equation (3), unit is meters); from9:00 5 January 2015 9:00 to 5 January the flag healthofflag the broadcast ephemeris showed that GEO-03 is in an unhealthy state; on 5 January 2015, theshowed broadcastthat ephemeris showed thatunhealthy GEO-03 is in an unhealthy state; 2015, on 5 January ephemeris GEO-03 is in an state; on 5 January among2015, the products among the products of precise orbits, notcalculate calculate coordinates of GEO-03; at among the products of precise orbits,GFZ GFZ did did not the the coordinates of GEO-03; at of precise orbits, GFZ did not calculate the coordinates of GEO-03; at station WUH1, the O-C station WUH1, the O-C of GEO-03 Figure 9 (X-axis represents unit is station WUH1, thevalues O-C values of GEO-03are areshown shown inin Figure 9 (X-axis represents time, unittime, is values of hours; GEO-03 arerepresents shown in Figure 9 (X-axisdistance represents time, unitdistance is hours; Y-axis represents Y-axis station-satellite observed of L3, unit hours; Y-axis represents that that station-satellite distanceminus minus observed distance of isL3, unit is that station-satellite distance minus observed distance of L3, unit is meters); meters); meters); 5.5

x 10

5

5

5.4 5.35 O-C/m

5.4

O-C/m

x 10

5.45

5.45

5.35

5.5

5.3 5.25

5.3

5.2

5.25

5.15 5.1

5.2

5.05 0

5.15

1

2

3

5.1 5.05 0

4

5

6

7

8

9

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Time/h

Figure 9. Results O-C of of GEO-03 at WUH1. Figure 9. Results ofofO-C GEO-03 at WUH1. 1

2

3

4

5

6

7

8

9

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

From Figure 8 and the other three verification methods for GEO-03, the testing results from 5 Time/h at 14:00–15:00, shown in Table 1. methods for GEO-03, the testing results From January Figure2015, 8 and the otherarethree verification

from

9. Results of O-C 5 January 2015, at 14:00–15:00,Figure are shown in Table 1.of GEO-03 at WUH1.

Table 1. Results of GEO-03 maneuvering detection at 14:00–15:00.

Method One Method Two (Conventional The Health Flag of Precise Orbit O-C Value at From Figure 8 and the three verification methods for GEO-03, the testing results Table 1. other Results of GEO-03 maneuvering detection at 14:00–15:00. (New Method) Method) Broadcast Ephemeris Products of GFZ WUH1 January 2015,Maneuvering at 14:00–15:00, are shown in Table 1. Unhealthy Maneuvering No GEO-03 There is a jump

Method One (New Method) Maneuvering

Method One (New Method) Maneuvering

Method Two The Health Flag of Precise Orbit (Conventional Method) Broadcast Ephemeris of GFZ Table 1. Results of GEO-03 maneuvering detectionProducts at 14:00–15:00. Maneuvering Unhealthy No GEO-03

Method Two (Conventional Method) Maneuvering

The Health Flag of Broadcast Ephemeris Unhealthy

Precise Orbit Products of GFZ No GEO-03

from 5

O-C Value at WUH1 There is a jump

O-C Value at WUH1 There is a jump

Sensors 2017, 17, 1129 Sensors 2017, 17, 1129 Sensors 2017, 17, 1129

1117 of 17 11 of 11 of 17

First-order difference rms8 T

5

10

-4 0 7

b1

4

x 10

15

20

25

10 15 method two(tref=10)

20

method two(tref=10)

7

T

rms8/m

1

0 0

rms8 T

2 1 5

10

0 0

5

15 Time/h 10 15 Time/h

20

25 20

3

x 10 a2

3

method one(tref=10)

5

x 10

method one(tref =10) difference rms First-order

5

8

T

a2

First-order difference rms8

2

T

2

1 0

5

rms8

b1

3

2

25

rms8/m

x 10

5

First-order difference rms8 /m

T

a1

-2

3 rms8/m

8

0

-2

4

method one(tref =10) difference rms First-order

7

2

0

-4 0

x 10

4

2

method one(tref=10)

7

a1

First-order difference rms8 /m

x 10

1 14 0

x 10 b2

x 10

4

5

3

4

2

3

1

2

0

1 14 0

25

16 14

18

20

22

16 18 20 method two(tref=10)

5

rms8/m

4

First-order difference rms8 /m

First-order difference rms8 /m

The results in Table 1 show that on 5 January 2015, at 14:00–15:00, both methods detect that The results in in Table 1 show on 5flag January 2015, atephemeris 14:00–15:00, detect that results boththat methods detect GEO-03The maneuvering may occur: thethat health of broadcast shows GEO-03 is in that a GEO-03 maneuvering may occur: the health flag of broadcast ephemeris shows that GEO-03 is GEO-03 state, maneuvering mayprevious occur: the health flag of health broadcast thatisGEO-03 inin a unhealthy but in the five hours, the flagephemeris shows thatshows GEO-03 also inisan a unhealthy state, but in the previous five hours, the health flag shows that GEO-03 is also in an 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 unhealthy state; GFZ might think that is is inin anan unhealthy state; and, there is a is jump initthe unhealthy GFZ might think thatGEO-03 GEO-03 unhealthy state; and, there a jump inO-C the O-C value at WUH1, which may represent that GEO-03 is in an unhealthy state. Furthermore, can value at WUH1, which may represent that GEO-03 is in an unhealthy state. Furthermore, it can be value atfrom WUH1, which may GEO-03 in an unhealthy state.isFurthermore, it can be O-C concluded Figure 8 that therepresent detectionthat result of theisconventional method affected by the concluded from Figure 8and that theresult detection result of the conventional is affected the selected be concluded from Figure 8the that the detection result of thefrom conventional method is affected by the selected reference epoch, of the epoch with far themethod reference epoch isbyobviously reference epoch, and the result of the epoch with far from the reference epoch is obviously affected selected epoch,and andtrend the result with far result from the reference epoch is obviously affected by reference the periodicity term, of sothe thatepoch the detection is not completely correct. But by the periodicity and trend term, so that the detection result is not completely correct. But the new 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: method not affected by these by factors: the new is method is not affected these factors: (2) Testing and analysis in IGSO-03(PRN8). The testing results of methods one and two are shown (2) and analysis in testing of one and two are Figure From 9 January 2015 14:00 toThe 9 January 2015 20:00, the health broadcast (2)in Testing Testing10. and analysis in IGSO-03(PRN8). IGSO-03(PRN8). The testingresults results of methods methods oneflag andof two are shown shown in 10. January to the flag of ephemeris showed that99IGSO-03 is in 14:00 an unhealthy state,2015 while20:00, at 11:00, the health is in a in Figure Figure 10. From From January2015 2015 14:00 to 99 January January 2015 20:00, the health health flagflag of broadcast broadcast ephemeris that isisin while at the health flag healthy state. showed On 9 January 2015, among theunhealthy products state, of precise orbits, GFZ did not calculate ephemeris showed thatIGSO-03 IGSO-03 inan an unhealthy state, while at 11:00, 11:00, the health flag is is in in aa healthy state. On 9 January 2015, among the products of precise orbits, GFZ did not calculate the healthy coordinate IGSO-03. At 2015, station BJF1,the theproducts O-C values of IGSO-03 shown in Figure 11the state.ofOn 9 January among of precise orbits,are GFZ did not calculate the coordinate oftime, IGSO-03. At hours; station BJF1, the O-C values ofstation-satellite IGSO-03 are shown in Figure 11 (X-axis represents unit is Y-axis represents distance minus coordinate of IGSO-03. At station BJF1, the O-C values ofthat IGSO-03 are shown in Figure 11 (X-axis (X-axis represents time, unit is hours; Y-axis represents that station-satellite distance minus observed distance B1, unit is meters); represents time,ofunit is hours; Y-axis represents that station-satellite distance minus observed observed distance ofisB1, unit is meters); distance of B1, unit meters);

24 22

method two(tref=10)

5

24 rms8 T

b2

rms8 T

16 14

18 20 Time/h 18 20 Time/h

16

22

24 22

24

Figure 10. (a1) The results for method one ( rms8 ( t ) ), (a2) The results for partial enlarged (a1), and Figure 10. 10. (a1) (a1)The Theresults results method one (∇rms results for partial enlarged (a1), 8 ( t )), (a2) forfor TheThe results for partial rms 8 ( t ) ), (a2) (b1)Figure The results for method two ( method (b2)( The results for partial enlarged (b1). enlarged (a1), and rms8 (t ) ),one and (b1) The results for method two (rms ( t ) ), (b2) The results for partial enlarged (b1). (b1) The results for method two ( rms (t ) ),8 (b2) The results for partial enlarged (b1). 8

4.5

3.5

4

3

3.5

2.5

3

2

2.5

O-C/m

O-C/m

4

6

1.5

2

1

1.5

0.5

1

0

0.5

-0.5 0

5

x 10 x 10

5

x 10 1

0 0 -1 -1 -2

-3

-2

-3 -4 -4

0 1 2 -0.5 0

1

6

O -C/m

x 10

O -C/m

4.5

3 1

4 2

5 3

6 4

7 5

8 6

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Time/h 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Time/h

-5

19 -5

Time/h 19

Time/h

20

21 20

21

(a) The O-C values of IGSO-03 at BJF1 (b) The results for partial enlarged (a) (a) The O-C values of IGSO-03 at BJF1 (b) The results for partial enlarged (a) Figure 11. Results of O-C of IGSO-03 at BJF1 (The graph on the right is a partial enlarged view of the Figure 11. Results of O-C of IGSO-03 at BJF1 (The graph on the right is a partial enlarged view of the left). Figure 11. Results of O-C of IGSO-03 at BJF1 (The graph on the right is a partial enlarged view of the left). left).

From Figure 10 and the other three verification methods for IGSO-03, from 9 January 2015, at Fromand Figure 10 and the verification methods for2 IGSO-03, from 9 January 2015, at 10:00–11:00 19:00–20:00, the other testingthree results are shown in Tables and 3, respectively. 10:00–11:00 and 19:00–20:00, the testing results are shown in Tables 2 and 3, respectively.

Sensors 2017, 17, 1129

12 of 17

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. Table 2. Results of IGSO-03 maneuvering detection at 10:00–11:00. Method One Method Two (New Method) (Conventional Method) Sensors 2017, 17, 1129 Anomalies

The Health Flag of Broadcast Ephemeris

Precise Orbit Products of GFZ

O-C Value at BJF1 12 of 17

Healthy

No IGSO-03

There is a jump

Anomalies

Table 2. Results of IGSO-03 maneuvering detection at 10:00–11:00.

Method One (New Method) Anomalies

Two The Health Flagdetection of Orbit Table 3.Method Results of IGSO-03 maneuvering atPrecise 19:00–20:00. (Conventional Method) Anomalies

Method One (New Method)

Method Two (Conventional Method)

Maneuvering

Maneuvering

Broadcast Ephemeris Healthy

O-C Value at BJF1 There is a jump

Products of GFZ No IGSO-03

The health flag of Broadcast Ephemeris

Precise Orbit Products of GFZ

Unhealthy

No IGSO-03

O-C Value at BJF1

Table 3. Results of IGSO-03 maneuvering detection at 19:00–20:00.

Method One (New Method) Maneuvering

Method Two (Conventional Method) Maneuvering

The health flag of Broadcast Ephemeris Unhealthy

There is a jump

Precise Orbit Products of GFZ No IGSO-03

O-C Value at BJF1 There is a jump

F i r s t - o r d e r d i f f e r e n c e r m s1 / m

From 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 From Figure 10, it can be seen that anomalies do not affect the detection of maneuvering. From IGSO-03 maneuvering occur: the health flag of broadcast ephemeris both shows that IGSO-03 is in an the results in Table 3,may it can be seen that on 9th January 2015, at 19:00–20:00, methods detect that unhealthy state, while in the previous five hours and the later one hour, the health flag shows IGSO-03 maneuvering may occur: the health flag of broadcast ephemeris shows that IGSO-03 is in an that IGSO-03 is alsostate, in anwhile unhealthy state; GFZ think IGSO-03 is in an unhealthy in the previous fivemay hours andthat the later one hour, theunhealthy health flag state; showsand, that there IGSO-03 is also an unhealthy state; GFZ think that is is in in anan unhealthy state; and, is a jump in the O-Cinvalue at BJF1, which maymay represent thatIGSO-03 IGSO-03 unhealthy state. there is a jump the O-C for value at detecting BJF1, whichability may represent IGSO-03 is in an state. Program two:inTesting the of BDS that orbit anomalies byunhealthy the robust method 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 in this paper. For GEO-01, the time was from 4th January 2015 00:00 to 10th January 2015 23:00.proposed The data comes from the BDS merged broadcast ephemeris provided by iGMAS. 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 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 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 orbit,GFZ GFZ calculated coordinates of GEO-01. in a healthy state. Among theproducts products of of precise precise orbit, calculated the the coordinates of GEO-01. At station BJF1,BJF1, the O-C values of GEO-01 are are shown in Figure 13 13 (X-axis represents At station the O-C values of GEO-01 shown in Figure (X-axis representstime, time,unit unitisishours; Y-axis represents that station-satellite distancedistance minus minus observed distance of B1, is meters). hours; Y-axis represents that station-satellite observed distance of unit B1, unit is meters). 1

x 10

method one(tref=84)

5

First-order difference rms1 T

0.5 0

-0.5 -1 0

6

x 10

20

40

60

80

100

120

140

160

180

120

140

160

180

method two(tref=84)

4

rms1 T

r m s1 /m

4

2

0 0

20

40

60

80

Time/h

100

Figure Theresults resultsfor formethod method one( two( rms1 (t ) ).(t)). 1 ( t )()/method Figure 12.12. The one(∇rms rms 1 t )/method two(rms 1

Sensors 2017, 17, 1129

13 of 17

Sensors 2017, 17, 1129

13 of 17 5

x 10 1.244 1.243 1.242

O -C /m

1.241 1.24

1.239 1.238 1.237 1.236 1.235

18

Time/h

19

Figure 13. Results of O-C of GEO-01 at BJF1.

Figure 13. Results of O-C of GEO-01 at BJF1.

From Figure 12 and other three verification ways for GEO-01, on 8 January 2015, at 18:00–19:00, 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.

the testing results are shown in Table 4. Table 4. Results of GEO-01 anomalies detection at 18:00–19:00. Method One (New Method) Anomalies Method One

TableMethod 4. Results detection at 18:00–19:00. Two of GEO-01 anomalies The Health Flag of Precise Orbit (Conventional Method) Anomalies Method Two

Broadcast Ephemeris Healthy The Health Flag of

Products of GFZ Have GEO-01 Precise Orbit

O-C Value at BJF1 There a jumpat O-CisValue

(New Method) (Conventional Method) Broadcast Ephemeris Products of GFZ BJF1 From the results in Table 4, it can be seen that on 8 January 2015, at 18:00–19:00, both methods Anomalies Anomalies Healthy Have GEO-01 There is a jump 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; 4, and, therebeis seen a jump in the value at BJF1,atwhich may represent that From the results in Table it can that on O-C 8 January 2015, 18:00–19:00, both methods GEO-01 in an unhealthy state.may Moreover, canare be seen fromby Figure 12 thaterrors the results at both ends the detected thatisGEO-01 anomalies occur itthat caused decoding or other reasons: of flag the detection window due to the impact the periodicity and trend termGFZ is clearly health of broadcast ephemeris show thatofGEO-01 is in a healthy state; may wrong. think that GEO-01 From Figure 10, it can be seen that maneuvering does not affect detection of anomalies. From is in a healthy state; and, there is a jump in the O-C value at BJF1, which may represent that GEO-01 the results in Table 2, it can be seen that on 9 January 2015, at 10:00–11:00, both methods detect that is in an unhealthy state. Moreover, it can be seen from Figure 12 that the results at both ends of the IGSO-03 may occur anomalies that are caused by decoding errors or other reasons: the health flag of detection window due to the impact of the periodicity and trend term is clearly wrong. broadcast ephemeris show that IGSO-03 is in a healthy state; GFZ may think that IGSO-03 is in an From Figure 10, it can beisseen that affect detection anomalies. unhealthy state; and, there a jump in maneuvering the O-C value atdoes BJF1,not which may representof that IGSO-03 isFrom the results in Table 2, it can be seen that on 9 January 2015, at 10:00–11:00, both methods detect that in an unhealthy state. IGSO-03 may occur anomalies are caused by decoding errorsmaneuvering or other reasons: the health flag of Program three: In orderthat to further test the effect of orbital information on the preciseephemeris orbit determination, the orbitis maneuvering information is used orbit broadcast show that IGSO-03 in a healthy state; GFZ may think for thatprecise IGSO-03 is in an determination. As there neither nor in WHU theatprecise the two satellites when these is in unhealthy state; and, isGFZ a jump the provides O-C value BJF1, orbits whichofmay represent that IGSO-03 two satellites maneuvered their orbits, orbit overlap comparison is used to analyze the precision of an unhealthy state. the orbit determination.

Program three: In order to further test the effect of orbital maneuvering information on the (1) orbit GEO-03 precise orbitthe maneuvering test. In thisinformation test, the BDSismeasured data of 20 stations from precise determination, orbit maneuvering used for precise orbit determination. MGEX and 10 stations from iGMAS from 3 to 6 January 2015 were used to obtain the final two As neither GFZ nor WHU provides the precise orbits of the two satellites when these precision orbit using three-day long arc POD. The comparison strategy is shown in Figure 14,of the satellites maneuvered their orbits, orbit overlap comparison is used to analyze the precision and the statistical results of orbit overlap comparison in radial direction, along direction, cross orbit determination. direction, 1D-RMS are shown in Table 5;

(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;

Sensors 2017, 17, 1129 Sensors 2017, 17, 1129

14 of 17 14 of 17

Sensors 2017, 17, 1129

14 of 17

Figure14. 14.Comparison ComparisonStrategy Strategyof ofGEO-03. GEO-03. Figure of GEO-03. TableFigure 5. RMS14. ofComparison orbit overlapStrategy comparison at 5 January (unit: m). Table 5. RMS of orbit overlap comparison at 5 January (unit: m). Satellite GEO-01 GEO-02 Satellite GEO-03 GEO-01 GEO-04 GEO-02 GEO-05 GEO-03 IGSO-01 GEO-04 IGSO-02 GEO-05 IGSO-03 IGSO-01 IGSO-04 IGSO-02 IGSO-05 IGSO-03

(2)

IGSO-04 IGSO-05 IGSO-03

PRN Radial Along Cross Table 5. RMS of orbit overlap comparison at 5 January (unit: m). 1 Satellite 0.1389 0.1695 PRN Radial 0.2016 Along Cross 1D-RMS 2 1.4875 1.7481 0.2402 PRN Cross GEO-01 Radial 1 0.1389 Along 0.2016 0.1695 0.2634 0.7614 1.8901 0.3392 13 GEO-02 0.1389 0.2016 0.1695 2 1.4875 1.7481 0.2402 2.3079 0.1435 0.3520 0.0743 24 0.2402 GEO-03 1.4875 3 0.7614 1.7481 1.8901 0.3392 2.0657 0.1576 0.4225 0.2662 35 0.3392 GEO-04 0.7614 4 0.1435 1.8901 0.3520 0.0743 0.3873 0.0828 0.1644 0.1258 46 0.1435 0.3520 0.0743 GEO-05 5 0.1576 0.4225 0.2662 0.5236 0.1170 0.2529 0.7386 57 0.1576 0.4225 0.2662 IGSO-01 6 0.0828 0.1644 0.1258 0.2230 0.1570 0.4986 0.4792 68 0.0828 0.1644 0.1258 IGSO-02 0.0787 7 0.1170 0.0916 0.2529 0.7386 0.7894 9 0.1428 7 0.1170 0.2529 0.7386 IGSO-03 8 0.1570 0.4986 0.4792 0.7091 10 0.0722 0.2861 0.6542 8 0.1570 0.4986 0.4792

IGSO-04

9

0.0787

0.0916

9 0.0787 0.0916 IGSO-05 10 0.0722 0.2861 10 0.0722 0.2861 precise orbit maneuvering test. In this test,

0.1428 0.6542

0.1870

0.1428 0.7177 0.6542 measured

1D-RMS 0.2634 2.3079 1D-RMS 2.0657 0.2634 0.3873 2.3079 0.5236 2.0657 0.2230 0.3873 0.7894 0.5236 0.7091 0.2230 0.1870 0.7894 0.7177 0.7091

0.1870 from0.7177 7 January

the BDS data to 10 2015 were used to obtain the final precision orbit. The comparison strategy is shown in (2) IGSO-03 precise orbitstatistical maneuvering test. In thistest, test,the theBDS BDSmeasured measured data from 7 10 January to (2) IGSO-03 precise maneuvering test. this data from 7 to January Figure 15, and the results ofIn orbit mutual difference of overlapping arc in radial 10 2015 were to obtain the final precision orbit. The comparison strategy is Figure shown15, in 2015 were usedused to obtain the final precision orbit. The strategy shown in direction, along direction, cross direction, 1D-RMS arecomparison shown in Table 6; is Figure 15, and the statistical results of orbit mutual difference of overlapping arc in radial and the statistical results of orbit mutual difference of overlapping arc in radial direction, along direction, cross direction, 1D-RMS are 6; shown in Table 6; direction,along cross direction, 1D-RMS are shown in Table

Figure 15. Comparison Strategy of IGSO-03. Figure Strategy Figure 15. Comparison Comparison Strategy of of IGSO-03. Table 6. RMS of 15. orbit overlap comparison atIGSO-03. 9 January (unit: m). Satellite GEO-01 GEO-02 Satellite GEO-03 GEO-01 GEO-04 GEO-02 GEO-05 GEO-03 IGSO-01 GEO-04 IGSO-02 GEO-05 IGSO-03 IGSO-01 IGSO-04 IGSO-02 IGSO-05 IGSO-03 IGSO-04 IGSO-05 From the

PRN Radial Along Cross Table Table6.6.RMS RMSof oforbit orbitoverlap overlapcomparison comparisonat at99January January(unit: (unit:m). m). 1 0.0554 0.0783 0.0458 2 0.0686 0.1772 0.1112 PRN Radial Along Cross Satellite 0.0940 PRN Radial 0.1030 Along Cross 0.0776 1D-RMS 13 0.0554 0.0783 0.0458 4 0.0288 0.1301 0.0370 2 GEO-01 0.0686 1 0.0554 0.1772 0.0783 0.0458 0.1112 0.1063 0.0565 0.0530 0.1104 35 GEO-02 0.0940 2 0.0686 0.1030 0.1772 0.1112 0.0776 0.2202 0.0328 0.1201 0.3887 46 GEO-03 0.0288 3 0.0940 0.1301 0.1030 0.0776 0.0370 0.1596 0.1074 0.1651 0.5064 57 GEO-04 0.0565 4 0.0288 0.0530 0.1301 0.0370 0.1104 0.1383 0.4460 0.3687 0.2069 68 GEO-05 0.0328 5 0.0565 0.1201 0.0530 0.1104 0.3887 0.1349 0.0775 0.1115 0.5155 79 IGSO-01 0.1074 6 0.0328 0.1651 0.1201 0.3887 0.5064 0.4082 10 0.0935 0.0888 0.0943 8 IGSO-02 0.4460 7 0.1074 0.3687 0.1651 0.5064 0.2069 0.5434 9

IGSO-03 0.0775 8

0.4460 0.1115 0.3687

0.2069

0.5155 0.6145

1D-RMS 0.1063 0.2202 1D-RMS 0.1596 0.1063 0.1383 0.2202 0.1349 0.1596 0.4082 0.1383 0.5434 0.1349 0.6145 0.4082 0.5331 0.5434 0.1598 0.6145 0.5331

10 of 0.0888 0.0943 0.1598 IGSO-04 9 IGSO-03 0.0775precise 0.1115orbit0.5155 0.5331tests, we can results GEO-030.0935 and maneuvering see that it is IGSO-05 10 0.0935 0.0888 0.0943 0.1598 difficult to provide precise orbit determination results for GEO-03 and IGSO-03 when there is no From themaneuvering results of GEO-03 and IGSO-03 precise orbit maneuvering tests,information we can see that it is satellite orbit information. However, when using the maneuver detected difficult to provide precise orbit determination results for GEO-03 and IGSO-03 when there is in this paper, the precise orbit determination results of these two satellites not only can no be 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

Sensors 2017, 17, 1129

15 of 17

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. 4. 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.

Sensors 2017, 17, 1129

16 of 17

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. Acknowledgments: This work was supported by the Collaborative Precision Positioning Project funded by the Ministry of Science and Technology of China (No. 2016YFB0501900) and China Natural Science Funds (No. 41231064, 41674022, 41574015). The authors would like to thank the IGS MGEX and Chinese iGMAS networks for providing the BeiDou data. Thanks to Baocheng Zhang and Lidan Zhou for giving careful guidance in modifying this article. Author Contributions: Fei Ye and Yunbin Yuan provided the initial idea for this study; Fei Ye and Bingfeng Tan conceived and designed the experiments; Fei Ye, Yunbin Yuan and Jikun Ou analyzed the experiment results; Fei Ye and Bingfeng Tan wrote the paper. Conflicts of Interest: The authors declare no conflict of interest.

References 1. 2. 3.

4.

5. 6. 7. 8. 9. 10. 11. 12.

13. 14. 15.

16.

Yang, Y.; Li, J.; Xu, J.; Tang, J.; Guo, H.; He, H. Contribution of the compass satellite navigation system to global PNT users. Chin. Sci. Bull. 2011, 56, 2813–2819. [CrossRef] Camacho-Lara, S. Current and Future GNSS and Their Augmentation Systems. In Handbook of Satellite Applications; Springer: New York, NY, USA, 2013; pp. 617–654. Montenbruck, O.; Hauschild, A.; Steigenberger, P.; Hugentobler, U.; Teunissen, P.; Nakamura, S. Initial assessment of the COMPASS/BeiDou-2 regional navigation satellite system. GPS Solut. 2013, 17, 211–222. [CrossRef] Li, X.; Zhang, X.; Ren, X.; Fritsche, M.; Wickert, J.; Schuh, H. Precise positioning with current multi-constellation global navigation satellite systems: GPS, GLONASS, Galileo and BeiDou. Sci. Rep. 2015, 5, 8328. [CrossRef] [PubMed] Tan, B.; Yuan, Y.; Zhang, B.; Hsu, H.Z.; Ou, J. A new analytical solar radiation pressure model for current BeiDou satellites: IGGBSPM. Sci. Rep. 2016, 6, 32967. [CrossRef] [PubMed] Shi, C.; Zhao, Q.; Hu, Z.; Liu, J. Precise relative positioning using real tracking data from COMPASS GEO and IGSO satellites. GPS Solut. 2013, 17, 103–119. [CrossRef] Cao, F.; Yang, X.; Li, Z.; Sun, B.; Kong, Y.; Chen, L.; Feng, C. Orbit determination and prediction of GEO satellite of BeiDou during repositioning maneuver. Adv. Space Res. 2014, 54, 1828–1837. [CrossRef] Wang, B.; Lou, Y.; Liu, J.; Zhao, Q.; Su, X. Analysis of BDS satellite clocks in orbit. GPS Solut. 2016, 20, 783–794. [CrossRef] Byun, S.H. Satellite orbit determination using triple-differenced GPS carrier phase in pure kinematic mode. J. Geod. 2003, 76, 569–585. [CrossRef] Steigenberger, P.; Hugentobler, U.; Hauschild, A.; Montenbruck, O. Orbit and clock analysis of Compass GEO and IGSO satellites. J. Geod. 2013, 87, 515–525. [CrossRef] Zhao, Q.; Wang, C.; Guo, J.; Liu, X. Assessment of the Contribution of BeiDou GEO, IGSO, and MEO Satellites to PPP in Asia—Pacific Region. Sensors 2015, 15, 29970–29983. [CrossRef] [PubMed] Geng, T.; Su, X.; Fang, R.; Xie, X.; Zhao, Q.; Liu, J. BDS Precise Point Positioning for Seismic Displacements Monitoring: Benefit from the High-Rate Satellite Clock Corrections. Sensors 2016, 16, 2192. [CrossRef] [PubMed] Liu, T.; Yuan, Y.; Zhang, B.; Wang, N.; Tan, B.; Chen, Y. Multi-GNSS precise point positioning (MGPPP) using raw observations. J. Geod. 2016. [CrossRef] Davis, T.; Baker, M.T.; Belchak, T.; Larsen, W. XSS-10 micro-satellite flight demonstration program. Proc. SPIE 2004, 5419, 16–25. Jacobovits, A.; Vaneck, T. AeroAstro’s Escort–a microsatellite for on-orbit inspection of space assets. In Proceedings of the 17th AIAA/USU Annual Conference on Small Satellites, SSC03-IV-7, Logan, UT, USA, 12 August 2003. Zhang, J.; Chen, F.; Fu, Q.; Xiao, H. Moving Target Recognition Based on HMM. J. Natl. Univ. Def. Technol. 2003, 25, 51–55.

Sensors 2017, 17, 1129

17. 18. 19. 20. 21. 22. 23. 24. 25.

17 of 17

Che, R.; Zhang, H. Relative Orbit Design of a Chaser Tracking a Non-cooperative Target in Space. Aerosap. Control 2006, 24, 40–45. Liu, Y.; Chang, Z.; He, F.; Guo, R. GEO Satellite Abnormity Recognition Based on Wavelet Analysis. J. Geod. Geodyn. 2010, 30, 73–76. Yan, X.; Huang, G.; Zhang, R.; Zhang, Q. A Method Based on Broadcast Ephemeris to Detect BDS Satellite Orbital Maneuver. J. Navig. Position 2015, 3, 35–38. Yuan, Y.; Ou, J. Auto-covariance estimation of variable samples (ACEVS) and its application for monitoring random ionospheric disturbances using GPS. J. Geod. 2001, 75, 438–447. [CrossRef] Stigler, S.M. Do robust estimators work with real data? Ann. Stat. 1977, 5, 1055–1098. [CrossRef] Rousseeuw, P.J. Least median of squares regression. J. Am. Stat. Assoc. 1984, 79, 871–880. [CrossRef] Rousseeuw, P.J.; Wagner, J. Robust regression with a distributed intercept using least median of squares. Comput. Stat. Data Anal. 1994, 17, 65–76. [CrossRef] Yang, Y. Robust estimation of geodetic datum transformation. J. Geod. 1999, 73, 268–274. [CrossRef] Yang, Y.; Xu, J. GNSS receiver autonomous integrity monitoring (RAIM) algorithm based on robust estimation. Geod. Geodyn. 2016, 7, 117–123. [CrossRef] © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).