A desertification ecological restoration quantification method and system based on a climate response model

Abstract

The present application relates to the technical field of environmental monitoring, and discloses a desertification ecological restoration quantification method and system based on a climate response model, which comprises the following steps: dividing a restoration area into a plurality of spatial units and identifying the treatment units in the spatial units, constructing a similarity graph, and determining an equivalent control candidate pool; fitting a climate response model for external driving, constructing a synthetic control weight, and calculating a control availability score; constructing counterfactual prediction values for the treatment units and estimating net effect values; calculating a standardized evidence quantity and making a restoration state determination; generating a determination evidence package and outputting sampling suggestions; the present application improves the comparability and stability of cross-period and cross-region controls by fitting a climate response model and extracting natural response fingerprints, performs counterfactual prediction and calculates a net effect sequence, reduces the interference of natural fluctuations and disaster impacts on acceptance conclusions, and improves the traceability, interpretability and engineering management usability of the conclusions.

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring technology, and more specifically, to a quantitative method and system for desertification ecological restoration based on a climate response model. Background Technology

[0002] Desertification ecological restoration projects primarily target desertified land, employing targeted measures to curb desertification and restore vegetation structure and ecosystem service functions. With the gradual strengthening of performance-oriented phased acceptance and long-term management systems, the management of desertification ecological restoration projects no longer requires only monitoring land ecological changes. It also necessitates quantifiable, explainable, and auditable assessments of the authenticity and effectiveness of desertification restoration, and whether the restoration results are caused by the control measures. This is crucial for supporting the allocation of desertification control resources, the iteration of restoration plans, and the determination of responsibility.

[0003] Chinese Patent CN119494471B discloses a method and system for monitoring and early warning of desert ecological environment. The method includes: determining the vegetation index and soil spectral characteristics of the monitoring area based on aerial images, and evaluating the monitoring score of the monitoring area; determining whether to conduct environmental monitoring of the monitoring area based on the relationship between the monitoring score and historical monitoring scores, wherein environmental monitoring is conducted when the monitoring score is lower than the historical monitoring score; determining the soil score and environmental score of the monitoring area based on soil and environmental data, and evaluating the desertification score of the monitoring area; obtaining historical desertification scores for adjacent time periods of the monitoring area, and determining the early warning level of the monitoring area based on the difference between the current desertification score and the historical desertification score.

[0004] However, existing technologies judge effectiveness solely based on the difference between current monitoring scores and historical scores, assuming that historical observations can serve as a reasonable comparison. They directly map changes in desertification ecological indicators under the combined effects of multiple factors to governance effects, lacking causal attribution mechanisms and evidence strength measurements. This can easily lead to misjudging natural climate fluctuations as contributions to desertification control or misjudging short-term disaster impacts as governance failures, resulting in biases in the acceptance and decision-making of desertification ecological restoration projects. This approach fails to meet the requirements of performance-oriented scenarios for assessing the net effect of desertification restoration, decomposing the contribution of measures, and ensuring cross-regional comparability. Summary of the Invention

[0005] The purpose of this invention is to provide a quantitative method and system for desertification ecological restoration based on a climate response model in order to solve the above-mentioned problems.

[0006] This invention provides a quantitative method for desertification ecological restoration based on a climate response model, comprising the following steps:

[0007] The remediation area is divided into multiple spatial units and the governance units within them are identified. The multi-source observation data of the spatial units are processed to obtain a weighted observation set and a similarity map is constructed. Based on the similarity map, an equivalent control candidate pool is determined.

[0008] For each spatial unit, a climate response model to external drivers is fitted to obtain natural response fingerprint parameters. Based on the equivalent control candidate pool and the natural response fingerprint parameters, a synthetic control weight is constructed and a control availability score is calculated.

[0009] Based on the climate response model and synthetic control weights, counterfactual predicted values ​​are constructed for the governance units and net effect values ​​are estimated.

[0010] Based on the control availability score, natural response fingerprint parameters, and net effect value, the standardized amount of evidence is calculated and the repair status is determined.

[0011] Based on the counterfactual predicted value and the net effect value, a judgment evidence package is generated and a sampling suggestion is output.

[0012] Further, determining the equivalent control candidate pool based on the similarity map includes:

[0013] The multi-source observation data collected at different times for each spatial unit is processed, and the multi-source observation data includes observation values ​​from multiple data channels;

[0014] The observations of each data channel are resampled along the time axis to be unified onto the same time axis, and the resampled observations are obtained and the observation weights are calculated. The weighted observation set includes the resampled observations and the observation weights.

[0015] For each spatial unit, extract its static attributes and long-term statistical features, and concatenate the static attributes and long-term statistical features to form a feature vector.

[0016] A similarity graph is constructed based on feature vectors. In the similarity graph, nodes are spatial units, edge weights are similarities, and the similarity is calculated using a Gaussian kernel function.

[0017] For each governance unit, multiple spatial units with the highest similarity and intervention intensity that meet the weak intervention condition are selected from the similarity map as equivalent control candidate pools. The equivalent control candidate pool includes multiple candidate units.

[0018] Furthermore, the data channels include remote sensing indices, soil moisture, surface temperature, and biological indicators;

[0019] The observation weight of the data channel of the spatial unit at time t is equal to the missing measurement mask of the data channel at time t divided by the sum of the estimated observation variance of the data channel at time t and the preset minimum normal number. The missing measurement mask is 1 when there is a valid observation of the data channel at time t, and 0 otherwise.

[0020] The static attributes of each spatial unit include geomorphic zoning, slope aspect, and soil type, while the long-term statistical characteristics include historical mean, seasonal range, and rainfall response coefficient.

[0021] The criterion for determining whether the intervention intensity meets the weak intervention condition is that the sum of the norms of the intervention vectors of the spatial unit at all times within the baseline time set is less than a preset weak intervention threshold. The intervention vectors record the intervention intensity of various governance measures received by the governance unit at different times.

[0022] Furthermore, constructing the synthetic control weights and calculating the control availability score includes:

[0023] Obtain an external driving sequence, which includes rainfall, temperature, wind speed, and evapotranspiration;

[0024] The multi-channel observations are weighted and fused according to the observation weights to obtain the composite value of the core indicators;

[0025] A climate response model is established for the composite of core indicators. The parameters are estimated using the least squares method or the weighted least squares method to obtain the natural response fingerprint parameters of each spatial unit, including the response coefficient vector, baseline bias constant, and residual variance estimate.

[0026] Calculate the fingerprint distance between the governance unit and each candidate unit in its candidate pool, calculate the synthetic control weight based on the fingerprint distance, and calculate the synthetic control sequence based on the synthetic control weight;

[0027] The control availability score is calculated based on the synthetic control sequence.

[0028] Furthermore, the climate response model represents the core index synthesis quantity of the spatial unit at time t as the inner product of the response coefficient vector and the external driving feature mapping vector plus the baseline bias constant, plus a noise term with a mean of 0. Here, the external driving feature mapping vector is the vector obtained by nonlinear feature expansion of the original external driving sequence at time t.

[0029] The fingerprint distance between the governance unit and the candidate unit is the square of the Euclidean distance between the difference of the response coefficient vectors after regularization of the two units, plus the preset balance coefficient multiplied by the square of the difference of the baseline bias constants of the two units.

[0030] The composite control weights of candidate units to governance units are calculated using a soft maximum normalization function;

[0031] The synthetic control sequence includes the synthetic control values ​​of each governance unit at multiple time points. The synthetic control value of the governance unit at time t is the result of the weighted sum of the core indicator synthetic values ​​of all candidate units in the equivalent control candidate pool of the governance unit at time t according to the synthetic control weight.

[0032] The control availability score is the geometric mean of the weighted concentration index, the baseline fit quality index, and the stability index.

[0033] Furthermore, constructing counterfactual predictive values ​​for governance units and estimating net effect values ​​includes:

[0034] Constructing counterfactual predictions: The counterfactual prediction of a governance unit at time t is equal to the weighted sum of the predicted values ​​of all candidate units in the equivalent control candidate pool at time t according to their climate response models, plus the mean of the baseline residuals of the synthetic control.

[0035] The net effect value is calculated as follows: the net effect value of the governance unit at time t is equal to the actual observed composite value of the core indicators of the governance unit at time t minus the counterfactual prediction value at time t.

[0036] For situations involving multiple overlapping measures, a measure attribution analysis is performed to decompose the intervention vector into components of each measure and establish a measure contribution model. In the measure contribution model, the net effect value is equal to the sum of the products of the marginal contribution coefficients of all measures and the intervention intensity of the corresponding measures at time t, plus the residual term.

