Abstract
Previous authors have reported on the morphology of GPS scintillations and irregularity zonal drift during the 2002 Conjugate Point Equatorial Experiment (COPEX) in Brazil. In this paper, we characterize the turbulent ionospheric medium that produced these scintillations. Using 10 Hz GPS carrier-to-noise measurements at Boa Vista (2.9°N, 60.7°W), Alta Floresta (9.9°S, 56.1°W), and Campo Grande (20.5°S, 54.7°W), we report on the variation of turbulent intensity, phase spectral index, and irregularity zonal drift as a function of latitude and local time for the evening of 1-2 November 2002. The method of analysis is new and, unlike analytical theories of scintillation based on the Born or Rytov approximations, it is valid when the scintillation index saturates due to multiple-scatter effects. Our principal findings are that (1) the strength of turbulence tended to be largest near the crests of the equatorial anomaly and at early postsunset local times, (2) the turbulent intensity was generally stronger and lasted two hours longer at Campo Grande than at Boa Vista, (3) the phase spectral index was similar at the three stations but increased from 2.5 to 4.5 with local time, and (4) our estimates of zonal irregularity drift are consistent with those provided by the spaced-receiver technique.
1. Introduction
The distribution of free electrons in the ionosphere is dictated by production from solar radiation, transport, and loss through chemical recombination. It is also subject to instability mechanisms that generate large-scale depletions and irregularities in the ambient electron density over a wide range of spatial scales (plasma turbulence). Radio waves that propagate through these irregularities experience scattering and diffraction, causing random fluctuations in amplitude and phase referred to as scintillations. The scintillation of satellite signals has been shown to severely degrade the performance of satellite communications systems such as AFSATCOM [1, 2], satellite global navigation satellite systems (GNSS) such as the Global Positioning System (GPS) [3–5], and space radars used to conduct cloud-free, day-and-night observations of the Earth’s surface [6–8]. Scintillations associated with irregularities in the equatorial ionosphere are generally the most intense encountered worldwide, and the occurrence morphology depends on local time, season, longitude, solar cycle, magnetic activity, and exhibits a high degree of night-to-night variability [9]. Ionospheric irregularities and scintillations constitute one of the most important space weather threats to modern technological systems which increasingly rely on transionospheric radio propagation. While a multitude of plasma instabilities may operate in the equatorial ionosphere after sunset, interchange instabilities which generate large-scale depletions in the ambient electron density, commonly referred to as equatorial plasma bubbles (EPB), are believed to be the dominant source of irregularities that cause scintillation of the GPS satellite signals at L band frequencies [2–5].
The Conjugate Point Equatorial Experiment (COPEX) was conducted in Brazil from October to December 2002. Its purpose was to explore the occurrence morphology of scintillations and the physical processes by which plasma instabilities occur in the equatorial ionosphere. In this experiment, a multitude of ionospheric monitoring instruments, including ionosondes, optical imagers, VHF receivers, and GPS receivers, were deployed to three sites nearly along the same magnetic meridian, one at the magnetic equator and the other two at magnetically conjugate points.
Several authors have reported on the morphology of GPS scintillations and irregularity zonal drift during this experiment, for example, Batista et al. [10], Muella et al. [11], Abdu et al. [12], Sobral et al. [13], and de Paula et al. [14]. In this paper, we characterize the turbulent ionospheric medium that produced these scintillations. For this purpose, we developed a new technique called Iterative Parameter Estimation (IPE) for inferring ionospheric turbulence parameters from a time series of scintillating intensity measurements resulting from radio propagation through the low latitude ionosphere. Unlike analytical theories of scintillation based on the Born or Rytov approximations, the IPE technique is valid for strong scatter when the scintillation index () saturates due to multiple-scatter effects. This generality is crucial for proper analysis of GPS scintillations during the COPEX campaign, since a large percentage of these observations were in the strong scatter regime, with frequently approaching or exceeding unity [11]. Weak scatter theory cannot be applied satisfactorily to analyze these observations. It is also beneficial that the IPE technique requires only the intensity fluctuations to infer the turbulence parameters. While the GPS phase measurements from the COPEX campaign are also available, their analysis is complicated by frequent cycle slips and loss of phase lock events which occur when the scintillation is strong. IPE employs the phase screen approximation [15–18] and a numerical inversion technique to estimate the parameters of the screen which are consistent with the ionospheric turbulence model employed and the scintillations observed. Due to the analytic intractability of the 4th moment equation governing fluctuations in the field we must determine the optimal parameters by numerical iteration, and therefore we cannot guarantee the uniqueness of this parametrization.
The purpose of this paper is to explore the latitudinal and local time variation of ionospheric turbulence characteristics during the COPEX campaign in Brazil. The IPE technique was developed to investigate these variations using ground-based observations rather than in situ measurements of the ionospheric turbulence, which are limited in several respects. First, the temporal coverage of in situ measurements within any geographic region is sparse, whereas ground-based observations of scintillation are continuous. Second, in situ observations of the density are usually made in the topside ionosphere, whereas the scattering of L band signals is generally believed to be concentrated near the F region peak where the density is highest. It is not obvious how to infer the level of density fluctuations near the F region peak from in situ measurements made in the topside (and it is clearly not possible to do so when the irregularities have not risen to the orbital altitude of the satellite making the in situ observations). An advantage of using ground-based GPS scintillation observations to characterize the morphology of ionospheric turbulence is that irregularities at all ionospheric altitudes are sampled by the satellite signals.
The organization of this paper is as follows. Section 2 describes the equipment and the measurements used in our analysis. Section 3 presents the methodology used to analyze the data. The exposition of the methodology is given in three parts. First, we describe a model for the correlation of phase fluctuations after the GPS satellite signals have penetrated the ionosphere. Second, we describe a technique that uses this correlation function to model the spectrum of intensity fluctuations on the ground. Third, we describe an iterative technique used to infer the turbulence parameters by fitting this model spectrum to the spectrum of measured intensity fluctuations. In Section 4, we present the results of our analysis and describe their significance. A summary of these results is provided in Section 5.
2. Experimental Techniques
During the COPEX experiment in 2002, Air Force Research Laboratory operated Ashtech μZ-CGRS model dual frequency GPS receivers at three locations in Brazil: Boa Vista (2.9°N, 60.7°W), Alta Floresta (9.9°S, 56.1°W), and Campo Grande (20.5°S, 54.7°W). Boa Vista is located near the northern crest of the equatorial anomaly, Alta Floresta is located near the magnetic equator, and Campo Grande is located near the southern crest of the equatorial anomaly. Figure 1 shows the location of the three stations and the local magnetic field geometry, while Table 1 provides their geographic and geomagnetic coordinates. All three stations lie on approximately the same magnetic meridian. The three AFRL GPS receivers recorded the carrier-to-noise ratio (C/No) of the C/A code on the L1 frequency of 1575.42 MHz at a rate of 10 samples per second. The data was collected in October through December of 2002, a period of high solar activity. Additional details of the COPEX experiment are summarized in [10–14] and the references therein. For validation purposes, we compare our results with the zonal irregularity drift measurements calculated by Muella et al. [11] using the GPS spaced-receiver technique, and also the 4-channel VHF receivers which were colocated with each of the GPS receivers. The spaced VHF receiver measurements of the zonal irregularity drift provided to us for this study were obtained by simple cross-correlation of the intensity fluctuations obtained along geostationary links from two receivers separated in the magnetic east-west direction. The geometrical corrections described in [19], which account for the look angles to the geostationary satellites and the magnetic dip angles at the penetration point locations, were not applied. The GPS spaced receiver measurements provided by Muella et al. [11] did account for these geometrical effects. Unfortunately, neither the VHF nor GPS spaced-receiver measurements of the zonal drift were provided to us with error bars. Table 2 gives the locations of the 350 km ionospheric penetration points and look angles for the VHF west links (receiver channels 1-2) and east links (receiver channels 3-4) for the three COPEX stations.
In this paper, we consider measurements collected during the evening of 1-2 November 2002, since during this period all three AFRL GPS receivers were operating, all three VHF spaced receivers were operating, and the spaced GPS receiver estimates of the zonal drift velocity were available for validation. This evening was typical of others during the COPEX experiment, with moderate magnetic activity (Kp ranged from 2° to 5−) and the 10.7 cm solar flux was 160 W/m2/Hz. Unfortunately, there were only a few other evenings during the COPEX campaign when all of these instruments were operating simultaneously, and most of these evenings were more magnetically active than the evening of 1-2 November 2002. Following [14], we restricted our analysis to evenings with low-to-moderate magnetic activity, since it is easier to compare different techniques for estimating the zonal irregularity drift when the actual drift is regular and changing slowly in the absence of storm-time perturbation electric fields.
For our analysis, we use the receiver reported C/No as a proxy for the signal intensity. From the time series of raw intensity fluctuations, we selected data segments which are approximately stationary and for which the scintillation intensity index, calculated as the standard deviation of normalized (by the mean) signal intensity, exceeded 0.3. This threshold criterion was applied to minimize the contribution of receiver noise to the scintillation statistics. A satellite elevation cutoff of 30° was used to avoid multipath and to reduce possibility that the radio wave may have propagated through multiple plasma bubbles with different turbulence characteristics. The stationary requirement is enforced to ensure the time series admits a spectral representation. We chose to use 4-minute data segments for the analysis, rather than 1-minute data segments which have become rather customary in scintillation studies, for example, [3]. The reason for this choice is that IPE analysis provides the most accurate results when all frequencies that contribute to the power spectrum of intensity fluctuations are resolved, including frequencies somewhat smaller than the Fresnel scale. Resolving smaller frequencies requires the analysis of a longer data segment. Over these relatively longer segment durations, however, the stationary requirement becomes more important to enforce. To confirm the (approximate) stationarity of the data, we calculated the index on a 1-minute cadence and also on a 4-minute cadence. Stationary segments were identified as segments for which 4 consecutive 1-minute values differed from the corresponding 4-minute by less than 0.1. The mean intensity was observed to vary slowly over a 4-minute timescale for satellites viewed at elevation angles above the 30° cutoff. The time to 50% decorrelation of intensity () was calculated from the segmented time series using the FFT technique. We computed the temporal spectrum of intensity fluctuations from the measured time series using Welch’s method [20], whereby each 4 minute data segment is subdivided into 1-minute subsegments. We apply a Hamming window to each of these 1-minute subsegments prior to FFT analysis and then average these spectra with 50% overlap. The averaging is performed to minimize noise in the spectra due to spectral leakage [20]. The resulting power spectral densities range from the minimum resolvable frequency of 1/60 sec = 0.0167 Hz to the Nyquist frequency of 10 Hz/2 = 5 Hz.
3. Method of Analysis
3.1. Ionospheric Turbulence Model
Following the development by Rino [16], a continuously displaced coordinate system is chosen in which the measurement plane follows the principal propagation direction. The propagation geometry is shown in Figure 2. The coordinate system is defined with origin at the ionospheric penetration point (IPP), (,, and ), which is located at the center of the ionospheric layer. At each point along the propagation path the , , and axes point toward geomagnetic north, geomagnetic east, and downward, respectively. The thickness of the ionospheric layer, assumed to contain homogenously distributed irregularities, is . The angle θ is the propagation (nadir) angle at the IPP, and φ is the magnetic meridian angle of the propagation vector.
As in [16], the 3D spatial spectrum of electron density fluctuations () is modeled as a power law with an outer scale as follows: In (1), is the strength of turbulence, is the irregularity spectral index, and is the isotropic (radial) wavenumber. The outer scale wavenumber, , is related to the outer scale of the turbulence, , as . In this model, is required for the variance of electron density fluctuations to remain finite. For irregularity scales such that this spectrum reduces to the pure power law form . In the equatorial region, electron density irregularities are highly elongated and aligned with the geomagnetic field. To account for irregularity anisotropy, Rino [16] employed a scaling and rotation of the coordinate system; this transformation is described in considerably more detail in [18]. Application of this transformation to (1), followed by straight-line integration through the ionospheric layer along the line of sight, followed by Fourier transformation, gives an expression for the correlation function of phase fluctuations just beneath the screen: In (2), is the classical electron radius, is the radio wavelength, is the separation distance in the transverse plane, is the modified Bessel function of order , and is Euler’s gamma function. Equation (2) above was first derived in [16] and appears in that paper as (11). G is a geometric enhancement factor which is given by Equation (3) also appeared in [16] but with a typographical error; the factor cos should be outside the radical as shown in (3). The scaling factors and in (3) elongate contours of constant phase correlation along and transverse to the magnetic field, respectively. The coefficients , , and depend on the direction of propagation and the orientation of the irregularity axes [16]: where In (5), is the magnetic inclination angle and is the angle at which irregularities are inclined from the xz plane (which we take to be zero). We use the International Geomagnetic Reference Field (IGRF) 2000 [21] to compute these parameters at the location of the ionospheric penetration points, using an assumed ionospheric shell height of 350 km.
It is convenient, at this point, to introduce the vertically integrated strength of turbulence at the 1 km scale [22]: It is common to evaluate the turbulence strength in this form since the layer thickness () and turbulent intensity appear together as a product and need not be known independently [5, 22–24]. Characterizing the strength of scatter using , rather than , is advantageous because the former is a property of the random medium alone, whereas the latter depends on the random medium, propagation geometry, and frequency of the radio wave.
Up to this point, the ionospheric turbulence model involves purely spatial quantities, whereas a time series of field fluctuations is measured by the receiver during an experiment. Assuming the random medium is invariant over the measurement interval (this is the Taylor hypothesis of frozen-in flow), spatial fluctuations and temporal fluctuations can be related by a model-dependent effective scan velocity. The effective scan velocity consistent with the Rino spectral model is [16, 18] Assuming the irregularity drift is purely in the zonal direction, the scan velocity in the continuously displaced coordinate system is given by where is the zonal irregularity drift velocity, and , , are the velocity components of the ionospheric penetration point (IPP) in the magnetic north, magnetic east, and down directions, respectively. In (8), the terms involving tan θ account for the horizontal translation of the continuously displaced coordinate system as it follows the propagation ray path. We note that the effective scan velocity has a particularly simple interpretation for the special case of normal propagation through infinitely long irregularities; in this case is equal to the zonal irregularity drift velocity minus the zonal component of the IPP velocity. The general case of oblique propagation and finite irregularity axial ratio is used for the calculations in this paper. Following the approach implemented by the Wideband Scintillation Model (WBMOD) [22] and also [25], the component of IPP velocity parallel to the line of sight has been removed prior to the calculation (8) in order to reduce spectral smearing. We compute the GPS satellite locations and velocities using the SGP4/SDP4 satellite propagator from Spacetrack [26].
In the next section, we use the model for the correlation function of phase fluctuations given in (2) and phase screen theory to compute the model spectrum of intensity fluctuations in the reception plane.
3.2. Calculation of the Model Intensity Spectrum
A widely employed simplifying assumption in transionospheric radio propagation problems is the phase screen approximation. In the phase screen approximation, the interaction of the wave with an extended irregularity layer is replaced by that with an equivalent thin phase-changing screen [15–18, 27]. This approximation neglects diffraction effects that develop within the irregularity layer, which can be considerable when the scattering is strong. Nevertheless, several authors (e.g., see [27] and the references therein) have demonstrated that an equivalent phase screen accurately reproduces the amplitude and phase scintillation predicted by a more complex formulation that accommodates diffraction within the scattering layer provided that the height of the phase screen is appropriately chosen. In this paper, we employ the phase screen approximation so that the ionospheric turbulence may be characterized using a small number of parameters (i.e., those parameters which specify the equivalent phase screen).
When a plane wave is incident on a thin phase-changing screen, the spectrum of intensity fluctuations at the reception plane, at some vertical distance beyond the screen, may be calculated using the algorithm proposed by Booker and MajidiAhi [28]. This algorithm provides an intensity spectrum that satisfies the differential equation governing the 4th moment of the field in the thin screen approximation (e.g., see Ishimaru [15, Section 20–14]). The Booker and MajidiAhi approach is not limited to weak scatter conditions, as are the analytical results relating the screen parameters to given in [16]. These analytical results are implemented in the Wideband Model [22, 29], which applies an empirical correction for strong scatter based on the questionable assumption of Rician statistics. By solving the 4th moment equation directly there is no need to assume Rician statistics when the scatter is strong. While the 4th moment equation has been solved in various other forms (e.g., [30–32]), Booker et al. [27] explains that the Booker and MajidiAhi formulation [28] provides the most accurate results for the same level of computational effort. In this section, we adapt the Booker and MajidiAhi approach so that it may be used with Rino’s model for the correlation function of phase fluctuations (2).
We begin by defining the Fresnel parameter, , in terms of the signal wavelength and the slant propagation distance, : Next, we form two functions that depend on the separation distance in the transverse plane, , and the spatial wavenumber, : Finally, the spectrum of intensity fluctuations is given by The factor of 2 appears in (11), instead of 4 as in the Booker and MajidiAhi paper [28], in order to change from a one-sided spectrum to a two-sided spectrum for consistency with Rino [16]. The integral in (11) is oscillatory and challenging to evaluate. This is especially true when the scatter is weak and Fresnel ringing is present in the spectrum. We use an adaptive bisection algorithm from QUADPACK [33] to evaluate (11). Once the intensity spectrum has been computed, the scintillation index may be calculated by integrating the spectrum over all wavenumbers: Once again, we have a formulation that involves spatial quantities, while a time series of intensity fluctuations are measured during an experiment. It can be shown by a change of variables that when temporal frequencies in the time series of intensity fluctuations and spatial wavenumbers are related as then temporal power spectral densities and spatial power spectral densities are related as Once the temporal intensity spectrum is known, the decorrelation time () corresponding to the model spectrum can be evaluated by Fourier transforming to obtain the intensity correlation function and then determining the time lag to 50% decorrelation.
3.3. Iterative Parameter Estimation
The preceding sections have described a technique to compute the temporal spectrum of intensity fluctuations corresponding to a given ionospheric screen. The free parameters of the model are the altitude of the phase screen (), the anisotropy parameters (, and ), the outer scale (), the turbulence strength (), irregularity spectral index (), and the zonal irregularity drift velocity (). All the other parameters in the model can be computed from the signal frequency and propagation geometry. In this paper we assume , , and = 10 km. The results of a sensitivity study revealed that the intensity spectrum so computed is insensitive to the value of once the ratio exceeds approximately 20. In general, we cannot unambiguously measure from the GPS observations alone, because is much larger than the Fresnel break scale. Scales larger than the Fresnel break scale are naturally filtered out by diffraction and therefore contribute little to the intensity fluctuations on the ground (except in the case of strong focusing). In a second sensitivity study it was determined that the phase spectral index, , as determined using the IPE technique, is insensitive to the value chosen for so long as the latter is much larger than the Fresnel scale. The effective scan velocity, , and therefore also our estimates of the zonal irregularity drift velocity (), depend on the assumed altitude of the phase screen. We assume the fixed value = 350 km for the calculations in this paper. With these assumptions, the remaining independent parameters are , , and ; these parameters specify the ionospheric screen and will be estimated from the scintillation observations by fitting the model intensity spectrum to the observed intensity spectrum as described next.
We define a metric to measure the difference between the modeled intensity spectrum () and measured intensity spectrum () in the log-log domain: Given an initial guess for the screen parameters , , and , we determine the optimal values of these parameters which provide the best fit of the model intensity spectrum to the measured intensity spectrum by minimizing . We use the downhill simplex method [34] to perform the multidimensional minimization. The factor 2 in (15) is to account for both positive and negative frequencies, and the frequency difference in the denominator is a normalization factor. The integration over frequency spans from to . We chose to be the smallest nonzero frequency after averaging the spectrum according to Welch’s method [20], which corresponds to 0.0167 Hz as explained in Section 2. The frequency is chosen as the maximum useable frequency before the noise floor in the PSD is reached. We estimate this frequency from the measured intensity spectrum by identifying frequency samples for which dB, and then taking as the frequency that separates the first 5% of these samples from the remaining 95%. The −35 dB cutoff is receiver specific and was determined by numerical experimentation.
We note that IPE analysis, as described above, is a computationally intensive procedure. The calculations shown in this paper required more than three days (wall clock time) to produce on a Pentium-class Quad Core workstation. A faster approach is possible, however, where the forward propagation calculation described in Section 3.2 is replaced by an approximate calculation, as follows. First, numerical realizations of the screen phase are generated which are characterized by the autocorrelation function given in (2). Next, a plane wave is propagated through the screen and through free space down to the ground (using, for example, the phase screen technique described in [8, 35] or [24]). Finally, this process is repeated many times and the results ensemble averaged to produce a smooth intensity spectrum, which is needed for least-squares fitting the measured intensity spectrum.
4. Results and Discussion
4.1. Example Application of the IPE Technique
The methodology described in Section 3 is a technique for inferring the parameters of the ionospheric screen from a time series of intensity fluctuations due to transionospheric propagation in the equatorial ionosphere. Figure 3 shows two example applications of the IPE technique using GPS measurements collected at Campo Grande on 1-2 November 2002. The time series of intensity fluctuations in Figure 3(a) correspond to a case of relatively weak scatter measured late in the evening (26:00 UT). The C/No fluctuates a few dB about the background level and the measured scintillation intensity index, , is relatively small. The rate of intensity fluctuations, quantified by the decorrelation time, = 1.2 sec, is moderate. Accounting for the propagation geometry, magnetic field configuration, and velocity of the GPS satellites at the time the measurement was taken, as described in Section 3, the IPE technique gives the ionospheric screen parameters as , , and m/s. The corresponding effective scan velocity in this case was m/s. The for the spectral fit was 0.052. Both the model intensity spectrum (Figure 3(b) red curve) and the measured spectrum (Figure 3(b) black curve) show evidence of Fresnel rings, which occur in weak scatter when the ionospheric scattering layer is relatively thin. Since the scatter is weak and the phase spectral index is relatively shallow, frequencies smaller than the Fresnel break frequency, , are suppressed by diffraction during free-space propagation [28]. For frequencies larger than , the intensity spectrum follows a power law with slope equal to the spectral index of the screen. The time series of intensity fluctuations in Figure 3(c) correspond to a case of strong scatter encountered earlier in the evening (23:34 UT). In this case, the C/No undergoes deep fades of up to 20 dB, and the scintillation index is , which is close to the saturation value of 1.0. The fades occur more rapidly with a decorrelation time of sec, which is partly due to multiple scatter effects which generate small scale features at the reception plane [28, 36], and also partly due to a faster effective scan velocity ( m/s) in this case [5]. Frequencies smaller than the Fresnel break scale are suppressed less effectively than in the weak scatter case due to the influence of refractive scatter [28]. For frequencies larger than but less than 1/, the intensity spectrum deviates substantially from a power law by broadening and steepening under the influence of refractive scatter. For frequencies larger than about 1/, the power law behavior of the intensity spectrum is restored, but at these frequencies the spectral power is often smaller than the noise floor of the receiver. In this case, it is not possible to measure the spectral index of the ionospheric screen directly from the slope of the intensity fluctuations when the scatter is strong. The IPE technique provides a way to circumvent this problem; the phase spectral index (and therefore also the irregularity spectral index) can be retrieved even when the measured intensity spectrum does not manifest a power law regime.
(a)
(b)
(c)
(d)
4.2. Directly Measured Parameters
We begin by presenting the latitudinal and local time morphology of those parameters which are measured directly from the GPS observations (i.e., without IPE analysis). Figure 4 shows the variation of vertical-equivalent total electron content (TEC), the scintillation index , and the decorrelation time (). The latitudes and local times for which the measurements have been assigned correspond to the latitudes and local times of the 350 km ionospheric penetration points. The total electron content was computed from the GPS observations using the two-frequency technique as described in [37]. For the TEC, all available data samples (with 60 second cadence) are shown in Figure 4(a). For the scintillation index and decorrelation time, the plots (Figures 4(b) and 4(c)) show data samples only for the statistically stationary segments of the intensity time series for which , as described in Section 2. Because of this threshold, fewer samples are shown for the equatorial station Alta Floresta than for the anomaly stations Boa Vista and Campo Grande, since the scintillations at Alta Floresta were much weaker.
(a)
(b)
(c)
We observe, from Figure 4(a), that a well-developed equatorial anomaly was present on this evening (1-2 November 2002). The northern crest of the anomaly was located a few degrees northward of Boa Vista. The southern crest of the anomaly was located a few degrees southward of Campo Grande. Multiple gaps in the TEC data are evident prior to local midnight when the scintillation activity caused loss of lock on either the GPS L1 or L2 carriers (or both). These TEC gaps preclude a direct quantitative correlation between TEC and during strong scintillations, but it is clear from Figure 4 that the largest values generally occur at or near the crests of the equatorial anomaly and during early postsunset local times. The decorrelation time tends to be short ( sec) when the scintillation is intense (i.e., larger than 0.5-0.6), but there can be exceptions to this rule. In general, the decorrelation time depends on both the strength of scatter and also the effective scan velocity [5, 36].
At this point we can address the validity of the IPE approach, as applied to this specific dataset. Booker et al. [27] state that the approach we have used to solve the 4th moment equation is limited to small-angle scatter such that , where is the signal wavelength (0.19 m for the GPS L1 signal) and is the correlation length of intensity fluctuations on the ground. The correlation length can be computed from the effective scan velocity and decorrelation time as . Since we have computed both and τ for every case in which the IPE technique has been applied, we can calculate the correlation lengths explicitly. The minimum correlation length we encountered in this dataset was 40 m. This corresponds to a maximum scattering angle of 8 × 10−4 radians, which is much smaller than 1, and therefore our propagation calculations using the IPE are valid even for the most strongly disturbed conditions that were observed during the experiment.
4.3. Parameters Inferred by IPE Analysis
As described in Section 3, application of the IPE technique provides the parameters of the ionospheric screen that result in a model intensity spectrum that best matches the measured intensity spectrum in a least-squares sense. Figure 5 shows the χ2 values from these least-squares fits for all applications of the IPE technique presented in this paper (from all three stations combined). We chose to discard the IPE results if the χ2 of the fit exceeded a threshold of 0.5, considering these as poor fits between the model and the measurements. As can be seen from the histogram, however, the number of IPE results rejected on this basis was very small. Figure 6 shows a comparison between the scintillation index and decorrelation time calculated from the model intensity spectrum with those calculated directly from the measurements. The close agreement between the measured and modeled scintillation index and decorrelation time indicate the accuracy of the least-squares fitting, and the ability of the ionospheric turbulence model to reproduce the observations.
(a)
(b)
Figure 7 shows the ionospheric screen parameters inferred from the IPE analysis of the GPS observations on this evening (1-2 November 2002). The calculated values for the vertically integrated turbulent strength () span roughly 2 decades, ranging from 1 × 1034 to 1 × 1036. Comparing Figures 7(a) and 4(a) reveals that the largest values of were generally encountered close to the crests of the anomaly, northward of Boa Vista and southward of Campo Grande. The largest values of CkL were observed by the receiver at Campo Grande, and the scintillations persisted 2 hours longer at Campo Grande than at Boa Vista. This is consistent with the results reported by de Paula et al. [14], which suggests the scintillations may be more intense at Campo Grande than at Boa Vista due to a combination of the effects of a larger postsunset vertical plasma drift and delayed collapse in peak electron density (over Campo Grande). Since the geometrical aspects of the scattering are not reflected in (it is a properly of the random medium alone), it is conjectured that this parameter may correlate better with the vertical TEC in the region than would the scintillation index , which depends on both the medium and the propagation geometry. We were not able to confirm this hypothesis in a quantitative manner using the GPS observations alone, due to the numerous loss of lock events caused by the scintillation itself (which prevented us from measuring the TEC). The difficulties associated with TEC estimation using satellite links that are strongly scintillating are widely known, for example, [38].
(a)
(b)
(c)
The values for the phase spectral index of the ionospheric screen (p), as inferred from the IPE analysis, generally ranged from 2.5 to 4.5. No significant differences in the phase spectral index were observed between the three COPEX stations. Figure 7(b) shows a clear increase in the spectral index with increasing local time, which the authors believe is a new result. The increase of with increasing local time may suggest the erosion (decay) of small scale features in the turbulence as it evolves in time. On the other hand, the values of (Figure 7(a)) or (Figure 4(b)) on this evening do not show much evidence of decay with increasing local time, except that after 25 LT there are no longer any data samples which meet the requirements of (and stationarity). The observation that the spectral index is steepening while the scintillation strength is unchanging may suggest that while the small scale features are eroding away, the large scales features near the Fresnel scale (on the order of 400 m) which contribute most strongly to the scintillations at L band, have been only modestly affected.
Figure 7(b) shows the irregularity zonal drift velocities () inferred from the IPE analysis. The irregularity zonal drift velocities were roughly the same at all three COPEX stations, generally ranging from 100 to 200 m/s, with the exception of those values computed for early local times. Prior to about 20:30 LT, the estimates of were anomalously small at Boa Vista (less than 100 m/s) and anomalously large at Campo Grande (larger than 200 m/s). At early local times the equatorial plasma bubbles that generate the plasma turbulence are still evolving and have a vertical component of velocity in addition to a zonal component. As pointed out by Rino et al. [16–18], the ionospheric turbulence model described in Section 2 is valid only for fully developed irregularities that do not evolve as they translate. Furthermore, we have assumed the irregularity motion to be purely zonal, which is clearly incorrect at early local times. Later in this paper we will show that after 20:30 LT the IPE inferred zonal irregularity drifts compare favorably with measurements of the zonal drift made using the spaced-receiver technique. It is interesting to note that it is possible, using the IPE technique, to infer the zonal drift using only a single receiver. This capability could provide new opportunities to measure the irregularity zonal drift velocity at sites where only a single receiver is available.
4.4. Characterization of Phase Scintillations
The parameters , , and , along with the assumed values we have used for the altitude , outer scale , and anisotropy ratio (), completely specify the ionospheric screen and its translational velocity. For the sake of completeness, we also show the latitudinal and local time variation of two additional parameters which characterize phase scintillations and can be calculated from these others, namely, the strength of the spectrum of phase fluctuations and absolute phase variance. Rino [16] shows that when the ionospheric screen is described as in Section 2, and if diffraction effects on the phase are neglected, the 1D temporal spectrum of phase fluctuations can be expressed as where is the strength of the temporal spectrum of phase at a scale of 1 Hz. The parameter can be expressed in terms of other parameters of the screen and the geometry which are known after the IPE analysis: Integrating over all frequencies gives an expression for the absolute phase variance: Note that the phase variance depends on the outer scale wavenumber, which we do not measure. Hence the values for the phase variance we report, like , must be considered as relative to our specific choice of outer scale (10 km). It is also important to note that the expressions (16) and (17) were derived while neglecting the effects of diffraction on the measured phase, which become important when the scatter is strong [8]. Nevertheless, the parameters and have been used extensively in the literature for characterizing the effect of scintillation on GPS tracking loop performance [39, 40] so we report these findings here.
Figure 8 shows the latitudinal and local time variations of the parameters and during the evening of 1-2 November 2002. The strength of the phase spectrum (T) ranged from −60 dB (weak) to −15 dB (very strong), with the largest values encountered at early local times at the anomaly stations Boa Vista and Campo Grande. The absolute phase variance reached up to 500 radians squared, which is well into the strong multiple-scatter regime [28]. The largest values of the absolute phase variance occurred near the anomaly crests, northward of Boa Vista and southward of Campo Grande. We note that using the IPE technique, it is possible to estimate the parameters and even when the GPS phase observations cannot be analyzed directly due to loss of phase lock. This may be helpful when studying the influence of scintillation on loss of lock and GPS positioning accuracy, where one of the principal challenges is to characterize the conditions under which loss of phase lock occurs in spite of the fact that loss of phase lock precludes direct measurement of and [5].
(a)
(b)
4.5. Local Time Variation at the Three Stations
The plots shown in Figures 4, 7, and 8 allow examination of the parameters as a function of latitude and local time. This format has been chosen in order to compare the parameters to the structure of the equatorial anomaly, which varies with latitude and local time. Nevertheless, the local time dependence of the screen parameters is perhaps more clearly shown in Figure 9. The plots shown in Figure 9 include all data, irrespective of the latitude at which they were measured. Figure 9(a) shows the dependence of with local time. The largest values of were encountered at the anomaly stations Boa Vista and Campo Grande. A clear latitudinal asymmetry in the turbulent intensity is evident in Figure 9(a), as the values were generally larger at Campo Grande than they were at Boa Vista. There was little change in with local time during this particular evening, but the scintillations exceeding our threshold of persisted 2 hours longer at Campo Grande than at Boa Vista. As commented earlier, the phase spectral index (Figure 9(b)) showed a marked change from 2.5 just after local sunset increasing to as large as 4.5 after local midnight. Again, this steepening of the spectrum may suggest that the small scale features of the turbulence are being eroded away. Figure 9(c) reveals that the strength of the phase fluctuation spectrum decreases with increasing local time. From (17), it is easy to show that if and all other parameters except are held constant, then will decrease approximately linearly with increasing on a logarithmic scale. Under these conditions the strength of phase scintillations decrease with increasing local time. We note that it would be very useful, for both academic interest and also modeling purposes, to know how the strength and spectral index of the plasma turbulence within the same equatorial plasma bubble change as a function of time. Since the measurements from the COPEX campaign are from receivers located along a fixed magnetic meridian, we cannot follow the evolution of a single bubble as it drifts zonally. It would be interesting to repeat this experiment using a chain of receivers distributed zonally to infer how the screen parameters vary in a frame of reference that follows the bubbles as they drift.
(a)
(b)
(c)
4.6. Comparison with Spaced-Receiver Measurements of the Irregularity Zonal Drift
The IPE analysis provides estimates for the parameters of the ionospheric screen which yield a temporal spectrum of intensity fluctuations that most closely match the observations in a least-squares sense. One of these screen parameters is its translational velocity in the zonal direction, . This is a model-inferred irregularity drift velocity, not a direct measurement of the rate at which the density irregularities sweep past two points in space as is provided by the spaced-receiver technique. Both the spaced-receiver technique and also IPE analysis are subject to errors if the irregularities evolve temporally as they drift. An excellent discussion of this issue and other potential sources of error when measuring the zonal irregularity drift using spaced GPS receivers may be found in [19]. A detailed analysis of zonal irregularity drift velocity using the spaced GPS receiver technique during the COPEX campaign is provided by Muella et al. [11] and de Paula et al. [14].
In this section, we compare the zonal drift estimates provided by the IPE technique with the estimates provided by the colocated 4-channel VHF receivers at Boa Vista, Alta Floresta, and Campo Grande, and also with the GPS spaced-receiver measurements provided by Muella et al., as described in [11]. The estimates from Muella were measured at Boa Vista, Cachimbo, and Campo Grande. We will compare our measurements at Alta Floresta (9.9°S, 56.1°W) to Muella’s measurements at Cachimbo (9.5°S, 54.8°W), since these stations are close to one another. For both the GPS spaced-receiver drift estimates and the IPE drift estimates, the geometry of propagation with respect to the magnetic field has been taken into account, and also the motion of the ionospheric penetration points. For the VHF spaced-receiver measurements, these geometrical considerations have not been accounted for, but these corrections are smaller when the satellites are viewed at high elevation angles. Therefore we only show results for the east links which are visible at higher elevation angles than the west links.
Figure 10 shows a comparison of the zonal drift velocity estimated using the IPE technique with estimates provided by the VHF east links (channels 3-4), and the GPS spaced-receiver measurements provided by Muella et al. [11]. We observe that after approximately 20:30 LT, the IPE estimates compare favorably with both the VHF spaced-receiver estimates and the GPS spaced-receiver estimates. All three techniques show considerable scatter at early local times, presumably because the plasma bubbles are still evolving and have significant vertical components of velocity.
(a)
(b)
(c)
The fact that the IPE inferred drifts agree with the VHF and GPS spaced-receiver drift estimates is encouraging and may lend confidence to the IPE estimates of and , parameters for which no independent source of validation data is readily available. We plan to repeat the IPE analysis for the remaining evenings of the COPEX campaign in order to validate the technique and also to study the average seasonal behavior of the ionospheric screen parameters during the duration of the experiment.
We should emphasize that the IPE technique is not intended as a replacement for the spaced-receiver technique. The spaced-receiver technique provides a direct measurement of the zonal irregularity drift once the necessary geometric corrections [19] have been applied, whereas IPE analysis provides a model-inferred drift estimate. For example, the IPE drift estimates would change if a different model for the irregularities were used (e.g., a two-component spectral model). As such, the IPE technique cannot be expected to provide as reliable estimates of the zonal drift as the spaced-receiver technique. Furthermore, the nonlinear nature of the inversion process makes it unclear how to provide error bounds for these drift estimates. Despite these limitations, however, the IPE technique should be useful for ground stations where zonal drift estimates are desired but only a single receiver is available.
5. Remarks and Conclusions
In this report, we introduce a new technique called Iterative Parameter Estimation (IPE) for inferring ionospheric turbulence parameters from a time series of scintillating intensity measurements resulting from propagation through the low latitude ionosphere. The method of analysis is new and, unlike analytical theories of scintillation based on the Born or Rytov approximations, it is valid for strong scatter when the scintillation index () saturates due to multiple-scatter effects. IPE employs the phase screen approximation and a numerical inversion technique to estimate the parameters of an ionospheric screen which are consistent with the ionospheric turbulence model employed and the scintillations observed. More specifically, IPE provides a set of screen parameters that produces a model intensity spectrum which best matches the measured intensity spectrum in a least-squares sense. Due to the analytic intractability of the 4th moment equation governing fluctuations in the field we must determine the optimal parameters by numerical iteration, and therefore we cannot guarantee the uniqueness of this parametrization.
The underlying ionospheric turbulence model we have used to fit the scintillation observations described in this paper is that developed by Rino [16]. This model assumes the electron density irregularities can be characterized by a slab of homogenous, anisotropic fluctuations with a single-slope power law spectral density function in three dimensions. Actual turbulence in the equatorial ionosphere, by contrast, consists of large-scale deterministic structure (equatorial plasma bubbles) with small-scale random structure embedded within these bubbles [17, 24]. Furthermore, there is evidence to suggest that a model based on a two-component power law is better able to characterize the plasma turbulence under some circumstances [41]. Despite these limitations, the Wideband Scintillation Model which is based on this formulation has been shown to produce a very satisfactory description of the scintillation of transionospheric signals at frequencies ranging from VHF to L band (e.g., see [22] and the references therein). The IPE technique itself is not constrained to use any particular turbulence model, however, and it would certainly be possible to extend the approach to accommodate a two-component spectrum, for example.
We use the IPE technique to investigate the latitudinal and local time variation of the plasma turbulence parameters from 10 Hz GPS C/No measurements during one evening (1-2 November 2002) of the Conjugate Point Equatorial Experiment (COPEX) in Brazil. We find that the strength of turbulence tends to be largest near the crests of the equatorial anomaly and at early postsunset local times. We note that the results of Muella et al. [42] suggest that the most intense scintillations do not occur exactly at the crests of the anomaly, but instead at the edges of the crests where TEC gradients are largest. Due to the limited number of observations considered in our analysis, and the gaps in our TEC estimates due to scintillation induced loss of lock, we can neither corroborate nor refute this claim. During the evening we considered, the strength of turbulence showed little variation as function of local time but was generally stronger and lasted two hours longer at Campo Grande than at Boa Vista. The phase spectral index was similar at the three stations but increased significantly from 2.5 to 4.5 with increasing local time. We hypothesize that this change in spectral slope may have occurred as small scale features of the turbulence were gradually eroded. The larger-scale sizes in the turbulence (near the Fresnel break scale), which contributed most to the intensity scintillations, may have decayed more slowly so that the turbulent intensity decreased more gradually. An alternative interpretation for the local time dependence of the spectral index is that the physical mechanisms which generated these irregularities may have differed depending on the time at which the bubbles were generated. With GPS receivers operating from fixed locations on the ground, it is not possible to distinguish between old bubbles which have been decaying for a long period of time and new bubbles which have developed at later local times. Additional work needs to be performed in order to resolve this issue. It is hoped that either in situ density observations or perhaps a longitudinally distributed chain of GPS receivers may enable us to determine how the spectral index varies as a function of the local time at which the bubbles developed. Even though we examined only the GPS intensity fluctuations and not the GPS phase measurements, we were able to infer the strength of phase fluctuations and the absolute phase variance based on theoretical considerations. From this analysis, we inferred that the absolute phase variance reached up to 500 rad2, which is well into the strong multiple-scatter regime. The strength of phase scintillations was inferred to decrease with increasing local time. We found that estimates of zonal irregularity drift using the IPE technique agreed favorably with those made using the VHF and GPS spaced-receiver techniques. This encouraging result suggests that IPE analysis may provide useful estimates of zonal irregularity drift at scintillation monitoring sites which are equipped with only a single receiver, so that the spaced-receiver technique cannot be employed.
Since the COPEX campaign was conducted back in 2002, the number of GPS receivers capable of providing high rate (10 Hz or faster) scintillation observations in South America has increased dramatically. The Low-Latitude Ionosphere Sensor Network (LISN) [43, 44] includes more than 35 such receivers and also will include five ionosondes distributed along a magnetic field line in a fashion similar to COPEX. The application of IPE analysis to GPS scintillation observations collected along this field line could enable routine and systematic investigation of the latitudinal and local time morphology of ionospheric turbulence that produces radio wave scintillations at low-latitudes. It may also be instructive to use a longitudinally distributed chain of GPS receivers to follow the evolution of individual plasma bubbles to determine how the turbulence evolves in a reference frame that follows the bubbles. Such an investigation might help better explain why the spectral slope was observed to increase significantly with local time.
Acknowledgments
The authors would like to thank Robert Livingston for providing the VHF spaced-receiver irregularity drift estimates and Marcio Muella for providing the GPS spaced-receiver irregularity drift estimates during the COPEX experiment in Brazil. This work was supported by AFRL contract FA8718-09-C0041 and also funding from AFOSR.