A Method for Identifying Key Influencing Factors of Ecological Benefits of Urban Blue-Green Spaces

Abstract

This invention discloses a method for identifying key influencing factors of urban blue-green space ecological benefits. First, a set of evaluation indicators for urban blue-green space ecological benefits is constructed. Second, the data for each indicator are calculated and standardized to establish a standardized data matrix of urban blue-green space ecological benefits. Third, the subjective weights of the indicators are determined using the analytic hierarchy process (AHP), and the objective weights are determined using the entropy weight method. The weights of each indicator are obtained through combined weighting. Then, an entropy weight-TOPSIS model is constructed to calculate the weighted standardized data matrix of urban blue-green space ecological benefits, completing the evaluation of urban blue-green space ecological benefits. Next, a set of influencing factors for urban blue-green space ecological benefits is established. Finally, based on the evaluation results, the GWR-XGBoost model is used to identify key influencing factors. This invention combines multi-criteria decision analysis, a geographically weighted regression model, and machine learning techniques to quickly and accurately obtain the key influencing factors of urban blue-green space ecological benefits.

Description

Technical Field

[0001] This 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 Technology

[0002] Urban blue-green spaces, comprised of both blue and green spaces, form the ecological foundation of a city and are a crucial aspect of high-quality urban development. Viewing the urban blue-green system as spatial infrastructure, and using it to coordinate ecological, economic, and social development, has become a broad global consensus for cities addressing the challenges of climate change and rapid urbanization. Currently, "ecological benefits" have been explicitly incorporated into the evaluation system for socio-economic development, attracting widespread attention from society. How to effectively utilize limited urban blue-green spaces to provide higher-quality ecological benefits has become an urgent problem to be solved.

[0003] However, existing ecological benefit assessments mostly focus on large-scale spaces, and their indicators and methods are not applicable to urban and zonal scales. Research on urban blue-green spaces is still in its early stages. Previous ecological benefit assessment methods have largely focused on single land use types or single objectives, lacking comprehensiveness and coordination. These methods are unsuitable for evaluating the ecological benefits of multi-objective blue-green spaces and cannot accurately identify the multi-objective factors influencing ecological benefits. Therefore, this invention aims to provide a method for identifying key influencing factors of urban blue-green space ecological benefits, which can quickly and accurately obtain these key influencing factors based on a comprehensive evaluation of urban blue-green space ecological benefits. Summary of the Invention

[0004] To address the aforementioned issues, this invention discloses a method for identifying key influencing factors of the ecological benefits of urban blue-green spaces. By combining multi-criteria decision analysis, geographically weighted regression models, and machine learning techniques, this method can quickly and accurately obtain key influencing factors of the ecological benefits of urban blue-green spaces, which is helpful in formulating strategies for the optimization and reorganization of blue-green spaces.

[0005] To achieve the above objectives, 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 spaces includes the following steps:

[0007] S1, Construct a set of evaluation indicators for the ecological benefits of urban blue-green spaces;

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

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

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

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

[0012] S6, based on the ecological benefit evaluation results of urban blue-green space in S4 and the influencing factor set in S5, uses the GWR (Geographically Weighted Regression) and XGBoost (Extreme Gradient Boosting) models in machine learning to identify key influencing factors.

[0013] Furthermore, in S1, the target layer of the evaluation index set for the ecological benefits of urban blue-green spaces is: ecological benefits; the criteria layer includes: stormwater resilience, cooling benefits, carbon sequestration benefits, and biodiversity; and the indicator layer includes: surface runoff, inundation risk, surface temperature, net primary productivity, carbon storage, and habitat quality.

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

[0015] S21, acquire remote sensing image data, land use data, normalized difference vegetation index (NDVI) data, vegetation type data, temperature data, precipitation data, solar radiation data, urban soil type data, and urban land cover data within the target area.

[0016] S22, based on the land use data, urban soil type data, urban surface cover data, and precipitation data obtained from S21, calculates the surface runoff of the target area 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 post-rain scenarios in the remote sensing image data of S21, 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 obtained by inversion using the atmospheric correction method;

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

[0020] S26, based on land use data from 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 uses the Habitat Quality module in the InVEST model to analyze and assess the habitat quality of the target area.

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

