Comparative Analysis of Non-Linear GNSS Geodetic Time Series Filtering Techniques: El Hierro Volcanic Process (2010–2014)

: GNSS geodetic time series analysis allows the study of the geodynamic behavior of a speciﬁc terrestrial area. These time series deﬁne the temporal evolution of the geocentric or topocentric coordinates obtained from geodetic stations, which are linear or non-linear depending, respectively, on the tectonic or volcanic–tectonic character of a region. Linear series are easily modeled but, for the study of nonlinear series, it is necessary to apply ﬁltering techniques that provide a more detailed analysis of their behavior. In this work, a comparative analysis is carried out between different ﬁltering techniques and non–linear GNSS time series analysis: 1sigma–2sigma ﬁlter, outlier ﬁlter, wavelet analysis, Kalman ﬁlter and CATS analysis (Create and Analyze Time Series). This comparative methodology is applied to the time series that describe the volcanic process of El Hierro island (2010–2014). Among them, the time series of the slope distance variation between FRON (El Hierro island) and LPAL (La Palma island) stations is studied, detecting and analyzing the different phases involved in the process.


Introduction
GNSS-GPS systems are very effective tools in the study of the geodynamic behavior of a region. From the processing of the GPS observations, geodetic time series of sub-centimetric accuracy are obtained. The analysis of these time series is essential to understand the geodynamic behavior, even distinguishing between tectonic and volcanic activities in those places where the geodynamics presents these or other complex situations.
In this work, a mathematical procedure is established for the study of nonlinear geodetic time series. The methodology consists of a pre-treatment of the time series to eliminate anomalous values that disturb the subsequent adjustments. This is performed using the Outlier R filter. To these filtered series, the analytical Kalman and wavelet filters, the statistical ARMA and ARIMA filters, and the CATS and STL techniques of linear fitting are applied.
This methodology is applied to the time series obtained from the FRON station, located on El Hierro island, in the Canary Islands. These series span from 2010 to 2014; therefore, they reflect the volcanic process that began on the island in July 2011. Due to the volcanic-tectonic geodynamic behavior of the region, the time series present non-linear characteristics. As a comparison, the time series of the IZAN station on Tenerife island are shown, which are not affected by volcanic activity. Therefore, these time series are linear.

Site Description
El Hierro island is located southwest of the Canary archipelago. The island's morphology has been interpreted as a triple volcanic rift: NE, NW and S rifts, with axes diverging about 120 • , Figure 1 [1]. In July 2011, an increase in surface deformation and seismicity on the island was detected. The climax of this unrest was a submarine eruption first detected on 10 October 2011 [2], and located at about 2 km SW of La Restinga, the southernmost village of the El Hierro island. The eruption ceased on 5 March 2012, although deformation and seismic activity did not cease after the eruption [1]. On the island, there was only one geodetic benchmark from which global navigation satellite systems (GNSS) provided continuous and publicly accessible data from the beginning of the volcanic unrest. This is the FRON station, located at the Frontera municipality in the El Golfo valley, and maintained by the Canarian Regional Government. Because of the ground deformation detected there through geodetic processing of global positioning system (GPS) data, other GNSS-GPS receivers were deployed by the Spanish National Geographic Institute (IGN) throughout a geodetic benchmarks' network designed by the Laboratory of Astronomy, Geodesy and Cartography of Cadiz University ( Figure 1). A first set of four benchmarks in the El Golfo valley was continuously observed near the end of July, forming an almost straight line (HI01, HI02, HI03, and HI04), but it was only later near the eruption's start that four other benchmarks were continuously observed, forming a three-tipped star covering all of the island spatially (HI00, HI05, HI08, and HI09). In addition, HI10 was continuously observed, while three others were only periodically observed (HI06, HI07 and HI11) [2,3].
In the rest of the islands of the Canary archipelago, there are other permanent GNSS-GPS stations managed by the IGN, among them IZAN on Tenerife island and LPAL on La Palma island. The LPAL station also belongs to the international IGS network (Figure 1).

Methodology
The evolution of the eruptive process has been studied from the analysis of the geodetic time series of the GNSS-GPS stations located on the islands. The GPS observations have been processed using the Bernese scientific software v5.0 [4]. The IGS LPAL station (La Palma island) has been used as a reference station, in sessions of 24 h and 30 s of sampling frequency, using the ITRF2008 reference frame [5]. For each observation session, geocentric coordinates (X, Y, Z) have been obtained with sub-centimeter accuracy. From a given initial time, the time series of the topocentric coordinates (east, north, up) and the time series of the distance variation between the reference station, LPAL, and the corresponding station have been built.
These time series can be linear or non-linear depending on the tectonic or volcanictectonic character of a region. Figure 2 shows the topocentric time series and distance variation time series between the LPAL-IZAN and LPAL-FRON stations from 2010 to 2014. The IZAN station, located on Tenerife island, is not affected by the volcanic process of El Hierro; therefore, it presents a linear behavior in all its components (Figure 2a). On the contrary, the FRON station is located on El Hierro island, so its time series are non-linear (Figure 2b). Linear series are easily modeled but, for the study of nonlinear series, it is necessary to apply filtering techniques that provide a more detailed analysis of their behavior. The methodology is summarized in Figure 3. For a better understanding, the distance variation time series between FRON and LPAL is shown with the result of each technique, but the next section shows the result for all the time series of FRON station.
First, a descriptive analysis of the raw data is carried out (Figure 4a). Thus, an initial visualization of the data is carried out, detecting errors due to different causes both physical and instrumental. Due to these anomalous data, it is necessary to treat the time series to eliminate outliers and noise that distort the subsequent analysis. For this reason, an initial filtering of the data is carried out using the Outlier R filter (Figure 4b). The objective of this filter is to ensure that the filtered series has linearity, homoscedasticity and follows a normal distribution. To achieve this objective, this filter uses the Box-Cox transform [6]. Subsequently, different analytical filtering techniques are applied to the filtered time series: Kalman and wavelet. The Kalman filter is an algorithm for updating, observation by observation, the linear projection of a system of variables on the set of available information, as new information becomes available. The Kalman filter makes it possible to easily calculate the likelihood of a linear, uni-equation or multi-equation dynamic model, estimating the parameters of the model, as well as obtaining predictions from these types of models [2,7]. The application of this filter to the time series is shown in Figure 4c.
Wavelet techniques allow to divide a complex function into simpler ones and study them separately. To apply the wavelet transform to a series of numerical data, it is necessary to implement the discrete wavelet transform (DWT) [8]. The objective of applying the DWT to a vector is to obtain a transformed vector that has in the middle, known as the high part (details), the same high-frequency information as the original vector and, in another half, known as the low part (approximations), the low-frequency information. Wavelet transforms comprise a large set of shapes. Over time, different versions of wavelets have been developed, which have given rise to families of wavelets [9]. In this work, the Coiflets family has been used, and specifically, the Coiflets of order 5. The result of this technique is shown in Figure 4d.
On the other hand, statistical filters are applied: ARMA and ARIMA. The ARMA model is given by the composition of autoregressive models (AR) and moving average models (MA). On the other hand, the ARIMA model results from the union of the autoregressive (AR), integrated and moving average (MA) models [10]. The results of both filters are shown in Figure 4e,f, respectively.
Finally, as adjustment and forecasting techniques, a linear adjustment, the CATS adjustment and the STL decomposition are performed. In order to carry out a linear fit and due to the non-linear characteristics of the time series, it is necessary to carry out this fit in parts (Figure 4g). The CATS adjustment (Create and Analyze Time Series) [11] consists of decomposing the time series in order to calculate the trend, the amplitudes of the sinusoidal terms and the magnitudes of the discontinuities that the series presents [8] ( Figure 4f). On the other hand, the STL decomposition (Seasonal and Trend decomposition procedure based on Loess) decomposes a time series into its three components: trend, seasonality and irregularities using local regression (loess) [12].  Figure 5 shows the results of applying the filters exposed in the methodology to the topocentric time series east (a), north (b) and height (c) of the FRON station, and to the distance variation between FRON-LPAL (d). The figures show: outlier R series (blue), wavelet series (red), Kalman series (pink), ARMA series (green), ARIMA series (orange), CATS series (light blue). The earthquakes of magnitude greater than 4 (light green color) that occurred in the region between 2011 and 2014 are also represented, obtained from the seismic catalog provided by the IGN (www.ign.es, accessed on 1 June 2021). The black line represents the earthquake of magnitude 5.1 that occurred on 27 December 2013. Data Availability Statement: Data supporting reported results can be found at http://www.ign.es, https://www.epncb.oma.be and http://www.grafcan.es.