Variation in Vegetation and Its Driving Force in the Middle Reaches of the Yangtze River in China

: The temporal and spatial characteristics of vegetation in the middle reaches of the Yangtze River (MRYR) were analyzed from 1999 to 2015 by trend analysis, co-integration analysis, partial correlation analysis, and spatial analysis using MODIS-NDVI time series remote sensing data. The average NDVI of the MRYR increased from 0.72 to 0.80, and nearly two-thirds of the vegetation showed a signiﬁcant trend of improvement. At the inter-annual scale, the relationship between NDVI and meteorological factors was not signiﬁcant in most areas. At the inter-monthly scale, NDVI was almost signiﬁcantly correlated with precipitation, relative humidity, and sunshine hours, and the effect of precipitation and sunshine hours on NDVI showed a pronounced lag. When the altitude was less than 2500 m, NDVI increased with elevation. NDVI increased gradually as the slope increased and decreased gradually as the slope aspect changed from north to south. NDVI decreased as the population density and per capita GDP increased and was signiﬁcantly positively correlated with afforestation policy. These ﬁndings provide new insights into the effects of climate change and human activities on vegetation growth.


Introduction
Natural ecosystems globally have become seriously threatened as climate change continues to intensify and the human population grows [1]. Vegetation, which generally refers to all of the plant communities covering the land surface, is one of the most important components of terrestrial ecosystems, as vegetation connects the atmosphere, soil, hydrology, and other ecological elements. The growth of vegetation is affected by both natural and human factors and is highly sensitive to environmental change [2,3]. Vegetation has long undergone dynamic changes; it thus plays an important role in the material cycle, energy flow, and ecosystem stability. Climate change and human activities have substantial effects on vegetation. Thus, study of the factors affecting changes in vegetation dynamics has implications for the maintenance of regional ecological balance [4,5]. The normalized difference vegetation index (NDVI) is an important indicator that reflects the growth status and dynamic changes in vegetation. It is highly sensitive to dynamic changes in surface vegetation (including vegetation coverage, biomass, and leaf area index) and has been widely used in vegetation monitoring; it has thus become a major focus of global change research [6][7][8].
The factors affecting vegetation cover change can be divided into natural factors and human factors. Previous studies have shown that climate factors may play a leading role at the regional scale, whereas human factors may play a leading role at the watershed scale [9]. In the Yangtze River Basin of China, precipitation and temperature are the most important factors affecting vegetation change [10]. Human activities can affect vegetation cover by modifying land use and land cover, which often forms a complex coupled humanenvironment system. In some traditional agricultural areas, the effect of human factors on vegetation cover requires consideration [11,12]. Vegetation cover change in the Yangtze River Basin is the result of the combined action of climate factors and human factors; however, human factors appear to be more important, and this is reflected in patterns of land use, related policies, and economic development [13].
Temporal and spatial patterns of vegetation change have been examined by various studies including changes in NDVI and the responses of NDVI to various factors [14][15][16]. For example, Coban determined that population, GDP, and topography have the greatest effect on vegetation change in the Western Black Sea region [17]. John and Ichii analyzed the relationship between NDVI and climate variables on a global scale and found that NDVI was significantly positively correlated with temperature in the middle and high latitudes of the northern region, but significantly negatively correlated with precipitation in the northern and southern semi-arid regions [15]. Tong et al. analyzed changes in vegetation cover from 2001 to 2013 and found that China's national and local governments implemented a series of ecological restoration projects to restore vegetation cover [18]. Zheng et al. noted that the ecological environment of the Yangtze River Delta urban agglomeration has improved over the past 20 years, but there are still significant differences between regions [19]. Yuan et al. and Tao et al. analyzed vegetation change in the Yangtze River Basin and found that NDVI may fluctuate in more than half of the study area in the future; they also found that altitude and annual average temperature were the main factors affecting changes in NDVI [20,21]. The responses of NDVI to different factors vary at different regional scales. Although several previous studies have conducted correlation analyses of vegetation change in the Yangtze River Basin, the changes and responses of vegetation to different factors in the middle reaches of the Yangtze River (MRYR) have not been well studied.

Study Area
The middle reaches of the Yangtze River (MRYR) include Hubei Province, Hunan Province, and Jiangxi Province (24 • 25 N-33 • 16 N, 108 • 24 E-118 • 23 E), with a total area of 564,600 km 2 ( Figure 1). This region features a subtropical monsoon climate, with annual precipitation of 1000-1600 mm·a −1 and annual average temperature of 16-18 • C [22]. The terrain of the MRYR mainly consists of mountainous hills and plains. The highest altitude is located in the mountainous area of Western Hubei, and the lowest is located in the central plain area. The mountainous hills account for more than 50% of the area, and the average altitude is 1497 m [23]. In 2018, the total population of the region was 174.38 million, accounting for 12.7% of the total population of the country; furthermore, the GDP was 9.8 trillion CNY (Chinese yuan), the urbanization level reached more than 50%, human economic activities were common, and regional disturbance was large [24]. The MRYR is rich in forest resources. In 2015, the forest coverage rates of Hubei Province, Hunan Province, and Jiangxi Province were 38.40%, 47.77%, and 60.01%, respectively [25][26][27].

Data Sources and Preprocessing
MODIS NDVI data from 1999 to 2015 derived from NASA's EOS/MODIS vegetation index product data MOD13A2 (https://wist.echo.nasa.gov/ap, accessed 17 April 2021) were used in this study. The spatial resolution of the data was 1 km × 1 km, and the temporal resolution was 16 d. ENVI5.3 was used to preprocess the original data, and the maximum NDVI of each year was obtained by the synthesis method. The meteorological data were derived from the National Meteorological Data Center (http://www.cma.gov.cn/2011qxfw/2011qsjgx/, accessed 25 July 2021). Monthly data (monthly average temperature, monthly precipitation, and monthly radiation) of 57 benchmark meteorological stations near the study area were from 1999 to 2015. In order to study the relationship between NDVI and terrain and socioeconomic factors, GDP and population, which are commonly used to represent socioeconomic factors as well as factors such as altitude, slope, and aspect used to represent terrain, were selected to analyze the relationship between NDVI and driving factors [28][29][30]. Terrain data (DEM data) and social data (population and GDP in 2015) were derived from the Resources and Environment Science and Data Center (http://www.resdc.cn, accessed 25 July 2021), which had a spatial resolution of 1 km × 1 km. Data for the artificial afforestation area were from the China Forestry Statistical Yearbook [31].

Theil-Sen Median Trend Analysis and Mann-Kendall Test
Theil-Sen median trend analysis was used to calculate the median slope between all n × (n − 1)/2 paired combinations of the time series data. The formula is as follows [32,33]: In Equation (1), β is the trend in vegetation change; NDVIj is the NDVI of the jth year; and NDVIi is the NDVI of the ith year. If β > 0, NDVI is increasing, indicating vegetation improvement or recovery. If β < 0, the NDVI is decreasing, indicating vegetation degradation.
The Mann-Kendall test can be used to assess the significance of the trend; this test has no specific requirements regarding the distribution of the samples, and the results are more accurate and not affected by outliers [34]. The formula is as follows: Set NDVIi, where i = 1999, 2000, …, 2015. The Z statistic can be defined as:

Data Sources and Preprocessing
MODIS NDVI data from 1999 to 2015 derived from NASA's EOS/MODIS vegetation index product data MOD13A2 (https://wist.echo.nasa.gov/ap, accessed on 17 April 2021) were used in this study. The spatial resolution of the data was 1 km × 1 km, and the temporal resolution was 16 d. ENVI5.3 was used to preprocess the original data, and the maximum NDVI of each year was obtained by the synthesis method. The meteorological data were derived from the National Meteorological Data Center (http://www.cma.gov.cn/2011qxfw/2011qsjgx/, accessed on 25 July 2021). Monthly data (monthly average temperature, monthly precipitation, and monthly radiation) of 57 benchmark meteorological stations near the study area were from 1999 to 2015. In order to study the relationship between NDVI and terrain and socioeconomic factors, GDP and population, which are commonly used to represent socioeconomic factors as well as factors such as altitude, slope, and aspect used to represent terrain, were selected to analyze the relationship between NDVI and driving factors [28][29][30]. Terrain data (DEM data) and social data (population and GDP in 2015) were derived from the Resources and Environment Science and Data Center (http://www.resdc.cn, accessed on 25 July 2021), which had a spatial resolution of 1 km × 1 km. Data for the artificial afforestation area were from the China Forestry Statistical Yearbook [31].

Theil-Sen Median Trend Analysis and Mann-Kendall Test
Theil-Sen median trend analysis was used to calculate the median slope between all n × (n − 1)/2 paired combinations of the time series data. The formula is as follows [32,33]: In Equation (1), β is the trend in vegetation change; NDVI j is the NDVI of the jth year; and NDVI i is the NDVI of the ith year. If β > 0, NDVI is increasing, indicating vegetation improvement or recovery. If β < 0, the NDVI is decreasing, indicating vegetation degradation.
The Mann-Kendall test can be used to assess the significance of the trend; this test has no specific requirements regarding the distribution of the samples, and the results are more accurate and not affected by outliers [34]. The formula is as follows: Set NDVI i , where i = 1999, 2000, . . . , 2015. The Z statistic can be defined as: In Equations (1)-(5), NDVI i and NDVI j represent the NDVI of pixel i and j years, respectively; n is the length of the time series; and sgn is a sign function. The value of the Z statistic is in the range of (−∞, +∞). When |Z| is greater than Z 1−α/2 at a given significance level α, there is a significant change in the time series at the α level. In this study, α = 0.05 was used to assess the significance of the pattern in regional NDVI change from 1999 to 2015 when |Z| > 1.96 at the confidence level of 0.05.

Stationarity Analysis and Co-Integration Analysis
If the time series factor does not pass the stationarity analysis and co-integration analysis, the situation of false regression might occur [35]. Therefore, we conducted a stationarity analysis for annual and monthly NDVI, temperature, precipitation, relative humidity, and sunshine hours, respectively. If the data were non-stationarity and each sequence was a single integral of the same order, the data would be analyzed for cointegration.

Stationarity Analysis
To analysis whether the time series being analyzed is stationarity, that is, to analysis whether there is a unit root, the commonly used method is ADF (Augmented Dickey-Fuller) analysis [35]. For the time series X = {x t }, the analysis equation is: where t is the time variable; α is the constant term; βt and β i are the trend items; ε t is the residual term; i is the lag order; and m is the maximum lag order. Assume that the null hypothesis H 0 : δ = 0, and the alternative hypothesis H 1 : δ = 0. The analysis was divided into three steps. In the first step, the analysis was carried out according to Equation (6). The second step was to delete the trend item for analysis. The third step was to delete the constant term and trend term for analysis. If H 0 is rejected in any step of the analysis, it means that there is no unit root in the time series X, that is, X is a stationarity time series, and the analysis can be stopped. Otherwise, the analysis should be continued until the end of the third step.

Co-Integration Analysis
The co-integration relation refers to the stationarity of the linear combination of two non-stationarity time series, so the co-integration analysis can be used to distinguish true regression from false regression [36]. The commonly used co-integration analysis method between two sequences is the E-G (Engle Granger) two-step method [37]. For the same order integral sequence X = {x t } and Y = {y t }, the analysis steps are as follows: The first step is to establish two regression equations of time series X and Y, and obtain the residual sequence {ε t } of the regression equation; and in step 2, the ADF analysis is performed on the residual sequence {ε t }. Assuming that the null hypothesis H 0 : δ = 0 and the alternative hypothesis H 1 : δ = 0, if the result rejects H 0 , it means that there is no unit root in the residual sequence {ε t }, that is, {ε t } is a stationarity time series, then the sequence X = {x t } and Y = {y t } can be considered to have a co-integration relationship.

Correlation Analysis
Correlation analysis and partial correlation analysis were used to analyze the relationship between NDVI and meteorological factors. The earth system is a complex system composed of multiple elements. Change in any element in the system can cause changes in other elements. Partial correlation analysis can effectively solve the above problems. When the correlation between two elements is studied separately, the effect of other elements is regarded as a constant by partial correlation [38]. The calculation formula is as follows: In Equation (6), r xy is the correlation coefficient between variables x and y; i is the number of samples; and x and y are the multi-year mean value of NDVI and the mean value of meteorological factors, respectively.
In Equation (7), r xy , r xz , and r yz are the correlation coefficients between elements x and y, x and z, and y and z, respectively. r xy·z is the partial correlation coefficient of factor X and y after factor Z is fixed. The t-test was used to test the significance of the partial correlation coefficients.

Relationship between NDVI and Driving Factors
The altitude gradient was divided into seven grades:

Temporal Variation in NDVI
In the time series of NDVI and meteorological factors, the first-order difference operator passed the stationarity analysis. There was a significant co-integration relationship between each variable and the change of time series. (Table A1, Table 1, Figures A1 and 2a).   T is temperature; P is precipitation; H is relative humidity; S is sunshine hours; (C, T, K) means that the analysis model contains intercept term, trend term, and lag order; Stationarity indicates that the variable passes the ADF stationarity analysis at a significance level of 5%.

Spatial Variation in NDVI
The spatial distribution of NDVI was based on the annual average NDVI data from 1999 to 2015 ( Figure 3). The western region is a mountainous area high in altitude. Natural forests, shrubs, and grasslands are widespread. Vegetation grew well, and NDVI was high. The northern region mainly consists of plains, as this is an important commodity grain production region in China, and NDVI was also high. Areas with low NDVI were mainly distributed in the capital cities of the three provinces (Changsha, Hunan Province; Wuhan, Hubei Province; and Nanchang, Jiangxi Province) and their adjacent urban areas. Most of these areas have experienced rapid economic development and substantial artificial disturbance. In 17 years, NDVI from 0.1-0.5 accounted for 1.40% of the area; NDVI from 0.5-0.7 accounted for 11.45%; NDVI from 0.7-0.8 accounted for 55.06%; and NDVI greater than 0.8 accounted for 32.09%.

Spatial Variation in NDVI
The spatial distribution of NDVI was based on the annual average NDVI data from 1999 to 2015 ( Figure 3). The western region is a mountainous area high in altitude. Natural forests, shrubs, and grasslands are widespread. Vegetation grew well, and NDVI was high. The northern region mainly consists of plains, as this is an important commodity grain production region in China, and NDVI was also high. Areas with low NDVI were mainly distributed in the capital cities of the three provinces (Changsha, Hunan Province; Wuhan, Hubei Province; and Nanchang, Jiangxi Province) and their adjacent urban areas. Most of these areas have experienced rapid economic development and substantial artificial disturbance. In 17 years, NDVI from 0.1-0.5 accounted for 1.40% of the area; NDVI from 0.5-0.7 accounted for 11.45%; NDVI from 0.7-0.8 accounted for 55.06%; and NDVI greater than 0.8 accounted for 32.09%.

Patterns of NDVI Change from 1999 to 2015
Several landscape patterns of NDVI in the MRYR from 1999 to 2015 were noted. Areas with pronounced vegetation improvement accounted for the largest proportion (67.39%) of the total area and were mainly distributed in the western and southeastern regions of the study area. Areas with slight improvement accounted for one-quarter (25.79%) of the total area and were mainly distributed in Jianghan Plain, Dongting Lake Plain, and Poyang Lake Plain. Areas with degradation accounted for 6.82% of the area. Areas with no change were denoted by stars and were scattered in the central area of major cities and outwards in the central area; the severely degraded locations were located in Wuhan, Changsha, and Nanchang, the capitals of the three provinces, and radiated outward, indicating that the expansion of urbanization had a significant impact on NDVI (Table 2 and Figure 4).

Patterns of NDVI Change from 1999 to 2015
Several landscape patterns of NDVI in the MRYR from 1999 to 2015 were noted. Areas with pronounced vegetation improvement accounted for the largest proportion (67.39%) of the total area and were mainly distributed in the western and southeastern regions of the study area. Areas with slight improvement accounted for one-quarter (25.79%) of the total area and were mainly distributed in Jianghan Plain, Dongting Lake Plain, and Poyang Lake Plain. Areas with degradation accounted for 6.82% of the area. Areas with no change were denoted by stars and were scattered in the central area of major cities and outwards in the central area; the severely degraded locations were located in Wuhan, Changsha, and Nanchang, the capitals of the three provinces, and radiated outward, indicating that the expansion of urbanization had a significant impact on NDVI (Table 2 and Figure 4).

Relationship between NDVI and Climate Change
In order to study the dynamic relationship between NDVI and meteorological factors, we needed to analyze the stationarity of NDVI, temperature, precipitation, relative humidity, and sunshine hours. All factors were stationarity after the first-order difference, so the cointegration could be measured (Table A1). Through the spatial analysis, there was a co integration relationship between the factors in some regions ( Table 3). The rejec-

Relationship between NDVI and Climate Change
In order to study the dynamic relationship between NDVI and meteorological factors, we needed to analyze the stationarity of NDVI, temperature, precipitation, relative humidity, and sunshine hours. All factors were stationarity after the first-order difference, so the cointegration could be measured (Table A1). Through the spatial analysis, there was a co integration relationship between the factors in some regions ( Table 3). The rejection ratio in stationarity analysis at the inter-monthly scale was 33.06%, 36.33%, 44.89%, and 45.60% higher than that at the inter-annual scale, respectively, indicating that there were more pixels at stationarity based on NDVI and meteorological factors at the inter-monthly scale ( Table 2). The rejection ratio of NDVI and temperature, NDVI and precipitation, NDVI and relative humidity, and NDVI and sunshine hours were higher than the acceptance ratio at the inter-monthly scale. The results indicated that there was a certain relationship between the NDVI and climate factors based on the vast majority of pixels in the MRYR at inter-monthly, and the correlation analysis could be conducted. Interannual variation in the correlation and partial correlations between NDVI and meteorological factors in the MRYR from 1999 to 2015 was analyzed (Table 4, Figure A1). NDVI was positively correlated with temperature and precipitation, with partial correlation coefficients of 0.488 and 0.277, respectively, and negatively correlated with relative humidity and sunshine hours, with partial correlation coefficients of −0.646 and −0.532, respectively. The correlation between NDVI and temperature was greater than that between NDVI and precipitation. The correlation coefficient and partial correlation coefficient between NDVI and climate factors were not significant (p > 0.05). The response of NDVI to climate change in the MRYR was not particularly strong. Note: T is temperature; P is precipitation; H is relative humidity; S is sunshine hours; R NDVI−T and R NDVI−T/PHS represents the correlation coefficient and partial correlation coefficient of NDVI and temperature, respectively, and so on.
Analysis between NDVI and climate across each of the 12 months from 1999 to 2015 (Table 5 and Figure A2) indicated that the growth of NDVI from January to April was significantly affected by precipitation and sunshine hours, and an obvious lag was noted in the effects of these variables. There was a significant correlation between NDVI in June and relative humidity in June. NDVI in July was affected by precipitation in May. There was a significant correlation between NDVI in August with temperature, relative humidity, and sunshine hours in August, and NDVI was most affected by climate factors. September was affected by the precipitation and relative humidity in September. There was no significant correlation between NDVI in November and various meteorological factors. NDVI in December was affected by the relative humidity and sunshine hours in December. Precipitation, relative humidity, and sunshine hours were the main climatic factors affecting seasonal variation in NDVI in the MRYR. Note: In the column of climate factors, T, P, H, and S are temperature, precipitation, relative humidity, and sunshine hours respectively. Jan., Feb., Mar., Apr., May., Jun., Jul., Aug., Sep., Oct., Nov. and Dec. are January, February, March, April, May, June, July, August, September, October, November, December. 0, 1, 2, and 3 denote the current month, the first month, the first two months, and the first three months, * denotes p < 0.05 and ** denotes p < 0.01.

Relationship between NDVI and Topographic Factors
NDVI increased from 0.7274 to 0.8512 (by 12%) at elevations less than 2500 m. At elevations greater than 2500 m, NDVI slightly decreased by 1.20% ( Figure 5 and Table 6). The proportion of areas with a slope less than 5 • in the study area was the highest (43.68%). The annual NDVI increased as each grade increased, but the rate of increase decreased as the grade increased, indicating that NDVI exhibited more pronounced changes at smaller slopes and that NDVI fluctuations were smaller at larger slopes. The study area was located in the Northern Hemisphere, and the north slope features both positive and negative slopes; NDVI was low (0.7437). The change in NDVI between other slope directions was not pronounced and remained at 0.76.
Water 2021, 13, x FOR PEER REVIEW 10 of 18 study area was located in the Northern Hemisphere, and the north slope features both positive and negative slopes; NDVI was low (0.7437). The change in NDVI between other slope directions was not pronounced and remained at 0.76.

Relationship of NDVI with Population and Economic Factors
NDVI decreased as the population density increased. When the population density increased from less than 500 persons/km 2 to more than 5000 people/km 2 , NDVI decreased from 0.8115 to 0.4071, a decrease of 49.83%. As the population density increased, the decrease in NDVI also increased gradually. When the population density reached 5000 persons/km 2 , the decline in NDVI was 39.24%. The overall results indicated that the population density had a significant negative effect on the vegetation cover in the study area, and  NDVI decreased as the population density increased. When the population density increased from less than 500 persons/km 2 to more than 5000 people/km 2 , NDVI decreased from 0.8115 to 0.4071, a decrease of 49.83%. As the population density increased, the decrease in NDVI also increased gradually. When the population density reached 5000 persons/km 2 , the decline in NDVI was 39.24%. The overall results indicated that the population density had a significant negative effect on the vegetation cover in the study area, and this effect strengthened as the population density increased ( Figure 6).
NDVI of the study area decreased as per capita GDP increased. When the per capita GDP increased from less than 500 CNY to more than 5000 CNY, NDVI decreased from 0.8422 to 0.6281, a decrease of 25.42%. NDVI decreased the most when the per capita GDP increased to more than 5000 CNY, indicating that economic development of the city center had a strong negative effect on vegetation. Per capita GDP also has an inhibitory effect on vegetation cover in the study area ( Figure 6 and Table 7).
Water 2021, 13, x FOR PEER REVIEW 11 of 18 NDVI of the study area decreased as per capita GDP increased. When the per capita GDP increased from less than 500 CNY to more than 5000 CNY, NDVI decreased from 0.8422 to 0.6281, a decrease of 25.42%. NDVI decreased the most when the per capita GDP increased to more than 5000 CNY, indicating that economic development of the city center had a strong negative effect on vegetation. Per capita GDP also has an inhibitory effect on vegetation cover in the study area ( Figure 6 and Table 7).

Relationship between NDVI and Policy Factors
During the study period, NDVI gradually increased as the artificial afforestation area increased, and there was a significant correlation between the two (r = 0.6828, n = 17, p < 0.001). In 2002 and 2003, the afforestation area was increased, and the growth rate of NDVI was high. Between 2004 and 2006, the afforestation area in the study area declined due to policy-oriented factors, but in general, afforestation policy has improved the coverage rate of forestland and increased the coverage area of forestland. During this period, some early afforestation was not appropriate, and the low efficiency plantations formed also appeared to have the degradation and death phenomenon [39,40]. Therefore, NDVI presented a relatively stable state with little change. The decrease in NDVI in 2011 stemmed from the effect of drought in that year. Overall, ecological projects such as returning farmland to forest have facilitated gradual increases in NDVI in the MRYR, and the condition of the vegetation has greatly improved (Figure 7).

Relationship between NDVI and Policy Factors
During the study period, NDVI gradually increased as the artificial afforestation area increased, and there was a significant correlation between the two (r = 0.6828, n = 17, p < 0.001). In 2002 and 2003, the afforestation area was increased, and the growth rate of NDVI was high. Between 2004 and 2006, the afforestation area in the study area declined due to policy-oriented factors, but in general, afforestation policy has improved the coverage rate of forestland and increased the coverage area of forestland. During this period, some early afforestation was not appropriate, and the low efficiency plantations formed also appeared to have the degradation and death phenomenon [39,40]. Therefore, NDVI presented a relatively stable state with little change. The decrease in NDVI in 2011 stemmed from the effect of drought in that year. Overall, ecological projects such as returning farmland to forest have facilitated gradual increases in NDVI in the MRYR, and the condition of the vegetation has greatly improved (Figure 7).

