A method for identifying spatial heterogeneity of ecosystem services based on driving factors

Abstract

The application discloses an ecosystem service spatial heterogeneity identification method based on driving factors, comprising the following steps: 1) fusing multi-source heterogeneous geographic spatial data to generate a basic quantitative data set; 2) constructing an ecosystem service driving factor database, and screening out dominant driving factors with significant explanatory power in combination with an improved random forest method; 3) based on the nonlinear response relationship of the dominant driving factors to the ecosystem service, a regional adaptive spatial error model is constructed to identify the differences of the dominant driving factors in the spatial heterogeneity of the ecosystem service; 4) the target region is divided into multiple heterogeneity control units by comprehensively considering the distribution of the dominant driving factors in the space, and corresponding optimization strategies are provided. The application solves the problem that the existing method lacks consideration of the comprehensive analysis of the spatial heterogeneity of the ecosystem service and driving factors, and is helpful to realize more accurate ecological protection and resource management, and provides a scientific basis for ecological environment protection decision-making.

Description

Technical Field

[0001] This invention relates to the field of ecological conservation technology, specifically a method for identifying spatial heterogeneity of ecosystem services based on driving factors. Background Technology

[0002] Ecosystem services, as the core functions of natural ecosystems in maintaining life, providing material supply, and regulating the environment, have spatial differentiation characteristics that directly affect the accuracy of ecological protection policies. Traditional methods mainly focus on the characteristics of the ecosystem itself or rely on a single evaluation indicator, depending on experience or rough regional delineation. The defined protected areas may have a certain degree of subjectivity, ignoring the spatial heterogeneity and spatiotemporal changes of ecosystem services. This makes the identification results lack adaptability in the face of complex and ever-changing geographical environments, unable to characterize the spatial differentiation characteristics of response intensity caused by geographical differences, and difficult to reveal the deep causes of the spatial pattern of ecosystem services. It is easy to cause a disconnect between the division of management units and the ecological process mechanism, resulting in low resource allocation efficiency and poor accuracy of the formulated ecological protection policies.

[0003] While existing studies have attempted to introduce a driving factor screening mechanism, they generally suffer from problems such as incomplete elimination of multicollinearity interference, coarse modeling of nonlinear response relationships, and lack the ability to quantitatively characterize the spatial nonstationary effects of driving factors. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies by providing a method for identifying spatial heterogeneity of ecosystem services based on driving factors. By combining multi-source spatial data with a quantitative assessment model of ecosystem services, a spatial distribution map of ecosystem services is constructed. The relationship between ecosystem services and driving factors is analyzed, and by comprehensively considering the spatial distribution of driving factors, the intrinsic mechanism of spatial heterogeneity of ecosystem services is systematically revealed. Heterogeneity control units are delineated, and targeted zoning optimization strategies are proposed. This method solves the problem that existing methods lack consideration for the spatial heterogeneity of ecosystem services and the comprehensive analysis of driving factors, contributing to more precise ecological protection and resource management and providing a scientific basis for ecological and environmental protection decision-making.

[0005] To achieve the above objectives, the present invention adopts the following technical solution.

[0006] A method for identifying spatial heterogeneity of ecosystem services based on driving factors includes the following steps: Step S1: Integrate multi-source heterogeneous geospatial data, including soil geological parameters, meteorological observation data, and multispectral vegetation indices, to generate a basic quantitative dataset that reflects the natural capital supply pattern of the target area. Step S2: Based on the basic quantitative dataset generated in Step S1, construct an ecosystem service driver factor database, and combine the improved random forest method to analyze the relationship between ecosystem services and driver factors, and screen out the dominant driver factors with significant explanatory power. The improved random forest method uses decision trees of different depths to construct a random forest model. The grid cells in the basic quantized dataset generated in step S1 are used as samples, and the constructed ecosystem service driving factor database is used as training data. The model is trained by assigning different weights to each sample based on the driving factors. After training, the final prediction model is obtained by weighted averaging. The feature importance is calculated by combining the weight of each tree in the prediction model and the contribution of each tree feature. This measures the driving ability of each driving factor in the database and selects the dominant driving factors with significant explanatory power. Step S3: For the selected dominant driving factors, based on the nonlinear response relationship of ecosystem services to the dominant driving factors, construct a regional adaptive spatial error model to obtain the spatial distribution of the dominant driving factors and identify the characteristics and differences of each dominant driving factor in the spatial heterogeneity of ecosystem services. Step S4: Based on the multi-driving factor spatiotemporal cluster partitioning method, obtain spatial clusters with the same dominant driving factor; The multi-driving factor spatiotemporal clustering method integrates single-year multi-driving factor clustering, multi-year clustering result hierarchical superposition and spatial contiguous correction method to identify and aggregate grid units with consistent driving characteristics, thereby obtaining spatial clusters with the same dominant driving factor. Step S5: Integrating the spatial differentiation characteristics of each dominant driving factor and the contiguous spatial clusters obtained by spatiotemporal clustering of multiple driving factors, the target region is divided into multiple heterogeneous control units by using spatial superposition and dominant factor priority assignment method, and targeted optimization strategies are provided.

