An ecosystem service and landscape ecological risk coupling partitioning method based on XGBoost-SHAP driving threshold identification

Abstract

The method for coupling partition of ecosystem services and landscape ecological risk based on XGBoost-SHAP driving threshold identification belongs to the technical field of ecological system comprehensive evaluation and land space fine management. The method solves the problems that the existing coupling partition boundary relies on experience division, lacks quantitative identification of driving threshold, and is difficult to depict the critical interval of coupling state transition. The technical scheme is as follows: preprocessing multi-source data to construct a grid database, calculating the ecosystem service value and the landscape ecological risk to obtain the coupling coordination degree grade; taking the coupling coordination degree as the dependent variable to construct an XGBoost regression model to identify the key driving factors, and extracting the critical threshold value through the SHAP method; combining the coupling coordination degree grade and the driving factor threshold interval to construct a two-dimensional partition rule, completing the ecological safety partition and integrating the patches. The method realizes the quantitative of the critical interval of the coupling state, improves the scientificity of the ecological safety partition, and is suitable for regional ecological evaluation and land space fine management.

Description

Technical Field

[0001] This invention belongs to the field of integrated ecosystem assessment and refined land space management technology, specifically involving an ecological security zoning method that integrates ecosystem service value assessment, landscape ecological risk characterization, coupling coordination degree calculation, and interpretable machine learning-driven threshold identification. Background Technology

[0002] Against the backdrop of rapid regional urbanization and continuous adjustments in land use structure, ecosystem service supply capacity and landscape ecological risk levels exhibit significant spatial differentiation characteristics. Ecosystem service value (ESV) measures the comprehensive contribution of ecosystems to human society, while landscape ecological risk (ERI) characterizes the vulnerability and stability of landscape patterns under disturbance. Analyzing the interaction between ESV and ERI by constructing a coupling coordination degree model has become an important technical approach for regional ecological security assessment. Existing technologies typically employ four-quadrant models or coupling coordination degree classification methods to categorize the relationship between ecosystem services and ecological risk. While these methods can intuitively reflect the coordination status of different regions, the zoning boundaries are mostly based on statistical means, quantiles, or grade intervals, primarily focusing on the static expression of the coupling pattern. They do not incorporate driving factors to identify the nonlinear response characteristics of coupling coordination degree for threshold identification, making it difficult to characterize the critical intervals where coupling states transition. This limits their effectiveness in determining zoning boundaries and constructing refined management rules. Therefore, it is necessary to introduce a driving threshold identification mechanism based on coupling coordination analysis to construct an ecological security zoning method capable of quantifying critical response intervals and optimizing zoning boundaries. Summary of the Invention

[0003] The purpose of this invention is to address the problems of existing ecosystem service and landscape ecological risk coupling zoning processes, which rely on experience-based boundary delineation and lack quantitative identification of driving thresholds. This invention provides an ecosystem service and landscape ecological risk coupling zoning method based on XGBoost-SHAP-driven threshold identification to achieve quantitative expression of the critical interval of the coupling state, thereby improving the scientific rigor and feasibility of ecological security zoning. To achieve the above objective, this invention adopts the following technical solution.

[0004] A coupled zoning method for ecosystem services and landscape ecological risk based on XGBoost-SHAP-driven threshold identification includes the following steps: Step S1: Determine the study area, study period and evaluation unit scale, collect land use data, ecological environment factor data and socio-economic factor data, perform projection unification, resampling and standardization processing on multi-source data, and construct a spatially consistent raster database. Step S2: Calculate the ecosystem service value of each evaluation unit based on the equivalent coefficient of ecosystem service value per unit area, and construct a landscape ecological risk index model to obtain the spatial distribution results of ecosystem services and ecological risks in each unit. Step S3: Standardize the ecosystem service value and landscape ecological risk, construct a coupling coordination degree model, and calculate the coupling coordination degree level of each evaluation unit; Step S4: Construct an XGBoost regression model with coupling coordination degree as the dependent variable and natural factors and human activity factors as independent variables to identify key driving factors affecting changes in coupling coordination. Step S5: Use the SHAP method to interpret the model output, calculate the contribution value of each driving factor and construct the feature dependency relationship, and identify the nonlinear response range and critical threshold of the key driving factors. Step S6: Combine the coupling coordination degree level with the threshold range of the driving factor to construct a two-dimensional partitioning rule, divide the study area into ecological security partitions, and obtain the final ecological partitioning results.

