## Abstract

Ionospheric anomalies may cause large differential range errors in Ground-Based Augmentation System (GBAS) users. To mitigate those integrity threats, worst-case ionosphere-induced position errors for potentially usable satellite geometries must be bounded by the GBAS ground facility. This mitigation method requires us to compute the worst-case range error for each satellite affected by a hypothetical ionospheric front. This paper presents a simulation-based method for deriving a closed-form expression of undetected ionosphere-induced range errors. Two types of ionospheric impact scenarios are defined in terms of the motion of an ionospheric front. Explicit expressions for outputs of the code-carrier smoothing filter and the code-carrier divergence monitor are derived to reduce the computational load of ionospheric impact simulations. An exhaustive search algorithm is applied to generate the worst undetected range error among all possible ionospheric impact conditions. Finally, a closed-form expression that bounds the maximum ionospheric range errors is determined as a linear function of the magnitude of gradient and the relative speed of the ionospheric front.

## 1 INTRODUCTION

The Ground-Based Augmentation System (GBAS) provides safe landing guidance to aircraft in local areas by enhancing measurements made by the Global Positioning System (GPS) (Enge, 1999). The GBAS ground facility generates and broadcasts differential corrections and integrity information to GBAS users within several tens of kilometers of a GBAS-equipped airport.

As most error sources, including tropospheric delay and ionospheric delay, are spatially correlated between the ground facility and airborne users, those common errors are eliminated via differential corrections. GBAS-equipped airborne users can compute error bounds on the true (but unknown) position error with an extremely high confidence level in real time by utilizing the integrity parameters broadcast from the ground facility.

However, very large ionospheric spatial variations within a local area, if undetected by GBAS integrity monitors, can introduce hazardous and misleading information to GBAS users. Previous research has discovered that ionospheric storms can possibly cause extremely large ionospheric spatial gradients.

During the latter part of Solar Cycle 23, observations from the network of Continuously Operating Reference Stations (CORS) showed that an ionospheric storm on November 20, 2003, created spatial gradients as large as 412 mm/km in slant ionospheric delay (Datta-Barua et al., 2010; Lee et al., 2011). Without mitigation, vertical position errors resulting from such a large ionospheric gradient (well over 100 times the typical value) could be a few tens of meters for GBAS users when combined with the worst-case airborne/satellite geometries and approach timing (Lee et al., 2011).

A strategy to mitigate such ionosphere anomalies to support Category I (CAT I) operations includes code-carrier divergence (CCD) monitoring in the GBAS ground facility (ICAO, 2018; Simili & Pervan, 2006). The time variation of ionospheric delay induces a divergence between code and carrier measurements that creates a bias in the carrier-smoothed code output, which in turn leads to errors in differential corrections which are small under nominal ionospheric conditions. The divergence rate is continuously estimated by the CCD monitor so that very large ionospheric spatial gradients are detected at the ground facility based on the temporal gradients that they typically induce.

However, a previous study reported that near-stationary fronts of ionospheric spatial change can produce large residual ionospheric delays without detection (Luo et al., 2005). This worst-case ionospheric front impact scenario occurs when the divergence rates caused by the near-stationary fronts are too small to trigger the CCD monitor before an aircraft reaches the CAT I decision height (DH) at the end of its approach.

In other words, there exists a potential integrity risk associated with ionosphere-induced undetected range errors due to the limitation of the CCD monitor. The GBAS ground facility also contains a smoothing filter innovation monitor, a carrier-phase-based acceleration monitor, and a B-value monitor (which compares candidate pseudo-range corrections among reference receivers) which can help mitigate ionospheric threats (Pullen, 2017). However, they add little to the effectiveness of the CCD monitor against worst-case threats.

One possible way to mitigate this integrity threat for CAT I GBAS users is to assume that the worst possible ionospheric anomaly always exists. The worst-case ionosphere-induced position errors for the satellite geometry that an aircraft might use for positioning can be calculated at the GBAS ground facility in near real time, as described in (Lee et al., 2011; Pullen et al., 2009). If there are satellite geometries that have acceptable vertical protection levels (VPLs) but have unacceptable potential position errors owing to potential ionospheric anomalies, these geometries are screened out by inflating (increasing) certain parameters broadcast from ground to aircraft in order to increase aircraft VPLs beyond the vertical alert limit (VAL) to which VPLs are compared for safety assurance. This mitigation technique is known as *geometry screening* and is illustrated by the functional flowchart in Figure 1 (Pullen et al., 2009).

