© Owned by the authors, Published by EDP Sciences 2011

1. Introduction

As a result of propagation through ionosphere electron density irregularities, transionospheric radio signals may experience amplitude and phase fluctuations. In equatorial regions, these signal fluctuations specially occur during equinoxes, after sunset and last for a few hours. They are more intense in periods of high solar activity. There is also a longitudinal dependency. Scintillations are more common in South America near the December solstice than at the equinoxes. These fluctuations result in signal degradation from VHF up to C band. They are a major issue for many systems including Global Navigation Satellite Systems (GNSS), telecommunications, remote sensing and earth observation systems.

The signal fluctuations, referred to as scintillations, are created by random fluctuations of the medium’s refractive index, which are caused by inhomogeneities inside the ionosphere. These inhomogeneities are substructures of bubbles, which may reach dimensions of several hundreds of kilometers as can be seen from radar observations (Costa et al. 2011). These bubbles present a patchy structure. They appear after sunset, when the sun ionization drops to zero. Instability processes develop inside these bubbles with creation of turbulences inside the medium. As a result, depletions of electron density appear. In the L band and for the distances usually considered, the diffracting pattern of inhomogeneities in the range of 1 km size is inside the first Fresnel zone and contributes to scintillation (Wernik et al. 1980).

Ionosphere scintillation is currently the object of many measurement campaigns, both at low and high latitudes. The main campaigns related to the areas covered by experimentation are the Low Latitude Ionospheric Sensor network in South America (Valladares 2009), the SCINDA network in America and Africa (Carrano 2006) and the CHAIN network in Canada (Jayachandran 2009). But there are many others worldwide. In this paper, we will refer to the PRIS measurement campaign conducted in the frame of an ESA/ESTEC contract during years 2005–2006, with measurements at low and high latitudes (Béniguel 2009).

A propagation model aimed at reproducing the signal fluctuations is presented in this paper. The calculation addresses the evaluation of the transmitted field and of its second-order moment. In most cases, as in GNSS applications, the knowledge of the transmitted field allows estimating the degradation of performances, due to both intensity and phase fluctuations. For radar observations, the knowledge of the second-order moment is also required. The signal coherence properties of the medium, both for time and frequency, are important in this case to assess the radar performances (Knepp 1989).

2. Scattered field calculation

2.1. Introduction

The model presented in this paper (Global Ionospheric Scintillation Propagation Model, GISM) uses the Multiple Phase Screen (MPS) technique (Knepp 1983; Béniguel 2002, 2004; Gherm 2005). The locations of transmitter and receiver are arbitrary. The incidence link angle is arbitrary regarding the ionosphere layers and the magnetic field vector orientation. It can cross the entire ionosphere or a small part of it. At each screen location along the line of sight, the parabolic equation (PE) is solved. GISM allows calculating mean errors and scintillations due to propagation through the ionosphere.

The mean errors are obtained using a ray technique solving the Haselgrove equations (Budden 1985). The ionosphere electron density at any point inside the medium, required for this calculation, is provided by the NeQuick model (Radicella 2009), which is included in the GISM.

The line of sight being determined, the fluctuations are calculated in a second stage using the multiple phase screen technique. The medium is divided into successive layers, each of them acting as a phase screen. In this technique, which is detailed hereafter, the field is scattered from one screen to the next one.

2.2. Theoretical formulation

The wave propagation is calculated solving the Helmholtz equation (Ishimaru 1978).(1)where

• is the local wave number.

• µ₀, ε₀ and k₀ are the free space permeability, permittivity and wave number.

• is one observation point inside the medium and z is the coordinate along the direction of propagation.

The dielectric permittivity along the main propagation axis z is written:(2)with being the random part of the relative dielectric permittivity.

Introducing the complex amplitude of the stochastic field(3)and assuming that the variation of the complex amplitude is mainly in the direction perpendicular to the main propagation axis (parabolic approximation), the stochastic PE for the complex amplitude can be written in the form(4)where is the transverse Laplacian.

2.3. Algorithm

To solve this equation, the medium is divided into series of successive layers (or screens) perpendicular to the main propagation axis, each one being characterized by local homogeneous statistical properties. The solution is then obtained by iterating successively scattering and propagation calculations as detailed hereafter.

