Modeling NDVI Using Joint Entropy Method Considering Hydro-Meteorological Driving Factors in the Middle Reaches of Hei River Basin

Abstract

Terrestrial vegetation dynamics are closely influenced by both hydrological process and climate change. This study investigated the relationships between vegetation pattern and hydro-meteorological elements. The joint entropy method was employed to evaluate the dependence between the normalized difference vegetation index (NDVI) and coupled variables in the middle reaches of the Hei River basin. Based on the spatial distribution of mutual information, the whole study area was divided into five sub-regions. In each sub-region, nested statistical models were applied to model the NDVI on the grid and regional scales, respectively. Results showed that the annual average NDVI increased at a rate of 0.005/a over the past 11 years. In the desert regions, the NDVI increased significantly with an increase in precipitation and temperature, and a high accuracy of retrieving NDVI model was obtained by coupling precipitation and temperature, especially in sub-region I. In the oasis regions, groundwater was also an important factor driving vegetation growth, and the rise of the groundwater level contributed to the growth of vegetation. However, the relationship was weaker in artificial oasis regions (sub-region III and sub-region V) due to the influence of human activities such as irrigation. The overall correlation coefficient between the observed NDVI and modeled NDVI was observed to be 0.97. The outcomes of this study are suitable for ecosystem monitoring, especially in the realm of climate change. Further studies are necessary and should consider more factors, such as runoff and irrigation.

1. Introduction

Terrestrial vegetation plays a key role in energy, water, and biogeochemical cycles, while variation in vegetation can significantly influence atmospheric and hydrological processes. Vegetation is a sensitive indicator of global change and can also be a natural link between the atmosphere, land surface, soil, and water [1,2]. Changes in vegetation cover are influenced by climate change, human activities, and the atmospheric CO₂ fertilization effect [3,4,5]. It is necessary to explore the response mechanism of vegetation dynamics to hydrological processes and climate change.

Vegetation emergence and senescence are closely related to characteristics of the lower atmosphere, including the annual cycle of weather pattern shifts, temperature and precipitation [6]. In particular, the amount and duration of precipitation and temperature play a significant role in controlling vegetation development [2,7,8,9]. During the last decades, the global average land surface air temperature has been increasing systematically in the northern middle and high latitude [10]. Elevated air temperature increased plant growth by lengthening the growing season [11], enhancing photosynthesis, and altering nitrogen availability by accelerating decomposition or mineralization. Precipitation is the dominant controlling climatic factor in water-limited, semi-arid and arid regions such as the northern China [10]. The increase in precipitation contributes to the growth of vegetation. Several studies have reported relationships between vegetation indices and precipitation and temperature [3,12,13,14,15,16,17]. However, most of these studies focused solely on either precipitation or temperature for vegetation models, thereby resulting in low precision and low correlation coefficients between vegetation indices and climatic variables [17]. For instance, Ichii et al. (2002) found a weaker relationship between vegetation growth and precipitation [18], while Nicholson et al. (1990) reported a weaker relationship between NDVI and precipitation, indicating a nonlinear overall relationship [19]. However, apart from the influence of climatic change, vegetation is also strongly related to the groundwater table level (GTL) in arid and semi-arid regions [2,20,21,22]. Groundwater is a key impact factor in sustaining the ecological environment, by supplying soil moisture through capillarity for maintaining the root system. The shallower the groundwater depth the more available soil moisture, and vice versa [22]. Many investigators have also reported the relationship between vegetation and groundwater availability [2,21]. For instance, Jin et al. (2016) found a linear correlation between NDVI and groundwater depth in the Qaidam basin [2]. The distribution of GTL is becoming more heterogeneous in the middle reaches of the Hei River as a result of increasing cultivated land area and runoff from Yingluoxia station [21]. The distribution of vegetation is also highly heterogeneous and depends upon multiple factors. In this study, different statistical tools, as well as the joint entropy method, were employed in order to understand the relationships between NDVI and precipitation, temperature, and GTL. This will be necessary for investigating the relevance of coupling variables.

The joint entropy expresses the total amount of uncertainty contained in the union of events [23] either in a linear or nonlinear system [24] and has been employed for various purposes. For instance, Li et al. (2012) employed joint entropy to design appropriate hydrometric networks [25], while Mishra et al. (2012) employed the entropy theory to overcome the nonlinear relationship between precipitation and vegetation, and further identified the downscaling method yielding higher mutual information [23]. Entropy theory has also been applied in ecological studies [5,17,26,27]. For instance, Sohoulande, Djebou and Singh (2015) employed joint entropy to retrieve vegetation growth patterns from climatic variables. However, this study emphasizes the mutual information of NDVI and the influencing factors.

In this paper, the mutual information of NDVI was investigated by coupling cumulative precipitation, the average temperature of growing season, and GTL. To that end, the study area was divided into 5 sub-regions, based on 3 spatial distributions of mutual information. In each sub-region, NDVI was modeled using nested statistical models at the grid and regional scales, respectively.

2. Methodology

2.1. Trend Analysis

Trends for each pixel were represented to reflect the characteristics of vegetation cover over different periods of time [28]. The rate of change of NDVI was calculated as [29].

$${\theta }_{slope}=\frac{n\times {\displaystyle \sum _{i=1}^{n}i\times NDV{I}_i}-{\displaystyle \sum _{i=1}^{n}i{\displaystyle \sum _{i=1}^{n}NDV{I}_i}}}{n\times {\displaystyle \sum _{i=1}^n{i}^2-({\displaystyle \sum _{i=1}^{n}i{)}^2}}}$$

(1)

where ${\theta }_{slope}$ is the slope of the trend line of annual NDVI; i represents the years; n is the time span; and NDVI_i is the NDVI value for the i-th year.

${\theta }_{slope}$ describes the trend in the annual NDVI within the study area. If ${\theta }_{slope}>0$, the NDVI increases; or else, the NDVI decreases (${\theta }_{slope}<0$) or remains constant (${\theta }_{slope}=0$).

2.2. Joint Entropy and Mutual Information

Mutual information was employed to explore the relationship between NDVI and influencing factors. Considering a random variable X, for each value of X, x_i represents an event with a corresponding probability of occurrence, p_i. The entropy H(X) can be expressed as:

$$H(X)=-{\displaystyle \sum _{i=1}^n{p}_i}{\log }_2(p_i)$$

(2)

Where i means the ith event. The logarithm is based on 2, because it is more convenient to use than logarithms based one or 10, however, the base can be taken as other than 2 without alteration [24]. The entropy is thereby measured in bits. The probability density function (PDF) of variable X is obtained using discrete intervals for the values of X. The discrete PDF is defined for n equal-width bins defined for the range of X [17,24,30]. The entropy H(X) is a measurement of information or uncertainty [24]. Similarly, the joint entropy can be computed for a joint probability distribution of two or more variables [24,25]. Specifically, with three variables, the joint PDF is computed based on a three-dimensional contingency analysis as illustrated in

Table 1. Illustration of a three-dimension contingency employed to compute the joint PDF. Consider 3 variables X, Y, and Z.

Variables			Discrete Probabilities
X = {X1,X2}	Y = {Y1,Y2}	Z = {Z1,Z2}
X = X1	Y = Y1	Z = Z1	p(X = X1, Y = Y1, Z = Z1)
		Z = Z2	p(X = X1, Y = Y1, Z = Z2)
	Y = Y2	Z = Z1	p(X = X1, Y = Y2, Z = Z1)
		Z = Z2	p(X = X1, Y = Y2, Z = Z2)
X = X2	Y = Y1	Z = Z1	p(X = X2, Y = Y1, Z = Z1)
		Z = Z2	p(X = X2, Y = Y1, Z = Z2)
	Y = Y2	Z = Z1	p(X = X2, Y = Y2, Z = Z1)
		Z = Z2	p(X = X2, Y = Y2, Z = Z2)

Note: Each variable is categorized into 2 classes (n = 2).

. Considering the three variables precipitation, temperature and NDVI, the process generates a discrete PDF made of n³ discrete probabilities such that:

$${\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{k=1}^{n}p(P=P_i}}},T=T_j,NDVI=NDV{I}_k)=1$$

(3)

The joint entropy of the three variables can be defined as:

$$H(P,T,NDVI)=-{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{k=1}^{n}p(P_i,T_j,NDV{I}_k){\log }_2[p(P_i,T_j,NDV{I}_k)]}}}$$

(4)

where n means the number of events.

The mutual information represents the amount of information common to both X and Y and provides a general measure of dependence between the random variables. It is superior to the Pearson correlation coefficient, because it captures both linear and nonlinear dependence while the Pearson correlation coefficient is only suitable for linear relationships [25]. It equals the difference between the sum of two marginal entropies and the total entropy:

$$M(X,Y)=H(X)+H(Y)-H(X,Y)$$

(5)

The two-dimensional cases of the amount of information transmitted can be extended to three or more dimensional cases [31]. Considering the three variables of precipitation, temperature, and NDVI, precipitation and temperature are considered inputs while NDVI is considered an output. Therefore:

$$M(P,T;NDVI)=H(P,T)+H(NDVI)-H(P,T,NDVI)$$

(6)

Moreover, applying the same methodology, the coupling of P and GTL; T and GTL along with NDVI were further explored. The ensemble E of candidate equations can be summarized as:

$$E=\left\{NDV{I}_{GTL,P}=f(GTL,P);NDV{I}_{GTL,T}=f(GTL,T);NDV{I}_{T,P}=f(T,P)\right\}$$

(7)

3. Study Area and Data

The study area occupied the middle reaches of the Hei River basin (38°30′–39°55′N, 98°55′–100°55′E; Figure 1), which is located in the middle of Hexi corridor of Gansu Province, northwest of China, and belongs to an arid climate zone. The total area is about 10,000 km². The mean annual temperature is around 8 °C, the mean annual precipitation is between 69 and 216 mm and is concentrated between June and September, while the mean annual potential evaporation is between 1453 and 2351 mm [32].

The Hei River is the second largest inland river in China with a total length of 821 km and a runoff of 16.0 × 10⁸ m³·year⁻¹. It flows through Yingluoxia and Zhengyixia stations, entering and exiting middle reaches. Runoff has shown an increasing/decreasing trend in the past decade. The annual runoff at Yingluoxai station increased to 16.0 × 10⁸ m³ in the 1990s from around 14.4 × 10⁸ m³ in the 1960s (Figure 2), while the runoff at Zhengyixai station decreased by 3 × 10⁸ m³ during the same time period (Figure 2). This situation accelerated desertification within the northern parts of the basin [33,34]. In order to mitigate deteriorating ecosystems in the Hei River basin, an Ecological Water Diversion Project (EWDP) was initiated by the Chinese government in 2000 to ensure the delivery of water to the lower reaches for ecological water needs [35]. After that, runoff at Zhengyixia station has increased to levels not recorded since the 1960s [33,36]. The implementation of the EWDP alleviated the ecological deterioration in the lower reaches of the Hei River but aggravated the water resources shortage in the middle reaches. In order to meet the increasing demand for water resources, groundwater exploration provided an option. After 2004, the recharge of groundwater also increased due to the increase in precipitation, irrigation and runoff replenishment. The temporal variation in groundwater levels of the observation wells over the period from2001 to 2011 is shown in Figure 3c which shows a rising trend in groundwater level along the Hei River [37]. However, the water level in the region away from the river showed a declining trend.