The box labeled *Worst-Case Range Error Determination* in Figure 1 requires a procedure to estimate the worst-case ionosphere-induced range error in the presence of a hypothetical ionospheric front. Real-time estimation of this error is currently infeasible due to both its computational demands and the undesirability of including exhaustive search routines (with a very large range of possible outputs) in safety-critical software. An extensive search is required because the differential range error is affected by numerous ionospheric and aircraft approach geometry and timing factors (e.g., at what point during an aircraft approach procedure does the leading edge of an ionospheric spatial gradient begin to affect it).

The method originally used in GBAS geometry screening analyses and implementations was to generate arrays of differential range errors ahead of time via extensive simulation, using these as look-up tables when doing real-time GBAS mitigation calculations (Lee et al., 2006). This approach remains feasible but cumbersome, as the range-domain simulations that generate these tables have to be repeated every time that the GBAS algorithms or operational parameters change.

Another solution to this difficulty is to generate bounding (conservative) closed-form expressions in advance that are flexible enough to cover all possibilities. This paper presents a derivation of a closed-form expression for undetected differential range errors due to ionospheric anomalies. A simplified set of expressions was developed by Enge (2007) but was not published or made widely available (Yoon et al., 2016).

This model was described in subsequent publications and was used to estimate ionosphere-induced range errors for GBAS mitigation simulations (Lee et al., 2011). However, the simplifications used in the original model abstracted out important features of the behavior of these errors. In particular, non-linear varying errors as a function of relative front speed were approximated to have a constant value. The resulting inaccuracies would become more significant as the maximum ionospheric gradient increased, as it does under equatorial ionospheric conditions (for users within ±20 degrees of the magnetic equator).

This motivated the derivation of improved closed-form expressions based on an empirical search for the worst-case error. The improved model was first used to express and bound ionospheric errors for a proposed GBAS installation in Brazil (Yoon et al., 2016, 2019). It is only slightly more complex than the original model and thus just as easy to implement in real-time GBAS ionospheric mitigation algorithms.

This paper is the first to fully explain the background and derivation of both the original and enhanced models. The remainder of this paper is organized as follows. Section 2 presents the ionospheric front model and impact scenarios of an ionospheric spatial gradient affecting a GBAS user.

Section 3 provides the background behind the calculation of undetected ionosphere-induced range error in GBAS. Explicit expressions for ground and airborne code-carrier smoothing filters and the ground CCD monitor are derived to reduce the computational burden of performing a vast number of simulations of all possible combinations of ionospheric threat parameters.

In Section 4, the revised method for deriving the closed-form expression is described in detail. An exhaustive search for the largest ionospheric differential range error from specific sets of ionospheric parameters is conducted, and a conservative linear model that bounds the worst-case errors is determined. The conclusion of this study is presented in Section 5.

## 2 IONOSPHERIC FRONT MODEL AND IMPACT SCENARIOS

### 2.1 Moving ionospheric front model

A semi-infinite wedge-shaped ionospheric front model was proposed in a previous study (Datta-Barua et al., 2010). Figure 2 illustrates this ionospheric front model described with four parameters: front gradient, *g* (mm/km); front width, *w* (km); front propagation speed relative to the Ionospheric Pierce Point (IPP) of the ground facility, Δ*υ* (m/s); and initial distance of the front from the ground facility, *D* (km).

The gradient (slope) represents a linear change in the slant ionospheric delay between high and low delay regions. Note that the product of gradient and width (*g* · *w*) yields a total delay change across the gradient. The model assumes that these parameters remain constant for at least the duration of a single aircraft approach.

