Entropy of GPS-measured Earth Tremor

Abstract

Modern GPS networks make it possible to study the tremors of the earth’s surface from the point of view of identifying anomalous areas. The use of the entropy of the distribution of wavelet coefficients provides a tool for highlighting the hidden and non-obvious properties of the earth’s surface tremors. The principal component method makes it possible to identify the most important general trends in the behavior of informative tremor statistics and determine areas of anomalous behavior. The application of these methods to the analysis of GPS data in California is presented. Particular attention is paid to time intervals and areas (clusters) with extreme entropy values. Periodicities in the occurrence of strong jumps in the average entropy of the entire region have been discovered, of which the period of 95 days is dominant. The trend of migration of areas of maximum entropy from the South to the North has been identified. As a result of the analysis, it was found that the area of minimum entropy values gravitates toward the San Andreas fault, and the vicinity of San Francisco has the selected properties of maximum information content and attracts low entropy trajectories.

1. Introduction

The fine structure of the time series of displacements of the earth’s surface measured using 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]. In works [5, 6, 7, 8], the maximum likelihood method was used to estimate the parameters of the shapes of power spectra and noise amplitude in various regions of the world. The noise structure and error estimates are considered in [7, 9, 10]. Errors in estimating the rate of displacement depending on the shape of the spectrum and the intensity of the noise component of the time series were analyzed in [11, 12]. Parametric models of GPS time series for the analysis of tectonically active areas, including the analysis of phase correlations, were used in [13, 14].

The structure of the high-frequency noise and low-frequency seasonal components of GPS time series in connection with the task of estimating the rates of displacement of tectonic plates was studied in [15, 16, 17, 18]. Seasonal effects in GPS signal series associated with hydrological load and tectonic movements of the earth’s crust were studied in [19, 20, 21, 22, 23, 24, 25, 26] using various methods, including a multivariate statistical approach. The influence of ionospheric and tropospheric delays on the sensitivity of GPS solutions was considered in [27]. Hypotheses about the occurrence of “anomalous” harmonics in GPS time series are considered in [28].

Multivariate statistics methods (principal components, empirical orthogonal functions, and singular spectrum estimation) in [29, 30] were used to identify common spatial and temporal components of GPS time series. A joint analysis of the noise component of time series and accelerometer readings was carried out in [31]. Non-stationary effects were studied in [32, 33, 34] to assess the mutual movements of blocks of the earth’s crust and the stability of station positioning.

The coherence of ground tremor was analyzed in [35, 36]. In [37, 38], coherence analysis of GPS time series was used to assess seismic hazards in Japan and California.

This chapter presents the further development of methods for analyzing ground tremors proposed in [39, 40] and analyzes the structure of high-frequency noise in time series of a network of GPS stations in the Western United States during the time period 2009–2022. The Western USA (mainly California) is one of the most tectonically active regions. A dense network of GPS stations makes it possible to study in detail the characteristics of vibrations of the earth’s surface and try to compare them with known regions of seismic activity. Previously, in [41], a study of the properties of seismic noise in Southern California was carried out in comparison with the seismic regime.

In this work, to analyze daily GPS time series, we use estimates of the entropy of the distribution of the squared wavelet coefficients of the expansion into orthogonal wavelets in a sliding time window. Since, due to the orthogonality of the wavelet transform, the squared modulus of each wavelet coefficient is equal to a certain “atom” of oscillation energy, the entropy of the energy distribution over the complete set of independent time-frequency intervals is actually studied.

The term tremor has many different meanings in geophysics. In seismology, tremor refers to slow movements of the earth’s crust as a result of the so-called quiet earthquakes. In volcanology, tremor is a shaking of the earth’s surface in the vicinity of volcanoes, usually associated with eruptions. In this chapter, tremor is the increments of daily GPS time series. The transition to increments is a well-known operation in time series analysis, which allows one to reduce the influence of strong low-frequency components and increase the stationarity of the analyzed signals.

2. Data used

Data on earth surface displacements measured by GPS were taken from the Nevada Geodetic Laboratory website at: http://geodesy.unr.edu/gps_timeseries/tenv3/IGS14/ [42].

Positions of 699 GPS stations for a rectangular domain ${32.5}^{\circ }\le \mathrm{Lat}\le {51}^{\circ }$, $-{128.5}^{\circ }\le \mathrm{Lon}\le -{115}^{\circ }$ which have daily time series of 5113 samples ranging from the beginning of 2009 to the end of 2022 (14 years), having a total number of gaps of less than 365 samples, and the longest gap of less than 182 samples, are presented at the left panel of Figure 1. Vertical components of earth displacement are investigated. Gaps in the GPS time series are filled using information from the left-hand and right-hand vicinities of the gap of the same length as the length of the gap [35].

Figure 1.

(a) It presents the positions of 699 GPS stations and their partition into 18 clusters. The numbered circles represent the centers of gravity of the clusters, the blue lines are the boundaries between Voronoi cells. (b) It shows a plot of pseudo-F statistics, which made it possible to select 18 as the number of clusters.

