A multi-scale and multi-scenario method for predicting and evaluating ecosystem service trade-offs and synergies

Abstract

This application relates to a multi-scale and multi-scenario ecosystem service trade-off-synergy prediction and assessment method. This method couples the SD, PLUS and InVEST models, comprehensively considers the impact of future climate change, economic development, population and policy implementation on ecosystem services, minimizes human subjectivity in the prediction process, and can more accurately predict multi-scenario ecosystem service trade-offs / synergies.

Description

Technical Field

[0001] The present application belongs to the technical field of ecosystem service assessment, and specifically relates to a multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and assessment method. Background Art

[0002] Ecosystem services are the benefits that ecosystems can provide to human life and production, covering all aspects, including primary categories such as provisioning services, regulating services, cultural services and supporting services. The assessment of ecosystem services can not only assess the quality of food, water and timber, but also reflect the service capacity of carbon storage, biodiversity and soil conservation. In the existing technology, the assessment of ecosystem services is mostly aimed at overall regional restoration, but lacks a regional ecological protection perspective based on the combination of current ecological services and future comprehensive ecological services. In addition, the current prediction of the trade-off-synergy relationship of ecosystem services mainly relies on the prediction of land use. It fails to provide a high-precision comprehensive ecosystem service assessment method suitable for the future multi-scenario climate framework, and does not consider the impact of changes in social factors on the future evolution of ecosystem services. Summary of the invention

[0003] The present invention provides a multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method, which can realize the evaluation and prediction of future ecosystem services and fully consider the impact of changes in social factors on ecosystem services.

[0004] The technical solution of the present invention is as follows:

[0005] 1. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method, characterized by comprising the following steps:

[0006] S1) Select administrative units, river basins or urban agglomerations as target areas for assessing ecosystem service trade-offs-synergies;

[0007] S2) Soil conservation services, water conservation services, habitat quality services, carbon storage services and food supply services are used as key evaluation indicators for assessing the target area;

[0008] S3) Selecting key timing nodes; the key timing nodes include historical timing nodes and future timing nodes;

[0009] S4) Obtain multi-source remote sensing data and social statistical data in the target area, input the multi-source remote sensing data and social statistical data into the InVEST model, and obtain the spatial differentiation characteristics of ecosystem service evaluation indicators and their key time series nodes;

[0010] S5) Obtain the historical change patterns of key evaluation indicators in the target area based on the local spatial autocorrelation index, land use products and landscape pattern index based on the common edge adjacency spatial weight;

[0011] S6) Establishing development scenarios, including natural development scenarios, ecological protection scenarios, and economic development scenarios;

[0012] The natural development scenario is based on the average annual temperature, average annual precipitation, population and GDP data under the future climate framework;

[0013] The ecological protection scenario is based on the ecological land area, annual average temperature, annual average precipitation, population and GDP data under the future climate framework;

[0014] The economic development scenario is based on the cultivated land red line, the annual average temperature, annual average precipitation, population and GDP data under the future climate framework;

[0015] S7) constructing a system dynamics model SD, and predicting the land use demand area of ​​the target area under different development scenarios according to the system dynamics model;

[0016] S8) constructing a PLUS model based on the land use demand area under each development scenario obtained in step S7), and coupling the SD-PLUS model to predict the land use spatial pattern of the target area under each development scenario at future time series nodes;

[0017] The PLUS model is established based on a land use simulation model of a cellular automaton model and a random patch generation mechanism;

[0018] S9) Based on different climate frameworks, the key ecological parameters for future ecosystem service prediction are revised, and the optimal InVEST model for prediction is obtained by combining the spatial pattern of land use under each development scenario obtained in step S8); the trade-off-synergy relationship between different scenarios at different scales in the region is obtained with the help of the optimal InVEST model and analysis method.

[0019] 2. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 1, characterized in that:

[0020] In step S2), the soil conservation service is calculated as follows:

[0021] ;

[0022] ;

[0023] ;

[0024] In the formula, For Grid soil conservation capacity; is the total potential soil erosion of land without management measures; is the actual amount of soil erosion after taking soil conservation measures; For Grid Sediment retention capacity; is the rainfall erosivity factor; is the soil erodibility factor; is the slope length factor; and represent soil conservation measures factor and vegetation cover management factor respectively;

[0025] The calculation formula for water conservation services is as follows:

[0026]

[0027] ;

[0028] In the formula, For Grid The water production; is the actual annual evapotranspiration of each grid cell; is the average annual precipitation for each grid; is the water conservation capacity of each grid unit; is the surface runoff velocity coefficient for each grid; It represents the saturated hydraulic conductivity of soil; is the topographic index; The water yield of the grid cell calculated by the water yield module in the InVEST model;

[0029] The calculation formula of habitat quality service is as follows:

[0030] ;

[0031] In the formula, Indicates land use type Habitat quality index of the grid cell; is habitat suitability; Indicates the default parameters in the model Next, the grid cell the extent of habitat degradation in is the half-saturation parameter;

[0032] The calculation formula for carbon storage services is as follows:

[0033] ;

[0034] In the formula, is the total carbon density of the land use type; 、 、 and represent aboveground, underground, soil and dead organic matter carbon pools respectively;

[0035] The calculation formula for food supply services is as follows:

[0036] ;

[0037] In the formula, For Region Inner Grid of food production; The grid for cultivated land type The normalized difference vegetation index; Indicates area The sum of the normalized difference vegetation index of all cultivated land units within; It is the total annual output of grain and cash crops.

[0038] 3. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 2, characterized in that:

[0039] In step S2), the rainfall erosivity factor, soil erodibility factor, slope length factor, soil conservation measure factor, vegetation cover management factor, annual average potential evapotranspiration, plant available water content, regional seasonal constant factor, soil saturated hydraulic conductivity, terrain index, habitat suitability and sensitivity, and regional forest carbon pool are optimized, including:

[0040] Based on the daily observation data of regional meteorological stations, the ANUSPLIN thin disk smoothing spline interpolation method is used to correct the precipitation erosivity;

[0041]

[0042] In the formula, and They represent the average monthly precipitation and the average annual precipitation respectively;

[0043] Optimize regional soil erodibility according to the erosion / productivity model based on the Chinese soil database;

[0044]

[0045]

[0046]

[0047]

[0048] In the formula, is the soil erodibility in the EPIC model; for the revised soil erodibility; , , and Respectively represent the contents of sand, silt, clay and organic carbon in the soil (%);

[0049] Correction of regional slope-slope length factor based on DEM products;

[0050]

[0051]

[0052]

[0053]

[0054] In the formula, and are slope length and slope factor respectively; represents the horizontal slope length; is the slope length index; is the slope;

[0055] Based on NDVI data, the vegetation cover management factors are optimized according to the relationship model between vegetation growth conditions and vegetation cover management factors;

[0056]

[0057] In the formula, and They represent vegetation cover management factor and normalized difference vegetation index respectively.

[0058] Based on DEM products, the soil conservation measures factors of different land use types are corrected according to the Wener formula;

[0059]

[0060] In the formula, Indicates slope.

[0061] Based on the daily observation data of the meteorological station, the annual average potential evapotranspiration is optimized according to the Hamon model;

[0062]

[0063]

[0064]

[0065] In the formula, represents the average annual potential evapotranspiration; is the average potential evapotranspiration for each month; For the Number of days in the month; Indicates the number of sunshine hours per day in a 12-hour clock; Indicates the saturated water vapor density for each month; Indicates Average monthly temperature for the month;

[0066] Based on the Chinese soil database, the plant available water content is corrected according to the soil distribution conditions;

[0067]

[0068]

[0069] In the formula, It indicates the amount of water available to plants; , , , They respectively represent the contents of sand, silt, clay and organic carbon in the surface soil layer.

[0070] Based on the daily observation data of the meteorological station, the regional seasonal constant factor is optimized according to the Zhang coefficient model;

[0071]

[0072] In the formula, is the average annual precipitation; The amount of water available to plants.

[0073] Based on the Chinese soil database, the SoilWater Characteristics module in the SPAW (Soil-Plant-Air-Water) simulation tool is used to correct the soil saturated hydraulic conductivity;

[0074] Based on the DEM product, the terrain index was corrected using the "Watershed Delineator" module in the SWAT (Soil and Water Assessment Tool, open source software) model;

[0075]

[0076] In the formula, is the catchment area of ​​the watershed, which is generated based on the “WatershedDelineator” module in the SWAT model; Indicates soil depth; Indicates slope.

[0077] Based on the special landform types in the study area, saline-alkali land and tidal flats were introduced to optimize habitat suitability and sensitivity;

[0078] Based on the statistical data of the study area and the ninth forest inventory report, the regional forest carbon pool was improved according to the IPCC biomass expansion factor method;

[0079]

[0080]

[0081]

[0082]

[0083] In the formula, for carbon storage; and represent total biomass and carbon content, respectively; is the aboveground biomass of the tree layer; Indicates tree species The amount of accumulation; Represents the biomass expansion factor of various tree species; is the wood density of the tree; The total number of major tree species in the forest system; is the underground biomass of the tree layer; For tree species Root-to-shoot ratio; Represents the aboveground biomass of the tree layer.

[0084] 4. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 3, characterized in that:

[0085] In step S3), the three years with a ten-digit number upward from the current year are taken as historical time series nodes, and the first year with a ten-digit number downward from the current year is taken as a future time series node.

[0086] 5. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 4, characterized in that:

[0087] Step S7) constructing the system dynamics model specifically includes:

[0088] S71) Determine the spatial boundary, time boundary and simulation step size of the system dynamics model;

[0089] S72) establishing subsystems of a system dynamics model, wherein the subsystems include a socio-economic subsystem, a population subsystem, a land use subsystem, and an environmental subsystem;

[0090] S73) Select indicator variables for the system dynamics model;

[0091] S74) designing multiple scenario variables according to the development scenarios, wherein the multiple scenario variables are used to control the simulation trend of the system dynamics model under each development scenario;

[0092] S75) determining a causal relationship expression of each indicator variable in the system network model, wherein the expression method includes linear fitting, WITH LOOKUP function, INTEG function and specific equation;

[0093] S76) Extract relevant data of historical time series nodes, use Vensim PLE software and adopt precision test parameter adjustment method to build the optimal SD model, and dynamically predict the land use demand in different scenarios in the future.

[0094] 6. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 5, characterized in that:

[0095] Step S8) specifically includes:

[0096] S81) Select driving factors that affect the change of land use spatial pattern;

[0097] S82) Based on the land expansion analysis strategy module and combined with data mining, the inertia probability of land use types is obtained, and the relationship model between land use types and various driving factors is obtained using the random forest algorithm;

[0098] S83) A PLUS model was established by coupling the system dynamics model and the cellular automaton model based on multi-class random patch seeds. Meanwhile, the neighborhood weight parameters and transfer cost matrix of the PLUS model were improved to simulate the spatial pattern of land use in the target area under various development scenarios at future time nodes.

[0099] 7. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 6, characterized in that:

[0100] The expression of the PLUS model is:

[0101] Land expansion analysis strategy module:

[0102]

[0103] In the formula, For the In the grid The development potential of this type of land; Represents the decision tree volume; is the vector composed of all driving factors; Represents the indicator function of the decision tree; for The vector is The land type predicted by the decision tree; Whether to allow other land types to Type conversion, the value is 1 or 0.

[0104] Cellular automaton model of multiple types of random patch seeds:

[0105]

[0106]

[0107]

[0108] In the formula, for At the moment Towards The comprehensive probability of land class conversion; For space Land use type The conversion probability of A random value between 0 and 1; is the threshold for generating new patches; represents the impact of future land use demand on type k of land use; is the neighborhood effect; is the number of cellular units; For land use type The total number of grid cells; is the number of iterations of the model during the dynamic simulation; Indicates the land use type; Represents the difference between the current stock and future demand for each type of land at the corresponding number of iterations.

[0109] 8. The multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 7, characterized in that:

[0110] Based on the landscape index, the total area change is used to improve the neighborhood weight parameter in the prediction process. The expression is:

[0111]

[0112] In the formula, For the Neighborhood weight parameter of the land use class; represents the expansion area of ​​each land use type; and They respectively represent the maximum and minimum areas of land type expansion.

[0113] 9. The multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 8, characterized in that:

[0114] In step S9), the trade-off-synergy relationship of each indicator under different scales and scenarios is demonstrated through Spearman correlation analysis, bivariate local spatial autocorrelation index and self-organizing map method.

[0115] Due to the adoption of the above technical solution, the beneficial effects of this application are as follows:

[0116] (1) This application has made regionalized improvements to the key ecological parameters of the InVEST model for assessing various ecosystem services, which can improve the applicability and accuracy of the quantitative assessment method of ecosystem services and effectively avoid the drawbacks of using existing product information as basic data, which may affect the accuracy of the assessment results.

[0117] (2) This application proposes various development scenarios coupled with future climate change, and couples the system dynamics model and the PLUS model based on the various development scenarios. The comprehensive model can predict the spatial pattern of land use in the target area under various development scenarios in a certain period in the future. The prediction results comprehensively consider the impact of the current social and economic development status, population changes, climate and environmental changes, and industrial structure characteristics on land use changes. The coupling of the system dynamics model and the PLUS model greatly compensates for the subjective defects in the current land use prediction process, making the prediction results more in line with the objective development situation and improving the accuracy, scientificity and applicability of the prediction.

[0118] (3) Based on the prediction results of the system dynamics model and the PLUS model, the present invention further couples the InVEST model, making use of the complementary advantages of the three models to make up for the limitations of a single model. At the same time, the key ecological parameters are improved in scenarios, and high-precision comprehensive ecosystem services suitable for the future multi-scenario climate framework are simulated. The trade-off / synergy relationship of ecosystem services under different scenarios is explored from a multi-scale perspective. The present invention not only provides a new simulation idea for predicting future ecosystem services, but also provides a multi-scale and multi-model vision basis for regional ecological restoration, ecological protection zone demarcation, and rational use of spatial land resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0119] The drawings described herein are used to provide a further understanding of the present application and constitute a part of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0120] Figure 1 A flow chart of a multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and assessment method provided for this application;

[0121] Figure 2 This is a map of the historical spatial distribution and change characteristics of the functions of each ecological service system in this application;

[0122] Figure 3 This is a distribution map of multiple ecosystem service hotspots in the embodiment of this application;

[0123] Figure 4 This is a composition diagram of multiple hotspots of ecosystem services in the embodiment of this application;

[0124] Figure 5 This is a multi-scenario ecosystem service prediction result diagram of an embodiment of the present application;

[0125] Figure 6 This is a result diagram of the trade-off / synergy relationship of ecosystem services in multiple scenarios at the global scale in the embodiment of this application;

[0126] Figure 7 Result diagram of trade-off / synergy relationship of ecosystem services in multiple scenarios at local scale in the embodiment of the present application;

[0127] Figure 8 Result diagram of watershed-scale multi-scenario ecosystem service trade-offs / synergies in the embodiments of the present application. DETAILED DESCRIPTION

[0128] Based on the background technology, this application aims to address the shortcomings of the existing technology in that it fails to provide a high-precision comprehensive ecosystem service assessment method suitable for predicting future multi-scenario climate frameworks, and does not consider the impact of changes in social factors on ecosystem services. Figure 1 As shown, the present application provides a multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method, comprising the following steps:

[0129] S1) Select administrative units, river basins or urban agglomerations as target areas for assessing ecosystem service trade-offs and synergies.

[0130] In the specific implementation, the selection of target areas needs to comprehensively consider the current state of social development, the challenges faced in the construction of ecological civilization, the important strategic significance of the country and the availability of some data.

[0131] In this embodiment, the determined target area is the urban agglomeration of Shandong Province, including Heze, Jining, Tai'an, Liaocheng, Jinan, Dezhou, Binzhou, Zibo and Dongying.

[0132] S2) Soil conservation services, water conservation services, habitat quality services, carbon storage services and food supply services are used as key evaluation indicators for assessing the target area.

[0133] The construction of the ecosystem service evaluation system should be closely related to the maintenance of ecological security and overall sustainable development in the study area, and can reflect the typical ecological problems in the region; the selected ecosystem service evaluation indicators should play the greatest role in the precise implementation of ecological restoration and ecological protection policies in the region.

[0134] In step S2), the soil conservation service is calculated as follows:

[0135] ;

[0136] ;

[0137] ;

[0138] In the formula, For Grid soil conservation capacity; is the total potential soil erosion of land without management measures; is the actual amount of soil erosion after taking soil conservation measures; For Grid Sediment retention capacity; is the rainfall erosivity factor; is the soil erodibility factor; is the slope length factor; and They represent soil conservation measures factor and vegetation cover management factor respectively.

[0139] The calculation formula for water conservation services is as follows:

[0140]

[0141] ;

[0142] In the formula, For Grid The water production; is the actual annual evapotranspiration of each grid cell in the urban agglomeration; is the average annual precipitation for each grid; is the water conservation capacity of each grid unit; is the surface runoff velocity coefficient for each grid; It represents the saturated hydraulic conductivity of soil; is the topographic index; The water yield of the grid cell calculated by the water yield module in the InVEST model.

