Improved Estimation of Aboveground Biomass of Disturbed Grassland through Including Bare Ground and Grazing Intensity

: Accurate approaches to aboveground biomass (AGB) estimation are required to support appraisal of the effectiveness of land use measures, which seek to protect grazing-adapted grasslands atop the Qinghai-Tibet Plateau (QTP). This methodological study assesses the effectiveness of one commonly used visible band vegetation index, Red Green Blue Vegetation Index (RGBVI), obtained from unmanned aerial vehicle (UAV), in estimating AGB timely and accurately at the local scale, seeking to improve the estimation accuracy by taking into account in situ collected information on disturbed grassland. Particular emphasis is placed upon the mapping and quantiﬁcation of areas disturbed by grazing (simulated via mowing) and plateau pika ( Ochotona curzoniae ) that have led to the emergence of bare ground. The initial model involving only RGBVI performed poorly in AGB estimation by underestimating high AGB by around 10% and overestimating low AGB by about 10%. The estimation model was modiﬁed by the mowing intensity ratio and bare ground metrics. The former almost doubled the estimation accuracy from R 2 = 0.44 to 0.81. However, this modiﬁcation caused the bare ground AGB to be overestimated by about 38 and 19 g m − 2 for 2018 and 2019, respectively. Although further modiﬁcation of the model by bare ground metrics improved the accuracy slightly to 0.88, it markedly reduced the overestimation of low AGB values. It is recommended that grazing intensity be incorporated into the micro-scale estimation of AGB, together with the bare ground modiﬁcation metrics, especially for severely disturbed meadows with a sizable portion of bare ground.


Introduction
Grasslands are vital ecosystems that occupy about 40% of the Earth's surface [1].
With an area of about 2.5 million km 2 , grasslands atop the Qinghai-Tibet Plateau (QTP) provide an important grazing resource that supports the livelihood of local pastoralists [2]. Maintenance of grazing-adapted ecosystems is an integral part of long-term environmental protection programmes in this region [3]. In recent decades, climate change, overgrazing and rodents threaten this ecosystem [4,5]. Around 40% of grassland areas in the Sanjiangyuan region of the QTP have been subjected to fragmentation and degradation [6][7][8]. Essentially, the prospect for sustainable practices reflects the balance between changes to grassland productivity, on the one hand (i.e., net primary productivity; NPP), and impacts of land use practices, especially stocking rates, on the other [9,10]. Given the vast expanse of the landscapes and ecosystems in this region, research to date has emphasized the use of remotely sensed applications to appraise spatial and temporal variability in the patterns and increase their value [62,63], thereby resulting in grassland AGB to be overestimated. At present, it remains unknown how the estimation accuracy can be improved via considering the effect of bare ground caused by pika disturbance.
In order to fill the identified knowledge gaps, this study aims to estimate alpine meadow AGB from multi-temporal drone images at a micro-scale and improve estimation accuracy in relation to two types of external disturbances (mowing-simulated grazing and rodents). The specific objectives are (a) to develop a method for alpine grassland AGB estimation from fine-resolution drone images without radiometric calibration; (b) to evaluate the improvement in AGB estimation accuracy by taking into account the meadow conditions caused by disturbances; and (c) to evaluate the influence of bare ground on the AGB estimation accuracy.

Study Area
The study area is located in the Henan Mongolian Autonomous County in the eastern Sanjiangyuan Region on the QTP (Figure 1). At an average elevation of 3600 m above sea level, this area has a plateau continental climate, with annual temperature ranging from 9.2 to 14.6 • C, and annual precipitation ranging from 597.1 to 615.5 mm. There are 5982.29 km 2 of available grassland in this county, 92.68% of which takes the form of typical alpine meadow. The grassland vegetation is dominated by Kobresia tibetica and has been affected by rodents and overgrazing [64,65], putting this area at the risk of degradation. The study site is located on a floodplain of the Yellow River with a gentle topography, and it suffers a moderate disturbance of pika (Ochotona curzoniae) and grazing.
Remote Sens. 2021, 13, x FOR PEER REVIEW and gnawing [58,59]. Such bare ground exerts significant influences on the spectra erties of vegetation canopies of arid grassland [60,61], thus degrading the estima curacy. Because VIs are sensitive to background factors, the soil background bri could increase their value [62,63], thereby resulting in grassland AGB to be overest At present, it remains unknown how the estimation accuracy can be improved v sidering the effect of bare ground caused by pika disturbance.
In order to fill the identified knowledge gaps, this study aims to estimate meadow AGB from multi-temporal drone images at a micro-scale and improve est accuracy in relation to two types of external disturbances (mowing-simulated graz rodents). The specific objectives are (a) to develop a method for alpine grasslan estimation from fine-resolution drone images without radiometric calibration; (b) uate the improvement in AGB estimation accuracy by taking into account the m conditions caused by disturbances; and (c) to evaluate the influence of bare ground AGB estimation accuracy.

