Urban blue-green space ecological benefit key influence factor identification method

Abstract

The invention discloses an urban blue-green space ecological benefit key influence factor identification method. Firstly, an urban blue-green space ecological benefit evaluation index set is constructed; secondly, various index data are calculated and standardized, and an urban blue-green space ecological benefit standardized data matrix is established; thirdly, determining subjective weights of the indexes by adopting an analytic hierarchy process, determining objective weights by adopting an entropy weight method, and obtaining weights of the indexes through combined weighting; and then, constructing an entropy weight-TOPSIS model to calculate a weighted urban blue-green space ecological benefit standardized data matrix, and completing urban blue-green space ecological benefit evaluation. Then, establishing an urban blue-green space ecological benefit influence factor set; and finally, based on an urban blue-green space ecological benefit evaluation result, identifying key influence factors by adopting a GWR-XGBoost model. According to the method, a multi-criterion decision analysis method, a geographical weighted regression model and a machine learning technology are combined, so that the key influence factors of the urban blue-green space ecological benefits can be quickly and accurately obtained.

Description

Technical Field

[0001] The invention belongs to the field of urban ecological environment management, and specifically relates to a method for identifying key influencing factors of the ecological benefits of urban blue-green spaces. Background Art

[0002] Urban blue-green space is composed of blue space and green space in the city. It is the ecological foundation of the city and an important part of the current high-quality urban development. It has become a broad consensus for global cities to cope with the challenges of climate change and rapid urbanization to regard the blue-green system of the city as spatial infrastructure to coordinate the development of ecology, economy, society and other aspects. At present, "ecological benefits" have been clearly included in the evaluation index system of social and economic development and have received widespread attention from the whole society. In the limited urban space, how to effectively utilize the limited blue-green space to provide better ecological benefits has become an urgent problem to be solved.

[0003] However, existing ecological benefit evaluations mostly focus on large-scale spaces, and indicators and methods are not applicable to urban and zoning scales. Research on urban blue-green spaces is still in its infancy. Ecological benefit evaluation methods constructed in previous studies mostly focus on a single land use type, or focus on a single goal, lacking comprehensiveness and coordination. The evaluation methods are not suitable for multi-objective blue-green space ecological benefit evaluation, and cannot accurately identify multi-objective factors that affect ecological benefits. To this end, the purpose of the present invention is to provide a method for identifying key influencing factors of the ecological benefits of urban blue-green spaces, which can quickly and accurately obtain the key influencing factors of the ecological benefits of urban blue-green spaces based on a comprehensive evaluation of the ecological benefits of urban blue-green spaces. Summary of the invention

[0004] To solve the above problems, the present invention discloses a method for identifying key influencing factors of the ecological benefits of urban blue-green spaces. By combining a multi-criteria decision analysis method, a geographically weighted regression model and machine learning technology, the key influencing factors of the ecological benefits of urban blue-green spaces can be quickly and accurately obtained, which is helpful for the formulation of strategies for the optimization and reorganization of blue-green spaces.

[0005] To achieve the above object, the technical solution of the present invention is as follows:

[0006] A method for identifying key influencing factors of the ecological benefits of urban blue-green space, comprising the following steps:

[0007] S1, construct an evaluation index set of ecological benefits of urban blue-green space;

[0008] S2, based on the indicator set constructed in S1, calculates and standardizes the indicator data in the target area to establish a standardized data matrix of the ecological benefits of urban blue-green space;

[0009] S3, based on the standardized data matrix of ecological benefits in S2, uses the analytic hierarchy process to determine the subjective weights of indicators, uses the entropy weight method to determine the objective weights, and obtains the weights of various indicators through combined weighting;

[0010] S4, construct the entropy weight-TOPSIS model to calculate the weighted standardized data matrix of the ecological benefits of urban blue-green space, thereby completing the ecological benefit evaluation of urban blue-green space;

[0011] S5, establish a set of influencing factors of the ecological benefits of urban blue-green space, including two categories: natural factors and spatial pattern factors;

[0012] S6, based on the evaluation results of the ecological benefits of urban blue-green spaces in S4 and the influencing factor set in S5, combined with the GWR model (Geographically weighted regression) and XGBoost model (Extreme Gradient Boosting) in machine learning to identify key influencing factors.

[0013] Furthermore, in S1, the target layer of the urban blue-green space ecological benefit evaluation index set is: ecological benefit; the criterion layer includes: stormwater resilience, cooling benefit, carbon sink benefit, and biodiversity; the indicator layer includes: surface runoff, flooding risk, surface temperature, net primary productivity, carbon storage, and habitat quality.

[0014] Furthermore, in S2, the steps of establishing a standardized data matrix of the ecological benefits of urban blue-green spaces include:

[0015] S21, obtaining remote sensing image data, land use data, normalized difference vegetation index NDVI data, vegetation type data, temperature data, precipitation data and daily radiation data, urban soil type data, and urban land cover data in the target area;

[0016] S22, based on the land use data, urban soil type data, urban surface cover data, and precipitation data obtained in S21, the surface runoff volume of the target area is calculated according to the SCS-CN model (Soil Conservation Service Curve Number, runoff curve model);

[0017] S23, based on the remote sensing images of continuous sunny days and after rain in S21 remote sensing image data, the density of waterlogging points is calculated to characterize the flooding risk in the study area;

[0018] S24, based on the spectral data in the remote sensing image data of S21, the surface temperature of the target area is inverted using the atmospheric correction method;

[0019] S25, using the land use data, normalized difference vegetation index NDVI data, vegetation type data, temperature data, precipitation data and daily radiation data obtained in S21, the CASA (Carnegie-Ames-Stanford Approach) model is used to calculate the net primary productivity of the target area;

[0020] S26, based on the land use data of S21, uses the carbon storage module of the InVEST model (Integrated Valuation of Ecosystem Services and Trade-offs) to generate the spatial distribution of carbon storage;

[0021] S27, the Habitat Quality module in the InVEST model was used for analysis to assess the habitat quality of the target area.

[0022] Furthermore, in S2, the specific calculation steps of the standardized data matrix of ecological benefits of urban blue-green space in the target area are:

[0023] S221, calculate the maximum possible water storage capacity based on the CN values ​​(Curve Numbers) of urban and residential land, and assign values ​​to soils of different urban land types in the target area in combination with land use data;

[0024] The soil cover composite number calculation formula is as follows:

[0025]

[0026] In the formula, S is the maximum possible water storage capacity, which indicates the maximum potential water storage capacity of the soil;

[0027] S222, combined with S21 precipitation data, calculates the surface runoff in the target area according to the SCS-CN runoff model calculation formula; the SCS-CN runoff model calculation formula is as follows:

[0028]

[0029] In the formula, Q represents the surface runoff, in millimeters; P represents the rainfall, also in millimeters; λ represents the initial loss rate, which is a dimensionless value, usually calculated using λ=0.2;

[0030] Secondly, in S23, the calculation steps of the flooding risk of the target area are:

[0031] S231, based on the remote sensing image data of S21, remote sensing images under the scenarios of continuous sunny days and after rain are selected to classify land use;

