Drought Forecasting with Vegetation Temperature Condition Index Using ARIMA Models in the Guanzhong Plain

This paper works on the agricultural drought forecasting in the Guanzhong Plain of China using Autoregressive Integrated Moving Average (ARIMA) models based on the time series of drought monitoring results of Vegetation Temperature Condition Index (VTCI). About 90 VTCI images derived from Advanced Very High Resolution Radiometer (AVHRR) data were selected to develop the ARIMA models from the erecting stage to the maturity stage of winter wheat (early March to late May in each year at a ten-day interval) of the years from 2000 to 2009. We take the study area overlying on the administration map around the study area, and divide the study area into 17 parts where at least one weather station is located in each part. The pixels where the 17 weather stations are located are firstly chosen and studied for their fitting models, and then the best models for all pixels of the whole area are determined. According to the procedures for the models’ development, the selected best models for the 17 pixels are identified and the forecast is done with three steps. The forecasting results of the ARIMA models were compared with the monitoring ones. The results show that with reference to the categorized VTCI drought monitoring results, the categorized forecasting results of the ARIMA models are in good agreement with the monitoring ones. The categorized drought forecasting results of the ARIMA models are more severity in the northeast of the Plain in April 2009, which are in good agreements with the monitoring ones. The absolute errors of the AR(1) models are lower than the SARIMA models, both in the frequency distributions and in the statistic results. However, the ability of SARIMA models to detect the changes of the drought situation is better than the AR(1) models. These results indicate that the ARIMA models can better forecast the category and extent of droughts and can be applied to forecast droughts in the Plain.


