Urban form ventilation potential assessment method and system based on rasterized spatial analysis and mathematical morphology

Abstract

The invention discloses an urban form ventilation potential assessment method and system based on rasterized spatial analysis and mathematical morphology, and belongs to the technical field of urban ventilation and thermal environment optimization. The method comprises the following steps: uniformly rasterizing a building and a road, and quickly calculating the windward area density by adopting a global rotation-block column extreme value-backfilling process; through multi-angle rotation and sliding window matching, pixel-level scoring of the downwind consistency of the road local is realized; extracting an opening width threshold value, a direction alignment degree and a shape factor of the semi-closed open space by using a multi-scale and anisotropic morphological opening operation; and superposing the indexes according to weights to generate a comprehensive ventilation potential map. According to the method, the calculation time consumption and memory occupation are effectively reduced while the calculation precision is ensured, the ventilation corridor, the high-resistance area and the local ventilation weakness can be accurately identified, and the method is suitable for urban planning, ventilation corridor special design, air quality simulation and microclimate optimization design.

Description

Technical Field

[0001] This invention belongs to the field of urban ventilation and thermal environment optimization technology, specifically relating to a method and system for assessing urban morphological ventilation potential based on gridded spatial analysis and mathematical morphology. Background Technology

[0002] Urban ventilation potential assessment is an important area of ​​research in urban climate and environment, playing a crucial role in mitigating the urban heat island effect, improving air quality, and enhancing human comfort. In urban planning and design, morphological parameters are typically used to indirectly assess the distribution characteristics of near-surface wind environments. Windward area density (FAD) is a core indicator for measuring spatial ventilation resistance, while the consistency between roads and wind direction, as well as the geometry of semi-enclosed open spaces, jointly influence the guidance and infiltration efficiency of local airflow. Currently, relevant technical guidelines and practices are gradually incorporating these morphological parameters into urban ventilation assessment systems to scientifically support the delineation of ventilation corridors and the optimization of spatial morphology.

[0003] Currently, there are two main technical paths for quantifying the aforementioned morphological parameters: the vector object method and the raster pixel method. The vector object method relies on geometric projection and topological operations, and its computational complexity increases sharply with the number of buildings and outline nodes, making it difficult to apply to large-scale, high-resolution scenes. Although the raster pixel method facilitates parallel computing through spatial discretization, it often requires pixel-by-pixel normal vector dot product to approximate the projection in FAD calculations, still facing bottlenecks in computational efficiency and memory usage. In addition, existing road downwind consistency assessments are mostly based on the global azimuth of the entire road segment, failing to accurately capture the wind-guiding effect of local straight segments in curved roads; and for semi-enclosed open spaces, commonly used indicators such as the sky visibility factor (SVF) are insufficient to effectively characterize the horizontal entrance width and wind direction alignment relationship, resulting in insufficient ability to identify ventilation weaknesses.