[0032] S232, selecting the water body classification data of the land use classification results under the continuous sunny day and after rain scenarios in S231, removing the permanent water area in the target area, obtaining the seasonal water body, and obtaining the waterlogged area after subtraction;

[0033] S233, after converting the waterlogged area into points, the point density was calculated using the “point density analysis” tool in the GIS software to obtain the waterlogged point density as the flooding risk of the waterlogged area;

[0034] Again, in S24, the steps for obtaining the surface temperature are:

[0035] S241, through the spectral data containing surface temperature information in the S21 remote sensing image data, the "FLAASH Atmospheric Correction" tool of ENVI software is used to perform atmospheric correction processing on the remote sensing image data based on the spectral band;

[0036] S242, applying a temperature inversion algorithm to convert the corrected remote sensing data into surface temperature, and using Planck's formula to calculate the true surface temperature;

[0037] Among them, the expression of the thermal infrared radiation brightness value received by the satellite sensor can be written as:

[0038] L λ =[ε·B(T)+(1-ε)L↓]·τ+L↑

[0039] Among them, the upward radiation brightness of the atmosphere L↑, ε is the surface emissivity, T is the true surface temperature, B(T) is the thermal radiation brightness of the black body at T derived from Planck's law, τ is the transmittance of the atmosphere in the thermal infrared band, then the radiation brightness B(T) of the black body with temperature T in the thermal infrared band is:

[0040]

[0041] The Planck calculation formula is as follows:

[0042]

[0043] Where K1 and K2 are Landsat8 sensor specific calibration constants.

[0044] Then, in S25, the steps to calculate the net primary productivity are:

[0045] S251, in GIS software, the monthly precipitation, temperature and solar radiation data of the target year are subjected to Kriging interpolation processing, and band synthesis is performed to obtain an annual average grid map; then, the grid map of precipitation, temperature and solar radiation data is obtained by clipping according to the target area range;

[0046] S252, using the CASA model, input precipitation, temperature and solar radiation data, normalized vegetation index NDVI map, vegetation data, set the parameters in the static file, including the maximum and minimum values ​​of the normalized vegetation index NDVI data, the maximum light energy utilization rate, and the vegetation index SR, and obtain the distribution of net primary productivity in the target area;

[0047] Secondly, in S26, the steps for obtaining the spatial distribution of carbon storage are:

[0048] S261, based on the S21 land use data, uses the carbon storage module of the InVEST model to calculate the total carbon storage of each land type based on the different land types and corresponding carbon pool data in the target area;

[0049] S262, classifying and summarizing the calculated carbon storage data by land type to generate a spatial distribution map of carbon storage;

[0050] Among them, the formula for calculating carbon storage through the InVEST model is divided into:

[0051] C t =C a +C b +C s +C d

[0052] Where: C t is the total carbon storage of the site (t / hm 2 );C a is the aboveground carbon storage (t / hm 2 );C b is the underground carbon storage (t / hm 2 );C s is soil carbon storage (t / hm 2 );C d is the dead organic carbon stock (t / hm 2 );

[0053] Finally, in S27, the steps for assessing habitat quality in the target area are:

[0054] S271, using the habitat quality module in the InVEST model for analysis, and assessing biodiversity based on the quality of the habitat obtained;

[0055] The formula for calculating habitat quality score is:

[0056] Q xj =H j (1-D xj )

[0057] Where Q xj is the habitat quality score of grid x in land use type j, H j is the habitat suitability score of land use type j, D xj is the threat level of grid x;

[0058] S272, the model assumes that areas with good habitat quality also have high biodiversity; threat source data are factors that threaten ecological land in land use types. GIS software is used to assign a value of 1 to the threatened area and a value of 0 to the non-threat area; the formula for calculating the grid threat level is:

[0059]

[0060] In the formula, t y is the threat intensity of grid y, ω r is the relative weight of threat source r, β x is the protection level of grid x, S jr is the sensitivity of land use type j to threat source r, i rxy is the impact level of grid y in threat source r on grid x;

[0061] Among them, the calculation formula of the impact level of grid y in threat source r on grid x is:

[0062]

[0063] Where, d xy is the Euclidean distance between grid x and grid y, d rmax is the maximum impact distance of the threat source r.

[0064] Furthermore, in S3, the steps for obtaining the weights of various indicators are:

[0065] S31, establish a number of sampling points, randomly extract the calculation results of the standardized data matrix of the ecological benefits of urban blue-green space in S2, and obtain the data set.

[0066] S32, the entropy value of each indicator is calculated using the entropy weight method, and then the objective weight of each indicator is obtained.

[0067] The specific steps of using the entropy weight method to calculate the entropy value of each indicator are as follows: First, data standardization is performed. Given n evaluation indicators and m evaluation objects, the original information matrix X = (x ij )m×n, normalization is used to perform data standardization, and its expression is:

[0068]

[0069] Where Y ijare the standardized indexes of the jth coordinate of the ith evaluation index; x ij is the original data of the jth coordinate of the ith evaluation index.

[0070] Secondly, determining the feature weight is a key step in the evaluation process. In order to calculate the weight of the jth indicator of the i-th evaluation object, we first need to define the weight p of the indicator. ij , obtained through S321 ij Perform calculations.

[0071]

[0072] Then, using the weight p of each indicator ij , calculate the information entropy, calculate the information entropy E of each evaluation index j .

[0073]

[0074] In the formula, m is the number of evaluation objects, E j is the information entropy of each evaluation index.

[0075] Finally, through the information entropy E j , use the entropy weight method to determine the indicator weight and calculate the objective weight β of the evaluation indicator j =β1,β2,…,β n , where β j is the weight corresponding to the jth indicator.

[0076]

[0077] In the formula, n is the number of evaluation indicators, β j is the entropy weight of the entropy weight method.

[0078] S33, in the AHP hierarchical analysis, the four ecological goals of stormwater resilience, cooling benefits, carbon sink benefits, and biodiversity are set to be equal weights. By comprehensively considering the calculation results of the AHP hierarchical analysis method and the entropy weight method, the final weights of each indicator are obtained. Among them, the final weight formula of each indicator is:

[0079]

[0080] In the formula, α j is the weight of the AHP hierarchical analysis method corresponding to the jth indicator, β j is the weight of the entropy weight method corresponding to the jth indicator.

[0081] Furthermore, in S4, the specific steps of constructing the entropy weight-TOPSIS model to calculate the weighted standardized data matrix of the ecological benefits of urban blue-green space are as follows:

[0082] S41, write a program in Matlab software, randomly select several sampling points, combine the indicator weight results of S3 with the calculation results of the standardized data matrix of the ecological benefits of urban blue-green space in S2, and calculate the ecological benefits of each point. According to the TOPSIS evaluation method, the positive and negative ideal solutions and the evaluation results of the ecological benefits corresponding to the sampling points are obtained.

[0083] Among them, the calculation steps of the TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution) evaluation method are:

[0084] First, in order to comprehensively evaluate each evaluation object, a weighted normalized decision matrix is ​​constructed according to the factor weights determined by S3. The matrix formula is:

[0085] V=X ij ·ω=[v ij ] m·n