[0143] The calculation formula of habitat quality service is as follows:

[0144] ;

[0145] In the formula, Indicates land use type Habitat quality index of the grid cell; is habitat suitability; Indicates the default parameters in the model Next, the grid cell the extent of habitat degradation in is the half-saturation parameter.

[0146] The calculation formula for carbon storage services is as follows:

[0147] ;

[0148] In the formula, is the total carbon density of the land use type; 、 、 and Represent aboveground, underground, soil and dead organic matter carbon pools respectively.

[0149] The calculation formula for food supply services is as follows:

[0150] ;

[0151] In the formula, For Region Inner Grid of food production; The grid for cultivated land type Normalized difference vegetation index; Indicates area The sum of the normalized difference vegetation index of all cultivated land units within; It is the total annual output of grain and cash crops.

[0152] In this step, the ecological parameters in the quantitative assessment process of various ecosystem services are further optimized and improved. Specifically, in step S2), the rainfall erosivity factor, soil erodibility factor, slope length factor, soil conservation measure factor, vegetation cover management factor, annual average potential evapotranspiration, plant available water content, regional seasonal constant factor, soil saturated hydraulic conductivity, terrain index, habitat suitability and sensitivity, and regional forest carbon pool are optimized, including:

[0153] Based on the daily observation data of regional meteorological stations, the ANUSPLIN thin disk smoothing spline interpolation method is used to correct the precipitation erosivity;

[0154]

[0155] In the formula, and They represent the average monthly precipitation and the average annual precipitation respectively;

[0156] Optimize regional soil erodibility according to the erosion / productivity model based on the Chinese soil database;

[0157]

[0158]

[0159]

[0160]

[0161] In the formula, is the soil erodibility in the EPIC model; for the revised soil erodibility; , , and Respectively represent the contents of sand, silt, clay and organic carbon in the soil (%);

[0162] Correction of regional slope-slope length factor based on DEM products;

[0163]

[0164]

[0165]

[0166]

[0167] In the formula, and are slope length and slope factor respectively; represents the horizontal slope length; is the slope length index; is the slope;

[0168] Based on NDVI data, the vegetation cover management factors are optimized according to the relationship model between vegetation growth conditions and vegetation cover management factors;

[0169]

[0170] In the formula, and They represent vegetation cover management factor and normalized difference vegetation index respectively.

[0171] Based on DEM products, the soil conservation measures factors of different land use types are corrected according to the Wener formula;

[0172]

[0173] In the formula, Indicates slope.

[0174] Based on the daily observation data of the meteorological station, the annual average potential evapotranspiration is optimized according to the Hamon model;

[0175]

[0176]

[0177]

[0178] In the formula, represents the average annual potential evapotranspiration; is the average potential evapotranspiration for each month; For the Number of days in the month; Indicates the number of sunshine hours per day in a 12-hour clock; Indicates the saturated water vapor density for each month; Indicates Average monthly temperature for the month;

[0179] Based on the Chinese soil database, the plant available water content is corrected according to the soil distribution conditions;

[0180]

[0181]

[0182] In the formula, It indicates the amount of water available to plants; , , , They respectively represent the contents of sand, silt, clay and organic carbon in the surface soil layer.

[0183] Based on the daily observation data of the meteorological station, the regional seasonal constant factor is optimized according to the Zhang coefficient model;

[0184]

[0185] In the formula, is the average annual precipitation; The amount of water available to plants.

[0186] Based on the Chinese soil database, the SoilWater Characteristics module in the SPAW (Soil-Plant-Air-Water) simulation tool is used to correct the soil saturated hydraulic conductivity;

[0187] Based on the DEM product, the terrain index was corrected using the "Watershed Delineator" module in the SWAT (Soil and Water Assessment Tool, open source software) model;

[0188]

[0189] In the formula, is the catchment area of ​​the watershed, which is generated based on the “WatershedDelineator” module in the SWAT model; Indicates soil depth; Indicates slope.

[0190] Based on the special landform types in the study area, saline-alkali land and tidal flats were introduced to optimize habitat suitability and sensitivity;

[0191] Based on the statistical data of the study area and the ninth forest inventory report, the regional forest carbon pool was improved according to the IPCC biomass expansion factor method;

[0192]

[0193]

[0194]

[0195]

[0196] In the formula, for carbon storage; and represent total biomass and carbon content, respectively; is the aboveground biomass of the tree layer; Indicates tree species The amount of accumulation; Represents the biomass expansion factor of various tree species; is the wood density of the tree; The total number of major tree species in the forest system; is the underground biomass of the tree layer; For tree species Root-to-shoot ratio; Represents the aboveground biomass of the tree layer.

[0197] S3) Selecting key timing nodes; the key timing nodes include historical timing nodes and future timing nodes.

[0198] In step S3), the three years with a ten-digit number upward from the current year are taken as historical time series nodes, and the first year with a ten-digit number downward from the current year is taken as a future time series node.

[0199] In this embodiment, the spatiotemporal differentiation characteristics of various ecosystem services in urban agglomerations are evaluated based on historical key time series nodes, and then the simulation time interval is determined based on the change trend of the current situation. In this embodiment, 2000, 2010 and 2020 are selected as historical key time nodes to quantitatively evaluate the spatial distribution of ecosystem services, and then the simulation interval is further determined to be 10 years based on the historical changes in ecosystem services, that is, 2030 is used as the key year for future simulation prediction.

[0200] S4) Obtain multi-source remote sensing data and social statistical data in the target area, input the multi-source remote sensing data and social statistical data into the InVEST model, and obtain the spatial differentiation characteristics of ecosystem service evaluation indicators and their key time series nodes.