The parabolic wave equation is split into two equations. The first one describes the phase change due to the presence of random fluctuations ε₁(r, z); r is the distance to the propagation direction main axis.(5)with solution(6)

The second equation describes propagation between two screens(7)

This equation is solved in the transform domain with solution(8)

The Fourier transform of the complex amplitude is calculated after the first step. It is as an input to the second step. Applying this two-step technique to each successive layer, the MPS solution of the PE is obtained (Yeh 1977; Knepp 1983). All these calculations can be performed using FFT techniques.

In most of the cases considered, the source point is very far away from the fluctuating medium. The incident field on the first layer is a plane wave and the initial value of the field complex amplitude on this screen is 1.

2.4. Phase synthesis on a phase screen

In the MPS technique, successive planes perpendicular to the direction of propagation are considered. On each one of these planes, a phase synthesis shall be performed. The field is then diffracted from one plane to the next one. The necessity to use 2D or 1D phase screens is addressed in this section.

In general, the medium’s Power Spectral Density (PSD) of phase fluctuations can be approximated by expression:(9)where

• C_p characterizes the turbulence strength. It is related to the variance of the electron density.

• q₀ = 2 < π/L₀ and L₀ is the outer scale of the inhomogeneities.

• p is the spectrum slope. The spectrum is consequently linear using Log-Log scales.

• The variables k and r are corresponding variables in the Fourier transform.

If we consider that the medium is a 2D isotropic medium, the phase autocorrelation function can be calculated using the integral given below (10) with corresponding result on the right-hand side.(10)with K the modified Bessel function.

If instead, 1D phase screens are used, the integral and corresponding solution is:(11)

The only difference is that the slope p shall be decreased by 1. One dimension phase screens can consequently be used with this slight modification with a significant simplification in the algorithm.

Models to reproduce the inhomogeneities development inside ionosphere can be built solving the plasma equations, namely the continuity and momentum equations (Béniguel 2011). Such models show that the inhomogeneities grow in a plane perpendicular to the terrestrial magnetic field and appear finally as tubes aligned with it. The dimensions of these tubes may reach hundreds of kilometers. The medium is consequently anisotropic (Rino 1979).

The PSD of phase fluctuations is given in this case by the more general expression:(12)where the coefficients a and b are the axial ratios of the irregularities and A, B and C depend on the orientation of the wave vector with respect to the irregularities principal axis. The wave number is a vector with components K_x⊥ and K_y⊥ in a plane perpendicular to the propagation direction.

Using an appropriate change of variables (Rino & Fremouw 1977; Rogers 2009), it can be shown that (11) is equivalent to (8) with an additional geometrical factor. As a consequence all the analysis can be done using a 2D geometry (1D phase screens). The axis reference system contains the direction of propagation (the line of sight) and the terrestrial magnetic field vector.

In the algorithm which has been implemented, at each phase screen, the phase synthesis is done taking the inverse Fourier transform of the product of a random numbers series with given uniform probability density by numbers with the required spectral density. This provides a series of random numbers with the required statistical properties (Picinbono 1968).

2.5. Results obtained

Inputs of the model are the transmitter and receiver locations, the time, day and year of observation, and the geophysical parameters. Based on the PRIS measurement campaign, experimental laws have been derived for the geographic and local time dependency. As mentioned previously the spectrum is characterized by three parameters: the slope, a typical dimension of inhomogeneities and the strength. Default values are respectively set to p = 3, L₀ = 1 km and. As geophysical parameters, only the average 10.7 cm solar flux number is considered. Its value is taken from the curves published by the National Oceanic and Atmospheric Administration (www.noaa.gov). The magnetic activity is ignored. It is not considered either in the NeQuick model.

For earth space links, the source point is the antenna on board the satellite and the observation point is on the ground. Once the line of sight is determined, phase screens are set along this line and statistical parameters are associated to each one of these screens. The algorithm provides the field U(r) at the observation point, where r is the dimension transverse to the direction of propagation. Dependency of the received field on the time can be obtained, considering both the source displacement and the medium drift velocity. A time series of the received field is obtained in that case.

