Characterization of ionospheric amplitude scintillations using wavelet entropy detrended GNSS data

Available online at www.sciencedirect.com

ScienceDirect Advances in Space Research 54 (2014) 2172–2183 www.elsevier.com/locate/asr

Characterization of ionospheric amplitude scintillations using wavelet entropy detrended GNSS data Yongqing Su, Hao Liu ⇑, Jiguang Yue, Yunfan Yang School of Electronics and Information Engineering, Tongji University, 4800 Caoan Road, Shanghai 201804, China Received 13 March 2014; received in revised form 11 August 2014; accepted 13 August 2014 Available online 23 August 2014

Abstract The extensive monitoring networks of Global Navigation Satellite System (GNSS) ionospheric scintillation have been established to continuously log observation data. Further, the amplitude scintillation index and the phase scintillation index, which are derived from scintillation observations, are anticipated to accommodate the accuracy requirement of both the user level and the monitoring station level. However, raw scintillation observations essentially measure superposed waveform impairments of GNSS signals propagating through ionosphere and troposphere. It implies that fluctuations of raw scintillation observations are caused by multiple factors from the entire radio propagation environment. Hence, it is crucial to characterize ionospheric scintillations from GNSS observation data. And the characterization is implemented through extracting fluctuations of raw observations merely induced by ionospheric scintillations. Designed to address this problem by means of Fourier filtering detrending, the present work investigates the influence of varying detrending cutoff frequencies on wavelet statistical energy and wavelet entropy distributions of scintillation data. It consequently derives criteria on the optimum detrending cutoff frequency for three types of raw amplitude scintillation data, which are classified by their wavelet energy distributions. Results of the present work verify that detrending with specific optimum cutoff frequencies rather than the fixed and universally applicable one renders the validity and credibility of characterizing ionospheric scintillations as the part of GNSS observation fluctuations purely induced by ionosphere electron density irregularities whose scale sizes are comparable with or smaller than the Fresnel scale. Ó 2014 COSPAR. Published by Elsevier Ltd. All rights reserved.

Keywords: Ionospheric scintillation monitoring; Optimum detrending; Wavelet analysis; Signal entropy

1. Introduction Radio signals of satellite-borne systems, such as GNSS and SAR, are impaired by ionosphere electron density irregularities of spatial and temporal random variations when propagating in the transionospheric links. Hence, signal waveforms in the plane of phase center of satelliteborne or ground-based receiver antennas present stochastically fluctuating amplitude, phase, angle of arrival etc., ⇑ Corresponding author.

E-mail addresses: [email protected] (Y. Su), 11liuhao@tongji. edu.cn (H. Liu), [email protected] (J. Yue), [email protected] (Y. Yang). http://dx.doi.org/10.1016/j.asr.2014.08.012 0273-1177/Ó 2014 COSPAR. Published by Elsevier Ltd. All rights reserved.

which are referred to as ionospheric scintillations. Severe ionospheric scintillations damage essential functions of satellite systems depending on transionospheric telecommunications. Even in moderate cases, stochastic fluctuations inevitably degrade the integrity, accuracy and validity of high-precision-oriented positioning and time service of GNSS to a significant degree. Ionospheric scintillations predominantly outbreak in the equatorial ionization anomaly (EIA), polar cap and auroral oval regions, rarely in mid-latitudes (Aarons, 1982; Forte et al., 2002; Kintner et al., 2007). Significant influences and widespread geographical distributions of this space weather phenomenon promote locally and globally successive investigations and establishments of monitoring

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

stations (Pi et al., 1997) in decades. From traditional in situ satellites to modern GNSS monitoring networks (Basu et al., 2002), ionospheric scintillation observations are obtained with wide spatial observation distributions and at high sample rates. Such precise observations contribute to facilitating the evolution analysis of ionospheric irregularities which cause scintillations, validating scintillation models, such as WBMOD (Secan et al., 1995) and its upgrade SCINTMOD (Secan et al., 1997), and providing crucial inputs of scintillation impact assessment models for distinguishing performance levels with regard to civil aviation, positioning and timing. Ionospheric scintillation indices, such as the amplitude scintillation index S 4 and the phase scintillation index r/ , which are derived from GNSS scintillation observations, measure random amplitude and phase fluctuations of GNSS signals respectively. Thus both indices could be roughly regarded as quantitative assessments of GNSS performance impacts under the condition of ionospheric scintillations. Further, in terms of the weak scatter theory (Rino, 1976, 1982), scintillation indices contribute to connecting observations with scintillation propagation models (Beach and Kintner, 1999). However, perils of scintillation indices (Beach, 2006; Mushini et al., 2012) lie in not only the discrepancy between observation-based indices and indices output by theoretical models, but also the overlooked fact that both observations and indices essentially measure entire fluctuation owing to ionospheric scintillations and other multi-factor patterns, including multipath, intersatellite interference (Beach and Baragona, 2007) and trends. All of the disturbance patterns superpose in the integrated GNSS link and they are indicated as the combined transmission impairment as a result of raw scintillation observations. For the purpose of accurate scintillation data and indices, it is fundamental to characterize ionospheric scintillations through extracting observation fluctuations merely induced by scintillations. With a fine antenna pattern located in the multipath-free and interference-free environment, it is proper to take no account of multipath and interference. However, trends are mainly caused by slow variations of background ionosphere and satellite motion, significantly superimpose extra fluctuations on those of ionospheric scintillations, and consequently lead to erroneous ionospheric scintillation indices. Hence, eliminating trends of GNSS ionospheric scintillation data, or referred to as detrending, has been continuously studied for years. The relevant investigation has gone through two stages. The initial detrending (Van Dierendonck et al., 1993; Van Dierendonck and Arbesser-Rastburg, 2004) adopted the filtering framework with the fixed and universally applicable cutoff frequency. The cutoff frequency was determined by abundant low-latitude observation data from in situ DNA Wideband satellites (Fremouw et al., 1978). Afterwards, the “phase without amplitude scintillation” phenomenon (Doherty et al., 2004) was discovered and drew attention to some research, including interpreting GNSS ionospheric scintil-

2173

