Performance analysis of SBAS ephemeris corrections and integrity algorithms in China region

Satellite Based Augmentation Systems (SBASs) improve the positioning accuracy and integrity by broadcasting to the civil aviation community the corrections and integrity parameters. A snapshot algorithm based on the minimum variance estimation is investigated in this study to calculate the satellite clock and orbit corrections. A chi-square test is performed on the remaining errors in the corrected ephemeris to guarantee the integrity. User Differential Range Error (UDRE) and scaling matrix contained in Message Type 28 are derived using the covariance information based on the assumption that one of the reference stations failed. A software package is developed and applied in the real data collected at 26 stations. International GNSS (Global Navigation Satellite System) Service (IGS) precise clock and orbit products are taken as the references to assess the accuracy of corrections. For both Global Positioning System (GPS) and BeiDou Navigation Satellite System (BDS), the range accuracy of 0.10 m can be achieved with the employment of the derived corrections. No obvious performance difference between GPS and BDS is found. UDREs for all visible satellites are generated with the maximum index of 12 and minimum index of 3. The geometric range differences calculated with IGS precise products and broadcast ephemeris are employed to assess the integrity of UDRE. It is found that the UDRE is able to bound the residuals with 99.9% confidence which meet the requirement of aviation users. With ionospheric delay corrected by Global Ionosphere Map (GIM), the positioning accuracy of 0.98 m with GPS corrections and 0.80 m with multi-constellation augmentation can be achieved which indicates a significant improvement of GPS standalone results.


Introduction
In order to meet the accuracy, continuity, availability, and especially integrity requirements of civil aviation navigation based on Global Navigation Satellite System (GNSS), various Satellite Based Augmentation Systems (SBASs) have been developed, such as Wide Area Augmentation System (WAAS) in the USA, European Geostationary Navigation Overlay Service (EGNOS) in Europe, GPS (Global Positioning System) Aided Geo Augmented Navigation (GAGAN) in India, and MTSAT (Multi-functional Transport Satellite) Satellite-based Augmentation System (MSAS) in Japan. BeiDou Satellite Based Augmentation System (BDSBAS) in China is also under development. Integrity refers to the reliability and trustworthiness of the information provided by a navigation system and to the system's ability to deliver timely warnings to users when the system should not be used for navigation because of signal corruption or some other errors or failures in the system (Sparks et al., 2011). SBAS broadcasts the corrections of GNSS satellite ephemeris and ionosphere delay along with the associated confidences (Jin et al, 2020). The confidences of the ephemeris corrections, User Differential Range Error (UDRE), and the scaling matrix contained in Message Type 28 (MT28)