Introduction
Earth experiences numerous extreme weather-induced disasters: droughts, floods, hurricanes, heat waves, bushfires, insect infestations and many others.Among them, drought is one of the most damaging environmental phenomena [1].Timely information about the onset of droughts, their extent, intensity, duration and impact is useful for alleviating the losses of life, human suffering and decreasing damages to economy and environment.In recent years, droughts have been occurring frequently, and their impacts are being aggravated by the rise in water demand and the variability in hydro-meteorological variables due to climate change [2].Because of this, droughts have been receiving much attention.The causes for the occurrence of droughts are complex, because they depend not only on atmospheric phenomena but also on hydrologic processes, which makes it difficult to forecast [3].Dozens of commonly used drought indices are proposed in recent decades.In general, a drought index is a prime variable for assessing the effects of a drought and defining different drought parameters, which include intensity, duration, severity and spatial extent.A drought variable should be able to quantify the drought at different time scales for which a long time series is essential.The time series of a drought index provides a framework for evaluating drought parameters of interest.A number of different indices are developed to quantify the drought, such as Standardize Precipitation Index (SPI) [4] and Palmer Drought Severity Index (PDSI) [5].Traditionally, these indices are based on data collected from the conventional ground measurements.
Satellite remote sensing provides a synoptic view of the land and a spatial context for estimating drought impacts over large areas timely and spatially [6].Many studies show that Normalized Difference Vegetation Index (NDVI) can be a useful index for studying vegetation and ecosystems in semi-arid environments where vegetation cover is less than 30% [7][8][9].Significant relationships between time series of NDVI and various vegetation indicators including green leaf area index, green biomass production, as well as rainfall or soil moisture in semi-arid environments have been reported [10][11][12][13][14][15].Consequently, Advanced Very High Resolution Radiometer (AVHRR) instruments on board National Oceanic and Atmospheric Administration (NOAA) derived NDVI and other related indices have been successfully used to identify and monitor areas affected by droughts at regional and local scales [16][17][18][19][20]. Tucker and Choudhury [17] found that NDVI could be used as a response variable to identify and quantify drought disturbance in semi-arid and arid lands since its low values correspond to stressed vegetation.Furthermore, Ji and Peters [20] found that NDVI is an effective indicator of vegetation response to drought in the Great Plains of the USA.Since NDVI has been proved to represent vegetation responses in a timely manner to climate variability as a normalized ratio, the Vegetation Condition Index (VCI), which is NDVI normalization for each pixel based on minimum and maximum NDVI values over time, was developed by Kogan [21] in order to relatively assess changes in the NDVI signal through time by reducing the influence of local ecosystem variables.In addition, the probability of the standardized NDVI anomaly, the Standardized Vegetation Index (SVI), has been used to monitor areas affected by droughts and vegetation conditions in terms of relative greenness at pixel level over time periods [22].Gutman [23] showed that the thermal data from polar orbiters might be useful for detecting the inter-annual changes in surface moisture, when the change in the vegetation index fails to produce a significant signal.However, NDVI has two main limitations for drought monitoring, first the apparent time lag between rainfall and the NDVI response and second, little influence of significant precipitation events later in the growing season (plant seed production period) on NDVI.Several authors used thermal (e.g., land surface temperature (LST)) products of the NOAA/AVHRR and MODIS to provide a more ecological and physical interpretation of remotely sensed data for examining drought conditions [24,25].
The combined application of the NDVI and LST for mapping and monitoring land surface provides a more mechanistic interpretation of the remotely sensed data.Goetz [26] reported that the negative correlation between LST and NDVI observed at several scales (25 m 2 to 1.2 km 2 ), was largely due to changes in vegetation cover and soil moisture, and indicted that the surface temperature can rise rapidly with water stress.Nemani et al [27] found the slope of LST versus NDVI to be negatively correlated to a crop-moisture index.Therefore, the ratio of LST/NDVI increases during times of drought.There are two methods currently being put forward.The first is a progression from the slop of LST versus NDVI approach, which describes the data as falling into a triangle [28][29][30][31].The second, the Vegetation Index/Temperature Trapezoid (VITT) [32,33], is an evolution of the Crop Water Stress Index (CWSI) [34], which promotes the ideas of data falling into a trapezoid.Based on AVHRR brightness temperature, temperature condition index (TCI) have been developed and used for monitoring drought [35].On the basis of the triangular space of LST and NDVI, the vegetation temperature condition index (VTCI) approach was proposed for monitoring drought occurrence at a regional level [36], and the temperature vegetation dryness index (TVDI) was developed for assessing soil surface moisture status [37].VTCI is time-dependent and usually region-specific, and is better used during plant growing seasons.It was determined that the VTCI can be used to monitor droughts in a near real time at regional level, to better indicate the severity of droughts in the Guanzhong Plain of China and the Great Plains of the USA [24,25].Lin et al. [38] validated that VTCI drought monitoring method at scale of 10 days is more suitable.Zhang et al. [39] classified and validated drought categories for the VTCI.Wang et al. [40] estimated soil water in northern China based on VTCI and determined that VTCI method has good precision and can monitor the drought conditions in northern China more accurately.
Drought forecasting is very important for effectively assisting local authorities to mitigate the drought impacts, and to reasonably use water resources.Nowadays, the time series forecasting has been widely applied and has become an important approach of drought forecasting.In a time series model, the past observations are analyzed to formulate a model describing the inner correlation among them.Then, the time series is extrapolated into the future according to the model.The Markov chain approach was used by Lohani and Loganathan [41] to develop an early warning tool.A non-homogeneous Markov chain formulation was adopted to derive drought characteristics and assess dry spells from long-term records of the PDSI in two climatic areas of Virginia, USA.Rao and Padmanahan [42] investigated the stochastic characteristics of yearly and monthly Palmer's Drought Index (PDI) to characterize those using valid stochastic models to forecast and simulate PDI series.Sen [43] predicted the possible critical drought durations that may result from any hydrologic phenomena during any future periods using the second order Markov chain.Kim and Valdes [44] used the PDSI as drought parameter to forecast droughts in the Conchos River Basin of Mexico.
Another most widely used time series model is the AutoRegressive Integrated Moving Average (ARIMA) model [45].However, it is basically a linear model assuming that time series data are stationary, and have a limited ability to capture non-stationarities and nonlinearities in the data.The main advantage of ARIMA model forecasting is that it only requires the time series data.Thus, in time series analyses of transportation freight and transportation demand, ARIMA models are mostly used [46].ARIMA modeling has been also implemented in runoff and inflow [47,48], electric grid [49].
ARIMA models effectively consider serial linear correlation among observations, whereas Seasonal AutoRegressive Integrated Moving Average (SARIMA) models can satisfactorily describe time series that exhibits non-stationarity both within and across seasons.The SARIMA models describe the seasonality and autocorrelation structure of the time series of VTCI more complexly and appear to be more suitable for assessing the relations within the series.In the AR(1) models' development [50], there is only one parameter to be estimated.If the estimated parameter is larger than the one, the forecasting results will be larger than the monitoring ones, and with the increase of the forecasting steps, the forecasting values are gradually increased.The SARIMA modeling approach has been promoted as the statistical method of choice in the case of data arising from observations collected over long periods of time.The main objective of the present study is to model time series of VTCI for drought forecasting by applying SARIMA models in the Guanzhong Plain of China.

Study Area
The study was conducted in the Guanzhong Plain, the study area having central geographic coordinates 34.65 • N and 108.75 • E, Shaanxi Province, China, which has low atmospheric water vapor content, low aerosol loading, and low probability of cloudiness to minimize atmospheric variation.The area is susceptible to spring and early summer droughts where the variability of annual precipitation is high.This study area has flat terra, fertile soil and relatively better natural conditions compared to other farmlands in the province.The cultivate index is about seventy percent and most of the study area is covered by winter wheat in spring and early summer.The study area includes irrigated farmlands as well as rainfed farmlands, which makes it suitable for using the VTCI approach for drought monitoring and forecasting.The drought forecasting for this area is very important to mitigate the effects of droughts and to enable water resource management.

