Introduction

Satellite-based augmentation system (SBAS) is a well-known augmentation technique in aeronautical navigation in many parts of the world (European GNSS Agency 2018; Li et al. 2020). It is particularly a vital part of enroute and possibly approach phases. The SBAS components consist of the reference network, master station, uplink station, geostationary SBAS satellite and users. The master station generates the SBAS corrections and integrity by using the observed data from the reference network, and then sends them to the SBAS satellite via the uplink station. The benefits of the SBAS technology include the corrections of the user positioning errors, the guarantee of the user’s availability, and the regional service. SBAS has been developed since 2003, when the global positioning system (GPS) satellites with L1 frequency are deployed such as the Wide Area Augmentation System (WAAS), ICAO SARPS (2006) and RTCA DO-229D (2006). Recently, the developments of the global navigation satellite system (GNSS) have matured, the dual-frequency multi-constellation (DFMC) SBAS based on multi-constellation satellites receive more attention, for example, DFMC SBAS based on GPS (ICAO 2018; Reid et al. 2013a, b), BeiDou SBAS (BDSBAS) and Quasi-Zenith Satellite System (QZSS) used to augment GPS, GLONASS, and Galileo constellations (Zhao et al. 2020; Takahashi et al. 2022). The L1 SBAS system broadcasts the correction and integrity information via the L1 frequency (1575.42 MHz) including the fast and long-term (pseudorange and orbit) corrections and the ionospheric grid points (IGP). Recently, the new DFMC SBAS standard using both L1 and L5 (1176.45 MHz) frequencies has been proposed; the advantage lies on the fact that it cancels the ionospheric effects, broadcasts only the satellite clock-ephemeris corrections and integrity information, and avoids the high dilution-of-precision (DOP) effects.

In DFMC SBAS, the satellite clock-ephemeris corrections can be generated by the residual errors in the pseudorange measurements (Shao et al. 2020). Typically, the ionosphere-free (IF) carrier smoothing code technique with the raw data of the reference stations is utilized to estimate the residual errors of the pseudorange measurements between a satellite and a receiver. The locally dispersed reference stations with the triple difference (TD) technique, Sickle and Dutton (2020), are used to approximate the residual satellite-clock offset errors. Although the IF combination is used to remove the first-order ionospheric delays (Hoque and Jakowski 2007), the higher-order ionospheric (HOI) delays still remain in the pseudorange residual errors (Subirana et al. 2013). Exactly, the HOI delays can increase the residual errors up to 10 cm (Liu et al. 2016). Therefore, the positioning accuracy and integrity in DFMC SBAS can be degraded by the satellite clock-ephemeris corrections with the HOI delays as well as long baseline (Sophan et al. 2022).

The HOI delays in various systems have recently been investigated (Qi et al. 2021; Guo et al. 2023; Xi et al. 2021; Yehun et al. 2020). A typical approach is to estimate the HOI delays from the total electron content (TEC) parameters (Marques et al. 2011; Hadas et al. 2017) based on the geometry-free (GF) combination (Krypiak-Gregorczyk and Wielgosz 2018). The HOI delays with the Gravity Recovery and Climate Experiment Follow-On (GRACE-FO) observation are estimated based on the satellite-borne GPS technique to improve the precise orbit determination of GRACE-FO (Qi et al. 2021). The maximum precision improvement is observed at the millimeter levels. The HOI delays with the low-earth orbit (LEO) satellite observation are approximated by applying the International Reference Ionosphere-2016 (IRI-2016) and the International Geomagnetic Reference Field: the 13th generation (IGRF-13) models (Guo et al. 2023). Moreover, the HOI delays with the Beidou Navigation Satellite (BDS) are analyzed by applying the TEC components of the IGRF-13 model and the Global Ionospheric Maps (GIM) model with the code pseudorange (Xi et al. 2021). The HOI effects can cause the user positioning errors in the ranges from 5 to 10 mm (Yehun et al. 2020).

To our knowledge, the local DFMC SBAS corrections with the HOI estimations or mitigations have not yet been considered in the previous research works. Therefore, in this work, we propose the procedures of local DFMC SBAS corrections together with the HOI delay mitigation. The HOI delays are computed from the local TEC values at stations in Thailand based on the Klobuchar model (Setti et al. 2019). In particular, the local DFMC SBAS corrections are generated by using the pseudorange residual errors of local reference stations with the minimum sum-variance and IF carrier smoothing code techniques. The system performances are evaluated on both quiet and disturbed days. Moreover, the availabilities of DFMC SBAS are evaluated for the Localizer Performance with Vertical guidance (LPV-200) and Category I precision approach (CAT-I) categories of the precise approach (PA) model (ICAO SARPS 2006).

In this article, the local DFMC SBAS corrections with the HOI delay mitigation have been derived in the background theory section. The methodology shows the simulation steps and data selections. In the simulation results, the DFMC SBAS corrections, user positioning errors, protection levels, and performances are illustrated in the results and discussions Section. Then, the conclusions are given at the end.

Background theories

User position estimations

The user coordinates and the receiver bias are typically calculated by using the single point positioning (SPP) algorithm based on the iterative weighted least-square estimation method (Takasu 2013). If the estimated user coordinate (\({\widehat{\mathbf{r}}}_{u}={\left[{\widehat{x}}_{u}, {\widehat{y}}_{u}, {\widehat{z}}_{u}\right]}^{T}\)) and the receiver bias (\({\widehat{b}}_{u}\)) are denoted as the column vector (\(\widehat{\mathbf{Z}}={\left[{\widehat{\mathbf{r}}}_{u}, {\widehat{b}}_{u}\right]}^{T}\)), the user position estimations of each time can be calculated by

