Satellite-Derived Estimation of Grassland Aboveground Biomass in the Three-River Headwaters Region of China during 1982–2018

Abstract

The long-term estimation of grassland aboveground biomass (AGB) is important for grassland resource management in the Three-River Headwaters Region (TRHR) of China. Due to the lack of reliable grassland AGB datasets since the 1980s, the long-term spatiotemporal variation in grassland AGB in the TRHR remains unclear. In this study, we estimated AGB in the grassland of 209,897 km² using advanced very high resolution radiometer (AVHRR), MODerate-resolution Imaging Spectroradiometer (MODIS), meteorological, ancillary data during 1982–2018, and 75 AGB ground observations in the growth period of 2009 in the TRHR. To enhance the spatial representativeness of ground observations, we firstly upscaled the grassland AGB using a gradient boosting regression tree (GBRT) model from ground observations to a 1 km spatial resolution via MODIS normalized difference vegetation index (NDVI), meteorological and ancillary data, and the model produced validation results with a coefficient of determination (R²) equal to 0.76, a relative mean square error (RMSE) equal to 88.8 g C m⁻², and a bias equal to −1.6 g C m⁻² between the ground-observed and MODIS-derived upscaled AGB. Then, we upscaled grassland AGB using the same model from a 1 km to 5 km spatial resolution via AVHRR NDVI and the same data as previously mentioned with the validation accuracy (R² = 0.74, RMSE = 57.8 g C m⁻², and bias = −0.1 g C m⁻²) between the MODIS-derived reference and AVHRR-derived upscaled AGB. The annual trend of grassland AGB in the TRHR increased by 0.37 g C m⁻² (p < 0.05) on average per year during 1982–2018, which was mainly caused by vegetation greening and increased precipitation. This study provided reliable long-term (1982–2018) grassland AGB datasets to monitor the spatiotemporal variation in grassland AGB in the TRHR.

1. Introduction

Grassland aboveground biomass (AGB) is defined as the organic matter resulting from grass photosynthetic activity [1,2,3]. It is an important indicator of carbon cycling over grassland ecosystems and provides technical and theoretical support for grassland utilization decisions and resource management [4,5,6,7]. In particular, grassland is a major component of the Three-River Headwaters Region (TRHR), which is an important nature reserve in the Qinghai-Tibet Plateau [8]. Spatiotemporal variation in AGB can reflect the carbon budget in grassland ecosystems [9,10]. Therefore, there is a considerable need to quantify the long-term spatiotemporal variation in grassland AGB in the TRHR to improve both the modeling and estimation of the carbon dynamics of grassland ecosystems [5].

Over the past 40 years, ground-based and satellite-based observations have been successfully used to support the efforts of estimating grassland AGB [11,12,13,14,15]. All ground-based methods and techniques only provide point-based observations that are limited at large spatial scales to show the heterogeneity of AGB [1,16]. Satellite-based observations can rapidly obtain data covering a wide region, so many studies have used various remotely sensed data to estimate grassland AGB at a regional scale [17,18,19,20]. One of the most common methods is to use satellite-based vegetation indices (VIs) [21,22,23,24,25]. Satellite-based VI data can directly indicate the condition of aboveground plant components [10]. Among the VI data, the normalized difference vegetation index (NDVI) is the most widely used index [26,27]. Therefore, the MODerate-resolution Imaging Spectroradiometer (MODIS) NDVI product has been widely used to estimate grassland AGB [28,29].

Table 1. Summary of relevant studies on estimating grassland AGB based on MODIS products in regions relative to TRHR or TP.

Study Domain	Time Horizon	Remotely Sensed Data	Models	Main Conclusions	References
TP	2005	MODIS NDVI product	Exponential model	Estimated precision of AGB in the TP is 80%.	[33]
TP	2001–2004	MODIS EVI product	EVI-based linear model	Approximately 54.2% of AGB variation in alpine grasslands could be explained by climatic variables and soil texture.	[34]
Golog Prefecture	2005.7–9	MODIS NDVI and EVI products	Exponential model	NDVI-based exponential model is optimal, and the precision was approximately 72%.	[43]
TP	2000–2012	MODIS NDVI product	Exponential model	The accuracy of model (AGB = 10.33e3.28NDVI) is R2 = 0.64.	[35]
TP	2011–2012	MODIS NDVI, EVI and reflectance products	General linear model	Large spatial variation of AGB in the TP was observed and AGB ranged from 1 to 136 g/m2.	[36]
TRHR	2003–2014	MODIS reflectance product	Multi-factors regression model	Multifactor model based on latitude, longitude, and grass cover is selected as the AGB model (R2 = 0.63).	[30]
TP	2000–2012	MODIS NDVI product	Linear and power models	For meadow grasslands, linear model is the best AGB estimation model (R2 = 0.51, RMSE = 139.5 kg/hm2). For steppe grasslands, power function is the best AGB estimation model (R2 = 0.63, RMSE = 146.5 kg/hm2).	[37]
NTP	2000–2013	MODIS FPAR/LAI product	Logarithm model	Aqua MODIS FPAR/LAI product may have the highest accuracy in AGB estimation.	[38]
TRHR	2001–2016	MODIS surface reflectance product	ANN model	ANN model based on grass cover, longitude, latitude, and NDVI has best accuracy (R2 = 0.62, RMSE = 589.9 kg/DW ha) in testing set.	[31]
Maqu	2016–2017	MODIS NDVI and EVI products	Linear model	R2 of the NDVI models were all greater than those of the EVI models.	[44]
NTP	2000–2014	MODIS NDVI product	logarithm and exponent models	The slopes of AGB in fenced grasslands and in livestock were −2.71 and −0.15 g m−2 year−1, respectively.	[39]
Shenzha	2000–2013	MODIS NDVI product	Exponential model	Annual average AGB was slightly increased and fluctuated between 1800 and 2700 kg/ha, with an average of 2272 kg/ha.	[45]
TP	2000–2014	MODIS NDVI and EVI products	RF model	Grassland AGB decreased from southeast to northwest with an average value of 77.12 g m−2 and increasing rate of 0.19 g m−2 year−1 in the TP.	[40]
TP	2000–2017	MODIS NDVI product	RF model	RF model with multiple factors performed best for AGB estimation (R2 = 0.75, RMSE = 346.06 kg/ha).	[41]
NTP	2004.9	MODIS NDVI and EVI products	NDVI-based exponential model	Optimal regression model for AGB estimation was NDVI-based exponential models (R2 = 0.63).	[42]
TRHR	2001–2019	MODIS NDVI product	HASM	HASM achieved the best results (R2 = 0.85, RMSE = 29 g m−2) and AGB increased by 1 g m−2 year−1.	[32]

Note: NTP: northern Tibetan Plateau; EVI: enhanced vegetation index; FPAR: fraction of photosynthetically active radiation; LAI: leaf area index; ANN: artificial neural network; HASM: high accuracy surface model; R²: coefficient of determination; RMSE: root mean square error; DW: dry weight; hm²/ha: hectare.

presents a summary of relevant studies focused on estimating grassland AGB based on MODIS products in the TRHR [30,31,32], Tibetan Plateau (TP) [33,34,35,36,37,38,39,40,41,42] or related regions [43,44,45]. These studies used many linear and nonlinear regression methods, including machine learning methods, to estimate grassland AGB with MODIS products. However, MODIS products cannot be used to estimate AGB before 2000 because of the lack of data. Solely relying on MODIS products has limited the monitoring of spatiotemporal variations in AGB in the TRHR prior to 2000.

Fortunately, advanced very high resolution radiometer (AVHRR) data, which have a relatively coarse spatial resolution (approximately 8 km × 8 km), have provided global satellite observations since 1981 [46]. Xiao et al. produced a new series of AVHRR NDVI product suites with a 5 km spatial resolution using the method of temporally continuous vegetation index-based land-surface reflectance reconstruction (VIRR) [47,48]. Therefore, AGB needs to be upscaled from ground observations to a 5 km spatial resolution. The upscaling methods consist of two fundamental steps [49,50,51,52]: (i) using statistical methods to establish a relationship between AGB ground observations and explanatory variables, including remotely sensed, meteorological, and ancillary data, and (ii) using gridded remotely sensed, meteorological, and ancillary data to obtain the spatiotemporal pattern of AGB by considering the established relationship. However, there is an obvious problem associated with the spatial representativeness of AGB upscaled from ground observations to a coarse spatial resolution directly [53]. A previous study [10] reported that the random forest (RF) model was used to upscale AGB directly from ground observations to an 8 km spatial resolution in the TP during 1982–2010. Due to the weak spatial representation of ground observations, this upscaling method leads to large uncertainties. Therefore, AGB is required to be upscaled from ground observations to a fine and medium spatial resolution and then from a fine and medium to coarse spatial resolution using multisource remotely sensed data with different spatial resolutions.

In this study, we upscaled grassland AGB from ground-based points to a 5 km spatial resolution using the MODIS NDVI product with a 1 km spatial resolution as a bridge to generate a long-term (1982–2018) grassland AGB product in the TRHR. We have three major objectives: (1) to upscale AGB from ground observations to a 1 km spatial resolution; (2) to upscale AGB from a 1 km to 5 km spatial resolution; and (3) to estimate long-term (1982–2018) grassland AGB of a 5 km spatial resolution in the TRHR.

2. Study Area and Datasets

2.1. Study Area and In Situ Data

The TRHR (89°31′~102°14′E, 31°38′~36°20′N) is situated in the backland of the Qinghai-Tibet Plateau, and it is the headwater source of the three famous rivers in China: the Yellow River, the Yangtze River, and the Lancang River (Figure 1). The total area is 302,500 km². The total grassland area is 209,897 km², while the area of usable grassland is 178,744 km². The area has an average elevation of over 4000 m and an annual temperature of −4.4 °C. Its annual precipitation is 262.2–772.8 mm and occurs mostly from June to September each year [54]. The grass growth period is less than 3 months when there is no absolute frost-free period. The major types of grassland include alpine meadow steppe and alpine steppe, which account for approximately 76% and 23% of the total grassland area, respectively [55].

Ground observations of grassland AGB were used to train the models and evaluate the accuracy of the AGB estimates. AGB ground observations were derived from the Qinghai Eco-Environment Monitoring Center, including a total of 125 effective AGB observations for the period of June to September, the period of grass growth in 2009. The 125 AGB observation sites were divided into two groups: 75 training sites and 50 validation sites. Each grassland sampling site consisted of a 1000 m² area with an approximate dominant grass type. Each site included 9 grass sample plots that were 1.0 m × 1.0 m in size. All the observation sites were all-natural grasslands with no human interaction during the study period. The data recorded from each sampling site included latitude, longitude, elevation, grass type, coverage, grass height, fresh weight, and dry weight (oven-dried at 65 °C for 48 h in the lab) [30]. The average of AGB observations from nine grass sample plots represented the AGB of each grassland sampling site.

2.2. Remotely Sensed Data

We used a series of AVHRR NDVI product suites with an 8-day temporal resolution and a 5 km spatial resolution produced by Xiao et al. [47]. The NDVI product representing the period of 1982–2018 was generated from AVHRR data by a refined VIRR method [48].

We also used the MODIS NDVI product (MOD13A2) with a 16-day temporal resolution and a 1 km spatial resolution [56] for 2009–2013 as a gridded level-3 product in the sinusoidal projection obtained from the Land Processes Distributed Active Archive Centre (https://lpdaac.usgs.gov, accessed on 15 September 2020). We extracted the NDVI values corresponding to the ground-based sampling sites, which can be used as an independent variable for training the models.

2.3. Meteorological Data

In this study, we used the China Meteorological Forcing Dataset (CMFD) containing air temperature (Ta), precipitation (P), and downward shortwave radiation (Rs) developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences [57]. These products cover the period of 1982–2018, with a spatial and temporal resolution of 10 km and 3 h. We resampled the dataset of the study area to a 1 km (5 km) spatial resolution using the nearest neighbor method. We also extracted the CMFD Ta, P, and Rs values corresponding to the ground-based sampling sites, which can be used as independent variables for training the models.

The CMFD Ta dataset was constructed by merging observations from 740 meteorological stations and the corresponding Princeton meteorological forcing data [58]. The CMFD P dataset was developed by combining three precipitation datasets including those of the same 740 meteorological stations, the Tropical Rainfall Measuring Mission (TRMM) 3B42 precipitation products [59], and the Asian Precipitation-Highly Resolved Observational Data Integration Towards Evaluation of the Water Resources project (APHRODITE) [60]. The CMFD Rs dataset was constructed by Global Energy and Global Energy and Water cycle Experiment-Surface Radiation Budget (GEWEX-SRB) radiation data and meteorological station measurements [61]. The 3-hourly Ta, P, and Rs were subsequently aggregated into daily means. Monthly P values were obtained by multiplying the average daily mean by 30 or 31 (days) [62].

2.4. Ancillary Data

We used digital elevation model (DEM) data to develop the upscaling model. The DEM data consist of Global Multiresolution Terrain Elevation Data 2010 with a 30 arcsecond spatial resolution in the GeoTIFF format, which can be obtained from ENVI 5.3. We resampled the elevation data from a 30 arcsecond to 1 km and 5 km spatial resolution using a simple average method.

We also used the UMD Land Cover Classification (UMD-LCC) dataset developed in 1998 by the Department of Geography at the University of Maryland [63]. The dataset was divided into 14 different land cover types. We used the UMD-LCC dataset with a 1 km spatial resolution and resampled the dataset of the study area to a 5 km spatial resolution using the nearest neighbor method.

3. Methodology

3.1. AGB Upscaled Procedure

Figure 2 shows the procedure and evaluation of upscaling AGB based on the gradient boosting regression trees (GBRT) model. The first step was to prepare the MODIS NDVI product (1 km), AVHRR NDVI product (5 km), CMDF Ta, P, Rs (1 km and 5 km), elevation dataset (1 km and 5 km), and ground AGB observation dataset. The second step was upscaling grassland AGB from ground observations to a 1 km spatial resolution based on GBRT model using all the input and output variables. Note that in the training process, the input variables were MODIS NDVI, P, Ta, Rs, and elevation at site scale, and the output was the ground AGB observation dataset collected at 75 sampling sites. Then, the developed model was applied to generate AGB with a 1 km spatial resolution from 2009 to 2013. Then, the upscaled AGB from MODIS NDVI, P, Ta, Rs, and elevation with a 1 km spatial resolution was evaluated using the ground AGB observations at the other 50 sampling sites as reference values. The third step was upscaling grassland AGB from a 1 km to 5 km spatial resolution using all the input and output variables based on the GBRT model. Note that in the training process, the input variables were AVHRR NDVI, P, Ta, Rs, and elevation with a 5 km spatial resolution, and the output was the MODIS-derived AGB with a 5 km spatial resolution from 2009 to 2011, which was aggregated from the MODIS-derived AGB with a 1 km spatial resolution using a simple average method. Then, the developed model was applied to generate AGB with a 5 km spatial resolution during 1982–2018. The upscaled AGB from AVHRR NDVI, P, Ta, Rs, and elevation with 5 km spatial resolution was evaluated using the MODIS-derived AGB with a 5 km spatial resolution during 2012–2013 as reference values. The final step was to detect the spatiotemporal variations in grassland AGB in the TRHR.

3.2. GBRT

In this study, we used the GBRT algorithm [64] to generate a long-term AGB dataset. The GBRT model, also referred to as the gradient boosting decision tree, is an ensemble learning approach that combines multiple decision trees to generate an estimation model in gradient boosting. It produces competitive, highly robust, interpretable procedures for regression, especially appropriate for mining less than clean data [64]. The GBRT model has been applied to AGB estimation in global forests [65] but not in grasslands in the TRHR. Therefore, we selected GBRT to generate a long-term AGB dataset in the TRHR. The approximation function according to GBRT can be expressed as [64]:

$$F(x;{\{{\beta }_m,a_m\}}_1^M)={\displaystyle \sum _{m=1}^M{\beta }_{m}h(x;a_m)}$$

(1)

where ${\beta }_m$ refers to the weight of a single decision tree $h(x;a_m)$. Each base learner is a J-terminal node regression tree. Each regression tree model itself has the additive form as follows [64]:

$$h(x;{\{b_j,R_j\}}_1^J)={\displaystyle \sum _{j=1}^J{b}_{j}I(x\in R_j)}$$

(2)

where ${\left\{R_j\right\}}_1^J$ are disjoint regions that collectively cover the space of all joint values of predictor variables x. These regions are represented by the terminal nodes of the corresponding tree. Indicator function I (·) has a value of 1 if its argument is true and zero otherwise. The “parameters” of this base learner (2) are the coefficients ${\left\{b_j\right\}}_1^J$ and the quantities that define the boundaries of the regions ${\left\{R_j\right\}}_1^J$. These are the splitting variables and the values of those variables that represent the splits at the nonterminal nodes of the tree. Because the regions are disjointed, (2) is equivalent to the prediction rule: if $x\in R_j$, then $h(x)=b_j$.

3.3. Other Machine Learning Methods

We replaced the GBRT in the AGB upscaling model by other machine learning algorithms, including random forest (RF) or extremely randomized tree (ERT) to compare the performance of models.

(1) RF: RF is a combination of tree predictors that every tree depends on the values of a random vector sampled independently and with the same distribution for all trees in the forest. [66,67]. The single tree predictor in the entirety is produced by stochastically choosing independent variables with the bootstrap sampling approach [68]. The RF method usually ensures optimal predictions through three main steps. First, original data are divided into subdivisions. Second, every subdivision is used to develop an individual decision tree representing a subregression model that provides its own results. Finally, outputs derived from the trees are averaged to estimate the optimal results [69]. The randomness added to the forests when coming into being decision trees may decouple estimated errors. Some errors can be removed by having an average of those estimations [70]. RF for regression are formed by growing trees depending on a random vector Θ such that the tree predictor h(x, Θ) takes on numerical values as opposed to class labels. The output values are numerical, and we assume that the training set is independently drawn from the distribution of the random vector Y, X. The mean-squared generalization error for any numerical predictor h(X) is:

$$E_{X,Y}{(Y - h(X\left)\right)}^2$$

(3)

The random forest predictor is formed by taking the average over k of the trees $\{h(x,{\Theta }_k)\}$ [66].

(2) ERT: ERT, which was raised by Geurts, is a tree-based machine learning ensemble model for supervised learning of classification and regression [71]. It consists of strongly randomizing both attribute and cut-point choices while splitting a tree node essentially. In the extreme case, it builds totally randomized trees with structures that are independent of the output values of the learning sample. The strength of randomization can be tuned to problem specifics by the appropriate choice of a parameter. ERT has three important parameters: K, $n_{\min }$, and M [71]. The parameter K is the number of random splits at every node to build extra trees. It may be selected from one to n, where n is the number of attributions. According to a given situation, the smaller the K is, the stronger the randomization of the trees and the weaker the dependence of their structure on the output values of the learning sample will be. Previous studies have demonstrated that the best value of K is K = n for regression [71]. The parameter $n_{\min }$ of original samples asked for splitting a node. Smaller trees, higher bias, and smaller variance result from the larger values of $n_{\min }$. Therefore, optimal value of $n_{\min }$ depends on the principle under the level of output noise in the dataset [71]. The parameter M refers to the amount of ensemble trees. It is well known that the behavior of estimation error is a monotonically decreasing function of M for randomization methods [66]. Relative variance reduction is used in regression. S_l and S_r refer to two subsets of cases from S corresponding to the two outcomes of a split s, then the score is defined as follows [71]:

$$Scor{e}_R(s,S)=\frac{\mathrm{var}\left\{y\right|S\}-\frac{\left|S_l\right|}{\left|S\right|}\mathrm{var}\{y\left|S_l\right\}-\frac{\left|S_r\right|}{\left|S\right|}\mathrm{var}\left\{y\right|S_r\}}{\mathrm{var}\left\{y\right|S\}}$$

(4)

where $\mathrm{var}\left\{y\right|S\}$ is the variance of the output y in the sample S.

3.4. Statistical Metrics

To evaluate the accuracy of the three machine learning upscaling methods, the coefficient of determination (R²), root mean square error (RMSE), and bias were selected as the evaluation metrics [72]. R², RMSE and bias can be written as:

$$R^2=\frac{{\displaystyle {(\sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y}))}^2}}{{\displaystyle \sum _{i=1}^N{(x_i-\overline{x})}^2\sum _{i=1}^N{(y_i-\overline{y})}^2}}$$

