A method for dynamically evaluating ecological system function damage value of natural resources

Abstract

This invention relates to the field of ecological environment monitoring and assessment technology, and particularly to a method for dynamic assessment of the value of ecosystem function damage oriented towards natural resources. The method includes the following steps: collecting data on morning dew formation and dissipation in a target area, and identifying land cover to obtain a land cover classification map; extracting areas where morning dew exists from the land cover classification map and performing environmental correlation to obtain a dynamic fingerprint map of morning dew; dividing the target area and calculating the micro-regional evaporation rate based on the dynamic fingerprint map of morning dew; performing evaporation trend analysis based on the micro-regional evaporation rate to obtain evaporation trend characteristics; extracting an evaporation spatial feature map from the micro-regional raster map based on the micro-regional evaporation rate; and comparing healthy land based on the evaporation spatial feature map and evaporation trend characteristics to obtain a soil vitality response curve. This invention provides more scientific and timely data support for ecological restoration decision-making through ecological environment monitoring and assessment technology.

Description

Technical Field

[0001] This invention relates to the field of ecological environment monitoring and assessment technology, and in particular to a method for dynamic assessment of the value of ecosystem function damage oriented towards natural resources. Background Technology

[0002] Traditional ecosystem function damage assessments primarily rely on field sampling and laboratory analysis. This static, point-based data collection method is not only costly and time-consuming but also struggles to cover large areas, resulting in insufficient spatiotemporal representativeness of the assessment results and failing to meet the needs of large-scale ecological management. Furthermore, the lack of a clear quantitative correlation mechanism between ecological status indicators and ecological function value often relies on expert judgment or simple linear models, leading to significant subjectivity and uncertainty in the assessment process and making it difficult to accurately reflect the true value of ecological damage. Existing assessment techniques focus primarily on the current state of the ecosystem, lacking scientific predictions of the ecosystem's self-repair capacity and dynamic recovery processes. This fails to provide a scientific basis for ecological restoration decisions in the temporal dimension, resulting in inefficient allocation of restoration resources.

[0003] In summary, existing technologies suffer from problems such as highly invasive assessment methods, limited coverage, ambiguous correlation between ecological functions and value, and lack of dynamic restoration assessment, which urgently need to be addressed. Summary of the Invention

[0004] Therefore, it is necessary to provide a dynamic assessment method for the value of ecosystem function damage oriented towards natural resources in order to solve at least one of the above-mentioned technical problems.

[0005] To achieve the above objectives, a dynamic assessment method for the value of ecosystem function damage oriented towards natural resources includes the following steps:

[0006] Step S1: Collect data on the formation and dissipation of morning dew in the target area, identify land cover, and obtain a land cover classification map; extract the areas where morning dew exists from the land cover classification map, and perform environmental correlation to obtain a dynamic fingerprint map of morning dew.

[0007] Step S2: Divide the target area into 10m × 10m grid units to obtain a micro-region grid map; calculate the micro-region evaporation rate based on the micro-region grid map and the morning dew dynamic fingerprint spectrum; perform evaporation trend analysis based on the micro-region evaporation rate to obtain evaporation trend characteristics; extract the evaporation spatial feature map from the micro-region grid map based on the micro-region evaporation rate; compare the evaporation spatial feature map and evaporation trend characteristics with healthy land to obtain the soil vitality response curve;

[0008] Step S3: Acquire and spatially align regional geographic information and micro-region raster maps to obtain a comprehensive geographic information layer; identify key ecological functions based on the comprehensive geographic information layer and land cover classification map to obtain ecological function indicators; perform water regulation function correlation analysis on soil vitality response curves to obtain a water regulation function map; quantify ecological function relationships based on the water regulation function map and ecological function indicators to obtain an ecological function integrity matrix.

[0009] Step S4: Assess the damage value of the ecological function integrity matrix to obtain the value decay curve; conduct ecological restoration intervention based on the value decay curve to obtain the spatiotemporal map of ecological damage value.

[0010] This invention utilizes UAVs equipped with high-resolution thermal imaging and visible light cameras for time-series data acquisition, achieving non-invasive, low-cost, and high-frequency monitoring of large-area natural resource ecological conditions. This completely overcomes the drawbacks of traditional ground sampling methods, such as high destructiveness, high cost, long processing time, and limited coverage. Specifically, by combining simple linear iterative clustering, random forest, and conditional random field models, the accuracy and spatial continuity of land cover type identification are significantly improved. Crucially, this step innovatively employs a dual verification mechanism of thermal imaging cold spots and high visible light reflectivity to accurately identify morning dew areas, effectively eliminating interference from single features such as shadows and water bodies, ensuring the high reliability of morning dew distribution data. The resulting dynamic fingerprint map of morning dew quantifies the entire process of morning dew formation and dissipation (start, duration, and dissipation time) for each pixel, forming an unprecedented high-dimensional spatiotemporal dataset capable of finely characterizing the hydrothermal dynamics of the surface microenvironment. This provides a solid and information-rich data foundation for subsequent in-depth analysis. By transforming raw morning dew observation data into scientific indicators that profoundly reveal the intrinsic health status of the soil, this invention achieves a leap from physical phenomena to ecological condition diagnosis. By standardizing the evaporation rate using saturated vapor pressure difference (VPD), the interference of real-time environmental meteorological conditions on the evaporation process was successfully eliminated, allowing the extracted evaporation characteristics to more accurately reflect the physical and biological properties of the land itself. Using a fourth-order polynomial to fit the evaporation rate sequence, the complete dynamic characteristics of the evaporation process (peak value, duration, acceleration, etc.) can be objectively and quantitatively captured, providing a far more comprehensive picture than a single rate value. Geostatistical and spatial analysis methods, such as standard deviation ellipses and semi-variograms, were employed to reveal the spatial clustering, directional trends, and autocorrelation structure of evaporation patterns, providing a macroscopic understanding of the patterns that traditional point-based analyses cannot achieve. Finally, by comparing multi-dimensional temporal and spatial characteristics with a health reference baseline and integrating them into a soil vitality response curve that includes response intensity, velocity, and stability, a multi-dimensional and comprehensive soil health diagnostic system was established, providing diagnostic results that are more profound and comprehensive than single-index assessments. By introducing and integrating geographic information such as digital elevation models and soil databases, and employing mature scientific models such as the Rosetta model, fuzzy logic inference systems, and random forest regression, a solid mechanistic support and computational basis were provided for each functional correlation analysis, ensuring the scientific rigor and objectivity of the functional assessment. The introduction of ground quadrat measured data and the correction of model results greatly improved the accuracy and reliability of the assessment results. The final generated ecological function integrity matrix not only accurately presents the spatial differentiation of various functions with a high spatial resolution of 10 meters × 10 meters, but also provides a comprehensive index that can comprehensively reflect the overall health status of the regional ecosystem through weighted analytic hierarchy process (AHP), providing clear and quantifiable spatial targets for subsequent value assessment and management decisions.By introducing an assessment of ecosystem self-repair capacity and a logistic recovery model, this method overcomes the limitations of traditional static assessments. It can dynamically predict the natural decay process of damage value over time, providing a scientific basis in the time dimension for formulating long-term restoration plans and evaluating the benefits of natural restoration. Using the social discount rate to calculate the present value of future damage makes cross-period cost-benefit comparisons possible, making decision-making more economically rational. By simulating the cost-benefit ratios of different restoration interventions, it provides decision-makers with a quantitative tool to optimize resource allocation and select the most efficient restoration solutions. The resulting spatiotemporal map of ecological damage value is an interactive decision support system integrating the current spatial distribution of damage value, future evolution paths, and intervention benefit analysis. It can present complex assessment results in an intuitive and dynamic way, greatly improving the scientific rigor, timeliness, and accuracy of ecological management and restoration decisions.

[0011] Therefore, this invention proposes a dynamic assessment method for the value of ecosystem function damage oriented towards natural resources. By collecting time-series data on the formation and dissipation of morning dew using drones, a dynamic fingerprint map of morning dew is constructed, and a quantitative correlation model between morning dew evaporation characteristics and ecological functions is established. Combined with the assessment of ecosystem self-repair capacity, this method achieves non-invasive and low-cost monitoring of the ecological status of large-area land resources, forming a logical closed loop from physical observation data to ecological function and then to economic value. At the same time, it provides comprehensive assessment results including spatial distribution and temporal evolution, providing more scientific and timely data support for ecological restoration decisions. Attached Figure Description

[0012] Figure 1 A schematic diagram illustrating the steps of a dynamic assessment method for the value of ecosystem function damage oriented towards natural resources;

[0013] Figure 2 A schematic diagram of the ecological function integrity recovery curve;

[0014] Figure 3 This is a schematic diagram of the damage value decay curve and the discounted present value.

[0015] The objectives, features, and advantages of this invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

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

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

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

[0019] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides a method for dynamic assessment of the value of ecosystem function damage for natural resources, comprising the following steps:

[0020] Step S1: Collect data on the formation and dissipation of morning dew in the target area, identify land cover, and obtain a land cover classification map; extract the areas where morning dew exists from the land cover classification map, and perform environmental correlation to obtain a dynamic fingerprint map of morning dew.

