Multiscale Analysis of GPS Time Series from Non-decimated Wavelet to Investigate the Effects of Ionospheric Scintillation^†

Abstracts

Due to the numerous application possibilities, the theory of wavelets has been applied in several areas of research. The Discrete Wavelet Transform is the most known version. However, the downsampling required for its calculation makes it sensitive to the origin, what is not ideal for some applications, mainly in time series. On the other hand, the Non-Decimated Discrete Wavelet Transform is shift invariant, because it considers all the elements of the sample, by eliminating the downsampling and, consequently, represents a time series with the same number of coefficients at each scale. In the present paper, the objective is to present the theorical aspects of the a multiscale/multiresolution analysis of non-stationary time series from non-decimated wavelets in terms of its implementation using the same pyramidal algorithm of the decimated wavelet transform. An application with real time series of the effect of the ionospheric scintillation on GPS satellite signals is investigated. Results obtained in this investigation of the time series of the S4 scintillation index indicate the presence of a pattern that is repeated in the series on consecutive days. Through decomposition in multiscale by wavelets, this periodic daily effect could be identified in smoother scales of the wavelet periodogram, estimated and corrected from the original time series. Therefore, this paper presents an innovative methodology for the analysis of the effects of ionospheric scintillation on GPS signals, and take a first step toward separating the effect of scintillation in relations the other effects that may influence its analysis.

Non-decimated wavelets; maximum overlap discrete wavelet transform; multiescale analysis; time series; ionospheric scintillation; GPS

Devido às inúmeras possibilidades de aplicação, a teoria de wavelets tem sido utilizada nas mais diversas áreas de pesquisa. A Transformada Wavelet Discreta é a versão mais conhecida. Porém, o processo de decimação necessário para seu cálculo, faz com que ela seja sensível à origem, o que para algumas aplicações não é o ideal, principalmente em séries temporais. Ao contrário, a Transformada Wavelet Discreta Não Decimada é invariante à translação, pois leva em consideração todos os elementos da amostra, de modo a representar uma série temporal com o mesmo número de coeficientes em cada escala, tendo o processo de decimação eliminado. Neste artigo pretende-se apresentar os aspectos teóricos da análise multiescala/multirresolução de séries temporais não estacionárias a partir de wavelets não decimadas utilizando, em termos de sua implementação, o mesmo algoritmo piramidal da transformada wavelet não decimada. Uma aplicação com a série temporal real do efeito da cintilação ionosférica sobre sinais de satélite GPS será investigada. Resultados desta investigação das séries obtidas dos índices S4 de cintilação indicam a presença de um padrão que se repete na série em dias consecutivos. Através da decomposição em multiescala por wavelets, esse efeito periódico diário pode ser identificado nas escalas mais suaves do periodograma wavelet, estimado e corrigido da série temporal original. Portanto, este artigo apresenta uma metodologia inovadora para a análise dos efeitos da cintilação ionosférica sobre sinais GPS, e dá um primeiro passo em direção a separar o efeito da cintilação de outros efeitos que podem influenciar sua análise.

wavelets não-decimadas; transformada wavelet não-decimada; análise multiescala; séries temporais; cintilação ionosférica; GPS

1 INTRODUCTION

The theory of wavelets has been quite widespread, providing major advances in several different areas of the science. There are numerous possibilities of aplication, being used for modeling,decomposition, removal of noise and/or unwanted effects, compression, detection of singularities, among other applications ⁶

In general, the decimated version (Discrete Wavelet Transform - DWT) is the most well known. It is calculated through an efficient algorithm developed by Mallat ², called Pyramidal Algorithm, which uses discrete filters and downsampling. This process eliminates the even or odd samples, so that at level j there is half of the coefficients of the previous level j ? 1. Applying an odd downsampling is similar to applying an even downsampling on data shifted by one unit. The DWT depends on the choice of downsampling, which corresponds to the choice of an origin in time or space ⁹

In this application, the multiscale analysis is applied to non-stationary time series (TS) from GPS satellite signals. Some previous applications have obtained good results applying DWT to improve GPS positioning methods ^{11 12 13}. However, to analyse TS it is preferable that the method is insensitive to the origin, which is possible with the shift invariant Non-Decimated Wavelet Transform (NDWT) ⁶.

