PCA is a linear dimensionality reduction technique that seeks projection of the data into uncorrelated (and orthogonal) directions of highest variability (Chatfield and Collins, 1980; Stone, 2004). It computes the principal components (PCs), which are the basis vectors of directions in decreasing order of variability. Thus, the first PC provides a basis vector for the direction of the highest detected variability. Then, the second PC provides the basis vector for the next direction orthogonal to the first PC, and so on follows the sub-component decomposition. Computation of PCs involves estimating the covariance matrix of the data, its eigenvalue decomposition, sorting of eigenvectors in the decreasing order of eigenvalues, and finally, a projection of the data into the new basis. Considering that our measured data consist of daily time series of a GNSS network, we can define X as the data matrix of mxn, where m is the number of stations of the network, and n the number of days of the time series. This matrix X contains only one of the three spatial components of the GNSS data (East, North or Up).

Thus, X _ij corresponds to the value of X for the i-th station (i = 1…m) and the j-th day (j = 1…n). Accordingly, we can define the mx1 vector x _j as the observations on the j-th day for all the m stations. To obtain the covariance matrix of X we must first compute the mean vector of the data as: x ¯ = 1 n ∑ j = 1 n x j (1) and the mxm covariance matrix C is defined by: C = 1 n − 1 ∑ j = 1 n ( x j − x ¯ ) ( x j − x ¯ ) T (2)

Let λ ₁, λ…λ _m, be the eigenvalues of C ordered so that λ ₁≥ λ ₂≥…≥ λ _m, and let v ₁, v ₂…v _m be the correspondingly ordered eigenvectors. We can define the mxm eigenvector matrix W as: W = [ v 1 T ⋮ v m T ] (3) then the spatial PCs for the j-th day are given by: y j = W ⋅ ( x j − x ¯ ) (4) where y _j is a mx1 vector representing the value of the PC for the j-th day at each GNSS station. The resulting PCs can be expressed for each day j = 1…n as a mxn matrix Y, where the first row corresponds to the first PC, the second row to the second PC, and so on, for each day with data in the GNSS time series. PCA is usually used to reduce the dimensionality of the data as well as for filtering noise, by discarding the PCs associated with the lowest eigenvalues. This can be done choosing a dimension d<m and projecting the data onto v ₁, v ₂…v _d , giving the approximation: W d = [ v 1 T ⋮ v d T ] ⇒ y ^ k = W d ⋅ ( x k − x ¯ ) (5) for k = 1…n, and y ^ k is a dx1 vector, resulting in a dxn matrix Y, where the rows correspond to the d first PCs. To partially recover the original data we can calculate: x ^ k = W d T y ^ k + x ¯ (6) where X ^ is the recovered filtered data, considering only d PCs. In this work, we use the estimated PCs as possible sources of transient information. Also, PCA was used to detrend the data, that is, recovering it from all PCs except the first component, which in most cases contains the secular trend components of GNSS time series, as observed in our analysis of Chilean GNSS data.

The PCA technique has been used in many fields, including GNSS time-series analysis (Dong et al., 2006; Kositsky and Avouac, 2010). However, PCA finds signals which are Gaussian and uncorrelated, which means that if the transients we are looking for are not Gaussians in nature, it is most likely that PCA will not be able to represent the transient with only one PC. This problem can be handled using a complementary technique, the ICA, capable of decomposing the original data into non Gaussian signals (Gualandi et al., 2016). In simpler words, while PCA decomposes the GNSS time series into a set of orthogonal spatial and temporal components, ICA does so into components that are statistically independent, without requiring mutual orthogonality. Generally, it is assumed that each measured signal depends on several distinct source signals. Moreover, each measured signal is essentially a linear mixture of these source signals. In many cases, these source signals are of primary interest, but they are interlaced within the set of measured signals, or signal mixtures. The problem of unmixing signals is known as blind source separation (BSS), and ICA (Comon, 1994; Hyvärinen and Oja, 2000; Stone, 2004) is a specific method for performing BSS. The solution of a BSS problem is a set of source signals that explain, using a linear mixture, the available measures. This problem refers to when neither the source signals nor the mixing structure are known, and can be written in matrix form as: X = A ⋅ S (7) where X is the matrix of measured signals x ₁…x _m; S is the matrix of source signals s ₁…s _m; and A the mixing matrix of the sources. Assuming that matrix A is invertible, we can write W = A ⁻¹, then Eq. 7 can be rewritten as: S = W ⋅ X (8)