Because both IPPs of the aircraft and ground facility (*IPP*_{AC} and *IPP*_{GF}) are spatially close to each other within a GBAS service area compared to the distance to a GPS satellite, the velocity components of both IPPs caused by satellite motion are common. Thus, the relative speed between two IPPs is approximated as the aircraft approach speed.

In this study, we assume the aircraft moves toward a runway (in line with that runway) with a constant speed of 70 m/s. The initial position of the aircraft IPP used in ionospheric impact scenarios varies depending on the relative motion of the aircraft to the ionospheric front, which is explained in the following section. The variable *x* denotes the distance of an evaluation point which is at the minimum 200-foot approach DH for CAT-I operations and is assumed to be located 6 km from the geographic centroid of the multiple reference receivers and antennas comprising the GBAS ground facility.

### 2.2 Ionospheric impact scenario

In terms of the motion of an ionospheric front relative to an aircraft, ionospheric impact scenarios are classified into two categories: *slow-moving* and *fast-moving*. The propagation speed of a *slow-moving* front (relative to that of the ground IPP for a particular satellite) is slower than that of the aircraft IPP for the same satellite relative to the ground IPP (Δ*υ* < 70 m/s). Alternatively, if the front relative to the ground IPP is moving faster than the aircraft IPP relative to the ground IPP (Δ*υ* > 70 m/s), it’s a *fast-moving* scenario.

In the slow-moving scenario, the initial position of the aircraft IPP is set at the rear (trailing) edge of the wave in order to ensure that the aircraft IPP experiences the entire ionospheric delay change immediately after the beginning of the simulation. For the same reason, the initial position of the aircraft IPP is set at the front (leading) edge of the front in the fast-moving front scenario. The initial position of the aircraft IPP is determined in order to deliberately simulate the situation that produces the largest differential range errors.

Figure 3(a) shows an example of the temporal change in ionospheric delay observed by the aircraft (solid line) and the ground facility (dashed line) for the slow-moving front scenario. Given that the aircraft IPP approaches from directly behind the ionospheric front, the ionospheric delay at the aircraft linearly decreases until the aircraft IPP entirely passes through the front. *t*_{DH} denotes the epoch when the aircraft reaches the decision height (DH), representing the end of CAT-I approach guidance.

Note that the aircraft IPP could be located in either the middle of the gradient or in the low-delay region when it arrives at the DH. On the other hand, the ionospheric delay observed from the ground facility remains low until the arrival of the ionospheric front at *t*_{GF} starts to increase. In Equation (1), *t*_{DH} and *t*_{GF} for the slow-moving front scenario is derived using the parameters described in Figure 2.

1

It is worth noting that *t*_{GF} can be larger (i.e., occuring later) than *t*_{DH} when Δ*v* is sufficiently small, indicating that the front cannot affect the ground facility at all until the aircraft reaches the DH. Furthermore, if Δ*υ* is near zero, the front would appear to be stationary with respect to the ground facility IPP. The worst differential range errors usually occur under these conditions because the ground monitors are not able to observe and detect the ionospheric front before the aircraft arrives at the DH (Luo et al., 2005).

The temporal ionospheric delay changes for the fast-moving front scenario as illustrated in Figure 3(b). The ionospheric delay observed from the aircraft linearly increases until the ionospheric front entirely overtakes the aircraft IPP. For the fast-moving front scenario, *t*_{DH} is derived according to Equation (2).

This is different from that of the slow-moving front scenario, while *t*_{GF} remains the same as Equation (1). Note that what differs between the two scenarios is the initial position of the aircraft, which only affects its arrival time at the DH:

2

## 3 DERIVATION OF UNDETECTED IONOSPHERE-INDUCED DIFFERENTIAL RANGE ERROR

The two key components of GBAS that are relevant to this threat, namely code-carrier smoothing and CCD monitoring, must be taken into account to derive the undetected ionosphere-induced differential range error. The GBAS Approach Service Type (GAST) C considered in this paper is intended to support precision approach operations down to CAT-I weather minima.

A GBAS facility which provides GAST-C service uses the code-carrier smoothing filter to mitigate the effects of receiver noise and multipath errors on GPS code (pseudo-range) measurements at both the aircraft and the ground facility. However, code-carrier divergence over time due to ionospheric delay adds up in the smoothing filter and introduces additional error to carrier-smoothed code measurements. The code-carrier divergence rate is continuously monitored by the CCD monitor implemented in the ground facility.

