Analyzing the Performance of GPS Data for Earthquake Prediction

Abstract

The results of earthquake prediction largely depend on the quality of data and the methods of their joint processing. At present, for a number of regions, it is possible, in addition to data from earthquake catalogs, to use space geodesy data obtained with the help of GPS. The purpose of our study is to evaluate the efficiency of using the time series of displacements of the Earth’s surface according to GPS data for the systematic prediction of earthquakes. The criterion of efficiency is the probability of successful prediction of an earthquake with a limited size of the alarm zone. We use a machine learning method, namely the method of the minimum area of alarm, to predict earthquakes with a magnitude greater than 6.0 and a hypocenter depth of up to 60 km, which occurred from 2016 to 2020 in Japan, and earthquakes with a magnitude greater than 5.5. and a hypocenter depth of up to 60 km, which happened from 2013 to 2020 in California. For each region, we compare the following results: random forecast of earthquakes, forecast obtained with the field of spatial density of earthquake epicenters, forecast obtained with spatio-temporal fields based on GPS data, based on seismological data, and based on combined GPS data and seismological data. The results confirm the effectiveness of using GPS data for the systematic prediction of earthquakes.

1. Introduction

In the immediate vicinity of the source of a future earthquake, anomalous changes in a number of processes occur. Abnormalities are instrumentally recorded in the frequency and intensity of seismic events, in the deformations of the Earth’s surface, in the chemical composition of fluids, in the groundwater level, in the time of passage of seismic waves, in the values of electric and geomagnetic fields, etc. [1,2,3,4]. Seismological monitoring data are currently the most widely available. Therefore, systematic earthquake prediction systems use only seismological data [5,6]. The use of additional information on spatial and temporal changes in a geological environment can provide more accurate predictions of earthquakes.

In recent years, the monitoring data of the Earth’s surface displacement obtained by Global Positioning System (GPS) have been published in real time for a number of seismically active regions. These data are used to study block models of the Earth’s crust and in earthquake prediction studies. In the articles [7,8,9,10], the displacements of the Earth’s surface were estimated from triplets of GPS ground network stations. Triangles with vertices at the stations’ locations can intersect different plates. It was shown that before strong earthquakes, there are changes in the area and other parameters of the triangles. In a number of later works, artificial intelligence methods were used to study the effectiveness of earthquake predictions based on GPS data [11,12,13,14]. In [15], GPS ground displacement data were used to predict 15 major earthquakes in western North America in 2007–2016. The forecast algorithm uses time series of horizontal and vertical displacements of GPS stations.

It is known that the accuracy of GPS measurements for vertical components of coordinates is 2–4 times less than for horizontal ones [16,17]. Therefore, in this paper, we considered only the horizontal components of the displacement velocities. It was theoretically shown in [18] that the field of horizontal velocities of the Earth’s surface displacements reflects the movements of the crustal blocks. Analysis of GPS data for the Bishkek geodynamic test area in the Tien Shan showed that seismicity increases in those areas where the horizontal component of the rate of change in the size and/or rotation of the Earth’s surface is close to its maximum values [19].

In this article, we use a statistical approach to assess the effectiveness of using GPS data for earthquake prediction. Our conclusion is based on the results of earthquake forecasting in two regions, which differ significantly in the seismotectonic and geodynamic regimes. Following a statistical approach, we consider GPS data to be effective if the result of a successful prediction of earthquakes based on these data is significantly better than the result of a random prediction of the same earthquakes or a prediction based on spatial seismicity data.

Obviously, the forecast result depends on both the GPS time series preprocessing methods and the forecast method. Currently, a number of works are underway using machine learning methods for predicting earthquakes. They include the study of mathematical models for earthquake prediction, machine learning methods and testing earthquake prediction algorithms [20,21,22,23,24,25,26,27,28,29]. Many articles are devoted to earthquake prediction using artificial neural networks [30,31] and more sophisticated methods, including hybrid neural network [32], and recurrent neural network [33], and a lot of other machine learning methods [34]. Despite the relatively positive results claimed in those studies, some simpler models can offer similar or better predictive powers [35]. This is due to the fact that the number of major seismic events is very low, making the forecasting of such events more difficult [36].

Here we use a machine learning prediction method that we call the method of the minimum area of alarm [37]. Forecast training is carried out according to the precedents of the occurred target earthquakes. The learning-based minimum alarm zone method systematically calculates a predictive alarm zone map in which a target earthquake with a magnitude greater than a predetermined magnitude is expected. The alarm zone is determined by the function of the grid spatio-temporal fields of features calculated on the basis of the initial data. For many seismically active regions, the quality of machine learning for predicting strong earthquakes is limited by a small number of precedents. Our method uses a mathematical model that allows, during training, to select anomalies in feature fields that precede target earthquakes, compare these anomalies with each other, and select other anomalies in feature fields similar to those selected during training. This model introduces constraints on the prediction rule that help compensate to some extent for a small number of use cases during training.

This article is divided into three sections. In Section 2, we briefly discuss the method of the minimum area of alarm. Section 3 and Section 4 present and discuss the results of modeling the forecast of strong earthquakes in Japan and California. For each region, we compare the following results: random forecast of earthquakes, forecast obtained with the field of spatial density of earthquake epicenters, forecast obtained with spatio-temporal fields based on GPS data, based on seismological data, and based on combined GPS data and seismological data.

2. Method

The method of the minimum area of alarm was designed to systematically predict strong earthquakes with a magnitude above a certain threshold [37]. The platform for systematic earthquake prediction regularly with a step $\Delta t$ calculates the alarm zone, in which the epicenter of the target earthquake is expected. A demo version of the platform has been available since 2018 at https://distcomp.ru/geo/prognosis/ (accessed on 5 March 2021).

We assume that strong earthquakes are preceded by abnormal manifestations, which can be represented in a spatio-temporal grid-based fields of features. The method of the minimum area of alarm is trained to calculate alarm zones based on retrospective observations consisting of the spatio-temporal fields of features, ${\mathbf{F}}_i,i=1,\dots ,I$, and a sample consisting of marked target abnormal objects and an unmarked mixture of target and normal objects. The abnormal objects here are the field values selected by the machine learning algorithm among the points preceding the epicenters of the target earthquakes $q=1,\dots ,Q$ with magnitudes $m\ge M$ (earthquake precursors). An unlabeled mixture of objects is represented by points with field values at all other grid nodes. This formulation in machine learning refers to one-class classification problems [38,39,40].

The fields ${\mathbf{F}}_i$ are set in a single coordinate grid. The values of the fields at the grid nodes $n=1,\dots ,N$ correspond to the points of the I-dimensional feature space ${\mathbf{f}}^{\left(n\right)}\in {\mathbb{R}}^I$. During training, the spatio-temporal field $\mathbf{\Phi }$, called the alarm field, is calculated. The values of the alarm field ${\varphi }^{\left(n\right)}\ge {\varphi }_0$ allocate a spatio-temporal area in it, called the alarm area. The temporal slice of the alarm area at the time of the forecast $t^{*}$ is exactly the alarm zone$S\left(t^{*}\right)$, in which the appearance of the epicenter of the target earthquake is predicted at the interval $[t^{*},t^{*}+\Delta t)$.

Figure 1 shows operations in the interval $[t^{*}-\Delta t,t^{*})$. The program downloads new data from remote servers and user computer, then preprocess them, calculates spatial and spatio-temporal fields of features, and selects data on target earthquakes from the catalog. Then, using the algorithm of the method of the minimum area of alarm, the program calculates the alarm field $\mathbf{\Phi }$ in the interval from the beginning of training to the moment of forecasting $t^{*}$. The time interval of the alarm field at the moment $t^{*}$ is a map of the alarm zone.

The target earthquake is predicted at the step $[t^{*},t^{*}+\Delta t)$ if its epicenter falls into the alarm zone. The larger the product $S\left(t^{*}\right)\Delta t$, the more successful the forecast. However, it is obvious that the value of this spatio-temporal area should be reasonably limited. The indicators of the forecast quality are the estimations of the probability of a successful forecast of events (the forecast probability) and the size of the alarm area in the training interval (volume of alarm). As a result of training, it is desirable to obtain a solution in which the maximum number of successful forecasts of target events is achieved for a given size of the volume of alarm.

The method of the minimum area of alarm is based on a model that introduces a measure of the anomality for points in the feature space. The algorithm uses this measure to select points that can be considered precursors of target earthquakes. The model contains two assumptions:

• Abnormality: In feature space, the target earthquakes are preceded by points (earthquake precursors) for which the values of some components (the values of some feature fields) are unlikely and close to the maximum (or minimum). To simplify the explanation, we assumed without loss of generality that the precursors refer only to the maximal values of the fields of features.
• Monotonicity: Feature space points that are, component-wise, larger than the earthquake precursor can also be precursors of similar target events.

For training, the earthquake prediction algorithm requires determining the precursor for each target event. Let the event q be preceded by a precursor ${\mathbf{f}}^{\left(q\right)}\in {\mathbb{R}}^I$. The precursor is associated with a feature space area whose points are, component-wise, greater than or equal to ${\mathbf{f}}^{\left(q\right)}$, i.e., the area $w^{\left(q\right)}=\{{\mathbf{f}}^{\left(n\right)}\in {\mathbb{R}}^I:f_i^{\left(n\right)}\ge f_i^{\left(q\right)}, i=1,2,\dots ,I\}$. We called this area the orthant $w^{\left(q\right)}$ with a vertex at ${\mathbf{f}}^{\left(q\right)}$ and the points ${\mathbf{f}}^{\left(n\right)}\in w^{\left(q\right)}$ are called the base points of the target event q.

According to the monotonicity condition, base points are also precursors of events that are similar to the event q. In geographical coordinates, each base point forms a cylinder of alarm with the radius R and the element T. The alarm cylinder of the base point ${\mathbf{f}}^{\left(n\right)}$ has a circular base center at the grid node n with coordinates $(x^{\left(n\right)},y^{\left(n\right)},t^{\left(n\right)})$, a radius of the base R, and the elements $[(x^{\left(n\right)},y^{\left(n\right)},t^{\left(n\right)}),(x^{\left(n\right)},y^{\left(n\right)},t^{\left(n\right)}+T)]$, where $t^{\left(n\right)}\in [t_0,t^{*}]$, in which $t_0$ is the start time of training and $t^{*}$ is the time of action on forecast. An earthquake can only be predicted if its epicenter falls into one of the alarm cylinders. The union of alarm cylinders formed by the base points of the orthant $w^{\left(q\right)}$ selects a set of grid nodes $W^{\left(q\right)}$, $\left|W^{\left(q\right)}\right|=L^{\left(q\right)}$. Accordingly, it follows that an earthquake q with the epicenter coordinates $(x^{\left(q\right)},y^{\left(q\right)},t^{\left(q\right)})$ can be predicted only if the cylinder with the center of the base of $(x^{\left(q\right)},y^{\left(q\right)},t^{\left(q\right)})$, a radius R, and the elements $[(x^{\left(q\right)},y^{\left(q\right)},t^{\left(q\right)}-T),(x^{\left(q\right)},y^{\left(q\right)},t^{\left(q\right)})]$ contains at least one base point. This cylinder is called a precursor cylinder. The precursor of event q is the point ${\mathbf{f}}^{\left(q\right)}\in {\mathbb{R}}^I$, which has the minimum value of the alarm volume $v^{\left(q\right)}=L^{\left(q\right)}/L$ among all points, corresponding to the grid nodes of the precursor cylinder, where L is the number of all grid nodes of the spatio-temporal area of analysis. The quantity $v^{\left(q\right)}$ (volume of the precursor) defines the measure of anomality for the precursor of event q. In the paper [41], the measure of the abnormality for target events was determined by the likelihood ratio estimate.

The forecast quality is determined by two indicators: (1) The probability of prediction U, equal to the fraction of correctly forecasted target events $Q^{*}$ to all target Q events, $U=Q^{*}/Q$, and (2) the volume of alarm V, equal to the fraction of the number of grid nodes of the alarm field $L^{*}$ with the values ${\varphi }^{\left(n\right)}\ge \varphi$ to the number of all grid nodes of the analyzed area L, $V=L^{*}/L$. It can be seen from the definition that the alarm volume V is equal to the probability of the detection of target events by random areas consisting of $L^{*}=VL$ grid nodes. Usually, the quality of the forecast is determined by the dependence $U\left(V\right)$, which is analogous to the error curve represented by the ROC (receiver operating characteristic) curve [42,43,44].

In classification problems, the decision rule is found by minimizing the function of losses from target detection errors and false alarms. The training algorithm is optimal if it calculates an alarm area with the volume $V_0$, which, for any value $V\le V_0$, provides the maximum value of U. In our case, this solution requires large calculations because the sets of base points associated with different precursors intersect. Therefore, we considered solutions that are close to optimal [37].

Consider the algorithm for constructing the alarm field.

• Generate a training sample set $\{{\mathbf{f}}^{\left(q\right)},v^{\left(q\right)}\}$. Arrange the precedents ${\mathbf{f}}^{\left(q\right)}$, $q=1,\dots ,Q$, in accordance with the increase in the alarm volume of the target events $v^{\left(1\right)}\le v^{\left(2\right)}\le \cdots \le v^{\left(q\right)}\le \cdots \le v^{\left(Q\right)}$.
• Calculate the alarm field $\mathbf{\Phi }$.

  (a)

  Assign to the nodes of the grid of the alarm field $\mathbf{\Phi }$ a value of 1.

  (b)

  Replace the value of 1 of the field $\mathbf{\Phi }$ with $V\left(1\right)=v^{\left(1\right)}$ in the set $W^{\left(1\right)}$ of the grid nodes corresponding to the base points of the precedent ${\mathbf{f}}^{\left(1\right)}$; replace the value of 1 with $V\left(2\right)=|W^{\left(1\right)}\cup W^{\left(2\right)}|/L$ at the grid nodes $W^{\left(2\right)}\backslash W^{\left(1\right)}$; replace the value of 1 with $V\left(3\right)=|W^{\left(1\right)}\cup W^{\left(2\right)}\cup W^{\left(3\right)}|/L$ at the grid nodes $W^{\left(3\right)}\backslash (W^{\left(1\right)}\cup W^{\left(2\right)})$ and then sequentially replace the values of 1 with the corresponding values of the alarm volumes. The resulting field $\mathbf{\Phi }$ takes the values $V\left(1\right)\le V\left(2\right)\le \cdots \le V\left(q\right)\le \cdots \le V\left(Q\right)$ or 1. The value $V\left(1\right)$ refers to grid nodes from the set $W^{\left(1\right)}$, $V\left(2\right)$ refers to grid nodes from the set $(W^{\left(2\right)}\backslash W^{\left(1\right)})$, $V\left(3\right)$ refers to grid nodes from the set $(W^{\left(3\right)}\backslash (W^{\left(1\right)}\cup W^{\left(2\right)}))$, etc.

Figure 2 shows a flowchart of the algorithm. The grid feature fields and the sample of target earthquakes are regularly updated with a step $\Delta t$. For each epicenter from the sample, the precursor cylinders are calculated. Then, among the points of the feature space corresponding to the grid nodes in these cylinders, precursors ${\mathbf{f}}^{\left(q\right)}$, orthants $w^{\left(q\right)}$, and alarm volumes of precursors $v^{\left(q\right)}$ are selected. Then the precursors and orthants are ordered according to the increase in the alarm volumes of the precursors. After that, the alarm area and the map of the alarm zone are calculated.

3. Modeling

3.1. Technology

