VISUAL ANALYSIS OF RECURRENCE OF TIME SERIES OF THE COORDINATES ENU IN THE GPS STATIONS

Análisis Visual de Recurrencia de las Series de Tiempo de las Coordenadas ENU en las Estaciones GPS

Abstract

The time series content information about the dynamic behavior of the system under study. This behavior could be complex, irregular and no lineal. For this reason, it is necessary to study new models that can solve this dynamic more satisfactorily. In this work a visual analysis of recurrence from time series of the coordinate’s variation ENU (East, North, Up) will be made. This analysis was obtained from nine continuous monitoring stations GPS (Global Navigation Satellite System); the intention is to study their behavior, they belong to the Equatorian GPS Network that materializes the reference system SIRGAS - ECUADOR. The presence of noise in the observations was reduced using digital low pass filters with Finite Impulse Response (FIR). For these series, the time delay was determined using the average mutual information, and for the minimum embedding dimension the False Nearest Neighbours (FNN) method was used; the purpose is to obtain the recurrent maps of each coordinates. The results of visual analysis show a strong tendency, especially in the East and North coordinates, while the Up coordinates indicate discontinued, symmetric and periodic behavior.

Keywords
Coordinates ENU; Time Series; Time Delay; Embedding Dimension; Recurrence

Resumen

Las series temporales contienen información acerca del comportamiento dinámico de un determinado sistema. Este comportamiento puede presentarse como complejo, irregular y no lineal. Por esta razón es necesario que se realicen estudios con nuevos modelos que satisfagan de mejor manera sus resultados. En este trabajo se hace un análisis en forma visual sobre la recurrencia de la variación de las coordenadas ENU (Este, Norte, Up-arriba) a partir de las series temporales obtenidas de las estaciones de monitoreo continuo GPS con la intención de ver que estructura dinámica presentan esta series. Para lo cual, fue reducido el ruido existente en las series usando filtros pasa-bajo con FIR (Finite Impulse Response). Con la finalidad de construir los mapas recurrentes de las coordenadas ENU fueron calculados el tiempo de retardo mediante el uso de la técnica ami (average mutual information), y la dimensión embebida usando el método de los falsos vecinos más cercanos (False Nearest Neighbours-FNN). Los resultados del análisis visual muestran que las coordenadas Este y Norte tienen una fuerte tendencia, mientras que la coordenada Up (arriba) indica la presencia de discontinuidades, además de un comportamiento simétrico y periódico.

Palabras-Clave
Coordenadas ENU; Series de tiempo; Tiempo de retardo; Dimensión Embebida; Recurrencia

1. Introduction

Generally, the analysis of time series of some phenomena that occurred in nature, has been made through linear models. Nevertheless, there are phenomena that have complex, irregular and non-linear characteristics. A non-linear deterministic system could generate series with aleatory appearance which make impossible to detect their deterministic nature with techniques of linear identification (Liu et al. 1993).

That is why they have to be analyzed with other techniques or theories, in order to know their dynamic behavior with more precision. One of these techniques is the dynamic system theory, which has been used in different fields of science. This theory allows the analysis of the behavior of these complex time series in order to determine if it is chaotic, deterministic or random (Gao et al. 2007).

The chaotic mathematics is an alternative in the statistic or probabilistic analysis of the phenomena that show a temporal evolution; apparently aleatory, irregular, aperiodic, and scarcely predictable. There are especially chaotic phenomena that are governed by deterministic and structurally simple models. The evidence of chaos, however weak it may be, is important and of interest to improve predictions about its future evolution (Scaffeta 2010).

There are several non-linear techniques or measures to determine the complexity of a time series, for example the correlation dimension, Lyapunov exponents, Poincare maps, Hurst coefficient, recurrence visual analysis, among others (Wolf et al. 1985, Nychka et al. 1992, Tierra et al. 2018). These technics allow the identification of non-linear presence and possibly chaotic behavior in these series. These techniques have been used in different fields as geodynamic (Tiwari and Krishnaveni 2015), economy (Carrasco et al.2015), vehicular traffic (Mancero et al.2017), and weather (Zhang et al. 2016), among others.

