Overbounding of Near Real-Time Estimated Ionospheric Gradient Slope in Low-Latitude Regions

3.1 Test Statistic and Threshold

First, the processing component of each station belonging to the network receives GNSS dual-frequency code and carrier-phase measurements and calculates the geometry-free combination of the carrier-phase measurements at the current epoch t as follows:

${\hat{I}}_r^j\left(t\right)=\frac{{\left(f_2\right)}^2}{{\left(f_1\right)}^2-{\left(f_2\right)}^2}·\left({\varphi }_{f_1,r}^j\left(t\right)-{\varphi }_{f_2,r}^j\left(t\right)\right)$ 1

where ${\varphi }_{f_i,r}^j\left(t\right)$ is the carrier-phase measurement in meters for frequency f_i (e.g., L1/E1 or L5/E5a), j is a single satellite from the GPS or other satellite constellation, and r is the station that calculates the geometry-free combination (i.e., estimate of the slant ionospheric delay with errors such as carrier-phase ambiguities or hardware biases).

The test statistics of the monitoring network are the rates of change of the estimated ionospheric delays corresponding to each satellite j and each station r, computed as follows:

$Ites{t}_r^j\left(t\right)=\frac{{\hat{I}}_r^j\left(t\right)-{\hat{I}}_r^j(t-1)}{\Delta T}$ 2

Here, ${\hat{I}}_r^j\left(t\right)$ is the slant ionospheric delay estimation (with errors such as carrier-phase ambiguities or hardware biases) calculated via Equation (1) for station r, satellite j, and epoch t. ${\hat{I}}_r^j(t-1)$ is the slant ionospheric delay estimation with errors for the previous epoch, and ΔT is the time difference between two consecutive epochs. Note that the errors in the slant ionospheric delays computed via Equation (1) remain constant over time if no cycle slips occur. As a result, these errors cancel out when the test statistic is computed via the first derivative over time (Equation (2)).

The test statistics are compared with predefined thresholds derived from historical data recorded in nominal ionospheric conditions from each station r and integrity requirements. These thresholds, which depend on satellite elevation, determine whether there is a significant perturbation affecting each of the satellite–station pairs in real time and are computed as follows:

$Th{r}_r\left({\theta }_m\right)={\mu }_{Ites{t}_r}\left({\theta }_m\right)\pm k_{fa}·I_{monitor,r}\left({\theta }_m\right)·{\sigma }_{Ites{t}_r}\left({\theta }_m\right)$ 3

where θ_m with m = 1, 2, … M represents the elevation bins, ${\mu }_{Ites{t}_r}\left({\theta }_m\right)$ is the mean of the $Ites{t}_r^j\left(t\right)$ distribution computed in each elevation bin, ${\sigma }_{Ites{t}_r}\left({\theta }_m\right)$ is the standard deviation, and I_monitor,r(θ_m) is the inflation factor calculated to overbound the tails of the distribution per elevation bin. The false alert multiplier, k_fa, is computed with the probability of false alert, P_fa. The detailed methodology for deriving the detection thresholds has been explained in the work by Caamano et al. (2021) and is not repeated in this paper. The historical data for nominal conditions used in this paper to derive the monitoring thresholds are described in Section 5.

3.2 Executive Monitoring

The detection information from each station observing a given satellite is shared within the network in real time. Here, the network distinguishes between three possible cases: (i) several monitoring stations have detected the same ionospheric gradient, (ii) none of the monitoring stations have detected an anomalous ionospheric gradient, and (iii) several monitoring stations have detected anomalous gradients, but the network is not able to ensure that this gradient is the same for all of the stations.

To determine whether several monitoring stations have detected the same ionospheric gradient, we compute the cross-correlation coefficient $\left({\alpha }_r^j\left(t\right)\right)$ in real time between test statistics belonging to the same satellite and different network receivers, as explained by Caamano et al. (2021). A cross-correlation coefficient of 1 between test statistics computed at two different stations indicates that the perturbation is the same, but delayed by a certain time interval $\left({\kappa }_r^j\left(t\right)\right)$ . A cross-correlation coefficient of 0 indicates that the perturbations at these two stations appear to be unrelated. Therefore, these cross-correlation coefficients ${\alpha }_r^j\left(t\right)$ are compared in each epoch with a minimum value α_min (e.g., 0.5), below which the perturbation occurring at station r is not considered to be the same as the perturbation occurring at the reference station (typically the first station that detects the gradient). Until ${\alpha }_r^j\left(t\right)$ reaches the value chosen for α_min, the network cannot determine whether both stations detected the same anomalous ionospheric gradient. Therefore, it declares a “Warning” to indicate to the GBAS stations in that area that they should use the “worst-case” threat model parameters. Once the algorithm converges, it outputs the cross-correlation coefficient ${\alpha }_r^j\left(t\right)$ and the time delay ${\kappa }_r^j\left(t\right)$ . While the convergence criteria are met, both parameters are very similar from one epoch to the next and can be considered constant. For the times when the cross-correlation coefficients no longer meet the convergence criteria but the gradient is still detected, the network again indicates the use of “worst-case” assumptions. The algorithms used to calculate the cross-correlation coefficients and time delays in real time and an explanation of the convergence criteria have been reported by Caamano et al. (2021).

