GBAS: fundamentals and availability analysis according to σ_vig

2.1 Pseudorange corrections computations

As a GNSS augmentation system, the main function of GBAS is to calculate the differential corrections from the reference GNSS stations and transmit them so that an aircraft equipped with an onboard receiver can perform more accurate positioning. It allows the use of the GNSS with more reliability and safety in precision approach and landing operations, even with poor visibility (EUROCAE, 2013). GBAS use smoothed pseudorange observable (P _CSC) through the Hatch filter:

where S is the GNSS data sampling, τ is the smoothing time constant for the Hatch Filter (100 s in this article), PD _t is raw pseudorange at time t (C/A) code modulated from L1 carrier phase), P CSC ( t − 1 ) is the previous smoothed pseudorange, λ is the wavelength at L1 frequency, and Φ is the carrier phase on L1.

The pseudorange correction is obtained by the difference of the geometric distance calculated by the Euclidean distance between the reference station and a satellite (i) with the pseudorange recorded by the reference station equation (2):

where PRC_CSC(i, j) is a pseudorange correction from satellite i and reference station j; R is a range between satellite i and reference station j; P CSC t is pseudorange smoothing for epoch t; c is the light velocity; and Δt _SV−GPS is the correction satellite clock error.

Then, the PRC_CSC is adjusted by removing the error from the receiver’s clock, which consists of an uncorrelated systematic effect, affecting each reference station differently (PULLEN, 2017). Equation (3) shows the removal of the receiver’s clock error in pseudorange correction:

where PRC_SCA is the pseudorange correction adjusted with the remotion of the clock receiver error; k _i is the weighting factor, with ∑ i ∈ S c k i = 1 ; S _c is the set of satellites tracked by all reference stations. One of the ways to weight the value of k can be by the sin function of the angle of elevation of the satellite (sin θ).

The broadcast pseudorange correction (PRC_TX) for the aircraft equation (4) is determined by the average of the PRC_SCA obtained for satellite i by the M reference stations with corrections valid for this satellite:

where S _i is the set of receivers with valid measures for satellite i.

According to EUROCAE (2013), the rate of change of broadcast pseudorange corrections (range rate corrections – RRC), equation (5) can be calculated by dividing the difference of the PRC_TX in the current time (t) and the previous time (t − 1) by the time difference between the two epochs.

2.2 Performance of the GFS parameters

Moreover, transmitting the pseudorange corrections and their variation rates, GFS also has some parameters that indicate the quality and consistency of the corrections for each reference station.

Among these parameters, B-values equation (6) indicates the consistency of each receiver and satellite measurements that contribute to the generation of the pseudorange correction. They are calculated by the difference between the transmitted pseudorange corrections, and the average of the corrections determined in equation (6), excluding the corrections of the receiver j, whose consistency will be verified (EUROCAE, 2013).

According to Pullen (2017); Lee (2005), the B-value can be interpreted as an error that would occur in the pseudorange correction of satellite i if the measurement of the receiver M were failed. If the receiver j has faults, the correction will be given by the average of the measurements of other receivers, without faults, that track satellite i.

The errors that occur in pseudorange correction are given by the B-values. They occur because, in GBAS, not all systematic errors can be corrected. It is the case of the multipath and noise of GNSS observations. Thus, an uncertainty called σ _{pr_gnd} that describes the contribution of the GBAS Ground Subsystem to errors in the pseudorange correction transmitted to the aircraft. This parameter is calculated using the B-values equations (7 and 8) and transmitted to the aircraft, and is used to determine the position as a component of the uncertainty associated with the pseudorange (EUROCAE, 2013).

where b represents a satellite elevation interval to be considered and N is the number of B-values in the satellite elevation range.

