A Multiple Algorithm Approach to the Analysis of GNSS Coordinate Time Series for Detecting Geohazards and Anomalies

In this study, a multiple algorithm approach to the analysis of GNSS coordinate time series for detecting geohazards and anomalies is proposed. This multiple algorithm approach includes the novel use of spatial and temporal analyses. In the spatial analysis algorithm, the spatial autoregressive model was used, assuming that the GNSS coordinate time series from a network of stations are spatially dependent. Whereas in the temporal analysis algorithm, it is assumed that the GNSS coordinate time series of a single station is temporally dependent and an artificial neural network is used to extract this dependency as a nonparametric model. This multiple algorithm approach was examined using (i) the BIGF network of GNSS stations in the British Isles and (ii) the GNSS stations of the GEONET network in Japan for the Tohoku‐Oki 2011 Mw9.0 earthquake. It was demonstrated in these case studies that this multiple algorithm approach can be used to detect the effect of a geohazard such as an earthquake on the GNSS network coordinate time series and to detect regional anomalies in the GNSS coordinate time series of a network. The spatial analysis algorithm seemed to be more suitable to detect coordinate offsets in the low‐frequency component and/or variations in the long‐term trends of the GNSS coordinate time series, while it is less reliable in detecting sudden large magnitude coordinate offsets due to earthquakes, as the effects at one station propagate to nearby stations. In contrast, the temporal analysis algorithm detects coordinate offsets in the high‐frequency component which makes it effective in detecting sudden large coordinate offsets in the GNSS coordinate time series such as those due to earthquakes. Thus, it was shown the complementary of the temporal and spatial analysis algorithms and their successful application for the magnitude and frequency content of the anomalies in the two case studies.


Introduction
Monitoring and early warning (EW) systems for the detection of natural hazards are usually based on a network of sensors and related to developments in geospatial engineering (Bhattacharya et al., 2012), which enable the protection of infrastructure and populations, and the mitigation of long-term consequences (Kubo et al., 2011). The developments of GPS technology, that is, sampling rate up to 100 Hz (Häberling et al., 2015;Zhou et al., 2018), the introduction of additional satellite systems as GLONASS, BeiDou, Galileo, and so forth (Msaewe et al., 2017;Teunissen et al., 2014), and the broad operation of permanent GNSS networks (Bock & Melgar, 2016) provide continuous time series of GNSS products, which can reflect potential ground deformation (Liu et al., 2017;Reilinger et al., 2006) and troposphere/ionosphere abnormalities (Wielgosz et al., 2005), all of which are related to geohazards. More specifically, GNSS coordinate time series are used in applications for the estimation of earthquake magnitude (Blewitt et al., 2006;Wright et al., 2012) and earthquake characteristics (Geng et al., 2013;Melgar et al., 2013;Psimoulis et al., 2014), tsunamis (Ohta et al., 2012), hydrological loadings (van Dam et al., 2001;Tregoning et al., 2009), vertical land movements and sea-level rise (Teferle et al., 2006), and ionospheric storms (Wielgosz et al., 2005).
The processing mode of the GNSS network data and the analysis of the resulting GNSS coordinate time series vary according to the network characteristics (e.g., size and density) and the characteristics of the monitored phenomenon (deformation magnitude, time scale, etc.). For instance, in the application of earthquake early warning systems, high-rate (i.e., 1 Hz or higher) GNSS coordinate time series usually derived from precise point positioning (PPP) mode are used in real time or near real time, to complement existing early warning systems in order to provide the amplitude of the coseismic displacement, ground motion characteristics (PGD, PGV, etc.), and the seismic wave detection (Pérez-Campos et al., 2013) or to predict the generation of tsunamis (Blewitt et al., 2009). A similar approach is also followed for the application of GNSS networks in the monitoring of landslides and the modeling of ground motion (Benoit et al., 2015;Zhou et al., 2018). Furthermore, recent studies have revealed the potential contribution of high-rate GNSS data in earthquake early warning systems by estimating the predominant period (Psimoulis, Houlié, & Behr, 2018) and the peak displacement (Crowell et al., 2016) of the P-waves, supplementing the current seismic data-based early warning techniques.
However, most of the abovementioned studies analyze the GNSS coordinate time series individually for each station, without considering the whole network in the analysis process. Only a few studies used the network mainly for studying wave patterns by visually showing the dynamics of ground motion (Grapenthin & Freymueller, 2011) or as a spatial check to filter out false alarms and outliers in early warning methods (Psimoulis, Houlié, Habboub, et al., 2018). On the other hand, some studies used the concept of spatial dependency within the surrounding stations to remove the common-mode error (CME) within a regional network of stations. The main adopted strategies were based on the method of stacking the residual time series of the stations  or the advanced method of weighted stacking (Bock et al., 2000) or spatial weighting (Tian & Shen, 2016). Hence, the spatial analysis of GNSS network data has not been extensively used in early warning applications yet. Therefore, the aim of the current study is to develop a multiple algorithm approach for the analysis of GNSS coordinate time series to detect potential geohazards and/or anomalies. The multiple algorithm approach is based on spatial and temporal analysis methods, which were used to develop two algorithms: the spatial analysis algorithm, which studies each station in relation to its surrounding stations, and the temporal analysis algorithm, which studies each station to its historical record. More specifically, the spatial autoregressive model is used for the spatial analysis algorithm, assuming that the GNSS coordinate time series of a network of stations are spatially dependent; consequently, any station that behaves differently to its surrounding stations will be considered anomalous. Whereas, for the temporal analysis algorithm, an artificial neural network is used to extract the temporal dependency of the GNSS coordinate time series for a single station; therefore, any station that behaves differently to its historical record will be considered anomalous.
In this paper, the main principles of neural networks and spatial autoregressive models are explained, and the development of the multiple algorithm approach is presented. The development of the multiple algorithm aims to enable this approach to detect potential rapid or low-frequency changes in the GNSS coordinate time series, covering the broader spectrum of geohazards, from earthquakes to slow land movements. To evaluate the multiple algorithm approach, it was applied in two different applications: (i) the analysis of the long-term, daily GNSS coordinate time series of the BIGF network in the British Isles, which have small coordinate offsets in long-term time series (greater than a few years) within a relatively small network, and (ii) the 1-Hz GNSS coordinate time series of the GEONET network in Japan, corresponding to the Tohoku-Oki 2011 Mw9.0 earthquake, where short-term (<300 s), large (up to few meters), and spatially correlated coordinate offsets occurred across the GNSS network. The performance of the multiple algorithm approach in detecting anomalies for the two case studies is assessed, and how the type (magnitude and frequency content) of the anomaly affects the detection is discussed.

