Performance improvement of the GAGAN satellite-based augmentation system based on local ionospheric delay estimation in Thailand

Satellite-Based Augmentation System (SBAS) is essential to support aircraft navigation. L1 SBAS operates on the L1 frequency (1575.42 MHz) and is currently still of interest since all GNSS satellites and receivers do not fully support additional frequencies such as L5 (1176.45 MHz). Although the Global Positioning System (GPS) aided Geo Augmented Navigation (GAGAN) SBAS is available, the performances are degraded due to the discrepancies of the ionospheric correction over Thailand and surrounding areas. Hence, in this work, we propose a new method based on the geometry-free ionospheric delay estimation with a single frequency (L1) and a single reference station requirement. The local ionospheric delays are estimated based on the proposed method with the observed GPS and GAGAN data in Thailand. Then the ionospheric corrections are obtained from the estimated local ionospheric delays. The analysis shows that using the estimated corrections, the positioning errors are reduced both on quiet days and locally disturbed days in 2019. More reductions in the positioning errors are found in September and December than other months. In addition, we perform a preliminary availability assessment of two critical phases of flights. The GAGAN performances with the proposed method for the APV-I and LPV-200 categories are improved up to 57% and 53%, respectively, in comparison with the baseline method of the IGP correction.


Introduction
Satellite-Based Augmentation System (SBAS) is crucial to improve services in aeronautical navigation based on Global Navigation Satellite System (GNSS). The augmented information for the SBAS-equipped GNSS receivers is broadcast by the SBAS geostationary satellites (European GSA 2020). Currently, there are two generations of SBAS: (1) the L1 SBAS operating on L1 frequency (1575.42 MHz) alone (RTCA 2006) and (2) the dual-frequency multi-constellation (DFMC) SBAS operating on the L1 and L5 (1176.45 MHz) frequencies (ICAO 2018b). However, since DFMC SBAS is still under the standardization process by the International Civil Aviation Organization (ICAO) and is not expected to be available until November 2022, the L1 SBAS system is currently still of much interest.
In the literature, most research works on the L1 SBAS focus on the assessment of SBAS availabilities and the improvement of augmented parameters. For example, in Bang et al. (2016), the Korea Augmentation Satellite System (KASS) was explored on the Korean peninsula and confirmed based on the two phase-of-flight categories. The LPV-200 and APV-I services (RTCA 2006) have availabilities of about 92% and 70%, respectively. Next, the GPS aided Geo Augmented Navigation (GAGAN) SBAS system over Sri Lanka was assessed by Dammalage et al. (2017), and it supports the APV-I service. Although SBAS services do not support the APV-II and CAT-I categories (the precision approach (PA) model based on the ICAO standard), some augmented parameters have been investigated for improvements. For instance, the improvements of the fast and long-term corrections for the Wide Area Augmentation System (WAAS) SBAS service were proposed by Chen et al. (2017) and Zheng et al. (2019), respectively. Moreover, Weng and Chen (2019) investigated the vertical integration of the SBAS service and a local monitoring station to reduce vertical protection levels (VPL). Dautermann et al. (2020) proposed to rebroadcast the corrected parameters of the SBAS service via the Ground-Based Augmentation System (GBAS) to the aircraft. Interestingly, based on the assumption that the DFMC SBAS system may not be operating soon, the modification of broadcast correction data frames is proposed by Reid et al. (2013). The L1 SBAS service from GAGAN and Multi-functional Satellite Augmentation System (MSAS) is available in Thailand and the surrounding areas. However, the ionospheric grid points (IGP) from GAGAN cover only some parts of Thailand (Pungpet et al. 2018) but still exceed those from MSAS.
The ionosphere extends from about 60 km to more than 2000 km above the earth. It is a partially ionized medium due to the X and UV rays of solar radiation and the incidence of charged particles. Importantly, the local ionospheric delays at the low latitude region are affected by the equatorial ionization anomaly (EIA) and equatorial plasma bubbles (EPB), (Balan et al. 2018;Timoçin et al. 2020); hence, the actual delays differ from the Klobuchar model (Mallika et al. 2020;Klobuchar 1987). To approximate the actual delays and the total electron content (TEC), most works focused on the estimation based on dual-frequency linear combinations. The minimum sum standard deviation of the TEC technique was proposed by Otsuka et al. (2002) and Ma et al. (2003) to remove the residual delays (combined ambiguity and receiver bias) within the estimated ionospheric delays. Sparks et al. (2011a, b) proposed the kriging method with GPS receiver networks in North America to upgrade the grid ionospheric vertical delay (GIVD) of WAAS. Meanwhile, regional ionospheric delays were mapped and predicted by a spherical cap harmonic function with past observations of the GPS data in China (Liu et al. 2011). Moreover, ionospheric irregularity threats that may evade detection by GNSS are constructed based on an undersampled ionospheric threat model (Bang et al. 2013) for Korean SBAS by using the GNSS station networks. Alternatively, the ionospheric estimations based on a single-frequency GNSS signal were proposed for users. The local ionospheric delays were estimated by the code-plus-carrier observation, which required the estimation of ambiguity parameters (Schüler et al. 2011) or known precise receiver location (Hein et al. 2016).
However, in single-frequency devices such as L1 SBAS receivers, dual-frequency data are not available.
Motivated by the above reasons, we propose a new local ionospheric delay estimation method based on the geometry-free ionospheric delay estimation (Krypiak-Gregorczyk and Wielgosz 2018) by using the observed data of GPS and GAGAN systems in Thailand. Afterward, we assess the GAGAN service quality for quiet and disturbed periods with and without the new method. Following, we give details of the fundamental theory and the proposed method. Then the methodology and experimental setup are described, followed by the results and discussions. Finally, the conclusions of this work are provided.