(5)

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{(x_i-y_i)}^2}}$$

(6)

$$bias=\frac{{\displaystyle \sum _{i=1}^N(x_i-y_i)}}{N}$$

(7)

where $N$ refers to the number of samples and $x_i$ and $y_i$ refer to the ith reference and estimated AGB, respectively. $\overline{x}$ and $\overline{y}$ refer to the average value of the reference and estimated AGB, respectively.

The sequential Mann-Kendall test was introduced by Sneyres [73] to discover abrupt potential turning points in long-term data. We used it to detect the turning point of long-term AGB variation. This test includes the calculation of two data series: one is progressive, and the other is retrogressive. When the two data series cross each other and deviate beyond the confidence limitation of 95%, then, that point will be a statistically significant turning point. It is calculated by comparing the values of x_j annual mean time series (j = 1,…, n) with x_i, (i = 1, …, j − 1) [74]. At the time of comparison, the cases where x_j > x_i is counted and denoted by n_j.

Then, test statistic t is calculated using the equation as follows:

$$t_j={\displaystyle \sum _1^j{n}_j}$$

(8)

The mean and variance of the test statistic are

$$E(t)=\frac{n(n-1)}{4}$$

(9)

And

$$Var(t_j)=[j(j-1)(2j+5)]/72$$

(10)

The sequential values of the statistic u(t) are

$$u(t)=\frac{t_j-E(t)}{\sqrt{Var(t_j)}}$$

(11)

The values of u′(t) retrogressive series are calculated similarly but starting from end of the series.

Linear trends were applied to quantify the interannual and pixel trends of the AGB and other variables. The equation can be written as [75]:

$$y_t=y_0+b{x}_t$$

(12)

where $y_t$ is the changed value; $y_0$ is the initial value; $x_t$ is the time interval; and b is the trend. Statistically significant differences were set as p < 0.05 unless otherwise stated.

4. Results

4.1. Evaluation of the AGB Upscaling Model from Ground Observations to a 1 KM Spatial Resolution

To upscale the AGB from ground observations to a 1 km spatial resolution using the GBRT, RF, and ERT models, the ground measurements collected at the 75 training sites were used to train the models. Figure 3a–c presents the scatter plots for the AGB observations and training results using the three models. GBRT yields the best training performance with the highest R² (0.95) and the lowest RMSE (39.3 g C m⁻²) of the AGB training results against ground observations among the three models, followed by RF with an R² of 0.87 and RMSE of 63.4 g C m⁻². The descending training performance order of the three models is GBRT, RF, and ERT.

To evaluate the three AGB upscaling models, the corresponding MODIS-derived AGB estimations with a 1 km spatial resolution were directly compared with AGB ground observations. Figure 3d–f shows the scatter plots for the AGB observations and estimations using the GBRT model and the other two models at the other 50 validation sites. The results show that the GBRT model yields the best validation performance with the highest R² (0.76), lowest RMSE (88.8 g C m⁻²) and lowest bias (−1.6 g C m⁻²) of the AGB estimations against ground observations among the three upscaling models. The accuracy of the RF model is lower than that of GBRT but higher than that of ERT, with an R² of 0.62, an RMSE of 112.0 g C m⁻², and a bias of −5.0 g C m⁻². The ERT model yields the lowest R² (0.62), the highest RMSE (114.3 g C m⁻²), and the highest bias (−6.0 g C m⁻²). The results indicate that the GBRT model performs best for estimating grassland AGB at a 1 km spatial resolution. Therefore, we selected the GBRT model to upscale AGB from ground observations to a 1 km spatial resolution. Then, we used the GBRT model driven by the MODIS NDVI product, elevation, Ta, P, and Rs with a 1 km spatial resolution to generate AGB datasets during 2009–2013.

4.2. Evaluation of the AGB Upscaling Model at a 1 km to 5 km Spatial Resolution

To upscale the AGB from a 1 km to 5 km spatial resolution using the GBRT, RF, and ERT models, the MODIS-derived reference AGB during 2009–2011 was aggregated from 1 km to 5 km using the simple average method and then used to train the models. Figure 4a–c presents the scatter plots for the MODIS-derived reference AGB and training results using the three models. RF yields the best training performance with the highest R² (0.76) and the lowest RMSE (54.7 g C m⁻²) of the AGB training results against the MODIS-derived reference AGB among the three models, followed by GBRT with an R² of 0.76 and RMSE of 55.3 g C m⁻². The ERT model yields the lowest R² (0.72) and the highest RMSE (59.7 g C m⁻²).

To evaluate the performance of the AGB upscaling model based on the GBRT model and the other two models, we validated the three models based on the MODIS-derived reference AGB estimation with a 5 km spatial resolution during 2012–2013. Figure 4d–f shows the scatter plots for the MODIS-derived reference AGB and AVHRR-derived estimated AGB by the three models with a 5 km spatial resolution. The results show that the GBRT model demonstrated the best performance with the highest R² (0.74) and lowest RMSE (57.8 g C m⁻²), and bias (−0.1 g C m⁻²) of the AVHRR-derived estimated AGB against the MODIS-derived reference AGB among the three models in the validation set. Therefore, we selected the GBRT model to predict the AGB based on AVHRR NDVI product, elevation, Ta, P and Rs with a 5 km spatial resolution to generate AGB datasets during 1982–2018.

Figure 5 shows the MODIS-derived reference and AVHRR-derived estimated average AGB maps with a 5 km spatial resolution determined by the GBRT model in the validation set during 2012–2013. The results confirmed that the spatial pattern of the estimated AGB was consistent with that of the reference AGB, so it could reflect the spatial distribution of actual grassland AGB, which could confirm the reliability of the estimated AGB datasets based on AVHRR data. Therefore, the GBRT model is selected to produce the AGB dataset during 1982–2018 using the AVHRR NDVI product, Ta, P, Rs, and elevation with a 5 km spatial resolution.

4.3. Importance of Input Variables to AGB Upscaling Model

Evaluation of the importance of different variables (i.e., NDVI, P, Rs, Ta, and elevation) to the upscaled AGB with a 5 km spatial resolution is important to understand the effect of variables on the performance of the GBRT models for AGB upscaling. The importance of variables was described by the increased percentage of RMSE relative to the baseline RMSE for which all the variables, including NDVI, P, Rs, Ta, and elevation, were used to develop the GBRT-based AGB upscaling model, indicating to what extent the RMSE increases when a variable is not involved in the model. The larger the value of the increased RMSE is, the more important the variable is to the AGB upscaling model. The GBRT models with five combinations of input variables and MODIS-derived reference AGB with a 5 km spatial resolution during 2009–2011 were trained independently to achieve the best model performance for each combination itself. RMSE values of the GBRT-based upscaled AGB relative to the MODIS-derived reference AGB with a 5 km spatial resolution during 2012–2013 for the five combinations were subsequently calculated to determine the relative importance of each variable. Figure 6 demonstrates that NDVI and P are the two most important variables (32.9% and 26.5%, respectively) in the AGB upscaling model. These results fully reflect the relationship between AGB and NDVI, as well as AGB and P in the TRHR, which is strongly consistent with most study results [30,31,32]. Zhu et al. [76] demonstrated that in 92.5% of the TRHR area, the vegetation was significantly affected by the severity of the drought during 2000–2018. Resulting from the soil water availability, the strong relationship between AGB and P is important in promoting plant production in dry grassland biomes [77]. In addition, Rs, Ta, and elevation show a relatively low (5.3%, 6.9%, and 4.6%, respectively) importance to the AGB upscaling model. Although the contribution of these variables is small, the best estimations can only be acquired using all five variables, because they are important for plant photosynthesis [78].

4.4. Spatiotemporal Variation in Grassland AGB during 1982–2018 in the TRHR

Figure 7 shows the overall spatial patterns of the annual average grassland AGB in the TRHR during 1982–2018. The annual average AGB decreases from the southeast to the northwest. A total of 3.2% of the total area in the grassland of the TRHR with a high AGB distribution (120–180 g C m⁻²) is mainly located in the far east region of the study area (Zeku, Henan and Jiuzhi). Figure 7 clearly shows a low AGB distribution (0–60 g C m⁻²), which occupies approximately 72.6% of the northern, western, and central areas of the TRHR (Qumalai, Maduo, Zhiduo, Tanggulashan, Zaduo and Chengduo). The southern and eastern areas of the TRHR (Yushu, Nangqian, Xinghai, Tongde, Maqin, Gande, Dari and Banma) have a relatively moderate AGB of approximately 60–120 g C m⁻² covering approximately 24.2% of the study area. The AGB distribution is mostly concentrated at 0–60 g C m⁻², and the number of AGB values larger than 100 g C m⁻² is relatively low.

Figure 8 presents the spatial patterns of the multiyear average seasonal AGB distribution. The AGB reached its highest value (320–360 g C m⁻²) in summer (JJA, June, July, and August). The change trend of the AGB spatial distribution in summer (Figure 8b) is similar to the multiyear patterns of AGB during 1982–2018 (Figure 7). The AGB is less than 80 g C m⁻² in spring (MAM, March, April, and May) and winter (DJF, December, January, and February). The seasonal AGB has overall moderate values of less than 200 g C m⁻² in autumn (SON, September, October, and November).

Figure 9 presents the turning point of AGB temporal variation was found in 1998 by the sequential Mann–Kendall test and the interannual variation in AGB, P, and NDVI over the period of 1982–2018. The annual AGB, P and NDVI in the TRHR increased over 1982–2018 with a linear fit with slopes of 0.37 g C m⁻² year⁻¹, 3.66 mm year⁻¹, and 0.9 × 10⁻³ year⁻¹, respectively (p < 0.05). An increase in annual AGB and NDVI occurred during 1982–1998, with linear trends of 0.16 g C m⁻² and 0.30 × 10⁻³ year⁻¹ (p < 0.05). Such an increase in the annual AGB is consistent with the NDVI increase in response to global warming. Coincident with the El Niño event in 1998, the positive trend in AGB seems to have accelerated (0.21 g C m⁻² year⁻¹, p < 0.05) during 1998–2018. The time slope of P was found to have changed from being negative (−2.81 mm year⁻¹, p < 0.05) during 1982–1998 to positive (5.02 mm year⁻¹, p < 0.05) during 1998–2018. However, since 1998, the increase in the annual NDVI seems to have stalled, with a 95% confidence interval trend of −0.50 × 10⁻³ year⁻¹ (p < 0.05).