The goal is then to estimate the coefficients W _ij from matrix W. A solution to this problem can be found by considering the statistical independence of the source signals. In this work, we apply the algorithm called FastICA (Hyvärinen and Oja, 2000), in which the key to the solution is the use of the central limit theorem. This theorem establishes that the distribution of a sum of independent random variables tends toward a Gaussian distribution. To estimate one of the independent components (ICs), let us consider a linear combination: y = w T ⋅ X = ∑ w i ⋅ x i (9) where w is a vector to be determined. If we make a change of variables z = A ^T ·w then: y = w T ⋅ X = w T ⋅ A ⋅ S = z T ⋅ S = ∑ z i ⋅ s i (10)

Thus, y is a linear combination of the sources s _i . Using the central limit theorem, we can observe that y (or z ^T ·S) is more Gaussian than any s _i and it becomes least Gaussian if it equals one of the s _i . Therefore, the task of ICA is to find a vector w that maximize the nongaussianity of y = w ^T ·X. This vector w will necessarily correspond to a vector z with only one nonzero component, that is, y = w ^T ·X= z ^T S corresponds to one of the independent components.

The main difference between ICA and PCA is that PCA decomposes a set of signal mixtures into a set of uncorrelated (orthogonal) signals, whereas ICA decomposes a set of signal mixtures into a set of statistically independent signals. Assuming that we are looking for transients that can be gaussians or non gaussians, we used in this work both techniques, in order to find the best performance for different kinds of transients, first with synthetic GNSS data and then with real data from different regions of Chile.

This detection algorithm searches for the principal or independent component - from PCA or ICA, respectively -that best represents the transient event that potentially exists within the GNSS time series. This search is performed by means of the temporal correlation between each time series and each estimated PC or IC based on the following hypothesis: the information of a transient, if it exists, is contained in at least one of the estimated PCA or ICA components. Therefore, the highest correlations between the time series and the PCA or ICA components are obtained when these components contain information of a transient event. A block diagram of the Correlation PCA/ICA-based algorithm is shown in Figure 2, which is explained in detail hereafter. The input data is composed by the daily GNSS network stations time series. This data is structured in a matrix X, with rows corresponding to the different stations in the network and the columns to the daily samples that span the full extents of the time series. The pre-processing stage (Stage 1 in Figure 2) is to apply PCA, calculate the first d PCs and recover the time series from these d PCs except the first one, that contains the secular trend or station velocity, obtaining the detrended time series for subsequent correlation with the PCs or ICs. Then, we apply ICA to obtain the corresponding ICs to correlate with the time series. As well as with PCA, we define d sources for the estimation in the ICA. We note that the correlations with the time series must be calculated using either PCA or ICA, so to obtain the results from both techniques, we apply the method twice, for PCA and for ICA, respectively.