Pseudorange measurement
In the GNSS system, the code pseudoranges are computed from the difference between the transmitted time from the i th satellite ( t s i ) and the received time at a receiver ( t r ). The code pseudorange of the i th satellite ( R i ) can be expressed as where i is the geodetic range, Δt r and Δt s i are the receiver and satellite clock offsets, c is the speed of light, i,ion is the ionospheric delay, i,tro is the tropospheric delay, and i is the measurement error. The geometric range (dashed line) and code pseudorange (blue line) are depicted in Fig. 1. The delays and the range errors in (1) are normally removed using prediction models or the error corrections of information processed by reference stations on the earth surface. (1)

GAGAN corrections
In L1 SBAS by GAGAN, the fast, long-term, and ionospheric corrections are broadcast to the users in the service areas, including Thailand (Sophan et al. 2020) via geostationary SBAS satellites. The information is composed of several message types (MT) as explained in RTCA (2006). The details of these corrections are reviewed as follows: (a) The fast correction (pseudorange correction: PRC ) is obtained from the MT2 to 4 for each available satellite. The corrected pseudorange R corr (t) at time t is computed by where R raw is the raw pseudoranges, and t of and t of ,pre are the time of applicability of the most recent and previous fast correction, respectively. where Δt s int is the interpolated satellite clock offset from the GPS navigation messages, and a f 0 and a f 1 are the clock offset and clock drift error corrections, respectively.
ionosphere model (Takasu 2013) is broadcast by MT26 and MT18. Then the slant ionospheric delays ( ion ) for each satellite are computed by where R e is the radius of the earth, El is the elevation angle between a receiver and a satellite, h I is the height between the earth surface and IPP, and vpp is the VIC at the IPP interpolated based on three-or four-point interpolation algorithm by using the latitude ( pp ) and longitude ( pp ) at the IPP, respectively. At present, a few IGPs of GAGAN (Sunda et al. 2015) cover a partial area of Thailand illustrated in Fig. 2.

SBAS integrity
The integrity parameters indicate the ability of a system, and they are used to compute the integrity bounds, the so-called protection levels (PL) RTCA (2006). Generally, the horizontal protection level (HPL) and vertical protection level (VPL) for the PA model of the aircraft are computed, respectively, by and where d 2 E and d 2 N are the variances of the true error distribution on the east and north axes, respectively, d 2 EN is the covariance of the true error distribution on the east-north axis, d U is the standard deviation of the true error distribution on the vertical axis, K H is the constant value of HPL for LNAV/VNAV, LP, and LPV, and K V is the constant value of VPL. Additionally, the model distribution matrix of the true error distribution is where is the geometry matrix, and is the weighting matrix. The i th row of the matrix is where Az is the azimuth angle. Meanwhile, the diagonal matrix is the total error (variance) of the i th satellite ( 2 i ), i.e., where 2 i,flt is the residual errors of the fast and long-term corrections characterized by the user differential range error (UDRE) parameters from the MT2 to 4, MT7, MT10, and MT28, 2 i,UIRE is the variances of the user's ionospheric range error (UIRE), 2 i,tro is the variance of tropospheric delays applied by the MOPS mode (RTCA 2006), and 2 i,air is the variance of airborne receiver errors assumed as the Class 2, 3, and 4 equipment. Additionally, 2 UIRE for each satellite is computed by where 2 UIVE is the variance of the grid ionospheric vertical error (GIVE) in the MT26 related to each IGP shown in Fig. 2.