Figure 10a shows the spatial distribution of AGB trends in the TRHR during 1982–2018. Our results illustrate that the AGB of the TRHR increased by 0–1.5 g C m⁻² year⁻¹ during 1982–2018. A total of 82.3% of the pixels in the grassland of the study area show an increasing trend in AGB. Except for the far west of the study area, the other regions have a significantly increasing trend (~79.4%, p < 0.05). The largest AGB increases occurred in the eastern TRHR (Figure 10a), where the annual AGB is also the highest, as shown in Figure 7. The overall increment in the annual AGB can largely be attributed to the upward trend of NDVI and P during 1982–2018 Figure 11a and Figure 12a).

We also analyzed changes in the annual AGB trends before and after 1998. Figure 10b shows that the AGB has the main distribution with a significantly increasing trend in the eastern, southern, and central parts of the study area (~43.9%, p < 0.05) between 1982 and 1998, which is relative to the regional vegetation (mainly grass) greening with the NDVI increase (Figure 11b). Figure 10c illustrates that the main regions of the TRHR show a significantly increasing AGB trend (~60.5%, p < 0.05), which is mainly attributed to the upward trend of P occurring in most areas of the TRHR during 1998–2018 (Figure 12c).

5. Discussion

5.1. Performance of AGB Upscaling Model

To meet various requirements of AGB estimation at grid scales in practice, coarser-spatial-resolution AGB estimates such as those at 1 km or 5 km resolution may be needed. In this study, to select the best AGB upscaling model, we used MODIS and AVHRR NDVI products and Ta, P, Rs, elevation, and AGB ground observations to establish the RF, ERT, and GBRT models for estimating the grassland AGB in the TRHR. The validation results illustrated that the GBRT model has the best performance among the three models, judging by the R² (0.74) and RMSE (57.8 g C m⁻²) of the AVHRR-derived estimated and MODIS-derived reference AGB with a 5 km spatial resolution (Figure 4d–f). Model validation demonstrated that the GBRT algorithm could improve the accuracy to prevent overfitting on the new data [64], which is corresponding to the reported results. Yang et al. [65] used GBRT as the optimal algorithm to estimate global forest AGB after comparing the performance of ANN and multivariate adaptive regression spline (MARS) algorithms. The validation results showed that the R² and RMSE between the estimated and reference AGB with 1 km spatial resolution were 0.90 and 35.87 Mg/ha, respectively. Bai et al. [78] used the GBRT algorithm to estimate the global monthly gross primary production (GPP) at a 5 km spatial resolution with the OCO-2 SIF product (GOSIF) and other explanatory variable data. Compared with the RF, MARS, light gradient boosting machine (LightGBM), and long short-term memory (LSTM) algorithms, GBRT was more suitable to estimate global GPP by GOSIF. The GBRT-based GPP estimations derived from GOSIF were validated by in situ measurements with highly accurate R² (0.58), RMSE (2.74 g C·m⁻²), and bias (–0.34 g C·m⁻²) values. Our conclusions are similar to previous findings. This may be because the GBRT algorithm is an iterative decision tree that utilizes the additive model and the forward distribution algorithm to combine several weak classifiers into a strong classifier in different proportions [64,79].

We compared estimated AGB errors to judge the performance of GBRT, RF, and ERT. Figure 13 illustrates the probability density distributions of the AVHRR-derived estimated AGB errors against the MODIS-derived reference AGB with a 5 km spatial resolution by the GBRT and the other two AGB upscaling models in the validation set. Figure 13 shows that all three predictive errors are closely centered on zero, even in the ERT model, which obtained the largest bias (0.6 g C m⁻²) (Figure 4f). The GBRT model shows a competitive and superior predictive performance compared to RF and ERT. This may be because the GBRT algorithm has the advantage of producing competitive, highly robust, interpretable procedures for regression [64]. Therefore, the GBRT algorithm can capture the AGB variance accurately and show the optimal estimation performance.

Training samples with spatial representativeness were essential for modeling the relationships between AGB and remotely sensed data with a coarse spatial resolution. To generate reliable long-term grassland AGB datasets, we proposed upscaling AGB data from ground observations to a dataset with a 1 km spatial resolution using MODIS data and then from a 1 km to 5 km spatial resolution using AVHRR data. A previous study confirmed this upscaling strategy [80]. Fractional vegetation cover (FVC) is a basic indicator of biomass production in ecosystems [81]. To enhance the spatial representativeness of FVC training samples for spatial matching with AVHRR reflectance data, Jia et al. [80] reprojected the Global LAnd Surface Satellite (GLASS)-MODIS FVC data from a spatial resolution of 500 m to 5 km as the training samples to generate a long-term GLASS-AVHRR FVC product since 1981 by MARS. Based on the use of field FVC measurements as the reference data, the validation results indicated that the performance of the GLASS-AVHRR FVC product (R² = 0.834, RMSE = 0.145) was slightly superior to that of the long-term GEOV1 FVC product (R² = 0.799, RMSE = 0.174). In this study, we used the established relationship between ground AGB observations and independent variables to generate a MODIS-derived AGB dataset with a 1 km spatial resolution. To enhance the spatial representativeness of the training set for generating long-term AGB datasets with a 5 km spatial resolution, we aggregated the MODIS-derived reference AGB from a 1 km to 5 km spatial resolution with a simple average method as the training set of the GBRT model. Compared with a previous study for upscaling AGB from ground observations to an 8 km spatial resolution directly [10], our upscaling strategy has the ability to generate a training set with stronger spatial representativeness to reduce the uncertainties [80].

