**Spatial–Temporal Evolution of Vegetation NDVI in Association with Climatic, Environmental and Anthropogenic Factors in the Loess Plateau, China during 2000–2015: Quantitative Analysis Based on Geographical Detector Model**

**Yi Dong 1, Dongqin Yin 1,2,3,\*, Xiang Li 4,5,6, Jianxi Huang 1,2, Wei Su 1,2, Xuecao Li 1,2 and Hongshuo Wang 1,2**


**Abstract:** In the Loess Plateau (LP) of China, the vegetation degradation and soil erosion problems have been shown to be curbed after the implementation of the Grain for Green program. In this study, the LP is divided into the northwestern semi-arid area and the southeastern semi-humid area using the 400 mm isohyet. The spatial–temporal evolution of the vegetation NDVI during 2000–2015 are analyzed, and the driving forces (including factors of climate, environment, and human activities) of the evolution are quantitatively identified using the geographical detector model (GDM). The results showed that the annual mean NDVI in the entire LP was 0.529, and it decreased from the semi-humid area (0.619) to the semi-arid area (0.346). The mean value of the coefficient of variation of the NDVI was 0.1406, and it increased from the semi-humid area (0.1165) to the semi-arid area (0.1926). The annual NDVI growth rate in the entire LP was 0.0079, with the NDVI growing faster in the semi-humid area (0.0093) than in the semi-arid area (0.0049). The largest increments of the NDVI were from grassland, farmland, and woodland. The GDM results revealed that changes in the spatial distribution of the NDVI could be primarily explained by the climatic and environmental factors in the semi-arid area, such as precipitation, soil type, and vegetation type, while the changes were mainly explained by the anthropogenic factors in the semi-humid area, such as the GDP density, land-use type, and population density. The interactive analysis showed that interactions between factors strengthened the impacts on the vegetation change compared with an individual factor. Furthermore, the ranges/types of factors suitable for vegetation growth were determined. The conclusions of this study have important implications for the formulation and implementation of ecological conservation and restoration strategies in different regions of the LP.

**Keywords:** Loess Plateau; China; normalized difference vegetation index (NDVI); spatial–temporal evolution; geographical detector model; driving forces

#### **1. Introduction**

Vegetation plays an indispensable role in regional terrestrial ecosystems, and constitutes an essential link between soil, water, and atmosphere. Moreover, vegetation is the material basis for the survival of surface organisms [1,2]. Vegetation coverage effectively

**Citation:** Dong, Y.; Yin, D.; Li, X.; Huang, J.; Su, W.; Li, X.; Wang, H. Spatial–Temporal Evolution of Vegetation NDVI in Association with Climatic, Environmental and Anthropogenic Factors in the Loess Plateau, China during 2000–2015: Quantitative Analysis Based on Geographical Detector Model. *Remote Sens.* **2021**, *13*, 4380. https://doi.org/ 10.3390/rs13214380

Academic Editors: Xuejia Wang, Tinghai Ou, Wenxin Zhang and Youhua Ran

Received: 7 October 2021 Accepted: 28 October 2021 Published: 30 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

reduces the surface soil erosion caused by exogenous forces such as wind, diminishes the splash erosion caused by raindrops, alleviates the hydrodynamic erosion of rivers, and improves the soil environment. Therefore, it is of crucial importance to explore the vegetation coverage changes and dynamics for the soil erosion prevention and control, ecological environmental protection, and sustainable social and economic development [3,4]

At present, the main monitoring method used in large-scale vegetation coverage change research is based on satellite remote sensing because of its wide spatial range and gradually improving resolution, which effectively makes up for the shortcomings of traditional monitoring methods [5]. The normalized difference vegetation index (NDVI) is strongly associated with the vegetation coverage, leaf area index, biomass, and land use, which can comprehensively reflect the vegetation's greenness, photosynthesis intensity, and vegetation metabolism intensity [6,7]. The NDVI can be used to quantitatively evaluate the regional vegetation coverage and growth, which is considered to be an effective indicator for monitoring terrestrial vegetation changes, and thus, has been widely used in research and management in various fields, such as agriculture and ecology [8].

The growth process of vegetation is affected by multiple factors. Temperature and precipitation are directly related to global climate change and are generally regarded as the most important natural factors affecting vegetation growth and changes in the longterm [9,10]. Since the 20th century, intense human activities, which are characterized by industrialization and urbanization, have brought about rapid population and economic growth, rapid changes in land use, and ecological and environmental problems such as vegetation degradation and soil erosion, resulting in significant impacts on vegetation growth and changes in the short-term [11–13].

The Loess Plateau (LP) of China is located in the semi-arid/humid zones, the local ecological environment of which is extremely fragile. Vegetation degradation and soil erosion have been particularly prominent in this region throughout history, making the LP main sediment source (nearly 90%) of the Yellow River (middle and lower reaches), once the most sand-laden river in the world [14]. Exploring the evolution of vegetation and its driving forces on the LP cannot only help formulate strategies on ecological restoration for this area, but also help predict the sediment and tackle the sediment related problems (such as reservoir sedimentation and flood control) in the Yellow River. Since the 1980s, particularly after 2000 when the Grain for Green program (GGP) was implemented in this area, the restoration of vegetation has been witnessed, and accordingly, the expansion trend of the soil erosion has been curbed [15]. In addition, a new stage of the GGP was launched in 2014 to consolidate the achievements [16].

Many scholars have studied the causes of the vegetation coverage changes on the LP. There are limitations in the existing research, which mainly include the following. (1) Many studies used NDVI data from original Global Inventory Modeling and Mapping Studies (GIMMS) or GIMMS NDVI data expanded from Moderate Resolution Imaging Spectroradiometer (MODIS) and Satellite Pour l'Observation de la Terre (SPOT) Vegetation (VGT) NDVI [6,17–19]. The spatial resolution of 8 km leads to serious mixed pixels and insufficient detail capture, making it difficult for the results to show the actual spatial– temporal changes in the vegetation coverage. (2) Several studies used statistical methods, such as correlation or regression, by assuming that the vegetation growth is linearly related to the potential factors with time, but they ignored the spatial heterogeneity [20,21]. Some evidence has suggested that there is no strict linear correlation between the factors and the vegetation growth under the climate change. A statistical linear correlation model may not be able to accurately describe the internal relationship [22,23]. (3) Several studies only considered the influences of climatic factors on vegetation coverage changes, and they neglected other environmental factors. Actually, certain environmental factors (e.g., altitude, slope, and slope aspect) also have key impacts on vegetation growth and changes. For instance, the altitude affects the temperature, precipitation, and soil, thus affecting the vertical distribution of the vegetation. The slope aspect affects the degree to which the slope receives various hydrothermal conditions, which has a certain impact on the type

and distribution of the vegetation. Moreover, the steep slopes mean that there may be serious soil erosion, which inhibits the vegetation growth. However, many environmental factors cannot be quantified using statistical methods, which limits the exploration of the effects of environmental factors on vegetation changes [24,25]. (4) Most previous studies merely evaluated the individual effects of the factors without quantitatively evaluating the interactive effects between multiple factors on vegetation changes. (5) The majority of previous studies focused on the entire area and local area of the LP (e.g., relevant research was carried out on Shannxi Province in the LP [26]), little comparative research has been conducted for different climatic regions, and the research conclusions are not conducive to guiding policymaking for different regions.

In this study, attempts were made to address the above research gaps. SPOT NDVI products were used to evaluate the vegetation coverage and growth. The spatial–temporal evolution of the NDVI on the LP during 2000–2015 was explained by the Geographical Detector Model (GDM) [27]. Based on the spatial analysis of variance, the GDM was shown to avoid linear assumption between variables and has a clear physical meaning, which reflects the explanatory power of an individual independent variable acting alone or of multiple independent variables interacting on a dependent variable. Moreover, the GDM can utilize a variety of type variables such as geomorphic type, soil type, vegetation type, and land-use type, and thus, is superior to the traditional statistical methods. The GDM is widely adopted in the exploration of the spatial differentiation characteristics of issues in natural and social sciences. For instance, the GDM was applied to study the explanatory powers of different factors on vegetation changes in Shannxi, Mongolia, Sichuan, and the Heihe River Basin [26,28–30].

The study begins by introducing the study area, method and data in Section 2. Then, the linear regression, coefficient of variation and transfer matrix models are used to identify the spatial–temporal evolution of the NDVI on the semi-arid/humid areas of the LP during 2000–2015 in Section 3. The influences of climatic (precipitation and temperature), environmental (altitude, slope, slope aspect, geomorphic type, soil type, and vegetation type), and anthropogenic factors (GDP density, population density, land-use type) on the spatial–temporal evolution of the vegetation NDVI on the semi-arid/humid areas of the LP are identified with the factor, interaction, risk and ecological detectors of the GDM. Section 4 further compares the results in the semi-arid and sub-humid areas, proposes suggestions for formulating and implementing ecological protection and restoration strategies in different areas, and discusses the connections and distinctions with other studies as well as the possible future work. Section 5 summarizes the main findings of the paper. The achievements gained from this study will benefit policy makers and administrative managers in their strategy development and implementation.

#### **2. Data and Methods**

#### *2.1. Study Area*

The LP is located north of central China, with a total area of 6.49 × 105 km<sup>2</sup> (107◦54 – 114◦33 E, 33◦41 –41◦16 N), as shown in Figure 1. It is the transitional zone from the coastal humid region to the inland arid region, and from forest to grassland. The terrain is generally low in the southeast and high in the northwest, including the mountain area, the loess hilly area, the loess tableland area, and the valley plain area. The LP has a temperate monsoon climate, with an annual mean temperature of 3.6–14.3 ◦C and precipitation of 200–800 mm. Both the temperature and precipitation increase from the northwest to the southeast. The annual and daily temperature ranges are large. Furthermore, the seasonal variation of precipitation is large. Heavy rains occur frequently during summer, resulting in the soil erosion, floods, landslides, debris flows and other disasters. The total annual radiation is 50.2 × 104 J/cm<sup>2</sup> to 67.0 × <sup>10</sup><sup>4</sup> J/cm2, with a long illumination time and high radiation. About 200 rivers have their headwaters on the LP, including the Tao River, Zuli River, Qingshui River, Huangfuchuan River, Kuye River, Wuding River, Beiluo River Wei River, Qin River, Fen River, etc. The study area includes the provinces (or autonomous regions) of Qinghai, Gansu, Ningxia, Inner Mongolia, Shaanxi, Shanxi, and Henan. Throughout history, the overall vegetation coverage has been low and soil erosion has been very serious on the LP due to the special natural and geographical environment, complex vegetation types, and severe impacts of human activities. Since the 1980s, a series of measures for soil and water conservation have been taken to control the soil erosion. Particularly, after 2000, when the GGP was implemented, the regional ecological environment has significantly improved.

**Figure 1.** Study area.

The study area was divided into semi-arid and semi-humid areas using the 400 mm isohyet (calculated using the annual mean precipitation data for the LP during 2000–2015) to carry out the comparative analysis. It should be noted that the 400 mm isohyet is of great significance to the physical geography and socio-economic regionalization in China, and it coincides with the Huhuanyong Line. In addition to dividing the climatic zones, this line divides the forest vegetation from the grassland vegetation and the densely populated area from the sparsely populated area [31,32].

#### *2.2. Geographical Detector Model (GDM)*

The GDM is proposed based on spatial differentiation theory and geographic information system (GIS) spatial analysis technique [27]. It is usually employed to study the factors affecting spatial hierarchical heterogeneities and the underlying mechanisms. The model assumes similar spatial pattern between the independent and dependent variables if these variables are spatially correlated [33]. The GDM consists of four detectors, including the factor, interaction, risk and ecological detectors.

(1) The factor detector can be utilized to analyze the spatial heterogeneity and to quantify the explanatory power (measured by the *q* value) of different impact factors *X* to dependent variable *Y*.

$$q = 1 - \frac{\sum\_{h=1}^{L} N\_h \sigma\_h^2}{N \sigma^2} \tag{1}$$

where *q* ∈ [0, 1], and with the increase in the *q* value, the explanatory power is expected to be stronger; *h* represents the stratum of variable *Y* (or factor *X*), *h* ∈ [1, L]; *Nh* is the unit number of stratum *h*; *N* is the unit number in all the strata; *σ*<sup>2</sup> *<sup>h</sup>* represents the variance of variable *Y* of stratum *h*; *σ*<sup>2</sup> is the variance of variable *Y* of all the strata.

(2) The interaction detector can be used to determine the explanatory power of interaction between two factors (say *X*<sup>1</sup> and *X*2) on the spatial heterogeneity of variable *Y*, and to judge whether the interactive effect on variable *Y* would be enhanced or weakened. The steps are as follows. (1) Compute the respective *q* value of factors *X*<sup>1</sup> and *X*<sup>2</sup> to variable *Y* (*qX*<sup>1</sup> and *qX*2). (2) Compute the *q* value of the interaction between *X*<sup>1</sup> and *X*<sup>2</sup> to *Y* (*qX*1&*X*2). (3) Compare *qX*1, *qX*2, and *qX*1&*X*2. If *max*(*qX*1,*qX*2) < *qX*1&*X*<sup>2</sup> < *qX*<sup>1</sup> + *qX*<sup>2</sup> is true, the factors *X*<sup>1</sup> and *X*<sup>2</sup> enhance each other (bi-enhance). If *qX*1&*X*<sup>2</sup> > *qX*<sup>1</sup> + *qX*<sup>2</sup> is true, the nonlinearity of two factors is enhanced (nonlinear enhancement). The interaction criteria are presented in Table 1.

**Table 1.** Definition of interaction detector.


(3) The risk detector is utilized to judge the difference of significance between the attribute mean values of two strata with the *t* statistic:

$$t\_{\overline{\mathcal{Y}}\_{h=1}} - y\_{h=2} = \frac{\overline{Y}\_{h=1} - \overline{Y}\_{h=2}}{\left[\frac{\text{Var}\left(\overline{\mathcal{T}}\_{h=1}\right)}{n\_{h=1}} + \frac{\text{Var}\left(\overline{\mathcal{T}}\_{h=2}\right)}{n\_{h=2}}\right]^{0.5}}\tag{2}$$

where *Yh* is the attribute mean value within stratum *h*; *nh* is the number of samples within stratum *h*; *Var* is the variance.

(4) The ecological detector is developed to compare the explanatory powers of two factors, *X*<sup>1</sup> and *X*2, on variable *Y* with the *F* statistic:

$$F = \frac{N\_{X1}(N\_{x2} - 1)SIV\_{X1}}{N\_{X2}(N\_{x1} - 1)SIV\_{X2}}\tag{3}$$

$$SIV\_{X1} = \sum\_{h=1}^{L1} N\_h \sigma\_{h'}^2\\ SIV\_{X2} = \sum\_{h=1}^{L2} N\_h \sigma\_h^2 \tag{4}$$

where *NX*<sup>1</sup> and *NX*<sup>2</sup> are the sample numbers of factors *X*<sup>1</sup> and *X*2, respectively; *SIVX*<sup>1</sup> and *SIVX*<sup>2</sup> are the sum of the internal variances of the strata from factors *X*<sup>1</sup> and *X*2, respectively; and *L*<sup>1</sup> and *L*<sup>2</sup> are the strata numbers of factors *X*<sup>1</sup> and *X*2, respectively.

#### *2.3. Data Description*

In this study, multi-source data (Table 2) were collected, mainly including the following.

(1) Vector data: The shape file data for the LP were downloaded from the Resource and Environment Science and Data Center (RESDC), Chinese Academy of Sciences (available from http://www.resdc.cn, accessed on 27 October 2021).

(2) Vegetation NDVI data: The SPOT NDVI data were used in this study (Figure 2), which were obtained from the RESDC. The time span is 2000–2015, the time step is 1 year, and the spatial resolution is 1 km × 1 km. It should be noted that the annual maximum synthesis method was used to obtain the annual value of NDVI, and therefore, there was no area with NDVI < 0.


**Table 2.** Climatic, environmental and anthropogenic factors.

**Figure 2.** Spatial distribution of NDVI on Loess Plateau during 2000–2015.

According to [34], the NDVI of the LP was divided into 5 intervals using the equal interval method, i.e., 0 to 0.2, 0.2 to 0.4, 0.4 to 0.6, 0.6 to 0.8, and 0.8 to 1.0, which represent 5 types of vegetation, i.e., low, relatively low, medium, relatively high, and high vegetation coverages, respectively.

(3) Meteorological data: The monthly mean precipitation and temperature from 101 meteorological stations on the LP and in its surrounding areas during 2000–2015 were collected for use, which were downloaded from the China Meteorological Data Service

Center (available from http://data.cma.cn/site/index.html, accessed on 27 October 2021). The annual value of the data from each meteorological station was calculated. Moreover, ANUSPLIN v4.3 was used to perform spatial interpolation of the meteorological station data, which were divided into 6 intervals using the equal interval method. Therefore, annual precipitation and temperature grid data (Figures 3 and 4) with the same projection mode, spatial resolution, and time series as the NDVI data were obtained. It should be noted that ANUSPLIN is a widely used software application for spatial interpolation, which was developed on the basis of the thin plate smoothing splines theory. The main feature of ANUSPLIN is that the topographic information can be used in the spatial interpolation of the meteorological data [35].

(4) DEM data: The DEM data were downloaded from the Shuttle Radar Topography Mission (available from http://srtm.csi.cgiar.org/, accessed on 27 October 2021). According to [28], the natural breakpoint method was used to divide the elevation data into 6 categories, including 90 to 790 m, 790 to 1228 m, 1228 to 1611 m, 1611 to 2136 m, 2136 to 2963 m, and 2963 to 4914 m (Figure 5a).

Based on the DEM data, the slope and slope aspect data for the study area were extracted using GIS spatial analysis tools. According to the policy of the Grain for Green [36], the slope was divided into 6 categories, including 0◦ to 5◦, 5◦ to 10◦, 10◦ to 15◦, 15◦ to 20◦, 20◦ to 25◦, >25◦ (Figure 5b).

The slope aspect was divided into 9 categories, which are denoted as no slope aspect (−1), east slope (67.5◦ to 112.5◦), west slope (247.5◦ to 292.5◦), south slope (157.5◦ to 202.5◦), north slope (0◦ to 22.5◦ and 337.5◦ to 360◦), southeast slope (112.5◦ to 157.5◦), northeast slope (22.5◦ to 67.5◦), southwest slope (202.5◦ to 247.5◦), and northwest slope (292.5◦ to 337.5◦) (Figure 5c).

**Figure 3.** Spatial distribution of precipitation on Loess Plateau during 2000–2015.

**Figure 4.** Spatial distribution of temperature on Loess Plateau during 2000–2015.

(5) Geomorphic type data: These data were obtained from the RESDC. Referring to the Geomorphological Atlas of the People's Republic of China (1:1 million) (2009), the geomorphology was divided into 6 types, including plain, platform, hill, small undulating mountain, medium undulating mountain, and large undulating mountain (Figure 5d).

(6) Soil type data: These data were downloaded from the RESDC. Referring to the 1:1 Million Soil Map of the People's Republic of China (1995), the soil was divided into 11 types, including alpine soil, anthropogenic soil, saline alkali soil, hydrogenetic soil, semi-hydrogenetic soil, primary soil, desert soil, arid soil, calcareous soil, semi-eluvial soil, and eluvial soil (Figure 5e).

(7) Vegetation type data: These data were obtained from the RESDC. Referring to the 1:1 Million Vegetation Atlas of China (2001), the vegetation was divided into 8 types, including cultivated vegetation, meadow, grass, grassland, desert, shrub, broad-leaved forest, and coniferous forest (Figure 5f).

(8) Socio-economic data: These data were obtained from the RESDC. The GDP density and population density data for China in 2000, 2005, 2010, and 2015 were obtained. The GDP density and population density were divided into 6 categories, respectively (Figures 6 and 7).

(9) Land-use type data: These data were downloaded from the RESDC. The spatial resolution is 1 km × 1 km. There are 6 land-use types, including unused land, construction land, water bodies, grassland, woodland, and farmland (Figure 8), with a total of 36 transfer combinations.

The influences of the various factors on the spatial–temporal evolution of the NDVI on the LP were calculated and analyzed with the GDM. The various factors were classified into climatic, environmental, and anthropogenic aspects. For the climatic factors, precipitation, and temperature were selected as proxy variables. For the environmental factors, altitude, slope, slope aspect, geomorphic type, soil type, and vegetation type were selected as proxy variables. For the anthropogenic factors, the proxy variables were the GDP density, population density, and land-use type.

**Figure 5.** Ranges/types of environmental factors on Loess Plateau. (**a**) Altitude. (**b**) Slope. (**c**) Slope aspect. (**d**) Geomorphic type. (**e**) Soil type. (**f**) Vegetation type.

**Figure 6.** Spatial distribution of GDP density on Loess Plateau during 2000–2015.

**Figure 7.** Spatial distribution of population density on Loess Plateau during 2000–2015.

**Figure 8.** Spatial distribution of land-use type on Loess Plateau during 2000–2015.

The following points should be noted: (1) Although environmental factors do not change in the short-term, these factors play a key role in vegetation growth and changes. Thus, they were also considered in this study according to [37]. (2) The data for each factor were discretized since the input of the GDM needs to be discrete data. (3) The data for the factors were converted into grid data with the 1 km × 1 km spatial resolution (same with the NDVI data) in order to facilitate calculation and analysis. (4) Due to the large area of the LP, the GDM's calculation capacity would be exceeded if all of the grid data were used. Therefore, the sampling method was used and 32,000 sampling points were evenly selected in space to carry out the calculations and analysis (Figure 1), and the consistency of the results of each calculation could be ensured. (5) Two types of inputs were considered in the analysis of the explanatory power of different factors for the spatial–temporal evolution of the NDVI on the LP using the GDM. First, the typical annual values of 2000, 2005, 2010, and 2015 were used to drive the GDM to analyze the influence of various factors on the spatial distribution state of the NDVI. Second, the differences between the annual mean values from 2000 to 2005 and those from 2010 to 2015 were calculated, which were used to drive the GDM to analyze the influence of various factors on the spatial distribution change of the NDVI.

#### **3. Results**

#### *3.1. Spatial–Temporal Evolution in NDVI*

(1) Grid scale area and area transfer of the NDVI and the driving forces:

The vegetation types on the LP are shown in Figure 5f. In the semi-arid area, the vegetation types were mainly grassland, cultivated vegetation, and desert, accounting for 15.9, 7.2 and 5.0% of the total area, respectively. While in the semi-humid area, cultivated vegetation, grassland, shrub, and broad-leaved forest were widely distributed, accounting for 37.8, 9.2, 6.7 and 5.4% of the total area, respectively.

The grid scale area transfer matrix of the NDVI on the LP during 2000–2015 is shown in Table 3. In 2000, over the entire LP, the coverage areas of the low, relatively low, medium,

relatively high, and high vegetation accounted for, respectively, 10.0, 31.5, 29.3, 28.1 and 1.3% of the total area; in the semi-arid area, these types of vegetation accounted for 9.4, 16.0, 5.1, 1.2 and 0.01%, respectively; and in the semi-humid area, they accounted for 0.5, 15.5, 24.2, 26.9 and 1.3%, respectively. In 2015, over the entire LP, the coverage areas of the low, relatively low, medium, relatively high, and high vegetation accounted for 8.7, 20.3, 22.9, 31.5 and 16.7%, respectively; in the semi-arid area, they accounted for 8.5, 14.5, 5.0, 3.3 and 0.4%, respectively; and in the semi-humid area, they accounted for 0.3, 5.8, 17.9, 28.1 and 16.3%, respectively. The proportion of the coverage areas with an NDVI of greater than 0.6 increased significantly, and the proportions of the relatively low vegetation converted to medium vegetation, medium vegetation converted to relatively high vegetation, and relatively high vegetation converted to high vegetation accounted for the largest proportions. The areas of these transitions accounted for, respectively, 10.8, 15.8 and 14.1% of the total area.


**Table 3.** Grid scale area transfer matrix of NDVI on Loess Plateau during 2000–2015 (km2).

The land-use type on the LP during 2000–2015 is shown in Figure 8, and the grid scale area transfer matrix is shown in Table 4. In 2000, over the entire LP, the proportions of the farmland, woodland, grassland, water bodies, construction land, and unused land were 33.3, 14.7, 41.5, 1.4, 2.5 and 6.7%, respectively; in the semi-arid area, they were 7.1, 1.2, 16.6, 0.7, 0.8 and 5.3%, respectively; and in the semi-humid area, they were 26.2, 13.5, 24.9, 0.8%, 1.7 and 1.4%, respectively. In 2015, over the entire LP, the proportions of farmland, woodland, grassland, water bodies, construction land, and unused land were 32.4, 15.1, 41.2, 1.5, 3.3 and 6.6%, respectively; in the semi-arid area, they were 6.9, 1.3, 16.4, 0.7, 1.1 and 5.2%, respectively; and in the semi-humid area, they were 25.4, 13.8, 24.8, 0.8, 2.1 and 1.4%, respectively. Overall, the area of the change in land-use types on the LP was not significant, and it was mainly in the 0.94% decrease in farmland, the 0.44% increase in woodland, and the 0.82% increase in construction land. To a certain extent, these changes reflect the impacts of human activities in the study area, such as the implementation of the GGP and urbanization-related development.

The difference in the NDVI of the different land-use type transfers on the LP during 2000–2015 is shown in Figure 9. It should be noted that Figure 9 was obtained by subtracting the NDVI data for 2000 from that for 2015 on the grid scale, and the differences in the NDVI of the different land-use type transfers were then counted. Figure 9a shows the mean difference in the NDVI of all grids of different land-use type transfers. Most of the transfers had a positive impact on the NDVI, except for the transfers of farmland and water bodies to construction land, which decreased the NDVI slightly (reduction rates of −0.0173 and −0.0016, respectively). Moreover, the largest increment of the NDVI was due to farmland transferred into woodland (increment of 0.1411), woodland transferred into farmland (increment of 0.1359), and unused land transferred into farmland (increment of 0.1338). Figure 9b shows the total difference in the NDVI of all grids of different land-use type transfers. It can be seen that the largest increments of the NDVI were due to grassland, farmland, and woodland, and these land-use types did not change. The increments were 28,539.55, 25,545.42, and 12,450.86, respectively. The probable reasons are as follows. On the LP, the precipitation is limited and there exists a gap between agricultural water demand and supply. In the past, there used to be large areas of farmland that could not be irrigated adequately and were greatly affected by the precipitation. With the implementation of the GGP, some infertile areas of farmland that could not be irrigated adequately were returned to the woodland or grassland, and what remained was adequately irrigated or fertile areas of farmland that were less affected by the precipitation. Moreover, with the development of agricultural technologies and the optimization and adjustment of planting structures, the crop yields on the LP were continuously increased.

**Table 4.** Grid scale area transfer matrix of land-use type on Loess Plateau during 2000–2015 (km2).


(2) Spatial changes in the NDVI on the grid scale:

The spatial changes in the NDVI on the grid scale on the LP during 2000–2015 are shown in Figure 10. Figure 10a indicates that the annual mean value of the NDVI on the LP was 0.529, with an uneven spatial distribution, i.e., decreasing from the southeastern to the northwestern areas. The high and low values were mainly distributed in the semi-humid and semi-arid areas, with annual mean values of 0.619 and 0.346, respectively.

Figure 10b shows that the mean value of the coefficient of variation of the NDVI was 0.1406 on the LP, with an uneven spatial distribution, i.e., increasing from the southeastern semi-humid area (mean value of 0.1165) to the northwestern semi-arid area (mean value of 0.1926). There was 7.0% of the total area with 0 < CV < 0.05, mainly distributed in the eastern and southern regions. The area with 0.05 ≤ CV < 0.1 accounted for 25.9% of the total area. The area with 0.1 ≤ CV < 0.15 accounted for 25.7% of the total area. Moreover, there was 41.4% of the total area with 0.15 ≤ CV < 0.2 and CV ≥ 0.2, mainly distributed in the semi-arid area as well as the junction of the semi-arid/humid areas.

(**b**)

**Figure 9.** Difference in NDVI of different land-use type transfers during 2000–2015. (**a**) Mean difference in NDVI of all grids of different land-use type transfers. (**b**) Total difference in NDVI of all grids of different land-use type transfers.

Figure 10c,d show that of the total area, the positive area accounted for 91.5% while the negative area accounted for 8.5%, which was obtained by subtracting the NDVI in 2000 from that in 2015. The NDVI in most areas of the LP (92.8% of the total area) increased from 2000 to 2015, indicating that the ecological environment in the region significantly improved. In addition, 70.0% of the area passed the significance test (*p* < 0.05); and 68.6% of the area increased significantly, while 1.4% of the area decreased significantly. The decrease in the NDVI was mainly concentrated in the semi-arid area as well as the urban area with rapid economic development and a large population concentration.

**Figure 10.** Spatial change in NDVI on grid scale on Loess Plateau during 2000–2015.

(3) Temporal changes in the NDVI on the regional scale:

The temporal changes in the NDVI on the regional scale on the LP during 2000–2015 are shown in Figure 11. The NDVI on the LP exhibited a fluctuating upward trend. The annual NDVI growth rates of the entire LP, the semi-arid area, and the semi-humid area were 0.0079, 0.0049, and 0.0093, respectively, indicating that the growth rate of the NDVI in the semi-humid area was higher than in the semi-arid area after the GGP was implemented.

**Figure 11.** Temporal change in NDVI on regional scale on Loess Plateau during 2000–2015.

#### *3.2. Individual Effects of Factors on NDVI*

(1) Results of typical years:

Using the typical annual values of 2000, 2005, 2010, and 2015 to drive the factor detector, the individual effects of the factors (represented by the *q* value) could be obtained, as shown in Figure 12.

For the entire LP, the order of the annual mean *q* values was precipitation (*q* = 0.5985) > vegetation type (*q* = 0.4790) > soil type (*q* = 0.3346) > land-use type (*q* = 0.2697) > geomorphic type (*q* = 0.2341) > temperature (*q* = 0.1469) > altitude (*q* = 0.1203) > population density (*q* = 0.1132) > slope (*q* = 0.0649) > GDP density (*q* = 0.0411) > slope aspect (*q* = 0.0012).

For the semi-arid area, the order of the annual mean *q* values was precipitation (*q* = 0.4796) > vegetation type (*q* = 0.3906) > soil type (*q* = 0.3409) > land-use type (*q* = 0.2368) > geomorphic type (*q* = 0.2013) > temperature (*q* = 0.1851) > altitude (*q* = 0.1695) > population density (*q* = 0.0954) > slope (*q* = 0.0774) > GDP density (*q* = 0.0110) > slope aspect (*q* = 0.0019).

For the semi-humid area, the order of the annual mean *q* values was geomorphic type (*q* = 0.2597) > vegetation type (*q* = 0.2466) > precipitation (*q* = 0.2250) > land-use type (*q* = 0.2053) > soil type (*q* = 0.1637) > temperature (*q* = 0.0965) > population density (*q* = 0.0783) > slope (*q* = 0.0541) > altitude (*q* = 0.0512) > GDP density (*q* = 0.0397) > slope aspect (*q* = 0.0086).

(2) Results of annual mean differences:

Using the differences between the annual mean values from 2000 to 2005 and those from 2010 to 2015 to drive the factor detector, the individual effects of the factors (represented by the *q* value) could be obtained, as shown in Figure 13.

For the entire LP, the order of the *q* values was soil type (*q* = 0.1287) > vegetation type (*q* = 0.1142) > land-use type (*q* = 0.0729) > temperature (*q* = 0.0680) > geomorphic type (*q* = 0.0532) > precipitation (*q* = 0.0508) > GDP density (*q* = 0.0463) > altitude (*q* = 0.0357) > population density (*q* = 0.0224) > slope (*q* = 0.0020) > slope aspect (*q* = 0.0010).

**Figure 12.** *Cont*.

**Figure 12.** Effects of individual factors derived from the factor detector with inputs of typical annual values. (**a**) Entire area. (**b**) Semi-arid area. (**c**) Semi-humid area.

**Figure 13.** Effects of individual factors derived from the factor detector with inputs of differences between annual mean values.

For the semi-arid area, the order of the *q* values was precipitation (*q* = 0.1476) > soil type (*q* = 0.1373) > vegetation type (*q* = 0.1353) > land-use type (*q* = 0.0864) > temperature (*q* = 0.0500) > geomorphic type (*q* = 0.0417) > altitude (*q* = 0.0358) > GDP density (*q* = 0.0202) > population density (*q* = 0.0195) > slope aspect (*q* = 0.0030) > slope (*q* = 0.0010).

For the semi-humid area, the order of the *q* values was GDP density (*q* = 0.1583) > geomorphic type (*q* = 0.1101) > soil type (*q* = 0.0999) > altitude (*q* = 0.0955) > land-use type (*q* = 0.0761) > temperature (*q* = 0.0705) > population density (*q* = 0.0520) > precipitation (*q* = 0.0452) > vegetation type (*q* = 0.0227) > slope aspect (*q* = 0.0049) > slope (*q* = 0.0040).

In summary, the explanatory powers of the factors on the spatial distribution state of the NDVI in typical years and the spatial distribution change of the NDVI during 2000–2015 were different. In addition, the explanatory powers were different for the entire LP, and for the semi-arid/humid areas. With respect to the spatial distribution state of the NDVI in typical years, the decisive climatic factor was precipitation, the decisive environmental factors were geomorphic type, soil type, and vegetation type, and the decisive anthropogenic factor was land-use type. With respect to the spatial distribution change of the NDVI during 2000–2015, in the semi-arid area, the climatic and environmental factors were the decisive factors, including precipitation, soil type, and vegetation type; the impacts of anthropogenic factors, such as the GDP density, land-use type, and population density, were more significant in the semi-humid area. *3.3. Interactive Effects between Factors on NDVI*

(1) Results of typical years:

Using the typical annual values of 2000, 2005, 2010, and 2015 to drive the interaction detector, the interactive effects of the factors (represented by the *q* value) could be obtained, as shown in Figure 14. For the entire LP and the semi-arid/humid areas, the interactive effects between factors were greater than their individual effects, indicating that none of the factors acted independently, but they had a certain enhancement effect, including nonlinear enhancement and bi-enhancement.

For the entire LP, 26.4% of the interactive factor combinations exhibited nonlinear enhancement and 73.6% exhibited bi-enhancement. The interactive effect between precipitation and vegetation type was the strongest (mean annual?*q* = 0.7034), followed by precipitation and soil type (mean annual?*q* = 0.6987).

For the semi-arid area, 32.3% of the interactive factor combinations exhibited nonlinear enhancement and 67.7% exhibited bi-enhancement. The interactive effect between precipitation and soil type was the strongest (mean annual?*q* = 0.6375), followed by precipitation and vegetation type (mean annual?*q* = 0.6235).

For the semi-humid area, 44.1% of the interactive factor combinations exhibited nonlinear enhancement and 55.9% exhibited bi-enhancement. The interactive effect between precipitation and geomorphic type was the strongest (mean annual?*q* = 0.3961), followed by precipitation and vegetation type (mean annual?*q* = 0.3701).

(2) Results of annual mean differences:

Using the differences between the annual mean values from 2000 to 2005 and those from 2010 to 2015 to drive the interaction detector, the interactive effects of the factors (represented by the *q* value) could be obtained, as shown in Figure 15.

(**a**)

(**b**)

**Figure 14.** *Cont*.

(**c**)

**Figure 14.** Effects of interactive factors derived from the interaction detector with inputs of typical annual values. (Note: \* indicates bi-enhance; \*\* indicates nonlinear enhance). (**a**) Entire area. (**b**) Semiarid area. (**c**) Semi-humid area.

For the entire LP, 90.9% of the interactive factor combinations exhibited nonlinear enhancement and 9.1% exhibited bi-enhancement. The interactive effect between soil type and temperature was the strongest (*q* = 0.2194), followed by soil type and vegetation type (*q* = 0.2164).

For the semi-arid area, 83.6% of the interactive factor combinations exhibited nonlinear enhancement and 16.4% exhibited bi-enhancement. The interactive effect between precipitation and soil type was the strongest (*q* = 0.2438), followed by soil type and land-use type (*q* = 0.2379).

For the semi-humid area, 61.8% of the interactive factor combinations exhibited nonlinear enhancement and 38.2% exhibited bi-enhancement. The interactive effect between GDP density and geomorphic type was the strongest (*q* = 0.2470), followed by GDP density and soil type (*q* = 0.2300).

(**a**)

(**b**)

**Figure 15.** *Cont*.

**Figure 15.** Effects of interactive factors derived from the interaction detector with inputs of differences between annual mean values. (**a**) Entire area. (**b**) Semi-arid area. (**c**) Semi-humid area.

#### *3.4. Ranges or Types of Factors for NDVI*

Using the typical annual values for 2000, 2005, 2010, and 2015 to drive the risk detector, the ranges or the types of factors for the NDVI could be obtained, as shown in Figure 16. The ranges or the types of factors for the NDVI were different in different years in the entire LP, and in the semi-arid/humid areas.

For the entire LP, the suitable precipitation range for the NDVI was greater than 650 mm (as the precipitation increased, the NDVI increased), the temperature range was smaller than 0 ◦C (as the temperature increased, the NDVI decreased first and then increased), the altitude ranged from 90 to 790 m or from 2963 to 4914 m (as the altitude increased, the NDVI decreased first and then increased), the slope ranged from 15◦ to 20◦, the slope aspect was no slope aspect, the geomorphic type was large undulating mountain, the soil type was eluvial soil, the vegetation type was broad-leaved forest, the GDP density ranged from 500 to 1000 CNY/km2, the population density ranged from 200 to 500 people/km2, and the land-use type was woodland.

For the semi-arid area, the suitable precipitation for the NDVI ranged from 450 to 550 mm, the temperature range was smaller than 0 ◦C, the altitude ranged from 2963 to 4914 m, the slope ranged from 20◦ to 25◦, the slope aspect was no slope aspect, the geomorphic type was large undulating mountain, the soil type was anthropogenic soil, the vegetation type was cultivated vegetation, the GDP density ranged from 500 to 1000 CNY/km2, the population density ranged from 200 to 500 people/km2, and the land-use type was farmland.

For the semi-humid area, the suitable precipitation range for the NDVI was greater than 650 mm, the temperature range was smaller than 0 ◦C, the altitude ranged from 2963 to 4914 m, the slope ranged from 15◦ to 20◦, the slope aspect was no slope aspect, the geomorphic type was large undulating mountain, the soil type was eluvial soil, the vegetation type was broad-leaved forest, the GDP density ranged from 3000 to 5000 CNY/km2, the population density ranged from 200 to 500 people/km2, and the land-use type was woodland.

**Figure 16.** *Cont*.

**Figure 16.** Ranges or types of factors for NDVI derived from the risk detector with inputs of typical annual values. (**a**) Entire area. (**b**) Semi-arid area. (**c**) Semi-humid area.

#### *3.5. Differences of Significance between Factors on NDVI*

Using the typical annual values for 2000, 2005, 2010, and 2015 to drive the ecological detector, the differences of significance between factors for the NDVI could be obtained, as shown in Table 5. For the entire LP, the differences were not significant between the precipitation and other factors in the different years (i.e., relatively certain). The significant differences between the soil type and the other factors did not vary in the different years (i.e., relatively certain). The soil type showed no significant differences with precipitation, GDP density, population density, and land-use type, but it exhibited significant differences with other factors. The differences of significance varied between the other pairs of factors in the different years (i.e., relatively uncertain).

**Table 5.** Statistical significance of factors derived from the ecological detector with inputs of typical annual values.


Note: Y indicates that there is a significant difference in the effect of two factors on NDVI (confidence is 95%); N indicates no significant difference; N/Y indicates that there is (or is no) significant difference across different years.

#### **4. Discussion**

#### *4.1. Comparison between Semi-Arid and Semi-Humid Areas*

Climatic factors influence the environmental conditions of the LP; while environmental factors determine the range and intensity of the human activities. The three types of factors are not independent of each other, and they interact with each other to influence the spatial distribution of the NDVI and its evolution over time on the LP. Based on the above results, there were significant spatial differences in the NDVI and its driving forces on both sides of the 400 mm isohyet on the LP.

For the semi-arid area to the northwest of the 400 mm isohyet, a number of observations were made. (1) The analysis of the spatial–temporal evolution revealed that the precipitation is relatively small and the temperature is relatively low, with annual mean values of 294 mm and 7.6 ◦C, respectively. The geomorphic types are mainly hills and plains (accounting for 12.3 and 11.0%, respectively). The soil types are mainly primary soil and arid soil (accounting for 11.1 and 8.0%, respectively). The main vegetation types are grassland, cultivated vegetation, and desert (accounting for 15.9, 7.2 and 5.0%, respectively). The grassland and farmland are the main land-use types (accounting for 16.4% and 6.9% in 2015, respectively). Moreover, the GDP density is relatively small (772 CNY/km2 in 2015). The population density is relatively small (96 people/km<sup>2</sup> in 2015). The vegetation coverage is relatively low (annual mean NDVI of 0.346), and the coefficient of variation is relatively large (mean value of 0.1926). The growth rate of the NDVI was 0.0049/a after 2000 when the GGP was implemented. (2) The analysis based on the GDM revealed the following: The individual factors that determine the changes in the spatial distribution of the NDVI are the climatic and environmental factors, such as precipitation (*q* = 0.1476), soil type (*q* = 0.1373), vegetation type (*q* = 0.1353), temperature (*q* = 0.0500), and geomorphic type (*q* = 0.0417). The main interactive factors are precipitation and soil type (*q* = 0.2438), and soil type and land-use type (*q* = 0.2379). The suitable precipitation for the NDVI ranged from 450 to 550 mm, the temperature range was smaller than 0 ◦C, the altitude ranged from 2963 to 4914 m, the slope ranged from 20◦ to 25◦, the slope aspect was no slope aspect, the geomorphic type was large undulating mountain, the soil type was anthropogenic soil, the vegetation type was cultivated vegetation, the GDP density ranged from 500 to 1000 CNY/km2, the population density ranged from 200 to 500 people/km2, and the land-use type was farmland.

Observations were also made regarding the semi-humid area to the southeast of the 400 mm isohyet. (1) The analysis of the spatial–temporal evolution revealed that both the precipitation and the temperature are relatively high, with annual mean values of 531 mm and 8.3 ◦C, respectively. The geomorphic types are mainly hills and medium undulating mountains (accounting for 19.4 and 13.3%, respectively). The soil types are mainly primary soil and semi-eluvial soil (accounting for 30.9 and 14.9%, respectively). The main vegetation types are cultivated vegetation and grassland (accounting for 37.8 and 9.2%, respectively). The farmland and grassland are the main land-use types (accounting for 25.4 and 24.8% in 2015, respectively). Moreover, the GDP density is relatively large (834 CNY/km2 in 2015). The population density is relatively large (212 people/km2 in 2015). The vegetation coverage is relatively high (annual mean NDVI of 0.619), and the coefficient of variation is relatively small (mean value of 0.1165). The growth rate of the NDVI was 0.0093/a after the implementation of the GGP in 2000. (2) The analysis based on the GDM revealed the following. The individual factors that determine the changes in the spatial distribution of the NDVI are the anthropogenic factors, such as GDP density (*q* = 0.1583), land-use type (*q* = 0.0761), and population density (*q* = 0.0520). The main interactive factors are GDP density and geomorphic type (*q* = 0.2470), and GDP density and soil type (*q* = 0.2300). The suitable precipitation range for the NDVI was greater than 650 mm, the temperature range was smaller than 0 ◦C, the altitude ranged from 2963 to 4914 m, the slope ranged from 15◦ to 20◦, the slope aspect was no slope aspect, the geomorphic type was large undulating mountain, the soil type was eluvial soil, the vegetation type was broad-leaved forest, the

GDP density ranged from 3000 to 5000 CNY/km2, the population density ranged from 200 to 500 people/km2, and the land-use type was woodland.

The conclusions above have important implications for ecological conservation and restoration in the different regions. The impacts of the climatic, environmental, and anthropogenic factors need to be comprehensively considered when formulating and implementing policies. The NDVI is sensitive to the climatic and environmental factors related to water changes in the semi-arid area where there is relatively little water. Therefore, appropriate vegetation types should be selected, regional vegetation structures should be optimized, and high water consumption vegetation should be gradually replaced by low water consumption vegetation such as grass and drought-tolerant plants. In the semihumid area, the regional water conditions basically meet the needs of vegetation growth because of the suitable climatic and natural conditions. Moreover, human activities have a relatively large disturbance effect on the NDVI. Thus, more positive human activities should be encouraged, including Grain for Green, afforestation, urban greening, and agricultural modernization projects. Furthermore, negative human activities should be controlled, including disordered urban expansion, population expansion, and overgrazing.

#### *4.2. Connections and Distinctions with Other Studies*

(1) Vegetation NDVI: In this study, it was found that the NDVI on the entire LP has increased significantly after 2000 when the GGP was implemented, consistent with the findings of [26,37,38].

(2) Climatic factors: In this study, it was found that among all factors, precipitation has the greatest influence on vegetation growth in the semi-arid area of the LP, consistent with the findings of [39–41]. In the semi-humid area, precipitation is not the decisive factor controlling vegetation growth because of the relatively good water conditions. Several studies have suggested that temperature is also a key factor affecting the vegetation growth. In general, when the temperature is low, the physiological activity of the vegetation is low. An increase in temperature promotes photosynthesis and vegetation growth, but as the temperature increases further, the growth of vegetation is inhibited by the accelerated evaporation and soil drying [42,43]. It was found in this study that conditions are suitable for vegetation growth when the temperature smaller than 0 ◦C on the LP. This is the case because the vegetation growth is influenced by many factors other than temperature. It can be seen from Figures 2–5 and 8 that the NDVI value is large in the southwest corner of semi-arid area of the LP, where the precipitation is large, the land-use types are mainly woodland and grassland, and the altitude ranges from 2963 m to 4914 m (the high altitude causes the low temperature). The finding is consistent with the results of [29]. It was also found that compared with other factors, temperature has relatively small explanatory power in terms of the spatial distribution and its change in the NDVI on the LP, which agrees with the findings and results of [39].

(3) Environmental factors: In this study, it was found that the environmental factors such as topographic types have an impact on the vegetation growth and restoration, consistent with the findings and results of [44]. The study revealed that the explanatory power of interactions between the factors would be increased for the spatial distribution and its change in the NDVI, which is consistent with the practical situation. For instance, the soil moisture is lower for steeper slopes, inhibiting vegetation growth. In this study, it was found that the NDVI value is large when the slope is between 15◦ and 20◦, while it decreases when the slope becomes greater.

(4) Human factors: Certain studies have suggested that the regional land-use types have changed since the implementation of the GGP [45,46], which is not consistent with the results of this study. The reason for this may be that the spatial resolutions of the remote sensing data used in these studies are different. In this study, it was found that the conservation and restoration of the original woodland and grassland have been strengthened since the implementation of the GGP on the LP, resulting in a significant increase in the NDVI throughout the entire region. Additionally, many studies have used the NDVI to

predict the grain yield because of the significant correlation between the grain yield and the NDVI [47,48]. The LP is a dry farming area, which has been affected by the scarcity of precipitation for a long time, resulting in a low grain yield [49]. With the agricultural technology developing recently, the optimization and adjustment of the crop planting structure has resulted in a continuous improvement in the grain yield on the LP. Therefore, the NDVI has increased significantly [19].

#### *4.3. Possible Future Work*

There is further research that can be carried out, including the following. (1) The factors affecting the spatial–temporal evolution of the NDVI are extremely complex. Due to data limitations, certain factors are not currently considered, such as agricultural fertilization, irrigation area, and CO2 concentration [50,51], but can be further improved with the acquisition of new data in the future. (2) The interaction detector of the GDM only considers the interactions between two factors. There are interactive effects involving more than two factors. The combined impacts of multiple factors would provide insights on the evolution of vegetation, which can be considered in the future by introducing other methods, such as the analytical hierarchy process (AHP) and principal component analysis (PCA) [52]. (3) It is necessary to perform further collection of data with a relatively higher resolution through field investigations and monitoring to determine the practical situation of the surface undulations in the study area. (4) The cumulative and time-lag effects of various factors and indicators (e.g., standardized precipitation and evapotranspiration index, SPEI) should also be evaluated in the future [53]. (5) It is possible for the GDM to identify the impacts of driving forces on the evolution of vegetation from the perspective of spatial analysis and mathematical statistics. In order to further explore the evolutionary mechanism of vegetation, some methods from the field of landscape studies can be introduced in the future [54].

#### **5. Conclusions**

The main conclusions of this paper include the following.

(1) The spatial–temporal evolution characteristics of the NDVI on the entire LP, and in the semi-arid/humid areas during 2000–2015 were analyzed via the linear regression, coefficient of variation, and transfer matrix models.


(2) The GDM was adopted to quantitatively identify the impact of multiple factors on the spatial–temporal evolution of the NDVI on the entire LP, and in the semi-arid/humid areas.

• Using the factor detector, it was found that in the semi-arid area, the climatic and environmental factors were the decisive factors influencing the spatial distribution changes of the NDVI during 2000–2015, including precipitation, soil type, and vegetation type. The impacts of anthropogenic factors, such as the GDP density, land-use type, and population density, were more significant in the semi-humid area.


The conclusions of this study have important implications for policy makers and administrative managers in terms of the formulation and implementation of ecological conservation and restoration strategies in the different regions. In the semi-arid area, appropriate vegetation types should be selected, regional vegetation structures should be optimized, and high water consumption vegetation should be gradually replaced by low water consumption vegetation such as grass and drought-tolerant plants. In the semihumid area, more positive human activities should be encouraged, including Grain for Green, afforestation, urban greening, and agricultural modernization projects. Furthermore, negative human activities should be controlled, including disordered urban expansion, population expansion, and overgrazing.

**Author Contributions:** Conceptualization, Data curation, Investigation, Methodology, Formal analysis, Writing—original draft: Y.D.; Conceptualization, Methodology, Formal analysis, Writing—review and editing, Supervision, Funding acquisition: D.Y.; Validation, Writing—review and editing: X.L. (Xiang Li); Project administration, Funding acquisition. J.H.; Supervision, Funding acquisition: W.S.; Validation, Writing—review and editing: X.L. (Xuecao Li); Supervision, Funding acquisition: H.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Natural Science Foundation of China (42171024, 51609122, 41701507), the National Key Research and Development Program of China (2018YFE0122700), the Open Research Fund Program of State Key Laboratory of Hydroscience and Engineering (sklhse-2020-A-07), and the Joint Open Research Fund Program of State key Laboratory of Hydroscience and Engineering and Tsinghua-Ningxia Yinchuan Joint Institute of Internet of Waters on Digital Water Governance (sklhse-2021-Iow06). Su Wei would like to thank the 2115 Talent Development Program of China Agricultural University.

**Data Availability Statement:** The shape file for the boundary of the LP, the SPOT NDVI data, and the influencing factors (including geomorphic type, soil type, vegetation type, socio-economic data, and land-use type) data were downloaded from the Resource and Environment Science and Data Center (RESDC), Chinese Academy of Sciences (available from http://www.resdc.cn, accessed on 6 October 2021). The DEM data were downloaded from the Shuttle Radar Topography Mission (available from http://srtm.csi.cgiar.org/, accessed on 6 October 2021). The precipitation and temperature data were from the China Meteorological Data Service Center (available from http: //data.cma.cn/site/index.html, accessed on 6 October 2021).

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

#### **References**