For the case in which several monitoring stations have detected the same ionospheric gradient for the same satellite and these stations are within a 200-km radius (to ensure that the gradient has not changed significantly during its propagation; see work by Caamano et al. (2021) for a detailed explanation), a central processor estimates the gradient parameters. The parameters are estimated by using the cross-correlation coefficients $\left({\alpha }_r^j\left(t\right)\right)$ and time delays $\left({\kappa }_r^j\left(t\right)\right)$ between the different stations involved in the cross-correlation process. This gradient parameter estimation process is performed for each satellite and requires that at least three stations have detected the same gradient. The first parameters that are estimated are the speed $\left({\hat{v}}_{iono}^j\right)$ and direction $\left({\hat{d}}_{iono}^j\right)$ of the gradient. These parameters can be estimated via a weighted-least-squares approach, as proposed by Caamano et al. (2021), with the following inputs: the position of the IPPs captured when the gradient is detected by each satellite–station pair, the time delays between stations, and the cross-correlation coefficients between test statistics computed in real time. Once the speed vector of the gradient is known, the slant slope of the gradient estimated by the considered satellite and station r is determined by the following geometrical relationship:

${\hat{g}}_{iono,r}^j\left(t\right)=\frac{Ites{t}_r^j\left(t\right)}{\Delta v_r^j\left(t\right)}$ 4

where $\Delta v_r^j\left(t\right)$ is the relative speed between the gradient and the IPP, projected in the direction of gradient propagation. In this way, the movement of the IPP relative to the gradient is corrected, and the gradients are not underestimated. The width of the gradient $\left({\hat{w}}_{iono,r}^j\right)$ is computed by multiplying the relative speed between the gradient and the IPPs by the time (in seconds) that has passed since the gradient was detected until the current time t. Note that the total width of the gradient can only be known when the whole gradient has passed and station r is no longer detecting the gradient.

If fewer than three stations have detected the gradient or the network is not able to ensure that the gradient detected is the same for all stations, the network triggers a “Warning” to indicate that the GBAS stations should use “worst-case” assumptions. If none of the monitoring stations have detected any anomalous ionospheric gradients, the network transmits the largest gradient that could be affecting the supported GBAS stations without being detected, or the “minimum detectable gradient” (MDG). Calculation of the MDG is part of future work.

Each GBAS station would be then responsible for using the network information to adjust the integrity parameters and support the already existing equipment in covering the GBAS approaches at each airport. Specifically, we propose to use a variable maximum gradient slope per satellite to calculate anomalous ionosphere-induced differential range errors within the PDGS algorithm (GAST C) instead of using the constant value from the threat model as is currently done. This variable maximum gradient slope would depend on the information received from the network. However, the estimated gradient slope derived from Equation (4) cannot be directly used in GBAS algorithms, as the associated estimation errors could compromise GBAS integrity. Therefore, in this work, we propose a methodology to overbound the estimated gradient slope in real time and meet GBAS integrity requirements. Incorporating this information into the PDGS and assessing its impact on system availability will be addressed in future work. Figure 2 presents the network decision logic described above.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Flow chart of the network decision logic.

4.1 Probability of Non-Bounded Errors

This work defines the probability of non-bounded errors as the probability of estimating error bounds erroneously in a way that compromises GBAS integrity. Because the goal of this work is to substitute the current conservative threat model for the network output, the network itself must cover the integrity budget allocated to the threat model. Although the ionospheric threat model is considered to always be true and, in principle, does not have a fixed associated integrity budget, Section D.3 of Eurocae (2013) (“GBAS integrity risk not covered by protection levels”) states that: “The GBAS Ground Subsystem integrity risk is less than 1.5 × 10⁻⁷ in any one approach. […] The integrity risk due to Ground Subsystem failures can be divided into risks associated with undetected events of environmental anomalies (Note 1) and Ground Subsystem failures resulting in erroneous GBAS messages. […] Anomalous ionosphere, anomalous troposphere, excessive radio-frequency interference (RFI) and excessive multipath are notably considered as environmental anomalies.” This description corresponds to the probability of loss of integrity allocated to faults different from individual satellite failures not bounded by any PL, which is 1 × 10 8 in the work by Pullen and Enge (2007). A part of this integrity budget should be sub-allocated to anomalous ionosphere. The GBAS manufacturer selects these allocations according to the characteristics of the GBAS station (e.g., detectability of the implemented monitors). As a simplification, this work considers a sub-allocated probability of loss of integrity per 150-s approach for abnormal ionospheric activity (P_{LOI,abnormal_iono}) of 1 × 10⁻⁸ from the total of 1.5 × 10⁻⁷ (integrity risk allocation for H2 conditions as reported by Yoon et al. (2019)). Sub-allocation for other anomalous conditions, such as anomalous troposphere, has not been considered. Therefore, we consider a probability of non-bounded errors (P_ne) equal to 1 × 10⁻⁸, assuming a value of 1 for the prior probability of an anomalous ionospheric gradient occurring.

4.2 Derivation of an Uncertainty Model for Estimated Gradient Parameters

The estimation process described in Section 3 is not free from errors that may cause the estimated values to deviate from the true values. These errors arise from several factors such as the presence of multipath and noise in the measurements, approximations and simplifications in the mathematical models, and insufficient time resolution of the measurements. In Section 7, we study these errors and their causes with simulated and real data.

To be able to substitute the conservative threat model by the gradient parameters estimated from the network, a statistical model of the estimation uncertainty is needed. To derive such a model, the gradient parameter estimation errors must first be determined. This work focuses specifically on the estimation errors of the gradient slope, as this parameter serves as the interface between the network and the GBAS stations. The estimation error of the slope for station r and satellite j can be computed as follows:

${\Delta }_{g_{iono,r}}^j\left(t\right)={\hat{g}}_{iono,r}^j\left(t\right)-g_{iono,r}^j\left(t\right)$ 5

where ${\hat{g}}_{iono,r}^j\left(t\right)$ is the maximum estimated slant slope for satellite j and station r and $g_{iono,r}^j\left(t\right)$ is the true value for the slant slope at the same epoch, for the same satellite and the same station. However, because the true values of $g_{iono,r}^j$ are unknown, Monte Carlo simulations are conducted by varying all of the gradient parameters shown in Figure 1 to calculate these errors in a controlled environment. In each simulation, $g_{iono,r}^j$ is known and varies within the range of values specified by the Brazilian threat model in Table 1. When the slope is downward, a minus sign is applied. The methodology described in Section 3 is applied to each simulation, allowing for an estimate of ${\hat{g}}_{iono,r}^j$ alongside the varied $g_{iono,r}^j$ and the other gradient parameters. This enables an offline computation of ${\Delta }_{g_{iono,r}}^j$ for all satellites, stations, and variations of all gradient parameters in the threat model. Estimation errors for all other gradient parameters can be similarly determined by subtracting the simulated value from the estimated value for each parameter. Section 6 provides a detailed description of the different simulation steps.

Once ${\Delta }_{g_{iono,r}}^j$ values have been computed for all satellites, stations, and variations of all gradient parameters, the algorithm sorts these error values in bins of estimated values in order to allow for variable behavior under different levels of ionospheric intensity while maintaining a sufficient number of samples per bin to obtain reliable statistics. Errors associated with both positive and negative estimated slant slopes of the same absolute value are grouped into the same bin. The number of samples in each bin and the interval of estimated values covered by each bin depend on the characteristics of the errors obtained. Section 7.2.3 explains, using simulated data, how the bin size is selected.

The estimation error uncertainty in each bin is modeled as a Gaussian distribution, which is commonly used in GNSS augmentation systems because of its simplicity for defining the error distribution with only two parameters (mean and standard deviation). Therefore, once the estimation errors are divided into bins, the estimation error samples (e.g., ${\Delta }_{g_{iono,r}}^j$ with $j\in 1,\dots ,J$ and $r\in 1,\dots ,ℛ$ ) are normalized by using the sample mean, ${\mu }_{g_{iono,m}}$ , and the sample standard deviation, ${\sigma }_{g_{iono,m}}$ , in each bin with $m=1,2,\dots ,ℳ$ representing the bins.

Typically, after normalization, the distribution is practically homogeneous, except for the bins corresponding to low elevations because they have less data owing to the manner in which the gradients are simulated. Section 7.2.3 explains this phenomenon with simulated data. Therefore, a single inflation factor, $I_{g_{iono}}$ , is calculated for the entire distribution to properly overbound the non-Gaussian tails of the normalized ${\Delta }_{g_{iono}}$ distribution formed by the normalized samples in all bins. This step is carried out following the methodology explained in Section 2.2 of the work by Xie (2004), which is one of the main references for the overbounding technique in the GBAS context. Thus, the inflated standard deviation per bin can be calculated as follows:

${\sigma }_{g_{iono,overbound,m}}=\left|{\mu }_{g_{iono,m}}\right|+I_{g_{iono\hspace{0.5em}}}·{\sigma }_{g_{iono,m}}$ 6

Then, when the algorithm is running in real time and the network estimates the value of a gradient parameter (e.g,. the slope), the previously determined inflation factors are applied to compute the appropriately overbounded estimated parameter as follows:

${\hat{g}}_{iono,overbound,r}^j\left(t\right)={\hat{g}}_{iono\text{,}r}^j\left(t\right)\pm k_{ne}·{\sigma }_{g_{iono,overbound,m}}$ 7

where k_ne is the scalar multiplier needed to meet the required probability of non-bounded errors described in Section 4.1, P_ne, computed from the inverse of the standard normal cumulative distribution, Q⁻¹(P_ne), and ${\sigma }_{g_{iono,overbound,m}}$ is the over-bounded standard deviation corresponding to the data bin m associated with the estimated slant slope ${\hat{g}}_{iono,r}^j\left(t\right)$ . Note that when the estimated slant slope ${\hat{g}}_{iono,r}^j\left(t\right)$ is positive, the second part of the formula is added to it, and when ${\hat{g}}_{iono,r}^j\left(t\right)$ is negative, the second part of the formula is subtracted from it.

Figure 3 depicts the process for calculating the overbounded estimated gradient slope. First, the real-time algorithm estimates the gradient slant slope $\left({\hat{g}}_{iono,r}^j\left(t\right)\right)$ . Then, the algorithm calculates the bin to which the estimated slant slope belongs (e.g., m = 3). This slope has an associated estimation error $\left({\Delta }_{g_{iono,r}}^j\right)$ that is unknown in the real case because the true slope value is also unknown. Therefore, the statistics derived from the ${\Delta }_{g_{iono,r}}^j$ values obtained from simulations, computed via $g_{iono,r}^j$ from the Monte Carlo simulations as per Equation (5), are used to overbound the estimated slant slope. This step is performed by using the mean and standard deviation of the bin to which the estimate with real data belongs (${\mu }_{g_{iono,3}}$ and ${\sigma }_{g_{iono,3}}$ in the example of Figure 3, respectively) to calculate the inflated sigma (${\sigma }_{g_{iono,overbound\text{,}3}}$ in Figure 3) (see Equation (6)). Finally, the previous elements and P_ne are used to compute the overbounded slant slope of the gradient (see Equation (7)). Section 7 presents an explanation of the above-described process with simulated and real data.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Process of calculating the overbounded estimated gradient parameters.

