Method for identifying key influencing factors on ecological benefits of urban blue-green spaces

Abstract

Disclosed in the present invention is a method for identifying key influencing factors on ecological benefits of urban blue-green spaces. The method comprises: first, constructing an evaluation indicator set for ecological benefits of urban blue-green spaces; then, calculating data of indicators and performing standardization processing, and establishing a standardized data matrix for the ecological benefits of the urban blue-green spaces; next, using an analytic hierarchy process to determine subjective weights of the indicators, using an entropy weight method to determine objective weights, and obtaining weights of the indicators by means of combination weighting; afterwards, constructing an entropy weighted TOPSIS model to calculate a weighted standardized data matrix for the ecological benefits of the urban blue-green spaces, so as to complete the evaluation of the ecological benefits of the urban blue-green spaces; subsequently, establishing an influencing factor set for the ecological benefits of the urban blue-green spaces; and finally, on the basis of evaluation results of the ecological benefits of the urban blue-green spaces, using a GWR-XGBoost model to identify key influencing factors. In the present invention, a multi-criteria decision analysis method, a geographically weighted regression model and a machine learning technique are combined, thereby enabling the rapid and accurate acquisition of key influencing factors on ecological benefits of urban blue-green spaces.

Description

A method for identifying key influencing factors of ecological benefits of urban blue-green space TECHNICAL FIELD

[0001] The present application belongs to the field of urban ecological environment management, and particularly relates to a method for identifying key influencing factors of ecological benefits of urban blue-green space. BACKGROUND

[0002] Urban blue-green space is composed of blue space and green space in the city, and is the ecological background of the city and an important content of the current high-quality development of the city. Treating the blue-green system of the city as spatial infrastructure to coordinate the development of ecology, economy, society and other aspects has become a widely accepted consensus of global cities in response to climate change and rapid urbanization challenges. 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 limited urban space, how to effectively use limited blue-green space to provide better ecological benefits has become a problem to be solved.

[0003] However, the existing evaluation of ecological benefits focuses on large-scale space, and the indicators and methods are not suitable for urban and sub-regional scales. The research on urban blue-green space is still in its infancy. The ecological benefit evaluation methods constructed in previous studies focus on a single land use type or a single target, lack comprehensiveness and overall planning, and the evaluation methods are not suitable for multi-target blue-green space ecological benefit evaluation and cannot accurately identify the multi-target factors affecting ecological benefits. Therefore, the present application aims to provide a method for identifying key influencing factors of ecological benefits of urban blue-green space, which can quickly and accurately obtain the key influencing factors of ecological benefits of urban blue-green space on the basis of comprehensive evaluation of the ecological benefits of urban blue-green space.

[0004] SUMMARY

[0005] To solve the above problems, the present application discloses a method for identifying key influencing factors of ecological benefits of urban blue-green space, which combines multi-criteria decision analysis method, geographic weighted regression model and machine learning technology, and can quickly and accurately obtain the key influencing factors of ecological benefits of urban blue-green space, which is helpful for the strategy making of blue-green space optimization and reorganization.

[0006] To achieve the above purpose, the technical scheme of the present application is as follows:

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

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

[0009] S2, based on the index set constructed in S1, calculate the index data in the target area and perform standardization processing to establish a standardized data matrix of urban blue-green space ecological benefits;

[0010] S3, based on the standardized data matrix of ecological benefits in S2, determine the subjective weight of the index by using the analytic hierarchy process, determine the objective weight by using the entropy weight method, and obtain the weight of each index by combining the weights;

[0011] S4, construct an entropy weight-TOPSIS model to calculate the weighted standardized data matrix of urban blue-green space ecological benefits, so as to complete the evaluation of urban blue-green space ecological benefits;

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

[0013] S6, based on the evaluation results of urban blue-green space ecological benefits in S4 and the set of influence factors in S5, combine the GWR model (Geographically weighted regression, Geographically weighted regression model) and the XGBoost model (Extreme Gradient Boosting) in machine learning to identify key influence factors.

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

[0015] Further, in S2, the steps for establishing the standardized data matrix of urban blue-green space ecological benefits include:

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

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

[0018] S23, according to the remote sensing images in S21 under the continuous sunny and post-rain scenarios, calculate the density of waterlogging points to represent the flooding risk in the study area;

[0019] S24, based on the spectral data in the remote sensing image data of S21, the atmospheric correction method is used to obtain the surface temperature of the target region;

[0020] S25, by using the CASA (Carnegie-Ames-Stanford Approach) model, the land use data, normalized vegetation index NDVI data, vegetation type data, air temperature data, precipitation data and daily radiation data obtained by S21 are used to calculate the net primary productivity of the target region;

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

[0022] S27, the Habitat Quality module in the InVEST model is used for analysis to evaluate the habitat quality of the target region.

[0023] Further, in S2, the specific calculation steps of the target region urban blue-green space ecological benefit standardized data matrix are:

[0024] S221, according to the CN value (Curve Numbers) of urban and residential land, the maximum possible water storage capacity is calculated, and the soil of different urban land types in the target region is valued in combination with the land use data;

[0025] Wherein, the soil cover composite number calculation formula is as follows:

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

[0027] S222, in combination with the precipitation data of S21, the surface runoff of the target region is calculated according to the SCS-CN runoff model calculation formula;

[0028] Wherein, the SCS-CN runoff model calculation formula is as follows:

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

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

[0031] S231, according to the remote sensing image data in S21, remote sensing images in continuous sunny days and after rain scenarios are selected respectively, and land use classification is carried out;

[0032] S232, water body classification data of land use classification results in continuous sunny days and after rain scenarios in S231 are selected, and the common water area in the target area is removed, so as to obtain seasonal water body, and the difference is obtained after the subtraction, so as to obtain the waterlogging area;

[0033] S233, after the waterlogging area is converted into points by using the "point density analysis" tool in the GIS software, the point density is calculated, so as to obtain the waterlogging point density as the flooding risk of the waterlogging area;

