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 xi (Takens 1981). This could be realized with the appropriate knowledge of the time delay (τ) and the embedded dimension (m) getting:
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):
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, MatLabTM ), 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, MatLabTM ). 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 (1 = 4(/1376 rad, and ((2 = 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.
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, MatLabTM), 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 (MatLabTM). 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.
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.
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
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.
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.