A Method and System for Grassland Drought Identification and Assessment Based on Multi-Source Remote Sensing Image Data

Abstract

This invention relates to the field of grassland drought management technology, and in particular to a method and system for grassland drought identification and assessment based on multi-source remote sensing image data. The method includes the following steps: acquiring multi-source remote sensing images and digital elevation models (DEMs) covering a target grassland area; performing terrain correction and band-weighted fusion on the multi-source remote sensing images based on the DEMs to obtain a terrain-corrected fused image; identifying and segmenting spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature based on the terrain-corrected fused image; fitting local wet and dry edges to the grassland homogeneous patch images and generating a grassland zoning drought grid; and implementing dynamic livestock carrying capacity management based on the drought intensity values ​​in the grassland zoning drought grid. This invention achieves refined grassland drought assessment by segmenting grassland homogeneous patches and assigning community sensitivity coefficients, thereby dynamically adjusting the livestock carrying capacity of pastoral areas.

Description

Technical Field

[0001] This invention relates to the field of grassland drought management technology, and in particular to a grassland drought identification and assessment method and system based on multi-source remote sensing image data. Background Technology

[0002] Remote sensing technology, with its macroscopic, periodic, and non-destructive characteristics, has become an important tool for monitoring drought in large-scale grasslands. Currently, grassland drought monitoring and assessment mainly rely on single remote sensing data sources or simple vegetation index threshold methods. Using a single data source for grassland drought assessment makes it difficult to comprehensively reflect the complex characteristics of drought. Grassland drought is a complex process involving the combined effects of multiple factors such as precipitation, temperature, soil moisture, and vegetation cover. A single remote sensing data source cannot capture this complexity, leading to biased and inaccurate assessment results. For example, relying solely on vegetation indices cannot effectively distinguish between vegetation degradation caused by drought and vegetation changes caused by other factors (such as overgrazing and insect pests), resulting in biased drought severity assessments. Furthermore, when grasslands experience drought, their carrying capacity declines rapidly. If livestock carrying capacity is not adjusted in time, it will lead to further grassland degradation, creating a vicious cycle and hindering the sustainable use of grassland resources. Summary of the Invention

[0003] Based on this, the present invention provides a grassland drought identification and assessment method and system based on multi-source remote sensing image data to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, a grassland drought identification and assessment method based on multi-source remote sensing image data includes the following steps:

[0005] Step S1: Acquire multi-source remote sensing images and digital elevation models covering the target grassland area, and perform terrain correction and band weighted fusion on the multi-source remote sensing images based on the slope and aspect in the digital elevation model to obtain a terrain-corrected fused image.

[0006] Step S2: Based on topographic correction, identify and segment spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature from the fused images; assign a drought sensitivity coefficient to the community type in the grassland homogeneous patch map, evaluate the drought intensity value of each patch based on the coefficient, and then fit the local wet and dry edges to generate a grassland zoning drought grid.

[0007] Step S3: Determine the dynamic carrying capacity of each grid cell under the current drought stress based on the drought intensity value in the grassland drought grid;

[0008] Step S4: The dynamic carrying capacity limit is compared with the real-time herd distribution data transmitted by the IoT device to determine the real-time carrying capacity. If the real-time herd distribution data exceeds the dynamic carrying capacity limit, it is determined to be an overgrazing state and dynamic carrying capacity management is implemented; otherwise, it is determined to be a safe carrying capacity state and the current grazing strategy is maintained.

[0009] The present invention also provides a grassland drought identification and assessment system based on multi-source remote sensing image data, which executes the grassland drought identification and assessment method based on multi-source remote sensing image data as described above. The grassland drought identification and assessment system based on multi-source remote sensing image data includes:

[0010] The terrain correction and fusion module is used to acquire multi-source remote sensing images and digital elevation models covering the target grassland area, and to perform terrain correction and band-weighted fusion on the multi-source remote sensing images based on the slope and aspect in the digital elevation model to obtain a terrain-corrected fused image.

[0011] The drought zoning assessment module is used to identify and segment spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature based on topographically corrected fused images. It assigns a drought sensitivity coefficient based on the community type in the grassland homogeneous patch map, determines the drought intensity value of each patch based on the coefficient, and then fits the local dry and wet edges to generate a grassland zoning drought grid.

[0012] The carrying capacity calculation module is used to determine the dynamic carrying capacity of each grid cell under the current drought stress based on the drought intensity value in the grassland drought grid.

[0013] The overgrazing monitoring and management module is used to determine the real-time carrying capacity by comparing the dynamic carrying capacity limit with the real-time distribution data of the herd transmitted by the Internet of Things devices. If the real-time distribution data of the herd exceeds the dynamic carrying capacity limit, it is determined to be an overgrazing state and dynamic carrying capacity management is implemented; otherwise, it is determined to be a safe carrying capacity state and the current grazing strategy is maintained.

[0014] The beneficial effects of this invention are as follows:

[0015] On the one hand, by simultaneously acquiring multi-source remote sensing images and digital elevation models, and using topographic factors such as slope and aspect to correct and weight the remote sensing images, this method can effectively eliminate the distortion of remote sensing signals caused by uneven solar radiation due to topographic undulations and differences in shady slopes. Compared with traditional methods that rely on a single, uncorrected data source, the topographically corrected fused image constructed in this invention serves as the basis for analysis, which can more realistically reflect the actual conditions of surface vegetation and temperature, ensuring the physical authenticity and data reliability of subsequent drought identification and assessment from the source, and significantly improving the accuracy of the assessment results.

[0016] On the other hand, by identifying and segmenting "homogeneous grassland patches" with vegetation and surface temperatures similar to each other, and innovatively introducing a "drought sensitivity coefficient" based on community type to assess drought intensity, this method leaps from traditional pixel-level analysis to assessment at the ecological functional unit level. This process can accurately distinguish the intrinsic response differences of different grassland types (such as meadow steppe and desert steppe) to drought stress, effectively solving the problem that existing technologies cannot distinguish between vegetation degradation caused by drought and other factors (such as overgrazing). This makes drought assessment no longer a single physical indicator judgment, but a comprehensive diagnosis of ecosystem vulnerability, greatly improving the scientific nature and pertinence of the assessment results.

[0017] On the other hand, this invention directly transforms the assessed "grassland zoning drought grid" into a quantified "dynamic carrying capacity quota," and links it with the real-time herd distribution data transmitted back by IoT devices to construct a closed-loop management system from remote sensing monitoring to precise control. This system achieves real-time and dynamic matching of grass-livestock relationships, can respond quickly when drought occurs, and avoids overgrazing through early warning and dynamic management instructions, providing strong technical support for the sustainable utilization of grassland resources and refined management of pastoral areas. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the steps in the grassland drought identification and assessment method based on multi-source remote sensing image data of the present invention;

[0019] Figure 2 This is a schematic diagram of the grassland drought identification and assessment system module based on multi-source remote sensing image data of the present invention;

[0020] Figure 3 Example remote sensing image of the target grassland area;

[0021] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0022] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0023] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0024] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated items listed.

[0025] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides a grassland drought identification and assessment method based on multi-source remote sensing image data, comprising the following steps:

[0026] Step S1: Acquire multi-source remote sensing images and digital elevation models covering the target grassland area, and perform terrain correction and band weighted fusion on the multi-source remote sensing images based on the slope and aspect in the digital elevation model to obtain a terrain-corrected fused image.

[0027] Step S2: Based on topographic correction, identify and segment spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature from the fused images; assign a drought sensitivity coefficient to the community type in the grassland homogeneous patch map, evaluate the drought intensity value of each patch based on the coefficient, and then fit the local wet and dry edges to generate a grassland zoning drought grid.

[0028] Step S3: Determine the dynamic carrying capacity of each grid cell under the current drought stress based on the drought intensity value in the grassland drought grid;

[0029] Step S4: The dynamic carrying capacity limit is compared with the real-time herd distribution data transmitted by the IoT device to determine the real-time carrying capacity. If the real-time herd distribution data exceeds the dynamic carrying capacity limit, it is determined to be an overgrazing state and dynamic carrying capacity management is implemented; otherwise, it is determined to be a safe carrying capacity state and the current grazing strategy is maintained.

[0030] In this embodiment of the invention, firstly, a digital elevation model (DEM) of the target grassland area (e.g., a grassland test area) is acquired with a spatial resolution of 30 meters. Simultaneously, Landsat 8OLI / TIRS satellite imagery covering the same area is acquired, selecting cloud-free images containing visible light (blue, green, and red bands), near-infrared bands, and thermal infrared bands.

[0031] In one implementation of this invention, DEM data is used to calculate the slope and aspect angle for each pixel in the image. Assume a pixel has a slope of 15° and an aspect angle of 135° (southeast slope). Then, based on the image's imaging time, the solar altitude angle at that moment is calculated to be 60° and the solar azimuth angle to be 150°. Based on these angular parameters, the solar incidence angle correction factor for that pixel is calculated to be 0.85.

[0032] In another implementation of this invention, it is assumed that the uncorrected brightness value of pixel P in the multi-source remote sensing image is 120 in the visible light band (taking the red light band B4 as an example), 180 in the near-infrared band (B5), and 90 in the thermal infrared band (B10). First, terrain correction is performed, dividing the brightness value of each band by the corresponding pixel's solar incidence angle correction coefficient to obtain the corrected brightness values: visible light band 120 / 0.85≈141.2, near-infrared band 180 / 0.85≈211.8, and thermal infrared band 90 / 0.85≈105.9. Then, weighted fusion is performed according to a weight ratio of 3:1:2, calculating the final fused value of the pixel in the terrain-corrected fused image as: (141.2×3+211.8×1+105.9×2) / (3+1+2)≈141.2. All pixels perform this operation, ultimately generating a terrain-corrected fused image that covers the entire target grassland area.

[0033] In this embodiment of the invention, the Normalized Difference Vegetation Index (NDVI) and Land Surface Temperature (LST) are retrieved from topographically corrected fused images. Assuming a pixel A in the image has an NDVI value of 0.65 and an LST value of 28℃, this pixel is selected as a seed point. Using this seed point as the core, a growing region is initiated from its 8 neighboring areas, with similarity conditions set as follows: the absolute value of the difference between the NDVI value of the pixel to be merged and the average NDVI value of the current growing patch is less than 0.05, and the absolute value of the difference between the LST value and the average LST value of the current growing patch is less than 2℃. Neighboring pixel B has an NDVI of 0.62 and an LST of 29℃, meeting the conditions and is merged. This process is repeated until no new pixels meet the merging conditions, forming a closed homogeneous vegetation patch, denoted as patch M. After all pixels have been processed, a grassland homogeneous patch map is generated.

[0034] In one implementation of this invention, patch M is classified into community types. The average NDVI of all pixels within patch M is calculated to be 0.63, and it is determined to be of the "typical grassland" type based on its canopy structure parameters (such as texture features). According to a preset sensitivity determination matrix, the drought sensitivity coefficient corresponding to "typical grassland" is 3 (range 1-5). Assuming that the original drought value of patch M is calculated to be 0.5 (range 0-1) using the Temperature Vegetation Drought Index (TVDI) model, the final drought intensity value of this patch is: Drought intensity value = Original drought value × Drought sensitivity coefficient = 0.5 × 3 = 1.5. The same calculation is performed on all patches, and the local dry and wet boundaries are fitted according to the spatial distribution of drought intensity values. Finally, the entire area is divided into multiple units with different drought intensity values, generating a grassland partition drought grid.

[0035] In this embodiment of the invention, a standard carrying capacity is set for the target grassland area, and this value is determined based on historical data and grassland type. It is assumed that for a "typical grassland" area, the standard carrying capacity under no drought stress is 1.5 standard sheep units / hectare. Based on the grassland zoning drought grid, the grid cell containing patch M is extracted, and its drought intensity value is 1.5. Using a preset loss model, the drought intensity value is converted into a grass yield loss coefficient. Assuming a linear relationship between the loss coefficient and the drought intensity value: loss coefficient = drought intensity value / 5 (where 5 is the maximum value of the drought sensitivity coefficient), then the loss coefficient of this grid cell is 1.5 / 5 = 0.3.

[0036] In another implementation of this invention, the dynamic carrying capacity of the grid cell is calculated. The calculation formula is: Dynamic carrying capacity = Standard carrying capacity × (1 - Loss coefficient); Dynamic carrying capacity = 1.5 × (1 - 0.3) = 1.05 standard sheep units / hectare. This means that under the current drought stress, the carrying capacity of the grid cell has decreased from 1.5 to 1.05 sheep units / hectare.

[0037] In one implementation of this invention, real-time herd distribution data is received from an Internet of Things (IoT) device (such as a sheep wearing a GPS collar). Assume that the grid cell containing patch M has an area of ​​10 hectares. At a certain moment, the IoT data platform shows that there are 12 standard sheep units within this grid cell. The theoretical total carrying capacity of this grid cell is calculated: Total carrying capacity = Dynamic carrying capacity × Cell area = 1.05 × 10 = 10.5 standard sheep units. Then, the real-time livestock number is compared with the total carrying capacity: 12 standard sheep units > 10.5 standard sheep units. Since the real-time herd distribution data exceeds the dynamic carrying capacity, the grid cell is determined to be overgrazing.

[0038] In another implementation of this invention, dynamic management of livestock carrying capacity is performed. When an overgrazing condition is determined, a warning message is immediately sent to the mobile terminal of the ranch manager, stating: "Warning: Grid G-07 is overloaded! The current carrying capacity is 12 sheep units, while the approved carrying capacity is 10.5 sheep units. It is recommended to immediately move at least 2 sheep to the adjacent safe carrying capacity area F-07." If the approved carrying capacity of another grid unit H is 8 sheep units, but the real-time data feedback shows 7 sheep units, then the system determines that it is in a safe carrying capacity state and maintains the current grazing strategy without sending any instructions.

[0039] Preferably, step S1 includes the following steps:

[0040] Step S11: Simultaneously acquire multi-source remote sensing images covering the same target grassland area, including visible light, near-infrared, and thermal infrared bands;

[0041] Step S12: Calculate the slope value and aspect angle of each pixel in the digital elevation model;

[0042] Step S13: Determine the solar incidence angle correction coefficient for the multi-source remote sensing image based on the slope value and aspect angle;

[0043] Step S14: Divide the brightness value of each band in the multi-source remote sensing image by the solar incidence angle correction coefficient of the corresponding pixel, and then perform weighted fusion of the multi-source remote sensing images of visible light band, near-infrared band and thermal infrared band according to a weight ratio of 3:1:2 to obtain the terrain-corrected fused image.

[0044] In one implementation of this invention, the target grassland area is first determined to be a 100-square-kilometer experimental pastoral area. Remote sensing images and elevation data covering this area are acquired simultaneously. Specifically, the remote sensing images are obtained from Landsat 9 satellite OLI-2 / TIRS-2 sensor data, which includes visible light band (Band 4 red light band as an example), near-infrared band (Band 5), and thermal infrared band (Band 10) with a spatial resolution of 30 meters. The digital elevation model (DEM) uses SRTM (Shuttle Radar Topography Mission) data products, which also have a spatial resolution of 30 meters, to ensure a one-to-one correspondence with the pixels of the remote sensing images.

