Issue 
J. Space Weather Space Clim.
Volume 13, 2023



Article Number  18  
Number of page(s)  17  
DOI  https://doi.org/10.1051/swsc/2023015  
Published online  16 June 2023 
Technical Article
A new approach for the generation of realtime GNSS lowlatitude ionospheric scintillation maps
^{1} Coordination for Applied Research and Technological Development (COPDT/INPE), National Institute for Space Research (INPE), Av. dos Astronautas, 1758, 12227010 São José dos Campos, São Paulo, Brazil
^{2}
Department of Heliophysics, Planetary Sciences and Aeronomy (DIHPA/INPE), National Institute for Space Research (INPE), Av. dos Astronautas, 1758, 12227010 São José dos Campos, São Paulo, Brazil
^{*} Corresponding author: andre.martinon@inpe.br
Received:
21
December
2022
Accepted:
15
May
2023
Ionospheric scintillation disturbs radio frequency signals affecting navigation based on Global Navigation System (GNSS), especially in Brazil, due to the large magnetic declination. The generation of realtime scintillation maps is an important way to provide scintillation monitoring. This work considers amplitude scintillation, given by the S4 index. Some existing and some proposed approaches for generating these maps are presented and tested, being each one a combination of an interpolation method with some existing and/or proposed sets of preprocessing options. These approaches are named after the related interpolation method as GRIDDATA, Inverse Distance Weighting, Radial Basis Functions and Gaussian Process Regression. The making of scintillation maps requires the interpolation of Ionospheric Pierce Point (IPP) samples, given by the S4 values for each IPP of each satellitestation link considering the set of GNSS stations of the given area and time interval. Some intervals of time that presented strong scintillation over Brazil were selected and the corresponding sets of IPP samples were used to obtain sequences of maps using all possible combinations of interpolation and preprocessing options. Furthermore, a fifth, more recent, approach was also included in the comparison. The quality of the resulting maps was assessed, concluding that the Gaussian Process Regression approach, with a specific set of preprocessing options, allows to generate the most accurate scintillation maps. The proposed map generation approach is part of a broader proposal being implemented to provide realtime scintillation maps covering the Brazilian territory.
Key words: Ionospheric scintillation map / Interpolation / Gaussian process regression / Realtime data acquisition / Realtime image processing
© A.R.F. Martinon et al., Published by EDP Sciences 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 satellitestation 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 preprocessing 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 realtime with a 1minute 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 preprocessing 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/scihome), 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 preprocessing options, and also an improved set of preprocessing 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 preprocessing options, allows the generation of the most accurate maps for the considered data. Each map integrates IPP samples of a 16minute interval (since there is 1minute 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 1minute 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 realtime scintillation maps covering the Brazilian territory, even with a 1minute 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 1minute 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 satellitestation 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 preprocessing options was originally employed in the IDW approach, requiring the following assumptions/steps:
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);