[0034] Thirdly, in S24, the step of obtaining the surface temperature is:

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

[0036] S242, the temperature inversion algorithm is applied to convert the corrected remote sensing data into the surface temperature, and the Planck formula is used to calculate the real surface temperature;

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

[0038] L λ =[epsilon cdot B(T) + (1 - epsilon) L downarrow] cdot tau + L uparrow

[0039] Wherein, the atmospheric upward radiation brightness L uparrow, epsilon is the surface emissivity, T is the real surface temperature, B(T) is the thermal radiation brightness of the black body at T derived by the Planck law, tau is the transmittance of the atmosphere in the thermal infrared band, and the radiation brightness B(T) of the black body at the thermal infrared band with temperature T is:

[0040] The Planck calculation formula is as follows:

[0041] In the formula, K1 and K2 are specific calibration constants of the Landsat8 sensor.

[0042] Then, in S25, the step of calculating the net primary productivity is:

[0043] S251, in the GIS software, the precipitation, air temperature and solar radiation data of the target year are processed by Kriging interpolation, and band synthesis is carried out, so as to obtain the annual average grid picture; then, the grid pictures of the precipitation, air temperature and solar radiation data are obtained by cutting according to the range of the target area;

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

[0045] Secondly, in S26, the step of obtaining the spatial distribution of carbon storage is:

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

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

[0048] Wherein, the formula for calculating carbon storage by InVEST model is divided into:

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

[0050] In the formula: 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 the soil carbon storage (t / hm 2 ); C d is the dead organic carbon storage (t / hm 2 );

[0051] Finally, in S27, the step of evaluating the habitat quality of the target area is:

[0052] S271, using the habitat quality module in the InVEST model for analysis, and evaluating the biological diversity according to the advantages and disadvantages of the obtained habitat quality;

[0053] Wherein, the formula for calculating habitat quality score is:

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

[0055] In the formula, 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;

[0056] S272, the model assumes that the area with good habitat quality has high biodiversity; the threat source data is a factor that threatens ecological land in the land use type, and the GIS software is used to assign a value of 1 to the threatened area and a value of 0 to the non-threatened area;

[0057] wherein the formula for calculating the threat level of the grid is:

[0058] In the formula, r y is the threat intensity of grid y, ω r is the relative weight of the 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 influence level of grid y in threat source r on grid x;

[0059] wherein the formula for calculating the influence level of grid y in threat source r on grid x is:

[0060] In the formula, d xy is the Euclidean distance between grid x and grid y, d rmax is the maximum influence distance of threat source r.

[0061] Further, in S3, the steps for obtaining the weights of each index are:

[0062] S31, a plurality of sampling points are established, and the calculation results of the urban blue-green space ecological benefit standardization data matrix in S2 are randomly extracted to obtain a data set.

[0063] S32, the entropy weight method is used to calculate the entropy values of each index, and the objective weights of each index are obtained.

[0064] Wherein, the specific steps for calculating the entropy values of each index using the entropy weight method are as follows: first, data standardization is performed, and under the condition that n evaluation indexes and m evaluation objects are given, the original information matrix X=(x ij )mxn is established, and data standardization is performed by using normalization, and the expression is:

[0065] In the formula, Y ij are the standardized indexes of the i-th evaluation index and the j-th coordinate, respectively; x ijThe original data of the jth coordinate of the ith evaluation index.

[0066] Secondly, determining the feature weight is a key link in the evaluation process. In order to calculate the proportion of the jth index of the ith evaluation object, the weight p ij of the index is first defined. ij The Y ij in S321 is used for calculation.

[0067] Then, the information entropy E j of each evaluation index is calculated using the index weight p ij .

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

[0069] Finally, the index weight is determined by the entropy weight method using the information entropy E j , and the objective weight β j of the evaluation index is calculated, where β n 1, β2, …, β j n, and β j j is the weight corresponding to the jth index.

[0070] In the formula, n is the number of evaluation indexes, and β j is the entropy weight value of the entropy weight method.

[0071] In S33, the rainwater resilience, cooling benefit, carbon sink benefit, and biodiversity are set as equal weights in the AHP hierarchical analysis. By comprehensively considering the calculation results of the AHP hierarchical analysis method and the entropy weight method, the final weight of each index is obtained. The formula for obtaining the final weight of each index is:

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

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

[0074] In S41, a program is written in Matlab software, a number of sampling points are randomly selected, the index weight result of S3 is combined with the calculation result of the standardized data matrix of the urban blue-green space ecological benefit of S2, and the ecological benefit of each point is calculated. According to the TOPSIS evaluation method, the positive and negative ideal solutions are obtained, and the evaluation result of the ecological benefit corresponding to the sampling points is obtained.

[0075] The calculation steps of the TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution) evaluation method are as follows:

[0076] First, in order to comprehensively evaluate each evaluation object, a weighted normalized decision matrix is constructed according to the factor weights determined in S3, and the matrix formula is as follows:

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

[0078] In the formula, V is the weighted matrix, X ij is the yth evaluation index.

[0079] Secondly, the positive and negative ideal solutions are determined by using the calculated weighted matrix V:

[0080] Then, according to the positive and negative ideal solutions V + , V - , the Euclidean distance between the positive and negative ideal solutions is calculated:

[0081] In the formula, is the distance of the positive ideal solution, is the distance of the negative ideal solution.

[0082] Finally, the calculated distances of the positive and negative ideal solutions are calculated to obtain the closeness degree:

[0083] According to the formula, the closeness degree T j to the ideal solution is calculated, and the value range of T j is 0-1. In the present application, the closer T j is to 1, the better the ecological benefit is.

[0084] S42, the index weight of the evaluation result being the positive ideal solution is selected for final calculation to obtain the final weight result; and then the urban blue-green space ecological benefit is calculated in combination with the urban blue-green space ecological benefit standardization data matrix in S2 to obtain the ecological benefit evaluation result.