Figure 4 illustrates a block diagram of the calculation of undetected range errors. The ionosphere-induced differential range error is computed by subtracting the code-carrier smoothing filter output of the ionospheric delay at the ground facility from the similar output at the aircraft. The time series of differential range errors and the output of the CCD monitor are obtained as the aircraft advances on its path toward the runway.

If the output of the CCD monitor exceeds its minimum detectable divergence rate (MDDR) before the aircraft reaches the DH, the differential range error for the corresponding scenario is discarded from the exhaustive search for the worst undetected range error (more details are in Section 4.1). Otherwise, the differential range error at the DH is recorded as the undetected range error for that scenario.

The recursive approach to computing the outputs of the two code-carrier smoothing filters and the ground CCD monitor (itself composed of two low-pass filters in series) at each epoch is computationally intensive. Moreover, a vast number of simulations need to be conducted to compute outputs from all possible combinations of ionospheric threat parameters and aircraft-ground facility ionospheric front geometries. Thus, we introduce explicit expressions for the outputs of code-carrier smoothing and the CCD monitor to compute undetected differential range errors efficiently in this section. The results of applying these equations will be used to search for the worst differential range error in Section 4.

### 3.1 Range error expression with carrier-smoothed code

The carrier-smoothed code measurement, denoted as , at the current epoch is expressed as follows:

3

where *ρ*(*t*) and *φ*(*t*) denote the single frequency code and carrier phase measurements, respectively, at a given epoch *t*. The carrier-smoothed code measurement at the previous epoch, *t* − *T*, is projected to the current epoch *t* by adding the change in the carrier phase measurement between the current and previous epochs (*T* represents the interval between samples). A recursive filter of length *M* weights this carrier-based measurement more heavily than the current code measurement.

The explicit expression for ionospheric delays is derived from Equation (3) using the simplified ionospheric front model of Section 2.1 whose ramp has a constant slope or gradient. The code-carrier smoothing filter outputs for the ramp inputs of increasing and decreasing delays are denoted as *C*_{inc}(*t*, *g*_{t}, *T*_{w}) and *C*_{dec}(*t*, *g*_{t},*T*_{w}), respectively, in Equation (4):

4

where *g*_{t} is the temporal gradient, and *τ* is the time constant of the code-carrier smoothing filter, which is equal to the sampling interval *T* times the filter length *M* (both ground and airborne receivers independently employ smoothing filters on their measurements with the same time constant).

In the ionospheric front impact model, the ionospheric delay increases or decreases between epochs zero and the temporal width *T*_{w}, the time required for the ionospheric front to completely pass through the relevant aircraft or ground IPP. The derivation of these functions is given in the Section A of the appendix.

As shown in Figure 3 for a ramp input of decreasing delay, the temporal change in ionospheric delay observed at the aircraft varies depending on the relative motion of the front with respect to the aircraft. Specifically, the aircraft experiences decreasing delays for the slow-moving front scenario or increasing delays for the fast-moving front scenario.

Figure 5 shows the smoothing filter output of ionospheric delays of the aircraft for slow-moving (Δ*υ* = 50 m/s) and fast-moving (Δ*υ* = 100 m/s) front cases when *g* = 500 mm/km and *w* = 25 km.

In contrast, the ground facility experiences ionospheric delays that increase for both scenarios as the ionospheric front passes. The differential range errors *ɛ* for both scenarios are then defined as a combination of smoothed range errors according to Equation (5):

5

where the temporal ionosphere parameters are:

6

Finally, the differential range error when the aircraft reaches the evaluation point, DH, can be computed by substituting *t*_{DH} [from Equation (1) and Equation (2)] for *t* in Equation (5).

Note that only *low-to-high* fronts as shown in Figures 2 and 3 (for which the ionospheric delay of the ground facility goes from low to high) are taken into account in this simulation search procedure.

