Applicability of a CEEMD–ARIMA Combined Model for Drought Forecasting: A Case Study in the Ningxia Hui Autonomous Region

Abstract

In the context of global warming, the increasing frequency of drought events has caused negative impacts on agricultural productivity and societal activities. However, the drought occurrences have not been well predicted by any single model, and precipitation may show nonstationary behavior. In this study, 60 years of monthly precipitation data from 1960 to 2019 for the Ningxia Hui Autonomous Region were analyzed. The standard precipitation index (SPI) was used to classify drought events. This study combined the strengths of autoregressive integrated moving average (ARIMA) and complementary ensemble empirical mode decomposition (CEEMD) to predict drought. First, based on the precipitation dataset, the SPI at timescales of 1, 3, 6, 9, 12, and 24 months was calculated. Then, each of these SPI time series was predicted using the ARIMA model and the CEEMD–ARIMA combined model. Finally, the models′ performance was compared using statistical metrics, namely, root-mean-square error (RMSE), mean absolute error (MAE), Kling–Gupta efficiency (KGE), Willmott index (WI), and Nash–Sutcliffe efficiency (NSE). The results show that the following: (1) Compared with the ARIMA forecast value, the prediction results of the CEEMD–ARIMA model were in good agreement with the SPI values, indicating that the combined model outperformed the single model. (2) Two different models obtained the lowest accuracy for the SPI1 prediction and the highest accuracy for the SPI24 prediction. (3) The CEEMD–ARIMA model achieved higher prediction accuracy than the ARIMA model at each time scale. The most precise model during the test phase was the CEEMD–ARIMA model at SPI24 at Xiji Station, with error measures of MAE = 0.076, RMSE = 0.100, NSE = 0.994, KGE = 0.993, and WI = 0.999. Such findings will be essential for government to make decisions.

1. Introduction

Drought, which frequently occurs around the world, causes tremendous losses to agricultural production and economic operation [1]. For instance, the 2014 California drought was a record-breaking event that cost the United States USD 2.2 billion [2]. Drought is one of the costliest disasters that humankind faces all over the world [3]. With climate changes and temperature increases, droughts are becoming more and more frequent. Quantitative studies on drought will help countries to avoid damage caused by climate disasters in the future. Improved drought-monitoring ability has obvious significance in city development, which could help in the creation of drought management agencies.

Droughts are generally categorized into four types: meteorological droughts, agricultural droughts, hydrological droughts, and socioeconomic droughts [4]. Meteorological drought initiates when precipitation presents with volumes below normal in a particular place. Such anomalies go on to affect agriculture and hydrology (agricultural and hydrological droughts, respectively) [5]. This study focuses on meteorological drought. An effective method of meteorological drought detection is to use drought indices [6]. Drought indices can be used to quantitatively evaluate the influence range of a drought. In recent years, a variety of indices have been used to study drought, such as the Palmer drought severity index (PDSI) [7], standardized precipitation index (SPI) [8], reconnaissance drought index (RDI) [9], and temperature vegetation drought index (TVDI) [10]. Among them, the SPI is widely used in drought research at home and abroad due to its variable time scale, and only precipitation data are used for calculation [6,11,12]. Based on these two major advantages of the SPI—and particularly its ability to describe drought on multiple time scales—Oliveira-Júnior et al. [13] evaluated the drought severity in northern and northwestern Rio de Janeiro State (SRJ) regions from 1967 to 2013. Wu et al. [14] and Xu et al. [15] used the SPI to investigate the characteristics of meteorological drought in China. Łabędzki [16] adopted the SPI to estimate meteorological drought frequency in the central part of Poland from 1861 to 2005. Therefore, the SPI was used in this study due to its wide acceptability and advantages in drought research. Drought prediction provides an early warning to decision makers in disaster management. However, it is often notably challenging to obtain proper forecasts, because of the complexity of measuring the precision of a time series [17,18]. In recent years, numerous models have been used in drought forecasting, such as the autoregressive integrated moving average (ARIMA), artificial neural network (ANN), and support-vector regression (SVR) [18,19,20,21,22]. Among them, ARIMA is widely used in drought prediction because of its flexibility and richer information on time-related changes [23]. Nevertheless, the forecasting accuracy of a single model cannot meet the needs of drought prediction. Moreover, the precipitation data have nonlinear and nonstationary characteristics. Therefore, hybrid models are used in drought research to improve the prediction accuracy. Empirical mode decomposition (EMD) has apparent advantages in the processing of nonlinear and nonstationary signal time–frequency sequences. Özger et al. [24] used EMD for decomposing self-calibrated Palmer drought severity index (sc-PDSI) time series into their sub-bands on drought prediction, but this decomposition method has the problem of mode aliasing. As a further improvement of EMD, ensemble empirical mode decomposition (EEMD) effectively reduces the occurrence of mode aliasing. Libanda et al. [25] used it to understand consecutive dry days, but in this decomposition method, Gaussian white noise was added to the original signal, and its influence on the results could not be ignored. Therefore, since CEEMD has achieved great results in many fields based on advantages in processing signals [26,27,28], and it can effectively reduce the residual white noise and process nonlinear and nonstationary signals, a new drought prediction method was proposed, combining the ARIMA model and complementary ensemble empirical mode decomposition (CEEMD).