GISM calculates a mean value of the scintillation indices both for intensity and phase fluctuations. The fluctuating medium is assumed to be statistically homogeneous. This could be improved however including a dependency of the statistical parameters on the altitude. The reality is different to some extent. The medium has a patchy structure and links meeting the geographic and time conditions may not be affected due to this patchy structure. Consequently a probability of occurrence should be given together with the mean value. This is not provided in the current version of the model. The corresponding probability shall be obtained from measurement analysis.

Two indices are defined to characterize the scintillations: the standard deviation of the normalized intensity, named S4, and the phase standard deviation. The scintillation event strength is defined with respect to the S4 value which is between 0 and 1. A value of 1 will correspond to about 35 dB peak to peak of intensity fluctuations. One example of the time series provided by the model is represented in Figure 1 for a strong fluctuation case. The scintillation indices are calculated from these 500 samples (intensity and phase). The scintillation strength is weak (S4 < 0.3), medium (0.3 < S4 < 0.6) or strong (S4 > 0.6) depending on the case. This usual classification refers to the fade levels and the resulting constraints on a navigation system, from −2 dB to +2 dB in the weak regime to more than 20 dB peak to peak for the strong regime.

Fig. 1. Intensity and phase time series / strong fluctuation event (S4 = 0.9).

In addition to the scintillation indices, the time series analysis enables one to estimate the intensity and phase probabilities, the fades statistics and the spectrum characteristics.

Figure 2 shows the correspondence between a Total Electron Content (TEC) map and a scintillation map. Those two maps were obtained by modeling using the NeQuick (Radicella 2009) model for the TEC and GISM for scintillations. They correspond to vertical links. The electron density is consequently integrated along a vertical line at each grid point on the map to get the TEC value. Slant observations may exhibit higher values. The propagation length inside the ionosphere would increase in that case, and by consequence the level of scintillation obtained would also increase.

Fig. 2. TEC (left panel) and scintillation map (right panel) obtained by modeling.

Figure 2 was obtained with the default spectral parameters as defined at the beginning of this section, including the electron density variance, and an average solar radio flux at 10.7 cm set to 150. This corresponds to a high value. Universal time is 12:00 p.m. At this time the peak values for the TEC occur in the Pacific Ocean area. For the scintillations, the time duration of the events is a few hours after sunset. Both plots exhibit peak values on each side of the magnetic equator near the crests of the equatorial anomaly. The values decrease with increasing latitude. For scintillations, the model calculates the effects at the equatorial regions. The high latitudes region, also affected by the scintillation, has not been considered in this example. The TEC maximum is 80 TEC units, which is a significant value. It is directly linked to the solar flux value. The peak value for the intensity RMS (S4 parameter) equal to 0.7 corresponds to strong fluctuations.

2.6. Comparison with measurements

The results reported hereafter are taken from the PRIS measurement campaign (Béniguel 2009) carried out under ESA/ESTEC contract N° 19530. For this study, a number of receivers were deployed both at low and high latitudes, in particular in Vietnam, Indonesia, Guiana, Cameroon, Chad and Sweden. These receivers were dedicated receivers, operating at 50 Hz. A data bank was constituted and the scintillation characteristics were derived from an extensive analysis of this data bank. Comparisons between measurements and results provided by the GISM model in the same conditions were performed both for the scintillation indices and on the spectrum.

2.6.1. Scintillation indices

One week of measurements at Cayenne, French Guiana was selected. The results are presented in Figure 3. The x axis corresponds to local time at receiver location. The 0 value has been set arbitrarily to saturday 19:00 PST, the week of observation.

Fig. 3. Intensity and phase scintillation indices measurements at Cayenne during GPS week N° 377. The PRN numbers refer to the 14 GPS satellites used in this figure.

2.6.2. Measurements

The local time of the x axis corresponds to hours in GPS time. Each point corresponds to a 1-min sample. Only points with a S4 value greater than 0.2 were retained in the analysis. A 5° mask elevation angle was taken when recording the data. In addition, multipath is rejected using the code carrier divergence algorithm recommended in the Novatel GSV 4004 user manual. As can be seen in Figure 3, the points are clustered every evening at post-sunset hours, typically 19:00–24:00. No average is taken on the data. The scintillation activity occurred quite regularly that week with comparable levels. The S4 average value is about 0.4. The flux number that week (GPS week N° 377, modulo 1024) was equal to 90.

