Summary
A fast Monte Carlo technique to simulate equatorial ionospheric scintillation on global navigation satellite system signals is proposed. The algorithm uses a single-layer phase-screen model of the ionosphere and the scintillation is expressed as a Huygens-Fresnel integral (HFI). By assuming a specially-tailored random phase screen, the HFI can be expressed in closed form as a combination of Fresnel integrals. We statistically characterize the amplitude and phase computed by the HFI for different values of the scintillation index S4. Results for the L1, L2, and L5 bands were obtained and compared with real data, showing good agreement. Some of the advantages of the proposed technique are: (a) the amplitude and phase of the scintillation process are simultaneously obtained; (b) arbitrarily long ionospheric scintillation time series with pre-defined stationary characteristics are synthesized; and (c) several scintillation time series corresponding to different carrier frequencies are generated using a common phase-screen model.
1 INTRODUCTION
Fast varying electron density irregularities cause random variations in the phase and amplitude of a wave propagating through the ionosphere. This effect, called ionospheric scintillation, generates deep intensity fadings and random phase shifts in global navigation satellite system (GNSS) received signals. Strong scintillation tends to increase the positioning errors and may even prevent the receiver operation (Conker et al., 2003; Humphreys et al., 2005). Ionospheric scintillation occurs, in general, in the polar regions and in the regions approximately 20 degrees north and south of the magnetic equator. The two types of scintillation have different characteristics as they are provoked by different physical phenomena (Ghafoori & Skone, 2015). Scintillation observed in the equatorial region is generally stronger, equally affecting the intensity and phase of the received signals with a maximum between sunset and midnight (ITU, 2019). In this contribution, we are mainly interested in simulating equatorial scintillation.
Equatorial amplitude scintillation affects both code and phase tracking, thus degrading pseudorange and carrier-phase measurements in GNSS receivers. Equatorial phase scintillation adversely affects the operation of the receiver’s phase-locked loop (PLL) and may lead to carrier cycle slips, navigation bit errors, and even the complete loss of carrier lock (Conker et al., 2003; Zhang et al., 2010). Several signal processing techniques can be used to improve carrier tracking robustness, including the frequency-lock loop assisted PLL and the adaptive-bandwidth PLL (Lin et al., 2014; Zhang et al., 2010). Scintillation is often characterized by a pair of indices: the amplitude scintillation index S4 (standard deviation of the received power normalized by its mean value) and the phase scintillation index σϕ (standard deviation of the de-trended carrier phase). Scintillation strength is typically classified as weak (S4 < 0.25), moderate (0.25 ≤ S4 < 0.7), or strong (S4 ≥ 0.7).
In spite of various campaigns that have been undertaken in the last decade to collect data from GNSS signals disturbed by scintillation, the amount of available data may not be sufficient or have the desired features to evaluate the performance of GNSS receivers. Recall that the peaks of solar activity (which lead to frequent events of strong ionospheric scintillation) are periodic phenomena spaced by approximately 11 years. Therefore, it is useful to have access to a Monte Carlo simulation tool that generates synthetic data simultaneously affected by amplitude and phase ionospheric disturbances in various equatorial scintillation scenarios. Note that amplitude is often easier to simulate than phase due to the lack of reliable phase measurements and models (Jiao et al., 2018), especially for intense scintillation in which the phase exhibits excursions outside the interval [−π, π] (Xu & Morton, 2018) and includes canonical fades, which are deep power fades accompanied by fast phase variations of approximately half of a cycle (Humphreys et al., 2010). Nevertheless, having reliable models for the phase induced by the ionosphere is crucial to test the robustness of GNSS receiver phase tracking loops (Conker et al., 2003; Ghafoori & Skone, 2015; Humphreys et al., 2005; Vilà-Valls et al., 2020).
There are essentially three different ways to produce such synthetic data (Vilà-Valls et al., 2015): physics-based models, statistical models, and phase-screen models. Physical models, such as the split-step (Jiao et al., 2018; Rino, 2011), although accurate, are in general complicated to use as they have to be tailored to each scenario. Statistical models like the Cornell Scintillation Model (Humphreys et al., 2009) are easy to apply but, since they are not supported by a physical model, the relation between the amplitude and the phase processes is not justified on physical grounds. Computer simulation of radio wave propagation through one-dimensional random phase screens has become an important tool in the study of scintillations as simulations provide values for the amplitude and phase of a monochromatic wave observed, for instance, on the ground by solving the Huygens-Fresnel integral (HFI; Beach & Lovelace, 1997). Moreover, they are able to reproduce the canonical fading phenomenon, which is considered the main effect that leads to the loss of lock in GNSS receivers (Humphreys et al., 2010). The goal of the phase screen model is to replace the interaction of the radio wave with an extended medium of scatterers by that with a thin phase-changing screen. Although being a simplification of the reality, the single phase screen tends to approximate quite well the behavior of multiple phase screens distributed through a 200-km thick ionospheric F-region (Carrano et al., 2013).
In general, the evaluation of the HFI requires numerical integration which leads to computationally heavy simulations. In Nunes et al. (2018), the authors proposed an analytic method to generate amplitude and phase samples of an ionospheric scintillation random process based on the phase-screen model. The method is computationally efficient as it requires solely the determination of sine and cosine Fresnel integrals. The goal of the present work is to statistically characterize the ionospheric scintillation process presented in Nunes et al. (2018). It is shown that the statistical distributions of the amplitude and phase vary widely with the scintillation index S4. The amplitude distribution is difficult to characterize, although, for strong scintillation, there is evidence that the amplitude is well approximated by a generalized gamma (or alpha-mu) distribution (Moraes et al., 2013; Nunes & Sousa, 2018). The phase distribution is monomodal for weak and moderate scintillation regimes and tends to be multimodal for the strong scintillation regime with S4 > 0.7 because the phase range typically exceeds 2π. A simple and robust phase unwrapper is developed to determine the phase from the modulo-2π version provided by the argument of the HFI.
The paper is organized as follows. Section 2 describes the proposed phase-screen model and presents expressions for the amplitude and phase of the received signal. Section 3 characterizes the amplitude and phase fading processes. Section 4 exhibits the most significant statistical results obtained with the proposed technique and Section 5 summarizes the main achievements of the contribution including their power spectra. Finally, the Appendix describes the adopted fast computation algorithm of the complex Fresnel integral.
2 PHASE-SCREEN MODEL
Plasma instabilities developing in the equatorial ionosphere after sunset can generate large-scale depletions in the electron density called equatorial plasma bubbles, which are accompanied by intense random density variations occurring over a broad range of spatial scales (plasma turbulence). A well-developed bubble structure is associated with irregularity scales ranging from dozens of centimeters to hundreds of kilometers, with hundreds of meter sizes responsible for very high frequency/ultra high frequency (VHF/UHF) scintillations (Abdu, 2005). The multiple-phase screen method is an efficient numerical method for simulating the propagation of radio waves through extended random media that allows both amplitude and phase fluctuations to accumulate within the medium (Carrano et al., 2011).
Herein, we consider a single phase-screen model for equatorial scintillation that assumes that the profile of the total electron content (TEC), which is the number of electrons in a tube of a 1 − m2 cross-section extending from the receiver to the satellite (Misra & Enge, 2006), is concentrated into a thin layer, typically at the altitude of the ionosphere peak electron density (z ≈ 350 km). The random medium can be thought to consist of a number of plasma density depletion areas (plasma bubbles). The scale size of the plasma depletion areas is approximately 10–100 km in the zonal direction and they are greatly elongated and aligned in the magnetic north-south direction (Kintner et al., 2004). The bubbles move from west to east (zonal direction) with typical drift velocities in the range of 25 to 200 m/s (Jiao et al., 2016; Kintner et al., 2004). When a bubble starts to grow or move upward, large electron density gradients on the bubble edges produce smaller irregularities (Ghafoori, 2012). The ionospheric irregularities that generate scintillation have a structure with a scale size of the Fresnel length 
, where λ is the wavelength of the incident wave (Kintner et al., 2004, 2007). For instance, for GPS L1 signals, the Fresnel length is rF ≈ 365 m. At time t, the carrier-phase advance, as compared with the received carrier phase in the absence of ionosphere, is given by (Klobuchar, 1996):
1
where c is the speed of light, f0 is the carrier frequency, and TEC(x, t) is expressed in units of electrons/m2. The Huygens-Fresnel diffraction theory permits us to determine the change effects in the amplitude and phase of a monochromatic wave as a function of horizontal position, x0, on an observing plane located at distance z from the one-dimensional phase screen for small-angle scattering. Assuming z ≫ λ, the wave complex amplitude is given by (Beach & Lovelace, 1997; Psiaki et al., 2007):
2
where A0 depends on the power of the satellite transmitted signal, k = ω0 / c = 2π / λ is the signal wave number, and ω0 = 2πf0. The quantity ϕsc(x, t) includes the phase introduced by the phase screen. The scintillation scenario is illustrated in Figure 1.
Geometry of the phase-screen scintillation scenario (adapted from Ghafoori [2012])
Consider for simplicity that the receiving antenna is at position x0 = 0; define the Huygens-Fresnel integral (HFI) with γ = k / (2z) = π / (λz), according to:
3
such that Equation (2) can be expressed as:
4
The infinite interval of integration in Equation (3) can be reduced to [−R, R], with R finite, without significant loss of accuracy, provided that γR2 ≫ ϕΔ, where ϕΔ = 80.6πΔTEC / (cf0) and ΔTEC stands for the maximum value of the varying part of TEC. This leads to 
. For instance, if z = 350 km, f0 = 1.572 GHz, and ΔTEC = 1 TECU, then R ≫ 337.9 m.
When the ionosphere is homogeneous with TEC(x, t) = TEC0, the HFI is simplified to:
5
with ϕ0 = 80.6πTEC0 / (cf0). The integral above was calculated in terms of the complex Fresnel integral (Abramowitz & Stegun, 1972):
6
where C(y) and S(y) are, respectively, the cosine and the sine Fresnel integrals. The plot of Ξ(y) in the complex plane (Cornu spiral), with Ξ(−y) = −Ξ(y), is shown in Figure 2(a) for |y| < 5. The curve converges to ±(1 + j) / 2 when |y| → ∞. The polynomial approximations described in the Appendix allow us to compute the Fresnel integrals efficiently.
Plots of (a) Ξ(y) = C(y) + jS(y) for |y| < 5 (Cornu spiral) and (b) |IA(a, b; γ)| for γ = 4.7 ×10−5
Consider the integral:
7
which may be expressed in terms of the complex Fresnel integral as:
8
with the asterisk denoting conjugate. Figure 2(b) sketches the modulus of IA(a, b; γ) for z = 350 km and the L1 band, yielding γ = 4.7 × 10−5.
For the homogeneous ionosphere, using Equation (4) and Equation (5) leads to:
9
which shows that, when TEC(x, t) is not varying with t, the amplitude and phase of Isc (0, t) + jQsc (0, t) remain constant and there is no scintillation effect. In the sequel, we model the ionosphere disturbance as a moving train of tubes (electron density irregularities) with triangular cross-sections according to:
10
where TEC0 corresponds to the absence of ionospheric electron density irregularities (homogeneous ionosphere) and ΔTEC is the maximum fluctuation of the ionospheric electron content. In Equation (10), {αi, i = 1, …, M} with 0 ≤ αi ≤ 1 is a set of independent random variables and the triangle function is defined as Λ(x) = 1−|x|, for |x| < 1 and 0, otherwise. Typically, the values of TEC0 vary from 1 to 103 TECU (Klobuchar, 1996), where the TEC unit (TECU) is equal to 1016 electrons/m2. The triangular shape of the tubes’ cross-section in Equation (10) may be considered a reasonable assumption as: (a) it is the simplest continuous shape for the electron density depletion regions and (b) it leads to a major simplification of the analytic expressions used by the HFI. Nevertheless, no physical evidence on the triangular shape was provided by the consulted bibliography.
We assume that the receiver is located at x0 = 0 and each triangle in Equation (10) drifts with a constant velocity vdrift, according to xi(t) = xi(0) + vdriftt. This model is a simplification of the case in reality because it considers that a frozen TEC(x) profile drifts eastward without any deformation (Yeh & Liu, 1982). Nonetheless, this frozen-in assumption is often used in simulations because it simplifies the calculations and yields results that keep a reasonable degree of realism (Psiaki et al., 2007). Notice that, besides the actual velocity of the plasma density irregularities, the effective velocity vdrift depends also on the components of the satellite and receiver velocities along the x-direction (Knight & Finn, 1998).
Using Equation (1) and Equation (10), we obtain for the phase induced by the diffraction at position x:
11
Therefore, only the varying part of the total electron content associated with ΔTEC is important to characterize the scintillation. The random nature of the induced phase is guaranteed by assigning independent random variables to αi, i = 1, …, M. Henceforth, for the sake of simplicity, a uniform distribution will be used with 0 ≤ αi ≤ 1, although other distributions can be assumed. Applying the model from Equation (11), a value of ΔTEC = 1.17 TECU will provoke a maximum phase change of 2π rad in an L1-band carrier.
We consider that the plasma tubes in Equation (10) are separated by a constant distance D, with D ≥ 2L. Figure 3 sketches an example of the phase ϕsc (x, t) for the fixed time t. According to the figure, the HFI for a finite train of M electron density irregularities may be determined from:
12
where xi(t) = x1(0) + vdriftt + (i–1)D, with i = 1, …, M, indicates the central position of the density irregularity i at time t.
Example of the induced phase versus the phase-screen diffraction position x
The integrals in Equation (12) are of two types: A and B. Type A is defined in Equation (7) and is computed using Equation (8) or, approximately, using expression (A2) of the Appendix. Type B is defined by:
13
and can be expressed in terms of IA as:
14
Thus, the integral IHF(t), with the plasma profile defined by Equation (10), can be written in a compact way as:
15
We define next the (normalized) ionospheric scintillation process as:
16
such that, in the absence of scintillation, Ψ(t) = 1 (see Equation [9]). The scintillation process is characterized by the amplitude:
17
and phase:
18
The scintillation index S4 is often used to characterize the strength of the scintillation activity (Rino, 2011):
19
where the 〈.〉 brackets indicate time average (with typical intervals of 60 s) and I(t) = A2(t) is the received signal intensity.
As an example, we consider a train of plasma electron density irregularities with vdrift = 100 m / s, L = D / 2 = 250 m, and z = 350 km. Figure 4 displays the normalized intensity and phase modulo-2π for ΔTEC = 2 TECU, which corresponds to a strong scintillation scenario (S4 ≈ 0.91).
Example of the normalized signal intensity for an infinite train of electron density irregularities with ΔTEC = 2 TECU and corresponding phase modulo-2π time sequence
Figure 5 displays the dependency of the scintillation index S4 versus ΔTEC for different electron density irregularity sizes with the L1 and L2 bands. The results were obtained using an observation interval of 300 s with vdrift = 100 m / s. Although the statistical distributions of amplitude are independent of the TEC(x) profile drift velocity, the variability of the fading process with time (e.g., autocorrelation function and power spectrum) depends on vdrift. The plots show that S4 grows almost linearly for small values of ΔTEC with the dependency becoming nonlinear for larger values. In the linear region, the rate S4 / ΔTEC decreases when the electron density irregularity size grows.
Scintillation indices S4 versus ΔTEC in TEC units for the L1 and L2 bands
One of the advantages of the proposed Monte Carlo simulation technique is the ability to generate simultaneously fading processes for different bands using a common ionospheric scintillation model with a plausible physical justification. Figure 6 illustrates this scenario for the GPS L1 (f0 = 1.5754 GHz) and L2 (f0 = 1.2276 GHz) bands with ΔTEC = 1.5 TECU, vdrift = 100 m/s, and L = D / 2 = 250 m. Notice that the deep fades tend to occur independently for signals transmitted in different bands in agreement with experimental data (Jiao et al., 2016; Seo et al., 2011).
Fading intensities generated by a common scintillation process for different carrier frequencies
3 CHARACTERIZATION OF THE SCINTILLATION PROCESS
3.1 Amplitude
The amplitude A(t) of the scintillation process is difficult to characterize analytically because it changes with ΔTEC in an unpredictable manner, as illustrated in the histograms of Figure 7. The simulations have shown that no simple amplitude distribution law can be applied for weak and average scintillation. In those cases, the histograms suggest that the amplitude distribution may be well described by a finite mixture of distributions with the probability density function taking the form (Mclachlan & Peel, 2000):
20
where the mixing proportions (or weights) πi are nonnegative and sum to one and where fi(a) are the component densities. These densities may belong to the same parametric family (for example Gaussian) or to different families. The selection of the appropriate densities and their parameters is, in general, a challenging task.
Normalized probability histograms of the amplitude for L1 band with L = D /2 = 250 m and vdrift = 100 m/s
For ΔTEC ≥ 5 TECU, it was found that the amplitude was well approximated by the generalized gamma (or alpha-mu) distribution. The probability density function (pdf) of the generalized gamma distribution is given by (Yacoub, 2007):
21
where Γ(·) is the gamma function. This distribution includes as a particular case the Nakagami-m distribution for α = 2. The n-th moment is determined by E{An} = ξn/2Γ(μ + n/α)/Γ(μ). Assuming normalization of the scintillation intensity, E{A2} = 1, then ξ = Γ(μ)/Γ(μ + 2/α) and the expression of the pdf in Equation (21) can be written solely in terms of parameters a and μ. The scintillation index is given by:
22
This formula shows that, for each value of S4, there is an infinite number of pairs (α, μ). The optimal pair may be selected by minimizing, for instance, the mean quadratic error between the pdf curve and the pdf histogram. Figure 7(d) exhibits the histogram and the probability density function (dash-dot curve) assuming a generalized gamma distribution with parameters α = 1.7 and μ = 1.04. Notice the good approximation to the histogram provided by the curve.
3.2 Phase
Unlike signal intensity, phase cannot be measured directly. A modulo-2π version of the phase is obtained from the four-quadrant arctangent function as atan2(Im{Ψ(t)}, Re{Ψ (t)}). However, if the phase excursions were to exceed the cycle, the arctangent would exhibit 2π-jumps which do not make sense physically. Hence, phase unwrapping is necessary to correct any phase changes over π radians between two adjacent phase samples (Jiao et al., 2018). To calculate the unwrapped phase ϕ(t), we resort to a phase unwrapping algorithm such as the one described next (see also the discussion in Rino et al. [2020]).
For k = 1, 2, … define the modulo-2π phase θk = atan2(Im{Ψk}, Re{Ψk}):
Step 1: Do k = 1 and ϕ1 = θ1.
Step 2: Increment k (k = k +1).
Step 3: Do
(where is the largest integer that does not exceed x).Step 4: Let Δϕk = ϕk – ϕk–1. If |Δϕ| < π, go to Step 2; else go to Step 5.
Step 5: If Δϕk ≤ 0, then do ϕk = ϕk + 2π; else do ϕk = ϕk – 2π. Go to Step 2.
This phase unwrapper eliminates the 2π phase jumps provided that the increments |ϕk+1 – ϕk| between consecutive phase samples are smaller than π
The severity of the scintillation is traditionally quantified by the indices S4 (amplitude scintillation index) and σϕ (phase scintillation index) which indicate the average intensity of the signal variations over the preceding minutes. However, they do not offer any insight into the correlation of the instantaneous phase and amplitude of scintillation observed on multiple frequencies (Vilà-Valls et al., 2020).
Figures 8 and 9, with ΔTEC = 2 TECU, illustrate the behavior of the fading intensity and unwrapped phase in the case of strong scintillation, with the polar plot of Figure 9 highlighting the relation between the scintillation intensity (square of the complex phasor magnitude) and the phase (polar angle of the phasor; Carrano et al., 2013). The interval A-B corresponds to a deep fade (herein defined as I (t) < −10 dB) with the threshold for the deep fade indicated by the dash-dotted line in Figure 8 and by the ellipse in Figure 9. During the A-B path, the unwrapped phase ϕ(t) experiences a variation of about π radians as the curve of the polar plot passes from the first to the third quadrant (canonical fade). Note that the large phase transitions are present even in the absence of thermal noise, thus suggesting they are a characteristic of the radio-frequency (RF) environment (due to diffraction) and not an effect generated by the receiver (Breitsch et al., 2020; Carrano et al., 2013). The absolute value of the phase rate is maximum during a deep fade and leads often to the loss of carrier tracking by the receiver, which may be temporary (cycle slips) or may correspond to the definitive loss of carrier phase. In that case, the receiver has to switch to the acquisition mode. Besides the existence of large values of the phase rate during a deep fade, the loss of performance of the receiver’s phase loop is aggravated by the low instantaneous carrier-to-noise ratio.
Excerpt of fading intensity I(t) and unwrapped phase ϕ(t) for ΔTEC = 2 TECU (S4 ≈ 0.89); the interval A-B corresponds to a canonical fade.
Polar plot for the scintillation process of Figure 8
The relation between the scintillation indices σϕ and S4 is displayed in Figure 10(a) for different values of ΔTEC between 0.05 TECU and 2 TECU, using the L1 band, L = D/2 = 250 m and vdrift = 100 m/s. For each value of ΔTEC, 10 independent pairs (S4, σϕ) were computed. Notice the existence of two distinct regions in the plot. For S4 ≤ 0.6 there is an approximately linear relation between parameters with σϕ ≈ (3/4)S4 (indicated by the dashed line); this corresponds to the region where the phase distribution is monomodal. When S4 > 0.6, there is a strong dispersion of σϕ which is linked to the multimodal nature of the phase. A similar behavior is observed in plots generated with the Global Ionospheric Scintillation Model (GISM), which is recommended by ITU to predict the intensity of the ionospheric scintillation between Earth and space (ITU, 2019). Figure 10(b) displays the plot of σϕ versus S4, produced by the GISM for a receiver location in Hanoi, Vietnam, on April 12, 2013, using all the visible GPS satellites between 16 h and 19 h UTC with a time step of 1 minute.
Relation between parameters σϕ and S4: (a) using the proposed method with 0.05 TECU ≤ ΔTEC ≤ 2 TECU; (b) using the GISM
3.3 Power Spectra
The spectral density function (SDF) of both phase and intensity scintillation follows an inverse power-law distribution of the form (Ghafoori, 2012; Knight & Finn, 1998):
23
where Φ0 is the spectral strength at 1 Hz (assuming f0 ≪ 1 Hz), f is the frequency of the fluctuations, f0 is a frequency corresponding to the maximum irregularity size, v is the power spectral index, and p = 2ν is the power spectrum slope.
For f ≫ f0, the SDF can be well approximated by Φ(f) ≈Φ0 f−p. The SDF of the intensity and phase scintillations follow similar power law relationships for high-fluctuation frequencies, but the intensity spectrum is heavily attenuated below a certain cutoff frequency. In the case of weak scintillation, the cutoff frequency of the intensity scintillation SDF (Fresnel cutoff frequency) is given by 
(Knight & Finn, 1998). For instance, with vdrift = 50 m/s, this frequency is approximately 0.13 Hz.
Figure 11 displays an example of intensity and phase power spectra for the L1 band obtained with vdrift = 50 m/s and sampling frequency fs = 50 Hz. The slope of the intensity and phase spectra is p ≈ 4.0 for weak scintillation (ΔTEC = 0.25 TECU) in the decaying region of the plots. For strong scintillation (ΔTEC = 5 TECU), the phase spectrum is shifted upward (meaning larger phase fluctuations), the slope decreases (larger high-frequency content), and the intensity spectrum becomes wider. The same effect was reported, for instance, in Humphreys et al. (2009).
Example of intensity and phase power spectra for the L1 band
Figure 12 illustrates the dependence of the intensity and phase spectra on the drift velocity of the electron density depletion regions. As expected, the curves experimented a frequency rightward shift when vdrift increased, which corresponds to a growth of the high-frequency content. The shift appears to have been approximately linear with the velocity increment.
Example of the dependence of the intensity and phase power spectra on vdrift for the L1 band
4 FADING STATISTICS
Figure 13 illustrates the main parameters that characterize the intensity and phase fading processes. For the intensity, we consider the fading duration, fading depth (in dB), and time separation between adjacent fades. Following Jiao et al. (2016), we define a deep fade event when the normalized signal intensity drops below −10 dB. For the phase, we define the absolute value of the peak-to-peak phase change during a fade and the maximum absolute phase rate during fading. These parameters are also considered, for instance, in Jiao et al. (2016, 2018). As pointed out in Jiao et al. (2015), the measurements of C / N0 are not suitable to estimate the fading characteristics due to the ionospheric scintillation. In fact, since the measurements of C / N0 are based on time averages to minimize the effect of thermal noise, they tend to underestimate the fading level and to overestimate the fading duration. Besides these parameters, it is also important to determine the correlation of the fade events in the L1 (f0 = 1.57542 GHz), L2 (f0 = 1.2276 GHz), and L5 (f0 = 1.17645 GHz) bands.
Definition of the intensity and phase fading parameters
Fading duration determines whether reacquisition is needed in case the receiver has lost lock of the signal as it would affect the accuracy of the positioning solutions if tracking is maintained (Jiao et al., 2015). The time separation between fades has an impact on signal reacquisition and on the carrier smoothing process of code measurements. The time separation between fades of different bands determines whether and how frequency diversity techniques can be employed to minimize the receiver performance degradation due to scintillation (Jiao et al., 2016). Figures 14, 15, and 16 display the histograms of fading duration, fading depth, separation between fades, absolute peak-to-peak phase change, and maximum phase rate for the L1 band, corresponding to a time interval of 60 minutes, using a sampling rate equal to 50 Hz. It was assumed that ΔTEC = 2.5 TECU, vdrift = 50 m/s, and L = D / 2 = 250 m.
Normalized probability histograms of (a) fading duration and (b) fading depth for the L1 band
Normalized probability histogram of separation between fades for the L1 band
Normalized probability histograms of (a) peak-to-peak phase change and (b) maximum phase change rate for the L1 band
The dashed-dotted curve in Figure 14(b) refers to the fitting exponential distribution with a pdf of f (x; μ) = exp[−(−x −10] / μ] / μ, where x is expressed in dB and μ = 8. The curve in Figure 16(b) refers to the fitting log-logistic distribution with pdf f(x; μ, σ) = exp(z)/[σx (1 + exp(z))2], where z = (ln x − μ) / σ. The parameters μ = 2.6 and σ = 0.37 were used. Both fittings show good agreement with those reported in Jiao et al. (2018).
Table 1 lists the values of means and standard deviations of the different intensity and phase fading parameters during the above-mentioned time interval for the L1, L2, and L5 bands. The scenario ΔTEC = 2.5 TECU, vdrift = 50 m/s, and L = D / 2 = 250 m was maintained for the three bands. It is noteworthy that the number of deep fades increased (the time separation between fades diminished) when the carrier frequency decreased. Except for the maximum phase rate, the other parameters presented similar values to those obtained in Jiao et al. (2018).
Means and Standard Deviations of Intensity and Phase Parameters
Table 2 presents the probability of deep fades (I(t) < −10 dB) for different configurations of the L1, L2, and L5 bands. As expected, these probabilities increased slightly when the carrier frequency was diminished. The probabilities of simultaneous deep fades for two and three bands are also indicated. Probabilities smaller than 0.02 were obtained for the pairs (L1, L2) and (L1, L5). For the three bands, the probability was smaller than 0.01. These results agree with the behavior exhibited by the experimental data analyzed in Jiao et al. (2015).
Probabilities of Deep Fades for Different Band Configurations
The results show that the deep fades in two largely frequency-separated bands tend to occur almost independently (note that this is not true for the pair [L2, L5]). Simultaneously using signals in the L1 and L2 (or L5) bands provides an efficient way to minimize the probability of loss of lock in GNSS receivers due to strong scintillation. That is, it is possible to utilize the tracking results of other bands to assist the tracking of the deep fading band during scintillation (Jiao et al., 2015).
The relations between the peak-to-peak phase change and the fading depth and between the maximum phase rate and the fading depth for the L1 band are displayed, respectively, in Figure 17(a) and Figure 17(b). The peak-to-peak phase change increased with the fading depth, although in a nonlinear form, the dispersion of values grew with the fading depth. The plot of the log10 of the maximum phase rate exhibits a highly correlated linear relationship. The linear fitting is depicted as a solid straight line defined by y = −0.0395x + 0.4276. The corresponding correlation coefficient is high in absolute value (0.82). For the sake of brevity, only the results relative to the L1 band have been presented. However, simulations carried out with the remaining bands indicate similar results.
Relations between (a) peak-to-peak phase change and fading depth and (b) between maximum phase rate and fading depth for the L1 band
In general, the results shown in Figures 14 through 17 reveal a good accordance with those presented, for instance, in Jiao et al. (2018; obtained with real data), leading us to believe that the proposed simulator is able to represent with realism the effects of ionospheric scintillation on GNSS received signals in low latitude scenarios.
In Figure 18(a), we display the intensity of an L1-band signal received by the ionospheric scintillation monitoring station at the Hanoi University of Science and Technology on April 12, 2013 (16h 00m and 03s UTC) for an observation of two minutes (left-side plot; data available on the JRC Scintillation Repository; for details, consult Curran et al. (2014, 2015). In this interval, the measured value of S4 is approximately equal to 1. In Figure 18(b), we exhibit a synthesized intensity time series using the proposed algorithm with the design parameters adjusted in order to achieve approximately the same S4 index, mean duration, and spacing between fades as in the real signal. The adjustment led to the following design parameters: ΔTEC = 3 TECU, vdrift = 15 m/s, L = 150 m, and D = 300 m.
(a) Real and (b) synthesized intensity time series with common S4 value and average deep fade duration and spacing
Recall that a significant advantage of the synthesized intensity time series over the experimental one is that a time series with constant scintillation statistical characteristics and arbitrarily long duration can be generated in contrast with real signals whose characteristics are typically non-stationary over intervals larger than a few minutes.
The flowchart of Figure 19 presents the operations to be carried out to synthesize ionospheric scintillation time series using a set of pre-defined (target) scintillation parameters S4, 
, and 
, with 
and 
indicating, respectively, the average deep fade duration and the average separation between adjacent deep fades. The design parameters to be adjusted in the proposed algorithm are: ΔTEC, vdrift, D, and L. For instance, the values of 
and 
may be adjusted by changing the TEC(x) profile velocity vdrift: the larger the velocity is, the smaller the values of 
and 
will be. This procedure was used to generate the synthesized time series of Figure 18. Note, however, that certain combinations of scintillation parameters may not be reachable using any set of design parameters. In general, large values of S4 (> 1.4) cannot be obtained with the single phase-screen model adopted by the proposed algorithm (see Figure 5). However, the admissible range of values (0 ≤ S4 ≤ 1.4) encompasses the scintillation indices typically found in experimental scenarios (ITU, 2019).
Generation of ionospheric scintillation time-series using pre-defined (target) scintillation parameters
5 CONCLUSION
The paper presents a new fast Monte Carlo technique to simulate the effect of equatorial ionospheric scintillation on signals transmitted by GNSS satellites, useful for assessing the performance of the receivers in scintillation scenarios. The method uses a single-layer phase-screen model of the ionosphere by which we assume a train of plasma TEC(x) profiles drifting eastward with constant velocity. The total electron content for each electron density irregularity was chosen randomly using, for instance, a uniform distribution generator. The shape of the electron density irregularities was made triangular in the x-direction, allowing us to analytically determine the ionospheric scintillation random process. This process is expressed by the Huygens-Fresnel complex integral, which can be computed by a sum of Fresnel cosine and sine integrals, thus avoiding the use of lengthy numerical integrations and making the computation fast.
We calculated the distributions for the amplitude and phase of the disturbed signals versus the scintillation index S4. It was shown that, in general, the amplitude of the fading process is difficult to characterize in terms of known distributions (it was suggested that it could be characterized as a finite mixture of distributions). However, for strong scintillation scenarios, it was verified that the amplitude is well approximated by the generalized gamma (or alpha-mu) distribution, as happens in real scenarios. The argument of the Huygens-Fresnel integral was utilized to generate unwrapped phase sequences that are useful to evaluate the performance of phase-locked loops in GNSS receivers. Several statistics were calculated for signals generated with the proposed technique, including fading duration and depth, separation between fades, peak-to-peak phase change, and maximum phase change rate. The results obtained for the L1, L2, and L5 bands show good agreement with the results provided by experimental data.
HOW TO CITE THIS ARTICLE
Nunes, F., Sousa, F., & Marcal, J. (2022) Multi-frequency simulation of ionospheric scintillation using a phase-screen model. NAVIGATION, 69(4). https://doi.org/10.33012/navi.545
ACKNOWLEDGMENTS
This work was partially supported by the Portuguese national funding agency FCT under project UIDB/EEA/50008/2020.
APPENDIX
Fast computation of the complex Fresnel integral Ξ(y) = C(y) + jS(y) can be obtained by resorting to the approximations presented in Chapter 7 of Abramowitz and Stegun (1972) for the sine and cosine Fresnel integrals. It can be shown that:
A1
where sign(y) = −1, if y < 0 and sign(y) = +1, otherwise. The auxiliary functions, f (x) and g(x), are given by the rational approximations:
A2
A3
The errors of these functions do not exceed 2 × 10−3. However, for more accurate approximations, the Boersma algorithm (Boersma, 1960) may be employed. Replacing the auxiliary functions in Equation (A1) of the Appendix and taking into account Equation (8) leads to:
A4
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.