[0021] In this embodiment of the invention, a hexacopter UAV equipped with a 640×512 thermal imaging camera and a 48-megapixel visible light camera conducts a baseline flight over a 1-square-kilometer target area two hours before sunrise. A visible light orthophoto with a ground resolution of 2 cm and a thermal orthophoto with a resolution of 15 cm are generated using a motion recovery structure algorithm, forming a baseline state map. Subsequently, from 45 minutes before sunrise to 120 minutes after sunrise, the same flight path is strictly followed 12 times at 15-minute intervals to collect a time-series raw image library. Using the set of visible light orthophotos with the best illumination, simple linear iterative clustering (SLIC) is used to segment it into superpixels. Spectral, texture, and geometric features are then extracted and input into a pre-trained random forest classifier for preliminary classification. Finally, a conditional random field (CRF) model is used for optimization to generate an accurate land cover classification map. Next, cold spots were identified by subtracting the thermal image from the baseline image at each time point and setting a threshold of -1.5 degrees Celsius. Simultaneously, the V channel of the visible light image was enhanced with CLAHE, and the top 5% of the brightness values ​​were extracted as high-reflectivity areas. A logical AND operation was performed on both images to accurately extract the morning dew distribution mask sequence. Finally, this mask sequence was correlated with synchronously recorded ambient temperature, humidity, and light intensity data. For each pixel, three indicators were calculated: the start time, duration, and dissipation time of the morning dew. These three indicators were stored as three bands in a multi-band raster dataset to generate the final dynamic fingerprint map of the morning dew.

[0022] Step S2: Divide the target area into 10m × 10m grid units to obtain a micro-region grid map; calculate the micro-region evaporation rate based on the micro-region grid map and the morning dew dynamic fingerprint spectrum; perform evaporation trend analysis based on the micro-region evaporation rate to obtain evaporation trend characteristics; extract the evaporation spatial feature map from the micro-region grid map based on the micro-region evaporation rate; compare the evaporation spatial feature map and evaporation trend characteristics with healthy land to obtain the soil vitality response curve;

[0023] In this embodiment of the invention, the target area is divided into 10m × 10m micro-region grids. Based on this grid, the morning dew distribution mask sequence is statistically analyzed by region to obtain a time series table of morning dew area for each micro-region, and the instantaneous evaporation rate is calculated using the central difference method. The saturated vapor pressure difference (VPD) is calculated using synchronously collected ambient temperature and humidity data. The instantaneous evaporation rate is standardized by dividing the VPD, and the peak value of the standardized rate is extracted as the characteristic evaporation rate of the micro-region. Next, the standardized evaporation rate time series of each micro-region is fitted with a fourth-order polynomial function, and four evaporation trend features—peak rate, peak occurrence time, effective evaporation duration, and maximum evaporation acceleration—are extracted from the fitted function. Simultaneously, the characteristic evaporation rate of the micro-region is assigned as an attribute to the grid center point to generate a spatial point set. By constructing a Queen adjacency space weight matrix, the Getis-OrdGi* statistic is calculated to generate a local hotspot distribution map. The standard deviation ellipse is calculated to obtain directional features, and a spherical model is fitted to the empirical semivariogram curve to extract three spatial dependence indicators: nugget, sill, and range, which together constitute the evaporation spatial feature map. Finally, based on the evaporation spatial feature map, the evaporation uniformity index and boundary clarity score are calculated. These two spatial indicators, along with four trend features, are compared with the preset healthy land reference database. The standardized deviation of the six indicators is calculated to form a health deviation matrix. These six deviations are then mapped and combined into quantitative scores for three dimensions: response intensity, response speed, and recovery stability. This three-dimensional feature vector constitutes the soil vitality response curve of the micro-region.

[0024] Step S3: Acquire and spatially align regional geographic information and micro-region raster maps to obtain a comprehensive geographic information layer; identify key ecological functions based on the comprehensive geographic information layer and land cover classification map to obtain ecological function indicators; perform water regulation function correlation analysis on soil vitality response curves to obtain a water regulation function map; quantify ecological function relationships based on the water regulation function map and ecological function indicators to obtain an ecological function integrity matrix.

[0025] In this embodiment of the invention, geographic information such as the digital elevation model (DEM), soil type, and vegetation cover of the study area is integrated and spatially aligned with a micro-regional raster map to form a comprehensive geographic information layer. Based on this, multiple ecological functions are quantitatively correlated on the soil vitality response curve. For the water regulation function, the "response speed" and "recovery stability" indices are extracted. Combined with topographic hydrological factors calculated from the DEM and saturated hydraulic conductivity and field capacity estimated from soil properties using the Rosetta model, infiltration, water holding, and water storage capacity index models are constructed respectively. Finally, a weighted summation is performed to obtain the water regulation function map. For the nutrient cycling function, a fuzzy logic inference system is used, with "response stability," "response intensity," and soil organic matter content as inputs, to output a nutrient cycling function index map. For the biodiversity support function, a weighted linear combination model is used to integrate the local Moran index, land use Shannon diversity index, and patch density to generate a biodiversity support function map. For the carbon sequestration function, a random forest regression model is used, with "evaporation duration," "response intensity," NDVI, and soil organic carbon as inputs, to predict net ecosystem productivity and generate a carbon sequestration function map. Subsequently, 30 quadrats were deployed in the study area for field measurements to verify and linearly correct the four functional maps, resulting in a corrected functional assessment set. Based on this, the percentage of impairment for each function was calculated by comparing it with the reference baseline value of a healthy ecosystem, forming a functional impairment table. Finally, the weights of the four functions were determined using the analytic hierarchy process (AHP), and a weighted comprehensive index of ecological function integrity for each microregion was calculated. Ultimately, a multidimensional ecological function integrity matrix was constructed, incorporating spatial location, the integrity of each individual function, and the comprehensive index.

[0026] Step S4: Assess the damage value of the ecological function integrity matrix to obtain the value decay curve; conduct ecological restoration intervention based on the value decay curve to obtain the spatiotemporal map of ecological damage value.

[0027] In this embodiment of the invention, a functional value reference table is constructed by combining literature research with local biomass and socioeconomic factors for correction. This table includes the annual economic value per unit area of ​​various functions under different land types. Based on this table and the ecological function integrity matrix, the annual economic loss caused by functional impairment in each micro-region is calculated and summarized into a basic damage value list. Next, a product model is used to integrate the "response stability" in the soil vitality response curve and the health status of the surrounding area to assess the natural restoration potential of each micro-region and generate a restoration potential distribution map. Then, the restoration potential is converted into the intrinsic growth rate r in the logistic recovery model to predict the time required for each function in each micro-region to recover to 95% level, forming a functional restoration time series table. Based on this recovery model, the function of damage value decay over time for each micro-region is derived, generating a set of value decay curves. A social discount rate of 3.5% is set to discount the damage value over the next 50 years into total present value, resulting in a discounted damage value table. Furthermore, an intervention measure library including costs and restoration acceleration coefficients is established to simulate the benefits of applying intervention measures to high-value damage areas, calculate the cost-benefit ratio, and form an intervention benefit evaluation table. Finally, in a geographic information system platform, the spatial distribution map of discounted damage value is used as the base map, and the value decay curve is integrated to achieve dynamic display over time. It also allows interactive querying of detailed attributes and intervention benefits of each micro-region, thereby generating a comprehensive spatiotemporal map of ecological damage value that integrates space, time, and information.

[0028] Preferably, step S1, which involves collecting data on the formation and dissipation of morning dew in the target area and identifying land cover, includes:

[0029] A baseline state map was obtained by collecting data on the target area before the formation of morning dew using a drone.

[0030] Based on the baseline state diagram, time-series data acquisition was performed during the morning dew period to obtain the original time-series image;

[0031] Land cover type identification is performed based on the original time-series images to obtain a land cover classification map.

[0032] In this embodiment of the invention, firstly, two hours before sunrise, a hexacopter industrial UAV equipped with a 640×512 pixel resolution thermal imaging camera and a 48-megapixel visible light camera is used to collect baseline data on a 1-square-kilometer forest-grassland ecotone with known geographical boundaries. The geographical boundary coordinates of the area are input into a ground control station, which automatically generates a serpentine scanning flight path covering the entire area, setting the flight altitude to 80 meters, the forward overlap rate to 75%, and the lateral overlap rate to 60%. The UAV flies autonomously according to this flight path, simultaneously acquiring raw thermal and visible light images covering the entire target area. After acquisition, all visible light images are input into an image stitching processing module. This module uses a structure from motion algorithm based on feature point matching and a multi-view stereo matching algorithm to generate an orthophoto image with a ground resolution of 2 centimeters. Similarly, all thermal images are stitched together to generate a thermal orthophoto image with a ground resolution of 15 centimeters. These two georegistered orthophoto images together constitute the baseline state map of the target area, providing an initial reference free from morning dew interference for subsequent temporal change analysis.

[0033] Based on the generated baseline state map and planned flight path, time-series data acquisition was conducted from 45 minutes before sunrise until the morning dew had largely dissipated 120 minutes after sunrise. The UAV strictly followed the same serpentine scanning path as during the baseline acquisition, repeating the flight at fixed 15-minute intervals. During each complete flight, the onboard thermal and visible light cameras recorded the entire process, acquiring image data covering the same spatial area as the baseline state map. This process involved a total of 12 flight sorties, each generating a set of hundreds of raw visible light and thermal images with precise timestamps and location information (POS data). These 12 sets of images acquired at different times together constitute a time-series raw image library recording the complete process of morning dew formation, development, and dissipation, with each image spatially corresponding precisely to the baseline state map.

[0034] Land cover type identification is performed using orthophoto maps generated by stitching together a set of visible light images from a time-series raw image database with optimal lighting conditions (e.g., images acquired 30 minutes after sunrise). The process begins by employing a Simple Linear Iterative Clustering (SLIC) algorithm to segment the orthophoto map into thousands of superpixel units with similar color and texture, using these superpixels as the basic analysis unit. Next, a multi-dimensional feature vector is extracted for each superpixel. This vector includes: spectral features composed of the mean and standard deviation of the R, G, and B channels; texture features such as contrast, correlation, energy, and homogeneity obtained by calculating the Gray-Level Co-occurrence Matrix (GLCM); and geometric features calculated from the superpixel area and perimeter. The extracted feature vectors are then input into a pre-trained random forest classifier model using known land cover samples (bare land, grassland, shrubs, woodland, and farmland). This model consists of multiple decision trees, which provide an initial land cover classification label for each superpixel through ensemble learning. To improve classification accuracy and ensure spatial continuity, the initial classification results are then input into a Conditional Random Field (CRF) model for optimization. The Conditional Random Field (CRF) model smooths the initial classification results by establishing a relational potential function between neighboring superpixel labels, correcting isolated misclassified points caused by spectral confusion. Ultimately, it generates a vectorized land cover classification map that accurately labels the boundaries and types of different land features, including bare land, grassland, shrubs, woodland, and farmland.