Temporal Dependency
Temporal dependency is mainly used in stochastic models, where the value of each observation (or measurement) at a time is a function of the values of previous observations within a level of uncertainty (Maybeck, 1982). The form of which these observations are ordered in relation to each other is called a time series (Fu, 2011). Since time series have an explicit order dependence between observations, researchers have used this to analyze the temporal dependency between observation values, where the goal is to develop a mapping function that expresses the relation between the input variables (e.g., past values in the time series) into output variables (e.g., future values in a time series). The structure of this function can either be considered known, parametric models (e.g., autoregressive and moving average models) (Cohen, 2014;Connor et al., 1994), or unknown, nonparametric models (Chen et al., 1997).
The nonlinear autoregressive artificial neural network (NAR-ANN) is one of the most common nonparametric machine learning models used in generic time series prediction (Ahmed et al., 2010;Fu, 2011) and in GNSS applications, such as GNSS/INS integration (Wang et al., 2007) or structural monitoring (Kaloop & Hu, 2015). It defines a complex, nonlinear form of the mapping function, with weights to fit the given data (Zhang, 2019).
Τhe main principle of a NAR-ANN is based on a model which consists of (i) the input observation values, y t − 1 , y t − 2 ,…,y t − p , (ii) the hidden layer(s) containing a number of neurons, and (iii) one output predicted value, b y t , which corresponds to the observed value, y t , and the weights associated with the input values that are to be optimized. Figure 1 represents the flow of the NAR-ANN for the case of two input observation values (p = 2) and a hidden layer with three neurons. Each single neuron works as a computational unit that adds up the weighted inputs and passes them to a so-called "nonlinear activation function" which, in turn, controls the neuron output (i.e., to be passed to the next neuron or not). The sigmoid function, which is considered the most common nonlinear activation function, was used in the temporal analysis algorithm, and it is expressed by the equations: (1) where T i is the output observation value of the function at neuron i, a i is the weighted sum at the neuron i from the input of the neurons in the previous layer, b i is the bias term of the neuron i, and g ij are the weights to be optimized connecting the p inputs to neuron i.
The cost function which represents the error in the estimation (i.e., the difference between the predicted value, b y t , and the observed value, y t ) shows how optimal the weights are given by the equation: where J(θ) is the cost function of the weight θ.
Accordingly, this error will be propagated backwards, in a process called "backpropagation" (Hecht-Nielsen, 1989), by apportioning it to each assigned weight according to the amount each weight is responsible for, in order to minimize this error through optimizing techniques such as the Bayesian regularization training function (Ali, 2009;Jiang et al., 2008) or the Levenberg-Marquardt (Kumar & Rajasekhar, 2017).
Since the performance of the NAR-ANN model depends on so many factors, such as the chosen architecture (i.e., the number of layers and the number of neurons), and the training functions (Zhang, 2019), it is always recommended to examine different NAR-ANN architectures in order to assess their normalized mean square error, which is expressed by the following equation: where N is the number of epochs within the time series data set andb y i is the predicted observation value of y i .
In this research, NAR-ANN is preferred to the parametric models (e.g., autoregressive and moving average models), as it does not require predetermination of the target function to model the time series while it is preferred to Kalman filter as the latter requires some prior models of the signal (e.g., linear or nonlinear) and noise type (e.g., Gaussian noise) to be addressed to optimally estimate the filter parameters. Moreover, NAR-ANN surpasses Kalman filter in the sense of its ability to extract the temporal dependency in both long-term (i.e., periodic signals) and short-term memories, depending on how many previous epochs are included, while Kalman filter practically provides predictions that reflect short-term memory.

