Detecting Thermal Anomalies of Earthquake Process Within Outgoing Longwave Radiation Using Time Series Forecasting Models

The catastrophic damages caused by the Jiuzhaigou earthquake in China of August 8, 2017 and the Mexico earthquake of September 20, 2017 have revealed some important weaknesses of currently operational earthquake-monitoring and forecasting systems. In this work, six time series forecasting models were applied to detect pre-earthquake anomalies within infrared outgoing longwave radiation. After comparing their prediction results using non-seismic time series data, the autoregressive integrated moving average (ARIMA) model was selected as the optimal model, and then a new prediction method based on this ARIMA model was proposed. The results show that the values observed on July 27 and August 5 before the Jiuzhaigou earthquake in China exceed the confidence interval for prediction and reaches the maximum on August 5, 2017. This indicates the infrared outgoing longwave radiation (IR-OLR) anomalies before the Jiuzhaigou earthquake in China. For the Mexico earthquake, pre-earthquake IR-OLR anomalies are detected on September 14, 18, and 19, and reaches the maximum on September 14, 2017. This demonstrates that the proposed time series forecasting model based on ARIMA could be an effective method for earthquake anomalies detection within infrared outgoing longwave radiation.

research in seismology, especially in earthquake prediction, this satellite makes it possible to solve the current major problems in earthquake prediction research and opens up new approaches to earthquake monitoring and prediction [Cheng et al., 2018].
A large number of earthquake cases have shown that there is a "thermal precursor" to earthquakes [Dai et al., 2016]. Although there is currently no consensus on the mechanism behind seismic infrared anomalies, it is widely believed that a large range of thermal anomalies could occur in the period preceding an earthquake. IR-OLR is the energy emitted from the Earth to space in the form of thermal radiation after the Earth's surface absorbs solar radiation. This raises two questions: (1) whether IR-OLR anomalies can be used as a precursor in earthquake monitoring, and (2) whether they are related to other factors [Prakash and Srivastava, 2017]. Based on the physical causes of IR-OLR pre-earthquake anomalies, Sun discussed the relationship between IR-OLR and earthquakes as early as 1994 [Sun, 1994], and introduced IR-OLR into the field of earthquake prediction for the first time. In recent decades, Chinese researchers have used the IR-OLR data provided by the National Oceanic and Atmospheres Administration of the United States (NOAA) to study the relationship between IR-OLR and earthquakes. They have identified anomalous IR-OLR changes preceding some typical earthquakes that occurred in China and other parts of the globe, such as the Ms 7.8 Tangshan earthquake in 1976, Ms 8.7 Sumatra earthquake in 2004, Ms 8.0 Wenchuan earthquake in 2008, and Ms 8.1 Nepal earthquake in 2015 [Liu et al., 1997;Liu and Kang, 2005;Kang et al., 2009;Zhou et al., 2016].
On this basis, the researchers put forward different methods for detecting IR-OLR anomalies, such as deviation from average, vorticity analysis, wavelet analysis, power spectrum analysis, and statistical analysis, and studied the mechanism responsible Dai et al., 2016;Qin, 2012;Ouzounov et al., 2007;Xiong and Shen, 2016;LU et al., 2017]. Research findings suggest that IR-OLR has a great potential for application in earthquake prediction.
Continuous advancements in scientific research have driven significant breakthroughs in earthquake forecasting.
A major manifestation is that a large number of forecasting methods have been applied to various fields and constantly improved in accuracy. Moreover, since the advent of the era of big data, data have gradually become a hot topic for research in various fields. Among them, time series is type of data with a very wide range of applications [Zhao and Chen, 2016]. Constrained by the physical periodicity of activity of the Sun-to-Earth, the IR-OLR emitted by the Earth exhibits periodic variations, especially seasonal variations [Jing et al., 2016;Sedlar et al., 2016], allowing IR-OLR to be used as a time series. Thanks to technological progress, people have gradually begun to apply the classic ARIMA, neural network and other models to preliminary predictions of earthquake magnitude, time, and location by means of time series forecasting [Fu, 2010]. Scholars have tried many methods to predict earthquakes in their research. In many domains, predictions are made based on a single model, such as ARIMA and BP neural network. For instance, Hikichi et al. [2017] used an ARIMA model to predict the growth in the number of ISO 14001 certifications in the United States and Zheng [2017] used an improved BP neural network model to predict the water quality of a river section in central Liaoning.
In this study, we discussed the method and theory for IR-OLR time series analysis and analyzed the predictive ability of different time series forecasting models. The 1-degree IR-OLR data (IR-OLR) provided by HIRS/4, an advanced sensing instrument operational on the NOAA-18, was used as the time series. Earthquake predictions were made using six prediction models, including an ARIMA model. Furthermore, the time series forecasting models were tested for validity, with the root mean square error (RMSE) of their predictions being used as a measure of accuracy. The results show that the pre-seismic IR-OLR anomalies predicted using the time series analysis based on the ARIMA model are obviously more accurate than the predictions from the Neural Network Auto-Regressive (NNAR), Seasonal decomposition of time series by loess-Exponential smoothing (STL-ETS), exponential smoothing state space model with Box-Cox transformation, ARMA errors, trend and seasonal component (TBATS), Naive, and Holt-Winters Models.
Therefore, the ARIMA model was used to predict the pre-seismic IR-OLR anomalies before the Jiuzhaigou in China and the Mexico earthquakes by time series analysis. The experimental results are satisfactory, suggesting that the IR-OLR time series analysis based on the ARIMA forecasting model built in this study has high predictive ability and considerable potential for wide application.
The paper is organized into six sections. Section 2 outlines the composition of HIRS/4 and the HIRS/4 OLR data products used in this paper. Section 3 describes the principles and process involved in the time series forecasting models and the procedure of data preprocessing for IR-OLR time series analysis followed by an evaluation of these models. Section 4 provides the results of the comparative analysis of the six time series forecasting models built and the results of pre-seismic IR-OLR time series analysis based on the ARIMA forecasting model. Section 5 discusses the findings of the study and section 6 presents the conclusions and expectations about future work in this domain. The HIRS/4 uses a single telescope and a rotating filter wheel containing 20 individual spectral filters to provide multi-spectral data from one visible channel (0.69 microns), seven shortwave channels (3.7 to 4.6 microns), and twelve longwave channels (6.7 to 15 microns). Each channel has an instantaneous field of view (IFOV) of approximately 7.0 degrees, which is from the spacecraft altitude of 870 km (470 micron) and covers a circular area of 10 km (6.2 microns) diameter at nadir on Earth. This is an improvement in resolution over the 20km (12.4 microns) HIRS/3 instrument that was flown on NOAA-KLM. An elliptical scan mirror provides a cross-track scan in 56 steps of 1.8 degrees per step. The mirror steps then rapidly stop and remain at their positions while the optical radiation that has passed through 20 spectral filters is sampled. Each Earth scan takes 6.4 seconds and covers 49.5 degrees from nadir. IR calibration of the HIRS/4 is achieved by views of space and the internal warm target, each viewed once every 38 Earth scans. The instrument measures scene radiance in the IR spectrum. The data is also used to determine ocean surface temperatures, total atmospheric ozone, precipitable water, cloud height and coverage, and surface radiance [Wang et al., 2006;Kelly, 2007]. Turner and Tett [2014] have indicated that the HIRS OLR product was developed in the late 1980s by Ellingson et al. [1994] and now provides a continuous data record spanning more than 3 decades [Lee et al., 2007]. It uses a multispectral regression technique to estimate OLR from radiances measured by 4 HIRS channels, including channels 3, 7, 10, and 12. After a change to HIRS/2 channel 10's spectral response function in the later HIRS/4,

