© A.R.F. Martinon et al., Published by EDP Sciences 2023

1 Introduction

During daylight, the equatorial ionospheric plasma is lifted due to eastward electric fields. After sunset and up to about 21 LT (local time) this upward plasma movement is intensified giving origin to a large vertical electron density gradient at the ionospheric F region base, which is one condition for plasma irregularities to be generated through the Rayleigh Taylor mechanism (Abdu, 2019). If the conditions are favorable these irregularities grow up and while moving upwards at the equatorial region they give origin to large structures known as ionospheric bubbles that map along the magnetic field lines reaching low latitudes. The bubbles can reach the Equatorial Ionization Anomaly (EIA) crests which are regions of high electron density around 20° north and south of the magnetic latitude. The EIA southern crest is located in the Brazilian region. Inside the bubbles smaller scale irregularities are generated by cascading processes and those having the scale size of the first Fresnel zone (about 400 m) cause amplitude and phase scintillation in the L band frequencies like those of the Global Navigation Satellite System (GNSS) signals. The ionospheric scintillation can cause loss of lock of the GNSS signal that crosses the bubbles, with a consequent decrease in the number of available satellites and a degradation in the GNSS receiver performance (Bandyopadhayay et al., 1997; Kintner et al., 2001) and the required accuracy to meet aerial navigation could not be reached. During scintillations the receiver signal can be affected by a fast phase variation that could cause a loss in the phase lock (Leick, 1995; Skone et al., 2001) and/or a rapid amplitude fluctuation, when the signal becomes too low to be sensed by the loop discriminator (Humphreys et al., 2004; Kintner et al., 2005). Augmentation systems like the Satellite Based Augmentation System (SBAS) (de Paula et al., 2007) and the Ground Based Augmentation System (GBAS) implementations are drastically affected during ionospheric scintillation and even making their installation unfeasible.

The making of scintillation maps requires the interpolation of Ionospheric Pierce Point (IPP) samples, given by the S4 values for each IPP of each satellite-station link considering the set of GNSS stations of the given area and time interval. In order to investigate ionospheric plasma irregularities and/or bubbles over the Brazilian territory, de Rezende et al. (2007) proposed to generate regional scintillation maps, i.e. maps of the S4 amplitude scintillation index. Vani (2018) employed a similar approach, with modified pre-processing options to generate these maps. Other works propose the generation of scintillation maps for different parts of the globe mainly by approaches based on kriging interpolation (Hamel et al., 2014; Kieft et al., 2014; Amabayo et al., 2021). All the cited approaches result in a temporal sequence of maps using or not a priori IPP samples from a previous map of the sequence. A more recent work (Koulouri, 2022) proposed the use of Kalman filtering for generating scintillation maps in real-time with a 1-minute resolution.

