Quantifying Grazing Intensity Using Remote Sensing in Alpine Meadows on Qinghai-Tibetan Plateau

Remote sensing data have been widely used in the study of large-scale vegetation activities, which have important significance in estimating grassland yields, determining grassland carrying capacity, and strengthening the scientific management of grasslands. Remote sensing data are also used for estimating grazing intensity. Unfortunately, the spatial distribution of grazing-induced degradation remains undocumented by field observation, and most previous studies on grazing intensity have been qualitative. In our study, we tried to quantify grazing intensity using remote sensing techniques. To achieve this goal, we conducted field experiments at Gansu Province, China, which included a meadow steppe and a semi-arid region. The correlation between a vegetation index and grazing intensity was simulated, and the results demonstrated that there was a significant negative correlation between NDVI and relative grazing intensity (p < 0.05). The relative grazing intensity increased with a decrease in NDVI, and when the relative grazing intensity reached a certain level, the response of NDVI to relative grazing intensity was no longer sensitive. This study shows that the NDVI model can illustrate the feasibility of using a vegetation index to monitor the grazing intensity of livestock in free-grazing mode. Notably, it is feasible to use the remote sensing vegetation index to obtain the thresholds of livestock grazing intensity.


Introduction
Livestock grazing is a dominant form of land use globally, as nearly one-quarter of the earth's land surface is used for grazing (22%) [1].The interaction between livestock and grassland is the key ecological process of the grazing ecosystem, and the grass-livestock relationship is often a subject of research.Most studies have shown that moderate grazing can effectively promote grassland productivity and improve biodiversity maintenance, corresponding with the intermediate disturbance hypothesis [2,3].Thus, plant diversity should be highest at intermediate levels of grazing intensity.However, overgrazing can lead to a large area of grassland degradation, as well as the deterioration of grassland primary productivity and soil's physical and chemical properties [4].In recent years, grasslands across the globe have experienced increasing pressure from increasing human populations, reduced areas with increasing livestock numbers, and declining terms of trade for livestock production [5].Unsustainable livestock farming is considered the greatest driver of land degradation, which is defined as a decrease in productivity [6].
The Qinghai-Tibetan Plateau (QTP, average elevation 4000 m) is the largest plateau in the world, but it is a relatively fragile system [7][8][9].Alpine meadow is the dominant vegetation type (more than 40%) and provides fundamental sources of livelihood for local residents and ecosystem services for the human population living downstream [10].The Tibetan Plateau ecosystem is particularly sensitive to global climate change [11], and related research on the plateau ecosystem plays a highly important role in research on global change [12].Degradation of the alpine grassland has been widely recorded for the QTP over recent decades [13][14][15].In some areas, the magnitude of degradation has even reached an alarming level [16].The degradation is not only damaging local herdsmen's livelihoods, but also threatening the unique flora and fauna of the plateau [17].Due to long-term over-grazing and other factors, the alpine meadow ecosystems clearly appear to have degraded [18].In conjunction with this recent degradation, grassland management and optimization, especially for the study of a grazing intensity model, is of great significance.
Grazing intensity refers to the amount of quantitative animal demand for forage placed upon the standing crop forage or forage mass, and to the resulting level of defoliation caused during grazing [19].Research on grazing intensity mainly focuses on the effects of different grazing intensities on grassland productivity, vegetation community characteristics, and diversity [20][21][22][23][24]. Information on the grazing intensity of grassland is currently mostly derived from field investigations, with the locations of herbivores obtained by the satellite-based global positioning system (GPS) [25].However, ground-based methods (field investigation) are not practical over extensive geographic areas.Remote sensing data have been widely used in the study of large-scale vegetation activities [26][27][28], which have important significance in realizing grassland yield estimation, determining grassland-carrying capacity, and strengthening the scientific management of grasslands [29].Remote sensing data are also used for grazing-intensity estimation.Unfortunately, the spatial distribution of grazing-induced degradation remains undocumented by field observation [30].Most previous studies on grazing intensity have been qualitative or controlled experiments [31][32][33].If we used a vegetation index to quantify grazing intensity, we could overcome the regional limits of grazing intensity quantification and monitor grazing behaviour in real time and adjust the range of livestock activities, to a certain extent, to ease the pressure of grassland degradation.In addition, we could conduct research on discipline fusion, using satellite remote sensing data analysis on a large scale and implementing vegetation surveys and experimental analysis on a small scale-thus integrating large-and small-scale research-as well as conduct research on remote sensing monitoring technology based on the grass-soil-livestock interaction.
For this reason, we first studied the spatial distribution characteristics of biomass in the study area under different grazing intensities and formulated the following general hypotheses: (1) Grass aboveground biomass stock and grass production are maximized under moderately grazed conditions.Grazing can promote an optimization process, increasing the productivity of defoliated species [34].In contrast, the total biomass production may increase or be maintained depending on the less-preferred species' response (i.e., higher abundance or not).This approach is supported by previous studies [35].In other words, under moderate grazing pressure, the aboveground biomass shows an increasing trend with an increase in grazing intensity.(2) With the strengthening of the grazing intensity, communities' height, aboveground biomass, Pielou evenness, and index of dominance have increased and then decreased [36].Finally, (3) a continuous and intensive grazing history reduces total biomass production due to intense and sustained defoliation [35], and the aboveground biomass decreases with increasing grazing intensity, which implies that there is overgrazing.Based on these assumptions, we studied the relationship between the aboveground biomass (2016-2017 field measurement data) and different vegetation indices, and simulated the correlation between each vegetation index and grazing intensity, thus illustrating the feasibility of using a vegetation index to monitor the grazing intensity of livestock in free-grazing mode.The business model of the Qinghai-Tibet Plateau is dominated by herd husbandry; therefore, it is of practical value to carry out a grazing experiment at the scale of herders to promote the sustainable use of alpine grassland ecosystems in the future.

Study Area
The study area was located in Maqu County, Gansu Province, China, on the eastern edge of the Qinghai-Tibet plateau (101 • 52 E, 33 • 40 N, 3650 m a.s.l.(Above Sea Level)) (Figure 1).The annual average temperature there is 1.2 • C, and the average annual rainfall is approximately 620 mm, which is mainly concentrated in May-September.The annual average sunshine duration is approximately 2580 h.The average annual number of frost days is over 270, and there is no absolute frost-free period, only a warm-cold season.The total area of grassland in the county is 858,700 hm 2 , accounting for 89.3% of the total land area, and the available grassland area is 830,300 hm 2 [37].The soil is sub-alpine meadow soil.The grassland type is alpine meadow, and vegetation is dominated by Kobresia capillifolia, Poa chalarantha, and Roegneria nutans, among other species.Forage generally turns green in late April and begins to yellow in late August.Historically, yaks (Bos gunense), as the dominant livestock, have played a crucial role in the study area.business model of the Qinghai-Tibet Plateau is dominated by herd husbandry; therefore, it is of practical value to carry out a grazing experiment at the scale of herders to promote the sustainable use of alpine grassland ecosystems in the future.

Study Area
The study area was located in Maqu County, Gansu Province, China, on the eastern edge of the Qinghai-Tibet plateau (101°52′E, 33°40′N, 3650 m a.s.l.( Above Sea Level)) (Figure 1).The annual average temperature there is 1.2 °C, and the average annual rainfall is approximately 620 mm, which is mainly concentrated in May-September.The annual average sunshine duration is approximately 2580 h.The average annual number of frost days is over 270, and there is no absolute frost-free period, only a warm-cold season.The total area of grassland in the county is 858,700 hm 2 , accounting for 89.3% of the total land area, and the available grassland area is 830,300 hm 2 [37].The soil is sub-alpine meadow soil.The grassland type is alpine meadow, and vegetation is dominated by Kobresia capillifolia, Poa chalarantha, and Roegneria nutans, among other species.Forage generally turns green in late April and begins to yellow in late August.Historically, yaks (Bos gunense), as the dominant livestock, have played a crucial role in the study area.

Ground-Measured Data
The experimental plots were mainly set up in the pastoral areas of the Aze Animal Husbandry Science and Technology Demonstration Zone.The grasslands belong to three herds of households and are similar in type (alpine meadow), with similar production and utilization patterns (yak grazing); moreover, the grasslands are at the same level of growth.That is, the theoretical edible pasture yields are equal in the grasslands and share the same growth cycle, with no other disturbances besides grazing.According to the number of herding yaks, the grassland area, the distance from the farmer's shed, the activity route of the yaks' selective feeding, and the intensity of herding activities during the frequency of stay, a transect sampling design along the grazing gradient was chosen to establish a grazing gradient with increasing distance from a poultry ring.The site surrounding each poultry ring was stratified into three concentric levels.The concentric levels represented the following: the first grazing gradient area (GI), near the entrance to the poultry ring; the second grazing gradient area (GII), located at a medium distance from the entrance to the poultry ring; and the third grazing gradient area (GIII) (Figure 2), located far from the entrance to the poultry ring.In addition, a long-term, non-grazing plot was chosen as a control (CK).Nine sampling plots measuring 0.5 m × 0.5 m were randomly selected in each grazing gradient area, and sampling

Ground-Measured Data
The experimental plots were mainly set up in the pastoral areas of the Aze Animal Husbandry Science and Technology Demonstration Zone.The grasslands belong to three herds of households and are similar in type (alpine meadow), with similar production and utilization patterns (yak grazing); moreover, the grasslands are at the same level of growth.That is, the theoretical edible pasture yields are equal in the grasslands and share the same growth cycle, with no other disturbances besides grazing.According to the number of herding yaks, the grassland area, the distance from the farmer's shed, the activity route of the yaks' selective feeding, and the intensity of herding activities during the frequency of stay, a transect sampling design along the grazing gradient was chosen to establish a grazing gradient with increasing distance from a poultry ring.The site surrounding each poultry ring was stratified into three concentric levels.The concentric levels represented the following: the first grazing gradient area (GI), near the entrance to the poultry ring; the second grazing gradient area (GII), located at a medium distance from the entrance to the poultry ring; and the third grazing gradient area (GIII) (Figure 2), located far from the entrance to the poultry ring.In addition, a long-term, non-grazing plot was chosen as a control (CK).Nine sampling plots measuring 0.5 m × 0.5 m were randomly selected in each grazing gradient area, and sampling points were spaced according to the spatial resolution of remote sensing data, enabling precise positioning in remote sensing images.Each sampling point was measured three times and mowed in late June, July, August, and September.During the vegetation survey, the names, coverings, and heights of the main species in the sample area were recorded.Then, the ground was mowed, and the fresh weight was measured and dried at 70 • C until the dry weight was reached.
Sustainability 2019, 11, x FOR PEER REVIEW 4 of 14 points were spaced according to the spatial resolution of remote sensing data, enabling precise positioning in remote sensing images.Each sampling point was measured three times and mowed in late June, July, August, and September.During the vegetation survey, the names, coverings, and heights of the main species in the sample area were recorded.Then, the ground was mowed, and the fresh weight was measured and dried at 70 °C until the dry weight was reached.

Vegetation Index Extraction and Calculation of Relative Grazing Intensity
The normalized vegetation index, which is usually the ratio of the difference between the reflectance in the red (R) and near-infrared (NIR) bands, was used to characterize the quantity distribution and quality of surface vegetation.The enhanced vegetation index (EVI) is an "optimized" vegetation index with improved sensitivity in high-biomass regions and improved vegetation monitoring through a de-coupling of the canopy background signal and a reduction in atmospheric effects.NDVI and EVI were calculated by the following formulas [38,39]: To maintain temporal consistency with the data measured on the ground, the 16-day synthesized MOD13Q1 product with a spatial resolution of 250 m from the NASA data centre in late June-September of 2016-2017 was selected to extract NDVI and EVI values.
According to the different grazing gradients, the relative grazing intensities of three herdsmen were calculated according to the following formula: Relative grazing intensity = number of yaks in each herding family/area of grazing gradient area (concentric area with different radius) The formula does not indicate the quantified grazing intensity, but rather, the intensity manifestation of grazing.In this study, assuming that each grazing gradient region is a circle with a different radius, the three gradient areas were three concentric circles centered on the ranch entrance.The area of the different gradient is the area of the circle, and it is assumed that the number of yaks in each gradient of the same herder is the same.Although this measure is the same as that in the traditional grazing intensity algorithm-namely, the number of heads of livestock divided by the area of grassland-based on the assumption that the number of yaks in each gradient is the same, only the area is changed.Therefore, the calculation result was set as the relative grazing intensity, which qualitatively indicates the degree of grazing.Selected pastures are grazed year-round, and livestock slaughter or death and other factors lead to changes in the number of livestock; thus, the two years of relative grazing intensity are different (Table 1).A comparison of the relative grazing

Vegetation Index Extraction and Calculation of Relative Grazing Intensity
The normalized vegetation index, which is usually the ratio of the difference between the reflectance in the red (R) and near-infrared (NIR) bands, was used to characterize the quantity distribution and quality of surface vegetation.The enhanced vegetation index (EVI) is an "optimized" vegetation index with improved sensitivity in high-biomass regions and improved vegetation monitoring through a de-coupling of the canopy background signal and a reduction in atmospheric effects.NDVI and EVI were calculated by the following formulas [38,39]: To maintain temporal consistency with the data measured on the ground, the 16-day synthesized MOD13Q1 product with a spatial resolution of 250 m from the NASA data centre in late June-September of 2016-2017 was selected to extract NDVI and EVI values.
According to the different grazing gradients, the relative grazing intensities of three herdsmen were calculated according to the following formula: Relative grazing intensity = number of yaks in each herding family/area of grazing gradient area (concentric area with different radius) The formula does not indicate the quantified grazing intensity, but rather, the intensity manifestation of grazing.In this study, assuming that each grazing gradient region is a circle with a different radius, the three gradient areas were three concentric circles centered on the ranch entrance.The area of the different gradient is the area of the circle, and it is assumed that the number of yaks in each gradient of the same herder is the same.Although this measure is the same as that in the traditional grazing intensity algorithm-namely, the number of heads of livestock divided by the area of grassland-based on the assumption that the number of yaks in each gradient is the same, only the area is changed.Therefore, the calculation result was set as the relative grazing intensity, which qualitatively indicates the degree of grazing.Selected pastures are grazed year-round, and livestock slaughter or death and other factors lead to changes in the number of livestock; thus, the two years of relative grazing intensity are different (Table 1).A comparison of the relative grazing intensities between the three herders in 2016 and 2017 showed that the overall grazing intensity in 2016 was stronger than that in 2017.In 2017, the number of yaks was reduced because of grazing-livestock slaughter, resulting in a decrease in grazing intensity compared with that in 2016.

Data Analysis
NDVI and EVI values were extracted from the downloaded images using MODIS (moderate-resolution imaging spectroradiometer) data reprojection software (MRT).The NDVI and EVI values corresponding to each square point were extracted by Arc-Map software, and the model was established with the measured biomass.In the field sampling process, due to the GPS fixed-point error and the errors in the weighing and drying processes (spillover or moisture absorption), outliers should be eliminated before establishing regression models.Excel and SPSS (Statistical Product and Service Solutions) software were used for the statistical analysis of relevance.
The determination coefficient (R 2 ) and mean absolute error (MAE) were used to evaluate the accuracy of models and the interrelationships of the models.R 2 and MAE were calculated by the following formulas: where f i is the estimated value, y i is the ground measured value, n is the sample number, and cov and std represent the covariance and the standard deviation, respectively.

Correlation Analysis of Measured Aboveground Biomass and Relative Grazing Intensity
The aboveground biomass decreased with the increase of grazing intensity in 2016-2017 (Figure 3).With the change of the growth rhythm of forage, the aboveground biomass showed a convexity parabola change in June-September, which showed a tendency of increasing first and then decreasing gradually, and reached the highest level in August.In June and August of 2016, there were significant differences in biomass under grazing treatment, except for GIII and CK.Compared with CK, biomass under GI and GII significantly decreased in July 2016 and June 2017.In September, there were significant differences between the grazing treatments, except for GI and GII.In August 2017, biomass decreased significantly with the increase of grazing intensity.In June 2016, biomass was significantly lower than other months at the GI and CK levels (p < 0.05).At the GII and GIII levels, biomass in July and August was significantly higher than in June and September.In June 2017, from GI to GIII, biomass was significantly lower than the other months.A regression model between aboveground biomass and relative grazing intensity was established for 2016-2017 (Figure 4).A significant negative correlation (p < 0.01) was found between aboveground biomass and grazing intensity-that is, the aboveground biomass decreased with an increase in grazing intensity.

