EXPERIENCE OF COMPLEXATION OF GLOBAL GEOPHYSICAL OBSERVATIONS

The technique of joint analysis of heterogeneous time series of geophysical monitoring systems for the pur‐ pose 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 electri‐ cal 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 pro‐ perties.


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 multivariate 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 dimensionless 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 nonstationary, 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.

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 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).Multielectrode 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, а).Iron electrodes (5005003 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 fourelectrode 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 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-Kamchatsky). 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.The initial records present signals sampled at onesecond 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 , subtracting mean values and dividing by , and clipping values which exceed thresholds 4 .These iterations are repeated for each channel until the values of 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.

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 [Martynov, 2008;Martynova, Martynov, 2006].These measurements are performed in Tula State University (Tula, Russia).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 µm; 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 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 photocurrent 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 roundthe-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[Martynov et al., , 2006b;;Shopin, 2009Shopin, , 2014Shopin, , 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.

SEISMIC NOISE PARAMETERS
In ours study, continuous seismic records from Fnet 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 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 [Lyubushin, 2009[Lyubushin, , 2010a[Lyubushin, , 2010b[Lyubushin, , 2010c[Lyubushin, , 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[Lyubushin, , 2012[Lyubushin, , 2013[Lyubushin, , 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]: where ∆ 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 ∆ 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[Lyubushin, , 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: where c k ,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 symlets 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[Lyubushin, , 2014b]].Multifractal parameters Δ and * .Multifractal singularity spectrum F(a) [Feder, 1988] of signal x(t) is defined as a fractal dimensionality of time moments   [Feder, 1988].If X(t) is a usual self-similar mono-fractal signal [Taqqu, 1988] with Hurst exponent value 0 1, then 1, 0∀ , 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: , , Δ , and * -an argument providing maximum to singularity spectra: * max .Parameter * is called a generalized Hurst exponent, and it gives the most typical value of Lipschitz-Holder exponent.Parameter Δ , singularity spectrum support width, can be regarded as a measure of variety of stochastic behavior.For removing scaledependent 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.
Another parameter of seismic noise is logarithm of variance, lg(Var).
Values lg(k), En, Δ , * 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, Δ , * 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.

METHODS OF ANALYSIS: MULTIPLE COHERENCE COEFFICIENT
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 .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 .Quantity 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: where is frequency, -spectral matrix of time series X(t) of size m×m, -cross-spectral matrix of e size m×n, , H is the sign of Hermitian conjugation.Quantity replaces the square of the modulus of the coherence spectrum in case of two multidimensional signals.
We now introduce the concept of component-bycomponent canonical coherences of q-dimensional time series Z(t) as squares of the moduli of the Рис. 6. Положение 78 широкополосных сейсмических станций на Японских островах.maximum canonical coherence in the situation when, in formula (1), the i-th scalar component of q-dimensional 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 .
By introducing the component-by-component canonical coherence, it becomes possible to determine one more frequency-dependent statistics [Lyubushin, 1999[Lyubushin, , 2007]], which characterizes the coherence of variations in all the components of vector series Z(t) at frequency : , 4 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 .Frequencydependent 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 be the time coordinate of the window having a length of L counts.Calculating the spectral matrices for the samples falling in time window , we obtain the two-parameter function, , .The bursts of the , 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 , with size q×q 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: , 5 where A k are matrices of autoregression parameters of size qq, p is autoregression order, e(t) is a q-dimensional 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 Φ 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, winsorization 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 A k and Φ, 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: where is a complex matrix: exp , 7 E is a unit matrix of size qq.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 , dependence.

METHODS OF ANALYSIS: FOURIER-AGGREGATED SIGNAL
An aggregated signal is constructed in two stages.At the 1 st 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 2 nd 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[Lyubushin, , 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 i (t) and filter (q-1) -dimensional series Z (i) (t) consisting of the all other components in such a way that scalar signal obtained at the filter output has the maximum coherence with extracted series Z i (t) at each frequency.In order to do this, components of the eigenvector of matrix (3), where Z i (t) appears as ( ) Y t and Z (i) (t) ap- pears as X(t), corresponding to its maximum eigenvalue (obviously equal to ), should be used as a multidimensional frequency filter [Lyubushin, 1998a[Lyubushin, , 1998b[Lyubushin, , 2007]].If component Z i (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 simply by its construction, and this is the meaning of such an operation.At the same time, series retains all the components of Z i (t) which are common for the other components of series Z(t), i.e., for signal Z (i) (t).
We now determine aggregated signal A Z (t) of multidimensional time series Z(t) as the first principal spectral component of multidimensional series C (Z) (t) composed of canonical components 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 A Z (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 A Z (t).Although the common components are extracted in the course of either filtering, aggregated signal A Z (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 periodogram 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 periodograms 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 socalled 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 , 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[Lyubushin, , 2007]].

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 2 nd 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 2 nd order aggregated signals for various numbers of stations in Kamchatka which have simultaneous measurements.Figure 9 shows results of the 2 nd order aggregation for combinations of different stations in Kamchatka.
Figure 10 illustrates the time-frequency structure of the 2 nd order aggregated signals.It is evident that the increase in the number of stations has an insignificant impact on the 2 nd order aggregated signal spectral structure.When choosing the number of stations included in the 2 nd 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 Рис. 8. Графики агрегированных сигналов от каждой станции электрических измерений, построенные после селекции каналов.Начальные моменты времени у сигналов различаются, но все они заканчиваются 10.05.2015 г. 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 2 nd 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 2 nd 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 7 th 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.

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 Рис. 9. Графики агрегированных сигналов 2-го порядка от различного числа станций электрических измерений (a-e) на Камчатке.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 7 th 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.

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 aggregated signal of daily time-series properties of seismic noise in Japan (see Fig. 7) and the 2 nd order aggregated signal of four electrical measurement stations in Kamchatka (see Fig. 9(b)) after transition to the oneday 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

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 Рис.12. Агрегированный сигнал приращений временных рядов крутильных маятников.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 phenomenological 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.

Fig. 8 .
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).

Fig. 9 .
Fig. 9. Graphs of the 2 nd order aggregated signals from different number of electrical measurements stations (a-e) in Kamchatka.

Fig. 12 .
Fig. 12. Aggregated signal for increments of time series from torsional pendulums.