CSES and DEMETER satellites were both launched specifically with the purpose of monitoring earthquakes. The DEMETER satellite was launched in 2004 by France and ceased data collection at the end of 2010, with an operating time of 6.5 years (Parrot, 2006). Data along more than 30,000 orbits were obtained, which provided a solid data foundation for the research in earthquake monitoring and ionospheric physics. During the operation of the DEMETER satellite, preliminary preparations for the CSES program were initiated and successfully launched on February 2, 2018 (Shen et al., 2020). However, compared to DEMETER, CSES is designed for a lower orbit altitude of 507 km, which is closer to the ionospheric peak region; the lifting and lowering nodes are at 14:00 and 02:00 local time, also known as day-side or night-side orbits, and the day-side observation time is when the peak ionospheric electron density occurs; the satellite revisiting time is set to 5 days, which is denser than the 16 days revisiting time of the DEMETER satellite.

Of the scientific payloads on CSES, the LAP is the in situ space plasma detection device (Liu et al., 2018). LAP can measure electron density (Ne) in the range of 5 × 10²–1×10⁷ cm⁻³, and electron temperature (Te) in the range of 500–10,000 K, with a relative accuracy of 10%, allowing the analysis of Ne and Te to study space plasma physical phenomena and ionospheric changes caused by earthquakes, magnetic storms, and other events. CSES data, including LAP electron density and electron temperature from August 1, 2018 to June 22, 2020, were utilized in this study. During this period, 6004 EQs with magnitudes ≥4.8 were recorded (USGS: http://www.usgs.gov). Due to the impact of the magnetic storm and the Sumatra earthquake from November to December 2004, we excluded the data from the DEMETER satellite in 2004 and started our study using the data from 2005 onwards. The DEMETER data used in this study included LAP electron density and electron temperature from 2005 to the end of 2010. During this period (about 6 years), 20,727 EQs with magnitudes ≥4.8 were recorded. The K _p index was referenced to avoid the effect of solar and magnetic activities (K _p > 3) in this paper.

To avoid mixing pre-and post-seismic effects, we deleted the aftershocks from the list of earthquakes in this study (Yan et al., 2017). In our study, all the seismic events in the area of 2° × 2° centered on the epicenter and within 15 days after the given earthquake were considered aftershocks. The choice of 15 days as the largest anticipation time of the pre-earthquake anomaly was made for convenience; otherwise, we could not exclude dependence on the impending earthquake magnitude (De Santis et al., 2019b). We first sorted the list of seismic events by time, selecting and removing the aftershocks for an earthquake (given earthquake) in the list in turn. After these operations, we dropped 3013 and 11,815 aftershocks for CSES and DEMETER data analysis. Finally, 2991 (CSES) and 8912 (DEMETER) independent earthquakes remained on the list, respectively. Additionally, we also deleted the data corresponding to the aftershocks.

To validate the dependability of machine learning technologies and to enhance their robustness, we generated the same number of artificial non-seismic events as real earthquakes, randomizing and shifting times and locations to avoid overlapping with the real earthquakes. We randomly sampled time, latitude, and longitude within the selected spatio-temporal range, following the given constraints: 1) the longitude or latitude is not within 10° before and after the longitude or latitude of a real earthquake, and 2) the time is not within 15 days before and after a real earthquake.

The operational modes of the Langmuir Probe onboard CSES include survey mode and burst mode (Yan et al., 2018a). As the sweeping period of the two modes is different, we used simple linear interpolation to interpolate the data in the survey mode so that the time resolution of the two modes was 1.5 s.

In our study, we combined a CNN and Bi-LSTM to train our proposed model SeqNetQuake (Figure 2). The model’s architecture is composed of one-dimensional convolutional layers, a one-dimensional Bi-LSTM structural layer, and a fully connected (FC) block. The SeqNetQuake model employs CNN layers to extract features from the input data and Bi-LSTMs for sequence prediction. SeqNetQuake reads subsequences of the main sequence as blocks, extracts features from each block, and then allows the LSTM layer to interpret the extract features. To enable the same CNN model to read each subsequence in the window, the whole CNN model is wrapped in a TimeDistributed layer. After flattening the retrieved features, they are forwarded to the Bi-LSTM layer for reading, and further features are extracted. Conceptually, the 1D convolutional layers are used to extract the data features, after which the Bi-LSTM structures optimize the feature extraction in the sequential data. Finally, an FC layer is employed to generate a classification probability. The loss function was categorical cross-entropy, and the optimization was performed using the Adam method (Kingma and Ba, 2014).

