Next Article in Journal
Distribution and Attribution of Terrestrial Snow Cover Phenology Changes over the Northern Hemisphere during 2001–2020
Next Article in Special Issue
Spatio-Temporal Distribution of Ground Deformation Due to 2018 Lombok Earthquake Series
Previous Article in Journal
The Intra-Tidal Characteristics of Tidal Front and Their Spring–Neap Tidal and Seasonal Variations in Bungo Channel, Japan
Previous Article in Special Issue
Analysis of Ocean Bottom Pressure Anomalies and Seismic Activities in the MedRidge Zone
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analyzing the Performance of GPS Data for Earthquake Prediction

The Institute for Information Transmission Problems, 127051 Moscow, Russia
*
Author to whom correspondence should be addressed.
Submission received: 11 March 2021 / Revised: 3 May 2021 / Accepted: 5 May 2021 / Published: 9 May 2021

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.

Graphical Abstract

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 Δ 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, F i , i = 1 , , 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 , , Q with magnitudes m 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 F i are set in a single coordinate grid. The values of the fields at the grid nodes n = 1 , , N correspond to the points of the I-dimensional feature space f ( n ) R I . During training, the spatio-temporal field Φ , called the alarm field, is calculated. The values of the alarm field ϕ ( n ) ϕ 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 ( t * ) , in which the appearance of the epicenter of the target earthquake is predicted at the interval [ t * , t * + Δ t ) .
Figure 1 shows operations in the interval [ t * Δ 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 Φ 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 * + Δ t ) if its epicenter falls into the alarm zone. The larger the product S ( t * ) Δ 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 f ( q ) R I . The precursor is associated with a feature space area whose points are, component-wise, greater than or equal to f ( q ) , i.e., the area w ( q ) = { f ( n ) R I : f i ( n ) f i ( q ) ,   i = 1 , 2 , , I } . We called this area the orthant w ( q ) with a vertex at f ( q ) and the points f ( n ) w ( q ) 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 f ( n ) has a circular base center at the grid node n with coordinates ( x ( n ) , y ( n ) , t ( n ) ) , a radius of the base R, and the elements [ ( x ( n ) , y ( n ) , t ( n ) ) , ( x ( n ) , y ( n ) , t ( n ) + T ) ] , where t ( n ) [ 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 ( q ) selects a set of grid nodes W ( q ) , | W ( q ) | = L ( q ) . Accordingly, it follows that an earthquake q with the epicenter coordinates ( x ( q ) , y ( q ) , t ( q ) ) can be predicted only if the cylinder with the center of the base of ( x ( q ) , y ( q ) , t ( q ) ) , a radius R, and the elements [ ( x ( q ) , y ( q ) , t ( q ) T ) , ( x ( q ) , y ( q ) , t ( q ) ) ] contains at least one base point. This cylinder is called a precursor cylinder. The precursor of event q is the point f ( q ) R I , which has the minimum value of the alarm volume v ( q ) = L ( q ) / 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 ( q ) (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 ϕ ( n ) ϕ 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 * = V L grid nodes. Usually, the quality of the forecast is determined by the dependence U ( V ) , 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 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 { f ( q ) , v ( q ) } . Arrange the precedents f ( q ) , q = 1 , , Q , in accordance with the increase in the alarm volume of the target events v ( 1 ) v ( 2 ) v ( q ) v ( Q ) .
  • Calculate the alarm field Φ .
    (a)
    Assign to the nodes of the grid of the alarm field Φ a value of 1.
    (b)
    Replace the value of 1 of the field Φ with V ( 1 ) = v ( 1 ) in the set W ( 1 ) of the grid nodes corresponding to the base points of the precedent f ( 1 ) ; replace the value of 1 with V ( 2 ) = | W ( 1 ) W ( 2 ) | / L at the grid nodes W ( 2 ) \ W ( 1 ) ; replace the value of 1 with V ( 3 ) = | W ( 1 ) W ( 2 ) W ( 3 ) | / L at the grid nodes W ( 3 ) \ ( W ( 1 ) W ( 2 ) ) and then sequentially replace the values of 1 with the corresponding values of the alarm volumes. The resulting field Φ takes the values V ( 1 ) V ( 2 ) V ( q ) V ( Q ) or 1. The value V ( 1 ) refers to grid nodes from the set W ( 1 ) , V ( 2 ) refers to grid nodes from the set ( W ( 2 ) \ W ( 1 ) ) , V ( 3 ) refers to grid nodes from the set ( W ( 3 ) \ ( W ( 1 ) W ( 2 ) ) ) , 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 Δ 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 f ( q ) , orthants w ( q ) , and alarm volumes of precursors v ( q ) 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 * + Δ 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 Φ V 0 and the alarm zone S ( t * ) 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 60 km for both regions, but with different magnitudes: m 6.0 in Japan and m 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 2 D 0.3–0.4. For this, in the time interval before the start of training, the 2D density field of earthquake epicenters 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 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 ( t ) and y ( t ) 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 ( t ) and g y ( t ) were determined by two coordinates of the GPS receiving station, spaced in time by the interval T 0 : g x ( t ) = ( x ( t ) x ( t T 0 ) ) / T 0 , g y ( t ) = ( y ( t ) 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) ×   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 V x and V y in the W–E and N–S directions. The fields were presented in the grid Δ x × Δ y × Δ t = 0 . 1 × 0 . 075 × 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 V x and 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 x n ( t ) = k = 1 K g x k ( t ) / r k p k = 1 K 1 / r k p ,
where V x n ( t ) 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 m a x , the values of which were used for interpolation, g x k ( t ) is the value of the W–E strain rate component for the station k , k = 1 , , K , at the time t, r k R m a x 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 m a x = 50 km, and p = 1 . If r k = 0 , then V x n ( t ) = g x k ( t ) . 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 N and 139 . 2 E for Japan: (A) coordinates of the receiving station, (B) daily velocities, and (C) values of the field components V x and 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 V x and 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 N and 123 . 2 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.
  • F 1 is the field of divergence of the strain rates:
    div V n = V x n x + V y n y
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.
  • F 2 is the field of rotor of the strain rates:
    rot V n = V x n y V y n x
The field values determine the direction and intensity of the field twisting around the vertical axis.
  • F 3 is the field of shear of the strain rates:
    sh V n = 1 2 V x n x V y n y 2 + V x n y + V y n x 2
Figure 8 and Figure 9 show a time series of fields F 1 , F 2 , and F 3 at the same points as in Figure 4, Figure 5, Figure 6 and Figure 7.
The fields of features F 4 , F 5 , and 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 Δ x × Δ y × Δ t = 0 . 1 × 0 . 075 × 30 days.
  • F 4 is the field of the temporal variations in the divergence strain rate.
The value of the field f 4 n ( t ) at time t is equal to the ratio of the difference ( div 2 n ¯ div 1 n ¯ ) between the mean values of the divergence in two consecutive intervals, namely, T 1 and T 2 , to the standard deviation of this difference σ n ( div ) , T 1 = T 2 = 361 days.
f 4 n ( t ) = ( div 2 n ¯ div 1 n ¯ ) / σ n ( div ) ,
where div 2 n ¯ is calculated from the values of field F 1 at the interval ( t T 2 , t ) , div 1 ¯ is calculated at the interval ( t T 2 T 1 , t T 2 ) .
  • F 5 is the field of the temporal variations in the rotor rate.
The values of field f 5 ( t ) are calculated similarly to the values of field F 4 ,
f 5 n ( t ) = ( rot 2 n ¯ rot 1 n ¯ ) / σ n ( rot ) .
  • F 6 is the field of the temporal variations in the shear deformation rate.
The values of field f 6 ( t ) are calculated similarly to the values of field F 4 ,
f 6 n ( t ) = ( sh 2 n ¯ sh 1 n ¯ ) / σ n ( sh ) ,
Figure 10 and Figure 11 show a time series of fields F 4 , F 5 , and F 6 at the same points as in Figure 4, Figure 5, Figure 6 and Figure 7.
Fields F 7 , F 8 , and F 9 represent spatial correlations of strain rate changes in a sliding window of 75 × 75 km 2 . With this window size, the correlation coefficients are estimated in approximately 70–80 grid points of the fields.
  • F 7 is the field of spatial correlations in fields F 4 and F 5 .
  • F 8 is the field of spatial correlations in fields F 4 and F 6 .
  • F 9 is the field of spatial correlations in fields F 5 and 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 F 7 , F 8 , and 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 F 7 , F 8 , and 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 F 10 field:
  • F 10 is the field of minimum values of the fields F 7 and F 9 :
    f 10 n ( t ) = min ( f 7 n ( t ) , f 9 n ( t ) )
The time series of field F 10 at point 35 . 45 N and 139 . 2 E in Japan and 39 . 49 N and 123 . 2 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 2.0 and a hypocenter depth of H 160 km.
  • S 1 is the 3D field of the density of earthquake epicenters.
  • S 2 is the 3D field of the mean earthquake magnitude.
  • S 3 is the 3D field of the negative temporal anomalies of the density of earthquake epicenters.
  • S 4 is the 3D field of the positive temporal anomalies of the density of earthquake epicenters.
  • S 5 is the 3D field of the negative temporal anomalies of the mean earthquake magnitude.
  • S 6 is the 3D field of the positive temporal anomalies of the mean earthquake magnitude.
  • 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.
  • 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 S 1 and S 2 was performed with the method of local kernel regression. The kernel function for the q-th earthquake has the form K q = [ c o s h 2 ( r q / R ) 2 c o s h 2 ( t q / T ) ] 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 S 1 and R = 100 km, and T = 730 days for S 2 . The field S 7 was calculated with the kernel function K q = [ c o s h 2 ( r q / R ) 2 ] 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 S 3 , S 4 , S 5 , and 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 S 3 , S 4 , S 5 , and 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.
  • S 9 = S 1 / ( S 8 + 0.001 ) is the field of ratios of the density values of the earthquake epicenters s 1 n 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 8 n + 0.001 ) .
  • 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.
  • 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 S 9 and 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 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 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 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 S 9 and 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 6.0 and 12 epicenters of target earthquakes in California with a magnitude of m 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 F 10 field and for California from the F 4 and F 6 fields. Table 1 and Table 2 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 0.20 . The final forecast results are summarized in Table 3 and Table 4. It can be seen that the forecast probabilities from GPS data with an alarm volume V 0.20 are significantly higher than the forecast probability for a random field and a spatial density field of earthquake epicenters 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 S 9 = S 1 / ( S 8 + 0.001 ) and field S 10 . The results of the forecasting based on seismological feature fields S 9 and S 10 , as well as based on fields S 9 , S 10 , and F 10 are shown in Table 5. Table 6 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 S 9 , S 10 , and F 12 practically do not differ from the forecast for fields S 9 and S 10 . Fields S 9 and 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 S 9 and S 11 , as well as based on fields S 9 , S 11 , F 4 , and F 6 are shown in Table 7 and Table 8.
Figure 16 shows the dependencies U ( V ) 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 ( V ) 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 F 10 field, and for California was given by the F 4 and F 6 fields. Field F 10 carries information about the consistency of the processes of changing the rates of various types of deformations in space. Fields F 4 and 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 F 10 for Japan and the fields F 4 and 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 F 10 for Japan, F 4 and 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 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 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 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.

Author Contributions

Conceptualization, V.G. and A.D.; methodology, V.G.; software, A.D. and K.P.; validation, V.G., A.D., and K.P.; formal analysis, V.G.; investigation, V.G. and K.P.; resources, A.D. and K.P.; data curation, A.D.; writing—original draft preparation, V.G.; writing—review and editing, A.D. and K.P.; visualization, A.D.; supervision, V.G.; project administration, V.G.; funding acquisition, V.G. All authors read and agreed to the published version of the manuscript.

Funding

This work was partially supported by the Russian Foundation for Basic Research, project 20-07-00445.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The GPS data were obtained from the Nevada Geodetic Laboratory (NGL), http://geodesy.unr.edu/about.php (accessed on 5 March 2021) at https://0-doi-org.brum.beds.ac.uk/10.1029/2018EO104623 [45]; the earthquake catalog data for Japan were obtained from the Japan Meteorological Agency (https://www.hinet.bosai.go.jp/) (accessed on 5 March 2021) [47,48]; the earthquake catalog data for California were obtained from the NEIC USGS catalog (https://earthquake.usgs.gov/earthquakes/feed/) (accessed on 5 March 2021) [55]. The earthquake catalog used in this study is produced by the Japan Meteorological Agency, in cooperation with the Ministry of Education, Culture, Sports, Science and Technology. The catalog is based on seismic data provided by the National Research Institute for Earth Science and Disaster Resilience, the Japan Meteorological Agency, Hokkaido University, Hirosaki University, Tohoku University, the University of Tokyo, Nagoya University, Kyoto University, Kochi University, Kyushu University, Kagoshima University, the National Institute of Advanced Industrial Science and Technology, the Geographical Survey Institute, Tokyo Metropolis, Shizuoka Prefecture, Hot Springs Research Institute of Kanagawa Prefecture, Yokohama City, and Japan Agency for Marine-Earth Science and Technology.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Sobolev, G.; Ponomarev, A. Earthquake Physics and Precursors; Publishing House Nauka: Moscow, Russia, 2003. [Google Scholar]
  2. Zavyalov, A. Medium-Term Earthquake Prediction: Principles, Methods, and Practice; Publishing House Nauka: Moscow, Russia, 2006. [Google Scholar]
  3. Lighthill, J. A Critical Review of VAN: Earthquake Prediction from Seismic Electrical Signals; World Scientific: Singapore, 1996. [Google Scholar]
  4. King, C.Y. Gas geochemistry applied to earthquake prediction: An overview. J. Geophys. Res. Solid Earth 1986, 91, 12269–12281. [Google Scholar] [CrossRef]
  5. Kosobokov, V.; Kejlis-Borok, V.I. Prognoz Zemletrjasenij: Osnovy, Realizacija, Perspektivy; GEOS: Moscow, Russia, 2005. [Google Scholar]
  6. Gitis, V.G.; Derendyaev, A.B. Web-Based GIS platform for automatic prediction of earthquakes. In Proceedings of the International Conference on Computational Science and Its Applications, Melbourne, Australia, 2–5 July 2018; Springer: Berlin/Heidelberg, Germany, 2018; pp. 268–283. [Google Scholar]
  7. Murai, S.; Araki, H. Earthquake Prediction Using GPS—A New Method Based on GPS Network Triangles. GIM Int. 2003, 17, 34–37. [Google Scholar]
  8. Murai, S.; Araki, H. Prediction of earthquake and volcanic eruption using GPS. Asian J. Geoinform. 2004, 4, 85–90. [Google Scholar]
  9. Murai, S. Can we predict earthquakes with GPS data? Int. J. Digit. Earth 2010, 3, 83–90. [Google Scholar] [CrossRef]
  10. Borghi, A.; Aoudia, A.; Riva, R.E.; Barzaghi, R. GPS monitoring and earthquake prediction: A success story towards a useful integration. Tectonophysics 2009, 465, 177–189. [Google Scholar] [CrossRef]
  11. Klein, E.; Vigny, C.; Fleitout, L.; Grandin, R.; Jolivet, R.; Rivera, E.; Métois, M. A comprehensive analysis of the Illapel 2015 Mw8. 3 earthquake from GPS and InSAR data. Earth Planet. Sci. Lett. 2017, 469, 123–134. [Google Scholar] [CrossRef] [Green Version]
  12. Liu, Y.; Ye, S.; Jiang, P.; Song, W.; Lou, Y. Combining GPS + GLONASS observations to improve the fixing percentage and precision of long baselines with limited data. Adv. Space Res. 2016, 57, 1258–1267. [Google Scholar] [CrossRef]
  13. Wang, Q.; Guo, Y.; Yu, L.; Li, P. Earthquake prediction based on spatio-temporal data mining: An LSTM network approach. IEEE Trans. Emerg. Top. Comput. 2017, 8, 148–158. [Google Scholar] [CrossRef]
  14. Chen, C.H.; Yeh, T.K.; Wen, S.; Meng, G.; Han, P.; Tang, C.C.; Liu, J.Y.; Wang, C.H. Unique pre-earthquake deformation patterns in the spatial domains from GPS in Taiwan. Remote Sens. 2020, 12, 366. [Google Scholar] [CrossRef] [Green Version]
  15. Li, N.; Kong, X.; Lin, L. Anomalies in continuous GPS data as precursors of 15 large earthquakes in Western North America during 2007–2016. Earth Sci. Inform. 2020, 13, 163–174. [Google Scholar] [CrossRef]
  16. Murray-Moraleda, J. GPS: Applications in crustal deformation monitoring. In Extreme Environmental Events; Meyers, R.A., Ed.; Springer: New York, NY, USA, 2011; pp. 589–622. [Google Scholar]
  17. Yuan, L.; Chao, B.F.; Ding, X.; Zhong, P. The tidal displacement field at Earth’s surface determined using global GPS observations. J. Geophys. Res. Solid Earth 2013, 118, 2618–2632. [Google Scholar] [CrossRef] [Green Version]
  18. Kuzikov, S.; Mukhamediev, S.A. Structure of the present-day velocity field of the crust in the area of the Central-Asian GPS network. Izv. Phys. Solid Earth 2010, 46, 584–601. [Google Scholar] [CrossRef]
  19. Sobolev, G.; Zakrzhevskaya, N.; Akatova, K.; Gitis, V.; Derendyaev, A.; Bragin, V.; Sycheva, N.; Kuzikov, S. Dynamics of interaction between fields of seismicity and surface deformations (Bishkek geodynamic test area). Izv. Phys. Solid Earth 2010, 46, 817–838. [Google Scholar] [CrossRef]
  20. Keilis-Borok, V.; Soloviev, A.A. Nonlinear Dynamics of the Lithosphere and Earthquake Prediction; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013. [Google Scholar]
  21. Kossobokov, V. User manual for M8. In Algorithms for Earthquake Statistics and Prediction; Healy, J.H., Keilis-Borok, V.I., Lee, W.H.K., Eds.; Springer: Berlin/Heidelberg, Germany, 1997; Volume 6, pp. 167–222. [Google Scholar]
  22. Kossobokov, V.; Shebalin, P. Earthquake prediction. In Nonlinear Dynamics of the Lithosphere and Earthquake Prediction; Springer: Berlin/Heidelberg, Germany, 2003; pp. 141–207. [Google Scholar]
  23. Corbi, F.; Sandri, L.; Bedford, J.; Funiciello, F.; Brizzi, S.; Rosenau, M.; Lallemand, S. Machine learning can predict the timing and size of analog earthquakes. Geophys. Res. Lett. 2019, 46, 1303–1311. [Google Scholar] [CrossRef] [Green Version]
  24. Zavyalov, A. Intermediate Term Earthquake Prediction; Nauka: Moscow, Russia, 2006; Volume 52, pp. 641–646. [Google Scholar]
  25. Shebalin, P.N.; Narteau, C.; Zechar, J.D.; Holschneider, M. Combining earthquake forecasts using differential probability gains. Earth Planets Space 2014, 66, 37. [Google Scholar] [CrossRef] [Green Version]
  26. Amei, A.; Fu, W.; Ho, C.H. Time series analysis for predicting the occurrences of large scale earthquakes. Int. J. Appl. Sci. Technol. 2012, 2, 64–75. [Google Scholar]
  27. Marzocchi, W.; Zechar, J.D. Earthquake forecasting and earthquake prediction: Different approaches for obtaining the best model. Seismol. Res. Lett. 2011, 82, 442–448. [Google Scholar] [CrossRef]
  28. Rhoades, D.A. Application of the EEPAS model to forecasting earthquakes of moderate magnitude in southern California. Seismol. Res. Lett. 2007, 78, 110–115. [Google Scholar] [CrossRef] [Green Version]
  29. Rhoades, D.A. Mixture models for improved earthquake forecasting with short-to-medium time horizons. Bull. Seismol. Soc. Am. 2013, 103, 2203–2215. [Google Scholar] [CrossRef]
  30. Alves, E. Earthquake forecasting using neural networks: Results and future work. Nonlinear Dyn. 2006, 44, 341–349. [Google Scholar] [CrossRef]
  31. Priambodo, B.; Mahmudy, W.; Rahman, M. Earthquake Magnitude and Grid-Based Location Prediction using Backpropagation Neural Network. Knowl. Eng. Data Sci. 2020, 3, 28–39. [Google Scholar] [CrossRef]
  32. Asim, K.; Idris, A.; Iqbal, T.; Martínez-Álvarez, F. Earthquake prediction model using support vector regressor and hybrid neural networks. PLoS ONE 2018, 13, e0199004. [Google Scholar] [CrossRef] [PubMed]
  33. Panakkat, A.; Adeli, H. Neural network models for earthquake magnitude prediction using multiple seismicity indicators. Int. J. Neural Syst. 2007, 17, 13–33. [Google Scholar] [CrossRef] [PubMed]
  34. Galkina, A.; Grafeeva, N. Machine learning methods for earthquake prediction: A survey. In Proceedings of the Fourth Conference on Software Engineering and Information Management (SEIM-2019), Saint Petersburg, Russia, 13 April 2019. [Google Scholar]
  35. Mignan, A.; Broccardo, M. Neural Network Applications in Earthquake Prediction (1994–2019): Meta-Analytic and Statistical Insights on their Limitations. Seismol. Res. Lett. 2020, 91, 2330–2342. [Google Scholar]
  36. Moustra, M.; Avraamides, M.; Christodoulou, C. Artificial neural networks for earthquake prediction using time series magnitude data or seismic electric signals. Expert Syst. Appl. 2011, 38, 15032–15039. [Google Scholar] [CrossRef]
  37. Gitis, V.G.; Derendyaev, A.B. Machine Learning Methods for Seismic Hazards Forecast. Geosciences 2019, 9, 308. [Google Scholar] [CrossRef] [Green Version]
  38. Bishop, C.M. Machine learning and pattern recognition. In Information Science and Statistics; Springer: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
  39. Kotsiantis, S.B.; Zaharakis, I.; Pintelas, P. Supervised machine learning: A review of classification techniques. Emerg. Artif. Intell. Appl. Comput. Eng. 2007, 160, 3–24. [Google Scholar]
  40. Khan, S.S.; Madden, M.G. A survey of recent trends in one class classification. In Proceedings of the Irish Conference on Artificial Intelligence and Cognitive Science, Dublin, Ireland, 19–21 August 2009; Springer: Berlin/Heidelberg, Germany, 2009; pp. 188–197. [Google Scholar]
  41. Gitis, V.; Derendyaev, A.; Petrov, K. A Method for Identification of Anomalous Geological Zones. J. Commun. Technol. Electron. 2020, 65, 1531–1541. [Google Scholar] [CrossRef]
  42. Metz, C.E. Basic principles of ROC analysis. In Seminars in Nuclear Medicine; Elsevier: Amsterdam, The Netherlands, 1978; Volume 8, pp. 283–298. [Google Scholar]
  43. Kim, D.S.; McCabe, C.J.; Yamasaki, B.L.; Louie, K.A.; King, K.M. Detecting random responders with infrequency scales using an error-balancing threshold. Behav. Res. Methods 2018, 50, 1960–1970. [Google Scholar] [CrossRef]
  44. Bradley, A.P. The use of the area under the ROC curve in the evaluation of machine learning algorithms. Pattern Recognit. 1997, 30, 1145–1159. [Google Scholar] [CrossRef] [Green Version]
  45. Blewitt, G.; Hammond, W.C.; Kreemer, C. Harnessing the GPS data explosion for interdisciplinary science. Eos 2018, 99. [Google Scholar] [CrossRef]
  46. Zadeh, L.A. Fuzzy logic. Computer 1988, 21, 83–93. [Google Scholar] [CrossRef]
  47. Okada, Y.; Kasahara, K.; Hori, S.; Obara, K.; Sekiguchi, S.; Fujiwara, H.; Yamamoto, A. Recent progress of seismic observation networks in Japan—Hi-net, F-net, K-NET and KiK-net—. Earth Planets Space 2004, 56, 15–28. [Google Scholar] [CrossRef] [Green Version]
  48. Obara, K.; Kasahara, K.; Hori, S.; Okada, Y. A densely distributed high-sensitivity seismograph network in Japan: Hi-net by national research institute for earth science and disasterprevention. Rev. Sci. Instrum. 2005, 76, 021301. [Google Scholar] [CrossRef]
  49. Barnhart, W.D.; Hayes, G.P.; Wald, D.J. Global earthquake response with imaging geodesy: Recent examples from the USGS NEIC. Remote Sens. 2019, 11, 1357. [Google Scholar] [CrossRef] [Green Version]
  50. Gitis, V.; Derendyaev, A. The Method of the Minimum Area of Alarm for Earthquake Magnitude Prediction. Front. Earth Sci. 2020, 8, 482. [Google Scholar] [CrossRef]
  51. Mjachkin, V.; Brace, W.F.; Sobolev, G.A.; Dieterich, J.H. Two models for earthquake forerunners. In Earthquake Prediction and Rock Mechanics; Springer: Berlin/Heidelberg, Germany, 1975; pp. 169–181. [Google Scholar]
  52. Gitis, V.; Derendyaev, A. Spatio-temporal analysis of Earth’s surface deformation by GPS and InSAR Data. In Proceedings of the International Conference on Computational Science and Its Applications, Guimaraes, Portugal, 30 June–3 July 2014; Springer: Berlin/Heidelberg, Germany, 2014; pp. 237–251. [Google Scholar]
  53. Garagash, I.; Bondur, V.; Gokhberg, M.; Steblov, G. Three-Year Experience of the Fortnight Forecast of Seismicity in Southern California on the Basis of Geomechanical Model and the Seismic Data. In Proceedings of the AGU Fall Meeting, San Francisco, CA, USA, 5–9 December 2011; Volume 2011. [Google Scholar]
  54. Lobkovsky, L.; Vladimirova, I.; Gabsatarov, Y.V.; Steblov, G. Seismotectonic Deformations Related to the 2011 Tohoku Earthquake at Different Stages of the Seismic Cycle, Based on Satellite Geodetic Observations. In Doklady Earth Sciences; Springer: Berlin/Heidelberg, Germany, 2018; Volume 481, pp. 1060–1065. [Google Scholar]
  55. Masse, R.; Needham, R. Earthquakes Volcanoes (USGS); NEIC: The National Earthquake Information Center: Golden, CO, USA, 1989; Volume 21, pp. 4–45.
Figure 1. Flowchart of operations in the intervals Δ t before the forecast.
Figure 1. Flowchart of operations in the intervals Δ t before the forecast.
Remotesensing 13 01842 g001
Figure 2. Flowchart of the algorithm.
Figure 2. Flowchart of the algorithm.
Remotesensing 13 01842 g002
Figure 3. Areas of analysis, Global Positioning System (GPS) ground receiving stations, and the epicenters of target earthquakes. (Left) Japan, epicenters of target earthquakes with magnitude m 6.0 in the interval 1 January 2011–26 July 2020; (Right) California, epicenters of target earthquakes with magnitude m 5.5 , in the interval 23 December 2009–14 November 2020. The epicenters of the target earthquakes for which the forecast was tested are highlighted in red and the epicenters on which the forecast was initial trained in yellow.
Figure 3. Areas of analysis, Global Positioning System (GPS) ground receiving stations, and the epicenters of target earthquakes. (Left) Japan, epicenters of target earthquakes with magnitude m 6.0 in the interval 1 January 2011–26 July 2020; (Right) California, epicenters of target earthquakes with magnitude m 5.5 , in the interval 23 December 2009–14 November 2020. The epicenters of the target earthquakes for which the forecast was tested are highlighted in red and the epicenters on which the forecast was initial trained in yellow.
Remotesensing 13 01842 g003
Figure 4. Time series of displacements in the W–E direction of the Japanese receiving station at 35 . 45 N and 139 . 2 E. (Left) Interval 1 January 2009–26 July 2020. (Right) Interval 1 January 2016–14 December 2016. (A,D) time series of coordinates x ( t ) of the station (mm); (B,E) time series of daily rates g x k ( t ) of the station (mm/day); (C,F) time series V x n ( t ) at the point n of the field V x (mm/day).
Figure 4. Time series of displacements in the W–E direction of the Japanese receiving station at 35 . 45 N and 139 . 2 E. (Left) Interval 1 January 2009–26 July 2020. (Right) Interval 1 January 2016–14 December 2016. (A,D) time series of coordinates x ( t ) of the station (mm); (B,E) time series of daily rates g x k ( t ) of the station (mm/day); (C,F) time series V x n ( t ) at the point n of the field V x (mm/day).
Remotesensing 13 01842 g004
Figure 5. Time series of displacements in the N–S direction of the Japanese receiving station at 35 . 45 N and 139 . 2 E. (Left) Interval 1 January 2009–26 July 2020. (Right) Interval 1 January 2016–14 December 2016. (A,D) time series of coordinates y ( t ) of the station (mm); (B,E) time series of daily rates g y k ( t ) of the station (mm/day); (C,F) time series V y n ( t ) at the point n of the field V y (mm/day).
Figure 5. Time series of displacements in the N–S direction of the Japanese receiving station at 35 . 45 N and 139 . 2 E. (Left) Interval 1 January 2009–26 July 2020. (Right) Interval 1 January 2016–14 December 2016. (A,D) time series of coordinates y ( t ) of the station (mm); (B,E) time series of daily rates g y k ( t ) of the station (mm/day); (C,F) time series V y n ( t ) at the point n of the field V y (mm/day).
Remotesensing 13 01842 g005
Figure 6. Time series of displacements in the W–E direction of the Californian receiving station at 39 . 49 N and 123 . 2 W. (Left) Interval 1 January 2008–14 November 2020. (Right) Interval 1 January 2013–1 October 2013. (A,D) time series of coordinates x ( t ) of the station (mm); (B,E) time series of daily rates g x k ( t ) of the station (mm/day); (C,F) time series V x n ( t ) at the point n of the field V x (mm/day).
Figure 6. Time series of displacements in the W–E direction of the Californian receiving station at 39 . 49 N and 123 . 2 W. (Left) Interval 1 January 2008–14 November 2020. (Right) Interval 1 January 2013–1 October 2013. (A,D) time series of coordinates x ( t ) of the station (mm); (B,E) time series of daily rates g x k ( t ) of the station (mm/day); (C,F) time series V x n ( t ) at the point n of the field V x (mm/day).
Remotesensing 13 01842 g006
Figure 7. Time series of displacements in the N–S direction of the Californian receiving station at 39 . 49 N and 123 . 2 W. (Left) Interval 1 January 2008–14 November 2020. (Right) Interval 1 January 2013–1 October 2013. (A,D) time series of coordinates y ( t ) of the station (mm); (B,E) time series of daily rates g y k ( t ) of the station (mm/day); (C,F) time series V y n ( t ) at the point n of the field V y (mm/day).
Figure 7. Time series of displacements in the N–S direction of the Californian receiving station at 39 . 49 N and 123 . 2 W. (Left) Interval 1 January 2008–14 November 2020. (Right) Interval 1 January 2013–1 October 2013. (A,D) time series of coordinates y ( t ) of the station (mm); (B,E) time series of daily rates g y k ( t ) of the station (mm/day); (C,F) time series V y n ( t ) at the point n of the field V y (mm/day).
Remotesensing 13 01842 g007
Figure 8. Japan. Time series at point 35 . 45 N and 139 . 2 E. (A) The divergence field F 1 ( 10 6 /day); (B) the rotor field F 2 ( 10 6 /day); (C) the shear field F 3 ( 10 6 /day).
Figure 8. Japan. Time series at point 35 . 45 N and 139 . 2 E. (A) The divergence field F 1 ( 10 6 /day); (B) the rotor field F 2 ( 10 6 /day); (C) the shear field F 3 ( 10 6 /day).
Remotesensing 13 01842 g008
Figure 9. California. Time series at point 39 . 49 N and 123 . 2 W. (A) The divergence field F 1 ( 10 6 /day); (B) the rotor field F 2 ( 10 6 /day); (C) the shear field F 3 ( 10 6 /day).
Figure 9. California. Time series at point 39 . 49 N and 123 . 2 W. (A) The divergence field F 1 ( 10 6 /day); (B) the rotor field F 2 ( 10 6 /day); (C) the shear field F 3 ( 10 6 /day).
Remotesensing 13 01842 g009
Figure 10. Japan. Time series at point 35 . 45 N and 139 . 2 E. (A) Field F 4 ; (B) field F 5 ; (C) field F 6 .
Figure 10. Japan. Time series at point 35 . 45 N and 139 . 2 E. (A) Field F 4 ; (B) field F 5 ; (C) field F 6 .
Remotesensing 13 01842 g010
Figure 11. California. Time series at point 39 . 49 N and 123 . 2 W. (A) Field F 4 ; (B) field F 5 ; (C) field F 6 .
Figure 11. California. Time series at point 39 . 49 N and 123 . 2 W. (A) Field F 4 ; (B) field F 5 ; (C) field F 6 .
Remotesensing 13 01842 g011
Figure 12. Japan. Time series at point 35 . 45 N and 139 . 2 E. (A) Field F 7 ; (B) field F 8 ; (C) field F 9 .
Figure 12. Japan. Time series at point 35 . 45 N and 139 . 2 E. (A) Field F 7 ; (B) field F 8 ; (C) field F 9 .
Remotesensing 13 01842 g012
Figure 13. California. Time series at point 39 . 49 N and 123 . 2 W. (A) Field F 7 ; (B) field F 8 ; (C) field F 9 .
Figure 13. California. Time series at point 39 . 49 N and 123 . 2 W. (A) Field F 7 ; (B) field F 8 ; (C) field F 9 .
Remotesensing 13 01842 g013
Figure 14. Japan. The time series of field F 10 at point 35 . 45 N and 139 . 2 .
Figure 14. Japan. The time series of field F 10 at point 35 . 45 N and 139 . 2 .
Remotesensing 13 01842 g014
Figure 15. California. The time series of field F 10 at point 39 . 49 N and 123 . 2 W.
Figure 15. California. The time series of field F 10 at point 39 . 49 N and 123 . 2 W.
Remotesensing 13 01842 g015
Figure 16. Dependences U ( V ) of the probability of a successful earthquake prediction U on the alarm volume V obtained with the different fields. (Left) Japan: (1) field S 7 ; (2) field F 10 ; (3) fields S 9 and S 10 ; (4) fields S 9 , S 10 , and F 10 . (Right) California: (1) field S 7 ; (2) fields F 4 and F 6 ; (3) fields S 9 and S 11 ; (4) fields S 9 , S 11 , F 4 , and F 6 .
Figure 16. Dependences U ( V ) of the probability of a successful earthquake prediction U on the alarm volume V obtained with the different fields. (Left) Japan: (1) field S 7 ; (2) field F 10 ; (3) fields S 9 and S 10 ; (4) fields S 9 , S 10 , and F 10 . (Right) California: (1) field S 7 ; (2) fields F 4 and F 6 ; (3) fields S 9 and S 11 ; (4) fields S 9 , S 11 , F 4 , and F 6 .
Remotesensing 13 01842 g016
Figure 17. Alarm area zones of prognosis for Japan and California. (Left) Japan alarm area zone on 14 December 2018 with the epicenter of target earthquake with magnitude m = 6.3 on 8 January 2019; (Right) California alarm area zone on 17 April 2020 with the epicenter of target earthquake with magnitude m = 6.5 on 15 May 2020. The epicenters of the target earthquake are highlighted in red. Alarm volume values are shown on palettes in the picture. Field values with V > 0.2 are not shown.
Figure 17. Alarm area zones of prognosis for Japan and California. (Left) Japan alarm area zone on 14 December 2018 with the epicenter of target earthquake with magnitude m = 6.3 on 8 January 2019; (Right) California alarm area zone on 17 April 2020 with the epicenter of target earthquake with magnitude m = 6.5 on 15 May 2020. The epicenters of the target earthquake are highlighted in red. Alarm volume values are shown on palettes in the picture. Field values with V > 0.2 are not shown.
Remotesensing 13 01842 g017
Table 1. Japan. Results of forecasting earthquakes with magnitudes m 6.0 from field F 10 and field S 7 .
Table 1. Japan. Results of forecasting earthquakes with magnitudes m 6.0 from field F 10 and field S 7 .
NDateTimeLong.Lat.Earthquake MagnitudeAlarm Volume for Prediction from Field F 10 Alarm Volume for Prediction from Field S 7
114 January 201606:25:33142.8041.976.70.460.24
214 April 201615:26:34130.8132.746.50.130.50
314 April 201618:03:46130.7832.706.40.130.50
415 April 201619:25:05130.7632.767.30.130.50
521 October 201608:07:22133.8635.386.60.321.00
624 November 201600:23:36141.3537.186.20.690.02
728 December 201615:38:49140.5736.726.30.040.02
88 April 201819:32:30132.5935.196.10.200.86
918 June 201801:58:34135.6234.846.10.450.76
107 July 201814:23:48140.5935.176.00.150.16
118 January 201915:39:30131.1730.576.00.100.28
1218 June 201916:22:19139.4838.616.70.681.00
1319 April 202023:39:05142.1038.896.20.200.03
1424 June 202023:39:05141.1135.556.10.770.04
Table 2. California. Results of forecasting earthquakes with magnitudes m 5.5 from fields F 4 and F 6 and field S 7 .
Table 2. California. Results of forecasting earthquakes with magnitudes m 5.5 from fields F 4 and F 6 and field S 7 .
NDateTimeLong.Lat.Earthquake MagnitudeAlarm Volume for Prediction from Fields F 4 and F 6 Alarm Volume for Prediction from Field S 7
124 May 201303:47:08121.0640.195.70.180.84
224 August 201410:20:44122.3138.216.00.170.62
38 December 201608:18:00118.9038.385.60.100.48
428 December 201608:22:12118.9038.405.60.100.48
528 December 201609:13:47118.9038.385.50.100.48
623 June 201903:53:02124.3040.285.60.470.30
74 July 201917:33:49117.5035.706.40.100.04
86 July 201903:19:53117.6035.777.10.100.04
96 July 201903:47:53117.7535.905.50.100.04
1015 May 202011:03:27117.8538.176.50.030.29
114 June 202001:32:11117.4335.615.50.480.04
1224 June 202017:40:49117.9836.455.80.030.29
Table 3. Japan. Probabilities of predicting earthquakes with magnitudes m 6.0 .
Table 3. Japan. Probabilities of predicting earthquakes with magnitudes m 6.0 .
Feature Field F 10 Spatial Field of Epicenter Density S 7
Alarm VolumeNumber of Successful PredictionsProbability of Successful PredictionsNumber of Successful PredictionsProbability of Successful Predictions
0.0100.0000.00
0.0510.0740.29
0.1020.1440.29
0.1550.3640.29
0.2080.5750.36
Table 4. California. Probabilities of predicting earthquakes with magnitudes m 5.5 .
Table 4. California. Probabilities of predicting earthquakes with magnitudes m 5.5 .
Feature Fields F 4 and F 6 Spatial Field of Epicenter Density S 7
Alarm VolumeNumber of Successful PredictionsProbability of Successful PredictionsNumber of Successful PredictionsProbability of Successful Predictions
0.0100.0000.00
0.0520.1740.33
0.1020.1740.33
0.1580.6740.33
0.20100.8340.33
Table 5. Japan. Results of forecasting earthquakes with magnitudes m 6.0 from seismological fields and from the same fields together with the GPS field.
Table 5. Japan. Results of forecasting earthquakes with magnitudes m 6.0 from seismological fields and from the same fields together with the GPS field.
NDateTimeLong.Lat.Earthquake MagnitudeAlarm Volume for Prediction from Fields S 9 and S 10 Alarm Volume for Prediction from Fields S 9 , S 10 , and F 10
114 January 201606:25:33142.8041.976.70.200.24
214 April 201615:26:34130.8132.746.50.170.15
314 April 201618:03:46130.7832.706.40.170.15
415 April 201619:25:05130.7632.767.30.170.15
521 October 201608:07:22133.8635.386.61.001.00
624 November 201600:23:36141.3537.186.20.010.01
728 December 201615:38:49140.5736.726.30.040.01
88 April 201819:32:30132.5935.196.10.840.38
918 June 201801:58:34135.6234.846.10.050.18
107 July 201814:23:48140.5935.176.00.190.09
118 January 201915:39:30131.1730.576.00.200.09
1218 June 201916:22:19139.4838.616.70.841.00
1319 April 202023:39:05142.1038.896.20.020.05
1424 June 202023:39:05141.1135.556.10.020.07
Table 6. Japan. Probabilities of earthquake prediction with magnitudes m 6.0 from seismological fields and from the same fields together with the GPS field.
Table 6. Japan. Probabilities of earthquake prediction with magnitudes m 6.0 from seismological fields and from the same fields together with the GPS field.
Feature Fields S 9 and S 10 Feature Fields S 9 , S 10 , and F 10
Alarm VolumeNumber of Successful PredictionsProbability of Successful PredictionsNumber of Successful PredictionsProbability of Successful Predictions
0.0110.0720.14
0.0550.3630.21
0.1050.3660.43
0.1550.3690.64
0.20100.71100.71
Table 7. California. Results of forecasting earthquakes with magnitudes m 5.5 from seismological fields and from the same fields together with GPS fields.
Table 7. California. Results of forecasting earthquakes with magnitudes m 5.5 from seismological fields and from the same fields together with GPS fields.
NDateTimeLong.Lat.Earthquake MagnitudeAlarm Volume for Prediction from Fields S 9 and S 11 Alarm Volume for Prediction from Fields S 9 , S 11 , F 4 , and F 6
124 May 201303:47:08121.0640.195.70.030.01
224 August 201410:20:44122.3138.216.00.050.02
38 December 201608:18:00118.9038.385.60.150.05
428 December 201608:22:12118.9038.405.60.150.05
528 December 201609:13:47118.9038.385.50.150.05
623 June 201903:53:02124.3040.285.60.200.11
74 July 201917:33:49117.5035.706.41.000.05
86 July 201903:19:53117.6035.777.10.200.05
96 July 201903:47:53117.7535.905.50.200.05
1015 May 202011:03:27117.8538.176.51.001.00
114 June 202001:32:11117.4335.615.50.110.04
1224 June 202017:40:49117.9836.455.80.440.07
Table 8. California. Probabilities of earthquake prediction with magnitudes m 5.5 from seismological fields and from the same fields together with GPS fields.
Table 8. California. Probabilities of earthquake prediction with magnitudes m 5.5 from seismological fields and from the same fields together with GPS fields.
Feature Fields S 9 and S 11 Feature Fields S 9 , S 11 , F 4 , and F 6
Alarm VolumeNumber of Successful PredictionsProbability of Successful PredictionsNumber of Successful PredictionsProbability of Successful Predictions
0.0100.0010.08
0.0510.0860.50
0.1020.17100.83
0.1530.25110.92
0.2090.75110.92
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Gitis, V.; Derendyaev, A.; Petrov, K. Analyzing the Performance of GPS Data for Earthquake Prediction. Remote Sens. 2021, 13, 1842. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13091842

AMA Style

Gitis V, Derendyaev A, Petrov K. Analyzing the Performance of GPS Data for Earthquake Prediction. Remote Sensing. 2021; 13(9):1842. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13091842

Chicago/Turabian Style

Gitis, Valeri, Alexander Derendyaev, and Konstantin Petrov. 2021. "Analyzing the Performance of GPS Data for Earthquake Prediction" Remote Sensing 13, no. 9: 1842. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13091842

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop