Methodology to estimate ionospheric scintillation risk maps and their contribution to position dilution of precision on the ground

Satellite-based communications, navigation systems and many scientific instruments rely on observations of trans-ionospheric signals. The quality of these signals can be deteriorated by ionospheric scintillation which can have detrimental effects on the mentioned applications. Therefore, monitoring of ionospheric scintillation and quantifying its effect on the ground are of significant interest. In this work, we develop a methodology which estimates the scintillation induced ionospheric uncertainties in the sky and translates their impact to the end-users on the ground. First, by using the risk concept from decision theory and by exploiting the intensity and duration of scintillation events (as measured by the S4 index), we estimate ionospheric risk maps that could readily give an initial impression on the effects of scintillation on the satellite-receiver communication. However, to better understand the influence of scintillation on the positioning accuracy on the ground, we formulate a new weighted dilution of precision (WPDOP) measure that incorporates the ionospheric scintillation risks as weighting factors for the given satellite-receiver constellations. These weights depend implicitly on scintillation intensity and duration thresholds which can be specified by the end-user based on the sensitivity of the application, for example. We demonstrate our methodology by using scintillation data from South America, and produce ionospheric risk maps which illustrate broad scintillation activity, especially at the equatorial anomaly. Moreover, we construct ground maps of WPDOP over a grid of hypothetical receivers which reveal that ionospheric scintillation can also affect such regions of the continent that are not exactly under the observed ionospheric scintillation structures. Particularly, this is evident in cases when only the Global Positioning System (GPS) is available.

amplitude scintillation (i.e. temporal fluctuation of the signal strength) and phase scintillation occur often at lower latitudes and can cause periods of reduced signal power at the receiver's antenna (fading). Therefore, GNSS signal-to-noise-ratio may drop below the receiver-tracking threshold [20] thereby leading to loss of positioning information, that can persist for a significant period of time [6].
Many characteristics of scintillation are well-studied [1,7,11,28,8,20,17,44]. It has been shown that scintillation activity varies with operating frequency, local time, season, geomagnetic activity level(planetary K-index, K p ), and the 11-year solar cycle [46,29,38,7,28,20]. Moreover, the scintillation effects vary according to the user's location (i.e. geomagnetic latitude). During high scintillation activity, these effects can be a serious problem for GNSS users in certain regions, such as polar and equatorial latitudes, while users in other regions (e.g. middle latitudes) are not affected to a large extent [48,2,20,31]. Equatorial regions experience strong scintillation mostly after sunset until a few hours after midnight local time [39,20,21]. In equatorial ionosphere, the most significant scintillation activity occurs with deep signal fades causing a GNSS receiver to lose lock and degrade positioning accuracy [47,6,36,33].
Limited case studies have been carried out related to the effects of scintillation for the end-user [60,36,49,58,41,3,31,15,18] and often global climatological models fail to describe extreme day-today variability of scintillation [57]. Moreover, because scintillation is frequency dependent and cannot be eliminated, for example by combining observable signals at different frequencies [32,54], the constant monitoring of scintillation activity using in-situ measurements alongside with the development of mathematical methodologies, which describe, quantify and interpret its impact on GNSS communications, is substantial to mitigate its effect. In this framework, we develop a mathematical methodology to estimate ionospheric scintillation risks based on scintillation measurements and decision theory.
These risks are given by the joint probability of scintillation intensity and duration estimated using information obtained from ground based scintillation monitors. Thus, in contrast to the commonly estimated statistics (e.g. averages or percentages of scintillation measurements above a threshold as in [48,28,49,19]), we use data collected from a network of scintillation monitors to construct ionospheric risk maps that describe and monitor scintillation behavior. While the ionospheric maps can readily be used to get an (initial) impression of the effects of scintillation on a satellite-receiver communication, in this work, we propose to incorporate the scintillation risk into the receiver's position dilution of precision (PDOP) [27].
Standard PDOP is a well-known way to quantify the reliability and integrity of a GNSS positioning system by using the geometry of the available satellite constellations [45,59,9,50,52]. Here, we show how the degradation of PDOP occurs relative to the location of ionization structures in the F-region due to scintillation. To do that, we propose to use a weighted position dilution of precision (WP-DOP) which is estimated by assigning different weights to receiver-satellite links using the ionospheric scintillation risks along the corresponding line-of-sights. Hence, in addition to the number and the geometry of the visible satellites, we take into account the effect of scintillation. Previous works have employed empirical error tracking models, e.g. Conker model [11], to implicitly introduce underlying scintillation activity by weighting receiver's measurements with the inverse of the variance of the tracking error. However, these approaches are often susceptible to model limitations and requirements of heuristically selected parameters and their applicability is limited due to the requirement of specific software/hardware that is available mainly in operating ground stations and not in standard GNSS receivers [5,13]. Instead, our approach utilizes weighing factors based on risks associated directly with monitored scintillation data. With the help of these weights, we are able to construct maps that describe the expected uncertainty in positioning results during ionospheric scintillation. To the best of our knowledge, this is the first time such a methodology is introduced.
The rest of the paper is organized as follows. In section 2, we present the theoretical background and mathematical expression of the risks and the proposed WPDOP. In section 3, we describe the formulae and technical details for the implementation of the proposed ionospheric risks and the corresponding ground WPDOP values. In section 4, we demonstrate our methodology by estimating ionospheric risk maps over South America using real scintillation data. Subsequently, we construct ground maps based on two different set-ups using either only GPS or GPS, GLONASS and Galileo.
Then, we discuss the main observations and further potentials of the proposed methodology. Finally, in section 5, we present a summary of our work and future aims.