Spatial Dependency
With spatial dependency, a variable is a function of space as its observation value in one location is related to the observation values of that variable in other locations (Griffith, 2003;Thanos et al., 2016). This spatial relationship between locations is expressed by the so-called "spatial weight matrix" which represents the spatial influence of the locations of observed values relative to each other. Spatial weights are based on the distances between the locations of these observation values, which is expressed as power distance weights: where w ij is the spatial weight which reflects the spatial influence of an observation at a location j on the desired location i, d ij is the Euclidean distance between location j and location i, and α is the distance exponent value which identifies the "distance decay function." The α takes the values of 2 or 3 (De Mesnard, 2013) to parabolic (Fann et al., 2012) or cubic (Bell, 2006) decay functions, respectively, and d is a cap distance which is used sometimes to exclude values farther than a certain distance, d.
In this study, the spatial regression approach (equation (6)) is preferred to spatial interpolation methods due to its simplicity and ability to be extended using exogenous variables (Manski, 1993). Furthermore, spatial regression models are preferred to the spatial autocorrelation index methods as the latter provides an overall index measuring the degree of spatial influence exerted by something over its neighbors (Griffith, 2003) with no site-specific estimation while the former provides both. In this study, the spatial autoregression (SAR) model is used among the spatial regression methods due to a lack of explanatory variables (i.e., external time series), and it is expressed by the following equation: where Y (n × 1 vector) represents the values of the low-frequency time series of the n stations for a specific time, ρ is the spatial autoregressive parameter to be estimated using maximum likelihood estimate (MLE) (LeSage & Pace, 2009;Elhorst, 2010), W is the hollow n × n spatial weight matrix, which is formed by the spatial weights (see equation (5)) and defines the spatial relation between the n stations, and ε is the spatial residuals when assumed to be white noise. Since W is a hollow matrix, the derived values of the Y vector (y i ) are equal to the weighted sum of the Y values of the rest of the stations (i.e., y j ≠ i ).

Methodology
The GNSS coordinate time series provided as North, East, Up (NEU) displacements were combined to form GNSS coordinate time series of 3-D displacements to be used as inputs for the two developed analysis algorithms ( Figure 2). Any gaps in the time series were allowed to remain, without applying interpolation so as not to alter the time series using artificial data and potentially increase the uncertainties.
The main concept of this approach is to model the normal and expected behavior of the temporal and spatial dependencies of GNSS time series and then to assess the difference between these modeled behaviors against the actual GNSS time series. The temporal and spatial analysis algorithms aim to detect potential anomalies occurring in the high-frequency and low-frequency component of the GNSS coordinate time series, respectively. To model the behavior of each GNSS coordinate time series, they were chronologically split into two data sets: (i) an algorithm training data set, using the beginning of the time series, and (ii) an algorithm testing data set, which follows the algorithm training data set as shown in Figure 3i. The training and testing data sets play different roles in the temporal and spatial analysis algorithms. For the temporal analysis, the algorithm training data set is used to train the NAR-ANNs and to define station-specific temporal analysis thresholds, while for the spatial analysis, it is only used for the determination of station-specific spatial analysis thresholds. For the testing data set, the trained NAR-ANNs, as well as the SAR models, are applied for the temporal and spatial analysis algorithms, respectively, and the predicted value is compared against the real sample to evaluate their difference with respect the corresponding predetermined thresholds (i.e., in the training data set).
The preliminary analysis in the algorithm training data set, common in both the temporal and spatial analysis algorithms, is to detrend by removing any long-term linear trend (i.e., tectonic motion) from each GNSS coordinate time series. Each detrended time series is then analyzed by using the discrete Fourier transform (DFT) to determine a cutoff frequency which separates the low frequencies, expressing the long-period behavior of the GNSS coordinate time series, from the white noise, represented by a flat spectral. Then a low-pass first-order Butterworth filter (Boore et al., 2002;Colombelli et al., 2013) is applied to filter out the frequencies higher than that cutoff frequency of the detrended time series, leading to a less noisy time series representing the expected long-period behavior of the GNSS station. This resultant time series will be referred to as the "low-frequency" time series and will be applied in both the temporal and spatial analyses ( Figure 3ii).

The Temporal Analysis Algorithm
The main hypothesis in the temporal analysis approach is that the GNSS coordinate time series of each station, i, consists of (i) a linear trend, (ii) a low-frequency component (i.e., without white noise) correlated in time, and (iii) a high-frequency component, including white noise and potential anomalies. In the temporal analysis algorithm (Figure 2a), the training algorithm data set was analyzed to train the NAR-ANNs by using the low-frequency time series of each GNSS station and model the corresponding low-frequency behavior. By applying the NAR-ANN model, the next epoch value of the low-frequency time series, y i t , is predicted by using the previous p epochs (Zhang, 2007). The predicted low-frequency time series, b y i (Figure 3a,iii), is subtracted from the detrended time series, x i (Figure 3a,ii), resulting in the temporal residuals, δ i (Figure 3a,v), which include the high-frequency component (i.e., white noise and any high-frequency anomaly) and any abnormalities in the low-frequency components that the NAR-ANN could not predict. The temporal analysis algorithm is expressed by the following equations: where y i t , x i t , and δ i t are the low-frequency, the detrended, and the temporal residual values at a time t and h is the mapping function.
The temporal residuals of the training data set for each GNSS station were used to compute the mean μ and the standard deviation σ and to define the corresponding temporal threshold, as μ ± 3σ range (Figure 3a,v). For the testing data set, if any temporal residual exceeds the corresponding temporal thresholds, it will be considered as a potential event (Figure 3a,v). To distinguish between any site-specific errors (e.g., outliers, multipath, and monumentation) and large-scale events, whenever a potential event is detected, a spatial search (i.e., buffer) is conducted to check if any of the surrounding stations also have potential events at that epoch (Psimoulis, Houlié, Habboub, et al., 2018). If this is the case, the algorithm considers that potential event as a geohazard; otherwise, it is considered as an outlier or site-specific anomaly if the detection remains for several samples (Figure 2a). The search buffer of the spatial search depends on the density of the network and the phenomenon to be monitored/detected. A large search buffer may be less effective as it may include stations that considered to be faraway whereas a small search buffer could isolate stations.
For the purpose of this study, different NAR-ANN architectures were examined, and their normalized mean square error (NMSE) were assessed. For the prediction of each low-frequency sample, the previous two epochs (p = 2) were used in order to make the algorithm more applicable for real-time applications and for GNSS time series with potential data gaps.