5.1 Data Description

As a representative monitoring network for the equatorial region, we selected eight stations situated in the state of São Paulo, Brazil. Figure 4 depicts the locations of these stations. Public data are available for this network, whose stations belong to the Brazilian network for continuous monitoring of GPS (RBMC), at a 15-s sampling rate for both L1 and L2 frequencies and GPS satellites (Instituto Brasileiro de Geografia e Estatística, 2014).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Locations of stations in Brazil.

As can be observed, the distances between most of the stations exceed 100 km, which can affect the estimation of the gradient parameters. Additionally, a sampling rate of 15 s is typically considered low for this application. However, the baselines and sampling rate required to successfully estimate the gradient parameters depend on the physical characteristics of the ionospheric perturbations. According to Caamano et al. (2021), the ionospheric perturbations in Alaska were small in magnitude but fast-changing and moving at high speeds, necessitating a sampling rate of least 1 s and baselines of approximately 5 km. These values were at the limit of what could be used to estimate the gradient parameters. However, as will be discussed in Section 7.3.2, ionospheric disturbances in the equatorial region, typically produced by EPBs, move at lower speeds (as reported by Yoon et al. (2017)) and change more slowly with distance and time. Under these conditions, baselines below 200 km and a 15-s sampling rate are suitable for gradient parameter estimation. Still, estimation errors can be expected, as baseline distances and sampling rates should ideally be adapted to the specific characteristics of the ionospheric disturbances in the region. For this reason, although all stations are used for detection, only stations “poli,” “mgin,” “sjsp,” “chpi,” and “uba1” are utilized for the estimation process.

Owing to the limited availability of data recorded at L1/L5 frequencies and other constellations (e.g., Galileo) during active ionospheric conditions, real measurements on L1 and L2 frequencies from only the GPS constellation were used for this study. Nowadays, more stations are available in this area with shorter baselines, and they provide data at a higher sampling rate. However, this work used the data available from 2014 because this year corresponded to a solar maximum and has been well studied in the literature (Yoon et al., 2017; Yoon, Kim et al., 2020).

Note that the network depicted in Figure 4 can support one or more GBAS stations located within its coverage area. For detection, assuming that a trigger in one station indicates a true anomalous gradient and not a false alert, the GBAS stations located within the polygon formed by “spbo,” “neia,” “uba1,” “chpi,” “mgin,” and “eesc” would be supported. For estimation, given that EPBs travel along MODIP lines (i.e., approximately 58° clockwise from the North Pole in this area), as explained in Section 2, and that at least three stations are required to estimate the gradient parameters (see Section 3), with baselines below 200 km, the GBAS stations situated within the polygon formed by “mgin,” “sjsp,” “uba1,” and “chpi” would benefit from the estimation of gradient parameters in addition to the detection feature.

5.2 Date Selection

The dates selected from all available data attend to two different purposes: (i) to study a real anomalous ionospheric gradient measured by the network depicted in Figure 4 and (ii) to derive a monitoring threshold for each station in the network. Both the active and nominal ionospheric days were determined based on the along-arc total electron content (TEC) rate (AATR), computed as reported by Juan et al. (2018) for one of the stations under study, “chpi.” Figures 5(a) and 5(b) show the specific AATR values for this station during active and nominal ionospheric periods, respectively.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

AATR values for reference station “chpi”. (a) Active conditions during 60 days of year 2014. (b) Nominal conditions during 20 days of year 2014.

The threshold of 0.25 TECU/min represents the value below which the ionosphere is considered “quiet” or nominal, and the threshold of 1 TECU/min is the value above which the ionosphere is considered very “active.” Figure 5(a) also shows the variation of the AATR values relative to the local time in Brazil during a selected active day.

As an AATR threshold to determine “active” days, 1 TEC unit (TECU) per minute was used, which corresponds to very high activity. Among the days for which the AATR exceeded 1 TECU/min, day 69 of year 2014, i.e., March 10, 2014, was selected. The zoomed-in part of Figure 5(a) shows the variation of the AATR relative to the local time in Brazil. This date falls within the time period associated with the maximum of Solar Cycle 24, when the largest anomalous ionospheric gradients were detected in Brazil (Yoon et al., 2017; Yoon, Kim et al., 2020).

For deriving the monitoring thresholds, both nominal days (i.e., 24 h in which no significant ionospheric activity was observed) and nominal hours of active days (i.e., 12 h during the daytime of active days) were selected. Note that the active hours of active days were not used for deriving the thresholds, because this high ionospheric activity is considered off-nominal and is what we aim to detect, even though this high activity is very frequent during the nighttime of active periods (see Figure 5(a)). Therefore, the first 10 “quiet” days were manually selected after the active period in March 2014. The calm period was selected after the active period and not before because the active period lasted from approximately December 2013 (Sanz et al., 2014) to April 2014, and in December 2013, the measurements from some of the stations in Figure 4 were not available. Then, 12 nominal hours belonging to 20 “active” days were selected to account for satellite geometry variability and to obtain more measurements. Specifically, the hours between 6 a.m. and 6 p.m. local time were assumed to be in nominal conditions, as can be observed in Figure 5(a). An AATR threshold of 0.25 TECU/min was used to distinguish nominal days/hours from days/hours with some ionospheric activity (see Figures 5(a) and 5(b)). The last step was to perform a visual inspection of the data recorded at the stations under study, depicted in Figure 4, during the days considered “nominal” and during the nominal hours of the “active” days. During this process, days with corrupted or missing measurements were discarded.

