Singular Spectrum Analysis of the Total Electron Content Changes Prior to M ≥ 6.0 Earthquakes in the Chinese Mainland During 1998–2013

Early studies have shown evidence of the seismo-ionospheric perturbations prior to large earthquakes. Due to dynamic complexity in the ionosphere, the identification of precursory ionospheric changes is quite challenging. In this study, we analyze the total electron content (TEC) in the global ionosphere map and investigate the TEC changes prior to M ≥ 6.0 earthquakes in the Chinese Mainland during 1998–2013 to identify possible seismo-ionospheric precursors. Singular spectrum analysis is applied to extract the trend and periodic variations including diurnal and semi-diurnal components, which are dominated by solar activities. The residual ΔTEC which is mainly composed of errors and possible perturbations induced by earthquakes and geomagnetic activities is further investigated, and the root-mean-square error is employed to detect anomalous changes. The F10.7 and Dst index is also used as criterion to rule out the anomalies when intense solar or geomagnetic activities occur. Our results are consistent with those of previous studies. It is confirmed that the negative anomalies are dominant 1–5 days before the earthquakes at the fixed point (35°N, 90°E) during 0600–1000 LT. The anomalies are more obvious near the epicenter area. The singular spectrum analysis method help to establish a more reliable variation background of TEC and thus may improve the identification of precursory ionospheric changes.

Due to the trend, long period and short period terms, such as annual, seasonal, and diurnal variations in the ionosphere, which might be much stronger than the seismo-ionospheric perturbations, the identification of precursory ionospheric changes is quite difficult. There are several methods to detect the ionospheric anomalies, including the median-interquartile range (IQR) method (Liu et al., 2006;Chen et al., 2015) and the average-standard deviation (STD) method (Kon et al., 2011). Both methods usually take previous some days as the reference background and compute the dynamic ranges of TEC variation. Recently, Guo et al. (2019) analyzed the TEC of two earthquake cases by singular spectrum analysis (SSA) and indicated the presence of obvious anomalous characteristics of the seismicionospheric coupling effect. The implementation of SSA could help to eliminate the influence of various periodical changes and might improve the detection of TEC anomalies. Chen et al. (2015) found the electron density decreased 1-5 days before M ≥ 6.0 earthquakes in China by using the previous 15-day as the reference background. To test possible relationship between ionospheric disturbances and earthquakes, we conduct similar analysis as Chen et al. (2015) and examine the changes in GIM-TEC for M ≥ 6.0 earthquakes in the Chinese mainland during 1998-2013. In this study, the SSA method is used to remove the trend and periodical components, particularly the 27-day (Pancheva and Lastovicka, 1989;Pancheva et al., 1991) and diurnal variations in TEC. The possible effect of solar and geomagnetic activity is further discussed, as the ionosphere is sensitive to solar and geomagnetic activities (Carter et al., 2013).

DATA TEC and Earthquake Data
GIM-TEC data are derived from GPS signals. There are seven ionosphere analysis centers (IAACs) all over the world, and each IAAC can calculate the global ionosphere map (GIM) with different computing strategies and the final GIM products can be obtained by weighting the GIM from each IAAC (Li et al., 2017). In this paper, we obtain TEC data which are generated on a daily basis from the Center for Orbit Determination in Europe (CODE), University of Berne, Switzerland. GIM-TEC covers ± 180°( longitude) and ± 87.5°(latitude) with a spatial resolution of 5.0°and 2.5°, respectively. The GIM-TEC data are released once a day with a 2-h resolution. We apply interpolation using a cubic spline function (de Boor, 1978) to obtain TEC data every 1 h. A list of 56 M ≥ 6.0 earthquakes in the Chinese Mainland during 1998-2013 are used in this study, the same as the earthquake catalog used by Chen et al. (2015), so that we can compare the results with previous studies. Figure 1 shows the locations of these earthquakes and the GIM-TEC used. Because it is not possible to know the GIM TEC over the "epicenter" before the earthquake occurs, therefore, for the practical application, we should examine GIM-TEC at a fixed location to analyze whether there are anomalies before the surrounding earthquake occurrence.

