The Temperature Vegetation Dryness Index (TVDI) Based on Bi-Parabolic NDVI-T_s Space and Gradient-Based Structural Similarity (GSSIM) for Long-Term Drought Assessment Across Shaanxi Province, China (2000–2016)

Abstract

Traditional NDVI-T_s space is triangular or trapezoidal, but Liu et al. (2015) discovered that the NDVI-T_s space was bi-parabolic when the study area was covered with low biomass vegetation. Moreover, the numerical value of the indicator was considered in most of the study when the drought conditions in the space domain were evaluated. In addition, quantitatively assessing the spatial-temporal changes of the drought was not enough. In this study, first, we used MODIS NDVI and T_s data to reexamine if the NDVI-T_s space with “time” and a single pixel domain is bi-parabolic in the Shaanxi province of China, which is vegetated with low biomass to high biomass. This is compared with the triangular NDVI-T_s space and one of the well-known drought indexes called the temperature-vegetation index (TVX). The results demonstrated that dry and wet edges exhibited a parabolic shape again in scatter plots of T_s and NDVI in the Shaanxi province, which was linear in the triangular NDVI-T_s space. The Temperature Vegetation Dryness Index (TVDI_c) was obtained from bi-parabolic NDVI-T_s andTVDI_t was obtained from the triangular NDVI-T_s space and TVX were compared with 10-cm depth relative soil moisture. By estimating the 10-cm depth soil moisture, TVDI_c was better than TVDI_t, which were all apparently better than TVX. Second, combined with MODIS data, the drought conditions of the study area were assessed by TVDI_c between 2000 to 2016. Spatially, the drought in the Shaanxi Province between 2000 to 2016 were mainly distributed in the northwest, North Shaanxi, and the North and East Guanzhong plain. The drought area of the Shaanxi province accounted for 31.95% in 2000 and 27.65% in 2016, respectively. Third, we quantitatively evaluated the variation of the drought status by using Gradient-based Structural Similarity (GSSIM) methods. The area of the drought conditions significantly changed and moderately changed at 5.34% and 40.22%, respectively, between 2000 and 2016. Finally, the possible reasons for drought change were discussed. The change of precipitation, temperature, irrigation, destruction or betterment of vegetation, and the enlargement of opening mining, etc., can lead to the variations of drought.

1. Introduction

Drought is a natural phenomenon that entails a lack of water. It is commonly divided into three types including meteorological drought (lack of precipitation), hydrological drought (decreasing of river flow, lake or reservoir capacity, and groundwater level), and agricultural drought (shortage of available soil moisture for the plant) [1,2]. The drought in this context is agricultural drought. As a popular and dangerous natural disaster, it can result in many hazards for farming, which makes the environment worse and potentially causes other natural disasters to happen. As such, it is very important to assess the spatial-temporal variations of droughts [3].

Remote sensing techniques can quickly monitor drought conditions in a large area [4]. Since the 1970s, many methods based on the vegetation index, land surface temperature (T_s), and albedo such as the Normalized Difference Vegetation Index (NDVI) [5], Perpendicular Vegetation Index (PVI) [6], the Vegetation Condition Index (VCI) [7,8,9,10], the Normalized Difference Water Index (NDWI) [11], Vegetation Condition Albedo Drought Index (VCADI) [12], and the Temperature Condition Index (TCI) [7,13] have been applied to evaluate the vegetation status and drought conditions. Most studies have pointed out that it is insufficient to monitor drought conditions with a single index. T_s and vegetation information have been combined for assessing droughts [1,2,4,14,15].

Goward et al. (1985), Hope et al. (1986), and Nemani et al. (1989) found that the slope of the T_s/NDVI curve could provide some soil moisture information [16,17,18]. In the study of Price (1990) and Carlson et al. (1994), a triangle was shaped in the scatterplots of NDVI and T_s [19,20]. A trapezoid was explored in scatterplots of NDVI and T_s in the research of Goward and Hope (1989) and Moran et al. (1994) [21,22]. The temperature vegetation dryness index (TVDI) [23] and the vegetation temperature condition index (VTCI) [24] obtained from the triangular NDVI-T_s space were developed to estimate drought or soil moisture conditions. The principle of VTCI is almost the same as TVDI and VTCI = 1-TVDI. These two indexes were only acquired from remotely sensed data and are widely applied to evaluate the drought status and assess its impacts [25,26,27,28,29,30].

In the triangular or trapezoidal NDVI-T_s space, the dry and wet edges are linear. However, the dry edge was not always linear. The dry edge was shaped in a parabola in the triangular NDVI-T_s space [19]. When the NDVI values were very low or high, the dry edge was demonstrated in nonlinear in the study of Carlson et al. (1994) [20]. A fourth-order polynomial was used to fit the dry edge in the NDVI-T_s space [31]. A bi-parabola was found in the NDVI-T_s space when the study area was covered with low biomass vegetation [32]. Moreover, most of studies only considered the numerical value of the indicator when evaluating drought conditions in the space domain, but quantitatively assessing the spatial-temporal changes of drought was not enough. Structural similarity (SSIM) was used to quantitatively and accurately describe the spatial distribution characteristics of the drought status in the Guanzhong plain of China [33]. However, gradient-based structural similarity (GSSIM) can better evaluate the quality of the image [34,35]. The structural information of the remote sensing image can reflect the spatial structure in the corresponding scene and the remote sensing images in the same area have the same or similar spatial structure. SSIM and GSSIM are used to measure the similarity of the structure information and structural characteristics of images [33,34,36].

In this study, a time-domain and pixel by pixel bi-parabolic NDVI-T_s space in which NDVI and T_s data from different time periods were mixed was re-examined in the Shaanxi province, China, which was vegetated from a low biomass to a high biomass. It was compared with the triangular NDVI-T_s space and one of the well-known drought index called the temperature-vegetation index (TVX) [37,38]. TVDI_c from the bi-parabolic NDVI-T_s space were applied to monitor distribution characteristics of drought between 2000 to 2016. Then, the GSSIM value of TVDI_c was calculated in 2000 to 2005, 2005 to 2010, 2010 to 2015, and 2000 to 2015, respectively, and analyzed the structural characteristics of TVDI_c in two phases. The spatial and temporal distribution features of the drought in the Shaanxi Province was quantitatively and accurately described based on GSSIM and compared with the change trend of drought based on the linear regression analysis method [39,40]. Finally, we discussed the possible reasons for drought changes.

2. Study Area and Data

2.1. Study Area

The Shaanxi province (105°29′–111°15′E, 34°42′N–39°35′N) is located in Northwest China and covers an area of 205,800 km². The elevation is higher in North and South Shaanxi when compared to the central region. From the north to the south, Shaanxi province are classified into three parts including the plateau areas in the north (called Shaanbei area), plain areas in the central region (called Guanzhong plain), and mountainous areas in the south (called Shaannan area) (Figure 1a). Shaanbei area riches in coal resources with windy desert soil lack water and vegetation cover. The Guanzhong plain is an important grain production base while Qinling-Dabashan Mountains are located in the Shaannan area with high vegetation cover (Figure 1b). It extends in three climatic zones from the north to the south, the middle temperate zone in the Shaanbei area, the warm temperate zone in the Guanzhong plain, and the subtropical zone in the Shaannan area. The annual average temperature and precipitation are about 13 °C and 580 mm, respectively, which decrease from the south to the north in the Shaanxi province [41].

2.2. Data

To produce dryness maps of the Shaanxi province, two tiles of moderate resolution imaging spectroradiometer (MODIS) data (h26v05, h27v05) acquired from 2000 to 2016 were used. The vegetation and land surface temperature (T_s) products with 1000 m resolution (MOD13A2 and MOD11A2, respectively) were downloaded from https://search.earthdata.nasa.gov/. Reflectance of red and NIR bands at 250 m resolution in an eight-day gridded level-3 product MOD09Q1 and MOD11A2 and filed measured soil moisture in April 2013 were used to verify the bi-parabolic NDVI-T_s space method. The projection system of these products are sinusoidal and geometrical. A radio-metrical correction has been made [1].

