Next Article in Journal
APCSMA: Adaptive Personalized Client-Selection and Model-Aggregation Algorithm for Federated Learning in Edge Computing Scenarios
Next Article in Special Issue
Visibility Graph Investigation of the Shallow Seismicity of Lai Chau Area (Vietnam)
Previous Article in Journal
Optimizing Distributions for Associated Entropic Vectors via Generative Convolutional Neural Networks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Prognostic Properties of Instantaneous Amplitudes Maxima of Earth Surface Tremor

Institute of Physics of the Earth RAS, Moscow 123242, Russia
*
Author to whom correspondence should be addressed.
Entropy 2024, 26(8), 710; https://doi.org/10.3390/e26080710
Submission received: 25 July 2024 / Revised: 19 August 2024 / Accepted: 20 August 2024 / Published: 21 August 2024
(This article belongs to the Special Issue Time Series Analysis in Earthquake Complex Networks)

Abstract

:
A method is proposed for analyzing the tremor of the earth’s surface, measured by GPS, in order to highlight prognostic effects. The method is applied to the analysis of daily time series of vertical displacements in Japan. The network of 1047 stations is divided into 15 clusters. The Huang Empirical Mode Decomposition (EMD) is applied to the time series of the principal components from the clusters, with subsequent calculation of instantaneous amplitudes using the Hilbert transform. To ensure the stability of estimates of the waveforms of the EMD decomposition, 1000 independent additive realizations of white noise of limited amplitude were averaged before the Hilbert transform. Using a parametric model of the intensities of point processes, we analyze the connections between the instants of sequences of times of the largest local maxima of instantaneous amplitudes, averaged over the number of clusters and the times of earthquakes in the vicinity of Japan with minimum magnitude thresholds of 5.5 for the time interval 2012–2023. It is shown that the sequence of the largest local maxima of instantaneous amplitudes significantly more often precedes the moments of time of earthquakes (roughly speaking, has an “influence”) than the reverse “influence” of earthquakes on the maxima of amplitudes.

1. Introduction

This article presents the further development of methods for analyzing ground surface tremor proposed in [1,2,3]. The coherence of the tremor of the earth’s surface was analyzed in [4,5]. In [6,7], coherent analysis of GPS time series was used to assess seismic hazard in Japan and California.
In this case, the main emphasis is on the use of the Hilbert–Huang expansion, which is well suited to take into account the effects of nonstationarity and nonlinearity in time series [8,9]. This method has been successfully used to analyze geodetic time series [10,11], when processing hydrological [12], financial [13], biological [14,15] and seismic [16] data.
The main purpose of the article is to clarify common hypotheses that movements of the earth’s crust recorded by GNSS may contain predictive information. That displacements recorded by geodetic methods respond to the effects of earthquakes is widely known and has been demonstrated many times. But extracting geodetic effects that predict seismic events is much more challenging. In our paper, we propose one method for detecting predictive effects in space geodesy data.
The works [17,18,19,20] analyzed the composition of GPS time series—both their high-frequency part and low-frequency seasonal components—in connection with the estimates of the velocities of tectonic plates. In [21,22,23,24,25,26,27], using multivariate statistical methods, the influence of hydrological loads on tectonic displacements of sections of the earth’s crust was studied. The influence of time delays on the sensitivity of GPS solutions due to the impact of ionospheric and tropospheric factors was analyzed in [28]. The causes of the occurrence of “anomalous harmonics” in the spectral decomposition of GPS time series were considered in [29].
The detailed structure of earth surface displacements presented in GPS time series has been analyzed in the works of a large number of scientists. One of the most popular approaches is the use of the maximum likelihood method for estimating the parameters of GPS time series models [30,31,32,33]. This method was used in [34,35,36] to estimate the parameters of the power spectrum shape and noise amplitude for data from different regions of the world, and the error estimates are discussed in [37,38]. The influence of the spectrum shape and noise amplitude on the errors of displacement velocity estimates was investigated in [39,40]. Phase correlations of GPS time series were studied in [41,42] using parametric models for a number of tectonically active regions.
In [43,44], principal component analysis, empirical orthogonal function analysis, and singular spectrum analysis were used to determine the most common spatial and temporal components of GPS time series. In [45], a joint analysis of accelerometer readings and the noise component of GPS time series was performed. The influence of non-stationary effects on the estimates of relative displacements of crustal blocks and station positions was studied in [46,47,48].

2. Data

GPS data of daily earth displacements are taken from the Nevada Geodetic Laboratory website [49] at: http://geodesy.unr.edu/gps_timeseries/tenv3/IGS14/, access at 19 August 2024.
The set of 1047 GPS stations within domain 30 Lat 46 , 128 Lon 146 , were chosen. These stations presented in Figure 1a have daily time series from the beginning of 2009 to the end of 2023 (15 years), with a total number of gaps of less than 365 samples, and the longest gap less than 182 samples. The vertical components of ground displacement are investigated. The choice of only the vertical component for Japan is due to the fact that it does not contain the catastrophic jump and subsequent long-term relaxation due to the impact of the mega-earthquake of 11 March 2011. Gaps in the GPS time series are filled in using information from the left and right neighborhoods of the gap of the same length as the gap length [4].
The set of stations was previously divided into 15 clusters (Figure 1). The number 15 was chosen as the number of clusters, which optimally splits their “cloud”. Let us split the set of station position vectors ζ into a given number q of clusters using the k-means clustering method [50]. Let C r , r = 1 , , q be clusters, let z r = ζ C r ζ / n r be the vector of the center of mass of the cluster C r , and let n r be the number of vectors in the cluster, r = 1 q n r = N . Vector ζ C r if the distance | ζ z r | is the minimum among the positions of all cluster centers. The k-means method minimizes the sum of squared distances:
G ( z 1 , , z q ) = r = 1 q ζ C r | ζ z r | 2 min z 1 , , z q
relative to the position of the cluster centers z r . Let Φ ( q ) = min z 1 , , z q G ( z 1 , , z q ) . We used a trial number of clusters in the range 2 q 50 . The problem of choosing the best number q * of clusters was solved using the maximum pseudo-F-statistic criterion [51]
P F S ( q ) = σ 1 2 ( q ) / σ 0 2 ( q ) max 2 q 50
where
σ 0 2 ( q ) = Φ ( q ) N q , σ 1 2 ( q ) = r = 1 q n r N | z r z 0 | 2 , z 0 = 1 N 1 N ζ
The plot in Figure 1b presents the pseudo-F-statistic values as a function of the trial number of clusters. The number 15 on the pseudo-F-statistics graph is the inflection point of the dependence on the trial number of clusters and realizes one of the largest local maxima for the number of clusters from 2 to 50. On the pseudo-F-statistics graph, they represent two local maxima with close values of the number of clusters 6 and 15. Of these two values, 15 was chosen as the largest in order to provide the most detailed breakdown of the set of stations. Figure 1a shows the division of a set of stations into 15 clusters along with Voronoi cells, which indicate that the stations belong to a particular cluster.
Clusters of stations are ordered by increasing latitude of the position of their centers of gravity. Table 1 shows for each cluster (first row) the number of stations in the cluster (second row).

