ARMA Prediction of SBAS Ephemeris and Clock Corrections for Low Earth Orbiting Satellites

For low earth orbit (LEO) satellite GPS receivers, space-based augmentation system (SBAS) ephemeris/clock corrections can be applied to improve positioning accuracy in real time.The SBAS correction is only available within its service area, and the prediction of the SBAS corrections during the outage period can extend the coverage area. Two time series forecasting models, autoregressive moving average (ARMA) and autoregressive (AR), are proposed to predict the corrections outside the service area. A simulated GPS satellite visibility condition is applied to the WAAS correction data, and the prediction accuracy degradation, along with the time, is investigated. Prediction results using the SBAS rate of change information are compared, and the ARMA method yields a better accuracy than the rate method.The error reductions of the ephemeris and clock by the ARMAmethod over the rate method are 37.8% and 38.5%, respectively.TheARmethod shows a slightly better orbit accuracy than the rate method, but its clock accuracy is even worse than the rate method. If the SBAS correction is sufficiently accurate comparing with the required ephemeris accuracy of a real-time navigation filter, then the predicted SBAS correction may improve orbit determination accuracy.


Introduction
A space-based augmentation system (SBAS) improves global navigation satellite system (GNSS) positioning accuracy by providing ephemeris/clock and ionospheric delay corrections.As with ground GNSS users, low earth orbit (LEO) satellites are able to use SBAS correction information if their GNSS receivers have the capability to receive SBAS signals.
The corrections are transmitted in real time along with GNSS signals.The SBAS correction transmission does not require extra data links and is suitable for space applications.
LEO satellites may use the SBAS corrections inside the system service area, but the corrections are not available outside the service area.The service area of an SBAS is determined by the geographical location of its monitoring stations.Due to this limitation, SBAS coverage is mainly limited to land area and does not cover most of the ocean area, which occupies 70% of the Earth.To maximize the usability of the SBAS for LEO satellites, it is necessary to extend the area where the SBAS corrections are available.Prediction of the SBAS corrections during the outage period can be a solution.
Among the SBAS corrections, the ionospheric correction is a position-specific value, and its prediction outside the service area is not useful.However, the ephemeris/clock correction is a GNSS satellite-specific value.Preceding researches discussed SBAS correction accuracy and its applications [1][2][3], but the prediction of the corrections has not been covered.A series of researches were performed for predicting differential GPS (DGPS) corrections [4,5], but the DGPS corrections are applied to the pseudoranges and their characteristics are entirely different from the SBAS corrections, which are applied to the three dimensional (3D) satellite coordinates.Other than the ephemeris correction, the author demonstrates that the SBAS ionosphere correction can be applied to the LEO GPS receiver with a proper scaling factor [6]. Prediction of the SBAS ephemeris/clock correction has not been attempted yet.If the SBAS ephemeris/clock correction is available for LEO satellites, the combined correction of SBAS ionosphere and ephemeris/clock can improve positioning accuracy of the LEO satellites.Recent advances in real-time LEO orbit determination filters enables submeter accuracy by fixing GPS carrier phase ambiguities [7].More accurate 2 International Journal of Aerospace Engineering ephemeris/clock information can improve the positioning accuracy and the filtering stability.A relative navigation between formation flying satellites can benefit from the SBAS corrections.
The SBAS message provides the rate of change data for the corrections to fill the time gap between data transmit intervals.The use of this data to extend the correction period is tested.However, the rate of change data was designed for a short interval only, and another method for long interval prediction is required.There are various time series forecasting models, for example, simple polynomial extrapolation, autoregressive moving average (ARMA), autoregressive (AR), and neural network.We propose two forecast models, an ARMA model and an AR model, to predict the ephemeris/clock corrections outside the service area.These models are selected because they are two of the most widely used models.Other prediction techniques have been tested but are not included in this paper: a polynomial model and a neural network model.The polynomial model yielded much lower prediction accuracy than the rate model and is dropped from consideration.A preliminary analysis by the neural network model showed low prediction accuracy and frequent divergence.Using the neural network requires more extensive work and is not included in this paper.The SBAS ephemeris/clock correction is predicted by using the ARMA and AR models by simulating the GPS satellite visibility condition.The ARMA/AR prediction accuracy is compared with the results using the SBAS rate of change data.