$${\widehat{\mathbf{Z}}}_{k+1}={\widehat{\mathbf{Z}}}_{k}+{\left({\mathbf{H}}^{T}\mathbf{W}\mathbf{H}\right)}^{-1}{\mathbf{H}}^{T}\mathbf{W}\left(\mathbf{Y}-h\left({\widehat{\mathbf{Z}}}_{k}\right)\right)$$
(1)

where k is the iterative index, and h(·) is the measurement vector function of errors, the coordinate unit matrix (\(\mathbf{H}\)) is

$$\mathbf{H}={\left[\begin{array}{cccc}-\frac{{{\varvec{\uprho}}}_{u}^{1}}{\Vert {{\varvec{\uprho}}}_{u}^{1}\Vert }& -\frac{{{\varvec{\uprho}}}_{u}^{2}}{\Vert {{\varvec{\uprho}}}_{u}^{2}\Vert }& \cdots & -\frac{{{\varvec{\uprho}}}_{u}^{{N}_{s}}}{\Vert {{\varvec{\uprho}}}_{u}^{{N}_{s}}\Vert }\\ 1& 1& \cdots & 1\end{array}\right]}_{\left(4\times {N}_{s}\right)}$$
(2)

where \({{\varvec{\uprho}}}_{u}^{i}\) is the column vectors of the equivalent geometric ranges between the satellite \(i\in \left\{1, 2, \dots , {N}_{s}\right\}\) and the user based on the earth centered earth fixed (ECEF) coordinates system of the world geodetic system 1984 (WGS-84), i.e.,

$${{\varvec{\uprho}}}_{u}^{i}={\left({\widehat{\mathbf{r}}}_{u}^{i}+{{\varvec{\updelta}}}^{i}\right)}^{T}-{\widehat{\mathbf{r}}}_{u}$$
(3)

\({\widehat{\mathbf{r}}}_{u}^{i}\) is the coordinate vectors (\({\left[{x}_{u}^{i}, {y}_{u}^{i}, {z}_{u}^{i}\right]}^{T}\)) of the \(i\) th satellite estimated by the pseudorange at the user location, \({{\varvec{\updelta}}}^{i}\) is the satellite position error correction vectors (\({\left[{\delta x}^{i}, {\delta y}^{i}, {\delta z}^{i}\right]}^{T}\)) of the DFMC SBAS corrections, and \({N}_{s}\) is the total number of available satellites.

The weighting matrix (\(\mathbf{W}\)) is a square matrix. The diagonal matrix equals the inversion of the total errors for all available satellites as written by

$$\mathbf{W}={\text{diag}}\left({\sigma }_{1}^{-2}, {\sigma }_{2}^{-2}, \dots , {\sigma }_{{N}_{s}}^{-2}\right)$$
(4)

where \({\sigma }_{i}^{2}\), \(i\in \left\{1, 2, \dots , { N}_{s}\right\},\) is the total variances of the available number of satellites \({N}_{s}\) with the user at each time expressed as

$${\sigma }_{i}^{2}={\sigma }_{i,DFRE}^{2}+{\sigma }_{i,UIRE}^{2}+{\sigma }_{i,tro}^{2}+{\sigma }_{i,air\_DF}^{2}$$
(5)

