Spatial-Temporal Characteristics and Driving Forces of Aboveground Biomass in Desert Steppes of Inner Mongolia, China in the Past 20 Years

Abstract

The desert steppe serves as a transitional zone between grasslands and deserts, and long-term monitoring of aboveground biomass (AGB) in the desert steppe is essential for understanding grassland changes. While AGB observation techniques based on multisource remote-sensing data and machine-learning algorithms have been widely applied, research on monitoring methods specifically for the desert steppe remains limited. In this study, we focused on the desert steppe of Inner Mongolia, China, as the study area and used field sampling data, MODIS data, MODIS-based vegetation indices (VI), and environmental factors (topography, climate, and soil) to compare the performance of four commonly used machine-learning algorithms: multiple linear regression (MLR), partial least-squares regression (PLS), random forest (RF), and support vector machine (SVM) in AGB estimation. Based on the optimal model, the spatial–temporal characteristics of AGB from 2000 to 2020 were calculated, and the driving forces of climate change and human activities on AGB changes were quantitatively analyzed using the random forest algorithm. The results are as follows: (1) RF demonstrated outstanding performance in terms of prediction accuracy and model robustness, making it suitable for AGB estimation in the desert steppe of Inner Mongolia; (2) VI contributed the most to the model, and no significant difference was found between soil-adjusted VIs and traditional VIs. Elevation, slope, precipitation, and temperature all had positive effects on the model; (3) from 2000 to 2020, the multiyear average AGB in the study area was 58.34 g/m², exhibiting a gradually increasing distribution pattern from the inner region to the outer region (from north to south); (4) from 2000 to 2020, the proportions of grassland with AGB slightly and significantly increasing trend in the study area were 87.08% and 5.13%, respectively, while the proportions of grassland with AGB slightly and significantly decreasing trend were 7.76% and 0.05%, respectively; and (5) over the past 20 years, climate change, particularly precipitation, has been the primary driving force behind AGB changes of the study area. This research holds reference value for improving desert steppe monitoring capabilities and the rational planning of grassland resources.

1. Introduction

The Desert steppe in Inner Mongolia is an important component of the temperate grasslands in Central Asia and plays a crucial role in regulating the regional climate, preserving soil and water resources, and providing habitats for numerous species [1,2,3]. However, climate change and inappropriate land use have put these areas under severe pressure [4]. Aboveground biomass (AGB) is a valuable indicator for assessing terrestrial ecosystem health and understanding its spatial–temporal characteristics and underlying factors is critical for effectively conserving and managing the desert steppe ecosystem [5,6]. Field observation and remote sensing are the two main methods for grassland AGB estimation. Compared with traditional field-based methods, remote-sensing-based techniques have emerged as a valuable tool for efficiently and accurately tracking AGB changes over large areas [5,7]. In recent decades, a multitude of remote-sensing-based AGB observation technologies has been developed and found to be suitable for vegetation observation problems at various scales [7,8,9]. However, constructing a fitting model that accurately captures changes in AGB depends heavily on the selection of both feature variables and appropriate modeling algorithms.

Feature selection is a critical step in remote-sensing-based AGB monitoring [10]. Currently, the main variables used to retrieve vegetation AGB are derived from spectral reflectance and vegetation indices (VI) [11,12]. As VIs are calculated based on the reflectance of multiple spectral bands, they exhibit high sensitivity to surface features, including vegetation, making them the preferred variables for biomass observation [13]. In addition to traditional VIs, such as the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI), modified indices, such as the soil-adjusted vegetation index (SAVI), the modified soil-adjusted vegetation index (MSAVI), the transformed soil-adjusted vegetation index (TSAVI), the optimized soil-adjusted vegetation index (OSAVI), and the adjusted transformed soil-adjusted vegetation index (ATSAVI), have been proposed and widely used in vegetation observation [14,15,16,17,18,19,20]. However, studies showed that some soil-adjusted VIs may not perform better than traditional indices in sparsely vegetated dry grasslands [21,22]. In addition, relying solely on VI in biomass monitoring models may result in regional differences in form and accuracy. To address this issue, some researchers integrated environmental variables such as geography, topography, and climate into single-factor models and developed multivariate regression models that significantly improve model-fitting accuracy [23,24,25,26]. Nevertheless, further research is necessary to determine which variables are suitable for estimating AGB in the desert steppe.

Modeling algorithm selection is a crucial aspect of the biomass estimation process [27]. Currently, there are numerous algorithms available for retrieving AGB, which can be broadly categorized into parametric and nonparametric algorithms based on whether they require assumptions about the model structure. Studies used VIs as the observed variable and constructed single-variable models, such as linear, exponential, and power to fit vegetation biomass [28,29,30]. However, the accuracy of these simple models is prone to be affected by soil and spectral saturation effects, particularly in heterogeneous regions [31,32]. To improve the generalization ability and accuracy, studies emphasized the establishment of multiple linear regression (MLR) models by increasing the feature variables [33,34,35]. In recent years, MLR models incorporating dimensionality reduction techniques such as the partial least squares regression (PLS) have also gained popularity [36,37]. However, as vegetation growth is a complex process, parametric models often lack the ability to fit this complex relationship, resulting in a lack of ecological significance in the entire operation process. In contrast, nonparametric models do not require assumptions about the form of the objective function and can thus fit complex relationships, making them the most popular machine-learning algorithms for estimating AGB [5,13]. Some studies showed that these algorithms, such as the random forest (RF), and support vector machine (SVM) perform better than traditional parametric models in estimating grassland AGB [38,39]. However, no indepth comparative study of AGB in Inner Mongolia’s desert steppe has been conducted. The desert steppe in Inner Mongolia, China, is characterized by sparse vegetation and complex environmental conditions, leading to uncertainty in the performance of these algorithms in fitting vegetation AGB. Therefore, further research is needed to evaluate the performance of these models in estimating AGB in this region.

Human activities and climate change are widely considered to be the main factors affecting grassland health, but quantifying their dominance is challenging [40,41,42,43]. Some studies used linear regression coefficients and variance ratios to assess the impact of external drivers on grasslands [43,44,45]. However, these methods are prone to overlooking ecological processes and cannot reflect complex ecological relationships [46]. Other studies attempted to establish quantitative analysis conditions from an ecological perspective by using ecological process models to explain driving factors [47,48,49]. However, this approach makes it difficult to provide data references for some specific driving forces, resulting in a lack of indepth analysis. Furthermore, some recent studies attempted to use nonparametric models such as RF [50] and boosted regression trees [51] to analyze the impacts of external driving forces on land changes. Similarly, we believe that this approach is also applicable to the topic of our research. The relationship between vegetation and external factors is complex and diverse, making nonparametric algorithms more advantageous in simulating this complex relationship. In addition, selecting appropriate algorithms and building grassland AGB inversion models that are suitable for the characteristics of the study area is key to accurately reflecting the spatial–temporal characteristics and driving forces of desert steppe in Inner Mongolia.

To fill the knowledge gaps mentioned above, this study focused on the desert steppe of Inner Mongolia, China, and carried out the following studies: (1) compared the performance of four commonly used machine-learning algorithms (MLR, PLS, RF, and SVM) on grassland AGB estimation, and generated a suitable modeling method; (2) extracted the spatial–temporal feature of grassland AGB in past 20 years (2000–2020); and (3) analyzed the driving forces of grassland AGB change in past 20 years.

2. Materials and Methods

2.1. Study Area

This research focuses on five regions (counties) in central Inner Mongolia, namely Sunite East Banner (SE), Sunite West Banner (SW), Erenhot City (ER), Siziwang Banner (SZW), and Damao Banner (DM) (Figure 1). The region is predominantly characterized by desert grasslands, covering an approximate area of 56,200 km², and represents the main body of the desert steppe in Inner Mongolia. The study area spans from 37°N to 46°N and 105°E to 115°E, with an elevation range of 800–2100 m and an average elevation of 1100 m. The climate in this area is characterized by a continental arid climate with an annual precipitation of 100–200 mm and an average annual temperature of 2–5 °C [2]. The grassland vegetation in this region is typified by a low diversity of species and a simple community structure. The drought-tolerant bunchgrasses are the dominant grass species, often accompanied by numerous drought-tolerant shrubs. The common herbaceous plants include Stipa klemenzii, Stipa breviflora, Artemisia frigida, and Cleistogenes songorica [2,52].

2.2. Data Used

2.2.1. Field Measurement Data

During the peak growing season (July and August) from 2018 to 2020, 97 field plots from various terrains and vegetation coverage were collected (Figure 1). The plot size was 500 m × 500 m, with each corner aligned as closely as possible with the corresponding MODIS pixel. Within each plot, 3–5 samples of 1 m × 1 m were selected along the diagonal, and the average value of these sample points was used as the AGB for that plot (Figure 1). The distance between adjacent plots was no less than 1 km, and the recorded information included fresh biomass, latitude, longitude, and elevation.

2.2.2. Remote-Sensing Data

• Multispectral Bands and Vegetation Indices

The MOD09GA product, which provides surface reflectance data, was acquired from the Google Earth Engine (GEE) platform. The time period covers from July 2000 to September 2020. The MOD09GA product records the surface spectral reflectance of MODIS bands 1 through 7 at a spatial resolution of 500 m, corrected for atmospheric conditions such as gases, aerosols, and Rayleigh scattering [53]. Specifically, the bands included in the MOD09GA are red (band 1), near infrared (band 2), blue (band 3), green (band 4), shortwave infrared 1 (band 6), and shortwave infrared 2 (band 7). In this study, in addition to selecting the above-mentioned six bands, seven MODIS-based VIs were considered. These VIs include two traditional VIs (NDVI, EVI) and five soil-adjusted vegetation indices (SAVI, MSAVI, TSAVI, OSAVI, and ATSAVI). The specific calculation formulas for these indices are given in

Table 1. Selected vegetation indices.

Vegetation Indices Type	Description
Traditional vegetation indices	$NDVI=\frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}}$
	$EVI=2.5\times \frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+6\times {\rho }_{RED}-7.5\times {\rho }_{BLUE}+1}$
Soil-adjusted vegetation indices	$SAVI=\left(1+L\right)\times \frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}+L}$
	$MSAVI=0.5\times [(2\times {\rho }_{NIR}+1)-\sqrt{{\left(2\times {\rho }_{NIR}+1\right)}^2-8\times \left({\rho }_{NIR}-{\rho }_{RED}\right)}]$
	$TSAVI=\frac{a\times ({\rho }_{NIR}-a\times {\rho }_{RED}-b)}{a\times {\rho }_{NIR}+{\rho }_{RED}+a\times b}$
	$OSAVI=1.16\times \frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}+0.16}$
	$ATSAVI=\frac{a\times ({\rho }_{NIR}-a\times {\rho }_{RED}-b)}{a\times {\rho }_{NIR}+{\rho }_{RED}-a\times b+0.08\times \left(1+a^2\right)}$

.

As some of these indices require correction parameters, 30 barren sample sites were selected in the study area to establish a soil line using the reflectance of soil in the near-infrared and red spectral bands, and obtain the correction parameters a and b.

• Terrain data

The terrain variables, elevation, slope, and aspect were from ASTER GDEM version 2, which is a global DEM generated through a collaboration between NASA and Japan’s Ministry of Economy, Trade, and Industry (METI). Despite variations in accuracy depending on terrain and other factors, the ASTER GDEM v2 is generally considered a high-quality DEM with a root mean square error (RMSE) of less than 10 m [54]. In this study, we applied the GEE platform to complete preprocessing operations such as projection transformation, resampling, and clipping.

• Meteorological data

Annual total precipitation data were obtained from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS), which is a quasi-global rainfall dataset spanning over 30 years. CHIRPS combines 0.05° resolution satellite imagery with in situ station data to generate gridded rainfall time series suitable for trend analysis and seasonal drought monitoring [55].

Annual mean temperature and annual mean wind data were sourced from the ERA5-Land dataset, the latest global meteorological reanalysis dataset released by the European Centre for Medium-Range Weather Forecasts (ECMWF). The dataset provides continuous meteorological records since 1950 at a pixel resolution of up to 9 km and encompasses all conventional meteorological variables [56].

Annual total evapotranspiration data were extracted from the MOD16A2 Version 6 product, which is an 8-day composite dataset produced at a 500 m pixel resolution. The algorithm utilized for the MOD16 data product collection is based on the Penman-Monteith equation and incorporates inputs of daily meteorological reanalysis data, along with MODIS remotely sensed data products such as vegetation property dynamics, albedo, and land cover [57]. In this study, we employed the GEE platform to preprocess the aforementioned factors, including projection conversion, resampling, cropping, and unit conversion. The datasets cover the time period from 2000 to 2020.

• Soil data

The soil factors, including soil bulk density, soil pH, clay content, and sand content at 10 cm depth, were obtained from the OpenLandMap soil dataset, which is an environmental data portal that provides global-scale data on land cover, vegetation, soil, climate, terrain, and other related information to the public [58]. The soil dataset in OpenLandMap is established through the integration of global soil profile data points and machine-learning algorithms and provides multiple soil properties data at six standard depths (0, 10, 30, 60, 100, and 200 cm) with a resolution of 250 m [59]. In this study, we utilized the GEE platform to extract each factor and performed preprocessing operations such as projection transformation, resampling, and cropping.

• Population data