SBAS Corrections
SBAS messages include GPS satellite ephemeris/clock corrections, pseudorange fast corrections, and ionospheric delay corrections.The message is classified into message type (MT), and one MT is transmitted every second.The ephemeris/clock corrections are mainly included in MT 25.The corrections and their rate of change are also included in MT 25 [8,9].The rates of change are optional and their presence is identified by a velocity code in MT 25.In the case of WAAS, the rate of changes is included most of the time.The ephemeris corrections are expressed in Earth centered Earth fixed (ECEF) coordinates and three values (, , and ) are included.The correction is calculated from the following equation: where  is a transmitted correction of ECEF , , , or clock and  is the transmitted rate of change. is a current time to compute the correction. 0 is the time of applicability that is included in the SBAS messages.Each MT 25 includes the issue of data (IOD).This IOD has to correspond with the IOD ephemeris (IODE) of the current GPS ephemeris for that satellite.If the IODE of the GPS broadcast ephemeris data does not match the IOD of the SBAS transmission, it is an indication that the GPS ephemeris data sets have changed.The user must continue to use the old matched data from the previous broadcasts, until a new MT 25 with matching IOD is broadcast for the new GPS ephemeris data set [8].
There is a difference between the SBAS coverage and service areas.The SBAS service area represents an area where SBAS provides precise differential corrections, while the SBAS coverage area represents an area where the SBAS signal covers.The service area is limited by ground based infrastructure such as GNSS monitoring stations.Since the SBAS signal is transmitted from a geostationary satellite in equatorial orbit, the coverage area is much wider than the service area.In the case of WAAS, both of North and South America are in the coverage area, but only a part of North America is in the service area.
An important aspect of the SBAS is the integrity function, which is given by providing the error bound of the corrections.The SBAS is primarily developed for aviation use, and the integrity function allows an aircraft user to estimate an error bound and to be alarmed on fault signals.The accuracy of the correction information degrades over time, and the error bound grows over time.MT 10 provides a degradation factor to adjust the integrity values.The SBAS equipment should only use data if it is current; that is, before it has timed out.The time-out intervals are different for each message type and are a function of the aviation flight phase.In case of MT 25, the time-out interval is different for the type of aircraft landing approach, 360 s for a nonprecision approach and 240 s for a precision approach.

Prediction Methodology
The AR model predicts a future output as a combination of past inputs.A linear AR model can be represented as follows: where  is the sampling time.The AR process () is a stationary stochastic time series, and the input or error () is a white noise.Parameter  is the AR order.  is the AR coefficient that can be determined by various method, for example, least square method and Yule-Walker equations.
The ARMA model is a combination of the AR model and the moving average (MA) model.The MA model predicts a future output as a combination of past errors.ARMA is appropriate when a system is a function of a series of unobserved noise as well as its own behavior.The general ARMA model can be expressed as [10,11] where   is the MA coefficient.The coefficients   and   are determined by using the past data.Parameter  is the AR order and  is the MA order, and they should be finely tuned to optimize the performance.The orders are usually set as  =  − 1.
Four corrections, ECEF , , , and clock, are predicted separately by applying an independent ARMA/AR process.If there are strong correlations among the prediction components, a combined prediction, that is, multivariate ARMA/AR, could be more efficient.However, the correlations among the components are not significant, which will be discussed later, and the independent prediction process is used.The rate method computes a linear prediction using the last available slope (rate of change) data.The ARMA method determines the model coefficients   and   using available data and then predicts ().The AR method determines   only.
Comparing with the rate model, the ARMA model may cause a divergence in its prediction.Unstable model coefficients may cause intermittent divergences of the prediction value.This is not a distinctive feature of the ARMA/AR model; any other prediction model may cause a divergence.Although the divergence rarely happens for the SBAS prediction, approximately less than once per day for each PRN, a counter measure should be prepared.In order to prevent the divergence problem, a simple monitoring algorithm is developed.The output of the ARMA/AR method is constantly compared with the output of the rate method for determining the divergence condition.If the magnitude of the ARMA/AR prediction output is greater than three times of the rate output, then the last output before the divergence is used instead.The factor three was determined, as a rule of thumb, from analyzing the prediction results.