Study Area
The study area is located in the Henan Mongolian Autonomous County in the Sanjiangyuan Region on the QTP (Figure 1). At an average elevation of 3600 m ab level, this area has a plateau continental climate, with annual temperature rangin 9.2 to 14.6 °C, and annual precipitation ranging from 597.1 to 615.5 mm. There are km 2 of available grassland in this county, 92.68% of which takes the form of typica meadow. The grassland vegetation is dominated by Kobresia tibetica and has been a by rodents and overgrazing [64,65], putting this area at the risk of degradation. Th site is located on a floodplain of the Yellow River with a gentle topography, and it a moderate disturbance of pika (Ochotona curzoniae) and grazing.

Setup of Experimental Plots
In total, 27 plots of 25 by 30 m were set up with a corridor of 5 m between neighbouring plots in a fenced site in June 2018 ( Figure 1). A wire mesh fence was set 50 cm above the ground and extended about 50 cm beneath the surface, deeper than the approximate depth (20 to 30 cm) of pika burrowing so as to prevent pika intrusion into and roaming among the plots ( Figure 2) [66,67]. All the plots were treated with three levels of disturbance (none, medium and high) for pika population and mowing intensity individually and jointly, resulting in nine treatments in total (Table 1). Each treatment was replicated three times. Spatially, the nine treatments of a given replica were randomly allocated to each of the three rows of plots ( Figure 1). In total, 27 plots of 25 by 30 m were set up with a corridor of 5 m between neighbouring plots in a fenced site in June 2018 ( Figure 1). A wire mesh fence was set 50 cm above the ground and extended about 50 cm beneath the surface, deeper than the approximate depth (20 to 30 cm) of pika burrowing so as to prevent pika intrusion into and roaming among the plots ( Figure 2) [66,67]. All the plots were treated with three levels of disturbance (none, medium and high) for pika population and mowing intensity individually and jointly, resulting in nine treatments in total (Table 1). Each treatment was replicated three times. Spatially, the nine treatments of a given replica were randomly allocated to each of the three rows of plots ( Figure 1).

Figure 2.
Images of the study site showing the experimental area with its weather station (a) and the above/below ground characteristics of the plot fencing (b). Plateau pika used in this study is a small mammal with a body weight of 130-200 g, widely scattered in the alpine grassland of the QTP [68]. It affects grassland not only through directly devouring forage but also via burrowing and digging [69,70]. The intensity of pika disturbance was set with reference to the literature, namely, 200 individual ha −1 or 14 in a plot for strong disturbance (Ph) and 7 for medium disturbance (Pm) [71,72]. No pika (Pn) plots served as the reference. In each plot, traps were placed to eradicate the existing pika population one month prior to field sampling. Afterwards, poisoning was applied to each plot to exterminate any residual pika. More traps were placed outside the plots to collect pika before the experiment. Later, they were released to the pre-designated plots in accordance with the pre-determined intensity at a male-female ratio of 1:1 for Ph and 3:4 for Pm [73]. The pika population and gender in each plot were re-assessed in July  Plateau pika used in this study is a small mammal with a body weight of 130-200 g, widely scattered in the alpine grassland of the QTP [68]. It affects grassland not only through directly devouring forage but also via burrowing and digging [69,70]. The intensity of pika disturbance was set with reference to the literature, namely, 200 individual ha −1 or 14 in a plot for strong disturbance (Ph) and 7 for medium disturbance (Pm) [71,72]. No pika (Pn) plots served as the reference. In each plot, traps were placed to eradicate the existing pika population one month prior to field sampling. Afterwards, poisoning was applied to each plot to exterminate any residual pika. More traps were placed outside the plots to collect pika before the experiment. Later, they were released to the pre-designated plots in accordance with the pre-determined intensity at a male-female ratio of 1:1 for Ph and 3:4 for Pm [73]. The pika population and gender in each plot were re-assessed in July 2019 to ensure that they still conformed to the required specifications. Any deviations from them were remedied by removing or supplementing pika to the affected plots as appropriate.
Livestock grazing was simulated via mowing twice a year in mid-July and mid-August. Mowing height could be used as a proxy for grazing intensity [74], as the same grazed pastures had a similar vegetation height [57]. So, mowing intensity was dictated by the mower blade height. Heavy grazing (Gh) was simulated by setting the blade at the extremely low level of 2 cm above the ground; the height was 12 cm for moderate grazing (Gm), which was about half of the plant height (24 cm). Gn was the reference treatment without mowing. The exact blade height was based on the one-third principle, namely, the clipped grass should not exceed one third of the plant height to reduce stress on the mown grass and to promote forage lateral growth [75].

