CC BY
37
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 9, ES2005, doi:10.2205/2007ES000220, 2007
Using modern seismological data to reveal earthquake precursors
G. A. Sobolev1,2, and A. A. Lyubushin1
Received 21 June 2007; accepted 16 September 2007; published 12 December 2007.
[1] Records obtained at the IRIS world-wide system of broadband seismic stations before strong earthquakes were investigated with a purpose of detecting hidden periodicities, multiple coherence effects and seeking for asymmetric impulses within a minute range of periods. The stations are located at different distances from the epicenters of earthquakes.
The initial realizations consisted of discrete measurements with a sampling rate of 20 Hz and the total volume of analyzed data exceeded 25 Gb. We used various programs of processing and analyzing time series: revealing of hidden periodicities in the sequences of peak values at a given level; wavelet analysis of microseisms flow; search for multiple coherence effects based on Fourier and wavelet approaches; estimates of spectral coherence measures evolution of variations of multi-fractal singularity indexes and others. Asymmetric pulses about 3-10 min long did arose several days before the Kronotskoe 05.12.1997 (M = 7.8), Neftegorskoe 27.05.1995 (M = 7.0), and Hokkaido 25.09.2003 (M = 8.5) earthquakes. Intervals of a stable manifestation of several periods (tens of minutes) of pulses were found before the Kronotskoe and Hokkaido events. Synchronization of microseismic oscillations at different stations was detected starting several days before Hokkaido and Sumatra 26.12.2004 (M =
9.2) earthquakes. Comparison of records obtained at different stations allows estimating the regional and local peculiarities of the anomalies. It is assumed that the nature of these phenomena is related to self-organization properties of the seismic process. The periodic vibrations, asymmetric pulses and synchronization intervals are indicators of the unstable state of a seismically active region and could be regarded as earthquake precursors.
INDEX TERMS: 7209 Seismology: Earthquake dynamics; 7223 Seismology: Earthquake interaction, forecasting, and prediction; 7294 Seismology: Seismic instruments and networks; KEYWORDS: IRIS world-wide system, microseismic oscillations, wavelet analysis, earthquake precursors.
Citation: Sobolev, G. A., and A. A. Lyubushin (2007), Using modern seismological data to reveal earthquake precursors, Russ. J. Earth. Sci, 9, ES2005,doi:10.2205/2007ES000220.
Introduction
[2] Numerous oscillating fields of different nature affect the earth in extremely wide range of periods. When this happens certain types of energy partially transform to other types. For example, magnetic waves energy coming to the earth from the outside causes elastic oscillations because of inverse piezoelectric and seismoelectric effects; elastic stresses in the earth appear with the coming heat owing to thermoelastic coupling coefficients and others. The intensity of external effects may be small as compared to the forces acting inside
1Schmidt Institute of the Physics of the Earth RAS, Moscow, Russia
2Geophysical Center RAS, Moscow, Russia
Copyright 2007 by the Russian Journal of Earth Sciences.
ISSN: 1681-1208 (online)
the earth but the extent of their influence depends on the energy saturation of rocks and cannot be explained by linear effects.
[3] The appearance of rhythms under the effects caused by internal and external sources that is synchronization was discussed in geophysical papers long ago. Solar activity, earth tides, and climate are known to influence seismicity.
[4] The problem remains to be solved of the external effect threshold that is sufficient to synchronize the process caused by more powerful forces. A system open in terms of energy and sensitive to minor external effects apparently is in metastable state. As the system approaches instability, the efficient external effect threshold lowers. The earth, however, is continuously affected by noise from natural and artificial sources. Therefore the threshold of effective influence that can be detected (including trigger mechanism) is apparently of finite value exceeding the noise level.
[5] The effect of hidden periodic oscillations in weak earth-
Figure 1. Position of the IRIS stations whose records were analyzed before the Kronotskoe and Neftegorskoe earthquakes. The epicenters of the earthquakes are shown by stars.
quakes and microseisms series discovered in [Sobolev, 2003, 2004] falls into the class of phenomena under discussion. In principle, it can be considered in terms of self-organized criticality concept (SOC) [Bak et al., 1989; Sornette and Sammis, 1995], which attaches much importance to the emergence of remote correlation of seismic events (collective behavior). However the physical mechanism of the possible remote correlation in the context of seismology is not clear yet; general theories of catastrophes and phase transitions in open-energy systems invite more detailed studies for heterogeneous environments.
[6] From the end of the 1990s and after the global system of wide-band seismic stations was established, seismic noise was studied in the range of 102-103 seconds. The authors [Tanimoto et al., 1998] believed that oscillations in the solid earth appeared under the influence of atmospheric pressure variations. The authors of alternative hypothesis [Kobayashi and Nishida, 1998] assume that oscillations are caused by numerous weak earthquakes of energy below seismic stations sensitivity. These studies as well as others show that minute range oscillations are actually permanent including time intervals free of strong earthquakes. However sequences of individual impulses divided with intervals of their missing evidently were not revealed. It may be related to the fact that most researchers used Fourier spectral analysis, which is not intended for research in non-stationary process comprising bursts of different amplitudes and duration. The use of the program complex described below including wavelet analysis appears to be more promising.
Initial Data
[7] We analyzed seismic records of the vertical component with discretization frequency of 20 Hz of wide-band stations IRIS before four strong earthquakes: Neftegorskoe earthquake of 27 May 1995 with coordinates [52.55°N, 142.75°E], M = 7.0; Kronotskoe earthquake of 5 December 1997 [54.64°N, 162.55°E], M = 7.8; Hokkaido earthquake of
25 September 2003 [41.81°N, 143.91°E], M = 8.3; Sumatra earthquake of26 December 2004 [3.32°N, 95.85°E], M = 9.2. The data were kindly given by the Geophysical Service RAS. Horisontal component records were analyzed incidentally. In the studies of Kronotskoe and Neftegorskoe earthquakes we used records of stations PET, MAG, YSS, YAK, ARU, and OBN; their location is shown in Figure 1. For Hokkaido earthquake the system comprised stations ERM, MAJ, INC, MDJ, BJT, PET, YSS shown in Figure 2 and OBN station (Figure 1). For Sumatra earthquake, the system of stations CHTO, KMI, XAN, COCO, PALK, MBWA, DGAR, DAV, QIZ is shown in Figure 3. Stations selected for the studies were located at distances ranging from 70 km to 7160 km from the above-mentioned earthquakes and in different seis-mogeological conditions. The total volume of analyzed data was more than 24 Gb.
[8] Amplitude-frequency characteristics of IRIS channels provide permanent sensitivity of recording the velocity of bases displacement in the period range of 0.3-357 seconds [Starovoit and Mishatkin, 2001]. Oscillations up to peaks of
Figure 2. Position of the IRIS stations whose records were analyzed before the Hokkaido earthquake. The epicenter of the earthquake is shown by star.
12 and 25 hours caused by earth tides are reliably registered in spite of the sensitivity drop in the range of longer periods. A typical spectrum is shown in Figure 4. We note three intervals of period values in the spectrum plot. In the short-period range up to 6 minutes, a drop in oscillation strength is observed that is a result of gradual decrease in the influence of microseisms of oceanic origin and weak earthquakes. In the period range of hundreds of minutes, the influence of earth tides is noted which is corroborated by peaks of 1440 and 720 minutes corresponding to diurnal and semidiurnal oscillations. We analyzed microseismic signals in the relatively low-noise minute range marked with an arrow in Figure 4.
Methods
[9] To obtain the results given in the paper we used four major methods to study observation series.
Method 1. Revealing Asymmetric Impulses
[10] To separate high-amplitude low-frequency impulses A. A. Lyubushin created a program [Sobolev and Lyubushin,
2006] that sequentially performs the following operations: signal aggregation by 20 times; removal of low-frequency Gaussian trend with scale parameter (averaging radius) of 1000 counts (seconds) to suppress oscillations caused by earth tides; calculations of Gaussian trend with parameter of
100 seconds to suppress oscillations of second range caused by microseisms of oceanic origin and earthquakes.
[11] The operations to calculate and to remove trends are as follows. Let X (t) be arbitrary limited integrated signal with continuous time. Let kernel averaging with scale parameter H > 0 be called the average value of X (t|H) in the time moment t calculated by formula:
x(t|H) = y x(t+h^) • V’m/y (1)
where ^(£) is arbitrary non-negative limited symmetric integrated function called averaging kernel [Hardle, 1989]. If ^(£) = exp(-£2), value X(t|H) is called Gaussian trend with parameter (radius) of averaging H [Hardle, 1989; Lyubushin,
2007].
[12] As a result of the above-described preliminary operations we obtain a signal with discretization interval of 1 sec with power spectrum in the period range approximately from 3 to 30 minutes The same result can be obtained with the common band Fourier filtration but Gaussian trends are more preferable because side effects caused by filter discrimination are missing and besides it is easier to control boundary effects caused by the finite character of the sampling being filtered.
[13] To analyze impulse sequence in the quantitative sense a program was developed to separate them automatically. We used Haar expansion by wavelets [Daubechies, 1992; Mallat, 1998]. After direct Haar wavelet transformation only a small part (1-a) of wavelet coefficients maximal in magnitude was left with a = 0.9995. Then reverse wavelet
Figure 3. Position of the IRIS stations whose records were analyzed before the Sumatra earthquake.
p, p € [0, 2n] and factor ^ > 0 (describing Poisson part of intensity) are the model parameters. Thus Poisson part of intensity is modeled by harmonic oscillations.
[15] Owing to considering an intensity model richer than for a random series of events and with harmonic component of given frequency w, the likelihood logarithmic function increment of point process [Cox and Lewis, 1966] is equal to [Lyubushin et al., 1998]:
A ln L(a, p| w) = ^2ln(1 + a cos(w^ + p))
(3)
+ N ln( wT/[wT + a( sin(wT + p) — sin(p))]) .
Here ti is the sequence of time moments of separated local maximums of the signal inside the window; N is the number of them; T is the length of time window. Suppose
R(w) = max A ln L(a, p| w),
a,v (4)
0 < a < 1, p € [0, 2n] .
[16] Function (4) may be considered as spectrum generalization to the sequence of events [Lyubushin et al., 1998]. The plot of the function shows how much more favorable
The epicenter of the earthquake is shown by star.
transformation was performed and as a result a sequence of impulses was separated with large amplitudes divided from one another with intervals of constant values that had been filled with noise before. In wavelet analysis, this operation is known as denoising [Daubechies, 1992; Mallat, 1998]. Haar wavelet was chosen for this operation because of simplicity of subsequent automated separation of rectangular impulses. The number of separated impulses and the extent of noise rejecting depend on the selected compression level a.
Method 2. Detecting Periodic Components in the Sequence of Events
[14] The method is intended for detecting periodic components in the sequence of events and was proposed in [Lyubushin et al., 1998]. Intensity model of events sequence was considered (in this case, time moments of essential local maximums, that is overshoots of microseisms time series), which supposedly contains harmonic component.
A(t) = ^ • (1 + a • cos(wt + p)) , (2)
where frequency w, amplitude a, 0 < a < 1, phase angle
the periodic model of intensity is as compared to a purely random model. Maximum values of function (4) indicate frequencies present in the event sequence.
[17] Let t be the time of the right edge of the moving time window of preset length TW. Expression (4) is actually a function of two arguments: R(w,t|Tw), which can be visualized in the form of 2-D map or 3-D relief on the plane of arguments (w,t). This frequency-time diagram allows us to study the dynamics of appearance and development of periodic components inside the event flow under investigation [Lyubushin, 2007; Sobolev, 2003, 2004; Sobolev et al., 2005].
Method 3. Initial Data Transformation Into Variations of Hurst Generalized Exponents
[18] Multifractal measure of synchronization or the evolution of the spectral measure of Hurst generalized exponent variations coherent behavior was studied with different station sets. We only touch upon the method briefly, referring for details to [Lyubushin, 2007; Lyubushin and Sobolev,
2006].
[19] Note that the analysis of multifractal characteristics of geophysical monitoring time series is one of promising lines of data studies in the solid earth physics. [Currenti et al., 2005; Lyubushin, 2007; Telesca et al., 2005]. It is determined by the fact that multifractal analysis is capable of investigating signals that are nothing more than white noise or Brownian motion in terms of covariation and spectrum theory.
[20] Let X(t) be a signal. We define the amplitude as variation measure ^(i,5) of the signal behavior X(t) at interval [t, t + 5] :
u(t,5) = max X(s) — min X(s) .
t<s<t+S t<s<t+S
(5)
ln(5)
M(5,q) = M{ (MM))*} .
(7)
K(q) = lim
5^0
ln M (5, q) ln(5)
(8)
[21] Holder-Lipschitz index h(t) for point t is defined as limit: ( )
ln( u(t, 5))
h(t) = lim \ ’ , (6)
that is in the vicinity of point t signal variation measure ^(t,5) with 5 ^ 0 descents by law 5h(t).
[22] Singularity spectrum F(a) is defined [Feder, 1989] as fractal dimensionality of a set of points t, for which h(t) = a (that is having the same Holder-Lipschitz index equal to a).
[23] The existence of singularity spectrum is ensured not for all signals but only for the so-called scale-invariant ones. If X(t) is a random process, let us calculate the average value of measures ^(t, 5) in exponent q:
[24] Random process is called scale-invariant if M (5, q) with 5 ^ 0 descents by law 5K(q), that is to say the limit exists:
[25] If the dependence K(q) is linear: K(q) = Hq, where H = const, 0 < H < 1, then the process is monofractal.
Figure 4. Power spectra estimate for micro-seismic oscillations at the station PET for time interval 20.11-05.12.1997 (strictly before Kronotskoe earthquake) after coming to 30-seconds sampling time interval. Double-headed red arrow indicates the main period range to be investigated.
Specifically for Brownian movement H = 0.5. Process X(t) is multifractal if dependence K(q) is non-linear.
[26] The idea of raising to different powers q in formula (7) implies that different weights may be given to time intervals with big and small measures of signal variability. If q > 0, then the major contribution into the average value M(5, q) is made by time intervals with great variability, whereas time intervals with small variability make maximum contribution with q < 0.
[27] If we estimate spectrum F(a) in moving time window its evolution may give information on the variation of the series random pulsations. Specifically the position and width of the spectrum carrier F(a) (values amin,amax and Aa = amax — amin and a* are given to function F(a) by maximum: F(a*) = max F(an are noise characteristics. Value a* can
a
be called generalized Hurst exponent. For monofractal signal the value of Aa is to be equal to zero, and a* = H. As for the value of F(a*), it is equal to fractal dimensionality of points in the vicinity of which scaling relationship (8) is fulfilled.
[28] Below to calculate singularity spectrum F(a) we applied the detrended fluctuation analysis [Kantelhardt et al., 2002] in programs given in detail in [Lyubushin, 2007; Lyubushin and Sobolev, 2006].
[29] Commonly F(a*) = 1, but in some windows we found F(a*) < 1. Recall that in the general case (not only for time series analysis) value F(a*) is equal to fractal dimensionality of multifractal measure support [Feder, 1989].
Figure 5. Noise components of 30-s discretized seismic records at the five stations after threshold filtering and removal of the first three detail levels (the resulting range of periods is 8-128 min).
Method 4. Calculations of Spectral Measure of Coherence
[30] To detect coherent elements of behavior that may have phase shift and be observed for several stations at one and the same time we used the method that uses canonical coherences estimate in the moving time window developed in [Lyubushin, 1998] to search for earthquake precursors from low-frequency geophysical monitoring data. In papers [Lyubushin et al., 2003, 2004], this method was applied to the analysis of multidimensional hydrological and oceanographic time series. In papers [Lyubushin and Sobolev, 2006; Sobolev and Lyubushin, 2007] this measure was used to analyze synchronization of microseismic oscillations. Technical details of its calculations may be found in [Lyubushin, 1998,
2007].
[31] Spectral measure of coherence A(t, w) is built as modulus of the product of canonical coherence component by component.
(9)
[32] Where q is the total number of time series analyzed together (dimensionality of multidimensional time series); w is frequency; t is time coordinate of the right edge of the moving time window comprising a certain amount of adjoining points; Vj (t, w) is canonical coherence of j-scalar time series, which describes the strength of coherence between this series and all other series. Value |vj (t, w)|2 is generalization of common square spectrum of coherence between two signals when the second signal is not scalar, but vector. Inequality 0 < |vj (t, w)| < 1 is fulfilled and the closer is value |vj (t, w)| to one, the stronger are linearly connected variations at frequency w in time window with coordinate t of j-series with analogous variations in all other series. Correspondingly value 0 < A(t, w) < 1 owing to its construction describes the effect of cumulative (synchronous, collective) behavior of all the signals.
[33] Note that, from the construction, the value of A(t, w) belongs to interval [0,1], and the closer are corresponding values to one the stronger is the connection between variations of the multidimensional time series components Z(t) at frequency w for time window with coordinate t. It should be emphasized that comparing absolute values of statistics
j=i
A(t, w) is only possible for one and the same number q of time series processed simultaneously since by formula (9) with the increase of q value A decreases, being product of q values that are below one.
[34] To implement this algorithm the spectral matrix assessment of the initial multidimensional series is required in each time window. In the following, preference is given to the use of vector autoregressive model of the 3rd order [Marple, 1987]. We took time window length equal to 109 values to obtain relationship A(t, w). Since each value of a* was obtained from time window of a length of 12 hours and the shift of those windows is 1 hour, then the length of time window to assess the spectral matrix is (109-1)-1+12 = 120 hours = 5 days.
Results
Asymmetric Impulses
[35] Method 1 was used to separate possible high-amplitude impulses in the minute range of periods. In Figure 5 the results are given of processing initial 20 Hz records of five stations (Figure 1) before Kronotskoe earthquake in which components were separated with periods ranging from 8 to 128 minutes. ARU data were not used because of the noise of technogenic origin in the day time. From the figure, it follows that in the time interval under investigation 84 hours before the earthquake no considerable variation of noise level was noted in any of the five stations. The amplitude of impulses at each of the stations only makes several percent of the level of common microseisms in the range of 2-6 seconds. We studied the variation of noise level with the period increase in all of the stations. It was revealed that microseisms have maximum amplitude in the range of 3-5 s, which is in agreement with the model of their oceanic origin. As the period lengthens the amplitude drops, so that with period of 1 minute it, on the average, decreases by 20 times. Then the decrease of amplitude slows down and after a period of 5-10 minutes it gradually starts growing. In the qualitative sense, the spectra of oscillations intensity have the same form as the spectrum for PET station that is shown in Figure 4. However, the farther is the station from PET, the slower is the decrease with lengthening of periods in the range from seconds to first minutes. Thus values of oscillations velocities shown in y-axis of plots in Figure 5 allow us to draw the following preliminary conclusion. The fact that their value gradually decreases as the distance of the station from Kamchatka (PET) increases and their level is higher for an order of magnitude in PET station as compared to other stations suggests that their source was located in the Pacific seismically active zone.
[36] The feature of PET station is asymmetric impulses of negative polarity. The dynamics of the impulses quantity change automatically separated by the program in 30-day interval before Kronotskoe earthquake is shown in Figure 6. Curve 1 shows the level of “common” microseisms of second range, where D is dispersion of initial record in successive
310 320 330 340
Days from 1 January 1997
Figure 6. Comparison of the microseismic level in the range of periods of the order of seconds (1) with the daily number of pulses of negative (2) and positive (3) polarities. The arrow shows the occurrence time of the Kronotskoe earthquake.
windows of length of 4 seconds (80 initial counts with sampling frequency of 20 Hz). Two upper curves show variations of the amount of impulses found by the program and having positive and negative polarity of minute range for 1 day with a step of 0.1 day. Gradual increase in the amount of impulses is noted especially of negative polarity N(—) before the earthquake.
[37] We compared oscillations of the amount of asymmetric impulses and the level of microseisms of second range of periods; no correlation between them was established. It means that impulses were not caused by storms in the seas and oceanic areas. The influence of more distant meteorological phenomena, for example in the Atlantic regions, is not excluded, but this possibility was not checked. Correlations between the quantity of impulses and semidiurnal, diurnal and two-week earth tides was not revealed either.
[38] Let us consider the structure of individual impulses and their sequence. One and the same section of PET station record of 3 December of about 3-hour duration is shown in plots 1 and 2 (Figure 7). The lower plot is the record with discreteness 1 s filled with microseism variations of the second range of periods. Plot 2 is the low-frequency constituent after high frequencies were rejected with smoothing by Gaussian kernels with averaging radius of 100 s. In
■2x10
0 2x103 4x103 6xl03 8X103 104s
Figure 7. Fragments of PET (1, 2) and OBN (3) records before the Kronotskoe earthquake: (1) initial record; (2, 3) variations within the range of periods of the order of minutes.
this plot you can see impulses with the following features: 1) duration of individual impulse is ~10—15 min, and intervals between them make ~40 min; 2) The impulse form changes from almost symmetrical to one-polar with gradual suppression of positive phase. Note that impulse amplitude in plot 2 is an order less than microseisms level in plot 1 and therefore they cannot be seen in the lower plot. A record fragment of the same time and processed in a similar way from Obninsk station located at the East European platform at a distance of 6500 km from Petropavlovsk station (see Figure 1) is shown in plot 3 (Figure 6) for comparison. Impulses of the type of plot 2 are not detected and the amplitude of oscillations is two less orders.
[39] Similar results were obtained before Neftegorskoe earthquake for Yuzhno-Sakhalinsk station (YSS), which was the nearest to the epicenter. In Figure 8, record fragments of microseisms (plot 1) and of low-frequency component (plot 2) of this station and of Obninsk station (plot 3) for comparative purposes are shown. In plot 2 both asymmetric and symmetric impulses of duration ~10-15 minutes are separated. Comparing Figure 8 and Figure 7 we note the following: durations of individual impulses (plots 2) coincide in both cases; as distinct from Kronotskoe earthquake, intervals between sequential impulses before Neftegorskoe earthquake are not regular; impulse amplitude in Figure 8 is approxi-
mately one twentieth as much as compared to Figure 7; microseisms amplitude (plot 1 Figure 8) is approximately one twenty-fifth as much as the one in Figure 7; oscillations amplitude at OBN station (plots 3) is comparable in the both cases (~5-7 nM s-1), and their structure differs from that of stations YSS and PET. Before Neftegorskoe earthquake the amount of negative polarity impulses increased in the last 5 days against the undisturbed background of second-range microseisms.
[40] Before Hokkaido earthquake the records of two stations shown in Figure 2 intense impulses were revealed with amplitudes above diurnal and semidiurnal tidal oscillations. It is demonstrated in Figure 9, where 4-day records (96 hours) are given for the period of 16-19 September 2003: plot ERM is of Erimo station located in the southeast of Hokkaido in subduction zone and actually in the epicentral zone of the earthquake; plot PET is of station Petropavlovsk on Kamchatka shore located in subduction zone; plot MDJ is of station Mudanjiang located in the continent in northeastern China. In ERM and PET plots, individual impulses of positive and negative polarity as well as series of impulses close in time can be seen against diurnal and semidiurnal tides. Positive polarity series is noted in the area of 1620 hours in ERM plot; negative polarity series in PET plot takes the interval of 82-96 hours. We note three features that in our opinion are significant: 1) impulses were only
11111 s"1 «
■800
0 2xl03 4xl03 6x103 8x103 1 04s
Figure 8. Fragments of YSS (1, 2) and OBN (3) records before the Neftegorskoe earthquake: (1) initial record; (2, 3) variations within the range of periods of the order of minutes.
manifested at stations ERM and PET located in subduction zone, 2) the times of both single impulses and series of impulses do not coincide in ERM and PET plots, 3) impulses amplitude in ERM plot is on the average greater as in magnitude and with respect to the span of tidal variations both as compared to plot PET. From data given in Figure 9, it is reasonable to assume that the sources of impulse oscillations are located in subduction zone. Since we do not know the actual sensitivity of a number of stations in the minute and hour range of periods, we give attention to relative amplitudes of oscillations in comparison with tides or microseisms of second range.
[41] More detailed structure of impulses of stations ERM and PET is shown in Figure 10. We changed impulses polarity in ERM plot artificially to inverse as compared to plot 9 so that comparing is more convenient. Low-frequency oscillations of hour range of periods caused by earth tides were rejected by deducting Gaussian trend with the radius of 1000 s. Intervals between sequential impulses make first thousands of seconds. From more detailed time base of the impulses 1, 2, 3, and 4, it can be seen that the time of impulse rise to extreme values is not the same and lies in the interval from 100 to 200 seconds. This interval is within standard frequency range of IRIS stations IRIS [Starovoit and Mishatkin, 2001]. From Figure 10 it follows that the comparative amplitude of impulses with respect to high-frequency noise of second-range microseisms is greater at ERM stations. Since BHZ measuring channels of IRIS stations record the vertical component of displacement velocity, it is believed that the presented impulses correspond to one-polar vertical step of the base movement. Comparison with horizontal component
Figure 9. Fragments of ERM, PET, MDJ records before the Hokkaido earthquake. The ordinates show the velocity of displacement V in arbitrary units.
Figure 10. Structure of asymmetric pulses in the minutes-range of periods recorded by stations ERM and PET before the Hokkaido earthquake. The ordinates show the velocity of displacement in arbitrary units.
record showed that horizontal component amplitudes of N-S and E-W impulses are almost one less order as compared to vertical components and are comparable to high-frequency noise amplitudes.
[42] Before Sumatra earthquake, quasi-sinusoidal symmetric oscillations in the minute range of periods were only detected in the records of stations shown in Figure 3.
Periodic Oscillations
[43] As a result of applying method 2 to the analysis of microseism variations before Kronotskoe earthquake of 5 December 1997 with M = 7.8 and coordinates [54.64° N, 162.55°E], spectrum-time diagrams of logarithmic likelihood function increment AlnL were obtained for six seismic stations IRIS: PET, YAK, OBN, MAG, YSS, and ARU that have similar characteristics and the location of which is shown in Figure 1. The stations are located in different seismic geological settings at considerable distances from each other. Station PET is located in subduction zone at a distance of A = 310 km from the epicenter. Station MAG second closest to the epicenter (A = 900 km) is located in
Time, h
Figure 11. Spectral time diagrams of the increment of the logarithmic likelihood function AlnL for the microseisms recorded at the PET and MAG stations. The ordinates show the spectral period. The arrows labeled as F, Fa, Fb mark the onset of foreshock activation and the two strongest foreshocks. The occurrence time of the Kronotskoe earthquake (M = 7.8) is shown by the thick arrow.
the north of the Sea of Okhotsk characterized by deep earthquakes. In the records of those two stations, periodic oscillations were revealed three hours before Kronotskoe earthquake (Figure 11). Red bands on the figure indicate emergence of periodic oscillations in the period range from 20 to 60 minutes. Increment AlnL was estimated for a sequence of time moments of seismogram maximums exceeding the level equal to the mean value in the window plus sample standard deviation in the same window. The count starts from 00 Greenwich time (UTC) 2 December 1997. Values AlnL were calculated in time window of duration AT = 3 hours with a shift AS = 1 hour, thus diagrams are presented in the figure, beginning at 3 o’clock 2 December. Last points in the time scale in diagrams (83 hours) correspond to 11 hours on 5 December (27 minutes before the moment of earthquake).
[44] The beginning of the period under investigation was chosen because 2 December 1997 was characterized by quiescent seismic conditions were noted in the zones where the enumerated-above stations are located and where local earthquakes of energy class K > 10 (M > 4) were not noted from Geophysical Service RAS data. Dramatic activation of the seismic process in the future Kronotskoe earthquake epicentral area started in Kamchatka in the middle of 3 December. That day three earthquakes of K > 10 and three earthquakes of K > 11 occurred there. The beginning of foreshock process is marked with an arrow and symbol F in the upper diagram of Figure 11. The activation went on with increasing number of shocks on the following day; seven events with K > 10, five events with K > 11 and one with K = 12.8 (M = 5.5) were recorded on 4 December. The latter is indicated with an arrow and symbol Fa in Figure 5. On
the day of Kronotskoe earthquake, 13 foreshocks K > 10, 12 foreshocks of K > 11 and 4 foreshocks of K > 12 including one of K = 12.5 (M = 5.3) were registered on 5 December before the earthquake moment. The latter is indicated with an arrow and symbol Fb in Figure 11. Altogether the foreshock series included 100 earthquakes of representative class K > 8.5 (M > 3).
[45] From Figure 11 it follows that the first series of periodic microseismic oscillations manifested itself at the end of 2 December - beginning 3 December. The second series started approximately at 2200 on 4 December (71 hours in diagrams) and was in progress until the time of the Kronotskoe earthquake main shock. Comparison was made with data of stations YAK, YSS, ARU, OBN, which are more remote from the source [Sobolev et al., 2005]. It was revealed that the first series of periodic oscillations was appeared in diagrams of AlnL of 2-3 December. In the second series, oscillations at those stations were noted immediately after foreshock Fa, but they were missing after foreshock Fb up to the earthquake moment. Since periodic oscillations 3 hours before the earthquake were only revealed at stations PET and MAG, which are the nearest to the epicenter, it is reasonable to assume their relation to the process in the sub-duction zone adjacent to the epicenter. The effect was most pronounced in PET station and maximum AlnL indicated a period of ~37 minutes 1 hour before the earthquake.
[46] Emergence of periodic oscillations in minute range of microseisms was revealed before Sumatra earthquake of
26 December 2004 with M = 9.2 and epicenter coordinates [3.32°N, 95.85°E]. Preliminary analysis of records made by stations in Figure 3 showed that station MBWA in Australia
M = 7.9 M= 9.2
Figure 12. Spectral time diagrams of the increment of the logarithmic likelihood function AlnL for the microseisms recorded at the KMI and CHTO stations. The ordinates show the spectral period. The occurrence time of the McQuary (M = 7.9) and Sumatra (M = 9.2) earthquakes are shown by the arrows.
had not been operated in the period of Sumatra earthquake, there had been failures and gaps in the records of DGAR and PALK stations and stations DAV and QIZ located in the Pacific region had a different structure of microseismic oscillations as compared to stations located in the Indian Ocean region. Therefore major studies were carried out from the data of stations CHTO, KMI, XAN, COCO, and partially PALK. We used records of vertical components with the exception of COCO station where this component had not been registered; in the latter case, horizontal component data were processed. The base of data of those stations that we used covered the interval from 6 to 26 December (00 hours б8 minutes) of that is to say until the moment of Sumatra earthquake.
[47] A singular feature was the fact that 2.б days before Sumatra earthquake in the southern hemisphere another strong earthquake had occurred of M = 7.9; the epicenter of the earthquake with coordinates [49.31°S, ^иЗб^] was located to the southwest of New Zealand (in Makkuori ridge area). Oscillations of this earthquake were hundreds of times as much as microseisms level in the above-mentioned stations.
[48] Time-frequency diagrams ДіпЬ shown in Figure 12 for stations KMI and CHTO were calculated with the use of method 2 described above. Arrows indicate the time when Sumatra earthquake (M = 9.2) and preceding Makkuori earthquake (M = 7.9) occurred. Periodic oscillations appeared after Makkuori and continued for one a day. Comparison with Figure ll suggests an effect similar to the one noted after the foreshock Fb of Kronotskoe earthquake. In the records of stations XAN, COCO and PALK periodic calculations were not revealed. Note that records of stations XAN and COCO were characterized with the noise increased level.
[49] In the last few days before Hokkaido earthquake of September 2003 with M = 8.3 and the epicenter coordinates [41.810N-143.910E] considerable foreshocks (M > б) or remote strong earthquakes were not noted. However the appearance of periodic oscillations was detected. We studied records of stations PET, YSS, OBN, ERM, MAJ, INC, MDJ, and BJT the location of which is shown in Figure 2. From calculations of parameter ДіпЬ it was revealed that at three stations PET, YSS and MDJ oscillations were noted 16 hours before the earthquake, which are presented in spectrum-time diagrams in Figure 13 (interval 4З00-б100 minutes). As distinct from Kronotskoe and Sumatra earthquakes, oscillations were revealed in more low-frequency range with periods of 120-160 min. One more group of oscillations was revealed fifty hours before the shock. Similar to Kronotskoe earthquake, periodic oscillations appeared in the records of stations close to the epicenter. Unfortunately ERM station located in the epicentral zone stopped recording 4 day before the earthquake.
Synchronization of Oscillations
[50] Methods 3 and 4 were used to search for oscillations synchronization effects at stations spaced apart.
Figure 13. Spectral time diagrams of the increment of the logarithmic likelihood function ДіпЬ for the microseisms recorded at the PET, YSS and MDJ stations. The ordinates show the spectral period. The occurrence time of the Hokkaido earthquake (M = 8.3) is shown by the arrow.
[51] In Figure 14, current values of Hurst generalized index а* are shown, which realizes the maximum of micro-seismic field singularity spectrum at б stations, the location of which is shown in Figure 1. We studied the interval of 30 days before Kronotskoe earthquake (309th-339th days in 1997). We used the window of 12-hour duration with a shift of one hour. With considerable variations of amplitude а* we failed to reliably separate intervals of synchronous oscillations at different stations. But coherence spectral measure calculation Л(т, w) as the modulus of canonical coherences product component by component (formula (9)) allows detecting them (Figure 1б). The major burst of coherence is centered in the vicinity of time markers of 40000-42000 minutes (several days before the earthquake). As the moving time window is approaching the earthquake moment, coherence level of а* -variations drops although it remains higher than background statistical fluctuations. Diagrams in Figure 1б testify to the increase of the time duration of low frequency “coherence spot” as the number of stations increases.
[52] To compare coherence level with different sets of stations their number should be constant according to formula (9). Sorting all possible combinations of 3 stations lead us to the conclusion that the set of stations PET, MAG, YAK have the highest value of а* = 0.6б and the distance from Kronotskoe earthquake epicenter to those stations is 3б0, 900 and 20б0 km; stations OBN, ARU and YAK have
M = 7.8
0 8000 16000 24000 32000 40000
Time, minutes from the beginning of 5 November 1997
Figure 14. Plots of the generalized Hurst exponent a* realizing the maximum of the singularity spectrum of the micro-seismic background at all stations with the estimation in a moving time window 12 h wide with a shift of 1 h. The coordinate of the right-hand end of the moving time window is plotted on the time axis.
the least value (0.32) and they are the farthest from the earthquake, located at 6800, 5900 and 2050 km. In all variants the coherence measure has a burst in the vicinity of the marker of 40000 minutes, which corresponds to observation time interval of 29.11-03.12.1997 that is 3-7 days before the shock. We note that on 03.12.1997 an intense series of foreshock activation started before Kronotskoe earthquake. Besides in this interval the increase of asymmetric impulses number was observed at PET station, which is the nearest to the earthquake epicenter (Figure 6). Just how random may be these coincidences is difficult to judge.
[53] In the analysis of the situation before Hokkaido earthquake, the largest number of stations participating in the computation was six: YSS, MDJ, INC, BJT, PET, OBN (Figure 2). Records of stations ERM and MAJ were not used because the former, as it was mentioned above, did not register microseisms in the last four days before the earthquake and the latter was not operated for 7 days two weeks before the earthquake. It was established that synchronization was manifested 2 days before the earthquake (interval 3300035000 minutes). It encompassed time periods from 3 hours (frequency 0.005 1 min-1) and longer ones. As the number of stations decreased, the amplitude A(t, w) increased in accordance with formula (9). It was significant that with complete sorting by 3 stations the most vivid effect was observed for stations nearest to the epicenter of Hokkaido earthquake. Spectrum-time diagram for such stations YSS,
MDJ, INC is shown in Figure 16. Three features may be noted: 1) synchronization with period ~3 hours (frequency ~0.005 1 min-1) started 9 days before the earthquake (23000 minutes); 2) most vividly and in a wide range of periods it was manifested 2 days before the earthquake (33000-35000 minutes); 3) a break in synchronization in the interval of 29000-31000 minutes was evidently associated with two remote strong earthquakes (shown with arrows) with magnitude 6.6. The first of them with the epicenter coordinates [19.72°N-95.46°E] occurred on 21 September and the second one with coordinates [21.16°N-71.67°W] occurred 10 hours later on 22 September. Of course, the arrival of seismic waves to stations at different times disturbed synchronization.
[54] In Figure 17, frequency-time diagram A(r, w) is given that was obtained from record processing of stations CHTO, KMI, XAN, COCO before Sumatra earthquake. Beginning with the time marker of 12800 minutes, coherence measure rise occurs with gradual lengthening of prevailing periods of oscillations from several minutes to tens of minutes. In the assumption of intra-terrestrial mechanism of the oscillations, they may be associated with resonance effects in the blocks increasing in scale and/or in the lithosphere layers and deeper layers of the earth. The analysis of microseism amplitude in the second range of periods showed that at all the above-mentioned stations in the interval of 16-26 December their level was practically stationary, which eliminates the atmosphere effects. Such phenomenon was previously noted in both the range of very long periods of the order of 1 year in seismic catalog studies and laboratory experiment of deforming and destructing a sample [Sobolev, 2003]. Apparently it is a fundamental characteristic of non-equilibrium system approaching instability.
Discussion
[55] We studied three types of effects: asymmetric impulses, oscillations periodicity, and noise synchronization in the range of periods from several minutes to tens of minutes from records of seismic stations of the same type in time intervals of the order of 1 month before 4 strong earthquakes.
[56] We show the results in Table 1 with the following
Table 1. Appearance of asymmetric impulses, oscillations periodicity and noise synchronization before the earthquakes
Earthquake I Iz P Pz S Sz
Neftegorskoe 27.05.1995, M = 7.0 + + - ? ?
Kronotskoe 05.12.1997, M = 7.8 + + + + -
Hokkaido 25.09.2003, M = 8.3 + + + + + +
Sumatra 26.12.2004, M = 9.2 - + ? + ?
Figure 15. Frequency-time diagrams of the evolution of the spectral measure of coherence of a* variation spectral series with the estimation in the moving time window 109 samples (5 days) long for a successively increasing number of simultaneously analyzed stations. Maximum values of the coherence measure are shown in each diagram after the codes of the stations analyzed.
symbols: in columns 2-7 I is for asymmetric impulses, P is periodicity, S is synchronization, Iz, Pz and Sz show location of the effect in the seismically active zone, where the earthquake occurred. Symbols (+), (-), (?) indicate appearance, non-appearance of the effect and uncertain result. Note that the table is based on our present knowledge based on limited amount of cases and relatively short time interval and may be improved in the future.
[57] Symbol (+) in columns Iz, Pz, Sz (location of the effect in the seismically active zone) means that the effect was more pronounced at stations that are located nearer to the earthquake epicenter. The uncertainty of results on synchronization before Neftegorskoe earthquake is explained by the fact that gaps in recording that did not coincide in time and lasted several days were left at most stations; interval
of 6 days suitable for joint analysis was too short to allow obtaining stable results. For Sumatra earthquake, symbol (?) is used in the columns of location in the seismically active zone, because now the data have not been analyzed that were obtained by stations located at distances considerably exceeding the rupture length (~1000 km) of this gigantic earthquake.
[58] We shall examine the feature of uniqueness of the effect appearance before the earthquake. If the effect appeared only once before the earthquake in the closing length of time interval under investigation, we shall conventionally consider it to be unique. Data having been obtained by the present do not allow us to conclude if the three types of effects mentioned above had the feature of uniqueness. It is established that the effects of oscillations periodicity
Right-hand end of the moving time window of the five-day length, minutes from the beginning of 1 September 2003
Figure 16. Frequency-time diagrams of the evolution of the spectral measure of coherence of a* variation spectral series with the estimation in the moving time window 109 samples (5 days) long for seismic stations YSS, MDJ, INC. Arrows indicate successively time moments of 2 remote earthquakes (M = 6.6) and of Hokkaido earthquake (the last one).
and noise synchronization may be repeated also after strong remote earthquakes and foreshocks. The possibility is not excluded that they might be repeated owing to trigger effect of tides and meteorological factors. The problem of uniqueness of asymmetric impulses appearance still remains to be solved.
[59] With account for the table we come to the conclusion
that the effects under investigation fall in the class of phenomena that are characteristic of non-equilibrium systems dynamics. The causes of their formation may be inside and outside the solid earth. Processes in the outer spheres of the earth (the atmosphere, ionosphere) are characterized by both random and quasi-periodic components. We proceed from the assumption that the dissipative system of seismi-
Figure
seismic
interval
17. Frequency-time diagrams of the evolution of the spectral measure of coherence A(t, w) for records of the stations XAN, KMI, CHTO, COCO after coming to 30-seconds sampling time Estimates were done within moving time window of the length 12 h with mutual shift 1 h.
cally active zone is in meta-stable state and processes going on in it have characteristics of determined chaos. Similar systems exist in different areas of the inner and outer spheres of the earth. If microseisms field in minute range of periods reflects space and time variations of different dynamic systems parameters and non-zero coefficients of relation between parameters of these systems exist, then mutual influence made by each system on another may be possible. It is well known that random systems show synchronization effects, especially in the attractors area. [Ott, 2002; Pykovski et al., 2003]. Synchronization of the systems dynamics may appear and be interrupted, and at some time intervals it may be stable [Gauthier and Bienfang, 1996].
[60] In applications, random systems are frequently encountered, in which the oscillations amplitude remaining finite changes in time irregularly from minimum to maximum and attractors are represented by cyclic orbits. [Rossler, 1976; Smirnov et al., 1997]. In such random systems, phase synchronization effects are manifested [Ott, 2002]. Characteristic curve of amplitude variation against time is shown in the upper part of Figure 18 (See also plots in Figures 5, 7, and 8). Let equation (10) describe random system that is affected by periodic oscillations.
dx/dt = F(x) + K * P(wt) (10)
[61] Suppose we deal with oscillations in the lithosphere and coefficient K shows the extent of influence of atmospheric pressure periodic disturbances made on them. Synchronization area in the frequency band w (Figure 18) is characterized by the following important characteristics [Ott, 2002]: it is not manifested if the relation coefficient K is less than threshold Ko; it expands as K increases. We can assume that as macro-instability (earthquake) approaches, the sensitivity of the lithosphere meta-stable area (value K) to the atmosphere pressure effect increases.
[62] In paper [Saltykov et al., 1997] facts were described of phase synchronization formation of high-frequency seismic noise (30 Hz) and tides. In our case, the rise of impulses in minute range of periods did not correlate with either the level of high-frequency microseisms of storm origin or the phases of earth tides. However, synchronization of stations records separated by thousands of kilometers and by tens of degrees in longitude suggests a common source. The rise of synchronization level when selecting stations located nearer to the earthquake epicenter may be associated with two features: the location of synchronization source in the corresponding seismically active zone and tensosensitivity rise of the zone to external effects caused by a remote source as the oncoming disaster approaches.
[63] Rhythms formation is a common phenomenon in nonequilibrium systems evolution [Nicolis and Prigogine, 1989]. It is important in this case that rhythmic oscillations may be of pulse form. Among impulses we revealed, there were symmetric and asymmetric impulses. By symmetry is meant approximately equal amplitude of positive and negative phases of oscillations. Figure 5 and curve 2 in Figure 7 show examples of such impulses. Asymmetric impulses locally appeared as series (Figure 9) and, as the moments of Kronotskoe and Neftegorskoe earthquakes approached, their number in-
dT = T(i+ 1) -7X0
Figure 18. Diagram illustrating the occurrence of the synchronization effects and pulses in a dynamic system containing chaotic and quasi-periodic components. The upper and lower curves are examples of temporal variations in the vibration amplitude (see explanations in the text).
creased (Figure 6). Such phenomena are known in the dynamics of random systems. The chart of impulse formation against oscillations of equation (10) type is presented in the lower plot (Figure 18). In paper [Gauthier and Bienfang, 1996] it is shown that with approaching macro instability in the stage of transitional bubbling conditions, synchronization intervals are interrupted by short-time bursts with large amplitude. With the development of macro instability the effect may be noted of bursts frequency increase [Ott, 2002].
P(dT) ~ (dT)-n , (11)
where P (dT) is probability of intervals formation of duration dT between bursts and n > 1, that is probability of short interval is higher. It was demonstrated in the studies of a relatively simple dynamic system when a particle placed in viscous liquid was affected by electric field and sinusoidal vibration. The situation inside the earth is much more complex, so this example only suggests that a phenomenon like this is possible in principle.
[64] Here it is worth noting one of the characteristics of impulse oscillations, which was only mentioned in the aforesaid. For the case of Kronotskoe, Neftegorskoe and Hokkaido earthquakes, the system of stations allowed us to compare records by stations located near the sources and by remote stations. For example, in the case of Kronotskoe earthquake, its epicenter was located at the following distances in kilometers from the stations participating in the studies: PET -350, MAG - 900, YSS - 1670, YAK - 2050, ARU - 5900, and OBN - 6800. If the impulse oscillations were elastic waves, then with duration (period) of T = 10 minutes the wavelength in the conditions of the upper lithosphere, the earth’s
crust would be from 2000 to 4000 km for surface and longitudinal body waves respectively. Then 2-3 stations closest to the epicenter would be within the wavelength (near zone) and we would expect that one and the same impulse would be recorded at those stations with a shift of several minutes. However attempts to identify clearly the same impulses at neighboring stations have not met with success; coincidence was observed only for strongest ones [Sobolev et al., 2005]. Therefore at the present stage we hold to the idea that the physical nature of the impulses is related to inelastic (quasi plastic) movements in the sources of earthquakes under preparation or fault zones located near them. In Kronotskoe and Hokkaido earthquakes it may have been subduction zone, in Neftegorskoe earthquake it may have been the regional fault extending along Sakhalin from the epicenter to YSS station. We draw attention to the fact that oscillations presented in Figures 5, 7 and 8 have dimensionality of displacement velocity (nM s-1). In case of one-polar impulse, they signify seismometer reaction to one step of ground displacement of several-minutes duration. The analogy of creep on the fault suggests itself.
Conclusions
[65] Microseismic field variations were studied in the minute range of periods from data of 20 wide-band IRIS stations before four earthquakes: Neftegorskoe with magnitude M = 7.0 in 1995; Kronotskoe earthquake (M = 7.8) in 1997; Hokkaido (M = 8.3) in 2003; and Sumatra (M = 9.2) in 2004.
[66] Several days before Kronotskoe, Neftegorskoe and Hokkaido earthquakes asymmetric impulse oscillations of several minutes duration were revealed that were broken by intervals of several tens of minutes. They were only manifested at stations located at the same seismically active zone where the earthquake occurred.
[67] Periodic oscillations were revealed before Kronotskoe earthquake with period of 20-60 minutes, Sumatra earthquake with period of 20-60 minutes and Hokkaido earthquake with period of 120-180 minutes. They appeared after foreshocks and remote strong earthquakes.
[68] Microseismic noise coherence was revealed at different stations before Kronotskoe earthquake with periods of more than 6 hours, Hokkaido with period more than 3 hours and Sumatra in the range of 2-60 minutes. Noise coherence measure increased if we chose stations located closer to the epicenter.
[69] These effects may appear and disappear several times in the process of earthquake preparation and fall in the category of phenomena characteristic of non-equilibrium systems dynamics.
[70] The nature of these effects is apparently associated with self-organization of seismic process in meta-stable lithosphere and synchronization of oscillations in the earth inner and outer spheres the dynamics of which comprise random and quasi-periodic components.
[71] To reveal their mechanism calls for further investigation with a greater length of observation series, an increased
amount of stations and greater number of earthquakes under investigation.
[72] Acknowledgments. The work was carried out with the support of the program “Electronic Earth” by the Presidium of the Russian Academy of Sciences, INTAS grant Ref. no. 05100008-7889 and RFBR grant no. 06-05-64625.
References
Bak, P., S. Tang, and K. Winsenfeld (1989), Earthquakes as selforganized critical phenomenon, J. Geophys. Res., 94, 15635.
Cox, D. R., and P. A. W. Lewis (1966), The Statistical
Analysis of Series of Events, 285 pp., Methuen, London.
Currenti, G., C. del Negro, V. Lapenna, and L. Telesca (2005), Multifractality in local geomagnetic field at Etna volcano, Sicily (southern Italy), Natural Hazards and Earth System Sciences, 5, 555.
Daubechies, I. (1992), Ten Lectures on Wavelets, no. 61 in
CBMS-NSF Series in Applied Mathematics, 449 pp., SIAM, Philadelphia.
Feder, J. (1989), Fractals, 247 pp., Plenum, New York, London.
Gauthier, D. J., and J. C. Bienfang (1996), Intermittent loss of synchronization in coupled chaotic oscillators: Towards a new criterion for high quality synchronization, Phys. Rev. Lett., 77, 1751, doi:10.1103/PhysRevLett.77.1751.
Hardle, W. (1989), Applied Nonparametric Regression,
333 pp., Cambridge University, Cambridge, New York, New Rochell, Melbourne, Sydney.
Kantelhardt, J. W., S. A. Zschiegner, E. Konscienly-Bunde, S. Havlin, A. Bunde, and H. E. Stanley (2002), Multifractal detrended fluctuation analysis of nonstationary time series, Physica A, 316, 87, doi:10.1016/S0378-4371(02)01383-3.
Kobayashi, N., and K. Nishida (1998), Continuous excitation of planetary free oscillations by atmospheric disturbances, Nature, 395, 357, doi:10.1038/26427.
Lyubushin, A. A. (1998), Analysis of canonical coherences in the problems of geophysical monitoring, Izv. Phys. Solid Earth, 34, 52.
Lyubushin, A. A. (2007), Geophysical and Ecological Monitoring Systems Data Analysis (in Russian), 228 pp., Nauka, Moscow.
Lyubushin, A. A., and G. A. Sobolev (2006), Multifractal measures of synchronization of microseismic oscillations in a minute range of periods, Izv. Phys. Solid Earth, 42(9), 734, doi:10.1134/S1069351306090035.
Lyubushin, A. A., V. Pisarenko, V. Ruzich, and V. Buddo (1998), A new method for identifying seismicity periodicities, Volcanology and Seismology, 20, 73.
Lyubushin, A., V. Pisarenko, M. Bolgov, andT. Rukavishnikova (2003), Study of general effects of rivers runoff Variations, Russian Meteorology and Hydrology, 7, 59.
Lyubushin, A. A., V. F. Pisarenko, M. V. Bolgov, M. Rodkin, and T. A. Rukavishnikova (2004), Synchronous variations in the Caspian Sea level from Coastal observations in 1977—1991, Atmospheric and Oceanic Physics, 40(6), 737.
Mallat, S. (1998), A Wavelet Tour of Signal Processing, 577 pp., Academic, San Diego, London, Boston, N.Y., Sydney, Tokyo, Toronto.
Marple, S. L., Jr. (1987), Digital Spectral Analysis With Applications, 492 pp., Prentice-Hall, Inc., Englewood Cliffs, New Jersey.
Nicolis, G., and I. Prigogine (1989), Exploring Complexity, An introduction, 328 pp., W. H. Freeman and Company, New York.
Ott, E. (2002), Chaos in Dynamic Systems, 478 pp., University, Cambridge.
Pykovsky, A., and al. et (2003), Synchronization, A Universal Concept in Nonlinear Science, 376 pp., University, Cambridge.
Rossler, O. E. (1976), An equation for continuous chaos, Phys. Lett., A 57, 397, doi:10.1016/0375-9601(76)90101-8.
Saltykov, V. A., and al. et (1997), Study of high-frequency seismic noise on the basis of regular observations on Kamchatka, Izv. Phys. Solid Earth (in Russian), 33(3), 39.
Smirnov, V. B. (1997), Experience of Estimating the Representativeness of Earthquake Catalog Data, Vulkanol. Seismol., (4), 93.
Sobolev, G. A. (2003), Evolution of periodic variations in the seismic intensity before strong earthquakes, Izv. Phys. Solid Earth, 39(11), 873.
Sobolev, G. A. (2004), Microseismic variations prior to a strong earthquake, Izv. Phys. Solid Earth, 40(6), 455.
Sobolev, G. A., A. A. Lyubushin, and N. A. Zakrzhevskaya(2005), Synchronization of microseismic variations within a minute range of periods, Izv. Phys. Solid Earth, 41(8), 599.
Sobolev, G. A., and A. A. Lyubushin (2006), Microseismic impulses as earthquake precursors, Izv. Phys. Solid Earth, 42(9), 721.
Sobolev, G. A., and A. A. Lyubushin (2007), Microseismic anomalies before the Sumatra earthquake of 26 December 2004, Izv.
Phys. Solid Earth, 43(5), 341, doi:10.1134/S1069351307050011.
Sornette, D., and C. G. Sammis (1995), Complex critical exponents from renormalization group theory of earthquakes: Implications for earthquake predictions, J. Phys. I. France, 5, 607, doi:10.1051/jp1:1995154.
Starovoit, O. E., and V. N. Mishatkin (2001), Seismic stations of Russian Academy of Sciences (in Russian), 86 pp., Geophysical Survey of Russian Academy of Sciences, Moscow-Obninsk.
Tanimoto, T., J. Um, K. Nishida, and N. Kobayashi (1998), Earth’s continuous oscillations observed on seismically quiet
days, Geophys. Res. Lett.,25(10),1553,doi:10.1029/98GL01223.
Telesca, L., G. Colangelo, and V. Lapenna (2005), Multifractal variability in geoelectrical signals and correlations with seismicity: a study case in southern Italy, Natural Hazards and Earth System Sciences, 5, 673.
G. A. Sobolev and A. A. Lyubushin, Schmidt Institute of the Physics of the Earth, Russian Academy of Sciences, Moscow, Russia