[0005] Specifically, step S1 involves preprocessing the multi-source data, including unified format conversion, spatial reference system unification, resampling, and standardization of data from different sources, to ensure spatial consistency and comparability among various types of data.

[0006] First, vector and raster data of different formats are converted into a unified raster data format and projected onto the same spatial coordinate system. Secondly, resampling is performed on data with different spatial resolutions, with a preferred resampling resolution of 1 km × 1 km. Subsequently, the continuous driving factor data is subjected to extreme value standardization, and the standardization formula is as follows: in, The original variable value, For the standardized variable values, and These are the minimum and maximum values ​​of the sample, respectively. Finally, a database of ecosystem services and driving factors with grid cells as the basic evaluation unit is constructed.

[0007] Specifically, step S2 involves constructing an ecosystem service value model and a landscape ecological risk index model.

[0008] The formula for calculating the value of ecosystem services is: in, Let k be the land use area of ​​the kth type within the i-th grid cell. Let m be the service value coefficient per unit area for the kth land use category, and m be the number of land use types. The service value equivalent per unit area is adjusted using the regional grain output value, and the adjustment formula is as follows: in, Let be the average price of the t-th crop. Let n be the yield per unit area, and n be the number of crop varieties.

[0009] The formula for calculating the landscape ecological risk index is: in, , For landscape disturbance degree, For landscape vulnerability, This represents the total area of ​​the unit.

[0010] Specifically, in step S3, a model for the coupling and coordination of ecosystem services and landscape ecological risks is constructed.

[0011] First, the ecosystem service value and landscape ecological risk index are standardized:

[0012] Secondly, calculate the coupling degree:

[0013] Subsequently, the comprehensive development index was calculated:

[0014] Finally, the formula for calculating the coupling coordination degree is:

[0015] The D value is used to reflect the level of coordinated development between ecosystem services and landscape ecological risks.

[0016] Specifically, in step S4, an XGBoost regression model is constructed, with the coupling coordination degree D as the dependent variable and natural factors and human activity factors as independent variables to establish a prediction model.

[0017] The prediction function is:

[0018] The objective function of the model is:

[0019] in, This is used to control model complexity.

[0020] Specifically, in step S5, the SHAP method is used to interpret the model output results and identify key driving factors and their critical threshold ranges.

[0021] The formula for calculating the SHAP value of feature i is:

[0022] Constructing feature dependency functions:

[0023] By calculating the second derivative of the function, when the following condition is met:

[0024] The corresponding variable interval is then defined as the driving threshold interval.

[0025] Specifically, in step S6, a two-dimensional ecological security zoning system is constructed based on the coupling coordination degree level and the threshold range of key driving factors.

[0026] First, the basic coordination level is divided according to the coupling coordination degree D value. When D < 0.35, it is divided into the severely unbalanced development zone; when 0.35 ≤ D < 0.5, it is divided into the relatively unbalanced development zone; when 0.5 ≤ D < 0.65, it is divided into the barely balanced development zone; when 0.65 ≤ D < 0.8, it is divided into the moderately balanced development zone; and when D ≥ 0.8, it is divided into the highly balanced development zone.

[0027] Secondly, based on the critical threshold of the key driving factors identified in step S5 Construct the sensitive region discrimination interval. Let the critical threshold of the k-th driving factor be... Then its sensitivity decision interval is defined as: That is, take the critical threshold as the center and expand the range by 10% in both the front and back.