UAV Images and Biomass Sampling
All experimental plots were photographed using an RGB digital camera (1 inch, 20 megapixel Complementary Metal Oxide Semiconductor sensor) aboard a Dajiang Phantom 4 UAV (Shenzhen Dajiang Baiwang Technology Co., Ltd., Shenzhen, China) in the auto mode in late August of 2018 and 2019. Acquired at a flying height of around 40 m above the ground, each drone image of a plot comprises 4864 × 3648 pixels (pixel size = 1 by 1 cm). The images were not radiometrically corrected, even though this can enhance the success of VI-based analyses [76]. This is permissible because VI can be reliably estimated from non-calibrated high-resolution images so long as they are obtained under the same radiometric conditions [77].
Within each plot, the location of the biomass-sampled quarter (1 by 1 m) was determined with reference to the plot fence, which was converted to local coordinates for quadrat demarcation on the drone images ( Figure 1). To synchronize with the drone photographing, the plant material in one random quarter of the quadrat (formed by the two diagonal axes) was collected by clipping all the grasses to the ground level as closely as possible in 2018, and one of the remaining three quarters was harvested randomly in 2019. In total, 216 samples (27 plots by 8 quadrats) were collected each year. After all the non-plant materials were removed, the clipped fresh grasses were oven-dried until weights stabilized, and this weight was regarded as the ground truth biomass. Analysis of variance (ANOVA) and least significant difference were performed to test the differences in AGB between the treatments.

Bare Ground Mapping and RGBVI Calculation
Bare ground devoid of any biomass in a plot was mapped using the Support Vector Machine classifier in ArcMap 10.7, together with grassy meadow. This classifier is a spectral classification method involving supervised learning. It is based on an iterative process of calculating the n-dimensional classification to find the optimum hyperplane boundary, which can classify data with unknown statistical distributions by using small training sets successfully [78,79]. As bare ground is visible on the drone images, all the training samples were randomly selected by visual inspection. Two indicators, overall accuracy and Kappa coefficient, were used to indicate classification accuracy based on 150 randomly selected evaluation points [80,81].
As the fence of each plot was clearly visible in the drone images, it served as the ground control for geo-referencing the UAV image of a plot to a dimension of 25 by 30 m in ArcMap ( Figure 1). The rectified image was clipped to 24.5 by 29.6 m to exclude the fence shadow, from which VI images were derived. These VIs were later used to estimate AGB from the 2018 and 2019 drone images separately to avoid any inconsistency in the illumination conditions between them. The Red Green Blue Vegetation Index (RGBVI) was used to estimate AGB in this study (Equation (1)); it has been successfully used for vegetation biomass estimation [82][83][84]. After the RGBVI calculation, pixels corresponding to the harvested quarter (0.25 m 2 ) in each quadrat were extracted (about 250 in total) from the geo-referenced VI images via their coordinates determined from the distance to the plot border.

AGB Modelling
Three types of model were built in this study ( Figure 3); separate regression models, both linear and second-order polynomial were built for the 2018 and 2019 images, and the best one was used to evaluate model performance. The initial model was constructed from the RGBVI values extracted from 216 pixels each year, during which the RGBVI value of the pixels falling into the same sampling quarter was averaged prior to its regression against the in situ measured AGB.
vegetation biomass estimation [82][83][84]. After the RGBVI calculation, pixels corresponding to the harvested quarter (0.25 m 2 ) in each quadrat were extracted (about 250 in total) from the geo-referenced VI images via their coordinates determined from the distance to the plot border.