7.1 Monitoring Thresholds

To compute the detection thresholds using Equation (3), we selected the false-alert probability (P_fa) to be the power of 10 immediately below the number of samples in each elevation bin, considering the following bin sizes: 2° between 5° and 55° of elevation, 5° between 55° and 70° of elevation, and 10° between 70° and 90° of elevation. These different bin sizes were used to account for the fewer number of samples available from high-elevation satellites and for the faster movement of satellites at lower elevations in comparison to their slow movement at high elevation. Therefore, assuming that all samples are statistically independent, a P_fa of 10⁻⁴ was empirically selected as an acceptable compromise for the tests with real measurements from Brazil. Nevertheless, this value is likely not optimal because of the limited data on which its selection was based. Optimization of this value will be investigated as part of future work.

Figure 7 shows that all stations in the network have similar performance for satellite elevations above 13°. Below this elevation, several satellites crossed the Appleton anomaly at different locations, causing different nominal delays. This phenomenon occurred for satellite elevations between 7° and 13°, where the thresholds are higher than between 5° and 7°. In addition, amongst the compared stations, the threshold for station “uba1” is slightly lower for satellite elevations below 13°. This result occurred because fewer measurements were available at these elevations for the nominal hours of the active days and, therefore, the maximum values of nominal ionospheric delays could not be reached. This fact resulted in a more sensitive performance of station “uba1,” which could trigger more false alerts than the others for satellites at lower elevations.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

Detection thresholds for all stations considered in Brazil derived from the data specified in Section 5.2 for P_fa = 10⁻⁴.

7.2 Evaluation of the Monitor with Simulated Ionospheric Gradients

This section analyzes the detection and estimation capabilities of the network of stations depicted in Figure 4.

7.2.1 Detection

Figure 8 shows an example of $Ites{t}_r^{G04}\left(t\right)$ calculated for five of the stations considered, satellite G04, and a simulated gradient with the following characteristics: vertical slope of 200 mm/km (downward and upward), speed of 100 m/s, width of 40 km (corresponding to each downward and upward slope), zonal length of 90 km, and movement along a direction of 58° measured from the North Pole. Here, the two gradients forming the simulated EPB are clearly distinguishable.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 8

Example of simulation test statistics for satellite G04, day 145 of year 2014, and the network of stations in Brazil.

The algorithm was able to distinguish both parts of the simulated EPBs; the part corresponding to downward slopes presented negative $Ites{t}_r^j$ values, while the part corresponding to upward slopes presented positive $Ites{t}_r^j$ values. Therefore, the algorithm treated both parts as two different gradients to compute their parameters, following the methodology in Section 3 (for more details, refer to Section 4 of the work by Caamano et al. (2021)). Because the EPB simulations in this work follow the threat model proposed by Wang et al. (2021) and Yoon et al. (2017) for GBASs in the equatorial region, which proposes to assume a constant ionospheric delay in the middle of both fronts, the test statistic for this intermediate part was below the detection threshold. Therefore, this intermediate region was not considered for the estimation process in the simulations. However, when real EPBs are observed, there are often ionospheric irregularities between the two ionospheric fronts large enough to be detectable and cause hazards if undetected. Section 7.3 shows the limitations and adaptations of the estimation method used to account for these irregularities between the two fronts in the study with real (actually observed) EPBs. Future work will address the impact of these intermediate irregularities in simulation using more complex EPB models.

7.2.3 Uncertainty Model of the Estimated Gradient Slope

Figure 9(a) shows the slope estimation errors (blue stars) considering variation of all gradient parameters within the Brazilian threat model. Positive errors refer to cases in which the gradients were overestimated (estimates of the slant slope were larger than the simulated slant slopes) while negative errors refer to cases in which the gradients were underestimated (estimates of the slant slope were lower than the simulated slant slopes).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 9

(a) Slant slope estimation errors versus estimated slant slope (blue stars: slant slope estimation errors; red dashed line: mean value per bin; solid light green lines: standard deviation per bin); (b) number of samples in each cell of 10 mm/km of estimated gradient slope × 10 mm/km of error.

As described in Section 4.2, the first step in deriving the uncertainty model is to order the estimation errors in bins to calculate the statistics. We selected a bin size of 10 mm/km in slant slope as a compromise, allowing for a sufficient number of data samples per bin while ensuring that negative errors (i.e., gradient underestimations) remained relatively consistent within each bin. However, it should be noted that because the slopes were simulated in the vertical domain, the number of positive samples in the first bins of estimated slant slope (e.g., 50–60 mm/km) was lower than that in the remaining bins. This result occurred because only the high-elevation satellites were able to observe low slant slopes, as the applied obliquity factors caused low-elevation satellites to observe steeper slant slopes than the simulated vertical slopes. Therefore, the positive samples in the first bins (e.g., 50–60 mm/km of estimated slant slope) belonged to simulated slant slopes between 50 mm/km and 60 mm/km that were slightly overestimated (e.g., 55 mm/km), but not enough to be sorted in other bins. The negative samples in the first bins belonged to simulated slant slopes greater than 60 mm/km (e.g., 62 mm/km) that were underestimated (e.g., 55 mm/km). Thus, these samples were sorted into a bin of estimated slope that was different from the corresponding simulated slope. The maximum values of the negative samples were higher than those of the positive samples in the first bins because they belonged to higher simulated slant slopes that were underestimated and because a greater absolute value of error is expected for high simulated slopes.