Variation in NDVI and Its Relationship with Climatic Factors
From 1999 to 2015, NDVI in the MRYR increased, which was consistent with the results of previous research in Hubei Province [41,42], Hunan Province [43], Jiangxi Province [44], and the Yangtze River Basin [20]. Average annual NDVI of the Yangtze River Delta, Pearl River Delta, and other coastal urban agglomerations were mostly between 0.3 and 0.5. Compared with these regions, the average NDVI in the MRYR (0.76) is relatively high (Figures 2 and 3). Approximately 67.39% of the regions showed a significant improvement in NDVI, and these regions were mainly distributed in non-urban areas. Stable and slightly improved areas accounted for one-quarter of the regional area (25.79%) and were mainly distributed in Jianghan Plain, Dongting Lake Plain, and Poyang Lake Plain. As these three plains are the main centers of commodity grain production in China (Table  1), the slight changes in NDVI observed are consistent with the slight changes observed in crops [45][46][47]. Degraded areas accounted for 6.82% of the total area, and Wuhan, Changsha, and Nanchang made up most of the degraded area; this indicates that the expansion of urbanization has a significant effect on vegetation coverage.
There was almost no significant relationship between NDVI and climate factors (Tables 3 and 4). On the monthly and annual scale, there was a 2-3-month lag in the effect of precipitation on vegetation growth, a 2-month lag for sunshine hours, and a 3-month lag for relative humidity, which indicates that the change in NDVI is also affected by soil water and nutrient availability (Tables 3 and 5). This is consistent with the results of global, regional, and other multi-scale studies; the effect of vegetation cover on climate change also shows a lag [36].

Relationship between Dynamic Change in NDVI and Terrain Factors
When the elevation was less than 2500 m, NDVI increased with elevation, but when the elevation was more than 2500 m, NDVI decreased ( Table 6). The explanation for this pattern is the greater intensity of human activities at low elevation. As the elevation increases, the effect of human activities on vegetation weakens, and NDVI increases [48][49][50]. As the elevational gradient increases, the rate of NDVI growth decreases, which may stem from changes in natural environmental factors such as temperature, precipitation, and light [51][52][53]. Higher elevation corresponds to lower soil temperature, higher humidity, and greater vulnerability to leaching, which leads to soil acidification and inhibits the growth of vegetation [54]. As the terrain slope increased, the NDVI of different slope grades gradually increased. This result is consistent with Lu showing that NDVI increases

Variation in NDVI and Its Relationship with Climatic Factors
From 1999 to 2015, NDVI in the MRYR increased, which was consistent with the results of previous research in Hubei Province [41,42], Hunan Province [43], Jiangxi Province [44], and the Yangtze River Basin [20]. Average annual NDVI of the Yangtze River Delta, Pearl River Delta, and other coastal urban agglomerations were mostly between 0.3 and 0.5. Compared with these regions, the average NDVI in the MRYR (0.76) is relatively high (Figures 2 and 3). Approximately 67.39% of the regions showed a significant improvement in NDVI, and these regions were mainly distributed in non-urban areas. Stable and slightly improved areas accounted for one-quarter of the regional area (25.79%) and were mainly distributed in Jianghan Plain, Dongting Lake Plain, and Poyang Lake Plain. As these three plains are the main centers of commodity grain production in China (Table 1), the slight changes in NDVI observed are consistent with the slight changes observed in crops [45][46][47]. Degraded areas accounted for 6.82% of the total area, and Wuhan, Changsha, and Nanchang made up most of the degraded area; this indicates that the expansion of urbanization has a significant effect on vegetation coverage.
There was almost no significant relationship between NDVI and climate factors (Tables 3 and 4). On the monthly and annual scale, there was a 2-3-month lag in the effect of precipitation on vegetation growth, a 2-month lag for sunshine hours, and a 3-month lag for relative humidity, which indicates that the change in NDVI is also affected by soil water and nutrient availability (Tables 3 and 5). This is consistent with the results of global, regional, and other multi-scale studies; the effect of vegetation cover on climate change also shows a lag [36].

Relationship between Dynamic Change in NDVI and Terrain Factors
When the elevation was less than 2500 m, NDVI increased with elevation, but when the elevation was more than 2500 m, NDVI decreased ( Table 6). The explanation for this pattern is the greater intensity of human activities at low elevation. As the elevation increases, the effect of human activities on vegetation weakens, and NDVI increases [48][49][50]. As the elevational gradient increases, the rate of NDVI growth decreases, which may stem from changes in natural environmental factors such as temperature, precipitation, and light [51][52][53]. Higher elevation corresponds to lower soil temperature, higher humidity, and greater vulnerability to leaching, which leads to soil acidification and inhibits the growth of vegetation [54]. As the terrain slope increased, the NDVI of different slope grades gradually increased. This result is consistent with Lu showing that NDVI increases as the slope increases in the Yangtze River Basin [55]. As the slope grade increased, the average annual NDVI of vegetation increased gradually, but the rate of increase reduced at higher slope grades, which indicated that the change in NDVI was more obvious when the slope was small and that the fluctuation in NDVI was smaller when the slope was large. NDVI was higher on sunny slopes than on shaded slopes.