First, we derived NDVI and T_s between 2000 to 2016 from a monthly composited product MYD13A2 and eight-day composited product MYD11A2, respectively, by using the MODIS Reprojection Tool (MRT) (https://lpdaac.usgs.gov/tools/modis_reprojection_tool). The reflectance of Red and NIR bands were obtained from MOD09Q1 and resampled into 1000 m resolution by the MRT tool. At the same time, the sinusoidal projection of MODIS was converted into albers equal area projection with WGS-84 datum. Second, the maximum value composite method (MVC) was applied to composite 12 scenes NDVI images and 46 scenes T_s images into one image in a year, respectively [42]. MVC can effectively reduce the influence of cloud, atmosphere, the sun height angle [43,44], and the study area was extracted using the vector boundary of the Shaanxi province. Third, taking an interval of 0.01, the java program language was used to obtain the maximum and minimum temperature corresponding to the same NDVI. Then dry and wet edges were acquired from the scatter plots of the NDVI-T_s space.

The filed measured soil moisture, the precipitation, annual average temperature, and temperature anomaly from 2000 to 2013 in the Shaanxi province (Figure 1a) were all obtained from the meteorological station (http://data.cma.cn/).

3. Methodology

3.1. The Technical Route

Figure 2 shows the flowchart of constructing bi-parabolic NDVI-T_s space in this study. We first derived yearly maximum NDVI and T_s images between 2000 to 2016 from MOD13A2 and MOD11A2, respectively (see in Section 2). We then constructed scatter plots of NDVI-T_s space and re-examined if NDVI-T_s space was shaped in bi-parabola in the Shaanxi province (see in Section 4.1.1). Then we constructed TVDI_c and compared it with TVX and TVDI_t (see Section 5.1). We used filed measured soil moisture in April 2013 to validate TVDI_c, TVX, and TVDI_t (see in Section 4.3). TVDI_c was used to evaluate the spatial and temporal evolution of drought from 2000 to 2016 in the study area (see Section 4.1.2). In Section 4.2 and Section 5.2, we quantitatively and accurately described drought characteristics based on GSSIM and compared it with the change trend of drought based on the linear regression analysis method. Finally, we analyzed the driving forces for drought variations (see Section 5.3).

3.2. TVDI in Bi-Parabolic NDVI-T_s Space

TVDI, which was acquired from the triangular NDVI-T_s space, was first developed by Sandholt et al. (2002) [23].

$$\mathrm{TVDI}=\frac{{\mathrm{T}}_{\mathrm{s}}-{\mathrm{T}}_{\mathrm{smin}}}{{\mathrm{T}}_{\mathrm{smax}}-{\mathrm{T}}_{\mathrm{smin}}},$$

(1)

where

$$\begin{array}{l}{\mathrm{T}}_{\mathrm{smax}}={\mathrm{a}}_1+{\mathrm{b}}_1\times \mathrm{NDVI}\\ {\mathrm{T}}_{\mathrm{smin}}={\mathrm{a}}_2+{\mathrm{b}}_2\times \mathrm{NDVI}\end{array},$$

(2)

where ${\mathrm{T}}_{\mathrm{s}}$ is the remotely sensed surface temperature, ${\mathrm{T}}_{\mathrm{smin}}$ formed the wet edge and represents the minimum surface temperature, and ${\mathrm{T}}_{\mathrm{smax}}$ formed the dry edge and represents the maximum surface temperature. ${\mathrm{a}}_1$ and ${\mathrm{b}}_1$ are the linear fitting coefficient of the dry edge. ${\mathrm{a}}_2$ and ${\mathrm{b}}_2$ are the linear fitting coefficient of the wet edge. A linear least squares estimation was carried out to estimate the coefficients for the dry and wet edges and these coefficients can be obtained from scatter plots of the NDVI-T_s space.

TVDI ranges from zero to one. The bigger TVDI indicates that T_s gets to the dry edge and the drought gets worse. Conversely, the smaller TVDI demonstrates that T_s gets to the wet edge and the drought eases.

However, Liu et al. (2015) found that the NDVI-T_s space was bi-parabolic and the dry or wet edges were not linear (Figure 3), which can be expressed by Equation [32].

$$\begin{array}{l}{\mathrm{T}}_{\mathrm{smax}}={\mathrm{a}}_1\times {\mathrm{NDVI}}^2+{\mathrm{b}}_1\times \mathrm{NDVI}+{\mathrm{c}}_1\\ {\mathrm{T}}_{\mathrm{smin}}={\mathrm{a}}_2\times {\mathrm{NDVI}}^2+{\mathrm{b}}_2\times \mathrm{NDVI}+{\mathrm{c}}_2\end{array},$$

(3)

where ${\mathrm{a}}_1$, ${\mathrm{b}}_1$, ${\mathrm{c}}_1$, ${\mathrm{a}}_2$, ${\mathrm{b}}_2$, and ${\mathrm{c}}_2$ define the dry and wet edges as the parabolic fit, which can be acquired from scatter plots of the NDVI-T_s space.

3.3. Gradient-Based Structural Similarity (GSSIM)

Based on the SSIM, Yang et al. (2006) proposed an improved image quality assessed method called gradient-based structural similarity (GSSIM) basing on the edge information, which is the most important image structure information [34]. GSSIM combines the information of luminance comparison l(x, y), contrast comparison c(x, y), and gradient-based structure comparison g(x, y) in the image, which is shown in the equation below.

$$\mathrm{GSSIM}(\mathrm{x},\mathrm{y})={\left[\mathrm{l}(\mathrm{x},\mathrm{y})\right]}^{\alpha }{\left[\mathrm{c}(\mathrm{x},\mathrm{y})\right]}^{\beta }{\left[\mathrm{g}(\mathrm{x},\mathrm{y})\right]}^{\gamma }$$

(4)

where

$$\mathrm{l}(\mathrm{x},\mathrm{y})=\frac{{2\mathsf{\mu }}_{\mathrm{x}}{\mathsf{\mu }}_{\mathrm{y}}+{\mathrm{c}}_1}{{\mathsf{\mu }}_{\mathrm{x}}^2{+\mathsf{\mu }}_{\mathrm{y}}^2{+\mathrm{c}}_1}$$

(5)

$$\mathrm{c}(\mathrm{x},\mathrm{y})=\frac{{2\mathsf{\sigma }}_{\mathrm{x}}{\mathsf{\sigma }}_{\mathrm{y}}{+\mathrm{c}}_2}{{\mathsf{\sigma }}_{\mathrm{x}}^2{+\mathsf{\sigma }}_{\mathrm{y}}^2{+\mathrm{c}}_2}$$

(6)

$$\mathrm{g}(\mathrm{x},\mathrm{y})=\frac{{2\text{∑}}_{\mathrm{j}}{\mathsf{\Sigma }}_{\mathrm{i}}{\mathrm{G}}_{\mathrm{x}}{(\mathrm{i},\mathrm{j})\mathrm{G}}_{\mathrm{y}}{(\mathrm{i},\mathrm{j})+\mathrm{c}}_3}{{\text{∑}}_{\mathrm{j}}{\mathsf{\Sigma }}_{\mathrm{i}}{\left[{\mathrm{G}}_{\mathrm{x}}(\mathrm{i},\mathrm{j})\right]}^2{+\text{∑}}_{\mathrm{j}}{\mathsf{\Sigma }}_{\mathrm{i}}{\left[{\mathrm{G}}_{\mathrm{y}}(\mathrm{i},\mathrm{j})\right]}^2{+\mathrm{c}}_3}$$

(7)

where µ_x and µ_y are the mean value of TVDI in two phases, respectively, which reflects the luminance comparison information. σ_x and σ_y are standard deviation of TVDI in two phases, respectively, which reflects the contrast comparison information. G_x(i, j) and G_y(i, j) represents the gradient value of the pixel in the row_i and column_j of TVDI imageries in two phases, respectively. c₁, c₂, and c₃ are small constants that prevent the denominator from being zero, respectively, and c₁ = (K₁L)², c₂ = (K₂L)² and c₃ = c₂/2 with K₁, K₂ ≤ 1 and L is the range of a gray level in the image. The parameter of α, β, and γ are greater than zero. In this paper, α = β = γ = 1 and c₁ = c₂ = 0.0001 and c₃ = 0.0005 [33]. The higher the GSSIM value is, the more similar the TVDI value in two periods will be. This indicates that the drought conditions are closer with TVDI in the former period and the soil moisture status will not change much more.

3.4. Linear Regression Analysis

The linear regression analysis method can acquire regression slopes of all pixels based on a time series of satellite data, which can describe the change trends of NDVI [39,40]. We replace the NDVI by TVDI and the regression slopes can be obtained by

$$\mathrm{slope}=\frac{\mathrm{n}\times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}\times {\mathrm{TVDI}}_{\mathrm{i}}-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}\times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{TVDI}}_{\mathrm{i}}}{\mathrm{n}\times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{i}}^2-{({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i})}^2}$$