Geomagnetic Dst and Solar Activity F 10.7 Index
Previous study suggested that the ionosphere was sensitive to solar and geomagnetic activities. Thus, it is necessary to analyze the geomagnetic indices and solar activity when identifying the earthquake-related TEC anomalies. The Dst index with a time resolution of 1-h provided by the World Data Center for Geomagnetism, Kyoto are used in this study. The geomagnetic storms can be classified into weak (−50 nT < Dst ≤ −30 nT), moderate (−100 nT < Dst ≤ −50 nT), and intense (Dst ≤ −100 nT) storms (Gonzalez et al., 1994). The solar radiation flux F 10.7 index provided by Space Physics Data Facility (SPDF) of NASA are used and F 10.7 is usually <100 SFU when solar activity is quiet (Guo et al., 2019).

METHODS
We apply the SSA method to study and decompose the TEC time series and use Correlation value to analyze main components. SSA is a model-free approach to analyze the periodic oscillation of time series, which can be used to extract information from the original Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 677163 signal containing different periodic components (Keppenne and Ghil, 1992;Robert and Pascal, 1992). First, we obtain the TEC dataset for each earthquake 31 days before and after its occurrence day. The timespan of each dataset is 63 days and centers on the earthquake day. The SSA results will be different if using different data length. We use the length of 63 days to do SSA and use 60 days to analyze so that it is consistent with the previous results. Some researcher used the different data length to do analyze, depending on their own demand (Shi et al., 2020;Zhang et al., 2021). Next, we create a trajectory matrix with a window length of 20 days for each dataset. We choose 20 days because that the window length (L) of SSA should be 2 < L < N/2 (Golyandina et al., 2001) and we use N/ 3, where N is the length of the time series. Then we obtain the reconstructed time series after embedding, SVD, grouping, and diagonal averaging. For two reconstructed series X (1) and X (2) , the inner product is defined by (Golyandina et al., 2001).
To measure the degree of approximate separability between two series X (1) and X (2) , the Correlation value is defined as The reconstructed components with large correlation values reflect that they are correlated from each other and usually composed of trend and periodic components, and the reconstructed components with small correlation values reflect that they are well separated from each other and contain different period components (Elsner, 2002). When extracting the reconstructed components, we select the first few principal components according to their energy, sorted by the corresponding eigenvalues λ i . When the energy components is lower than 0.12%, we stop and take the rest components as the noises part. We select 0.12% as cut-off energy because this threshold value can ensure the TEC data of 56 earthquakes to exclude diurnal and semi-diurnal variations. Finally, we remove the semi-diurnal, diurnal, and long-term variations and obtain the residual denoted by ΔTEC. ΔTEC is mainly composed of computing errors and possible perturbations induced by earthquakes and geomagnetic activities. Note that here we use the TEC principal components for each earthquake 30 days before and after its occurrence day to avoid the boundary effects.

Extraction of Main Components by SSA
We start by looking at an example of TEC data before and after the 2008 Wenchuan earthquake (Ms8.0) to show the extraction of main components by SSA. Figure 2 is an example of the Wenchuan earthquake from April 12, 2008 to June 11, 2008. We extract the first seven components based on the cut-off energy. Figure 2A is the correlation result of the first 20 reconstruction components by using SSA. It can be found the first seven components are separable and the eighth to later components are mixed with each other. Figure 2B shows the first seven reconstructed components. The first component is the trend component. The second to the fifth components are the diurnal variation and semi-diurnal variation, respectively. The sixth to the seventh components are approximately the 27-day variation and its half cycle. The ionosphere is an integrated product of the interaction between the Sun and the Earth's environment. The TEC changes are mainly controlled by the intensities of solar electromagnetic radiation (He et al., 2012). Analyzing the correlation between the TEC component derived from SSA including trend and periodical components and geophysical indices including F 10.7 and Dst may help us to know the ionospheric effect from solar activity. The comparison of reconstructed components and F 10.7 /Dst are as shown in Figure 3, and their Pearson correlation analysis are calculated. Because the F 10.7 data are of the 1-day resolution, the mean GIM-TEC of each day are used. The daily geomagnetic disturbance is measured by the minimum of the Dst indices on the day. As shown in Figure 3A, it is evident that the variation of the reconstructed components (RCs) 1-7 is similar to that of the F 10.7 , and the RC 1 shows similar trend with F 10.7 , suggesting a high correlation. The Pearson correlation analysis is demonstrated in Figure 3B, the correlation coefficient between RCs 1-7 and F 10.7 reaches 0.8, and the correlation coefficient between RC 1 and F 10.7 is 0.63. By the correlation analysis, it could be found that the influence from solar radiation at the daily timescale and above can be efficiently removed through the SSA method. In contrast, in Figure 3C, the correlation coefficient between RCs 1-7 and Dst is −0.22, between RC 1 and Dst is −0.26, showing a weak negative correlation.
After using SSA to extract the trend and periodic components which are mainly from the effects of solar activities, we subtract them from TEC data and obtain the residual ΔTEC. To give a more panoramic view, Figure 4A shows the original and reconstructed TEC results of the Wenchuan earthquake. Figure 4B shows the ΔTEC of the Wenchuan earthquake and its upper/lower bounds with ±RMS. The RMS is the root-meansquare error provided by the CODE on each grid point, which reflects the precision of the TEC data. Similar to (Guo et al., 2019), we use the RMS as the criterion for TEC anomaly detection. The ΔTEC value exceeding the upper/lower bounds implies that the ionosphere might be affected by solar activities or earthquakes. The ΔTEC at the fixed point (35°N, 90°E) for Wenchuan earthquake agrees with the results of (Chen et al., 2015). The ΔTEC decreased 1-4 days (from May 8 to May 12, 2008) before the earthquake and increased from May 3 to 5 May, as shown in Figure 4B.