In addition, Figure 9(a) shows the mean value per bin (red dashed line) and the standard deviation per bin (solid light green line). The mean values were positive but close to zero for all bins of estimated slant slope except the first bin (i.e., 50–60 mm/km of estimated slant slope), and the values of the standard deviation were also low.

Figure 9(b) shows the number of samples in each cell, defined by a size of 10 mm/km of estimated gradient slope × 10 mm/km of error. It can be observed that although some estimation errors are high, the number of samples in cells corresponding to these high errors is close to zero.

Figure 10 shows the estimated maximum value (blue stars along a straight line) that could occur inside each bin (e.g., 60 mm/km for the bin of 50–60 mm/km depicted in the center of the bin) and the overbounded value considering Equation (7) (black triangles). In addition, the figure shows the overbounded value if instead of assuming a probability of non-bounded errors of 1 ×10⁻⁸ (see Section 4.1), a value of 1 ×10⁻⁵ was assumed (green triangles), considering a prior probability of occurrence of an anomalous ionospheric gradient of 1 × 10⁻³, as in the work by Yoon, Lee et al. (2020) for the Conterminous United States (CONUS) and Yoon et al. (2019) for Brazil. Note that the prior probability of 1 × 10⁻³ was chosen to be conservative for CONUS, but as can be observed in Figure 5(a), it would not be conservative for Brazil at local nighttime. For an estimated gradient slope of up to 600 mm/km for the most conservative case (no prior probability considered) or 660 mm/km when assuming a prior probability, a benefit is obtained by using this methodology in Brazil. Above these values, the overbounded slant slopes would be higher than the Brazilian threat model, and in this case, the threat model would be used as-is instead.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 10

Worst-case overbounded gradient slant slopes in comparison to the Brazilian threat model (solid red line) (blue stars along a straight line: maximum possible estimated slant slope inside each data bin; black triangles: worst-case overbounded slant slope in each bin; green triangles: overbounded slant slopes in each bin when a prior probability of occurrence of an ionospheric gradient of 1 × 10⁻³ is assumed).

7.3 Evaluation of the Monitor with a Real Ionospheric Gradient

This section evaluates the methodology with a real anomalous ionospheric gradient measured by the network in Brazil. Unlike the simulated gradients, the real gradient was accompanied by many cycle slips that made parameter estimation more difficult. However, it should be noted that in a real implementation, receivers with a higher bandwidth, which are more robust against scintillation, could be installed. Furthermore, these receivers should have measurements available on the L5 frequency, which is more robust than L2 for this problem. To detect cycle slips, we used the cycle slip detector described in Section 4.3.1.1 of the work by Sanz et al. (2013), using a constant threshold of 4 m. In addition, after the application of the cycle slip detector, the cycle slips were visually inspected to avoid discarding excessive EPB measurements when the ionosphere variations were large but not due to a cycle slip.

The network was able to detect all affected satellites. Furthermore, as will be demonstrated in the following sections, the algorithm could estimate the selected real gradient, both in post-processing and real time, by using the measurements from satellite G31 and the stations in the network under study on day 68 (local time) and 69 (GPS time) of year 2014, i.e., during the peak of Solar Cycle 24. The post-processed results were used for comparison with those obtained in real time to evaluate the feasibility of the real-time concept.

7.3.1 Detection

Figure 11(a) shows slant ionospheric delays corresponding to the epochs when the test statistics exceeded the threshold for each station. During the local times depicted, satellite G31 had an elevation between 33° and 53°. As can be seen, there is a vertical offset in the estimated slant ionospheric delays, due to the ambiguities present in the carrier-phase measurements in Equation (1).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 11

(a) Slant ionospheric delays corresponding to the epochs when the test statistics exceeded the threshold for each station; (b) absolute value of the test statistics, zoomed in to highlight the main depletion for the studied gradient experienced by satellite G31 (on day 68 local time and day 69 GPS time of 2014 in Brazil).

The solid orange line in Figure 11(b) represents the detection threshold for all stations, as they all exhibited similar values during this period.

Figure 11(b) shows the absolute value of the test statistics for the same satellite, $Ites{t}_r^{G31}\left(t\right)$ , zoomed in to highlight the main depletion for the studied gradient. Here, the previously mentioned vertical offset is not present because the test statistics are computed as the first derivative of the slant ionospheric delays with respect to time. The threshold for all stations is represented by a solid orange line for simplicity, as the thresholds for all stations exhibited very similar values during this period.

As can be observed, the algorithm achieved nearly continuous gradient detection for all stations. Note that the main EPB was accompanied by earlier or later minor EPBs (e.g., station “chpi”) or by intermediate periods between the two main EPB slopes with many oscillations (e.g., station “chpi,” where the upward part was not observed). For this reason, in highly active periods, where the EPBs come in sequences and are accompanied by excessive scintillation, the network would be constantly alerting and the EPB parameters would be especially difficult to estimate. In these cases, the network would indicate the use of the “worst-case” threat model for all times when the gradient parameters could not be estimated. To avoid on-and-off alerting under these conditions, a minimum alert period could be implemented (e.g., 2 min). During this period, network alerts would be temporarily deactivated, and the conservative parameters of the current threat model would be applied. Normal network operation would only resume after this period, ensuring stability before alerts are reactivated. A similar approach has been proposed by Marini-Pereira et al. (2024).