The suggested model was developed using the Keras (v 2.3.0) interface with TensorFlow 2.0 (Abadi et al., 2016). To facilitate fast training, all models were built on a server equipped with two Intel Xeon E5-2650 v4 CPU processors, 128 GB of RAM, and an NVIDIA GeForce RTX 2080 Ti graphics processing unit (GPU) (Oh and Jung, 2004). As the proposed method is sensitive to the selected parameters, Bayesian hyperparameter tuning (Snoek et al., 2012) was utilized to determine the parameters that provide the best performance, and the Hyperopt Python package (Bergstra et al., 2013) was used to implement it. In this process, the negative of the F-measure (F1) was used as the return value (loss) of the objective function. The process selects the most promising hyperparameters based on their ability to minimize an objective function by building a probability model based on the past evaluation results. Thus, this study can better perform with fewer iterations than the required efforts to conduct random and grid searches. Supplementary Table S2 provides the information on the search space for the significant parameters of SeqNetQuake. The maximum number of iterations for each model was set to 100. Supplemental Dataset S1 provides the hyperparametric optimization trial results for all the datasets passed through SeqNetQuake and other benchmarking classifiers when we train the model.

Transfer learning describes when a complex model that has been trained using a large dataset for a specific task is further trained for a related task with limited data (Pan and Yang, 2010). We used the training data from the large DEMETER dataset (TR0) for initial training, the test data of the large DEMETER dataset (TS0) was used for testing the initial model, and then transfer learning was employed using the small CSES training dataset (TR1) (Figure 2). All the weights/parameters that were learned using TR0 were retrained for TR1. We trained all the model weights without freezing any layer during transfer learning, as this approach provides better results (Hanson et al., 2020). The small CSES test dataset (TS1) was used as an independent test set during transfer learning, and the hyperparameters used were the same as those employed when we train the models on the TS0 dataset.

The DEMETER and CSES datasets are often class imbalanced where the number of the samples representing non-seismic class is much higher than the number of the samples in the other classes (Japkowicz and Stephen, 2002). In this case, a trivial classifier that predicts each sample as the majority class can obtain very high accuracy, and therefore, the overall classification accuracy is not accurate to evaluate the performance. Therefore, we use the F-measure (F1) evaluation model, which considers the correct classification of each class to be equally important. The F1-score is a metric that combines precision and recall. It is usually described as the harmonic mean of both. Thus, the class imbalance is countered by weighting different classes according to their sample proportions.

Beyond the metric mentioned above, which emphasize the positives, the Matthews correlation coefficient (MCC) (Matthews, 1975) is also adopted.

Further, receiver operating characteristic (ROC) curves, plots of the true positive rate (TPR) against the false-positive rate (FPR), were used to evaluate the classifier’s output quality in this study. ROC curves are typically used in binary classification contexts to evaluate the output of a classifier. To extend the ROC curve and ROC area for multiclass classification, the output is binarized, and one ROC curve can be drawn and used to evaluate classifier quality per class. In addition, we calculated the area under the ROC curve, termed as AUC, which is used to distinguish different models. Higher AUC values were considered to be indicative of superior methods for the identification of pre-earthquake ionospheric perturbations.

Finally, to visually demonstrate the classification performance of each class, ternary probability diagrams and confusion matrixes were used to show the probability distributions for each input class of the test data and the distribution of the predictions and actual values.

Five state-of-the-art methods were benchmarked for the study task: the gradient boosting machine (GBM) (Friedman, 2001), deep neural network (DNN) (LeCun et al., 2015), random forest (RF) (Geurts et al., 2006), CNN (Krizhevsky et al., 2012), and LSTM (Hochreiter and Schmidhuber, 1997) models. These methods were implemented in Python (v 3.6) with scikit-learn (v 0.20.0) and Keras (v 2.3.0). As the investigated models are sensitive to parameter selection, we chose parameters that yielded the best performance using Bayesian hyperparameter tuning, as described above. After optimal parameters were determined for each method, the performances of the different methods were compared.

We first directly trained the proposed SeqNetQuake model using the CSES dataset, which was further contiguously divided into a training set (TR1) and a test set (TS1). Given that there is no universal standard for the lengths of the input sequence and the spatial window, we initially configured the data with 10 consecutive observations as the input sequence length, a spatial window centered at the epicenter, a deviation of 3°, and nighttime data in the initial configuration (DataSet 01 in Supplementary Table S1).

