Singular points of the tremor of the earth's surface

1. Introduction

The thin structure of the time series of displacements of the earth’s surface, measured by means of GPS, is the subject of research by a large number of specialists. The maximum likelihood approach for estimating parameters of GPS time series models was used in [1,2,3,4]. The structure of the noise and estimates of errors were considered in [5,6,7]. Parametric models of GPS time series for the analysis of phase correlations were used in [8]. Seasonal effects in GPS tine series connected with hydrological loading and tectonic movements of the crust were studied in [9,10,11,12,13,14,15,16] by different methods including multidimensional statistic approach. Influence of ionospheric and tropospheric delays on the sensitivity of the GPS solutions were considered in [17]. Hypotheses about origin of “abnormal” harmonics in GPS time series were regarded in [18].

The coherence of the earth’s surface tremor was analyzed in [19,20]. In [21,22], fractal analysis of high-frequency GPS time series was used to assess seismic hazard in Japan and California.

This article presents a further development of the methods for analyzing the earth’s surface tremor, proposed in [23], and analyzes the high-frequency noise structure in the time series of the GPS station network in the Western United States during the time interval 2009-2022. The US West is one of the most seismically active regions. A dense network of GPS stations makes it possible to study in detail the features of the earth’s surface oscillations and try to compare them with global seismic activity. Note that earlier in [24], a study was made of the properties of seismic noise in Southern California in comparison with the seismic regime.

2. Materials

GPS ground displacement data taken from the Nevada Geodetic Laboratory website at: http://geodesy.unr.edu/gps_timeseries/tenv3/IGS14/.

These data represent displacements of the earth’s surface in three orthogonal directions, which we will further denote as E, N and U (Up) for displacements in the west-east, north-south and vertical directions, with a time step of 1 day [25].

In what follows, we will consider a regular grid of 100×150 nodes, covering an area with latitude from 32.5°N to 51°N and longitude from 115°W to 128.5°W. Figure 1(a) shows the positions of GPS points for measuring displacements of the earth’s surface. Figure 1(b) shows the distribution density of GPS stations obtained using the Gaussian kernel estimate [26] at the nodes of this regular grid with a smoothing radius of 0.25° (about 28 km).

We will consider estimates of the properties of GPS time series in a moving time window 365 days long with a shift of 7 days. For each time window, operable stations are determined from the condition that the number of gaps in the considered time window does not exceed 5% of the time window length, that is, that this number is not more than 18 values. Missing values are filled in with artificial data, which are found from the value to the left and right of the missing time interval of the same length as the length of the gap [19]. Figure 2 shows a graph of the number of working stations within a sliding time window of 365 days for 14 years of observations, 2009-2022.

3. Description of the time series statistics used

For each station in each time window, 4 surface tremor statistics are computed for the 3 ground displacement components.

Spectral exponent. Let $x(t), t=1,,,N$ be a sample of the GPS time series within a window of sample lengths $N$. Denote by $S_{xx}({\omega }_j)$ the estimate of the signal $x(t)$ power spectrum ${\omega }_j$, - a sequence of discrete values of frequencies with a step $\Delta \omega$, $\Delta \omega =2\pi /(N_f\cdot \Delta s)$, $j=1,\mathrm{...},N_f$, $\Delta s$, - a time step. Here $N_f$ is the minimum length, which is greater than or equal to the sample length $N$: $N_f=2^n\ge N$, the index $j$ in formula (6) varies up to $N_f/2$ due to the symmetry of the power spectrum of real signals with respect to the Nyquist frequency ${\omega }_{N_f/2+1}$: $S_{xx}({\omega }_{N_f/2+1-k})=S_{xx}({\omega }_{N_f/2+1+k}),k=1,\mathrm{...},N_f/2-1$. We introduce the spectral exponent $\alpha$ by the formula:

$$\log (S_{xx}(T_j))=\alpha \cdot \log (T_j))+d+{\eta }_j,j=1,\mathrm{...},N_f/2$$

(1)

where period $T_j=2\pi /{\omega }_j$, the value $\alpha$ is determined by the least squares method: ${\displaystyle {\sum }_{j=1}^{N_f/2}{\eta }_j^2}\to \underset{d,\alpha }{\min }$, ${\eta }_j$ is a sequence of random variables with zero mean. In each time window, the power spectrum was computed after transition to increments using an autoregressive model (maximum entropy estimate) [27] with an autoregressive order of $N/10$, i.e., 36 for a window of length 365 samples. The transition to increments before calculating the power spectrum aims to free from the dominant effect of low frequencies in the daily time series of GPS.

Wavelet-based entropy and wavelet-based spectral slope. The minimum normalized wavelet information entropy for the sample $x(t), t=1,,,N$ is determined by the formula:

$$En=-{\displaystyle \sum _{k=1}^N{p}_k\cdot \log (p_k)}/\log (N), p_k=c_k^2/{\displaystyle \sum _{j=1}^N{c}_j^2}$$

(2)

Here $c_k$ denoted the coefficients of the orthogonal wavelet decomposition. Daubechies bases were used with the number of vanishing moments from 1 to 10 [28]. We chose such a basis from this set for which the value (2) is minimal for a given time window and for the considered component of the displacement of the earth’s surface. By formula (2) $0\le En\le 1$.

For the orthogonal wavelet which is found from minimum of entropy, the wavelet power spectrum can be calculated as the average values of the squares of the wavelet coefficients at each level of detail of the decomposition:

$$S_k={\displaystyle \sum _{j=1}^{M_k}|c_j^{(k)}{|}^2}/M_k$$

(3)

In formula (3) $c_j^{(k)}$ are the wavelet coefficients at the detail level with number $k, k= 1,\mathrm{...}{k}_{\max }$, index $j= 1,\mathrm{...}{M}_k$, where $M_k$ is the number of wavelet coefficients at the detail level with number $k$, which corresponds to the frequency band $[\frac{1}{2^{(k+1)}\Delta s}, \frac{1}{2^k\Delta s}]$ with the central period $T_k=\frac{2\Delta s}{2^{-k}+2^{-(k+1)}}$ [28]. The maximum detail level number $k_{\max }$ is chosen such that it contains at least a given number $L_{\min }$ of wavelet coefficients. Recall that the number of wavelet coefficients at the level of detail decreases by a factor of 2 when the level number increases by one. The value $L_{\min }=4$ was used in the calculations, which gives the value $k_{\max }=6$ for the window length $N=365$. Values $S_k=S(T_k)$, are similar to conventional power spectra. The regression model is used to calculate the wavelet-based spectral slope:

$$\log (S(T_k))=\beta \cdot \log (T_k))+c+{\epsilon }_k,k=1,\mathrm{...},k_{\max }$$

(4)

where ${\epsilon }_k$ is a sequence of random variables with zero mean. The parameter $\beta$ in formula (4) is the wavelet-based spectral exponent, the value of which can be found by the least squares method: ${\displaystyle {\sum }_{k=1}^{k_{\max }}{\epsilon }_k^2}\to \underset{\beta ,c}{\min }$.

Donoho-Johnstone index (DJ-index). The threshold $T_{DJ}$ is defined by the formula [28,29]:

$$T_{DJ}=\sigma \sqrt{2\cdot \ln N}$$

(5)

The threshold (5) separates rather large (informative) in their absolute values wavelet coefficients from other coefficients which are considered to be noisy. Thus, we can consider the dimensionless signal characteristics $\gamma$, $0\le \gamma \le 1$, as the ratio of the number of the largest wavelet coefficients, for which the inequality $\left|c_k\right|> T_{DJ}$ is fulfilled, to the number N of all wavelet coefficients.

The value $\sigma$ in the formula (5) is the noise standard deviation estimate under the assumption that the noise is most concentrated in the 1st detail level of orthogonal wavelet decomposition. The estimate of standard deviation $\sigma$ should be robust with respect to outliers in the values of the coefficients at the first level. To provide this, a median estimate of the standard deviation for a normal random variable is used:

$$\sigma =med\left\{\left|c_k^{(1)}\right| , k=1,\mathrm{...},N/2\right\}/0.6745$$

(6)

For each station in each time window, 4 surface tremor statistics are computed for the 3 ground displacement components where с_k⁽¹⁾ are wavelet coefficients at the first level of detail. The estimate of the standard deviation σ from formula (6) determines the value (5) as a threshold for extracting noise wavelet coefficients. The quantity (5) is known in wavelet analysis as the Donoho–Johnstone threshold, and the expression for this quantity is based on the formula for the asymptotic probability of the maximum deviations of Gaussian white noise [28,29]. The DJ index $\gamma$ values can be interpreted as a measure of the non-stationarity of seismic noise. For stationary Gaussian white noise, the index $\gamma$ is zero. In [30,31,32] DJ-index was used as one of the main statistics for investigating properties of global seismic noise.

Summarizing, we will briefly describe in the list the characteristics of the values of the tremor statistics for stationary white Gaussian noise, which is considered as a standard of a completely random oscillation:

1. The “usual” spectral exponent $\alpha$ is the coefficient of regression between the logarithm of the power spectrum values and the logarithm of the period (“spectral slope”). For white noise it is equal to zero (white noise has a “flat spectrum”).

2. Wavelet-based spectral exponent $\beta$ - the regression coefficient between the logarithm of the mean values of the squared wavelet coefficients for each level of detail and the logarithm of the period of the central frequency corresponding to the level of detail of the wavelet decomposition (“spectral slope”). For white noise is zero.

