Method for identifying patches of grassland sustainability decline anomaly based on time series grid

Abstract

The application discloses a kind of grassland sustainability decline abnormal patch identification method based on timing grid, belong to image data processing field.To solve the problem that single index representation is not complete in existing method, fixed threshold is poor in adaptability, abnormal result is broken, it is difficult to accurately identify the abnormal area under the condition of relative historical state decline, the grassland sustainability index of fusion multidimensional index is constructed, to detect year before continuous years index is constructed year by year historical baseline, the difference of grassland sustainability decline is calculated, and maximum interclass variance method is used to determine annual adaptive threshold, after spatial filtering and patching processing, abnormal patch and attribute data are generated.The application can be used for grassland monitoring and degradation positioning monitoring in ecological environment industry, and can effectively improve the adaptability and spatial integrity of abnormal patch identification.

Description

Technical Field

[0001] This invention relates to the field of image data processing, and more specifically to a method for identifying abnormal patches of grassland sustainability degradation based on temporal raster. Background Technology

[0002] Grassland ecosystems are an important component of terrestrial ecosystems, serving multiple functions including ecological barrier, livestock production, and biodiversity maintenance. Their operational status is jointly regulated by climate fluctuations, vegetation productivity, grazing intensity, and grassland carrying capacity. Due to significant spatial heterogeneity in hydrothermal conditions, grassland types, and utilization patterns across different regions, the same region can exhibit marked state fluctuations in different years due to factors such as interannual climate fluctuations and changes in human utilization intensity. Grassland degradation processes can exhibit both long-term evolutionary trends and localized, sudden, phased declines. Therefore, accurately identifying abnormal areas of degradation that deviate from their historical normal state from monitoring data spanning multiple years is a key technical problem that urgently needs to be solved in the field of grassland ecological remote sensing monitoring.

[0003] Currently, grassland ecological monitoring primarily uses single vegetation indicators such as Normalized Difference Vegetation Index (NDVI), vegetation cover, aboveground biomass, or net primary productivity (NPMP) to characterize grassland growth. While these indicators can directly reflect the growth level of grassland vegetation, their dynamic changes are not entirely equivalent to changes in the sustainable use status of grassland. On the one hand, in areas where grassland productivity remains at a high level, if grazing pressure consistently approaches or even exceeds the grassland's carrying capacity, the apparent vegetation growth may remain relatively stable, but the risk of overexploitation of grassland resources has been accumulating continuously, and a single vegetation indicator cannot identify such hidden degradation risks. On the other hand, a short-term decline in grassland productivity may only be caused by abnormal climate factors such as annual precipitation and temperature, and does not necessarily represent a continuous deterioration in the sustainable status of grassland. Relying solely on a single vegetation indicator cannot simultaneously cover multiple dimensions such as productivity, utilization pressure, interannual stability, and long-term overgrazing risk, and therefore cannot comprehensively depict the true state of grassland sustainability.

[0004] To improve the comprehensiveness of grassland status assessment, existing technologies typically standardize multiple ecological indicators and then weight them to synthesize a comprehensive evaluation index. Grassland status levels are then classified based on preset fixed thresholds. While this method can reflect spatial differences in grassland status across different regions, its evaluation results tend to be static and do not fully consider the differences in background conditions among different spatial pixels in the image. Different geographic units have different productivity bases and grazing use backgrounds. A uniform fixed threshold is difficult to adapt to the historical change characteristics of each pixel. For grasslands that have been at a low productivity level for a long time but whose status is relatively stable, a fixed threshold easily leads to their continued classification as low-value degraded areas. Conversely, for grasslands that were originally in good condition but have recently experienced rapid decline, a fixed threshold cannot promptly identify their relative degradation anomalies. Classification based solely on absolute values ​​cannot accurately locate abnormal areas that have significantly declined relative to their historical status.

[0005] In terms of temporal change detection, existing grassland dynamic analysis methods mostly use the adjacent year difference method, long-term trend analysis method, or unified reference period mean method to determine changes. Among them, the adjacent year difference method is easily affected by single-year climate anomalies and data noise, resulting in a high false detection rate. The long-term trend analysis method is more suitable for depicting the overall evolution direction over multiple years, but it is insufficient for responding to local sudden declines in a single year. Although the unified reference period mean method can provide a unified comparison benchmark, it does not fully consider the differences in the number of effective observations and historical value distribution of each pixel, making it difficult to take into account the heterogeneity of the historical state of different spatial pixels and unable to conduct relative anomaly determination for the historical state of each pixel.

[0006] Furthermore, anomaly results extracted directly based on numerical thresholds are typically composed of discrete, independent pixels. These are easily contaminated with data noise, isolated anomaly pixels formed by local data gaps and boundary disturbances. If such discrete pixels are directly output as anomaly regions, the results often exhibit highly fragmented characteristics, making it difficult to form monitoring objects with clear spatial boundaries and thus unable to directly support grassland management applications. Simply using connected component labeling methods also retains a large number of fragmented areas composed of a small number of pixels, which similarly lack practical management significance and are difficult to use as effective targets for grassland degradation monitoring and control.

[0007] In summary, existing technologies are insufficient to simultaneously achieve multi-dimensional comprehensive characterization of grassland sustainability, differentiated construction of pixel-by-pixel historical baselines, adaptive determination of annual anomaly thresholds, and spatial patch aggregation of discrete anomaly pixels. This fails to meet the operational needs for accurate identification and spatial monitoring of grassland sustainability degradation anomalies. Therefore, a new identification method is urgently needed to overcome these shortcomings. Summary of the Invention

[0008] To alleviate or partially alleviate the above-mentioned technical problems, the solution of the present invention is as follows:

[0009] A method for identifying anomalous patches of grassland sustainability decline based on temporal raster data includes:

[0010] The aboveground biomass raster data and grazing intensity raster data of the area to be detected over many years are obtained. After preprocessing with coordinate system one and spatial alignment, productivity status index, grazing pressure index, time stability index and overload frequency index of each pixel in each year are constructed. The productivity status index, grazing pressure index, time stability index and overload frequency index are weighted and fused to generate grassland sustainability index raster for each year.

[0011] A pixel-by-pixel historical baseline is constructed using the grassland sustainability index raster of several consecutive years prior to the year to be detected. The grassland sustainability index of the year to be detected is compared pixel-by-pixel with the aforementioned pixel-by-pixel historical baseline to calculate the difference in grassland sustainability decline in each pixel.

[0012] The maximum inter-class variance method is used to determine the adaptive anomaly threshold for the year based on the distribution of valid values ​​in the grassland sustainability decline difference raster of the year to be detected, and the pixels whose grassland sustainability decline difference reaches or exceeds the adaptive anomaly threshold are marked as initial anomaly pixels.

[0013] Spatial neighborhood consistency filtering is performed on the initial abnormal pixels, connected region marking is performed on the filtered abnormal pixels, spatially connected abnormal pixels are aggregated into connected regions, connected regions are filtered according to the minimum patch size, and connected regions that reach the minimum patch size are retained to generate abnormal patches.

[0014] Output the spatial distribution raster and attribute data of the abnormal patches.

[0015] Furthermore, the methods for constructing the aforementioned productivity status indicators, grazing pressure indicators, time stability indicators, and overload frequency indicators include:

[0016] Based on the aboveground biomass raster data, the productivity status index was constructed using the quantile normalization method.

[0017] Based on the grazing intensity raster data, with a grazing intensity ratio of 1 as the grass-livestock balance benchmark, transition intervals are set on both sides of the benchmark, and the grazing pressure index is constructed using a piecewise linear response function.

[0018] Based on the aboveground biomass raster data, the interannual variation coefficient of each pixel is calculated using a rolling time window, and the time stability index is constructed by inverse standardization.

[0019] Based on the grazing intensity raster data, a rolling time window corresponding to the time stability index is used to count the proportion of years in which the grazing intensity ratio of each pixel exceeds the grassland carrying capacity within the rolling time window, and the overload frequency index is constructed by reverse transformation.

[0020] Furthermore, the lower limit of the transition interval is 0.8 times the grass-livestock balance benchmark, and the upper limit is 1.2 times the grass-livestock balance benchmark. The piecewise linear response function decreases linearly with the increase of the livestock intensity ratio within the transition interval.

[0021] Furthermore, before generating the annual grassland sustainability index raster by weighted fusion of productivity status indicators, grazing stress indicators, time stability indicators, and overgrazing frequency indicators, the following steps are also included:

[0022] An objective weighting method based on the indicator removal effect is adopted to determine the fusion weights of the productivity status indicator, grazing pressure indicator, time stability indicator, and overload frequency indicator.

[0023] Furthermore, the construction of a pixel-by-pixel historical baseline using grassland sustainability index rasters from consecutive years prior to the year to be detected includes:

[0024] The grassland sustainability index grid for the five consecutive years preceding the year to be tested is used as a historical window, which is rolled over year by year as the year to be tested.

[0025] For each pixel, count the number of valid grassland sustainability index values ​​within the historical window;

[0026] When the number of valid values ​​is not less than 4, the median of the valid values ​​is used as the historical baseline of the pixel.

[0027] When the number of valid values ​​is less than 4, the historical baseline of the pixel is set to invalid.

[0028] Furthermore, the calculation of the difference in grassland sustainability decline for each pixel includes:

[0029] When the historical baseline of a pixel is higher than the grassland sustainability index of the year to be detected, the difference between the historical baseline and the grassland sustainability index of the year to be detected is used as the difference in grassland sustainability decline of that pixel.

[0030] When the grassland sustainability index of the year to be tested is not lower than the historical baseline, the difference in grassland sustainability decline of that cell is set to 0;

[0031] When the grassland sustainability index of the historical baseline or the year to be tested is invalid, the difference in grassland sustainability decline of that cell is set to invalid.

[0032] Furthermore, the spatial neighborhood consistency filtering includes:

[0033] For each pixel, a 3×3 neighborhood window is constructed, and the number of abnormal pixels within the neighborhood window is counted.

[0034] If the center cell is the initial abnormal cell and the number of abnormal cells in the neighborhood window is not less than 4, retain the abnormal label of the center cell; otherwise, remove the abnormal label of the center cell.

[0035] Furthermore, during the neighborhood count process, invalid pixels are temporarily included in the calculation with a value of 0; after the neighborhood filtering is completed, invalid areas are restored based on the effective area mask of the initial abnormal pixel raster to avoid unobserved areas being judged as non-abnormal areas.

[0036] Furthermore, the connected region marking adopts the eight-neighbor connectivity rule, and adjacent abnormal pixels in the horizontal, vertical and diagonal directions are all determined to have spatial connectivity; the minimum patch size is 4 pixels.

[0037] Furthermore, the attribute data includes: the patch area and average degree of decline of each abnormal patch, as well as the average changes of the productivity status index, grazing pressure index, time stability index, and overload frequency index within the range of each abnormal patch.

[0038] The technical solution of the present invention has one or more of the following beneficial technical effects:

[0039] (1) By integrating four sub-indicators—productivity status, grazing pressure, time stability, and overgrazing frequency—a grassland sustainability index can be constructed. This index can simultaneously reflect grassland production supply capacity, grazing utilization pressure, production function fluctuation characteristics, and long-term overgrazing risk. It can identify the hidden degradation risk under the apparent stability of vegetation, distinguish between temporary decline and continuous degradation caused by climate fluctuations, and achieve a multi-dimensional comprehensive representation of grassland sustainability status, effectively making up for the representational defects of a single vegetation indicator.

