Deep Learning of Detecting Ionospheric Precursors Associated With M ≥ 6.0 Earthquakes in Taiwan

Abstract

A short‐term (30 days before an earthquake) prediction of an earthquake is a big challenge in seismology. As a first step, we apply deep learning to the ionospheric total electron content (TEC) data between 2003 and 2014 to detect the seismo‐ionospheric precursors of M ≥ 6.0 earthquakes in Taiwan. The bidirectional Long Short‐Term Memory (Bi‐LSTM) network is employed to use observed input data (features) to obtain the sequential TEC variations. The five input features are sequential vectors of TEC, the geomagnetic index Dst, the solar activity index F10.7, sunspot number (SSN), and solar emission index Lyman‐α. The daily values of F10.7, SSN, and Lyman‐α are converted into hourly values, depending on the solar elevation angle. The calculated hourly TEC variations can be more precisely predicted with this data conversion. We calculate the normalized difference of errors between two 15‐day adjacent stages as the “relative error”. Three trained models with the best discrimination between the relative errors of earthquake and no‐earthquake cases are chosen as classifiers. These three classifiers are then used to have a majority vote to declare whether the 30‐day period is related to the preparation of an earthquake or not. The results show that all 22 positive cases (earthquakes) are successfully predicted, giving a true positive rate of 100%. Among the 19 negative cases (normal cases), 10 of them are true negative.” Overall, a high accuracy of 78.05% is obtained.

Introduction

The ionosphere is a region of the upper atmosphere above 60 km altitude. This region is partially ionized by Extreme Ultra-Violet (EUV) and X-rays of the solar radiation, and by collisions with energetic particles that penetrate the atmosphere. The electron density variation with altitude displays a similar basic structure at all latitudes. The ionospheric total electron content (TEC) is defined as the total number of electrons per unit ground area. The Vertical TEC (VTEC) is calculated along a vertical path from the Global Position System (GPS) satellite to a ground receiver or a Global Navigation Satellite System (GNSS) receiver. The Slant TEC is estimated along the line from the GPS satellite to a ground receiver or GNSS receiver. TEC is measured in units of 10¹⁶ electrons per square meter, where 10¹⁶ electrons/m² = 1 TEC unit (TECu) (Hofmann-Wellenhof et al., 2001). The TEC is modified by changing solar EUV radiation, geomagnetic storms, and the atmospheric waves that propagate up from the lower atmosphere. The TEC will therefore depend on local time, latitude, longitude, season, geomagnetic conditions, solar cycle and activity, and troposphere conditions.

Leonard and Barnes (1965) reported the ionospheric disturbances after the 1964 Alaska earthquake by four ionosonde stations. The seismo-ionospheric anomalies have been observed before and after large earthquakes in the following decades. Molchanov et al. (1998) discovered a significant change of the VLF signal which passed through the epicenter of the 1995 Kobe earthquake. They suggested that the signal change was caused by the decrease of the VLF reflection height of the ionosphere. Pulinets (1998) found that seismo-ionospheric variations have their individual features, different from variations induced by magnetic disturbances and other kinds of external activity. Liu et al. (2000) analyzed the critical frequency foF2 by the Chung-Li ionosonde before M ≥ 6.0 earthquakes during 1994–1999. They found that the seismo-ionospheric signatures appeared 1–6 days prior to the earthquakes. Liu et al. (2001) further examined the GPS-TEC before the 1999 Chi-Chi earthquake. The results showed that TEC decreases significantly 1, 3, 4 days before the earthquake and around the epicenter. Pulinets and Boyarchuk (2005) described seismo-ionospheric coupling, which is associated with the electric field penetration into the ionosphere. Liu et al. (2006) reported the pre-earthquake ionospheric anomalies (PEIAs) by statistically investigating the relationship between variations of the foF2 and 184 earthquakes with magnitudes M ≥ 5.0 during 1994–1999 in the Taiwan area. These results showed strong statistical evidence of ionospheric precursors. Jhuang et al. (2010) examined the TEC obtained from Global Ionosphere Maps (GIM) before 2008 Wenchuan earthquake. They found that remarkable TEC reductions appear locally around the epicenter and its conjugate points 13, 6, 5 days before the earthquake. They suggested the conjugate signature implies the seismo-generated electric field is essential. Le et al. (2011) presented a statistical study of GIM TEC before 736 M ≥ 6.0 earthquakes during 2002–2010. The results showed that the occurrence rates of abnormal days before the earthquakes depend on the magnitude, the depth, and the number of days prior to the earthquake. Parrot (2012) and Li and Parrot (2013) showed statistical analyses of ion density recorded by DEMETER satellite before large earthquakes during 2004–2010. They found that the intensities of the ionospheric perturbations increase with the earthquake magnitudes. He and Heki (2017) explored preseismic TEC changes for 32 magnitude M7.0–M8.0 earthquakes. They found the preseismic changes emerged when the background TEC was large (>50 TECU), and the precursor time duration was relatively short (<20 min). Ho et al. (2018) found both ion density and ion velocity observed by DEMETER increased anomalously before M ≥ 6.5 earthquakes in 2010. Furthermore, the density variations have a positive correlation with the perpendicular velocity. De Santis et al. (2019) analyzed the electron density and magnetic field data from Swarm three-satellite during 2014–2018. They found that the electron density and magnetic anomalies occurred from more than 2 months to some days before earthquakes.

Possible mechanisms for the generation of PEIAs have been proposed. These include (a) atmospheric acoustic and acoustic gravity waves (AGWs), which propagate through the atmosphere, reach ionospheric altitudes, and resulting in the generation of electric field disturbances and modulation of charged particles density (Liperovsky et al., 2008), (b) radon ionization, charged aerosols, and change of load resistance in the global electric circuit (Fu et al., 2019, 2020; Ouzounov et al., 2020; Pulinets & Ouzounov, 2011; Pulinets et al., 2000; Sorokin & Hayakawa, 2013), (c) a stressed rock-atmosphere–ionosphere coupling model (Kuo et al., 2011, 2014, 2015, 2018) based on experimental evidence of stressed rock currents (Freund, 2000; Freund et al., 2009), and (d) ionosphere dynamics with imposed zonal (west-east) electric field (Namgaladze et al., 2012; Zolotov et al., 2011, 2012). Other reported earthquake precursors include anomalies in geomagnetic field intensity (Chen et al., 2013; Hattori & Han, 2018), groundwater level (Igarashi et al., 1995; Yeh et al., 2015), geo-electrical activity (Uyeda et al., 2009), and seismicity (Holliday et al., 2005).

