INSIGHT INTO SPATIAL ANALYSIS OF GPS BROADCAST ORBIT ERROR TOWARDS ORBITAL IMPROVEMENT

Broadcast orbits are compared against final orbit to get the error of broadcast orbit. The errors are analysed by presenting the error over space, especially longitude. The satellite trajectory is divided into three sector namely northern, southern, and transitional sectors. Spatial analysis show that the error is correlated with the latitude and longitude. Some consistency pattern can be observed from the distribution of the error in the spatial analysis. Standard deviation (SD) is used to quantify the consistency, providing more quantitative insights into the spatial analysis. Four patterns can be observed in the error distribution, namely consistency in northern and southern sector, consistency of transitional sector, changes after transitional sector, and correlation between ∆X component and ∆Y component. The spatial analysis shows potential to be used in broadcast orbit error estimation and prediction. A model that uses this predicted broadcast orbit error as a correction will be designed in the future to improve the broadcast orbit accuracy.


INTRODUCTION
Global Navigation Satellite System (GNSS) are satellite-based positioning systems that have played playing major roles in positioning, navigation, and timing. Advances in GNSS technology have diversify its use from daily activities such as ehailing, food delivery and navigation to various fields such as surveying, agriculture, maritime and aviation (Balafoutis et al., 2017;Kealy & Moore, 2017;Kin-Chung, 2020;Rizos, 2017). It has also seen its share of application in scientific fields such as meteorology, geodynamics and space weather (Langely et al., 2017). Since the full deployment of United States' Navstar Global Positioning System (GPS) in 1993, there are many GNSS that have been developed and providing active service including Russia's Global Navigation Satellite System (GLONASS), European Union Galileo, and China's BeiDou Navigation Satellite System (BDS). GPS is the earliest and one of the most commonly used GNSS in the world (Hein, 2020).
GPS positioning works using trilateration technique. The position of the user is determined from signals originating from a minimum of four satellites (Hofmann-Wellenhof et al., 2012). The position and clock bias of the satellite are used to determine the geometric range between the receiver and the satellite. A rootmean-squared (RMS) of 2 m error can be induced in the GPS measurement due to the orbital error (Karaim et al., 2018). Therefore, it is important to know the accurate position of the satellites during the emission of the GPS signals. This information is recorded and transmitted in ephemerides.
Currently, there are many types of ephemerides available, namely broadcast, ultra-rapid, rapid and final ephemerides (Johnston et al., 2017). Broadcast ephemerides are provided in real-time with an RMS accuracy of ±100 cm. On the other hand, final ephemerides provide the orbit position of highest accuracy at RMS of ±2.5 cm but with a latency of 14 to 18 days.
Final ephemerides are often taken as the 'true' position of the GPS satellites (Kim & Kim, 2015;Tusat & Ozyuksel, 2018). Error of the broadcast orbit can be computed by comparing it against the final orbit. Broadcast orbit error has been shown to portray variation over time (Jia et al., 2014;Kim & Kim, 2015) and satellite (Tusat & Ozyuksel, 2018). Some research has also considered into the signal in space range error (SISRE) (Montenbruck et al., 2014) and reference frame (Nicolini & Caporali, 2018) using broadcast orbit error.
This paper aims to analyse the broadcast orbit error from a spatial perspective. Analysing the spatial trend of the broadcast orbit error would provide some insight in modelling the broadcast orbit error as a correction to improve the accuracy of broadcast orbit.
The paper is divided into four sections. First section briefly introduces GPS and the research made on broadcast orbit error. The second section presents the methodology of the research, covering broadcast orbit error computation and spatial analysis. The third section presents the result and analysis of spatial analysis. The paper is then concluded with future recommendation in the last section.

METHODOLOGY
This section will present the methodology of broadcast orbit error computation and spatial analysis. The dataset sampled are GPS orbit from 3 rd January 2021 to 16 th January 2021, for a total of 14 days. All 32 GPS satellites were operational during this period.

