Principles of Organizing Earthquake Forecasting Based on Multiparameter Sensor-WEB Monitoring Data

The paper describes an approach that allows, basing on the data of multiparameter monitoring of atmospheric and ionospheric parameters and using ground-based and satellite measurements, to select from the data stream a time interval indicating the beginning of the final stage of earthquake preparation, and finally using intelligent data processing to carry out a short-term forecast for a time interval of 2 weeks to 1 day before the main shock. Based on the physical model of the lithosphere-atmosphericionospheric coupling, the precursors are selected, the ensemble of which is observed only during the precursory periods, and their identification is based on morphological features determined by the physical mechanism of their generation, and not on amplitude selection based on statistical data processing. Basing on the developed maquette of the automatic processing service, the possibility of real-time monitoring of the situation in a seismically active region will be demonstrated using the territory of the Kamchatka region and the Kuril Islands. 1 Precursory period identification In [1], short-term precursors of earthquakes are listed, the generation of which is due to the development of a cascade process of the Earth's outer geo-shells coupling and which we use to create a short-term forecast (if data on these parameters are available). As shown in [2]. the processes that cause the generation of these precursors are genetically related to processes in the earth's crust and used in seismology as precursors, in particular, a decrease in the value of the parameter b in the Gutenberg-Richter relation [3], foreshocks [4], etc. At the same time, continuously there are discussions about the problem of correlation between the probabilistic approach to forecasting and determinism conditioned by the presence of precursors. In comparison with the uncompromising position of the seismologists of the 90s of the last century [5], the positions of the convergence of seismologists and scientists engaged in the study of physical precursors of earthquakes have emerged [6]. It seems to us * Corresponding author: pulse@rssi.ru © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). E3S Web of Conferences 196, 03004 (2020) https://doi.org/10.1051/e3sconf/202019603004


Precursory period identification
In [1], short-term precursors of earthquakes are listed, the generation of which is due to the development of a cascade process of the Earth's outer geo-shells coupling and which we use to create a short-term forecast (if data on these parameters are available). As shown in [2]. the processes that cause the generation of these precursors are genetically related to processes in the earth's crust and used in seismology as precursors, in particular, a decrease in the value of the parameter b in the Gutenberg-Richter relation [3], foreshocks [4], etc. At the same time, continuously there are discussions about the problem of correlation between the probabilistic approach to forecasting and determinism conditioned by the presence of precursors. In comparison with the uncompromising position of the seismologists of the 90s of the last century [5], the positions of the convergence of seismologists and scientists engaged in the study of physical precursors of earthquakes have emerged [6]. It seems to us that the combination of both approaches is quite acceptable, moreover, it is natural. Here's a simple analogy. An egg lying on the edge of a table is always in danger of falling to the floor. We do not know when and for what reason this will happen, whether from a sharp stream of air when opening the door of the apartment, or whether a person passing by will touch him with the hollow clothes, or a small child will be interested in a white ball lying on the table. But when this happens, the end is inevitable, the egg will fall to the floor and break. So, the period of time until the egg flies from the edge of the table to the floor, we can call the precursor period, meaning the presence of short-term precursors of earthquakes. In this case, the main task is the timely detection and identification of those precursors that are characteristic to a period of time when a return to a state of equilibrium is impossible, a state of irreversibility. In fact, all the studies of physical precursors in recent years, conducted by our group, were aimed at identifying the most reliable precursors (and their combinations), which would unequivocally indicate the approach of an earthquake within the next few days from the moment the precursor was discovered. At the same time, efforts were made not only to develop technology for the unambiguous identification of precursors, but also to automate this process. That would allow real-time precursor monitoring [7]. It should be noted that the procedure for identifying precursors is based not on the assessment of the amplitude of the deviation of the ionospheric parameters from the unperturbed value, but on the recognition of their morphological characteristics due to the physical mechanism of their generation. The process actually boils down to the recognition of the precursor image, where, for example, in the case of ionospheric precursors, an image was developed  the "mask" of the precursor, which manifests itself in the form of a strong positive variation of the electron concentration over the earthquake preparation zone [8]. This is the only parameter for the interpretation of which human intervention is required, since it is the image that is recognized in the form of a "stalactite-stalagmite" structure -vertical structures on the "mask" image corresponding to the evening and morning positive variations (Fig. 1). In fact, this is an ideal picture, and characteristic, only for this region, because nighttime positive variations can appear not only a day before the earthquake, but also several days (up to 10), and not only during one night, but several nights in a row. In Fig. 2 shows the mask of ionospheric variations at Tsukuba station during the preparation of the M9.1 megaearthquake in Japan on March 11, 2011.