[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 the soil of different urban land types in the target area in combination with land use data;

[0024] The formula for calculating the soil cover composite number is as follows:

[0025]

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

[0027] Based on precipitation data from S222 and S21, the surface runoff in the target area is calculated using the SCS-CN runoff generation model formula. The SCS-CN runoff generation model formula is as follows:

[0028]

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

[0030] Secondly, in S23, the steps for calculating the flooding risk of the target area are as follows:

[0031] S231. Based on the remote sensing image data of S21, land use classification is carried out by selecting remote sensing images under continuous sunny days and after rain scenarios respectively.

[0032] S232: Select water body classification data from the land use classification results under the continuous sunny day and post-rain scenarios in S231, remove the perennial water areas in the target area to obtain seasonal water bodies, and subtract them to obtain the waterlogged areas.

[0033] S233: The "point density analysis" tool in GIS software is used to convert the waterlogged area into points, and the point density is calculated to obtain the waterlogged point density, which is used as the flood risk of the waterlogged area.

[0034] Furthermore, in S24, the steps for obtaining the surface temperature are as follows:

[0035] S241: Using the spectral data containing surface temperature information in the S21 remote sensing image data, atmospheric correction processing of the remote sensing image data is performed based on the spectral band using the "FLAASH Atmospheric Correction" tool of ENVI software.

[0036] S242, the temperature inversion algorithm is used to convert the corrected remote sensing data into surface temperature, and Planck's formula is used to calculate the true surface temperature.

[0037] The expression for the thermal infrared radiation brightness value received by the satellite sensor can be written as:

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

[0039] Where L↑ is the upward atmospheric radiance, ε is the surface emissivity, T is the true surface temperature, B(T) is the thermal radiance of a blackbody at temperature T derived from Planck's law, and τ is the atmospheric transmittance in the thermal infrared band, then the radiance B(T) of a blackbody at temperature T in the thermal infrared band is:

[0040]

[0041] Planck's calculation formula is as follows:

[0042]

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

[0044] Then, in S25, the steps for calculating net primary productivity are as follows:

[0045] S251. In GIS software, monthly precipitation, temperature and solar radiation data for the target year are processed by Kriging interpolation and band synthesis is performed to obtain an annual average raster map; then, the raster map of precipitation, temperature and solar radiation data is obtained by cropping according to the target area.

[0046] S252 uses the CASA model, inputting precipitation, temperature and solar radiation data, normalized difference vegetation index (NDVI) map, and vegetation data, and setting parameters in the static file, including the maximum and minimum values ​​of the normalized difference vegetation index (NDVI) data, maximum light energy utilization rate, and vegetation index (SR), to 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 as follows:

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

[0049] S262, the calculated carbon storage data are classified and summarized according to land type to generate a spatial distribution map of carbon storage;

[0050] Among them, the formulas for calculating carbon storage using the InVEST model are divided into:

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

[0052] In the formula: C t Total carbon storage of the site (t / hm) 2 ); C a Carbon storage in the aboveground portion (t / hm) 2 ); C b underground carbon storage (t / hm) 2 ); C s Soil carbon storage (t / hm) 2 ); C d Dead organic carbon storage (t / hm) 2 );

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

[0054] S271 uses the habitat quality module in the InVEST model for analysis and assesses biodiversity based on the quality of the obtained habitat.

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

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

[0057] In the formula Q xj H represents the habitat quality score of raster x in land use type j. j D represents the habitat suitability score for land use type j. xj The threat level of grid x;

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

[0059]

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

[0061] The formula for calculating the influence level of grid y in threat source r on grid x is as follows:

[0062]

[0063] In the formula, d xy Let d be the Euclidean distance between grid x and grid y. rmax This represents the maximum influence distance of the threat source r.

[0064] Furthermore, in S3, the steps to obtain the weights of each indicator are as follows:

[0065] S31. Establish several sampling points and randomly extract the calculation results of the standardized data matrix of urban blue-green space ecological benefits from S2 to obtain the dataset.

[0066] S32 uses the entropy weight method to calculate the entropy value of each indicator, thereby obtaining the objective weight of each indicator.

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

[0068]

[0069] In the formula, Y ijThese are the standardized indices for the i-th evaluation index and the j-th coordinate, respectively; x ij This is the original data for the j-th coordinate of the i-th evaluation index.

[0070] Secondly, determining the 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 the calculation.

[0071]

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

[0073]

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

[0075] Finally, through information entropy E 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 Let be the weight corresponding to the j-th indicator.

[0076]

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

[0078] In the AHP (Analytic Hierarchy Process) analysis (S33), the four ecological objectives—rainfall resilience, cooling benefits, carbon sequestration benefits, and biodiversity—are assigned equal weights. By comprehensively considering the calculation results of the AHP and entropy weight method, the final weights of each indicator are derived. The formulas for deriving the final weights of each indicator are as follows:

[0079]

[0080] In the formula, α j Let β be the weight of the AHP (Analog-Philippines Hierarchical Analysis) method corresponding to the j-th indicator. j Let be the weight of the entropy weight method corresponding to the j-th index.

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

[0082] S41. Write a program in Matlab to randomly select several sampling points. Combine the index weight results from S3 with the calculation results of the standardized data matrix of urban blue-green space ecological benefits from S2 to calculate the ecological benefits of each point. Based on the TOPSIS evaluation method, obtain the positive and negative ideal solutions and the evaluation results of the ecological benefits corresponding to the sampling points.

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

[0084] First, in order to comprehensively evaluate the various evaluation objects, a weighted normalized decision matrix is ​​constructed based on the factor weights determined in S3. The matrix formula is as follows:

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

[0086] In the formula, V is the weighting matrix, and X... ij Let be the i-th evaluation index.

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

[0088]

[0089] Then, based on the positive and negative ideal solutions V + V - The distance between the positive and negative ideal solutions is calculated using Euclidean distance:

[0090]

[0091] In the formula, The distance of the ideal solution. It is the distance to the negative ideal solution.

[0092] Finally, the distances between the positive and negative ideal solutions are calculated, and the closeness is calculated:

[0093]

[0094] The approximation T between the calculated formula and the ideal solution is... j T j The value range of T is 0 to 1. In this invention, T... j The closer the value is to 1, the better the ecological benefits.

