Research ArticleOriginal Article
Open Access

Integrity risk evaluation of impact of ionospheric anomalies on GAST D GBAS

Moonseok Yoon, Jiyun Lee and Sam Pullen
NAVIGATION: Journal of the Institute of Navigation June 2020, 67 (2) 223-234; DOI: https://doi.org/10.1002/navi.339

Abstract

This study develops a three-step Monte Carlo method to evaluate the worst possible integrity risk of ionospheric spatial gradients for GAST D GBAS. Impact simulation parameters are classified into two groups, “worst-case” and “average,” based on the underlying integrity requirements. Unlike “average” parameters, “worst-case” parameters are those for which a clear basis for averaging could not be established because the probabilistic distribution of these parameters cannot be developed with sufficient confidence due to the lack of observation data. In calculating integrity risk, these “worst-case” parameters use the worst-case value that maximizes the integrity risk. Each step of the randomized search narrows down the parameter ranges in sequence and identifies the two worst-case parameter sets (based on the largest position error and the largest missed-detection probability) for each worst-case region identified in the initial step. The resulting integrity risk values are well below 10−9, showing that the GAST D SARPs integrity requirement is met.

1 INTRODUCTION

Ground-Based Augmentation System (GBAS) Approach Service Type D (GAST D) is designed to support Category (CAT) II/III precision approach operations for single frequency GBAS users. One notable change from GAST C is that responsibility for mitigating ionospheric gradients is shared between airborne users as well as the ground system. As a result, several monitoring and geometry screening techniques, including airborne monitors, have been implemented.1 In particular, ground-based monitoring of ionospheric spatial gradients using measurements across reference receiver baselines (known as Ionospheric Gradient Monitoring or IGM) has been added to detect severe gradients affecting the ground system that may not be detectable in GAST C (eg, because they are slow-moving).2-6 Screening to remove potentially unsafe satellite geometries (in the presence of the worst undetectable ionospheric gradients) has been moved from ground to the airborne system, which reduces conservatism because the airborne receiver knows the geometry it is applying for positioning.7 Code-carrier divergence (CCD) monitoring similar to that utilized in the GAST C ground system is also applied in the airborne subsystem to detect gradients before they reach the ground system.8,9

A new requirement to achieve a missed detection probability (PMD) of 10−9 or lower for 30-s differential range errors exceeding 2.75 m at the Landing Threshold Point (LTP) due to anomalous ionospheric gradients was codified in the GAST D SARPs (Standards and Recommended Practices).10 This requirement was previously validated in Kim et al.11 using a Monte Carlo simulation approach that enables probabilistic analysis of the effects of anomalous ionospheric gradients on GAST D GBAS approach and landing operations. In this approach, all simulation parameters, which form each ionospheric gradient impact scenario, were randomly selected from their selected probabilistic distributions using Monte Carlo sampling. Then, the differential range error at the LTP and the probabilities of missed detections (PMD) from GBAS ground and airborne monitors were computed as part of the simulation of each sampled scenario. These monitors include the ground CCD monitor, the airborne range-domain dual-smoothing ionosphere gradient monitor algorithm (DSIGMA) monitor, and the ground ionosphere gradient monitor (IGM).12 More than 109 sampled events were then combined to compute the probability of Hazardously Misleading Information (PHMI), which is defined as the sum of PMDs of all events (in which the differential error at the LTP exceeds 2.75 m) divided by the total number of sampled scenarios.

This method is useful for verifying that the integrity requirement is met in a probabilistic sense using the calculated PHMI values. One concern raised by this approach, however, is the lack of sufficient independent observation data used to model ionospheric anomalies that results in insufficient confidence in the probability distributions of the simulation parameters used for sampling. This can lead to improper averaging when calculating the PHMI values. To address that concern, this paper defines a rationale to classify the simulation parameters into two groups: “worst-case (WC)” and “average.” Based on those two definitions, we have developed a three-step worst-case parameter search and integrity risk evaluation method. In each step, this method refines the search for the worst-case set of threat parameters and narrows down the threat parameter ranges in sequence. By doing this, single-valued worst-case sets can finally be identified from different parameter ranges. As a result, we are able to calculate the worst possible PHMI values.

This paper first reviews a detailed description of the ionospheric anomaly threat model parameters and geometric parameters that represent the anomalous ionospheric gradient impact scenario for GBAS-supported precision approaches. It then details the rationale and the selection of WC or average for specific threat parameters. Detailed algorithms of the three-step worst-case parameter search and integrity risk evaluation are presented in the following section. This paper also analyzes the simulation results from the three-step method. The final section reports the conclusions of this paper.

2 REVIEW OF IONOSPHERIC GRADIENT IMPACT SCENARIO AND PARAMETERS

In this section, we review the ionospheric anomaly threat model and geometric parameters (previously defined in Belabbas and Meurer4), which were used to analyze the effects of anomalous ionospheric gradients on users of GAST D GBAS. Figure 1 shows an overhead view of a GAST D ionospheric threat scenario where a hypothetical aircraft is impacted by an anomalous ionospheric wave front. The ground facility is assumed to be located LLTP (5 km in this paper) away from the LTP with an orientation angle α. With the predetermined speed profiles defined in Working Paper 3,13 the aircraft is configured to run into the moving wedge-shape ionospheric wave front model in a relative approach direction to the front of β (which is calculated using the front moving direction fdir and the aircraft approach direction acdir). The predetermined speed profiles as a function of remaining time before reaching LTP (randomly chosen in the simulation) are shown in Figure 2. The profiles that are labeled “fast” (solid line), “moderate” (dashed line), and “slow” (dashed dotted line) are characterized with landing ground speed, time at landing speed, deceleration rate, and ground speed at start of deceleration. The threat model captures the mid-latitude ionospheric storms observed in the conterminous U.S. (CONUS).14 The slope of the front (g) represents a linear change in ionospheric delay from high to low. Note that the product of the gradient (g) and the width of the ionospheric front (W) yield complete delay depletion, which is also constrained by the threat model. The ranges of the anomalous gradient model parameters are summarized in Table 1 in the next section.

FIGURE 1
FIGURE 1

Anomalous ionosphere impact scenario parameters for Monte Carlo simulations

FIGURE 2
FIGURE 2

Aircraft speed profiles

View this table:
TABLE 1

Ionospheric and geometric scenario parameters: WC vs. Average classification and rationale

The detection capability of the GAST D monitors mainly relies on the speed (ΔvΔv) of the ionospheric front relative to the speed of a satellite ionospheric pierce point (IPP) projected onto a line perpendicular to the front. Thus, the simulation assumes that the front moves with a constant speed of Δv, which is calculated from the front speed (vfrvfr), front moving direction (fdir), and IPP geometry (based on satellite geometry). The parameter γ is defined as a relative angle between the LTP orientation angle α and the front approach angle relative to the aircraft approach direction β. For a given set of simulation parameters, the position (ΔW) of the ionospheric front relative to the aircraft at the end of the simulation is selected. The simulation period (tsim) and the corresponding initial position of the aircraft are then determined by rewinding it back until both the aircraft and the ground facility are outside the ionospheric front, at the beginning of the simulation. In some cases, tsim could be abnormally large when Δv is very small or close to the approaching speed of the aircraft. Thus, tsim tsim is constrained to be less than 40 minutes.

Monte Carlo sampling is employed (see Kim et al.11) to select the threat model and aircraft parameters for a single simulated aircraft approach. Once a given set of the threat model and the aircraft parameters is sampled, an approach simulation using those parameters starts with the aircraft moving toward the LTP of a runway and ends when the aircraft reaches the evaluation point at the LTP. Differential range errors and PMD from GBAS ground and airborne monitors at each epoch are produced under an “in-between” scenario where the aircraft remains within the ionospheric front at the end of the simulation (ie, when the aircraft reaches the LTP).

3 RATIONALE FOR “WORST-CASE” VS. “AVERAGE” PARAMETERS

A combination of ionospheric threat model parameters (eg, gradient, width, and speed) and geometric parameters from Figure 1 influences the differential error and PMD obtained at the LTP. A key question for a Monte Carlo simulation that uses these parameters as inputs is how to regard the distributions of each parameter when the final probability of hazardously misleading information (ie, PHMI—essentially, the probability of violating the GAST D integrity requirements) is calculated. Parameters treated as WC must use the worst-case value of that parameter in calculating PHMI, meaning the value of the parameter that maximizes PHMI, even though the WC value typically has a low probability of occurring (even given that an anomaly occurs). On the other hand, parameters treated as average can use a sampled distribution (usually from a uniform distribution) between the minimum and maximum values of that parameter as the basis for calculating PHMI. In other words, for average parameters, the effect of each sampled value of that parameter is evaluated in the simulation and given a weight corresponding to its sampled probability. This is much less conservative than assuming that the WC value (which usually must be found via search) always occurs.

Table 1 identifies each of these parameters and shows whether they are treated as WC or average in the calculation of PHMI. This table gives the minimum and maximum ranges for each parameter, its status as WC or average, and a brief rationale for why WC or average was chosen. Parameters chosen to be WC are those for which a clear basis for averaging could not be established. For example, the gradient magnitude (g) is chosen to be WC because few ionospheric storms with gradients in the range of concern (above 200 mm/km) have been observed, and the resulting data is insufficient to establish the probabilistic distribution of these gradients (ie, what would be needed to implement averaging) with sufficient confidence. The data that is available suggests that, as expected, the probability of observing an anomalous gradient of a given magnitude decreases as the magnitude increases, which, in turn, suggests that averaging based upon a uniform distribution of gradients between 200 and 500 mm/km would be acceptable. However, it was decided that more independent observations of severe ionospheric storms (as opposed to multiple observations of the same storm from different ground stations and look angles) would be needed to validate this or any other averaging approach. Similar rationales apply to several other key ionospheric front parameters (velocity vfr and width W).

