Multifractal and Entropic Properties of Seismic Noise in the Japanese Islands

1. Introduction

The article is devoted to the analysis of seismic noise data recorded in the F-net network of stations on the Japanese islands for 29 years, 1997-2025. During this period, on March 11, 2011, a mega-earthquake with a magnitude of 9.1 occurred in Japan. Processes within the Earth’s crust, including the preparation of strong seismic events, are reflected in changes in the properties of seismic noise. The availability of the F-net seismic network with free access to data, provided by the National Research Institute for Earth Science and Disaster Resilience (NIED), makes it possible to test various hypotheses about how changes in the properties of seismic noise reflect the preparation of a strong seismic event. In works [1,2,3], it was shown that the processes of earthquake preparation are preceded by changes in the statistical structure of seismic noise, consisting in simplification of noise, namely, an increase in entropy and a loss of multifractality. This article is a continuation of works [4,5] on the study of the properties of seismic noise in Japan. The possibility of a catastrophic earthquake in the deep-sea Nankai Trough region south of Tokyo has long been a subject of discussion among Japanese seismologists [6,7]. The Tohoku mega-earthquake of March 11, 2011, revived interest in assessing the likelihood of a repeat mega-earthquake in this area [8,9]. A hypothesis was put forward in [10] that a new mega-earthquake could even exceed the Tohoku event in energy and its magnitude could reach 10. The emergence of new data makes it possible to assess the possibility of a new strong earthquake. This paper shows that significant changes in the spatiotemporal structure of seismic noise in Japan occurred in 2020-2025, which may indicate an acceleration of the preparation processes for a new mega-earthquake.

2. Initial Data

The analysis was performed using vertical seismic oscillation data with a sampling frequency of 1 Hz [11]. These data are available at the address for 78 seismic stations of the F-net network in Japan. The analysis was performed over the period from the beginning of 1997 to December 31, 2025. Seismic data with a sampling frequency of 1 Hz were normalized to a time step of 1 minute by calculating the average values in adjacent time intervals of 60 samples.

Figure 1. Purple dots represent the locations of F-net seismic stations on the Japanese islands. Gray circles are drawn with a radius of 250 km centered on each station. The union of these circles represents the network’s “area of influence,” which is being explored further.

Figure 1. Purple dots represent the locations of F-net seismic stations on the Japanese islands. Gray circles are drawn with a radius of 250 km centered on each station. The union of these circles represents the network’s “area of influence,” which is being explored further.

Figure 2 shows a daily graph of the number of working stations. A seismic station is considered operational for the current day if there are no gaps during that day.

3. Seismic Noise Statistics

The minimum entropy of a time series is determined by the formula $En=-{\displaystyle \sum _k{p}_k\cdot \log (}{p}_k)/\log (N)$, where $p_k=c_k^2/{\displaystyle {\sum }_j{c}_j^2}$, $c_k$ are the wavelet coefficients of the signal decomposition, and $N$ is the total number of coefficients $c_k$. Seventeen orthogonal Daubechies wavelets were used: 10 ordinary minimum support bases with a number of zeroed moments from 1 to 10 and 7 so-called Daubechies symlets [13] with a number of vanishing moments from 4 to 10. In each time window, the wavelet for which the value $En$ is minimal is selected.

The Donoho-Johnston (DJ) wavelet-based index $\gamma$ is the ratio of “large” wavelet coefficients in absolute value to their total number. By definition $0\le \gamma \le 1$. The threshold separating “large” wavelet coefficients is $\sigma \sqrt{2\cdot \ln N}$, where $\sigma =med\left\{\left|c_k^{(1)}\right|\hspace{0.33em},\hspace{0.33em}k=1,\mathrm{...},N/2\right\}/0.6745$ is a robust estimate of the standard deviation of the normal distribution, $c_k^{(1)}$ are the wavelet coefficients at the first detail level of decomposition, and $N/2$ is the number of such coefficients [13,14].

The singularity spectrum support width $\Delta \alpha$ is considered as a measure of the diversity of the stochastic behavior of the signal $u(t)$. It is defined as $\Delta \alpha ={\alpha }_{\max }-{\alpha }_{\min }$, where ${\alpha }_{\min }$ and ${\alpha }_{\max }$ are the minimum and maximum values of the Hölder-Lipschitz exponent [15], which governs the behavior of the signal in the vicinity of the time instant $t$: $|u(t+\frac{\delta }{2})-u(t-\frac{\delta }{2})|\hspace{0.33em}\sim \hspace{0.17em}|\delta {|}^{\hspace{0.17em}\alpha },\delta \to 0$. For a monofractal signal, the exponent $\alpha$ is the same for all time instants $t$. If this exponent differs, then the signal is multifractal [16].