The set of stations was previously divided into clusters. To divide the network of stations into clusters, 18 reference points were selected (Figure 1). The number 18 was chosen as the optimal number of station clusters, which split their “cloud” using the k-means method. Let us partition the set of station position vectors $\overrightarrow{\zeta }$ into a given trial number $q$ of clusters using the popular k-means clustering method [43]. Denote by $C_r,r=1,\dots ,q$ clusters, let ${\overrightarrow{z}}_r={\sum }_{\overrightarrow{\zeta }\in C_r}\overrightarrow{\zeta }/n_r$ be the center of mass vector of the cluster $C_r$, $n_r$ be the number of vectors in the cluster $C_r$, ${\sum }_{r=1}^q{n}_r=N$. A vector $\overrightarrow{\zeta }\in C_r$ if the distance $\mid \overrightarrow{\zeta }-{\overrightarrow{z}}_r\mid$ is minimal among the positions of all cluster centers. The k-means method minimizes the sum of distances:

$$G\left({\overrightarrow{z}}_1,\dots ,{\overrightarrow{z}}_q\right)=\sum _{r=1}^q\sum _{\overrightarrow{\zeta }\in C_r}{\left|\overrightarrow{\zeta }-{\overrightarrow{z}}_r\right|}^2\to \underset{{\overrightarrow{z}}_1,\dots ,{\overrightarrow{z}}_q}{min}$$E1

with respect to the position of cluster centers ${\overrightarrow{z}}_r$. Let $\mathrm{\Phi }\left(q\right)=\underset{{\overrightarrow{z}}_1,\dots ,{\overrightarrow{z}}_q}{min}G\left({\overrightarrow{z}}_1,\dots ,{\overrightarrow{z}}_q\right)$. We used a trial number of clusters in the range $2\le q\le 50$. The problem of choosing the best number of clusters $q^{\ast }$ was solved using the pseudo-F-statistic maximum criterion [44]

$$\mathit{PFS}\left(q\right)={\sigma }_1^2\left(q\right)/{\sigma }_0^2\left(q\right)\to \underset{2\le q\le 50}{max}$$E2

where

$${\sigma }_0^2\left(q\right)=\frac{\mathrm{\Phi }\left(q\right)}{N-q},{\sigma }_1^2\left(q\right)=\sum _{r=1}^q\frac{n_r}{N}\cdot {\left|{\overrightarrow{z}}_r-{\overrightarrow{z}}_0\right|}^2,{\overrightarrow{z}}_0=\frac{1}{N}\sum _1^N\overrightarrow{\zeta }$$E3

The right-hand graph in Figure 1 presents the values of the pseudo-F-statistic as a function of the trial number of clusters. The number 18 on the graph of the pseudo-F-statistics is the breakpoint of the dependence on the trial number of clusters and realizes the maximum for the number of clusters from 2 to 50. On the graph of the pseudo-F statistics, there are two local maxima with close values for the number of clusters 18 and 21. Of these two values, the number 18 was chosen as the smaller one. In Figure 1, the division of a set of stations into 18 clusters is shown along with Voronoi cells, which clearly indicate that the stations belong to one or another cluster.

Clusters of stations are ordered by increasing latitude of the positions of their centers of gravity. Table 1 shows for each cluster (first line) the number of stations in the cluster (second line).

Clust#	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18
Nsta	65	97	46	66	50	52	57	32	14	50	15	20	11	17	39	19	35	14

Table 1.

Number of N_sta stations in each Clust# cluster.

3. Principal components of increments in a sliding time window

Since the goal is to study the tremor of the earth’s surface, that is, the high-frequency part of the earth’s surface displacements, the analysis was carried out for increments of time series. The transition to increments reduces the dominant influence of low frequencies in the daily GPS time series and ensures stationarity of time series fragments within the 365-day time windows that are used further.

The division of stations into 18 clusters is used for subsequent application of the principal component method [45]. For each cluster of stations, the first principal component of the time series of increments of vertical displacements of the earth’s surface was calculated in a sliding adaptation time window of 365 days in length.

Let us consider several time series $P\left(t\right)={\left(P_1\left(t\right),\dots ,P_m\left(t\right)\right)}^T$ of common dimensionality $m$. In our case $m=18$. It is necessary to estimate the first principal component in a sliding time window of long $L$ samples. To do this, consider samples with time indexes $t$ under the condition $s-L+1\le t\le s$, $s$ - the right end of the time window. The correlation matrix $\mathrm{\Phi }\left(s\right)$ of the size $m\times m$ is calculated using the formulas:

$$\mathrm{\Phi }\left(s\right)=\left({\phi }_{\mathit{ab}}^{\left(s\right)}\right),{\phi }_{\mathit{ab}}^{\left(s\right)}=\sum _{t=s-L+1}^s{q}_a^{\left(s\right)}\left(t\right)q_b^{\left(s\right)}\left(t\right)/L,a,b=1,\dots ,m$$E4

