A quantitative measurement method and system for the three-dimensional pattern of impermeable surfaces

Abstract

This invention discloses a quantitative measurement method and system for the three-dimensional pattern of impervious surfaces, comprising the following steps: reconstructing spatially continuous Landsat imagery; extracting the feature set with the highest accuracy in measuring impervious surfaces as the optimal feature set; estimating the abundance information of impervious surfaces using fully constrained least squares; gridding the city using a three-dimensional spatial grid; constructing quantitative indicators using building vector data and DEM data, and determining the weight of each indicator in the measurement of the three-dimensional pattern of impervious surfaces using the entropy method, calculating the development score of the three-dimensional pattern of impervious surfaces, and describing the three-dimensional pattern of impervious surfaces. This invention quantitatively describes the three-dimensional pattern of impervious surfaces, which helps to understand the three-dimensional spatial distribution of cities and provides a reference for urbanization patterns.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, specifically relating to a quantitative measurement method and system for the three-dimensional pattern of an impermeable surface. Background Technology

[0002] Cities are strategic core functional areas for socio-economic development. Against the backdrop of global urbanization, rapid urbanization has placed enormous pressure on the carrying capacity of urban ecological environments, and increasingly severe urban ecological problems are hindering sustainable urban development. How cities can achieve high-quality and effective development within the carrying capacity of their ecological environments, and how the ecological environment can better provide high-quality spatial carriers for urban development, are key issues for regional sustainable development.

[0003] Impermeable surfaces, comprised of building rooftops, roads, parking lots, plazas, and other impermeable landforms, are a key factor in measuring the level of urban development. Rapid urban development leads to the expansion of impermeable surfaces, a reduction in urban ecological land, and exacerbates the deterioration of the urban thermal environment. The urban thermal environment is influenced by both urban surface conditions and human activities, and is closely related to urban energy and resource consumption, air pollution, and the health of human living environments. It is a crucial component of the urban ecological environment and has become a hot topic in regional climate and environmental research. Studying impermeable surface patterns is of great significance for creating ecologically livable cities. Therefore, to meet the needs of global urbanization development strategies, it is urgent to focus on the spatial patterns of impermeable surfaces, thereby providing a scientific basis for sustainable urban development.

[0004] Current analyses of impermeable surface patterns are limited to two-dimensional space, neglecting the three-dimensional spatial information of impermeable surfaces. The height, spacing, and elevation differences of impermeable surfaces affect urban airflow and circulation, altering the city's internal energy balance and air movement. Three-dimensional spatial information about impermeable surfaces can more precisely reflect the city's internal energy balance. Therefore, a quantitative measurement method for the three-dimensional pattern of impermeable surfaces is urgently needed. Summary of the Invention

[0005] To address the problem that existing impermeable surface patterns are limited to two-dimensional spatial analysis and cannot accurately represent three-dimensional spatial information, this invention provides a quantitative measurement method and system for three-dimensional impermeable surface patterns, which can quantitatively describe urban three-dimensional space.

[0006] To achieve the above objectives, the technical solution of the present invention is a quantitative measurement method for the three-dimensional pattern of impermeable surfaces, comprising the following steps:

[0007] Step 1: Using the Landsat remote sensing dataset of the study area as the basic data source and Sentinel-2 imagery as a supplementary data source, reconstruct spatially continuous Landsat imagery.

[0008] Step 1.1: Preprocess the Landsat and Sentinel-2 images of the study area;

[0009] Step 1.2: Using a three-dimensional spatiotemporal interpolation method, the spatial and temporal information of pixels is fully utilized to fill in missing data and construct a spatially continuous Landsat image.

[0010] Step 2: Extract the feature set from the Landsat image reconstructed in Step 1, and select the feature set with the highest accuracy for impermeable surfaces as the optimal feature set.

[0011] Step 2.1: Use a selective integration strategy to integrate local feature selectors, combining the advantages of each local feature selector;

[0012] Step 2.2: For the features selected in Step 2.1, the Naive Bayes algorithm is used to evaluate the feature subset generated by the forward search of the sequence. The termination criterion is the decrease in the accuracy of the impermeable surface. Multiple feature sets are compared, and the feature subset with the highest accuracy of the impermeable surface is selected as the optimal feature set.

[0013] Step 3: Based on the optimal feature set, estimate the abundance information of impermeable surfaces using the fully constrained least squares algorithm;

[0014] Step 3.1: Based on the optimal feature set, extract the abundance of impermeable surfaces using the fully constrained least squares method;

[0015] Step 3.2: Samples are extracted from Landsat images using stratified random sampling, and the accuracy is evaluated using ten-fold cross-validation.

[0016] Step 4: Grid the study area using a three-dimensional spatial grid;

