A grazing intensity inversion method based on statistical analysis and deep learning method

Abstract

The present application belongs to the technical field of grassland ecological remote sensing monitoring and data inversion, and provides a grazing intensity inversion method based on statistical analysis and deep learning method, which comprises: acquiring and preprocessing dynamic variable data and static variable data required by an LSTM model based on a remote sensing cloud platform; constructing an NDVI dataset by using a maximum synthesis method; calculating the time stability and spatial stability of each period of NDVI data, and screening the non-grazing window; inverting the potential NDVI of each pixel; quantifying the residual signal between each pixel potential NDVI and actual NDVI as grassland utilization rate; according to the standard sheep unit daily forage quota, the residual signal is mapped into a standardized carrying capacity with physical meaning, and finally a high-precision long-time-series grazing intensity dataset is obtained.

Description

Technical Field

[0001] This invention belongs to the field of grassland ecological remote sensing monitoring and data inversion technology, specifically involving a method for inverting grazing intensity based on statistical analysis and deep learning. Background Technology

[0002] Grasslands, as the dominant vegetation type in arid regions, play an irreplaceable and vital role in regional ecosystems. They not only provide crucial ecosystem services such as carbon sequestration, water regulation, and soil conservation, but also serve as a vital guarantee for economic development and food security in pastoral areas. However, against the backdrop of global climate change, rising temperatures, altered precipitation patterns, and the impacts of human activities have placed grassland ecosystems under increasingly complex pressures. Unreasonable grazing practices have brought severe challenges to grassland ecosystems and livestock production.

[0003] Grazing is one of the most important and continuous human activities in grassland pastoral areas. Its intensity is an important indicator for evaluating the impact of grazing activities on grasslands. Changes in grazing intensity directly affect key ecological processes such as regional grassland productivity, vegetation resilience, and the balance between grass and livestock. Accurate and scientific inversion of grazing intensity and quantification of its spatial distribution characteristics are of great significance for grassland ecological security assessment.

[0004] Remote sensing technology, with its advantages of large scale, strong continuity, and rapid update rate, has provided important support for the inversion of regional grazing intensity. Numerous studies have constructed grassland utilization indicators based on Normalized Difference Vegetation Index (NDVI), Net Primary Productivity (NPP), and aboveground biomass, and combined these with machine learning, linear regression, and other methods to characterize the degree of grassland disturbance caused by grazing, yielding a series of research results. With the rapid development of deep learning methods, their advantages in modeling complex nonlinear relationships are becoming increasingly prominent, providing a new technical path for grazing intensity inversion and potentially further improving the accuracy and reliability of the inversion results.

[0005] Existing technical solutions include: (1) Traditional manual monitoring methods: mainly relying on field quadrat surveys, herder statistical records, and empirical estimation to assess grassland utilization; (2) Methods based on remote sensing indices and statistical models: using low-to-medium resolution satellite data (such as MODIS) or medium-to-high resolution satellite data (such as Sentinel-2, Landsat) to construct time-series vegetation indices (such as NDVI), and estimating grazing intensity through rate of change analysis, threshold discrimination, or linear regression models; (3) Data-driven methods based on machine learning: combining multi-source remote sensing data and environmental variables, using machine learning models such as random forests and support vector machines to invert grazing intensity. The existing technology has the following drawbacks: While grazing intensity monitoring technology has a certain application foundation in grassland resource assessment and management, it still faces several limitations. Traditional manual monitoring methods mainly rely on quadrat surveys, herder statistical records, and empirical parameter estimations, playing a crucial role in early grassland ecological research. However, these methods are highly dependent on manual operation for data acquisition, resulting in high costs and susceptibility to factors such as the subjective experience of surveyors, the accuracy of herder reports, and insufficient representativeness of sampling points. This leads to significant uncertainties in the accuracy and consistency of monitoring results. Furthermore, their observation methods are mostly intermittent surveys with low temporal resolution, making it difficult to capture continuous dynamic changes in grazing activities and effectively characterize the spatial heterogeneity of grazing pressure at the regional scale.