Define the IPP for each satellitereceiver link, as shown in Figure 1, as being the intersection of the lineofsight receiversatellite (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 satelliteground station lineofsight depicting the Ionospheric Pierce Point.
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);
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 preprocessing options as well as the newly proposed ones were applied to the four mapgeneration approaches: (1) GDA, (2) IDW, (3) RBF, and (4) GPR. The GDA is a triangulationbased 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 thinplate 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 greatcircle 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 preprocessing options. The first modification was not to consider the S4 measured along the receiversatellite 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 GNSSNavAer network (Vani et al., 2017) stations at a 1minute rate, while for the stations of the LowLatitude 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 1minute 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 preprocessing 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, threeletter abbreviations were adopted to denote the preprocessing 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 preprocessing options of the former approaches are SMR or SAR. The newlyproposed 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 IDWinterpolation approach with the preprocessing options SMR.
Therefore, each map generation approach was tested with each possible set of preprocessing 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 multiconstellation. 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).
LISN GPS ionospheric monitoring stations coordinates.
INCT GNSS ionospheric monitoring stations coordinates.
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).
Unified S4 data file format for each S4 observation.
Scintillation data was acquired for 10 selected nonconsecutive onehour intervals that mostly presented moderate and strong scintillation events. Each onehour interval was taken from a different night between November 2013 and November 2014. Each one of the 10 onehour 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 1minute overlap between consecutive maps. Figure 3 presents the distribution of the S4 values classified according to the scintillation severity for each onehour 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 Lband 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 onehour nonconsecutive selected intervals of time. Each bar chart depicts the date and starting time of a onehour interval that is further divided into four 16minute intervals for generating the maps (with a 1minute overlap between them). The corresponding values of the geomagnetic index Kp and the adjusted F10.7 solar flux index are also included for each onehour interval. 
4 Error and correlation results
The presented 12 preprocessing options and 4 interpolation methods were evaluated for the set of 40 sixteenminute intervals extracted from the 10 nonconsecutive onehour 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 onehour dataset is expanded by an extra minute, being split into 4 sixteenminute intervals with 1minute 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 distancebased 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 mutuallyexclusive sets of IPP samples are obtained from the original set of IPP samples described in the previous section (10 onehour sets), which is split according to a validation scheme.
Two crossvalidation 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 10fold scheme was chosen, thus the IPP samples in a given 16minute 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 16minute 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 sixteenminute intervals, the 10fold SSS scheme results in 400 scintillation maps, for each of the 4 approaches and each of the 12 preprocessing options, yielding a total of 19,200 maps (40 × 10 × 4 × 12). Each iteration of the 10fold scheme (for each approach and each set of preprocessing 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 bestconsidered set of preprocessing options (VQI) and 5 LOGO options, according to the leftout 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 16minute interval, there are 16 samples of a fixedposition 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 preprocessing 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 preprocessing 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 preprocessing 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 3rdquartile reductions (SQR and SQI) consistently presented higher errors and lower correlations than the remaining VQR/VQI, SAR/SAI, and VAR/VAI, with the corresponding verticallyprojected 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 preprocessing 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.
GPR interpolation error metrics and correlation results for all preprocessing options using the SSS crossvalidation scheme for the main test set of 40 maps.
Concerning the reduction function of the S4 values of the IPP samples, the Q preprocessing 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 preprocessing options, it can be noticed that the V preprocessing 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 crossvalidation results are shown in Table 5 shows the results for the 4 interpolation approaches and only for the VQI set of preprocessing 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 GNSSstation density areas, the GPR and IDW results are similar, and better than those of GDA and RBF. For lowdensity areas, the GPR approach performed better, proving its robustness.
Interpolation error metrics and correlation of all four approaches for the VQI set of preprocessing options using the LOGO crossvalidation scheme (left column shows station or set of stations left out from the interpolation, but used for evaluation), for the main test set of 40 maps.
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 1minute 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 crossvalidation 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: 1minute, 2minute, and also 16minute intervals. The first one (1 min) implies not using any a priori information, the second one (2 min) uses 1minute prior information, and the last one uses 15minute 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 crossvalidation scheme presented similar error results (Table 6) for the three integration intervals, but the correlation was better for the 16minute interval, despite larger over and underestimation errors. In the case of the LOGO crossvalidation 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.
GPR(VQI) error metrics and correlation for the three integration intervals using the crossvalidation schemes SSS (“All”) and LOGO (left column shows station or set of stations left out from the interpolation, but used for evaluation), for the additional test set of 475 maps.
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 10hour dataset were generated in sequences of 16 min for each of the 4 approaches and each of the 12 preprocessing options, and are available in the Zenodo repository, as described in the “Data availability statement” section, in the file “Scintillation maps for the 10hour dataset.zip”.
The comparison of the scintillation maps generated by the four approaches (GDA, IDW, RBF, and GPR), each one with three sets of preprocessing options (SAR, SMR and VQI), appears in Figures 4, 5, and 6, for a moderate, a moderatetostrong, and a strong scintillation event, respectively. The existing approaches are GDA and IDW, and the proposed ones are RBF and GPR. Above, the baseline preprocessing options are SAR and SMR, and the proposed one is VQI. These three sets of preprocessing 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 preprocessing options for the same 16minute interval corresponding to a moderate scintillation event starting at 1:16 UT of 20131128, 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 preprocessing options for the same 16minute interval corresponding to a moderatetostrong scintillation event starting at 00:31 UT of 20141026, 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 preprocessing options for the same 16minute interval corresponding to a strong scintillation event starting at 00:31 UT of 20141124, included in the main test set of 40 maps. 
A comparison plot for all the 40 scintillation maps generated using the 10hour 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 preprocessing options (SAR, SMR, and VQI), allow to conclude that the SMR set of preprocessing 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 preprocessing 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 8hour comparison dataset.zip”, but using three different time intervals to integrate the IPP samples: 1minute (as in that work), 2minute and 16minute, 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 8hour comparison dataset.zip”, for the three resolutions (1, 2, and 15minute). 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 15minute 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) 1minute 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) 2minute 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) 16minute 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 preprocessing 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 satellitestation 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 1hour 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 realtime maps that can be actually useful for users and researchers in terms of monitoring and/or mitigating scintillation. This realtime 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, preprocessing, interpolation, and image rendering using a web server.
It was possible to generate in realtime a sequence of scintillation maps, each one integrating 16 min of IPP samples, using a 1minute sliding window that sweeps these samples, resulting in a 1minute 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 16minute 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 preprocessing 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 preprocessing option that uses the slant S4 value. The reason is that, for each lineofsight satellitereceiver, 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 preprocessing 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 crossvalidation, 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 precalculation of variograms, and thus precluding the generation of maps in realtime. Despite its wide use for interpolation, kriging is just a particular case of Gaussian Process Regression. The proposed GPRbased interpolation embeds an implicit error minimization of the set of interpolation functions, which is automatically performed, making feasible the generation of realtime 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 realtime 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 preprocessing 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 scintillationprone 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 realtime 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 realtime maps proposed in the manuscript.
The proposed approach GPR(VQI) is simple, robust, accurate, and wellproven 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 realtime generation of scintillation maps, even considering a 1minute resolution. As soon as all related GNSS stations are configured for such realtime 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 multiinstitutional 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/20206. The GNSSNavAer INCT Project is supported by grants CNPq 465648/20142, FAPESP 2017/501150, and CAPES 88887.137186/201700. Author Eurico R. de Paula acknowledges the support from CNPq under process number 302531/20190. The authors acknowledge the LISN (Lowlatitude 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 preprocessing 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 10hour dataset was generated in sequences of 16 min, for each of the 4 approaches, and each of the 12 preprocessing options, is also available in the same Zenodo repository, in the file “Scintillation maps for the 10hour dataset.zip”. Also considering Section 5, the comparison plot for all the 40 scintillation maps generated using the 10hour 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 8hour comparison dataset.zip”, as well as the corresponding animations of the sequence of maps, in the file “Animations of GPR(VQI) maps for the 8hour comparison dataset.zip”, for the three specified resolutions.
Appendix
Acronyms (following order of appearance in the text).
References
 Abdu MA. 2019. Daytoday and shortterm variabilities in the equatorial plasma bubble/spread F irregularity seeding and development. Prog Earth Planet Sci 6: 11. https://doi.org/10.1186/s4064501902581. [CrossRef] [Google Scholar]
 Alken P, Thébault E, Beggan CD, Amit H, Aubert J, et al. 2021. International geomagnetic reference field: The thirteenth generation. Earth Planets Space 73: 49. https://doi.org/10.1186/s4062302001288x. [NASA ADS] [CrossRef] [Google Scholar]
 Amabayo EB, Andima G, Ssenyunzi RC. 2021. Instantaneous ionospheric scintillation mapping over the east African region by use of GPS derived amplitude scintillation proxy. Asian J Res Rev Phys 4(2): 6–20. https://doi.org/10.9734/ajr2p/2021/v4i230138. [CrossRef] [Google Scholar]
 Bandyopadhayay T, Guha A, DasGupta A, Banerjee P, Bose A. 1997. Degradation of navigational accuracy with global positioning system during periods of scintillation at equatorial latitudes. Electron Lett 33(12): 1010–1011. https://doi.org/10.1049/el:19970692. [CrossRef] [Google Scholar]
 Broomhead DH, Lowe D. 1988. Multivariable functional interpolation and adaptive networks. Complex Systems 2: 321–355. [Google Scholar]
 de Paula ER, Kherani EA, Abdu MA, Batista IS, Sobral JHA, et al. 2007. Characteristics of the ionospheric Fregion plasma irregularities over Brazilian longitudinal sector. Indian J Radio Space Phys 36(4): 268–277. http://nopr.niscpr.res.in/handle/123456789/4706. [Google Scholar]
 de Paula ER, Martinon ARF, Moraes AO, Carrano C, Neto AC, Doherty P, et al. 2021. Performance of 6 different global navigation satellite system receivers at low latitude under moderate and strong scintillation. Earth Space Sci 8: e2020EA001314. https://doi.org/10.1029/2020EA001314. [CrossRef] [Google Scholar]
 de Rezende L, de Paula ER, Kantor I, Kintner P. 2007. Mapping and survey of plasma bubbles over brazilian territory. J Navig 60(1): 69–81. https://doi.org/10.1017/S0373463307004006. [CrossRef] [Google Scholar]
 Hamel P, Sambou DC, Darces M, Beniguel Y, Hélier M. 2014. Kriging method to perform scintillation maps based on measurement and GISM model. Radio Sci 49: 746–752. https://doi.org/10.1002/2014RS005470. [CrossRef] [Google Scholar]
 Humphreys TE, Ledvina BM, Psiaki ML, Cerruti AP, Kintner PM. 2004. Analysis of ionospheric scintillations using wideband GPS L1 C/A signal data. In: Proceedings of the 17th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2004), Long Beach, CA, September 2004, pp. 399–407. [Google Scholar]
 Koulouri A. 2022. Realtime ionospheric imaging of S_{4} scintillation from limited data with parallel Kalman filters and smoothness. IEEE Trans Geosci Remote Sens 60: 4106012. https://doi.org/10.1109/TGRS.2022.3140600. [CrossRef] [Google Scholar]
 Kieft P, Aquino M, Dodson A. 2014. Using ordinary kriging for the creation of scintillation maps. In: Mitigation of ionospheric threats to GNSS, Notarpietro R, Dovis F, Franceschi GD, Aquino M, (Eds.) IntechOpen, Rijeka. https://doi.org/10.5772/58781. [Google Scholar]
 Kintner PM, Kil H, Beach TL, de Paula ER. 2001. Fading timescales associated with GPS signals and potential consequences. Radio Sci 36(4): 731–743. https://doi.org/10.1029/1999RS002310. [CrossRef] [Google Scholar]
 Kintner PM, Ledvina BM, de Paula ER. 2005. An amplitude scintillation test pattern for evaluating GPS receiver performance. Space Weather 3: S03002. https://doi.org/10.1029/2003SW000025. [Google Scholar]
 Leick A. 1995. GPS satellite surveying, 2nd edn, John Wiley and Sons, USA, ISBN10. ISBN10: 0471306266 / ISBN13: 9780471306269. [Google Scholar]
 Martinon ARF, Stephany S, de Paula ER. 2022a. Challenges in realtime generation of scintillation index maps. Colóquio Brasileiro de Ciências Geodésicas/V Simpósio Brasileiro de Geomática, Curitiba, PR, Brazil, November 2022. Extended abstract. https://cbcg.ufpr.br/wpcontent/uploads/2023/02/AnaisCBCG_SBG2022VersaoFinalCompactada.pdf. [Google Scholar]
 Martinon ARF, Stephany S, de Paula ER. 2022b. Ionospheric scintillation awareness for GPS and other constellations by probability maps. Proceedings of the XXIII Brazilian Symposium on GeoInformatics 23: 275–280. ISSN 21794847. São José dos Campos, SP, Brazil, November 2022. Short paper. http://urlib.net/ibi/8JMKD3MGPDW34P/487MADB [Google Scholar]
 Prol FS, Camargo PO, Muella MTAH. 2017. Comparative study of methods for calculating ionospheric points and describing the GNSS signal path. Bol Ciênc Geod 23(4): 669–683. https://doi.org/10.1590/S198221702017000400044. [CrossRef] [Google Scholar]
 Rasmussen CE, Williams CKI. 2016. Gaussian processes for machine learning. The MIT Press, USA. ISBN: 026218253X. [Google Scholar]
 Rino CL. 1979. A power law phase screen model for ionospheric scintillation: 1. Weak scatter. Radio Sci 14(6): 1135–1145. https://doi.org/10.1029/RS014i006p01135. [CrossRef] [Google Scholar]
 Skone S, Knudsen K, Jong M. 2001. Limitations in GPS receiver tracking performance under ionospheric scintillation conditions. Phys Chem Earth 26: 613–621. https://doi.org/10.1016/S14641895(01)001107. [CrossRef] [Google Scholar]
 Spogli L, Alfonsi L, de Franceschi G, Romano V, Aquino MHO, Dodson A. 2009. Climatology of GPS ionospheric scintillations over high and midlatitude European regions. Ann Geophys 27: 3429–3437. https://doi.org/10.5194/angeo2734292009. [CrossRef] [Google Scholar]
 Spogli L, Alfonsi L, Romano V, de Franceschi G, Monico G, Shimabukuro M, Bougar B, Aquino M. 2013. Assessing the GNSS scintillation climate over Brazil under increasing solar activity. J Atmos Sol Terr Phys 105–106: 199–206. https://doi.org/10.1016/j.jastp.2013.10.003. [CrossRef] [Google Scholar]
 Valladares CE, Chau JL. 2012. The lowlatitude ionosphere sensor network: Initial results. Radio Sci 47(04): 1–18. https://doi.org/10.1029/2011RS004978. [CrossRef] [Google Scholar]
 van Dierendonck AJ, Klobuchar J, Hua Q. 1993. Ionospheric scintillation monitoring using commercial single frequency C/A code receivers. In Proceedings of the 6th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GPS 1993), Salt Lake City, UT, September 1993, 1333–1342. [Google Scholar]
 Vani BC, Shimabukuro MH, Monico JFG. 2017. Visual exploration and analysis of ionospheric scintillation monitoring data: the ISMR query tool. Comput Geosci 104: 125–134. https://doi.org/10.1016/j.cageo.2016.08.022. [CrossRef] [Google Scholar]
 Vani BC. 2018. Investigações sobre modelagem, mitigação e predição de cintilação ionosférica na região brasileira, Tese (doutorado). Universidade Estadual Paulista. Faculdade de Ciências e Tecnologia, Presidente Prudente. Available at: http://hdl.handle.net/11449/153701. [Google Scholar]
