Method for automatic monitoring of regional desertification based on remote sensing cloud computing

Abstract

The application discloses a kind of based on remote sensing cloud computing's regional desertification automatic monitoring method, comprising: obtaining target area time series multispectral remote sensing image, extracting vegetation index and surface albedo;Based on composite ground object identification rule, eliminate non-target ground object pixel, construct basic sampling pixel set;Based on spatial homogeneity verification mechanism, sample extraction effective fitting sample set;The sample set is segmented and is statistically processed and is regressed, extracts two-dimensional characteristic space's dry side equation and wet side equation;Calculate the spatial geometric relative distance of each reserved pixel to dry and wet double sides, combined with the out-of-bound physical constraint truncation mechanism processing obtains relative desertification difference index;Based on the index, divide desertification grade.The present application overcomes the defect that extreme value fitting is fragile and threshold generalization fails, and realizes high-precision automatic monitoring.

Description

Technical Field

[0001] This invention relates to the field of remote sensing monitoring and geographic information technology, and in particular to an automatic monitoring method for regional desertification based on remote sensing cloud computing. Background Technology

[0002] Desertification monitoring is crucial for quantitatively characterizing surface water stress and ecological succession. Constructing a two-dimensional feature space using multispectral imagery, integrating vegetation indices and surface albedo, effectively decouples moisture signals from vegetation and soil background. Accurately retrieving boundary parameters of this feature space and extracting drought indices provides core technical support for the objective quantification of large-scale desertification spatial distribution.

[0003] Current spatial monitoring of desertification characteristics mostly adopts an offline processing paradigm. In a local computing environment, the absolute maximum and minimum values ​​of albedo within a predetermined vegetation index interval are extracted as extreme points. Ordinary least squares method is used to fit and generate dry-side equations and wet-side equations representing the extreme moisture state. Only the slope of the dry-side equation is used to construct the absolute difference index of a single boundary, and a fixed empirical numerical threshold is configured to discretize the desertification level of the target area.

[0004] However, existing technologies, when evolving towards large-scale fully automated batch processing, suffer from three major deep-seated technical defects: first, boundary fitting is extremely fragile; second, the exponent lacks cross-domain generalization ability; and third, the underlying computation is prone to algebraic deadlock. Extracting absolute extrema is susceptible to the leverage effect of heterogeneous noise such as residual thin clouds or saline-alkali land, and ordinary regression fails to distinguish confidence differences in sample density across local intervals, leading to a significant deviation of the fitting boundary from the true physical limit. The lack of wet-edge normalization constraints in unilateral absolute exponents causes their numerical domain to drift drastically with different climatic backgrounds, making it impossible to reuse fixed classification thresholds across time and space. In highly homogeneous terrain or high-vegetation convergence areas, the lack of adaptive sampling and overflow physical truncation mechanisms easily leads to missing or inverted feature boundaries, causing frequent system crashes due to division by zero in the denominator or overflow. Summary of the Invention

[0005] The purpose of this invention is to provide an automatic monitoring method for regional desertification based on remote sensing cloud computing, in order to solve at least one of the aforementioned problems in the existing technology.

[0006] Technical solution: A method for automatic monitoring of regional desertification based on remote sensing cloud computing, comprising:

[0007] Based on the acquired time-series multispectral remote sensing images of the target area, the vegetation index and surface albedo are extracted pixel by pixel and combined into a global sampling set;

[0008] Non-target ground object pixels in time-series multispectral remote sensing images are removed based on preset composite ground object identification rules to obtain a retained pixel sampling set;

[0009] Based on a preset spatial homogeneity verification mechanism, spatial sampling is performed on the retained pixel sampling set to extract an effective fitting sample set;

[0010] The effective fitted sample set is subjected to segmented statistical and regression processing to extract the dry edge equation and wet edge equation of the two-dimensional feature space composed of vegetation index and surface albedo.

[0011] Calculate the spatial geometric relative distances of each retained pixel to the dry and wet edges;

[0012] By combining the pre-set boundary-crossing physical constraint truncation mechanism with the spatial geometric relative distance, the relative desertification difference index of each retained pixel is obtained. Based on this, the desertification level of the target area is classified, and the desertification monitoring results are obtained.

[0013] Beneficial effects: This invention overcomes the defects of fragile extreme value fitting and threshold generalization failure, and achieves high-precision automated monitoring. Attached Figure Description

[0014] Figure 1 This is a schematic diagram of the overall process of an automatic monitoring method for regional desertification based on remote sensing cloud computing provided in an embodiment of this application.

[0015] Figure 2 This is a schematic diagram of the process of removing non-target land cover pixels from time-series multispectral remote sensing images based on preset composite land cover identification rules, and obtaining a basic sampled pixel set composed of vegetation index and surface albedo of the retained pixels, provided by an embodiment of this application.

[0016] Figure 3 This is a schematic diagram of the process for removing built-up area pixels based on preset built-up area identification rules provided in this application embodiment.

[0017] Figure 4 This is a schematic diagram of the NDVI-Albedo feature space provided in the embodiments of this application.

[0018] Figure 5 This is a scatter plot of the 1985 NDVI-Albedo feature space provided in the embodiments of this application. Detailed Implementation

[0019] Example 1 details the overall process of an automatic monitoring method for regional desertification based on remote sensing cloud computing, such as... Figure 1 As shown.

[0020] Step 101: Acquire time-series multispectral remote sensing images of the target area.

[0021] In this embodiment, the target area is the predetermined geographic spatial range for desertification monitoring. Data acquisition and processing are performed using a remote sensing cloud computing platform with massive distributed computing capabilities. A full-stack automated computation graph is constructed using the platform's server-side lazy evaluation mechanism, avoiding local downloading and transmission of intermediate data. Specifically, the time-series multispectral remote sensing imagery uses the 8th generation Landsat 2nd Edition Level 1 surface reflectance dataset, and is limited to continuous observation data within the vegetation growing season. After image acquisition, the system performs bitwise AND cloud masking based on the image's built-in pixel quality assessment bands, and removes abnormal observations indicating sensor saturation using the radiation saturation assessment bands, providing a clean observational base map for subsequent calculations.

[0022] Furthermore, to reduce the dimensionality of time-series multispectral remote sensing images into a single effective base map, the system employs the maximum vegetation index synthesis method to perform pixel-level synthesis of multiple observation data within the time series. A single-period composite image is generated by extracting the spectral reflectance of each pixel at a relatively vigorous vegetation growth point within the time series. This image serves as the base data for subsequent calculations, ensuring that the reflectance of each band in each pixel originates from the actual physical state.

[0023] In some alternative implementations, when the computing power of the target area is limited or it is not the vegetation growing season, the median composite method can be used as an alternative image composite strategy, that is, the median of the time series is extracted independently band by band to form a composite image.

[0024] In some alternative implementations, time-series multispectral remote sensing imagery can also be derived from Sentinel-2 satellite multispectral imagery, medium-resolution imaging spectrometer products, or other multispectral satellite data sources with continuous spectral coverage from visible to shortwave infrared. When using different satellite data sources, only the corresponding calibration scaling formulas and band weight coefficients in the empirical weighted model of surface albedo need to be adjusted according to the band configuration and radiometric calibration parameters of the data source; the rest of the data processing procedures remain consistent.

[0025] Step 102: Extract the vegetation index and surface albedo per pixel based on time-series multispectral remote sensing images.

[0026] In this embodiment, the original remote sensing image is typically stored as a digital quantized value. A pre-configured scaling factor needs to be applied to convert the quantized value into the true surface reflectance to avoid errors in subsequent calculations. The scaling formula is as follows:

[0027] ρ = DN * 0.0000275 - 0.2;

[0028] Where ρ is the true surface reflectance, DN is the original image pixel value, 0.0000275 is the pre-configured multiplication scaling factor, and 0.2 is the pre-configured additive offset.

[0029] After obtaining the true surface reflectance, vegetation indices are used to characterize the abundance and growth status of surface vegetation. Specifically, the following standard normalized vegetation index formula is used to process the near-infrared and red bands in time-series multispectral remote sensing images to calculate the vegetation index for each pixel:

[0030] NDVI=(ρ NIR -ρ Red ) / (ρ NIR +ρ Red );

[0031] Wherein, NDVI is the vegetation index, ρ NIR ρ represents the true surface reflectance in the near-infrared band. Red This represents the true surface reflectance in the red light band.

[0032] Meanwhile, surface albedo is used to characterize the energy balance of the Earth's surface. Specifically, an empirical weighted model for surface albedo is used to sum the actual surface reflectance in the blue, red, near-infrared, first shortwave infrared, and second shortwave infrared bands to obtain the surface albedo for each pixel. The formula for the empirical weighted model is as follows:

[0033] Albedo=0.356*ρ Blue +0.130*ρ Red +0.373*ρ NIR +0.085*ρ SWIR1 +0.072*ρ SWIR2 -0.0018;

[0034] Where Albedo is the surface albedo, ρ Blue ρ represents the true surface reflectance in the blue light band. Red ρ represents the true surface reflectance in the red band. NIR ρ represents the true surface reflectance in the near-infrared band. SWIR1 ρ represents the true surface reflectance in the first shortwave infrared band. SWIR2 The value represents the true surface reflectance of the second shortwave infrared band, while the other values ​​represent the band weight coefficients and fitting intercepts in the empirical weighted model of surface albedo.

[0035] Step 103: Based on the preset composite land cover identification rules, remove non-target land cover pixels from the time series multispectral remote sensing image to obtain a basic sampled pixel set composed of the vegetation index and surface albedo of the retained pixels. The basic sampled pixel set is the retained pixel sample set.

[0036] In other words, this step can also be to remove non-target land cover pixels from time-series multispectral remote sensing images based on preset composite land cover identification rules, and obtain retained pixels, with the vegetation index and surface albedo corresponding to the retained pixels forming a basic sampling pixel set.

[0037] The composite land cover identification rules are a set of logical criteria used to identify and filter out land cover types that interfere with the accuracy of desertification monitoring. Specifically, non-target land cover pixels mainly include water bodies, snow and ice, and built-up areas, which do not belong to the continuous physical succession gradient of vegetation and bare soil. By applying masking rules to remove noisy pixels with abnormally high albedo or low vegetation indices, pixels participating in subsequent feature space construction can have rigorous physical meaning. The pixels remaining after masking are the retained pixels. The system extracts the vegetation index and land albedo corresponding to all retained pixels and combines them into a basic sample pixel set for subsequent fitting. This set is specifically represented in the cloud system as a vector set structure containing numerous elements, each carrying clear spatial coordinate information and corresponding two-dimensional attribute values.

[0038] Step 104: Based on the preset spatial homogeneity verification mechanism, spatial sampling is performed on the retained pixel sampling set to extract the effective fitting sample set.

[0039] A spatial homogeneity verification mechanism is used to pre-quantify the richness of the ecological gradient within the target area. Since the fitting of the feature space is highly dependent on the complete transition zone from heavily bare soil to dense vegetation, if the area exhibits extreme homogeneity, fitting will fail. The system pre-verifies the attribute range of the base sampled pixel set. After the homogeneity verification passes or the corresponding compensation mechanism is triggered, adaptive random sampling is performed to extract a statistically representative set of effective fitted samples while controlling computational memory consumption.

[0040] Step 105: Perform piecewise statistical and regression processing on the effective fitted sample set to extract the dry edge equation and wet edge equation of the two-dimensional feature space composed of vegetation index and surface albedo.

[0041] The two-dimensional feature space is a scatter plotted space constructed with vegetation index as the horizontal axis and surface albedo as the vertical axis. In this space, the dry edge equation reflects the extreme drought state with the lowest surface water content under the same vegetation cover, while the wet edge equation represents the ideal state with sufficient soil moisture. The effective fitted sample set is segmented statistically, i.e., divided into multiple intervals with a fixed step size, and noise-resistant boundary feature points are extracted within each interval. Using the number of effective samples in each interval as weights, a weighted least squares regression with a constant term is performed to calculate the dry edge equation and wet edge equation representing the upper and lower boundaries. The automated cloud-based fitting mechanism avoids the reliance on manual point selection in traditional methods, improving the robustness of parameter calculation.

[0042] Step 106: Calculate the spatial geometric relative distance of each preserved pixel in the two-dimensional feature space with respect to the dry edge corresponding to the dry edge equation and the wet edge corresponding to the wet edge equation.

[0043] Spatial geometric relative distance is used to quantitatively measure the relative position of a predetermined pixel's true surface albedo between extreme drought and extreme wetness under current vegetation cover conditions. Specifically, based on the aforementioned dry-edge and wet-edge equations obtained as spatial benchmarks, the span from each preserved pixel to the dry edge and to the wet edge is determined. By using relative distance measurement, the dimensional distortion and incomparability of absolute values ​​under different regional environmental backgrounds are eliminated, giving this distance index a basis for cross-temporal and spatial generalization.

[0044] Step 107: Combine the preset cross-boundary physical constraint truncation mechanism to process the spatial geometric relative distance and obtain the relative desertification difference index of each retained pixel.

[0045] The out-of-bounds physical constraint truncation mechanism is a mandatory error correction fallback logic designed to address the outlier overflow inevitably caused by statistical boundary fitting and the risk of dividing the mathematical denominator by zero. In actual calculations, the surface albedo of some retained pixels may exceed the theoretical dry-wet envelope. Through this mechanism, the system identifies these outliers that exceed the limits of physical norms, as well as pixels with equation crossover failures, and forcibly assigns or masks off their spatial geometric relative distances. After the above physical constraint correction, the final output value is the relative desertification difference index strictly normalized to the standard interval.

[0046] Step 108: Based on the relative desertification difference index of each preserved pixel, the desertification level of the target area is classified to obtain the desertification monitoring results.

[0047] Desertification monitoring results are presented as raster imagery with a clear spatial projection system and discrete numerical classifications. Specifically, since the relative desertification difference index has been uniformly mapped to a fixed numerical range, the system can directly divide it into multiple standard desertification levels according to a preset equal-interval step size. Through quantitative mapping and division, the system generates a final operational output that intuitively reflects the spatial distribution pattern of desertification in the region.

[0048] It should be noted that the above order can be adjusted. For example, non-target feature pixels can be removed first, and then the vegetation index and surface albedo of the retained pixels can be extracted to form the basic sampled pixel set. In other words, non-target feature pixels in time-series multispectral remote sensing images are removed based on preset composite feature identification rules, and the data corresponding to the retained pixels are selected from the extracted pixel-by-pixel vegetation index and surface albedo to obtain the basic sampled pixel set. The step numbers are not used to limit the execution flow, but are only for the convenience of description.

[0049] According to one aspect of this application, a method for automatic monitoring of regional desertification based on remote sensing cloud computing may further include the following steps:

[0050] Based on the acquired time-series multispectral remote sensing images of the target area, the vegetation index and surface albedo are extracted pixel by pixel and combined into a global sampling set;

[0051] Non-target ground object pixels in time-series multispectral remote sensing images are removed based on preset composite ground object identification rules to obtain a retained pixel sampling set;

[0052] Based on a preset spatial homogeneity verification mechanism, spatial sampling is performed on the retained pixel sampling set to extract an effective fitting sample set;

[0053] The effective fitted sample set is subjected to segmented statistical and regression processing to extract the dry edge equation and wet edge equation of the two-dimensional feature space composed of vegetation index and surface albedo.

[0054] Calculate the spatial geometric relative distances of each retained pixel to the dry and wet edges;

[0055] By combining the pre-set boundary-crossing physical constraint truncation mechanism with the spatial geometric relative distance, the relative desertification difference index of each retained pixel is obtained. Based on this, the desertification level of the target area is classified, and the desertification monitoring results are obtained.

[0056] Example 2 details how, before constructing a two-dimensional feature space, various non-target ground object pixels that interfere with the physical meaning of the desertification feature space can be accurately identified and eliminated through predetermined band combinations and external high-precision data sources. Figure 2 As shown.

[0057] Step 201: Non-target ground feature pixels include water body pixels, ice and snow pixels, and built-up area pixels.

[0058] In this embodiment, the premise for constructing a two-dimensional feature space of vegetation index and surface albedo is that the pixels involved in the fitting must be able to represent a continuous ecological succession gradient from completely bare soil to higher vegetation cover. Non-target land cover pixels refer to land cover types that exhibit abnormal reflectance in spectral features and whose physical properties are unrelated to desertification water stress processes.

[0059] Specifically, water body pixels typically exhibit low near-infrared reflectance and surface albedo; if not removed, they will lower the wet edge baseline in the feature space. Ice and snow pixels have high visible light reflectance, causing an abnormally positive shift in the slope of the dry edge equation. Built-up area pixels, due to the high albedo of artificial impermeable surfaces, are easily confused with severely desertified, arid, bare soil. By explicitly defining these three typical interference sources as non-target ground feature pixels, a foundation is laid for subsequent high-precision composite masking operations.

[0060] Step 202: Extract the green band reflectance, shortwave infrared 1 band reflectance, and shortwave infrared 2 band reflectance from the time-series multispectral remote sensing image.

[0061] The time-series multispectral remote sensing imagery contains multiple optical observation channels within predetermined wavelength ranges. The system selectively extracts corresponding band data for exponential calculations based on the differences in spectral absorption and reflectance between water bodies and ice / snow. Specifically, taking the eighth-generation Landsat dataset as an example, the green band reflectance corresponds to the third observation channel with a wavelength range of 0.525 μm to 0.600 μm. The shortwave infrared 1 band reflectance corresponds to the sixth observation channel with a wavelength range of 1.560 μm to 1.660 μm. The shortwave infrared 2 band reflectance corresponds to the seventh observation channel with a wavelength range of 2.100 μm to 2.300 μm. All of the above band reflectance data have undergone radiometric calibration and atmospheric correction in previous steps, representing the true physical reflectance characteristics of the Earth's surface.

[0062] Step 203: Calculate the corrected normalized water index based on the green light band reflectance and the shortwave infrared 1 band reflectance, and remove water pixels by combining the preset water threshold.

[0063] To accurately identify surface water bodies, an existing index model with enhanced sensitivity for water body identification is employed. This model leverages the spectral difference between water bodies, which exhibit high reflectivity in the green light band and strong absorption characteristics in the shortwave infrared band. Specifically, based on the reflectivity in the green light band and the shortwave infrared band 1, the standard modified normalized water index formula is used to calculate the index value for each pixel. The formula for calculating the modified normalized water index is as follows:

[0064] MNDWI=(ρ Green -ρ SWIR1 ) / (ρ Green +ρ SWIR1 );

[0065] Where MNDWI is the modified normalized water index, ρ Green ρ is the reflectivity in the green light band. SWIR1 The reflectance is for the shortwave infrared band 1.

[0066] The value is compared pixel-by-pixel with a preset water threshold. The preset water threshold can be set to 0. Pixels with a value greater than 0 are identified as water pixels and are masked out from the computation process.

[0067] Step 204: Calculate the improved normalized snow index based on the green light band reflectance and the shortwave infrared 2 band reflectance, and remove ice and snow pixels by combining the preset ice and snow threshold.

[0068] For seasonal snow cover or alpine glaciers commonly found in high-altitude areas of arid and semi-arid regions, the system employs an improved normalized snow index with fine-tuned band combinations. Traditional snow indices typically use the green light band and shortwave infrared band 1; in this embodiment, the longer-wavelength shortwave infrared band 2 is selected to replace shortwave infrared band 1. The absorption depth of ice crystal particles near 2.1 micrometers is significantly greater than their absorption depth near 1.6 micrometers, allowing snow and ice pixels to obtain a larger positive response value in the normalized difference model constructed based on shortwave infrared band 2. The calculation formula for the improved normalized snow index is as follows:

[0069] NDSI star =(ρ Green -ρ SWIR2 ) / (ρ Green +ρ SWIR2 );

[0070] Among them, NDSI star For the improved normalized snow cover index, ρ Green ρ is the reflectivity in the green light band. SWIR2 The reflectance is for the shortwave infrared band 2.

[0071] After calculation, the obtained improved normalized snow cover index is compared with a preset snow and ice threshold. The preset snow and ice threshold can be set to 0.4. The system classifies pixels with a calculation result greater than 0.4 as snow and ice pixels and removes them from the subsequent feature space construction dataset. This step effectively prevents snow and ice with high albedo from being mistaken for extremely arid bare soil or saline-alkali land, ensuring the scientific validity of the dry edge equation fitting.

[0072] Step 205: Remove built-up area pixels based on preset built-up area identification rules.

[0073] Because there is a high degree of overlap in optical reflectivity between man-made structures and natural bare soil, this embodiment introduces a more flexible mechanism to process built-up area pixels. Built-up area identification rules are automatically invoked by the system or manually pre-configured based on the specific drought level and surface background characteristics of the target area, adapting to the monitoring needs of different regions.

[0074] In a preferred embodiment, for extremely arid areas with sparse vegetation and large areas of exposed bare soil, built-up area pixels are removed based on preset built-up area identification rules, such as... Figure 3 As shown, this is achieved in the following way:

[0075] Step 2051: Obtain the pre-configured global land cover classification product.

[0076] In extremely arid regions, natural sandy land and gravelly Gobi deserts are almost indistinguishable from urban concrete and asphalt roads using conventional spectral indices. Forcibly using spectral indices often leads to the most severely desertified areas being mistakenly identified as urban areas. This preferred embodiment directly obtains publicly available third-party land cover products with high interpretation accuracy. For example, a pre-configured global land cover classification product can be obtained from the 10-meter spatial resolution world land cover dataset published by the European Space Agency.

[0077] Step 2052: Extract predefined category pixels from the pre-configured global land cover classification product based on preset category labels.

[0078] The system extracts a set of pixels with clear semantic attributes based on preset category labels. For example, when using the European Space Agency's data products, the numerical code representing built-up areas in its internal data dictionary is 50. Using this value as the preset category label, the corresponding grid is extracted to form pixels of the predetermined category.

[0079] Step 2053: Pixels in the time-series multispectral remote sensing image whose spatial location overlaps with that of a predetermined category of pixels are removed as built-up area pixels.

[0080] The system performs spatial projection alignment and resampling operations, spatially overlaying the currently processed time-series multispectral remote sensing image with the mask layer containing pixels of predetermined categories. Image pixels falling within the predetermined category's pixel range on the geospatial coordinates are marked as built-up area pixels and removed using a hard mask. This direct import of spatial prior knowledge solves the technical problem of confusing artificial structures with severely desertified bare soil in extremely arid areas.

[0081] In addition, in another optional implementation, for semi-arid or semi-humid areas with good vegetation cover, the removal of built-up area pixels based on preset built-up area identification rules is achieved in the following way:

[0082] Step 2054: Extract near-infrared reflectance based on time-series multispectral remote sensing images.

[0083] The system extracts predetermined near-infrared band true surface reflectance data from time-series multispectral remote sensing images that have undergone cloud masking and radiometric calibration. For the eighth-generation Landsat dataset, the near-infrared band reflectance corresponds to the fifth band observation channel with a wavelength range between 0.845 micrometers and 0.885 micrometers.

[0084] Step 2055: Calculate the normalized building index based on the shortwave infrared band 1 reflectance and the near-infrared band reflectance.

[0085] The system reuses the shortwave infrared band 1 reflectance extracted in the previous steps, and combines it with the near-infrared band reflectance, directly using the existing normalized building index formula for calculation. The calculation formula is as follows:

[0086] NDBI=(ρ SWIR1 -ρ NIR ) / (ρ SWIR1 +ρ NIR );

[0087] Wherein, NDBI is the Normalized Building Index, ρ SWIR1 ρ is the reflectivity of the shortwave infrared band 1. NIR This refers to the reflectivity in the near-infrared band.

[0088] Because in areas with relatively dense vegetation, the short-wave infrared reflectance of buildings is significantly higher than their near-infrared reflectance, while the surrounding vegetation exhibits the opposite characteristics, this classical spectral index is sufficient to distinguish non-target features.

[0089] Step 2056: Determine the normalized building index based on the preset building threshold and remove built-up area pixels.

[0090] The system compares the calculated normalized building index with a preset building threshold. The preset building threshold can be set to 0. Pixels with a value greater than 0 will be identified as built-up area pixels and removed by the system. This solution enables end-to-end data looping during cloud computing, providing an effective parallel alternative when external high-precision classification products are not required.

[0091] Step 206: The remaining pixels after the removal process are taken as retained pixels to obtain the retained pixel sampling set.

[0092] After sequentially or in parallel identifying and masking water body pixels, snow and ice pixels, and built-up area pixels, contaminated or disturbed non-target features in the time-series multispectral remote sensing imagery have been filtered out. The system uniformly marks the remaining clean pixels not covered by the mask as retained pixels. These pixels accurately and rigorously reflect the physical succession gradient of the interaction between surface vegetation and soil moisture. The system extracts and summarizes the spatial geometric coordinates, vegetation index values, and surface albedo values ​​of all retained pixels to jointly construct the basic sampled pixel set used to support the subsequent two-dimensional feature space core fitting algorithm.

[0093] Example 3 elaborates on how, before extracting feature space fitting samples, ecological heterogeneity of the target region is pre-examined, and adaptive large-scale random sampling is performed in normal regions, or a buffer resampling mechanism is triggered to extend outward in extremely homogeneous regions.

[0094] Step 301: Based on the preset spatial homogeneity verification mechanism, spatial sampling is performed on the retained pixel sampling set to extract the effective fitting sample set.

[0095] In this embodiment, the scientific fitting of the dry and wet dual boundaries of the two-dimensional feature space is highly dependent on the objective physical laws that a complete transition zone from heavily bare soil to dense vegetation must exist simultaneously within the target area. If the input target area is extremely small, or located in a region with strong homogeneous surface features such as a purely natural desert environment, the spectral features of all pixels within it will exhibit extreme homogeneity. In this case, forcibly performing binning regression on the pixels within the region will lead to a mathematical collapse of the feature space, causing system computation to crash or output absurd boundaries. Before performing formal spatial sampling, the system introduces a spatial homogeneity verification mechanism as a safety valve, providing a decision basis for subsequently adaptively extracting a valid fitting sample set with strong statistical representativeness.

[0096] Step 302: Extract the maximum and minimum vegetation index values ​​of the retained pixels in the basic sampled pixel set.

[0097] A global attribute traversal is performed on the base sampled pixel set after preliminary composite masking purification. Specifically, using the distributed statistical reduction operator of the remote sensing cloud computing platform, the maximum vegetation index representing the densest vegetation cover in the current area and the minimum vegetation index representing the most extremely bare surface are retrieved and extracted from all retained pixels. These extreme value data objectively reflect the upper and lower limits of vegetation cover at the current monitoring scale. In practice, to avoid extreme value artifacts caused by a few residual noise points, the 99th percentile can be used to replace the absolute maximum vegetation index, and the 1st percentile can be used to replace the absolute minimum vegetation index, thereby improving the robustness of the statistics.

[0098] Step 303: Calculate the difference between the maximum and minimum vegetation index values ​​to obtain the vegetation index range.

[0099] The macro-ecological heterogeneity within the current target region is quantitatively assessed by taking the algebraic difference between the two extracted global extrema. The calculation formula is as follows:

[0100] ΔNDVI=NDVI max -NDVI min ;

[0101] Where ΔNDVI is the vegetation index range, NDVI max The maximum value of the vegetation index, NDVI min This represents the minimum value of the vegetation index.

[0102] Based on this, the vegetation index range is used as the core hub for subsequent determination of whether the buffer fallback mechanism is triggered, which directly and sensitively reflects the tolerance of vegetation gradient in the target area.

[0103] Step 304: Based on the comparison results between the vegetation index range and the preset homogeneity threshold, execute the corresponding spatial sampling strategy to extract the effective fitted sample set.

[0104] The system acquires a pre-configured judgment criterion, which involves comparing the preset homogeneity threshold with the calculated vegetation index range. The preset homogeneity threshold represents the minimum vegetation cover gradient that the system can tolerate in fitting a statistically significant two-dimensional feature space. The preset homogeneity threshold can be set to 0.15. When the vegetation index range is greater than or equal to the preset homogeneity threshold, it indicates that the region exhibits a reasonable succession from bare soil to vegetation; conversely, it indicates that the region is in an abnormal state of high homogeneity. Based on the mutual exclusivity of the comparison results, the system will automatically switch to and execute one of two parallel spatial sampling strategies to ensure that the effective fitting sample set ultimately delivered to the fitting engine always possesses sufficient feature diversity.

[0105] Furthermore, in a typical implementation for heterogeneous landforms, based on the comparison between the vegetation index range and a preset homogeneity threshold, a corresponding spatial sampling strategy is executed to extract an effective fitted sample set. This is achieved in the following way:

[0106] Step 305: If the vegetation index range is greater than or equal to the homogeneity threshold, then random sampling is performed on the retained pixel sampling set based on the pre-configured sampling size parameter to extract the effective fitted sample set.

[0107] When the system determines that the vegetation index range meets or exceeds a preset homogeneity threshold, it indicates that the target area itself possesses an extremely rich ecological gradient. For example, in a typical grassland area, the system calculates that its maximum vegetation index is 0.85 and its minimum vegetation index is 0.10, corresponding to a vegetation index range of 0.75, which is much greater than 0.15. In this case, the system does not need to rely on external data and directly performs large-scale random sampling within the current basic sampled pixel set.

[0108] Specifically, to prevent memory overflow exceptions when performing global reduction on tens of millions of pixels in the cloud, a pre-configured sampling size parameter is used to scale and truncate the overall data. This pre-configured sampling size parameter can be set to 50,000 pixels. The system calls a stratified random sampling algorithm to strictly limit the total sample size within a safe range, while simultaneously declaring a memory cap coefficient for expanding distributed slicing operations to the underlying operators. Through these parameter configurations, the system randomly extracts and combines samples from a massive amount of raw pixels to form an effective fitted sample set, ensuring the representativeness of the fitted samples while avoiding the computing power bottleneck of the cloud computing platform.

[0109] Furthermore, in an alternative implementation, when the computing power of the cloud computing platform is extremely limited, the system can use fixed random point sampling with a lower sample size as an alternative solution. In this solution, the system directly generates one thousand fixed spatial coordinate points randomly within the target area, and extracts only the vegetation index and surface albedo at these points to form an effective fitting sample set, achieving rapid and coarse monitoring while sacrificing some fitting robustness.

[0110] In another fault-tolerant implementation for extremely homogeneous terrain, based on the comparison between the vegetation index range and a preset homogeneity threshold, a corresponding spatial sampling strategy is executed to extract an effective fitted sample set. This is achieved in the following way:

[0111] Step 306: If the vegetation index range is less than the homogeneity threshold, the spatial boundary of the target area is extended outward based on the preset buffer extension distance to generate an extended buffer area.

[0112] When the system determines that the vegetation index range fails to reach the preset homogeneity threshold, the system will immediately trigger a safety fallback mechanism to prevent feature space collapse. For example, if the input target area is a very small block located in the heart of a desert, the system calculates that the vegetation index range of its internal pixels is only 0.04, far less than 0.15.

[0113] At this point, the system automatically acquires the preset buffer extension distance and, based on the spatial vector boundary of the original target area, performs polygonal expansion buffering in all outward directions. The preset buffer extension distance can be set to fifty kilometers. This geometrical extension operation forcibly includes the geographical environment surrounding the target area, artificially reconstructing the lost ecological gradient on a macro scale, forming a broader extended buffer zone.

[0114] Step 307: Extract the vegetation index and surface albedo of the pixels located in the extended buffer area based on the time series multispectral remote sensing image. Remove non-target feature pixels based on the composite feature identification rules. Construct a candidate buffer pixel set by combining the vegetation index and surface albedo of the remaining pixels.

[0115] After obtaining the outer buffer zone, it is essential to ensure that the newly introduced surrounding environmental pixels maintain absolute consistency with the original core area pixels in terms of physical attributes and data quality. Therefore, the system directly reuses the underlying data source of the time-series multispectral remote sensing image and only extracts pixels falling within the spatial range of the outer buffer zone.

[0116] The system rigorously reuses all the purification logic defined above, and based on the same composite feature recognition rules, performs strict masking operations on the newly extracted extension pixels, removing all interfering non-target feature pixels. After purification and filtering to the same standard as the core area, the vegetation index and surface albedo of the remaining high-quality pure extension pixels are integrated by the system to construct a candidate buffer pixel set independent of the original set.

[0117] Step 308: Merge the candidate buffer cell set with the retained cell sampling set, and randomly sample the merged cell set based on the pre-configured sampling size parameter to extract the effective fitted sample set.

[0118] The system seamlessly merges and stitches together a set of candidate buffer pixels rich in gradient information with the original, highly homogeneous set of base sample pixels in cloud memory. Through the borrowing and compensation of spatial scale, the merged new pixel pool breaks free from the extreme spectral constraints of a single landform, forming a complete succession sequence from severe desertification to healthy vegetation.

[0119] After completing the merging and reconstruction, the system again applies the aforementioned pre-configured sampling scale parameters to perform large-cardinality random sampling on the large and heterogeneous merged pixel pool, ultimately extracting an effective set of fitted samples that can support the normal operation of the underlying regression engine. Through a multi-track parallel adaptive sampling design, this scheme breaks through the stringent constraints of traditional feature space models on the natural topographical environment of the target area.

[0120] Example 4 details how, after obtaining effective fitted samples, a series of automated mathematical processes on the cloud server, such as equal-width binning, percentile noise-resistant extraction, adaptive sample filtering, and underlying matrix reconstruction, are used to solve the highly robust two-dimensional feature space boundary equation.

[0121] Step 401: Perform piecewise statistical and regression processing on the effective fitted sample set to extract the dry edge equation and wet edge equation of the two-dimensional feature space composed of vegetation index and surface albedo.

[0122] In traditional desertification remote sensing monitoring methods, the dry and wet boundaries in a two-dimensional feature space often rely on researchers exporting data locally and manually fitting the data using extreme value lines from scatter plots or subjective bounding boxes. This operation not only disrupts the automated workflow but also makes manual extreme values ​​susceptible to interference from high-reflectivity anomalies such as residual thin clouds and saline-alkali land, causing the fitted line to deviate significantly from the true physical boundary. This invention constructs a black-box bilateral parameter calculation engine on a cloud computing platform server that requires no manual intervention. By transforming the discrete distribution structure of all valid fitted sample sets into a standardized mathematical matrix system, it performs highly robust piecewise statistics and adaptive weighted solving, accurately outputting the dry-side equation representing the region's extreme water stress state and the wet-side equation representing the optimal water supply state.

[0123] Step 402: Divide the vegetation indices in the effective fitted sample set into intervals based on a preset partitioning step size to obtain multiple vegetation index boxes.

[0124] The system first obtains a preset partitioning step size and performs a regular spatial slicing operation in the two-dimensional feature space, using the vegetation index as the independent variable on the horizontal axis. The preset partitioning step size determines the precision of subsequent searches for local extrema or percentiles on the horizontal axis. In actual system configurations, the preset partitioning step size can be set to 0.01. Based on the minimum and maximum values ​​of the vegetation index within the extracted effective fitted sample set, a continuous equidistant scale is constructed using a cloud-based sequence generation function, dividing the entire continuous vegetation index span into dozens or hundreds of non-overlapping sub-intervals. Each independent sub-interval is instantiated in the system as a structured data aggregation unit, i.e., a vegetation index bin. Through this operation, the originally disordered distribution of massive scattered points is systematically placed into the corresponding vegetation index bins, providing a discretized basic framework for subsequent local feature extraction.

[0125] Step 403: Count the number of pixels in each vegetation index box to obtain the number of samples in the box.

[0126] After completing the basic binning and placement, the system initiates a local statistical reduction process for each vegetation index bin. Specifically, the system iterates through all the feature objects contained within each independent vegetation index bin and calculates the length of its feature set, i.e., the number of all pixels falling within the predetermined vegetation index span. The statistically obtained value is the sample count within the bin. In natural ecological succession systems, the vast majority of pixels are concentrated in bins representing moderate or dominant ground cover types, while edge bins in extreme desertification or extreme wetland states often contain only extremely sparse pixels. The sample count within a bin is not just a simple count value; it is a core benchmark representing the reliability and confidence level of the statistical characteristics of each local bin in subsequent weighted regression calculations.

[0127] Step 404: Extract boundary features for each vegetation index box and generate corresponding dry edge candidate points and wet edge candidate points.

[0128] For each completed vegetation index box, the system needs to search for points on the vertical axis, i.e., the direction of surface albedo, that can represent the physical meaning of extreme drought and extreme wetness under that local vegetation cover. These extracted key scattered points constitute the dry-edge candidate points and wet-edge candidate points used for subsequent linear fitting trajectory anchoring. Connecting the discrete candidate point sets together initially outlines the nonlinear envelope contour of the two-dimensional feature space.

[0129] Furthermore, addressing the technical weakness of traditional line fitting based on interval maximum and minimum values, which is prone to failure, the following method is used to extract boundary features for each vegetation index box and generate corresponding dry and wet edge candidate points:

[0130] Step 405: Extract the 95th percentile and 5th percentile of the surface albedo of each pixel within the vegetation index box.

[0131] The system does not extract absolute extreme values, but employs percentile statistics, which are highly resistant to random noise. Within each independent vegetation index box, the surface albedo of all pixels is sorted and distributed according to their numerical values. To mask the heterogeneous noise at the tail end of the distribution, such as abnormally low albedo due to residual shadows or abnormally high albedo due to incompletely masked salt crusts, values ​​at predetermined positions in the distribution sequence are directly extracted. Specifically, the values ​​at the 95th percentile, arranged from smallest to largest, are extracted as the 95th percentile, and the values ​​at the 5th percentile are extracted as the 5th percentile.

[0132] Step 406: Extract the center value of vegetation index for each vegetation index box.

[0133] After anchoring the percentiles on the ordinate, a unique representative coordinate for the horizontal axis needs to be determined for this vegetation index box. The lower limit of the vegetation index interval for this box is obtained, and half of the preset division step size is added to calculate the central value of the vegetation index at the exact center of the interval. By employing strict center alignment positioning, the horizontal axis position offset error that may be caused by the skewed distribution of samples within the box is eliminated.

[0134] Step 407: Combine the central value of the vegetation index with the corresponding 95th percentile to form a candidate dry edge point.

[0135] The obtained vegetation index center value is used as the x-axis component, and the previously extracted 95th percentile with noise resistance is used as the y-axis component. Both are assembled into a two-dimensional coordinate point object, namely, a dry edge candidate point, in the cloud data structure. This candidate point physically represents the robust extreme drought state under the current predetermined vegetation cover level, excluding occasional outlier disturbances.

[0136] Step 408: Combine the central value of the vegetation index with the corresponding 5th percentile to form wet edge candidate points.

[0137] Similarly, the system again uses the central value of the vegetation index as the abscissa component and the corresponding 5th percentile as the ordinate component to assemble another two-dimensional coordinate point object, namely the wet edge candidate point. This candidate point represents the robust optimal water supply state after excluding the interference of abnormal dark target shadows.

[0138] Step 409: Using the number of samples within each vegetation index box as the regression weight, perform weighted regression processing on the dry edge candidate points and wet edge candidate points respectively to obtain the dry edge equation and wet edge equation.

[0139] After acquiring all robust boundary point pairs contributed by each bin, the system uses a linear function shape to optimally approximate the envelope trend based on statistical principles. Traditional feature space fitting not only uses extreme points but also typically performs indiscriminate ordinary least squares regression. This results in extremely sparse bin candidate points, containing only two or three random pixels, located at the edge of the vegetation gradient, having the same leverage effect when pulling the fitted straight line as dense bin candidate points, located in the ecological core zone, containing tens of thousands of pixels, which is statistically dangerous. This embodiment uses the number of samples within each bin obtained from previous statistics as a forced constraint, giving densely packed bins a greater straight-line traction force, so that the final output dry-edge equation and wet-edge equation can fit the vast majority of high-confidence core point pair sequences to the greatest extent.

[0140] Furthermore, to address the percentile degradation problem caused by lower sample sizes and to correct the numerical singularity and boundary overflow defects caused by the forced lack of intercept constraint in the linear regression operator of the cloud computing platform at the underlying mathematical structure, weighted regression processing was performed on dry and wet candidate points using the number of samples within each vegetation index box as the regression weight. This yielded the dry and wet equations, which were achieved through the following adaptive modified algebra construction strategy:

[0141] Step 410: Count the total number of pixels in the effective fitted sample set to obtain the total sample size, and count the total number of vegetation index boxes to obtain the total number of boxes. Based on the total sample size of the effective fitted sample set and the total number of vegetation index boxes, calculate the adaptive sample lower limit threshold.

[0142] The system avoids the generalization dilemma caused by hard-coding fixed values ​​and dynamically infers a reasonable in-bin elimination baseline based on the overall abundance of the current input dataset. If a fixed constant is forcibly used, such as requiring all bins to have more than 10 pixels, this constant will either lead to the complete elimination of valid samples or render the noise reduction meaningless due to an excessively low threshold when dealing with high-resolution small regions or low-resolution large regions. The system calculates the adaptive sample lower limit threshold based on the cumulative sum of all pixels in the valid fitted sample set (i.e., the total sample size) and the number of vegetation index bins (i.e., the total number of bins), using a dynamic lower limit inference formula.

[0143] n limit =max(10,N total / (k*50));

[0144] Where, n limit The calculated adaptive sample lower limit threshold is given by N, where max is a function operator that takes the larger of the two values. total To effectively fit the sample set, the total sample size is defined by k, which is the total number of vegetation index boxes currently divided, 50 is the preset empirical dilution coefficient constant, and 10 is the minimum hard bottom line to prevent percentile statistics from failing.

[0145] Step 411: Remove vegetation index boxes and their corresponding dry and wet edge candidate points if the number of samples in the box is less than the adaptive sample lower limit threshold, and extract the set of retained boxes.

[0146] The system iterates through all vegetation index boxes, comparing the number of samples within each box with the dynamically calculated adaptive lower limit threshold. Boxes with a sample count below this dynamic threshold are deemed to have excessively sparse samples, whose calculated percentiles are not representative and are at risk of degenerating into one-sided extrema. The system directly and forcibly removes and destroys these unqualified boxes, along with their generated dry-edge and wet-edge candidate points. After this rigorous adaptive purification and filtering mechanism, all box units with sufficient sample support constitute a set of retained boxes with a high level of statistical confidence.

[0147] Step 412 introduces a constant term column vector to construct a joint independent variable matrix for the retained box set.

[0148] In the servers of mainstream remote sensing cloud computing platforms, the built-in linear regression operator, when receiving only a single independent variable sequence—for example, only inputting a column of vegetation index center values—is forced to limit its underlying matrix operation logic to solving a direct proportional function model passing through the origin and excluding the intercept. This means it forcibly fits Y=a'X, where Y is the dependent variable, X is the independent variable, and a' is the slope in the linear regression model. However, the actual wet-dry boundary does not always pass through the origin of the coordinate system. If this underlying constraint is not forcibly broken during the algebraic input stage, when the system subsequently attempts to extract the non-existent intercept parameter, it will immediately trigger memory errors and deadlocks such as array out-of-bounds errors.

[0149] To legally solve the complete linear equation, before the system inputs the data into the regression engine, a one-dimensional array structure—the constant term column vector—is manually constructed in cloud memory. This array structure has the same length as the number of elements in the retained bin set, and all its elements have a constant value of 1. The system uses this constant term column vector as the zeroth dimension feature on the horizontal axis and merges it vertically with the vegetation index center value sequence, which serves as the first dimension feature. This generates a two-dimensional joint independent variable matrix in memory, containing two columns of independent variable features.

[0150] Step 413: Using the joint independent variable matrix as the independent variable feature and the number of samples in the bin corresponding to the retained bin set as the regression weight, perform weighted least squares regression to obtain the dry side equation and wet side equation containing the intercept parameter.