lation index distortions, defining trends of observation fluctuations and modifying detrending methods. Subsequently, investigations into the phenomenon are concentrated on the detrending filtering cutoff frequency (Forte et al., 2002), the geometrical control of ionospheric scintillation indices (Forte and Radicella, 2004; Forte, 2007) and the weak scatter approximation (Doherty et al., 2004). Finally, there is a consensus on the interpretation of the phenomenon. It is indiscriminately applying the fix and universally applicable cutoff frequency originally fitting for observations of low latitudes to those of high latitudes that induces abnormal correlations between amplitude and phase scintillation indices of high latitude scintillation events. The inherent correlation is essentially dominated by the definition of GNSS ionospheric scintillation. To take the Fresnel scale as a reference, ionospheric irregularities induce radio observation fluctuations subject to diverse propagation mechanisms in terms of scales (Knepp and Wittwer, 1984). Diffraction or angle scattering effects owing to smaller scale irregularities dominate variations of amplitude fluctuations, or referred to as amplitude scintillations. Whereas larger scale irregularities, whose equivalent focal lengths are larger than the propagation distance of radio signals, contribute to significant phase fluctuations with negligible amplitude fluctuations. Hence, only phase fluctuations owing to smaller scale irregularities, or referred to as phase scintillations, correlate with amplitude scintillations. Comparisons of power spectral density between amplitude and phase scintillations demonstrate that detrending with the fixed and universally applicable cutoff frequency leads to an overestimate of phase scintillations in the polar region (Forte et al., 2002). To address this situation, an appropriately enlarged cutoff frequency eliminates the residual trend and dramatically weakens fluctuations of phase observations owing to larger scale irregularities, which represent slowly varying trends (Forte, 2005; Forte et al., 2011). The determination of a proper filtering cutoff frequency is crucial for detrending raw scintillation data. For amplitude scintillation data, an adaptive determination of the optimum cutoff frequency is presented in Materassi and Mitchell (2007). Moreover, running mean and polynomial methods are utilized to fit trends, and even time length of raw scintillation data is used to determine the cutoff frequency (Materassi et al., 2009). It is their respective equivalent cutoff frequencies that distinguish these techniques. Then the wavelet-filtering-based detrending presented in Mushini et al. (2012), essentially utilizes wavelet transform reconstruction to detrend and denoise raw observations. However, its lower and upper cutoff scales, similar to the filtering cutoff frequency, are acquired by surveying scalograms of 400 ionospheric scintillation events in high latitudes. Thus, cutoff scales in that case, as a priori information, still adhere to the fixed and universally applicable cutoff frequency approach, which probably leads to ionospheric scintillation index distortions for raw observations in other latitudes.

2174

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

Regardless of its inherent under-detrending, the adaptive detrending emphasizes that determination of the optimum cutoff frequency depends on the assessment of detrended results with varying cutoff frequencies. In that case, the S 4 index is regarded as a traditional and convenient assessment approach. However, shortening time window length of S 4 calculation generates nondistinctive assessments for different detrending methods (Materassi et al., 2009). Since the first order statistic index intrinsically smoothes fluctuations within the time window, the assessment approach consequently induces a resolution limit, not to mention the implied stationarity assumption. Apart from scintillation indices, fading timescale (Kintner et al., 2001) and wavelet transform analyze scintillations in time domain and time–frequency domain respectively. As a more sophisticated assessment approach, scalograms obtained from wavelet transform, illustratively present detrended results without the stationarity assumption (Materassi and Mitchell, 2007). Note that the qualitative analysis approach is incapable of deriving the optimum cutoff frequency. The optimum detrending cutoff frequency delimits ionospheric scintillation components of raw observations. Detrending with the optimum cutoff frequency completely preserves scintillation components as much as possible and adequately eliminates trend components. Consequently, aiming at quantitative criteria on optimum cutoff frequencies for diverse raw observations, influences of varying cutoff frequencies on wavelet statistical energy and wavelet entropy distributions of scintillation data are investigated in this paper. The detrending approach with the optimum cutoff frequency highlights diversity of ionospheric scintillation observations and provides reliable data for ionospheric scintillation and ionospheric irregularities investigations. Section 2 describes the raw amplitude scintillation data and the detrending filtering framework. Section 3 analyzes effects of varying cutoff frequencies on wavelet statistical energy of scintillation data and derives the optimum cutoff frequency criteria in terms of wavelet entropy. Then optimum detrended results and the discussion are presented in Section 4. Finally, the conclusion is drawn in Section 5. 2. Scintillation data and detrending GNSS ionospheric scintillation observation data is provided by the Canadian high Arctic Ionospheric Network (CHAIN) (Jayachandran et al., 2009). The ground-based ionospheric monitoring network consists of 10 stations distributed from the auroral oval to the polar cap regions. GSV4004B ionospheric scintillation and TEC monitors (GISTM) (Van Dierendonck et al., 1993; Van Dierendonck and Arbesser-Rastburg, 2004), part of CHAIN instruments, record ionospheric scintillation indices and output raw power time series of GPS L1 signals at 50 Hz sample rate, which matches the data message rate of GPS signals. The time series measuring GNSS signal amplitude

fluctuations is a function of in-phase and quadrature signal components output by the channel correlator. It is a proper approximation of carrier to noise ratio (Beach and Kintner, 2001; Forte, 2005). The raw power time series (Van Dierendonck and Arbesser-Rastburg, 2004) is denoted as: xðtÞ ¼ F ðNBP  WBP Þ

ð1Þ

F ; the scaling factor, becomes meaningless due to being removed by subsequent detrending and S 4 index calculation. NBP , the narrow band power, is an estimator of signal strength, and it is defined as: !2 !2 N N X X NBP ¼ In þ Qn ð2Þ n¼1

n¼1

WBP ; the wide band power, another signal strength estimator, which is defined as: WBP ¼

N X

ðI 2n þ Q2n Þ

ð3Þ

n¼1