For *high-to-low* fronts, which are the reverse cases, the delay of the ground facility decreases for both slow-moving and fast-moving scenarios. Since the high-to-low front cases give symmetric results as those in Equation (5), no searching among these scenarios is required.

### 3.2 Code-carrier divergence monitor

The derivation for the test statistic of the CCD monitor is described as follows. The amount of code-carrier divergence Δ is simply expressed as:

7

Because the ionospheric delay has an opposite sign in code and carrier measurements, the amount of code-carrier divergence is twice the ionospheric delay *I*(*t*):

8

While this analysis does not directly represent receiver noise and multipath errors on single frequency code and carrier measurements, they are accounted for as nominal errors in the derivation of detection statistics to follow.

With the sampling interval denoted by *T* as before, the divergence rate *d*Δ is defined as follows:

9

Common errors such as tropospheric delays, satellite and receiver clock offsets, and ephemeris errors are removed by subtracting divergences between two consecutive epochs. Assuming that there is no cycle slip within the interval *T*, the carrier-phase integer ambiguity is eliminated in Equation (9).

A typical CCD monitor for GBAS is composed of a series of two cascaded first-order low-pass filters with the divergence rate from Equation (9) as the input (Simili & Pervan, 2006):

10

where *Z*_{i}(*t*) represents the output of each filter at a corresponding epoch with a filter weight of . The filter time constant of the CCD monitor *τccd* is different from the time constant of the ground and airborne code-carrier smoothing filters.

An explicit expression for the test statistic of the CCD monitor is given by Equation (11) using the unit step function *u*. The derivation of this expression is provided in Section B of the appendix.

Since the ionospheric front starts to affect the ground facility after *t*_{GF}, the test statistic can be computed by substituting *t* − *t*_{GF} for *t* in Equation (11):

11

In order to account for the effect of measurement noise in Equation (11), the minimum detectable divergence rate *MDDR*_{ccd} is used to guarantee a sufficiently low probability of missed detection according to Equation (12):

12

where *σ*_{d} is the fault-free standard deviation of the output of the CCD monitor *Z*_{2}, which is the test statistic. The value of *σ*_{d} was estimated as 0.00399 m/s in (Simili & Pervan, 2006). The probability of fault-free detection allocation *P*_{ffd} is 1 × 10^{− 7} per 15 seconds, and the probability of missed detection is 1 × 10^{−4} per approach for all satellites.

By assuming a maximum of 10 critical satellites, this study uses a P_{ffd} of 1 × 10^{−8} and a P_{md} of 1 × 10^{−5} per satellite. A *k*_{ffd} of 5.73 and a *k*_{md} of 4.26 are constant multipliers selected to meet the allocated continuity and integrity requirements respectively, and *k*_{ffd} + *k*_{md} is set to 10. Any CCD output above *MDDRccd* is assumed to be detected with the probability represented by *k*_{md}, assuming the monitor threshold has been set as *k*_{ffd}*σ*_{d}.

Figure 6 shows the filter output of the CCD monitor as a function of time when the ionospheric threat parameters are assumed to be Δ*υ* = 100 m/s, *g* = 500 mm/km, and *w* = 25 km. The CCD monitor time constant is set to be 30 s (Simili & Pervan, 2006). Detection with a probability of at least one minus the required missed-detection probability (P_{md}; from which *k*_{md} is derived) is assumed to occur when the output of the CCD monitor crosses *MDDRccd*.

Figure 7 shows the CCD monitor outputs with different Δ*υ*. For slow-moving front cases, the CCD output tends to be smaller, while the period in which the output is affected by the front tends to be longer.

## 4 GENERATION OF CLOSED-FORM EXPRESSION

### 4.1 Exhaustive search for the largest undetected ionosphere-induced differential range error

A particular ionospheric impact scenario includes specific values for each of the four threat parameters, which leads to a large number of permutations. An exhaustive search algorithm is conducted offline to search for the largest undetected range error among all possible ionospheric anomaly impact scenarios. The goal is to derive a bounding closed-form expression for the worst-case undetected ionosphere-induced differential range error as a function of threat parameters Δ*υ* and *g*.

