Enhancing GAST D Availability by Using an Ionospheric Field Monitor

2.1.1 Simulation Design

To evaluate the GAST D performance, we first computed the residual ionospheric range error (E_IG) with a selected set of integrity monitor thresholds by repeating a massive number of simulations with various parameters, including ionospheric delay gradient, GBAS ground station location, aircraft (approach speed and direction), and GNSS satellite (speed and direction of ionospheric pierce point [IPP]) parameters. In each simulation, the residual ionospheric delay errors and their missed-detection probabilities are computed. The maximum error with a missed-detection probability higher than a certain tolerable level is defined as E_IG. The improvement in availability was evaluated based on E_IG values obtained with and without the IFM.

Figure 1(a) shows the simulation geometries of a standard GAST D. An IFM station was introduced for a specified direction and distance from the runway threshold (Figure 1(b)) as described by Pullen et al. (2017). In each simulation run, ionospheric delays of a GNSS satellite observed by the aircraft, GBAS reference receiver, and IFM station were calculated at each epoch to obtain the ionospheric delay gradient. Instead of computing the IPP velocity from the orbit, the IPP velocity with respect to the ground was sampled within a possible range of IPP velocities.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 1

(a) Simulation geometry with the standard GAST D configuration; (b) simulation geometry with the IFM; rxs: receivers

The aircraft speed profile was obtained from the International Civil Aviation Organization (ICAO) (2023). The ionospheric gradient parameters were adapted from values for the Asia and Pacific regions (Saito et al., 2017; Saito & Yoshihara, 2017) and from values used in the validation of GAST D standards (ICAO, 2023), as summarized in Table 1. The effect of GNSS satellite motion is represented by the motion of the IPP at which the satellite-to-receiver path crosses a certain altitude. The IPP altitude was assumed to be 350 km, which is the typical altitude for an ionospheric density peak. The satellite IPP motion depends on the GBAS ground station location. Figure 2 shows the IPP speed and direction at an altitude of 350 km at Tokyo (36˚N, 139˚E), the North Pole, and the equator, calculated by using the GPS standard 24-satellite constellation (RTCA, 2020). The IPP speed range was almost equivalent for different latitudes, at 50–600 m/s. However, the IPP direction varied based on the satellite inclination and the geographic latitude of the location. At the equator, no IPP moves westward; thus, the direction was confined to a range of 0˚–180˚ (clockwise from the north). In contrast, at the North and South Poles, the IPP direction was distributed across all angles. At Tokyo, which is a mid-latitude region, there were narrow ranges of directions that were not occupied by satellites. To be conservative, all theoretically possible IPP locations and directions were considered in this study. The satellite IPP motion parameters are listed in Table 2.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 2

GPS satellite IPP velocity and direction at an altitude of 350 km at (a) Tokyo, (b) the North Pole, and (c) the equator; CW: clockwise

Based on the calculated ionospheric delays for the aircraft and ground reference receivers, carrier smoothing using the Hatch filter (Hatch, 1983) was applied with two time constants, 30 and 100 s (ICAO, 2023), to describe the accumulation of time-varying ionospheric delays associated with the Hatch filter. The code-carrier divergence (CCD) associated with the ionospheric delay variation was also calculated for the aircraft and ground reference receiver. For a 30-s smoothed ionospheric delay, the residual ionospheric delay errors (i.e., the difference between the 30-s smoothed ionospheric delays of the aircraft and ground reference receiver) were calculated at the time when the aircraft arrived at the runway threshold. Based on the ionospheric delays (raw, 30-s smoothed, and 100-s smoothed) and CCD values, the missed-detection probability of the error (P_md) was estimated for the integrity monitors utilized by the GAST D system. If P_md was lower than 10^–9, the error was assumed to be well detected and excluded by the integrity monitors; in contrast, if the P_md was greater than 10^–9, the error was considered to pose an unacceptable risk.

2.1.2 Integrity Monitors

For a minimally compliant standard GAST D system, the ground CCD monitor, ground ionospheric gradient monitor (IGM), airborne dual smoothing pseudorange ionospheric gradient monitoring algorithm (DSIGMA), and airborne CCD monitor are used for ionosphere anomaly monitoring (ICAO, 2023; RTCA, 2019). The CCD monitor detects different behaviors in the code and carrier. Ionospheric delay gradients are one of the causes of CCD. IGM monitors the ionospheric delay gradient among the ground reference receivers. The IGM test statistic is the ionospheric delay gradient at the GBAS reference point estimated by the simulation. For simplicity, the spatial distribution of the individual ground reference receivers is not taken into account. DSIGMA monitors the difference between the 100-s and 30-s pseudoranges between the ground and air and detects ionospheric delay gradients due to different accumulated errors with different smoothing time constants. The ionospheric delay is the main cause of increased DSIGMA test statistics. DISGMA test statistics are computed with 30-s smoothed and 100-s smoothed ionospheric delays for the ground and air. The total missed-detection probability (P_{md, tot}) is taken as a function of the missed-detection probability of each monitor (P_{md, CCDgnd} for the ground CCD, P_{md, IGM} for the IGM and P_{md, DSIGMA} for DSIGMA). The airborne CCD monitor is not used in computing P_{md, tot} because it is highly correlated with DSIGMA (Pullen et al., 2017). The total missed-detection probability in the absence of the IFM is defined as follows:

$P_{md,tot}=\min \left(P_{md,IGM}\cdot P_{md,DSIGMA},P_{md,IGM}\cdot P_{md,CCDgnd}\right)$ 1

because a combination of the IGM and DSIGMA or that of the IGM and ground CCD monitor is sufficient to detect an anomaly.

For simulations using the IFM, IFM-based gradient monitoring was considered as an additional integrity monitor that works in parallel with the IGM; that is, the P_md value for both IGM and IFM is the smaller value of P_md, _IGM or P_md, _IFM. The total missed-detection probability with the IFM is described as follows:

$P_{md,tot}=\min \left(\min \left(P_{md,IGM},P_{md,IFM}\right)\cdot P_{md,DSIGMA},\min \left(P_{md,IGM},P_{md,IFM}\right)\cdot P_{md,CCDgnd}\right)$ 2

The integrity monitor parameters used in this study are summarized in Table 3. P_md for each integrity monitor with test statistics x and its standard deviation for integrity (σ_int)are computed as follows:

$P_{md}=\frac{1}{2}\left[1+erf\left(\frac{x_{th}-x}{\sqrt{2{\sigma }_{\mathrm{int}}^2}}\right)\right]$ 3

where erf is the error function. The threshold value (x_th) is given by the following:

$x_{th}=K_{FFD}\cdot {\sigma }_{cont}$ 4

The parameters for the IGM, ground CCD, and DSIGMA were obtained from the work by Pullen et al. (2017). The parameters for the IFM are assumed to be the same as those of the IGM, which are based on the algorithm utilized for the IFM (Saito et al., 2012).

The minimum detectable error (MDE) is given as follows:

$MDE=K_{md}\cdot {\sigma }_{\mathrm{int}}+K_{FFD}\cdot {\sigma }_{cont}$ 5

When the monitor test statistic exceeds the MDE, an anomaly is detected by the monitor with a missed-detection probability of less than 10^–9. Excluded satellites may be re-admitted when the monitor test statistic is reduced to a lower value. Pullen et al. (2017) suggested that satellites could be re-admitted when the monitor test statistic becomes lower than (x_th – 2*σ_int). In this study, satellites are considered for re-admission at any time for the ground CCD and DSIGMA, which are based on smoothed pseudoranges and have filter delay effects. This conservative assumption results in higher residual errors. For the IGM and IFM, which are based on instantaneous carrier-phase measurements, we assume that satellites can be re-admitted 200 s after the test statistic becomes lower than the MDE. A waiting time of 200 s was chosen to ensure that the smoothing filter for the pseudorange with a time constant of 100 s reaches a steady state after initialization.

2.2.1 GAST D Protection Level

The first condition is given as follows:

$VPL<VAL$ 6

and:

$VPL=K_{ffmd}{\sigma }_{vert}+D_V$ 7

where:

$K_{ffmd}=5.85$ 8

and:

${\sigma }_{vert}=\sqrt{\sum _{i=i}^N{S}_{Apr\_vert,i}^2\times {\sigma }_{i_x}^2}$ 9

Here, S_{Apr_vert, i} is the projection of range to vertical errors with respect to the approach path, with a glide path angle (GPA) of θ_GPA given by the following:

$S_{Apr\_vert,i}=S_{zi}+S_{x,i}\tan {\theta }_{GPA}$10

S is the projection matrix from the range to position given as follows:

$\begin{array}{cc}S& ={(G^{T}WG)}^{-1}{G}^{T}W\\ & =\left[\begin{array}{cccc}{S}_{x,1}& S_{x,2}& \mathrm{...}& S_{x,N}\\ S_{y,1}& S_{y,2}& \mathrm{...}& S_{y,N}\\ S_{z,1}& S_{z,2}& \mathrm{...}& S_{z,N}\\ S_{t,1}& S_{t,2}& \mathrm{...}& S_{t,N}\end{array}\right]\end{array}$ 11

where G is the geometry design matrix consisting of the azimuth and elevation angles (Az and El, respectively) of the i-th satellite:

$G_i=[\cos E{l}_i\cos A{z}_i,-\cos E{l}_i\sin A{z}_i-\sin E{l}_i,1]$ 12

The weight matrix W is given by the following:

$\begin{array}{cc}{W}^{-1}& \left[\begin{array}{cccc}{\sigma }_{1_x}^2& 0& \cdots & 0\\ 0& {\sigma }_{2_x}^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }_{N_x}^2\end{array}\right]\end{array}$13

σ_{i_x} is the uncertainty of the 100-s smoothed pseudorange:

${\sigma }_{i_x}^2={\sigma }_{pr\_gnd\_x,i}^2+{\sigma }_{tropo,i}^2+{\sigma }_{pr\_air,i}^2+{\sigma }_{iono\_x,i}^2$ 14

