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

Traditional NDVI-Ts space is triangular or trapezoidal, but Liu et al. (2015) discovered that the NDVI-Ts 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 Ts data to reexamine if the NDVI-Ts 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-Ts 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 Ts and NDVI in the Shaanxi province, which was linear in the triangular NDVI-Ts space. The Temperature Vegetation Dryness Index (TVDIc) was obtained from bi-parabolic NDVI-Ts andTVDIt was obtained from the triangular NDVI-Ts space and TVX were compared with 10-cm depth relative soil moisture. By estimating the 10-cm depth soil moisture, TVDIc was better than TVDIt, which were all apparently better than TVX. Second, combined with MODIS data, the drought conditions of the study area were assessed by TVDIc 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.


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

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 2 . 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].

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 (Ts) 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-Ts space method. The projection system of these products are sinusoidal and geometrical. A radiometrical correction has been made [1].
First, we derived NDVI and Ts 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 Ts 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

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/). 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 Sections 4.2 and 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).

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]. where where T s is the remotely sensed surface temperature, T smin formed the wet edge and represents the minimum surface temperature, and T smax formed the dry edge and represents the maximum surface temperature. a 1 and b 1 are the linear fitting coefficient of the dry edge. a 2 and 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].
where a 1 , b 1 , c 1 , a 2 , b 2 , and 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.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 22 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-Ts space. TVDI ranges from zero to one. The bigger TVDI indicates that Ts gets to the dry edge and the drought gets worse. Conversely, the smaller TVDI demonstrates that Ts gets to the wet edge and the drought eases.
However, Liu et al. (2015) found that the NDVI-Ts space was bi-parabolic and the dry or wet edges were not linear ( Figure 3), which can be expressed by Equation [32]. T =a NDVI +b NDVI+c where a , b , c , a , b , and c define the dry and wet edges as the parabolic fit, which can be acquired from scatter plots of the NDVI-Ts space.

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.

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.
where l(x, y) = 2µ x µ y + c 1 µ 2 x +µ 2 y +c 1 (5) c(x, y) = 2σ x σ y +c 2 σ 2 x +σ 2 y +c 2 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 1 , c 2 , and c 3 are small constants that prevent the denominator from being zero, respectively, and c 1 = (K 1 L) 2 , c 2 = (K 2 L) 2 and c 3 = c 2 /2 with K 1 , K 2 ≤ 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 1 = c 2 = 0.0001 and c 3 = 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.

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 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.

Pearson Correlation Analysis
The Pearson correlation can assess the correlation between two parameters or images [45,46], which is expressed below.
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.

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 2 > 0.74) and R 2 of most wet-edges is above 0.50, which indicates that the second-order polynomial can describe the trend of dry and wet edges.

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 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

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 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%.  [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.

Quantitative Analysis of Spatial Distribution of Drought
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 22 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.    To analyze the changes of GSSIM in different periods quantitatively, we obtained the area and proportion of drought variation.  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 (Table 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 2 = 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 2 = 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].

Scatter Plots of Bi-Parabolic and Triangular NDVI-Ts Space
It shows that the NDVI-Ts space is bi-parabolic during three periods in April 2013 (Figure 10a,c,e), the coefficient of determination R 2 > 0.80 in the dry edge, and R 2 > 0.38 in the wet edge. In contrast, NDVI-Ts space becomes triangular during three periods in April 2013 when parts of NDVI < 0.15 are omitted from the bi-parabolic NDVI-Ts space (Figure 10b,d,f), with R 2 > 0.85 in the dry edge and R 2 > 0.21 in the wet edge. Although R 2 in the dry edge in triangular NDVI-Ts space is higher than that of bi-parabolic NDVI-Ts space, R 2 between TVDIc and 10 cm soil moisture during these three periods is obviously higher than that of TVDIt and TVX (Figure 9). The results show that the bi-parabolic NDVI-Ts space is better than the triangular NDVI-Ts space for monitoring drought conditions and TVDIc is better than TVDIt and TVX for assessing 10-cm-depth soil moisture. The NDVI-Ts 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 Ts when NDVI < 0.15 and a negative linear relationship between NDVI and Ts when NDVI > 0.15 [32].

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 2 > 0.80 in the dry edge, and R 2 > 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 2 > 0.85 in the dry edge and R 2 > 0.21 in the wet edge. Although R 2 in the dry edge in triangular NDVI-T s space is higher than that of bi-parabolic NDVI-T s space, R 2 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].
Ts images into one image in a year. The value of every pixel in yearly composite NDVI and Ts images is the maximum value of these 12 scenes monthly NDVI and 46 scenes 8-day Ts 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 Ts 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]. 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  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 Figure 11 shows drought conditions during three periods in April 2013 from TVDIc, TVDIt, and TVX, respectively. Drought status assessed by TVDIc was very similar to that of TVDIt. However, there were some differences in the Shaanbei area where drought conditions were more severe by TVDIc than TVDIt. Drought conditions monitored by TVX were very different from TVDIc and TVDIt, 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

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.

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.

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. (b) Significant linear regression slope values for trends derived from TVDI observations during 2000-2016 (reliable at 5% significance level).

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.

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.