Analysis of Swarm Satellite Magnetic Field Data Before the 2016 Ecuador (Mw = 7.8) Earthquake Based on Non-negative Matrix Factorization

In this paper, based on non-negative matrix factorization (NMF), we analyzed the ionosphere magnetic field data of the Swarm Alpha satellite before the 2016 (Mw = 7. 8) Ecuador earthquake (April 16, 0.35°N, 79.93°W), including the whole data collected under quiet and disturbed geomagnetic conditions. The data from each track were decomposed into basis features and their corresponding weights. We found that the energy and entropy of one of the weight components were more concentrated inside the earthquake-sensitive area, which meant that this weight component was more likely to reflect the activity inside the earthquake-sensitive area. We focused on this weight component and used five times the root mean square (RMS) to extract the anomalies. We found that for this weight component, the cumulative number of tracks, which had anomalies inside the earthquake-sensitive area, showed accelerated growth before the Ecuador earthquake and recovered to linear growth after the earthquake. To verify that the accelerated cumulative anomaly was possibly associated with the earthquake, we excluded the influence of the geomagnetic activity and plasma bubble. Through the random earthquake study and low-seismicity period study, we found that the accelerated cumulative anomaly was not obtained by chance. Moreover, we observed that the cumulative Benioff strain S, which reflected the lithosphere activity, had acceleration behavior similar to the accelerated cumulative anomaly of the ionosphere magnetic field, which suggested that the anomaly that we obtained was possibly associated with the Ecuador earthquake and could be described by one of the Lithosphere–Atmosphere–Ionosphere Coupling (LAIC) models.