Fig. 2.
Ionospheric precursor mask for the Tohoku earthquake, Japan, March 11, 2011 according to the Tsukuba GPS receiver As one can see, positive nighttime variations begin to appear a week before the earthquake and continue until the event itself and the next night after it. The question arises as to what caused the positive variations during the day (we see them between the morning and evening variations). The thing is that the whole week before the Tohoku earthquake and on the day of the earthquake itself was characterized by very high solar and geomagnetic activity. From February 27 (50 DOY) to March 8 (67 DOY), the solar radio flux index F10.7 increased from 90 to 155, and, as is known, there is a very high correlation between the total electron content and F10.7, with the electron content lagging by about 2 days in relation to F10.7 [9]. In addition, two moderate magnetic storms took place on March 3 and 11. These two factors contributed to the daytime TEC disturbances. Variations of F10.7 and the global equatorial geomagnetic index Dst are shown in Fig. 3. It was mentioned above that the precursor period is more reliably distinguished using a combination of precursors. As such a combination, we use variations in the total electron content (or the critical frequency foF2 in the presence of an ionosonde inside the earthquake preparation zone), anomalous flux of outgoing infrared radiation in the long-wavelength frequency range OLR [10] and corrections for the chemical potential of water vapor in the surface layer of the atmosphere [11]. Both of these parameters are associated with the process of ionization of the surface layer of the atmosphere by alpha particles emitted by radon during the intensification of its release from the earth's crust before an earthquake. Condensation of water vapor occurs on the new formed ions, and latent heat is released, which is the source of OLR. Hydration of ions leads to a decrease in relative humidity and an increase in air temperature, causing a failure of the equilibrium state of the functions of temperature and pressure in the atmosphere [12]. Fig. 4 shows the dynamics of the spatial distribution of anomalies of the outgoing flow of long-wave infrared radiation OLR during the preparation of the Tohoku earthquake [13]. Now let's turn to the variations of the chemical potential correction. Unfortunately, satellite data with a sufficiently high time resolution (3 hours) appeared only in September 2011, so we will use the data of the Ishinomaki ground based meteorological station (Fig. 5). Quite a distinct maximum is observed on March 8 (67 DOY).
One of the well-proven precursors revealed by multi-year analysis of precursors is the cross-correlation coefficient between ionospheric monitoring data (ground-based ionosondes or GPS receivers) for stations located at different distances from the epicenter [14]. For the Tohoku case, we used data from two Japanese ionosondes: Kokubunji near the epicenter and Yamagawa in the very south of Japan, Kyushu island. As an earthquake approaches, the cross-correlation coefficient according to [14] should decrease, which we observe in Fig. 6. We see that the minimum is also reached on March 8 (67 DOY).  Let us summarize. To determine the precursor period, we have data from four sources: variations in the GPS TEC (63-71 DOY), variations in the OLR E_index (66-71 DOY), variations in the chemical potential (67 DOY), and variations in the cross-correlation coefficient (67 DOY). We can say with confidence that the main day giving grounds for a short-term forecast is March 8 (3 days before the earthquake), and we can determine the precursor period (for reliability, we take the coincidence of OLR and GPS TEC) days from March 7 to 11, 2011. The precursor period is a reference point for determining the date of an earthquake. To improve the accuracy, a multi-year analysis of short-term precursors in the selected region is required, from which we choose the interval between the day of the earthquake and the first day of the precursor period. As we can see, for Greece this period is very short -1 day, and for Italy it is 5-6 days. For different regions, this interval can be from 10 to 1 day. Therefore, when making a forecast, a period of time typical for a given region is indicated, and is determined from multi-year observation data.