[0004] Therefore, existing technologies have not yet achieved a good balance between computational efficiency, accuracy, and the completeness of the indicator system. Developing an assessment method that can take into account rapid calculation, pixel-level fine evaluation, and multi-dimensional ventilation feature extraction has become an urgent need to improve the practical ability of urban ventilation potential analysis and support scientific planning. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention provides a rapid and robust quantitative assessment of wind resistance, road / open space downwind consistency, and ventilation characteristic parameters of semi-enclosed open spaces at multiple scales—regional, district, and street block—providing a methodological and data foundation for land spatial planning, ventilation corridor projects, microclimate-friendly urban design, and air quality pre-assessment. This invention proposes a method for assessing the ventilation potential of urban morphology based on gridded spatial analysis and mathematical morphology, comprising the following steps: Step 1: Convert the 3D vector data of urban buildings into a building height raster matrix. The resolution is r meters, non-building areas are assigned a value of 0, and building areas are assigned the building height. The building height raster matrix of the study area is divided into multiple m×m pixel blocks, and a block label matrix is ​​generated. ; arranging multi-segment vector data of the road into angle sequences After rotation, it is rasterized into a binary road matrix with grid size g. A set; Step 2: Calculate the windward area and windward area density of the grid blocks under any prevailing wind direction a, and obtain the windward area density distribution of each grid block under each prevailing wind direction in the study area. Step 3: Obtain the binary road matrix Medium pixel at all rotation angles Directional response intensity The average index is obtained by taking the average of the indices of the maximum and minimum rotation angles corresponding to the maximum directional response intensity value, and then obtaining the difference between the rotation direction corresponding to the average index and the prevailing wind direction 'a'. The weights of the directional response intensity corresponding to the average index are determined based on the weights of the directional response intensity of each pixel and the difference. The corresponding direction matching score is used to calculate the road tailwind consistency score; Step 4: Based on the height threshold Binarize the building height raster into an open area matrix. Employing a multi-scale structuring element set For open area matrix Perform morphological opening operations, where at least some structural elements are related to the prevailing wind direction. Aligned anisotropic structuring elements; Identification of semi-closed open space connected regions based on opening operation results. And extract its inlet width threshold. Orientation alignment and shape factor parameters Composite ventilation index ; Step 5: Normalize the windward area density and then invert it; simultaneously, calculate the road-wind consistency score and composite ventilation index. The data are normalized separately and then weighted and overlaid to generate a comprehensive ventilation potential map.

[0006] Further, step 2 includes: Taking due north as the wind direction, the building height raster matrix Bh and the block label matrix Lb are rotated counterclockwise by an angle a around the image center to obtain the rotation matrix: , In multi-wind direction analysis, wind direction angles are set as equally spaced discrete sets. Step length ; Extract the submatrix of the grid belonging to block grid i from the building height raster matrix. Find the maximum height from the height value sequence along the column direction. , This represents the maximum building height in the j-th column within the i-th grid. It is an indicator function, if If the value is equal to i, the value is 1; otherwise, it is 0. The column direction is the direction of the incoming wind. Calculated column sum Multiplying by the resolution r yields an approximate value for the windward area. Then, the windward area density of each grid segment is obtained. ; The windward area density of each grid segment Backfill to the positions corresponding to the block mesh in the rotation matrix. Then, by reversing the rotation, it is mapped back to the coordinate system before rotation to form a windward area density distribution map.

[0007] Calculate the maximum value sequence in the column direction At that time, an envelope correction strategy on adjacent columns is introduced, which uses the local maximum of the current column and its adjacent columns for smoothing to suppress isolated spikes caused by supertall buildings.

[0008] Furthermore, in step 2, the rotation operation uses bilinear interpolation, and the reverse rotation operation uses nearest neighbor interpolation.

[0009] Further, step 3 includes: for each rotation angle binary road matrix Using a 3×1 convolution kernel, sliding window convolution calculations are performed to obtain the pixel at all rotation angles. directional response intensity Composed of sequences , and sequence The index sequence of response intensities in all directions; select the indices of the first and last occurrences of the maximum directional response intensity from the index sequence. , ,by , The average value is used as the average index. The rotation angle corresponding to the average index As the local dominant direction of this pixel ;calculate With the target wind direction minimum angle difference According to the formula Calculate the road tailwind consistency score for each pixel; in, Weights for directional response strength Angle sequence Rotation angle interval in It is the minimum angle difference; The response intensity weighting coefficient in step 3 The value can be: when , ,when , ,when , ,when , .

[0010] Furthermore, in step 4, the anisotropic structural element is related to the prevailing wind direction. Aligned ellipses, their major and minor axis ratios can be calculated as follows: Alternatively, the settings can be based on the seasonal differences in prevailing winds given by the wind rose diagram; in, and These are the wind speed components in the parallel and vertical directions, respectively, used to reflect anisotropic ventilation sensitivity.

[0011] Furthermore, in step 4, an area threshold is applied to the opening operation result. To remove excessively large background open areas and excessively small noise, among which , The statistical quantiles were set based on the area of ​​the open zone in the study region.