[0201] SC (soil conservation service), WC (water conservation service), HQ (habitat quality service), CF (carbon storage service) and FS (food supply service) were obtained through the InVEST model assessment. The spatial heterogeneity characteristics of SC, WC, HQ, CF and FS were quantitatively assessed from the perspective of the entire region and the county in the historical time series, and the results were verified and the hotspot area assessment and analysis were carried out. Among them, the scope of hotspot areas with one or more types of ecosystem services is defined as type I to V hotspot area types. As shown in the attached figure Figure 2 As shown in Figure 1, at the global scale, soil conservation capacity first increased and then decreased, water conservation and food supply services increased year by year, and habitat quality and carbon sequestration capacity decreased year by year. At the county scale, the counties in the central mountainous areas of the urban agglomeration (Zhangqiu District, Laiwu District, Licheng District, Yiyuan County, etc.) and the southern mountainous areas (Sishui County, Xintai City, Ningyang County, etc.) carry a variety of high-value ecosystem service areas, and the growth of ecosystem services is concentrated in the above counties. As shown in the attached Figure 3 and attached Figure 4As shown in the figure, the area of ​​multiple ecosystem service hotspots first increases and then decreases. The Class I hotspot is the most widely distributed, and the total area of ​​Class III, IV and V increases year by year, mainly concentrated in the core areas of the Shandong Mountains and southern hills. Cultivated land, water areas and construction land are the main contributors to Class I and Class II hotspots, and typical ecological land such as cultivated land, forest land and grassland are the main contributors to Class III and above hotspots. At the same time, 87.99%, 81.7% and 35.21% of unused land, construction land and cultivated land do not significantly provide any ecosystem services.

[0202] S5) The local spatial autocorrelation index, land use product and landscape pattern index based on the common edge adjacency spatial weights are used to obtain the historical change patterns of key evaluation indicators in the target area.

[0203] The historical change law analysis is about the changing trend of various indicators in the historical period (for example: increasing or decreasing year by year, and to what extent), which is objective, because it is followed by the predicted future time series indicators. In order to see the changing law of indicators from 2000 to 2030, it is necessary to use the historical time series as a reference period to analyze the changing trend of future indicators.

[0204] S6) Establishing development scenarios, including natural development scenarios, ecological protection scenarios, and economic development scenarios;

[0205] The natural development scenario is based on the average annual temperature, average annual precipitation, population and GDP data under the future climate framework;

[0206] The ecological protection scenario is based on the ecological land area, annual average temperature, annual average precipitation, population and GDP data under the future climate framework;

[0207] The economic development scenario is based on the cultivated land red line, the annual average temperature, annual average precipitation, population and GDP data under the future climate framework.

[0208] In this embodiment, the detailed designs of different scenarios are as follows:

[0209] Natural development scenario (NDS). Following the historical land use development law, without considering the impact of policy planning on land use, the future land use transfer characteristics are basically consistent with that of this period, and no additional restrictions are added. This scenario couples the average annual temperature, average annual precipitation and population under the SSP245 (SSP2-RCP4.5) scenario, the GDP under the SSP4 scenario, and the moderate social development speed, which generally represents a medium social vulnerability and medium radiative forcing scenario.

[0210] Ecological protection scenario (EPS). With ecological protection as the main goal, restrict the conversion of ecological land such as forest, grassland and water to other land types, and promote the conversion of unused land and construction land to ecological land. Coupled with the average annual temperature, average annual precipitation and population under the SSP126 (SSP1-RCP2.6) scenario, GDP under the SSP1 pathway. This scenario strictly implements the requirements of ecological protection policies, and generally represents the low vulnerability and low radiative forcing scenarios of CIMP6.

[0211] Economic development scenario (EDS). Based on the natural development scenario, combined with the reality of the rapid economic development stage, while adhering to the 1.8 billion mu of cultivated land red line. Restrict the transfer of cultivated land to land types other than construction land, and couple the average annual temperature, average annual precipitation and population under the SSP585 (SSP5-RCP8.5) scenario, and the GDP under the SSP5 path. It represents the CIMP6 high-forcing scenario, with the goal of promoting rapid urban development at the expense of fossil energy combustion.

[0212] S7) constructing a system dynamics model SD, and predicting the land use demand area of ​​the target area under different development scenarios based on the system dynamics model.

[0213] In this embodiment, the land use classification area of ​​the urban agglomeration in 2010 is extracted and the relevant statistical data from 2010 to 2020 are obtained. The causal relationship equation of the entire SD system network is constructed with the help of Vensim PLE to simulate the land use demand area in 2020, and then the actual land use classification area in 2020 is extracted to verify the prediction accuracy of the SD model. After verification, except for unused land, the relative errors of the remaining land types are within 5%, with high accuracy. Then, the numerical values ​​of relevant variables and multi-scenario variables are adjusted to dynamically simulate the land use demand area of ​​the urban agglomeration under the NDS, EPS and EDS scenarios in 2030.

[0214] S8) constructing a PLUS model based on the land use demand area under each development scenario obtained in step S7), and coupling the SD-PLUS model to predict the land use spatial pattern of the target area under each development scenario at future time series nodes;

[0215] The PLUS model is established based on a land use simulation model of a cellular automaton model and a random patch generation mechanism.

[0216] Step S7) constructing the system dynamics model specifically includes:

[0217] S71) Determine the spatial boundary, time boundary and simulation step size of the system dynamics model;

[0218] S72) establishing subsystems of a system dynamics model, wherein the subsystems include a socio-economic subsystem, a population subsystem, a land use subsystem, and an environmental subsystem;

[0219] S73) Select indicator variables for the system dynamics model;

[0220] S74) designing multiple scenario variables according to the development scenarios, wherein the multiple scenario variables are used to control the simulation trend of the system dynamics model under each development scenario;