[0151] After completing the reconstruction matrix of the independent variable features, the system sends a final computation request to the underlying regression engine. At this point, the system uses the joint independent variable matrix as the input to the design matrix and the albedo values ​​of the corresponding candidate points as the input to the target matrix. The system explicitly specifies the number of samples within each bin as the external weighting array. The underlying operator executes a weighted matrix partial derivative solving algorithm. Since the independent variable design matrix contains a constant term column vector consisting entirely of 1s, the matrix operation fully allocates the systematic constant deviation of the dependent variable to the regression coefficients of the constant column vector. This ensures that the operator safely and legally returns an array containing two coefficient results. The system extracts the slope parameter corresponding to the first dimension and the intercept parameter corresponding to the zeroth dimension from this array, constructing complete dry-side and wet-side equations that accurately characterize the physical properties of extreme moisture values ​​in the current remote sensing image under different vegetation covers.

[0152] Example 5 elaborates on how, after obtaining high-precision dry and wet edge equations, the absolute coordinate positions of pixels in the two-dimensional feature space are transformed into cross-regional universal spatial geometric relative distances that can offset environmental background differences by constructing multi-track parallel geometric and algebraic ranging logic. At the same time, the absolute difference index method as a downgrade alternative is disclosed.

[0153] Step 501: Calculate the spatial geometric relative distance of each preserved pixel in the two-dimensional feature space relative to the dry edge corresponding to the dry edge equation and the wet edge corresponding to the wet edge equation. That is, calculate the spatial geometric relative distance of each preserved pixel to the dry edge and the wet edge.

[0154] After obtaining the dry and wet dual boundaries describing the baseline of the current study area's limiting moisture conditions, this embodiment abandons the traditional approach of calculating the absolute numerical index of a single side solely using the slope of the dry side. Because the physical distance between the dry and wet sides varies drastically in the two-dimensional feature space of different climate zones and years, the numerical range of the single-side absolute index lacks comparability. This embodiment constructs a more rigorous bilateral constraint framework. By measuring the physical intervals of each preserved pixel to the two edges of the drought and wet limits in a scattered coordinate system, and by constructing a ratio model, the dimensional and absolute value distortions caused by different environmental backgrounds in the study area are eliminated. Finally, a scalar theoretically normalized to a predetermined interval is obtained, which is the spatial geometric relative distance.

[0155] Step 502: Calculate the first distance from each retained pixel to the dry edge.

[0156] The system iterates through all retained pixels involved in the calculation, using the previously fitted dry edge equation as a rigid mathematical benchmark. For each predetermined pixel coordinate, the system quantifies the physical distance difference between the pixel's current location and its most severe drought condition by calling a preset spatial ranging model, obtaining a first distance characterizing its deviation from the drought limit. The smaller this distance, the closer the pixel's true surface albedo is to the theoretical limit under that vegetation cover condition, indicating a more severe degree of desertification.

[0157] Step 503: Calculate the second distance from each retained pixel to the wet edge.

[0158] Similarly, the system takes the current coordinates of the pixel as the starting point and the wet edge corresponding to the wet edge equation representing the most humid state as the ending point. It then performs the predetermined ranging model calculation again to quantify the interval between the current position and the theoretically most humid state, obtaining the second distance. The magnitude of this distance directly reflects the gap between the pixel's current moisture supply and its potential maximum.

[0159] Step 504: Calculate the sum of the first distance and the second distance to obtain the total dry and wet distance.

[0160] After extracting the independent measures from the pixel to the two side boundaries, the first distance and the second distance are algebraically summed. This operation reconstructs the total dry-wet tolerance span of the entire two-dimensional feature space at the predetermined vegetation index position of the pixel, which is the total dry-wet distance. The total dry-wet distance is a dynamic variable that tends to converge or diverge as the vegetation index on the horizontal axis increases. This non-parallel attribute reflects the differentiated sensitivity of different vegetation cover to changes in soil background moisture.

[0161] Step 505: Calculate the ratio of the first distance to the total dry and wet distance to obtain the spatial geometric relative distance.

[0162] After obtaining the absolute value of a pixel's deviation from the drought limit and the system's overall tolerance at that location, the system combines the two into a dimensionless fractional structure. Specifically, using the first distance as the numerator parameter and the total dry-wet distance as the denominator scale benchmark, the spatial geometric relative distance is obtained by dividing the two. Through relative ratio calculation, the drought level of a pixel is naturally projected and mapped to a relative range of 0 to 1, providing a scientific basis for horizontal comparison and classification of monitoring results from different environmental backgrounds on the same absolute numerical scale.

[0163] Furthermore, in a preferred embodiment that can effectively compensate for dimensional mixing distortion, calculating a first distance from each retained pixel to the dry edge and calculating a second distance from each retained pixel to the wet edge includes the following steps:

[0164] Step 506: Substitute the vegetation index of each preserved pixel into the dry edge equation and the wet edge equation respectively to obtain the corresponding theoretical dry edge albedo and theoretical wet edge albedo.

[0165] In this preferred ranging branch, the system keeps the x-coordinate of each pixel constant. For a pixel whose x-coordinate is a given vegetation index value, its vegetation index value is substituted into the dry-side linear function expression with known slope and intercept to calculate the theoretical maximum albedo that the surface may exhibit under the predetermined vegetation cover. This maximum theoretical value is the theoretical dry-side albedo. Simultaneously, it is substituted into the wet-side linear function expression to calculate the theoretical minimum albedo that the surface may exhibit under the same vegetation cover. This minimum theoretical value is the theoretical wet-side albedo.

[0166] Step 507: Calculate the absolute value of the difference between the surface albedo of each preserved pixel and the corresponding theoretical dry edge albedo to obtain the first distance.

[0167] Extract the actual observed vertical coordinate value of the pixel, i.e., the surface albedo. Calculate the algebraic difference between this surface albedo and the theoretical dry-edge albedo derived above, and take the absolute value to obtain the linear spacing along the vertical axis. This linear spacing is the first distance. This calculation method ensures that the calculated distance is a pure difference in reflectance in terms of physical dimensions.

[0168] Step 508: Calculate the absolute value of the difference between the surface albedo of each preserved pixel and the corresponding theoretical wet edge albedo to obtain the second distance.

[0169] Similarly, by calculating the absolute value of the algebraic difference between the surface albedo and the corresponding theoretical wet edge albedo, the length of the longitudinal line segment representing the pixel's deviation from the wet potential limit is obtained; this length is the second distance. When using the longitudinal distance to construct the relative distance ratio, both the numerator and denominator have the same physical unit of reflectance, and dividing them eliminates the unit difference.

[0170] Based on this, the preferred formula for longitudinal algebraic calculation of spatial geometric relative distance is as follows:

[0171] RDDI i =(Albedo dry_i -Albedo i ) / (Albedo dry_i -Albedo wet_i );

[0172] Among them, RDDI i For the spatial geometric relative distance of each preserved cell, Albedo dry_i For the corresponding theoretical dry edge albedo, Albedo i For true surface albedo, Albedowet_i This represents the corresponding theoretical wet edge albedo.

[0173] The following numerical example illustrates the complete calculation process of the aforementioned spatial geometric relative distance. Assume a pixel has a vegetation index of 0.30 and a surface albedo of 0.25. Assume the dry edge equation obtained after the preceding fitting steps is Albedo... dry =-0.20×NDVI+0.45, the wet edge equation is Albedo wet =0.15×NDVI+0.10. Substituting the vegetation index of 0.30 for this pixel into the dry edge equation, we get the theoretical dry edge albedo as -0.20×0.30+0.45=0.39; substituting it into the wet edge equation, we get the theoretical wet edge albedo as -0.15×0.30+0.10=0.055. The first distance is |0.39-0.25|=0.14, and the second distance is |0.25-0.055|=0.195. The total dry-wet distance is 0.14+0.195=0.335. The spatial geometric relative distance, i.e., the relative desertification difference index, is 0.14 / 0.335≈0.418. Since 0.418 is within the equidistant sub-interval of 0.4 to 0.6, this pixel is classified as moderate desertification. The above calculations show that the surface albedo of this pixel is slightly biased towards the drought side but does not reach the level of severe drought, and the judgment result is in line with physical expectations.

[0174] Furthermore, in another alternative embodiment, to prevent design loopholes in the geometric formula, the present invention also provides a Euclidean vertical distance calculation scheme based on a pure geometric coordinate system. Calculating the first distance from each retained cell to the dry edge, and calculating the second distance from each retained cell to the wet edge, includes the following steps:

[0175] Step 509: Convert the dry edge equation and the wet edge equation into the standard straight line general form equations respectively.

[0176] The system performs algebraic identity transformations on linear boundary functions with slope intercepts. By rearranging terms to align the dependent variable and constant terms on the same side of the equation, a standard linear equation in the form of a general algebraic polynomial summing to zero is constructed. This equation structure is a prerequisite mathematical building block for performing geometric projection distance measurements from arbitrary points to lines in space.

[0177] Step 510: Extract the two-dimensional coordinate points formed by the vegetation index and the surface albedo of each preserved pixel.

[0178] The system no longer processes the two attributes of a pixel separately and in isolation. Instead, it instantiates a geometric coordinate entity, i.e. a two-dimensional coordinate point, representing the true distribution location of the pixel, by using the vegetation index as the horizontal component and the surface albedo as the vertical component.

[0179] Step 511: Calculate the first Euclidean perpendicular distance from each two-dimensional coordinate point to the general equation of the standard straight line representing the dry side, and use it as the first distance.

[0180] The system invokes the classic analytical geometry formula for the distance from a point to a line. Specifically, it substitutes the x and y coordinate components of the cell into the left-hand polynomial of the general equation of the dry side and takes the absolute value, then divides it by the square root of the sum of the square of the slope parameter and 1. The calculation formula is as follows:

[0181] d dry =|a dry *NDVI i -Albedo i +b dry | / sqrt(a dry 2 +1);

[0182] Where, d dry For the calculated first Euclidean perpendicular distance, a dry NDVI represents the slope parameter of the general form equation of the standard straight line representing the dry side. i For the vegetation index components of two-dimensional coordinate points, Albedo i Let b be the surface albedo component of a two-dimensional coordinate point. dry Let be the intercept constant of the standard straight line equation representing the dry side, sqrt be the square root function, and a be the intercept constant. dry 2 It is the square of the slope parameter.

[0183] The distance calculated above is directly used as the first distance of this pixel.

[0184] Step 512: Calculate the second Euclidean perpendicular distance from each two-dimensional coordinate point to the general equation of the standard straight line representing the wet edge, and use it as the second distance.

[0185] Using the same analytical geometric logic, the corresponding parameters of the general equation of the standard straight line representing the wet side are substituted into the above ranging model to calculate the geometrically shortest line segment length projected from the pixel coordinate point onto the wet side along the normal direction perpendicular to the wet side line. This calculation result is adopted as the second distance and subsequently fed back into the ratio model to participate in the construction of the relative distance. This pure geometric ranging branch scheme can output equivalent substitution results highly consistent with the longitudinal ranging method in certain predetermined climatic zones where the slopes of the dry and wet sides are similar.

[0186] As a downgraded alternative to the aforementioned dual-track relative ranging scheme, in certain unique local scenarios where cross-regional lateral comparisons are not required, this scheme can strip away the relative ranging logic and directly degenerate into the absolute difference index method. In this alternative, the system only needs to rely on the slope parameter of the dry-side equation to directly construct the difference calculation model:

[0187] DDI i =a dry *NDVI i -Albedo i ;

[0188] Among them, DDI i For the extracted absolute difference index, a dry NDVI is the slope parameter of the dry-side equation. i For the vegetation index of pixels, Albedo i denoted as the surface albedo of a pixel.

[0189] This alternative is computationally simple, but because it lacks the ability to normalize for environmental background tolerance, its absolute output value is prone to drastic fluctuations with small changes in image phase and boundary slope. It cannot be applied with a uniform classification threshold and is only retained as a low-level downgrade alternative when computing power is extremely limited.

[0190] Example 6 details how, after calculating the spatial geometric relative distance, a mandatory error correction fallback algorithm based on physical principles is introduced to prevent system crashes caused by underlying mathematical deadlocks such as division by zero in the denominator and numerical overflow during large-scale automated batch processing in the cloud. Through boundary checks and multi-segment mathematical truncation assignment, the calculated relative desertification difference index is strictly constrained within a physically interpretable normalized interval.

[0191] Step 601: The spatial geometric relative distances are processed using a pre-defined boundary-crossing physical constraint truncation mechanism to obtain the relative desertification difference index for each retained pixel. This includes the following steps:

[0192] In the automated feature space construction of remote sensing big data, the use of percentiles as candidate points for boundary fitting in the preliminary steps indicates that, according to objective statistical laws, some samples will inevitably be distributed outside the fitted line segment. Simultaneously, in nature, the wet-dry line in high-vegetation-coverage areas exhibits a very small mathematical convergence angle. This objective situation inevitably leads to some retained pixels having actual surface albedo exceeding the theoretical fitted envelope. If the system directly substitutes these outliers into the aforementioned spatial geometric relative distance formula, it will generate an abnormal overflow index greater than 1 or less than 0, and may even produce a fatal error with a denominator of 0 at the predetermined vegetation index where the wet-dry edges mathematically intersect, causing the entire cloud computing graph to collapse. This embodiment constructs a pre-emptive physical constraint truncation mechanism as a safety net to forcibly intercept and correct abnormal pixels that may cause underlying errors or invalid physical meanings, outputting a safe, legal, and normalized relative desertification difference index.

[0193] Step 602: Substitute the vegetation index of each preserved pixel into the dry edge equation and the wet edge equation respectively to calculate the corresponding theoretical dry edge albedo and theoretical wet edge albedo.

[0194] For each predetermined pixel being calculated, its unique abscissa component, the vegetation index, is extracted again. To establish a local benchmark, this value is independently substituted into the previously globally fitted dry-edge and wet-edge equations. Through independent derivation and calculation of two linear equations, the system obtains the theoretically achievable maximum and minimum albedo limits for the land surface under the predetermined vegetation cover. These two local limit values ​​are the theoretical dry-edge albedo and the theoretical wet-edge albedo, respectively. Through independent recalculation, the system establishes the basic geometric reference for subsequent triple cross-comparison and boundary-crossing detection.

[0195] Step 603: If the corresponding theoretical dry edge albedo is less than or equal to the theoretical wet edge albedo, then the cross-failure masking mechanism is triggered to remove the corresponding retained pixels.

[0196] After extracting the local upper and lower boundary limits, a mathematical validity pre-interception based on common sense is performed. Physically, the surface albedo in arid conditions must always be higher than that in wet conditions with the same vegetation cover. However, in the abscissa edge areas where vegetation is extremely dense or sparse, due to the natural convergence of sample points, the dry edge line fitted by linear regression and the wet edge line may intersect or even be inverted and cross each other in the mathematical coordinate system.

[0197] Specifically, the system compares the calculated theoretical dry edge albedo with the theoretical wet edge albedo. If the theoretical dry edge albedo is found to be numerically less than the theoretical wet edge albedo, it indicates that the denominator in subsequent calculations will be zero or negative. For example, in areas with higher vegetation cover, the system calculates a theoretical dry edge albedo of 0.10, while the theoretical wet edge albedo is 0.12. The system detects this anomaly and immediately triggers a cross-failure masking mechanism. The value of the corresponding reserved pixel is forcibly set to a transparent state with no data. By executing this mechanism, anomalous pixels that are in a mathematical deadlock state and have no physical significance for drought monitoring are removed from the computational flow queue, preventing the underlying operator from crashing.

[0198] Step 604: If the surface albedo of each preserved pixel is greater than or equal to the corresponding theoretical dry edge albedo, the extreme drought cutoff mechanism is triggered, and the relative desertification difference index is assigned a value of 0.

[0199] For normal region pixels that have passed the above intersection test, the system compares their actual observed attributes with the theoretical upper limit. The true surface albedo of the pixel is extracted, and it is determined whether it has exceeded the driest ceiling defined by the fitted envelope. If the actual surface albedo of a pixel is numerically not less than the theoretical dry edge albedo under its corresponding vegetation conditions, it indicates that the actual soil moisture deficiency of that pixel is even more severe than the statistically extreme dry edge point.

[0200] For example, the actual surface albedo of a pixel might be 0.50, while the theoretical dry-edge albedo calculated by the system is only 0.42. Without intervention, substituting this into the pre-ratio formula would produce a negative abnormal relative index. In this case, the system immediately intercepts the pixel and activates the extreme drought cutoff mechanism. The system bypasses the normal ratio division operation and directly uses the hard-coded constant assignment function to forcibly assign the pixel's final relative desertification difference index a value of 0. In relative distance calculations, 0 represents a zero deviation from the dry edge, i.e., the most severe drought level. This forced zeroing operation, while preserving the extreme drought physical properties, effectively suppresses algebraic overflow errors caused by quantile exclusion.

[0201] Step 605: If the surface albedo of each preserved pixel is less than or equal to the corresponding theoretical wet edge albedo, the healthy vegetation cutoff mechanism is triggered, and the relative desertification difference index is assigned a value of 1.

[0202] Similarly, the system performs symmetrical over-limit capture and interception for the lower limit of surface albedo. If the actual surface albedo of a pixel is abnormally low, numerically not greater than the theoretically derived wet edge albedo, it indicates that the environment in which the pixel is located is extremely humid or that there is residual water shadow on the surface. At this time, the system triggers the corresponding healthy vegetation cutoff mechanism. The ratio ranging calculation is forcibly stopped, and the final relative desertification difference index of the predetermined pixel is directly assigned a constant value of 1. In relative distance, 1 represents coincidence with the wet edge, i.e., the optimal water supply state.

[0203] Step 606: For each preserved pixel that has not triggered the cross-failure masking mechanism, the extreme drought cutoff mechanism, or the healthy vegetation cutoff mechanism, the corresponding spatial geometric relative distance is directly used as the relative desertification difference index.

[0204] Only normal pixels whose actual surface albedo strictly falls within the envelope of the theoretical dry-edge albedo and theoretical wet-edge albedo can be exempted from the aforementioned triple extreme capture mechanism. For compliant pixel groups, such as normal pixels with an actual albedo of 0.20 and no boundary inversion, the system shifts and formally assigns the reasonable relative distance value within the open interval of 0 to 1 obtained from the previous steps through conventional ratio division to the relative desertification difference index.

[0205] To accomplish the aforementioned complex conditional branching and parallel decision interception using a single system-level codebase, this embodiment provides an algorithmic structure that transforms textual logic into a purely multi-segment mathematical piecewise function. This multi-segment mathematical truncation assignment function is expressed as follows:

[0206] RDDI i =Piecewise(masked,Albedo dry_i <=Albedo wet_i ;0,Albedo i >=Albedo dry_i ;1,Albedo i <=Albedo wet_i ;D rel_i (otherwise);

[0207] Among them, RDDI i The final output is the relative desertification difference index. Piecewise is a multi-stage conditional branching function operator. Albedo i To preserve the true surface albedo of each pixel, Albedo dry_i For the corresponding theoretical dry edge albedo, Albedo wet_i D represents the corresponding theoretical wet edge albedo, masked is the data-free variable representing the triggering cross-failure masking mechanism, and 0 is the constant representing the most severe drought assignment that triggers the extreme drought cutoff mechanism. rel_i "Otherwise" represents the spatial geometric relative distance of normal pixels that have not exceeded the boundary.

[0208] Example 7 elaborates on how, after obtaining the relative desertification difference index that has undergone physical boundary truncation and strict normalization to a predetermined numerical range, a universally applicable equidistant mapping standard across time and space is constructed. This eliminates the drawbacks of relying on manual experience and blindly trying to fix classification thresholds in traditional remote sensing monitoring, and enables effective automatic interpretation and visualization of desertification levels in large-scale regions.

[0209] Step 701: Based on the relative desertification difference index of each preserved pixel, the desertification level of the target area is classified to obtain the desertification monitoring results, specifically:

[0210] In the final engineering implementation stage of remote sensing monitoring big data, accurately converting continuous floating-point index values ​​into discrete classification labels with clear ecological guidance is crucial for business decision-making. Due to the dual processing in Examples 5 and 6, the relative desertification difference index output by this invention is constrained by physical mechanisms and naturally projected into a closed interval of 0 to 1. This allows the index to not only characterize the succession position relative to extremely arid and extremely wet states, but also, by eliminating underlying environmental differences and dimensional interference, achieve a generalization ability that maintains its relative physical meaning across different ecologically vulnerable areas and different climatic years. This embodiment constructs an automatic classification and engineering visualization mapping method that is independent of subjective human labeling, transforming complex ecological response gradients into intuitive desertification monitoring results.

[0211] Step 702: Obtain the pre-configured equidistant division step size.

[0212] The system retrieves pre-defined scalar parameters from memory to segment the continuous exponential space. In general operational scenarios for desertification monitoring, since the soil moisture content and vegetation cover, represented by the relative desertification difference index, exhibit an approximately linear positive correlation within the dry-wet boundary, and given the lack of a universal monitoring framework with numerous ground-based control points for nonlinear function calibration, using equidistant linear segmentation is statistically the most unbiased and robust prior strategy. The pre-configured equidistant segmentation step size determines the granularity of the severity levels in the final output monitoring results. In actual engineering configurations, the pre-configured equidistant segmentation step size can be uniformly set to 0.2.

[0213] Step 703: Based on the equidistant division step size, the normalized interval with values ​​from 0 to 1 is divided into multiple equidistant sub-intervals, and a desertification level label is assigned to each equidistant sub-interval.

[0214] The system utilizes the obtained equidistant division step size to perform an equal division operation on the normalized interval, which is strictly locked by the physical extrema to a value between 0 and 1. Specifically, taking an equidistant division step size of 0.2 as an example, the continuous normalized interval is uniformly divided, generating five non-overlapping, progressively equidistant sub-intervals in the mathematical dimension. The five equidistant sub-intervals are arranged from the arid limit to the wet limit as follows: [0,0.2), [0.2,0.4), [0.4,0.6), [0.6,0.8), and [0.8,1].

[0215] After completing the equidistant partitioning of the mathematical space, five purely algebraic intervals were bound to specific ecological business semantics. Following the decreasing trend of water stress, each equidistant sub-interval was assigned a corresponding desertification level label. Specifically, the first equidistant sub-interval was assigned the label of extreme drought or severe desertification, indicating that the pixel's surface albedo was extremely close to the drought envelope and soil moisture was extremely scarce; the second equidistant sub-interval was assigned the label of severe drought or severe desertification, indicating that soil moisture was significantly below the regional median level; the third equidistant sub-interval was assigned the label of moderate drought or moderate desertification, representing that the pixel was in an intermediate transitional state between the dry and wet envelopes; the fourth equidistant sub-interval was assigned the label of mild drought or mild desertification, indicating that moisture was slightly below the adequacy line; and the fifth equidistant sub-interval was assigned the label of no drought or no desertification, indicating that the area had abundant moisture and healthy vegetation. This pre-defined label mapping system, thanks to the anti-interference capability of relative distance calculation, can maintain good consistency in continuous automated batch processing across a wide spatial scale, avoiding the heavy workload of manually calibrating the grading threshold for each area.

[0216] Step 704: Determine the equidistant sub-interval to which the relative desertification difference index of each preserved pixel belongs, and extract the corresponding desertification level label to obtain the desertification monitoring results.

[0217] The system iterates through each retained pixel participating in the calculation and obtains its relative desertification difference index value after physical truncation correction. This continuous floating-point value is compared with the upper and lower thresholds of the previously defined equidistant sub-intervals to determine its unique numerical segment. Once it is confirmed that the value falls within a predetermined equidistant sub-interval, the discrete desertification level label bound to that sub-interval is immediately extracted and used as the final desertification interpretation conclusion for that pixel, thus obtaining the desertification monitoring result.

[0218] In the final visualization stage of cloud-based engineering export, desertification monitoring results are specifically represented as single-band raster image layers with a clearly defined spatial reference coordinate system. For example, the system calls the image export interface, forcibly setting the output results to use the globally universal EPSG:4326 (the spatial reference identifier code corresponding to the WGS84 geographic coordinate system) spatial reference coordinate system, and strictly limiting the spatial resampling grid resolution to 30 meters. The system configures different pseudo-color rendering mapping relationships for five different discrete level labels according to business requirements. For example, continuous visualization color bands of red, orange, yellow, light green, and dark green can be configured sequentially from extremely severe to non-desertification. Simultaneously, when pushing large-scale remote sensing monitoring results to external cloud disks, the system sets the maximum output pixel limit parameter to 10. 13 This level ensures the smooth transfer and distribution of high-definition spatial distribution maps of desertification in large-scale areas.

[0219] In some alternative implementations, when the system performs local single-scene monitoring without cross-regional comparison using the alternative absolute difference index method in Example 5, an alternative set of empirical grading standards is provided to replace the above-mentioned normalized equidistant division method. When applying the absolute difference index, since its value is not distributed within the range of 0 to 1, but rather floats freely with the slope of the dry edge, a fixed empirical classification threshold table must be used for matching. In this alternative implementation, the system determines whether the absolute difference index of each pixel is less than a predetermined empirical constant threshold and assigns a corresponding label. Specifically, if the absolute difference index is less than -0.3, a label of extremely severe desertification is assigned; if the absolute difference index is greater than or equal to -0.3 and less than -0.1, a label of severe desertification is assigned; if the absolute difference index is greater than or equal to -0.1 and less than 0.1, a label of moderate desertification is assigned; if the absolute difference index is greater than or equal to 0.1 and less than 0.3, a label of mild desertification is assigned; and if the absolute difference index is greater than or equal to 0.3, a label of non-desertification is assigned. This empirical threshold is heavily dependent on the temporal and regional background parameters of the predetermined image, and has weak generalization ability. It is only used as a downgrade alternative when the computing power of the cloud computing platform is extremely scarce and cannot support bilateral relative constraint inference.

[0220] Example 8 provides a supplementary explanation of the aforementioned online cloud-based large-scale automated monitoring process. It elaborates on how, after the system executes the cloud-based black-box dual-boundary automatic fitting engine, the key intermediate process data structures are exported to the offline local computing environment. Then, using visualized scatter plot rendering and rigorous statistical accuracy calculations, the high scientific validity and reliability of the output parameters of the cloud-based adaptive noise-resistant bilateral robust fitting engine are demonstrated from a post-evaluation perspective.

[0221] Step 801: Export the extracted effective fitting sample set and the parameters of the dry and wet edge equations in the two-dimensional feature space to the offline local computing environment.

[0222] In the actual operational deployment of remote sensing monitoring big data, although this invention proposes a black-box dual-boundary parameter calculation engine that requires no manual intervention on the cloud computing platform server, in some high-risk or ecological red line assessment areas with high requirements for model absolute accuracy, it is necessary to perform white-box accuracy review of the dry-edge equations and wet-edge equations output from the cloud. This embodiment reclaims the data export mechanism and local offline verification process as the basis for supporting the verifiability of the system of this invention. Specifically, after extracting the effective fitting sample set containing pixel vegetation index and surface albedo from the cloud, the system provides a data downlink interface. The export task is called to serialize the spectral attribute pairs of each retained pixel into a structured comma-separated value data file CSV. This offline file is pushed to an external storage medium. At the same time, key parameters such as the slope and intercept of the dry-edge equation and the slope and intercept of the wet-edge equation calculated by the cloud fitting engine are also exported along with the data stream, preparing complete prerequisites for offline reproduction of the local environment.

[0223] Step 802: Obtain the exported valid fitted sample set in the offline local computing environment, and draw a two-dimensional feature space scatter plot and a one-to-one reference line.

[0224] The system utilizes a data processing framework supporting scientific visualization within the local computing environment for data parsing and reconstruction. It acquires a local comma-separated data file, extracts the vegetation index sequence as the x-axis array, and extracts the land surface albedo sequence as the y-axis observation array. An offline plotting and rendering library is then used to generate a high-resolution scatter plot. To prevent severe visual overlap caused by a large number of sample points on the map, the rendering size of individual scatter points is reduced, and a transparency attenuation coefficient is forcibly applied. Transparency overlay rendering technology makes the densely packed ecological core areas appear as deep, brightly colored centers, reproducing the natural distribution of the effectively fitted sample set. Simultaneously, to intuitively assess the deviation between actual observations and theoretical expectations, an additional diagonal line with strictly equal proportions on the x and y axes, passing through the origin, is overlaid on the scatter plot as a one-to-one reference line for evaluating data convergence characteristics.

[0225] Step 803: Overlay and render the theoretical dry edge line and the theoretical wet edge line in the two-dimensional feature space scatter plot.

[0226] To verify the accuracy of the cloud-based calculation parameters, the synchronously exported dry and wet boundary constant parameters were extracted. Using the slope and intercept of the dry edge equation, a series of continuous theoretical limiting surface albedo coordinate points were calculated within a given vegetation index abscissa range, and these points were connected to form a smooth straight line running through the upper edge of the scatter plot. Similarly, using the slope and intercept of the wet edge equation, a corresponding straight line was generated at the lower edge of the scatter plot. By forcibly overlaying the two physical boundary lines onto the previously rendered real scatter plot, it is possible to visually determine whether the straight line fitted using percentile noise-resistant extraction and dynamic weighting of sample size conforms to the confidence boundary of the vast majority of sample distributions.

[0227] Step 804: Calculate various regression statistical evaluation indicators using a comparative evaluation model to quantify the approximation degree of the cloud-based fitting parameters.

[0228] In addition to the geometrically intuitive verification described above, this embodiment introduces a model accuracy evaluation system based on statistical indicators in an offline local environment. The system substitutes the x-coordinate values ​​of all sample points in the offline data matrix into the derived dry-edge and wet-edge equations to calculate a series of theoretical surface albedo predictions derived from the model. The actual surface albedo observation values ​​carried by these sample points are extracted and used as the prediction input array and actual input array for evaluating the model, respectively.

[0229] Calculate the coefficient of determination between the predicted input array and the actual input array. The coefficient of determination quantifies the proportion of the variance of the actual sample distribution explained by the linear boundary fitted from the cloud. The closer the value is to 1, the more accurately the boundary fitted from the cloud captures the scattering pattern under extreme water stress. The calculation formula is as follows:

[0230] R square =1-(SS res / SS tot );

[0231] Among them, R square As the coefficient of determination, SS res SS is the sum of squared residuals between the actual observed value series and the theoretically predicted target value series. tot It is the sum of squares of the deviations between the true observation sequence and its mean.

[0232] The system calculates the root mean square error (RMSE) between the predicted input array and the actual input array. The RMSE is extremely sensitive to extreme outliers where the sample deviates from the straight line. By calculating the square root of the mean square of the differences between the observed and predicted values, the system can effectively verify whether the percentile noise reduction mechanism used in the cloud truly eliminates residual outliers. The calculation formula is as follows:

[0233] RMSE=sqrt(∑((Y pred_i-Y actual_i ) 2 ) / N sample );

[0234] Where RMSE is the root mean square error, sqrt is the square root function, ∑ is the summation operator, and Y... pred_i Y is the theoretical predicted surface albedo value for the i-th sample. actual_i Let N be the true surface albedo observation value for the i-th sample. sample This represents the total number of samples participating in the evaluation.

[0235] Furthermore, the mean absolute error (MAE) between the predicted input array and the actual input array is calculated. The MAE provides a robust error estimate unaffected by the squared amplification of extreme outliers, directly reflecting the average absolute distance between the model's predicted values ​​and the actual observed values. The calculation formula is as follows:

[0236] MAE=∑(|Y pred_i -Y actual_i |) / N sample ;

[0237] Where MAE is the mean absolute error, Y pred_i Y is the theoretical predicted surface albedo value for the i-th sample. actual_i Let N be the true surface albedo observation value for the i-th sample. sample This represents the total number of samples participating in the evaluation.

[0238] For example, in a test in a typical arid and semi-arid grassland transition zone, the dry-edge equation output by the cloud-based adaptive noise-resistant fitting engine achieved a coefficient of determination of 0.91, a root mean square error of 0.018, and a mean absolute error of 0.014; the wet-edge equation achieved a coefficient of determination of 0.87, a root mean square error of 0.012, and a mean absolute error of 0.009. These indicators show that the dual-boundary equation fitted using percentile noise-resistant extraction and dynamic weighting of in-bin sample numbers can accurately approximate the true physical boundary of the effective fitted sample set.

[0239] By calculating the three core accuracy metrics in an offline local environment, the system can use rigorous mathematical statistics to verify the high reliability of the cloud-based weighted regression engine with a large sample size. This auxiliary evaluation mechanism provides an indispensable quality inspection loop for the algorithm development and verification phases and for important decision-making applications.

[0240] The experimental results of the above embodiments demonstrate that the proposed dual-boundary fitting method based on percentile noise-resistant extraction and dynamic weighting of in-bin sample numbers significantly improves boundary fitting robustness compared to traditional absolute extremum extraction and ordinary least squares fitting methods. Furthermore, the relative desertification difference index constructed based on bilateral relative distances shows significantly improved cross-regional generalization ability compared to the traditional one-sided absolute difference index. Those skilled in the art will understand that the specific improvement may vary depending on the landform type of the study area, the image time phase, and the data quality.

[0241] According to one aspect of this application, a method for monitoring regional desertification based on remote sensing cloud computing may further include:

[0242] Step 901: Define the study area and prepare data.

[0243] Import the vector boundary of the study area into GEE, such as Shapefile, GeoJSON, etc.

[0244] Load and display the study area boundary, import by file name, and select the applicable Landsat image dataset, such as Landsat 8 SR Tier 1;

[0245] Set a time frame, such as the vegetation growing season: June to September, and apply cloud masking and radiometric calibration pretreatment.

[0246] Step 902: Calculate the vegetation index NDVI and the surface albedo. The specific calculation formula has been explained in the previous embodiments and will not be described in detail here.

[0247] Step 903: Construct the Albedo-NDVI feature space.

[0248] Several sampling points are randomly generated within the study area, such as 1000, and the NDVI and Albedo values ​​of each point are extracted;

[0249] ; ;

[0250] Where N and A are the normalized NDVI and Albedo, respectively; NDVI min Represents the minimum value of NDVI, Albedo max and Albedo min This represents the maximum and minimum values ​​of Albedo.

[0251] Export as a CSV file for subsequent scatter plotting and dry edge fitting.

[0252] Step 904: Dry edge fitting and parameter a determination.

[0253] Read the sampled data in the Python environment and plot the Albedo-NDVI scatter plot;

[0254] Albedo = a × NDVI + b;

[0255] Where 'a' is a coefficient and 'b' is a parameter. As desertification intensifies, the NDVI value gradually decreases, while the Albedo value gradually increases. The Albedo-NDVI two-dimensional feature space can intuitively reflect the desertification process and its surface features, such as... Figure 4 As shown.

[0256] Use Python to obtain a scatter plot, such as Figure 5 As shown in the figure, Original represents the initial 1000 randomly selected points; Trend represents the trend line; R 2 The coefficient of determination, or coefficient of determination, is a core indicator in statistics used to measure the goodness of fit of a regression model to observed data. Its value ranges from 0 to 1, and can also be expressed as 0% to 100%. The closer the value is to 1, the better the model fits the data, and the stronger the explanatory power of the independent variable on the variation of the dependent variable.

[0257] The dry-side equation is fitted using linear regression or binning statistics to obtain the slope parameter, which serves as a key parameter for DDI calculation.

[0258] Step 905: Calculate the Desertification Differential Index (DDI).

[0259] In GEE, the DDI value for each cell is calculated using the following formula:

[0260] ;

[0261] in, The slope of the fitted line between Albedo and NDVI is given. The larger the DDI, the less desertification; conversely, the lower the DDI, the more severe the desertification. A higher DDI value indicates a less severe degree of desertification, while a lower DDI value indicates a more severe degree. Based on the natural breakpoint of the DDI value or a preset threshold, pixels are divided into several desertification levels.

[0262] Step 906: Visualize and output the results.

[0263] Visualize desertification level maps in GEE; download them to the cloud drive;

[0264] Export DDI raster or classification results to Google Drive (a commercial cloud storage service) or a GIS platform for further analysis.

[0265] This application employs binning to extract percentiles instead of absolute extrema for extracting boundary candidate points, and performs weighted least squares regression with the number of samples within each bin as the weight. This effectively shields heterogeneous noise interference from residual thin clouds or saline-alkali land, suppresses the erroneous leverage effect of sparse edge samples, and improves the physical robustness of the double-boundary solution. It also solves the problem of fragile boundary fitting.

[0266] This application constructs a ratio framework based on the relative geometric distance between pixels and the dry / wet dual-boundary space. By normalizing the wet edge, it effectively offsets the dimensional and absolute numerical distortions caused by differences in environmental background in different climate zones, naturally projecting the desertification index into a closed interval of 0 to 1. This achieves cross-temporal and spatial universality and direct reuse of equidistant fixed classification thresholds, solving the problem of the index's lack of cross-domain generalization ability.

[0267] This application introduces a dual safeguard at both the spatial and algebraic levels. On one hand, it mitigates feature space collapse in highly homogeneous terrains through spatial homogeneity range verification and buffered extrapolation sampling. On the other hand, it proposes a multi-segmented physical constraint truncation mechanism, performing mask removal on inverted areas where dry and wet areas intersect and forcibly assigning values ​​to overflowing pixels. This approach effectively avoids system crashes caused by denominator division by zero and algebraic overflow in automated batch processing, thus resolving the problem of algebraic deadlock in low-level computation.

[0268] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for automatic monitoring of regional desertification based on remote sensing cloud computing, characterized in that, include: Based on the acquired time-series multispectral remote sensing images of the target area, the vegetation index and surface albedo are extracted pixel by pixel and combined into a global sampling set; Non-target ground object pixels in time-series multispectral remote sensing images are removed based on preset composite ground object identification rules to obtain a retained pixel sampling set; Based on a preset spatial homogeneity verification mechanism, spatial sampling is performed on the retained pixel sampling set to extract an effective fitting sample set; The effective fitted sample set is subjected to segmented statistical and regression processing to extract the dry edge equation and wet edge equation of the two-dimensional feature space composed of vegetation index and surface albedo. Calculate the spatial geometric relative distances of each retained pixel to the dry and wet edges; By combining the pre-set cross-boundary physical constraint truncation mechanism to process the spatial geometric relative distance, the relative desertification difference index of each retained pixel is obtained. Based on this, the desertification level of the target area is classified, and the desertification monitoring results are obtained. The effective fitted sample set is subjected to piecewise statistical and regression processing to extract the dry edge equation and wet edge equation of the two-dimensional feature space composed of vegetation index and surface albedo, including: Based on the preset partitioning step size, the vegetation indices in the effective fitted sample set are divided into intervals to obtain multiple vegetation index boxes. The number of pixels within each vegetation index box is counted to obtain the number of samples within the box; Boundary features are extracted for each vegetation index box to generate corresponding dry edge candidate points and wet edge candidate points; Using the number of samples within each vegetation index box as the regression weight, weighted regression processing was performed on dry edge candidate points and wet edge candidate points respectively to obtain dry edge equations and wet edge equations. By combining a pre-defined boundary-crossing physical constraint truncation mechanism with the spatial geometric relative distances, the relative desertification difference index of each preserved pixel is obtained, including: If the corresponding theoretical dry edge albedo is less than or equal to the theoretical wet edge albedo, the cross failure masking mechanism is triggered, and the corresponding retained pixels are removed. If the surface albedo of each preserved pixel is greater than or equal to the corresponding theoretical dry edge albedo, the extreme drought cutoff mechanism is triggered, and the relative desertification difference index is assigned a value of 0. If the surface albedo of each preserved pixel is less than or equal to the corresponding theoretical wet edge albedo, the healthy vegetation cutoff mechanism is triggered, and the relative desertification difference index is assigned a value of 1. For each preserved pixel that did not trigger the cross-failure masking mechanism, the extreme drought cutoff mechanism, or the healthy vegetation cutoff mechanism, the corresponding spatial geometric relative distance was directly used as the relative desertification difference index.

2. The method according to claim 1, characterized in that, Non-target ground feature pixels include water body pixels, ice and snow pixels, and built-up area pixels; Based on preset composite land cover identification rules, non-target land cover pixels are removed from time-series multispectral remote sensing images to obtain a retained pixel sample set, including: Reflectance of green band, shortwave infrared 1 band, and shortwave infrared 2 band was extracted from time-series multispectral remote sensing images. The normalized water index is calculated based on the green light band reflectance and the shortwave infrared 1 band reflectance, and water pixels are removed by combining the preset water body threshold. An improved normalized snow cover index is calculated based on the reflectivity of the green light band and the reflectivity of the short-wave infrared band 2, and snow and ice pixels are removed by combining the preset snow and ice threshold. Pixels in built-up areas are removed based on preset built-up area identification rules; The remaining pixels after the removal process are selected as the retained pixels, resulting in the retained pixel sampling set.

3. The method according to claim 2, characterized in that, Based on preset built-up area identification rules, built-up area pixels are removed, including: Obtain pre-configured global land cover classification products; Based on preset category labels, predefined category pixels are extracted from pre-configured global land cover classification products; Pixels in the time-series multispectral remote sensing image whose spatial location overlaps with that of pixels of a predetermined category are treated as built-up area pixels and removed.

4. The method according to claim 1, characterized in that, Based on a pre-defined spatial homogeneity verification mechanism, spatial sampling is performed on the retained pixel sampling set to extract a valid fitted sample set, including: Extract the maximum and minimum vegetation index values ​​of the retained pixels in the retained pixel sampling set; The difference between the maximum and minimum vegetation index values ​​is calculated to obtain the vegetation index range. Based on the comparison results between the vegetation index range and the preset homogeneity threshold, the corresponding spatial sampling strategy is executed to extract the effective fitted sample set.

5. The method according to claim 4, characterized in that, Based on the comparison results between the vegetation index range and the preset homogeneity threshold, the corresponding spatial sampling strategy is executed to extract the effective fitted sample set, including: If the vegetation index range is greater than or equal to the homogeneity threshold, then random sampling is performed on the retained pixel sampling set based on the pre-configured sampling size parameter to extract the effective fitted sample set; otherwise: Based on a preset buffer extension distance, the spatial boundary of the target area is extended outward to generate an extended buffer area; Based on time-series multispectral remote sensing images, the vegetation index and surface albedo of pixels located in the extended buffer area are extracted. Non-target feature pixels are eliminated based on composite land cover identification rules, and the vegetation index and surface albedo of the remaining pixels are used to form a candidate buffer pixel set. The candidate buffer cell set is merged with the reserved cell sampling set, and the merged cell set is randomly sampled based on the pre-configured sampling size parameter to extract the effective fitted sample set.

6. The method according to claim 1, characterized in that, For each vegetation index box, boundary features are extracted to generate corresponding dry edge candidate points and wet edge candidate points, including: Extract the 95th and 5th percentiles of the surface albedo of pixels within each vegetation index box; Extract the central value of vegetation index for each vegetation index box; The central value of the vegetation index is combined with the corresponding 95th percentile to form dry edge candidate points; The central value of the vegetation index is combined with the corresponding 5th percentile to form wet edge candidate points.

7. The method according to claim 1, characterized in that, Calculate the spatial geometric relative distances of each preserved pixel to the dry and wet edges, including: Calculate the first distance from each retained pixel to the dry edge; Calculate the second distance from each retained pixel to the wet edge; Calculate the sum of the first distance and the second distance to obtain the total dry and wet distance; Calculate the ratio of the first distance to the total dry and wet distance to obtain the spatial geometric relative distance.

8. The method according to claim 7, characterized in that, Calculate the first distance from each retained cell to the dry edge, and calculate the second distance from each retained cell to the wet edge, including: Substitute the vegetation index of each preserved pixel into the dry edge equation and the wet edge equation respectively to obtain the corresponding theoretical dry edge albedo and theoretical wet edge albedo. The absolute value of the difference between the surface albedo of each preserved pixel and the corresponding theoretical dry edge albedo is used to obtain the first distance; The second distance is obtained by calculating the absolute value of the difference between the surface albedo of each preserved pixel and the corresponding theoretical wet edge albedo.

Images