Parameters chosen to be average are those which are known to be essentially random from the point of view of a GBAS precision approach impacted by an ionospheric front. The parameter ΔW, for example, is based on the relative timing of the approaching aircraft and ionospheric front, which is variable and random for each simulated aircraft approach. Similarly, the front direction (of approach) parameter fdir is also treated as random for mid-latitude locations where there is no reason to expect that one approach direction is more or less likely than another. Note that, in Figure FIGURE 1, the angle β between the front direction of approach fdir and the aircraft approach direction are directly related, so it is equivalent to treat either fdir or β as the random parameter. Also note that low-latitude regions affected by equatorial plasma bubbles at nighttime observe the resulting gradients mostly in a west-to-east direction, so the “uniform likelihood of direction” model applied here (for mid-latitudes) would not apply there (at night).

4 THREE-STEP METHOD FOR WORST-CASE PARAMETER SEARCH AND INTEGRITY RISK EVALUATION

4.1 Overview of three-step method

The three-step worst-case parameter and integrity risk evaluation method shown in Figure 3 searches for and identifies the worst-case sets of ionospheric threat parameters that are designated as WC and the impact of averaging over the parameters that can be averaged. Note that Ntot, Nave, and NWC denote the total number of samples, the number of average variables, and the number of WC variables, respectively. Ntot can be obtained by multiplying Nave with NWC (ie, Ntot = Nave × NWC). The first step performs a randomized search over all of the parameters in Table 1 with a very large total number of samples (ideally using a sample size of 109). The goal of this step (Step 1) is to identify one or more regions of the threat space (combinations of threat parameters) where the worst-case result (in terms of the maximum PHMI) may lie.

FIGURE 3
FIGURE 3

Three-step worst-case parameter search and integrity risk evaluation

Once one or more (perhaps as many as five but more commonly two or three) potential worst-case regions have been identified, each identified region is re-simulated in Step 2 with much narrower parameter ranges designed to focus the sampled parameters on the region in question. The goal of Step 2 is to refine the search for the worst-case set of threat parameters so that a single worst-case set of parameters can be selected for each region. This is typically a manual process as the results of the first Step simulation for a given region may suggest altering and further tightening the parameter ranges to search more closely in a follow-up Step 2 simulation. Note that the set of threat-model parameters referred to here only includes those that must be set at worst-case values according to Table 1 as the other two averaged parameters (ΔW and β) are always allowed to vary.

Once each identified region has a single set of worst-case threat-model parameters, separate Step 3 simulations can be conducted for each region to estimate PHMI, which is what is compared with the anomalous ionosphere requirements given in ICAO State Letter.10 In each Step 3 simulation, the worst-case set of threat-model parameters is set as fixed and is not varied; thus, no sampling within worst-case parameter ranges occurs. Only the averaged parameters (two in this case) are varied, and the results should represent the worst possible performance within each region identified after Step 1. For each Step 3 result, PHMI is computed as follows:

Embedded Image 1

where Pprior represents the prior probability of the anomalous ionospheric events sampled in the Monte Carlo simulation, i ∈ φ represents the subset of all sampled events that have differential errors at the LTP that exceed a certain critical error value E, which is defined to be 2.75 m by the GAST D requirements, and PC,MD is the combined PMD over the three independent monitors simulated for each gradient event. This paper allows for a prior probability of 10−3 against ionospheric gradients in mid-latitudes, and thus Pprior is set to be 10−3. (The derivation of Pprior will be given in the next subsection.)

Equation (1) sums the PC,MD values of all events in which the differential error at the LTP exceeds 2.75 m, and this sum of probabilities is divided by the total number of sampled scenarios, which can be averaged over. In other words, Equation (1) collects all combinations of worst-case and averaged (sampled) parameters for which the simulated aircraft approach has a differential range error exceeding 2.75 m at the LTP, sums the PMDs of these cases, and then divides by the total number of simulated scenarios (ie, Ntot). Note that Ntot includes all simulated scenarios, not only those that generated differential range errors exceeding 2.75 m at the LTP. The prior probability of being in the simulated anomalous ionospheric state (Pprior, which was determined to be 10−3 for mid-latitude regions in the next subsection) is then included. Each resulting PHMI (for each region) simulated must be below 10−9 per worst-case event for the requirements to be met. If Pprior = 10−3 is separated from the calculation of PHMI, as when computing an overall PMD for a particular set of worst-case parameters, the requirement becomes 10−6 instead.