The current work analyzes four approaches for the generation of scintillation maps covering Brazil, two existing and two new approaches proposed here. Each one of these approaches involves an interpolation method and a set of pre-processing options. All approaches are named after the corresponding interpolation method. The two existing approaches are (1) GDA (Matlab’s GRIDDATA), employed at INPE by the space weather monitoring program (EMBRACE, available at http://www2.inpe.br/climaespacial/portal/sci-home), and (2) IDW (Inverse Distance Weighting), as proposed by de Rezende et al. (2007) and Vani (2018). The two new approaches are (3) based on RBF (Radial Basis Function) interpolation, described in Broomhead & Lowe (1988), and (4) GPR (Gaussian Process Regression), described in Rasmussen & Williams (2016). The analysis of these four approaches includes a baseline set of pre-processing options, and also an improved set of pre-processing options proposed here.

The four approaches were evaluated for intervals of time that presented strong scintillation events. The quality of the maps was assessed using a set of metrics and evaluating the smoothness and absence of artifacts of the maps. It can be concluded that the GPR interpolation, with a specific set of pre-processing options, allows the generation of the most accurate maps for the considered data. Each map integrates IPP samples of a 16-minute interval (since there is 1-minute overlap, the actual resolution is 15 min) that is suitable for the monitoring of scintillation and for studying the evolution of plasma irregularities and/or bubbles. This resolution was chosen to cope with the high temporal and spatial variability of the scintillation, but not long enough to hide its temporal trend. However, a comparison was also performed with a fifth approach, the one proposed by Koulouri (2022), using a 1-minute resolution, showing also the corresponding animations of the sequence of scintillation maps.

The proposed approach for generating more accurate scintillation maps has a low computational cost, allowing the making of real-time scintillation maps covering the Brazilian territory, even with a 1-minute resolution in a standard personal computer. Such capability is already being implemented at the National Institute for Space Research (INPE), Brazil.

2 Proposed approach for ionospheric scintillation maps

Amplitude scintillation is given by the S4 index, estimated from radio frequency signals acquired at 50 Hz by receivers of GNSS monitoring stations. Typically, the signal intensity of 3000 samples of a 1-minute interval is used to calculate the S4 index value given by the intensity standard deviation normalized by the mean intensity (van Dierendonck et al., 1993). The S4 index is calculated for each GNSS satellite-station link from one or more of the several constellations currently available (GPS, GLONASS, GALILEO, BEIDOU, SBAS, etc.). As a consequence, the number of GNSS satellites that are visible for a station can exceed 30 in the same minute. In addition, for a given constellation, calculations of the index S4 are compatible between GNSS stations equipped with different receivers (de Paula et al., 2021). The present work employs all available GNSS constellations for the considered time interval to produce the scintillation maps, except for SBAS satellites which are only used to assess errors of the maps.

Independently from the adopted approach, the generation of the scintillation maps employs a large number of IPP samples (in this case, more than 3000), thus requiring to aggregate groups of IPP samples into “interpolation samples” in order to make the interpolation feasible. In this work, for any map generation approach, it is adopted an interpolation grid with the same spatial resolution of the map (0.25° × 0.25°), covering latitudes −39° to 9° and longitudes −78° to −30°, and composing a regular grid of 193 × 193 elements. Besides the interpolation grid, it is also adopted a coarser “aggregation grid” (1.0° × 1.0°) that defines square cells of the same size (1.0° × 1.0°) centered at each grid point. All IPP samples inside a cell are then assigned to the corresponding grid point with a reduced S4 value, composing an “interpolation sample”. Some tests were performed showing that a finer aggregation grid did not improve the quality of the maps, but increased processing time.

The baseline set of pre-processing options was originally employed in the IDW approach, requiring the following assumptions/steps:

1. Model the ionosphere as a thin spherical shell over the Earth at the mean altitude of 350 km, which corresponds to the maximum observed TEC (Total Electron Content);

2. Define the IPP for each satellite-receiver link, as shown in Figure 1, as being the intersection of the line-of-sight receiver-satellite (the slant path) with the ionosphere, using the receiver coordinates and azimuth and elevation angles of the satellite (Prol et al., 2017);

   Fig. 1 Geometry of a satellite-ground station line-of-sight depicting the Ionospheric Pierce Point.

3. Define an aggregation grid with 1.0° × 1.0° resolution with square cells of the same size centered at each grid point, and reduce the IPP samples inside a cell by the S4 maximum or average value to an “interpolation sample” (cells without IPP samples are discarded);

4. Perform the interpolation using the “interpolation samples”, each one with a reduced value at the aggregation grid point, obtaining the interpolated values for each point of the considered interpolation grid, i.e. obtaining the scintillation map (de Rezende et al., 2007).

The baseline set of pre-processing options as well as the newly proposed ones were applied to the four map-generation approaches: (1) GDA, (2) IDW, (3) RBF, and (4) GPR. The GDA is a triangulation-based cubic interpolation, implemented here using the Python equivalent of Matlab’s GRIDDATA. The IDW interpolation is a modified Shepard’s method, which includes an influence radius R that determines the neighboring interpolation IPP samples that will be considered in the weight definition, and weights each “interpolation sample” by the inverse of the square of the distance between its position and the interpolation grid point. The RBF method allows the use of a set of kernel functions that determines the way the interpolation will be carried out, in this case, the thin-plate kernel was chosen. Finally, the GPR method also employs a set of kernel functions to modify the interpolation result, in this case, a Rational Quadratic kernel. Except for GDA, a distance metric must be defined, usually the Euclidean distance, but the Haversine distance was used instead. The Haversine distance determines the great-circle distance between two points on a sphere given their longitudes and latitudes.

In order to generate more accurate scintillation maps, some modifications were made to the baseline pre-processing options. The first modification was not to consider the S4 measured along the receiver-satellite slant path, but its vertical projection. Such projection takes into account the geometrical effects on the measurements made at different elevation angles (Spogli et al., 2009). The vertical projection is achieved using the spectral slopes of the phase PSD (Power Spectral Density) for each observation, which are provided by the INCT GNSS-NavAer network (Vani et al., 2017) stations at a 1-minute rate, while for the stations of the Low-Latitude Ionospheric Sensor Network (LISN) (Valladares & Chau, 2012), the spectral slope value is assumed as 2.6 (Rino, 1979, Spogli et al., 2013). The same assumption is made for some 1-minute INCT observations that do not have the spectral slope value. The second modification was the introduction of a new reduction function, defined by the average of S4 values above the 3rd quartile (top 25% values), being the reduced value assigned to the corresponding grid point, composing an “interpolation sample”. A third modification was the location of the “interpolation samples”, given by a centroid instead of the regular grid point. Considering the set of IPP samples that are related to an “interpolation sample”, each centroid is defined according to the reduction function, as being: (i) the location of the IPP sample with the maximum S4 value, (ii) the average location of all IPP samples or (iii) the average location of the IPP samples with S4 values above the 3rd quartile.

The set of pre-processing options is based on a regular aggregation grid with the corresponding cells centered at each grid point. The IPP samples of each cell are grouped/aggregated with a reduced S4 value given by a reduction function, composing an “interpolation sample”. In order to simplify the presentation of the results, three-letter abbreviations were adopted to denote the pre-processing options that define the “interpolation sample” (each one has a reduced value and a location), as below:

• [S/V] – S (slant S4 value) or V (S4 value projected vertically) of each IPP sample;

• [M/A/Q] – Reduction function for the S4 values of the IPP samples of a cell, given by M (maximum), A (average) or Q (average of S4 values above the 3rd quartile);

• [R/I] – Location of the “interpolation sample” for each cell, being R (the regular grid point of the cell) or I (a centroid defined according to the choice of the reduction function).

The baseline pre-processing options of the former approaches are SMR or SAR. The newly-proposed options for all approaches are VMR, VMI, SMI, VAR, VAI, SAI, VQR, VQI, SQR, and SQI. Henceforth, each of the 4 approaches will be denoted with the corresponding set of preprocessing options, for instance, IDW(SMR) to refer to the IDW-interpolation approach with the pre-processing options SMR.

Therefore, each map generation approach was tested with each possible set of pre-processing options, being evaluated by a set of error and correlation metrics in relation to the S4 values of the test set of IPP samples. Considering any particular map, the set of IPP samples employed to generate the map (interpolation IPP samples) and the set of IPP samples employed to test the accuracy of the map (test IPP samples), are mutually exclusive, being defined according to a validation/partition scheme, as shown in Section 4.

3 Data description and filtering

The GNSS data used in this work was obtained from two ionospheric monitoring networks, the LISN monitoring stations covering most of South America, and the INCT stations, which cover only the Brazilian territory. LISN receivers are only compatible with GPS and SBAS satellites, while INCT receivers are multi-constellation. These networks provided a total of 45 ionospheric monitoring stations, with the corresponding locations shown in Figure 2. Not all the 45 stations of these networks have been continuously operational, with the total number of available stations ranging between 26 and 34.

Fig. 2 Geographic location of the 45 GNSS stations of the LISN and INCT networks that provided scintillation data for this work.

Table 1 presents the geographic locations and the magnetic dip latitudes for the 32 LISN monitoring stations, while Table 2 presents the same information for the 13 INCT monitoring stations. Since the intervals of time selected for this work are between November 2013 and November 2014, the magnetic dip latitude curves in Figure 1 and the values presented in these tables were derived using the year 2014 as reference and were calculated using the IGRF 13 model (Alken et al., 2021).

Since the format of the data provided by the two monitoring networks is different, a unified data format was adopted, as shown in Table 3. The “i_lat” and “i_lon” values correspond to the IPP coordinates. The “s4v” value corresponds to the vertical projection of the slant S4 value, as proposed by Spogli et al. (2009) and the spectral slope of the phase PSD (p) was defined according to Rino (1979) for LISN receivers (p = 2.6), and Spogli et al. (2013), for INCT receivers (p being measured every minute).

Scintillation data was acquired for 10 selected non-consecutive one-hour intervals that mostly presented moderate and strong scintillation events. Each one-hour interval was taken from a different night between November 2013 and November 2014. Each one of the 10 one-hour intervals of time was divided into four consecutive intervals, each one integrating 16 min of IPP data to generate the corresponding sequence of maps using a sliding window of 15 min. Therefore, temporal resolution is 15 min and there is a 1-minute overlap between consecutive maps. Figure 3 presents the distribution of the S4 values classified according to the scintillation severity for each one-hour interval using the following thresholds: null [0.00–0.15], weak (0.15–0.30], moderate (0.30–0.70], or strong (0.70–1.4]. S4 values above 1.4 were considered outliers, being rounded to that value. In order to filter out data of L-band links affected by ground interference and the related multipath reflections, only satellite elevations higher than 30° were considered. In addition, for each interval, the corresponding values of the planetary Kp geomagnetic index and the adjusted F10.7 index (solar flux at the 10.7 cm wavelength) are presented, showing the absence of geomagnetic storms.

Fig. 3 Bar charts showing the distribution of the IPP samples according to the scintillation class for the 10 one-hour non-consecutive selected intervals of time. Each bar chart depicts the date and starting time of a one-hour interval that is further divided into four 16-minute intervals for generating the maps (with a 1-minute overlap between them). The corresponding values of the geomagnetic index Kp and the adjusted F10.7 solar flux index are also included for each one-hour interval.

4 Error and correlation results

The presented 12 pre-processing options and 4 interpolation methods were evaluated for the set of 40 sixteen-minute intervals extracted from the 10 non-consecutive one-hour datasets described in the previous section. This is the main test set of 40 maps. Furthermore, an additional test set of maps is described in Section 4.1, which is specific for a comparison with a recent work (Koulouri, 2022).

Concerning the main test set of 40 maps, for compatibility with the receiver models, each one-hour dataset is expanded by an extra minute, being split into 4 sixteen-minute intervals with 1-minute overlaps, as previously mentioned, thus resulting in a resolution of 15 min. Scintillation maps were generated using a 0.25° × 0.25° latitude and longitude resolution, covering latitudes −39° to 9° and longitudes −78° to −30°, which corresponds to a regular grid of 193 × 193 elements. The evaluation of the maps employs each original IPP sample comparing their S4 value with the one obtained from the interpolated map. This is achieved by performing a distance-based weighted average of the four grid points of the map that are closer to the IPP sample position. It is important to distinguish between the interpolation points (obtained from the IPP samples) used to interpolate the map, and the test points (IPP samples used to evaluate the map). These two mutually-exclusive sets of IPP samples are obtained from the original set of IPP samples described in the previous section (10 one-hour sets), which is split according to a validation scheme.

Two cross-validation schemes were chosen for the evaluation tests: Stratified Shuffle Split (SSS) and Leave One Group Out (LOGO). The SSS scheme splits IPP samples into interpolation/test sets using stratified randomized folds, with each fold preserving the percentage of samples for each scintillation class (null, weak, moderate, and strong). A usual 10-fold scheme was chosen, thus the IPP samples in a given 16-minute interval were split into 10 sets, using 9 sets for interpolation and the remaining set, for evaluation. This scheme is repeated 10 times, swapping the evaluation set. The LOGO scheme also splits the IPP samples into interpolation/test sets, but each test set refers to IPP samples of a given GNSS station left out of the interpolation. Therefore, IPP samples of all but one station are used to generate the map for the 16-minute interval, and the map is evaluated using the IPP samples of the remaining station. This scheme is repeated swapping the station to be tested. In particular, in one case, a set of three stations in the same city were left out.

Considering the 40 sixteen-minute intervals, the 10-fold SSS scheme results in 400 scintillation maps, for each of the 4 approaches and each of the 12 pre-processing options, yielding a total of 19,200 maps (40 × 10 × 4 × 12). Each iteration of the 10-fold scheme (for each approach and each set of pre-processing options) employed 100,470 IPP samples for interpolation, 11,138 for evaluation/test, both from the GPS, GLONASS, and GALILEO constellations and further 15,413 IPP samples for test from the SBAS satellites. Therefore, the IPP samples of the test set in each iteration amount to about 21% of the total IPP samples.

In the case of the LOGO scheme, the 4 approaches were only tested with the best-considered set of pre-processing options (VQI) and 5 LOGO options, according to the left-out station, resulting in 800 maps (40 × 5 × 4 × 1). Each of the 5 LOGO options employed an average of 105,740 IPP samples for interpolation and 5,445 IPP samples for evaluation/test, both from the GPS, GLONASS, and GALILEO constellations and (in average) further 694 IPP samples for test from the SBAS satellites. The IPP samples of each test set amount to 5.5% of the total IPP samples.

As already mentioned, the IPP samples for interpolation and test sets were obtained from the GPS, GLONASS, and GALILEO constellations. The IPP samples of the SBAS constellation were only used for the test sets, and since the SBAS satellites are geostationary, for a 16-minute interval, there are 16 samples of a fixed-position IPP. A single S4 value is then obtained for each SBAS station using the same reduction function employed in the generation of the map.

In order to compare the 4 interpolation methods with the 12 sets of pre-processing options, 7 different error and correlation metrics were employed in the evaluation of the maps using the SSS scheme: overall Mean Absolute Error (MAE), overall Root Mean Square Error (RMSE), Mean of Maximum Absolute Errors (MXAE), Mean of Minimum Errors (MMIN), Mean of Maximum Errors (MMAX), overall Standard Deviation of Absolute Error (STDA) and the overall Pearson’s Correlation Coefficient (CORR). Therefore, MXAE, MMIN, and MMAX are calculated from the corresponding values of each map.

The complete error evaluation table for all the 12 pre-processing options employed for each interpolation method is available in the Zenodo repository, as described in the “Data availability statement” section, in the file “Complete table of interpolation errors and correlation.csv”. This table allows to conclude that the proposed GPR approach performs better than the remaining three approaches. Therefore, for the sake of simplicity, Table 4 shows only the results for the GPR interpolation method and the 12 pre-processing options. Considering the same table, it is difficult to choose the best set of options. Clearly, options with maximum value reduction (SMI, SMR, VMI, VMR) and also the slant S4 3rd-quartile reductions (SQR and SQI) consistently presented higher errors and lower correlations than the remaining VQR/VQI, SAR/SAI, and VAR/VAI, with the corresponding vertically-projected options performing better. In addition, options SAR/SAI and VAR/VAI performed better than VQR/VQI. However, these seven error and correlation metrics are not the only criteria, and VQI was chosen as the best overall set of pre-processing options based on the smoothness and absence of artifacts of the generated map. The latter feature shows how well the interpolation approach extrapolates grid points in regions with low density of GNSS stations.

Concerning the reduction function of the S4 values of the IPP samples, the Q pre-processing option (average of the 3rd quartile) resulted in errors that are bound by those of option M (maximum value) and option A (average value). Despite the lower MAE and RMSE errors of the option A, the corresponding MMIN errors were worse, thus underestimating strong scintillation. In addition, option M presented the worst MMAX errors, thus overestimating strong scintillation. Furthermore, for any approach and set of pre-processing options, it can be noticed that the V pre-processing option (vertical projection of S4 value) results in scintillation maps with lower errors than the corresponding S options (slant S4 value). Finally, changing the location of the “interpolation samples”, given by a centroid instead of the regular grid point lowered the error metrics for the IDW, RBF, and GPR approaches, but not for the GDA approach, which present lower errors with the “interpolation samples” located at regular grid points.

The LOGO cross-validation results are shown in Table 5 shows the results for the 4 interpolation approaches and only for the VQI set of pre-processing options. It is possible to check the performance for regions of high density (SJCE/LSJK/SJCU and INCO) and low density (LBSB, LDOU, and LNTA) of GNSS monitoring stations. For high GNSS-station density areas, the GPR and IDW results are similar, and better than those of GDA and RBF. For low-density areas, the GPR approach performed better, proving its robustness.

4.1 Error and correlation results for the additional test set

Furthermore, the GPR(VQI) approach was employed to generate maps using the same dataset evaluated by Koulouri (2022), which covers the period from December 1, 2014 21:30 UT to December 2, 2014 05:24 UT with a temporal sliding window of 1-minute resolution, resulting in 475 scintillation maps, with the same interpolation grid (193 × 193 grid points) and spatial resolution (0.25° × 0.25°). This evaluation used the same set of seven different error and correlation metrics for the cross-validation schemes already shown. Since Koulouri (2022) stated the use of prior S4 information to generate the scintillation maps, this comparison was performed using three different time intervals to integrate the IPP samples: 1-minute, 2-minute, and also 16-minute intervals. The first one (1 min) implies not using any a priori information, the second one (2 min) uses 1-minute prior information, and the last one uses 15-minute prior information. However, in all the three cases above, the temporal sliding window always advances/shifts 1 min, and thus the resolution is 1 min.

The GPR(VQI) with the SSS cross-validation scheme presented similar error results (Table 6) for the three integration intervals, but the correlation was better for the 16-minute interval, despite larger over and underestimation errors. In the case of the LOGO cross-validation scheme, areas with high density (SJCE/LSJK/SJCU, SJCE, and INCO) have similar results for the different integration intervals, as well as areas with low density (LNTA and LBOA), but with larger errors and significantly lower correlations, as it would be expected.

Koulouri (2022) obtained RMSE = 0.058 and CORR = 0.8573 for the SJCE station not using it in the interpolation, and RMSE = 0.195 and CORR = 0.464 for the LBOA station, not using it in the interpolation. The GPR(VQI) approach proposed here presents lower errors, mainly in areas with low density of stations, besides evaluating the errors for all satellites visible to the test stations, and not just the four satellites of the SJCE station and two of the LBOA station, as shown in that work. The LSJK and SJCU stations are in the same city as the SJCE station. Therefore, tests shown above were conducted leaving out these three stations (SJCE/LSJK/SJCU). Obviously, the results leaving out only the SJCE station were better but were included for the sake of comparison with that work.

5 Scintillation maps generated for the test cases

Considering the main test set of 40 maps, a total of 1920 maps for the 10-hour dataset were generated in sequences of 16 min for each of the 4 approaches and each of the 12 pre-processing options, and are available in the Zenodo repository, as described in the “Data availability statement” section, in the file “Scintillation maps for the 10-hour dataset.zip”.

The comparison of the scintillation maps generated by the four approaches (GDA, IDW, RBF, and GPR), each one with three sets of pre-processing options (SAR, SMR and VQI), appears in Figures 4, 5, and 6, for a moderate, a moderate-to-strong, and a strong scintillation event, respectively. The existing approaches are GDA and IDW, and the proposed ones are RBF and GPR. Above, the baseline pre-processing options are SAR and SMR, and the proposed one is VQI. These three sets of pre-processing options were chosen since SAR and SMR are commonly applied in other works and VQI resulted in the best maps for any of the four approaches. Each map results from the interpolation of IPP samples for an interval of 16 min.

Fig. 4 Scintillation maps generated by the GDA, IDW, RBF, and GPR approaches with each of the SAR, SMR, and VQI pre-processing options for the same 16-minute interval corresponding to a moderate scintillation event starting at 1:16 UT of 2013-11-28, included in the main test set of 40 maps.

Fig. 5 Scintillation maps generated by the GDA, IDW, RBF, and GPR approaches with each of the SAR, SMR, and VQI pre-processing options for the same 16-minute interval corresponding to a moderate-to-strong scintillation event starting at 00:31 UT of 2014-10-26, included in the main test set of 40 maps.

Fig. 6 Scintillation maps generated by the GDA, IDW, RBF, and GPR approaches with each of the SAR, SMR, and VQI pre-processing options for the same 16-minute interval corresponding to a strong scintillation event starting at 00:31 UT of 2014-11-24, included in the main test set of 40 maps.

A comparison plot for all the 40 scintillation maps generated using the 10-hour dataset is available in the Zenodo repository, as described in the “Data availability statement” section, in the file “Set of maps for all 4 approaches with the SAR, SMR, and VQI sets of options.zip”. The generated maps allow to conclude that the GDA approach generated maps with many artifacts and areas with saturation of the S4 values. On the other hand, the RBF approach generated maps with few artifacts, but with areas of saturation of the S4 values. The IDW performed better than GDA and RBF concerning the saturation of S4 values, but presents many artifacts that resulted from undue extrapolation, mainly in areas with low density of GNSS stations that implied in less “interpolation samples”. The maps generated by the GPR approach presented a good balance of smoothness, absence of artifacts, and low extrapolation in areas of low density of GNSS stations. It can be concluded that the GPR approach generated better maps than the other three approaches.

The generated maps and the error and correlations metrics for all four approaches, but only related to the three best pre-processing options (SAR, SMR, and VQI), allow to conclude that the SMR set of pre-processing options tends to saturate the S4 values. Considering the same map, the SMR set generates higher S4 values, followed by SAR and VQI. However, it is worth to note that VQI implies performing the vertical projection of the slant S4 values, which are then attenuated.

The above results for the test set of 40 maps show that the GPR(VQI) is the best approach with the best set of pre-processing options, rendering maps with better smoothness, absence of artifacts, no saturation of S4 values, lower undue extrapolation in areas with low density of GNSS stations, and also lower errors and better correlation. Therefore, GPR(VQI) is the approach proposed in this work.

5.1 Scintillation maps for the additional test set

Additional maps for the same dataset evaluated by Koulouri (2022) were generated using GPR(VQI), being available in the Zenodo repository, as described in the “Data availability statement” section, in the file “Scintillation maps for the 8-hour comparison dataset.zip”, but using three different time intervals to integrate the IPP samples: 1-minute (as in that work), 2-minute and 16-minute, corresponding to the resolutions of 1, 2, and 15 min, respectively. Besides the maps, the corresponding animations of the sequence of maps corresponding to that dataset are available in the Zenodo repository, as described in the “Data availability statement” section, in the file “Animations of GPR(VQI) maps for the 8-hour comparison dataset.zip”, for the three resolutions (1, 2, and 15-minute). Figures 7, 8, and 9 show samples of a GPR(VQI)-generated sequence of maps in the same time intervals of Koulouri (2022), with 1, 2, and 15-minute resolutions, respectively. These maps allow to observe the evolution of ionospheric irregularities on the night of 1st to 2nd of December 2014. These maps are smoother, with the absence of artifacts and absence of areas with the saturation of S4 values, in comparison with the maps of that work, mainly considering the map boundaries and regions with low density of GNSS stations.

Fig. 7 Snapshots ordered in time of GPR(VQI) 1-minute integration interval scintillation maps of the additional test set, showing the evolution of ionospheric irregularities on the night of 1st to 2nd of December 2014.

Fig. 8 Snapshots ordered in time of GPR(VQI) 2-minute integration interval scintillation maps of the additional test set, showing the evolution of ionospheric irregularities on the night of 1st to 2nd of December 2014.

Fig. 9 Snapshots ordered in time of GPR(VQI) 16-minute integration interval scintillation maps of the additional test set, showing the evolution of ionospheric irregularities on the night of 1st to 2nd of December 2014.

6 Discussion and conclusions

This work analyzed two existing and two proposed approaches for generating maps of the S4 amplitude scintillation index, each one being a combination of an interpolation method with some existing and/or proposed sets of pre-processing options. These approaches are named after the related interpolation method as GDA, IDW, RBF, and GPR. The making of scintillation maps requires the interpolation of the IPP samples of each satellite-station link considering the set of GNSS stations. The IPP samples of this work are related to GNSS stations covering most of South America for 10 sets of 1-hour intervals, which were selected in order to present strong scintillation events. Each IPP sample is composed of the IPP location and its S4 value.

The quality of the scintillation maps was assessed using three main criteria: (i) error and correlation metrics, (ii) the smoothness of the map, and also (iii) the absence of artifacts. These artifacts usually appear in the boundaries of the map due to interpolation problems or appear in areas with low density of GNSS stations due to undue extrapolation. The goal is to provide real-time maps that can be actually useful for users and researchers in terms of monitoring and/or mitigating scintillation. This real-time capability was already demonstrated in Martinon et al. (2022a), which proposed and implemented four data processing steps: S4 data acquisition from LISN and INCT networks of GNSS stations, pre-processing, interpolation, and image rendering using a web server.

It was possible to generate in real-time a sequence of scintillation maps, each one integrating 16 min of IPP samples, using a 1-minute sliding window that sweeps these samples, resulting in a 1-minute resolution. The high temporal and spatial variability of the ionospheric scintillation precludes the use of low map integration interval for its effective monitoring. Therefore, the 16-minute integration interval of IPP samples was adopted for the generation of the maps. The GPR approach consistently generated the more accurate maps with any set of pre-processing options, but the GPR(VQI) had the best overall performance concerning the above criteria, and also error and correlation metrics. Considering the set of IPP samples of the cell of a given grid point, it employs the vertical projection of the S4 values, the reduction of these values by the average of S4 values above the 3rd quartile, and the location of the resulting “interpolation sample” given by the centroid of these values.

As extensively demonstrated by the maps generated in this work, the vertical projection renders maps with a more accurate and smoother distribution of S4 values than the equivalent pre-processing option that uses the slant S4 value. The reason is that, for each line-of-sight satellite-receiver, the GNSS signal travels a longer path within ionospheric irregularities than it does for a vertical path, depending on the elevation angle. Thus, the vertical projection normalizes all S4 values of observations with different elevation angles, assuming a 90° line of sight in the IPP.

It is worth to note that an extensive set of tests was performed for the 4 approaches, each one with 12 pre-processing options, using all IPP samples and two validation schemes, resulting in over 20,000 maps for IPP samples of a total of 10 h. IPP samples were divided into mutually exclusive interpolation and test sets. For instance, a single test case using cross-validation, employed almost 18,000 IPP samples to generate a scintillation map, and almost 6000 IPP samples to test/evaluate it. Most of the recently proposed approaches for generating scintillation maps only evaluated a few time series of IPP samples from some PRNs (satellite identifier number), like Amabayo et al. (2021), Hamel et al. (2014), Kieft et al. (2014), and Koulouri (2022). The exception was Vani (2018), but even this one did not perform such an extensive set of evaluations.

Furthermore, except for Koulouri (2022) and Vani (2018), the other approaches employ interpolation by kriging, requiring the pre-calculation of variograms, and thus precluding the generation of maps in real-time. Despite its wide use for interpolation, kriging is just a particular case of Gaussian Process Regression. The proposed GPR-based interpolation embeds an implicit error minimization of the set of interpolation functions, which is automatically performed, making feasible the generation of real-time scintillation maps. This feature was already demonstrated in an actual scenario with GNSS stations, with a personal computer acting as a hub and generating scintillation maps in real time (Martinon et al., 2022a). For instance, Koulouri (2022), performed a real-time download of past scintillation data from an online database, not from actual GNSS stations.

It is worth to note that all evaluations performed here are based on the scintillation maps generated by the interpolation approach with the set of pre-processing options, and not the interpolator itself that is derived from the interpolation points for each map. In other words, the evaluations considered the actual maps that a user would download from a website in order to monitor scintillation-prone areas of the map. Therefore, only S4 values at grid points of the generated map are considered. The evaluation of the S4 error of the generated map is made by comparing the S4 of each IPP sample of the test set with the corresponding S4 value in the map. However, since the map contains S4 values only at grid points, such comparison requires performing an interpolation of the S4 values at the four neighboring grid points of the IPP sample location.

The availability of real-time scintillation maps with better quality would provide better data for new approaches of scintillation prediction. However, there is another future work cited in Martinon et al. (2022b), the generation of scintillation probability maps, which is being currently implemented and tested, in order to refine specific parameters of its algorithm. The generation of scintillation probability maps will complement the current real-time maps proposed in the manuscript.

The proposed approach GPR(VQI) is simple, robust, accurate, and well-proven for the considered set of data for the South American territory. The related computer code will be made available soon in a GitHub repository. Its implementation conforms to the real-time generation of scintillation maps, even considering a 1-minute resolution. As soon as all related GNSS stations are configured for such real-time requirements, it is expected to have better monitoring of scintillation over Brazil.

Acknowledgments

Authors thank the National Laboratory for Scientific Computing (MCTIC/LNCC, Brazil) for the use of the Santos Dumont supercomputer, as part of the project “Implementation of machine learning approaches that demand supercomputing for ionospheric scintillation prediction and meteorological forecast of severe convective events”. Authors take part and acknowledge the Brazilian multi-institutional project INCT “GNSS Technology for Supporting Aerial Navigation”. Author André R. F. Martinon acknowledges the National Council for Scientific and Technological Development (CNPq, Brazil) for the doctorate grant 140562/2020-6. The GNSS-NavAer INCT Project is supported by grants CNPq 465648/2014-2, FAPESP 2017/50115-0, and CAPES 88887.137186/2017-00. Author Eurico R. de Paula acknowledges the support from CNPq under process number 302531/2019-0. The authors acknowledge the LISN (Low-latitude Ionospheric Sensor Network) for the data employed in this work. The authors thank the Brazilian Ministry of Science, Technology, and Innovation, and the Brazilian Space Agency. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Data availability statement

Supplementary data pointed out in Sections 4 and 5 which includes the complete error evaluation, additional maps, and animations of sequences of maps, is available in a specific Zenodo repository at https://doi.org/10.5281/zenodo.7438643. The complete error evaluation table for all the 12 pre-processing options employed for each interpolation method, mentioned in Section 4, appears in the file “Complete table of interpolation errors and correlation.csv” of the Zenodo repository. Considering the main test set of 40 maps of Section 5, the total set of 1920 maps for the 10-hour dataset was generated in sequences of 16 min, for each of the 4 approaches, and each of the 12 pre-processing options, is also available in the same Zenodo repository, in the file “Scintillation maps for the 10-hour dataset.zip”. Also considering Section 5, the comparison plot for all the 40 scintillation maps generated using the 10-hour dataset is also available in the same repository, in the file “Set of maps for all 4 approaches with the SAR, SMR and VQI sets of options.zip”. Section 5.1 states that additional maps for the same dataset evaluated by Koulouri (2022) were generated using GPR(VQI), appearing in the Zenodo repository, in the file “Scintillation maps for the 8-hour comparison dataset.zip”, as well as the corresponding animations of the sequence of maps, in the file “Animations of GPR(VQI) maps for the 8-hour comparison dataset.zip”, for the three specified resolutions.

Appendix

References

All Tables

All Figures

	Fig. 1 Geometry of a satellite-ground station line-of-sight depicting the Ionospheric Pierce Point.
In the text

	Fig. 2 Geographic location of the 45 GNSS stations of the LISN and INCT networks that provided scintillation data for this work.
In the text

Fig. 3 Bar charts showing the distribution of the IPP samples according to the scintillation class for the 10 one-hour non-consecutive selected intervals of time. Each bar chart depicts the date and starting time of a one-hour interval that is further divided into four 16-minute intervals for generating the maps (with a 1-minute overlap between them). The corresponding values of the geomagnetic index Kp and the adjusted F10.7 solar flux index are also included for each one-hour interval.

In the text

	Fig. 4 Scintillation maps generated by the GDA, IDW, RBF, and GPR approaches with each of the SAR, SMR, and VQI pre-processing options for the same 16-minute interval corresponding to a moderate scintillation event starting at 1:16 UT of 2013-11-28, included in the main test set of 40 maps.
In the text

	Fig. 5 Scintillation maps generated by the GDA, IDW, RBF, and GPR approaches with each of the SAR, SMR, and VQI pre-processing options for the same 16-minute interval corresponding to a moderate-to-strong scintillation event starting at 00:31 UT of 2014-10-26, included in the main test set of 40 maps.
In the text

	Fig. 6 Scintillation maps generated by the GDA, IDW, RBF, and GPR approaches with each of the SAR, SMR, and VQI pre-processing options for the same 16-minute interval corresponding to a strong scintillation event starting at 00:31 UT of 2014-11-24, included in the main test set of 40 maps.
In the text

	Fig. 7 Snapshots ordered in time of GPR(VQI) 1-minute integration interval scintillation maps of the additional test set, showing the evolution of ionospheric irregularities on the night of 1st to 2nd of December 2014.
In the text

	Fig. 8 Snapshots ordered in time of GPR(VQI) 2-minute integration interval scintillation maps of the additional test set, showing the evolution of ionospheric irregularities on the night of 1st to 2nd of December 2014.
In the text

	Fig. 9 Snapshots ordered in time of GPR(VQI) 16-minute integration interval scintillation maps of the additional test set, showing the evolution of ionospheric irregularities on the night of 1st to 2nd of December 2014.
In the text