A method for inverting urban greenbelt cooling service function by fusing three-dimensional greenbelt structure and dynamic population exposure

Abstract

The application discloses a kind of fusion three-dimensional green space structure and dynamic population exposure city green space cooling service function inversion method, including obtaining research area multi-source input data, remote sensing image pre-processing and heat month identification, construct spatial partition, construct three-dimensional green space structure comprehensive representation value G3D, partition inversion three-dimensional correction cooling efficiency CE3D, calculate three-dimensional correction cooling capacity CC3D, construct dynamic population exposure weight, calculate dynamic cooling benefit DCB, generate result layer and statistical table and result check and storage archive;The application no longer only with NDVI or green area representation green space cooling effect, but tree height, crown cover, building shadow and sky visible factor are included in G3D, so that cooling efficiency can reflect the difference of vegetation stereo structure, shading process and radiation exchange.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing inversion, urban thermal environment monitoring, three-dimensional urban morphology analysis, geographic information system spatial analysis, dynamic population exposure assessment and urban ecosystem service function assessment, specifically a method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure. Background Technology

[0002] Urbanization alters land cover structure, building morphology, and surface energy exchange processes, leading to an enhanced urban heat island effect. Urban green spaces, through processes such as shading, evapotranspiration, airflow regulation, and surface energy balance regulation, can reduce surrounding surface temperatures and are crucial for urban thermal environment management and natural solutions. Existing urban green space cooling studies typically utilize remote sensing imagery to retrieve surface temperature and combine it with indicators such as vegetation index, green space area, patch shape, and buffer zone temperature difference to evaluate the cooling effect. Current methods often use two-dimensional vegetation cover or surface temperature difference as the primary evaluation criteria, which can reflect the physical cooling effect of urban green spaces to some extent. However, the cooling process of urban green spaces is not only related to vegetation cover but is also influenced by three-dimensional structural factors such as tree height, canopy cover, sky visibility, building shadows, and street valley morphology. Furthermore, the actual access residents receive to cooling services varies depending on daytime, nighttime, weekdays, holidays, and different population activity patterns; static population grids cannot accurately reflect these dynamic exposure differences.

[0003] Existing technologies closely related to this application mainly fall into three categories: First, green space cooling evaluation methods based on remote sensing surface temperature and vegetation index. These methods typically acquire thermal infrared remote sensing images of the target city, retrieve the LST (Local Temperature Stress), extract NDVI (Natural Density Variable Index) or green space coverage, and then statistically analyze the average temperature difference according to green space patches, buffer zones, or administrative units, using the temperature reduction magnitude to characterize the green space cooling effect. This approach can reflect the local temperature differences between green spaces and the surrounding built environment, but it primarily remains at the level of two-dimensional surface temperature changes. Second, urban thermal environment analysis methods based on land cover and building morphology zoning. These methods can divide the city into different surface types or local climate zones and statistically analyze variables such as surface temperature, vegetation coverage, and building height within each zone. However, they typically do not further incorporate three-dimensional interaction characteristics between green spaces and the built environment, such as tree height, canopy coverage, sky visibility factor, and building shadows, nor do they form a three-dimensional structure-corrected cooling efficiency inversion process. Third, green space cooling benefit evaluation methods based on population grids. These methods can incorporate population density into assessments of cooling capacity or green space accessibility to characterize the spatial relationship between cooling services and population distribution. However, these methods typically use static resident population or annual average population data, making it difficult to distinguish between daytime working population, nighttime resident population, weekday commuter population, holiday activity population, and the differences in heat exposure among sensitive groups such as the elderly and children at different times.

[0004] Existing technologies suffer from at least the following drawbacks: First, existing methods often use NDVI, green space area, or local temperature difference as cooling evaluation indicators, neglecting the impact of three-dimensional structures such as tree height, canopy coverage, sky visibility factor, and building shadows on shading, radiation exchange, and evapotranspiration processes. This makes the estimation of green space cooling efficiency easily limited by insufficient two-dimensional vegetation indicators. Second, existing methods typically conduct statistics within administrative units, green space patches, or simple buffer zones, failing to fully consider the differences in the relationship between LST and the three-dimensional structure of green spaces under different surface types, building heights, and canopy structures. This easily leads to mixed estimation errors of green space cooling capacity under different spatial backgrounds. Third, existing technologies usually use static population rasters when calculating cooling benefits, failing to reflect the dynamic changes in population exposure under different time periods, date types, and population types. This results in insufficient spatial coupling between cooling service supply and actual population demand. Fourth, in existing technologies, LST inversion, three-dimensional structural feature construction, cooling efficiency inversion, cooling capacity estimation, and dynamic population benefit calculation are often implemented in a decentralized manner, lacking a reproducible and transferable automated technology chain, which affects the comparability between different cities, different years, and different time periods.

[0005] Therefore, the technical problem to be solved by this application is to provide a method for inverting the urban green space cooling service function by integrating multi-source remote sensing data, urban three-dimensional green space structure data and dynamic population exposure data, inverting the three-dimensional corrected cooling efficiency according to urban spatial zoning, and further constructing a method for inverting the urban green space cooling service function by integrating the three-dimensional corrected cooling capacity and the dynamic population exposure weighted cooling benefit, so as to improve the structural integrity, time-time adaptability and spatial expression accuracy of urban green space cooling function identification. Summary of the Invention