4.2 GAST D prior probability of threatening anomalies

A prior probability of 10−3 against anomalous events sampled by Monte Carlo simulation can be justified based on previous CONUS data collected from 1999 from 2005. First, we assume that potentially anomalous (stormy) mid-latitude ionospheric conditions are limited to periods with Kp indices equal to or larger than 5, while anomalies on periods with Kp indices less than 5 are extremely unlikely to create threatening gradients for GBAS. From the record of Kp indices (at 3-h intervals), the equivalent number of days with Kp index equal to or larger than 5 from 1999 to 2005 was calculated to 197.75. From the analysis in Kim et al.,11 gradients below 200 mm/km were found to be unable to create potentially hazardous errors.

Using these criteria, we found four threatening days (with periods of gradients above 200 mm/km) out of 191.75 observed days from 1999 to 2005. These were all caused by the extreme solar coronal mass ejection (CME) in October 2003: 10/29/2003, 10/30/2003, 10/31/2003, and 11/20/2003. One further assumption made from CONUS ionospheric observations was that no more than 1 hour of a threatening day would experience gradients larger than 200 mm/km (on any or all satellites) as observed by a single site. The resulting prior probability is calculated as follows:

Embedded Image 2

We introduced further conservatism to this value by rounding up to 0.001 or 10−3 because it is not possible to prove from the measurements available to us that only four threatening days occurred during this period.

While several assumptions were needed to produce the prior probability in Equation (2), we believe that these assumptions are valid based on our years of studying mid-latitude ionospheric gradients and that the result is conservative when applied to mid-latitude regions. First, although only data from CONUS was used, data from other mid-latitude regions in Europe,15,16 Japan,17 and Australia18 show that gradients above 200 mm/km appear to be even rarer in these places, except for ionospheric pierce points that fall outside mid-latitude regions. For example, the largest spatial gradient in the Australian study18 was only 148 mm/km. In the largest European study conducted to date,15 gradients as high as 403 mm/km were observed in equatorial locations (the Canary Islands), but within the mid-latitudes, only one low-elevation gradient reaching 343 mm/km was observed (from Sweden)—the rest were all below 200 mm/km.

Second, in computing the probability of hazardous gradients from the CONUS data, the four days with such gradients were treated as separate instances, even though they were caused by the same physical event, as explained above. However, the CME that caused these extreme gradients did not do so uniformly. Gradients above 200 mm/km were only observed over specific localized regions on specific days and times, while most of the mid-latitude regions were not so severely affected. Given this, it seemed reasonable to treat these occurrences as separate, although they are only likely to occur following a very strong solar-magnetic disturbance.

Finally, the assumption that gradients above 200 mm/km would persist at a single location for no more than one hour is based on results showing the passage of ionospheric fronts (as described by the GBAS ionospheric threat model explained above) over fixed reference receivers (eg, see Datta-Barua et al.,19). In most cases, satellites are affected by such large gradients for less than 15 to 30 minutes before the front passes by. A conservative value of 60 minutes (1 hour) was used in Equation (2) to also cover “slow-moving” ionospheric front scenarios (ie, small magnitude of Δv due to a small difference between front velocity and satellite IPP velocity with respect to the ground).

4.3 Nested loop approach to implement three-step method

In order to implement the proposed method, a nested loop approach is employed to allow averaging only over relative front approach direction (β) and aircraft/front phasing (ΔW) instead of treating all parameters as random variables. As shown in Figure 4, the Monte Carlo sampling simulation that exists outside of each approach simulation has a loop structure to sample satellite geometry, time epoch, threat model, and aircraft approach parameters. The outermost loop of the simulation updates GPS satellite geometries (Nsat satellites per geometry) every minute, which in turn yields 1440 epochs (Nepochs) in a day. The middle loop performs Nave samples of the average variables in Table 1, which are β and ΔW. Note that Nave includes Nsat and Nepochs since the time-varying satellite geometry (or position) is treated as an average parameter. The inner loop performs NWC samples of the WC variables in Table 1, which are the ionospheric front gradient (g), front width (W), front speed (vfr), front direction (fdir), aircraft speed (vac), and LTP orientation angle (α). The parameters are randomly generated from the uniform distribution of threat model values given in Table 1. The total number of samples (Ntot) is set to be approximately 109 (ie, Ntot = Nave × NWC).

FIGURE 4
FIGURE 4

Nested loop approach for three-step method