Open Access
Satellite Navigation https://satellite-navigation.springeropen.com/ *Correspondence: jinbiao366788@126.com 1 Space Star Technology Co., Ltd., Beijing 100194, China Full list of author information is available at the end of the article are transmitted to users for computing the positioning protection levels (RTCA, 2013).
Several methods have been developed to estimate the errors in broadcast ephemeris or to determine the precise GNSS clocks and orbits (Griffiths & Ray, 2009;Dach et al., 2009;Guo et al., 2016). For the implementation of SBAS, detailed snapshot methods to solve for the clock and orbit corrections were proposed in (Enge et al., 1996;Tsai, 1999). Theoretical descriptions of ephemeris corrections calculated by a minimum mean square estimator were introduced in Wu & Peck, 2002). A correction generation model by combining pseudo range and phase data was presented in (Chen, Yang, et al., 2017), where the phase smoothed code observation was applied to define the datum and epoch differenced phase was adopted to define epoch wise variation of the satellite clock and orbit corrections. Chen, Huang, et al. (2017) computed clock corrections by combining prior knowledge of satellite clock prediction and parameter-fitting errors with real time measurements, which could reduce the User Range Error (URE) in comparison with applying WAAS corrections. Zhao et al. (2020) reported that it was difficult to obtain the orbit that is more accurate than the broadcast orbit using regional monitoring stations, and therefore suggested a SBAS correction generation approach with only clock corrections calculated. Satellite dynamic model was introduced by Hughes company to determine the WAAS ephemeris corrections (Ceva, 1997). A new module for GNSS satellite orbit determination and clock synchronization with the ability of detecting and reacting to a satellite maneuver was developed for EGNOS (Rouanet et al., 2016). The dynamic model method provides higher accuracy and can separate satellite orbit and clock errors (Ceva, 1997), but its main problem lies in the computational complexity.
SBAS provides UDRE for each satellite to inform users the accuracy level of the ephemeris corrections (Walter et al., 2018). UDRE is a scalar value and it is conservative to take the largest value to protect all users within the service volume. Instead, the relative covariance of the orbit and clock corrections contained in MT28 was proposed and used along with the UDRE. Users can reconstruct their location specific error bound rather than applying the largest one (Authié et al., 2017;Walter et al., 2001). UDRE generation based on the residuals between the computed distances and observations was introduced in (Peck & Tekawy, 1997). UDRE and MT28 calculation based on the covariance solution with the assumption that one of the monitoring stations failed was proposed in (Wu & Peck, 2002). A detailed derivation of UDRE and MT28 for dual frequency SBAS was presented in  considering the range error induced by antenna bias and the difference between estimated and broadcast corrections. Blanch et al. (2013) evaluated the potential of changes in WAAS monitoring algorithms among which a new clock and ephemeris algorithm was engaged to mitigate the effect of a degraded GPS constellation. A more implementable algorithm to derive the UDRE was described in (Blanch et al., 2014) by taking into account the message broadcast interval and the lower limit on UDRE. The modifications and details of the improved service and integrity performance including the UDRE algorithm were recorded in (Walter et al., 2018) which provided a comprehensive understanding of the WAAS. A Dual Frequency Range Error (DFRE) estimation method based on the maximal projection of ephemeris-clock covariance matrix to form an envelope for the corrected error was introduced in (Shao et al., 2020). Finally, Rho and Richard (2007) and Heßelbarth and Wanninger (2013) examined the accuracy of WAAS orbits and clocks as well as investigated the possibility to use these corrections in the Precise Point Positioning (PPP). PPP implemented using GPS broadcast navigation messages and SBAS corrections in conjunction with integrity algorithms to produce protection levels were proposed in (Gunning et al., 2019).
The previous investigations focus on the generation of ephemeris corrections or integrity parameters, while this paper gives comprehensive description to calculate the SBAS message and performance analysis of the derived contents. Moreover, many studies focused mainly on the augmentation of GPS. With the development of satellite navigation systems, more constellations and signals are available to users, leading to smaller positioning errors and better availability. The SBAS Interoperability Working Group (IWG) has proposed a standard for next generation satellite based augmentation system, namely, Dual-Frequency Multi-Constellation (DFMC) SBAS (Zhao et al., 2020). All of 30 BeiDou-3 Navigation Satellite System (BDS-3) satellites have been successfully lunched as of June, 2020 and BDSBAS is also under developing to provide the navigation service for civil aviation users. It is meaningful to analyze the accuracy of satellite ephemeris corrections and the integrity performance of both GPS and BDS-3. In this research, a snapshot method based on the prior information of broadcast ephemeris is applied to estimate the satellite clock and orbit corrections, and the UDRE and MT28 derivation method based on the covariance matrix is investigated. The data collected at 26 stations in China are employed to perform the analysis.
This paper is comprised of 4 sections. Following the introduction, section "Snapshot ephemeris error estimation" derives the mathematical model to calculate the ephemeris corrections and integrity parameters. Section "Performance analysis" introduces the data and products used to carry out the research and discusses quantitatively the results. Conclusions are drawn in the last section.

Snapshot ephemeris error estimation
Snapshot algorithm uses a geometric approach to get the ephemeris errors, but it is sensitive to the satellite-station observation geometry and data loss. The minimum variance estimation based on the prior information of the state is employed in this analysis to overcome the under determined problem. The mathematical models to estimate the corrections and integrity parameters are presented in this section.

Calculation of ephemeris corrections
GNSS dual frequency ionosphere free observation can be expressed as: where ρ j i is the ionosphere free observation, ρ 0 is the geometric distance between satellite j and station i, G is the unit vector from station to satellite, (x sat , y sat , z sat ) is satellite position, (x, y, z) is station position, r sat is the vector from earth center to satellite, r is the position vector of station, x, y, z is satellite orbit error, dt i and dt j are station and satellite clock errors, respectively, c is the speed of light, d trop is the troposphere delay and corrected with Hopfield model (Hopfield, 1969), and ε is the combination of measurement noise, multipath and residual errors. Assuming only satellite and station clock errors are involved, and other errors neglected, then Eq. (1) can be rewritten as: A certain station is selected as reference to calculate the clock errors of satellites and other stations using the least squares estimation. It can be deduced that the satellite clock errors will absorb the radial errors of the satellite orbit. The resulted satellite clock errors are treated as the long-term clock corrections. Only orbit errors are contained in the observation residuals after removing the clock errors. The minimum variance estimation is applied to determine the orbit errors Tsai, 1999): X MV the estimated orbit error, and Λ is the a priori covariance matrix based on the uncertainty of broadcast ephemeris. The satellite orbit uncertainty here are 2.61 m, 5.45 m and 13.25 m and the satellite clock uncertainty is 2.61 m (Jefferson & Bar-Sever, 2000). Note that the above sigma values are defined in the radial, cross-track and along-track coordinate system, therefore, Λ is first transformed to the Earth Centered Earth Fixed (ECEF) coordinate system. W the measurement weight matrix, z is the range residual acquired with the initial values of the estimation parameters, and P MV is the posterior covariance. The orbit errors obtained in this step are served as long-term orbit corrections. The range residuals after removing the long-term corrections are used to estimate the clock errors again and will be taken as the fast corrections.