[0012] Furthermore, the composite ventilation index in step 4 Calculated using the following formula: in, It is a connected region; 0-1 normalization relative to the entrance width; For orientation alignment, take a non-negative value. Indicates the dominant direction of the connected component. The prevailing wind direction used in the assessment The minimum included angle between them; The shape ventilation efficiency factor; α, β, and γ are weighting coefficients and AR represents the area of ​​the connected region. The solidity of the connected region is the ratio of the area of ​​the connected region to the area of ​​the circumscribed convex hull.

[0013] Furthermore, the weighted summation in step 5 is achieved using the following formula: in, , , , The weighting coefficients are λ1+λ2+λ3=1.

[0014] The present invention also relates to an urban morphology ventilation potential assessment system based on rasterized spatial analysis and mathematical morphology, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements an urban morphology ventilation potential assessment method based on rasterized spatial analysis and mathematical morphology.

[0015] Compared with the prior art, the beneficial effects of the present invention are: For calculating windward area density (FAD), this invention employs an innovative process of "global rotation—block-by-block column extremum—backfilling," transforming complex geometric projection calculations into efficient matrix column direction extremum and summation operations, achieving near-linear computational complexity. In large-scale, high-resolution scenarios, compared to traditional vector projection methods or pixel-by-pixel normal vector dot product methods, this invention can significantly reduce computation time by ≥50% and memory usage by ≥60% while maintaining a relative error of ≤5%, demonstrating excellent engineering deployment and parallel computing potential.

[0016] By employing a combination of multi-angle rotating rasterization and sliding window convolution, this invention achieves pixel-level local road orientation recognition and consistency scoring. This method effectively overcomes the misjudgment of curved roads in traditional road azimuth-based evaluations, accurately highlighting the wind-guiding value of long, straight sections. Furthermore, by introducing a first-to-last maximum averaging strategy and response intensity weights, it effectively suppresses the bimodal directional effect and noise interference, improving the robustness and continuity of the scoring results.

[0017] This invention innovatively applies multi-scale and anisotropic morphological opening operations to automatically and robustly extract key ventilation characteristic parameters such as the opening width threshold, orientation alignment, and shape factor of semi-enclosed spaces. This index system compensates for the shortcomings of the Sky Visibility Factor (SVF) in characterizing horizontal channels, enhances the ability to identify ventilation weaknesses such as U / Π-shaped enclosed clusters, and provides accurate data support for targeted renovation of open spaces.

[0018] This invention integrates the core assessment process into a unified grid computing framework, natively supporting multi-wind direction analysis, adjustable index weights for comprehensive ventilation potential mapping, and easy coupling with CFD, climate models, and planning decision-making systems. The method features a clear process and well-defined parameters, significantly improving the automation level and practicality of urban ventilation assessment results, and strongly supporting the transformation from scientific research to engineering practice. Attached Figure Description

[0019] Figure 1 This is a schematic diagram showing the rotation of the building height grid matrix according to the direction of the incoming wind. Figure 2 This is a schematic diagram of the calculation process for FAD; Figure 3 This is an example of a problem in calculating road-wind direction consistency based on vector data; Figure 4 This is the multi-scale extraction result for a semi-open space. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of this disclosure clearer, the technical solutions of the embodiments of this disclosure will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this disclosure. All other embodiments obtained by those skilled in the art based on the described embodiments of this disclosure without creative effort are within the scope of protection of this disclosure.

[0021] Example: A method and system for assessing urban morphological ventilation potential based on rasterized spatial analysis and mathematical morphology. This embodiment takes the main urban area of ​​a city (approximately 6.2km × 12.0km) as an example, but the implementation of this invention is not limited to this area.

