Spatial Correlations of Global Seismic Noise Properties

: A study of global seismic noise during 1997–2022 was carried out. A property of waveforms known as the Donoho–Johnston (DJ) index was used, which separates the values of the wavelet coefﬁcients into “small” and “large”. For each reference point in an auxiliary network of 50 points, a time series was calculated with a time step of one day for the median of the values at the ﬁve nearest stations. In a moving time window of 365 days, correlations between the index values at the reference points were calculated. A decrease in the average values of the DJ-index and an increase in correlations were interpreted as a sign of an increase in global seismic danger. After 2011, there was a sharp increase in the maximum distances between reference points with large correlations. The high amplitude of the response of the DJ-index to the length of the day for 2020–2022 could predict a strong earthquake in the second half of 2023. The purpose of this study was to improve the mathematical apparatus for assessing the current seismic hazard according to the properties of seismic noise.


Introduction
Information about seismic noise makes it possible to study the processes preceding strong earthquakes [1][2][3]. Atmosphere and ocean (cyclone movement and the impact of ocean waves on the shelf) are the main sources for the energy of seismic noise [4][5][6][7][8][9][10]. At the same time, processes inside the earth's crust are reflected in changes in the properties of seismic noise and studying these properties helps in investigating the structural features of the earth's crust [11][12][13].
The Donoho-Johnston index (DJ-index) of a random signal can be defined as the ratio of the number of the coefficients of orthogonal wavelet decomposition, which in absolute value exceed the threshold introduced in [14], to the total number of wavelet coefficients. The threshold that separates "large" wavelet coefficients from others was defined first in [14] and is used to shrink noise from signals and images using wavelets. It has turned out that the DJ-index has the ability to most clearly highlight the effects of the spatial correlation of seismic noise, and it outperforms other noise statistics (such as entropy and the width of the multi-fractal spectrum of the singularity) in terms of highlighting predictive effects [15][16][17][18]. This article presents a detailed study of the DJ-index of global seismic noise, both the change in time and space of the correlation properties of noise and its response to the irregularity of the Earth's rotation, as well as the possible use of this response to assess the current seismic hazard.

Wavelet-Based Measure of Time-Series Non-Stationarity
Let us consider a random signal ( ) x t where 1,..., t N = is an integer discrete-time index. For a finite sample of the signal, the entropy [20] could be defined as: Here, k c 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 2 ln , where σ is the standard deviation of wavelet coefficients of the first

Wavelet-Based Measure of Time-Series Non-Stationarity
Let us consider a random signal x(t) where t = 1, . . . , N is an integer discrete-time index. For a finite sample of the signal, the entropy [20] could be defined as: Here, c k 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 c (1) k . The robust estimate was proposed: where N/2 is the number of wavelet coefficients on the first level, and med 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 γ, 0 < γ < 1, as the ratio of the number of the most informative wavelet coefficients, for which the inequality |c k | > T DJ 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 γ = 0. 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).
where N/2 is the number of wavelet coefficients on the first level, and med 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  , 01   , as the ratio of the number of the most informative wavelet coefficients, for which the inequality k DJ cT  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 0 =  . 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).     [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.  [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.

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 Here, h is the distance, which corresponds to an arc of 15 angular degrees, which is approximately equal to 1700 km, and ( ) P t is a divisor that provides the normalization

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 ζ k , k = 1, . . . , m = 50 coordinates of the reference points. Let Z k be the values of any property of the noise at the reference point with the number k. We will be interested in maps of the distribution of probability densities of extreme values (minimums or maxima) of values Z 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 Z, the values of which will be estimated in time windows of various lengths. For each such window, we define the time index t, which corresponds to its right end. For each time index t, find the reference point with the number k * t , at which the extremum of one or another property Z is realized, and let ζ * t be the coordinate of the reference point with the number k * t . Let ρ(ζ k , r ij ) be the distance between the points ζ k and the coordinates of the nodes of the regular grid r ij . Then, the value of the probability density of extreme values for the quantity Z at the node r ij of the regular grid can be calculated according to Gaussian distribution: Here, h is the distance, which corresponds to an arc of 15 angular degrees, which is approximately equal to 1700 km, and P(t) 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 Z extrema distribution for each current window with a time index t. This makes it possible to calculate the average probability density for a given time interval from t 0 to t 1 : where n(t 0 , t 1 ) is the number of time windows labeled from t 0 to t 1 .
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 Z 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].
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].

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.

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

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.

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. 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 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 Figures 5-7, we also studied the synchronization of seismicnoise 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.
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 Figures 5-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" 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. 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.

Conclusions
The results of a detailed analysis of the properties of seismic noise based on the use of the Donoho-Johnston index, defined as the proportion of the maximum values of the Mean correlation Figure 13. The average value of the correlation coefficients between the seismic-noise index γ responses per LOD for each of the 50 control points with responses at other points. 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.

Conclusions
The results of a detailed analysis of the properties of seismic noise based on the use of the Donoho-Johnston index, defined as the proportion of the maximum values of the moduli of the wavelet coefficients separating the "noise" coefficients, was presented. The main results of the conducted research were as follows: (1) Three characteristic time intervals for the behavior of seismic noise were identified: 1997-2003, 2004-2015 and 2016-2022. (2) For these time intervals, there was a parallel decrease in the average values of the DJ-index (interpreted as a simplification of the structure of seismic noise and an approximation of its properties to white noise) and an increase in average spatial correlations. According to general ideas about the behavior of complex systems, this means that the system is approaching a critical state, i.e., an increase in the probability of strong earthquakes. (3) After 2011 (after the Tohoku mega-earthquake in Japan on 11 March 2011), there was an explosive increase in the maximum distances at which strong correlations occur. Such behavior is known in the theory of critical phenomena as "an increase in the radius of fluctuation correlations" [26], which also heralds a sharp change in the state of complex systems. (4) The estimation of the correlation function between the logarithm of the released seismic energy and the average response of seismic noise to the irregularity of the Earth's rotation shifted towards negative time delays, which corresponds to the shift of the coherence maxima towards the peaks of the seismic-energy release. The maximum correlation fell on a time shift of 534 days. This estimate assumes that the burst of LOD response corresponding to the data analysis in 2020-2022 precedes a strong earthquake with an average delay of 1.5 years. (5) From the estimation of the change in time of the synchronization measure of the seismic-noise response to the irregularity of the Earth's rotation (Figure 13), one can also see a growth of correlations, similar to the growth of the initial correlations in Figure 6.
It should be noted that the noise level in real observations of the seismic background is inevitably extremely high. Transition to the low-frequency part of the spectrum allows the elimination of most of the man-made noise. Moreover, this transition makes it possible to study that part of the spectrum that is in the intermediate, little-studied region of periods between traditional seismology and gravimetry. The author's previous studies give hope that it is in this region of periods that the effects preceding strong earthquakes can be detected. This work is presented as a further development of these methods. The complexity of the data and processes underlying the formation of low-frequency seismic noise in the Earth's crust does not allow direct use of conventional statistical procedures for estimating confidence intervals, standard deviations and other measures of uncertainty. It is possible that further study of the statistical structure of low-frequency seismic noise will enable us to build such models that will make it possible to estimate the measures of uncertainty of the obtained statistical inferences, for example, using the bootstrap technique.
Funding: This research received no external funding.