Grassland type intelligent identification and space-time monitoring method for grassland ecology

Abstract

The invention relates to the technical field of mobile robot environment perception and autonomous navigation, in particular to a grassland type intelligent identification and space-time monitoring method for grassland ecology, which comprises the steps of data preparation, feature processing, sample construction, model training and application and monitoring result analysis. According to the method, autonomous extraction of high-order topographic features such as surface roughness and effective integration of the high-order topographic features and multi-source remote sensing features such as spectrum, elevation and texture are realized, and the sensitivity and recognition capability of the system to grassland type differences under different landform backgrounds are enhanced; a channel attention mechanism and a multi-scale residual structure are applied to grassland remote sensing classification and segmentation, the automatic attention capability of a model on key features is improved, meanwhile, the accuracy of deep space information expression and boundary discrimination is enhanced, and through the multi-temporal earth surface remote sensing monitoring and change detection process of the system, the grassland remote sensing classification and segmentation efficiency is improved. Spatio-temporal evolution results such as annual transfer, degradation and recovery of grassland types can be automatically output.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of mobile robot environment perception and autonomous navigation, and particularly relates to a grassland type intelligent identification and spatio-temporal monitoring method for grassland ecology. BACKGROUND

[0002] Remote sensing image is a film, photo or digital image recorded by sensor, mainly including aerial photograph and satellite photograph, which is widely used in grassland resource monitoring, degradation diagnosis and dynamic evaluation in the field of ecological management. The current grassland type intelligent identification and spatio-temporal monitoring method mainly includes remote sensing classification method based on single spectral information or vegetation index (such as NDVI) as the core criterion, and technical path of applying machine learning algorithm (such as support vector machine SVM, random forest RF) to identify grassland type. Convolutional neural network represented by UNet is gradually applied to remote sensing image feature segmentation, usually taking multi-channel spectral band or vegetation index as model input to realize automatic classification of grassland type.

[0003] The above technical solutions mostly only use the data characteristics of multi-spectral remote sensing image, rely on single spectrum or vegetation index characteristics, and are difficult to effectively fuse multi-source heterogeneous information such as elevation, terrain and texture, which limits the classification accuracy in areas with complex ground types and strong terrain heterogeneity. For large-scale, complex terrain and strong heterogeneity grassland-desert areas, there are still problems of single feature expression, limited model generalization ability and insufficient time series dynamic monitoring ability. In addition, due to the limited feature channel expression ability of deep neural network structure such as original UNet model, it is difficult to automatically focus on important ground areas, and the spatial resolution is gradually lost in the down-sampling process, which affects the fine positioning of ground object boundary, and it is difficult to realize high-precision, automatic, multi-temporal grassland coverage type identification and dynamic change tracking, which cannot provide strong data support for actual needs such as grassland degradation diagnosis and ecological restoration management. Therefore, a grassland type intelligent identification and spatio-temporal monitoring method for grassland ecology is proposed. SUMMARY

[0004] In order to solve the above technical problems in the prior art, the present application provides a grassland type intelligent identification and spatio-temporal monitoring method for grassland ecology.

[0005] In order to solve the above technical problems, the present application provides the following technical scheme: a grassland type intelligent identification and spatio-temporal monitoring method for grassland ecology, the grassland type intelligent identification and spatio-temporal monitoring method comprising the following steps:

[0006] S1: Collect remote sensing image data of Landsat8 in the study area during the cloudless or extremely low cloud coverage period of each growing season, collect global digital elevation model data, number and record the metadata of the collected remote sensing image data, and select representative annotation samples;

[0007] S2: Perform radiation calibration, image registration and atmospheric correction preprocessing on the collected original remote sensing image data, convert the original remote sensing image data value to actual ground reflectivity, eliminate atmospheric and instrument interference, then crop the image within the vector boundary of the study area, remove irrelevant areas, unify the resolution, and perform cloud removal processing;

[0008] S3: Perform spectral, ground ruggedness and texture feature extraction, and normalize and dimensionless process all extracted features, then correspondingly register in space to form a multi-channel feature image group, preparing for subsequent model input;

[0009] S4: Establish a segmentation model and perform model training, use the encoding-decoding network structure of UNet as the backbone of the segmentation model, optimize the characteristics of the ground objects in the study area within the resolution range of the remote sensing image, and introduce a dynamic sampling mechanism in the decoding stage to restore the features, and then perform feature fusion after restoration;

[0010] S5: After the training of SE-ResUNet is completed, use the optimized SE-ResUNet deep segmentation model to perform end-to-end automatic inference and segmentation on the multi-temporal and multi-modal remote sensing image data of the entire area, and monitor and analyze the dynamic changes of the grassland.

[0011] Preferably, the step S1 specifically comprises the following contents:

[0012] S11: The record content includes the acquisition date, the center coordinates, the image range, the sensor parameters and the cloud coverage;

[0013] S12: Select an image or elevation data with strict overlap of all multi-temporal and multi-source data as a spatial reference, and define it as a representative annotation sample.

[0014] Preferably, the step S2 comprises the following processing in sequence:

[0015] S21, radiation calibration:

[0016] The metadata in the MTL file attached to the original remote sensing image data of Landsat8 is converted into physical radiance and ground reflectivity, and the calibration formula is as follows:

[0017] L λ =Μ L ·Q cal +Α L

[0018] wherein, L λ denotes the calibrated radiance, M L denotes the band calibration coefficient, Q cal denotes the original DN value, A L denotes the addition correction term;

[0019] S22, atmospheric correction:

[0020] The top layer radiance is corrected to the actual ground reflectance by the fast atmospheric correction method, and the remote sensing professional software is used as the operation platform to maximize the elimination of signal confusion caused by atmospheric scattering and absorption signals;

[0021] S23, image registration:

[0022] The geometric correction is performed on the remote sensing image data of each period of Landsat8 based on DEM elevation data or regional high-precision maps, 30 or more uniformly distributed ground control points are selected, the registration accuracy is controlled within 1 pixel, and the output is unified to the WGS-84 UTM projection system;

[0023] S24, image cropping:

[0024] The spatial cropping is performed on all the images after the above pretreatment based on the shp format of the research area boundary vector;

[0025] S25, resampling:

[0026] The multi-source data is uniformly resampled and standardized to 30-meter resolution;

[0027] S26, cloud removal:

[0028] The Fmask algorithm is used to detect and mask the cloud or shadow of the Landsat8 image, the detection is performed by visible light, short-wave infrared band threshold judgment, and the Fmask is based on the rule-threshold process of multispectral and thermal infrared combination: first, the cloud is identified by using the high reflectivity of visible light and SWIR and the low brightness temperature of TIRS, then the cloud shadow is identified by using the solar geometric projection and NIR / SWIR darkening characteristics, and finally the stable cloud or shadow mask is formed through morphological processing, thereby generating a cloud mask layer, the cloud and cloud shadow area is set to invalid value, and is removed in subsequent feature statistics, after the removal, the hole filling is performed by using other phase images of the same year.

[0029] Preferably, the step S3 specifically comprises the following contents:

[0030] S31: extracting the ground vegetation coverage area in the image by using the spectral characteristics:

[0031] Information from each band of Landsat 8 was extracted, and the Enhanced Vegetation Index (EVI) for the study area was calculated using ENVI software. The specific calculation formula is as follows:

[0032]

[0033] Where, ρ NIR ρ R and ρ Blue These represent the surface reflectance in the near-infrared, red, and blue light bands, respectively.

[0034] By introducing the blue light band and utilizing its sensitivity to atmospheric aerosols, atmospheric effects are corrected by adding a gain coefficient, typically 2.5.

[0035] S32: Select a center point, extract the elevation range and average slope of the area near the center point, and introduce pixel descriptions to quantitatively express the influence of terrain conditions on vegetation distribution and type classification, thereby improving the classification and recognition ability under heterogeneous terrain. The extraction range is: with the target pixel as the center, statistics are performed within a multi-scale window, such as r∈{3×3, 5×5, 7×7, 11×11}, with 5×5 as the main scale. After parallel calculation of multi-scale features, the data is then stitched together. The specific calculation formula is as follows:

[0036] The actual coverage area S of the target object on the ground plane Surface The calculation formula is:

[0037]

[0038] The area S of the object projected onto the ground along the sun's direction or a specified direction of the sensor's line of sight. Projection The calculation formula is:

[0039] S Projection =(d-1) 2 ·r

[0040] The area S of the pixel in the i-th row and j-th column ij The calculation formula is:

[0041] S ij =S ΔABD +S ΔACD

[0042]

[0043]

[0044] Among them, S Surface S represents the surface area of ​​the extracted window T. ij S represents the area of ​​the (i,j)th pixel. ProjectionLet (i,j) represent the surface projection area of ​​the extraction window T, (i,j) represent the cell row and column indices, d represent the extraction window parameters, r represent the data resolution, and S represent the surface projection area of ​​the extraction window T. ΔABD and S ΔABD P1 and P2 represent the areas of the triangles formed by points A, B, D, and A, C, D, respectively. P1 and P2 represent weighting coefficients used for linear weighted fusion of multiple quantities. AB l BD l AD l AC l CD and l XY Both represent the Euclidean distance between the two points, h XY h AB h BD h AC h CD and h AD Both represent the relative degree between two points, t 11 The self-maintenance of the representative period t1, t 12 This represents a transition from t1 to t2, where t 21 This represents a rotation from t2 to t1, where t 22 The self-maintained state of the representative period t2;

[0045] S33: Extract image texture features:

[0046] Statistical methods were used to extract image texture features at different window scales. The required features specifically included mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment of angle, and correlation. The specific calculation formulas are as follows:

[0047] The mean is:

[0048] Mean = ∑ i ∑ j p(i,j)*i;

[0049] The variance is:

[0050] Variance=∑ i Σ j p(i,j)*(i-Mean) 2 ;

[0051] Homogeneity is:

[0052]

[0053] Contrast is:

[0054] Contrast=∑ i ∑ j p(i,j)*(ij)2 ;

[0055] Dissimilarity is:

[0056] Dissimilarity=∑ i ∑ j p(i,j)*|ij|;

[0057] Entropy is:

[0058] Entropy = ∑ i ∑ j p(i,j)*lnΡ(i,j);

[0059] The second moment of the angle, ASM, is:

[0060] ASM=∑ i ∑ j p(i,j) 2 ;

[0061] The correlation is:

[0062]

[0063] Where (i,j) represents the index of the quantized gray level, with a value range of {0, 1, ..., G-1}, G represents the total number of quantized levels, and P(i,j) represents the joint probability. P(i,j) is a normalized GLCM element calculated under a given window, orientation, and pixel spacing, representing the joint probability of gray level pair (i,j) occurring, satisfying ∑ i ∑ j P(i,j), where i represents the gray level value in the row direction.

[0064] Preferably, step S4 specifically includes the following:

[0065] Step S41: In each encoding and decoding module of the UNet encoder-decoder network structure, embed the SE channel attention mechanism to extract global information and recalibrate the features, so that the segmentation model automatically focuses on the channel layer that is most representative of the land surface type distinction.

[0066] Step S42: Introduce multi-layer residual connection modules into the feature pathways at each level in the segmentation model to alleviate gradient vanishing and information decay that occur during the forward and backward propagation of deep features;

[0067] Step S43: Align all preprocessed and normalized features pixel by pixel in space, and fuse Landsat8 spectral bands, EVI, DEM elevation, surface ruggedness, image textures and statistical features at all levels in a channel stitching manner. Channel stitching connects multi-source, multi-scale or multi-type features in the channel dimension as input or intermediate layer features of the neural network. After fusion, the features are stitched into different channels of the input layer to form a multi-channel input tensor.

[0068] Step S44: Combining known field sampling data, existing grassland or land cover thematic databases, and expert manual interpretation results, construct a high-precision supervised classification and labeling sample set, label the main land surface types, and all samples strictly follow the unified interpretation standards and spatial consistency.

[0069] Step S45: A dynamic sampling mechanism is introduced in the decoding stage to achieve accurate feature restoration. Through the skipconnection mechanism, the shallow spatial details at the encoding end are directly fused with the high-level semantic features at the decoding end.

[0070] Preferably, the automated reasoning and segmentation process in step S5 includes the following:

[0071] Step S51: Input the integrated multi-channel remote sensing and auxiliary terrain feature images into the trained segmentation model in the form of a standardized three-dimensional tensor. Specific features include Landsat8 multi-band, EVI, DEM elevation, surface ruggedness and texture features.

[0072] Step S52: The segmentation model automatically identifies the spatial feature combination of each input cell, and extracts, fuses and discriminates features layer by layer through the backbone network and attention mechanism, and outputs a cell-level grassland classification result that is consistent with the original input space. Cell-level grassland classification refers to the category determination result of each raster cell, that is, the output is a classification map that is consistent with the input raster in spatial size and assigns values ​​to each cell.

[0073] Step S53: Align all processed annual category raster layer data spatially with the initial image in chronological order to continuously reconstruct the spatial change process of grassland types, and then perform automated change detection on grassland cover types between years.

[0074] Preferably, the automated change detection of grassland cover type between years in step S53 specifically includes the following:

[0075] Step S531: Using the classified annual raster layer data, statistically analyze the changes in grassland type for the same pixel across different years:

[0076] Step S532: Construct a grassland type transfer matrix to quantify the relationship between the increase, loss and transformation of different grassland types, and achieve accurate measurement of the annual changes in the area of ​​different types;

[0077] Step S533: Automatically generate a spatial distribution map of grassland type changes, which visually shows the degradation path of high-coverage grassland to low-coverage or bare land and sandy land, as well as the restoration and evolution process of some areas.

[0078] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0079] 1. This invention achieves autonomous extraction of high-order terrain features such as surface ruggedness and effective integration with multi-source remote sensing features such as spectrum, elevation, and texture. It enhances the system's sensitivity and recognition ability to grassland type differences under different terrain and landform backgrounds. Through the improved SE-ResUNet deep learning framework, channel attention mechanism and multi-scale residual structure are applied to grassland remote sensing classification and segmentation, improving the model's ability to automatically focus on key features. At the same time, it also enhances the accuracy of deep spatial information expression and boundary discrimination. Through the system's multi-temporal surface remote sensing monitoring and change detection process, it can automatically output the spatiotemporal evolution results of annual grassland type transfer, degradation, and recovery.

[0080] 2. This invention achieves automatic tracking and quantitative analysis of annual dynamic changes in grassland cover types, degradation and restoration directions by deeply integrating remote sensing data from multiple periods. It also provides multi-dimensional information on grassland degradation, desertification expansion, and ecological restoration responses in a spatiotemporal manner, providing scientific, accurate, open and transparent data support for decision-making and management in regional ecological protection, grassland management, and grassland resource evaluation. Furthermore, it can automatically generate land cover type transfer maps and statistical reports, significantly improving the intelligence and efficiency of grassland resource surveys and reducing the workload and subjective errors of manual interpretation.

[0081] 3. This invention utilizes a novel method for intelligent identification and change analysis of grassland types by integrating multimodal remote sensing data, deep feature attention mechanisms, and dynamic monitoring analysis. This method achieves high-precision, robust, and scalable intelligent remote sensing monitoring of grassland resources. Furthermore, through end-to-end integration of multimodal features and network structure optimization, the distinction between different categories such as grassland, bare land, and sandy land in terms of spatial boundaries and texture structure becomes clearer, significantly improving accuracy compared to existing methods. Moreover, the monitoring data can be imported into GIS, mapping, or reporting systems to assist users in practical work such as grassland ecological monitoring, degradation early warning, and governance assessment. This method can be extended to other types of land cover monitoring, ecological environment analysis, and even large-scale national land resource dynamic surveys. Attached Figure Description

[0082] Figure 1 This is a flowchart illustrating an embodiment of the present invention;

[0083] Figure 2 This is a schematic diagram illustrating the ruggedness extraction principle of the present invention;

[0084] Figure 3 This is a diagram of the SE-ResUNet network architecture of the present invention;

[0085] Figure 4 This is a diagram of the UNet network architecture of the present invention. Detailed Implementation

[0086] The present invention will be further described below with reference to the accompanying drawings and embodiments, which illustrate the above and other technical features and advantages of the present invention. However, the following embodiments are merely preferred embodiments of the present invention and are not exhaustive.

[0087] Example:

[0088] like Figures 1-4 As shown, this invention provides a method for intelligent identification and spatiotemporal monitoring of grassland types for grassland ecology. The method includes the following steps:

[0089] S1: Collect multi-temporal Landsat 8 remote sensing image data of the study area during the growing season from June to August each year when there is no cloud or the cloud cover is extremely low, and collect global digital elevation model data. Number the collected remote sensing image data and record metadata, and select representative labeled samples.

[0090] S2: Perform radiometric calibration, image registration and atmospheric correction preprocessing on the acquired raw remote sensing image data, convert the raw remote sensing image data values ​​into actual surface reflectance, eliminate atmospheric and instrument interference, then cut the silhouette image with the vector boundary of the study area as the boundary, remove irrelevant areas, unify the resolution, and perform cloud removal processing.

[0091] S3: Extract spectral, surface ruggedness and texture features, normalize and dimensionlessly process all extracted features, and then spatially register them to form a multi-channel feature image group to prepare for subsequent model input.

[0092] S4: Establish a segmentation model and train the model. Use the UNet encoder-decoder network structure as the backbone of the segmentation model. Optimize the structure of the features of the study area within the resolution range of the remote sensing image. In the decoding stage, introduce a dynamic sampling mechanism to restore features. After restoration, perform feature fusion.

[0093] S5: After SE-ResUNet is trained, the optimized SE-ResUNet deep segmentation model is used to perform end-to-end automated inference and segmentation of multi-temporal and multi-modal remote sensing image data of the entire region, while monitoring and analyzing the dynamic changes of grassland.

[0094] In areas with strong geomorphic heterogeneity, the synergy of topography, texture, and spectrum can significantly reduce class confusion (such as low-coverage grassland and sparse shrubland, sandy boundaries), enhance feature restoration and boundary discrimination, and improve boundary localization accuracy; SE attention enhances automatic attention to key channels, multi-layer residuals alleviate gradient vanishing and achieve steady-state convergence, dynamic sampling improves upsampling detail reconstruction, and comprehensively achieves higher overall accuracy; the automated generation of annual alignment and transition matrix improves business availability and spatiotemporal analysis capabilities.

[0095] In this embodiment, step S1 specifically includes the following:

[0096] S11: The recorded information includes the acquisition date, center coordinates, image range, sensor parameters, and cloud cover.

[0097] S12: Select imagery or elevation data that strictly overlaps across all multi-temporal and multi-source data in a given period as the spatial benchmark and designate it as the representative annotation sample. The specific selection criteria and parameter requirements for the annotation sample are as follows:

[0098] Category coverage: Covers all target categories such as "high-coverage grassland, medium-coverage grassland, low-coverage grassland, bare land, and sandy land", and the number of samples in each category is balanced (no less than N=1000 pixels or ≥50 ROI patches in each category, adjusted according to the area and heterogeneity of the study area);

[0099] Spatiotemporal consistency: Strictly overlaps with the selected "spatial reference period", and the sampled ROI is free of clouds or cloud shadows on all annual images;

[0100] Geomorphological representativeness: stratified sampling covers different altitudes, slopes, aspects, topographic relief zones, and major parent material / soil types;

[0101] Consistency of time phase: Try to choose the peak period in the middle to late part of the growing season from June to August and pass over on the same day / near the same day to reduce phenological differences;

[0102] Quality control and rejection: The proportion of outlier pixels in the sample is <5%, and cloud / thin cloud / shadow pixels are rejected by quality band or Fmask.

[0103] In this embodiment, the following processes are performed sequentially in step S2:

[0104] S21: Metadata in the MTL file accompanying the raw Landsat 8 remote sensing image data. The numerical values ​​DN in the metadata are converted into physical emissivity and surface reflectivity. The calibration formula is as follows:

[0105] L λ =M L ·Q cal +Α L

[0106] Among them, L λ M represents the calibrated emissivity. L Q represents the band calibration coefficient. cal Represents the original DN value, A L Indicates the addition correction term.

[0107] The metadata is specifically in the Landsat8MTL file:

[0108] Absolute scaling coefficients: RADIANCE_MULT_BAND_x, RADIANCE_ADD_BAND_x, or

[0109] REFLECTANCE_MULT_BAND_x, REFLECTANCE_ADD_BAND_x

[0110] Solar altitude angle / zenith angle: SUN_ELEVATION / SUN_AZIMUTH

[0111] Quantization bit depth: QUANTIZE_CAL_MAX_BAND_x / QUANTIZE_CAL_MIN_BAND_x

[0112] Image acquisition time, row and column number, sensor information, cloud cover, and quality assessment indicator (BQA).

[0113] This achieves the conversion of DN → radiance → apparent / surface reflectance;

[0114] S22, Atmospheric Correction:

[0115] The top-level radiance is corrected to the actual surface reflectance using a rapid atmospheric correction method. Remote sensing software is used as the operating platform to minimize signal confusion caused by atmospheric scattering and absorption. The operating platform can be ENVI, ArcGIS, or SNAP.

[0116] S23, Image Registration:

[0117] Based on DEM elevation data or regional high-precision maps, geometric correction is performed on Landsat 8 remote sensing image data from each period. Affine or polynomial transformation is used during registration. 30 or more uniformly distributed ground control points are selected, and the registration accuracy is controlled within 1 pixel (30 meters). The output is uniformly converted to the WGS-84UTM projection system.

[0118] S24, Image Cropping:

[0119] Based on the boundary vector of the study area in shp format, all images that have undergone the above preprocessing are spatially cropped to ensure that the spatial range of multi-period, multi-source data within the analysis range is highly consistent.

[0120] S25, Resampling:

[0121] Multi-source data are uniformly resampled using the nearest neighbor method or bilinear interpolation algorithm, and all are standardized to a 30-meter resolution to ensure that each data channel pixel corresponds one-to-one and the corresponding latitude and longitude coordinates are completely consistent.

[0122] S26, Cloud-free processing:

[0123] The Fmask algorithm was used to detect and mask clouds or shadows in Landsat 8 images. Detection was performed using thresholding in the visible light and shortwave infrared bands. The Fmask rule-thresholding process, based on a combination of multispectral and thermal infrared methods, works as follows: First, clouds are identified using high reflectivity of visible light and SWIR, and low brightness temperature of TIRS. Then, cloud shadows are identified using solar geometric projection and NIR / SWIR darkening features. Finally, morphological processing is used to form a stable cloud or shadow mask, thus generating a cloud mask layer. Cloud and shadow areas are set to invalid values ​​and removed during subsequent feature statistics. After removal, holes are filled using other temporal images from the same year. The specific method for hole filling is as follows:

[0124] (1) Prioritize the synthesis of neighboring time phases in the same year: For pixels obscured by clouds / shadows, perform time-series synthesis of clear pixels in different time phases of the same year (such as taking the median / quantile or the nearest date value based on the phenological window);

[0125] (2) If no value is available for the same year, the value is backfilled with the average / time series fitted value of the same month over multiple years, and the quality indicator is marked.

[0126] (3) Pixels are filled independently by band, and the original unshaded pixels are preserved to avoid introducing system bias across years.

[0127] In this embodiment, step S3 specifically includes the following:

[0128] S31: Extracting surface vegetation cover areas from images using spectral features:

[0129] Information from each band of Landsat 8 was extracted, and the Enhanced Vegetation Index (EVI) for the study area was calculated using ENVI software. The specific calculation formula is as follows:

[0130]

[0131] Where, ρ NIR ρ R and ρ Blue These represent the surface reflectance in the near-infrared, red, and blue light bands, respectively.

[0132] By introducing the blue light band and utilizing its sensitivity to atmospheric aerosols, atmospheric effects are corrected. At the same time, a gain coefficient, typically 2.5, is added to improve the sensitivity of vegetation cover images, especially in areas with high vegetation cover.

[0133] S32: Extracting surface ruggedness features from the image:

[0134] like Figure 2 As shown, a center point is selected, and the elevation range and average slope of the area near the center point are extracted. Pixel descriptions are introduced to quantitatively express the influence of topographic conditions on vegetation distribution and type classification, improving the classification and identification ability under heterogeneous topography. The extraction range is: statistical analysis is performed within a multi-scale window centered on the target pixel, such as r∈{3×3,5×5,7×7,11×11}, with the unit being pixels, using 5×5 as the main scale. After parallel calculation of multi-scale features, the data is stitched together. The specific calculation formula is as follows:

[0135] The actual coverage area S of the target object on the ground plane Surface The calculation formula is:

[0136]

[0137] The area S of the object projected onto the ground along the sun's direction or a specified direction of the sensor's line of sight. Projection The calculation formula is:

[0138] S Projection =(d-1) 2 ·r

[0139] The area S of the pixel in the i-th row and j-th column ij The calculation formula is:

[0140] S ij =S ΔABD +S ΔACD

[0141]

[0142] Among them, S Surface S represents the land surface area of ​​the extraction window T, used to characterize the "true size" of the object. In ecological / resource accounting, it is used for area statistics, land cover ratio, patch size, area change, etc. ij S represents the area of ​​the (i,j)th cell, used to convert cell counts into actual areas, or for weighted accumulation in partition statistics. ProjectionThe surface projection area of ​​the extraction window T represents the area used for geometric calculation of cloud and mountain shadows, radiation or energy illumination analysis, and direction-dependent visibility or occlusion assessment. (i,j) represents the cell row and column index, serving as the coordinate anchor point for all cell-level quantities to ensure channel alignment and cell-by-cell calculation. d represents the extraction window parameters, and r represents the data resolution, used to determine the spatial scale of local statistics such as ruggedness and texture, reflecting the multi-scale representation capability from detail to coarse scale. ΔABD and S ΔABD P1 and P2 represent the areas of the triangles formed by points A, B, and D, and A, C, and D, respectively. These are used to decompose complex polygons into triangles for accurate calculation of area, centroid, and centroid moments, or as partitioned areas in local geometric derivations. P1 and P2 represent weighting coefficients used for linear weighted fusion of multiple quantities. AB l BD l AD l AC l CD and l XY Both represent the Euclidean distance between two points, providing a basis for indicators such as shape index, fragmentation, and boundary complexity. XY h AB h BD h AC h CD and h AD Both represent the relative degree between two points, used to characterize terrain undulation, slope / volume estimation, cloud shadow geometry, and can also serve as a basic quantity for ruggedness, t 11 The self-maintenance of the representative period t1, t 12 This represents a transition from t1 to t2, where t 21 This represents a rotation from t2 to t1, where t 22 The self-maintained state of period t2 is used to support the K×K dimensional transition matrix M∈R. K×K The construction and interpretation of these data facilitate the identification of statistical area transfer, net increase / loss, degradation / recovery, and the generation of annual change maps.

[0143] S33: Extract image texture features:

[0144] Statistical methods are used to extract image texture features at different window scales. The extraction method can simultaneously employ the gray-level co-occurrence matrix method and the local variance method to balance roughness and directionality representation; when computational resources are limited or for comparative experiments, either method can be chosen. Specific features required include mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment of angle, and correlation, further reflecting the spatial structural differences of grassland, bare land, shrubs, and other land features. The specific calculation formulas are as follows:

[0145] The mean represents the expectation of gray levels based on the edge distribution of the GLCM (Gray Co-occurrence Matrix). The specific calculation formula is as follows:

[0146] Mean = ∑ i ∑ j p(i,j)*i;

[0147] Variance represents the degree of dispersion of gray levels relative to the mean, and the specific calculation formula is as follows:

[0148] Variance=∑ i Σ j p(i,j)*(i-Mean) 2 ;

[0149] Homogeneity indicates that the closer the gray levels of neighboring areas are (i.e., the smaller the difference), the greater the weight. A higher value indicates a smoother texture. The specific calculation formula is as follows:

[0150]

[0151] Contrast emphasizes the difference in gray levels and edges. The square of the difference magnifies the overall difference. The specific calculation formula is as follows:

[0152] Contrast=Σ i Σ j p(i,j)*(ij) 2 ;

[0153] Dissimilarity represents a linear measure of the frequency of gray-level differences. Compared to contrast, it amplifies large differences more gently. The specific calculation formula is as follows:

[0154] Dissimilarity=Σ i Σ j p(i,j)*|ij|;

[0155] Entropy represents the randomness or disorder of a texture, and its specific calculation formula is as follows:

[0156] Entropy = Σ i Σ j p(i,j)*lnΡ(i,j);

[0157] The second moment of the angle (ASM) represents the regularity or repetition of the texture. The specific calculation formula is as follows:

[0158] ASM=∑ i ∑ j p(i,j) 2 ;

[0159] Correlation measures the degree of synchronous shift of a pair of adjacent gray levels (i,j) relative to their respective means, and normalizes it to a comparable range of approximately [-1,1] using variance. The specific calculation formula is as follows:

[0160]

[0161] Where (i,j) represents the index of the quantized gray level, with a value range of {0,1,…,G-1}, G represents the total number of quantized levels, and P(i,j) represents the joint probability. P(i,j) is a normalized GLCM element calculated under a given window, orientation, and pixel spacing, representing the joint probability of gray level pair (i,j) occurring, satisfying ∑ i ∑ j P(i,j), where i represents the gray level value in the row direction;

[0162] Mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment of angle, and correlation are all calculated from each input channel of the remote sensing image (such as a single band or its derived index) at different window scales. Among them, the mean and variance are directly based on the gray-level statistics of pixels within the window, and the remaining features are calculated based on the gray-level co-occurrence matrix (GLCM) of the window. Multi-scale textures are obtained by changing the window size (such as 3×3, 5×5, 7×7, etc.).

[0163] The functions of each feature are as follows:

[0164] The mean reflects the local brightness level, while the variance describes the local contrast and fluctuations.

[0165] Homogeneity measures the similarity of gray levels in neighboring areas; a high value indicates smooth texture. Contrast emphasizes strong edges and bands.

[0166] Dissimilarity emphasizes the frequent occurrence of grayscale differences; entropy measures randomness and disorder, with high values ​​usually corresponding to complex textures or noise.

[0167] The angular second moment (ASM / energy) represents the regularity / uniformity of the texture; a high value indicates high structural repeatability and low noise.

[0168] Correlation characterizes the linear dependence and directionality of neighborhood gray levels.

[0169] The above features can improve the separability of land cover categories in the texture dimension, and are particularly effective in distinguishing and delineating easily confused categories such as woodland / grassland, bare land / arable land, and urban construction / roads.

[0170] In this embodiment, step S4 specifically includes the following:

[0171] Step S41: In each encoding and decoding module of the UNet encoder-decoder network structure, embed the SE channel attention mechanism to extract global information and relabel the features. During deep learning, this enables the segmentation model to automatically focus on the channel layer that is most representative of the land surface type.

[0172] UNet's structure resembles the letter U, consisting of an encoder (shrinking path) and a decoder (expanding path), with skip connections fusing multi-scale features. The encoder borrows from the concept of Convolutional Neural Networks (CNNs), using a series of convolution and pooling operations to gradually extract image features. The encoder comprises multiple convolutional modules, each consisting of two 3×3 convolutional layers. Following each convolutional layer is a ReLU activation function, enhancing the model's non-linear expressive power and enabling it to learn more complex image features. After each convolutional block, a 2×2 max pooling layer with a stride of 2 is applied to reduce the spatial resolution of the image while expanding the receptive field, allowing the model to capture more comprehensive image information. As the network depth increases, the size of the feature maps gradually decreases, while the number of channels gradually increases. Therefore, the model can learn more representative image features at lower resolutions.

[0173] The decoder section has a symmetrical structure with the encoder. It uses deconvolution and upsampling to restore the low-resolution feature map to the original image size to achieve the image segmentation prediction task. The decoder also includes multiple deconvolution modules. Each deconvolution module consists of a deconvolution layer with a stride of 2 and a size of 2×2, and two 3×3 convolution layers. After the convolution layers, a ReLU activation function is connected to enhance the nonlinear mapping capability of the model. This allows the decoder to effectively extract and utilize feature information during the process of restoring the image size and accurately complete the image segmentation prediction. The deconvolution layer is responsible for upsampling the feature map to gradually increase its size, while the convolution layers are used to further fuse features and refine the segmentation results.

[0174] In this way, the decoder gradually recovers the detailed information of the image and finally outputs a segmentation mask with the same size as the input image.

[0175] The UNet network concatenates feature maps from different levels in the encoder with corresponding feature maps from the decoder, thereby fusing feature information at different scales. Specifically, the output feature map of a certain layer in the encoder is directly connected to a network layer in the decoder with the same or similar resolution. For example, in a certain layer of the decoder, the feature map obtained from the deconvolution of that layer is concatenated with the feature map from the encoder after the same number of pooling operations along the channel dimension, and then subsequent convolution operations are performed. This allows the decoder to not only utilize the features obtained from its own upsampling and convolution operations when recovering image details, but also to leverage the richer, lower-level feature information extracted from the corresponding layer in the encoder. This enables it to better capture edge positions and small objects in the image, improving the accuracy and precision of segmentation. The structure of UNet is as follows: Figure 4 As shown;

[0176] Among them, the SE channel attention mechanism obtains the channel description through global average pooling, and generates the channel weight vector through a two-layer fully connected network of compression and activation. This is used to scale the feature channels one by one, realize the recalibration of the features, that is, change the relative weight of each channel in the subsequent convolution.

[0177] SE specific operation steps and the meaning of recalibration:

[0178] For the input feature map X∈R ^ {H×W×C}:

[0179] Where R represents that the values ​​of each element of the tensor belong to the set of real numbers, that is, X is a real tensor with shape H×W×C; C represents the number of channels, such as the number of channels after stitching multi-band, EVI, DEM, ruggedness and texture, etc., and each extracted one-dimensional feature is counted as one channel; H represents the height of the input feature image; W represents the width of the input feature image;

[0180] Squeeze: Uses global average pooling, zc = (1 / HW)∑h∑wX(h,w,c);

[0181] Excitation: Dimensionality reduction of the fully connected layer (C→C / r, r is usually 8 or 16) + ReLU, then dimension increase of the fully connected layer C / r→C + Sigmoid, to obtain the channel weight vector s∈(0,1)^C;

[0182] The recalibration is: Y(h,w,c)=s(c)·X(h,w,c), which means scaling each channel by channel with s(c), which is essentially changing the importance of the channels by redistributing the weights.

[0183] Step S42: Introduce multi-layer residual connection modules into each level of feature pathways in the segmentation model. The specific model structure parameters are as follows:

[0184] The convolution kernel size is uniformly 3×3 with a stride of 1. The activation function is ReLU. The entire encoder-decoder network is set to a four-level structure. The specific number of layers and other parameters can be flexibly adjusted according to the size of the input data space. This alleviates the gradient vanishing and information decay that occur during the forward and backward propagation of deep features, while accelerating the model training convergence speed, improving network stability, and thus improving the final segmentation accuracy.

[0185] The multi-layer residual connection module consists of 2–3 stacked standard residual blocks. Each standard residual block contains “3×3 convolution-BN-ReLU-3×3 convolution-BN” plus an identity shorting. If necessary, 1×1 convolution is used to match the channels. The module is deployed in each level of feature pathway to alleviate gradient vanishing and improve training stability and expressive power.

[0186] Identity shorting means that the output y = F(x) + x is the input x is directly added to the output of the convolutional sub-network F(x) through a jump connection. This achieves a fast path for identity mapping, which facilitates gradient backpropagation, alleviates gradient vanishing, and stabilizes training.

[0187] When the number of input and output channels or spatial dimensions do not match, the projection shorting of a 1×1 convolution is used to replace the identity to match the dimensions.

[0188] Step S43: Align all preprocessed and normalized features pixel by pixel in space, and fuse Landsat8 spectral bands, EVI, DEM elevation, surface ruggedness, image textures and statistical features at all levels in a channel stitching manner. Channel stitching connects multi-source, multi-scale or multi-type features in the channel dimension as input or intermediate layer features of the neural network. After fusion, the features are stitched into different channels of the input layer to form a multi-channel input tensor.

[0189] During pixel alignment: a unified spatial reference and pixel size are used to strictly register each data source. The same anchor point grid is used to resample data of different resolutions to a unified grid. This ensures that the channel values ​​of each pixel correspond to the same row and column index positions, thereby achieving pixel-by-pixel alignment.

[0190] Step S44: Combining known field sampling data, existing grassland or land cover thematic databases, and expert manual interpretation results, construct a high-precision supervised classification annotation sample set. The annotation categories cover high-coverage grassland, medium-coverage grassland, low-coverage grassland, bare land, sandy land, and other major land surface types. All samples strictly follow unified interpretation standards and spatial consistency. The annotation samples are divided into training set and validation set in an 8:2 ratio. The number of samples in each category is controlled between 1,000 and 10,000 typical patches or pixels, covering the distribution characteristics of major landforms, climates, and vegetation types. This ensures both the sufficiency of model training and meets the needs of generalization evaluation.

[0191] In sample construction, a unified interpretation standard and spatial consistency requirements are implemented: First, based on a closed and mutually exclusive land surface category system (such as cultivated land, forest land, grassland, water bodies, construction land, and bare land) and clear inclusion / exclusion rules, combined with multi-source evidence and fixed-caliber thresholds (such as NDVI, MNDWI, NDBI, and DEM slope / elevation based on contemporaneous windowing), manual / semi-automatic interpretation is performed, and easily confused categories are handled according to the minimum mapping unit and neighborhood consistency rules; Second, all data are unified to the same projection coordinate system and resolution to ensure that the geometric registration error is less than 0.5 pixels. Samples are assigned values ​​using the pixel center or coverage threshold method and are independently interpreted by two people + arbitration review to eliminate anomalies caused by cloud / snow / shadow and temporal inconsistencies, ensuring that the labels correspond one-to-one in space and are consistent in time sequence, thereby forming a reproducible high-precision supervised classification label sample set.

[0192] Step S45: A dynamic sampling mechanism is introduced in the decoding stage to enhance feature restoration and boundary discrimination, and achieve accurate feature restoration. Through the skipconnection mechanism, the shallow spatial details at the encoding end are directly fused with the high-level semantic features at the decoding end.

[0193] The dynamic sampling mechanism works as follows: First, the decoder obtains a preliminary category probability map and calculates uncertainty or importance indices, such as pixel entropy U = Σ. k P k ·logP k Based on the boundary response or residual, sampling points with high uncertainty and boundary neighborhoods are adaptively selected according to the quantile threshold or a fixed ratio. Then, at the sampling point location, a continuous sampling grid is constructed based on the sampling offset of the lightweight prediction. Features are resampled from the multi-scale feature / original image using methods such as bilinear interpolation, and locally refined by a small decoding head (1×1 / 3×3 convolution + normalization + activation) to update the class probability at the corresponding position. During the training phase, difficult samples are given higher weights or repeated sampling. During the inference phase, it can be iterated once or multiple times until the confidence level is raised to the threshold or the iteration limit is reached. This mechanism is executed in parallel / alternately with conventional global decoding, which significantly improves the classification accuracy of boundary and easily confused regions with a slight increase in computation.

[0194] It fully preserves multi-scale information such as spatial texture, boundaries and land cover morphology in the input images, improves the model's ability to identify and segment the continuity of small patches, narrow grasslands and edge ecological types, and ultimately significantly improves the overall accuracy and practical application value of grassland cover type segmentation.

[0195] In this embodiment, the specific automated reasoning and segmentation process in step S5 includes the following:

[0196] Step S51: Input the multi-channel remote sensing and auxiliary terrain feature images, which are integrated according to time sequence and spatial range, into the trained segmentation model in the form of standardized three-dimensional tensors. The specific features include Landsat8 multiband, EVI, DEM elevation, surface ruggedness and texture features. Align and fuse the Landsat8 multiband, EVI, DEM elevation, surface ruggedness and texture feature channels by pixel. After inference, align by year to generate "transition matrix + change spatial distribution map".

[0197] The standardized three-dimensional tensor form (H×W×C) is a common data organization method for deep learning image segmentation.

[0198] Input tensor size: H×W×C, where C is the number of channels (e.g., the number of channels after stitching together multi-band, EVI, DEM, ruggedness, and texture, etc.), and each channel is normalized to [0,1] or standardized to zero mean and unit variance;

[0199] The specific features correspond to the following:

[0200] Landsat-8 multiband: corresponds to OLI / SR reflectivity channels

[0201] Common combinations: B2 (Blue), B3 (Green), B4 (Red), B5 (NIR), B6 ​​(SWIR1), B7 (SWIR2)

[0202] Optional: B1 (Coastal), B9 (Cirrus), B10 / 11 (TIRS brightness temperature) will be treated as independent channels if they are used in modeling;

[0203] EVI (Enhanced Vegetation Index): Corresponds to the vegetation index channel.

[0204] Calculation parameters: EVI = 2.5 * (NIR - Red) / (NIR + 6 * Red - 7.5 * Blue + 1),

[0205] Where NIR = B5, Red = B4, Blue = B2 (or the corresponding band of the sensor);

[0206] DEM elevation: corresponds to the absolute elevation channel.

[0207] Parameter source: Digital elevation model (unit: m), resampled at the same resolution / projection as the image.

[0208] Surface ruggedness: corresponds to the intensity of terrain undulation.

[0209] Texture features: Corresponding to a multi-channel parameter set based on grayscale statistics / GLCM.

[0210] Step S52: The segmentation model automatically identifies the spatial feature combination of each input pixel, and extracts, fuses and discriminates features layer by layer through the backbone network and attention mechanism, outputting a pixel-level grassland classification result consistent with the original input space. Pixel-level grassland classification refers to the category determination result of each raster pixel, that is, the output is a classification map with the same spatial size as the input raster and assigned values ​​pixel by pixel.

[0211] Step S53: Spatially align all processed annual category raster layer data with the initial image in chronological order to continuously reconstruct the spatial change process of grassland types. Then, perform automated change detection of grassland cover types between years. The specific automated implementation method is as follows:

[0212] The classification results for each year are resampled / reprojected onto a grid with the same spatial reference, resolution, range, and origin to ensure that year t and t+1 are completely consistent in the raster rows and columns. A "baseline period" raster is fixed as a template in the production chain, and alignment transformation is performed on other years.

[0213] The specific implementation method of automated monitoring is as follows:

[0214] Pixel-by-pixel comparison of the category codes for year t and t+1 generates a pixel-level transition coding map (e.g., code=C). t ×K+C t+1 (K is the total number of categories), and accumulates to the full time series. This automatically calculates the transfer count and area.

[0215] In this embodiment, step S53, which involves automatically detecting changes in grassland cover type between years, specifically includes the following:

[0216] Step S531: Using the classified annual raster layer data, statistically analyze the changes in grassland type for the same pixel across different years:

[0217] Step S532: Construct a grassland type transfer matrix to quantify the area increase, loss, and transformation relationships among different grassland types such as high-coverage grassland, medium-coverage grassland, low-coverage grassland, bare land, and sandy land, achieving accurate measurement of annual changes in area by type. The transfer matrix quantification process is as follows:

[0218] Let the set of categories be L = {1...K}, and the transition matrix be M ∈ R. ^ {K×K}, element M ij The pixel count for transitioning from category i (year t) to category j (year t+1); area quantization A ij =M ij ×A px , where A pxThis is the area of ​​a single pixel (e.g., 30m × 30m); further calculations can be made of indicators such as annual net increase / loss, conversion rate, and stability.

[0219] Step S533: Automatically generate a spatial distribution map of grassland type changes. The specific process for generating the spatial distribution map is as follows:

[0220] Based on the pixel-level transfer coding map, color palette rendering is set according to the type of transfer of interest (e.g., high cover → medium cover, grassland → bare land / sand); or rules are used to classify "degradation / restoration / stabilization", for example:

[0221] Degradation: Coverage level declines or turns into bare land / sand.

[0222] Restoration: Coverage level increases or bare / sandy land transforms into grassland.

[0223] Stable: Category remains unchanged;

[0224] This visually demonstrates the degradation path of high-cover grasslands to low-cover or bare land and sandy land, as well as the restoration and evolution process of some areas.

[0225] In practical applications, the model can process large-scale remote sensing images in batches, enabling fully automatic segmentation of grassland cover in the entire study area by year, season, or multiple periods without human intervention. Each pixel is classified into a certain category, such as high-coverage grassland, medium-coverage grassland, low-coverage grassland, bare land, sandy land, etc., achieving efficient surveys and dynamic tracking of long-term grassland evolution trends, providing timely and high-precision spatial information support for grassland resource management and ecological status monitoring.

[0226] The above description is merely a preferred embodiment of the present invention and is illustrative rather than restrictive. Those skilled in the art will understand that many changes, modifications, and even equivalents can be made within the spirit and scope defined by the claims of the present invention, all of which will fall within the protection scope of the present invention.

Claims

1. A method for intelligent identification and spatiotemporal monitoring of grassland types for grassland ecology, characterized in that, The intelligent grassland type identification and spatiotemporal monitoring method includes the following steps: S1: Collect multi-temporal Landsat 8 remote sensing image data of the study area during the growing season each year when there is no cloud or the cloud cover is extremely low, and collect global digital elevation model data. Number the collected remote sensing image data and record the metadata, and select representative labeled samples. S2: Perform radiometric calibration, image registration and atmospheric correction preprocessing on the acquired raw remote sensing image data, convert the raw remote sensing image data values ​​into actual surface reflectance, eliminate atmospheric and instrument interference, then cut the silhouette image with the vector boundary of the study area as the boundary, remove irrelevant areas, unify the resolution, and perform cloud removal processing. S3: Extract spectral, surface ruggedness and texture features, normalize and dimensionlessly process all extracted features, and then spatially register them to form a multi-channel feature image group to prepare for subsequent model input. S4: Establish a segmentation model and train the model. Use the UNet encoder-decoder network structure as the backbone of the segmentation model. Optimize the structure of the features of the study area within the resolution range of the remote sensing image. In the decoding stage, introduce a dynamic sampling mechanism to restore features. After restoration, perform feature fusion. S5: After SE-ResUNet is trained, the optimized SE-ResUNet deep segmentation model is used to perform end-to-end automated inference and segmentation of multi-temporal and multi-modal remote sensing image data of the entire region, while monitoring and analyzing the dynamic changes of grassland.

2. The method for intelligent identification and spatiotemporal monitoring of grassland types for grassland ecology as described in claim 1, characterized in that, Step S1 specifically includes the following: S11: The recorded information includes the acquisition date, center coordinates, image range, sensor parameters, and cloud cover. S12: Select images or elevation data that strictly overlap with all multi-temporal and multi-source data in a period as the spatial reference and designate them as representative annotation samples.

3. The method for intelligent identification and spatiotemporal monitoring of grassland types for grassland ecology as described in claim 1, characterized in that, The following processes are performed sequentially in step S2: S21, Radiation Calibration: The metadata in the MTL file accompanying the raw Landsat 8 remote sensing image data is used to convert the numerical values ​​DN in the metadata into physical emissivity and surface reflectivity. The calibration formula is as follows: L λ =M L ·Q cal +A L Among them, L λ M represents the calibrated emissivity. L Q represents the band calibration coefficient. cal Represents the original DN value, A L Indicates the addition correction term; S22, Atmospheric Correction: By using a rapid atmospheric correction method, the top-level radiance is corrected to the actual surface reflectance. Using remote sensing software as the operating platform, signal confusion caused by atmospheric scattering and absorption is minimized. S23, Image Registration: Based on DEM elevation data or regional high-precision maps, geometric correction is performed on Landsat 8 remote sensing image data from each period. 30 or more evenly distributed ground control points are selected, and the registration accuracy is controlled within 1 pixel. The data is then output to the WGS-84UTM projection system. S24, Image Cropping: Based on the boundary vector of the study area in shp format, spatial cropping is performed on all images that have undergone the above preprocessing. S25, Resampling: Multi-source data are spatially resampled and standardized to a 30-meter resolution. S26, Cloud-free processing: The Fmask algorithm is used to detect and mask clouds or shadows in Landsat 8 images. Detection is performed using thresholds in the visible light and shortwave infrared bands. The Fmask rule-threshold process based on the combination of multispectral and thermal infrared is as follows: First, clouds are identified using the high reflectivity of visible light and SWIR and the low brightness temperature of TIRS. Then, cloud shadows are identified using solar geometric projection and NIR / SWIR darkening features. Finally, a stable cloud or shadow mask is formed through morphological processing, thereby generating a cloud mask layer. Cloud and cloud shadow areas are set to invalid values ​​and removed during subsequent feature statistics. After removal, holes are filled using other temporal images from the same year.

4. The method for intelligent identification and spatiotemporal monitoring of grassland types for grassland ecology as described in claim 1, characterized in that, Step S3 specifically includes the following: S31: Extracting surface vegetation cover areas from images using spectral features: Information from each band of Landsat 8 was extracted, and the Enhanced Vegetation Index (EVI) for the study area was calculated using ENVI software. The specific calculation formula is as follows: Where, ρ NIR ρ R and ρ Blue These represent the surface reflectance in the near-infrared, red, and blue light bands, respectively. By introducing the blue light band and utilizing its sensitivity to atmospheric aerosols, atmospheric effects are corrected by adding a gain coefficient, typically 2.

5. S32: Select a center point, extract the elevation range and average slope of the area near the center point, and introduce pixel descriptions to quantitatively express the influence of terrain conditions on vegetation distribution and type classification, thereby improving the classification and recognition ability under heterogeneous terrain. The extraction range is: with the target pixel as the center, statistics are performed within a multi-scale window, such as r∈{3×3, 5×5, 7×7, 11×11}, with 5×5 as the main scale. After parallel calculation of multi-scale features, the data is then stitched together. The specific calculation formula is as follows: The actual coverage area S of the target object on the ground plane Surface The calculation formula is: The area S of the object projected onto the ground along the sun's direction or a specified direction of the sensor's line of sight. Projection The calculation formula is: S Projection =(d-1) 2 ·r The area S of the pixel in the i-th row and j-th column ij The calculation formula is: S ij =S ΔABD +S ΔACD Among them, S Surface S represents the surface area of ​​the extracted window T. ij S represents the area of ​​the (i,j)th pixel. Projection Let (i,j) represent the surface projection area of ​​the extraction window T, (i,j) represent the cell row and column indices, d represent the extraction window parameters, r represent the data resolution, and S represent the surface projection area of ​​the extraction window T. ΔABD and S ΔABD P1 and P2 represent the areas of the triangles formed by points A, B, D, and A, C, D, respectively. P1 and P2 represent weighting coefficients used for linear weighted fusion of multiple quantities. AB l BD l AD l AC l CD and l XY Both represent the Euclidean distance between the two points, h XY h AB h BD h AC h CD and h AD Both represent the relative degree between two points, t 11 The self-maintenance of the representative period t1, t 12 This represents a transition from t1 to t2, where t 21 This represents a rotation from t2 to t1, where t 22 The self-maintained state of the representative period t2; S33: Extract image texture features: Statistical methods were used to extract image texture features at different window scales. The required features specifically included mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment of angle, and correlation. The specific calculation formulas are as follows: The mean is: Mean=∑ i ∑ j p(i,j)*i; The variance is: Variance=∑ i Σ j p(i,j)*(i-Mean) 2 ; Homogeneity is: Contrast is: Contrast=∑ i ∑ j p(i,j)*(i-j) 2 ; Dissimilarity is: Dissimilarity=∑ i ∑ j p(i,j)*|i-j|; Entropy is: Entropy=∑ i ∑ j p(i,j)*lnΡ(i,j); The second moment of the angle, ASM, is: ASM=∑ i ∑ j p(i,j) 2 ; The correlation is: Where (i,j) represents the index of the quantized gray level, with a value range of {0, 1, ..., G-1}, G represents the total number of quantized levels, and P(i,j) represents the joint probability. P(i,j) is a normalized GLCM element calculated under a given window, orientation, and pixel spacing, representing the joint probability of gray level pair (i,j) occurring, satisfying ∑ i ∑ j P(i,j), where i represents the gray level value in the row direction.

5. The method for intelligent identification and spatiotemporal monitoring of grassland types for grassland ecology as described in claim 1, characterized in that, Step S4 specifically includes the following: Step S41: In each encoding and decoding module of the UNet encoder-decoder network structure, embed the SE channel attention mechanism to extract global information and recalibrate the features, so that the segmentation model automatically focuses on the channel layer that is most representative of the land surface type distinction. Step S42: Introduce multi-layer residual connection modules into the feature pathways at each level in the segmentation model to alleviate gradient vanishing and information decay that occur during the forward and backward propagation of deep features; Step S43: Align all preprocessed and normalized features pixel by pixel in space, and fuse Landsat8 spectral bands, EVI, DEM elevation, surface ruggedness, image textures and statistical features at all levels in a channel stitching manner. Channel stitching connects multi-source, multi-scale or multi-type features in the channel dimension as input or intermediate layer features of the neural network. After fusion, the features are stitched into different channels of the input layer to form a multi-channel input tensor. Step S44: Combining known field sampling data, existing grassland or land cover thematic databases, and expert manual interpretation results, construct a high-precision supervised classification and labeling sample set, label the main land surface types, and all samples strictly follow the unified interpretation standards and spatial consistency. Step S45: A dynamic sampling mechanism is introduced in the decoding stage to achieve accurate feature restoration. Through the skipconnection mechanism, the shallow spatial details at the encoding end are directly fused with the high-level semantic features at the decoding end.

6. The method for intelligent identification and spatiotemporal monitoring of grassland types for grassland ecology as described in claim 1, characterized in that, The specific automated reasoning and segmentation process in step S5 includes the following: Step S51: Input the integrated multi-channel remote sensing and auxiliary terrain feature images into the trained segmentation model in the form of a standardized three-dimensional tensor. Specific features include Landsat8 multi-band, EVI, DEM elevation, surface ruggedness and texture features. Step S52: The segmentation model automatically identifies the spatial feature combination of each input cell, and extracts, fuses and discriminates features layer by layer through the backbone network and attention mechanism, and outputs a cell-level grassland classification result that is consistent with the original input space. Cell-level grassland classification refers to the category determination result of each raster cell, that is, the output is a classification map that is consistent with the input raster in spatial size and assigns values ​​to each cell. Step S53: Align all processed annual category raster layer data spatially with the initial image in chronological order to continuously reconstruct the spatial change process of grassland types, and then perform automated change detection on grassland cover types between years.

7. The method for intelligent identification and spatiotemporal monitoring of grassland types for grassland ecology as described in claim 6, characterized in that, The automated change detection of grassland cover type between years in step S53 specifically includes the following: Step S531: Using the classified annual raster layer data, statistically analyze the changes in grassland type for the same pixel across different years: Step S532: Construct a grassland type transfer matrix to quantify the relationship between the increase, loss and transformation of different grassland types, and achieve accurate measurement of the annual changes in the area of ​​different types; Step S533: Automatically generate a spatial distribution map of grassland type changes, which visually shows the degradation path of high-coverage grassland to low-coverage or bare land and sandy land, as well as the restoration and evolution process of some areas.

Images