[0017] Step 5: Using building vector data and DEM data, construct four quantitative indicators for the three-dimensional spatial structure, three-dimensional morphology, three-dimensional compactness, and three-dimensional shape index of impermeable surfaces. Use the entropy method to determine the weight of each indicator in the measurement of the three-dimensional pattern of impermeable surfaces and calculate the development score of the three-dimensional pattern of impermeable surfaces.

[0018] Step 5.1: Use a geographic detector, geographic weighted regression, and clustering to describe the three-dimensional spatial structure of the impermeable surface;

[0019] Step 5.2: Use quantitative indicators to describe the three-dimensional morphology of the impermeable surface;

[0020] Step 5.3: Calculate the three-dimensional compactness of the impermeable surface;

[0021] Step 5.4, calculate the three-dimensional shape index of the impermeable surface;

[0022] Step 5.5: Use the entropy method to determine the weight of each index in the three-dimensional pattern measurement of impermeable surfaces, and calculate the development score of the three-dimensional pattern of impermeable surfaces.

[0023] Furthermore, in step 1.1, the Sentinel-2 image is resampled to the same spatial resolution as the Landsat image, and preprocessing operations such as geometric fine correction, BRDF correction, and spectral channel correction are performed on the Landsat image and the Sentinel-2 image.

[0024] Furthermore, in step 1.2, a multi-temporal image A = {A1, A2, ..., A...} is constructed using Landsat imagery and resampled adjacent temporal Sentinel-2 imagery. n}, where A i Let A represent an image, and n be the number of images. i The missing pixel P is estimated using M1 neighboring pixels within a 3×3×3 spatiotemporal window of pixel P. The 3D spatiotemporal interpolation method aims to find missing pixels. Make the distance-weighted least squares function Minimize, that is:

[0025]

[0026]

[0027] In the formula, P j This represents the neighboring pixels of the missing pixel P. W(d) represents the value of the missing pixel P. j ) is the weighting function, d j Let be the distance between the j-th neighboring pixel and the missing pixel P in the spatiotemporal domain. It represents the average distance between neighboring pixels and the missing pixel P in the spatiotemporal domain.

[0028] Furthermore, the final stability score of each local feature selector in step 2.1 is:

[0029]

[0030] In the formula, σ score σ is the score of the local feature selector. i Let k be the score for the i-th feature set, and k be the number of feature sets extracted by the local feature selector from the reconstructed Landsat image.

[0031] Based on the score of each local feature selector, all local feature selectors are sorted in descending order. A selective ensemble strategy is then used to select some local feature selectors for ensemble, that is, the top δ local feature selectors are selected for ensemble, and feature selectors with poor performance after the δ-th place are removed.

[0032] Furthermore, in step 2.2, the feature set is composed of features extracted by the top δ local feature selectors. N1 is the number of features in the feature set. Let the initial optimal feature set be... Suppose that the feature subset F currently contains p features has been selected. p The N1-p features that were not selected constitute the feature subset F. q (q=1,2,…,N1-p), the specific operation for selecting the optimal feature set is as follows: ① First, calculate F one by one. q Internal characteristics and F p Accuracy of the combined impermeable surface like but Select the feature set, the new feature set is ②Then calculate like but Select the feature set, the new feature set is ③ Repeat the above judgment and calculation process continuously until C calculated in the current iteration is less than C calculated in the previous iteration, at which point the iteration ends; ④ After the iteration ends, multiple feature subsets are obtained. The subset with the highest accuracy C for impermeable surfaces is selected as the optimal feature set F0. The accuracy of impermeable surfaces is described by weighted comprehensive analysis of four quantitative indicators: overall accuracy, precision, recall, and F-value, using the entropy method. The accuracy of impermeable surfaces in the feature subset can be used to determine whether the feature subset can effectively identify impermeable surfaces.

[0033] Furthermore, the FCLS model in step 3.1 includes two constraints: the sum of the abundance values ​​within a single pixel is 1, and the abundance value of each endmember is between 0 and 1, specifically expressed as follows:

[0034]

[0035]

[0036]

[0037] In the formula, r is the feature value of the sample pixel, r i (0≤i≤N2) is the endmember reflectivity, f i (0≤i≤N2) represents the endmember abundance, and N2 represents the number of endmembers;

[0038] The abundance value that satisfies the requirement is the product of reflectance and abundance of all endmembers in a pixel that differs from the original pixel the least.

[0039] Furthermore, in step 3.2, the evaluation indicators selected are root mean square error, mean absolute error, and systematic error. Among them, root mean square error and mean absolute error are used to evaluate the accuracy of impervious surface abundance estimation, and systematic error is used to quantify the deviation of impervious surface abundance estimation. The calculation method is as follows:

[0040]