\({\sigma }_{i,DFRE}^{2}\) is the dual frequency range error (DFRE) variance of the \(i\) th satellite, \({\sigma }_{i,UIRE}^{2}\) is the ionospheric-free residual variance (user ionospheric range error: UIRE based on ICAO SARPs (2018) of the \(i\) th satellite with the user, \({\sigma }_{i,tro}^{2}\) is the variance of the tropospheric delay applied based on the minimum operational performance specification (MOPS), RTCA DO-229D (2006), \({\sigma }_{i,air\_DF}^{2}\) is the variance of the airborne receiver errors for the dual frequency (L1-L5) modified from the single frequency of MOPS (Jie et al. 2018).

In addition, the observed vectors (\(\mathbf{Y}\)) of all available satellites with the user equal

$$\mathbf{Y}={\left[\begin{array}{c}{R}_{u}^{1}+{c\widehat{\Delta t}}_{u}^{1}-{\widehat{\delta T}}_{u}^{1}-{\widehat{\delta \Delta I}}_{u}^{1}\\ {R}_{u}^{2}+{c\widehat{\Delta t}}_{u}^{2}-{\widehat{\delta T}}_{u}^{2}-{\widehat{\delta \Delta I}}_{u}^{2}\\ \vdots \\ {R}_{u}^{{N}_{s}}+{c\widehat{\Delta t}}_{u}^{{N}_{s}}-{\widehat{\delta T}}_{u}^{{N}_{s}}-{\widehat{\delta \Delta I}}_{u}^{{N}_{s}}\end{array}\right]}_{\left({N}_{s}\times 1\right)}$$
(6)

where \({R}_{u}^{i}\) is the smoothed pseudorange measurements (in meter) between the u-th user and the i-th satellites, where \(i\in \left\{1, 2, \dots ,{N}_{s}\right\}\) based on ICAO SARPs (ICAO 2018), \(c\) is the speed of light in the vacuum (in meters/sec), \({\widehat{\Delta t}}_{u}^{i}\) is the \(i\) th satellite-clock offset (in sec) estimated by the code pseudorange at the user, \({\widehat{\delta T}}_{u}^{i}\) is the tropospheric delay (in meter) of the \(i\) th satellite estimated based on MOPS, RTCA DO-229D (2006), with the user, and \({\widehat{\delta \Delta I}}_{u}^{i}\) is the HOI delays (in meter) for the \(i\) th satellite with the user estimated from the local TECs based on the Klobuchar model (Setti et al. 2019). The HOI delays are derived by Marques et al. (2011), Hadas et al. (2017), and Qi et al. (2021). Additionally, the local DFMC SBAS corrections without and with the HOI mitigations are investigated. In the case without HOI mitigations, the \({\widehat{\delta \Delta I}}_{u}^{i}\) term for the available satellites in (6) will be omitted.

Local DFMC SBAS corrections

The DFMC SBAS corrections based on ICAO (2018) consist of the error corrections of satellite-position and satellite-clock errors. Typically, the satellite-position error correction is utilized to correct the satellite coordinate vectors in (3) whereas the satellite-clock error correction is used to reduce the estimated satellite-clock offset errors in (6). These corrections are broadcast via the geostationary SBAS satellite. The information is composed of several message types (MT) explained in ICAO (2018). The MT31-32 and MT34-37 are designed to support the DFMC SBAS service. In this work, we evaluate the performance without the SBAS satellites, the DFMC SBAS corrections have not been encoded into the message format.

The proposed local DFMC SBAS corrections with the HOI mitigations are described in the four steps as follows.

Step 1: Given the precise coordinate vectors of the receiver, the residual errors between the \(i\) th satellite and the \(j\) th receiver (\({\Delta R}_{j}^{i}\)) can be expressed below.

Method A (baseline): The residual errors (\({\Delta R}_{j,A}^{i}\)) at each time including the HOI effects are computed from

$$\Delta R_{{j,A}}^{i} = R_{j}^{i} - \left\| {{\mathbf{r}}_{j}^{i} - {\mathbf{r}}_{j} } \right\| + c\widehat{{\Delta t}}_{j}^{i} - \widehat{{\delta T}}_{j}^{i} ~ = \underbrace {{c\Delta t_{j} }}_{{\beta _{j} }} + \underbrace {{\left( {\delta \Delta I_{j}^{i} - c\delta \Delta t_{j}^{i} + \varepsilon _{j}^{i} } \right)}}_{{\beta _{{j,A}}^{i} }}$$
(7)

where \({R}_{j}^{i}\) is the smoothed pseudorange measurements between the \(i\) th satellite and the \(j\) th receiver based on ICAO SARPs (ICAO 2018), \({\mathbf{r}}_{j}^{i}\) is the coordinate vectors (in meter) of the \(i\) th satellite estimated by the code pseudorange at the \(j\) th receiver (\(\left[{x}_{j}^{i}, {y}_{j}^{i}, {z}_{j}^{i}\right]\)), and \({\mathbf{r}}_{j}\) is the precise coordinate vectors of the \(j\) th receiver (\(\left[{x}_{j}, {y}_{j}, {z}_{j}\right]\)), \({\widehat{\Delta t}}_{j}^{i}\) is the \(i\) th satellite clock offset (in sec) estimated by the code pseudorange at the \(j\) th receiver, \({\widehat{\delta T}}_{j}^{i}\) is the tropospheric delay of the \(i\) th satellite estimated based on MOPS (RTCA DO-229D, 2006) with the \(j\) th receiver, \({\beta }_{j}\) is the \(j\) th receiver biases (in meters), \({\delta \Delta t}_{j}^{i}\) is the residual satellite-clock offset errors of the \(i\) th satellite at the \(j\) th receiver (in sec), \({\delta \Delta I}_{j}^{i}\) is the HOI delays (in meter) for the \(i\) th satellite with the \(j\) th receiver, \({\varepsilon }_{j}^{i}\) is the total errors (in meter) of the estimated parameters, difference code biases, and smoothed errors, and \({\beta }_{j,A}^{i}\) is the residual errors of the \({\delta \Delta I}_{j}^{i}\), \({\delta \Delta t}_{j}^{i}\), and \({\varepsilon }_{j}^{i}\) parameters for the \(i\) th satellite with the \(j\) th receiver. Note that \({\text{t}}\) he \({c\delta \Delta t}^{i}\) term for the \(i\) th satellite is equivalent to the satellite-clock error correction of DFMC SBAS.

Method B (proposed): The residual errors (\({\Delta R}_{j,B}^{i}\)) of each time excluding the HOI effects are calculated by

$$\Delta R_{j,B}^{i} = \Delta R_{j,A}^{i} - \widehat{\delta \Delta I}_{j}^{i} = \underbrace {{c\Delta t_{j} }}_{{\beta_{j} }} - \underbrace {{\left( {c\delta \Delta t_{j}^{i} - \varepsilon_{j}^{i} } \right)}}_{{\beta_{j,B}^{i} }}$$
(8)

where \({\widehat{\delta \Delta I}}_{j}^{i}\) is the estimated HOI delays for the \(i\) th satellite with the \(j\) th receiver estimated from the local TECs based on the Klobuchar model (Setti et al. 2019), and \({\beta }_{j,B}^{i}\) is the \(i\) th satellite biases with the \(j\) th receiver excluding the HOI effects (in meters).

Step 2: From step 1, the \({\beta }_{j}\) term should be approximated before the \({\beta }_{j,A}^{i}\) term since the precise receiver coordinates can be calculated by historical data. Empirically, the \({\beta }_{j}\) should be the constant value for a short period (i.e., 5 min) due to the fixed receiver. In this work, the minimum sum-variance technique is applied to the TD observations (between-satellites, between-receivers, and between-epochs). Then each receiver \(j\) th (\({\widehat{\beta }}_{j}\)) at time \(t\) can be approximated by

$${\widehat{\beta }}_{j}\left(t\right)\approx \underset{{X}_{\mathit{min}}\le {b}_{j}\le {X}_{\mathit{max}}}{{\text{min}}}\left(\sum_{i=1}^{{N}_{s}}\sum_{\tau =1}^{{N}_{\tau }}\frac{{\left({\Delta R}_{j,A}^{i}\left(t-\tau \right)-{b}_{j}\right)}^{2}}{{N}_{s}+{N}_{\tau }}\right)$$
(9)

where \({X}_{min}\) and \({X}_{max}\) are the minimum and maximum values of the \(\Delta {R}_{j,A}^{i}\) at the time in {t − 1, t − 2, …, t − \({N}_{\tau }\)}, \({b}_{j}\) is an arbitrary mean value selected from min-to-max values, and \({N}_{\tau }\) is the time window size. The \({\beta }_{j}\) value is approximated depending on the \({N}_{\tau }\) and \({N}_{s}\) parameters. The \({N}_{s}\) parameter requires at least three satellites since the optimal mean can be computed. Additionally, the satellite and receiver coordinates are based on the WGS-84 ECEF coordinate system, Takasu (2013).

Step 3: Once the \({\beta }_{j}\) in (9) is obtained, the remaining term of \({\beta }_{j,A}^{i}\) corresponding the \(i\) th satellite with the reference receiver \(j\) th are computed. Since the baselines between the reference receivers for the SBAS ground segment are generally more than 100 km, the \({\beta }_{j,A}^{i}\) are estimated by the different values. Hence, in this work, the minimum sum-variance technique is also applied to the TD observations. The mean values of each satellite \(i\) th (\({\widehat{\beta }}^{i}\)) at time \(t\) are estimated from

$${\widehat{\beta }}^{i}\left(t\right)\approx \underset{{Q}_{\mathit{min}}\le {b}^{i}\le {Q}_{\mathit{max}}}{{\text{min}}}\left(\sum_{j=1}^{{N}_{r}}\sum_{\tau =1}^{{N}_{\tau }}\frac{{\left({\widehat{\beta }}_{j,A}^{i}\left(t-\tau \right)-{b}^{i}\right)}^{2}}{{N}_{r}+{N}_{\tau }}\right)$$
(10)

where \({\widehat{\beta }}_{j,A}^{i}\) is the remaining term of \({\Delta R}_{j,A}^{i}-{\widehat{\beta }}_{j}\) at time \(t\), \({Q}_{min}\) and \({Q}_{max}\) are the minimum and maximum values of the \({\widehat{\beta }}_{j}^{i}(t)\) where \(j\in \left\{1, 2, \dots ,{N}_{r}\right\}\), \({N}_{r}\) is the total number of available receivers, and \({b}^{i}\) is the arbitrary mean value selected from the min-to-max values. Additionally, the method B is also computed similar to the method A by changing the \({\Delta R}_{j,A}^{i}\) to \({\Delta R}_{j,B}^{i}\) and \({\beta }_{j,A}^{i}\) to \({\beta }_{j,B}^{i}\) instead.

Equations (9) is utilized to observe the \({\widehat{\beta }}_{j}^{i}\) of the \(i\) th satellites (multi-constellation) at the receivers \(j\) th, and (10) is applied to optimize the \({\widehat{\beta }}^{i}\) values by using the \({\widehat{\beta }}_{j}^{i}\) of many reference receivers with the long baselines. Thus, the \({\varepsilon }_{j}^{i}\) parameter should be reduced by the minimum sum-variance technique with the TD observations.

Step 4: The \({\widehat{\beta }}^{i}\) values, equivalent to the residual satellite-clock offset errors of \(i\) th satellite, cause the errors of the satellite coordinate vector interpolations (based on by the code pseudorange at each receiver). Hence, the local DFMC SBAS corrections, the \(i\) th satellite position error correction vectors (\({{\varvec{\updelta}}}^{i}\)) can be generated from the \({\widehat{\beta }}^{i}\) values. They are projected to the (X,Y,Z) coordinate by using the unit vectors of the virtual satellite coordinates. The \({{\varvec{\updelta}}}^{i}\) vectors of the \(i\) th satellite at each time can be computed from

$${{\varvec{\updelta}}}^{i}=\frac{{\mathbf{r}}^{i}}{\Vert {\mathbf{r}}^{i}\Vert }\times {\widehat{\beta }}^{i}$$
(11)

where \({\mathbf{r}}^{i}\) is the virtual coordinate vectors of the \(i\) th satellite (\(\left[{x}^{i}, {y}^{i}, {z}^{i}\right]\)) approximated by the code pseudorange at the master station.

DFMC SBAS integrity

The integrity parameters are utilized to determine the availabilities of the SBAS operation (RTCA DO-229D 2006). In practice, the horizontal protection level (HPL) and the vertical protection level (VPL) for the PA model at each time are, respectively, computed 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}}}$$
(12)

