Method for measuring nutrient balance of natural grassland grazing system

Abstract

This invention discloses a method for calculating the nutrient balance of natural grassland grazing systems, belonging to the field of agricultural ecology and grassland management technology. To address the problems of low accuracy and difficulty in guiding precise management in existing nutrient balance calculations for grazing systems, this invention divides spatial units by coupling vegetation spatial heterogeneity index and livestock utilization intensity index; determines the diet composition and deducts vegetation mechanical damage to calculate net nutrient output; constructs a spatiotemporal dynamic layer of excrement deposition; uses a soil animal promotion coefficient-corrected decomposition model to calculate effective nutrient return; and performs spatial overlay analysis of nutrient output and effective nutrient return to calculate the nutrient balance index of each spatial unit, generating decisions for rest grazing, rotational grazing, or supplemental feeding. This invention enables accurate calculation of the nutrient balance of grazing systems, providing a scientific basis for the sustainable utilization and precise management of natural grasslands.

Description

Technical Field

[0001] This invention relates to the field of agricultural ecology and grassland management technology. More specifically, this invention relates to a method for calculating the nutrient balance of a natural grassland grazing system. Background Technology

[0002] The nutrient balance of natural grassland grazing systems directly affects grassland productivity and ecological stability. Accurately measuring nutrient input and output during grazing is crucial for developing reasonable grazing management measures. However, due to the complexity and openness of natural grassland grazing systems, precise calculation of their nutrient balance has long faced numerous technical challenges.

[0003] Spatially, natural grasslands exhibit significant spatial heterogeneity, with vegetation type, soil conditions, and topographic features varying considerably across different locations. This results in inconsistent nutrient supply capacity and livestock utilization intensity at different points within the same grazing unit. Traditional nutrient balance calculation methods typically treat the entire grazing area as a homogeneous unit, but this extensive spatial approach fails to reflect the true situation, leading to insufficient representativeness of the calculation results.

[0004] Existing methods for estimating nutrient output have the following main shortcomings: Firstly, the measurement of livestock feed intake often relies on short-term, fixed-point observations or estimation models based on physiological needs, making it difficult to obtain continuous and accurate feed intake data during actual grazing. This is particularly true in natural grasslands where livestock have large activity ranges and dispersed grazing behaviors, where the coverage and timeliness of traditional observation methods are limited. Secondly, while the composition of livestock diets is a key factor influencing nutrient output, the vegetation composition of natural grasslands is complex, and livestock selectively consume plant species under different seasons and habitat conditions, making it difficult to accurately determine the actual plant species and proportions ingested by livestock using traditional botanical surveys and microscopic analysis. Furthermore, in addition to grazing, livestock also cause mechanical damage to vegetation through trampling and resting during grazing. Although this vegetation is not ingested, it still results in nutrient loss, which traditional methods typically overlook, leading to systematic biases in nutrient output estimation.

[0005] Existing technologies also face significant limitations in measuring nutrient return. Excrement is the primary carrier of nutrient return in grazing systems, but its spatial distribution on grasslands is extremely uneven, exhibiting clear spatial aggregation characteristics influenced by factors such as livestock behavior, topography, and water source location. Traditional methods often estimate nutrient return based on average excrement per unit area, neglecting the spatial heterogeneity of actual excrement distribution, leading to significant discrepancies between calculated and actual results. More importantly, not all nutrients in excrement can promptly enter the soil-plant cycle and be converted into usable nutrients for plants; their decomposition and transformation process is influenced by a combination of meteorological conditions such as temperature and humidity, as well as soil biological activity. Existing studies typically use fixed conversion coefficients or simple empirical formulas for estimation, failing to fully consider the dynamic changes in meteorological conditions and the influence of biological factors such as soil animals, making it difficult to guarantee the timeliness and accuracy of the calculation results.

[0006] In summary, existing methods have shortcomings in handling spatial heterogeneity, accurately monitoring livestock behavior, precisely determining diet composition, deducting vegetation mechanical damage, characterizing the spatial distribution of excrement, and estimating the effective conversion rate of nutrients. They are unable to truly reflect the actual nutrient balance of grazing systems and thus restrict the scientific nature of grazing management decisions. Summary of the Invention

[0007] One object of the present invention is to solve at least the above-mentioned problems and to provide at least the advantages that will be described later.

[0008] Another objective of this invention is to provide a method for calculating the nutrient balance of a natural grassland grazing system. This method aims to solve the technical problems of existing technologies, which fail to accurately calculate nutrient output and return due to neglecting the spatial heterogeneity of natural grasslands and the dynamic nature of livestock grazing behavior. Furthermore, it is difficult to dynamically correct the process of excrement nutrients actually being converted into available soil nutrients, thus failing to truly reflect the nutrient balance of the grazing system and making it difficult to support precise management decisions.

[0009] To achieve these objectives and other advantages of the present invention, a method for calculating the nutrient balance of a natural grassland grazing system is provided, comprising the following steps: S1. Acquire multi-source remote sensing data and grazing trajectory data of the target area, and divide the natural grassland into multiple spatial units by coupling the vegetation spatial heterogeneity index and the livestock utilization intensity index. S2. By identifying feeding behavior through GPS locators and jaw movement sensors worn by livestock, combined with the diet composition determined by DNA macrobarcode sequencing, and deducting the amount of mechanical damage to vegetation caused by trampling and resting, the total net nutrient intake of livestock in each spatial unit is calculated as the nutrient output. S3. Based on the location and time of excretion events recorded by GPS locators and mandibular motion sensors, and combined with the spatial distribution of excrement patches identified by UAV thermal imaging, a spatiotemporal dynamic layer of excrement deposition is constructed. Representative excrement samples were collected for decomposition bag experiments, and the changes in nitrogen and oxygen stable isotopes were measured simultaneously to determine the baseline conversion coefficient of excrement nutrients to soil available nutrients during the grazing cycle. At the same time, the population density of soil animals was investigated, and the promoting coefficient of soil animal activity on the decomposition rate of excrement was fitted through indoor controlled experiments. The baseline conversion coefficient is input into the decomposition model coupled with micrometeorological data, and the decomposition rate in the model is dynamically corrected using the soil animal promotion coefficient. Combined with the spatiotemporal dynamic layer of excrement deposition, the effective nutrient return amount that actually enters the soil-plant cycle in each spatial unit during the grazing cycle is calculated. S4. Perform spatial overlay analysis on the nutrient output of S2 and the effective nutrient return of S3, calculate the nutrient balance index of each spatial unit, identify nutrient surplus or deficit areas, and generate control decisions for rest grazing, rotational grazing, or supplementary feeding.

[0010] Preferably, in the method for calculating the nutrient balance of a natural grassland grazing system, step S1 involves dividing the natural grassland into multiple spatial units, specifically including: The Normalized Differential Vegetation Index (NDVI) of forage grass was obtained by using multispectral data from UAVs or satellite remote sensing data. A regression model between NDVI and forage grass biomass was established by combining field quadrat surveys, and a spatial distribution map of forage grass biomass was generated. Slope data was extracted based on the digital elevation model (DEM), and a comprehensive grazing suitability index was constructed by combining the distance to water source and the type of livestock. Using GIS spatial analysis tools, the spatial distribution map of pasture biomass was overlaid with the comprehensive index of grazing suitability for analysis, and the study area was divided into multiple homogeneous spatial units according to the set grid size.

[0011] Preferably, in the method for calculating the nutrient balance of the natural grassland grazing system, step S2, calculating the total net nutrient intake of livestock in each spatial unit, specifically includes: S21. Synchronize the data from the GPS locator worn by livestock with the data from the mandibular motion sensor, identify feeding, walking, and resting behaviors through the characteristics of the mandibular motion waveform, and count the feeding time and walking frequency of each spatial unit. S22. Set up representative quadrats in the study area, collect samples of major forage grasses to determine the nutrient content of crude protein, neutral detergent fiber and other nutrients, and establish a forage grass nutrient database. S23. Differentiate between foraging loss and trampling loss, and calculate the net intake separately: Feed loss ,in Let j be the feeding time of the j-th livestock in unit i. Let j be the number of feed intakes per unit time for the j-th livestock. The average nutrient content of pasture grass consumed by livestock; Trampling damage ,in S step,ij Let j be the walking cadence of the j-th livestock in the i-th unit. L step For the average step size, W trample To trample on the effective width, B i Let i be the pasture biomass of the i-th unit. α i The trampling damage coefficient is dynamically corrected using synchronously monitored soil moisture content. α i = α 0×(1+ k × θ i ), α 0 is the base damage coefficient. θ i Let be the soil volumetric water content of the i-th unit. k Moisture sensitivity coefficient; Net nutrient output i = F i - T i .

[0012] Preferably, in the method for calculating the nutrient balance of the natural grassland grazing system, when calculating the effective nutrient return in step S3, the decomposition and conversion coefficient of excrement is further dynamically corrected by combining meteorological data during the grazing cycle, specifically including: S31. Based on the excretion events recorded by the GPS locator and the mandibular motion sensor, generate the time series and location information of the excretion events in each spatial unit, and combine the area of ​​excrement patches identified by UAV thermal imaging to establish a spatiotemporal dynamic layer of excrement deposition. S32. Collect representative fresh excrement samples, determine the content of total nitrogen, total phosphorus, total potassium, and organic matter, and calculate the nutrient content per unit area of ​​excrement. C dung ; S33. Establish a baseline conversion coefficient through a decomposition bag experiment over at least one grazing cycle. R 0Nylon biodegradable bags containing a known mass of excrement were placed in representative plots. Residual nutrient content was measured at multiple time points, and nitrogen and oxygen stable isotope analyses were performed simultaneously. The proportion of nutrients actually entering the soil-plant cycle was calculated based on isotope mass balance and used as the baseline conversion coefficient. R 0 ; S34. Simultaneously record daily meteorological data during the decomposition experiment, including accumulated temperature. T sum and effective precipitation P eff Based on the historical multi-year average meteorological data of the region, a dynamic correction model for the decomposition transformation coefficient and meteorological factors was established: ; in , This refers to the multi-year average accumulated temperature and effective precipitation. a , b These are the empirical coefficients fitted through regression analysis; S35. Based on the spatiotemporal dynamic layer of excrement deposition, for each excretion event, apply the corresponding dynamic decomposition and transformation coefficient R according to its deposition time, and sum to obtain the effective nutrient return of each spatial unit: ; in A patch,ik Let be the area of ​​the k-th excrement patch in the i-th spatial unit. D dung It refers to the dry matter weight of excrement per unit area.

[0013] Preferably, in the method for calculating the nutrient balance of the natural grassland grazing system, step S4, generating a decision on rest grazing, rotational grazing, or supplemental feeding, specifically includes: S41. Calculate the nutrient balance index of each spatial unit according to the following formula: B i =Output i / (Supply soil,i +Supply deposition,i ); Output i The nutrient output of the i-th spatial unit is calculated in step S2; Supply soil,i For the available nutrients in the i-th spatial unit, the available nutrient content is determined by collecting soil samples and calculating the available nutrient content in combination with soil bulk density; Supply deposition,i The total amount of this nutrient input to atmospheric dry and wet deposition was obtained by actual measurement using deposition collectors deployed in the study area or by referencing data from nearby meteorological stations. S42. Set the nutrient balance index threshold: when Bi When the range is ∈ [0.8, 1.2], it is determined to be a nutrient balance zone, B i >1.2 indicates a nutrient deficiency zone, B i <0.8 indicates a nutrient surplus zone; S43. For nutrient-deficient areas, generate a grazing rest area map by combining the remaining biomass of pasture based on remote sensing inversion and topographic conditions; for nutrient-surplus areas, generate a priority grazing area map for the next grazing cycle by combining water source distribution and slope conditions; for areas with continuous nutrient deficiency that livestock must pass through, generate a map of recommended locations for supplementary feeding points.

[0014] Preferably, the method for calculating the nutrient balance of the natural grassland grazing system is reviewed and corrected after each grazing cycle. A correlation analysis is performed between the current nutrient balance calculation results and the livestock production performance data of the same period. If the calculation results of two consecutive grazing cycles are contrary to the actual production performance trend, then the soil animal promotion coefficient in the decomposition model described in step S3, or the water sensitivity coefficient in step S23, will be adjusted. k Or the dynamic correction of model parameters in step S34 a , b Make corrections.

[0015] Preferably, in the method for calculating the nutrient balance of the natural grassland grazing system, when dividing spatial units in step S1, a livestock spatial memory preference factor is further introduced for correction, specifically including: Historical GPS trajectory data from at least one complete grazing cycle was collected, and foraging hotspot distribution maps were extracted using kernel density analysis to generate a livestock spatial memory preference index. M i ; Livestock Spatial Memory Preference Index M i As a correction weight, the comprehensive index of grazing suitability in step S1 is weighted and corrected: ; in I suit,i The index represents the overall grazing suitability index before revision. I suit,i ′ The corrected index is λ, which is the preference intensity coefficient, obtained by regression analysis of historical grazing trajectory data and vegetation spatial distribution. Spatial units were divided using a modified comprehensive index of grazing suitability.

[0016] Preferably, in the nutrient balance calculation method for the natural grassland grazing system, when calculating the trampling loss in step S23, the spatial heterogeneity of trampling damage is further considered. Specifically, this includes: connecting the GPS trajectory points of livestock to generate a walking path, using the path as the center line, dynamically adjusting the buffer zone width according to the livestock walking speed, establishing an inner buffer zone, a middle buffer zone, and an outer buffer zone respectively, and assigning different damage weights w1, w2, and w3, where w1>w2>w3; Calculate the weighted trampling damage coefficient for each spatial unit based on the proportion of areas falling into different buffer zones: ; in α i The coefficient of damage caused by trampling. A ik Let be the area within the i-th spatial cell that falls into the k-th level buffer. A i Let i be the total area of ​​the i-th spatial unit. w k Let be the damage weight of the k-th level buffer; Weighted trampling damage coefficient α i ′ Calculate the amount of damage caused by trampling.

[0017] The present invention has at least the following beneficial effects: First, this invention divides natural grasslands into multiple spatial units by coupling the vegetation spatial heterogeneity index and the livestock utilization intensity index, thus solving the problem of insufficient representativeness of measurement results caused by traditional methods that treat grazing areas as homogeneous units and ignore spatial heterogeneity. By combining GPS trajectory and grazing behavior recognition with jaw movement sensors, it can accurately distinguish between grazing and trampling losses, significantly improving the accuracy of nutrient output measurement and providing a reliable data foundation for the refined management of nutrient balance in grazing systems.

[0018] Secondly, this invention introduces DNA macrobarcoding technology to accurately determine the composition of livestock diets, and combines it with a vegetation mechanical damage deduction mechanism to overcome the systematic bias of traditional methods that rely on empirical estimation and ignore trampling losses. By constructing a two-factor calculation model of "feeding loss + trampling loss" and introducing a dynamic soil moisture correction coefficient for trampling damage, it achieves dynamic and accurate quantification of actual nutrient intake by livestock, improving the accuracy and adaptability of nutrient output measurement.

[0019] Third, this invention constructs a dynamic layer of excrement deposition based on the spatiotemporal distribution of excretion events and UAV thermal imaging identification, solving the problem of traditional methods neglecting the spatial heterogeneity of excrement. Simultaneously, by combining decomposition bag experiments, nitrogen and oxygen stable isotope analysis, soil animal promotion coefficient correction, and a dynamic correction decomposition model based on meteorological data, it achieves for the first time a dynamic simulation of the entire process of effective nutrient return from excrement, significantly improving the timeliness and accuracy of nutrient return measurement.

[0020] Fourth, this invention accurately identifies nutrient surplus and deficit areas through spatial overlay analysis and nutrient balance index calculation, and generates control decisions such as rest grazing, rotational grazing, or supplementary feeding, solving the problem that existing technologies cannot support precise management decisions. A verification and correction mechanism is introduced to dynamically adjust model parameters based on livestock production performance, ensuring consistency between the calculation results and actual production trends, and improving the method's adaptability, robustness, and application value.

[0021] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Detailed Implementation

[0022] The present invention will be further described in detail below with reference to embodiments, so that those skilled in the art can implement it based on the description.

[0023] It should be understood that terms such as “having,” “comprising,” and “including” as used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0024] It should be noted that, unless otherwise specified, the experimental methods described in the following implementation plan are all conventional methods, and the reagents and materials described are all commercially available unless otherwise specified.

[0025] This invention provides a method for calculating the nutrient balance of a natural grassland grazing system, comprising the following steps: S1. Acquire multi-source remote sensing data and grazing trajectory data of the target area, and divide the natural grassland into multiple spatial units by coupling the vegetation spatial heterogeneity index and the livestock utilization intensity index. S2. By identifying feeding behavior through GPS locators and jaw movement sensors worn by livestock, combined with the diet composition determined by DNA macrobarcode sequencing, and deducting the amount of mechanical damage to vegetation caused by trampling and resting, the total net nutrient intake of livestock in each spatial unit is calculated as the nutrient output. S3. Based on the location and time of excretion events recorded by GPS locators and mandibular motion sensors, and combined with the spatial distribution of excrement patches identified by UAV thermal imaging, a spatiotemporal dynamic layer of excrement deposition is constructed. Representative excrement samples were collected for decomposition bag experiments, and the changes in nitrogen and oxygen stable isotopes were measured simultaneously to determine the baseline conversion coefficient of excrement nutrients to soil available nutrients during the grazing cycle. At the same time, the population density of soil animals was investigated, and the promoting coefficient of soil animal activity on the decomposition rate of excrement was fitted through indoor controlled experiments. The baseline conversion coefficient is input into the decomposition model coupled with micrometeorological data, and the decomposition rate in the model is dynamically corrected using the soil animal promotion coefficient. Combined with the spatiotemporal dynamic layer of excrement deposition, the effective nutrient return amount that actually enters the soil-plant cycle in each spatial unit during the grazing cycle is calculated. S4. Perform spatial overlay analysis on the nutrient output of S2 and the effective nutrient return of S3, calculate the nutrient balance index of each spatial unit, identify nutrient surplus or deficit areas, and generate control decisions for rest grazing, rotational grazing, or supplementary feeding.

[0026] Specifically, the multi-source remote sensing data mentioned in step S1 includes UAV multispectral imagery or satellite remote sensing data, and grazing trajectory data is collected by GPS trackers worn by livestock. First, the Normalized Differential Vegetation Index (NDVI) of pasture is obtained through UAV multispectral imagery (flight altitude 100-200m, spatial resolution 0.1-0.5m) or satellite remote sensing data (such as Sentinel-2, resolution 10m). Combined with field quadrat surveys (at least 3 quadrats per hectare, quadrat area 1m×1m), the aboveground biomass of pasture within the quadrats is measured, and a regression model between NDVI and pasture biomass is established, thereby generating a spatial distribution map of pasture biomass for the entire study area, with biomass units of kg / m³. 2 or g / m 2 Secondly, slope data was extracted based on the Digital Elevation Model (DEM), and combined with distance to water sources (calculated through GIS buffer analysis; water sources include rivers, lakes, and artificial drinking water points) and livestock species to construct a comprehensive grazing suitability index. The comprehensive grazing suitability index can be calculated using a weighted overlay method: I suit =w1 × normalized biomass value + w2 × normalized slope value + w3 × normalized water source distance value, where the weights w1, w2, and w3 are determined through the analytic hierarchy process (AHP) or expert experience, for example, 0.5, 0.3, and 0.2 respectively, and the normalized values ​​are all converted to the range of 0-1. Finally, using GIS spatial analysis tools, the spatial distribution map of pasture biomass is overlaid with the comprehensive index of grazing suitability for analysis. The study area is divided into multiple homogeneous spatial units according to the set grid size (e.g., 10m × 10m). The selection of the grid size needs to comprehensively consider the grassland type, terrain complexity, and monitoring accuracy requirements, and is usually set to 5-20m. In practice, if the terrain is flat and the vegetation is uniform, the grid size can be appropriately increased (e.g., 20m); if the terrain is fragmented and the vegetation heterogeneity is high, the grid size should be reduced (e.g., 5m) to ensure homogeneity within the unit.

[0027] In step S2, the total net nutrient intake of livestock in each spatial unit is calculated. First, the data from the GPS locator worn by the livestock (sampling frequency 1Hz, recording latitude and longitude coordinates and time) is synchronized with the data from the mandibular motion sensor (sampling frequency 10-50Hz, recording triaxial acceleration). The mandibular motion sensor uses a pre-set machine learning classification model to identify behaviors such as feeding, walking, and resting, with an accuracy rate exceeding 95%. The feeding time (minutes) and walking cadence (steps / minute) of each livestock in each spatial unit are recorded. Next, representative quadrats are established in the study area, and samples of major forage grasses are collected to determine the nutrient content, including crude protein and neutral detergent fiber. Near-infrared spectroscopy or conventional chemical analysis methods are used to establish a forage nutrient database. Then, feeding loss and trampling loss are distinguished, and the net intake is calculated separately: Feeding Loss F. i =∑ j t ij ×r j ×c, where t ij r represents the feeding time (in minutes) of the j-th livestock in unit i. j denoted as , the number of bites (bites / minute) of the j-th livestock, obtained by counting chewing movements using a mandibular sensor; denoted as , the average nutrient content of the forage consumed by the livestock (g / bite), determined by collecting forage samples from areas marked by the livestock's grazing. Trampling loss T. i =∑ j S step,ij ×L step ×W trample ×B i ×α i S step,ij Let L be the walking frequency (steps / minute) of the j-th livestock in the i-th unit. step The average stride length (meters per step) can be estimated from the animal's size and GPS trajectory; for cattle, it is generally 0.6-1.0m, and for sheep, it is 0.4-0.6m. W trample The effective trampling width (in meters), i.e., the lateral range within which livestock cause damage to vegetation when walking, was determined through field observations. For cattle, it is generally 0.4-0.6 meters, and for sheep, it is 0.3-0.4 meters. i The pasture biomass of unit i (kg / m³) 2 ), obtained from the inversion in step S1, α i The trampling damage coefficient is dynamically corrected using synchronously monitored soil moisture content: α i =α0×(1+k×θ i ), where α0 is the basic damage coefficient, determined through fence comparison experiments, and is generally taken as 0.1-0.3, θ iThe soil volumetric water content (%) of unit i is obtained through soil moisture sensor or remote sensing inversion, and k is the moisture sensitivity coefficient, which is experimentally fitted and ranges from 0.5 to 2.0. Net nutrient uptake output. i =F i -T i The unit is kg or g, and can be calculated separately according to different nutrients (N, P, K).

[0028] In step S3, when calculating the effective nutrient return, the time series and location information of excretion events within each spatial unit are first generated based on excretion events recorded by the GPS locator and mandibular motion sensor. Simultaneously, UAV thermal imaging flights are conducted periodically during the grazing cycle. Because fresh excrement generates heat due to microbial fermentation, it appears as distinct hotspots in the thermal imaging. The outlines of excrement patches are extracted using an image segmentation algorithm, and the area is calculated to establish a spatiotemporal dynamic layer of excrement deposition. Each excretion event records the deposition time, location coordinates, and patch area. Secondly, representative fresh excrement samples are collected, and the contents of total nitrogen, total phosphorus, total potassium, and organic matter are measured to calculate the nutrient content C per unit area of ​​excrement. dung (g / m 2 Simultaneously, the dry matter weight D of the excrement was measured. dung (kg / m 2 Then, a baseline conversion coefficient R0 was established through a decomposition bag experiment over at least one grazing cycle: nylon decomposition bags containing a known mass of excrement were placed in representative plots, and residual nutrient content was measured at multiple time points. Simultaneously, nitrogen and oxygen stable isotope analyses were performed on the samples. The proportion of nutrients actually entering the soil-plant cycle was calculated based on isotope mass balance and used as the baseline conversion coefficient R0. Daily meteorological data, including accumulated temperature T, were recorded concurrently during the decomposition experiment. sum and effective precipitation P eff Based on the historical multi-year average meteorological data of the region, a dynamic correction model for the decomposition transformation coefficient and meteorological factors was established: ,in , The values ​​represent the multi-year average accumulated temperature and effective precipitation, with a and b being empirical coefficients fitted through regression analysis. Simultaneously, soil animal population density was investigated, and the promoting coefficient β of soil animal activity on the decomposition rate of excrement was fitted through indoor controlled experiments. β was multiplied into the decomposition rate term in the decomposition model to correct the baseline transformation coefficient. Finally, based on the spatiotemporal dynamic layer of excrement deposition, for each excretion event, the corresponding dynamic decomposition transformation coefficient R was applied according to its deposition time, and the effective nutrient return of each spatial unit was summed to obtain: A patch,ik The area (m²) of the k-th excrement patch in the i-th spatial unit. 2 ).

[0029] In step S4, the decision to implement rest grazing, rotational grazing, or supplemental feeding is generated. First, the nutrient balance index of each spatial unit is calculated according to the following formula: , where output i The nutrient output of the i-th spatial unit, Supply soil,i For the available nutrients in the i-th spatial unit, the available nutrient content is determined by collecting soil samples and calculating the available nutrient content in combination with soil bulk density. deposition,i The total amount of this nutrient input for atmospheric dry and wet deposition was obtained through actual measurements using deposition collectors deployed in the study area or by referencing data from nearby meteorological stations. Secondly, a nutrient balance index threshold was set: when B... i When the range is ∈ [0.8, 1.2], it is determined to be a nutrient balance zone, B i >1.2 indicates a nutrient deficiency zone, B i Areas with a nutrient surplus of less than 0.8 are considered nutrient-deficient areas. Finally, for nutrient-deficient areas, a grazing rest area map is generated by combining the remaining biomass of pasture grasses based on remote sensing inversion and topographic conditions; for nutrient-surplus areas, a priority grazing area map for the next grazing cycle is generated by combining water source distribution and slope conditions; and for areas that are continuously deficient and must be traversed by livestock, a map of recommended locations for supplementary feeding points is generated.

[0030] In the above technical solution, spatial units are divided by coupling the vegetation spatial heterogeneity index and the livestock utilization intensity index, which solves the problem of measurement bias caused by treating grazing areas as homogeneous units in traditional methods. By identifying grazing behavior using GPS and mandible sensors, determining dietary composition using DNA macrobarcoding, and deducting mechanical damage, the problems of inaccurate feed intake estimation and systematic bias are solved. By constructing a spatiotemporal dynamic layer of excrement deposition, combined with decomposition bag experiments, stable isotope analysis, soil animal promotion coefficient correction, and meteorological dynamic correction, the problems of uneven spatial distribution of excrement and difficulty in dynamically correcting nutrient transformation are solved. By generating regulatory decisions through spatial overlay analysis and nutrient balance index calculation, the problem of measurement results being insufficient to support precise management is solved. This invention achieves accurate measurement of nutrient balance in grazing systems, providing a scientific basis for the sustainable utilization and precise management of natural grasslands.

[0031] In another technical solution, step S1 divides the natural grassland into multiple spatial units, specifically including: The Normalized Differential Vegetation Index (NDVI) of forage grass was obtained by using multispectral data from UAVs or satellite remote sensing data. A regression model between NDVI and forage grass biomass was established by combining field quadrat surveys, and a spatial distribution map of forage grass biomass was generated. Slope data was extracted based on the digital elevation model (DEM), and a comprehensive grazing suitability index was constructed by combining the distance to water source and the type of livestock. Using GIS spatial analysis tools, the spatial distribution map of pasture biomass was overlaid with the comprehensive index of grazing suitability for analysis, and the study area was divided into multiple homogeneous spatial units according to the set grid size.

[0032] Specifically, the Normalized Differential Vegetation Index (NDVI) of pasture is first acquired using UAV multispectral or satellite remote sensing data. If using UAV multispectral imagery, the flight altitude is typically set to 100-120m, with a forward overlap of at least 80% and a lateral overlap of at least 70% to ensure image stitching quality. If using satellite remote sensing data, Sentinel-2 satellite imagery can be selected, with a spatial resolution of 10m and a revisit period of 5 days, which can meet the dynamic monitoring needs of the grazing cycle. The formula for calculating NDVI is NDVI = ( ρ NIR - ρ R ) / ( ρ NIR + ρ R ),in ρ NIR For near-infrared reflectivity, ρ R This refers to the reflectivity in the red light band.

[0033] Simultaneously with acquiring NDVI data, a field quadrat survey was conducted. Stratified random sampling was used for quadrat placement. Based on the preliminary NDVI classification results, the study area was divided into high, medium, and low cover layers, with quadrats randomly placed within each layer. Each quadrat had an area of ​​1m × 1m, and the center coordinates were recorded using differential GPS. All forage within the quadrats was harvested at ground level and dried to constant weight in a 65℃ oven. The dry matter biomass (g / m²) was recorded. The biomass data was paired with the corresponding NDVI values ​​to establish a regression model, such as a linear model. B = a ×NDVI+ b or exponential model B = a × eb ×NDVI, select the model with the highest coefficient of determination R². Apply the model to the NDVI raster layer of the entire study area to generate a spatial distribution map of pasture biomass.

[0034] Secondly, slope data is extracted based on the Digital Elevation Model (DEM). The DEM can be derived from SRTM (30m resolution) or a high-precision DEM generated by UAV photogrammetry. Slope is extracted using slope analysis tools in GIS software, expressed in degrees or percentages. Simultaneously, the distance to water sources, including rivers, lakes, and artificial drinking water points, is calculated using the Euclidean distance tool, generating a distance raster map.

[0035] A comprehensive grazing suitability index was constructed based on livestock species. Different livestock species showed varying preferences for slope and water source: cattle were highly dependent on water sources and preferred gentle slopes; sheep had moderate water dependence and could adapt to gentle slopes; goats were less dependent on water sources and could adapt to steep slopes. The comprehensive grazing suitability index was calculated using a weighted summation method. I suit = w 1×B norm + w 2×S norm + w 3×D norm ; Among them B norm S norm D norm These are the normalized values ​​(range 0-1) for biomass, slope, and distance to water source, respectively. Slope and distance to water source need to be inversely normalized (higher values ​​indicate better suitability). Weight w 1, w 2, w 3. The values ​​are determined using the analytic hierarchy process (AHP). For example, for a cattle grazing system, values ​​of 0.6, 0.2, and 0.2 can be used.

[0036] Finally, GIS spatial analysis tools were used to overlay the spatial distribution map of pasture biomass with the comprehensive grazing suitability index, dividing the study area into multiple homogeneous spatial units according to a set grid size. The grid size was typically set to 5-20m: 15-20m for typical grasslands with uniform vegetation and flat terrain; and 5-10m for meadows with high vegetation heterogeneity and fragmented terrain. The average pasture biomass and grazing suitability index within each grid unit were calculated using a zoning statistical tool and used as the attribute value for that unit, thus completing the spatial unit division.

[0037] In the above technical solution, by coupling the spatial distribution of vegetation biomass with the grazing suitability index, a refined characterization of the spatial heterogeneity of natural grasslands is achieved, solving the problem of measurement bias caused by treating the entire grazing area as a homogeneous unit in traditional methods. A biomass distribution map is generated by combining multi-source remote sensing data with field surveys, ensuring both spatial coverage and accuracy. By constructing a comprehensive grazing suitability index that considers livestock species preferences, the unit division is made more consistent with actual grazing utilization patterns. The gridded overlay analysis method makes the nutrient supply capacity and utilization intensity of each spatial unit closer to reality, providing a reasonable unit basis for the subsequent accurate calculation of nutrient output and return.

[0038] In another technical solution, step S2 involves calculating the total net nutrient intake of livestock in each spatial unit, specifically including: S21. Synchronize the data from the GPS locator worn by livestock with the data from the mandibular motion sensor, identify feeding, walking, and resting behaviors through the characteristics of the mandibular motion waveform, and count the feeding time and walking frequency of each spatial unit. S22. Set up representative quadrats in the study area, collect samples of major forage grasses to determine the nutrient content of crude protein, neutral detergent fiber and other nutrients, and establish a forage grass nutrient database. S23. Differentiate between foraging loss and trampling loss, and calculate the net intake separately: Feed loss ,in Let j be the feeding time of the j-th livestock in unit i. Let j be the number of feed intakes per unit time for the j-th livestock. The average nutrient content of pasture grass consumed by livestock; Trampling damage ,in S step,ij Let j be the walking cadence of the j-th livestock in the i-th unit. L step For the average step size, W trample To trample on the effective width, B i Let i be the pasture biomass of the i-th unit. α i The trampling damage coefficient is dynamically corrected using synchronously monitored soil moisture content. α i = α 0×(1+ k × θ i ), α 0 is the base damage coefficient. θ i Let be the soil volumetric water content of the i-th unit. k Moisture sensitivity coefficient; Net nutrient output i = F i - T i .

[0039] Specifically, in step S21, the GPS locator worn by the livestock is set to a sampling frequency of 1Hz to record the livestock's latitude and longitude coordinates and timestamps in real time. The mandibular motion sensor uses a triaxial accelerometer with a sampling frequency of 20-50Hz, worn on the livestock's mandible to record the acceleration waveforms of chewing, swallowing, and rumination. During data preprocessing, the GPS data and mandibular sensor data are synchronized using a unified time reference, with a time registration accuracy required to be within ±0.1 seconds.

[0040] Behavior recognition is achieved through machine learning algorithms. First, manually observed and labeled behavioral data is collected as training samples, including typical behaviors such as feeding (continuous chewing movements with regular high-frequency oscillations in the waveform), walking (regular head swaying coupled with gait frequency), and resting (head lowered with a flat waveform). Temporal features (mean, variance, peak value), frequency domain features (spectral energy distribution), and wavelet features are extracted. A classification model is constructed using random forest or support vector machine, achieving an accuracy of over 95%. The trained model is applied to continuously collected sensor data to identify the real-time behavioral status of livestock and associate it with GPS location by timestamp. Through GIS spatial analysis, the cumulative feeding time (t) for each livestock within each spatial unit is statistically analyzed. ij (minutes) and walking cadence S step,ij (Steps per minute, i.e., the number of steps per minute while walking).

[0041] In step S22, the collection of forage samples must cover the main vegetation types and topographic features of the study area. Before the start of the grazing cycle, stratified random sampling is used to establish quadrats, each 1m × 1m in size, with no fewer than 10 quadrats for each vegetation type. The aboveground parts of the main forage species within each quadrat are collected, placed in cloth bags, and brought back to the laboratory. The samples are dried to constant weight in a 65℃ oven, pulverized, and sieved through a 1mm sieve. The Kjeldahl nitrogen determination method is used to determine the crude protein content, and the Van der Waals detergent fiber analysis method is used to determine the neutral detergent fiber, acid detergent fiber, and lignin content. All samples are measured three times, and the average value is taken. A forage nutrient database is established, including forage species, sampling locations, sampling times, crude protein, neutral detergent fiber, and other indicators. The average nutrient content c of livestock feed needs to be weighted and calculated in conjunction with the dietary composition determined by DNA metabarcoding in step S2: c = ∑ m p m ×c m , where p m Let c be the mass proportion of the m-th type of forage in the diet. m This refers to the nutrient content (such as crude protein content) of this type of forage.

[0042] In step S23, the calculation of feed loss is based on feeding time and feeding rate. The number of bites per unit time (rjrj) is directly counted using mandibular sensors: during feeding behavior, each chewing action is recorded as one bite, and the number of chews per minute is calculated. Feeding rates vary among different livestock; cattle typically have 50-70 bites / minute, and sheep 40-60 bites / minute, requiring individual calibration. The average step length (L) is also considered. stepGPS trajectory estimation: Select continuous straight-line walking trajectory segments, divide the trajectory length by the number of steps (obtained by multiplying the step frequency and time during walking behavior), and take the average of multiple measurements. The average stride length for cattle is generally 0.6-1.0m, and for sheep it is 0.4-0.6m. Effective trampling width (Wtrample): Determined through field observation: Randomly select sampling points along the grazing path, measure the lateral range of vegetation damage caused by livestock walking, and take the average. For cattle, it is generally 0.4-0.6m, and for sheep it is 0.3-0.4m. Pasture biomass (B) i Extracted from the spatial distribution map of pasture biomass generated in step S1, in units of kg / m³. 2 or g / m 2 .

[0043] Trampling damage coefficient α i Dynamic correction is one of the key innovations of this technical solution. The basic damage coefficient α0 is determined through a fenced comparative experiment: multiple fenced closed plots are set up in the study area (to exclude grazing interference), while grazing plots are set up at corresponding locations outside the fences. Vegetation cover and biomass changes are monitored regularly, and the proportion of additional damage caused by grazing is calculated. The experiment needs to cover different seasons and soil moisture conditions to obtain the basic damage coefficient α0, which generally ranges from 0.1 to 0.3. Soil volumetric water content θ i Soil moisture is acquired through multiple sources: soil moisture sensors (such as TDR probes) are deployed at key locations for continuous monitoring. Combined with Sentinel-1 SAR satellite data, regional soil moisture distribution is retrieved to generate the average soil moisture content (%) within spatial units. The moisture sensitivity coefficient k is obtained through experimental fitting: simulated trampling experiments are conducted under different soil moisture conditions to determine the degree of vegetation damage. A regression relationship between the damage coefficient and soil moisture is established. The value of k generally ranges from 0.5 to 2.0. The higher the soil clay content and the finer the texture, the larger the k value, indicating a more significant amplification effect of moisture on trampling damage.

[0044] Net nutrient output i =F i -T i The calculations need to differentiate between different nutrient types. For nitrogen, c is calculated using crude protein content (crude protein contains approximately 16% nitrogen); for other nutrients such as phosphorus and potassium, corresponding nutrient databases and calculation formulas need to be established. The final result is the nutrient output of each spatial unit during the grazing cycle, expressed in kg or g.

[0045] In the above technical solution, high-precision behavior recognition and spatial positioning enable the separation and quantification of livestock feed intake and trampling losses. By employing GPS and a mandible sensor working in tandem, the spatiotemporal trajectory of livestock behavior is accurately correlated with the spatial distribution of vegetation for the first time. Machine learning algorithms identify behavioral states such as feeding, walking, and resting, avoiding the subjectivity and discontinuity of manual observation. Dynamic soil moisture correction of the trampling damage coefficient addresses the problem of traditional methods neglecting the influence of environmental conditions on the trampling effect. Combining the dietary composition determined by DNA macrobarcoding to calculate the weighted average nutrient content makes the calculation of feed intake losses more closely reflect the actual feed consumption of livestock. These combined technical approaches significantly improve the accuracy and dynamic adaptability of nutrient output measurement, providing reliable output data for nutrient balance assessment of grazing systems.

[0046] In another technical solution, when calculating the effective nutrient return in step S3, the decomposition and conversion coefficient of excrement is further dynamically adjusted by combining meteorological data during the grazing cycle, specifically including: S31. Based on the excretion events recorded by the GPS locator and the mandibular motion sensor, generate the time series and location information of the excretion events in each spatial unit, and combine the area of ​​excrement patches identified by UAV thermal imaging to establish a spatiotemporal dynamic layer of excrement deposition. S32. Collect representative fresh excrement samples, determine the content of total nitrogen, total phosphorus, total potassium, and organic matter, and calculate the nutrient content per unit area of ​​excrement. C dung ; S33. Establish a baseline conversion coefficient through a decomposition bag experiment over at least one grazing cycle. R 0 Nylon biodegradable bags containing a known mass of excrement were placed in representative plots. Residual nutrient content was measured at multiple time points, and nitrogen and oxygen stable isotope analyses were performed simultaneously. The proportion of nutrients actually entering the soil-plant cycle was calculated based on isotope mass balance and used as the baseline conversion coefficient. R 0 ; S34. Simultaneously record daily meteorological data during the decomposition experiment, including accumulated temperature. T sum and effective precipitation P eff Based on the historical multi-year average meteorological data of the region, a dynamic correction model for the decomposition transformation coefficient and meteorological factors was established: ; in , This refers to the multi-year average accumulated temperature and effective precipitation. a , b These are the empirical coefficients fitted through regression analysis; S35. Based on the spatiotemporal dynamic layer of excrement deposition, for each excretion event, apply the corresponding dynamic decomposition and transformation coefficient R according to its deposition time, and sum to obtain the effective nutrient return of each spatial unit: ; in A patch,ik Let be the area of ​​the k-th excrement patch in the i-th spatial unit. D dung It refers to the dry matter weight of excrement per unit area.

[0047] Specifically, in step S31, the identification of defecation events is based on the coordinated operation of a mandibular motion sensor worn by the livestock and a GPS locator. The mandibular motion sensor identifies defecation actions through acceleration waveform characteristics—when the sensor detects a specific combination of tail movement and a stationary state, it determines it as a defecation event and automatically triggers the GPS to record the current time and location coordinates, with positioning accuracy reaching sub-meter level. This process requires no manual intervention and can be continuously recorded around the clock.

[0048] Drone thermal imaging flights are conducted regularly during the grazing cycle, with the frequency determined by the grazing intensity, generally 1-2 times per week. Flights are typically conducted in the early morning or late afternoon when the surface temperature difference is smaller, making the thermal characteristics of excrement more apparent. The drones are equipped with thermal imaging cameras, flying at altitudes of 50-80 meters, achieving a thermal imaging resolution of 640×512 pixels and a ground resolution of 0.1-0.2 meters. Fresh excrement generates heat due to microbial fermentation, appearing as distinct high-temperature patches in the thermal images, creating a significant temperature difference with the surrounding background, typically 3-8°C.

[0049] After acquiring thermal infrared images, image processing software was used for stitching and radiometric calibration, and a threshold segmentation method was employed to extract excrement patches. Specifically, the temperature distribution of the background region in the image was statistically analyzed, and the average background temperature T was calculated. 背景 And standard deviation σ, set temperature threshold T th =T 背景 +2σ identifies pixels above a threshold as candidate excrement patches. Connectivity analysis merges adjacent pixels into independent patches, eliminating those with excessively small areas (e.g., <0.1m). 2 The area of ​​each valid patch is calculated for the noise patches identified by thermal imaging. The patches identified by thermal imaging are spatiotemporally matched with excretion events recorded by GPS: the time window is set to ±1 day, and the distance window is set to ±5m. Successfully matched excretion events yield accurate patch area information. Finally, a spatiotemporal dynamic layer of excrement deposition is established, with layer attributes including: excretion event ID, livestock ID, deposition time, center coordinates, patch area, and corresponding spatial unit number.

[0050] In step S32, the collection of excrement samples must be sufficiently representative. Stratified sampling should be conducted based on the type of grazing livestock (cattle, sheep, horses, etc.) and season (spring, summer, autumn), with at least 10 fresh excrement samples collected from each combination (within 24 hours of excretion). Samples should be immediately refrigerated and transported back to the laboratory after collection. One portion will be used for dry matter content determination (dried at 105°C to constant weight), and the other portion will be frozen for nutrient analysis.

[0051] Nutrient determination methods are as follows: Total nitrogen was determined using the Kjeldahl method. After digestion with sulfuric acid and hydrogen peroxide, the sample was titrated by distillation using a Kjeldahl apparatus. Total phosphorus was determined using the molybdenum-antimony colorimetric method. After the digestion solution reacted with the molybdenum-antimony colorimetric reagent, the colorimetric concentration was measured at a wavelength of 700 nm. Total potassium was determined using a flame photometer. Organic matter was determined using the potassium dichromate titration method. Under external heating conditions, organic carbon was oxidized with potassium dichromate-sulfuric acid solution, and the remaining potassium dichromate was titrated with ferrous sulfate. Each sample was measured three times, and the average value was taken.

[0052] The method for calculating the nutrient content of excrement per unit area is as follows: D dung The dry matter weight of excrement per unit area (kg / m²) 2 D dung The method was determined through field measurements: a sampling frame of known area (e.g., 0.25m) was placed on the excrement patch. 2 or 0.5m 2 Collect all excrement in the collection box, weigh the fresh weight, and bring it back to the laboratory to dry to constant weight. Calculate the dry matter weight per unit area. For different livestock and different seasons, D dung The values ​​may differ and need to be measured separately.

[0053] In step S33, the decomposition bag experiment is a key technique for determining the baseline transformation coefficient. Decomposition bags are deployed in representative plots (covering major vegetation types and soil conditions), with no fewer than 30 bags at each sampling point to meet the needs of multiple sampling. The decomposition bags are made of nylon mesh, with the mesh size determined according to the experimental purpose: 0.05 mm mesh is used to exclude the influence of large soil animals; 2 mm mesh is used to allow small soil animals to enter, considering the effects of soil animals under natural conditions. Each bag contains 100g of fresh excrement sample (accurately recording the initial dry weight and nutrient content), is sealed, laid flat on the ground, and secured with metal nails to prevent displacement.

[0054] Sampling time points are set according to the length of the grazing cycle, generally 7-8 time points, such as days 7, 15, 30, 45, 60, 75, and 90. At each time point, 3-5 decomposition bags are randomly collected, surface adhering substances are carefully removed, and the residual dry weight and nutrient content are measured to calculate the nutrient residue rate. The decay curve of the residue rate over time reflects the decomposition dynamics of excrement.

[0055] Simultaneous nitrogen and oxygen stable isotope analysis was performed on both the original and decomposed samples. After grinding and sieving, the samples were analyzed using a stable isotope mass spectrometer to determine the δ¹⁸O values. 15 N and δ 18 O value. Based on the principle of isotope fractionation, lighter isotopes are preferentially released into the environment during the decomposition of excrement, leading to a relative enrichment of heavier isotopes in the residue. The proportion of nutrients actually entering the soil-plant cycle can be calculated using an isotope mass balance model. The simplified model is as follows: ; Where Δ is the empirical fractionation coefficient, calibrated through prior controlled experiments. This method can distinguish between leaching loss and mineralization transformation, accurately obtaining the proportion of nutrients truly converted into plant-available forms. A baseline dataset is established by correlating the R0 values ​​at each time point with the corresponding cumulative temperature and effective precipitation.

[0056] In step S34, meteorological data acquisition and model building are crucial. Small automatic weather stations are deployed at representative locations in the study area to record daily data such as temperature, precipitation, humidity, and wind speed. Accumulated temperature T... sum The calculation uses the daily accumulated temperature value where the daily average temperature is greater than 0℃, i.e., T. sum =∑max(T mean ,0). Effective precipitation P eff The actual infiltration rate, calculated by subtracting runoff loss and deep seepage from daily precipitation, can be expressed using the empirical formula P for grassland systems. eff =P×(1−α), where α is the runoff coefficient, which is determined based on the slope and vegetation cover. The greater the slope and the lower the vegetation cover, the greater the runoff coefficient, which is generally taken as 0.1-0.3.

[0057] Multi-year average accumulated temperature and effective precipitation The data were obtained from nearly 30 years of statistical data from nearby meteorological stations or calculated from long-term location monitoring data of the study area.

[0058] The dynamic correction model was constructed using multiple regression analysis. The actual conversion coefficient R0 at each time point obtained from the decomposition experiment was used as the dependent variable, and... and Using the least squares method to solve for the coefficients a and b as independent variables, and applying the constraint a + b = 1, the physical meaning of the model is guaranteed. For example, in a decomposition experiment in a typical grassland area of ​​Inner Mongolia, the fitted values ​​were a = 0.62 and b = 0.38, indicating that temperature has a slightly greater influence on the decomposition rate than precipitation in this region. The goodness of fit R0... 2 The model has a good predictive ability, reaching a value of 0.85.

[0059] In step S35, the calculation of effective nutrient return integrates all the above data results. For each discharge event, the grazing cycle stage is determined based on its deposition time, and the corresponding accumulated temperature and effective precipitation data are retrieved to calculate the dynamic correction coefficient R. If soil fauna survey data is available within this time period, it also needs to be multiplied by the soil fauna promotion coefficient β to obtain the final conversion coefficient R′ (the method for determining the soil fauna promotion coefficient is detailed in the relevant description in step S3).

[0060] The effective nutrient return of a single excrement plaque is calculated as follows: plaque area A patch,ik Multiply by the dry matter weight per unit area, D dung The total dry matter mass of the patch is obtained, then multiplied by the nutrient content C. dung The total nutrient input of the plaque is obtained, and then multiplied by the conversion coefficient R(t). k The actual effective nutrient return entering the soil-plant cycle is obtained. Summing over all excretion events yields the total effective nutrient return for each spatial unit throughout the entire grazing cycle: ; When performing calculations, ensure that the units are consistent and that the final result is in kg / ha or g / m³. 2 This indicates that it is convenient to compare with nutrient output later.

[0061] In the aforementioned technical solution, a multi-dimensional approach, including spatiotemporal dynamic layers of excrement, decomposition bag experiments, stable isotope tracing, and meteorological dynamic correction, was used to achieve, for the first time, a dynamic simulation of the entire process of effective nutrient return from excrement. Based on GPS and jaw sensors, the spatiotemporal location of excretion events was accurately recorded, and combined with UAV thermal imaging, the area of ​​excrement patches was precisely identified, solving the problem of traditional methods neglecting the spatial heterogeneity of excrement. A baseline conversion coefficient was established through decomposition bag experiments and stable isotope analysis, avoiding the arbitrariness of empirical estimations. A dynamic correction model for meteorological factors was introduced, allowing the conversion coefficient to be adjusted in real time according to the grazing process, improving the timeliness and accuracy of the calculation. Finally, the effective return of each unit was summed to provide reliable input data for the nutrient balance assessment of the grazing system.

[0062] In another technical solution, step S4 generates a decision on rest grazing, rotational grazing, or supplemental feeding, specifically including: S41. Calculate the nutrient balance index of each spatial unit according to the following formula: ; Output i The nutrient output of the i-th spatial unit is calculated in step S2; Supply soil,iFor the available nutrients in the i-th spatial unit, the available nutrient content is determined by collecting soil samples and calculating the available nutrient content in combination with soil bulk density; Supply deposition,i The total amount of this nutrient input to atmospheric dry and wet deposition was obtained by actual measurement using deposition collectors deployed in the study area or by referencing data from nearby meteorological stations. S42. Set the nutrient balance index threshold: when B i When the range is ∈ [0.8, 1.2], it is determined to be a nutrient balance zone, B i >1.2 indicates a nutrient deficiency zone, B i <0.8 indicates a nutrient surplus zone; S43. For nutrient-deficient areas, generate a grazing rest area map by combining the remaining biomass of pasture based on remote sensing inversion and topographic conditions; for nutrient-surplus areas, generate a priority grazing area map for the next grazing cycle by combining water source distribution and slope conditions; for areas with continuous nutrient deficiency that livestock must pass through, generate a map of recommended locations for supplementary feeding points.

[0063] Specifically, in step S41, the nutrient balance index is calculated based on the key nutrient (usually nitrogen), but the balance indices of nitrogen, phosphorus, and potassium can also be calculated separately as needed. i Obtained from step S2, in units of kg / ha or g / m³ 2 Soil nutrient supply soil,i The determination method is as follows: Before the start of the grazing cycle, 5-10 sampling points are set up in an "S" shape or checkerboard pattern within each spatial unit to collect mixed soil samples from the 0-20cm soil layer. After air-drying, the samples are sieved through a 2mm sieve to determine the contents of available nitrogen (alkaline available nitrogen), available phosphorus (Olsen-P or Bray-P), and available potassium (ammonium acetate extraction). Simultaneously, the soil bulk density is determined using the ring sampler method, and the soil mass per unit area (t / ha) is calculated. Soil nutrient availability = available nutrient content (mg / kg) × soil mass (t / ha) × unit conversion factor. For example, if the alkaline available nitrogen content is 50mg / kg and the bulk density is 1.2g / cm³, the soil nutrient availability is 1.2g / cm³. 3 Therefore, the soil mass of the 0-20cm soil layer is 1.2×10000×0.2=2400t / ha, and the available nitrogen amount is 50×2400 / 1000=120kg / ha.

[0064] Atmospheric dry and wet deposition Supply deposition,iData acquisition methods include: deploying precipitation and dust collectors at representative locations within the study area to collect samples monthly and determine nitrogen, phosphorus, and potassium inputs; or using long-term monitoring data from nearby meteorological / ecological stations. For wet deposition, precipitation is collected using rain gauges, and the concentrations of ammonium nitrogen, nitrate nitrogen, total phosphorus, and total potassium are measured, multiplied by the precipitation amount to obtain the input; dry deposition is collected using dust collection tanks, and the soluble nutrient content is measured. If no measured data is available within a spatial unit, Kriging interpolation can be used to extend point data to an area.

[0065] In step S42, the threshold setting is based on the nutrient balance theory of grassland ecosystems and practical management experience. The range of 0.8-1.2 takes into account measurement errors and natural fluctuations. When the index is below 0.8, it indicates that the input is much greater than the output, which may lead to soil nutrient accumulation and the risk of eutrophication of water bodies; when it is above 1.2, it indicates that the output is greater than the input, which will lead to a decline in soil fertility in the long term. The threshold can be fine-tuned according to grassland type and management objectives. For example, in degraded grassland restoration areas, the deficit threshold can be appropriately increased (e.g., 1.5) to promote vegetation restoration, while in water source protection areas, the surplus threshold can be decreased (e.g., 0.7) to prevent non-point source pollution.

[0066] In step S43, the generation of the grazing rest area map needs to comprehensively consider the distribution of deficient areas, remaining pasture biomass, and topographic conditions. Remaining pasture biomass can be obtained through inversion from remote sensing images at the end of the grazing cycle (using the same method as step S1). Combined with topographic slope data (DEM extraction), areas within the deficient area with remaining biomass below the grazing warning line (e.g., 300 kg / ha) and a slope greater than 15° are prioritized for grazing rest areas, as grazing on steep slopes easily exacerbates soil erosion. Patches meeting the criteria are extracted and merged in the GIS to generate a vector map of the grazing rest area, with suggested start and end times for the grazing rest period indicated.

[0067] The generation of priority grazing area maps targets nutrient surplus areas, taking into account water source distance and slope conditions. Areas within the surplus area that are close to water sources (<1km) and have gentle slopes (<10°) are suitable as priority grazing areas for the next cycle to promote nutrient recycling. Through buffer zone analysis and slope classification, areas meeting the criteria are marked as priority grazing areas, and the output is a recommended grazing sequence map.

[0068] The generation of supplementary feeding point recommendation maps targets areas with persistent deficiencies that livestock must traverse, such as watering routes and transhumance paths. These areas cannot be left to graze even when deficient, requiring supplementary feeding to compensate for nutrient deficiencies. First, units with deficiencies for two consecutive grazing cycles are identified. Then, areas with high livestock GPS track density (kernel density analysis) are selected, and the center points of these intersection areas or points evenly distributed along the paths are designated as recommended supplementary feeding points. Point attributes include coordinates, recommended supplementary feeding time (e.g., lambing season, dry season), and suggested supplementary feeding amount (calculated based on the degree of deficiency).

[0069] In the aforementioned technical solution, surplus and deficit areas are accurately identified through the nutrient balance index, and an intuitive control decision map is generated by combining spatial information, providing scientific and operable guidance for grazing managers. Input items are calculated using measured soil data and atmospheric deposition data, avoiding the bias of relying solely on model estimations; threshold settings consider ecological security and management needs, exhibiting flexibility and adaptability; the generated grazing rest, rotational grazing, and supplementary feeding maps can be directly imported into mobile GIS devices or printed as operational maps, facilitating implementation by frontline managers. This invention solves the problem that traditional methods struggle to support precise management decisions, achieving sustainable grassland utilization and optimized allocation of nutrient resources.

[0070] In another technical solution, a review and correction are performed after each grazing cycle: A correlation analysis is performed between the current nutrient balance calculation results and the livestock production performance data of the same period. If the calculation results of two consecutive grazing cycles are contrary to the actual production performance trend, then the soil animal promotion coefficient in the decomposition model described in step S3, or the water sensitivity coefficient in step S23, will be adjusted. k Or the dynamic correction of model parameters in step S34 a , b Make corrections.

[0071] Specifically, the verification and correction mechanism is a key element in achieving adaptive optimization and ensuring long-term reliability in this method. At the end of each grazing cycle (e.g., at the end of each winter or grazing season), two data collection tasks are carried out simultaneously: first, the nutrient balance index and summary data of each spatial unit calculated through steps S1-S4 during this cycle; and second, the production performance data of livestock during the same period.

[0072] The collection of livestock production performance data requires determining core indicators based on the type of grazing system. For meat livestock, core indicators include average daily weight gain (kg / day), slaughter weight (kg), and feed conversion ratio; for wool livestock, these include wool yield (kg / animal) and wool fiber diameter; and for breeding livestock, these include reproductive rate, lamb survival rate, and weaning weight. Data acquisition methods include: regular weighing records (weighing marked livestock monthly or quarterly), review of production records, and slaughter measurement data. To ensure data representativeness, each grazing unit or livestock group should have at least 30 valid samples.

[0073] Correlation analysis uses statistical methods to test the consistency between the measured results and production performance. The specific procedure is as follows: the study area is divided into several analysis units (which can be aggregated by spatial units or divided by livestock groups), and two variables are obtained from each unit—the nutrient balance index B. i The average value and livestock production performance index P i (e.g., average daily weight gain). Calculate the Pearson correlation coefficient r between the two: ; In theory, nutrient balance should be positively correlated with production performance: in nutrient-deficient areas (B>1.2), livestock face high feed pressure and limited nutrient availability, resulting in lower production performance; in nutrient-surplus areas (B<0.8), vegetation grows vigorously and feeds well, leading to higher production performance. If the calculation results conform to this rule, r should be positive and statistically significant (p<0.05).

[0074] The criteria for determining "opposite trends" need to consider both direction and degree. If one of the following two conditions occurs, and two consecutive grazing cycles persist, then the trend is considered opposite: The first scenario is a negative correlation: the calculation results show that the production performance in the nutrient surplus area is actually lower than that in the deficit area, with a negative correlation coefficient r and an absolute value greater than 0.3. For example, a certain area is calculated to be a nutrient surplus area (B=0.6), but the average daily weight gain of livestock is only 0.5 kg / day, lower than the 0.7 kg / day in the deficit area, contrary to expectations.

[0075] The second scenario is no correlation or weak correlation: the calculation results have no significant correlation with production performance, the correlation coefficient |r| < 0.1, and professional judgment indicates a significant contradiction. For example, long-term calculations show a persistent nutrient deficiency in a certain region, but livestock condition is good and production performance is stable, indicating that the actual nutrient supply may be underestimated.

[0076] Parameter correction is the core operation of verification and calibration, involving targeted adjustments to key parameters based on possible causes: If the trend is reversed, it may be related to a deviation in the calculation of nutrient return from excrement. In this case, the soil animal promotion coefficient β in the decomposition model in step S3 should be corrected first. The correction method is to conduct supplementary soil animal surveys and decomposition experiments in the next grazing cycle and refit the promotion coefficient. For example, if the original model uses β=1.2, but the supplementary experiment shows a significant increase in earthworm density and a faster decomposition rate, then β should be adjusted to 1.3-1.5. If supplementary experimental conditions are lacking, an optimization algorithm can be used to find the optimal value: with production performance data as the target, search for the β value within a certain range (e.g., 1.0-2.0) that best matches the calculated result with the actual trend.

[0077] If the trend is reversed, it may be related to a bias in the estimation of trampling loss. In this case, the moisture sensitivity coefficient k in step S23 should be corrected first. The correction method is as follows: add trampling damage observation plots under different soil moisture conditions, measure the relationship between the damage coefficient and soil moisture, and refit the k value. For example, the original model used k=1.2, but the actual measurement shows that the soil texture in this area is heavy and the amplification effect of moisture on trampling damage is more significant. Therefore, k is adjusted to 1.5-1.8. After adjustment, historical data needs to be used for verification to ensure that the corrected model can better fit the observation results of previous years.

[0078] If the trend is reversed, it may be related to a bias in the response of excrement decomposition to meteorological conditions. In this case, the dynamic correction model parameters a and b in step S34 should be adjusted first. The correction method is as follows: collect meteorological data and corresponding decomposition experiment results for multiple grazing cycles, perform regression analysis again, and fit new coefficients a and b. For example, if the original model fits a=0.6 and b=0.4, but data from two consecutive years show that precipitation has a greater impact on the decomposition rate than temperature, then adjust a=0.4 and b=0.6. The constraint condition a+b=1 must be ensured during the correction.

[0079] After parameter correction, backtesting using historical grazing cycle data is required: substitute the corrected parameters into the model, recalculate the nutrient balance index for the historical cycle, and examine whether the correlation with livestock production performance during the same period has improved. If the improvement is significant, the new parameters are used for subsequent grazing cycle calculations; if there is no improvement, other possible sources of bias need to be further investigated, such as remote sensing data inversion errors or insufficient representativeness of soil sampling.

[0080] In the above technical solution, a verification and correction mechanism is introduced to achieve dynamic optimization and adaptive adjustment of the calculation model. The nutrient balance calculation results are correlated with livestock production performance—the ultimate verification indicator—establishing a closed-loop feedback system that allows the model to continuously improve based on actual production feedback. Targeted correction paths are designed for key parameters that may cause deviations (soil animal promotion coefficient, water sensitivity coefficient, and meteorological correction coefficient), avoiding blind adjustments. By determining the opposite trends over two consecutive periods, the influence of random fluctuations is eliminated, ensuring the reliability of correction decisions. This mechanism enables the method to adapt to differences in different years, climatic conditions, and grassland types, significantly improving the method's adaptability, robustness, and long-term application value, ensuring that the calculation results remain highly consistent with actual conditions.

[0081] In another technical solution, when dividing spatial units in step S1, a livestock spatial memory preference factor is further introduced for correction, specifically including: Historical GPS trajectory data from at least one complete grazing cycle was collected, and foraging hotspot distribution maps were extracted using kernel density analysis to generate a livestock spatial memory preference index. M i ; Livestock Spatial Memory Preference Index M i As a correction weight, the comprehensive index of grazing suitability in step S1 is weighted and corrected: ; in I suit,i The index represents the overall grazing suitability index before revision. I suit,i′ The corrected index is λ, which is the preference intensity coefficient, obtained by regression analysis of historical grazing trajectory data and vegetation spatial distribution. Spatial units were divided using a modified comprehensive index of grazing suitability.

[0082] Specifically, the introduction of livestock spatial memory preference factors is based on an important finding in behavioral ecology: livestock develop spatial memory during grazing and tend to repeatedly visit areas with previously favorable foraging conditions, even if the vegetation conditions in these areas have changed. This preference behavior leads to a discrepancy between actual grazing utilization intensity and the suitability index calculated based on static environmental factors (vegetation, topography, water resources). Therefore, by mining livestock spatial preferences through historical trajectory data and correcting the suitability index, spatial unit division can more closely reflect actual grazing patterns.

[0083] Historical GPS trajectory data collection must cover at least one complete grazing cycle to ensure the capture of livestock activity patterns under different seasons and weather conditions. If possible, collecting data from 2-3 cycles will yield better results. The GPS locator's sampling frequency should be set to 1Hz or higher, and the recorded information should include timestamps, latitude and longitude coordinates, instantaneous speed, and positioning accuracy. During data collection, attention should be paid to the device's battery life and data storage capacity; typically, each livestock can generate approximately 8-10MB of trajectory data per day.

[0084] Data preprocessing is fundamental to kernel density analysis. First, data cleaning is performed to remove location drift points (such as outliers exceeding the livestock's maximum speed), duplicate records, and locations with an accuracy less than 10 meters. Then, trajectory segmentation is performed, dividing continuous trajectories by day or by grazing period to facilitate subsequent analysis. For livestock groups with multiple grazing individuals, individual trajectory data need to be merged for analysis to reflect spatial preferences at the group level.

[0085] Kernel density analysis is a core method for generating foraging hotspot distribution maps. Kernel density estimation generates a continuous probability density surface by calculating the density of trajectory points within a certain range around each spatial point, reflecting the spatial clustering of livestock activity. Its mathematical expression is: ; Where (x,y) are the coordinates of the point to be estimated, n is the total number of trajectory points, h is the bandwidth (search radius), K is the kernel function (usually a Gaussian kernel or a fourth kernel), and d i Let be the distance from the point to be estimated to the i-th trajectory point. The choice of bandwidth h has a significant impact on the analysis results: too small a bandwidth will produce many isolated hotspots, while too large a bandwidth will cause the hotspot regions to become overly smoothed. Generally, an empirical formula is used based on the daily activity range of livestock and the area of ​​the study area. The calculation is performed, where σ is the standard deviation and IQR is the interquartile range. In practice, multiple bandwidth values ​​can be tried, and the result that clearly reflects the feeding aggregation characteristics and is consistent with field observations can be selected.

[0086] Kernel density analysis is performed in GIS software (such as ArcGIS's Kernel Density tool and QGIS's Heatmap plugin). The input layer is the preprocessed trajectory points. The cell size (usually consistent with the spatial cell grid size, such as 10m), search radius (determined based on the scale of livestock activity, such as 200-500m), and kernel function type are set. The output is a raster layer, where each cell's value is a kernel density estimate, reflecting the intensity of livestock activity.

[0087] The generation of a feeding hotspot distribution map requires the integration of behavior recognition results. GPS tracks alone cannot distinguish feeding behavior from other behaviors (such as walking and resting), therefore, it needs to be coordinated with the mandibular motion sensor data from step S2. Through time synchronization, behavior labels are added to the track points, and track points in the feeding state are selected for kernel density analysis to generate a feeding hotspot distribution map. If behavior data is unavailable, all track points can be analyzed, but it should be noted that high-density points on walking paths may be misidentified as feeding hotspots.

[0088] Livestock Spatial Memory Preference Index M i The generation process employs normalization. The raster layer obtained from kernel density analysis is then normalized by range normalization: ; KDE i For the kernel density estimate of the i-th spatial cell (or raster pixel), KDE min and KDE max These represent the minimum and maximum kernel density values ​​for the entire study area, respectively. After normalization, M... i The value ranges from 0 to 1, with a larger value indicating a higher preference for that area among livestock. For spatial units divided by a grid, M within each unit is calculated using a partitioning statistical tool. i The average value is used as the preference index for that unit.

[0089] Determining the preference intensity coefficient λ is a crucial step, as it reflects the weight of spatial memory preference in grazing suitability evaluation. The value of λ typically ranges from 0 to 1. When λ=0, the correction is ineffective; when λ=1, the influence of preference factors is comparable to the original suitability index; and when λ>1, preference factors play a dominant role.

[0090] The fitting of λ was performed using regression analysis, based on the relationship between historical grazing trajectory data and vegetation spatial distribution. The specific steps were as follows: The first step is to extract vegetation environmental features from foraging sites along historical grazing routes. For each foraging site, environmental factors such as pasture biomass B, slope S, and distance to water source D are extracted to form a dataset {(B j ,S j D j )}.

[0091] The second step is to construct a proxy variable for grazing utilization intensity. The density of trajectory points within a certain range (e.g., 50m) around each grazing point is used as an approximation of the actual utilization intensity, denoted as U. j .

[0092] The third step is to establish a multiple regression model. This model uses the actual utilization intensity U... j As the dependent variable, environmental factors (B) j ,S j D j ) and preference index M j Using as the independent variable, a linear regression model is adopted: U j =β0+β1B j +β2S j +β3D j +β4M j +ε; Where β1, β2, and β3 are the coefficients of environmental factors, and β4 is the coefficient of the preference index. The preference intensity coefficient λλ can be defined as the ratio of the combined influence of β4 to β1, β2, and β3, or the normalized value of β4 can be directly adopted based on professional judgment.

[0093] In practical applications, if sufficient foraging data is lacking, a simplified method can be used to determine λ. For example, through comparative analysis, the study area can be divided into two categories—a high-density historical grazing area (M>0.7) and a low-density area (M<0.3). The average pasture biomass and grazing suitability index can be calculated for each category, and the ratio between the two can be analyzed to estimate the influence of preference factors. A general empirical value is 0.3-0.5. For pastures with a long history of grazing and stable livestock populations, λ can be set higher; for new pastures or pastures with frequent livestock population changes, λ should be set lower.

[0094] The spatial suitability index is corrected using a weighted formula: I suit,i =I suit,i ×(1+λ×M i ) Where I suit,i This refers to the original comprehensive index of grazing suitability calculated in step S1. The principle of correction is: in the area where livestock prefer grazing (M... iEven if the environmental factors indicate a moderate level of suitability, the actual utilization intensity may still be high, thus requiring an upward adjustment of the suitability index. Conversely, in areas with low preference, even with favorable environmental conditions, livestock may utilize the land less, necessitating maintaining the suitability index or slightly lowering it. This is achieved by multiplying by (1+λM). i The factor enables non-linear adjustment of the original index.

[0095] The modified grazing suitability composite index I was adopted. suit,i The spatial units are re-divided. The specific method is the same as in step S1: the spatial distribution map of pasture biomass is overlaid with the revised comprehensive index of grazing suitability, and the study area is divided into multiple homogeneous spatial units according to the set grid size. The boundaries of the revised units may differ from the original division; areas favored by livestock will be more finely divided, while areas with low preference may be merged.

[0096] The effectiveness of the correction can be verified using cross-validation: historical trajectory data are divided into a training set (e.g., 70%) and a validation set (30%). The training set is used to fit λ and generate the corrected fitness index. Then, the proportion of actual foraging points in the validation set falling within high fitness areas is examined to see if it has significantly increased. If the correction is effective, the density of foraging points in high fitness areas (e.g., the upper 20th percentile) should be significantly higher than before the correction.

[0097] In the aforementioned technical solution, by introducing a livestock spatial memory preference factor, livestock behavioral habits are integrated into spatial unit division, making unit boundaries more consistent with actual grazing patterns. Kernel density analysis using historical GPS trajectory data objectively quantifies livestock spatial preferences, avoiding the arbitrariness of subjective judgment. Regression analysis is used to fit preference intensity coefficients, ensuring the correction degree is data-supported and adjustable. A weighted correction formula is employed, preserving the physical meaning of the original suitability index while incorporating behavioral information. This improvement is particularly suitable for pastures with a long history of grazing, significantly improving the accuracy of nutrient output measurement, making spatial unit division more realistic, and providing a more reliable unit basis for subsequent nutrient balance assessments.

[0098] In another technical solution, when calculating the trampling loss in step S23, the spatial heterogeneity of trampling damage is further considered, specifically including: The GPS track points of livestock are connected to generate a walking path. The width of the buffer zone is dynamically adjusted according to the walking speed of the livestock, with the path as the center line. Inner buffer zone, middle buffer zone and outer buffer zone are established respectively, and different damage weights w1, w2 and w3 are assigned, where w1>w2>w3. Calculate the weighted trampling damage coefficient for each spatial unit based on the proportion of areas falling into different buffer zones: ; in α i The coefficient of damage caused by trampling. A ik Let be the area within the i-th spatial cell that falls into the k-th level buffer. A i Let i be the total area of ​​the i-th spatial unit. w k Let be the damage weight of the k-th level buffer; Weighted trampling damage coefficient α i ′ Calculate the amount of damage caused by trampling.

[0099] Specifically, traditional trampling damage calculations typically employ the assumption of uniform damage, assuming that the trampling effect is the same at all locations within a cell. However, actual observations show that trampling damage is most severe along the livestock's walking path, with the degree of damage decreasing with increasing distance on both sides of the path. Therefore, this technical solution significantly improves calculation accuracy by refining the spatial distribution of trampling damage.

[0100] The walking path is generated based on continuous trajectory points recorded by GPS trackers worn by livestock. First, the raw GPS trajectory data is preprocessed to remove drifting points (those with speeds greater than the livestock's maximum speed, e.g., cattle > 2.5 m / s, sheep > 2.0 m / s) and stationary points (those with speeds < 0.1 m / s and a duration exceeding 5 minutes). Then, the cleaned trajectory points are connected in chronological order to generate a continuous linear path. For adjacent points with time intervals exceeding 30 seconds, path interruptions are considered, and they are not connected. Path generation can be achieved using the "point set to line" tool in GIS software.

[0101] The dynamic adjustment of the buffer zone width is based on the walking speed of livestock. Studies have shown that the faster the livestock walks, the larger their stride, the stronger the impact force, and the more significant the crushing and tearing effects on vegetation, as well as the wider the damage area. Therefore, this scheme introduces a speed-dependent buffer zone width setting method. For each trajectory segment (the line segment between two adjacent trajectory points), the buffer zone width is calculated based on the average walking speed v (unit: m / s) within that segment: Inner buffer half-width = a × v + b; The coefficients a and b were determined through field observations. For example, observational data on cattle showed that when the walking speed v = 0.5 m / s, the half-width of the inner buffer zone was approximately 0.5 m; when v = 1.5 m / s, the half-width of the inner buffer zone was approximately 1.0 m. The fitted values ​​yielded a = 0.5 and b = 0.25. The widths of the middle and outer buffer zones can be proportionally increased based on the width of the inner buffer zone; for example, the half-width of the middle buffer zone could be twice that of the inner buffer zone, and the half-width of the outer buffer zone could be three times that of the inner buffer zone. While fixed width gradations can be used to simplify calculations, dynamic speed adjustment better reflects the actual situation.

[0102] The weights assigned to buffer zone damage levels should reflect the attenuation of damage with distance. The inner buffer zone (the innermost area on either side of the path centerline) is where livestock hooves directly trample, resulting in the most severe vegetation damage, and is assigned the highest weight w1. The middle buffer zone (the area outside the inner buffer zone) is mainly affected by mud splashed up by hooves and indirect crushing, resulting in moderate damage, and is assigned a weight w2. The outer buffer zone (the outermost area) is mainly affected by friction from livestock bodies and sporadic coverage by feces, resulting in less damage, and is assigned a weight w3. The weight values ​​must satisfy w1 > w2 > w3 and be normalized (i.e., the sum of the three weights does not necessarily equal 1, but their relative magnitudes reflect the proportion of damage). Based on field experiments, the recommended value ranges are: w1 = 0.8–1.0, w2 = 0.3–0.5, w3 = 0.1–0.2. For example, in a calibration experiment in a typical grassland area, the vegetation biomass in the inner buffer zone decreased by 80%, the middle buffer zone by 30%, and the outer buffer zone by 10%. Therefore, w1=0.8, w2=0.3, and w3=0.1 were selected.

[0103] The area ratio calculation is performed in GIS. The generated buffer layer (containing inner, middle, and outer buffer surface features) is overlaid and analyzed with the spatial cell grid layer. Using the "Intersection" or "Identifier" tool, the area A within each spatial cell that falls into each level of buffer is calculated. i1 A i2 A i3 Simultaneously obtain the total area A of this spatial unit. i (i.e., grid area). Calculate the area ratio p of each level of buffer zone. ik =A ik / A i .

[0104] The weighted trampling injury coefficient is calculated using the formula described above. Note that α here... i This is the basic damage coefficient obtained in step S23 after soil moisture correction. It reflects the proportion of vegetation damage caused per unit step frequency under the assumption of uniform trampling. After introducing spatial heterogeneity correction, the actual damage coefficient α i ′ equals the base coefficient multiplied by the area-weighted average damage weight of each buffer level. That is: ; The physical meaning of this formula is: for the unit damage coefficient αi under the assumption of uniform trampling, when only some areas within a spatial unit are subjected to trampling of varying intensities, the actual average damage level should be adjusted through a weighted average. For example, if the inner buffer zone accounts for 20%, the middle buffer zone 30%, and the outer buffer zone 50% within a unit, with weights of 0.8, 0.3, and 0.1 respectively, then the weighted sum is 0.2 × 0.8 + 0.3 × 0.3 + 0.5 × 0.1 = 0.3, and the actual damage coefficient αi is... i =0.3α i The coefficient is lower than the original coefficient because most of the area is in the less damaged outer buffer zone.

[0105] Finally, the corrected damage coefficient α is used. i Recalculate trampling losses: ; This calculation will more accurately reflect the actual damage to vegetation caused by livestock walking paths, avoiding overestimation or underestimation due to the uniformity assumption.

[0106] To further verify the correction effect, fixed monitoring quadrats can be set up within the sample plots to compare the differences between the calculated values ​​under the uniformity assumption, the weighted corrected calculated values, and the measured vegetation damage (such as the percentage reduction in biomass). Measured data can be obtained by comparing the mowing of quadrats before and after grazing. If the weighted corrected calculated values ​​match the measured values ​​better (error <10%), then the method is proven to be effective.

[0107] In the aforementioned technical solution, a refined characterization of the spatial heterogeneity of trampling damage is achieved for the first time by introducing multi-level buffer weighted correction. Walking paths are generated based on GPS trajectories, and the buffer width is dynamically adjusted according to walking speed, making the identification of damaged areas more realistic. By assigning different damage weights to each level of buffer, the natural law of damage attenuation with distance is reflected. Area ratio weighted calculation makes the average damage coefficient within spatial units more accurate. This improvement solves the systematic bias caused by the traditional method treating trampling damage as a homogeneous distribution, significantly improving the calculation accuracy of trampling loss, and thus providing a more reliable basis for the accurate measurement of nutrient output. At the same time, this method fully utilizes the potential of existing GPS trajectory data, requires no additional hardware investment, and has the advantages of low cost and easy promotion.

[0108] 1. Verify excretion parameters (refer to Table 1): When calculating the nutrient content of excrement in step S3, the excretion data of "sheep / goats" in Table 1 (1.5-3.5 kg / d) is introduced to verify whether the total excretion of the whole herd calculated by GPS and thermal imaging methods is within a reasonable range, thereby proving the accuracy of the method of the present invention.

[0109] Table 1. Emissions and excrement decomposition rates of grazing livestock (based on International Livestock Units) livestock species Fresh fecal excretion (kg / d) Urine excretion (kg / d) Total excrement (fresh weight kg / d) dairy cow 30-50 15-25 45-75 beef cattle 20-35 10-18 30-53 Sheep / Goats 1.5-3.5 0.6-1.5 2.1-5 horse 15-25 5-10 20-35 2. Verify the decomposition rate (refer to Table 2): After obtaining the baseline conversion coefficient through the decomposition bag experiment in step S3, the decomposition rate range of "sheep manure" in Table 2 (0.3-2 years) is introduced to verify whether the decomposition rate (or conversion coefficient) obtained from the experiment is consistent with common knowledge in the field, and to further enhance the reliability of the experimental results.

[0110] Table 2. Decomposition rate range of various organic materials in natural grasslands Organic types Decomposition rate (70% decomposition) Key features and descriptions Alfalfa (legume) 0.3-1 year It decomposes fastest. Its high nitrogen content and low C / N ratio provide an abundant nitrogen source for microorganisms. Natural grassland litter 0.5-2 years It adapts to local microorganisms, has a moderate C / N ratio, and its decomposition rate varies greatly along the precipitation gradient. sheep dung 0.3-2 years The particles are small, resulting in a large contact area. They decompose quickly under humid conditions but harden more slowly under dry conditions. cow dung 0.6-2 years Decomposition is highly dependent on dung beetle activity. Decomposition is significantly accelerated in ecosystems with dung beetles. horse manure 0.6-3 years It has a high fiber content, a high C / N ratio, and a coarse texture, making it one of the slowest-decomposing livestock manures. oat straw 0.8-2 years It is softer than wheat and corn stalks, has a lower C / N ratio, and decomposes slightly faster. wheat straw 1-2.5 years It has a high degree of lignification, a high C / N ratio, and decomposes slowly on grasslands. corn stalks 1.5-3 years The slowest decomposing straw. The stalks are thick and hard, with extremely high lignin content, resulting in very slow initial decomposition. Camel dung 1-3 years The slowest decomposing type of feces. Extremely dry, in the form of hard cakes, with a very high C / N ratio, it mainly decomposes through physical weathering. urine Hours to days It is not a solid substance; it quickly seeps into the soil and is absorbed by plants or volatilized or leached away. The technical solution of the present invention will be further described in detail below with reference to specific application examples, so that those skilled in the art can implement it accordingly. This embodiment takes the Tan sheep grazing system in the typical grassland area of ​​Ningxia, my country as an example to fully illustrate the method of the present invention, but the scope of protection of the present invention is not limited to this specific example.

[0111] Implementation scenario: A natural grassland pasture in the arid central region of Ningxia, with a total area of ​​approximately 1,200 hectares, was selected. The grazing livestock were Tan sheep, totaling 500 head, for a complete grazing cycle (June to October of the same year, a total of 150 days).

[0112] S1. Couple multi-source data and divide into refined spatial units. This step aims to address the spatial heterogeneity of natural grasslands by dividing the study area into relatively homogeneous spatial units.

[0113] Step 1: Data Acquisition and Biomass Inversion. Before the start of the grazing cycle (end of May), an aerial survey of the study area was conducted using a DJI Phantom 4 multispectral drone. The flight altitude was set at 100 meters, with 85% forward overlap and 75% lateral overlap, acquiring multispectral images with a spatial resolution of 5 cm. Through preprocessing and stitching, orthophotos and NDVI layers of the study area were generated. Simultaneously, 30 1m × 1m ground plots were established in the study area, and the center point coordinates were recorded. Forage was cut at ground level, dried at 65℃, and weighed to obtain the aboveground biomass dry weight (g / m³). 2 A linear regression model was established between NDVI and measured biomass: Biomass = 375 × NDVI (R0) 2 =0.88). The model was applied to the NDVI layer to generate a spatial distribution map of pasture biomass in the study area, in kg / ha.

[0114] Step 2: Constructing a comprehensive grazing suitability index. Slope data was extracted based on a 12.5m resolution digital elevation model (DEM) of the study area. The distances to the five existing water sources were calculated using the Euclidean distance tool in GIS. Considering the ecological habits of the Tan sheep (dry-loving, hardy, and moderately dependent on water sources), the weights were determined using the analytic hierarchy process (AHP) to construct the comprehensive grazing suitability index. I suit =0.5×(Biomass normalized value)+0.2×(1-slope normalized value)+0.3×(1-water source distance normalized value). All factors are normalized to the 0-1 range.

[0115] Step 3: Spatial Unit Division. Taking into account the topography and vegetation heterogeneity of the study area, a 50m × 50m grid was established, dividing the entire study area into approximately 4800 spatial units. Through GIS zoning statistics, an average pasture biomass B was assigned to each unit. i and average grazing suitability index I suit,i .

[0116] Step 4 (Optimal Solution): Introducing a livestock spatial memory preference factor for correction. GPS trajectory data of 500 Tan sheep from the same period of the previous year were collected. Kernel density analysis was used to generate a grazing hotspot distribution map of the Tan sheep. After normalization, the livestock spatial memory preference index M for each unit was obtained. i Regression analysis with historical grazing point data yielded a preference intensity coefficient λ = 0.4. The grazing suitability index was then adjusted as follows: I suit,i '=I suit,i ×(1+0.4×M i ).

[0117] Using the modified exponent I suit,i 'As the basis for the final spatial unit division and evaluation.'

[0118] S2. Precisely quantify the dual losses from "feeding + trampling" and calculate the net nutrient output. This step aims to accurately calculate the total net nutrients that livestock ingest from each spatial unit.

[0119] Step 1: Livestock Behavior Monitoring. Thirty representative Tan sheep (6% of the flock) were fitted with collars integrating GPS locators and mandibular motion sensors. The GPS sampling frequency was 1Hz, and the mandibular sensor sampling frequency was 25Hz. Using machine learning algorithms, based on the characteristics of mandibular motion waveforms, behavioral states such as feeding, walking, and resting were accurately identified. During the grazing cycle, the spatiotemporal trajectory and behavioral data of each sheep were continuously recorded. Finally, the total feeding time t of all monitored sheep within each spatial unit i was calculated. ij and average walking cadence S step,ij .

[0120] Step 2: Establish a forage nutrient database. Samples of major forage species (Stipa benthamiana, Artemisia frigida, Artemisia scoparia, etc.) were collected in the study area, and their nutrient contents, including crude protein, neutral detergent fiber, and acid detergent fiber, were determined. Using DNA macrobarcoding technology to analyze Tan sheep manure samples, the diet composition was determined to be: 60% Stipa benthamiana, 25% Artemisia frigida, and 15% other miscellaneous grasses. Based on the dietary composition ratio, the weighted average crude protein content (c) of the forage consumed by Tan sheep was calculated to be 10.5%.

[0121] Step 3: Calculate the losses from foraging and trampling to obtain the net output.

[0122] Calculation of feed loss: F i =Σ j (t ij ×r j ×c) Among them, the average feeding rate r of Tan sheep was measured by a jaw sensor. j The flow rate is 55 mouthfuls / minute. Taking nitrogen (N) as the calculation object (crude protein nitrogen content 16%), then c(N) = 10.5% × 16% = 1.68%. Assume t... ij This represents the total feeding time of all monitored sheep within a certain unit. By extrapolating the model to the entire flock (500 sheep), the nitrogen output F caused by feeding in that unit is finally calculated. i .

[0123] Calculation of trampling loss: T i =Σ j S step,ij ×L step ×W trample ×B i ×α i ' Among them, the average stride length of Tan sheep is L step Take 0.5m as the effective trampling width W trample Take 0.3m. B i Provided by step S1. The basic damage coefficient α0 was determined to be 0.15 through fencing experiments. The average soil volumetric water content θ of each unit during grazing was obtained through remote sensing inversion. i Given a moisture sensitivity coefficient k of 1.0 and a base trampling damage coefficient α of 15%, the base trampling damage coefficient is... i =0.15×(1+1.0×15%)=0.1725.

[0124] Trampling Space Heterogeneity Correction (Preferred Solution): GPS track points are used to generate walking paths, and buffer zones are dynamically established. Three levels of buffer zones are defined based on the walking speed of the Tan sheep and assigned weights: inner buffer zone (damage weight w1=0.8), middle buffer zone (w2=0.4), and outer buffer zone (w3=0.1). Calculations show that if the area ratios of the inner, middle, and outer buffer zones within a given unit are 10%, 20%, and 70%, respectively, then the weighted trampling damage coefficient is: α i =0.1725×(0.1×0.8 + 0.2×0.4 + 0.7×0.1)=0.1725×0.23=0.0397. Use this corrected α i 'Recalculate T' i The result was about 77% lower than before the correction, which is more in line with the actual situation that the activities of Tan sheep are scattered and the impact of trampling is limited.

[0125] Net Output: Output i =F i -T i Summing over the entire grazing cycle and all units in the entire study area yielded a total nutrient (nitrogen) output of approximately 9800 kg.

[0126] S3. Dynamically simulate the entire process of "excretion-decomposition-transformation" to calculate the effective nutrient return. This step aims to accurately calculate the effective nutrients in livestock excrement that actually enter the soil-plant cycle.

[0127] Step 1: Constructing a dynamic spatiotemporal layer of excrement deposition. Excretion events were automatically identified using jaw movement sensors monitoring sheep, recording the time and location. During the grazing cycle, a UAV thermal imaging flight (in the early morning) was conducted weekly to identify fresh excrement patches. These thermal imaging patches were then spatiotemporally matched with GPS excrement events. The result was a dynamic layer containing approximately 8500 excrement events (patches), recording the deposition time, location coordinates, and area (average area approximately 0.15m²) for each patch. 2 The spatial unit to which they belong is also considered. Based on this, it is estimated that the total daily excrement of 500 Tan sheep is approximately (8500 patches / 150 days) × 0.15m. 2 / patch × 0.6kg / m 2 The result is approximately 5.1 kg / d (fresh weight), which is basically consistent with the range of sheep fresh feces excretion (1.5-3.5 kg / d) listed in Table 1. Considering the actual moisture content and the estimation error of the patches, the calculation results of this invention have high reliability.

[0128] Step 2: Determine basic parameters of excrement. Collect fresh sheep feces samples and determine their dry matter content (40%) and dry matter weight per unit area (D). dung Approximately 0.6 kg / m 2The total nitrogen content of the fecal sample was determined to be 2.0% (dry matter basis). The nitrogen content C per unit area of ​​excrement was calculated. dung =0.6kg / m 2 ×2.0%=12g / m 2 .

[0129] Step 3: Establishing the baseline conversion coefficient. In-situ decomposition bag experiments were conducted in the study area. Nylon mesh bags (2mm mesh) filled with fresh sheep manure were placed in the plots, and samples were taken on days 15, 30, 60, 90, and 150 to determine the total nitrogen content of the residues. Simultaneously, nitrogen and oxygen stable isotope analyses were performed. Using an isotope mass balance model, the baseline ratio R0 of excrement nitrogen actually converted into available soil nutrients was calculated at different time points. For example, at the end of a 150-day grazing cycle, R0 was 0.45. Based on the decomposition experiment data, it was estimated that the time required for 70% of the organic matter in sheep manure to decompose was approximately 0.8-1.5 years. This is highly consistent with the sheep manure decomposition rate range (0.3-2 years) listed in Table 2, verifying the representativeness of the experimental results.

[0130] Step 4: Dynamically adjust the decomposition and transformation coefficients.

[0131] Meteorological factor correction: Small weather stations were set up at the decomposition experiment sites to record daily temperature and precipitation. The cumulative accumulated temperature T during the grazing cycle was calculated. sum The effective precipitation is 3200℃·d. eff It is 180mm. According to local meteorological data, the multi-year average accumulated temperature T... sum The average annual effective precipitation is 3000℃·d. eff The value is 160mm. Regression analysis fitted empirical coefficients a=0.58, b=0.42. Therefore, the meteorologically corrected conversion coefficient is: R0 weather =0.45×(0.58×(3200 / 3000)+0.42×(180 / 160))=0.45×(0.62+0.47)=0.49.

[0132] Soil animal promotion coefficient correction: The survey found that the density of soil animals such as dung beetles was high in the study area. The promotion coefficient β was fitted to 1.2 through indoor controlled experiments.

[0133] Final effective conversion factor: R = R weather ×β = 0.49 × 1.2 = 0.588.

[0134] Step 5: Calculate the effective nutrient return. For each excrement patch, based on its deposition time, apply the corresponding dynamic conversion coefficient `R` (for patches with earlier deposition times, the decomposition time is longer, and the conversion coefficient can be slightly higher than the average; for simplicity, the average coefficient of 0.588 over the entire period is used here) to calculate its effective nutrient return, and sum them to obtain the effective return for each spatial unit: Supply deposition,i =Σ k [A patch,ik ×D dung ×C dung ×R].

[0135] Summing over the entire study area, the total effective nutrients (nitrogen) returned by excrement during the entire grazing cycle is approximately 5900 kg.

[0136] S4. Spatial overlay analysis to generate precise control decisions. This step aims to assess nutrient balance by taking into account soil baseline supply and to develop an actionable grazing management plan.

[0137] Step 1: Calculate the nutrient balance index for each spatial unit. After the grazing cycle ends, soil samples from 0-20cm depth are collected from each spatial unit to determine the available nitrogen content. Combined with soil bulk density (1.3g / cm³), the available soil nutrients are calculated. soil,i The average nitrogen input from dry and wet atmospheric deposition across the entire study area is approximately 65 kg / ha. atmo,i Referencing data from nearby weather stations, a uniform nutrient balance of 5 kg / ha was used. The nutrient balance index for each unit was calculated: B i =Output i / (Supply soil,i +Supply deposition,i ).

[0138] Step 2: Identify nutrient surplus and deficit zones. Set a threshold: B i ≤0.8 indicates a nutrient surplus zone; 0.8 i <1.2 is the equilibrium zone; B i A value ≥1.2 indicates a nutrient deficit area. GIS analysis results show: Nutrient surplus areas are mainly distributed within 500 meters of the water source, accounting for 18% of the total study area. Nutrient balanced areas account for 65% of the total study area. Nutrient deficit areas are mainly distributed on slope tops and gentle slopes far from the water source, accounting for 17% of the total study area.

[0139] Step 3: Generate the regulation and control decision map.

[0140] ​Resting grazing area map: For nutrient-deficient areas, the remaining biomass of pasture (areas below 300 kg / ha) and topographic slope (areas greater than 15°) obtained by remote sensing inversion at the end of grazing are combined to designate the overlapping areas that meet the conditions as resting grazing areas at the beginning of the next grazing cycle. The vector map of the resting grazing area and the resting grazing suggestion (such as resting grazing until the pasture turns green again the following year) are output.

[0141] Priority grazing area map: For nutrient surplus areas, based on the distance to water source (less than 1.5km) and slope (less than 10°), suitable areas for priority utilization are selected to generate a priority grazing area map for the next grazing cycle.

[0142] Recommended location map for supplementary feeding points: For areas with persistent deficiencies and located on livestock migration routes, it is recommended to set up 3 supplementary feeding points, and estimate the recommended amount of supplementary feed for each feeding point during the dry season based on the degree of deficiency (e.g., 0.2 kg of concentrate per sheep per day).

[0143] Verification and Correction: After the grazing cycle ends, a correlation analysis is performed between the average nutrient balance index of the entire farm calculated by this method and the average daily weight gain of Tan sheep during the same period (measured at 0.15 kg / day). The results are consistent with historical experience trends, indicating that the model parameters are reliable and no correction is needed. If the trend is opposite, the verification and correction method will be used after each grazing cycle to correct k, a, b, or the soil animal promotion coefficient.

[0144] In summary, this embodiment, through the implementation of the method of the present invention, successfully achieved accurate calculation of the nutrient balance of the typical grassland sheep grazing system in Ningxia, and generated specific suggestions for rest grazing, rotational grazing, and supplementary feeding, providing a solid scientific basis for the sustainable utilization and refined management of grasslands.

[0145] The number of devices and processing scale described herein are for the purpose of simplifying the description of the invention. Applications, modifications, and variations of the invention will be readily apparent to those skilled in the art.

[0146] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details.

Claims

1. A method for calculating the nutrient balance of a natural grassland grazing system, characterized in that, Includes the following steps: S1. Acquire multi-source remote sensing data and grazing trajectory data of the target area, and divide the natural grassland into multiple spatial units by coupling the vegetation spatial heterogeneity index and the livestock utilization intensity index. S2. By identifying feeding behavior through GPS locators and jaw movement sensors worn by livestock, combined with the diet composition determined by DNA macrobarcode sequencing, and deducting the amount of mechanical damage to vegetation caused by trampling and resting, the total net nutrient intake of livestock in each spatial unit is calculated as the nutrient output. S3. Based on the location and time of excretion events recorded by GPS locators and mandibular motion sensors, and combined with the spatial distribution of excrement patches identified by UAV thermal imaging, a spatiotemporal dynamic layer of excrement deposition is constructed. Representative excrement samples were collected for decomposition bag experiments, and the changes in nitrogen and oxygen stable isotopes were measured simultaneously to determine the baseline conversion coefficient of excrement nutrients to soil available nutrients during the grazing cycle. At the same time, the population density of soil animals was investigated, and the promoting coefficient of soil animal activity on the decomposition rate of excrement was fitted through indoor controlled experiments. The baseline conversion coefficient is input into the decomposition model coupled with micrometeorological data, and the decomposition rate in the model is dynamically corrected using the soil animal promotion coefficient. Combined with the spatiotemporal dynamic layer of excrement deposition, the effective nutrient return amount that actually enters the soil-plant cycle in each spatial unit during the grazing cycle is calculated. S4. Perform spatial overlay analysis on the nutrient output of S2 and the effective nutrient return of S3, calculate the nutrient balance index of each spatial unit, identify nutrient surplus or deficit areas, and generate control decisions for rest grazing, rotational grazing, or supplementary feeding.

2. The method for calculating the nutrient balance of a natural grassland grazing system as described in claim 1, characterized in that, Step S1 divides the natural grassland into multiple spatial units, specifically including: The Normalized Differential Vegetation Index (NDVI) of forage grass was obtained by using multispectral data from UAVs or satellite remote sensing data. A regression model between NDVI and forage grass biomass was established by combining field quadrat surveys, and a spatial distribution map of forage grass biomass was generated. Slope data was extracted based on the digital elevation model (DEM), and a comprehensive grazing suitability index was constructed by combining the distance to water source and the type of livestock. Using GIS spatial analysis tools, the spatial distribution map of pasture biomass was overlaid with the comprehensive index of grazing suitability for analysis, and the study area was divided into multiple homogeneous spatial units according to the set grid size.

3. The method for calculating the nutrient balance of a natural grassland grazing system as described in claim 1, characterized in that, Step S2 calculates the total net nutrient intake of livestock in each spatial unit, specifically including: S21. Synchronize the data from the GPS locator worn by livestock with the data from the mandibular motion sensor, identify feeding, walking, and resting behaviors through the characteristics of the mandibular motion waveform, and count the feeding time and walking frequency of each spatial unit. S22. Set up representative quadrats in the study area, collect samples of major forage grasses to determine the nutrient content of crude protein, neutral detergent fiber and other nutrients, and establish a forage grass nutrient database. S23. Differentiate between foraging loss and trampling loss, and calculate the net intake separately: Feed loss ,in Let be the feeding time of the j-th livestock in unit i. Let j be the number of feed intakes per unit time for the j-th livestock. The average nutrient content of pasture grass consumed by livestock; Trampling damage ,in S step,ij Let j be the walking cadence of the j-th livestock in the i-th unit. L step For the average step size, W trample To trample on the effective width, B i Let i be the pasture biomass of the i-th unit. α i The trampling damage coefficient is dynamically corrected using synchronously monitored soil moisture content. α i = α 0×(1+ k × θ i ), α 0 is the base damage coefficient. θ i Let be the soil volumetric water content of the i-th unit. k Moisture sensitivity coefficient; Net nutrient output i = F i - T i .

4. The method for calculating the nutrient balance of a natural grassland grazing system as described in claim 1, characterized in that, In step S3, when calculating the effective nutrient return, the decomposition and conversion coefficient of excrement is dynamically adjusted by combining meteorological data during the grazing cycle. Specifically, this includes: S31. Based on the excretion events recorded by the GPS locator and the mandibular motion sensor, generate the time series and location information of the excretion events in each spatial unit, and combine the area of ​​excrement patches identified by UAV thermal imaging to establish a spatiotemporal dynamic layer of excrement deposition. S32. Collect representative fresh excrement samples, determine the content of total nitrogen, total phosphorus, total potassium, and organic matter, and calculate the nutrient content per unit area of ​​excrement. C dung ; S33. Establish a baseline conversion coefficient through a decomposition bag experiment over at least one grazing cycle. R 0 Nylon biodegradable bags containing a known mass of excrement were placed in representative plots. Residual nutrient content was measured at multiple time points, and nitrogen and oxygen stable isotope analyses were performed simultaneously. The proportion of nutrients actually entering the soil-plant cycle was calculated based on isotope mass balance and used as the baseline conversion coefficient. R 0 ; S34. Simultaneously record daily meteorological data during the decomposition experiment, including accumulated temperature. T sum and effective precipitation P eff Based on the historical multi-year average meteorological data of the region, a dynamic correction model for the decomposition transformation coefficient and meteorological factors was established: ； in , This refers to the multi-year average accumulated temperature and effective precipitation. a , b These are the empirical coefficients fitted through regression analysis; S35. Based on the spatiotemporal dynamic layer of excrement deposition, for each excretion event, apply the corresponding dynamic decomposition and transformation coefficient R according to its deposition time, and sum to obtain the effective nutrient return of each spatial unit: ； in A patch,ik Let be the area of ​​the k-th excrement patch in the i-th spatial unit. D dung It refers to the dry matter weight of excrement per unit area.

5. The method for calculating the nutrient balance of a natural grassland grazing system as described in claim 1, characterized in that, Step S4 generates a decision on rest grazing, rotational grazing, or supplemental feeding, specifically including: S41. Calculate the nutrient balance index of each spatial unit according to the following formula: B i =Output i / (Supply soil,i +Supply deposition,i ); Output i The nutrient output of the i-th spatial unit is calculated in step S2; Supply soil,i For the available nutrients in the i-th spatial unit, the available nutrient content is determined by collecting soil samples and calculating the available nutrient content in combination with soil bulk density; Supply deposition,i The total amount of this nutrient input to atmospheric dry and wet deposition was obtained by actual measurement using deposition collectors deployed in the study area or by referencing data from nearby meteorological stations. S42. Set the nutrient balance index threshold: when B i When the range is ∈ [0.8, 1.2], it is determined to be a nutrient balance zone, B i >1.2 indicates a nutrient deficiency zone, B i <0.8 indicates a nutrient surplus zone; S43. For nutrient-deficient areas, generate a grazing rest area map by combining the remaining biomass of pasture based on remote sensing inversion and topographic conditions; for nutrient-surplus areas, generate a priority grazing area map for the next grazing cycle by combining water source distribution and slope conditions; for areas with continuous nutrient deficiency that livestock must pass through, generate a map of recommended locations for supplementary feeding points.

6. The method for calculating the nutrient balance of a natural grassland grazing system as described in claim 1, characterized in that, A review and correction are conducted at the end of each grazing cycle: A correlation analysis is performed between the current nutrient balance calculation results and the livestock production performance data of the same period. If the calculation results of two consecutive grazing cycles are contrary to the actual production performance trend, then the soil animal promotion coefficient in the decomposition model described in step S3, or the water sensitivity coefficient in step S23, will be adjusted. k Or the dynamic correction of model parameters in step S34 a , b Make corrections.

7. The method for calculating the nutrient balance of a natural grassland grazing system as described in claim 2, characterized in that, In step S1, when dividing the spatial units, a livestock spatial memory preference factor is further introduced for correction, specifically including: Historical GPS trajectory data from at least one complete grazing cycle was collected, and foraging hotspot distribution maps were extracted using kernel density analysis to generate a livestock spatial memory preference index. M i ; Livestock Spatial Memory Preference Index M i As a correction weight, the comprehensive index of grazing suitability in step S1 is weighted and corrected: ； in I suit,i The index represents the overall grazing suitability index before revision. I suit,i ′ The corrected index is λ, which is the preference intensity coefficient, obtained by regression analysis of historical grazing trajectory data and vegetation spatial distribution. Spatial units were divided using a modified comprehensive index of grazing suitability.

8. The method for calculating the nutrient balance of a natural grassland grazing system as described in claim 3, characterized in that, When calculating the trampling loss in step S23, the spatial heterogeneity of trampling damage is further considered. Specifically, this includes: connecting the GPS trajectory points of livestock to generate a walking path, using the path as the center line, dynamically adjusting the buffer width according to the walking speed of the livestock, establishing an inner buffer, a middle buffer, and an outer buffer respectively, and assigning different damage weights w1, w2, and w3, where w1>w2>w3. Calculate the weighted trampling damage coefficient for each spatial unit based on the proportion of areas falling into different buffer zones: ； in α i The coefficient of damage caused by trampling. A ik Let be the area within the i-th spatial cell that falls into the k-th level buffer. A i Let i be the total area of ​​the i-th spatial unit. w k Let be the damage weight of the k-th level buffer; Weighted trampling damage coefficient α i ′ Calculate the amount of damage caused by trampling.

Images