[0041]

[0042]

[0043] In the formula, RMSE is the root mean square error, MAE is the squared absolute error, SE is the systematic error, and h i h represents the impermeability abundance information of the i-th pixel in the reconstructed Landsat image sample extracted using the method in step 3.1. i N3 is the number of samples, representing the impermeability abundance information of the i-th pixel in the same sample of a remote sensing image with a resolution higher than Landsat image, obtained using visual interpretation methods.

[0044] Furthermore, in step 4, the center of gravity of the study area is calculated using a centroid model, and the study area is meshed using a three-dimensional spatial coordinate system with the centroid as the origin. The centroid model calculation is as follows:

[0045]

[0046]

[0047] In the formula, X and Y represent the longitude and latitude of the city's center of gravity, respectively, and x i y i P represents the longitude and latitude of the i-th pixel, respectively. i N represents the value of the i-th cell, and N4 is the number of cells.

[0048] Furthermore, the three-dimensional spatial structure of the impermeable surface in step 5.1 includes point structures, single-center structures, single-center axis structures, multi-center cluster structures, and multi-center network structures. The point structure represents a single, compact block-like form of the impermeable surface; the single-center structure represents the outward extension of the core urban area along transportation routes; the single-center axis structure represents the outward development of the core urban area in a strip-like pattern; the multi-center cluster structure represents the block-like form of multi-center clusters after urban expansion; and the multi-center network structure represents the networked, multi-center, multi-axis form of the impermeable surface evolution. Geographic detectors, geographic weighted regression, and clustering are used to describe the three-dimensional spatial structure of the impermeable surface.

[0049] Furthermore, step 5.2 describes the three-dimensional morphology of the impermeable surface from both global and local perspectives. The global perspective analyzes the impermeable surface morphology within the grid based on building area ratio, floor area ratio, and maximum building surface area. The local perspective uses kernel density estimation and three-dimensional fractal dimension to intuitively express the three-dimensional morphological characteristics of the impermeable surface. The kernel density distribution is characterized by statistically analyzing the kernel density of building area ratio, shape index, floor area ratio, or clustering density to represent the three-dimensional building pattern distribution. The three-dimensional fractal dimension is as follows:

[0050] lgF(s)=α-δlgs (12)

[0051] In the formula, F(s) is the fractal dimension curve with s as the variable, lgF(s) represents the degree of irregularity of the three-dimensional shape of the impermeable surface, s is the scale, δ is the three-dimensional fractal dimension, and α is a constant.

[0052] Furthermore, the calculation method for the three-dimensional compactness of the impermeable surface in step 5.3 is as follows:

[0053]

[0054] In the formula, u i and u j Let d(i,j) be the volume of the impermeable surface within grids i and j, and d(i,j) be the geometric distance between grids i and j. k and u l Let the volume of the impermeable surface of the equivalent cone be the volume within grids k and l. M1 represents the geometric distance between grid k and grid l, M2 represents the number of impermeable surface grids in the study area, and N5 represents the number of grids in the equivalent cone.

[0055] Furthermore, the method for calculating the three-dimensional shape index of the impermeable surface in step 5.4 is as follows:

[0056]

[0057] In the formula, Shape 3d S is the ratio of the surface area of ​​an impermeable surface to the surface area of ​​a sphere of the same volume. i The impermeable surface area is calculated based on the number of pixels; P i The perimeter of the impermeable surface is calculated from the building's vector data; H i The height of the impermeable surface is obtained from the DEM data; g is the volume of the impermeable surface within the grid; and n6 is the number of grids.

[0058] Furthermore, in step 5.5, the entropy method is used to determine the weights of each index in the three-dimensional pattern measurement of impermeable surfaces, and the development score of the three-dimensional pattern of impermeable surfaces is calculated. The specific calculation formula is as follows:

[0059]

[0060]

[0061]

[0062] In the formula, M3 is the grid scale, N7 is the total number of indicators, and P ij For indicator X j The proportion of the i-th grid, X j w is any one of the four categories of indicators: three-dimensional spatial structure, three-dimensional morphology, three-dimensional compactness, and three-dimensional shape index. ij For indicator X ij The weight of the i-th grid, Y ij For indicator X j The overall development score in the i-th grid.

[0063] The present invention also provides a quantitative measurement system for the three-dimensional pattern of impermeable surfaces, for implementing the quantitative measurement method for the three-dimensional pattern of impermeable surfaces as described above.

[0064] Furthermore, it includes a processor and a memory, the memory being used to store program instructions, and the processor being used to call the stored instructions in the memory to execute a quantitative measurement method for the three-dimensional pattern of an impermeable surface as described above.

[0065] Compared with the prior art, the present invention has the following advantages:

[0066] This invention utilizes a three-dimensional spatiotemporal interpolation method to compensate for missing data and obtain spatially continuous remote sensing imagery. It employs a selective ensemble strategy to optimize feature sets and maximize the feature differences among different land cover types. Fully constrained least squares method is used to estimate impervious surface abundance information. By establishing a three-dimensional spatial grid, quantitative indicators describing the three-dimensional pattern of impervious surfaces are constructed based on building vector data and DEM data. This invention provides a quantitative description of the three-dimensional pattern of impervious surfaces, which helps in understanding the three-dimensional spatial distribution of cities and provides a reference for urbanization patterns. Attached Figure Description

[0067] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0068] Figure 1 This is a flowchart of the method according to an embodiment of the present invention. Detailed Implementation

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

[0070] Example 1

[0071] like Figure 1 As shown, this embodiment of the invention provides a quantitative measurement method for the three-dimensional pattern of impermeable surfaces, comprising the following steps:

[0072] Step 1: Using the Landsat remote sensing dataset of the study area as the primary data source and Sentinel-2 imagery as a supplementary data source, reconstruct spatially continuous Landsat imagery.

[0073] Step 1.1: Preprocess the Landsat and Sentinel-2 images of the study area.

[0074] In this embodiment, atmospherically corrected Landsat L2 level surface reflectance products are downloaded via USGS, and atmospherically corrected Sentinel-2 L2A level surface reflectance data are downloaded via ESA. Due to registration errors and spatial resolution inconsistencies between Landsat and Sentinel-2 images, to ensure spatial consistency, the Sentinel-2 images are resampled to the same spatial resolution as the Landsat images. Furthermore, geometrical correction, BRDF correction, and spectral channel correction are preprocessed on both Landsat and Sentinel-2 images.

[0075] Step 1.2: Using the three-dimensional spatiotemporal interpolation method, the spatial and temporal information of pixels is fully utilized to fill in the missing data and construct a spatially continuous Landsat image.

[0076] A multi-temporal image A = {A1, A2, ..., A...} was constructed using Landsat imagery and adjacent temporal Sentinel-2 imagery resampled to 30 meters. n}, where A i Let A represent an image, and n be the number of images. i The missing pixel P is estimated using M1 neighboring pixels within a 3×3×3 spatiotemporal window of pixel P. The 3D spatiotemporal interpolation method aims to find missing pixels. Make the distance-weighted least squares function Minimize, that is:

[0077]

[0078]

[0079] In the formula, P j This represents the neighboring pixels of the missing pixel P. W(d) represents the value of the missing pixel P. j ) is the weighting function, d j Let be the distance between the j-th neighboring pixel and the missing pixel P in the spatiotemporal domain. It represents the average distance between neighboring pixels and the missing pixel P in the spatiotemporal domain.

[0080] Step 2: Extract the feature set from the Landsat image reconstructed in Step 1, and select the feature set with high accuracy as the optimal feature set.

[0081] Step 2.1: Use a selective integration strategy to integrate local feature selectors, combining the advantages of each local feature selector.

[0082] This embodiment selects five local feature selectors: relevance, information gain, gain ratio, statistics, and consistency measure. Each local feature selector contains selectors based on multiple algorithms. The final stability score of each local feature selector is:

[0083]

[0084] In the formula, σ score σ is the score of the local feature selector. i The score is given for the i-th feature set, and k is the number of feature sets extracted by the local feature selector from the reconstructed Landsat image.

[0085] Based on the score of each local feature selector, all local feature selectors are sorted in descending order. A selective ensemble strategy is used to select some local feature selectors for ensemble. In this embodiment, the top 5 local feature selectors are selected for ensemble, and feature selectors with poor performance after the top 5 are removed.

[0086] Step 2.2: For the features selected in Step 2.1, the Naive Bayes algorithm is used to evaluate the feature subset generated by the forward search of the sequence. The termination criterion is the decrease in the accuracy of the impermeable surface. Multiple feature sets are compared, and the feature subset with the highest accuracy of the impermeable surface is selected as the optimal feature set.

[0087] Let the feature set consist of the features extracted by the top 5 local feature selectors. N1 is the number of features in the feature set. Let the initial optimal feature set be... Suppose that the feature subset F currently contains p features has been selected.p The N1-p features that were not selected constitute the feature subset F. q (q=1,2,…,N1-p), the specific operation for selecting the optimal feature set is as follows: ① First, calculate F one by one. q Internal characteristics and F p Accuracy of the combined impermeable surface like but Select the feature set, the new feature set is ②Then calculate like but Select the feature set, the new feature set is ③ Repeat the above judgment and calculation process continuously until C calculated in the current iteration is less than C calculated in the previous iteration, at which point the iteration ends; ④ After the iteration ends, multiple feature subsets are obtained. The subset with the highest accuracy C for impermeable surfaces is selected as the optimal feature set F0. The accuracy of impermeable surfaces is described using a weighted comprehensive method based on four quantitative indicators: overall accuracy, precision, recall, and F-value. The accuracy of impermeable surfaces within a feature subset determines whether that subset can effectively identify impermeable surfaces. The time complexity of this algorithm is ≤ N1(N1-1) / 2.