[0045] In one implementation of this invention, a surface analysis tool in a geographic information system (GIS) software is invoked to process the acquired DEM data. Specifically, taking any pixel P(i, j) in the DEM raster data as an example, its slope value and aspect angle are determined by calculating the rate of change of elevation between this pixel and its eight neighboring pixels. Assuming that after calculation, the slope value of pixel P(i, j) is determined to be 25°, and its aspect angle is determined to be 225° (i.e., southwest slope). This calculation process is applied to every pixel in the DEM data, thereby generating a slope layer and aspect layer with the exact same dimensions as the original DEM.

[0046] In one implementation of this invention, the calculation of the solar incidence angle correction coefficient depends not only on the terrain but also on the sun's position at the time of imaging. First, the solar zenith angle at the time of imaging is read from the Landsat 9 image's metadata file as 35° and the solar azimuth angle as 150°.

[0047] It should be noted that the solar incidence angle correction factor is essentially the cosine of the angle between the sunlight and the Earth's surface normal, cos(θ), and its mathematical expression is:

[0048] cos(θ)=cos(S)×cos(Z)+sin(S)×sin(Z)×cos(A_s-A_p);

[0049] Where S is the pixel slope, Z is the solar zenith angle, A_s is the solar azimuth angle, and A_p is the pixel aspect angle.

[0050] Specifically, taking the aforementioned pixel P(i,j) as an example, its slope S = 25° and aspect angle A_p = 225°. Substituting these parameters, the calculation is: cos(θ) = cos(25°) × cos(35°) + sin(25°) × sin(35°) × cos(150°―225°) ≈ 0.805. Therefore, the solar incidence angle correction factor for pixel P(i,j) is determined to be 0.805.

[0051] In one implementation of this invention, the original brightness values ​​(Digital Number, DN) of each band are first subjected to terrain correction. Assume that the original DN values ​​of pixel P(i,j) in the visible (red) band, near-infrared band, and thermal infrared band are 95, 170, and 105, respectively. The terrain-corrected DN values ​​are calculated as follows: Corrected visible light DN = 95 / 0.805 ≈ 117.9; Corrected near-infrared DN = 170 / 0.805 ≈ 211.2; Corrected thermal infrared DN = 105 / 0.805 ≈ 130.4.

[0052] In one implementation of this invention, the corrected DN values ​​of each band are weighted and fused. This fusion process aims to generate a single image layer that comprehensively reflects vegetation, moisture, and thermal conditions. The formula for calculating the fusion value is: Fusion value = (Corrected visible light DN × 3 + Corrected near-infrared DN × 1 + Corrected thermal infrared DN × 2) / (3 + 1 + 2); specifically, the final fusion value of pixel P(i, j) is: Fusion value = (117.9 × 3 + 211.2 × 1 + 130.4 × 2) / 6 ≈ 137.6. This topographic correction and weighted fusion process is applied to all pixels within the target grassland area, ultimately generating a complete topographically corrected and fused image that eliminates topographic influences and integrates multi-band information.

[0053] In one implementation of this invention, each pixel in the digital elevation model is classified into a slope grid and a slope aspect grid based on the slope value and aspect angle. The slope grid is reclassified, with a slope of 0°-5° defined as gentle, 5°-15° as sloping, and greater than 15° as steep. The aspect grid is also reclassified, with an aspect angle of 135°-225° defined as sunny slope and 315°-45° as shady slope.

[0054] Preferably, before identifying and segmenting spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature based on terrain-corrected fused images in step S2, the method further includes:

[0055] The normalized vegetation index is calculated based on the remote sensing images in the near-infrared and red bands of the topographically corrected fused image, and the original grassland surface temperature is obtained by inversion using the thermal infrared band of the topographically corrected fused image.

[0056] Determine the normalized vegetation index (NDI) of each pixel in the digital elevation model. If it is greater than 0.75, the pixel is marked as a pure vegetation pixel; if it is less than 0.2, it is marked as a pure bare soil pixel; if it is between 0.2 and 0.75, it is marked as a mixed pixel.

[0057] Centered on the mixed pixel, search for all pixels marked as pure bare soil pixels within a preset semi-circular neighborhood, and calculate the arithmetic mean of their surface temperatures based on the original grassland surface temperature. Use this average value as the bare soil background temperature of the area.

[0058] The original grassland surface temperature corresponding to the mixed pixel is corrected based on the bare soil background temperature to obtain the corrected grassland surface temperature.

[0059] The digital elevation models labeled with pure vegetation pixels, pure bare soil pixels, and mixed pixels were mapped to the corresponding modified grassland surface temperature and normalized vegetation index to obtain grassland vegetation distribution data.

[0060] In one implementation of this invention, reflectance layers for the near-infrared (NIR) and red bands are separated from the topographically corrected fused image. Specifically, taking pixel A in the image as an example, its corrected red band reflectance is 0.08, and its near-infrared band reflectance is 0.52. The Normalized Difference Vegetation Index (NDVI) is calculated using the following mathematical expression: NDVI = (NIR - Red) / (NIR + Red). Substituting the data for pixel A: NDVI = (0.52 - 0.08) / (0.52 + 0.08) = 0.44 / 0.60 ≈ 0.733. Simultaneously, the surface temperature is retrieved using the thermal infrared band layer from the topographically corrected fused image through a single-window algorithm. This algorithm requires input parameters such as surface emissivity and atmospheric transmittance, which can be preset based on the NDVI value and an atmospheric parameter database. After the inversion calculation, the original grassland surface temperature (LST) of pixel A is found to be 31.5℃.

[0061] In one implementation of this invention, each pixel is classified and labeled based on the NDVI layer generated in the previous step. Specifically, taking pixel A as an example, its NDVI value is 0.733. Since this value is determined to be between a preset range of 0.2 and 0.75, pixel A is labeled as a mixed pixel.

[0062] In another implementation of this invention, it is assumed that pixel B has an NDVI value of 0.81, which is greater than 0.75, and is therefore labeled as a pure vegetation pixel. It is also assumed that pixel C has an NDVI value of 0.15, which is less than 0.2, and is therefore labeled as a pure bare soil pixel. This classification and labeling process will traverse the entire study area, generating a pixel type classification map.

[0063] In one implementation of this invention, a neighborhood search is performed for each pixel marked as a hybrid pixel.

[0064] It should be noted that the "semi-circular neighborhood" here is preset for more accurately obtaining representative bare soil background temperatures. Its range can be set to a radius of 20 pixels (i.e., 600 meters), oriented in the opposite direction of the solar azimuth, to prioritize searching for bare soil in shady or drier areas. Specifically, the search is performed within the preset semi-circular neighborhood centered on the mixed pixel A. Assuming five pixels are found labeled as pure bare soil pixels, their corresponding original grassland surface temperatures are 42.1℃, 43.5℃, 41.8℃, 42.9℃, and 43.2℃, respectively. The arithmetic mean of these five temperature values ​​yields the bare soil background temperature for that area: Bare soil background temperature = (42.1 + 43.5 + 41.8 + 42.9 + 43.2) / 5 = 213.5 / 5 = 42.7℃.

[0065] In one implementation of this invention, the correction process aims to eliminate the effect of bare soil in a mixed pixel on the overall pixel temperature. Specifically, taking mixed pixel A as an example, its original grassland surface temperature is 31.5℃, and its calculated bare soil background temperature is 42.7℃. The correction process is adjusted based on the pixel's vegetation cover (which can be estimated by NDVI). Assuming that the vegetation cover of pixel A is estimated to be 78% based on its NDVI value (0.733), a correction model is applied that considers the pixel temperature to be a weighted average of the vegetation and bare soil temperatures. By introducing a localized bare soil background temperature, the true temperature of the vegetation portion can be estimated more accurately. After correction calculation, the corrected grassland surface temperature of pixel A is determined to be 29.8℃.

[0066] It should be noted that for pixels marked as pure vegetation pixels and pure bare soil pixels, their surface temperature does not need to be corrected; the corrected grassland surface temperature is the same as their original grassland surface temperature.