AGB Modelling
Three types of model were built in this study ( Figure 3); separate regression models, both linear and second-order polynomial were built for the 2018 and 2019 images, and the best one was used to evaluate model performance. The initial model was constructed from the RGBVI values extracted from 216 pixels each year, during which the RGBVI value of the pixels falling into the same sampling quarter was averaged prior to its regression against the in situ measured AGB. The initial model was modified by two ratios. The first is the mowing ratio (F) based on mowing intensity that was treated as a proxy for vegetation biomass. The mowing effect on the biomass removal was proportional to the removed biomass corresponding to the reference plots. It is calculated as: where n is the number of sample plots having the same mowing intensity; B is the averaged in situ measured AGB of the sampled quarter in the reference plots (PnGn), and b is the measured AGB of the sampled quarter in the plots of the same mowing intensity. F was used to modify the RGBVI value: The initial model was modified by two ratios. The first is the mowing ratio (F) based on mowing intensity that was treated as a proxy for vegetation biomass. The mowing effect on the biomass removal was proportional to the removed biomass corresponding to the reference plots. It is calculated as: where n is the number of sample plots having the same mowing intensity; B is the averaged in situ measured AGB of the sampled quarter in the reference plots (PnGn), and b is the measured AGB of the sampled quarter in the plots of the same mowing intensity. F was used to modify the RGBVI value: where t denotes one of the three mowing intensities; (1 − F t ) k is the biomass density factor, which can reflect the vertical properties of the grassland; k is a constant determined via non-linear regression (k = 4 produced the best result for the study area but could be slightly different for other areas). The second ratio is the bare ground metrics. This modification was implemented by extracting the corresponding VI of bare ground pixels in a newly created mask layer and assigning a value of zero to them in the output metrics. These bare ground modification metrics were then used as the ground truth dataset, which can provide more modelling samples in model construction and validation. This metric was used as the supplementary training dataset based on the MF model to further improve the AGB estimation accuracy (Figure 3).
All models were validated using the "leave one out cross validation" (LOOCV) method. LOOCV can minimize the overfitting problem and provide an accurate assessment of model prediction strength, using each sample for validation and the rest for training [85,86]. The coefficient of determination (R 2 ), root mean square error (RMSE) and mean error were employed to evaluate the over-or under-estimated biomass. The last was used to assess the model's prediction errors to represent overestimations (positive number) and underestimations (negative number). The closer the mean error to zero, the more accurate the estimate.
where y i and Y i are the measured and modelled AGB, respectively, and n is the total number of samples. However, R 2 has an inherent problem in explaining the additional input variables relating to the model outputs. The changes in the training dataset may cause the overfitting problem in the model. So, adjusted R 2 is introduced to check whether the bare ground modification metrics improved model accuracy.

Measured AGB and Mowing Modification Ratio (F)
The field measured AGB is statistically summarized in Table 2. Although the year of 2018 had more AGB than the year of 2019 (178.20 vs. 145.64 g m −2 on average), the control plots and intensively mowed plots had the maximum and minimum AGB, respectively, in both years. In the mowed plots, AGB dropped along with the mowing intensity, but pika disturbance had no such effect as the medium disturbance group had the lowest AGB. Moreover, the AGB of pika treated plots had a narrow range from 166.47 to 188.46 g m −2 than the mowed plots from 140.52 to 200.24 g m −2 in 2018. Furthermore, the latter varied significantly between disturbance severities while the former did not. This finding still held true in 2019. Figure 4 displayed the mowing intensity modification ratio (F) based on the field collected data, which had a value less than 0.4 in both 2018 and 2019. F of the reference plots was close to zero, and the moderately mowed plots had a slightly higher value of 0.11 and 0.14 in both years, respectively. The severely mowed plots had a value of 0.36 on average. The average value at each severity was used for model modification. SD: standard deviations. LSD: least significant difference (α < 0.05). Lowercase letters (a-c) represent the significant difference in the two treatments separately.