Statistical Results
We select the most significant time 0600-1000 LT and calculate the percentages of the earthquakes with negative and positive anomalies that appear 30 days before and after the 56 earthquakes Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 677163 (following Chen et al., 2015). If one-third of the time in the 0600-1000 LT time span (3 h or more) occurred negative or positive anomalies, we count it a negative or positive anomaly day (following Chen et al., 2015). Figure 5 displays the percentages of the earthquakes with negative and positive precursory days and gives which kind of precursor is dominated. There are obvious negative anomalies 1-5 days before earthquakes, especially for M ≥ 6.8 and M ≥ 7.0. In general, the observed negative or positive precursors are consistent with the result of Chen et al. (2015). It should be noticed there are also higher positive abnormal proportion on 7-19 days before earthquakes, implying that the observed TEC tend to be relatively larger during this period. Therefore, if using the previous 15 days as the background, the upper and lower Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 677163 4 bounds of the dynamic ranges for TEC anomaly detection will become high. This might lead to more negative (below the lower bound) and less positive (above the upper bound) in the subsequent days. The proportion of anomalies overall in (Chen et al., 2015) is higher than that in our study, implying that the TEC anomalies reduce after applying SSA. This is consistent with the results of (He et al., 2012) that the ranges of TEC anomalies decreased on the global scale after removing background variations produced by solar radiations.
Even if the results in Figure 5 can be deemed to remove the influence of solar activity after SSA method, we redefine an anomaly hour when the TEC exceeds the upper or lower bound and both geomagnetic and solar activity are quiet (Dst index is more than −30 nT and F 10.7 is <100 SFU). If 3 h or more in the 0600-1000 LT occurred negative or positive anomalies, it is counted as a negative or positive anomaly day. Figure 6 shows results after ruling out the anomalies when geomagnetic or solar activities are not quiet. Negative anomalies are still dominant in 1-5 days before earthquakes, but the proportion decreased compared to that in Figure 5.