Composition of the NDVI and LST products
Data collected from AVHRR instruments onboard National Oceanic and Atmospheric Administration (NOAA) satellites were used.The AVHRR remotely sensed data acquired during the daytime at nadir and near-nadir views under clear-sky conditions were selected for the study.The maximum value compositing (MVC) approach was applied to generate the 10-day maximum value composite NDVI and LST products for each year at 10-day intervals.We also used MVC approach to generate the 10-day NDVI and LST products for the study years at the given period; these are known as the multiyear maximum value composite NDVI and LST products.When combining the NDVI and LST, the purpose of using the MVC technique was to reduce the effects of the solar zenith angle, the satellite's view angle, the orbit drift and cloud contamination.These effects cannot be minimized by the minimum value compositing approach.The minimum value compositing approach is used to generate the multiyear minimum value composite LST products using each year's maximum value composite LST products for the 10 days; this is called the multiyear maximum-minimum value composite LST product [25].

VTCI calculations
Vegetation temperature condition index is defined as [24]: where LST NDV I i min and LST NDV I i max are the minimum LST and maximum LST of the pixels that have the same NDV I i value in a study region, respectively, and LST NDV I i denotes the LST of one pixel whose NDVI is NDV I i , Coefficients a, b, a and b can be determined from an area large enough where the soil moisture in the surface layer should span from wilting point to field capacity at the pixel level.The coefficients of Equations ( 2) and (3) were estimated from the scatter plots of NDVI and LST using different composite NDVI and LST products.
VTCI can be physically interpreted as the ratio of LST differences among the pixels with a specific NDVI value in an area large enough to provide wide ranges of NDVI and soil moisture at surface layer.The VTCI drought monitoring approach integrates the remotely sensed land surface reflectance and thermal properties, and gives the emphasis on changes in both LST and NDVI over a region.The VTCI drought monitoring results at different time scales show that the method at the scale of 10 days is more accurate and practical [35].The time series of VTCI in the Guanzhong Plain is processed at 10-day intervals in this study.

SARIMA Model
ARIMA models are one of the classes of stochastic models for describing time series.When the process remains in equilibrium about a constant mean level, which is so called stationary, the AutoRegressive (AR) models, Moving Average (MA) models or mixed AutoRegressive Moving Average (ARMA) models are implied.When the time series exhibit non-stationary behavior and in particular do not vary about a fixed mean, the ARIMA and SARIMA models are implied.
In this study, we mainly discuss the SARIMA models on our VTCI dataset at the ten-day intervals.The seasonal effect indicates that an observation of a particular ten-day, for example the first ten days of April 2009, is related to the observations of the first ten days of April 2008.The form of a SARIMA model is where z t is the observed series with period s, and α t is the white noise process.The subscripts p, P, q and Q are the orders of the models, and the resulting multiplicative process is known to be of order (p, d, q) (P, D, Q) s , whereas d and D are the degrees of simple differencing and seasonal differencing, and B is the backward shift operator which is defined by Bz t = z t−1 .Hence, B s z t = z t−s .∇ is another important operator of backward difference operator, which can be written in terms of B, since In turn, ∇ has for its inverse the summation operator s given by ∇ and the moving average operator of q is defined by ) and Θ (B s ) are operators in B s of degrees of P and Q.

SARIMA Modeling Procedure
A SARIMA model is usually best achieved by a three stage iterative procedure based on identification, estimation, and diagnostic checking [45].Identification means the use of data, and any information on how the series is generated; estimation means the efficient use of the data to make inferences about parameters conditional on the adequacy of the entertained model; and diagnostic checking means checking the fitted model in its relation to the data with an intention to reveal any inadequacy of the model and to perform a model improvement.

Identification
The autocorrelation and partial autocorrelation functions are often used in the identification.The autocorrelation function (ACF) of series z t at lag k (ρ k ) is where µ is the mean value of the series, and the variance σ 2 z of the stochastic process can be estimated by The partial autocorrelation function (PACF) of series z t at lag k (φ kk ) is where The ACF is used to identify the degree of differencing in the time series.In the case of a stationary model, the ACF will quickly "die out" for moderate and large k.Generally, if the ACF tails off, and is significantly different from zero after the lag 3, the time series is considered non-stationarity, and the periodic components reflect in very large correlations at lags 1s, 2s and 3s.Conversely, the time series is deemed to be stationarity, if the time series contains both general and periodic components, making simple differencing and seasonal differencing for z t as many times it needs to produce stationarity.
Having tentatively decided d and D, we further identify the orders of the SARIMA process and the principal tools of ACF and PACF as well.Studying the general appearances of the ACF and PACF of the appropriately differenced series provides clues about the choice of the orders of p, P, q and Q for the autoregressive along with moving average operators.To do so, recall the characteristic behaviors of theoretical ACF and the theoretical PACF for the moving average, autoregressive and mixed processes.
Briefly, whereas the ACF of an autoregressive process of order p tails off, its PACF has a cutoff after lag p. Conversely, the ACF of a moving average process of order q has a cutoff after lag q, while its PACF tails off.If both the ACF and PACF tail off, a mixed process is suggested [51].The identification of P and Q are similar to p and q, just consider 1s, 2s and more integral period lags.