7.3.2 Estimation of Gradient Parameters

Figure 11(a) shows that the slant ionospheric delays presented a similar shape when being measured by all affected stations and the same satellite. The figure also suggests that the gradient changed slightly during its propagation; the drop in ionospheric delay corresponding to the depletion was smaller when the gradient reached the last stations impacted by it. However, the slant slope for each depletion and the test statistics were more similar than the drop in meters. Therefore, as long as the cross-correlation coefficients computed with the test statistics for the same satellite and different stations were greater than 0.5, we considered the gradient to be the same for these stations. As for the estimation with simulated gradients, real data from stations “poli,” “sjsp,” “mgin,” “uba1,” and “chpi” were used for the estimation of gradient parameters with the real gradient. The algorithm was able to estimate the gradient parameters both in post-processing and in real time for satellite G31.

Estimation of Gradient Parameters in Post-Processing with Satellite G31

The first step of the algorithm searched for the reference station. All stations presented high cross-correlations (above 0.9) with each other, meaning that the gradient did not change substantially while propagating from the viewpoint of the test statistic, as mentioned previously. Therefore, the first station that detected the gradient, “poli,” was selected as the reference station.

Initially, the gradient parameters observed by satellite G31 were calculated in post-processing by using a spline interpolation of the data at 10 Hz to obtain the best estimation possible and to evaluate whether there are any problems with the time resolution of the adopted data set. According to the post-processing results, the measured ionospheric gradient was an EPB that traveled with an estimated direction of 44° (from southwest to northeast) and a low–medium speed, 165.8 m/s. The slant slopes estimated by each network station were as follows: –481.6 mm/ km for “poli,” –507.6 mm/km for “sjsp,” –543.7 mm/km for “mgin,” –379.1 mm/km for “uba1,” and –474.6 mm/km for “chpi.” As can be observed in Figure 11(a), only the results corresponding to the downward slopes were available owing to cycle slips affecting the upward slopes, and therefore, the sign of the slope estimates was negative. The widths estimated by each network station were as follows: 59.8 km for “poli,” 41.8 km for “sjsp,” 38.2 km for “mgin,” 35.8 km for “uba1,” and 33.5 km for “chpi.” These results are compatible with the Brazilian threat model introduced in Section 2.

Estimation of Gradient Parameters in Real Time with Satellite G31

The application of the algorithm in real time with data recorded at a 15-s sampling rate presents the same limitations as in the simulation: it is necessary to have sufficient time for the cross-correlation estimate to converge, sufficient temporal resolution, adequate distance between stations to find a sufficient correlation between the test statistics of gradients that are the same, and test statistics that are not excessively noisy between the two slopes (downward and upward) and during the occurrence of the gradients.

Figure 12(a) shows results regarding the first limitation, where the convergence requirement (the difference in cross-correlation coefficients at the current epoch and previous epochs is less than 1 × 10⁻² for the last three samples, as in Equation (10) of the work by Caamano et al. (2021)) was not met until 21:54:14 for station “sjsp,” 21:56:44 for station “uba1,” and 22:03:29 for station “chpi” in local time. Station “mgin” did not meet the convergence criteria and thus was not used to determine the real-time gradient propagation parameters. Furthermore, Figure 12(a) shows that even once the algorithm converged, the cross-correlation coefficients did not remain constant. Instead, the cross-correlation coefficients experienced several jumps, and their value decreased as more samples from the test statistics were considered in their calculation. This result occurred because at the beginning of the detection at each station, the samples that were considered for the cross-correlation calculation corresponded to the first of the gradient slopes (downward). As this slope is similar at almost all stations, the algorithm converged to high cross-correlation coefficients and time delays that were comparable to the post-processed results.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 12

(a) Cross-correlation coefficients and (b) time delays calculated in real time for the test statistics depicted in Figure 11(b), considering “poli” as the reference station.

Once the first main slope had passed, the test statistic samples belonging to the intermediate period between slopes (the EPB depletion part) were considered. In the simulation, we assumed that the test statistics would quickly decrease below the threshold during the depletion part, allowing for a clear separation of the downward and upward slopes to be considered as two different gradients. However, when real measurements were used, the test statistics oscillated and experienced numerous cycle slips, causing the test statistics to be considerably different from station to station and, therefore, triggering the thresholds. In these cases, the separation in two gradients would be produced by a cycle slip, by a sufficient number of test statistic samples below the corresponding threshold, or by the cross-correlation coefficients finally decreasing below 0.5. Figure 12(b) shows that although the value of the cross-correlation coefficients changed, the time delays remained practically constant during the time period under study until the signal was lost or there was a cycle slip.