INTRODUCTION
Earthquakes are the most energetic phenomena in the lithosphere and on the Earth (Bolt, 1999). During the long-term earthquake preparation processing and at the largest energy release moment of earthquakes, some seismic-linked anomalies could occur in the lithosphere, atmosphere, and ionosphere (Parrot, 1995;Freund, 2000;Liu et al., 2015;Hattori and Han, 2018;Piscini et al., 2019).
Swarm is a satellite mission of the European Space Agency (ESA). The mission consists of three satellites (Alpha, Bravo, and Charlie) that are devoted to studying the geomagnetic field and its dynamics (Olsen et al., 2006;Bouffard et al., 2019). The magnetic field data of Swarm satellites are widely used in fields, such as studying the ionospheric current systems (Alken, 2016) and the magnetic storms (Wang et al., 2019), and constructing the magnetic field models (Finlay et al., 2015(Finlay et al., , 2016. Beyond this, the high-precision magnetic field data are also applied to study the ionosphere anomalies, which are possibly related to earthquakes Marchetti et al., 2019a). After the launch of the Swarm satellites at the end of 2013, some large earthquakes occurred. We focused on the 7.8 Mw earthquake that occurred in Ecuador (0.35 • N, 79.93 • W) on April 16, 2016.
Since satellite magnetic field data are affected by geomagnetic activity (Zhima et al., 2014;Perrone et al., 2018), in the routine analysis of the anomalies of earthquakes, the data under strong geomagnetic conditions are usually deleted. Zhang et al. (2009) studied the 2007 Chile M7.9 earthquake using the DEMETER satellite magnetic field data with Kp < 3 and Dst > −20 nT and found that the low frequency electromagnetic disturbances increased 1 week before the earthquake. Akhoondzadeh et al. (2019) used the magnetic field data collected under quiet geomagnetic conditions to study the 2017 Sarpole Zahab (Mw = 7.3) earthquake and observed ionospheric magnetic anomalies between 8 and 11 days prior to the occurrence of the event. De  used the Swarm satellite magnetic field data with |Dst| ≤ 20 nT and a p ≤ 10 nT to study 12 earthquakes with magnitudes from 6.1 to 8.3, and they observed some ionospheric magnetic field anomalies before most of the events (Yan et al., 2013;Marchetti and Akhoondzadeh, 2018;De Santis et al., 2019b;Marchetti et al., 2019a,b). Hattori et al. (2004) applied the principal component analysis to decompose the magnetic field data of three ground-based stations (including the data under strong geomagnetic activity). The first and the second principal components were found to be associated with the geomagnetic variation and man-made activity, respectively, and the residual third component contained the earthquake precursory signature. Then, in 2019, the principal component analysis was applied to the satellite magnetic field data by Zhu et al. (2019); they removed the component, which has a high correlation coefficient with the a p index, and extracted earthquake-related anomalies using the residual component.
Non-negative matrix factorization (Lee and Seung, 1999) is a parts-based matrix decomposition method that has the advantage of obtaining local features from the whole data, which differs from the principal component analysis. NMF has been successfully applied in various fields. Lee and Seung (1999) proposed NMF and applied the method to decompose face images into basis images, such as the mouth, the nose, and other facial parts. Smaragdis and Brown (2003) applied NMF to decompose complex piano music and obtained the spectral bases of each note and its weight in time. Mouri et al. (2009) applied NMF to decompose electromagnetic data with extremely low-frequency bands from multiple ground stations and obtained local signals (which are emitted by regional electromagnetic radiation sources) from the extremely low-frequency data. Since the impact of earthquakes occurs only near the epicenters. Based on the advantage of NMF that could obtain local features from the overall data, in this paper, we used the NMF to analyze the ionosphere magnetic field data of Swarm Alpha satellite, including the data collected under quiet geomagnetic conditions and those under strong geomagnetic conditions, to explore the 2016 Ecuador earthquake (April 16, 0.35 • N, 79.93 • W).
First, we briefly described the Swarm satellite magnetic field data and the NMF method. Next, we performed NMF to decompose the short-time Fourier transform (STFT) magnitude spectra of each track and obtained three basis features and their corresponding weights. We calculated the energy-entropy ratio of the three weight components. For the weight component with the highest energy-entropy ratio, we computed the cumulative number of tracks, which had anomalies inside the earthquakesensitive area. In addition, we analyzed the influence of geomagnetic activity and plasma bubbles on the cumulative anomaly. From the random earthquake study and low-seismicity period study, we explored whether the cumulative anomaly was obtained by chance. Finally, we studied the correlation between the ionosphere magnetic anomaly and the lithosphere activity.

Swarm Satellite Magnetic Field Data
The Swarm mission consisted of three identical satellites (Alpha, Bravo, and Charlie), which had been launched into nearpolar orbits on November 22, 2013 (Bouffard et al., 2019). The satellites, Alpha and Charlie, flew almost side-by-side (longitudinal separation of 1.4 • at the equator) at an altitude near 462 km (initial altitude) and the Bravo flew higher, near 511 km (initial altitude). The three satellites completed one track in about 90 min, accomplished nearly 15 day-and-night passes every 24 h, and drifted about 1 h in local time (LT), every 10 days. The objective of the mission was to provide the best survey and the temporal evolution of the geomagnetic field, obtain a spacetime characterization of the internal field sources in the Earth, while improving the understanding of the interior of the Earth and the Geospace environment (Friis-Christensen et al., 2006Olsen et al., 2006).
Each satellite of the Swarm mission was equipped with seven instruments. The main sensors included two magnetometers, an Absolute Scalar Magnetometer that provided measurements of the field intensity, and a Vector Field Magnetometer that measured the magnetic field from three different directions. These sensors made high-precision and high-resolution measurements of the strength, direction, and variation of the magnetic field. The observations are provided as Level-1b data, which are then calibrated and formatted time series (Olsen et al., 2013).
The observation data of vector magnetic field in L1b products are shown in two reference systems: (1) the instrument frame and (2) the North (X) East (Y) Center (Vertical, Z) frame. The magnetic field of the Y-East component could be affected by lithospheric activity and is less influenced by external magnetic disturbances (Pinheiro et al., 2011). In this study, we analyzed the Y-East component magnetic field data (Level-1b 1 Hz data) of the Swarm Alpha satellite to study the 2016 Ecuador earthquake.

Non-negative Matrix Factorization
The non-negative matrix factorization (Paatero and Tapper, 1994;Lee and Seung, 2000;Smaragdis, 2004;Cichocki et al., 2006) is a matrix factorization algorithm that approximates a given non-negative matrix as a product of two other non-negative matrices. The two decomposed non-negative matrices consist of a few basis vectors which contain the potential structures of the given non-negative matrix. This algorithm can obtain partsbased representations of non-negative data and make no further assumptions about their statistical dependencies. Thus, NMF could obtain local features from the overall data.
The goal of NMF is to approximate the non-negative matrix V as a product of two non-negative matrices, W and H. The approximate factorization can be written as in Equation (1).
where, V ∈ R ≥0, m×n is the original matrix, W ∈ R ≥0,m×r is the basis matrix, and H ∈ R ≥0,r×n is the weight matrix. W * 1 , . . . , W * r are the columns of matrix W and can be interpreted as basis vectors of matrix V. H 1 * , . . . , H r * are the rows of matrix H, which have a one-to-one correspondence with the columns of W and can be interpreted as the weights of these basis vectors. The rank r corresponds to the number of basis vectors, and is generally selected, such that (n + m)r < nm.
To minimize the error of reconstruction of V by WH, we use the cost function based on the generalized Kullback-Leibler divergence. We constrained the solutions of NMF problems by the minimum determinant constraint (Schachtner et al., 2009), which could achieve unique and optimal solutions in a general setting. The cost function is shown in Equation (2).
where, α is the balance parameter. When α remains small enough, the reconstruction error does not increase during an iteration step, and very satisfactory results are obtained.
To optimize this function, we used an efficient multiplicative update algorithm introduced by Lee and Seung (2000); the iterative updated rules of basis matrix, W and weight matrix, H are shown in Equations (3) and (4).
We stopped the iteration when the D detKL (V|W, H)is was smaller than the threshold Th or the number of iterations reached the upper limit, N iter . The formula for calculating Th is shown in Equation (5).
where, ε is the error factor. This parameter should be a small value and set according to the actual needs.
Overall, the NMF algorithm can be summarized as follows: For the initialization (Step 1), we used the Gauss random. According to NMF, each element in matrix V can be described by the corresponding row of basis matrix W and the column of weight matrix H, as shown in Equation (6).
where, Vij is the element in i row and j column of matrix V. From this formula, we can see that the element Vij in matrix V can be calculated by the corresponding i row of W and j column of H. Since NMF uses a small amount of data to represent the original data, matrix W contains the main vertical features of matrix V and matrix H contains the horizontal information of these features. When there exists a local feature in the original matrix V, which has a special vertical feature and concentrates in a certain local position in the horizontal direction, this feature can be represented by one of the basis vectors of W and its corresponding weight of H. We expect that the anomalies possibly produced by earthquakes are special, in terms of frequency and/or pattern, and occur not so far from the epicenters. Thus, some of these anomalies could have local features, and through NMF, we can potentially obtain the components which represent the anomalies possibly related to the earthquake.

Ecuador Earthquake and Data Selection
On April 16, 2016, at 23:58:36 Coordinated Universal Time (UTC), a strong earthquake with a magnitude of Mw = 7.8, occurred ∼27 km south-southeast of Muisne, in the province of Esmeraldas, offshore of the west coast of northern Ecuador. The Ecuador earthquake was "the result of shallow thrust faulting on or near the plate boundary between the Nazca and South American plates. At the location of the earthquake, the Nazca plate subducted eastward, beneath the South American plate at a velocity of 61 mm/year" [from the United States Geological Survey (USGS)]. The location of the Ecuador earthquake is marked as a red star (0.35 • N, 79.93 • W, depth 21.0 km) in Figure 1A.
We studied the tracks that crossed the rectangular region, with R = 2259.4 km (half the side length), which is shown as a black rectangle in Figure 1A. The distance, R is the Dobrovolsky radius strain (Dobrovolsky et al., 1979); R = 10 0.43M (M is the magnitude of the earthquakes) estimates the size of the effective precursor manifestation zone of an earthquake. We assumed that the rectangular region was affected by the Ecuador earthquake, and referred to it as the "earthquake-sensitive area, " while the outside area was not influenced by the earthquake. Moreover, the geomagnetic latitudes of the tracks that we investigated ranged from −50 • to + 50 • (we avoided aurora and polar interferences). Taking April 13-15, 2016, as an example, the tracks of the Swarm A satellites flying over the earthquake-sensitive area are shown in Figure 1B.
Since we expected the influence of earthquakes to occur only near the epicenters, and due to the advantage of NMF to obtain local features from the overall data, we applied NMF to all the observed satellite magnetic field data, and did not delete the data under strong geomagnetic activity. We analyzed 302 effective tracks from February 16 to May 16, 2016.

Swarm Alpha Satellite Magnetic Field Data Decomposition
In this part, we applied NMF to decompose the Y-East component magnetic field data of the Swarm Alpha satellite. First, we preprocessed the magnetic field data. For each track, the predicted value from the CHAOS-6 geomagnetic field model was subtracted from the magnetic field measurement to remove the core and the static crustal magnetic fields (Finlay et al., 2016). Then, we calculated the first difference of the residuals, along the tracks. The preprocessed magnetic field data for one track data is represented by the vector X as reported in Equation (7).
where, L is the number of samples along the track. After preprocessing, we performed STFT to X, setting the window length to be equal to 64 and the step length to be equal to 16 (with 1 Hz sampling rate of the data, the spatial resolution would be about 122 km). Then the STFT magnitude spectrum V was decomposed by NMF to three basis vectors, namely W * 1 , W * 2, and W * 3 , and their corresponding weights , H 1 * , H 2 * , and H 3 * , as shown in Equation 8, with the number of basis vector, r = 3. For convenience, we refer below the three basis vectors as W1, W2, and W3, and the three corresponding weights as H1, H2, and H3. In addition, we set the times of the iterative update threshold, N iter = 1000, the balance parameter α = 0.001 in Equations (2) and (4), and the error factor, ε = 0.0001 in Equation (5).
The magnitude spectrum, V describes the amplitudefrequency characteristics at different locations along the track. The W1, W2, and W3 describe the basis features (vertical structure) of V. H1, H2, and H3 describe the weight coefficients of the three basis features over location and could show where the features W1, W2, and W3 appeared along the track. If there exists a basis feature, Wx (x = 1, 2, or 3), the high value of the corresponding weight Hx concentrates inside the earthquake-sensitive area; this feature reflects the local feature of the area and it is possibly associated with the Ecuador earthquake. We regarded it as a kind of anomaly inside the earthquakesensitive area.
Taking one track on April 9, 2016 as an example, the NMF decomposition results of the Y-East component magnetic field data for this track is shown in Figure 2. The preprocessed magnetic field data of the track is shown in Figure 2A. The STFT magnitude spectrum V of the data is shown in Figure 2B. The three basis vectors, W1, W2, and W3 are shown in Figure 2C. The three weight components, H1, H2, and H3 are shown in Figure 2D.
As shown in Figures 2A,B, before applying NMF, the features of the preprocessed data occurred inside and outside the earthquake-sensitive area. After decomposition, in Figures 2C,D, for feature W1, the corresponding weight, H1 has high values inside the earthquake-sensitive area at ∼5 • N latitude [with the highest value of 0.5560 which is over 5 * RMS (Root Mean Square) of H1], and it has low values outside the area. For feature W2, the corresponding weight, H2 has slightly high values outside the earthquake-sensitive area (with the highest value of 0.4015, which is 0.78 of 5 * RMS of H2). For feature W3, the corresponding weight, H3 has a very high value at approximately 22 • S latitude outside the earthquake-sensitive area (with the highest value of 0.6831, which is over the 5 * RMS of H3) and slightly higher weights inside the area (with the highest value of 0.28244, which exceeds the 5 * RMS of H3 by half of the same value). From these decomposition results, for the example track, the high values of weight, H1 is concentrated inside the earthquakesensitive area. Thus, the basis feature of W1 reflects the local feature of this area and is possibly associated with the Ecuador earthquake. The high value of H1 is an anomaly inside the earthquake-sensitive area.
The basis vectors describe features and their corresponding weights describe where these features are located along the track. We utilized the weight H components to detect whether anomalies existed inside the earthquake-sensitive area, which are possibly related to earthquakes. The higher the concentration of the H components concentrated inside the earthquake-sensitive area, the more likely this component would have anomalies inside the earthquake-sensitive area. We evaluated the concentration of energy and entropy inside the earthquake-sensitive area for three H components of all the studied tracks, by the energy-entropy ratio. The energy-entropy ratio shows the energy and entropy of the earthquake-sensitive area over those of the whole track, as shown in Equation (9).
where, the points between P s and P e are located inside the earthquake-sensitive area, as shown in Figure 2D. Xsn is the time-domain reconstruction signal for the Hn component. Z sXsn is the Shannon entropy (Shannon, 2001) of the whole track for Xsn, and Z inXsn is the entropy inside the earthquake-sensitive area for Xsn. The average and standard deviation of the largest energyentropy ratio, the second largest energy-entropy ratio, and the smallest energy-entropy ratio for the three decomposed H components of all the studied tracks are 0.571 ± 0.106, 0.456 ± 0.114, and 0.340 ±0.103, respectively. The largest average value is 25.22% higher than the second largest average value, and 67.94% higher than the smallest average value. From this statistical result, we could see that there existed an H component among the three decomposed H components, whose energy and entropy is much more concentrated inside the earthquake-sensitive area than in the other two H components. Therefore, this H component is more likely to have anomalies inside the earthquake-sensitive area, which is possibly related to the Ecuador earthquake. We refer to this component as Hs1 (the H component with the largest energy-entropy ratio), and the remaining two H components as Hs2 and Hs3, by the descending order of their energyentropy ratio.

Definition of Anomalous Tracks
Here, we defined the anomalous tracks for the H components. We calculated the RMS of the Hsn ( n = 1, 2, or 3) component and set the threshold, K * RMS (K = 5). The values that exceeded the threshold are the anomalies of this Hsn component, and a track with one or more anomalies is an anomalous track for this Hsn component. The track is an inside anomalous track for this Hsn component, if one or more anomalies occurred inside FIGURE 3 | The cumulative number of inside anomalous tracks and those of outside anomalous tracks for components Hs1, Hs2, and Hs3, from 60 days before to 30 days after the earthquake. (A) The cumulative number of inside anomalous tracks for components Hs1, Hs2 and Hs3. The red curve is the cubic spline with 5-knots points using all points. The solid line indicates that we use the points during this period to do the linear fit. The dashed line indicates that we do not use the points during this period to do the linear fit. The day of the earthquake is represented as a red vertical dotted line. (B) The cumulative number of outside anomalous tracks for components, Hs1, Hs2, and Hs3. The day of the earthquake is represented as a red vertical dotted line.
the earthquake-sensitive area. The track is an outside anomalous track for this Hsn component, if one or more anomalies occurred outside the earthquake-sensitive area. For example, the track in Figure 2 is an inside anomalous track for Hs1 component, but an outside anomalous track for Hs3 component.

Results
In this part, we calculated the cumulative number of inside anomalous tracks for Hs1 component, from 60 days before the earthquake to 30 days after the earthquake, to study the Ecuador earthquake. In addition, the cumulative number of inside anomalous tracks for Hs2 and Hs3 components and that of the outside anomalous tracks for the three H components were also computed for comparison. Finally, we obtained six cumulative results over time, which are presented in Figure 3. The three inside cumulative results are shown in Figure 3A, and the three outside cumulative results are shown in Figure 3B.
Moreover, if there are nx anomalous tracks on the t th day, the cumulative sequence Nx(t) is increased by nx on the previous day. For the cumulative results without any value on the first day (day −60), we made 0 as the value of the first day (day −60). For the cumulative results without any value on the last day of the cumulative (day +30), we made the last value of the result as the value of this day (day +30). Thus, all the cumulative results For the inside cumulative result of Hs1 component, we performed sigmoid fit, 3-degree polynomial fit, 6-degree polynomial fit, cubic spline with 2-knot points, and cubic spline with 5-knot points (D'Errico, 2021), using all of the blue points. The two cubic spline fits are constrained to be monotonical. We compared the mean squared error, adjusted the coefficient of determination and akaike information criterion, and found that the cubic spline with 5-knot points has the best fitting effect. Thus, for the inside cumulative result of Hs1 component, we performed a cubic spline with 5-knot points using all of the points, as shown in the red curve in Figure 3A. Based on this fit, we selected the points from day −60 to day −24 and from day +23 to day +30, and performed a linear fit to show the anomaly part of the cumulative result and we referred to it as background linear fit of the Hs1 component.
In Figure 3A, the inside cumulative result of Hs1 component shows a clear acceleration before the earthquake, deviating from the previous linear growth, and then recovers after the event. Compared to this result, the inside cumulative result of Hs2 exhibits approximately linear growth, and the Hs3 component does not have an inside anomalous track, as shown in Figure 3A. Figure 3B, Hs1 component has only one outside anomalous track. The outside cumulative result of the Hs2 component presents roughly a linear growth. There are two clear increase groups in the cumulative result. We checked the Dst, a p , and AE indices of the anomalous tracks for these two groups. The group from-14 days to −11 days has three anomalous tracks, two of them with |Dst| > 20 nT and a p > 10 nT. The group from +10 days to +14 days has five anomalous tracks; the tracks in day +11, day +12, and day +14 present values of the AE index higher than 100 nT, 80 nT, and 300 nT, respectively. Moreover, the AE index is 299 nT 2 h before the anomalous track in day +11 and 3 h before the anomalous track in day +12. From these results, we infer that the anomalous tracks of these two groups are likely to be affected by the geomagnetic activity. The outside cumulative result of Hs3 component shows a quasi-linear trend.

As shown in
Moreover, on comparing the cumulative number of inside cumulative results with those outside for the three H components, we found that almost all anomalies of Hs1 component occur inside the earthquake-sensitive area. However, for the Hs2 component, some anomalies are located inside the earthquake-sensitive area, some anomalies are located outside this area, and all anomalies of Hs3 component are outside the earthquake-sensitive area. From these results, we found that the inside cumulative result of Hs1 component is more likely to reflect the anomalies of the Ecuador earthquake.
Then, we further analyzed the inside result of the Hs1 component. From Figure 3A, until 20 days before the main shock, the cumulative number of Hs1 has exhibited linear growth and then begins to show accelerated growth. The increased speed reaches its maximum 11 days before the main shock, which is shown as the maximum acceleration anomaly, A1 in Figure 3A. Meanwhile, the cumulative number deviates from the previous linear growth, and the deviation extent reaches its highest in 2 days after the Ecuador earthquake, which is shown as the maximum deviation anomaly, A2 in Figure 3A. Eight days after the main shock, the cumulative number stops increasing temporarily, which is consistent with the deceleration time of the aftershock cumulative seismicity for the northern and southern patches (Agurto-Detzel et al., 2019). Thus, we speculated that the gap between +8 days to +22 days is possibly affected by the decrease of the seismic activity. Then, 23 days after the main shock, the cumulative number recovers to its linear growth. The abnormal phenomenon of the inside cumulative result for Hs1 component that the cumulative number accelerated before the earthquake and recovered after it is consistent with the studies of De Santis et al. (2017, 2019c). These studies indicated that the cumulative number of anomalies for a critical system would accelerate when approaching its critical time and recover after a large event , that is, the 7.8 Mw earthquake.
We also found that the accelerated cumulative anomaly of Hs1 component is similar to the results of the other two studies on the Ecuador earthquake, which also used the magnetic field data of the Swarm satellites but excluded the influences of the geomagnetic activity. The maximum acceleration anomaly, A1 in Figure 3A corresponds to a study by Akhoondzadeh et al. (2018). Their study shows that the cumulative number of anomalous tracks accelerates ∼9 days before the Ecuador earthquake, when the magnetic field data collected under quiet geomagnetic conditions are considered. The maximum deviation anomaly, A2 in Figure 3A is consistent with the study of Zhu et al. (2019). Their study indicated that the increased speed of the cumulative number of anomalous tracks reaches its maximum around the time of the Ecuador earthquake, after using principal component analysis to remove the component affected by geomagnetic activity.
In addition, we performed some confutation analysis to prove that the accelerated cumulative anomaly is possibly related to the Ecuador earthquake rather than caused by other disturbances.

CONFUTATION ANALYSIS Geomagnetic Activity Influence Study
To study whether the accelerated cumulative anomaly of Hs1 component is affected by geomagnetic activity, we performed two types of research. First, we analyzed the correlation between the time-domain reconstruction signal of Hs1 component and the geomagnetic index and compared the results with those of the original time domain signal and the other two H components.
For each investigated track (from 30 days before the earthquake to 30 days after the earthquake), we computed the median value of the time-domain signal X before the NMF decomposition and those of the time-domain reconstruction signal for the H components, Xs1, Xs2, Xs3 and Xs2 + Xs3. We compared the variation of the median value for the signals with their corresponding a p index [from the International Service of Geomagnetic Indices (ISGI)], as shown in Figure 4.
In Figure 4, Xs1 is not significantly correlated to the a p index, but the X, Xs3 and Xs2 + Xs3 seem to be adequately correlated with the geomagnetic a p index. Then, we calculated the correlation coefficient between the a p index and X, Xs1, Xs2, Xs3, and Xs2 + Xs3, where the results show the values, a −0.76, −0.20, −0.40, −0.82, and −0.84, respectively. These results reflect that the time-domain signals, X, Xs3, and Xs2, + Xs3 have strong correlations with the a p index; the absolute values of these correlation coefficients are over 0.8. In contrast, the correlation coefficient between Xs1 and the a p index is very low; the absolute value of the correlation coefficient is 0.20. Thus, after the NMF decomposition, the signal influenced by the geomagnetic activity seems to be included in the Xs2 and Xs3 components. However, the Xs1 component is not likely to relate to the geomagnetic activity.
Second, we examined the influence of the geomagnetic activity on the anomalous tracks of Hs1 component, and compared the results with those of the other two components. We set the Dst index (from ISGI) classification standard and divided the spatial conditions into four levels, from L0 toL3. Then, we calculated the numbers and percentages of anomalous tracks for Hs1, Hs2, and Hs3 components, before the earthquake, at different Dst levels, as shown in Table 1.
From Table 1, with the increase of the Dst level, the percentage of anomalous tracks for the Hs1 component shows a slight fluctuation; it first decreases by 20% and then increases by 35%. The percentage of anomalous tracks for the Hs2 and Hs3 components continuously increase with stronger geomagnetic activity. The percentage of anomalous tracks for Hs2 increases by 85% and then increases by 103%. At the L0 level, the percentage of anomalous tracks for the Hs3 is similar to Hs1, but at the L 2 and L 1 levels, the results of Hs3 are approximately twice as large as those of Hs1. Therefore, the anomalous tracks of Hs2 and Hs3 components are likely to be influenced by the geomagnetic activity. However, the anomalous tracks of Hs1 component are almost not affected by the geomagnetic activity.
From these two studies, we come to know that the accelerated cumulative anomaly of Hs1 component is hardly influenced by the geomagnetic activity.
We also assess the solar activity index, namely the F10.7 index within the study time range. Although the solar activity index near the earthquake is high [max: 107 solar flux unit (sfu)], no sudden change is evident. In addition, a previous study showed that the solar radio flux parameter probably did not affect the seismo-ionospheric anomalies detected around the date of Ecuador earthquake . Therefore, we did not consider the influence of the solar activity index, F10.7 on the ionosphere magnetic anomalies.

Plasma Bubbles Influence Study
The ionospheric magnetic field data might be affected by smallscale ionospheric irregularities, namely plasma bubbles. To explore whether the accelerated cumulative anomaly of Hs1 component is affected by plasma bubbles, we checked the Bubble index and the Bubble probability of the anomalous tracks (Park et al., 2013). Both the Bubble index and the Bubble probability of the points inside the earthquake-sensitive area of the inside anomalous tracks for Hs1 component, from day−14 to day +2, at the critical accelerated phase are 0, which means that the data are not affected by the plasma bubbles. So, we can speculate that the accelerated cumulative anomaly of Hs1 component is not probably caused by the plasma bubble.

Random Earthquake Study and Low-Seismicity Period Study
In this section, we performed a random earthquake study proposed by Parrot (2011) and low-seismicity period study to analyze the relationship between the Ecuador earthquake and the accelerated cumulative anomaly of the magnetic field data.
We randomly selected four pseudo earthquake epicenters around the world, which cover different latitude positions (excluding high latitudes), and studied them with the same time duration and the same size of earthquake-sensitive area as those of the Ecuador earthquake, without earthquakes over 5.5, as shown in Table 2 and Figure 5. In particular, we noted that the latitude of the Random 1 epicenter is similar to that of the Ecuador earthquake epicenter and the Random 2 epicenter is in a region that sometime has active ionospheric irregularities (Yizengaw and Groves, 2018). The LT of the four random pseudo earthquakes is almost the same as that of the Ecuador earthquake.
We also studied two low-seismicity periods with the same location and earthquake-sensitive area as those of the Ecuador earthquake; no earthquake with a magnitude over 5.8 occurred in these two periods, as shown in Table 2. Although the LT of the data for the two cases are not exactly same with those for the Ecuador earthquake, they cover the LT (sweeps from 21 to 20 LT for nighttime tracks and from 9 to 8 LT for daytime tracks) at the increased phase of the cumulative results.
Likewise, we calculated the cumulative number of the inside anomalous tracks of Hs1 component for the four random pseudo earthquakes and two low-seismicity periods, and then compared them with the cumulative results of the Ecuador earthquake, as shown in Figures 6, 7.
In Figures 6, 7, the cumulative number of the inside anomalous tracks for four random pseudo earthquakes and two low-seismicity period cases show approximately linear growth. To avoid the influence of LT, we checked the cumulative results of Low 1 and Low 2 from 21 (9) to 20 (8) LT, as shown by the green curves in Figure 7, and we found that they did not show accelerated increase. Thus, in non-earthquake regions or during non-earthquake periods, the cumulative results exhibited linear increase trend. The results are consistent with the "typical random process, " proposed by De Santis  2017); for a random process, the accumulated value is expected to show a statistically linear increase. However, for the actual Ecuador earthquake, as the earthquake approaches, the cumulative number of anomalous tracks accelerates growth and then recovers after the event, confirming that this trend seems to be related to the seismic event.
Moreover, from Figure 6 before the day of the earthquake, the average number of anomalous tracks for the four random earthquakes is ∼10, while the number of anomalous tracks of the real Ecuador earthquake is ∼25. This result shows that the number of anomalies for non-earthquake region cases is less than that of the real earthquake, which possibly supports the earthquake source of at least a subset of the anomalies (15 during the 2 months before the Ecuador earthquake), which is also consistent with the study of Parrot (2011).
These results further verify that the accelerated cumulative anomaly of Hs1 component of the ionosphere magnetic field is not obtained by simple chance or influenced by the study time period, geographical location, and LT, and it is possibly associated with the Ecuador earthquake.

Cumulative Benioff Strain S Study
In this section, we studied the correlation between the accelerated cumulative anomaly of the ionosphere magnetic field data and the lithosphere activities. Considering the cumulative effect of a series of N earthquakes before a large earthquake, we calculated the cumulative Benioff strain S (Benioff, 1949;De Santis et al., 2019a), which could estimate the strain-rebound increment by the earthquake energy, to study the lithosphere activity before the Ecuador earthquake, and compared it with the cumulative result of Hs1 component for the ionosphere magnetic field data.
Because of the limitations in the detection capability of the seismograph network, some weak earthquakes are not recorded. We used the maximum curvature technique (Xie et al., 2019) to estimate the smallest magnitude Mc that could be completely detected in the seismic-relevant area (with half large size of the earthquake-sensitive area of the Ecuador earthquake). The estimated complete magnitude is Mc = 4.3 by the seismic events (from USGS) in 2016. Then we selected the seismic events with magnitudes ≥4.3 and hypocentral depths ≤50 km, from February 16, 2016 until the moment of the Ecuador earthquake, in the seismic-relevant area to calculate the cumulative Benioff strain S.
For each selected seismic event, we calculated the released energy as explained in Equation (10) (Han et al., 2014).
where, Mi is the magnitude of the ith earthquake. Then we computed the cumulative Benioff strain S as explained in Equation (11)  .
where, nc is the total number of the earthquakes. The cumulative Benioff strain S around the epicenter of the Ecuador earthquake and the cumulative number of the inside anomalous tracks of Hs1 component for ionospheric magnetic field data are shown in Figure 8.
In Figure 8, until 8 days before the main shock, the cumulative Benioff strain S shows an approximate linear increase. After this day, the cumulative result shows a trend of accelerated growth and it continues until the moment of occurrence before the Ecuador earthquake. The clustering of earthquakes before a large earthquake might be a result of several physical mechanisms operating in the seismogenic crust (Shcherbakov et al., 2019). Therefore, this acceleration phenomenon may correspond to the energy release before the Ecuador earthquake and may be connected to the lithosphere activities before the earthquake, that is, the seismic activation of the fault system . Comparing the cumulative Benioff strain S and the cumulative results of Hs1 component of the ionosphere magnetic field data, we found that they have similar accelerate growth trend before the earthquake. And, as the earthquake approached, the anomaly phenomena in both the ionosphere and lithosphere became more and more obvious. Meanwhile, the acceleration day (8 days before the earthquake) of the cumulative Benioff strain S is near the day (11 days before the earthquake) when the increased speed of the cumulative number of the magnetic field data reached its maximum. This consistency indicates that the accelerated cumulative anomaly of Hs1 component for the ionosphere magnetic field data are possibly associated with the lithosphere activities, which are related to the Ecuador earthquake.
The correspondence between the accelerated cumulative anomaly of the ionosphere magnetic field data and the lithosphere activities is consistent with the study on the Nepal earthquake by De . This study indicates that the cumulative number of ionosphere magnetic anomalies and those of M4+ earthquakes have similar trends before and after the earthquake. By De Santis, this similar behavior between them supports the LAIC and the hypothesis that "the noticed magnetic anomalies in Swarm data are mostly of internal origin, due to a LAI coupling." According to LAIC models, the anomalies, which occurred in different layers before strong earthquakes, could be explained as a synergy between the lithosphere, atmosphere (included the Earth's surface), and ionosphere processes  and anomalous variations which are usually named as medium/short-term earthquake precursors. From the LAIC model, before the earthquake, fault activation releases some positive charged holes (Freund, 2011) or leads to gas migration including radon emanation (Pulinets and Ouzounov, 2011). Then the release of radon induces the ionization of the atmosphere. Formation of large ion clusters led to variations in the atmospheric electricity which is the main source of ionospheric anomalies over seismically active areas including the electromagnetic anomaly before the earthquake. Other coupling mechanisms have also been hypnotized; for example, a complete electric coupling induced by a change in the ground resistivity due to the variation of the strain on the fault (Kuo et al., 2014), or even a surface warming could produce an acoustic gravity wave before the occurrence of the earthquake (Hayakawa, 2011).
This correspondence between the anomalous behavior in the ionosphere and that in the lithosphere is consistent with the synergy among the processes of the different layers, and their anomalies are explained by one of the LAIC models. Therefore, our result supports the LAIC effects and we suggest that the Ecuador earthquake possibly involved a physical coupling between the lithosphere and the ionosphere in its preparation phase, and the accelerated cumulative anomaly of Hs1 component for the ionosphere magnetic field data is likely to originate in the lithosphere.

Investigation of Anomalies During the Accelerated Increase Phase
We investigated the anomalies extracted during the accelerated increase phase, from 26 days before the earthquake to the day of the earthquake (the anomalies in these days are more possibly related to the earthquake).
We checked the extracted anomalies one by one and found that the anomaly at track 5 on April 2, 2016, 14 days before the earthquake, is the closest anomaly to the epicenter of the Ecuador earthquake (<500 km away from the epicenter). The preprocessed East component of magnetic field data X for this anomalous track, the time-domain reconstruction signals, Xs1, Xs2, and Xs3 of the three decomposition components of this track, the feature Ws1, the weight, Hs1 and the location of this anomaly are shown in Figure 9. From Figure 9, it is clear that in Xs1, the anomaly lasted from −7 • latitude to +7 • latitude and the maximum value of this anomaly is almost at the same latitude of the epicenter. According to Ws1, the frequency with the highest amplitude of this anomaly is around 0.36 Hz. According to Hs1, the weight of this anomaly is surely concentrated inside the earthquake-sensitive area. Compared with Xs1, the original data X has two obvious anomalies. The northern one centered about at 0 • latitude is found also in Xs1, while a double pattern is found in Xs2 and Xs3. Although these two anomalies in X are nearer and symmetrical to the geomagnetic equator (around −10 • geographic latitude), the LT is 21 o'clock, they are not directly to be affected by the sunset. Moreover, the NFM technique underlines that the northern anomaly is different from the standard double pattern, which could be a residual of the daily interaction between the Sun and the ionosphere. From these results, the anomaly of Xs1 is more likely related to the earthquake.
In addition, we analyzed the relationship between the amplitude of anomalies and their distance from the epicenters, from 26 days before the earthquake to the day of the earthquake. The smallest distance from the epicenter is 465 km and the largest distance from the epicenter is 2,792 km. We divided the anomalies into three groups, according to their distance from the epicenter. For each group, we calculated the average energy of the anomalies. The relationship between the average energy of the anomalies and their distances from the epicenter is shown in Figure 10A. From Figure 10A, we can see that as the location of the anomaly is farther from the epicenter of the earthquake, the average energy of the anomaly decreased. This result is consistent with one of the standards for excellent earthquake precursors proposed by the International Association of Seismology and Physics of Earth's Interior (IASPEI) in 1991 and 1997 (Wyss, 1991(Wyss, , 1997Wyss and Booth, 1997). This standard observes that "the amplitude of the observed anomaly should bear a relation to the distance from the eventual mainshock, " and it is also affirmed by Rikitake and Yamazaki (1985).
The relationship between the amplitude and the time of occurrence of the anomalies was also investigated. The anomalies were divided into three groups, according to their time of occurrence. The average energy of the anomalies for each group was computed. The relationship between the average energy of the anomalies and their occurrence time is shown in Figure 10B. According to Figure 10B, at the advent of earthquake, the average energy of the anomalies increased. This result is reasonable, since generally with the advent of the earthquake, the anomalies should be stronger.
These two analyses further prove that the anomalies we extracted at the accelerated increase phase before the earthquake are likely to be related to the earthquake.

CONCLUSIONS
In this paper, we used the NMF to analyze the magnetic field data of Swarm Alpha satellite and explored the 2016 Mw 7.8 Ecuador earthquake, including all the observation data, regardless of the strength of the geomagnetic activity. After decomposition, we obtained three H components, and found that almost all the anomalies of one of the H components, that is the Hs1 component, occurred inside the earthquake-sensitive area, and the inside cumulative result of this component showed a clear acceleration before the Ecuador earthquake and recovered after it, which obeys the power-law behavior of a critical system . By different analyses, we excluded the influence of the geomagnetic activity and the plasma bubbles on the accelerated cumulative anomaly, and proved that this anomaly was not obtained by simple chance. Moreover, we found that the accelerated cumulative anomaly for the ionosphere magnetic field data was possibly associated with the lithosphere activities before the Ecuador earthquake, which also corresponds to the LAIC effect. By investigating the anomalies, we further found that as the location of the anomaly is farther from the epicenter of the earthquake, the average energy of the anomaly decreased, and with the advent of the earthquake, the average energy of the anomalies increased. Thus, the anomalies that we extracted are possibly related to the Ecuador earthquake.
These abovementioned results show that the NMF method has the capacity to detect the local feature of the earthquakesensitive area. For the Ecuador earthquake, based on NMF decomposition, by using all the observed magnetic field data without considering the geomagnetic activity, we obtained a weight component, whose cumulative result has an accelerated anomaly that is possibly related to the earthquake and not affected by the geomagnetic activity. This type of method could introduce a new perspective and analyze as much observation data as possible to study the earthquakes.
We suggest that our paper promotes the development from the search of anomalies that are likely to be related to earthquakes at the quiet geomagnetic conditions (which is still a difficult task) to the search of the anomalies that are likely to be related to earthquakes without considering the geomagnetic activity. We propose to address this problem with the use of NMF, but we do not exclude other methods that can more efficiently achieve the same goal.
We remind that the objective of these present analyses is not a retrospective prediction of an earthquake, but to emphasize that by NMF, we can use more data to analyze earthquakes. We would expect that using the satellite magnetic field data without considering the geomagnetic condition might be a better and more comprehensive method to understand earthquakes than using the data at quiet geomagnetic conditions alone.
In addition, further study should analyze the regular patterns among the local feature of the earthquake-sensitive area for different tracks. More case studies should be undertaken especially the investigation of other large earthquakes will be an important matter for future work.

AUTHOR CONTRIBUTIONS
KZ conceived and designed the project, guided the writing of the manuscript, and revised the manuscript. MF worked on the data analysis, the interpretation of the data, and wrote the manuscript. XH and KL contributed to data acquisition and revised the manuscript. DM contributed to data interpretation, gave some guiding opinions, and revised the manuscript. ZY, CC, HS, and YC contributed to the data interpretation and revised the manuscript. All authors contributed to the revision of the manuscript, read, and approved the submitted version.