After switching to a 1-minute sampling time step, the seismic records from each station were divided into 1-day-long time fragments (1440 samples), and the parameters $(En,\gamma ,\Delta \alpha )$ of daily seismic noise signals were calculated for each fragment. The methods for estimating the parameters $En$, $\gamma$ and $\Delta \alpha$ for seismic noise records in a sliding time window are described in detail in [4]. The DFA method [17] was used to calculate the values of the singularity spectrum carrier width. Detrending of seismic noise signals using an 8th-order polynomial was used before calculating the entropy and index in each daily time window. Thus, time series of values with a 1-day time step were obtained for each of the seismic stations. The graphs of these quantities are presented in Figure 4 (a1, b1, c1).

The entropy value $En$ used in this paper, defined through the coefficients of the orthogonal wavelet decomposition, has features in common with multiscale entropy used, for example, in biology [18,19]. In [20,21], entropy is used within the framework of the natural time approach for seismic data analysis. In [22,23], non-extensive Tsallis entropy is used for processing seismic noise data.

The use of multifractal analysis to study the behavior of various complex systems has a long history. Particular attention is paid to the effect of parameter $\Delta \alpha$ reduction (loss of multifractality) preceding changes in system properties. In medicine, a decrease in the $\Delta \alpha$ value of various parameters accompanies age-related changes [24,25,26]. In [27,28], multifractal analysis is used to analyze geoelectric signals and wind speed. Within the framework of the natural time approach, multifractal analysis is used to study both seismicity and other time sequences [29,30].

In the future, we will be interested in the positions of the time points of a certain number of the largest local maxima or the smallest local minima of daily seismic noise properties in comparison with the times of earthquakes with a magnitude of at least 6 (Figure 3). To eliminate the influence of low-frequency components of changes in statistical values on the determination of the times of local extrema, the time series of noise properties were subjected to the operation of removing low frequencies using Gaussian kernel smoothing. Let $u(t)$ be a time series with discrete time $t$. Gaussian kernel averaging of a time series $u(t)$ with radius (scale parameter) $h>0$ at time $t$, is calculated using the formula [31]:

$$\overline{u}(t|h)={\displaystyle \sum _{s}u(s)\cdot \exp (-{\left(\frac{t-s}{h}\right)}^2}/{\displaystyle \sum _s\exp (-{\left(\frac{t-s}{h}\right)}^2}$$

(1)

Next, the average values of the time series for an averaging radius $h$ of 2 days were subtracted from the time series of statistical changes, and the most significant local extremum points were found for the residuals. These operations are illustrated by the graphs in Figure 4 (a2, b2, c2).

Figure 3. Time sequence of 352 earthquakes with magnitude of at least 6 in the rectangular vicinity of the Japanese islands (Figure 1) for the time period 1997-2025. Data source [12].

Figure 3. Time sequence of 352 earthquakes with magnitude of at least 6 in the rectangular vicinity of the Japanese islands (Figure 1) for the time period 1997-2025. Data source [12].

Figure 4. Figures (a1), (b1), (c1) present the graphs of daily median values of the Donoho-Johnston index $\gamma$, the singularity spectrum support width $\Delta \alpha$, and the entropy $En$, respectively. The green lines represent the moving averages in a 57-day window. The gray lines in figures (a2), (b2), (c2) represent the result of removing local trends with a Gaussian smoothing window of radius 2 days for the curves presented in figures (a1), (b1), (c1). The blue dots in figures (a2) and (b2) represent the 352 smallest local minima of the values $\gamma$and $\Delta \alpha$ after removing local trends; the red dots in figure (c2) represent the 352 largest local maxima of the entropy values $En$ after removing local trends. The number 352 is the number of earthquakes with a magnitude of at least 6 in the rectangular area shown in Figure 1, 1997-2025 (Figure 3).

Figure 4. Figures (a1), (b1), (c1) present the graphs of daily median values of the Donoho-Johnston index $\gamma$, the singularity spectrum support width $\Delta \alpha$, and the entropy $En$, respectively. The green lines represent the moving averages in a 57-day window. The gray lines in figures (a2), (b2), (c2) represent the result of removing local trends with a Gaussian smoothing window of radius 2 days for the curves presented in figures (a1), (b1), (c1). The blue dots in figures (a2) and (b2) represent the 352 smallest local minima of the values $\gamma$and $\Delta \alpha$ after removing local trends; the red dots in figure (c2) represent the 352 largest local maxima of the entropy values $En$ after removing local trends. The number 352 is the number of earthquakes with a magnitude of at least 6 in the rectangular area shown in Figure 1, 1997-2025 (Figure 3).

Let’s consider how the time points of the most significant local extrema of the statistics $En$, $\gamma$ and $\Delta \alpha$ are distributed. To do this, we combine all the extremum points (352 smallest local minima $\gamma$, $\Delta \alpha$ and 352 largest local maxima of entropy $En$) into a single time sequence, sort it in ascending order, and calculate the empirical probabilities of time values falling within successive six-month intervals; that is, we calculate a probability histogram of time points. Since the total length of observations is 29 years, the number of such intervals is 58. Figure 5 shows a graph of this histogram. It shows that the last six-month interval has the highest empirical probability of containing local extremum points.

4. Lead Times Between Local Extremes of Seismic Noise Properties and Earthquakes

To quantitatively estimate the relative positions of time points of local extrema of seismic noise properties and earthquake times, a measure of the mutual advance of point processes, proposed in [32] for the multivariate case, is used. For the particular case of joint processing of two point processes, this method was used in [5,34,35,36,37,38] to analyze the relationships between seismic noise properties, fluctuations of the global magnetic field, GPS tremor of the Earth’s surface, properties of the solar proton flux, and anomalies of meteorological parameters with the seismic process.

The lead measure is calculated using a parametric model of the intensity of two interacting point processes. Let $t_{\hspace{0.17em}j}^{(\alpha )},\hspace{0.33em}j=1,\mathrm{...},N_{\alpha };\hspace{0.33em}\alpha =1,2$ represent the time instants of two event streams. Let their intensities be represented as:

$$\hspace{0.33em}{\lambda }^{(\alpha )}(t)=b_{\hspace{0.17em}0}^{(\alpha )}\hspace{0.17em}+b_{\hspace{0.17em}1}^{(\alpha )}{g}^{(1)}(t)+b_{\hspace{0.17em}2}^{(\alpha )}{g}^{(2)}(t)\hspace{0.17em}$$

(2)

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

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

(3)

According to formula (3), the weight of event number $j$ becomes nonzero for times $t\hspace{0.17em}>t_j^{(\beta )}$ and decays with characteristic time $\tau$. The parameter $b_{\beta }^{(\alpha )}$ determines the degree of influence of the event stream $\beta$ on the stream $\alpha$. The parameter $b_{\alpha }^{(\alpha )}$ determines the degree of influence of the stream $\alpha$ on itself (self-excitation), and the parameter $b_{\hspace{0.17em}0}^{(\alpha )}$ reflects the purely random (Poisson) component of intensity. 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 to over the time interval $[0,T]$ [39]:

$$\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}},\alpha =1,2$$

(4)

It is necessary to find the maximum of functions (4) with respect to parameters $b_{\hspace{0.17em}0}^{(\alpha )},b_{\beta }^{(\alpha )}$. It can be shown [32,37] that problem (4) is reduced to a maximization problem

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

(5)

where $\Delta g^{(\beta )}(t)=g^{(\beta )}(t)-{\overline{g}}^{(\beta )}$, ${\overline{g}}^{(\beta )}={\displaystyle {\int }_0^T{g}^{(\beta )}(s)ds}/T$, ${\lambda }_{\hspace{0.17em}0}^{(\alpha )}=N_{\alpha }/T$, under restrictions:

$$b_{\hspace{0.33em}1}^{(\alpha )}\ge 0,b_{\hspace{0.33em}2}^{(\alpha )}\ge 0,b_{\hspace{0.33em}1}^{(\alpha )}{\overline{g}}^{(1)}+b_{\hspace{0.33em}2}^{(\alpha )}{\overline{g}}^{(2)}\le {\lambda }_{\hspace{0.17em}0}^{(\alpha )}$$

(6)

Having solved the problem (5-6) numerically for a given $\tau$, we can introduce the elements of the influence matrix ${\kappa }_{\beta }^{(\alpha )},\alpha =1,2;\hspace{0.33em}\beta =0,1,2$ according to the formulas:

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

(7)

The quantity ${\kappa }_0^{(\alpha )}$ is the portion of the average intensity ${\lambda }_{\hspace{0.17em}0}^{(\alpha )}$ of the process with index $\alpha$, which is purely stochastic, the portion ${\kappa }_{\alpha }^{(\alpha )}$ caused by self-excitation $\alpha \to \alpha$, and the portion ${\kappa }_{\beta }^{(\alpha )},\hspace{0.33em}\beta \ne \alpha$ caused by external influence $\beta \to \alpha$. As a result, the influence matrix can be determined:

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

(8)

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