Integrity check of corrections
The range residuals obtained with both the fast and long-term corrections will be applied to derive the residual errors remaining in the satellite ephemeris with Eq. (5) (Wu & Peck, 2002): where e is the residual or remaining error in the corrected ephemeris, and P is the a priori variance of the corrected ephemeris. The difference between Λ in formula (4) and P here is that Λ represents the variance of original ephemeris and P is the variance of corrected ephemeris. P is the a posteriori variance of e , and z c is the measurement residual corrected by both the fast and long-term corrections which is also different from z in Eq. (4). Afterward, the resulted residual errors in ephemeris and their associated covariance are utilized to guarantee the integrity of corrections: k state is the chi-square statistic of residual errors remaining in the ephemeris, and k FA is a test threshold that sets the false alert rate. Theoretically, if the state e has a Gaussian distribution and has a covariance matrix specified by P , k state will follow the chi-square distribution, so the constant k FA can be determined. For example, a value of 4.3 for k FA will give a false alert probability of 1 × 10 -3 (Wu & Peck, 2002). The test is a simple comparison of k state with the constant k FA such that if k state < k FA then the satellite passes the test, otherwise the satellite fails the test and will be flagged as unusable.

Generation of UDRE and MT28
Integrity parameters UDRE and MT28 are broadcasted along with corrections to users for calculating the positioning protection levels. MT28 is used to adjust the UDRE to reconstruct location specific error bounds for users in different places (Rho & Richard, 2007;Walter et al., 2001). For each satellite, the SBAS receiver forms the standard deviation corresponding to the clock and orbit errors from the UDRE index and the 4 by 4 upper triangular matrix R included in MT28. The UDRE index provides the value of σ UDRE and R gives the value of the scaling matrix C MT28 . The uncertainty of the fast and long-term corrections is given by: σ flt is the standard deviation that characterizes the residual error associated with the fast and long-term corrections. σ flt is applied to weight the satellite observation and to compute the protection levels in SBAS positioning. σ UDRE is related to the UDRE index (RTCA, 2013) and is the minimum projected value of the ephemeris corrections covariance. The first three elements in vector u are related to the line-of-sight vector from a user in the service area to satellite, and the fourth is one. σ UDRE and C MT28 can be computed from the covariance matrix with formula (8) and (9): where P broadcast is the scaled covariance matrix and will be broadcasted to users through UDRE index and MT28. k md is determined by missed detection probability corresponding to Gaussian distribution and is taken as 6.13 here for a missed detection probability of 4.5 × 10 −10 . k FA is related to the false alert rate and is the same as in Eq. (6). P is the full covariance matrix calculated using the measurements at all stations viewing the satellite, and P is a function of satellite location, station observational geometry, and measurement confidence. The SBAS threat model should protect against a worse case of error measurement at any reference stations, so the P broadcast should meet the requirement: where e i is the residual error in ephemeris and P i is the subset covariance matrix constructed by removing the measurements from the i-th reference station in the least squares estimation process. Theoretically computing the best P broadcast is a complicated mathematical problem and different methods have been proposed (Blanch et al., 2014;Wu & Peck, 2002). The first method reported in (Wu & Peck, 2002) is used in this research and a brief introduction is given below. A theoretical UDRE of user j at the one sigma level can be defined as: where σ j is the standard deviation for user j in the service volume and is the reference value used to compute the broadcast covariance matrix. The line-of-sight vector for the j-th user augmented by one in the fourth component is represented by u j . Then a scaling factor F 0 can be derived with F 0 is the maximum value over all possible user j in the service area. Then the broadcast covariance matrix P broadcast can be determined with The UDRE and R matrix in MT28 are calculated as follows: where U is the upper triangular decomposition matrix of P broadcast and U 4,4 will be the UDRE.