As noted above, the simulation of each aircraft approach includes GBAS ground and airborne monitors, which operate in the presence of the sampled anomalous ionospheric gradient parameters to fully analyze the threat that they pose. These monitors include the ground CCD monitor, the airborne range-domain DSIGMA monitor, and the ground IGM. The ground CCD monitor (which is defined as a series of two cascading first-order low-pass filters) currently implemented in CAT I GBAS detects ionospheric gradients by observing CCD rates.20 A CCD monitor is also present in the airborne system,21 but the airborne DSIGMA monitor is very highly correlated with it and detects more reliably in every case. Therefore, separate modeling of the airborne CCD monitor is not necessary and would not change the results. DSIGMA monitors the discrepancy between two carrier-smoothed code measurements with different time constants (100 and 30 seconds) to detect a smoothing filter lag difference caused by ionospheric gradients.1,21 IGM directly measures ionospheric gradients by utilizing double-difference carrier-phase measurements observed from IGM receivers.2-6 In this paper, we assume that four IGM receivers (which usually also serve as GBAS reference stations) are located at vertices of a square constructed around the centroid of the GBAS ground facility. From each pair of IGM receivers, the magnitude of the gradient is measured based on observed ionospheric delays. The maximum value of the gradient measured over the receiver pairs is taken as the IGM test statistic.

More details regarding the implementation of these monitors are described in Kim et al.11 The monitor performance characteristics assumed in the simulations conducted in this paper are the ones given in Kim et al.,11 except that several parameters were modified to better represent specific GAST D GBAS ground systems now under development. The exact values of these parameters affect the results somewhat, but they do not affect the simulation approach or methodology described in this paper.

For each of these three monitors, missed detection probabilities (PMDs) are computed at every second to search for the minimum PMD values for each monitor over the approach simulation time (tsim). The three minimum PMDs for each monitor are then combined as follows:

Embedded Image 3

The combined PMD (or PC,MD) is used to derive the value of PHMI. The differential range error at the LTP and the value of PC,MD are stored for evaluation from each simulation run. After all loops and samples have completed, the results of all Ntot scenarios are combined to estimate the overall PHMI (representing an unsafe condition) using Equation (1).

5 RESULTS FROM THREE-STEP METHOD AND INTEGRITY RISK EVALUATION

The standard RTCA 24-satellite GPS constellation (see RTCA SC-1598) and the GBAS ground facility at Newark Airport (Newark, New Jersey, near New York City) were chosen for the simulation. The simulation also included a re-admittance constraint, which is related to the recovery process in which GBAS monitors allow excluded satellites to be included again after passing tightened thresholds during (or following) a defined exclusion period. For each monitor, re-admittance was permitted only if the test statistic remained below twice the bounding test statistic sigma for 60 consecutive seconds after monitor initialization.

The Monte Carlo simulation using the parameter ranges shown in Table 1 represents only the first step among the three steps of this method as an Ntot of about 109 samples is not sufficient to examine multiple worst-case regions or families. Thus, the first step results are used here to identify two worst-case regions where events with PC,MD greater than 10−7 exist for Step 2 simulation. This constraint (PC,MD ≥ 10−7) was chosen because 10−7 represents 10% of the derived requirement on overall PMD (10−6) once the prior probability of 10−3 is taken into account. Events with PC,MD above 10−7 suggest that nearby (unsampled) events with PC,MD exceeding 10−6 might exist and should be investigated in Step 2.

The Step 2 Monte Carlo simulations use the same Nave and NWC values but greatly narrow the ranges of worst-case parameter values to focus on each region of interest. Then the single-valued worst-case parameter vectors are selected from the Step 2 results to run in Step 3 Monte Carlo simulations where NWC is equal to one and Ntot (or Nave) is about 109. In Step 3, PHMI is computed using Ntot.

Figure 5 shows a scatter plot of the resulting 30-s differential range errors (or 30-s smoothed corrected pseudorange errors at the LTP; along the y-axis) at the LTP as a function of the speed of ionospheric front relative to the IPP observed from by ground facility (along the x-axis) and the ionospheric gradient magnitude (gray scale). Again, events with a combined PC,MD less than 10−7 are not shown in this plot. Three different peaks (above 2.6 m) denoted by Regions 1 and 2 and Region 2 (representing gradients around 330 and 500 mm/km, respectively) are shown in Figure 5. The first peak in Region 1 is mainly due to the limited observability of IGM (ie, when IGM observes gradients below about 300 mm/km, its PMD is higher than 10−9). The second peak of range errors (in Region 2), whose gradients exceed 400 mm/km, is shown in a range of relative speed between 0 and +30 m/s. In this region, the anomalous ionospheric gradients are difficult to detect by the ground monitors until the aircraft reaches the LTP, and thus only the airborne monitor is responsible for detecting anomalous gradients but could not observe them in all of these cases. The third peak (Region 3) between +70 and +170 m/s occurred due to the weakness of the airborne DSIGMA monitor at these speeds because its detection capability, which relies on the CCD rate over time, is reduced at these speeds. The three peaks were defined as Region 1, Region 2, and Region 3, respectively, and we used them to specify the worst-case regions for Step 2.