(8)

where the slope is the regression slopes, TVDI_i represents the TVDI image in a year, the parameter i is the number of the year, and i = 1, 2, 3, 4, … 17.

3.5. Pearson Correlation Analysis

The Pearson correlation can assess the correlation between two parameters or images [45,46], which is expressed below.

$$\mathrm{r}=\frac{\mathrm{N}\times \sum {\mathrm{x}}_{\mathrm{i}}{\mathrm{y}}_{\mathrm{i}}-\sum {\mathrm{x}}_{\mathrm{i}}\sum {\mathrm{y}}_{\mathrm{i}}}{\sqrt{\mathrm{N}\times \sum {\mathrm{x}}_{\mathrm{i}}^2-{\left(\sum {\mathrm{x}}_{\mathrm{i}}\right)}^2}\times \sqrt{\mathrm{N}\times \sum {\mathrm{y}}_{\mathrm{i}}^2-{\left(\sum {\mathrm{y}}_{\mathrm{i}}\right)}^2}}$$

(9)

where r is the Pearson correlation coefficient of two parameters, which ranges from −1 to 1. Higher values correspond to a better coherence between the inter-compared data sets. x_i represents the TVDI image in a year and i is the same as in Equation (8). y_i represents the annual average precipitation, the annual average temperature, and temperature anomaly images, respectively, which can be acquired by the inverse distance weighted method combined in GIS software.

4. Results

4.1. Spatio-Temporal Dynamic Changes of Drought