I n and Qn are respectively in-phase and quadrature signal component samples. Here N ¼ 20 is chosen to match the number of C/A code periods contained in one GPS navigation message data bit. This paper presents three raw signal power time series from CHAIN and they were observed in the year 2013. The lengths of them are identically half an hour. The station involved in this work is illustrated in Table 1 and details of these observations are described as follows: Time series I was observed on the link from the satellite PRN28 starting at 6:00 UT on March 20, 2013 (see Fig. 1). Time series II was observed on the link from the satellite PRN17 starting at 15:20 UT on November 4, 2013 (see Fig. 2). Time series III was observed on the link from the satellite PRN21 starting at 00:05 UT on September 19, 2013 (see Fig. 3). The raw amplitude scintillation data is converted into the canonical form of the decibel unit. These time series display diverse patterns of observation fluctuations. Time series I possesses several big bursts, and time series II shows a giant peak. As a comparison, time series III exhibits obvious trends with negligible scintillations. These features of interest are characterized as highlights of corresponding scalograms in Section 3. Note that antenna phase centers of CHAIN scintillation monitors are configured to operate in a fine line of sight and the observation data involved in this work was obtained above an elevation angle of 30°. The elevation angle mask is a conventional technique to mitigate multipath. It is based on the general observation that signals at low elevation angles are more likely to suffer from multipath. However, a previous study (Beach and Kintner, 2001) has indicated that multipath can occur at high elevation angles as well. It is generally caused by improper antenna installations and nearby buildings or structures extending to the field view. Amplitude fluctuation patterns owing to multipath repeat from day to day with the

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

2175

Table 1 Geographic and corrected geomagnetic latitude coordinates of the field of view of the CHAIN site at 350 km. Station

Geographic latitude

Geographic longitude

Corrected geomagnetic latitude

Corrected geomagnetic longitude

Cambridge Bay

69.12°N

105.03°W

76.49°N

46.50°W

Fig. 3. The amplitude scintillation event was observed from the satellite PRN21 by the receiver in the Cambridge Bay station between 00:05 UT and 00:35 UT on September 19, 2013. The plot is counted in minutes from the beginning of the observation at 00:05 UT and the raw power is given in the dB unit. The amplitude scintillation observation is referred to as time series III.

Fig. 1. The amplitude scintillation event was observed from the satellite PRN28 by the receiver in the Cambridge Bay station between 6:00 UT and 6:30 UT on March 20, 2013. The plot is counted in minutes from the beginning of the observation at 6:00 UT and the raw power is given in the dB unit. The amplitude scintillation observation is referred to as time series I.

Fig. 2. The amplitude scintillation event was observed from the satellite PRN17 by the receiver in the Cambridge Bay station between 15:20 UT and 15:50 UT on November 4, 2013. The plot is counted in minutes from the beginning of the observation at 15:20 UT and the raw power is given in the dB unit. The amplitude scintillation observation is referred to as time series II.

sidereal time shift in the invariant line of sight environment of antennas. Thus, for permanent monitoring stations, it is not proper to output observations contaminated by multipath at high elevation angles. The anomaly should be avoided by relocating stations and re-installing antennas. If the anomaly persists, it can be detected through a comparison among observations of successive days at the same time (with nearly a 4 min time shift each day). All of the raw time series involved in this work go through multipath

check. Hence, post-processing raw time series is of tolerance to neglect fluctuations due to multipath. The initial detrending procedure for Wideband satellite and GNSS ionospheric scintillation data begins with a low-pass filter deploying a specific cutoff frequency. Then filtering results, low frequency components of raw time series, are treated as trends. Finally, raw time series divided by trend time series are detrended time series (Fremouw et al., 1978; Van Dierendonck et al., 1993). An additional normalization prior to detrending renders signal power time series fluctuating around unity (Materassi and Mitchell, 2007). Since then Fourier-filtering-based method has become the basic framework of detrending GNSS ionospheric scintillation data. Other detrending approaches also depend on the exact definition of trends, such as running mean instead of low-pass filtering. A novel detrending framework should be given cautious consideration. The novel framework replacing division with subtraction (Forte and Radicella, 2002; Forte et al., 2011), satisfies constraints of phase screen theory on peak patterns of amplitude scintillation power spectral density. However, note that in the case of slowly varying time series, subtracting trends may induce approximate zero values, which lead to abnormally huge S 4 index values. Therefore, the division detrending is utilized in this paper, and the formulas are defined as: xN ðtÞ ¼

xðtÞ xðtÞ

xtrend ðt; f~ Þ ¼ ½xN ðtÞlpf ;f~ xdet ðt; f~ Þ ¼

xN ðtÞ xtrend ðt; f~ Þ

ð4Þ ð5Þ ð6Þ

xN ðtÞ; the normalized amplitude scintillation time series, is filtered by a Fourier low-pass filter to fit the trend time series xtrend ðt; f~ Þ. The detrended time series xdet ðt; f~ Þ is obtained by normalized time series dividing trend time series point by point. It is obvious that detrending with one specific

2176

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

filtering cutoff frequency obtains one specific trend time series and one detrended time series as well. To obtain the optimum detrended scintillation data, the optimum cutoff frequency is determined in the next section. 3. Methods The framework of Fourier-filtering-based detrending characterizes trends and scintillations of raw observations as low-frequency and high-frequency components of raw time series respectively. However, the optimum filtering cutoff frequency strictly delimiting them can only be derived in terms of specific raw time series. Since taking no account of diverse ionospheric scintillation observation patterns, detrending with the non-optimum cutoff frequency inevitably suffers from over-detrending or underdetrending. These two defective results respectively indicate excessively detrending with a too large cutoff frequency and inadequately detrending with a too small one. In this paper wavelet energy and wavelet entropy are utilized to evaluate detrended results and then derive optimum cutoff frequency criteria. 3.1. Wavelet statistical energy Wavelet transform is capable of analyzing locally slight variations in time–frequency domain of time series with fine frequency domain resolutions in low frequency band and fine time domain resolutions in high frequency band. Moreover, scalograms, the graphical representation of wavelet statistical energy distribution (Materassi and Mitchell, 2007; Materassi et al., 2009; Mushini et al., 2012), are used to analyze nonstationary time series consisting of multi-frequency components. This visual analysis tool highlights scintillation components therein.Continuous wavelet transform of detrended time series is defined as: Z 1 xðs; tÞ ¼ pffiffi w ðtÞxdet ðt; f~ Þdt ð7Þ s I s;s And Morlet wavelet selected as mother wavelet is defined as:   1 2 1=4 w0 ðtÞ ¼ p exp ix0 t  t ð8Þ 2 After time translation  and wavelet scale variation, w0 ðtÞ and wavelet function is denoted is transformed into w ts s as: t  s 1=2 ws;s ðtÞ ¼ jsj w ð9Þ s The superscript of ws;s ðtÞ in Eq. (7) denotes complex conjugation operation. Owing to Morlet wavelet employed here, wavelet transform involved in this paper belongs to nonorthogonal wavelet transform. The nonorthogonal wavelet transform induces redundance at large scales due to highly correlated wavelet spectrum at adjacent times.