[0067] In one implementation of this invention, a new multi-layer raster dataset, namely grassland vegetation distribution data, is created. Each layer of this dataset stores one type of attribute information. Specifically, the operations are as follows: First layer: stores the markers (pure vegetation, pure bare soil, mixed pixels) in the land cover type classification map. Second layer: stores the corrected grassland surface temperature values ​​for all pixels (including pure vegetation, pure bare soil, and corrected mixed pixels). Third layer: stores the Normalized Difference Vegetation Index (NDVI) values ​​for all pixels. By precisely aligning these three layers of data spatially, each pixel location simultaneously possesses three key attributes: type, temperature, and vegetation index.

[0068] Preferably, step S2, which involves identifying and segmenting spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature based on terrain-corrected fused images, includes:

[0069] Based on grassland vegetation distribution data, local extreme high and low value points of vegetation index are identified, and the corresponding pixels are defined as growth seed points for lush vegetation patches and sparse vegetation patches, respectively.

[0070] Taking the seed point as the core, the region that initiates growth is used as a growth patch, and the vegetation index and surface temperature value of the neighboring pixels are obtained.

[0071] Based on the vegetation index and surface temperature values ​​of neighboring pixels, neighboring pixels that meet the similarity conditions of vegetation index and surface temperature are identified as pixels to be merged and merged into the current growing patch until no new pixels meet the merging conditions, forming a closed boundary of homogeneous vegetation patches, and finally generating a grassland homogeneous patch map.

[0072] In one implementation of this invention, a Normalized Difference Vegetation Index (NDVI) layer is extracted from grassland vegetation distribution data. To identify local extrema, a 3×3 moving window is used to traverse the entire NDVI layer.

[0073] Specifically, when the center of the moving window is located at cell X, the NDVI value of the center cell X is compared with the NDVI values ​​of its eight neighboring cells. If the NDVI value of cell X is greater than the values ​​of all eight neighboring cells, then cell X is identified as a local extreme point and defined as a seed point for the growth of lush vegetation patches.

[0074] In another implementation of this invention, the same method is used to identify local minimum values ​​(LNUs). If the NDVI value of pixel Y is less than the values ​​of all eight of its neighboring pixels, then pixel Y is identified as a LNU and defined as a seed point for the growth of sparse vegetation patches. For example, assuming pixel X has an NDVI of 0.82 and the NDVI values ​​of its eight surrounding pixels are all between 0.75 and 0.80, then pixel X is marked as a seed point for lush vegetation. Assuming pixel Y has an NDVI of 0.21 and the NDVI values ​​of its eight surrounding pixels are all between 0.23 and 0.28, then pixel Y is marked as a seed point for sparse vegetation.

[0075] In one implementation of this invention, an unprocessed seed point (e.g., a lush vegetation seed point X) is selected from the seed point set generated in the previous step as the starting core to initiate the regional growth process. The initial growth patch contains only the seed point X itself. Specifically, the eight spatial neighboring pixels of seed point X (i.e., the eight pixels directly adjacent to X) are identified as the neighborhood to be checked. For each neighboring pixel to be checked, its corresponding NDVI value and corrected grassland surface temperature (LST) value are retrieved from the grassland vegetation distribution data. For example, if a neighboring pixel of seed point X is P1, the NDVI value of P1 is 0.80 and the LST value is 26.1℃, as read from the grassland vegetation distribution data.

[0076] In one implementation of this invention, the determination of similarity conditions is the core of region growing. These conditions are preset as follows: the absolute value of the difference between the NDVI value of the pixel to be inspected and the average NDVI value of all pixels within the current growing patch is less than a first threshold T_NDVI (e.g., T_NDVI = 0.05); and the absolute value of the difference between the LST value of the pixel to be inspected and the average LST value of all pixels within the current growing patch is less than a second threshold T_LST (e.g., T_LST = 1.5℃). Only when a neighboring pixel simultaneously meets both of these conditions is it determined to be a pixel to be merged.

[0077] Specifically, taking the seed point X and neighboring pixel P1 as an example, the initial growth patch contains only X, with an average NDVI of 0.82 and an average LST of 25.8℃. The neighboring pixel P1 is evaluated as follows: NDVI difference = |0.80 - 0.82| = 0.02. Since 0.02 < 0.05 (T_NDVI), condition one is satisfied. LST difference = |26.1 - 25.8| = 0.3℃. Since 0.3 < 1.5 (T_LST), condition two is satisfied. Because P1 satisfies both conditions, it is identified as a pixel to be merged and is merged into the current growth patch.

[0078] It is important to note that the average NDVI and average LST values ​​of the growth patch are recalculated whenever a new pixel is merged. For example, after P1 is merged, the new growth patch contains X and P1, and its new average NDVI is (0.82+0.80) / 2 = 0.81, and its new average LST is (25.8+26.1) / 2 = 25.95℃.

[0079] In another implementation of this invention, the merging process continues. Unchecked pixels in the neighborhood of the newly merged pixel P1 are added to a queue to be checked. Then, the next pixel to be checked is taken from the queue and compared with the updated patch average. This iterative process repeats until the queue is empty, meaning all neighboring pixels have been checked and no new pixels meet the merging conditions. At this point, a closed, homogeneous patch boundary is formed, with internal vegetation growth highly similar to surface temperature.

[0080] Preferably, determining the similarity conditions satisfied by the pixels to be merged includes:

[0081] The absolute value of the difference between the vegetation index value of the pixel to be merged and the average vegetation index value of the current growing patch is used to obtain the pixel to be merged.

[0082] The absolute value of the difference between the surface temperature value of the pixel to be merged and the average surface temperature value of the current growing patch is calculated to obtain the surface temperature difference value.

[0083] If the difference in vegetation index is less than a preset first threshold and the difference in surface temperature is less than a preset second threshold, then the pixel is determined to meet the merging conditions.

[0084] In one implementation of this invention, this determination process is the core of the region growing algorithm, used to determine whether a neighboring cell should be incorporated into a growing homogeneous patch.

[0085] It should be noted that the first threshold (T_NDVI) and the second threshold (T_LST) are preset according to the ecological characteristics of the target grassland area and the quality of remote sensing data. For example, for a typical grassland with a relatively single vegetation type, T_NDVI can be set to 0.05 and T_LST can be set to 1.5°C. These thresholds directly control the internal homogeneity degree of the finally generated homogeneous patches.

[0086] In an implementation manner of the embodiment of the present invention, assume that a regional growth process is in progress, and the current growing patch consists of 10 pixels. The average vegetation index (NDVI_avg) of this patch has been calculated to be 0.68, and the average land surface temperature (LST_avg) is 29.5°C. Now, it is necessary to determine the similarity conditions for a neighboring pixel P_neighbor of this patch. The vegetation index value (NDVI_neighbor) of P_neighbor is read from the grassland vegetation distribution data as 0.71, and the land surface temperature value (LST_neighbor) is 28.8°C.

[0087] In another implementation manner of the embodiment of the present invention, consider the situation of another neighboring pixel Q_neighbor. Assume that its vegetation index value is 0.75 and the land surface temperature value is 30.1°C. Calculate the absolute value of the vegetation index difference: ΔNDVI = |0.75 - 0.68| = 0.07. Calculate the absolute value of the land surface temperature difference: ΔLST = |30.1°C - 29.5°C| = 0.6°C. Conduct a logical judgment: Judgment 1 (vegetation index): ΔNDVI (0.07) < T_NDVI (0.05)? The result is false. Judgment 2 (land surface temperature): ΔLST (0.6°C) < T_LST (1.5°C)? The result is true.

[0088] It should be noted that the merging condition requires that both judgments must be true. Since the vegetation index difference of pixel Q_neighbor exceeds the preset first threshold, even though its land surface temperature difference is within the allowable range, it is finally determined that this pixel does not meet the merging condition.