As mentioned above, with the nonlinear and nonstationary characteristics of precipitation data, it is important to accurately predict the occurrence of drought. The main objectives of the present study are as follows: to (1) quantify the precipitation situation by multi-timescale SPI, (2) develop the ARIMA model and then propose the hybrid model by combining the strengths of ARIMA with CEEMD, and (3) evaluate the efficiency of the ARIMA model and the CEEMD–ARIMA model according to the evaluated indices.

2. Study Area

The Ningxia Hui Autonomous Region extends from 104°17′ E to 107°39′ E and from 35°14′ N to 39°23′ N, with altitudes mostly above 1000 m (Figure 1). Helan Mountain is located in the northwest of the province. As a natural barrier, this mountain reduces intrusion from Tengger Desert quicksand and cold northwest winds into the Ningxia Hui Autonomous Region. The southern part of the province is the Liupan mountainous area. As the wettest region in Ningxia Hui Autonomous Region, it has a humid climate and dense jungles. The area from the Helan Mountains to the Weining Plain has an arid climate, and the area from the Weining Plain to the Liupan mountainous area has a semiarid climate. The Liupan mountainous area consists of jungle, wet, and pleasantly cool areas. The climate characteristics of Ningxia vary widely from north to south on spatial and temporal scales. The annual average temperature decreases from the north to the south in Ningxia, while the annual precipitation shows the opposite tendency. The mean annual temperature is between 5.3 and 9.9 °C, with the southern part below 7 °C, the central part above 7 °C, and the northern part above 8 °C. The mean annual precipitation is between 150 mm and 600 mm, and the average annual water surface evaporation in Ningxia is 1250 mm.

3. Materials and Methods

3.1. Data Sources