[0037] Furthermore, the mean of the baseline residuals of the synthetic control is the weighted sum of the climate response residuals of the candidate units according to the synthetic control weights. The climate response residual is the difference between the synthetic value of the core indicator actually observed in the candidate unit and the predicted value of its climate response model.

[0038] Furthermore, calculating the standardized amount of evidence and determining the repair status includes:

[0039] The standardized evidence quantity of a governance unit is equal to the time average of the net effect value within the time set of the evaluation period of the governance unit divided by the standard error of the net effect, and then multiplied by the control availability score of the governance unit.

[0040] The standardized evidence quantity is mapped to the confidence level. The mapping adopts a sigma-shaped function, and the input is the product of a preset calibration coefficient and the standardized evidence quantity of the governance unit.

[0041] The repair status is determined by setting a set of repair statuses, including no improvement, uncertain improvement, significant improvement, and risk of regression. The determination is made jointly based on the average net effect value during the assessment period and the confidence level of the determination.

[0042] Furthermore, the rule for determining the repair status is as follows:

[0043] When the average net effect value during the evaluation period is greater than or equal to the preset positive threshold, and the confidence level is greater than or equal to the preset high confidence threshold, the repair status will be judged as significant improvement.

[0044] When the absolute value of the average net effect value during the assessment period is less than the preset positive threshold and the confidence level is less than the preset low confidence threshold, the repair status will be judged as improvement uncertainty.

[0045] When the average net effect value during the assessment period is less than or equal to the preset negative threshold and the confidence level is greater than or equal to the preset high confidence threshold, the repair status is judged as a risk of regression. The preset negative threshold is the negative value of the preset positive threshold.

[0046] In all other cases, the repair status will be judged as no improvement.

[0047] This invention provides a quantitative system for desertification ecological restoration based on a climate response model, used to implement the aforementioned quantitative method for desertification ecological restoration based on a climate response model. The system includes:

[0048] The data processing module divides the repair area into multiple spatial units and identifies the governance units within them. It processes the multi-source observation data of the spatial units to obtain a weighted observation set and constructs a similarity map. Based on the similarity map, it determines the equivalent control candidate pool.

[0049] The response model module fits a climate response model to external drivers for each spatial unit to obtain natural response fingerprint parameters. Based on the equivalent control candidate pool and the natural response fingerprint parameters, it constructs synthetic control weights and calculates control availability scores.

[0050] The net effect estimation module constructs counterfactual predicted values ​​for the governance units and estimates net effect values ​​based on the climate response model and synthetic control weights.

[0051] The status determination module calculates the standardized amount of evidence and determines the repair status based on the control availability score, natural response fingerprint parameters, and net effect value.

[0052] The evidence output module generates a judgment evidence package and outputs sampling suggestions based on the counterfactual prediction value and net effect value.

[0053] The beneficial effects of this invention are as follows: By dividing the remediation area into spatial units and constructing a similarity map based on static attributes and long-term statistical characteristics, this invention forms an equivalent control candidate pool with controlled intervention intensity, avoiding misjudgments caused by using only historical means as a comparison. By fitting a climate response model at the baseline period and extracting natural response fingerprints, the synthetic control weights are determined by fingerprint distance, improving the comparability and stability of cross-time and cross-regional comparisons. During the evaluation period, counterfactual predictions are made based on climate drivers, and the net effect sequence is calculated, reducing the interference of natural fluctuations and disaster impacts on the acceptance conclusions. By introducing observation weights and uncertainty propagation, and combining the control usability score to construct evidence strength and judgment credibility, this invention generates an evidence package containing data lineage, model parameter hashes, control weights, and net effects, and provides verification sampling suggestions, improving the traceability, interpretability, and engineering management usability of the conclusions. Attached Figure Description

[0054] Figure 1 This is a flowchart illustrating a quantitative method for desertification ecological restoration based on a climate response model, according to the present invention.

[0055] Figure 2 This is an example graph showing the calculated usability score of the present invention;

[0056] Figure 3 This is an example diagram illustrating the repair status determination of the present invention;

[0057] Figure 4 This is a module example diagram of a desertification ecological restoration quantitative system based on a climate response model according to the present invention. Detailed Implementation

[0058] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0059] Example 1:

[0060] A quantitative method for desertification ecological restoration based on climate response models, such as Figure 1 As shown, it includes the following steps:

[0061] Step 100: Divide the remediation area into several spatial units and identify the governance units within them. Process the multi-source observation data of the spatial units to obtain a weighted observation set and construct a similarity map. Determine the equivalent control candidate pool based on the similarity map.

[0062] First, the restoration area is divided into several spatial units. These spatial units refer to basic geographic units within the restoration area, defined by a certain spatial resolution. Each spatial unit corresponds to an area with a clearly defined geographic boundary. Using spatial unit division transforms continuous geographic space into discrete analytical units. This discretization provides a unified spatial benchmark for subsequent control matching and effect estimation, enabling quantitative comparisons of ecological states between different areas. The spatial unit division method involves spatially discretizing the restoration area according to a regular grid or natural geographic boundaries, based on the spatial resolution of remote sensing imagery and the actual needs of ecological restoration management, forming several spatially non-overlapping basic units. The spatial resolution of remote sensing imagery is chosen as the basis for division because remote sensing data is the primary observation data source for this method. Dividing according to its spatial resolution ensures that the observation data within each spatial unit is sufficiently representative, avoiding data sparsity due to excessively small units or spatial heterogeneity due to excessively large units. The static attributes of each spatial unit include geomorphic zoning, slope aspect, and soil type, which remain relatively stable within the observation timeframe. The purpose of extracting static attributes is to provide a stable feature basis for subsequent similarity calculations, because spatial units with similar static attributes are often closer in ecological response characteristics, thus serving as more reliable control units.

[0063] After completing the spatial unit division, it is necessary to identify the governance units within it. The governance unit refers to the spatial unit within the restoration area that actually received governance interventions, and is the main object of this method for determining the restoration status. The purpose of identifying governance units is to clarify which spatial units received governance interventions, thereby dividing the spatial units within the restoration area into two categories: governance units and potential control units, laying the foundation for subsequent control matching and effect estimation. The method for obtaining governance units is as follows: based on the implementation records of the ecological restoration project, including spatial distribution maps of governance measures, construction scope vector data, and project management ledgers, spatial overlay analysis is performed on the spatial scope of these governance measures and the spatial units to identify spatial units that received governance measures such as enclosure, reseeding, sand fixation, water replenishment, or soil improvement after the baseline period. These spatial units that received governance measures are marked as governance units. Each governance unit has a corresponding intervention vector, which records the intervention intensity of various governance measures received by the governance unit at different times. Each component of the intervention vector corresponds to a different type of governance measure, and the value of the component represents the intervention intensity of the measure at that time.

[0064] The method for obtaining intervention vectors involves first extracting the types and implementation information of governance measures received by each governance unit at various times from the implementation records of the ecological restoration project. Then, the intervention intensity of each governance measure is quantified. The reason for quantifying intervention intensity is that the implementation methods and intensities of different governance measures vary greatly; only by converting them into a unified numerical representation can quantitative comparisons and effect decomposition be performed in subsequent measure attribution analysis. Different quantitative indicators are used for quantifying intervention intensity based on the type of governance measure. The reason for using differentiated quantitative indicators is that different types of governance measures have different mechanisms of action and implementation characteristics. For example, enclosure measures mainly work through management intensity, while reseeding measures reflect intervention intensity through seeding quantity. Therefore, it is necessary to select the quantitative indicator that best reflects the intervention intensity for each measure to ensure the accuracy and comparability of the quantification results. For enclosure measures, the intervention intensity is quantified as the enclosure management intensity level. This level is comprehensively evaluated based on factors such as enclosure fence density, patrol frequency, and management personnel configuration, with a value range of 0 to 1, where 0 represents no enclosure measures and 1 represents the highest intensity of fully enclosed management. For reseeding, the intervention intensity is quantified as the ratio of the seeding rate per unit area to the standard seeding rate for that ecological type. This ratio reflects the level of investment in reseeding. When the actual seeding rate equals the standard seeding rate, the intervention intensity is set to 1; when the actual seeding rate exceeds the standard seeding rate, the intervention intensity is greater than 1. For sand fixation, the intervention intensity is quantified as the coverage density of sand fixation facilities. This density is calculated by dividing the actual coverage area of ​​sand fixation grids, sand barriers, or vegetation sand fixation belts by the total area of ​​the spatial unit, with a value ranging from 0 to 1. For water replenishment, the intervention intensity is quantified as the ratio of the replenished water volume to the average annual natural precipitation for that ecological type. This ratio reflects the degree to which artificial water replenishment enhances natural moisture conditions. For soil improvement, the intervention intensity is quantified as the ratio of the soil conditioner application rate to the recommended application rate, or the ratio of soil tillage depth to standard tillage depth. The appropriate quantitative indicator is selected based on the specific improvement method. The quantitative indicators of intervention intensity for all the above-mentioned measures were calculated using data such as construction parameters, material usage, and work area recorded in the project implementation records. When a certain type of measure was not implemented in a certain treatment unit at a certain moment, the intervention intensity corresponding to that measure was set to 0. The intervention intensity of various treatment measures for a certain treatment unit at a certain moment was arranged in the order of enclosure, reseeding, sand fixation, water replenishment, and soil improvement to form the intervention vector of that treatment unit at that moment. The dimension of the intervention vector is equal to the total number of treatment measure types.