[0089] Preferably, before assigning a drought sensitivity coefficient to each patch in the grassland homogeneous patch map in step S2, determining the drought intensity value of each patch according to this coefficient, and then fitting the local wet-dry edge and generating a grassland partition drought grid, it further includes:

[0090] According to each patch in the grassland homogeneous patch map, extract the vegetation index values of all pixels corresponding to it in the grassland vegetation distribution data;

[0091] Calculate the mean value of the vegetation index values of all pixels within each patch as a parameter representing the vegetation coverage;

[0092] Calculate the ratio of the near-infrared band to the visible light band in the region corresponding to the pure vegetation pixel within each patch, and determine the parameters characterizing the complexity of the canopy structure based on the ratio;

[0093] Based on parameters characterizing vegetation cover and parameters characterizing canopy structure complexity, patches are classified into grassland community distribution types, and each patch is assigned a drought sensitivity coefficient ranging from 1 to 5 through a preset sensitivity determination matrix.

[0094] Each patch is spatially correlated with its assigned drought sensitivity coefficient value to generate a patch stress sensitivity map.

[0095] In one implementation of this invention, taking a patch (denoted as patch A) in a grassland homogeneous patch map as an example, this patch consists of 50 pixels. Using the spatial boundary of patch A as a mask, the Normalized Difference Vegetation Index (NDVI) values ​​corresponding to these 50 pixels are accurately extracted from the grassland vegetation distribution data. Specifically, a spatial query operation is performed, returning a list containing 50 NDVI values, for example: [0.72, 0.68, 0.75, ..., 0.71].

[0096] In one implementation of this invention, the arithmetic mean of the 50 NDVI values ​​extracted in the previous step is calculated. Specifically, the calculation process is as follows: Average NDVI = (0.72 + 0.68 + 0.75 + ... + 0.71) / 50; assuming the calculation result is 0.73. Therefore, the parameter representing vegetation cover for patch A is determined to be 0.73.

[0097] In one implementation of this invention, firstly, from 50 pixels in patch A, all pixels labeled as "pure vegetation pixels" are selected based on pixel type markers in the grassland vegetation distribution data. Assume 40 pure vegetation pixels are selected. Next, from the terrain-corrected fused image dataset, the reflectance values ​​in the near-infrared (NIR) and visible red (Red) bands corresponding to these 40 pure vegetation pixels are extracted. Specifically, for each of these 40 pixels, its ratio (Ratio = NIR / Red) is calculated. For example, if one pixel has an NIR of 0.55 and a Red of 0.10, its ratio is 5.5. The ratios of all 40 pure vegetation pixels are calculated, and their average value is taken. Assume the final average ratio is 4.1. Then, according to preset rules, this average ratio is converted into canopy structure parameters: if the average ratio > 3.5, the parameter is "high canopy" (corresponding to a canopy height of 20-40 cm); if the average ratio is between 2.0 and 3.5, the parameter is "medium canopy" (corresponding to a canopy height of 10-20 cm); if the average ratio < 2.0, the parameter is "low canopy" (corresponding to a canopy height of less than 10 cm). Since the calculated average ratio is 4.1, which is greater than 3.5, the parameter representing the canopy structure complexity of patch A is determined to be "high canopy".

[0098] In one implementation of this invention, a preset sensitivity determination matrix constructed based on field survey data is used. For example:

[0099] Vegetation coverage parameters Canopy structural parameters Grassland community distribution types drought sensitivity coefficient >0.7 High Crown meadow grassland 5 0.4-0.7 Middle canopy Typical grassland 3 <0.4 Low canopy desert steppe 1

[0100] Specifically, the two parameters of patch A were substituted into the matrix for matching: the parameter representing vegetation cover (0.73) and the parameter representing canopy structure complexity ("high canopy"). The matrix query revealed that 0.73 > 0.7, and the structure was "high canopy," perfectly matching the first row of rules. Therefore, patch A was classified as "meadow steppe" and assigned a drought sensitivity coefficient of 5.

[0101] It should be noted that this sensitivity determination matrix can be customized and expanded according to grassland types in different geographical regions to improve the accuracy of classification and sensitivity assessment.

[0102] In one implementation of this invention, the calculated non-spatial attribute (drought sensitivity coefficient) is assigned to the spatial entity (patch). Specifically, a new raster layer is created, whose spatial extent and resolution are completely identical to the grassland homogeneous patch map. Then, all 50 pixels belonging to patch A in the grassland homogeneous patch map are assigned a value of 5 at their corresponding positions in the new layer.

[0103] In another implementation of this invention, suppose there is another patch B, which, after the above calculation, is classified as "typical grassland" with a drought sensitivity coefficient of 3. Then, in the new layer, all pixel positions belonging to patch B will be assigned a value of 3. After this process traverses all patches, a complete patch stress sensitivity map is finally generated.

[0104] Preferably, the determination of parameters characterizing the complexity of the canopy structure further includes:

[0105] Calculate vegetation cover for mixed pixels;

[0106] The standard deviation of the normalized vegetation index within the mixed pixels was calculated using a 3×3 pixel window.

[0107] When the standard deviation is greater than the preset difference threshold, for areas with vegetation coverage greater than 60%, the first derivative of the reflectance in the near-infrared band of the topographically corrected fused image is extracted.

[0108] When the first derivative value is in the range of 0.02-0.05, it is identified as high canopy density grassland; when the first derivative value is less than 0.02, it is identified as low and sparse grassland.

[0109] In one implementation of this invention, a more refined canopy structure analysis is performed on pixels labeled as "mixed pixels." First, the fractional vegetation coverage (FVC) of these mixed pixels needs to be calculated. Specifically, the FVC is calculated using a pixel-based binary model, and its mathematical expression is:

[0110] FVC=(NDVI―NDVI_soil) / (NDVI_veg―NDVI_soil);

[0111] Wherein, NDVI is the normalized vegetation index value of the mixed pixel; NDVI_soil is the NDVI value of the pure bare soil pixel; and NDVI_veg is the NDVI value of the pure vegetation pixel.

[0112] It should be noted that the values ​​of NDVI_soil and NDVI_veg can be obtained from the NDVI frequency histogram of the entire study area. Typically, the NDVI value corresponding to 5% of the cumulative frequency is taken as NDVI_soil, and the NDVI value corresponding to 95% is taken as NDVI_veg. Assume that NDVI_soil = 0.08 and NDVI_veg = 0.85 are determined using this method. Taking a pixel P labeled as a "mixed pixel" as an example, its NDVI value is 0.58. Substituting into the formula to calculate its vegetation cover: FVC = (0.58 - 0.08) / (0.85 - 0.08) = 0.50 / 0.77 ≈ 0.649, or 64.9%.

[0113] In one implementation of this invention, a 3×3 pixel window is constructed centered on the mixed pixel P. Then, the NDVI values ​​of nine pixels within this window are extracted from the NDVI layer. Specifically, assume the nine extracted NDVI values ​​are: [0.55, 0.59, 0.60, 0.56, 0.58, 0.61, 0.54, 0.57, 0.62]. The standard deviation of these nine values ​​is calculated. The standard deviation is first calculated by taking the mean: (0.55 + ... + 0.62) / 9 ≈ 0.58. Then, the sum of the squares of the differences between each value and the mean is calculated, divided by the sample size minus one, and finally the square root is taken. After calculation, the standard deviation of the NDVI values ​​within this window is 0.028.

[0114] In one implementation of this invention, the standard deviation calculated in the previous step is compared with a preset difference threshold (e.g., T_std = 0.025). Simultaneously, the vegetation cover of the pixel is compared with a preset coverage threshold (60%). Specifically, for pixel P: Judgment 1 (Standard Deviation): Standard deviation (0.028) > difference threshold (0.025)? The result is true. Judgment 2 (Vegetation Coverage): Vegetation cover (64.9%) > coverage threshold (60%)? The result is true. Since both conditions are true, the next step is performed on the pixel, namely, extracting the first derivative value of the near-infrared band reflectance.