Another advantage of using the NDWT is its easy implementation. It is possible to implement the NDWT applying the DWT pyramidal algorithm twice, first in the original vector and then in the shifted vector, eliminating the downsampling and causing that all the elements are considered at each stage of the algorithm.

The growing influence of GPS (Global Positioning System) in the navigation and remote sensing, is quite evident nowadays. However, until satellite signal gets to the GPS receiver, effects, often severe, influence the signal, causing errors or even the loss of the satellite signal by the receiver.

The ionosphere is a plasmatic region of the atmosphere that contains ions and electrons in sufficient quantity to affect signals of frequency band between 3 and 30 KHz that passes through it. The ionospheric scintillations, which are irregularities in the density of ions in the ionosphere is even worse. Its effects are more severe at high latitudes, where their occurrence is related to periods of high solar activity, magnetic storms and other extreme activities, and in equatorial and low latitudes regions, where the equatorial ionization anomaly (EIA) occurs ¹. The use of wavelets in the study of the scintillation have contributed for its better characterization in high latitude regions ^{3 5 16}. However, for equatorial and low latitude regions, in which is included Brazil, there is still no work in that direction.

For high-precision applications it is necessary the correct understanding of the effect of the ionospheric scintillation on GPS signals, enabling the characterization of this effect to allow the development of methodologies to reduce their interference.

In this paper it is intended to investigate the behavior of the effect of the ionospheric scintillation on GPS signals from a TS approach that goes beyond of the classical methods, which assume stationarity. Due to the irregularities and scattering of the ionosphere, the effects of scintillation in the satellite signals may not be stationary. The use of the NDWT for a complete analysis in multiscale will make information and hidden patterns that could not be detected in the time or frequency domain can be explained in the time/frequency domain ⁷. Even if some Fourier methods could be applied to non-stationary TS ^{8 10}, when some effects need to be investigated and separated in different scales at once, the wavelets methods became advantageous. This is because the wavelets can deal naturally with the severe non-stationarity in time caused by spikes of different magnitudes, identifying their frequency of occurrence, localization in time, and making a reliable approximation of magnitude of this effects.

2 THEORETICAL ASPECTS