[0095] S42, select the index weights that have positive ideal solutions for evaluation results and perform final calculations to obtain the final weight results; then combine the standardized data matrix of urban blue-green space ecological benefits in S2 to calculate the ecological benefits of urban blue-green space and obtain the ecological benefit evaluation results.

[0096] Furthermore, the ecological benefit influencing factors of urban blue-green space in S5 includes two categories: "spatial pattern factors" and "natural factors." Spatial pattern factors are further divided into patch layer factors and type layer factors. 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, clustering degree, separation degree, cohesion index, and fragmentation degree. Natural factors consist of elevation, slope, and aspect.

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

[0098] 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.

[0099] S62 uses the global Moran index tool to calculate its spatial autocorrelation.

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

[0101]

[0102] 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.

[0103] S63. Write a VIF calculation program in Python software to calculate the multicollinearity of each influence factor in dataset 1 of S61.

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

[0105]

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

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

[0108] The GWR model calculation formula is as follows:

[0109]

[0110] In the formula, Y represents the explanatory value of the ecological benefits of city i, (u i ,v i X represents the geographic coordinates of city i; ik β represents the explanatory value of the independent variable for city i, and β represents the impact factor. k (u i ,v i ) Unit i region centroid (u i ,v i The regression parameters at ); ε i This is the random error term.

[0111] The formula for bandwidth determination using the AIC (Akaike Information Criterion) method is as follows:

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

[0113] Where n represents the number of sample points, RSS represents the sum of squared residuals, and p is 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, uses factors with 5 ≤ VIF < 10 as 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, constructing a new dataset 3, which is then imported into the XGBoost model for training. Bayesian optimization is used to adjust the hyperparameters of the XGBoost model.

[0118] S66. Model accuracy is validated by selecting mean absolute error, root mean square error, and coefficient of determination as evaluation metrics for prediction performance, and using ten-fold cross-validation to assess the model's generalization performance.

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

[0120] S68, In the fitting results of the GWR model, select R. 2 Regression results with high values ​​are considered key influencing factors.

[0121] To further avoid overfitting of the model, a new dataset was constructed by ranking the important features in S67, including the urban blue-green spatial pattern factors with high influence and their corresponding ecological benefits. After training the model, the dataset was interpreted to obtain the ranking of the importance of the influencing factors.

[0122] The beneficial effects of this invention are as follows:

[0123] This invention provides a method for identifying key influencing factors of the ecological benefits of urban blue-green spaces, applicable to ecological benefit assessment and influencing factor identification at the urban scale. The method comprehensively evaluates and obtains the weights of each ecological benefit indicator through the AHP method, entropy weight method, and TOPSIS method, which is more objective than traditional weight determination methods and avoids interference from subjective human factors. The introduction of machine learning, an advanced method, accurately and quickly identifies key factors affecting the ecological benefits of urban blue-green spaces, contributing to the formulation of strategies for blue-green space optimization and restructuring. Attached Figure Description

[0124] Figure 1 This is a flowchart of the method for constructing an evaluation system for the ecological benefits of urban blue-green spaces according to the present invention;

[0125] Figure 2 This is a specific remote sensing data area map of an embodiment of the present invention;

[0126] Figure 3 This is a surface runoff map of an embodiment of the present invention;

[0127] Figure 4 This is a flooding risk diagram of an embodiment of the present invention;

[0128] Figure 5 This is a surface temperature inversion diagram from an embodiment of the present invention;

[0129] Figure 6 This is a diagram showing the carbon solids content of an embodiment of the present invention;

[0130] Figure 7 This is a carbon storage diagram of an embodiment of the present invention;

[0131] Figure 8This is a biodiversity diagram of an embodiment of the present invention;

[0132] Figure 9 This is a diagram showing the ecological benefit evaluation results of an embodiment of the present invention;

[0133] Figure 10 This is an elevation, slope, and aspect diagram of an embodiment of the present invention;

[0134] Figure 11 This is a calculation diagram of the patch layer index in an embodiment of the present invention;

[0135] Figure 12 This is a calculation diagram of the indicator type layer of the implementation case of this invention;

[0136] Figure 13 This is a flowchart of the GWR-XGBoost model construction process in an embodiment of the present invention;

[0137] Figure 14 This is a graph showing the results of spatial autocorrelation analysis of the embodiments of this invention;

[0138] Figure 15 This is a ranking diagram of important features of spatial pattern factors in an embodiment of the present invention;

[0139] Figure 16 This is a ranking diagram of key spatial pattern factors and important features in an embodiment of the present invention. Detailed Implementation

[0140] The present invention will be further illustrated below with reference to the accompanying drawings and specific embodiments. It should be understood that the following specific embodiments are for illustrative purposes only and are not intended to limit the scope of the invention.

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

[0142] S1 constructs a set of evaluation indicators for the ecological benefits of urban blue-green spaces. The target layer is ecological benefits; the criteria layer includes stormwater resilience, cooling benefits, carbon sequestration benefits, and biodiversity; and the indicator layer includes surface runoff, inundation 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 Spaces

[0144]

[0145] S2 selects Nanjing City as the case area. Based on the indicator set constructed in S1, it calculates and standardizes the data of various indicators in the target area, establishes a standardized data matrix of urban blue-green space ecological benefits, and obtains surface runoff, stormwater 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. Acquire remote sensing imagery, land use data, Normalized Difference Vegetation Index (NDVI) data, vegetation type data, temperature data, precipitation data, solar radiation data, urban soil type data, and urban land cover data within the target area. The data sources for this case are shown in Table 2, and the selected remote sensing data area for Nanjing City is as follows: Figure 2 The specific process is as follows:

[0147] Table 2 Data Sources

[0148]

[0149]

[0150] First, land use data: Based on the Landsat 8 (2020) remote sensing image of Nanjing urban area downloaded from the geospatial cloud platform, supervised classification and manual visual interpretation correction were performed 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, temperature, precipitation, and solar radiation data: data from 13 meteorological stations in Nanjing were obtained from the Geographic Remote Sensing Ecosystem Network;

[0152] Secondly, vegetation type data: vegetation type data for Nanjing was obtained from the geographic remote sensing ecological network, with an accuracy of 30m×30m; similarly, normalized vegetation index (NDVI) data: in Google Earth Engine, based on Landsat 8 remote sensing data, after correcting and declouding the remote sensing images, the normalized vegetation index (NDVI) data for Nanjing was obtained according to the calculation formula of the normalized vegetation index (NDVI), with an accuracy of 30m×30m.

[0153] The formula for calculating the Normalized Difference Vegetation Index (NDVI) is as follows:

[0154]

[0155] In the formula, NIR is the reflectance in the near-infrared band, and RED is the reflectance in the red band.

[0156] Finally, soil data: Download soil classification data from the U.S. Department of Agriculture classification system in Google Earth Engine, with a resolution of 30m x 30m.

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

[0158] Table 3 Soil Data

[0159]

[0160]

[0161] First, based on the CN values ​​of urban and residential land, as shown in Table 3, the maximum possible water storage capacity is calculated. Then, soil values ​​are assigned to different urban land use types within the target area, combined with land use data. The formula for calculating the soil cover composite number 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 soil of different urban land use types in the study area was assigned values ​​using the "grid calculator" tool in the GIS software;

[0165] Finally, combining the precipitation data of Nanjing City (S21), and based on the calculation formula of the SCS-CN runoff model, the surface runoff in the study area was calculated using the "raster calculator" tool in GIS software. The calculation formula of the SCS-CN runoff model is as follows:

[0166]

[0167] In the formula, Q represents surface runoff in millimeters (mm); P represents rainfall in millimeters (mm); and λ represents initial loss rate, which is a dimensionless value and is usually calculated using λ = 0.2.

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

[0169] First, based on the S21 remote sensing image data, land use classification was carried out by selecting remote sensing images under continuous sunny days and after rain scenarios.

[0170] Then, water body classification data from land use classification results under continuous sunny and post-rain scenarios in S231 are selected, and the perennial water areas within the target area are removed to obtain seasonal water bodies. The difference is then used to obtain the waterlogged areas.

[0171] Finally, after converting the flooded area into points using GIS tools, the point density is calculated using the "point density analysis" tool to obtain the flood point density, thereby characterizing the flood risk of the flooded area.

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

[0173] First, using the spectral data containing surface temperature information in the S21 remote sensing image data, atmospheric correction was performed on the remote sensing image data based on spectral bands using the "FLAASH Atmospheric Correction" tool of the ENVI software platform.

[0174] Secondly, a temperature inversion algorithm was applied to convert the corrected remote sensing data into land surface temperature, and Planck's formula was used to calculate the land surface temperature. Then, the accuracy of the inverted land surface temperature data was verified using ground observation data, and the results were analyzed and applied.

[0175] The expression for the thermal infrared radiation brightness value received by the satellite sensor can be written as:

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

[0177] Where L↑ is the upward atmospheric radiance, ε is the surface emissivity, T is the true surface temperature, B(T) is the thermal radiance of a blackbody at temperature T derived from Planck's law, and τ is the atmospheric transmittance in the thermal infrared band, then the radiance B(T) of a blackbody at temperature T in the thermal infrared band is:

[0178]

[0179] Planck's calculation formula is as follows:

[0180]

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

[0182] S25, using land use data, Normalized Differential Vegetation Index (NDVI) data, vegetation type data, temperature data, precipitation data, and solar radiation data obtained from S21, calculates the net primary productivity of the target area using the CASA model. The results are as follows: Figure 6 ;

[0183] First, using GIS software, the meteorological data of S21 was processed, 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 bands of the target year were synthesized into an annual average raster, and resampled to a spatial resolution of 30m×30m. Raster maps of precipitation, temperature, and solar radiation data for the target year were then created by cropping the data according to the study area.

[0185] Finally, in the ENVI software platform, using the CASA model, we input precipitation, temperature, solar radiation data, Normalized Difference Vegetation Index (NDVI) data, and vegetation data. We set the parameters in the static file, including the maximum and minimum values ​​of the NDVI data, maximum light energy utilization, and vegetation index SR, to obtain the spatial distribution map of net primary productivity within the actual case study, as shown below. Figure 6 .

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

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

[0188] In the formula: C t Total carbon storage of the site (t / hm) 2 ); C a Carbon storage in the aboveground portion (t / hm) 2 ); C b underground carbon storage (t / hm) 2 ); C s Soil carbon storage (t / hm) 2 ); Cd Dead organic carbon storage (t / hm) 2 ).