[0022] This invention, within a unified grid computing framework, sequentially completes rapid calculation of windward area density, local scoring of road downwind consistency, extraction of ventilation indicators for semi-enclosed and open spaces, and outputs a comprehensive ventilation potential map, specifically including: Step 1: Data Preparation and Rasterization First, we collect 3D vector data of buildings, polyline vector data of roads, and prevailing wind direction information. Then, we convert the 3D vector data of buildings into a 2.5D building height raster matrix. The resolution is set to r = 2 meters. For each pixel (x, y), if the pixel is located inside a building, the element value is... The height is equal to the building's height; if a pixel is in a non-building area, it is assigned a value of 0. The prevailing wind direction is represented by angle 'a', which is the angle formed by rotating clockwise from due north to the prevailing wind direction. The 2.5D building height raster matrix of the study area is divided into M×N sub-grids, each sub-grid being a regular grid consisting of m×m pixels. Let m=100, and each sub-grid be approximately 200 meters × 200 meters in size. A sub-label matrix is ​​then generated. Each grid cell has a unique number. For each road segment vector data, sort the data by angle sequence. ( They are rotated and rasterized into multiple binary road matrices. The set, with a grid size g of 10 meters, where the angle The reference direction is due north, and the angle between the clockwise direction and due north is calculated. All data are uniformly based on the CGCS2000 coordinate system, and quality control measures such as geometric cleaning, outlier removal, and range alignment are performed.

[0023] Step 2: Rapid Calculation of Windward Area Density This step uses a "global rotation—block column extreme value—backfill" method to calculate the city's frontal area density (FAD) under any wind direction. It includes: Taking due north as the wind direction, the building height raster matrix Bh and the block label matrix Lb are rotated counterclockwise by an angle a around the image center to obtain the rotation matrix: , In multi-wind direction analysis, wind direction angles are set as equally spaced discrete sets. Step length For densely built-up areas with tall, thin buildings, sub-pixel resampling (i.e., super-resolution interpolation) is performed before rotation. The rotation interpolation uses a bilinear interpolation algorithm, which is a type of nearest-neighbor interpolation. Boundary regions are filled with 0, meaning that 0 represents a grid without buildings. During this operation, because the grid is a square grid, if the rotation angle is not an integer multiple of 90 degrees, it is impossible to ensure a one-to-one correspondence between the rotated and unrotated grids. Therefore, interpolation must be used to fill each grid.

[0024] Extract the submatrix of the grid belonging to block grid i from the building height raster matrix. Find the maximum height from the height value sequence along the column direction. , This represents the maximum building height in the j-th column within the i-th grid. It is an indicator function, if If the value is equal to i, the value is 1; otherwise, it is 0. The column direction is the direction of the incoming wind. Calculated column sum Multiplying by the resolution r yields an approximate value for the windward area. Then, the windward area density of each grid segment is obtained. .

[0025] The windward area density of each grid segment Backfill to the positions corresponding to the block mesh in the rotation matrix. Then, by inverse rotation, the coordinates are mapped back to the original coordinate system (preferably nearest neighbors to avoid cumulative secondary interpolation) to form the FAD distribution map corresponding to the prevailing wind direction 'a'. This method controls the computational complexity to the O(MN) level, with a relative error of no more than 5%, and supports multi-threaded parallel computation.

[0026] To handle isolated spikes caused by data anomalies or supertall buildings, this step introduces a "peripheral envelope on adjacent columns" correction strategy. This involves using the local maximum values ​​of the current column and its adjacent columns to correct the spikes. Smoothing is performed to effectively suppress unreasonable extreme values.

