Analysis of the Spatial and Temporal Evolution Patterns of Grassland Health and Its Driving Factors in Xilingol

Abstract

The combination of natural environment changes and human activities affects the growth of grasslands. In order to quantitatively assess the causes of spatial and temporal variation of grasslands in Xilingol, this study assessed the spatial and temporal evolution patterns of grassland health based on MOD13A1 long time series Normalized Difference Vegetation Index (NDVI) data from 2000–2019 using trend analysis. The geodetector model was used to explore the dominant drivers of spatial variation in grassland NDVI, combined with 34 factors covering natural environmental changes and human disturbances over the same period. The results show that the grasslands of Xilingol showed an overall recovery trend from 2000 to 2019, with an average annual NDVI growth rate of 0.0028/a, a monthly increasing rate of 0.0005/month, and 68.06% of the grassland at an average recovery level. Moisture-dominated natural climate change factors, such as Growing Season Precipitation (Prep2), Annual Mean Water Vapor Pressure (WVP), and Annual Mean Relative Humidity (RH), were the underlying cause of grassland health changes during the study period, with the highest explanatory factor being growing season precipitation (q value of 0.59 on a multi-year average). The influence of primary production value among human activities was greater, and the explanatory factor of tertiary production value showed an increasing trend. The interactions among natural and anthropogenic factors significantly enhances their explanatory credibility for NDVI, with the type of interaction dominated by the two-factor enhancement. Risk detection of the top 10 dominant drivers in terms of q statistic were carried out to obtain the threshold range of each driver in the high zone of grassland NDVI, which can provide a scientific reference for the sustainable restoration of grassland.

1. Introduction

As one of the most important ecosystems, grasslands not only contain huge natural resources, but also participate heavily in the global carbon cycle, and are closely related to human production and life. Therefore, it is of great significance for regional sustainable development to study the spatial and temporal patterns of grassland health and its driving mechanisms [1,2]. Normalized Difference Vegetation Index (NDVI) is an important evaluation index to respond to the growth status of grassland vegetation [3], and its annual maximum value can effectively respond to the above-ground biomass [4,5] and other grassland health indicators. Previous studies investigated vegetation cover changes based on time series NDVI data using Sen’s slope and Mann–Kendall test methods [6,7], and analyzed the response of NDVI to natural environment conditions using statistical methods. The results showed that climatic change explained vegetation growth better, whereas for grassland, a fragile vegetation type, in addition to natural factors such as precipitation and temperature, the impact of human activities such as sheep stocking, large livestock stocking, and livestock density on grassland health also needs to be considered. [8,9,10].

The geodetector model was developed by Wang et al. [11,12,13] based on Spatial Stratified Heterogeneity (SSH), which is a tool to measure Spatial Heterogeneity (SH) and to explore the determinants of SH. It can quantify the explanatory power of factors and the power of interactions. In recent years, Tian et al. and Zhou et al. have quantitatively evaluated the drivers of grassland health, but mainly analyzed from a single perspective, such as meteorological factors [14,15] and human activity impacts [16], without considering the integration of driving factors on grassland NDVI. The geodetector model can aptly solve the above problems simply and quickly [17,18,19].

The Xilingol grassland is a typical area of arid semi-arid temperate grassland in northern China and a major component of the Inner Mongolia grassland. The region is currently the largest natural reserve of the grassland and meadow ecosystem type in China, and holds an important spatial position and international influence in the conservation of grassland biodiversity [20]. The Xilingol grassland is also the closest grassland pasture and source of wind and sand to Beijing, Tianjin, and Hebei, and is an essential ecological barrier for a vast inland region [21]. However, with overuse of grasslands due to increasing population pressure, grasslands are also facing a series of ecological and environmental health problems, such as productivity decline, area reduction, and land degradation, which have huge impacts on the balance of grassland ecosystems and the sustainable development of grassland grazing [22,23], as well as bringing risks to human health, such as the spread of plague [24].

Therefore, our study selects Xilingol in Inner Mongolia as the study area, and explores the spatial and temporal evolution patterns of grassland health in Xilingol League based on MOD13A1 NDVI time series data from 2000–2019. There were 34 factors covering climate, natural conditions, population, economy, etc., selected from both natural and human activities, and the impact of these factors on NDVI changes was then quantitatively assessed, and the suitable range of factors for healthy grassland was obtained through multi-year analysis, so as to provide scientific support for grassland ecological project construction.

2. Materials and Methods

2.1. Study Area and Data Source

2.1.1. Study Area

The Xilingol grassland is the hinterland of the Eurasian steppe belt, located in the central part of the Inner Mongolia Autonomous Region, with a geographical location of 115°13′–117°06′E and 43°02′–44°52′N (Figure 1). The total land area of the study area is 202,580 km², with a grassland coverage of more than 95% [25]. The elevation of the study area is between 761–1927 m, with an average elevation of 1104 m, and the topography is mainly plateau plains, accompanied by hills, shallow mountains, and other landscape types. The climate type is typical temperate arid semi-arid continental climate. The average annual precipitation is 200–300 mm, mainly concentrated in the vegetation growth season, and the precipitation gradually decreases from east to west, with obvious regional differences. The average annual temperature of the whole region is 1–2 °C, the average minimum temperature is −20 °C, and the average maximum temperature is 21 °C. The soil types are black calcium soil, chestnut calcium soil, weathered soil, and brown calcium soil. The study area has three types of grasslands: meadow grassland, typical grassland, and desert grassland, which are distributed from east to west. Meadow grassland is mainly distributed in the eastern mountainous region, with elevation above 1098 m. The region is rich in precipitation, with annual precipitation above 350 mm, and the annual average temperature is between 0.7–5.1 °C. Typical grassland is widely distributed with undulating terrain, ranging from 740–1440 m, with annual precipitation between 240–350 mm, and an annual average temperature between 3.1–4.5 °C. Desert grassland is more flat compared to the other grassland types, with elevation ranging from 924–1098 m, average annual precipitation less than 240 mm, and an average annual temperature greater than 4.5 °C.

2.1.2. Data Source

The data used in this study include remote sensing data and statistical yearbook data. Among them, remote sensing data are mainly NDVI, and also include soil type data, land use data, etc.

MODIS NDVI data (MOD13A1: h25v04 and h26v04) for the period 2000–2019 May to September each year were downloaded from the Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/(accessed on 10 July 2022)), with 500 m spatial resolution and 16d temporal resolution. The NDVI data were geometrically corrected, stitched, projected, and formatted using MRT; after which, the data were cropped according to the study area. Pixel quality layers of MOD13A1 were used to eliminate invalid values and poor-quality pixels; after which, monthly and yearly NDVI data were extracted in ArcGIS by MVC (maximum value composite) for further analysis.

Meteorological data and DEM were provided by the Resource and Environment Science Data Centre of the Chinese Academy of Sciences (http://www.resdc.cn (accessed on 15 July 2022)); soil type and land use data were provided by National Tibetan Plateau Data Center (http://data.tpdc.ac.cn (accessed on 1 July 2022)); population density data were downloaded from the Worldpop (https://www.worldpop.org (accessed on 20 July 2022)); economic and livestock data were obtained from the statistical yearbooks provided by the Bureau of Statistics of Inner Mongolia (http://tj.nmg.gov.cn (accessed on 20 June 2022)) and Xilingol (http://tjj.xlgl.gov.cn (accessed on 25 June 2022)); and the corresponding data for 2000, 2005, 2010, 2015, and 2019 were downloaded. The DEM was adopted from the SRTMDEM with a spatial resolution of 90 m; the slope and aspect were obtained on this basis. Soil type data were obtained according to the soil classification method of SYMBOL90 of the Harmonized World Soil Database (HWSD) [26], with 31 soil types and a spatial resolution of 1 km.

In order to present the spatial heterogeneity of anthropogenic activities, socio-economic and livestock data at the banner and county scales in the study area were gridded from statistical yearbooks. Land cover data were also used to calculate the human activity intensity of land surface (HALS) to determine the intensity of human utilization, modification, and exploitation on the natural land surface [27]. As the majority of the data used in the study were at 1 km spatial resolution, and existing studies have taken the same 1 km spatial resolution in the region [28], to reduce the cost of operations, the most common nearest neighbor method was used to resample the above data to 1 km, using the Asia_North_Albers_Equal_Area_Conic projection. Due to the method used in the study requiring the input of category data, the above data needed to be reclassified into 8 categories from lowest to highest by data value field, except for soil type and land use type, which maintained their original categories, and the segmentation method used was the natural breaks (Jenks) method.

2.2. Materials

Xilingol grassland is ecologically vulnerable to variations of the natural environment and disturbances from human activities [29]. In this study, the spatial distribution and temporal changes of grassland health in Xilingol were analyzed based on long time series NDVI data, and the trends of grassland changes over the past 20 years since 2000 were calculated using trend analysis to obtain the spatial distribution of degraded and restored areas. NDVI data for 2000, 2005, 2010, 2015, and 2019 were selected as dependent variables, and the three independent variables of climate change, natural environment, and anthropogenic activity were combined to analyze the driving factors of grassland health using the geodetector model. The workflow of this research is shown in Figure 2.

2.2.1. Analysis of Vegetation Variation Trends

NDVI trends can reflect changes such as the degradation/restoration of land in the study area. Xilingol has more than 95% grass cover, so the analysis of NDVI can reflect the health of the grass. Stow D et al. found that when linear fit analysis of NDVI time series data was performed, the slope of the fitted trend line can reflect the vegetation trend over a long time series [30]. The NDVI time series data set obtained by smoothing with the maximum synthesis method was used to calculate the NDVI trends in Xilinhot City on an image-by-image basis over a 20-year period. The calculation method is as follows:

$$\mathrm{Slope}=\frac{\mathrm{n}\times {\sum }_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}\times {\mathrm{NDVI}}_{\mathrm{i},\mathrm{xy}}-{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}\times {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{NDVI}}_{\mathrm{i},\mathrm{xy}}}{\mathrm{n}\times {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{i}}^2-{\left({\sum }_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}\right)}^2}$$