[0088] Step 3: Use the fully constrained least squares algorithm to estimate the abundance information of impermeable surfaces.

[0089] Step 3.1: Based on the optimal feature set, extract the abundance of impermeable surfaces using the fully constrained least squares (FLCS) method.

[0090] The FCLS model includes two constraints: the sum of abundance values ​​within a single pixel must be 1, and the abundance value of each endmember must be between 0 and 1. Specifically:

[0091]

[0092]

[0093]

[0094] In the formula, r is the feature value of the sample pixel, r i (0≤i≤N2) is the endmember reflectivity, f i (0≤i≤N2) represents the endmember abundance, and N2 represents the number of endmembers.

[0095] In remote sensing imagery, a pixel contains multiple endmembers. The abundance of an endmember is the percentage of that endmember's area within the pixel, representing the proportion of impermeable surfaces within that pixel. Assuming there are N pixels containing impermeable surface abundance values, the impermeable surface value required for subsequent impermeability measurement is the sum of the impermeable surface abundance values ​​of these N pixels. The abundance value that minimizes the difference between the product of reflectance and abundance of all endmembers within a pixel and the original pixel is the required abundance value.

[0096] Step 3.2: Samples are extracted from Landsat images using stratified random sampling, and the accuracy is evaluated using 10-fold cross-validation.

[0097] The evaluation metrics selected are root mean square error, mean absolute error, and systematic error. Root mean square error and mean absolute error are used to evaluate the accuracy of impervious surface abundance estimation, while systematic error is used to quantify the deviation in impervious surface abundance estimation. The calculation methods are shown below:

[0098]

[0099]

[0100]

[0101] In the formula, RMSE is the root mean square error, MAE is the squared absolute error, SE is the systematic error, and h i ′ represents the impermeability abundance information of the i-th pixel in the reconstructed Landsat image sample extracted using the method in step 3.1, k i N3 is the number of samples, representing the impermeability abundance information of the i-th pixel in the same sample of a remote sensing image with a resolution higher than Landsat image, obtained using visual interpretation methods.

[0102] Step 4: Grid the study area using a three-dimensional spatial grid.

[0103] The centroid of the study area was calculated using a centroid model, and the area was meshed using a 3D spatial mesh with the centroid as the origin of the coordinate system. The centroid model calculation is as follows:

[0104]

[0105]

[0106] In the formula, X and Y represent the longitude and latitude of the city's center of gravity, respectively, and x i y i P represents the longitude and latitude of the i-th pixel, respectively. i N represents the value of the i-th cell, and N4 is the number of cells.

[0107] Step 5: Using building vector data and DEM data, construct four quantitative indicators for the three-dimensional spatial structure, three-dimensional morphology, three-dimensional compactness, and three-dimensional shape index of the impermeable surface. Use the entropy method to determine the weight of each indicator in the measurement of the three-dimensional pattern of the impermeable surface, and calculate the development score of the three-dimensional pattern of the impermeable surface.

[0108] Step 5.1: Use a geographic detector, geographic weighted regression, and clustering to describe the three-dimensional spatial structure of the impermeable surface.

[0109] The three-dimensional spatial structure of impermeable surfaces includes point structures, single-center structures, single-center axis structures, multi-center cluster structures, and multi-center network structures. Point structures represent a single, compact block-like form of the impermeable surface; single-center structures represent the outward extension of the core urban area along transportation routes; single-center axis structures represent the outward strip-like development of the core urban area; multi-center cluster structures represent the block-like form of multi-center clusters after urban expansion; and multi-center network structures represent the networked, multi-center, multi-axis form of the impermeable surface evolution. Geographic detectors, geographic weighted regression, and clustering are used to describe the three-dimensional spatial structure of impermeable surfaces.

[0110] Step 5.2: Use quantitative indicators to describe the three-dimensional morphology of the impermeable surface.

[0111] This paper describes the three-dimensional morphology of impermeable surfaces from both global and local perspectives. The global perspective analyzes the morphology of impermeable surfaces within a grid based on data such as building area ratio, floor area ratio, and maximum building surface area. The local perspective uses kernel density estimation methods and three-dimensional fractal dimension to intuitively express the three-dimensional morphological characteristics of impermeable surfaces. Kernel density distribution can be used to characterize the distribution of three-dimensional building patterns. Specifically, it can be characterized by statistically analyzing the kernel density of building area ratio, shape index, floor area ratio, or clustering degree.