[0006] The purpose of this invention is to provide a method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure, so as to solve the problems existing in the prior art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: a method for inverting the cooling service function of urban green spaces by integrating three-dimensional green space structure and dynamic population exposure, comprising the following steps:

[0008] S1 acquires multi-source input data within the target city area, including at least remote sensing image data, land cover data, building height data, digital elevation model data, distance from water body data, tree height data, canopy coverage data, sky visibility factor data, building shadow data, dynamic population raster data, and study area boundary data.

[0009] S2 preprocesses remote sensing image data, calculates the Normalized Difference Vegetation Index (NDVI) and Land Surface Temperature (LST), and identifies the hot month or heat wave period of the target year.

[0010] S3 divides the target city into multiple spatial zones based on land cover data, building height data, and tree canopy cover structure;

[0011] S4 unifies the spatial resolution and coordinate system of NDVI, tree height, canopy coverage, building shadow ratio and sky visibility factor, and after normalization, weighted fusion is used to construct the three-dimensional green space structure comprehensive characterization value G3D.

[0012] S5 uses each spatial zone as an independent calculation unit to construct a multiple regression model with LST as the dependent variable, G3D as the core independent variable, and introduces terrain and building environment covariates to obtain the three-dimensional corrected cooling efficiency CE3D of each zone.

[0013] S6 determines the minimum value of G3D within the target city area, G3Dmin, and calculates the three-dimensional corrected cooling capacity CC3D at the zonal and city scales by combining the intensity difference between CE3D of each zone and G3D of the grid.

[0014] S7 constructs a grid-scale dynamic population exposure weight WPop for different population exposure scenarios;

[0015] S8 combines three-dimensional corrected cooling capacity with dynamic population exposure weights to calculate dynamic cooling benefit DCB at the zonal and city scales.

[0016] S9 generates spatial layers and statistical results for urban green space cooling service functions;

[0017] S10 verifies, stores, and updates the output results.

[0018] Furthermore, in step S1, multi-source input data is acquired within the target city area. All data must cover the entire study area, and the timeframes should ideally match the hottest periods. The multi-source input data includes at least:

[0019] Remote sensing image data: including optical and thermal infrared bands, used to retrieve land surface temperature and calculate vegetation indices. Optical bands are used to extract vegetation indices and land cover information, while thermal infrared bands are used for land surface temperature retrieval.

[0020] Land cover data: used to identify green spaces, water bodies, built-up land and other land surface types, including clear land use codes and boundary information.

[0021] Building height data: used to characterize the three-dimensional form of urban buildings, supporting spatial zoning, sky visibility factor calculation, and building shadow simulation.

[0022] Digital Elevation Model (DEM) data is used to characterize terrain undulations, serves as a covariate for cooling efficiency inversion, and also participates in the calculation of sky visibility factor and building shadows.

[0023] Distance data to water bodies: used to characterize the moderating effect of water bodies on surface temperature, serving as a covariate for cooling efficiency inversion; based on the water body extent in the land cover data, the straight-line distance from each grid cell to the nearest water body is calculated through Euclidean distance analysis.

[0024] Tree height data: used to characterize the vertical structural features of arbor vegetation, and is a core component of three-dimensional green space structure.

[0025] Canopy coverage data: used to characterize the horizontal coverage of tree canopies, reflecting the shading range of vegetation canopies.

[0026] Sky visibility factor data (SVF data) is used to characterize the proportion of the sky hemisphere visible from a location on the ground, reflecting the degree of sky obstruction by buildings and tree canopies.

[0027] Building shadow data: used to characterize the spatial distribution range of shadows cast by buildings and tree canopies, reflecting the actual area affected by shading.

[0028] Dynamic population raster data: used to calculate the cooling benefits under different time periods, date types, or population types.

[0029] Study area boundary data: namely, the administrative boundary or study area boundary vector data of the target city, used for cropping and unifying the scope of all spatial data.

[0030] Furthermore, in step S2, the remote sensing image data undergoes systematic preprocessing to ensure the accuracy of subsequent index calculations and temperature inversion. This specifically includes:

[0031] Cloud and cloud shadow masking: Using the quality assessment bands built into remote sensing images, identify contaminated pixels such as clouds, cloud shadows, and snow, generate mask files, remove contaminated pixels, and ensure that all pixels involved in the calculation are in a clear, cloudless state.

[0032] Radiometric calibration: Using the radiometric calibration coefficients in the image header file, the digital quantization value of the image is converted into an apparent radiance value, eliminating the radiometric error caused by the sensor's own gain and offset.

[0033] Scale factor conversion: For surface temperature or reflectivity products with scale factors, multiply by the corresponding scale coefficient according to the product instructions to convert to the actual physical quantity value.

[0034] Atmospheric correction: Atmospheric correction is performed on optical band images using a radiative transfer model to eliminate the effects of atmospheric molecular scattering, aerosol scattering and absorption on surface reflectivity and obtain accurate surface reflectivity data.

[0035] Projection transformation: Convert all input data into the same projected coordinate system to ensure that all spatial data have a consistent spatial reference and avoid spatial misalignment.