[0028] Construct a driver-sensitive discrimination function: When a certain evaluation unit has at least one key driving factor that satisfies When this occurs, the unit is determined to be in a drive threshold sensitive state.

[0029] Subsequently, the coupling coordination level and the driving threshold sensitive state are superimposed for analysis: (1) When D < 0.35, it is classified as an ecologically fragile area; (2) When D≥0.8 and all key driving factors are not in the sensitive interval, it is designated as a coordinated conservation area; (3) When 0.35 ≤ D < 0.65 and there exists at least one key driving factor that satisfies At that time, it was divided into sensitive areas; Further, the data is further subdivided based on the dominant driving factor category that enters the sensitive region: If the sensitive factor belongs to the natural ecological factor, it is classified as a vegetation-dominated sensitive area; If the sensitive factor belongs to the socio-economic category, it is classified as a socially dominant sensitive area; (4) The remaining evaluation units are divided into stable transition zones.

[0030] Finally, connectivity analysis and patch integration are performed on spatially adjacent grid units of the same type to eliminate isolated patches and generate the final ecological security zoning map. Attached Figure Description

[0031] Figure 1 This is a flowchart illustrating the overall technical process of the ecosystem service and landscape ecological risk coupled zoning method based on XGBoost-SHAP-driven threshold identification described in this invention. Detailed Implementation

[0032] The present invention will be further described below with reference to the accompanying drawings and specific embodiments. It should be understood that the following embodiments are for illustrative purposes only and are not intended to limit the scope of protection of the present invention. Equivalent substitutions or modifications made by those skilled in the art without departing from the concept of the present invention should fall within the scope of protection of the present invention. Example 1

[0033] This embodiment takes Shandong Province as the research area to illustrate a method for coupled zoning of ecosystem services and landscape ecological risks based on XGBoost-SHAP-driven threshold identification.

[0034] The study area is Shandong Province, and the study period is 2000, 2005, 2010, 2015, and 2020, covering five time phases. The evaluation unit uses a 1 km × 1 km raster cell. The data includes land use data, digital elevation model data, normalized difference vegetation index (NDV) data, precipitation data, biomass data, habitat quality data, population density data, human activity footprint data, road data, and statistical data on grain yield and price. The land use data is used for calculating ecosystem service value and constructing a landscape ecological risk index, while the natural ecological factors and socioeconomic factors are used for analyzing the coupling coordination mechanism and threshold identification.

[0035] Specifically, in step S1, multi-source data is preprocessed. First, vector and raster data from different sources are converted to a unified format and projected onto the same spatial coordinate system. Second, data with different spatial resolutions are resampled, with a preferred resampling resolution of 1 km × 1 km. Subsequently, continuous driving factor data are normalized to eliminate dimensional differences and improve the comparability between different factors. Finally, a spatially consistent database of ecosystem services, landscape ecological risks, and driving factors is constructed using 1 km × 1 km raster units as the basic evaluation units.

[0036] Specifically, in step S2, an ecosystem service value model and a landscape ecological risk index model are constructed. The ecosystem service value is calculated using the equivalent factor method of ecosystem service value per unit area. Based on land use types such as cultivated land, forest land, grassland, water area, wetland, construction land, and unused land, the equivalent value per unit area is corrected by incorporating regional grain output value, thereby calculating the ecosystem service value of each evaluation unit. This method has been used in the calculation of five periods of data from Shandong Province from 2000 to 2020, and the overall changes in ecosystem service value are presented. Simultaneously, the landscape ecological risk index is constructed from landscape disturbance degree, landscape vulnerability degree, and landscape loss degree, used to characterize the comprehensive risk level of landscape pattern within the evaluation unit due to external disturbances and internal structural instability. By calculating the area proportion of different landscape types and their corresponding risk parameters within each evaluation unit, the spatial distribution results of the landscape ecological risk index for each grid unit are obtained.