Next, optimization is performed with respect to the parameters $L$, $\tau$ with the aim of maximizing the average value of the difference $\Delta \kappa ={\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}-{\kappa }_{\hspace{0.33em}1}^{(2)}\to \underset{\tau ,\hspace{0.33em}L}{\max }$, where ${\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}$ is the portion of the intensity of the earthquake sequence (the first process) caused by the “influence” of the sequence of local extrema (the second process), and ${\kappa }_{\hspace{0.33em}1}^{(2)}$ is the portion of the intensity of the sequence of local extrema caused by the leading influence of seismic events. We will agree to call ${\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}$ the “direct” lead, and the ${\kappa }_{\hspace{0.33em}1}^{(2)}$ “reverse” lead. The greater the difference $\Delta \kappa$, the greater the effect of the leading of the local extrema by the time instants relative to the time instants of the earthquakes. The problem of finding the maximum is solved within an a priori range of parameter variations: $L_{\min }\le L\le L_{\max }$, ${\tau }_{\min }\le \tau \le {\tau }_{\max }$. The values $L_{\min }=0.5$, $L_{\max }=5$, ${\tau }_{\min }=0.05$, ${\tau }_{\max }=5$, were used in the calculations; all time units are in years. The shift $\Delta L$ of the time windows was taken equal to 0.05 years. The Table 1 presents parameters of the model (2-3) which were obtained after solving the problem $\Delta \kappa \to \underset{\tau ,\hspace{0.33em}L}{\max }$.

times.

Figure 6 shows graphs of changes in the values of ${\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}$ and ${\kappa }_{\hspace{0.33em}1}^{(2)}$ for points of local extrema of the analyzed properties of seismic noise with the optimal choice of parameters.

Let us denote by ${\overline{\kappa }}_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}$ the result of averaging the measures ${\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}$ of “forward” advance for various properties of seismic noise, presented by the graphs in Figure 6 ((a1), (b1), (c1)). This averaging is realized by calculating the Gaussian means of the union of these dependencies using formula (1) for an averaging radius $h$ equal to 10 days and for a uniform sequence of time points with a step of 10 days. Next, we calculate the same result ${\overline{\kappa }}_{\hspace{0.33em}1}^{(2)}$ of averaging the measures ${\kappa }_{\hspace{0.33em}1}^{(2)}$ of “reverse” advance, presented by the graphs in Figure 6 ((a2), (b2), (c2)). Figure 7 shows their difference $\Delta \overline{\kappa }={\overline{\kappa }}_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}-{\overline{\kappa }}_{\hspace{0.33em}1}^{(2)}$ as a function of the position of the right end of successive time intervals of 10 days.

The value $\Delta \overline{\kappa }$ is an integrated measure of the lead time of local extrema in seismic noise properties relative to the times of earthquakes with magnitudes of at least 6. The vertical blue line in Figure 7 represents the time of the Tohoku mega-earthquake of March 11, 2011. The graph in Figure 7 shows that this measure is mostly positive; that is, for most times, the points of local extrema in seismic noise properties lead the times of earthquakes. An exception is a time interval of just under 4 years, adjacent to the time of the Tohoku event on the right and extending to the 2015 time stamp, when the value $\Delta \overline{\kappa }$ is negative. The graph in Figure 7 also shows other time segments when $\Delta \overline{\kappa }<0$, but they are not as long as after the Tohoku earthquake. The time intervals $\Delta \overline{\kappa }<0$ correspond to a situation when, on average, the attainment of extreme values of properties $\gamma$, $\Delta \alpha$ and $En$ is a consequence of earthquakes. If $\Delta \overline{\kappa }>0$, then this means that local extremes of noise properties foreshadow seismic events.

For two long time periods before the Tohoku event and after the 2015 timestamp, quasi-periodic fluctuations in the measure $\Delta \overline{\kappa }$ values are noticeable. We estimate their periods by fitting harmonic oscillations (shown as green lines in Figure 7) using the least-squares method with unknown periods, which we find by minimizing the residual variance. It turns out that these periods are quite close to each other, equal to 3.5 years before the Tohoku event and 3.2 years after 2015.

5. The First Principal Components of the Amplitudes of the EEMD Envelope Decompositions of Seismic Noise Properties

Let us consider the issue of identifying the general fluctuations of the median values of the noise properties $\gamma$, $\Delta \alpha$ and $En$, the graphs of which are presented in Figure 4 ((a1), (b1), (c1)). To solve this problem, we will use the Hilbert-Huang method [40,41] of decomposing signals into sequences of intrinsic mode functions (IMF), also known as empirical mode decomposition (EMD). The use of EMD decomposition has a wide range of applications. In [29], this method was used to analyze seismicity on various time scales in combination with multifractal analysis. In [35,37], EEMD decomposition was used to analyze the prognostic properties of GPS tremor of the earth’s surface and anomalies of meteorological time series.