Monthly precipitation datasets from 10 meteorological stations were used in this study. The datasets from January 1960 to December 2019 were obtained from the Ningxia Hui Autonomous Region weather station in the National Meteorological Data Center (http://data.cma.cn/ accessed on 13 March 2020). The elevation data were obtained from the Geospatial Data Cloud (http://www.gscloud.cn/search accessed on 28 June 2021). Due to the terrain inclines from the southwest to the northeast in the Ningxia Hui Autonomous Region, there is a large spatial difference in climate characteristics from the south to the north. Since this paper focuses on the applicability of the combined model in drought prediction, three meteorological stations were selected as representative stations of the southern, northern, and central Ningxia Hui Autonomous Region.

Table 1. Information about meteorological stations in the sample.

Station Number	Station Name	Longitude/°E	Latitude/°N	Altitude/m
53519	Huinong	106.46	39.13	1092.5
53810	Tongxin	105.54	36.58	1339.3
53903	Xiji	105.43	35.58	1916.5

shows the information of three representative stations. The observed and predicted values of 10 sites are visualized by the empirical Bayesian kriging interpolation method of ArcGIS.

3.2. Research Methods

3.2.1. SPI

The SPI was developed by McKee et al. [29] to quantify precipitation on different time scales. It can be calculated based solely on precipitation. The time scale of the SPI is variable. Short-timescale SPI could reflect the water supply of crops. On a 1-month time scale, the SPI can be used to reflect short-term precipitation conditions. The SPI on a 3-month time scale can be used to analyze seasonal changes in precipitation. The calculated SPI data for February, May, August, and November reflect winter, spring, summer, and autumn drought conditions, respectively. The 6-month time scale SPI can be used to reflect mid-term precipitation conditions. On a 9-month time scale, the SPI can be used to characterize groundwater level changes over a longer period of time. The 12-month time scale SPI and 24-month time scale SPI are good indicators of long-term drought conditions [30]. The SPIs of the 1-, 3-, 6-, 9-, 12-, and 24-month time scales are denoted as SPI1, SPI3, SPI6, SPI9, SPI12, and SPI24, respectively. The computation procedure of SPI followed the method of Lloyd-Hughes and Saunders [31]:

$$\mathrm{SPI}=B\left(t-\frac{c_0+c_{1}t+c_2{t}^2}{1+d_{1}t+d_2{t}^2+d_3{t}^3}\right)$$

(1)

where B is the positive and negative coefficient of probability density; for $B=-1$, $t=\sqrt{ln\frac{1}{G{\left(x\right)}^2}},$ and for $B=1$, $t=\sqrt{ln\frac{1}{{\left(1-G\left(x\right)\right)}^2}}$, where $G\left(x\right)$ is a cumulative probability [32]. The constants $c_0$ = 2.515517, $c_1$ = 0.802853, $c_2$ = 0.010328, $d_1$ = 1.432788, $d_2$ = 0.189269, and $d_3$ = 0.001308.

The classification of drought based on SPI values is shown in

Table 2. Drought classification based on SPI.

SPI Value	Category
SPI > −0.5	No drought
−1.0 < SPI ≤ −0.5	Mild drought
−1.5 < SPI ≤ −1.0	Moderate drought
−2.0 < SPI ≤ −1.5	Severe drought
SPI ≤ −2.0	Extreme drought

[31].

3.2.2. ARIMA Model

The ARIMA model developed by Box and Jenkins [33] includes three basic types: autoregressive (AR) models, moving average (MA) models, and the combined AR and MA (ARMA) models. AR, MA, and ARMA can be used when the data are stationary [19]. Nonstationary and nonwhite noise sequences can be predicted by the ARIMA model. First, a stationary time series is obtained from a nonstationary series through d-order difference. Then, the ARMA model is used for prediction. The formula of the ARIMA (p, d, q) model is written as follows:

$$X_t={\omega }_1{X}_{t-1}+{\omega }_2{X}_{t-2}+\cdots +{\omega }_p{X}_{t-p}+h_t-{\theta }_1{h}_{t-1}-{\theta }_2{h}_{t-2}-\cdots -{\theta }_q{h}_{t-q}$$

(2)

where $X_t$ is a time-series value, and ${\omega }_i\left(i=1,2,\cdots ,p\right)$ and ${\theta }_j\left(j=1,2,\cdots ,q\right)$ are the autoregressive coefficient and moving average coefficient, respectively. $h_t$ is a white noise sequence, $h_t~N\left(0,{\sigma }^2\right)$.

The modeling process of the ARIMA model is as follows:

First, each time series should go through stationary testing. In this paper, the augmented Dickey–Fuller test (ADF) is used for judgment. If it is a nonstationary time series, the d-order difference of the series is determined. Then, the value range of the model order should be determined. The value range of P and Q is determined according to the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the data. After that, the Akaike information criterion (AIC) and Bayesian information criterion (BIC) are used to determine the order of the model. The formulae of the AIC and BIC are as follows:

$$AIC\left(p,q\right)=Nln{\sigma }^2\left(p,q\right)+2\left(p+q+1\right)$$

(3)

$$BIC\left(p,q\right)=Nln{\sigma }^2\left(p,q\right)+\left(p+q+1\right)lnN$$

(4)

where N is the number of parameters. In different combinations of p and q, the combination corresponding to the minimum value of AIC and BIC is selected to obtain the optimal ARIMA model. The datasets were divided by grid search and cross-validation; 80% of the data were selected as the training set for model prediction, and 20% of the data were selected as the test set.

3.2.3. CEEMD

As proposed by Yeh et al. [34], CEEMD has apparent advantages in the processing of nonlinear and nonstationary signal time–frequency sequences. It can adaptively decompose the original sequence into several intrinsic mode function (IMF) components, with different scales that are mutually independent and a residual trend quantity. The steps are as follows:

A group of white noise includes positive noise and negative noise. P(t) is the original sequence, with n groups of auxiliary white noise added to the positive noise sequence Y₁ and negative noise sequence Y₂. Now, the total number of sequences obtained is 2n,

$$\left[\begin{array}{c}{Y}_1\\ Y_2\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{c}P\\ N\end{array}\right]$$

(5)

where N is an auxiliary sequence. The obtained sequences are decomposed by EMD to obtain m IMF components, and each group of components is denoted as $C_{ij}^{+}\left(t\right)$ and $C_{ij}^{-}\left(t\right)$ (i = 1, …, n, j = 1, …, m). $C_{ij}^{+}\left(t\right)$ and $C_{ij}^{-}\left(t\right)$ of the IMF components in each group are averaged to obtain the ${\mathrm{IMF}}_j$ value.

$${\mathrm{IMF}}_j=\frac{1}{2n}{\displaystyle \sum }_{i=1}^n\left(C_{ij}^{+}\left(t\right)+C_{ij}^{-}\left(t\right)\right)$$

(6)

Take the decomposed IMF components as the final result. The original sequence is decomposed into:

$$P_{\left(t\right)}={\displaystyle \sum }_{j=1}^m{\mathrm{IMF}}_{j\left(t\right)}+R_{\left(t\right)}$$

(7)

where $R_{\left(t\right)}$ is a residual trend quantity.

3.2.4. CEEMD–ARIMA Combined Model

The original sequences with strong volatility are decomposed by CEEMD to obtain a set of IMF components with low volatility, improving the predictability of the sequences predicted by the ARIMA model. The CEEMD and ARIMA models are combined using Python to form the CEEMD–ARIMA model. The steps are as follows:

Step 1: The SPI sequence is imported into the CEEMD for decomposition to obtain IMF1, IMF2…, IMFn, and residual trends from high frequency to low frequency.

Step 2: The sequences decomposed by CEEMD are imported into the ARIMA model. The stationarity of each component is tested through the ARIMA model. After the order and prediction are determined, the prediction result is obtained. The predicted results are denoted as P1, P2, and Pn + 1.

Step 3: Finally, sum P1, P2, and…, Pn + 1 as the combined model’s predicted results.

$$P={\displaystyle \sum }_{i=1}^{n+1}{P}_i$$

(8)

The modeling process of the CEEMD–ARIMA model is as shown in Figure 2.

3.2.5. Evaluation Index

To evaluate the performance of the ARIMA model and the CEEMD–ARIMA model, statistical criteria such as root-mean-square error (RMSE), mean absolute error (MAE), Kling–Gupta efficiency (KGE), Willmott index (WI), and Nash–Sutcliffe efficiency (NSE) were used. A model with the lowest RMSE and MAE and the highest KGE, WI, and NSE was selected and proposed as an appropriate model. The formulae are as follows [35,36,37,38]:

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left(y_i-\widehat{y_i}\right)}^2}$$

(9)

$$MAE=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\left|y_i-\widehat{y_i}\right|$$

(10)

$$KGE=1-\sqrt{{\left(CC-1\right)}^2+{\left(\alpha -1\right)}^2+{\left(\beta -1\right)}^2}$$

(11)

$$WI=\left|1-\left[\frac{{\sum }_{i=1}^N{\left(y_i-\widehat{y_i}\right)}^2}{{\sum }_{i=1}^N{\left(\left|y_i-\overline{y}\right|+\left|\widehat{y_i}-\overline{y}\right|\right)}^2}\right]\right|$$

(12)

$$NSE=1-\frac{{\sum }_{i=1}^N{\left(y_i-\widehat{y_i}\right)}^2}{{\sum }_{i=1}^N{\left(y_i-\overline{y}\right)}^2}$$

(13)

where $y_i$ is the observed value, $\overline{y}$ is the average value of $y_i$, $\widehat{y_i}$ is the forecasted value, and $N$ is the total data size of $y_i$. $CC$, $\alpha$, and $\beta$ in the KGE index illustrate the correlation coefficient, the standard deviation ratio, and the average ratio of the observed and forecasted data, respectively.

4. Results

The calculation of the SPI and the fitting of the ARIMA model were both accomplished on the Python 3.7 platform.

4.1. SPI Values at Different Time Scales

The research applicability of the CEEMD–ARIMA model in drought is mainly through the prediction of the SPI on the time scales of 1, 3, 6, 9, 12, and 24 months. Monthly precipitation data from 10 meteorological stations in the Ningxia Hui Autonomous Region during 1960–2019 were used to calculate the SPI. The calculated SPI results characterize drought conditions, as shown in Table 2. Huinong, Tongxin, and Xiji were taken as examples to demonstrate multiple-timescale SPIs, and the calculated SPIs of the sample stations are shown in Figure 3. Through the Mann–Kendall trend test, the SPI12 and SPI24 sequences of Tongxin Station and the SPI9, SPI12, and SPI24 sequences of Xiji Station have a decreasing trend. The other sequences have no trend.

4.2. The ARIMA Modeling and Prediction

The first 80% of the calculated SPI data were used as observation training data, and the last 20% of the data were used as prediction comparison data. That is, the data from 1960 to 2007 were used as the training set, and the data from 2008 to 2019 were used as the test set. The stability of the 80% training data should be judged before prediction. If the data are a stable series, then $d=0$ in the ARIMA model can be used for prediction; if not, then $d\ne 0$. Through the ADF test, the p-values of all SPI sequences of the sample sites were less than 0.05 (

Table 3. Unit root test of SPI original sequence.

Example Stations	SPI Series	ADF	Critical Value			p-Value
			1%	5%	10%
Huinong	SPI1	−20.0550	−3.4418	−2.8666	−2.5694	0.0000
	SPI3	−9.6732	−3.4419	−2.8666	−2.5695	1.2610 × 10−16
	SPI6	−6.9028	−3.4420	−2.8667	−2.5695	1.2693 × 10−9
	SPI9	−5.3241	−3.4423	−2.8668	−2.5696	4.8806 × 10−6
	SPI12	−4.7455	−3.4423	−2.8668	−2.5696	6.9075 × 10−5
	SPI24	−4.1882	−3.4423	−2.8668	−2.5696	0.0007
Tongxin	SPI1	−21.6155	−3.4418	−2.8666	−2.5694	0.0000
	SPI3	−9.6077	−3.4419	−2.8666	−2.5695	1.8469 × 10−16
	SPI6	−6.7922	−3.4420	−2.8667	−2.5695	2.3486 × 10−9
	SPI9	−4.9104	−3.4423	−2.8668	−2.5696	3.3288 × 10−5
	SPI12	−4.4071	−3.4423	−2.8668	−2.5696	0.0003
	SPI24	−3.7087	−3.4423	−2.8668	−2.5696	0.0040
Xiji	SPI1	−22.0945	−3.4418	−2.8666	−2.5694	0.0000
	SPI3	−10.7739	−3.4419	−2.8666	−2.5695	2.3469 × 10−19
	SPI6	−7.3216	−3.4420	−2.8667	−2.5695	1.1900 × 10−10
	SPI9	−4.1113	−3.4423	−2.8668	−2.5696	0.0009
	SPI12	−3.4578	−3.4422	−2.8668	−2.5696	0.0091
	SPI24	−3.3257	−3.4423	−2.8668	−2.5696	0.0138

). Therefore, all SPI sequences were stationary time series.