4.1.1. Scatter Plot of NDVI-T_s Space

We obtained the scatter plots of NDVI-T_s space between 2000 to 2016 by the java program language and the dry and wet edges fitting equations were defined by Excel. Figure 4 shows that the dry edge is clearly parabolic and the wet edge is characterized as a downward parabola. It indicates that NDVI-T_s space is roughly bi-parabolic. The coefficient of determination of most dry edges is high (R² > 0.74) and R² of most wet-edges is above 0.50, which indicates that the second-order polynomial can describe the trend of dry and wet edges.

4.1.2. Temporal and Spatial Evolution of Drought

Acquiring the fitting equation of dry and wet edges from Figure 4, TVDI images between 2000 to 2016 were obtained from Equation (1) and were mapped in GIS software. With an interval of 0.2, drought was classified into five categories [32]: very wet, wet, no dry, dry, and very dry. Figure 5 and

Table 1. The proportion of drought types in the Shaanxi province between 2000 to 2016.

	Drought Types and Proportion (%)
	Very Wet	Wet	No Dry	Dry	Very Dry	Drought Affected Area
2000	1.31	12.01	43.32	35.75	7.19	42.94
2001	1.24	14.65	40.72	34.83	8.13	42.96
2002	0.66	13.16	56.52	24.27	5.05	29.32
2003	1.95	30.72	44.42	18.38	4.34	22.72
2004	1.15	13.65	53.03	26.78	5.07	31.84
2005	3.00	22.38	36.52	29.10	8.55	37.65
2006	1.66	30.65	35.94	25.83	5.55	31.38
2007	2.31	18.63	50.30	23.96	4.25	28.20
2008	2.20	29.94	47.53	15.62	4.29	19.91
2009	2.31	24.83	36.44	30.55	5.44	35.99
2010	3.35	22.80	36.02	33.08	4.58	37.66
2011	3.74	28.60	39.70	23.32	4.22	27.55
2012	2.42	28.90	37.55	27.37	3.37	30.75
2013	2.68	22.87	45.33	22.97	5.93	28.90
2014	1.95	28.69	52.42	14.66	1.85	16.50
2015	0.58	9.14	59.65	26.79	3.51	30.29
2016	1.28	22.22	48.56	24.25	3.28	27.54

show the spatial distribution of TVDI and the proportion of drought types in the Shaanxi province between 2000 to 2016, respectively.

In space, the drought in the Shaanxi province mainly occurred in the Shaanbei area, central and eastern regions of the Guanzhong plain, and eastern parts of the Shaannan area. There was no drought in most areas of Shaannan and southern and northern parts of the Guanzhong plain in the past seventeen years.

(i) The proportion of drought affected areas in the Shaanxi province fluctuated between 16.50% and 42.96% in the past 17 years and the mean value of proportion was 30.71%. The drought area was the smallest in 2014, which was 16.50% and the drought area was up to 42.96% in 2001.

(ii) Drought area accounted for 42.94% in 2000, which was mainly concentrated in the Shaanbei area, central and eastern parts of the Guanzhong plain, and some parts in the eastern region of Shaannan area. The drought situation in 2001 was more severe than of the one in 2000. It accounted for 42.96%, which mainly occurred in the Shaanbei area and the Guanzhong plain. In the following three years, the drought conditions were significantly released in the Shaanbei area, but droughts were distributed in the Yulin city and central parts of the Guanzhong plain, which accounted for 29.31%, 22.71%, and 31.84%, respectively. It aggravated in 2005, which was 37.65%. The drought was eased year by year, which accounted for 31.38%, 28.20%, and 19.91%, respectively, between 2006 to 2008. The drought was aggravated in the Shaanbei area the following two years, which were 35.99% and 37.66%, respectively. The drought condition eased in the Shaanbei area, which was more severe than in the Guanzhong plain from 2011 to 2013. Drought eased in the whole province, which accounted for 16.50% in 2014 and was the smallest drought in the past 17 years. In 2015, the wet area decreased in Yan’an city and the Shaannan area and the area of the drought enlarged in Shaanbei, which was 30.29%. The wet area was increased and the drought conditions were slightly eased in the Shaanbei area and the central parts of the Guanzhong plain in 2016, which accounted for 27.54%.

4.2. Quantitative Analysis of Spatial Distribution of Drought

Based on Equation (4), GSSIM imageries were obtained by Matlab software between two years including 2000–2005, 2005–2011, 2011–2016 and 2000–2016 and mapped in GIS software. GSSIM was classified into three categories: (i) mutation (0 < GSSIM ≤ 0.25) [33], which indicated that the value of TVDI in two years has mutated and the drought status has changed significantly; (ii) Moderate change (0.25 < GSSIM ≤ 0.65), which showed that the value of TVDI in two years has been altered and drought conditions had changed moderately, (iii) Low change (0.65 < GSSIM ≤ 1), which indiated that value of TVDI in two years are similar and the drought status is almost unchanged.