The decomposition of an arbitrary signal into empirical oscillation modes is given by the formula:

$$x(t)={\displaystyle \sum _{j=1}^n{h}_j(t)+r_n(t)}$$

(9)

where $h_j(t)$ is the j-th empirical mode, $r_n(t)$ is the remainder, and n is the number of empirical modes.

The algorithm for decomposing into a sequence of empirical modes is iterative for each level $j$. We denote by $m,\hspace{0.33em}m=0,1,\mathrm{...},M_j$ the iteration index, where $M_j$ is the maximum number of iterations for level $j$. The iterations are described by the formula

$$h_{\hspace{0.33em}\hspace{0.33em}j}^{(m+1)}(t)=h_{\hspace{0.33em}\hspace{0.33em}j}^{(m)}(t)-z_{\hspace{0.33em}j}^{(m)}(t)$$

(10)

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

Iterations (10) are initialized with a zero step for the first level ($j=1$) of the decomposition $h_{\hspace{0.33em}\hspace{0.33em}1}^{(0)}(t)=x(t)$. Next, the upper $p_{\hspace{0.33em}1}^{(0)}(t)$ and lower $q_{\hspace{0.33em}1}^{(0)}(t)$ envelopes and are found, the mean line $z_{\hspace{0.33em}1}^{(0)}(t)$ is calculated, and $h_{\hspace{0.33em}\hspace{0.33em}1}^{(1)}(t)$ is found using formula (10). For $h_{\hspace{0.33em}\hspace{0.33em}1}^{(1)}(t)$, the upper $p_{\hspace{0.33em}1}^{(1)}(t)$ and lower $q_{\hspace{0.33em}1}^{(1)}(t)$ envelopes and are determined, and the mean line $z_{\hspace{0.33em}1}^{(1)}(t)$ are found, and so on, up to the last iteration index $M_1$, after which the first empirical mode $h_1(t)$ is considered to have been found.

The condition for stopping the iterations is usually chosen in the form of the inequality:

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

(11)

where $\delta$ is some small number, such as 0.01. Once the mode $h_j(t)$ is found, an iterative process is started to determine the empirical mode $h_{j+1}(t)$ of the next level. This process is initialized by the formula for the initial iteration index $m=0$:

$$h_{j+1}^{(0)}(t)=x(t)-h_j(t)$$

(12)

According to formula (12), the high-frequency component is subtracted from the original signal, and the resulting lower-frequency signal is considered a new signal for subsequent decomposition. The construction of empirical oscillation modes continues until the number of local extrema becomes too small to be used to construct envelopes.

As the empirical mode level number $j$ increases, the signals $h_j(t)$ become increasingly low-frequency and tend toward an unchanging form. The sequence $h_1(t),\hspace{0.33em}{h}_2(t),\mathrm{...},h_n(t)$ is constructed such that its sum yields an approximation of the original signal, which can be represented as (9).

Empirical oscillation modes, known as Intrinsic Mode Functions (IMFs), are orthogonal to each other, thus forming the empirical basis for decomposing the original signal. From here on, these decomposition levels will be referred to as IMF levels.

One of the drawbacks of the EMD method is the problem of mode mixing, which occurs when a single empirical mode includes signals of different scales or when signals of the same scale are distributed across different empirical modes. For example, if a signal exhibits “intermittency,” meaning that short-lived sections of higher-frequency behavior appear against a smooth background, then EMD decomposition results in a mixing of modes with different frequencies, as the relatively rare local extrema of the smooth behavior are interspersed with significantly more frequent local extrema of the high-frequency component.

To combat this effect, the ensemble empirical mode decomposition (EEMD) method was proposed in [41]. It is a regularization of the EMD method, in which finite-amplitude white noise is added to the original data. This allows the true empirical modes to be determined as the average value over an ensemble of trials, each of which is the sum of the signal and white noise.

The EEMD algorithm includes the following steps:

• Adding a white noise realization to the original data.
• Decomposing the data with the added white noise into empirical modes.
• Repeating steps 1 and 2 a sufficiently large number of times with independent white noise realizations.
• Obtaining the ensemble mean for the corresponding empirical modes.

Thus, numerous artificial observations are simulated in real time series:

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

(13)

where ${\epsilon }_i(t)$ is the i-th realization of white noise. In our calculations, we used 1000 realizations of independent white noise with a standard deviation equal to 0.1 of the standard deviation of the original signals.