[0007] Specifically, step S1 generates a basic quantitative dataset reflecting the natural capital supply pattern of the target region, including the following: Step S11: Acquisition and preprocessing of multi-source heterogeneous geospatial data; Soil geological parameters, meteorological observation data, and multispectral vegetation indices were collected, and the collected data were standardized in format, calibrated in coordinates, and outlier removed. First, data from different sources are converted to a unified format and standardized, with standardized field naming. Then, coordinate calibration is performed, converting all data into latitude and longitude coordinates. Finally, a time-weighted Gaussian mixture model (GMM) is constructed to calculate the probability density. To eliminate outliers, the probability density The calculation formula is as follows: ; ; In the above formula, yes The probability density at time t. For time indexing; For the first The weight of each component; The number of Gaussian components; It is a weighting factor that changes over time; express dimensional parameter vector, , for The One component; for 3D mean vector; for Covariance matrix; Let be the probability density function of a multivariate Gaussian normal distribution, which is read as given parameters and Down, dimensional parameter vector The probability of occurrence; Represents an exponential function; This is the time decay factor; Finally, the extremely low probability quantile is chosen as the threshold. If at a certain time t, If the threshold is less than 1, then a decision is made. Data at any given time that is outlier will be removed. The selection of extremely low probability quantiles as thresholds specifically involves taking the 1st to 5th percentile of the probability density values ​​of all samples as the threshold. This threshold range can effectively eliminate outliers caused by sensor errors, abnormal data entry, etc., while avoiding excessive deletion of valid samples located at the distribution boundary.

[0008] Step S12: Data integration and dataset generation; Data is integrated using spatial fusion technology, and a basic quantitative dataset reflecting the regional natural capital supply pattern is generated through quantitative calculations.

[0009] Specifically, step S2 identifies the dominant driving factors with significant explanatory power, including the following: Step S21: Integrate decision trees of different depths; A random forest model is constructed using multi-level decision trees, where each tree in the constructed random forest model has a different maximum depth; Step S22: Train and build a prediction model; Using the grid cells in the basic quantized dataset generated in step S1 as samples, the model is trained by assigning different weights to each sample based on the driving factors. After training, the final prediction model expression is obtained through a weighted average method: ; In the above formula, This is the predicted output for the entire tree; For the first The predicted output for each tree; For the first The weight of each tree; This represents the total number of decision trees in the random forest. The input feature vector corresponds to the values ​​of each driving factor; Combine the weight of each tree The weighted representation of feature importance is obtained by considering the contribution of each tree feature: ; In the above formula, Representation of features Weighted feature importance; Indicates the first Features of a tree The importance of; For the first The weight of each tree; Finally, based on the calculated feature importance, the driving power of each driving factor is measured, and a threshold is selected or the dominant driving factors with significant explanatory power are screened out by sorting according to their importance.

[0010] Specifically, in step S3, a regional adaptive spatial error model is constructed to obtain the distribution of the dominant driving factors in each space. The process is as follows: Based on the traditional spatial error model, and considering the characteristics of the correlation strength and influence range of spatial errors varying with region, a regional adaptive spatial error model is constructed. The model's formula is expressed as follows: ; ; In the above formula, Serving the ecosystem in ecological areas The vector of observed values; For the first A matrix of driving factors; For the first A vector of driving factor coefficients for each ecological region; The spatial error term is divided into... Subset ,correspond One ecological region; For the first The spatial weight matrix of each ecological region is constructed based on geographical proximity or ecological similarity; For the first Spatial autocorrelation parameters of an ecological region; The error term is an independent error term that follows a normal distribution. , For the error variance, It is the identity matrix; Will Plotting this on a map reveals the spatial differences in ecosystem service-driven factors. Furthermore, the aforementioned The map contains the following details: First, the maximum likelihood estimation method is used to solve the parameters of the regional adaptive spatial error model: the basic quantization dataset generated in step S1 is grouped according to ecological regions, with ecosystem services in ecological regions. Observation vector As the dependent variable, the driving factor matrix As the independent variable, combined with the spatial weight matrix The regression coefficient vector is estimated simultaneously through an iterative optimization algorithm. Spatial autocorrelation coefficient and error variance To maximize the likelihood function, we obtain the first... Driving factor coefficient vector of each ecological region ; Then the resulting driving factor coefficient vector As an attribute field, it is associated with the corresponding ecological region in ArcGIS software. In the spatial surface feature attribute table; for each driving factor, its location in the ecological region is selected. The corresponding coefficient components (such as the first) The coefficient components of each driving factor Ecological regions are classified using the natural discontinuity grading method. The coefficient estimates are divided into 5 to 7 levels and rendered with gradient colors to output a spatial distribution map of the driving effect of the driving factor on ecosystem services. A positive coefficient estimate indicates that the driving factor has a positive promoting effect on ecosystem services, while a negative coefficient estimate indicates a negative inhibiting effect. The larger the absolute value of the coefficient, the higher the driving intensity of the driving factor on ecosystem services. By analogy, the spatial distribution maps of all dominant driving factors are completed one by one. By comparing and analyzing the spatial differentiation patterns of each driving factor, the characteristics and differences of each dominant driving factor in the spatial heterogeneity of ecosystem services are identified.

