Abstract
Large ionospheric gradients acting between a Ground Based Augmentation System (GBAS) reference station and an aircraft on approach could lead to hazardous position errors if undetected. Current GBAS stations provide solutions against this threat that rely on the use of “worst-case” conservative threat models, which could limit the availability of the system.
This paper presents a methodology capable of detecting ionospheric gradients in real time and estimating the actual threat model parameters based on a network of dual-frequency and multi-constellation GNSS monitoring stations. First, we evaluate the performance of our algorithm with synthetic gradients that are simulated over the nominal measurements recorded by a reference network in Alaska. Afterwards, we also assess it with one real ionospheric gradient measured by the same network.
Results with both simulated gradients and a real gradient show the potential to support GBAS by detecting and estimating these gradients instead of always using “worst-case” models.
1 INTRODUCTION
The Ground Based Augmentation System (GBAS) is a local-area, airport-based augmentation of Global Navigation Satellite Systems (GNSS). Its main purpose is enhancing GNSS performance in terms of integrity, continuity, accuracy, and availability. A GBAS reference station broadcasts differential corrections along with integrity parameters. The differential corrections enable an aircraft approaching an airport to correct the navigation signals from the satellites by removing the spatially correlated errors between the ground station and the aircraft. Thus, the aircraft is able to improve the accuracy of its position estimation. Additionally, the integrity parameters enable the airborne system to calculate bounds on the residual position errors and ensure safety of the operation.
Single-frequency (L1) single-constellation (GPS) differential corrections, as they are provided today by GBAS, enable airborne users to correct most of the GPS ranging errors, especially the ones caused by the ionospheric delay. Nevertheless, nominal residual errors still remain and are overbounded by the so-called protection levels (PLs). The protection levels are calculated by the avionics using error models and the integrity parameters received from the ground station. Then, they are compared to the alert limits (ALs), the maximum allowable bounds, to determine safety of the operation for each user. If the PLs exceed the ALs, GBAS is unavailable. The component of the PLs that describes and overbounds the remaining residual errors due to the nominal ionospheric decorrelation between the reference station and the user is called 𝜎𝑣𝑖𝑔 (Chang et al., 2019; Lee et al., 2007; Mayer et al., 2009). This is an integrity parameter broadcast by the ground station that ensures the safety of the operation only as long as the state of the ionosphere is nominal.
However, ionospheric anomalies, like large ionospheric gradients, might produce a significant difference between the ionospheric error experienced by the GBAS reference station and the aircraft on approach. This ionospheric delay difference could lead to hazardous position errors if undetected, since it is not correctly overbounded by the integrity parameters and therefore results in misleading integrity information. For that reason, the GBAS Approach Service Types (GAST) C (RTCA, 2008) and D (RTCA, 2017) provide solutions to mitigate the ionospheric gradient threat, but the methods employed still face challenges by limiting availability in certain cases.
A GAST-C ground station, developed to support CAT-I operations, inflates the integrity parameters in order to exclude potentially usable satellite geometries that could produce unacceptably large position errors if affected by a “worst-case” gradient in a process called “geometry screening” (Lee et al., 2011; Seo et al., 2012). This “worst-case” gradient is modeled as a wave front of a certain magnitude (represented by the slope and width) moving with constant speed and direction (Figure 1). Furthermore, the “worst-case” gradient belongs to a threat model that is derived based on the worst-ever-experienced ionospheric gradients measured in the relevant region and defines a range of values for the gradient parameters that combined could harm GBAS. In Table 1, we show the main parameters of the threat model for the CONUS region (Datta-Barua et al., 2010), Germany (Mayer et al., 2009), and Brazil (Yoon et al., 2017). However, the main problem is that this methodology assumes that the “worst-case” gradient is always present, which is very unlikely. This strong assumption, together with the high values of the threat model in regions with severe ionospheric conditions like Brazil, can result in satellite geometries being excluded even in days where no perturbation is present, leading to a loss of availability (Yoon et al., 2019). Moreover, the threat model is derived based on historical data; thus, this procedure cannot protect users against a larger gradient that might occur at some point in time in the future.
GBAS ionospheric threat model
GBAS ionospheric gradient parameters
The GAST-D concept, developed to support CAT-II/III operations, contains additional monitors to mitigate this threat in mid-latitudes (Pullen et al., 2017). The main problem in this case is that the baseline concept of GAST D assumes that the prior probability of occurrence of an anomalous ionospheric gradient is one. This assumption, as previously mentioned, is conservative and causes the monitors designed to protect integrity in GAST D to be very sensitive and trigger false alerts. Recent studies (Yoon et al., 2019; Yoon et al., 2020) already suggest reducing this prior probability of occurrence to 10−3, which would significantly relax the monitoring thresholds making the implementation of GAST D easier in regions where its availability is degraded. However, the derivation of the prior probabilities of occurrence of an anomalous ionospheric gradient relies on a statistical approach that uses historical data. Thus, the lack of sufficient historical data collected during large ionospheric events in certain regions might limit the use of this statistical approach. Moreover, one of the monitors designed to protect GAST D against ionospheric gradients (the Ionospheric Gradient Monitor, IGM (Khanafseh et al., 2012)) has associated siting constraints of the GBAS reference receivers to ensure that all potential gradients can be detected. This makes it difficult to deploy in certain airports with space constraints. While trying to adapt the GAST D concept to other regions with more severe ionospheric conditions than those present in mid-latitudes, the previously mentioned issues could have a negative impact on GBAS availability (ICAO, 2017).
In previous work (Caamano et al., 2017), we proposed a real-time ionospheric monitoring approach that could reduce the conservative assumptions currently applied in GAST C and D and thus improve their availability under active ionospheric conditions. Moreover, we described how the monitor could work under nominal conditions and evaluated the impact of noise and multipath on the monitoring concept.
In this paper, we build on Caamano et al. (2017) to propose an algorithm that addresses the detection of ionospheric gradients in real time and estimates the gradient parameters in near real time. Furthermore, we evaluate our algorithm with simulated ionospheric gradients, assessing the differences between the known simulated gradient parameters and the parameters estimated by our algorithm. These simulations correspond to regional ionospheric disturbances (i.e., beyond the GBAS local scale), which are based on actual perturbations. Additionally, we also evaluate our algorithm with a real ionospheric gradient measured by monitoring stations in Alaska to show the differences between using simulated gradients and real gradients.
This paper is divided into six main sections: Section 2 reviews the real-time ionospheric monitoring concept, Section 3 introduces the data used for our evaluations, Section 4 describes our algorithm for detection of the gradients and estimation of the gradient parameters, Section 5 introduces the simulation setup and the ionospheric gradient simulator, Section 6 presents the results of the evaluation of the monitor with simulated gradients and a real gradient observed by a reference network, and Section 7 discusses the applicability of the proposed methodology to GBAS.
2 DUAL-FREQUENCY MONITORING CONCEPT
Our proposed monitor provides protection to all GBAS users by using a wide area network of dual-frequency and multi-constellation GNSS stations situated in carefully surveyed locations. The network would consist of mainly the GBAS stations installed in the relevant region that would need dual-frequency capable receivers or have a second dual-frequency receiver installed for monitoring purposes only. Since it is foreseen that the GBAS stations will be installed principally at important airports, the network might need other external dual-frequency monitoring stations to enhance the coverage of the GBAS stations while monitoring for gradients. These external monitoring stations could be newly installed and/or already existing dual-frequency reference stations like (e.g.) Satellite Based Augmentation System (SBAS) stations. Note that the quality of the measurements provided by existing reference stations must be sufficient to guarantee integrity in order to be considered as part of the network. The dual-frequency measurements coming from both types of stations (GBAS and external) would be used to estimate the ionospheric delay reliably in order to detect ionospheric anomalies in real-time. Furthermore, the utilization of multiple constellations would provide improved sampling of the ionosphere. Additional ionospheric measurements coming from other constellations could also be used to support the single-frequency and single-constellation GBAS service types C and D by extending their knowledge of the ionospheric state beyond what would be available from a single constellation. The applicability of the dual-frequency monitoring concept to GAST-C and GAST-D is further discussed in Section 7.
The functionality of the proposed monitoring network is described in the following:
The processing component of each of the stations receives GNSS dual-frequency and multi-constellation code and carrier-phase measurements and calculates the slant ionospheric delay at current time 𝑡 as in Equation (1) (Misra & Enge, 2006):
1where
is the carrier-phase measurement in meters for frequency 𝑓1 (L1/E1) and the carrier-phase measurement in meters for frequency 𝑓2 (L5/E5a). A single satellite from the GPS or other satellite constellations is represented by 𝑗 and identified with its system name and its number within its system (e.g., G03 is satellite number 3 in the GPS constellation). The station that calculates the slant ionospheric delay is 𝑟, which is identified with a certain station name. Moreover, each station 𝑟 belongs to the monitoring network of stations denoted as . Note that the slant ionospheric delay is calculated at so-called Ionospheric Pierce Points (IPPs), the points of intersection between the satellite-receiver line of sight and the ionosphere modeled as a “thin shell” located at 350 km above the Earth’s surface (Klobuchar, 1987). Since in this work we consider only the slant ionospheric delay calculated for frequency (L1/E1), we omit the frequency subscripts in the following.The rate of change of the estimated ionospheric delay corresponding to each satellite 𝑗 and each station 𝑟 is compared with a predefined threshold derived with real measurements from this station 𝑟 and integrity requirements. This threshold decides whether there is a perturbation affecting this satellite-station pair in real time. In Sections 4.1 and 4.2, we give more details about the threshold derivation and the detection algorithm, respectively.
The detection information from each of the stations observing a given satellite is shared within the network in real time. Here, the network distinguishes between two possible cases: i) several monitoring stations have detected the same ionospheric gradient, and ii) none of the monitoring stations have detected anomalous ionospheric gradients. In the case that several monitoring stations have detected the same ionospheric gradient, a central processor estimates its parameters. This gradient parameter estimation process is done per satellite, and it requires that at least three stations have detected the gradient (see Section 4.2.2). If fewer than three stations have detected the gradient, the network is not able to estimate the actual gradient parameters, and it triggers a “Warning” to indicate that the GBAS stations should use “worst-case” assumptions. If none of the monitoring stations have detected any anomalous ionospheric gradients, the network calculates the largest gradient that could be affecting the supported GBAS stations without being detected.
Each GBAS station would then be responsible of using the network information to support the already existing equipment in covering the GBAS approaches at each airport.
3 DATA
3.1 Data description
In this section, we describe the data that we used to evaluate the detection and estimation algorithm proposed in this work. As the monitoring network, we selected five reference stations situated in Alaska and depicted in Figure 2 with the coordinates of each station in Table 2. Public data is available for this network at a 1 Hz sampling rate both for L1 and L2 frequencies and GPS satellites (https://www.unavco.org).
Station coordinates.
Station coordinates
Note that, due to the limited availability of data recorded with L1/L5 frequencies and other constellations (e.g., Galileo) during active ionospheric conditions, we used for our tests real measurements on L1 and L2 frequencies and only the GPS constellation. It is expected that the performance of our methodology evaluated with measurements from multiple satellite constellations broadcasting L1/E1 and L5/E5a will be better than the performance achieved in the present paper with only GPS and L1/L2 frequencies. The reason why this is expected is that measurements recorded on L5 are typically less noisy than on L2, and more measurements available from more satellites will provide better sampling of the ionosphere.
3.2 Date selection
The dates selected from all the data available attend to two different purposes: i) the study of a real anomalous ionospheric gradient measured by the network depicted in Figure 2, and ii) the derivation of the monitoring threshold for each of the stations in the network.
As an “active” ionospheric day, we selected the geomagnetic storm that occurred on the 17th of March of 2015, St. Patrick’s Day, which has been extensively studied in the literature (Béniguel et al., 2017; Jacobsen & Andalsvik, 2016). For the threshold derivation, we selected manually ten “quiet” days prior to the “St. Patrick’s Day Storm.”
We determined both the active and the nominal days based on an ionospheric activity index, the Along Arc TEC Rate (AATR) (Juan et al., 2018b), calculated for one of the stations under study, av17. The specific AATR values for this station can be found in Figure 3. As can be observed, the AATR for the day of year 76, i.e., St. Patrick’s Day, is above 1 TECU/min, which can be considered as very high activity (Juan et al., 2018b). Furthermore, “quiet” ionospheric conditions refer to AATR values below 0.2 TECU/min. As a last step, we performed a visual inspection of the data recorded in the stations under study, ac59, av17, av16, av20, and av01 during the days considered “nominal.” During this process, the days with corrupted or missing measurements were discarded.
Along Arc TEC Rate (AATR) values for reference station “av17” during 30 days of year 2015. The AATR values are calculated every 5 minutes. The threshold of 0.2 TECU/min represents the value below which the ionosphere is considered “quiet” or nominal and the threshold of 1 TECU/min is the value above which the ionosphere is considered very “active” or abnormal.
According to this criteria, the nominal days used to derive the detection thresholds for the stations av16, av17, ac59, and av20 were: 52, 56, 57, 58, 63, 64, 68, 69, 72, and 73 of year 2015. In the case of station av01, the measurements presented certain discontinuities on the days 63 and 64. Thus, we substituted these corrupted days with days 67 and 70 of year 2015.
4 METHODS
With the aim of monitoring for possible gradients, we developed an algorithm that is able to detect anomalous ionospheric gradients in real time and estimate the gradient parameters like slope, width, speed, and direction in near real time. However, before the real-time operation of the algorithm, a preprocessing part is needed in order to derive detection thresholds while taking into account the characteristics of each of the stations in the monitoring network. In the following, we explain both parts: the derivation of the monitoring thresholds and the real-time operation of the algorithm.
4.1 Derivation of the monitoring thresholds
A correct monitoring threshold derivation is the key part in the algorithmic chain since this threshold determines whether we measure an anomalous ionospheric gradient or not. Therefore, it is important to study the expected performance of each of the monitoring stations inside the network in days when the state of the ionosphere was considered nominal. For that purpose, we select manually for each of the monitoring stations ten days of dual-frequency GNSS measurements recorded under nominal ionospheric conditions as explained in Section 3.2.
Note that the amount of nominal data that is selected for the monitoring threshold derivation should be large enough to cover all possible satellite geometries used while avoiding false alerts due to environmental features (e.g., multipath).
Nevertheless, even in nominal conditions, the carrierphase measurements might contain cycle slips that could result in false gradient detections; thus, it is necessary to “clean” the data before calculating the threshold.
This data cleaning process consists of the following two stages:
A cycle slip detector explained in Sanz et al. (2013) searches for jumps in the slant ionospheric delay estimation
computed as in Equation (1). Combining measurements from two different carrier phases removes the geometry influence, including clocks and all non-dispersive effects. In non-scintillation conditions, the slant ionospheric delays experience smooth changes between consecutive epochs, even if the receivers have a sampling rate of 30 s. This cycle slip detector predicts by performing a polynomial fitting of the last 𝑁𝐼 samples, where 𝑁𝐼 is the size of the window used in the polynomial fitting (e.g., 10 samples in Sanz et al., 2013). Considering a sampling rate of 1 Hz, a cycle slip is declared when the actual and its predicted value differ more than 3.18 cm (see Sanz et al., 2013 for an explanation of the cycle slip detector threshold value). However, since this test signal is still driven by frequency-dependent effects in regions where “nominal conditions” include some activity in the ionosphere, some less evident cycle slips could remain undetected. Examples of this problem can be found in Juan et al. (2018a) and Juan et al. (2017). As a consequence, the computed thresholds using only this cycle slip detector could suffer from an increase of their values degrading the detection performance of the monitoring network.For this reason, in this work, a manual check is applied after the cycle slip detector to remove possible remaining cycle slips and outliers. Another possibility in this step could be to use a second more complex and sensitive cycle slip detector, which would be preferable when working with greater amounts of data. One example of these more sensitive cycle slip detectors can be found in Juan et al. (2017), where the authors propose a method similar to the GBAS acceleration-ramp-step monitor applied to the ionospheric-free combination of carrier phases (i.e., unaffected by ionosphere). Using the knowledge of the precise coordinates of the receivers on the ground, the precise orbits, and the clock corrections, they remove most of the physical effects present on the ionospheric-free combination. The resulting values are completely independent of the ionospheric effects, are very accurate, and can be used to detect these smaller cycle slips.
After data cleaning, we derive the detection threshold for each of the monitoring stations 𝑟. First, we compute the test statistic 
for each epoch 𝑡, each satellite 𝑗, and each station 𝑟, as the first derivative or rate of the cleaned ionospheric delay in order to remove the unresolved ambiguities in the carrier-phase measurements.
2
Here, 
is the ionospheric delay estimation for station 𝑟, satellite 𝑗, and epoch 𝑡, 
is the ionospheric delay estimation for station 𝑟, satellite 𝑗, and the previous epoch, and 𝛿𝑡 is the time difference between two consecutive epochs. Note that any cycle slip would produce a high rate value if it is not detected and removed from the slant ionospheric delay estimation.
After that, we sort these 
computed for all the satellites with elevations 
, all epochs 𝑡 during ten days, and one of the stations 𝑟 in elevation bins. The size of these elevation bins depends also on the elevation: 2◦ between 5◦ and 25◦ of elevation, 5◦ between 25◦ and 50◦ of elevation, and 10◦ between 50◦ and 90◦ of elevation. This different binning size is used to account for the fewer amount of samples available from high elevation satellites. Then, given an acceptable false alert probability, 𝑃𝑓𝑎, a threshold for each station can be defined as
3
where 𝜎𝐼𝑡𝑒𝑠𝑡𝑟 (𝜃𝑚) is the standard deviation of the 𝐼𝑡𝑒𝑠𝑡𝑟(𝜃𝑚) distribution composed of the 
samples for all satellites and epochs arranged into elevation bins, with 𝜃𝑚 and 𝑚 = 1, 2,.., 𝑀 representing the elevation bins. Moreover, 𝑘𝑓𝑎 is the false alert multiplier computed from the inverse of the standard normal cumulative distribution, 
.
However, this methodology only applies if 𝐼𝑡𝑒𝑠𝑡𝑟(𝜃𝑚) is Gaussian. In the case this probability distribution presents non-Gaussian behavior, we calculate a Gaussian overbound of the tails of 𝐼𝑡𝑒𝑠𝑡𝑟(𝜃𝑚) per bin following the approach in Shively and Braff (2000) and Xiong (2015).
For that purpose, first, we normalize the data by substracting from each 
sample the mean computed with the samples inside each elevation bin, 𝜇𝐼𝑡𝑒𝑠𝑡𝑟 (𝜃𝑚), and dividing the result by the standard deviation of each bin, 𝜎𝐼𝑡𝑒𝑠𝑡𝑟 (𝜃𝑚). The probability distribution composed of these normalized data samples is denoted as 𝐼𝑡𝑒𝑠𝑡𝑟,𝑛𝑜𝑟𝑚(𝜃𝑚).
Then, for each elevation bin, we calculate the cumulative distribution function (CDF) of 𝐼𝑡𝑒𝑠𝑡𝑟,𝑛𝑜𝑟𝑚(𝜃𝑚) and its complement (1-CDF). Also, the CDF and its complement are calculated for a zero-mean Gaussian distribution 
with an inflated standard deviation 𝑖𝑚 and 𝑚 = 1, 2, .., 𝑀 representing the elevation bins. These two functions are used to overbound the tails of 𝐼𝑡𝑒𝑠𝑡𝑟,𝑛𝑜𝑟𝑚(𝜃𝑚) per bin. Therefore, the standard deviation inflation factor 𝑖𝑚 is calculated such as 𝐶𝐷𝐹(𝐼𝑡𝑒𝑠𝑡𝑟,𝑛𝑜𝑟𝑚(𝜃𝑚)) < 𝐶𝐷𝐹(𝐺𝑚) and 1−𝐶𝐷𝐹(𝐼𝑡𝑒𝑠𝑡𝑟,𝑛𝑜𝑟𝑚(𝜃𝑚)) < 1 − 𝐶𝐷𝐹(𝐺𝑚) per elevation bin. Examples of the tail overbounding process can be found in Shively and Braff (2000) and Xiong (2015).
The detection threshold for each elevation bin at station 𝑟 is then defined as
4
Figure 4 shows an example of different test statistic values (in blue) and different values for the detection threshold depending on the different probabilities of false alert considered. The data used in this example corresponds to all visible GPS satellites for the ten days specified for station “av17” in Section 3.2. In this case, the probability that a certain 
value exceeds the threshold is smaller than the 𝑃𝑓𝑎 utilized to calculate the 𝑘𝑓𝑎 multiplier due to the overbounding performed, which means that the requirement of false alert probability of the monitor is satisfied. Furthermore, we can observe that the threshold values depend greatly on satellite elevation due to the higher effects of noise and multipath on low elevation satellites, which also make the inflation factors required for overbounding larger (above 1.5 for satellite elevations below 30◦ but below 1.2 for satellite elevations above 35◦). Concerning the different probabilities of false alert, when we allow a larger number of false alerts, the threshold gets more restrictive, and when we allow fewer false alerts, the threshold is relaxed, but then some gradients could be missed. Therefore, we selected a 𝑃𝑓𝑎 of 10−6 as an acceptable compromise for our tests with real measurements. It should be noted that a false alert in a certain station of the network does not automatically lead to a loss of continuity but instead requires assuming the “worst-case” ionospheric threat model parameters on that area of the network. Thus, the selected 𝑃𝑓𝑎 is likely not optimal, and the optimization of its value will be investigated as part of future work.
Example of different detection thresholds for different𝑓𝑎 considered. The blue dots show the 𝐼𝑡𝑒𝑠𝑡𝑎𝑣17(𝜃𝑚) distribution calculated for 10 nominal days of year 2015 for all elevation bins at station “av17.” All elevation bins contain more than 350,000 samples with the exception of the 80◦ to 90◦ bin, which contains 121,493. The black continuous curve shows the threshold per elevation bin for 𝑃𝑓𝑎 = 10−6 and the red dashed curve the threshold per elevation bin for𝑓𝑎 = 10−9.
4.2 Real-time operation of the algorithm
The real-time operation of the algorithm can also be divided into two main parts: the detection step and the estimation step.
4.2.1 Detection step
As previously stated in Section 2, the detection step is performed individually per station 𝑟, satellite 𝑗, and epoch 𝑡, monitoring for ionospheric anomalies in the rate of the estimated slant ionospheric delays. The detection algorithm receives the GNSS dual-frequency and multi-constellation carrier-phase measurements and calculates the slant ionospheric delay as in Equation (1). Then, these measurements undergo a cleaning process intended to remove possible cycle slips in the data. Once the data is cleaned, we compute the test statistic 
(Equation (2)). Given a predefined threshold curve for a certain monitoring station 𝑟 and considering its value for the elevation of a satellite 𝑗 at a certain epoch 𝑡, 
, the condition for detecting the ionospheric gradients is
5
The algorithm, presented in Algorithm 1, outputs detection information that is shared over the network in real time. This detection information consists of the 
values, a detection flag 
that becomes one when a gradient is detected, and a signal flag 
that indicates that valid measurements in both frequencies exist at the current time 𝑡.
Detection algorithm for a single monitor station
Note that the detection capabilities of each of the individual reference stations are just as important as those of the network. In the case that the reference stations forming the network are so far apart that the use of multiple constellations is not sufficient to monitor the whole area, ionospheric irregularities such as small plasma bubbles might not be detected, as mentioned in Yoon et al. (2017). This issue should be carefully investigated when assessing the detection capabilities of the network, and a criteria for a maximum distance between monitoring stations should be established. However, the study of the “desired” value of the station separation for ideal detection of these small irregularities is not addressed in this paper and is part of future work.
Furthermore, in the real-time operation of the algorithm, we perform only stage 1 of the data cleaning process presented in Section 4.1. The main reason is that this cycle slip detector can be implemented in real time, whereas a manual check or another more complex cycle slip detector (such as the one proposed in Juan et al., 2017) is not yet real-time capable, although studies on how to adapt these more complex cycle slip detectors to real time are ongoing. Since we designed our monitoring thresholds using both stages in order to guarantee integrity, in the very rare occasions where an undetected cycle slip occurs, our monitoring thresholds will trigger false alerts. In these cases, the monitoring network tries to estimate the ionospheric gradient parameters as if the false alerts were real gradients. In Section 4.2.2, we explain this process and discuss the behavior of the algorithm during false alert events. Moreover, our monitoring network would consist of either GBAS stations with very high quality and reliable measurements or external stations that are also expected to provide high-quality measurements, therefore minimizing the likelihood of this problem.
4.2.2 Estimation step
The estimation step of the algorithm collects in a central processor the detection information shared in the network per station 𝑟 in real time. First, this central processor groups the information coming from different stations for the same satellite. Then, it estimates the gradient parameters explained in Section 1 per satellite. All the formulas in the following are therefore expressed for a single satellite 𝑗, and thus we omit the 𝑗 superscript.
Speed and direction of the gradient
In order to calculate the speed of the ionospheric gradient, we need to track its spatial evolution with time or, in other words, the time delay 𝜏𝑟 between detections in two stations that are spatially separated, a station 𝑟, and a station that we choose as reference. In principle, we select as the reference station the first station that detects the gradient, and we distinguish it by substituting the subscript 𝑟 with zero. Additionally, we assume that the ionospheric gradient is local, maintaining its characteristics of magnitude (slope and width) and propagation (speed and direction) constant only over a certain time and distance. This assumption, which we refer to as the “locality principle” in this work, is verified at the end of this subsection. Furthermore, we also consider that the ionospheric disturbance propagates as a planar wave that moves with a certain speed and impacts the different IPPs corresponding to the same satellite and different stations 𝑟 at different times. Note, that the assumption of a planar wave propagation is reasonable for station baselines of a few kilometers, as considered in this work (Figure 2).
However, since several ionospheric gradients can occur in a short period of time at each ground station, it is necessary to identify the same gradient occurring at different ground stations. Under the assumption that the gradient maintains its characteristics of magnitude during a certain period of time, this means that we need to identify the gradient with the same shape at different stations. To this end, we compute the cross-correlation between the test statistic values calculated in different stations for the same satellite.
Nevertheless, due to the real-time constraint in GBAS, it is necessary to perform the cross-correlation between the test statistics as fast as possible instead of using the complete set of measurements over a day. Thus, we calculate the cross-correlation between 𝐁_𝐈𝐭𝐞𝐬𝐭𝐫, a buffer containing a history of test statistic values until current time 𝑡 for any of the stations 𝑟 ≠ 0 ∈ 
, and 𝐁_𝐈𝐭𝐞𝐬𝐭𝟎, the buffer for the station of reference.
The two buffers involved in the cross-correlation process are defined as
6
where 𝑡𝑑0 is the epoch when the station of reference detected the anomalous ionospheric gradient for the first time, 𝑁𝐵 is the size of both buffers, and 𝑑𝑤 is a window of time designed to capture the part of the test statistic that starts to increase when a gradient begins but is still not sufficiently large to trigger the detection thresholds. In this work, we chose a size of 30 seconds for 𝑑𝑤 taking into consideration the characteristics of the data that we had available. As can be observed in Equation (6), the buffers belonging to all non-reference stations begin to fill at the same moment the station of reference detected an anomalous gradient, 𝑡𝑑0. Therefore, all buffers are initialized with the test statistic values from the stations to which they belong from time 𝑡𝑑0 −𝑑𝑤 to 𝑡𝑑0 and continue to receive and store the test statistic values until current time 𝑡, 𝐼𝑡𝑒𝑠𝑡𝑟(𝑡). Thus, all buffers have the same length. The maximum length of these buffers is another design parameter, and its determination should take into account the range of values of the expected anomalous gradient parameters.
We describe the cross-correlation process in more detail in Algorithm 2. For the times when 𝑑𝑡𝑐𝑟(𝑡) is equal to one, the algorithm checks if there is any sample of the test statistics missing in any of the two buffers involved in the cross-correlation process. These missing data are often due to a cycle slip, a data gap, or when the satellite is no longer visible. Since we cannot ensure integrity if we have a gap in the data, if this occurs, a “Warning” is declared to indicate to the network that there has been a detection, but the gradient parameters cannot be determined. In this case, the network should use the “worst-case” threat model parameters on that area. After this initial check, Algorithm 2 computes the cross-correlation between the two buffers, 𝐁_𝐈𝐭𝐞𝐬𝐭𝟎 of the reference station and 𝐁_𝐈𝐭𝐞𝐬𝐭𝐫,as
7
Real-time cross-correlation algorithm for a single satellite
For each time 𝑡 in which the cross-correlation is computed, we find the maximum of the cross-correlation and the index of the cross-correlation vector where the maximum occurs 𝑝𝑚𝑎𝑥(𝑡). Then, the time delay between the two stations at time 𝑡 is computed as
8
Once the time delay is found, the cross-correlation coefficient, 𝛼𝑟(𝑡), is calculated between 𝐁𝟎 and 𝐁𝐫, two buffers containing the relevant parts of the signal from 𝐁_𝐈𝐭𝐞𝐬𝐭𝟎 and 𝐁_𝐈𝐭𝐞𝐬𝐭𝐫. This cross-correlation coefficient, also known as the Pearson correlation coefficient, is determined by
9
where 𝑐𝑜𝑣 represents the covariance of the two vectors and 𝜎 the standard deviation.
A cross-correlation coefficient of one between the tests statistics computed at two different stations means that the perturbation is the same, but delayed by a certain time interval, which we use to estimate the propagation parameters. A cross-correlation coefficient of zero means that the perturbations at these two stations are totally different. Therefore, these cross-correlation coefficients 𝛼𝑟(𝑡) are compared in each epoch with a minimum value 𝛼𝑚𝑖𝑛, below which we do not consider the perturbation occurring at the station 𝑟 to be the same as the one occurring at the station of reference. In this work, we use for 𝛼𝑚𝑖𝑛 a value of 0.5, which is a common value used in signal processing. Until 𝛼𝑟(𝑡) reaches the value chosen for 𝛼𝑚𝑖𝑛, the network does not know if both stations detected the same anomalous ionospheric gradient. Therefore, it declares a “Warning” to indicate to the GBAS stations in that area that they should use the “worst-case” threat model parameters.
When 𝛼𝑟(𝑡) is above 𝛼𝑚𝑖𝑛, the network considers that the anomalous gradient measured at both stations, 𝑟 and the reference, are the same. However, when the difference between current epoch 𝑡 and time of detection at a station different from the reference, 𝑡𝑑𝑟, is low, the real-time cross-correlation algorithm finds poor correlations due to the noise and multipath present in the 𝐼𝑡𝑒𝑠𝑡𝑟(𝑡) samples and the very limited amount of data corresponding to the gradient stored in the respective buffers. Therefore, the algorithm needs a certain time to converge. We consider that the algorithm converges when the difference between the cross-correlation coefficients 𝛼𝑟 at the current epoch and at the previous epoch is below 1𝑥10−2 for the last 𝑁𝑐 samples (e.g., three samples). The convergence condition is summarized in Equation (10).
10
The amount of time that the algorithm needs to converge depends on the characteristics of the ionospheric gradient and the level of noise and multipath present on the measurements, and it is further discussed in the results (Section 6). Until the convergence criterion is met, the two stations involved in Algorithm 2, which protect a certain part of the area of coverage of the network, are treated as having detected an anomalous ionospheric gradient, but they still do not have any information of the size and propagation of the gradient. Therefore, a “Warning” is issued to warn the network that there is an anomalous ionospheric gradient in that area that cannot yet be estimated. Thus, the areas of the network that those stations cover should use the “worst-case” threat model.
Once the algorithm converges, it outputs the cross-correlation coefficient 𝛼𝑟(𝑡) and 𝜏𝑟(𝑡). While the convergence criteria is fulfilled, both parameters are very similar from one epoch to the next and can be considered constant.
For the times when the cross-correlation coefficients no longer fulfill the convergence criteria but 𝑑𝑡𝑐𝑟(𝑡) is one, the network indicates again the use of “worst-case” assumptions. Algorithm 2 continues searching for high cross-correlations until the gradient is no longer detected (𝑑𝑡𝑐𝑟(𝑡) changes from one to zero), and the buffer for station 𝑟 is reset to empty to wait for another gradient to come.
As stated in Section 4.2.1, when false alerts from different stations occur, they are treated as “real alerts,” and they also go through the cross-correlation procedure. However, since the source of these alerts is not real anomalous ionospheric gradients, Algorithm 2 does not find sufficient correlation between the test statistics from the different stations and the same satellite, and the estimation of the gradient parameters is not calculated. In these cases, Algorithm 2 triggers a “Warning” indicating that the GBAS stations affected by these false alerts should use more conservative approaches because the network cannot guarantee the integrity otherwise.
When at least three stations are impacted by the same gradient, we can already estimate the speed vector of the gradient, 𝐯. For this purpose, we start by applying the work accomplished in Mayer et al. (2009) and Hernández-Pajares et al. (2006). Under the assumption that the ionospheric disturbance propagates as a planar wave, the planar wave phase 𝜙 can be expressed as
11
where 𝐤 is the angular wave number vector, 𝐱 is the position vector, 𝜔 is the angular frequency, and 𝜙𝑖𝑛 is the initial phase. Note that all vectors are defined as row vectors unless it is specified differently and their transposed form is indicated with a 𝑇 superscript.
Given an ionospheric gradient affecting the observations of a certain GNSS satellite, the condition to obtain the same phase between the corresponding perturbation observed from two stations of the network (for example, the reference station and another station 𝑟) can be written as
12
where 𝐱𝐫 and 𝐱𝟎 are the position vectors of the IPPs at the times of detection of the gradient. These are expressed in a local reference coordinate system. That is, 𝐱𝐫 = [𝑥𝑟,𝐸𝑎𝑠𝑡,𝑥𝑟,𝑁𝑜𝑟𝑡ℎ]. Note that the variables that follow in the manuscript are also expressed in this coordinate system. Defining the vector slowness as 
as in Hernández-Pajares et al. (2006), we can express Equation (12) in terms of the vector difference of the position between IPPs 𝐱𝐫,𝟎, and the time delay 𝜏𝑟 as
13
Now, we consider Equation (13) for the different stations 𝑟 belonging to the network 
that have detected the gradient until the current time 𝑡 denoted as 𝑁. Thus, we define matrix 𝐗 and vector 𝐳 as
14
Here, 𝐗 contains the IPP position vectors and 𝐳 the time delays between the different stations 𝑟 and the station of reference. Note that the difference position vectors composing 𝐗 and the time delays in 𝐳 appear in the order in which the different stations detect the gradient. This order replaces the name of the stations represented by subscript 𝑟. The dimensions of 𝐗 and 𝐳 are [N-1 x 2] and [N-1 x 1], respectively, because the first station that detects the gradient is used as reference for calculating the cross-correlation coefficients and the time delays.
Then, isolating 𝐬 from Equation (13):
15
Since 𝐗−𝟏 in Equation (15) is not invertible, we use the Moore-Penrose pseudoinverse, and thus we can estimate 𝐬 as in Equation (16), which can be solved by least squares:
16
However, since we also have the information of the level of cross-correlation between the test statistics (𝛼𝑟), we use it to estimate 
. In this way, we introduce in the estimation of 
information about the trust that we can give to the measurements coming from the different stations. Thus, our estimation of 
can be solved with a weighted least squares as
17
where 𝐖 is a diagonal matrix containing the weights of the weighted least squares (based on the cross-correlation results) as
18
Then, the speed of the gradient 𝐯 is the Samelson inverse of the vector 
:
19
Nevertheless, as previously stated, we assume that the perturbation is local, which means that the result of comparing the measurements from two stations with a long distance between them (e.g., 200 km) might be inaccurate since the perturbation might change during its propagation from one station to the next. Moreover, the geometry of the IPPs in the sky should not be aligned in order to avoid singularities while calculating the speed vector.
Therefore, before the estimation of the gradient parameters, we group the IPPs corresponding to the different stations in clusters, and we validate these clusters to check if they are suitable to calculate a reliable estimation of the speed vector. At the time of detection in each of the stations, we compute the position vectors 𝐱𝐫 for each of the IPPs corresponding to the considered satellite and different stations 𝑟 impacted by the same gradient in a local coordinate system. Then, we form a cluster with these position vectors, which is validated attending to two different criteria: i) a radius around the central point of the local coordinate system to guarantee the locality principle, and ii) a geometry index that ensures that the geometry of the IPPs in the sky is acceptable. For the first criteria, we choose a validation radius of 200 km. This value is based on the maximum station separation distance considered in the literature that addresses the ionospheric threat model derivation in the CONUS region (Datta-Barua et al., 2010) and in Korea (Kim et al., 2015). For the second criteria, we define the geometry index (GI) as
20
This 𝐺𝐼 evaluates the impact of the geometry of the IPPs in the solution of Equation (17) and plays a similar role as the Geometry Dilution of Precision (GDOP) in the estimation of the position solution: a high GI means a bad geometry, whereas a low GI means a good geometry from the point of view of the resolution of Equation (17). The values that the GI can get are further discussed in Section 6.
Slope and width of the gradient
When the speed vector of the gradient is known, the slant slope of the gradient estimated by the considered satellite and station 𝑟 is determined by the following geometrical relationship:
21
where Δ𝑣𝑟(𝑡) is the relative speed between the gradient, 𝐯, and the IPP, 𝐯𝐈𝐏𝐏𝐫 (𝑡), projected in the direction of propagation of the gradient calculated as
22
As we can observe in Equation (21), the slope is calculated per epoch 𝑡, and it is not considered constant during the propagation of the gradient as the speed vector. The reason behind this is that we look for the “worst-case” or highest slope and therefore, we need to consider all the values of the 𝐼𝑡𝑒𝑠𝑡𝑟(𝑡) until we find the maximum. Once the maximum slope for this specific satellite-station pair is found, we assume that it is constant until the perturbation is not detected any more. Another possibility to calculate the slope of the gradient is the station-pair method, which calculates the difference of the slant ionospheric delay between two stations and divides by the distance between the stations (Datta-Barua et al., 2010). However, this methodology requires very accurate integer ambiguity resolution, and its estimation accuracy is often degraded by any remaining erroneous systematic offset. In contrast, the methodology used in this work, widely known as a time-step method, is not highly sensitive to remaining biases on ionospheric delay estimates, since the integer ambiguity is removed when calculating the first derivative of the ionospheric delay. Note that gradient slopes estimated by the time-step method are conservative because they include both spatial and temporal changes in the ionospheric delay when only the spatial changes are desired. However, this is acceptable when searching for the “worst-case” gradient estimations in real time.
The width of the gradient estimated by the considered satellite and station 𝑟 is calculated as
23
where 𝑇𝑊(𝑡) is the time in seconds that have passed since the gradient was detected until current time 𝑡. Note that the total width of the gradient can only be known when the whole gradient has passed and the station 𝑟 is not detecting the gradient any more.
All the satellites affected by the gradient separately estimate the gradient parameters and provide the current state of the ionosphere in the monitored area. Then, it is the task of each GBAS ground station to use this information to update the threat model to be used at a certain time epoch.
5 SIMULATION SETUP
In this section, we introduce the simulation setup that we use to evaluate the differences between known simulated gradient parameters and the parameters estimated by our algorithm. The simulation setup consists of the measurements recorded by the real network of stations introduced in Section 3 on the day 73 of year 2015 (one of the “quiet” days selected for the threshold derivation in all stations).
On top of the slant ionospheric delays calculated with these measurements, we simulate synthetic ionospheric perturbations designed to be representative of the GBAS threat model shown in Figure 1 occurring at different moments of the day. For that purpose, we define the synthetic perturbation as a planar wave front depicted in Figure 5. This figure shows an ionospheric gradient or change in the ionospheric delay values between the dark blue area and the dark red area, which are the areas where the ionospheric delay values are constant and no gradient is present. In this case, the ionospheric delay values are expressed in the vertical domain to be independent of the elevation of each of the satellites, which allows us to simulate the same ionospheric gradient for all the satellites. These vertical ionospheric delays are simulated by the utilization of two different simulation gradient parameters: the vertical slope in millimeters per kilometer and the width in kilometers. Therefore, when an IPP corresponding to a certain satellite and station moves into the region with the gradient, we compute the distance from the IPP to the point where the perturbation starts in kilometers and we multiply it with the vertical slope to get the vertical ionospheric delay. This vertical delay is then translated into a slant delay by multiplying with an obliquity factor (see the equation 5.28 in Misra & Enge (2006)) that depends on the elevation of the satellite. Furthermore, this ionospheric gradient moves with a constant speed (meters per second) and direction (degrees) over a “thin shell” layer at a height of 350 km above the Earth’s surface. Here, by “speed” and “direction” of the simulated gradient, we refer to the magnitude and direction of the speed vector 𝐯, measured in the clockwise direction from the North Pole. Notice that although the perturbation is simulated as a planetary planar wave, which is unrealistic, the measurement of the front is made by IPPs separated by only a few kilometers. Therefore, since an actual front size would be several times larger than the baselines of the network used, the approximation of a planar wave front is applicable to this work.
Example of one synthetic perturbation.
For our studies, we vary all the already defined simulation gradient parameters within their ranges in the CONUS and German threat models up to a maximum differential delay of 50 meters. In Table 3, we present the parameter bounds for our simulation. Note that the slope of the gradient is simulated in the vertical domain until the maximum value defined in slant for the CONUS threat model. The reason behind this is that a high elevation satellite presents a similar value for the slant delay as for the vertical delay, and therefore we need to consider the vertical values up to the maximum. This means that in the case of low elevation satellites much higher slopes than the maximum for the CONUS threat model are simulated as long as the differential delay achieved does not exceed 50 meters.
Simulation parameters
6 RESULTS AND DISCUSSION
6.1 Simulated gradients
6.1.1 Detection
In this section, we analyze the detection capabilities of the network of stations depicted in Figure 2. For this purpose, we initially introduce the performance of each of the monitoring stations inside the network by means of their detection thresholds (Figure 6). Then, we relate the minimum ionospheric rates that each of the stations can detect with the minimum ionospheric rates that the simulation produces when considering ionospheric gradients with different parameters (Figure 7). In this way, we are able to identify which of the simulated ionospheric gradients results in 100% detectability depending on the performance of each of the stations.
Detection thresholds for all stations considered derived with data from 10 nominal days of year 2015 and p𝑓𝑎 = 10−6.
Minimum ionospheric rate generated by a gradient with different simulated speeds and slant slopes that should be detected by the stations.
The detection thresholds depicted in Figure 6 show that all the stations of the network have similar performance with the exception of av01, which has a higher threshold overall. This fact is translated into worse performance of station av01 because the simulated ionospheric gradients have to generate higher ionospheric rates in order to be detected.
However, it might happen that the characteristics of an ionospheric gradient or the propagation of the satellite measuring it generates an ionospheric rate that is not high enough to trigger the detection thresholds. To illustrate this issue, in Figure 7, we show the minimum ionospheric rate in millimeters per second that is generated by a gradient with a certain slope and a certain speed considering all the GPS satellites, values of the rest of the parameters (width and direction), and time of occurrence of the gradient during the day. Note that the values shown in Figure 7 are calculated considering only the effects introduced by the synthetic gradients. That is, the noise, multipath, and nominal ionosphere present in the real nominal measurements on which these synthetic perturbations are simulated are not taken into account for this study. The reason behind this is that we seek to find the ideal performance that would be needed by the stations in order to have 100% detectability of all simulated gradients.
Therefore, each pixel of Figure 7 contains the minimum ionospheric rate generated by a gradient of a certain speed and a slant slope in a range between a lower limit plus 25 mm/km. As an example, the minimum ionospheric rate that was generated by a gradient of 900 m/s of speed and a slope between 675 and 700 mm/km was 260 mm/s (Figure 7). This value is low in comparison to the rate that we would expect for such a large gradient, but it corresponds to a satellite that was moving almost perpendicularly with respect to the direction of propagation of the gradient. This means that the propagation of the satellite also has an impact on how we observe the gradients, and all satellite propagations must be taken into account to avoid integrity issues.
Hence, the minimum rates showed in Figure 7 are the lowest that should be detected by all stations in the network in order to achieve 100% detectability for a gradient with these characteristics. Still, from both Figures 6 and 7, we can conclude that the detectability is not 100% for all simulated gradients. All stations present 100% detectability for all gradients that produce ionospheric rates above 45 mm/s (see Figure 6) independently of the elevation of the satellite or speed and propagation characteristics of the IPPs. These are all simulated gradients above the black dashed curve in Figure 7. For all the other simulated gradients that fall below this curve, achieving 100% detectability depends highly on the performance of each of the stations in the network.
Thus, there is a trade-off between the quality of the measurements that the monitoring stations should have and 100% detectability of harmful gradients for GBAS. For example, all stations with the exception of av01 are able to detect 100% of the gradients that produce ionospheric rates above 10 mm/s with satellites above 35 degrees of elevation. This corresponds to all gradients above the black continuous curve in Figure 7. This means that if stations ac59, av17, av16, and av20 do not detect any gradients with satellites above 35 degrees of elevation, the probability that there is an anomalous ionospheric gradient above 100 mm/km of slope and 100 m/s of speed in the area that they support is the same as the probability of missed detection of the network. Therefore, it would be possible to lower the maximum values of the threat model to be used in the supported GBAS stations in that area of coverage down to 100 mm/km of slope and 100 m/s of speed for satellites with elevations above 35 degrees.
For all other ionospheric rate values, the minimum detectable gradients could also be calculated, and these values could be used instead of the more conservative threat models as long as their values are lower than the conservative ones.
Notice that the calculation of the probability of missed detection of the network and the application of its value in the supported GBAS stations is out of the scope of this paper and will be addressed as part of future work.
6.1.2 Estimation
In the case of the estimation capabilities of the network, we analyze two different aspects: the estimation errors associated with each of the gradient parameters and the real-time performance. The results of the parameter estimation errors were calculated after the algorithm converged.
Gradient parameter estimation errors
First, we start with the evaluation of the gradient parameter estimation errors. These errors are calculated by substracting from the estimated value of the parameter the simulated value of it. In Table 4, we show the summarized mean and maximum estimation errors (maximum positive and maximum negative) for all simulations considered. Each row of Table 4 considers the change of one of the parameters among the values in Table 3 while maintaining all other parameters constant at values for which the simulation performed well. These values are 100 m/s for the speed, 180 degrees direction (propagation from North to South), 100 km of width, and 200 mm/km for the vertical slope. In this case, we see the impact of the change of one of the parameters on the estimation error of all the other parameters. As an example, the first row of Table 4 presents the mean and maximum estimation errors of all the gradient parameters when changing only the speed of the gradient from 0 m/s to 1200 m/s.
Summarized general simulation results
The results in Table 4 show that the errors in the columns “Mean” stay with a very low value for all the gradient parameters with the exception of the slope determination. Thus, the algorithm in general appears to work well under the simulation conditions examined. The particular case of the slope is further discussed in the subsection “Slope estimation error.”
If we consider the general results in the columns “Maximum” of Table 4, the following conclusions can be reached:
The greatest impact on the maximum estimation errors for all the parameters is due to the change of the speed (first row in Table 4), which also impacts at most the speed determination. This point is further discussed in the subsection “Speed estimation error.”
The maximum errors for the estimation of the slope are high for the change of each of the gradient parameters (sixth column in Table 4). This point is further discussed in the subsection “Slope estimation error.”
The maximum estimation errors for the direction and the width parameters are much lower than for the speed and the slope and present comparable values for all changes of the gradient parameters.
The change of the direction and the width parameters does not present any impact on the estimation errors. The change of the slope has an effect only when its value generates ionospheric rates that are close to the detection thresholds because the gradient might not be completely detected during all its duration.
Now, we study in detail the speed and the slope estimation errors.
Speed estimation error
As commented previously, the greatest impact on the speed estimation errors is the change of the speed itself. For a direction of the gradient that is not aligned with all or most of the stations on the ground, this mainly has two reasons: i) the alignment of the IPPs due to low elevation satellites (bad geometries) and ii) an insufficient time resolution (1s) that can cause appreciable errors in the detection and cross-correlation steps.
To study the impact of the IPP geometries on the estimation of the gradient parameters, we use the Geometry Index (GI) introduced in Equation (20). As explained in Section 4.2.2, this GI can give an idea of the suitability of the IPP geometry at the moment of calculating the gradient parameters. The higher GI gets, the worse the mentioned IPP geometry becomes. In Figure 8, we show the estimation error in speed versus the GI to evaluate the influence of the IPP geometry in the errors that we estimated. Here, the GI stays between 2𝑥10−4 and 4𝑥10−4 in most cases. These values stay within the limits for suitable IPP geometries, since the worst cases in the simulation showed values up to 9.3𝑥10−4 for very low elevation satellites. Moreover, the largest speed estimation errors correspond to the highest simulated speeds (green stars in Figure 8) while presenting average GIs. Therefore, the majority of these errors can be associated with the time resolution problem.
Error in speed estimation over Geometry Index for different simulated speeds of the gradient. Each marker represents a simulated speed and each point the GI calculated at the moment of the estimation.
Regarding the impact of the time resolution, we show in Figure 9 the mean (blue circles) and maximum speed estimation errors (up red triangles and down black triangles) when changing the simulated speeds of the gradient. As we can see, the mean estimation errors stay low, and we only notice a slight change for speeds of 0 m/s (the gradient is not moving) and 1200 m/s. In the first case, when the gradient is not moving, the errors that we observe are due to the geometrical approximations used in order to calculate the IPP locations at a specific shell height. This error is not very important in comparison with the error of the movement of the gradient and is hidden by it as the gradient moves. In the case of the gradient moving at very high speeds, the time resolution problem appears. In this case, the detection is not triggered exactly when it happens, and it could be delayed up to 1s, which in this case could mean up to 1,200 meters of error. Since some of the stations in the considered network (av01, av16, av17, and av20) are very close together (around 5 km), this can translate into large speed estimation errors. In the case of the maximum errors, this behaviour is more evident, and the more we increase the speed of the gradient, the worse the estimation gets.
Mean (blue circles), maximum positive (red up triangles), and maximum negative (black down triangles) values for speed estimation errors with respect to the simulated speed of the gradient.
Slope estimation error
In the case of the slope or spatial gradient magnitude, one factor that influences both the mean and the maximum errors is that we consider as true slope only the one generated by the simulator (i.e., the synthetic gradient). This means that the nominal ionosphere on top of which these synthetic gradients are simulated is not considered as part of the ground truth since it was not possible to separate it from the multipath and noise coming from the real measurements. However, the “nominal” ionospheric gradients present in the measurements represent a small part of the ionospheric rates used as test statistics.
The largest impact on the slope estimation is the one caused by the phase noise and multipath when calculating the ionospheric rates every second (upper part of Equation (21)). This effect is depicted in Figure 10 with respect to the change of the speed. As we can see, the maximum values become large due to the fact that the test statistics are noisy and the fact that we divide their values by small relative speeds between the gradient and the IPPs. However, the more we increase the speed of the gradient, the more this error is smoothed by dividing by larger relative speeds. The problem in this case is that we could underestimate the slope of the gradient (for speeds above 450 m/s) since we get up to a -25 mm/km error. This issue should be taken into account by calculating the “maximum possible slope” from the “estimated slope” when determining the integrity of the methodology, which is part of future work.
Mean (blue circles), maximum positive (red up triangles), and maximum negative (black down triangles) error in slope estimation for different simulated speeds of the gradient.
Real-time performance
Finally, we analyze the real-time capability of our algorithm. In Figure 11, we show an example of the 
calculated for all stations considered and satellite G03. Here, the test statistic values are noisy, and when the algorithm applies the cross-correlation algorithm in real time considering station ac59 as reference, at the beginning, it only finds noise. Therefore, the algorithm in real time needs a certain time to converge. In Figure 12, we present the cross-correlation coefficients for the example depicted in Figure 11. As we can observe, the algorithm needs at least 100 seconds to satisfy the convergence condition presented in Equation (10).
Example of simulation test statistics and thresholds for satellite G03 and day 73 of year 2015.
Cross-correlation coefficients calculated in real-time for the example of Figure 11.
Note that these maximum estimation errors and necessary waiting times for the convergence of the algorithm would be considered by the network when calculating the “worst-case” real-time gradients that would be transmitted to the GBAS stations.
6.2 Real gradients
Additionally, we also evaluate our algorithm with a real gradient measured by the network in Alaska (Figure 2). The selected real gradient was measured by satellite G03 and the stations in the network under study on day 76 of year 2015, i.e., during the St. Patrick’s Day Storm.
6.2.1 Detection
After applying the thresholds derived as explained in Section 4.1 and depicted in Figure 6, we show the detected gradients in Figure 13 represented with the respective markers for each station. The detected gradients coincide with the steepest slopes in slant ionospheric delays. Thus, we conclude that the gradients are adequately detected. Note that, if the monitoring stations have poor performance and their thresholds are high, some anomalous gradients could remain undetected. However, this does not imply an integrity issue for the GBAS stations supported by the network. If a monitoring station does not detect any gradients, it does not automatically assume nominal conditions, but instead it considers that there might be a gradient present that is equal to its minimum detectable gradient (introduced in Section 6.1). In this way, the network ensures that the supported GBAS stations are always considering the largest undetectable gradient that could be affecting them.
Slant ionospheric delays for the studied gradient measured by satellite G03 (day 76 of year 2015, Alaska). The lines with the respective markers for each station show the times when the gradient was detected.
In Figure 14, we also show the test statistics, 
, for all the stations in the times where they exceeded the corresponding thresholds.
Slant ionospheric rates of the studied gradient measured by satellite G03 (day 76 of year 2015, Alaska). The values of the test statistics are highlighted when the thresholds were exceeded.
6.2.2 Estimation
Figure 13 shows that the gradient presents a similar shape when being measured by all affected stations. However, in Figure 14, the test statistics suggest that the gradient is narrower and steeper when it reaches the last stations impacted by it, which implies that our previous assumption of a non-changing ionospheric gradient does not hold for the farthest station. Note that this assumption is used to simplify the analysis of the problem, but it is expected to not be true in general for real gradients measured by stations separated by more than a few kilometers.
Nevertheless, the algorithm is still able to estimate the gradient parameters. First, it assumes as the station of reference the first station impacted by the gradient (ac59). However, as we can observe in Figure 14, the shape of 
(black stars) is different from all the other test statistics. Thus, the cross-correlation coefficients calculated with respect to this station result in low values. This issue can be seen in Table 5, where we show the cross-correlation coefficients calculated in post-processing considering as reference each of the different stations in the network. The stations in Table 5 appear in the order they detected the real ionospheric gradient.
Cross-correlation coefficients calculated in post-processing considering the different stations as reference
Thus, it is necessary to adapt the requirements for a better selection of the reference station. Therefore, the algorithm searches for the first impacted station whose test statistic presents above 0.9 cross-correlation with the test statistics from at least two other stations of the network. A cross-correlation coefficient of one between the tests statistics computed at two different stations would mean that the perturbation is the same, but delayed by a certain time interval, as stated in Section 4.2.2. However, since the real anomalous gradients are not ideal, we have lowered that value to 0.9 to allow enough margin for the algorithm to find the station considered as “reference.” Note that, this value was tested only with the real gradient under study depicted in Figure 13 since in the simulations it was not necessary because all cross-correlation coefficients were very close to one due to the simulated ideal gradients.
Thus, according to Table 5, station av16 is considered as the reference from now on. However, the algorithm still keeps the information from all past stations, since it applies weighting according to the cross-correlation coefficients (see Equation (17)). In this process, the information coming from stations with a cross-correlation coefficient below 0.5, i.e., 𝛼𝑚𝑖𝑛 in Section 4.2.2, is excluded. As explained in Section 4.2.2, below this cross-correlation value, the ionospheric gradient should not be considered the same, and the inclusion of this information could make the estimation of the gradient parameters worse.
Thus, assuming station av16 as reference, we initially calculate the gradient parameters in post-processing to get the best estimation possible as ground truth. Then, we compare the post-processing results with the ones obtained in real time to evaluate the feasibility of the real-time concept.
Estimation of the gradient parameters in post-processing
Since the gradient under study propagates with a high speed, as we can already observe in Figure 14, where three of the stations detect it almost at the same time, we use a spline interpolation of the data at 10 Hz to achieve better accuracy in our ground truth and avoid the problem of time resolution.
The results in post-processing are summarized in Table 6. The ionospheric gradient measured is a short-duration perturbation that travels with an estimated direction of 203.3◦ (from the Northeast to Southwest), a high speed, 2,473.3 m/s, and not a very steep slope, up to 55.8 mm/km. These results are compatible with a propagation of the perturbation following the magnetic field lines that in this region have a declination of 15◦ and also agree with the short time delays calculated between the stations. Moreover, the slope and width estimations indicate that the high values of the test statistics in Figure 14 can be attributed more to the high speed of the gradient than to its size. These kinds of ionospheric perturbations, local and generated in the auroral region, are often the source of Large Scale Travelling Ionospheric Disturbances (LSTIDs) at mid-latitudes (Borries et al., 2009; Wang et al., 2007).
Estimated real gradient parameters in post-processing
Estimation of the gradient parameters in real time
Applying the algorithm in real time with the data recorded at 1 Hz presents the same limitations as in the simulation: we need a convergence time for the cross-correlation and sufficient time resolution. The first limitation can be seen in Figure 15, where the convergence requirement (Equation (10)) is not met until 09:20:11 (local time). Moreover, once the algorithm converges, it outputs time delays of 0 seconds for two of the stations (av01 and av20 in Figure 16) because the time resolution is insufficient. These limitations can also be recognized in the estimation of the speed in real time (Figure 17), where we get high errors when the algorithm did not converge and no output when the time resolution is insufficient. Note that the values of the time delays before convergence would not be used by the algorithm. Therefore, the speed estimations before convergence of the algorithm are only shown in Figure 17 for better understanding.
Cross-correlation coefficients calculated in real time for the test statistics depicted in Figure 14 considering as station of reference av16.
Time delays calculated in real time for the test statistics depicted in Figure 14 considering as station of reference av16.
Estimated speed in real time of the real gradient depicted in Figure 13.
After 21 seconds of continuous tracking by all stations in the network, the algorithm converges, and it has enough resolution to output results comparable to the ground truth calculated in post-processing, 2,083.9 m/s for the speed estimation and 212.9◦ for the direction.
Finally, the slope and width parameters are also estimated. In the case of slope determination, two possibilities are considered. The first possibility is calculating it in real time. However, the maximum slope that we find is 34.6 mm/km, because the “worst-case” gradient occurs before the algorithm converges. The second possibility is recomputing the slope backwards to find the “worst-case” that could be transmitted to the GBAS stations. This second solution is depicted in Figure 18. Here, the slope estimations are higher than in Table 6 because the speed estimation used (2,083.9 m/s) is lower. Thus, the same ionospheric rate (Figure 14) is translated into higher slopes.
Estimated slope of the gradient depicted in Figure 13 using speed estimation after convergence.
The width parameter is determined after the gradient has finished affecting all stations. The results for the estimation of this parameter per station are: 50.0 km for av17, 60.4 km for av16, 58.3 km for av01, and 58.4 km for av20. In this case, the width is shorter than in the post-processed results because the estimation of the speed is lower. Thus, the same time duration of the gradient 𝑇𝑊 in Equation (23) is attributed to a shorter width.
7 APPLICABILITY OF THE DUAL-FREQUENCY IONOSPHERIC MONITORING CONCEPT TO GBAS
In this section, we discuss the applicability of the proposed concept to the current GAST-C and GAST-D systems.
A GAST-C ground station would use a variable maximum gradient slope, 𝑔, in the calculation of anomalous ionosphere-induced differential range errors in the geometry screening algorithm instead of using the constant value from the threat model as done today. This variable maximum gradient slope would depend on the information received from the network, which could be:
i. The minimum detectable gradient for the coverage area of the network when the network does not detect any anomalous ionospheric activity. This value would be used as the maximum anomalous gradient that can currently occur instead of the maximum “worst-case” gradient defined by the regional ionospheric threat model.
ii. The “maximum possible slope” estimated and overbounded in real time for the coverage area of the network when the network detects an anomalous ionospheric gradient and the algorithms converge. This value would be larger than the minimum detectable gradient from (i) but typically smaller than the maximum “worst-case” gradient from the ionospheric threat model.
iii. A “Warning” that implies the use of the “worst-case” ionospheric threat model applicable to GBAS for that area when the network detects an anomalous ionospheric gradient but the algorithms are not yet able to reliably estimate the gradient parameters.
Therefore, the integrity parameters transmitted by the GAST-C station would require less inflation most of the time. Thus, availability would increase in regions where conservative ionospheric threat models have to be applied to protect integrity.
A GAST D ground station would receive the probability of missed detection from the network. This probability of missed detection could be used as a prior probability of occurrence of an anomalous ionospheric gradient in certain monitors (such as the Ionospheric Gradient Monitor, IGM), which would significantly relax the detection thresholds and thus improve the false alert rate from these monitors. This concept differs from a similar strategy applied in Yoon et al. (2020) by using real-time or near real-time observations instead of a statistical approach from historical data. Moreover, it could help relax the stringent siting criteria of the IGM and enable the use of GAST-D in areas with more challenging ionospheric conditions.
A methodology for deriving both the “worst-case possible estimated gradient slope in real time” and the “probability of missed detection of the network” is part of future work.
8 CONCLUSIONS AND FUTURE WORK
In this work, we present a method capable of detecting ionospheric gradients in real time and estimating their parameters in near real time based on a wide area network of dual-frequency and multi-constellation GNSS monitoring stations.
First, we explain the derivation of detection thresholds for each of the monitoring stations. Then, we describe the methodology during real-time operation of the algorithm. Here, we introduce weighting within the equations for the estimation of the gradient parameters that takes into consideration the cross-correlation coefficients between the test statistics coming from different pairs of stations. Additionally, we define a geometry index to check if the clusters formed by the IPPs to estimate the gradient parameters are suitable to limit geometrical approximation errors in the methodology.
Finally, we evaluate our algorithm with simulated gradients and a real gradient measured by a reference network in Alaska. Simulation results show promising performance of the algorithm since the mean estimation errors are low. Therefore, our algorithm shows potential to support GBAS by detecting gradients and estimating the gradient parameters instead of always assuming “worst-case” conditions. In our study, the largest estimation errors occur when varying the speed of the gradient. For high simulated speeds, the highest influence on the estimation errors is the time resolution of the measurements, which causes the gradient to not be detected exactly at the moment it affects a station. Additionally, the real-time capabilities of the algorithm are analyzed. Here, the algorithm needs time to converge that depends on the characteristics of the ionospheric gradient and the performance of the monitoring network.
Results evaluating the performance with a real gradient show the need to adapt the algorithm to the characteristics of the ionospheric perturbations in the studied area, especially if they are fast traveling, have a short duration, and are fast changing with their propagation. For the gradient studied, this is translated into adding minimum limits for the cross-correlation coefficients used as weights in the resolution of the speed vector estimation. We require that the cross-correlation coefficients of a certain station are at least 0.9 with two other stations to select it as reference and 0.5 to consider that a station is measuring the same perturbation as the reference. In this way, we could obtain promising estimations of the gradient parameters.
Future work will continue improving the methodology further to reduce the impact of carrier phase noise and multipath on the detection and estimation of the gradient parameters and to provide robustness against more complex simulated gradients. Additionally, the methodology for translating the network’s outputs into GBAS station actions and for evaluating the integrity of the presented algorithm will be developed.
HOW TO CITE THIS ARTICLE
Caamano M, Juan JM, Felux M, Gerbeth D, González-Casado G, Sanz J. Network-based ionospheric gradient monitoring to support GBAS. NAVIGATION. 2021;68:135–156. https://doi.org/10.1002/navi.411
ACKNOWLEDGMENTS
Open access funding enabled and organized by Projekt DEAL.
- Received April 3, 2020.
- Revision received December 1, 2020.
- Accepted December 13, 2020.
- © 2021 The Authors. NAVIGATION published by Wiley Periodicals LLC on behalf of 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.
Jump to section
Related Articles
Cited By...
- No citing articles found.