CC BY
44
GEODYNAMICS & TECTONOPHYSICS
PUBLISHED BY THE INSTITUTE OF THE EARTH'S CRUST SIBERIAN BRANCH OF RUSSIAN ACADEMY OF SCIENCES
2016 VOLUME 7 ISSUE 1 PAGES 1-21
ISSN 2078-502X
http://dx.doi.org/10.5800/GT-2016-7-1 -0194
Experience of complexation of global geophysical observations
A. A. Lyubushin1, V. S. Bobrovskiy2, S. A. Shopin3
1 O.Yu. Schmidt Institute of Physics of the Earth of RAS, Moscow, Russia
2 Distant School 'Cosmic-Meteo-Tectonics', Petropavlovsk-Kamchatsky, Russia
3 Tula Geological Exploration Crew LLC, Tula, Russia
Abstract: The technique of joint analysis of heterogeneous time series of geophysical monitoring systems for the purpose of detecting time intervals and specific periods of bursts of synchronous behavior is presented. The technique is based on the use of the Fourier-aggregated signals and spectral measures of coherent behavior of multivariate time series, estimated in moving time windows. The article presents results of the analysis of data of underground electrical surveys at stations located in Kamchatka, Altai and Italy; the data were analysed together with torsion pendulum movements in Tula (Russia) and the time series of seismic noise parameters at the Japanese islands for the interval 2012-2015. The analysis identified a number of significant bursts of coherent behavior for these observations, some of which are presumably due to the strong mantle Okhotsk Sea earthquake of 24 May 2013. The coherent behavior of various geophysical fields before and after strong earthquakes is interpreted as a manifestation of the general pattern of increasing synchronization of fluctuations of complex systems at their approach to the rapid changes in their properties.
Key words: complex geophysical monitoring; multidimensional time series; coherence; aggregated signals; synchronization
Recommended by V.V. Ruzhich
For citation: Lyubushin A.A., Bobrovskiy V.S., Shopin S.A. 2016. Experience of complexation of global geophysical observations. Geodynamics & Tectonophysics 7 (1), 1-21. doi:10.5800/GT-2016-7-1-0194.
Опыт комплексирования глобальных геофизических
наблюдений
А. А. Любушин1, В. С. Бобровский2, С. А. Шопин3
1 Институт физики Земли им. О.Ю. Шмидта РАН, Москва, Россия
2 Дистантная школа 'Космо-Метео-Тектоника', Петропавловск-Камчатский, Россия
3 ООО «Тульская геологоразведочная партия», Тула, Россия
Аннотация: Излагается методика совместного анализа разнородных временных рядов систем геофизического мониторинга с целью выделения временных интервалов и характерных периодов всплесков их синхронного поведения. Методика основана на использовании Фурье-агрегированных сигналов и спектральных мер когерентного поведения многомерных временных рядов, оцениваемых в скользящих временных окнах. В качестве примера рассматриваются данные подземных электрических наблюдений на станциях, расположенных на Камчатке, Алтае и Италии, совместно с показаниями крутильных маятников в Туле и временными рядами изменения параметров сейсмического шума на Японских островах для интервала наблюдений 2012-2015 гг. В результате анализа выделен ряд значимых всплесков когерентного поведения полей в рассмотренном ряду наблюдений, часть из которых предположительно связана с сильнейшим мантийным
Охотоморским землетрясением 24.05.2013 г. Когерентное поведение различных геофизических полей до и после сильных землетрясений интерпретируется как проявление общей закономерности увеличения синхронизации флуктуаций сложных систем при их приближении к резким изменениям своих свойств.
Ключевые слова: комплексный геофизический мониторинг; многомерные временные ряды; когерентность; агрегированные сигналы; синхронизация
1. Introduction
One of the fundamental problems of geophysical monitoring is 'complexation' of different observations and measurements. This term means joint analysis of observations of different geophysical fields or observations of the same field, but at different measurement points, or both simultaneously. The idea of complexation is based on the hypothesis that using of a large number of monitored parameters can help extracting a weak common signal, which 'drowns' in the strong noise when each individual measurement is considered separately. The main feature of such a common signal must be its coherence (correlation) in a variety of observations, the use of which allows one to detect the very existence of the common components, despite the fact that the frequency content of the common signal may coincide with the frequency content of the strong local noise. Thus, the idea of complexation envisages joint analysis of geophysical characteristics, rather than their simultaneous measurements.
Identification of precursors of earthquakes or other geological catastrophes is among the most challenging problems of complexing measurements. In such a problem, a weak common signal is related to earthquake preparation processes, such as consolidation of the Earth's crust matter in a future earthquake focus and around it [Sobolev, Ponomarev, 2003]. The search for precursors of catastrophes as the occurrence of synchronous components in a variety of observations is the general idea about increasing the correlation radius of the random fluctuations of the parameters of a complex system as it approaches a sharp change in its properties as a result of its own dynamics [Gilmore, 1981; Nicolis, Prigogine, 1989].
The idea of complexation requires using methods of analysis of multivariate time series. In [Lyubushin, 2007], a set of algorithms for the analysis of multivari-ate time series monitoring systems is presented, which is an effective tool for discovering the hidden relations between processes, including those of different nature and structure. An important part of the developed algorithms is a preliminary analysis of time series with different scales with a purpose to extract the dimension-less characteristics that are independent in the physical
meaning of time series within successive adjacent time intervals of short length. Analysis of the noise is often neglected, although statistical regularities of the noise structure are an important source of hidden information about upcoming sharp changes in the properties of the objects under consideration. The methods are based on the analysis of canonical coherences, multidimensional spectral matrices and canonical correlation coefficients of the wavelet decomposition of signals in moving time windows, with all available samples (the so-called method of aggregated signals). The purpose of these algorithms is the detection of a very weak non-stationary signals of common origin, having both harmonic oscillation behavior or sharply non-stationary, wavelet character, within multivariate time series monitoring to identify their specific periods (or time scales).
This article presents an attempt of the joint complex analysis by applying software units to heterogeneous time series derived from several very distant from each other measuring points located in Kamchatka, Altai, Tula (a city in the middle of the European part of Russia), and Italy. Besides, time series of changing mean properties of seismic noise in the Japanese islands are used in the multivariate analysis. This sort of a global observation network covers a significant part of the Northern Hemisphere, and our analysis includes the data for the period from the beginning of 2012 to early May 2015.
An important result of this study is extracting the time interval of strong coherence before and after the powerful mantle earthquake (M=8.3) that occurred in the Okhotsk Sea on 24 May 2013.
2. Underground electrical measurements
The first group of data are results of underground electric measurements obtained from stations at different geographic locations. Technically, underground electric measurements can be traced back to classical telluric measurements being performed over a long period of time in seismic-hazardous areas (Japan, Greece, Kamchatka) [Varotsos et al., 1993; Uyeda et al., 2000; Lyubushin, Kopylova, 2004]. Potential differences
Fig. 1. a - measuring hole: 1, 2, 3, 4 - electrodes; b - measured potential differences. Stations automatically collect and transmit data to the central server: http://www.cosmetecor.org.
Рис. 1. а - измерительный шурф: 1, 2, 3, 4 - электроды; b - измеряемые разности потенциалов. Станции работают в автономном режиме и осуществляют передачу данных на центральный сервер, доступный в сети Интернет по адресу http://www.cosmetecor.org.
between electrodes are measured values for both methods. Using multielectrode systems in shallow holes at the tectonosphere-atmosphere boundary in case of underground electric measurements is the principal difference since telluric measurements are performed using long measuring lines (up to kilometers). Multi-electrode systems have been in use in studies of the Kamchatka region since 1989. The measurement technology was originally developed by Dr. D.A. Kuznetsov [Kuznetsov, 1991].
Typically, an underground electric measurement site has three shallow holes, south-west (SW), central (C) and north-east (NE), spaced by 5-10 m and 1.5 to 3.0 m deep (Fig. 1, a). Iron electrodes (500x500x3 mm) are placed horizontally in the hole and separated from each other by a tamped soil layer with thickness of ~300-500 mm. Each electrode is connected to a coaxial cable. The electrode-cable connection is protected from corrosion by the resin sealing. Typically, each hole has four measurement electrodes, but more electrodes (six, eight, etc.) can also be used.
Measured values are potential differences between electrodes in a single hole (electrode-electrode scheme), between electrodes in two different holes (subhorizontal scheme) and between each electrode and the local electrical ground (GND). A general scheme of measurements is shown in Fig. 1, b.
Arrowed lines show measured potential differences. Electrodes in holes are marked by increasing sequential numbers from uppermost to lowermost. A power main neutral wire or a local earthing system or a steel water pipeline, having electrical contact with ground, are used as the local electrical ground.
The full set of measured parameters for the four-electrode measurement scheme includes:
• Twelve underground electromotive forces between each electrode and GND in the GND scheme (shown by solid arrows between each electrode and GND);
• Nine underground electromotive forces in the electrode-electrode scheme (shown by dotted line arrows between each electrode for every hole);
• Eight underground electromotive forces in the subhorizontal scheme (shown by dotted line arrows between the second electrode in the central hole and each electrode in the NE and SW holes). Therefore, 29 underground electromotive forces
(EMF) are measured.
Currently, measurements are performed in the frequency band of 0...4 kHz. For each channel, AC and DC components of underground EMF are measured. The AC component is measured as a half-period average voltage.
Measurements are performed using standard external USB DAQ module E14-140M by L-Card LLC (http://www.lcard.ru), which is certified and approved for measuring usage.
Electrode signals are digitized with the sample rate of 12.5 or 14 kHz. Calculation of AC and DC components is performed by the software. The signal processing scheme is presented in Fig. 2.
The DC component is calculated by the low-pass filtering of input signal (cut-off frequency is 10 mHz) followed by the decimation to 1 Hz sampling rate.
The AC component is calculated by the high-pass filtering of input signal (cut-off frequency is 4 Hz), followed by rectification, low-pass filtering (cut-off
Fig. 2. Scheme of processing signals from electrodes. Рис. 2. Схема обработки сигналов электродов.
frequency is 10 mHz) and decimation to 1 Hz sampling frequency. All used digital filters are designed as the 6th order IIR Bessel filters and implemented as cascaded biquadratic structures. The above-described processing scheme is similar to the measurement of AC and DC components by the regular multimeter. It is used to provide for comparability of new data with the data sets obtained prior to 2012 by manual multimeter measurements.
The current data acquisition network consists of nine stations, which parameters are shown in Table [Bobrovskiy, Kuznetsov, 2011; Bobrovskiy, 2011] (in the table, abbreviation P-K means Petropavlovsk-Kamchat-sky). Stations S1_PK andS1_UZ are parts of the single station that is technically subdivided in two parts because of the constraint on the input channel count of the measuring equipment. In Table 1, only Station S1_PK + S1_UZ has all 29 measuring channels shown in Fig. 1, b. Channel count for other stations is decreased with respect to different considerations.
Each station data is the multichannel record, containing potential differences between pairs of electrodes, recessed at different depths and in different holes. Channel count for stations varies between 26 and 32.
EMF measuring stations
Измерительные станции
Station name Location EMF count
S1_PK P-K 16
S1_UZ P-K 13
S1_IMFSET P-K 16
S2 P-K 14
S3 P-K 16
S7 P-K 16
S4 Esso, Kamchatka 14
S5 Gorno-Altaisk 16
S6 Chiety, Italy 16
S8 Crimea 16
Time series of stations contain the values of AC and DC components sampled at one-second time interval. The number of time series for each station equals the doubled number of station EMFs. For instance, time series of Station S1_PK has 32 values of EMFs. First 16 of them are DC components, and last 16 are AC components of corresponding EMFs.
The initial records present signals sampled at one-second time interval. Such time step is too small for the analysis of synchronization and coherence effects between both different channels and records at different stations. That is why a transfer to one-hour sampling time step is performed by computing mean values within adjacent time intervals for one hour. Then, increments of the signals are calculated to suppress dominated low-frequency components and implement the winsorization operation [Huber, 1981] which is an iteration procedure of computing mean values and standard deviations a, subtracting mean values and dividing by a, and clipping values which exceed thresholds ±4c. These iterations are repeated for each channel until the values of a will not be changed (usually 10 iterations are enough). Winsorization is a simple tool which provides robustness of estimates with respect to outliers within increments of the signals [Huber, 1981].
Fig. 3 presents graphics of all 32 channels of increments after transition to one-hour step for Station S1_PK for the time period from 10 August 2012 to 10 May 2015.
3. THE DATA OF TORSIONAL PENDULUM SYSTEMS
The second group of the source data includes measurement results of special pendulum systems similar to the asymmetric horizontal torsional pendulum [Mar-tynov, 2008; Martynova, Martynov, 2006]. These measurements are performed in Tula State University (Tula, Russia).
40 0 -40
40 0 -40
8 4 0 -4 -8 4 2 0 -2 -4 -6
4
0 -4
2 1 0 -1 -2 4 2 0 -2 -4 4 2 0 -2 -4
02
IliMllini
■UP
03
04
05
UNWHf
06
ФШ
07
08
РШНШ
10000 20000 30000
10000 20000 30000
10000 20000 30000
25
26
&
27
28
29
30
31
32
10000 20000 30000
Fig. 3. Graphs of time series increments after coming to the one-hour sampling time step at each of 32 channels of Station S1_PK from 10 August 2012 to 10 May 2015. Time marks at the X-axis correspond to number of hours since the beginning of 2012.
Рис. 3. Графики приращений временных рядов показаний на каждом из 32 каналов станции S1_PK после перехода к шагу по времени 1 час для временного фрагмента 10.08.2012 г. - 10.05.2015 г. На оси абсцисс отложены временные метки в часах от начала 2012 г.
One variant of the construction of the torsion system is shown in Fig. 4, a. A 15-20 m long beam is suspended with a thread at the center of gravity. A tungsten wire (diameter of 50-100 [im; length of 0.4-1.5 m) is used as a suspension thread. Two weights are fastened at the beam's ends. Their typical mass values are about 10 g. One of the weights has a complex shape
providing asymmetry for the whole torsion system. Elements of the torsion system are made from nonmagnetic materials. Thread suspension point O is stationary.
The measured value is the twist angle of the thread of the torsion system, which is the same as the angle of rotation of the beam. Every instrument has several
Fig. 4. a - torsion system of the instrument: O - suspension point of the torsion system; 1 - suspension thread; 2 - beam; 3 - weight having a complex shape; 4 - weight-counterbalance; b - instrument case construction: 1 - working volume of the case; 2 - subassembly for fastening and adjustment of the torsion system; 3 - instrument base; 4 - cover of access hole; 5 - duct; 6 - duct support; c - method of rotation angle measurement: 1 - mirror; 2 - light source; 3 - photodiode; 4 - beam of light; 5 - suspension thread.
Рис. 4. а - крутильная система прибора: O - точка подвеса крутильной системы; 1 - нить подвеса; 2 - коромысло; 3 - груз сложной формы; 4 - груз-противовес; b - конструкция корпуса прибора: 1 - рабочий объем корпуса; 2 - узел крепления и регулировки крутильной системы; 3 - основание прибора; 4 - крышка рабочего монтажного окна; 5 - штанга; 6 - опора штанги; c - способ измерения угла поворота: 1 - зеркальце; 2 - источник света; 3 - фотодиод; 4 - луч света; 5 - нить подвеса.
torsion systems placed inside the grounded metal screen case (Fig. 4, b). Every torsion system is associated with a separate measurement channel of the instrument. Beams with weights are located inside the working volume of the case made from thick steel (thickness is about 20 mm). Torsion systems suspension threads are located inside the ducts. Fastening of suspension threads and setting of torsion system equilibrium positions is performed using the subassembly.
Every torsion system has the optoelectronic system (Fig. 4, c) measuring the beam angle of rotation and transmitting the acquired information to the PC. Measurement is performed by the registration of the pho-tocurrent of the photodiode illuminated by the light spot reflected by the mirror that is rigidly connected to the beam. A light emitting diode or a laser diode is used as the light source. Components of angle-of-rotation sensors are located in the duct supports (Fig. 4, b).
Measurements are taken in the automatic round-the-clock mode. Since 1993, a significant data base has been accumulated. The construction of the instruments and the method of angle of rotation measurement are described in more details in [Martynov et al., 2006a, 2006b; Shopin, 2009, 2014, 2015]. Investigation of similar torsion systems are performed by I. Kalinnikov and
colleagues in the Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences [Kalinnikov et al., 2011].
Time series of WBG-3, a three-channel instrument, are used in the present work. These are series of angular positions of three torsion systems of instrument in arbitrary units. The torsion systems have the same construction, but different suspension thread lengths and spatial orientation. The time series are sampled with a 0.5 Hz sample rate. Graphs of instrument data are shown in Fig. 5 after transition to the one-hour time step by calculating mean values in the adjacent one-hour time windows.
4. Seismic noise parameters
In ours study, continuous seismic records from F-net broadband seismic network in Japan are used as a source of seismic data. This network was established in 1997, and its data are available for free downloading at http://www.fnet.bosai.go.jp/faq/?LANG=en.
On 11 March 2011, a mega-earthquake with magnitude 9 occurred in Japan. This earthquake is remarkable by the fact that it was predicted in advance. This
Time, hours from the beginning of 2012
Fig. 5. Readings from three torsional pendulums of differ-rent orientations (conventional units) from the beginning of 2012 to 11 March 2015 after coming to the one-hour sampling time step.
Рис. 5. Показания 3-х крутильных маятников (условные единицы) различной ориентации за интервал времени с начала 2012 г. по 11.03.2015 г. после перехода к шагу по времени 1 час.
prediction was based on the analysis of seismic noise properties from F-net . The hypothesis of a great seismic catastrophe approaching Japan was formulated at the middle of 2008. Later on, as the new data were coming from the network, the estimates of the future event became more precise. In April 2010, an estimate was made that the middle of 2010 could be regarded as the beginning of waiting time for earthquake with magnitude 8.5-9.0. This prediction was voiced at international conferences before the seismic catastrophe and published in several papers and abstracts [Lyu-bushin, 2009, 2010a, 2010b, 2010c, 2011a], The paper [Lyubushin, 2011a] was published after the event but it was submitted in April 2010.
After the earthquake, the experience of its prediction was published in a number of scientific papers, and it was shown that the available software technique provides for estimating a place of a future event as well, but it was done retrospectively [Lyubushin, 2011b, 2012, 2013, 2014a],
Positions of F-net stations are shown in Fig. 6. The vertical seismic records with sampling rate 1 Hz (LHZ
records) were downloaded from the site and transformed to the one-minute time step by calculating mean values within adjacent time intervals for 60 samples.
Seismic noise waveforms statistics for further use is described below.
Logarithm of kurtosis, lg(k). Kurtosis k is defined by the formula [Cramer, 1999]:
k = ((Ax)4)/((Ax)2)2 -3, (1)
where Ax is deflection of seismic noise waveform from its trend, (...) is the symbol of sample estimate of mean value. Kurtosis characterizes the sharpness of probability distribution forms and gives a measure of deflection of Ax from normal distribution for which k=0. The values in equation (1) are computed within adjacent time intervals (without overlapping) of the length 1440 neighbor one-minute samples (one day). Before calculating (1) within each daily time window, the trend is removed by orthogonal polynomial of the 8-th order. Removing the trend eliminates deterministic trends which are caused by the influence of tidal and temperature deformations of the Earth's crust and impact the investigation of the noise properties. For seismic noise, kurtosis is positive, and usually k=1 - that is why we will consider its logarithm, lg(k). Kurtosis was used for investigating properties of seismic noise in [Lyubushin, 2014b, 2015].
Minimum normalized entropy En of squared wavelet coefficients. Let x(t) be the finite sample of the signal t=1,..., N - index, numerating the counts. The normalized entropy is defined by the formula:
N
V- log(Pfc)
N
Pk = <En< 1, (2)
7 = 1
where ck,k=1,N are the orthogonal wavelet coefficients calculated from the minimized value (2). We try 17 orthogonal wavelets [Mallat, 1998]: 10 usual wavelets of Daubechies (number of vanishing moments equals to integer numbers from 1 up to 10) and 7 so-called sym-lets with numbers of vanishing moments varying from 4 up to 10. For low-frequency noise, parameters En were estimated within adjacent time windows of length N=1440, i.e. one day, after removing the trend by polynomial of the 8-th order. Minimum normalized entropy En was suggested in [Lyubushin, 2012] and used for investigating seismic noise properties in [Lyubushin, 2013,2014b],
Multifractal parameters Aa and a*. Multifractal singularity spectrum F(a) [Feder, 1988] of signal x(t) is defined as a fractal dimensionality of time moments
128 132 136 140 144 148
E, deg
I Fig. 6. Positions of 78 broadband seismic stations in Japan islands.
I Рис. 6. Положение 78 широкополосных сейсмических станций на Японских островах.
All details of the multi-fractal properties computation methods used here are described in [Lyubushin, 2007, 2009, 2010b, 2011b], Multi-fractal characteristics of seismic noise are widely used to detect precursors of earthquakes and dynamic estimates of seismic danger, and the results of such studies were presented in [Lyubushin, 2009, 2010a, 2010b, 2010c, 2011a, 2011b, 2012, 2013, 2014a, 2014b, 2015].
Another parameter of seismic noise is logarithm of variance, lg(Var).
Values lg(k), En, Aa, a*and lg(Var) are performed in the adjacent time windows for 1440 samples with the one-minute time step which correspond to one day. Before computing within each time window, the trend is removed by polynomial of the 8-th order. Thus, for each station, we obtain time series of values of lg(k), En, Aa, a*and lg(Var) with the one-day time step. For each day, we calculate median values of the seismic noise parameters by values from all stations of the network. Daily median values for the time period from the beginning of 2012 to 31 July 2015 are shown in Fig. 7.
5. Methods of analysis: multiple
COHERENCE COEFFICIENT
ta which have the same value of local Lipschitz-Holder exponent: h(t) = lim5^0(ln(^(t, 5))/ln(5)), i.e. h(ta) = a, where ^(t, 5) = (max x(s) — minxes)), maximum and minimum values are taken for argument t — 5/2 < s<t + 5/2. Value ^(t, 5) is a measure of signal variability in the vicinity of time moment t [Feder, 1988]. If X(t) is a usual self-similar mono-fractal signal [Taqqu, 1988] with Hurst exponent value 0< H < 1, then F(H) = 1, F(a) = 0Va ^H, but the finite sample estimate of singularity spectrum does not obey these rigorous theoretical conditions.
Practically, the most convenient method for estimating singularity spectrum is a method of multifractal detrended fluctuations analysis (DFA) [Kantelhardt et al., 2002] which is used here. Function F(a) can be characterized by the following parameters: amin, amax , Aa = amax -amin, and a* - an argument providing maximum to singularity spectra: F(a*) = = maxa F(a). Parameter a* is called a generalized Hurst exponent, and it gives the most typical value of Lipschitz-Holder exponent. Parameter Aa, singularity spectrum support width, can be regarded as a measure of variety of stochastic behavior. For removing scale-dependent trends (which are mostly caused by tidal variations) in multifractal DFA-method of singularity spectrums, estimates a local polynomials of the 8-th order are used.
An ordinary spectrum of the coherence of two processes can be non-rigorously defined as the square of the correlation coefficient for these processes at the frequency <d. The canonical coherences represent a generalization of the concept of the coherence spectrum to the situation when, instead of a pair of scalar time series, it is necessary to investigate the relation between two vector time series: the m-dimensional series X(t) and the n-dimensional series Y(t), at different frequencies. Without loss of generality, we consider that n<m
Quantity p2 (w) called a square of the modulus of the first (maximum) canonical coherence of the series X(t)and Y(t), which replaces an ordinary coherence spectrum in this case, is calculated as the maximum eigenvalue of the following matrix:
t/(w) Sxx SXySyy SyXSxx , (3)
where ^ is frequency, Sxx(w) - spectral matrix of time series X(t) of size m*m, Sxy(w) - cross-spectral matrix of e size mxn, Syx(ui) = Sxy(w), H is the sign of Hermi-tian conjugation. Quantityp2(w) replaces the square of the modulus of the coherence spectrum in case of two multidimensional signals.
We now introduce the concept of component-by-component canonical coherences of q-dimen-
sional time series Z(t) as squares of the moduli of the
3
2
1
0
-1 1
0.8 0.6 0.4 0.2
0
3
2.5 2 1.5 1
0.5
0 0.6
0.5
0.4
0.3
0.2
0.1 8
6
4 2 0
igM
En
Да
or
lg(Var]
ll
liu
2012
2013
2014
2015
I Fig. 7. Graphs of medians of five daily low-frequency seismic noise properties calculated from the data obtained by F-net seismic network in Japan from the beginning of 2012 to 31 July 2015.
I Рис. 7. Графики медиан ежесуточных значений 5 свойств низкочастотного сейсмического шума, вычисленных по данным сейсмической сети F-net в Японии с начала 2012 г. по 31.07.2015. г.
maximum canonical coherence in the situation when, in formula (1), the i-th scalar component of q-dimen-sional series Z(t) is assumed to be series Y(t), and (q-1)-dimensional series consisting of other components is assumed to be X(t) series. Therefore, quantity characterizes the coherence of variations in the i-th component with variations in the set of all the other components at frequency w.
By introducing the component-by-component canonical coherence, it becomes possible to determine
one more frequency-dependent statistics k(a>) [Lyu-bushin, 1999, 2007], which characterizes the coherence of variations in all the components of vector series Z(t) at frequency w:
k(a)) = ^Vi(a)), (4)
¿=i
where q is a dimensionality (number of scalar components) of time series Z(t) in our study.
Note that, by virtue of the construction, the value of belongs to interval [0.1], and the closer is the corresponding value to one, the stronger is the relation between variations in the components of multidimensional time series Z(t) at frequency w. Frequency-dependent quantity (2) can be called the spectral measure of the coherence of a multidimensional time series or multiple coherence coefficient. It should be noted that values of can be compared for the same values of dimensionality q only, because value (4) is defined as the product of quantities which are less than 1. If q=2, measure is a usual squared coherence
spectrum.
Thus, estimating the multiple coherence of multidimensional time series is the problem of calculating statistics and investigating its peculiarities.
To estimate the temporal variability of the interaction of recorded processes, it is necessary to perform calculations in the moving time window of a specified length. Let r be the time coordinate of the window having a length of L counts. Calculating the spectral matrices for the samples falling in time window t, we obtain the two-parameter function, k(r, w). The bursts of the k(r, w) value will determine the frequency bands and time intervals in which the collective behavior of jointly analyzed processes is enhanced.
To realize this algorithm, it is necessary to have the estimate of spectral matrix Szz(j,m) with size qxq in each time window. Below, we prefer to use the model of vector autoregression [Marple, 1987]. The method
consists in the estimation of model parameters:
p
Z(t) + YjAk-Z(t-k)=eit), (5)
k=1
where Ak are matrices of autoregression parameters of size qxq, p is autoregression order, e(t) is a q-dimen-sional time series of identification residuals which is assumed to be the sequence of independent Gaussian vectors with a zero mean and an unknown covariance matrix O = M{e(t)eT(t)} which is considered independent on the time index. It is important to note that model (5) was constructed after the preliminary operations of eliminating the general linear trend, winsoriza-tion and normalizing each scalar component to the unit variance. These operations were performed independently in each time window of processing and for each scalar component of the multidimensional series. Their meaning consists in eliminating the influence of diversification in scale in the series processed and provides spectral matrix estimate robustness.
To estimate the matrices Ak and O, the Durbin-Levinson recurrence procedure [Marple, 1987] is used, for which the sampling estimates of the covariance matrices must be preliminarily calculated.
The estimate of spectral matrix is calculated by the formula:
Szz(^) = F-10)-O-F-"(^), (6)
where F(^) is a complex matrix:
p
F(o)) = E+ -exp(-i<j)k), (7)
k=1
E is a unit matrix of size qxq.
Estimate (6) has a rather high resolution in frequency for short samples and, therefore, is more preferable for estimations in a moving time window than, for example, nonparametric estimates in terms of the averaging of multidimensional periodograms. There are no formalized procedures for choosing autoregression order p. In the calculations, p was chosen by the trial method as the minimum value, which further increase does not lead to any substantial change in the main elements of the behavior of k(r, w) dependence.
6. Methods of analysis: Fourier-aggregated signal
An aggregated signal is constructed in two stages. At the 1st stage, initial multidimensional time series is substituted by time series of the same dimensionality but composed of so-called canonical components. The canonical components accumulate signals that are common for all initial components and are free of local variations that are specific for individual scalar time series only. At the 2nd stage, the common signals are amplified by constructing a single scalar time series, their first principal component. Thus, an aggregated signal is the first principal component of canonical components [Lyubushin, 1999,2007].
An aggregated signal is intended for extracting common harmonic stationary oscillations and is based on the analysis not in moving time windows of rather short length but on using information from the whole length of available sample or from its given part. Let's consider q -dimensional time series Z(t).
We will extract the i-th scalar component Z(t) and filter (q-1) - dimensional series Z(i)(t) consisting of the all other components in such a way that scalar signal Cf(t) obtained at the filter output has the maximum coherence with extracted series Z(t) at each frequency. In order to do this, components of the eigenvector of matrix (3), where Z(t) appears as Y(t) and Z(i)(t) appears as X(t), corresponding to its maximum eigenvalue (obviously equal to should be used as a multidimensional frequency filter [Lyubushin, 1998a, 1998b, 2007]. If component Z(t) contains the noise that is characteristic solely for this series and is absent in
the other components of the series Z(t), i.e. in Z(i)(t), it is absent in signal Cf(t) simply by its construction, and this is the meaning of such an operation. At the same time, series Cf(t) retains all the components of Z,(t) which are common for the other components of series Z(t), i.e., for signal ZW(t).
We now determine aggregated signal AZ(t) of multidimensional time series Z(t) as the first principal spectral component of multidimensional series C(Z)(t) composed of canonical components Cf(t) of each scalar time series forming initial series Z(t). Recall that the main spectral component is the projection of the vector of Fourier transforms on the eigenvector of the spectral matrix corresponding to its maximum value [Brillinger, 1975; Hannan, 1970]. It should be emphasized that series AZ(t) differs from the simple first principal component. In either of the cases, the series are determined by multidimensional filtering, in which the eigenvectors of the spectral matrices corresponding to their maximum eigenvalues are assumed to be multidimensional frequency filters. However, for the ordinary first principal component, the matrix of initial time series Z(t) is such a spectral matrix, whereas the spectral matrix of series C(Z)(t) is the matrix of series AZ(t). Although the common components are extracted in the course of either filtering, aggregated signal AZ(t) is preferable, because individual noises are completely eliminated in the course of its construction.
Contrary to the estimation in a moving time window, nonparametric estimation by the frequency averaging of periodograms and cross-periodograms [Brillinger, 1975; Hannan, 1970; Marple, 1987] was used to estimate the spectral matrix required for the construction of the aggregated signal. Such a choice was related to a higher structural stability of the classical periodo-gram estimates of power spectra for long time series compared to parametric autoregression estimates of spectral matrices (6), which are advantageous for short samples. We used a deep averaging (smoothing) of pe-riodograms in the frequency window with the length equal to 1/32 part of the total number of discrete frequency values.
Before computing spectral matrixes, the data are processed by a number of operations, including linear trend removing, winsorization and smoothing at the end parts of time interval by cosine window (the so-called tapering operation) [Brillinger, 1975; Hannan, 1970]. The length of the end parts is taken 1/8 share of the whole length. These preliminary operations are fulfilled before computing spectral matrix only. Further on, for computing canonical components Cf(t), data are transformed to frequency image by the fast discrete Fourier transform, but without winsorization and cosine tapering, and the resulting multidimensional Fourier transforms are projected to spectral matrix eigenvectors of the matrix (3). The same procedure is re-
peated for computing aggregated signal: at first, computing spectral matrix of canonical components, then determining eigenvector corresponding maximum eigenvalue of this spectral matrix and projection of frequency Fourier images of canonical components on these vectors. All these operations are done for each frequency value. Final inverse Fourier transfer of the results of such projections provides temporal realization of the aggregated signal.
Note that the aggregated signal has no physical dimensionality. Since it is constructed after the sequence of operations aimed at the normalization of the initial data, its meaning consists solely in the formal extraction of the most common harmonic variations. All the elements of the computational technology are described in detail in [Lyubushin, 1999,2007].
7. Analysis of underground electrical measurements
Before the analysis of the electrical data, channel selection was done at each station. For this purpose, tables of pairwise correlation coefficients were computed, and in pairs of the channels which have absolute values of correlations exceeding 0.9, one channel was excluded from the analysis. The number of excluded channels may be very significant because many channels are strongly correlated with each other, and the correlation coefficient amounts to 0.999-1 (it is noticeable even from a purely visual analysis in Fig. 3). For example, for Station S1_PK 14, channels 1, 3, 12, 14, 15, 19, 20, 21, 22, 23, 24, 25, 27, and 28excluded from the analysis.
Then, for each station, an aggregated signal was computed for the channels remaining after selection. Such signals for the electrical measurements at the stations are called stations' aggregated signals (Fig. 8).
After building the stations' aggregated signals, it becomes possible to do the aggregation of the 2nd stage, i.e. calculation of aggregate signals from different sets of stations' aggregated signals. Because among 10 measuring stations , seven are located in the Kamchatka Peninsula, it is reasonable to build the 2nd order aggregated signals for various numbers of stations in Kamchatka which have simultaneous measurements. Figure 9 shows results of the 2nd order aggregation for combinations of different stations in Kamchatka.
Figure 10 illustrates the time-frequency structure of the 2nd order aggregated signals. It is evident that the increase in the number of stations has an insignificant impact on the 2nd order aggregated signal spectral structure. When choosing the number of stations included in the 2nd order aggregation, a compromise is necessary: the number of stations involved in the aggregation should be large enough to suppress the noise and to identify patterns and, at the same time, if the
S1 PK, P-K
S1_UZ P-K
S1JMFSET P-K
S2, P-K
S4, Kamchatka |
i iilltiiii ° iiihi rillli ......1w ...... ppfpfipf
Time, hours from the beginning of 2012
5000
10000
15000
20000
25000
I Fig. 8. Graphs of aggregated signals from each electrical measurement station, which are built after selecting the channels. Start time moments are different whereas the end time moments are the same (10 May 2015).
I Рис. 8. Графики агрегированных сигналов от каждой станции электрических измерений, построенные после селекции каналов. Начальные моменты времени у сигналов различаются, но все они заканчиваются 10.05.2015 г.
Geodynamics & Tectonophysics 2016 Volume 7 Issue 1 Pages 1-21 (a) 3 stations: S1_PK, S2, S4
(b) 4 stations: S1_PK, S1_UZ, S2, S4
(C) 5 stations: S1_PK, S1_UZ, S2, S4, S7
(d) 6 stations: S1_PK, S1_UZ, S2, S3, S4, S7
(e) 7 stations: S1_PK, S1_UZ, S1JMFSET, S2, S3, S4, S7
Time, hours from the beginning of 2012
5000
10000
15000
20000
25000
I Fig. 9. Graphs of the 2nd order aggregated signals from different number of electrical measurements stations (a-e) in Kamchatka.
I Рис. 9. Графики агрегированных сигналов 2-го порядка от различного числа станций электрических измерений (a-e) на Камчатке.
number stations is too large, the time interval of joint measurements is too short. From sorting options, it appears that the most appropriate is the 2nd order aggregated signal for four stations (Fig. 9, b).
Figure 11 shows the time-frequency diagram of the evolution of the coefficients of coherence: in Figure 11, a, - pairwise coherence quadratic coefficient, and Figure 11, b, - 11, g, - the coefficients of multiple coherence between the aggregated signals from stations in Altai and Italy and the 2nd order aggregated signals from different combinations of stations in Kamchatka. Estimation was done in a moving time window of the length 672 hours or 28 days (lunar month) with mutual shift of windows 24 hours, using the model of autoregression of the 7th order.
Figure 11 implies the existence of an essential burst of coherence for oscillations with periods of about 24 hours in a time fragment of 12000-16000 hours from the beginning of 2012, which corresponds to the posi-
tions of time windows of the length 28 days from mid-April to the end of October of 2013. Time mark 15000 hours from the beginning of 2012 corresponds to the position of the right-hand end of the time window to 16 September 2013.
8. Joint analysis of electrical observations data
AND DATA OF TORSION PENDULUMS
In order to include the data of torsional pendulums in Tula, we must build their aggregated signal. Figure 12 shows a graph of the aggregated signal of increments of time series of the torsion pendulums.
Information from the torsional pendulums was included into analysis of the coherent behavior of aggregated signals from four different geographical points, Altai, Kamchatka, Italy, and Tula. For this purpose, the time-frequency diagram for the multiple coherence
B000 10000 12000 14000 10000 13000 20000 22000 24000 26000 28000
Time, hours from the beginning of 2012, right-hand end of moving time window of the length 672 hours with mutual shift 24 hours.
10000 20000 22000 24000 26000 Time, hours from the beginning of 2012, right-hand end of moving time window of the length 672 hours with mutual shift 24 hours.
Fig. 10. a' and b' - time-frequency diagrams of power spectra logarithms of the 2nd order aggregated signals for three and six stations in Kamchatka (Fig. 9, a, and 9, d); a" and b" - graphs of correspondent averaged spectra from all time windows.
Рис. 10. а' и Ь' - частотно-временные диаграммы эволюции логарифмов спектров мощности агрегированных сигналов 2-го порядка для трех и для шести станций на Камчатке (рис. 9, а, и 9, а'' и Ь" - графики соответствующих усредненных спектров от всех окон.
coefficient was computed (Fig. 13). As for previous cases, the time window of the length of 672 hours was used, which was taken with mutual shift 24 hours for the autoregression model of the 7th order.
Figure 13 (compare with Figure 11) shows again a burst of coherence in the range of time marks of the right-hand end of the moving time window of the length of 28 days for the time marks from 12000 up to 16000 hours from the beginning of 2012, i.e. for the time interval from mid-April to the end of October 2013.
9. Joint analysis with variations of seismic
NOISE PARAMETERS
One of the hypotheses to explain the occurrence of a burst of coherence in the time-frequency diagram in Figures 11 and 13 for the time windows within the time period of April - October 2013 (in the vicinity of
the time mark of 15000 hours from the beginning of 2012), is the global synchronization of the geophysical fields shortly before and shortly after the largest mantle earthquake (M=8.3) in the Okhotsk Sea on 24 May 2013 [Chebrov et al., 2013], the response to which felt across almost the entire territory of the Russian Federation, including its European part.
It should be noted that the very distant correlation of geophysical fields in connection to major earthquakes have been repeatedly noted by many researchers., For instance, such effects concerning the fields of seismic noise are described in [Sobolev, 2015].
In order to investigate the coherence in connection to the Okhotsk Sea earthquake, let us consider two of the regions which are closest to the focus of the seismic event, namely, Kamchatka and the Japanese islands.
It should be noticed that the underground electrical measurements were recorded with the sampling time step of 1 hour, whereas the time series of seismic noise parameters are daily. Figure 14 shows graphs of the
10000
15000
20000
25000
15000
20000
25000
10000
15000
20000
25000
254
1.5
й It
(d)
H
0.5 I
15000
20000
25000
0.70 ■0.65 0.60 -0.55 ■0.50 '0.45 0.« -0.35 0.30 0.25 -0.20 -0.15 -0.10 -0.05
Time, hours from the beginning of 2012, right-hand end of moving time window of the length 672 hours with mutual shift 24 hours.
Fig. 11. Time-frequency diagrams for spectral measure of coherence: a - between two signals: aggregated signal of Station S5 (Altai), and the 2nd order aggregated signal of four stations in Kamchatka (Fig. 9, b); b - between three signals: aggregated signal of Station S5 (Altai), aggregated signal of Station S6 (Italy), and the 2nd order aggregated signal of four stations in Kamchatka; с and d - the same combination of signals as in (b) but the 2nd order aggregated signals from Kamchatka are built for five and six stations (Fig. 9, c, and 9, d).
Рис. 11. Частотно-временные диаграммы эволюции спектральной меры когерентности: а - между двумя сигналами: агрегированным сигналом станции S5 (Алтай) и агрегированным сигналом 2-го порядка для четырех станций на Камчатке (рис. 9, b); b - между тремя сигналами: агрегированным сигналом станции S5 (Алтай), агрегированным сигналом станции S6 (Италия) и агрегированным сигналом 2-го порядка для четырех станций на Камчатке (рис. 9, b); c и d - такая же комбинация сигналов, что и (b), но агрегированные сигналы 2-го порядка на Камчатке строятся для пяти и шести станций (рис. 9, с, и 9, d).
aggregated signal of daily time-series properties of seismic noise in Japan (see Fig. 7) and the 2nd order aggregated signal of four electrical measurement stations in Kamchatka (see Fig. 9(b)) after transition to the one-day sampling time step by calculating the average values in successive time windows of 24 hours. The aggregated signal of the seismic noise parameters dropped the initial time intervals in order both series have the same starting point, 12 December 2012. Further on, in all graphs, time marks are counted from that date.
Results of our analysis aimed at studying effects of
coherence between the aggregated signal on Kamchatka and the aggregated signal in Japan, are shown in Figure 15. The time-frequency diagram shows the evolution of the squared spectrum of the coherence between the two aggregated signals calculated in the moving time window of the length of 90 days. It shows five bursts of coherence with gradually increasing characteristic periods. At the end of the investigated time interval, we have the maximum peak of coherence corresponding to the maximum value of the period (the last position of the right-hand end of the time
Time, hours from the beginning of 2012
-1-
10000
—I— 20000
Fig. 12. Aggregated signal for increments of time series
from torsional pendulums.
Рис. 12. Агрегированный сигнал приращений временных рядов крутильных маятников.
window - 12 May 2015). It is noteworthy interesting that the first burst of coherence corresponds to the right-hand end of the time window 120-150 days from 12 December 2012 and has the characteristic period of about four days - it falls in the time interval before the mantle Okhotsk Sea earthquake of 24 May 2013 (the
163th day from 12 December 2012). Periods of five coherence maximums are approximately equal to 4, 6, 10-16 and 30-90 days, which means that the synchronization process is observed at more and more lower frequencies.
10. Conclusion
Based on the multivariate time series analysis, the following facts are established:
- Underground electrical measurements in the three different regions of the globe (Kamchatka, Altai, Italy) show a significant burst of joint coherence for oscillations with periods of about 24 hours in the time fragment 12000-16000 hours from the beginning of 2012, which corresponds to the positions of the moving time window of the length of 28 days from mid-April to the end of October 2013;
- The joint analysis of the electrical measurements and the torsional pendulum data from the four regions (Kamchatka, Altai, Italy, and Tula) shows a burst of coherence in the same time interval, from mid-April to the end of October 2013;
- The joint analysis of the electrical measurements
10000 12000 14000 16000 18000 20000 22000 24000 26000
Time, hours from the beginning of 2012, right-hand end of moving time window of the length 672 hours with mutual shift 24 hours.
Fig. 13. Estimate of the multiple coherence coefficient evolution between four signals: (1) aggregated signal of Station S5 (Altai); (2) aggregated signal of Station S6 (Italy); (3) the 2nd order aggregated signal of four stations in Kamchatka (Fig. 9(b)); (4) aggregated signal for increments of three torsional pendulums (Tula).
Рис. 13. Оценка эволюции множественной меры когерентности между четырьмя сигналами: 1) агрегированным сигналом на станции S5 (Алтай); 2) агрегированным сигналом на станции S6 (Италия); 3) агрегированным сигналом 2-го порядка для четырех станций на Камчатке (рис.9(з)); 4) агрегированным сигналом приращений трех крутильных маятников (Тула).
Time, days from 12 Dec 2012
200
400
600
800
Fig. 14. a - the 2nd order aggregated signal of four stations in Kamchatka (Fig. 9, b) after coming to the one-day sampling time step; b - aggregated signal of median values for five daily seismic noise properties according to data from F-net broadband seismic network in Japan.
Рис. 14. а - агрегированный сигнал 2-го порядка для четырех станций на Камчатке (рис. 9, Ь) после перехода к шагу по времени в одни сутки; Ь - агрегированный сигнал значений медиан пяти ежесуточных значений свойств сейсмического шума на сети широкополосных сейсмических станций F-net в Японии.
100 200 300 400 500 600 700 S00
Time, days from 12 Dec 2012, right-hand end of moving time window of the length 90 days
Fig. 15. Evolution of the squared coherence spectrum between two signals: (1) the 2nd order aggregated signal of four stations in Kamchatka (Fig. 9, b) after coming to the one-day sampling time step, and (2) aggregated signal of median values for five daily seismic noise properties according to data from F-net broadband seismic network in Japan.
Рис. 15. Эволюция квадратичного спектра когерентности между двумя сигналами: 1) агрегированным сигналом 2-го порядка для четырех станций на Камчатке (рис. 9, b) после перехода к шагу по времени в одни сутки; 2) агрегированным сигналом значений медиан пяти ежесуточных значений свойств сейсмического шума на сети широкополосных сейсмических станций F-net в Японии.
obtained in Kamchatka and the time series variation of the seismic noise parameters in the Japanese islands shows five bursts of coherence with the gradually increasing characteristic period with the maximum coherence peak at the end of the investigated time interval (May 2015).
Our hypothesis is that the burst of coherence in April - October 2013 is related to the global synchronization of the geophysical fields, which was associated with the great mantle Okhotsk Sea earthquake (M=8.3) that occurred on 24 May 2013.
In studies of such a complex multi-component system as the Earth's crust, it is highly challenging to identify a set of the main deterministic reasons, which would define all the features of the global seismic regime, particularly those which control long-term changes in the intensity of potential seismic events. Solving this problem may be facilitated by the phenomeno-logical approach based on the use of coherent noise generated by the system in the course of its evolution. For the Earth's crust, the ambient noise is a product of
its 'life'. Coherence (or synchronization) of the behavior of characteristics of a complex system, described by data of different nature and structure, is an important feature for assessments of its approach to rapid changes in the condition, which are often referred to as a 'catastrophe'. Searching for precursors of catastrophes, which may be manifested by the occurrence of synchronous components in a variety of observations, is the general idea for increasing the correlation radius of random fluctuations of parameters of a complex system as it approaches a sharp change in its properties as a result of its own dynamics [Gilmore, 1981; Nicolis, Prigogine, 1989]. This property of coherence of ambient noise of the Earth is investigated in this article.
11. Acknowledgements
The work was financially supported by the Ministry of Education and Science of the Russian Federation (Contract No. 14.577.21.0109, Project ID RFMEFI57714X0109).
12. References
Bobrovskiy V.S., 2011. The results of subterranean electric measurements on Kamchatka as global effects of proton tectogenesis: damaging earthquakes in Indonesia and China. In: P. Guarnieri (Ed.), Recent progress on earthquake geology. Nova Science Publishers, New York, p. 189-248.
Bobrovskiy V.S, Kuznetsov D.A., 2011. Cosmo-meteo-tectonics. Chapters 01-10, Cosmo-Meteo-Tectonics Distant School, Petropavlovsk-Kamchatsky, 294 p. VINITI Deponent No. 82-V2011, Feb. 2011 (in Russian) [Бобровский В.С., Кузнецов Д.А. «КосмоМетеоТектоника». Главы 01-10. Петропавловск-Камчатский: Дистантная школа «КосмоМетеоТектоника», 2011. 294 с. Деп. ВИНИТИ 24.02.2011 №82-В2011].
Brillinger D.R., 1975. Time Series. Data Analysis and Theory. Holt, Rinehart and Winston, Inc., N.Y., Chicago, San Francisco, 540 p. [Русский перевод: Бриллинджер Д. Временные ряды. Обработка данных и теория. М.: Мир, 1980. 536 с.].
Chebrov V.N., Kugaenko Yu.A., Vikulina S.A., Kravchenko N.M., Matveyenko E.A., Mityushkina S.V., Raevskaya A.A., Saltykov V.A., Chebrov D.V., Lander A.V., 2013. Deep Mw=8.3 Sea of Okhotsk Earthquake on May 24, 2013 - the largest event ever recorded near Kamchatka Peninsula for the whole period of regional seismic network operation (in Russian). Bulletin of Kamchatka Regional Association "Training and Research Center", Series: Earth Sciences (1), 17-24 (in Russian) [Чебров В.Н., Кугаенко Ю.А., Викулина С.А., Кравченко Ю.А., Матвеенко Е.А., Матюшкина С.В., Раевская А.А., Салтыков В.А., Чебров Д.В., Ландер А.В. Глубокое Охотоморское землетрясение 24.05.2013 г. с магнитудой Мw=8.3 - сильнейшее сейсмическое событие у берегов Камчатки за период детальных сейсмологических наблюдений // Вестник КРАУНЦ. Серия Науки о Земле. 2013. № 1. С. 17-24].
Cramer H., 1999. Mathematical Methods of Statistics. Princeton University Press, 575 p.
Feder J., 1988. Fractals. Plenum Press, New York, London, 284 p. [Русский перевод: Федер Е. Фракталы. М.: Мир, 1991. 254 с.].
Gilmore R., 1981. Catastrophe Theory for Scientists and Engineers. John Wiley and Sons, Inc., New York, 666 p. [Русский перевод: Гилмор Р. Прикладная теория катастроф: в 2-х книгах. М.: Мир, 1984. Кн. 1, 350 с.; Кн. 2, 285 с.].
Hannan E.J., 1970. Multiple Time Series. John Wiley and Sons, Inc., New York, London, Sydney, Toronto, 540 p. [Русский перевод: Хеннан Э. Многомерные временные ряды. М.: Мир, 1974. 575 с.].
Huber P.J., 1981. Robust Statistics. John Wiley and Sons, New York, Chichester, Brisbane, Toronto, 380 p. [Русский перевод: Хьюбер П. Робастность в статистике. М.: Мир, 1984. 303 с.].
Kalinnikov 1.1., Manukin A.B., Koneshov V.N., Matyunin V.P., Karagioz O.V., Volfson G.B., 2011. Investigation of variable gravitational gradients and specific features of the microseismic background with a torsion balance. lzvestiya, Physics of the Solid Earth 47 (5), 456-463. http://dx.doi.org/10.1134/S1069351311040033.
Kantelhardt J.W., Zschiegner S.A., Konscienly-Bunde E., Havlin S., Bunde A., Stanley H.E., 2002. Multifractal detrended fluctuation analysis of nonstationary time series. Physica A: Statistical Mechanics and its Applications 316 (1-4), 87-114. http://dx.doi.org/10.1016/S0378-4371(02)01383-3.
Kuznetsov D.A. 1991. Practice of short-term prediction of earthquakes: astro-, cosmo-geophysical impulses of Vernad-sky-Vlasov-Vorobjev-Prigozhin on the vertical sequence of underground electrodes at Pedlnstitute fault at the magnetic meridian of Petropavlovsk-Kamchatsky. KGPI, Petropavlovsk-Kamchtasky, 9 p., VINITI Deponent No. 3256-V91, July 1991 (in Russian) [Кузнецов Д.А. Практика краткосрочного прогноза землетрясений: аст-ро-космо-геофизические импульсы Вернадского-Власова-Воробьева-Пригожина на вертикальной последовательности подземных электродов в разломе «Пединститутский» на магнитном меридиане Петропав-ловска-Камчатского. Петропавловск-Камчатский: КГПИ, 1991. 9 с. Деп. в ВИНИТИ 30.07.91. №3256-В91].
Lyubushin A.A., 1998a. An aggregated signal of low-frequency geophysical monitoring systems. lzvestiya, Physics of the Solid Earth 34 (3), 238-243.
Lyubushin A.A., 1998b. Analysis of low-frequency multidimensional time series for geophysical monitoring and earthquake prediction. Journal of Earthquake Prediction Research 7 (4), 496-509.
Lyubushin A.A., 1999. Analysis of multidimensional geophysical monitoring time series for earthquake prediction. Annali di Geofisica 42 (5), 927-937. http://dx.doi.org/10.4401/ag-3757.
Lyubushin A.A, 2007. Analysis of Geophysical and Ecological Monitoring Systems Data. Nauka, Moscow, 228 p. (in Russian) [Любушин А.А. Анализ данных систем геофизического и экологического мониторинга. М.: Наука, 2007. 228 с.].
Lyubushin A.A., 2009. Synchronization trends and rhythms of multifractal parameters of the field of low-frequency microseisms. lzvestiya, Physics of the Solid Earth 45 (5), 381-394. http://dx.doi.org/10.1134/S10693513090 50024.
Lyubushin A.A., 2010a. The statistics of the time segments of low-frequency microseisms: trends and synchronization. lzvestiya, Physics of the Solid Earth 46 (6), 544-554. http://dx.doi.org/10.1134/S1069351310060091.
Lyubushin A.A., 2010b. Multifractal parameters of low-frequency microseisms. Chapter 15. In: V. de Rubeis, Z. Czechowski, R. Teisseyre (Eds.), Synchronization and Triggering: from Fracture to Earthquake Processes. Springer-Verlag, Berlin, Heidelberg, p. 253-272.http://dx.doi.org/10.1007/978-3-642-12300-9_15.
Lyubushin A.A., 2010c. Synchronization of multifractal parameters of regional and global low-frequency microseisms. Geophysical Research Abstracts 12, EGU2010-696 (EGU General Assembly). Available from: http:// meetingorganizer.copernicus.org/EGU2010/EGU2010-696.pdf.
Lyubushin A.A., 2011a. Cluster analysis of low-frequency microseismic noise. lzvestiya, Physics of the Solid Earth 47 (6), 488-495. http://dx.doi.org/10.1134/S1069351311040057.
Lyubushin A.A., 2011b. Seismic catastrophe in Japan on March 11, 2011: Long-term prediction on the basis of low-frequency microseisms. lzvestiya, Atmospheric and Oceanic Physics 46 (8), 904-921 [Любушин А.А. Сейсмическая катастрофа в Японии 11 марта 2011 года. Долгосрочный прогноз по низкочастотным микросейсмам // Геофизические процессы и биосфера. 2011. Т. 10. № 1. С. 9-35]. http://dx.doi.org/10.1134/S00014 33811080056.
Lyubushin A., 2012. Prognostic properties of low-frequency seismic noise. Natural Science 4 (8A), 659-666. http://dx.doi.org/10.4236/ns.2012.428087.
Lyubushin A., 2013. How soon would the next mega-earthquake occur in Japan? // Natural Science 5 (8A1), 1-7. http://dx.doi.org/10.4236/ns.2013.58A1001.
Lyubushin A.A., 2014a. Dynamic estimate of seismic danger based on multifractal properties of low-frequency seismic noise. Natural Hazards 70 (1), 471-483. http://dx.doi.org/10.1007/s11069-013-0823-7.
Lyubushin A.A., 2014b. Analysis of coherence in global seismic noise for 1997-2012. lzvestiya, Physics of the Solid Earth 50 (3), 325-333. http://dx.doi.org/10.1134/S1069351314030069.
Lyubushin A.A., 2015. Wavelet-based coherence measures of global seismic noise properties. Journal of Seismology 19 (2), 329-340. http://dx.doi.org/10.1007/s10950-014-9468-6.
Lyubushin A.A., Kopylova G.N., 2004. Multidimensional wavelet analysis of time series of electrotelluric observations in Kamchatka. lzvestiya, Physics of the Solid Earth 40 (2), 163-175.
Mallat S., 1998. A wavelet tour of signal processing. Academic Press, San Diego, London, Boston, New York, Sydney, Tokyo, Toronto, 577 p. [Русский перевод: Малла С. Вэйвлеты в обработке сигналов. М.: Мир, 2005. 671 c.].
Marple S.L. (Jr.), 1987. Digital Spectral Analysis with Applications. Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 492 p. [Русский перевод: Марпл С.Л. Цифровой спектральный анализ и его приложения. М.: Мир, 1990. 584 с.].
Martynov O.V., 2008. The natural catastrophes forecasting system concept and practical results obtained from nonlinear physics, mathematics and system data. Nelineinyi mir (Nonlinear World) 6 (10), 579-615 (in Russian) [Мартынов О.В. Концепция системы прогноза природных катастроф и практические результаты, полученные на основе аппарата нелинейной физики, математики и данных системы // Нелинейный мир. 2008. Т. 6. № 10. С. 579-615].
Martynov O. V., Semenov L.L., Kurotchenko S.P., Legkov A. V., Parshutin R. V., 2006b. System of reading, storage and information processing of four channel wideband gradiometer WBGM-4. Proceedings of Tula State University. Geodynamics, Physics, Mathematics, Thermodynamics, Geoecology Series (3), 68-77 (in Russian) [Мартынов О.В.,
СеменовЛ.Л., Куротченко С.П., Легков А.В., Паршутин Р.В., 2006. Система считывания, хранения и обработки информации четырехканального широкополосного градиентометра ШГММ-4 // Известия Тульского государственного университета. Серия геодинамика, физика, математика, термодинамика, геоэкология. 2006. Вып. 3. С. 68-77].
Martynov O.V., Semenov L.L., Surkov A.V., Shopin S.A., 2006a. System of reading, storage and information processing of three channel wideband gradiometer. Proceedings of Tula State University. Geodynamics, Physics, Mathematics, Thermodynamics, Geoecology Series (3), 59-68 (in Russian) [Мартынов О.В., Семенов Л.Л., Сурков А.В., Шопин С.А. Система считывания, хранения и обработки информации трехканального широкополосного градиентометра // Известия Тульского государственного университета. Серия геодинамика, физика, математика, термодинамика, геоэкология. 2006. Вып. 3. С. 59-68].
Martynova M., Martynov O, 2006. Physico-mathematical base and monitoring system for the earthquake forecast (M>6), its place and time. In: Proceedings of First European Conference on Earthquake Engineering and Seismology (3-6 September 2006). Geneva, Switzerland, p. 3294-3300.
Nicolis G., Prigogine 1., 1989. Exploring Complexity, an Introduction. W.H. Freedman and Co., New York, 328 p. [Русский перевод: Николис Г., Пригожин И. Познание сложного. М.: Мир, 1990. 344 с.].
Shopin S.A., 2009. Instrumentation system for registration of ultra low frequency gravitational field disturbances. In: Proceedings of International Conference on Ecology, Energy, Economy Security in a Non Linear World (3E-Security) (Belgrade, Serbia, 6-7 November 2009). Swiss Association "NON-LINEARITE, Geneva, p. 86-100. http://dx.doi.org/10.13140/RG.2.1.4757.6161.
Shopin S.A. 2014. Influence of microseisms and variations of atmospheric pressure on the instrumentation systems based on horizontal torsion balance. Proceedings of Tula State University. Natural Sciences 1 (1), 249-263 (in Russian) [Шопин С.А. Влияние микросейсм и вариаций атмосферного давления на измерительные системы на основе горизонтальных крутильных весов // Известия Тульского государственного университета. Естественные науки. 2014. Вып. 1. Ч. 1. С. 249-263].
Shopin S.A., 2015. On the works of O.V. Martynov on earthquake prediction. In: A.E. Fedorov (Ed.), Planet Earth System: 200 years of Holy Alliance. LENAND, Moscow, p. 102-120 (in Russian) [Шопин С.А. О работах О.В. Мартынова по прогнозу землетрясений // Система "Планета Земля": 200 лет Священному союзу / Ред. А.Е. Федоров. М.: ЛЕНАНД, 2015. С. 102-120].
Sobolev G.A., 2015. Seismic Noise. Nauka i Obrazovanie, Moscow, 272 p. (in Russian) [Соболев Г.А. Сейсмический шум. М.: ООО «Наука и образование», 2014. 272 с.].
Sobolev G.A, Ponomarev V.A, 2003. Physics of Earthquakes and Precursors. Nauka, Moscow, 270. (in Russian) [Соболев Г.А., Пономарев А.В. Физика землетрясений и предвестники. М.: Наука, 2003. 270 с.].
Taqqu M.S., 1988. Self-similar processes. In: Encyclopedia of statistical sciences, vol. 8. Wiley, New York, p. 352-357.
Uyeda S., Nagao T., Orihara Y., Yamaguchi T., Takahashi 1., 2000. Geoelectric potential changes: Possible precursors to earthquakes in Japan. Proceedings of the National Academy of Sciences of the United States of America 97 (9), 4561-4566. http://dx.doi.org/10.1073/pnas.97.9.4561.
Varotsos P., Alexopoulos K., Lazaridou M., 1993. Latest aspects of earthquake prediction in Greece based on seismic electric signals, II. Tectonophysics 224 (1-3), 1-37. http://dx.doi.org/10.1016/0040-1951(93)90055-0.
Lyubushin, Alexey A., Doctor of Physics and Mathematics, Head of laboratory O.Yu. Schmidt Institute of Physics of the Earth of RAS 10 Bol'shaya Gruzinskaya street, Moscow D-242 123242, GSP-5, Russia И e-mail: [email protected]
Любушин Алексей Александрович, докт. физ.-мат. наук, заведующий лабораторией
Институт физики Земли им. О.Ю. Шмидта РАН
123242, ГСП-5, Москва Д-242, ул. Большая Грузинская, 10, Россия
И e-mail: [email protected]
Bobrovskiy, Vadim S., Director Distant School 'Cosmic-Meteo-Tectonics'
15/2 Prospect 50 Let Oktyabrya, office 16, Petropavlovsk-Kamchatsky 683024, Russia e-mail: [email protected]
Бобровский Вадим Сергеевич, директор Дистантная школа 'Космо-Метео-Тектоника'
683024, Петропавловск-Камчатский, проспект 50 Лет Октября, 15/2, офис 16, Россия e-mail: [email protected]
Shopin, Sergey А., chief specialist Tula Geological Exploration Crew LLC 13 Smidovich street, Tula 300028, Russia e-mail: [email protected]
Шопин Сергей Александрович, главный специалист ООО «Тульская геологоразведочная партия» 300028, Тула, ул. Смидович, 13, Россия e-mail: [email protected]