Performance analysis
In order to verify the performance of the proposed algorithms, GPS and BDS-3 satellite corrections as well as UDRE and MT28 are estimated with the data collected at 26 stations on 1st to 7th July 2020. The station distribution is given in Fig. 1.
In section "Corrections with precise ephemeris", the corrections related to the IGS Multi-GNSS Experiment (MGEX) precise clock and orbit products are derived. With the high accuracy of IGS products, the satellite clock and orbit corrections should be very small and can be used to assess the precision of the results. The corrections along with the integrity results related to broadcast ephemeris are presented in section "Corrections and integrity results with broadcast ephemeris". Section "Positioning results" shows the positioning improvement with the employment of the corrections.

Corrections with precise ephemeris
IGS MGEX precise clock and orbit products are used to determine the corrections. Ionosphere free range residuals can be acquired with the original observations, station coordinates, precise satellite orbit and clock information.
A cutoff elevation angle of 15° is selected in the process. The errors induced by the IGS products can be neglected, for the associated corrections should be very small theoretically. Figure 2 illustrates the long-term clock corrections for GPS and BDS. Figure 2a shows the corrections for GPS PRN01 satellite while the lower panel Fig. 2b for BDS PRN19. The x-axis represents time in second of week, and y-axis the magnitude of corrections in nanosecond. The red dots are the clock corrections while the cyan dots are the average satellite elevation angle. The clock corrections are statistically (− 0.28 ± 0.44) ns and (− 0.26 ± 0.26) ns for GPS and BDS, respectively, which show a good consistency with the IGS products considering only the code observations are used. A slip in GPS corrections at 4.93 × 10 5 s is because of the change of the reference station clock. In the beginning and end of the satellite track arc, the accuracy of the results is lower because of the observations are more noisy and fewer stations have a common view. Figure 3 presents orbit corrections for the same satellite. The red dots in Fig. 3 are the calculated orbit errors in X component, the green and blue dots are the errors in Y and Z components, respectively. The results show that the average three-dimensional accuracy of orbit corrections is better than 0.4 m for both GPS PRN01 and BDS PRN19. The clock and orbit corrections obtained in this step will be used as the long-term corrections. In order to investigate the combined effect of the long-term corrections, clock error is recalculated and handled as the fast correction. The long-term corrections are generated every second in this research, so the error caused by broadcast delay is not introduced. The fast corrections of the two satellites are displayed in Fig. 4, with accuracy better than 0.15 ns.
Due to the usage of code observations, the clock and orbit corrections can be determined with single epoch data which makes the algorithm react rapidly to a satellite maneuver or clock fault.
Root Mean Square (RMS) of the long-term and fast corrections for all satellites with 7 days' data are plotted in Fig. 5. Figure 5a is correction statistics of GPS while Fig. 5b is of BDS satellites. The red bar in Fig. 5 is for long-term clock corrections, green bar is for orbit corrections, and blue one is for fast clock corrections. The unit for both clock and orbit corrections is meter. The average RMSs of different type corrections are shown in the figure. The upper panel shows the corrections of all GPS satellites and the lower part that of BDS satellites. After corrected the long-term errors, the remaining errors as the fast corrections are smaller than 0.1 m for GPS and BDS.