Calculation of non-seismic data values using time series forecasting models
Time series analysis is a statistical analysis method for analyzing and exploring a set of time-ordered data in order to predict its future trend. Variations in IR-OLR are strongly associated with geographical location, season, total electron content (TEC), and the activities taking place in the interior of the Sun and the Earth. These factors contribute to the obvious periodicity, especially seasonality, in IR-OLR [Yun and Jie, 2005]. Application of time series analysis to this series allows us not only to take its deterministic components into account, but also to calculate the uncertain components of IR-OLR changes, so as to get more accurate non-seismic data value of the IR-OLR. The process of modeling and calculating the non-seismic data value of IR-OLR is shown in Figure 3. It involves four major steps. The first step is to perform a stationary test (difference methods) to stabilize the IR-OLR series. Then a periodic model of the stationary series is constructed. The next step is to select the optimal model by using the RMSE accuracy criterion. The optimal model is then utilized to predict non-seismic data values of IR-OLR. The specific steps are as follows: Firstly, a differential operation is applied to the IR-OLR series, with the cycle length s being used as the step length. The cycle length (s) is distinguished based on the observation timing chart and the autocorrelation function of the observation series, and the periodic difference is given by formula (1): Earthquake anomalies detecting within OLR In the formula, is a backward shift operator and -= , ∇ is a difference operator. If a trend still exists in the time series after the periodic difference, the normal difference ∇∇ is then applied to the series (t): After the difference, the time series forecasting models are used to fit the stationary series. In general, the original series can be transformed into a stationary time series by a periodic difference and a normal difference. If not, continue the difference until it becomes stationary.
After modeling is completed, a significance test of the models is carried out. Then, the models that pass the significance test are evaluated based on RMSE in order to select the optimal model. The optimal model is used to build model and predict the IR-OLR non-seismic series, so as to obtain the non-seismic data values of the IR-OLR to be predicted.

Building time series-based forecasting models
In addition to the forecast non-seismic data values, the prediction of the pre-earthquake IR-OLR anomalies also requires a reasonable error range for non-seismic data values. In this paper, a more reasonable bound determination strategy is proposed for determining the upper and lower bounds of the non-seismic data value range in the time series analysis based on the time series forecasting models. The specific steps are as follows: Dulin Zhai et al.
6 The first step is to select the IR-OLR data that are closest to the epicenter and more than one year before earthquake, which can ensure the absence of seismic disturbance. The data cover the period from the 15 th to the 27 th months before the earthquakes and are spaced at one-month intervals were used as training datasets, and the 5-days spaced data for the 13 th months before earthquakes was used as the test set, which starts at the day one year before the earthquake and terminates at the 30 th day backwards.
The reasons for choosing this temporal window in this paper are as follows: (1) Dai et al. [2009] used IR-OLR to extract the earthquake-related anomaly information from 15 different earthquake cases based on wavelet packet method. The results show that different degrees of IR-OLR earthquake anomaly information are extracted during 10 months before the shocks, and there is a significant correlation between the duration of energy attenuation and earthquake magnitude. Meanwhile, Jing et al. [2009] used IR-OLR daily grid data (2.5°) to analyze the mid-term infrared emission anomalies related to medium-term and large earthquakes in Sichuan-Yunnan area using wavelet time-frequency analysis method. It was found that the anomalies are mainly occurred on half a year to one year before earthquakes. Therefore, considering that both the Jiuzhaigou earthquake in China and the Mexico earthquake are large earthquakes, and the magnitude of the earthquake is relatively high. In this paper, the period from the 15 th to the 27 th months before the earthquakes and are spaced at one-month intervals were used as training datasets, which can effectively avoid a certain degree of interference caused by the anomaly of IR-OLR before the earthquake.
(2) In recent decades, a large number of experts and scholars have studied IR-OLR anomalies such as Lu et al.
[2014]and Liu et al. [1997] The results show that the IR-OLR anomaly occurs within a few days to one month before the earthquake. Therefore, this paper chooses the 5-days spaced data for the 13 th months before earthquakes was used as the test set, which starts at the day one year before the earthquake and terminates at the 30 th day backwards, which can improve the detection results of IR-OLR anomalies.
(3) Before using RST method to analyze TIR at the time of the Abruzzo 6 April 2009 earthquake,  and Pergola et al. [2010] verified that the absence of TIR anomalies in a relatively seismically unperturbed period, the confutation phases has been performed by considering the same period (15 March-15 April) but in a different year (2008). Therefore, in this paper, the same approach is adopted, that is, in order to verify that the absence of IR-OLR anomalies in a relatively seismically unperturbed period, the confutation phases has been performed by considering the same period (the 5-days spaced data for the 13 th months before earthquakes was used as the test set, which starts at the day one year before the earthquake and terminates at the 30 th day backwards) but in a different year. By predicting each modeling forecasting models with this temporal window, and the accuracy of each building time series-based forecasting models is obtained successively, and then the optimal time series-based forecasting models is selected to improve the accuracy of the pre-earthquake IR-OLR anomaly detection.
Finally, the differences between the non-seismic data values for this period and the actual observations were calculated. The residual value with 95% confidence level, Δ, was determined by statistical analysis. Then the upper and lower bounds of the non-seismic data range can be calculated from the residual value determined and the non-seismic data value for the normal period using the equations below: (3) In the formula, is the non-seismic data value predicted by the time series analysis based on a time series forecasting model; Δ is the residual value with 95% confidence level for the predicted non-seismic IR-OLR data (without disturbance) that is closest to the epicentre; and are respectively the upper and lower bounds of the non-seismic data range, the bounds beyond which is considered as anomaly. The overall IR-OLR anomaly prediction process is illustrated in Figure 4.

Data Preprocessing
To ensure the reliability of the model, the IR-OLR data is used to build model is more than one year before the selected earthquakes, which were not affected by seismic disturbance. At the same time, in order to verify the proposed model, the earthquakes with magnitude of 7.0 or higher are selected. Therefore, the Ms 7.0 Jiuzhaigou earthquake in China on August 8, 2017 and the Ms 7.1 Mexico earthquake on September 20, 2017 were selected as the case study.
In accordance with the process steps shown in Figure 1 and 2, the night time IR-OLR data points that were closest to the epicenter and more than one year before the Jiuzhaigou in China and Mexico earthquakes were selected for modeling. This ensures that the selected data set has not been affected by seismic disturbance. For more reliable predictive modeling, the selected IR-OLR series were subjected to stationarity and white noise tests. The unit root test was used for the stationarity test in the study. Later, the results of the stationarity test using the IR-OLR series of the Jiuzhaigou in China and Mexico earthquakes were obtained.
The stationarity test results show that the p-value of the unit root of the IR-OLR series for the Jiuzhaigou earthquake in China is less than 0.05. This indicates that the IR-OLR series is stationary. However, the p-value of the IR-OLR series for the Mexico earthquake being modeled exceeds 0.05, indicating that this series is nonstationary.
Therefore, it is necessary to make the IR-OLR series for the Mexico earthquake stationary by the first-order difference.
After the first-order difference of the IR-OLR series for the Mexico earthquake, the series is uniformly distributed around 0, indicating general stationarity. The later unit root test showed that the p-value is less than 0.05. This confirms that the IR-OLR series for the Mexico earthquake becomes stationary. After the stationarity test, the pvalue of the results of white noise test is greater than 0.05 for both earthquakes, indicating that the residual error of the fitted model can be detected as white noise.

Time Series Forecasting Models
This study analyzed the prediction principles and process involved in the time series forecasting models, including ARIMA, NNAR, STL-ETS, TBATS, Naïve and Holt-Winters.
The ARIMA model, also known as the Box-Jenkins Model or the B-J model, is a method of time series forecasting [Wang and Gu, 2000]. It has the following form: where represents time series data, is related to -1 (i = 1,2, ... , p); denotes the residual term, which is related to -(j = 1,2, ... , q); is a delay operator, and satisfies; = -; p represents the autoregressive order; q is the moving average order; is the difference order; ∇ is the difference operator, ∇ = (1-); Φ( ) represents the polynomial of the autoregressive coefficient; Θ( ) represents the polynomial of moving average coefficient; and is a white noise series that is independent ofand -.
The model that satisfies equation (4) is ARIMA (p, d, q), an autoregressive integrated moving average model which combines the ARMA (p, q) model and the differential operation organically. It is able to provide high-precision short-term predictions and does not require highly structured data. To provide a good fit, this model only needs to find a regular pattern in the data itself [Mini et al., 2015]. However, as the IR-OLR series is univariate and contains a seasonal component, this study used the ARIMA (p, d, q) (P, D, Q) m model proposed by Hyndman and Khandakar [2008] to obtain the best model parameter values through automatic predictions.
The neural network auto-regressive (NNAR) model can be used to model the complex nonlinear relationships between input and output variables. In the case of the NNAR model, the lagged value of time series is used as the model input, and the output is the predicted value of the time series. Due to the presence of seasonal components Φ ( )∇ = Θ ( ) ( ) = 0, ( ) = 2 , ( ) = 0, ≠ ( ) = 0, ⩝ < in the IR-OLR series, this study used the NNAR (p, P, k) m model proposed by Hyndman and Athanasopoulos [2014].
Single seasonal exponential smoothing methods are among the most widely used forecasting procedures in practice [Snyder et al., 2002;Makridakis et al., 2010]. Seasonal decomposition of time series by loess-exponential smoothing (STL-ETS) model is a time series decomposition method based on robust locally weighted regression. Before fitting the STL-ETS model to this time series, smoothing was first performed to explore its overall trend and decomposes it to see if the time series is affected by seasonal factors. In this study, the seasonal IR-OLR series data were decomposed into trend factors, seasonal factors and random factors [Cleveland et al., 1990;.
These models require estimation of the initial seasonal values of 2(k 1 + k 2 ++ k T ), which is expected to greatly reduce the number of parameters required compared to BATS models. The use of trigonometric functions also allows modeling of non-integer seasonal frequencies [De Livera et al., 2009]. Additionally, due to presence of a seasonal component in the IR-OLR series, the value predicted by the Naive model is equal to the last value of the same season ( +ℎ| = +ℎ-, where m is equal to the length of the seasonal period, and k = +1). In general, this model has better prediction performance for seasonal time series [Makridakis and Hibon 2004;Barot, 2004]. The Holt-Winters model is an improvement and development of the moving average method for exponential smoothing. Therefore, the Holt-Winters' three-parameter exponential smoothing model is essentially a high-level exponential smoothing model capable of simultaneously dealing with trend and seasonal variations. Moreover, it can properly smooth out random fluctuations and predict future values of time series that exhibits both a long-term trend and a seasonal pattern [Hyndman and Athanasopoulos, 2018]. In view of the fact that the patterns in the data analyzed in this study are consistent with the additive model, an additive model was selected for seasonality and trend and the random terms were smoothed [Mini et al., 2015;Elmunim et al., 2015].

Construction and evaluation of time series forecasting models
In order to make the determination of the parameters in the aforementioned prediction model more convenient and more optimal, we utilized existing time series forecasting models in combination with a powerful and widely applicable automatic prediction algorithm to simplify the process of determining specific parameters [Hyndman et al., 2000]. The steps involved are summarized below.
First, the time series forecasting models were applied to the appropriate IR-OLR data, and the model parameters were optimized to best fit the data in each case. Then, the optimal model was selected according to the RMSE.
Forecasts were later produced by the optimal model with optimized parameters. By assuming a standard normal distribution of the non-seismic data covering 13 months, we obtained the prediction intervals and then estimated the 95% confidence intervals for all forecasting models.
Before the comparative analysis based on RMSE, modeling and prediction processing were performed using time series cross-validation. The specific method and principle are as follows: First, after verifying the stability and periodicity of the series, prediction was made by using IR-OLR data unaffected by seismic disturbance as training datasets. The data cover the period from the 15 th to the 27 th months before the earthquakes and are spaced at one-month intervals were used as training datasets.
Next, the 5-days spaced data for the 13 th months before each earthquake was used as the test set, which starts at the day one year before the earthquake and terminates at the 30 th day backwards. represents the number of predicting days covered by the test set more than one year before the earthquake, and "Training-Months" represents the number of months covered by the training sets more than one year before the earthquake. For example, Training = 3 means the training set covers the 3 months from the 13 th to the 15 th month before the earthquake). As described above, in order to identify the superior performance of each prediction model and validate the models' predictions, RMSE was used as a performance measure to analyze and compare the accuracy of these models.

RMSE
RMSE can be calculated using the equation below [Mini et al., 2015]: where represents the difference between the observed value and the fitted value of the time series at time t.
The smaller the RMSE value, the higher the accuracy of the prediction. Conversely, a higher RMSE indicates lower prediction accuracy.

Analysis of time series forecast results
The time series forecasting models described above were applied to the N-OLR series used as the training sets (data covering from the 15 th to the 27 th month before the Jiuzhaigou in China and Mexico earthquakes). Then the total RMSE of the models' predictions for different training set and test set intervals were obtained for the purpose of comparative analysis based on the aforementioned time series cross-validation method and RMSE. More prediction results from the models are detailed form Table 1 to Table 12 in  As shown in Figure 5 and Figure 6, for both earthquakes and across all time intervals, the traditional ARIMA model's predictions have the lowest total RMSE compared to other five models' predictions. This suggests that the traditional ARIMA method is superior to other prediction models in terms of prediction accuracy.
Based on the total RMSE curves for the models, the traditional ARIMA model was selected as the optimal model for predictive modeling of the Jiuzhaigou in China and Mexico earthquakes. To verify this model, further forecasting was performed using the same principle. During verification, the IR-OLR data spaced at one-month intervals for the period from the 3 rd to the 15 th month before the earthquakes were used as the training sets.
Then the data spaced at 5-days intervals for the month preceding each earthquake was used as the test set, which starts at the day before earthquake and terminates at the 30th day backwards. Then the optimal model selected was used to predict and the corresponding predictions were obtained. The total RMSE of the ARIMA model's predictions for different training set and test set intervals are provided in Tables 3 and 4. Table 3 shows that the prediction for the 20 days before the Jiuzhaigou earthquake in China using 12 months series data before the earthquake as the training set has the lowest RMSE, at 34.55936. This suggests that the model prediction for the Jiuzhaigou earthquake in China has the highest prediction accuracy when Training time = 12 months and Test time = 20 days. For the Mexico earthquake (Table 4), the prediction for the 30 days before the earthquake made 12 months before the earthquake (Training time = 12 months) has the lowest RMSE, at 50.86853, demonstrating that the model prediction for this earthquake has the optimal prediction accuracy when Training time = 12 months and Test time = 30 days. Therefore, the selected Training and Test sets that can ensure optimal prediction accuracy were applied to the prediction using the ARIMA model. The corresponding predicted values and residual value Δ with 95% confidence level were obtained for the two earthquakes.  Table 3. RMSE using the ARIMA forecasting model for IR-OLR series of the Jiuzhaigou earthquake in China ("Testing-Days" represents the number of predicting days covered by the test set before the earthquake, and "Training-Months" represents the number of months covered by the training set before the earthquake. For example, Training = 3 means the training set covers the 3 months from the 3 rd to the 5 th month before the earthquake).

RMSE
Testing-Days  Table 4. RMSE using the ARIMA forecasting model for IR-OLR series of the Mexico earthquake ("Testing-Days" represents the number of predicting days covered by the test set before the earthquake, and "Training-Months" represents the number of months covered by the training set before the earthquake. For example, Training = 3 means the training set covers the 3 months from the 3 rd to the 5 th month before the earthquake).

Analysis of pre-seismic IR-OLR anomalies predicted by the ARIMA model
In order to detect pre-seismic anomalies from the time series of IR-OLR data, a confidence interval was determined by calculating the upper and lower bounds based on the predictions and residual value with 95% confidence level obtained. Finally, the differences between the calculated upper and lower bounds L and their actual values P were obtained as the predictions of IR-OLR anomalies. The prediction results for the Jiuzhaigou in China and Mexico earthquakes are shown in Figure 7 and Figure 8, respectively.
As can be seen from Figure 7 and Figure 8, both the Jiuzhaigou in China and Mexico earthquakes were effectively predicted by the ARIMA model. This demonstrates that IR-OLR data can be used as valid, reasonable and reliable time series in earthquake prediction and has certain practical value for relevant research. Moreover, the observed values for July 27 and August 5 before the Jiuzhaigou earthquake in China exceed the estimated confidence interval for prediction, and the maximum (39.44 W/m 2 ) exceed bounds occurred on August 5, 2017. So, we inferred that IR-OLR anomalies existed on July 27 and August 5, twelve and three days before the Jiuzhaigou earthquake in China respectively, and the intensity of IR-OLR anomalies reached the maximum 3 days before the earthquake on August 5, 2017.
As can be seen from Figure 8, the observed values for September 14, 18, and 19 before the Mexico earthquake are also exceed the predicted confidence interval and the exceed bounds reaches the maximum (50.83W/m 2 ) on September 14, 2017. Therefore, we inferred that IR-OLR anomalies occurred one, two and six days before the Mexico earthquake respectively, and the intensity of IR-OLR anomalies reached the maximum during six days before the earthquake (September 14, 2017).

Dulin Zhai et al.
14 Figure 7. Predictions of IR-OLR anomalies before the Jiuzhaigou earthquake in China.

Discussion
In the past few decades, many researchers, such as Chen et al. [2010], have discussed the relationship between IR-OLR and earthquakes, and introduced IR-OLR into the field of earthquake predictions. They presented the potential relationship between outgoing long-wave radiation anomalies and earthquakes and revealed anomalies prior to major earthquakes. Nevertheless, there has been no consensus on the causes of pre-earthquake IR-OLR anomalies so far. We believe that the near surface air ionization and latent heat exchange due to humidity variations [Pulinets et al., 2006] may cause interactions between the lithosphere-hydrosphere and atmosphere, which is related to changes in the near-surface electrical field and gas composition.
In addition, researchers from China and other countries have also used different processing methods to extract IR-OLR data for pre-seismic IR-OLR anomaly detection. For instance, Xiong et al. [2010] proposed a spatial-temporal analysis method based on eddy field calculation mean and wavelet maxima for detecting seismic anomalies within the outgoing longwave radiation data. Zhou et al. [2016] used a standard deviation threshold to detect earthquake anomalies in outgoing longwave radiation data from NOAA. And Sun et al. [2017] applied RST and anomaly index algorithms to analyze the characteristics of IR-OLR from satellites, to detect seismic anomalies. Though all the methods have produced good results, the results for a given earthquake vary between different methods of the same earthquake. Besides, these methods are only applicable to detection of limited pre-seismic IR-OLR anomalies and can't effectively explain the anomalies.
Considering the limitations of the above-mentioned methods in detecting pre-seismic anomalies, we designed a new prediction method based on time series forecasting models (such as the ARIMA, NNAR, STL-ETS, TBATS, Naive and Holt-Winters models). This method uses N-OLR data as time series. Later, the prediction results produced by the six models based on the time-series cross-validation method were compared against the RMSE, and then an optimal model was selected for earthquake prediction. This method can predict pre-earthquake IR-OLR anomalies while avoid the limitations of the wavelet-based and other detection methods mentioned above. It has great advantages in terms of time, performance and methodology. In terms of time, this method allows us to make effective predictions before an earthquake occurs. In terms of performance, the ARIMA model is a well-established classical prediction model with wide applications in many areas, and it significantly outperform other models analyzed in this study. In terms of methodology, the use of cross-validation and RMSE for accuracy evaluation enables the proposed prediction model to deliver more reliable and valid predictions of pre-earthquake IR-OLR anomalies. Therefore, the ARIMA based prediction model has the potential to be a useful research tool for prediction of pre-earthquake IR-OLR anomalies.

Conclusions
A new method for predicting pre-earthquake infrared radiation anomalies was proposed based on time series forecasting models in this paper. The Ms 7.0 Jiuzhaigou earthquake in China on August 8, 2017 and the Ms 7.1 Mexico earthquake on September 20, 2017 are analyzed using the proposed method. The following conclusions can be reached: (1) Among the six basis time series models, the ARIMA model delivered more valid and reliable predictions than the other five models, as indicated by the accuracy evaluation based on RMSE.
(2) Considering the probability of IR-OLR anomalies, the ARIMA model was used to predict IR-OLR anomalies for the 20 days (Test time = 20 days) preceding the Jiuzhaigou earthquake in China and 30 days (Test time = 30 days) preceding the Mexico earthquake. The series covering the 12 months before the earthquakes were selected as the time series for prediction. Figures 8 and 9 demonstrate that the ARIMA model has good (3) This paper reveals that obvious IR-OLR anomalies did exist the twelve and three days before the Jiuzhaigou earthquake in China, and the intensity of IR-OLR anomalies reached the maximum during three days before the earthquake (August 5, 2017). Moreover, this paper also reveals that obvious IR-OLR anomalies did exist the one, two and six days before the Mexico earthquake, and the intensity of IR-OLR anomalies reached the maximum during six days before the earthquake (September 14, 2017). The consistency between the observations and predictions demonstrates that anomalous IR-OLR changes will probably occur before earthquakes, and it is reasonable to expect that predicting pre-earthquake IR-OLR anomalies by using the time series analysis based on the ARIMA model will be a useful approach to imminent earthquake prediction.
Author Contributions. Dulin Zhai wrote the majority of the paper along with Pan Xiong. Xuemin Zhang contributed to the writing of the paper and interpretation of the results.