FIGURE 5
FIGURE 5

Step 1 simulation result: Selecting worst-case regions

We selected Region 2 first for integrity risk evaluation since the threat scenarios in this region are considered to be the most hazardous to GBAS users. Among the threat points in Region 2, a common set of results with differential error (E) greater than 2.6 m and PC,MD larger than 10−7 were selected, and the narrow ranges of threat parameters of the points in this common set were then specified based on the largest error and the largest PC,MD. Note that a common set is defined as a set of ionospheric scenarios in which the resulting differential error (E) is larger than 2.6 m and PC,MD is larger than 10−7. The ionospheric model and geometric parameters of the selected scenarios are represented in Figure 6 and the resulting ranges of parameters are shown in Table 2. In Figure 6, the ionospheric impact scenarios selected with the largest error are denoted by a rectangle, and those selected with the largest PC,MD are denoted by a triangle. The magnitudes of error and PC,MD are represented with the gray scales, respectively. The ranges of gradient and width of the threat points in the common sets are from 486 to 499 mm/km and from 26 to 52 km, respectively, as shown in Figure 6A. As the magnitude of gradient increases, the size of error also increases. In Figure 6B, the β of 75 to 83 degrees indicates that the airborne DSIGMA monitor could hardly detect the anomalous ionospheric gradient due to the low time-variation of CCD since the ionospheric front moves in a direction nearly perpendicular to the approach path of the aircraft. Furthermore, the stationary front with a Δv of 0 to 11 m/s is hardly observable by the ground CCD monitor. Thus, β and Δv in these ranges lead to an increase in PC,MD values. As γ approaches 0°, the direction of the baseline between the ground facility and the LTP becomes perpendicular to the front. Thus the differential range error is maximized when γ gets closer to zero degrees, as shown in Figure 6C. Figure 6D indicates that fast speed profiles caused the aircraft to experience large temporal CCD rates, and the resulting carrier-smoothed code differential errors tend to be larger than those from slower speed profiles. These worst-case parameter ranges that maximize the error magnitude and PC,MDs were used as sampling limits within the Step 2 simulation.

FIGURE 6
FIGURE 6

Step 1 simulation result: Ionospheric model and geometric parameters of the scenarios selected with the largest PC,MD and the largest error. A, Front spatial structure; B, Front motion relative to AC; GF-LTP direction relative to front; D, ICAO AC speed profiles

View this table:
TABLE 2

Step 1 simulation result: Region 2 WC parameter ranges from the common sets of ionospheric gradient impact scenario

Figure 7 shows the results of the Step 2 simulation using the narrow ranges of threat parameter values determined in Step 1. Among the threat points with error greater than 2.75 m, two worst-case sets of parameter values (vectors) were determined based on the worst error (called WC A) and the worst PC,MD (called WC B), respectively. The resulting worst-case parameter vectors are presented in Table 3. The difference between the two cases is that β of WC B is slightly larger than that of WC A (ie, aircraft approach direction is a little more perpendicular to the front moving direction). Thus, compared with WC A, the approaching aircraft experiences a smaller CCD rate over time, and the resulting PMD of the airborne DSIGMA becomes larger. Both single-valued worst-case parameter vectors were then fixed and used to run two Step 3 simulations, where β (analogous to fdir) and ΔW were averaged.

FIGURE 7
FIGURE 7

Step 2 simulation result of Region 2: Selecting two individual sets of worst-case parameters

View this table:
TABLE 3

Step 2 simulation result: Region 2 WC parameter vectors

Figure 8 shows the histogram of PHMI results binned by differential range error at the LTP for the Region 2 WC B implementation of Step 3. The sample size to generate this histogram was approximately 109. The y-axis represents the PHMI on a logarithmic scale. Each binned range error value represents the critical error value for the computation of PHMI for that bin. The vertical magenta lines indicate an error of 2.75 m that is deemed to be hazardous for the GAST D GBAS service. The resulting PHMI value at an error of 2.75 m using Equation (1) was 8.46 × 10−11, which is well below the ICAO SARPs requirement of 10−9. As the binned range errors approach the required bound of 2.75 m shown by the vertical magenta lines, the PHMI for the events with larger errors shrinks rapidly because (a) the set of events meeting this criterion becomes very small relative to Ntot and (b) the combined missed detection probability PC,MD for the severe events meeting this criterion is also very small. The PHMI value for Region 2 WC A was computed as 6.73 × 10−11 and is also well below 10−9.

FIGURE 8
FIGURE 8

Step 3 simulation result: PHMI histogram for Region 2 WC B