As shown in Supplementary Figure S2, ROC curves are adopted as a performance metric, as they depict relative trade-offs between true positives (benefits) and false positives (costs) for each class; the model’s performance using nighttime data is displayed for three classes. Supplementary Figures S2A–C demonstrates that the AUC values of classes 0 and 1 are higher than 0.7, indicating that the model can roughly identify the time series related to earthquakes and non-seismic events, but the AUC of class 2 is only 0.6128, indicating that the model’s accuracy in identifying space weather such as magnetic storms is not high. This may be attributed to the smaller number of samples used to train class 2, which could have led the model to fail to extract the features of this class. Supplementary Figure S2D shows the bar plot curves of MCC, F1 score, and accuracy, which reflect the model’s overall performance. The results are similar to those implied by the ROC curves, indicating that the model can distinguish, to some extent, earthquakes, non-seismic, and space events. In general, the performance of the model based on the initial configuration was reasonable but relatively weak. Therefore, we explored whether combining datasets with different temporal and spatial features or different models may offer better performance.

A DEMETER dataset was generated based on the optimal spatio-temporal feature configuration for the CSES data (DataSet 11 in Supplementary Table S1). In this study, deep transfer learning, although it is first applied to DEMETER to be transferred to CSES, considers the properties and features of CSES data. In other words, it is shaped/adapted to CSES data features.

We first trained the SeqNetQuake model using the full DEMETER dataset, which was contiguously divided into a training set (TR0) and test set (TS0). Using TR0 for training, TS0 for testing, and the sliding sequences as model input, we trained deep-learning models with its architecture shown in the left panel of Figure 2. The performance of the best model using TS0 is shown in Table 1. The F1 score of 0.7002 and MCC of 0.3763 for TS0 suggest the robustness of the trained model.

The model obtained using the DEMETER data was further trained using the CSES data employing TR1 (training set) and TS1 (test set), which corresponds to a transfer-learning process. The TS1 set was independent of the training data (TR0 and TR1). Table 1 and Figure 6D further illustrate the performance of the initial training and transfer training. The F1 score improved by 8% from 0.7002 to 0.7607, and MCC increased from 0.3763 to 0.5523 using TR0 for training, confirming the robustness of the model trained using the larger DEMETER dataset. The ROC curves in Figures 6A–C provide a visual performance comparison of each class and further demonstrate the result.

We selected the model with the best performance in direct learning, the SeqNetQuake model trained on Dataset 07 (called SeqNetQuake-DT), and compared it with the model trained using transfer learning (called SeqNetQuake-TL). Table 2 compares the performance of the two approaches, and Figure 6D shows the bar plot of the performance metrics. Figures 6A–C compares the ROC curves yielded by direct training and transfer learning for the independent test set TS1.

The performance of direct learning was significantly lower than that of transfer learning, yielding a 10% reduction in the F1 score and a 23% reduction of the MCC value. These results confirm the difficulty of using small training data sets (TR1) for direct learning and the need to use large data sets that can effectively exploit the capabilities of the deep-learning networks.

To further confirm the performance of transfer learning, ternary probability diagrams and a confusion matrix were used to indicate the distribution of the predicted and true values and allow more profound insight into the classification performance of the models. Supplementary Figure S3 and Figure 7 show the ternary probability diagrams and the confusion matrix for the three classes obtained from SeqNetQuake-DT and SeqNetQuake-TL, respectively.

The ternary probability diagrams allow a qualitative evaluation of the classification results. The number of the correctly classified samples in class 0 increases significantly, as evidenced from Supplementary Figure S3A and Figure 7A; the same trend also appears in the qualitative comparison of class 1 in Supplementary Figure S3B and Figure 7B. For class 2, the direct training results (Supplementary Figure S3B) show that the probability that the model would discriminate the true class is generally only around 0.4. With transfer learning, the probability is significantly improved and is typically higher than 0.5 (Figure 7C). As class 2 only represents 10.5% of samples, far less than the proportions of the other two classes (which are 45.3 and 44.1%, respectively), the model likely fails to fully learn the features of this class, causing a higher error rate.