Estimation
For a complete understanding of the estimation situation, it is necessary to make a thorough analytical and graphical study of the likelihood function, or in the Bayesian framework, the posterior distribution of the parameters, which, in the situations we consider, are dominated by the likelihood.
In the examples thus far, for moderate and large samples, the log-likelihood function will be unimodal and could be adequately approximated over a sufficiently extensive region near the maximum by a quadratic function.In such cases, the log-likelihood function can be described by its maximum and its second derivatives' at the maximum.The values of the parameters that maximize the likelihood function, or equivalently the log-likelihood function, are called maximum likelihood estimates.The second derivatives of the log likelihood function provide measures of "spread" of the likelihood function and can be used to calculate approximate standard errors for the estimates.
In addition, as the ACF and the PACF determine more than one model, the Akaike Information Criteria (AIC) is used to identify the best fitted model among them [52,53].
where ω is the number of estimated parameters, and β is the maximum likelihood function values.AIC consists of two parts.The first item reflects the model precision.It decreases along with the increase of the order number.The second marks the number of model parameters, which presents positive relation with the order number.These indicate that the model is a compromise with the precision and the number of parameters.

Diagnostic Checking
The model has been identified and the parameters estimated; the diagnostic check is then applied to the fitted model, as the diagnostic check is based on the analysis of residuals.We refer to the quantities as residuals.When the parameters are well estimated, the tentative model accuracy is validated by examining the autocorrelations ρ k ( α) of the residuals in order to simulate the white noise process.Furthermore, the Q-statistics test is applied to verify the tentative adequateness of the model.
where m is the specified delay lags, and n is the length of the residuals.If the calculated value of Q exceeds the critical value of χ 2 with m degrees of freedom obtained from the chi-square tables, the tentative model is tuned as inadequate; otherwise, the model is adequate.

Data Processing
At regional scale there are abundant of literatures documenting the use of NOAA/AVHRR data archives to perform time series analysis in regions all over the world and addressing applicability in wide and vibrant range.The remotely sensed NOAA/AVHRR's data at 1.1 km geometric resolution are used in our study.The acquired time series is from March to May of the years 2000 to 2009.The AVHRR's remotely sensed data acquired at daytime at nadir and near nadir views under clear sky conditions are selected for the study.The atmospheric correction, radiometric calibration and geometric correction of AVHRR's data are carried out using the methods described by Wan and Sun [24,25].The time series of VTCI in the Guanzhong Plain is processed at ten-day intervals (e.g., the first ten days of March as from 1 March to 10 March of each year) using our proposed approaches [24,25].Due to the cloud contamination at some ten-day intervals, there are some missing data in the time series of VTCI.In order to obtain an equally spaced time series, it is necessary to complement the missed ones.Thus far, the cubic spline interpolation method is chosen [54] Considering the study area is overlaid on the administration map around the study area before modeling the time series, divide the study area into 17 parts where at least one weather station is located in each part, as shown in Figure 1.The grey shading in Figure 1 is the study area, and the closed lines are borders of the administration counties.The selected 17 weather stations have a relatively uniform distribution in the study area.The pixels where the 17 weather stations are located are chosen at first and studied for their fitting models, and then considered the best models for all pixels of the whole area to be determined.According to the above-mentioned procedures for the models' development, the selected best models for the 17 pixels are identified and forecasted in three steps.By comparing the forecasting VTCI values with the monitoring ones (original VTCI drought monitoring results), the absolute errors (the difference between the forecasting results and the monitoring results) are calculated pixel by pixel for each forecasting image and the frequency distributions of absolute errors are the distributions of the absolute errors for each forecasting image.The average, root mean square error (RMSE) and histogram of the absolute errors are used for the accuracy analysis.The RMSE of the pair wise differences of the two datasets can serve as a measure of how far the average error is from zero, and is defined as: where x is the forecasting results, x is the monitoring results, and n is the number of pixels.

