Precursor Analysis Associated With the Ecuador Earthquake Using Swarm A and C Satellite Magnetic Data Based on PCA

A large magnitude 7.8 earthquakes occurred in Ecuador at 23:58:36 Coordinated Universal Time (UTC) on April 16, 2016. In this paper, we revisit this earthquake by simultaneously analyzing the magnetic field data of Swarm satellite A and satellite C based on principal component analysis (PCA), and the eigenvalues and principal components are calculated throughout 2016. We find that the first principal component mainly contains the signal originating from solar-terrestrial effects such as geomagnetic activity since the first eigenvalue and the geomagnetic index are highly correlated. Therefore, the second principal component is used to extract the anomalies associated with the Ecuador earthquake in terms of skewness and kurtosis. The anomalous tracks of the S-K (Skewness-Kurtosis) coefficient are accumulated from 90 days before the event to 30 days after. The cumulative number follows an accelerating power-law behavior before the earthquake and decelerating recovery behavior after the earthquake; moreover, the inflection point of the sigmoidal fitting curve is close to the time of the earthquake. The cumulative number of anomalous tracks in the random regions shows a linear increase, further verifying the correlation between anomalies extracted and the Ecuador earthquake. This phenomenon could be related to the preparation phase and the aftershock phase of the Ecuador earthquake.


I. INTRODUCTION
Short-term earthquake prediction is one of the most challenging issues in the world. Electromagnetic disturbances associated with earthquakes were discussed a few decades ago [1]- [3] and are probably the most promising candidates as possible precursor signals. As ground-based observations are limited by the study area and events, satellite electromagnetic observations have been developed since the 1980s, and ionospheric electromagnetic precursors based on satellites have been investigated [4]- [8]. Electromagnetic perturbations related to earthquakes have a fairly wide frequency range, from direct current to very high frequency [9]- [11]. Since the ULF (ultralow frequency) electromagnetic waves may penetrate through the lithosphere to the surface and propagate into the upper ionosphere and magnetosphere with small attenuation [12]- [14], the ULF (f<10 Hz) disturbance The associate editor coordinating the review of this manuscript and approving it for publication was Nilanjan Dey. is considered one of the most promising earthquake precursors [15]- [17]. This potential has also been verified by the ULF-band electromagnetic observations of DEMETER (Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions) [18]- [20].
The Swarm satellites measure the Earth's magnetic field with unprecedented precision, and their global high spatial coverage and short revisit time provide powerful data and technical support for the accurate separation and modeling of the magnetic fields in each layer. To study the core magnetic field and its long-term changes, various geomagnetic field models have been investigated, such as the IGRF (International Geomagnetic Reference Field) model, the CHAOS model, and the CM (Comprehensive Model), among others [21]- [23]. The lithospheric magnetic field is of great significance for geophysical explorations and tectonic interpretations. Lithosphere magnetic field inversions based on Swarm magnetic fields have been carried out [24], [25]. CI (comprehensive inversion) [26] and VOLUME 7, 2019 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ SI (sequence inversion) [27] are two main methods for lithospheric field inversion modeling. Precision orbit determination is being studied to monitor Earth's gravity field with high precision via the on-board GPS (Global Positioning System) tracking of Swarm [28], [29]. Field-aligned currents have also been studied to understand the solar wind-magnetosphereionosphere coupling mechanism [30]- [32].
In addition, based on the high-precision magnetic field of the Swarm satellites, ionospheric anomalies that occur before earthquakes have been studied. Studies on the magnetic field in the ULF-band recorded by the Swarm satellites before the 2015 Nepal earthquake further demonstrate the possible coupling between the lithosphere and the ionosphere [33]. As there are three satellites in the Swarm constellation operating at different orbits and altitudes, anomalies observed in each satellite can be checked against the data from the others. For examples, the magnetic field anomalies observed in the Swarm A satellite also appeared in the Swarm C satellite in an analysis of parameters around the 2016 Ecuador earthquake due to the short distance between the two satellites [34]. Moreover, the electron density differences between Swarm A and C were calculated to detect anomalies since these two satellites fly side by side at almost the same altitude [35].
Multi-satellite monitoring is similar to multi-station monitoring on the ground. For ground-based stations, PCA (principal component analysis) has been applied to ULF magnetic field data from three stations to separate signals generated by different sources and extract the main anomalous features of an earthquake [36]. PCA has also been applied to 3-component borehole strain data to extract the anomalies in the eigenvalues and eigenvectors [37].
Since the earthquake in Ecuador in 2016, many studies have focused on anomalous analysis before the earthquake. Based on VTEC data, He and Heki (2017) identified anomalous signal almost immediately (∼17 min) before the onset of the earthquake through the use of the reference curve method, and Carlos Sotomayor-Beltran (2019) observed a negative ionospheric anomaly 10 days before the incident by using a statistical method [38], [39]. An unusual increase in the environmental radiation level a few hours before the earthquake highlights the importance of the precursor role of environmental radiation in the precise location of earthquakes in Ecuador [40], [41]. Twenty-four GPS stations distributed throughout Ecuador recorded significant changes in their usual locations a few minutes before the earthquake, which can be used for forecasting in the medium and short terms along active continental margins [41].
In this paper, we explore a possible precursor before 2016 Ecuador earthquake by analyzing Swarm magnetic field data. First, we briefly introduced the Swarm satellite data and the method of application. In the following, we regard Swarm A and C satellites as two moving stations and apply PCA to the magnetic field data. Then, we calculate an S-K coefficient of the second PC (principal component) based on skewness and kurtosis to extract the anomalies associated with the 2016 M7.8 Ecuador earthquake, and the accumulation of anomalous tracks is calculated over time to study the preparation phase and aftershock phase of the earthquake. Finally, we determine the correlation between the extraction anomalies and the earthquake by comparing the cumulative number of anomalous tracks in the random regions with that in the actual seismic region.