Cite this article as: Martinon ARF, Stephany S & de Paula ER 2023. A new approach for the generation of realtime GNSS lowlatitude ionospheric scintillation maps. J. Space Weather Space Clim. 13, 18. https://doi.org/10.1051/swsc/2023015
All Tables
GPR interpolation error metrics and correlation results for all preprocessing options using the SSS crossvalidation scheme for the main test set of 40 maps.
Interpolation error metrics and correlation of all four approaches for the VQI set of preprocessing options using the LOGO crossvalidation scheme (left column shows station or set of stations left out from the interpolation, but used for evaluation), for the main test set of 40 maps.
GPR(VQI) error metrics and correlation for the three integration intervals using the crossvalidation schemes SSS (“All”) and LOGO (left column shows station or set of stations left out from the interpolation, but used for evaluation), for the additional test set of 475 maps.
All Figures
Fig. 1 Geometry of a satelliteground station lineofsight 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 onehour nonconsecutive selected intervals of time. Each bar chart depicts the date and starting time of a onehour interval that is further divided into four 16minute intervals for generating the maps (with a 1minute overlap between them). The corresponding values of the geomagnetic index Kp and the adjusted F10.7 solar flux index are also included for each onehour 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 preprocessing options for the same 16minute interval corresponding to a moderate scintillation event starting at 1:16 UT of 20131128, 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 preprocessing options for the same 16minute interval corresponding to a moderatetostrong scintillation event starting at 00:31 UT of 20141026, 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 preprocessing options for the same 16minute interval corresponding to a strong scintillation event starting at 00:31 UT of 20141124, included in the main test set of 40 maps. 

In the text 
Fig. 7 Snapshots ordered in time of GPR(VQI) 1minute 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) 2minute 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) 16minute 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.