In the case of time series of the East, North, Up (ENU) coordinates obtained from the GPS continuous monitoring, it is hard to justify that they have a non-linear behavior. It is necessary to introduce new concepts and techniques related with nonlinearity and to determinate their dynamic behavior. This will allow the series to be modeled with more reliability and precision, and to have more accurate predictions.

In this analysis, a first study of the time series of the coordinates ENU is made; obtained from nine continuous monitoring stations which belong to the Ecuadorian GPS Network. The recurrence visual analysis technic will allow the construction of the respective recurrent two-dimension maps of the ENU coordinates. With the visual analysis of these maps, it will be possible to identify the presence of patterns or structures that indicate their dynamic behavior.

2. Visual Analysis of Recurrence

The VAR (Visual Analysis of Recurrence) technique uses recurrence graphics to detect patterns and structural changes hidden in the data, in order to find hints of the chaotic behavior in the time series. The recurrence maps are a set of points, equally spaced in a square matrix of NxN (N is the total data of time series), where the axes represent chronologic sequence of the vectors in the reconstructed space. When the recurrence map shows an identifiable structure or pattern, the series indicate a sign of determinism; otherwise, the series indicate a sign of randomness (Ekmann et al. 1987). This technic allows investigating the trajectory (system dynamics) in an m-dimensional phase space through a two-dimension representation of its recurrences. This recurrence at state of time i at different time j is visualized within a two-dimension square matrix (recurrence matrix) with black and white dots, where black dots mark a recurrence and the two axes represent the time (Iwanski and Bradley 1998, Marwan et al. 2007). The recurrence graph is a qualitative process that allows observing hints of system behavior in order to localize hidden recurrent patterns.

According to Takens’s theorem an image can be topologically reconstructed equivalent to an original multidimensional system from a scalar time series x_i (Takens 1981). This could be realized with the appropriate knowledge of the time delay (τ) and the embedded dimension (m) getting:

$$x_i^m=\left(x_i,x_{i+\tau },x_{i+2\tau },\dots ,x_{i+\left(m-1\right)\tau }\right)$$ (1)

The time delay can be computed by the average mutual information function (AMI) and the first minimum is chosen. With this time delay, the embedding dimension can be calculated by false nearest neighbors method (FNN), and the value is chosen when the minimum rate is achieved (Cao 1992, Kennel et al. 1992, Frazer and Swinney 1986). Both AMI and FNN were used in this study, because they are standard approaches for finding optimal parameters for time delay and embedded dimension respectively. The recurrence image can be expressed by matrix according to eq(2) (Marwan et al. 2007):

$$R_{ij}\left(\epsilon \right)=H\left(\epsilon -‖x_i-x_j‖\right),i,j=1,\dots N$$ (2)

Where N is the total number points of the time series, ε is a threshold distance, and H is the Heaviside function.

Attention should be paid on how to choose the threshold ε., it is recommended to choose an optimal threshold. Suggestions from literature show that this threshold should not exceed 10% the average of the maximum diameter of the phase space (Zbilut and Weber 1992; Koebbe and Myare-Kress 1992; Mindlin and Gilmor 1992).

The time series, depending of their evolutionary process, may present different structures, which indicate the existence of periodicity, homogeneity, tendency, disorganized, stationary, signal with white noise, non-stationary persistent, chaotic, black bands, among others (Zbilut and Weber 1992).

3. Times series of the ENU coordinates in GPS station

The determination of time series can be developed with coordinates and their uncertainties obtained from the processing of GPS data. This data belongs to 45 sites of the Equatorian GPS Network (green triangles). The Equatorian GPS Network materializes the reference system SIRGAS - ECUADOR in the country, and its sites are distributed all around the country. In this study, time series of ENU (East, North, Up) coordinates of nine stations were used, between latitude 1° N - 3.25°S and longitude 77°W-80°W (Figure 1). The sites are continuously monitoring and capturing data for 24 hours a day and 365 days per year. In this way, these sites provide necessary information for GPS processing considering the variation of coordinates within the time. The data processing was processed everyday using GAMIT-GLOBK scientific software (Herring et al. 2015) following the standards of SIRGAS. The ENU coordinates in the ITRF2008 from year 2012.925 to 2016.999 were calculated. The displacement in millimeters with respect to an a priori position was obtained, and differences de, dn, and du in each coordinate were generated.