[0006] Furthermore, existing methods for grazing intensity inversion based on remote sensing and machine learning still face the problem of insufficient sample support. Grazing activities are highly spatially dispersed, temporally dynamic, and subject to variations in human management. Actual studies often lack clearly defined grazing and non-grazing samples that have undergone rigorous field validation, making it difficult to provide high-quality, distinguishable, and generalizable training data for machine learning models. Especially in large-scale grassland monitoring scenarios, the boundaries between grazing and non-grazing areas are often unclear, and some areas may exhibit complex situations such as light grazing, seasonal grazing, or rotational grazing, leading to uncertainty and confounding in sample labeling. Insufficient training samples, unclear category classification, and weak sample representativeness directly affect the model's ability to learn grazing disturbance characteristics, thereby reducing the reliability of grazing intensity identification and estimation results.

[0007] Simple machine learning methods and empirical statistical regression models are often insufficient to fully reveal the complex nonlinear relationships between vegetation indices and environmental variables such as climate, topography, soil, water conditions, and human grazing activities, making it impossible to accurately quantify the degree of disturbance to grassland vegetation caused by grazing activities. Existing studies on grazing intensity are mostly limited to annual scales, failing to reflect the monthly dynamic changes in grazing activities during the growing season, thus hindering timely guidance and regulation of herders' grazing behavior. Due to the significant lag in monitoring results, it is difficult to promptly detect the risk of overgrazing or grassland degradation in local areas, potentially exacerbating grassland ecosystem degradation. Summary of the Invention

[0008] To address the aforementioned technical problems, this invention provides a grazing intensity inversion method based on statistical analysis and deep learning, thereby resolving the issues in the prior art. The technical solution adopted by this invention is as follows: A method for inverting grazing intensity based on statistical analysis and deep learning includes: Step S1: Based on the remote sensing cloud platform, acquire long-term satellite remote sensing data covering the target area as the original image, along with the corresponding meteorological data and background environmental data; after preprocessing the original image, synthesize it into NDVI data, and then construct a continuous monthly NDVI time series dataset. Step S2: Based on the continuous monthly NDVI time series dataset, screen potential ungrazed windows with high vegetation growth, stable temporal performance, and continuous spatial distribution. Step S3: Construct a convolutional long short-term memory network model. Based on the continuous monthly NDVI time series dataset obtained in step S1 and the potential ungrazed window obtained in step S2 as input, simulate the potential NDVI of the target area after iterative optimization. Step S4: Perform grazing intensity inversion. Based on potential NDVI and actual NDVI data, convert NDVI loss into vegetation dry weight foraging to complete the quantitative inversion of grazing intensity. Step S5: Perform grazing intensity inversion pixel by pixel, construct a spatial dataset of grazing intensity in the target area, and perform spatial mapping of grazing intensity in the target area.

[0009] Furthermore, in step S1, the preprocessing of the original image includes spatial cropping and quality control processing to remove cloud contamination and abnormal pixels; the maximum value synthesis method is used to synthesize NDVI data month by month.

[0010] Furthermore, step S2 includes: Step S2-1: Divide the continuous monthly NDVI time series dataset of the target area into several spatial windows according to the preset window size, and define stable high NDVI pixels. Step S2-2: Calculate the temporal stability of the spatial window based on the stable high NDVI pixels identified in Step S2-1; Step S2-3: Calculate the spatial clustering of the spatial window based on the stable high NDVI pixels identified in Step S2-1; Step S2-4: For each valid time phase, when the temporal stability and spatial clustering of its corresponding spatial window simultaneously meet the preset threshold, the spatial window is considered to have passed the stability test in that time phase and is regarded as a potential ungrazed window.

[0011] Furthermore, in step S2-1, defining stable high NDVI pixels includes: For any spatial window, the NDVI pixel values ​​within the spatial window are read sequentially, and invalid or abnormally low-value pixels are removed; in each time phase, the high quantile threshold of the valid NDVI values ​​of the spatial window is calculated, and pixels not lower than the high quantile threshold are identified as high NDVI pixels of that time phase. The frequency of each pixel being identified as a high NDVI pixel in multiple effective time phases is counted. If the frequency of any pixel being repeatedly identified as a high NDVI pixel in an effective time phase is not lower than the preset core frequency threshold, it is classified into the long-term stable core area and regarded as a stable high NDVI pixel.