Results
The As can be seen in Figure 2, the series contains non-stationary characteristics, since it represents "uncontrolled" behaviors.The dataset is formed by ten years with nine data in each year, which indicates the periodicity, as s equals to nine, as shown in Figure 2a.The estimated ACF of the series indicates that at least one differencing is necessary in Figure 2b.Therefore, considering first simple differencing (∇) and first seasonal differencing (∇ 9 ) of z t , the ACF and PACF of ∇∇ 9 z t are plotted in Figure 3.The ACF of ∇∇ 9 z t are small after the second lag, which shows the series ∇∇ 9 z t is stationary, and the seasonal and trend components within the time series of VTCI are eliminated.The next step is to identify the values of p, P, q and Q.The ACF and PACF both tail off exponentially in Figure 3, which indicates a combined autoregressive moving average (ARIMA) process that contains a pth order autoregressive component and a qth order moving average component.p is smaller than 3 (Figure 3b), q is smaller than 2 (Figure 3a), and P and Q are both equal to zero.In order to get the fitting model among the different combination models, the AIC values of each model for Tongchuan weather station are calculated using Equation ( 6), as shown in Table 1.The ACF and PACF both tail off exponentially in Figure 3, which indicates a combined autoregressive moving average (ARIMA) process that contains a pth order autoregressive component and a qth order moving average component.p is smaller than 3 (Figure 3b), q is smaller than 2 (Figure 3a), and P and Q are both equal to zero.In order to get the fitting model among the different combination models, the AIC values of each model for Tongchuan weather station are calculated using Equation ( 6), as shown in Table 1.
In Table 1, the values of AIC reach their minimum when p = 1 and q = 1, and therefore, ARIMA(1, 1, 1)(0, 1, 0) 9 is considered the suitable model for Tongchuan weather station.After the estimation steps, the diagnostic check is necessary.The autocorrelations of the ARIMA(1, 1, 1)(0, 1, 0) 9 model residuals are shown in Table 2.The diagnostic check is provided by the quantity Q test (Equation ( 7)), where n = 84, m = 25, which is approximately distributed as χ 2 with 25 degrees of freedom.The observed value of Q = 84 × 0.2314 = 19.4 is smaller than χ 2 0.95 (25) = 37.6, and the residual is white noise.The described diagnostic check shows that the model is adequate.
The same process is applied to determine the parameters of VTCI dataset in the other 16 weather stations.The ACF of the datasets does not die out rapidly, suggesting that the time series of VTCIs of the dataset are non-stationary.Therefore, the first simple differencing and first seasonal differencing are used to process the datasets to stationary as ∇∇ 9 z t .From the ACFs and PACFs of the 16 time series of VTCIs, the values of P and Q are both equal to zero, while the values of p and q of the models vary, as seen in Table 3.Most values of p are equal to one or two, and values of q are all equal to one except for one model in which q equals to zero.Although only the model structures of Baoji (ARIMA(2, 1, 0)(0, 1, 0) 9 ) and Liquan (ARIMA(0, 1, 1)(0, 1, 0) 9 ) areas in particular are mentioned in Table 3, observe that the models of the 17 parts can be mainly divided into two model types: ARIMA(1, 1, 1) (0, 1, 0) 9 and ARIMA(2, 1, 1)(0, 1, 0) 9 .Integration of ground observed data in the entire area as well as in the specified areas shows that most of the parts identified in ARIMA(2, 1, 1)(0, 1, 0) 9 models are irrigated farmlands, while the ARIMA(1, 1, 1)(0, 1, 0) 9 types are rainfed farmlands.These results show that the irrigated farmlands of the current ten-day conditions are associated with the past first and second ten-day intervals, while in the rainfed farmlands, the current ten-day conditions are only relative with the previous ten-day interval.We can see the whole area is relatively uniform.Furthermore, the ARIMA(2, 1, 1)(0, 1, 0) 9 models are gathered in the middle part of the Guanzhong Plain (Table 3 and Figure 1), and ARIMA(1, 1, 1)(0, 1, 0) 9 models are in the east part of the Plain (Table 3 and Figure 1).Baishui is the only mountain county in Weinan City, differing from its surrounding counties, which leads to the different model structure.In this respect, the identified models are reasonable for the Plain.
Regarding the values of p and q identified in each part of Table 3, the SARIMA approach has been applied to forecast droughts in the whole area of the Guanzhong Plain.Our former developed AR(1) models were also used for the forecasting [50].Employing the results shown in Figure 4, the forecasting and monitoring results for the ten days periods of the SARIMA and the AR(1) models are compared in order to determine the better performing model.The forecasting results of the SARIMA models show that the east of the Guanzhong Plain has drought conditions in the second and last ten days of April 2009, and the VTCI's values in the west of the study area are higher than in the east, which shows that the drought conditions are not severe in the west of the Plain.The AR(1) models results show that the values of VTCI are relatively high and uniform compared to the results of the SARIMA models in the whole study area, and from the Figure 4 (Figure 4g-i) we cannot find obvious drought conditions from the ten days of April, 2009.In general, the monitoring results (Figure 4) show that there is drought occurrence in the northeast of the Plain in the middle and last ten days of April 2009.The forecasting results of the SARIMA models are in good agreement with the monitoring ones up to a certain extent, while those of the AR(1) models approaches are not in good agreement with the actual drought conditions.