[0011] Traditional spatial error models typically output regression results in numerical tables, which are difficult to intuitively reflect the spatial differentiation patterns of the influence intensity of driving factors. This invention spatializes the coefficient vectors of driving factors in each ecological region, enabling a direct visual representation of the positive and negative directions, intensity levels, and spatial distribution patterns of the impact of driving factors on ecosystem services. This provides an intuitive spatial basis for subsequent steps S4, which involves spatiotemporal clustering based on the similarity of driving factors, and S5, which involves delineating heterogeneous control units.

[0012] Specifically, in step S4, based on the multi-driving-factor spatiotemporal cluster partitioning method, spatial clusters with the same dominant driving factor are obtained. The specific process is as follows: Step S41: Single-year multi-driving factor clustering; Based on the multi-driving factor ensemble data of the region adaptive spatial error model, the comprehensive distance of each raster to other rasters is calculated: ; In the above formula, Indicates the first The grid to the first The combined distance of each grid cell; The total number of driving factors, Indicates the driving factor index; This is the value of the s-th raster after normalization to the q-th driving factor; This is the value of the c-th raster after normalization to the q-th driving factor; ; ; In the above formula, Let be the initial value of the s-th raster in the q-th driving factor; This is the initial value of the c-th raster in the q-th driving factor; The function is for finding the minimum value; This is a function to find the maximum value. Let q be the set of the q-th driving factor across all grids; Grouping factors with a combined distance less than a preset threshold into the same cluster, we finally obtain the driving factor combination clusters for each year; ; In the above formula, F is the set of multiple driving factors for each year; For the first The set of multiple driving factors for a year; m is the total number of years; Step S42: Layered overlay of multi-year clustering results; Based on all the single-year combination clusters in step S41, the most frequently occurring driving combination clusters of each grid over many years are statistically analyzed to determine the long-term dominant driving factor combination. ; In the above formula, The combination of driving factors corresponding to the cluster that appears most frequently; For the first The number of times class-driven factor combinations occur in all years; The total number of combinations of driving factors; Step S43: Spatial alignment correction; Based on the multi-year clustering overlay results of step S42, further spatial clustering corrections are made to ensure spatial cluster contiguousness, ultimately obtaining the spatial clusters of the dominant driving factors. ; In the above formula, For the first Spatial clusters of dominant driving factors after grid correction; For the first The long-term dominant driving factor combination of each grid; For window Within, the long-term dominant driving factor combination is as follows The number of grid cells; For window The most frequently occurring combination of driving factors, For the set of all driving factor combination categories; Therefore, the first A 3×3 neighborhood window centered on a grid.

[0013] Specifically, in step S5, the target region is divided into multiple heterogeneous control units using spatial overlay and dominant factor priority assignment methods, and corresponding optimization strategies are provided, as follows: Step S51: Based on the spatial distribution of the dominant driving factors and the spatial cluster boundaries, perform spatial superposition to clarify the type, intensity and combination pattern of the dominant driving factors within each spatial cluster; Step S52: Assign priority weights to different types of dominant driving factors based on the ecological management objectives of the target area; Step S53: Based on the criteria of consistent dominant driving factor type, intensity difference ≤10%, and continuous and unbroken spatial boundary within the same control unit, the target area is divided into multiple heterogeneous control units. Targeted regulation and optimization strategies are proposed for the driving factor characteristics and ecological problems of each unit. The targeted regulation and optimization strategies include: For areas dominated by human activity, the optimization strategy focuses on controlling the ecological impact of human activities, including improving land use patterns, enhancing urban greening, and reducing pollution emissions; in agricultural areas, precision agriculture technologies are adopted to reduce over-cultivation and irrational irrigation; in urbanized areas, the focus is on urban ecological restoration, including adding green roofs and parks. For regions dominated by ecological characteristics, the optimization strategy focuses on ecological protection and restoration, including strengthening the management of nature reserves, returning farmland to forest and restoring wetlands, and protecting biodiversity; in forest and grassland ecological areas, optimization measures include the construction of ecological corridors and the protection of species habitats. For regions dominated by meteorological factors, the optimization strategy focuses on climate adaptation and disaster prevention, including implementing water resource management, strengthening drought-resistant planting, and improving soil water retention capacity in arid areas; and strengthening wind disaster protection and flood control in storm-prone areas.

[0014] At the same time, the sustainability of ecosystem services in response to climate change can be ensured by setting up climate forecasting and emergency response mechanisms.

[0015] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention takes the core driving characteristics of ecosystem services as its starting point, uses an improved random forest method to screen the dominant driving factors, and combines a regional adaptive spatial error model to quantify the spatial heterogeneity response characteristics of each driving factor. This can reveal the differentiation patterns of ecosystem services in different regions and at different times, fully present their spatiotemporal change mechanisms, and provide a scientific basis for clarifying differentiated ecological protection strategies.

[0016] 2. This invention is based on a multi-driving factor spatiotemporal clustering method, which integrates single-year clustering, multi-year hierarchical superposition, and spatial contiguous correction to identify contiguous spatial clusters with consistent driving characteristics. On this basis, the target area is divided into heterogeneous control units by combining the spatial distribution and priority weight of the dominant driving factors. This avoids the problems of traditional administrative divisions being disconnected from ecological processes, relying on experience for division, and being highly subjective, making the definition of control units more ecologically reasonable.

[0017] 3. The zoning scheme based on driving factors in this invention starts from the root of the formation and evolution of ecosystem services, ensuring the homogeneity of ecological processes and service supply within the unit. Case analysis shows that the standard deviation of the trade-offs and synergies of ecosystem services under the driving factor zoning is significantly smaller than that under the administrative division, which verifies the superiority of this zoning scheme in matching the ecological background characteristics and reducing internal variation.

[0018] 4. This invention proposes differentiated regulation and optimization strategies for different types of dominant driving factors, promoting the transformation of ecological protection policies from "one-size-fits-all" to "precision and refinement," which helps improve resource allocation efficiency and enhance the pertinence and effectiveness of ecological restoration and management. Attached Figure Description

[0019] Figure 1 This is a flowchart of a method for identifying spatial heterogeneity of ecosystem services based on driving factors, according to the present invention. Detailed Implementation

[0020] To facilitate understanding and implementation of the present invention by those skilled in the art, the various steps of the method proposed in this invention are described in detail below. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. Furthermore, it should be understood that after reading the teachings of this invention, those skilled in the art can make various modifications or alterations to the invention, and these equivalent forms also fall within the scope defined by the appended claims.

[0021] Example like Figure 1 As shown, this invention discloses a method for identifying spatial heterogeneity of ecosystem services based on driving factors, comprising the following steps: Step S1: Integrate multi-source heterogeneous geospatial data, including soil geological parameters, meteorological observation data, and multispectral vegetation indices, to generate a basic quantitative dataset that reflects the natural capital supply pattern of the target area. In this embodiment, the target area is the Yangtze River Economic Belt, and the research period is from 2000 to 2020. Since the Yangtze River Economic Belt has multiple needs such as water source protection, ecological security and carbon emission reduction, water conservation, soil conservation and carbon sink are selected as core ecosystem services for quantitative calculation. Step S2: Based on the basic quantitative dataset generated in Step S1, construct an ecosystem service driver factor database, and combine the improved random forest method to analyze the relationship between ecosystem services and driver factors, and screen out the dominant driver factors with significant explanatory power. The improved random forest method uses decision trees of different depths to construct a random forest model. The grid cells in the basic quantized dataset generated in step S1 are used as samples, and the constructed ecosystem service driving factor database is used as training data. The model is trained by assigning different weights to each sample based on the driving factors. After training, the final prediction model is obtained by weighted averaging. The feature importance is calculated by combining the weight of each tree in the prediction model and the contribution of each tree feature. This measures the driving ability of each driving factor in the database and selects the dominant driving factors with significant explanatory power. Step S3: For the selected dominant driving factors, based on the nonlinear response relationship of ecosystem services to the dominant driving factors, construct a regional adaptive spatial error model to obtain the spatial distribution of the dominant driving factors and identify the characteristics and differences of each dominant driving factor in the spatial heterogeneity of ecosystem services. Step S4: Based on the multi-driving factor spatiotemporal cluster partitioning method, obtain spatial clusters with the same dominant driving factor; The multi-driving factor spatiotemporal clustering method integrates single-year multi-driving factor clustering, multi-year clustering result hierarchical superposition and spatial contiguous correction method to identify and aggregate grid units with consistent driving characteristics, thereby obtaining spatial clusters with the same dominant driving factor. Step S5: Integrating the spatial differentiation characteristics of each dominant driving factor and the contiguous spatial clusters obtained by spatiotemporal clustering of multiple driving factors, the target region is divided into multiple heterogeneous control units by using spatial superposition and dominant factor priority assignment method, and targeted optimization strategies are provided.

[0022] Specifically, step S1 generates a basic quantitative dataset reflecting the natural capital supply pattern of the target region, including the following: Step S11: Acquisition and preprocessing of multi-source heterogeneous geospatial data; Soil geological parameters, meteorological observation data, and multispectral vegetation indices were collected, and the collected data were standardized in format, calibrated in coordinates, and outlier removed. First, data from different sources are converted to a unified format and standardized, with standardized field naming. Then, coordinate calibration is performed, converting all data into latitude and longitude coordinates. Finally, a time-weighted Gaussian mixture model (GMM) is constructed to calculate the probability density. To eliminate outliers, the probability density The calculation formula is as follows: ; ; In the above formula, yes The probability density at time t. For time indexing; For the first The weight of each component; The number of Gaussian components; It is a weighting factor that changes over time; express dimensional parameter vector, , for The One component; for 3D mean vector; for Covariance matrix; Let be the probability density function of a multivariate Gaussian normal distribution, which is read as given parameters and Down, dimensional parameter vector The probability of occurrence; Represents an exponential function; This is the time decay factor; Finally, the extremely low probability quantile is chosen as the threshold. If at a certain time t, If the threshold is less than 1, then a decision is made. Data at any given time that is outlier will be removed. In this embodiment, the extremely low probability quantile is selected by taking the 1% quantile of the probability density values ​​of all samples as the threshold. This threshold can effectively remove outliers caused by sensor errors, abnormal data entry, etc., and can also avoid excessive deletion of valid samples located at the distribution boundary.