A. SWARM A AND C SATELLITE DATA
Swarm, a satellite mission of the ESA (European Space Agency), was launched into a near-circular polar orbit on November 22, 2013, to accurately measure the different magnetic signals that originate in the Earth's core, mantle, crust, oceans, ionosphere and magnetosphere; respectively, the signals from these sources form the Earth's magnetic field. The Swarm satellite mission consists of three satellites, namely, A (Alpha), B (Bravo), and C (Charlie). The Swarm A and C satellites fly almost side by side at a constant speed of ∼7.5 km/s in low orbit at an altitude of approximately 460 km and an inclination of 87.4 • ; in the east-west direction, they have a longitudinal separation of 1.4 • at the Equator, and in the north-south direction, they have a maximum differential delay of approximately 10 s. Swarm B flies at a higher orbit with an altitude of 510 km and inclination of 88 • . Swarm B drifts longitudinally from Swarm A and C by 20 • every year. The three satellites complete one orbit in ∼90 min and hence accomplish nearly 15 day and night passes every 24 h. Because they drift slowly by ∼1h in LT (local time) every 11 days, Swarm satellites A and C need approximately 133 days to cover all 24 LTs and Swarm B needs approximately 141 days for covering all 24 LTs [42]- [43]. The Swarm constellation and its orbital properties enable the accurate measurement of geomagnetic field signals, especially those linked to the lithosphere.
Each satellite is equipped with the following set of identical instruments: an ASM (absolute scalar magnetometer), a VFM (vector field magnetometer), an STR (star tracker), an EFI (electric field instrument), a GPSR (GPS receiver), an LRR (laser retro-reflector) and an ACC (accelerometer). Among them, there are two magnetic field measuring instruments: the ASM, which measures the magnetic field intensity and provides a scalar magnetic field to calibrate the VFM, and the VFM, which provides a high-precision vector magnetic field.
Since Swarm A and C are very close together and fly at almost the same height, the magnetic fields measured by these two satellites at the same time have close similarities. In this study, we employ Level 1b vector magnetic field data with a 1-Hz sampling rate in the NEC (north east center) frame recorded by the VFM of the Swarm A and C satellites (https://swarm-diss.eo.esa.int).

B. PRINCIPAL COMPONENT ANALYSIS
PCA [44][45][46][47] is a multivariate statistical method that projects the original time series (from strong components to weak components) onto a new set of spatial orthogonal bases, thereby transforming multiple original indicators into a few comprehensive indicators that are independent of each other. The goals of PCA are to extract the most important information form the data table, compress the size of the data set by keeping only this important information to simplify the description of the data set and analyze the structure of the observations and variables. The stronger the correlation in the original data is, the better the feature extraction effect is. In this way, the problem of identifying relatively weak signals against strong background interference can be solved.
The original data X are expressed as where m is the number of samples, and n is the dimension of the data. First, we calculate the covariance matrix R X (n×n), and the elements r uv can be calculated as follows: where x iu and x iv are the uth and vth columns, respectively, of the ith row of data, andX u andX v are the averages of the uth and vth columns of data, respectively.
We next perform eigenvalue decomposition on the covariance matrix R X : where is the eigenvalue matrix (λ 1 , λ 2 , . . . , λ n ), and V is the corresponding eigenvector matrix, the columns of which are v 1 , v 2 , . . . , v n . Through PCA, we project the observed signals onto projection axes that are orthogonal to each other. The projection variance of X on v 1 is the largest (occupies the most energy) and is thus called the first PC. The eigenvector v 1 is considered to be the most dominant or intense signal subspace, and the corresponding eigenvalue λ 1 reflects the amount of energy in this subspace. The projections of the observed signals X on v 2 , . . . , v n are called the second PC, third PC and so on, and λ 2 , . . . , λ n correspondingly reflect the energy occupied by the second PC, third PC and so on. Each obtained eigenvector indicates the basis function of the corresponding signal subspace, and its eigenvalue reveals the power within that subspace.
The principal components are calculated as follows: where contains all the information of the original data and is arranged according to the variance from largest to smallest. An n-dimensional signal can be divided into n unrelated principal components at maximum.

A. PCA OF THE SWARM A AND C MAGNETIC FIELD DATA
In this part, we apply the PCA to the magnetic field data (1 Hz) measured in 2016 by the VFM instruments onboard the Swarm A and C satellites at geomagnetic latitudes between 50 • and −50 • to avoid polar disturbances. We consider the Y-component of the magnetic field (the east component in the NEC frame) since the clear anomalies generally occur in the Y-component [33], [34]. Since the ionospheric variation caused by the sun during the daytime may overwhelm the small disturbance caused by an earthquake, only satellite passes during the local nighttime (18:00-06:00 LT) are considered [48], [49]. Then the main IGRF [21] model is subtracted from the magnetic field measured by Swarm A and C satellites. The magnetic field data after subtracting the IGRF model of a Swarm A satellite track whose geomagnetic latitude is ±50 • are represented as follows: Similarly, the data of Swarm C satellite are represented as follows: where m indicates the length of tracks. The magnetic field of the A satellite is the same as the magnetic field of the C satellite at the UTC time. The data of Swarm A and Swarm C are stored as columns as follows: ].
Finally, the PCA is performed on the data matrix X = [B T A , B T C ] and the eigenvalues (λ 1 λ 2 (λ 1 > λ 2 )) and PCs are calculated through (3) and (4), respectively, when n is 2.
The eigenvalues λ 1 and λ 2 over 2016 are shown in Fig. 1, where λ 1 and λ 2 correspond to the first and second PCs, respectively. These eigenvalues indicate the energy of the two PCs. VOLUME 7, 2019 Through the PCA, we find that the amplitude of the first eigenvalue is much larger than that of the second eigenvalue. The contribution rate of the first PC is above 90%, which indicates that the first PC occupies the main part of the signal.
Swarm mainly measures the geomagnetic field [50]. The geomagnetic field is composed mainly of the internal magnetic field and the external magnetic field, of which the former occupies approximately 99% of the geomagnetic field. The internal source field is composed mainly of the core field (main magnetic field, ∼95%) and the crustal magnetic field (lithospheric magnetic field, ∼4%), which is stable. In contrast, the external source field is composed mainly of the ionospheric magnetic field, the magnetospheric magnetic field and the induced magnetic field, which is complex and varies with time, primarily due to diurnal variation, magnetic storms, solar rotation cycles and seasonal changes [51]. Earthquakes occur in the lithosphere; accordingly, earthquakes will lead to changes in the thermal state and stress state of the lithosphere, resulting in changes in its conductivity and magnetic susceptibility. The difference in the lithospheric magnetic field before and after an earthquake is enormously helpful for extracting seismic-related information. Therefore, the key issue in the extraction of seismo-ionospheric anomalies is how to separate seismic-related signals from other signals.
We have removed the main magnetic field by subtracting the IGRF model. The external source field is a changing magnetic field, which is divided into a calming magnetic field and a disturbing magnetic field. The latter originates mainly from various short-lived current systems formed by solar winds in the near-Earth space and is affected primarily by magnetic storms and geomagnetic pulsations. The spatial characteristics of the external source are therefore global in nature. Since the Swarm A and C satellites are very close together at almost the same altitude, the effects of global-scale interferences such as magnetic storms and geomagnetic pulsations on these two satellites should be similar. However, the time and distance required for signals associated with earthquakes to propagate to the magnetic sensors of the two satellites vary, and thus, the abnormalities recorded by the two satellites differ.
PCA projects a signal into a new coordinate system, and the original signal is divided into several unrelated principal components according to the variance from largest to smallest. Therefore, through PCA, the signal with the strongest correlation and significant features between the magnetic field recorded the by Swarm A and C satellites is separated into the first PC, and the weaker signals that are less correlated in the original signal are projected into the remaining principal components. In order to distinguish seismo-ionospheric perturbations from geomagnetic disturbances, geomagnetic indices should be considered [52], [53]. In Fig. 2, we illustrate the variations in the eigenvalues and in the corresponding geomagnetic ap index and Dst index. Fig. 2(a), Fig. 2(c) and Fig. 2(d) show that when the first eigenvalue increases, the corresponding geomagnetic index also increases. The variation in λ 1 seems to be adequately correlated with the geomagnetic index, and there are a number of simultaneous peaks. In contrast, the variation in the second eigenvalue λ 2 is not significantly similar to the geomagnetic index.
To quantitatively describe the correlations between the principal components and the geomagnetic activity, we use the cross-correlation function to calculate the correlation coefficients between eigenvalues λ 1 and λ 2 and the ap index. The formula for calculating the correlation coefficients of the two sequences x(t) and y(t) is as follows: where σ x and σ y are the standard deviations of x(t) and y(t), respectively; x and y are the mean values of x(t) and y(t), respectively; C xy (k) is the covariance of the two time series under the time delay k; r xy (k) is the cross-correlation coefficient of the two time series under the time delay k; and n is the length of the time series. The crosscorrelation coefficients between eigenvalues λ 1 and λ 2 and the ap index under a time delay of 0 are calculated as follows: where p denotes the number of tracks. We calculated the correlation coefficient between the first eigenvalue and the ap index in a sliding window. The result is shown in Fig.3.
As shown in Fig. 3, the correlation coefficient is mostly concentrated above 0.8-the median value and mean value of the correlation coefficient are 0.90 and 0.84, respectively, with confidence interval above 95%. These findings reflect the strong correlation between the geomagnetic index and the first eigenvalue and it is sufficient to show that PC1 mainly contains the influence of geomagnetic activity. In contrast, the correlation coefficient between the second eigenvalue and the ap index is approximately 0.2, which is very low. Therefore, the correlation of PC2 with the geomagnetic activity is not as obvious and PC2 can be seen to be hardly affected by geomagnetic activity.
Since the correlation between PC1 and the geomagnetic index is high, the data affected most by geomagnetic activity are separated mainly into PC1 through PCA. Next, to extract seismo-ionospheric anomalies associated with the earthquake, we will analyze PC2.

