Pre-Seismic Anomaly Detection from Multichannel Infrared Images of FY-4A Satellite

Abstract

Research of seismic infrared remote sensing has been undertaken for several decades, but there is no stable and effective earthquake prediction method. A new algorithm combining the long short-term memory and the density-based spatial clustering of applications with noise models is proposed to extract the anomalies from the multichannel infrared remote sensing images of the Fengyun-4 satellites. A statistical analysis is used to validate the correlation between the anomalies and earthquakes. The results show that the correlation rate is 64.29%, the hit rate is 68.75%, and the probability gain is about 1.91. In the Madoi and YangBi earthquake cases, the infrared anomaly detected in this paper is correlated with the TEC anomaly found in the previous research. This indicates that it is feasible to combine multi-source data to improve the accuracy of earthquake prediction in future studies.

1. Introduction

Earthquake prediction is a complicated issue [1]. The application of multi-source data may benefit to collect more detailed surface and atmospheric information and interpret pre-seismic anomalies [2]. Researchers have found multiple types of pre-seismic anomalies, including crustal deformation [3], infrared radiation [4], temperature [5,6], humidity [7], electromagnetic field [8,9], atmospheric composition [10], and so on. Huang [11] found changes in the seismicity pattern before a strong earthquake.

The previous research on seismic infrared remote sensing includes three aspects: the mechanism of anomaly [12,13,14], anomaly detection methods [15], and correlation analysis [16]. According to experiments and observations, changes in the earth’s crust, its atmosphere, and the earth’s electric field may cause anomalies of infrared radiation [14,17]. The data from the NOAA [18], Terra [19], Aqua [20], and Fengyun satellites [4,21,22] could capture infrared images continuously over a large area for a long time and have been applied to analyze earthquake precursors. The earth’s system is complex and infrared radiation is affected by several factors [23]. The changes in infrared radiation caused by earthquakes are difficult to observe directly. It is necessary to eliminate the background field and high-frequency noise before anomaly detection. Many signal processing methods have also been used to detect this pre-seismic anomaly, such as robust satellite techniques (RST) [24,25], wavelet transforms [26], power spectrum [4], and other classical signal processing methods. Machine learning has also been used to detect anomalies related to earthquakes. Akhoondzadeh et al. [27] compared the classical methods and artificial intelligence in the pre-seismic anomaly detection from infrared data. Researchers have used the Kalman filter [28], auto-regressive integrated moving average [29], support vector machines [30], neural networks, and genetic algorithm [31] to predict the time series of infrared or total electron content (TEC) data for pre-seismic anomaly detection. Zhai et al. [32] used the ARIMA model to detect the infrared anomalies before two earthquakes. Jing et al. [33] found anomalies in microwave brightness temperature and outgoing longwave radiation before the Madoi earthquake. Yang et al. [34] found the thermal infrared brightness temperature anomalies before this event.

Most of the infrared data used for earthquake prediction are long-wave infrared images, and a few are medium-wave infrared images [35]. In this paper, the multichannel infrared remote sensing images from the FY-4A satellite were applied, which included long-wave and medium-wave infrared images, and a new algorithm based on the long short-term memory (LSTM) and the density-based spatial clustering of applications with noise (DBSAN) models was proposed to detect the anomalies. The anomalies and earthquakes were analyzed statistically to validate the performance of this method in earthquake prediction.

2. Materials and Methods

2.1. Data and Study Area

The Fengyun-4 satellites are China’s latest generation of geostationary-orbit meteorological satellites. The first one (FY-4A) among the Fengyun-4 satellites was launched in December 2016 and has been providing data since 2018. The FY-4A satellite has two main optical detectors. One is the advanced geostationary radiation imager (AGRI), and the other is the geostationary interferometric infrared sounder (GIIRS). The AGRI of FY-4A could obtain images from 14 channels, including visible and infrared light [36]. The wavelengths of these channels are shown in

Table 1. The wavelengths of 14 channels of the advanced geostationary radiation imager (AGRI).