Risk concept
While it may be straightforward to extract statistics (e.g. averages or standard deviations) to describe scintillation [30], nevertheless it is important to interpret those statistics for an end-user, in terms of the costs expected for a particular application when scintillation is experienced. For this purpose, we use risk analysis from decision theory [37]. Risk can be defined as the expected loss incurred when an activity is performed (for example satellite-receiver communication). To estimate the risk we need i) to find features that describe the condition in which we perform this activity (e.g. communication under scintillation can be described through the measured S 4 value), ii) the frequency of occurrence of these features and iii) how injurious the condition described by these features is to the successful outcome of the activity. These are respectively expressed by i) a feature vector ii) a probability and iii) a loss function.
In this section, we describe the background theory and definitions about loss function and risk before making it application-specific in the next section. Here, we denote random variables by capital letters and their realizations by lowercase letters. We start with the measurable feature Z ∈ Z , where Z represents the set of all possible feature vectors (e.g. scintillation index, etc).
Let us introduce a 0-1 loss function which implies there is no loss when values of z are less than a specified threshold, but otherwise losses are equally injurious. Although this is a simple loss function, it nevertheless is instructive for this particular problem without increasing computational complexity, and as will be shown it permits a straightforward interpretation.
The risk is then defined as the expected value of the loss function given by where E[.] denotes the expectation value and p(z) is the probability density funciton of the measurable feature Z. From (eq. 1) and (eq. 2), we have that where π(.) denotes probability. Thus, the risk has a simple interpretation in terms of probability.

Risk associated with ionospheric scintillation
In this section, we apply the previous analysis to estimate risk maps using scintillation features.
Let X denote a scintillation index (e.g. S 4 value). Moreover, we define as a scintillation event a time series of scintillation index values above a scintillation threshold denoted by x th over a duration limit d th (see Fig. 1). The duration of scintillation events can be treated as a random variable denoted by D. Based on the analysis of section 2.1, we have that (X, D) ≡ Z.
Using the insight that the vulnerability of a link depends on the value of the scintillation index (e.g. strong scintillation when S 4 > 0.6) and duration of a scintillation event, we define the loss function where x th and d th are the intensity and duration thresholds respectively. These thresholds can be chosen to reflect the specifications of the receiver and the needs of safety critical applications.
From (eq. 3), the risk due to scintillation features is Hence, the risk of the underlying activity, depends on the joint probability of the duration and intensity of scintillation events and takes values between 0 and 1.

Effect of scintillation on the positioning accuracy
The proposed risk can readily be used to give an initial impression of the effect of scintillation in a satellite-receiver communication. However, to better understand the influence of scintillation on the positioning accuracy on the ground, a well-known quantity as the dilution of precision can be estimated [27].
Towards this aim, this section introduces a weighted position dilution of precision (WPDOP) as a measure of uncertainty in estimating the receiver's position. To put our analysis in perspective, we first discuss the theoretical basis of the standard position dilution of precision (PDOP) from the Bayesian point of view. Then, we explain how in addition to the availability and geometry of the satellites, the proposed WPDOP can include information about the effect of ionospheric scintillation.