B. DETECTION OF ANOMALIES ASSOCIATED WITH THE 2016 ECUADOR EARTHQUAKE USING SKEWNESS AND KURTOSIS
A large earthquake with a magnitude of 7.8 occurred on April 16, 2016, at 23:58:36 UTC; the epicenter was near the coast of Ecuador, and the depth was 20.6 km. The earthquake ruptured the Ecuador-Colombia subduction boundary between the Nazca and South American plates, which is the area of the crustal deformation caused by the subduction of the Carnegie Ridge. The border between the Nazca and South American plates ranges from the southern coast of Chile to the Panamanian fault zone, covering almost the entire west coast of South America, with a total length of approximately 7,000 km, where several historical megathrust earthquakes have occurred [54], [55]. The earthquake tectonics and plates of 2016 Mw7.8 Ecuador are shown in Fig. 4. Dobrovolsky et al. (1979) showed that the area affected by an earthquake could be estimated by the following equation: R = 10 0.43×M , where R is the radius in km of the  affected area, and M is the magnitude of the earthquake [56]. To extract seismo-ionospheric anomalies associated with the Ecuador earthquake, we calculate PC2 of the Y-component magnetic field data acquired along tracks in a circular region defined by R = 2259.4km from the Ecuador earthquake epicenter from January 17 to May 16, 2016 (90 days before the earthquake and 30 days after the earthquake). The most significant earthquakes that occurred in this time interval are shown in Table 1.
Taking April 13-15, 2016, as an example, the epicenter, the Dobrovolsky circular and the tracks of the Swarm A and C satellites flying over the epicenter are shown in Fig.5.
Probability density functions (PDFs) have been suggested to be useful for understanding the phenomena. PDFs indicate whether a given phenomenon has a random character following a Gaussian distribution or whether it is intermittent and asymmetric, suggesting turbulence within the system [57]. The most useful parameters for evaluating turbulence are skewness and kurtosis. Skewness and kurtosis have been used for seismo-ionospheric anomaly extraction based on data recorded by DEMETER [57], [58].
Skewness, which describes the symmetry of a distribution, is the third central moment of a measured physical value normalized by the variance: where x represents the measured values and X denotes their mean value. Positive skewness indicates that the PDF has a longer tail for x − X > 0 than for x − X < 0. Hence, positive VOLUME 7, 2019  skewness means that the variable x is more likely to take large positive values than negative values. Kurtosis, which describes the thickness of the tail of a distribution, is defined as the fourth central moment of a measured physical value normalized by the variance: where x represents the measured values and X denotes their mean value. A time series in which most values are clustered around the average has low kurtosis, while a time series dominated by intermittent extreme events has high kurtosis. First, we calculate the first difference of the PC2 along the satellites' tracks. The results are very sensitive to any abrupt change by the first derivative. Then, the skewness and kurtosis of each track are calculated, as shown in Fig. 6.
It can be seen that the skewness and kurtosis exhibit little variation for a long time before the earthquake while closer to the earthquake, the amplitudes of the skewness and kurtosis are much higher than usual. After the earthquake, the skewness and kurtosis amplitudes are still slightly larger due to the influence of aftershocks, but they eventually return to normal levels. Thus, to evaluate the background value, we calculate the skewness and kurtosis of the PC2 in the study region for the two years from 2015 to 2016. We find that the skewness and kurtosis values are very low in most cases, with means of 0.4574 and 10.6719, respectively. Therefore, we define the S-K coefficient γ as in (12): where S is skewness, K is kurtosis, i is the ith number and N is the length of data. The S-K coefficient γ is sensitive to changes in the skewness and kurtosis, and thus can be used to extract anomalies in the skewness and kurtosis values. The curve of the S-K coefficient is shown in Fig. 7. At the beginning of 2016, the amplitude of the S-K coefficient γ is small on most tracks. Upon nearing the earthquake, the S-K coefficient γ exceeds the threshold (µ + 2 × σ ) on approximately 10 tracks and then recovers after a period of time following the earthquake. Therefore, this abnormality that changes with the earthquake may be related to the Ecuador earthquake. In addition, since we choose the tracks of the local night-time for analysis, 90% of the tracks are descending (from north to south), and it is very interesting that all abnormal tracks are descending.
Usually, for a random process, the accumulation of anomalies will show a linear increase. However, if a critical phenomenon occurs, anomalies will appear more frequently approaching the critical point, and the frequency will decrease subsequently. This behavior manifests as an accelerating power law before the event and as recovery afterward.
To study the preparation phase and aftershock phase of the earthquake, we calculate the accumulation of anomalous tracks over time. When the S-K coefficient γ of the track is greater than the threshold, we regard the track as an anomalous track. Moreover, if there are m anomaly tracks on the tth day, the cumulative sequence N(t) is increased by m on the previous day. The cumulative number of anomalous tracks over time is shown in Fig. 8. Fig. 8 clearly illustrates that before the earthquake, the number of anomalous tracks slowly increased. Approaching the earthquake, the number of anomalous tracks rapidly increased and then recovered to slow growth after the earthquake. We use a sigmoid function to fit the cumulative number of anomalous tracks. The sigmoid function is expressed as follows: where A1, A2, x0, and dx are calculated by fitting and x0 is the inflection point of the function. From the fitting results shown in Fig. 8, we can see that there is a good fit between the sigmoid function and the cumulative number of anomalous tracks. The variance is very small close to 0, and the correlation coefficient is close to 1. The fitting curve is concave upward before the earthquake and concave downward after the earthquake. As the earthquake approaches, the slope of the curve increases, reaching its maximum near the earthquake, and then slowly decreasing. Also, it is interesting that the value of x0 obtained in the fitting result is very close to the time of the Ecuador earthquake, here fixed as 0.
Moreover, this result is very similar to the results of a precursory study on an earthquake in Nepal [33]. The anomaly accumulation of the Swarm A satellite magnetic field data one month before and after the Nepal earthquake without consideration of night or day, quiet or disturbed times shows that the cumulative number follows the power law, showing an accelerated growth before the earthquake and decelerating recovery after the earthquake. In addition, the inflection point of the sigmoid fitting curve is close to the time of the earthquake.