Because the SPI series of the three sites were stationary time series, the ARMA model was selected for the prediction. ACF and PACF were used to rank the ARMA model, and p- and q-values corresponding to the minimum AIC and BIC values were selected. The model ranking results of each sequence are shown in

Table 4. Model order based on SPI values of six time scales.

Example Stations	SPI Series	Model Select	AIC	BIC	Model Order Estimation
Huinong	SPI1	ARMA	1826.071	1839.804	ARMA (1, 0)
	SPI3	ARMA	1631.778	1650.079	ARMA (0, 2)
	SPI6	ARMA	1398.692	1412.404	ARMA (1, 0)
	SPI9	ARMA	1026.739	1045.006	ARMA (1, 0)
	SPI12	ARMA	538.884	579.946	ARMA (5, 2)
	SPI24	ARMA	64.999	87.725	ARMA (3, 0)
Tongxin	SPI1	ARMA	1937.225	1950.959	ARMA (1, 0)
	SPI3	ARMA	1593.929	1612.230	ARMA (0, 2)
	SPI6	ARMA	1302.638	1343.776	ARMA (5, 2)
	SPI9	ARMA	957.282	970.982	ARMA (1, 0)
	SPI12	ARMA	536.069	586.256	ARMA (7, 2)
	SPI24	ARMA	43.954	62.136	ARMA (2, 0)
Xiji	SPI1	ARMA	2012.614	2026.347	ARMA (0, 1)
	SPI3	ARMA	1628.778	1647.078	ARMA (0, 2)
	SPI6	ARMA	1453.959	1472.242	ARMA (2, 0)
	SPI9	ARMA	1061.371	1075.071	ARMA (1, 0)
	SPI12	ARMA	575.482	616.544	ARMA (5, 2)
	SPI24	ARMA	31.131	62.949	ARMA (3, 2)

. The optimal model of the SPI series at various time scales was applied to predict the SPI series of the three stations from 2008 to 2019.

4.3. The CEEMD–ARIMA Combined Model

Multiscale SPI was decomposed by CEEMD. After several parameters were modified and compared, when the Gaussian white noise logarithm was 100, the total number of modes (not including the trend) was 7, and the standard deviation of the original time series was multiplied by 0.2. CEEMD had the best decomposition effect. Seven IMF components and one trend item were obtained by the CEEMD decomposition of the SPI3 sequence of Xiji Station (Figure 4). The trend term represents the general trend of a sequence over time. As shown in Figure 4, the fluctuation range of the IMF component obtained by decomposition is smaller than that of the original sequence, and with the gradual progress of decomposition, the fluctuation of the component tends to be smooth. Therefore, the predictability of the subsequence obtained after decomposition is higher than that of the original sequence.

A stationary test was carried out for the components decomposed by CEEMD. IMF1~IMF6 were stationary time series; therefore, the ARMA model was selected for prediction. IMF7 and RES were nonstationary time series, and the ARIMA model was used for prediction after the stationary series was obtained by data difference. AIC and BIC were used to determine the order of the model, and the ranking results are shown in

Table 5. Model order based on IMFs of the SPI3 sequence at Xiji Station.

SPI Series	Decompose Results	Model Select	Model Order Estimation
SPI3	IMF1	ARMA	ARMA (1, 1)
	IMF2	ARMA	ARMA (2, 5)
	IMF3	ARMA	ARMA (4, 2)
	IMF4	ARMA	ARMA (4, 5)
	IMF5	ARMA	ARMA (4, 6)
	IMF6	ARMA	ARMA (2, 1)
	IMF7	ARIMA	ARIMA (4, 1, 1)
	Res	ARIMA	ARIMA (3, 1, 1)

.

In this paper, the data from 1960 to 2007 were used as observation training data. Therefore, only this part of the SPI sequence was decomposed by CEEMD and then predicted by the ARIMA model, and the sum of the predicted results of each component was used as the final prediction result of the SPI sequence. The prediction comparison plot includes the actual calculated SPI values and the predicted SPI values of the CEEMD–ARIMA model and the ARIMA model, as shown in Figure 5, Figure 6 and Figure 7 for Huinong, Tongxin, and Xiji sample stations, respectively.