LEO GPS and SBAS Data Processing
A LEO satellite orbit is used to simulate the GPS and SBAS signal availability.The NASA/DLR Gravity Recovery and Climate Experiment (GRACE) satellite's orbit is selected [12].The GRACE mission is a gravity mapping mission using two low Earth satellites.The satellites, launched in 2002, are currently operational and are measuring static and timevariable Earth gravities.The initial altitude of the GRACE satellites was 480 km, and its inclination is 89.5 deg.Each GRACE satellite carries a high precision GPS receiver, and the GPS satellite visibility information is obtained from the actual GPS measurements.Any LEO satellite orbit can be used for the simulation, but the GRACE orbit is selected due to its orbit condition, that is, a near polar orbit covering most of the area.Another reason is its high precision GPS receiver.The GPS satellite visibility may be simulated by using the LEO and GPS satellite orbits only, but the actual measurements are used to simulate more realistic conditions.Between the two GRACE satellites' data, we used GRACE-B data because it has fewer outlier signals than GRACE-A data.Its sampling rate is 10 seconds.
The first step of the simulation is to determine the GPS satellite visibility condition for the LEO satellite from the actual GPS measurements.This process determines which GPS satellite signal is available for the LEO satellite at each epoch.The second step is matching SBAS corrections to the GPS satellite at each epoch.The third step is to determine the SBAS signal availability using the LEO satellite location and the SBAS coverage area.This SBAS information is the input of the prediction process and the SBAS information of the second step is the reference information to match the prediction output.The reference SBAS signal is generated by applying the rate of change of the corrections up to 360 s, which is the time-out interval of the nonprecision approach.The data interval is 10 s as the GRACE GPS measurements.
Among the SBAS, that is, WAAS, MSAS, EGNOS, and GAGAN, one SBAS is used for the analysis.Combined use of different SBAS corrections can extend the signal availability, but this option is not adopted for this analysis due to different levels of correction characteristics and accuracies.The correction accuracy varies significantly with the systems [13,14], and it is not possible to design a single optimal prediction model for multiple systems.WAAS is selected for the prediction analysis.This is because WAAS has many ground stations in a large area, and WAAS has the largest service area and has good correction accuracy.Another candidate is EGNOS, which has correction accuracy comparable to WAAS.Since the EGNOS coverage area is nearly on the opposite side of the WAAS coverage area, a mixed use of WAAS and EGNOS extends the coverage area, and it fulfills most of the flight times except in the south hemisphere.GAGAN can also be considered for its geographical location, close to the equator, and it may be helpful to enhance the prediction accuracy in the south hemisphere.After the optimization of the prediction algorithm for single SBAS, the prediction algorithm can be tuned for other SBAS.A combined prediction algorithm which consists of multiple system-specific algorithms can be used to extend the coverage.
Figure 1 represents a simulated SBAS message availability plot using GRACE satellite trajectory on May 1, 2014.At the time three WAAS satellites are working, PRN 133 at 98 ∘ W, PRN 135 at 133 ∘ W, and PRN 138 at 107 ∘ W [15].The three satellites transmit nearly identical corrections and could be treated as the same satellite.The signal coverage is the sum of the three satellite coverages.The mean of the GRACE orbit altitude was 432.0 km and the inclination was 89.0 ∘ .The color of the trajectory represents the number of available SBAS messages for visible GPS satellites.In North America the SBAS messages are available for most of visible GPS satellites.The number of visible satellites ranges from 4 to 10.In South America, the number is decreased because the WAAS messages are generated from the GPS satellites visible to the WAAS monitoring stations in North America.In Polar International Journal of Aerospace Engineering regions, the WAAS signal is not available, but propagation using the rate of change information up to 360 s enables the use of the WAAS corrections.In North Pole, the availability is different for ascending and descending paths in North America.During the ascending path, moving south to north, the WAAS correction is propagated.During the descending path, no WAAS signals are available.Figure 2 represents the number of visible GPS satellites and the number of WAAS-available GPS satellites per each epoch.The time represents the second of the day on May 1, 2014.The number of SBAS corrections is presented as the color code in Figure 1.When GRACE satellite is within the WAAS service area, that is, from 31000 s to 32000 s, the number of GPS signals is equal to the number of SBAS signals.
As GRACE moves out of the service area but within the coverage area, from 32000 s to 33000 s, the number of SBAS signals drops rapidly.During this period, many GPS satellites are visible to GRACE but not to WAAS monitoring stations.
Figure 3 shows the time series of the SBAS corrections on May 1, 2014, for PRN 22.The clock correction unit is converted from seconds to meters by multiplying by the speed of light.The correction data is not available when the satellite is not in the WAAS service area.The time series variation does not show a simple pattern and it may not be possible to model the time series with a simple function, for example, polynomials.The magnitude of the ephemeris corrections, norm of , , and  corrections, has a mean of 3.183 m for all PRNs on this date.The RMS of the ephemeris and clock corrections are 4.276 m and 2.436 m, respectively.The data interval, when the SBAS signal is available for data acquisition, ranges from 370 s to 2690 s, while the prediction interval (outage of available corrections) ranges from 10 s to 3270 s.The mean of the data acquisition and prediction intervals is 1868 s and 1914 s, respectively.