[0023] Step S12: Data integration and dataset generation; Data is integrated using spatial fusion technology, and a basic quantitative dataset reflecting the regional natural capital supply pattern is generated through quantitative calculations.

[0024] Specifically, step S2 involves screening out the dominant driving factors with significant explanatory power, including: Step S21: Integrate decision trees of different depths; A random forest model is constructed using multi-level decision trees, where each tree in the constructed random forest model has a different maximum depth; Step S22: Train and build a prediction model; Using the grid cells in the basic quantized dataset generated in step S1 as samples, the model is trained by assigning different weights to each sample based on the driving factors. After training, the final prediction model expression is obtained through a weighted average method: ; In the above formula, This is the predicted output for the entire tree; For the first The predicted output for each tree; For the first The weight of each tree; This represents the total number of decision trees in the random forest. The input feature vector corresponds to the values ​​of each driving factor; Combine the weight of each tree The weighted representation of feature importance is obtained by considering the contribution of each tree feature: ; In the above formula, Representation of features Weighted feature importance; Indicates the first Features of trees The importance of; For the first The weight of each tree.

[0025] Finally, based on the calculated feature importance, the driving power of each driving factor is measured, and a threshold is selected or the dominant driving factors with significant explanatory power are screened out by sorting according to their importance.

[0026] Specifically, in step S3, a regional adaptive spatial error model is constructed to obtain the distribution of the dominant driving factors in each space, including the following: Based on the traditional spatial error model, and considering the characteristics of the correlation strength and influence range of spatial errors varying with region, a regional adaptive spatial error model is constructed. The model's formula is expressed as follows: ; ; In the above formula, Serving the ecosystem in ecological areas The vector of observed values; For the first A matrix of driving factors; For the first A vector of driving factor coefficients for each ecological region; The spatial error term is divided into... Subset ,correspond One ecological region; For the first The spatial weight matrix of each ecological region is constructed based on geographical proximity or ecological similarity; For the first Spatial autocorrelation parameters of an ecological region; The error term is an independent error term that follows a normal distribution. , For the error variance, It is the identity matrix; Using ArcGIS software By plotting on a map, we can obtain a spatial distribution map of the differences in the impact of different driving factors on ecosystem services.

[0027] In this embodiment, the step of... The map contains the following details: First, the maximum likelihood estimation method is used to solve the parameters of the regional adaptive spatial error model: the basic quantization dataset generated in step S1 is grouped according to ecological regions, with ecosystem services in ecological regions. Observation vector As the dependent variable, the driving factor matrix As the independent variable, combined with the spatial weight matrix The regression coefficient vector is estimated simultaneously through an iterative optimization algorithm. Spatial autocorrelation coefficient and error variance To maximize the likelihood function, we obtain the first... Driving factor coefficient vector of each ecological region ; Then the resulting driving factor coefficient vector As an attribute field, it is associated with the corresponding ecological region in ArcGIS software. In the spatial surface feature attribute table; for each driving factor, its location in the ecological region is selected. The corresponding coefficient components (such as the first) The coefficient components of each driving factor Ecological regions are classified using the natural discontinuity grading method. The coefficient estimates are divided into 5 to 7 levels and rendered with gradient colors to output a spatial distribution map of the driving effect of the driving factor on ecosystem services. A positive coefficient estimate indicates that the driving factor has a positive promoting effect on ecosystem services, while a negative coefficient estimate indicates a negative inhibiting effect. The larger the absolute value of the coefficient, the higher the driving intensity of the driving factor on ecosystem services. By analogy, the spatial distribution maps of all dominant driving factors are completed one by one. By comparing and analyzing the spatial differentiation patterns of each driving factor, the characteristics and differences of each dominant driving factor in the spatial heterogeneity of ecosystem services are identified.

[0028] Specifically, in step S4, based on the multi-driving-factor spatiotemporal cluster partitioning method, spatial clusters with the same dominant driving factor are obtained. The specific process is as follows: Step S41: Single-year multi-driving factor clustering; Based on the multi-driving factor set data of the regional adaptive spatial error model, the comprehensive distance of each grid cell to other grid cells is calculated. ; In the above formula, Indicates the first The grid to the first The combined distance of each grid cell; The total number of driving factors, Indicates the driving factor index; This is the value of the s-th raster after normalization to the q-th driving factor; This is the value of the c-th raster after normalization to the q-th driving factor; ; ; In the above formula, Let be the initial value of the s-th raster in the q-th driving factor; This is the initial value of the c-th raster in the q-th driving factor; The function is for finding the minimum value; This is a function to find the maximum value. Let q be the set of the q-th driving factor across all grids; Grouping factors with a combined distance less than a preset threshold into the same cluster, we finally obtain the driving factor combination clusters for each year; ; In the above formula, F is the set of multiple driving factors for each year; For the first The set of multiple driving factors for a year; m is the total number of years; Step S42: Layered overlay of multi-year clustering results; Based on all the single-year combination clusters in step S41, the most frequently occurring driving combination clusters of each grid over many years are statistically analyzed to determine the long-term dominant driving factor combination. ; In the above formula, The combination of driving factors corresponding to the cluster that appears most frequently; For the first The number of times class-driven factor combinations occur in all years; The total number of driving factor combinations; Step S43: Spatial alignment correction; Based on the multi-year clustering overlay results of step S42, further spatial clustering corrections are made to ensure spatial cluster contiguousness, ultimately obtaining the spatial clusters of the dominant driving factors. ; In the above formula, For the first Spatial clusters of dominant driving factors after grid correction; For the first The long-term dominant driving factor combination of each grid; For window Within, the long-term dominant driving factor combination is as follows The number of grid cells; For window The most frequently occurring combination of driving factors For the set of all driving factor combination categories; Therefore, the first A 3×3 neighborhood window centered on a grid.

[0029] Specifically, in step S5, the target region is divided into multiple heterogeneous control units using spatial overlay and dominant factor priority assignment methods, and corresponding optimization strategies are provided, as follows: Step S51: Based on the spatial distribution of the dominant driving factors and the spatial cluster boundaries, perform spatial superposition to clarify the type, intensity and combination pattern of the dominant driving factors within each spatial cluster; Step S52: Assign priority weights to different types of dominant driving factors based on the ecological management objectives of the target area; Step S53: Based on the criteria of consistent dominant driving factor type, intensity difference ≤10%, and continuous and unbroken spatial boundary within the same control unit, the target area is divided into multiple heterogeneous control units. Targeted regulation and optimization strategies are proposed for the driving factor characteristics and ecological problems of each unit. The targeted regulation and optimization strategies include: For areas dominated by human activity, the optimization strategy focuses on controlling the ecological impact of human activities, including improving land use patterns, enhancing urban greening, and reducing pollution emissions; in agricultural areas, precision agriculture technologies are adopted to reduce over-cultivation and irrational irrigation; in urbanized areas, the focus is on urban ecological restoration, including adding green roofs and parks. For regions dominated by ecological characteristics, the optimization strategy focuses on ecological protection and restoration, including strengthening the management of nature reserves, returning farmland to forest and restoring wetlands, and protecting biodiversity; in forest and grassland ecological areas, optimization measures include the construction of ecological corridors and the protection of species habitats. For regions dominated by meteorological factors, the optimization strategy focuses on climate adaptation and disaster prevention, including implementing water resource management, strengthening drought-resistant planting, and improving soil water retention capacity in arid areas; and strengthening wind disaster protection and flood control in storm-prone areas.

[0030] The technical effects of the method of the present invention will be further illustrated below through an example of spatial heterogeneity identification in the Yangtze River Economic Belt region.

[0031] This study focuses on the Yangtze River Economic Belt as the target region, covering the period from 2000 to 2020. Data on three types of ecosystem services—water conservation, soil conservation, and carbon sequestration—were collected within the target region. ArcGIS software was used for preprocessing the multi-source data, including format standardization, coordinate calibration, and outlier removal. RStudio software was used for statistical analysis. Origin software was used for visualization, resulting in the calculated results of ecosystem service trade-offs and synergies under different zoning schemes, as shown in Table 1 below.

[0032] Table 1. Calculation results of ecosystem service trade-offs and synergies under different zoning schemes ; As shown in Table 1, in the administrative zoning scheme, the mean values ​​of ecosystem services for water conservation and carbon sink, carbon sink and soil conservation, and water conservation and soil conservation are 0.52, 0.53, and 0.63, respectively; while in the driving factor zoning scheme, the corresponding mean values ​​are 0.41, 0.40, and 0.48, respectively. Further analysis using standard deviation to quantify the degree of variation in the trade-offs and synergies of ecosystem services within each group under each zoning scheme shows that the trade-offs and synergies of ecosystem services within each group are more consistent under the driving factor zoning scheme, and the within-group standard deviation is significantly smaller than that of the administrative zoning scheme.

[0033] The above results indicate that the method of this invention, based on the core influencing factors of ecosystem service formation and evolution (such as climate, topography, land use, and intensity of human activities), ensures the homogeneity of ecological processes and service supply within the region from the root, thus resulting in smaller fluctuations in intra-group trade-offs and synergies. In contrast, the traditional administrative zoning method uses administrative boundaries as the dividing standard, without considering the continuity and intrinsic connections of ecosystems, which can easily lead to significant differences in ecological background within the same administrative unit, resulting in higher intra-group standard deviations.

[0034] In summary, the driving factor zoning scheme provided by this invention has significant advantages in matching the background characteristics of ecosystems and reducing internal heterogeneity. The low standard deviation within groups directly confirms the rationality and effectiveness of the zoning scheme, indicating that the method of this invention can more scientifically reveal the spatial differentiation patterns of ecosystem services.

[0035] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for identifying spatial heterogeneity of ecosystem services based on driving factors, characterized in that, Includes the following steps: Step S1: Integrate multi-source heterogeneous geospatial data, including soil geological parameters, meteorological observation data, and multispectral vegetation indices, to generate a basic quantitative dataset that reflects the natural capital supply pattern of the target area. Step S2: Based on the basic quantitative dataset generated in Step S1, construct an ecosystem service driver factor database, and combine the improved random forest method to analyze the relationship between ecosystem services and driver factors, and screen out the dominant driver factors with significant explanatory power. The improved random forest method uses decision trees of different depths to construct a random forest model. The grid cells in the basic quantized dataset generated in step S1 are used as samples, and the constructed ecosystem service driving factor database is used as training data. The model is trained by assigning different weights to each sample based on the driving factors. After training, the final prediction model is obtained by weighted averaging. The feature importance is calculated by combining the weight of each tree in the prediction model and the contribution of each tree feature. This measures the driving ability of each driving factor in the database and selects the dominant driving factors with significant explanatory power. Step S3: For the selected dominant driving factors, based on the nonlinear response relationship of ecosystem services to the dominant driving factors, construct a regional adaptive spatial error model to obtain the spatial distribution of the dominant driving factors and identify the characteristics and differences of each dominant driving factor in the spatial heterogeneity of ecosystem services. Step S4: Based on the multi-driving factor spatiotemporal cluster partitioning method, obtain spatial clusters with the same dominant driving factor; The multi-driving factor spatiotemporal clustering method integrates single-year multi-driving factor clustering, multi-year clustering result hierarchical superposition and spatial contiguous correction method to identify and aggregate grid units with consistent driving characteristics, thereby obtaining spatial clusters with the same dominant driving factor. Step S5: Integrating the spatial differentiation characteristics of each dominant driving factor and the contiguous spatial clusters obtained by spatiotemporal clustering of multiple driving factors, the target region is divided into multiple heterogeneous control units by using spatial superposition and dominant factor priority assignment method, and targeted optimization strategies are provided.

2. The method for identifying spatial heterogeneity of ecosystem services based on driving factors according to claim 1, characterized in that, Step S1 generates a basic quantitative dataset reflecting the natural capital supply pattern of the target region, which includes: Step S11: Acquisition and preprocessing of multi-source heterogeneous geospatial data; Soil geological parameters, meteorological observation data, and multispectral vegetation indices were collected, and the collected data were standardized in format, calibrated in coordinates, and outlier removed. First, data from different sources are converted to a unified format and standardized, with standardized field naming. Then, coordinate calibration is performed, converting all data into latitude and longitude coordinates. Finally, a time-weighted Gaussian mixture model (GMM) is constructed to calculate the probability density. To eliminate outliers, the probability density The calculation formula is as follows: ； ； In the above formula, yes The probability density at time t. For time indexing; For the first The weight of each component; The number of Gaussian components; It is a weighting factor that changes over time; express dimensional parameter vector, , for The One component; for 3D mean vector; for Covariance matrix; Let be the probability density function of a multivariate Gaussian normal distribution, which is read as given parameters and Down, dimensional parameter vector The probability of occurrence; Represents an exponential function; This is the time decay factor; Finally, the extremely low probability quantile is chosen as the threshold. If at a certain time t, If the threshold is less than 1, then a decision is made. Data at any given time that is an outlier will be removed. Step S12: Data integration and dataset generation; Data is integrated using spatial fusion technology, and a basic quantitative dataset reflecting the regional natural capital supply pattern is generated through quantitative calculations.

3. The method for identifying spatial heterogeneity of ecosystem services based on driving factors according to claim 2, characterized in that, Step S2 involves screening out the dominant driving factors with significant explanatory power, specifically including: Step S21: Integrate decision trees of different depths; A random forest model is constructed using multi-level decision trees, where each tree in the constructed random forest model has a different maximum depth; Step S22: Train and build a prediction model; Using the grid cells in the basic quantized dataset generated in step S1 as samples, the model is trained by assigning different weights to each sample based on the driving factors. After training, the final prediction model expression is obtained through a weighted average method: ； In the above formula, This is the predicted output for the entire tree; For the first The predicted output for each tree; For the first The weight of each tree; This represents the total number of decision trees in the random forest. The input feature vector corresponds to the values ​​of each driving factor; Combine the weight of each tree The weighted representation of feature importance is obtained by considering the contribution of each tree feature: ； In the above formula, Representation of features Weighted feature importance; Indicates the first Features of trees The importance of; For the first The weight of each tree; Finally, based on the calculated feature importance, the driving power of each driving factor is measured, and a threshold is selected or the dominant driving factors with significant explanatory power are screened out by sorting according to their importance.