Figure 5a, Figure 6a and Figure 7a show that there is a big difference between the predicted value of the ARIMA model and the actual value of SPI at the 1-month time scale, and the predicted result of the CEEMD–ARIMA combination model at the SPI1 time scale is better than that of the ARIMA model. The poor stationarity of the 1-month time series resulted in bad prediction results of the ARIMA model. As the time scale increases and the data stationarity improves, the ARIMA prediction results become closer and closer to the actual situation. The prediction of the ARIMA model in the combined model relies on the stable basis provided by the CEEMD. At the time scales of 1 month and 3 months, the predicted value of the CEEMD–ARIMA model is different from the actual value. At the time scales of 6, 9, 12, and 24 months, the predicted values of the combined model are less different. According to the comparison diagram of SPI3 and SPI6 in Figure 5, Figure 6 and Figure 7, the prediction results of the combined model in the extreme drought year (SPI ≤ −2) are closer to the actual situation, indicating that the combined model is more suitable for the study of extreme drought than the single model. With increasing time scale, the difference between the predicted value of the model and the actual calculated value tends to decrease. The two models′ predicted values of SPI12 and SPI24 were very close to the actual values, and the combined model was closer than the single model.

In 2009, due to high temperature and drier weather, a severe large-scale drought occurred in Ningxia, which spread from the central region. From the predicted values of the combined model shown in Figure 5, Figure 6 and Figure 7, drought first appeared near Tongxin Station, and then around Huinong Station and Xiji Station, which is consistent with the record in the China Meteorological Network. The combined model predicts that there would be no drought or mild drought at Huinong Station and Tongxin Station in 2016, but drought at Xiji Station, which is also consistent with the record. In July 2017, large-scale precipitation occurred in Ningxia, with heavy rain in some areas. In 2019, the precipitation in the whole region was relatively high, with annual precipitation of 341.7 mm, and the southern mountainous region was rainy for eight consecutive months. The prediction results of the combined model for drought in these periods were nearly consistent with the actual situation. This indicates that the prediction of the combined model has a high consistency with the actual situation, illustrating that the combined model is suitable for studying drought prediction.