[0065] The baseline period refers to the historical observation period prior to the large-scale implementation of the governance measures. It is used to establish a baseline of the ecosystem's natural state under no-intervention conditions. The timeframe of the baseline period is determined based on the project implementation time and the availability of historical data, typically selecting observation data from several consecutive years prior to the implementation of the governance measures. The reason for setting a baseline period is to establish a baseline of the ecosystem's natural state before the implementation of the governance measures. This baseline is used to verify the fitting quality of the synthetic control sequence during the no-intervention period and to extract long-term statistical characteristics and climate response patterns of spatial units, thus providing a reliable reference for subsequent control matching and counterfactual prediction. The evaluation period refers to the period after the implementation of the governance measures, during which the ecosystem response is observed. The timeframe of the evaluation period is determined based on project acceptance requirements and the ecological response cycle, typically selecting observation data from several consecutive years after the implementation of the governance measures. The reason for setting an evaluation period is to continuously observe changes in the ecosystem's response over a period after the implementation of the governance measures. By comparing the actual observations during the evaluation period with counterfactual predictions, the net effect is calculated, thereby quantifying the actual effectiveness of the governance measures.

[0066] Multi-source observation data collected at different times for each spatial unit are processed. This multi-source data includes observations from multiple data channels, such as remote sensing indices, soil moisture, surface temperature, and biological indicators. The reason for using multi-source observation data instead of a single data source is that ecosystem status is a multi-dimensional and comprehensive concept. A single data channel often only reflects one aspect of the ecosystem's characteristics. For example, vegetation indices mainly reflect vegetation cover and growth, soil moisture reflects water conditions, and surface temperature reflects thermal environment conditions. Only by integrating multiple data channels can the overall state of the ecosystem be comprehensively depicted, thereby improving the accuracy and comprehensiveness of the restoration effect assessment.

[0067] Time-axis resampling is performed on the observations of each data channel to unify the observations from different channels onto the same time axis, resulting in resampled observations. The reason for time-axis resampling is that multi-source observation data often have different time sampling frequencies. For example, optical remote sensing data is affected by cloud cover, leading to irregular observation times, while meteorological station data is recorded at fixed time intervals. Only by unifying these data onto the same time axis can subsequent fusion analysis and time-series modeling be performed. Missing data is handled by retaining the missing data and propagating uncertainty, rather than directly interpolating. This approach, instead of traditional interpolation, is used because interpolation artificially introduces non-existent observation information, potentially masking data quality issues and underestimating uncertainty. Retaining missing data and propagating uncertainty accurately reflects the true usability of the data. By appropriately reducing the weight of missing data periods in subsequent calculations, the problem of missing data is addressed, ensuring a more accurate reliability assessment of the final results.

[0068] Observation weights are constructed based on data quality information. The multi-source observation data includes observations from multiple data channels. The purpose of constructing observation weights is to reasonably weight each data channel according to its quality differences during the multi-source data fusion process. This ensures that higher-quality observations have a greater weight in the fusion result, while the influence of lower-quality or more uncertain observations is appropriately reduced, thereby improving the reliability and representativeness of the synthesized core indicators after fusion. Specifically, the observation weight of a data channel of a spatial unit at a certain time is equal to the missing measurement mask of that data channel at that time divided by the sum of the estimated observation variance of that data channel at that time and a preset minimum positive constant. The reciprocal of the estimated observation variance is used as the basis for weight calculation because, in statistics, the reciprocal of variance naturally represents the observation precision. The smaller the observation variance, the more precise the observation, and the greater its weight should be. This weighting method conforms to the theoretical basis of optimal linear unbiased estimation and can minimize the estimation variance during the fusion process. The missing measurement mask is used to identify the validity of the observation data. When there is a valid observation for that data channel at that time, the missing measurement mask is set to 1; otherwise, it is set to 0. The observation variance estimate is generated from information such as sensor self-checks, cloud cover, and image quality scores. It reflects the degree of uncertainty in the observation data, and its dimension is consistent with the variance of the observed values. The preset minimum normal value is used to prevent calculation anomalies caused by a denominator of 0. Its default value is one ten-thousandth of the typical value of the observation variance estimate. This value is determined through numerical stability testing to ensure numerical stability even when the observation variance estimate is close to 0, while also preventing excessive concentration of weights in low-variance data channels. The dimension of the observation weights is the reciprocal of the variance; normalization is used in subsequent weighted fusion to ensure dimensional consistency.

[0069] For each spatial unit, its static attributes and long-term statistical features are extracted. These long-term statistical features include historical mean, seasonal amplitude, and rainfall response coefficient. The static attributes and long-term statistical features are concatenated to form a feature vector. The reason for extracting long-term statistical features is that these features can characterize the ecological behavior patterns of spatial units over long time scales. The historical mean reflects the basic ecological level of the unit, the seasonal amplitude reflects its response strength to seasonal changes, and the rainfall response coefficient reflects its sensitivity to water conditions. These long-term features, together with the static attributes, constitute the comprehensive ecological fingerprint of the spatial unit, enabling similarity calculations to simultaneously consider the inherent attributes and dynamic response characteristics of the spatial unit.

[0070] A similarity graph is constructed based on feature vectors. Nodes in the similarity graph are spatial units, and edge weights represent similarity scores. The similarity is calculated using a Gaussian kernel function. The purpose of constructing the similarity graph is to organize the similarity relationships between spatial units in a graph structure. This graph structure can efficiently support subsequent candidate unit retrieval and neighborhood analysis. The Gaussian kernel function is used to calculate similarity because it has good mathematical properties; its output value exhibits a smooth exponential decay as the feature distance increases. This naturally transforms the Euclidean distance in the feature space into a similarity metric ranging from 0 to 1. Furthermore, the Gaussian kernel function is robust to outliers, preventing extreme differences in individual feature dimensions from having an excessive impact on the similarity calculation. Specifically, the similarity between any two spatial units is equal to the negative power of the natural constant. The exponent of this negative power is the square of the Euclidean distance between the feature vectors of the two spatial units divided by the square of a preset bandwidth parameter. The Euclidean distance between the feature vectors measures the degree of difference between the two spatial units in the feature space. The preset bandwidth parameter is used to control the rate at which similarity decays with distance. The larger the preset bandwidth parameter, the slower the similarity decays with distance. Its value ranges from 0.5 times to 2 times the median of the Euclidean distance of the feature vectors. The specific value is determined by the cross-validation method, so that the constructed similarity map can preserve the local neighborhood structure while avoiding excessive sparsity or excessive density.

[0071] For each governance unit, the top few spatial units with the highest similarity from the similarity map and whose intervention intensity meets the weak intervention condition are selected as equivalent control candidate pools. The equivalent control candidate pool includes multiple candidate units, with a preset candidate pool size of 20 spatial units (default). This value is determined by balancing control quality and computational efficiency. A candidate pool that is too small may lead to insufficient representativeness, while a candidate pool that is too large will increase computational burden and introduce low-similarity candidate units. The reason for selecting the spatial units with the highest similarity as candidate units is that units with high similarity are closer to the governance unit in terms of ecological characteristics and response patterns, and can more accurately simulate the natural evolution trajectory of the governance unit under no-intervention conditions. Simultaneously, the intervention intensity of the candidate units is required to meet the weak intervention condition because the core idea of ​​equivalent controls is to find units similar to the governance unit but not subject to governance measures to construct counterfactual benchmarks. If the candidate units themselves have also received strong governance intervention, they cannot truly reflect the natural state under no-intervention conditions, leading to biased net effect estimation. The criterion for meeting the weak intervention condition is that the sum of the first norms of the intervention vectors of the candidate spatial units at all times within the baseline time set is less than a preset weak intervention threshold. Here, the intervention vector represents the intensity of various governance measures received by the spatial unit at a certain moment, and the norm of the intervention vector represents the sum of the absolute values ​​of all components in the vector, used to measure the total intensity of the governance measures at that moment. The default value of the preset weak intervention threshold is 0.1 times the sum of the norms of the intervention intensities of the governance units. This threshold is determined by analyzing the intervention distribution during the baseline period to ensure that the degree of intervention received by the candidate unit during the baseline period is significantly lower than that of the governance unit, thereby meeting the basic requirements of equivalent control.