Figure 6 shows that the distribution of the mutated area is similar in four periods. It is mainly located in the north and east of Yulin city including both sides of the boundary between Yulin and Yan’an cities, the middle and southern parts of Yan’an city, the northern and central parts of Weinan city, both sides of the boundary between Hanzhong and Baoji cities, the southeast of Shangluo city. The drought in these areas changes significantly. Figure 6 also indicates that the moderately changed area surrounds the mutated area and the drought status changes moderately. However, the slightly changed area surrounds the moderately changed area and the drought status is almost unchanged.

Taking GSSIM between 2008 and 2015 as an example, samples of mutation A, B, C, and D and samples of moderate change E, F, and G were selected to analyze the variation of drought, respectively (Figure 7). High-resolution images obtained from Google Earth™ in 2008 and 2015 combined with the difference in the TVDI image between 2008 and 2015 were used to verify drought changes. Figure 8 demonstrates that the river changes into a barren area in mutated sample A and vegetation covered sample B becomes the rural area. Figure 8 shows that farmland in the mutated sample C moves into ditches and the vegetation coverage obviously increases in the mutated sample D. Combined with the differences in the TVDI image between 2008 and 2015, the drought conditions in mutated samples A and B were aggravated while it was significantly relieved in mutated samples C and D. Figure 8 shows that the coal yard enlarges from 2008 to 2015 in moderate change sample E and some vegetation cover parts in moderate change sample F become the open coal mining area. A region of the open coal mining area is covered with vegetation in the moderate change sample G (Figure 8). It indicated that the underlies were moderately changed and the drought status was moderately aggravated in moderate change samples E and F while it was moderately eased in moderate change sample G.

To analyze the changes of GSSIM in different periods quantitatively, we obtained the area and proportion of drought variation.

Table 2. The proportion of GSSIM in the Shaanxi province between 2000–2016 (%).

	Times	Drought Status	2000–2005	2005–2011	2011–2016	2008–2015	2000–2016
GSSIM
GSSIM ≤ 0.25 (mutation)		eased	3.20	5.11	2.68	1.34	3.53
		aggravated	2.22	2.87	3.04	3.44	1.81
0.25~0.65 (moderate change)		eased	25.58	27.00	15.77	7.86	28.45
		aggravated	14.27	13.29
GSSIM > 0.65 (low change)		changed a little	54.73	51.73	54.37	54.69	54.44

describes that the area of drought eased in the mutated area and moderate change areas increased from 3.20%, 25.58% during 2000–2005 to 5.11%, 27.00% during 2005–2011 and then decreases to 2.68% and 15.77% during 2011–2016, respectively. We can see that the area of drought aggravated in the moderate change region slightly decreases from 14.27% between 2000 to 2005 to 13.27% between 2005 to 2011 and increases to 24.15% between 2011 to 2016. Table 2 indicates that the area of drought aggravated in the mutated region slightly increases in these three phases and more than half of Shaanxi province’s drought conditions has no obvious change in five phases. Table 2 also shows that the area of drought relived is larger than that of aggravated between 2000–2005, 2005–2011, and 2000–2016 and the drought are relieved while the area of drought affected is enlarged between 2008–2015 and 2011–2016, which coincide with the change of the mean value of TVDI in these years (

Table 3. The mean value of TVDI in the Shaanxi province.

2000	2005	2008	2011	2015	2016
0.5708	0.5372	0.4793	0.4910	0.5577	0.5117

).

4.3. Relations between Filed Measured Soil Moisture, TVDI, and TVX

Soil moisture (SM) can be used as a direct indicator for the drought status. We applied filed measured soil moisture to validate the performance of TVDI_c compared with TVX and TVDI_t. Soil moisture at 10 cm depth was available from 26 agricultural meteorology-observing stations across the whole region for the 8, 18, and 28 days of each month. However, it is often difficult to get cloud-free images on these days and the vegetation over semi-arid regions, which needs roughly 5 days to respond to a soil moisture change [47]. The reflectance of red and NIR band in 8-day gridded level from MOD09Q1 and T_s from MOD11A2 in April 2013 were applied to acquire TVDI_c, TVDI_t, and TVX, respectively. Figure 9 indicates that there are significant negative correlations (reliable at 1% significance level) between TVDI_c, TVDI_t, and 10-cm-depth soil moisture during three periods in April 2013, respectively. The coefficient of determination R² = 0.37~0.43 between TVDI_c and 10-cm-depth soil moisture is obviously higher than TVDI_t. In contrast, the TVX is also negatively related to soil moisture with R² = 0.005~0.0625, which is obviously lower than TVDI_c and TVDI_t (Figure 9c,f,i). Although TVX is very simple and integrates both the reflective and thermal information of remote sensing data, it is easily influenced by the variation of land cover [48], atmospheric impacts, cloud and sensor drift, etc. [49].