[0115] It should be noted that a large standard deviation usually means that the surface cover in that local area is highly heterogeneous, with vegetation patches and bare soil interspersed.

[0116] In another implementation of this invention, if the standard deviation of another mixed cell Q is 0.018 (less than 0.025), then condition one is not met, and the first derivative calculation will not be performed on the cell. Its canopy structure parameters will be determined by other methods (such as the aforementioned ratio method) or default values.

[0117] In one implementation of this invention, the calculation of the first derivative is used to reflect the rate of change of near-infrared reflectance in space, which is closely related to the density and structure of the canopy. The first derivative is estimated by analyzing the near-infrared reflectance values ​​of pixel P and its neighboring pixels. Specifically, assuming the near-infrared reflectance of pixel P is 0.45, and the near-infrared reflectances of its neighboring pixels (e.g., along the east-west direction) are 0.42 and 0.49, respectively. The first derivative value can be approximated as the change in reflectance divided by the pixel distance (30 meters). After calculation, the first derivative value at pixel P is approximately 0.04. Then, this first derivative value is compared with a preset interval: Is the first derivative value (0.04) within the range of 0.02-0.05? The result is true. Therefore, the canopy structure parameter of pixel P is identified as "high canopy density grassland". This parameter will serve as supplementary or alternative information characterizing the complexity of the canopy structure for subsequent grassland community type classification and drought sensitivity coefficient assignment.

[0118] In another implementation of this invention, if the first derivative of another pixel R is calculated to be 0.015, its canopy structure parameter will be identified as "low and sparse grassland" because this value is less than 0.02.

[0119] Preferably, step S2, which involves assigning a drought sensitivity coefficient based on the community type in the grassland homogeneous patch map, determining the drought intensity value of each patch based on this coefficient, and then fitting the local wet-dry edges to generate a grassland zoning drought grid, includes:

[0120] For each patch in the patch stress sensitivity map, the mean surface temperature and mean vegetation index are extracted from the grassland vegetation distribution data.

[0121] The original drought degree of the current water stress state in the patch is calculated based on the average surface temperature and the average vegetation index.

[0122] The original drought intensity is multiplied by the drought sensitivity coefficient corresponding to the patch in the patch stress sensitivity map to obtain the final drought intensity value of the patch.

[0123] The local wet-dry edges of the target grassland region are fitted based on the drought intensity values, and the regions of the local wet-dry edges are closed and connected to generate a grassland zoning drought grid.

[0124] In one implementation of this invention, a patch (denoted as patch A) in a patch stress sensitivity map is used as an example. First, using the spatial boundary of patch A as a mask, the corrected grassland surface temperature (LST) and normalized difference vegetation index (NDVI) values ​​corresponding to all pixels within the patch are extracted from the grassland vegetation distribution data. Specifically, assuming patch A consists of 50 pixels, two lists containing 50 values ​​are extracted. Then, the arithmetic mean of these two lists is calculated. For example, the calculated average surface temperature (LST_avg) of patch A is 30.2℃; the average vegetation index (NDVI_avg) of patch A is 0.73.

[0125] In one implementation of this invention, the original drought degree is calculated using the Temperature Vegetation Drought Index (TVDI) model. This model constructs a feature space consisting of surface temperature and vegetation index.

[0126] It should be noted that constructing a TVDI model requires first determining the "dry edge" and "wet edge" of the study area. The "dry edge" represents the surface temperature under specific vegetation cover when water stress is most severe, while the "wet edge" represents the surface temperature when water is plentiful. These two edges are typically obtained by linear regression fitting of the LST-NDVI scatter plot of the entire study area. Assuming that the fitting yields the equation for the dry edge as LST_dry = -15 × NDVI + 48 and the equation for the wet edge as LST_wet = -10 × NDVI + 32, the mathematical expression for TVDI is: TVDI = (LST_avg - LST_wet) / (LST_dry - LST_wet); specifically, for patch A, its NDVI_avg = 0.73. First, calculate the dry and wet edge temperatures corresponding to this NDVI: LST_wet_A = 24.7℃; LST_dry_A = 37.05℃; then substitute LST_avg (30.2℃) of patch A into the TVDI formula. Therefore, the original drought degree of patch A representing the current water stress state is determined to be 0.445.

[0127] In one implementation of this invention, the drought sensitivity coefficient corresponding to patch A is queried from the patch stress sensitivity map. Assume that in the previous steps, patch A was classified as "meadow steppe" and assigned a drought sensitivity coefficient of 5. Then, a multiplication operation is performed: final drought intensity value = original drought intensity × drought sensitivity coefficient; final drought intensity value = 0.445 × 5 = 2.225; therefore, the final drought intensity value of patch A is determined to be 2.225.

[0128] In one implementation of this invention, the final drought intensity value of all patches is assigned to their corresponding pixels to generate a drought intensity distribution map. Then, the map is segmented according to a preset drought level threshold. Specifically, the drought level thresholds are assumed to be set as follows: mild drought: drought intensity value 1.0-2.0; moderate drought: drought intensity value 2.0-3.0; severe drought: drought intensity value > 3.0. The drought intensity distribution map is reclassified, merging adjacent pixels whose values ​​fall within the same interval to form regional patches representing different drought levels. The boundaries of these patches are the local wet / dry edges. These closed drought level regional patches are then vectorized or gridded to generate the final grassland zoning drought grid.

[0129] Of particular importance, the original drought degree for calculating the current water stress state in patches based on mean surface temperature and mean vegetation index can also include:

[0130] Using grassland homogeneous patch images as masks, remote sensing images of land surface temperature and vegetation index at the same time phase are partitioned.

[0131] For each microenvironment zone represented by the identification code, the average surface temperature and average vegetation index of all pixels within it are calculated independently, and these two values ​​are combined to generate the temperature-vegetation baseline pair for that zone.

[0132] By traversing each pixel in the study area, the true surface temperature value and vegetation index value are obtained, and the corresponding value in the temperature-vegetation reference pair corresponding to the microenvironment zone to which the pixel belongs is subtracted to obtain the temperature anomaly value and vegetation anomaly value.

[0133] Create dual-band raster data, write the temperature anomaly values ​​of all pixels into the first band, and write the vegetation anomaly values ​​into the second band to generate a relative stress field.

[0134] The original drought degree is calculated based on the relative stress field.

[0135] In one implementation of this invention, each patch in the grassland homogeneous patch image has a unique identifier, representing an independent microenvironment partition. Specifically, taking patch A (identifier code 1, representing the "meadow steppe" microenvironment) as an example, its spatial boundary is used as a mask to filter out the LST and NDVI values ​​of all pixels belonging to that patch from the grassland vegetation distribution data. Similarly, the same operation is performed on patch B (identifier code 2, representing the "typical grassland" microenvironment), thereby logically segmenting the continuous LST and NDVI images into array sets that correspond one-to-one with the microenvironment partitions.

[0136] In one implementation of this invention, statistical calculations are performed independently on each array set obtained from the previous partitioning step. Specifically, for patch A with identifier code 1, assuming it contains 50 pixels, its average LST is calculated to be 28.5℃ and its average NDVI is 0.75. These two averages are combined to form the temperature-planting baseline pair for this patch, denoted as B_A = (28.5, 0.75). For patch B with identifier code 2, assuming it contains 120 pixels, its average LST is calculated to be 31.0℃ and its average NDVI is 0.60. These two averages are combined to form the temperature-planting baseline pair for this patch, denoted as B_B = (31.0, 0.60).