Prediction Results
SBAS ephemeris/clock corrections are predicted with the ARMA and AR models.To determine an optimal ARMA order, a series of tests has been performed by changing the order.The orbit error differences among the different orders are less than 5 cm, far below the SBAS accuracy, and higher orders (>10) caused intermittent divergence.Differences among the , , , and clock are not significant.Since a higher ARMA order yields a lower error in general, the ARMA order is selected as 7, a stable highest order.With the same procedure, the AR order is selected as 5.The AR model causes more divergence than the ARMA model, and its order is set to be lower than the ARMA.A prediction using the rate of change is performed for comparison.Figure 4 compares the predicted corrections by the ARMA, AR, and the rate of change methods on May 1, 2014.ECEF Z components of PRN 22 from 30000 s to 42000 s are presented where the time interval is the same as Figure 2.During the first prediction interval, correction data from 31320 s to 33170 s are used to build the ARMA and AR models, and the data from 33180 s to 36200 s are predicted.The IOD was changed on 36200 s, and the prediction is performed until the same IOD as the data acquisition interval is used.The rate method caused a linear deviation, the slope of which is the last available rate of change during the data interval.During the first prediction interval, the ARMA and AR methods return the true value after a short time interval, but the rate method yields a large deviation.During the second data acquisition interval (38800 s∼41710 s), the last rate of change of the previous data acquisition interval is zero, and the rate method yields a constant value prediction and its accuracy is better than the ARMA for 650 s.After this short interval, the true value is changed and the ARMA method follows this change while the rate method remains as a constant.The AR method shows a mixture of the two methods, similar to the ARMA during the first prediction interval but similar to the rate during the second interval.
Figures 5, 6, and 7 represent the prediction errors of PRN 22 on May 1, 2014.The correction errors grow over the prediction time, and the ARMA and AR methods show a substantially smaller error than the rate method.In all methods, the maximum error occurs in  corrections, and the maximum values are 4.56 m, 2.25 m, and 3.86 m for the rate, ARMA, and AR methods, respectively.The difference in the maximum value is reduced for the clock error: 1.09 m, 0.75 m, and 0.99 m for the rate, ARMA, and AR methods, respectively.Correlations among the , , , and clock errors are not significant.The constant level of the clock errors, shown at approximately 40000 s and 63000 s, correspond to the constant level of the clock corrections in Figure 3.The constant corrections are caused by the zero rate of change information.
Theoretically, the residuals from the ARMA or AR should be close to white noise if the model fits the time series perfectly.The fitting residuals from the ARMA and AR are analyzed in the frequency domain, but they are not shown as a figure.All the residual power spectral density (PSD), that is, , , , or clock component, shows a similar linear decrease in the log-log plot with a slope of −0.485log()/log(Hz).In the frequency range from 0.03 Hz to 0.05 Hz, the PSD is  close to white noise.Existence of the linear slope implies that the prediction/fitting model can be improved by removing a long term trend.A preliminary experiment was performed by using an autoregressive integrated moving average (ARIMA) model that is useful to reduce the long term effect.However, no improvement was made.The SBAS correction is a piecewise signal separated by IOD changes, and the fitting interval is relatively short when compared with the number of data points.For this reason, detrending the long term effect requires further research.
Figures 8 and 9 represent the prediction RMS errors for all 30 PRNs on May 1, 2014.The 3D orbit correction error is represented as a norm of , , and  RMS errors, and the clock error is computed separately.Mean of the 3D orbit RMS errors for all PRNs is 1.049 m, 1.670, and 1.744 m for the ARMA, AR, and rate methods, respectively; the reduction ratio by the ARMA method over the rate method is 40.0%.In Figure 3, the mean RMS of the SBAS ephemeris  correction magnitude was 4.276 m.The ARMA RMS error of 1.744 m is substantially lower than the signal magnitude, which demonstrates the usefulness of the prediction results.An exception is PRN 26, where the rate error is smaller than the ARMA error.However, both of the two models show a high error, and it is not meaningful to compare the error level for PRN 26.The standard deviation among the PRNs is 0.822 m, 1.237 m, and 1.099 m for the ARMA, AR, and rate methods, respectively.The error magnitude variation among the PRNs is similar to all three methods because the error level mainly depends on common factors: length of prediction intervals, length of data acquisition period, and so forth.The mean of the clock RMS errors for all PRNs is 0.610 m, 1.321 m, and 0.977 m for the ARMA, AR, and rate methods, respectively.The reduction ratio by the ARMA is 37.6% and close to the orbit reduction ratio.
The prediction accuracy depends on the length of the prediction time.Figure 10 shows the orbit and clock predicted accuracies for different time intervals on May 1, 2014.The axis value 600 s represents the prediction time from 10 s to 600 s, and the value 1200 s represents the time from 610 s to 1200 s, and so forth.The total number of prediction data is 60688 and the data percentage ranges from 29.4% (0 s∼600 s) to 0.6% (3010 s∼3600 s).The rate method shows a linear increase of the error in both orbit and clock.The ARMA method shows a linear increase in orbit but the slope is lower than the rate method, and the difference grows as the time increases.The difference between the rate and ARMA methods is magnified in the clock error.The clock error by the ARMA method is significantly lower.At 2410 s∼3000 s interval, the RMS error reductions by the ARMA method over the rate method are 46.0%(orbit) and 42.2% (clock).
The AR method yields a lower or similar accuracy than the rate method in short intervals up to 1200 s.The clock accuracy of the AR method is even worse than the rate method.The stochastic nature of the clock time series does not have a significant correlation over time, and the AR is not appropriate for the clock prediction.As mentioned earlier, the RMS of the correction magnitude is 4.276 m for the orbit and 2.436 m for the clock on this date.If we set a maximum allowable prediction error less than a half of the correction magnitude, they become approximately 2 m for the orbit and 1 m for the clock.With this justification, the ARMA prediction up to 1800 s can be feasible for use.If a single SBAS is used for the corrections, the 1800 s is not sufficient for covering most of the flight time.If multiple SBAS corrections, for example, EGNOS, MSAS, and GAGAN, are used, the 1800 s is sufficient for most of the flight time except in the South Pole region.
The prediction accuracy may depend on the correction data characteristics, and it is reasonable to analyze other days of data.Five days from 2014 are chosen to evaluate the accuracy, from January 1 to September 1 with two month interval.The latest 2014 dates are selected because the WAAS correction data accuracy and characteristics changed with the improvement of the system [16].Figure 11 shows the prediction RMS error on five different days in 2014.The prediction errors in all time intervals from 10 s to 3600 s are computed.At all dates, the ARMA method results in a better accuracy than the rate method.The mean of the daily orbit RMS errors is 1.12 m (ARMA) and 1.77 m (rate) with an average of a 37.8% reduction.The mean of the daily clock RMS errors is 0.64 m (ARMA) and 1.04 m (rate), and the reduction percentage is 38.5%, close to the orbit results.For the prediction interval from 1210 s to 1800 s, the five day mean of RMS errors are 1.164 m for the orbit and 0.614 m for the clock.The AR method shows a slightly better orbit accuracy than the rate method, but its clock accuracy is even lower than the rate method.
The impact of the SBAS correction on the user position accuracy depends on the various effects, for example, navigation filters and other errors.In general, a high precision filter using carrier phase measurements can benefit from the SBAS [3].Another aspect is the SBAS integrity function.In order to evaluate SBAS accuracy, proper integrity information should be used to determine the use of the correction information.Regardless of these facts, it is clear that improved ephemeris information is helpful to improve the user position.Evaluation of SBAS accuracy is a very complicated process due to the integrity.At this time the integrity is suited for aircraft use, and it should not be applied to spacecraft.Without applying the integrity information, it is not appropriate to compare the positioning accuracy.For this reason the positioning accuracy is not covered in this paper.