The population data was sourced from the WorldPop project, which is an interdisciplinary research group dedicated to improving the evidence base for spatial demographics and utilizing the data for health and development applications [60,61]. The project takes recent population counts from censuses and breaks them down into grid cells of approximately 100 m × 100 m using machine-learning techniques. These techniques exploit the relationships between population densities and different geospatial covariate layers. In this study, we employed the GEE platform to carry out preprocessing tasks. The datasets used encompass the period from 2000 to 2020.

2.2.3. Statistical Data

The statistical data involved in this study include annual livestock quantities, agricultural gross domestic product (GDP), industrial GDP, and consumer GDP for each county from 2000 to 2020. These data were obtained from local statistical yearbooks. To capture the multiyear trends of these statistical variables, we employed the least squares method and assigned the resulting trends to the corresponding administrative boundary polygons.

2.3. Analysis Methods

The technical route of this study can be divided into four parts: data preprocessing, optimal model construction, extraction of AGB spatial–temporal features, and analysis of driving forces (Figure 2). As mentioned in the previous section, the data used in this study mainly consist of field-sample data, remote-sensing data, and statistical data. First, extract the raster values corresponding to the coordinates of the sample plots for the year 2020 to build a sample dataset. Then, split the dataset into training and testing datasets using a 70:30 ratio. Next, the training dataset is used in combination with leave-one-out cross-validation (LOOCV) to train four different models: MLR, PLSR, RF, and SVM. The performance of these models is compared on the testing dataset, and the best model and variable combination are selected. Then, using the selected model, a spatial–temporal dataset of AGB in the desert steppe of the study area is generated for the period from 2000 to 2020. Following this, change detection and trend analysis methods are employed to characterize the changes in desert steppe AGB over the past 20 years. Lastly, the RF algorithm is utilized to quantitatively analyze the potential driving forces behind the AGB changes.

2.3.1. Modeling Algorithms

MLR is the most used parametric statistical algorithm for relating linear relationships between various predictors and independent variables [62]. Stepwise regression is a statistical approach used to build an MLR model by selecting the most significant predictors from a larger set. The algorithm adds or removes variables from the model in a stepwise manner based on predetermined criteria until the optimal subset of variables is achieved. The final model consists of the most significant predictors, and the coefficients represent the strength and direction of the relationship between predictors and the outcome variable [63].

PLS is a statistical method that is similar to principal components regression in that it finds a linear regression model by projecting the predicted variables and the observable variables to a new space rather than finding hyperplanes of maximum variance between the response and independent variables [64]. This algorithm is especially useful when the predictor matrix contains more variables than observations and there is multicollinearity among the variables.

RF is a powerful machine-learning algorithm used for regression analysis that falls under the category of ensemble learning algorithms. The algorithm works by creating multiple decision trees on randomly sampled subsets of the training data and averaging their predictions to obtain the final output [65]. Generally, RF involves two important hyperparameters: the number of decision trees (n_estimators) and the maximum number of features considered when splitting a node in each decision tree (max_features) [65]. The n_estimators refer to the number of decision trees used to construct the random forest, and the max_features refer to the maximum number of features considered when splitting a node in each decision tree.

SVM is a popular machine-learning technique for predicting continuous output values. The algorithm aims to find the best hyperplane that minimizes the error between the predicted and true output values, with a penalty for the weight magnitude [66]. It can handle nonlinear data using kernel functions such as the radial basis function (RBF), which maps data into higher dimensional feature spaces. The SVM-RBF algorithm is capable of modeling nonlinear relationships between input and output variables and is robust to noise and outliers. Two hyperparameters need to be adjusted: penalty factor (c) and kernel parameter (gamma). The C parameter controls the trade-off between achieving a low error on the training data and minimizing the weights’ mean, while the gamma parameter controls the distance of influence of a single training point.

Feature selection has always been an important step in model training. In this study, the recursive feature elimination (RFE) method was employed to iteratively select feature variables for both RF and SVM. RFE is a feature selection technique that aims to iteratively eliminate features with lower predictive contributions to the model. In this method, the model is repeatedly trained, and less informative features are removed until a subset of the most predictive features is obtained [67].

2.3.2. Model Assessment

In this study, three validation standards, namely R², RMSE, and RMSEr, were selected to evaluate model accuracy, with the following specific formulas:

$$R^2=1-\frac{\sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{\sum _{i=1}^n{(y_i-\stackrel{-}{y})}^2}$$

(1)

$$RMSE=\sqrt{\sum _{i=1}^n\frac{{({\widehat{y}}_i-y_i)}^2}{n}}$$

(2)

$$RMSEr=\frac{RMSE}{\stackrel{-}{y}}\times 100$$

(3)

Here, ${\widehat{y}}_i$ represents the measured AGB value of the i-th sample, $y_i$ represents the predicted AGB value of the i-th sample, $\stackrel{-}{y}$ represents the average AGB value of the samples, and n represents the total number of samples.

2.3.3. Change Point Detection

The Pettitt test was used in this study to extract the abrupt change point of grassland AGB time-series data. The Pettitt test is a nonparametric test based on Mann–Whitney two independent samples and its hypothesis is that H₀ states that if two samples have no difference, then their central positions will be the same, while H₁ states that if there is a difference, then their central positions will not be the same [68].

The nonparametric statistic is defined as follows:

$$\widehat{U}=\mathrm{m}\mathrm{a}\mathrm{x} \left|U_k\right|$$

(4)

$$U_k=2\sum _{i=1}^k{r}_i-k(n+1),k=1,2,\dots ,n$$

(5)

Here, $\widehat{U}$ is the test statistic, and its value is equal to the maximum absolute value of the vector $U_k$; $r_i$ is the order of the series data, and k is the possible change year.

The formula for the p-value of Pettitt’s test:

$$p=2\mathrm{e}\mathrm{x}\mathrm{p}(\frac{-6{\widehat{U}}^2}{n^3+n^2})$$

(6)

Here, if p is less than a given significance level β, reject the null hypothesis, namely series data has a change point; On the contrary, if p is greater than a given significance level β, accept the null hypothesis, namely the series data without change point.

2.3.4. Trend Analysis

In this study, the grassland AGB changing trend was analyzed using Sen’s slope method. The Mann–Kendall test was employed to test the significance of the sequence trend. The Mann–Kendall test was preferred as it does not make any assumptions about the distribution of the sequence and is robust to outliers [69].

Sen’s slope formula is:

$$\beta =Median(\frac{{AGB}_j-{AGB}_i}{j-i})$$

(7)

Here, Median(x) represents the median function, and β is the Sen’s slope value. A positive value of β indicates an increasing trend of AGB, while a negative value of β indicates a decreasing trend.

The Mann–Kendall test involves testing two hypotheses: H₀ assumes that the sequence data are randomly arranged and there is no significant change trend, while H₁ assumes that the sequence data exhibit an upward or downward trend. The formula for the test statistic S is:

$$S=\sum _{i=1}^{n-1}\sum _{j=i+1}^{n}sgn({AGB}_j-{AGB}_i)$$

(8)

$$sgn\left({AGB}_j-{AGB}_i\right)=\left\{\begin{array}{c}1\\ 0\\ -1\end{array}\right.\begin{array}{c}{AGB}_j-{AGB}_i>0\\ {AGB}_j-{AGB}_i=0\\ {AGB}_j-{AGB}_i<0\end{array}$$

(9)

Here, sgn(x) is the sign function that returns 1, −1, or 0, depending on whether x is positive, negative, or zero, respectively.

The values of Sen’s slope (β) and the Mann–Kendall test statistic (p) were used to determine the levels of the AGB change trend, which included significantly increasing, slightly increasing, significantly decreasing, and slightly decreasing (

Table 2. Levels of AGB change the trend of grassland.

AGB Change Trend and Test Results	Levels of AGB Change Trend
ꞵ > 0, p < 0.05	significantly increasing
ꞵ > 0, p > 0.05	slightly increasing
ꞵ < 0, p < 0.05	significantly decreasing
ꞵ < 0, p > 0.05	slightly decreasing

).

2.3.5. Driving Factors Analysis

The RF algorithm was employed in this study to analyze the driving factors behind changes in grassland AGB quantitatively. The following steps were taken in the analysis: (1) nine potential driving forces, from climate and human activity categories, were selected as independent variables to investigate their impact on grassland AGB changes in the study area (

Table 3. Selected potential driving forces of AGB changes.

Factor Type	Factor Name	Description	No.	Unit
Climate change	Precipitation trends	Annual precipitation trends from 2000 to 2020	D1	mm
	Temperature trends	Annual temperature trends from 2000 to 2020	D2	°C
	Evapotranspiration trends	Annual evapotranspiration trends from 2000 to 2020	D3	mm
	Wind trends	Annual wind trends from 2000 to 2020	D4	m
Human activity	Population density trends	Annual population density trends from 2000 to 2020	D5	%
	Livestock growth rate	Annual average livestock growth rate	D6	%
	Primary industry GDP growth rate	Averaged annual growth rate of primary industry GDP from 2000 to 2020	D7	%
	Secondary industry GDP growth rate	Averaged annual growth rate of secondary industry GDP from 2000 to 2020	D8	%
	Tertiary industry GDP growth rate	Averaged annual growth rate of tertiary industry GDP from 2000 to 2020	D9	%

); (2) 300 sample points were randomly selected, and the corresponding independent variable and the AGB change trend (dependent variable) values were extracted for each point; (3) the sample points were divided into training and testing datasets in a 70:30 ratio, and an RF regression model was optimized using the training dataset and a 10-fold cross-validation method; (4) the feature importance was extracted to quantitatively evaluate the driving factors.

RF provides two metrics for measuring the importance of features: mean decrease impurity (MDI) and mean decrease accuracy (MDA). The MDI metric measures the number of times a feature is selected in the decision tree and its contribution to the impurity reduction, which can evaluate the contribution of a feature to classification. On the other hand, the MDA metric calculates the change in model accuracy on a test set by randomly shuffling the values of a feature, where a greater decrease in accuracy indicates a higher feature importance. It should be noted that although both metrics can be used to measure feature importance, they are not always consistent.

3. Results

3.1. The Performance of the Selected Algorithms

Figure 3 presents the relationship between algorithm performances and the number of variables (components). The cross-validation accuracy of these algorithms on the training dataset was unstable, indicating that an excessive number of variables does not improve the performance of the model. MLR retained eight variables for optimal performance after stepwise feature testing, indicating the presence of multicollinearity among the selected features. The optimal performance of PLS was achieved by retaining only one principal component. After evaluating each feature using the RFE method, RF and SVM chose 17 and 22 effective variables, respectively, to achieve the best performance.

Table 4. Selected effective variables for models.

Algorithms	Selected Variables	No. of Variables
MLR	MIR, SWIR, MSAVI, ATSAVI, soil bulk density, clay content, precipitation, temperature	8
PLS	All	22
RF	Blue, green, red, NIR, MIR, SWIR, NDVI, EVI, SAVI, TSAVI, OSAVI, MSAVI, ATSAVI, elevation, slope, precipitation, temperature	17
SVM	All	22

lists the detailed features selected by the modeling algorithms. MLR selected the feature variables MIR, SWIR, MSAVI, ATSAVI, soil bulk density, clay content, precipitation, and temperature to construct the AGB regression model, expressed as AGB = 35.14 × MIR − 33.69 × SWIR − 200.18 × MSAVI + 229.79 × ATSAVI + 18.44 × soil bulk density − 21.27 × clay content + 22.94 × precipitation − 12.82 × temperature + 73.41. PLS retained all variables since it integrates dimensionality reduction techniques. RF identified the effective variables as blue, green, red, NIR, MIR, SWIR, NDVI, EVI, SAVI, TSAVI, OSAVI, MSAVI, ATSAVI, elevation, slope, precipitation, and temperature, while SVM selected all variables.

Figure 4 displays the model performance of RF and SVM under different hyperparameter combinations at the optimal variable conditions. For the RF, the model achieved high accuracy when the n_estimators was set between 500 and 700, and the best model performance was obtained when the n_estimators and max_features were set to 700 and 4, respectively. As for the SVM, the model achieved high prediction accuracy when gamma was set between 0.02 and 0.05, and the best model performance was achieved when gamma and cost were set to 0.03 and 8, respectively.

Figure 5 compares the prediction accuracy of selected algorithms on the training and testing datasets under effective feature and hyperparameter conditions. The results indicate that nonparametric models outperformed linear models significantly. On the testing dataset, RF achieved the highest prediction accuracy (R² = 0.66, RMSE = 39.85, RMSEr = 39.12), followed by SVM (R² = 0.60, RMSE = 52.16, RMSEr = 51.21), PLS (R² = 0.59, RMSE = 44.47, RMSEr = 43.67), and MLR (R² = 0.52, RMSE = 48.85, RMSEr = 47.96). These findings suggest that RF has the best predictive ability for unseen data, which is consistent with the inversion requirements of AGB in Inner Mongolia’s desert steppe.

3.2. Spatial–Temporal Feature of Desert Steppe AGB in Inner Mongolia during the Past 20 Years

As seen in Figure 6, from 2000 to 2020, the vegetation cover of the desert steppe in the study area exhibited a gradually increasing distribution pattern from north to south, with an average AGB of 58.34 g/m². Among them, the SZW had the highest average AGB of 68.91 g/m², followed by SZ (60.90 g/m²), DM (56.52 g/m²), SY (49.82 g/m²), and ER (35.32 g/m²). Within the entire study area, the AGB was highest in the southern parts of SZW and DM, exceeding 120 g/m², while the vegetation in most parts of ER, the northern part of SZW, and the western part of SE was sparse, with an AGB ranging from 20–50 g/m².

Figure 7a,b demonstrates that while AGB changed in most of the study area, only a small number of regions experienced significant changes (p < 0.05). This suggests that over the past 20 years, the AGB changes in the study area have been comparatively stable. Figure 7c shows that the areas with significant changes over the past 20 years were mainly distributed in the southern SZW and central DM, and the years with significant changes mainly occurred in 2010, 2011, and 2012.

Figure 8 displays the spatiotemporal variation of AGB in the study area from 2000 to 2020. Overall, AGB in desert steppe increased with an average annual growth rate of 1.33 g/m² during the past two decades (Figure 8a,d). According to the significance testing results (Figure 8b), grassland with a significant increase in AGB accounted for 5.13%, primarily found in the central and southern SZW and central DM, while grassland with a slight increase in AGB accounted for 87.08%, found throughout the study area. The proportion of grassland with a significant decrease in AGB was less than 0.5% and the proportion of grassland with a slight decrease was only 7.76%, with the majority concentrated in the eastern SW, central SE, and northeastern ER (Figure 8).

At the county level, SZW had the most significant restoration, with an annual increase in grassland AGB of 1.84 g/m², and the proportion of grassland with a significant increase in AGB was 9.95%, while the proportion of grassland with a slight increase in AGB was 87.17%. DM followed, with an annual increase in grassland AGB of 1.44 g/m², with 94.21% of grassland showing an increase in AGB. In SE, the annual increase in grassland AGB was 1.30 g/m² and the proportion of grassland with an increase in AGB was 92.38%. In ER, the annual increase in grassland AGB was 0.89 g/m², with 90.09% of grassland showing an increase in AGB. In SW, the annual increase in grassland AGB was 0.84 g/m² and the proportion of grassland with an increase in AGB was 85.57% (Figure 8c,d).

3.3. Driving Factors for Desert Steppe AGB in Inner Mongolia during the Past 20 Years

Figure 9 displays the drivers behind the desert steppe AGB changes in the study area from 2000 to 2020. Both indicators provide similar rankings. In the MDA ranking, precipitation trends were the most important, followed by evapotranspiration trends > wind trends > temperature trends > livestock growth rate > secondary industry GDP growth rate > population growth rate > tertiary industry GDP growth rate > primary industry GDP growth rate. Similarly, in the MDI, precipitation trends were the most important, followed by temperature trends > wind trends > population growth rate > livestock growth rate > secondary industry GDP growth rate > tertiary industry GDP growth rate > primary industry GDP growth rate. Combining the two indicators, it can be concluded that climate change is the primary factor influencing the changes in AGB in the study area’s desert steppe and human activities have a much smaller impact compared to climate change.

4. Discussion

4.1. Modeling Algorithm and Feature Selection for Estimating AGB in the Desert Steppe of Inner Mongolia

The desert steppe of Inner Mongolian is characterized by sparse vegetation and simple structure, which poses challenges for reflecting the spatial and temporal characteristics of AGB in the region. In this study, the performance of four popular machine-learning algorithms, MLR, PLS, RF, and SVM, were compared for AGB estimation. From the perspective of prediction accuracy, RF performed the best. Moreover, RF also outperformed other algorithms in terms of noise resistance. This fully demonstrates that RF is suitable for AGB inversion tasks in the Inner Mongolian desert steppe. Similarly, some other studies have observed the same results [24,70]. RF predicts by combining multiple decision-tree models, and due to the ensemble effect of the model, RF has high accuracy and good handling capabilities for high-dimensional and missing data, as well as advantages in interpretability and noise resistance [65]. This makes it one of the commonly used algorithms in machine learning. However, these nonparametric models also have disadvantages, such as low interpretability and reduced model portability due to the influence of spatial heterogeneity [27]. This uncertainty raises the issue of whether the trained RF model is applicable to other desert steppe areas outside the study area, which requires further research for generalization. In comparison, parametric models such as MLR and PLS have advantages in model interpretability, flexibility, and ease of operation. However, based on the results of this study, these models are easily affected by multicollinearity, leading to the exclusion of some key variables from the modeling process. In conclusion, RF is the optimal model for obtaining AGB in the Inner Mongolian desert steppe.

Among the tested variables, VIs contributed the most to AGB inversion in the desert steppe, particularly in nonparametric models, where these indices were all selected as effective variables. However, in the MLR model, only MSAVI and ATSAVI were chosen, indicating the presence of multicollinearity among these VIs. A Spearman correlation test revealed significant correlations among the chosen vegetation indices (Figure A1). This suggests that soil-adjusted VIs do not significantly enhance the spatial heterogeneity representation of AGB in the desert steppe at MODIS spatial resolution scale. Similarly, Baghi and Oldeland found comparable results when comparing soil-adjusted VIs and standard VIs for AGB inversion in Northern Iran’s semiarid grasslands using MODIS data [22]. Ren et al. emphasized that the performance of soil-adjusted VIs based on MODIS data was inferior to traditional VIs [21]. This could be attributed to the image resolution, as most of these indices were developed using Landsat imagery, which may affect their performance on coarser spatial-resolution images. On the other hand, soil–vegetation indices might lose their ability to describe the green vegetation information in areas with low green-vegetation cover [21]. In terms of environmental variables, elevation, slope, precipitation, and temperature all contributed positively to AGB inversion in the desert steppe. This implies that the selected samples exhibit spatial differences in these features. The selection of soil factors should be further determined based on the modeling algorithm.

4.2. Spatial–Temporal Feature and Driving Forces of AGB Changes in the Desert Steppe of Inner Mongolia

Over the past 20 years, the multiyear average AGB of the desert steppe in the study area was 58.34 g/m² and it exhibited a gradually decreasing trend from the outer region towards the inner region. This pattern is consistent with the distribution characteristics of annual precipitation in this area. Similarly, some previous studies have also observed the same distribution patterns [24,71]. This indicates that precipitation is the primary factor determining the vegetation distribution characteristics in the region. In terms of the change trend, the AGB exhibited an increasing tendency, with an average annual growth of 1.33 g/m², and most regions have not experienced abrupt shifts. That is, the desert steppe in Inner Mongolia has undergone a stable recovery during these 20 years. Human activities and climate change have long been considered the primary factors driving changes in grasslands [44]. However, quantitatively analyzing their impacts on grasslands presents challenges. This study considered that the interaction mechanisms between external factors and grassland vegetation are complex and diverse. Therefore, employing nonparametric algorithms to address these intricate relationships is more consistent with the actual situation. Based on the feature importance ranking from the RF model, climate change emerges as the dominant factor influencing AGB changes in the grasslands of the study area, particularly the contributions of precipitation and evapotranspiration. During the past 20 years, most of the regions in the study area have experienced an upward trend in precipitation, with annual precipitation surpassing annual evapotranspiration. This implies that the water required for vegetation growth has been adequately replenished, promoting vegetation productivity. The desert steppe is primarily dominated by xerophytic, perennial bunchgrasses, but during wet years, annual hydrophilic vegetation can easily form a dominant layer, positively affecting community biomass [46]. Moreover, results show that temperature and wind speed are also major factors causing AGB changes. Nevertheless, the variations in these two factors over the past 20 years have been relatively minor, implying that their associations with AGB are more prominently manifested in spatial distribution patterns. Compared with climatic factors, human activities exert a less pronounced impact on AGB within the desert steppe of the study area. However, it is worth noting that quantifying social statistical indicators such as GDP and livestock numbers for spatial use, particularly on smaller scales, may be challenging and could lead to uncertainties in the results. Therefore, the study of spatial transformation methods for socioeconomic indicators is a key issue for future quantitative analysis of driving forces.

5. Conclusions

This study employed field-survey data and multisource remote-sensing data to compare the performance of four commonly used machine-learning algorithms, including linear and nonparametric methods (MLR, PLS, RF, and SVM), in estimating desert steppe AGB in Inner Mongolia, China. The optimal model was used to analyze the spatiotemporal characteristics of AGB from 2000 to 2020, and the driving forces, including climate change and human activities, on AGB change were quantitatively assessed based on the RF algorithm. The main conclusions are as follows: (1) RF outperformed other algorithms regarding prediction accuracy and model robustness, rendering it suitable for AGB estimation in the desert steppe of Inner Mongolia; (2) among the tested variables, VI contributed the most to the model and no significant difference was found between soil-adjusted VIs and traditional VIs. Elevation, slope, precipitation, and the temperature had positive effects on model construction; (3) from 2000 to 2020, the multiyear average AGB in the study area was 58.34 g/m² and it exhibited a gradually decreasing trend from the outer region towards the inner region; (4) from 2000 to 2020, grassland with a slight increase in AGB accounted for 87.08%, significantly increased accounted for 5.13%, significantly decreased accounted for only 0.05%, and slightly decreased accounted for 7.76% in the study area; (5) over the past two decades, climate change, particularly precipitation, is the dominant driving force behind AGB change in the desert steppe of the study area.