The procedure of the exhaustive search algorithm is described as follows. Recall that we derived explicit expressions for the smoothing filter and CCD monitor outputs in Equation (4) and Equation (11).

For a particular set of ionospheric front parameters, we calculate both the differential range error and the output of the CCD monitor only at the arrival time of the aircraft at the DH, rather than during each epoch of the approach, which greatly reduces the computational load. For given values of *w*, Δ*υ*, and *g*, a series of differential range errors and CCD monitor outputs at *t*_{DH} can be obtained across the possible values of *D*, as illustrated in Figure 8.

*D*_{CCD} is the initial distance of the front from the ground facility for which the CCD monitor output at *t*_{DH} is *MDDRccd*. For the slow-moving front scenario, cases with *D* ≤ *D*_{CCD} are regarded as detected cases. Since the front moves slower than the aircraft IPP, the period of the ground IPP being affected by the front gets longer for a shorter *D*, which causes a larger CCD monitor output. Thus, detection occurs before the aircraft arrives at the DH in these cases.

Conversely, cases with *D* ≥ *DCCD* are regarded as detected cases for the fast-moving front scenario. The worst-case undetected range error is then chosen as the maximum of the differential range errors among the undetected cases for all possible values of *D* and *w*. This process is repeated for each set of values for Δ*υ* and *g*.

The ranges for the parameters and step sizes used in the exhaustive search are specified in Table 1. The bounds for the search space associated with the ionospheric front are determined based on the ionospheric threat model for the Conterminous United States (CONUS) (Lee et al., 2017).