The prediction results of the two models were evaluated using MAE, RMSE, NSE, KGE, and WI. With increasing time scale, the MAE and RMSE values of the two models decreased, while the NSE, KGE, and WI generally increased (

Table 6. The statistical criteria of the ARIMA and CEEMD–ARIMA models.

Example Stations	SPI Series	Model	Training					Testing
			MAE	RMSE	NSE	KGE	WI	MAE	RMSE	NSE	KGE	WI
Huinong	SPI1	ARIMA	0.634	0.850	−31.759	−3.881	0.204	0.667	0.892	−48.453	−4.992	0.150
		CEEMD–ARIMA	0.459	0.580	−0.020	0.440	0.830	0.465	0.596	−0.058	0.420	0.817
	SPI3	ARIMA	0.535	0.708	−0.544	0.284	0.750	0.549	0.723	−0.663	−0.250	0.730
		CEEMD–ARIMA	0.393	0.502	0.497	0.654	0.894	0.407	0.526	0.448	0.632	0.886
	SPI6	ARIMA	0.429	0.609	0.363	0.643	0.867	0.440	0.618	−0.013	0.452	0.783
		CEEMD–ARIMA	0.244	0.312	0.860	0.886	0.962	0.250	0.321	0.808	0.861	0.954
	SPI9	ARIMA	0.304	0.434	0.711	0.816	0.934	0.315	0.460	0.384	0.671	0.850
		CEEMD–ARIMA	0.143	0.188	0.927	0.893	0.981	0.150	0.199	0.906	0.876	0.977
	SPI12	ARIMA	0.219	0.348	0.883	0.921	0.972	0.226	0.363	0.604	0.783	0.896
		CEEMD–ARIMA	0.125	0.186	0.925	0.927	0.982	0.129	0.194	0.884	0.923	0.972
	SPI24	ARIMA	0.149	0.233	0.939	0.953	0.985	0.157	0.248	0.670	0.831	0.911
		CEEMD–ARIMA	0.067	0.087	0.957	0.978	0.990	0.069	0.090	0.954	0.972	0.989
Tongxin	SPI1	ARIMA	0.711	0.909	−87.660	−7.274	0.127	0.724	0.918	−100.523	−8.116	0.115
		CEEMD–ARIMA	0.452	0.557	0.415	0.133	0.879	0.466	0.574	0.374	0.130	0.868
	SPI3	ARIMA	0.578	0.729	−0.286	0.360	0.783	0.606	0.740	−0.395	0.133	0.758
		CEEMD–ARIMA	0.343	0.416	0.787	0.377	0.952	0.349	0.424	0.750	0.369	0.944
	SPI6	ARIMA	0.437	0.588	0.489	0.704	0.890	0.467	0.626	0.357	0.499	0.859
		CEEMD–ARIMA	0.207	0.275	0.934	0.553	0.985	0.224	0.296	0.894	0.541	0.974
	SPI9	ARIMA	0.323	0.472	0.731	0.791	0.938	0.325	0.482	0.632	0.783	0.915
		CEEMD–ARIMA	0.138	0.181	0.960	0.804	0.991	0.142	0.187	0.952	0.797	0.988
	SPI12	ARIMA	0.235	0.336	0.873	0.916	0.969	0.239	0.341	0.823	0.853	0.957
		CEEMD–ARIMA	0.090	0.122	0.984	0.967	0.996	0.096	0.130	0.976	0.962	0.994
	SPI24	ARIMA	0.159	0.247	0.937	0.956	0.985	0.172	0.253	0.921	0.944	0.980
		CEEMD–ARIMA	0.062	0.079	0.996	0.975	0.999	0.065	0.083	0.992	0.972	0.998
Xiji	SPI1	ARIMA	0.782	0.961	−116.898	−10.640	0.237	0.825	1.036	−126.675	−36.326	0.224
		CEEMD–ARIMA	0.570	0.706	0.269	0.205	0.846	0.584	0.739	0.256	0.182	0.831
	SPI3	ARIMA	0.574	0.731	−0.313	0.370	0.774	0.649	0.820	−0.487	−0.528	0.752
		CEEMD–ARIMA	0.391	0.481	0.717	0.794	0.939	0.407	0.508	0.689	0.776	0.930
	SPI6	ARIMA	0.492	0.657	0.332	0.529	0.877	0.547	0.670	0.313	0.262	0.869
		CEEMD–ARIMA	0.235	0.297	0.930	0.842	0.984	0.247	0.309	0.923	0.835	0.981
	SPI9	ARIMA	0.346	0.490	0.711	0.768	0.936	0.412	0.576	0.696	0.482	0.933
		CEEMD–ARIMA	0.211	0.279	0.948	0.934	0.988	0.221	0.291	0.940	0.923	0.985
	SPI12	ARIMA	0.229	0.354	0.890	0.888	0.980	0.245	0.377	0.890	0.625	0.974
		CEEMD–ARIMA	0.102	0.136	0.987	0.937	0.997	0.107	0.141	0.987	0.921	0.997
	SPI24	ARIMA	0.158	0.233	0.950	0.753	0.989	0.188	0.285	0.949	0.514	0.988
		CEEMD–ARIMA	0.069	0.087	0.995	0.994	0.999	0.076	0.100	0.994	0.993	0.999

), indicating that the prediction accuracy of the two models gradually improved with increasing time scale, and reached the maximum at the 24-month time scale. For example, the ARIMA model implemented at SPI1 had a WI = 0.150, and at SPI24 had a WI = 0.911, at Huinong Station. The evaluation index values of the two models were compared. At all time scales, the MAE and RMSE values of the CEEMD–ARIMA model were lower than those of the ARIMA model, and the NSE, KGE, and WI values were higher than those of the ARIMA model, illustrating the higher prediction accuracy of the combined model, which is more suitable for the prediction of multiscale SPI. At the 1-month time scale, the prediction accuracy of the combined model was much higher than that of the single model, and the prediction accuracy of SPI9, SPI12, and SPI24 was slightly higher than that of the single model. The most precise model during the test phase was the CEEMD–ARIMA model at SPI24 at Xiji Station with MAE = 0.076 RMSE = 0.100, NSE = 0.994, KGE = 0.993, and WI = 0.999. With increasing time scale, the fluctuation of the SPI series tended to be flat, and the fitting degree of the ARIMA model to the calculated SPI value was gradually improved.

The actual calculated SPI values, the ARIMA-predicted values, and the CEEMD–ARIMA-predicted values of 10 sites in 2019 were visualized by the empirical Bayesian kriging interpolation method in ArcGIS. The SPI at different time scales is suitable for different analyses. In this paper, SPI3 was selected to show the drought situation of spring, summer, autumn, and winter in the Ningxia Hui Autonomous Region, which can be used to analyze the seasonal variation in drought. As shown in Figure 8, the prediction of the CEEMD–ARIMA model was closer to the actual situation than that of the ARIMA model, and the predicted results were consistent with the actual approximation. In the summer of 2019, the precipitation in southern Ningxia Hui Autonomous Region was abnormal and excessive. The instability of precipitation data led to a big difference between the ARIMA-predicted values and the observed SPI values. Based on the advantage of CEEMD in nonstationary signal processing, the precision of the combined model is good, and is consistent with actual states.

5. Discussion

Recent studies have highlighted the superiority of EMD, EEMD, and CEEMD in forecasting aspects [39,40,41,42]. Ali et al. [43] forecasted rainfall at a monthly time scale. They resolved the non-stationarity challenges faced by rainfall forecasting models via CEEMD. Their study, which used a hybrid model for forecasting rainfall, achieved a WI value of 0.966. Their findings indicated that the CEEMD effectively avoids the non-stationarity in rainfall forecasting. Previous studies in drought prediction have used single models. For example, using the ARIMA model, Shatanawi et al. [20] predicted 3 out of 4 actual moderate droughts at Amman and Mafraq Stations. Similarly, Liu et al. [6] used the ARMA model to predict SPI9; the results showed that the prediction results of ARMA model were closer to the observed values in Longkou Station (the average relative error is 20.39%). However, when considering the entire stations, the highest standard error was as much as 43.69%. This result showed that more emphasis should be given to studying larger areas, which will be essential in regional drought management to make decisions. In this study, data stability is considered, along with the applicability of the CEEMD–ARIMA model over a large area.