[0072] When the number of candidate units meeting the weak intervention criteria is insufficient to constitute the preset candidate pool size, a degradation strategy is adopted. The reason for adopting this strategy is that in practical applications, there may be situations where the governance measures cover a wide range or the spatial heterogeneity within the remediation area is significant, resulting in an insufficient number of candidate units that meet the strict weak intervention criteria and have high similarity. Directly abandoning the analysis of these governance units would reduce the applicability of the method. However, by gradually relaxing the constraints and recording control quality warnings, the limitations of control quality can be accurately reflected while ensuring the applicability of the method. First, the preset weak intervention threshold is gradually relaxed, increasing in steps of 0.2, 0.3, and 0.5, until the number of candidate units reaches the minimum requirement of the preset candidate pool size. The default minimum requirement is 5 spatial units. If the minimum requirement is still not met after relaxing to 0.5 times the threshold, the candidate search range is further expanded, selecting the second most similar spatial units from the similarity map and including them in the candidate pool, while simultaneously recording a control quality warning. If the candidate pool is still empty after expanding the search scope, the governance unit is marked as having no available control. In subsequent steps, counterfactual prediction and net effect estimation are not performed. Instead, the judgment result of insufficient data is directly output, and the governance unit is listed as the highest priority on-site verification object in the sampling recommendation.

[0073] Step 200: Fit a climate response model to external drivers to each spatial unit to obtain natural response fingerprint parameters. Construct synthetic control weights based on the equivalent control candidate pool and the natural response fingerprint parameters, and calculate the control availability score, as detailed below. Figure 2 As shown.

[0074] External driving sequences are acquired, including meteorological and hydrological data such as rainfall, temperature, wind speed, and evapotranspiration, as well as markers of major events such as dust storms and floods. The reason for introducing external driving sequences is that ecosystem state changes are not only influenced by governance measures but also strongly driven by climate and environmental factors. Without considering these external driving factors, it is impossible to accurately distinguish how much of the observed ecological change stems from governance measures and how much from natural fluctuations in climate conditions. By establishing a climate response model, the response pattern of each spatial unit to external driving forces can be quantified, thereby accurately separating climate effects from governance effects in subsequent counterfactual predictions.

[0075] The core indicator composite quantity is obtained by weighting and fusing multi-channel observations according to their observation weights. Specifically, the core indicator composite quantity of a spatial unit at a certain time is equal to the sum of the products of the observation values ​​of each data channel and their corresponding observation weights at that time, divided by the sum of the observation weights of all data channels at that time. This weighted fusion method ensures that the core indicator composite quantity comprehensively reflects the information from multiple data channels, and at the same time, it reasonably weights the data according to the observation quality of each channel, making the fusion result both comprehensive and reliable. As a single comprehensive ecological state indicator, the core indicator composite quantity has normalized dimensions and its value range is usually between 0 and 1, with a larger value indicating a better ecological state.

[0076] Then, a climate response model is established for the composite quantities of core indicators. The purpose of establishing the climate response model is to characterize the response relationship of the ecological state of a spatial unit to external climate drivers. This response relationship constitutes the natural response fingerprint parameters of that spatial unit. The climate response model represents the composite quantity of the core indicators of a spatial unit at a certain time as the inner product of the response coefficient vector and the external driving feature mapping vector, plus the baseline bias constant, plus a noise term with a mean of 0. The linear model form is adopted because linear models have the advantages of stable parameter estimation, high computational efficiency, and strong interpretability. At the same time, the feature mapping method can capture nonlinear response relationships while maintaining the linear structure of the model.

[0077] The input to the climate response model includes an external driving feature mapping vector, which is a high-dimensional feature representation obtained by feature expansion of the original external driving sequence. The original external driving sequence includes observed values ​​of meteorological and hydrological driving variables such as rainfall, temperature, wind speed, and evapotranspiration at the current moment and several moments prior. Introducing driving variables from the previous few moments is to capture the lagged effects of climate driving on the ecosystem, as the ecosystem's response to climate change often has a time lag; for example, vegetation growth after rainfall requires a certain amount of time to be reflected in remote sensing observations. The length of the time window for the previous few moments is determined according to the ecological type and response cycle: the default value is 3 months for grassland ecosystems, 6 months for desert ecosystems, and 12 months for forest ecosystems.

[0078] The external driving feature mapping vector is a vector obtained by nonlinearly expanding the original external driving sequence. This vector includes the original driving variables and their predefined polynomial terms and interaction terms. The reason for nonlinear feature expansion is that ecosystem responses to climate drivers often exhibit nonlinear characteristics. For example, vegetation growth tends to saturate in response to rainfall when water is plentiful, and there is an optimal temperature range for the response to temperature. Furthermore, there are interactions between different climate factors. Introducing polynomial terms and interaction terms can effectively capture these nonlinear and interactive effects within a linear model framework. Specifically, for each original driving variable, its first-order and second-order terms are calculated as predefined polynomial terms, such as the first-order and squared-order terms of rainfall, and the first-order and squared-order terms of temperature. The default order of the predefined polynomial terms is 2, which is determined through cross-validation to capture nonlinear response relationships while avoiding overfitting caused by excessively high feature dimensions. For the predefined interaction terms, pairwise products of all original driving variables are calculated, such as the product of temperature and rainfall, and the product of wind speed and evapotranspiration. The preset interaction terms include pairwise products of all original driving variables, used to capture synergistic or antagonistic effects between different climate factors. Furthermore, for major event markers such as dust storms and floods, they are included as binary features in the feature mapping vector, taking a value of 1 when the corresponding event occurs and 0 otherwise. All the above features are arranged in the order of linear terms of original driving variables, quadratic terms of original driving variables, pairwise interaction terms, and major event markers to form the external driving feature mapping vector. The dimension of this feature mapping vector is equal to the number of original driving variables multiplied by 2, plus the number of pairwise combinations of the original driving variables, plus the number of major event types. For example, when the original driving variables include four variables (rainfall, temperature, wind speed, and evapotranspiration), and the major events include two types (dust storms and floods), the dimension of the feature mapping vector is 4 times 2 plus the number of pairwise combinations of 4 (6) plus 2, equaling 16 dimensions.

[0079] The output of the climate response model is the predicted value of the composite quantity of the core indicators for the spatial unit at that moment. This predicted value equals the inner product of the response coefficient vector and the external driving feature mapping vector, plus the baseline bias constant. The response coefficient vector characterizes the spatial unit's sensitivity to each feature mapping component. The dimension of the response coefficient vector is the same as that of the external driving feature mapping vector. Each component of the response coefficient vector corresponds to a component of the feature mapping vector, representing the expected change in the composite quantity of the core indicators when that feature component changes by one unit. The baseline bias constant represents the benchmark indicator value for the spatial unit when the external driving force is at its average level. This constant reflects the inherent ecological level of the spatial unit and does not change with short-term fluctuations in the external driving force. The noise term represents random fluctuations that the model cannot explain. This noise term is assumed to follow a normal distribution with a mean of 0, and its variance is the residual variance, used to quantify the uncertainty of the model fit.

[0080] The parameters of a climate response model include the response coefficient vector, baseline bias constant, and residual variance estimate; these parameters are collectively referred to as natural response fingerprint parameters. The parameters are set through the model training process, which uses observational data from the baseline period. Baseline period data is chosen for training because governance measures have not yet been implemented on a large scale during this period, and changes in the ecological state of spatial units are mainly determined by climate-driven factors and natural evolution. Therefore, baseline period data can accurately reflect the natural response patterns of spatial units to climate-driven factors, without being affected by governance interventions.

