The Dynamic Change of Vegetation Cover and Associated Driving Forces in Nanxiong Basin , China

Natural climate change and human activities are the main driving forces associated with vegetation coverage change. Nanxiong Basin is a key ecosystem-service area at the national level with a dense population and highly representative of red-bed basins, which are considered as fragile ecological units in humid regions. In this study, the authors aimed to determine the trends in vegetation cover change over past two decades and the associated driving forces in this study area. The Normalized Difference Vegetation Index (NDVI) of 2000–2015, derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) remote sensing dataset along with the application of statistical methods and GIS (geographic information system) techniques were used to quantify vegetation cover change. The results show that human-induced factors can explain most variations at sites with significant cover change. That is to say that human activities are the main drivers of vegetation dynamics in this study area, which shows a significant reduction trend in vegetation cover during the industrialization and urbanization processes of the study period and noticeable recovery trend in 2000–2015 under the plantation and enclosed forest policy.


Introduction
Vegetation has a considerable impact on surficial energy exchange processes and acts as an interface between land and atmosphere.It affects local and regional climate [1,2] and the hydrologic balance at the surface.Vegetation dynamics have been recognized to be of primary importance in the global change of terrestrial ecosystems [3,4].
Plants take part in maintaining a natural ecosystem's health.They are greatly sensitive to changes in environmental elements, and vegetation is always considered as an important indicator for measuring local environmental conditions and detecting environmental change [5].Thus, vegetation cover change is gaining recognition as a key component of global environment change [6].Changes in vegetation coverage rank as one of the key components of vegetation dynamics, which can help to assess the environmental effects of vegetation cover change more clearly [7].In recent years, with global climate change along with increases in population and resource consumption, ecological degradation has become a major environmental problem that influences human production and quality of life.Thus, the monitoring of vegetation change is of great significance.
Ground-based measurements are difficult to use for detecting and predicting vegetation dynamics on a regional to global scale because such data are collected at small spatial and temporal scales and vary with respect to type and reliability [8].The remotely sensed Normalized Difference Vegetation Index (NDVI) can indicate growth status, phenology and vegetation cover over large spatial extents and at a finer temporal resolution [9][10][11].It has been widely used to monitor vegetation dynamics and vegetation responses to human factors [12,13].
Due to fine temporal resolutions and timely updates, the time series of remotely sensed data are widely used in environmental-change detection.In recent decades, with the application of remote sensing and geographic information science (GIS) techniques, numerous case studies have been conducted in order to monitor and map changes in vegetation coverage at various spatial scales, such as global [14][15][16][17], national [13,[18][19][20][21], provincial [22][23][24][25][26][27], and also county scales [7,28].However, few studies have been conducted for a complete natural geographical unit, i.e., a basin.In addition, long-term change monitoring may average out multiple short-term changes [26], and, therefore fails to detect the change processes at fine temporal intervals [7].On the contrary, short-term change detection in recent decades can clearly convey environmental change.
China is facing many challenges that affect the sustainable development of society, including an ever-growing demand on the environment and natural resources.Past governmental policies have favored economic growth over environmental protection, which did not improve until the 1998 floods stimulated large-scale conservation programs (e.g., natural forest conservation and grain-to-green programs).With the implementation of a series of nation-wide environmental protection measurements, it is important to assess their effectiveness, especially in areas with dense population where the ecological environment is more likely to be affected by human behavior.
Nanxiong Basin is a red-bed basin with a severe soil erosion problem due to its dominant purple soil texture (Calcaric Regosols in FAO taxonomy) and the fact that red clay-rich sediments of the parent rock are mainly distributed in the central part of the basin.This basin is representative of the purple soil regions in many aspects, including water and soil erosion conditions, land degradation, material circulation of farmland ecosystems, regulation and slope runoff.However, little is known about the overall ecosystem of these areas because of limited and fragmented data.Therefore, we choose Nanxiong Basin as our study area and analyze its spatiotemporal variations in vegetation coverage in order to assess ecological change and further understand the exerted influence of human factors on these typical ecosystems of South China.
Overall, the goals of this study are to monitor and map the spatiotemporal dynamics of vegetation coverage in Nanxiong Basin from 2000-2015.The vegetation coverage is indicated by the remotely sensed NDVI dataset, and its temporal trend is detected by statistical methods including trend analysis, overlay analysis and the correlation analysis.In addition, GIS techniques were applied to detect the associated driving forces of vegetation coverage change in combination with the statistical methods.

Description of the Study Area
Nanxiong Basin (24 • 33 -25 • 24 N, 113 • 52 -114 • 45 E) is located in the north of Guangdong Province, with an area of 8848 km 2 (Figure 1).The elevation ranges from 48 to 1421 m above sea level (ASL).A subtropical monsoon climate prevails, with long hot summers and short winters.Based on the meteorological observations collected from Nanxiong Station (1956-2010), the mean annual temperature of the study basin is 19.6  2b) mainly include woodland (60.4%), shrubland (13.1%), cultivated land (14.8%), construction land (3.1%), garden plots (1.9%) and grassland (1.4%).Purple soil accompanies the red-bed parent material distributed in the central part of the basin.Nanxiong City functions as a key ecosystem-service area at the national level and is also one of the state-level poverty counties.