[0137] In one implementation of this invention, each pixel within the study area is traversed, and anomaly calculations are performed. Specifically, taking pixel P as an example, the microenvironmental zone to which it belongs is first determined by querying the grassland homogeneous patch map. Assume pixel P belongs to patch A, identified by code 1. Then, the true surface temperature of pixel P is read from the grassland vegetation distribution data as 30.0℃, and the vegetation index is 0.71. Next, the corresponding values ​​in the temperature-vegetation baseline pair B_A(28.5, 0.75) corresponding to the zone to which pixel P belongs are subtracted: Temperature anomaly = True LST - Baseline LST = 30.0℃ - 28.5℃ = +1.5℃. Vegetation anomaly = True NDVI - Baseline NDVI = 0.71 - 0.75 = -0.04. A positive temperature anomaly indicates that the point is hotter than the average state of its microenvironment, while a negative vegetation anomaly indicates that its vegetation growth is worse than the average state; both indicate the presence of water stress.

[0138] In one implementation of this invention, a new dual-band raster data file with the exact same size and resolution as the original remote sensing image is created.

[0139] In one implementation of this invention, the original drought degree is calculated based on two anomalies in the relative stress field. To eliminate the influence of dimensions, the anomalies of the two bands first need to be normalized, mapping them to the range of 0-1. Specifically, the maximum-minimum normalization method is used. Assume the temperature anomaly range for the entire study area is [-5℃, +10℃], and the vegetation anomaly range is [-0.2, +0.1]. After normalization, the temperature anomaly of pixel P (+1.5) becomes 0.43, and the vegetation anomaly (-0.04) becomes 0.53 (it should be noted that the smaller the vegetation anomaly, the more severe the stress; therefore, the normalization process needs to be reversed). Then, the original drought degree is calculated by weighted summation, and its mathematical expression is:

[0140] Original drought degree = W_T × (normalized temperature anomaly) + W_V × (normalized vegetation anomaly);

[0141] Where W_T and W_V are weighting coefficients, and W_T + W_V = 1. Since temperature responds more rapidly to drought, W_T = 0.6 and W_V = 0.4 can be set. Substituting the normalized anomaly value of pixel P: raw drought degree = 0.47. Therefore, the raw drought degree of pixel P is finally determined to be 0.47.

[0142] Preferably, determining the solar incidence angle correction coefficient for multi-source remote sensing images based on slope value and aspect angle includes:

[0143] Extract the solar azimuth angle of each pixel in the digital elevation model corresponding to the image imaging time from multi-source remote sensing images;

[0144] Based on the range of slope values, the solar incidence angle correction factor is dynamically calculated using the slope aspect angle and solar azimuth angle.

[0145] In one implementation of this invention, the solar azimuth describes the sun's horizontal position in the sky, and this information is typically recorded in the metadata file of a remote sensing image. Specifically, the metadata file (usually in MTL.txt format) of the multi-source remote sensing image used (e.g., Landsat 9 image) is first located and parsed. In this file, fields related to the sun's geometric position are searched. For example, the line "SUN_AZIMUTH=152.50" is read. This means that at the time the image was captured, the sun's azimuth was 152.5 degrees (measured clockwise from true north).

[0146] It should be noted that this solar azimuth angle value is for the entire image scene and can be considered a constant within a small study area. Therefore, this value will be applied to the subsequent calculations for each pixel within that area.

[0147] Of particular importance, the solar incidence angle correction factor can be calculated as follows:

[0148] When the slope value is within the range of 15°-30°, the solar incidence angle correction factor is set to 1.2×cos(slope aspect angle - solar azimuth angle);

[0149] When the slope is greater than 30°, the solar incidence angle correction factor is set to 1.5×cos(slope angle - solar azimuth angle).

[0150] In another implementation of this invention, it is assumed that the solar azimuth angle obtained from the metadata is 145°. Each cell is traversed, and a conditional judgment is performed based on its slope value. Different formulas are applied to calculate its solar incidence angle correction coefficient. Calculations are performed for cells with slope values ​​within the range of 15°-30°. Assume the system processes cell A, whose value in the slope layer is 20° and its value in the aspect layer is 180° (due south). First, it is determined that its slope value of 20° falls within the preset range of 15°-30°, thus triggering the first calculation rule.

[0151] Specifically, the solar incidence angle correction factor is calculated as follows: correction factor = 1.2 × cos(slope angle - solar azimuth angle); according to the trigonometric function value, cos(35°) ≈ 0.819, then: correction factor = 1.2 × 0.819 ≈ 0.983; therefore, the solar incidence angle correction factor for pixel A is determined to be 0.983, and this value will be used to correct the brightness value of pixel A in the future.

[0152] In another implementation of this invention, calculations are performed for pixels with slope values ​​greater than 30°. Assume the system processes pixel B, whose value in the slope layer is 35° and its value in the aspect layer is 270° (due west slope). The slope value of 35° is determined to be greater than the preset threshold of 30°, thus triggering the second calculation rule. Specifically, the solar incidence angle correction coefficient is calculated as follows: Correction coefficient = 1.5 × cos(aspect angle - solar azimuth angle) ≈ -0.861;

[0153] It should be noted that when the calculated solar incidence angle correction factor is negative or zero, it usually indicates that the pixel is in the self-shadow zone of the terrain and cannot receive direct sunlight.

[0154] See Figure 2 The present invention also provides a grassland drought identification and assessment system based on multi-source remote sensing image data, which executes the grassland drought identification and assessment method based on multi-source remote sensing image data as described above. The grassland drought identification and assessment system based on multi-source remote sensing image data includes:

[0155] S101: Terrain Correction and Fusion Module, used to acquire multi-source remote sensing images and digital elevation models covering the target grassland area, and perform terrain correction and band-weighted fusion on the multi-source remote sensing images based on the slope and aspect in the digital elevation model to obtain a terrain-corrected fused image.

[0156] S102: Arid Zoning Assessment Module, used to identify and segment spatially continuous grassland homogeneous patches based on topographically corrected fused images and where vegetation growth and surface temperature are similar; combine the community type in the grassland homogeneous patch map to assign a drought sensitivity coefficient, determine the drought intensity value of each patch based on the coefficient, and then fit the local dry and wet edges to generate a grassland zoning drought grid.

[0157] S103: Carrying capacity calculation module, used to determine the dynamic carrying capacity of each grid cell under the current drought stress based on the drought intensity value in the grassland zoning drought grid;

[0158] S104: Overgrazing monitoring and management module, used to determine the real-time carrying capacity by comparing the dynamic carrying capacity limit with the real-time distribution data of the herd transmitted by the Internet of Things device. If the real-time distribution data of the herd exceeds the dynamic carrying capacity limit, it is determined to be an overgrazing state and dynamic carrying capacity management is implemented; otherwise, it is determined to be a safe carrying capacity state and the current grazing strategy is maintained.

[0159] Therefore, the embodiments should be considered exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0160] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for identifying and assessing grassland drought based on multi-source remote sensing image data, characterized in that, Includes the following steps: Step S1: Acquire multi-source remote sensing images and digital elevation models covering the target grassland area, and perform terrain correction and band weighted fusion on the multi-source remote sensing images based on the slope and aspect in the digital elevation model to obtain a terrain-corrected fused image. Step S2: Based on topographic correction and image fusion, identify and segment spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature; assign drought sensitivity coefficients based on community types in the grassland homogeneous patches, evaluate the drought intensity of each patch based on these coefficients, and then fit local wet-dry edges to generate a grassland zoning drought grid; before identifying and segmenting spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature based on topographic correction and image fusion in step S2, the following steps are also included: The normalized vegetation index is calculated based on the remote sensing images in the near-infrared and red bands of the topographically corrected fused image, and the original grassland surface temperature is obtained by inversion using the thermal infrared band of the topographically corrected fused image. Determine the normalized vegetation index (NDI) of each pixel in the digital elevation model. If it is greater than 0.75, the pixel is marked as a pure vegetation pixel; if it is less than 0.2, it is marked as a pure bare soil pixel; if it is between 0.2 and 0.75, it is marked as a mixed pixel. Centered on the mixed pixel, search for all pixels marked as pure bare soil pixels within a preset semi-circular neighborhood, and calculate the arithmetic mean of their surface temperatures based on the original grassland surface temperature. Use this average value as the bare soil background temperature of the area. The original grassland surface temperature corresponding to the mixed pixel is corrected based on the bare soil background temperature to obtain the corrected grassland surface temperature. The digital elevation models of pure vegetation pixels, pure bare soil pixels, and mixed pixel markers are mapped to the corresponding modified grassland surface temperature and normalized vegetation index to obtain grassland vegetation distribution data. Step S3: Determine the dynamic carrying capacity of each grid cell under the current drought stress based on the drought intensity value in the grassland drought grid; Step S4: The dynamic carrying capacity limit is compared with the real-time herd distribution data transmitted by the IoT device to determine the real-time carrying capacity. If the real-time herd distribution data exceeds the dynamic carrying capacity limit, it is determined to be an overgrazing state and dynamic carrying capacity management is implemented; otherwise, it is determined to be a safe carrying capacity state and the current grazing strategy is maintained.