Corrections and integrity results with broadcast ephemeris
Broadcast ephemeris are used to obtain the corrections in this part. Figures 6 and 7 illustrate the long-term items of GPS PRN01 and BDS PRN19 satellites.
The long-term clock corrections of GPS in Fig. 6a and BDS in Fig. 6b are less than 2.0 ns. Orbit corrections in Fig. 7a, b for both systems are less than 2.0 m. The Because of no broadcast delay the fast term can be used to assess the accuracy of the long-term clock and orbit corrections. Figure 8 gives the RMS of the fast corrections for all satellites. The mean RMS is better than 0.30 ns. It can be concluded that the range error caused by the corrected ephemeris will be smaller than 0.10 m, which is sufficient to meet the system availability requirement.
Satellite UDRE is generated according to the method described in section "Generation of UDRE and MT28". Equation (6) is first applied to guarantee the integrity of the generated corrections and the threshold is set 4.3. The results of k state for all satellites observed at stations are shown in Fig. 9. The red line is the threshold, and the green dots are the statistics of the corrections. It can be seen that all satellites pass the test which means the preliminary integrity of the corrections.
As described in (RTCA, 2013), a total of 16 UDRE indices from 0 to 15 will be broadcasted to users. UDRE index (UDREI) of 15 indicates a detected fault or integrity risk of satellite, and index of 14 denotes that the satellite is not viewed by the monitoring stations or the satellite fails the test. The satellites with these two UDRE indices cannot be employed in the positioning solution. Figure 10 presents the UDRE of the GPS satellites viewed at the stations. The colorful lines in the figure are the ground tracks of satellites. Different colors represent different UDRE indices as indicated in the color bar. It can be found that UDRE at the start and end of the satellite track arc possess bigger indices because of the fewer viewed stations and lower accuracy observations. The maximum UDRE index in this analysis is 12 corresponding to a range accuracy of 15.2 m, and the minimum UDRE index of 3 implies a range uncertainty of 0.53 m. The minimum UDRE indices are over the China, which is reasonable due to the station locations. MT28 is also broadcasted to provide the relative covariance matrix for clock and ephemeris errors. It is an expansion of the UDRE which specifies the correction confidence as a function of user location (Walter et al., 2001). Figure 11 shows the projected covariance of a satellite, as an example, for all surface users with a viewing angle greater than 5°.
The satellite is located in 12°N and 78°E, represented as a red star in the figure. In this scene, the minimum UDRE is in China with a value of 1.15 m and the maximum UDRE occurs in the southwest direction with a value about 7.50 m. With the original covariance matrix displayed in Fig. 11, the scaling matrix C MT28 is derived using Eq. (14). The discretization error is not considered here. Figure 12 shows the projection of C MT28 onto the user space. Note it is similar in shape to Fig. 11. As expected, the minimum projected value of this normalized matrix is 1.0. With this message, users can reconstruct their location specific error bound rather than applying the largest bound in the service volume. Thus, both the availability within the service volume and the integrity in the region outside can be improved (Walter et al., 2001).
By the definition, UDRE should bound 99.9% of the combined fast and long-term correction errors (RTCA, 2013). The geometric ranges derived with IGS precise products and corrected ephemeris are used to evaluate the UDRE integrity. The difference between these two ranges is served as the remaining errors in the corrected ephemeris. Two issues need to be considered: one is the satellite antenna phase center correction when using the precise orbit product, and the other is the different time reference implied in precise and broadcast ephemeris clocks. Figure 13 plots the GPS range differences of all 7 days' data which exhibit a systematic bias of 0.6 m possibly due to the different time reference. The bias is neglected here so the result is somewhat conservative. Figure 14 presents the UDRE and associated range residuals of GPS and BDS satellites viewed by the station selected as time reference. The results show that the residuals can be bounded by UDRE with high confidence which verifies the validity of the UDRE generation method.

Positioning results
Three positioning tests are performed using the broadcast ephemeris and associated corrections generated in section "Corrections and integrity results with broadcast ephemeris". The first case is the positioning with only GPS ephemeris and corrections while the second case is the positioning with only the BDS data. GPS and BDS products are the combined together in the third case to evaluate the multi-constellation performance. The results are compared with those of the GPS standalone positioning. Global Ionosphere Map (GIM) generated by Center for Orbit Determination in Europe (CODE) was applied to correct the ionosphere delay for both positioning modes. Figure 15 gives the positioning errors of 26 stations. The x-axis shows the station number and y-axis the position error. The red bar represents the GPS standalone results, the green and blue ones are the results of BDS and GPS with the ephemeris and associated corrections and will be donated as BDS SBAS and GPS SBAS positioning. The cyan one is the results of multi constellation. The RMSs of GPS standalone, BDS SBAS and GPS SBAS position errors are 1.70, 1.08 and 0.98 m, respectively, while it is 0.80 m for the multi constellation case.  Page 12 of 14 Jin et al. Satell Navig (2021) 2:15 Conclusions SBAS is a civil aviation safety critical system that provides the correction and integrity information as a complementary data to the primary GNSS constellation for performance and reliability enhancement. The minimum variance estimation is applied to calculate the broadcast ephemeris corrections and a covariance based UDRE and MT28 derivation method is implemented in this research. GPS and BDS corrections together with integrity parameters are derived with real data. The results indicate that range accuracy of 0.10 m can be obtained with the application of the generated corrections. Satellite integrity parameters are calculated based on the a posteriori covariance. A clear dependence of UDRE on satellite elevation angle is observed. The relationship between the covariance and the scaling matrix C MT28 is presented to illustrate the advantage of employing MT28. The geometric range differences between the precise and broadcast ephemeris are utilized to verify the integrity of UDRE which shows the ability to bound the 99.9% correction errors. Performance improvements show the validity of the method proposed in this investigation which is helpful to the development of BDSBAS and Dual Frequency Multi Constellation (DFMC) SBAS.