where

$$\begin{array}{ll}{q}_a^{\left(s\right)}\left(t\right)& =\left(P_a\left(t\right)-{\overline{P}}_a^{\left(s\right)}\right)/{\sigma }_a^{\left(s\right)},{\overline{P}}_a^{\left(s\right)}=\sum _{t=s-L+1}^s{P}_a\left(t\right)/L,{\left({\sigma }_a^{\left(s\right)}\right)}^2=\\ & =\sum _{t=s-L+1}^s{\left(P_a\left(t\right)-{\overline{P}}_a^{\left(s\right)}\right)}^2/\left(L-1\right),a=1,\dots ,m\end{array}$$E5

Let ${\theta }^{\left(s\right)}={\left({\theta }_1^{\left(s\right)},\dots ,{\theta }_m^{\left(s\right)}\right)}^T$ be the eigenvector of the matrix $\mathrm{\Phi }\left(s\right)$ with the maximum eigenvalue. Let us calculate the value:

$${\psi }^{\left(s\right)}\left(t\right)=\sum _{\alpha =1}^m{\theta }_a^{\left(s\right)}\cdot q_a^{\left(s\right)}\left(t\right)$$E6

and define the scalar time series $\psi \left(t\right)$ of the adaptive first principal component in a sliding time window in accordance with the formula:

$$\psi \left(t\right)=\left\{\begin{array}{c}{\psi }^{\left(L-1\right)}\left(t\right),0\le t\le \left(L-1\right)\\ {\psi }^{\left(t\right)}\left(t\right),t\ge L\end{array}\right.$$E7

Eq. (4)–(7) are applied independently in each time window. According to them, in the first time window, a series $\psi \left(t\right)$ consisting of $L$ values, the values of the principal component are calculated according to Eqs. (6)–(7). In all subsequent time windows, the time series values $\psi \left(t\right)$ correspond to the single estimates at the far right ends of the respective windows. Thus, outside the first time window, the values of the series $\psi \left(t\right)$ depend only on past values $P\left(t\right)$. The calculations of the principal components were carried out in a sliding time window of 365 days with a step of 3 days. The length of the adaptation time window of 365 days is equal to 1 year and was chosen for the reason that it coincides with the natural period of variations in GPS time series. After calculating the first principal components, information is compressed and the dimension of the analyzed data is reduced from 699 to 18. Figure 2 shows the graphs of the first principal components for 4 clusters as an example.

Figure 2.

Graphs of first principal components within moving time window of the length 365 days for 4 clusters.

4. Entropy of the first principal components

The minimum normalized wavelet-based entropy for the sample $x\left(t\right),t=1,\mathrm{..},N$ is determined by the formula:

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

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 [46]. We chose such a basis from this set for which the value Eq. (8) is minimal for a given time window and for the considered component of the displacement of the earth’s surface. By Eq. (8) $0\le \mathit{En}\le 1$. We will consider estimates of the entropy of principal components of GPS time series increments in a sliding time window 365 days long with a shift of 3 days.

Entropy Eq. (8) has been repeatedly used in [41, 47] to analyze the variability of seismic noise in different regions of the Earth and at the global level as an independent tool. By its design, entropy Eq. (8) is also multiscale, like the entropy proposed in [48] for studying the properties of random signals. In [49], the author used the non-extensive Tsallis entropy to analyze seismic noise. The “natural time” approach to the analysis of random data uses a related definition of entropy in [50].

Figure 3 shows graphs of entropy values for the first principal components of the increments of the vertical components of the displacement of the earth’s surface, estimated in a sliding time window of 365 days in length with a displacement of 3 days for 4 clusters of stations, the same as those presented in the previous Figure 2.

Figure 3.

Graphs of minimum normalized wavelet-based entropy for first principal components within moving time window of the length 365 days with mutual shift 3 days for 4 clusters.

5. Mean entropy and its wavelet transform maximum modulus stepwise approximation

Entropy is a measure of the complexity of the behavior of a random signal. A simple signal which is close by its properties to white noise has maximum entropy. On the contrary, the entropy of a signal, saturated with non-stationary components and sharp changes in its amplitude, reaches minimum values. Therefore, of particular interest are those time fragments and those areas of the region under study where entropy reaches extreme values, minimums or maximums. Next, to isolate such time fragments, a smoothing method is used, the purpose of which is to approximate the signal with a piecewise step function. The scale-dependent piecewise step approximation makes it possible to formally determine the moments in time of a significant change in the current average signal values. These moments in time can further be considered as a kind of event in the dynamics of a random time series. To construct such an approximation, continuous wavelet transforms with a kernel in the form of derivatives of a Gaussian function are used, which makes it possible to formalize the selection of time points at which a scale-dependent significant change in the average value of the noisy signal occurs [46, 51].

Let $Y\left(t\right)$ be an arbitrary signal and consider the corresponding smoothed signal obtained by using a scale-dependent kernel:

$$c_0\left(t,a\right)=\underset{-\infty }{\overset{+\infty }{\int }}Y\left(t+\mathit{av}\right)\cdot {\psi }_0\left(v\right)\mathit{dv}/\underset{-\infty }{\overset{+\infty }{\int }}{\psi }_0\left(v\right)\mathit{dv}$$E9

where $a>0$ is a time scale and ${\psi }_0\left(t\right)$ is some function decaying sufficiently fast on both sides of its single maximum; further on, ${\psi }_0\left(t\right)=exp\left(-t^2\right)$ shall be used. The wavelet function is defined as: ${\psi }_n\left(t\right)={\left(-1\right)}^n\cdot \frac{d^n{\psi }_0\left(t\right)}{{\mathit{dt}}^n}\equiv {\left(-1\right)}^n\cdot {\psi }_0^{\left(n\right)}\left(t\right)$. Using integration by part and exploiting the fast decay of the function ${\psi }_0\left(t\right)$ at $t\to \pm \infty$, the following formula for the Taylor’s coefficients (the $n$-th derivative of the smoothed signal, divided by $n!$) of the smoothed signal is obtained:

$$c_n\left(t,a\right)\equiv \frac{1}{n!}\frac{d^n{c}_0\left(t,a\right)}{{\mathit{dt}}^n}=\underset{-\infty }{\overset{+\infty }{\int }}Y\left(t+\mathit{av}\right){\psi }_n\left(v\right)\mathit{dv}/a^n\underset{-\infty }{\overset{+\infty }{\int }}{v}^n{\psi }_n\left(v\right)\mathit{dv}$$E10

The WTMM-point $\left(t,a\right)$ for $n\ge 1$ is defined as the point for which $\mid c_n\left(t,a\right)\mid$ exhibits a local maximum with respect to time $t$ for a given time scale $a$. For $n=0$, the WTMM-points are defined as points of local extremes (maxima or minima) of the smoothed signal $c_0\left(t,a\right)$, that may be joined to form chains. The set of all chains forms a WTMM-skeleton of the signal. If ${\psi }_0\left(t\right)$ is a Gaussian, then a given WTMM-skeleton chain does not end when the scale $a$ is decreased [52]. The WTMM-points for the 1st-order derivative $c_1\left(t,a\right)$ indicate time points of the maximum trend (positive or negative) of the smoothed signal $c_0\left(t,a\right)$ for the given scale. This can be used to extract points of large and abrupt changes of the mean value. These times correspond to the occurrence of rather long WTMM-chains that begin to grow. Let us define a WTMM stepwise approximation $S_Y\left(t|a_{\ast }\right)$ for the signal $Y\left(t\right)$, as a function which is equal to the sequence of constant values $g_k$ in the successive intervals $t\in \left[{\tau }_k\left(a_{\ast }\right),{\tau }_{k+1}\left(a_{\ast }\right)\right]$. Here ${\tau }_k\left(a_{\ast }\right)$ are the beginnings of the WTMM-chains for $c_1\left(t,a\right)$ which exceed the threshold time scale $a_{\ast }$ and $g_k={\sum }_{t={\tau }_k\left(a_{\ast }\right)}^{{\tau }_{k+1}\left(a_{\ast }\right)}Y\left(t\right)/\left({\tau }_{k+1}\left(a_{\ast }\right)-{\tau }_k\left(a_{\ast }\right)+1\right)$ is equal to the mean value of $Y\left(t\right)$ within the time interval $\left[{\tau }_k\left(a_{\ast }\right),{\tau }_{k+1}\left(a_{\ast }\right)\right]$.

Let us calculate the weighted average entropy of all clusters. As weights we use the squares of the eigenvector components of the correlation matrix of the entropy values of the first principal components of all clusters. The calculation of such an average is the following application of the principal component method. In Figure 4b the black line shows the graph of the weighted average entropy, and in Figure 4c the values of the weights. From graph 4(c), it can be seen that the weight of cluster #7 stands out significantly compared to the weights of other clusters. Thus, cluster #7 could be regarded as the most informative. The WTMM piecewise stepwise approximation for time variations of the weighted average entropy of the principal components from all clusters is represented in Figure 4b by the red line. To construct this approximation, a sampling scale $a_{\ast }=10$ was used. Since the entropy values were calculated in a sliding time window of length 365 days with an offset of 3 days, the dimensional value of this scale is equal to 30 days. The absolute values of the increments of the average entropy for the piecewise-step approximation in Figure 4a form a sequence of 56 “events” of sufficiently large changes in the average entropy.

Figure 4.

(a) – absolute increments of average weighted entropy (b, black line) for piecewise stepwise approximation (b, red line). The blue line in plot (b) represents the linear trend of the average entropy. Figure (c) represents the weight values of each cluster as the squares of the eigenvector components corresponding to the maximum eigenvalue of the correlation matrix.

The graphs in Figure 4 provide an integral assessment of the change in entropy of the entire region and allow us to identify a close to periodic regime (with a period of about 3 years) of changes in the statistical structure of the earth’s surface tremor. In addition, a general trend of decreasing entropy over the entire observation period is visible. A decrease in average entropy means that the composition of GPS time series is enriched with non-stationary components, which reflect an increase in the intensity of mutual movements of small blocks of the earth’s crust.

Let us try to find out whether there is a statistically significant periodicity in the occurrence of jumps in the average entropy value, presented in Figure 4a. To test this assumption, let us estimate the periodic components of the shock intensity using the parametric model proposed in [53].

Let $t_i,i=1,\dots ,N$ be the times of the sequence of events observed on the interval $\left(0,T\right]$. Consider the following intensity model containing a periodic component:

$$\lambda \left(t\right)=\mu \cdot \left(1+acos\left(\mathit{\omega t}+\phi \right)\right)$$E11

where frequency $\omega$, amplitude $a,0\le a\le 1$, phase angle $\phi$, $\phi \in \left[0,2\pi \right]$ and multiplier $\mu >0$ (describing the Poisson part of the intensity) are the model parameters. Thus, the Poisson part of the intensity is modulated by a harmonic oscillation. Let us fix some frequency value $\omega$. The log likelihood function [54] in this case for a series of observed events is equal to:

$$\begin{array}{ll}lnL\left(\mu ,a,\phi |\omega \right)& =\sum _{t_i}ln\left(\lambda \left(t_i\right)\right)-\underset{0}{\overset{T}{\int }}\lambda \left(u\right)\mathit{du}=\\ & =Nln\left(\mu \right)+\sum _{t_i}ln\left(1+acos\left(\omega t_i+\phi \right)\right)-\mathit{\mu T}-\frac{\mathit{\mu a}}{\omega }\left[sin\left(\mathit{\omega T}+\phi \right)-sin\left(\phi \right)\right]\end{array}$$E12

Taking the maximum of expression Eq. (12) with respect to the parameter $\mu$, it is easy to find that:

$$\mu =\widehat{\mu }\left(a,\phi |\omega \right)=\frac{N}{T+a\left(sin\left(\mathit{\omega T}+\phi \right)-sin\left(\phi \right)\right)/\omega }$$E13

Substituting Eq. (13) into formula Eq. (12) we get

$$ln\left(L\left(\widehat{\mu },a,\phi |\omega \right)\right)=\sum _{t_i}ln\left(1+acos\left(\omega t_i+\phi \right)\right)+N\cdot ln\left(\widehat{\mu }\left(a,\phi |\omega \right)\right)-N$$E14

It should be noted that the expression $\widehat{\mu }\left(a=0,\phi |\omega \right)\equiv {\widehat{\mu }}_0=N/T$ is an estimate of the intensity of the process, provided that it is homogeneous Poisson (purely random). Thus, the increment of the log-likelihood function due to consideration of a richer intensity model with a harmonic component with a given frequency $\omega$ than for a purely random flow of events is equal to:

$$\mathrm{\Delta }lnL\left(a,\phi |\omega \right)=\sum _{t_i}ln\left(1+acos\left(\omega t_i+\phi \right)\right)+N\cdot ln\left(\widehat{\mu }\left(a,\phi |\omega \right)/{\widehat{\mu }}_0\right)$$E15

Let

$$R\left(\omega \right)=\underset{a,\phi }{max}\mathrm{\Delta }lnL\left(a,\phi |\omega \right),0\le a\le 1,\phi \in \left[0,2\pi \right]$$E16

Function Eq. (16) can be considered as a generalization of the spectrum for a sequence of events. The graph of this function shows how “more profitable” the periodic intensity model is compared to a purely random model. The maximum values of function Eq. (16) highlight the frequencies present in the stream of events.

An important issue in applying this method to real data is determining the statistical significance of the obtained peak values of statistics $R\left(\omega \right)$. To solve it, one can apply Wilks’ asymptotic theory [55]. Let us consider 2 hypotheses for the same data set consisting of independent observations:

1. hypothesis $H_0$: $X^{\left(N\right)}$distributed according to density $p_0\left(X^{\left(N\right)}|{\theta }_0\right)$;

2. hypothesis $H_1$: $X^{\left(N\right)}$ distributed according to density $p_1\left(X^{\left(N\right)}|{\theta }_1\right)$.

where ${\theta }_0$ and ${\theta }_1$ are vectors of unknown parameters having dimensions $m_0$ and $m_1$, and the hypothesis $H_1$ is more “rich”: $m_1>m_0$, and the vector of parameters ${\theta }_1$ completely includes the components of the vector ${\theta }_0$. Let us consider the difference between the log likelihoods for these two hypotheses, provided that the parameter vectors are taken from their maximum likelihood estimates:

$$\mathrm{\Delta }lnL\left(X^{\left(N\right)}\right)=ln\left(\underset{{\theta }_1}{max}{p}_1\left(X^{\left(N\right)}|{\theta }_1\right)\right)-ln\left(\underset{{\theta }_0}{max}{p}_0\left(X^{\left(N\right)}|{\theta }_0\right)\right)$$E17

It is obvious that $\mathrm{\Delta }lnL\left(X^{\left(N\right)}\right)\ge 0$. According to Wilks’ theorem, if the hypothesis $H_0$ is valid, quantity Eq. (17) has an asymptotic distribution:

$$\mathrm{\Delta }lnL\left(X^{\left(N\right)}\right)\sim \frac{{\chi }_m^2}{2},m=m_1-m_0,N\to \infty$$E18

In our case $m=2$, and therefore, the doubled value Eq. (18) has an asymptotic distribution ${\chi }_2^2$ density equal to $e^{-x/2}/2$, and the value Eq. (18) itself is distributed asymptotically with an exponential density $e^{-x}$ or

$$Pr\left\{R\left(\omega \right)<x\right\}=1-e^{-x},N\to \infty$$E19

provided that the analyzed sequence of time points is distributed according to the Poisson law with constant intensity. Expression Eq. (19) allows us to give thresholds for statistics $R\left(\omega \right)$, which allow us to state that only when exceeded, the sequence of time instants differs from the Poisson sequence with a given probability. In particular, setting the threshold probability equal to 0.99, we obtain that the threshold is equal to ln(100) = 4.6.

Figure 5 shows graphs of function estimates $R\left(\omega \right)$ for 1000 frequency values with periods from 10 to 1000 days, the values of which are taken uniformly on a logarithmic scale.

Figure 5.

Graph of increments $R\left(\omega \right)$ of the maximum value of the log-likelihood function of a periodic model of the intensity of mean entropy jumps compared to a purely random Poisson process model.

If we choose the value 4.6 as the statistical $R\left(\omega \right)$ threshold, above which model (11) has an advantage over the Poisson model with a probability exceeding 0.99, then 5 selected period values arise: 11.4, 18.9, 27.8, 28.7 and 95 days, of which the last period is clearly dominant.

6. Dynamics of extreme entropy values

Figure 6(a1,b1) shows graphs of changes in cluster numbers in which the minimum (a1) and maximum (b1) entropy values were realized for each position of a time window of 365 days in length. These graphs show that there are quite long periods of time when extreme entropy values are delayed in close clusters of stations. Having constructed histograms of the distribution of cluster numbers in which extreme entropy values are realized (Figure 6(a2,b2)), one can notice that there are clusters that “attract” these extreme values to themselves. In particular, entropy minima are most often realized in clusters #4 and #7, and maximum values are in clusters #15 and #13. Such differences highlight “characteristic areas” of earth tremor [39, 40], which reflect differences in the dynamics of the behavior of crustal blocks.

Figure 6.

Graphs (a1) and (a2) represent sequences of cluster numbers in which minima (graph (a1)) and maxima (graph (a2)) of entropy values were realized in each time window of 365 days with an offset of 3 days. Graphs (b1) and (b2) represent histograms of the distribution of cluster numbers, and which had minima (b1) and maxima (b2) of entropy.

To represent the spatiotemporal dynamics of the behavior of extreme entropy values, we calculate a histogram of the distribution of cluster numbers in which minima and maxima are realized in a sliding time window, which includes a certain number of neighboring extreme entropy values, presented in the graphs of Figure 6(a1,b1). Since each entropy value is calculated in a “small” sliding time window of length 365 days with an offset of 3 days, then if we take $M$ neighboring entropy values, then the dimensional length of the “large” time window in which the histograms are calculated will be equal to $365+3\cdot \left(M-1\right)$ days. With a value of $M=244$, the length of the “large” time window of the histogram will be equal to 1094 days, which is approximately equal to 3 years. Histograms of extreme entropy values in a 3-year time window are presented in Figure 7 in the form of two-dimensional diagrams. A similar method for presenting changes in the spatiotemporal properties of extreme values of global seismic noise statistics was previously used in [56].

Figure 7.

Histograms of clusters numbers where minima (a) and maxima (b) of first principal component entropy are realized within moving time window of the length 3 years.

Note that in Figure 7b, there is a noticeable trend in the space of cluster numbers for the evolution of the maximum values of the probability of the distribution of entropy maximums for the timestamps of the right end of the “large” 3-year time window from 2012 to 2019 from cluster #4 to cluster #17. Considering that the clusters are numbered in ascending latitude, physically this change corresponds to the migration of entropy maxima from south to north for the time interval 2009–2019. (taking into account that the time stamps in Figure 7 correspond to the right end of the 3-year time window).

Let $p_t^{\left(j\right)}$ be the probabilities of histograms of the distribution of extreme entropy values (minimums or maximums) for the moment in time of the right end $t$ of a “large” time window of 3 years in length, which are presented in Figure 7, ${\sum }_j{p}_t^{\left(j\right)}=1$. Let us denote

$$j^{\ast }\left(t\right)=\underset{1\le j\le 18}{argmax}{p}_t^{\left(j\right)}$$E20

the number of the cluster in which the extreme entropy value was realized in the 3-year time window with the time stamp $t$ of the right end of the window. Let $\left(X_j,Y_j\right)$ be the longitude and latitude of the centers of gravity of station clusters, $j=1,\dots ,18$. Then the jump-like trajectory of changes in the centers of gravity of those cluster numbers for which the maximum probability of extreme entropy values is realized will be a sequence of points $\left(X_{j^{\ast }\left(t\right)},Y_{j^{\ast }\left(t\right)}\right)$. In addition, you can calculate the continuous weighted average trajectory of the distribution of extreme entropy values using the formulas:

$$X_t=\sum _j{p}_t^{\left(j\right)}\cdot X_j,Y_t=\sum _j{p}_t^{\left(j\right)}\cdot Y_j$$E21

In Figure 8, the dynamics of changes in the probability distribution of extreme entropy values is transferred from the space of cluster numbers, as in Figure 7, to the plane of geographic coordinates. In Figure 8a, the blue line shows the jump trajectory $\left(X_{j^{\ast }\left(t\right)},Y_{j^{\ast }\left(t\right)}\right)$ for the distribution of minimum entropy values. In addition, in Figure 8a, the continuous blue line represents the weighted average trajectory calculated according to formula Eq. (21). Figure 8b, similar to Figure 8a, shows the trajectories of changes in the probabilities of the distribution of maximum entropy values. The jump-like red line corresponds to a change in time in the sequence of cluster centers $\left(X_{j^{\ast }\left(t\right)},Y_{j^{\ast }\left(t\right)}\right)$, in which for each position of the 3-year time window the maximum probability value $p_t^{\left(j\right)}$ of the histogram of the distribution of entropy maxima among all clusters of stations is realized. The continuous pink line in Figure 8b depicts the weighted average trajectory Eq. (21) of the maximum entropy values, which is similar to the blue line in Figure 8a for the minimum entropy values. Note that those clusters for which the probabilities in the histograms in Figure 6(a2) and Figure 6(b2) are small are not visited by the jump-like trajectories in Figure 8.

Figure 8.

In plot (a), the blue lines show the trajectory of changes in the numbers of clusters in which the maximum probability was realized, and the cyan line shows the corresponding weighted average trajectory of the distribution of the minimum entropy values over all clusters (Figure 5a). In plot (b), the red lines show the trajectory of changes in the numbers of clusters in which the maximum probability was realized, and the pink line shows the corresponding weighted average trajectory of the distribution of maximum entropy values over all clusters (Figure 5b). The thick brown line in (a) and (b) shows the san-Andreas fault. The largest cities in California, Los Angeles and San Francisco, are located in clusters #2 and #7, respectively.

From a comparison of Figure 8a and b, it is noticeable that the region of minimum entropy values gravitates toward the Pacific coast, in particular, toward the San Andreas fault, which is apparently due to more intense movements of crustal blocks in this region.

Another way to present trajectories of changes in the probabilities of extreme entropy values is presented in Figure 9, where the trajectories are plotted separately for latitude and longitude. For jump-like trajectories of changes in maximum probabilities along piecewise-constant graphs, numerical labels are placed in the appropriate color to indicate membership in one or another cluster of stations. In Figure 9a, on a continuous graph of changes in the weighted average probability (pink line), one can see the trend of migration of the probability of maximum entropy values to the north (in the direction of increasing latitude) for time stamps of the right end of the 3-year window until 2019, which is noticeable in Figure 7b. Some clusters of stations are visited by both the trajectory of maximum probabilities of minimum entropy values (blue line in Figure 8a) and the trajectory of maximum probabilities of maximum entropy values (red line in Figure 8b). However, these trajectories never coincide in time, as follows from the graph in Figure 9c, which shows the distance between these trajectories.

Figure 9.

In plots (a) and (b) by blue and cyan colors for minimum entropy values and by red and pink colors for maximum entropy values the changes in latitude (a) and longitude (b) are presented for maximum and weighted mean probabilities. Along the lines of changes in the coordinates of the maximum probabilities of extreme entropy values (blue for the minimum and red for the maximum), the corresponding colors mark the belonging of the trajectory sections to one or another cluster. Plot (c) presents distance between trajectories of maximum probabilities of extremum entropy values.

Based on the interpretation of the entropy of GPS time series as a measure of non-stationarity, from the graphs in Figure 9 we can conclude that for the last time fragment of 3 years in length, 2018–2023, the most “restless”, with the lowest entropy, is cluster #7, in which San-Francisco is located.

7. Spatiotemporal correlations of entropy values in station clusters

In a 3-year sliding time window, we calculate absolute values of all pairwise correlation coefficients between entropy values in each cluster. For 18 clusters, the number of such pairwise correlations is 153. Figure 10a shows a graph of the change in the average value of absolute pairwise correlations depending on the position of the right end of the sliding 3-year window. This graph highlights the 2017–2019 time interval as the time slice of maximum average correlations. Let us calculate the number of pairs of clusters between which an absolute correlation has arisen that exceeds the threshold of 0.6. The change in this number of pairs of clusters with the strongest correlations of entropy values is presented in Figure 10b, which also highlights the time interval 2017–2019 with the largest number of pairwise strong correlations. Finally, Figure 10c shows a plot of changes in the maximum distance between cluster centers for which the absolute entropy correlation exceeded 0.6. This graph contains two bursts of maximum distance, one of which (the second in a row) again falls on the time interval 2017–2019.

Figure 10.

(a) – mean value of absolute pair correlations between values of entropy of clusters within moving time window of the length 3 years; (b) – number of pairs of clusters for which absolute correlations exceed the threshold 0.6; (с) – maximum distance between clusters centers for which absolute correlations exceed 0.6.

In Figure 11, to visualize strong spatial correlations, straight lines connect the centers of those clusters for which the absolute values of correlations between entropy values in a sliding time window of 3 years length exceed the threshold of 0.6. The connection of cluster centers is carried out for three intervals of timestamp values at the right end of the 3-year time window. The time stamp boundaries are chosen so that each graph in Figure 11 corresponds to a time interval of the same length. Figure 11a and c correspond to the beginning and end of the analyzed time period. The central Figure 11b presents the strong correlations for the time interval corresponding to the maximum point of the average correlations and the numbers of pairs of clusters with strong correlations in Figure 10a and b. It can be seen from Figure 11b that strong correlations form a network that covers the entire region, while for weak average correlations in Figure 11a and c strong pairwise correlations are observed only between the southern clusters of stations and those located along the San-Andreas fault.

Figure 11.

Purple straight lines connect the centers of those clusters for which pairwise absolute correlations between entropy values in sliding time windows of 3 years in length exceed 0.6 for three intervals of timestamp values at the right end of the time windows.

Figure 12 provides a detailed view of the time evolution of correlations between clusters. For each cluster, Figure 12a shows graphs of changes in the numbers of pairs of other clusters for which the correlation between entropy values in a sliding time window of 3 years exceeded 0.6. The maximum number of such pairwise correlations is 3 - therefore, all Y-axes in these graphs have the same scale of 3. For each cluster, Figure 12b shows the probability that the pairwise correlation coefficient for entropy values with other clusters will exceed 0.6. From this graph, it is clear that clusters #9 and #11 are the least correlated with other clusters and act as a kind of border dividing the entire region into northern and southern parts. Figures 11 and 12 show that strong correlations are most often observed for cluster pairs located along the San Andreas Fault and on the Pacific Coast.

Figure 12.

Figure (a) for each cluster of stations shows graphs of the numbers of pairwise correlations of the principal components of this cluster with the principal components of other clusters exceeding the threshold of 0.6 in a sliding time window of 3 years. Figure (b) plots sample estimates of the probabilities that the correlations of the principal component of a given cluster with other clusters exceed the 0.6 threshold.

8. Conclusion

A method for joint analysis of a large number of synchronous monitoring time series is proposed using the example of data on daily vertical displacements of the earth’s surface in the Western USA (California) for the observation interval 2009–2023. The analysis is based on dividing the initial data (699 GPS time series) into blocks (clusters of stations) and sequentially applying the principal component method. The minimum normalized entropy of the distribution of squares of orthogonal wavelet coefficients, which is a measure of the information content of the distribution of oscillation energy over a system of independent time-frequency intervals, was chosen as the main research tool. Particular attention is paid to time intervals and areas (clusters) with extreme entropy values. Periodicities in the occurrence of strong jumps in the average entropy of the entire region have been discovered, of which the period of 95 days is dominant. The trend of migration of areas of maximum entropy from the South to the North has been identified. As a result of the analysis, it was found that the area of minimum entropy values gravitates toward the San Andreas fault, and the vicinity of San Francisco (cluster #7) has the selected properties of maximum information content and the last 3-year time fragment (2021–2023) of low entropy trajectories belongs to this cluster.

Further development of methods for analyzing the fine structure of GPS time series, going beyond the traditional analysis of the field of displacements of the earth’s surface, involves, in particular, a more detailed analysis of spectral properties. In particular, the transition from studying the field of maximum traditional Pearson correlations to the field of frequency-dependent correlations, that is, to the analysis of maximum coherences, seems promising. In this direction, it is also possible to use Hilbert-Huang expansions [57] in terms of a system of independent intrinsic vibration modes, which provides additional opportunities for studying nonlinear nonstationary effects.

Acknowledgments

The work was carried out within the framework of the state assignment of the Institute of Physics of the Earth of the Russian Academy of Sciences (topic FMWU-2022-0018).

Conflict of interest

The authors declare no conflict of interest.

Funding

This research received no external funding.