Evaluation of Remote Sensing Inversion Error for the Above-Ground Biomass of Alpine Meadow Grassland Based on Multi-Source Satellite Data

Abstract

It is not yet clear whether there is any difference in using remote sensing data of different spatial resolutions and filtering methods to improve the above-ground biomass (AGB) estimation accuracy of alpine meadow grassland. In this study, field measurements of AGB and spectral data at Sangke Town, Gansu Province, China, in three years (2013–2015) are combined to construct AGB estimation models of alpine meadow grassland based on these different remotely-sensed NDVI data: MODIS, HJ-1B CCD of China and Landsat 8 OLI (denoted as NDVI_MOD, NDVI_CCD and NDVI_OLI, respectively). This study aims to investigate the estimation errors of AGB from the three satellite sensors, to examine the influence of different filtering methods on MODIS NDVI for the estimation accuracy of AGB and to evaluate the feasibility of large-scale models applied to a small area. The results showed that: (1) filtering the MODIS NDVI using the Savitzky–Golay (SG), logistic and Gaussian approaches can reduce the AGB estimation error; in particular, the SG method performs the best, with the smallest errors at both the sample plot scale (250 m × 250 m) and the entire study area (33.9% and 34.9%, respectively); (2) the optimum estimation model of grassland AGB in the study area is the exponential model based on NDVI_OLI, with estimation errors of 29.1% and 30.7% at the sample plot and the study area scales, respectively; and (3) the estimation errors of grassland AGB models previously constructed at different spatial scales (the Tibetan Plateau, Gannan Prefecture and Xiahe County) are higher than those directly constructed based on the small area of this study by 11.9%–36.4% and 5.3%–29.6% at the sample plot and study area scales, respectively. This study presents an improved monitoring algorithm of alpine natural grassland AGB estimation and provides a clear direction for future improvement of the grassland AGB estimation and grassland productivity from remote sensing technology.

1. Introduction

As the largest terrestrial biome on the Earth’s surface [1], the grassland biome occupies approximately 40% of the total land area [2], and its net primary productivity accounts for approximately 20% of the production capability of the entire land biome [3]. In China, the grassland biome accounts for approximately 1/3 of the total national territory. Grasslands play a critical role not only in animal husbandry, but also in energy exchange, carbon sequestration and the biogeochemical cycle. In particular, the effect of the grassland biome is more prominent in the vast central and western regions of China [4,5].

Grassland vegetation is defined as permanent vegetation composed of an herbaceous plant community [2]; the above-ground biomass (AGB) of grassland provides the basis for estimating the net primary productivity of grassland [6]. Overall, AGB can be used to directly estimate grassland productivity [7,8]. In general, AGB is defined as the weight of dry grass in the above-ground portion within one unit area [9]. Monitoring using remote sensing data is the most effective method for collecting continuous spatial and temporal data on the regional or global scale [10,11], because satellite remote sensing can provide large-scale, frequent, low cost and massive information [12]. Therefore, it has gradually replaced traditional methods of ground biomass monitoring, which are inefficient and expensive. Since NDVI was first applied to study natural grasslands in the 1970s, research on the linkage between vegetation indices and AGB has had a history extending over several decades [13,14,15,16,17]. Because grassland AGB estimations can directly guide livestock production, evaluating the accuracy of estimation models is highly important. Numerous studies have indicated that factors such as the representativeness of ground sampling sites, the temporal and spatial resolutions of satellite images, the types of sensors and the methods of remote sensing image processing [18] are responsible for the large differences between the grassland AGB models developed by different scholars even for the same region or the same type of grassland. These differences exist not only in the form of the models, but also in the estimation errors. At present, grassland AGB estimation is often based on the statistical models on a specific remote sensing vegetation index [19,20,21,22]. For alpine meadow grassland, which type of remote sensing estimation model is the most suitable is currently unknown. Existing studies indicate that using both high spatial resolution imagery and filtering the vegetation index can increase the accuracy of grassland AGB estimation [18,23,24,25]. However, it is not yet clear whether there is any difference in using remote sensing data of different spatial resolutions and filtering methods to improve the AGB estimation accuracy of alpine meadow grassland. Therefore, further comparisons and studies must be conducted.

In this study, considering the factors discussed above, a region of alpine meadow grassland in the Tibetan Plateau area is used to perform the following investigations: (1) by comparing and analyzing the AGB estimation models and their accuracies based on the NDVI of multi-source remote sensing data, the influences of different remote sensing data and filtering methods on the error of grassland AGB estimation are revealed; (2) the data from the sample plots observed inside the study area are used to validate the applicability of previous grassland models based on MODIS data and to investigate the reasons for errors from different models; and (3) based on the above research results, we propose a method to improve the accuracy of grassland AGB estimation.

2. Data and Methods

2.1. Study $A=\pi r^2$

The study area (102°23′–102°26′ E, 35°5′–35°7′ N) is located in the Yangji Community of Sangke Town in Xiahe County, Gansu Province, with a size of approximately 3.86 km (N–S) × 2.77 km (E–W) and a mean elevation of 3050 m (Figure 1). The natural grassland type in the study area is alpine meadow. The dominant pasture plants include alpine Kobresia pygmaea C B Clarke, Elymus nutans Griseb, Festuca ovina L, Poa annua L, Koeleria litwinowii Domin var and Astragalus membranaceus Bunge, while the primary malignant weeds include Ligularia virgaurea Mattf, Leontopodium japonicum, Potentilla chinensis Ser and Pedicularis resupinata L. The dominant grazing livestock are yak and Ganjia sheep. The study area is cold and wet throughout the year, and it belongs to the continental monsoon climate of the temperate plateau. The annual average temperature is 2.1 °C; the annual average rainfall is 580 mm.

The study area consists of five fenced experimental sample areas (A–E) of natural alpine meadow grassland (Figure 1), with a total area of 161.36 ha. Specifically, Area A (19.38 ha) was used to conduct natural grassland reseeding tests; Area B (16.06 ha) was used for grazing utilization experiments; Area C (7.52 ha) acted as a non-treatment testing control area; Area D (19.30 ha) was an experimental area for grassland enclosure; and Area E (99.10 ha, accounting for 61.55% of the total experimental area) was used for artificial fertilization experiments. The grassland growth in these five experimental areas differed because of the above different usages.

2.2. Sampling Strategy and Data Collection

A total of 13 sample plots of 250 m × 250 m were set up inside the five experimental sample areas. In each sample plot, a 30 m × 30 m sub-plot was set up for data collection. In each sub-plot, five quadrats (1.5 m × 1.5 m each) were set up as shown in Figure 2, taking the central point and the four corner points to represent the entire sub-plot. To reduce the artificial sampling error of biomass measurements, a strategy of 9 sub-quadrats was also used for avoiding repeat sampling in each year (Figure 2). The locations of the 13 sample plots were selected based on the following factors: (1) the growth status of the grassland was relatively uniform and was spatially representative; and (2) each sample plot (250 m × 250 m) corresponded to one MODIS pixel. However, in each sample plot, the sub-plot of 30 m× 30 m was randomly selected, corresponding to a Landsat 8 OLI and HJ-1B CCD pixel.

During the grassland growth seasons from 2013 to 2015, observations were performed approximately every 10 d, and a total of 20 field investigations was conducted over the three years. Due to factors such as weather conditions and satellite image time phases, only five times of field measurement data were matched well with satellites and were selected. They were August 2013, July 2014 and July, August and October 2015; with a total of 325 quadrat observations (see the details in Section 2.4) used to construct and analyze the biomass model.

Grass samples in each sub-quadrat (0.5 m × 0.5 m) were collected and recorded in the same way through the years. It was cut to the ground level using shears and then packed and weighed, after removing all residual litter and other non-plant materials. The sampling record included latitude, longitude, elevation, grassland types, dominant species, grassland vegetation coverage, grass height, the fresh weight and dry weight (after being dried in the lab at 64 °C until the weight remained constant).

In the field, we also collected spectral data of grassland, in a total of 18 out of the 20 field campaigns at 1170 sub-quadrats. These data were used to analyze the influence of the three filtering methods on MODIS NDVI. The spectral data were acquired using an AvaField-3 portable spectrometer (made by Holland Avantes Company). The spectral range of AvaField-3 is 300–2500 nm, particularly with a resolution of 1.4 nm and a sampling interval of 0.6 nm in the 300–1100 nm and a resolution of 15 nm and a sampling interval of 6 nm in 1100–2500 nm. The fieldwork was conducted on several consecutive sunny days. Measurements were taken on clear sunny days between 11:30 a.m. and 2:00 p.m. The fiber optic sensor, with a field view of 25°, was pointed on the target at nadir position from about a 1-m height above the ground surface. In order to increase the stability and precision of the instrument as advised in the operational manual, we preheated the spectrometer for half an hour before measurements. For each target measurement, the downwelling radiance was first measured by pointing to the white reference panel. The spectral reflectivity was directly recorded by using the upwelling radiance (i.e., reflective radiance) to the target divided by the downwelling radiance. Ten spectral reflectivity curves for each target were collected for later processing in the lab. An example of the mean spectrum of 13 sample plots on 28–29 July 2014 is shown in Figure 2. The artifacts seen (in the Figure 3) around 1150 nm could be mostly due to the mixed effect of different objects in the sample plots, including grass, soil, gravel and dead grass.

In addition, a portable GPS device for 20 ground control points (GCP) was used in the study to record the longitude, latitude and elevation, for later precise geometric corrections of the Landsat and HJ-1B satellite images.

2.3. Preprocessing of MODIS Vegetation Index Data

The MODIS vegetation index data were selected from the MODIS 16-d maximum composite NDVI vegetation index product (MOD13Q1) of the United States National Aeronautics and Space Administration. In total, 69 images spanning the years 2013–2015 taken with a spatial resolution of 250 m from orbit number h26v05 were selected. Three filtering methods, the asymmetric Gaussian function (GA), double-logistic curve (LO) and Savitzky–Golay (SG), were used to reconstruct the time series data of NDVI. The four vegetation indices were denoted as NDVI_MOD (MOD13Q1 NDVI), NDVI_GA (after GA filtering), NDVI_LO (after LO filtering) and NDVI_SG (after SG filtering). The performance of the three filtering methods was examined using the NDVI values derived from spectrum measurements in the 18 out of 20 field campaigns, as shown in

Table 1. The dates between MODIS images (h26v05) and spectral measurements in 13 sample plots in the study area.

Date of MODIS	Measurement Time
2013.08.30–09.14	2013.09.12–09.13
2014.05.26–06.10	2014.05.30–05.31
2014.06.11–06.26	2014.06.14–06.16
2014.06.27–07.12	2014.06.28–06.29
2014.07.13–07.28	2014.07.11–07.13
2014.07.13–07.28	2014.07.26–07.28
2014.08.14–08.29	2014.08.14–08.15
2014.08.30–09.14	2014.09.01–09.02
2014.09.15–09.30	2014.09.26–09.28
2014.10.17–11.01	2014.10.20–10.22
2015.05.10–05.25	2015.05.20–05.22
2015.07.13–07.28	2015.07.14–07.15
2015.07.13–07.28	2015.07.24–07.25
2015.07.29–08.13	2015.08.10–08.11
2015.08.14–08.29	2015.08.20–08.23
2015.08.30–09.14	2015.09.11–09.13
2015.10.01–10.16	2015.10.10–10.11
2015.10.17–11.01	2015.10.20–10.22

.

For the period of field sampling in August 2013, July 2014 and July, August and September 2015, the NDVI_MOD, NDVI_GA, NDVI_LO and NDVI_SG were used to construct the grassland AGB models, and the performance of the models was assessed using the field measurements.

2.4. Data Processing of Landsat 8 OLI and HJ-1B CCD and Calculation of the Vegetation Index

The Landsat 8 OLI data were obtained from the website of the United States Geological Survey (USGS). The OLI imager includes eight multi-spectral bands with a spatial resolution of 30 m and one panchromatic band with a spatial resolution of 15 m. The imaging range is 185 km × 185 km, and the revisit period is 16 d. In this study, five scenes of Landsat 8 OLI satellite images with no cloud cover in the study area in Xiahe County, Gansu Province, during 2013–2015 were downloaded. The HJ-1B CCD data were acquired from the China Center for Resource Satellite Data and Applications (http://www.cresda.com/n16/index.html). The Environment and Disaster Monitoring Satellite B (HJ-1B) carries two CCD cameras. Each CCD creates four multi-band images with a spatial resolution of 30 m; the frame width of a single CCD image is 360 km. Together, the frame width of the two CCD images spans 700 km, and the revisit period is 4 d. Based on the time of the ground surveys, the revisiting periods of the Landsat 8 and the HJ-1B satellites and cloud cover in the images, a total of five scenes of the HJ-1B CCD images and five scenes of Landsat 8 images close to the field investigation times that cover the entire study area without clouds (

Table 2. The date between satellite images and field measurements in the study area.

Date of Satellite Images	Satellite	Sensor Type	Path	Row	Sampling Time
2013.08.08	Landsat8	OLI	131	36	2013.08.06–09
2013.08.09	HJ-1B	CCD2	12	73	2013.08.06–09
2013.07.29–08.13	MODIS	Terra	26	5	2013.08.06–08.09
2014.07.26	Landsat8	OLI	131	36	2014.07.27–31
2014.07.28	HJ-1B	CCD2	13	72	2014.07.27–31
2014.07.29–08.13	MODIS	Terra	26	5	2014.07.27–07.31
2015.07.13	Landsat8	OLI	131	36	2015.07.11–17
2015.07.13	HJ-1B	CCD1	20	72	2015.07.11–17
2015.07.13–07.28	MODIS	Terra	26	5	2015.07.11–07.17
2015.08.14	Landsat8	OLI	131	36	2015.08.10–11
2015.08.12	HJ-1B	CCD2	16	72	2015.08.10–11
2015.07.29–08.13	MODIS	Terra	26	5	2015.08.10–08.11
2015.09.15	Landsat8	OLI	131	36	2015.09.14–18
2015.09.14	HJ-1B	CCD1	19	72	2015.09.14–18
2015.09.15–09.30	MODIS	Terra	26	5	2015.09.14–09.18

) were selected.

The OLI and CCD data were both processed using ENVI 5.0 software, and the Radiometric Calibration module, the Image to Image module in the registration and the FLAASH Atmospheric Correction module were used for converting the original DN value to atmospheric surface reflectance, precise geometric correction and atmospheric correction of OLI and CCD images, respectively. The image projection was defined as WGS_1984_UTM_ZONE_47N [26]. Then, the Band Math module was used to calculate the NDVI values for the Landsat 8 OLI and HJ-1B CCD images. The NDVI values were extracted from the Landsat 8 OLI and HJ-1B CCD images in the 13 sample sub-plots (250 m × 250 m). These results were used as the NDVI values (namely, NDVI_CCD and NDVI_OLI) corresponding to the ground sampling sites.

2.5. Spectral Data Processing and Accuracy Evaluation of MODIS NDVI

The spectral data were processed by using the Viewer 7 software of AvaField-3. First, we previewed the spectral curves in the Viewer and removed the abnormal curves. Then, we calculated the mean value of the normal curves for each plot. Third, according to the spectral response function of the MODIS sensor, the reflectance of the red and near-infrared bands corresponding to the MODIS sensor (near-infrared band: 841–876 nm; red band: 620–670 nm) was calculated (Equation (1)). Finally, the NDVI was calculated as Equation (2).

$${\mathsf{\rho }}_i=\frac{{{\displaystyle \sum }}_{{\mathsf{\lambda }}_{ai}}^{{\mathsf{\lambda }}_{bi}}\mathsf{\phi }\left({\mathsf{\lambda }}_i\right)\mathsf{\psi }\left({\mathsf{\lambda }}_i\right)}{{{\displaystyle \sum }}_{{\mathsf{\lambda }}_{ai}}^{{\mathsf{\lambda }}_{bi}}\mathsf{\phi }\left({\mathsf{\lambda }}_i\right)}$$

(1)

where ${\mathsf{\rho }}_i$ represents the reflectivity of band i, ${\mathsf{\lambda }}_{ai}$ represents the starting wavelength of band i, ${\mathsf{\lambda }}_{bi}$ represents the termination wavelength of band i, $\mathsf{\psi }\left({\mathsf{\lambda }}_i\right)$ represents the reflectivity value at wavelength $\mathsf{\lambda }$ and $\mathsf{\phi }\left({\mathsf{\lambda }}_i\right)$ represents the spectral response factor at wavelength $\mathsf{\lambda }$ of band i.

$$\mathrm{NDVI}=（{\mathsf{\rho }}_{Nir}-{\mathsf{\rho }}_{Red}）/（{\mathsf{\rho }}_{Nir}+{\mathsf{\rho }}_{Red}）$$

(2)

where ${\mathsf{\rho }}_{Nir}$ represents the reflectance of the near-infrared band and ${\mathsf{\rho }}_{Red}$ represents the reflectance of the red band.

According to the measured spectrum of grassland in the sub-quadrat, the NDVI values of the grassland were calculated by taking the average NDVI of the five quadrats as the ground-measured NDVI value for that sub-plot and then calculating the root mean square error (RMSE) with NDVI_MOD, NDVI_GA, NDVI_LO and NDVI_SG (Equation (3)) and their correlation coefficients (r). The smaller the RMSE and the larger the r value are, the better the vegetation index.

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\mathrm{y}}_i-{\mathrm{y}}_i^{\prime }\right)}^2}{\mathrm{n}}}$$

(3)

where ${\mathrm{y}}_i$ represents the ground-measured value, ${\mathrm{y}}_i^{\prime }$ is the estimated value, $i$ represents a sampling plot and $\mathrm{n}$ stands for the number of sample plots.

2.6. Construction of Grassland Biomass Monitoring Model and Accuracy Evaluation

The average dry weight of AGB in all of the sub-quadrats of each sub-plot was estimated as the biomass for that sample plot. Taking the grassland biomass of different sample plots as the dependent variable and taking NDVI_MOD, NDVI_GA, NDVI_LO, NDVI_SG, NDVI_CCD and NDVI_OLI as the independent variables, the leave-one-out cross validation (LOOCV) method, RMSE and relative estimate error (REE) (Equation (4)) were used to evaluate the accuracy and estimation errors of the linear, exponential, logarithmic and power regression models for the 6 vegetation indices at the sample plot level.

$$\mathrm{REE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left[\left({\mathrm{y}}_i-{\mathrm{y}}_i^{\prime }\right)/{\mathrm{y}}_i^{\prime }\right]}^2}{\mathrm{n}}}$$

(4)

where ${\mathrm{y}}_i$ represents the ground-measured value, ${\mathrm{y}}_i^{\prime }$ is the estimated value, $i$ represents a sampling plot and $\mathrm{n}$ stands for the number of sample plots.

In addition, to study the relative estimation error of the alpine meadow grassland AGB monitoring model on the regional scale (i.e., the experimental area of this study), the average AGB dry weight for all of the sample plots in each sample area was adopted as the ground-measured biomass in this sample area. Then, the RMSE and REE were used to assess the performance of each model for estimating total yield and yield per unit area.

3. Results and Analysis

3.1. Statistical Analysis of Ground Observation AGB and the Corresponding Multi-Source Satellite NDVI

Table 3. Descriptive statistics of grassland above-ground biomass (AGB) and corresponding multi-source satellite NDVI during July–September of 2013–2015 in the study area, Xiahe County.

Index	Statistics	Plot
		E1	E2	E3	E4	E5	E6	E7	D8	D9	C10	B11	A12	A13	All
Biomass (103 kg/ha)	Maximum	2.20	2.28	2.70	2.08	2.52	2.94	3.20	2.47	1.77	2.03	2.41	3.96	2.67	3.96
	Minimum	1.12	1.15	0.83	0.86	0.79	1.13	0.97	0.94	1.18	0.95	0.93	0.83	0.75	0.75
	Average	1.77	1.89	1.86	1.29	1.48	1.71	1.75	1.82	1.51	1.37	1.44	1.68	1.28	1.81
	Standard deviation	0.47	0.52	0.81	0.534	0.74	0.84	1.00	0.68	0.25	0.46	0.67	1.52	0.93	0.85
	Cv	0.27	0.27	0.44	0.42	0.50	0.50	0.57	0.37	0.16	0.34	0.47	0.12	0.73	0.47
	n	25	25	25	25	25	25	25	25	25	25	25	25	25	325
HJ-CCD NDVI	Maximum	0.68	0.74	0.80	0.76	0.80	0.78	0.73	0.74	0.68	0.71	0.71	0.67	0.67	0.80
	Minimum	0.57	0.64	0.67	0.62	0.68	0.64	0.65	0.65	0.62	0.58	0.57	0.53	0.51	0.51
	Average	0.65	0.71	0.75	0.67	0.73	0.72	0.70	0.70	0.65	0.64	0.64	0.61	0.59	0.68
	Standard deviation	0.04	0.05	0.05	0.07	0.05	0.06	0.04	0.04	0.03	0.05	0.06	0.06	0.07	0.07
	Cv	0.07	0.07	0.06	0.10	0.07	0.08	0.06	0.05	0.05	0.08	0.09	0.09	0.12	0.10
	n	5	5	5	5	5	5	5	5	5	5	5	5	5	65
Landsat-8 OLI DVI	Maximum	0.77	0.86	0.90	0.86	0.87	0.86	0.85	0.87	0.82	0.78	0.82	0.79	0.71	0.97
	Minimum	0.60	0.70	0.74	0.65	0.73	0.73	0.75	0.76	0.63	0.66	0.69	0.64	0.55	0.55
	Average	0.70	0.79	0.84	0.755	0.81	0.80	0.81	0.82	0.72	0.72	0.75	0.70	0.62	0.78
	Standard deviation	0.07	0.07	0.07	0.10	0.06	0.06	0.05	0.05	0.08	0.06	0.06	0.07	0.07	0.10
	Cv	0.10	0.09	0.08	0.13	0.08	0.08	0.06	0.06	0.11	0.08	0.08	0.10	0.11	0.13
	n	5	5	5	5	5	5	5	5	5	5	5	5	5	65
MOD13Q1 NDVI	Maximum	0.82	0.84	0.82	0.82	0.81	0.81	0.82	0.79	0.73	0.75	0.79	0.69	0.67	0.84
	Minimum	0.66	0.69	0.68	0.65	0.61	0.64	0.64	0.58	0.62	0.46	0.49	0.49	0.44	0.44
	Average	0.73	0.74	0.73	0.72	0.70	0.71	0.72	0.69	0.68	0.60	0.60	0.57	0.58	0.69
	Standard deviation	0.06	0.07	0.07	0.07	0.08	0.07	0.08	0.08	0.05	0.12	0.13	0.09	0.10	0.10
	Cv	0.09	0.09	0.09	0.10	0.12	0.10	0.11	0.12	0.07	0.20	0.21	0.15	0.18	0.15
	n	5	5	5	5	5	5	5	5	5	5	5	5	5	65

Note: Cv depicts the coefficient of variation, equal to the standard deviation divided by the mean; n depicts the number of observations.

shows the statistical results of AGB in the surveyed sample plots and the NDVI calculated using the corresponding remote sensing data from 2013–2015. There are considerable differences in the biomass during the grass growing season in the 13 sample plots of the experimental area. The average values range from 1280–1887 kg DW/ha over the past three years, while the coefficient of variation ranges from 0.16–0.73. The largest biomass from the 13 sample plots is 3963 kg DW/ha, and the minimum is 745.5 kg DW/ha. The magnitude of the variation in the NDVI_CCD, NDVI_OLI and NDVI_MOD values is relatively small, and their dispersion degree is similar; the standard deviation and coefficient of variation range from 0.03–0.13 and 0.05–0.18, respectively. In terms of mean NDVI value ranges, they are very similar, from 0.59–0.75 and 0.57–0.74, respectively, for NDVI_CCD and NDVI_MOD and 0.62–0.83 for NDVI_OLI, slightly larger than the former two.

For the entire experimental area, the average total biomass in August 2013–2015 is 285 × 10³ kg, and the average grass yield is 1813.5 kg DW/ha. In the experimental area, the average biomass of the five times in July–September of 2013–2015 is 1808.4 kg DW/ha, with the standard deviation of 847.3 kg DW/ha and the variation coefficient of 0.47. Similar to the sample plots level, the NDVI_CCD and NDVI_MOD means are very close (0.68–0.69), slightly smaller than the 0.78 of NDVI_OLI.

3.2. Influence of Different Filtering Methods on MODIS NDVI

The NDVI calculated using the measured spectrum data on the ground is treated as the ground-measured value to examine the accuracy of NDVI_MOD and the three NDVI time series after the filters are applied, i.e., NDVI_GA, NDVI_LO and NDVI_SG. The results are shown in

Table 4. The result of NDVI_MOD, NDVI_GA, NDVI_LO and NDVI_SG as compared with the field measurements in 13 sample plots of the study area.

Vegetation Index	Index	Plot
		E1	E2	E3	E4	E5	E6	E7	D8	DB9	C10	B11	A12	A13	Average
NDVIMOD	RMSE	0.035	0.083	0.081	0.079	0.087	0.089	0.105	0.115	0.110	0.087	0.051	0.120	0.111	0.102
	R	0.83 **	0.85 **	0.86 **	0.83 **	0.81 **	0.81 **	0.78 **	0.70	0.74 **	0.82 **	0.87 **	0.82	0.77 **	0.77 **
NDVIGA	RMSE	0.029	0.078	0.077	0.062	0.077	0.070	0.087	0.094	0.094	0.079	0.051	0.113	0.107	0.091
	R	0.90 **	0.89 **	0.88 **	0.90 **	0.86 **	0.89 **	0.86 **	0.82 *	0.81 **	0.86 **	0.89 **	0.85 **	0.79 **	0.83 **
NDVILO	RMSE	0.028	0.078	0.076	0.062	0.077	0.072	0.088	0.095	0.095	0.080	0.050	0.112	0.107	0.091
	R	0.90 **	0.89 **	0.88 **	0.90 **	0.86 **	0.88 **	0.86 **	0.81 *	0.80 **	0.86 **	0.89 **	0.86 **	0.78 **	0.83 **
NDVISG	RMSE	0.028	0.079	0.080	0.060	0.075	0.070	0.085	0.093	0.095	0.078	0.051	0.098	0.107	0.090
	R	0.91 **	0.88 **	0.87 **	0.91 **	0.87 *	0.89 **	0.86 **	0.82 *	0.81 **	0.88 **	0.88 **	0.91 **	0.78 **	0.84 **

** p < 0.01, * p < 0.05.

, indicating that the three filtering methods perform better overall than the original MOD13Q1 NDVI, and the SG filtering is the best. At the sample plot level, the SG filtering method achieves the best results for Sample Plots 1, 4–8, 10, 12 and 13, with RMSE ranging from 0.028–0.098 and r values ranging from 0.78–0.91. The logistic filtering method achieved the best results for Sample Plots 3 and 11; with RMSE of 0.050–0.076 and r values of 0.88–0.89. The Gaussian filtering method achieved the best results for Sample Plots 2 and 9; with RMSE of 0.078–0.094 and r values of 0.81–0.89. Among all of the sample plots, all three filtering methods achieved relatively good results at Plot 1, with RMSE of 0.028–0.029 and r of 0.90–0.91. In contrast, all three filtering methods achieve relatively poor results at Plot 13, with RMSE values ranging from 0.098–0.113 and r values from 0.85–0.91.

3.3. Grassland Biomass Monitoring Model in the Study Area and Evaluation of Its Accuracy at the Sample Plot Level

The results of the accuracy evaluation validated by LOOCV for the grassland biomass models constructed on the three satellites, NDVI_MOD, NDVI_CCD and NDVI_OLI, as well as the three filtering methods applied to MODIS NDVI only are listed in

Table 5. The validation results by leave-one-out cross validation for the grassland biomass models based on multi-source satellite data. REE, relative estimate error.

Vegetation Index	Model	Accuracy Evaluation
		RMSE (kg/ha)	REE (%)
NDVIMOD	Linear	594.5	47.8
	Exponential	574.6	35.3
	Logarithm	619.2	61.4
	Power	598.8	36.7
NDVISG	Linear	573.3	45.0
	Exponential	549.7	33.9
	Logarithm	594.6	57.7
	Power	571.4	35.0
NDVILO	Linear	581.1	44.1
	Exponential	560.4	34.1
	Logarithm	602.4	53.8
	Power	582.0	35.2
NDVIGA	Linear	583.8	44.2
	Exponential	562.0	34.1
	Logarithm	605.8	53.9
	Power	585.2	35.3
NDVICCD	Linear	557.1	45.2
	Exponential	548.4	31.6
	Logarithm	567.6	74.0
	Power	552.6	32.1
NDVIOLI	Linear	516.8	33.7
	Exponential	511.6	29.1
	Logarithm	528.9	41.2
	Power	512.1	29.4

. It is found that in the four types of grassland AGB estimation models (linear, logarithmic, power and exponential), the exponential model performs the best, with the smallest RMSE and REE, while the linear or power model performs the second best. Among the six vegetation indices, the exponential model based on NDVI_OLI is the best, with the smallest RMSE of 511.6 kg DW/ha and REE of 29.1% and the highest estimation accuracy of the AGB (Figure 4d). The exponential model based on NDVI_CCD is the second best (Figure 4c). The exponential model based on NDVI_MOD is worse (Figure 4a). All of the filtering methods perform better than the NDVI_MOD, with NDVI_SG the best (Figure 4b).

Table 6. The results of model fits with the optimum inversion models based on multi-source satellite data.

Vegetation Index	Parameter Estimation and T Test			Regression Significance Test
	Parameter	Estimated Value	T	R2	F
NDVIMOD	b	3.0149	6.817 **	0.46	46.478 **
	a	193.442	3.229 **
NDVISG	b	3.4496	7.400 **	0.50	54.759 **
	a	144.265	3.078 **
NDVILO	b	3.496	7.030 **	0.47	49.413 **
	a	140.404	2.888 **
NDVIGA	b	3.487	6.889 **	0.46	47.464 **
	a	141.383	2.839 **
NDVICCD	b	5.0715	8.581 **	0.57	73.634 **
	a	48.325	2.457 **
NDVIOLI	b	3.6787	10.017 **	0.65	100.341 **
	a	85.916	3.427 **

Note: ** represents p < 0.01; a and b represent the constant and exponential term of the model, respectively; T and F are the significant values according to the T and F tests.

shows the parameters of the best fit model for each of the six vegetation indices. As shown, all six models pass the T test and F test at a significance level of 0.01, with NDVI_OLI the best, followed by NDVI_CCD, NDVI_SG and NDVI_MOD.

4. Discussion

4.1. Influence of Different Remote Sensing Data on the Estimation Error of Grassland Biomass

In this study, NDVI_MOD, NDVI_CCD and NDVI_OLI data were used in combination with ground observation data during 2013–2015 to construct a grassland biomass inversion model. The errors based on NDVI_MOD, NDVI_CCD and NDVI_OLI reached 35.3%, 31.6% and 29.1%, respectively, at the sample plot level (Table 5). Our results indicate that the yield per unit area estimated using the exponential model based on NDVI_OLI (1518 kg DW/ha) is closest to the ground-measured value, and its estimation error is the lowest (30.7%), followed by NDVI_CCD (32.4%) and NDVI_MOD (39.6%) (

Table 7. The estimation error of the grassland AGB evaluated by different remote sensing data and filtering methods in study area.

Data Type	Vegetation Index	Formula	Yield (AGB) (kg DW/ha)	REE (%)
Different remote sensing data	NDVIOLI	y = 85.916e3.6787x	1518.0	30.7
	NDVICCD	y = 48.325e5.0715x	1472.7	32.4
	NDVIMOD	y = 193.442e3.0149x	1431.0	39.6
Different filtering methods	NDVISG	y = 132.146e3.584x	1564.1	34.9
	NDVILO	y = 140.404e3.496x	1422.3	38.6
	NDVIGA	y = 141.383e3.478x	1408.1	39.3

). This result suggests that the inverse models based on 30-m resolution (namely, NDVI_CCD and NDVI_OLI) are better than that based on MODIS data with 250-m resolution, and the model based on NDVI_OLI is better than that on NDVI_CCD.

There are multiple reasons for the above results. Existing studies indicate that the accuracy of the grassland biomass inversion model is not only affected by the temporal and spatial resolutions of the satellite images, sensor type and method of image processing [18,27], but is also subject to the influence of the size, number and representativeness of ground sampling sites. In this study, the five-point method is adopted to determine the sub-plot (30 m × 30 m), which matches the spatial resolution of the Landsat 8 OLI and HJ-1B CCD images. Although the Landsat or HJ-1B pixel could not always exactly match a sub-plot (30 m × 30 m) on the ground due to the orbit drift, we consider this a minor problem, since the relatively uniform distribution of grass in each 250 m × 250 m sample plot was one of the selection criteria. A slight mismatch between sub-plot and satellite pixel does not really make much difference in our case. The average grassland biomass based on five quadrats can reasonably represent the variation of grassland biomass within a full 30 m × 30 m sub-plot; therefore, the accuracy is relatively high for the grassland biomass monitoring model constructed using the remote sensing image data with 30-m resolution (namely, NDVI_OLI and NDVI_CCD). However, the spatial matching with the MODIS vegetation index, which uses a 250-m spatial resolution, has some limitations; consequently, the accuracy of the simulation model is relatively low. In addition, there is an obvious difference between the estimation errors of the biomass models based on NDVI_OLI and NDVI_CCD, which have the same spatial resolution (30 m). The error of the AGB estimation model based on NDVI_OLI is smaller than NDVI_CCD. This difference occurred mainly because the near-infrared (NIR) and red bands (RED) used to calculate their NDVI values come from different sensors that have different band ranges. Specifically, the NIR and RED band ranges used to calculate NDVI_CCD are 0.76–0.90 μm and 0.63–0.69 μm, respectively, while for NDVI_OLI, they are 0.845–0.885 μm and 0.630–0.680 μm, respectively. The narrow band receives less signal disturbances. In addition, the OLI NIR also excludes the influence of water vapor absorption at 0.825 μm, and therefore, the OLI NIR’s contribution to retrieve grassland information is more prominent [28].

4.2. Influence of Three Filtering Methods on the Error of Grassland AGB Estimation Based on MODIS NDVI

The exponential models (based on the accuracy of the best inversion mode) for grassland AGB of the three filtering methods including NDVI_SG, NDVI_LO and NDVI_GA, can better simulate the grassland AGB than the NDVI_MOD, in which the error decreased by approximately 1.40, 1.14 and 1.13%, respectively (Table 5). The yield per unit area for grassland AGB estimated by NDVI_SG (1564.1 kg DW/ha) is the closest to the ground-measured value, and its estimation error was the lowest (34.9%), followed by NDVI_LO and NDVI_GA (Table 7). This result indicates that the SG filtering method is better at eliminating abnormal values from the MOD13Q1 NDVI; consequently, it yields biomass estimates closer to the ground-measured value. This finding is similar to the research results of Chen et al. (2004) and Geng et al. (2014) Their studies in Southeast Asia indicated that SG filtering can partially eliminate atmospheric interference and the influence of mixed pixels on the grassland vegetation index and can improve the representativeness of MODIS NDVI for grassland vegetation [23,29].

4.3. Assessment of Previously-Established Biomass Inversion Models Based on the MODIS Vegetation Index over the Tibetan Plateau

Many scholars have used the MODIS vegetation index to conduct numerous successful studies on grassland biomass in the Tibetan Plateau area. However, the existing studies show considerable differences in model errors for alpine grassland in the Tibetan Plateau area at the sample plot scale (

Table 8. The existing AGB estimation models for alpine grassland in the Tibetan Plateau region.

Model	Study Area	Area (104 ha)	MODIS	Formula	R2	Literature
I	Tibetan Plateau	25,724	NDVI	y = 225.42 × e4.4368x	0.75	[19]
II	Tibetan Plateau	25,724	NDVI	y = 268.810 × e2.398x	0.49	[30]
III	The northeast of Tibetan Plateau (Gannan Prefecture)	380	EVI	y = 3738.073x1.553	0.63	[31]
IV	The northeast of Tibetan Plateau (Gannan Prefecture)	380	EVI	y = 5320.7x1.9776	0.62	[32]
V	The northwest of Gannan Prefecture	62.74	EVI	y = 1719.1x2.2588	0.63	[33]
	(Xiahe County)

). For example, the grassland biomass inversion model constructed by Xu et al. (2007) used field measurement data in the Tibetan Plateau area and the 10-d maximum composite MODIS NDVI (with a 1-km resolution) from the end of July to the end of September in 2007, to construct an exponential function (Model I), with R² of 0.75 and overall estimation accuracy up to 80% [19]. In the grassland biomass regression models constructed by Feng et al. (2011), which combined measurement data from the Tibetan Plateau during July–August of 2005 and 2006 with the monthly maximum composition 1-km resolution MODIS NDVI, the exponential function model is the optimum inversion model (Model II). However, its accuracy is lower than that of Model I; with R² of 0.49 and RMSE up to 671.8 kg DW/ha [30]. Cui et al. (2011), Wang et al. (2010) and Bao et al. (2010) used the MODIS vegetation index at 500-m, 1-km and 250-m resolution, respectively, to construct biomass regression models for the Gannan prefecture. Their results are similar. The optimum inversion model (in the order of Models III, IV and V) is a power function based on MODIS EVI, with the R² values of 0.63, 0.62 and 0.63, respectively [31,32,33]. Clearly, there are large differences among models. Here, we intend to examine the feasibility of these models all developed at large scales to be applied to our small study area with detailed field data.

Table 9. The estimation error for the grassland AGB based on existing inversion models at the sample plot and the entire study area scales.

Model	Sample Plot Scale	Study Area Scale
	REE (%)	Yield (kg DW/ha)	REE (%)
I	71.7	5748.8	68.9
II	58.6	1352.9	48.3
III	47.8	1397.8	46.2
IV	47.2	1470.8	48.0
V	47.3	1551.3	44.6

shows the results of statistical analysis using the biomass data measured on the ground surface during 2013–2015 in our study area to examine the error of the aforementioned five models. As shown, the Model IV error is the smallest (47.2%), followed by V, III, II and I. Among the five models, Models I and II based on the entire Tibetan Plateau perform the worse. Although Model I achieved the lowest error in the Tibetan Plateau (Table 8), its estimation error of grassland AGB in the study area is highest (71.7%). Models III, V and IV are built on the regional area (Xiahe and Gannan) of the Tibetan Plateau; their estimation errors are much smaller, around 47%.

By analyzing the differences in biomass estimation among Models I–V in the study area (Table 8), we can see that the yield per unit area of grassland biomass (1551.3 kg DW/ha) estimated by Model V in Xiahe County is closest to the ground-measured value of the grassland AGB in the experimental area, and its estimation error for yield per unit area is the lowest among the five models (44.6%). The yield per unit area of grassland AGB estimated by Model I for the experimental area (5748.8 kg DW/ha) is much larger than the ground-measured value, and its estimation error is also the largest (68.9%); the error of the yield per unit area estimated by the other three models ranked from small to large as III, IV and II.

Compared with this study, either at the sample plot scale or the regional scale (namely, the study area of Xiahe County), the estimation accuracy of grassland AGB based on NDVI_MOD is the highest, and its estimation error for the yield of grassland AGB per unit area is the smallest in the experimental area (39.3%) (Table 7). Although the previous inversion Models I–V can reflect the overall variation trend of grassland biomass, their errors are higher when applied to our study area. The reason is that in large areas (the Tibetan Plateau area or an entire prefecture or county), the types of natural grassland pasture are complicated, the geographical distribution is wide, the spatial heterogeneity is strong and the vegetation index value is subject to influences from many features and factors (e.g., the ecological environment of the grassland). Therefore, the models built based on a large scale, when applied to small local areas (such as this experimental area), would result in large errors and low accuracy, poor stability and large spatial variations. This conclusion is similar to the results of other studies in other areas [16,34,35,36,37].

4.4. Limitations and Prospects of Remote Sensing Monitoring Biomass

In this study, although the optimum models for the remote sensing monitoring of grassland biomass in the study area are determined, due to the limitations of factors, such as the duration of sampling on the ground (only the growing seasons of the grassland from 2013–2015 were sampled), there is still some error and uncertainty for these inversion models, especially the grassland biomass inversion model based on the MODIS NDVI data. There are clearly two areas that we may improve in the future. One is the NDVI saturation problem, as is clearly seen in Figure 4 when the biomass is larger than 3000 kg DW/ha. One way to solve this problem is to use the wide dynamic range vegetation index (WDRVI) [38,39,40], which was developed mostly to simulate the biomass for crops, usually with very high biomass. It is worthy to introduce WDRVI for the grassland biomass study in the future, especially to explore its sensitivity to the high end of grassland biomass. The other one is the bidirectional reflectance distribution function (BRDF) effects on all remotely-sensed NDVI. As the BRDF effects are already considered in the MODIS NDVI products, our derived NDVI from the Landsat and HJ-1B are not. Although it is not the intention of this study to evaluate the methods for deriving BRDF from these two satellites, we notice that quite a few studies have already reported their methods [41,42,43,44]. For example, one study [45] found that there are 0.02–0.06 differences for Landsat images between regular reflectivity and BRDF, and it is worthy to explore and use BRDF-derived NDVI for modeling purposes in our future efforts.

The area of grassland resources is distributed broadly worldwide (accounting for approximately 40% of the global land biome) [2], and the variation on the temporal and spatial scales is dramatic [46]. Satellite data using the relatively high resolutions of Landsat MSS, TM and SPOT show relatively high accuracy in the monitoring of grassland AGB over extended periods; however, because of influences from the transit period, frame width and cloud cover, it is difficult to obtain long time series of high-quality images. Although MODIS has relatively low spatial resolution and relatively large estimation error of grassland AGB, it is still particularly suitable for monitoring widely distributed grassland AGB, mostly due to its high temporal resolution (daily) and large spatial coverage (2330 km). Therefore, it is important to explore new research approaches to grassland AGB monitoring based on MODIS data to improve the accuracy of grassland remote sensing inversion over large regions in the future. These approaches mainly involve the following aspects: (1) enhance the spatial representativeness of ground sampling sites, for example by increasing their number and area and improving the range observed by the ground sampling site and the corresponding spatial matching problem given the size of satellite image pixels; (2) improve the temporal matching between ground sampling sites and remote sensing data, for example by better scheduling the times of field investigations to reduce the time differences between ground surveys and satellite image acquisitions; (3) incorporate new remote sensing observation techniques (e.g., hyperspectral imagery and the UAV remote sensing technique) and strengthen research on the spectral characteristics of the grassland vegetation community and the applications of the narrow band remote sensing vegetation index in monitoring grassland AGB [47,48]; and (4) construct multi-factor grassland AGB estimation models based on statistical analysis and machine learning techniques. These multiple factors include climatic factors (e.g., sunlight, temperature and rainfall), soil factors (e.g., soil nutrients, soil structure and fertility), biological factors (e.g., grassland type, species richness and distribution of malignant weeds) and management factors (e.g., pasture, fencing enclosures and rotational grazing) [49,50,51]. For example, Li et al. (2013) used neural networks to build an AGB model based on multiple MODIS vegetation indices and showed higher accuracy than a model based on a single index, by decreased RMSE of 433 kg/ha and increased R² of 0.35 [49]. A multi-factor model by Liang et al. (2016) showed decreased RMSE by 14.5% as compared with the optimum single-factor model [50]. Diouf et al. (2016) studied the semi-arid grassland in the Sahel region and indicated that a combined photosynthetic radiation and meteorological data model had better performance (R² = 0.69 and RMSE = 483 kg DW/ha) than the single-factor model of photosynthetic radiation or meteorological data (R² = 0.63 and 0.55 and RMSE = 550 kg DW/ha and 585 kg DW/ha, respectively) [51].

5. Conclusions

In this study, based on MOD13Q1, HJ-1B CCD and Landsat 8 OLI remote sensing data, grassland observation data in the Sangke grassland of Xiahe County during 2013–2015 are combined to construct a grassland AGB estimation model based on different remote sensing data, and the influence of different filtering approaches for MODIS NDVI on the biomass estimation error of alpine meadow grassland is investigated. The simulation errors of several grassland AGB models in the study area are compared and analyzed, and the applicability of the models is evaluated. The following primary conclusions have been reached.

(1)

There is a significant difference in the estimation errors of alpine meadow grassland AGB using remote sensing data from the Chinese HJ-1B CCD, Terra MODIS and Landsat 8 OLI. In this study, the grassland AGB optimum inversion model of the experimental area is the exponential model based on NDVI_MOD, NDVI_OLI and NDVI_CCD, but different models show considerable differences in the error of grassland AGB inversion. The errors for the estimation of grassland AGB for the optimum models based on NDVI_MOD, NDVI_CCD and NDVI_OLI at the sample plot level are 35.3%, 31.6% and 29.1%, respectively. Their yield per unit area estimations for grassland AGB in the experimental area indicate that the exponential model based on NDVI_OLI yielded values closest to the ground-measured value; its estimation error for yield per unit area is the smallest (30.7%). The estimation error for yield per unit area for the experimental area with the optimum AGB inversion model based on NDVI_OLI decrease by eight and two percentage points, respectively, compared to the optimum inversion models based on NDVI_MOD and NDVI_CCD.

(2)

The filtering and de-noising processing of MOD13Q1 NDVI are key for reducing the AGB inversion error of alpine meadow grassland based on MODIS data. At the sample plot level, the estimation errors of the AGB estimation models based on NDVI_SG, NDVI_LO and NDVI_GA decreased by 1.40, 1.14 and 1.13 percentage points, respectively, compared to the AGB estimation model based on NDVI_MOD. On the study area scale (161.36 ha), the estimation errors for the yield per unit area of grassland AGB based on NDVI_SG, NDVI_LO and NDVI_GA decreased by 4.48, 0.95 and 0.22, respectively, compared to that based on NDVI_MOD.

(3)

The feasibility study on previous models (I and II, III and IV and V) developed (on MODIS indices) at broad scales to apply to our small study area suggests that the estimation error of these models is higher than that of the NDVI_MOD model constructed in this study by 11.9%–36.4% at the sample plot scale and 5.3%–29.6% at the study area scale. Models V, IV and III based on Xiahe County and Gannan Prefecture do not show considerable difference on the estimation error of AGB, ranging from 47.2%–47.8% at the sample plot level and 44.6%–48.0% of the yield per unit area at the study area level. However, Models I and II based on the Tibetan Plateau scale show much larger estimation error, up to 71.7% and 58.6%, respectively, at the sample plot scale and 68.9% and 48.3% of the yield per unit area at the study area scale.