[0035] Preferably, step S1 involves extracting the areas where morning dew exists from the land cover classification map and performing environmental correlation, including:

[0036] Extracting temperature difference map sequences from land cover classification maps;

[0037] Temperature cold spots were identified from the temperature difference map sequence to obtain a cold spot marker map sequence;

[0038] The land cover classification map is subjected to reflective enhancement processing to obtain a reflective enhancement map;

[0039] Extracting high-reflectivity mask sequence from reflectivity enhancement images;

[0040] The morning dew area was determined based on the high reflectivity mask sequence and the cold spot marker sequence, and the morning dew distribution mask sequence was obtained.

[0041] Environmental factors were correlated with the distribution mask sequence of morning dew to obtain environmental-morning dew correlation data;

[0042] A dynamic fingerprint map of morning dew is generated based on the environmental-morning dew correlation data.

[0043] In this embodiment of the invention, thermal orthophoto image sequences and visible light orthophoto image sequences from a time-series raw image library are utilized, combined with a land cover classification map. First, to obtain a temperature difference map sequence, a pixel-by-pixel difference operation is performed between the thermal orthophoto image at each time point in the time-series raw image library and the thermal image in the baseline state map. That is, the temperature value at each time point is subtracted from the baseline temperature value, resulting in a series of temperature difference maps reflecting the change in surface temperature relative to the dew-free state. Next, for temperature cold spot identification, a fixed temperature threshold, such as -1.5 degrees Celsius, is applied to each image in the temperature difference map sequence, and all pixels with temperature changes below this threshold are marked as potential cold spot areas. Subsequently, a morphological opening operation (erosion followed by dilation) is performed on the marked binary image using a 3×3 pixel structuring element. This operation aims to remove isolated small-area cold spots caused by sensor noise, ultimately generating an accurate cold spot marking map sequence. Simultaneously, the visible light images are processed. To obtain the reflectivity enhancement images, the visible light orthophoto images at each time point in the time-series original image library are converted from the RGB color space to the HSV color space. A contrast-limited adaptive histogram equalization (CLAHE) algorithm is applied to the V (luminance) channel. This algorithm redistributes luminance by calculating the histogram of local image regions, thereby significantly enhancing the contrast between the specular reflection areas caused by morning dew and the surrounding background. Then, to extract the high reflectivity mask sequence, the pixel value histogram of the V channel of each reflectivity enhancement image is calculated, and pixels with luminance values ​​above the 95th percentile are identified as high reflectivity regions, generating corresponding binarized high reflectivity masks. After completing the above preparations, the morning dew area is determined by performing a pixel-level logical AND operation on the cold spot marker map and the high reflectivity mask map at the same time point. That is, a pixel is finally confirmed as a morning dew-covered area only when it is marked (value 1) in both images. This dual verification mechanism eliminates the interference of single features (such as low temperature caused by shadows or high reflectivity caused by wet surfaces), thereby generating a high-precision morning dew distribution mask sequence. Subsequently, environmental factor association processing is performed. The environmental temperature, humidity, and light intensity data synchronously recorded by the UAV in each time-series acquisition are precisely matched with the timestamp of each morning dew distribution mask. For environmental data during flight intervals, linear interpolation is used for estimation, forming an environment-morning dew association dataset containing fields such as time, mask, temperature, humidity, and light intensity.Finally, a dynamic fingerprint map of morning dew is generated based on the mask sequence of morning dew distribution and the environment-morning dew correlation data. This process involves: creating a multi-band raster data structure with the same spatial range and resolution as the target area; traversing the state of each pixel throughout the entire time series; and calculating and filling in three core dynamic indicators for that pixel: 1. The time point at which morning dew begins to appear; 2. The total duration of morning dew (obtained by multiplying the number of masks with a pixel value of 1 by a 15-minute time interval); 3. The time point at which morning dew completely disappears. These three indicators are stored as three independent bands in the raster data structure. The resulting multi-band data product, containing full-temporal information on the morning dew formation and disappearance process for each pixel, is the dynamic fingerprint map of morning dew.

[0044] Preferably, the calculation of the micro-region evaporation rate based on the micro-region grid map and the morning dew dynamic fingerprint spectrum in step S2 includes:

[0045] Extract the time series table of morning dew area from the dynamic fingerprint of morning dew based on the micro-region raster map;

[0046] Calculate the area change rate sequence based on the morning dew area time series table;

[0047] The instantaneous evaporation rate is calculated from the area change rate sequence to obtain the evaporation rate sequence;

[0048] Environmental factor sequences were extracted from the dynamic fingerprint of morning dew.

[0049] Temperature and humidity were standardized based on the environmental factor sequence and evaporation rate sequence to obtain the standardized evaporation rate.

[0050] The evaporation region correlation is performed between the standardized evaporation rate and the micro-region grid map to obtain the micro-region evaporation rate.