Artificial intelligence (AI) has been tried to use for predicting seismic events in recent years. The seismic indicators based on variation of parameters of Gutenberg-Richter law or possible seismic electric signals have been used to predict the magnitude of the upcoming M ≥ 4.5 earthquakes by several algorithms of machine learning (Asencio-Cortés et al., 2017; Asim et al., 2016; Panakkat & Adeli, 2007). The results showed good accuracies ≥70%. Asim et al. (2018) used 6 seismic indicators as the input to run Support Vector Regressor neural networks and Hybrid Neural Network and predicted the occurrence of earthquakes with an accuracy of up to 80%. DeVries et al. (2018) predicted the spatial distributions of aftershocks by using deep learning. The ROC (receiver operating characteristic) analysis reveals that the neural network forecast can predict aftershock locations better than the classic Coulomb failure stress change. In addition, Aji et al. (2017) applied the N-Model Neural Network Model by using the values of the TEC and Dst index at midnight (00:00 LT) to detect TEC earthquake precursors in Sumatra from December 2004 to March 2005. Akhoondzadeh Hanzaei (2018) used a Multi-Layer Perceptron (MLP) neural network to estimate TEC variations induced by the 2017 Mexico earthquake. However, the authors mentioned that it is difficult to conclude that the irregular TEC variations are related to the Mexico earthquake due to high geomagnetic and solar activities.

Sun, Xu, Huang, Zhang, Yuan, Chen, and Yan (2017) applied Long Short-Term Memory (LSTM) to forecast the ionospheric TEC for the following 24 hr based on earlier sequence data of TEC, F10.7, and ap index. They then demonstrated that bidirectional LSTM (Bi-LSTM) (Graves & Schmidhuber, 2005; Hochreiter & Schmidhuber, 1997; Sun, Xu, Huang, Zhang, Yuan, & Yan, 2017) can capture the cyclic change of TEC and gives a better TEC forecast in comparison with LSTM and multi-LSTM.

In this study, we apply deep learning to the GIM TEC data between 2003 and 2014 to identify the seismo-ionospheric precursors (SIP)of the 22 earthquakes with M ≥ 6.0 and 19 cases without the occurrence of M ≥ 5.3 earthquakes in Taiwan. In the Bi-LSTM network, the 7-day (168-hr) observational data of the feature inputs are used to predict the next hour TEC. The feature inputs in the Bi-LSTM include sequence vectors of TEC, Dst, F10.7, sunspot number (SSN), and solar emission index Lyman-α. By comparing the root-mean-square errors (RMSE) and mean absolute errors (MAE) between the predicted and observed TEC data, the anomalous TEC variation can be identified as a precursor of the upcoming earthquake.

The deep learning model is an empirical model, which employs a non-linear regression from the observed data. The predicted TEC from deep learning can be regarded as normal TEC variation, in which the earthquake causes little effect. If the observed TEC is different from the predicted TEC, the observed TEC is regarded as an abnormal variation, which is mainly caused by earthquakes in our assumption.

Furthermore, the TEC variations during the 19 periods, in which there is no earthquake of M ≥ 5.3, are also analyzed. The true and false predictions (or not) are then used to construct a confusion matrix. The accuracy will be calculated from the confusion matrix. The formal definition of confusion matrix and accuracy will be given in Section 5 and Table 2.

Spatial and Temporal Distributions of Earthquakes for Study