[0081] The model training process employs least squares or weighted least squares to estimate the parameters of the climate response model. Specifically, for a given spatial unit, the observed values ​​of the composite core indicator and the corresponding external driving feature mapping vectors are collected for all times during the baseline period, forming the training dataset. The number of samples in the training dataset equals the number of time points during the baseline period, and each sample includes one observed value of the composite core indicator and one external driving feature mapping vector. When the number of time points during the baseline period is less than twice the dimension of the feature mapping vector, a regularization constraint needs to be applied to the response coefficient vector to avoid underdetermined parameter estimation. The specific method will be explained in the subsequent regularization processing section.

[0082] The reason for using the least squares method for parameter estimation is that it is a classic estimation method for linear regression models, providing unbiased parameter estimates with minimal variance when the noise term satisfies the independent and identically distributed assumption. The goal of the least squares method is to find a set of response coefficient vectors and baseline bias constants that minimizes the sum of squares of the differences between the observed values ​​and model predictions of the composite core indices at all times during the baseline period. Analytical solutions for the response coefficient vectors and baseline bias constants can be obtained by taking the derivative of the objective function and setting it to zero. Specifically, the design matrix is ​​formed by arranging the external driving feature mapping vectors at all times during the baseline period, and the observation vector is formed by arranging the observed values ​​of the composite core indices at all times during the baseline period. The response coefficient vectors and baseline bias constants can then be obtained by multiplying the pseudo-inverse of the design matrix by the observation vectors. In practical calculations, to improve numerical stability, numerical methods such as singular value decomposition or QR decomposition are used to solve for the pseudo-inverse.

[0083] When the observation quality varies across different times, weighted least squares can differentiate the fitting error at different times by weighting the observations according to their respective observation weights. This allows times with higher observation quality to play a greater role in parameter estimation, thereby improving the accuracy and stability of parameter estimation. The goal of weighted least squares is to find a set of response coefficient vectors and baseline bias constants such that the sum of the squares of the differences between the observed and predicted values ​​of the composite core indicator at all times during the baseline period is minimized by weighting the composite core indicator values ​​according to their observation weights at that time. These observation weights are derived from the observation weights calculated during the fusion of multi-channel observations in step 100. The observation variance of the composite core indicator is calculated using the variance propagation formula, and the observation weights are taken as the reciprocals of the observation variance. The solution method for weighted least squares is similar to that of ordinary least squares. It simply involves weighting the design matrix and observation vectors according to the square roots of the observation weights, and then applying the solution formula for ordinary least squares.

[0084] After estimating the response coefficient vector and baseline bias constant, the residual variance estimate is calculated. Specifically, for each time point in the baseline period, the difference between the observed value of the composite value of the core indicators and the predicted value of the climate response model is calculated to obtain the residual at that time. The residual variance estimate is equal to the sum of squares of the residuals at all times in the baseline period divided by the degrees of freedom, where the degrees of freedom equal the number of time points in the baseline period minus the number of parameters, and the number of parameters equals the dimension of the response coefficient vector plus one. The residual variance estimate is used to quantify the fitting uncertainty of the climate response model, and this estimate serves as one of the sources of model fitting uncertainty in the subsequent calculation of the net effect standard error.

[0085] The estimated natural response fingerprint parameters characterize the response features of a spatial unit to climate-driven factors. These parameters are used in subsequent control matching to measure the similarity of response patterns among different spatial units. The response coefficient vector reflects the sensitivity and response direction of the spatial unit to various climate factors, the baseline bias constant reflects the inherent ecological level of the spatial unit, and the residual variance estimate reflects the natural fluctuation of the ecological state of the spatial unit. Spatial units with similar natural response fingerprint parameters will exhibit similar ecological response trajectories under the same climate-driven conditions, and therefore can serve as control units for each other.

[0086] Based on the equivalent control candidate pool obtained in step 100 and the aforementioned natural response fingerprint parameters, a synthetic control sequence is constructed for each governance unit. Before constructing the synthetic control sequence, the response coefficient vector of the climate response model is regularized to improve the stability and repeatability of the fingerprint parameters. Regularization is performed because when the feature mapping dimension is high, ordinary least squares is prone to overfitting, causing the response coefficient vector to be overly sensitive to noise in the training data, resulting in poor repeatability of the estimated fingerprint parameters across different datasets. The regularization process uses ridge regression, adding the square of the L2 norm of the response coefficient vector as a penalty term to the least squares objective function. The weight coefficient of the penalty term is determined through cross-validation, achieving an optimal balance between prediction error and parameter complexity during the baseline period. Ridge regression, by penalizing the magnitude of the response coefficients, effectively suppresses the variance of parameter estimation, making the fingerprint parameters more stable and thus improving the reliability of subsequent control matching and counterfactual prediction. For cases with high feature mapping dimensionality, an elastic network regularization method is further adopted, and a first-norm penalty term is introduced to achieve feature selection, eliminating feature mapping components that contribute less to the response, thereby reducing the risk of overfitting and improving the interpretability of fingerprint parameters.

[0087] The fingerprint distance between the governance unit and each candidate unit in its candidate pool is calculated. Specifically, the fingerprint distance between a governance unit and a candidate unit is equal to the square of the Euclidean distance between the differences in the regularized response coefficient vectors of the two units, plus a preset balance coefficient multiplied by the square of the differences in the baseline bias constants of the two units. The purpose of calculating the fingerprint distance is to measure the similarity between the governance unit and the candidate units in the natural response fingerprint parameter space. The smaller the fingerprint distance, the more similar the response patterns of the two units to climate-driven events, and the closer their baseline ecological levels. Such candidate units can more accurately simulate the evolution trajectory of the governance unit under no-intervention conditions. The response coefficient distance and the baseline bias distance are combined to calculate the comprehensive fingerprint distance because considering only the similarity of response patterns may lead to significant differences between the selected candidate units and the governance unit at the baseline level, while considering only the closeness of the baseline level may ignore the differences in response patterns. Combining both aspects ensures that the candidate units are consistent with the governance units in multiple dimensions. The square of the Euclidean distance between the response coefficient vector differences is used to measure the degree of difference in the response patterns of the two units to external driving events. The preset balance coefficient is used to adjust the relative weight of the response coefficient distance and the baseline bias distance, so that the two types of distances have a reasonable contribution ratio in the comprehensive fingerprint distance. Its default value is 0.5. This value is determined by sensitivity analysis of the baseline matching quality, so that the similarity of response patterns and the closeness of baseline levels are balanced in the selection of controls.

[0088] The synthetic control weight is calculated based on fingerprint distance. Specifically, the synthetic control weight of a candidate unit relative to a governance unit is calculated using a soft-maximum normalization function. This synthetic control weight is equal to the negative power of the natural constant divided by the sum of the negative powers of the natural constants of all candidate units in the equivalent control candidate pool for that governance unit. The exponent of each negative power is the negative value obtained by dividing the fingerprint distance between the governance unit and the corresponding candidate unit by a preset temperature parameter. The reason for using a soft-maximum normalization function to calculate the weight is that this function can transform the fingerprint distance into a normalized weight distribution, allowing candidate units with smaller fingerprint distances to receive larger weights, while ensuring that the sum of the weights of all candidate units is 1, satisfying the mathematical requirements of weighted average. Compared to directly selecting the single candidate unit with the smallest fingerprint distance, using a weighted combination of multiple candidate units can reduce the impact of observational noise and specific fluctuations of a single candidate unit, improving the stability and representativeness of the synthetic control sequence. The preset temperature parameter is used to control the concentration of the synthetic control weight distribution. The smaller the preset temperature parameter, the more concentrated the synthetic control weight is on the few candidate units with the closest fingerprint distance. The larger the preset temperature parameter, the more uniform the synthetic control weight distribution. Its value range is 0.1 to 1 times the standard deviation of fingerprint distance in the candidate pool. The specific value is determined by minimizing the baseline fitting error, so that the synthetic control sequence has the best fitting effect with the actual observed core indicators of the treatment unit at the baseline.

[0089] The synthetic control sequence is calculated based on the synthetic control weights. This sequence includes the synthetic control value for each governance unit at each time point during the baseline period. Specifically, the synthetic control value for a governance unit at a given time point during the baseline period is equal to the weighted sum of the composite values ​​of the core indicators of all candidate units in the equivalent control candidate pool at that time, weighted by the synthetic control weights. The purpose of calculating the synthetic control sequence is to construct a virtual control unit through a weighted combination of candidate units. This virtual control unit should have a highly similar evolutionary trajectory to the governance unit during the baseline period. The synthetic control sequence serves to verify the fit between the actual observations of the governance unit during the baseline period. By comparing the degree of fit between the composite values ​​of the core indicators of the actual observations of the governance unit and the synthetic control values ​​during the baseline period, the rationality of the selected candidate units and their weight configurations is evaluated, thus providing a reliable guarantee for counterfactual predictions in the subsequent evaluation period. If the synthetic control sequence fits well with the actual observations of the governance unit during the baseline period, it indicates that the selected candidate units and their weight configurations can effectively simulate the natural evolution of the governance unit. Therefore, using the same candidate units and weights to construct counterfactual predictions during the evaluation period has high credibility.

[0090] Simultaneously, a control availability score is calculated. This score comprehensively considers three dimensions: concentration of synthetic controls, baseline fit quality, and stability. The purpose of calculating the control availability score is to quantitatively assess the reliability of the constructed synthetic control sequence. Since not all governance units can find high-quality control units, the control availability score reflects the extent to which synthetic controls can serve as credible counterfactual benchmarks, thus providing a reliable basis for subsequent net effect estimation and repair status determination. First, the weight concentration index is calculated, which is equal to the maximum value of the synthetic control weights of all candidate units in the equivalent control candidate pool for that governance unit. Then, the baseline fit quality index is calculated, which is equal to 1 minus the ratio of the root mean square error between the actual observed core indicator synthetic quantity and the synthetic control value within the baseline time set of the governance unit, divided by the standard deviation of the actual observed core indicator synthetic quantity of the governance unit. This ratio reflects the goodness of fit of the synthetic control sequence to the governance unit during the baseline period. Next, stability indices are calculated using a placebo test. Specifically, several points within the baseline period are randomly selected as pseudo-evaluation periods, and the timeframes following these points are considered pseudo-evaluation periods. The difference between the actual observed composite values ​​of core indicators in the treatment unit and the composite control values ​​during these pseudo-evaluation periods is calculated as the pseudo-net effect. The stability index is equal to 1 minus the standard deviation of the pseudo-net effect divided by the standard deviation of the baseline observations. A smaller ratio indicates less fluctuation in the composite control sequence without actual intervention, suggesting better stability of the constructed composite control. The control availability score for a given treatment unit is equal to the geometric mean of the weighted concentration index, the baseline fit quality index, and the stability index, ranging from 0 to 1. A lower control availability score indicates weaker reliability of the composite control for that treatment unit.

[0091] Step 300: Based on the climate response model and the synthetic control weights, construct counterfactual predicted values ​​for the governance units and estimate the net effect values.

[0092] A counterfactual prediction sequence is constructed, which is the expected observation sequence of the governance unit during the assessment period under the assumption of no governance intervention. To avoid double-counting climate effects, a residual correction method is used to construct the counterfactual prediction sequence. The reason for using the residual correction method is that if the climate response model predictions of the candidate units are directly used for weighted summation, although the climate driving changes during the assessment period can be reflected, the systematic differences between the candidate units and the governance units during the baseline period will be lost. These systematic differences may originate from local environmental factors or ecological processes that are not captured by the climate response model. Without correction, the counterfactual prediction will be biased. By introducing the mean of the baseline residuals for correction, the systematic differences between the candidate units and the governance units during the baseline period can be transferred to the counterfactual predictions during the assessment period while retaining the climate driving effects, thereby improving the accuracy of the counterfactual predictions. Specifically, the climate response residuals of the candidate units during the baseline period are first calculated, which is the difference between the composite quantity of the core indicators actually observed in the candidate units and their climate response model predictions. This residual reflects the systematic biases in the candidate units that are not explained by climate driving forces. Then, the climate response residuals of the candidate units are weighted and summed according to the composite control weights to obtain the mean baseline residuals of the composite control. The counterfactual predicted value of a governance unit at a certain moment in the assessment period is equal to the weighted sum of the climate response model predicted values ​​of all candidate units in the equivalent control candidate pool of that governance unit at that moment, plus the mean baseline residuals of the composite control. The climate response model predicted value of a candidate unit at that moment is equal to the inner product of the candidate unit's response coefficient vector and the external driving feature mapping vector at that moment, plus the candidate unit's baseline bias constant. This construction method ensures that the counterfactual prediction sequence reflects both the climate driving changes in the assessment period and preserves the systematic differences between candidate units and governance units in the baseline period, while avoiding the double inclusion of climate effects. It should be noted that the counterfactual prediction sequence and the synthetic control sequence in step 200 differ in their construction methods. The synthetic control sequence directly uses the actual observed values ​​of the candidate units to sum up in a weighted manner and is only used for baseline period fitting verification. In contrast, the counterfactual prediction sequence uses the climate response model prediction values ​​of the candidate units to sum up in a weighted manner and adds baseline residual correction, and is used for prediction of the no-intervention scenario during the assessment period. This differentiated treatment ensures that the climate-driving effects and governance measures effects can be accurately separated during the assessment period.

[0093] The net effect sequence, i.e., the sequence of gains brought about by governance measures, is calculated, and includes multiple net effect values. Specifically, the net effect value of a governance unit at a certain point in the assessment period is equal to the actual observed composite value of the core indicators of that governance unit at that point minus the counterfactual predicted value at that point. The reason for defining the net effect using the difference between the actual observed value and the counterfactual predicted value is that the counterfactual predicted value represents the expected ecological state of the governance unit under the assumption of no governance intervention, while the actual observed value is the actual ecological state after the implementation of governance measures. The difference between the two is precisely the net gain or net loss brought about by the governance measures. This difference method can effectively eliminate the influence of climate-driven and natural evolution, thereby accurately separating the true effect of the governance measures. A positive net effect value indicates that the governance measures have produced positive effects, and a negative net effect value indicates that the actual performance after governance is lower than the expected value under the condition of no intervention. The net effect sequence includes the net effect values ​​of the governance unit at all points in the assessment period, used to characterize the evolution of the governance effect over time.

[0094] For situations involving multiple overlapping measures, further attribution analysis is conducted. The intervention vector is decomposed into individual measure components, including enclosure, reseeding, sand fixation, water replenishment, and soil improvement, and a measure contribution model is established. The reason for conducting measure attribution analysis is that in actual ecological restoration projects, multiple governance measures are often implemented simultaneously. Simply calculating the overall net effect cannot answer how much each measure contributes individually. Measure attribution analysis can decompose the overall net effect into specific measures, thus providing a quantitative basis for optimizing governance schemes and resource allocation. This measure contribution model represents the net effect value of a governance unit at a certain time as the weighted sum of the marginal contributions of each measure plus a residual term. The measure contribution model is established using a linear weighted sum, based on the assumption of marginal effect additivity, which assumes that the contributions of each measure can be approximately separated. This assumption provides reasonable attribution results when the interaction between measures is not significant, while maintaining the simplicity and interpretability of the model. Specifically, the net effect value equals the sum of the products of the marginal contribution coefficients of all measures and the intervention intensity of the corresponding measure at that time, plus a residual term. The total number of measure types represents the number of different governance measures involved in the analysis. The intervention intensity of each measure represents the magnitude of the effect of implementing that measure at that moment. The marginal contribution coefficient represents the net effect gain brought about by a unit of intervention intensity, characterizing the marginal utility of the measure. The residual term represents the net effect fluctuations that the measure contribution model cannot explain.

[0095] In parameter estimation, a segmented constraint strategy is adopted based on the different characteristics of the measure type and the assessment period. The reason for using this strategy is that different types of governance measures have different ecological response cycles and effect characteristics. Applying uniform constraints to all measures and all time periods might lead to estimation results that violate ecological common sense or produce unreasonable negative effects. For long-term benefit measures such as enclosure and soil improvement, a non-negative constraint is applied to their marginal contribution coefficient in the latter half of the assessment period, reflecting that the positive benefits of these measures mainly manifest in the medium to long term. For measures that may cause short-term disturbances, such as water replenishment and sand fixation, no sign constraint is applied in the first half of the assessment period, allowing the marginal contribution coefficient to be negative to capture short-term negative effects. A non-negative constraint is applied in the latter half of the assessment period to reflect the positive benefits after the disturbance subsides. For measures with moderate response cycles, such as reseeding, a weak non-negative constraint is applied throughout the assessment period, allowing the marginal contribution coefficient to be negative within a certain range but imposing a penalty term. The weight of the penalty term is determined through calibration using historical project data. The above-mentioned segmented constraint strategy is implemented through a constraint optimization algorithm to ensure that the estimated marginal contribution coefficient minimizes the fitting error under the constraint conditions, while improving the rationality and credibility of the attribution of measures.

[0096] Step 400: Based on the control availability score, natural response fingerprint parameters, and net effect value, calculate the standardized amount of evidence and determine the repair status, specifically as follows: Figure 3 As shown.