After obtaining EEMD realizations of daily time series of noise properties $\gamma$, $\Delta \alpha$ and $En$ for a sequence of IMF levels, we calculate the amplitudes of their envelopes using the Hilbert transform, which is efficiently implemented using the fast Fourier transform [42]. For each IMF level, we have 3 time functions of the envelope amplitudes corresponding to different properties of seismic noise. To isolate their common properties at each level, we calculate their first principal component after centering and normalization to unit variance [43,44]. The graphs of the first principal components of the envelope amplitudes are presented in Figure 8 for the first 8 IMF levels of the EEMD decomposition. Each IMF level approximately corresponds to a frequency band with boundary periods from $2^j\Delta t$ to $2^{(j+1)}\Delta t$, where $j$ is the IMF level number, $\Delta t$ is the time sampling step (in our case, 1 day).

Among the graphs in Figure 8, the low-frequency IMF levels 5-8 are particularly striking. At IMF levels 5-7, peaks in the first principal component centered on the time of the Tohoku mega-earthquake of March 11, 2011, are visible, indicating a post-seismic response. However, the IMF level 8 graph (periods from 256 to 512 days) displays a significant peak for the time interval from mid-2024 to late 2025, which may indicate preparation for the next strong earthquake.

6. Probability Densities of Extreme Values of Seismic Noise Properties

Let us study the variability of the spatial distribution of extreme values of seismic noise properties. For this purpose, we consider a regular grid of 30×30 nodes, covering an area with latitude from 28°N to 46°N and longitude from 128°E to 146°E (Figure 1). Let $W$ be any value of $\gamma$, $\Delta \alpha$ or $En$. For each grid node $(i,j)$ and for each day with number $t$, we find the 5 nearest working seismic stations, which yields 5 values of $W$. Let us denote by $W_{ij}^{(t)}$ the median value of these 5 properties at a node $(i,j)$ on a day with number $t$. The values $W_{ij}^{(t)}$ in the set of nodes of the regular grid can be considered as an elementary daily map. For each daily elementary map with a discrete time index $t$, we find the coordinates of the nodes ${\zeta }_{\hspace{0.17em}mn}^{(t)}=(x_m^{(t)},y_n^{(t)})$ at which a given number $n_m$ of extreme values of $W$ is reached relative to all other nodes of the regular grid. Minima are found for $W=\Delta \alpha$ or $W=\gamma$, and maxima are found for $W=En$. Further, $n_m=10$ extreme values will be used. The set of two-dimensional vectors ${\zeta }_{\hspace{0.17em}mn}^{(t)}$ considered within the time interval $t\in [t_0,t_1]$ forms a random set. We estimate the two-dimensional probability distribution function of the vectors ${\zeta }_{\hspace{0.17em}mn}^{(t)}$ for each node of the regular grid. For each node $(i,j)$, the probability density function of the distribution of extreme values of the property $W$ is calculated according to the Parzen-Rosenblatt estimate with the Gaussian kernel function [44]:

$$p({\zeta }_{ij}|t_0,t_1)={\displaystyle \frac{1}{2\pi n_m{H}^2{N}_{t_0,t_1}}}{\displaystyle \sum _{t=t_0}^{t_1}{\displaystyle \sum _{mn}\exp \left(-\frac{|{\zeta }_{ij}-{\zeta }_{\hspace{0.17em}mn}^{(t)}{|}^2}{2H^2}\right)}}$$

(14)

Here $H$ is the kernel averaging radius, $t_0,t_1$ are integer indices that number the daily elementary maps, and $N_{t_0,t_1}=(t_1-t_0+1)$ is the number of daily maps in the time interval under consideration. A smoothing bandwidth of $H=1^{\circ }$ was used.

Let us calculate the 2-dimensional distribution densities (14) of extreme values in successive time windows of 10 days in length ($N_{t_0,t_1}=10$). Since we have three such distributions, we will calculate their weighted average by taking as weights the squares of the components of the eigenvector of the correlation matrix of probability densities corresponding to the maximum eigenvalue of the matrix [43]. By construction, the sum of the squares of the components is equal to 1. Figure 9 shows graphs of the weighting coefficients depending on the position of the right end of time windows of 10 days in length.

To quantitatively describe the sequence of values of the weighted probability density function (PDF) for the distribution of extreme seismic noise properties in 10-day time windows, we consider the entropy of these densities. Let $M$ be the number of grid nodes within the region formed by the union of circles with a radius of 250 km, centered at each F-net station (Figure 1). We calculate the histogram of the weighted probability density functions for each 10-day interval $j$ as a sequence of values:

$$h_{\hspace{0.33em}k}^{(j)}=m_{\hspace{0.33em}k}^{(j)}/M,k=1,\mathrm{...},M_b$$