The TEC data are obtained from the GIM, which are retrieved from the Center for Orbit Determination Europe (CODE), University of Bern, Switzerland (http://ftp.aiub.unibe.ch/CODE/). CODE uses the GPS data from more than 200 International GNSS Service (IGS) receivers to produce the global TEC map. The global TEC data are provided with a 2.5° × 5° in latitudinal and longitudinal resolutions. The time resolution of GIM TEC data are 2 hr (1 hr) before (after) 19 October 2014.

The TEC data point closest to Taiwan is located at 120°E, 22.5°N. We then choose the earthquakes around this location within the region of 117–123°E and 19.5–25.5°N for this study. In the earthquake catalog from 2003 to 2014 in the local time frame obtained from the Central Weather Bureau of Taiwan, there are 32 earthquakes with M ≥ 6.0 and 110 earthquakes with 5.3 ≤ M < 6.0 (Figure 1). In this study, 22 cases with the presences of M ≥ 6.0 earthquakes and 19 cases without an earthquake of M ≥ 5.3 are chosen for identifying the SIP. The detail of case selection will be discussed in Section 3.2.

Figure 1. The distribution of M ≥ 6.0 earthquakes (red dots) and 5.3 ≤ M [less than] 6.0 earthquakes (green dots) during 2003–2014 in the Taiwan region. The blue cross denotes the data point of Global Ionosphere Maps total electron content in this study.

Features and Selection of Positive, Training, and Negative Periods

In this paper, the solar and ionosphere-related parameters are set as features into the network Bi-LSTM to predict TEC variations (Figure 2). The feature inputs in the Bi-LSTM include sequence vectors of TEC, SSN, F10.7, Lyman-α, and Dst index for the past 7 days. Various hyperparameters are chosen for the Bi-LSTM: the neurons are 50 and 100, batch sizes are 45, 60, 72, and 120, and iteration is 200. In Section 3.1, the daily values of the features, F10.7, SSN, and Lyman-α, are converted into hourly values, depending on the solar elevation angle. The 2-hr TEC data are interpolated to obtain the hourly TEC data. In Section 3.2, three types of time periods are chosen for study: positive periods, training periods, and negative periods. The criteria of selecting periods for study are also given. In Section 3.3, we use the relative errors of two stages to identify the seismo-ionospheric precursors.

Figure 2. The input and output features of bidirectional Long Short-Term Memory in this study.

Features: TEC, SSN, F10.7, Lyman-α, and Dst

The data of SSN, F10.7, and Lyman-α are provided with a temporal resolution of 1-day. On the other hand, the solar energy impact on the ionosphere depends on the local time. To precisely predict the hourly variation of TEC, the daily values of the feature inputs of F10.7, SSN and, Lyman-α are converted into hourly values, which depend on the solar elevation angle. The Dst data already have an hourly value. Figure 3 is a sketch of the incident solar energy flux density on the Earth's surface. Let the solar energy flux density (per unit area) in the solar radial direction be F and the elevation angle between the incident solar ray and Earth surface be β. The unit area at P is projected to Q with an area of $$1/\mathrm{sin}\,\beta $$. The incident solar energy flux density on the Earth's surface is $$F\,\mathrm{sin}\,\beta $$.

Figure 3. F is the solar energy flux density (per unit area) in the solar radial direction and β is the elevation angle between the incident solar ray and Earth's surface. The unit area at P is projected to Q with an area of 1/sinβ $1/\mathrm{sin}\,\beta $. The solar energy flux density incident on the Earth's surface is then Fsinβ $F\,\mathrm{sin}\,\beta $.

The solar elevation angle (Cooper, 1969) β can be written as: [Image Omitted. See PDF]

Here L is the geographic latitude, and the declination angle δ in degrees for a given day in Taiwan is [Image Omitted. See PDF]where d is the day of the year. The hour angle H in Taiwan is [Image Omitted. See PDF]where LT is the local hour. Then we define the parameter of solar effect as γ: [Image Omitted. See PDF]

In Equation (4), we neglect the effect of ionosphere height and the residual of the daytime activity on the TEC variation, which should be considered. The hourly data of SSN, F10.7 and Lyman-α value can be estimated as [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

Figure 4 shows a sketch for (a) the input layer and (b) the output target. Let n be the length of the sliding window in units of hours. We assume that the TEC at time $${t}_{n+1}$$ is related to the features from time $${t}_{1}$$ to time $${t}_{n}$$. Therefore the TEC, SSN_h, F10.7_h, Lyman-α_h, and Dst index from $${t}_{1}$$ to $${t}_{n}$$ are set as the input layer and the TEC at $${t}_{n+1}$$ is set as the output. For example, the features of 168 hr (7 days) are used to predict the TEC at the 169th hour. Then the input window shifts 1 hour to obtain the TEC output at $${t}_{n+2}$$. All data in the training period are run through to obtain the necessary informations for the trained model.

Figure 4. A sketch for the input layer and the output target. (a) The sequence of the input layer with time length n (n = 168 hr, 7 days). (b) A sketch of the predicted output (target) from the deep learning the time series of data in the input layer.

Using Stages for Selection of Positive, Training, and Negative Periods

Abdu et al. (2006) analyzed the power spectrum of the equatorial electro-jet (EEJ) strength from March 15 to 26 April 1994. The spectrum showed prominent spectral peaks near the periodicities at ∼6.5 and ∼14 days. In addition, Wu et al. (2020) reported that the appearance time of equatorial ionization anomaly (EIA) reveals a prominent semimonthly lunar tide of 14.77 days, which shows the existence of a nearly 15 days periodic variation in the Earth's ionosphere. In our study, we calculate the difference of errors between two neighboring stages; each stage consists of TEC data for 15 days. We then divide the error difference by the error of the first stage to obtain the relative error (RE) between these two neighboring stages. The relative errors (REs) are used for detecting earthquake precursors. With the choice of two 15-day neighboring stages, the effects of the semimonthly lunar tide on the obtained relative errors will be significantly reduced. The seasonal effects and the solar activity effects can also be reduced by the use of the “relative errors (REs)” between two neighboring stages for prediction.

After preparing the features of input, we select three types of periods for the present study: positive periods, training periods, and negative periods (Figure 5a). The data length of a positive (negative) period consists of 45(60) days, which are divided into three(four) stages, each with 15 days. The relative errors between two stages are calculated. The starting date of Stage 1 and ending date of Stage 2 for positive periods and negative periods are listed respectively in Tables S3 and S4 of the Supporting Information S1. Table S3 in Supporting Information S1 also lists the earthquake magnitude in the positive period. Table S5 in Supporting Information S1 lists the starting and ending dates of Stage 1 for the training periods.

Figure 5. (a) A sketch of analyzing periods. Red, blue, and green blocks represent positive, training, and negative periods, respectively. The positive periods are 1–30 days before M ≥ 6.0 earthquake and without other M ≥ 6.0 earthquakes. The white block before the positive period shows the 15-day period without M ≥ 6.0 earthquake. The training and negative periods are 45 and 30 days without M ≥ 5.3 earthquakes, respectively. The white blocks before and after the training and negative periods show the periods without M ≥ 5.3 earthquakes. The yellow areas are the event day (Day 0). The number under the positive and negative blocks shows the days to the event day. (b) The distribution of positive periods (red horizontal line, 1–30 days before M ≥ 6.0 earthquakes), training model periods (blue horizontal line), and negative periods (green horizontal line) in this study. The black bars denote the occurrence of M ≥ 6.0 earthquakes, and the gray bars denote the occurrence of 5.3 ≤ M [less than] 6.0 earthquakes.

Figure 5b shows the distribution of positive periods (red horizontal bars), training periods (blue horizontal bars), and negative periods (green horizontal bars). We first choose the positive periods. Each positive period consists of 45 days (15 days for each of Stage 0, Stage 1, and Stage 2) before an M ≥ 6.0 earthquake (black bars in Figure 5b); there is no M ≥ 6.0 earthquake during the 45 days. A training period of 75 days for constructing a predicting model is then chosen, in which there is no M ≥ 5.3 earthquake (15 days for Stage 0, 45 days for Stage 1, and 15 days for Stage 3). The training period is designed for a longer data length with 45 days in Stage 1 in order to have a more comprehensive deep learning. The training period has a total length of 75 days, while the negative period has 60 days in total. We select all training periods first and the negative period after, because it is harder to find a 75-day slot than a 60-day slot. During a negative period of 60 days (15 days for each of Stage 0, Stage 1, Stage 2, and Stage 3), it has been verified that there is also no M ≥ 5.3 earthquake. Figure 5b shows the temporal distribution of time periods of the three types. To avoid mutual interference among the periods, a buffer interval of 15-day slot (Stage 0 or Stage 3) is included. There is no temporal overlap between two training periods, between two positive periods, or between two negative periods, but the time overlap between one positive period and one negative period is allowed.

For a training period, only the TEC data in Stage 1 (45 days) are used to train a model. The trained models are then used to predict the TEC in Stage 1 and Stage 2 for a positive or negative period.

Table 1 lists the numbers of positive periods, negative periods, and training models as a function of earthquake magnitude based on the criteria of the associated day length of the period. To guarantee the absence of the PEIA in the training and negative periods, the magnitude of an earthquake is required as small as possible. In addition, the number of negative periods needs to be comparable to the number of positive periods to avoid data imbalance (He & Garcia, 2009). If we choose M < 6.0 as the criteria of negative periods, the number of negative periods is 67, which is three times the number of positive periods (22). These two groups of periods with a skewed proportion are called imbalanced data, which is not suitable for using accuracy to check the performance of a model. However, the number of negative periods for M < 5.0 is only three. It will lead to imbalanced data again if we use M < 5.0 as the criteria of negative periods. Therefore, we choose M < 5.3 as the criterion for the training and negative periods. Moreover, the influence of M < 5.3 earthquakes on the ionosphere is small enough for the use of generating the trained models. The influence of the M ≥ 6.0 earthquakes on the ionosphere is large enough to be detected by the trained models. The 22 cases with the presences of M ≥ 6.0 earthquakes are first chosen. Second, 13 periods without an earthquake of M ≥ 5.3 are chosen as the training period. Finally, 19 cases without an earthquake of M ≥ 5.3 are chosen for identifying the SIP (Figure 5b).

Table 1 The Numbers of Positive Periods, Negative Periods, and Training Models as a Function of Earthquake Magnitude Based on the Selected Criteria as Discussed in Figure 5a and Section 3.2

Magnitude	Positive periods	Negative periods	Training models
5.0	16	3	3
5.1	23	14	7
5.2	28	18	10
5.3	34	19	13
5.4	33	22	19
5.5	32	31	24
5.6	31	34	26
5.7	33	42	29
5.8	32	50	33
5.9	26	61	39
6.0	22	67	39

The predicted TEC from the trained model can be regarded as normal TEC variation. If the observed TEC is different from the predicted TEC, the observed TEC is regarded as an abnormal variation, which could be mainly caused by an earthquake in our assumption. Two types of error metrics are taken into account: RMSE and MAE. The relative errors (REs) between the errors at Stage 1 (RMSE₁, MAE₁) and the errors at Stage 2 (RMSE₂, MAE₂) can be expressed as: [Image Omitted. See PDF]and [Image Omitted. See PDF]

For illustration purpose, let RMSE_sea2 be the RMSE associated with the seasonal effect at Stage 2, RMSE_sea1 be the part of RMSE associated with the seasonal effect at Stage 1, and RMSE_eqk2 be the RMSE associated with the earthquake precursor at Stage 2. The Equation 6(a) becomes [Image Omitted. See PDF]where $${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{s}\mathrm{e}\mathrm{a}2}\approx {\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{s}\mathrm{e}\mathrm{a}1}$$. The seasonal effect is then highly reduced and only the term $${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{e}\mathrm{q}\mathrm{k}2}$$ remains.

We compare the relative error for Stage 1 and Stage 2 of a positive or negative period. If the relative error exceeds the given threshold, the trained model will declare that the previous 30 days before Day 0 as “ positive”. Otherwise, the trained model will declare the 30 days before Day 0 as “negative”.

The Relative Errors for Identification of Earthquakes Precursors

This study aims to detect the seismo-ionospheric precursors by calculating the normalized difference (relative error) between the predicted TEC variation obtained from the trained models and the observed TEC data.

A pair of negative and positive periods are presented in Figure 6. For the negative period N3 from 16 August 2004 to 14 September 2004 shown in Figure 6a, the blue line is the predicted TEC value by the trained model T5. The green line is the observed TEC, and the black shadow is the absolute deviation between predicted TEC and observed TEC above the horizontal axis. The RMSE in the 30 days period is 0.0740. The associated error metric RMSE (MAE) at Stage 1 and 2 are 0.0690 and 0.0787 (0.0525 and 0.0596), respectively. The relative error $${RE}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=0.1408\,\left({RE}_{\mathrm{M}\mathrm{A}\mathrm{E}}=0.1351\right)$$. Figure 6b shows a predicted and observed TEC variation in a positive period from 15 September to 14 October 2004 (P4). An earthquake M7.1 occurred at longitude 121.85°E and latitude 24.16°N with depth 91 km on 15 October 2004 (Day 0). The blue line is the predicted TEC value by the trained model T5. The red line is the observed TEC, and the black shadow is the absolute deviation between predicted TEC and observed TEC. The RMSE in the 30 days period is 0.0964. The associated error metric RMSE (MAE) in stage 1 and 2 are 0.0720 and 0.1158 (0.0552 and 0.0903), respectively. The relative error $${RE}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=0.6085\,\left({RE}_{\mathrm{M}\mathrm{A}\mathrm{E}}=0.6347\right)$$. The relative error of the positive period is four times larger than that of the negative period shown in Figure 6. Stage 2 of Figure 6b shows a higher deviation and persists longer than stage 1. Compared with the deviations at stages 1 or 2 in Figure 6a, this deviation is also higher than the deviations at stages 1 or 2 in Figure 6a. Therefore, the relative error for the positive period in Figure 6b is much higher than that for the negative period in Figure 6a.

Figure 6. A pair of examples of total electron content (TEC) data processing results. (a) TEC data in normal condition (negative period) from 16 August to 14 September 2004. The blue line shows the predicted TEC from trained model T5, the green curve shows the observed TEC, and the black shadow is the absolute deviation between predicted TEC and observed TEC. (b) TEC data before the occurrence of M7.1 earthquake (positive period) from 15 September to 14 October 2004. The blue curve shows the predicted TEC from trained model T5, the red curve shows the observed TEC, and the black shadow is the absolute deviation between predicted TEC and observed TEC.

In Figure 6, the peaks and dips are due to the daily periodic TEC variation. It is interesting to note that the maximum differences between the predicted and the observed TEC data often occur around these peaks and dips. We call the difference between the predicted and the observed TEC data as “the model error” associated with the trained model. The model error is affected by the semimonthly variation, seasonal variation, and solar activity variation of the TEC data, and is not appropriate for detecting earthquake precursors. Therefore, we introduce the relative (normalized) error in Equations 6a and 6b, which can reduce the periodic effects.

Effect on Detecting Ionospheric Precursors Associated With M ≥ 6.0 Earthquakes by the Ionospheric Disturbances From a Geomagnetic Storm, Typhoon, and Solar Activity

There are sixteen (nine) cases with geomagnetic storms and four cases with typhoons during the negative (positive) periods. Among the 16 negative cases with geomagnetic storms, nine of them are with long-term geomagnetic storms (lasting more than 7 days). Figures 7 and 8 show the predicted TEC from the trained model T5 and observed TEC under the influences of a geomagnetic storm. Figure 9 shows the predicted TEC from the trained model T5 and observed TEC under the influences of a typhoon. In addition, we show two cases in Figure 10 with different solar activity (different years) but in the same season (Fall) to discuss the prediction capability of the trained models.

Figure 7. An example of predicted and observed total electron content (TEC) variations from 15 August to 13 September 2004 with a geomagnetic storm. (a) The predicted TEC data from trained model T5 are produced by the bidirectional Long Short-Term Memory models with the hourly features of TEC, SSNh, F10.7h, Lyman-αh, and Dst. The blue line is the predicted TEC from model T5, the green curve shows the observed TEC, and the black shadow is the absolute deviation between predicted TEC and observed TEC. (b) The Dst index variations from 15 August to 13 September 2004. The root-mean-square errors (RMSE) and the relative error RERMSE are 0.0715 and 0.1191, respectively. The presence of a geomagnetic storm from 30 August to 4 September 2004 does not affect the capability of predicting the TEC variations.

Figure 8. An example of predicted and observed total electron content (TEC) variations from 1 October to 30 October 2013 with the presence of two geomagnetic storms. (a) The predicted TEC data from trained model T5 are produced by the bidirectional Long Short-Term Memory models with the hourly features of TEC, SSNh, F10.7h, Lyman-αh, and Dst. The blue line is the predicted TEC from model T5, the red curve shows the observed TEC, and the black shadow is the absolute deviation between predicted TEC and observed TEC. (b) The Dst index variations from 1 to 30 October 2013. The root-mean-square errors (RMSE) and the relative error RERMSE are 0.1292 and 0.1869, respectively. The presence of two geomagnetic storms are from 2 to 7 October (Storm 1) and from 9 to 12 October (Storm 2), 2013.

Figure 9. An example of total electron content (TEC) processing results from 15 days before to 15 days after a typhoon arrived Taiwan from 29 August to 27 September 2008. (a) The blue line shows the predicted TEC from the trained model T5, the green curve shows the observed TEC, and the black shadow is the absolute deviation between predicted TEC and observed TEC. (b) The Dst index variation from 29 August to 27 September 2008. The root-mean-square errors (RMSE) and the relative error RERMSE are 0.0873 and 0.0515 from input features of TEC, SSNh, F10.7h, Lyman-αh, and Dst.

Figure 10. (a) The observed (green line), predicted (blue) total electron content (TEC) data, and their absolute deviation (black shadow) for negative period N3 from 16 August to 14 September 2004. (b) The observed and predicted TEC data, and their absolute deviation for period N10 from 15 August to 13 September 2009. The predicted TEC values are obtained from the trained model T5. The solar activity is high for N3 in 2004 and quiet for N10 in 2009.

Geomagnetic Storm

Figure 7a shows the predicted (blue line) and observed (green line) TECs profiles during the negative period N3 from August 15 to 13 September 2004. The absolute deviation between predicted and observed TECs are plotted as black shadow above the horizontal axis. Figure 7b shows the Dst variation during the period. A geomagnetic storm occurred from 30 August to 4 September, as indicated by the red line. The hourly input features of the trained model T5 in Figure 7a are TEC, SSN_h, F10.7_h, Lyman-α_h, and Dst, and the RMSE is 0.0715. The daily maximum of the observed TEC ranges from 25 to 55 TECu from 30 August to 4 September, while the geomagnetic storm may affect the TEC variation. By excluding this period, the daily maximum of the observed TEC ranges from 35 to 80 TECu. Meanwhile, the predicted TEC also has a similar variation. Here, the RMSE₁ is 0.0674 at Stage 1, and the RMSE₂ is 0.0754 at Stage 2. The relative error RE_RMSE (=0.1191) is lower than 0.15, and model T5 calls these 30 days “Negative.” Likewise, the RE_RMSE obtained from T7 and T10 is 0.0792 and 0.0429, respectively. T7 and T10 also call this period as “Negative.” Therefore, the majority votes “No” in this case. The presence of a geomagnetic storm from 30 August to 4 September 2004 does not affect the capability of predicting the TEC variations.

Similarly, Figure 8a shows the predicted (blue line) and observed (red line) TECs profiles during the positive period P22 from 1 to 30 October, 2013. An M6.4 earthquake occurred at 121.4°E and latitude 23.6°N on 31 October 2013, with a depth of 15 km. The absolute deviation between predicted TEC and observed TEC is plotted as black shadow. Figure 8b shows the Dst variation during the period. Two geomagnetic storms (Storm1 and Storm 2) occurred at Stage 1 from 2 to 6 October and from 9 to 13 October, as indicated by the red line. The hourly input features of the trained model T5 in Figure 8a are TEC, SSN_h, F10.7_h, Lyman-α_h, and Dst. The deviation (black shadow) in Storm 1 ranges from 0.2 TECu to 15.1 TECu. The deviation in Storms 2 ranges from 0.1 TECu to 12.2 TECu. For each storm, the most significant deviation occurs in the geomagnetic storm's main phase. The deviation at Stage 2 are large and can be more than 20 TECu. The absolute values of deviation at Stage 1, which are likely associated with the geomagnetic storms, are less than the absolute values of deviation at Stage 2, which are likely associated with the earthquake. The RMSE₁ is 0.1178 at Stage 1, and the RMSE₂ is 0.1398 at Stage 2. The RMSE from October 1 to 30 is 0.1292. The relative error RE_RMSE (= 0.1869) is higher than 0.15, and model T5 calls these 30 days “Positive.” Likewise, the RE_RMSE obtained from T7 and T10 is 0.1501 and 0.3728, respectively. T7 and T10 also call this period as “Positive.” Therefore, the majority votes “Yes” in this case. The presence of two geomagnetic storms at Stage 1 from 2 to 7 and from 9 to 12 October, 2013 seems to not affect the capability of predicting the TEC variations for earthquake precursors.

Typhoon

Figure 9 shows the predicted and observed TECs from 29 August to 27 September 2008. The hourly input features of the trained model T5 in Figure 9a are TEC, SSN_h, F10.7_h, Lyman-α_h, and Dst during the passage of a typhoon. Chou et al. (2017) reported the observations of the concentric traveling ionospheric disturbances triggered by a typhoon. An intense intensity typhoon, “Sinaku,” was formed at 125.7°E, 16.7°N, approached Taiwan in midnight (0150LT) on 14 September 2008, and caused strong wind and heavy rain in Taiwan until the early morning (0800LT) of 15 September 2008. Here, the RMSE₁ is 0.0896 at Stage 1, and the RMSE₂ is 0.0849 at Stage 2. The relative error RE_RMSE (= 0.0515) is lower than 0.15, and model T5 calls these 30 days “Negative.” The RE_RMSE obtained from T7 and T10 is 0.2004 and 0.1243, respectively. T7 (T10) calls this period “Positive (Negative).” Therefore, the majority vote this case “No”. The results show that the presence of a typhoon from September 11 to 15, 2008 seems to not affect the capability of predicting the TEC variations associated with earthquakes.

Solar Activity

Figure 10 shows the observed TEC (green line) and the predicted TEC (blue line) profiles during (a) the negative period N3 from 16 August to 14 September in 2004, and (b) the negative period N10 from 15 August to 13 September in 2009. The absolute deviations between predicted and observed TECs are plotted as black shadows in each panel. The negative case N3 (N10) occurs during the solar high (quiet) period. The predicted TECs for both cases are obtained from the trained model T5. The observed TEC data range from 10 to 75 TECu in the solar active period (N3), and from 5 to 30 TECu in the solar quiet period (N10). The predicted TECs show a range from 10 to 80 TECu in 2004 (solar active period) and 5 to 32 TECu in 2009 (solar quiet period). The trained model predicts a larger (smaller) TEC variation during solar active (quiet) period. The RE_RMSE are 0.1408 and 0.0942 in Figures 10a and 10b, respectively. These small relative errors show that the trained model is stable during both periods, even with different solar activities.

Figure 11 shows the relative errors RE_RMSE obtained from three classifiers, T5, T7, and T10, for the four negative cases in Figures 7, 9, and 10. Only the relative error obtained from T7 for the typhoon case exceeds the threshold of 0.15. For the rest of the six negative cases with short-term geomagnetic storms and three cases with typhoons, all majorities from three classifiers vote “No” to these cases.

Figure 11. The relative errors REsRMSE are obtained from three classifiers, T5, T7, and T10 for cases in Figures 7, 9, and 10. The black circle denotes the case with a geomagnetic storm in Figure 7. The blue diamond denotes the case with a typhoon in Figure 9. The triangles denote the cases for solar active period (red) and for solar quiet period (gray) in Figure 10. The black dashed line is the threshold 0.15.

In summary, the trained models (three classifiers) can identify whether the TEC variations are related to an earthquake or not, even under the influence of a geomagnetic storm, the possible disturbance of the gravity waves from a typhoon, and different solar activities.

Ensemble Learning: Classifiers, Vote and Random Seeds

Figure 12 shows the distribution of the relative errors in (a) RMSE and (b) MAE from all positive and negative periods obtained from the 13 trained models. The relative errors range from $$1.14\times {10}^{-4}$$ to 6.83. Some trained models (e.g., T1, T2, T12, and T13) have wide distributions in associated relative error metrics. These broad distributions can underline a distorted TEC prediction. Therefore, we apply ensemble learning to determine the trained models for this study. Three classifiers will be selected based on the following discussion and will be used for a majority vote. The majority vote is used to identify the TEC variations for having a “positive” or “negative” association with an M ≥ 6.0 earthquake. In this Section, we present the procedure of selecting classifiers and obtaining the majority vote. The hyperparameters are chosen as follows: the number of neurons is 100, batch size is 45, and iteration is 200.

Figure 12. The distribution of the relative errors in (a) root-mean-square errors, and (b) mean absolute errors from all positive and negative periods obtained from 13 trained models. The upper and lower boundaries of the boxes denote the upper and lower quartiles. The red lines in the boxes are the median. The red crosses are the outliers.

Selection of Three Classifiers

Figure 13 shows the process of calculating the relative errors of the example of model T1 and the other models have been evaluated using the same scheme. We calculate the error metrics between predicted and observed TEC of each trained model for each period. For the events before October 2014, each trained model has 22 pairs of RE_RMSE and RE_MAE for positive periods and 19 pairs of RE_RMSE and RE_MAE for negative periods. The median $${\tilde{RE}}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$$ ($${\tilde{RE}}_{\mathrm{M}\mathrm{A}\mathrm{E}}$$) of the 22 RE_RMSE (RE_MAE) of positive periods and of the 19 RE_RMSE (RE_MAE) of negative periods are obtained. The difference $${{\increment}\tilde{RE}}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$$ ($${\increment}{\tilde{RE}}_{\mathrm{M}\mathrm{A}\mathrm{E}}$$) between positive medium $${\tilde{RE}}_{\mathrm{R}{\mathrm{M}\mathrm{S}\mathrm{E}}_{P}}$$ ($${\tilde{RE}}_{{\mathrm{M}\mathrm{A}\mathrm{E}}_{P}}$$) and negative medium $${\tilde{RE}}_{{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{N}}$$ ($${\tilde{RE}}_{{\mathrm{M}\mathrm{A}\mathrm{E}}_{N}}$$) is given by: [Image Omitted. See PDF] [Image Omitted. See PDF]

Figure 13. A sketch of model validation for training period 1 (T1). There are 22 positive periods (P1–P22) and 19 negative periods (N1–N19). RE is the relative error in percentage as shown in equations 6a and 6b, RE∼ $\tilde{RE}$ is the median of RE, ΔRE∼ ${\Delta}\tilde{RE}$ is the median difference between positive and negative RE∼ $\tilde{RE}$, and ΔRE∼total ${\Delta}{\tilde{RE}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$ is the summation of ΔRE∼ ${\Delta}\tilde{RE}$ of root-mean-square errors (RMSE) and mean absolute errors (MAE). The subscript RMSE and MAE show the types of error, P and N denote the positive and negative periods, respectively.

The subscript P is for a positive period and N is for a negative period. The total error of each model is the sum of $${\increment}{\tilde{RE}}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$$ and $${\increment}{\tilde{RE}}_{\mathrm{M}\mathrm{A}\mathrm{E}}$$: [Image Omitted. See PDF]

Three trained models with the largest $${\increment}{\tilde{RE}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$$ are selected as the classifiers for the ensemble learning process. T5, T10, and T7 are selected as the classifiers before October 2014.

Vote by Selected Trained Models

Three trained models with the best discrimination between the relative errors of positive period and of negative period are selected as the classifiers. These three selected classifiers declare that each period is positive or negative, and then a majority voting is adopted. Three sets of predicted TEC and its RE_RMSE and RE_MAE are obtained by the three classifiers for each period. A set of the threshold from 5% to 30% is used to test the results. If the RE_RMSE or RE_MAE of one period is larger than a given threshold, this classifier declares that this period is “positive”. If two or three classifiers announcing “positive” for a given period, the final decision is positive, or “YES”, that is, the classifiers predict that an M ≥ 6.0 earthquake will occur from Day 0 to Day 14. Otherwise, the final decision is negative, or “NO”, that is, the classifiers predict that no M ≥ 5.3 earthquake will occur on Day 0. The predicted final decisions for all positive and negative periods are then examined whether they are true or false in comparison with the truly observed earthquakes. The highest accuracy is obtained with the threshold of 15% for the period from January 2003 to October 2014. In addition, the effect on detecting ionospheric precursors associated with M ≥ 6.0 earthquakes by the ionospheric disturbances from a geomagnetic storm, typhoon, and solar activity are discussed in the Supporting Information S1.

Selection of Length of Sliding Windows

Bi-LSTM predicts data based on the features of the prior sequence window. In this Section, we survey the influence of window length from 1 day to 8 days in the input layer of the Bi-LSTM models. Figure 14a shows the relative error for the same periods N3 (green bars) and P4 (red bars) shown in Figure 6 and the same hyperparameters but with different lengths of sliding windows. For the length of sliding window with 1–3 days, the relative error of N3 is larger than the positive period P4. However, when the length of the sliding window is longer than 3 days, the relative error of P4 becomes larger than that of the negative period N3. The figure shows that the 7-day sliding window length has the best discrimination between the negative period N3 and the positive period P4.

Figure 14. (a) The relative error are plotted as a function of length of sliding window from 1 to 8 days for the negative period N3 (green bar) and for the positive period P4 (red bar), those are discussed in Figure 6. The red bars are the RERMSE for P4, and the green bars are the RERMSE for N3. (b) The distribution of accuracy with different hyperparameters in different window lengths from January 2003 to October 2014. The upper and lower boundaries of the boxes denote the upper and lower quartiles. The red lines in the boxes are median and the green diamonds are the average.

On the other hand, we check the accuracy distribution as a function of window length with various hyperparameters. For a given window length, we input the hyperparameters of 50 and 100 neurons per layer with batch sizes of 45, 60, 72, and 120 to Bi-LSTM. There are 8 sets of hyperparameters for each sequence window. The accuracy is calculated from the confusion matrix of true and false prediction for all positive and negative periods. Figure 14b shows the accuracy distribution of these 8 sets of hyperparameters. In general, most of the accuracy is distributed in 60%–75% for the data before October 2014. Moreover, the interquartile range (IQR) for each window length ranges from 4% to 12%. For narrower IQRs, the data of the distribution are more concentrated, which suggests that the Bi-LSTM models with various hyperparameters in this study are more stable. The average accuracies of the models with the window length of 7 days (72.62%) are higher than those with other window length (64.29%–68.15%) before October 2014. Therefore, the 7-day window is chosen in the input layer of Bi-LSTM for later calculations. As discussed in Section 3.2, Abdu et al. (2006) observed one of the prominent spectral peaks near ∼6.5 days of EEJ strength. This can provide a possible physical reason for 7 days as the best length of sliding window. In addition, the set of other hyperparameters with the best accuracy is a neuron of 100, batch size of 45, iteration of 200, and activation function tanh.

Effect of Random Seeds and Stability of Prediction Scheme

As mentioned before, this study attempts to identify the SIP by using the normalized difference (relative error) between the predicted TEC from trained models and the observed TEC data. Figure 15 illustrates the distribution of the accuracies obtained for different error metrics ($${RE}_{\mathrm{e}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}},{RE}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}},{RE}_{\mathrm{M}\mathrm{A}\mathrm{E}}$$) with same three classifiers, same set of hyperparameters, and different random seeds for the period from January 2003 to October 2014. Each box is estimated from the same set of hyperparameters and input features discussed in Section 4.3 but with 10 different random seed numbers. The seed value of the random number has an impact on the accuracy as indicated by the vertical spread of accuracy in the Figure 15.

Figure 15. The distribution of accuracies obtained for different error metrics (REeither,RERMSE,REMAE ${RE}_{\mathrm{e}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}},{RE}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}},{RE}_{\mathrm{M}\mathrm{A}\mathrm{E}}$) with the same three classifiers, the same set of hyperparameters, and different random seeds for the period from January 2003 to October 2014. Each box-and-whisker represents the distribution of the accuracies with the same hyperparameters but 10 different random seeds. Different boxes are estimated by different error metrics (REeither,RERMSE,andREMAE ${RE}_{\mathrm{e}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}},\,{RE}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}},\,\mathrm{a}\mathrm{n}\mathrm{d}\,{RE}_{\mathrm{M}\mathrm{A}\mathrm{E}}$). The red lines and green diamonds in the boxes show the median and average of each data set.