Bare Ground Modification Metrics
In order to assess the effect of bare ground on AGB estimation, bare ground in ea plot has been mapped to an overall accuracy higher than 94.0% and a Kappa over 0.89. 2018, the RGBVI of the sampled quarter was mainly positive ranging from −0.03 to 0. (Figure 5a). Only few negative RGBVI values were found, which still corresponded to th field measured AGB values (above zero), while bare ground had smaller RGBVI value which were all negative numbers, without overlapping with the former. Although th sampled quarters had overlapping bare ground in 2019 (Figure 5b), their interquarti ranges were distinct.

Disturbance-Specific Models
The performance of the disturbance-specific models was evaluated based on the li ear model for each disturbance severity (none, medium and high) using RGBVI and th measured AGB in both 2018 and 2019 ( Figure 6). RGBVI had a close relationship with th measured AGB at all mowing severities, with an R 2 around 0.70 in both years, while th reference groups had a lower accuracy less than 0.70 and 0.68, respectively (Figure 6a, However, the models of the pika-disturbed groups were much less accurate with R 2 < 0.3 It is noteworthy that the lesser the mowing severity, the higher the measured AGB, b the rank of model reliability remains unchanged in both years. However, there were mo negative VIs in the medium and severe groups (joint disturbances by mowing and pik in 2019. Overall, mowing had significant impacts on the regression models, but not pik disturbance. This indicates the possibilities of improving the AGB estimation accura based on mowing intensity.

Bare Ground Modification Metrics
In order to assess the effect of bare ground on AGB estimation, bare ground in each plot has been mapped to an overall accuracy higher than 94.0% and a Kappa over 0.89. In 2018, the RGBVI of the sampled quarter was mainly positive ranging from −0.03 to 0.47 (Figure 5a). Only few negative RGBVI values were found, which still corresponded to the field measured AGB values (above zero), while bare ground had smaller RGBVI values, which were all negative numbers, without overlapping with the former. Although the sampled quarters had overlapping bare ground in 2019 (Figure 5b), their interquartile ranges were distinct.

Disturbance-Specific Models
The performance of the disturbance-specific models was evaluated based on the linear model for each disturbance severity (none, medium and high) using RGBVI and the measured AGB in both 2018 and 2019 ( Figure 6). RGBVI had a close relationship with the measured AGB at all mowing severities, with an R 2 around 0.70 in both years, while the reference groups had a lower accuracy less than 0.70 and 0.68, respectively (Figure 6a,c). However, the models of the pika-disturbed groups were much less accurate with R 2 < 0.30. It is noteworthy that the lesser the mowing severity, the higher the measured AGB, but the rank of model reliability remains unchanged in both years. However, there were more negative VIs in the medium and severe groups (joint disturbances by mowing and pika) in 2019. Overall, mowing had significant impacts on the regression models, but not pika disturbance. This indicates the possibilities of improving the AGB estimation accuracy based on mowing intensity.

AGB Models and Accuracy
Both linear and polynomial models were evaluated for their accuracy in estimating AGB (Figure 7). The best models were judged from the scatterplots in Figure 8. Generally, the polynomial models outperformed the linear ones except the initial models. Their accuracy in both 2018 and 2019 was rather low (R 2 < 0.45). However, model accuracy was highly improved if the model was modified with the mowing factor. Mowing modified models achieved a remarkable accuracy of R 2 around 0.70 (linear) and 0.82 (polynomial) in both years. The combination of both mowing and bare ground modifications further improved model accuracy by about 0.04 for 2018 and 0.08 for 2019 (Figure 7). In particular, the polynomial models achieved the highest R 2 of 0.86 for 2018, and 0.88 for 2019 ( Figure 8). It was worth noting that RMSE increased from 22.9 to 24.3 for 2018 and from 18.2 to 19.9 for 2019 after bare ground metrics were incorporated into the mowing-modified models. However, this metric helped to improve model reliability. models achieved a remarkable accuracy of R 2 around 0.70 (linear) and 0.82 (polynomial) in both years. The combination of both mowing and bare ground modifications further improved model accuracy by about 0.04 for 2018 and 0.08 for 2019 (Figure 7). In particular, the polynomial models achieved the highest R 2 of 0.86 for 2018, and 0.88 for 2019 ( Figure  8). It was worth noting that RMSE increased from 22.9 to 24.3 for 2018 and from 18.2 to 19.9 for 2019 after bare ground metrics were incorporated into the mowing-modified models. However, this metric helped to improve model reliability.   models achieved a remarkable accuracy of R 2 around 0.70 (linear) and 0.82 (polynomial) in both years. The combination of both mowing and bare ground modifications further improved model accuracy by about 0.04 for 2018 and 0.08 for 2019 (Figure 7). In particular, the polynomial models achieved the highest R 2 of 0.86 for 2018, and 0.88 for 2019 ( Figure  8). It was worth noting that RMSE increased from 22.9 to 24.3 for 2018 and from 18.2 to 19.9 for 2019 after bare ground metrics were incorporated into the mowing-modified models. However, this metric helped to improve model reliability.

Effects of Disturbance Severity and Bare Ground on AGB Estimation
The influence of the severity of sole and joint disturbances on AGB estimation was assessed in terms of mean error (Figure 9). The mowing-modified models underestimated AGB for the non-mowed groups (reference and pika) with an error around −10% in 2018, and a further −15% for the reference group of 2019. The mowed groups (mowing and joint) were all overestimated with an error about 10% in both years. On average, AGB was overestimated by these models with an error around 4.0% in both years. However, if modified by both mowing and bare ground, absolute mean error of estimation was reduced to about 1.0% for the reference groups and less than 3.0% for the other groups in both years. The mean error of all groups was found to be around 2.0% and 1.0% for 2018 and 2019, respectively. Although all the groups had little improvement after both mowing and bare ground modifications in both years, the average mean error of estimation was reduced to close to zero.
In order to assess the performance of bare ground metrics on improving AGB estimation, the two modified models were compared in terms of the predicted AGB in bare ground areas (Figure 10). The mowing-modified model overestimated bare ground AGB by 32-50 g m −2 for 2018 and 12-30 g m −2 for 2019. This estimation error was effectively reduced to a range from −5 to 18 g m −2 for 2018 and from −8 to 13 g m −2 for 2019 after bare ground metrics were incorporated into the mowing-modified models, although underestimation occurred in some treatments. On average, the consideration of bare ground modification metrics reduced the overestimation from 38 to 4 g m −2 for 2018 and from 19 to 2 g m −2 for 2019. Therefore, both modifications can help to reduce estimation errors caused by the soil background influence.

Effects of Disturbance Severity and Bare Ground on AGB Estimation
The influence of the severity of sole and joint disturbances on AGB estimation was assessed in terms of mean error (Figure 9). The mowing-modified models underestimated AGB for the non-mowed groups (reference and pika) with an error around −10% in 2018, and a further −15% for the reference group of 2019. The mowed groups (mowing and joint) were all overestimated with an error about 10% in both years. On average, AGB was overestimated by these models with an error around 4.0% in both years. However, if modified by both mowing and bare ground, absolute mean error of estimation was reduced to about 1.0% for the reference groups and less than 3.0% for the other groups in both years. The mean error of all groups was found to be around 2.0% and 1.0% for 2018 and 2019, respectively. Although all the groups had little improvement after both mowing and bare ground modifications in both years, the average mean error of estimation was reduced to close to zero. In order to assess the performance of bare ground metrics on improving AGB estimation, the two modified models were compared in terms of the predicted AGB in bare ground areas ( Figure 10). The mowing-modified model overestimated bare ground AGB by 32-50 g m −2 for 2018 and 12-30 g m −2 for 2019. This estimation error was effectively reduced to a range from −5 to 18 g m −2 for 2018 and from −8 to 13 g m −2 for 2019 after bare ground metrics were incorporated into the mowing-modified models, although underestimation occurred in some treatments. On average, the consideration of bare ground modification metrics reduced the overestimation from 38 to 4 g m −2 for 2018 and from 19 to 2 g m −2 for 2019. Therefore, both modifications can help to reduce estimation errors caused by the soil background influence.

Comparison of Modifications
Of the two modifications, mowing boosted the estimation accuracy to a level similar to that achieved using canopy volume [87][88][89]. However, the estimation accuracy was improved only slightly if the model was further modified by the bare ground metrics, alt-

Comparison of Modifications
Of the two modifications, mowing boosted the estimation accuracy to a level similar to that achieved using canopy volume [87][88][89]. However, the estimation accuracy was improved only slightly if the model was further modified by the bare ground metrics, although both modifications jointly improved the estimation accuracy to the highest level possible for both years (Figure 8). Mowing modification was more effective at improving the estimation accuracy than bare ground modification because it was applied universally to the entire experiment plots. In comparison, bare ground was much less extensive spatially in all the plots. It can be predicted that this modification would improve the estimation accuracy more had there been more denudated patches inside the selected experiment plots. Furthermore, after these modifications, the best estimation model changed its nature from linear to polynomial. This change is explained by the pyramidal distribution of grassland biomass vertically [90].
In this study, grazing was simulated via mowing. Studies elsewhere have shown that livestock grazing generally influences grassland AGB differently from artificial grazing simulated via mowing, but the differences are rather small [91]. Grazing effectively influences the vertical structure of biomass by removing biomass equivalently and producing a uniform, smooth horizontal biomass distribution [92]. In comparison, livestock preferentially graze palatable species [93], resulting in an uneven height of plants, while mowing always leads to a uniform height of plants, at least immediately following the mowing. This equivalent height of mowed plants is conducive to achieving a higher AGB estimation accuracy. For instance, the mowed plots had an R 2 over 0.70, slightly higher than that of the reference plots ( Figure 6), because the latter had an uneven plant height due to natural growth. Grass height affects the accuracy of biomass estimation [94]. The lack of consideration of vegetation height caused high AGB to be underestimated but low AGB to be overestimated by the initial models ( Figure 8). Such inaccurate estimations are attributed to the fact that the RGBVI derived from the drone images can only show the aerial view of the grass from above [95].

Implications
The achievement of a high estimation accuracy requires substantial training samples representative of all the variable vegetation characteristics [96]. An inspection reveals huge differences in RGBVI between the sampled and bare ground areas. Thus, the use of the AGB samples collected from the experimental site may not be sufficiently representative to achieve an accuracy higher than 0.88 because of the under-representation of low AGB. The models constructed from these samples typically overestimate AGB, especially in severely disturbed plots ( Figure 10). Therefore, there is a need to build bare ground metrics for expanding the training samples to reduce the estimation inaccuracy of bare ground AGB.
This study has demonstrated the feasibility of using local-scale drone images in estimation grassland AGB based on RGBVI at a rather high accuracy. Naturally, the accuracy will be lower (e.g., less than 0.70) if the estimation is implemented on a broad scale due to much wider spatial variability in vegetation conditions [97,98]. The proposed modification method can be easily replicated elsewhere. In the absence of mowing intensity data, the estimation accuracy can still be improved by substituting the mowing intensity by the stocking rate that can be regarded as a proxy for grazing intensity. This substitution is permissible, as grassland AGB bears a strong relationship with stocking rate and grazing intensity (R 2 > 0.72) [55,56]. Thus, the proposed modification methods have the potential to improve the AGB estimation accuracy even for livestock grazed grassland. If applied at a broad scale, the estimation of grassland AGB should be achieved from satellite images with near infrared and shortwave infrared bands. Such images enable other VIs to be derived. It will be interesting to investigate whether the same level of improvement in estimation accuracy can be achieved at a broad scale by both modifications.

Conclusions
This study proposes an efficient method to improve the accuracy of estimating grassland AGB from high-resolution drone images. The commonly adopted RGBVI estimation models were further modified by grazing (simulated via mowing) intensity and the quantity of bare ground mapped from the UAV images using support vector machine. Without such modifications, the estimation was achieved at a rather low accuracy (R 2 around 0.44) due to overestimation of low AGB values but underestimation of high AGB values, both by the same amount of 10%. However, the estimation accuracy was significantly improved to 0.81 after the model was modified by the mowing intensity that was considered as the proxy for the quantity of biomass on the ground. Nevertheless, further modification by the amount of bare ground improved the estimation only marginally (by 0.04 up to 0.08) to 0.88, an accuracy comparable to that achieved at a broad scale. The proposed mowing modification is more effective than bare ground modification in improving estimation accuracy, as it was applied universally to the entire experimental plots and is able to take into account the vertical structure of grasses, a key indicator of biomass in this instance. In comparison, bare grounds were of limited spatial presence in the study area. Since this modification markedly reduced the overestimation of low AGB values, it can be anticipated that the estimation accuracy would have been improved much more if there had been more bare ground enclosed in the experimental plots. The proposed methods can be used to assess the grassland carrying capacity quickly, thereby supporting more systematic and precise determination of the impacts of stocking rates and grazing practices.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

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