As mentioned earlier, GBAS is based on the principle of spatiotemporal correlation errors, mainly atmospheric errors due to the ionosphere and troposphere. However, in regions where these layers are irregular, the correlation decreases, and there are many errors in the pseudorange corrections. These errors can be mitigated by two uncertainties related to the effects of the ionosphere (σ _iono) and troposphere (σ _tropo) transmitted for the composition of errors associated with pseudorange. According to EUROCAE (2013), the uncertainty associated with the tropospheric error can be obtained by equation (9).

where σ _N is the refractivity uncertainty; h ₀ is the height in meters of the troposphere layer used to determine tropospheric error; θ _i is the elevation angle of satellite i; ∆h is the height difference of the aircraft (meters) and the GBAS reference point. As σ _N and h ₀ are parameters transmitted by GFS, they have to be evaluated prior to the operationalization of GBAS.

In the case of the ionosphere, σ _vig is the parameter transmitted by GFS that describes the variability of the vertical ionospheric gradient. The σ _vig is used to determine σ _iono that represent uncertainty associated with the ionospheric error. Also, the σ _iono parameter is crucial in the calculation of the protection levels and aircraft position.

where F _PP is the slant factor for converting the standard deviation of the ionospheric error from the vertical to satellite–receiver direction, X _aircarft is the distance between aircraft and GBAS reference point (GFS), τ is the time constant used in the smoothing filter (100 s for GBAS Approach Service Type [GAST] C), and V _aircarft is the horizontal speed of the aircraft, given in m/s.

In equatorial and low-latitude regions, the effects of the ionosphere on GNSS signals are one of the main obstacles to a GBAS certification since they represent a threat to accuracy and especially to integrity. According to Pereira et al. (2021), the ionospheric gradient and its uncertainty are parameters of the Ionospheric Threat Model that must be determined to implement GBAS. For such computation, a time series of GNSS data representing the ionospheric irregularity events have to be evaluated to identify its variability. It provides information on the uncertainty of the ionospheric gradient (Figure 2) for the GBAS coverage region.

Figure 2

Representation of variation of the ionospheric gradient for GBAS.

The ionospheric gradient computation can be performed by two methods: Stations pairs and time step. In the studies by Lee et al. (2007); Pereira et al. (2017), details of these methods for determining the ionospheric gradient can be found.

The σ _vig used results from the analysis of the probability density function of the vertical ionospheric gradient calculated for days with nominal ionospheric conditions and days with ionospheric anomalies. The transmitted σ _{vig_inf} is inflated to contemplate the conditions of more significant variations in the ionospheric gradient, for which the mean σ _vig value (μ _vig), the normalized σ _vig value, and an inflation factor (f) are considered.

An example of σ _vig determination (Figure 3) can be found in the study by Lee et al. (2007), in which the sigma overbound method was developed based on the ionospheric gradient probability density function. This method was used to assess extreme or anomalous conditions from the CONUS (CONterminous United States) threat model.

Figure 3

The probability density function of normalized vertical ionospheric gradients on July 2, 2000 (Lee et al., 2007).

2.3 Aircraft position computation

After receiving the corrections and integrity parameters transmitted by the GFS, the aircraft GNSS receiver will calculate the corrected pseudorange at the aircraft level. But, before being applied, some information is checked, such as the received signal strength, latency time valid, and the compatibility of the identical satellites with corrections. The accepted measures are applied in equation (12) of the aircraft’s pseudorange correction (PRC_{corr_air}) (Pullen, 2017):

where P CSC_air ( t ) is the smoothed pseudorange recorded by satellite i at the aircraft receiver at time t before applying correction; PRC_TX is the pseudorange correction transmitted by GFS, being valid for time t _{z_count}; RRC is the rate of change of transmitted pseudorange corrections valid for time t _{z_count}; TC is the differential correction of tropospheric delay; c is the speed of light in a vacuum; ( Δ t ( i ) ) SV corresponds to the satellite clock error applied at time t.

Then, the three-dimensional components of aircraft position (X, Y, and Z) are estimated in the WGS84, as well as the receiver’s clock applying the weighted least squares method, according to EUROCAE (2013); Pullen (2017):

where x ˆ is the vector to be estimated, with the three-dimensional components of position X, Y, and Z and the clock error; A is the design (Jacobian) satellite–receiver matrix, with a dimension of N rows and four columns; W is the weight matrix, with the variance of each measurement by satellite. The Jacobian matrix A is composed of the partial derivatives of the satellite-aircraft line-of-sight vector augmented by 1 for the clock.

where X _si , Y _si , and Z _si are satellite positions; X _a , Y _a , and Z _a are approximated aircraft position; ρ _si is the computed geometric distance from the aircraft to the satellite i (with i = 1, 2, 3, …, n).

However, geocentric coordinates are not most adequate to evaluate the positioning accuracy. According to EUROCAE (2013), Jacobian matrix A can be transformed to the Local Cartesian Frame, North, East, and Down, with origin fixed in the runway on landing threshold point/fictitious threshold point using azimuth (Az) and elevation angle (θ) of each satellite (ranging from 1 to n). In such a case, one has:

The weight matrix W of the satellite measurements consists of a diagonal matrix composed of the inverse of the fault-free variance ( σ P 2 ). This variance term contains all the elements of the ground and aircraft subsystems seen previously:

where airborne error contribution (σ _{pr_air}) is defined in RTCA DO-245A (2004) as the combination of two error sources, receiver noise, and interference and multipath:

where σ multipath ( i ) = 0.13 + 0.53 ⁎ e θ i 10 ; and σ noise ( i ) = a 0 + a 1 ⁎ e θ i θ 0 . The σ _noise term is calculated according to the performance defined by aircraft accuracy designator (AAD) for two levels A and B (Table 1).

The σ _tropo and σ _iono were presented, respectively, in equations (9 and 10). As aircraft was in static mode, the σ iono ( i ) = F pp ( i ) ⋅ σ vig ⋅ X aircarft .

Thus, one can then determine the matrix of W:

According to EUROCAE (2013); Felux et al. (2017), the term (A ^T PA)⁻¹ A ^T P corresponds to the projection matrix, and making S = (A ^T PA)⁻¹ A ^T P, it can determine the uncertainties of the position errors, according to equations (19 and 20). Another way to calculate standard deviations is to express them for errors in the horizontal and vertical components:

where S _horiz. (i) and S _vert. (i) are, respectively, the horizontal and vertical components in the projection matrix determined by equations (21 and 22) (Felux et al., 2017):

The elements S ₁(i), S ₂(i), and S ₃(i) represent, respectively, in the first, second, and third columns of row i of the projection matrix S. The term (θ)_GPA is the angle of the approach path to the runway (3 degrees is generally used).

It can also be derived from the covariance propagation in equation (23), and one obtains the covariance matrix of the estimate parameters (coordinates and clock error):

2.4 Protection levels for precision approach

The protection levels (PL) represent limits for the position errors of each aircraft, based on extremely low probabilities of loss integrity that are required according to GBAS precision approach categories (Pullen, 2017) (Figure 4).

Figure 4

Protection Levels and Alert Limit representation (Lee, 2005).

According to Pullen (2017), the protection levels are built by extrapolating probability distributions limiting the errors expected for the user (aircraft). PL calculation is determined in the vertical (VPL) and horizontal (HPL) components. For the precision approach (PA) service, HPL and VPL are calculated as the maximum value that position errors can achieve, considering hypothesis tests. For the first one, H0 is the nominal case, free of faults, while the alternative H1 assumes a single fault of the reference receiver (Dautermann et al., 2011).

For the nominal case, HPL and VPL are obtained from equations (24 and 25):

where K _ffmd is a multiplication factor for fault detection, the value of which is defined according to the number of reference stations used to generate pseudorange corrections, in this case with 4 receivers K _ffmd = 5.847 m; D _ffmd. and D _vert. are, respectively, the magnitude of the horizontal and vertical projections, considering the difference between the smoothed position by the 30 and 100 second on smoothing pseudorange filter. For the GAST C, the values of D _horiz. and D _vert. are equal to zero, since only the 100 s smoothing filter is used (EUROCAE, 2013).