The accuracy is in the range of 61.90%–78.05%. The IQRs of accuracy for different error metrics $${RE}_{\mathrm{e}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}},\,{RE}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}},\,{RE}_{\mathrm{M}\mathrm{A}\mathrm{E}}$$ are 4.76%, 4.76%, and 7.14% respectively. This suggests that the trained models by $${RE}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$$ are more stable. The best accuracy reaches 78.05% before October 2014 as estimated by 𝑅𝐸_{𝑀𝐴𝐸} even if it shows a larger IQR equal to 7.14%. The accuracies estimated by $${RE}_{\mathrm{e}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}}$$ are more stable than the results by $${RE}_{\mathrm{M}\mathrm{A}\mathrm{E}}$$, and have higher accuracy values than the results by $${RE}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$$. Figure 15 shows that our scheme with the chosen hyperparameters and the use of $${RE}_{\mathrm{e}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}}$$ is more stable and suitable for the prediction of TEC values and earthquake occurrence.

Results and Discussions: Confusion Matrix, Accuracy, TPR, and NPV

In Table 2, we list the prediction results before October 2014 in the form of a confusion matrix. A confusion matrix is a contingency table with “actual” and “predicted” numbers. In the table, “EQ” and “no EQ” indicate whether or not the earthquake occurs in observation; “Predict YES” and “Predict NO” indicate the model predictions. TP (true positive) stands for successful predictions of the occurrence of earthquakes; TN (true negative) for successful prediction of no earthquake; FP (false positive) for false prediction of earthquake occurrence but earthquakes do not occur as predicted; FN (false negative) for false prediction of no earthquake but earthquakes occur. In the matrix, from January 2003 to October 2014, 22 positive (TP) and 9 negative (FP) periods are announced as positive (predict “YES”). In contrast, there are 0 positive (FN) and 10 negative (TN) periods being announced as negative (predict “NO”). The accuracy ($$\mathrm{A}\mathrm{C}\mathrm{C}=\frac{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}+\mathrm{F}\mathrm{N}}$$) is 78.05%.