MODIS Time Series of NDVI
The NDVI time series from 18 February 2000 to 17 November 2015 were acquired by a moderate resolution imaging spectrometer (MODIS) Terra Satellite (MOD13Q1, https://reverb.echo.nasa.gov/reverb/).This 16-day composite time series has a spatial resolution of 250 m.To ensure that only reliable observations were used in the analysis, quality assurance (QA) datasets with 0 in Rank Key pixel reliability were used in order to exclude observations with cloudy and/or snow/ice conditions.In addition, in order to be consistent with other input data, the NDVI dataset was reprojected into the Elipse Lambert Azimuthal Equal Area coordinate system.
The MODIS sensor has a spatial and temporal resolution that allows for monitoring and mapping of vegetation at a regional level.With minimal variations associated with external influences (atmosphere, view and sun angles), the MODIS vegetation index product provides consistent spatial and temporal comparisons of vegetation conditions [30].Although these types of images do not have a high spatial resolution, it is possible to have cloud-free images because of the high revisit frequency of the satellite (1-2 days).Higher spatial resolution satellites are usually only available a few times per year for a given location [31].This is due to the lower revisit frequency and cloud cover frequency [32,33].NDVI data derived from Landsat sensors (e.g., TM/ETM + OLI with 30-m resolution at visible and near infrared bands) can successfully capture the spatial details and alleviate the spectral mixture issue to some extent [34,35].However, longer observation cycles of these satellites and frequent cloud contamination in some regions limit their application for detecting rapid changes in ecosystems [36,37].The authors, therefore, did not use Landsat images in this study.

Rainfall
The daily time series of rainfall from a total of 16 National Meteorological Observing Stations, covering January 1980 to December 2013, were obtained from the National Meteorological Information Center, China Meteorological Administration [29].Nanxiong Station is one of the national-level stations located inside the study area with coordinates of 25°08′N and 114°19′E at 1338 m ASL.To obtain the spatial distribution of rainfall within the study basin, monthly precipitation data were interpolated using the Kriging method at the same resolution as NDVI dataset.

MODIS Time Series of NDVI
The NDVI time series from 18 February 2000 to 17 November 2015 were acquired by a moderate resolution imaging spectrometer (MODIS) Terra Satellite (MOD13Q1, https://reverb.echo.nasa.gov/reverb/).This 16-day composite time series has a spatial resolution of 250 m.To ensure that only reliable observations were used in the analysis, quality assurance (QA) datasets with 0 in Rank Key pixel reliability were used in order to exclude observations with cloudy and/or snow/ice conditions.In addition, in order to be consistent with other input data, the NDVI dataset was reprojected into the Elipse Lambert Azimuthal Equal Area coordinate system.
The MODIS sensor has a spatial and temporal resolution that allows for monitoring and mapping of vegetation at a regional level.With minimal variations associated with external influences (atmosphere, view and sun angles), the MODIS vegetation index product provides consistent spatial and temporal comparisons of vegetation conditions [30].Although these types of images do not have a high spatial resolution, it is possible to have cloud-free images because of the high revisit frequency of the satellite (1-2 days).Higher spatial resolution satellites are usually only available a few times per year for a given location [31].This is due to the lower revisit frequency and cloud cover frequency [32,33].NDVI data derived from Landsat sensors (e.g., TM/ETM + OLI with 30-m resolution at visible and near infrared bands) can successfully capture the spatial details and alleviate the spectral mixture issue to some extent [34,35].However, longer observation cycles of these satellites and frequent cloud contamination in some regions limit their application for detecting rapid changes in ecosystems [36,37].The authors, therefore, did not use Landsat images in this study.

Rainfall
The daily time series of rainfall from a total of 16 National Meteorological Observing Stations, covering January 1980 to December 2013, were obtained from the National Meteorological Information Center, China Meteorological Administration [29].Nanxiong Station is one of the national-level stations located inside the study area with coordinates of 25 • 08 N and 114 • 19 E at 1338 m ASL.To obtain the spatial distribution of rainfall within the study basin, monthly precipitation data were interpolated using the Kriging method at the same resolution as NDVI dataset.

Quantifying Trends in NDVI
The authors used MATLAB ® to calculate the following statistics from the MOD13Q1 NDVI time series per pixel: a monotonic trend based on the nonparametric M-K τ correlation coefficient (ranging between −1 and 1, which is simply the relative frequency of increases minus the relative frequency of decreases) and M-K significance, which calculates the z and p values of the significance of the monotonic trend.In this study, the default significance level was set as 0.05.In addition, in order to reduce the impact of outliers on the results, the authors used nonlinear statistical estimates of trends, as NDVI trends may not be linear [38].These statistics were calculated for the entire time series covering 16 years (with a total of 363 images).
The aim was to examine the drivers of vegetation change, so the authors digitized a sample of areas throughout Nanxiong Basin that were subsequently analyzed in detail.These regions of interest (ROIs) were manually chosen with the aim of representing locations with negative and positive trends in NDVI.In addition, the distribution of the selected ROIs should cover the entire area of Nanxiong Basin as random as possible.To this end, the authors digitized 82 ROIs, where statistically significant trends in NDVI were identified in these ROIs (Figure 2d).The authors had no prior knowledge regarding the type of land cover changes driving NDVI trends in the chosen ROIs (Figure 2c).An average nearest neighbor analysis (a global Moran's I spatial autocorrelation test) indicated that the spatial pattern of these 82 ROIs corresponded to a clustered pattern (z = −1.69,p = 0.09), which was controlled by the clustered pattern of pixels with significant change.The area of these ROIs ranged between three and 103 pixels, with a median area size of six pixels (i.e., between 0.19 km 2 and 6.44 km 2 , with a median area size of 0.38 km 2 ).The overall area of these ROIs covered 53.44 km 2 , among which about 31.4% of the area showed statistical significance in M-K τ over 16 years (i.e., about 0.6% of the entire study area).It is believed that these ROIs are adequate for generalizing the study outputs of the entire Nanxiong Basin.
The cumulative sum (CUSUM) analysis is a sequential analysis technique that is typically used for monitoring change detection [39].The authors used CUSUM analysis of NDVI values to identify the years when an event, typically a change on average, took place (i.e., when a change in NDVI values took place in each of the 82 ROIs).The principle behind the CUSUM analysis is that it displays a running total of the differences from an average [40].Changes in direction of the obtained line from CUSUM analysis indicate that an event has occurred.

Statistical Analyses
For each of the 82 ROIs, the authors calculated the minimum, maximum, standard deviation, mean and coefficient of variation (CV) values of NDVI over the entire 16-year time series.In addition, the authors identified the land cover/land-use types in 2000 and 2015 respectively (i.e., the beginning and end of the NDVI time series).The data used include available images on Google Earth, orthophotos, and land-use cover data of 2000 and 2009, representing urbanized areas, agricultural areas, planted forests and water bodies (from the Department of Land resources of Guangdong Province).Following this, fine spatial resolution imagery and a detailed land-use map (instead of classifying Landsat images for accuracy consideration) were linked and visually examined.For most ROIs, the observed changes in NDVI can be explained and further categorized as: (1) plantation and recovery; (2) forest-cutting; (3) construction; and (4) forest-to-terrace.Due to the saturation phenomenon in NDVI [41], mistakes may occur in the identification of the forest-cutting category when the canopy cover is greater than 80%.However, deforestation generally caused huge decreases in NDVI values and as a result, there was no impact of NDVI saturation on the identification of the forest-cutting category.
With the goal of detecting the cause behind observed changes in NDVI from 2000 to 2015, the authors recorded forest gain and loss and ignored any subsequent afforestation that occurred during the study period.Forest-cutting as a deforestation event was defined in the study period.The authors used three-year moving windows to adjust deforestation events only if the mean NDVI of the first year of a window was higher than the next two years by 20% (considering that NDVI saturation occurs when canopy cover is >80% [41]) and no conversion of forestland to another type of land use was recorded in the 2000 and 2015 layers.For plantation and recovery, a three-year moving window was also used to detect plantation and recovery only if all the mean NDVI values of the third year were not lower than the previous two years, and no conversion of forestland to another type of land use was recorded in 2000 and 2015.The third and fourth categories (i.e., construction and forest-to-terrace) were easy to detect by comparing land-use types between 2000 and 2015.For arid areas where standard deviation of NDVI and mean NDVI are very low (<0.01 and <0.13 respectively), the precision of obtained NDVI dataset is within error bounds, and thus the NDVI trends could not be reliably determined in these areas [42,43]; however, in this study area with humid climate, this did not need to be taken into consideration.
In order to examine the correspondence between the NDVI and rainfall time series (of the 82 ROIs), the authors calculated three types of correlations.The first type of correlation is between monthly rainfall values and monthly NDVI values.The second correlation is between cumulative monthly rainfall values (ranging between 1 and 6 months) and monthly NDVI values.This allowed the authors to account for the time lag in the response of vegetation growth to rainfall [37].The third correlation is between annual rainfall values and mean annual NDVI values, which removed the seasonal influence of rainfall on plant phenology.

Vegetation Dynamics Observed from Multiyear Mean NDVI
Across the entire study area, the average value of the multiyear mean NDVI (NDVI m ) from 2000 to 2015 was 0.65 (Figure 2c).The areas with an NDVI m > 0.7 account for 24% of the basin and are mainly comprised of woodland cover.
The areas with NDVI m < 0.6 make up 14% of the basin and are distributed mainly in the central parts.The average value of CV for the annual mean NDVI was 0.056.In addition, the spatial patterns of CV and NDVI m show a negative correlation.That is to say, areas with lower NDVI m have a higher CV value and vice versa (Figure 2c,d).Furthermore, the higher the NDVI value, the lower the fluctuation of vegetation greenness.Most of the areas with relatively lower NDVI were around the city due to urban expansion in the recent decade.Moreover, since irrigation could not be guaranteed in the study area, NDVI fluctuation was also more significant for cultivated land and garden plots around Nanxiong Basin (Figure 3).was recorded in the 2000 and 2015 layers.For plantation and recovery, a three-year moving window was also used to detect plantation and recovery only if all the mean NDVI values of the third year were not lower than the previous two years, and no conversion of forestland to another type of land use was recorded in 2000 and 2015.The third and fourth categories (i.e., construction and forest-toterrace) were easy to detect by comparing land-use types between 2000 and 2015.For arid areas where standard deviation of NDVI and mean NDVI are very low (<0.01 and <0.13 respectively), the precision of obtained NDVI dataset is within error bounds, and thus the NDVI trends could not be reliably determined in these areas [42,43]; however, in this study area with humid climate, this did not need to be taken into consideration.
In order to examine the correspondence between the NDVI and rainfall time series (of the 82 ROIs), the authors calculated three types of correlations.The first type of correlation is between monthly rainfall values and monthly NDVI values.The second correlation is between cumulative monthly rainfall values (ranging between 1 and 6 months) and monthly NDVI values.This allowed the authors to account for the time lag in the response of vegetation growth to rainfall [37].The third correlation is between annual rainfall values and mean annual NDVI values, which removed the seasonal influence of rainfall on plant phenology.

Vegetation Dynamics Observed from Multiyear Mean NDVI
Across the entire study area, the average value of the multiyear mean NDVI (NDVIm) from 2000 to 2015 was 0.65 (Figure 2c).The areas with an NDVIm > 0.7 account for 24% of the basin and are mainly comprised of woodland cover.
The areas with NDVIm < 0.6 make up 14% of the basin and are distributed mainly in the central parts.The average value of CV for the annual mean NDVI was 0.056.In addition, the spatial patterns of CV and NDVIm show a negative correlation.That is to say, areas with lower NDVIm have a higher CV value and vice versa (Figure 2c,d).Furthermore, the higher the NDVI value, the lower the fluctuation of vegetation greenness.Most of the areas with relatively lower NDVI were around the city due to urban expansion in the recent decade.Moreover, since irrigation could not be guaranteed in the study area, NDVI fluctuation was also more significant for cultivated land and garden plots around Nanxiong Basin (Figure 3).

Vegetation Cover Change in ROIs and Driving Forces
Over the study period of 16 years, statistically significant (p = 0.05) trends in NDVI values were found in 76% of study area using M-K statistical analysis, with 9% of the study area exhibiting negative trends and 67% positive trends (Figure 2a).Out of the selected 82 ROIs, 35  Previous research suggests that a lag effect exists in the response of monthly NDVI to rainfall [44].For most ROIs, maximum correlation between monthly rainfall amounts and monthly mean NDVI were found when using cumulative rainfall over a period of 4 months.The average correlation coefficient was 0.60 compared to 0.39 when the time lag response of NDVI values to rainfall was not considered (Figure 4).When examining the correlation between annual rainfall amounts and mean annual NDVI, no positive correlation values above 0.5 were found (Figure 5).

Vegetation Cover Change in ROIs and Driving Forces
Over the study period of 16 years, statistically significant (p = 0.05) trends in NDVI values were found in 76% of study area using M-K statistical analysis, with 9% of the study area exhibiting negative trends and 67% positive trends (Figure 2a).Out of the selected 82 ROIs, 35  Previous research suggests that a lag effect exists in the response of monthly NDVI to rainfall [44].For most ROIs, maximum correlation between monthly rainfall amounts and monthly mean NDVI were found when using cumulative rainfall over a period of 4 months.The average correlation coefficient was 0.60 compared to 0.39 when the time lag response of NDVI values to rainfall was not considered (Figure 4).When examining the correlation between annual rainfall amounts and mean annual NDVI, no positive correlation values above 0.5 were found (Figure 5).The four categories of land cover conversion for the 82 ROIs were: plantation and recovery (35); forest-cutting (21); construction (13); and forest-to-terrace (13).The factors driving land cover changes in the classes of agriculture and construction dynamics were associated with both positive and negative trends in NDVI values.Trends of NDVI in ROIs derived from M-K τ were not correlated with rainfall, mean NDVI and standard deviation or range of NDVI values (Figure 4).Based on the CUSUM analysis, NDVI change can be identified in 39 of the 82 ROIs (48%), and the change mostly occurred between 2008 and 2014 (Figure 6).
NDVI trends driven by agricultural processes were identified in 62 ROIs, of which 35 were associated with an increase in NDVI values and 27 exhibited a decrease in NDVI values.NDVI values within agricultural areas (plantation and recovery, forest-to-terrace and forest-cutting) were often associated with relatively high values of mean, standard deviation and range (Figure 5).The oscillations related to agricultural practices belong to three major subclasses: (i) Conversion of a shrubland or woodland into a terrace or fruit forest or tea yard, resulting in a decrease in NDVI values and, consequently, low standard deviation or range of NDVI (Figure 5).(ii) The recovery and planting of shrubbery and forest resulted in an increase in NDVI values.A negative correlation between standard deviation and the M-K τ value can be found.(iii) The cutting of forests resulted in decreasing NDVI values, as shown in Figure 7.The four categories of land cover conversion for the 82 ROIs were: plantation and recovery (35); forest-cutting (21); construction (13); and forest-to-terrace (13).The factors driving land cover changes in the classes of agriculture and construction dynamics were associated with both positive and negative trends in NDVI values.Trends of NDVI in ROIs derived from M-K τ were not correlated with rainfall, mean NDVI and standard deviation or range of NDVI values (Figure 4).Based on the CUSUM analysis, NDVI change can be identified in 39 of the 82 ROIs (48%), and the change mostly occurred between 2008 and 2014 (Figure 6).
NDVI trends driven by agricultural processes were identified in 62 ROIs, of which 35 were associated with an increase in NDVI values and 27 exhibited a decrease in NDVI values.NDVI values within agricultural areas (plantation and recovery, forest-to-terrace and forest-cutting) were often associated with relatively high values of mean, standard deviation and range (Figure 5).The oscillations related to agricultural practices belong to three major subclasses: (i) Conversion of a shrubland or woodland into a terrace or fruit forest or tea yard, resulting in a decrease in NDVI values and, consequently, low standard deviation or range of NDVI (Figure 5).(ii) The recovery and planting of shrubbery and forest resulted in an increase in NDVI values.
A negative correlation between standard deviation and the M-K τ value can be found.(iii) The cutting of forests resulted in decreasing NDVI values, as shown in Figure 7.The ROIs of forest-cutting have a scattered distribution.However, due to immediate reforestation, all decreases recovered in the next few years, and high mean NDVI values in ROIs of forest-cutting can be found.The second most common factor driving land cover change was construction (numbers of ROIs is 13), with all 13 ROIs exhibiting a sudden decrease in NDVI down to a very low level.The correlation between annual rainfall and NDVI is low in all ROIs; more specifically, the plantation and recovery class has positive values and most of the forest-cutting, construction and forest-to-terrace classes show negative values.There are two major types of trends identified in construction areas: a sudden decrease in NDVI that cannot be recovered to the original state and, thus, results in low mean NDVI values (Figure 7).Generally, precipitation is the limiting factor for plant growth in arid and semi-arid areas across the globe [38,45], while temperature is the limiting factor in the high latitudes of the northern hemisphere and Qinghai-Tibet Plateau [46,47].Thus, the impact of temperature on NDVI change in the lowlatitude study area was not taken into consideration.
Furthermore, as shown in Figure 8, different classes show different characteristics in the relationship between average NDVI and its corresponding CV.Specifically, most of the ROIs in the forest-cutting class exhibit low CV and relatively high NDVI due to rapid reforestation while woodland exhibits high NDVI.ROIs in the plantation and recovery class have a large range in average NDVI value because the change in ROIs in this class was caused by plantations on barren land or by recovery in closed forests; as a result, average NDVI values can be high or low, respectively.Figure 7 shows the impact of construction on NDVI, which can explain the fact that ROIs in the construction class have a large range in CV.In the construction class, there is a sudden decrease in NDVI values that cannot recover to the original state, resulting in low mean NDVI values (Figure 7).

Driving Forces of Vegetation Cover Change in the Study Area
It is convincing that the general ecological environment reflected by NDVI trend is continuously improving during the past 16 years.According to the NDVI change derived by the M-K τ method in different land-use types (Figure 9), it was found that during the study period, areas with M-K τ > 0.2 account for the largest amount (>61.2%) of the area in woodlands, a certain area in shrubland and cultivated land, and a small proportion (<2%) of construction land, garden plots and grassland.This also reflects that in all types of land use, except construction land, the area with M-K τ > 0.2 is much larger than the area with M-K τ < −0.2.As for the ratio of area with M-K τ > 0.2 and its distribution in different land-use types, 18% of shrubland underwent an increase in NDVI, compared with 10% for There are two major types of trends identified in construction areas: a sudden decrease in NDVI that cannot be recovered to the original state and, thus, results in low mean NDVI values (Figure 7).Generally, precipitation is the limiting factor for plant growth in arid and semi-arid areas across the globe [38,45], while temperature is the limiting factor in the high latitudes of the northern hemisphere and Qinghai-Tibet Plateau [46,47].Thus, the impact of temperature on NDVI change in the low-latitude study area was not taken into consideration.
Furthermore, as shown in Figure 8, different classes show different characteristics in the relationship between average NDVI and its corresponding CV.Specifically, most of the ROIs in the forest-cutting class exhibit low CV and relatively high NDVI due to rapid reforestation while woodland exhibits high NDVI.ROIs in the plantation and recovery class have a large range in average NDVI value because the change in ROIs in this class was caused by plantations on barren land or by recovery in closed forests; as a result, average NDVI values can be high or low, respectively.Figure 7 shows the impact of construction on NDVI, which can explain the fact that ROIs in the construction class have a large range in CV.In the construction class, there is a sudden decrease in NDVI values that cannot recover to the original state, resulting in low mean NDVI values (Figure 7).There are two major types of trends identified in construction areas: a sudden decrease in NDVI that cannot be recovered to the original state and, thus, results in low mean NDVI values (Figure 7).Generally, precipitation is the limiting factor for plant growth in arid and semi-arid areas across the globe [38,45], while temperature is the limiting factor in the high latitudes of the northern hemisphere and Qinghai-Tibet Plateau [46,47].Thus, the impact of temperature on NDVI change in the lowlatitude study area was not taken into consideration.
Furthermore, as shown in Figure 8, different classes show different characteristics in the relationship between average NDVI and its corresponding CV.Specifically, most of the ROIs in the forest-cutting class exhibit low CV and relatively high NDVI due to rapid reforestation while woodland exhibits high NDVI.ROIs in the plantation and recovery class have a large range in average NDVI value because the change in ROIs in this class was caused by plantations on barren land or by recovery in closed forests; as a result, average NDVI values can be high or low, respectively.Figure 7 shows the impact of construction on NDVI, which can explain the fact that ROIs in the construction class have a large range in CV.In the construction class, there is a sudden decrease in NDVI values that cannot recover to the original state, resulting in low mean NDVI values (Figure 7).

Driving Forces of Vegetation Cover Change in the Study Area
It is convincing that the general ecological environment reflected by NDVI trend is continuously improving during the past 16 years.According to the NDVI change derived by the M-K τ method in different land-use types (Figure 9), it was found that during the study period, areas with M-K τ > 0.2 account for the largest amount (>61.2%) of the area in woodlands, a certain area in shrubland and cultivated land, and a small proportion (<2%) of construction land, garden plots and grassland.This also reflects that in all types of land use, except construction land, the area with M-K τ > 0.2 is much larger than the area with M-K τ < −0.2.As for the ratio of area with M-K τ > 0.2 and its distribution in different land-use types, 18% of shrubland underwent an increase in NDVI, compared with 10% for

Driving Forces of Vegetation Cover Change in the Study Area
It is convincing that the general ecological environment reflected by NDVI trend is continuously improving during the past 16 years.According to the NDVI change derived by the M-K τ method in different land-use types (Figure 9), it was found that during the study period, areas with M-K τ > 0.2 account for the largest amount (>61.2%) of the area in woodlands, a certain area in shrubland and cultivated land, and a small proportion (<2%) of construction land, garden plots and grassland.This also reflects that in all types of land use, except construction land, the area with M-K τ > 0.2 is much larger than the area with M-K τ < −0.2.As for the ratio of area with M-K τ > 0.2 and its distribution in different land-use types, 18% of shrubland underwent an increase in NDVI, compared with 10% for gardens, 0.45%-0.75%for grasslands, cultivated land, woodland and construction land.However, for the ratio of area with M-K τ < −0.2 and its distribution in different land-use types, only construction land has a proportion close to 5%; the ratio of the remaining land-use types is <2.5%.As NDVI saturation will occur with a high canopy cover, and the main land cover type in the study area is forest with a high coverage rate [41], it stands to reason that when calculating the area of woodland with M-K τ > 0.2 or <−0.2, saturation may cause the calculated area to be slightly smaller than it should be.
Sustainability 2017, 9, 443 10 of 15 gardens, 0.45%-0.75%for grasslands, cultivated land, woodland and construction land.However, for the ratio of area with M-K τ < −0.2 and its distribution in different land-use types, only construction land has a proportion close to 5%; the ratio of the remaining land-use types is <2.5%.As NDVI saturation will occur with a high canopy cover, and the main land cover type in the study area is forest with a high coverage rate [41], it stands to reason that when calculating the area of woodland with M-K τ > 0.2 or <−0.2, saturation may cause the calculated area to be slightly smaller than it should be.By comparing the total area ratio of M-K τ > 0.2 and M-K τ < −0.2 (Figure 9), it is evident that in the study period, there were equal areas of construction land that showed increasing and decreasing Although the total amount of construction land with significant vegetation change is not large, construction is considered to be held responsible for NDVI changes in 13 of the 82 (16%) ROIs.The authors found the main drivers for a decrease in NDVI value in construction land to be urbanization and industrialization.The authors used two measures to respectively indicate urban development and industrialization: (1) the area of residential buildings under construction and (2) the numbers of industrial enterprises with annual revenue of 20 million ¥.It is shown in Figure 10 that there was a noticeable increase in these two indicators since 2007.Development of industrialization and urbanization mean more infrastructure, civilian construction and industrial parks, which further leads to sharp decreases in vegetation cover, after which vegetation cover increases very slowly (Figure 7).As the largest developing country in the world, China has undergone significant urbanization since the 1970s [48].Consequently, rapid economic growth rates and urban expansion have led to widespread land use and land cover (LULC) changes [49].In the previous studies investigating LULC change in humid mountainous regions, such as northwestern Yunnan of China, land-use change was also identified to be responsible for vegetation coverage change, which was reflected by vegetation coverage grades [7].For the areas with M-K τ > 0.2, the increase in vegetation change can be attributed to government-sponsored greening in urbanized areas.Urban green coverage was 35% until 2015 in Nanxiong City [50] and according to the Shaoguan Statistical Yearbook of 2015, the "Ten-thousand-green program" has resulted in a total of 0.2 million trees being planted by organized voluntary tree-planting events [50].By comparing the total area ratio of M-K τ > 0.2 and M-K τ < −0.2 (Figure 9), it is evident that in the study period, there were equal areas of construction land that showed increasing and decreasing trends.Although the total amount of construction land with significant vegetation change is not large, construction is considered to be held responsible for NDVI changes in 13 of the 82 (16%) ROIs.The authors found the main drivers for a decrease in NDVI value in construction land to be urbanization and industrialization.The authors used two measures to respectively indicate urban development and industrialization: (1) the area of residential buildings under construction and (2) the numbers of industrial enterprises with annual revenue of 20 million ¥.It is shown in Figure 10 that there was a noticeable increase in these two indicators since 2007.Development of industrialization and urbanization mean more infrastructure, civilian construction and industrial parks, which further leads to sharp decreases in vegetation cover, after which vegetation cover increases very slowly (Figure 7).As the largest developing country in the world, China has undergone significant urbanization since the 1970s [48].Consequently, rapid economic growth rates and urban expansion have led to widespread land use and land cover (LULC) changes [49].In the previous studies investigating LULC change in humid mountainous regions, such as northwestern Yunnan of China, land-use change was also identified to be responsible for vegetation coverage change, which was reflected by vegetation coverage grades [7].For the areas with M-K τ > 0.2, the increase in vegetation change can be attributed to government-sponsored greening in urbanized areas.Urban green coverage was 35% until 2015 in Nanxiong City  As shown in Figure 2, a total of 35 ROIs of plantation and recovery accounts for 43% of all selected ROIs.Meanwhile in Figure 9, almost 20% of woodland shows an increase in vegetation.This significant increase is also driven by government-sponsored afforestation.Figure 11a shows that accumulative afforestation covers 139.8 km 2 (not including reforestation areas during the study period) in 2000-2013 (data from 2003 and 2010 were missing), which accounts for 2.6% of the total woodland area.Furthermore, the average area and maximum area fenced off for afforestation were 861 km 2 and 1899 km 2 , which accounts for 15.2% and 33.3% of the woodland area, respectively.The increased area fenced off for afforestation was also responsible for the increase in vegetation cover change.In 1999, the Chinese government initiated the Natural Forest Resources Conservation Program and the Restoring Farmland into Forest Program.Logging has been prohibited in most natural forests, and cultivated land on areas with slopes steeper than 25° must be restored to forests or grasslands [51].Consequently, the increase in vegetation cover can be attributed to the plantation and afforestation policy.Figure 2d shows that 21 of the 82 ROIs fall into the forest-cutting class, suggesting that forestcutting is a significant factor in vegetation dynamics change.Due to rapid reforestation each year and the relative small ratio between average deforestation volume (177,252 m 3 ) and average growth volume (725,683 m 3 ), the NDVI in timber harvest ROIs recovered rapidly (Figure 12).In addition, it As shown in Figure 2, a total of 35 ROIs of plantation and recovery accounts for 43% of all selected ROIs.Meanwhile in Figure 9, almost 20% of woodland shows an increase in vegetation.This significant increase is also driven by government-sponsored afforestation.Figure 11a shows that accumulative afforestation covers 139.8 km 2 (not including reforestation areas during the study period) in 2000-2013 (data from 2003 and 2010 were missing), which accounts for 2.6% of the total woodland area.Furthermore, the average area and maximum area fenced off for afforestation were 861 km 2 and 1899 km 2 , which accounts for 15.2% and 33.3% of the woodland area, respectively.The increased area fenced off for afforestation was also responsible for the increase in vegetation cover change.In 1999, the Chinese government initiated the Natural Forest Resources Conservation Program and the Restoring Farmland into Forest Program.Logging has been prohibited in most natural forests, and cultivated land on areas with slopes steeper than 25 • must be restored to forests or grasslands [51].Consequently, the increase in vegetation cover can be attributed to the plantation and afforestation policy.As shown in Figure 2, a total of 35 ROIs of plantation and recovery accounts for 43% of all selected ROIs.Meanwhile in Figure 9, almost 20% of woodland shows an increase in vegetation.This significant increase is also driven by government-sponsored afforestation.Figure 11a shows that accumulative afforestation covers 139.8 km 2 (not including reforestation areas during the study period) in 2000-2013 (data from 2003 and 2010 were missing), which accounts for 2.6% of the total woodland area.Furthermore, the average area and maximum area fenced off for afforestation were 861 km 2 and 1899 km 2 , which accounts for 15.2% and 33.3% of the woodland area, respectively.The increased area fenced off for afforestation was also responsible for the increase in vegetation cover change.In 1999, the Chinese government initiated the Natural Forest Resources Conservation Program and the Restoring Farmland into Forest Program.Logging has been prohibited in most natural forests, and cultivated land on areas with slopes steeper than 25° must be restored to forests or grasslands [51].Consequently, the increase in vegetation cover can be attributed to the plantation and afforestation policy.Figure 2d shows that 21 of the 82 ROIs fall into the forest-cutting class, suggesting that forestcutting is a significant factor in vegetation dynamics change.Due to rapid reforestation each year and the relative small ratio between average deforestation volume (177,252 m 3 ) and average growth volume (725,683 m 3 ), the NDVI in timber harvest ROIs recovered rapidly (Figure 12).In addition, it Figure 2d shows that 21 of the 82 ROIs fall into the forest-cutting class, suggesting that forestcutting is a significant factor in vegetation dynamics change.Due to rapid reforestation each year and the relative small ratio between average deforestation volume (177,252 m 3 ) and average growth volume (725,683 m 3 ), the NDVI in timber harvest ROIs recovered rapidly (Figure 12).In addition, it has been demonstrated that controlled deforestation and rapid plant replacement is environmentally acceptable in humid regions such as Malaysia [53].
Sustainability 2017, 9, 443 12 of 15 has been demonstrated that controlled deforestation and rapid plant replacement is environmentally acceptable in humid regions such as Malaysia [53].Vegetation variation in cultivated areas is the direct result of agricultural tillage.Previous studies have found that vegetation on different land-use types will show different responses to climate change and human activities [54][55][56][57].In addition, natural environmental change and human activities are the main two driving forces associated with vegetation cover change.In this case, four few-noticed but important driving forces, i.e., plantation, construction, forest-cutting, forest-toterrace, are representatives of human impact.By analyzing the interactions between forest above ground biomass (AGB) and climate data as well as forest disturbance in northern Guangdong of China, Shen et al. [58] found that human activities make a greater contribution to short-term AGB dynamics than climate change.
Globally, rainfall was identified as a major driver for NDVI variability, explaining more than 40% of variability in temperate to tropical water-limited herbaceous systems, whereas, temperature was a major driver in northern mid-to high-latitudes in the spring and early summer, explaining up to 70% of NDVI variance [59].However, in contrast to the aforementioned global studies, climate variability was found to be a relatively minor driver of changes in NDVI values in the study area at Nanxiong Basin.

Conclusions
Taking Nanxiong Basin in Guangdong Province of China as the case study area, vegetation coverage change and associated driving forces in mountainous areas were analyzed with the application of GIS techniques and statistical methods.
The results showed that vegetation coverage (indicated by NDVI derived from the MODIS remote sensing dataset) for the whole basin conveyed an overall increasing trend during 2000-2015, and the significant change in vegetation coverage was characterized by a gradual increase in all types of land-use (except for construction land), which results from environmental policy implementation.The general ecological environment indicated by NDVI has been improved during the past 16 years.
For construction land, the increase of vegetation coverage was mainly caused by governmentsponsored greening in urbanized areas, while the decrease of vegetation in the same land cover was caused by urbanization and industrialization.The total area of woodland with an increase in vegetation cover is much higher than that with a decrease in vegetation cover, and the former was mainly driven by afforestation and accounted for by afforestation policy.In conclusion, most of the dynamics in NDVI during the study period can be attributed to anthropogenic effects in Nanxiong Basin.Vegetation variation in cultivated areas is the direct result of agricultural tillage.Previous studies have found that vegetation on different land-use types will show different responses to climate change and human activities [54][55][56][57].In addition, natural environmental change and human activities are the main two driving forces associated with vegetation cover change.In this case, four few-noticed but important driving forces, i.e., plantation, construction, forest-cutting, forest-to-terrace, are representatives of human impact.By analyzing the interactions between forest above ground biomass (AGB) and climate data as well as forest disturbance in northern Guangdong of China, Shen et al. [58] found that human activities make a greater contribution to short-term AGB dynamics than climate change.
Globally, rainfall was identified as a major driver for NDVI variability, explaining more than 40% of variability in temperate to tropical water-limited herbaceous systems, whereas, temperature was a major driver in northern mid-to high-latitudes in the spring and early summer, explaining up to 70% of NDVI variance [59].However, in contrast to the aforementioned global studies, climate variability was found to be a relatively minor driver of changes in NDVI values in the study area at Nanxiong Basin.

Conclusions
Taking Nanxiong Basin in Guangdong Province of China as the case study area, vegetation coverage change and associated driving forces in mountainous areas were analyzed with the application of GIS techniques and statistical methods.
The results showed that vegetation coverage (indicated by NDVI derived from the MODIS remote sensing dataset) for the whole basin conveyed an overall increasing trend during 2000-2015, and the significant change in vegetation coverage was characterized by a gradual increase in all types of land-use (except for construction land), which results from environmental policy implementation.The general ecological environment indicated by NDVI has been improved during the past 16 years.
For construction land, the increase of vegetation coverage was mainly caused by governmentsponsored greening in urbanized areas, while the decrease of vegetation in the same land cover was caused by urbanization and industrialization.The total area of woodland with an increase in vegetation cover is much higher than that with a decrease in vegetation cover, and the former was mainly driven by afforestation and accounted for by afforestation policy.In conclusion, most of the dynamics in NDVI during the study period can be attributed to anthropogenic effects in Nanxiong Basin.

Figure 1 .
Figure 1.Elevation of Nanxiong Basin; the subset figure in the upper right and upper left show the location of the study area in Guangdong Province and in China.

Figure 1 .
Figure 1.Elevation of Nanxiong Basin; the subset figure in the upper right and upper left show the location of the study area in Guangdong Province and in China.

Figure 1 .
Figure 1.Elevation of Nanxiong Basin; the subset figure in the upper right and upper left show the location of the study area in Guangdong Province and in China.

Figure 2 .
Figure 2. Nonparametric Mann-Kendall τ (M-K τ) correlation coefficient calculated from the NDVI time series covering the study period (only statistically significant trends are shown) (a); land-use map of Nanxiong Basin in 2010 (b); the coefficient of variation (CV) in the NDVI time series and (c); the mean NDVI over the study period and spatial distribution of regions of interest (ROIs) (d).

Figure 2 .
Figure 2. Nonparametric Mann-Kendall τ (M-K τ) correlation coefficient calculated from the NDVI time series covering the study period (only statistically significant trends are shown) (a); land-use map of Nanxiong Basin in 2010 (b); the coefficient of variation (CV) in the NDVI time series and (c); the mean NDVI over the study period and spatial distribution of regions of interest (ROIs) (d).

Figure 3 .
Figure 3.The relationship between average annual mean NDVI and CV of annual mean NDVI in six main types of land use: (a) construction land; (b) cultivated land; (c) garden plots; (d) grassland; (e) shrubland; and (f) woodland.

Figure 3 .
Figure 3.The relationship between average annual mean NDVI and CV of annual mean NDVI in six main types of land use: (a) construction land; (b) cultivated land; (c) garden plots; (d) grassland; (e) shrubland; and (f) woodland.
of them (covering 14.25 km 2 ) represented areas with positive trends in NDVI, and 47 ROIs (covering 39.19 km 2 ) represented areas with a negative trend in NDVI.A global Moran's I spatial autocorrelation test on the M-K τ values of the 82 ROIs found a nearest neighbor ratio of −0.11 (z = −1.69,p = 0.09), indicating a clustered spatial pattern of observed NDVI trends in the 82 ROIs.
of them (covering 14.25 km 2 ) represented areas with positive trends in NDVI, and 47 ROIs (covering 39.19 km 2 ) represented areas with a negative trend in NDVI.A global Moran's I spatial autocorrelation test on the M-K τ values of the 82 ROIs found a nearest neighbor ratio of −0.11 (z = −1.69,p = 0.09), indicating a clustered spatial pattern of observed NDVI trends in the 82 ROIs.

Figure 4 .
Figure 4. Correlation coefficients for the 82 regions of interest showing annual rainfall NDVI values as a function of the correlation coefficients between cumulative rainfall (over four months) and mean monthly NDVI values.

Figure 4 .
Figure 4. Correlation coefficients for the 82 regions of interest showing annual rainfall NDVI values as a function of the correlation coefficients between cumulative rainfall (over four months) and mean monthly NDVI values.

Figure 5 .
Figure 5. M-K τ correlation coefficient of the 82 ROIs as a function of annual rainfall (a); standard deviation of NDVI values (b); range of NDVI values (c); and mean NDVI values (d).All of the trends in the 82 regions of interest were statistically significant.

Figure 6 .
Figure 6.Number of ROIs for the years when change in NDVI values took place.

Figure 5 . 15 Figure 5 .
Figure 5. M-K τ correlation coefficient of the 82 ROIs as a function of annual rainfall (a); standard deviation of NDVI values (b); range of NDVI values (c); and mean NDVI values (d).All of the trends in the 82 regions of interest were statistically significant.

Figure 6 .
Figure 6.Number of ROIs for the years when change in NDVI values took place.

Figure 6 .
Figure 6.Number of ROIs for the years when change in NDVI values took place.

Figure 7 .
Figure 7.Typical trends of NDVI values in four types of ROIs: (a) plantation and recovery; (b) construction; (c) forest-cutting; and (d) forest-to-terrace.

Figure 8 .
Figure 8.The relationship between the average value and CV of annual mean NDVI.

Figure 7 .
Figure 7.Typical trends of NDVI values in four types of ROIs: (a) plantation and recovery; (b) construction; (c) forest-cutting; and (d) forest-to-terrace.

Figure 7 .
Figure 7.Typical trends of NDVI values in four types of ROIs: (a) plantation and recovery; (b) construction; (c) forest-cutting; and (d) forest-to-terrace.

Figure 8 .
Figure 8.The relationship between the average value and CV of annual mean NDVI.

Figure 8 .
Figure 8.The relationship between the average value and CV of annual mean NDVI.

Figure 9 .
Figure 9.The areas with M-K τ > 0.2 and <−0.2 in each land-use type (a); and the ratio of areas with M-K τ > 0.2 and <−0.2 and total area in each land-use type (b).

Figure 9 .
Figure 9.The areas with M-K τ > 0.2 and <−0.2 in each land-use type (a); and the ratio of areas with M-K τ > 0.2 and <−0.2 and total area in each land-use type (b).
[50]  and according to the Shaoguan Statistical Yearbook of 2015, the "Ten-thousand-green program" has resulted in a total of 0.2 million trees being planted by organized voluntary tree-planting events [50].

Figure 10 .
Figure 10.The areas of residential buildings under construction each year and the number of industrial enterprises with an annual revenue of 20 million ¥.Data were obtained from the Shaoguan Statistical Yearbook [50].

Figure 11 .
Figure 11.Annual plantation area and cumulative plantation area, excluding the reforestation area each year (a), data were obtained from Shaoguan Statistical Yearbook [50]; The area of closed forest for recovery each year (b), data were collected from Nanxiong Forestry Bureau [52].

Figure 10 .
Figure 10.The areas of residential buildings under construction each year and the number of industrial enterprises with an annual revenue of 20 million ¥.Data were obtained from the Shaoguan Statistical Yearbook [50].

Figure 10 .
Figure 10.The areas of residential buildings under construction each year and the number of industrial enterprises with an annual revenue of 20 million ¥.Data were obtained from the Shaoguan Statistical Yearbook [50].

Figure 11 .
Figure 11.Annual plantation area and cumulative plantation area, excluding the reforestation area each year (a), data were obtained from Shaoguan Statistical Yearbook [50]; The area of closed forest for recovery each year (b), data were collected from Nanxiong Forestry Bureau [52].

Figure 11 .
Figure 11.Annual plantation area and cumulative plantation area, excluding the reforestation area each year (a), data were obtained from Shaoguan Statistical Yearbook [50]; The area of closed forest for recovery each year (b), data were collected from Nanxiong Forestry Bureau [52].

Figure 12 .
Figure 12.Annual forest-cutting volume and reforestation area during 2000-2013 (a), data were obtained from the Shaoguan Statistical Yearbook [50]; The annual growth of tree volume in 2003-2013 (b), data were collected from Nanxiong Forestry Bureau [52].

Figure 12 .
Figure 12.Annual forest-cutting volume and reforestation area during 2000-2013 (a), data were obtained from the Shaoguan Statistical Yearbook [50]; The annual growth of tree volume in 2003-2013 (b), data were collected from Nanxiong Forestry Bureau [52].
• C; the annual mean precipitation and potential evaporation are 1555.1 mm and 1678.7 mm, respectively [29].The rainy season is from March to August.According to the 2009 1:5000 land-use map from the Shaoguan Municipal Bureau of Land and Resources, land-use types (Figure