The upper bound on *D* (see Figure 2) is purposely set to be an unrealistically large value for CAT-I precision approach procedures in order to produce the worst possible undetected range error for a particular satellite when Δ*υ* approaches the aircraft speed. A gradient greater than 50/*w* is discarded by the search process because the largest change in slant ionospheric delay among the observed ionospheric anomalies in CONUS was no greater than 50 m (Datta-Barua et al., 2010).

The step sizes of the parameters, which directly relate to computation time, were determined to be sufficiently dense such that the ionospheric impact scenario producing the largest undetected error would not be missed. For the ranges specified in Table 1, the number of possible ionospheric impact scenarios is approximately 1.6 × 10^{11}.

The offline computation time of the exhaustive search over all these cases was measured using MATLAB’s built-in functions. Using a Windows PC with an Intel Core i7-10700 CPU, the entire search procedure with MATLAB scripts took about 2.5 hours. However, as Figure 8 illustrates, an upper bound on *D* of 100 km (instead of 100,000 km) is sufficient to generate the largest undetected error in most cases. The unrealistically large value of *D* = 100,000 km is only needed for cases where Δ*υ* is nearly equal to *υ*_{ac}.

The computation time is significantly reduced to less than 1 minute by applying an upper bound on *D* of 100 km. If we specify targeted bounds on *D* with respect to Δ*υ*, the computation time can be greatly reduced without omitting the worst possible undetected errors for the entire range of values for Δ*υ*.

Figure 9 shows the worst (largest-magnitude) ionosphere-induced differential range errors *ɛ* over all combinations of *w* and *D* as a function of Δ*υ* for three different values of *g*. One can observe that the curves of *ɛ* have similar trends, which can be divided into three regions in terms of Δ*υ*.

First, in the region of small Δ*υ*, where the error has a constant value, the worst differential range error is computed as its maximum possible value. In this small Δ*υ* case, the CCD monitor does not capture the ionospheric front in at least one permutation among all possible *w* and *D* combinations, because the change of ionospheric delay over time is too small to be observed.

In such a case, the maximum undetected range error could occur as a combination of the physical separation *x* and the synthetic separation due to the memory of the gradient remaining in the code-carrier smoothing filter. Thus, the maximum range error can be expressed as a product of the magnitude of the gradient *g* and the effective separation between the aircraft and ground facility, which is a sum of the physical and synthetic separation.

Second, in the region of moderate Δ*υ*, the worst differential range error starts to decrease considerably. This is because the ionospheric delay changes in time quickly enough to be observed by the ground facility. In every permutation of this region, the CCD monitor detects the front before the differential range error reaches its maximum value. As Δ*υ* becomes larger, the range error continuously decreases until the minimum value of the undetected error is reached.

Third, in the region of fast Δ*υ*, a transition of the error state takes place. The worst differential range error occurs among the cases in which the ionospheric front reaches the ground and is detected right after the aircraft arrives at the evaluation point. In these cases, the differential range errors are only caused by the code-carrier smoothed ionospheric delays of the aircraft, not by those of the ground facility.

Since carrier measurements dominate in the code-carrier smoothing filter directly after the impact of the front, the initial response of the carrier-smoothed code has an overshoot in the negative direction (Luo et al., 2003). The depth of this overshoot is proportional to the temporal gradient according to Equation (4).

As the effect of negative errors in carrier-smoothed code measurements of the aircraft increases in the region of fast Δ*υ*, the worst differential range error increases. When Δ*υ* is sufficiently large, the ionospheric front passes the aircraft so quickly that the worst differential range error is totally due to the effect of physical separation. Thus, the worst differential range error converges to *g**x*, while the memory effect within the code-carrier smoothing filter decreases to zero.

### 4.2 Linear closed-form expression

This section explains the procedure for developing a linear closed-form expression for the worst-case ionospheric differential range error. Because it is difficult to find an analytical function of the worst undetected range error generated by the exhaustive search, a conservative linear function (i.e., one that upper-bounds all the simulated range errors) is created for simplicity.

Note that, since the results of future exhaustive searches may change when the simulation parameters (such as the threat parameters, filter time constants, monitor thresholds, aircraft speeds, DH distance with a variation of more than a couple of kilometers, and/or runway configurations) change, the linear function may need to be redesigned to upper-bound the results of updated searches.

Figure 10 shows an example result of the worst-case undetected differential range errors induced by the impact of ionospheric fronts with a gradient of 500 mm/km. Simulated range errors as a function of Δ*υ* were obtained by using the exhaustive search algorithm (red dashed curve). A linear closed-form expression (black solid line) was determined to fit the simulation result.

In this process, two transition points of the linear expression are defined. Point *a* indicates the minimum Δ*υ* value beyond which the CCD monitor is capable of detecting the ionospheric front with the required probability as determined by the results of the exhaustive search. For a particular gradient *g*, point *a* can be derived analytically from the mathematical relationship of the gradient of the ionospheric front and the minimum detectable divergence rate:

13

Point *b* denotes a Δ*υ* for which the worst undetected range error is the smallest value. We performed a curve fitting method to model *b* as a function of the gradient *g* from 200 mm/km to 500 mm/km.

In Figure 11, the simulated *b* values are marked with blue dots, and the fitted curve is represented by the red curve. The fitted curve of point *b* is expressed as follows:

14

Using transition points *a* and *b*, a bounding linear closed-form expression for the worst undetected range error can be expressed in Equation (15) as a function of Δ*υ* and *g*:

15

## 5 CONCLUSION

This study develops a closed-form expression for the undetected differential range error for CAT-I Ground-Based Augmentation System (GBAS) users under anomalous ionospheric conditions. Using the simplified moving front model, ionospheric impact scenarios are defined to help formulate a description of the resulting undetected range errors.

The key steps in calculating these errors include the derivation of expressions for the outputs of the code-carrier smoothing filter and the ground code-carrier divergence (CCD) monitor. We then conducted exhaustive simulations to search for the worst-case undetected differential range error among all possible combinations of ionospheric threat parameters and aircraft approach geometries. Finally, a closed-form expression is derived as a linear function that conservatively bounds the worst-case undetected range errors for all combinations of parameters.

While not all possible approach geometries were considered, a similar technique could be used to derive closed-form solutions for other GBAS service types, such as GAST-D for CAT II/III approaches, which employ aircraft monitoring as well. For example, it is known that varying aircraft speeds and runways not aligned with the vector between the GBAS reference receiver centroid and the evaluation point (DH) can result in different worst-case errors.

The proposed closed-form expression for the maximum undetected differential range error can be utilized to calculate worst-case ionosphere-induced position errors within the GBAS ground facility to screen out potentially unsafe satellite geometries in GAST-C operations. This work can also be applied to the development of GBAS for not only CONUS but also other regions that have different ionospheric threat models.

However, the exhaustive search must be repeated to recalculate a model for the closed-form expression when the simulation parameters are changed. Those include the ranges of ionospheric threat parameters, time constants of the code-carrier smoothing filter and CCD monitor, CCD monitor threshold, and distance from the centroid of GBAS reference receivers to the evaluation point. The current closed-form expression can be further improved by simulating additional aircraft profiles or considering aircraft speed to be an additional simulation parameter when modeling differential range errors.

## HOW TO CITE THIS ARTICLE

Kim, D., Yoon, M., Pullen, S., & Lee, J. (2021). Closed-form analysis of undetected range errors due to ionospheric impacts for GBAS category I operations. *NAVIGATION*, *68,* 507–519. https://doi.org/10.1002/navi.442

## ACKNOWLEDGMENTS

As explained in the Introduction, the origin of this research was the work of Prof. Per Enge in deriving the original simplified set of closed-form ionosphere-induced range error equations. The work reported here would not have been possible without his contributions. Sadly, Prof. Enge passed away in 2018 and was unable to review this research. We all miss him greatly.

We would also like to thank the FAA Satellite Navigation Program Office and the FAA William J. Hughes Technical Center for their support of this research over many years. Many people within the FAA supported this work, including Leo Eldredge, John Warburton, Randy Key, Barbara Clark, and Shelly Beauchamp. In addition, many FAA Key Technical Advisors and GBAS industry members (too many to name) provided helpful comments and recommendations.

Dongwoo Kim was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2020R1A2C1011745) and the Korea Polar Research Institute (PE19900).

## APPENDIX A: CARRIER-SMOOTHED-CODE MEASUREMENT

If we only consider the effect of ionospheric delay on pseudo-range error; then the carrier smoothed code measurements in Equation (3) are simply expressed as *ρ*(*t*) = *I*(*t*) and *φ*(*t*) = − I(*t*) yields :

The expression for the carrier-smoothed code measurement in the Laplace domain is as follows:

According to the Padé approximation, exponential terms can be expanded as a series of polynomial functions. Here, we approximate the exponential function with the first-order Padé approximation:

Using the above approximation, the transfer function of the smoothing filter can be derived in Equation (A4). Both the aircraft and the ground facility use a smoothing time constant *τ* of 100 s, which is equal to a sampling interval *T* of 0.5 s times a filter length *M* of 200 (RTCA, 2017).

We have two types of ionospheric delay changes: increasing and decreasing. For the temporal gradient of *g*_{t} and the temporal width of *T*_{w}, the increasing type of I(*t*) can be simply expressed with unit step functions:

By taking the Laplace transform of I(*t*):

Then, an explicit expression for *C*(*t*) can be obtained by taking the inverse Laplace transform of *H*(*s*)*I*(*s*):

Finally, we obtain an explicit expression of *C*(*t*) for the increasing input:

For the decreasing input case, the equation is derived using the same steps.

## APPENDIX B: OUTPUT OF THE CODE-CARRIER DIVERGENCE MONITOR

By taking the Laplace transform of Equation (10), the output of the filter can be expressed in the *s* domain:

Because the CCD monitor is present in the ground facility, the ionospheric delay will increase linearly as the front passes. Similar to Equation (A6), the increasing input of the ground facility for the temporal parameters *g*_{GF} and *T*_{w} in s domain is derived as:

Then, *Z*_{2}(*s*) yields:

The inverse Laplace transform function *Z*_{2}(*t*) can be expressed as:

In the case of decreasing I(*t*) (i.e., the *high-to-low* front scenario), the output of the CCD monitor will be symmetric to that of the increasing type of I(*t*) in (B4) and thus have an opposite sign.

## Footnotes

**Funding information**National Research Foundation of Korea, Grant/Award Number: 2020R1A2C1011745; Korea Polar Research Institute, Grant/Award Number: PE19900

- Received August 19, 2020.
- Revision received July 26, 2021.
- Accepted July 31, 2021.

- © 2021 Institute of Navigation

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.