Correlation Analysis between Vegetation Index and Aboveground Biomass
Several fitting models (i.e., exponential, linear, polynomial, power) between vegetation index and aboveground biomass were established, and the results showed that the linear model had the best effect, so the linear regression models of two vegetation indices (NDVI and EVI) and aboveground biomass in 2016 and 2017 were established (Figure 5a-h).The results showed that the fitting between different vegetation indices and aboveground biomass exhibited the same tendency, as the biomass increased with the vegetation index.All the models satisfied the significant relationship between aboveground biomass and vegetation index (p < 0.01), and the R-squared values of the linear models of the NDVI exponent were all greater than those of the EVI exponential models for 2016 and 2017.A regression model between aboveground biomass and relative grazing intensity was established for 2016-2017 (Figure 4).A significant negative correlation (p < 0.01) was found between aboveground biomass and grazing intensity-that is, the aboveground biomass decreased with an increase in grazing intensity.A regression model between aboveground biomass and relative grazing intensity was established for 2016-2017 (Figure 4).A significant negative correlation (p < 0.01) was found between aboveground biomass and grazing intensity-that is, the aboveground biomass decreased with an increase in grazing intensity.

Correlation Analysis between Vegetation Index and Aboveground Biomass
Several fitting models (i.e., exponential, linear, polynomial, power) between vegetation index and aboveground biomass were established, and the results showed that the linear model had the best effect, so the linear regression models of two vegetation indices (NDVI and EVI) and aboveground biomass in 2016 and 2017 were established (Figure 5a-h).The results showed that the fitting between different vegetation indices and aboveground biomass exhibited the same tendency, as the biomass increased with the vegetation index.All the models satisfied the significant relationship between aboveground biomass and vegetation index (p < 0.01), and the R-squared values of the linear models of the NDVI exponent were all greater than those of the EVI exponential models for 2016 and 2017.

Correlation Analysis between Vegetation Index and Aboveground Biomass
Several fitting models (i.e., exponential, linear, polynomial, power) between vegetation index and aboveground biomass were established, and the results showed that the linear model had the best effect, so the linear regression models of two vegetation indices (NDVI and EVI) and aboveground biomass in 2016 and 2017 were established (Figure 5a-h).The results showed that the fitting between different vegetation indices and aboveground biomass exhibited the same tendency, as the biomass increased with the vegetation index.All the models satisfied the significant relationship between aboveground biomass and vegetation index (p < 0.01), and the R-squared values of the linear models of the NDVI exponent were all greater than those of the EVI exponential models for 2016 and 2017.

Correlation Analysis of Relative Grazing Intensity and Vegetation Index
A power regression model was established between the relative grazing intensities and the two vegetation indices on the corresponding gradients (Figure 6), and the results show that there is a significant negative correlation between them.With the decrease in NDVI, the relative grazing intensity increased, and when the relative grazing intensity reached a certain level (approximately 9 head•ha −1 ), the response of NDVI to the relative grazing intensity was no longer sensitive.The model fitting effect indices, such as the R-squared, mean absolute error, and significance values, were calculated (Table 2).Table 2 shows the simulation effect between the NDVI and relative grazing intensity is better than between the EVI and relative grazing intensity, and that the R-squared and significance values are greater than those of the EVI simulation.

Correlation Analysis of Relative Grazing Intensity and Vegetation Index
A power regression model was established between the relative grazing intensities and the two vegetation indices on the corresponding gradients (Figure 6), and the results show that there is a significant negative correlation between them.With the decrease in NDVI, the relative grazing intensity increased, and when the relative grazing intensity reached a certain level (approximately 9 head•ha −1 ), the response of NDVI to the relative grazing intensity was no longer sensitive.The model fitting effect indices, such as the R-squared, mean absolute error, and significance values, were calculated (Table 2).Table 2 shows the simulation effect between the NDVI and relative grazing intensity is better than between the EVI and relative grazing intensity, and that the R-squared and significance values are greater than those of the EVI simulation.

Aboveground Biomass under Different Grazing Intensities
A significant negative correlation was observed between total aboveground biomass and grazing intensity (p < 0.01)-namely, the aboveground biomass decreased with an increase in grazing intensity.This finding agrees with previous studies showing that with an increase in grazing intensity, the height, cover, and aboveground biomass exhibit a significant decreasing trend, such as in the Kobresia pima alpine meadow community of the eastern part of the Qinghai-Tibetan Plateau [40] and most grassland regions of northern China [41][42][43][44].Yan [45] found that grazing has negative effects on grassland biomass, and the grazing effects change with environmental conditions and grazing intensity.In addition, grazing increased the non-grass component of the vegetation, and

Aboveground Biomass under Different Grazing Intensities
A significant negative correlation was observed between total aboveground biomass and grazing intensity (p < 0.01)-namely, the aboveground biomass decreased with an increase in grazing intensity.This finding agrees with previous studies showing that with an increase in grazing intensity, the height, cover, and aboveground biomass exhibit a significant decreasing trend, such as in the Kobresia pima alpine meadow community of the eastern part of the Qinghai-Tibetan Plateau [40] and most grassland regions of northern China [41][42][43][44].Yan [45] found that grazing has negative effects on grassland biomass, and the grazing effects change with environmental conditions and grazing intensity.In addition, grazing increased the non-grass component of the vegetation, and species richness and diversity increased in the temperate grassland on the Mongolian Plateau.These results provide evidence backing the third hypothesis, as continuous and intensive grazing reduced total biomass production, which suggests a tendency towards overgrazing throughout the present study area.

Correlation between Different Vegetation Indices and Aboveground Biomass
Different vegetation indices were used to simulate aboveground biomass, and the results showed that the EVI was better than the NDVI for simulating aboveground biomass from June to August in 2016.However, the MODIS-NDVI-based linear model could better simulate aboveground biomass in the Maqu alpine meadow, and the NDVI increased with the aboveground biomass, which was consistent with the trend observed for the NDVI and biomass models for grazing pasture [46,47].In alpine grasslands of the Tibetan Plateau [48,49] and the Three-River-Source Region, Du [50] demonstrated that the NDVI was more closely related to AGB (aboveground biomass) than was the EVI.Liang [51] also found the NDVI-based AGB model to be the best (46%) among all linear remote sensing models they tested in the pastoral area of the southern Qinghai Province of China.Zhang [52] found that when the biomass is between 370-720 g•m −2 , the simulation results of the linear model for yield estimation was very good.In this study, the maximum biomass was 402.76 g•m −2 , and the linear model was simulated well by comparison, which is consistent with the previous results.
The comparison of aboveground biomass and vegetation index values in different months showed that the biomass and the vegetation index in September decreased.In general, the forage began to revert to green in May, and the aboveground biomass gradually increased from the beginning of the growing season to the maximum in August, and then entered the non-growing season.The grass began to turn yellow, and the compensatory growth of biomass lagged behind that of grazing; thus, grazing affected its community composition and decreased the biomass of grasses [29].Therefore, the aboveground biomass decreased in September, but the decrease was not significant.The red band is the main band for the photosynthesis of green plants.In September, yellowing of the aboveground biomass led to weakened photosynthesis and increased reflection in the red band.Therefore, at a certain biomass, the vegetation index in September rapidly declined, and the vegetation index was obviously smaller than that in June, July, and August.When months 6-9 were integrated, two independent parts were observed (obviously separated in September and June-August), and the overall fitting effect was not good.Therefore, September was modelled separately, and the simulation effect was determined to be better.In addition, the correlation coefficient of the model was relatively low, and may have been affected by the soil and climate conditions in the study area.In addition, the spatial resolution of the selected MODIS product data was limited, which also affects the accuracy of the model.

Correlation between Grazing Intensity and Vegetation Index
The comparison of the regression models between different vegetation indices and relative grazing intensity demonstrated that the simulation effect of the NDVI was better than that of the EVI (Table 2), and the results showed that there was a significant negative correlation between vegetation index and relative grazing intensity (p < 0.05).This result showed that there is a negative correlation between animal density and NDVI in the study area.However, some studies have shown that there is a positive correlation between animal density and NDVI.Ma [53] found that grasshopper density had a significant positive correlation with NDVI under different grazing intensities (p < 0.05).Guo [54] investigated the relationship between pikas' burrow density and the NDVI and found a significant positive correlation between the two.Pikas' burrow density is low when the NDVI is less than 0.2, indicating that Plateau pikas avoid dwelling in degraded grassland.Behnke [55] found that if larger water points were associated with more livestock, it is reasonable to assume that larger water points are also associated with increased grazing pressure and higher levels of forage depletion.High NDVI values near water points can examine this relationship and indirectly confirm the positive correlation between NDVI and animal density.However, Mongolian gazelles were shown to prefer an intermediate range of vegetation productivity, presumably facing quality/quantity trade-offs, where areas with a low NDVI are limited by low ingestion rates and areas with a high NDVI are limited by the low digestibility of mature forage [56].High levels of herbage mass with moderate grazing pressure enable animals to select better-quality diets, and when Mongolian gazelles prefer short or intermediate grasslands, according to an energy-maximizing strategy, higher NDVI values might not indicate better habitat quality for gazelles [57].Several factors may explain why gazelles did not use areas where the NDVI was higher.For instance, a higher NDVI may not directly reflect a better site for gazelles.For grazing ruminants, the daily energy intake rate will be the highest on intermediate biomass swards, where the grazer can optimize their dry-matter intake rate and energy digestibility [58,59].The above-mentioned studies all indicate that herbivores tend to choose a higher-productivity environment for feed intake over a wide range of areas.However, the results of our study indicated just the opposite.The NDVI is negatively correlated with the density of livestock, indicating that under limited resources, biomass will be reduced in areas where livestock concentrate their feeding, which implies overgrazing in the study area.
A threshold value of the relative grazing intensity can be found by analyzing the relationship between relative grazing intensity and NDVI.Under this value, while the response of NDVI to relative grazing intensity is no longer sensitive and can be used as a warning value for measuring the intensity of grazing degree, it also can be an important reference value for the reasonable grazing intensity.Of course, this value will vary depending on the study area or other circumstances.We know that the direct harvest estimation of grassland aboveground biomass requires considerable time and labor, as well as a sufficient number of samples to ensure that the samples meet statistical requirements (a limited number of samples leads to limited estimation accuracy) [60].In addition, the direct harvesting method will destroy vegetation in the survey site.Such destruction will affect the vegetation restoration in certain ecologically fragile areas, and thus, restoration will not be carried out over a large area; moreover, coupled to the activity of livestock is the process of spatial dynamic change.If the correlation between grazing intensity and aboveground biomass is used to determine the threshold of grazing intensity, heavy workloads, as well as spatial constraints, will have a strong effect on the research results.However, compared with traditional methods based on quadrat biomass monitoring, remote sensing estimation is not only more time-efficient but also has incomparable advantages over traditional methods, such as the monitoring of change in unmanned alpine areas [61].Therefore, it is feasible to use the remote sensing vegetation index to obtain the thresholds of livestock grazing intensity.However, the current study is only a preliminary exploration and there is still a lack of practical application, which will be supplemented in our following study.
Livestock feed species that are not resistant to grazing trampled plants, such as Potentilla anserina, Festuca arundinacea, Elymus nutans, and Poa pratensis, were shown to be greatly affected by grazing disturbance, and could not fully grow and develop; furthermore, their coverage and richness declined [62].That is, the growth rate of plants was lower than the feeding rate of domestic animals, which may be one of the reasons why the NDVI value decreased with the increase in grazing intensity.In addition, from the perspective of the whole community, grazing significantly reduced the aboveground biomass, mainly due to the reduction in leaf photosynthetic area and the reduction in organic matter accumulation efficiency [63].Cao [64] selected typical meadow grasslands in different degradation stages as their research object.Through a study of the mineralization ability and supply of nitrogen in the soil and the demand for nitrogen of forage grass, it was found that with an increase in grazing intensity, the supply and demand of soil-forage nutrition were imbalanced, and soil nutrients were biologically fixed and thus rendered ineffective by the excessive absorption of underground roots [43], which is also one of the major causes of biomass reduction and may lead to the decline of the vegetation index in the study area.An NDVI value of zero indicates the presence of rock or bare soil and that NIR and R are approximately equal; positive values indicate the presence of vegetation cover and a positive correlation with coverage.An unfavorable effect of overgrazing on grassland is that livestock excessively feed on the branches and leaves of plants, causing the leaf area of plants to continuously decrease and the photosynthetic capacity to decrease, thereby affecting the plant's growth and development and reproductive capacity [65] and reducing vegetation coverage [66][67][68]; this result may affect the vegetation index.The threshold NDVI in this study was within a positive range, indicating that the grazing intensity threshold did not result in exposure of the land and avoided the problem of vegetation damage at lower levels due to overestimation of the threshold.