Broadcast orbit error computation
Before computing the broadcast orbit error, it is important to extract the satellite position from the broadcast and final ephemerides. The satellite position in final ephemerides is given in an Earth-Centred, Earth-Fixed (ECEF) coordinate system. The final orbit is published in ITRF2014 geodetic datum with 15 minutes interval. Final orbit is also referenced to the (CoM) of the satellite (Kouba & Héroux, 2001).
On the other hand, broadcast orbit is also published in Keplerian elements. These elements are updated every 2 hours. Users will need to compute the satellite position in ECEF or other relevant coordinate system. Broadcast orbit is published in WGS84 geodetic datum. The coordinates are referenced to the APC of the satellite (Dunn & DISL, 2012). These differences must be considered during the computation of broadcast orbit error. It is necessary to put the broadcast and final orbit in a same coordinate system before comparing them. In addition, antenna phase centre offset will need to be applied to the broadcast orbit for it to be referenced to the CoM (Montenbruck et al., 2014). In terms of geodetic datum, WGS84 is consistent with ITRF2014, with a difference of ±10 cm (Qinsy, 2020). This difference will be absorbed as part of the measurement error in the broadcast orbit error.
The steps to compute broadcast orbit can be found directly in the Interface Specification 200 (IS-GPS-200) (Dunn & DISL, 2012). Borre (2003) also published a set of MATLAB codes to compute the satellite position. Table 2 shows the Keplerian elements distributed in the broadcast ephemerides.

Parameter Symbol
Amplitude of cosine harmonic correction term to the argument of latitude  (Pestana, 2015) Derivation of equation 1 to 16 uses the symbols in Table 2.
Computation of the broadcast orbit is done for each satellite on a per epoch basis. The computation starts with computing the semimajor axis (a) of the reference epoch.
The eccentric anomaly (E) is computed from the mean anomaly (M) using Newton's iteration. Initial approximation of eccentric anomaly takes mean anomaly as the eccentric anomaly.
Computation of the true anomaly (v) of the GPS satellite is performed followed by the unperturbed argument of latitude (u ̅).
The unperturbed argument of latitude will then be used to compute the periodic corrections of the radius, argument of latitude and inclination. These periodic corrections will be applied to the parameters to get the perturbed values.
= cos( ) (11) = sin( ) The corrected longitude of ascending node ( ) is required to convert the satellite position from the orbital plane to an ECEF coordinate system with the rotation of the earth (̇) at the constant value of 7.2921151467x10 7 rad/sec.
Lastly, the satellite position is outputted in ECEF coordinate system.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVI-4/W3-2021 Joint International Conference Geospatial Asia-Europe 2021 and GeoAdvances 2021, 5-6 October 2021, online Antenna phase centre correction must be applied to translate from APC (XYZAPC,BRDC) to CoM (XYZCoM,BRDC). The value can be taken from IGS14.atx in the antenna exchange (ANTEX) format.
Lastly, broadcast orbit error (∆X, ∆Y, ∆Z) is computed by comparing the final orbit with the CoM broadcast orbit. The comparison is made every 15 minutes, in consistent with final orbit's interval, in all three components. The broadcast orbit error for a given satellite is described as a function of epoch. [

Spatial analysis
While broadcast orbit error is commonly presented as a time series, it could also be presented over spatial changes. GPS satellites has a medium Earth orbit (MEO) with a period of 23 hours 56 minutes. It will visit the same location approximately once per sidereal day. Figure 1 shows the 3D orbit trajectory of a satellite with pseudo-random number (PRN) 12.

Figure 1. 3D orbit trajectory of PRN12
Spatial analysis of broadcast orbit error analyses the changes of the broadcast orbit error over space, especially longitude. The latitude and longitude are calculated from the coordinates of the broadcast orbit using the ellipsoid model for WGS84. The spatial information is then rearranged to match the broadcast orbit error of the corresponding epoch. Figure 2 shows the 2D satellite trajectory of PRN12.

Figure 2. 2D orbit trajectory of PRN12
The trajectory of the satellite will be analysed along with the broadcast orbit error by dividing them into three categories according to their latitude, namely northern sector (above 0.7 rad), southern sector (below -0.7 rad), and transitional sector (between 0.7 rad to -0.7 rad).
Consistency of the spatial trend in each individual sector will also be analysed using standard deviation (SD) that can be calculated using the equation below (Streiner, 1996):

