Seismic Reﬂection Analysis of AETA Electromagnetic Signals

: Acoustic and electromagnetics to artiﬁcial intelligence (AETA) is a system used to predict seismic events through monitoring of electromagnetic and geoacoustic signals. It is widely deployed in the Sichuan–Yunnan region (22 ◦ N–34 ◦ N, 98 ◦ E–107 ◦ E) of China. Generally, the electromagnetic signals of AETA stations near the epicenter have abnormal disturbances before an earthquake. When a signiﬁcant decrease or increase in the signal is observed, it is difﬁcult to quantify this change using only visual observation and conﬁrm that it is related to an upcoming large earthquake. Considering that the AETA data comprise a typical time series, current work has analyzed the anomalism of AETA electromagnetic signals using the long short-term memory (LSTM) autoencoder method to prove that the electromagnetic anomaly of the AETA station can be regarded as an earthquake precursor. The results show that there are 2–4% anomalous points and some outliers exceeding 0.7 (after normalization) in the AETA stations within 200 km of the epicenter of the Jiuzaigou earthquake (M. 7.0) and the Yibin earthquake (M. 6.0) half a month before the earthquakes. Therefore, the AETA electromagnetic disturbance signal can be used as an earthquake precursor and for further earthquake prediction.


Introduction
Before a large earthquake (M. > 6.0) occurs, there are many abnormal phenomena related to the earthquake which are called earthquake precursors. The International Association of Seismology and Physics of the Earth's Interior (IASPEI) has a subcommittee that developed a list of such precursors, namely seismic quiescence, a decrease in radon concentration, other geochemical phenomenon, and changes in ground water levels [1]. Precursor phenomena not included in the IASPEI list have also been confirmed as potential earthquake prediction tools, such as anomalous electromagnetic field changes, abnormal animal behavior, and fault creep or continuous strain. With the improvement of computer computing performance, many data analysts have attempted to utilize earthquake-related signals with machine learning or deep learning models for earthquake prediction or warning. However, it is necessary to establish that the signal anomaly is related to seismic activities before using algorithms to analyze the signal.
The correlation analysis between seismic precursors and earthquake activities has always attracted the attention of geophysicists, geologists, and data analysts. Precursor signals with great correlation have stronger interpretability and even better performance in downstream tasks such as earthquake early warning. Some animals have evolved to develop a vibration-triggered early warning response, as outlined by Kirschvink, which acts in the short time interval between the arrival of P and S waves [2]. Kulahci utilized time series analysis to determine that stress accumulation preceding the earthquake could cause the radon anomaly [3]. Han, P. et al., who used a method named superposed epoch analysis (SEA), observed that the ULF geomagnetic anomalies are more likely seen 6-15 days before sizeable earthquake events (Es ≥ 108) [4]. The moving median method was used in total electron content (TEC), land surface temperature (LST), and aerosol optical depth (AOD) analysis, and the result in Sekertekin's study confirmed that those data are important potential precursors for earthquake prediction [5]. Genzano, N. executed a correlation analysis between significant sequences of thermal anomalies (SSTAs) and M. ≥ 6 Japanese earthquakes by using robust satellite techniques (RSTs), which address the constraints concerning space, time, and magnitude [6].
Many such methods have been used in the correlation analysis of seismic precursors. The LSTM autoencoder (LAE) is an unsupervised neural network that analyzes the anomaly information of the original data by reconstructing it, and this approach is often used for anomaly detection of temporal series. Feng et al. used the autoencoder method to explore the latent space distribution without anomaly analysis from past seismic events, which was used as input for a subsequent earthquake prediction model [7]. Paul, S. et al. developed an unsupervised model using the LAE algorithm to mimic the daily activities of employees to detect anomalous employee behavior of a community [8]. The algorithm was also used to evaluate 12 representative types of anomalies from 1555 robot-assisted feeding executions by Park's team, and, compared with five baseline detectors, the LAE achieved the best performance [9]. Hashimoto, S. et al. utilized the LAE method to prevent industrial accidents by taking the sequential skeleton map as the input [10]. LSTM autoencoder and one-class support vector machines were combined by researchers in synthetic data and real data of fashion retail, which comprised a multivariate time series, to detect anomalies; the obtained results led to better performance compared to other machine learning methods [11]. Khan et al. proposed a hybrid intelligent intrusion detection system, which was developed using the Spark MLlib and an LSTM autoencoder deep learning approach; the simulation results significantly outperformed existing intrusion detection approaches for the ISCX-UNB dataset [12]. A combination of LSTM and variational AE was used in Kim's research, which can diagnose abnormal network situations as well as ensure that the credibility of the diagnosis is estimated via the anomaly score-based negative log-likelihood [13]. Due to the timing of AETA electromagnetic signals, the aim of the present study was to confirm, using the LAE algorithm, that the AETA station has an abnormal electromagnetic disturbance before the earthquake, which can be considered as an earthquake precursor. The signals have the potential to be used as a tool for researchers to predict earthquakes.
In Section 2, we introduce the detailed parameters of the AETA equipment developed by our laboratory and the deployment of AETA stations in Sichuan and Yunnan. In Section 3, the method of LAE and interquartile range (IQR) and the construction of input samples are described. Section 4 describes the analysis of specific earthquake cases. Afterwards, the research work of this article is summarized in the Conclusions.