[0037] Specifically, in step S3, the ecosystem service value and landscape ecological risk index are standardized, and a coupling coordination degree model is constructed. First, the coupling degree C is calculated, then the comprehensive development index T is calculated, and finally, the coupling coordination degree D is calculated to characterize the coordinated development level between ecosystem service supply capacity and landscape ecological risk. In this embodiment, ecosystem service value and landscape ecological risk have the same weight in the comprehensive development index. Based on the coupling coordination degree D value, the evaluation unit is divided into five categories: severely imbalanced development area, relatively imbalanced development area, barely balanced development area, moderately balanced development area, and highly balanced development area. Specifically: when 0 ≤ D < 0.35, it is a severely imbalanced development area; when 0.35 ≤ D < 0.5, it is a relatively imbalanced development area; when 0.5 ≤ D < 0.65, it is a barely balanced development area; when 0.65 ≤ D < 0.8, it is a moderately balanced development area; and when 0.8 ≤ D < 1, it is a highly balanced development area.

[0038] Specifically, in step S4, an XGBoost regression model is constructed with coupling coordination degree D as the dependent variable and natural ecological factors and human activity factors as independent variables. The natural ecological factors include at least DEM, precipitation, NDVI, biomass, and habitat quality; the human activity factors include at least construction density, population density, human activity footprint, distance from roads, and distance from water systems. All samples are divided into training and validation sets, preferably using a 70% training set and 30% validation set split, to train and validate the model, identifying key driving factors influencing changes in coupling coordination. The corresponding results in the paper show that the model achieves high fitting accuracy on both the training and validation sets, and can effectively reveal the nonlinear relationship between coupling coordination degree and multiple driving factors. Specifically, in step S5, the SHAP method is used to interpret the output of the XGBoost model.

[0039] By calculating the SHAP value of each driving factor, the contribution intensity and direction of influence of each factor on the coupling coordination degree are obtained. Furthermore, feature dependency curves are constructed, and the nonlinear response intervals and critical thresholds of key driving factors are identified by combining the curve slope changes. In this embodiment, the identified key driving factors mainly include construction density, population density, NDVI, and DEM, among which construction density has the most significant impact on coupling coordination degree. The threshold characteristics of each key factor are further identified: the construction density threshold is approximately 24.69%, and the population density threshold is approximately 1468 person-km². -2 The NDVI threshold is approximately 0.71, and the DEM threshold is approximately 110 m.

[0040] Specifically, in step S6, the coupling coordination level is superimposed with the threshold interval of key driving factors to construct a two-dimensional ecological security zoning rule. First, based on the critical thresholds of the key driving factors identified in step S5, a sensitive area discrimination interval is constructed. Preferably, the threshold sensitivity interval of the corresponding factor is formed by expanding forward and backward by 10% with the critical threshold of each driving factor as the center.

[0041] When at least one key driving factor of an evaluation unit falls within its corresponding sensitive interval, the evaluation unit is determined to be in a driving threshold sensitive state. Subsequently, the coupling coordination level and the driving threshold sensitive state are overlaid for analysis, with the specific partitioning rules as follows: (1) When the coupling coordination degree of the evaluation unit is in a state of severe developmental imbalance, it is classified as an ecologically fragile area; (2) When the evaluation unit coupling coordination is in a highly balanced development state and all key driving factors are not in the sensitive range, it is divided into a coordinated conservation area; (3) When the coupling coordination degree of the evaluation unit is not in a state of severe developmental imbalance, and at least one key driving factor is in the sensitive interval, it is classified as a sensitive area; Furthermore, if the dominant sensitive factor belongs to the natural ecological factor, it is classified as a vegetation-dominated sensitive area; if the dominant sensitive factor belongs to the socio-economic factor, it is classified as a socially dominated sensitive area. (4) The remaining evaluation units are divided into stable transition zones.