IV. CONFUTATION ANALYSIS
To study the correlation between extracted anomalies and earthquakes, Parrot proposed a method for random earthquake distribution in 2011 [59]. The location of the earthquake epicenter is randomly changed while the size of the study area and the study time range are kept unchanged. The anomalous variations are compared between the random region and the actual seismic region, and the correlation between the anomalies and the earthquake is determined. In order to verify the relationship between the PC2 anomalies extracted and the Ecuador earthquake, we did a random earthquake distribution study. The distribution of earthquakes above Mw 6 worldwide from January 17 to May 16, 2016 is shown in Fig. 9.
We selected three locations worldwide as the epicenters of random earthquakes. In Fig. 9, the red, pink, and green asterisks are the epicenter positions of the three random earthquakes. The circles are random regions. To avoid anomalies being affected by other earthquakes, there are no earthquakes in the three random regions. To be more representative, our three locations are selected at the same latitude as the actual earthquake and to the south and north of the actual earthquake. Information on the random earthquakes is shown in Table 2.
To compare the anomalous variations in the actual seismic region with those in the random seismic regions, we calculated the S-K coefficients for PC2 in the random seismic regions and accumulated the anomalous tracks. The results are shown in Fig. 10.
As shown in Fig. 10, the cumulative number of anomalous tracks in random regions increases almost linearly with time. Before and after the earthquake, abnormal changes are not observed in the cumulative number of abnormal tracks. However, in the actual seismic research area, as the earthquake approaches, the cumulative number of abnormal tracks increases rapidly and recovers to a slow growth after the earthquake. This pattern is also consistent with the ''typical random process'' and ''typical critical system'' proposed by De Santis et al.: ''The cumulated value of a typical random process has a statistical linear increase. In particular, in case of critical phenomena, we would expect more frequent anomalies when they approach the critical point, and less frequent anomalies after'' [33]. This phenomenon further  verifies the correlation between the anomaly extracted by PC2 and the Ecuador earthquake.