3. Principal Components of Increments in a Moving 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. Switching 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 a set of stations into 15 clusters is used for the subsequent application of the principal component method [52]. 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 there be a p-dimensional cloud of similar N-dimensional signals { y j ( k ) , k = 1 , , N } , j = 1 , , p . Let us choose the size of the sliding window w and center the signals,
x j ( k ) = y j ( k ) y j ( k ) ¯ , j = 1 , , p ,   k = w , , N ,
y j ( k ) ¯ = 1 w i = 1 w y j ( k w + i ) , k = w , , N .
The next step is to calculate the sample estimate of the covariance ( p × p ) -dimensional matrix in a sliding window:
r X X ( m , n ) = 1 w i = 1 w x m ( i w + k ) x n ( i w + k ) , k , j = 1 , , p ,   k = w , , N .
Let ϕ ( s ) = ( ϕ 1 ( s ) , , ϕ p ( s ) ) T be the eigenvector of this matrix corresponding to the maximum eigenvalue. Let us put
ξ s ( k ) = j = 1 p ( ϕ j ( k ) ) 2 x j ( k ) ,   k = w , , N .
Having generalized Formulas (4)–(7) with understandable changes to the case k = 1 , , w 1 , let us determine the weighted average in a sliding time window of length w using the formula
ξ ( k ) = { ξ w 1 ( k ) , k < w , ξ k ( k ) , k w .
Thus, Formulas (7) and (8) determine the values of the weighted average increments of vertical time series of displacements of the earth’s surface. The squared values of the eigenvector of the covariance matrix in the sliding time window corresponding to the largest eigenvalue are taken as weights. The sum of these weights is equal to one.
Within each of the 15 clusters, a transition to a weighted average was made using the method described above; the length w of the sliding time window was taken as equal to 365 samples, i.e., 1 year. At the same time, in order to eliminate the influence of large outliers, before calculating the weighted average, the so-called winsorization procedure was carried out [53], which consists in eliminating outliers falling beyond the level μ ± 4 σ by cutting off the values of the time series in a sliding time window ( μ and σ are sample estimates of the mathematical expectation and standard deviation for the current time window). The procedure is repeated iteratively until the values μ and σ stop changing.
Figure 2 shows graphs of the first principal components of increments (in the form of weighted averages) of vertical displacements of the earth’s surface in each of the identified 15 clusters.

4. Empirical Mode Decomposition

Let x ( k ) be the analyzed discrete signal. Empirical mode decomposition (EMD) [8,9] represents the decomposition of the signal into modes of oscillation:
x ( k ) = j = 1 n h j ( k ) + r n ( k )
where h j ( k ) is the j-th empirical mode, r n ( k ) is the remainder, n is the number of empirical modes.
The algorithm for decomposing into a sequence of empirical modes is iterative for each level j . Let us denote by m , m = 0 , 1 , , M j the index of iterations, where M j is the maximum number of iterations for level j . The iterations are described by the formula
h j ( m + 1 ) ( k ) = h j ( m ) ( k ) z j ( m ) ( k ) ,
Here z j ( m ) ( k ) = ( p j ( m ) ( k ) + q j ( m ) ( k ) ) / 2 , where p j ( m ) ( k ) and q j ( m ) ( k ) are both the upper and lower envelopes for the signal, which are constructed using spline interpolation (usually a third order spline) over all local maxima and minima of the signal h j ( m ) ( k ) .
Iterations (10) are initialized with step zero for the first level ( j = 1 ) of the expansion h 1 ( 0 ) ( k ) = x ( k ) . Next, the upper p 1 ( 0 ) ( k ) and lower q 1 ( 0 ) ( k ) envelopes are found, and the middle line z 1 ( 0 ) ( k ) is calculated and found h 1 ( 2 ) ( k ) using Formula (10). For h 1 ( 1 ) ( k ) , the upper p 1 ( 1 ) ( k ) and lower q 1 ( 1 ) ( k ) envelopes are determined, and the middle line z 1 ( 1 ) ( k ) is found, and so on, until the last iteration index M 1 , after which it is considered that the first empirical mode h 1 ( k ) has been found.
The condition for stopping iterations is usually chosen in the form of the following inequality:
k ( h j ( m + 1 ) ( k ) h j ( m ) ( k ) ) 2 / k ( h j ( m ) ( k ) ) 2 δ ,
where δ is some small number, for example 0.01. After the mode h j ( k ) is found, the iterative process of determining the empirical mode h j + 1 ( k ) of the next level is started. This process is initialized by the formula for the initial iteration index m = 0 :
h j + 1 ( 0 ) ( k ) = x ( k ) h j ( k ) ,
According to formula (12), the high-frequency part is subtracted from the original signal and a new, lower-frequency signal is considered as a new signal for subsequent decomposition. The construction of empirical oscillation modes continues until the number of local extrema becomes too small for envelopes to be constructed from them. As the number j of the empirical mode level increases, the signals become increasingly low-frequency and tend towards an unchanging form. The sequence h 1 ( k ) , h 2 ( k ) , , h n ( k ) is constructed in such a way that its sum gives an approximation to the original signal x ( k ) , which can be represented in formula (9) [8,9]. Empirical modes are orthogonal to each other, thus constituting a certain empirical basis for the decomposition of the original signal.
In the practical implementation of the method, technical difficulties arise due to edge effects, since the continuation of the envelopes beyond the first and last points of local extrema is ambiguous. To overcome this difficulty, there are several approaches, in particular mirror continuation [8,13] of the analyzed sample back and forth for a sufficiently long period of time. It was this approach that was used in this work.

5. Ensemble Empirical Mode Decomposition

One of the key disadvantages of the EMD method is the problem of mode mixing, which occurs when one empirical mode includes signals of different scales or when signals of the same scale are distributed over different empirical modes. For example, if “intermittency” is observed in the signal, that is, against the background of a smooth signal, short-term sections of higher frequency behavior appear, then with EMD decomposition behavior modes with different frequencies are mixed, since relatively rare points of local extrema of smooth behavior are interspersed with much more frequent points of local extrema of the high-frequency component.
To combat this effect, ref. [9] proposed the ensemble empirical mode decomposition (EEMD) method. This is a regularization of the EMD method in which white noise of finite amplitude is added to the original data. This allows the true empirical modes to be determined as the average over an ensemble of trials, each of which is the sum of signal and white noise.
The EEMD algorithm includes the following steps:
  • Add a white noise implementation to the original data.
  • Decomposition of data with the addition of white noise into empirical modes.
  • Repeat steps 1 and 2 quite a large number of times with different implementations of white noise.
  • Obtain the ensemble average for the corresponding empirical modes.
Thus, numerous “artificial” observations are simulated:
x ( i ) ( k ) = x ( k ) + ε i ( k ) ,
where ε i ( k ) is the i-th realization of white noise.
The true component, according to the EEMD definition, for a sequence of all levels of empirical modes is calculated as the average value of the expansions of the noisy modes (13). It is important to note that EEMD largely eliminates this mixing problem [9]. Adding independent white noise to the sample has a regularizing effect, since it simplifies the construction of envelopes (after adding a small amount of white noise, many local extrema are immediately created). The operation of averaging over a sufficiently large number of independent implementations of white noise makes it possible to eliminate the influence of the noise component and to isolate the true internal modes of oscillations.
For each of the resulting 15 time series, EEMD waveforms of principal components were calculated. EEMD waveforms are obtained by averaging 1000 decompositions of the original signals, to which are added independent Gaussian white noises with a standard deviation of 0.1 from the standard deviation of the weighted average from each cluster. Figure 3 shows EEMD waveform plots for the first 6 expansion levels for three of the 15 clusters.

6. Hilbert Transform

The Hilbert transform of the signal is determined by the formula [54]:
H X ( t ) = + X ( u ) π ( t u ) d u = X ( t ) ( 1 π t )
where f ( t ) g ( t ) = + f ( u ) g ( t u ) d u is the convolution of two functions. If f ˜ ( ω ) and g ˜ ( ω ) are the Fourier transforms of convolutional functions f ˜ ( ω ) = + f ( t ) e i ω t d t , then, as is known, the Fourier transform of convolution is equal to the product of the Fourier transforms of convolutional functions. Fourier transform from 1 / ( π t ) equals:
+ e i ω t π t d t = i s i g n ( ω ) = { i , ω > 0 0 , ω = 0 i , ω < 0
Thus, if X ˜ ( ω ) there is a Fourier transform of X ( t ) , then
H ˜ X ( ω ) = i s i g n ( ω ) X ˜ ( ω ) = { i X ˜ ( ω ) , ω > 0 0 , ω = 0 i X ˜ ( ω ) , ω < 0
If you present X ˜ ( ω ) = | X ˜ ( ω ) | e i φ ( ω ) , then
H ˜ X ( ω ) = { | X ˜ ( ω ) | e i ( φ ( ω ) + π / 2 ) , ω > 0 0 , ω = 0 | X ˜ ( ω ) | e i ( φ ( ω ) π / 2 ) , ω < 0
In practice, it is more convenient to calculate the Hilbert transform using the concepts of an analytical signal:
Z X ( t ) = X ( t ) + i H X ( t ) = | Z X ( t ) | e i ϑ ( t ) A X ( t ) e i ϑ ( t )
where A X ( t ) = X 2 ( t ) + H X 2 ( t ) are the amplitudes of the signal X ( t ) envelope, and ϑ ( t ) is the instantaneous phase. The derivative ν ( t ) = d ϑ ( t ) / d t is called instantaneous frequency. The Fourier transform of the analytical signal Z ˜ X ( ω ) = X ˜ ( ω ) + i H ˜ X ( ω ) = X ˜ ( ω ) ( 1 + s i g n ( ω ) ) or:
Z ˜ X ( ω ) = { 2 X ˜ ( ω ) , ω > 0 X ˜ ( 0 ) , ω = 0 0 , ω < 0
after which the Hilbert transform is equal to the imaginary part of the result of the inverse Fourier transform of Z ˜ X ( ω )
H X ( t ) = Im ( 1 2 π + Z ˜ X ( ω ) e i ω t d ω )
For each of the resulting 15 time series, EEMD waveforms of principal components were calculated. EEMD waveforms are obtained by averaging 1000 decompositions of the original signals, to which are added independent Gaussian white noises with a standard deviation of 0.1 from the standard deviation of the weighted average from each cluster. Figure 4 shows EEMD waveform plots for the first 6 expansion levels for three of the 15 clusters.
For a discrete-time signal X ( t ) , t = 0 , , ( N 1 ) , this transformation can be calculated using the discrete Fourier transform:
d X ( N ) ( ω k ) = t = 0 N 1 X ( t ) exp ( i ω k t ) , ω k = 2 π N ( k 1 ) , k = 0 , 1 , , ( N 1 )
after which the second part of the Fourier coefficients (corresponding to negative frequencies) should be reset to zero: h X ( N ) ( ω k ) = 0 , k = N / 2 + 1 , , ( N 1 ) while the 1st part should be doubled: h X ( N ) ( ω k ) = 2 d X ( N ) ( ω k ) , k = 1 , , N / 2 . The Hilbert transform is then calculated as the imaginary part of the inverse discrete Fourier transform:
H X ( t ) = Im ( k = 0 N 1 h X ( N ) ( ω k ) exp ( i ω k t ) / N ) , t = 0 , 1 , , ( N 1 )
Thus, after determining the instantaneous amplitudes and frequencies of the EEMD, the Hilbert–Huang decomposition can be represented as follows:
x ( t ) = j = 1 n h j ( t ) + r n ( t ) , h j ( t ) = Re { A j ( t ) exp ( i ν j ( s ) d s ) }

7. Influence Matrix

To solve the problem of estimating the connection between two sequences of random events, a parametric intensity model is used. In [55,56,57] this method was used to test the hypotheses that local extrema of the average values of certain properties of seismic noise and magnetic field precede the instants of strong earthquakes.
Let t j ( α ) , j = 1 , , N α ; α = 1 , 2 represent the moments in time of 2 streams of events.
In our case it is:
(1)
a sequence of moments in time corresponding to the largest local maxima of the amplitudes of the envelopes at certain levels of the EEMD Huang decomposition
(2)
a sequence of times of seismic events with a magnitude not less than a given value.
Let us present their intensities in the form:
λ ( α ) ( t ) = b 0 ( α ) + β = 1 2 b β ( α ) g ( β ) ( t )
where b 0 ( α ) 0 , b β ( α ) 0 are parameters, and g ( β ) ( t ) is the influence function of flow events t j ( β ) with number β :
g ( β ) ( t ) = t j ( β ) < t exp ( ( t t j ( β ) ) / τ )
According to formula (25), the weight of an event with number j becomes non-zero for times t > t j ( β ) and decays with a characteristic time τ . The parameter b β ( α ) determines the degree of influence of the sequence β on the sequence α . The parameter b α ( α ) determines the degree of self-excitation to which the flow α influences itself, and the parameter b 0 ( α ) reflects the purely random (Poisson) intensity component. Let us fix the parameter τ and consider the problem of estimating the parameters b 0 ( α ) , b β ( α ) .
The log-likelihood function for a non-stationary Poisson process within the time interval [ 0 , T ] equals [58]:
ln ( L α ) = j = 1 N α ln ( λ ( α ) ( t j ( α ) ) ) 0 T λ ( α ) ( s ) d s , α = 1 , 2
We need to find parameters b 0 ( α ) , b β ( α ) from maximum of functions (26). It is easy to obtain the following expression:
b 0 ( α ) ln ( L α ) b 0 ( α ) + β = 1 2 b β ( α ) ln ( L α ) b β ( α ) = N α 0 T λ ( α ) ( s ) d s
The parameters b 0 ( α ) , b β ( α ) are non-negative. It means that each term on the left side of (27) equals to zero at point of maximum of function (26)—either due to the necessary conditions for the extremum (if the parameters are positive), or, if the maximum is reached at the boundary, then the parameters themselves are equal to zero. Therefore, at the maximum point of the likelihood function the equality holds:
0 T λ ( α ) ( s ) d s = N α
Let us substitute the expression g ( β ) ( t ) from (28) into (27) and divide by T . Then we get another form of formula (28):
b 0 ( α ) + β = 1 m b β ( α ) g ¯ ( β ) = λ 0 ( α ) N α / T
where
g ¯ ( β ) = 0 T g ( β ) ( s ) d s / T
average value of the influence function. Substituting b 0 ( α ) from (29) into (26), we obtain the following maximum problem:
Φ ( α ) ( b 1 ( α ) , b 2 ( α ) ) = j = 1 N α ln ( λ 0 ( α ) + β = 1 2 b β ( α ) Δ g ( β ) ( t j ( α ) ) ) max
where Δ g ( β ) ( t ) = g ( β ) ( t ) g ¯ ( β ) , under restrictions:
b 1 ( α ) 0 , b 2 ( α ) 0 , β = 1 2 b β ( α ) g ¯ ( β ) λ 0 ( α )
It could be shown that function (31) is convex with a negative definite Hessian. Therefore, problem (31) and (32) has a unique solution. The problem (31) and (32) is solved numerically for a given relaxation parameter τ . After this step we can define the influence matrix with elements κ β ( α ) , α = 1 , 2 ; β = 0 , 1 , 2 using the formulas:
κ 0 ( α ) = b 0 ( α ) λ 0 ( α ) 0 , κ β ( α ) = b β ( α ) g ¯ ( β ) λ 0 ( α ) 0
The value κ 0 ( α ) is the share of the average intensity λ 0 ( α ) of the process with number α , which is purely random, part κ α ( α ) is caused by the influence of self-excitation α α , and κ β ( α ) , β α is due to external influence β α . From Formula (29) the normalization condition follows:
κ 0 ( α ) + β = 1 2 κ β ( α ) = 1 , α = 1 , 2
As a result, the influence matrix can be determined:
( | κ 0 ( 1 ) κ 0 ( 2 ) | | κ 1 ( 1 ) κ 2 ( 1 ) κ 1 ( 2 ) κ 2 ( 2 ) | )
The first column of matrix (35) is composed of Poisson shares of average intensities. The diagonal elements of the right submatrix of size 2 × 2 consist of self-exciting elements of mean intensity, while the off-diagonal elements correspond to mutual excitation. The sums of the component rows of the influence matrix (34) equal 1. The influence matrices are estimated in a certain moving time window of length L with an offset Δ L and with a given value of the relaxation parameter τ .

8. Estimation of Connections between the Times of Local Amplitude Maxima and Seismic Events

The further plan of the article is to use the apparatus of influence matrices to assess the relationship between the times of maximum average amplitudes of the envelopes and the times of sufficiently strong earthquakes. A magnitude threshold of 5.5 was chosen. There were 673 such seismic events in the vicinity of the Japanese islands during the period of time 2009–2023—see Figure 5a. However, the mega-earthquake of 11 March 2011 with a magnitude of 9 gave rise to a surge in aftershock activity, as a result of which, if we consider the time interval 2012–2023, when the intensity of aftershocks had already decreased, the number of earthquakes with a magnitude of at least 5.5 will decrease by almost 2 times to 349—see Figure 5b. It should be noted that accurately estimating the time of the end of the Tohoku earthquake aftershocks is a difficult task [59] and in this case we used a rough estimate based on the visual perception of the intensity of seismic events.
The working hypothesis is that for certain levels of EEMD decomposition, the times of the largest local maxima of the average amplitudes of the envelopes precedes the times of earthquakes. For a correct comparison of two streams of events, it is necessary that their average intensity be approximately equal. This means that the number of the largest local maxima of the amplitudes of the envelopes during the time period under study should be equal to the number of earthquakes. From these considerations, it becomes clear that for a correct analysis of the connections between the time instants of local maximum amplitudes of envelopes and seismic events, the time interval of aftershock activity must be excluded. Therefore, further analysis is carried out for the time period 2012–2023 lasting 12 years.
Figure 6 shows the distribution of epicenters of earthquakes with a magnitude of at least 5.5.
It should be taken into account that as the number of the decomposition level increases, both the waveforms of the levels themselves and the amplitudes of their envelopes become increasingly low-frequency. As a result, it is possible to select the 349 largest local maxima of average amplitudes in the interval 2012–2023 only for a certain number of lower expansion levels. For the time interval 2012–2023, only the first two levels are suitable for selecting 349 local maxima of amplitudes, since the number of local extrema already at the third level is 242, that is, less than 349. From the first two levels, it was decided to choose the second as the lower frequency, and for which the results of mutual influence assessments are more expressive.
In Figure 7, the red dots represent the selected 349 largest local maxima of the average amplitude of the envelopes at the second level of decomposition.
Let us call the “direct” influence of the moments of time of earthquakes on the occurrence of local maxima of the average amplitudes of the envelopes, and the “reverse” —correspondingly, the advance of the moments of time of local maxima of amplitudes relative to the times of earthquakes. Figure 8 shows graphs of changes in the components of the matrix of “direct” and “reverse” influence for level 2 when assessed in a sliding window of 2 years.
Of these graphs, the pair (a1, a2) is of greatest importance: a1 represents the change in the components of the direct influence of seismic events on the positions of local amplitude maxima, while a2 represents the inverse influence of the times of local amplitude maxima on seismic events. From a comparison of these two graphs, it is clear that the reverse influence significantly exceeds the direct influence; that is, there is a delay effect of seismic events relative to the maximum amplitudes. In other words, there is a predictive effect. Graphs (b1, b2) represent changes in the self-exciting component of average intensities, while graphs (c1, c2) represent changes in the purely random (Poisson) component. Finally, the pair of plots (d1, d2) represents the change in the number of jointly analyzed time points in each time window. Let us recall once again that the sum of the components (a1, b1, c1) and (a2, b2, c2) is equal to 1 for any position of the time window.
Figure 8 presents the results of estimates of the mutual influence of two sequences of events only for a time window of 2 years. In order to increase the representativeness of this result, we will carry out similar estimates for a sufficiently large set of time lengths varying within specified limits. In this case, for each value of the length of the time window, we will identify local maxima of the components of the influence matrix, which are responsible for the mutual influence of sequences of events when the time windows are shifted. Let us describe a method for constructing a set of maximum components of mutual influence matrices in the form of numbered points.
  • The minimum L min and maximum L max lengths of time windows and N L —the number of lengths of time windows in this interval are selected. Thus, the lengths of the time windows took on the values L k = L min + ( k 1 ) Δ L , k = 1 , , N L , Δ L = ( L max L min ) / ( N L 1 ) . In our calculations, we took L min as equal to 1 year, and L max —3 years, N L = 100 .
  • Each time window of length L k was shifted from left to right along the time axis with some offset Δ t . Let us denote by t j ( L k ) , j = 1 , , M ( L k ) the sequence of moments in time of the positions of the right windows with length L k . The number M ( L k ) of time windows in length L k is determined by their time offset Δ t . We used a time window offset Δ t of 0.01 year.
  • For each position of a time window of length L k , the elements of the influence matrix (35) are estimated for a given relaxation time τ of the model (26–27), corresponding to the mutual influence of the two processes being analyzed. We took a value τ equal to 0.1 year. For definiteness, we will consider one influence, for example, of the first process on the second. As a result of such estimates, we obtain their values in the form ( t j ( L k ) , c j ( L k ) ) , where c j ( L k ) is the corresponding element of the influence matrix for a position with a time window number j of length L k .
  • In the sequence ( t j ( L k ) , c j ( L k ) ) , we select elements ( t j * ( L k ) , c j * ( L k ) ) corresponding to local maxima of values c j ( L k ) , that is, from the condition c j 1 ( L k ) < c j * ( L k ) < c j + 1 ( L k ) . Let us present each element ( t j * ( L k ) , c j * ( L k ) ) as a vertical segment of length c j * ( L k ) located at a time point t j * ( L k ) . The combination of such vertical graphic elements for all k = 1 , , N L , j = 1 , , M ( L k ) visualizes the “strength” of the mutual influence of processes on each other.
So, the full set of free parameters of the method: τ , L min , L max , N L , Δ t . The result of such estimates is presented in Figure 9.
The results presented in Figure 9 confirm the conclusions made earlier based on the graphs in Figure 8: the “reverse” influence of the time instants of local maxima of envelope amplitudes at the second level on the time instants of earthquakes is significantly greater (Figure 9b) than the “direct” influence of earthquakes on the occurrence of local maxima in the amplitudes of the envelopes (Figure 9a).

9. Conclusions

Traditional methods of analyzing data on crustal movements obtained using space geodesy are focused on identifying systematic low-frequency components that can be interpreted as manifestations of slow tectonic movements. The high-frequency component of these time series, which can be called the “tremor” of the earth’s surface, is most often interpreted as a manifestation of noise arising from atmospheric and ionospheric fluctuations. Our point is that, despite the presence of this process noise, it is in the high-frequency component of GPS data that there is hidden prognostic information. In this article, the joint use of cluster analysis, principal component analysis, Hilbert–Huang decomposition and evaluation of parametric models of the mutual influence of sequences of events allowed us to obtain a result confirming the presence of predictive information in high-frequency tremor of the earth’s surface.
The authors used software which was written by them using programming languages Fortran and Python.

Author Contributions

A.L.—idea, programming, paper preparation; E.R.—programming, data preparation. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

The open access data from the source: http://geodesy.unr.edu/gps_timeseries/tenv3/IGS14/ (accessed on 15 January 2024) were used.

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).

Conflicts of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

  1. Lyubushin, A. Identification of Areas of Anomalous Tremor of the Earth’s Surface on the Japanese Islands According to GPS Data. Appl. Sci. 2022, 12, 7297. [Google Scholar] [CrossRef]
  2. Lyubushin, A. Singular Points of the Tremor of the Earth’s Surface. Appl. Sci. 2023, 13, 10060. [Google Scholar] [CrossRef]
  3. Lyubushin, A. Entropy of GPS-measured Earth tremor. In Revolutionizing Earth Observation—New Technologies and Insights; Abdalla, R.M., Ed.; IntechOpen: Rijeka, Croatia, 2024. [Google Scholar] [CrossRef]
  4. Lyubushin, A. Global coherence of GPS-measured high-frequency surface tremor motions. GPS Solut. 2018, 22, 116. [Google Scholar] [CrossRef]
  5. Lyubushin, A. Field of coherence of GPS-measured earth tremors. GPS Solut. 2019, 23, 120. [Google Scholar] [CrossRef]
  6. Filatov, D.M.; Lyubushin, A.A. Fractal analysis of GPS time series for early detection of disastrous seismic events. Phys. A Stat. Mech. Its Appl. 2017, 469, 718–730. [Google Scholar] [CrossRef]
  7. Filatov, D.M.; Lyubushin, A.A. Precursory Analysis of GPS Time Series for Seismic Hazard Assessment. Pure Appl. Geophys. 2019, 177, 509–530. [Google Scholar] [CrossRef]
  8. Huang, N.E.; Shen, Z.; Long, S.R.; Wu, V.C.; Shih, H.H.; Zheng, Q.; Yen, N.C.; Tung, C.C.; Liv, H.H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. Ser. A 1998, 454, 903–995. [Google Scholar] [CrossRef]
  9. Huang, N.E.; Wu, Z. A review on Hilbert-Huang transform: Method and its applications to geophysical studies. Rev. Geophys. 2008, 46, RG2006. [Google Scholar] [CrossRef]
  10. Pan, Y.; Shen, W.-B.; Ding, H.; Hwang, C.; Li, J.; Zhang, T. The Quasi-Biennial Vertical Oscillations at Ghlobal GPS Stations: Identification by Ensemble Empirical Mode Decomposition. Sensors 2015, 15, 26096–26114. [Google Scholar] [CrossRef]
  11. Li, W.; Guo, J. Extraction of periodic signals in Global Navigation Satellite System (GNSS) vertical coordinate time series using the adaptive ensemble empirical modal decomposition method. Nonlin. Process. Geophys. 2024, 31, 99–113. [Google Scholar] [CrossRef]
  12. Huang, Y.; Schmitt, F.G.; Lu, Z.; Liu, Y. Analysis of daily river flow fluctuations using empirical mode decomposition and arbitrary order Hilbert spectral analysis. J. Hydrol. 2009, 373, 103–111. [Google Scholar] [CrossRef]
  13. Huang, N.E.; Wu, M.; Qu, W.; Long, S.R.; Shen, S.S.P. Applications of Hilbert–Huang transform to non-stationary financial time series analysis. Appl. Stoch. Models Bus. Ind. 2003, 19, 245–268. [Google Scholar] [CrossRef]
  14. Li, H.; Kwong, S.; Yang, L.; Huang, D.; Xiao, D. Hilbert-Huang Transform for Analysis of Heart Rate Variability in Cardiac Health. IEEE/ACM Trans. Comput. Biol. Bioinform. 2011, 8, 1557–1567. [Google Scholar] [CrossRef]
  15. Wei, H.-C.; Xiao, M.-X.; Chen, H.-Y.; Li, Y.-Q.; Wu, H.-T.; Sun, C.-K. Instantaneous frequency from Hilbert-Huang transformation of digital volume pulse as indicator of diabetes and arterial stiffness in upper-middle-aged subjects. Sci. Rep. 2018, 8, 15771. [Google Scholar] [CrossRef] [PubMed]
  16. Sarlis, N.V.; Skordas, E.S.; Mintzelas, A.; Papadopoulou, K.A. Micro-scale, mid-scale, and macro-scale in global seismicity identified by empirical mode decomposition and their multifractal characteristics. Sci. Rep. 2018, 8, 9206. [Google Scholar] [CrossRef]
  17. Beavan, J. Noise properties of continuous GPS data from concrete pillar geodetic monuments in New Zealand and comparison with data from U.S. deep drilled braced monuments. J. Geophys. Res. 2005, 110, B08. [Google Scholar] [CrossRef]
  18. Langbein, J. Noise in GPS displacement measurements from Southern California and Southern Nevada. J. Geophys. Res. 2008, 113, B05405. [Google Scholar] [CrossRef]
  19. Blewitt, G.; Lavallee, D. Effects of annual signal on geodetic velocity. J. Geophys. Res. 2002, 107, 2145. [Google Scholar] [CrossRef]
  20. Bos, M.S.; Fernandes, R.M.S.; Williams, S.D.P.; Bastos, L. Fast error analysis of continuous GPS observations. J. Geod. 2008, 82, 157–166. [Google Scholar] [CrossRef]
  21. Liu, B.; Xing, X.; Tan, J.; Xia, Q. Modeling Seasonal Variations in Vertical GPS Coordinate Time Series Using Independent Component Analysis and Varying Coefficient Regression. Sensors 2020, 20, 5627. [Google Scholar] [CrossRef]
  22. Liu, B.; Dai, W.; Liu, N. Extracting seasonal deformations of the Nepal Himalaya region from vertical GPS position time series using Independent Component Analysis. Adv. Space Res. 2017, 60, 2910–2917. [Google Scholar] [CrossRef]
  23. Tesmer, V.; Steigenberger, P.; van Dam, T.; Mayer-Gurr, T. Vertical deformations from homogeneously processed GRACE and global GPS long-term series. J. Geod. 2011, 85, 291–310. [Google Scholar] [CrossRef]
  24. Yan, J.; Dong, D.; Burgmann, R.; Materna, K.; Tan, W.; Peng, Y.; Chen, J. Separation of Sources of Seasonal Uplift in China Using Independent Component Analysis of GNSS Time Series. J. Geophys. Res. Solid Earth 2019, 124, 11951–11971. [Google Scholar] [CrossRef]
  25. Santamaría-Gómez, A.; Memin, A. Geodetic secular velocity errors due to interannual surface loading deformation. Geophys. J. Int. 2015, 202, 763–767. [Google Scholar] [CrossRef]
  26. Fu, Y.; Argus, D.F.; Freymueller, J.T.; Heflin, M.B. Horizontal motion in elastic response to seasonal loading of rain water in the Amazon Basin and monsoon water in Southeast Asia observed by GPS and inferred from GRACE. Geophys. Res. Lett. 2013, 40, 6048–6053. [Google Scholar] [CrossRef]
  27. Chanard, K.; Avouac, J.P.; Ramillien, G.; Genrich, J. Modeling deformation induced by seasonal variations of continental water in the Himalaya region: Sensitivity to Earth elastic structure. J. Geophys. Res. Solid Earth 2014, 119, 5097–5113. [Google Scholar] [CrossRef]
  28. He, X.; Montillet, J.P.; Fernandes, R.; Bos, M.; Yu, K.; Hua, X.; Jiang, W. Review of current GPS methodologies for producing accurate time series and their error sources. J. Geodyn. 2017, 106, 12–29. [Google Scholar] [CrossRef]
  29. Ray, J.; Altamimi, Z.; Collilieux, X.; van Dam, T. Anomalous harmonics in the spectra of GPS position estimates. GPS Solut. 2008, 12, 55–64. [Google Scholar] [CrossRef]
  30. Roncagliolo, P.A.; García, J.G.; Mercader, P.I.; Fuhrmann, D.R.; Muravchik, C.H. Maximum-likelihood attitude estimation using GPS signals. Digit. Signal Process. 2007, 17, 1089–1100. [Google Scholar] [CrossRef]
  31. Wang, F.; Li, H.; Lu, M. GNSS Spoofing Detection and Mitigation Based on Maximum Likelihood Estimation. Sensors 2017, 17, 1532. [Google Scholar] [CrossRef]
  32. Parkinson, B.W. Global Positioning System: Theory and Applications; AIAA: Reston, VA, USA, 1996; p. 781. [Google Scholar]
  33. Xu, J.; He, J.; Zhang, Y.; Xu, F.; Cai, F. A Distance-Based Maximum Likelihood Estimation Method for Sensor Localization in Wireless Sensor Networks. Int. J. Distrib. Sens. Netw. 2016, 12, 2080536. [Google Scholar] [CrossRef]
  34. Langbein, J.; Johnson, H. Correlated errors in geodetic time series, Implications for time-dependent deformation. J. Geophys. Res. 1997, 102, 591–603. [Google Scholar] [CrossRef]
  35. Williams, S.D.P.; Bock, Y.; Fang, P.; Jamason, P.; Nikolaidis, R.M.; Prawirodirdjo, L.; Miller, M.; Johnson, D.J. Error analysis of continuous GPS time series. J. Geophys. Res. 2004, 109, B03412. [Google Scholar] [CrossRef]
  36. Wang, W.; Zhao, B.; Wang, Q.; Yang, S. Noise analysis of continuous GPS coordinate time series for CMONOC. Adv. Space Res. 2012, 49, 943–956. [Google Scholar] [CrossRef]
  37. Agnew, D. The time domain behavior of power law noises. Geophys. Res. Lett. 1992, 19, 333–336. [Google Scholar] [CrossRef]
  38. Amiri-Simkooei, A.R.; Tiberius, C.C.J.M.; Teunissen, P.J.G. Teunissen Assessment of noise in GPS coordinate time series: Methodology and results. J. Geophys. Res. 2007, 112, B07413. [Google Scholar] [CrossRef]
  39. Caporali, A. Average strain rate in the Italian crust inferred from a permanent GPS network—I. Statistical analysis of the time-series of permanent GPS stations. Geophys. J. Int. 2003, 155, 241–253. [Google Scholar] [CrossRef]
  40. Zhang, J.; Bock, Y.; Johnson, H.; Fang, P.; Williams, S.; Genrich, J.; Wdowinski, S.; Behr, J. Southern California permanent GPS geodetic array: Error analysis of daily position estimates and site velocities. J. Geophys. Res. 1997, 102, 18035–18055. [Google Scholar] [CrossRef]
  41. Li, J.; Miyashita, K.; Kato, T.; Miyazaki, S. GPS time series modeling by autoregressive moving average method, Application to the crustal deformation in central Japan. Earth Planets Space 2000, 52, 155–162. [Google Scholar] [CrossRef]
  42. Kermarrec, G.; Schon, S. On modelling GPS phase correlations: A parametric model. Acta Geod. Geophys. 2018, 53, 139–156. [Google Scholar] [CrossRef]
  43. Teferle, F.N.; Williams, S.D.P.; Kierulf, H.P.; Bingley, R.M.; Plag, H.P. A continuous GPS coordinate time series analysis strategy for high-accuracy vertical land movements. Phys. Chem. Earth Parts A/B/C 2008, 33, 205–216. [Google Scholar] [CrossRef]
  44. Chen, Q.; van Dam, T.; Sneeuw, N.; Collilieux, X.; Weigelt, M.; Rebischung, P. Singular spectrum analysis for modeling seasonal signals from GPS time series. J. Geodyn. 2013, 72, 25–35. [Google Scholar] [CrossRef]
  45. Bock, Y.; Melgar, D.; Crowell, B.W. Real-Time Strong-Motion Broadband Displacements from Collocated GPS and Accelerometers. Bull. Seismol. Soc. Am. 2011, 101, 2904–2925. [Google Scholar] [CrossRef]
  46. Goudarzi, M.A.; Cocard, M.; Santerre, R.; Woldai, T. GPS interactive time series analysis software. GPS Solut. 2013, 17, 595–603. [Google Scholar] [CrossRef]
  47. Hackl, M.; Malservisi, R.; Hugentobler, U.; Jiang, Y. Velocity covariance in the presence of anisotropic time correlated noise and transient events in GPS time series. J. Geodyn. 2013, 72, 36–45. [Google Scholar] [CrossRef]
  48. Khelif, S.; Kahlouche, S.; Belbachir, M.F. Analysis of position time series of GPS-DORIS co-located stations. Int. J. Appl. Earth Observ. Geoinf. 2013, 20, 67–76. [Google Scholar] [CrossRef]
  49. Blewitt, G.; Hammond, W.C.; Kreemer, C. Harnessing the GPS data explosion for interdisciplinary science. Eos 2018, 99, 485. [Google Scholar] [CrossRef]
  50. Duda, R.O.; Hart, P.E.; Stork, D.G. Pattern Classification; Wiley: Hoboken, NJ, USA, 2000. [Google Scholar]
  51. Vogel, M.A.; Wong, A.K.C. PFS clustering method. IEEE Trans. Pattern Anal. Mach. Intell. 1979, 1, 237–245. [Google Scholar] [CrossRef] [PubMed]
  52. Jolliffe, I.T. Principal Component Analysis; Springer: Berlin/Heidelberg, Germany, 1986. [Google Scholar]
  53. Huber, P.J. Robust Statistics; Wiley: Toronto, ON, Canada; Chichester, UK; New York, NY, USA, 1981. [Google Scholar]
  54. Bendat, J.S.; Piersol, A.G. Random Data. Analysis and Measurement Procedures, 4th ed.; Wiley & Sons: Hoboken, NJ, USA, 2010. [Google Scholar]
  55. Lyubushin, A. Investigation of the Global Seismic Noise Properties in Connection to Strong Earthquakes. Front. Earth Sci. 2022, 10, 905663. [Google Scholar] [CrossRef]
  56. Lyubushin, A. Seismic Hazard Indicators in Japan based on Seismic Noise Properties. J. Earth Environ. Sci. Res. 2023, 5, 1–8. [Google Scholar] [CrossRef]
  57. Lyubushin, A.; Rodionov, E. Wavelet-based correlations of the global magnetic field in connection to strongest earthquakes. Adv. Space Res. 2024. [Google Scholar] [CrossRef]
  58. Cox, D.R.; Lewis, P.A.W. The Statistical Analysis of Series of Events; Methuen: London, UK, 1966. [Google Scholar]
  59. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S.; Toshiyasu, N.; Masashi, K. The unusual case of the ultra-deep 2015 Ogasawara earthquake (Mw7.9): Natural time analysis. EPL 2021, 135, 49002. [Google Scholar] [CrossRef]
Figure 1. (a) shows the positions of 1047 GPS stations and their division into 15 clusters. The numbered circles indicate the centers of gravity of the clusters, and the blue lines indicate the boundaries between Voronoi cells. The blue star shows the position of the center of mass of all cluster centers. (b) shows a plot of the pseudo-F-statistic that allowed us to select 15 as the number of clusters.
Figure 1. (a) shows the positions of 1047 GPS stations and their division into 15 clusters. The numbered circles indicate the centers of gravity of the clusters, and the blue lines indicate the boundaries between Voronoi cells. The blue star shows the position of the center of mass of all cluster centers. (b) shows a plot of the pseudo-F-statistic that allowed us to select 15 as the number of clusters.
Entropy 26 00710 g001
Figure 2. Graphs of weighted average vertical displacements of the earth’s surface in each of the selected 15 clusters in a sliding time window of 365 days. The Y axes represent displacement increments in mm.
Figure 2. Graphs of weighted average vertical displacements of the earth’s surface in each of the selected 15 clusters in a sliding time window of 365 days. The Y axes represent displacement increments in mm.
Entropy 26 00710 g002
Figure 3. EEMD waveform plots for the first 6 decomposition levels for 3 clusters (numbers 1, 9 and 15). Decomposition level numbers are indicated on the left.
Figure 3. EEMD waveform plots for the first 6 decomposition levels for 3 clusters (numbers 1, 9 and 15). Decomposition level numbers are indicated on the left.
Entropy 26 00710 g003
Figure 4. Plots of mean EEMD waveforms and mean instantaneous amplitudes averaged over all 15 station clusters for the first 6 decomposition levels. The left side of the figure, separated by a vertical blue line, shows pairs of graphs (waveforms—their amplitudes) for levels 1–3, and the right side shows graphs for levels 4–6. The decomposition level numbers are indicated between the waveform and amplitude graphs.
Figure 4. Plots of mean EEMD waveforms and mean instantaneous amplitudes averaged over all 15 station clusters for the first 6 decomposition levels. The left side of the figure, separated by a vertical blue line, shows pairs of graphs (waveforms—their amplitudes) for levels 1–3, and the right side shows graphs for levels 4–6. The decomposition level numbers are indicated between the waveform and amplitude graphs.
Entropy 26 00710 g004
Figure 5. Time sequence of earthquakes with a magnitude of at least 5.5 in the vicinity of the Japanese Islands: (a) in the time period 2009–2023; (b) in the period of time 2012–2023.
Figure 5. Time sequence of earthquakes with a magnitude of at least 5.5 in the vicinity of the Japanese Islands: (a) in the time period 2009–2023; (b) in the period of time 2012–2023.
Entropy 26 00710 g005
Figure 6. Distribution of epicenters of 349 earthquakes with a magnitude of at least 5.5 in the vicinity of the Japanese Islands in the time period 2012–2023. The red asterisk marks the center of gravity of the centers of 15 clusters of GPS stations (Figure 1a), which is chosen as the center of a circle with a radius of 1500 km.
Figure 6. Distribution of epicenters of 349 earthquakes with a magnitude of at least 5.5 in the vicinity of the Japanese Islands in the time period 2012–2023. The red asterisk marks the center of gravity of the centers of 15 clusters of GPS stations (Figure 1a), which is chosen as the center of a circle with a radius of 1500 km.
Entropy 26 00710 g006
Figure 7. Average instantaneous amplitudes of the second level of EEMD decomposition and the 349 largest local maxima (red dots).
Figure 7. Average instantaneous amplitudes of the second level of EEMD decomposition and the 349 largest local maxima (red dots).
Entropy 26 00710 g007
Figure 8. Graphs of changes in the components of the influence matrix between sequences of seismic events with a magnitude of at least 5.5 and the sequence of moments in time of 349 of the largest local maxima of the average amplitudes at the second level of EEMD decomposition. The estimates were made in a sliding time window of length 2 years with a shift of 0.01 year for a relaxation time τ of the model (24, 25) of 0.1 year. The graphs (a1c1) refer to the components of the influence matrix (35), which refers to the intensity fractions of the sequence of the largest local amplitude maxima, while the graphs (a2c2) refer to the intensity fractions of the sequence of seismic events. Plots (d1) and (d2) present the numbers of local maxima of amplitudes (d1) and the number of seismic events (d2) within moving time window. Other explanations are in the text.
Figure 8. Graphs of changes in the components of the influence matrix between sequences of seismic events with a magnitude of at least 5.5 and the sequence of moments in time of 349 of the largest local maxima of the average amplitudes at the second level of EEMD decomposition. The estimates were made in a sliding time window of length 2 years with a shift of 0.01 year for a relaxation time τ of the model (24, 25) of 0.1 year. The graphs (a1c1) refer to the components of the influence matrix (35), which refers to the intensity fractions of the sequence of the largest local amplitude maxima, while the graphs (a2c2) refer to the intensity fractions of the sequence of seismic events. Plots (d1) and (d2) present the numbers of local maxima of amplitudes (d1) and the number of seismic events (d2) within moving time window. Other explanations are in the text.
Entropy 26 00710 g008
Figure 9. Maximum values of the elements of the influence matrices in a sequence of 100 time windows of length from 1 to 3 years, taken with an offset of 0.01 year: (a) for the “direct” influence of the time points of seismic events on the positions of the largest local maxima of average amplitudes on the second EEMD level of decomposition; (b) for the “reverse” influence of the positions of amplitude maxima on earthquakes. The relaxation time of the model is 0.1 year.
Figure 9. Maximum values of the elements of the influence matrices in a sequence of 100 time windows of length from 1 to 3 years, taken with an offset of 0.01 year: (a) for the “direct” influence of the time points of seismic events on the positions of the largest local maxima of average amplitudes on the second EEMD level of decomposition; (b) for the “reverse” influence of the positions of amplitude maxima on earthquakes. The relaxation time of the model is 0.1 year.
Entropy 26 00710 g009
Table 1. Number of Nsta stations in each Clust# cluster.
Table 1. Number of Nsta stations in each Clust# cluster.
Clust#.123456789101112131415
Nsta575654836961787791765795488857
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Lyubushin, A.; Rodionov, E. Prognostic Properties of Instantaneous Amplitudes Maxima of Earth Surface Tremor. Entropy 2024, 26, 710. https://doi.org/10.3390/e26080710

AMA Style

Lyubushin A, Rodionov E. Prognostic Properties of Instantaneous Amplitudes Maxima of Earth Surface Tremor. Entropy. 2024; 26(8):710. https://doi.org/10.3390/e26080710

Chicago/Turabian Style

Lyubushin, Alexey, and Eugeny Rodionov. 2024. "Prognostic Properties of Instantaneous Amplitudes Maxima of Earth Surface Tremor" Entropy 26, no. 8: 710. https://doi.org/10.3390/e26080710

APA Style

Lyubushin, A., & Rodionov, E. (2024). Prognostic Properties of Instantaneous Amplitudes Maxima of Earth Surface Tremor. Entropy, 26(8), 710. https://doi.org/10.3390/e26080710

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop