Analysis and Prediction of the SARIMA Model for a Time Interval of Earthquakes in the Longmenshan Fault Zone

Abstract

Based on the catalog data of earthquakes with Ms ≥ 2.5 in the Longmenshan fault zone from January 2012 to September 2021, we establish an earthquake time interval series grouped by earthquake magnitude and then use the SARIMA model to predict the series in different periods. By analyzing the fitting effect of the models, the optimal model parameters of different magnitude sequences and the corresponding period values are obtained. Among them, the adjusted R2 values of each model with Ms ≥ 2.5 and Ms ≥ 3.0 sequences are more than 0.86, up to 0.911; the short-time prediction effects are good, and the values of predicted RMSE are 10.686 and 8.800. The prediction results of the models show that the overall trend of the subsequent earthquake time interval in the Longmenshan fault zone is stable, and the prediction results of the Ms ≥ 3.0 sequence have a weak fluctuating growth trend; that is, the number of earthquakes with the Ms ≥ 3.0 in this area will decrease slightly, and the seismicity will decrease in a period of time. The analysis results and method can provide a scientific basis for earthquake risk management and a feasible way to predict earthquake occurrence times.

1. Introduction

The Longmenshan fault zone is located on the eastern edge of the Qinghai-Tibet Plateau and the northwestern edge of the Sichuan Basin. It extends in a NE-SW direction, reaches the Qinling fault zone in the north, and ends at the Xianshuihe-Xiaojiang fault zone in the south. It is approximately 500 km long and 50 km wide and is mainly composed of four main faults, including the Houshan fault, Central fault, Qianshan fault, piedmont hidden fault, and the thrust nappe controlled by them [1, 2]. Dozens of important cities and towns, such as Dujiangyan, are distributed within its scope. In a southeastern direction, this large fault at the front of the main mountain range is close to the densely populated Chengdu Plain area and the Chengdu-Chongqing economic circle. Exploring the regularity of earthquake occurrence times along the Longmenshan fault zone can provide a scientific basis for earthquake management and decision-making for earthquake prediction in this area.

Currently, earthquake research is an interdisciplinary and comprehensive field. For earthquake time series data, many researchers have used statistical and probabilistic methods to predict important information, such as earthquake occurrence time. Some scholars have focused on the exploration of magnitude time series [3, 4, 5]. There have been studies on the time interval sequence of earthquake occurrence, most of which have analyzed the probability distribution of the earthquake time interval [6, 7, 8, 9, 10]. There have also been studies on the prediction of time intervals. Guo & Xu [11] made a confirmatory prediction of the Ms8.0 earthquake in northwestern China using a method of multiplied periods with golden section, which appeared to be in good agreement with the testing results. Mariani and Tweneboah [12] proposed applying the stochastic differential equation based on the superposition of independent Ornstein-Uhlenbeck processes to earthquake research and used this model to fit an earthquake sequence in South America. In addition, some scholars have expanded the prediction approach, using SVM, LSTM and regression algorithms to predict earthquake times [13, 14, 15, 16].

The seasonal autoregressive moving average model can be used not only to predict natural phenomena such as river flow and rainfall [17, 18, 19, 20] but also to predict social and economic phenomena such as transportation passenger flow and trade volume [21, 22, 23, 24, 25]. In the field of earthquake prediction, some scholars have used ARIMA to identify earthquake precursor anomalies [26, 27, 28, 29]. However, few ARIMA or SARIMA models have been used to fit earthquake occurrence time data.

In this paper, we selected earthquake events in the Longmenshan fault zone in Sichuan Province in China since 2012 as the research object. We used the SARIMA model to fit the time interval sequence according to magnitude to predict the development trend of the earthquake occurrence time interval and the time of the next earthquake occurrence.

2. SARIMA model

2.1 Model theory

The ARIMA model was first proposed by Box and Jenkins in the 1970s. Its full name is the differential autoregressive moving average model, denoted as ARIMA (p, d, q), where p and q represent autoregressive and moving average orders, respectively; d is the difference order, and the basic model structure is as follows:

$$\varnothing \left(B\right){\nabla }^d{x}_t=\Theta \left(B\right){\epsilon }_t$$E1

In the formula, $\varnothing \left(\mathrm{B}\right)=1-{\varnothing }_1\mathrm{B}-\dots -{\varnothing }_{\mathrm{p}}{\mathrm{B}}^{\mathrm{p}}{\nabla }^{\mathrm{d}}={\left(1-\mathrm{B}\right)}^{\mathrm{d}}$ is the autoregressive coefficient polynomial of order p, where B is the delay operator, and Bx_t = x_t-1; ${\nabla }^d$ is the d-th order difference; and $\Theta \left(B\right)=1-{\theta }_{1}B-\dots -{\theta }_q{B}^q$ is the moving average coefficient polynomial of order q [30].

ARIMA can model nonstationary time series without seasonal effects. However, in real life, many time series have a certain periodicity. For the series that contains both seasonal effects and long-term trend effects and a complex interaction between them, we can use the SARIMA model.

The general expression of the SARIMA model was originally proposed by Wang et al. (Wang et al., 2018) and is recorded as $\text{ARIMA}\left(\mathrm{p},\mathrm{d},\mathrm{q}\right)\times \left(\mathrm{P},\mathrm{D},\mathrm{Q}\right)\mathrm{s}$, where P and Q are seasonal autoregressive and seasonal moving average orders, respectively; D is the order of seasonal difference; and s is the number of seasonal cycles. As an extension of the ARIMA model, this model extracts seasonal effect information from the series by seasonal difference.

For the time series x_t, the SARIMA model expression is as follows:

$${\nabla }^d{\nabla }_s^D{x}_t=\frac{\theta \left(B\right){\theta }_s\left(B\right)}{\varnothing \left(B\right){\varnothing }_s\left(B\right)}{\epsilon }_t$$E2

where ${\nabla }_s^D$ is the D seasonal difference in s steps; ${\mathrm{\epsilon }}_{\mathrm{t}}$ is the random interference of the error term at time t; ${\varnothing }_{\mathrm{s}}\left(\mathrm{B}\right)$ and ${\mathrm{\theta }}_{\mathrm{s}}\left(\mathrm{B}\right)$ represent the P-order seasonal autoregressive coefficient polynomial and Q-order seasonal moving average coefficient polynomial, respectively. The meanings of other variables are shown above.

2.2 Forecast model construction

For the prediction of small sample data, the traditional time series prediction model can fully extract the information and can also avoid overfitting. After analyzing the data characteristics, this paper chooses to use SARIMA to model the earthquake time series. As a traditional statistical prediction model, the SARIMA model has a good fitting degree to small sample data. It has high prediction accuracy and fast training speed and can reflect the dynamic changes in data. While capturing the trend of the series, the model can also extract its periodic fluctuations [17, 31].

To obtain the time interval sequence, we take the difference from the earthquake sequence. Then, we use SARIMA to model the data that are grouped according to the magnitude and obtain models with different parameters, model-fitting value series ${\mathrm{Y}}_{\mathrm{t}}^{\prime }$ and predicted value series ${\mathrm{Y}}_{\mathrm{t}}^{\prime }$ . To ensure that the model can extract short-term and long-term information and improve the accuracy of prediction in each period, long, medium, and short periods with the best modeling effect are selected, expressed as ${\mathrm{s}}_{\mathrm{a}}$, ${\mathrm{s}}_{\mathrm{b}}$, and ${\mathrm{s}}_{\mathrm{c}}$. The model analysis flowchart is shown in Figure 1.

Figure 1.

SARIMA model analysis flow chart.

3. Study on the time interval of earthquakes in the Longmenshan fault zone

3.1 Data and samples

This paper uses earthquake occurrence time data from January 9, 2012, to September 24, 2021, in the earthquake catalog of the China Seismological Network (http://www.ceic.ac.cn). The earthquake occurrence time, geographic coordinates, magnitude, and other factors of the earthquake event time series are characterized by volatility and nonlinearity. Therefore, we preprocess this information into data in a specific space and a certain magnitude; that is, the time series with an Ms2.4 and above in the longitude and latitude range of the Longmenshan fault zone (29.5° N-33.5° N, 102° E-107° E). By further screening whether the earthquake location is in the fault zone, 437 experimental sample data are obtained, and the earthquake occurrence time is from January 9, 2012, to August 2021, Ms. ∈ (2.5,7).

The earthquake time interval sequence is obtained by the difference operation of the selected samples. The sample has the following characteristics:

1. The sample size is relatively small. The monitoring data in the early years are not sufficiently complete, which can lead to a significant reduction in the accuracy of calculating the time interval values. Therefore, an earthquake sequence of nearly 10 years is selected to effectively avoid this problem, but it also reduces the amount of data.

2. Large time span. The data are from January 2012 to August 2021, including the earthquake events of the Longmenshan fault zone in the past 10 years.

3. Nonstationary and strong volatility. There are many factors affecting the data, there may be complex interactions, and the earthquake sequence has complex periodic characteristics.

In addition, the 58th earthquake recorded in the sequence was the 2013 Ya’an Lushan Ms. 7.0 earthquake, followed by several very close aftershocks, resulting in a short time interval between the 58th and 203th earthquakes in the monitoring data, and the sequence approached zero in a period of time. The original sequence is represented with X_t, and the sequence diagram is shown in Figure 2.

Figure 2.

Sequence diagram of earthquake time intervals.

3.2 Model determination

The characteristics of the data are analyzed, and its trend, seasonality, and residual series are extracted through seasonal decomposition. The results are shown in Figure 3, indicating that the series has a certain trend and seasonality. It is preliminarily speculated that the minimum period is approximately 5 days.

Figure 3.

Trend and seasonal decomposition of earthquake time interval series.

The SARIMA model requires the sequence to meet the precondition of stationarity or post-difference stationarity. After the stationarity test of X_t, we determine the number of differences d = 2. The sequence after the difference is marked as X_d, and its sequence diagram is shown in Figure 4. The ADF test is used to determine whether the data after the difference have a trend. The test principle is to check whether there is a unit root in the sequence: if the sequence is stable, there is no unit root. After testing, the ADF test statistic of sequence X_d is −15.4634, and the p value is less than 0.05, indicating that the sequence is stable.

Figure 4.

Sequence diagram of the earthquake time interval after secondary difference.

The ratio of the model training set and test set of sequence Xd is 9 to 1, and the order of the model is determined according to the ACF and PACF of the data. ACF is used to measure the correlation between the current value of the series and its lagged terms. PACF is used to measure the correlation between the current value and its lagged terms after removing the influence explained by the previous lag. The nonseasonal order of the model is preliminarily determined from the ACF and PACF charts in Figure 5. The PACF diagram can be regarded as a trailing or ninth order truncation.

Figure 5.

ACF and PACF diagrams of the earthquake time interval sequence.

The optimal order of the model in combination with the adjusted R² of the model is selected. The fitting effect of ARIMA models with different orders is shown in Table 1. From this, the parameters p = 9 and q = 1 in the model are determined.

ARIMA model	Adjusted R^2	Parameter test	LB test(p value)
ARIMA (0, 2, 1)	0.643	not significant	0.0001
ARIMA (7, 2, 1)	0.803	pass	0.031
ARIMA (9, 2, 1)	0.813	pass	0.395
ARIMA (10, 2, 1)	0.813	not significant	0.305

Table 1.

Comparison of ARIMA models with different autoregressive orders.

Then, we analyze the ACF and PACF of integral multiple orders of s after a one-order, s-step difference of X. Taking s = 22 as an example, as shown in Figure 6, the ACF and PACF coefficients lagging 22 orders are outside the range of two times the standard deviation and then gradually converge. It is preliminarily assumed that P = 1 and Q = 1 in the model. We further carry out parameter tests and LB tests for the hypothetical model. A parameter test is used to verify the validity of the model; the LB test is used to check whether the residual sequence after model fitting is a white noise sequence.

Figure 6.

ACF and PACF diagrams of the earthquake time interval sequence after secondary difference.

The test results are shown in Table 2, and the above model contains an insignificant parameter (first-order seasonal autoregressive coefficient). Through comparison, the best seasonal order of the model, P = 0, and Q = 1 are determined.

SARIMA model	Adjusted R^2	Parameter test	LB test(p value)
ARIMA(9, 2, 1)x(0, 1, 0)22	0.817	pass	0.048
ARIMA(9, 2, 1)x(1, 1, 0) 22	0.846	pass	0.009
ARIMA(9, 2, 1)x(0, 1, 1) 22	0.869	pass	0.175
ARIMA(9, 2, 1)x(1, 1, 1) 22	0.869	not significant	0.115

Table 2.

Comparison of SARIMA models with different seasonal orders.

3.3 Model evaluation

The earthquake sequences used in the study have complex overlapping of different cycle periods. Therefore, when analyzing the periodicity of data, different seasonal periods are selected for short-, medium-, and long-term prediction. The center moving average of sequence X_d with span n is calculated, and the resulting sequence is marked as X_n, n ∈ (4,45). The difference between X_n and X_d is compared, the seasonal elimination degree of each cycle is observed, and the models of s ∈ (4,45) are fitted, respectively. To better compare the prediction effects between different models, this paper uses the RMSE to represent the fitting error of the model:

$$\text{RMSE}=\sqrt{\frac{1}{\mathrm{N}}\sum _{\mathrm{t}=1}^{\mathrm{N}}{\left({\mathrm{x}}_{\mathrm{t}}-{\mathrm{y}}_{\mathrm{t}}\right)}^2}$$E3

where x_t is the true value, and y_t is the fitted value. For the model passing the parameter test and LB test, we analyze its adjustment R² value and fitting RMSE to determine the fitting effect. The AIC is also commonly used to measure the complexity and goodness of fit of statistical models; the smaller the AIC value is, the better the model performance. The principle of model selection in this paper is to prioritize the model with a high R² value. When the difference between the R² values is not significant, the fitting RMSE of the model is compared to determine the optimal model, and the AIC value of the model is taken as an auxiliary reference. Table 3 provides the fitting of some models with different periods that have passed the parameter test and LB test.

Period	SARIMA model	Adjusted R^2	LB test (p value)	Fitting RMSE	AIC
Short period	ARIMA(9, 2, 1)x(0, 1, 1)5	0.870	0.194	12.049	2924.767
	ARIMA(9, 2, 1)x(0, 1, 1)6	0.911	0.215	11.785	2914.680
Medium period	ARIMA(9, 2, 1)x(0, 1, 1)22	0.869	0.175	13.083	2705.475
Long period	ARIMA(9, 2, 1)x(0, 1, 1)32	0.872	0.655	13.133	2579.307
	ARIMA(9, 2, 1)x(0, 1, 1)40	0.871	0.369	13.791	2473.724
	ARIMA(9, 2, 1)x(0, 1, 1)42	0.881	0.298	13.212	2435.224

Table 3.

Comparison of SARIMA models with different periods.

Therefore, we determine that the short-, medium-, and long-term prediction models of seismic data are ARIMA (9, 2, 1) $\times$ (0, 1, 1)₆，ARIMA(9, 2, 1) $\times$ (0, 1, 1)₂₂，and ARIMA(9, 2, 1) $\times$ (0, 1, 1)₄₂. The model-fitting effect is as follows (Figure 7):

Figure 7.

Fitting effect of models with different periods.

4. Magnitude grouping

By grouping the original sequence by magnitude, 390 points of data for the Ms. ≥ 3.0 sequence and 58 points of data for the Ms. ≥ 4.0 sequence are obtained. After analyzing and modeling each sequence according to the above steps, we consider fitting X_d3 with model ARIMA(7,2,1) $\times$ (0,1,1)_s and X_d4 with model ARIMA(3,2,1) $\times$ (0,1,0)_s. The comparison of each group’s models in different periods is shown in Tables 4 and 5.

Period	SARIMA model	Adjusted R^2	LB test (p value)	Fitting RMSE	AIC
Short period	ARIMA(7, 2, 1)x(0, 1, 1)7	0.872	0.058	13.901	2705.256
Medium period	ARIMA(7, 2, 1)x(0, 1, 1)14	0.881	0.061	14.059	2604.135
	ARIMA(7, 2, 1)x(0, 1, 1)16	0.891	0.201	13.970	2577.105
	ARIMA(7, 2, 1)x(0, 1, 1)23	0.849	0.179	17.243	2442.787
Long period	ARIMA(7, 2, 1)x(0, 1, 1)32	0.867	0.383	15.112	2349.765
	ARIMA(7, 2, 1)x(0, 1, 1)35	0.863	0.074	15.041	2294.26

Table 4.

Comparison of different periodic models of sequence X_d3.

Period	SARIMA model	Adjusted R^2	LB test (p value)	Fitting RMSE	AIC
Short period	ARIMA(3, 2, 0)x(0, 1, 0)7	0.809	0.108	73.242	472.971
Medium period	ARIMA(3, 2, 0)x(0, 1, 0)10	0.708	0.106	76.754	443.281
	ARIMA(3, 2, 0)x(0, 1, 0)11	0.775	0.672	60.594	414.368
	ARIMA(3, 2, 0)x(0, 1, 0)15	0.77	0.165	59.875	372.062
Long period	ARIMA(3, 2, 0)x(0, 1, 0)20	0.781	0.321	59.076	318.484

Table 5.

Comparison of different periodic models of sequence X_d4.

According to the above principles, we use s = 7, s = 16, and s = 32 to make short-term, medium-term, and long-term predictions for Xd3, respectively. The prediction model is ARIMA(7,2,1) $\times$ (0,1,1)₇，ARIMA(7,2,1) $\times$ (0,1,1)₁₆，and ARIMA(7,2,1) $\times$ (0,1,1)₃₂. For X_d4, due to the small amount of earthquake sequence data with an Ms. ≥ 4.0, it is difficult to capture the long period changing trend, and the fitting effect of those long period models is poor. In addition, the fitting effect of the short-period models is unsatisfactory. The X_d4 sequence only shows an obvious medium periodic. Therefore, the medium-term prediction of X_d4 is made with s = 11, and the prediction model is ARIMA (3,2,0) × (0,1,0)11. The fitting effects of different models of sequences are shown in Figures 8 and 9.

Figure 8.

Fitting effect of the optimal model with different periods of sequence X_d3. A. ARIMA(7,2,1) $\times$ (0,1,1)₇. B. ARIMA(7,2,1) $\times$ (0,1,1)₁₆. C. ARIMA(7,2,1) $\times$ (0,1,1)₃₂.

Figure 9.

Fitting effect of the optimal model of sequence X_d4.

The predicted results are compared with the true values, as shown in Figures 10 and 11. The long-term prediction of the X_d and X_d3 series shows that the model can capture their periodic changes, but the overall prediction value is higher than the real value.

Figure 10.

Comparison of prediction effects of optimal models with different periods of sequences X_d and X_d3. A. Prediction of X_d. B. Prediction of X_d3.

Figure 11.

Prediction effect of the optimal model of sequence X_d4.

Short-term prediction of X_d4 shows that ARIMA(3,2,0) $\times$ (0,1,0)₁₁ can predict well the development trend of this series in recent few times and does not show a high prediction result. The predicted RMSE of each model is shown in Table 6.

Series	Optimal model	Predicted RMSE
${\mathrm{X}}_{\mathrm{d}}$	ARIMA(9, 2, 1)x(0, 1, 1)6	35.546
	ARIMA(9, 2, 1)x(0, 1, 1)22	35.437
	ARIMA(9, 2, 1)x(0, 1, 1)42	35.206
${\mathrm{X}}_{\mathrm{d}3}$	ARIMA(7, 2, 1)x(0, 1, 1)42	37.385
	ARIMA(7, 2, 1)x(0, 1, 1)16	37.571
	ARIMA(7, 2, 1)x(0, 1, 1)32	37.060
${\mathrm{X}}_{\mathrm{d}4}$	ARIMA(3, 2, 0)x(0, 1, 0)11	56.091

Table 6.

Prediction RMSE of different periodic models for each magnitude series.

The average level of the true and predicted values of the series X_d and X_d3 is calculated. The average error of the models is 17.578 days at the lowest and 32.967 days at the highest. With reference to the general accuracy of the prediction of the earthquake occurrence time, the prediction error can be considered to be within an acceptable range. We use the predicted value to subtract the average error to correct the series predicted value. Taking ARIMA(9, 2, 1) x (0, 1, 1)₂₂ for X_d and ARIMA(7, 2, 1) x (0, 1, 1)₃₂ for X_d3 as examples, the correction effect is shown in Figure 12, which shows that the models have good prediction performance for trend and periodicity.

Figure 12.

Prediction effects after correction of sequences X_d and X_d3.

Since the prediction step of sequences X_d and X_d3 is more than 40, which may reduce the prediction accuracy, we consider making 15-step predictions for sequence X_d and 5-step predictions for sequence X_d3, which has a relatively small data volume. Taking the above two models as examples, the prediction results are shown in Figure 13. The predicted RMSEs are 10.6860 and 8.8009, respectively. There is no obvious trend of the earthquake occurrence interval predicted by sequence X_d, while the prediction sequence of X_d3 has a slightly increasing trend.

Figure 13.

Prediction effects of sequences X_d and X_d3 after reducing prediction times.

Two new datasets for an earthquake sequence with an Ms. ≥ 4.0 in the Longmenshan fault zone as of March 6, 2022 are added, and the ARIMA(3, 2, 0) x (0, 1, 0)₁₁ model is used to predict the time interval series. Comparing the predicted value with the true value, as shown in Figure 12, the predicted RSME is 55.7112. It is found that the trend is still captured well, and the model accuracy can be further improved on this basis (Figure 14).

Figure 14.

Prediction results of new data (Ms ≥ 4.0) in the Longmenshan fault zone in 2022.

5. Conclusion

This paper proposes to use SARIMA to model the time series from the perspective of the time interval of the Longmenshan fault zone, analyze the hidden information of the earthquake time series, and predict the next earthquake occurrence time. According to the model analysis and prediction results, the following conclusions are drawn: (1) The SARIMA model is applicable to the analysis and prediction of earthquake time interval series. The optimal model adjusted R² value of each series is above 0.86, up to 0.911. (2) For long-term prediction, the models of series X_d and X_d3 have higher prediction values than the true values, and the prediction performance for lower values (time interval approaching 0) is relatively poor. (3) In short-term prediction, the optimal models of sequences X_d and X_d3 have good prediction effects and can predict the sequence periodicity well. The prediction result of X_d3 shows a slight fluctuation growth trend; that is, the number of earthquakes with an Ms. ≥ 3.0 in the Longmenshan fault zone decreases slightly. The periodicity of sequence X_d4 is obvious, and in short-term prediction, the model can capture well its development. The prediction trend is consistent with the real situation. The prediction accuracy of the model can be further improved.