Evaluating the Models Using Statistical Analysis
The frequency distributions of the absolute errors of the SARIMA models and the AR(1) models in Figure 5 are nearly symmetric, and the ranges of the absolute errors are increased with the

Evaluating the Models Using Statistical Analysis
The frequency distributions of the absolute errors of the SARIMA models and the AR(1) models in Figure 5 are nearly symmetric, and the ranges of the absolute errors are increased with the increase of the forecasting time periods for the two kind of models.The absolute errors of the SARIMA models and the AR(1) models of the first step prediction are mainly distributed in the range from −0.1 to 0.1, and the ranges of the two models are nearly the same.However, the ranges of the second and third step predictions for the SARIMA models are wider than those of the AR(1) models.Considering this, the precisions of the AR(1) models might be better than those of the SARIMA models, especially for the predictions of the second and the third steps.The averages, RMSEs and ranges of the absolute errors between the forecasting VTCIs and actual monitoring ones are calculated and used for accuracy analysis developed models, as listed in Table 4.The averages and RMSEs of the first step prediction (first ten days of April 2009) of the SARIMA models and the AR(1) models are 0.0243 and 0.0818 and −0.0047 and 0.0513, and those of the second step (middle ten days of April 2009) prediction and the third step (last ten days of April 2009) prediction are 0.1194 and 0.1745 and 0.0412 and 0.0697, and 0.0961 and 0.2004 and 0.0730 and 0.1414, respectively.For both the SARIMA models and the AR(1) models, the errors of the first step prediction is the lowest and the errors are increased with the increase of the forecasting time periods.These results indicate the forecasting abilities of the models decrease with the increase of the forecasting steps.In addition, comparing the AR(1) and SARIMA models, we can find that the absolute errors of the AR(1) models are lower than the SARIMA models, both in the frequency distributions figures and in the statistic results.The principles of the AR(1) models, which predict VTCI values that are close to the mean values of the time series, they also leads to the small deviations from the monitoring ones in most pixels.However, in some particular areas, the AR(1) models cannot detect the changes of the drought situation, and we can evaluate the two models in another aspect, as discussed in the following section.