In the same manner as Region 2, we evaluated the integrity risk for Regions 1 and 3. The resulting PHMI values computed using the worst-case vectors of WC A (selected based on the worst error) and WC B (selected based on the worst PC,MD) for Regions 1 and 3 are summarized in Table 4 with the results for Region 2. The PHMI values of both WC A and WC B for Regions 1 and 3 were also well below the ICAO SARPs requirement of 10−9. This analysis demonstrates that all combinations of the worst-case parameters specified in Regions 1, 2, and 3 were not hazardous to GAST D GBAS users using the revised probabilistic definition of the requirement proposed in this paper (combining worst-case and averaged parameters).

View this table:
TABLE 4

Step 3 simulation result: PHMI values for Regions 1, 2, and 3

6 CONCLUSIONS

This study develops a three-step method for worst-case parameter search and integrity evaluation under the impact of ionospheric anomalies on GAST D GBAS. It first defines the rationale for worst-case (WC) and average parameters and then classified the GAST D ionospheric anomaly simulation parameters into these two groups based on a specific-risk-based integrity rationale that determines when sufficient knowledge of randomness exists to allow averaging. Then, the three-step method identifies the worst-case sets of ionospheric threat parameters and evaluates the worst possible PHMI. More specifically, Step 1 performs a randomized search over the allowed ranges of all parameters with a total number of samples of about 109 and identifies one or more regions or “families” where PHMI might be high enough to exceed the overall integrity risk requirement. Step 2 determines two worst-case sets of parameter vectors for each defined worst-case region based on the largest differential range error at the LTP and the largest PMD value over all simulated approaches, respectively. For each region, both single-valued worst-case parameter vectors are then fixed and used to run two Step 3 simulations, where average parameters were allowed to vary. For the mid-latitude example considered in this paper, the resulting PHMI values for all regions were well below 10−9, showing that the GAST D SARPs integrity requirements can be met based on this combined worst-case/average definition of the anomalous ionospheric threat.

This work is a basis for research on more-detailed system design and ionospheric mitigation techniques for GAST D GBAS. More detailed studies will be performed in future work, including integrity risk evaluations using multiple monitor parameter sets and various ionospheric anomaly threat models, ionospheric monitor designs, examination of GAST D differential range errors beyond the LTP, and feasibility analysis of GAST D GBAS operation in low-latitude ionospheric environments.

HOW TO CITE THIS ARTICLE

Yoon M, Lee J, Pullen S. Integrity risk evaluation of impact of ionospheric anomalies on GAST D GBAS. NAVIGATION-US. 2020,67:223–234. https://doi.org/10.1002/navi.339

ACKNOWLEDGMENTS

The authors gratefully acknowledge Shelly Beauchamp of the Federal Aviation Administration (FAA) William J. Hughes Technical Center and her team for their support. We would like to thank the FAA Key Technical Advisors for their help and advice during the development of this methodology, including Boris Pervan of IIT, Rick Cassell of Engility Corp., and Barbara Clark, John Warburton, Randy Key, and Jed Dennis of the FAA. We would also like to thank Oliver Jeannot and Dieter Guenter of Tetra Tech AMT and Bruce Johnson, Mats Brenner, Doug Weed, and Lucas Cypriano of Honeywell for their support of this work. Moonseok Yoon was supported by the Space Core Technology Development Program of the National Research Foundation (NRF) funded by the Ministry of Science, ICT & Future Planning (NRF-2014M1A3A3A02034937) and the Korea Polar Research Institute (KOPRI, PE17900).