The Spatial Analysis Algorithm
In the spatial analysis algorithm (Figure 2b), the spatial weight matrix of the SAR model (equations (5) and (6)) was formed from the geometry of the GNSS network stations; two distance decay exponent values were used: α = 2 and α = 3 to examine the sensitivity of this parameter. For the training algorithm data set, the SAR model was applied independently for each epoch of each GNSS station time series to compute the low-frequency component (Figure 3b,iv), then to compare it to the corresponding GNSS low-frequency time series, and finally to compute the spatial residuals ε. The mean μ and standard deviation σ of the spatial residual time series were computed to define the corresponding spatial thresholds, as μ ± 3σ range, for each GNSS station. Then the same process was followed for the testing algorithm data set, to compute the spatial residuals for each GNSS station for each epoch and compare them against the corresponding spatial thresholds. If any spatial residual exceeds these thresholds, it will be considered as a potential event (Figure 3b,v). Likewise to the temporal analysis, a spatial search is used to distinguish between a geohazard and sitespecific anomalies.

Case Studies and Evaluation
The multiple algorithm approach was applied in two case studies, which cover different aspects: the size of the network (e.g., small vs. large network), the length of the time series (e.g., long-term vs. short-term time series), and the temporal resolution of the time series (daily vs. 1-Hz time series), as well as the dynamics (e.g., high-frequency vs. low-frequency), the size (large vs. small), and the type (spatially correlated vs. site-specific) of the chosen phenomena.
The first case study is based on the analysis of the long-term, daily GPS coordinate time series of the BIGF network of stations in the British Isles and was seen as a relatively small network (22 stations) with longterm (>5 years) daily time series. Due to the fact that these GPS coordinate time series were provided by BIGF with all within-outliers were removed as well as with all documented changes were provided, it was seen to be useful to test the ability of the multiple algorithm to detect these known documented changes and also to examine whether the multiple algorithm can detect any other anomalous behaviors in the time series. The BIGF network data were also used to evaluate the different NAR-ANN architectures.
The second case study is based on the analysis of the 1-Hz GPS coordinate time series of the GEONET stations in Japan for the Tohoku-Oki 2011 Mw9.0 earthquake, which was seen as a relatively large network (>800 stations) with short-term (<1 day), 1-Hz time series with a large and spatially correlated phenomenon to be monitored (e.g., the Tohoku-Oki 2011 Mw9.0 earthquake).
To evaluate the performance of the multiple algorithm approach, the detected anomalies were compared against those that were detected by using a method which is based on considering as anomalies (i) the values exceeding the range μ ± 3σ, where μ and σ are the mean average and standard deviation of a moving window, and (ii) fulfill the spatial search. This method, which will be called the "conventional method," was adopted to develop the RT-SHAKE algorithm and successfully applied in Psimoulis, Houlié, Habboub, et al. (2018), to detect the ground motion using GPS network data for the Tohoku 2011 Mw9.0 earthquake. In that study, it was proved that the detection using GPS data was consistent with that of strong-motion data. Similar approaches to the conventional method have also been adopted by other studies and current early warning systems (Allen & Ziv, 2011;Colombelli et al., 2013). However, in the current study, the conventional method is only used as reference and guidance for the evaluation of the reliability and validity of the detection by using the temporal and spatial analysis algorithms.

Case Study: The BIGF Network of GNSS Stations in the British Isles
The British Isles continuous GNSS Facility (BIGF) holds data and products dating back to 1997 (BIGF, n.d.; Hansen et al., 2012). For some stations, the long-term, daily GNSS coordinate time series have coordinate offsets, which can be due to documented changes (i.e., receiver and/or antenna equipment change) or undocumented changes. The cumulative effect of even small coordinate offsets can significantly alter velocity estimates that are used for the determination of long-term land movement. In the following sections, therefore, the multiple algorithm approach was evaluated by examining its ability to detect any site-specific anomalies due to documented and undocumented changes.

Available Data and Preprocessing
From a processing using Bernese software version 5.2 and GPS only to output Cartesian (XYZ) coordinate estimates in daily SINEX files, data were provided by BIGF as long-term, daily GNSS coordinate time series of North, East, Up (NEU), expressed as displacements with respect to reference coordinates. A total of approximately 150 stations in the British Isles was included in the BIGF processing along with some 200+ IGS stations. For the processing the daily coordinate estimates were created using (i) C13 (CODE repro2/repro_2013) reanalyzed satellite orbit and earth orientation parameter products; (ii) a global network of reference stations and the IGb08 reference frame; (iii) mitigation of firstand higher-order (second and third order and ray bending) ionospheric effects; (iv) a priori modeling of troposphere effects using VMF1G and mitigation using zenith path delay and gradient parameters; (v) I08.ATX models for satellite and receiver antenna phase center offsets and variations; (vi) ambiguities fixed, where possible; and (vii) models for Solid Earth tides, ocean tidal loading and atmospheric tidal loading.
As part of this process, any days that have one or more of the NEU displacements identified as an outlier at this stage are removed; an outlier is identified as a residual greater than three times the weighted RMS of a detrended time series with coordinate offsets due to documented changes accounted for.
All of the available daily GNSS coordinate time series were then considered in order to eliminate any stations with short time series and any stations with a single gap of 10% of the data, which leads to better trend extraction, spectral analysis, and NAR-ANN learning. Furthermore, the stations were preferred to have coordinate offset due to documented changes in the period of the algorithm testing data set for validation purposes. However, since the majority of the documented changes occurred in 2009 (antenna changes), the algorithm training data set was chosen to be from 2000 to 2008 and the algorithm testing data set from 2008 to 2010. As a result, a total of 22 stations of the BIGF network was included in the analysis using the multiple algorithm approach (Figure 4). To examine the performance of the algorithm, it was applied twice (i.e., two runs): (i) after correcting coordinate offsets due to documented changes in the algorithm training data set only; the two analysis algorithms were applied to detect any anomalous behaviors in the algorithm testing data set due to documented changes (supporting information Figure S1); and (ii) after correcting the coordinate offsets