Proposed local ionospheric delay estimation method
As the IGP locations in the GAGAN system do not fully cover Thailand, therefore, the IGP edge may not represent the actual ionospheric delays. Hence, we propose a new method based on the geometry-free ionospheric delay estimation with a single frequency (L1) signal and a single reference station requirement. The advantage of the single-frequency method is that there is no need for differential code bias (DCB) parameters and avoiding the errors existing at lower frequency (such as scintillation in L5) and the multifrequency method (e.g., L1-L2 or L1-L5). Moreover, if the precise receiver position is known, the L1 GPS receivers can be employed to estimate the observed ionospheric delays. The right-hand side terms ( i , Δt s i , and i,tro ) of (1) are obtained in three steps as follows. First, if the precise receiver position ( r = x, y, z T ) is known then the parameter i is computed based on the Pythagorean theory to calculate the geometric distance between the i th GPS satellite position corrected by GAGAN ( s i,corr ) in (4) and the precise receiver position, i.e., Afterward, the Δt s i and i,tro parameters from the GPS navigation message are corrected by equation (5) and the MOPS model, respectively. Therefore, i , Δt s i , and i,tro can be canceled, i.e., Moreover, the slant-factor effects for each elevation angle for the i th GPS satellite at time t ( El i,t ) are corrected based on the single-layer ionosphere model (Takasu 2013). The residual error function ( X i ) can be expressed as where k is a positive real number for estimating cΔt r t + i,t . Next, we propose the bias estimations ( B ) of cΔt r t + i,t based on the minimum sum standard deviation technique, i.e., where N s is the number of GPS satellites, T is the period for estimating the parameter B, and i,k is the matrix of X i (t, k) for t = 0, … , T . However, the underestimated ionospheric delays in period T can occur due to measurement errors. Thus, we use the zero-TEC approach to obtain the estimated receiver bias cΔt r Figure 3 shows the overall process of estimating the local ionospheric correction and the variances of UIVE using the L1 GPS data at the reference station. The procedures are described by four steps as follows. In the process, the observed GPS and GAGAN data such as the code and carrier pseudorange, the ephemeris, and the fast and long-term corrections are fed input. First, the noise effects in the code pseudoranges are reduced based on the carrier smoothed pseudorange (RTCA 2006) to achieve the smoothed pseudoranges. Second, the smoothed pseudorange and satellite position errors are corrected based on the fast and long-term corrections of GAGAN according to equation (2) and equation (4-5), respectively. Next, the VICs of each satellite are estimated by equation (18). Finally, the 2 UIVE values are computed by (19) in the next section. Importantly, a user then applies the VIC and 2 UIVE parameters based on the satellite-by-satellite basis and uses the fast and long-term corrections of GAGAN based on the SBAS standard (RTCA 2006) to obtain the estimated position.

UIVE variance estimation
Generally, the planar fit technique with several reference stations is applied to compute the errors at each IGP. The ionospheric vertical error (UIVE) is often used to indicate positioning errors related to the ionospheric correction.
and are the formal error on the estimations and the observed matrix of i th IPPs based on the planar fit technique as explained in Fig. 4, 2 decorr and R 2 irreg represent the decorrelation function and the inflation factor defined based on the previous work (Walter et al. 2001), respectively. It should be noted that in (20), the W matrix is ignored because the new IGP point is still not created and not interpolated to each IPP. Additionally, the W matrix created from the measurement variances of the vertical ionospheric delays requires multiple reference stations (Walter et al. 2001). One suggestion is to construct the W matrix by logging the local vertical ionospheric delays based on the proposed method with a specific window size (for example, 5 min).
In Fig. 4, d IPP,Ref is the distance between the IPP and Ref points, and Ê and N are the unit vectors of the local coordinate system in the east and north axes, respectively, (Takasu 2013) as computed by (20)

Experimental setup
In this work, the KMIT and STFD stations are the user position and reference station, respectively. Their locations illustrated in Fig. 5 are around 12 km apart. The ionospheric correction of GAGAN is computed based on the three-point interpolation algorithm (RTCA 2006). Three IGPs in the latitude and longitude directions as shown in Fig. 2

Dataset
The parameters in the experiment and data processing are set as shown in Table 1. The precise locations of KMIT and STFD stations are obtained from the AUSPOS Online GPS Processing Service on the website https:// gnss. ga. gov. au/ auspos using the GPS data on March 17, 2019. Five-disturbed days in March (March D ) and five-quiet days in each month for March, September (equinox season), and June and December (solstice season) in 2019 are analyzed to show the improvements. Additionally, the selected period for estimating b is from 7:00 to 10:00 local time (LTC) because the fluctuation of the ionosphere is relatively lower than at other times and, importantly, EPB events are not observed during this period. The reason we do not choose the nighttime is that the EPBs can occur, particularly at post-sunset, postmidnight hours to pre-sunrise hours.
In this work, we identify the locally disturbed or quiet days, particularly, due to EPB events using the rate of TEC index (ROTI) in the unit of TECU/min (Schaer 1999), which is obtained based on the geodetic-free combination technique by using the carrier phase pseudoranges of L1 and L2 GPS data (of the STFD station). In Fig. 6, the ROTI plots are in the universal time (UTC) on the dates considered in Table 1. The time window of 5 min and the elevation masks of 30° are used, whereby distinct colors represent the PRNs of GPS satellites. In Fig. 6a, for March D , the ROTIs after sunset are above 0.5 TECU/min, but no geomagnetic storm occurs as confirmed by the Dst index (Nose et al. 2015). Therefore, the local disturbances  are possibly due to EPB. In Fig. 6b-e, for locally quiet days, the ROTIs are relatively lower than those in Fig. 6a.

Results and discussions
The performances of GAGAN are evaluated based on the proposed ionospheric correction. We specifically consider two phase-of-flights: the APV-I and CAT-I categories.

Comparisons between the IGP and local ionospheric corrections
In Fig. 7, the VICs of three 'IGPs' interpolated to the KMIT station are compared with the 'local' VICs estimated based on the proposed method using the observed GPS and GAGAN data of the KMIT and STFD stations. Note that the VICs of the proposed method still include the remaining errors (ionospheric delays + remaining fast and long-term errors), we can compute the estimated actual ionospheric delays ( i,vion ) from the difference between the VICs and remaining errors at the reference station ( t r ), i.e., where the t r term can be viewed as the 'receiver bias,' which can be estimated based on the SPP algorithm using the iterative least-square estimations (Takasu 2013) of the observed GPS and GAGAN data at the reference station. In addition, the KMIT (IGPs) results are computed using the three-point interpolation algorithm. The KMIT (Local), KMIT* (Local), and STFD* (Local) results are calculated based on the medians of all the visible satellites with the elevation mask of 30° and the moving average of a 4-min interval, without and with (*) the remaining errors, respectively. Figure 7a, b and d shows the results of the fifteen days in March and September, respectively. From the figures, the ionospheric delays at KMIT (IGPs) increase during daytime but decrease during nighttime exhibiting the typical shape of ionospheric delays. Similarly, the KMIT (Local) follows a similar trend to the KMIT(IGPs) except for some results after sunset on the 2nd and 5th of March D . On the other hand, the VIC at both KMIT* (Local) and STFD* (Local) are markedly different from the KMIT(IGPs) and KMIT (Local) since both are not the actual ionospheric delays, but they include the remaining errors in fast and long-term corrections (Chen et al. 2020). These discrepancies are common because the error corrections are different at each location, and Thailand is at the edge of the GAGAN service. Additionally, the overestimated or underestimated ionospheric values are dependent on the selected period of estimated b .

User positioning errors
Here, the horizontal positioning errors (HPE) and vertical positioning errors (VPE) are evaluated based on the SPP algorithm at the KMIT station illustrated in Fig. 8 (middle and bottom), respectively. The 'baseline' plots represent the performances of the GAGAN system, whereas the 'proposed' method refers to the local ionospheric correction of the STFD station. Figure 8 (top) shows the comparison of the ionospheric delays between the IGPs (GAGAN) and the proposed method in March 2019 (disturbed and quiet days, left panel) and June 2019 (quiet days only, right panel), respectively.
From the HPE and VPE results, the proposed method outperforms the GAGAN service. Higher VPEs appear in March than June 2019. Significantly, most positioning error improvements are made based on the ionospheric fluctuations, meaning that the proposed method can reduce the positioning errors always. Additionally, the occasional spikes are caused based on the multipath effects of the low elevation angle (5°) and the poor distribution of satellites (high dilution of precision). Alternatively, the positioning errors are compared based on the histogram plot illustrated in Fig. 9. The HPEs and VPEs tend to have longer tails in March than June so the equinox season (March) shows more errors than the solstice season (June). Likewise, the results of the proposed method in Fig. 10 show comparable positioning errors in either month. Moreover, HPEs and VPEs in Fig. 10 are lower than in Fig. 9. Figure 11 shows the 95% user positioning errors in 2019. The proposed method can decrease position errors during disturbed and quiet days. Significantly, the HPE and VPE results can reach below 0.5 and 1 m, respectively. In practice, the user can be far from the reference station. Thus, the interpolated approach is needed to obtain the more precise ionospheric correction.

Protection levels
Generally, for the aircraft positioning and navigation, at each altitude, when protection levels are above the alert limit (AL), the system will warn the users. Here, the HPL and VPL values are computed based on (7) and (8), where 2 i,UIRE is calculated based on (19) for the 'proposed' method. The K H and K V are respectively 6 and 5.33 (SBAS standard). Figure 12 shows the HPL and VPL results in March 2019 and indicates the 'proposed' method generally gives lower errors than the 'baseline' method. Moreover, these results are presented based on the 95% errors illustrated in Fig. 13. The proposed method reduces the errors of HPL by 52.5 and 57.6 m in September and December, respectively. On the other hand, error reductions of VPL by 87.2 and 98.0 m are observed in September and December, respectively. The best improvements in December are caused by the high levels and fluctuations of 2 UIVE values based on the baseline method on days: 2, 7, and 10.

Performance evaluations
Finally, we evaluate the performances of the improved GAGAN using the proposed local ionospheric estimation in March 2019 for the PA phase of the ICAO standard. The Stanford chart of each case on the five days (17, 19, 22, 26, and 27) of data in March are produced based on the vertical direction. The vertical alert limits (VAL) for the APV-I and CAT-I equivalent (LPV-200) in the SBAS system are 50 and 35 m (ICAO 2018a), respectively. The chart consists of four conditions: available, misleading information (MI), hazardous misleading information (HMI), and unavailable conditions. In Fig. 14, we plot the vertical performances of APV-I and  in March 2019 at the KMIT station, then show that the 'available' percentages are 92% and 66%, respectively. The HMI events in Fig. 14 are due to two factors: overestimation of VPE and underestimation of VPL. Overestimation of VPE is due to the overestimation of local VIC (Vertical ionospheric correction which include remaining errors). It is caused by the cycle slip occurrences at low elevation angles below 15°. The VIC is dependent on 2 i,air and 2 i,UIRE . The 2 i,air term is based on the modeled multipath errors in the airborne environment. Because the user station is on the ground, there may be cases with larger multipath errors than in the airborne environment.
The HMIs associated with this type of error could be removed by excluding satellites from position solution when  (20), the W matrix for local ionosphere is not used here. As previously described, it is constructed by the local vertical ionospheric delays in a specific window size. Although not shown in paper, with the elevation angles of 15° (rather than 5°) and the W matrix inclusion, the HMI events on the Stanford chart are reduced to zero.
In Table 2, we compare the availabilities of the PA model between the baseline and proposed methods in 2019. Then the available percentages of APV-I and LPV-200 are improved by 42% and 47% in September 2019, respectively. Additionally, the available percentages in December 2019 are also improved by 57% and 53%, respectively.

Conclusions
In this work, we propose a new local ionospheric delay estimation method by using the single-frequency L1 GPS data with a reference station requirement to improve the performance of GAGAN system in Thailand. The user positioning errors during locally disturbed and quiet days in 2019 are improved using the estimated local ionospheric corrections. The HPE and VPE are reduced to below 0.5 and 1 m, respectively. Additionally, the available GAGAN service is assessed based on APV-I and LPV-200 categories by the proposed method showing improvements of about 57% and 53%, respectively.
As the proposed method will aid the locations at the edge of a service area, a few solutions are possible (upon future discussions), for example, to coordinate with the SBAS ground control and include VIC and 2 UIVE into the MT18 and MT26 blocks, particularly for the missing grids, or to transmit these values from a local airport or local ground station to the approaching aircraft through designated wireless communication systems.

Data availability
The datasets analyzed during the current study are available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.