Footnotes

  • Funding information National Research Foundation, Grant/Award Number: NRF-2014M1A3A3A02034937; Korea Polar Research Institute, Grant/Award Number: PE17900

  • Received February 24, 2019.
  • Revision received July 22, 2019.
  • Accepted August 30, 2019.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

  1. 1.↵
    1. Murphy T,
    2. Harris M
    . More ionosphere anomaly mitigation considerations for Category II/III GBAS. In: Proceedings of the 20th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2007). Fort Worth, TX; September 2007:438-452.
  2. 2.↵
    1. Khanafseh S,
    2. Pullen S,
    3. Warburton J
    . Carrier phase ionospheric gradient ground monitor for GBAS with experimental validation. NAVIGATION. 2012;59:51-60.
  3. 3.
    1. Reuter R,
    2. Weed D,
    3. Brenner M
    . Ionosphere gradient detection for CAT III GBAS. In: Proceedings of the 25th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2012). Nashville, TN; September 2012: 2175-2183.
  4. 4.↵
    1. Belabbas B,
    2. Meurer M
    . Carrier phase and code based absolute slant ionosphere gradient monitor for GBAS. In: Proceedings of the 25th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2012). Nashville, TN; September 2012:2201-2208.
  5. 5.
    1. Jing J,
    2. Khanafseh S,
    3. Langel S,
    4. Pervan B
    . Detection and isolation of ionospheric fronts for GBAS. Proceedings of the 27th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2014). Tampa, FL; September 2014:3526–3531.
  6. 6.↵
    1. Patel J,
    2. Khanafseh S,
    3. Pervan B
    . Spatial gradient monitor for GBAS. Proceedings of IEEE/ION PLANS 2018. Monterey, CA; April 2018:571–582.
  7. 7.↵
    1. Azoulai L,
    2. Souchon F,
    3. Brunet C,
    4. Lipp A,
    5. Duchet D
    . GBAS GAST D parametric study and availability simulations to support Category III operations. In: Proceedings of the 25th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2012). Nashville, TN; September 2012: 2123-2140.
  8. 8.↵
    1. RTCA SC-159
    . Minimum operational performance standards for global positioning system/wide area augmentation system airborne equipment. RTCA Document Number DO-229E. December 2016.
  9. 9.↵
    1. Dvorska J,
    2. Podivin L,
    3. Zaviralova L,
    4. Musil M,
    5. Kren M
    . GBAS CAT II/III operational validation input for business aircraft – flight tests and simulations. In: Proceedings of the ION 2015 Pacific PNT Meeting. Honolulu, HI; April 2015:779-798.
  10. 10.↵
    1. ICAO State Letter
    . SL 17–18-AX10-Vol I-GBAS-SBAS-Approach Landing. March 30, 2017.
  11. 11.↵
    1. Kim D,
    2. Yoon M,
    3. Lee J,
    4. Pullen S,
    5. Weed D
    . Monte Carlo simulation for impact of anomalous ionospheric gradient on GAST-D GBAS. In: Proceedings of the ION 2017 Pacific PNT Meeting. Honolulu, HI; May 2017:47-55.
  12. 12.↵
    1. Pullen S,
    2. Cassell R,
    3. Johnson B,
    4. Brenner M,
    5. Weed D, et al
    . Impact of ionospheric anomalies on GBAS GAST D service and validation of relevant ICAO SARPs requirements. In: Proceedings of the 30th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2017). Portlaned, OR; September 2017:2085-2105.
  13. 13.↵
    1. Working Paper 3
    . Proposed Amendments to Annex 10: Ground-based Augmentation System (GBAS) Provisions. NSP/CSG Meeting. Presented by CAT II/III Subgroup (CSG) Rapporteur. April 2015.
  14. 14.↵
    1. Lee J,
    2. Seo J,
    3. Park YS,
    4. Pullen S,
    5. Enge P
    . Ionospheric threat mitigation by geometry screening in ground-based augmentation systems. AIAA J Aircraft. 2011;48(4):1422-1433.
  15. 15.↵
    1. Robert E,
    2. Jonas P,
    3. Vuillaume J,
    4. Salos D,
    5. Hecker L,
    6. Yaya P
    . Development of a European ionosphere threat model in support of GBAS deployment. In: Proceedings of IEEE/ION PLANS 2018. Monterey, CA; April 2018:1181-1190.
  16. 16.↵
    1. Mayer C,
    2. Belabbas B,
    3. Jakowski N,
    4. Meurer M,
    5. Dunkel W
    . Ionosphere threat space model assessment for GBAS. In: Proceedings of the 22nd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2009). Savannah, GA; September 2009:1091-1099.
  17. 17.↵
    1. Yoshihara T,
    2. Saito S,
    3. Fujii N
    . A study of nominal ionospheric gradient for GBAS (Ground-based Augmentation System) in Japan. In: Proceedings of the 23rd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2010). Portland, OR; September 2010:2689-2694.
  18. 18.↵
    1. Terkildsen M
    . GBAS ionospheric threat model evaluation: Continental Australian region. IPS Radio and Space Services. April 2010; Sydney, NSW, Australia:IPS-CR-09-01-P.
  19. 19.↵
    1. Datta-Barua S,
    2. Lee J,
    3. Pullen S, et al
    . Ionospheric threat parameterization for local area GPS-based aircraft landing systems. AIAA J Aircraft. July-August 2010;47(4):1141-1151.
  20. 20.↵
    1. Simili DV,
    2. Pervan B
    . Code-carrier divergence monitoring for the GPS Local Area Augmentation System. In: Proceedings of IEEE/ION PLANS 2006. San Diego, CA; April 2006:483-493.
  21. 21.↵
    1. RTCA SC-159
    . Minimum operational performance standards for GPS local area augmentation system airborne equipment. RTCA Document Number DO-253D. July 2017.
Back to top

In this issue

Print
Download PDF
Article Alerts
Email Article
Citation Tools
Integrity risk evaluation of impact of ionospheric anomalies on GAST D GBAS
Moonseok Yoon, Jiyun Lee, Sam Pullen
NAVIGATION: Journal of the Institute of Navigation Jun 2020, 67 (2) 223-234; DOI: 10.1002/navi.339
Twitter logo Facebook logo Mendeley logo
Bookmark this article

Jump to section