The non-stationarity of data affects the drought prediction results of the model. The prediction result of the ARIMA model for SPI1 was significantly different from the actual situation. The ARIMA model had lower prediction accuracy at short time scales and higher prediction accuracy at long time scales, depending on the characteristics of the ARIMA model. As an overall linear autoregressive model, the prediction of the ARIMA model tends to become stable gradually with increasing test set time. In this study, the data volume of the 1-month time scale was larger than that of the 3-, 6-, 9-, 12-, and 24-month timescales, and the data series tended to be strictly stationary (i.e., the sequence distribution structure does not change over time). Therefore, the lowest prediction accuracy of the ARIMA model was obtained at the 1-month time scale. With increasing time scale, the amount of time-series data decreased, and the data series tended to be weakly stationary (i.e., the expectation, variance, and covariance of the stochastic process were constant; that is, the future value was related to the past value). Therefore, the fitting accuracy of the ARIMA model gradually improves as the time scale increases. At the same time, with increasing time scale, the SPI sequence obtains more information from the original sequence, and the fit of the predicted value with the actual calculated value becomes increasingly better.

Some studies have compared and analyzed the signal decomposition methods of EMD, EEMD, and CEEMD, and the analysis results show that the effects of different decomposition methods are very good, and that CEEMD can control residual noise at a relatively low level [34,44,45]. The large residual auxiliary noise of the defective EEMD cannot be avoided, and influences the experimental results. However, the influence from CEEMD can be ignored. Therefore, to stabilize the SPI sequence, in this study, CEEMD was used to extract the local features of the original sequence at different scales. Drought prediction based on CEEMD decomposition provides a stable premise for the ARIMA model. Therefore, the CEEMD–ARIMA model has a high prediction accuracy. If the data stationarity is poor, the prediction accuracy of the CEEMD–ARIMA is be reduced, as determined by the characteristics of the ARIMA model. In this case, although the prediction accuracy of the CEEMD–ARIMA model has a significant improvement over the single model, as shown in Figure 5a, Figure 6a and Figure 7a, its predictive effect is still poor. In August 2016, local rainstorms and short-term heavy precipitation occurred in the central and northern parts of Ningxia, precipitation in the southern mountainous area was rare, and the whole region presented a rare flood in the north and drought in the south. However, the combined model predicted conditions consistent with the record, indicating that CEEMD–ARIMA is suitable for use in drought research.

One limitation of this study is the diversity of factors that contribute to drought. In the arid zone of central Ningxia, the average annual precipitation is only 183.1 mm, and it is concentrated between June and August. However, with low precipitation, transpiration is high—close to 2000 mm in arid areas. Therefore, considering the influence of evapotranspiration is of great significance for accurately judging the arid zones′ drought status. In a subsequent study on the drought situation in Ningxia Hui Autonomous Region, it was necessary not only to judge the drought situation based on precipitation, but also to consider the influence of evapotranspiration. This can be achieved by dividing the study area and then selecting different drought indices to analyze the drought in each area.

6. Conclusions

In this paper, the multiscale SPI was calculated based on precipitation data from 10 stations in the Ningxia Hui Autonomous Region. Combining CEEMD in the signal processing field and the ARIMA model in the machine learning field to predict the SPI, through comparative analysis of the prediction results, the following conclusions were obtained:

(1)

As an effective nonlinear and nonstationary time-series decomposition method, CEEMD can extract the change trend of the SPI series and describe the characteristics of drought trends under climate change. Using CEEMD to decompose the SPI sequence of the Ningxia Hui Autonomous Region, seven IMF components and one trend item were obtained. The fluctuation of the component quantity became smoother than that of the original sequence, providing a basis for model prediction.

(2)

The ARIMA model had the lowest prediction accuracy on the 1-month time scale and the highest on the 24-month time scale. At the same time scales, the prediction accuracy of the CEEMD–ARIMA model was higher than that of the ARIMA model. According to the visual display of the forecast results of the 3-month time scale, in the seasons of spring, summer, autumn, and winter, the drought conditions predicted by CEEMD–ARIMA were more consistent with the actual conditions.

(3)

The drought prediction of CEEMD–ARIMA was approximately consistent with the China Meteorological Network records, indicating that the combined model is suitable for drought prediction. The original sequence was decomposed by CEEMD, and then the decomposed sequence was predicted by the ARIMA model. Finally, the predicted values of each component were added together to obtain the final prediction result. The final prediction result had high precision. According to the prediction results, the CEEMD–ARIMA model obtains higher prediction accuracy than the ARIMA model at multiple time scales, meaning that the combined model can better fit the SPI sequence at different time scales.