Conclusions
A regression model between aboveground biomass and relative grazing intensity shows that the aboveground biomass decreased with an increase in grazing intensity.This result provided evidence backing the third hypothesis, as continuous and intensive grazing reduced total biomass production, which suggests a tendency towards overgrazing throughout the present study area.
This research applied two quantitative models-EVI and NDVI-to simulate the correlation with relative grazing intensity, and found that NDVI could better simulate the correlation and showed a significant negative correlation (p < 0.05), indicating that NDVI can be used to evaluate grazing intensity in the study area.The result in the model indicates that when the relative grazing intensity increased to a certain extent (approximately 9 head ha −1 ), NDVI was insensitive to it; thus, this grazing intensity value can be used as a warning value for measuring the intensity of grazing degree, and it also has significance for determining the appropriate grazing intensity under free-grazing modes.

Figure 1 .
Figure 1.The location of Aze station in Maqu county.

Figure 1 .
Figure 1.The location of Aze station in Maqu county.

Sustainability 2019 , 14 Figure 3 .
Figure 3. Aboveground biomass under different grazing intensity in 2016 (a) and 2017 (b).Note: Capital letters indicate differences between groups, and lowercase letters indicate intra-group differences.Note: Raw data have uploaded in the supplementary materials.

Figure 4 .
Figure 4.The power regression model of aboveground biomass and relative grazing intensity in 2016-2017.Note: Raw data have uploaded in the supplementary materials.

Figure 3 .
Figure 3. Aboveground biomass under different grazing intensity in 2016 (a) and 2017 (b).Note: Capital letters indicate differences between groups, and lowercase letters indicate intra-group differences.Note: Raw data have uploaded in the supplementary materials.

Sustainability 2019 , 14 Figure 3 .
Figure 3. Aboveground biomass under different grazing intensity in 2016 (a) and 2017 (b).Note: Capital letters indicate differences between groups, and lowercase letters indicate intra-group differences.Note: Raw data have uploaded in the supplementary materials.

Figure 4 .
Figure 4.The power regression model of aboveground biomass and relative grazing intensity in 2016-2017.Note: Raw data have uploaded in the supplementary materials.

Figure 4 .
Figure 4.The power regression model of aboveground biomass and relative grazing intensity in 2016-2017.Note: Raw data have uploaded in the supplementary materials.

Figure 5 .
Figure 5.Comparison of linear models of aboveground biomass and different vegetation indexes in 2016-2017.(a) Correlation between NDVI and aboveground biomass in June-August 2016; (b) Correlation between EVI and aboveground biomass in June-August 2016; (c) Correlation between

Figure 5 .
Figure 5.Comparison of linear models of aboveground biomass and different vegetation indexes in 2016-2017.(a) Correlation between NDVI and aboveground biomass in June-August 2016; (b) Correlation between EVI and aboveground biomass in June-August 2016; (c) Correlation between NDVI and aboveground biomass in September 2016; (d) Correlation between EVI and aboveground biomass in September 2016; (e) Correlation between NDVI and aboveground biomass in June-August 2017; (f) Correlation between EVI and aboveground biomass in June-August 2017; (g) Correlation between NDVI and aboveground biomass in September 2017; (h) Correlation between EVI and aboveground biomass in September 2017.Note: Raw data have uploaded in the supplementary materials.
NDVI and aboveground biomass in September 2016; (d) Correlation between EVI and aboveground biomass in September 2016; (e) Correlation between NDVI and aboveground biomass in June-August 2017; (f) Correlation between EVI and aboveground biomass in June-August 2017; (g) Correlation between NDVI and aboveground biomass in September 2017; (h) Correlation between EVI and aboveground biomass in September 2017.Note: Raw data have uploaded in the supplementary materials.

Figure 6 .Table 2 .
Figure 6.Comparison of power regression model of relative grazing intensity and different vegetation indexes in 2016-2017.Note: Raw data have uploaded in the supplementary materials.

Figure 6 .Table 2 .
Figure 6.Comparison of power regression model of relative grazing intensity and different vegetation indexes in 2016-2017.Note: Raw data have uploaded in the supplementary materials.Table 2. The results of the relative grazing intensity and two vegetation indices regression models.The R-squared value is the coefficient of determination, MAE is the mean absolute error.Vegetation Index Model R-Squared MAE Sig.NDVI y = 0.0544x −11.8 0.558 2.246 0.005 EVI y = 0.0833 −4.924 0.233 2.831 0.112

Table 1 .
The relative grazing intensities of the three surveyed herdsmen from 2016 to 2017.