[0012] Furthermore, in step S2-2, the formula for calculating the temporal stability of the spatial window is: in, This represents the time stability score for the t-th time phase; C indicates the high-value area in the current period; C indicates the long-term stable core area.

[0013] Furthermore, in steps S2-3, the spatial clustering calculation process of the spatial window includes: Determine whether there are other high NDVI pixels in the neighborhood of each high NDVI pixel, and count the proportion of high NDVI pixels that meet the neighborhood clustering condition to obtain the spatial clustering index. The formula is: in, This represents the spatial clustering score of the t-th time phase. This indicates the total number of high NDVI pixels at that time phase. This indicates the number of high NDVI pixels with neighborhood clustering relationships.

[0014] Furthermore, step S4 includes: Step S4-1: Based on the potential NDVI obtained in step S3 and the actual NDVI obtained in step S1, calculate the residual between them pixel by pixel: in, This indicates NDVI in a state of no grazing or potential. Indicates the actual observed NDVI. This indicates the NDVI loss signal caused by grazing interference; Step S4-2: Combining the grazing gradient experiment data from the stations, the NDVI loss signal is converted into actual vegetation dry weight. Based on the standard daily feed intake quota for sheep, a mapping function is constructed: Where C represents the dry weight of vegetation consumed; f (·) represents the mapping function obtained by fitting the experimental data of grazing gradient at the station.

[0015] The present invention has the following beneficial effects: (1) Improve the objectivity and reliability of identifying potential ungrazed samples; Existing methods for grazing intensity inversion typically rely on manual sampling, herder statistics, or empirical judgment to determine grazing and ungrazing samples. This results in high sample acquisition costs and is susceptible to subjective human factors, insufficient sample representativeness, and uncertainty in sample labels. This invention, based on multi-temporal NDVI data, utilizes temporal and spatial stability indicators to screen potential ungrazing windows. It can identify long-term stable high-value areas of vegetation based on the characteristics of remote sensing data itself, thereby reducing reliance on manual samples and improving the objectivity and consistency of ungrazing sample screening.

[0016] (2) Improve the accuracy of potential NDVI simulation; This invention introduces a Convolutional Long Short-Term Memory (ConvLSTM) model for simulating potential NDVI. This model can not only learn the time lag relationship and cumulative effect between climatic factors such as temperature and precipitation and vegetation growth, but also extract spatial correlation features between neighboring pixels through its convolutional structure. Compared to traditional linear regression, empirical statistical models, and general machine learning methods, this invention can better characterize the nonlinear relationships and spatiotemporal coupling features in the vegetation growth process, thereby improving the estimation accuracy of potential NDVI.

[0017] (3) Achieve effective separation between grazing interference signals and climate-driven signals; In existing technologies, NDVI changes are often simultaneously influenced by multiple factors such as climate change, soil moisture, topographic conditions, and human grazing activities, making it difficult to directly distinguish grazing disturbance signals. This invention constructs a potential NDVI under environmental constraints and performs difference analysis with the actual observed NDVI. The residual signal between the potential and actual states is used as a characterization of grazing disturbance intensity, thereby effectively reducing the interference of climate fluctuations on grazing intensity inversion and enhancing the ecological interpretability of the inversion results.

[0018] (4) Improve the spatiotemporal resolution of grazing intensity monitoring; Existing grazing intensity assessments primarily rely on administrative region statistics or annual grassland utilization rates, which fail to reflect the dynamic changes in grazing activities during the growing season and are insufficient to characterize the spatial heterogeneity of local areas. This invention, based on long-term remote sensing data and pixel-scale modeling methods, enables dynamic inversion of grazing intensity at the monthly and pixel-scale levels. It reveals more precisely the temporal changes and spatial variations in grazing activities, providing technical support for identifying local overgrazing areas, seasonal grazing pressure changes, and grassland degradation risks.