10.1029/2019JB018104
Journal of Geophysical Research: Solid Earth due to documented changes in both the algorithm training and testing data sets; the two analysis algorithms were applied to detect any further anomalous behaviors in the algorithm testing data set, including those due to undocumented changes.
To correct a coordinate offset due to a documented change in the algorithm training data set, the weighted mean method was applied (Montillet et al., 2015;Williams, 2003a). Then, if there are any offsets due to documented changes in the testing data set, the estimated linear trend from the training data set is extended to detrend windows of n days before and after the epoch of the documented change, and then the mean values are subtracted from each other. This technique is based on the conventional averaging method (Blewitt, 1993;Bock et al., 1997;Montillet et al., 2015). The width of these windows, n, is subjective to the data availability and how close the analysis is to a real-time scenario. The wider this window, the better as it then considers including seasonal variations (Montillet et al., 2015). However, a very wide window will be less suitable for real-time analysis as well as introduce more complications in a case where two offsets appear within one window frame. In this research, a 30-day window was found to be optimal for the used data set. Finally, following a spatial search (last step in Figure 1), the anomalous behaviors in this case study were all considered to be site-specific, and none were considered as geohazards, which is due in part to the small number of stations as well as the large distances between them.

NAR-ANN Architectures
NAR-ANN architectures of different (i) number of layers, (ii) number of neurons per layer, and (iii) for two training functions (Bayesian regularization and Levenberg-Marquardt) were assessed for the 22 GPS BIGF time series, to define the optimum NAR-ANN algorithm, using as criterion the normalized mean square error (NMSE; Table 1). The overall performance of the Bayesian regularization is better than Levenberg-Marquardt as the latter gives higher NMSE. Moreover, the difference between the training data set and testing data set is larger using Levenberg-Marquardt which might be due to the fact that Levenberg-Marquardt is more vulnerable to "overfitting problem" than the Bayesian regularization as the latter uses only the nontrivial weights between neurons in the trained NAR-ANN oppose to using all weights in a NAR-ANN Figure 5. Long-term 3-D displacements of the selected 22 stations of the BIGF case study, after correcting for any coordinate offsets due to documented changes in the training data set only, and detected potential events in the testing data set at each station using the spatial (red) and temporal (blue) analysis algorithms. The dashed blue line separates the training from the testing data set. The documented changes are marked with black circles.

10.1029/2019JB018104
Journal of Geophysical Research: Solid Earth (Burden & Winkler, 2008). Overfitting can happen when the NAR-ANN learns the details and noise in the data to the extent that it negatively impacts the performance of the model on new data which in turn leads to higher NMSE in the testing than the training data set.
Furthermore, even though the NMSE values of different architectures using the Bayesian regularization method are almost the same, using a two-layer NAR-ANN slightly enhanced the training process as reflected by the lower NMSE. However, a two-layer network with 10 neurons per layer was preferred as opposed to 20 neurons per layer as the former gives a slightly better result on the testing data set and to reduce the complexity of the NAR-ANN.

The Multiple Algorithm
The application of the multiple algorithms led to the detection of potential events for each station. Figure 5 shows the first run, where the coordinate offsets due to documented changes were corrected in the algorithm training data set and 11 out of 14 offsets due to documented offsets were detected in the testing data set by at least one of the two algorithms. Figure 6 shows the second run (i.e., after correcting for coordinate offsets due to documented changes in both the algorithm training and testing data sets) which illustrates the different nature of the detected anomalies in both the temporal and spatial analysis algorithms. Table 2 presents the number of the detected anomalous behaviors at each station for the two different analysis algorithms (temporal and spatial). However, to check the sensitivity of distance decay exponent value (equation (5)) on the spatial analysis algorithm results, two distance decay exponent values, α = 2 and α = 3, were used. The table also shows the temporal and spatial thresholds values for each station (i.e., derived from the algorithm training data set residuals), the standard deviation of the detrended time series of the algorithm training data set, and the calculated documented offsets in the algorithm testing data set.  Figure 6. Long-term 3-D displacements of the selected 22 stations of the BIGF case study, after correcting for any coordinate offsets due to documented changes, in both training and testing data sets, and detected potential events in the testing data set at each station using the spatial (red) and temporal (blue) analysis algorithms. The dashed blue line separates the training from the testing data set. The documented changes are marked with black circles.

Analysis
In the first run ( Figure 5), 78.5% of the offsets due to documented changes in the algorithm testing data set were detected along with other anomalous behaviors. The magnitude of these documented offsets is shown in Table 2.