Relationship between Dynamic Change in NDVI and Socio-Economic Factors
Dynamic changes in vegetation in the MRYR are not only affected by natural factors such as climate and topography, but also by human activities [56,57]. This region is the cradle of human civilization, as it consists of rivers and lakes and a fertile alluvial plain, making it a center of human agricultural production and social life [58]. During the study period, NDVI significantly decreased as the population density and per capita GDP increased, and changes in NDVI with population density were more pronounced. When the population density increased from less than 500 people/km 2 to more than 5000 people/km 2 , NDVI decreased by 39.24%, and when the per capita GDP increased from less than 500 CNY to more than 5000 CNY, NDVI decreased by 16.60% (Table 7). With the rapid growth of the population and GDP, many people from rural and suburban regions have moved to cities. To accommodate a larger population, cities have further expanded and occupied several high-quality land types such as cultivated land, forest land, and unused land, which leads to a gradual reduction in regional vegetation coverage and has a substantial effect on the ecological environment. This is consistent with the results of other studies [59][60][61].
To cope with global climate change and problems associated with the ecological environment, a series of ecological construction projects have been implemented at the national level over the past 30 years such as the natural forest resources protection project, conversion of farmland to forest project, national desertification control, protection forest construction project in the middle and upper reaches of the Yangtze River, and nature reserve construction [62][63][64][65]. The MRYR plays an important ecological role, as it is China's economic center and main grain production area [66,67]. At the regional and national levels, there is much focus on regional vegetation construction, and several ecological projects and national ecological policies have been implemented. These projects have played a major role in regional vegetation restoration. As the artificial afforestation area increased, NDVI gradually increased, and there was a significant correlation between NDVI and artificial afforestation area (r = 0.6828, n = 17, p < 0.001) (Figure 7). In 2002 and 2003, the afforestation area and the growth rate of NDVI were high. The construction of urban buildings, roads, and other structures that accompany economic development inevitably damages vegetation. However, this development is also accompanied by efforts to beautify the environment, as reflected by afforested areas, the conversion of farmland to forests, and the construction of green space, parks, and ecological vegetation in the city. Both of these reflect the effect of human activities on vegetation cover in the study area, and dynamic changes in vegetation are becoming increasingly affected by human production and life.

Conclusions
NDVI increased from 1999 to 2015, and nearly two-thirds of the vegetation showed a significant trend of improvement in the MRYR over this period. At the inter-annual scale, the relationship between NDVI and meteorological factors was not significant in most areas. At the inter-monthly scale, precipitation, relative humidity, and sunshine hours were the main factors affecting changes in NDVI, and the effect of precipitation and sunshine hours on NDVI showed a pronounced lag. NDVI increased with elevation when the altitude was less than 2500 m and decreased when the altitude was greater than 2500 m. As the slope increased, NDVI gradually increased. NDVI gradually decreased when the slope aspect changed from south to north. Per capita GDP and population density had a significant negative effect on NDVI. There was a significant correlation between ecological projects (afforestation area) and NDVI. The ecological afforestation project implemented by the state improved the regional vegetation. The dynamic changes in NDVI in the study area are thus jointly affected by climate change and human impacts.  Data Availability Statement: Some or all data, models, or code generated or used during the study are available from the author by request (yiyang0307@sjtu.edu.cn).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Water 2021, 13, x FOR PEER REVIEW 14 of 18 a significant negative effect on NDVI. There was a significant correlation between ecological projects (afforestation area) and NDVI. The ecological afforestation project implemented by the state improved the regional vegetation. The dynamic changes in NDVI in the study area are thus jointly affected by climate change and human impacts. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Some or all data, models, or code generated or used during the study are available from the author by request (yiyang0307@sjtu.edu.cn).

Conflicts of Interest:
The authors declare no conflicts of interest.
Appendix A Figure A1. Trends in temperature changes in the middle reaches of the Yangtze River from 1999 to 2015.   Yes Note: T is temperature; P is precipitation; H is relative humidity; S is sunshine hours; ∆ represents the first difference of the variables; (C, T, K) means that the analysis model contains intercept term, trend term, and lag order; Stationarity indicates that the variable passes the ADF stationarity analysis at a significance level of 5%.