Conclusions
The position accuracy of a GPS receiver is limited by the broadcast ephemeris/clock errors.For LEO satellite GPS receivers, SBAS ephemeris/clock corrections can be an alternative solution for real-time use.LEO satellites move at a high speed and the duration for LEO satellites to be inside the SBAS service area is relatively short.Prediction of the SBAS corrections during the outage period can extend the service area.Two widely used forecasting models, ARMA and AR models, are proposed to predict the ephemeris/clock corrections outside the coverage area.The simulated GPS satellite visibility is applied to the WAAS correction data, and the prediction accuracy degradation along with the time is investigated.A comparison with the prediction results using the SBAS rate of change information is performed.The ARMA method shows a greater accuracy than the rate method.The error reductions of the ephemeris and clock by the ARMA method over the rate method are 37.8% and 38.5%, respectively.For the prediction interval from 20 min.to 30 min., the RMS errors are 1.164 m for the orbit and 0.614 m for the clock.The AR method shows a slightly better orbit accuracy than the rate method, but its clock accuracy is even worse than the rate method.The positioning accuracy improvement by the raw or predicted SBAS correction is not covered in this paper.Impact of the SBAS corrections on real-time orbit determination mainly depends on the SBAS correction accuracy level.If the SBAS correction is sufficiently accurate comparing with the required ephemeris accuracy of a real-time navigation filter, then the raw or predicted SBAS correction may improve orbit determination accuracy.

Figure 2 :
Figure 2: Number of total and WAAS-available GPS satellites observed at GRACE satellite orbit (May 1, 2014).

Figure 10 :
Figure 10: RMS of orbit and clock prediction errors for different time intervals (May 1, 2014).

Figure 11 :
Figure 11: Orbit and clock prediction errors for different days in 2014.