[0221] S75) determining a causal relationship expression of each indicator variable in the system network model, wherein the expression method includes linear fitting, WITH LOOKUP function, INTEG function and specific equation;

[0222] S76) Extract relevant data of historical time series nodes, use Vensim PLE software and adopt precision test parameter adjustment method to build the optimal SD model, and dynamically predict the land use demand in different scenarios in the future.

[0223] In this embodiment, the expansion potential of each land use type in two time series was explored through the LEAS module of the PLUS model based on the land use in 2010 and 2020. A total of 16 driving factors, including terrain factors, soil factors, climate and environmental factors, and socio-economic factors, which have the greatest impact on the land use change of urban agglomerations, were selected as driving factors for exploring the potential for land use expansion. The spatial distribution of land use in urban agglomerations in 2020 was simulated through the CARS module, and compared with the actual land use results for verification. Verification shows that the Kappa coefficient is 0.784, the overall accuracy reaches 0.831, and the Fom index is 0.113, with a high overall accuracy. After the accuracy meets the requirements, the base period land use is adjusted to 2020, and the multi-scenario land use demand predicted by the SD model is coupled. The land use spatial distribution of urban agglomerations in 2030 under the NDS, EPS and EDS scenarios is dynamically simulated through the LEAS and CARS modules of the PLUS model. Among them, the neighborhood weight parameter is improved from the perspective of landscape index, the transfer cost matrix is ​​improved, the transfer degree between different land types is improved, and the nature reserve restriction is set.

[0224] Step S8) specifically includes:

[0225] S81) Select driving factors that affect the change of land use spatial pattern;

[0226] S82) Based on the land expansion analysis strategy module and combined with data mining, the inertia probability of land use types is obtained, and the relationship model between land use types and various driving factors is obtained using the random forest algorithm;

[0227] S83) A PLUS model was established by coupling the system dynamics model and the cellular automaton model based on multi-class random patch seeds. At the same time, the domain weight parameters and transfer cost matrix of the PLUS model were improved to simulate the spatial pattern of land use in the target area under various development scenarios at future time nodes.

[0228] The expression of the PLUS model is:

[0229] Land expansion analysis strategy module:

[0230]

[0231] In the formula, For the In the grid The development potential of this type of land; Represents the decision tree volume; is the vector composed of all driving factors; Represents the indicator function of the decision tree; for The vector is The land type predicted by the decision tree; Whether to allow other land types to Type conversion, the value is 1 or 0.

[0232] Cellular automaton model of multiple types of random patch seeds:

[0233]

[0234]

[0235]

[0236] In the formula, for At the moment Towards The comprehensive probability of land class conversion; For space Land use type The conversion probability of A random value between 0 and 1; is the threshold for generating new patches; represents the impact of future land use demand on type k of land use; is the neighborhood effect; is the number of cellular units; For land use type The total number of grid cells; is the number of iterations of the model during the dynamic simulation; Indicates the land use type; Represents the difference between the current stock and future demand for each type of land at the corresponding number of iterations.

[0237] Based on the landscape index, the total area change is used to improve the neighborhood weight parameter in the prediction process. The expression is:

[0238]

[0239] In the formula, For the Neighborhood weight parameter of the land use class; represents the expansion area of ​​each land use type; and Respectively represent the maximum and minimum areas of land type expansion

[0240] S9) Based on different climate frameworks, the key ecological parameters for future ecosystem service prediction are revised, and the optimal InVEST model for prediction is obtained by combining the spatial pattern of land use under each development scenario obtained in step S8); the trade-off-synergy relationship between different scenarios at different scales in the region is obtained with the help of the optimal InVEST model and analysis method.

[0241] In step S9), the trade-off-synergy relationship of each indicator under different scales and scenarios is demonstrated through Spearman correlation analysis, bivariate local spatial autocorrelation index and self-organizing map method.

[0242] Based on the multi-scenario design under the climate framework, the key ecological parameters of various ecosystem services in the future are improved, including precipitation erosivity (R), runoff velocity coefficient (V), annual average potential evapotranspiration (ET) and precipitation (PRE). Then, the multi-scenario land use predicted by the SD model and the PLUS model are coupled with the InVEST model to predict the spatial distribution of key ecosystem services in the NDS, EPS and EDS scenarios of the urban agglomeration in 2030, such as Figure 5 shown.

[0243] The results show that there are obvious differences between ecosystem services under various scenarios, and different scenarios show a downward or upward trend compared to 2020, which mainly depends on the law and degree of transformation between ecological land and non-ecological land. According to the heterogeneity analysis of spatial differentiation pattern, the central mountains and southern mountains have made significant contributions to maintaining various ecosystem service functions under different scenarios, and the development orientation of service levels is significantly different due to the interaction between climate and topography.

[0244] Figure 6The trade-off / synergy relationship of ecosystem services under the global scale of multiple scenarios is shown. Under the NDS scenario, there is a high degree of synergy between habitat quality and water conservation, carbon sequestration, and water conservation and carbon sequestration, with correlation coefficients of 0.60***, 0.86***, and 0.80***, respectively. Soil conservation and food supply services show a slight trade-off relationship (-0.084***). The correlation coefficient between food supply and other ecosystem services is generally low, with an average of only 0.27***. Under the EPS scenario, the strong synergy is in the order of habitat quality and carbon storage (0.89***) > water conservation and carbon storage (0.75***) > soil conservation and water conservation (0.71***) > water conservation and habitat quality (0.59***) > soil conservation and carbon storage (0.56***). The trade-off relationship between food supply and soil conservation has been improved, from the negative correlation in the NDS scenario to a positive correlation (0.27***). Secondly, soil conservation and habitat quality, water conservation and food supply, and habitat quality and food supply all have synergistic relationships. Under the EDS scenario, water conservation and carbon storage, habitat quality and water conservation, and carbon storage and habitat quality have strong synergistic relationships, with correlation coefficients of 0.82***, 0.60***, and 0.85***, respectively, but food supply and soil conservation have a moderate trade-off relationship (-0.14***).