[0036] Spatial cropping: Using the boundary vector data of the study area, all remote sensing images and spatial data are cropped to retain only the data within the target city area.

[0037] Resolution unification: A resampling method is used to unify all spatial data to the same spatial resolution, ensuring that the pixels in subsequent raster calculations correspond one-to-one.

[0038] After preprocessing, the Normalized Difference Vegetation Index (NDVI) is calculated based on the corrected surface reflectance, and the Land Temperature (LST) is obtained based on thermal infrared bands or temperature products. LST and NDVI images are synthesized monthly or daily according to the target year. The average LST for each month or each high-temperature day within the target city area is statistically analyzed. The month with the highest average LST is selected as the target hot month, or consecutive high-temperature days are selected to constitute a heat wave period. Water body pixels and pixels with NDVI less than 0 are masked to reduce the interference of low-temperature effects from water bodies and anomalous pixels on the results.

[0039] The preprocessing includes cloud masking, cloud shadow masking, radiometric calibration, scale factor conversion, projection conversion, spatial cropping, and resolution unification; masking is performed on water body pixels and pixels with NDVI less than 0; the identification method for hot months or heat waves is as follows: synthesize LST images month by month or day by day, statistically analyze the average LST for each time period within the target city area, select the month with the highest average LST as the hot month, or select consecutive high-temperature days to constitute a heat wave period.

[0040] Furthermore, in step S3, based on land cover data, building height data, and tree canopy cover structure, the target city is divided into multiple spatial zones, which serve as the basic spatial units for subsequent zone regression and indicator calculation. The specific division rules are as follows:

[0041] Initial land use classification: Based on land cover data, two major land cover categories are extracted: building land and vegetation cover.

[0042] Building Classification: Within the building site, the area is divided according to building height data. A building height threshold is set, and areas with an average building height below the threshold are classified as low-rise building areas, while areas with an average building height above or equal to the threshold are classified as mid-rise building areas.

[0043] Vegetation Classification: Within the vegetation-covered area, tree vegetation and non-tree vegetation are first distinguished by combining tree height data, canopy cover data, and land cover sub-type. Vegetation areas with tree height below a threshold and no obvious tree canopy are classified as non-tree vegetation areas. For tree vegetation areas, further subdivision is based on canopy cover. A canopy cover threshold is set, and areas with canopy cover below the threshold are classified as open canopy areas, while areas with canopy cover above or equal to the threshold are classified as closed canopy areas.

[0044] Ultimately, five spatial zones are formed: non-tree vegetation zone, low-rise building zone, mid-to-high-rise building zone, open canopy zone, and closed canopy zone. Each zone has a relatively consistent underlying surface background and three-dimensional structural features.

[0045] The spatial zoning includes five categories: non-tree vegetation zone, low-rise building zone, mid-to-high-rise building zone, open canopy zone, and closed canopy zone. The zoning rules are as follows: vegetation-covered zones and building zones are distinguished based on land cover data; within the building zone, low-rise building zones and mid-to-high-rise building zones are distinguished based on building height thresholds; within the vegetation-covered zone, non-tree vegetation and tree vegetation are distinguished; within the tree vegetation zone, open canopy zones and closed canopy zones are distinguished based on canopy coverage thresholds.

[0046] Furthermore, in step S4, NDVI, tree height, canopy coverage, building shadow ratio, and sky visibility factor are unified to the same spatial resolution and coordinate system. After normalization of each variable, they are weighted and fused to obtain G3D, which is used to characterize the three-dimensional structure and shading capacity of green space.

[0047] The calculation equation for the comprehensive characterization value G3D of the three-dimensional green space structure is as follows:

[0048]

[0049] In the formula, This is the comprehensive three-dimensional green space structure representation value of the i-th grid cell; The normalized vegetation index; The normalized tree height; This represents the normalized canopy coverage. The normalized building shadow ratio; The normalized sky visibility factor; , , , , These are the weighting coefficients.

[0050] Furthermore, the normalization process employs a minimum-maximum normalization method to map each variable to the [0,1] interval, eliminating dimensional differences. Since a higher sky visibility factor value indicates weaker shading and is negatively correlated with the cooling effect of shading, the normalized value is subtracted from 1 after normalization, making its value positively correlated with the cooling contribution of the three-dimensional structure. The weighting coefficients are determined using any one of the following methods: equal weighting, entropy weighting, expert weighting, regression coefficient normalization, or machine learning feature contribution. A higher G3D value indicates a more complete three-dimensional structure of the green space in the grid, and stronger shading, evapotranspiration, and radiation regulation capabilities.

[0051] Furthermore, in step S5, a multiple regression model is constructed using each spatial partition as an independent calculation unit. The model uses LST as the dependent variable, G3D as the core explanatory variable, and introduces elevation, building height, and distance from the water body as covariates to control for interference from environmental background factors. The calculation equation of the multiple regression model is as follows:

[0052]

[0053] In the formula, For the intercept term; , , , These are the regression coefficients; The elevation of the i-th grid cell; Let i be the building height of the i-th grid cell; Let be the Euclidean distance from the i-th grid cell to the nearest body of water; This is the error term;

[0054] The calculation equation for the three-dimensional corrected cooling efficiency CE3D is as follows:

[0055]

[0056] In the formula, The three-dimensional corrected cooling efficiency for the lcz-th spatial partition; The regression coefficients of G3D on LST within the corresponding partition;

[0057] Before regression modeling, the variance inflation factor (VIF) of the covariates is calculated. When the VIF is greater than or equal to a preset threshold, the corresponding variable is removed to reduce multicollinearity interference. The above regression model is executed for each spatial partition to obtain the regression coefficients of G3D on LST within the corresponding partition. Since enhanced three-dimensional green space structures typically correspond to a decrease in surface temperature, [the following is taken]: The absolute value of the value is used as the three-dimensional corrected cooling efficiency for this spatial partition. .

[0058] Furthermore, in step S6, the minimum value of G3D within the target city area is obtained. As a baseline, it represents the worst pixel level of the 3D green space structure. For each spatial partition, the G3D of each raster cell within that partition is compared with... The difference between them is taken as the increment of the three-dimensional green space structure intensity, and multiplied by the CE3D corresponding to the zone to obtain the grid-scale three-dimensional corrected cooling capacity contribution value. The calculation equation for the zone-scale three-dimensional corrected cooling capacity CC3D is as follows:

[0059]

[0060] Divide the space into zones The urban-scale 3D corrected cooling capacity CC3D is obtained by weighting the corresponding number of grid cells by area. The calculation equation for the urban-scale 3D corrected cooling capacity CC3D is as follows:

[0061]

[0062] In the formula, The minimum G3D value within the target city area; This represents the number of grid cells within the corresponding spatial partition; This represents the number of raster cells within the lcz-th spatial partition. This represents the total number of grid cells involved in the calculation within the target city.

[0063] Furthermore, in step S7, to reflect the characteristics of population exposure changes over time and activity patterns, dynamic population raster data is acquired or constructed. Dynamic population raster data can include daytime population, nighttime population, weekday population, holiday population, commuter population, active population, and vulnerable populations such as the elderly and children. It can be generated by fusing data from mobile phone signaling, internet location service data, transportation data, POI functional distribution, building function type, residential population raster, and employment distribution. For each time scenario t, the average population density within the target city area is calculated. And the population of the i-th grid cell Divide by To obtain dynamic population exposure weights The dynamic population exposure scenarios include daytime, nighttime, weekdays, holidays, and sensitive populations; the calculation equation for the dynamic population exposure weight is:

[0064]

[0065] In the formula, The dynamic population exposure weight of the i-th grid cell in scene t; Let be the population of the i-th grid cell in scene t; Let be the average population density of the target city in scenario t.

[0066] Furthermore, in step S8, the three-dimensional corrected cooling capacity is combined with the dynamic population exposure weight to calculate the dynamic population exposure weighted cooling benefit. The calculation equation for the zonal-scale dynamic cooling benefit DCB is as follows:

[0067]

[0068] Divide the space into zones The dynamic cooling benefit of the target city at the city scale in scenario t is obtained by area weighting according to the corresponding number of grid cells. The calculation equation for the DCB (Dynamic Cooling Benefit) at the city scale is as follows:

[0069]

[0070] In the formula, The dynamic cooling of the lcz-type spatial partition in scenario t benefits; The target city benefits from dynamic cooling at the city scale in scenario t.

[0071] Further, step S9 specifically involves generating an urban green space cooling service function layer for the target city based on the CE3D, CC3D, and DCB results obtained from S5 to S8; the layer includes at least: a 3D corrected cooling efficiency CE3D spatial distribution layer, a 3D corrected cooling capacity CC3D spatial distribution layer, a dynamic cooling benefit DCB spatial distribution layer, a dynamic population exposure weight layer, a spatial zoning layer, and a zoning statistics table.

[0072] Step S10 specifically involves: checking the output results for valid value range, null values, water mask, partition coverage integrity, raster quantity consistency, 3D structural variable integrity, and dynamic population weight outlier. When abnormally high or low values ​​or missing partitions are found, the corresponding step is returned to check the input data, mask range, regression model, 3D structural weight, or dynamic population weight. After verification, the data source, target year, hot month or heat wave period, time scenario, resolution, regression parameters, CE3D / CC3D / DCB results, and statistical tables used in this inversion are written to the results database. For different years, different cities, or different time scenarios, S1 to S10 can be repeated.

[0073] Compared with the prior art, the beneficial effects of the present invention are:

[0074] This application no longer uses NDVI or green area to characterize the cooling effect of green space, but incorporates tree height, canopy coverage, building shadow and sky visibility factors into G3D, so that the cooling efficiency can reflect the differences in vegetation three-dimensional structure, shading process and radiation exchange.

[0075] This application establishes LST-G3D regression relationships within local climate zones or surface functional zones, and controls for covariates such as elevation, building height, and distance from water bodies, which can reduce estimation bias caused by the mixing of different spatial backgrounds.

[0076] This application calculates CC3D based on G3D intensity difference rather than a single NDVI difference, which can distinguish the differences in cooling capacity of different green space structures such as low grass, open canopy and closed canopy.

[0077] This application extends static population weights to dynamic population exposure weights such as daytime, nighttime, weekdays, holidays, or sensitive groups, which can identify mismatch areas between cooling service supply and population demand in different time scenarios.