4. The method for identifying spatial heterogeneity of ecosystem services based on driving factors according to claim 3, characterized in that, Step S3 involves constructing a regional adaptive spatial error model to obtain the spatial distribution of the dominant driving factors. Specific details include: Based on the traditional spatial error model, and considering the characteristics of the correlation strength and influence range of spatial errors varying with region, a regional adaptive spatial error model is constructed. The model's formula is expressed as follows: ； ； In the above formula, Serving the ecosystem in ecological areas The vector of observed values; For the first A matrix of driving factors; For the first A vector of driving factor coefficients for each ecological region; The spatial error term is divided into... Subset ,correspond One ecological region; For the first The spatial weight matrix of each ecological region is constructed based on geographical proximity or ecological similarity; For the first Spatial autocorrelation parameters of an ecological region; The error term is an independent error term that follows a normal distribution. , For the error variance, It is the identity matrix; Will Plotting this on a map reveals the spatial differences in ecosystem service-driven factors.

5. The method for identifying spatial heterogeneity of ecosystem services based on driving factors according to claim 1, characterized in that, In step S4, the spatial clusters with the same dominant driving factor are obtained based on the multi-driving-factor spatiotemporal cluster partitioning method. The specific process is as follows: Step S41: Single-year multi-driving factor clustering; Based on the multi-driving factor ensemble data of the region adaptive spatial error model, the comprehensive distance of each raster to other rasters is calculated: ； In the above formula, Indicates the first The grid to the first The combined distance of each grid cell; The total number of driving factors, Indicates the driving factor index; This is the value of the s-th raster after normalization to the q-th driving factor; This is the value of the c-th raster after normalization to the q-th driving factor; ； ； In the above formula, Let be the initial value of the s-th raster in the q-th driving factor; This is the initial value of the c-th raster in the q-th driving factor; The function is for finding the minimum value; This is a function to find the maximum value. Let q be the set of the q-th driving factor across all grids; Grouping factors with a combined distance less than a preset threshold into the same cluster, we finally obtain the driving factor combination clusters for each year; ； In the above formula, F is the set of multiple driving factors for each year; For the first A set of multiple driving factors for the year; m is the total number of years; Step S42: Layered overlay of multi-year clustering results; Based on all the single-year combination clusters in step S41, the most frequently occurring driving combination clusters of each grid over many years are statistically analyzed to determine the long-term dominant driving factor combination. ； In the above formula, The combination of driving factors corresponding to the cluster that appears most frequently; For the first The number of times class-driven factor combinations occur in all years; The total number of driving factor combinations; Step S43: Spatial alignment correction; Based on the multi-year clustering overlay results of step S42, further spatial clustering corrections are made to ensure spatial cluster contiguousness, ultimately obtaining the spatial clusters of the dominant driving factors. ； In the above formula, For the first Spatial clusters of dominant driving factors after grid correction; For the first The long-term dominant driving factor combination of each grid; For window Within, the long-term dominant driving factor combination is as follows The number of grid cells; For window The most frequently occurring combination of driving factors, For the set of all driving factor combination categories; Therefore, the first A 3×3 neighborhood window centered on a grid.

6. The method for identifying spatial heterogeneity of ecosystem services based on driving factors according to claim 1, characterized in that, In step S5, the target region is divided into multiple heterogeneous control units using spatial overlay and dominant factor priority assignment methods, and corresponding optimization strategies are provided, as detailed below: Step S51: Based on the spatial distribution of the dominant driving factors and the spatial cluster boundaries, perform spatial superposition to clarify the type, intensity and combination pattern of the dominant driving factors within each spatial cluster; Step S52: Assign priority weights to different types of dominant driving factors based on the ecological management objectives of the target area; Step S53: Based on the criteria of consistent dominant driving factor type, intensity difference ≤10%, and continuous and unbroken spatial boundary within the same control unit, the target area is divided into multiple heterogeneous control units. Targeted regulation and optimization strategies are proposed for the driving factor characteristics and ecological problems of each unit. The targeted regulation and optimization strategies include: For areas dominated by human activity, the optimization strategy focuses on controlling the ecological impact of human activities, including improving land use patterns, enhancing urban greening, and reducing pollution emissions; in agricultural areas, precision agriculture technologies are adopted to reduce over-cultivation and irrational irrigation; in urbanized areas, the focus is on urban ecological restoration, including adding green roofs and parks. For regions dominated by ecological characteristics, the optimization strategy focuses on ecological protection and restoration, including strengthening the management of nature reserves, returning farmland to forest and restoring wetlands, and protecting biodiversity; in forest and grassland ecological areas, optimization measures include the construction of ecological corridors and the protection of species habitats. For regions dominated by meteorological factors, the optimization strategy focuses on climate adaptation and disaster prevention, including implementing water resource management, strengthening drought-resistant planting, and improving soil water retention capacity in arid areas; and strengthening wind disaster protection and flood control in storm-prone areas.

Images