(1)

where Slope represents the vegetation change trend over a 20-year period from 2000–2019; n represents the length of the time series from 2000–2019 with a value of 20; and i represents the year number, 1–20 from 2000–2019, respectively. NDVI_(i, xy) denotes the image element value of the raster image element in row x and column y on the maximum synthetic NDVI image in year i. When Slope < −0.1%, it indicates that the vegetation of this image element shows a degradation trend in the study time series, and the smaller the value, the more serious the degradation; when Slope is between ±0.1%, it indicates that there is no significant change in this image element in the study time series, and the vegetation growth in the region is relatively stable; when Slope > 0.1%, it indicates that this image element shows a recovery trend in the study time series, and the larger the value, the more obvious the recovery is.

2.2.2. Geodetector Model

Geodetector is a popular geostatistical model that analyzes spatial variations and reveals the driving factors behind them. Geodetector consists of four subdetectors: factor detector, interaction detector, risk detector, and ecological detector. In this study, the former three detectors are used to investigate the driving mechanisms behind NDVI change [11,12,31].

(1)

Factor detector

The factor detector is calculated by the following q statistic:

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

(2)

where 0 ≤ $q$ ≤ 1, and the larger the value, the greater the explanatory power of the factor. When the $q$ value is 0, it means factor has no relationship with NDVI. $h$ is the number of strata for variables or factors, $N$ represents the number of units in stratum $h$, and ${\sigma }_h^2$ and ${\sigma }^2$ denote the variance in the stratum $h$ and entire study area, respectively. $SSW$ and $SST$ are the sum of squares within the data and the total sum of squares, respectively.

(2)

Interaction detector

The interaction detector was used to identify the explanatory power of the different factors for NDVI when they were interacted with each other. The new variable can be obtained as $q\left(X1{\displaystyle \cap }X2\right)$, which indicates the explanatory power of the of the interaction of variables X1 and X2 on Y. By comparing the relationship among $q\left(X1{\displaystyle \cap }X2\right)$, $q\left(X1\right)$, and $q\left(X2\right)$, interactions can be classified into five types, as shown in

Table 1. Interaction detector types.

$\mathit{q} \mathrm{Value\; Comparison}$	Interaction
$q\left(X\mathbf{1}{\displaystyle \cap }X\mathbf{2}\mathbf{\right)}<Min\left(q\left(X\mathbf{1}\right),q\left(X\mathbf{2}\right)\right)$	Non-linear weakening
$Min\left(q\left(X\mathbf{1}\right),q\left(X\mathbf{2}\right)\right)<q\left(X\mathbf{1}{\displaystyle \cap }X\mathbf{2}\right)<Max\left(q\left(X\mathbf{1}\right),q\left(X\mathbf{2}\right)\right)$	Single-factor nonlinear attenuation
$q\left(X\mathbf{1}{\displaystyle \cap }X\mathbf{2}\right)>Max\left(q\left(X\mathbf{1}\right),q\left(X\mathbf{2}\right)\right)$	Two-factor enhancement
$q\left(X\mathbf{1}{\displaystyle \cap }X\mathbf{2}\right)=q\left(X\mathbf{1}\right)+q\left(X\mathbf{2}\right)$	Independent
$q\left(X\mathbf{1}{\displaystyle \cap }X\mathbf{2}\right)>q\left(X\mathbf{1}\right)+q\left(X\mathbf{2}\right)$	Non-linear enhancement

.

(3)

Risk detector

The mean NDVI values of each influence factor within different sub-regions are calculated, and t-statistics are used to determine whether the differences between the sub-regions are significant. Therefore, the risk detector is also commonly used to find the appropriate range and type of NDVI drivers.

(4)

Driving factor selection

In the study, the NDVI of grassland was used as the dependent variable, and 34 factors were selectedas independent variables (

Table 2. Grassland change driving factors.

Type	Factors	Abbreviation	Unit
Climate Change	Annual Mean Temperature	Temp1	°C
	Annual Mean Maximum Temperature	Temp2	°C
	Annual Mean Minimum Temperature	Temp3	°C
	Annual Maximum Temperature	Temp4	°C
	Annual Minimum Temperature	Temp5	°C
	Annual Precipitation	Prep1	mm
	Annual Mean Air Pressure	Pres1	hpa
	Annual Maximum Air Pressure	Pres2	hpa
	Annual Minimum Air Pressure	Pres3	hpa
	Annual Mean Water Vapor Pressure	WVP	hpa
	Annual Mean Relative Humidity	RH	%
	Annual Sunshine Hours	Sun	h
	Annual Mean Wind Speed	Wind	m/s
	Growing Season Mean Temperature	Temp6	°C
	Growing Season Maximum Temperature	Temp7	°C
	Growing Season Minimum Temperature	Temp8	°C
	Growing Season Precipitation	Prep2	mm
Fundamental Natural Environment	Elevation	Elev	m
	Slope	Slope	degree
	Aspect	Aspect	°
	Soil Type	Soil	categorical
Anthropogenic Activity	Land-use Type	LUT	categorical
	Population Density	Pop	people/km2
	Gross Domestic Product	GDP	10,000 yuan
	Primary Industry	PI	10,000 yuan
	Secondary Industry	SI	10,000 yuan
	Tertiary Industry	TI	10,000 yuan
	Per Capita GDP	PGDP	yuan
	Large Animals	LA	10,000 heads
	Sheep and Goats	SG	10,000 heads
	Livestock Density	LD	heads/km2
	Human Activity Intensity of Land Surface	HALS	%
	Cultivated Area	CA	%
	Road Density	RD	km/km2

) from climate change, fundamental natural environment, and anthropogenic activities, taking as reference the selection of factors in Xilingol [32,33,34]. Seventeen indicators from meteorological stations were selected for the climate change factor, including thirteen annual indicators, mainly selected for temperature, precipitation, humidity, barometric pressure, wind speed, sunshine hours, etc. The temperature indicators considered the annual average temperature, annual average high temperature, and annual average low temperature, which were used to represent the annual average of temperature, whereas the maximum and minimum temperatures considered the influence of extreme weather conditions on grass growth. Four growing season indicators were selected, mainly in terms of precipitation and air temperature. Four indicators were selected as fundamental natural environment factors, as there were differences among the natural habitats in terms of elevation, slope, aspect, and soil type. For anthropogenic activity indicators, 13 indicators, such as population density, economic and industrial structure, road density, and others, were considered, whereas for Xilingol, as a typical grazing area, unreasonable grazing activities also bring risks to the healthy development of grassland, so grazing-related indicators, such as livestock density, were chosen [35].

The indicators selected for this study were widely chosen to adequately consider the factors affecting NDVI changes, so that potential influences would not be ignored.

3. Results

3.1. Spatial and Temporal Variation Patterns of NDVI in Xilingol

3.1.1. Spatial Distribution Characteristics of NDVI

The spatial distribution of NDVI of the Xilingol grassland from 2000–2019 has obvious spatial heterogeneity (Figure 3). In general, the variation of NDVI ranged from −0.24 to 0.90, with a mean value of 0.39, showing a trend of low in the east and high in the west. The NDVI in the east and northeast is high, with a mean value of 0.68. This area is a meadow grassland, which is a transition zone from forest to grassland. The NDVI in the central area is in the middle area, with a mean value of 0.39, and the grassland type is a typical grassland, which is the main part of the Xilingol grassland. The NDVI in the west is low, with a mean value of 0.19, and the grassland type in this area is a desert grassland, with a dry climate, little rainfall, and sparse vegetation distribution. In the south, Taipusi Banner and Duolun County have high NDVI values, with a mean value of 0.47. This area is in the transition zone between grassland and farmland. There is one NDVI low value area in Dongwuzhumuqin Banner, which is the high wall of Ulagai, which is in the catchment area of the northeastern Inner Mongolia Plateau, so its NDVI value is low. Another NDVI low value area is in the southern part of Abaga Banner, which is a sandy grassland in the Hunsandak Sands running through the central part of Xilingol [10]. In this way, one can observe how the spatial differentiation of grassland types in Xilingol determines the spatial variability of NDVI.

3.1.2. Temporal Variation Characteristics of NDVI in Xilingol

