Prognostic Properties of Instantaneous Amplitudes Maxima of Earth Surface Tremor

1. Introduction

This article presents the further development of methods for analyzing ground surface tremor proposed in [1,2,3]. The coherence of 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] and biological [14,15] 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 structure of high-frequency noise and low-frequency seasonal component of GPS time series in connection with the task of estimating the rates of displacement of tectonic plates was studied in [16,17,18,19]. Seasonal effects in GPS signal series associated with hydrological load and tectonic movements of the earth's crust were studied in [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].

The fine structure of 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 [29,30,31,32]. In [33,34,35], the maximum likelihood method was used to estimate the shape parameters of power spectra and noise amplitude in different regions of the world. The noise structure and error estimates are discussed in [36,37]. In [38,39] analyzed the errors in estimating the displacement rate depending on the shape of the spectrum and the intensity of the noise component of the time series. Parametric models of GPS time series for the analysis of tectonically active areas, including the analysis of phase correlations, were used in [40,41].

Methods of multivariate statistics (principal components, empirical orthogonal functions and singular spectrum estimation) in [42,43] 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 [44]. Non-stationary effects were studied in [45,46,47] to assess the mutual movements of crustal blocks and the stability of station positions.

2. Data

Data on ground displacements measured using GPS are taken from the Nevada Geodetic Laboratory website [48] at:

http://geodesy.unr.edu/gps_timeseries/tenv3/IGS14/

The positions of 1047 GPS stations for the rectangular area ${30}^{\circ }\le \mathrm{Lat}\le {46}^{\circ }$, ${128}^{\circ }\le \mathrm{Lon}\le {146}^{\circ }$, which have daily time series from the beginning of 2009 to the end of 2023 (15 years), for which the total number of gaps is less than 365 samples and the longest gap is less than 182 samples, are presented in Figure 1(a). The vertical components of ground displacement are investigated. Gaps in the GPS time series are filled 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 clusters. To divide the network of stations into clusters, 15 reference points were selected (Figure 1). The number 15 was chosen as the optimal number of station clusters, which splits their “cloud” using the k-means method. Let's split the set of station position vectors $\overrightarrow{\zeta }$ into a given sample number $q$ of clusters using the popular k-means clustering method [49]. Let us denote by $C_r,\hspace{0.17em}r=1,\mathrm{...},q$ clusters, let ${\overrightarrow{z}}_r={\displaystyle {\sum }_{\overrightarrow{\zeta }\in \hspace{0.17em}{C}_r}\overrightarrow{\zeta }}/n_r$ be the vector of the center of mass of the cluster $C_r$, $n_r$ be the number of vectors in the cluster, ${\displaystyle {\sum }_{r=1}^{\hspace{0.17em}q}{n}_r}=N$. Vector $\overrightarrow{\zeta }\in C_r$ if the distance $|\overrightarrow{\zeta }-{\overrightarrow{z}}_r|$ is the minimum among the positions of all cluster centers. The k-means method minimizes the sum of squared distances:

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

(1)

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

$$PFS(q)={\sigma }_1^2(q)/{\sigma }_0^2(q)\to \underset{2\hspace{0.17em}\hspace{0.17em}\le \hspace{0.17em}\hspace{0.17em}q\hspace{0.17em}\hspace{0.17em}\le \hspace{0.17em}\hspace{0.17em}50}{\max }$$

(2)

where

$${\sigma }_0^2(q)={\displaystyle \frac{\Phi (q)}{N-q}},\hspace{0.17em}{\sigma }_1^2(q)={\displaystyle \sum _{r=1}^q\frac{n_r}{N}\cdot |{\overrightarrow{z}}_r-{\overrightarrow{z}}_0{|}^2},\hspace{0.17em}{\overrightarrow{z}}_0={\displaystyle \frac{1}{N}}{\displaystyle \sum _1^N\overrightarrow{\zeta }}$$

(3)

The plot in Figure 1(b) 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 1(a) 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 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. 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 [51]. 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),\hspace{0.17em}\hspace{0.17em}k=1,\mathrm{...},N\},\hspace{0.17em}j=1,\mathrm{...},p$. Let us choose the size of the sliding window $w$ and center the signals,