The Temporal Analysis Algorithm
The temporal analysis algorithm was able to detect the documented offsets at SCAO, KING, INVR, and MALX stations as the offset values are relatively large (i.e., −10 to −19 mm), even compared to their standard deviation (Table 2), but it missed all the smaller offsets. This is due to the fact that NAR-ANN predicts the next epoch of the low-frequency time series using the previous two epochs; therefore, if the offset appears after a missing observation (i.e., gap), it won't be detected as the first two epochs after the gap will be the inputs. This explains why the temporal analysis could not detect the documented changes such as SHAM and NCAS. Similarly, if an offset happened gradually (i.e., the difference between the two successive epochs is small), the temporal analysis may miss the offset as the NAR-ANN will cope with such offsets especially if they are of small magnitude.
After correcting the coordinate offsets in the algorithm testing data set (i.e., the second run), the temporal analysis algorithm was only able to detect a few anomalies for the stations (i.e., ≤13 days per station).

The Spatial Analysis Algorithm
The spatial analysis algorithm shows a better performance in detecting the documented offsets, by detecting those at BLAP, CARL, DARE, EDIN, INVR, KING, MALX, NCAS, NEWL, and SCAO stations. However, the coordinate offsets due to documented changes were not detected at GLAS, LEED, and SHAM stations. This might be due to the relatively small amplitude of these offsets with respect to the corresponding standard deviation (see GLAS and SHAM; Table 2) or due to the impact of the nearby stations with large offsets which masks the offset of the examined station. For instance, the offset of LEED station with magnitude −3.97 mm

Journal of Geophysical Research: Solid Earth
was not detected, as LEED is surrounded by stations with large negative offsets (e.g., SCAO: −19.05 mm and DARE: −5.73 mm) which masked the anomalous behavior at LEED station.
After correcting for the coordinate offsets due to documented changes in the algorithm testing data set (i.e., the second run), long-term trend problems were still observed at LOWE station ( Figure S2), which justifies the high number of anomalous days in Table 2, as the detected potential events happened at certain times and continue through time. This behavior is a result of either variations in the long-term trend or of a gradual/small coordinate offset due to an undocumented change as the NAR-ANN could not detect this change.
The GNSS stations with temporary anomalous behaviors, however, have a relatively low number of detected anomalous days (e.g., BLAP and IESG stations). It is worthy to compare the spatial analysis results of CAMB station, as it does not have any documented changes in the algorithm testing data set, but there is a difference in the detected events between Figures 5 and 6. This is due to the fact that the offset in the testing data set at NEWL station of Figure 5 was left uncorrected and it affected the estimation of the CAMB station through the spatial regression analysis.
Regarding the impact of the distance decay exponent value, α, it was insignificant, as the number of detected anomalous days was similar for different values of α, mainly due to the relatively small amplitude of the offsets and large distances between the stations. Furthermore, the conventional method failed to detect any of the potential events for the 22 stations which were detected by the temporal and/or spatial analysis algorithm. This is due to the fact that the GNSS coordinate time series were free of outliers and the potential anomalies were of too small a magnitude to be detected by the conventional method (i.e., μ ± 3σ).

Case Study: The GEONET Network of GNSS Stations in Japan in Relation to the Tohoku-Oki 2011 Mw9.0 Earthquake
The Tohoku-Oki Mw9.0 earthquake occurred at 05:46:23 UTC on 11 March 2011, and it was recorded by the GEONET network of GNSS stations operated by the Geospatial Information Authority of Japan (GSI). Numerous studies have used the GNSS network data to estimate the coseismic and transient displacement, in order to model the fault rupture and earthquake characteristics (Ohta et al., 2012;Wright et al., 2012). In this case study, the performance of the multiple algorithm approach was examined for relatively short GNSS coordinate time series (i.e.,~6 hr) with rapid large magnitude coordinate offsets (i.e., up to 5-6 m) in a very short time period (i.e., <300 s).