[0040] (2) By constructing a rolling historical median baseline for each pixel, anomaly detection is transformed from a global absolute value comparison to a comparison of each pixel relative to its own historical state. Simultaneously, the maximum inter-class variance method is used to determine the adaptive anomaly threshold year by year, which can adapt to the differences in grassland background and the overall annual fluctuation characteristics in different regions, reducing the probability of misjudgment in low-productivity stable areas and missed detection in high-productivity rapidly declining areas. This significantly improves the adaptive capability and spatial adaptability of decline anomaly identification, avoiding the limitations of fixed thresholds and global benchmarks.

[0041] (3) Through three-level spatial processing of spatial neighborhood consistency filtering, eight-neighbor connected region marking and minimum patch size constraint, isolated noise pixels are eliminated, discrete abnormal pixels are aggregated into abnormal patches with clear boundaries and unique numbers, and attribute data containing patch area, degree of decline and changes in sub-indicators are generated simultaneously, which can directly support the accurate positioning, monitoring and early warning and management decision-making of grassland degradation areas.

[0042] Furthermore, other beneficial effects of the present invention will be mentioned in the specific embodiments. Attached Figure Description

[0043] Figure 1 This is an overall flowchart of an embodiment of the present invention;

[0044] Figure 2 This is a schematic diagram of the grassland sustainability index construction process according to an embodiment of the present invention;

[0045] Figure 3 This is a schematic diagram of the abnormal patch identification process according to an embodiment of the present invention;

[0046] Figure 4 This is a schematic diagram of the system structure according to an embodiment of the present invention. Detailed Implementation

[0047] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0048] To facilitate a clear description of the technical solutions in the embodiments of the present invention, the terms "first" and "second" are used to distinguish identical or similar items with essentially the same function and effect. Those skilled in the art will understand that the terms "first" and "second" do not limit the quantity or execution order.

[0049] The term "grassland sustainability index grid" refers to the basic spatial data characterizing the sustainable use status of grassland in this invention. This grid is generated by weighted fusion of productivity status indicators, grazing pressure indicators, temporal stability indicators, and overgrazing frequency indicators. Each pixel value reflects the overall sustainable use level of the grassland at that location. Higher values ​​indicate higher relative productivity, lower grazing pressure, stronger productivity stability, and lower overgrazing risk, indicating a better sustainable grassland use status. Lower values ​​suggest that the spatial pixel may have problems such as low productivity, strong interannual fluctuations, high grazing pressure, or frequent overgrazing.

[0050] The term "pixel-by-pixel historical baseline" refers to the spatialized data used in this invention to characterize the historical normal state of grassland sustainability for each spatial pixel. For each year to be detected, a historical window is formed by assembling grassland sustainability index rasters from several consecutive years preceding the year to be detected. For each spatial pixel, the central tendency (e.g., median) of all valid annual values ​​within the historical window is statistically analyzed, serving as the historical reference benchmark for that pixel in the year to be detected. By constructing a historical baseline independently for each pixel, anomaly detection is transformed from traditional global absolute value comparison to a relative comparison of each pixel relative to its own historical state, thereby eliminating the influence of differences in grassland background conditions across different regions on anomaly identification.

[0051] The term "adaptive anomaly threshold" is an annual segmentation criterion used in this invention to determine whether the decline in grassland sustainability has reached an "abnormal" level. For each year to be tested, a segmentation threshold adapted to the data distribution of that year is automatically calculated using the maximum inter-class variance method, based on the numerical distribution characteristics of all valid positive samples in the decline difference raster for that year. This threshold dynamically changes with the decline difference distribution of different years, rather than being a fixed value preset by humans. Its core function is to adaptively separate the relatively small values, which may belong to normal fluctuations, from the relatively large values, which have anomaly significance, in the decline difference samples of that year.

[0052] The term "minimum patch size" refers to the numerical standard used in this invention to constrain the minimum spatial range of anomalous patches, expressed in pixels. After connected region labeling, each connected region contains a different number of anomalous pixels. The minimum patch size is the lowest pixel threshold required to determine whether a connected region can be retained as a valid anomalous patch. Only connected regions with anomalous pixels reaching or exceeding this size are considered anomalous patches with practical management significance and are retained. Small connected regions that do not reach this size are considered scattered pixel combinations without significant spatial distribution characteristics and are discarded. This size value is related to the spatial resolution of the raster data; the corresponding pixel number standard can be adjusted for data of different resolutions.

[0053] The term "spatial neighborhood consistency filtering" is a technical means used in this invention to perform spatial constraint processing on initial anomalous pixels. Its basic approach involves constructing a sliding neighborhood window of a preset size (e.g., a 3×3 pixel window) for each pair of pixels, and counting the number of anomalous pixels within the window. If the central pixel is the initial anomalous pixel and the number of anomalous pixels within the window meets a preset condition, its anomalous label is retained; if the number of anomalous pixels within the window is insufficient, the central pixel is determined to be an isolated or noise-induced anomalous pixel, and its anomalous label is removed. This filtering method utilizes the characteristic that grassland anomalies are usually spatially continuous, using neighborhood consistency checks to eliminate isolated anomalous pixels formed by data noise, local missing data, and boundary perturbations.

[0054] The term "connected region labeling" refers to the processing method used in this invention to group spatially adjacent anomalous pixels into the same spatial unit. Based on preset pixel spatial adjacency relationships (such as the eight-neighbor connectivity rule, i.e., horizontal, vertical, and diagonal adjacency are all considered connected), all anomalous pixels are traversed. Anomalous pixels that satisfy the adjacency relationship are labeled as the same connected region, while spatially disconnected anomalous pixels are labeled as different connected regions. Through connected region labeling, discrete anomalous pixels are aggregated into spatially connected objects with clear boundaries, providing basic spatial units for subsequent patch size filtering and attribute statistics.

[0055] This embodiment takes the grassland region of the Qinghai-Tibet Plateau as the research object and provides an exemplary description of the grassland sustainability degradation abnormal patch identification method based on temporal grids described in this invention. Figure 1 This is an overall flowchart of an embodiment of the present invention. Figure 1 As shown, the overall process includes eight steps: basic data acquisition and preprocessing, grassland sustainability index construction, pixel rolling historical baseline construction, grassland sustainability decline difference calculation, adaptive threshold determination and initial abnormal pixel extraction, spatial neighborhood consistency filtering, connected region labeling and minimum patch filtering, and abnormal patch result generation and output. It can realize adaptive identification of grassland sustainability decline anomalies and spatialized patch output.

[0056] Step S1 involves basic data acquisition and preprocessing. In this embodiment, aboveground biomass (AGB) data and grazing intensity data for the study area from 2001 to 2019 are acquired. Both types of data have a spatial resolution of 250m. The aboveground biomass data is sourced from the National Tibetan Plateau Data Center, and the grazing intensity data is sourced from the Figshare database. Both types of data are raster data with time-series attributes. After acquiring the raw data, the two datasets are subjected to coordinate system one and spatial alignment processing to ensure that all raster data have a consistent spatial reference, spatial resolution, spatial extent, pixel position, and number of rows and columns, providing a unified spatial benchmark for subsequent pixel-by-pixel calculations. The term "coordinate system one" refers to the unified conversion to the metric projected coordinate system in WGS-84 (a geocentric coordinate system) to ensure a standard spatial resolution of 250 meters. This is to avoid the confusing use of geographic coordinates and projected coordinates, which is a common practice in this field.

[0057] Step S2 is to construct the Grassland Sustainability Index (GSI). Figure 2 This is a schematic diagram illustrating the grassland sustainability index construction process according to an embodiment of the present invention. Figure 2 As shown, this step constructs sub-indicators from four dimensions: productivity status, grazing pressure, time stability, and overgrazing frequency. After objective weighting, these indicators are integrated to generate an annual grassland sustainability index grid, which comprehensively represents the sustainable utilization status of grassland.

[0058] A productivity status index was constructed based on preprocessed multi-year aboveground biomass time series. To reduce the interference of extreme high and low values ​​on the index calculation results and improve the comparability of productivity status across different years and regions, quantile normalization was used to standardize the AGB (Above Ground Biomass). The calculation formula is as follows:

[0059] ;

[0060] Among them, AGB i,t AGB represents the aboveground biomass of the i-th spatial cell at time t; 5% The 5th percentile of all effective values ​​for the multi-year aboveground biomass grid in the study area; AGB 95% P represents the 95th percentile of all valid values ​​of the multi-year aboveground biomass raster in the study area; clip is a cutoff function used to restrict the calculation results to the interval [0,1], with values ​​outside the interval taking the corresponding boundary value; P i,trepresents the productivity state index of the i-th spatial pixel in the t-th year, and a larger value thereof indicates a higher relative productivity level of the corresponding spatial pixel. This normalization adopts the global quantile of multi-year above-ground biomass in the study area as a benchmark, which can ensure that the productivity state indexes of different years and different regions are horizontally comparable.

[0061] constructing a grazing pressure index, which is based on the preprocessed multi-year grazing intensity time series, and converts the corresponding relationship between the actual grazing intensity and grassland carrying capacity into a grazing pressure sustainability score according to the grass-livestock balance concept and ecological threshold theory. To reflect the buffering capacity of the grassland ecosystem to moderate grazing disturbance and avoid sudden changes in evaluation results caused by a single carrying threshold, in this embodiment, the grazing intensity R i,t =1 is taken as the grass-livestock balance benchmark (that is, the ratio of the actual grazing intensity to the grassland carrying capacity is 1), a 20% transition interval, that is, 0.8 to 1.2, is set above and below the carrying threshold, and a piecewise linear response function is adopted to construct the grazing pressure index G i,t , the calculation method is:

[0062] ;

[0063] wherein R represents the livestock intensity ratio, that is, the ratio of the actual grazing intensity to the grassland carrying capacity; when R ≤ 0.8, the grassland is in a relatively low-pressure utilization state; when R ≥ 1.2, it indicates that grazing utilization significantly exceeds the grassland carrying capacity; when 0.8 < R < 1.2, it indicates that the grassland transitions from moderate utilization to overloading risk. A larger value of the grazing pressure index G i,t indicates that the grazing pressure of the corresponding spatial pixel is lower, and the sustainable utilization state of the grassland is better.

[0064] constructing a temporal stability index, which is based on the preprocessed multi-year AGB time series, and uses a rolling time window of a preset length (that is, the index calculation rolling window) to calculate the interannual fluctuation degree of the grassland production function of each spatial pixel year by year, so as to construct a temporal stability index time series. In this embodiment, a 5-year rolling time window is adopted, specifically, the AGB data from 2001 to 2005 is used to calculate the temporal stability index of 2005, the data from 2002 to 2006 is used to calculate the temporal stability index of 2006, and so on, rolling backward year by year, until the data from 2015 to 2019 is used to calculate the temporal stability index of 2019.

[0065] For the index year t, five consecutive years of AGB data ending in that year are selected, and the coefficient of variation (CV) of the i-th spatial pixel within the rolling time window is calculated. To reduce the influence of extreme values ​​on the index calculation results, the quantile normalization method is used to inversely standardize the CV values, constructing a time stability index S, calculated using the following formula:

[0066] ;

[0067] Among them, CV i,t CV represents the coefficient of variation of ground biomass in the i-th spatial cell within a rolling window ending in year t; 5%,t The 5th percentile of the coefficient of variation of all spatial pixels in the study area in year t; CV 95%,t is the 95th percentile of the coefficient of variation of all spatial pixels in the study area in year t; clip is the cutoff function used to restrict the calculation results to the interval [0,1]; S i,t This represents the temporal stability index of the i-th spatial pixel in year t. The larger the value, the smaller the interannual fluctuation of the grassland productivity of the corresponding spatial pixel, and the stronger its temporal stability. This normalization uses the quantile of the coefficient of variation of the study area for the corresponding year as the benchmark, which can accurately reflect the relative level of productivity stability of each pixel in that year.

[0068] An overgrazing frequency index is constructed based on a preprocessed multi-year grazing intensity time series. A rolling time window corresponding to the time stability index is used to statistically analyze the number of times grazing pressure exceeds grassland carrying capacity for each spatial pixel within the rolling time window each year, thus constructing a time series for the sustainability index of overgrazing frequency. This embodiment also uses a 5-year rolling window. The overgrazing frequency index for 2005 is calculated using grazing intensity data from 2001 to 2005, and the overgrazing frequency index for 2006 is calculated using data from 2002 to 2006, rolling forward year by year until 2019. For the index year t, five consecutive years of livestock intensity data ending in that year are selected to determine the frequency of overgrazing for that spatial pixel within those five years. The formula for calculating the sustainability index of overgrazing frequency is:

[0069] ;

[0070] Where #(·) is a counting function used to count the number of years that meet the conditions; This represents the number of times within 5 years that the livestock intensity ratio exceeds the carrying capacity, where k ranges from 0 to 4, corresponding to 5 consecutive years ending in year t; F i,tThe value represents the overgrazing frequency index of the i-th spatial cell in year t. The larger the value, the lower the frequency of overgrazing and the closer the long-term grazing utilization is to sustainability; conversely, the smaller the value, the higher the risk of overgrazing for the corresponding spatial cell.

[0071] After constructing the four sub-indicators, the weights of each indicator were determined. An evaluation matrix was formed by combining the productivity status indicator, grazing stress indicator, temporal stability indicator, and overgrazing frequency indicator. The weights of each indicator were determined using the Method based on the Removal Effects of Criteria (MEREC). This method measures the importance of an indicator by calculating the change in the overall evaluation result before and after its removal. The greater the change in the evaluation result caused by the removal of an indicator, the greater its contribution to the grassland sustainability evaluation, and the higher its weight. The weights of the productivity status indicator, grazing stress indicator, temporal stability indicator, and overgrazing frequency indicator were calculated using this method and denoted as α, β, θ, and ω, respectively, satisfying α + β + θ + ω = 1.

[0072] Finally, based on the four sub-indicators and their corresponding weights, a weighted comprehensive method was used to calculate the grassland sustainability index from 2005 to 2019 to characterize the sustainable use status of grassland in different spatial pixels. The formula for calculating the grassland sustainability index of the i-th spatial pixel in year t is:

[0073] ;

[0074] Among them, GSI i,t P represents the grassland sustainability index of the i-th spatial pixel in year t; i,t G represents the productivity status index of the i-th spatial cell in year t; i,t S represents the grazing pressure index of the i-th spatial cell in year t; i,t F represents the temporal stability index of the i-th spatial pixel in year t; i,t Let represent the overload frequency index of the i-th spatial pixel in year t; α, β, θ and ω represent the weights of the productivity status index, grazing pressure index, time stability index and overload frequency index, respectively.

[0075] A higher GSI value indicates that the corresponding spatial pixel has a higher relative productivity level, lower grazing pressure, stronger production stability, and lower overgrazing risk, indicating a better sustainable grassland utilization status. Conversely, a lower GSI value suggests that the spatial pixel may have problems such as low productivity, strong interannual fluctuations, high grazing pressure, or frequent overgrazing. Based on the calculation results of each spatial pixel, spatial distribution data of grassland sustainability index is generated.

[0076] Step S3 is to construct the historical baseline for cell rolling. Figure 3 This is a schematic diagram of the abnormal patch identification process according to an embodiment of the present invention. Figure 3 As shown, this step uses a 5-year rolling historical window (i.e., a historical baseline rolling window) combined with the median to construct the historical baseline, which can reduce the impact of individual anomalous years and extreme values ​​on the representation of grassland historical state. Following the annual order, the GSI grids of the five consecutive years preceding the year to be detected are used as the historical window, and the historical window rolls annually with the year to be detected. Since the historical baseline needs to rely on the index sequence of the five consecutive years preceding the year to be detected, in this embodiment, when the input year range is 2005 to 2019, the range of anomalous detection years is 2010 to 2019. Specifically, when 2010 is the year to be detected, the historical window is 2005 to 2009; when 2011 is the year to be detected, the historical window is 2006 to 2010, and so on, rolling forward year by year until 2014 to 2018 is used as the historical window for 2019.

[0077] For the i-th pixel, count the number of valid GSI values ​​within the historical window. When there are at least 4 valid annual values, use the median of the valid annual values ​​as the historical baseline of the pixel. When there are fewer than 4 valid annual values, set the historical baseline of the pixel in the year to be detected to an invalid value. Finally, generate a pixel-by-pixel historical baseline raster corresponding to the year to be detected.

[0078] Step S4 involves calculating the difference in grassland sustainability degradation. The GSI raster for the year in question is compared pixel-by-pixel with its corresponding historical baseline raster to calculate the difference in grassland sustainability degradation for each pixel. The calculation formula is as follows:

[0079] ;

[0080] Among them, D i,t For the difference in grassland sustainability decline in the i-th pixel in year t, B i,t GSI is the historical baseline of the i-th pixel in year t. i,t This refers to the grassland sustainability index (GSI) for the year to be tested. During the calculation, when the GSI for the year to be tested is lower than the historical baseline, the difference between the historical baseline and the current year's index is retained to characterize the degree of decline of that pixel relative to its historical state. When the GSI for the year to be tested is not lower than the historical baseline, it indicates that the current state has not declined, and the corresponding decline difference value is set to 0. When the historical baseline or the current year's index is invalid, the corresponding decline difference is set to invalid. A grassland sustainability decline difference raster for the year to be tested is generated according to the above rules, so that it only retains the negative deviation information of the GSI relative to the historical baseline.

[0081] Step S5 involves adaptive threshold determination and initial outlier pixel extraction. The grassland sustainability degradation difference raster generated in step S4 is read, and pixels with values ​​greater than 0 and valid are selected as the adaptive threshold calculation objects. To reduce the computational load of large-scale raster data, regular sampling is performed according to the raster spatial location, with a target sampling quantity of 50,000 pixels; when the total number of raster pixels is less than 50,000, the actual number of pixels is used as the sampling quantity.

[0082] When the number of valid samples obtained by regular sampling is less than 1000, random sampling is used to extract a maximum of 10000 valid pixels from the descent difference raster. After sampling, invalid values, non-finite values ​​(including infinite values, null values, and other outliers), and samples less than or equal to 0 are further deleted, retaining only positive value samples that can characterize the decline in grassland sustainability. The obtained descent difference samples are divided into 128 histogram intervals, and the number of samples in the k-th histogram interval is denoted as h. k Then the probability of this interval is:

[0083] ;

[0084] Among them, h k h is the number of decreasing difference samples contained in the k-th histogram interval; K is the total number of histogram intervals, h j p represents the number of decreasing difference samples contained in the j-th histogram interval; j is the traversal index of the histogram interval, ranging from 1 to K; in this embodiment, K=128; k denoted as the sample probability corresponding to the k-th histogram interval.

[0085] Let the center value of the k-th histogram interval be c. k For a candidate segmentation position q, calculate the cumulative probability of the first q intervals:

[0086] ;

[0087] Where q is the candidate segmentation position, corresponding to the interval index of the histogram; ω q Let p be the cumulative probability of the first q intervals. k denoted as the sample probability corresponding to the k-th histogram interval.

[0088] And cumulative average:

[0089] ;

[0090] Where q is the candidate segmentation position, c k μ is the center value of the k-th histogram interval. q Let p be the cumulative mean of the first q intervals. k denoted as the sample probability corresponding to the k-th histogram interval.

[0091] The population mean of all samples is:

[0092] ;

[0093] Where K is the total number of histogram intervals, μ T c is the population mean of all samples showing a decrease in difference. k p is the center value of the k-th histogram interval; k denoted as the sample probability corresponding to the k-th histogram interval.

[0094] Given that both sides of the candidate segmentation position contain valid samples, calculate the corresponding inter-class variance:

[0095] ;

[0096] in, ω represents the inter-class variance corresponding to the candidate segmentation position q; q μ is the cumulative probability of the first q intervals; q μ is the cumulative mean of the first q intervals; T This represents the population mean of all samples.

[0097] Iterate through all candidate segmentation positions and select the position with the largest inter-class variance:

[0098] ;

[0099] Where, q * The optimal candidate segmentation position that maximizes the inter-class variance; q represents the inter-class variance corresponding to the candidate split position q; arg max represents the value of the independent variable that maximizes the objective function, where the independent variable is the candidate split position q.

[0100] The boundary value of the histogram interval corresponding to the optimal segmentation position is determined as the adaptive anomaly threshold for the year t to be detected.

[0101] ;

[0102] Where, τ t Let b be the adaptive anomaly threshold for year t; b is the sequence of boundary values ​​for the histogram interval. q*+1 The optimal candidate segmentation position q * The corresponding left boundary value of the next histogram interval. This method essentially utilizes the histogram distribution of the decreasing difference samples to adaptively separate smaller decreasing values ​​from larger decreasing values ​​through the maximum inter-class variance. Different years determine their thresholds based on their respective decreasing difference grids, thereby avoiding the use of the same fixed threshold for all years.

[0103] This step also sets up multiple fallback mechanisms to ensure the stability of the results. When all valid sample values ​​are the same, the sample value is directly used as the segmentation threshold. When the maximum effective inter-class variance cannot be obtained, the 95th percentile of the decreasing difference sample is used as the backup threshold. When the number of valid positive samples is less than the preset minimum number, the year is set as the threshold for no valid outliers, and the valid region is marked as non-outlier to avoid unstable segmentation caused by a small number of samples.

[0104] After obtaining the adaptive anomaly threshold, the decreasing difference raster is judged pixel by pixel. When the i-th pixel shows a decrease in grassland sustainability D in the year t to be detected, the threshold is determined. i,t A valid value that is greater than 0 and reaches or exceeds the adaptive anomaly threshold τ for that year. t At that time, the initial anomaly label A of the pixel is set. i,t A value of 1 indicates that the pixel is identified as an initial anomalous pixel. When the decrease difference is a valid value but less than the adaptive anomaly threshold, or when the decrease difference is equal to 0, the initial anomaly marker is assigned a value of 0, indicating that the pixel is a non-anomaly pixel. When the decrease difference is an invalid value, the corresponding initial anomaly marker is set to an invalid value and does not participate in subsequent anomaly identification. Initial anomaly rasters corresponding to the spatial location of grassland sustainability decrease difference rasters are generated according to the above rules.

[0105] Step S6 is spatial neighborhood consistency filtering. The initial anomalous cell raster generated in step S5 is read, and a 3×3 sliding neighborhood window is constructed centered on each cell. Focus statistics are used to calculate the number of initial anomalous cells within this window. The calculation formula for the number of anomalous neighboring cells of the i-th cell in the year t to be detected is:

[0106] ;

[0107] Among them, M i,t Ω represents the number of abnormal pixels in a 3×3 neighborhood centered at the i-th pixel, where j is the traversal index of pixels within the neighborhood, used to refer to any pixel within the neighborhood; i Let i be the spatial neighborhood centered on the i-th pixel, including the central pixel and its 8 surrounding adjacent pixels. These are binary pixel values ​​used for neighborhood statistics, where the initial abnormal pixels are set to 1, and the non-abnormal and invalid pixels are set to 0 during the neighborhood number statistics process.

[0108] Based on the initial abnormal state of the central pixel and the number of abnormal pixels in its neighborhood, the initial abnormal pixels are filtered. If the central pixel is the initial abnormal pixel and the number of abnormal pixels in its 3×3 neighborhood is not less than 4, the abnormal label of the pixel is retained; otherwise, its abnormal label is updated to 0.

[0109] During neighborhood statistics, invalid pixels are only temporarily counted as having zero entries. After neighborhood filtering is completed, invalid regions are restored based on the effective region mask of the initial anomaly pixel raster to avoid misidentifying unobserved regions as non-anomaly regions. The filtering results are written into a binary raster with the same georegistration information as the initial anomaly pixel raster to generate a spatial neighborhood-filtered anomaly raster, which serves as input for subsequent connected region labeling and anomaly patch generation.

[0110] This step improves the spatial continuity and reliability of grassland sustainability degradation anomaly identification results by utilizing the continuity of grassland anomalies between adjacent spatial pixels to remove isolated anomalies and scattered noise that lack support from surrounding anomaly pixels, while retaining anomaly regions with certain spatial clustering characteristics.

[0111] Step S7 involves connected region labeling and minimum patch filtering. The spatial neighborhood filtering anomaly raster generated in step S6 is read, and the anomaly pixels are used as connected region identification objects. Non-anomaly pixels and invalid pixels do not participate in connected region labeling.

[0112] Furthermore, this embodiment employs an eight-neighbor connectivity rule to traverse anomalous pixels. When two anomalous pixels are adjacent in the horizontal, vertical, or diagonal direction, they are determined to have a spatial connectivity relationship, and the interconnected anomalous pixels are grouped into the same connected region. Let C be the c-th connected region in the year t to be detected. c,t The number of anomalous pixels it contains is:

[0113] ;

[0114] Where, n c,t Let be the number of anomalous pixels contained in the c-th connected region in year t; i is the traversal index of the raster pixels in the study area, used to refer to any spatial pixel; C c,t Let I be the c-th connected region in year t; I(·) is the indicator function, which takes the value 1 when the i-th cell belongs to the connected region, and takes the value 0 otherwise.

[0115] Minimum patch filtering is performed based on the number of pixels in each connected region. In this embodiment, the minimum patch pixel count is set to 4. A connected region containing at least 4 anomalous pixels is retained, while a region containing fewer than 4 anomalous pixels is deleted. For data with a spatial resolution of 250m, the nominal area corresponding to 4 pixels is approximately 25 hectares.

[0116] This step organizes adjacent anomalous cells into spatially continuous patches using eight-neighbor connectivity markers and removes small, scattered anomalous areas through minimum scale filtering, thereby improving the spatial integrity of the anomalous results and forming grassland sustainability degradation anomalous patches with clear locations, boundaries, and numbers.

[0117] Step S8 generates and outputs abnormal patch results. Based on the connected region markers and minimum patch filtering results obtained in step S7, a raster of abnormal patch numbers for grassland sustainability decline is generated. Abnormal pixels belonging to the same connected region are assigned the same integer patch number, while spatially disconnected abnormal patches are assigned different patch numbers to characterize the spatial location, boundary range, and mutual distinction of each abnormal patch.

[0118] Furthermore, this embodiment can simultaneously generate a binary anomaly raster based on the anomaly patch number raster, where anomaly patch pixels are assigned a value of 1, non-anomaly valid pixels are assigned a value of 0, and pixels without valid observations are retained as invalid values. The anomaly patch number raster serves as the main output of this invention, while the binary anomaly raster serves as an auxiliary result. Using the anomaly patch number as a partition identifier, the GSI raster, historical baseline raster, decline difference raster, and productivity status index P, grazing pressure index G, time stability index S, and overload frequency index F are partitioned and statistically analyzed to generate an anomaly patch attribute table. The attribute data includes patch area, average decline degree, and the average change of the four component indices: productivity status, grazing pressure, time stability, and overload frequency.

[0119] Furthermore, the present invention can also be implemented in the form of a system. Figure 4 This is a schematic diagram of the system structure according to an embodiment of the present invention. Figure 4 As shown, the electronic device 310 (e.g., a computer server with program execution capabilities) includes at least one processor 311, a power supply 314, and a memory 312 and an input / output interface 313 communicatively connected to the at least one processor 311. The memory 312 stores instructions executable by the at least one processor 311, which are executed by the at least one processor 311 to enable the at least one processor 311 to perform the methods disclosed in any of the foregoing embodiments. The input / output interface 313 may include a display, keyboard, mouse, and USB interface for inputting and outputting data. The power supply 314 provides power to the electronic device 310.

[0120] Those skilled in the art will understand that all or part of the steps of the above method embodiments can be implemented by hardware related to program instructions. The aforementioned program can be stored in a computer-readable storage medium. When the program is executed, it performs the steps of the above method embodiments. The aforementioned storage medium includes various media that can store program code, such as mobile storage devices, read-only memory, magnetic disks, or optical disks.

[0121] When the integrated units of this invention are implemented as software functional units and sold or used as independent products, they can also be stored in a computer-readable storage medium. Based on this understanding, the technical solutions of the embodiments of this invention, or the parts that contribute to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as mobile storage devices, read-only memories, magnetic disks, or optical disks.

[0122] It should be noted that the parameters such as the historical window length, number of valid years, number of adaptive threshold samples, neighborhood window size, neighborhood retention conditions, connectivity direction, and minimum patch size in the above embodiments can be adjusted according to the data time span, spatial resolution, and characteristics of the study area. The adjustment methods are all conventional choices made by those skilled in the art under the technical concept of this invention and do not depart from the protection scope of this invention.

[0123] To better illustrate the present invention, numerous specific details have been provided in the detailed embodiments described above. Those skilled in the art should understand that the present invention can be practiced even without certain specific details. In some instances, methods, means, components, and circuits well known to those skilled in the art have not been described in detail in order to highlight the main points of the present invention.

[0124] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for identifying anomalous patches of declining grassland sustainability based on temporal grids, characterized in that, include: The aboveground biomass raster data and grazing intensity raster data of the area to be detected over many years are obtained. After preprocessing with coordinate system one and spatial alignment, productivity status index, grazing pressure index, time stability index and overload frequency index of each pixel in each year are constructed. The productivity status index, grazing pressure index, time stability index and overload frequency index are weighted and fused to generate grassland sustainability index raster for each year. A pixel-by-pixel historical baseline is constructed using the grassland sustainability index raster of several consecutive years prior to the year to be detected. The grassland sustainability index of the year to be detected is compared pixel-by-pixel with the aforementioned pixel-by-pixel historical baseline to calculate the difference in grassland sustainability decline in each pixel. The maximum inter-class variance method is used to determine the adaptive anomaly threshold for the year based on the distribution of valid values ​​in the grassland sustainability decline difference raster of the year to be detected, and the pixels whose grassland sustainability decline difference reaches or exceeds the adaptive anomaly threshold are marked as initial anomaly pixels. Spatial neighborhood consistency filtering is performed on the initial abnormal pixels, connected region marking is performed on the filtered abnormal pixels, spatially connected abnormal pixels are aggregated into connected regions, connected regions are filtered according to the minimum patch size, and connected regions that reach the minimum patch size are retained to generate abnormal patches. Output the spatial distribution raster and attribute data of the abnormal patches.

2. The method for identifying anomalous patches of declining grassland sustainability based on temporal grids according to claim 1, characterized in that, The methods for constructing the aforementioned productivity status indicators, grazing stress indicators, time stability indicators, and overload frequency indicators include: Based on the aboveground biomass raster data, the productivity status index was constructed using the quantile normalization method. Based on the grazing intensity raster data, with a grazing intensity ratio of 1 as the grass-livestock balance benchmark, transition intervals are set on both sides of the benchmark, and the grazing pressure index is constructed using a piecewise linear response function. Based on the aboveground biomass raster data, the interannual variation coefficient of each pixel is calculated using a rolling time window, and the time stability index is constructed by inverse standardization. Based on the grazing intensity raster data, a rolling time window corresponding to the time stability index is used to count the proportion of years in which the grazing intensity ratio of each pixel exceeds the grassland carrying capacity within the rolling time window, and the overload frequency index is constructed by reverse transformation.

3. The method for identifying abnormal patches of declining grassland sustainability based on temporal grids according to claim 2, characterized in that: The lower limit of the transition interval is 0.8 times the grass-livestock balance benchmark, and the upper limit is 1.2 times the grass-livestock balance benchmark. The piecewise linear response function decreases linearly with the increase of the livestock intensity ratio within the transition interval.

4. The method for identifying anomalous patches of declining grassland sustainability based on temporal grids according to claim 1, characterized in that, Before generating the annual grassland sustainability index raster by weighted fusion of productivity status indicators, grazing stress indicators, time stability indicators, and overgrazing frequency indicators, the following steps are also included: An objective weighting method based on the indicator removal effect is adopted to determine the fusion weights of the productivity status indicator, grazing pressure indicator, time stability indicator, and overload frequency indicator.

5. The method for identifying abnormal patches of declining grassland sustainability based on temporal grids according to claim 1, characterized in that, The construction of a pixel-by-pixel historical baseline using grassland sustainability index rasters from consecutive years prior to the year to be detected includes: The grassland sustainability index grid for the five consecutive years preceding the year to be tested is used as a historical window, which is rolled over year by year as the year to be tested. For each pixel, count the number of valid grassland sustainability index values ​​within the historical window; When the number of valid values ​​is not less than 4, the median of the valid values ​​is used as the historical baseline of the pixel. When the number of valid values ​​is less than 4, the historical baseline of the pixel is set to invalid.

6. The method for identifying anomalous patches of declining grassland sustainability based on temporal grids according to claim 1, characterized in that, The calculation of the differences in grassland sustainability decline for each pixel includes: When the historical baseline of a pixel is higher than the grassland sustainability index of the year to be detected, the difference between the historical baseline and the grassland sustainability index of the year to be detected is used as the difference in grassland sustainability decline of that pixel. When the grassland sustainability index of the year to be tested is not lower than the historical baseline, the difference in grassland sustainability decline of that cell is set to 0; When the grassland sustainability index of the historical baseline or the year to be tested is invalid, the difference in grassland sustainability decline of that cell is set to invalid.

7. The method for identifying abnormal patches of declining grassland sustainability based on temporal grids according to claim 1, characterized in that, The spatial neighborhood consistency filtering includes: For each pixel, a 3×3 neighborhood window is constructed, and the number of abnormal pixels within the neighborhood window is counted. If the center cell is the initial abnormal cell and the number of abnormal cells in the neighborhood window is not less than 4, retain the abnormal label of the center cell; otherwise, remove the abnormal label of the center cell.

8. The method for identifying anomalous patches of declining grassland sustainability based on temporal grids according to claim 7, characterized in that: During the neighborhood count, invalid pixels are temporarily assigned a value of 0 for calculation. After neighborhood filtering is completed, invalid regions are restored based on the effective region mask of the initial abnormal pixel raster to prevent unobserved regions from being judged as non-abnormal regions.

9. The method for identifying abnormal patches of declining grassland sustainability based on temporal grids according to claim 1, characterized in that: The connected region marking adopts the eight-neighbor connectivity rule, and adjacent abnormal pixels in the horizontal, vertical and diagonal directions are all determined to have spatial connectivity; the minimum patch size is 4 pixels.

10. The method for identifying anomalous patches of grassland sustainability decline based on temporal grids according to claim 1, characterized in that, The attribute data includes: The patch area and average degree of decline of each abnormal patch, as well as the average changes of the productivity status index, grazing pressure index, time stability index and overload frequency index within the range of each abnormal patch.

Images