AETA System
Acoustic and electromagnetics to artificial intelligence (AETA) is a system designed by the IMS Laboratory of Peking University for earthquake monitoring and prediction. It consists of an underground signal acquisition probe, a ground data processing terminal, a monitoring cloud platform, and a data analysis system, as shown in Figure 1. The underground part contains electromagnetic and geoacoustic sensor probes, which are utilized to collect electromagnetic and geoacoustic signals in real time. First, the collected data are transmitted to the terminal with the function of data processing on the ground. Processed data are then received by the monitoring cloud platform through the network for further analysis by the researchers [14].
This article focuses on the anomaly analysis of an impending earthquake based on AETA electromagnetic data. The design of the AETA electromagnetic sensor is similar to that of the search-coil magnetometer instrument (IMSC) on the French DEMETER and that on the Chinese electromagnetic satellite, based on Faraday's electromagnetic theory. The induced electromotive force obtained from the magnetic field changing in the vertical Appl. Sci. 2021, 11, 5869 3 of 16 direction is inputted into the front-end signal processing circuit for amplification, and is then filtered and processed via digital-to-analog conversion (DAC). The processed signal passes through the back-end acquisition circuit, which has a strong processing speed, STM32F407 as the control unit, and a W5300 high-speed Ethernet communication interface chip as the data transmission module to meet the needs of real-time detection. The parameters of the electromagnetic sensor are as follows: the electromagnetic probe monitors 0.1 Hz~10 kHz, covering very low frequency and ultra-low frequency electromagnetic bands with a wide dynamic range of 0.1~1000 nT. The device has an 18-bit resolution, a sensitivity of >20 mV/nT@ 0.1 Hz-10 kHz, and a noise level of 0.1-0.2 pT/Hz@ (10 Hz-1 kHz). It has a low frequency sampling rate of 500 Hz and a full frequency sampling rate of 30 kHz [15]. This article focuses on the anomaly analysis of an impending earthquake based on AETA electromagnetic data. The design of the AETA electromagnetic sensor is similar to that of the search-coil magnetometer instrument (IMSC) on the French DEMETER and that on the Chinese electromagnetic satellite, based on Faraday's electromagnetic theory. The induced electromotive force obtained from the magnetic field changing in the vertical direction is inputted into the front-end signal processing circuit for amplification, and is then filtered and processed via digital-to-analog conversion (DAC). The processed signal passes through the back-end acquisition circuit, which has a strong processing speed, STM32F407 as the control unit, and a W5300 high-speed Ethernet communication interface chip as the data transmission module to meet the needs of real-time detection. The parameters of the electromagnetic sensor are as follows: the electromagnetic probe monitors 0.1 Hz~10 kHz, covering very low frequency and ultra-low frequency electromagnetic bands with a wide dynamic range of 0.1~1000 nT. The device has an 18-bit resolution, a sensitivity of >20 mV/nT@0.1 Hz-10 kHz, and a noise level of 0.1-0.2 pT/Hz@ (10 Hz-1 kHz). It has a low frequency sampling rate of 500 Hz and a full frequency sampling rate of 30 kHz [15].
So far, the IMS laboratory has installed more than 200 AETA devices nationwide, and 150 of them are in the regions of Sichuan and Yunnan (22° N-34° N, 98° E-107° E). The reason for this is that this area has high incidence of earthquakes, and more than 200 earthquakes with a magnitude of 3.5 or higher have occurred there since 2017, which requires further investigation. The distribution of all stations in Sichuan and Yunnan is shown in Figure 2. So far, the IMS laboratory has installed more than 200 AETA devices nationwide, and 150 of them are in the regions of Sichuan and Yunnan (22 The reason for this is that this area has high incidence of earthquakes, and more than 200 earthquakes with a magnitude of 3.5 or higher have occurred there since 2017, which requires further investigation. The distribution of all stations in Sichuan and Yunnan is shown in Figure 2.

Anomalies of Electromagnetic Signals before the Earthquake
The analyzed cases were selected from earthquakes (M. > 6.0) in the target area (22 • N-34 • N, 98 • E-107 • E) in recent years: 1. On 8 August 2017, a magnitude 7.0 earthquake occurred in Jiuzhaigou County, Aba Prefecture, Northern Sichuan Province. Figure 3 shows the changes of the electromagnetic disturbance signal of the AETA stations near the epicenter; 2. On 17 June 2019, a magnitude 6.0 earthquake occurred in Changning County, Yibin City, Sichuan Province. Figure 4 shows the changes in electromagnetic disturbance signals of the AETA stations near the epicenter. The red line at the top of the figure indicates the main shock and the aftershocks that followed.
It can be seen from Figures 3 and 4 that there are some signal anomalies before the earthquakes, and we argue that the abnormal change is related to the occurrence of earthquakes. It is difficult to quantify such anomalies directly from manual observations, and there may be other undetectable changes in addition to this obvious information. These subtle anomalies mixed with abnormal disturbances are difficult to detect with the naked eye. Therefore, a lot of potentially abnormal information will be lost directly through simple observation. In order to discover the hidden information, we chose the LSTM autoencoder model to capture anomalies.

LAE and IQR Method
The autoencoder (AE) is an unsupervised neural network that aims to learn the reconstruction close to its original input. In general, it consists of two parts, an encoder network and a decoder network. The encoder reduces the dimension of the input vector x ∈ R m and compresses it to the hidden vector h ∈ R n (n < m), and then the decoder raises the dimension of the obtained hidden vector h and maps it back to the original input space to obtain the reconstructed vectorx ∈ R m . The difference between the original input vector x and the reconstructed vectorx is called the reconstruction error L(x,x) [16]. The autoencoder updates the weight of the network by minimizing the reconstruction error L: The role of the autoencoder is not just to simply copy the input to the output. By constraining the latent space to have a smaller dimension than the input, i.e., n < m, the autoencoder is forced to learn the most salient features of the training data [17]. In other words, the autoencoder has a significant function in that it can find the regular or general pattern of the input data by reconstructing it.
Common autoencoders include the vanilla autoencoder [18], the convolutional autoencoder [19], the regularized autoencoder [20], and the LSTM autoencoder (LAE). As a form of sensor monitoring data, the AETA electromagnetic signal is a typical time series. The ability of LSTM to learn the patterns of long data makes it suitable for time series prediction or anomaly detection, and it is often used to capture the time dependence in multivariate data and determine whether data are out of distribution in a long time series. [21] In an LAE, both the encoder and decoder are LSTM networks. On the one hand, the fixed number of hidden units limit the model to learning the trivial mapping of the input; on the other hand, the same LSTM operation is performed recursively on the learned representation at any stage of decoding, which further hinders the model from learning the original mapping. This also explains why the LAE algorithm can achieve a reconstructed signal that presents almost all the input information but with a representation no better than the input [22]. The structure of the AE and LAE is shown as Figure 5.
The LAE created for AETA electromagnetic anomaly monitoring had two units corresponding to the encoding and decoding parts. The encoder consists of an LSTM layer with a size of 100 and a repeating module that is used to increase the dimension of the vector output by the former. There is also an LSTM layer and a flatten layer, which are used to transform the hidden layer vector from the encoder into reconstruction sequences with the same dimension as the input. The activation function applied for the LAE model was ReLU, which is based on a trial and error validation. We chose the Adam optimization as the optimizer and the mean square error (MSE) shown in Formula (1) as the loss function. The training data were shuffled for every epoch, and each station was trained for a sufficient number of epochs, about 200. The learning rate was initialized to 0.001. The weight decay was set to 0.005 and the dropout at 0.5, for all layers. The model stopped training data when the loss of validation set did not decrease for 50 consecutive generations.
Quartiles were always used to draw box plots in statistics; that is, all values were arranged from small to large and divided into four equal parts. The value at the position of the three division points is the quartile. Interquartile range (IQR) is a method in descriptive statistics that represents a set of data arranged in order, as well as the difference between the upper quartile Q 3 and the lower quartile Q 1 [23]. When the IQR method is utilized for the monitoring of abnormal phenomena, the outliers are usually defined as values less than Q lower or greater than Q upper , as the following formula: Q upper = Q 3 + 1.5IQR (4) words, the autoencoder has a significant function in that it can find the regular or general pattern of the input data by reconstructing it. Common autoencoders include the vanilla autoencoder [18], the convolutional autoencoder [19], the regularized autoencoder [20], and the LSTM autoencoder (LAE). As a form of sensor monitoring data, the AETA electromagnetic signal is a typical time series. The ability of LSTM to learn the patterns of long data makes it suitable for time series prediction or anomaly detection, and it is often used to capture the time dependence in multivariate data and determine whether data are out of distribution in a long time series. [21] In an LAE, both the encoder and decoder are LSTM networks. On the one hand, the fixed number of hidden units limit the model to learning the trivial mapping of the input; on the other hand, the same LSTM operation is performed recursively on the learned representation at any stage of decoding, which further hinders the model from learning the original mapping. This also explains why the LAE algorithm can achieve a reconstructed signal that presents almost all the input information but with a representation no better than the input [22]. The structure of the AE and LAE is shown as Figure 5. The LAE created for AETA electromagnetic anomaly monitoring had two units corresponding to the encoding and decoding parts. The encoder consists of an LSTM layer with a size of 100 and a repeating module that is used to increase the dimension of the vector output by the former. There is also an LSTM layer and a flatten layer, which are used to transform the hidden layer vector from the encoder into reconstruction sequences with the same dimension as the input. The activation function applied for the LAE model was ReLU, which is based on a trial and error validation. We chose the Adam optimization as the optimizer and the mean square error (MSE) shown in Formula (1) as the loss function. The training data were shuffled for every epoch, and each station was trained for a sufficient number of epochs, about 200. The learning rate was initialized to 0.001. The weight decay was set to 0.005 and the dropout at 0.5, for all layers. The model stopped The reason for the robustness of IQR is that up to 25% of the data can become arbitrarily distanced without greatly disturbing Q 3 or Q 1 so that outliers cannot affect this standard [24].

Sample Construction
The AETA signal is highly correlated with time, the quantity of data is large, and the information contained in a single point is insufficient. It is difficult for analysts who directly observe signals or utilize machine learning methods to process data to obtain valuable information, and they are likely to fall into a dilemma that they cannot resolve.
n represents the total number of AETA electromagnetic data for half an hour. We selected the electromagnetic disturbance data from 10 June 2017 to 10 December 2020, and, through the above formula, we obtained the average value of electromagnetic granularity data for 30 min. After that, the sliding window of the feature data processed from the AETA electromagnetic sensor was used as the input sequence, and the label had the same nature as the input with the size of 48, which was used to train the model more accurately, expand the training samples, and seek the abnormal information of the electromagnetic signal. Generally, since there are mostly samples with small magnitudes and too few instances with high magnitudes, it is necessary to consider how to handle class imbalance when performing earthquake prediction or magnitude classification tasks. However, the current study is an anomaly analysis of the AETA signal before the large earthquake, and there is no problem of sample imbalance. The sliding window method is shown in Figure 6, and a total of 58,302 samples were obtained. Appl. Sci. 2021, 11, x FOR PEER REVIEW 8 of 16 Figure 6. The sliding window method used to obtain samples.

Seismic Case Analysis
During the period from 10 June 2017 to 10 December 2020, 92 earthquakes with magnitudes above 4.0 occurred in the target area (22° N-34° N, 98° E-107° E) selected in this study. Two earthquakes with a magnitude over 6.0 were chosen for this study through analyzing the outliers of the difference between reconstructed and original series to verify that the abnormal AETA data can reflect an impending large earthquake. The epicenter of two major earthquakes and nearby stations is shown in Figure 7.  Table 1 shows the location and installation date of stations within 200 km of the earthquakes. Figure 8 shows the fitting effect of the LAE on the electromagnetic data of JZG and SP stations half a month before and

Seismic Case Analysis
During the period from 10 June 2017 to 10 December 2020, 92 earthquakes with magnitudes above 4.0 occurred in the target area (22 • N-34 • N, 98 • E-107 • E) selected in this study. Two earthquakes with a magnitude over 6.0 were chosen for this study through analyzing the outliers of the difference between reconstructed and original series to verify that the abnormal AETA data can reflect an impending large earthquake. The epicenter of two major earthquakes and nearby stations is shown in Figure 7 Table 1 shows the location and installation date of stations within 200 km of the earthquakes. Figure 8 shows the fitting effect of the LAE on the electromagnetic data of JZG and SP stations half a month before and one week after the earthquake. Figure 9 shows the curve fitted by the model on a specific day with or without anomalies, which can show the fitting ability of the model more clearly.
From Figures 8 and 9, it can be seen that the LAE algorithm has an excellent ability to reconstruct time series data, which can fit the rising or falling trend of original data well. Figure 8 shows that both stations showed anomalous phenomena before the earthquake day (8 August). The anomalies gradually disappeared, and the deviation value tended to be stable after that day. In particular, the JZG station showed some abnormalities on 1 August and 3 August, and the SP stations showed strong abnormalities on 31 July and 1 August. At other times, the difference between the reconstructed and the raw value was within the normal range. From the detailed diagram of a single day, the model shows a strong fitting ability and can reconstruct the input data well, while it can no longer find the trend of the data, showing a larger deviation from the original input.
In order to enhance the robustness of the conclusion, we counted the number of outliers in all six stations before and after the earthquake, so that the anomalies of each station appear more intuitively in mathematical statistics. The abnormal points are defined as follows: 0, i f |y −ŷ| ≤ y upper |y −ŷ|, i f |y −ŷ| > y upper (6) where y andŷ represent the input and reconstructed data point, respectively. The deviation is ∆y, and the IQR threshold value is y upper . When ∆y > 0, the point belongs to an abnormal point. Otherwise, we consider the point as a normal point. If the point ∆y > 0.7, the point belongs to "outliers over 0.7" at the same time. The statistical results are shown in Table 2. It can be seen in Table 2 that, before the earthquake, about 2% of the points had anomalies. Among them, MX stations had the most total outliers, and, for the class of outliers exceeding 0.7, SP showed the most with two points. There were only six abnormal points in QCX with no value exceeding 0.7. In summary, before the Jiuzhaigou earthquake, the electromagnetic data of all the AETA stations showed abnormal disturbances, and that the number of anomalous points was irrelevant to the distance from the station to the epicenter. The reason may be that the electromagnetic waves received by different stations contain different seismic information that depends on geological structure when the electromagnetic data are propagated underground.  When the Yibin M. 6.0 earthquake occurred, there were more than 30 AETA monitoring stations within 200 km of the epicenter, 27 of which were operating normally and had robust data. Table 3 shows the location and installation date of stations within 200 km of the earthquakes. Figure 10 shows the abnormal values after LAE fitting, which can more intuitively show the abnormality before and after the earthquake.  Figure 10 shows the statistical abnormal points at the PSFR, GAX and MBXB stations. Among them, IQR − Diff refers to ∆y in Formula (6), and IQR − Prop is the abnormal point obtained after we normalized the outlier and recalculated the IQR anomaly. We can see that some stations exhibited anomalies before the earthquake in the first and last pictures in Figure 10, and the middle section shows that some stations have abnormal data after the earthquake. The number or amplitude of the outliers obtained by Diff and Prop are similar, which further demonstrates the robustness of the method of extracting anomalies from the results of LAE processing. The IQR-Prop calculation method is defined as follows: y max = max(x n ) (8) y = ∆y − y min y max − y min (9) ∆y Prop = 0, i f y ≤ y upper y, i f y > y upper (10) where x n is the deviation series of LAE input and output, y is the normalized result of x n , y upper is the IQR anomaly threshold calculated after the sequence is standardized, and the anomaly value is ∆y Prop . The number of abnormal points at each station is shown in Table 4.

Comparison with Principal Component Analysis (PCA) Method
PCA is a method widely used in the field of data mining and data analysis. It has the advantages of separating different signals and effectively identifying relatively weak signals and can be used for pre-earthquake electromagnetic disturbance analysis and total electron content (TEC) anomaly detection [25][26][27]. Therefore, we compared the results of the anomaly analysis of the Jiuzhaigou earthquake with those of the LAE using the PCA method. The results of the PCA are shown in Figure 11. Figure 11 shows the thermogram of the abnormal values of each monitoring station from 1 August 2017 to 31 August 2017. The abscissa and the ordinate represent date and time, respectively. It can be seen from the figure that 5 days before the earthquake a continuous anomaly band with an increasing level of abnormality appeared at the JZG monitoring station at 5 ∼ 6 o'clock, and that the value of the abnormal points exceeded two in the 2 days before the earthquake, which was significantly higher than the other historical abnormal points in the background. The outlier point appeared at 21 ∼ 22 o'clock on the day of the earthquake, which means that the anomaly degree of the electromagnetic data at this time is greater than that at other times of the day. After the earthquake, the anomaly degree weakened generally, and the abnormal band appeared again at 5 ∼ 6 o'clock and lasted 15 days. However, other monitoring stations near the seismic epicenter have a relatively small degree of anomaly, and their thermograms do not constitute a continuous band anomaly. In short, compared with the PCA method, the LAE can be used to observe the anomaly time and amplitude of the station more intuitively and conveniently. At the same time, in terms of the ability to identify seismic anomalies, the LAE method is robust. For example, for stations such as SP, it can capture anomalies not found by the PCA method. Therefore, the LAE is more suitable for anomaly analysis of AETA data before earthquakes.

Conclusions
In this study, the LAE and IQR methods were utilized to analyze the electromagnetic signals of the AETA stations near the Jiuzhaigou and Yibin epicenter. The results show that the electromagnetic signals of the AETA stations are obviously abnormal before the two earthquakes. First, the time, value, and number of abnormal points in stations at different geographical locations are not consistent, which is probably due to the difference of geological structure and the distance from the epicenter. The stations that were further away from the epicenter may show stronger anomalies than those closer. In general, there was obvious abnormal information about half a month before the earthquake. Secondly, the stations near the seismic epicenter generally exhibited outliers accounting for 2 ∼ 4% of the total data, and some stations had outliers with anomaly values exceeding 0.7 before the occurrence of a large earthquake, which shows that the electromagnetic anomaly of the AETA station has a certain correlation with the earthquake. Finally, comparing the methods of LAE and PCA through the case of the Jiuzhaigou earthquake, we found that the abnormal performance of the station is more common based on the LAE method, and that the ability to respond to abnormalities is also stronger and more intuitive. While the PCA method has good earthquake reflection performance at the JZG station, the performance of the other stations is not effective, which further proves the robustness of the LAE method used in this study.
In summary, the experiment in this study proves that the electromagnetic disturbance anomalies of the AETA station could imply an upcoming large earthquake, and that AETA electromagnetic signals can be regarded as a type of earthquake precursor. According to the experimental results, we believe that if 2% of abnormal points appear in most monitoring stations within half a month, it means that a certain area has a higher risk of strong earthquakes.
In the future, we consider using a Generative Adversarial Networks (GAN) approach to generate several representative anomalous sequences. When the original signal matches these anomalous sequences, the system can automatically alert, which indicates that a certain area is facing the risk of high-level earthquakes. In addition, due to the complexity of actual seismic activity, it is difficult for us to directly estimate the geographical location of the epicenter from statistical values such as the number and amplitude of anomalous points. The use of AETA electromagnetic signals for earthquake prediction or early warning is awaiting further research in the future.