(15)

where $M_b$ is the number of intervals dividing the range of density variation within a given 10-day window into intervals of equal length. According to recommendations [31,44], we choose the value $M_b=\sqrt{M}$. Obviously ${\displaystyle {\sum }_k{h}_{\hspace{0.33em}k}^{(j)}=1}$, the values of the histogram (15) also have the properties of a discrete probability distribution. We define the normalized entropy of this distribution by the formula:

$$Q_j=-{\displaystyle \sum _{k=1}^{M_b}{h}_{\hspace{0.33em}k}^{(j)}\log (h_{\hspace{0.33em}k}^{(j)})/\log (M_b),0\le Q_j\le 1}$$

(16)

In Figure 10, the gray line represents the graph of entropy change (16) depending on the position of the right end of 10-day time windows. The red line represents the result of smoothing with a Gaussian kernel (formula (1)) with a smoothing kernel width $h$ of 400 days.

The blue line represents a stepwise approximation using so-called long WTMM (Wavelet Transform Modulus Maxima) chains of the signal skeleton. Scale-dependent stepwise approximation allows one to formally determine the moments in time of significant changes in the mean entropy values. This approximation is constructed using continuous wavelet transforms with a kernel in the form of derivatives of a Gaussian function, which allows one to formalize the selection of moments in time at which a significant, scale-dependent change in the mean value of a noisy signal occurs [13,45].

To describe the construction of the stepwise WTMM approximation, we will use continuous-time formulas, which are easily adapted to the discrete-time case but are more compact. The smoothed entropy value is calculated using the formula:

$$Q_0(t,a)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}Q(t+av)\cdot {\psi }_0(v)\hspace{0.17em}dv}/{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\psi }_0(v)\hspace{0.17em}dv},\hspace{0.33em}{\psi }_0(v)=\exp (-v^2)$$

(17)

where $a\hspace{0.33em}>0$ is the smoothing time scale. For the $n$-th derivative of the smoothed signal, divided by $n!$ (the Taylor coefficient), the following expression holds:

$$Q_n(t,a)\equiv {\displaystyle \frac{1}{n!}}{\displaystyle \frac{d^n{Q}_0(t,a)}{d{t}^n}}={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}Q(t+av){\psi }_n}(v)\hspace{0.17em}dv/a^n{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{v}^n{\psi }_n(v)\hspace{0.17em}dv}$$

(18)

where ${\psi }_n(\nu )={(-1)}^n\cdot {\displaystyle \frac{d^n{\psi }_0(\nu )}{d{\nu }^n}}\equiv {(-1)}^n\cdot {\psi }_{\hspace{0.33em}\hspace{0.17em}0}^{\hspace{0.17em}(n)}(\nu )$. A two-dimensional WTMM point $(t,a)$ for $n\ge 0$ is defined as a point for which $|Q_n(t,a)|$ exhibits a local maximum in time $t$ for a given time scale $a$. For $n=0$, WTMM points are defined as points of local extrema (maxima or minima) of the smoothed signal $Q_0(t,a)$, which can be combined into chains. The set of all chains forms the WTMM skeleton of the signal. If the kernel ${\psi }_0(\nu )$ is Gaussian, then a given chain of the WTMM skeleton does not terminate (does not “hang in the air”) when decreasing the scale and touches the time axis [46]. The WTMM points for the 1st-order derivative $Q_1(t,a)$ indicate the time points of maximum trend (positive or negative) of the smoothed signal $Q_0(t,a)$ for a given scale $a$. This can be used to formally identify points of large changes in the mean. These times correspond to the roots of long WTMM chains reaching given values of the scale $a^{*}$. We define a stepwise WTMM approximation for a signal $Q(t)$ as a function that is equal to a sequence of constant values ${\overline{Q}}_k$ in successive intervals $[{\tau }_k(a^{*}),{\tau }_{k+1}(a^{*})]$. Here ${\tau }_k(a^{*})$ denote the beginnings of WTMM chains for $Q_1(t,a)$, which exceed a threshold time scale $a^{*}$ and stepwise approximation is equal to the sequence of average values ${\overline{Q}}_k$ within the time intervals $[{\tau }_k(a^{*}),{\tau }_{k+1}(a^{*})]$.

In Figure 10, the blue line represents the stepwise WTMM approximation for the parameter $a^{*}=450$ days. This threshold time scale ensures that the overall time interval is divided into 12 fragments with different average entropy values, which are then interpreted as stages in the development of seismic noise field features in the Japanese islands. A comparison of the red and blue lines in Figure 10 shows that both averaging methods exhibit roughly the same behavior in the most general trends. In particular, the final time interval from the end of 2020 is characterized by elevated average entropy values.

