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_{i}^{{{s}}}\)) and the received time at a receiver (\(t^{{{r}}}\)). The code pseudorange of the \(i\)th satellite (\(R_{i}\)) can be expressed as

$$R_{i} = \rho_{i} + c\left( {\Delta t^{{{r}}} - \Delta t_{i}^{{{s}}} } \right) + \delta_{{i,{{ion}}}} + \delta_{{i,{{tro}}}} + \varepsilon_{i}$$
(1)

where \(\rho_{i}\) is the geodetic range, \(\Delta t^{{{r}}}\) and \(\Delta t_{i}^{{{s}}}\) are the receiver and satellite clock offsets, c is the speed of light, \(\delta_{{i,{{ion}}}}\) is the ionospheric delay, \(\delta_{{i, {{tro}}}}\) is the tropospheric delay, and \(\varepsilon_{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.

Fig. 1
figure 1

Comparison between the geodetic range \({\rho }_{i}\) and the code pseudorange \({R}_{i}\)

Full size image

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:

  1. (a)

    The fast correction (pseudorange correction: \({{PRC}}\)) is obtained from the MT2 to 4 for each available satellite. The corrected pseudorange \(R_{{{{corr}}}} \left( t \right)\) at time \(t\) is computed by

    $$R_{{{{corr}}}} \left( t \right) = R_{{{{raw}}}} \left( t \right) + PRC\left( {t_{{{{of}}}}} \right) + \frac{{PRC\left({t_{{{{of}}}}} \right) - PRC\left( {t_{{{of,pre}}}} \right)}}{{t_{{{{of}}}} - t_{{{of,pre}}}}} \times \left({t - t_{{{{of}}}}} \right)$$
    (2)

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.

  1. (b)

    The long-term correction (interpolated satellite position error) is included in the MT25. The actual satellite position error vector for each satellite (\(\mathbf{\delta r}^{{{s}}}\)) can be computed from

    $$\mathbf{\delta r}^{{{s}}} \left( t \right) = \left[ {\begin{array}{*{20}c} {\delta x} \\ {\delta y} \\ {\delta z} \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} {\delta \dot{x}} \\ {\delta \dot{y}} \\ {\delta \dot{z}} \\ \end{array} } \right] \times \left( {t - t_{{{o}}} } \right)$$
    (3)

where \(\left[ {\delta x, \delta y, \delta z} \right]^{{{T}}}\) is the broadcast satellite position errors, \(\left[ {\delta \dot{x}, \delta \dot{y}, \delta \dot{z}} \right]^{{{T}}}\) is the rate-of-change of the satellite position errors, and \(t_{o}\) is the time-of-day applicability. Afterward, the corrected satellite position error vector (\({\mathbf{r}}_{{{{corr}}}}^{{{s}}}\)) is calculated by

$${\mathbf{r}}_{{{{corr}}}}^{{{s}}} \left( t \right) = {\mathbf{r}}_{{{{int}}}}^{{{s}}} \left( t \right) + {\mathbf{\delta r}}^{{{s}}} \left( t \right)$$
(4)

where \({\mathbf{r}}_{{{{int}}}}^{{{s}}} = \left[ {x, y, z} \right]^{{{T}}}\) is the interpolated satellite position coordinate vector based on the Earth-Center Earth-Fixed (ECEF) coordinate of the World Geodetic System 1984 (WGS-84) computed from the GPS navigation message. Additionally, the corrected clock offset for each satellite (\(\Delta t_{{{{corr}}}}^{{{s}}}\)) at the received time are obtained from MT25, i.e.,

$$\Delta t_{{{{corr}}}}^{{{s}}} \left( t \right) = \Delta t_{{{{int}}}}^{{{s}}} \left( t \right) + \delta a_{{{{f}}0}} \left( t \right) + \delta a_{{{{f}}1}} \times \left( {t - t_{{{o}}} } \right)$$
(5)

where \(\Delta t_{{{{int}}}}^{{{s}}}\) is the interpolated satellite clock offset from the GPS navigation messages, and \(\delta a_{{{{f}}0}}\) and \(\delta a_{{{{f}}1}}\) are the clock offset and clock drift error corrections, respectively.

  1. (c)

    The vertical ionospheric correction (VIC) of each ionospheric pierce point (IPP) based on the single-layer ionosphere model (Takasu 2013) is broadcast by MT26 and MT18. Then the slant ionospheric delays (\(\delta_{{{{ion}}}}\)) for each satellite are computed by

    $$\delta_{{{{ion}}}} = - \left( {\sqrt {1 - \left( {\frac{{R_{{{e}}} \cos \left( {El} \right)}}{{R_{{{e}}} + h_{{{I}}} }}} \right)^{2} } } \right)^{ - 1} \times \tau_{{{{vpp}}}} \left( {\varphi_{{{{pp}}}} ,\lambda_{{{{pp}}}} } \right)$$
    (6)

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 \(\tau_{{{{vpp}}}}\) is the VIC at the IPP interpolated based on three- or four-point interpolation algorithm by using the latitude (\(\varphi_{{{{pp}}}}\)) and longitude (\(\lambda_{{{{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.

Fig. 2
figure 2

Locations of some IGPs drawn on the geodetic map

Full size image

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

$${{HPL}} = K_{{{H}}} \times \sqrt {\frac{{d_{{{E}}}^{2} + d_{{{N}}}^{2} }}{2} + \sqrt {\left( {\frac{{d_{{{E}}}^{2} - d_{{{N}}}^{2} }}{2}} \right)^{2} + d_{{{{EN}}}}^{2} } }$$
(7)

and

$${{VPL}} = K_{{{V}}} \times d_{{{U}}}$$
(8)

where \(d_{{{E}}}^{2}\) and \(d_{{{N}}}^{2}\) are the variances of the true error distribution on the east and north axes, respectively, \(d_{{{{EN}}}}^{2}\) 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

$$\left[ {\begin{array}{*{20}c} {d_{{{E}}}^{2} } & {d_{{{{EN}}}} } & {d_{{{{EU}}}} } & {d_{{{{ET}}}} } \\ {d_{{{{EN}}}} } & {d_{{{N}}}^{2} } & {d_{{{{NU}}}} } & {d_{{{{NT}}}} } \\ {d_{{{{EU}}}} } & {d_{{{{NU}}}} } & {d_{{{U}}}^{2} } & {d_{{{{UT}}}} } \\ {d_{{{{ET}}}} } & {d_{{{{NT}}}} } & {d_{{{{UT}}}} } & {d_{{{T}}}^{2} } \\ \end{array} } \right] = \left( {{\mathbf{G}}^{{{T}}} \cdot {\mathbf{W}} \cdot {\mathbf{G}}} \right)^{ - 1}$$
(9)

where \({\mathbf{G}}\) is the geometry matrix, and \({\mathbf{W}}\) is the weighting matrix. The \(i\)th row of the matrix \({\mathbf{G}}\) is

$${\mathbf{G}}_{i} = \left[ {\begin{array}{*{20}c} { - \cos \left( {{{El}}_{i} } \right)\sin \left( {{{Az}}_{i} } \right)} \\ { - \cos \left( {{{El}}_{i} } \right)\cos \left( {{{Az}}_{i} } \right)} \\ { - \sin \left( {{{El}}_{i} } \right)} \\ 1 \\ \end{array} } \right]^{T}$$
(10)

where \({{Az}}\) is the azimuth angle. Meanwhile, the diagonal \({\mathbf{W}}\) matrix is the total error (variance) of the \(i\)th satellite (\(\sigma_{i}^{2}\)), i.e.,

$$\sigma_{i}^{2} = \sigma_{{i,{{flt}}}}^{2} + \sigma_{{i,{{UIRE}}}}^{2} + \sigma_{{i,{{tro}}}}^{2} + \sigma_{{i,{{air}}}}^{2}$$
(11)

where \(\sigma_{i,flt}^{2}\) 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, \(\sigma_{i,UIRE}^{2}\) is the variances of the user’s ionospheric range error (UIRE), \(\sigma_{i,tro}^{2}\) is the variance of tropospheric delays applied by the MOPS mode (RTCA 2006), and \(\sigma_{i,air}^{2}\) is the variance of airborne receiver errors assumed as the Class 2, 3, and 4 equipment. Additionally, \(\sigma_{UIRE}^{2}\) for each satellite is computed by

$$\sigma_{UIRE}^{2} = \left( {1 - \left( {\frac{{R_{e} \cos \left( {El} \right)}}{{R_{e} + h_{I} }}} \right)^{2} } \right)^{ - 1} \times \sigma_{{{{UIVE}}}}^{2}$$
(12)

where \(\sigma_{{{{UIVE}}}}^{2}\) 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 multi-frequency 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 (\(\rho_{i}\), \(\Delta t_{i}^{{{s}}}\), and \(\delta_{{i,{{tro}}}}\)) of (1) are obtained in three steps as follows. First, if the precise receiver position (\({\mathbf{r}}^{{{r}}} = \left[ {x,y,z} \right]^{T}\)) is known then the parameter \(\rho_{i}\) is computed based on the Pythagorean theory to calculate the geometric distance between the \(i\)th GPS satellite position corrected by GAGAN (\({\mathbf{r}}_{{i,{{corr}}}}^{{{s}}}\)) in (4) and the precise receiver position, i.e.,

$$\rho_{i} = ||{\mathbf{r}}_{{i,{{corr}}}}^{{{s}}} - {\mathbf{r}}^{{{r}}} ||= \sqrt {\left( {x_{i} - x} \right)^{2} + \left( {y_{i} - y} \right)^{2} + \left( {z_{i} - z} \right)^{2} }$$
(13)

Afterward, the \(\Delta t_{i}^{{{s}}}\) and \(\delta_{{i,{{tro}}}}\) parameters from the GPS navigation message are corrected by equation (5) and the MOPS model, respectively. Therefore, \(\rho_{i}\), \(\Delta t_{i}^{{{s}}}\), and \(\delta_{{i,{{tro}}}}\) can be canceled, i.e.,

$$\delta_{{i, {{t}},{{ion}}}} + c\Delta t_{{{t}}}^{{{r}}} + \varepsilon_{{i,{{t}}}} = R_{{i,{{t}}}} - \rho_{{i,{{t}}}} + c\Delta t_{{i,{{t}},{{corr}}}}^{{{s}}} - \delta_{{i,{{t}},{{tro}}}}$$
(14)

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

$$X_{i} \left( {t,k} \right) = \sqrt {1 - \left( {\frac{{R_{{{e}}} \cos \left( {El_{{i,{{t}}}} } \right)}}{{R_{{{e}}} + h_{I} }}} \right)^{2} } \times \left( {\delta_{{i, {{t}},{{ion}}}} + c\Delta t_{{{t}}}^{{{r}}} + \varepsilon_{{i,{{t}}}} - k} \right)$$
(15)

where \(k\) is a positive real number for estimating \(c\Delta t_{{{t}}}^{{{r}}} + \varepsilon_{{i,{{t}}}}\). Next, we propose the bias estimations (\(B\)) of \(c\Delta t_{{{t}}}^{{{r}}} + \varepsilon_{{i,{{t}}}}\) based on the minimum sum standard deviation technique, i.e.,

$$B = \mathop {\min }\limits_{0 \le k} \left( {\mathop \sum \limits_{i = 1}^{{{N^s}}} \sqrt {\mathop \sum \limits_{t = 0}^{{{T}}} \frac{{\left( {X_{i} \left( {t,k} \right) - {\overline{\mathbf{X}}}_{i,k} } \right)^{2} }}{T}} } \right)$$
(16)

where \(N^s\) is the number of GPS satellites, \(T\) is the period for estimating the parameter \(B\), and \({\mathbf{X}}_{i,k}\) is the matrix of \(X_{i} \left( {t,k} \right)\) for \(t = 0, \ldots , 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\Delta t_{{{t}}}^{{{r}}} + \varepsilon_{{i,{{t}}}}\) values (\(\hat{b}\)), i.e.,

$$\hat{b} \approx B + \min \left( {{\mathbf{X}}_{{i,{{B}}}} } \right)$$
(17)

where \({\mathbf{X}}_{{i,{{B}}}}\) is the matrix of \(X_{i} \left( {t,B} \right)\) for \(t = 0, \ldots , T\) and \(i = 1, \ldots , N^s\). Finally, the estimated VIC of the \(i\)th GPS satellite at time \(t\) (\({{VIC}}_{{i,{{t}}}}\)) can be computed by

$${{VIC}}_{{i,{{t}}}} = \sqrt {1 - \left( {\frac{{R_{{{e}}} \cos \left( {El_{{i,{{t}}}} } \right)}}{{R_{{{e}}} + h_{{{I}}} }}} \right)^{2} } \times \left( {\delta_{{i, {{t}},{{ion}}}} + c\Delta t_{{{t}}}^{{{r}}} + \varepsilon_{{i,{{t}}}} - \hat{b}} \right)$$
(18)

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 \(\sigma_{{{{UIVE}}}}^{2}\) values are computed by (19) in the next section. Importantly, a user then applies the \({{VIC}}\) and \(\sigma_{{{{UIVE}}}}^{2}\) 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.

Fig. 3
figure 3

Block diagram for estimating the VICs and the UIVE variances

Full size image

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. The variance of a user’s UIVE, \(\sigma_{{{{UIVE}}}}^{2}\), can be computed from

$$\sigma_{{{{UIVE}}}}^{2} = \left( {\sigma_{{{{VIC}}_{{{{Ref}}}} }}^{2} + \sigma_{{{{decorr}}}}^{2} } \right) \cdot R_{{{{irreg}}}}^{2}$$
(19)

where

$$\sigma_{{{{VIC}}_{{{{Ref}}}} }}^{2} = \left( {\left[ {\begin{array}{*{20}c} 1 \\ 0 \\ 0 \\ \end{array} } \right]^{T} \cdot \left( {{\mathbf{D}}^{{{T}}} {\mathbf{D}}} \right)^{ - 1} \cdot \left[ {\begin{array}{*{20}c} 1 \\ 0 \\ 0 \\ \end{array} } \right]} \right)$$
(20)

and

$${\mathbf{D}} = \left[ {\begin{array}{*{20}c} 1 & {d_{{{{IPP}}_{1} ,{{Ref}}}} \cdot \hat{E}} & {d_{{{{IPP}}_{1} ,{{Ref}}}} \cdot \hat{N}} \\ 1 & {d_{{{{IPP}}_{2} ,{{Ref}}}} \cdot \hat{E}} & {d_{{{{IPP}}_{2} ,{{Ref}}}} \cdot \hat{N}} \\ \vdots & \vdots & \vdots \\ 1 & {d_{{{{IPP}}_{N^s} ,{{Ref}}}} \cdot \hat{E}} & {d_{{{{IPP}}_{N^s} ,{{Ref}}}} \cdot \hat{N}} \\ \end{array} } \right]_{N^s \times 3}$$
(21)

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, \(\sigma_{{{{decorr}}}}^{2}\) and \(R_{{{{irreg}}}}^{2}\) 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).

Fig. 4
figure 4

Estimation of ionospheric vertical errors by using the positions of IPPs and the reference station

Full size image

In Fig. 4, \(d_{{{IPP,Ref}}}\) is the distance between the IPP and Ref points, and \(\hat{E}\) and \(\hat{N}\) are the unit vectors of the local coordinate system in the east and north axes, respectively, (Takasu 2013) as computed by

$$\left[ {\begin{array}{*{20}c} {\hat{E}} \\ {\hat{N}} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} { - \sin \left( {\lambda_{{{{Ref}}}} } \right)} & { - \sin \left( {\varphi_{{{{Ref}}}} } \right)\cos \left( {\lambda_{{{{Ref}}}} } \right)} \\ {\cos \left( {\lambda_{{{{Ref}}}} } \right)} & { - \sin \left( {\varphi_{{{{Ref}}}} } \right)\sin \left( {\lambda_{{{{Ref}}}} } \right)} \\ 0 & {\cos \left( {\varphi_{{{{Ref}}}} } \right)} \\ \end{array} } \right]^{T} \cdot \left( {\frac{{{\mathbf{r}}_{{{{corr}}}}^{{{s}}} - {\mathbf{r}}^{{{{Ref}}}} }}{||{{\mathbf{r}}_{{{{corr}}}}^{{{s}}} - {\mathbf{r}}^{{{{Ref}}}}|| }}} \right)$$
(22)

where \(\lambda_{{{{Ref}}}}\) and \(\varphi_{{{{Ref}}}}\) are, respectively, the longitude and latitude of the reference station, and \({\mathbf{r}}^{{{{Ref}}}} = \left[ {x,y,z} \right]^{{{T}}}\) is the reference position vector.

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 consist of [15° N, 100° E], [10° N, 100° E], and [15° N, 105° E], respectively. Likewise, the ionospheric correction of the reference station is estimated based on the proposed method. The user positions applying the ionospheric correction of the GAGAN (IGPs) or proposed (Local) are calculated based on the single-point positioning (SPP) algorithm.

Fig. 5
figure 5

Locations of user position (KMIT) and reference station (STFD)

Full size image

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 (MarchD) 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 \(\hat{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, post-midnight hours to pre-sunrise hours.

Table 1 Parameters in the experiment and data processing
Full size table

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 MarchD, 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.

Fig. 6
figure 6

ROTI plots at the STFD station during a MarchD, b March, c June, d September, and e December 2019

Full size image

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 (\(\delta_{{i,{{vion}}}}\)) from the difference between the VICs and remaining errors at the reference station (\(\delta t^{{{r}}}\)), i.e.,

$$\delta_{{i,{{vion}}}} \left( t \right) = VIC_{i} \left( t \right) - \sqrt {1 - \left( {\frac{{R_{{{e}}} \cos \left( {El_{{i,{{t}}}} } \right)}}{{R_{{{e}}} + h_{{{I}}} }}} \right)^{2} } \times \delta t^{{{r}}} \left( t \right)$$
(23)
Fig. 7
figure 7

Comparisons of the IGP-VIC at KMIT station and the local VICs at KMIT and STFD stations during a MarchD, b March, c June, d September, and e December 2019

Full size image

where the \(\delta 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 MarchD. 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 \(\hat{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.

Fig. 8
figure 8

Comparison between ionospheric delays (top) and user positioning errors in terms of HPE (middle) and VPE (bottom) in MarchD and March 2019 (equinox season, left panel) and in June 2019 (solstice season, right panel)

Full size image

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.

Fig. 9
figure 9

Positioning error comparison of the baseline GAGAN method between March and June 2019

Full size image
Fig. 10
figure 10

Positioning error comparison of the improved GAGAN method with proposed local ionospheric estimation between March and June 2019

Full size image

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.

Fig. 11
figure 11

Comparison of the 95% user positioning errors between the proposed and baseline methods in terms of HPE (top) and VPE (bottom) in 2019

Full size image

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 \(\sigma_{{i,{{UIRE}}}}^{2}\) 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 \(\sigma_{{{{UIVE}}}}^{2}\) values based on the baseline method on days: 2, 7, and 10.

Fig. 12
figure 12

Comparisons of HPL (top) and VPL (bottom) in March 2019

Full size image
Fig. 13
figure 13

Comparison of the protection levels between the proposed and baseline methods in terms of HPL (top) and VPL (bottom) in 2019

Full size image

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 LPV-200 (CAT-I equivalent) 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 \(\sigma_{{i,{{air}}}}^{2}\) and \(\sigma_{{i,{{UIRE}}}}^{2}\). The \(\sigma_{{i,{{air}}}}^{2}\) 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.

Fig. 14
figure 14

Vertical performances of APV-I (top), and LPV-200 (bottom) in March 2019

Full size image

The HMIs associated with this type of error could be removed by excluding satellites from position solution when the multipath errors (estimated by code-minus-carrier) is large. Based on equation (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.

Table 2 Percentages of availability at the KMIT station for APV-I and LPV-200
Full size table

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 \(\sigma_{{{{UIVE}}}}^{2}\) 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.