where σ_{pr_gnd_x,i}, σ_tropo,i, σ_{pr_air,i}, and σ_{iono_x,i} are the nominal uncertainties for ground pseudorange noise and multipath, tropospheric delay, airborne pseudorange noise and multipath, and ionospheric delay, respectively. These parameters are obtained from models, as described later.

D_V is the difference between the vertical components of position solutions with 100-s and 30-s smoothed pseudoranges. Here, we used the empirical equation given by Murphy et al. (2010):

$D_V=0.15\times 2\times {\sigma }_{vert}$ 15

2.2.2 Error Models

To evaluate the VPL, each error component in Equation (14) must be modeled. σ_{pr_gnd_x,i} is retrieved from the ground accuracy designator C model, as described by the ICAO (2023):

${\sigma }_{pr\_gnd\_x,i}(El<{35}^{\circ })=\sqrt{\frac{{0.24}^2}{N_{RR}}+{0.04}^2}$ 16

${\sigma }_{pr\_gnd\_x,i}(El\ge {35}^{\circ })=\sqrt{\frac{{(0.15+0.84\cdot e^{-\frac{E{l}_i}{15.5}})}^2}{N_{RR}}+{0.04}^2}$ 17

where N_RR is the number of ground reference receivers, which is assumed to be four in this study.

σ_{tropo, i} is given by the ICAO (2023) as follows:

${\sigma }_{tropo,i}={\sigma }_n{h}_0\frac{{10}^{-6}}{\sqrt{0.002+{\sin }^{2}E{l}_i}}(1-e^{-\frac{\Delta h}{h_0}})$18

where σ_n, h₀, and Δh are the uncertainties in the refractive index, tropospheric scale height, and height of the aircraft above the GBAS reference point (GRP), respectively. σ_n and h₀ are set to 40 and 8.1 km, respectively.

σ_{iono_x, i} is given as follows:

${\sigma }_{iono\_x,i}=F_{pp,i}{\sigma }_{vig}(x_{air}+2\cdot 100\cdot v_{air})$ 19

where σ_vig is the nominal vertical ionosphere gradient parameter broadcast in the GBAS message type 2, x_air is the geometric distance between the aircraft and GRP, v_air is the aircraft speed, and F_{pp, i} is the mapping function from the vertical to slant ionospheric delay. σ_vig is set to 9.5 mm/km (Saito, 2021). The second term in the parentheses accounts for the smoothing filter (100-s) build-up effect resulting from the CCD caused by the ionosphere. The mapping function is defined as a function of the satellite elevation angle as follows:

$F_{pp,i}={(1-{\left(\frac{R_e\cos E{l}_i}{R_e+h_{shell}}\right)}^2)}^{-\frac{1}{2}}$ 20

where R_e and h_shell are the Earth’s radius (6378.1363 km) and ionospheric shell height (350 km), respectively.

Airborne noise and multipath are given as follows:

${\sigma }_{pr\_air,i}^2={\sigma }_{noise,i}^2+{\sigma }_{multipath,i}^2+{\sigma }_{div,i}^2$ 21

where:

${\sigma }_{multipath,i}^2=0.13+0.53\cdot e^{-E{l}_i/10}$ 22

and:

${\sigma }_{noise,i}=0.11\sim 0.15$ 23

Here, σ_{div, i} represents the steady-state value after smoothing filter convergence, which is assumed to be zero (RTCA, 2019).

2.2.3 Airborne Geometry Screening

Airborne geometry screening was performed with the two largest projection coefficients, $\left|\max S_{Apr\_vert}1\right|$ and $\left|\max S_{Apr\_vert}2\right|$:

$\left|\max S_{Ap{r}_{vert}}1\right|+\left|\max S_{Ap{r}_{vert}}2\right|<|limit{S}_{vert}2|$ 24

The threshold limitS_vert 2 is determined based on E_ig and maxE_V, which are aircraft-type-specific parameters that represent the maximum allowed vertical position error for safe landing. In this study, maxE_V was assumed to be 10 m.

The formula for calculating the projection matrix is the same as that for the protection level, except that the smoothing time constant is 30 s instead of 100 s (RTCA, 2019]. The change in smoothing time constant is represented in the weighting matrix as follows:

$W^{-1}=\left[\begin{array}{cccc}{\sigma }_{1-y}^2& 0& \cdots & 0\\ 0& {\sigma }_{2-y}^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }_{N-y}^2\end{array}\right]$ 25

where:

${\sigma }_{i_y}^2={\sigma }_{pr\_gnd\_y,i}^2+{\sigma }_{tropo,i}^2+{\sigma }_{pr\_air,i}^2+{\sigma }_{iono\_y,i}^2$ 26

σ_{pr_gnd_y,i} and _{pr_air, i} were scaled by a factor of $\sqrt{100/30};$ however, this scaling factor is questionable, as the multipath error may not necessarily be represented as Gaussian white noise (Murphy et al., 2023). The uncertainty in tropospheric delay is independent of the smoothing time constant. σ_{iono_y, i} is given as follows:

${\sigma }_{iono\_y,i}=F_{pp,i}{\sigma }_{vig}(x_{air}+2\cdot 30\cdot v_{air})$ 27