Wavelet is the function (t) ( L ² (ℝ) that satisfies the following properties ⁹:

where (t) is function of time t ( ℝ.

There are two known functions that perform a primary role in the wavelet analysis, the scaling function φ called father wavelet and the wavelet φ, also called mother wavelet.

The father wavelet, discretized form, in can be expressed as

where h_k is a low-pass filter, called scaling filter, and k represents the scale covered by the filter.

On the other hand, in the relationship between the filters h_k and g_k ,

g_k = (-1)^kh_1-k

called quadrature mirror filter relation, it is considerated that the mother wavelet φ can be defined as

where g_k is a high-pass filter called wavelet filter.

Through appropriate dilations and translations for φ and φ, described by

(_j,k (t) = 2^-j/2(2^{− j} t − k) and (_j,k(t) = 2⁻ j/ ²(2⁻ j t − k), j, k ( ℤ

it is possible to construct bases for various function spaces. This is the main motivation for the development of the Wavelet Transform (WT), whose goal is to expand f ( L ² (ℝ) using a set of wavelet functions that have time/frequency location.

Thus, it is considered that any signal f(t) ( L ² (ℝ) can be represented as ⁴

where

and j ₀ represents the lowest level of resolution.

The coefficients c, d described in (2.1), called smooth (or scaling) and detail (or wavelet) coefficients, respectively, are result of the filtering of the function f(t) ( L ² (ℝ), by scaling and wavelets parent filters.

Aiming at practical application with sampled signals the following notation of DWT, in particular the NDWT, will be used to facilitate its implementation.

2.1 NDWT

The great motivation for the use of the NDWT in this work is to use a transform that operates similarly to DWT, but does not suffer sensitivity to the origin choice in time for a TS, i.e., it is shift invariant. This sensitivity of the DWT is due entirely to dowsample the output from the scaling and wavelet filters at each scale (or level) of the pyramid algorithm, where for each two outputs of the filter, one output is discarded.

The idea of the NDWT is to eliminate the downsampling, considering all the elements, and representing a TS with the same number of coefficients at each scale.

The outputs of the filters that are discarded at the first level of the DWT pyramid algorithm, can be obtained by applying the DWT pyramid algorithm in the shifted vector T X instead of X, where X is a sequence of observations (or a TS). This suggests that the NDWT can be obtained by applying the usual DWT pyramidal algorithm twice, once to X and the other to T X, and after merging the two sets of DWT coefficients together ⁹.

The NDWT scaling and wavelet coefficients, of level j, are calculed using the scaling coefficients of level j ? 1, as in DWT. However, downsampling does not occur and the filters are modified at each scale.

The NDWT scaling and wavelet coeffients can be seen as the result of the filtering of a TS X with the NDWT scaling and wavelet filters, which are presented below.

The NDWT scaling and wavelet K-length filters, must satisfy

and

respectively, for a nonzero integer n.

So, given a TS {X_t : t = 0, ... N - 1}, the result of the filtering of {X_t } with the NDWT scaling and wavelet filters is given, respectively, by

These two sequences given by equation (2.2) represent the NDWT of level j, and stipulate that the elements of , and are obtained by circularly filtering {X_t } with the filters , and , respectively. The term "mod" of this equation allows a circularly filtering making X to be represented with the same numbers of coefficients at each scale.

Moreover, it can be shown (see section 5.5 (⁹)) that it is possible to obtain and by filtering of by the following equation

The equations in (2.3) represent the NDWT pyramid algorithm and they can be written in matrix form as

where the rows of contain circularly shifted versions of after it has been upsampled to width 2j ^-1 (k - 1) + 1 (that consists of inserting 2j ^-1 zeros between each of the k values of the original filter) and then periodized to lenght N, and with a similar construction for Ãj based upon ⁹.

The NDWT also allows to reconstruct from and . The inverse DWT can be calculated via inverse pyramidal algorithm described by the following equation

or, in matrix form,

Through the wavelet periodogram, the energy of the process (or TS) {Xt} is decomposed indifferent scales and with time/frequency location. It is calculated from the NDWT wavelet coefficients given in (2.2) and can be described by (⁴)

where is calculated by the equation (2.2).

The periodogram provides a good description of where the significant changes are located in the function. Abrupt changes in the function correspond proportionately large magnitude of the coefficients.

From of the wavelet periodogram, a decomposition of variances is obtained in multiscale of the TS, i.e., in different levels of resolution j. Based on these results, it is possible to calculate the variance (energy) at each scale, with the goal to know what scale has the most energy. This task can be executed by summing all coefficients at each scale, according to the following equation

and we have the global wavelet spectrum or power spectral density (PSD).

3 APPLICATION AND RESULTS

In March 2010, the CIGALA (Concept for Ionospheric Scintillation Mitigation for Professional GNSS in Latin America) project started to investigate the causes of the ionospheric scintillation on low latitudes and in the equatorial regions, as well as its effects on technologies under development. The choice of Latin America, especially Brazil is due to two reasons: it is in one of the most critical regions under effect of ionospheric irregularities and it has academic and commercial institutions of relevance to implementation and development of the project, such as UNESP (http://is-cigala-calibra.fct.unesp.br/is/).

With this project, 12 stations with high sampling rate receivers were installed in strategic places, as illustrated in Figure 1(a), so that the collections, maintenance and re-occupations are possible for validation of possible methodologies and algorithms developed. Thus, important data for the study of the ionospheric scintillation are available since February 14, 2011.

Figure 1:
(a) Stations of CIGALA network, highlightingstation PRU1 used in this paper. Source: http://is-cigala-calibra.fct.unesp.br; (b) Cycle 24 - period of maximum solar activity in 2013. Source: http://solarscience.msfc.nasa.gov/predict.shtml.

In order to investigate the TS of the effect of the ionospheric scintillation on GPS satellite signals, and evaluate the proposed methodology, the data of the CIGALA Project available for the entire year of 2012 were used. The PRU1 station located at UNESP, Presidente Prudente, was chosen and the GPS satellite 11 was selected due to its highest scintillation index. The receiver, a Septentrio PolaRxS, represents the state of the art on receivers for tracking multi constellations and triple-frequency in the monitoring of the ionosphere. It has technology for extremely low noise and allows gathering intervals up to 100 Hz. It is possible to obtain the S4 scintillation index that can be used to characterize and help to understand the ionospheric irregularities that cause scintillation, as well as establish strategies for forecasting effects of the scintillation that can make the GPS ineffective ¹⁴.

The index S4 is defined by

where the operator represents the average (expectation E) of the intensity (I) of ionospheric scintillation. The S4 index numerator represents the variance of I (E[I ²] - E ²[I]) and the denominator, the squared average of this intensity. Thus, S4 index is a unitless measure whose square root is like the coeficient of variation in statistics. Although S4 amplitude index does not measure only scintilation effects, it has been used as an indicator of this effect in the literature. If S4 ³ 1.0 the level of scintillation is classified as strong, when 0.5 ? S4 ? 1.0 the scintillation is moderate, and for 0 ? S4 ? 0.5 the scintillation is weak ¹⁵.

The ionosphere is a very unstable layer that undergoes variations on several time scales - during the day, throughout the seasons and solar cycles that occur at each 11 years (Fig. 1(b)). In 2012, the highest index S4 occured in March, as shown in Figure 2, where the index S4 exceeds 2.1.

Figure 2:
Scintillation index S4 from 01/01/2012 to 12/31/2012.

With a "zoom" on period from 03/12/2012 to 03/14/2012, it can be seen more clearly in Figure 3 (top) that the TS of scintillation index presents gaps caused by the lack of data when the satellites are not being tracked (i.e, daily gaps that occur because the satellite is invisible, below the horizon) and periodic behavior at each day when there are data.

Figure 3:
Scintillation index S4 from 03/12/2012 to 03/14/2012 (period of strong scintillation) with elevation angle of the satellite 11 plotted in the graph below.

This periodic behavior is related to effects that repeats about daily, such as multipath effect due to reflections of satellite signals on surfaces nearby of the receptor. If the receiver environment remains unchanged, the multipath effect changes according to the movement of the satellite, which has daily repeatability. It is noted also that this periodic behavior has a shape of "U". Whereas the elevation angle between the satellite pointing direction and the local horizontal plane follows a shape of inverted "U"(bottom of the Figure 3), it is expected that the reflections (multipath) and noises are more expressive to lower angles.

It is important to eliminate or separate such daily behavior that does not characterize scintillation because it may influence the analysis of ionospheric scintillation index S4.

3.1 Results

The NDWT was applied to obtain the wavelet and scaling coefficients, given by 2.2, considering the wavelet family Symmlets with 10 vanishing moments (SYM10).

The Figure 4(a) presents the TS of the index S4 of the period from 06/12/2012 to 06/14/2012 that caracterizes weak ionospheric scintillation, as preconized by the literature. Extracting information in multiscale, as theory presented in section 2, we obtain in (b) the daily behavior of the TS in Figure 3 evidenced by repeatability at their three smoothest scales. So, we propose in this paper to estimate this behavior, which is evidenced at smoother scales of the wavelet periodogram and subtract it from the TS. Thus, in (c), the three smoothest scales are reconstructed to obtain a estimation of such behavior (or effect). For a more detailed analysis noindent multiscale, it was plotted (Fig. 4(d)) the wavelet spectrum estimated by the wavelet periodogram, described by the equation 2.4, and the global spectrum, which represents the energy or the total variance at each level of resolution. We can note in the global wavelet spectrum that the scale of level of resolution 9 is the greatest one. By interpretation of the frequency of resolution levels, the 9^thlevel is related to the effects of 512 to 1024 minutes. This demonstrates that this effect is the closest to the daily behavior, which reinforces the existence of daily variations, which include the multipath, on scintillation index S4.

Figure 4:
Estimation of the periodic behavior presented in the ST of scintillation index S4.

It is important to highlight that the wavelet multiscale decomposition is essential to identify in which scales the daily repeatability occurs. This methodology has advantages over methods in time domain or only in frequency domain, even those can deal with missing data, such as Fourier SLEX analysis or Lomb-Scargle periodograms ^{8 10}. Wavelets can deal naturally with the severe non-stationarity in time caused by sudden spikes that change part of the U-shapes: sometimes in the beginning, middle and/or end of the U-shape, with different magnitudes. Hence, the difficulty is not to identify the frequency of occurrence of the periodicity, but to estimate only this effect without an overestimation due to mainly the spikes in edges of the U-shape.

Having extracted the daily behavior of the TS of index S4 of the period of strong scintillation (03/12/2012 to 03/14/2012), a TS corrected from the cyclical effect is plotted in Figure 5.

Figure 5:
Index S4 of the period with strong scintillation, with smoother scales (estimated in the period without scintillation) removed.

Decomposing in multiscale the TS plotted in Figure 5, we can observe in Figure 6 that the scale relative to level of resolution 6 of the periodogram is most evident, i.e., it concentrates higher energy. This scale is closest to the effects of 1 hour. Having removed the daily effect of the serie, the periodic effects are not expressive in the spectrum anymore, not influencing the analysis of the scintillation index TS. It is also important to notice that the peak of the scintillation was very evident in the wavelet analysis, showing the behavior of the most abrupt change in the series of index S4 and where it occurs.

Figure 6:
Periodogram of the scintillation index with smoother scales removed.

4 CONCLUSIONS

In this article a new methodology is presented for the investigation of the ionospheric scintillation through the wavelet spectrum of the TS of the index S4 on GPS signals. It was possible to give a first step in relation to separate scintillation effect from other daily effects (mainly multipath) that may influence the S4 index.

After the estimation and removal of the daily effect, corrected series were obtained, making possible the investigation of an index that represents better the ionospheric scintillation. It is noteworthy that for effective analysis of the scintillation it is necessary to analyze a longer period of data, all available satellites and the other stations.

REFERENCES

• ¹
  M.C. Kelley. The earths ionosphere: plasma physics and electrodynamics. San Diego: Academic Press (1989).
• ²
  S.G. Mallat. A theory for multiresolution signal decomposition: The wavelet representation. IEEE Trans. on Pattern Anal. and Mach. Intell., 11 (1989), 674-693.
• ³
  M. Materassi et al. Detrend effect on the scalograms of GPS power scintillation. Advances in Space Research, 43(11) (2009), 1740-1748.
• ⁴
  P.A. Morettin. Ondas e Ondaletas: da análise de Fourier à análise de ondaletas. EDUSP, São Paulo (1999).
• ⁵
  S.C. Mushini et al. Improved amplitude- and phase-scintillation indices derived from wavelet detrended high-latitude GPS data. GPS Solution, 16 (2012), 363-373.
• ⁶
  G.P. Nason. Wavelet Methods in Statistics with R. Springer, New York (2008).
• ⁷
  G.P. Nason & B.W. Silverman. The stationary wavelet transform and some statistical applications. Lecture Notes in Statistics, 103 Springer-Verlag (1995), 281-299.
• ⁸
  H. Ombao, R. von Sachs & W. Guo. SLEX analysis of multivariate non-stationary time series. J. of the Am. Stat. Assoc., 100(470) (Jun., 2005), 519-531.
• ⁹
  D.B. Percival & A.T. Walden. Wavelets Methods for Time Series Analysis. Cambridge University Press, Cambridge, England (2000).
• ¹⁰
  J.D. Scargle. Statistical aspects of spectral analysis of unevenly spaced data. Astrophys. J., 263 (1992), 835-853.
• ¹¹
  E.M. de Souza et al. Processamento de Sinais com Wavelets: Redução do Efeito do Multicaminho do Sinal GPS. TEMA - Tend. Mat. Apl. Comput., 6(1) (2005), 177-186.
• ¹²
  E.M. de Souza et al. Comparação das Bases de Wavelets Ortonormais e Biortogonais: Implementação, Vantagens e Desvantagens no Posicionamento com GPS. TEMA - Tend. Mat. Apl. Comput., 8 (2007), 149-158.
• ¹³
  E.M. de Souza et al. Método dos Mínimos Quadrados com Regressão Wavelets para Redução de Erros no Processamento de Sinais GPS. TEMA - Tend. Mat. Apl. Comput., 8 (2007), 139-147.
• ¹⁴
  J.H. Strangeways. Determining scintillation effects on GPS receivers. Radio Science, 44 (2009), doi:10.1029/2008RS004076.
  » https://doi.org/10.1029/2008RS004076
• ¹⁵
  R. Tiwari et al. WBMod Assisted PLL GPS Software Receiverfor Mitigating Scintillation Affect in High Latitude Region. URSI - Int. Union of Radio Sci., IEEE (2011).
• ¹⁶
  R. Tiwari et al. Investigation of ionospheric irregularities and scintillation using TEC at high latitude. Advances in Space Research, 52(6) (2013), 1111-1124.

• †
  Thanks to Foundation for Research Support of the State of São Paulo - FAPESP, for the financial support; Paper presented at CMAC - Southeast in 2013.

Publication Dates

• Publication in this collection
  Aug 2015

History

• Received
  30 Nov 2013
• Accepted
  15 Apr 2015