And it is suitable for processing time series with smoothly and continuously varying wavelet amplitudes. However, the orthogonal wavelet transform is not utilized here because it induces a different wavelet spectrum due to the aperiodic shift of raw time series (Torrence and Compo, 1998). The definition of the scales in terms of fractional powers of two is denoted as: sn ¼ s0 2nj ;

n ¼ 0; 1; . . . ; J

1

J ¼ j log2 ðT =s0 Þ

ð10Þ ð11Þ

s0 is the smallest resolvable scale, and j determines the largest number of the resolvable scale. T is the time length of data, denoted as: T ¼ N  dt

ð12Þ

N is total sample point number, and dt is the sample interval of time series. Then, wavelet statistical energy eðs; tÞ ¼ jxðs; tÞj

2

ð13Þ

or referred to as wavelet power spectrum, characterizes wavelet energy distribution of time series at scales. To render the comparability of wavelet statistical energy among scales for a single or multiple time series, a normalization for wavelet function at each scale is necessary, which is defined as: t  s  s 1=2 t  s wNs;s ws;s ¼ ð14Þ s dt s And the conversion of wavelet function energy into unity is denoted as: 2 E½eðs; tÞ ¼ N  E½j^xdet ðt; f~ Þj 

ð15Þ

Here, E½ is the expectation operator and the Fourier transform of xdet ðt; f~ Þ is ^xdet ðt; f~ Þ. Therefore, Fourier transform of time series dominates its wavelet statistical energy distribution. And scalograms are obtained by plotting wavelet statistical energy distribution in decibels (depicted in Figs. 4–6). The structured “scale band” feature (Materassi et al., 2009) of scalograms distinguishes fluctuation patterns due to scintillation, trend and noise components of raw time series. Large scales corresponding to trends indicate smooth and slow fluctuation patterns. Note that trend components of the three time series occupy most of their wavelet statistical energy. A typical scintillation scale band (scales > 0.09 s and <4.7 s) of scalograms is indicated in the previous study (Mushini et al., 2012). In terms of the scintillation scale band, diverse scintillation patterns pertaining to the three time series are described as follows. The scalogram of time series I depicts a concentrated distribution of wavelet statistical energy in the scale range of 2.5–5.7 s, which indicates narrow scale band scintillations (see Fig. 4). On the contrary, the scalogram of time series II displays a feature that a main peak extends in the scale range of 0.1–10 s due to scintillation components of significant wide scale band (see Fig. 5). Since wavelet

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

2177

any appreciable scintillations in the scale range of interest (see Fig. 6). Beyond the scintillation scale band, all the scalograms possess more or less fibrous high-frequency components due to noise or interference. Nevertheless, the noise scale band is not as significant as that of scintillations in terms of wavelet energy distribution. And the aforementioned wavelet-reconstruction-based filtering could suppress the high-frequency components. Hence, the discrimination between scintillations and trends is primarily concerned here and dealing with the noise scale band is beyond the scope of this article. Fig. 4. Scalogram of time series I as shown in Fig. 1 counted from minute 0 to minute 30 at the beginning of 6:00 UT. The plot shows a narrow scintillation band arising from minute 0 to minute 30 in the wavelet scale range of 1–10 s.

3.2. Cutoff frequency effect The cutoff frequency of filtering-based detrending, increased from zero to the upper limit, has an influence on wavelet statistical energy distribution of raw time series, which is denoted as: 2 Z 1 N ~ ~ eðs; t; f Þ ¼ pffiffi ws;s ðtÞxdet ðt; f Þdt ð16Þ s I

Fig. 5. Scalogram of time series II as shown in Fig. 2 counted from minute 0 to minute 30 at the beginning of 15:20 UT. The plot shows a wide scintillation band arising from minute 0 to minute 13 in the wavelet scale range of 0.1–10 s.

Fig. 6. Scalogram of time series III as shown in Fig. 3 counted from minute 0 to minute 30 at the beginning of 00:05 UT. The plot shows negligible scintillation in the wavelet scale range of 1–10 s.

transform provides fine time–frequency analysis, scalogram highlights of scintillation components keep pace with scintillation patterns of raw time series in terms of time evolution. For time series I, the most violent burst in Fig. 1 corresponds to the scintillation components in Fig. 4 from the 13th to 18th minute. In the case of time series II, both peaks of Figs. 2 and 5 approximately appear at the 11th minute. As a comparison, the scalogram of time series III presents noticeable trend components without

As a result, low-frequency trend components of raw time series are eliminated more and more adequately with the gradually enlarged cutoff frequency. Owing to the normalization of wavelet statistical energy, elimination of trend components renders scintillation components, if any, occupying more wavelet energy at scales of interest. Namely, wavelet statistical energy quantitatively emphasizes the scintillation scale band because of detrending. To investigate the influence of varying cutoff frequencies on wavelet statistical energy scale distribution of scintillation time series, for a certain cutoff frequency f~ 0 ; wavelet energy marginal scale spectrum density is defined as: R eðs; t; f~ 0 Þdt ~ qðs; f 0 Þ ¼ R R ð17Þ eðs; t; f~ 0 Þdtds The quantity indicates normalized wavelet relative energy distribution. Namely, it measures wavelet energy distribution in terms of scales. Hence, wavelet energy marginal scale distribution of scintillation components is quantified as wavelet energy accumulated throughout the entire time evolution of detrended time series at a specific scale within the scintillation scale band. And the inherent normalization renders the comparability of scintillation component energy distributions with different cutoff frequencies. A previous comparison of wavelet energy between three different detrending methods for the same raw time series is presented in Materassi et al. (2009). Those methods are essentially distinguished by their equivalent cutoff frequencies. However, there exist consistent peaks of wavelet energy scale spectra, referred to as the total wavelet power in that paper, at the same scintillation scale. The occurrence of the consistent peaks demonstrates scintillation components dominate wavelet energy due to the elimination of trends. And the peak scale indicates

