Analysis of geoacoustic emission and electromagnetic radiation signals accompanying earthquake with magnitude Mw = 7.5

The paper is devoted to the analysis of frequency spectra and pulse waveform variety of the geoacoustic and electromagnetic signals recorded on Kamchatka Peninsula at “Karymshina” site during seismically calm and active periods. Signal pre-processing includes pulse detection and their waveforms reconstruction. A frequency spectrum is analyzed using the Adaptive Matching Pursuit algorithm. To study a variety of waveforms, each pulse is encoded by a special descriptive matrix. Then pulse classification based on similarity of the descriptive matrices is performed. Thus, a signal alphabet is formed. The authors analyzed the geophysical signals recorded before, during and after the earthquake with the magnitude Mw = 7.5 dated March 25, 2020. The obtained estimates of frequency spectra and signal alphabets are compared with the analysis results of signal recoded during the seismically calm period of March 22, 2020.


Introduction
The results of earlier studies [1][2][3] show presence of a pre-seismic response in geoacoustic emission and VLF electromagnetic radiation signals. In some cases, this response is a characteristic change in a signal frequency [4]. Study and analysis of dynamics of the geoacoustic and electromagnetic pulse streams are promising tasks in view of identifying complex earthquake precursors.
Registration of geophysical signals is carried out on Kamchatka Peninsula at "Karymshina" site (52.83º N, 158.13º E). Monitoring of electromagnetic signals is carried out using a multichannel VLF recorder. This recorder allows us to receive weak noise signals and consists of three antennas, pre-amplifiers, a filter device, and an analog-todigital converter. Geoacoustic signals are observed using piezoceramic receivers located near the bottom of artificial reservoirs [5]. Electromagnetic and geoacoustic signals are sampled at the rate of 48 kHz. Both types of the signals are pulse streams. A single electromagnetic signal of natural origin propagating in the Earth-ionosphere waveguide is called "sferic". A single pulse of geoacoustic signal is called "geoacoustic pulse". This article is devoted to the analysis of frequency spectrum and waveforms of the geoacoustic pulses and electrical component of sferics recorded under good weather conditions (no strong wind and precipitation) during the seismically calm period on March 22, 2020 (3665 geoacoustic pulses and 27436 sferics) and during the earthquake dated March 25, 2020 (3050 geoacoustic pulses and 12074 sferics). The earthquake with the magnitude M w = 7.5 was registered on March 25, 2020 at 02:49:21 UTC at the depth of 57.8 km, at the epicentral distance of 430 km. The epicenter coordinates are 48.964° N and 157.696° E (NEIC, https://earthquake.usgs.gov/earthquake). According to the Earthquakes Catalogue for Kamchatka and the Commander Islands (1962 -present) provided by KB GS RAS the earthquake local magnitude M L is 7.7 (http://sdis.emsd.ru/info/earthquakes/catalogue.php). This is one of the strongest regional earthquakes in the last 20 years.

Experiment
In general, the process of geophysical signal analysis can be represented as the diagram shown in Fig. 1. Let's consider the analysis stages in more detail. Firstly, signal pre-processing is performed. It includes removing a trendline by applying the high-pass Butterworth filter with the cutoff frequency of 1 kHz, centering and normalizing.
Then search and detection of pulses are performed. The adaptive threshold scheme is used. It adjusts to the noise level of a processed signal. A signal fragment detected as a pulse is checked for integrity of the sequence of local maxima and minima recorded along its length. If the length of the continuous sequence of local extrema belongs to the given interval from N min to N max (in the experiment, N min = 3, N max = 200) then this signal fragment is a pulse. And vice versa, if the number of local extrema continuously following one after another lies outside the specified interval, the detected signal fragment is not a pulse. Detailed description of the pulse detection procedure is presented in [6,7].
Distortions of a geophysical pulse waveform negatively affect quality of the signal processing and analysis; therefore, the next stage is waveform reconstruction of each detected pulse using wavelet denoising. The family of fourth-order Symlet wavelets were used for forward and inverse wavelet transforms. Detail coefficient modification obtained by forward wavelet transform was carried out according to the posterior median rule. Threshold values were calculated using the Empirical Bayes method [8]. A waveform reconstruction example of some geoacoustic pulses and sferics is shown in Fig. 2. Next, frequency spectrum and waveform of the reconstructed pulses are analyzed. The features of time-frequency analysis of sferics and geoacoustic pulses using sparse approximation are described in detail in [9,10]. To determine a time-frequency pulse structure, the Adaptive Matching Pursuit algorithm (AMP) is proposed. After applying the algorithm, the analyzed pulse is represented as a linear combination of basic functions.
The results presented in [10] show that geoacoustic pulses and sferics are described accurately by the modulated Berlage functions: where A is the pulse amplitude; t is the time, 0 ≤ t ≤ t end ; p max is the pulse envelope maximum position defined relating to t end ; n(p max ) is the minimum value of the parameter affecting pulse envelope slope chosen so that g(p max · t end ) ≤ 0.05 g(t end ); Δ is the coefficient of pulse envelope slope variation, the more Δ value is, the steeper the envelope is, Δ > 1; f is the function frequency; φ 0 is the initial phase. During the computational experiment, 21 fifteen-minute records of geoacoustic signals were processed and analyzed (16 records made under good weather conditions during the seismically calm period on March 22, 2020; 5 records made in the time interval ±30 minutes near the earthquake registered on March 25, 2020 at 02:49:21 UTC). 22 fifteen-minute records of electromagnetic signals were processed and analyzed (17 records made under good weather conditions during the seismically calm period; 5 records made in the time interval ±30 minutes near the earthquake). A total of 6715 geoacoustic pulses and 39510 sferics were detected.
For each detected pulse, the set consisting of 216 Berlage functions was generated. The functions have following parameters:  8 values of f from 1 kHz to 20 kHz;  3 values of t end (at a percentage of a maximum possible function duration, the maximum function duration coincides with the analyzed pulse duration) from 50% to 100%;  3 values of p max (at a percentage of the actual function duration) from 5% to 25%;  3 values of  from 1.1 to 5;  initial phase φ 0 is 0. The pulses were analyzed using the AMP algorithm with the following parameters:  maximum error level is 5%,  maximum number of refinement iterations is 5,  refinement accuracy is 5 Hz, 2%, 2%, and 0.1 for the parameters f, t end , p max , and  respectively. To visualize signal sparse approximation on time-frequency plane, the Wigner-Ville transform is used. Fig. 3 shows the results of applying the AMP algorithm to a geoacoustic pulse and a sferic.  The detected geoacoustic pulses and sferics have a wide variety of waveforms. To identify and classify pulse waveforms, the structural description method is proposed. It is to encode a pulse with special binary matrix called "descriptive matrix" [7]. To construct the descriptive matrix, amplitudes of pulse local extrema and time intervals between them are compared:   Fig. 4 shows a pulse and a corresponding descriptive matrix. Next, the described pulses are classified according to similarity of their descriptive matrices. The pulses with descriptive matrices that coincide by more than 80% belong to the same class. Each class can be considered as a symbol, and a set of symbols forms a signal alphabet. It is proposed to visualize a signal alphabet as a histogram. Orders of descriptive matrices are plotted along the horizontal axis; and probabilities of pulses described by corresponding symbol are plotted along the vertical axis. A visualization example of fifteen-minute geoacoustic signal alphabet is shown in Fig. 5.  Table 1 lists main characteristics of the processed geoacoustic signals. The signal records were taken from the database of Acoustic research laboratory of IKIR FEB RAS (September 2016 -present).  Table 1 shows that pulse repetition rate is higher in the time interval ±30 minutes near the earthquake than in the seismically calm period, but pulse duration and amplitude differ slightly.

Analysis results of geoacoustic emission signals
For each fifteen-minute signal record, a frequency histogram of the function with maximum absolute value of the approximation coefficients (hereinafter "first decomposition function") was constructed. The analysis results of the histograms show that pulses with the frequencies below 3 kHz prevail in the time interval near the earthquake. So, the histograms have a pronounced peak at the corresponding frequencies. Fig. 6 shows histograms plotted for the geoacoustic signals recorded on March 22, 2020 at 06:45, 17:45 and 18:30 UTC. There is no peak at the frequencies below 3 kHz on the histograms.  The analysis results obtained by the structural description method show that the alphabets of background signals and signals accompanying the earthquake also differ. The alphabets of signals registered on March 25, 2020 contains a symbol encoding a prevalent waveform (more than 20% of the detected pulses). Therefore, a single peak corresponding to the most common symbol appears on alphabet histograms. The alphabet histograms of signals registered on March 22, 2020 at 06:45, 17:45, and 18:30 UTC are shown in Fig. 8. The background pulses have no prevalent waveform. Fig. 9 illustrates alphabet histograms of signals registered on March 25, 2020 at 02:15, 02:45 and 03:00 UTC. There is a peak corresponding to the most common symbol on these histograms.
The authors assume that the listed differences are associated with the pulse generation peculiarities during the earthquake development. These peculiarities consist in increasing the number of pulses with a certain frequency and waveform.   Table 2 lists main characteristics of the processed electromagnetic signals. The signal records were provided by Electromagnetic radiation laboratory of IKIR FEB RAS.  Table 2 shows that the sferics recorded before and after the earthquake are have more duration, amplitude, and repetition rate than the sferics recorded during the seismically calm period. It should be noted that the fifteen-minute electromagnetic signals recorded on March 25, 2020 are very different from the background signals recorded on March 22, 2020. Fig. 10 shows graphs of two electromagnetic signals recorded on March 22, 2020 and March 25, 2020. The pre-seismic signal contains more pulses and has higher average amplitude.

Analysis results of VLF electromagnetic radiation signals
Analysis by sparse approximation method shows that frequency spectra of the background signals and signals accompanying the earthquake differ. The first decomposition function frequency is distributed from 4 to 18 kHz for the background sferics. But it belongs to the range from 0 to 2 kHz for the pre-and post-seismic sferics. Fig. 11 shows frequency histograms plotted for the electromagnetic signals recorded on March 22, 2020 at 06:39, 17:39 and 18:24 UTC.  Due to peculiarities of signal propagation in the Earth-ionosphere waveguide, there should be no frequencies below the waveguide critical frequency in sferics (according to calculations the waveguide critical frequency is 1709 Hz for the studied signals). Thus, a classical sferic does not contain frequencies bellow 1709 Hz. This is confirmed by the results of analysis of the signals recorded on March 22, 2020. All sferics having the timefrequency structure different from the declared one can be considered as anomalies. Such differences are caused by the presence of inhomogeneities (areas with an increased electron concentration) in the Earth-ionosphere waveguide. Appearance of these inhomogeneities is associated with the earthquake preparation [11].
Alphabets of the background signals and signals accompanying the earthquake also have differences, but less obvious. The alphabets of pre-and post-seismic signals contain the symbols described by matrices of small orders from 3 to 5. Fig. 13-14 show the alphabet histograms of the electromagnetic signals recorded in the seismically calm and active periods.

Conclusions
The presented work is devoted to complex analysis of the background, pre-and postseismic pulse signals of geoacoustic emission and VLF electromagnetic radiation. The signals recorded under good weather conditions during the seismically calm period of March 22, 2020 and during the time interval ±30 minutes near the earthquake registered on March 25, 2020 at 02:49:21 UTC (the depth is 57.8 km; the epicentral distance is 430 km; the epicenter coordinates are 48.964° N and 157.696° E; the magnitude M w is 7.5) were processed and analyzed.
It was found that frequency spectra and waveforms of the pulses accompanying the earthquake and pulses registered in the background period differ. The analysis was carried out by sparse approximation and the structural description method. The first method is used to analyze a signal frequency spectrum. The second one is used to analyze a signal waveform.
The frequencies below 3 kHz are present in spectra of many geoacoustic pulses recorded during the time interval ±30 minutes near the earthquake. More than 20% of the pulses have a similar waveform. The authors suggest that this is associated with increasing the number of pulses with a certain frequency and waveform during the earthquake development.
Sferics with frequencies below the critical frequency of the Earth-ionosphere waveguide appear in the pre-and post-seismic electromagnetic signals. These sferics are encoded by symbols with descriptive matrices of small orders. Such symbols correspond to the pulse waveforms containing 3-5 local extrema. These peculiarities are caused by the presence of inhomogeneities at the waveguide boundaries. The inhomogeneity appearance is associated with the earthquake preparation.
The obtained results are consistent with the results of the previous study devoted to analysis of the signals recorded before the Zhupanov earthquake registered on January 30, 2016 at 03:25 UTC [4]. In the future, the presented work will be useful for identifying geophysical short-term precursors of strong Kamchatka earthquakes. A study of sferics and geoacoustic pulses preceding other Kamchatka seismic events is planned.