[0086] Where V is the weight matrix, X ij is the ith evaluation index.

[0087] Secondly, using the calculated weighting matrix V, determine the positive and negative ideal solutions:

[0088]

[0089] Then, according to the positive and negative ideal solution V + 、V - , use Euclidean distance to calculate the distance between positive and negative ideal solutions:

[0090]

[0091] In the formula, Positive ideal solution distance, is the negative ideal solution distance.

[0092] Finally, the calculated positive and negative ideal solution distances are used to calculate the closeness:

[0093]

[0094] According to the formula, calculate the closeness T to the ideal solution j , T j The value range of T is 0 to 1. j The closer it is to 1, the better the ecological benefit.

[0095] S42, select the indicator weights whose evaluation results are positive ideal solutions for final calculation to obtain the final weight results; then combine the standardized data matrix of the ecological benefits of the urban blue-green space in S2 to calculate the ecological benefits of the urban blue-green space to obtain the ecological benefit evaluation results.

[0096] Furthermore, the influencing factors of the ecological benefits of urban blue-green space in S5 include two categories: "spatial pattern factors" and "natural factors". The spatial pattern factors are composed of patch layer factors and type layer factors. Among them, patch layer factors include area, perimeter, fractal dimension, near circle index, proximity index, and patch distance; type layer factors include area proportion, maximum patch index, edge density, landscape shape index, similarity adjacency ratio, aggregation, separation amount, cohesion index, and fragmentation. Natural factors are composed of elevation, slope, and slope aspect.

[0097] Furthermore, in S6, the specific steps of combining the GWR model and the XGBoost model in machine learning to identify key influencing factors are as follows:

[0098] S61, construct a data set based on the natural factors and pattern factors of the urban blue-green space in S5 and the corresponding ecological benefit evaluation results calculated in S4. In the GIS software, a random sampling tool is used to create random sampling points, extract the influencing factors and ecological benefit evaluation results to the sampling points, and eliminate the points that do not belong to the urban blue-green space to form data set 1.

[0099] S62, the spatial autocorrelation was calculated using the global Moran's index tool.

[0100] The calculation formula of Moran Index is as follows:

[0101]

[0102] In the formula, n is the total number of spatial units, y i and j Respectively represent the attribute values ​​of the i-th spatial unit and the j-th spatial unit, is the mean of all spatial unit attributes, w ij is the spatial weight value.

[0103] S63, writing a VIF calculation program in Python software to calculate the multicollinearity of each influencing factor of the data set 1 in S61.

[0104] The formula for calculating VIF is:

[0105]

[0106] In the formula, R 2 is the coefficient of determination of the model.

[0107] S64, constructing data set 2 with variables with VIF < 5 as independent variables, and constructing data set 2 with the dependent variable as the ecological benefit evaluation results, introducing the GWR model for regression analysis, using the Gaussian function as the weight function of the model, and using the AIC method to select a bandwidth that can effectively balance the bias and variance and avoid overfitting.

[0108] Among them, the GWR model calculation formula is:

[0109]

[0110] Where Y is the explanation value of the ecological benefit of city i, (u i ,v i ) is the geographical coordinates of city i; X ik is the explanatory value of the independent variable of city i, is the influencing factor; β k (u i ,v i ) Unit i area centroid (u i ,v i ) at the regression parameter; ε i is a random error term.

[0111] The calculation formula for determining the bandwidth using the AIC (Akaike Information Criterion) method is:

[0112] AIC = nlog e (RSS)+2p

[0113] Among them, n represents the number of sample points, RSS represents the residual sum of squares, and p represents the number of unknown parameters. Based on this, the AIC calculation method for fusion with the GWR model is as follows:

[0114]

[0115] in, Maximum likelihood variance estimation for the GWR model:

[0116]

[0117] S65, based on S63, takes the factor of 5≤VIF<10 as the independent variable. The dependent variable is the combination of the ecological benefit evaluation result of S4 and the residual of the GWR model in S64. A new data set 3 is constructed and imported into the XGBoost model for training. The Bayesian optimization method is used to adjust the hyperparameters of the XGBoost model.

[0118] S66, to verify the accuracy of the model, the mean absolute error, root mean square error and determination coefficient were selected as the evaluation indicators of the prediction effect, and the ten-fold cross-validation technique was used to evaluate the generalization performance of the model.

[0119] S67, the positive and negative impact of each feature was calculated using the SHAP model in Python software.

[0120] S68, in the GWR model fitting results, select R 2 The regression results with high values ​​are taken as key influencing factors.

[0121] S69, in order to further avoid the overfitting of the model, the urban blue-green spatial pattern factors with high influence and the corresponding ecological benefits after the important features in S67 are ranked are constructed into a new data set, which is interpreted after model training to obtain the importance ranking of the influencing factors.

[0122] The beneficial effects of the present invention are:

[0123] The invention provides a method for identifying key influencing factors of the ecological benefits of urban blue-green spaces, which is suitable for ecological benefit assessment and influencing factor identification at the urban scale. The weights of various ecological benefit indicators are obtained through comprehensive evaluation by the AHP method, the entropy weight method, and the TOPSIS method, which is more objective than the traditional method of determining weights and avoids the interference of human subjective factors. The introduction of machine learning, an advanced method, can accurately and quickly identify the key factors affecting the ecological benefits of urban blue-green spaces, which is helpful for the formulation of strategies for the optimization and reorganization of blue-green spaces. BRIEF DESCRIPTION OF THE DRAWINGS

[0124] Figure 1 It is a flow chart of the method for constructing the urban blue-green space ecological benefit index evaluation system of the present invention;

[0125] Figure 2 It is a specific remote sensing data area map of the implementation case of the present invention;

[0126] Figure 3 It is a surface runoff diagram of an implementation case of the present invention;

[0127] Figure 4 It is a flooding risk map of an implementation case of the present invention;

[0128] Figure 5 It is a surface temperature inversion diagram of an implementation case of the present invention;

[0129] Figure 6 is a carbon sequestration diagram of an embodiment of the present invention;

[0130] Figure 7 It is a carbon storage diagram of an implementation case of the present invention;

[0131] Figure 8It is a biodiversity diagram of an implementation case of the present invention;

[0132] Fig. 9 This is a diagram of the ecological benefit evaluation results of the implementation case of the present invention;

[0133] Fig.10 It is an elevation, slope and aspect map of the implementation case of the present invention;

[0134] Fig.11 It is a plaque layer index calculation diagram of an implementation case of the present invention;

[0135] Fig.12 It is a type layer index calculation diagram of the implementation case of the present invention;

[0136] Fig.13 This is a flow chart of the GWR-XGBoost model construction in the implementation case of the present invention;

[0137] Fig.14 This is a result diagram of spatial autocorrelation analysis of an implementation case of the present invention;

[0138] Fig.15 It is a ranking diagram of important characteristics of spatial pattern factors of implementation cases of the present invention;

[0139] Fig.16 It is a ranking diagram of important features of key spatial pattern factors in implementation cases of the present invention. DETAILED DESCRIPTION

[0140] The present invention will be further explained below in conjunction with the accompanying drawings and specific embodiments. It should be understood that the following specific embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention.