Determination of the magnitude and position of the epicenter of an earthquake
As the data of long-term researchers of various groups working in the field of monitoring of the physical precursors of earthquake show, the assessment of the earthquake preparation zone obtained back in 1979 [15] has been brilliantly confirmed in practice. Surprisingly, it was obtained in those times when remote sensing data did not exist, and an empirical estimate was made from ground-based point measurements and an estimate of the zone of elastic deformation at the level10 -8 . Using the Dobrovolsky formula, we can determine the magnitude of the future earthquake as M = [log (R)]/0.43, where R is the radius of the preparation zone in kilometers.
It turned out that using remote sensing data, this area can literally be seen. There is a technology for detecting surface thermal anomalies, which allows building a map of their distribution in real time [16]. In Fig. 7 are shown two similar distributions for two earthquakes with different magnitudes: the left panel is the earthquake in Aquila M6.3 on April 6, 2009, the right panel is the earthquake in Gujarat, India M7.7 on January 26, 2001 [1]. In the figure, the preparation zone is marked in blue according to the Dobrovolsky formula, R = 10 0.43M , and in red, the so-called activation zone according to [17]  As one can see, the anomalies practically do not go beyond the boundaries of the designated zones, although they do not completely fill the circle (which corresponds to the definition of a preparation zone). However, the disadvantage of this technology is that anomalies of this type are visible only in cloudless weather, which imposes serious restrictions on its application. At the same time, this is a great opportunity to visualize and prove the reality of the concept of an earthquake preparation zone.
As it turned out, the scale of ionospheric variations before earthquakes also has an order of magnitude determined by the Dobrovolsky formula. Maps of the distribution of ionospheric anomalies can be built both according to the data of local networks of stationary GPS/GLONASS receivers (in the case of a moderate earthquake magnitude), and according to the IONEX index, which present ready-made maps of the global distribution of the total electron content, but with low spatial resolution and restrictions, superimposed by heterogeneous distribution of IGS network receivers.
Despite the low resolution of IONEX maps (also called GPS GIM), today they are the most suitable option in terms of data availability and efficiency. IGS data in IONEX format is a matrix whose elements are TEC values multiplied by 10. The matrix resolution is 2.5 degrees in latitude and 5 degrees in longitude. TEC values are calculated by IGS every 2 hours (at present, the transition to a time resolution of 1 hour is underway). The calculation and construction of differential maps of the global TEC TECGIM, which represent the deviation of the current TECGIM of the TEC values from the background TECGIMA, is performed according to the formula: TECGIM = TECGIM -TECGIMA. The deviation from the background value is expressed in TECU units.
As an example, Fig. 8 shows a difference map obtained on March 8, 2011 (the maximum of ionospheric deviations before the Tohoku earthquake). For a magnitude of 9.1, the preparation zone radius will be ~ 8200 km. In fact, we are talking about global variation in the ionosphere, which occupies half of the earth's global longitude range. The third parameter that makes it possible to estimate the size of the earthquake preparation zone is the spatial distribution of the chemical potential [2]. In the Fig. 9 an example is shown from the publication [2] where the estimation of the earthquake preparation zone for the M6.4 earthquake on March 20, 2016 off the eastern coast of Kamchatka is given.
However, it is possible to assess the magnitude of an earthquake not by the geometric factors of the registered precursors, but based on the physical analysis of the dynamics of the precursor before the earthquake. This is the OLR heat flow. The rapid increase in radiation and a transient change in OLR were recorded at the top of the atmosphere over seismically active regions. In the first approximation, we can define the atmospheric anomaly in the Euler frame of reference by subtracting the mean value. The mean can be defined as the average for the same day of the year, local time, and location over more than 12 years (larger than 11 years, one solar cycle). The advantage of this approach is its simplicity and effectiveness with the availability of long-time satellite observations. Following this, the OLR anomalous variations were defined as an E_index [18].
Where: t=1, K days, S *(xij, yij ,t) is the current OLR and is the computed mean of the OLR field, defined for multiple years of observations over the same location and same local time, ti,j is the standard deviation.
The Magnitude estimation is based on the calculation of the speed of change of E_index. We use the Magnitude Assessment a modified version of E_Index (1), representing the regional calibration of E_Index estimates we apply for the Kamchatka region and defines as: Where A and B are regional calibration coefficients; A -is a mask, mainly defined by the regional seismo-tectonic patterns and frequency of appearance of OLR anomalies for the historical events and B -regional normalization of thermal flux energy, normalizes each of NOAA satellites 15 and 18 -time series of OLR data to the same time coverage (twelve years or more). The Kamchatka region A values vary around 0.5-0.7; B -range 1-3. The Magnitude assessment data have been computed with a resolution of 2.5° × 2.5°. The final step includes additional preprocessing to avoid aliasing short wavelengths and spatial filtering based on a "minimum curvature" algorithm [10].
At present, the accuracy of determining the magnitude is 0.5 on the Richter scale. Experience shows that determining the position of the epicenter of a future earthquake is the most difficult task. Actually, we use two methods for determining the position of the epicenter: finding the coordinates of the center of the identified earthquake preparation zone and the position of the OLR anomaly in its maximum development.
The distribution of the chemical potential correction shown in Fig. 9 is the exception rather than the rule, and this is natural. Not necessarily the maximum radon flux before an earthquake is observed at a point coinciding with the epicenter of a future earthquake. Then, not always observed anomalies both in the ionosphere and near the earth's surface are symmetrical with respect to the position of the epicenter. The low resolution of the IONEX index (2.5 in latitude and 5 in longitude) should also be noted, which greatly coarsens the accuracy of determining the position of the epicenter from global TEC maps.
As one can see from Fig. 4, the OLR anomaly drifts along active tectonic faults and tectonic plate boundaries, but its position of the OLR anomaly at its maximum value is always close to the epicenter position, and the deviation is no more than 2.5-3. In principle, the accuracy of determining the epicenter position from the data on the distribution of the chemical potential also lies within these limits.
Summing up, we can say that at the moment the accuracy of our short-term forecasting technology is: • Time of the earthquake (10-1 days) • Magnitude of an earthquake ± 0.5 on the Richter scale • Position of epicenter 2.5-3 in latitude and longitude It should be noted that OLR has two periods of leading time: middle term 20-5 days, and short term 5-1 days.