With all the signals preprocessed and estimated, we calculate the correlation coefficient between them (Stage 2 in Figure 2). Let X ^ be the mxn matrix whose rows correspond to the detrended GNSS time series, and columns to each day of the given time span. Also, we can define a dxn matrix S, which contains the d PCs or ICs over the n days of the time series. The correlation will be calculated for sliding windows of length R days, so that: ρ k l ( p , q ) = 1 R − 1 ∑ j = 1 R ( S k p ( j ) − μ ( S k p ) σ ( S k p ) ) ( X ^ l q ( j ) − μ ( X ^ l q ) σ ( X ^ l q ) ) (11) where ρ _kl (p,q) is the correlation coefficient between time windows S _kp and X ^ _lq . Time indexes p = 1…n-R and q = 1…n-R, thus, we have a (n-R)x(n-R) matrix ρ _kl of correlations for every kl combination. The window S _kp corresponds to the k-th PC or IC (k = 1…d) in time window from p to p + R; and X ^ _lq to the l-th station time serie (l = 1…m) in time window from q to q + R. Also, μ and σ are the mean and the standard deviations of time windows S _kp and X ^ _lq . To summarize all this new data, a new dxm matrix H _ρ is generated (Stage 3 in Figure 2), containing only the highest correlations from every ρ _kl . Also, two new dxm matrices are obtained: a dxm array containing the p value (time index) for the highest correlation from every ρ _kl , and another dxm array containing the q value (time index) for the highest correlation. With this information, a search for transient information is performed for each PC or IC analyzed, that is, for each row of matrix H _ρ (Stage 4 in Figure 2).

The search procedure consists of:

1. Search for the maximum value of the matrix row and set a threshold (relative to this maximum value) to consider all the correlation values above this threshold in the row. All stations with correlations below it are discarded.

The resulting cluster is composed of stations that are spatially consistent, with a time index that corresponds approximately to the onset time of the transient, and with a duration related to the sliding window length R. This method was tested on synthetic GNSS data and then applied on real and documented slow slip events in Chile.

To evaluate performance of the method, different values of R must be applied, in order to obtain a set of different detection results that can be evaluated through performance indexes. We use known a priori information that includes the onset time of the transient and the group of stations containing significant transient information in the region of the simulated slip. Thus, we can evaluated the error between the real onset time and the detected onset time, as well as the inclusion of the detected stations in the predefined group of stations. A normalized time error index is defined by: T e r r = | t i − t i d | n (12) and two known indexes are also used as spatial performance evaluation: Pr e c i s i o n = T p T p + F p         Re c a l l = T p T p + F n (13) where t _i is the onset time, t _id is the onset time detected and n the number of days of the time series. The T _err index evaluates the error in the transient onset time detection, normalized by the total number of days in the time series. Precision and Recall are used in binary classification algorithms to evaluate the performance of the results. Since we are developing a method to detect (or not detect) transients in time and space, we can see the problem as a binary classification, so it is possible to evaluate the performance of our method using these two indexes. Here, T _p is the number of true positive cases, that is, the number of stations that contain transient information and, consistently, are correctly detected; T _n is the number of true negative cases, namely, the number of stations that do not contain transient information and are not detected. Similarly, F _p is the number of false positive cases, in other words, the number of stations that do not contain transient information but are detected; and F _n the number of false negative cases, that is, the number of stations that contain transient information but are not detected. Thus, Precision evaluates the ability to detect true positive cases over all the positive cases detected, and Recall evaluates the ability to detect true positive cases over all the real positive cases.