7. A Sequence of Distribution Maps of Average-Weighted Probability Densities of Extreme Values of Noise Properties

Figure 11 shows a sequence of distributions of average-weighted probability densities of extreme values of seismic noise properties in the region of the Japanese islands in the union of circles with a radius of 250 km, constructed around each seismic station for 12 time intervals allocated by the stepwise WTMM approximation of entropy (16).

Within the framework of the considered hypothesis, the areas of concentration of the weighted probability density are interpreted as “spots of increased seismic hazard.” The time intervals identified by the stepwise approximation of entropy (16) can be represented as a sequence of 12 stages of precursor development based on the properties of seismic noise. The first six distribution maps (α01-α06) correspond to the time intervals preceding the Tohoku event, the epicenter of which is shown by a black asterisk. In the first two maps (α01, α02), it is noticeable how a single area of concentration of the probability density at α01 splits into two parts, one of which at α02 closely approaches the epicenter of the Tohoku event. Stage α03 lasts for approximately six months, and during this short period of time, for which the average entropy value (16) declines, a region of probability density concentration formed in the region of 30°N≤Lat≤34°N, 136°E≤Lon≤140°E, which remains stable for all later stages. A clear precursor to the future Tohoku mega-earthquake of March 11, 2011, formed during stages (α04-α06). Note the significant duration of the time interval corresponding to stage α06.

Five stages (α07-α11) demonstrate a relatively calm evolution of “seismic hazard spots” with periodic moderate concentration of probability density in the Tohoku aftershock region, the bursts of which correspond to an increase in the average entropy values (16) in Figure 10. Finally, the last stage α12 is characterized by a significant rearrangement of the configuration of seismic hazard spots. Concentration of probability density is observed along the entire length of the Pacific coast of Japan (subduction region), while the area of the highest density values remains low-lying in the region 30°N≤Lat≤34°N, 136°E≤Lon≤140°E. During stage α12, the Aomori earthquake occurred on December 8, 2025, M=7.6, the epicenter of which is shown by a black asterisk on the α12 map.

8. Conclusions

A method for analyzing the evolution of low-frequency seismic noise field properties in a seismically active region is proposed. This method is based on studying the temporal and spatial distribution of extreme values of noise properties, which describe the complexity of its statistical structure. An increase in seismic hazard is interpreted as a simplification of the noise structure and its approach to the properties of white noise, as well as the synchronization of temporal variations in these properties. A method for identifying stages in the development of a seismic noise field is proposed based on a stepwise WTMM approximation of the change in the entropy of the distribution of weighted average densities of the distribution of extreme values of seismic noise properties. The method is applied to the analysis of noise properties on the Japanese islands over 29 years of observations, 1997-2025.

The main conclusions are that the spatial and temporal behavior of low-frequency seismic noise properties in Japan since 2020, and particularly in the last two years, reveal a number of indicators that point to the preparation of the next strong earthquake. Figure 5 shows a spike in the probability histogram of the distribution of time points when the values of the properties under consideration reach extreme values over the last six months. Figure 7 shows a pattern of quasi-periodic oscillations in the average integrated measure of the lead time of seismic noise property extreme points relative to the time points of strong earthquakes with a magnitude of at least 6 since 2015, with a peak value at the end of 2025, following a four-year interval when this pattern disappeared following the Tohoku mega-earthquake of March 11, 2011. Figure 8 shows a spike in the first principal component of the envelope amplitudes for the EEMD Hilbert-Huang decomposition of the daily median values of three seismic noise properties for the lowest-frequency IMF level number 8 (characteristic periods from 256 to 512 days) for the 2024-2025 time interval. This spike indicates synchronization of low-frequency variations in seismic noise properties. Figure 10 shows that the average entropy values of the weighted average probability density distribution of noise property extremes reach their maximum values from the end of 2020 to the end of 2025. Figure 11, which visualizes 12 stages of seismic noise field development on the Japanese islands, shows that the situation before the Tohoku mega-earthquake (stages (α04-α06)) is repeated at the final stage, α12, in that new “seismic hazard” spots appear, covering the entire Pacific coast of Japan. At the same time, the maximum probability density in the southern part of the studied area of the Japanese islands, which formed at stage α03, remains stable.

Taken together, these features are interpreted as indicators of increased current seismic hazard, with particular attention to the region 30°N≤Lat≤34°N, 136°E ≤Lon≤140°E. Naturally, this is not a prediction of the timing of the next strong earthquake; the essence of this study is to detect a trend of increasing seismic hazard.