Table 2 Confusion Matrix of Prediction for the Period Between January 2003 and October 2014

Earthquake
Prediction	EQ	No EQ
Predict YES	22 (TP)	9 (FP)
Predict NO	0 (FN)	10 (TN)

TP is true positive; FP is false positive; FN is false negative; TN is true negative. $$\mathrm{A}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}=\frac{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}+\mathrm{F}\mathrm{N}}$$

We study the nine FP cases from January 2003 to October 2014. In all these cases, there are geomagnetic storms or substorms lasting for more than 7 days (the sliding window for the training models). We refer to the hourly Dst data and the daily average of Ap geomagnetic indices data from the International Service of Geomagnetic Indices (ISGI) to define a geomagnetic storm. A geomagnetic storm is defined as a local minimum value of the Dst index below −30 nT (Kamide et al., 1998; Loewe & Prölss, 1997). The lasting time of a geomagnetic storm or a substorm is defined from when the value of the Dst index is smaller than zero and back to zero.

The reason for FP in these cases could be that the training models cannot properly resolve the TEC variations under the influence of storms or substorms lasting for more than the sliding window. An example of an FP case (7 May to 5 June 2007) is shown in Figure 16. In this case, the relative errors given by both the classifiers T5 and T10 exceed the threshold, leading to a “positive” prediction. Figure 16a shows the observed TEC (green line) and predicted TEC from the trained model T10 (blue line). At Stage 1, the RMSE₁ and MAE₁ are 0.0842 and 0.0677, respectively. At Stage 2, the RMSE₂ and MAE₂ are 0.1097 and 0.0905, respectively. The relative errors RE_RMSE = 0.3028 and RE_MAE = 0.3375 are greater than the threshold (15%). Figure 12b shows the geomagnetic index Dst variation. At Stage 2 (22 May to 6 June 2007), a geomagnetic storm lasted for 8 days beginning from 22–30 May 2007. The large relative error, RE_RMSE = 0.3028, is caused by the long duration (8 days) geomagnetic storm. If all the cases with long-term storms or substorms are removed, the number of FP cases becomes 0 before October 2014. The accuracy becomes 100%.

Figure 16. The variations of (a) predicted and observed total electron content (TEC), and (b) geomagnetic index Dst during the negative period N8 from 7 May to 5 June 2007. The blue line is the predicted TEC from trained model T10. The green line is the observed TEC. The black shadow in (a) is the absolute deviation between predicted TEC and observed TEC. The large relative error, RERMSE = 0.3028, is caused by the long duration (8 days) geomagnetic storm.

Table 3 lists the statistical indices: accuracy (ACC), true positive rate (TPR) (hitting rate, H), false positive rates (F), the rate of random agreement ($${p}_{e}$$), and Cohen's kappa coefficient (κ) (Cohen, 1960). The “corrected” in Table 3 indicates that the cases with long-term storms or substorms are excluded. Cohen's kappa coefficient is an estimate of the performance of an algorithm in classifying an object compared to a classification due to chance. The formula is given as: [Image Omitted. See PDF]where $${p}_{0}$$ is the observed proportion agreement ($${p}_{0}$$ equals to ACC), and $${p}_{e}$$ is the rate of random agreement. As listed in Table 3, after excluding cases with the long-term storms or substorms, the false alarm rate is reduced from 47.37% to 0 for the cases before October 2014.

Table 3 The Statistical Indices: Overall Accuracy (ACC), True Positive Rate (Hit Rate, H), False Alarms Rate (F), the Rate of Random Agreement (p_e), and Cohen's Kappa Coefficient (κ) of All Events Before October 2014