5.2. Spatiotemporal Variation in Grassland AGB during 1982–2018

In this study, the spatial distribution of grassland AGB confirmed previous findings that the grassland AGB decreased gradually from southeast to northwest in the TRHR [31]. We found a strong spatial correlation between the distribution of Ta and P in the TRHR. This was also documented by Bei et al. [82]. Figure 9b presents that the grassland AGB increased 0.37 g C m⁻² per year during 1982–2018 in the TRHR. However, the main cause of the interannual increasing trend of grassland AGB is still controversial [37,40,41]. Zeng et al. [40] suggested that the interannual increasing trend of grassland AGB in the TP was significantly and positively correlated with temperature during 2000–2014. Gao et al. [41] suggested that the national conservation policies of returning cultivated land to forests or grasslands since 1999 and grazing ban laws since 2004 may have promoted grassland restoration [83], which caused the increasing trend of grassland AGB in the TP during 2000–2017. Liu et al. [37] presented that the primary factors influencing the increasing trend of grassland AGB in the TP were anthropogenic activities and climate warming during 2000–2012. Reflected by our results, NDVI and P are the two primary factors influencing the interannual variation in grassland AGB.

Our study also discussed the turning point in the long-term temporal dynamics of grassland AGB in the TRHR during 1982–2018 (Figure 9). The turning point of AGB was identified as occurring in 1998 by the sequential Mann–Kendall test (Figure 9a). Before 1998, the increasing trend of AGB was mainly attributed to the increasing NDVI. The positive relationship between NDVI and AGB in the TRHR was also documented by Liang et al. [30]. The AGB increased rapidly in 1998, coincident with the El Niño event, when the rising P was matched by an increasing available soil water content [82,84]. After 1998, the increasing AGB coincided with an increasing P. The annual average AGB was relatively higher in the years with higher mean annual P (e.g., 2012 and 2014) than in the surrounding years (Figure 9b). Previous studies demonstrated that AGB increased with P in the TRHR or TP [30,34], which provides evidence for our results. Therefore, P was the main contributor to the increasing trend of AGB during 1998–2018.

5.3. Limitations and Outlooks for Future Study

Our strategy was based on the GBRT model, a machine learning method relying on training sample datasets. The performance of the data-driven model usually depends on the quantity of the training set [85]. The number and distribution uniformity of sampling sites also affected the accuracy of estimation [31]. Due to limited transportation conditions, many depopulated zones distributed in the western region of the TRHR are inaccessible. There were very few field surveys taken in these regions, which might result in lower accuracy compared to other regions [41]. Our ground observations were captured only during the growth period of grassland in 2009. The lack of AGB ground observations increases the uncertainty of the estimated AGB results in the non-growing season. In future studies, the AGB observation data of non-growing season can be properly supplemented to reduce the uncertainty of the estimation results. It will be helpful to monitor the long-term spatiotemporal variation of AGB in the TRHR.

Although the GBRT model achieves a higher accuracy, the physical mechanism of the GBRT to estimate AGB is unclear. Undoubtedly, ecosystem and physical models can better explain the spatial and temporal heterogeneity of ecosystem biomass variation, although they still have large uncertainties [86,87]. In future work, we will focus on combining machine learning models, such as GBRT, with physical-based methods to establish a better grassland AGB inversion model for upscaling the ground observations to the target grid scale, which could be explained by the internal physical mechanisms, and obtain better accuracy.

There was a time difference between satellite and ground observations, even if we made great efforts to select the MODIS data close to the field measurement time [30]. The spatial scale mismatch between the AGB observations and satellite-based remotely sensed data would also introduce uncertainties to upscaled AGB results.

6. Conclusions

Based on the AGB upscaling model from GBRT, we generated long-term (1982–2018) grassland AGB datasets with a 5 km spatial resolution by remotely sensed data with different spatial resolutions (MODIS and AVHRR NDVI products), meteorological data (Ta, P, and Rs), ancillary data (elevation), and ground AGB observations in the TRHR. Compared with previous studies, our study can effectively reflect the long-term spatiotemporal variation in AGB during 1982–2018 in the TRHR. Major conclusions are given as follows.

(1)

MODIS-derived upscaled AGB with a 1 km spatial resolution based on the GBRT model was validated by AGB ground observations with the best performance (R² = 0.76, RMSE = 88.8 g C m⁻², and bias = −1.6 g C m⁻²).

(2)

The GBRT model also showed the best validation performance between the AVHRR-derived upscaled and MODIS-derived reference AGB with a 5 km spatial resolution with an R² of 0.74, an RMSE of 57.8 g C m⁻², and a bias of −0.1 g C m⁻².

(3)

The annual trends of grassland AGB increased by 0.37 g C m⁻² year⁻¹ on average during 1982–2018, which was mainly attributed to vegetation greening and increasing precipitation.

In summary, our study successfully generated reliable long-term grassland AGB datasets in the TRHR covering the period of 1982–2018, thereby meeting the need for monitoring the long-term spatiotemporal variation in grassland AGB for grassland utilization decisions and resource management in the TRHR.