Figure 13 shows the real-time estimates of the speed and direction parameters before the convergence criteria of the algorithm were met (black dots) and the estimates calculated by the algorithm after they were met (red dots). As can be observed, the estimation errors were initially high for both the speed and direction, when the algorithm did not converge and the cross-correlation coefficients were low. However, the algorithm consistently estimated the speed and direction parameters after convergence. The estimated values changed smoothly when the same set of stations was used for their estimation, and there were no changes in the estimated time delays. However, this was not the case when the set of stations used changed or there was a variation in the time delay estimates. The two jumps that can be observed after the algorithm first converged in Figure 13 between 21:58:04 and 22:06:24 local time were due to a small variation in the time delay estimate for station “sjsp” at 22:00:04 local time and the use of stations “poli,” “uba1,” and “chpi” instead of “poli,” “sjsp,” and “uba1” for the estimation process from 22:02:44 local time. The set of stations used was different because, after 22:02:44 local time, the data from station “sjsp” had a gap during some epochs and the algorithm concluded that the gradient had ended for this station. At the same time, the algorithm converged for station “chpi,” and thus, this station could be used in the estimation process. It can also be seen that the algorithm correctly estimated the propagation parameters at some points, but these points are not highlighted in red because the convergence criteria were not met. This result was due to the aforementioned jumps in the cross-correlation coefficients that did not affect the estimation of time delays between stations. This finding indicates that the convergence criteria used might be overly conservative for the case of Brazil. Nevertheless, optimization of the convergence criteria for different data sets was not addressed in this paper and is part of future work.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 13

(a) Estimated speed and (b) direction of the real gradient depicted in Figure 11(a) (red points: estimations that were calculated after the algorithm converged; black points: estimations that were calculated before convergence).

The real-time speed estimates varied slightly between 114 m/s and 155 m/s, while the direction estimates varied between 37.5° and 48.5°. These values, especially the last points in Figure 13 when the stations used were “poli,” “uba1,” and “chpi,” 143.6 m/s and 41.3°, are similar to the 165-m/s speed and 44° direction obtained in post-processing with the data interpolated at 10 Hz, indicating that the final real-time results converged to the post-processed values when using three of the five stations that the algorithm considered to obtain the post-processed results.

The variations in the above parameters can be explained as follows. Firstly, because the stations used are separated by 100 km or more, it is possible that the EPB varied slightly in its propagation and shape, with the different stations observing slightly different parts of the front. In addition, the elevation of satellite G31 also changed from the time of detection at “poli” to the time of detection at “chpi,” which had an impact on the IPP speed and, therefore, on the estimation of the gradient propagation parameters. Secondly, the algorithm used two different sets of stations for the estimation, as discussed previously. Because the center of the local coordinate system for estimating the parameters can vary depending on the position of the IPPs at the time of detection (see Section 4.2.2 of the work by Caamano et al. (2021)), the center was different for calculations with different sets of stations. Moreover, as the IPPs moved farther from the center of the local reference coordinate system, the approximation errors at each location increased with increasing distance between the IPPs and the center of the local reference coordinate system. Thirdly, because only three of the five available stations were used for real-time estimation, there was no redundancy, as the minimum number of stations for the estimation process is three. Therefore, any errors (e.g., due to an insufficient sampling rate) in the estimation process had a greater impact in this case than when a larger number of stations was used for estimation.

Finally, the real-time algorithm estimated the slant slope and width parameters. Because the algorithm was able to estimate the gradient speed and direction after the maximum absolute value of $Ites{t}_r^{G31}\left(t\right)$ was reached for stations “poli,” “sjsp,” “mgin,” and “uba1,” it recalculated the gradient slope backwards to find the “worst case” that could be transmitted to the GBAS stations. To this end, the algorithm used the first speed and direction estimates after convergence. If real-time estimates of speed and direction were available while the slope was being determined, the algorithm used these real-time estimates to calculate the slope at those times. Figure 14 shows the estimated slant slope obtained through this approach. Here, the slant slope estimates were slightly lower but close to those calculated in post-processing because the relative speed between the gradient and the IPPs was slightly higher. Thus, the same ionospheric delay rates resulted in lower slopes. Note that the gradient slopes are negative when the gradient is downward and positive when it is upward, unless absolute values are taken.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 14

Estimated slant slope of the gradient depicted in Figure 11(a) using speed estimation after convergence.

The results from the estimated width parameter, which was determined after the gradient had affected all stations, were 62.7 km for “poli,” 43.4 km for “sjsp,” 40.7 km for “mgin,” 37.8 km for “uba1,” and 35.1 km for “chpi.” In this case, the width was larger than in the post-processed results because the relative speed estimate between the gradient and the IPPs was larger. Note that the determination of the width was based only on the part of the test statistics that was continuously above the threshold, which corresponded to the downward slope. For this reason, the algorithm estimated a single width per station and satellite instead of two separate widths corresponding to the downward and upward slopes.

7.3.3 Overbound of the Estimated Slant Slope in Real Time in Brazil

The gradient parameters estimated for this satellite were within the Brazilian threat model, and thus, the uncertainty model for the slant slope estimation determined in Section 7.2.3 was valid. Figure 15 shows the overbounded gradient slant slopes for satellite G31 and each station for which estimation was possible. The overbounding method was not applied to the values between ± 50 mm/km of estimated slant slopes in Figure 14, as these gradients are considered to be nominal. As can be observed, the worst case for the estimated and overbounded gradient slant slopes for satellite G31 were -662.62 mm/km for “poli,” -677.89 mm/km for “sjsp,” -727.71 mm/km for “mgin,” -561.44 mm/km for “uba1,” and -563.22 mm/km for “chpi.” These values were calculated conservatively, considering a probability of non-bounded errors of 1 × 10⁻⁸ and, therefore, k_ne = 5.61 in Equation (7). It can be observed that even after overbounding has been applied, the resulting gradient slopes are still below the upper bound of the Brazilian threat model (orange line in Figure 15); therefore, increased availability and continuity can be expected when using this method. An evaluation of how this method improves availability and continuity will be addressed in future work.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 15

Overbounded slant slopes of the gradient for each station.