and

$$VPL={K}_{V}\times {d}_{U}$$
(13)

where \({K}_{H}\) and \({K}_{V}\) are the constant values of HPL and VPL, \({d}_{E}^{2}\) and \({d}_{N}^{2}\) are the variances of the true error distribution in the east and north axis, respectively, \({d}_{EN}^{2}\) is the covariance of the true error distribution in the east-north axis, \({d}_{U}\) is the standard deviation of the true error distribution in the up-axis. The true error distribution can be calculated by utilizing the inversion of the geometry matrix and weighting matrix derived in RTCA DO-229D (2006).

From (4)-(5), the weighting matrix \(\mathbf{W}\) is computed from the total variances including the DFRE values (in meters). They can be calculated based on the projection method with the observed data of the reference stations derived by Shao et al. (2020). This technique can bound the satellite correction errors as 99.9% of the time; it is also suitable for use in the multi-constellation. Hence, the DFRE values can be computed from

$${\sigma }_{i,DFRE}=\sqrt{{\left(\sqrt{{\lambda }_{1}}+\sqrt{{\mathcal{C}}_{i}}\right)}^{2}-{\mathcal{C}}_{i}^{2}+{C}_{i,t}}$$
(14)

where \({\lambda }_{1}\) is the maximal distance from the point on the ellipsoidal surface to the origin, \({\mathcal{C}}_{i}\) equals \({\mathbf{C}}_{i,ot}^{T}{\mathbf{C}}_{i,o}^{-1}{\mathbf{C}}_{i,ot}\), and \({C}_{i,t}\) is the inversion of the \(i\) th satellite clock variance. The covariance matrix (\({\mathbf{C}}_{i,o,(3\times 3)}\), \({\mathbf{C}}_{i,ot,(3\times 1)}\), and \({C}_{i,t,}\)) of the \(i\) th satellite orbit coordinate (\(o\)) with the clock error (\(t\)) is

