3. Wavelet-Based Measure of Time-Series Non-Stationarity
Let us consider a random signal
where
is an integer discrete-time index. For a finite sample of the signal, the entropy [
20] could be defined as:
Here,
are the coefficients of orthogonal wavelet decomposition. Within a finite set of Daubechies wavelet bases [
20] with the number of vanishing moments from 1 to 10, the optimal basis is chosen from the minimum of (1).
The DJ-index was introduced in the problem of noise reduction and compression of information by setting all “small” wavelet coefficients to zero and performing an inverse wavelet transform [
14,
20]. This operation of noise reduction is performed for the optimal wavelet basis which was found from the minimum of entropy. The problem is how to define the threshold which separates “large” wavelet coefficients from “small” ones. In [
14], this problem was solved using a formula of the asymptotic probability of the maximum deviations of Gaussian white noise. The DJ threshold is given by the expression
, where
is the standard deviation of wavelet coefficients of the first detail level of wavelet decomposition. The first level is the most high frequency and is considered as the most noisy. Thus, in order to estimate the noise standard deviation, it is possible to calculate the standard deviation of the first-level wavelet coefficients
. The robust estimate was proposed:
where
N/2 is the number of wavelet coefficients on the first level, and
is the median. Formula (2) gives the relationship between the standard deviation and the median for the Gaussian random variable.
Finally, it is possible to define the dimensionless signal characteristic , , as the ratio of the number of the most informative wavelet coefficients, for which the inequality is satisfied, to the total number N of all wavelet coefficients. For stationary Gaussian white noise, which is a kind of ideal stationary random signal, the index . That is why the DJ-index could be called a measure of non-stationarity.
For each reference point, the index value
is calculated daily as a median of the values in the five nearest operational stations. Thus, continuous time series are obtained with a time sampling of one day at 50 reference points.
Figure 2 shows the graphs of the index
for 9 reference points (point numbers are indicated next to each graph).
Figure 3 shows a graph of average daily values of index
calculated for all 50 reference points. The red lines represent the linear trends of the average values calculated for the time fragments 1997–2003, 2004–2015 and 2016–2022. The boundaries of time fragments were chosen so that the continuity of the broken line at the union of linear trends was approximately observed. For the time intervals 1997–2003 and 2016–2022, linear trends coincided with the average values and it can be seen that during the intermediate time interval 2004–2015 there had been a significant decrease in the global average
. According to the results published in articles [
15,
17,
18], a decrease in the average value of the index
, which is a sign of the simplification of the statistical structure of noise (approximation of its properties to the properties of white noise), was associated with an increase in seismic hazard. Thus, the graph of average values in
Figure 3 can be interpreted as a transition to a higher level of global seismic hazard during 2004–2015.
4. Probability Densities of Extreme Values of Seismic-Noise Properties
Further analysis of the properties of seismic noise was based on the identification of regions that are characterized by the most frequent occurrences of extreme values for certain seismic-noise statistics. To do this, it was necessary to estimate the probability densities of the distribution over the space of extreme values by their values in the network of reference points. The following process was used. Denote by coordinates of the reference points. Let be the values of any property of the noise at the reference point with the number . We will be interested in maps of the distribution of probability densities of extreme values (minimums or maxima) of values over the Earth’s surface. To construct such maps, consider a regular grid of 50 nodes in latitude and 100 nodes in longitude covering the entire earth’s surface. Next, maps will be built for various seismic-noise statistics , the values of which will be estimated in time windows of various lengths. For each such window, we define the time index , which corresponds to its right end. For each time index , find the reference point with the number , at which the extremum of one or another property is realized, and let be the coordinate of the reference point with the number .
Let
be the distance between the points
and the coordinates of the nodes of the regular grid
. Then, the value of the probability density of extreme values for the quantity
at the node
of the regular grid can be calculated according to Gaussian distribution:
Here,
is the distance, which corresponds to an arc of 15 angular degrees, which is approximately equal to 1700 km, and
is a divisor that provides the normalization condition so that the numerical integral of function (3) over the Earth’s surface is equal to 1. Formula (3) gives an “elementary map” of the probability density of the property
extrema distribution for each current window with a time index
. This makes it possible to calculate the average probability density for a given time interval from
to
:
where
is the number of time windows labeled from
to
.
Another way to estimate the variability of the probability density of extreme properties of seismic noise is to construct histograms of the distribution of the numbers of control points, in which the statistics extrema are realized. The construction of such histograms in a moving time window, the length of which should include a sufficiently large number of initial time windows within which the noise properties are estimated, makes it possible to determine the reference points at which the minima or maxima of the seismic-noise properties are realized most often.
This method has the disadvantage that it does not give a direct distribution of the places of the most probable realization of the noise-property extrema in the form of a geographical map, although each reference point has geographical coordinates. On the other hand, this method makes it possible to compactly visualize the temporal dynamics and, in particular, to determine the time intervals when the places of concentration of extreme values for the noise statistics change sharply. To solve a similar problem of visualizing temporal dynamics using averaged probability densities (4), one would have to build a laborious and long sequence of maps for time fragments consisting of the same number of time windows.
Therefore, in the future, the probability densities of extreme values of noise properties will be presented simultaneously in two forms: as averaged densities according to Formula (4) for two time intervals, before and after 2011, and as a temporal histogram of the distribution of numbers of reference points in which extremes are realized. In this case, the number of base intervals (bins) of histograms is taken equal to the number of reference points, that is, 50, which makes it possible to visualize the dynamics of the emergence and disappearance of bursts of the probability of extreme values at each reference point.
Figure 4a,b show maps of the distribution of the minimum values of the DJ-index, calculated by Formulas (3) and (4) for windows with time stamps before and after 2011. The choice of 2011 as a boundary time mark is connected, as will be seen from the following presentation, with the fact that the correlation properties of seismic noise changed dramatically after two mega-earthquakes: on 27 February 2010,
M = 8.8 in Chile, and on 11 March 2011,
M = 9.1 in Japan. Previously, this change has been already detected in [
3,
18]. Increased attention to the areas where the minimum values
are most often realized, as was already mentioned above when discussing
Figure 3, is due to the fact that an increase in seismic hazard is associated with a simplification of the statistical noise structure and a decrease in
. For regions that are quite densely covered with seismic networks, for example, Japan, this property of statistics
makes it possible to estimate the location of a possible strong earthquake [
17].
However, the global network of seismic stations is sparse and therefore the selection of places where the probability of occurrence of small values is increased does not necessarily coincide with the known places of the occurrence of strong earthquakes. In
Figure 4a,b, the maxima of the probability-density distribution of the minimum values
are located in the vicinity of reference points with numbers 23, 25 and 37, 38 in the western part of the Pacific Ocean, where the epicenters of strong earthquakes and volcanic manifestations are really concentrated, including the largest volcanic eruption, Hunga Tonga-Hunga Ha’apai on 15 January 2022 [
21]. Note that in
Figure 4c, which shows a diagram of the change in the histogram of the numbers of reference points, in which the minimum value
was realized, one can see the switching of the histogram maximum from the reference point number 23 to the reference point number 38 just after 2011. As for the concentration of probabilities of minimum values
in the vicinity of reference points with number 46 in the southern Indian Ocean (in the vicinity of Kerguelen Island) and in the vicinity of point 40 in the Atlantic Ocean (Saint Helena) after 2011, these features can be associated with a mantle-plume activation regime [
22].
5. Estimation of Spatial Long-Range Correlations of Seismic-Noise Properties
The simplification of the statistical structure of seismic noise manifests itself in a decrease in the average value of the index
(
Figure 3), and it is also reflected in an increase in entropy (1) and a decrease in the width of the multi-fractal singularity spectrum support width (“loss of multifractality”). These changes in noise properties are indicators of the transition of a seismically active region (including the whole world) to a critical state [
3,
17].
Another indication of an increase in seismic hazard is an increase in the average correlation as well as the radius of strong correlations between noise properties in different parts of the system. In order to find out how the spatial scale of correlations changes with time, for each reference point we calculated the correlation coefficient between the index
values at this point and at all other reference points. We performed these calculations in a sliding time window 365 days long with a shift of 3 days. Next, for each reference point, we calculated the average value of the correlation coefficient with the values of
at all other reference points. Since the result of calculating such an average value depends on the time window, we obtained graphs of changes in average correlations for all reference points.
Figure 5 shows examples of graphs of changes in correlations for nine reference points.
From the graphs in
Figure 5, a general trend towards an increase in the average correlations for each reference point is noticeable. As a result of averaging all the graphs of the type shown in
Figure 5 for all 50 reference points, we obtained a general graph of average correlations, as shown in
Figure 6.
It can be noticed from the graph in
Figure 6 that the time interval 2004–2016, in which the index
decreases (
Figure 3), corresponds to an increase in average correlations between the properties of noise at various reference points and that there is a transition from small correlations before 2004 to fairly large correlations after 2016.
In order to assess the scale of strong correlations between seismic-noise properties, for each reference point we defined another point that had a maximum correlation with this point, provided that this maximum exceeded the threshold of 0.8, and we calculated the distance between these points.
Figure 7 shows the graphs of the changes in the average and maximum values of the distances between points for which the correlation maxima exceeded 0.8.
It can be seen from the graph in
Figure 7b that after 2011 there was a sharp increase in the maximum distance at which strong maximum correlations occur. Thus, the growth of average correlations, which is visible in
Figure 6, was supplemented by an explosive growth of the maximum scale of strong correlations.
It is of interest to identify those reference points in the vicinity of which strong correlations most often occur. As before, for this purpose we used estimates (3–4) of the probability density of the distribution over space of the maximum correlation values for the time intervals before and after 2011, as well as the distribution histograms of the numbers of control points in which the maximum correlation values occurred.
To select the length of the time window for calculating histograms, we recall that the correlation coefficients between the values at each reference point with the values of the DJ-index at other reference points are calculated in windows of 365 days with a shift of 3 days. Let us choose the length of the window for calculating histograms so that the dimensional length of the window is equal to five years. It is easy to calculate that this length will be 488 values of “short” time windows with a length of 365 days with an offset of 3 days. In this case, the histogram is calculated in a window of length 365 + 3·(488 − 1) = 1826 days, which is approximately equal to five years, given that every fourth year is a leap year. As for the number of base intervals (bins) of histograms, in this case there will not be 50, as in
Figure 4c, but 44, since for reference points with numbers 1–6 there are no realizations of correlation maxima.
Estimates of the position and temporal dynamics of places with an increased probability of maximum correlations are presented in
Figure 8. Comparison of
Figure 8a,b shows that the places of the concentration of maximum correlations have changed significantly since 2011. In particular, there has been a rapid shift of such a center from the Southwest Pacific to Central Eurasia, which can be seen in the histogram diagram in
Figure 8c.
6. Seismic-Noise Response to the Irregularity of the Earth’s Rotation
The connection between the irregular rotation of the Earth and seismicity was investigated in [
23]. The influence of strong earthquakes on the Earth’s rotation was studied in [
24,
25].
Figure 9 shows a graph of the time series of the length of the day for the time interval 1997–2022.
Further, the LOD time series is used as a kind of “probing signal” for the properties of the seismic process [
3,
15,
16,
17,
18]. To study the relationship between the properties of the seismic process and the irregularity of the Earth’s rotation, the squared coherences maxima between the LOD and the properties of seismic noise at each reference point were estimated. The coherences were estimated in moving time windows of 365 days with a shift of 3 days using a vector (two-dimensional) fifth-order autoregressive model [
15,
16,
17,
18]. In what follows, these maximum coherences will be referred to as the seismic-noise responses to the LOD.
Figure 8 shows examples of the LOD response graphs at nine reference points.
We define the integrated response of global seismic noise to the LOD as the average of the responses of the type shown in
Figure 10 from all 50 reference points. Let us compare the behavior of the average response to the LOD with the amount of released seismic energy. To do this, we calculate the logarithm of the released energy as a result of all earthquakes in a moving time window 365 days long with a shift of 3 days. In
Figure 11a,b, there are graphs of the logarithm of the released energy and the average response to the LOD. It is visually noticeable that bursts of the average maximum coherence in
Figure 11b precede the energy spikes in
Figure 11a. In order to quantitatively describe this delay, we calculated the correlation function between the curves in
Figure 11a,b. The graph of this correlation function for time shifts of ±1200 days is shown in
Figure 11c.
The estimate of the correlation function in
Figure 11c between the logarithm of the released seismic energy in
Figure 11a and the average response of seismic noise to the LOD had a significant asymmetry and was shifted to the region of negative time shifts, which corresponded to the coherence maxima being ahead of the maxima of seismic energy emissions. The correlation maximum corresponded to a negative time shift of 534 days. This estimate of cross-correlations produced a foundation to propose that the burst shown in
Figure 11b by the magenta line may precede a major earthquake with an average delay of 1.5 years.
The curve in
Figure 11b can be split into three sections. The first two intervals with time marks of the right ends of windows before and after 2012 differ significantly from each other by their mean values presented by red lines.
The third time interval in
Figure 11b is presented by the magenta line and refers to the processing of the most recent data which correspond to right-hand end windows after 2021, i.e., referring to the time span of 2020–2022. A short third segment was identified based on the results of data processing for 2021–2022 due to an unusually large spike in the response of noise properties to the LOD. The maximum value of the response in
Figure 11b was reached at the time corresponding to the date 9 May 2022. The “naive” forecast, in which 534 days are added to this date, gives the most probable date for the future strong event as being 24 October 2023. A more realistic prediction would be to blur this date with some uncertainty interval. However, it is not yet possible to give such an interval due to the short history of using the seismic-noise response to the LOD as a predictor.
Similar to how it was done to identify the places of concentration of the maxima of spatial correlations of seismic noise in
Figure 8, we can search for places where the maximum response of seismic-noise properties to the uneven rotation of the Earth is most likely. The results of such a search are shown in
Figure 12. Comparison of
Figure 12a,b shows that the maximum response to the LOD is concentrated in the subarctic regions.
In addition to the issue of the synchronization of seismic-noise properties, which are reflected in the graphs in
Figure 5,
Figure 6 and
Figure 7, we also studied the synchronization of seismic-noise responses to the irregularity of the Earth’s rotation. Previously, this issue has been considered for the synchronization of seismic-noise responses to LODs in Japan and California [
16].
For each reference point, we estimated the correlation of the response at this point with the responses at all other points. It should be taken into account that the correlated responses themselves were calculated in sliding time windows of 365 days with a shift of 3 days. These were “short” time windows. In this case, we needed to take a certain number of consecutive “short” windows and form a “long window” from them, such that it contained a sufficiently large number of estimates of the noise responses to the LOD in the “short” time windows. Next, 250 adjacent “short” windows were selected to form a “long” window. Correlation coefficients between the responses to the LOD at various reference points were calculated in a sliding “long” window with a shift of 1, that is, in a time window with a length of 365 + 3·(250 − 1) = 1112 days, which is approximately 3 years. A shift of 1 count meant a real shift of 3 days.
The goal of such calculations for each reference point was the average correlations of responses over all other reference points. As a result of such averaging, 50 dependences of the same type as shown in
Figure 5 were obtained. In order to obtain an integrated measure of the correlation of global noise responses to the LOD, the average dependences for all reference points needed to be averaged secondarily. As a result, we obtained the graph shown in
Figure 13.
Figure 13 shows the change in time of the measure of synchronization of the response of seismic noise to the irregularity of the Earth’s rotation. This graph also shows the growth of correlations, similar to the growth of the initial correlations in
Figure 6.