The annual NDVI of each pixel was calculated by the annual maximum value of NDVI in order to eliminate the noise from the cloud and solar altitude. The Moderate-resolution imaging spectroradiometer (MODIS)-based NDVI products(MYD13A2 and MOD13A2)used for this study were acquired from the National Natural Science Foundation of China “Environmental and Ecological Science Data Center in the West of China (http://westdc.westgis.ac.cn/)” with a spatial resolution of 1 km and a temporal resolution of 1 day, groundwater depths of 53 wells were obtained from the National Natural Science Foundation of China “Hei River Project Data Management Center (http://westdc.westgis.ac.cn/)” while monthly precipitation and temperature datasets were acquired from the China Meteorological Network (http://data.cma.cn/), with a 0.5° spatial resolution, interpolated according to the observed data of 2747 stations in China.

The ordinary kriging method was employed to interpolate the observed values of ground water depths. For consistency regarding the NDVI spatial resolution, the interpolated groundwater table levels, precipitation, and temperature datasets were rescaled to a 1 km spatial resolution.

4. Results and Discussion

4.1. Variation of Hydrometeorological Variables and NDVI

During the period from2001 to 2011, both the average temperature and cumulative precipitation increased at the three hydrologic stations (Figure 3a,b). The temperature increases significantly at Gaotai station (p < 0.05). However, the trend of temperature increase is not obvious at another two stations. The distributions of cumulative precipitation and average temperature are shown in Figure 4a,b. Temperature gradually rises from south to north, while precipitation decreases were also observed. This is mainly determined by topography.

Figure 5a illustrates the annual maximum NDVI in the middle reaches of the Hei River basin from 2000 to 2011. The slope of the linear regression trend-line was 0.005. The NDVI time series corresponded to the average of all the pixels in the study area and reflected the overall trend. Figure 5b shows the spatial distribution of the linear regression slope (${\theta }_{slope}$). In most cases, the area in which NDVI increased occupied about 94% of the total area. It can be seen that NDVI changes obviously. In order to explore the influencing factors on NDVI variation, the whole region was divided into 5 sub-regions based on the joint entropy method, and different regression models were constructed considering hydro-meteorological variables.

4.2. Mutual Information of NDVI with Coupling of Variables

To understand the relationship between NDVI and hydro-meteorological variables, mutual information was employed to measure the dependence between the NDVI time series and each coupled variable. The greater the values of mutual information the higher the correlation. The spatial distribution of mutual information between NDVI and coupled variables are shown in Figure 6a–c, respectively. It can be observed from the figures that the mutual information of NDVI and coupled variables of precipitation and GTL generally decreased from southeast to northwest. The mutual information values of NDVI and the coupled variables of precipitation and temperature were dispersedly distributed, increasing from the center to the edge of regions, while the mutual information of NDVI and the coupled variables of temperature and groundwater table levels were greater in the northwestern and southeastern regions. In order to appropriately model the NDVI using dependent variables, the study region was artificially divided into 5 sub-regions (Figure 6d), based on three spatial distributions of mutual information, topographic and land use maps. In order to simply construct models, as far as possible each sub-region was made to have the same characteristics. For instance, sub-region I and sub-region IV are mainly mountainous and M (P, T; NDVI) is higher than M (T, GTL; NDVI) and M (P, GTL; NDVI) in most areas of these sub-regions. As the same time, the land use types are desert. Likewise, in most areas of sub-region II and sub-region V, M (GTL, P; NDVI) is higher than M (P, T; NDVI) and M (T, GTL; NDVI), and these areas are distributed as artificial oases.

4.3. Modelling NDVI

In sub-regions I and IV, precipitation and temperature were used to model the NDVI, since the mutual information of NDVI and the variables of precipitation and temperature were relatively higher. GTL and precipitation were used to model the NDVI in sub-regions II and V, while temperature and GTL were selected to model the NDVI in sub-region III.

Nested statistical models were developed for modeling NDVI in each sub-region. In most cases, the linear fitting of the NDVI time series to precipitation, temperature, and GTL resulted in poor correlation coefficients and high root means square error [17]. The vegetation growth may not be monotonic in relation to the atmospheric variables. For example, the temperature sensitivity of photosynthesis increases up to an optimum and later decreases as temperature gets higher [17]. A two-dimensional fitting was used to model NDVI using each of the hydro-meteorological variables as a factor. In order to select appropriate functional forms for modeling NDVI, four simple functional forms including linear function, quadratic function, exponent function and logarithm function were compared using Akaike information criterion (AIC) criteria. AIC is a measure of the relative quality of statistical models for a given set of data. The comparison of AIC values between different functional forms are as shown in

Table 2. The comparison of AIC values between different functional forms on 5 sub-regions; AIC in bold indicates most appropriate functional form.

NDVI and Temperature
	Region	Region I	Region II	Region III	Region IV	Region V
Function
liner		−3.6153		−4.7407	−6.1732
quadratic		−3.8733		−5.4122	−6.5508
exponent		−0.6537		−1.3973	−2.9187
logarithm		−3.7049		−5.2447	−6.5427
NDVI and Groundwater
	Region	Region I	Region II	Region III	Region IV	Region V
Function
liner			−1.8812	−3.8162		−2.7076
quadratic			−3.8876	−5.1042		−4.3212
exponent			−3.9708	−2.2641		−4.3329
logarithm			0.3744	−5.1401		0.4571
NDVI and Precipitation
	Region	Region I	Region II	Region III	Region IV	Region V
Function
liner		−3.5757	−1.7001		−4.9864	−2.8516
quadratic		−4.5808	−4.714		−7.1646	−4.4654
exponent		−4.5530	−4.5093		−7.1620	−4.4443
logarithm		−4.7402	−4.7974		−7.2847	−2.5677

. The smaller the AIC value, the more appropriate the form of the regression function. It was discovered that a logarithmic transformation of precipitation produced a better correlation with NDVI in some sub-regions (I, II, IV). Considering GTL and temperature, a negative exponential and a quadratic model showed a better relationship with NDVI. However, the logarithm regression model exhibited a better relationship between GTL and NDVI in sub-region III.

Different candidate equations were used to fit NDVI by coupling pair-wise precipitation, temperature, and GTL. The aim of the fitting procedures was to maximize the correlation coefficient and minimize the sum of squared residuals. However, the main importance of the fitting procedure was to maximize the proportion of the area in which the correlation coefficient met with significance testing.

In each sub-region, analysis was performed at the grid scale. In sub-regions I and IV, the coupling of precipitation and temperature resulted in a model of NDVI time series with an average correlation coefficient of 0.68 and 0.58, respectively. Li et al. (2016) also reported that in the desert regions in middle reaches of the Hei River, precipitation is a dominant influencing factor on annual NDVI variation [38]. Zhou et al. (2013) also demonstrated that precipitation and temperature are important factors driving the increase of NDVI in the study area [39]. However, the percentage of the area in which the correlation coefficient met with significance testing was only 42% (p < 0.05) and 15% (p < 0.01) in sub-region IV. This was mainly because the temperature was relatively high, thereby increasing evapotranspiration and water demand for vegetation growth that precipitation was unable to satisfy. In sub-regions II and V, the best fitting models were obtained by combining pair-wise negative exponentially transformed GTL respectively with the logarithmically and linearly transformed precipitation. In these regions, groundwater depths became an important influencing factor. Zhao et al. (2014) have also reported that in an oasis ecosystem, groundwater plays an important role in the formation of vegetation productivity [40]. In sub-region V, the percentage area in which the correlation coefficient met with the significance testing was only 50% (p < 0.05) and 16% (p < 0.01). There was a large area of arable land and artificial grassland in sub-region V (Figure 1). In sub-regions III and V, the growth of vegetation was strongly influenced by human activities, e.g., irrigation. Therefore, irrigation was the main influencing factor for vegetation growth. Beyond that, there was a considerable amount of water in sub-region III, thereby influencing prediction accuracy. After various fittings in different sub-regions, the correlation coefficients were combined for the entire study area (Figure 7). Results showed that there was about 53% of the area in which relationships were significant.

Analysis was also performed on a regional scale. Using the average NDVI time series of different sub-regions as predicted and each of the average hydro-meteorological variables as a predictor, a two-dimensional fitting was employed in the five sub-regions, respectively. This prediction also employed the same formulas used for the prediction on a grid scale (

Table 3. Fitting models, correlation coefficients and area proportion meeting significance testing in five sub-regions.

Sub-Region No.	Fitting Formula	Average Correlation Coefficient	Area Proportion	Area Proportion
			(p < 0.05)	(p < 0.01)
I	$NDV{I}_{T,P}=a\cdot \ln (b\cdot P)+c\cdot T^2+d\cdot T+e$	0.68	0.75	0.40
II	$NDV{I}_{GTL,P}=a\cdot e^{-(b\cdot GTL)}+c\cdot \ln (d\cdot P)+e$	0.66	0.72	0.42
III	$NDV{I}_{GTL,T}=a\cdot \ln (b\cdot GTL)+c\cdot T^2+d\cdot T+e$	0.47	0.36	0.16
IV	$NDV{I}_{T,P}=a\cdot \ln (b\cdot P)+c\cdot T^2+d\cdot T+e$	0.58	0.42	0.15
V	$NDV{I}_{GTL,P}=a\cdot GT{L}^2+b\cdot GTL+c\cdot P^2+d\cdot P+e$	0.57	0.50	0.16

). Figure 8a–e showed the results of the three-dimensional scatter-plot couplings. Comparisons of these figures demonstrated that sub-regions I and II had a better accuracy of fitting. Figure 8f presents the overall relationship between observed NDVI and modeled NDVI for the whole region, where the correlation coefficient was observed to be 0.97, which indicated a high accuracy. This showed that this method was suitable for NDVI modeling.

5. Conclusions

In this study, the spatio-temporal variation of vegetation cover in the middle reaches of theHei River basin and the main driving factors during the period from 2001 to 2011 were analyzed using MODIS NDVI data and hydro-meteorological data. Considering the nonlinear relationship between hydro-meteorological variables, the entropy theory was employed to calculate the mutual information between the NDVI time series and coupled hydro meteorological variables. Based on the spatial distribution of mutual information, the whole region was further divided into 5 sub-regions, and nested statistical models were employed to simulate NDVI in each sub-region. The main conclusions are as follows:

(1)

The average annual NDVI increased at a rate of 0.005/a over the past 11 years in the middle reaches of the Hei River. The percentage area in which NDVI increased occupied 94% of the total area.

(2)

In the desert sub-regions (I and IV), temperature and precipitation are the main driving factors for vegetation growth. In sub-region I, NDVI is consistent with the trend of temperature and precipitation (Figure 9a). However, in sub-region IV, the trend of temperature change is not obvious and the change in NDVI is mainly due to the increase in precipitation (Figure 9b). In the oasis regions (sub-region II and sub-region III), groundwater was an important factor for vegetation growth.

(3)

In coupling hydro-meteorological variables, a nested statistical model was proposed for modeling NDVI on a regional scale. The overall correlation coefficient between observed NDVI and modeled NDVI was observed to be 0.97. This high simulation accuracy further proves the suitability of this method.

(4)

Due to the influence of human activities, the modeling accuracy was not effective within the artificial oasis (sub-region III and sub-region V). For instance, in irrigation areas, vegetation can absorb water from irrigation to sustain growth but in non-irrigation areas the over-exploitation of groundwater caused by the increase in the amount of irrigation breaks the natural ecological balance, affecting the growth of natural vegetation. Additionally, the shortage of the dataset may also be a factor influencing the modeling accuracy. For instance, the scarcity of temperature and precipitation data may cause nearby cells to have similar values of T and P. Therefore, further studies are necessary for modeling NDVI that consider additional factors such as runoff and irrigation using long-term and higher resolution datasets.

The results can provide some scientific guidance for water resources management and the reasonable exploitation of groundwater, especially in the artificial oasis areas. The outcomes demonstrated the ability to successfully employ the joint entropy method in vegetation dynamics analysis in relation to hydro-meteorological variables.