$${\left[\begin{array}{cc}{\mathbf{C}}_{i,o}& {\mathbf{C}}_{i,ot}\\ {\mathbf{C}}_{i,ot}^{T}& {C}_{i,t}\end{array}\right]}_{\left(4\times 4\right)}={\left({\mathbf{A}}_{i}^{T}{{\varvec{\Lambda}}}_{i}{\mathbf{A}}_{i}\right)}^{-1}$$
(15)

\({\mathbf{A}}_{i}\) is the unit vector matrix (\(4\times {N}_{r}\)) of the \(i\) th satellite position error correction vectors with the \(j\) th receiver (\({\left[{\delta x}_{j}^{i}, \delta {y}_{j}^{i}, \delta {z}_{j}^{i}, 1\right]}^{T}\)), and \({{\varvec{\Lambda}}}_{i}\) is the weighting matrix where the diagonal matrix is the inversion of variances in the \(i\) th satellite-clock offset errors depending on the \({N}_{r}\) and \({N}_{\tau }\) parameters.

Methodology

In this work, the local DFMC SBAS corrections are generated based on the method A and B. We evaluate the system with and without the HOI delays on both ionospheric quiet and disturbed days.

System simulations

The assessments of DFMC SBAS are demonstrated in Fig. 1 consisting of six steps. First, the raw observations are fed as input. The first-order ionospheric delays are removed by the IF combinations. Second, the multipath noise is typically removed by the elevation mask component, whereby the satellites below the elevation mask are cut-off. Third, the IF carrier smoothing code technique is utilized to approximate the pseudorange measurements. Fourth, this step, both method A and B are tested. For method A, the satellite position error correction vectors and the integrity parameters are generated (Eq. 11,14) based on the residual errors including the HOI effects (Eq. 7), however, for method B, they are based on the residual errors with the HOI effects removed (Eq. 8). Fifth, the user coordinate vectors are calculated based on the SPP algorithm in (1), then the horizontal positioning errors (HPE) and vertical positioning errors (VPE) are computed based on the Haversine formula. Finally, the available percentages of DFMC SBAS are assessed by comparing the HPE and VPE values with the HPL and VPL values, respectively, on the Stanford diagram (Ogier 2021).

Fig. 1
figure 1

Block diagram of the DFMC SBAS assessments (Method A represents a baseline method, whereas Method B is the proposed method with additional HOI mitigation)

Full size image

In the simulations, the important parameters are described in Table 1. To investigate the HOI mitigation on DFMC SBAS, the quiet and disturbed days are selected based on the rate of TEC index (ROTI) values without the geomagnetic storms (the Dst above -30, Nose et al. 2015, and the Kp below 3, Jürgen et al. 2021). Then the observed periods of four days (non-consecutive) are considered explaining in Subsection-dataset. Normally, an observed sample is defined as 1 s, the SBAS service updated at least every 6 s (updating time) can be done. The elevation mask is defined as 15°, which is suitable for the multi-constellation case. Additionally, the \({N}_{\tau }\) of 300 samples (5 min) is a suitable sample for estimating the \({\widehat{\beta }}_{j}\) values in the near real-time processing.

Table 1 Important parameters for the DFMC SBAS simulation
Full size table

Figure 2 shows the selected stations covering Thailand. The reference stations include CHAN, CHMA, DPT9, NKSW, PJRK, SISK, SOKA, and UDON. The remaining three stations (CNBR, UTTD, and SRTN) are considered as the users. Since the DPT9 station is selected as the master station, distances between the master and user (CNBR, UTTD, and SRTN) stations are about 60, 432, and 529 km, respectively.