[0189] Table 4 Carbon density of different land use types in the case study area

[0190]

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

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

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

[0194] In the formula Q xj H represents the habitat quality score of raster x in land use type j. j D represents the habitat suitability score for land use type j. xj The threat level of grid x;

[0195] The formula for calculating the threat level of a grid is as follows:

[0196]

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

[0198] The formula for calculating the influence level of grid y in threat source r on grid x is as follows:

[0199]

[0200] In the formula, d xy Let d be the Euclidean distance between grid x and grid y. rmax This represents the maximum influence distance of the threat source r.

[0201] Table 5 Threat Factors in Habitat Quality Assessment

[0202]

[0203] 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.

[0204] 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;

[0205] 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.

[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 For m×n, the data is standardized using normalization, and its expression is:

[0207]

[0208] 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;

[0209] 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;

[0210]

[0211] 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;

[0212]

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

[0214] S324, using the information entropy Ej from S323, determines the index weights using the entropy weight method, and calculates the objective weight β of the evaluation index. j =β1, β2, ...,βn , where β j The weight corresponding to the j-th indicator;

[0215]

[0216] In the formula, n is the number of evaluation indicators, and β j The entropy weight is the value of the entropy weight method.

[0217] In the AHP (Analytic Hierarchy Process) analysis, the four ecological objectives—rainfall resilience, cooling benefits, carbon sequestration benefits, and biodiversity—were assigned equal weights. By comprehensively considering the calculation results of the AHP and entropy weight methods, the final weights of each indicator were derived. The formulas for the final weights of each indicator are shown in Table 6.

[0218] The formula for the final weight is:

[0219]

[0220] In the formula, α j Let β be the weight of the AHP (Analog-Philippines Hierarchical Analysis) method corresponding to the j-th indicator. j Let be the weight of the entropy weight method corresponding to the j-th index.

[0221] Table 6 Calculation Results of Evaluation Index Weights

[0222]

[0223] S4, based on the quantitative assessment results of S2 and the weight calculation results of S3, constructs an entropy-weighted TOPSIS model to calculate a weighted standardized data matrix of urban blue-green space ecological benefits, thereby completing the evaluation of urban blue-green space ecological benefits. Figure 9 .

[0224] S41. Write a program in Matlab to randomly select several sampling points. Combine the index weight results from S3 with the calculation results of the standardized data matrix of urban blue-green space ecological benefits from S2, as shown in Table 7, to calculate the ecological benefits of each point. Based on the TOPSIS evaluation method, obtain the positive and negative ideal solutions, and the evaluation results of the ecological benefits corresponding to the sampling points.

[0225] Furthermore, in order to comprehensively evaluate the various evaluation objects, a weighted normalized decision matrix is ​​constructed based on the factor weights determined in S3. The matrix formula is as follows:

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

[0227] In the formula, V is the weighting matrix, and X... ij Let i be the i-th evaluation index;

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

[0229]

[0230] Then, based on the obtained positive and negative ideal solutions V + V - The distance between the positive and negative ideal solutions is calculated using Euclidean distance:

[0231]

[0232] In the formula, The distance of the ideal solution. It is the distance to the negative ideal solution;

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

[0234]

[0235] The approximation T between the calculated formula and the ideal solution is... j T j The value of is in the range of 0 to 1. T j The closer the value is to 1, the better the ecological benefits.

[0236] Table 7 Calculation Results of Positive and Negative Ideal Solutions

[0237]

[0238] S42, the weights of the indicators whose evaluation results are positive ideal solutions are selected for final calculation to obtain the final weight results; the values ​​of the above calculation results are correlated with the values ​​of each point in the geographic space, Kriging interpolation is used, and then combined with the standardized data matrix of urban blue-green space ecological benefits in S2 to calculate the ecological benefits of urban blue-green space, so as to obtain the ecological benefit evaluation results, such as... Figure 9 .

[0239] S5 establishes a set of influencing factors for the ecological benefits of urban blue-green spaces, comprising two main categories: natural factors and spatial pattern factors. The influencing factors for the ecological benefits of urban blue-green spaces include both "spatial pattern factors" and "natural factors." Spatial pattern factors are further divided into patch layer factors and type layer factors. 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, clustering degree, separation degree, cohesion index, and fragmentation degree. Natural factors consist of elevation, slope, and aspect, as shown in Table 8.

[0240] Table 8. Set of Influencing Factors of Ecological Benefits of Urban Blue-Green Spaces

[0241]

[0242]

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

[0244] Furthermore, using Fragstats software, a quantitative analysis of the urban blue-green spatial pattern characteristics for selected years was conducted, and the patch layer calculation results are as follows: Figure 11 The type layer calculation results are as follows Figure 12 .

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

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

[0247] 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.

[0248] S62, in GIS software, the global Moran's index tool is used to process the ecological benefit evaluation results and calculate their spatial autocorrelation. The results are as follows: Figure 14 The evaluation criteria are shown in Table 9. The spatial autocorrelation Moran's index of the ecological benefit evaluation results within the target area of ​​the case study 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 target area of ​​the case study is not random and there is a significant positive correlation with the spatial distribution. Therefore, the standard regression model is not applicable in this case, and the GWR model is more suitable.

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

[0250]

[0251] In the formula, n is the total number of spatial units, y i and y jLet 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;

[0252] Table 9. Spatial Autocorrelation Test Criteria

[0253]

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

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

[0256]

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

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

[0259] The GWR model calculation formula is as follows:

[0260]

[0261] In the formula, Y i Let i be the explanatory value for the ecological benefits of city i, (u i ,v i X represents the geographic coordinates of city i; ik β represents the explanatory value of the independent variable for city i, and β represents the impact factor. k (u i ,v i ) is the centroid of region i (u i ,v i The regression parameters at ); ε i This is the random error term;

[0262] The formula for bandwidth determination using the AIC method is as follows:

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

[0264] Where n represents the number of sample points, RSS represents the sum of squared residuals, and p is 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, uses factors with 5 ≤ VIF < 10 as 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, constructing a new dataset 3, which is then imported into the XGBoost model for training. Bayesian optimization is used to adjust the hyperparameters of the XGBoost model.

[0269] First, the dataset was standardized. To avoid overfitting, the dataset was then regularized. 80% of the data was used as the training set and 20% as the test set for model validation.

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

[0271] Specifically, 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. The objective function defines the model's evaluation metrics, such as accuracy or loss value.

[0272] 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 Hyperopt will try. The optimal hyperparameter combination is obtained by combining the previously defined parameter space `space_xgb` with the result returned by the `fmin` function.

[0273] Table 10 Model Accuracy Verification

[0274]

[0275] S66. Model accuracy is verified by selecting mean absolute error, root mean square error and coefficient of determination as evaluation indicators of prediction performance, and using ten-fold cross-validation to evaluate the generalization performance of the model.

[0276] S67, using the SHAP model in Python to calculate the positive and negative influence of each feature, the results are as follows: Figure 15 ;

[0277] This includes setting up the environment for model interpretation using the SHAP library in Python, and ensuring that the drawn graphs support Chinese character display and the correct display of plus and minus signs;

[0278] S68. Identify key influencing factors and select R from the GWR model fitting results. 2 The regression results with high values ​​are shown in Table 11. The results indicate that among the three natural geographical factors, slope has the highest impact on ecological benefits, followed by elevation, while aspect has a relatively small impact.

[0279] Based on the results of S67, S69 shows that the urban blue-green spatial pattern factors with high influence after ranking important features are AREA, PLAND, SPLIT, CONTIG, ENN, FRAC, CIRCLE, PLADJ, and ED. To further avoid overfitting, a new dataset was constructed by combining the above features with their corresponding ecological benefits. After training the model, the dataset was interpreted to obtain... Figure 16 Among the patch layer features of blue-green spatial patterns, the top five most important are AREA, CONTIG, ENN, FRAC, and CIRCLE, with the area index of blue-green patches being the most important. Among the type layer features of blue-green spatial patterns, the top four most important are PLAD, SPLIT, ED, and PLADJ, with the proportion of blue-green spatial area being the most important.

[0280] Table 11 GWR Model Fitting Results

[0281]

[0282] It should be noted that the above content merely illustrates the technical concept of the present invention and should not be construed as limiting the scope of protection of the present invention. For those skilled in the art, various improvements and modifications can be made without departing from the principle of the present invention, and all such improvements and modifications fall within the scope of protection of the claims of the present invention.

Claims

1. A method for identifying key influencing factors of the ecological benefits of urban blue-green spaces, characterized in that, Includes the following steps: S1, Construct a set of evaluation indicators for the ecological benefits of urban blue-green spaces; S2, based on the indicator set constructed by S1, calculates the data of various indicators in the target area and performs standardization processing to establish a standardized data matrix of ecological benefits of urban blue-green space; S3, based on the ecological benefit standardized data matrix of S2, uses the analytic hierarchy process to determine the subjective weights of the indicators, uses the entropy weight method to determine the objective weights, and obtains the weights of each indicator through combined weighting. S4. Construct the entropy weight-TOPSIS model to calculate the weighted standardized data matrix of urban blue-green space ecological benefits, thereby completing the evaluation of urban blue-green space ecological benefits; S5, establish a set of factors influencing the ecological benefits of urban blue-green spaces, including two major 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 set of influencing factors in S5, identifies key influencing factors by combining the GWR model and XGBoost model in machine learning. The specific steps are as follows: S61 is a dataset constructed 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 GIS software, a random sampling tool is used to create random sampling points, and each influencing factor and ecological benefit evaluation result is extracted to the sampling points. Points that do not belong to urban blue-green space are removed to form dataset 1. S62 uses the global Moran index tool to calculate its spatial autocorrelation; The Moran index is calculated using the following formula: 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; S63, Write a VIF calculation program in Python software to calculate the multicollinearity of each influence factor in dataset 1 of S61; The formula for calculating VIF is as follows: In the formula, R 2 These are the coefficients of determination for the model; S64. Dataset 2 was constructed using variables with VIF < 5 as independent variables and its ecological benefit evaluation results as dependent variables. The GWR model was introduced for regression analysis, and the Gaussian function was used as the weight function of the model. The AIC method was used to select a bandwidth that can effectively balance bias and variance and avoid overfitting. The GWR model calculation formula is as follows: In the formula, Y i Let i be the explanatory value for the ecological benefits of city i, (u i ,v i X represents the geographic coordinates of city i; ik β represents the explanatory value of the independent variable for city i, and β represents the impact factor. k (u i ,v i ) Unit i region centroid (u i ,v i The regression parameters at ); ε i This is the random error term; The formula for bandwidth determination using the AIC method is as follows: AIC = log e (RSS)+2p Where n represents the number of sample points, RSS represents the sum of squared residuals, and p is the number of unknown parameters; based on this, 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, uses factors with 5 ≤ VIF < 10 as 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 to construct a new dataset 3, which is then imported into the XGBoost model for training; the hyperparameters of the XGBoost model are adjusted using the Bayesian optimization method. S66. Model accuracy is verified by selecting mean absolute error, root mean square error and coefficient of determination as evaluation indicators of prediction performance, and using ten-fold cross-validation to evaluate the generalization performance of the model. S67, calculate the positive and negative impact of each feature using the SHAP model in Python software; S68, In the fitting results of the GWR model, select R. 2 Regression results with high values ​​are considered key influencing factors. To further avoid overfitting of the model, a new dataset was constructed by ranking the important features in S67 and constructing the urban blue-green spatial pattern factors with high influence and their corresponding ecological benefits. After training the model, the dataset was interpreted to obtain the ranking of the importance of the influencing factors.

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

3. The method for identifying key influencing factors of urban blue-green space ecological benefits according to claim 1, characterized in that, The steps for establishing a standardized data matrix of the ecological benefits of urban blue-green spaces in S2 include: S21, acquire remote sensing image data, land use data, normalized difference vegetation index (NDVI) data, vegetation type data, temperature data, precipitation data, solar radiation data, urban soil type data, and urban land cover data within the target area. S22, based on the land use data, urban soil type data, urban surface cover data, and precipitation data obtained from S21, calculates the surface runoff of the target area according to the SCS-CN runoff generation model; S23. Based on the remote sensing images of continuous sunny days and post-rain scenarios in the remote sensing image data of S21, 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 obtained by inversion using the atmospheric correction method; S25 uses the land use data, normalized difference vegetation index (NDVI) data, vegetation type data, temperature data, precipitation data, and solar radiation data obtained from S21 to calculate the net primary productivity of the target area using the CASA model. S26, a carbon storage module using the InVEST model, generates the spatial distribution of carbon storage. S27 uses the habitat quality module in the InVEST model for analysis to assess the habitat quality of the target area.

4. The method for identifying key influencing factors of urban blue-green space ecological benefits according to claim 3, characterized in that, In S22, the calculation steps for surface runoff in the target area are as follows: S221, calculate the maximum possible water storage capacity based on the CN value of urban and residential land, and assign values ​​to the soil of different urban land types in the target area in combination with land use data; The formula for calculating the soil cover composite number is as follows: In the formula, S is the maximum possible water storage, which represents the maximum potential water storage capacity of the soil; Based on precipitation data from S222 and S21, the surface runoff in the target area is calculated using the SCS-CN runoff generation model formula. The SCS-CN runoff generation model formula is as follows: In the formula, Q represents surface runoff in millimeters; P represents rainfall in millimeters; and λ represents initial loss rate, which is a dimensionless value and is usually calculated using λ = 0.

2. Secondly, the steps for calculating the flooding risk of the target area in S23 are as follows: S231. Based on the remote sensing image data of S21, land use classification is carried out by selecting remote sensing images under continuous sunny days and after rain scenarios respectively. S232: Select water body classification data from the land use classification results under the continuous sunny day and post-rain scenarios in S231, remove the perennial water areas in the target area to obtain seasonal water bodies, and subtract them to obtain the waterlogged areas. S233: The "point density analysis" tool in GIS software is used to convert the waterlogged area into points, and the point density is calculated to obtain the waterlogged point density, which is used to characterize the flood risk of the waterlogged area. Furthermore, in S24, the steps for obtaining the surface temperature are as follows: S241: Using the spectral data containing surface temperature information in the S21 remote sensing image data, atmospheric correction processing of the remote sensing image data is performed based on the spectral band using the "FLAASH Atmospheric Correction" tool of ENVI software. S242, the temperature inversion algorithm is used to convert the corrected remote sensing data into surface temperature, and Planck's formula is used to calculate the true surface temperature. The expression for the thermal infrared radiation brightness value received by the satellite sensor can be written as: L λ =[ε·B(T)+(1-ε)L↓]·τ+L↑ Where L↑ is the upward atmospheric radiance, ε is the surface emissivity, T is the true surface temperature, B(T) is the thermal radiance of a blackbody at temperature T derived from Planck's law, and τ is the atmospheric transmittance in the thermal infrared band, then the radiance B(T) of a blackbody at temperature T in the thermal infrared band is: Secondly, Planck's formula is as follows: In the formula, K1 and K2 are calibration constants specific to the Landsat sensor; Then, in S25, the steps for calculating net primary productivity are as follows: S251. In GIS software, monthly precipitation, temperature and solar radiation data for the target year are processed by Kriging interpolation and band synthesis is performed to obtain an annual average raster map; then, the raster map of precipitation, temperature and solar radiation data is obtained by cropping according to the target area. S252 uses the CASA model, inputting precipitation, temperature and solar radiation data, normalized difference vegetation index (NDVI) map, and vegetation data, and setting parameters in the static file, including the maximum and minimum values ​​of the normalized difference vegetation index (NDVI) data, maximum light energy utilization rate, and vegetation index (SR), to 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 as follows: S261, based on S21 land use data, uses the InVEST model carbon storage module to calculate the total carbon storage of each land type based on different land types and corresponding carbon pool data in the target area. S262, the calculated carbon storage data are classified and summarized according to land type to generate a spatial distribution map of carbon storage; Among them, the formulas for calculating carbon storage using the InVEST model are divided into: C t =C a +C b +C s +C d In the formula: C t Total carbon storage of the site (t / hm) 2 ); C a Carbon storage in the aboveground portion (t / hm) 2 ); C b underground carbon storage (t / hm) 2 ); C s Soil carbon storage (t / hm) 2 ); C d Dead organic carbon storage (t / hm) 2 ); Finally, in S27, the steps for assessing the habitat quality of the target area are as follows: S271 uses the habitat quality module in the InVEST model for analysis and assesses biodiversity based on the quality of the obtained habitat. The formula for calculating habitat quality score is as follows: Q xj =H j (1-D xj ) In the formula Q xj H represents the habitat quality score of raster x in land use type j. j D represents the habitat suitability score for land use type j. xj The threat level of grid x; S272, Threat source data refers to factors that pose a threat to ecological land in land use types. Using GIS software, the threatened areas are assigned a value of 1, and the non-threatened areas are assigned a value of 0. The formula for calculating the threat level of a grid is as follows: In the formula, r y For the threat intensity of grid y, ω r β represents the relative weight of threat source r. x S represents the protection level of raster x. jr For the sensitivity of land use type j to threat source r, i rxy The level of influence of grid y in threat source r on grid x; The formula for calculating the influence level of grid y in threat source r on grid x is as follows: In the formula, d xy Let d be the Euclidean distance between grid x and grid y. r max This represents the maximum influence distance of the threat source r.

5. The method for identifying key influencing factors of urban blue-green space ecological benefits according to claim 1, characterized in that, In S3, the steps to obtain the weights of each indicator are as follows: S31. Establish several sampling points and randomly extract the calculation results of the standardized data matrix of urban blue-green space ecological benefits in S2 to obtain the dataset; S32, the entropy weight method is used to calculate the entropy value of each indicator in the dataset, thereby obtaining the objective weight of each indicator; the specific steps of calculating the entropy value of each indicator using the entropy weight method are as follows: First, data standardization is performed. Given n evaluation indicators and m evaluation objects, an original information matrix X = (x ij For m×n, the data is standardized using normalization, and its expression is: In the formula, Y xj 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; Secondly, determining the 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 Through Y ij Perform calculations; Then, using the weights p of each indicator ij Calculate the information entropy, and calculate the information entropy E of each evaluation index. j ; In the formula, m is the number of evaluation objects, E j The information entropy of each evaluation indicator; Finally, through information entropy E 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; In the formula, n is the number of evaluation indicators, and β j The entropy weight is the value of the entropy weight method. S33, in the AHP (Analytic Hierarchy Process) analysis, the four ecological objectives of stormwater resilience, cooling benefits, carbon sequestration benefits, and biodiversity are assigned equal weights. By comprehensively considering the calculation results of the AHP analysis and the entropy weight method, the final weights of each indicator are derived. The formulas for deriving the final weights of each indicator are as follows: In the formula, α j Let β be the weight of the AHP (Analog-Philippines Hierarchical Organization) method corresponding to the j-th indicator. j Let be the weight of the entropy weight method corresponding to the j-th index.

6. The method for identifying key influencing factors of urban blue-green space ecological benefits according to claim 1, characterized in that, In S4, the specific steps for constructing the entropy-weighted TOPSIS model to calculate the weighted standardized data matrix of urban blue-green space ecological benefits are as follows: S41. Write a program in Matlab software to randomly select several sampling points. Combine the index weight results of S3 with the calculation results of the standardized data matrix of urban blue-green space ecological benefits of S2. According to the TOPSIS evaluation method, obtain the positive and negative ideal solutions and the evaluation results of the ecological benefits corresponding to the sampling points. The calculation steps for the TOPSIS evaluation method are as follows: First, construct the weighted normalized decision matrix based on the factor weights determined in S3. The matrix formula is as follows: V=X ij ·ω=[v ij ] m·n In the formula, V is the weighting matrix, and X... ij Let i be the i-th evaluation index; Secondly, using the calculated weighted matrix V, the positive and negative ideal solutions are determined: Then, based on the positive and negative ideal solutions V + V - The distance between the positive and negative ideal solutions is calculated using Euclidean distance: In the formula, The distance of the ideal solution. It is the distance to the negative ideal solution; Finally, the distances between the positive and negative ideal solutions are calculated, and the closeness is calculated: The approximation T between the calculated formula and the ideal solution 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. S42, select the index weights that have positive ideal solutions for evaluation results and perform final calculations to obtain the final weight results; then combine the standardized data matrix of urban blue-green space ecological benefits in S2 to calculate the ecological benefits of urban blue-green space and obtain the ecological benefit evaluation results.

7. The method for identifying key influencing factors of urban blue-green space ecological benefits according to claim 1, characterized in that, The ecological benefit influencing factors of urban blue-green space in S5 include two categories: "spatial pattern factors" and "natural factors." 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 degree, separation degree, cohesion index, and fragmentation degree; natural factors consist of elevation, slope, and aspect.

Images