Evaluating the Models Using Drought Categories
The proposed VTCI is a quantitative drought monitoring approach has two ways to emphasize drought studies: first, a drought monitoring index that is defined quantitatively, and second an index that is based on categorizes, such as normal or wet, and mild, moderate, severe and extreme droughts.The purpose of drought forecasting is to detect droughts and abnormality; if the forecasting results do not show drought conditions, they are of little significance.According to our research on VTCI drought categories in Table 5 [41], the VTCI's forecasting and monitoring has three periods, as categorized in Figure 6.The category distribution and extent of droughts are useful to evaluate the drought forecasting accuracies.Based on the categorized drought monitoring results in Figure 6a-c, most of the Guanzhong Plain has no drought conditions (normal or wet categorized in Table 4) in the first ten days of April 2009, and mild droughts appear in the northeast of the Plain in the second ten days of April 2009.Furthermore, the drought situation of the northeast is more severe in the last ten days of April 2009.The categorized drought forecasting results of the SARIMA models in Figure 6d-f coincide in coordination with the monitoring ones, as drought conditions in the northeast becomes more severe.Most of the categorized drought forecasting results of the AR(1) models in Figure 6g-i show no drought grades in the first, second and last ten days of April 2009.Based on the analysis, it is concluded that the forecasting results of the SARIMA models could better forecast the category and extent of droughts.
Severe drought <0.37 Based on the categorized drought monitoring results in Figure 6a-c, most of the Guanzhong Plain has no drought conditions (normal or wet categorized in Table 4) in the first ten days of April 2009, and mild droughts appear in the northeast of the Plain in the second ten days of April 2009.Furthermore, the drought situation of the northeast is more severe in the last ten days of April 2009.The categorized drought forecasting results of the SARIMA models in Figure 6d-f coincide in coordination with the monitoring ones, as drought conditions in the northeast becomes more severe.Most of the categorized drought forecasting results of the AR(1) models in Figure 6g-i show no drought grades in the first, second and last ten days of April 2009.Based on the analysis, it is concluded that the forecasting results of the SARIMA models could better forecast the category and extent of droughts.Moreover, the categorized errors and their statistical results in Figure 7 show that both SARIMA models and AR(1) models have good ability for the prediction of no drought conditions.However, the key point of drought forecasting is to predict the drought conditions, and no drought conditions are not of particular concern.From this perspective, the SARIMA models are more reliable than AR(1) models.In the first forecasting step, as seen in Figure 7a-d, most of the drought categories are correct in the first ten days of April 2009 (the categorized error pixels of the SARIMA models and the AR(1) models are 4.856% and 1.723%, respectively) as most of the study area of the monitoring image has no drought conditions (Figure 7a).For the second step forecasting (Figure 7b,e), numerous forecasting drought categorized errors of the SARIMA models and the AR(1) models in the second ten days of April 2009 start to appear, and the categorized error pixels of the SARIMA models and the AR(1) models are 25.293% and 21.025%, respectively.The SARIMA models forecast has no drought conditions to mild droughts, and mild droughts to moderate droughts (Figure 7b) in the second ten days of April 2009, while the AR(1) models forecast almost all kinds of the drought grades to no drought conditions (Figure 7e).These results indicate that the AR(1) models cannot predict the drought categories accurately in the second step forecasting.For the third step forecasting (Figure 7c,f), more categorized errors of the SARIMA models and the AR(1) models in the last ten days of April 2009 appear, and the categorized error pixels of the SARIMA models and the AR(1) models are 48.424% and 46.444%, respectively.Figure 7c shows various error conditions.The total error pixels of categorizing other drought grades to no drought grades are 29.041%, while in Figure 7f the total numbers of error pixels are 45.910%.The so-called categorizing errors of the AR(1) models are mainly reflected in categorizing various drought levels into no drought grades.Only 0.533% pixels are other types of categorizing errors, in which the drought types of prediction results have less severity than those of the monitoring ones.From this analysis, we came to the conclusion that the AR(1) models have low ability to forecast drought conditions in the third step forecasting, and the SARIMA models are better to forecast drought conditions.Considering the three steps forecasting, the performance of the SARIMA models is very promising.Moreover, the categorized errors and their statistical results in Figure 7 show that both SARIMA models and AR(1) models have good ability for the prediction of no drought conditions.However, the key point of drought forecasting is to predict the drought conditions, and no drought conditions are not of particular concern.From this perspective, the SARIMA models are more reliable than AR(1) models.In the first forecasting step, as seen in Figure 7a-d, most of the drought categories are correct in the first ten days of April 2009 (the categorized error pixels of the SARIMA models and the AR(1) models are 4.856% and 1.723%, respectively) as most of the study area of the monitoring image has no drought conditions (Figure 7a).For the second step forecasting (Figure 7b,e), numerous forecasting drought categorized errors of the SARIMA models and the AR(1) models in the second ten days of April 2009 start to appear, and the categorized error pixels of the SARIMA models and the AR(1) models are 25.293% and 21.025%, respectively.The SARIMA models forecast has no drought conditions to mild droughts, and mild droughts to moderate droughts (Figure 7b) in the second ten days of April 2009, while the AR(1) models forecast almost all kinds of the drought grades to no drought conditions (Figure 7e).These results indicate that the AR(1) models cannot predict the drought categories accurately in the second step forecasting.For the third step forecasting (Figure 7c,f), more categorized errors of the SARIMA models and the AR(1) models in the last ten days of April 2009 appear, and the categorized error pixels of the SARIMA models and the AR(1) models are 48.424% and 46.444%, respectively.Figure 7c shows various error conditions.The total error pixels of categorizing other drought grades to no drought grades are 29.041%, while in Figure 7f the total numbers of error pixels are 45.910%.The so-called categorizing errors of the AR(1) models are mainly reflected in categorizing various drought levels into no drought grades.Only 0.533% pixels are other types of categorizing errors, in which the drought types of prediction results have less severity than those of the monitoring ones.From this analysis, we came to the conclusion that the AR(1) models have low ability to forecast drought conditions in the third step forecasting, and the SARIMA models are better to forecast drought conditions.Considering the three steps forecasting, the performance of the SARIMA models is very promising.

Conclusions
The Guanzhong Plain is one of the areas subjected to droughts.In this study, the seasonal ARIMA models are developed to forecast droughts using the quantitative drought monitoring index called VTCI, and compared with Han's AR(1) models for better performance of the models.The