[0051] In this embodiment of the invention, based on the 10m × 10m cell boundary defined by the micro-region grid map, the morning dew distribution mask sequence is partitioned and statistically analyzed to extract the morning dew area time series table. Specifically, this process involves traversing each micro-region cell and, at each acquisition time point (e.g., ... , ,..., The total number of pixels marked as 1 by the morning dew distribution mask within the cell is calculated, and this total number of pixels is multiplied by the actual geographic area of ​​a single pixel (e.g., 0.04 square meters / pixel) to obtain the total morning dew coverage area of ​​the micro-region at that time. The morning dew area data of all micro-regions at all time points are summarized to form a morning dew area time series matrix with micro-region IDs as rows and time points as columns. Next, based on this time series matrix, the area change rate sequence is calculated. For the area time series A(t) of each micro-region, the central difference method is used to calculate its change rate at each time point, i.e., at time point... The rate of change of area is equal to (A( -1)-A( +1)) / (2×Δt), where Δt is the 15-minute data collection interval, thus obtaining a sequence characterizing the rate of decrease in morning dew area. Subsequently, the instantaneous evaporation rate is calculated from the area change rate sequence. The area change rate calculated in the previous process is directly regarded as the instantaneous evaporation rate of the micro-region, in square meters per minute, thereby generating a curve of evaporation rate versus time for each micro-region, i.e., an evaporation rate sequence. Simultaneously, ambient temperature (T) and relative humidity (RH) data synchronized with each collection time point are extracted from the environment-morning dew association dataset to form an environmental factor sequence. Temperature and humidity standardization is performed based on the environmental factor sequence and the evaporation rate sequence. The logic of this operation is to remove the influence of the background environment on the evaporation process to highlight the characteristics of the land itself. Specifically, the saturated vapor pressure difference (VPD) at each time point is calculated using temperature and humidity data. This VPD value characterizes the driving force for the atmosphere to absorb moisture from the land surface. Then, the instantaneous evaporation rate value of each micro-region at each time point is divided by the corresponding VPD value, i.e., standardized evaporation rate = instantaneous evaporation rate / VPD. This standardization process unifies evaporation rates under different environmental conditions to a comparable benchmark. Finally, the standardized evaporation rates and micro-region raster maps are correlated to determine the evaporation region. To assign a single, representative evaporation rate index to each micro-region, the peak value—the maximum standardized evaporation rate throughout the entire evaporation process—is extracted from the standardized evaporation rate time series of each micro-region. This peak value is used as the final evaporation rate feature value for that micro-region and assigned to the corresponding cell in the micro-region raster map. After completing the calculations and assignments for all micro-regions, a new raster map is generated. The value of each pixel in the map represents the characteristic evaporation rate of its corresponding 10m × 10m micro-region; this map is the micro-region evaporation rate map.

[0052] Preferably, step S2, which involves analyzing the evaporation trend based on the micro-region evaporation rate, includes:

[0053] The evaporation rate of the micro-region is serialized to obtain the evaporation rate sequence.

[0054] By fitting the trend of the evaporation rate sequence, the evaporation fitting function is obtained.

[0055] Extract peak rate features from the evaporation fitting function;

[0056] The evaporation duration is obtained by calculating the evaporation duration using the evaporation fitting function.

[0057] Evaporation acceleration is estimated by fitting the evaporation function to obtain the evaporation acceleration characteristics;

[0058] Based on the peak rate characteristics, evaporation duration, and evaporation acceleration characteristics, key time points are calibrated to obtain evaporation trend characteristics.

[0059] In this embodiment of the invention, for each 10m × 10m micro-region, the standardized evaporation rate values ​​at all collection time points are arranged in chronological order to form a discrete time-series data point set, which is the evaporation rate sequence. Next, to transform the discrete data points into a continuous, analyzable mathematical expression, a trend fitting operation is performed on the evaporation rate sequence of each micro-region. Specifically, the least squares method is used to fit the sequence to a fourth-order polynomial function: ,in express Evaporation rate at time, coefficient This is the unique set of parameters calculated by the algorithm that minimizes the sum of squared errors between the fitted curve and the original data points. The resulting continuous function is the evaporation fitting function for this micro-region. Subsequently, feature extraction is performed based on this function. To obtain the peak rate feature, the evaporation fitting function is... Find the first derivative and order Solving for the real roots of this equation reveals the moment t_peak when the rate reaches its extreme value. Substituting t_peak into the original function R(t) yields the peak evaporation rate R_max, and this (t_peak, R_max) coordinate pair represents the peak rate characteristic. Then, the evaporation duration is calculated, setting a rate threshold of 5% of the peak rate R_max. This is achieved by solving the equation... =0.05×R_max yields two time points, t_start and t_end, which represent the effective start and end of the evaporation process, respectively. The difference between them, (t_end-t_start), is the evaporation duration. Next, the evaporation acceleration is estimated by taking the second derivative of the evaporation fitting function R(t). The second derivative represents the acceleration of the change in evaporation rate. (Calculation) The maximum value within the time interval [t_start, t_peak] is the evaporation acceleration characteristic, which quantifies the rapid increase in the evaporation rate. Finally, based on the peak rate characteristic, evaporation duration, and evaporation acceleration characteristic, key time points are calibrated. The four core parameters calculated for each micro-region—peak rate R_max, peak occurrence time t_peak, effective evaporation duration, and maximum evaporation acceleration—are combined into a four-dimensional vector. This vector comprehensively and quantitatively describes the dynamic characteristics of the morning dew evaporation process in that micro-region. This vector set is the final generated evaporation trend characteristic set.

[0060] Preferably, step S2, extracting the evaporation spatial feature map from the micro-region grid map based on the micro-region evaporation rate, includes:

[0061] Spatial data processing is performed on the micro-region evaporation rate and the micro-region raster map to obtain the spatial point set of evaporation rate;

[0062] By constructing neighborhood relationships for the spatial point set of evaporation rates, a spatial weight matrix is ​​obtained;

[0063] Local hotspot analysis is performed based on the spatial point set of evaporation rate and the spatial weight matrix to obtain a local hotspot distribution map;

[0064] A directional characteristic table is obtained by performing a directional analysis on the spatial point set of evaporation rates.

[0065] Calculate the semivariogram curve of the spatial point set of evaporation rate;

[0066] Spatial dependence was assessed on the semivariogram curve to obtain a spatial dependence index;

[0067] An evaporation spatial feature map is generated based on spatial dependence indicators, local hotspot distribution maps, directional feature tables, and spatial dependence indicators.

[0068] In this embodiment of the invention, spatial data processing is performed on the micro-region evaporation rate and the micro-region grid map. The geometric center of each 10m × 10m grid cell is extracted as a point feature, and the micro-region evaporation rate value of that cell is used as the attribute of that point feature. This generates a spatial point set of evaporation rates containing all micro-region center points and their evaporation rate values. Next, neighborhood relationships are constructed for this spatial point set using the Queen adjacency criterion. That is, if two micro-region grid cells share an edge or a vertex, they are defined as neighbors. Based on this, an N×N spatial weight matrix W is generated, where N is the total number of micro-regions. If cells i and j are neighbors, then the matrix element Wij = 1; otherwise, Wij = 0. Then, local hotspot analysis is performed based on the spatial point set of evaporation rates and the spatial weight matrix. Specifically, the Getis-Ord Gi statistic is applied, and a Gi value is calculated for each point i. This value measures the deviation of the sum of the evaporation rate values ​​of point i and its neighborhood from the sum of all points. By comparing the calculated Gi* value with the expected value using Z-score and p-value tests, each micro-region is classified into statistically significant high-value clusters (hotspots), low-value clusters (cold spots), or insignificant clusters, ultimately generating a local hotspot distribution map. Simultaneously, directional analysis is performed on the spatial point set of evaporation rates. This is achieved by calculating the Standard Deviation Ellipse, an ellipse whose major axis indicates the primary distribution direction of high or low evaporation rate clusters, its minor axis indicates the secondary direction, and its area reflects the dispersion of the distribution. The center coordinates, rotation angle, semi-major axis length, and semi-minor axis length of the ellipse are calculated and recorded; these parameters together constitute a directional feature table. Furthermore, to quantify spatial autocorrelation, an empirical semivariogram curve of the spatial point set of evaporation rates was calculated. This process grouped all point pairs according to their distance h and calculated the semivariogram γ(h) within each distance group. Plotting γ(h) against h yielded the empirical semivariogram. A spherical model was then fitted onto this empirical graph. Three key spatial dependence indicators were extracted from the fitted theoretical model: the nugget value, representing random error and spatial variation smaller than the sampling scale; the sill value, representing the total variance of the data; and the range, representing the maximum distance at which spatial autocorrelation exists.Finally, an evaporation spatial feature map is generated based on spatial dependence indicators, local hotspot distribution maps, and directional feature tables. This map is a comprehensive geographic data product containing multiple data layers, specifically composed of the following layers: 1. The local hotspot distribution map is used as the base map; 2. The standard deviation ellipse from the directional feature table is overlaid as a vector layer; 3. The three spatial dependence indicators—nuclear value, sill value, and range—are stored as global attributes for the entire study area. This multi-layered dataset comprehensively reveals the spatial structure, clustering patterns, and directional trends of evaporation rates.

[0069] Of particular importance is the extraction of the evaporation spatial feature map from the micro-region raster map based on the micro-region evaporation rate, specifically as follows:

[0070] Spatial data processing is performed on the micro-region evaporation rate and the micro-region raster map to obtain the spatial point set of evaporation rate;

[0071] By constructing neighborhood relationships for the spatial point set of evaporation rates, a spatial weight matrix is ​​obtained;

[0072] Local hotspot analysis is performed based on the spatial point set of evaporation rate and the spatial weight matrix to obtain a local hotspot distribution map;

[0073] A directional characteristic table is obtained by performing a directional analysis on the spatial point set of evaporation rates.

[0074] Calculate the semivariogram curve of the spatial point set of evaporation rate;

[0075] Spatial dependence was assessed on the semivariogram curve to obtain a spatial dependence index;

[0076] An evaporation space feature map is generated based on the local hotspot distribution map of the spatial dependence index, the directional feature table, and the spatial dependence index.

[0077] In this embodiment of the invention, a geographic information processing script is used to traverse each 10m × 10m raster cell in the micro-region evaporation rate map. For each raster cell, its pixel value is first read, which represents the characteristic evaporation rate of the micro-region. Next, the geometric center coordinates (X, Y) of the raster cell are calculated. Then, a new point feature is created, and the calculated center coordinates (X, Y) are assigned to the geometric attribute of the point feature, while the read evaporation rate value is used as the numerical attribute of the point feature. This operation is repeated for all raster cells in the micro-region evaporation rate map until all cells have been transformed. Finally, all generated point features and their attribute data are stored in a new vector data file (e.g., Shapefile format), which is the spatial point set of evaporation rates. Each point in this point set uniquely represents a micro-region and carries its geographical location and evaporation rate information.

[0078] The Queen adjacency criterion, based on contiguity, is used to define neighborhoods. This criterion states that two micro-regions are neighbors if they share an edge or a vertex in the original raster graph. Based on this criterion, an N×N spatial weight matrix W is constructed, where N is the total number of micro-regions (i.e., the number of points in the spatial point set). The elements of this matrix... Assign values ​​according to the following rules: If the point... and points If a matrix element is determined to be a neighbor, then the matrix element... The value is 1; if they are not neighbors, then The value is 0. Also, by convention, a point itself is not considered its neighbor, therefore the elements on the matrix diagonal are... The sum of all elements in a row is always 0. For ease of calculation, this matrix is ​​typically row-normalized, meaning each element in a row is divided by the sum of the elements in that row (i.e., the total number of neighbors of that point), so that the sum of all elements in each row equals 1. This calculated and normalized spatial weight matrix W accurately and quantitatively describes the spatial adjacency relationships between all micro-regions within the entire study area.

[0079] The Getis-OrdGi* (pronounced Gi-star) statistic is used for calculation. This is performed for each point in the spatial set of evaporation rates. ,That The formula for calculating the statistic is: =[Σ( × )-X_bar×Σ ] / [S× ],in It is a point The evaporation rate value, In the spatial weight matrix and The weights between points are: X_bar is the average evaporation rate across all points, S is the standard deviation of the evaporation rate across all points, and n is the total number of points. The logic of this formula is to compare points... The difference between the local sum and the global sum of W and its neighborhood (defined by W). The calculated... The value is a Z-score, which represents the degree of deviation between the local observed value and the expected value. According to statistical principles, if a point has a significantly positive Z-score and a small p-value (e.g., p < 0.05), then the point and its neighborhood form a high-value cluster, i.e., a "hot spot"; if the Z-score is significantly negative and the p-value is small, then it forms a low-value cluster, i.e., a "cold spot"; if the Z-score is not significant, then it is a non-clustered area. After calculating and statistically testing all points, each point is assigned a classification label (hot spot, cold spot, insignificant), and then rendered on a map using different colors (e.g., red for hot spots, blue for cold spots), ultimately generating a local hot spot distribution map that visually displays the spatial clustering pattern of evaporation rates.

[0080] The Standard Deviation Ellipse (SDE) method is employed. SDE constructs an ellipse by calculating the standard deviation of all point coordinates. This ellipse summarizes the central tendency, distribution direction, and shape of the point set. The calculation process includes: 1. Calculating the mean center of all points (…). , 1. The center of the ellipse is determined by the following steps: 2. Calculate the standard deviations σx and σy of the x and y coordinates. 3. Calculate the covariance of the point coordinates to determine the rotation direction of the ellipse. The major axis of the ellipse indicates the primary direction of the point data distribution, while the minor axis indicates the secondary distribution direction. The area of ​​the ellipse reflects the concentration or dispersion of the point distribution; a smaller area indicates a more concentrated distribution. After calculation, extract and record the five core parameters of the ellipse: the X-coordinate of the center point, the Y-coordinate of the center point, the rotation angle (the angle between the major axis and true north), the length of the semi-major axis, and the length of the semi-minor axis. Organize these five parameters into a table, which is the directional characteristic table. It quantitatively describes the global geometric characteristics of the evaporation rate distribution throughout the entire study area on a macroscopic scale.

[0081] Geostatistical methods are used to reveal the spatial autocorrelation of evaporation rates, that is, how the similarity between the evaporation rate value of a point and its neighboring points varies with distance. The semivariogram γ(h) measures the similarity between two points h apart (…). , The expected value of the square of the difference between the attribute values ​​of +h is half of the expected value, and its formula is: γ(h)=(1 / 2N(h))×Σ[z( )-z( +h)] 2 , where z( ) is a point The evaporation rate value is given by N(h), which is the total number of point pairs at a distance of h. The specific steps are as follows: 1. Calculate the distance between all point pairs in the spatial set of evaporation rates and the square of the difference in evaporation rates. 2. Group these point pairs according to distance h, for example, with a distance step size of 5 meters (lag distance). 3. Calculate the mean semivariogram within each distance group, obtaining a series of (h, γ(h)) data points. 4. Plot these data points on a two-dimensional coordinate system, with the horizontal axis representing distance h and the vertical axis representing the semivariogram γ(h). These scatter plots constitute the empirical semivariogram curve. This curve visually shows how the difference in evaporation rate changes with increasing distance.

[0082] To extract quantitative structured information from the empirical semivariogram curve, a theoretical model is needed to fit it. A spherical model is chosen for fitting, with the functional expression: when h ≤ a, γ(h) = +C×(1.5×(h / a)-0.5×(h / a) When h > a, γ(h) = +C. The model parameters were adjusted using nonlinear least squares to achieve the highest possible fit between the theoretical curve and the empirical scatter plot. After successful fitting, three key spatial dependence indicators were extracted from the theoretical model: 1. Nugget value (Nugget, 1. The intercept of the model at h=0, representing the random component caused by measurement error and spatial variability smaller than the sampling scale. 2. Sill value (Sill, +C): The plateau value reached by the semivariogram as distance increases; it represents the total variance of the sample data. 3. Range (a): The distance at which the semivariogram reaches the sill value; it represents the maximum distance at which spatial autocorrelation exists, and points beyond this distance are statistically considered independent. These three indicators collectively quantify the spatial structural characteristics of the evaporation rate: the degree of randomness, the degree of total variability, and the range of spatial correlation.

[0083] The results of all preceding spatial analyses are integrated to form a multi-layered comprehensive geographic data product that comprehensively describes the spatial characteristics of evaporation rates. The specific construction method is as follows: 1. A local hotspot distribution map is used as the base map of this comprehensive map, providing information on the aggregation of evaporation rates at the microscale. 2. The standard deviation ellipse calculated from the directional feature table is used as an independent vector layer, overlaid on the hotspot distribution map, which displays the global distribution direction and dispersion trend at the macroscale. 3. Three spatial dependence indicators (nuclear value, sill value, and range) are stored and labeled as metadata or global attributes for the entire study area, providing key parameters regarding spatial autocorrelation structure. These three layers (hotspot map, ellipse map, and global parameters) are integrated into a single geographic information system project file, collectively forming an evaporation spatial feature map. This map, through the combination of data at different levels, achieves a comprehensive and three-dimensional display from micro-aggregation patterns to macro-directional trends and spatial structure parameters.

[0084] Preferably, step S2, which compares healthy land based on the evaporation spatial characteristic map and evaporation trend characteristics, includes:

[0085] Calculate the evaporation uniformity index based on the evaporation space characteristic map and the evaporation rate of the micro-region;

[0086] Based on the evaporation space feature map, the boundary clarity of the micro-region raster map is evaluated to obtain the boundary feature score map;

[0087] A comparative analysis of healthy land was conducted on the evaporation trend characteristics, evaporation uniformity index map, and boundary feature scoring map to obtain the health deviation matrix.

[0088] Soil vitality response curves are generated based on the health deviation matrix.

[0089] In this embodiment of the invention, to calculate the evaporation uniformity index map, a 5×5 moving window is used to traverse the micro-region evaporation rate map. For each micro-region in the map, the evaporation rate values ​​of itself and its 24 neighboring micro-regions are extracted. The standard deviation and arithmetic mean of these 25 values ​​are calculated. Then, the coefficient of variation of the central micro-region is obtained using the formula: Uniformity Index = (Standard Deviation / Arithmetic Mean) × 100%. This coefficient is the evaporation uniformity index; the lower the value, the more uniform the local evaporation. The calculated index is assigned to the corresponding central micro-region. After completing the calculation for all micro-regions, an evaporation uniformity index map covering the entire study area is generated. Next, to evaluate boundary sharpness, a 3×3 Sobel operator is applied to the micro-region evaporation rate map. This operator quantifies the sharpness of the boundary by calculating the gradient of the evaporation rate values ​​between each micro-region and its eight neighbors. Regions with larger gradient magnitudes indicate drastic changes in evaporation rate and clear boundaries. The calculated gradient magnitude map is used as the boundary feature scoring map. Subsequently, a comparative analysis of healthy land was conducted. This process incorporated a pre-established healthy land reference database, which defined a set of ideal characteristic parameter standard intervals for each land cover type (e.g., healthy grassland), including: peak rate, peak occurrence time, evaporation duration, evaporation acceleration, evaporation evenness index, and boundary feature score. For each micro-region within the study area, the calculated six characteristic values ​​were compared with the standard intervals for the corresponding land cover type in the database. The standardized deviation of each characteristic was calculated using the formula: Deviation = |Measured value - Reference interval median| / (Reference interval width / 2). This value standardized the deviation of all characteristics to a range of 0 to 1. The deviations of the six characteristics of all micro-regions were integrated into an N×6 dimensional health deviation matrix, where N is the total number of micro-regions, each row of the matrix represents a micro-region, and each column represents the deviation of a characteristic. Finally, a soil vitality response curve is generated based on the health deviation matrix. The six deviation indicators in the health deviation matrix are recombined and mapped to three core soil vitality dimensions: 1. Response intensity, directly determined by the deviation of peak rate; 2. Response speed, a weighted combination of the deviations of peak occurrence time and evaporation acceleration; 3. Recovery stability, comprehensively evaluated by the deviations of evaporation duration, evaporation evenness index, and boundary feature score. A comprehensive score for these three dimensions is calculated for each microregion. This set of data, containing the three quantitative indicators of response intensity, speed, and stability, constitutes the soil vitality response curve for that microregion. It is not a two-dimensional graph, but a multi-dimensional feature vector that quantitatively describes the soil's ecological response characteristics.

[0090] Preferably, step S3, which involves performing a water regulation function correlation analysis on the soil vitality response curve, includes:

[0091] Extracting moisture characteristic parameters from soil vitality response curves;

[0092] Topographic and hydrological factors are calculated based on the integrated geographic information layer;

[0093] Soil physical properties are integrated from the comprehensive geographic information layer to obtain soil moisture characteristics;

[0094] Based on soil moisture characteristics and moisture feature parameters, the morning dew-infiltration relationship was constructed to obtain an infiltration capacity index map.

[0095] Based on soil moisture characteristics and moisture feature parameters, the relationship between morning dew and water holding capacity was constructed to obtain a water holding capacity index map.

[0096] A water storage capacity index map is constructed based on water characteristic parameters and topographic hydrological factors;

[0097] Functional regulation analysis was performed on the infiltration capacity index, water holding capacity index, and water storage capacity index to obtain the water regulation function diagram.

[0098] In this embodiment of the invention, to extract water characteristic parameters, two indicators most directly related to water dynamics, "response speed" and "recovery stability," are separated from the multidimensional feature vector of the soil vitality response curve of each micro-region. Response speed reflects the residence time of water on the surface, while recovery stability indirectly characterizes the water-holding capacity of the soil structure. Simultaneously, based on the digital elevation model (DEM) in the integrated geographic information layer, hydrological analysis tools are used to calculate the slope (unit: degrees) and runoff accumulation (the number of grid cells flowing into the unit from upstream) for each 10m × 10m micro-region. These two indicators together constitute topographic hydrological factors. Further, the integrated geographic information layer integrates soil physical properties, extracting soil texture (percentage of sand, silt, and clay) and organic matter content data for each micro-region from the integrated soil database. This data is then input into a pedotransfer function model called Rosetta, which estimates the saturated hydraulic conductivity (Ksat) and field capacity (FieldCapacity) based on basic soil physical properties. These two estimates constitute the quantified soil water characteristics. Next, a morning dew-infiltration relationship was constructed based on soil moisture characteristics and parameters. The logic is that a slower evaporation response indicates more time for water to infiltrate. Combined with the soil's high saturated hydraulic conductivity, this signifies strong infiltration capacity. Therefore, the infiltration capacity index is defined as the product of the standardized (1 - response speed score) and the standardized saturated hydraulic conductivity, thus generating an infiltration capacity index map. Similarly, a morning dew-water holding capacity relationship was constructed. Higher recovery stability indicates a more effective soil structure in retaining moisture. Combined with high field capacity, this signifies strong water holding capacity. Therefore, the water holding capacity index is defined as the product of the standardized recovery stability score and the standardized field capacity, generating a water holding capacity index map. Subsequently, a water regulation capacity index map is constructed based on water characteristic parameters and topographic hydrological factors. The logic is that areas with gentle terrain (low slope), located along runoff paths (high runoff accumulation), and stable evaporation processes (high recovery stability) possess stronger surface runoff regulation capacity. This index is calculated using a product model integrating (1-standardized slope), standardized runoff accumulation, and standardized recovery stability scores. Finally, the infiltration capacity index map, water holding capacity index map, and water regulation capacity index map undergo functional adjustment analysis. By assigning preset weighting coefficients to these three index maps (e.g., setting W_infiltration = 0.4, W_water holding = 0.4, W_regulation = 0.2 based on regional ecological characteristics), the three index values ​​for each micro-region are weighted and summed to obtain a comprehensive water regulation function score. This comprehensive score is assigned to the corresponding cell in the micro-region raster map, ultimately generating a spatial distribution map that accurately quantifies the strength of water regulation capacity in each micro-region—this is the water regulation function map.

[0099] Preferably, step S3, which involves quantifying the ecological function relationship based on the water regulation function map and ecological function indicators, includes:

[0100] Nutrient cycling function correlation analysis was performed on the soil vitality response curve to obtain the nutrient cycling function map;

[0101] A biodiversity support function correlation analysis was performed on the soil vitality response curve to obtain a biodiversity support function map;

[0102] A carbon fixation function correlation analysis was performed on the soil activity response curve to obtain a carbon fixation function map;

[0103] Functional correlations of the water regulation function map, nutrient cycling function map, biodiversity support function map, and carbon fixation function map were validated to obtain a corrected functional assessment set.

[0104] Calculate the functional impairment scale based on the corrective function assessment set;

[0105] An ecological function integrity matrix is ​​generated based on the functional impairment table and ecological function indicators.

[0106] In this embodiment of the invention, a nutrient cycling potential assessment model based on soil vitality and soil biochemical properties is established. The core logic of this model is that soil microbial activity is a key driver of nutrient cycling, and microbial activity is directly affected by the stability and suitability of soil hydrothermal conditions. Two indicators, "response stability" and "response intensity," are extracted from the soil vitality response curve of each micro-region. A high "response stability" score indicates a smooth morning dew evaporation process, small fluctuations in the soil hydrothermal environment, and is conducive to microbial community stability; a moderate "response intensity" score indicates that the evaporation rate is neither too fast, leading to drought stress, nor too slow, leading to an anaerobic environment, placing it within the optimal range for microbial metabolism. Simultaneously, soil organic matter content (SOM) and carbon-to-nitrogen ratio (C:N) data for the corresponding micro-regions are extracted from the comprehensive geographic information layer. A three-input, single-output fuzzy logic inference system is constructed to quantify nutrient cycling function. The input variables of this system are: 1. Response stability score; 2. Response intensity score; 3. Soil organic matter content. The output variable is: nutrient cycling function index. Three fuzzy sets—low, medium, and high—are defined for each input variable, and their membership functions are set, such as using a Gaussian membership function. The system's rule base contains a series of "IF-THEN" rules, such as: "IF response stability is high AND response intensity is medium AND organic matter content is high THEN nutrient cycling function index is high"; "IF response stability is low AND response intensity is high THEN nutrient cycling function index is low". The three input variable values ​​for each micro-region are input into the system, and a nutrient cycling function index between 0 and 1 is calculated through fuzzification, rule inference, and defuzzification (e.g., using the centroid method). This index is assigned to each micro-region grid cell, ultimately generating a nutrient cycling function map characterizing the spatial differentiation of nutrient cycling function in the study area.

[0107] Indicators for quantifying spatial heterogeneity were extracted from the evaporation spatial feature map, specifically the absolute value of the Local Moran's I for each micro-region. This value reflects the degree of difference between the evaporation rate and its neighborhood; a higher value indicates stronger local heterogeneity. Secondly, based on the land cover classification map, a 50m × 50m moving window was used to calculate the Shannon Diversity Index (SHDI) and Patch Density (PD) for each land cover type within the window. These two landscape pattern indices quantify the richness and fragmentation of land cover types, respectively; high SHDI and moderate PD are generally associated with higher biodiversity. The calculated absolute values ​​of the Local Moran's I, Shannon Diversity Index, and Patch Density were then normalized using a minimum-maximum method to unify their numerical range to between 0 and 1. Finally, a weighted linear combination model was used to calculate the Biodiversity Support Function Index, expressed as: Biodiversity Support Function Index = × Normalized local Moran index + W2 × Normalized Shannon diversity index + W3 × Normalized patch density. Where, the weighting coefficients are... , , The regression relationship between historical biodiversity survey data and these three indicators was analyzed to determine, for example, that... =0.4、 =0.4、 =0.2. Assign values ​​to the indexes calculated for all micro-regions to generate a biodiversity support function map that quantifies the ability of each region to support biodiversity.

[0108] Two indicators, "evaporation duration" and "response intensity," were extracted from the soil vitality response curves of each micro-region. A longer "evaporation duration" indicates that the soil can provide more sustained effective moisture to vegetation roots in the early morning, which is beneficial for photosynthesis; an excessively high "response intensity" (i.e., peak evaporation rate) indicates poor soil water retention, which can easily lead to water stress in vegetation. Simultaneously, the Normalized Difference Vegetation Index (NDVI) was calculated using visible light orthophotos collected by UAVs as a direct representation of vegetation growth status and biomass. Furthermore, soil organic carbon (SOC) content data were extracted from the integrated geographic information layer. A Random Forest Regressor model was constructed to predict carbon sequestration function. This model, ensembled from hundreds of decision trees, is capable of handling complex nonlinear relationships. The model's input feature vector includes four variables: evaporation duration, response intensity, NDVI value, and SOC content. The model's output target variable is net ecosystem productivity (NEP), expressed in grams of carbon per square meter per day. The model was trained and validated using NEP data measured by eddy covariance analyzers or box cascade methods within the study area or similar ecosystems. After training, the four input feature values ​​of each micro-region within the study area are fed into the trained random forest model for prediction, yielding an estimated carbon fixation function for each micro-region. These estimated values ​​are then assigned to the corresponding grid cells to generate a carbon fixation function map.

[0109] To ensure the accuracy of the model evaluation results, 30 ground quadrats were established within the study area, covering different land cover types and functional evaluation value gradients. Within each quadrat, ground measurements corresponding to the four functions were performed. The moisture regulation function was validated by continuously monitoring soil moisture content changes using a time-domain reflectometer (TDR) probe; the nutrient cycling function was determined by collecting soil samples and conducting a net nitrogen mineralization rate incubation experiment in the laboratory; the biodiversity support function was calculated by recording the species and quantity of all plant species within the quadrat to determine the species richness index; and the carbon fixation function was evaluated by measuring the net photosynthetic rate of dominant plants using a portable photosynthesis measurement system. The measured functional values ​​from the 30 quadrats were mapped one-to-one with the predicted values ​​at corresponding locations on the four functional maps generated by the model, forming 30 sets of "predicted value-measured value" data pairs. Pearson correlation analysis was performed on the data pairs for each function, calculating the correlation coefficient r and root mean square error (RMSE) to verify the model accuracy. If the correlation is significant (e.g., r > 0.7, p < 0.05), a linear regression equation is established between the predicted and measured values: Measured value = a × Predicted value + b. This regression equation is then used to perform pixel-by-pixel linear correction on each of the originally generated function maps. That is, the predicted value of each pixel in the map is transformed using this equation, thereby eliminating systematic bias. The four corrected function maps together constitute the corrected function evaluation set.

[0110] Based on regional ecological functional zoning or referencing similar functional assessment values ​​of undisturbed ecosystems (e.g., the core area of ​​a nature reserve), a "health reference baseline" is determined for each ecological function (water regulation, nutrient cycling, biodiversity support, carbon sequestration) under each land cover type. This baseline represents the optimal level of the function under healthy conditions. Then, each micro-region grid cell in each functional map of the corrected functional assessment set is traversed, and the degree of impairment of a function in that cell is calculated using the following formula: Functional Impairment (%) = (Health Reference Baseline - Corrected Functional Value of the Cell) / Health Reference Baseline × 100%. The result of this formula represents the percentage loss of the current functional value compared to the ideal state. If the result is negative, it is recorded as 0%, indicating that the function is unimpaired or better than the reference state. This calculation is repeated for all micro-regions and all four functions, ultimately generating a functional impairment table. This table is a structured data table, with rows representing unique IDs for each microregion and columns including the microregion ID and the percentage of damage to four functions: water regulation damage, nutrient cycling damage, biodiversity support damage, and carbon fixation damage.

[0111] The Analytic Hierarchy Process (AHP) was used to determine the relative importance weights of four ecological functions (water regulation, nutrient cycling, biodiversity support, and carbon sequestration) in the regional ecosystem. Specifically, experts in ecology, economics, and management constructed a 4×4 judgment matrix, comparing each pair of functions and assigning them a scale value of 1-9 based on their relative importance. The weight coefficients for each function were obtained by calculating the largest eigenvalue and its corresponding eigenvector of the judgment matrix and normalizing the eigenvectors. (Σ) =1). Next, calculate the individual functional integrity of each micro-region using the formula: Functional integrity = 1 - Functional impairment / 100%. Then, for each micro-region, calculate its comprehensive ecological functional integrity index using a weighted summation method, using the formula: Comprehensive index = ×Moisture Integrity+ ×Nutrient integrity+ ×Biological integrity+ × Carbon fixation integrity. Finally, an ecological function integrity matrix is ​​constructed, which is an N×7 dimensional geospatial data structure, where N is the total number of microregions. Each row of the matrix corresponds to a microregion, and its columns are: microregion ID, center point X coordinate, center point Y coordinate, water regulation function integrity, nutrient cycling function integrity, biodiversity support function integrity, carbon fixation function integrity, and comprehensive ecological function integrity index. This matrix comprehensively quantifies the overall health status and internal structure of the ecosystem within the study area in a spatial form.

[0112] Preferably, step S4 includes the following steps:

[0113] Step S41: Obtain regional ecological and economic data, and construct a regional value reference system by combining the ecological function integrity matrix to obtain a functional value reference table;

[0114] Step S42: Calculate the basic damage value list of the ecological function integrity matrix based on the functional value reference table;

[0115] Step S43: Assess restoration capacity based on the ecological function integrity matrix to obtain a restoration potential distribution map;

[0116] Step S44: Based on the functional value reference table and the repair potential distribution map, predict the repair time to obtain the functional repair time sequence table;

[0117] Step S45: Calculate the dynamic changes in value based on the functional repair timeline and the basic damage value list to obtain the value decay curve;

[0118] Step S46: Perform present value conversion on the value decay curve according to the preset social discount rate to obtain the discounted damage value table;

[0119] Step S47: Based on the discounted damage value table and the restoration potential distribution map, conduct an intervention application and benefit assessment to obtain an intervention benefit assessment table;

[0120] Step S48: Generate a spatiotemporal map of ecological damage value based on the intervention benefit assessment table, the discounted damage value table, and the value decay curve.

[0121] In this embodiment of the invention, the ecological service value assessment results of the province or similar climate and economic zone where the study area is located are collected through literature research and public database queries. The focus is on extracting the annual value per unit area of ​​four functions: water regulation, nutrient cycling, biodiversity support, and carbon sequestration, expressed in yuan / hectare / year. For example, the water regulation value of a forest ecosystem in a certain province is found to be 5000 yuan / hectare / year from the "China Ecosystem Service Value Assessment Report". Next, to make these macro-value data more consistent with the actual situation of the study area, a value equivalent factor correction method is used for localization. Specifically: 1. Biomass factor correction: The average NDVI value of the study area is compared with the average NDVI value of the value source areas in the literature, and the ratio is calculated as a correction coefficient; 2. Socioeconomic factor correction: The per capita GDP and population density of the cities and counties where the study area is located are collected and compared with the national average, and the ratio is calculated. The average of the two ratios is used as a correction coefficient for social willingness to pay. The original value benchmark is multiplied by the biomass correction coefficient and the socioeconomic correction coefficient to obtain the localized functional unit value. This process is repeated for all four functions, and different unit values ​​are assigned to different land types (such as forest land and grassland) based on the land cover classification map, ultimately forming a functional value reference table. This table uses land cover type and function type as indexes and clearly lists the annual economic value per unit area for each combination. For example, the value of "forest land - water regulation function" is 5,500 yuan / hectare / year.

[0122] The process iterates through each 10m × 10m microregion in the ecological function integrity matrix. For each microregion, it obtains the integrity scores for its land cover type and four functions: water regulation, nutrient cycling, biodiversity support, and carbon sequestration. Then, it looks up the annual value per unit area for each function under the corresponding land cover type in the functional value reference table for that microregion. The damage value of a single function is calculated using the formula: Damage value per unit area = Annual value per unit area × (1 - Function integrity score) × Microregion area. For example, a microregion is forest land (0.01 hectares), with a water regulation function integrity score of 0.7. The table shows that the unit value for forest land water regulation function is 5500 yuan / hectare / year. Therefore, the damage value for the water regulation function of this microregion is 5500 × (1 - 0.7) × 0.01 = 16.5 yuan / year. This calculation is repeated for all four functions of the microregion, and the four damage values ​​are summed to obtain the total damage value of the microregion. The same calculation process is performed for all microregions within the study area, ultimately compiling a basic damage value list. The list is a detailed table, with each row representing a micro-region. The columns include the micro-region ID, land cover type, integrity of each function, damage value of each function, and total damage value.

[0123] The "response stability" index was extracted from the soil vitality response curve. This index reflects the stability of the soil microenvironment and is a core representation of intrinsic resilience. Secondly, based on the ecological function integrity matrix, the average comprehensive ecological function integrity index of all micro-regions within a 500-meter radius of each micro-region was calculated. This average quantifies the neighborhood health status at the landscape scale, representing the "source" effect from the surrounding ecosystem. Both the response stability and neighborhood health indices were normalized using a minimum-maximum normalization process. Then, a restoration potential index was calculated using a product model: Restoration Potential Index = Normalized Response Stability × Normalized Neighborhood Health. This model implies that a region's restoration potential is highest only when it possesses inherent resilience and its surrounding environment is healthy. The calculated restoration potential index was assigned to each micro-region, generating a restoration potential distribution map. The value of each pixel (between 0 and 1) in the map visually represents the magnitude of the natural restoration potential at that location.

[0124] This study predicts the time required for each damaged microregion to recover to a healthy state. The prediction is based on a logistic growth model, which assumes that the recovery process of ecological functions follows an "S"-shaped curve: slow initially, accelerating in the middle stage, and leveling off in the later stage. The model expression is: V(t) = K / (1 + ((K-V0) / V0) × exp(-r × t)). Here, V(t) is the functional value at time t, K is the total value of the function in a healthy state (obtained from a functional value reference table), V0 is the current value of the damaged function (obtained from a basic damage value list), and r is the intrinsic growth rate, representing the restoration rate. The key here is to determine the intrinsic growth rate r of each microregion and establish a direct linear relationship with the restoration potential index: r = r_max × restoration potential index, where r_max is a preset maximum restoration rate, determined by referring to literature data on similar ecosystem restoration cases. For example, the maximum annual restoration rate of a certain type of grassland ecosystem is 0.1. For each damaged function in each microregion, K, V0, and the calculated r are substituted into the logistic model. Then, a recovery target is set, such as restoring to 95% of the K value. The equation 0.95×K=K / (1+((K-V0) / V0)×exp(-r×t)) is solved to obtain the time t. This t is the predicted repair time for that function in that microregion. This calculation is performed for all damaged microregions and all functions, ultimately generating a function repair time series table. This table records the number of years required for each function in each microregion to recover to the target level.

[0125] Based on the logistic recovery model used in step S44, a function for the change of damage value over time can be derived. The damage value at time t, D(t), is equal to the total value K minus the recovered value V(t), i.e., D(t) = KV(t). Substituting the logistic model, the decay function of damage value is obtained: D(t) = KK / (1 + ((K-V0) / V0) × exp(-r × t)). For each damaged function in each micro-region of the basic damage value list, its value decay function can be completely determined using its initial damage value D0 (i.e., K-V0), total value K, and intrinsic growth rate r. Using this function, the remaining damage value of the micro-region in each year of the next 50 years from the current time (t=0) can be calculated. Connecting the damage value data points of all micro-regions over the next 50 years forms a set of curves, which is the set of value decay curves. Each curve represents the path of the evolution of the damage value of a function in a micro-region over time.

[0126] This method unifies the value of damage occurring at different points in the future to the value at the current point in time, enabling cross-period comparisons and decision-making. The Social Discount Rate (SDR) reflects society's valuation of future benefits or costs. First, a social discount rate suitable for long-term ecological projects in the region is determined; for example, based on the recommendation of the National Development and Reform Commission, it is set to 3.5%. Then, each curve in the value decay curve set is traversed, representing the future damage value sequence D(t) for each function in each micro-region. The present value calculation formula is applied to the damage value D(t) in the t-th year: PV(t) = D(t) / Calculate the present value of the damage for each year from t=1 to t=50. Sum the present values ​​of all years to obtain the Total Present Value of Damage (TPVD), i.e., TPVD = Σ[D(t) / The summation interval is t=1 to 50. This calculation is repeated for all damaged functions in all micro-regions, ultimately generating a discounted damage value table. The structure of this table is similar to the basic damage value list, but its value columns are replaced with the total present value of each functional damage and a summed present value of the total damage to the micro-regions.

[0127] An intervention library is established, containing specific remediation techniques for different types of damage (such as water imbalance and nutrient deficiency). For example, for areas with severely impaired water regulation, the intervention is "terracing"; for areas with poor nutrient cycling, the intervention is "applying more organic fertilizer." Each intervention is associated with an implementation cost per unit area (C) and a remediation acceleration coefficient (α>1). This coefficient indicates that after the intervention, the remediation rate r will become α×r. Then, several "hotspot" micro-regions with the highest total present value of damage in the discounted damage value table are identified as priority intervention targets. One or more interventions are simulated and applied to these hotspot regions. After the intervention is applied, the remediation rate r is recalculated. =α×r, and based on the new rate r Recalculate its functional recovery time and value decay curve, and then calculate the total present value of damage (TPVD) after intervention. The benefit (B) of the intervention is the reduction in the present value of the damage, B = TPVD - TPVD Calculate the cost-benefit ratio (BCR = B / C) of the intervention. Simulations are performed on different combinations of interventions, ultimately generating an intervention benefit assessment table. This table lists different intervention options, their applicable micro-regions, expected costs, expected benefits, and cost-benefit ratios, providing decision-makers with a quantitative basis for selecting the optimal remediation strategy.

[0128] All analytical results are integrated into a multi-dimensional, interactive, visualized decision support product. This map is a comprehensive data system built on a Geographic Information System (GIS) platform. Its spatial dimension visualizes the total present value of damage in micro-regions from the discounted damage value table, generating a spatial distribution base map of the current ecological damage value. Different color shades represent value levels, and the map can be overlaid with a restoration potential distribution map and priority intervention areas from the intervention benefit assessment table. Its temporal dimension is achieved through a time slider. As the user drags the time slider, the map dynamically displays the spatial distribution of damage value in any future year under a "no intervention" scenario, based on data from the value decay curve set. Its information dimension allows users to click on any micro-region, popping up an information window that details the micro-region's land cover type, the completeness of various functions, the composition of damage value, restoration potential score, predicted natural restoration time, and cost-benefit analysis results under different intervention measures. This dynamic visualization system, integrating spatial patterns, temporal evolution, and multi-level information, constitutes the final spatiotemporal map of ecological damage value.

[0129] Please see Figure 2 This is a schematic diagram of the ecological function integrity restoration curve, showing the S-shaped logistic growth pattern under three restoration potential levels:

[0130] Pattern 1: High repair potential area (r=0.15), reaching an inflection point at point A (year 5);

[0131] Pattern 2: Medium repair potential area (r=0.08), inflection point is reached at point B (year 12);

[0132] Pattern 3: Low repair potential area (r=0.04), inflection point is reached at point C (year 25).

[0133] The figure shows nonlinear recovery characteristics, with functional recovery exhibiting a sluggish pattern: slow in the early stages, accelerated in the middle stages, and gradual decline in the later stages; the intrinsic growth rate r directly affects the recovery speed and the decay of damage value.

[0134] Please see Figure 3 This diagram illustrates the damage value decay curve and its discounted present value. The solid line represents the original damage value decay curve, while the dashed line represents the present value curve after applying a 3.5% social discount rate. The diagram clearly shows the differences in value decay paths under different restoration potentials. For areas with long-term recovery potential, the discount rate significantly impacts the current economic value assessment; areas with high restoration potential should be prioritized for intervention due to their superior cost-effectiveness.

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

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

Claims

1. A method for dynamic assessment of the value of ecosystem function damage oriented towards natural resources, characterized in that, Includes the following steps: Step S1: Collect data on the formation and dissipation of morning dew in the target area, identify land cover, and obtain a land cover classification map; extract the areas where morning dew exists from the land cover classification map, and perform environmental correlation to obtain a dynamic fingerprint map of morning dew. Step S2: Divide the target area into 10m × 10m grid units to obtain a micro-region grid map; The evaporation rate of a micro-region is calculated based on the micro-region raster map and the dynamic fingerprint of morning dew; the evaporation trend is analyzed based on the micro-region evaporation rate to obtain the evaporation trend characteristics; Evaporation spatial feature maps are extracted from micro-region raster maps based on micro-region evaporation rates; healthy land is compared based on evaporation spatial feature maps and evaporation trend characteristics to obtain soil vitality response curves; Step S3: Acquire and spatially align regional geographic information and micro-regional raster maps to obtain a comprehensive geographic information layer; identify key ecological functions based on the comprehensive geographic information layer and land cover classification map to obtain ecological function indicators; perform water regulation function correlation analysis on soil vitality response curves to obtain a water regulation function map; quantify ecological function relationships based on the water regulation function map and ecological function indicators to obtain an ecological function integrity matrix; wherein, the quantification of ecological function relationships based on the water regulation function map and ecological function indicators in step S3 includes: Nutrient cycling function correlation analysis was performed on the soil vitality response curve to obtain the nutrient cycling function map; A biodiversity support function correlation analysis was performed on the soil vitality response curve to obtain a biodiversity support function map; A carbon fixation function correlation analysis was performed on the soil activity response curve to obtain a carbon fixation function map; Functional correlations of the water regulation function map, nutrient cycling function map, biodiversity support function map, and carbon fixation function map were validated to obtain a corrected functional assessment set. Calculate the functional impairment scale based on the corrective function assessment set; An ecological function integrity matrix is ​​generated based on the functional impairment scale and ecological function indicators. Step S4: Assess the damage value of the ecological function integrity matrix to obtain the value decay curve; conduct ecological restoration intervention based on the value decay curve to obtain the spatiotemporal map of ecological damage value.

2. The method for dynamic assessment of ecosystem function damage value oriented towards natural resources according to claim 1, characterized in that, Step S1 involves collecting data on the formation and dissipation of morning dew in the target area and identifying land cover, including: A baseline state map was obtained by collecting data on the target area before the formation of morning dew using a drone. Based on the baseline state diagram, time-series data acquisition was performed during the morning dew period to obtain the original time-series image; Land cover type identification is performed based on the original time-series images to obtain a land cover classification map.

3. The method for dynamic assessment of ecosystem function damage value oriented towards natural resources according to claim 1, characterized in that, Step S1 involves extracting the areas where morning dew exists from the land cover classification map and performing environmental correlation, including: Extracting temperature difference map sequences from land cover classification maps; Temperature cold spots were identified from the temperature difference map sequence to obtain a cold spot marker map sequence; The land cover classification map is subjected to reflective enhancement processing to obtain a reflective enhancement map; Extracting high-reflectivity mask sequence from reflectivity enhancement images; The morning dew area was determined based on the high reflectivity mask sequence and the cold spot marker sequence, and the morning dew distribution mask sequence was obtained. Environmental factors were correlated with the distribution mask sequence of morning dew to obtain environmental-morning dew correlation data; A dynamic fingerprint map of morning dew is generated based on the environmental-morning dew correlation data.

4. The method for dynamic assessment of ecosystem function damage value oriented towards natural resources according to claim 1, characterized in that, Step S2, which calculates the micro-region evaporation rate based on the micro-region raster map and the morning dew dynamic fingerprint spectrum, includes: Extract the time series table of morning dew area from the dynamic fingerprint of morning dew based on the micro-region raster map; Calculate the area change rate sequence based on the morning dew area time series table; The instantaneous evaporation rate is calculated from the area change rate sequence to obtain the evaporation rate sequence; Environmental factor sequences were extracted from the dynamic fingerprint of morning dew. Temperature and humidity were standardized based on the environmental factor sequence and evaporation rate sequence to obtain the standardized evaporation rate. The evaporation region correlation is performed between the standardized evaporation rate and the micro-region grid map to obtain the micro-region evaporation rate.

5. The method for dynamic assessment of ecosystem function damage value oriented towards natural resources according to claim 1, characterized in that, Step S2, which involves analyzing evaporation trends based on micro-region evaporation rates, includes: The evaporation rate of the micro-region is serialized to obtain the evaporation rate sequence. By fitting the trend of the evaporation rate sequence, the evaporation fitting function is obtained. Extract peak rate features from the evaporation fitting function; The evaporation duration is obtained by calculating the evaporation duration using the evaporation fitting function. Evaporation acceleration is estimated by fitting the evaporation function to obtain the evaporation acceleration characteristics; Based on the peak rate characteristics, evaporation duration, and evaporation acceleration characteristics, key time points are calibrated to obtain evaporation trend characteristics.

6. The method for dynamic assessment of ecosystem function damage value oriented towards natural resources according to claim 1, characterized in that, Step S2, which involves extracting the evaporation spatial feature map from the micro-region raster map based on the micro-region evaporation rate, includes: Spatial data processing is performed on the micro-region evaporation rate and the micro-region raster map to obtain the spatial point set of evaporation rate; By constructing neighborhood relationships for the spatial point set of evaporation rates, a spatial weight matrix is ​​obtained; Local hotspot analysis is performed based on the spatial point set of evaporation rate and the spatial weight matrix to obtain a local hotspot distribution map; A directional characteristic table is obtained by performing a directional analysis on the spatial point set of evaporation rates. Calculate the semivariogram curve of the spatial point set of evaporation rate; Spatial dependence was assessed on the semivariogram curve to obtain a spatial dependence index; An evaporation spatial feature map is generated based on spatial dependence indicators, local hotspot distribution maps, directional feature tables, and spatial dependence indicators.

7. The method for dynamic assessment of ecosystem function damage value oriented towards natural resources according to claim 1, characterized in that, Step S2, which compares healthy land based on the spatial evaporation characteristic map and evaporation trend characteristics, includes: Calculate the evaporation uniformity index based on the evaporation space characteristic map and the evaporation rate of the micro-region; Based on the evaporation space feature map, the boundary clarity of the micro-region raster map is evaluated to obtain the boundary feature score map; A comparative analysis of healthy land was conducted on the evaporation trend characteristics, evaporation uniformity index map, and boundary feature scoring map to obtain the health deviation matrix. Soil vitality response curves are generated based on the health deviation matrix.

8. The method for dynamic assessment of ecosystem function damage value oriented towards natural resources according to claim 1, characterized in that, Step S3 involves performing a water regulation function correlation analysis on the soil vitality response curve, including: Extracting moisture characteristic parameters from soil vitality response curves; Topographic and hydrological factors are calculated based on the integrated geographic information layer; Soil physical properties are integrated from the comprehensive geographic information layer to obtain soil moisture characteristics; Based on soil moisture characteristics and moisture feature parameters, the morning dew-infiltration relationship was constructed to obtain an infiltration capacity index map. Based on soil moisture characteristics and moisture feature parameters, the relationship between morning dew and water holding capacity was constructed to obtain a water holding capacity index map. A water storage capacity index map is constructed based on water characteristic parameters and topographic hydrological factors; Functional regulation analysis was performed on the infiltration capacity index, water holding capacity index, and water storage capacity index to obtain the water regulation function diagram.

9. The method for dynamic assessment of ecosystem function damage value oriented towards natural resources according to claim 1, characterized in that, Step S4 includes the following steps: Step S41: Obtain regional ecological and economic data, and construct a regional value reference system by combining the ecological function integrity matrix to obtain a functional value reference table; Step S42: Calculate the basic damage value list of the ecological function integrity matrix based on the functional value reference table; Step S43: Assess restoration capacity based on the ecological function integrity matrix to obtain a restoration potential distribution map; Step S44: Based on the functional value reference table and the repair potential distribution map, predict the repair time to obtain the functional repair time sequence table; Step S45: Calculate the dynamic changes in value based on the functional repair timeline and the basic damage value list to obtain the value decay curve; Step S46: Perform present value conversion on the value decay curve according to the preset social discount rate to obtain the discounted damage value table; Step S47: Based on the discounted damage value table and the restoration potential distribution map, conduct an intervention application and benefit assessment to obtain an intervention benefit assessment table; Step S48: Generate a spatiotemporal map of ecological damage value based on the intervention benefit assessment table, the discounted damage value table, and the value decay curve.

Images