Channel	$\mathbf{Wavelength} \left(\mathsf{\mu }\mathbf{m}\right)$
1	0.45~0.49
2	0.55~0.75
3	0.75~0.90
4	1.36~1.39
5	1.58~1.64
6	2.10~2.35
7	3.50~4.00
8	3.50~4.00
9	5.80~6.70
10	6.90~7.30
11	8.00~9.00
12	10.30~11.30
13	11.50~12.50
14	13.20~13.80

. In this paper, the data of the FY-4A satellite from 2019 to 2021 were used to analyze the multichannel infrared anomalies before earthquakes. The former two-year (2019–2020) data were used to train the model, and the last one-year (2021) data were used to detect anomalies. The radiations emitted by the earth’s surface are mainly medium-wave infrared and long-wave infrared light, so only nighttime data from the seventh to the fourteenth channel were applied to avoid the influence of sunshine.

China is one of the countries with frequent earthquakes. The China earthquake networks center (www.ceic.ac.cn (accessed on 11 August 2022)) provided the earthquake list. There were 106 earthquakes with a magnitude over 5 in China and its surrounding area (from 70°E to 140°E, from 0°N to 60°N) in 2021. The seismic information is collected in Supplementary Materials. The geographic locations of the epicenters are shown in Figure 1. The strongest earthquake in the study area is the M 7.4 Madoi earthquake on 22 May 2021 (UTC+8). The Yangbi earthquake with a magnitude of 6.4 occurred about 4 h before the Madoi earthquake. There may be a series of earthquakes in the same region during a short term. The strongest earthquake in the earthquake sequence is the mainshock, and multiple foreshocks or aftershocks may occur before and after it. Taking the Madoi M7.4 and Yangbi M6.4 earthquakes in May 2021 as examples, these two earthquakes and their foreshock and aftershocks are shown in

Table 2. Seismic list around Madoi M7.4 and Yangbi M6.4 earthquakes.

Time (UTC+8)	Magnitude	Latitude (°N)	Longitude (°E)	Depth (km)	Type
22 May 2021, 10:29:34	5.1	34.85	97.5	10	Aftershock
22 May 2021, 02:04:11	7.4	34.59	98.34	17	Mainshock
21 May 2021, 22:31:10	5.2	25.59	99.97	8	Aftershock
21 May 2021, 21:55:28	5.0	25.67	99.89	8	Aftershock
21 May 2021, 21:48:34	6.4	25.67	99.87	8	Mainshock
21 May 2021, 21:21:25	5.6	25.63	99.92	10	Foreshock

. The Yangbi earthquake with a magnitude of 6.4 occurred at 21:48:34 on 21 May 2021 (UTC+8). The epicenter was at 99.87°E, 25.67°N. The depth of the hypocenter was 8 km. There was one foreshock and two aftershocks with a magnitude of over five around it. The Madoi earthquake with a magnitude of 7.4 occurred at 02:04:11 on 22 May 2021, (UTC+8). The epicenter was at 98.34°E, 34.59°N. The depth of the hypocenter was 17 km. There was one aftershock with a magnitude of over five around it. The study area in this paper is the region in the red box in Figure 1 (from 90°E to 110°E, from 17°N to 42°N), and there were 20 earthquakes with a magnitude of over 5, including 16 mainshocks.

2.2. Anomaly Detection

2.2.1. Long Short-Term Memory (LSTM)

The data with a spatial resolution of 4 km are used in this paper. The changes in infrared radiation in the nearby region are similar, and the seismic infrared anomalies commonly covered a large region [37]. The original data were down-sampled with an interval of 5 pixels to improve the efficiency of data processing. The brightness temperatures were obtained by looking up the calibration table, and the everyday average was calculated with the brightness temperature during nighttime (from 23:00 to 04:00 UTC+8). The preprocessed data are four-dimension data (two-dimensional geographic coordinates, time, and channel).

In this paper, the model based on the LSTM is used to predict the multichannel infrared data. The model comprised two layers (an LSTM layer and a fully connected layer). The structure of the LSTM layer is shown in Figure 2 [38]. $X_{t,m,n}$ is the input variable comprised of three-day and eight-channel data, as Equations (1) and (2) show [39]:

$$X_{t,m,n}=\left[\begin{array}{c}{S}_{t-3,m,n}\\ S_{t-2,m,n}\\ S_{t-1,m,n}\end{array}\right],$$

(1)

$$S_{t-k,m,n}=\left[B{T}_{t-k,m,n,7},B{T}_{t-k,m,n,8},\dots ,B{T}_{t-k,m,n,14}\right],$$

(2)

where $B{T}_{t-k,m,n,p}$ denotes the brightness temperature of the data at the mth line and the nth column at the $t-k$ moment in the pth channel. In Figure 2, $H_{t,m,n}$ is the output of the LSTM layer at the t moment. It was transferred to the state at the $t+1$ moment, and it was also inputted into the next layer to calculate the predicted value ${\widehat{Y}}_{t,m,n}$, which has the same dimension as the output variables $Y_{t,m,n}$, as Equations (3) and (4) show:

$${\widehat{Y}}_{t,m,n}=\left[{\widehat{BT}}_{t,m,n,7},{\widehat{BT}}_{t,m,n,8},\dots ,{\widehat{BT}}_{t,m,n,14}\right],$$

(3)

$$Y_{t,m,n}=\left[B{T}_{t,m,n,7},B{T}_{t,m,n,8},\dots ,B{T}_{t,m,n,14}\right],$$

(4)

where $B{T}_{t,m,n,p}$ denotes the brightness temperature of the data at the mth line and the nth column at the t moment in the pth channel; ${\widehat{BT}}_{t,m,n,p}$ denotes the predicted value at the mth line and the nth column at the t moment in the pth channel.

The root-mean-square error of the predicted and observed values is shown in Equation (5) [32]:

$$RMS{E}_P=\sqrt{\frac{{{\displaystyle \sum }}_{m=1}^M{{\displaystyle \sum }}_{n=1}^N{{\displaystyle \sum }}_{t=1}^T{\left({\widehat{BT}}_{t,m,n,p}-B{T}_{t,m,n,p}\right)}^2}{M\times N\times T}},$$

(5)

where $RMS{E}_P$ is the root-mean-square error in the pth channel; M and N are the numbers of the lines and columns of the image; and T is the time length of the data.

2.2.2. Density-Based Spatial Clustering of Application with Noise (DBSCAN)

The difference between the observed and predicted values is shown in Equation (6):

$$\begin{array}{c}\Delta Y_{t,m,n}=Y_{t,m,n}-{\widehat{Y}}_{t,m,n}\hfill \\ =\left[B{T}_{t,m,n,7}-{\widehat{BT}}_{t,m,n,7},B{T}_{t,m,n,8}-{\widehat{BT}}_{t,m,n,8},\dots ,B{T}_{t,m,n,14}-{\widehat{BT}}_{t,m,n,14}\right]\hfill \\ =\left[\Delta B{T}_{t,m,n,7},\Delta B{T}_{t,m,n,8},\dots ,\Delta B{T}_{t,m,n,14}\right]\hfill \end{array}$$

(6)

where $\Delta B{T}_{t,m,n,p}$ is the difference between the observed and predicted values at the mth line and the nth column at the t moment in the pth channel. The data in different channels were normalized, as Equations (7)–(10) show:

$$R_{t,m,n}=\left[R_{t,m,n,7},R_{t,m,n,8},\dots ,R_{t,m,n,14}\right],$$

(7)

$$R_{t,m,n,p}=\frac{\Delta B{T}_{t,m,n,8}-{\mu }_p}{{\sigma }_p},$$

(8)

$${\mu }_p=\frac{{{\displaystyle \sum }}_{m=1}^M{{\displaystyle \sum }}_{n=1}^N{{\displaystyle \sum }}_{t=1}^T\Delta B{T}_{t,m,n,p}}{M\times N\times T},$$

(9)