Fig. 2
figure 2

Locations of the reference (Δ) and user (□) stations in Thailand (https://gnss-portal.rtsd.mi.th/)

Full size image

Dataset

In the simulations, the GPS, Galileo, and QZSS satellites are utilized in the band of L1 and L5 frequencies. The GPS and Galileo satellites are in a medium earth orbit (MEO), however, the QZSS satellite is in a quasi-zenith orbit (QZO). Currently, the BDS-3 satellites with same bands are available, however, the BeiDou data are not available from the selected reference stations. Since not all GPS satellites broadcast the L5 frequency at present, fewer than four satellites are available at times. For example, at 12:00 and 22:30 h (local time: LT), three core systems with both L1 and L5 are illustrated in Fig. 3. The distinct colors indicate each constellation, i.e., GPS (G, O), Galileo (E, □), and QZSS (J, ◇) satellites. Normally, the GPS satellites during 12:00 LT have six satellites, but during 22:30 LT, only one satellite is available. Hence, the Galileo and QZSS satellites are utilized to avoid the positioning outage as well as the high DOP effects.

Fig. 3
figure 3

Sky plots of the GPS (G, O), Galileo (E, □), and QZSS (J, ◇) satellites with both L1 and L5 at 12:00 LT (left) and 22:30 LT (right)

Full size image

Figure 4 shows the ROTI values of the CNBR station. The top figure (quiet days) consists of the DOY of the 36, 38, 45, and 50 in 2022 labeled as the day numbers of the 1 to 4, respectively. Each ROTI line is computed by the time window of 5 min and the elevation masks of 30°. The data gap is due to missing observational data. Similarly, the bottom figure (disturbed days) shows the DOY of the 60, 62, 74, and 75 in 2022. Normally, the ROTIs on the quiet days are below 0.5 TECU/min with low fluctuation, whereas on the disturbed days are above 0.5 TECU/min after sunset (Schaer 1999). Although there is no geomagnetic storm on the disturbed days during the studied period, the ROTIs after sunset indicate local ionospheric disturbances due to the equatorial plasma bubbles (EPB).

Fig. 4
figure 4

ROTI plots of the CNBR station on the quiet (top) and disturbed (bottom) days (Day 1 to 4 represents DOY 36, 38, 45, and 50 in 2022)

Full size image

Similarly, Fig. 5 shows the ROTI plots at UTTD and SRTN stations. Mostly, the ROTIs with the large fluctuation are observed after sunset. At the SRTN station, relatively lower ROTI levels are seen on the first day. Interestingly, the SRTN station generally experiences lower ROTIs than the UTTD station. These results indicate that the ionospheric disturbances vary by location.

Fig. 5
figure 5

ROTI plots of the UTTD (top) and SRTN (bottom) stations on the disturbed days (Day 1 to 4 represents DOY 36, 38, 45, and 50 in 2022)

Full size image

Results and discussions

The DFMC SBAS availabilities with the quiet and disturbed events are evaluated. The local DFMC SBAS correction and integrity with and without the HOI delays are utilized for the user segments such as the CNBR, SRTN, and UTTD stations in Fig. 2. Then the simulation results are demonstrated and discussed as follows.

DFMC SBAS parameters

Figure 6 shows the satellite position error correction (\({\delta x}^{i}, {\delta y}^{i}, {\delta z}^{i}\)) and integrity parameter (\({\sigma }_{i,DFRE}\)) from the method B on the quiet days. The six reference stations (CHAN, CHMA, DPT9, NKSW, PJRK, and UDON) are utilized for the computation of the master station processing. The distinct color lines indicate each satellite number of the GPS (G04 and G10), Galileo (E02 and E03), and QZSS (J02 and J03) constellation. The satellite position error correction vectors (\({{\varvec{\updelta}}}^{i}\)) are computed by (11) with the virtual coordinate vectors of the \(i\) th satellite approximated by the code pseudorange at DPT9, while the integrity parameters (\({\sigma }_{i,DFRE}\)) are calculated by (14).

Fig. 6
figure 6

Satellite position error correction and integrity parameter with the method B for each satellite on the quiet days (Day 1 to 4 represents DOY 36, 38, 45, and 50 in 2022)

Full size image

In Fig. 6, the \({\delta x}^{i}\), \({\delta y}^{i}\), and \({\delta z}^{i}\) values are within ± 10 m. The position error correction of J02 and J03 show the least fluctuations (± 2 m), possibly due to the trajectory of the quasi-zenith orbit (QZO) with the local quiet event. The \({\sigma }_{i,DFRE}\) values are more than 3 m are caused from some spikes or outliers of the position error correction. For example, \({\delta y}^{i}\) and \({\delta z}^{i}\) vectors of E02 and E03 are on the first day. Hence, the \({\sigma }_{i,DFRE}\) parameter is utilized as the error indicator related to the local DFMC SBAS corrections.

Figure 7 shows the 95th percentile of \({\sigma }_{i,DFRE}\) versus the number of reference stations (4, 6 and 8) for each satellite. The distinct color bars are computed on four days, whereby the dotted bars are the averages of all color bars. On the quiet days, the \({\sigma }_{i,DFRE}\) values of each satellite are obviously decreased as the number of reference stations increases showing the averaged \({\sigma }_{i,DFRE}\) of 1.08, 0.26, and 0.12 m for 4, 6, and 8 reference stations, respectively. On the disturbed days, the similar trend of \({\sigma }_{i,DFRE}\) are still seen, however, the averaged \({\sigma }_{i,DFRE}\) are relatively larger compared with the quiet days. In the following analysis, we select the number of the reference stations to be 6.

Fig. 7
figure 7

The 95th percentile of \({\sigma }_{i,DFRE}\) and averages (Avg) based on the method B versus the number of reference stations on the quiet (top) and disturbed (bottom) days

Full size image

User positioning errors

The preliminary DFMC SBAS operation can be evaluated by the user positioning errors (user processing). The precise coordinate vectors of the reference station are defined by the AUSPOS Online GPS Processing Service on the website https://gnss.ga.gov.au/auspos. This website is administered by Geoscience Australia; the users need to submit the observed data with the receiver independent exchange format (RINEX). The raw GPS pseudorange values are utilized with the reference stations of the international GNSS service (IGS) and the Asia Pacific Reference Frame (APREF) to compute the precise receiver coordinate based on the double-difference strategy (relative positioning); the accuracy is in the centimeter levels. Then an AUSPOS report will be emailed to users with the Geocentric Datum of Australia 2020 (GDA2020), Geocentric Datum of Australia 1994 (GDA94), and International Terrestrial Reference Frame (ITRF) coordinates. Further information is described by Malys and Slater (1994). Figure 8 shows the number of satellites, HPE, and VPE at the user station (CNBR) processing on the quiet and disturbed days. The ‘SPS’ legend indicates the standard positioning system (SPS) based on the GPS L1 frequency and the Klobuchar model, whereas the ‘DFMC SBAS-A’ and ‘DFMC SBAS-B’ legends represent the method A and B, respectively. The SPS method is used to validate the DFMC SBAS algorithm. On the quiet days, the number of available satellites with both L1 and L5 tends to decrease from after sunset to midnight. The HPEs are typically lower than 5 m, but the VPEs are typically 10 m or less. As expected, the VPEs are higher than HPEs by about 1.5 m due to the geography and ionospheric delays. Both HPE and VPE increase significantly as the number of satellites decreases. Importantly, the DFMC SBAS-B method can better improve the errors than the DFMC SBAS-A.

Fig. 8
figure 8

Number of satellites (top panel) and user positioning errors of the SPS, DFMC SBAS-A, and DFMC SBAS-B methods in terms of HPE (middle panel) and VPE (bottom panel) on the quiet (left panel) and disturbed (right panel) days at CNBR station (Day 1 to 4 represents DOY 36, 38, 45, and 50 in 2022)

Full size image

For ease of comparison, the 95th percentile of HPE and VPE at three stations (CNBR, UTTD, SRTN) are summarized in Fig. 9. Each graph is obtained from the positioning errors for four days. On the quiet days, both the HPEs and VPEs from the method B are relatively lower than the method A. Compared with DFMC SBAS-A, the DFMC SBAS-B method can improve the HPEs at the CNBR, UTTD, and SRTN stations by 0.14 (13%), 0.07 (11%), and 0.05 (5%) m, and the VPEs about 0.21 (8%), 0.15 (12%), and 0.16 (10%) m, respectively. The improved percentage is computed by (1 – (DFMC SBAS-B/DFMC SBAS-A)) × 100%. Similarly, on the disturbed days, the HPEs and VPEs of each user station are also improved, but the improvements are slightly less. The HPEs of the DFMC SBAS-B over DFMC SBAS-A results are about 0.12 (10%), 0.07 (10%), and 0.05 (5%) m at the CNBR, UTTD, and SRTN stations, respectively. While the VPEs at the CNBR, UTTD, and SRTN stations are about 0.25 (8%), 0.13 (8%), and 0.15 (7%) m, respectively. The improvement due to Method B illustrates that when the HOI effects are mitigated, the positioning errors are reduced by 5–20 cm.

Fig. 9
figure 9

User positioning errors in the 95.th percentile HPE (top panel) and VPE (bottom panel) on the quiet (left panel) and disturbed (right panel) days for three user stations (CNBR, UTTD, and SRTN)

Full size image

Here, we provide some discussions on the results in Fig. 9, the DFMC SBAS corrections with the local HOI mitigation can significantly enhance the user positioning accuracy on both quiet and disturbed days in the range of 7%-13%, but the improvements on the disturbed days are slightly smaller than those on the quiet days by about 2%. It is clear from the results in Fig. 7 (disturbed), larger HOI delays are observed for any number of reference stations. Therefore, the estimated HOI delays on the disturbed days are less reliable than on the quiet days leading to the degradation on the precision of the DFMC SBAS corrections. Since we estimate the HOI values from Klobuchar model, more accurate HOI can be obtained from the actual pseudorange observations.

Protection levels

In practice, the HPL and VPL, rather than HPE and VPE values are estimated by the airborne receivers since actual positions are not known. These values can be utilized to assess the availability of the SBAS operation. Figure 10 shows the number of satellites, HPL, and VPL results of the user station (CNBR) computed by (12) and (13) with \({K}_{H}=6.0\) and \({K}_{V}=5.33\), respectively. Similar trends in HPL and VPL are also seen here, showing the increased protection levels up to 40 and almost 60 m, respectively, when the number of satellites is decreased below 10. Similar to Fig. 9, the protection levels (PL) based on the 95th percentile are summarized in Fig. 11.

Fig. 10
figure 10

Number of satellites (top) and the user protection levels in terms of HPL (middle) and VPL (bottom) on the disturbed days at the CNBR station (Day 1 to 4 represents DOY 36, 38, 45, and 50 in 2022)

Full size image
Fig. 11
figure 11

95.th percentile HPL (top panel) and VPL (bottom panel) on the quiet (left panel) and disturbed (right panel) days for three user stations (CNBR, UTTD, and SRTN)

Full size image

In Fig. 11, we investigate the HPLs and VPLs similarly to Fig. 10 for the three stations. On the quiet days, the HPLs and VPLs with the HOI effects for all user stations can be reduced by the DFMC SBAS-B method by about 21% and 16%, respectively. In contrast, on the disturbed days, the improved HPLs with DFMC SBAS-B for the CNBR, UTTD, and SRTN stations are about 25%, 24%, and 26%, while the improved VPLs are about 19%, 20%, and 21%, respectively. Similarly, the protection levels follow the trends of the user positioning errors in Fig. 9, but the improvements on the disturbed days are slightly larger than those on the quiet days. It is reasonable that the large HOI delays on the disturbed days caused the overestimated protection levels with DFMC SBAS-A.

DFMC SBAS performances for LPV-200 and CAT-I phases

Finally, the ICAO requirements of the PA model (ICAO SARPS 2006) are considered. For the LPV-200 (equivalent CAT-I) and CAT-I categories, the vertical alert limit (AL) bounds are determined as 35 and 15 m, respectively. In the assessments of DFMC SBAS, both the positioning errors and protection levels are plotted in the Stanford diagram.

Figure 12 compares the Stanford diagrams between the DFMC SBAS-A and DFMC SBAS-B methods at the user station (CNBR) on four quiet days. The MI and HMI stand for misleading information and the hazard misleading information, respectively. The color bar indicates the values of the histogram for four days. The availabilities of DFMC SBAS-A in the LPV-200 and CAT-I categories are about 99.77% and 97.72%, and of DFMC SBAS-B in both categories are 99.98% and 99.46%, respectively. Evidently, the DFMC SBAS-B method can improve the DFMC SBAS-A method by 0.21% and 1.74% with the LPV-200 and CAT-I categories, respectively.

Fig. 12
figure 12

Stanford diagrams of the DFMC SBAS-A (top) and DFMC SBAS-B (bottom) methods at the CNBR station on the four quiet days

Full size image

Table 2 shows the summarized availabilities of DFMC SBAS for three users. As expected, the local ionospheric disturbances decrease the availability at each user station. On the disturbed days, the available percentages are reduced depending on the baselines between the users and the master station. For example, the available (%)/baselines (km) in CAT-I are about 97.75/60, 97.42/432, and 97.20/529 at the CNBR, UTTD, and SRTN stations, respectively. Note that the baselines can reduce the availability below 1%. Additionally, if we assess the available percentages in LPV-200 and CAT-I with the ICAO requirements of 99.99%, the available DFMC SBAS with the method B in Thailand cannot yet be achieved.

Table 2 Available percentages of the DFMC SBAS-B/DFMC SBAS-A methods at the CNBR, UTTD, and SRTN stations for the LPV-200 and CAT-I categories, respectively
Full size table

We can roughly divide the correction of the HOI delay problem into: receiver-end (user) processing of satellite signals and the base (master) station or control center processing of satellite signals. At the master station, the HOIs need to be removed from the residual errors, then the DFMC SBAS corrections (satellite orbit errors, satellite clock errors) are generated. In the process, the estimated HOI delays of each satellite need to be computed at each reference station. On the other hand, on the receiver side, the HOIs need to be estimated and removed from the pseudorange observations. As space weather plays an important role in the variation of ionospheric delays, HOI included, both types of processing need to consider the influencing factors of the HOIs such as solar activity, geomagnetic activity and other space weather parameters. The annual, seasonal and temporal variation of HOIs need to studied and, importantly, HOI models should be developed. Location (latitude and longitude) factors are also of importance. During the upcoming solar maximum, the ionospheric delays tend to increase and during geomagnetic storms typically as a result of solar flare events, the ionospheric irregularity (EPB events) often occurs, hence, more data collection and analysis are continually required. The master and receiver sides will benefit from such HOI models. With disturbance detection, both the receiver and the master station can adapt to more accurate HOI values.

Conclusion

In this work, we evaluate the DFMC SBAS with and without the HOI mitigation in Thailand. The procedures for local DFMC SBAS corrections and integrity parameters are proposed and then generated by using the local GNSS reference stations. The estimated local HOI delays are used in the HOI mitigation showing that the improvements of user positioning errors are up to in 13%, but the improvements on the disturbed days are slightly less than those on the quiet days. The results on the protection levels follow the trends of positioning errors. The availability of DFMC SBAS with the local DFMC SBAS corrections for both LPV-200 and CAT-I categories are 99.98% and 99.46%, respectively. It is worth noting that the suitable SBAS operation with the customer requirements must be below 99%. The HOI mitigation is hence beneficial to DFMC SBAS operation. Further studies and analyzes of DFMC SBAS during the upcoming solar maximum need to be made together with advanced algorithms and multi-constellations to achieve availability levels according to ICAO standards.

In addition, the integrity parameter (UIRE) needs to be created for the (new) regional model based on the actual observations and the local ionospheric disturbances. The regional model will improve the performance of DFMC SBAS particularly on a local scale.