Index
Events	ACC (%)	H (%)	F (%)	pe	κ
All events before October 2014	78.05	100	47.37	0.52	0.54
Selected eventsa	100	100	0	1	1

^aThe cases with long-term storms or substorms are excluded.

All the selected M ≥ 6.0 earthquakes between 2003 and 2014 are successfully predicted, so the TPR ($$\mathrm{T}\mathrm{P}\mathrm{R}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}$$) is 100%. On the other hand, the results show no occurrence of M ≥ 5.3 (also no M ≥ 6.0) earthquakes on Day 0 when the trained models predict “NO”. The negative predictive value (NPV) ($$\mathrm{N}\mathrm{P}\mathrm{V}=\frac{\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{N}}$$) is 100% in our study. However, 9 of 19 negative cases are predicted “YES”, that is, there are nine “false positive (FP)” cases. The false prediction may be caused by the geomagnetic storms lasting more than 7 days during these periods.

Summary

In this paper, we apply deep learning to the ionospheric TEC data between January 2003 and October 2014 to detect the SIP of M ≥ 6.0 earthquakes in Taiwan based on the hourly data of TEC, Dst, F10.7_h, SSN_h, and Lyman-α_h. The impact on the ionosphere from the solar activities depends on the solar elevation angle, which varies with local time. The sun-related input features, F10.7, SSN, and Lyman-α, are corrected by using the solar elevation angle to obtain the hourly values.