3. Donoho-Johnstone threshold (DJ-index) $\gamma$ - the ratio of the number of “large” in absolute magnitude orthogonal wavelet coefficients to the total number of coefficients. For white noise is zero.

4. Minimum entropy $En$ of the distribution of squared orthogonal wavelet coefficients in the enumeration of wavelet basis functions in the class of Daubechies wavelets with the number of vanishing moments from 1 to 10. For white noise, the entropy is maximum.

Since we use 4 statistics for 3 components of the ground displacement vector, we get 12 statistics for each type of tremor. Let’s put them in one matrix:

$$\left[\begin{array}{cc}\left(\begin{array}{c}{\alpha }_{E,\max }\\ {\beta }_{E,\max }\\ {\gamma }_{E,\max }\\ E{n}_{E,\min }\end{array}\begin{array}{c}{\alpha }_{N,\max }\\ {\beta }_{N,\max }\\ {\gamma }_{N,\max }\\ E{n}_{N,\min }\end{array}\begin{array}{c}{\alpha }_{U,\max }\\ {\beta }_{U,\max }\\ {\gamma }_{U,\max }\\ E{n}_{U,\min }\end{array}\right)& \left(\begin{array}{c}{\alpha }_{E,\min }\\ {\beta }_{E,\min }\\ {\gamma }_{E,\min }\\ E{n}_{E,\max }\end{array}\begin{array}{c}{\alpha }_{N,\min }\\ {\beta }_{N,\min }\\ {\gamma }_{N,\min }\\ E{n}_{N,\max }\end{array}\begin{array}{c}{\alpha }_{U,\min }\\ {\beta }_{U,\min }\\ {\gamma }_{U,\min }\\ E{n}_{U,\max }\end{array}\right)\end{array}\right]$$

(7)

On the left side of the matrix (7) there are the designations of the characteristics of active tremor, on the right side of the matrix (7) there are the designations of the characteristics of passive tremor. Recall that we are considering a regular grid of nodes with a size of 100 nodes in longitude and 150 nodes in latitude. For each time window with a length of 365 days with a shift of 7 days, the values of all 12 sequences of the distribution of tremor properties at the nodes of the regular grid are determined.

4. Extreme value probability density maps

Let $U$ there be any value from the matrix (7), the estimates of which were obtained in a moving time window. For each grid node $(i,j)$ and for each time window labeled $t$ at the right end of the window, we find the 10 closest seismic working stations, which gives 10 values of $U$. Let’s take their median value $U_{ij}^{(t)}$. The values $U_{ij}^{(t)}$ form an “elementary” map corresponding to a time window of 365 days. Consider the values of the parameter as a function of two-dimensional vectors $z_{ij}=(x_i,y_j)$ of longitudes and latitudes of nodes $(i,j)$ in an explicit form: $U_{ij}^{(t)}\equiv U^{(t)}(z_{ij})$. For each “elementary map” with a discrete time index $t$, we find the coordinates of the nodes $z_{ mn}^{(t)}=(x_m^{(t)},y_n^{(t)})$ at which the extreme value $U$ is reached with respect to all other nodes of the regular grid.

As noted above, the choice of the minimum for entropies and the maximum for spectral slopes is due to the considerations that the regions of active tremor should be characterized by the values of statistics that reflect the maximum deviation from the properties of white noise. For passive tremor, the opposite is true: maxima for entropies and minima for spectral slopes. A cloud of 2D vectors $z_{ mn}^{(t)}$ considered within a certain time interval $t\in [t_0,t_1]$ composes some random set. Let us estimate two-dimensional probability density function for each node of the regular grid. For this, the Parzen–Rosenblatt estimate with the Gaussian kernel function [26] will be applied:

$$p(z_{ij}|t_0,t_1)=\frac{1}{2\pi h^2(t_1-t_0+1)}{\displaystyle \sum _{t=t_0}^{t_1}{\displaystyle \sum _{mn}\exp \left(-\frac{|z_{ij}-z_{ mn}^{(t)}{|}^2}{2h^2}\right)}}$$

(8)

Here, $h$ is the averaging radius, integer indices $t_0,t_1$ enumerate the maps of each time window. The $(t_1-t_0+1)$ equals the number of time windows in the considered time interval. The width of the smoothing band $h={0.25}^{\circ }$ was used, which is approximately equal to 28 km.

In our analysis, time windows of 2 lengths are used: a “short” window of 365 days in length, which is used to obtain estimates of tremor parameters, and a “long” window of 1095 days (3 years) in length, in which a sequence of estimates of the distribution densities of extreme values of statistics is obtained tremor. Estimates of characteristics from matrix (7) are calculated using “short” time windows 365 days long. Next, the probability densities of extreme values are calculated from the values that fall within the current “long” time window of 1095 days, which are also shifted by 7 days. Each “long” window includes 105 “short” time windows.