The confusion matrixes shown in Supplementary Figure S3D and Figure 7D show significant improvements in precision regarding class 2, from 75.3% before transfer learning to 95.2% after. Further, transfer learning achieves 9 and 13.06% improvements in precision for classes 0 and 1, respectively. Therefore the misclassification rate of each class drops significantly after transfer learning. Transfer learning achieves satisfactory results in the classification of classes 0 and 1, and the probability distributions for each class of the test data (Figure 7D) show that 77.9 and 85.1% of the input samples, respectively, are correctly classified (recall).

Table 2 and Figure 8D report the performance of our transfer-learning model (SeqNetQuake-TL) with five other benchmarking classifiers for the independent test set TS1-07. The performance of the existing methods ranges from F1 = 0.5757 to 0.7607 and MCC = 0.2365 to 0.552. However, SeqNetQuake-TL offers the best performance, improving MCC by 64.2% MCC and F1 by about 24.5% over those for the next-best CNN model. Figures 8A–C compares the ROC curves obtained for the SeqNetQuake-TL model with those of the five other classifiers, and SeqNetQuake-TL again demonstrates the best performance with a 14.5% improvement in AUC for class 0 and a 12.2% improvement for class 1 over the second-best CNN model and a 4.4% improvement in AUC for class 2 over the second-best LSTM model.

Excellent neural network architecture may explain why SeqNetQuake outperforms all the other predictors. First, SeqNetQuake’s CNN network layers allows us to extract local parallel features. Then, the Bi-LSTM layer (composed of multiple memory modules with a two-cell topological structure) extracts long-distance dependent features and performs sequence learning, thereby mitigating the influence of the front and back features of each attribute feature point. Finally, the model obtains the classification results through the FC layer and the Softmax classifier, improving SeqNetQuake’s accuracy and reducing false positives. Notably, the DNN, CNN, and LSTM models performed better than the RF and GBM models, suggesting that the deeper neural network models are much more efficient in computation and the number of the parameters.

Lithospheric-ionospheric coupling is a topic that has been qualitatively discussed in several papers. Several studies have investigated the physical mechanisms of ionospheric pre-earthquake perturbations (Hayakawa et al., 2010; Pulinets and Ouzounov, 2011; Wu et al., 2012; Ouzounov et al., 2018; De Santis et al., 2019a; De Santis et al., 2020; Freund et al., 2021). Hypotheses regarding these mechanisms were presented by Pulinets et al. (2015) and Kuo et al. (2014), who proposed complex lithosphere–atmosphere–ionosphere coupling as the physical basis of the generation of short-term earthquake precursors. Recently, a detailed mathematical model providing a quantitative description of the Magnetospheric-Ionospheric-Lithospheric-Coupling (MILC) was introduced and successfully tested with data from individual earthquakes (Piersanti et al., 2020a). A key feature of the MILC model is based on the development of an Acoustic Gravity Wave (AGW) (Carbone et al., 2021) interacting mechanically with the ionosphere and then electromagnetically with the magnetosphere. The AGW coupling mechanism is quite general and can provide for both co-seismic and precursor lithospheric-ionospheric coupling. In the case of co-seismic coupling, the AGW is generated by ground surface motion (solid, liquid). In contrast, in the case of precursor coupling, the AGW can be generated by phenomena which modifying the temperature or the electrical properties of the atmospheric column above the EQ preparation zone.

One example of such mechanisms is radon release following by its radioactive decay, emitting alpha particles of 5.2 MeV, a large value used in comparison to the ionization energy required for the dissociation of an atmospheric molecule (32 eV). One alpha particle is sufficient to generate 150,000 pairs of positive and negative ions, thereby creating an excess of positive airborne ions near the Earth’s surface, which can influence the ionosphere. Radon gas accumulates with time in the lower crust and upper mantle and can exhibit large-scale distribution in the rocks. During the preparation phase of an earthquake, these gas domains can become hydrostatically unstable and force their way upward through the lithosphere. The rapid lithospheric degassing followed by radon atoms radioactive decays would trigger several atmospheric processes near Earth’s surface, leading to changes in air conductivity and temperature and therefore changes to the near-ground atmospheric electric field (Ouzounov et al., 2018; Pulinets and Ouzounov, 2018).

Crustal strain measurements typically fail to detect any unusual changes before earthquakes (Soter, 1999), while radon outgassing has been observed before the earthquake (Gold and Soter, 1985; Riggio and Santulin, 2015; Fu et al., 2019), and they are considered as a possible EQ precursor.