Three of 13 trained models are selected to predict the TEC data for 22 positive periods. Each positive period consists of 45 days before an earthquake of M ≥ 6.0. Each of the 19 negative period consists of 60 days without any earthquake of M ≥ 5.3. The results show that all positive cases are successfully predicted, giving a TPR of 100%, while 9 of 19 negative cases are predicted “YES”. The nine false predictions may be caused by the geomagnetic storms lasting more than 7 days during these periods. Overall, a high accuracy of 78.05% is obtained. Meanwhile, there is no occurrence of M ≥ 5.3 (so, also no M ≥ 6.0) earthquake when the trained models predict “NO”; the NPV is also 100%. In the next paper, the current procedures of detecting SIP will be used to examine the earthquakes after 2015 in Taiwan, in which the GIM data have a 1-hr time resolution. The current study cannot predict the exact occurrence day of an M ≥ 6.0 earthquake. However, this study can provide a clue about an impending earthquake. In the future, we plan to integrate the multiple observations of ground-based data from Gamma-ray and Radon (Fu et al., 2019, 2020; Ouzounov et al., 2020) and anomalies in geomagnetic field intensity (Chen et al., 2013; Hattori & Han, 2018) to determine the spatial range of epicenter of an impending earthquake in Taiwan.

Acknowledgments

This research is supported by the Ministry of Science and Technology (Grant No. MOST 110-2111-M002-017, MOST 110-2740-M-001-001, MOST 109-2111-M-002-001, MOST 108-2111-M-492-001 and MOST 109-2116-M-002-031) and the Central Weather Bureau (Grant No. MOTC-CWB-109-E-01). We thank the National Center for High-performance Computing (NCHC) for providing computational and storage resources. We thank Dr. MC Liang for the helpful suggestions.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

The Global Ionosphere Map data are provided by the Center for Orbit Determination Europe (CODE, http://ftp.aiub.unibe.ch/CODE/), University of Bern, Switzerland. The data of SSN, F10.7, Lyman-α, and Dst index are obtained from Goddard Space Flight Center Space Physics Data Facility (https://omniweb.gsfc.nasa.gov/form/dx1.html?fbclid=IwAR0iPdGDQFk_pIhs9BcLV9BiGBOCTvk9mvYP474LFqEnCVvguWmYecNUFDM).