V. DISCUSSION AND CONCLUSION
This paper attempts to use the advantages of Swarm satellite constellation monitoring to simultaneously analyze the magnetic field data of Swarm A and C satellites to extract the anomalies associated with the Ecuador earthquake. Recently, pre-earthquake anomaly analyses have been performed for the Ecuador earthquake using the Swarm satellite magnetic field data. The analysis of magnetic field data of daytime recorded by Swarm A shows that when only the tracks under quiet magnetic conditions are considered, the accumulation of anomalous tracks is significantly accelerated on the 9th day before the Ecuador earthquake [34].
However, in our results, the accumulation of anomalous tracks is significantly accelerated around the earthquake day. We speculate that the difference between our results and those previous studies may be due to two reasons. First, since the ionospheric variation caused by the sun during the daytime may overwhelm the small disturbance caused by an earthquake, we select data during local nighttime hours to avoid the unclear complexity and diversity of the data recorded during the daytime. Moreover, we apply PCA to the magnetic fields measured by Swarm A and C and separate the signal mainly affected by geomagnetic activity into PC1 through PCA and the signal that is mostly unaffected by geomagnetic activity into PC2 to extract the pre-earthquake anomaly without having to consider quiet or disturbed times. Through the accumulation of anomalous tracks, we also find that the inflection point of the sigmoidal fitting curve is surprisingly close to the time of the earthquake, which is similar to the results of De Santis et al [33]. In addition, the cumulative number of abnormal tracks in other random regions showed a linear increase, further supporting the correlation between the extracted anomalies and earthquakes.
Before the earthquake, fault activation leads to gas migration, and the resulting particles ionize the air. Then the variations of air conductivity lead to the variations of atmospheric electricity and finally induce variations in ionosphere including the electromagnetic before the earthquake. Due to the coupling mechanism of the LAIC (lithosphere-Atmosphere-Ionosphere Coupling), earthquake anomalies originate from the lithosphere and are then transferred to the atmosphere and the ionosphere [60]. The Swarm satellite is at the height of the ionosphere, Due to the LAIC mechanism, the Swarm satellite can detect the ionospheric magnetic field anomalies before the earthquake. Our results also suggest a physical coupling between the lithosphere and ionosphere.
In this paper, we consider Swarm A and C as two moving stations, and apply the PCA method to the data of two satellites to separate seismic-related signals from other signals. The results suggest that the multi-satellite combination analysis approach may represent a promising method for seismic precursor extraction.
KAIGUANG ZHU received the Ph.D. degree in earth exploration and information technology from Jilin University, Changchun, China, in 2002. Her current research interests include signal processing involved with airborne electromagnetic survey, and earthquake precursor analysis from both ground-based and satellite observations. KAIYAN LI received the bachelor's degree from the College of Instrumentation and Electrical Engineering, Jilin University, Changchun, China, in 2017, where she is currently pursuing the master's degree. Her current research interests include seismic precursor anomaly detection and satellite data processing.
MENGXUAN FAN received the bachelor's degree from the College of Instrumentation and Electrical Engineering, Jilin University, Changchun, China, in 2017, where she is currently pursuing the master's degree. Her current research interests include signal processing, seismic precursor anomaly detection, and satellites data processing.
CHENGQUAN CHI received the B.S. degree in automation and the M.S. degree in control theory and control engineering from Northeast Forestry University, Harbin, China, in 2009 and 2013, respectively. He is currently pursuing the Ph.D. degree in measurement technology and instrumentation with Jilin University, Changchun, China. His research interests include blind source separation, geophysics, and signal processing, especially earthquake precursor data.
ZINING YU received the bachelor's degree from the College of Instrumentation and Electrical Engineering, Jilin University, Changchun, China, in 2016, where she is currently pursuing the Ph.D. degree. Her current research interests include signal processing, seismic observation data analysis, and seismic precursor analysis.