The Quantitative Estimation of Grazing Intensity on the Zoige Plateau Based on the Space-Air-Ground Integrated Monitoring Technology

Grazing intensity (GI) is an important indicator for grazing situations in pastoral areas. However, it has been difficult to be observed directly in the field, due to the randomness and dynamics of the grazing behavior of livestock. Consequently, the lack of actual GI information has become a common issue in studies on quantitatively estimating GI. In this paper, a novel quantitative estimation method is proposed based on the Space-Air-Ground integrated monitoring technology. It systematically integrates GPS tracking technology, Unmanned Aerial Vehicle (UAV) observation technology, and satellite remote sensing technology. Taking Xiangdong Village on the Zoige Plateau as a study area, the trajectory data and UAV images were acquired by the GPS tracking experiments and UAV observation experiments, respectively. The GI at paddock scale (PGI) was then generated with the Kernel Density Estimation (KDE) algorithm and the above data. Taking the generated PGI as training data, an estimation model of GI at region scale (RGI) was constructed by using the time-series satellite remote sensing images and random forest regression algorithm. Finally, the time-series RGI data with a spatial resolution of 10 m in Xiangdong Village were produced by the above model. The accuracy assessment demonstrated that the generated time-series RGI data could reflect the spatial-temporal heterogeneity of actual GI, with a mean absolute error of 0.9301 and r2 of 0. 8573. The proposed method provides a new idea for generating the actual GI on the ground and the time-series RGI data. This study also highlights the feasibility and potential of using the Space-Air-Ground integrated monitoring technology to generate time-series RGI data with high spatial resolution. The generated time-series RGI data would provide data support for the formulation of policies and plans related to the sustainable development of animal husbandry.


Introduction
Animal husbandry is the pillar industry in pastoral areas. The benign and healthy development of animal husbandry is of considerable significance to the sustainable development of society and economy in these areas. Moderate grazing is beneficial not only to the productivity of grassland but also to the income of herdsmen [1,2]. In other words, it achieves a win-win situation of ecological and economic benefits. Overgrazing can yield temporary economic benefits. However, in the long term, it will destroy the grassland ecosystem, cause irreversible harm to animal husbandry, and finally lead to the double damages of ecological and economic benefits [2,3]. At present, an increasing number of pastoral areas are suffering from overgrazing due to the combined influences of multiple factors, such as tourism and the rising price of domestic animal products [4][5][6]. These phenomena have also been prominent on the Zoige Plateau on the eastern margin of the Qinghai-Tibetan Plateau [7]. They also lead to a series of eco-environmental issues, such as grassland desertification, weed or exotic invasion, and wetland degradation [8]. It is urgent to restore and repair the damaged eco-environment. More importantly, a series of scientific and reasonable control policies and development plans for animal husbandry must be formulated [9]. Understanding the grazing situations and quantitatively estimating the grazing intensity (GI) on the Zoige Plateau would provide data support for the formulation of corresponding policies and plans [10].
GI is the cumulative effect of grazing livestock on pastures during a particular period [11]. The livestock on the Zoige Plateau is usually fenced in each paddock. Therefore, the GI of a single paddock can be simply calculated as the ratio of livestock inventory to paddock area [2]; here, we call such a method the statistical method. However, the GI of a single paddock calculated by the statistical method is a constant in both space and time dimensions. In fact, due to the comprehensive effects of micro-topography, weed or exotic invasion, grassland degradation, and so on, the growth status of grass varies significantly within a paddock. At the same time, the discrepancies of phenology and growth environment for various grass species result in differences in their growth process. Coupled with the interferences of herdsmen in grazing behaviors, the actual GI in a paddock shows obvious heterogeneity in both time and space dimensions [12]. Obviously, the GI estimated with the statistical method ignores the spatial-temporal heterogeneity of the actual GI in a paddock.
Grazing reduces the vegetation coverage and aboveground biomass of grassland [13]. Many scholars have quantitatively estimated GI by using various remote sensing-based vegetation indices (such as NDVI, EVI) [8,[14][15][16][17][18]. The prominent merit of these methods is that a time-series GI dataset at a regional scale can be generated [2,8], because satellite remote sensing provides a unique capability for regular, repeated observations of the entire globe or any large areas. The actual GI information is necessary for these methods to train the estimation model and determine its parameters. However, the actual GI was difficult to observe directly in the field [2], because the grazing behavior of livestock in a paddock is random and dynamic [19]. As a result, the effects of various GIs on the spectral curve of grassland, vegetation index, biomass, or other indicators have been measured by limited fence experiments [20,21], and the experimental results have been taken as the ground truth for building the estimation model of GI. Considering the heterogeneity of actual GI in both time and space dimensions, it is urgent to find a new and more direct method to obtain the actual GI information in some local regions, such as a paddock, to serve the construction of GI estimation model.
A GPS collar is a device that automatically records the geographical location of the monitored object in real-time or at fixed intervals through navigation satellites (such as GPS or Beidou) or mobile base stations [22]. Such devices have been widely used in many studies; for example, wildlife monitoring and animal behaviors [23][24][25]. Cattle, sheep, and other livestock are social animals. The daily trajectories of these herds are easily tracked by a GPS collar attached to one or more livestock. These trajectory data can represent the grazing behaviors of monitored herds [26,27]. If the tracking points are concentrated in a region, it means that the grass in that region is more likely to be grazed by the monitored herds. In other words, the GI of that region is high. More importantly, the spatiotemporal dynamics of GI can also be identified by the trajectory data over a long time span. Therefore, the actual GI at paddock scale (PGI) might be quantitatively estimated based on these trajectory data recorded by GPS collars. Moreover, this generated PGI data can be taken as the ground truth for building the estimation model of GI at regional scale (RGI).
The livestock inventory of each paddock and its changes in the monitoring period are the critical coefficients of the PGI estimation model based on trajectory data. Field investigation is a common method to obtain this information. However, this approach faces many problems, such as inaccurate livestock inventory and its time-consuming nature [28]. Remote sensing provides an alternative approach to obtain surface information objectively and quickly. Unfortunately, it is difficult to directly identify individual livestock from satellite images, even if commercial satellite images with a resolution of less than one meter are adopted [29], because the size of individual livestock is too small for the pixel size of the prevailing satellite image [29]. The emergence and rapid development of Unmanned Aerial Vehicle (UAV) technology have made it possible to acquire higher spatial resolution images (less than one decimeter) quickly and flexibly [30,31]. At present, UAV images have been successfully used for identifying domestic and wild animals [28,30]. Apparently, it can also be used to identify livestock, which seems to offer a new opportunity for obtaining livestock inventory objectively and quickly.
On the whole, this study attempts to propose a novel method for estimating GI on the Zoige Plateau based on GPS tracking technology (Ground), UAV observation technology (Air), and satellite remote sensing technology (Space), i.e., Space-Air-Ground integrated monitoring technology. The specific objectives of this research are to (1) generate actual GI data at paddock scale (PGI) with spatial-temporal heterogeneity based on GPS tracking technology and UAV observation technology; (2) build an estimation model of GI at regional scale (RGI) using the generated PGI data and satellite remote sensing technology; and (3) produce time-series GI datasets.

Study Area
Xiangdong Village, Zoige County, located on the eastern margin of the Qinghai-Tibetan Plateau (Figure 1), was selected as a study area. Its average altitude is approximately 3500 amsl. The climate of this area is typical of alpine climates, with low temperatures and high humidity. The annual average temperature was 0.7-3.3 °C, and the annual accumulated precipitation was approximately 648.5 mm, mainly occurring from June to September [32]. Animal husbandry is the dominant industry, and grazing is the leading lifestyle of residents. In 1982, the implementation of the "household contraction of collective livestock" significantly increased the enthusiasm of herdsmen [33]. Subsequently, animal husbandry started to develop rapidly, and the livestock inventory increased obviously [34]. At the same time, with the increase in population, the per-capita pasture area began to decrease. To meet the needs of life, herdsmen started to expand the livestock inventory in their paddocks [34]. As a result, some problems, such as overgrazing and grassland degradation, began to worsen [7]. The native grasslands were divided into summer paddocks, winter paddocks, and all-year paddocks by local residents. In the summer paddocks, continuous grazing occurs from June to October. The grazing paddocks from November to May of the next year are defined as the winter paddocks [8].

Materials and Methods
A novel method for quantitatively estimating RGI on the Zoige Plateau based on the Space-Air-Ground integrated monitoring technology is proposed in this paper. The overall flowchart of the proposed method is shown in Figure 2. First, the trajectory data and UAV images were acquired, relying on Space-Air-Ground integrated monitoring technology. The PGI of the monitored paddocks were then produced by using the Kernel Density Estimation (KDE) algorithm and trajectory data. Taking the estimated PGI as training data, an estimation model of RGI was constructed by using the time-series satellite remote sensing images and random forest regression algorithm. Finally, the time-series PGI in Xiangdong Village was generated based on the above models. Details are given in the following sections.

Trajectory Data of Monitored Paddocks
In July-December 2018, a long-term tracking experiment using GPS collars on 10 randomly selected paddocks in Xiangdong Village was conducted. Because it is necessary to attach the customized GPS collars on the yaks, the selected paddocks are a compromise between the randomness of the spatial distribution and the support of the paddock's owner. The geographical location of these selected paddocks is shown by the green dot in Figure 1. For each monitored paddock, two adult male yaks were randomly selected to attach the same GPS collar. The collar has two kinds of positioning modes: one is GPS and Beidou dual-mode positioning, and the other is LBS (Location Based Services) positioning. When the GPS or Beidou signal is active, the dual-mode positioning is adopted, which is also the default positioning mode. Otherwise, the LBS positioning mode is used. Battery capacity is one of the limiting factors for the GPS collar [26]. Once the battery is exhausted, the GPS collar needs to be retrieved from the yak. It can then be reattached after the battery is fully charged. To reduce the time of battery charging and increase the time span of trajectory data, the time interval of the GPS collars was set to 1 h according to test reports of the manufacturer. The trajectory data recorded by the GPS collars were uploaded directly to the data server through the mobile communication network. The researcher can download these trajectory data from the data server at any time without retrieving the GPS collars.
From 20 July to 2 December, 2018, with the help of local herdsmen, the trajectory data of 10 monitored paddocks was collected. The preliminary analysis of these trajectory data ( Figure 3) showed that most of the trajectories of herds from August to November in eight paddocks (Paddocks #2, #3, #4, #5, #6, #7, #8, and #10) were entirely recorded. Unfortunately, for Paddocks #1 and #9, only a small amount of trajectory data was recorded because their collars were damaged in the early stages of monitoring. The trajectory data of Paddocks #1 and #9 were excluded in the following analysis. To verify the positioning accuracy of GPS collars, spatial overlay analysis was adopted. The trajectory data was overlaid with satellite images and UAV images to judge whether the tracking points fall within the monitored paddocks. There were apparent position offsets between the tracking points obtained by LBS mode and the boundary of monitored paddocks. Almost all of the tracking points obtained by the GPS and Beidou dual-mode positioning were located within the boundary of monitored paddocks. Therefore, only tracking points obtained by the dual-mode positioning and located within the boundary of monitored paddocks were considered in the following analysis.

UAV Images
Three aerial observation experiments (20 July, 28 September, and 2 December, 2018) were conducted during the monitoring period, which were synchronized with the attaching and retrieving of the GPS collars. All regions covered by the 10 monitored paddocks were observed. The first experiment was used to calculate the area and livestock inventory of each paddock ( Table 1). The main purpose of the second and third experiments was to find the changes in the livestock inventory of each paddock during the monitoring period. A DJI Phantom 4 Pro UAV was adopted for the aerial observation experiment. It was equipped with a visible-light camera with an 8.8 mm focal length and 20 megapixels. In the UAV experiments, the flight altitude was set at 400 m, and the corresponding ground spatial resolution was approximately 8 cm. The livestock and fence could be directly identified from these UAV images. The forward overlap and sideways overlap were set to 80% and 65%, respectively, to seamlessly mosaic the acquired UAV images. The entire UAV flight campaign was conducted from 9:00-14:00, with no wind. The flight time of the fully charged UAV was approximately 20 min, and the corresponding coverage area was approximately 1 km 2 . To ensure that all regions covered by each paddock were observed, the spatial extension of UAV flight was far larger than the boundary of each paddock. Therefore, for each paddock, at least two flight campaigns were conducted. During the UAV observation, a handheld Trimber GPS receiver was used to simultaneously record the geographic location of some common reference objects, such as the corners of fences and the intersections of the road. These objects have been used for the geometric registration of UAV images.
Several UAV observation experiments have been conducted by the research team on the Zoige Plateau, and the geometric offset of the post-processing UAV images was less than 0.5 m [35]. The same preprocessing methods were adopted in this study to ensure the geometric accuracy of UAV images. After a series of preprocessing steps, including geometric rectification and image mosaicking, UAV images of each monitored paddock were generated.

Satellite Images
MODIS image with a spatial resolution of 250 m is the most commonly used data source for estimating GI due to the high temporal resolution [8,36]. Field investigations have found that the width of most paddocks on the Zoige Plateau is approximately 100-300 m and that their length is approximately 3-8 km. When MODIS image was adopted, the generated GI data was difficult to characterize the spatial heterogeneity of actual GI in paddock scale. Therefore, GI with a finer spatial resolution was needed. Landsat-8 OLI images and Sentinel-2A/B images are ideal data sources for estimating GI at regional scale. Their spatial resolutions are between 10 m and 30 m, which can depict the spatial heterogeneity of GI in a paddock. The revisit cycle of Landsat-8 is 16 days, while that of Sentinel-2A/B is 5 days [37]. Generally, the frequency of cloud cover is high on the Zoige Plateau during the experimental period. Therefore, Sentinel-2A/B images were chosen. During the monitoring period, four Sentinel-2A/B Level 1C products [37] (top of atmosphere (TOA) reflectance data and cloud coverage ≤ 5%) were selected in this study ( Table 2). All Sentinel-2A/B images were downloaded from the United States Geological Survey (USGS) EarthExplorer platform. Bands 1, 9, and 10 have a relatively coarse spatial resolution (60 m) and are primarily designed for radiometric correction purposes (cirrus detection and estimation of water vapor content). Thus, only bands 2 to 8, 8A, 11, and 12 were used in GI estimation. The bands with a spatial resolution of 20 m (bands 5-7, 8A, 11, and 12) were resampled to 10 m.

Land Cover Data
Grazing mainly occurred on grassland. Other land cover classes needed to be excluded when estimating RGI. A 30-m-resolution land cover product (ChinaCover-SW) [38] generated by the research team was adopted to extract the spatial distribution of grassland. The accuracy of this product reached 87.14% [38]. It contained 38 land cover classes, but only the spatial distributions of meadow, steppe, tussock, and sparse grassland were used in this study. To eliminate the potential impacts of classification errors and spatial resolution inconsistency on the estimation results, the land cover product was overlaid with a Sentinel-2A image acquired on 23 August, 2018 (the appropriate period to distinguish between grassland and non-grassland) to identify omission errors and commission errors of grassland. For omission errors, new grassland regions were manually added, while for commission errors, the misclassified grassland regions were removed.

Livestock and Fence Identification Based on UAV Images
The vast majority of yaks on the Zoige Plateau are black ( Figure 4d) and appear as black clumps in the acquired UAV-based RGB images (Figure 4b). The grassland shows as light to dark green in RGB images depending on their coverage (Figure 4b, c). These spectral differences between livestock and grassland are obvious in the acquired UAV-based RGB images, which can be used to visually identify individual livestock. The livestock inventory of each paddock can be obtained by counting the number of livestock identified. Moreover, the texture and spectral features of adjacent paddocks are also noticeable (Figure 4c) due to the differences in grazing intensity, grazing time, and so on. Therefore, the fences can be easily identified and vectorized from these UAV-based RGB images [12], and the area of each paddock can be obtained by calculating the area of these identified fences.

Estimation of PGI Based on Kernel Density Estimation and Trajectory Data
The estimation of PGI based on the hourly trajectory data could be regarded as the issue of estimating overall characteristics based on samples, because only 24 tracking points were recorded by GPS collars each day. Kernel density estimation (KDE) is just a method for estimating the overall probability density function based on limited samples [39,40]. Therefore, KDE was adopted in this study to generate PGI. It posits that every location in a region has a measurable event intensity, and the event intensity of that location can be estimated by all event points within a certain distance [41]. Specifically, it takes each target point x as the center and calculates the density contribution of each tracking point within a specified distance (a circle with bandwidth h) by using a kernel function [40]. For each tracking point, the closer it is to the target point, the higher its density contribution. GI is a time-independent quantity. Therefore, the effects of monitoring duration and sampling frequency need to be eliminated. The KDE value at target point x is a dimensionless probability value [39] that needs to be multiplied by the livestock inventory in a paddock to be converted to GI. The basic principle of estimating PGI based on KDE is given in Formula 1: where GI(x) stands for the GI at the target point x; T represents the monitoring duration in days; f refers to the sampling frequency of the GPS collar within a day; N is the livestock inventory in a single paddock, and the unit is Sheep Unit (SU); h stands for bandwidth, which can be quantitatively described by the radius of the spatial distribution zone of monitored herds; n is the number of tracking points that fall within the bandwidth of the target point x, and a is the adjustment factor; K represents the kernel function, and its value decreases with the increase in the distance between the target point x and the tracking point xi. In this paper, the Gaussian kernel was chosen as the kernel function, and its calculation formula is as follows: The spatial resolution of the selected satellite images is 10 m. To conveniently extend the PGI to RGI without other additional operations, the spatial resolution of the generated PGI was also set to 10 m. N and h are two critical parameters for the proposed PGI estimation method. The livestock inventory was obtained by counting the number of livestock identified in each paddock. Generally, the type of livestock raised in each paddock might be different. To make the generated GI comparable across paddocks, the numbers of various livestock were uniformly converted into SU in this paper. The conversion coefficient used by Wang et al. [8] was adopted in this study. Among them, 1 sheep = 1 SU and 1 yak = 4 SU. The simplest way to obtain bandwidth h was to measure the largest circular envelope of each monitored herd in space. However, the spatial distribution of each herd acquired by a single UAV observation was not sufficient to represent their general conditions. Therefore, a linear regression relationship between the livestock inventory and bandwidth was established with the least square method, based on the statistical data of livestock inventory and bandwidth of herds that came from the interpretation results of UAV images. The bandwidth of each herd was calculated by using the constructed linear relationship. Since the bandwidth of some tracking points near the fence exceeded the boundary of the paddock, some pixels in adjacent areas outside the paddock might also be assigned a GI value when the KDE algorithm was adopted to estimate PGI, which might underestimate the GI for a single paddock. Therefore, the adjustment factor a was adopted to eliminate the systematic underestimation of PGI. The value of the adjustment factor was determined by the degree of deviation between the statistical GI of each paddock and their initial estimated PGI.
The GI is a time-independent quantity. Therefore, the GI at any period can be estimated. Because the PGI data would be taken as training data to build the RGI estimation model, the time span of the generated PGI should be consistent with the acquisition time of the selected time-series remote sensing images used to retrieve RGI, which means producing PGI data from August 23 to November 1 (PGI-PA). To verify the performances of the proposed RGI model in producing time-series RGI datasets, the time-series PGI datasets were also produced simultaneously, i.e., the PGI data from 23 August to 22 September (PGI-P1), from 22 September to 17 October (PGI-P2), and from 17 October to 1 November (PGI-P3).

Estimation of RGI Based on the Random Forest Regression Algorithm
Grazing by livestock changes the surface reflectance of grassland [42]. The relationship between GI and the changes of surface reflectance of remote sensing images during the monitoring period was constructed with a nonlinear function (f) using the random forest regression method: where ΔSR x represents the change values of TOA reflectance of band x in two Sentinel-2A/B images at the beginning and end of the monitoring period. A total of 10 bands (bands 2 to 8, 8A, 11, and 12) were adopted to construct the RGI estimation model in this study. Because the NDVI was sensitive to grazing [19,43], the change value of NDVI (ΔNDVI) in the same period was an optional dependent variable to construct the estimation model of RGI. Considering that vegetation phenology is also an important factor causing NDVI changes, except grazing, it is necessary to minimize the impacts of vegetation phenology on GI estimation when ΔNDVI is adopted in model. Usually, NDVI changes in ungrazing regions were mainly caused by vegetation phenology. Therefore, the average of ΔNDVI in ungrazing regions (MΔNDVIug) was calculated to represent the NDVI changes caused by the vegetation phenology in this study. Accordingly, the ΔNDVI used in the estimation model is the change values of NDVI of two corresponding Sentinel-2A/B images subtracted to the MΔNDVIug in the same period. The random forest regression algorithm (RF) was adopted to construct the estimation model of RGI. In RF, each tree is built using a deterministic algorithm by selecting a random set of variables and a random sample from the training dataset [44]. The final RF predictor is formed by taking the average over all trees. To assess the explanatory power of each variable, Gini importance was adopted and taken as an important indicator of feature relevance. A higher Gini importance value indicates that the variable is relatively more important [45]. This strategy increases the diversity between trees to avoid over-fitting, increases the robustness of the model, and has the ability to handle high data dimensionality and multicollinearity [45,46]. Currently, RF has been widely used in ecological parameter inversion, impervious surface estimation, and so on [46,47]. In this study, the number of variables extracted at each split (mtry) was set to 3, and the number of trees (ntree) was set to 500, based on the results of a large number of experiments and referring to the parameter values in similar studies [46,47].
The generated PGI was used as the actual GI to train the RGI estimation model. It was further divided into training samples and validation samples according to the following three criteria: (1) the homogeneity of spatial distribution, (2) the representativeness of the GI value, and (3) randomness. Specifically, the generated PGI data of each paddock was first divided into several sub-data at intervals of 1 SU/ha. The total number of pixels in each sub-data was then counted, and the size of training samples to be extracted from each sub-data was calculated according to the proportion of 70%. A pseudo-random number function was used to determine the location of training samples in each sub-data [47], and the remaining data were validation samples.

Accuracy Assessment of GI
The accuracies of both PGI and RGI were assessed in this study. For a single paddock, the GI can be simply calculated by the ratio of livestock inventory to paddock area, which we call the statistical GI (SGI). In this study, the SGI values of eight selected paddocks were calculated. Simultaneously, the average values of PGI and RGI of each selected paddock were counted. By comparing the PGI average with SGI, as well as the RGI average with SGI, the reliability of the generated PGI and RGI could be assessed, respectively. The independent validation samples were used to quantitatively evaluate the accuracy of RGI and the spatiotemporal scale expansion ability of the RGI model. The mean absolute error (MAE), root mean square error (RMSE), and coefficient of determination (r 2 ) were calculated [47] as follows: where Mi represents the RGI at pixel i; Ri stands for the corresponding PGI at pixel i; cov(M-R) 2 is the covariance of the RGI and PGI; and var(M) and var (R) are the variances of RGI and PGI, respectively.

The Trajectory of the Monitored Paddocks
Based on the GPS collars, hourly trajectory data of eight monitored herds over the whole monitoring period (21 July -1 December) were acquired completely. Figure 5 shows the hourly trajectories of two yaks recorded by GPS collars in Paddock #4 during three consecutive days. There exists a certain degree of similarity between the trajectories of two monitored yaks. They usually left their habitat and grazed freely on the paddock at approximately 7:00 a.m., and returned to the habitat at approximately 7:00 p.m. The daily routes of grazing were significantly different, even on two adjacent days ( Figure 5). The average grazing time of livestock was approximately 13 h. Therefore, only the tracking points between 7:00 a.m. and 7:00 p.m. could be used to produce PGI. In other words, only 13 tracking points were adopted for every day.

Livestock and Fence of Each Paddock
The livestock and fence of each monitored paddock were identified from the UAV images ( Figure 6). Simultaneously, the area of each monitored paddock and its livestock inventory for three UAV observation experiments were counted (Table 1). During the period of the first and second UAV observations, there was no significant change in the livestock inventory for each paddock, and only one yak was slaughtered for daily subsistence in some paddocks (Paddocks #2, #4, #5, #8, and #9). However, the livestock inventory of six monitored paddocks (Paddocks #1, #4, #5, #7, #8, and #10) changed significantly when the results of the third UAV observation were compared with those of the previous two UAV observations. The change in livestock inventory was related to the fact that adult yaks began to be sold during this period. The field investigation also confirmed that the date for livestock sales in 2018 was concentrated around 20 November. However, the monitoring period of this study was from 23 August to 1 November 2018. The changes in livestock inventory have no effects on the estimation results.  The linear regression relationship between the livestock inventory and bandwidth was constructed using the least-squares method, as shown in Figure 7. The spatial aggregation of the herd (bandwidth) has an obvious positive linear correlation with the size of the herd (livestock inventory). The spatial distribution pattern of livestock in a paddock showed an obvious aggregation effect, which was confirmed by the identified results and field investigation. That is, the grazing behavior of individual livestock was greatly affected by the grazing behavior of the whole herd. Therefore, it also indirectly confirmed the hypothesis of the GPS tracking experiment that the trajectory of individual livestock can represent the trajectory of the whole herd.

The PGI of the Monitored Paddocks
Based on the estimation method described in Section 3.6, the PGI-PA and time-series PGI data (PGI-P1, PGI-P2, and PGI-P3) of eight monitored paddocks with complete trajectory data were generated. The spatial distribution of the generated PGI data for four of them (Paddocks #4, #5, #8, and #10) are shown in Figure 8. From August 23 to September 22, the livestock grazed on the summer paddock (Figure 8a-1, 8b-1). The livestock migrated gradually from the summer paddock to the winter paddock between September 22 and October 17 (Figure 8a-2, 8b-2). By comparing the trajectory data and asking herdsmen of these paddocks, the seasonal migration from summer paddock to winter paddock was found to occur in early October. From October 17 to November 1, the livestock grazed on the winter paddock (Figure 8a-3, 8b-3). For Paddocks #8 and #10, the obvious migration did not occur during the monitoring period (Figure 8c, 8d), but there are some differences in the spatial pattern of PGI.

The Accuracy of PGI
To assess the accuracy of PGI, a comparison between the average value of PGI and corresponding SGI of the eight monitored paddocks was conducted (Figure 9). The average value of PGI was very close to the corresponding SGI for all monitoring periods, with all r 2 greater than 0.95. On the whole, the generated PGI data was reliable and could be taken as the actual GI for the training of the RGI estimation model. The advantage of the generated PGI data was that it could describe the spatial-temporal heterogeneity of the actual GI of a single paddock and enriched the information of traditional SGI.

The RGI of Xiangdong Village
The RGI of the entire monitoring period (RGI-PA) with a spatial resolution of 10 m in Xiangdong Village has been produced based on the PGI-PA data and the estimation model described in Section 3.7. Referring to Wang's grading method of GI on the Zoige County [8], the spatial distribution of graded RGI-PA in Xiangdong Village is shown in Figure 10d. Moreover, the constructed model has also been directly used to estimate the time-series RGI in Xiangdong Village, that is, the RGI data from 23 August, 2018 to 22 September, 2018 (RGI-P1), from 22 September, 2018 to 17 October, 2018 (RGI-P2), and from 17 October, 2018 to 1 November, 2018 (RGI-P3), which are shown in Figure 10a, 10b, and 10c, respectively. For the RGI-PA, light grazing (1-5 SU/ha) was dominant in Xiangdong Village, and its proportion reached 77.30%. However, heavy grazing (10-20 SU/ha) still existed in some areas, and its proportion was approximately 0.60%. For the RGI-P1 data, approximately 1.06% of GI consisted of heavy grazing, which exceeded its proportion in the RGI-PA. The proportion of heavy grazing decreased slightly in the RGI-P2 data because a large number of livestock began to migrate from summer paddocks to winter paddocks. The grazing regions of livestock increased, and the GI per unit area decreased. However, extremely heavy grazing (> 20 SU/ha) existed in some areas in this period, with a proportion of 0.18%. For the RGI-P3 data, the proportion of heavy grazing was increased obviously, which was mainly related to two reasons. On the one hand, the monitoring period was relatively short. On the other hand, to ensure that the whole herd has plentiful forage throughout the winter, herdsmen would artificially limit the grazing on a sub-paddock for a period. Frequent grazing in a small region during a short time would inevitably increase the GI in the grazing regions.

The Accuracy of RGI
Two strategies were adopted to assess the accuracy of the RGI in Xiangdong Village. One relied on the independent validation samples. The other was consistent with that adopted by the PGI. Figure 11 shows the fitting and predicted results of the training samples ( Figure 11a) and independent validation samples (Figure 11b). For the fitting performance of the estimation model of RGI, the MAE was 0.5298, with an r 2 of 0.9684, indicating that the models can fit the training samples well by the selected input variables. The independent validation samples did not take part in the model construction, which only served to reveal the predictability of the constructed model. According to the statistics of error indicators of the independent validation samples, the RMSE was 0.01546, with an MAE of 0.9301 and r 2 of 0.8573. However, there was a phenomenon that cannot be ignored: some low values were overestimated, while some high values were underestimated. For regions with a low GI value, the natural growth of grass might enhance the changes in surface reflectance during the monitoring period, which would lead to the overestimation of GI. For regions with a high GI value, the compensatory growth of grass would reduce the changes in surface reflectance of the monitoring period.   To assess the accuracy of RGI, a comparison between the RGI average and corresponding SGI of the eight monitored paddocks was also conducted ( Figure 13). In general, the RGI-PA averages of these paddocks were close to their SGI, with an r 2 of 0.8531, indicating that the RGI estimated by the proposed method could reflect the actual grazing of paddocks. The estimation performance of the model for producing the RGI during the whole monitoring period was better than that for producing time-series RGI.

The Role of Space-Air-Ground Integrated Monitoring for the Estimation of GI at Regional Scale
The specific contribution of ground monitoring is to collect the actual grazing trajectory of livestock, which is the major challenge for obtaining actual GI information, because the grazing behavior of livestock is random and dynamic [19]. Generally, a single or limited repeated observation does not provide a complete picture of the actual GI. Long-term continuous observation with the GPS tracking technology, such as the one adopted in this study, is an appropriate strategy. The grazing trajectories of livestock in eight monitored paddocks over a long monitoring period were recorded by GPS collar in this study, which further proved the availability and performance of the long-term observation. The advantage of this technology is that it can record the trajectory information of livestock automatically and regularly, with no human operation, which ensures the objectivity and authenticity of the data to some extent.
The role of UAV observation was to obtain the accurate livestock inventory and area of the monitored paddocks in this study, which are the critical parameters for the estimation model of PGI. Usually, a household survey is the primary method to obtain these parameters. However, this approach faces many problems, such as inaccurate livestock inventory and its time-consuming nature [28]. Considering the size of individual livestock and the width of the fence (Figure 4), they are difficult to be directly identified from satellite images [29]. The UAV observation simply provides an effective means to solve the above problems. Generally, a UAV image has a high spatial resolution (< 0.1 m), which could accurately identify yak, sheep, and other livestock in the paddock [48]. At the same time, the flexibility of the UAV allowed it to acquire images of any area at any time, and the frequency of UAV observation could be adjusted flexibly according to the requirements and weather conditions [49]. Therefore, it was very suitable for detecting the changes in livestock inventory.
The main contribution of the satellite monitoring to the estimation of GI was to realize the spatial-temporal scale expansion of GI [50]. MODIS images with a spatial resolution of 250 m have been the most commonly used data for estimating GI. Such images could reflect the differences in GI among paddocks, but there is difficulty in describing the spatial heterogeneity of GI within a paddock. The spatial resolution of GI obtained in this study reached 10 m. It could adequately express the spatial-temporal heterogeneity of GI within a paddock.
In this study, the Space-Air-Ground integrated monitoring was conducted in the fenced grazing areas. For free-grazing areas, this proposed method is also suitable. As the proposed method is related to the livestock inventory and grazing area, it is preferable that the monitoring regions are closed. In other words, the livestock inventory in the monitoring area is fixed, or the changes in livestock inventory can be monitored. For public pastures without fences, the performance of the proposed method needs to be further assessed because the livestock inventories and their changes are difficult to monitor. The standard sheep unit (SU) was taken as the unit of GI in this study. Therefore, this method is also suitable for the coexistence of various livestock in the same paddocks, which extends the application scope of the proposed method.

The Effects of Phenology on Estimation of RGI
The quantitative estimation of GI using remote sensing technology usually depends on the changes of TOA reflectance or NDVI of satellite images in different periods. In fact, both grazing and vegetation phenology can lead to changes in TOA reflectance or NDVI. Therefore, it is necessary to reduce the interference of phenology to the RGI estimation result. The GI estimation models considering the phenology and not considering the phenology were constructed in this study. The precision indicators (Table 3) show that the MAE of the model considering the phenology is smaller and r 2 is higher than that of the model not considering the phenology. The average GI in Xiangdong Village estimated by the model considering the phenology was 0.5650 less than that estimated by the model not considering the phenology. The effect of phenology is considered in model construction, which contributes to improving the estimation accuracy of the time-series GI datasets.

Implications and Future Work
GI is a significant indicator to reflect the degree of grassland utilization [2]. Quantifying GI will contribute to understanding the overall grazing situation in a region, in particular identifying the core regions where overgrazing occurs. More importantly, the decision-makers can formulate various regulation policies according to the situation of GI in different regions to serve the sustainable development of animal husbandry. The effects of grazing on grassland productivity, ecosystem services, vegetation community characteristics, and diversity have been studied [3,8,51,52]. The spatially explicit PGI data is particularly indispensable for these studies. For example, net primary productivity (NPP) is a critical parameter in vegetation growth and terrestrial ecosystem processes. If the effects of grazing on NPP (Net Primary Productivity) are ignored, the NPP estimated by the remote sensing approach may underestimate the actual NPP of grassland to some extent [8]. In order to accurately estimate the NPP at regional scale, GI information is needed.
This study proposes a new method for the estimation of the time-series GI with a finer spatial resolution based on the Space-Air-Ground integrated monitoring technology. The RGI average of the selected paddocks are close to their SGI (Figure 13), which indicates that the RGI estimated by the proposed method is statistically credible. The spatial statistical result shows that the average GI in Xiangdong Village was 2.9465 SU/ha, which was similar to the average GI (2.745 SU/ha) estimated by Wang et al. in Zoige County [8]. Since Xiangdong Village is located in the heart of the grazing area in Zoige County, it is normal for the average GI in Xiangdong Village to be slightly higher than the average GI of this county. Therefore, it is also reasonable to believe that the GI estimated in this paper is convincing. Furthermore, the spatial resolution of the GI estimated in this study is finer than that estimated by MODIS images, and can better describe the spatial heterogeneity of GI.
This study explored the estimation method of GI, but there is still a long way to go to generate accurate and practical time-series GI data at regional scale. Subsequent efforts need to be attempted in future research. (1) For ground monitoring, the battery capacity of the collars needs to be enhanced. At the same time, the time interval of the GPS collars needs to be shortened in the critical stages of vegetation growth because the GI at the stage of grass growth determines the yield and carrying capacity of the paddocks. More factors, such as unstable weather and terrain occlusion, would affect the positioning capability of GPS collars, resulting in the loss of some trajectory data. These factors need to be considered when applying to other regions, such as complex mountainous areas. New navigation technology will improve the acquisition ability and quality in the future. Additionally, corresponding algorithms for trajectory data with missing values need to be proposed. (2) For aerial monitoring, automatic recognition methods of the livestock and fences based on UAV images need to be explored, which would provide technical support for the commercial monitoring of GI. As the livestock inventory varied dynamically within one year due to the effects of calving, slaughtering, and selling, the frequency of UAV observation needs to be increased in these key periods. Multispectral or hyperspectral remote sensing images could be acquired directly based on the UAV platform, with the development of UAV technology, which would provide the possibility for directly quantitatively estimating GI based on UAV images. (3) For satellite monitoring, time-series satellite datasets with equal interval could be generated by using the multi-source data fusion method [53,54]. The time-series GI with equal intervals might better support the related application studies. (4) The proposed method was validated as a very good contribution to the estimation of GI with a finer spatial resolution in Zoige County. For estimating the GI at a regional scale, more details will be considered in the future, such as how to optimize the selection of sampling paddock to achieve the trade-off between cost and effectiveness and whether the experiment results of this study can be directly applied to other similar regions.

Conclusion
A novel estimation method for GI on the Zoige Plateau is proposed in this paper based on the Space-Air-Ground integrated monitoring technology. The time-series GI data in Xiangdong Village were produced by the proposed method. For ground monitoring, the complete trajectory data of eight monitored herds were recorded by using GPS collars in this study. Three UAV observation experiments were conducted, synchronized with ground experiments. The livestock inventory and area of the monitored paddocks were extracted from the UAV images. Based on the KDE algorithm, the trajectory data, and interpretation data of UAV images, the estimation model of PGI was constructed. The PGI of the eight monitored paddocks were produced by the proposed model. Taking these produced PGI as training data, an estimation model of RGI was constructed by using the time-series satellite remote sensing images and random forest regression algorithm. The time-series GI data in Xiangdong Village with a spatial resolution of 10 m was produced based on the RGI estimation model. According to the accuracy assessment of RGI, the generated GI data can reflect the spatial-temporal heterogeneity of actual GI, with an MAE of 0.9301 and r 2 of 0. 8573. More importantly, the proposed method provides a new idea for obtaining the actual GI on the ground and the time-series GI data. In the future, more efforts need to be conducted, such as the generation of time-series satellite datasets, developing an automatic recognition method for UAV images, and acquisition of trajectory data for the whole year. The generated time-series GI data would provide data support for the formulation of policies and plans related to the sustainable development of animal husbandry.