Available Data and Preprocessing
The 1-Hz GNSS coordinate time series of North, East, Up (NEU) displacements for 847 GNSS stations from the GEONET network, covering the period of 5 hr and 47 min before and 1 hr after the earthquake, were available. The GNSS data, which were processed in Bernese software version 5.2, using PPP processing mode have been successfully applied for the definition of the earthquake characteristics (Psimoulis et al., 2014(Psimoulis et al., , 2015Michel et al., 2017) and in the application of RT-SHAKE for ground motion detection (Psimoulis, Houlié, Habboub, et al., 2018). More details about the GNSS network data processing can be found in Psimoulis et al. (2014Psimoulis et al. ( , 2015.
For the application of the analysis algorithms, the 5 hr and 40 min of the 1-Hz GNSS coordinate time series were used as the algorithm training data set, followed by 12 min for the algorithm testing data set, corresponding to 7 and 5 min before and after the earthquake, respectively.
Due to the fact that the GEONET stations form a dense network and the offsets in the 1-Hz GPS coordinate time series of 3-D displacements are significantly large in magnitude (i.e.,~6 m), apart from using two distance exponent values (i.e., α = 2 and α = 3) in spatial analysis, different cap distances were also introduced. More specifically, the GNSS stations that are located out of the radius defined by the corresponding cap distance are ignored (equation (5)).
In this study, three cases of spatial weight matrix were applied: one with a cap distance of 100 km; one with a cap distance of 200 km; and one using all stations (i.e., with no cap distance applied). Moreover, since the aim of this case study is detecting geohazards, rather than site-specific anomalous behaviors, a spatial search was conducted to filter out site-specific outliers. This search buffer is chosen to be 50 km, as for that distance 96%

10.1029/2019JB018104
Journal of Geophysical Research: Solid Earth of the GEONET stations are surrounded by at least three other stations, which is consistent with Crowell et al. (2009) who used three stations or more to initiate the hypocenter estimation. A cutoff frequency of 0.1 Hz was identified from the DFT analysis and used for the filtering of the low-frequency component.

Results
From the application of the two algorithms (temporal and spatial) along with the three cases of spatial weight matrix, the arrival time of the seismic motion at each GNSS station was detected and analyzed with respect their corresponding distance from the epicenter and compared with the corresponding time arrival detection of the conventional method ( Figure 7). Furthermore, the performance of the algorithms with respect the conventional method was also examined by calculating the percent difference between the corresponding time arrival, as expressed by the following equation (Figure 8): where ATCM is the detected Arrival Time using the Conventional Method, ATDA is the detected Arrival Time using the Developed Algorithm, and n is the number of GNSS stations. The negative difference means that the developed algorithm detected the ground motion later than the conventional method, while the positive difference means that the developed algorithms detected the ground motion earlier than the conventional method. Figure 9 presents the spatial distribution of the GNSS stations where the ground motion due to the seismic waves was firstly detected by the temporal and spatial analyses as well as the conventional method.

Analysis 5.3.1. The Temporal Analysis Algorithm
The temporal analysis algorithm performs similar to the conventional method (Figure 7), as in both cases each station is analyzed separately, with no interaction between the stations in the network. However, the temporal analysis algorithm detects potential events earlier than the conventional method (Figure 8), due Figure 7. Distance from the epicenter versus the detected arrival time of the seismic motion by using the conventional method (black), the temporal analysis algorithm (blue), and the spatial analysis algorithm (red), referenced to the origin earthquake time. The latter was examined for different distance exponents (α = 2 at the left and α = 3 at the right) and different cap distances (100 km, 200 km, or no cap distance). The P-waves velocity was assumed about 6 km/s (Psimoulis, Houlié, Habboub, et al., 2018).

10.1029/2019JB018104
Journal of Geophysical Research: Solid Earth to its higher sensitivity, through the NAR-ANN modeling and the analysis of the high-frequency GNSS coordinate time series. This led the temporal analysis algorithm to obtain an elapsed time of the first detection of the geohazard 3 s earlier than the conventional method for the closest GNSS stations to the earthquake epicenter ( Figure 9a). It is also expected that the temporal analysis algorithm would have performed better than the conventional method if the GPS coordinate time series had been of longer time period, where the impact of periodic signals and linear trends would be more evident and their modeling would benefit even more by the NAR-ANN. In Figure 9, it is also shown that the temporal analysis algorithm led to the detection of the ground motion in a broader area than that of the spatial analysis algorithm. This may be due to the averaging in the spatial regression through using the spatial weight matrix to compare the relative behaviors. Thus, the estimation at a GNSS station by the spatial regression model is formed as a weighted average of the surrounding stations; for a GNSS station with larger displacement to its neighbors, the estimated displacement will be smaller than the true displacement of the GNSS station.

The Spatial Analysis Algorithm
The application of the spatial analysis algorithm with α = 2 for the spatial weight matrix and all the available GNSS stations (i.e., no cap distance) leads to false detection of the ground motion (Figure 7), as the corresponding arrival time corresponds to a propagation velocity significantly faster than that of the P-wave (i.e., typical value of 6 km/s), with the latter being also estimated using strong-motion data in Psimoulis, Houlié, Habboub, et al. (2018). Hence, this configuration was not further analyzed; the large displacement (reaching up to 6 m) at a station close to the epicenter affected the estimates at all stations, especially when the decay function is not steep enough to dissipate this propagation (i.e., comparing with using α = 3).
For all of the other combinations of α and cap distance, the overall performance of the spatial analysis algorithm seems similar, with the general trend being that the detected arrival time of the seismic motion is slower than the P-wave and slower than the conventional method (Figure 7). This is also reflected numerically (Figure 8) in the percent difference values which are all negative, meaning that the spatial analysis algorithm detected the ground motion later than the conventional method. From all the examined configuration of the spatial analysis algorithm (combinations of α and cap distance), the one with α = 2 and 200 km cap distance seems to be the closest to the conventional method ( Figure 8).
Furthermore, it is observed that, for a given cap distance for the spatial regression, the smaller the α used for the spatial weight matrix, the more the impact of the seismic motion from far stations, which makes the algorithm more susceptible to false detection due to large displacements from distant stations, as in the case of a large earthquake. Likewise for a given α, the larger the cap distance, the more susceptible the spatial analysis algorithm is to false detection due to the impact from distant stations. Therefore, the combination of the cap distance and α should be optimized for the spatial analysis algorithm, according to the geometry of the GNSS network and the expected coordinate offset due to the monitored geohazard, in order to reduce the occurrence of displacements in the spatial residual time series, which are artifacts created by large offsets of distant GNSS stations. However, at the beginning of the detection of the geohazard, the characteristics of the detection (stations and time detection) were very similar regardless of the values of α and cap distance, indicating that these parameters affect the detection of the geohazard for distant stations (Figure 9).
In Figure 9, it is also shown that potential events are firstly detected at stations 171, 176, and 550 for the temporal analysis algorithm and 550 and 1179 for the spatial analysis algorithm (hollow circles). However, due to the limited spatial buffer (i.e., 50 km) of the spatial search, they were flagged but not considered as geohazards, which indicates the potential impact of the spatial buffer on the detection of the geohazard.

Discussion and Conclusions
The multiple algorithm approach was applied in two different GNSS networks of different data rates and density of GNSS stations (see GEONET and BIGF characteristics) to detect different type of anomalies in the GNSS coordinate time series: (i) rapid large magnitude changes in the time series (Tohoku Mw9.0 Figure 9. Spatial distribution of the GNSS stations where the ground motion was detected and flagged as geohazard (marked by filled circles) as well as the GNSS stations with detected motion flagged as potential event (marked by hollow circles) through a time window of 6 s using the (a) temporal (blue) and (b) spatial (red) analysis algorithms, along with (c) the conventional method (black). The distribution of the GNSS stations for the spatial analysis algorithm expresses the outcome of the detection regardless the value of α or cap distance. earthquake) and (ii) slow and small magnitude changes in the time series (BIGF network). The two presented case studies showed that the multiple algorithm approach can be applied to detect displacement related to geohazards and anomalies, of different magnitude and frequency content, thanks to the complementary use of the temporal and spatial analysis.
More specifically, the spatial analysis algorithm seemed to be sensitive in detecting slow and small magnitude changes in the GNSS coordinate time series for a station or a small cluster of stations, based on the comparison with surrounding stations, as in the case of BIGF stations. This makes the spatial analysis approach effective in detecting site-specific anomalies, due to GNSS station problems (e.g., unstable monumentation; Hackl et al., 2011) or regional anomalies (i.e., vertical land movement; Wöppelmann and Marcos, 2016). Also, it is less vulnerable to outliers, since outliers in the high-frequencies content are filtered, leading to narrow-spread spatial residuals which is effective in detecting small long-period offsets. Characteristic is the example of LOWE station, where it was possible to detect an anomaly of magnitude 7 mm slowly developed in a period of 1 year ( Figure S2), which could represent potential malfunction of a station (i.e., corrosion of GNSS antenna) or geophysical processes (e.g., postseismic relaxation; Savage et al., 2005;Broerse et al., 2015). However, for rapid large magnitude changes in the GNSS coordinate time series, as in the case of Tohoku-Oki 2011 earthquake, the spatial algorithm proved ineffective, as large anomalies in specific stations bias the GNSS coordinate time series of stations which do not appear to have an anomaly at the corresponding epoch; the latter is known as "spatial spillovers" (LeSage & Pace, 2009;Elhorst, 2010) and can result in a false triggering of the detection of the algorithm at those stations. The adjustment of the parameters (cap distance and α) of the spatial analysis algorithm limits these biases, but still, the spatial analysis algorithm proved less effective than the temporal analysis algorithm for rapid large magnitude changes in the GNSS coordinate time series (e.g., earthquakes).
In contrast, the temporal analysis algorithm was proved to be more effective in the application of the Tohoku-Oki 2011 earthquake (i.e., rapid large magnitude changes in the GNSS coordinate time series) as the temporal residuals are dominated by high-frequency components, which contain the rapid displacement signal of the ground motion. However, the high-frequency components of the temporal residuals are also characterized by white noise resulting in a large threshold (i.e., μ ± 3σ), which makes the temporal analysis (i) less effective in detecting low-frequency offsets of small magnitude and (ii) more susceptible to outliers. The latter makes necessary the spatial check criterion in order to avoid potential outliers and false alarms. Based on the examined case studies, an anomaly detected primarily by the temporal analysis algorithm indicates a relatively rapidly developed phenomenon (high-frequency content), which if it is not observed in surrounding stations, should correspond to a site-specific effect.
Thus, the complementary function of the temporal and spatial analysis algorithms, which could be implemented even for real-time applications, reveals their potential to support the operation of existing GNSS networks and contribute either to the reliable operation of the GNSS networks, by detecting site-specific anomalies at GNSS stations and/or inconsistencies across the GNSS network, or to supplement existing early warning systems (e.g., seismic sensor networks) for the detection of ground motion (e.g., earthquake and landslide). The main requirement of the multiple algorithm is a training data set to use as reference to train the NAR-ANN and define the thresholds, and then the multiple algorithm can be applied in real time, by predicting epoch by epoch and comparing it with the corresponding real GNSS coordinate value.
Even though our multiple algorithm approach proved reliable in detecting small amplitude anomalies in GNSS time series (i.e., 3-to 4-mm level of anomalies in BIGF time series), further investigation is needed to enhance the performance of the developed method and detect anomalies close to the noise floor of the GNSS time series (i.e., 1-2 mm). For example, it is needed (i) to assess the performance (i.e., robustness and sensitivity) of more advanced techniques for the temporal analysis (e.g., Long Short-Term Memory Neural Network model) and the spatial analysis (e.g., Manski model) and (ii) to adopt more sophisticated techniques to model the noise level and the long-period characteristics of the GNSS time series (e.g., Generalized Gauss Markov; Langbein, 2004;He et al., 2020;Langbein & Svarc, 2019) and to define the threshold criteria (e.g., Pierce's criterion and modified Thompson tau test). Furthermore, the multiple algorithm approach should also be applied for different type of geohazards (e.g., landslide and land subsidence), different data sets (e.g., ionosphere and troposphere), and different sizes of monitoring networks in order to evaluate its sensitivity for different data sets.

10.1029/2019JB018104
Journal of Geophysical Research: Solid Earth