[0042] Finally, connectivity analysis and patch integration are performed on spatially adjacent grid units of the same type to eliminate isolated small patches and generate the final ecological security zoning map. This embodiment verifies that the present invention can achieve spatial identification of the coupling state of ecosystem services and landscape ecological risks, quantitative extraction of key driving factor thresholds, and sensitive constraint expression of zoning boundaries, thereby improving the scientific rigor, relevance, and operability of ecological security zoning results.

Claims

1. A method for coupled zoning of ecosystem services and landscape ecological risk based on XGBoost-SHAP-driven threshold identification, characterized in that, Includes the following steps: S1. Determine the study area, time period, and evaluation unit scale; collect multi-source data and perform projection unification, resampling, and standardization to construct a spatially consistent raster database. S2. Calculate the ecosystem service value of each evaluation unit, construct a landscape ecological risk index model, and obtain its spatial distribution results. S3. After standardizing both, construct a coupling coordination degree model and calculate the coupling coordination degree level of each evaluation unit. S4. Construct an XGBoost regression model with coupling coordination degree as the dependent variable and natural and human activity factors as independent variables to identify key driving factors. S5. Calculate the contribution value of driving factors using the SHAP method, construct feature dependencies, and identify their nonlinear response intervals and critical thresholds. S6. Combine the coupling coordination degree level and the driving factor threshold intervals to construct a two-dimensional partitioning rule, complete the ecological security partitioning of the study area, and integrate patches to obtain the final results.

2. The method according to claim 1, characterized in that, Step S1, multi-source data preprocessing, includes: converting vector and raster data into a unified raster format and uniform projection; resampling the data to a preferred 1 km × 1 km resolution; and applying extreme value normalization to continuous driving factors. The normalization formula is as follows: Ultimately, a database of ecosystem services and driving factors with grids as the evaluation unit will be constructed.

3. The method according to claim 1, characterized in that, The formula for calculating the value of ecosystem services in step S2 is as follows: The service value equivalent per unit area, adjusted for regional grain output value, is calculated using the following formula: The formula for calculating the landscape ecological risk index is: ,in , D k For landscape disturbance degree, V k Landscape vulnerability.

4. The method according to claim 1, characterized in that, In step S3, the ecosystem service value and landscape ecological risk index are first standardized to obtain U1 and U2, and then the coupling degree is calculated. Comprehensive Development Index T=0.5U 1 + 0.5U 2 The final coupling coordination degree is obtained. Based on the D value, five coordination levels are defined: D<0.35 is the severely imbalanced zone, 0.35≤D<0.5 is the relatively imbalanced zone, 0.5≤D<0.65 is the barely balanced zone, 0.65≤D<0.8 is the moderately balanced zone, and D≥0.8 is the highly balanced zone.

5. The method according to claim 1, characterized in that, The prediction function of the XGBoost regression model in step S4 is: The objective function is: The formula for controlling model complexity is: .

6. The method according to claim 1, characterized in that, The formula for calculating the SHAP value of feature i in step S5 is as follows: ; Constructing feature dependency functions When satisfied When the variable interval is equal to the driving threshold interval, the corresponding variable interval is equal to the driving threshold interval.

7. The method according to claim 1, characterized in that, The sensitive region discrimination interval for the driving factor in step S6 is: Construct a sensitivity discrimination function F k ( x k ),when x k ∈ I k hour F k ( x k If )=1, it is determined to be a sensitive state; The superimposed coupling coordination level and sensitivity status are used to divide the ecological zones into four categories: D<0.35 is the ecologically fragile zone, D≥0.8 and no sensitive factors is the coordinated conservation zone, 0.35≤D<0.65 and there are sensitive factors is the sensitive zone, which is further subdivided according to factor type, and the rest are stable transition zones. Finally, connectivity analysis and patch integration are performed on the same type of raster.

Images