The grassland NDVI in Xilingol showed a trend of a fluctuating increase on the time scale. The temporal variation of grassland NDVI from 2000 to 2019 is shown in Figure 4. Both the monthly and annual NDVI showed an increase, with an increasing rate of 0.0005/a and 0.0028/a. The mean value of grassland NDVI during the growing season over 20 years was 0.29, and the high value of NDVI occurred in July and August. The highest value of the monthly average NDVI occurred in August 2018, and the lowest value occurred in May 2006. Combining the hydrothermal conditions at the corresponding times revealed that the precipitation and the average temperature were both high in July and August 2018, whereas the hydrothermal conditions were poor in May 2006. This indicates that grassland NDVI is influenced by meteorological factors.

Based on the fluctuating increasing trend in grassland NDVI during the 20-year period, we divided this period into four phases, the first phase being 2000–2003. During this phase, grassland NDVI showed a significant increase, with an increment rate of 0.037/a. The second phase was 2003–2009. During this phase, NDVI showed a decrease, with a reduction rate of 0.013/a. The third phase was in 2009–2012. During this phase, NDVI showed an extremely significant increase, with an increment rate of 0.051/a. The last phase was 2012–2019, where NDVI showed a decrease again, with a reduction rate of 0.003/a.

3.1.3. Variation Trends of NDVI in Xilingol

The spatial distribution of grassland restoration and degradation over 20 years can be shown by the variation trends of NDVI. As shown in Figure 5a, the variation trend ranges from −0.0472 to 0.0371, with a mean value of 0.0028, indicating that grasslands within the study area of Xilingol showed an overall recovery trend from 2000 to 2019, but there is obvious spatial heterogeneity in grassland variation trends [36]. The areas with NDVI increase are concentrated in the northeast, central, and southern parts of the study area, and the areas with NDVI decrease are concentrated in the northern, eastern, and western parts of the study area.

Variation trends were divided into five types by percentage, with less than −1% as significant degradation, −1% to −0.1% as degradation, −0.1% to 0.1% as no significant change, 0.1% to 1% as restoration, and greater than 1% as significant restoration. The reclassification results of grassland variation trends are shown in Figure 5b. Combined with the statistical results of different grades of variation trends within each administrative division (

Table 3. Percentage of grassland variation trends in different zones.

	Significant Degradation	General Degradation	No Significant Change	General Restoration	Significant Restoration
Dongwuzhumuqin Banner Xiwuzhumuqin Banner Xilinhot	0.09%	7.32%	14.21%	75.01%	3.37%
	0.44%	7.40%	13.06%	77.52%	1.58%
	0.22%	8.07%	12.91%	75.50%	3.30%
Abaga Banner Sunitezuo Banner	0.04%	10.72%	18.44%	70.57%	0.23%
	0.00%	4.11%	34.00%	61.84%	0.05%
Suniteyou Banner Erlianhot Xianghuang Banner Zhengxiangbai Banner Zhenglan Banner	0.00%	6.74%	43.76%	49.35%	0.14%
	0.00%	1.22%	21.30%	77.20%	0.27%
	0.00%	16.01%	31.91%	52.00%	0.08%
	0.02%	14.49%	22.68%	62.58%	0.23%
	0.03%	6.33%	13.57%	75.19%	4.88%
Taipusi Banner Duolan County	0.07%	5.49%	7.57%	72.09%	14.77%
	0.11%	2.87%	3.44%	65.71%	27.87%

), it was found that the highest percentage of general restoration areas in the study area was 68.06% from 2000 to 2019, and the cities and counties were also heavily dominated by restoration, which further illustrates the overall recovery trend of the Xilingol grassland. Significant degradation areas accounted for 0.10%, mainly concentrated in the eastern part of Xiwuzhumuqin Banner and the center of Xilinhot. The degraded areas in Xiwuzhimuqin Banner are concentrated in the NDVI high-value area, which is influenced by human activities such as overgrazing [37]. The degraded areas in the Xilinhot city center are mainly affected by the expansion of urban areas. The general degradation areas are mainly distributed in the central area of desert grassland in Suniteyou Banner and the high-altitude area in the north of Abaga Banner. Significant restoration areas accounted for 2.32%, mainly concentrated in Duolun County, Taipusi Banner, Zhenglan Banner, and the north-central part of Xilinhot City, indicating that many ecological projects such as the return of cultivated land to forest and grass, the Green Wall of China, and the treatment of Beijing–Tianjin wind and sand sources over the years have achieved certain results, and it should be noted that the NDVI in the high-walled waters of Dongwuzhumuqin Banner in Ulagai also showed an increasing trend, which may be due to the shrinking water area [38].

3.2. The Relationship between NDVI and Driving Factors in Xilingol

3.2.1. Detection of Differences in the Explanatory Credibility of Different Factors

The geodetector model was applied to identify the impacts of driving factors on NDVI distribution. The results of factor detection (q statistic) reflect the explanatory credibility of each factor on NDVI. Figure 6 shows the average q statistic of each driving factor from 2000 to 2019. Prep2 has the strongest credibility, with a multi-year average q statistic of 0.59, followed by Soil > WVP > PI > RH, and the q statistics of the above factors on NDVI are greater than 0.5. In addition to natural conditions, the impact of the primary industry on the spatial distribution of NDVI should not be neglected. The stronger factors were LD > Wind > GDP > Prep1 > LA > SG > CA > Temp1 > PGDP > Temp2 > SI > TI> Temp6 > Temp3, with the q statistics of all the above factors being greater than 0.4. Temp8, HALS, Pres1, RD, Pres2, Pres3, Temp7, Temp4, Sun, Temp5, and LUT had fair explanatory credibility in affecting NDVI, with q values ranging from 0.2–0.4. Slope, Elev, Pop, and Aspect had insignificant effects on NDVI, with q values less than 0.12.

The influence of natural conditions dominated the spatial distribution of NDVI, with Prep2 having the strongest explanatory credibility. The interannual variation of 19 factors with q statistics greater than 0.4 was counted, and the results are shown in Figure 7. The interannual variation of the explanatory credibility of most meteorological factors for NDVI is relatively consistent, showing a trend of fluctuating changes. Q statistics of meteorological factors such as Temp, WVP, RH, and Wind are in the high-value zone in 2005 and 2015, in the low-value zone in 2000 and 2010, and in the medium-value zone in 2019. Among them, the explanatory power of Soil for NDVI is high and stable from year to year, fluctuating between 0.50 and 0.63. The explanatory credibility of Prep2 for NDVI showed an increasing trend, from 0.33 in 2000 to 0.70 in 2019. RH had the strongest explanatory power for NDVI in 2005 and 2015. Compared with Prep2, the explanatory credibility of temperature factors such as Temp1, Temp2, Temp3, and Temp6 was weaker for NDVI. In 2010, the q statistic of Prep2 decreased significantly, and the rest of the temperature factors changed to a lesser extent.

The explanatory credibility of human activity factors on the spatial distribution of NDVI is weaker compared with natural factors, and the interannual variation is larger. From 2000–2005, the explanatory credibility of all factors increased significantly; from 2005–2010, the explanatory power of all factors decreased significantly, except for the SI and the TI; from 2010–2015, the q statistics of SG, LD, and PGDP decreased, and the remaining factors rose; and from 2015–2019, the q statistics of GDP, PI, SI, and LA decreased and the remaining factors increased. The explanatory credibility of LD for NDVI was greater than 0.5 in 2005, 2015, and 2019, and its trend was more consistent with SG. In terms of economic structure, PI had the highest q statistic from 2000 to 2005, when PI remained in the high-value area, but had a decreasing trend. SI increased significantly and surpassed PI in 2010, and from 2015 to 2019, PI and SI q statistics all decreased, whereas TI was the largest. The explanatory credibility of TI on NDVI has been increasing during the 20 years. The explanatory credibility of the arable land area share on NDVI shows a fluctuating upward trend.

3.2.2. Interaction Detection of Different Driving Factors

An interaction detector was used to assess the explanatory credibility of different driving factors when they acted together on the spatial distribution of NDVI. The results show that the q statistics after the interaction of any two factors are larger than any single factor, and the interaction types are two-factor enhancement and non-linear enhancement, without weakening or mutual independence, thus indicating that different factors have stronger effects on NDVI when they act together.

Figure 8 shows the interaction detection results of different driving factors. The type of interaction is dominated by two-factor enhancement. In 2000, the explanatory factor of the interactions of annual mean water vapor pressure (WVP) with other natural and anthropogenic activity factors, and the interactions of soil type (Soil) with anthropogenic activity factors on the spatial distribution of NDVI are dominant, with both being greater than 0.62. The strongest explanatory power is annual mean water vapor pressure (WVP) ∩ soil type (Soil) (q = 0.68). In 2005, the explanatory credibility of the interaction between growing season precipitation (Prep2), annual precipitation (Prep1), annual mean water vapor pressure (WVP), and the remaining factors is higher, and the strongest explanatory factor is growing season precipitation (Prep2) ∩ Annual Minimum Temperature (Temp5) (q = 0.68). In 2010, the explanatory credibility of the interaction of growing season temperature (Temp6, Temp7, Temp8), annual mean temperature (Temp1), and soil type (Siol) with anthropogenic activity factors is higher, all greater than 0.7, and the strongest explanatory factor is growing season mean temperature (Temp6) ∩ primary industry output (PI) (q = 0.73). In 2015, annual precipitation (Prep1), growing season precipitation (Prep2), and annual mean water vapor pressure (WVP) interacting with the rest of the factors have higher explanatory credibility, all greater than 0.78, and the strongest explanatory factor is annual precipitation (Prep1) ∩ annual mean air pressure (Press1) (q = 0.82). In 2019, annual precipitation (Prep1) and growing season precipitation (Perp2) interacting with the rest of the factors have higher explanatory credibility, all greater than 0.75, and the strongest explanatory factor is annual precipitation (Perp1) ∩ gross domestic product (GDP) (q = 0.80). All anthropogenic activity factors have significantly higher explanatory credibility after interacting with other factors.

The non-linear enhancement was mainly concentrated in the interaction of single factors that have low q statistics with other factors. The interaction of elevation (Elev), slope (Slope), and population density (Pop) with other factors in all years; the interaction of annual sunshine hours (Sun) with other factors in 2000 and 2005; the interaction of annual maximum temperature (Temp4) with other factors in 2015; the interaction of road density (RD) with natural factors in 2019; and the interaction of temperature factors with air pressure factors in some years were all non-linear enhancement.

3.2.3. Relationship between the Dominant Driving Factors and NDVI

The risk detector can determine whether the differences between the different types of each factor are significant. The results found that, except for slope, there were significant differences between the types of the driving factors, and all of them passed the t-test, indicating that the spatial partitioning of the driving factors in this study is reasonable. Based on the results of the factor detector and interaction detector, it was found that growing season precipitation (Prep2), soil type (Soil), primary industry output (PI), annual mean water vapor pressure (WVP), annual mean relative humidity (RH), annual precipitation (Perp1), annual mean temperature (Temp1), annual mean wind speed (Wind), livestock density (LD), and tertiary industry output (TI) could better explain the spatial distribution of grassland NDVI in Xilingol, and were considered as the dominant drivers of grassland health. Figure 9 shows the statistics of the risk detection results for the above 10 factors, where the horizontal coordinates are the subinterval category amounts for each factor, and the vertical coordinates are the corresponding mean NDVI values.

Among the natural environment factors, the NDVI of grassland tends to increase with increasing growing season precipitation (Perp2), average water vapor pressure (WVP), annual mean relative humidity (RH), and annual precipitation (Perp1). The highest NDVI value is found when the category value is 8, which is the highest value in the range of the factors, indicating that the precipitation-based moisture factor plays a role in promoting grassland health. Additionally, the growing season precipitation was higher than the mean NDVI values of all other factors, which also proves its dominant role in the spatial distribution of NDVI. The trend of NDVI with annual precipitation in 2010 is very different from other years, with NDVI reaching its highest value at the precipitation category, which is 5. The spatial distribution of annual precipitation shows that the precipitation in that year is concentrated in the southern part of the study area where agriculture and livestock are intermingled, whereas the high-NDVI areas in the study area are mainly concentrated in the eastern region, thus causing anomalous changes in the curve. As the mean wind speed and mean annual temperature increased, the mean NDVI values of grassland showed a decreasing trend, indicating that wind speed and temperature had a suppressive effect on the growth of grassland. Among the soil types, the highest mean NDVI values were found for grassland in the soil types Haplic Greyzems, Haplic Luvisols, and Chernozems, but the most widely distributed soil type in the study area was that of Kastanozems, with the average value of NDVI being 0.48.

Among the human activity factors, NDVI showed a fluctuating upward trend as the output value of the primary industry increased, indicating that the mean NDVI value was higher in areas with a developed primary industry. As the value of the tertiary industry increases, NDVI also shows a fluctuating upward trend, but when the value of the tertiary industry reaches the fifth category, NDVI remains stable or shows a decreasing trend, indicating that the development of the tertiary industry on the one hand will promote the healthy growth of the grassland, but may hinder growth when the excessive exploitation of tourism resources causes damage to the health of the grassland. The NDVI increases with the increase in livestock density until the category amount is 4; after which, the NDVI remains stable.

Based on the NDVI high-value areas for the dominant drivers analyzed above, the appropriate range for each driver factor at a high value of NDVI can be inferred.

Table 4. Suitable threshold range of dominant driving factors.

Factor	Range	Unit
Prep2	314.06–412.64	mm
Soil	Greyzems Luvisols Chernozems	−
PI	9.11–12.87	100 million yuan
WVP	6.39–6.65	hpa
RH	61.49–64.35	%
Prep1	413.36–490.91	mm
Temp1	0.64–1.61	°C
Wind	2.51–2.83	m/s
LD	170–263	heads/km2
TI	15.37–16.27	100 million yuan