5. Discussion

5.1. Comparison between the TVDI_c, TVDI_t and TVX

5.1.1. Scatter Plots of Bi-Parabolic and Triangular NDVI-T_s Space

It shows that the NDVI-T_s space is bi-parabolic during three periods in April 2013 (Figure 10a,c,e), the coefficient of determination R² > 0.80 in the dry edge, and R² > 0.38 in the wet edge. In contrast, NDVI-T_s space becomes triangular during three periods in April 2013 when parts of NDVI < 0.15 are omitted from the bi-parabolic NDVI-T_s space (Figure 10b,d,f), with R² > 0.85 in the dry edge and R² > 0.21 in the wet edge. Although R² in the dry edge in triangular NDVI-T_s space is higher than that of bi-parabolic NDVI-T_s space, R² between TVDI_c and 10 cm soil moisture during these three periods is obviously higher than that of TVDI_t and TVX (Figure 9). The results show that the bi-parabolic NDVI-T_s space is better than the triangular NDVI-T_s space for monitoring drought conditions and TVDI_c is better than TVDI_t and TVX for assessing 10-cm-depth soil moisture. The NDVI-T_s space should include the area with NDVI < 0.15, which should not be omitted [32]. Figure 3 shows that the ability of evaporation is strengthened from water to the bare surface and the partial vegetation area. Therefore, there is a positive linear relationship between NDVI and T_s when NDVI < 0.15 and a negative linear relationship between NDVI and T_s when NDVI > 0.15 [32].

The main hypothesis of T_s/NDVI approaches is that atmospheric conditions (including elevation effects) are relatively constant or uniform over the whole image [50,51]. In the triangular NDVI-T_s space, the dry and wet edges are derived from scatter plots of NDVI and T_s images in the same time collected in a given spatial domain [51]. TVDI_c is based on a pixel by pixel bi-parabolic NDVI-T_s space, which is obtained by 8-day (yearly) composite NDVI and T_s images. The proposed bi-parabolic NDVI-T_s space mixes time points together from periods with potentially different atmospheric conditions. In this space, the “time” is the data domain and a single pixel is the spatial domain.

As previously discussed, we can see that combining data from different time periods together with potentially different atmospheric conditions goes against the theoretical basis of the T_s/NDVI approach, but that doing so appears to be empirically valid based on the results in Figure 10. The potential reasons may be that the air temperature was relatively constant in theses 8-day periods when atmospheric conditions are likely to be relatively stable. We used MVC to composite 46 scenes 8-day T_s images into one image in a year. The value of every pixel in yearly composite NDVI and T_s images is the maximum value of these 12 scenes monthly NDVI and 46 scenes 8-day T_s images, respectively, and the maximum NDVI and air temperature in the Shaanxi province often occurs in July/August in a particular year. In this case, the atmospheric conditions are likely to be also relatively stable. In the research of Liu et al. (2017), Cao et al. (2017), and Khan et al. (2018), they also used MVC to composite 8-day T_s and 16-day NDVI images into monthly or yearly images and it obtained dry and wet edges from time series MDOIS data [2,52,53].

5.1.2. Dryness Map of TVDI_c, TVDI_t, and TVX

Figure 11 shows drought conditions during three periods in April 2013 from TVDI_c, TVDI_t, and TVX, respectively. Drought status assessed by TVDI_c was very similar to that of TVDI_t. However, there were some differences in the Shaanbei area where drought conditions were more severe by TVDI_c than TVDI_t. Drought conditions monitored by TVX were very different from TVDI_c and TVDI_t, which was more severe in the Shaanbei area and the soil water content was higher in the Shaannan area. As discussed above, TVX is often influenced by other factors, which may cause more disturbances to the relationship between the soil water condition and TVX. Figure 11 also indicates that the drought was aggravated in the Yulin city and most areas of Shaannan from 7 April, 2013 to 14 April, 2013 to 15 April, 2013 to 22 April, 2013, which was eased in most areas of Shaannan between 23 April, 2013 to 30 April, 2013.

5.2. Comparison between the GSSIM and Linear Regression Analysis