2. The grassland drought identification and assessment method based on multi-source remote sensing image data according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: Simultaneously acquire multi-source remote sensing images covering the same target grassland area, including visible light, near-infrared, and thermal infrared bands; Step S12: Calculate the slope value and aspect angle of each pixel in the digital elevation model; Step S13: Determine the solar incidence angle correction coefficient for the multi-source remote sensing image based on the slope value and aspect angle; Step S14: Divide the brightness value of each band in the multi-source remote sensing image by the solar incidence angle correction coefficient of the corresponding pixel, and then perform weighted fusion of the multi-source remote sensing images of visible light band, near-infrared band and thermal infrared band according to a weight ratio of 3:1:2 to obtain the terrain-corrected fused image.

3. The grassland drought identification and assessment method based on multi-source remote sensing image data according to claim 1, characterized in that, Step S2 involves identifying and segmenting spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature based on topographically corrected fused images, including: Based on grassland vegetation distribution data, local extreme high and low value points of vegetation index are identified, and the corresponding pixels are defined as growth seed points for lush vegetation patches and sparse vegetation patches, respectively. Taking the seed point as the core, the region that initiates growth is used as a growth patch, and the vegetation index and surface temperature value of the neighboring pixels are obtained. Based on the vegetation index and surface temperature values ​​of neighboring pixels, neighboring pixels that meet the similarity conditions of vegetation index and surface temperature are identified as pixels to be merged and merged into the current growing patch until no new pixels meet the merging conditions, forming a closed boundary of homogeneous vegetation patches, and finally generating a grassland homogeneous patch map.

4. The grassland drought identification and assessment method based on multi-source remote sensing image data according to claim 3, characterized in that, The determination of the similarity conditions satisfied by the pixels to be merged includes: The absolute value of the difference between the vegetation index value of the pixel to be merged and the average vegetation index value of the current growing patch is used to obtain the pixel to be merged. The absolute value of the difference between the surface temperature value of the pixel to be merged and the average surface temperature value of the current growing patch is calculated to obtain the surface temperature difference value. If the difference in vegetation index is less than a preset first threshold and the difference in surface temperature is less than a preset second threshold, then the pixel is determined to meet the merging conditions.

5. The grassland drought identification and assessment method based on multi-source remote sensing image data according to claim 3, characterized in that, Step S2, which assigns a drought sensitivity coefficient based on the community type in the grassland homogeneous patch map, determines the drought intensity value of each patch based on this coefficient, and then fits the local wet-dry edges to generate a grassland regional drought grid, also includes: Based on each patch in the grassland homogeneous patch map, extract the vegetation index value of all pixels corresponding to it in the grassland vegetation distribution data. The mean value of the vegetation index of all pixels in each patch is calculated as a parameter characterizing vegetation cover. Calculate the ratio of the near-infrared band to the visible light band in the region corresponding to the pure vegetation pixel within each patch, and determine the parameters characterizing the complexity of the canopy structure based on the ratio; Based on parameters characterizing vegetation cover and parameters characterizing canopy structure complexity, patches are classified into grassland community distribution types, and each patch is assigned a drought sensitivity coefficient ranging from 1 to 5 through a preset sensitivity determination matrix. Each patch is spatially correlated with its assigned drought sensitivity coefficient value to generate a patch stress sensitivity map.

6. The grassland drought identification and assessment method based on multi-source remote sensing image data according to claim 5, characterized in that, The determination of parameters characterizing the complexity of canopy structure also includes: Calculate vegetation cover for mixed pixels; The standard deviation of the normalized vegetation index within the mixed pixels was calculated using a 3×3 pixel window. When the standard deviation is greater than the preset difference threshold, for areas with vegetation coverage greater than 60%, the first derivative of the reflectance in the near-infrared band of the topographically corrected fused image is extracted. When the first derivative value is in the range of 0.02-0.05, it is identified as high canopy density grassland; when the first derivative value is less than 0.02, it is identified as low and sparse grassland.

7. The grassland drought identification and assessment method based on multi-source remote sensing image data according to claim 1, characterized in that, Step S2 involves assigning a drought sensitivity coefficient based on the community type in the grassland homogeneous patch map, determining the drought intensity value of each patch based on this coefficient, and then fitting the local wet-dry edges to generate a grassland regional drought grid, including: For each patch in the patch stress sensitivity map, the mean surface temperature and mean vegetation index are extracted from the grassland vegetation distribution data. The original drought degree of the current water stress state in the patch is calculated based on the average surface temperature and the average vegetation index. The original drought intensity is multiplied by the drought sensitivity coefficient corresponding to the patch in the patch stress sensitivity map to obtain the final drought intensity value of the patch. The local wet-dry edges of the target grassland region are fitted based on the drought intensity values, and the regions of the local wet-dry edges are closed and connected to generate a grassland zoning drought grid.

8. The grassland drought identification and assessment method based on multi-source remote sensing image data according to claim 2, characterized in that, The solar incidence angle correction factor for multi-source remote sensing images, determined based on slope and aspect angle, includes: Extract the solar azimuth angle of each pixel in the digital elevation model corresponding to the image imaging time from multi-source remote sensing images; Based on the range of slope values, the solar incidence angle correction factor is dynamically calculated using the slope aspect angle and solar azimuth angle.

9. A grassland drought identification and assessment system based on multi-source remote sensing image data, characterized in that, For executing the grassland drought identification and assessment method based on multi-source remote sensing image data as described in claim 1, the grassland drought identification and assessment system based on multi-source remote sensing image data includes: The terrain correction and fusion module is used to acquire multi-source remote sensing images and digital elevation models covering the target grassland area, and to perform terrain correction and band-weighted fusion on the multi-source remote sensing images based on the slope and aspect in the digital elevation model to obtain a terrain-corrected fused image. The drought zoning assessment module is used to identify and segment spatially continuous grassland homogeneous patches with similar vegetation growth and surface temperature based on topographically corrected fused images. It assigns a drought sensitivity coefficient based on the community type in the grassland homogeneous patch map, determines the drought intensity value of each patch based on the coefficient, and then fits the local dry and wet edges to generate a grassland zoning drought grid. The carrying capacity calculation module is used to determine the dynamic carrying capacity of each grid unit under the current drought stress based on the drought intensity value in the grassland zoning drought grid. The overgrazing monitoring and management module is used to determine the real-time carrying capacity by comparing the dynamic carrying capacity limit with the real-time distribution data of the herd transmitted by the Internet of Things devices. If the real-time distribution data of the herd exceeds the dynamic carrying capacity limit, it is determined to be an overgrazing state and dynamic carrying capacity management is implemented; otherwise, it is determined to be a safe carrying capacity state and the current grazing strategy is maintained.

Images