2178

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

where scintillation components of the detrended time series distribute the maximum wavelet energy, independent of detrending methods. To exactly delimit wavelet energy distributions of scintillation and trend components, the extreme conditions of under-detrending and over-detrending should be avoided. Detrended results with an excessively small low-pass filtering cutoff frequency retain evident trend components, which still dominate wavelet energy distribution. Hence, the expected peak occurs within trend scale band or none of peaks occurs in under-detrending case. On the contrary, detrending with an excessively large cutoff frequency not only rigorously eliminates trend components, but also impairs part of scintillation components. As a result, the peak within scintillation scale band degenerates into a cliff whose edge is towards large scales. Thus, the optimum detrending method presented in this paper adopts a varying cutoff frequency from zero to the upper limit to avoid overdetrending. The effects of varying cutoff frequencies on wavelet statistical energy distributions of the three time series are illustrated in Figs. 7–9, respectively. With regard to time series I and II, both scintillation and trend components are

Fig. 9. Variation of wavelet energy marginal scale spectrum density of time series III with varying cutoff frequencies from 0 to 0.5 Hz. The plot shows very slight wavelet energy distribution of scintillation components corresponding with the scalogram signature shown in Fig. 6.

evidently distinguished by varying cutoff frequencies in terms of wavelet energy scale spectrum where low-frequency trends are rapidly attenuated and high-frequency scintillations are gradually revealed in the direction of increasing the cutoff frequency. Nevertheless, the distinction between these two time series is wavelet energy distributions of their respective scintillation components. The concentrated scale distribution is shown in Fig. 7 for time series I and the wide scale distribution of time series II is depicted in Fig. 8. Note that both scale distribution patterns are consistent with their respective scalogram patterns of scintillation components (see Figs. 4 and 5). In addition, the evident trends and very slight scintillations of time series III are also demonstrated by wavelet energy scale distribution with varying cutoff frequencies (see Fig. 9). 3.3. Wavelet entropy and the optimum cutoff frequency

Fig. 7. Variation of wavelet energy marginal scale spectrum density of time series I with varying cutoff frequencies from 0 to 0.5 Hz. The plot shows concentrated wavelet energy distribution of scintillation components corresponding with the scalogram signature shown in Fig. 4.

Wavelet entropy (Rosso et al., 2001) combining wavelet statistical energy and Shannon entropy can be used to derive criteria for analyzing and comparing wavelet statistical energy distribution of raw time series with varying cutoff frequencies. It is defined as: X S WT ðf~ Þ ¼  qðsj ; f~ Þ  ln qðsj ; f~ Þ ð18Þ j

Fig. 8. Variation of wavelet energy marginal scale spectrum density of time series II with varying cutoff frequencies from 0 to 0.5 Hz. The plot shows wide wavelet energy distribution of scintillation components corresponding with the scalogram signature shown in Fig. 5.

Wavelet entropy evaluates the degree of order (disorder) of raw time series. In terms of wavelet statistical energy, a concentrated (wide) wavelet energy distribution at scales, which obtains a relatively small (large) wavelet entropy value, consequently demonstrates a high degree of order (disorder). And in Eq. (18), each term of the summation constituting the total wavelet entropy is the wavelet entropy at a specific scale for a specific cutoff frequency, which is denoted as: S WT ðsj ; f~ Þ ¼ qðsj ; f~ Þ  ln qðsj ; f~ Þ

ð19Þ

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

The quantity measures wavelet entropy scale distribution (see Figs. 10–12), which indicates contribution of trends and scintillations on wavelet entropy variation. Aiming at the optimum cutoff frequency, the effects of varying cutoff frequencies on wavelet entropy distribution are investigated to derive the optimization criteria. When varying cutoff frequencies are much smaller than the optimum value, trend components occupy most of the wavelet energy spectrum density and thus overwhelmingly dominate the variation of wavelet entropy scale distribution. When varying cutoff frequencies gradually approach the optimum value, the elimination of trend components is accompanied with rapidly decreasing wavelet energy spectrum density and wavelet entropy distributed at large scales. Meanwhile, increasing wavelet energy spectrum density of scintillation components gradually dominates wavelet entropy scale distribution. Then in the case of detrending with a cutoff frequency within or even in excess of the scintillation band, impaired scintillation components possess a weak and smooth wavelet energy distribution whose higher degree of disorder shapes the total wavelet entropy curve into a plateau (see Fig. 13). Both scintillations and trends are multiple scale components, and especially the discrete multi-scale feature of trends is demonstrated by isolated peaks in their wavelet energy and wavelet entropy distributions. Consequently the total wavelet entropy may fluctuate within the phase of trends (see Figs. 13–15). Hence, to derive the optimum cutoff frequency delimiting trends and scintillations, discrete multiscale features of trends should be taken into account. Based on the specific analysis of wavelet entropy with varying cutoff frequencies, this paper presents the optimum cutoff frequency criteria for the three raw time series respectively. For the concentrated wavelet energy distribution (see Fig. 7), trend and scintillation components dominate the variation of wavelet entropy (see Fig. 10) one after the other with varying cutoff frequencies. Scintillation components with higher degree of order replace the disordered trend components and render wavelet entropy continu-

2179

Fig. 11. Variation of wavelet entropy scale distribution of time series II with varying cutoff frequencies from 0 to 0.5 Hz. The plot shows wide wavelet entropy distribution of scintillation components corresponding with the wavelet energy distribution signature shown in Fig. 8.

Fig. 12. Variation of wavelet entropy scale distribution of time series III with varying cutoff frequencies from 0 to 0.5 Hz. The plot shows very slight wavelet entropy distribution of scintillation components corresponding with the wavelet energy distribution signature shown in Fig. 9.

Fig. 13. Variation of total wavelet entropy of time series I with varying cutoff frequencies from 0 to 0.5 Hz. The curve is derived from the wavelet entropy scale distribution shown in Fig. 10.