Figure 1:
Map of the Ecuadorian GPS Network. The nine GPS stations used to generate the time series are indicated with the red triangles.

a. Preprocessing Temporal Sequences

Temporal sequences analyses, using digital signal processing techniques in the time and frequency domains, have evolved significantly in the last 40 years (León 2014, Box and Jenkins 1976, Oppenheim and Schafer 2010). The data of the displacements experienced by well-known monitoring stations distributed across the country, with respect to a fixed point of the last four years can be preprocessed with simple techniques of low computational cost and easy to interpret. For this purpose, the database was built with a ratio of one sample per day.

b. Preprocessing in the Frequency Domain

Before preprocessing the temporal sequences with filtering techniques in the time domain, we investigated the spectral content of the sequences (de, dn and du) of each station, obtaining its respective SPECTROGRAM (Rabin and Schafer 2007, MatLab^TM ), with a duration of the processing window of 512 days and an overlap of 1 day. The results illustrated in Figure 2 show that the highest energetic content of the temporal sequences (de, dn and du) of each station is in very low frequency (below of ((/54.5) rad/sample). Similarly, Figure 3 illustrates the Fast Fourier Transform (FFT) of the temporal sequences du (León 2014, MatLab^TM ). In adition, at the highlighted window, in the center of Figure 3, it is showing that the sequence du have, fundamentally, two very low frequencies (₁ = 4(/1376 rad, and ((₂ = 12(/1376 rad) with high energy. The behavior of the other two components shown similar behavior, i.e. the highest energetic content of the temporal sequence is concentrated in very low frequencies.

Figure 2:
An example spectrogram of the temporal sequences (data = du).

Figure 3:
Fast Fourier Transform of the temporal sequences (data = du).

c. Preprocessing in the Time Domain

Although the results shown in item b, previous section, it was obtained by processing the du time sequence, our tests were applied in all three sequences. Therefore, according to our analysis, the time sequence can be modeled as the sum of very low frequency signals with Additive White Gaussian Noise (AWGN) (Van Trees et al. 2013). We projected five temporary filters: first with four known filters as the Moving Average with windows of duration of 25, 50, 75 and 100 days (Oppenheim and Schafer 2010, MatLab^TM), whose frequency response is equivalent to a low pass filter. Nevertheless, low computational cost and an excellent monitor of the instantaneous average value of the process. Second, a low pass filter, Finite Impulse Response (FIR), with cut-off frequency equal to ((/54.5 rad/second), product of analysis of results of b) an order 300 (MatLab^TM). Note that the computational cost of low pass FIR filter is higher than the Moving Average one but the first is more selective. The results of the application of the five filters to the temporal sequences (de, dn and du) of each station are illustrated in Figure 4. We can observe that all the filters smoothen the temporal sequences, mitigating what we have so far called AWGN, whose analysis will be the object of future work.

Figure 4:
Results of filtering the temporal sequences (data = du) with five filters.

From the results obtained, we will proceed with the analysis of the temporal sequences preprocessed only by the Moving Average filter of window duration equal to 50, considering that it is the one that best trade off is obtained between computational cost, delay produced and noise mitigation.

4. Results and Discussion

In order to reduce noise in the observations, the time series were pre-processed according to the sections 3.a and 3.b. As an example, Figure5a shows series in ALEC station, where EAST coordinate is illustrated in blue, and filtered series in red. Figure 5b shows the recurrence map obtained from the original series, and the presence of noise can be observed, making difficult to identify the patterns or structures in series. Figure5c shows an image of the recurrence map (without threshold) constructed with the noise filtered series, being clearer, and easier to distinguish patterns in the series.

Figure 5:
Coordinate East of the station ALEC a) Original time series in blue and filtered series in red, b) Recurrence map with noise, c) Recurrence map without noise

The remaining recurrence maps (without threshold) of ENU coordinates of the GPS stations can be observed in Table2. For this filtered series the time delay and the embedding dimension were calculated and the calculated values are given in the Table 1.

The time when the first minimum is reached was chosen to determine the optimum time delay. As example the Figure 6 shows that the value for the North coordinate of the FOEC station is 3.

Figure 6:
Time Delay of FOEC (North) station

Later, the embedding dimension was determined from previously calculated time delay, and its dimension was obtained when the minimum was reached. Figure 7 indicates that the minimum embedding dimension of North coordinate from FOEC station is 3

Figure 7:
Embedding dimension of FOEC (North) station with a time delay of 3.

Table 1:
Time delay (τ) and embedding dimension(m)

Station	East		North		Up
	τ	m	τ	m	τ	m
ALEC	5	4	6	5	4	4
COEC	6	4	10	4	5	3
EPEC	7	5	6	3	3	3
ESMR	4	3	2	3	4	3
FOEC	6	3	3	3	3	3
GUEC	6	4	4	3	3	4
LPEC	4	4	3	3	3	3
QUEM	6	4	3	4	3	4
SIEC	5	3	5	4	6	5

In Table 2, column 1, the respective recurrence maps of the time series obtained with values of the time delay and embedding dimension (calculated in the Table 1) are shown. To obtain a better view of the patterns in the series, the recurrence maps were graphed in black and white, and the 10% of the maximus radius was used. This radius was determinate from the maximum Euclidean distance between the vectors of the recurrence matrix (eq.2). The column 2 of the Table 2, shows the respective recurrence maps.

Table 2:
Recurrence maps without threshold (color) and with threshold (black and white)

Making a visual analysis of the black and white images it is possible to observe that around the diagonal E, N, coordinates, black spots are concentrated; likewise, it can be observed the presence of arrows. This behavior indicates the presence of a strong tendency in those time series. Also, it is easy to identify how the areas change (with black spots) in most of the stations, indicating the presence of discontinuities. Furthermore, the stations have different structures grouped in grid form. This indicates that they have different behaviors and it would help in the detection of structural changes.

Instead, in the Up coordinate recurrence maps present various structures. For the QUEM station, the recurrence spots form parallel lines to the main diagonal until approximately 1200 time, given hints of a possible deterministic structure, and with a period of approximately 300 days. For ALEC it can be perceived that spots form parallel and perpendicular lines to the main diagonal, indicating the presence of a symmetric and periodic structure trough the time, with a period of 300 days approximately, the same happens for COEC, EPEC, ESMR stations but at different time intervals. FOEC, GUEC, and SIEC stations present the same structure in a disorganized way, which indicate that the series presents abrupt changes in its dynamic.

From the visual analysis, the different structures presented by the recurrence maps of the time series, give hints of their dynamic behavior. Nevertheless, in order to determine their behavior in a more objective and accurate way, it is necessary to use new techniques. For this reason future works will be realized with the quantitative analysis of this recurrence maps using different metrics of the Recurrence Quantification Analysis- RQA.

5. Conclusion

According to the results obtained from the visual analysis of recurrence of the time series from the ENU coordinates, the following can be indicated

• The recurrence maps from the E,N coordinates, present structures that indicate that the process contains a tendency

• In the U coordinate, it was possible to determine that the structures present different behaviors such as periodic process, symmetrical and disorganized.

Acknowledgment

To Massachusetts Institute of Technology - MIT trough Prof. King for providing the GAMIT/GLOBK scientific software for the processing of GPS data

References

• Box, G. and Jenkins, G. 1976. Time Series Analysis: Forecasting and Control Holden-Day, San Francisco.
• Cao, L. 1992. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D 58, pp. 43-50.
• Carrasco, R., Vargas, M. , Alfaro, M., Soto, I. and Fuertes, G. 2015. Cooper Metal Price Using Chaotic Time Series Forecasting. IEEE Latin America Transactions, 13(6), pp. 1961-1965.
• Eckmann, J., Kamphorst, S., and Ruelle, D. 1987. Recurrence Plots of Dynamical Systems. Europhysics Letters, 4 (9), pp. 973-977.
• Frazer, A. and Swinney, H. 1986. Independent coordinates for strange attractors from mutual information. Physical Review A 33, pp. 1134-1140.
• Gao, J., Cao, Y., Tung,W. and Hu, J. 2007. Multiscale analysis of complex time series: Integration of chaos and random fractal theory, and Beyond Wiley & Sons, New Jersey, USA.
• Kennel, MB., Brown, R. and Abarbanel, HD. 1992. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys Rev A, 45(6), pp.3403-3411.
• Koebbe, M. and Mayre-Kress, G. 1992. Use of Recurrence Plots in the Analysis of Time-Series Data. In: Casdagli, M. Eubank, S. (Eds.), Proceedings of SFI Studies in Science of Complexity Vol. XXI. Addison-Wesley, Redwood City, pp. 361-378.
• Herring, T., King, R., Floyd, M. and McClusky, -. 2015. Introduction to GAMIT/GLOBK: Release 10.6 Massachusetts Institute of Technology, USA.
• Iwanski, J. and Bradley, E. 1998. Recurrence plot analysis: To embed or not to embed?. Chaos, 8, pp. 861-871.
• León, R. 2014. Introducción al Procesamiento Digital de Señales Estocásticas ISBN: 978-9978-301-28-9, Universidad de las Fuerzas Armadas, Ecuador.
• Liu, C., Granger, J. and Heller, W. 1993. Using de correlation exponent to decide whether an economic series is chaotic. Non-linear Dynamics, Chaos and Econometrics J. Wiley, New York.
• Mancero, H., Morales, B. and Tierra, A. 2017. Hurst coeficient estimation of traffic time series in urban areas by Range_Rescaled. Revista Geoespacial, 14(1), pp. 103-120.
• Marwan, N., Romano, M., Thiel, M. and Kurths, J. 2007. Recurrence plots for the analysis of complex systems. Physics Reports, 438, pp. 237-329.
• MatLab, MathWorks, Apple Hill Drive, Natick, MA 01760-2098, USA, Available at: <Available at: http://www.mathworks.com > [Accessed 15 September 2017].
  » http://www.mathworks.com
• Mindlin, G. and Gilmor, R. 1992. Topological Analysis and Synthesis of Chaotic Time Series. Physica D 58 (1-4), pp. 229-242.
• Nychka, D., Ellner, S., Gallant, A. and McCafrey, D. 1992. Finding Chaos in Noisy Systems. J.R. Statist. Soc. B, 54(2), pp. 399-426.
• Oppenheim, A. and Schafer, R. 2010. Discrete - Time Signal Processing, 3rd ed. ISBN: 978-0-13-198842-2, Pearson.
• Rabiner, L. and Schafer, R. 2007. Introduction to Digital Speech Processing, ISBN: 978-1-60198-070-0, NOW the essence of knowledge.
• Scaffeta, N. 2010. Fractal and Diffusion Entropy Analysis of Time Series: Theory, concepts, applications and computer codes for studying fractal noises and Lévy walk signals VDM Verlag. USA.
• Takens, F. 1981. Detecting Strange Attractors in Turbulence. Vol. 898 of Lecture Notes in Mathematics Springer, Berlin, 1981. pp. 366-381.
• Tierra, A., Luna, M., Staller, A., Pilapanta, C., Romero, R. and Porras, L. 2018. Hurst Coefficient Estimation by Rescaled Range and Wavelet of the ENU Coordinates Time Series in GPS Network. IEEE Latin America Transaction, Vol. 16, N°. 4, pp. 1064-1069.
• Tiwari, R. and Krishnaveni, P. 2015. Evidence of self-organization in Sumatra earthquakes: recurrence time and its geodynamical implications. Nat. Hazards, 77:S51-S63.
• Van Trees, H., Bell, H. and Tia, Z. 2013. Detection Estimation and Modulation Theory, Part I: Detection, Estimation and Filtering Theory 2nd ed. ISBN: 978-0-470-54296-5, John Wiley and Sons.
• Wolf, A., Swuift, J., Awinney , H. and Vastano, J. 1985. Determining Lyapunov Exponents from a Time Series. Physica D, 16, pp. 285-317.
• Zhang, W., Feng, G., and Liu, Q.2016. The application of recurrence quantification analysis in detection of abrupt climate change. Discrete Dynamics in Nature and Society, 2016, 7 pag, article ID 2689429.
• Zbilut, J., and Webber, -. 1992. Embeddings and delays as derived from quantification of recurrence plots. Physics Letters A, 171, pp. 199-203.

Publication Dates

• Publication in this collection
  Oct-Dec 2018
• Date of issue
  Dec 2018

History

• Received
  06 Oct 2017
• Accepted
  27 May 2018