$${\sigma }_p=\sqrt{\frac{{{\displaystyle \sum }}_{m=1}^M{{\displaystyle \sum }}_{n=1}^N{{\displaystyle \sum }}_{t=1}^T{\left(\Delta B{T}_{t,m,n,p}-{\mu }_p\right)}^2}{M\times N\times T}},$$

(10)

where $R_{t,m,n}$ is the normalized data; the ${\mu }_p$ and ${\sigma }_p$ are the average and standard deviation of $\Delta B{T}_{t,m,n,p}$. The DBSCAN model could be used to detect anomalies. The Euclidean distance $D_{t_1,m_1,n_1,t_2,m_2,n_2}$ between two points ($P_1\left(t_1,m_1,n_1\right)$ and $P_2\left(t_2,m_2,n_2\right)$) is shown in Equation (11):

$$D_{t_1,m_1,n_1,t_2,m_2,n_2}=\sqrt{{{\displaystyle \sum }}_{p=7}^{14}{\left(R_{t_1,m_1,n_1,p}-R_{t_2,m_2,n_2,p}\right)}^2},$$

(11)

For a given radius ($E_0$) and minimum integer ($M_0$), a point $P_c\left(t_c,m_c,n_c\right)$ is a core point if there are more than $M_0$ points within $E_0$ around it. It is denoted as $f_c\left(t_c,m_c,n_c\right)=1$, as Equations (12) and (13) show [40]:

$$f_c\left(t_c,m_c,n_c\right)=\{\begin{array}{c}1,if {\displaystyle \sum }_{m=1}^M{\displaystyle \sum }_{n=1}^N{\displaystyle \sum }_{t=1}^{T}h(D_{t_c,m_c,n_c,t,m,n})>M_0\\ 0,otherwise\end{array},$$

(12)

$$h\left(x\right)=\{\begin{array}{c}1,if x\le E_0\\ 0,otherwise\end{array}.$$

(13)

A point $P_b\left(t_b,m_b,n_b\right)$ is a border point if it is not a core point and there is any core point within $E_0$ around it. It is denoted as $f_b\left(t_b,m_b,n_b\right)=1$, as Equation (14) shows [40]:

$$\begin{array}{c}\hfill f_b\left(t_b,m_b,n_b\right)=\hfill \\ \hfill \{\begin{array}{c}1,if f_c\left(t_b,m_b,n_b\right)=0 and {\displaystyle \sum }_{m=1}^M{\displaystyle \sum }_{n=1}^N{\displaystyle \sum }_{t=1}^T{f}_c(t,m,n)\times h(D_{t_b,m_b,n_b,t,m,n})>0\\ 0,otherwise\end{array}.\hfill \end{array}$$

(14)

Other points $P_n\left(t_n,m_n,n_n\right)$ are noise points. In this paper, the noise points are regarded as abnormal points. It is denoted as $f_n\left(t_n,m_n,n_n\right)=1$, as Equation (15) shows [40]:

$$f_n\left(t_n,m_n,n_n\right)=\{\begin{array}{c}1,if f_c\left(t_n,m_n,n_n\right)=0 and f_b\left(t_n,m_n,n_n\right)=0\\ 0,otherwise\end{array}.$$

(15)

2.3. Statistical Method

An anomaly detection algorithm needs to be validated statistically. We assume that the pre-seismic infrared anomalies commonly occurred near the epicenter. For the given predicted radius ($R_0$) and predicted time window ($W_0$), the region within $R_0$ around the anomaly and the time within $W_0$ after the anomaly determine the predicted range. The hit rate (HR) is defined in Equation (16) [41]:

$$HR=\frac{N_P}{N_E},$$

(16)

where $N_P$ is the number of earthquakes within the predicted range, and $N_E$ is the number of all earthquakes within the study region. The two points ($P_1\left(t_1,m_1,n_1\right)$ and $P_2\left(t_2,m_2,n_2\right)$) belong to an anomaly sample, if they follow the condition of Equations (17)–(21):

$$f_n\left(t_1,m_1,n_1\right)=1,$$

(17)

$$f_n\left(t_2,m_2,n_2\right)=1,$$

(18)

$$\left|t_1-t_2\right|\le 1,$$

(19)

$$\left|m_1-m_2\right|\le 1,$$

(20)

$$\left|n_1-n_2\right|\le 1,$$

(21)

The correlation rate (CR) is defined in Equation (22) [42]:

$$CR=\frac{N_Y}{N_A},$$

(22)

where $N_Y$ is the number of anomaly samples related to earthquakes, and $N_A$ is the number of all anomaly samples. The probability gain (Gain) is defined in Equation (23) [43]:

$$Gain=\frac{HR}{\rho },$$

(23)

where ρ is the spatiotemporal occupancy of the predicted range.

3. Results

3.1. The Deviation of the LSTM Model

Both the pre-processed and predicted data are four-dimensional data. The data at 99.13°E, 29.68°N are shown in Figure 3. The pre-processed data are complicated and varies irregularly. The brightness temperatures in the 9th, 10th, and 14th channels are lower than those in the other channels because of the absorption of water vapors in the atmosphere. The deviation data reveal the difference between the observed and predicted data. As Figure 4 shows, it fluctuates around zero. The blue line is the deviation for the training dataset, and the orange line is the deviation for the test dataset. The images on a randomly selected date (30 January 2021) are used as examples to compare the original data with the deviation images. The images of the brightness temperature on 30 January 2021 are shown in Figure 5, and the deviation images are shown in Figure 6, where the digits denote the terrain in kilometers. The geospatial distribution of brightness temperature is related to the latitude and the topography. There were lower temperatures in the region with a high latitude and a high altitude. The spatiotemporal difference near the absorption band of water vapor was relatively less pronounced. The deviation images show the independence of the latitude. This indicates that the LSTM model could eliminate the background field of infrared brightness temperature, especially the difference caused by the latitude.

The training and testing scores of the LSTM model are shown in

Table 3. The training and testing scores of the LSTM model.

Channel	$\mathbf{Training} \mathbf{Score} \left(\mathsf{\mu }\mathbf{m}\right)$	$\mathbf{Testing} \mathbf{Score} \left(\mathsf{\mu }\mathbf{m}\right)$
7	7.4826	7.8116
8	7.2683	7.3682
9	4.3431	4.4512
10	4.9965	5.0628
11	7.7708	7.9001
12	8.4068	8.5414
13	8.624	8.7911
14	5.1109	5.4133

. A smaller score, which means a smaller prediction deviation, indicates a better performance of the model. The RMSEs in the channels around the water vapor absorption band (channel 9, 10, and 14) are lower. Among the other channels, the RMSEs in the long-wave infrared channels (channels 11, 12, and 13) are slightly higher than those in the medium-wave infrared channels (channels 7 and 8). The frequency distributions of the deviation values are shown in Figure 7. Most of the deviations are around zero.

3.2. Clustering Results

The clustering results of the DBSCAN model depend on its two parameters ($E_0$ and $M_0$). The abnormal rates for different parameters are shown in

Table 4. The abnormal rate for different radii and minimum integers.

$\mathbf{Radius} {\mathit{E}}_0$	$\mathbf{Minimum} \mathbf{Integer} {\mathit{M}}_0$	Abnormal Rate (%)
1	10	14.17
1	102	29.23
1	103	82.01
2	10	2.54
2	102	4.02
2	103	6.72
3	10	1.12
3	102	1.89
3	103	2.24

. These are the proportions of abnormal points in the whole study area. For the same radius, the larger the minimum integer is, the larger the abnormal rate is. A larger radius increases the amount of computation. Considering the algorithm efficiency and anomaly rate, the radius of 2 and the minimum integer of 100 are chosen in this study. For the chosen parameters ($E_0=2, M_0=100$), the deviation data were classified into two classes. About 4% of all the points were abnormal, and the others were normal.

All anomalies within one month before the Madoi and Yangbi earthquakes in the study area are recorded in

Table 5. The information on anomalies from 20 April 2021 to 20 May 2021 (UTC+8).

Start Time (Day)	Duration (Day)	Distance 1 (km)	Distance 2 (km)	Coverage Area (Pixel)
−29	2	1218.514	522.3128	335
−27	2	1614.612	1082.72	76
−27	3	694.1449	147.2815	233
−27	5	357.3739	1338.696	1310
−22	2	256.0721	989.042	174
−21	4	416.9308	308.354	777
−19	2	1818.656	871.9	65
−19	2	774.8023	1100.769	295
−17	4	11.84356	769.3796	809
−16	2	633.2236	515.8192	73
−14	2	622.0482	1264.55	138
−13	6	420.0436	388.9196	1157
−12	2	1438.045	704.7843	40
−11	4	297.1712	573.796	1261
−10	3	415.1837	1417.791	48
−10	2	629.9462	1580.574	85
−8	2	1046.354	79.24198	545
−6	2	682.544	1520.562	160
−5	5	412.239	1011.247	562
−4	2	557.2012	554.8826	349

. The start time is the number of days relative to the Yangbi Earthquake time (21 May 2021, UTC+8). The duration is the number of days from the anomaly starting to the anomaly ending. Distance 1 is the distance between the anomaly and Yangbi epicenter, and Distance 2 is the distance between the anomaly and Madoi epicenter. The coverage area is the number of pixels covered by the anomaly. The longest anomaly occurred from 8 to 13 May 2021. It is about 420 km away from the Yangbi epicenter and 389 km away from the Madoi epicenter. The anomaly images are shown in Figure 8. The black points represent the abnormal points marked by the DBSCAN model. The two red stars represent the locations of the two mainshocks. The coverage area is the largest on 12 May 2021, which is 1157 pixels.

Figure 9 and Figure 10 show the anomaly information of the Yangbi and Madoi earthquakes, respectively. The size of the scatter denotes the anomaly area (S), while the color denotes the anomaly duration. The coordinates denote the spatiotemporal interval between the anomaly and earthquake. The red dotted line is the radius of the seismogenic region calculated according to the so-called Dobrovolsky area, defined as Equation (24) [44]:

$$R={10}^{0.43M},$$

(24)

where M is the magnitude of the earthquake. The large-area and long-duration anomalies are within the seismogenic zone (the region below the red dotted line in Figure 9 and Figure 10).

3.3. Statistical Results

For the given predicted time window ($W_0$) of 14 days,

Table 6. The correlation rate, the hit rate, and the probability gain with different parameters.

Predicted Radius (km)	Correlation Rate	Hit Rate	Probability Gain
1000	0.6429	0.6875	1.6193
800	0.6429	0.6875	1.7358
600	0.6429	0.6875	1.9137
400	0.5714	0.5625	1.7971
200	0.3571	0.375	1.4576

shows the correlation rate, the hit rate, and the probability gain with different predicted radii ($R_0$). The anomaly samples with a duration of over 5 days and a coverage area over 1000 pixels are used to determine the predicted range. There are 14 anomaly samples and 20 earthquakes with a magnitude over 5 in the study area, including 16 mainshocks, 1 foreshock, and 3 aftershocks. The magnitudes of the foreshock and aftershocks are no more than 5.6, as Table 2 shows. The foreshocks and aftershocks are excluded because they occupy a similar spatiotemporal area to the main shock. The highest probability gain is 1.9137, which is calculated with the predicted radius of 600 km. This means that multichannel infrared data could reduce the uncertainty of earthquake prediction. Meanwhile, the correlation is 64.29%, and the hit rate is 68.75%. This means that most of the anomaly samples are related to the earthquakes, and the method we proposed could detect the anomalies from multichannel infrared images before most of the earthquakes. The Molchan diagram is shown in Figure 11. The horizontal axis is the spatiotemporal occupancy of the predicted range, and the vertical axis is the miss rate (1-HR) [45]. The slope of the line through the star in Figure 12 is negative and is equal to the negative of the probability gain [46]. There are both the lowest miss rate and highest probability gain when the predicted radius is 600 km.

In this study, the anomalies were detected within 14 days before 11 earthquakes within 600 km around the epicenters. The information on these 11 earthquakes is listed in

Table 7. The considered earthquakes with a pre-seismic anomaly.

Time (UTC+8)	Magnitude	Latitude (°N)	Longitude (°E)	Depth (km)
16 September 2021, 04:33:31	6	29.2	105.34	10
26 August 2021, 07:38:18	5.5	38.88	95.5	15
13 August 2021, 12:21:35	5.8	34.58	97.54	8
29 July 2021, 16:39:27	5.7	22.7	96.04	20
7 July 2021, 14:43:48	5.2	19.65	101.2	10
16 June 2021, 16:48:58	5.8	38.14	93.81	10
12 June 2021, 18:00:46	5	24.96	97.89	16
10 June 2021, 19:46:07	5.1	24.34	101.91	8
22 May 2021, 02:04:11	7.4	34.59	98.34	17
21 May 2021, 21:48:34	6.4	25.67	99.87	8
19 March 2021, 14:11:26	6.1	31.94	92.74	10

. The anomaly times and the distances between the anomaly and the epicenter are shown in Figure 12. Two earthquakes occurred after the anomaly disappeared. Two earthquakes occurred on the same day when the anomaly disappeared. In the other seven cases, the anomaly continued for several days after the earthquake.

4. Discussion

In this paper, the one-year infrared anomalies are detected and their correlations with the earthquakes are analyzed statistically. The maximum probability gain was 1.91, which was similar to the application result of the total electron content (TEC) anomaly [47]. Probability gains are influenced by the data used, the algorithms, and the regions studied. Filizzola et al. [48] studied the thermal anomalies in Turkey using robust satellite techniques, the probability gain for earthquakes with a magnitude over 5 was 2.2 in their study. A probability gain greater than one means that the prediction method is better than a random guess. Jiao et al. [48] compared true earthquakes and random synthetic earthquakes to verify that the recognition ability of their method was higher than the random guess. At the same time, they proposed that it is expected to improve earthquake prediction performance by combining two–five kinds of anomalies.

As Figure 13 shows, the anomaly occurred in the southeast of the Madoi epicenter from 24 to 26 April and was 147 km away from the epicenter. It covered 233 pixels. The red stars indicate the locations of the epicenters, while the green points indicate the locations of the TEC anomaly. In our previous studies [47], there were TEC anomalies near the region on 23 April 2021. The abnormal region of TEC on 23 April 2021 overlaps with that of infrared data on 25 April 2021. In addition, TEC anomalies were observed near the two epicenters on 18 May 2021, and the infrared anomaly lasted from 17 May 2021 to 18 May 2021, as Figure 14 shows. During this period, the two types of anomalies were closest in location on 18 May 2021. This indicates that there is a certain spatial and temporal correlation between infrared and TEC anomalies before earthquakes.

Although the correlation rate of 64.29% and the hit rate of 68.75% were higher than that of the TEC anomalies, the probability gain was not improved. It may be that the spatiotemporal occupancy of the prediction range was too high because the abnormal coverage area in this paper is large, the prediction area radius is 600 km, and the prediction time window is 14 days. It is hard to pinpoint the exact location of an earthquake. Future research needs to combine multiple data sources to enhance the accuracy of earthquake prediction. Since both LSTM and DBSCAN models are suitable for multi-dimensional feature vectors, the method proposed in this paper can be used for pre-earthquake anomaly extraction from multi-source data, but normalizations of the spatiotemporal scale and amplitude are necessary for data pre-processing.

5. Conclusions

A new algorithm combining the LSTM and DBSCAN models was proposed to extract pre-seismic anomalies from the multichannel data of the Fengyun-4A satellite. The probability gain of 1.91 shows that this method is better than a random guess and can reduce the uncertainty of earthquake prediction. The correlation rate is 64.29%, and the hit rate is 68.75%. In this seismic case analysis, infrared and TEC anomalies exhibited a spatiotemporal correlation. This indicates that the method could be used to extract pre-seismic anomalies from multi-source data in future research for better prediction accuracy.