[0112] The three-dimensional fractal dimension is as follows:

[0113] lg F(s)=α-δlg s (12)

[0114] In the formula, F(s) is the fractal dimension curve with s as the variable, s is the scale, δ is the three-dimensional fractal dimension, and α is a constant with a value of 3.

[0115] lg F(s) represents the degree of irregularity in the three-dimensional shape of the impermeable surface. An increase in the degree of irregularity indicates that the increase in the three-dimensional volume of the impermeable surface is mainly due to expansion; a decrease in the degree of irregularity indicates that the increase in the three-dimensional volume of the impermeable surface is mainly due to edge filling; and no change in the degree of irregularity indicates that urban development has entered a stable period. The larger the volume of the impermeable surface within the grid, the larger the value of the three-dimensional fractal dimension.

[0116] Step 5.3: Calculate the three-dimensional compactness of the impermeable surface.

[0117] The three-dimensional compactness of impermeable surfaces reflects the degree of compactness and dispersion of impermeable surfaces in three-dimensional space, and is an important indicator for measuring the three-dimensional spatial agglomeration of urban land use. The calculation method for the three-dimensional compactness of impermeable surfaces is as follows:

[0118]

[0119] In the formula, u i and u j Let d(i,j) be the volume of the impermeable surface within grids i and j, and d(i,j) be the geometric distance between grids i and j. k and u l Let the volume of the impermeable surface of the equivalent cone be the volume within grids k and l. M1 represents the geometric distance between grid k and grid l, M2 represents the number of impermeable surface grids in the study area, and N5 represents the number of grids in the equivalent cone.

[0120] TCI 3d The larger the value, the more compact the impermeable surface space pattern.

[0121] Step 5.4: Calculate the three-dimensional shape index of the impermeable surface.

[0122] The three-dimensional shape index of an impermeable surface reflects its shape characteristics in three-dimensional space and characterizes its spatial complexity. The calculation method for the three-dimensional shape index of an impermeable surface is as follows:

[0123]

[0124] In the formula, Shape 3d S is the ratio of the surface area of ​​an impermeable surface to the surface area of ​​a sphere of the same volume. i The impermeable surface area is calculated based on the number of pixels; P i The perimeter of the impermeable surface is calculated from the building's vector data; H i The height of the impermeable surface is obtained from the DEM data; g is the volume of the impermeable surface within the grid; and N6 is the number of grids.

[0125] Step 5.5: Use the entropy method to determine the weight of each index in the three-dimensional pattern measurement of impermeable surfaces, and calculate the development score of the three-dimensional pattern of impermeable surfaces.

[0126] The specific calculation formula is as follows:

[0127]

[0128]

[0129]

[0130] In the formula, M3 is the grid scale, N7 is the total number of indicators, and P ij For indicator X j The proportion of the i-th grid, X j w is any one of the four categories of indicators: three-dimensional spatial structure, three-dimensional morphology, three-dimensional compactness, and three-dimensional shape index. ij For indicator X ij The weight of the i-th grid, Y ij For indicator X j The overall development score for the i-th grid. A higher score indicates that the three-dimensional layout of the area is more suitable for human habitation.

[0131] Example 2

[0132] Based on the same inventive concept, the present invention also provides a quantitative measurement system for the three-dimensional pattern of an impermeable surface, including a processor and a memory. The memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute the quantitative measurement method for the three-dimensional pattern of an impermeable surface as described above.

[0133] In specific implementation, the method proposed in the technical solution of this invention can be automatically executed by those skilled in the art using computer software technology. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of this invention and computer equipment including the computer program running the corresponding computer program, should also be within the protection scope of this invention.