We use synthetic and observed GNSS data in order to develop and test the PICCA algorithm. Synthetic GNSS time-series were generated based on the model proposed by Carr Agnew (2013). This model defines the position for the i-th station (i = 1…m) and the j-th day (j = 1…n) as: X i ( t j ) = [ ( t j − t 0 ) v i + r i ( t j ) + y i ( t j ) + c i ( t j ) + ∑ l = 1 L j α l i H ( t j − t l i ) + ∑ ε = 1 Ε G i ε s ε ( t j ) ] ⋅ D i ( t j ) (14) where v is the velocity, r is the individual noise per station, y the seasonal component and c is the common mode noise. The first sum corresponds to the jumps or offsets in the data, however, for this work they were not considered in the analysis to simplify the complexity in the use of the proposed method. The second sum corresponds to tectonic displacements due to slips at Ε different fault patches. We generated Green’s functions for the upper boundary of the slab using the Slab 2.0 geometry (Hayes et al., 2018). We use a geometry in Cartesian coordinates for the fault extending from 18°S to 45°S. This may distort the geometry due to the sphericity of the earth. However, the modeled slow earthquakes extend over a fault zone ∼300 km long, so at this scale the use of Cartesian coordinates is appropriate. The mesh was discretized into 1,667 triangles and the surface displacements (at each station in the continuous GNSS network) were estimated using the TDdispHS, which is a triangular dislocation analytical code (Nikkhoo and Walter, 2015). These static displacements were temporally distributed to simulate slow earthquake events of varying duration and magnitude. The function used for this temporal distribution is given by: s ( τ j ) = ( 1 + τ j a ) b − 1 (15) where we choose a = 2, b = −0.5 and the time τ _j varies from j = 1…h, where h is the transient duration in days. In these synthetic tests, only one transient event occurring in a given time window was used and its daily offsets were added to the other components (Eq. 14) to generate the time series. Finally, D corresponds to a binary function that indicates whether i-th data is available or not (1 or 0 respectively). In order to simplify the application of the PCA and ICA separation methods, this function was considered to always have a value of 1. Thus, the model proposed in this work for the synthetic series is reduced to the following form: X i ( t j ) = ( t j − t 0 ) v i + r i ( t j ) + y i ( t j ) + c i ( t j ) + G i s i ( t j ) (16)

An example of this model is shown in Figure 3, where a synthetic signal is generated from the five components presented in Eq. 16. The sets of synthetic data generated in this work are based on 150 GNSS stations (Figure 4) from the Nevada Geodetic Laboratory, and the three components of displacements (East, North and Up) were used. We created three synthetic datasets of daily observations, in which we used different configurations for the duration and magnitude of a slow earthquake event. The first data set (Transient No. 1) consists of a 3 years span (1,095 days) time series, with a tectonic slip of 0.05 m and a magnitude of M _W 7.2, distributed in 150 days and starting on day 501 (Figures 3, 4). The second data set (Transient No. 2), with time series of 3 years span (1,095 days), a fault slip of 0.03 m and a magnitude of M _W 7.0, distributed in 14 days and starting on day 1,082, simulating a precursor activity of an earthquake; and in the third data set (Transient No. 3), time series of 1-year span (365 days), with a slip of 0.03 m and a magnitude of M _W 7.0, distributed in 3 days and starting on day 363, simulating a “fast” precursor activity. The complexity of the velocity field together with the noise and seasonal motions causes the 150 days-long transient signal to be hidden in the 3 years time window (Figure 4C), even for stations close to the transient event. So, even an obvious transient might not be apparent with a cumulative displacement map.

We applied the Maximum Correlation PCA/ICA-based algorithm to analyze the three synthetic datasets. For both PCA and ICA, 10 components were estimated (Supplementary Figure S2 in the Supplementary Material). We estimated the correlation matrix ρ _kl (Eq. 11) between each PC or IC time window and the original time series for each station. Figure 5 shows an example of the matrix ρ _kl of correlations for 2 cases of kl combination, stations IMCH and MAUL, using synthetic data. In Figure 5A, the highest correlation is obtained from PC 2 and the station IMCH time series window that starts on day 487. However, in Figure 5B, there are no matches between time indexes for the highest correlations, which also are lower than the values found in Figure 5A. The application of the method is shown in Supplementary Figures S3A, S3B in the Supplementary Material. In this figure, results for Transient No. 1, stage by stage, are shown based on the block diagram of Figure 2.

To analyze the results of the PICCA method on synthetic data, different values of R were considered, from 10 to 200 days, with a step of 10 days, The minimum number of days of R was selected considering the frequency of noise in the data. In Figure 6, two plots are shown for each case: the left plot is a two y-axis graph, where the left y-axis correspond to the values of the three indexes, Precision, Recall and Time error, from 0 to 1, while the right y-axis correspond to the component number, from 1 to 10. The values of the indices were calculated from the spatial and temporal results, and the plotted component is the resulting component for each value of R. The right plot corresponds to an histogram of these components for all the values of R. Table 1 shows the components from PCA and ICA that contain observable transient information for every case, showing the ability of the method to choose the component that better represents the transient analyzed.