[0097] The standardized evidence quantity (SEF) is calculated by comprehensively considering the statistical significance of the net effect and the usability of the control. The purpose of calculating the SEF is to integrate the three dimensions of net effect estimation—size, uncertainty, and control quality—into a unified strength of evidence index. This index reflects the reliability of the net effect estimate and provides a quantitative basis for subsequent remediation status determination. Specifically, the SEF for a given governance unit is equal to the time average of the net effect values ​​within the assessment period of that governance unit, divided by the standard error of the net effect, and then multiplied by the control usability score of that governance unit. Using the net effect mean divided by the standard error is similar to the test statistic in statistical testing, simultaneously considering effect size and estimation precision. When the net effect mean is large and the standard error is small, the SEF is large, indicating that the observed net effect has strong statistical significance. Multiplying by the control usability score further considers the quality of the synthetic control. When the control quality is low, even if the net effect is statistically significant, its strength of evidence should be reduced to avoid overly definitive conclusions based on low-quality controls. The time average of the net effect values ​​is the arithmetic mean of the net effect values ​​at all times during the assessment period. The calculation of the net effect standard error needs to consider the propagation of multiple sources of uncertainty, including observational uncertainty, model fitting uncertainty, and synthetic control construction uncertainty. Observational uncertainty originates from the observation weights of each data channel in step 100, and the observational variance of the composite core indicator is calculated using the variance propagation formula for weighted fusion. Model fitting uncertainty originates from the residual variance estimation of the climate response model in step 200. Synthetic control construction uncertainty is calculated by weighting and summing the observation variances of candidate units according to the squares of the synthetic control weights. The net effect standard error is equal to the square root of the sum of the variances of the above three types of uncertainty divided by the number of time points in the assessment period, representing the overall uncertainty of the net effect mean estimate. The control availability score, derived from the calculation results of step 200, is used to weight and adjust the standardized evidence quantity, reflecting the reliability of the synthetic control. It should be noted that the standardized evidence quantity is not a strict statistical test quantity, but rather a scoring index of the overall evidence strength. Its value reflects the reliability and significance of the net effect estimate, but does not directly correspond to a specific significance level or confidence interval.

[0098] Then, the standardized evidence quantity is mapped to the confidence level. Specifically, the confidence level of a governance unit is equal to 1 divided by 1 plus a negative power of the natural constant, where the exponent of the negative power is the negative product of the preset calibration coefficient and the standardized evidence quantity of that governance unit. This mapping uses a sigma function to map the range of standardized evidence quantity from negative infinity to positive infinity to the confidence level range of 0 to 1. The reason for using a sigma function is that it has the characteristic of smooth monotonically increasing, which can transform the unbounded statistical quantity of standardized evidence quantity into a bounded probabilistic index, making the confidence level range uniform between 0 and 1, which is convenient for subsequent threshold determination and result interpretation. At the same time, the sigma function has a large slope when the standardized evidence quantity is close to 0, which can effectively distinguish small differences in evidence strength, while it tends to saturate when the absolute value of the standardized evidence quantity is large, avoiding extreme values ​​from having an excessive impact on the confidence level. The preset calibration coefficient is used to adjust the steepness of the mapping curve. Its default value is 1. This value can be calibrated and adjusted by historical acceptance samples or expert annotations. The specific adjustment method is to determine the optimal preset calibration coefficient by maximizing the area under the receiver operating characteristic curve that maximizes the consistency between the judgment confidence and the actual acceptance results, so that the judgment confidence and the actual acceptance results have good consistency.

[0099] Next, the restoration status is determined. A set of restoration statuses is established, including no improvement, uncertain improvement, significant improvement, and risk of regression. The reason for setting four restoration statuses is that ecological restoration effect assessment not only needs to determine whether improvement has occurred, but also needs to distinguish the significance of the improvement and the credibility of the determination, while also identifying areas at risk of ecological degradation. A joint determination is made based on the average net effect value during the assessment period and the credibility of the determination. A joint determination method is used instead of relying solely on the net effect value because the size of the net effect value only reflects the direction and magnitude of the effect, while the credibility of the determination reflects the reliability of this estimate. Only by considering both dimensions simultaneously can a scientific and reasonable determination be made, avoiding the use of unreliable estimates as the basis for judgment. The determination rules are as follows:

[0100] When the average net effect value during the evaluation period is greater than or equal to the preset positive threshold, and the confidence level is greater than or equal to the preset high confidence threshold, the repair status is judged as significantly improved.

[0101] When the absolute value of the average net effect value during the assessment period is less than the preset positive threshold and the confidence level is less than the preset low confidence threshold, the repair status will be judged as improvement uncertainty.

[0102] When the average net effect value during the assessment period is less than or equal to a preset negative threshold, and the confidence level is greater than or equal to a preset high confidence threshold, the repair status is judged as a risk of regression. The preset negative threshold is the negative value of a preset positive threshold.

[0103] In all other cases, the repair status will be judged as no improvement.

[0104] The preset positive threshold is the preset net effect threshold, with a default value of 0.2 times the baseline standard deviation of the core indicator's composite quantity. This value is determined by analyzing the minimum detectable ecological improvement in historical projects, ensuring that the determined significant improvement status has practical ecological significance. The default value for the preset low confidence threshold is 0.3, and the default value for the preset high confidence threshold is 0.7. These confidence thresholds are determined by analyzing the ternary loci of the confidence distribution, ensuring reasonable differentiation between different restoration statuses in terms of confidence. All of the above preset thresholds are configurable parameters and can be set differently according to ecological type. Different ecological types, such as grasslands, deserts, and wetlands, can have different combinations of preset thresholds set according to their ecological response characteristics and management needs.

[0105] Step 500: Generate a judgment evidence package based on the counterfactual prediction value and the net effect value, and output sampling suggestions.

[0106] A judgment evidence package is generated, outputting for each governance unit the actual observation sequence, counterfactual prediction sequence, net effect sequence, synthetic control weight, control availability score, standardized evidence quantity, judgment confidence, remediation status, and marginal contribution coefficient of the measures, as well as data lineage metadata, including data source identifier, time window, model version, and parameter hash value. The purpose of generating the judgment evidence package is to provide a complete chain of evidence and traceable information for the judgment results of each governance unit, enabling the judgment results to not only provide the final remediation status conclusion but also demonstrate the entire intermediate process and basis for reaching that conclusion, thereby supporting the verification, review, and response to challenges. Outputting data lineage metadata ensures the reproducibility of the judgment results, clearly distinguishing between different versions of data or models, and avoiding confusion caused by inconsistencies in results due to data or method updates. Among them, the actual observation sequence is the time series of the composite quantity of the core indicators of the governance unit during the baseline period and the evaluation period, the counterfactual prediction sequence is the time series of the counterfactual prediction value of the governance unit during the evaluation period, the net effect sequence is the time series of the net effect value of the governance unit during the evaluation period, and the composite control weight is the weight vector of each candidate unit in the equivalent control candidate pool of the governance unit. These pieces of information together constitute a complete chain of evidence for judgment, ensuring the traceability and verifiability of the judgment results.

[0107] Output spatially distributed repair status layer, project-level summary net effect and judgment confidence, and ranking of measures' contributions.

[0108] Sampling recommendations are generated for areas requiring further verification. For governance units deemed uncertain in terms of improvement and with low confidence levels, as well as those deemed at risk of regression, sampling priority scores are calculated. The purpose of generating sampling recommendations is to combine remote sensing assessments with ground verification to form a closed-loop validation mechanism. For areas with low confidence levels or ecological risks, targeted on-site sampling improves the reliability of assessments or allows for timely problem detection, thereby enhancing the accuracy and credibility of the overall assessment results. Specifically, the sampling priority score for a governance unit is equal to 1 minus the confidence level of that unit's assessment, multiplied by the unit's preset risk sensitivity index. The confidence level is derived from the calculation result in step 400. The preset risk sensitivity index is calculated by comprehensively considering factors such as wind erosion sensitivity, soil erosion sensitivity, and population exposure. It reflects the potential harm caused by ecological degradation in the governance unit, with a value ranging from 0 to 1. The specific calculation method involves normalizing each sensitivity factor and then summing them according to preset weights. The default preset weights for wind erosion sensitivity, soil erosion sensitivity, and population exposure are 0.4 and 0.2, respectively. These preset weights are determined based on the regional ecological risk characteristics and management priorities. A higher sampling priority score indicates that the governance unit requires priority for on-site verification to improve the reliability of the assessment.

