Abstract
This paper presents a physics-based, multifrequency, strong scintillation simulator that requires only two input parameters: the S4 index and decorrelation time τ0. This simulator is developed based on the two-component power-law phase screen model, which is specified by five parameters and suitable for representation of strong equatorial scintillation. By setting three of the model parameters to representative values, numerical mappings from the user input parameter set (S4 and τ0) to the remaining model parameter subset of two parameters can be established through numerical evaluation. The numerical mappings are then used to drive the scintillation model to generate time series of scintillation amplitude and phase. A large group of strong scintillation data (with S4 > 0.6) from two equatorial sites are processed to obtain the representative values of the three model default parameters and to validate the numerical parameter mappings. A MATLAB implementation of the scintillation simulator has been made available.
1 INTRODUCTION
Ionospheric plasma irregularities along the propagation path of radio wave signals cause temporal fluctuations in signal phase and amplitude that are referred to as ionospheric scintillation.1 Ionospheric scintillation most frequently occurs in equatorial and high latitude areas, while the strongest effects are often observed in equatorial areas during the postsunset and premidnight period characterized with simultaneous deep amplitude fading and fast phase fluctuations.2-4 The scintillation effects can severely degrade receiver performance5,6 and impact a variety of GNSS applications, including remote sensing,7-10 high accuracy positioning,11-13 applications that rely on Satellite-Based Augmentation Systems (SBAS),14-16 etc. It is therefore of great interest to develop and evaluate tracking algorithms that can mitigate these adverse effects and improve robustness for GNSS receivers encountering ionospheric scintillation.
Although in recent years, there have been increased availability of real scintillation intermediate-frequency (IF) data collected globally,17,18 a scintillation simulator is still a necessary tool for the scientific, engineering, and GNSS application community to study the physical mechanism and effects of scintillation and to develop receivers that can mitigate these effects. For users focused on the latter objective, it is desirable that the simulator be not only capable of capturing realistic scintillation effects, but also convenient and intuitive to configure and use. The objective of this paper is to present an ionospheric scintillation simulator developed to address this issue. The simulator presented in this paper is a physics-based, multifrequency, strong scintillation simulator that only requires two user input parameters: the scintillation index S4 and the intensity decorrelation time τ0.
Scintillation models that have been studied and implemented for GNSS signal simulation include statistical models and physics-based models. Statistical models simplify the scintillation effects as a stochastic process parameterized by a small number of scintillation indicators. The corresponding probability distribution functions (PDF) and spectral density functions (SDF) are obtained from empirical scintillation data through statistical analysis and spectrum shaping.19-21 The statistical modeling approaches are straightforward and simple to configure. For example, the widely used Cornell scintillation model only requires the expected S4 index and complex scintillation signal decorrelation time as input.21 However, these models do not preserve the information on the ionospheric irregularity structures and propagation geometry that are responsible for causing scintillation effects observed in the empirical data.22 In addition, due to the lack of strong scintillation data, existing statistical models have not been able to simulate realistic strong scintillation signals. Moreover, they are unable to retain the correlated effects among multifrequency scintillation, and they cannot simulate scintillation effects experienced on dynamic platforms.
Physics-based scintillation models, on the other hand, are based on theory of radio wave propagation through the ionospheric plasma structures. The most commonly used are the power-law phase screen models that employ phase screens to abstract irregularity structures in a statistical sense and allow oblique signal propagation in an anisotropic medium. The phase screens are spectrally characterized with a power-law model that is typically implemented as one-component or two-component power laws. The one-component power-law model was implemented and demonstrated suitable for weak to moderate scintillation simulation in Deshpande et al23, Chartier et al24 and Ghafoori and Skone25 while strong equatorial scintillation effects are better captured by the two-component power-law phase screen model (TPPSM) as studied in Rino et al26,27 and Carrano and Rino.28 The specification of physics-based models previously required knowledge of a large number of parameters involved in the calculation of propagation geometry and the propagation signal’s ionospheric piercing point (IPP) scan velocity, which made them cumbersome to work with.
A compact TPPSM was presented in Rino et al22,29 where a space-to-time scaling parameter (ρF/veff) was introduced that can absorb all dependency on the propagation geometry and IPP scan velocity. This scaling parameter, along with a strength parameter (U) and three spectral parameters (p1,p2,μ0), forms the set of parameters that defines the scintillation model. Given a segment of real scintillation data, an irregularity parameter estimation (IPE) method developed in Carrano and Rino28 and Carrano et al30 extracts the values of the entire parameter set to generate statistically equivalent scintillation effects. No explicit specification of the propagation geometry is needed since the geometry is already embedded in the extracted value of ρF/veff. This model also enabled the study of scintillation effects on GNSS signals observed on dynamic platforms such as aircraft and low-Earth-orbit (LEO) satellites as presented in Jiao et al31 and Xu et al,32 respectively. In both references, {U, p1, p2, μ0} were extracted from real scintillation data segments using the IPE to specify the ionospheric structure realization, while ρF/veff was calculated for various user platform dynamics to investigate the platform dynamics’ impact on scintillation signal propagation.
The simulator in this paper is developed with the motivation to relate the TPPSM parameters to a smaller set of scintillation indicators on the receiver observation plane. Only two parameters will be used: the S4 index and the intensity decorrelation time τ0. Both parameters are widely used indicators that directly reflect the amplitude fluctuation severity and the temporal variability of the scintillation effects on the signals, respectively. They also capture the stress of a receiver tracking loop under scintillation conditions, which make this simulator more practical for receiver algorithm testing. In Carrano and Rino,28 {U,ρF/veff} were shown numerically to be more closely related to the temporal characteristics of the scintillation signals (ie, {S4,τ0}) than the three spectral parameters {p1, p2, μ0}. Therefore, in this paper, by defaulting the three spectral parameters to representative values during numerical evaluation of the model, numerical mappings are established between the model parameter subset {U,ρF/veff} and the user input parameter set {S4,τ0}.
To obtain the representative values of the spectral parameters and to validate the numerical mappings, we first apply the IPE method to a large group of strong scintillation data (with S4 > 0.6) from two equatorial sites (Ascension Island and Hong Kong) to establish the profiles of the five model parameters for a strong equatorial scintillation scenario. The threshold of 0.6 was heuristically selected for strong scintillation based on experience processing a large amount of scintillation data. The representative values for the spectral content parameters are then determined based on their respective profiles. Based on the estimates of {U,ρF/veff} and {S4,τ0} obtained from the real scintillation data, the numerical mappings are demonstrated to represent the most typical relationships between the two parameter sets observed in real data. The results also confirmed that the reduction in the input parameters does not undermine the TPPSM’s effectiveness in generating realistic scintillation effects and hence the simulator’s performance. In addition, we validated the IPE method’s performance in terms of triple-frequency consistency using the real scintillation data.
Based on these numerical mappings, the simulator can then be controlled by specifying the expected S4 and τ0 values for a ground-based stationary receiver to perform the TPPSM-based structure realization and wave propagation. Statistically equivalent realizations of scintillation effects are therefore generated. In addition, this simulator can generate realistic scintillation effects that are observed on platforms with user-defined dynamics under the ionospheric condition responsible for producing the ground received scintillation, as in the work presented in Jiao et al31 and Xu et al.32 A MATLAB implementation of this two-parameter simulator is made available for public access. It should be mentioned, the numerical mappings that relate the scintillation indicators to model parameters can also be established in the same fashion for weak and moderate scintillation, where the phase screen model reduces to a one-component power law. This paper is intended to focus on the strong scintillation scenario which is more threatening for the navigation community.
The rest of this paper is organized as follows. Section 2 presents the architecture of this two-parameter scintillation simulator. Section 3 summarizes the mathematics in the TPPSM-based phase screen realization and wave propagation. Section 4 describes the IPE method and the calculation of S4 and τ0. Section 5 describes the real scintillation data set used in this study. Section 6 discusses real data analysis, including the evaluation of IPE consistency across the three GPS carrier frequency bands, distributions of the model spectral parameters obtained from the data set, and the parameter mappings and validation using the real data. Section 7 concludes the new simulator model approach. Access to the MATLAB implementation of this simulator are given in the Appendix, accompanied with its configuration guide.
2 TWO-PARAMETER SCINTILLATION SIMULATOR ROUTINE DESCRIPTION
This section gives a description of the two-parameter scintillation simulator architecture. Figure 1 depicts the flow chart of the simulation routine. The parallelograms on the left side depict the user inputs, and the rectangular components show the different steps of simulation.
Flow chart of the scintillation simulator. The parallelograms on the left side depict the parameters to be specified by users, while the rectangular components are the different steps within the scintillation simulation. The simulator outputs are the scintillation amplitude (δA) and phase (δϕ) time series
As can be seen in Figure 1, the parameters to be specified by users include two categories: ground observed scintillation indicators and propagation geometric parameters. The static ground observed scintillation indicators are the S4 index and intensity decorrelation time τ0. The S4 index is the standard deviation of the normalized signal intensity. In this study, we shall focus on strong scintillation where S4>0.6. The intensity decorrelation time τ0 is defined as the time delay where the normalized autocorrelation function of the signal intensity decreases21 to e−1.
The user-input {S4,τ0} are first converted to 
in the parameter mapping step. The 
is the time scaling factor of the static ground observed scintillation and will be used in the propagation geometry calculation step. As mentioned earlier, the three spectral parameters {p1, p2, μ0} required to completely specify the phase screen realization are defaulted to their typical values, which will be determined in detail in Section 6.2. Based on the mapped value of U and the defaulted spectral parameters, a realization of phase screen is then generated to represent the ionospheric plasma structure responsible for ground observed scintillation effects. A brief description of the phase screen realization will be given in the next section.
For stationary platforms, the users only need to specify the two ground observed scintillation indicators, and the simulator will generate scintillation effects with S4 and τ0 values that match the user input. To simulate scintillation effects experienced on dynamic platforms, the users will need to specify the propagation geometric parameters to enable the calculation of ρF/veff for the user-defined dynamic platform. The propagation geometric parameters include the platform position and velocity (PV), date and time, and satellite PRN number. Based on the user specified PRN, date, and time, the satellite orbit is calculated using the corresponding ephemeris, which is automatically downloaded from Crustal Dynamics Data Information System (CDDIS) of NASA.
All the geometric and dynamic dependencies are encapsulated in the calculation of the scaling factor ρF/veff in the compact TPPSM.31 The quantities involved in this calculation include the angle between the signal propagation direction and the geomagnetic field vector, the satellite scan velocity and receiver scan velocity at the phase screen, and the drift velocity of the ionospheric irregularities (vdrift).33 Among these quantities, the only unknown is vdrift, while the others can be calculated with the user-input propagation geometric parameters (satellite orbit and receiver PV) and existing models (such as the IGRF model). In this simulator, in order to obtain the ρF/veff value for user-defined dynamic platform (ρF/veffdyn), the value of vdrift is first numerically solved using the ρF/veffstatic value from ground observed scintillation.
Finally, a plane wave is propagated through the phase screen realization, following the user specified geometry and dynamics embedded in ρF/veffdyn (or ρF/veffstatic for stationary-platform scintillation). This gives simulated scintillating wave fields at the receiver, ie, the scintillation induced phase (δϕ) and amplitude (δA) time series. A description of the wave propagation step will also be given in the next section. In this simulator, the sampling rate for the scintillation effects time series is 100 Hz, which adequately captures the scintillation effects spectrum and is consistent with the sampling rate of the receiver measurements used for numerous scintillation studies. The simulated scintillation amplitude and phase can be further interpolated and modulated onto nominal GNSS baseband samples to generate GNSS scintillation signals for receiver processing. Mathematical representations for the scintillation signal modulation can be found in Rino et al.29
3 PHASE SCREEN REALIZATION AND WAVE PROPAGATION
The propagation theory behind the scintillation generator depicted in Figure 1 was mainly developed in previous study33 and updated in Rino et al.22,29 This section summarizes the fundamental mathematics for the phase screen realization step and wave propagation step in this simulator. Readers are referred to Rino et al29,33 for more details.
In the equatorial region, ionospheric irregularities are typically highly elongated along the geomagnetic field lines.34,35 The phase screen structure is simplified to be two-dimensional, namely, the only variant crossing the field-aligned direction at a given height. Based on this approximation, the complex field ψ at a propagation distance x from the phase screen is generated by the following forward and inverse discrete Fourier transforms (DFTs)29:
1
2
where y is the direction perpendicular to the propagation plane, ϕ(y) is the path-integrated phase structure, q is the wavenumber along the y direction with Δq being its resolution, k is the signal carrier’s wavenumber, and N is the size of the DFTs.
The path-integrated phase structure is characterized by a one-dimensional, two-component power-law SDF as follows:
3
where Cp is the turbulence strength, p1 and p2 are the spectral indices, and q0 is the break wavenumber.
Simplification is achieved by scaling the wavenumber q with the Fresnel scale 
to a normalized unit such that μ = qρF. After such normalization, the SDF of the phase screen can be rewritten as
4
where 
. Thus, the universal scattering strength U can be defined as
5
U is essentially the normalized phase spectral power at the Fresnel scale as U ≡ P(μ = 1).28
A statistically equivalent phase screen realization can then be generated by imposing the above desired SDF on white noise:
6
where ηn is a zero-mean Gaussian random process with the Hermitian property.
Applying the phase screen realization 
and substituting x with the Fresnel scale, the propagation Equations (1) and (2) from the phase screen to the observation plane can be implemented in the simulation as
7
8
To convert the complex field ψ from space to time domain, an effective scan velocity veff is used such that y = vefft. The calculation of veff involves knowledge of the ionosphere anisotropy, propagation geometry including different angles formed between the line-of-sight signal and the geomagnetic field, and velocities including the effective scan velocities of satellite and receiver at the phase screen and the drift velocity of the ionosphere irregularities vdrift. For detailed veff calculation procedure, readers are referred to a previous study.33
The conversion from Doppler frequency to normalized wavenumber is
9
where fD is the Doppler frequency. A sampled phase screen constructed with P(μn)2πΔfD/Δμ = P(μn)ρF/veff where ΔfD = 1/(NΔt) and μn = nΔμ will generate a realization of time series (with a sampling interval Δt) that are statistically equivalent to the scintillation defined by the phase screen structure. As mentioned earlier, the sampling rate for the simulator is 100 Hz, corresponding to a Δt of 0.01 s. Intuitively, the ρF/veff directly determines how much compression or decompression the scintillation variations have over the temporal domain and is therefore tightly related to the decorrelation time of the scintillation intensity τ0, which will be shown in Section 6.3.
The time series of the scintillation induced phase and amplitude are then obtained from the propagated complex field as
10
11
To summarize, a GNSS complex-field scintillation realization can be generated by specifying the TPPSM parameters {U,p1,p2,μ0, and ρF/veff} and the sampling parameters Δt and N. Among the TPPSM parameters, U, μ0, and ρF/veff are dependent on the signal carrier frequency, whereas p1 and p2 remain the same for different frequencies. To construct the same phase screen for different frequencies, these three parameters are scaled from one frequency (denoted as the reference frequency fr) to another (denoted as the frequency with scaled parameters fs) using the following equations:
12
13
14
Using TPPSM parameters of different frequencies that satisfy the relationships in Equations (12) through (14) and the same white noise realization ηn in Equation (6), the simulator generates realistic multifrequency scintillation in terms of consistent scintillation level and correlated scintillation effects. This interfrequency consistency of the simulated scintillation is essential for the evaluation of multifrequency receiver algorithms36 and the study of multifrequency scintillation characteristics.29 To assess the performance of the model and the parameter estimation method, the parameters estimated from real multifrequency scintillation data are used to validate the relationships represented in Equations (12) through (14). This portion of the work will be presented in Section 6.1.
4 PARAMETER ESTIMATION
4.1 IPE
As mentioned earlier, the IPE technique is used in this study to extract the TPPSM parameters from real scintillation data and establish their distributions under a strong equatorial scintillation scenario. The IPE technique was developed in Carrano et al,28,30 while a brief introduction is provided in this section for completeness.
The IPE technique is essentially an iterative fitting procedure to obtain the TPPSM parameter estimates that yield the best match to the intensity SDF of the real scintillation data under the least square error (LSE) criteria. The IPE technique can also be implemented under the maximum likelihood criteria as presented in Carrano.37 In this paper, the LSE approach was chosen because of easier implementation and clarity of interpretation. It should also be mentioned that the IPE technique is applied to the intensity SDF rather than the phase SDF due to the higher quality obtained in the signal intensity measurements than that in the phase measurements during strong scintillation.38
The intensity SDF I(μ) observed on the receiver plane from a given phase screen can be theoretically evaluated as follows28:
15
where γ(·) is the structure interaction function defined in Carrano28 as
16
As can be seen in Equation (16), γ(·) is related to the phase screen SDF P(μ) defined in Equation (4) and, therefore, is also parameterized by the model parameter subset {p1, p2, μ0, U}.
Conversion for the intensity SDF from the normalized wavenumber domain (in μ) in Equation (15) to the temporal frequency domain (in fD) is achieved through
17
where ΦI,model is the intensity SDF in the temporal frequency domain.
The metric used to measure the fitting error between the model SDF ΦI,model and the measured SDF 
(the periodogram estimate of the real intensity measurement) is then defined as follows30:
18
where fmin and fmax indicate the frequency range [fmin, fmax] of the SDF over which the IPE fitting is performed, which is limited by the length N and update interval Δt of the intensity measurement by: 1/(NΔt) < fmin < fmax < 1/(2Δt).
A highly efficient algorithm for numerically computing I(μ) using Equation (15) as a function of {p1,p2,μ0,U} is presented in Carrano.28 Using Equation (17), the model intensity SDF in the time frequency domain ΦI,model can then be evaluated as a function of all five model parameters. Given an initial guess for the intensity SDF, the IPE program then iterates to obtain the optimal estimates of these parameters which provide the best fit of ΦI,model to 
by minimizing χ2. This multidimensional minimization is performed with the Nelder-Mead simplex method conveniently implemented in the Matlab function fminsearch.
The current IPE obtains the TPPSM parameter estimates from intensity measurements of different frequencies separately. This provides a means to evaluate the accuracy of the IPE method in terms of interfrequency consistency, which are defined by Equations (12) through (14) with real multifrequency data. This part of the work will be presented in Section 6.1.
4.2 S4 and τ0 calculation
As mentioned earlier, S4 and τ0 are the input for the simulator. To validate the parameter mappings, S4 and τ0 for both model and real data can be calculated from the intensity SDF as29
19
20
where ΦI = ΦI,model is the model intensity SDF (can be obtained using Equation (17)), and 
is the real measurement intensity SDF. RSI(·) is the normalized autocorrelation function (ACF) of the signal intensity measurements, which can be obtained as the DFT of ΦI.
5 DATA DESCRIPTION
The real scintillation GPS data sets used in this paper were collected at two equatorial sites: Hong Kong (geographic: 22.2° N, 114.3° E; geomagnetic: 16° N, 187.0° E) and Ascension Island (geographic: 7.9° S, 14.4° W; geomagnetic: 12.3° S, 57.0° E) using Septentrio PolaRxS ionospheric scintillation monitoring (ISM) receivers. The receivers are part of the global GNSS network deployed by University of Colorado Boulder for ISM and studies.17 For GPS civilian signals, the ISM receiver tracks L1 C/A signal, L2 CL signal, and L5Q signal and outputs various types of measurements including 100-Hz update rate carrier phase and signal intensity. Both receiver locations are close to the equatorial anomaly crest (around geomagnetic latitude 15°) where scintillation is known to be frequent and strong.6 Using a machine learning-based scintillation detection approach, we were able to identify a large amount of strong scintillation events collected from both sites.39
The signal intensity measurements are first divided into 5-min segments and then selected to extract the TPPSM parameters using IPE and {S4, τ0} values using Equations (19) and (20) for later analysis. The segment length of 5 minutes is heuristically chosen based on our experience working with a large amount of scintillation data. It is short enough to assume stationarity in scintillation effects, but also long enough to generate a reasonable intensity SDF estimate. It should be mentioned that the current simulator generates scintillation effects as a stationary process defined by the user input parameters during the simulation duration. In order to simulate a time-varying scintillation scenario resulting from the IPP scanning through different structures, an intuitive approach is to simulate the scintillation effects from different structures (defined by different parameter sets) separately and then combine the resulting time series. The authors are currently working on an updated version that incorporates such capability. Figure 2 shows 50-min data on 5 March 2014 from Hong Kong with strong scintillation observed on all three frequencies on GPS PRN 1. Strong scintillation with frequent deep fading exceeding 40 dB can be observed in the signal intensity on the top panel. The center of each 5-min segment is marked with Segment Numbers 1 to 9 at the bottom of the panel. The average S4 index of each segment is plotted in the middle panel, while the bottom panel shows the elevation of PRN 1.
50-min processing results on 5 March 2014 in Hong Kong with strong scintillation observed on all three frequencies on GPS PRN 1, starting at 13:27:00 UTC. The triangles mark the center of each 5-min segment that the intensity results are divided into
To ensure the measurement fidelity for IPE processing and scintillation strength, we imposed three criteria in data segment selection: (a) the segment average elevation is above 25°; (b) the average S4 index is over 0.6 on the L1 signal; and (c) the segment contains consistent signal fading to meet a good stationarity assumption. For example, among the 9 segments shown in Figure 2, it can be seen that Segments 3, 6, 7, and 9 meet all three criteria.
Based on these three criteria, the following three groups of scintillation intensity segments (forming a total of 174 segments) are selected from our scintillation data archives for IPE processing:
1. Hong Kong spring data: selected from 11 days during March 2014, containing 85 segments from seven satellites;
2. Hong Kong fall data: selected from 8 days during September through early November 2014, containing 55 segments from four satellites;
3. Ascension Island spring data: selected during 8 days of March 2013, containing 34 segments from four satellites.
It can be seen that all three groups of data were collected during the maximum of solar cycle 24. These three groups of data were chosen for comparison so that potential seasonal and longitudinal dependency (the two sites reside at similar geomagnetic latitudes) of the model spectral parameters can be observed. A subset of the data containing 45 segments with triple-frequency scintillation from groups 1 and 2 (on PRN 1, 3, 24, and 25) were selected to evaluate the IPE’s performance, which will be presented in the following subsection. The group 3 data did not contain any triple-frequency scintillation events from satellites above a 25 degree elevation.
6 REAL DATA ANALYSIS
6.1 IPE triple-frequency consistency evaluation
Figure 3 shows an example where the triple-frequency signal intensity from Segment 7 in Figure 2 have been fitted using the IPE method. The top panels plot the data intensity SDFs (blue) and the fitted SDFs (red). The intensity spectrum is computed by first breaking the 300-second segment into three consecutive 120-second segments with 40-second overlaps in between and then averaging their periodograms. Only the midfrequency sections of the data intensity SDFs are used for the fitting, as indicated in the figure with the dashed, double-sided arrows. The lowest frequencies may be distorted by large-scale departures from stationarity, and the highest frequencies may be contaminated by receiver noise.28 Both are therefore excluded from the fitting. The S4 and τ0 values estimated from the data and fitted model are also listed. The bottom panels plot the data intensity ACFs and the fitted ACFs.
SDF and ACF for intensity scintillations observed on GPS triple-frequency signals from segment 7 in Figure 2. The SDFs and AFCs for the measurements are shown in blue, while those for the fitted model are shown in red. The dashed, double-sided arrow on each top panel indicates the frequency range over which the IPE fitting was performed
From Figure 3, it can be seen that the shapes of the IPE-derived SDF and ACF are both very good fits to those of the data for all three frequencies, and the S4 and τ0 values of the IPE results and of the data are also all in good agreement. This example shows that IPE is effective in fitting the observed data’s SDF with TPPSM.
The agreement in the spectra fitting and the resulting S4 and τ0 values is representative throughout the IPE processing of the three selected data groups. In order to quantitively evaluate the IPE method’s performance in estimating these model parameters, the consistency across the carriers in the IPE-generated results were assessed using the triple-frequency data set. As mentioned in Section 3 with Equations (12) through (14), the five model parameters have different theoretical interfrequency relationships that can be summarized as follows:
1. p1 and p2: These two parameters are not frequency dependent. In order to evaluate their frequency independence, three 2-D scatter plots can be formed with IPE estimates of three pairs of frequencies (L1-L2, L1-L5, and L2-L5). A linear relationship is then fitted with slope a and vertical offset b. This linear relationship is compared against the theoretical linear function with slope 1 and offset 0. The resulting deviation can be described by the differences in the slopes (Δa = a − 1) and vertical offset (Δb = b);
2. μ0 and ρF/veff: According to Equations (12) and (13), these two parameters have the same frequency dependency of
, with g = μ0 or ρF/veff. Again, 2-D scatter plots are shown for each of the three pairs of frequencies for μ0 and ρF/veff. A linear relationship can be fitted for each scatter plot. For a certain pair Li-Lj, the deviation between the fitted relationship and the theoretical relationship can be expressed as and Δb = b;3. U: The interfrequency relationship of U is the most complicated one according to Equation (14). It has four different forms depending on the μ0 values on both frequencies in question. For the triple-frequency data subset used in this paper, the IPE-estimated μ0 values were all smaller than 1 for all three frequencies. This result limits U’s theoretical interfrequency relationship to the fourth form in Equation (14) as
. This relationship is dependent on the frequencies and the p2 value of the phase screen spectrum, which obviously varies among different data segments. Let us denote the IPE estimates of U from real data segments of a certain frequency Li as . As a result of the p2 dependency, the interfrequency estimate pairs from different data segments of the set cannot be directly combined into a scatter plot to fit a linear relationship for L1-L2, L1-L5, and L2-L5.
In order to assess the interfrequency consistency of U estimation in a similar fashion as the other four parameters, for a certain frequency pair Li − Lj, a scatter plot is shown with 
still assigned to the vertical axis. For the horizontal axis, instead of directly using 
, a 
-based prediction of U(Lj) (denoted as 
) is employed. 
is yielded using 
, where 
is the p2 estimate from each data segment. 
and 
should then be theoretically equivalent, and a linear relationship can be fitted from the resulting scatter plot. The corresponding correlation coefficient and deviation (as Δa = a − 1 and Δb = b) in turn reflect how close the 
relationship agrees with theory. It should be mentioned that 
used in calculating 
scaling is the mean value of the triple-frequency p2 estimates, in order to minimize the p2 estimation error introduced into the scaling.
In Figure 4, the interfrequency scatter plots summarized above for all five parameters are plotted with blue markers, each with three pairs (L1-L2 in the left panel, L1-L5 in the middle panel, and L2-L5 in the right panel). Within each panel, the least-square-fitted linear relationship and the model-defined relationships are plotted in red and black, respectively, with the correlation coefficients and the deviation Δa and Δb between the linear-fitted and the model relationships listed.
The interfrequency relationships for all five parameters, each with three pairs (L1-L2, L1-L5, and L2-L5). Real data results are plotted in blue markers, and the least-square-fitted and the model-defined ones are in red and black, respectively. The correlation coefficients, the differences in the slope (Δa) and offset (Δb) between the linear-fitted, and the model relationships are also listed
As can be seen in Figure 4, all five parameters have generally high correlation coefficients among all three pairs of carriers results and small deviations (Δa and Δb) from the model, which validate the model-defined linear relationships and IPE’s performance. Among the results of three different frequency pairs, the L2-L5 pair generates the best results in both correlation coefficients and model deviations for all five parameters. This is most likely because L2 and L5 have stronger scintillation effects than L1, and, therefore, lead to a more dominant contribution in their own intensity SDFs over other error sources, which improves the IPE fitting accuracy. Among the results of the five parameters, the ρF/veff showed the best results in both correlation coefficients and model deviations, suggesting the highest sensitivity in fitting. This finding sheds light on the promising usage of the IPE as a way of inferring ionospheric background information such as vdrift for a standalone receiver, as has been mentioned in previous studies.28-30
6.2 Spectral parameter characterization
This section presents the distributions of the L1 signal spectral parameters {p1,p2,μ0} obtained by IPE from the three groups of data sets to determine their most representative values and potential seasonal and longitudinal dependency. Figure 5 plots the distributions of p1 (left panel), μ0 (middle panel), and p2 (right panels) from Group 1 (denoted as G1 in green), Group 2 (as G2 in blue), and Group 3 (as G3 in red) with their mean values listed in the corresponding legends.
Distributions of the IPE-estimated spectral parameters from the three data groups (Group 1 in green, Group 2 in blue, and Group 3 in red), with the mean values listed in the legends
By examining Figure 5, there are no considerable discrepancies among the distributions of the parameters from the three data groups for all three parameters, suggesting no apparent seasonal or longitudinal dependency of these parameters. The ranges for p1 and μ0 lie from 1.4 to 3 and 0.3 to 0.7, respectively, while p2 generally ranges between 3 and 4.6. The results from all three groups are, therefore, combined to calculate the mean for each parameter as their representative values for the general strong equatorial scenario (p1 = 2.45, μ0 = 0.55, and p2=3.70).
It should be mentioned that the authors have investigated the correlation between parameter pairs for all five parameters in all three data groups, and no apparent correlations have been observed (the highest correlation coeffect was 0.41, which was between p1 and U estimates in Group 3).
6.3 Parameter mappings and real data evaluation
After defaulting the spectral parameters to their representative values obtained from the previous subsection, the mappings between {S4, τ0} and {U, ρF/veff} can be established by numerically evaluating Equations (15), (17), (19), and (20).
By examining Equations (15), (17), and (19), the effect of different parameters from the whole model parameter set {U,p1,p2,μ0,ρF/veff} on the S4 index calculation can be broken down as follows:
1. According to Equation (19), the S4 value is directly determined by the shape and magnitude of ΦI(f);
2. By looking at Equation (17), ΦI(f) is obtained by linearly mapping I(μ) into the time frequency domain with the space-to-time scaling factor of ρF/veff;
3. {p1,p2,μ0,U} then affect the calculation of S4 as the shape of I(μ) is jointly determined by {p1,p2,μ0,U}, while its magnitude is directly determined by U, as demonstrated in Carrano28;
4. The value of ρF/veff affects the calculation of S4 as it determines which part of the I(μ) corresponds to the part of ΦI(f) being integrated over the range of [fmin, fmax] in Equation (19). This is indeed the case for scintillation effects on platforms with high dynamics and in an orbit close to the ionosphere (such as LEO satellites), where the value ρF/veff may have a small magnitude of 10−3 and become the dominant factor within Equation (19).32 However, since these numerical mappings are meant to relate stationary ground-observed scintillation indicators to the model parameters, the range of ρF/veff is limited to 0.5 to 2.0 as observed in the three data groups. Such limited variability does not affect the I(μ) to ΦI(f) mapping sufficiently to impact the outcome of Equation (19). Therefore, in the numerical mapping established for this simulator, the S4 index is independent of ρF/veff and only related to U.
Figure 6 shows the numerically evaluated relationship (denoted as model, in red) and the scatter plot of the S4 – U estimate pairs from the combined scintillation data set (denoted as data, in black markers). In addition, an iterative fitting is performed in search for the optimal {p1,p2,μ0} values that can achieve the LSE with respect to the real data’s S4 – U pairs. The resulting S4 – U numerical relationship is also plotted in Figure 6 in blue, and the optimal {p1,p2,μ0} values under LSE are listed in the legend.
The numerically evaluated relationship (denoted as model, in red) and the scatter plot of the S4 index and U estimates from the combined scintillation data set (denoted as data, in black markers), as well as the relationship that resulted from fitting the p1,p2,μ0 to optimal values that satisfy the LSE criteria with respect to the data (denoted as fitted, in blue)
As can be seen in Figure 6, the LSE fitted values for {p1,p2,μ0} are very close to their corresponding representative values obtained in the previous subsection, and the resulting two relationships are in close agreement with each other and appear to lie in the center of the variation range of the real data’s S4-U pairs. This validates the numerical mapping between the S4 index and U.
As for τ0, the numerical evaluation shows that the ρF/veff is linearly correlated with τ0, and the slope of this linear relationship varies considerably with respect to the variations of U within its typical range (usually 1-3 for strong L1 scintillation). Figure 7 plots the numerically evaluated relationships between τ0 and ρF/veff color coded with respect to different U values ranging from 1 to 20 with a spacing of 1.
The numerically evaluated relationships between τ0 and ρF/veff color coded with respect to the corresponding U values
As can be seen in Figure 7, τ0 and ρF/veff are linearly correlated with a slope that is determined by the value of U. As U increases, the slope becomes smaller, which means that under the same ρF/veff value, more severe scintillation will result in faster decorrelation. In order to validate these linear relationships between τ0 and ρF/veff, the real data set is divided into four categories according to their ranges in U: [0.6-1], [1.0-1.5], [1.5-2.0], and [2.0-]. For each category, a numerical relationship between τ0 and ρF/veff is obtained under the mean U value of the data samples within the corresponding category.
Figure 8 plots the scatter plot of the τ0 versus ρF/veff estimates from the four categories of the real data distinguished by color and marker shapes. The numerical relationship under the mean U value of each data category is plotted as a straight line in the same color as the scatter plot of the corresponding category. The U range and its mean value of each data category are listed in the legend. Figure 8 clearly shows close agreement between the model evaluated relationship and the τ0 and ρF/veff relationship in real data, despite the variability of U within each data category.
The scatter plot of τ0 versus ρF/veff estimates from the four categories of the real data distinguished by color and marker shapes. The numerical relationship under the mean U value of each data category is plotted as a straight line in the same color as the corresponding category. The U range and its mean value of each data category are listed in the legend
After obtaining the representative values to default the model spectral parameter subset {p1,p2,μ0} and the numerical mappings between {S4, τ0} and {U, ρF/veff}, the simulator can now be controlled by only specifying the expected {S4, τ0} values. An example of triple-frequency signal intensity is simulated using the same {S4, τ0} values as estimated from a 5-min real data segment collected on PRN 24 in Hong Kong starting from 12:25:00 UTC on 5 October 2013. The intensity and phase measurements from both the real and simulated data are shown in Figure 9 with the corresponding {S4, τ0} values for each signal listed in the legend. The simulated signal intensity and phase measurements do not contain any low-frequency trend and, therefore, do not need detrending operations, while the detrending for the real data intensity and carrier was performed with a wavelet method as described in Jiao et al.40 It should be mentioned that the carrier phase measurements recorded by commercial receivers during strong scintillation usually experience numerous large cycle slips that require advanced cycle slip repair procedures, which is outside the scope of this paper.41
Comparison of the signal intensity and phase in the simulated data and real data
As can be seen in Figure 9, the simulation data and real data show very similar characteristics in a statistical sense on both the signal intensity and carrier phase. The {S4, τ0} values between the real and simulator-generated data are also in good agreement for signals of all three frequencies. This validates the effectiveness of the two-parameter scintillation generator in simulating realistic, coherent triple-frequency scintillation according to the users’ expected scintillation condition. A more comprehensive comparison between the simulated and real data have been conducted in terms of scintillation characteristics, including the amplitude fading duration and separation and the multifrequency fading consistency. The results further confirmed the simulator’s effectiveness in generating realistic scintillation signals. This comprehensive validation is the topic of a new manuscript which is ready for journal submission.
7 CONCLUSIONS
This paper presented a TPPSM-based multifrequency strong scintillation simulator that requires only the expected scintillation index S4 and the intensity decorrelation time τ0. This simulator was developed by defaulting three of the TPPSM spectral parameters {p1, p2,μ0} to representative values and therefore obtaining numerical mappings from the user input parameter set {S4,τ0} to the remaining parameter subset {U,ρF/veff}. The numerical evaluation shows that the S4 has a one-to-one mapping relationship with U, whereas τ0 is linearly correlated with ρF/veff with a slope dependent on U’s value. Based on these numerical mappings, the scintillation model can then be controlled by specifying the expected S4 and τ0 values from scintillation measurements obtained from a stationary ground-based receiver. The simulator can generate statistically equivalent realizations of the scintillation effects. Using this basic two-parameter set, the simulator can also generate realistic scintillation effects that are observed on platforms with user-defined dynamics under the same ionospheric phase screen responsible for producing the ground received scintillation. A MATLAB implementation has been made available for users, and the downloading information and configuration guide are provided in the appendix.
A total of 174 5-min segments of strong scintillation data (with S4 > 0.6) from two equatorial sites were processed using IPE to establish the profiles of the model parameters. Based on the profiles of the three spectral parameters from data collected from different locations and seasons, the distributions of the three parameters indicate that there are no clear longitudinal and seasonal dependencies. Therefore, the profiles are then combined to yield the defaulting values of the three spectral parameters {p1,p2,μ0} used in obtaining the parameter numerical mappings, which are eligible to represent the most typical equatorial strong scintillation scenario. Based on the estimates of {U,ρF/veff} and {S4,τ0} obtained from this multisite real scintillation data set, numerical mappings were validated to be generally accurate to represent the case in observed strong equatorial scintillation. In addition, this paper also evaluated the IPE method in terms of triple-frequency consistency using a subset of the real scintillation data set with triple-frequency scintillation, which confirmed the IPE’s effectiveness.
HOW TO CITE THIS ARTICLE
Xu D, Morton YTJ, Rino CL, Carrano CS, Jiao Y. A two-parameter multifrequency GPS signal simulator for strong equatorial ionospheric scintillation: modeling and parameter characterization. NAVIGATION. 2020; 67:181–195. https://doi.org/10.1002/navi.350
ACKNOWLEDGEMENTS
The receiver deployed in Ascension Island was funded by a grant through AFRL (FA8650-08-D-1451), Wright Patterson Air Force Base. The authors appreciate the support of Dr Todd Pederson from the Air Force Research Laboratory at Kirtland AFB that made this deployment happen. The receiver deployed in Hong Kong was funded by the Consortium of Ohio Universities on Navigation and Timekeeping (COUNT), and the authors would like to thank Dr Zhizhao Liu from Hong Kong Polytechnic University for hosting the receiver. The work conducted in this paper was funded by Colorado State University, University of Colorado Boulder, and a NASA Grant (NNX15AT54G).
APPENDIX A. SIMULATOR PROGRAM CONFIGURATION AND EXECUTION
A MATLAB implementation of this two-parameter scintillation simulator can be downloaded from https://github.com/cu-sense-lab/gnss-scintillation-simulator_2-param. The program folder contains the.m file of function Main and the folder named “Libraries.” The “Libraries” subfolder contains the supporting functions for function Main.
The user-input parameters to be specified are all included in Main.m:
1. Simulation start date and time in UTC;
2. Simulation length in seconds;
3. Satellite PRN number;
4. S4 index and τ0;
5. Platform position in latitude-longitude-height (LLH) coordinates and velocity vector in east-north-up (ENU) coordinates to generate stationary scintillation data, simply assign the velocity vector to all zeros;
6. Number of frequencies to be simulated: 1—GPS L1 only||2—L1 and L2||3—L1, L2, and L5; and,
7. Figure plotting switch (if on, figures to be plotted include satellite sky plot, ionosphere piercing point [IPP] LLH 3-D plot, simulated scintillation intensity and phase).
The variables for each input parameter are named in function Main in a self-explanatory manner and stored in the struct “userInput.” Comments are also provided above the value assignment of each variable to explain each variable and the correct input format. After specifying each parameter, run the Main function, and the scintillation amplitude and phase time series of specified length with an update rate of 100 Hz will be recorded in the variables “Scin_amp” and “Scin_phi.”
For Windows or Mac users, if the program reports an error during the function ExtractRINEXeph.m which downloads the RINEX ephemeris file and extracts ephemeris parameters for satellite orbit computation, simply go to the simulator program folder and unzip the downloaded brdc****.**n ephemeris file by hand and then rerun the program.
Footnotes
Funding information AFRL, Grant/Award Number: FA8650-08-D-1451; Colorado State University; Consortium of Ohio Universities on Navigation and Timekeeping (COUNT); National Aeronautics and Space Administration, Grant/Award Number: NNX15AT54G; University of Colorado Boulder
- Received June 5, 2019.
- Revision received August 28, 2019.
- Accepted October 17, 2019.
- © 2020 Institute of Navigation
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
- Article
- Abstract
- 1 INTRODUCTION
- 2 TWO-PARAMETER SCINTILLATION SIMULATOR ROUTINE DESCRIPTION
- 3 PHASE SCREEN REALIZATION AND WAVE PROPAGATION
- 4 PARAMETER ESTIMATION
- 5 DATA DESCRIPTION
- 6 REAL DATA ANALYSIS
- 7 CONCLUSIONS
- HOW TO CITE THIS ARTICLE
- ACKNOWLEDGEMENTS
- APPENDIX A. SIMULATOR PROGRAM CONFIGURATION AND EXECUTION
- Footnotes
- REFERENCES
- Figures & Data
- References
- Info & Metrics
Related Articles
Cited By...
- No citing articles found.