The efficiency of using the information on the Earth’s surface displacements according to GPS data for predicting earthquakes largely depends on the methods of data processing and analysis. For earthquake prediction, we used the method of the minimum area of alarm. Modeling simulates the operation of a demo platform (https://distcomp.ru/geo/prognosis/) (accessed on 5 March 2021). We performed modeling on the GIS GeoTime 3 (http://geo.iitp.ru/GT3/) (accessed on 5 March 2021).

The modeling technology, similarly to the forecast, is as follows. The area of analysis is selected. Point fields, time series, and rasters can be used as input. All initial data are converted into uniform grid-based spatial and spatio-temporal fields, which are used to calculate the zones of expected epicenters of target earthquakes at each moment of the forecast. The forecast is carried out systematically with a constant step. At each step n, training is performed using the data available from the start of training to the moment of forecasting $t^{*}$, and the decision rule computed during training is checked at the interval $[t^{*},t^{*}+\Delta t)$. In the process of training, according to the monitoring data, new values of the forecast feature grid fields are calculated, the training sample of target earthquake epicenters is supplemented, the alarm area $\mathbf{\Phi }\le V_0$ and the alarm zone $S\left(t^{*}\right)$ are calculated. Testing verifies whether new epicenters of target earthquakes are in the calculated alarm zone.

3.2. Area of Analysis

The regions under study are in Japan within the boundaries of 130–145.5° E and 30–43° N, and California (USA) within the boundaries of 125–114° W and 32–43° N. We forecasted earthquakes with a hypocenter depth of $H\le 60$ km for both regions, but with different magnitudes: $m\ge 6.0$ in Japan and $m\ge 5.5$ in California. The analysis of the depths of the hypocenters characteristic of the foci of tectonic earthquakes and the zones of their preparation in each of the regions under consideration was not carried out in this work. The area of analysis was determined by the intersection of two areas, A and B. Area A was defined as the territory for which the distance from any point to GPS ground receiving stations does not exceed 50 km. Area B was calculated from seismological data. It was selected such that, with an alarm volume of $V=0.2$, the probability of predicting target earthquakes from the field of spatial (2D) density of earthquake epicenters $U_{2D}\approx$ 0.3–0.4. For this, in the time interval before the start of training, the 2D density field of earthquake epicenters $\mathbf{S}$ with an averaging radius of 50 km was calculated. The boundary of the area of analysis was determined by the values of the field $s_n\ge C$. If the forecast probability for the threshold C was greater than 0.3–0.4, then the threshold C had increased, and the procedure was repeated. Furthermore, the area of analysis was found as the intersection of areas A and B.

3.3. GPS Data Preprocessing

We analyzed the time series of daily horizontal displacements of the Earth’s surface at the intervals 1 January 2009–26 July 2020 for Japan and 1 January 2008–14 November 2020 for California. The data were obtained from the Nevada Geodetic Laboratory (NGL), http://geodesy.unr.edu/about.php (accessed on 5 March 2021) [45]. There are 1420 and 1803 GPS receiving stations in Japan and California, respectively, and the analysis areas contain 1229 and 1204 stations, respectively. The networks of GPS receiving stations, the areas of analysis, and the epicenters of target earthquakes are shown in Figure 3. The stations evenly cover the analysis area. The average minimum distance between stations is 12.8 km for Japan and 9.38 km for California, with standard deviations of 5.4 and 5.74 km, respectively.

The calculation of the feature fields used for the earthquake prediction based on GPS data was performed in two stages. The purpose of the first stage was to extract a useful signal from the time series of coordinates of the receiving stations. The purpose of the second stage was to calculate spatio-temporal fields of forecast features.

3.3.1. Time Series of the Earth’s Surface Displacement Velocities

The initial data were daily time series of coordinates $x\left(t\right)$ and $y\left(t\right)$ of GPS ground receiving stations in the W–E and N–S directions in the intervals 1 January 2009–26 July 2020 for Japan and 1 January 2008–14 November 2020 for California. The daily horizontal velocities of the Earth’s surface displacements $g_x\left(t\right)$ and $g_y\left(t\right)$ were determined by two coordinates of the GPS receiving station, spaced in time by the interval $T_0$: $g_x\left(t\right)=(x\left(t\right)-x(t-T_0))/T_0$, $g_y\left(t\right)=(y\left(t\right)-y(t-T_0))/T_0$. There were discontinuities (gaps) in the time series. In our case, for each coordinate, there were 23,682 gaps and 168,218 days of missed measurements for Japan, and 29,977 gaps and 357,067 days of missed measurements for California. Since the displacement rate estimates are ahead of the time of the values of the first station coordinates $x(t-T_0)$ and $y(t-T_0)$ by $T_0$ days, each gap in the time series of the station coordinates increases the number of missing values in the time series of velocities by $T_0$ days. With a large number of gaps, the number of missing velocity values can significantly exceed the number of missing coordinate values. To limit the number of missed velocity values, we linearly interpolated the coordinate values in the gaps less than or equal to $T_0$. For gaps of more than $T_0$ days, we ended the calculation of the speed at the last value of the station coordinate before the start of the rupture and re-estimated the rates, starting from the first value of the station coordinate after the rupture.

To calculate the daily rates, the interval $T_0=30$ days was selected. For this interval, the station movement values were comparable to the noise values of the daily measurements. For 30-day intervals, there were 23,119 gaps in each coordinate of the area in Japan, which was 97.62% of all gaps in the time series, and 50,537 blanks in the measurements (30.04%). For California, there were 27,990 measurement gaps (93.37%) and 88,387 blank measurements (24.75%) over a 30-day period for each coordinate. At the same time, the number of missing speed values in each of the W–E and N–S directions increased for Japan by (23,682 − 23,119) $\times 30=16,890(10.04\%)$, and for California by (29,977 − 27,990) × 30 = 59,610 (16.69%).

The first stage was completed by calculating the spatio-temporal fields of the rate components ${\mathbf{V}}_x$ and ${\mathbf{V}}_y$ in the W–E and N–S directions. The fields were presented in the grid $\Delta x\times \Delta y\times \Delta t=0.1^{\circ }\times 0.{075}^{\circ }\times 1$ day. The calculation of the fields was carried out using an interpolation technique known as inverse distance weighting. During interpolation, the gaps in the values of the time series were not filled, but they were taken into account as the absence of the receiving station. The values of the fields ${\mathbf{V}}_x$ and ${\mathbf{V}}_y$ at the grid points for each time slice of the field of the velocity component W–E were calculated by the formula:

$$V_{xn}\left(t\right)=\frac{{\sum }_{k=1}^K{g}_{xk}^{\left(t\right)}/r_k^p}{{\sum }_{k=1}^{K}1/r_k^p},$$

(1)

where $V_{xn}\left(t\right)$ is the value of the field of the W–E strain rate component at the grid node n at the moment t, K is the maximum number of stations closest to the node n in the circle of radius $R_{max}$, the values of which were used for interpolation, $g_{xk}\left(t\right)$ is the value of the W–E strain rate component for the station $k,k=1,\dots ,K$, at the time t, $r_k\le R_{max}$ is the distance from the k-th station to the grid node n, and p is the degree that determines the dependence of the station weight on its distance to the grid node. The interpolation parameters were $K=5$, $R_{max}=50$ km, and $p=1$. If $r_k=0$, then $V_{xn}\left(t\right)=g_{xk}\left(t\right)$. The calculations of the field of the N–S strain rate components were similar.

Figure 4 and Figure 5 show an example of the time series at the point with coordinates $35.{45}^{\circ }$N and $139.2^{\circ }$E for Japan: (A) coordinates of the receiving station, (B) daily velocities, and (C) values of the field components ${\mathbf{V}}_x$ and ${\mathbf{V}}_y$. The time series in the figure on the left refer to the interval 1 January 2009–26 July 2020, and on the right, to the interval 1 January 2016–14 December 2016. The jump in the time series of the station coordinates and the ejection of the velocity series was caused by the Tohoku earthquake (11 April 2011). It can be seen that the gaps in the values of the coordinates of the stations cause discontinuities in the time series of the daily rates of the displacement of the Earth’s surface. However, when calculating ${\mathbf{V}}_x$ and ${\mathbf{V}}_y$ strain rate fields, the missing velocity values were filled in by interpolating data from other GPS receiving stations. Figure 6 and Figure 7 show similar example of the time series at the point with coordinates $39.{49}^{\circ }$N and $123.2^{\circ }$W for California for the intervals 1 January 2008–14 November 2020 and 1 January 2013–1 October 2013.

3.3.2. Spatio-Temporal Fields of Features

We assume that strong earthquakes are preceded by spatio-temporal anomalous changes in the regime of various deformations of the Earth’s surface. Therefore, we looked for fields containing information about the anomalous values of the change in the deformation mode. The basis of the considered fields of features was the following invariants of the strain rate fields.

• ${\mathbf{F}}_1$ is the field of divergence of the strain rates:

  $$\mathrm{div}{V}_n=\frac{\partial V_{xn}}{\partial x}+\frac{\partial V_{yn}}{\partial y}$$

  (2)

The maximum and minimum values of the divergence field refer to places where there is a relative contraction or expansion of the size of a small horizontal area.

• ${\mathbf{F}}_2$ is the field of rotor of the strain rates:

  $$\mathrm{rot}{V}_n=\frac{\partial V_{xn}}{\partial y}-\frac{\partial V_{yn}}{\partial x}$$

  (3)

The field values determine the direction and intensity of the field twisting around the vertical axis.

• ${\mathbf{F}}_3$ is the field of shear of the strain rates:

  $$\mathrm{sh}{V}_n=\frac{1}{2}\sqrt{{\left(\frac{\partial V_{xn}}{\partial x}-\frac{\partial V_{yn}}{\partial y}\right)}^2+{\left(\frac{\partial V_{xn}}{\partial y}+\frac{\partial V_{yn}}{\partial x}\right)}^2}$$

  (4)

Figure 8 and Figure 9 show a time series of fields ${\mathbf{F}}_1$, ${\mathbf{F}}_2$, and ${\mathbf{F}}_3$ at the same points as in Figure 4, Figure 5, Figure 6 and Figure 7.

The fields of features ${\mathbf{F}}_4$, ${\mathbf{F}}_5$, and ${\mathbf{F}}_6$ represent changes in the fields of the strain rate invariants over time. They are equal to the ratios of the mean values of the invariants in two consecutive intervals to the standard deviation of this difference. The values of the fields are converted into the grid $\Delta x\times \Delta y\times \Delta t=0.1^{\circ }\times 0.{075}^{\circ }\times 30$ days.

• ${\mathbf{F}}_4$ is the field of the temporal variations in the divergence strain rate.

The value of the field $f_{4n}\left(t\right)$ at time t is equal to the ratio of the difference $(\overline{{\mathrm{div}}_{2n}}-\overline{{\mathrm{div}}_{1n}})$ between the mean values of the divergence in two consecutive intervals, namely, $T_1$ and $T_2$, to the standard deviation of this difference ${\sigma }_n\left(\mathrm{div}\right)$, $T_1=T_2=361$ days.

$$f_{4n}\left(t\right)=(\overline{{\mathrm{div}}_{2n}}-\overline{{\mathrm{div}}_{1n}})/{\sigma }_n\left(\mathrm{div}\right),$$

(5)

where $\overline{{\mathrm{div}}_{2n}}$ is calculated from the values of field ${\mathbf{F}}_1$ at the interval $(t-T_2,t)$, $\overline{{\mathrm{div}}_1}$ is calculated at the interval $(t-T_2-T_1,t-T_2)$.

• ${\mathbf{F}}_5$ is the field of the temporal variations in the rotor rate.

The values of field $f_5\left(t\right)$ are calculated similarly to the values of field ${\mathbf{F}}_4$,

$$f_{5n}\left(t\right)=(\overline{{\mathrm{rot}}_{2n}}-\overline{{\mathrm{rot}}_{1n}})/{\sigma }_n\left(\mathrm{rot}\right).$$

(6)

• ${\mathbf{F}}_6$ is the field of the temporal variations in the shear deformation rate.

The values of field $f_6\left(t\right)$ are calculated similarly to the values of field ${\mathbf{F}}_4$,

$$f_{6n}\left(t\right)=(\overline{{\mathrm{sh}}_{2n}}-\overline{{\mathrm{sh}}_{1n}})/{\sigma }_n\left(\mathrm{sh}\right),$$

(7)

Figure 10 and Figure 11 show a time series of fields ${\mathbf{F}}_4$, ${\mathbf{F}}_5$, and ${\mathbf{F}}_6$ at the same points as in Figure 4, Figure 5, Figure 6 and Figure 7.

Fields ${\mathbf{F}}_7$, ${\mathbf{F}}_8$, and ${\mathbf{F}}_9$ represent spatial correlations of strain rate changes in a sliding window of $75\times 75$ km${}^2$. With this window size, the correlation coefficients are estimated in approximately 70–80 grid points of the fields.

• ${\mathbf{F}}_7$ is the field of spatial correlations in fields ${\mathbf{F}}_4$ and ${\mathbf{F}}_5$.
• ${\mathbf{F}}_8$ is the field of spatial correlations in fields ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$.
• ${\mathbf{F}}_9$ is the field of spatial correlations in fields ${\mathbf{F}}_5$ and ${\mathbf{F}}_6$.

The time series of the correlation fields at the same points as in Figure 4, Figure 5, Figure 6 and Figure 7 are shown in Figure 12 and Figure 13.

Correlation fields ${\mathbf{F}}_7$, ${\mathbf{F}}_8$, and ${\mathbf{F}}_9$ carry information about the spatial relationship between the values of the change in the rate of different pairs of deformation types. The minimum or maximum fields of the correlation fields combine this information.

To combine information on the spatial relationship between the values of the change in the rate of different pairs of deformation types, you can use the field of maximum or minimum values of the correlation fields ${\mathbf{F}}_7$, ${\mathbf{F}}_8$, and ${\mathbf{F}}_9$. This operation can be interpreted in terms of fuzzy logic [46]. For Japan, the most successful prediction of target earthquakes was obtained using the ${\mathbf{F}}_{10}$ field:

• ${\mathbf{F}}_{10}$ is the field of minimum values of the fields ${\mathbf{F}}_7$ and ${\mathbf{F}}_9$:

  $$f_{10n}\left(t\right)=min(f_{7n}\left(t\right),f_{9n}\left(t\right))$$

  (8)

The time series of field ${\mathbf{F}}_{10}$ at point $35.{45}^{\circ }$N and $139.2^{\circ }$E in Japan and $39.{49}^{\circ }$N and $123.2^{\circ }$W in California are shown in Figure 14 and Figure 15.

3.4. Seismological Data Preprocessing

Seismological data for Japan and California were taken from the Japan Meteorological Agency earthquake catalogs [47,48] and the National Earthquake Information Center (NEIC) [49] at intervals 02 June 2002–26 July 2020 and 1 January 1995–20 December 2020, represented by earthquakes with a magnitude of $m\ge 2.0$ and a hypocenter depth of $H\le 160$ km.

• ${\mathbf{S}}_1$ is the 3D field of the density of earthquake epicenters.
• ${\mathbf{S}}_2$ is the 3D field of the mean earthquake magnitude.
• ${\mathbf{S}}_3$ is the 3D field of the negative temporal anomalies of the density of earthquake epicenters.
• ${\mathbf{S}}_4$ is the 3D field of the positive temporal anomalies of the density of earthquake epicenters.
• ${\mathbf{S}}_5$ is the 3D field of the negative temporal anomalies of the mean earthquake magnitude.
• ${\mathbf{S}}_6$ is the 3D field of the positive temporal anomalies of the mean earthquake magnitude.
• ${\mathbf{S}}_7$ is the 2D field of the density of earthquake epicenters: Kernel smoothing with the parameter $R=50$ km in the interval from the beginning of the analysis to the start of training.
• ${\mathbf{S}}_8$ is the 3D field of quantiles of the background density of earthquake epicenters, calculated using the interval from the beginning of the analysis to the start of training, which corresponds to the density values of earthquake epicenters at the current time.

The estimation of 3D fields ${\mathbf{S}}_1$ and ${\mathbf{S}}_2$ was performed with the method of local kernel regression. The kernel function for the q-th earthquake has the form $K_q={\left[cos{h}^2{(r_q/R)}^{2}cos{h}^2(t_q/T)\right]}^{-1}$, where $r_q<Rϵ$ and $t_q<Tϵ$ are the distance and time interval between the q-th epicenter of the earthquake and the node of the 3D grid of the field, $ϵ=2$, $R=50$ km, $T=100$ days for ${\mathbf{S}}_1$ and R = 100 km, and T = 730 days for ${\mathbf{S}}_2$. The field ${\mathbf{S}}_7$ was calculated with the kernel function $K_q={\left[cos{h}^2{(r_q/R)}^2\right]}^{-1}$. The parameters for evaluating the fields, the radii R, and the interval T were chosen empirically, considering the step of the network of fields and the approximate number of events in the evaluation window. To calculate the fields ${\mathbf{S}}_3$, ${\mathbf{S}}_4$, ${\mathbf{S}}_5$, and ${\mathbf{S}}_6$, Student’s t-test was used. This t-statistic was determined for each grid node as the ratio of the difference in the average values of the current interval $T_2$ and the background interval $T_1$ to the standard deviation of this difference. Positive values of the t-statistic correspond to an increase in the value on the test interval.

We also analyzed fields similar to the fields ${\mathbf{S}}_3$, ${\mathbf{S}}_4$, ${\mathbf{S}}_5$, and ${\mathbf{S}}_6$, but with different values of the $T_2$ and $T_1$ intervals.

When predicting from seismological data, the following three fields turned out to be the most informative.

• ${\mathbf{S}}_9={\mathbf{S}}_1/({\mathbf{S}}_8+0.001)$ is the field of ratios of the density values of the earthquake epicenters $s_{1n}$ to the values of the quantiles of the density of the epicenters calculated on the interval from the beginning of the analysis to the start of training, which corresponds to the density values of earthquake epicenters at the current time $(s_{8n}+0.001)$.
• ${\mathbf{S}}_{10}$ is the field of negative anomalies of Student’s t-statistic of the density of earthquake epicenters with the intervals $T_1=1095$ and $T_2=365$ days.
• ${\mathbf{S}}_{11}$ is the field of negative anomalies of Student’s t-statistic of the mean earthquake magnitude with the intervals $T_1=1095$ and $T_2=730$ days.

For Japan, the most informative were the fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{10}$. Both of them previously proved to be the most effective in predicting earthquakes and their magnitudes in Kamchatka and the Aegean region [50]. The anomalous values of the ${\mathbf{S}}_9$ field correspond to areas of the seismic process in which the density values of earthquake epicenters are quite high but significantly less than the average values of the density of epicenters in the interval from the beginning of the analysis to the start of training. The anomalous values of the ${\mathbf{S}}_{10}$ field correspond to the spatio-temporal regions of the seismic process, in which the average values of the density of earthquake epicenters in the $T_2$ interval are significantly lower than the average field values in the $T_1$ interval. These changes highlight anomalous areas in which a quiescence sets in after the activation of the seismic process. The time series of the ${\mathbf{S}}_{10}$ field simulates the preparation of strong earthquakes proposed by the IPE model (named after the Institute of Physics of the Earth) proposed in [51]. For California, the most informative were the fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{11}$.

3.5. Earthquake Forecast

The initial training intervals started for Japan and California on 1 January 2011 and 23 December 2009 and ended before the next forecast, starting on 20 November 2015 and 19 January 2013, respectively. The testing intervals were 20 November 2015–10 September 2020 for Japan and 19 January 2013–14 November 2020 for California. The areas of analysis at the testing intervals contained 14 epicenters of target earthquakes in Japan with a magnitude of $m\ge 6.0$ and 12 epicenters of target earthquakes in California with a magnitude of $m\ge 5.5$. The alarm cylinder parameters were the cylinder radius $R=16$ km and its generatrix $T=91$ days for Japan and $R=18$ km and its generatrix $T=61$ days for California. The best forecast of target earthquakes according to GPS data for Japan was obtained from the ${\mathbf{F}}_{10}$ field and for California from the ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$ fields.

Table 1. Japan. Results of forecasting earthquakes with magnitudes $m\ge 6.0$ from field ${\mathbf{F}}_{10}$ and field ${\mathbf{S}}_7$.

N	Date	Time	Long.	Lat.	Earthquake Magnitude	Alarm Volume for Prediction from Field ${\mathbf{F}}_{10}$	Alarm Volume for Prediction from Field ${\mathbf{S}}_7$
1	14 January 2016	06:25:33	142.80	41.97	6.7	0.46	0.24
2	14 April 2016	15:26:34	130.81	32.74	6.5	0.13	0.50
3	14 April 2016	18:03:46	130.78	32.70	6.4	0.13	0.50
4	15 April 2016	19:25:05	130.76	32.76	7.3	0.13	0.50
5	21 October 2016	08:07:22	133.86	35.38	6.6	0.32	1.00
6	24 November 2016	00:23:36	141.35	37.18	6.2	0.69	0.02
7	28 December 2016	15:38:49	140.57	36.72	6.3	0.04	0.02
8	8 April 2018	19:32:30	132.59	35.19	6.1	0.20	0.86
9	18 June 2018	01:58:34	135.62	34.84	6.1	0.45	0.76
10	7 July 2018	14:23:48	140.59	35.17	6.0	0.15	0.16
11	8 January 2019	15:39:30	131.17	30.57	6.0	0.10	0.28
12	18 June 2019	16:22:19	139.48	38.61	6.7	0.68	1.00
13	19 April 2020	23:39:05	142.10	38.89	6.2	0.20	0.03
14	24 June 2020	23:39:05	141.11	35.55	6.1	0.77	0.04

and

Table 2. California. Results of forecasting earthquakes with magnitudes $m\ge 5.5$ from fields ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$ and field ${\mathbf{S}}_7$.

N	Date	Time	Long.	Lat.	Earthquake Magnitude	Alarm Volume for Prediction from Fields ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$	Alarm Volume for Prediction from Field ${\mathbf{S}}_7$
1	24 May 2013	03:47:08	121.06	40.19	5.7	0.18	0.84
2	24 August 2014	10:20:44	122.31	38.21	6.0	0.17	0.62
3	8 December 2016	08:18:00	118.90	38.38	5.6	0.10	0.48
4	28 December 2016	08:22:12	118.90	38.40	5.6	0.10	0.48
5	28 December 2016	09:13:47	118.90	38.38	5.5	0.10	0.48
6	23 June 2019	03:53:02	124.30	40.28	5.6	0.47	0.30
7	4 July 2019	17:33:49	117.50	35.70	6.4	0.10	0.04
8	6 July 2019	03:19:53	117.60	35.77	7.1	0.10	0.04
9	6 July 2019	03:47:53	117.75	35.90	5.5	0.10	0.04
10	15 May 2020	11:03:27	117.85	38.17	6.5	0.03	0.29
11	4 June 2020	01:32:11	117.43	35.61	5.5	0.48	0.04
12	24 June 2020	17:40:49	117.98	36.45	5.8	0.03	0.29

show the forecast results for each target event in Japan and California. The forecast is considered successful if the epicenter of the target event falls into the alarm zone with the alarm volume value $V\le 0.20$. The final forecast results are summarized in

Table 3. Japan. Probabilities of predicting earthquakes with magnitudes $m\ge 6.0$.

	Feature Field ${\mathbf{F}}_{10}$		Spatial Field of Epicenter Density ${\mathbf{S}}_7$
Alarm Volume	Number of Successful Predictions	Probability of Successful Predictions	Number of Successful Predictions	Probability of Successful Predictions
0.01	0	0.00	0	0.00
0.05	1	0.07	4	0.29
0.10	2	0.14	4	0.29
0.15	5	0.36	4	0.29
0.20	8	0.57	5	0.36

and

Table 4. California. Probabilities of predicting earthquakes with magnitudes $m\ge 5.5$.

	Feature Fields ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$		Spatial Field of Epicenter Density ${\mathbf{S}}_7$
Alarm Volume	Number of Successful Predictions	Probability of Successful Predictions	Number of Successful Predictions	Probability of Successful Predictions
0.01	0	0.00	0	0.00
0.05	2	0.17	4	0.33
0.10	2	0.17	4	0.33
0.15	8	0.67	4	0.33
0.20	10	0.83	4	0.33

. It can be seen that the forecast probabilities from GPS data with an alarm volume $V\le 0.20$ are significantly higher than the forecast probability for a random field and a spatial density field of earthquake epicenters ${\mathbf{S}}_7$.

Let us consider the prediction of target earthquakes based on feature fields calculated from the earthquake catalog. The alarm fields for Japan and California were calculated in the same grid and with the same parameters of the alarm cylinders as for the forecasted results from the GPS fields.

Two fields of features were selected for Japan: Field ${\mathbf{S}}_9={\mathbf{S}}_1/({\mathbf{S}}_8+0.001)$ and field ${\mathbf{S}}_{10}$. The results of the forecasting based on seismological feature fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{10}$, as well as based on fields ${\mathbf{S}}_9$, ${\mathbf{S}}_{10}$, and ${\mathbf{F}}_{10}$ are shown in

Table 5. Japan. Results of forecasting earthquakes with magnitudes $m\ge 6.0$ from seismological fields and from the same fields together with the GPS field.

N	Date	Time	Long.	Lat.	Earthquake Magnitude	Alarm Volume for Prediction from Fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{10}$	Alarm Volume for Prediction from Fields ${\mathbf{S}}_9$, ${\mathbf{S}}_{10}$, and ${\mathbf{F}}_{10}$
1	14 January 2016	06:25:33	142.80	41.97	6.7	0.20	0.24
2	14 April 2016	15:26:34	130.81	32.74	6.5	0.17	0.15
3	14 April 2016	18:03:46	130.78	32.70	6.4	0.17	0.15
4	15 April 2016	19:25:05	130.76	32.76	7.3	0.17	0.15
5	21 October 2016	08:07:22	133.86	35.38	6.6	1.00	1.00
6	24 November 2016	00:23:36	141.35	37.18	6.2	0.01	0.01
7	28 December 2016	15:38:49	140.57	36.72	6.3	0.04	0.01
8	8 April 2018	19:32:30	132.59	35.19	6.1	0.84	0.38
9	18 June 2018	01:58:34	135.62	34.84	6.1	0.05	0.18
10	7 July 2018	14:23:48	140.59	35.17	6.0	0.19	0.09
11	8 January 2019	15:39:30	131.17	30.57	6.0	0.20	0.09
12	18 June 2019	16:22:19	139.48	38.61	6.7	0.84	1.00
13	19 April 2020	23:39:05	142.10	38.89	6.2	0.02	0.05
14	24 June 2020	23:39:05	141.11	35.55	6.1	0.02	0.07

.

Table 6. Japan. Probabilities of earthquake prediction with magnitudes $m\ge 6.0$ from seismological fields and from the same fields together with the GPS field.

	Feature Fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{10}$		Feature Fields ${\mathbf{S}}_9$, ${\mathbf{S}}_{10}$, and ${\mathbf{F}}_{10}$
Alarm Volume	Number of Successful Predictions	Probability of Successful Predictions	Number of Successful Predictions	Probability of Successful Predictions
0.01	1	0.07	2	0.14
0.05	5	0.36	3	0.21
0.10	5	0.36	6	0.43
0.15	5	0.36	9	0.64
0.20	10	0.71	10	0.71

shows the probabilities of a successful prediction for several alarm volumes. From Table 5 and Table 6, it can be seen that the forecast results for fields ${\mathbf{S}}_9$, ${\mathbf{S}}_{10}$, and ${\mathbf{F}}_{12}$ practically do not differ from the forecast for fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{10}$. Fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{11}$ proved to be the most informative for predicting earthquakes in California. The results of the forecasting and the probabilities of a successful prediction based on seismological feature fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{11}$, as well as based on fields ${\mathbf{S}}_9$, ${\mathbf{S}}_{11}$, ${\mathbf{F}}_4$, and ${\mathbf{F}}_6$ are shown in

Table 7. California. Results of forecasting earthquakes with magnitudes $m\ge 5.5$ from seismological fields and from the same fields together with GPS fields.

N	Date	Time	Long.	Lat.	Earthquake Magnitude	Alarm Volume for Prediction from Fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{11}$	Alarm Volume for Prediction from Fields ${\mathbf{S}}_9$, ${\mathbf{S}}_{11}$, ${\mathbf{F}}_4$, and ${\mathbf{F}}_6$
1	24 May 2013	03:47:08	121.06	40.19	5.7	0.03	0.01
2	24 August 2014	10:20:44	122.31	38.21	6.0	0.05	0.02
3	8 December 2016	08:18:00	118.90	38.38	5.6	0.15	0.05
4	28 December 2016	08:22:12	118.90	38.40	5.6	0.15	0.05
5	28 December 2016	09:13:47	118.90	38.38	5.5	0.15	0.05
6	23 June 2019	03:53:02	124.30	40.28	5.6	0.20	0.11
7	4 July 2019	17:33:49	117.50	35.70	6.4	1.00	0.05
8	6 July 2019	03:19:53	117.60	35.77	7.1	0.20	0.05
9	6 July 2019	03:47:53	117.75	35.90	5.5	0.20	0.05
10	15 May 2020	11:03:27	117.85	38.17	6.5	1.00	1.00
11	4 June 2020	01:32:11	117.43	35.61	5.5	0.11	0.04
12	24 June 2020	17:40:49	117.98	36.45	5.8	0.44	0.07

and

Table 8. California. Probabilities of earthquake prediction with magnitudes $m\ge 5.5$ from seismological fields and from the same fields together with GPS fields.

	Feature Fields ${\mathbf{S}}_9$ and ${\mathbf{S}}_{11}$		Feature Fields ${\mathbf{S}}_9$, ${\mathbf{S}}_{11}$, ${\mathbf{F}}_4$, and ${\mathbf{F}}_6$
Alarm Volume	Number of Successful Predictions	Probability of Successful Predictions	Number of Successful Predictions	Probability of Successful Predictions
0.01	0	0.00	1	0.08
0.05	1	0.08	6	0.50
0.10	2	0.17	10	0.83
0.15	3	0.25	11	0.92
0.20	9	0.75	11	0.92

.

Figure 16 shows the dependencies $U\left(V\right)$ obtained for Japan and for California from the GPS and seismological fields.

4. Discussion

The results of forecasting earthquakes in Japan and California based on the data of space geodesy were considered herein. The criterion for the effectiveness of forecasting according to GPS data was the probability of successful forecasting of strong earthquakes, obtained by the method of the minimum area of alarm. We can see from Figure 16 that for alarm volume $V=0.2$, the probability of forecasting target earthquakes according to GPS data was significantly higher than the probability of forecast based on random data and forecast based on the field of spatial density of earthquake epicenters. However, the $U\left(V\right)$ dependencies, the plots of which are shown in Figure 16, were found for a small number of target events. More convincing results about the effectiveness of GPS data in predicting strong earthquakes can be obtained using statistical hypothesis testing models. To do this, we compared the results of the prediction of target earthquakes using GPS data with the forecast results using a random field. If the target events were independent, we used the binomial distribution statistics.

Considering the results of earthquake forecasts in Japan. The test data contained 14 events. Table 1 shows that events 2, 3, and 4 were dependent. We will consider them as one event. Then, our sample consisted of $Q=12$ independent target events. Table 1 shows that $Q^{*}=6$ events were predicted in a sample of 12 events according to the GPS data. For the forecast, we chose the alarm volume $V=0.2$. We interpreted the value of the alarm volume V as the probability p of a successful prediction of events with a purely random distribution of alarm locations with a total alarm volume V in the analysis area. The observed frequency of a successful forecast $p^{*}=Q^{*}/Q=0.5$ exceeded the probability $p=0.2$ of successful forecast for a random field. Thus, the question arises: Can the observed difference between these two values be explained by the random nature of the sample? In terms of hypothesis testing, this means that the null hypothesis $p=0.2$ requires testing against the right-sided alternative $p^{*}>0.2$. Taking a binomial distribution model with $Q=12$ number of trials and a $p=0.2$ probability of a successful event, we found from the binomial distribution tables that the probability of obtaining $Q^{*}=6$ or more successful predictions was 0.0194. Since this value is significantly less than the generally accepted significance levels of 0.05, the null hypothesis was rejected.

A number of transformations were required to calculate feature fields based on GPS data. These included interpolation of time series at relatively small time intervals of discontinuity in the operation of GPS receiving stations, calculation of time series of components of the rates of horizontal displacements of stations, calculation of spatio-temporal fields of components of the rate of the Earth’s surface deformations, calculation of fields of invariants of velocities and fields of variation of invariants of rates in time, calculation of spatial correlation fields, and calculation of minimum and maximum values of correlations. A number of parameters were used in the algorithms for calculating these transformations: the time interval for estimating the daily displacement rates, the sizes of the spatial and temporal smoothing windows, and the windows for estimating the spatial correlation coefficients, as well as the time intervals for calculating the field of invariants of the strain rates. The GPS fields for Japan and California were calculated with the same parameters. The choice of transformation parameters, as well as the choice of the feature fields themselves, requires special study. In this work, such a study was not carried out. The types of transformations of the initial data into the fields of features and transformation parameters were selected based on qualitative considerations about the methods of cleaning signals from noise, recovering missing values, and disclosing information about the spatio-temporal properties of geodynamic processes.

The best forecast for Japan was given by the ${\mathbf{F}}_{10}$ field, and for California was given by the ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$ fields. Field ${\mathbf{F}}_{10}$ carries information about the consistency of the processes of changing the rates of various types of deformations in space. Fields ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$ carry information about the positive values of the change in the mean values of the divergence and shift of the strain rate with an interval of 361 days. When simulating earthquake predictions, we analyzed several types of such fields. A number of fields containing information on the change in strain rate invariants provided similar results for a successful prediction. We considered the field ${\mathbf{F}}_{10}$ for Japan and the fields ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$ for California to be more efficient, since, for them, the probability of a successful forecast practically did not change with small changes in the parameters of the alarm field: With the radius of the signal cylinder, R changed from 14 to 18 km when its generatrix T changed from 61 to 91 days.

Figure 17 shows maps of alarm zones when forecasting earthquakes using GPS data. Alarm zones highlight subdomains in which the values of the fields ${\mathbf{F}}_{10}$ for Japan, ${\mathbf{F}}_4$ and ${\mathbf{F}}_6$ for California take on anomalous values. Shaded alarm zones indicate the volume of the alarm. The level of shading of the alarm zones shows the value of the alarm volume. On alarm zone maps, these values are approximately equal to fractions of the alarm area relative to the analysis area. On the maps, target events fell into subdomains: for Japan, the proportion of the alarm area is about 15%, and for California, about 3%.

In the modeling, we used the experience of our previous studies to assess the relationship between data on horizontal displacements of the Earth’s surface with seismicity [19,52] and on the application of the method of the minimum area of alarm to earthquake prediction [50].

5. Conclusions

A number of seismically active regions are equipped with a rather dense network of GPS receiving stations that track the movements of the Earth’s surface. In our study, we tried to answer two questions: (1) Are space geodesy data effective for systematic earthquake prediction?; (2) is the earthquake forecast improved if seismological data are supplemented with space geodesy data? Obviously, the answers to these questions depend on the spatial density of the network of receiving stations, on the parameters of the time series of GPS measurements, on the method of preprocessing of GPS data, and on the method for forecasting earthquakes.

The results in this article were obtained from data from Japan and California. Space geodesy data are represented by daily time series of horizontal displacements of the Earth’s surface. The heart of the processing of GPS time series is occupied by the calculation of spatio-temporal fields of changes in the invariants of the seismic strain rate. To predict earthquakes, we used our developed machine learning method that we called the method of the minimum area of alarm. For Japan, it was shown that with an alarm volume of $V=0.2$, the probability of forecasting earthquakes with magnitudes $m\ge 6.0$ according to GPS data was statistically significantly higher than a forecast based on random data. It can be seen from the tables and Figure 16 that with the same volume of alarms, the probability of predicting earthquakes with magnitudes of $m\ge 5.5$ according to GPS data in California significantly exceeded the probability of forecasting based on random data and spatial density field earthquakes. Figure 16 also shows that adding GPS data to seismological data increased the probability of predicting an earthquake at an alarm volume of $V\le 0.2$.

The method of the minimum area of alarm is universal for various types of initial data since, for forecasting, all data are converted into uniform spatial and spatio-temporal grid fields. Modeling earthquake predictions in the regions of Japan and California showed that the anomalous field values that determine the change in the rates of deformations of the Earth’s surface, the change in the characteristics of the seismic regime, and the spatial correlation of these processes quite well distinguish the spatio-temporal areas preceding the appearance of the epicenters of strong earthquakes. Taking into account that this result was obtained for regions that differ significantly in seismotectonic and geodynamic regimes [53,54], it can be concluded that fields reflecting anomalous changes in seismotectonic and geodynamic processes are the most informative for predicting earthquakes.

The obtained results require additional verification, since the analysis material was limited by a relatively short observation time, two investigated regions, and a small number of target earthquakes.