[0027] Note: Summing the maximum values ​​in a column is equivalent to a discrete approximation of the height envelope by line integration along the wind direction. The core parameters for FAD calculation are resolution *r* and block size *m*, which have a clear impact on performance and accuracy. Decreasing *r* (e.g., from 2m to 1m) can improve accuracy and reduce error, but will significantly increase computation time and memory usage; increasing *m* (e.g., from 100 to 150 pixels) can reduce computation time, but in areas with drastic fluctuations in building height, it will introduce slightly higher errors due to the increased approximation. When the relative error is sufficiently small (r≤2m) and the rotational interpolation is controllable, denote the upper bound of the relative error. ( )satisfy (Right now (related to r / typical projected side length and rotation interpolation error), where Typical projected side length, This represents the equivalent error of rotational interpolation. In practice, a trade-off must be struck between these parameters depending on the scale of the study area and computational resources.

[0028] Step 3: Local scoring of road tailwind consistency Pixel-level road direction consistency evaluation is achieved through multi-angle rotation and sliding window matching.

[0029] Use 3×1 convolution kernel (i.e., a window of 3 consecutive pixels), in each rotated binary road matrix Slide upwards to obtain the local triple pixel directional response intensity at the pixel being slid to. Obtain pixels at all rotation angles directional response intensity Composed of sequences , and sequence An index sequence of response intensities in all directions. , The larger the value, the higher the value at that rotation angle. The stronger the evidence, the closer the pixel is to the one that forms a "three-pixel straight line"; Select the indices from the index sequence where the maximum directional response intensity occurs for the first and last time. , ,by , The average value is used as the average index, i.e. The rotation angle corresponding to the average index As a local direction estimate ,calculate With the prevailing wind direction minimum angle difference .

[0030] Finally, the tailwind severity of each pixel is weighted and scored by combining the directional response intensity and the minimum angular difference, resulting in a road tailwind score map. The scoring formula is as follows: In the formula, To address the response intensity weighting factor, this embodiment introduces a weight. Weighting, i.e. , , , ; Represents a pixel; It is the minimum angle difference; The directional response intensity of the pixel; This is the angle index corresponding to the maximum directional response intensity.

[0031] Note that the key parameters for road tailwind scoring are the angular interval Δθ and the grid size g. Reducing Δθ (e.g., by 10°) can improve angular resolution, but will increase computation by approximately 50%. The grid size g affects spatial accuracy and computational efficiency; g=5m is suitable for detailed analysis in densely roaded old urban areas, while g=15m is more suitable for rapid screening over a large area. If necessary, a 3×3 median filter can be applied to dense road network areas to suppress isolated noise. The angle index corresponding to the maximum directional response intensity in step 3... The determination of the angle using the average of the first and last maximum value positions is to alleviate the directional ambiguity caused by the bimodal distribution (i.e., to suppress the influence of the bimodal effect on score calculation), which can improve the stability of angle estimation. Furthermore, the rotation and convolution operations in this step are both performed based on global matrix transformations, avoiding object-level geometric topology calculations.

[0032] Step 4: Extraction of ventilation index for semi-enclosed open spaces The system automatically identifies "semi-enclosed but ventilated" open spaces (such as courtyards enclosed by buildings but with entrances) from the building height grid, quantifies their ventilation capacity, and outputs a comprehensive index, ESI.

[0033] First, pixels in the building height raster matrix with heights below the height threshold h0 are designated as passable open zones and assigned a value of 0; otherwise, they are assigned a value of 1, thus obtaining the open zone matrix. , rice.

[0034] Employing multi-scale structural elements For open area matrix Perform morphological opening operations; Among them, structural elements It includes anisotropic and isotropic structural elements. The anisotropic structural elements adopt a major axis that is aligned with the wind direction. An ellipse aligned with its minor axis perpendicular to the wind direction, with a major-to-minor axis ratio of Alternatively, the setting can be dynamically adjusted based on the seasonal prevailing wind differences given by the wind rose diagram (e.g., 1.5~2.0). and These represent the wind speed components in the parallel and vertical wind directions, respectively, used to reflect anisotropic ventilation sensitivity. The size of the minor axis b is selected from the set... , For diameter, The meter is taken every 5 to 10 meters, and the isotropic structural element is a disk.

[0035] Identifying Connected Regions in Semi-Enclosed Open Spaces Based on Opening Operation Results The result of the opening operation at each scale Perform area filtering, retaining only areas within the specified range. The semi-enclosed open space connecting areas within the scope; , The statistical quantile can be set according to the area of ​​the open zone in the study region. The area is set to the 95th percentile of all open areas to effectively exclude large background open areas such as rivers, lakes, and large squares. Set it to the 5th percentile of all open area areas to effectively eliminate excessively small noise.

[0036] By tracing the connecting areas of semi-enclosed open spaces Changes in connectivity at different scales, when The scale first led to connected regions When disconnected from external connectivity, determine its inlet width threshold. , exist The condition is defined as "disconnected from the outside." Here, the criterion for "disconnected from the outside" is that the connected region becomes an independent connected component in the binary image, no longer connected to the image boundary or other external regions, and the region orientation of the connected region is calculated. (The direction of this region refers to the direction of the major axis of the circumscribed ellipse.) , (Refers to the lengths of the major and minor axes of the circumscribed ellipse of the connected region, respectively), Solidity (the solidity of the connected region, i.e., the ratio of the area of ​​the connected region to the area of ​​the circumscribed convex hull), AR (the area of ​​the connected region), and other geometric parameters, which are used to construct a composite ventilation index. : In the formula, (α=β=γ=1 / 3); 0-1 normalization relative to the entrance width; Let be the orientation alignment, and take the non-negative value, where Indicates the region direction of the connected region The prevailing wind direction used in the assessment The minimum included angle between them; It is the shape ventilation efficiency factor.

[0037] Engineering Explanation and Supplement: The opening operation is equivalent to a "doorway filter," which automatically assigns an inlet width threshold when the inlet is narrower than the structural element diameter, thus disconnecting it from the outside. Anisotropic structural elements introduce wind direction weights to enhance sensitivity to differences in monsoon dominance. ESI focuses on horizontal channels and inlet geometry, while SVF focuses on upper hemisphere openness; combining the two provides a more comprehensive understanding.

[0038] Step 5: Comprehensive Ventilation Potential Mapping and Output The comprehensive ventilation potential is generated by weighting and superimposing normalized and weighted windward area density, road-tailwind consistency score, and composite ventilation index. .

[0039] Specifically, the ventilation resistance index is obtained by normalizing FAD and then inverting it. The road tailwind consistency index, after normalization, is: The normalized composite ventilation index is . norm(X)=(X−min) / (max−min); Using weights The final ventilation potential value is obtained by linearly weighting (λ1+λ2+λ3=1). .

[0040] in The values ​​are inverted to reflect that "the smaller the FAD, the lower the resistance and the greater the ventilation potential." The results are output in GeoTIFF format and support PNG visualization and CSV table export for easy urban planning applications.

[0041] To objectively verify the comprehensive performance of the method of this invention, a typical area containing 18,511 buildings (grid resolution r=2m) was selected in the study area. The method of this invention (denoted as the "rotation-column extremum-backfill" method) was systematically compared with two existing mainstream methods (i.e., the vector projection method based on building-by-building geometric projection and topological intersection calculation, and the grid method based on pixel-by-pixel normal vector dot product approximate projection). The experiment used multi-dimensional quantitative indicators for evaluation: in terms of efficiency, the total computation time and peak memory usage were examined; in terms of accuracy, the median (P50) and 95th percentile (P95) of the relative error were calculated by comparing the calculation results of each method with the high-precision geometric analytical solution; in terms of stability, the variance of the results under different combinations of resolution (r) and block size (m) was analyzed for evaluation.

[0042] Table 1 shows the performance comparison results of the three methods in typical scenarios. Experimental results show that the method of this invention has significant advantages in computational efficiency and memory usage, reducing the total time consumption by approximately 63% and memory usage by approximately 71% compared to the vector projection method. In terms of computational accuracy, the median error of the method of this invention is comparable to that of the pixel-by-pixel method, while the 95% quantile error is lower, demonstrating superior error control and stability. The underlying mechanism for its efficiency improvement lies in: reducing the complex three-dimensional projection problem to the column direction extremum calculation of a two-dimensional grid through global rotation; reducing the data scan range in a single operation through block processing; and allowing the global rotation results to be reused in multi-wind direction analysis, avoiding redundant calculations.

[0043] To further verify the necessity of each technical component in this invention, the following ablation experiments were conducted: Elimination of the "average of first and last maximum values" strategy in road scoring: If this strategy is canceled and the direction of the first maximum value is taken directly, a falsely high consistency score will be generated at the bends, resulting in an increase in the overall score noise level of about 8-12%.

[0044] Impact of sliding window kernel size: The convolution kernel in the road consistency analysis was changed from the recommended 3×1 to 5×1. Although the response of long straight road segments was enhanced, the ability to detect continuity in narrow streets and alleys decreased, resulting in breaks in effective ventilation paths.

[0045] The influence of isotropic morphological structural elements: In the extraction of indices in semi-enclosed spaces, if only isotropic disk structural elements are used and anisotropic elliptical kernels are abandoned, the sensitivity of the obtained ESI index to the difference in the dominant monsoon direction will be significantly reduced, weakening its assessment ability under complex wind field conditions.

[0046] The above ablation experiments demonstrate that the key technical components used in this invention, such as the "average of the first and last maximum values", the 3×1 sliding window kernel, and anisotropic morphological operations, play an irreplaceable role in ensuring the accuracy, continuity, and sensitivity to wind direction of the evaluation results, and together constitute the advanced nature and robustness of the technical solution of this invention.

[0047] The present invention also provides an urban morphology ventilation potential assessment system based on rasterized spatial analysis and mathematical morphology, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements an urban morphology ventilation potential assessment method based on rasterized spatial analysis and mathematical morphology.

[0048] This invention's system is deployed on a computer with an 8-core, 16-thread or higher CPU and at least 32GB of RAM. An optional GPU with at least 8GB of video memory can be added to accelerate computation. It supports Windows / Linux operating systems, is developed using Python or MATLAB, and uses NumPy and Numba / CuPy for matrix operations, and GDAL / OGR for reading and writing geographic data.

[0049] The system adopts a data flow-driven sequential execution mode. After startup, each core module executes in turn: starting with data input and preprocessing, through rasterization, FAD calculation, road consistency analysis, and morphological index extraction, finally completing the comprehensive evaluation and visualization output. To enhance the system's flexibility and integrability, the system also provides a set of RESTful application programming interfaces (APIs), allowing users to call core functional modules or dynamically set parameters through specific HTTP requests.

[0050] Based on the above assessment results, the following targeted applications can be carried out in urban planning and management: Ventilation corridor identification and construction: high consistency in identifying road downwind routes ( High value area) and low building ventilation resistance ( In areas with low air quality, these sections will be prioritized for connection and integration as potential ventilation corridors, forming key pathways to promote urban air circulation.

[0051] Optimize ventilation in residential clusters: Identify ventilation bottlenecks in clusters with low ESI (Energy Stability Index) values ​​for semi-enclosed open spaces. Prioritize planning interventions in these areas, such as adding connecting passageways, widening entrances, or creating linear green spaces to improve local wind conditions.

[0052] Development intensity and form guidance: High density of building windward area ( For neighborhoods with high (high) values ​​but low overall ventilation potential (V), stricter controls should be implemented. It is recommended to limit the height and volume of new buildings, or to improve spatial openness through optimized layout to reduce their obstruction of urban airflow.

[0053] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of this invention. It should be understood that the above descriptions are merely specific embodiments of this invention and are not intended to limit the scope of protection of this invention. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A method for assessing the ventilation potential of urban morphology based on rasterized spatial analysis and mathematical morphology, characterized in that, Includes the following steps: Step 1: Convert the 3D vector data of urban buildings into a building height raster matrix. The resolution is r meters, non-building areas are assigned a value of 0, and building areas are assigned the building height. The building height raster matrix of the study area is divided into multiple m×m pixel blocks, and a block label matrix is ​​generated. ; arranging multi-segment vector data of the road into angle sequences After rotation, it is rasterized into a binary road matrix with grid size g. A set; Step 2: Calculate the windward area and windward area density of the grid blocks under any prevailing wind direction a, and obtain the windward area density distribution of each grid block under each prevailing wind direction in the study area. Step 3: Obtain the binary road matrix Medium pixel at all rotation angles Directional response intensity The average index is obtained by taking the average of the indices of the maximum and minimum rotation angles corresponding to the maximum directional response intensity value, and then obtaining the difference between the rotation direction corresponding to the average index and the prevailing wind direction 'a'. The weights of the directional response intensity corresponding to the average index are determined based on the weights of the directional response intensity of each pixel and the difference. The corresponding direction matching score is used to calculate the road tailwind consistency score; Step 4: Based on the height threshold Binarize the building height raster into an open area matrix. Employing a multi-scale structuring element set For open area matrix Perform morphological opening operations, where at least some structural elements are related to the prevailing wind direction. Aligned anisotropic structural elements; Identifying Connected Regions in Semi-Enclosed Open Spaces Based on Opening Operation Results And extract its inlet width threshold. Orientation alignment and shape factor parameters Composite ventilation index ; Step 5: Normalize the windward area density and then invert it; simultaneously, calculate the road-wind consistency score and composite ventilation index. The data are normalized separately and then weighted and overlaid to generate a comprehensive ventilation potential map.

2. The method according to claim 1, characterized in that, Step 2 includes: Taking due north as the wind direction, the building height raster matrix Bh and the block label matrix Lb are rotated counterclockwise by an angle a around the image center to obtain the rotation matrix: , In multi-wind direction analysis, wind direction angles are set as equally spaced discrete sets. Step length ; Extract the submatrix of the grid belonging to block grid i from the building height raster matrix. Find the maximum height from the height value sequence along the column direction. , This represents the maximum building height in the j-th column within the i-th grid. It is an indicator function, if If the value is equal to i, the value is 1; otherwise, it is 0. The column direction is the direction of the incoming wind. Calculated column sum Multiplying by the resolution r yields an approximate value for the windward area. Then, the windward area density of each grid segment is obtained. ; The windward area density of each grid segment Backfill to the positions corresponding to the block mesh in the rotation matrix. Then, by reversing the rotation, it is mapped back to the coordinate system before rotation to form a windward area density distribution map.

3. The method according to claim 2, characterized in that, Calculate the maximum value sequence in the column direction At that time, an envelope correction strategy on adjacent columns is introduced, which uses the local maximum of the current column and its adjacent columns for smoothing to suppress isolated spikes caused by supertall buildings.

4. The method according to claim 2, characterized in that, In step 2, the rotation operation uses bilinear interpolation, and the reverse rotation operation uses nearest neighbor interpolation.

5. The method according to claim 1, characterized in that, Step 3 includes: for each rotation angle binary road matrix Using a 3×1 convolution kernel, sliding window convolution calculations are performed to obtain the pixel at all rotation angles. directional response intensity Composed of sequences , and sequence The index sequence of response intensities in all directions; select the indices of the first and last occurrences of the maximum directional response intensity from the index sequence. , ,by , The average value is used as the average index. The rotation angle corresponding to the average index As the local dominant direction of this pixel ;calculate With the target wind direction minimum angle difference According to the formula Calculate the road tailwind consistency score for each pixel; in, Weights for directional response strength Angle sequence Rotation angle interval in It is the minimum angle difference; The response intensity weighting coefficient in step 3 The value can be: when , ,when , ,when , ,when , .

6. The method according to claim 1, characterized in that, In step 4, the anisotropic structural element is related to the target wind direction. Aligned ellipses, their major and minor axis ratios can be calculated as follows: Alternatively, the settings can be based on the seasonal differences in prevailing winds given by the wind rose diagram; in, and These are the wind speed components in the parallel and vertical directions, respectively, used to reflect anisotropic ventilation sensitivity.

7. The method according to claim 1, characterized in that, In step 4, an area threshold is applied to the opening operation result. To remove excessively large background open areas and excessively small noise, among which , The statistical quantiles were set based on the area of ​​the open zone in the study region.

8. The method according to claim 1, characterized in that, The composite ventilation index in step 4 Calculated using the following formula: in, It is a connected region; 0-1 normalization relative to the entrance width; For orientation alignment, take a non-negative value. Indicates the dominant direction of the connected component. The prevailing wind direction used in the assessment The minimum included angle between them; The shape ventilation efficiency factor; α, β, and γ are weighting coefficients and AR represents the area of ​​the connected region. The solidity of the connected region is the ratio of the area of ​​the connected region to the area of ​​the circumscribed convex hull.

9. The method according to claim 1, characterized in that, The weighted summation in step 5 is achieved using the following formula: in, , , , The weighting coefficients are λ1+λ2+λ3=1.

10. A system for assessing the ventilation potential of urban morphology based on rasterized spatial analysis and mathematical morphology, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the method as described in any one of claims 1 to 9.

Images