[0134] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to replace them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A quantitative measurement method for the three-dimensional pattern of impermeable surfaces, characterized in that, Includes the following steps: Step 1: Using the Landsat remote sensing image dataset of the study area as the basic data source and Sentinel-2 imagery as a supplementary data source, reconstruct spatially continuous Landsat imagery. Step 2: Extract the feature set from the Landsat image reconstructed in Step 1, and select the feature set with the highest accuracy for impermeable surfaces as the optimal feature set. Step 3: Based on the optimal feature set, estimate the abundance information of impermeable surfaces using the fully constrained least squares algorithm; Step 4: Grid the study area using a three-dimensional spatial grid; Step 5: Using building vector data and DEM data, construct four quantitative indicators for the three-dimensional spatial structure, three-dimensional morphology, three-dimensional compactness, and three-dimensional shape index of impermeable surfaces. Use the entropy method to determine the weight of each indicator in the measurement of the three-dimensional pattern of impermeable surfaces and calculate the development score of the three-dimensional pattern of impermeable surfaces. The three-dimensional spatial structure of impermeable surfaces includes point structures, single-center structures, single-center axis structures, multi-center cluster structures, and multi-center network structures. Point structures represent a single, compact block-like form of impermeable surfaces; single-center structures represent the outward extension of the core urban area along transportation routes; single-center axis structures represent the outward development of the core urban area in a strip-like pattern; multi-center cluster structures represent the block-like form of multi-center clusters after urban expansion; and multi-center network structures represent the networked, multi-center, multi-axis form of impermeable surface evolution. Geographic detectors, geographic weighted regression, and clustering are used to describe the three-dimensional spatial structure of impermeable surfaces. The three-dimensional morphology of impermeable surfaces is described from both global and local perspectives. Globally, the morphology of impermeable surfaces within a grid is analyzed based on building area ratio, floor area ratio, and maximum building surface area. Locally, kernel density estimation and three-dimensional fractal dimension are used to intuitively express the three-dimensional morphological characteristics of impermeable surfaces. Kernel density distribution is characterized by statistically analyzing the kernel density of building area ratio, shape index, floor area ratio, or clustering degree to represent the distribution of the three-dimensional building pattern. The three-dimensional fractal dimension is as follows: lgF(s)=α-δlgs (12) In the formula, F(s) is the fractal dimension curve with s as the variable, lgF(s) represents the degree of irregularity of the three-dimensional shape of the impermeable surface, s is the scale, δ is the three-dimensional fractal dimension, and α is a constant; The calculation method for the three-dimensional compactness of impermeable surfaces is as follows: In the formula, u i and u j Let d(i,j) be the volume of the impermeable surface within grids i and j, and d(i,j) be the geometric distance between grids i and j. k and u l Let the volume of the impermeable surface of the equivalent cone be the volume within grids k and l. M1 represents the geometric distance between grid k and grid l, M2 represents the number of grids on the impermeable surface of the study area, and N5 represents the number of grids on the equivalent cone. The method for calculating the three-dimensional shape index of impermeable surfaces is as follows: In the formula, Shape 3d S is the ratio of the surface area of ​​an impermeable surface to the surface area of ​​a sphere of the same volume. i The impermeable surface area is calculated based on the number of pixels; P i The perimeter of the impermeable surface is calculated from the building's vector data; H i The height of the impermeable surface is obtained from the DEM data; g is the volume of the impermeable surface within the grid; and N6 is the number of grids.

2. The quantitative measurement method for the three-dimensional pattern of impermeable surfaces as described in claim 1, characterized in that: Step 1 includes the following steps: Step 1.1: Preprocess the Landsat and Sentinel-2 images of the study area; The Sentinel-2 image was resampled to the same spatial resolution as the Landsat image, and preprocessing operations such as geometric fine correction, BRDF correction, and spectral channel correction were performed on the Landsat image and the Sentinel-2 image. Step 1.2: Using a three-dimensional spatiotemporal interpolation method, the spatial and temporal information of pixels is fully utilized to fill in missing data and construct a spatially continuous Landsat image. A multi-temporal image A = {A1, A2, ..., A2} is constructed using Landsat imagery and resampled adjacent temporal Sentinel-2 imagery. n }, where A i Let A represent an image, and n be the number of images; for image A i The missing pixel P is estimated using M1 neighboring pixels within a 3×3×3 spatiotemporal window of pixel P; the three-dimensional spatiotemporal interpolation method aims to find the missing pixel. Make the distance-weighted least squares function Minimize, that is: In the formula, P j This represents the neighboring pixels of the missing pixel P. W(d) represents the value of the missing pixel P. j ) is the weighting function, d j Let be the distance between the j-th neighboring pixel and the missing pixel P in the spatiotemporal domain. It represents the average distance between neighboring pixels and the missing pixel P in the spatiotemporal domain.

3. The quantitative measurement method for the three-dimensional pattern of impermeable surfaces as described in claim 1, characterized in that: Step 2 includes the following steps: Step 2.1: Use a selective integration strategy to integrate local feature selectors, combining the advantages of each local feature selector; The final stability score for each local feature selector is: In the formula, σ score σ is the score of the local feature selector. i The score is given for the i-th feature set, and k is the number of feature sets extracted by the local feature selector on the reconstructed Landsat image. Based on the score of each local feature selector, all local feature selectors are sorted in descending order. A selective ensemble strategy is used to select some local feature selectors for ensemble, that is, the top δ local feature selectors are selected for ensemble, and feature selectors with poor performance after the δth place are removed. Step 2.2: For the features selected in Step 2.1, the Naive Bayes algorithm is used to evaluate the feature subset generated by the forward search of the sequence. The termination criterion is the decrease in the accuracy of the impermeable surface. Multiple feature sets are compared, and the feature subset with the highest accuracy of the impermeable surface is selected as the optimal feature set.