[0078] This application outputs spatial layers such as CE3D, CC3D, and DCB, as well as zoning statistical tables, which can support the addition of green spaces at the block scale, canopy enhancement, shading optimization, and the delineation of priority areas for thermal risk management.

[0079] This application organizes thermomoon identification, image preprocessing, 3D structural feature construction, CE3D inversion, CC3D calculation, dynamic population weight construction, DCB calculation, and result output into a continuous process, reducing the differences in caliber and error propagation caused by manual decentralized processing. Attached Figure Description

[0080] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0081] The technical solutions of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0082] See Figure 1 The present invention adopts the following technical solution: a method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure.

[0083] Example:

[0084] This embodiment takes the central urban area of ​​a large subtropical city in southern China as the research scope, and carries out the inversion of the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure. It fully implements all the technical processes described in this invention, and the specific implementation steps are as follows:

[0085] Step 1: Acquire multi-source input data for the study area

[0086] Four images from Sentinel-2A / B satellite remote sensing data collected between July and August 2023 were selected, showing clear skies and no clouds. These images included 10m resolution optical band data and corresponding land surface temperature products for the respective time periods. Land cover data was derived from the city's 2022 land survey, with a resolution of 10m, covering land types such as cultivated land, orchards, forest land, grassland, construction land, and water areas. Building height data was obtained from the city's 2022 building census vector data, including the number of floors, base extent, and building height information for each building. The digital elevation model (DEM) used SRTM 30m resolution DEM data, resampled to 10m resolution via bilinear interpolation. Distances to water bodies were calculated using Euclidean distance analysis tools based on the extent of rivers, lakes, and reservoirs within the land cover, generating 10m resolution distance grids. Tree height data was obtained from 10m resolution tree height data retrieved from 2022 airborne LiDAR point clouds. High-resolution images cover the entire central urban area; canopy coverage is obtained at 10m resolution based on Sentinel-2 imagery using a pixel-based binary model, reflecting the horizontal coverage ratio of tree canopies; sky visibility factor is calculated using a 10m resolution viewpoint analysis algorithm based on building height, tree height, and DEM data, representing the sky visibility ratio of each pixel; building shadows are generated at 10m resolution using a shadow simulation algorithm based on the noon solar altitude and azimuth angles on July 15th, combined with building height and tree height data, and then the shadow ratio around each pixel is calculated using a 100m moving window; the dynamic population raster is a 100m resolution population raster generated by aggregating mobile signaling data from the city in 2023, including three typical scenarios: weekday daytime, weekday nighttime, and holiday daytime; the study area boundary is the administrative boundary vector data of the city's central urban area, clarifying the spatial scope of this inversion.

[0087] Step 2: Remote Sensing Image Preprocessing and Thermomoon Identification

[0088] Cloud masking was performed on Sentinel-2 imagery, and thick clouds and cloud image pixels were identified and removed using the QA60 band, retaining clear-sky surface pixels. Radiometric calibration was performed on the optical bands, converting DN values ​​into apparent radiance, and then atmospheric correction was completed using the Sen2Cor atmospheric correction tool to obtain surface reflectance data. Surface temperature products were scaled to convert quantized values ​​into true surface temperature in degrees Celsius. All imagery and spatial data were uniformly projected onto the UTM50N coordinate system, and all data were spatially cropped using the administrative boundaries of the central urban area. Data with different resolutions, such as DEM and population raster, were resampled to a 10m spatial resolution to ensure that the number of rows and columns of all raster data were consistent and their spatial positions were strictly aligned.

[0089] The Normalized Difference Vegetation Index (NDVI) was calculated based on the corrected red and near-infrared bands to characterize vegetation cover. Spatial statistics were performed on the monthly average land surface temperature from June to September 2023. The results showed that the average land surface temperature in July was the highest of the year, so July was determined as the target hot month for this inversion. The water body range was extracted based on the land cover data, and a uniform mask file was generated by combining the abnormal pixels with NDVI less than 0. All subsequent calculation steps removed such pixels to avoid interfering with the model fitting effect.

[0090] Step 3: Construct spatial partitions

[0091] Based on land cover data, two main categories of areas were extracted: built-up land and vegetation cover. Built-up land includes artificially paved surfaces such as residential, commercial, industrial, and road areas, while vegetation cover includes woodland, gardens, grasslands, and parkland. Within built-up land, an average building height of 10m was used as a threshold, and the average building height of each plot was statistically calculated through zoning. Areas with an average height less than 10m were classified as low-rise building areas, and areas with an average height of 10m or more were classified as mid- to high-rise building areas. Within vegetation cover areas, tree height data was used for further classification. Areas with tree height less than 2m were classified as non-tree vegetation areas, including lawns, low shrubs, and herbaceous flowers; areas with tree height greater than or equal to 2m were classified as tree vegetation areas. A canopy coverage of 40% was then used as a threshold, with areas having less than 40% canopy coverage classified as open canopy areas, and areas having 40% or more canopy coverage classified as closed canopy areas. Ultimately, five spatial zones were formed: non-tree vegetation zone, low-rise building zone, mid-to-high-rise building zone, open canopy zone, and closed canopy zone. The zones have continuous and complete boundaries with no overlap or omissions, covering the entire central urban area.