Fig. 10. Variation of wavelet entropy scale distribution of time series I with varying cutoff frequencies from 0 to 0.5 Hz. The plot shows concentrated wavelet entropy distribution of scintillation components corresponding with the wavelet energy distribution signature shown in Fig. 7.

ously decreasing towards a high degree of order. Then detrending with larger cutoff frequencies filters scintillation components. And wavelet entropy steeply increases and rapidly becomes flat, indicating a disordered post-processing result. Note that during the initial phase, wavelet entropy rapidly fluctuates. The fluctuation demonstrates the multi-scale feature of trends. As for the main decline phase of wavelet entropy curve, it is the result of a mix of residual trends and gradually impaired scintillations in

2180

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

Fig. 14. Variation of total wavelet entropy of time series II with varying cutoff frequencies from 0 to 0.5 Hz. The curve is derived from the wavelet entropy scale distribution shown in Fig. 11.

Fig. 15. Variation of total wavelet entropy of time series III with varying cutoff frequencies from 0 to 0.5 Hz. The curve is derived from the wavelet entropy scale distribution shown in Fig. 12.

response to the varying cutoff frequencies. Hence, there exist rough subsection points (see Fig. 13) indicating abrupt alternation of wavelet entropy distribution from trend components to scintillation components and undesirable elimination of large scale scintillation components (see Fig. 10). Addressing this type of raw time series, it is proper to determine the frequency corresponding to the starting point of the main decline section as the optimum detrending cutoff frequency. The total wavelet entropy curve for wide wavelet energy distribution of disordered scintillation components is shown in Fig. 14. In view of wavelet entropy scale contribution (see Fig. 11), the wavelet entropy maximum corresponds to the cutoff frequency with which trends are mainly eliminated and the part of scintillation components adjacent to trends begins to be filtered. As the concentrated wavelet energy distribution is analyzed above, the optimum cutoff frequency for time series II is determined simply as the starting point of the main decline section, namely the wavelet entropy maximum. Consequently, the criterion completely preserves scintillation components, especially those components adjacent to trends at scales, and adequately eliminates trends. Compared with previous two types of scintillation time series, time series III, as the most likely type of GNSS

ionospheric scintillation observations, presents barely detectable scintillations submerged in evident trends. Its total wavelet entropy is shown in Fig. 15. In this case, trend components dominate nearly all the variation of wavelet energy distribution (see Fig. 9), which equivalently evaluates the degree of order/disorder for trends irrespective of very slight scintillations. When wavelet entropy scale distribution of trend components varies due to the increasing cutoff frequency, various peaks at different scales (see Fig. 12) dominate the variation of wavelet entropy successively. And the main peak corresponds to upward steep of wavelet entropy (see Fig. 15). Then disordered wavelet energy distribution of trend components becomes smooth due to large cutoff frequencies, but still suppresses that of scintillation components. Not to mention that a relatively small cutoff frequency is incapable of highlighting weak scintillations from evident trends. Therefore, a compromise between adequate elimination of trends and complete preservation of scintillations is essential for the criterion. In this case, the optimum cutoff frequency is determined as the frequency corresponding to the starting point of steady decline in the total wavelet entropy curve, namely the wavelet entropy maximum. The optimum cutoff frequency criteria depend on features of wavelet energy distribution and wavelet entropy for the three types of raw scintillation time series and they take account of the variety and diversity of ionospheric scintillation observations. Compared with the fixed cutoff frequency, such as 0.1 Hz, 0.3 Hz and even 0.5 Hz, the optimum cutoff frequency is capable of eliminating trends adequately and preserving scintillations completely for each type of scintillation time series, which is denoted as: f o : S WT ðf o Þ P S WT ðf~ Þ;

f~ < f U

ð20Þ

It follows that the optimum cutoff frequency make the total wavelet entropy maximum for all cutoff frequencies smaller than the upper limit f U . The upper limit is determined as a large value to avoid an erroneous global maximum corresponding to an excessively large cutoff frequency. Note that for time series I, the optimum cutoff frequency corresponds to the local maximum in the anterior segment of the total wavelet entropy curve, and the optimum cutoff frequencies of the other two correspond to the global wavelet entropy maxima. In the case of time series I, determining the local maximum rather than the global maximum, or the higher degree of disorder owing to the elimination of all scintillation components, is subject to the upper limit. 4. Results and discussion According to the optimum detrending cutoff frequency criteria, the optimum cutoff frequencies for the three raw scintillation time series are respectively derived as: f o1 ¼ 0:0170 Hz

ð21Þ

f o2 ¼ 0:0545 Hz

ð22Þ

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

f o3 ¼ 0:1560 Hz

2181

ð23Þ

The detrended results of power time series according to the optimum cutoff frequencies are shown in Figs. 16–18. Compared with plots of raw time series (see Figs. 1–3), the successive bursts of time series I become more evident and similar to each other (see Fig. 16). The peak pattern of time series II appears unchanged only with the thorough removal of trends after the optimum detrending (see Fig. 17), as well as the case of time series III. Note that all trends of these raw time series are eliminated adequately with a relatively smaller cutoff frequency than the fixed one so as to completely preserve scintillation components. Scalograms of detrended power time series are shown in Figs. 19–21. In comparison with scalograms of raw time series (see Figs. 4–6), both concentrated and wide wavelet energy distributions of scintillation components are enhanced due to detrending (see Figs. 19 and 20). Adequate detrending renders the wavelet energy of time series III distributed within the scintillation scale band despite negligible scintillations (see Fig. 21). Although edge artifacts restraining is adopted in the process of wavelet transform, part of edge effects, or referred to as cone of influence, is left in full-scale scalograms due to limited sample length of time series. However, the inherent imperfection due to assumptive signal periodicity of Fourier transform during wavelet transform has no explicit influence on the determination of the optimum cutoff frequencies. Hence, corresponding to the total length of time series, complete scalograms containing the maximum scale 1800 s are shown in Figs. 19–21. The detrending methods for scintillation observations from GNSS links show some limitations on the significant deviation in scintillation statistical indices, which are related to scintillation patterns and geophysical conditions of observation stations. To fairly evaluate these detrending methods, comparing the variation of scintillation statistical indices according to different detrending methods is intuitive and easy to operate (Forte, 2005). Based on wavelet analysis of scintillation time series (Materassi and Mitchell, 2007), scalograms, as a qualitative evaluation

Fig. 17. Detrended time series II compared with the raw power as shown in Fig. 2. The plot is counted in minutes from the beginning of the observation at 15:20 UT and the detrended power fluctuates around 1 with its optimum cutoff frequency f o2 ¼ 0:0545 Hz.

Fig. 18. Detrended time series III compared with the raw power as shown in Fig. 3. The plot is counted in minutes from the beginning of the observation at 00:05 UT and the detrended power fluctuates around 1 with its optimum cutoff frequency f o3 ¼ 0:1560 Hz.

Fig. 19. Scalogram of detrended time series I compared with the raw scalogram as shown in Fig. 4. The plot is counted in minutes from the beginning of the observation at 6:00 UT with its optimum cutoff frequency f o1 ¼ 0:0170 Hz.

Fig. 16. Detrended time series I compared with the raw power as shown in Fig. 1. The plot is counted in minutes from the beginning of the observation at 6:00 UT and the detrended power fluctuates around 1 with its optimum cutoff frequency f o1 ¼ 0:0170 Hz.

tool of higher resolution, display detrended results precisely. However, scalogram results retain consistency with those of statistical indices. Then the frequency spectrum of the S 4 index is introduced to evaluate detrended results (Forte et al., 2011). A novel evaluation depending on correlation coefficient between amplitude and phase scintillation index is presented in Mushini et al. (2012). And the evaluation satisfies constraints of phase screen theory on

2182

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

Fig. 20. Scalogram of detrended time series II compared with the raw scalogram as shown in Fig. 5. The plot is counted in minutes from the beginning of the observation at 15:20 UT with its optimum cutoff frequency f o2 ¼ 0:0545 Hz.

Fig. 21. Scalogram of detrended time series III compared with the raw scalogram as shown in Fig. 6. The plot is counted in minutes from the beginning of the observation at 00:05 UT with its optimum cutoff frequency f o3 ¼ 0:1560 Hz.

scintillation statistical indices, which implicitly assumes strong phase scintillations accompany with strong amplitude scintillations. However, there exists an implicit limitation related to geophysical conditions of scintillation observations. For instance, amplitude scintillations isolated from phase scintillations frequently erupt in the equatorial anomaly region. The absence of phase scintillations leads to a small correlation coefficient value, which is against the expectation of the detrending evaluation. Hence, the evaluation based purely on correlation of scintillation indices may obscure the optimum cutoff frequency in some cases. A proper evaluation of detrending methods is derived from the emphasis that the aim of the post-processing or semi-real time processing is to extract fluctuations of observations due to only scintillations. It is concerned with the definition of ionospheric scintillations. Amplitude scintillations are generated by diffraction and focusing mechanism jointly affecting GNSS signals. Therefore, the frequency spectrum of amplitude scintillations is restricted to lower bands than that of phase scintillations, exactly around the Fresnel frequency. The Fresnel filtering of amplitude scintillations and the spectrum spread of phase fluctuations require detrending raw scintillation time series in terms of

amplitude scintillation observations. The traditional scintillation indices measure fluctuations of observations and provide an intuitive and convenient evaluation for detrended results, which is especially applicable to general users of low precision. In this paper, a novel evaluation approach combining the cutoff frequency optimization is presented. The detrending process with the cutoff frequency optimization implicitly depends on determining the optimum detrended results, namely the optimum evaluation. As a result, scintillation components of raw time series, even weak scintillations or those at large scales are preserved and trend components are adequately eliminated. Hence, the comparison and correlation of scintillation indices before and after detrending are not presented in this paper. Note that the optimum cutoff frequency is derived through the evolution from under-detrending to over-detrending. And when investigating the influence of varying cutoff frequencies on wavelet energy and wavelet entropy distribution of scintillation data, the degree of disorder/order pertaining to scintillation components is mainly analyzed and that of trend components at large scales is reasonably simplified. Considering that amplitude scintillation data involving in this paper was observed in high latitudes, the application of the detrending method to equatorial scintillation cases should cautiously take account of the fade stretching (Kintner et al., 2001). In terms of the fading timescale, the fade stretching phenomenon indicates significantly reduced fading rates of equatorial scintillation data, which might induce an overlap between scintillation component band and frequency content of trend components. Thus special attention should be paid to fade stretching and appropriate detrending for the equatorial scintillation study.

5. Conclusions Based on detrending raw ionospheric scintillation time series with the Fourier filter, this paper investigated the variation of wavelet energy and wavelet entropy distribution of scintillation data due to varying cutoff frequencies, presented the optimum cutoff frequency criteria, and characterized ionospheric scintillation observations for three types of scintillation time series, namely concentrated wavelet energy distribution, wide wavelet energy distribution and evident trends with negligible scintillations. With the optimum cutoff frequency derived in this paper, the detrending can extract the amplitude fluctuations purely induced by scintillations from impaired GNSS signal amplitude during propagation. Exactly, the aim of detrending ionospheric scintillation observations is to characterize the fluctuations only related with scintillations. And the optimum cutoff frequency is derived for the specific time series rather than roughly utilizing the fixed and universally applicable cutoff frequency for detrending. Hence, the optimal detrended results in spite of some

Y. Su et al. / Advances in Space Research 54 (2014) 2172–2183