The regression slopes of all pixel locations for TVDI changes trends during 2000–2016 (reliable at 5% significance level), which were derived from Equation (8) and, therefore, categorized into “drought eased not significantly,” “drought eased significantly,” “drought aggravated significantly,” and “drought aggravated not significantly” and then mapped accordingly. Figure 12b demonstrates that drought conditions accounting for 69.31% are not significantly alleviated, which is mainly distributed in most areas of Shaannan, southern and northern parts of Guanzhong plain, northern and southern area of Yulin city, and southern, northern, and central parts of Yan’an city. In other words, the drought conditions in these areas have not changed much more between 2000 to 2016. At the same time, drought conditions have been significantly alleviated, which accounted for 14.45% and were mainly located in northern parts of Yulin city, central parts of Yan’an city, central area of Guanzhong plain, and some parts in Shaannan. Figure 12a and Table 2 show that the area of drought eased in the mutation area accounts for 3.53%, which is lower than that of linear regression slope and the distribution is similar to the regression slope. The drought is not significantly aggravated in southern and northwestern parts of the Yulin city, the northern and southern part of the Guanzhong plain, and the southern parts of Shannan, which accounts for 15.35%. The drought, significantly aggravated in theses area, only accounts for 0.88% while the drought conditions aggravated area accounts for 1.81% by GSSIM. There are some apparent differences between GSSIM and the linear regression slope in northern regions of Shaanbei area and some parts of Shaannan area. It should be noted that only two TVDI images were used to obtain GSSIM while 17 scenes of TVDI were used to acquire the linear regression slope. This may lead to some uncertainties in the results.

5.3. The Reasons for Drought Status Variations

The Pearson correlation coefficients of all pixel locations for assessing the relationship between TVDI and meteorological factors from 2000 to 2013 (reliable at 10% significance level) were derived from Equation (9) and, therefore, categorized into “negative significant,” “negative not significant,” “positive significant,” and “positive not significant” and then mapping accordingly (Figure 13). Figure 13a shows that there is a significant negative correlation between TVDI and annual average precipitation in most parts of Yulin city, central parts of Yan’an city, north of Ankang and Weinan cities, eastern area of Shangluo city, western and northern parts of Baoji area, which accounted for 23.74%. It is indicated that the drought in these areas is mainly affected by annual rainfall. When the rainfall increases, TVDI decreases and the drought can be relieved. The correlation between TVDI and annual average precipitation in the rest area of the Shaanxi province has not passed the significance test, which indicates that annual precipitation is not a dominant factor for drought variations in these regions.

It describes in Figure 13b that a significant positive relationship exists between the annual average temperature and TVDI in central and southern parts of Baoji city and west and northwest areas of Hanzhong city, which accounts for 8.21%. It is expressed that drought conditions in these areas are mainly influenced by the annual average temperature. As the temperature increases, TVDI increases and the drought is aggravated. There is a significant negative relationship between the annual average temperature and TVDI in central and southern parts of Yan’an city, north of Weinan city, most area of Shangluo city and east of Xi’an city, which accounts for 14.11%. It illustrates that drought in these regions is mainly influenced by the annual average temperature. When the temperature increases, TVDI decreases and the drought eases. Figure 13c demonstrates that there is a positive correlation between TVDI and temperature anomaly, which accounts for 58.70% while 35.62% of the whole province are negatively related with a temperature anomaly. However, the Person correlation coefficient is not reliable at a 10% significance level and indicates that the temperature is not the main factor for the change of drought status in the Shaanxi province.

As discussed in Section 4.2, the changes of drought were not only influenced by precipitation and temperature, but the change of underlies, the vegetation improvement, and the enlargement of open mining and buildings and irrigation.

6. Conclusions

MODIS NDVI and T_s time series data were used to reexamine the time-domain bi-parabolic NDVI-T_s space when the area was vegetated from low biomass to high biomass and compared with the triangular NDVI-T_s space and TVX. Then the GSSIM was used to quantitatively evaluate the variation of drought in two periods. Dry and wet edges in this research exhibited a parabolic shape again in a scatter plot of T_s and NDVI, which was linear in the triangular NDVI-T_s space. TVDI_c, TVDI_t, and TVX were compared with 10-cm depth relative soil moisture. It was shown that TVDI_c was better than TVDI_t, which were all apparently better than TVX in revealing 10-cm depth soil moisture. The bi-parabolic NDVI-T_s space is the complement and enlargement of the triangular NDVI-T_s space. In the past 17 years, there was a slightly fluctuating trend of the drought status and the results of GSSIM indicated that the drought in more than half of the Shaanxi province changed a little and was relatively stable. The change of meteorological factors and underlies are the reasons for drought variations. We will do more work to examine whether the NDVI-T_s space with “time” and a single pixel spatial domain is bi-parabolic with other study areas and verify if the atmospheric conditions remain constant with yearly composite NDVI-T_s space in the future.