[0245] Figure 7 The trade-off / synergy relationship of ecosystem services at the local scale of multiple scenarios is shown. At the local scale of the three scenarios, the interrelationships of most ecosystem services are mainly low-low synergy, and high-high synergy mostly appears in the mountains of central Shandong and the southern mountains, and is highly clustered. Therefore, the high-level comprehensive services that the ecosystem in the study area can provide have certain distance limitations. Only the trade-off / synergy relationship under the ecological protection scenario is optimal. Compared with other scenarios, the high-high synergy areas near rural settlements have been effectively protected, but the high-high synergy areas in the other scenarios still show a decreasing trend to varying degrees, and some have turned into low-low synergy or trade-off states.

[0246] Figure 8 The trade-off / synergy relationship of ecosystem services at the watershed scale under multiple scenarios is shown. At the watershed scale, ecosystem service clusters under different scenarios have different ecosystem service combination characteristics and service levels. The EPS scenario is the best overall, with a significantly improved trade-off / synergy relationship and the highest service level. The ecological source area has made an important contribution to promoting the ecological development and ecological regulation of the entire watershed of the Shandong Yellow River urban agglomeration. Compared with the SC, WC, HQ and CF of NDS and EDS, EPS has stronger service capabilities and high-high synergy relationships.

[0247] Anything not described in this application can be achieved by adopting or drawing on existing technologies.

[0248] The above is only an embodiment of the present application and is not intended to limit the present application. For those skilled in the art, the present application may have various changes and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the scope of the claims of the present application.

Claims

1. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and assessment method, characterized in that: The following steps are involved: S1) Select administrative units, river basins or urban agglomerations as target areas for assessing ecosystem service trade-offs-synergies; S2) Soil conservation services, water conservation services, habitat quality services, carbon storage services and food supply services are used as key evaluation indicators for assessing the target area; S3) Selecting key timing nodes; the key timing nodes include historical timing nodes and future timing nodes; S4) Obtain multi-source remote sensing data and social statistical data in the target area, input the multi-source remote sensing data and social statistical data into the InVEST model, and obtain the spatial differentiation characteristics of ecosystem service evaluation indicators and their key time series nodes; S5) Obtain the historical change patterns of key evaluation indicators in the target area based on the local spatial autocorrelation index, land use products and landscape pattern index based on the common edge adjacency spatial weight; S6) Establishing development scenarios, including natural development scenarios, ecological protection scenarios, and economic development scenarios; The natural development scenario is based on the average annual temperature, average annual precipitation, population and GDP data under the future climate framework; The ecological protection scenario is based on the ecological land area, annual average temperature, annual average precipitation, population and GDP data under the future climate framework; The economic development scenario is based on the cultivated land red line, the annual average temperature, annual average precipitation, population and GDP data under the future climate framework; S7) constructing a system dynamics model SD, and predicting the land use demand area of ​​the target area under different development scenarios according to the system dynamics model; S8) constructing a PLUS model based on the land use demand area under each development scenario obtained in step S7), and coupling the SD-PLUS model to predict the land use spatial pattern of the target area under each development scenario at future time series nodes; The PLUS model is established based on a land use simulation model of a cellular automaton model and a random patch generation mechanism; S9) Correct the key ecological parameters for future ecosystem service prediction based on different climate frameworks, and combine the land use spatial pattern under each development scenario obtained in step S8) to obtain the optimal InVEST model for prediction; Use the optimal InVEST model and analysis method to obtain the trade-off-synergy relationship of different scenarios at different scales in the region; In step S2), the soil conservation service is calculated as follows: ； ； ； In the formula, For Grid soil conservation capacity; is the total potential soil erosion of land without management measures; is the actual amount of soil erosion after taking soil conservation measures; For Grid Sediment retention capacity; is the rainfall erosivity factor; is the soil erodibility factor; is the slope length factor; and represent soil conservation measures factor and vegetation cover management factor respectively; The calculation formula for water conservation services is as follows: ； In the formula, For Grid The water production; is the actual annual evapotranspiration of each grid cell; is the average annual precipitation for each grid; is the water conservation capacity of each grid unit; is the surface runoff velocity coefficient for each grid; It represents the saturated hydraulic conductivity of soil; is the topographic index; The water yield of the grid cell calculated by the water yield module in the InVEST model; The calculation formula of habitat quality service is as follows: ； In the formula, Indicates land use type Habitat quality index of the grid cell; is habitat suitability; Indicates the default parameters in the model Next, the grid cell the extent of habitat degradation in is the half-saturation parameter; The calculation formula for carbon storage services is as follows: ； In the formula, is the total carbon density of the land use type; 、 、 and represent aboveground, underground, soil and dead organic matter carbon pools respectively; The calculation formula for food supply services is as follows: ； In the formula, For Region Inner Grid of food production; The grid for cultivated land type Normalized difference vegetation index; Indicates area The sum of the normalized difference vegetation index of all cultivated land units within; It is the total annual output of grain and cash crops; In step S2), the rainfall erosivity factor, soil erodibility factor, slope length factor, soil conservation measure factor, vegetation cover management factor, annual average potential evapotranspiration, plant available water content, regional seasonal constant factor, soil saturated hydraulic conductivity, terrain index, habitat suitability and sensitivity, and regional forest carbon pool are optimized, including: Based on the daily observation data of regional meteorological stations, the ANUSPLIN thin disk smoothing spline interpolation method is used to correct the precipitation erosivity; In the formula, and They represent the average monthly precipitation and the average annual precipitation respectively; Optimize regional soil erodibility according to the erosion / productivity model based on the Chinese soil database; In the formula, is the soil erodibility in the EPIC model; for the revised soil erodibility; , , and Respectively represent the contents of sand, silt, clay and organic carbon in the soil; Correction of regional slope-slope length factor based on DEM products; In the formula, and are slope length and slope factor respectively; represents the horizontal slope length; is the slope length index; is the slope; Based on NDVI data, the vegetation cover management factors are optimized according to the relationship model between vegetation growth conditions and vegetation cover management factors; In the formula, and represent vegetation cover management factor and normalized difference vegetation index respectively; Based on DEM products, the soil conservation measures factors of different land use types are corrected according to the Wener formula; In the formula, Indicates slope; Based on the daily observation data of the meteorological station, the annual average potential evapotranspiration is optimized according to the Hamon model; In the formula, represents the average annual potential evapotranspiration; is the average potential evapotranspiration for each month; For the Number of days in the month; Indicates the number of sunshine hours per day in a 12-hour clock; Indicates the saturated water vapor density for each month; Indicates Average monthly temperature for the month; Based on the Chinese soil database, the plant available water content is corrected according to the soil distribution conditions; In the formula, It indicates the amount of water available to plants; , , , Respectively represent the contents of sand, silt, clay and organic carbon in the surface soil layer; Based on the daily observation data of the meteorological station, the regional seasonal constant factor is optimized according to the Zhang coefficient model; In the formula, is the average annual precipitation; is the amount of water available to plants; Based on the Chinese soil database, the Soil Water Characteristics module in the SPAW (Soil-Plant-Air-Water) simulation tool is used to correct the soil saturated hydraulic conductivity; Based on the DEM product, the terrain index is corrected using the "Watershed Delineator" module in the SWAT (Soil and Water Assessment Tool, open source software) model; In the formula, is the catchment area of ​​the watershed, which is generated based on the "WatershedDelineator" module in the SWAT model; Indicates soil depth; Indicates slope; Based on the special landform types in the study area, saline-alkali land and tidal flats were introduced to optimize habitat suitability and sensitivity; Based on the statistical data of the study area and the ninth forest inventory report, the regional forest carbon pool was improved according to the IPCC biomass expansion factor method; In the formula, for carbon storage; and represent total biomass and carbon content, respectively; is the aboveground biomass of the tree layer; Indicates tree species The amount of accumulation; Represents the biomass expansion factor of various tree species; is the wood density of the tree; The total number of major tree species in the forest system; is the underground biomass of the tree layer; For tree species Root-to-shoot ratio; Represents the aboveground biomass of the tree layer.

2. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 1, characterized in that: In step S3), the three years with a ten-digit number upward from the current year are taken as historical time series nodes, and the first year with a ten-digit number downward from the current year is taken as a future time series node.

3. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 2, characterized in that: Step S7) constructing the system dynamics model specifically includes: S71) Determine the spatial boundary, time boundary and simulation step size of the system dynamics model; S72) establishing subsystems of a system dynamics model, wherein the subsystems include a socio-economic subsystem, a population subsystem, a land use subsystem, and an environmental subsystem; S73) Select indicator variables for the system dynamics model; S74) designing multiple scenario variables according to the development scenarios, wherein the multiple scenario variables are used to control the simulation trend of the system dynamics model under each development scenario; S75) determining a causal relationship expression of each indicator variable in the system network model, wherein the expression method includes linear fitting, WITH LOOKUP function, INTEG function and specific equation; S76) Extract relevant data of historical time series nodes, use Vensim PLE software and precision test parameter adjustment method to build the optimal SD model, and dynamically predict the land use demand in different scenarios in the future.

4. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 3, characterized in that: Step S8) specifically includes: S81) Select driving factors that affect the change of land use spatial pattern; S82) Based on the land expansion analysis strategy module and combined with data mining, the inertia probability of land use types is obtained, and the relationship model between land use types and various driving factors is obtained using the random forest algorithm; S83) A PLUS model was established by coupling the system dynamics model and the cellular automaton model based on multi-class random patch seeds. Meanwhile, the neighborhood weight parameters and transfer cost matrix of the PLUS model were improved to simulate the spatial pattern of land use in the target area under various development scenarios at future time nodes.

5. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 4, characterized in that: The expression of the PLUS model is: Land expansion analysis strategy module: In the formula, For the In the grid The development potential of this type of land; Represents the decision tree volume; is the vector composed of all driving factors; Represents the indicator function of the decision tree; for The vector in The land type predicted by the decision tree; Whether to allow other land types to Type conversion, the value is 1 or 0; Cellular automaton model of multiple types of random patch seeds: In the formula, for At the moment Towards The comprehensive probability of land class conversion; For space Land use type The conversion probability of A random value between 0 and 1; is the threshold for generating new patches; represents the impact of future land use demand on type k of land use; is the neighborhood effect; is the number of cellular units; For land use type The total number of grid cells; is the number of iterations of the model during the dynamic simulation; Indicates the land use type; Represents the difference between the current stock and future demand for each type of land at the corresponding number of iterations.

6. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 5, characterized in that: Based on the landscape index, the total area change is used to improve the neighborhood weight parameter in the prediction process. The expression is: In the formula, For the Neighborhood weight parameter of the land use class; represents the expansion area of ​​each land use type; and They respectively represent the maximum and minimum areas of land type expansion.

7. A multi-scale and multi-scenario ecosystem service trade-off-synergy relationship prediction and evaluation method according to claim 6, characterized in that: In step S9), the trade-off-synergy relationship of each indicator under different scales and scenarios is demonstrated through Spearman correlation analysis, bivariate local spatial autocorrelation index and self-organizing map method.

Images