[0141] As shown in the figure, the method for identifying key influencing factors of the ecological benefits of urban blue-green space described in the present invention specifically comprises the following steps:

[0142] S1, construct an evaluation index set of ecological benefits of urban blue-green space, the target layer is ecological benefits; the criterion layer includes stormwater resilience, cooling benefits, carbon sink benefits, and biodiversity; the indicator layer includes surface runoff, flooding risk, surface temperature, carbon sequestration, carbon storage, and habitat quality. See Table 1.

[0143] Table 1 Evaluation index set of ecological benefits of urban blue-green space

[0144]

[0145] S2, select Nanjing as the case area, calculate and standardize the indicator data in the target area based on the indicator set constructed in S1, establish the standardized data matrix of the ecological benefits of urban blue-green space, and obtain the surface runoff, rain and flood inundation range, surface temperature inversion map, net primary productivity spatial distribution map, carbon storage spatial distribution map, and biodiversity distribution map of the study area. The specific steps include:

[0146] S21, obtain remote sensing image data, land use data, normalized difference vegetation index NDVI data, vegetation type data, temperature data, precipitation data and daily radiation data, urban soil type data, urban surface cover data in the target area. The case data source is shown in Table 2, and the selected Nanjing remote sensing data range is as follows: Figure 2 , the specific process is:

[0147] Table 2 Data sources

[0148]

[0149]

[0150] First, land use data: Based on the Landsat8 (2020) remote sensing images of Nanjing urban area downloaded from the geospatial cloud platform, supervised classification and manual visual interpretation and correction were carried out on the ENVI software platform. After accuracy verification, relatively accurate land use data in the study area has been obtained, with an accuracy of 30m×30m;

[0151] Then, the temperature, precipitation and daily radiation data were obtained from 13 meteorological stations in Nanjing from the Geographic Remote Sensing Ecological Network;

[0152] Secondly, vegetation type data: the vegetation type data in Nanjing was obtained from the Geographic Remote Sensing Ecological Network, with an accuracy of 30m×30m; similarly, the Normalized Difference Vegetation Index (NDVI) data: in Google Earth Engine, based on Landsat8 remote sensing data, the remote sensing images were calibrated and pre-processed to remove clouds, and the Normalized Difference Vegetation Index (NDVI) data in Nanjing was obtained according to the calculation formula of the Normalized Difference Vegetation Index (NDVI), with an accuracy of 30m×30m;

[0153] Among them, the calculation formula of normalized vegetation index NDVI is:

[0154]

[0155] Where NIR is the reflectivity of the near infrared band, and RED is the reflectivity of the red band;

[0156] Finally, soil data: In Google Earth Engine, download the soil classification data under the U.S. Department of Agriculture classification system with a data accuracy of 30m×30m.

[0157] S22, based on the urban soil type data, urban surface cover data, land use data, and meteorological data obtained in S21, the surface runoff in the case is calculated according to the SCS-CN runoff model calculation formula. The results are as follows: Figure 3 As shown, the specific process is:

[0158] Table 3 Soil data

[0159]

[0160]

[0161] First, according to the CN values ​​of urban and residential land, as shown in Table 3, the maximum possible water storage capacity is calculated, and the soil of different urban land types in the target area is assigned a value in combination with land use data. The soil cover composite number calculation formula is as follows:

[0162]

[0163] In the formula, S is the maximum possible water storage capacity (mm) and represents the maximum potential water storage capacity of the soil;

[0164] Then, based on the CN value, the soils of different urban land types in the study area were assigned values ​​using the “raster calculator” tool in the GIS software;

[0165] Finally, combined with the S21 Nanjing precipitation data, according to the SCS-CN runoff model calculation formula, the surface runoff in the study area was calculated using the "raster calculator" tool in the GIS software. The SCS-CN runoff model calculation formula is as follows:

[0166]

[0167] In the formula, Q represents the surface runoff, the unit is millimeter (mm); P represents the rainfall, the unit is also millimeter (mm); λ represents the initial loss rate, which is a dimensionless value, usually λ = 0.2 is used for calculation.

[0168] S23, based on the remote sensing images of S31 remote sensing image data under continuous sunny days and after rain, the density of waterlogging points is calculated to characterize the flooding risk in the study area. The results are as follows: Figure 4 , the specific process is:

[0169] Firstly, based on the S21 remote sensing image data, remote sensing images under continuous sunny days and after rain were selected to classify land use;

[0170] Then, the water body classification data of the land use classification results under the continuous sunny day and after-rainfall scenarios in S231 are selected, and the permanent water area in the target area is removed to obtain the seasonal water body, and the waterlogged area is obtained after subtraction;

[0171] Finally, the waterlogged areas were converted into points using GIS tools, and the point density was calculated using the “point density analysis” tool to obtain the waterlogged point density, which was used to characterize the flooding risk in the waterlogged areas.

[0172] S24, based on the spectral data in the remote sensing image data of S21, uses the atmospheric correction method to invert the surface temperature of the target area. The results are as follows: Figure 5 , the specific process is:

[0173] Firstly, the spectral data containing surface temperature information in the S21 remote sensing image data was used, and the “FLAASH Atmospheric Correction” tool of the ENVI software platform was used to perform atmospheric correction on the remote sensing image data based on spectral bands.

[0174] Secondly, the temperature inversion algorithm is used to convert the corrected remote sensing data into surface temperature, and the surface temperature is calculated using the Planck formula. Then, the accuracy of the inverted surface temperature data is verified through ground observation data, and the results are analyzed and applied;

[0175] Among them, the expression of the thermal infrared radiation brightness value received by the satellite sensor can be written as:

[0176] L λ =[ε·B(T)+(1-ε)L↓]·τ+L↑

[0177] Among them, the upward radiation brightness of the atmosphere L↑, ε is the surface emissivity, T is the true surface temperature, B(T) is the thermal radiation brightness of the black body at T derived from Planck's law, τ is the transmittance of the atmosphere in the thermal infrared band, then the radiation brightness B(T) of the black body with temperature T in the thermal infrared band is:

[0178]

[0179] The Planck calculation formula is as follows:

[0180]

[0181] Where K1 and K2 are Landsat8 sensor specific calibration constants.

[0182] S25, through the land use data, normalized difference vegetation index NDVI data, vegetation type data, temperature data, precipitation data and daily radiation data obtained in S21, the net primary productivity of the target area is calculated using the CASA model. The results are as follows: Figure 6 ;

[0183] First, the meteorological data of S21 were processed using GIS software, and Kriging interpolation was performed on the monthly precipitation, temperature, and solar radiation data of the target year;

[0184] Then, based on the above data, the 12-month raster map bands of the target year were synthesized into an annual average raster map, and resampled to a spatial resolution of 30m×30m. According to the research scope, the precipitation, temperature and solar radiation data raster maps of the target year were produced;

[0185] Finally, in the ENVI software platform, using the CASA model, input precipitation, temperature, solar radiation data, normalized vegetation index NDVI data, vegetation data, set the parameters in the static file, including the maximum and minimum values ​​of the normalized vegetation index NDVI data, the maximum light energy utilization rate, and the vegetation index SR, and obtain the spatial distribution map of net primary productivity in the actual case, such as Figure 6 .

[0186] S26, based on the land use data of S21, uses the carbon storage module of the InVEST model (Integrated Valuation of Ecosystem Services and Trade-offs) to calculate the total carbon storage of each land type based on the different land types and corresponding carbon pool data in the case, as shown in Table 4. The calculated carbon storage data are classified and summarized by land type to generate a spatial distribution map of carbon storage, as shown in Figure 7 The formula for calculating carbon storage using the InVEST model is as follows:

[0187] C t =C a +C b +C s +C d

[0188] Where: C t is the total carbon storage of the site (t / hm 2 );C a is the aboveground carbon storage (t / hm 2 );C b is the underground carbon storage (t / hm 2 );C s is soil carbon storage (t / hm 2 );Cd is the dead organic carbon stock (t / hm 2 ).

[0189] Table 4 Carbon density of each land use type in the case area

[0190]

[0191] S27, using the habitat quality module in the InVEST model for analysis, and assessing biodiversity based on the quality of the habitat obtained. The model assumes that areas with good habitat quality also have high biodiversity; threat source data are factors that threaten ecological land in land use types, and the values ​​are as shown in Table 5. Using GIS software, the threat area is assigned a value of 1 and the non-threat area is assigned a value of 0. The results are as follows Figure 8 , the specific calculation steps are:

[0192] First, the formula for calculating habitat quality score is:

[0193] Q xj =H j (1-D xj )

[0194] Where Q xj is the habitat quality score of grid x in land use type j, H j is the habitat suitability score of land use type j, D xj is the threat level of grid x;

[0195] The calculation formula of grid threat level is:

[0196]

[0197] In the formula, r y is the threat intensity of grid y, ω r is the relative weight of threat source r, β x is the protection level of grid x, S jr is the sensitivity of land use type j to threat source r, i rxy is the impact level of grid y in threat source r on grid x;

[0198] Among them, the calculation formula of the impact level of grid y in threat source r on grid x is:

[0199]

[0200] Where, d xy is the Euclidean distance between grid x and grid y, d rmax is the maximum impact distance of the threat source r.

[0201] Table 5 Threat factors for habitat quality assessment

[0202]

[0203] S3, based on the standardized data matrix of ecological benefits of S2, uses the hierarchical analysis method to determine the subjective weights of indicators, uses the entropy weight method to determine the objective weights, and obtains the weights of each indicator through combined weighting.

[0204] S31, in the GIS software, create 20,000 random points, and use the “Extract Multi-Values ​​to Points” tool to extract the calculation results of the above six indicators to points to obtain a data set;

[0205] S32, using Matlab software to write relevant programs, by using the entropy weight method to calculate the entropy value of each indicator, and then obtain the objective weight of each indicator;

[0206] S321, data standardization is an important data processing step. Given n evaluation indicators and m evaluation objects, the original information matrix X = (x ij )m×n, normalization is used to perform data standardization, and its expression is:

[0207]

[0208] Where Y ij are the standardized indexes of the jth coordinate of the ith evaluation index; x ij is the original data of the jth coordinate of the ith evaluation index;

[0209] S322, determining the feature weight is a key step in the evaluation process. In order to calculate the weight of the jth indicator of the i-th evaluation object, it is first necessary to define the weight p of the indicator. ij , obtained through S321 ij Perform calculations;

[0210]

[0211] S323, using the weights p of each indicator obtained in S322 ij , calculate the information entropy, calculate the information entropy Ej of each evaluation index;

[0212]

[0213] In the formula, m is the number of evaluation objects, E j is the information entropy of each evaluation index;

[0214] S324, through the information entropy Ej of S323, use the entropy weight method to determine the indicator weight and calculate the objective weight β of the evaluation indicator j =β1,β2,…,βn , where β j is the weight corresponding to the jth indicator;

[0215]

[0216] In the formula, n is the number of evaluation indicators, β j is the entropy weight of the entropy weight method;

[0217] S33, in the AHP hierarchical analysis, the four ecological goals of stormwater resilience, cooling benefits, carbon sink benefits, and biodiversity are set to be equal weights. By comprehensively considering the calculation results of the AHP hierarchical analysis method and the entropy weight method, the final weights of each indicator are obtained. Among them, the formula for the final weight of each indicator is obtained, and the results are shown in Table 6;

[0218] The final weight formula is:

[0219]

[0220] In the formula, α j is the weight of the AHP hierarchical analysis method corresponding to the jth indicator, β j is the weight of the entropy weight method corresponding to the jth indicator.

[0221] Table 6 Evaluation index weight calculation results

[0222]

[0223] S4, based on the quantitative evaluation results of S2 and the weight calculation results of S3, construct the entropy weight-TOPSIS model to calculate the weighted standardized data matrix of the ecological benefits of urban blue-green space, thereby completing the ecological benefit evaluation of urban blue-green space, such as Fig. 9 .

[0224] S41, write a program in Matlab, randomly select several sampling points, and calculate the ecological benefit of each point by combining the indicator weight results of S3 with the calculation results of the standardized data matrix of the ecological benefits of urban blue-green space in S2, as shown in Table 7. According to the TOPSIS evaluation method, the positive and negative ideal solutions and the evaluation results of the ecological benefits corresponding to the sampling points are obtained;

[0225] Furthermore, in order to comprehensively evaluate each evaluation object, a weighted normalized decision matrix is ​​constructed according to the factor weights determined by S3. The matrix formula is:

[0226] V=X ij ·ω=[v ij ] m·n

[0227] Where V is the weight matrix, X ij is the i-th evaluation index;

[0228] Secondly, using the calculated weighting matrix V, determine the positive and negative ideal solutions:

[0229]

[0230] Then, according to the obtained positive and negative ideal solutions V + 、V - , use Euclidean distance to calculate the distance between positive and negative ideal solutions:

[0231]

[0232] In the formula, Positive ideal solution distance, is the negative ideal solution distance;

[0233] Finally, the closeness is calculated based on the calculated positive and negative ideal solution distances:

[0234]

[0235] According to the formula, calculate the closeness T to the ideal solution j , T j The value range of T is 0 to 1. j The closer it is to 1, the better the ecological benefit;

[0236] Table 7 Calculation results of positive and negative ideal solutions

[0237]

[0238] S42, select the indicator weights whose evaluation results are positive ideal solutions for final calculation to obtain the final weight results; associate the values ​​of the above calculation results with the values ​​of each point in the geographic space, use Kriging interpolation, and then calculate the ecological benefits of the urban blue-green space in combination with the standardized data matrix of the ecological benefits of the urban blue-green space in S2 to obtain the ecological benefit evaluation results, such as Fig. 9 .

[0239] S5, establish a set of influencing factors of the ecological benefits of urban blue-green space, including two categories: natural factors and spatial pattern factors. Among them, the influencing factors of the ecological benefits of urban blue-green space include "spatial pattern factors" and "natural factors". The spatial pattern factors are composed of patch layer factors and type layer factors. Among them, patch layer factors include area, perimeter, fractal dimension, near circle index, proximity index, and patch distance; type layer factors include area proportion, maximum patch index, edge density, landscape shape index, similarity adjacency ratio, aggregation, separation amount, cohesion index, and fragmentation. Natural factors are composed of elevation, slope and slope aspect, as shown in Table 8.

[0240] Table 8 Set of factors affecting the ecological benefits of urban blue-green space

[0241]

[0242]

[0243] Furthermore, the elevation data was downloaded from the geospatial data cloud, and based on the DEM data, the slope and aspect data were obtained after analysis and processing using the slope and aspect tools in the GIS software. The data processing results are shown in Figure 2. Fig.10 As shown;

[0244] Furthermore, the Fragstats software was used to quantitatively analyze the characteristics of urban blue-green spatial patterns in the selected years. The patch layer calculation results are shown in the figure below. Fig.11 , the type layer calculation results are as follows Fig.12 .

[0245] S6, based on the urban blue-green space ecological benefit evaluation results of S4 and the influencing factor set of S5, combined with the GWR model and XGBoost model in machine learning to identify key influencing factors. The steps are as follows: Fig.13 .

[0246] The specific steps of combining the GWR model and the XGBoost model in machine learning to identify key influencing factors are:

[0247] S61, constructing data sets based on the natural factors and pattern factors of urban blue-green space in S5 and the corresponding ecological benefit evaluation results calculated in S4. In the GIS software, using random sampling tools, creating random sampling points, extracting each influencing factor and ecological benefit evaluation results to the sampling points, and eliminating points that do not belong to urban blue-green space to form data set 1;

[0248] S62, in GIS software, the global Moran index tool is used to process the ecological benefit evaluation results and calculate their spatial autocorrelation. The results are as follows Fig.14 The evaluation criteria are shown in Table 9. The spatial autocorrelation Moran index of the ecological benefit evaluation results in the case target area is 0.78, P = 0, and the absolute value of the Z-score is much greater than 2.58. This shows that the spatial distribution of ecological benefits in the case target area is not random, and there is a significant positive correlation with the spatial distribution. Therefore, the standard regression model is not applicable at this time, and the GWR model is more appropriate;

[0249] The calculation formula of Moran Index is as follows:

[0250]

[0251] In the formula, n is the total number of spatial units, y i and jRespectively represent the attribute values ​​of the i-th spatial unit and the j-th spatial unit, is the mean of all spatial unit attributes, w ij is the spatial weight value;

[0252] Table 9 Spatial autocorrelation test standards

[0253]

[0254] S63, write a VIF calculation program in Python to calculate the multicollinearity of each influencing factor of data set 1 in S61;

[0255] The formula for calculating VIF is:

[0256]

[0257] In the formula, R 2 is the coefficient of determination of the model.

[0258] S64, constructing dataset 2 with variables with VIF < 5 as independent variables and constructing dataset 2 with the dependent variable as the ecological benefit evaluation results, introducing the GWR model for regression analysis, using the Gaussian function as the weight function of the model, and using the AIC method to select a bandwidth that can effectively balance the bias and variance and avoid overfitting;

[0259] Among them, the GWR model calculation formula is:

[0260]

[0261] Where Y i is the explained value of the ecological benefit of city i, (u i ,v i ) is the geographical coordinates of city i; X ik is the explanatory value of the independent variable of city i, is the influencing factor; β k (u i ,v i ) is the centroid of the unit i region (u i ,v i ) at the regression parameter; ε i is the random error term;

[0262] The calculation formula for determining the bandwidth using the AIC method is:

[0263] AIC = nlog e (RSS)+2p

[0264] Among them, n represents the number of sample points, RSS represents the residual sum of squares, and p represents the number of unknown parameters. Based on this, the AIC calculation method for fusion with the GWR model is as follows:

[0265]

[0266] in, Maximum likelihood variance estimation for the GWR model:

[0267]

[0268] S65, based on S63, the factor of 5≤VIF<10 is used as the independent variable. The dependent variable is the combination of the ecological benefit evaluation results of S4 and the residual of the GWR model in S64. A new data set 3 is constructed and imported into the XGBoost model for training. The Bayesian optimization method is used to adjust the hyperparameters of the XGBoost model;

[0269] First, the data set is standardized. To avoid overfitting of the model, the data set is regularized, and 80% of the data is used as a training set and 20% as a test set for model verification.

[0270] Secondly, the Bayesian optimization method is used to adjust the XGBoost model hyperparameters, and the results are shown in Table 10;

[0271] The fmin function is used to find the optimal parameter combination from the defined parameter space space_xgb to minimize the return value of the objective function objective. The objective function defines the evaluation indicators of the model, such as accuracy or loss value;

[0272] Secondly, the tree-based Parzen estimator is specified as the search algorithm. max_evals = 50 indicates the maximum number of evaluations, that is, the maximum number of times Hyperopt will try different hyperparameter combinations. The best hyperparameter combination is obtained by combining the previously defined parameter space space_xgb and the result returned by the fmin function;

[0273] Table 10 Model accuracy verification

[0274]

[0275] S66, verify the accuracy of the model, select mean absolute error, root mean square error and determination coefficient as the evaluation indicators of prediction effect, and use ten-fold cross validation technology to evaluate the generalization performance of the model;

[0276] S67, the positive and negative influence of each feature is calculated using the SHAP model in Python. The results are as follows Fig.15 ;

[0277] Among them, the environment settings when using the SHAP library in Python to interpret the model, and ensure that the drawn graphics support Chinese display and the correct display of positive and negative signs;

[0278] S68, identify key influencing factors, select R in the GWR model fitting results. 2 The regression results with high values ​​are shown in Table 11. The results show that among the three natural geographical factors, slope has the greatest impact on ecological benefits, followed by elevation, while aspect has a smaller impact;

[0279] S69, based on the results of S67, it can be seen that the urban blue-green spatial pattern factors with high influence after the important features are sorted are AREA, PLAND, SPLIT, CONTIG, ENN, FRAC, CIRCLE, PLADJ, and ED. In order to further avoid the overfitting phenomenon of the model, the above features and the corresponding ecological benefits are constructed into a new data set, and the model is trained and interpreted to obtain Fig.16 The top five important features of the blue-green spatial pattern patch layer are AREA, CONTIG, ENN, FRAC, and CIRCLE, among which the highest ranking is the area index of blue-green patches; the top four important features of the blue-green spatial pattern type layer are PLAND, SPLIT, ED, and PLADJ, among which the highest ranking is the proportion of blue-green space area.

[0280] Table 11 GWR model fitting results

[0281]

[0282] It should be noted that the above content merely illustrates the technical idea of ​​the present invention and cannot be used to limit the protection scope of the present invention. For ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications all fall within the protection scope of the claims of the present invention.

Claims

1. A method for identifying key influencing factors of the ecological benefits of urban blue-green space, characterized in that: The following steps are involved: S1, construct an evaluation index set of ecological benefits of urban blue-green space; S2, based on the indicator set constructed in S1, calculates and standardizes the indicator data in the target area to establish a standardized data matrix of the ecological benefits of urban blue-green space; S3, based on the standardized data matrix of ecological benefits in S2, uses the analytic hierarchy process to determine the subjective weights of indicators, uses the entropy weight method to determine the objective weights, and obtains the weights of various indicators through combined weighting; S4, construct the entropy weight-TOPSIS model to calculate the weighted standardized data matrix of the ecological benefits of urban blue-green space, thereby completing the ecological benefit evaluation of urban blue-green space; S5, establish a set of influencing factors of the ecological benefits of urban blue-green space, including two categories: natural factors and spatial pattern factors; S6, based on the evaluation results of the ecological benefits of urban blue-green space in S4 and the influencing factor set in S5, combined with the GWR model and XGBoost model in machine learning to identify key influencing factors.

2. A method for identifying key influencing factors of the ecological benefits of urban blue-green space according to claim 1, characterized in that: In S1, the target layer of the urban blue-green space ecological benefit evaluation index set is ecological benefit; the criterion layer includes stormwater resilience, cooling benefit, carbon sink benefit, and biodiversity; the indicator layer includes surface runoff, flooding risk, surface temperature, net primary productivity, carbon storage, and habitat quality.

3. The method for identifying key influencing factors of the ecological benefits of urban blue-green space according to claim 1 is characterized in that: The steps for establishing a standardized data matrix of the ecological benefits of urban blue-green spaces in S2 include: S21, obtaining remote sensing image data, land use data, normalized difference vegetation index NDVI data, vegetation type data, temperature data, precipitation data and daily radiation data, urban soil type data, and urban land cover data in the target area; S22, based on the land use data, urban soil type data, urban surface cover data, and precipitation data obtained in S21, the surface runoff of the target area is calculated according to the SCS-CN runoff generation model; S23, based on the remote sensing images of continuous sunny days and after rain in S21 remote sensing image data, the density of waterlogging points is calculated to characterize the flooding risk in the study area; S24, based on the spectral data in the remote sensing image data of S21, the surface temperature of the target area is inverted using the atmospheric correction method; S25, using the land use data, normalized difference vegetation index NDVI data, vegetation type data, temperature data, precipitation data and daily radiation data obtained in S21, and using the CASA model, calculate the net primary productivity of the target area; S26, using the carbon storage module of the InVEST model to generate the spatial distribution of carbon storage; S27, the habitat quality module in the InVEST model was used for analysis to assess the habitat quality in the target area.

4. A method for identifying key influencing factors of the ecological benefits of urban blue-green space according to claim 3, characterized in that: In S22, the calculation steps of surface runoff in the target area are: S221, calculate the maximum possible water storage capacity based on the CN values ​​of urban and residential land, and assign values ​​to soils of different urban land types in the target area in combination with land use data; The soil cover composite number calculation formula is as follows: In the formula, S is the maximum possible water storage capacity, which indicates the maximum potential water storage capacity of the soil; S222, combined with S21 precipitation data, calculates the surface runoff in the target area according to the SCS-CN runoff model calculation formula; the SCS-CN runoff model calculation formula is as follows: In the formula, Q represents the surface runoff, in millimeters; P represents the rainfall, also in millimeters; λ represents the initial loss rate, which is a dimensionless value, usually calculated using λ=0.2; Secondly, the calculation steps of flooding risk of target area in S23 are as follows: S231, based on the remote sensing image data of S21, remote sensing images under the scenarios of continuous sunny days and after rain are selected to classify land use; S232, selecting the water body classification data of the land use classification results under the continuous sunny day and after rain scenarios in S231, removing the permanent water area in the target area, obtaining the seasonal water body, and obtaining the waterlogged area after subtraction; S233, after converting the waterlogged area into points, the point density was calculated using the "point density analysis" tool in the GIS software to obtain the waterlogged point density, which was used to characterize the flooding risk of the waterlogged area; Again, in S24, the steps for obtaining the surface temperature are: S241, through the spectral data containing surface temperature information in the S21 remote sensing image data, the "FLAASH Atmospheric Correction" tool of ENVI software is used to perform atmospheric correction processing on the remote sensing image data based on the spectral band; S242, applying a temperature inversion algorithm to convert the corrected remote sensing data into surface temperature, and using Planck's formula to calculate the true surface temperature; Among them, the expression of the thermal infrared radiation brightness value received by the satellite sensor can be written as: L λ =[ε·B(T)+(1-ε)L↓]·τ+L↑ Among them, the upward radiation brightness of the atmosphere L↑, ε is the surface emissivity, T is the true surface temperature, B(T) is the thermal radiation brightness of the black body at T derived from Planck's law, τ is the transmittance of the atmosphere in the thermal infrared band, then the radiation brightness B(T) of the black body with temperature T in the thermal infrared band is: Secondly, the Planck calculation formula is as follows: Where K1 and K2 are Landsat sensor specific calibration constants; Then, in S25, the steps to calculate the net primary productivity are: S251, in GIS software, the monthly precipitation, temperature and solar radiation data of the target year are subjected to Kriging interpolation processing, and band synthesis is performed to obtain an annual average grid map; then, the grid map of precipitation, temperature and solar radiation data is obtained by clipping according to the target area range; S252, using the CASA model, input precipitation, temperature and solar radiation data, normalized vegetation index NDVI map, vegetation data, set the parameters in the static file, including the maximum and minimum values ​​of the normalized vegetation index NDVI data, the maximum light energy utilization rate, and the vegetation index SR, and obtain the distribution of net primary productivity in the target area; Secondly, in S26, the steps for obtaining the spatial distribution of carbon storage are: S261, based on the S21 land use data, uses the carbon storage module of the InVEST model to calculate the total carbon storage of each land type based on the different land types and corresponding carbon pool data in the target area; S262, classifying and summarizing the calculated carbon storage data by land type to generate a spatial distribution map of carbon storage; Among them, the formula for calculating carbon storage through the InVEST model is divided into: C t =C a +C b +C s +C d Where: C t is the total carbon storage of the site (t / hm 2 );C a is the aboveground carbon storage (t / hm 2 );C b is the underground carbon storage (t / hm 2 );C s is soil carbon storage (t / hm 2 );C d is the dead organic carbon stock (t / hm 2 ); Finally, in S27, the steps for assessing habitat quality in the target area are: S271, using the habitat quality module in the InVEST model for analysis, and assessing biodiversity based on the quality of the habitat obtained; The formula for calculating habitat quality score is: Q xj =H j (1-D xj ) Where Q xj is the habitat quality score of grid x in land use type j, H j is the habitat suitability score of land use type j, D xj is the threat level of grid x; S272, threat source data refers to the factors that threaten the ecological land in the land use type. GIS software is used to assign a value of 1 to the threat area and a value of 0 to the non-threat area; The calculation formula of grid threat level is: In the formula, r y is the threat intensity of grid y, ω r is the relative weight of threat source r, β x is the protection level of grid x, S jr is the sensitivity of land use type j to threat source r, i rxy is the impact level of grid y in threat source r on grid x; Among them, the calculation formula of the impact level of grid y in threat source r on grid x is: Where, d xy is the Euclidean distance between grid x and grid y, d rmax is the maximum impact distance of the threat source r.

5. The method for identifying key influencing factors of the ecological benefits of urban blue-green space according to claim 1 is characterized in that: In S3, the steps to obtain the weights of each indicator are: S31, establish a number of sampling points, randomly extract the calculation results of the standardized data matrix of the ecological benefits of urban blue-green space in S2, and obtain a data set; S32, using the entropy weight method to calculate the entropy value of each indicator in the data set, and then obtain the objective weight of each indicator; wherein, the specific steps of using the entropy weight method to calculate the entropy value of each indicator are: first, perform data standardization, and establish the original information matrix X = (x ij )m×n, normalization is used to perform data standardization, and its expression is: Where Y xj are the standardized indexes of the jth coordinate of the ith evaluation index; x ij is the original data of the jth coordinate of the ith evaluation index; Secondly, determining the feature weight is a key step in the evaluation process; in order to calculate the weight of the jth indicator of the i-th evaluation object, we first need to define the weight p of the indicator. ij , through Y ij Perform calculations; Then, using the weight p of each indicator ij , calculate the information entropy, calculate the information entropy E of each evaluation index j ; In the formula, m is the number of evaluation objects, E j is the information entropy of each evaluation index; Finally, through the information entropy E j , use the entropy weight method to determine the indicator weight and calculate the objective weight β of the evaluation indicator j =β1,β2,…,β n , where β j is the weight corresponding to the jth indicator; In the formula, n is the number of evaluation indicators, β j is the entropy weight of the entropy weight method; S33, in the AHP hierarchical analysis, the four ecological goals of stormwater resilience, cooling benefits, carbon sink benefits, and biodiversity are set to be equal weights; the final weights of each indicator are obtained by comprehensively considering the calculation results of the AHP hierarchical analysis method and the entropy weight method; among them, the final weight formula of each indicator is obtained as follows: In the formula, α j is the weight of the AHP hierarchical analysis method corresponding to the jth indicator, β j is the weight of the entropy weight method corresponding to the jth indicator.

6. The method for identifying key influencing factors of the ecological benefits of urban blue-green space according to claim 1 is characterized in that: In S4, the specific steps of constructing the entropy weight-TOPSIS model to calculate the weighted standardized data matrix of the ecological benefits of urban blue-green space are as follows: S41, write a program in Matlab software, randomly select a number of sampling points, combine the indicator weight results of S3 with the calculation results of the standardized data matrix of the ecological benefits of urban blue-green space in S2, and obtain the positive and negative ideal solutions and the evaluation results of the ecological benefits corresponding to the sampling points according to the TOPSIS evaluation method; Among them, the calculation steps of the TOPSIS evaluation method are: First, a weighted normalized decision matrix is ​​constructed based on the factor weights determined by S3. The matrix formula is: V=X ij ·ω=[v ij ] m·n Where V is the weight matrix, X ij is the i-th evaluation index; Secondly, using the calculated weighting matrix V, determine the positive and negative ideal solutions: Then, according to the positive and negative ideal solution V + 、V - , use Euclidean distance to calculate the distance between positive and negative ideal solutions: In the formula, Positive ideal solution distance, is the negative ideal solution distance; Finally, the calculated positive and negative ideal solution distances are used to calculate the closeness: According to the formula, calculate the closeness T to the ideal solution j , T j The value range of T is 0 to 1; j The closer it is to 1, the better the ecological benefit; S42, select the indicator weights whose evaluation results are positive ideal solutions for final calculation to obtain the final weight results; then combine the standardized data matrix of the ecological benefits of the urban blue-green space in S2 to calculate the ecological benefits of the urban blue-green space to obtain the ecological benefit evaluation results.

7. The method for identifying key influencing factors of the ecological benefits of urban blue-green space according to claim 1 is characterized in that: The influencing factors of the ecological benefits of urban blue-green space in S5 include two categories: "spatial pattern factors" and "natural factors"; the spatial pattern factors are composed of patch layer factors and type layer factors; among them, patch layer factors include area, perimeter, fractal dimension, near-circle index, proximity index, and patch distance; type layer factors include area proportion, maximum patch index, edge density, landscape shape index, similarity adjacency ratio, aggregation, separation amount, cohesion index, and fragmentation; natural factors are composed of elevation, slope, and slope aspect.

8. The method for identifying key influencing factors of the ecological benefits of urban blue-green space according to claim 1 is characterized in that: The specific steps in S6 to identify key influencing factors by combining the GWR model and the XGBoost model in machine learning are: S61, constructing data sets based on the natural factors and pattern factors of urban blue-green space in S5 and the corresponding ecological benefit evaluation results calculated in S4; using random sampling tools in GIS software, creating random sampling points, extracting various influencing factors and ecological benefit evaluation results to the sampling points, and eliminating points that do not belong to urban blue-green space to form data set 1; S62, the spatial autocorrelation was calculated using the global Moran's index tool; The calculation formula of Moran Index is as follows: In the formula, n is the total number of spatial units, y i and j Respectively represent the attribute values ​​of the i-th spatial unit and the j-th spatial unit, is the mean of all spatial unit attributes, w ij is the spatial weight value; S63, writing a VIF calculation program in Python software to calculate the multicollinearity of each influencing factor of data set 1 in S61; The formula for calculating VIF is: In the formula, R 2 is the coefficient of determination of the model; S64, constructing dataset 2 with variables with VIF < 5 as independent variables and constructing dataset 2 with the dependent variable as the ecological benefit evaluation results, introducing the GWR model for regression analysis, using the Gaussian function as the weight function of the model, and using the AIC method to select a bandwidth that can effectively balance the bias and variance and avoid overfitting; Among them, the GWR model calculation formula is: Where Y i is the explained value of the ecological benefit of city i, (u i ,v i ) is the geographical coordinates of city i; X ik is the explanatory value of the independent variable of city i, is the influencing factor; β k (u i ,v i ) Unit i area centroid (u i ,v i ) at the regression parameter; ε i is the random error term; The calculation formula for determining the bandwidth using the AIC method is: AIC = log e (RSS)+2p Among them, n represents the number of sample points, RSS represents the residual sum of squares, and p represents the number of unknown parameters. On this basis, the AIC calculation method for fusion with the GWR model is as follows: in, Maximum likelihood variance estimation for the GWR model: S65, based on S63, takes the factor of 5≤VIF<10 as the independent variable; the dependent variable is the combination of the ecological benefit evaluation result of S4 and the residual of the GWR model in S64, constructs a new data set 3, imports the XGBoost model for training; uses the Bayesian optimization method to adjust the hyperparameters of the XGBoost model; S66, verify the accuracy of the model, select mean absolute error, root mean square error and determination coefficient as the evaluation indicators of prediction effect, and use ten-fold cross validation technology to evaluate the generalization performance of the model; S67, the positive and negative influence of each feature was calculated using the SHAP model in Python software; S68, in the GWR model fitting results, select R 2 The regression results with high values ​​are taken as key influencing factors; S69, in order to further avoid the overfitting of the model, the urban blue-green spatial pattern factors with high influence and the corresponding ecological benefits after the important features in S67 are ranked are constructed into a new data set, which is interpreted after model training to obtain the importance ranking of the influencing factors.

Images