The phase fluctuations are plotted concurrently. The mean value is about 0.2, consequently lower than the S4 value. This observation is quite common. A few points exhibit high values. Deep fades occur concurrently to these high values. In the case of very small values this creates phase jumps. As a consequence the phase and intensity standard deviations are no longer related.

The scintillation characteristics, both indices and spectral parameters, have been calculated using 1-min samples. This calculation brings no particular difficulty for the intensity, which practically does not change during one minute. The calculation of the phase parameters is more difficult. A high-pass 6th-order filter is used to remove the low-frequency components of the signal, due to the satellite motion on its trajectory.

2.6.3. Modeling

The scenario was replayed using the corresponding Yuma files for one particular day of the week (cf. Fig. 4). A different day will not bring significant differences considering that the geophysical parameters would have been quite identical. The flux number, input to GISM, has been set to 90. As mentioned previously, the model provides a mean value. It overestimates the number of affected links due to the fact that the probability of occurrence is not considered. Only the mean values can be compared. The scintillation intensity index mean value is about 0.4, corresponding to the measurements. The scintillation phase index mean value is slightly greater than the one recorded in the measurements. In both cases the phase RMS is lower than the intensity RMS, and in both cases some points exhibit high values due to the phase jumps.

Fig. 4. Intensity and phase scintillation indices on day 314, GPS week N° 377, obtained by modelling. Different GPS satellites are identified by different numbers, symbols and colours.

2.6.4. Spectrum of received signal

The measurement spectral parameters have been obtained using a periodogram analysis. One typical result for a medium value of the indices (S4 = 0.5) is presented in the right panel of Figure 5. The low-frequency level (below 0.1 Hz) is arbitrary and is meaningless. In both cases the meaningful part of the spectrum is between 0.1 and 1 Hz. This frequency window has been selected to calculate the slope and strength of scintillations. The phase signal is significantly affected by noise in that case, the signal being corrupted by the receiver noise.

Fig. 5. Spectrum comparisons: typical result on 1 sample: modeling (left panel) and measurements (right panel).

The left panel of Figure 5 shows the corresponding PSD obtained with GISM model for comparable scintillation indices values. The intensity slopes are in good agreement. This is not the case for the phase due to the noise level.

This study will be complemented in the future using a wavelet analysis to fit this window with more appropriate values if necessary. This might be the case in particular at high latitudes for which the spectral components are expected to be different.

The signal spectrum deduced from the measurements has been approximated by expression Tf^−P which applies to most of the cases, where T, defined as the value at 1 Hz, is related to the turbulence strength. The GPS week N° 377 (modulo 1024) was selected again to derive the parameters p and T. Samples corresponding to S4 values greater than 0.2 were only selected in this analysis in order to diminish the effect on the noise on the results. The slope value, plotted in Figure 6, has a median value equal to 2.7. This result is in agreement with what is usually considered in the literature (Wernik 2007). The slope value decreases with time after sunset corresponding to the fact that the inhomogeneities sizes decrease with time after sunset. It should be noticed however that the PRIS measurement campaign was done in a year close to solar minimum. High solar activity values might be different, in particular for the strength.

Fig. 6. Slope value deduced from measurements in Cayenne (1 week of 50 Hz raw data files analysis).

3. Second-order moment of the field

For a radar application, the coherence properties of the transmitted field are required. The mutual coherence function (MCF), noted Γ, characterizes the coherence properties of the transmitted field and its determination is required for a radar application.(13)

The Γ function is obtained starting from equation (4) written for two frequencies and two positions, then combining the different equations (Yeh & Liu 1977; Ishimaru 1978). The final equation is written below:(14)where

• B_ε (z, ρ) = 〈ε(r₁) ε(r₂)〉 is the autocorrelation of permittivity fluctuations

• ξ = σ_Ne / Ne and Ne is the electron density

• A_ξ (ρ) = ∫ B_ξ(z, ρ) dz

• are the Laplacians with respect to r₁ and r₂ and ρ = r₁ − r₂.