Conclusions
The Guanzhong Plain is one of the areas subjected to droughts.In this study, the seasonal ARIMA models are developed to forecast droughts using the quantitative drought monitoring index called VTCI, and compared with Han's AR(1) models for better performance of the models.The structures of the models are identified in 17 pixels of the selected area where the weather stations are located, and the results of model identification show that the whole area can be divided into two types of model structures, which are related to irrigated farmlands and rainfed farmlands.Moreover, the absolute errors of the AR(1) models are less than those of the SARIMA models, while the forecasted drought conditions of the AR(1) models are not in good agreement with the monitoring ones.Based on the categorized drought forecasting results, the AR(1) models can predict no drought grades very well, and the SARIMA models can not only predict the no drought grades but also the drought grades, which indicates that the performance of the SARIMA models is better than that of the AR(1) models.Result shows that the selected SARIMA models are appropriate for drought forecasting.
Further, the SARIMA models describe the seasonality and autocorrelation structure of the time series of VTCI more complexly and appear to be more suitable for assessing the relations within the series.The SARIMA modeling approach has been promoted as the statistical method of choice in the case of data arising from observations collected over long periods of time.Compared with the AR(1) modeling approach, the SARIMA modeling approach has several advantages, in particular, its forecasting capability and its richer information on time-related changes.The VTCI time series of this study are constructed by 10 years, and has seasonal variations and fluctuations, Overall, the SARIMA models are better than the AR(1) models.The forecasting abilities of the SARIMA models decrease with the increase of the forecasting steps.
This research work contributes in developing a feasible approach for drought forecasting in the Guanzhong Plain.Results demonstrate that the SARIMA models can be used for drought forecasting; however, the absolute errors of the SARIMA models are still large and challenging.Further studies to improve the precision of the models by employing other methods to identify seductive models are being considered.
. After the processing, about 90 VTCI images are acquired at the ten-day intervals in the Guanzhong Plain from the erecting stage to the maturity stage of winter wheat during 2000-2009.The time series of VTCI from the first ten days of March 2000 to the last ten days of March 2009 is used for the model development, and the dataset from the first ten days of April to the last ten days of April 2009 is used for the model validation.
acquired time series is from the first ten days of March to the last ten days of May during 2000-2009.The dataset from the first ten days of March of the year 2000 to the last ten days of March of the year 2009 was used for model development and for obtaining the fitted model for each weather station.First, start with model and analyze the time series of VTCI in each weather station, and take Tongchuan station as an example to present the subsequent modeling process.The graph and ACF of the time series VTCI in Tongchuan station are shown in Figure 2.

Figure 2 .Figure 3 .
Figure 2. The time series of VTCI in Tongchuan weather station (a); and its autocorrelation function (b).

Figure 3 .
Figure 3.The autocorrelation (a) and partial autocorrelation functions (b) at differenced time series of VTCI in Tongchuan weather station.

Figure 4 .
Figure 4.The drought forecasting and monitoring results of April 2009: (a-c) the monitoring results in the first, second and last ten days of April 2009, respectively; (d-f) the forecasting results of the SARIMA models; and (g-i) the forecasting results of the AR(1) models.

Figure 4 .
Figure 4.The drought forecasting and monitoring results of April 2009: (a-c) the monitoring results in the first, second and last ten days of April 2009, respectively; (d-f) the forecasting results of the SARIMA models; and (g-i) the forecasting results of the AR(1) models.

Figure 5 .
Figure 5.The absolute errors of the forecasting results of the SARIMA models (a) and the AR(1) models (b), and their frequency distributions of the first, second and last ten days of April 2009.

Figure 5 .
Figure 5.The absolute errors of the forecasting results of the SARIMA models (a) and the AR(1) models (b), and their frequency distributions of the first, second and last ten days of April 2009.

Figure 6 .
Figure 6.The drought forecasting and monitoring categories results of April 2009: (a-c) the monitoring categories in the first, second and last ten days of April 2009, respectively; (d-f) the forecasting categories of the SARIMA models; and (g-i) the forecasting categories of the AR(1) models.

Figure 6 .
Figure 6.The drought forecasting and monitoring categories results of April 2009: (a-c) the monitoring categories in the first, second and last ten days of April 2009, respectively; (d-f) the forecasting categories of the SARIMA models; and (g-i) the forecasting categories of the AR(1) models.

Figure 7 .
Figure 7.The images of categorized errors and the statistical results that are calculated by the differences of the categorized drought monitoring and forecasting results: (a-c) the results of the SARIMA models in the first, second and last ten days of April 2009, respectively; and (d-f) the results of the AR(1) models

Figure 7 .
Figure 7.The images of categorized errors and the statistical results that are calculated by the differences of the categorized drought monitoring and forecasting results: (a-c) the results of the SARIMA models in the first, second and last ten days of April 2009, respectively; and (d-f) the results of the AR(1) models

Table 1 .
The AIC values of the identified SARIMA models of Tongchuan weather station.

Table 2 .
Estimated autocorrelations of the fitting model residuals of Tongchuan weather station.

Table 3 .
The values of p and q of the time series of VTCIs in the 17 weather stations.

Table 4 .
The statistics of the absolute errors between the forecasting VTCIs and the monitoring ones.

Table 5 .
Vegetation temperature condition index drought categories.