DISCUSSION AND CONCLUSIONS
A number of studies have shown the ionospheric TEC changes prior to great/mega earthquakes worldwide (e.g., Heki, 2011;Guo et al., 2015;He and Heki, 2016;Iwata and Umeno, 2016;He and Heki, 2017;Liu et al., 2018). As suggested by He et al. (2012), one of the key points is the elimination of ionospheric effect from solar radiation. In this study, we demonstrate another possible approach of SSA to eliminate such solar effects. By using SSA, the GIM-TEC anomalies prior to 56 M ≥ 6.0 earthquakes in the Chinese mainland during 1998-2013 have been revealed. Our findings agree with previous studies that negative anomalies are dominant during 1-5 days before earthquakes at the fixed point (Chen et al., 2015). The results persist after removing the data during intense geomagnetic and solar activity. As shown in both Figures 5 and 6, overall, the proportion of precursory days increases with the rise of earthquake magnitude, implying that the GIM-TEC anomalies are more likely to cooperate with stronger events.
Due to the ionospheric dynamic complexity, the earthquakerelated perturbations might be far less than the seasonal and diurnal variations in the TEC data. Most studies took the previous 15/30 days as a reference and built the bounds for the TEC variation on the target day. As shown in Figure 2B, there are long-term trends in the TEC data. It may overestimate the background using the previous 15-day data if the TEC is on a downswing. On the contrary, the background may be underestimated if the TEC is in an uptrend. These may cause bias in detecting anomalies in the TEC data. As an attempt, we apply SSA to extract long-term, diurnal, and semi-diurnal variations. The residual ΔTEC contains errors, and possible perturbations induced by earthquakes and geomagnetic activities. By using the RMS as a threshold, we could remove the influences of computing errors to a certain degree, as the RMS gives the precision of the TEC data. The approach based on SSA and RMS have an advantage in identifying anomalous changes in TEC variation (Guo et al., 2019). It should be noticed that some earthquakes' epicenters may be far from the point at the GIM map that we used to analyze preseismic TEC variations, especially for the M6-M6.5 earthquakes, some of them are outside the Dobrovolsky earthquake preparation zone, leading to unobvious ionosphere reaction on the GIM selected point. Thus, we also take into account the preearthquake anomalies over the epicenter of the 56 earthquakes. Figure 7 shows the ΔTEC of the point which is closest to epicenter for the Wenchuan earthquake and its upper/lower bounds with ±RMS, here we use the same method above to extract the main components. There are more obvious positive and negative anomalies compared to Figure 4B, especially positive abnormal enhancement on May 9, which are in line with the results of Liu et al. (2009).
In our statistical results, we focus on the time of 0600-1000 LT, if one-third of the time in the 0600-1000 LT (3 h or more) occurred negative or positive anomalies, we count it a negative or positive anomaly day (following Chen . For 1-5 days before M ≥ 6.0 earthquakes, the highest proportion of both positive and negative anomalies is about 23%, as shown in Figure 5. After ruling out the anomalies when intense geomagnetic or solar activities occur, the highest proportion of negative and positive anomalies are 19 and 13%, respectively, as seen in Figure 6. The proportion decreases compared to that in Figure 5 but negative anomalies are still dominant on 1-5 days before earthquakes. Figure 8 displays the percentages of the earthquakes with negative and positive precursory days over epicenter and gives which kind of precursor is dominated. In Figure 8, when using the TEC data close to the epicenter, the highest proportion of negative and positive anomalies increase to 28 and 34% on 1-5 days before M ≥ 6.0 earthquakes. The percentage is higher than those based on the TEC data in a fixed area in Figure 5, and the proportion increases with the magnitude. The proportions reach up to 100% (negative anomalies) and 33% (positive anomalies) for M ≥ 7.0 earthquakes. The results display that the anomalous variations are more significantly enhanced near the epicenter area. After ruling out the anomalies when intense geomagnetic or solar activities occur, the highest proportion of negative and positive anomalies are 15 and 18% respectively on 1-5 days before M ≥ 6.0 earthquakes, as shown in Figure 9, the proportion decreases compared to that in Figure 8 but negative anomalies are still dominant on 5 days before. The positive anomalies become obvious on 2-4 days before M ≥ 7.0 earthquakes, showing agreement with the positive abnormal enhancement on May 9, 3 days before the Wenchuan earthquake (Liu et al., 2009).
The characteristics of pre-earthquake TEC changes were quite different in different areas. For example, the negative anomalies are found to be dominant 1-6 days before earthquakes in China (Liu et al., 2006;Chen et al., 2015); while the positive anomalies were considered prominent 1-5 days before earthquakes in Japan (Kon et al., 2011). If the earthquake samples are analyzed together, one may not find clear results as the correlation between either negative or positive anomalies and earthquakes will decrease. This might be the reason why no significant anomalies for global earthquakes were detected in some studies (Zhu and Jiang, 2020). In addition, the RMS of GIM-TEC data is different at different places due to the uneven distribution of GPS stations. The level of uncertainty is high if using the GIM-TEC with large RMS, so the data accuracy should be also considered when doing global statistical analysis.
As for the physical mechanism of ionospheric anomalies associated earthquake, there are several ideas have been proposed by some scholars in the previous studies. The internal gravity waves (IGWs) mechanism has been provided by Pertsev and Shalimov (1996), which origin from seismogenic zone with period of 1-3 h due to inflowing of lithosphere gases into the atmosphere before earthquake. Liu et al. (2016) indicated that the triggered acoustic and/or gravity waves in the atmosphere near the surface can be activated by vertical ground motions, which is considered coseismic phenomenon of Tohoku earthquake . Some scholars described the formation mechanism of earthquake ionospheric precursors by seismogenic electric field with amplitude (Sorokin and Chmyrev, 1999;Pulinets and Boyarchuk, 2004), which charged aerosols FIGURE 5 | Percentages of the earthquakes with negative (black dot) and positive (gray dot) precursory days that appear 30 days before and after the earthquakes in China during 1998-2013 in the 0600-1000 LT. The black bar represents the amount of percentage in which negative anomaly is over positive anomaly, while the gray bar denotes the amount of percentage in which positive anomaly is over negative anomaly.
Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 677163 FIGURE 6 | Percentages of the earthquakes with negative (black dot) and positive (gray dot) precursory days after ruling out the anomalies when geomagnetic or solar activities are not quiet. The black bar represents the amount of percentage in which negative anomaly is over positive anomaly, while the gray bar denotes the amount of percentage in which positive anomaly is over negative anomaly. Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 677163 8 FIGURE 8 | Percentages of the earthquakes with negative (black dot) and positive (gray dot) precursory days over epicenter that appear 30 days before and after the earthquakes in China during 1998-2013 in the 0600-1000 LT. The black bar represents the amount of percentage in which negative anomaly is over positive anomaly, while the gray bar denotes the amount of percentage in which positive anomaly is over negative anomaly.
FIGURE 9 | Percentages of the earthquakes with negative (black dot) and positive (gray dot) precursory days over epicenter after ruling out the anomalies when geomagnetic or solar activities are not quiet. The black bar represents the amount of percentage in which negative anomaly is over positive anomaly, while the gray bar denotes the amount of percentage in which positive anomaly is over negative anomaly.
Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 677163 9 injected in the atmosphere and air ionization by radon followed by formation of large ion clusters of aerosol size during preparation. Namgaladze et al. (2009) proposed the most probable formation was the vertical transport of F2-region ionosphere plasma under action of zonal electric field (Namgaladze et al., 2009), some phenomena including geomagnetic conjugacy of the ionosphere precursors are the strong arguments in favor of this hypothesis. One more possible formation mechanism about the abnormal electromagnetic fields and emissions has been offered (Hayakawa and Molchanov, 2002), which originated from the ground due to certain mechanisms such as piezo-electric effects (Finkelstein et al., 1973), tribo-electric effect (Gershenzon et al., 1993) and positive holes (P-holes) (Freund, 2000), although this mechanism is found to be insufficient because of the weak intensity of lithosphere radio emissions (Hayakawa, 2007). Among the possible mechanisms, air ionization by natural ground radioactivity is considered more credible by some scholars (Pulinets, 2012;Guo et al., 2015), which is conform to the ionospheric characteristics of Chile and Wenchuan earthquakes (Pulinets and Davidenko, 2014).
In summary, by applying the SSA method, we extract the trend and periodical background components of GIM-TEC to analyze the pre-earthquake ionospheric anomalies associated with M ≥ 6.0 earthquakes in the Chinese Mainland during 1998-2013. The results are consistent with those of Chen et al. (2015) in general. It is confirmed that the negative anomalies are dominant 1-5 days before the earthquake at the fixed point (35°N, 90°E) during 0600-1000 LT. The anomalies are more obvious near the epicenter area. As there are considerable errors in the GIM-TEC data, it is worthwhile to analyze highprecision satellite observation data to obtain more accurate results.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: "The TEC data are from the Center for Orbit Determination in Europe (ftp://ftp.unibe.ch/aiub/CODE), the solar radiation flux F 10.7 index are from Space Physics Data Facility (SPDF) of NASA (https://cdaweb.gsfc.nasa.gov/index. html/) and the Dst data are from the World Data Center for Geomagnetism (http://wdc.kugi.kyoto-u.ac.jp/)."

AUTHOR CONTRIBUTIONS
Conceptualization, PH and HC; methodology, HC and PH; software, HC; validation, HC, MM, YC, QW, XS, KH, and PH; formal analysis, HC; investigation, HC, QW, XS, KH, and PH; resources, HC, MM, and YC; data curation, QW, XS, KH, and PH; writing-original draft preparation, HC; writing-review and editing, XS, KH, and PH; visualization, HC; supervision, PH; All authors have read and agreed to the published version of the manuscript.