SPATIAL ANALYSIS
This section presents the spatial analysis of broadcast orbit error. Figure 3 shows the spatial analysis of PRN12. Blue, green, and red are used to represent northern, southern, and transitional sector respectively. The three rows in Figure 3 show the broadcast orbit error over longitude of ∆X(n), ∆Y(n), and ∆Z(n) respectively. The value below each sector is the SD of the error in that sector and is calculated using equation 18.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVI-4/W3-2021 Joint International Conference Geospatial Asia-Europe 2021 and GeoAdvances 2021, 5-6 October 2021, online Consistency of the broadcast orbit error can be categorised into stable or dispersive. A stable consistency is defined as having an error distribution which has small variation and is more constant.
On the other hand, the error distribution of dispersive consistency has a larger variation and is more spread out. Consistency can be compared relatively in terms of SD.
Four observation was made from the spatial analysis and its SD analysis, namely (i) consistency in northern and southern sector, (ii) consistency of transitional sector, (iii) changes after transitional sector, and (v) correlation between ∆X and ∆Y.
The consistency of the ∆X and ∆Y in the northern and southern sector interchange between stable and dispersive. Take the ∆Y in Figure 3 as example, the error in -3.14 rad to -2 rad is relatively dispersive (SD: ±0.705 m) while it was relatively stable (SD: ±0.489 m) from -1.8 rad to -0.4 rad. The consistency of the ∆Z is mostly stable during the northern or southern sector, with SD less than ±0.409 m.
The consistency of the transitional sector for ∆X and ∆Y are dependent on the consistency of the sector that follows. When the following sector is stable, the consistency of that transitional sector is relatively stable. For example, the first transitional sector of ∆X is relatively dispersive (SD: ±0.597 m) since it is followed by a southern sector which is dispersive (SD: ±0.831 m). On the other hand, the second transitional sector of ∆X is relatively stable (SD: ±0.314 m), as the following sector was also relatively stable (SD: ±0.416 m). In the case of ∆Z, the transitional sector is mostly dispersive, having an SD of more ±0.685 m to ±1.020 m.
In addition, a change in the error distribution can be observed after every transitional trajectory. Taking the same example of ∆Y in Figure 3, the properties of consistency changed from dispersive to stable or vice versa after transitional trajectory.
Another interesting point to note is the relationship between the consistency properties of ∆X and ∆Y. The consistency shows opposite properties between ∆X and ∆Y. For example, ∆X was relatively stable (SD: ±0.449 m) while ∆Y was relatively dispersive (SD: ±0.705 m) during -3.14 rad to -2 rad. On the other hand, when the satellite entered -1.8 rad to -0.4 rad, ∆X became relatively dispersive (SD: ±0.831 m) while ∆Y became relatively stable (SD: ±0.489 m).
Similar properties were also observed in all available 32 satellites. In general, the errors can be dispersive or stable in the spatial analysis of all satellites. Some satellite displayed dispersive or stable consistency properties only in its trajectory. The first and second observation is only a general categorisation of the error. The third observation on the changes after every transitional trajectory is consistent in all satellites, with only some differences in the offset value. The fourth observation on the correlation on ∆X and ∆Y is also valid in all satellites, given that there were dispersive to stable changes in the error on the satellites.

CONCLUSION AND RECOMMENDATION
Spatial analysis and its SD shown some interesting information whereby stable or dispersive consistency trend can be observed in the northern and southern sector. The consistency of transitional sector also showed dependency on the consistency of northern or southern sector, especially for the ∆X and ∆Y components. The transitional sector on ∆Z displayed dispersive trend without dependency on the following sector. There are also changes in the characteristics and magnitude of broadcast orbit error after transitional sector. Lastly, the error in the ∆X and ∆Y component opposes each other where when one is dispersive, the other would most likely be stable in the same sector.
Analysing the broadcast orbit error trend from spatial perspective shed some new light in designing a model for broadcast orbit improvement. Models that could predict the broadcast orbit error from a spatial perspective can be designed. The predicted errors can then be applied to real-time broadcast orbit to minimise the orbital error. At this point, the spatial analysis is still insufficient for broadcast orbit error prediction. An in-depth investigation on the trend of the spatial variation still needs to be conducted and this will be covered in future research.