[0109] Based on the sampling priority scores, the samples are sorted from high to low, and a list of suggested UAV routes and ground retest points is output for subsequent closed-loop improvement of judgment reliability.

[0110] Example 2

[0111] See Figure 4 As shown, a system for quantifying desertification ecological restoration based on a climate response model is provided. This system stores computer-readable instructions, which, when read, can execute the aforementioned method for quantifying desertification ecological restoration based on a climate response model. The system includes:

[0112] The data processing module 101 divides the repair area into several spatial units and identifies the governance units within them. It processes the multi-source observation data of the spatial units to obtain a weighted observation set and constructs a similarity map. Based on the similarity map, it determines the equivalent control candidate pool.

[0113] The response model module 102 fits a climate response model to external drivers for each spatial unit to obtain natural response fingerprint parameters, constructs synthetic control weights based on the equivalent control candidate pool and the natural response fingerprint parameters, and calculates the control availability score.

[0114] Net effect estimation module 103 constructs counterfactual prediction values ​​for the governance unit and estimates net effect values ​​based on the climate response model and synthetic control weights;

[0115] The status determination module 104 calculates the standardized amount of evidence and determines the repair status based on the control availability score, natural response fingerprint parameters and net effect value.

[0116] The evidence output module 105 generates a judgment evidence package and outputs sampling suggestions based on the counterfactual prediction value and net effect value.

[0117] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A quantitative method for desertification ecological restoration based on a climate response model, characterized in that, Includes the following steps: The restoration area is divided into multiple spatial units and the governance units within them are identified. The multi-source observation data of the spatial units are resampled along the time axis and the observation weights are calculated. Static attributes and long-term statistical features are extracted to form feature vectors. A similarity graph is constructed based on the feature vectors. The nodes in the similarity graph are spatial units, and the edge weights are similarities. The similarity is calculated using a Gaussian kernel function. Multiple spatial units with the highest similarity and intervention intensity that meet the weak intervention condition are selected from the similarity graph as equivalent control candidate pools. External driving sequences are acquired, and multi-channel observations are weighted and fused according to observation weights to obtain the composite quantity of core indicators. A climate response model is established for the composite quantity of core indicators of each spatial unit, and parameters are estimated using the least squares method or weighted least squares method to obtain natural response fingerprint parameters. The natural response fingerprint parameters include response coefficient vector, baseline bias constant, and residual variance estimate. Based on the equivalent control candidate pool and the natural response fingerprint parameters, the fingerprint distance between the governance unit and each candidate unit in its candidate pool is calculated. The composite control weight is calculated based on the fingerprint distance. The composite control sequence is calculated based on the composite control weight. The control availability score is calculated based on the composite control sequence. Based on the climate response model and the composite control weights, counterfactual predicted values ​​are constructed for the governance units and net effect values ​​are estimated. The counterfactual predicted value of the governance unit at time t is equal to the weighted sum of the predicted values ​​of all candidate units in the equivalent control candidate pool at time t according to their climate response model, plus the mean of the baseline residuals of the composite control. The net effect value of the governance unit at time t is equal to the actual observed composite amount of the core indicators of the governance unit at time t minus the counterfactual predicted value at time t. Based on the control availability score, natural response fingerprint parameters, and net effect value, a standardized amount of evidence is calculated and a repair status is determined. The standardized amount of evidence for a governance unit is equal to the time average of the net effect value within the evaluation period of the governance unit divided by the standard error of the net effect, and then multiplied by the control availability score of the governance unit. The standardized amount of evidence is mapped to a decision confidence level. The repair status is determined based on the average net effect value during the evaluation period and the decision confidence level to obtain repair statuses including no improvement, uncertain improvement, significant improvement, and risk of regression. Based on the counterfactual predicted value and the net effect value, a judgment evidence package is generated and a sampling suggestion is output.

2. The method for quantifying desertification ecological restoration based on a climate response model according to claim 1, characterized in that, The data channels include remote sensing indices, soil moisture, surface temperature, and biological indicators; The observation weight of the data channel of the spatial unit at time t is equal to the missing measurement mask of the data channel at time t divided by the sum of the estimated observation variance of the data channel at time t and the preset minimum normal number. The missing measurement mask is 1 when there is a valid observation of the data channel at time t, and 0 otherwise. The static attributes of each spatial unit include geomorphic zoning, slope aspect, and soil type, while the long-term statistical characteristics include historical mean, seasonal range, and rainfall response coefficient. The criterion for determining whether the intervention intensity meets the weak intervention condition is that the sum of the norms of the intervention vectors of the spatial unit at all times within the baseline time set is less than a preset weak intervention threshold. The intervention vectors record the intervention intensity of various governance measures received by the governance unit at different times.

3. The method for quantifying desertification ecological restoration based on a climate response model according to claim 1, characterized in that, The climate response model represents the core index synthesis quantity of the spatial unit at time t as the inner product of the response coefficient vector and the external driving feature mapping vector, plus the baseline bias constant, plus a noise term with a mean of 0. The external driving feature mapping vector is the vector obtained by nonlinear feature expansion of the original external driving sequence at time t. The fingerprint distance between the governance unit and the candidate unit is the square of the Euclidean distance between the difference of the response coefficient vectors after regularization of the two units, plus the preset balance coefficient multiplied by the square of the difference of the baseline bias constants of the two units. The composite control weights of candidate units to governance units are calculated using a soft maximum normalization function; The synthetic control sequence includes the synthetic control values ​​of each governance unit at multiple time points. The synthetic control value of the governance unit at time t is the result of the weighted sum of the core indicator synthetic values ​​of all candidate units in the equivalent control candidate pool of the governance unit at time t according to the synthetic control weight. The control availability score is the geometric mean of the weighted concentration index, the baseline fit quality index, and the stability index.

4. The method for quantifying desertification ecological restoration based on a climate response model according to claim 1, characterized in that, For situations involving multiple overlapping measures, a measure attribution analysis is performed to decompose the intervention vector into components of each measure and establish a measure contribution model. In the measure contribution model, the net effect value is equal to the sum of the products of the marginal contribution coefficients of all measures and the intervention intensity of the corresponding measures at time t, plus the residual term.

5. The method for quantifying desertification ecological restoration based on a climate response model according to claim 1, characterized in that, The mean baseline residual of the synthetic control is the weighted sum of the climate response residuals of the candidate units according to the synthetic control weights. The climate response residual is the difference between the synthetic value of the core indicator actually observed in the candidate unit and the predicted value of its climate response model.

6. The method for quantifying desertification ecological restoration based on a climate response model according to claim 1, characterized in that, The standardized evidence quantity is mapped to the confidence level. The mapping adopts a sigma-shaped function, and the input is the product of a preset calibration coefficient and the standardized evidence quantity of the governance unit.

7. The method for quantifying desertification ecological restoration based on a climate response model according to claim 1, characterized in that, The rules for determining the repair status are as follows: When the average net effect value during the evaluation period is greater than or equal to the preset positive threshold, and the confidence level is greater than or equal to the preset high confidence threshold, the repair status will be judged as significant improvement. When the absolute value of the average net effect value during the assessment period is less than the preset positive threshold and the confidence level is less than the preset low confidence threshold, the repair status will be judged as improvement uncertainty. When the average net effect value during the assessment period is less than or equal to the preset negative threshold and the confidence level is greater than or equal to the preset high confidence threshold, the repair status is judged as a risk of regression. The preset negative threshold is the negative value of the preset positive threshold. In all other cases, the repair status will be judged as no improvement.

8. A quantitative system for desertification ecological restoration based on a climate response model, used to implement the quantitative method for desertification ecological restoration based on a climate response model as described in any one of claims 1-7, characterized in that, The system includes: The data processing module divides the repair area into multiple spatial units and identifies the governance units within them. It processes the multi-source observation data of the spatial units to obtain a weighted observation set and constructs a similarity map. Based on the similarity map, it determines the equivalent control candidate pool. The response model module fits a climate response model to external drivers for each spatial unit to obtain natural response fingerprint parameters. Based on the equivalent control candidate pool and the natural response fingerprint parameters, it constructs synthetic control weights and calculates control availability scores. The net effect estimation module constructs counterfactual predicted values ​​for the governance units and estimates net effect values ​​based on the climate response model and synthetic control weights. The status determination module calculates the standardized amount of evidence and determines the repair status based on the control availability score, natural response fingerprint parameters, and net effect value. The evidence output module generates a judgment evidence package and outputs sampling suggestions based on the counterfactual prediction value and net effect value.

Images