[0092] Step 4: Construct the comprehensive characterization value G3D of the three-dimensional green space structure

[0093] Five types of 10m resolution raster data were extracted: NDVI, tree height, canopy coverage, shade ratio, and sky visibility factor. Linear transformations were performed using the min-max normalization method to map the values ​​of each variable to the [0,1] interval. For the sky visibility factor, since a higher value indicates greater sky openness and weaker shading effect, the normalized value was subtracted from the normalized value after normalization, ensuring a positive correlation between the transformed value and the contribution of green space to cooling.

[0094] The entropy weight method is used to calculate the weight coefficients of five categories of variables. First, the information entropy of each category of variables is calculated. The smaller the information entropy, the higher the discrimination of the variable, and the greater the weight should be. The weight coefficients are then calculated. to The values ​​were 0.22, 0.24, 0.21, 0.18, and 0.15, respectively. Weighted summation was performed pixel-by-pixel according to the G3D calculation formula to obtain the spatial distribution data of the comprehensive three-dimensional green space structure characterization value (G3D) at a resolution of 10m for the study area. High G3D value areas were mainly concentrated in densely wooded areas such as urban parks and mountainous woodlands, while low value areas were concentrated in high-density built-up areas and bare land, consistent with the actual characteristics of urban green space distribution.

[0095] Step 5: Partition Inversion 3D Correction of Cooling Efficiency CE3D

[0096] Five spatial zones were used as independent calculation units, and multiple linear regression models were constructed for each. Each model used land surface temperature (LST) as the dependent variable and the three-dimensional green space structure comprehensive characterization value (G3D) as the core explanatory variable. At the same time, three environmental covariates, namely elevation, building height, and distance from water bodies, were introduced to control for the interference of topography, building form, and water cooling effect on land surface temperature, and to isolate the independent influence of G3D on LST.

[0097] Before performing regression modeling, a multicollinearity test was performed on the independent variables for each partition, and the variance inflation factor (VIF) of each independent variable was calculated. The test showed that the VIF values ​​of all independent variables within all partitions were less than 5, indicating no severe multicollinearity; therefore, all variables were retained for inclusion in the model. Ordinary least squares regression was then performed on each partition, and the regression coefficients corresponding to G3D regression were extracted. Because enhanced green space structures lower surface temperature, all zones... All values ​​are negative; the absolute value is taken to obtain the three-dimensional corrected cooling efficiency (CE3D) of the corresponding partition.

[0098] The results showed that the CE3D value was the highest in the closed canopy area, indicating that the cooling effect brought about by the improvement of the unit three-dimensional structure was the most significant; the open canopy area was the second highest; the CE3D value was the lowest in the non-tree vegetation area, which is consistent with the difference in cooling effect of different vegetation types; the CE3D value in the building area showed different numerical characteristics due to the interference of building form, which reflects the advantage of zonal regression in capturing spatial heterogeneity.

[0099] Step Six: Calculate the 3D corrected cooling capacity CC3D

[0100] First, the G3D values ​​of all valid rasters within the study area are calculated, and the minimum value is extracted. As the baseline state of the three-dimensional green space structure, it represents the pixel level with the worst green space structure and weakest cooling capacity in the study area. Its G3D value is calculated pixel by pixel and compared with the baseline value. The difference represents the increment of the three-dimensional green space structure of the grid relative to the worst baseline; multiplying this increment by the three-dimensional corrected cooling efficiency (CE3D) corresponding to the partition to which the grid belongs, we obtain the cooling capacity contribution value of each grid cell, that is, the additional cooling range that the grid can provide relative to the baseline state.

[0101] The average three-dimensional corrected cooling capacity of each zone is obtained by averaging the cooling capacity contribution values ​​of all grids within the five zones. The average cooling capacity of each zone is then multiplied by the number of effective grid cells in that zone, summed, and divided by the total number of effective grid cells in the city. This area-weighted calculation yields the overall 3D corrected cooling capacity (CC3D) for the central urban area. Finally, a 10m resolution CC3D spatial distribution layer is generated, visually representing the spatial differences in the cooling capacity of urban green spaces.

[0102] Step 7: Construct dynamic population exposure weights

[0103] Dynamic population exposure weighting was calculated for three typical scenarios: weekday daytime, weekday nighttime, and holiday daytime. First, the total population of the central urban area was calculated for each scenario, and then divided by the total area of ​​the central urban area to obtain the citywide average population density for that scenario. Count the population in the corresponding scene grid by grid. Divide by the city's average population density to obtain the dynamic population exposure weight of the grid in the corresponding scenario. .

[0104] The results show that during weekdays, high-weighted areas are concentrated in central business districts, industrial parks, and other areas with high employment concentrations, while low-weighted areas are concentrated in suburban residential areas and ecological green spaces. During weekday nights, high-weighted areas are concentrated in large residential areas, while low-weighted areas are concentrated in commercial and office areas. During holidays, high-weighted areas are concentrated around core business districts, large parks, and transportation hubs, consistent with the spatiotemporal patterns of urban population activities. All population data has undergone spatial aggregation and anonymization processing and does not contain any personally identifiable information.

[0105] Step 8: Calculate the dynamic cooling benefit DCB

[0106] By combining the 3D corrected cooling capacity with the dynamic population exposure weights of the three scenarios, the dynamic cooling benefit (DCB) is calculated. For each partition and each scenario, the cooling capacity contribution of each grid cell within the partition is multiplied by the population of the corresponding scenario. The sum of these values ​​is then divided by the product of the city's average population density for that scenario and the number of grid cells in the partition, yielding the average dynamic cooling benefit for that partition under the corresponding scenario. Then, multiply the dynamic cooling benefit of each zone by the number of grid cells in that zone, sum them up, and divide by the total number of grid cells in the city to obtain the overall dynamic cooling benefit of the central urban area in the corresponding scenario. .

[0107] The final result is a 10m resolution DCB spatial distribution layer for three scenarios. The results show that during weekdays, high-value DCB areas are mainly distributed in parks and street green spaces around business districts; these areas have strong cooling capabilities and are densely populated during the day. During weekdays at night, high-value DCB areas are concentrated in the green spaces within large residential areas, serving the nighttime population. On holidays, high-value DCB areas are concentrated in urban core parks and green spaces around commercial districts, corresponding to areas where leisure travelers congregate. The DCB layers clearly identify the degree of matching between cooling service supply and population demand, pinpointing key areas of supply-demand mismatch.

[0108] Step 9: Generate result layers and statistical tables

[0109] All calculation results were compiled to generate a standardized outcome system. For spatial layers, six raster layers were output: 3D corrected cooling efficiency (CE3D), 3D corrected cooling capacity (CC3D), dynamic cooling benefit (DCB) for three scenarios, spatial zoning, and dynamic population exposure weight. These were all in GeoTIFF format, maintaining a 10m resolution and UTM50N projection, and a unified invalid value encoding was applied. For statistical tables, a zoning statistics table and a city comprehensive index table were generated in CSV format. The zoning statistics table calculated the average, maximum, minimum, and standard deviation of CE3D, CC3D, and the three types of DCB, as well as the area percentage of each zoning category. The city comprehensive index table summarized the city's weighted average CE3D, average CC3D, and average DCB for each scenario, among other core indicators. Simultaneously, tiered thematic maps were created for outcome presentation and planning decision support.

[0110] Step 10: Result Verification and Storage Archiving

[0111] All results were subjected to quality checks according to the preset verification rules: First, the effective value range was checked to confirm that the core indicators such as LST, G3D, CE3D, CC3D, and DCB were all within a reasonable value range, and CE3D was positive with no abnormal extreme values; Second, null values ​​were checked to confirm that null values ​​were only distributed in the water body and boundary areas of the mask, with no continuous abnormal null values; Third, the water body mask was checked to confirm that all calculation results had removed water body pixels and had not been mixed with the low temperature effect of water bodies; Fourth, the zonal coverage integrity was checked to confirm that the total area of ​​the five types of zones was consistent with the area of ​​the study area, with no omissions or overlaps; Fifth, the raster number consistency was checked to confirm that the number of rows and columns of all raster layers was completely consistent with the total number of raster cells; Sixth, the 3D structural variable integrity was checked to confirm that there were no large-area missing G3D input variables and that the normalized results were all within the [0,1] interval; Seventh, the dynamic population weight outlier check was checked to confirm that the population weight distribution conformed to the urban activity pattern and had no extreme outliers.

[0112] After verification, all indicators met the quality requirements, and the results were reasonable and reliable. The metadata from this inversion and all result data were archived and stored together. The metadata records key information such as the data source used, target year, hottest month period, time scenario, spatial resolution, regression model parameters, and weight settings; the result data includes all spatial layers and statistical tables. This inversion fully verified the feasibility and stability of the method and can be extended to the evaluation of green space cooling service functions in other cities, other years, and other time periods.

[0113] Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for inverting the cooling service function of urban green spaces by integrating three-dimensional green space structure and dynamic population exposure, characterized in that: Includes the following steps: S1 acquires multi-source input data within the target city area, including at least remote sensing image data, land cover data, building height data, digital elevation model data, distance from water body data, tree height data, canopy coverage data, sky visibility factor data, building shadow data, dynamic population raster data, and study area boundary data. S2 preprocesses remote sensing image data, calculates the Normalized Difference Vegetation Index (NDVI) and Land Surface Temperature (LST), and identifies the hot month or heat wave period of the target year. S3 divides the target city into multiple spatial zones based on land cover data, building height data, and tree canopy cover structure; S4 unifies the spatial resolution and coordinate system of NDVI, tree height, canopy coverage, building shadow ratio and sky visibility factor, and after normalization, weighted fusion is used to construct the three-dimensional green space structure comprehensive characterization value G3D. S5 uses each spatial zone as an independent calculation unit to construct a multiple regression model with LST as the dependent variable, G3D as the core independent variable, and introduces terrain and building environment covariates to obtain the three-dimensional corrected cooling efficiency CE3D of each zone. S6 determines the minimum value of G3D within the target city area, G3Dmin, and calculates the three-dimensional corrected cooling capacity CC3D at the zonal and city scales by combining the intensity difference between CE3D of each zone and G3D of the grid. S7 constructs a grid-scale dynamic population exposure weight WPop for different population exposure scenarios; S8 combines three-dimensional corrected cooling capacity with dynamic population exposure weights to calculate dynamic cooling benefit DCB at the zonal and city scales. S9 generates spatial layers and statistical results for urban green space cooling service functions; S10 verifies, stores, and updates the output results.

2. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure as described in claim 1, characterized in that: In step S2, the preprocessing includes cloud masking, cloud shadow masking, radiometric calibration, scale factor conversion, projection conversion, spatial cropping, and resolution unification; masking is performed on water body pixels and pixels with NDVI less than 0; the identification method for the hot month or heat wave period is: synthesizing LST images month by month or day by day, statistically analyzing the average LST for each time period within the target city area, selecting the month with the highest average LST as the hot month, or selecting consecutive high-temperature days to constitute the heat wave period.

3. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure as described in claim 1, characterized in that: In step S3, the spatial partitioning includes five categories: non-tree vegetation area, low-rise building area, mid-to-high-rise building area, open canopy area, and closed canopy area. The partitioning rules are as follows: vegetation coverage area and building area are distinguished based on land cover data; low-rise building area and mid-to-high-rise building area are divided based on building height threshold within the building area; non-tree vegetation area and tree vegetation are distinguished within the vegetation coverage area; and open canopy area and closed canopy area are divided based on canopy coverage threshold within the tree vegetation area.

4. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure as described in claim 1, characterized in that: In step S4, the calculation equation for the comprehensive characterization value G3D of the three-dimensional green space structure is as follows: In the formula, This is the comprehensive three-dimensional green space structure representation value of the i-th grid cell; The normalized vegetation index; The normalized tree height; This represents the normalized canopy coverage. The normalized building shadow ratio; The normalized sky visibility factor; , , , , These are the weighting coefficients.

5. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure as described in claim 4, characterized in that: The normalization process uses the min-max normalization method to map each variable to the [0,1] interval; the weight coefficients are determined using any one of the following methods: equal weight, entropy weighting, expert weighting, regression coefficient normalization, or machine learning feature contribution.

6. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure as described in claim 1, characterized in that: In step S5, the calculation equation of the multiple regression model is: In the formula, For the intercept term; , , , These are the regression coefficients; The elevation of the i-th grid cell; Let i be the building height of the i-th grid cell; Let be the Euclidean distance from the i-th grid cell to the nearest body of water; This is the error term; The calculation equation for the three-dimensional corrected cooling efficiency CE3D is as follows: In the formula, The three-dimensional corrected cooling efficiency for the lcz-th spatial partition; The regression coefficients of G3D on LST within the corresponding partition; Before regression modeling, the variance inflation factor (VIF) of the covariates is calculated. When the VIF is greater than or equal to a preset threshold, the corresponding variable is removed to reduce multicollinearity interference.

7. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure as described in claim 1, characterized in that: In step S6, the calculation equation for the partition-scale three-dimensional corrected cooling capacity CC3D is as follows: The calculation equation for the three-dimensional corrected cooling capacity (CC3D) at the city scale is as follows: In the formula, The minimum G3D value within the target city area; This represents the number of grid cells within the corresponding spatial partition; This represents the number of raster cells within the lcz-th spatial partition. This represents the total number of grid cells involved in the calculation within the target city.

8. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure as described in claim 1, characterized in that: In step S7, the dynamic population exposure scenarios include daytime, nighttime, weekdays, holidays, and sensitive populations; the calculation equation for the dynamic population exposure weight is: In the formula, The dynamic population exposure weight of the i-th grid cell in scene t; Let be the population of the i-th grid cell in scene t; Let be the average population density of the target city in scenario t.

9. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure as described in claim 1, characterized in that: In step S8, the calculation equation for the DCB benefit of dynamic cooling at the zonal scale is as follows: The equation for calculating the dynamic cooling benefit DCB at the city scale is as follows: In the formula, The dynamic cooling of the lcz-type spatial partition in scenario t benefits; The target city benefits from dynamic cooling at the city scale in scenario t.

10. The method for inverting the urban green space cooling service function by integrating three-dimensional green space structure and dynamic population exposure according to claim 1, characterized in that: Step S9 specifically involves generating an urban green space cooling service function layer for the target city based on the CE3D, CC3D, and DCB results obtained from S5 to S8. The layer includes at least the following: a 3D corrected cooling efficiency CE3D spatial distribution layer, a 3D corrected cooling capacity CC3D spatial distribution layer, a dynamic cooling benefit DCB spatial distribution layer, a dynamic population exposure weight layer, a spatial zoning layer, and a zoning statistics table. Step S10 specifically involves: checking the output results for valid value range, null values, water mask, partition coverage integrity, raster quantity consistency, 3D structural variable integrity, and dynamic population weight outlier. When abnormally high or low values ​​or missing partitions are found, the corresponding step is returned to check the input data, mask range, regression model, 3D structural weight, or dynamic population weight. After verification, the data source, target year, hot month or heat wave period, time scenario, resolution, regression parameters, CE3D / CC3D / DCB results, and statistical tables used in this inversion are written to the results database. For different years, different cities, or different time scenarios, S1 to S10 can be repeated.

Images