$$x_j(k)=y_j(k)-\hspace{0.17em}\overline{y_j(k)},\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},p, k=w,\mathrm{...},N,$$

(4)

$$\overline{y_j(k)}={\displaystyle \frac{1}{w}}{\displaystyle \sum _{i=1}^w{y}_j(k-w+i)},\hspace{0.17em}\hspace{0.17em}k=w,\mathrm{...},N$$

(5)

The next step is to calculate the sample estimate of the covariance $(p\times p)$-dimensional matrix in a sliding window:

$$r_{XX}^{(m,n)}={\displaystyle \frac{1}{w}}{\displaystyle \sum _{i=1}^w{x}_m(i-w+k)x_n(i-w+k)},\hspace{0.17em}\hspace{0.17em}k,j=1,\mathrm{...},p, k=w,\mathrm{...},N.$$

(6)

Let

${\varphi }^{(s)}={({\varphi }_1^{(s)},\mathrm{...},{\varphi }_p^{(s)})}^T$ be the eigenvector of this matrix corresponding to the maximum eigenvalue. Let's put

$${\xi }_s(k)\hspace{0.17em}=\hspace{0.17em}{\displaystyle \sum _{j=1}^p{\left({\varphi }_j^{(k)}\right)}^2{x}_j(k)}, k=w,\mathrm{...},N.$$

(7)

Having generalized formulas (4) – (7) with understandable changes to the case $k=1,\mathrm{...},w-1$, let us determine the weighted average in a sliding time window of length $w$ using the formula

$$\xi (k)=\{\begin{array}{l}{\xi }_{w-1}(k),\hspace{0.17em}k<w,\\ {\xi }_k(k),\hspace{0.17em}k\ge w.\end{array}\hspace{0.17em}$$

(8)

Thus, formulas (7-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 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 [52], which consists in eliminating outliers falling beyond the level $\mu \pm 4\sigma$ by cutting off the values ​​of the time series in a sliding time window ($\mu$ and $\sigma$ are sample estimates of the mathematical expectation and standard deviation for the current time window). The procedure is repeated iteratively until the values $\mu$ ​​and $\sigma$ 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 Modes 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)={\displaystyle \sum _{j=1}^n{h}_j(k)+r_n(k)}$$

(9)

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,\hspace{0.17em}m=0,1,\mathrm{...},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_{\hspace{0.17em}\hspace{0.17em}j}^{(m+1)}(k)=h_{\hspace{0.17em}\hspace{0.17em}j}^{(m)}(k)-z_{\hspace{0.17em}j}^{(m)}(k)$$

(10)

Here $z_{\hspace{0.17em}j}^{(m)}(k)=\left(p_{\hspace{0.17em}j}^{(m)}(k)+q_{\hspace{0.17em}j}^{(m)}(k)\right)/2$, where $p_{\hspace{0.17em}j}^{(m)}(k)$ and $q_{\hspace{0.17em}j}^{(m)}(k)$ are both the upper and lower envelopes for the signal, which are constructed using spline interpolation (usually a 3rd order spline) over all local maxima and minimums of the signal $h_{\hspace{0.17em}\hspace{0.17em}j}^{(m)}(k)$.

Iterations (10) are initialized with step zero for the first level ($j=1$) of the expansion $h_{\hspace{0.17em}\hspace{0.17em}1}^{(0)}(k)=x(k)$. Next, the upper $p_{\hspace{0.17em}1}^{(0)}(k)$ and lower $q_{\hspace{0.17em}1}^{(0)}(k)$ envelopes and are found, the middle line $z_{\hspace{0.17em}1}^{(0)}(k)$ is calculated and found $h_{\hspace{0.17em}\hspace{0.17em}1}^{(2)}(k)$ using formula (10). For $h_{\hspace{0.17em}\hspace{0.17em}1}^{(1)}(k)$, the upper $p_{\hspace{0.17em}1}^{(1)}(k)$ and lower $q_{\hspace{0.17em}1}^{(1)}(k)$ envelopes are determined, and the middle line $z_{\hspace{0.17em}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:

$${\displaystyle {\sum }_k{(h_{\hspace{0.17em}\hspace{0.17em}j}^{(m+1)}(k)-h_{\hspace{0.17em}\hspace{0.17em}j}^{(m)}(k))}^2}/{\displaystyle {\sum }_k{(h_{\hspace{0.17em}\hspace{0.17em}j}^{(m)}(k))}^2}\hspace{0.17em}\le \hspace{0.17em}\hspace{0.17em}\delta$$

(11)

where $\delta$ 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)$$

(12)

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),\hspace{0.17em}{h}_2(k),\mathrm{...},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 the form (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 of the analyzed sample back and forth for a sufficiently long period of time. It was this one that was used in this work.

5. Ensemble Empirical Modes 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, [9] proposed the ensemble empirical mode decomposition (EEMD) method. It 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:

• Adding 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.
• Obtaining the ensemble average for the corresponding empirical modes.

Thus, numerous “artificial” observations are simulated:

$$x^{(i)}(k)=x(k)+{\epsilon }_i(k)$$

(13)

where ${\epsilon }_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 small white noise, there immediately become a lot of local extrema). The operation of averaging over a sufficiently large number of independent implementations of white noise makes it possible to get rid of 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 [53]:

$$H_X(t)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{X(u)}{\pi (t-u)}}\hspace{0.17em}du=X(t)\ast \left({\displaystyle \frac{1}{\pi t}}\right)$$

(14)

where $f(t)\ast g(t)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}f(u)\cdot g(t-u)\hspace{0.17em}du}$ is the convolution of two functions. If $\tilde{f}(\omega )$ and $\tilde{g}(\omega )$ are the Fourier transforms of convolutional functions $\tilde{f}(\omega )={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}f(t)\cdot e^{-i\omega t}\hspace{0.17em}dt}$, 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/(\pi t)$ equals:

$${\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{e^{-i\omega t}}{\pi t}}\hspace{0.17em}dt=-i\cdot sign(\omega )=\{\begin{array}{c}-i,\hspace{1em}\omega >0\\ 0,\hspace{1em}\omega =0\\ i,\hspace{1em}\omega <0\end{array}$$

(15)

Thus, if $\tilde{X}(\omega )$ there is a Fourier transform of $X(t)$, then

$${\tilde{H}}_X(\omega )=-i\cdot sign(\omega )\cdot \tilde{X}(\omega )=\{\begin{array}{c}-i\cdot \tilde{X}(\omega ),\hspace{1em}\omega >0\\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\omega =0\\ i\cdot \tilde{X}(\omega ),\hspace{1em}\omega <0\end{array}$$

(16)

If you present $\tilde{X}(\omega )=|\tilde{X}(\omega )|e^{-i\phi (\omega )}$, then

$${\tilde{H}}_X(\omega )=\{\begin{array}{c}|\tilde{X}(\omega )|e^{-i(\phi (\omega )+\pi /2)},\hspace{1em}\omega >0\\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\omega =0\\ |\tilde{X}(\omega )|e^{-i(\phi (\omega )-\pi /2)},\hspace{1em}\omega <0\end{array}$$

(17)

In practice, it is more convenient to calculate the Hilbert transform using the concepts of an analytical signal:

$$Z_X(t)=X(t)+i\cdot H_X(t)=|Z_X(t)|\cdot e^{i\vartheta (t)}\equiv A_X(t)\cdot e^{i\vartheta (t)}$$

(18)

where $A_X(t)=\sqrt{X^2(t)+H_X^2(t)}$ are the amplitudes of the signal $X(t)$ envelope, and $\vartheta (t)$ is the instantaneous phase. The derivative $\nu (t)=d\vartheta (t)/dt$ is called instantaneous frequency. Fourier transform of the analytical signal ${\tilde{Z}}_X(\omega )=\tilde{X}(\omega )+i\cdot {\tilde{H}}_X(\omega )=\tilde{X}(\omega )\cdot (1+sign(\omega ))$ or:

$${\tilde{Z}}_X(\omega )=\{\begin{array}{c}2\tilde{X}(\omega ),\hspace{1em}\omega >0\\ \tilde{X}(0),\hspace{1em}\hspace{1em}\omega =0\\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\omega <0\end{array}$$

(19)

after which the Hilbert transform is equal to the imaginary part of the result of the inverse Fourier transform of ${\tilde{Z}}_X(\omega )$

$$H_X(t)=\mathrm{Im}\left(\frac{1}{2\pi }{ʃ}_{-\infty }^{+\infty }{\tilde{Z}}_X(\omega )e^{i\omega t}d\omega \right)$$

(20)

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.

For a discrete-time signal $X(t),\hspace{0.17em}\hspace{0.17em}t=0,\mathrm{...},(N-1)$, this transformation can be calculated using the discrete Fourier transform:

$$d_X^{(N)}\left({\omega }_k\right)=\sum _{t=0}^{N-1}X(t)·\exp \left(-i{\omega }_{k}t\right),\hspace{1em}{\omega }_k=\frac{2\pi }{N}(k-1),\hspace{1em}k=0,1,\dots ,(N-1)$$

(21)

after which the 2nd part of the Fourier coefficients (corresponding to negative frequencies) should be reset to zero: $h_{\hspace{0.17em}X}^{(N)}({\omega }_k)=0,\hspace{1em}k=N/2+1,\mathrm{...},(N-1)$ while the 1st part should be doubled: $h_{\hspace{0.17em}X}^{(N)}({\omega }_k)=\hspace{0.17em}\hspace{0.17em}2\cdot d_X^{(N)}({\omega }_k),\hspace{1em}k=1,\mathrm{...},N/2$. The Hilbert transform is then calculated as the imaginary part of the inverse discrete Fourier transform:

$$H_X(t)=\mathrm{Im}\hspace{0.17em}\left({\displaystyle \sum _{k=0}^{N-1}{h}_{\hspace{0.17em}X}^{(N)}({\omega }_k)\cdot \exp (i{\omega }_{k}t)/N}\right),\hspace{1em}t=0,1,\mathrm{...},(N-1)$$

(22)

Thus, after determining the instantaneous amplitudes and frequencies of the EEMD, the Hilbert-Huang decomposition can be represented as follows:

$$x(t)={\displaystyle \sum _{j=1}^n{h}_j(t)+r_n(t)},\hspace{1em}{h}_j(t)=\mathrm{Re}\hspace{0.17em}\left\{A_{\hspace{0.17em}j}(t)\cdot \exp \left(i\cdot {\displaystyle \int {\nu }_j(s)\hspace{0.17em}ds}\right)\right\}$$

(23)

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.

7. Influence Matrix

The next stage of the analysis is to study the connection between two sequences of random events. To solve this problem, a parametric intensity model is used. In [54,55,56] 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_{\hspace{0.17em}j}^{(\alpha )},\hspace{0.17em}j=1,\mathrm{...},N_{\alpha };\hspace{0.17em}\alpha =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's present their intensities in the form:

$$\hspace{0.17em}{\lambda }^{(\alpha )}(t)=b_{\hspace{0.17em}0}^{(\alpha )}\hspace{0.17em}+\hspace{0.17em}{\displaystyle \sum _{\beta =1}^2{b}_{\beta }^{(\alpha )}\cdot g^{(\beta )}(t)}$$

(24)

where $b_{\hspace{0.17em}0}^{(\alpha )}\ge 0,b_{\beta }^{(\alpha )}\ge 0$ are parameters, $g^{(\beta )}(t)\hspace{0.17em}$ is the influence function of flow events $t_j^{(\beta )}$ with number $\beta$:

$$g^{(\beta )}(t)={\displaystyle {\sum }_{t_j^{(\beta )}<t}\exp (-(t-t_j^{(\beta )})/\tau \hspace{0.17em})}$$

(25)

According to formula (25), the weight of an event with number $j$ becomes non-zero for times $t\hspace{0.17em}>t_j^{(\beta )}$ and decays with a characteristic time $\tau$. The parameter $b_{\beta }^{(\alpha )}$ determines the degree of influence of the stream $\beta$ on the stream $\alpha$. The parameter $b_{\alpha }^{(\alpha )}$ determines the degree to which the flow $\alpha$ influences itself (self-excitation), and the parameter $b_{\hspace{0.17em}0}^{(\alpha )}$ reflects the purely random (Poisson) intensity component. Let us fix the parameter $\tau$ and consider the problem of determining the parameters $b_{\hspace{0.17em}0}^{(\alpha )},b_{\beta }^{(\alpha )}$.

The log-likelihood function for a non-stationary Poisson process is equal over the time interval $[0,T]$ [57]:

$$\ln (L_{\alpha })={\displaystyle \sum _{j=1}^{N_{\alpha }}\ln ({\lambda }^{(\alpha )}(t_j^{(\alpha )}))-{\displaystyle \underset{0}{\overset{T}{\int }}{\lambda }^{(\alpha )}(s)ds}},\hspace{1em}\alpha =1,2$$

(26)

It is necessary to find the maximum of functions (26) with respect to the parameters $b_{\hspace{0.17em}0}^{(\alpha )},b_{\beta }^{(\alpha )}$. It is easy to obtain the following expression:

$$b_{\hspace{0.17em}0}^{(\alpha )}{\displaystyle \frac{\partial \ln (L_{\alpha })}{\partial b_{\hspace{0.17em}0}^{(\alpha )}}}+{\displaystyle \sum _{\beta =1}^2{b}_{\beta }^{(\alpha )}\frac{\partial \ln (L_{\alpha })}{\partial b_{\beta }^{(\alpha )}}}=N_{\alpha }-{\displaystyle \underset{0}{\overset{T}{\int }}{\lambda }^{(\alpha )}(s)ds}$$

(27)

Since the parameters $b_{\hspace{0.17em}0}^{(\alpha )},b_{\beta }^{(\alpha )}$ are non-negative, then each term on the left side of this formula is equal to zero at the maximum point 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:

$${\displaystyle \underset{0}{\overset{T}{\int }}{\lambda }^{(\alpha )}(s)ds=N_{\alpha }}$$

(28)

Let's substitute the expression $g^{(\beta )}(t)$ from (28) into (27) and divide by $T$. Then we get another form of formula (28):

$$b_{\hspace{0.17em}0}^{(\alpha )}+{\displaystyle \sum _{\beta =1}^m{b}_{\beta }^{(\alpha )}\cdot {\overline{g}}^{(\beta )}={\lambda }_{\hspace{0.17em}0}^{(\alpha )}\equiv }{N}_{\alpha }/T$$

(29)

where

$${\overline{g}}^{(\beta )}={\displaystyle {\int }_0^T{g}^{(\beta )}(s)ds}/T$$

(30)

- average value of the influence function. Substituting $b_{\hspace{0.17em}0}^{(\alpha )}$ from (29) into (26), we obtain the following maximum problem:

$${\Phi }^{(\alpha )}(b_1^{(\alpha )},b_2^{(\alpha )})={\displaystyle \sum _{j=1}^{N_{\alpha }}\ln ({\lambda }_{\hspace{0.17em}0}^{(\alpha )}+{\displaystyle \sum _{\beta =1}^2{b}_{\beta }^{(\alpha )}\cdot \Delta g^{(\beta )}(t_j^{(\alpha )}))\to \max }}$$

(31)

where $\Delta g^{(\beta )}(t)=g^{(\beta )}(t)-{\overline{g}}^{(\beta )}$, under restrictions:

$$b_1^{(\alpha )}\ge 0,b_2^{(\alpha )}\ge 0,{\displaystyle \sum _{\beta =1}^2{b}_{\beta }^{(\alpha )}{\overline{g}}^{(\beta )}}\le {\lambda }_{\hspace{0.17em}0}^{(\alpha )}$$

(32)

Function (31) is convex with a negative definite Hessian and, therefore, problem (31-32) has a unique solution. Having solved problem (31-32) numerically for a given $\tau$, you can enter the elements of the influence matrix ${\kappa }_{\beta }^{(\alpha )},\alpha =1,2;\hspace{0.17em}\beta =0,1,2$ according to the formulas:

$${\kappa }_0^{(\alpha )}={\displaystyle \frac{b_{\hspace{0.17em}0}^{(\alpha )}}{{\lambda }_{\hspace{0.17em}0}^{(\alpha )}}}\ge 0,\hspace{1em}{\kappa }_{\beta }^{(\alpha )}={\displaystyle \frac{b_{\beta }^{(\alpha )}\cdot {\overline{g}}^{(\beta )}}{{\lambda }_{\hspace{0.17em}0}^{(\alpha )}}}\ge 0$$

(33)

The value ${\kappa }_0^{(\alpha )}$ is part of the average intensity ${\lambda }_{\hspace{0.17em}0}^{(\alpha )}$ of the process with number $\alpha$, which is purely stochastic, part ${\kappa }_{\alpha }^{(\alpha )}$ is caused by the influence of self-excitation $\alpha \to \alpha$ and ${\kappa }_{\beta }^{(\alpha )},\hspace{0.17em}\beta \ne \alpha$ is due to external influence $\beta \to \alpha$. From formula (29) the normalization condition follows:

$${\kappa }_0^{(\alpha )}+{\displaystyle \sum _{\beta =1}^2{\kappa }_{\beta }^{(\alpha )}=1},\hspace{1em}\alpha =1,2$$

(34)

As a result, the influence matrix can be determined:

$$\left(\left|\begin{array}{l}{\kappa }_{\hspace{0.17em}\hspace{0.17em}0}^{(1)}\\ {\kappa }_{\hspace{0.17em}\hspace{0.17em}0}^{(2)}\end{array}\right|\left|\begin{array}{l}{\kappa }_{\hspace{0.17em}\hspace{0.17em}1}^{(1)}\hspace{0.17em}{\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}\\ {\kappa }_{\hspace{0.17em}1}^{(2)}\hspace{0.17em}{\kappa }_{\hspace{0.17em}2}^{(2)}\end{array}\right|\right)$$

(35)

The first column of matrix (35) is composed of Poisson fractions of average intensities. The diagonal elements of the right submatrix of size 2×2 consist of self-exciting elements of medium intensity, while the off-diagonal elements correspond to mutual excitation. The sums of the component rows of the influence matrix (34) are equal to 1. The influence matrices are estimated in a certain sliding time window of length $L$ with an offset $\Delta L$ and with a given value of the attenuation parameter $\tau$.

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. Such seismic events occurred in the vicinity of the Japanese islands during the period of time 2009-2023 it turned out to be 673 events - see Figure 5(a). However, the mega-earthquake of March 11, 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 sharply decreased, the number of earthquakes with a magnitude of at least 5.5 will decrease by almost 2 times and become equal 349 – see Figure 5(b).

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 2 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 2 levels are suitable for selecting 349 local maxima of amplitudes, since the number of local extrema already at the 3rd level is 242, that is, less than 349. From the first 2 levels, it was decided to choose the 2nd, as 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 2nd 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 the model (23, 24) of 0.1 year. The graphs (a1, b1, c1) refer to the components of the influence matrix (34), which refers to the intensity fractions of the sequence of the largest local amplitude maxima, while the graphs (a2, b2, c2) refer to the intensity fractions of the sequence of seismic events. Other explanations in the text.

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 one, 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)\Delta L$, $k=1,\mathrm{...},N_L$, $\Delta L=(L_{\max }-L_{\min })/(N_L-1)$. In our calculations, we took $L_{\min }$ 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 $\Delta t$. Let us denote by $t_j(L_k)$, $j=1,\mathrm{...},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 $\Delta t$. We used a time window offset $\Delta 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 $\tau$ of the model (26-27), corresponding to the mutual influence of the two processes being analyzed. We took a value $\tau$ 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),\hspace{0.17em}{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),\hspace{0.17em}{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's present each element $(t_j^{*}(L_k),\hspace{0.17em}{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,\mathrm{...},N_L$, $j=1,\mathrm{...},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: $\tau$, $L_{\min }$, $L_{\max }$, $N_L$, $\Delta 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 maximums of envelope amplitudes at the 2^nd level on the time instants of earthquakes is significantly greater (Figure 9(b)) than the “direct” influence of earthquakes on the occurrence of local maxima in the amplitudes of the envelopes (Figure 9(a)).

9. Conclusion

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.