Transient No. 1 (Figure 6A) presents a maximum Precision and Recall for values of R of 90 and 100 days, and between 130 and 170 days. The histogram shows that PC 2 is the most selected component over the values of R, and is consistent to the Precision values, demonstrating that the component selection is correct. Also, we can observe that the time error is very low for that component selection. In Transient No. 2 the results for PCA are more ambiguous (Figure 6B), since PC 2 (incorrect component) had as many matches as PC 5 (correct component). However, the method is still capable of obtaining correct results for low values of R, where performance indexes present good results. For ICA (Figure 6C), the results are more clear, the most selected component is IC 6 and the best performances are for an R between 20 and 70 days. For Transient No. 3, similar results as those from transient No. 2 are obtained: PCA-based method (Figure 6D) results in more matches with the incorrect component (PC 3), while in ICA-based method (Figure 6E) the most selected component is IC 2, in a correct match, specially for lower values of R, consistently with the performance indexes. This suggests that PCA and ICA have different performance depending on the duration of the transient event.

From the performance evaluation, we chose the results with best performance and summarized them in Table 2. In this table we present details of transient event characteristics (duration, time analyzed) and detection results such as event onset time, duration and stations detected, for the PCA and ICA methods. For Transient No. 1, the correlations in matrix H _ρ for every PC and every station (Figure 7A) show that the time series of stations IMCH and PECL have the highest correlations with the PC 2, while time series from station MAUL has low correlation (as shown in the time plot of Figure 7B). It is noted that from matrix H _ρ it is quite difficult to visually establish the PC that effectively contains the transient. However, matrix H _ρ is calculated only in Stage 3 of the method (see Figure 2), with 3 processing stages remaining to obtain the final result, namely, the clustering stage to detect the groups of similar time indexes, the component selection stage and the final spatial consistency stage (Supplementary Figure S3B in the Supplementary Material shows an example where these remaining stages can be more clearly understood). Time series with high correlation (Figure 7C) correspond to stations closer to the area where the simulated slip has its maximum magnitude (Figure 4B). The same analysis using ICA fails to detect correlations that would indicate the stations close to the transient event as well as its onset time. For Transient No. 2, results for the ICA-based method are shown in Figure 8. Clearly, the highest correlations are obtained from IC 6 and time series from stations IMCH, PECL, and others. Time series from station QTAY has low correlation, which is evident in the time plot in Figure 8B. The highest correlation was for the window that starts on day 1,073. In Figure 8C, time series with highest correlation clearly corresponds to the stations in the area where the simulated slip occurs. In this dataset, the PCA analysis provides similar results to ICA. The onset time of the transient event is found as the beginning of the window with the highest correlation. However, to estimate the duration of the event it is necessary to find the optimal width of the time window R, which cannot be estimated a priori. Results for Transient No. 3 are shown in Supplementary Figure S4 in Supplementary Material.

Once our method was calibrated and tested with synthetic data, we used GNSS time series that recorded the three main transient events in Chile to assess the ability of the PICCA algorithm to detect these events automatically. In doing so, we select three time-windows and locations: Iquique area, spanning between February 4, 2013 to April 1, 2014. We focused on detecting precursory transient movements in the days before the April 1, 2014 earthquakes, which is the most prominent transient signal prior to this event (Schurr et al., 2014); The second location is Valparaíso, spanning between March 25, 2016 to April 24, 2017. We attempt to detect the intense precursory activity documented 2 days before (e.g., Ruiz et al., 2017; Ruiz et al., 2018; Caballero et al., 2021) the M _W 6.9 earthquake in Valparaíso. The third location is Copiapó, spanning between January 1, 2014 to September 12, 2015, where a long-term transient was detected between 2014 and 2016 (Klein et al., 2018). In the analyses, we use only the time series measured in the EW direction. Only in the case of the Copiapo event, we use the vertical signal. However, within the results shown in the following section, the inferred values for the directions not analyzed (i.e., N-S for Iquique and Valparaiso, and additionally E-W for Copiapo) are extrapolated from our analysis, taking the time windows detected by PICCA and extracting the detrended displacement in that period.

As well as with synthetic data, we apply the method for different values of R, as we can see in Figure 9. Performance indexes and the number of the resulting component selected were plotted for all the values of R, and histograms of the selected components are shown for each data set. Also, Table 1 is used to validate the results.

For the Iquique data set, the PCA-based method is not capable of generating consistent results (see Figure 9A). The component selection failed for most of the values of R, with only one correct result obtained for R = 20 days, which is consistent with the transient duration. Although this result is correct and consistent with the transient duration, the method is not reliable using PCA. On the other hand, the ICA-based method presents a high consistency between the correctly selected component (IC 2) and the performance indexes (Figure 9B). The results have a high performance for values of R from 10 to 100 days, similarly to the results from the ICA-based method for synthetic data of 15 days slip (Transient No. 2). In the second data set, corresponding to Valparaíso, the ICA-based method does not perform stable in function of R (Figure 9C). Only two values of R were able to select the correct component (IC 2), with only one of these values, R = 20 days, with high performance. We also can observe that the most selected component was IC 10, for larger values of R, however, the onset time estimated differs substantially from the real onset time of the transient. Finally, for the Copiapó event, most of the component selection were correct (PC 3), obtaining a high performance, specifically for higher values of R, consistently with the several-months duration of the transient (Figure 9D).

From these performance results, we summarize the best results in Table 3. For the Iquique region (see all the principal and independent components from Iquique dataset in Supplementary Figure S5 in the Supplementary Material), we can observe that the highest correlations are given by IC 2 and time series from stations ATJN, CGTC and IQQE (Figure 10A). The corresponding time series (red lines in Figure 10B) show the detection window, from day 392 (March 2, 2014), including the precursor activity that starts on day 406 (March 16, 2014). As we can note in Figure 10C, these time series with high correlation correspond to the stations closer to the 2014 earthquake rupture zone and possibly where the pre-earthquake slow slip occurred. GNSS-derived cumulative displacements show a consistent trenchward pattern (with magnitudes>10 mm) for the transient event window. The direction of motion associated with this SEE is consistent with the distribution of seismicity (extracted from NEIC, Masse and Needham, 1989) in the time window analyzed, indicating that this SEE is related to this seismic activity (Bedford et al., 2015).

For the Valparaíso region, Figure 11A shows high correlations between IC 2 and time series from stations BN05, VALN and QTAY. Because of the size R of the window in relation to the duration of the transient, the onset time detected is far before the real onset time, however, the transient is included in the window, thus, correctly detected. The cumulative movement in the time window reaches 5 mm, with a consistent pattern of motions towards a zone with high seismicity activity in the analyzed time period. On the third real dataset, in the Copiapó region, the Up component was analyzed, following the finding of Klein et al. (2018). This transient is a long-term deformation (possibly several months) that has not an accurate onset time, so we established that the day 200 of the time series was the onset time, only as a reference for the calculation of the time error and compare results. In Figure 12A we can observe that the highest correlations are given by PC 3 and time series from stations BN03 and COPO. Despite the large amount of noise and gap in the data, as we can see in Figure 12B, the method was capable of detecting 2 stations from this region with displacements consistent with the previous studies (Figure 12C). During the analyzed period, no seismic activity was recorded that could be related to the occurrence of an SSE. The horizontal displacements are of very low magnitude, and it remains to be seen whether this is a phenomenon related to a slow earthquake or an unusual deformation affecting the vertical component. Whatever its origin, our method is able to detect this signal, which affected the Copiapó area.

Summarizing, similar results can be observed in both synthetic and real data. For long-term transients (150 days or more), the PCA-based method was able to find the components with displacement information, while the ICA-based method performs better for shorter duration events (15 days for example).