Equation (13) can be re-organized as written below:(15)This new expression is again a PE. It can be solved using the same technique as the one presented in Section 2. It is two dimensions with respect to distance and frequency separation. In the transform domain, it provides the medium scattering function dependency with respect to Doppler frequency and delay.

The algorithm is similar to the one used for equation (4), alternating scattering and propagation calculations. If a quadratic approximation of the phase structure is used, which can be demonstrated whatever the spectrum slope value is, most of the calculations can be performed analytically (Nickisch 1992; Knepp & Nickisch 2009).

Despite the fact that the Γ function depends on two variables, the Fourier transforms reduce to a 1D FFT, the second transformation being done analytically. The two results presented in Figure 7 have been obtained respectively in HF (left panel) and P band (right panel). A spread factor, named Q, related to the medium parameters can be defined and depending on its value, the shape of the scattering function may change significantly. The inhomogeneities sizes, the frequency and the distances have a strong influence on the signal spreading.

Fig. 7. Mutual coherence function; one single screen. HF propagation (left panel, large Q value) and P band propagation (right panel, Q value around 1).

Solving equation (14) also provides the space and frequency coherence of the medium. The space and time (assuming a displacement with velocity V) coherences are given by (Knepp 1989):(16)This coherence time may be a limitation for some applications, in particular for Space-Based Synthetic Aperture Radars, considering the satellite displacement velocity. On the contrary, the frequency coherence is quite large, and its decrease due to ionosphere turbulence appears to be less critical for most applications.

In the case of one single screen the whole calculation can be performed analytically. The solution is given by the expression below where S, B and P are functions of the phase variance, the frequency and the medium parameters as introduced by (Nickisch 1992)(17)with(18)and(19)z_i is the coordinate along the propagation direction.

This expression allows conducting a parametric analysis study.

4. Conclusion

The GISM model, presented in this paper, uses a classical phase screen technique algorithm. Submodels have been included to estimate some specific parameters and take the geophysical dependencies into account. This concerns in particular the local time and seasonal dependency, the spectrum parameters, the inhomogeneities dimensions and their correlation distance. They were derived from measurement campaign results.

When compared to measurements, the modeling results show a relatively good agreement as presented in Section 2.6. Some work still needs to be done with respect to the phase characterization and to the calculation of its spectral parameters. No results were presented for high latitudes. Data from this region will be collected in the framework of the Monitor campaign. Monitor is a new ESA promising measurement campaign, currently ongoing, with a higher level of requirements (Prieto Cerdeira 2011). The signals from stations mostly located at low and high latitudes will be received in quasi-real time. The measurements will be multi-frequency and will use in some places co-located receivers in order to derive the medium drift velocity and correlation distance. In addition, in case of high scintillations a bitgrabber will be activated for post-processing analysis. This will allow a better characterization of the signal for extreme events which are of particular interest for GNSS applications.

The last section is focused on the calculation of the transmitted field MCF. One analytical solution has been derived by assimilating the medium to one single phase screen. Comparisons for this point were made only with respect to published results. The characterization of this function is of particular interest for radar observations and remote sensing applications.

Acknowledgments

The authors thank the reviewers for many helpful remarks. This study was supported by ESA / ESTEC and by SOlar –TERrestrial Investigations and Archives (SOTERIA) project, under FP7 Grant Agreement N°. 218816.

References

All Figures

	Fig. 1. Intensity and phase time series / strong fluctuation event (S4 = 0.9).
In the text

	Fig. 2. TEC (left panel) and scintillation map (right panel) obtained by modeling.
In the text

	Fig. 3. Intensity and phase scintillation indices measurements at Cayenne during GPS week N° 377. The PRN numbers refer to the 14 GPS satellites used in this figure.
In the text

	Fig. 4. Intensity and phase scintillation indices on day 314, GPS week N° 377, obtained by modelling. Different GPS satellites are identified by different numbers, symbols and colours.
In the text

	Fig. 5. Spectrum comparisons: typical result on 1 sample: modeling (left panel) and measurements (right panel).
In the text

	Fig. 6. Slope value deduced from measurements in Cayenne (1 week of 50 Hz raw data files analysis).
In the text

	Fig. 7. Mutual coherence function; one single screen. HF propagation (left panel, large Q value) and P band propagation (right panel, Q value around 1).
In the text