Standard position dilution of precision
Following similar analysis as in [27], for a ground receiver at location r = (r x , r y , r z ) T , the positioning error at time t, denoted by ∆r = (∆r x , ∆r y , ∆r z ) T , can be expressed through the linear where b ∈ R S is a vector with the differences between the measured and modelled pseudorange values, where n s is a unit column vector pointing from the modelled (approximated) receiver position to the satellite and S is the total number of visible satellites at a time instance t.
Error ε ∈ R S represents the measurement noise plus model errors and ionospheric effects (including scintillation). Further details about the derivation of the aforementioned linear system are given in In the PDOP standard probabilistic analysis [27,35], the noise ε is modelled as i.i.
where p(∆r|b) is the posterior distribution of ∆r given by p(∆r|b) ∝ p(b|∆r) π(∆r) based on the } is the likelihood and p(∆r) is the prior distribution [23].
With an uninformative prior p(∆r), the positioning error can be obtained by ∆r := arg max ∆r p(b|∆r) which is the maximum likelihood (ML) estimate. The ML estimate is and it is equal to the least squares solution of system (eq. 6). The solution (eq. 8) exists as long as there are 3 or more satellites located in distinct directions in the sky. Quantitatively, the rank of (A T A) can be used as an indicator whether the problem is well posed or not. For example, if rank(A T A) ≥ 3, then we can deduce that there are enough satellites for the position estimation. However, the accuracy of the solution (eq. 8) depends on the inverse matrix (A T A) −1 (whether it is well or ill-conditioned).
From the Bayesian point of view, this inverse matrix multiplied by error variance γ is the posterior covariance of p(∆r|b) given by Γ = γ(A T A) −1 and its diagonal elements (posterior variances) can help to quantify the expected accuracy of the estimate ∆r [23]. Large posterior variances may indicate that there is possibly an overlapping between satellites (eclipse phenomenon), collinearity or inappropriate geometric constellation. The square root of the trace of the posterior Γ normalized by the constant variance γ is the well-known position dilution of precision (PDOP) [27]. In particular, for a ground 1 We note here that the clock offset is omitted (as a variable to be estimated) from the current formulation. Matrix A, ∆r, b and ε are time-varying variables but for simplicity in the notation, time t has not been used in this context. Also, linear system (eq. 6) helps us to introduce the different dilution of precision formulae. The estimation of vector b is out of the scope of the current work, the interested reader is referred to for example [53]. receiver at location r and time t, the PDOP is given by where tr() denotes the trace of matrix given by tr ( Hence, PDOP can be used to quantify the expected uncertainty in the estimate (eq. 8).

Weighted position dilution of precision
Now, by considering the scintillation effect independently from all the other errors, the noise term in the observation model (eq. 6) can be decomposed into two uncorrelated errors The error due to ionospheric scintillation can be approximated by a Gaussian distribution given by ε sc ∼ N (0, Γ εsc ) [35] where its covariance is a diagonal matrix Γ εsc = diag(γ (s) sc ) with the variances depending on the scintillation risks along the available satellite-receiver ray-paths . Therefore, the error variances due to scintillation can be expressed as a function of the ionospheric risk r Function h(r  The remaining error ε rem , that encloses all the other modelling and measurement related uncertainties, is modelled as Gaussian, i.e. ε rem ∼ N (0, γI S×S ) with γ constant as in the standard PDOP analysis 2 [27]. We notice that in the current analysis the errors between different links are considered uncorrelated which is a plausible assumption since there are usually large distances between different ray-paths 3 .
Then, the covariance matrix of total error ε given by Γ ε = γI S×S + Γ εsc is a diagonal matrix The positioning error estimate is now when S ≥ 3 and rank(A T Γ −1 ε A) ≥ 3. The posterior covariance of ∆r is (14) and the weighted position dilution of precision (WPDOP) for a ground receiver at location r and time t is defined (in a similar way as in eq. 9) as the square root of the trace of the posterior Γ normalized by the constant γ and it is given by where tr ( where 0 ≤ w s ≤ 1 and γ constant. A rigorous derivation of (eq. 15) is given in appendix Appendix B.
Based on equation (eq. 16) we have that w s = 1 when r The shape of h() in (eq. 16) may be selected based on prior knowledge or application requirements. In general, the proposed WPDOP uses both information about the number of available satellites, their geometry in the sky and the risks associated with scintillation along the visible ray-paths. If rank(A T W A) < 3 i.e. the number of available satellites is less than 3 then WPDOP r,t is undefined.

Single shell ionospheric model
The F-region (especially around 350 km of altitude) is of main concern for GNSS users due to propagation errors caused by the dynamics of the large electron density. To construct surface ionospheric maps it is standard to consider the ionospheric thin shell approximation at 350 km altitude, based on the assumption that the entirety of the electron content (which affects a link) is acting at the point where a ray-path of a satellite-receiver intersects (or pierces) the ionospheric shell at that altitude [14]. Hence, an ionospheric risk map can be estimated using ground measurements that are directly projected on the ionospheric shell at the corresponding ionospheric pierce points (IPPs) 4 (similarly as in [55]). The projected (scintillation) measurements at the 350km ionospheric shell will be referred to as IPP data in the following text.

Ionospheric computational domain and data accumulation
For the construction of an ionospheric risk map, the thin ionospheric shell was discretized using a uniform grid where each pixel got its own risk value. In the results section, we estimated scintillation risk maps over South America and as a reliable index of scintillation at low-latitudes, we used S 4 data which is the ratio of the standard deviation of the received signal intensity (measured by a 50 Hz sampling rate) to the averaged intensity during an 1 minute time interval. A risk was estimated at each ionospheric pixel v based on the joint probability (eq. 5) using IPP data accumulated over a time interval τ . We note that there is an implicit assumption of statistical stationarity i.e. the joint probability does not change over time.

A scintillation event
A scintillation event was introduced briefly in section 2. Here, we explain more analytically the properties of such events. Particularly, duration of a scintillation event is the time interval where the scintillation values for a specific (satellite-scintillation monitor) link are greater than a given threshold.
In the schematic example of Figure 1, we show how different scintillation events are determined using sequences of IPP data samples (depicted by small circles) for 3 different links imagined that their ray-paths are passing through the same ionospheric pixel v. The horizontal dashed line indicates a selected intensity threshold x th . All the links experience scintillation events since there are data points above the threshold. For link 1 (in blue), we can observe that there are two distinctive scintillation events between [t 2 − t 6 ] and [t 11 − t 14 ]. However, for links 2 (in yellow) and 3 (in black) we have to deal with data gaps related either to strong scintillation or other reasons. In these cases, we have to decide whether the available data and the corresponding gaps are considered as a single scintillation event or not. In the current implementation, we consider that there is a single scintillation event if the gap is small (i.e. less than 4 minutes). For a longer gap, we consider that two separate scintillation events of shorter durations take place. Extra conditions based on the values of the scintillation samples before and after the gap may also be applied to differentiate single or multiple events. For example very high values of the S 4 samples on either side of a large gap (e.g. S 4 ≥ 0.6) can imply a single event rather than two distinct events.

Risk estimation
In this section we describe the mathematical formulae for the estimation of the ionospheric risk maps. As we mentioned earlier, risk maps are constructed by estimating pixel-wise the joint probability (eq. 5) using the accumulated IPP data at each pixel of the discretized ionospheric domain over a period of time τ . Based on that and (eq. 5) we can mathematically express a risk in a pixel as the conditional probability given pixel v and time window τ , i.e.
Moreover, we introduce the following notation for clarity. A scintillation measurement is given by is the geographic location (at the ionospheric pierce point) at time t. Super-script s is used to identify that the measured index corresponds to a particular satellitemonitor link. Then, a scintillation event is a time series of scintillation samples Here, ∆ is the sampling period and d th is an integer denoting the minimum number of samples required to qualify for the event (a.k.a. duration threshold). Thus, each IPP sample can be represented pairwise ) that includes the scintillation value and the number of samples in the scintillation event to which x s (v(t)) belongs. For example, d s (t) = 0 when x s (v(t)) < x th which means that sample x s (v(t)) does not belong to a scintillation event. On the other hand when x s (v(t)) ≥ x th then sample x s (v(t)) belongs to a scintillation event and then d s (t) equals the number of consecutive samples (before and after(and including) x s (v(t)) at time t) that exceed threshold x th .
In the following analysis, for simplicity in the probability expressions, we write measurement is the value of the scintillation index at time t j for satellite-receiver link s and the corresponding IPP (ionospheric location) of The joint probability of X and D (eq. 17) can be re-written as The previous equation can be further expanded, i.e.
where d k is an integer denoting the number of samples ( Based on equation (eq. 19), risk (eq. 17) is estimated by first calculating probability π(X ≥ x th |v, τ ) according to where U is a step function defined as index S is the total number of available links, α is the intensity threshold and N s the number of the available scintillation samples (i.e. total number of measurements) for link s over a time window τ .
Furthermore, the conditional probability where function H d is defined as and gives either 1 or 0 based on the sequence Finally by estimating pixel by pixel a risk (eq. 17) at each v, we can obtain a map in matrix form denoted by

WPDOP using scintillation risks
To estimate WPDOP (eq. 15) first we compute matrix A ∈ R S×3 which includes the relative positions between a ground receiver and S available satellites following the details of appendix Appendix A.
Regarding the weights, without further knowledge about the precise shape of function h() (the properties of which we anticipate may be gleaned from either statistical or physical information, or alternatively selected based on the requirements of a specific application), we approximate w s (eq. 16) by where r (s) ion is the ionospheric risk value at the intersection point between the line-of-sight of (link) s and the ionospheric risk map at 350km. Weight (eq. 25) satisfies the limiting values for w s (eq. 16) as risk varies and parameter k ≥ 0 penalizes our trust to a ray-path based on the scintillation risk value.
For the estimation of WPDOP in the result section, we set k = 2 to ensure that relatively high trust is given to ray-paths that are associated with low scintillation risks.
So, given a hypothetical ground receiver and S available satellites, different weights are assigned to the ray-paths (between the receiver and the satellites) based on the points the ray-paths pierce the ionospheric map. Then, WPDOP (eq. 15) is estimated using the weighting matrix Figure 3.b explains schematically how to assigns a weight to each each ray-path using the ionospheric risk map.

Data, map construction, results and discussion
In this section, we present ionospheric risk maps and corresponding WPDOP ground maps over the area of South America at geographic latitude between ∼ −40 • N and ∼ 10 • N and longitude between ∼ −90 • E and ∼ −35 • E. This region experiences significant scintillation activity that causes deep signal fades inducing a GNSS receiver to lose lock of one or more satellite signals and includes the equatorial anomaly [4].

S 4 data
Scintillation activity maximizes mostly after sunset until a few hours after midnight local time during the equinoctial months in the equatorial ionosphere [16]. For this reason in this study we used available S 4 data between 23.00 and 03.00 Universal Time (UT). The S 4 data (that was used to produce a risk map) was measured from a network of 38 ISMR scintillation monitors using all the available GPS, GLONNASS and Galileo satellites over the 4 hour period.
In particular, the S 4 measurements ( with sampling period ∆ = 1 minute) and their corresponding IPPs at 350km were provided from the CIGALA/CALIBRA network -UNESP web server [10].

Construction of ionospheric maps
For the ionospheric risk maps, a uniform standard grid of spatial resolution 2 • × 2 • was used for the estimation of the risk in each ionospheric pixel (as suggested also by [51]). A scintillation risk (joint probability eq. 17) was estimated in each pixel using accumulated S 4 data over the interval 23.00-03.00 (UT) (for example Figure 2.b shows the available data over the 4 hours interval for a set of nine pixels). Ionospheric pixels that did not have any IPP data were left blank in the following ionospheric maps.
For clarity in the subsequent text, the night time interval (23.00-03.00 UT) crossing from 02 to 03 November 2014 will be simply referenced as that of 02 November 2014, and similarly the night time from 03 to 04 November 2014 will be referenced as that of 03 November 2014.

Construction of ground maps
The spatial resolution of the ground maps was selected as 1 • × 1 • (in latitude and longitude) by placing hypothetical receivers (denoted by green triangles) in the centre of each pixel as shown in Figure 3.c. To calculate the WPDOP and PDOP for a set of hypothetical ground receivers, the positions of the satellites in view have to be estimated. There are different ways that this can be done using for example orbital information. In our study, the locations of the available satellites were computed geometrically at time t using downloaded data from [10].  Moreover, we can see that there are two bands of scintillation i.e. north and south of the geomagnetic equator. This corresponds to the well known equatorial anomaly [24]. The risk of scintillation in the north band seems to disappear for increasing thresholds relatively drastically (see second and third row in Figure 4) which indicates that scintillation of shorter duration and lower intensity prevails in this ionospheric region. We note that this could also be related to the number of the ground scintillation monitors which are less in the north band compared to the south band (so the south band may give us more reliable estimates).

Ionospheric risk maps over South America
In Figure 5 scintillation structures remain similar as in Figure 4, however, we can notice that the scintillation in Figure 4 covers a broader area of the south band compared to Figure 5 where high scintillation activity is observed predominantly over Brazil. Overall, we can say that scintillation of low intensity (0.3 ≤ S 4 < 0.5) affects a broad area of South America with risks approaching 0.8 whereas scintillation risk drops drastically for S 4 ≥ 0.7 during both days.

Ground maps over South America
We illustrate the effect of scintillation on the ground by estimating the proposed weighted dilution of precision in a set of hypothetical ground receivers (small green triangles in Figure 3 In Figure 6 (and similarly in Figure 7), along column 1 we can observe the ionospheric risk maps; along column 2 the maps of WPDOP estimated based on the risk maps (of column 1); along column 3 the map of PDOP (we note that maps (c) and (g) are identical and have been repeated to ease the comparisons); and along column 4 the contribution of scintillation to the total dilution of precision (in a percentage form). In Figure 6.i, PDOP values in the ground map (c) and (g) are low and increase slightly towards the west coast; in particular in the area marked with the red ellipse the PDOP values are higher than in the corresponding area in Figure 6.ii (map (c) and (g)) due to less satellite coverage. and (f)) the WPDOP values are relatively low even though the risk due to scintillation over that region is significant. This is most likely because the ray-paths of the satellites in view (that were used to estimate the WPDOP at 01.00.09 (UT)) do not intersect the scintillation active region at this WPDOP ground map (f) in Figure 6.i (column 2 and second row) reveals that high S 4 threshold (i.e S 4 ≥ 0.5) reduces the effect of scintillation on the WPDOP values. This is because the strong scintillation activity is more focal based on the corresponding risk map (e) in Figure 6.i. However, by comparing the percentage maps (d) and (h) along the column 4 of Figure 6.i, we can observe that the loss in dilution of precision is localized similarly both when S 4 ≥ 0.3 and S 4 ≥ 0.5. Both (WPDOP-PDOP)/WPDOP maps of column 4 (in Figure 6.i) give a percentage of at least 30% contribution from scintillation which can be observed for longitude between ∼ −60 • E and ∼ −75 • E and latitude ∼ −20 • N and ∼ −25 • N (area marked with yellow ellipses in map (d) and (h)). This suggests that scintillation has a significant influence in the position accuracy and has to be taken into account. Furthermore, we can notice that further away from the ionospheric scintillation structures, the scintillation contribution to WPDOP drops to zero.
With the introduction of extra constellations, the effect of scintillation is reduced broadly with a small exception in the north part of the continent as we can observe in WPDOP maps (b) and (f) along column 2 of Figure 6.ii (area marked with an ellipse in red color). This happens due to scintillation and not due to poor geometry since the PDOP values of the PDOP map (c) and (g) in column 3 of Furthermore, test case of Figure 7 demonstrates that the effect of scintillation activity can yield different ground structures. This is happening because the geometry of the available constellations and the risk maps are different between Figure 6 and Figure 7. The maps (d) and (h) along column 4 in Figure 7 reveal "shadowing effects" on the Earth due to scintillation around the ionospheric active region, with the contribution decreasing with increasing distance from the scintillation structures.
Similarly as in the example of Figure 7 we can conclude that the effect of scintillation on the ground is reduced when more satellites are in use and when the high risk regions in the ionospheric map are more localized.
(i) Ionospheric risk map estimated using S 4 data over 23.00-03.00 (UT) and ground maps using only GPS constellation at 01.00.09 (UT) (ii) Ionospheric risk map using S 4 data over 23.00-03.00 (UT) and ground maps using GPS, GLONASS and Galileo constellation at 01.00.09 (UT) Figure 6: Risk maps and ground maps for different sets of satellite constellations for 02 November 2014. Ionospheric static risk maps were estimated using data over the interval 23.00-03.00 (UT). The geometry of the constellations for the ground maps was that of time instant 01.00.09 (UT). Within (i) or (ii) each row corresponds to different thresholds in S4. From left to right, column 1 shows ionospheric risk maps for thresholds in S4 of 0.3 and 0.5, and scintillation duration of 5 minutes, column 2 illustrates WPDOP ground maps which explicitly incorporate the scintillation risks as weights from ionospheric map in column 1, column 3 gives the PDOP ground maps (with no scintillation being included), column 4 give us the difference between WPDOP and PDOP, expressed as a percentage of WPDOP. The regions that are described in the text are marked with colored ellipses.
(i) Ionospheric risk map estimated using S 4 data over 23.00-03.00 (UT) and ground maps using only GPS constellation at precisely 01.00.09 (UT) (ii) Ionospheric risk map using S 4 data over 23.00-03.00 (UT) and ground maps using GPS, GLONASS and Galileo constellation at 01.00.09 (UT) Figure 7: The results are presented as in Figure 6. Ionospheric risk maps and ground maps for different sets of satellite constellations for 03 November 2014.

Current implementation and intrinsic limitations
In the current application of our methodology, we were restricted by the number of available ground scintillation monitors (unfortunately, at the moment the ground monitors are very few and unevenly spread over the continent) and the fixed number of satellites. So, in effect these limitations forced us to produce ionospheric scintillation risk maps using data collected over a relatively long time window. The selection of the 4 hour time window (period of the highest scintillation activity) was used to capture the average scintillation activity. Therefore, the estimated static risk maps described the average scintillation activity during the specific time window.
A WPDOP value was estimated for a specific locations using the number of available satellites at that moment and the risk map, as shown in Figure 3.b. The weights were assigned to different ray-paths based on the ionospheric scintillation risk maps and were directly related to statistical information about scintillation since the risks that were used to estimate the weights (eq. 26) had been computed from scintillation data (i.e. scintillation intensity and duration which has been extracted from S 4 measurements). Therefore, the proposed WPDOP included the impact of the average scintillation activity.
In the future, if more scintillation monitors are built on the ground then our methodology can produce ionospheric scintillation maps with higher temporal resolution (say, every 30 or 15 minutes).
This will allow to update the weights in the WPDOP constantly and to observe the impact of the dynamic scintillation activity. Currently, we can say that a moving window and/or a shorter time window for the data collection (to construct ionospheric risk maps) could be implemented only in areas that are covered by a dense network of scintillation monitors. By using shorter data collection intervals and considering for example lower scintillation thresholds, we could produce scintillation risk maps useful in aircraft navigation. These risk maps can update the weights in the PDOP and provide information to the pilot whether to trust the navigation system.

Comparing WPDOP with PDOP
The dilution of precision is used as a quality measure in positioning systems since it expresses the expected uncertainty in the positioning estimates. Roughly speaking, a low dilution of precision implies high confidence level in the positioning results, and a high dilution of precision indicates low confidence level.
In general, we can notice that the standard PDOP considers implicitly the scintillation effects in the sense that strong scintillation can lead to loss of one or more satellites resulting in high measured PDOP values (due to poor satellite geometry of the remaining satellites). However, by using the standard PDOP it is not possible to quantify the impact that low-to-high scintillation activity has on the accuracy of a positioning system when losses of lock do not occur. In such cases, the degradation of the GNSS positioning may be significant even though the satellite constellation geometry would appear promising. The concept of WPDOP allows to estimate the impact that scintillation induced errors can have on PDOP by assigning the proposed scintillation weights in all the available satellite-receiver links.
Overall, we observed that the more satellites were used in the estimation of the dilution of precision, the lower were the values both for PDOP and WPDOP and hence lower position errors were expected 5 .
So, the combination of GPS, GLONASS and Galileo can improve precision. However, there were still some areas which were affected by scintillation (especially for weak scintillation of duration 5 minutes or more) even when all the available constellations were used. This suggests that by increasing the number of satellites in view, the positioning precision can be further improved in the case of strong localized scintillation. This outcome is in line with a recent work by [34] and the review article [26].

Novelty and transferability
In this work, instead of using basic statistical analysis, we exploited the concept of risk from decision theory to quantify the effects that the ionospheric scintillation can have on potential enduser activities. In particular, we derived risks considering that the quality of a satellite-receiver communication can be characterized only by scintillation measurements and then we constructed ionospheric risk maps using S 4 data provided by the UNESP server in South America. As we showed in section 4.4 this methodology can be used for the construction of sequences of ionospheric risk maps ( Figure. 4and 5). Such maps could be further used as a-priori information for ionospheric stochastic models, for enriching scintillation climatology data-bases, or for improving understanding of scintillation phenomena over multiple years [42,22]. Moreover, by explicitly incorporating scintillation risks into dilution of precision, one can overcome intrinsic limitations of empirical models and receiver specifications. Finally, the proposed methodology allows to produce high resolution ground map which could potentially be used to help mitigate the effects of ionospheric scintillation on GNSS positioning.
The proposed methodology for the construction of ionospheric risk and ground maps is reusable and transferable. Notably, if one wants to tailor risk to specific applications, then this methodology has the advantage that its parameters can be optimized and validated for specific applications. For example the loss function (eq. 4) allows us to select scintillation intensity and duration thresholds, x th and d th respectively, based on the specifications of the receiver and the needs of safety critical applications. By employing a 0-1 loss function, we showed that risks can be easily interpreted in terms of probability of scintillation events. However, not only the thresholds, but also alternative loss functions to equation (eq. 4) can be adopted. Moreover, our risk formalism allows to consider extra (measurable) inputs in feature Z (of section 2.1) in addition to scintillation data X and durations of events D, such as information about data gaps and other ionospheric quantities (e.g. total electron content (TEC)). Therefore, the risks (eq.5) (as joint probabilities of multiple parameters) can be estimated to illustrate an overall ionospheric activity. Regarding the weighting factors w s (eq. 16), different shapes of the error variance function h(.) (eq. 11) can be adopted based on the outcomes of future research without any changes in the proposed methodology of section 2.3.
Alternatively, the construction of multiple risk maps using different sets of inputs (measurable features), such as scintillation values measured in different frequencies which reflect how small-scale plasma density irregularities in F region affect different propagating signals and gradients of the TEC which can capture large scale irregularities, can offer versatile and complimentary information to 5 For further details about the link between WPDO/PDOP and position error see Appendix Appendix C the users. As a further step, these sets of risk maps can infer implicit statistical information about different (multi-frequency) errors which can be considered in the position error model (eq. 6) and the estimation of the weights in the WPDOP. Hence, the concept of WPDOP can be extended to include other ionospheric-related errors to inform how different propagation errors arising from ionospheric structures can impact the positioning on the ground.
In particular, our formalism enables the adaptation of the WPDOP given the presence of scintillation or any other ionospheric-related errors (e.g. in the equatorial F-region). For example the error term in the observation model (eq. 6) could be decomposed into more uncorrelated error terms e.g. ε = ε sc + ε ∆TEC + ε rem and the corresponding ionospheric risks could be used to provide information about the different error variances (expressed as weights) for the WPDOP, e.g.
where g(r (s) ∆TEC ) is the error covariance due to ∆T EC ionospheric effects. Therefore, based on the available statistical knowledge, a WPDOP calculated by any hypothetical receiver could depend upon the combination of scintillation levels on multiple frequencies (utilized for positioning), the presence of losses of lock on some or all the frequencies associated with a given link, other ionospheric errors (e.g. TEC) and the geometry of the available links. Then by comparing WPDOP against PDOP, we can infer information about the impact that different ionospheric irregularities can have on the positioning accuracy on the ground.

Conclusions
This paper proposed a methodology which connects the ionospheric uncertainties in the sky originating from small-scale irregularities in the F-region (e.g. in the equatorial ionosphere) with their impact on the positioning on the ground. In particular, first we exploited decision theory to develop risks defined as expected losses that incur during a satellite communication activity which is assumed to be fully characterized by ionospheric measurable features, in our case scintillation intensity and duration. The proposed scintillation risks relied on a 0-1 loss function and thus they were equal to the joint probability of scintillation intensity and duration above specified thresholds. The derived risk formulation was used to estimate pixel-wise risk values in an uniform ionospheric grid at 350 km and the estimated risk maps reflected the ionospheric plasma density inhomogeneities causing scintillation during the observed time window.
We demonstrated our methodology by estimating risk maps in the sky using accumulated data obtained from GPS, GLONASS and Galileo scintillation monitors in South America during local nighttime hours and using different scintillation intensity and duration thresholds. Overall, we observed that the two bands of scintillation i.e. north and south of the geomagnetic equator were visible in cases where the intensity threshold did not exceed value 0.5 and when the duration of the scintillation events was up to 10 minutes. Subsequently, to understand the effect of scintillation on the ground, we proposed to use a weighted position dilution of precision (WPDOP) which was estimated by assigning different weights to receiver measurements from different satellites using the ionospheric scintillation risks along the line-of-sights (LoSs) of the corresponding links. Thus, the dilution of precision on the ground depended upon the combination of the geometry of the available links and the scintillationinduced observable errors (quantified though the scintillation risks) on the same available links.
Finally, we constructed instantaneous WPDOP ground maps combining information from the ionospheric risk maps and the available link geometries. The WPDOP maps revealed those areas on the ground which were more affected by ionospheric scintillation. We noted that the scintillation effects on the ground receiver were more prominent when only GPS was available. However, the effects were not totally eliminated with the addition of extra constellations.
The present work focused on the effect of ionospheric scintillation to the dilution of precision. In the future, other sources of errors or modelling uncertainties can be taken into account in the dilution of precision estimations (i.e. superposition of risks associated, not only with scintillation but also with other ionospheric measurable features, can be considered). Also, we expect to further study how errors due to scintillation limit the accuracy in position estimation both from theoretical and applied perspective and how statistical information about these errors can be incorporated into existing precise point positioning softwares.

Author contribution statement
Methodology for ionospheric maps: AK, NDS and BCV with support from the other co-authors; Theory for ground maps: AK, NDS and VR; Implementation of ionospheric maps: AK and BCV; Implementation of ground maps: AK; Analyzing the results: AK and NDS; Writing the paper: AK as main author in collaboration with NDS and VR and with editorial support from BCV, IA and BF; Scoping the research and ongoing discussion: BF and IA.

Data availability statement
The S 4 data used to estimate the ionospheric maps in section 4.4 was downloaded from http:// is-cigala-calibra.fct.unesp.br; The MATLABa R codes for the estimation of the ionospheric and ground maps can be downloaded from the link https://github.com/AlexandraKoulouri/Ionospheric-Scintillation-Maps-and-PDOP. .  and ε ∈ R S . Here, we assume that even in the absence of a loss of lock, scintillation can induce higher-order errors in the observables that do not cancel out in dual frequency combinations [5,54].

Appendix B. Defining the position dilution of precision
As we will show in the following text, the position dilution of precision is a measure of the uncertainty of the estimates ∆r and it is parameterized according to the statistical characterization of the errors ε and the linear model used for the position estimates (eq. 6). In particular, by denoting Γ xx , Γ yy and Γ zz the diagonal elements of the general covariance Γ = (A T Γ −1 ε A) −1 (eq. 14), the standard deviations of ∆r (eq. 13) in X,Y and Z direction are σ x = √ Γ xx , σ y = Γ yy and σ z = √ Γ zz respectively. The dilution of precision is defined as the normalized (with a common scaling factor κ) root mean square (RMS) of the standard deviations given by RMS = σ 2 x + σ 2 y + σ 2 z = Γ xx + Γ yy + Γ zz = tr(Γ) = tr((A T Γ −1 ε A) where W = Γε κ −1 . Based on that ratio, we can see that the weights can be derived based on the prior knowledge about the error statistics. Furthermore, from the expression of the PDOP XY Z we can understand that the dilution of precision is not the actual positioning error, but a measure of the uncertainty of the estimates ∆r. Therefore, the key feature of the dilution of precision is that it is a single value that can be used as an indicator about the reliability of a positioning system. Moreover, we can notice that as matrix A and W do not depend on the measurements, but only on the geometry and the weighting scheme, dilution of precision can be computed from the satellite orbital information without needing real measurements. Currently, there are empirical tables that relate PDOP values with the expected accuracy of a positioning system (see for example table 1 in [56]).
In the following text, we will see that the scaling factor κ is defined based on the model of Therefore, the common scaling factor is κ = γ and W = I S×S . Then, we have the well-known standard position dilution of precision given by (eq. 9).
For an observation error modelled as in (eq. 10 of section 2.3.2), the error covariance has two terms, i.e. as in (eq. 16). Thefore, based on the previous error analysis, PDOP is given by

Appendix C. Average position error and position dilution of precision
Let us first clarify that ∆r = r tr − r 0 is the position difference between the true XYZ position vector, r tr = (r xtr , r ytr , r ztr ), and the computed (approximate) one, r 0 = (r x 0 , r y 0 , r z 0 ), and ∆r =r − r 0 is the position vector difference between the estimated position vectorr = (r x ,r y ,r z ) and r 0 .
Based on the linear system (eq. 6) and the observation error modelled using a Gaussian distribution i.e. ε ∼ N (0, Γ ε ), the estimated position is given bŷ