[0019] (5) To realize the physical dimension expression of grazing intensity; Traditional remote sensing inversion results often remain at the level of vegetation index differences, grassland utilization rates, or relative intensity levels, lacking direct physical meaning for grassland management. This invention combines grazing gradient experimental data from stations to map the NDVI residual signal to the dry weight of vegetation consumed per unit area, and further introduces the standard daily feed intake quota per sheep unit to convert it into a standardized carrying capacity (SSR). Therefore, the inversion results can not only represent the intensity of grazing disturbance but also be expressed dimensionally in terms of "sheep units / square meter / month," facilitating grassland-livestock balance accounting, grassland carrying capacity assessment, and management policy formulation.

[0020] (6) Reduce the cost of large-scale monitoring and improve application efficiency; This invention primarily relies on remote sensing data, cloud computing platforms, and deep learning models to invert grazing intensity, eliminating the need for large-scale, high-frequency manual surveys and significantly reducing manpower, time, and financial costs. Furthermore, this method possesses good automation capabilities and regional scalability, making it suitable for rapid monitoring, long-term assessment, and dynamic updating of grazing intensity on large-scale grasslands. Attached Figure Description

[0021] Figure 1 This is a schematic diagram of the process of the present invention; Figure 2 Example image for identifying potential ungrazed windows; Figure 3 Example diagram of potential NDVI simulation distribution; Figure 4 This is an example figure showing the fitting of potential NDVI and actual NDVI in 90 ungrazed sample points. Detailed Implementation

[0022] The following will be described in conjunction with embodiments of the present invention. Figures 1-4 The technical solutions in the embodiments of the present invention will be clearly and completely described. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Unless otherwise specified, the technical means used in the embodiments are conventional means well known to those skilled in the art.

[0023] like Figures 1-4 This invention proposes a method for inverting grazing intensity based on statistical analysis and deep learning, comprising the following steps: Step S1, Monthly NDVI synthesis and variable data acquisition and processing: Leveraging remote sensing cloud platforms (such as Google Earth Engine, GEE), long-term time-series satellite remote sensing data (such as MODIS data) covering the target area are acquired, along with corresponding meteorological data (dynamic variables such as temperature and precipitation) and background environmental data (static variables such as soil type, grassland type, and elevation). Spatial cropping and quality control processing are performed on the original images to remove cloud contamination and anomalous pixels. Maximum Value Compositing (MVC) is used to synthesize NDVI data monthly, constructing a continuous monthly NDVI time-series dataset as the core input variable for the model. After unifying the spatial resolution of dynamic and static variables, they are aligned with the NDVI dataset to jointly constitute the multi-source input feature set for subsequent models.

[0024] In practice, other remote sensing indices that can reflect vegetation growth or changes in surface biomass can be used to replace NDVI, such as Enhanced Vegetation Index (EVI), Leaf Area Index (LAI), and Net Primary Productivity (NPP); however, NDVI data is preferred in this invention.

[0025] Step S2, Identification of potential ungrazed windows: In this embodiment, the potential ungrazed window identification module is used to screen potential ungrazed windows with high vegetation growth, stable temporal performance and continuous spatial distribution from the study area based on the continuous monthly NDVI time series dataset.

[0026] Step S2 includes: Step S2-1, Long-term stable high-value pixel statistics: First, the continuous monthly NDVI time series dataset of the target area is divided into several spatial windows according to a preset window size. In a specific embodiment, the window size is set to 25×25 pixels, and each window is an independent evaluation unit. For any window, the NDVI pixel values ​​within the window are read time-by-time, and invalid or abnormally low-value pixels are removed. Subsequently, the high quantile threshold of the effective NDVI values ​​of the window is calculated in each time phase. The 95th percentile is used as the high value threshold, and pixels not lower than this threshold are determined as high NDVI pixels of that time phase. The expression is: Where H_t represents the set of high NDVI pixels in the t-th time phase, NDVI_t(p) represents the NDVI value of pixel p in the t-th time phase, and P95(NDVI_t) represents the 95th percentile of the effective NDVI values ​​of the window in the t-th time phase. To calculate whether the high NDVI region within the window has spatiotemporal stability, the frequency of each pixel being identified as a high NDVI pixel in multiple effective time phases is statistically analyzed, and the long-term stable core region is determined accordingly. If the frequency of a pixel being repeatedly identified as a high NDVI pixel in an effective time phase is not lower than a preset core frequency threshold, it is included in the long-term stable core region. In a specific embodiment, the core frequency threshold is 0.40. The long-term stable core region can be represented as: Where C represents the long-term stable core region, F(p) represents the frequency at which pixel p becomes a high NDVI pixel in multiple effective time phases, and γ represents the core frequency threshold.

[0027] Step S2-2, Temporal Stability Calculation: Based on the stable high-value pixels identified in S2-1, calculate the temporal stability for each period window. The temporal stability index characterizes the degree of overlap between the current high NDVI pixels and the long-term stable core area. Its calculation method is as follows: ,in, The value represents the high-value region in the current period, and TS_t represents the temporal stability score of the t-th time phase. The larger the value, the more concentrated the high NDVI region of that time phase appears in the long-term stable position, indicating that the window has strong temporal stability.

[0028] Step S2-3, Spatial Stability Calculation: Based on the stable high-value pixels identified in S2-1, calculate the spatial clustering for each period window. The spatial stability index characterizes whether high NDVI pixels in the current time phase exhibit a spatially continuous distribution. Specifically, determine whether there are other high NDVI pixels in the neighborhood of each high NDVI pixel, and statistically analyze the proportion of high NDVI pixels that satisfy the neighborhood clustering condition to obtain the spatial clustering index. in, This represents the spatial clustering score of the t-th time phase. This represents the total number of high NDVI pixels at that time phase, where N_cluster,t represents the number of high NDVI pixels with neighborhood clustering relationships. The larger this value, the more likely high NDVI pixels are to form continuous patches rather than isolated noise points.

[0029] Step S2-4, Screening of Potential Ungrazed Windows: For each valid time phase, if both its temporal stability index and spatial clustering index simultaneously meet preset thresholds, the window is considered to have passed the stability test for that time phase. In a specific embodiment, the temporal stability threshold is set to 0.50, and the spatial clustering threshold is set to 0.60, that is: When the proportion of the number of time phases that pass the stability test to the number of valid time phases is not less than 0.7, the window is considered to be a potential ungrazed window.

[0030] Step S3, Potential NDVI Simulation: Construct a Convolutional Long Short-Term Memory (ConvLSTM) network model. This network structure includes a spatial convolution module (extracting shared spatial features of adjacent habitat pixels) and an LSTM temporal module (extracting the lag and cumulative effects of climate and vegetation growth in long-term series). Based on the input variables obtained in Step S1 and the potential ungrazed window obtained in Step S2 as data input samples, the model is iteratively optimized to finally simulate the pixel-by-pixel potential NDVI of the target area.

[0031] The Convolutional Long Short-Term Memory (ConvLSTM) network model constructed in this invention comprises a spatial convolution module (extracting shared spatial features of adjacent habitat pixels) and an LSTM temporal module (extracting the lag and cumulative effects of climate and vegetation growth in long-term series). This model uses variables such as climate factors, topographic factors, soil factors, and historical NDVI sequences obtained in step S1, as well as potential ungrazing window samples determined in step S2, as input data. Through joint learning of spatiotemporal features, it establishes a nonlinear mapping relationship between NDVI and environmental variables under grazing-free conditions.

[0032] The specific processing steps include: First, the input variables are reconstructed into a time series to form a spatiotemporal sample tensor centered on the target pixel and containing information about its neighboring pixels. Second, the samples are input into a spatial convolution module to extract spatial correlation and habitat similarity features between adjacent pixels through convolution operations. Then, the convolution features are input into a ConvLSTM time series module to extract the lag and cumulative effects of factors such as climate change, soil moisture, temperature, and precipitation on vegetation growth. Furthermore, a spatiotemporal attention module is introduced at the output of the ConvLSTM to assign adaptive weights to different time steps and different environmental factors to enhance the model's ability to identify key growing seasons, extreme climate periods, and major limiting factors. Finally, the potential NDVI of each pixel in the target area under an ungrazing scenario is output through a fully connected mapping layer.

[0033] The state update process of the ConvLSTM unit is as follows: In the formula, i and f t o t These represent the input gate, forget gate, and output gate, respectively. Indicates candidate memory state; C t and H t represents the cell state and hidden state, respectively; * represents convolution operation; ⊙ represents element-wise multiplication; σ represents the sigmoid activation function; tanh represents the hyperbolic tangent activation function; W and b are the model's learnable parameters and bias terms, respectively. The input gate i is determined by the static entity variable X. s Controlled and unchanging over time within the same pixel or window, it is used to characterize the inherent differences in the responses of different grassland entities to external climate drives under long-term natural backgrounds. (Forget gate f) t and output gate o t Controlled by dynamic climate variables and the hidden state of the previous time step, it is used to characterize the impact of current and previous climate conditions on grassland growth memory retention, renewal, and output. Step S4, grazing intensity inversion, includes: Step S4-1, Residual Signal Quantization: Based on the potential NDVI obtained in step S3 and the actual NDVI obtained in step S1, calculate the residual between them pixel by pixel. ,in, This indicates NDVI in a state of no grazing or potential. Indicates the actual observed NDVI; This represents the NDVI loss signal caused by grazing interference, i.e., the residual between the potential NDVI and the actual NDVI.

[0034] Step S4-2, Physical Dimension Mapping: Combining the grazing gradient experimental data from the stations, NDVI loss is converted into actual vegetation dry weight. Based on the introduced standard sheep daily feed quota (such as physical / biological parameters like the daily dry matter intake required per standard sheep), a mapping function is constructed: Where C represents the dry matter mass of vegetation consumed per unit area, i.e., the dry weight of vegetation consumed, and the unit can be kg / m². 2 or g / m 2 Δt represents the residual between the potential NDVI and the actual NDVI; f(·) represents the mapping function fitted based on the experimental data of grazing gradient at the station. This mapping function can be constructed using linear, nonlinear, or piecewise functions based on the characteristics of the measured data.

[0035] Furthermore, a standard daily feed intake quota for sheep is introduced, converting the aforementioned dry weight of vegetation per unit area into a standardized carrying capacity. Let the standard daily dry matter intake for sheep be... If the statistical period is L, then the standardized carrying rate (SSR, unit: sheep units / square meter / month) can be expressed as: .

[0036] Step S5, Spatial mapping of grazing intensity: Grazing intensity is inverted pixel by pixel to construct a high-precision long-term target area grazing intensity spatial dataset, and spatial mapping of grazing intensity in the target area is performed.

[0037] Specifically, based on the pixel-by-pixel vegetation grazing dry weight and / or standardized carrying capacity obtained in step S4, the grazing intensity results for each time phase are spatially rasterized and mapped to construct a long-term spatial dataset of grazing intensity in the target area. The standardized carrying capacity (SSR) obtained in step S4 is used as a quantitative indicator of grazing intensity. SSR represents the standard sheep carrying capacity per unit area of ​​grassland within the statistical period. A higher SSR value indicates higher vegetation grazing and greater grazing intensity in the area corresponding to that pixel during that period; a lower SSR value indicates weaker grazing pressure in the area corresponding to that pixel. For any pixel (p) and any time phase (t) within the target area, the standardized carrying capacity calculated in step S4 is denoted as... It is defined as the grazing intensity of that pixel in that time phase. ,Right now: in, This represents the grazing intensity of pixel (p) at time phase (t). This represents the standardized carrying capacity of pixel (p) at time phase (t), calculated from the dry weight of vegetation grazing. Thus, the standardized carrying capacity, which has physical dimensions, is directly converted into a pixel-scale grazing intensity value.

[0038] Furthermore, based on the unified spatial reference, spatial resolution, and temporal scale of the target area, the (GI(p,t)) value of each pixel is written to the corresponding raster position, forming a monthly-scale grazing intensity raster dataset. For grazing intensity at the growing season or annual scale, the comprehensive grazing intensity can be obtained by summarizing the grazing intensity results from multiple time phases: in This represents the cumulative grazing intensity of pixel p during the statistical period.

[0039] The above embodiments are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Any modifications, alterations, alterations, or substitutions made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for inverting grazing intensity based on statistical analysis and deep learning, characterized in that, include: Step S1: Based on the remote sensing cloud platform, acquire long-term satellite remote sensing data covering the target area as the original image, along with the corresponding meteorological data and background environmental data; after preprocessing the original image, synthesize it into NDVI data, and then construct a continuous monthly NDVI time series dataset. Step S2: Based on the continuous monthly NDVI time series dataset, screen potential ungrazed windows with high vegetation growth, stable temporal performance, and continuous spatial distribution. Step S3: Construct a convolutional long short-term memory network model. Based on the continuous monthly NDVI time series dataset obtained in step S1 and the potential ungrazed window obtained in step S2 as input, simulate the potential NDVI of the target area after iterative optimization. Step S4: Perform grazing intensity inversion. Based on potential NDVI and actual NDVI data, convert NDVI loss into vegetation dry weight foraging to complete the quantitative inversion of grazing intensity. Step S5: Perform grazing intensity inversion pixel by pixel, construct a spatial dataset of grazing intensity in the target area, and perform spatial mapping of grazing intensity in the target area.

2. The grazing intensity inversion method based on statistical analysis and deep learning as described in claim 1, characterized in that, In step S1, the preprocessing of the original image includes spatial cropping and quality control processing to remove cloud contamination and abnormal pixels; the maximum value synthesis method is used to synthesize NDVI data month by month.

3. The grazing intensity inversion method based on statistical analysis and deep learning as described in claim 1, characterized in that, Step S2 includes: Step S2-1: Divide the continuous monthly NDVI time series dataset of the target area into several spatial windows according to the preset window size, and define stable high NDVI pixels. Step S2-2: Calculate the temporal stability of the spatial window based on the stable high NDVI pixels identified in Step S2-1; Step S2-3: Calculate the spatial clustering of the spatial window based on the stable high NDVI pixels identified in Step S2-1; Step S2-4: For each valid time phase, when the temporal stability and spatial clustering of its corresponding spatial window simultaneously meet the preset threshold, the spatial window is considered to have passed the stability test in that time phase and is regarded as a potential ungrazed window.

4. The grazing intensity inversion method based on statistical analysis and deep learning as described in claim 3, characterized in that, In step S2-1, stable high NDVI pixels are defined as follows: For any spatial window, the NDVI pixel values ​​within the spatial window are read sequentially, and invalid or abnormally low-value pixels are removed; in each time phase, the high quantile threshold of the valid NDVI values ​​of the spatial window is calculated, and pixels not lower than the high quantile threshold are identified as high NDVI pixels of that time phase. The frequency of each pixel being identified as a high NDVI pixel in multiple effective time phases is counted. If the frequency of any pixel being repeatedly identified as a high NDVI pixel in an effective time phase is not lower than the preset core frequency threshold, it is classified into the long-term stable core area and regarded as a stable high NDVI pixel.

5. The grazing intensity inversion method based on statistical analysis and deep learning as described in claim 3, characterized in that, In step S2-2, the formula for calculating the temporal stability of the spatial window is: in, This represents the time stability score for the t-th time phase; C indicates the high-value area in the current period; C indicates the long-term stable core area.

6. The grazing intensity inversion method based on statistical analysis and deep learning as described in claim 3, characterized in that, In steps S2-3, the spatial clustering calculation process of the spatial window includes: Determine whether there are other high NDVI pixels in the neighborhood of each high NDVI pixel, and count the proportion of high NDVI pixels that meet the neighborhood clustering condition to obtain the spatial clustering index. The formula is: in, This represents the spatial clustering score for the t-th time phase. This indicates the total number of high NDVI pixels at that time phase. This indicates the number of high NDVI pixels with neighborhood clustering relationships.

7. The grazing intensity inversion method based on statistical analysis and deep learning as described in claim 1, characterized in that, Step S4 includes: Step S4-1: Based on the potential NDVI obtained in step S3 and the actual NDVI obtained in step S1, calculate the residual between them pixel by pixel: in, This indicates NDVI in a state of no grazing or potential. Indicates the actual observed NDVI. This indicates the NDVI loss signal caused by grazing interference; Step S4-2: Combining the grazing gradient experiment data from the stations, the NDVI loss signal is converted into actual vegetation dry weight. Based on the standard daily feed intake quota for sheep, a mapping function is constructed: Where C represents the dry weight of vegetation consumed; f (·) represents the mapping function obtained by fitting the experimental data of grazing gradient at the station.

Images