[0085] Further, the urban blue-green space ecological benefit influencing factors in S5 include two categories of "spatial pattern factors" and "natural factors". The spatial pattern factors are divided into patch layer factors and type layer factors. The patch layer factors include area, perimeter, fractal dimension, near-circular index, proximity index, and patch distance. The type layer factors include area proportion, largest patch index, edge density, landscape shape index, similar adjacency ratio, aggregation degree, separation degree, cohesion index, and fragmentation. The natural factors are composed of elevation, slope, and aspect.

[0086] Further, in S6, the specific steps of identifying key influencing factors by combining the GWR model and the XGBoost model in machine learning are as follows:

[0087] S61, based on the urban blue-green space natural factors and pattern factors in S5 and the corresponding ecological benefit evaluation results calculated in S4, a data set is constructed. In the GIS software, a random sampling tool is used to create random sampling points, and each influencing factor and the ecological benefit evaluation result are extracted to the sampling points. Points not belonging to the urban blue-green space are removed to form data set 1.

[0088] S62, the global Moran's I tool is used to calculate the spatial autocorrelation.

[0089] The calculation formula of Moran's I is as follows:

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

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

[0092] The formula for calculating VIF is as follows:

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

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

[0095] The GWR model calculation formula is:

[0096] In the formula, Y i is the explanatory value of the ecological benefit of city i, (u i ,v i ) is the geographical coordinate of city i; X ik is the explanatory value of the independent variable of city i, and is the influence factor; β k (u i ,v i ) is the regression parameter of the centroid (u i ,v i ) of unit i; and ε i is a random error term.

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

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

[0099] where n represents the number of sample points, RSS represents the residual sum of squares, and p is the number of unknown parameters. On this basis, the AIC calculation method for fusion with the GWR model is as follows:

[0100] where, is the maximum likelihood variance estimation of the GWR model:

[0101] S65, based on S63, the factors with 5≤VIF<10 are taken as the independent variables. The dependent variable is the combination of the ecological benefit evaluation results of S4 and the residuals of the GWR model in S64, and a new data set 3 is constructed and imported into the XGBoost model for training. The Bayesian optimization method is used to adjust the XGBoost model hyperparameters.

[0102] S66, model accuracy verification is performed, the mean absolute error, root mean square error and coefficient of determination are selected as the evaluation indexes of the prediction effect, and the ten-fold cross-validation technology is used to evaluate the generalization performance of the model.

[0103] S67, the SHAP model in Python software is used to calculate the positive and negative influence degree of each feature.

[0104] S68, in the fitting results of the GWR model, the regression results with high R 2 value are selected as the key influence factors.

[0105] S69, in order to further avoid the overfitting phenomenon of the model, the city blue-green space pattern factor with high influence degree after the important feature ranking in S67 is combined with the corresponding ecological benefit to build a new data set, and the model is trained and interpreted to obtain the importance ranking of the influence factors.

[0106] The present application has the following advantages:

[0107] The city blue-green space ecological benefit key influence factor identification method provided by the present application is suitable for ecological benefit evaluation and influence factor identification at the city scale. The AHP method, entropy weight method and TOPSIS method are used for comprehensive evaluation to obtain the weight of each ecological benefit index, which is more objective than the traditional method for determining the weight, and avoids the interference of subjective factors. The advanced method of machine learning is introduced to accurately and quickly identify the key factors affecting the ecological benefit of the city blue-green space, which is helpful for the strategy making of the optimization and reorganization of the blue-green space. BRIEF DESCRIPTION OF DRAWINGS

[0108] Fig. 1 is a flow chart of the city blue-green space ecological benefit index evaluation system construction method of the present application;

[0109] Fig. 2 is a specific remote sensing data area graph of the embodiment of the present application;

[0110] Fig. 3 is a surface runoff graph of the embodiment of the present application;

[0111] Fig. 4 is a submergence risk graph of the embodiment of the present application;

[0112] Fig. 5 is a ground surface temperature inversion graph of the embodiment of the present application;

[0113] Fig. 6 is a carbon storage amount graph of the embodiment of the present application;

[0114] Fig. 7 is a carbon storage amount graph of the embodiment of the present application;

[0115] Fig. 8 is a biological diversity graph of the embodiment of the present application;

[0116] Fig. 9 is an ecological benefit evaluation result graph of the embodiment of the present application;

[0117] Fig. 10 is an elevation, slope and slope direction graph of the embodiment of the present application;

[0118] Fig. 11 is a patch layer index calculation graph of the embodiment of the present application;

[0119] Fig. 12 is a type layer index calculation graph of the embodiment of the present application;

[0120] Fig. 13 is a GWR-XGBoost model construction flow chart of the embodiment of the present application;

[0121] Fig. 14 is a spatial autocorrelation analysis result graph of the embodiment of the present application;

[0122] Figure 15 is a plot of the important features of the spatial pattern factor in the embodiment of the application;

[0123] Figure 16 is a plot of the important features of the key spatial pattern factor in the embodiment of the application. DETAILED DESCRIPTION

[0124] The application will be further clarified by the following description and specific embodiments, which should be understood as merely illustrating the application and not limiting the scope of the application.

[0125] As shown in the figure, the specific steps of the method for identifying key influencing factors of urban blue-green space ecological benefits according to the application include:

[0126] S1, construct an evaluation index set of urban blue-green space ecological benefits, with ecological benefits as the target layer, rainwater resilience, cooling benefit, carbon sink benefit, and biodiversity as the criterion layer, and surface runoff, flooding risk, surface temperature, carbon storage amount, carbon storage amount, and habitat quality as the index layer. As shown in Table 1.

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

[0128] S2, select Nanjing as the case area, calculate the index data in the target area based on the index set constructed in S1, and perform standardization processing to establish a standardized data matrix of urban blue-green space ecological benefits, and obtain the surface runoff, rainwater flooding range, surface temperature inversion graph, net primary productivity spatial distribution graph, carbon storage spatial distribution graph, and biodiversity distribution graph of the study area. The specific steps include:

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

[0130] Table 2 Data sources

[0131] First, land use data: based on the Landsat8 (2020) remote sensing image of Nanjing city downloaded from the geographic spatial cloud platform, supervised classification and manual visual interpretation correction are performed on the ENVI software platform. After accuracy verification, the land use data in the study area has been obtained with an accuracy of 30m x 30m;

[0132] Then, air temperature, precipitation, and daily radiation data: obtain the data of 13 weather stations in Nanjing from the geographic remote sensing ecological network;

[0133] Secondly, vegetation type data: obtain the vegetation type data in Nanjing from the geographic remote sensing ecological network, with a precision of 30m x 30m;

[0134] Similarly, normalized vegetation index NDVI data: in the Google Earth Engine, based on the Landsat8 remote sensing data, after the correction and cloud removal preprocessing of the remote sensing image, according to the calculation formula of the normalized vegetation index NDVI, the normalized vegetation index NDVI data in Nanjing is obtained, with a precision of 30m x 30m;

[0135] The calculation formula of the normalized vegetation index NDVI is as follows:

[0136] In the formula, NIR is the reflectivity of the near-infrared band, and RED is the reflectivity of the red light band;

[0137] Finally, soil data: in the Google Earth Engine, download the soil classification data under the classification system of the United States Department of Agriculture, with a precision of 30m x 30m.

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

[0139] Table 3 soil data

[0140] Firstly, according to the CN value of the city and residential land, the maximum possible water storage capacity is calculated, and the soil of different city land types in the target area is valued in combination with the land use data. The soil cover composite number calculation formula is as follows:

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

[0142] Then, based on the CN value, the soil of different city land types in the study area is valued through the "raster calculator" tool in the GIS software;

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

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

[0145] S23, according to S31 remote sensing image data in the continuous sunny and after rain scene, calculate the density of waterlogging points to represent the flooding risk in the study area, the results as shown in Figure 4, the specific process is:

[0146] First, according to S21 remote sensing image data, respectively select the continuous sunny and after rain scene, land use classification is carried out on remote sensing image;

[0147] Then, select S231 in the continuous sunny and after rain scene land use classification results of water classification data, remove the target area of the water area, get the seasonal water body, subtraction after getting the waterlogging area;

[0148] Finally, the waterlogging area is converted into points by using GIS tool, and the point density is calculated by using "point density analysis" tool, and the waterlogging point density is obtained, which represents the flooding risk of waterlogging area.

[0149] S24, based on the spectral data in S21 remote sensing image data, the atmospheric correction method is used to retrieve the surface temperature of the target area, the results as shown in Figure 5, the specific process is:

[0150] First, through the spectral data containing surface temperature information in S21 remote sensing image data, the "FLAASH Atmospheric Correction" tool of ENVI software platform is used to carry out atmospheric correction processing on remote sensing image data based on spectral band.

[0151] Second, the corrected remote sensing data is converted into surface temperature by using temperature inversion algorithm, and the surface temperature is calculated by using Planck formula. Then, the accuracy of the retrieved surface temperature data is verified by ground observation data, and the results are analyzed and applied;

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

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

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

[0155] The Plank calculation formula is as follows:

[0156] In the formula, K1 and K2 are calibration constants specific to the Landsat8 sensor.

[0157] S25, land use data, normalized vegetation index NDVI data, vegetation type data, air temperature data, precipitation data and daily radiation data obtained by S21 are used to calculate the net primary productivity of the target area by using the CASA model, and the result is as shown in Fig. 6.

[0158] Firstly, the meteorological data of S21 is processed by using the GIS software, and the Kriging interpolation processing is performed on the monthly precipitation, air temperature and solar radiation data of the target year;

[0159] Then, based on the above data, the 12-month grid map bands of the target year are combined into an annual average grid map, and it is resampled to a spatial resolution of 30m x 30m. The precipitation, air temperature and solar radiation data grid maps of the target year are cut according to the research range and made;

[0160] Finally, in the ENVI software platform, the CASA model is used to input the precipitation, air temperature and solar radiation data, normalized vegetation index NDVI data and vegetation data, and 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, to obtain the actual case of the net primary productivity spatial distribution map, as shown in Fig. 6.

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

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

[0163] In the formula: C t is the total carbon storage of the site (t / hm 2 ); C a is the carbon storage of the aboveground part (t / hm 2 ); Cb Carbon storage of belowground part (t / hm2); C 2 ) ; C s Carbon storage of soil (t / hm2); C 2 ) ; C d Carbon storage of dead organic matter (t / hm2); C 2 ).

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

[0165] S27, the habitat quality module in the InVEST model is used for analysis, and the advantages and disadvantages of the obtained habitat quality are used to evaluate the biodiversity. The model assumes that the area with good habitat quality has high biodiversity; the threat source data is the factor threatening the ecological land in the land use type, and the assignment is as shown in Table 5. The GIS software is used to assign the threat area as 1 and the non-threat area as 0, and the result is as shown in Fig. 8. The specific calculation steps are as follows:

[0166] Firstly, the habitat quality score formula is calculated as follows:

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

[0168] In the formula, 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, and D xj is the threat level of grid x.

[0169] Among them, the calculation formula of the threat level of the grid is as follows:

[0170] 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, and i rxy is the influence level of grid y in threat source r to grid x.

[0171] Among them, the calculation formula of the influence level of grid y in threat source r to grid x is as follows:

[0172] In the formula, d xy is the Euclidean distance between grid x and grid y, and d rmax is the maximum influence distance of threat source r.

[0173] Table 5 Threat factor for habitat quality evaluation

[0174] S3, based on the ecological benefit standardized data matrix of S2, uses the analytic hierarchy process (AHP) to determine the subjective weights of the indicators and the entropy weight method to determine the objective weights. The weights of each indicator are obtained through combined weighting.

[0175] S31. In GIS software, create 20,000 random points, and use the "Multi-value Extraction to Point" tool to extract the calculation results of the above 6 indicators to the points to obtain the dataset;

[0176] S32 uses Matlab software to write relevant programs, and calculates the entropy value of each indicator by using the entropy weight method, thereby obtaining the objective weight of each indicator.

[0177] S321, Data standardization is an important data processing step. Given n evaluation indicators and m evaluation objects, the original information matrix X = (x ij For m×n, the data is standardized using normalization, and its expression is:

[0178] In the formula, Y ij These are the standardized indices for the i-th evaluation index and the j-th coordinate, respectively; x ij This refers to the original data for the j-th coordinate of the i-th evaluation index;

[0179] S322, Determining feature weights is a crucial step in the evaluation process. To calculate the weight of the j-th indicator for the i-th evaluation object, it is first necessary to define the weight p of that indicator. ij Y obtained from S321 ij Perform calculations;

[0180] S323, using the weights p of each indicator derived in S322. ij Calculate the information entropy, and calculate the information entropy Ej for each evaluation index;

[0181] In the formula, m is the number of evaluation objects, E j The information entropy of each evaluation indicator;

[0182] S324, through the information entropy E of S323 j The entropy weight method is used to determine the index weights, and the objective weight β of the evaluation index is calculated. j =β1, β2, ...,β n , where β j The weight corresponding to the j-th indicator;

[0183] wherein n is the number of evaluation indexes, β j is the entropy weight value of the entropy weight method;

[0184] S33, in the AHP hierarchical analysis, rainwater resilience, cooling benefit, carbon sink benefit, and biodiversity are set as equal weights. Through comprehensive consideration of the calculation results of the AHP hierarchical analysis method and the entropy weight method, the final weight of each index is obtained. The formula for obtaining the final weight of each index is as follows: the results are shown in Table 6;

[0185] wherein the formula for the final weight is:

[0186] wherein a j is the weight of the jth index corresponding to the AHP hierarchical analysis method, β j is the weight of the jth index corresponding to the entropy weight method.

[0187] Table 6 Calculation results of evaluation index weight

[0188] S4, based on the quantitative evaluation results of S2 and the weight calculation results of S3, an entropy weight-TOPSIS model is constructed to calculate the weighted urban blue-green space ecological benefit standardized data matrix, so as to complete the evaluation of urban blue-green space ecological benefit, as shown in FIG. 9.

[0189] S41, a program is written in Matlab, a plurality of sampling points are randomly selected, and the index weight results of S3 and the calculation results of the urban blue-green space ecological benefit standardized data matrix of S2 are combined, as shown in Table 7, to calculate the ecological benefit of each point. According to the TOPSIS evaluation method, the positive and negative ideal solutions and the evaluation results of the ecological benefit corresponding to the sampling points are obtained;

[0190] Further, in order to comprehensively evaluate each evaluation object, a weighted normalized decision matrix is constructed according to the factor weight determined in S3, and the matrix formula is:

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

[0192] wherein V is the weighted matrix, X ij is the ith evaluation index;

[0193] Secondly, the positive and negative ideal solutions are determined by using the calculated weighted matrix V:

[0194] Then, according to the obtained positive and negative ideal solutions V + , V - , the Euclidean distance between the positive and negative ideal solutions is calculated:

[0195] wherein, positive ideal solution distance, negative ideal solution distance;

[0196] Finally, the closeness degree is calculated according to the calculated positive and negative ideal solution distances:

[0197] The closeness degree T with the ideal solution is calculated according to the formula j , T j The value range of T is 0-1. The closer T is to 1, the better the ecological benefit is; j

[0198] Table 7 calculation results of positive and negative ideal solutions

[0199] S42, the index weight of the evaluation result as the positive ideal solution is selected for final calculation to obtain the final weight result; the value of the above calculation result is associated with the value of each point in the geographic space, Kriging interpolation is adopted, and then the ecological benefit of the urban blue-green space is calculated in combination with the standardized data matrix of the urban blue-green space ecological benefit in S2 to obtain the ecological benefit evaluation result, as shown in FIG. 9.

[0200] S5, an urban blue-green space ecological benefit influence factor set is established, including natural factors and spatial pattern factors. The urban blue-green space ecological benefit influence factors include “spatial pattern factors” and “natural factors”. The spatial pattern factors are divided into patch layer factors and type layer factors. The patch layer factors include area, perimeter, fractal dimension, near-circle index, proximity index, and patch distance; the type layer factors include area ratio, maximum patch index, edge density, landscape shape index, similar adjacency ratio, aggregation degree, separation degree, cohesion index, and fragmentation. The natural factors include elevation, slope, and aspect, as shown in Table 8.

[0201] Table 8 urban blue-green space ecological benefit influence factor set

[0202] Further, the elevation data is downloaded from the geographic space data cloud, the slope and aspect data are obtained by analyzing and processing the slope and aspect tools in the GIS software based on the DEM data. The data processing result is shown in FIG. 10;

[0203] Further, the Fragstats software is used to quantitatively analyze the urban blue-green space pattern characteristics in the selected year. The patch layer calculation result is shown in FIG. 11, and the type layer calculation result is shown in FIG. 12.

[0204] ​S6, based on the evaluation results of the ecological benefits of urban blue-green space in S4 and the set of influencing factors in S5, identifies key influencing factors by combining the GWR model and XGBoost model in machine learning. The steps are shown in Figure 13.

[0205] The specific steps for identifying key influencing factors by combining the GWR model and XGBoost model in machine learning are as follows:

[0206] S61 is a dataset constructed based on the natural and pattern factors of urban blue-green space in S5 and the corresponding ecological benefit evaluation results calculated in S4. In GIS software, a random sampling tool is used to create random sampling points, and the influencing factors and ecological benefit evaluation results are extracted to the sampling points. Points that do not belong to urban blue-green space are removed to form dataset 1.

[0207] S62. In GIS software, the Global Moran's Index (GWR) tool was used to process the ecological benefit evaluation results and calculate their spatial autocorrelation, as shown in Figure 14. The evaluation criteria are shown in Table 9. The spatial autocorrelation Moran's Index of the ecological benefit evaluation results within the case target area is 0.78, P = 0, and the absolute value of the Z-score is much greater than 2.58. This indicates that the spatial distribution of ecological benefits in the case target area is not random and has a significant positive correlation with the spatial distribution. Therefore, the standard regression model is not applicable in this case, while the GWR model is more suitable.

[0208] The Moran index is calculated using the following formula:

[0209] In the formula, n is the total number of spatial units, y i and y j Let represent the attribute values ​​of the i-th and j-th spatial units, respectively. w is the mean of all spatial unit properties. ij This represents the spatial weight value;

[0210] Table 9. Spatial Autocorrelation Test Criteria

[0211] S63, Write a VIF calculation program in Python to calculate the multicollinearity of each influence factor in dataset 1 of S61;

[0212] The formula for calculating VIF is as follows:

[0213] In the formula, R 2 The coefficient of determination for this model is denoted as .

[0214] S64, with VIF < 5, the variable construction dataset 2 as the independent variable, the dependent variable is the ecological benefit evaluation results of the construction of data set 2, the introduction of GWR model regression analysis, using Gaussian function as the weight function of the model, with AIC method to select a bandwidth that can effectively balance the deviation and variance, avoid over fitting;

[0215] The GWR model calculation formula is:

[0216] In the formula, Y i is the explanatory 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, which is the influencing factor; β k (u i , v i ) is the regression parameter at the centroid (u i , v i ) of unit i; ε i is the random error term;

[0217] The calculation formula for bandwidth determination using AIC method is:

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

[0219] Where n represents the number of sample points, RSS represents the residual sum of squares, and p is the number of unknown parameters. On the basis of GWR model, the AIC calculation method is as follows:

[0220] , which is the maximum likelihood variance estimation of GWR model:

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

[0222] Firstly, the data set is standardized to avoid the over fitting phenomenon of the model, and the data set is further regularized. 80% of the data is taken as the training set and 20% as the test set for model verification;

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

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

[0225] Secondly, the tree-based Parzen estimator is specified as the search algorithm. max_evals = 50 represents the maximum number of evaluations, i.e. the maximum number of different hyperparameter combinations that Hyperopt will try. Combined with the previously defined parameter space space_xgb and the return value of the fmin function to obtain the best hyperparameter combination;

[0226] Table 10 Model Accuracy Verification

[0227] S66, model accuracy verification is performed, and mean absolute error, root mean square error and coefficient of determination are selected as evaluation indicators of prediction effect, and ten-fold cross-validation technology is used to evaluate the generalization performance of the model;

[0228] S67, the positive and negative influence degree of each feature is calculated by SHAP model in Python, and the result is shown in Figure 15;

[0229] Among them, the environment setting when using SHAP library for model interpretation in Python, and ensuring that the drawn graphics support Chinese display and correct display of positive and negative signs;

[0230] S68, identify key influencing factors, in the fitting results of GWR model, select the regression results with high R 2 value, such as Table 11. The results show that among the three natural geographical factors, the influence of slope on ecological benefit is the highest, followed by elevation, and the influence of slope direction is smaller;

[0231] S69, based on the results of S67, the important feature ranking of urban blue-green space pattern factor is AREA, PLAND, SPLIT, CONTIG, ENN, FRAC, CIRCLE, PLADJ, and ED. To further avoid the overfitting phenomenon of the model, the above features and corresponding ecological benefits are constructed into a new data set, and the model is trained and interpreted to obtain Figure 16. The top five important features in the blue-green space pattern patch layer are AREA, CONTIG, ENN, FRAC, and CIRCLE, among which the highest ranking is the area index of blue-green patch; the top four important features in the blue-green space pattern type layer are PLAND, SPLIT, ED, and PLADJ, among which the highest ranking is the proportion of blue-green space area.

[0232] Table 11 GWR Model Fitting Results

[0233] It should be noted that the above merely illustrates the technical idea of the present application, and cannot be used to limit the protection scope of the present application. Those skilled in the art can make several improvements and refinements without departing from the principle of the present application, and these improvements and refinements all fall within the protection scope of the present application.

Claims

1. A method for identifying key influence factors of ecological benefits of urban blue-green space, characterized in that, The method comprises the following steps: S1, constructing an urban blue-green space ecological benefit evaluation index set; S2, based on the index set constructed in S1, calculating the index data in the target region and performing standardization processing to establish an urban blue-green space ecological benefit standardized data matrix; S3, based on the ecological benefit standardized data matrix in S2, determining the subjective weight of the index by using the analytic hierarchy process, determining the objective weight by using the entropy weight method, and obtaining the weight of each index by combination weighting; S4, constructing an entropy weight-TOPSIS model to calculate the weighted urban blue-green space ecological benefit standardized data matrix, so as to complete the evaluation of the urban blue-green space ecological benefit; S5, establishing an urban blue-green space ecological benefit influence factor set, including two categories of natural factors and spatial pattern factors; S6, based on the urban blue-green space ecological benefit evaluation result in S4 and the influence factor set in S5, combining the GWR model and the XGBoost model in machine learning to identify the key influence factors.

2. The method according to claim 1, wherein, In S1, the target layer of the urban blue-green space ecological benefit evaluation index set is ecological benefit; the criterion layer includes rainwater resilience, cooling benefit, carbon sink benefit and biodiversity; the index layer includes surface runoff, flooding risk, surface temperature, net primary productivity, carbon storage and habitat quality.

3. The method of claim 1, wherein the method is characterized by, The steps for establishing the urban blue-green space ecological benefit standardized data matrix in S2 include: S21, obtaining remote sensing image data, land use data, normalized vegetation index NDVI data, vegetation type data, air temperature data, precipitation data and daily radiation data, urban soil type data and urban land cover data in the target region; S22, based on the land use data, urban soil type data, urban land cover data and precipitation data obtained in S21, the surface runoff of the target region is calculated according to the SCS-CN runoff model; S23, according to the remote sensing image data in S21, the remote sensing image under the continuous sunny and rainy scenes is calculated to represent the flooding risk in the study area; S24, based on the spectral data in the remote sensing image data in S21, the surface temperature of the target region is obtained by using the atmospheric correction method; S25, by using the land use data, normalized vegetation index NDVI data, vegetation type data, air temperature data, precipitation data and daily radiation data obtained in S21, the net primary productivity of the target region is calculated by using the CASA model; S26, the spatial distribution of carbon storage is generated by using the carbon storage module of the InVEST model; S27, the habitat quality of the target region is evaluated by using the habitat quality module in the InVEST model.

4. The method according to claim 3, wherein, In S22, the calculation steps of the surface runoff of the target region are as follows: S221, the maximum possible water storage capacity is calculated according to the CN value of the urban and residential land, and the soil of different urban land types in the target region is valued combined with the land use data; The soil cover composite number is calculated according to the following formula: In the formula, S is the maximum possible water storage capacity, which represents the maximum potential water storage capacity of the soil; S222, combined with the precipitation data in S21, the surface runoff of the target region is calculated according to the SCS-CN runoff model calculation formula; SCS-CN model calculation formula is as follows: In the formula, Q represents the surface runoff, the unit is millimeter; P represents the rainfall, the unit is also millimeter; λ indicates the initial loss rate, which is a dimensionless number, and λ=0.2 is usually used for calculation; Secondly, the step of calculating the flooding risk of the target area in S23 is: S231, according to the remote sensing image data in S21, select the remote sensing image under the continuous sunny day and after the rain respectively, and carry out land use classification; S232, select the water body classification data of the land use classification results in S231 under the continuous sunny day and after the rain, remove the constant water area in the target area, obtain the seasonal water body, and obtain the waterlogging area after subtraction; S233, after the waterlogging area is converted into points by using the "point density analysis" tool in the GIS software, the point density is calculated to obtain the waterlogging point density, which represents the flooding risk of the waterlogging area; Thirdly, in S24, the step of obtaining the surface temperature is: S241, through the spectral data containing surface temperature information in the remote sensing image data in S21, the "FLAASH Atmospheric Correction" tool of ENVI software is used to carry out atmospheric correction processing on the remote sensing image data based on spectral band; S242, the corrected remote sensing data is converted into surface temperature by using temperature inversion algorithm, and the real temperature of the surface is calculated by using Planck formula; Wherein, the expression of the thermal infrared radiation brightness value received by the satellite sensor can be written as: L λ = [ε·B(T) + (1 - ε)L↓]·τ + L↑ Wherein, the atmospheric upward radiation brightness L↑, ε is the surface emissivity, T is the surface true temperature, B(T) is the thermal radiation brightness of the blackbody at T derived from Planck's law, τ is the transmittance of the atmosphere in the thermal infrared waveband, the radiation brightness B(T) of the blackbody at T in the thermal infrared waveband is: Secondly, the Planck calculation formula is as follows: In the formula, K1 and K2 are specific calibration constants of Landsat sensor; Then, in S25, the step of calculating the net primary productivity is: S251, in the GIS software, the precipitation, temperature and solar radiation data of the target year are processed by Kriging interpolation, and the band synthesis is carried out to obtain the annual average grid picture; Then, the precipitation, temperature and solar radiation data grid picture are obtained by cutting according to the range of the target area; S252, by using CASA model, inputting the precipitation, temperature and solar radiation data, normalized vegetation index NDVI picture and vegetation data, setting the parameters in the static file, including the maximum and minimum values of normalized vegetation index NDVI data, maximum light energy utilization rate and vegetation index SR, the distribution of net primary productivity in the target area is obtained; Secondly, in S26, the step of obtaining the spatial distribution of carbon storage is: S261, based on the land use data in S21, the carbon storage module of InVEST model is used to calculate the total carbon storage of each type of land based on different land types and corresponding carbon storage data in the target area; S262, the calculated carbon storage data is classified and summarized according to the land type to generate the spatial distribution map of carbon storage; Wherein, the formula for calculating the carbon storage by InVEST model is: C t = C a + C b + C s + C d wherein: 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 belowground carbon storage (t / hm 2 ); C s is the soil carbon storage (t / hm 2 ); C d is the dead organic carbon storage (t / hm 2 ); Finally, in S27, the step of evaluating the habitat quality of the target area is: S271, the habitat quality module in InVEST model is used for analysis, and the biological diversity is evaluated according to the advantages and disadvantages of the obtained habitat quality; Wherein, the formula for calculating the habitat quality score is: Q xj = H j (1-D xj ) where Q xj is the habitat quality score for grid x in land use type j, H j is the habitat suitability score for land use type j, D xj is the threat level for grid x; S272, the threat source data is a factor threatening the ecological land in the land use type, the threat area is valued as 1 by using the GIS software, and the non-threat area is valued as 0; In the formula, the grid threat level calculation formula is: where r y is the threat intensity of the grid y, ω r is the relative weight of the threat source r, β x is the protected level of the grid x, S jr is the sensitivity of the land use type j to the threat source r, i rxy is the impact level of the grid y in the threat source r on the grid x; In the formula, the calculation formula of the influence level of the grid y in the threat source r on the grid x is as follows: In the formula, d xy is the Euclidean distance of grid x and grid y, d rmax is the maximum influence distance of threat source r.

5. The method of claim 1, wherein, In S3, the step of obtaining the weight of each index is: In S31, a plurality of sampling points are established, and the calculation results of the urban blue-green space ecological benefit standardization data matrix in S2 are randomly extracted to obtain a data set; S32, the entropy weight method is used to calculate the entropy value of each index in the data set, and then the objective weight of each index is obtained; wherein the specific steps of using the entropy weight method to calculate the entropy value of each index are as follows: first, data standardization, under the condition that n evaluation indexes and m evaluation objects are given, the original information matrix X=(x ij )m×n is established, and data standardization is carried out by using normalization, and the expression is as follows: In the formula, Y xj Xi,j is the standardized index of the i-th evaluation index and the j-th coordinate; x ij Xi,j is the original data of the i-th evaluation index and the j-th coordinate; Secondly, determining the feature weight is a key link in the evaluation process; in order to calculate the proportion of the jth index of the ith evaluation object, the weight p of the index needs to be defined first ij , and the calculation is carried out through Y ij ; Then, the information entropy E is calculated by using the weight p of each index ij j ; ​ In the formulae, m number of evaluation objects, E j is the information entropy of each evaluation index. Finally, by information entropy E j , the entropy weight method is used to determine the index weight, and the objective weight β j of the evaluation index is calculated n , β j is the weight corresponding to the jth index; In the formula, n is the number of evaluation indexes, β j is the entropy weight value of the entropy weight method; S33, in the AHP hierarchical analysis, the rain flood resilience, cooling benefit, carbon sink benefit, and biodiversity are set as equal weight; by comprehensively considering the calculation results of the AHP hierarchical analysis method and the entropy weight method, the final weight of each index is obtained; wherein the formula for obtaining the final weight of each index is: In the formula, α j is the weight of the AHP hierarchical analysis method corresponding to the jth index, and β j is the weight of the entropy weight method corresponding to the jth index.

6. The method of claim 1, wherein, In S4, the specific steps of constructing the entropy weight-TOPSIS model to calculate the weighted urban blue-green space ecological benefit standardization data matrix are: In S41, a program is written in the Matlab software, a plurality of sampling points are randomly selected, the index weight result of S3 is combined with the calculation result of the urban blue-green space ecological benefit standardization data matrix of S2, and 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; The calculation steps of the TOPSIS evaluation method are: Firstly, the weighted normalized decision matrix is constructed according to the factor weight determined in S3, and the matrix formula is: V = X ij • ω = [v ij ] m·n In the formula, V is a weighted matrix, X ij is the ith evaluation index; Second, the positive and negative ideal solutions are determined using the calculated weighted matrix V: Then, according to the positive and negative ideal solutions V + , V - , the Euclidean distance is used to calculate the distance between the positive and negative ideal solutions: where D j + D is the positive ideal solution distance j - is the negative ideal solution distance Finally, the positive and negative ideal solution distances are calculated, and the closeness is calculated: The approximation T between the calculated formula and the ideal solution is... j T j The value range of T is 0 to 1; j The closer the value is to 1, the better the ecological benefits. In S42, the index weight of the evaluation result as the positive ideal solution is selected for final calculation to obtain the final weight result; and then the urban blue-green space ecological benefit is calculated in combination with the urban blue-green space ecological benefit standardization data matrix in S2 to obtain the ecological benefit evaluation result.

7. The method of claim 1, wherein the method is characterized by, In S5, the urban blue-green space ecological benefit influence factors include "space pattern factors" and "natural factors"; the space pattern factors are divided into patch layer factors and type layer factors; wherein, the patch layer factors include area, perimeter, fractal dimension, near-circle index, adjacent index, patch distance; the type layer factors include area ratio, maximum patch index, edge density, landscape shape index, similar adjacency ratio, aggregation degree, separation degree, cohesion index, and fragmentation; the natural factors are composed of elevation, slope and slope direction.

8. The method of claim 1, wherein the method is characterized by, In S6, the specific steps of identifying the key influence factors by combining the GWR model and the XGBoost model in machine learning are: In S61, data sets are constructed based on the urban blue-green space natural factors and pattern factors in S5 and the corresponding ecological benefit evaluation results calculated in S4; in the GIS software, random sampling points are created by using the random sampling tool, and each influence factor and the ecological benefit evaluation result are extracted to the sampling points to form data set 1 by excluding points not belonging to the urban blue-green space; In S62, the global Moran index tool is used to calculate the spatial autocorrelation; The formula for calculating the Moran's index is as follows: In the formulae, n is the total number of spatial units, y i and y j respectively represent the attribute value of the i-th spatial unit and the j-th spatial unit, w for the mean of all spatial unit attributes ij for spatial weight values; In S63, a VIF calculation program is written in the Python software to calculate the multicollinearity of each influence factor in data set 1 in S61; wherein the formula for calculating VIF is: wherein R 2 is the coefficient of determination for the model; In S64, the variables with VIF<5 are used to construct data set 2 as independent variables, data set 2 is constructed with the ecological benefit evaluation results as dependent variables, the GWR model is introduced for regression analysis, the Gaussian function is used as the weight function of the model, and the AIC method is used to select a bandwidth that can effectively balance the bias and variance and avoid overfitting. GWR model calculation formula is: where Y i is the explained value of the urban ecological benefit, (u i ,v i ) is the geographical coordinate of city i; X ik is the independent variable explained value of city i, is the influence factor; β k (u i ,v i ) is the regression parameter of unit i at the region centroid (u i ,v i ); and ε i is the random error term. In the formula, the calculation formula for bandwidth determination using the AIC method is: Wherein, n represents the sample point number, RSS represents the residual sum of squares, and p is the number of unknown parameters; on this basis, the AIC calculation mode fused with the GWR model is as follows: wherein For GWR model maximum likelihood variance estimation: S65, based on S63, the factor of 5≤VIF<10 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 the XGBoost model is introduced for training; the Bayesian optimization method is used to adjust the XGBoost model hyperparameters; S66, the model accuracy is verified, the mean absolute error, the root mean square error and the determination coefficient are selected as the evaluation indexes of the prediction effect, and the ten-fold cross-validation technology is used to evaluate the generalization performance of the model; S67, the positive and negative influence degree of each feature is calculated by the SHAP model in Python software; S68, in the fitting results of the GWR model, select R 2 regression results with high values as key influencing factors; S69, in order to further avoid the overfitting phenomenon of the model, the city blue-green space pattern factor with high influence degree after the important feature sorting in S67 is combined with the corresponding ecological benefit to construct a new data set, and the model is trained and interpreted to obtain the importance degree sorting of the influence factors.

Images