As a result of such assessments, 24 maps of the distribution density of extreme values of tremor parameters were constructed, presented in Figure 3, on which 12 maps for active tremor are in the upper half, and 12 maps for passive tremor are in the lower half.

Purely visually, a correlation of probability density maps of the distribution of extreme values of tremor statistics is noticeable. The principal component method [33] is using to extract the most common spatial characteristics. In each “long” 3-year time window, a correlation matrix of the size 12×12 between the probability densities of the distribution of extreme values of tremor statistics is calculated. The eigenvector corresponding to the maximum eigenvalue of the correlation matrix is found and the squares of its components were used as weights to calculate the average weighted probability density in each “long” window. The coordinates of the point that realizes the maxima of average density were found. These operations performed makes it possible to obtain the trajectory of geographic coordinates of points that realize the maximum values of the average probability density. These points are called singular tremor points.

The results of weighting the probability density maps of the distribution of extreme properties of tremor by the method of principal components are shown in Figure 4. On the map 4(a), red circles indicate the positions of singular points of active tremor, which are combined into 5 clusters. The clusters are numbered in descending order of the latitude of their centers of mass; the cluster numbers are also shown. On the map 4(b), blue circles indicate the positions of singular points of passive tremor, which form 3 clusters. On Figure 4(c) and (d) there are plots of histograms of the distribution of the number of singular points in each cluster for each type of tremor.

Figure 5 shows graphs of changes in the squares of the components of sequences of correlation matrices for probability density maps of the distribution of extreme values, respectively, for active (12 upper graphs) and passive (12 lower graphs). As noted above, these values are used as weight functions to obtain weighted average probability density maps. The result of such a weighing operation using the principal component method is shown in Figure 4. The graphs in Figure 5 show the average values of the weights calculated for all “long” time windows. For each tremor variant, 4 maximum average values of weights are highlighted in red - they select those components from matrix (7) that make the greatest contribution to obtaining the weighted average sum of probability densities.

The trajectories of changes in geographic latitudes and longitudes of singular points over time are shown in Figure 6(a) and (b). Along the trajectories of change of coordinates there are marks of belonging to clusters of stations. Figure 6(c) shows how the distances between the singular points of active and passive tremor change as the time window of the length 3 years moves from left to right. From this graph it can be seen that the concentration of singular points of both types of tremor in those regions where their positions are very close to each other, namely, cluster #2 of active singular points and cluster #1 of passive singular points mainly occurs in different time windows. The exception is the time intervals with labels of the right ends of the “long” 3-year windows 2012-2012.5 and 2015.5-2016, when the distances from singular points of different types are less than 100 km, but do not reach zero values.

Let ${\tilde{\lambda }}_1$ be the normalized maximum eigenvalue of the correlation matrix:

$${\tilde{\lambda }}_1={\lambda }_1/{\displaystyle \sum _{k=1}^q{\lambda }_k}$$

(9)

In formula (9) ${\lambda }_1,\mathrm{...},{\lambda }_q$ are the eigenvalues of the correlation matrix, sorted in descending order: ${\lambda }_k\ge {\lambda }_{k+1}, {\displaystyle \sum _{k=1}^q{\lambda }_k}=q$. According to the principal component method [33], the value ${\tilde{\lambda }}_1$ is equal to the share of the total variance attributable to the first principal component and can be interpreted as a measure of synchronization of changes in the scalar components of the multivariate time series of the probability density values of the tremor statistics from the matrix in formula (7) in nodes regular grid covering the region under study.

Figure 7 shows graphs of changes in the eigenvalues of the correlation matrix for active and passive tremor. It can be seen from these plots that the active tremor (Figure 7(a)) is much more synchronized than the passive tremor (Figure 7(b)). Moreover, for active tremor, there is a very strong synchronization for the positions of the right end of the 3-year time window from 2013 to 2015. Given the length of the window, this synchronization takes a period of time from 2010 to 2015.

As was shown in [23,24], after a couple of the strongest earthquakes in the world, Maule in Chile on 2010.02.27 and Tohoku in Japan on 2011.11.03, there was an explosive increase in the radius of spatial correlations of global seismic noise properties. It is possible that the increase in synchronization of active tremor in the Western United States in 2010-2015 is also a response to this pair of mega earthquakes. Note that a sharp jump in the magnitude of global correlations of the earth’s surface tremor in 2010-2011 was shown in [19], and [34] showed that before the 2011 Tohoku earthquake in Japan, long-term correlations arose between seismic noise in California and Japan.