4. The quantitative measurement method for the three-dimensional pattern of impermeable surfaces as described in claim 3, characterized in that: In step 2.2, let the feature set consist of the features extracted by the top δ local feature selectors. N1 is the number of features in the feature set. Let the initial optimal feature set be... Suppose that the feature subset F currently contains p features has been selected. p The N1-p features that were not selected constitute the feature subset F. q For q = 1, 2, ..., N1-p, the specific steps for selecting the optimal feature set are as follows: ① First, calculate F one by one. q Internal characteristics and F p Accuracy of the combined impermeable surface like but Select the feature set, the new feature set is ②Then calculate like but Select the feature set, the new feature set is ③ Repeat the above judgment and calculation process continuously until C calculated in the current iteration is less than C calculated in the previous iteration, at which point the iteration ends; ④ After the iteration ends, multiple feature subsets F are obtained. p , The subset with the highest accuracy C for impermeable surfaces is selected as the optimal feature set F0. The accuracy of impermeable surfaces is described by weighted comprehensive analysis of four quantitative indicators: overall accuracy, precision, recall, and F-value, using the entropy method. The accuracy of impermeable surfaces in the feature subset is used to determine whether the feature subset can effectively identify impermeable surfaces.

5. The quantitative measurement method for the three-dimensional pattern of impermeable surfaces as described in claim 1, characterized in that: Step 3 includes the following steps: Step 3.1: Based on the optimal feature set, extract the abundance of impermeable surfaces using the fully constrained least squares method; The FCLS model includes two constraints: the sum of abundance values ​​within a single pixel must be 1, and the abundance value of each endmember must be between 0 and 1. Specifically: In the formula, r is the feature value of the sample pixel, r i f is the end-member reflectivity. i The term is the endmember abundance, where 0 ≤ i ≤ N2, and N2 is the number of endmembers; The abundance value that satisfies the requirement is the product of reflectance and abundance of all endmembers in a pixel that differs from the original pixel the least. Step 3.2: Samples are extracted from Landsat images using stratified random sampling, and the accuracy is evaluated using 10-fold cross-validation.

6. The quantitative measurement method for the three-dimensional pattern of impermeable surfaces as described in claim 5, characterized in that: In step 3.2, the evaluation indicators selected are root mean square error, mean absolute error, and systematic error. Root mean square error and mean absolute error are used to evaluate the accuracy of impervious surface abundance estimation, while systematic error is used to quantify the deviation in impervious surface abundance estimation. The calculation methods are shown below: In the formula, RMSE is the root mean square error, MAE is the squared absolute error, SE is the systematic error, and h i h represents the impermeability abundance information of the i-th pixel in the reconstructed Landsat image sample extracted using the method in step 3.

1. i N3 is the number of samples, representing the impermeability abundance information of the i-th pixel in the same sample of a remote sensing image with a resolution higher than Landsat image, obtained using visual interpretation methods.

7. The quantitative measurement method for the three-dimensional pattern of impermeable surfaces as described in claim 1, characterized in that: Step 4 uses a centroid model to calculate the centroid of the study area, and uses the centroid as the origin of the three-dimensional spatial coordinate system to mesh the study area using a three-dimensional spatial mesh. The centroid model calculation is as follows: In the formula, X and Y represent the longitude and latitude of the city's center of gravity, respectively, and x i y i P represents the longitude and latitude of the i-th pixel, respectively. i N represents the value of the i-th cell, and N4 is the number of cells.

8. The quantitative measurement method for the three-dimensional pattern of impermeable surfaces as described in claim 1, characterized in that: In step 5, the entropy method is used to determine the weights of each index in the three-dimensional pattern measurement of impermeable surfaces, and the development score of the three-dimensional pattern of impermeable surfaces is calculated. The specific calculation formula is as follows: In the formula, M3 is the grid scale, N7 is the total number of indicators, and P ij For indicator X j The proportion of the i-th grid, X j w is any one of the four categories of indicators: three-dimensional spatial structure, three-dimensional morphology, three-dimensional compactness, and three-dimensional shape index. ij For indicator X ij The weight of the i-th grid, Y ij For indicator X j The overall development score in the i-th grid; The higher the score, the more suitable the three-dimensional layout of the area is for human habitation.

9. A quantitative measurement system for the three-dimensional pattern of impermeable surfaces, characterized in that, It includes a processor and a memory, the memory being used to store program instructions, and the processor being used to call the program instructions in the memory to execute a quantitative measurement method for the three-dimensional pattern of an impermeable surface as described in any one of claims 1-8.

Images