inherent limitations provide reliable observations and model parameters for ionospheric scintillation investigations. Recent analyses of ionospheric scintillation time series are mainly based on wavelet. Whereas the non-stationarity of scintillation data attracts intense interest and discussions. It is of concern to apply some novel analysis tools for ionospheric scintillation observations, such as Empirical Mode Decomposition (EMD) and multifractal analysis in order to overcome some inherent imperfections of wavelet, including minimum resolution limits and edge artifacts. Furthermore, establishing a relation between the detrended observations and ionospheric scintillation models, such as the propagation model and ionospheric irregularity parameter model is another study the authors of this paper attempt to investigate. Acknowledgments This work was supported by the National High Technology Research and Development Program (863 Program) of China (No. 2011AA120503). Infrastructure funding for CHAIN is provided by the Canada Foundation for Innovation and the New Brunswick Innovation Foundation. CHAIN operation is conducted in collaboration with the Canadian Space Agency. Wavelet software was provided by C. Torrence and G. Compo, and is available at URL: http://paos.colorado.edu/research/wavelets/. The authors thank the two anonymous referees for their comments and suggestions improving this paper substantially. References Aarons, J., 1982. Global morphology of ionospheric scintillations. Proc. IEEE 70 (4), 360–378. Basu, S., Groves, K.M., Basu, Su, et al., 2002. Specification and forecasting of scintillations in communication/navigation links: current status and future plans. J. Atmos. Sol. Terr. Phys. 64, 1745–1754. Beach, T.L., 2006. Perils of the GPS phase scintillation index (rU ). Radio Sci. 41, RS5S31. http://dx.doi.org/10.1029/2005RS003356. Beach, T.L., Baragona, C.A., 2007. Quasiperiodic scintillation and data interpretation: nongeophysical GPS amplitude fluctuations due to intersatellite interference. Radio Sci. 42, RS3010. http://dx.doi.org/ 10.1029/2006RS003532. Beach, T.L., Kintner, P.M., 1999. Simultaneous global positioning system observations of equatorial scintillations and total electron content fluctuations. J. Geophys. Res. 104 (A10), 22553–22565. Beach, T.L., Kintner, P.M., 2001. Development and use of a GPS ionospheric scintillation monitor. IEEE Trans. Geosci. Remote Sens. 39 (5), 918–928. Doherty, P.H., Delay, S.H., Valladares, C.E., et al., 2004. Ionospheric scintillation effects on GPS in the equatorial and auroral regions. J. Inst. Navig. 50 (4), 235–246, Winter 2003–2004. Forte, B., 2005. Optimum detrending of raw GPS data for scintillation measurements at auroral latitudes. J. Atmos. Sol. Terr. Phys. 67, 1100– 1109.

2183

Forte, B., 2007. On the relationship between the geometrical control of scintillation indices and the data detrending problems observed at high latitudes. Ann. Geophys. 50 (6), 699–706. Forte, B., Radicella, S.M., 2002. Problems in data treatment for ionospheric scintillation measurements. Radio Sci. 37 (6), 1096. http://dx.doi.org/10.1029/2001RS002508. Forte, B., Radicella, S.M., 2004. Geometrical control of scintillation indices: what happens for GPS satellites. Radio Sci. 39, RS5014. http://dx.doi.org/10.1029/2002RS002852. Forte, B., Radicella, S.M., Ezquer, R.G., 2002. A different approach to the analysis of GPS scintillation data. Ann. Geophys. 45, 551–561. Forte, B., Materassi, M., Alfonsi, L., et al., 2011. Optimum parameter for estimating phase fluctuations on transionospheric signals at high latitudes. Adv. Space Res. 47, 2188–2193. Fremouw, E.J., Leadabrand, R.L., Livingston, R.C., et al., 1978. Early results from the DNA Wideband satellite experiment—complex-signal scintillation. Radio Sci. 13 (1), 167–187. Jayachandran, P.T., Langley, R.B., MacDougall, J.W., et al., 2009. Canadian high arctic ionospheric network (CHAIN). Radio Sci. 44, RS0A03. http://dx.doi.org/10.1029/2008RS004046. Kintner, P.M., Kil, H., Beach, T.L., et al., 2001. Fading timescales associated with GPS signals and potential consequences. Radio Sci. 36 (4), 731–743. Kintner, P.M., Ledvina, B.M., De Paula, E.R., 2007. GPS and ionospheric scintillations. Space Weather 5, S09003. http://dx.doi.org/ 10.1029/2006SW000260. Knepp, D.L., Wittwer, L.A., 1984. Simulation of wide bandwidth signals that have propagated through random media. Radio Sci. 19 (1), 303–318. Materassi, M., Mitchell, C.N., 2007. Wavelet analysis of GPS amplitude scintillation: a case study. Radio Sci. 42, RS1004. http://dx.doi.org/ 10.1029/2005RS003415. Materassi, M., Alfonsi, L., De Franceschi, G., et al., 2009. Detrend effect on the scalograms of GPS power scintillation. Adv. Space Res. 43, 1740–1748. Mushini, S.C., Jayachandran, P.T., Langley, R.B., et al., 2012. Improved amplitude- and phase-scintillation indices derived from wavelet detrended high-latitude GPS data. GPS Solution 16, 363–373. Pi, X., Mannucci, A.J., Lindqwister, U.J., Ho, C.M., 1997. Monitoring of global ionospheric irregularities using the worldwide GPS network. Geophys. Res. Lett. 24 (18), 2283–2286. Rino, C.L., 1976. Ionospheric scintillation theory-a mini-review. IEEE Trans. Antennas Propag. 24 (6), 912–915. Rino, C.L., 1982. On the application of phase screen models to the interpretation of ionospheric scintillation data. Radio Sci. 17 (4), 855–867. Rosso, O.A., Blanco, S., Yordanova, J., et al., 2001. Wavelet entropy: a new tool for analysis of short duration brain electrical signals. J. Neurosci. Methods 105, 65–75. Secan, J.A., Bussey, R.M., Fremouw, E.J., et al., 1995. An improved model of equatorial scintillation. Radio Sci. 30 (3), 607–617. Secan, J.A., Bussey, R.M., Fremouw, E.J., et al., 1997. High-latitude upgrade to the Wideband ionospheric scintillation model. Radio Sci. 32 (4), 1567–1574. Torrence, C., Compo, G.P., 1998. A practical guide to wavelet analysis. Bull. Am. Meteorol. Soc. 79 (1), 61–78 (Boston). Van Dierendonck, A.J., Arbesser-Rastburg, B., 2004. Measuring Ionospheric Scintillation in the Equatorial Region over Africa, Including Measurements from SBAS Geostationary Satellite Signals, Paper Presented at Institute of Navigation GNSS 2004, Inst. of Navig., Long Beach, California. Van Dierendonck, A.J., Klobuchar, J., Hua, Q., 1993. Ionospheric Scintillation Monitoring using Commercial Single Frequency C/A Code Receivers, Paper Presented at Institute of Navigation GPS-93, Inst. of Navig., Arlington, VA, September 1993.