shows the appropriate threshold ranges for each of the dominant drivers for healthy grassland.

4. Discussion

4.1. Dominant Driving Factors of Grassland Health in Xilingol

The fundamental cause of changes in grassland health in Xilingol over the years has primarily been natural environment changes. The results of the factor detector show that precipitation, soil type, mean water vapor pressure, and mean relative humidity have a high explanatory power for NDVI, with moisture-dominated factors, including precipitation, mean water vapor pressure, and mean relative humidity, dominating the spatial distribution of grassland in the study area [39,40]. The influence of growing season precipitation was significantly greater than that of annual precipitation [41]. The moisture factor played a positive role in NDVI, whereas temperature played a negative role. The results of interaction detection showed that in 2000, 2005, 2015, and 2019, precipitation and mean water vapor pressure had the highest q values after interaction with other factors, but in 2010, temperature had the highest q values after interaction with other factors, whereas the q values of precipitation decreased significantly, and the mean NDVI values were lower in the same years. Combining the actual precipitation and temperature data revealed that precipitation dominated the spatial distribution of NDVI when precipitation was sufficient, and temperature became the dominant driving force affecting the distribution of NDVI when precipitation was insufficient. When water and heat conditions are favorable, soil type is less influential, whereas when water and heat conditions are poor, suitable soil types are more favorable for grassland growth.

4.2. Impact of Human Activities on Grassland Health

Among the anthropogenic disturbances, the output value of the primary industry has the strongest explanatory power for the spatial distribution of NDVI. Livestock density is consistent with the change in q value of primary industry output, indicating that grazing is the main pillar of primary industry in the study area. As a typical pastoral area, the population distribution in the study area is relatively sparse, so the influence of population density on the distribution of grassland NDVI is relatively weak. The opposite change in livestock density and the q value of the output value of the secondary industry also indicates that grazing accounts for a smaller share in areas where the industry is more developed. The q value of tertiary industry output is on the increase and is gradually replacing primary industry output as the most important driver of the human activity factors. The q value of livestock density is greater than 0.5 in 2005, 2010, and 2019, and the main contribution of this value comes from the density of sheep in the study area. Industrial restructuring had a large impact on the spatial distribution of NDVI. It gradually changed from primary-industry-dominated to primary-and-secondary-industry-jointly-influenced, and finally to tertiary-industry-dominated in 2000–2019.

4.3. The Differences among Interactions of the Natrual and Anthropogenic Factors

The results of the interaction detector are divided into three parts: the first part is the interaction among natural environmental factors, the second part is the interaction between natural environmental factors and the anthropogenic disturbance, and the third part is the interaction among anthropogenic factors. Over the years, the q statistics were higher after the interaction of the first two components, whereas the q statistics of the third component were consistently lower and less varied. This is because the spatial scales of the anthropogenic disturbance factors used in this study are mostly at the banner and county scales, and the administrative areas of Inner Mongolia are quite large, thus, to some extent, blurring the differences between the different spatial scales, and resulting in a small difference in q statistics between the anthropogenic factors.

Further work can take this as a starting point to improve the spatial resolution of anthropogenic factors using remote sensing data, and, thus, evaluate the impact of anthropogenic activities on changes in grassland health in a more detailed way.

5. Conclusions

This study uses trend analysis methods and geodetector methods to examine the spatial and temporal evolutionary patterns of grasslands and their drivers from 2000–2019 in Xilingol, leading to the following conclusions.

The NDVI of Xilingol grassland shows a spatial distribution pattern of a high value in the east and low value in the west, and a mean NDVI value of 0.39. The temporal pattern shows a fluctuating upward trend, with a monthly increasing rate of 0.0005/month and an annual increasing rate of 0.0028/a. Spatially, it also shows a trend of recovery, with 68.06% of the grassland in the region at an average recovery level.

Using geodetector methods, the evaluation of the dominant drivers can be completed quickly and effectively. The results show that moisture-dependent natural environmental factors, soil type, and primary industry output are the dominant drivers of grassland health in Xilingol. The essential causes of changes in grassland health are changes in natural climatic conditions, with moisture factors playing a facilitating role, and temperature and wind speed factors playing a suppressing role. Among the anthropogenic activities’ disturbances, the adjustment of industrial structure has the greatest influence on the spatial distribution of NDVI. As a typical grazing area, livestock density also plays an important role, but due to the limited spatial scale, the explanatory credibility varies less between the factors. Different types of drivers can increase the explanatory credibility when they act together. When the dominant driver is within the range of multi-year average thresholds, it can contribute to the sustained recovery of grasslands.