Creation of a model of automated monitoring and short-term forecasting of strong earthquakes in the Kamchatka region
To monitor the state of the atmosphere and ionosphere and identify the precursors of earthquakes in the Kamchatka region, it is required to continuously process data from selected receivers of signals from global navigation satellite systems (GNSS) and a groundbased vertical sounding ionosonde in Petropavlovsk-Kamchatsky simultaneously with data on solar and geomagnetic activity, data from monitoring atmospheric parameters , and only their joint analysis allows us to make an informed conclusion about the possibility of a seismic event in a certain place, at a certain time and with a certain magnitude [1]. The logistics of data processing automation can be conditionally divided into three stages: primary processing, data validation and intelligent processing, as a result of which, in fact, precursors are recognized. This division does not necessarily mean the same sequence in time (some operations can be performed simultaneously), but it helps to understand the process of processing itself.
In addition to ionospheric data for the analysis of the heliogeomagnetic conditions, it is planned to use data on the geomagnetic conditions -the Ap, Kp and Dst indices (https://www.gfz-potsdam.de; http://wdc.kugi.kyoto-u.ac.jp /dst_realtime/index.html) and solar radio flux data F10.7 (ftp://ftp.swpc.noaa.gov/pub/indices/old_indices/), freely available on the Internet and updated in real time.
Simultaneously with ionospheric data for the implementation of a synergistic approach in accordance with a complex model of lithosphere-atmospheric-ionosphere-magnetospheric coupling [11,19], it is proposed to use meteorological data (temperature and humidity of atmospheric air), for example, from the site https://meteoinfo.ru/ archive-pogoda, by which a generalized parameter is calculated, called the correction of the chemical potential of water vapor in the atmosphere [11]. If there are other sources of data on the impending earthquake in this region, for example, remote sensing data on thermal anomalies [18], they are also included in the process of intelligent processing.
The primary processing phase is the most expensive in terms of machine time, since it involves downloading large amounts of data from various sources simultaneously. Typically, data on source servers is stored in a packed form, so it needs to be unpacked, read in the source format, and converted to a format convenient for analysis. In this case, it is necessary to bring data from different sources to a common time scale: data from various geophysical indices can be provided with a time resolution of one hour, three hours, etc. Vertical sounding data can also be both with a resolution of 15 minutes and 1 hour, and etc. At the same time, TEC can be calculated both with a resolution of 15 s, 30 s, and with any given resolution at large time intervals.
Equally important is the validation phase. Unfortunately, not all receivers work flawlessly. There are whole periods when corrupted data comes in. The receiver may not work at all for some time, there may be outliers of values far beyond the permissible limits. Sometimes TEC processing algorithms give negative values, which is physically impossible. If all this is put into the recognition system, then we will face severe disappointment. Therefore, already at this stage, it is necessary to use sufficiently intelligent algorithms for recognizing failures in the data stream.
Only after bringing the data into proper form we can start machine processing. And, if automatic processing processes have been successfully used for a long time, then machine analysis systems for data containing complex patterns have become popular in geophysics only recently with the development of machine learning methods. At present, the accuracy of data analysis using artificial intelligence methods is not inferior to classical methods [20,21], while significantly exceeding them in terms of speed and usability. Among the numerous machine learning methods for recognizing anomalies, convolutional neural networks (CNNs) have found the greatest application. The application of these methods is widespread in geological exploration problems [22]. A number of similar machine learning methods have also been used to recognize ionospheric earthquake precursors based on two-dimensional TEC distribution maps [23].
For automated processing of ionospheric monitoring data and analysis of ionospheric precursors, a machine data processing system is proposed (see Fig. 10), which is described below.
Processing includes proven algorithms, including: 1. Analysis of the TEC (or foF2) data sets using the pattern recognition method  for the correspondence of the ionospheric precursor mask to the current changes in the ionosphere over the seismically active region [8].
2. Correlation analysis of arrays of daily TEC values (or critical frequency foF2) between a pair of adjacent GPS/GLONASS receivers (or ground stations for vertical radio sounding of the ionosphere) [14].
3. In the presence of a dense local network of stationary GPS/GLONASS receivers, the calculation of the coefficient of regional variability of the ionosphere can be provided [24]. 4. Calculation and construction of difference maps of the global TEC TECGIM in order to determine the position of the epicenter of the future earthquake and its magnitude) [1]. If there is a dense local network of fixed GPS/GLONASS receivers, differential maps can be calculated using local data rather than GPSGIM.
5. Comparison of variations of the global TEC [25] with the local TEC with reference to the solar activity index F10.7 [1].
6. Calculation of the correction of the chemical potential of water vapor according to the local temperature and relative humidity data to determine the time of a seismic event [11].
7. Construction of maps of distribution of the correction of the chemical potential according to the data of local temperature and relative humidity to determine the position of the epicenter of the future earthquake and estimate its magnitude [1].
8. In the case of low-latitude earthquakes, analysis of the dynamics of the equatorial anomaly (EA) in order to detect the absolute anomaly and the longitudinal effect in EA [1]. 9. Multiparameter analysis using operational data on other physical precursors, if any (radon activity, crustal conductivity, OLR, anomalous cloud structures) [1].
As a result of the analysis, the predicted values of the earthquake magnitude for a given region and their probability are estimated. Based on the forecasts, the final assessment of the probability of an earthquake is made based on machine learning data. This does not exclude the operator's expert judgment. Thus, the system carries out a multivariate analysis of the state of the ionosphere, capable of recognizing a unique image of an earthquake precursor.
More detailed description of the process of multiparameter data processing and interpretation is provided in [7].

Conclusion
The proposed approach is based on multi-years experience in analysing both ionospheric precursors of earthquakes and other physical precursors, which together create a generalized image of the final stage of preparation of strong earthquakes. The creation of the proposed system for processing ionospheric monitoring data and analysing ionospheric precursors will make it possible to create an earthquake prediction service capable of recognizing and identifying geophysical variations that are precursors of earthquakes in an automatic mode.
The authors thank NASA's Goddard Space Flight Center and UNAVCO for the GPS (Rinex) and GPS (IONEX) data, NOAA's Climate Prediction Center (USA), NASA Goddard Earth Science Center (GES