One additional protection level is defined to cover failures not detected by the ground system. In this case, HPL_PA,H1 and VPL_PA,H1 are obtained by equation (26 and 27, according to EUROCAE (2013); Pullen (2017):

where K _md is the failure detection multiplier, which also depends on the number of reference receivers used: for 2 receivers K _md = 2.935 m; for 3 K _md receivers = 2.898 m; and using 4 the K _md receivers = 2.878 m.

3.1 Generation of GBAS messages

The GBAS simulation carried out met the category of precision approach CAT I or GAST C (GBAS Approach Service Type) that requires the generation and transmission to the aircraft GBAS Message Types 1, 2, and 4. The analysis period of the GBAS simulation was 9 h of data on 03/16/2019, starting at 05 h 08 min 25 s and ending at 14 h 08 min 24 s (Universal Time). All receivers tracked pseudorange (code C/A) and carrier phase on GPS L1 frequency data, with a sampling interval of 0.5 s. The smoothing filter initialization constant was 100 s to obtain smoothed pseudorange, as provided by Standard Recommended and Practices for CAT I or GAST C (GBAS Approach Service Type).

Part of the GGC module configuration is done with the insertion of text files with parameters referring to the geodetic coordinates of the stations; the satellite’s elevation weight function to correct the receiver’s clock error; and in this case, the sen² (θ _i ) function. Other relevant configurations were GAD 2 (ground accuracy designator); minimum elevation angle of the satellites equal to 5°; refractivity of 763; refractivity uncertainty of 20; GCID 1 (GBAS Continuity and Integrity Designator); 42 km for Dmax; kmPOSGPS of 6.0 m; kmCAT1GPS of 5.0 m; vertical alert limit of 10 m; horizontal alert limit of 40 m; GPA (glide path approach) equal to 3°; 49 m for threshold crossing height; and 105 m course width. A more detailed description of all these parameters can be found at EUROCAE (2013).

Another parameter configured in the GGC module is σ _vig, which is part of the content of Type 2 messages. Although it is not a variable used in pseudorange correction generation algorithms and ground subsystem integrity parameters, σ _vig is transmitted as a parameter of the region covered by the airport. For this work, five scenarios with σ _vig values were evaluated: 4, 8, 12, 16, and 24 mm/km. As stated before, such selection aims to evaluate the impact of the ionospheric gradients on the availability of GBAS.

3.2 Positioning and integrity parameters

After the generation of GBAS Message content by the GGC module, the simulated aircraft data were used in the GNSS Solution module to obtain estimates of coordinates, protection levels, position errors, and other indicators inherent to GNSS positioning.

The GPS receiver representing the aircraft (simulated) was in Guaratinguetá – SP (Figure 5b) and had a similar setting to the receivers of the reference stations mentioned in Section 4.1. The AAD parameter was selected as level A. The AAD parameter is used to determine the pseudorange uncertainty calculated by the aircraft (σ _{pr_air}) and is based on the multipath and noise at the aircraft’s receiver.

As the experiment was carried out in static mode and as the position of the simulated aircraft is known, they were inserted to analyze the horizontal and vertical positioning errors (HPE and VPE).

4.1 Results in the position domain

Now, GBAS simulations at the aircraft level will be carried out, with the estimates in the positions and their errors and the integrity indicators and performance analysis.

Thus, Figure 10 shows the time series of the VPE and VPL calculated for the GAST C precision approach. In this scenario, a σ _vig value of 4 mm/km was applied.

Figure 10

Performance of positioning errors and protection levels for GBAS simulation considering σ _vig of 4 mm/km.

In the horizontal component, the HPE was lower than the HPL during the whole period analyzed for all scenarios of the σ _vig. The horizontal alert limit is 40 m and was not exceeded in any evaluated scenarios.

On the other hand, in the vertical component, 3.5% of the time, the VPE was higher than the VPL, indicating the occurrence of Misleading Information (MI). In this case, the system is declared available because VPE and VPL did not exceed the vertical alert limit. The protection level represents the maximum estimated error in positioning for the precision approach to occur safely.

The second scenario evaluated considered the σ _vig value equal to 8 mm/km. In Figure 11, it is possible to observe the performance of the horizontal and vertical components for this scenario.

Figure 11

Performance of positioning errors and protection levels for GBAS simulation considering σ _vig of 8 mm/km.

The VPE had only 9 epochs between 6 am and 7 am, where there was a slight superiority of the VPE over VPL. The σ _vig of 8 mm/km was the scenario with the highest availability of 99.9% and the lowest occurrence of MI events. However, there was an increase in the protection levels compared to the previous scenario.

For the 12 mm/km σ_vig scenario (Figure 12), HPE and VPE were lower than the HPL and VPL. However, there was an increase in the values of the protection levels, and 2.2% times, the VPL exceeded the alert limit that indicates loss of availability for the CAT-I approach.

Figure 12

Performance of positioning errors and protection levels for GBAS simulation considering σ _vig of 12 mm/km.

This condition of the vertical protection levels exceeding the vertical alert limit can also be seen in the other tested scenarios with a σ _vig of 16 mm/km and 20 mm/km, as shown in Figures 13 and 14. For σ _vig of 16 mm/km was 17.3% time and σ _vig of 20 mm/km was 38.4%, the VPL was higher than VAL. Another interpretation of protection levels is that they represent the maximum estimated error that position errors can achieve. Therefore, Figures 13 and 14 show the reduced availability of GBAS due to the estimated error (VPL) being more extensive than the VAL. This reduction decreases with increasing σ _vig, although the position error is still below the alert limit.

Figure 13

Performance of positioning errors and protection levels for GBAS simulation considering σ _vig of 16 mm/km.

Figure 14

Performance of positioning errors and protection levels for GBAS simulation considering σ _vig of 20 mm/km.

A critical factor observed was the presence of peaks and discontinuities in the protection levels, increasing and decreasing the PL values. This occurrence coincides with the change in satellite geometry, as shown in Figure 15, which is related to PDOP (precision dilution of position) and the number of satellites visible (NSV) at each positioning epoch.

Figure 15

Comparison of VPL according to the adopted σ _vig and the influence of satellite geometry.

The exclusion of satellites can be due to several reasons such as no corrections, decorrelation of ephemeris, and the restart of the pseudorange smoothing process due to cycle slip. However, a detailed analysis was not carried out about each occurrence, only the identification of the behavior of the VPL with the variations of PDOP and NSV.

An anomalous ionospheric gradient occurrence may be one of the causes for the reduction of satellites observed in Figure 15, so Table 2 shows the impact of σ _vig on the residuals errors uncertainty from the ionosphere (σ _iono) calculated according to equation 10 and the approximate distance of 31 km between the reference stations and simulated aircraft, as well as elevation angle: 5, 10, 15 and 20°. As the aircraft is in static mode, the speed is at zero, and the σ _iono is the result of the multiplication of the terms F _pp, σ _vig, and X _aircraft.

The increase in residual uncertainty due to the ionosphere can be observed as the σ _vig, and the elevation angle of the satellite increases, representing, therefore, a multiplicative factor of uncertainty as the pseudorange measure to be applied in the positioning of the aircraft and the calculation of the protection levels as shown in Figure 14.

It is important to emphasize that in the future, the tendency is for this degradation caused by the geometry and availability of satellites to be reduced, as there will be a more significant number of satellites, as well as the possibility of new GNSS signals that will provide linear combinations that aim to minimize effects systematic as well as the variation of the ionospheric gradient.