Urban underlying surface type extraction method based on high-resolution image

Abstract

The invention discloses an urban underlying surface type extraction method based on a high-resolution remote sensing image, and aims to solve the problems that a traditional segmentation method depends on artificial parameters, shadow interference is significant and classification model precision is insufficient. A non-dominated sorting genetic algorithm is adopted to realize automatic selection of an optimal segmentation scale, a multi-scale segmentation method is adopted to segment a high-resolution remote sensing image according to the optimal segmentation scale, a JM distance optimized by a Gaussian kernel function is introduced to dynamically merge object optimization segmentation results, and object-oriented data is obtained; and on the basis of the object-oriented data, performing feature optimization and data dimension reduction by adopting non-negative matrix factorization (NMF), and inputting the processed data into a Light GBM base classifier for image classification so as to complete urban underlying surface type extraction. According to the method, the difficult problem of complex city ground feature form recognition is effectively solved, the shadow interference influence is eliminated, and high-precision automatic city underlying surface type extraction is realized.

Description

Technical Field

[0001] The present invention relates to the field of surface cover classification, and in particular to a method for extracting urban underlying surface types based on high-resolution images. Background Art

[0002] As global warming continues, the frequency of urban heat events is increasing. These events create heat loads, increase morbidity and energy consumption, and impact human health. Accurately capturing urban thermal zones facilitates timely disaster warnings and the precise implementation of disaster prevention measures. During high-temperature events, urban underlying surface type is a core parameter in thermal environment simulations, and its accurate extraction directly impacts the effectiveness of disaster warnings and urban heat island mitigation measures.

[0003] Currently, urban underlying surface types are extracted using low- to medium-resolution remote sensing data such as Landsat and MODIS. While these data have large-scale observation capabilities, their insufficient spatial resolution allows them to only simplify cities into a single feature category and cannot identify detailed underlying surface types such as buildings, roads, and green spaces. While high-resolution remote sensing images can capture the geometric texture characteristics of objects, they face the following technical bottlenecks: Limitations of segmentation methods: Existing segmentation methods (such as trial and error method and index optimization method) rely on manually set parameters, resulting in highly subjective segmentation results, blurred boundaries, and difficulty in adapting to complex urban landforms.

[0004] Shadow interference problem: Dynamic changes such as building shadows and vegetation shadows will significantly affect the accuracy of underlying surface recognition, and traditional methods lack a dynamic correction mechanism.

[0005] Classification model defects: Single classifiers (such as support vector machines and random forests) have limited ability to express high-dimensional and heterogeneous features. Summary of the Invention

[0006] The purpose of this invention is to provide a method for extracting urban underlying surface types based on high-resolution images, which achieves high-precision and automated urban underlying surface classification using an integrated classification model by dynamically optimizing the segmentation scale and improving the feature optimization strategy.

[0007] To achieve the above object, the present invention provides a method for extracting urban underlying surface types based on high-resolution remote sensing images, comprising: Determine the optimal segmentation scale of high-resolution remote sensing images, segment the high-resolution remote sensing images according to the optimal segmentation scale through multi-scale segmentation methods, and optimize the segmentation results with probabilistic distance to obtain object-oriented data; Based on the object-oriented data, multiple training subsets are generated by random sampling without replacement. Non-negative matrix factorization (NMF) is used to perform feature optimization and data dimensionality reduction on each training subset to obtain the processed training subset. The processed training subset is input into the LightGBM base classifier for image classification to complete the extraction of urban underlying surface types.

[0008] Furthermore, the method for obtaining object-oriented data specifically includes: A comprehensive optimization objective function is constructed by using a homogeneity evaluation index and a heterogeneity evaluation index, and the optimal segmentation scale of the high-resolution remote sensing image is determined according to the comprehensive optimization objective function by using the RMNE method; Segmenting the high-resolution remote sensing image according to the optimal segmentation scale using a multi-scale segmentation method to generate an initial segmentation result, which includes multiple initial image objects; Based on the initial segmentation result, the JM distance between two adjacent initial image objects is calculated. When the JM distance is less than the target threshold T, the two adjacent initial image objects are merged. Otherwise, the segmentation result is retained to obtain the final segmentation result. According to the final segmentation results, object-oriented data is constructed.

[0009] Furthermore, the RMNE method determines the optimal segmentation scale of the high-resolution remote sensing image according to the comprehensive optimization objective function, including: Generate a set of candidate segmentation scales for high-resolution remote sensing images , each candidate segmentation scale Corresponding to different spatial texture information entropy and spectral heterogeneity index , the spatial texture information entropy As an evaluation index of internal homogeneity of objects, the calculation formula is: (1) In formula (1), The first Probability distribution of texture-like features; The spectral heterogeneity index As an evaluation index of heterogeneity between adjacent objects, the calculation formula is: (2) In formula (2), For the current object Band mean, is the mean value of the corresponding band of adjacent objects, is the spectral standard deviation, is the weight coefficient and satisfies ; According to formula (1), the internal homogeneity evaluation index of each candidate segmentation scale is obtained, and according to formula (2), the heterogeneity evaluation index between adjacent objects corresponding to each candidate segmentation scale is obtained; Constructing a comprehensive optimization objective function based on the homogeneity evaluation index and the heterogeneity evaluation index , the comprehensive optimization objective function The calculation formula is: (3) In formula (3), is the maximum value of spatial texture information entropy in the candidate segmentation scale set, is the minimum value of spectral heterogeneity in the set of candidate segmentation scales, is the weight coefficient and satisfies ; The non-dominated sorting genetic algorithm is used to perform multi-objective optimization on each candidate segmentation scale to minimize the comprehensive optimization objective function. As the goal, generate the Pareto front solution set; According to the comprehensive optimization objective function The value of is selected from the Pareto front solution set to make the comprehensive optimization objective function Minimum candidate segmentation scale as the optimal segmentation scale.

[0010] Furthermore, the JM distance is calculated using a Gaussian kernel function for probability density estimation, and the calculation formula is: (4) In formula (4), is the Gaussian kernel bandwidth parameter, for In the The difference in feature dimensions; is the probability density function based on Gaussian distribution.

[0011] Furthermore, the method for determining the target threshold T includes: Based on the initial segmentation results, an initial object set is constructed and the JM distance between two adjacent initial image objects is calculated to form a sample set. ; The sample set is constructed using Gaussian kernel density estimation method The probability distribution model is used to determine the peak point of the probability distribution model by maximum likelihood estimation. ; Peak point As the center, generate candidate threshold intervals : (5) In formula (5), For the sample set The standard deviation of is the adaptive coefficient, The value of is based on the sample set Dynamic adjustment of skewness; Through grid search in the candidate threshold interval Step length Calculate the value for each candidate The proportion of merged objects corresponding to the value; Select the candidate that meets the screening conditions from all candidate T values Value as target threshold , the screening conditions include: the proportion of merged objects of the target threshold T is within the statistical confidence interval of the natural merging rate of the sample set D; and the average JM distance between objects merged through the target threshold T is the largest; When there is no valid target threshold that meets the screening conditions in the candidate threshold interval When , the sample set is reconstructed by Box-Cox transformation The distribution is recalculated; when the recalculated candidate threshold interval still has no valid target threshold When as the target threshold T.

[0012] Furthermore, the method for extracting the urban underlying surface type specifically includes: Construct a feature training set based on object-oriented data, and set the number L of LightGBM base classifiers and the number K of feature subsets; randomly dividing the feature training set into K mutually exclusive training subsets according to the number K of the feature subsets, resampling each training subset, and applying a non-negative matrix factorization (NMF) method to the object-oriented data in the resampled training subset to generate rotation coefficients; The rotation coefficients are reorganized in the order corresponding to the corresponding training subsets and then multiplied with the feature training set to construct an enhanced feature set; The enhanced feature set is trained in parallel by L LightGBM base classifiers to obtain L trained LightGBM base classifiers; For the image objects to be classified, the L trained LightGBM base classifier models are used to perform predictions respectively to obtain L preliminary classification results; A weighted vote is performed on the L preliminary classification results based on the classification confidence, and the urban underlying surface type with the highest vote is output as the final extracted urban underlying surface type.

[0013] Furthermore, the division of the training subsets satisfies: (6) In formula (6), is the total number of features in the feature training set, is the number of preset feature subsets, To obey the uniform distribution , generate K mutually exclusive training subsets through the Fisher-Yates shuffle algorithm, each training subset contains Features.

[0014] Furthermore, the object-oriented data in the training subset is subjected to a non-negative matrix factorization (NMF) method to generate rotation coefficients, including: Construct a non-negative data matrix for the object-oriented data in each training subset , the non-negative data matrix is ​​transformed into Decompose into basis matrix Coefficient matrix ; Define the objective function of non-negative matrix factorization (NMF): (7) In formula (7), is a non-negative data matrix, is the basis matrix, is the coefficient matrix, is the sparse regularization coefficient, is the number of features in the training subset, is the feature dimension of the training subset, is the preset implicit feature dimension; The objective function is minimized by alternating least squares method, and the basis matrix is ​​iteratively updated. and coefficient matrix ; The error of non-negative matrix decomposition after each iteration is measured by calculating KL divergence. When the error change is less than the preset threshold or the maximum number of iterations is reached, the decomposition of the non-negative data matrix is ​​terminated and the final converged basis matrix is ​​extracted. as the rotation factor.

[0015] The present invention also protects an electronic device comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the method for extracting urban underlying surface types based on high-resolution remote sensing images is implemented.

[0016] The present invention also protects a computer-readable storage medium storing a computer program, wherein the computer program enables a computer to execute the above-mentioned method for extracting urban underlying surface types based on high-resolution remote sensing images.

[0017] The beneficial effects of the present invention are: By dynamically balancing the homogeneity and heterogeneity within the object through a multi-objective optimization method and combining it with a non-dominated sorting genetic algorithm to generate the optimal segmentation scale, the problem of traditional segmentation scale relying on manual experience is solved, and the error of segmentation boundaries is greatly reduced. By dynamically generating the judgment threshold of merging objects based on the JM distance of the Gaussian kernel, the segmentation results are optimized, the boundaries of homogeneous objects are effectively distinguished, and the merging misjudgment rate is reduced. It is especially suitable for complex scenes at the junction of building shadows and permeable surfaces. In the improved rotation forest framework, the non-negative matrix factorization (NMF) method is used instead of the traditional PCA method to construct an enhanced feature set. This method preserves the physical meaning of features during dimensionality reduction, avoids feature confusion caused by principal component direction shift, and improves the interpretability of classification features. Combining the non-negative matrix factorization (NMF) method and the LightGBM feature importance screening method to compress feature dimensions, effectively reduce computational complexity, and enable large-scale image processing. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 It is a flow chart of the method for extracting urban underlying surface types based on high-resolution remote sensing images of the present invention; Figure 2 It is a schematic diagram of the process of obtaining object-oriented data according to the present invention; Figure 3 Schematic diagram of the process of determining the optimal segmentation scale of the present invention; Figure 4 Schematic diagram of the process of determining the custom threshold T of the present invention; Figure 5 It is a schematic diagram of the process of image classification of the present invention; Figure 6 This is a flowchart of the non-negative matrix factorization (NMF) algorithm. DETAILED DESCRIPTION

[0019] The following will provide a clear and complete description of the technical method of the present invention in conjunction with the accompanying drawings. Obviously, the embodiments described are only some of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making any creative efforts shall fall within the scope of protection of the present invention.

[0020] Example 1 See also Figure 1 ,The urban underlying surface type extraction method based on high-resolution remote sensing images includes: Determine the optimal segmentation scale of high-resolution remote sensing images, segment the high-resolution remote sensing images according to the optimal segmentation scale through multi-scale segmentation methods, and optimize the segmentation results with probabilistic distance to obtain object-oriented data; Based on the object-oriented data, multiple training subsets are generated by random sampling without replacement. Non-negative matrix factorization (NMF) is used to perform feature optimization and data dimensionality reduction on each training subset to obtain the processed training subset. The processed training subset is input into the LightGBM base classifier for image classification to complete the extraction of urban underlying surface types.

[0021] Specifically, high-resolution remote sensing images from Sentinel-2 are used to cover the target urban area. Sensor noise is eliminated through radiometric calibration and converted into reflectance / radiance values. Atmospheric correction is performed using the FLAASH model to eliminate atmospheric interference. Ground control points (GCPs) are used to achieve multi-temporal image geometric registration, and the NDVI vegetation index is calculated to assist in ground object classification.

[0022] See also Figure 2 As a preferred embodiment of the present invention, the method for obtaining object-oriented data specifically includes: A comprehensive optimization objective function is constructed by using a homogeneity evaluation index and a heterogeneity evaluation index, and the optimal segmentation scale of the high-resolution remote sensing image is determined according to the comprehensive optimization objective function by using the RMNE method; Segmenting the high-resolution remote sensing image according to the optimal segmentation scale using a multi-scale segmentation method to generate an initial segmentation result, which includes multiple initial image objects; Based on the initial segmentation result, the JM distance between two adjacent initial image objects is calculated. When the JM distance is less than the target threshold T, the two adjacent initial image objects are merged. Otherwise, the segmentation result is retained to obtain the final segmentation result. According to the final segmentation results, object-oriented data is constructed.

[0023] See also Figure 4 As a preferred embodiment of the present invention, the RMNE method determines the optimal segmentation scale of the high-resolution remote sensing image according to the comprehensive optimization objective function, including: Generate a set of candidate segmentation scales for high-resolution remote sensing images , each candidate segmentation scale Corresponding to different spatial texture information entropy and spectral heterogeneity index , the spatial texture information entropy As an evaluation index of internal homogeneity of objects, the calculation formula is: (1) In formula (1), The first Probability distribution of texture-like features; The spectral heterogeneity index As an evaluation index of heterogeneity between adjacent objects, the calculation formula is: (2) In formula (2), For the current object Band mean, is the mean value of the corresponding band of adjacent objects, is the spectral standard deviation, is the weight coefficient and satisfies ; According to formula (1), the internal homogeneity evaluation index of each candidate segmentation scale is obtained, and according to formula (2), the heterogeneity evaluation index between adjacent objects corresponding to each candidate segmentation scale is obtained; Specifically, a set of candidate segmentation scales is set based on urban characteristics (densely built-up areas, farmland patches, water areas, etc.), and the scale range covers typical urban elements.

[0024] Spatial texture information entropy calculate: Calculate 8 types of texture features (contrast, correlation, energy, homogeneity, etc.) based on GLCM (Gray Level Co-occurrence Matrix). For each segmented object, count the 8 types of texture feature values ​​and normalize them into probability distribution. ,The spatial texture information entropy is calculated by formula (1), which reflects the complexity of the internal texture of the object. The higher the entropy value, the worse the homogeneity.

[0025] Spectral heterogeneity index calculate: Select four key bands B2 (blue), B3 (green), B4 (red), and B8 (near infrared) that are sensitive to vegetation and buildings, calculate the absolute value of the mean difference between the current object and the adjacent objects in the selected bands, and calculate the standard deviation of the adjacent objects in band c. Calculate the spectral heterogeneity index using formula (2) and adjust the weights , balancing spectral continuity and diversity.

[0026] Constructing a comprehensive optimization objective function based on the homogeneity evaluation index and the heterogeneity evaluation index , the comprehensive optimization objective function The calculation formula is: (3) In formula (3), is the maximum value of spatial texture information entropy in the candidate segmentation scale set, is the minimum value of spectral heterogeneity in the set of candidate segmentation scales, is the weight coefficient and satisfies ; Specifically, the weights are determined through expert experience or experimental data Specific values, balancing homogeneity and heterogeneity.

[0027] The non-dominated sorting genetic algorithm is used to perform multi-objective optimization on each candidate segmentation scale to minimize the comprehensive optimization objective function. As the goal, generate the Pareto front solution set; According to the comprehensive optimization objective function The value of is selected from the Pareto front solution set to make the comprehensive optimization objective function Minimum candidate segmentation scale as the optimal segmentation scale.

[0028] Specifically, the segmentation scale The encoding is a real chromosome, the population size is set to 50, the maximum number of iterations is 100 generations, and the multi-objective optimization function is constructed using formula (3). Through dynamic weighting strategy adjustment, the weight of texture features is increased in dense urban areas; the spectral weight is enhanced in suburban farmland; congestion distance sorting is used to retain the top 10% of non-dominated solutions as candidate optimal scales; the non-dominated sorting genetic algorithm (NSGA-II) is used to generate a non-dominated solution set, which can ensure that the segmentation scale is optimal among multiple objectives. The classification accuracy at different scales is evaluated through cross-validation to ensure robustness.

[0029] As a preferred embodiment of the present invention, the JM distance is calculated using a Gaussian kernel function to perform probability density estimation, and the calculation formula is: (4) In formula (4), is the Gaussian kernel bandwidth parameter, for In the The difference in feature dimensions; is the probability density function based on Gaussian distribution.

[0030] See also Figure 5 As a preferred embodiment of the present invention, the method for determining the target threshold value T includes: Based on the initial segmentation results, an initial object set is constructed and the JM distance between two adjacent initial image objects is calculated to form a sample set. ; The sample set is constructed using Gaussian kernel density estimation method The probability distribution model is used to determine the peak point of the probability distribution model by maximum likelihood estimation. ; Peak point As the center, generate candidate threshold intervals : (5) In formula (5), For the sample set The standard deviation of is the adaptive coefficient, The value of is based on the sample set Dynamic adjustment of skewness; Through grid search in the candidate threshold interval Step length Traversing candidates Calculate the value for each candidate The proportion of merged objects corresponding to the value; Select the candidate that meets the screening conditions from all candidate T values Value as target threshold , the screening conditions include: the proportion of merged objects of the target threshold T is within the statistical confidence interval of the natural merging rate of the sample set D; and the average JM distance between objects merged through the target threshold T is the largest; When there is no target threshold that meets the screening conditions in the candidate threshold range When , the sample set is reconstructed by Box-Cox transformation distribution after recalculation; when the recalculated candidate threshold interval still has no target threshold When as the target threshold T.

[0031] Specifically, for the initial object set after segmentation, traverse all adjacent object pairs and calculate based on formula (4): Distance, store the JM distance values ​​of all adjacent object pairs into a collection , forming a length of Sample set , Gaussian kernel function is used to perform probability density estimation to construct the sample set The probability distribution model is used to determine the peak point of the probability distribution model by maximum likelihood estimation. , generate the interval according to formula (5) , with a step length Traverse all candidate T values ​​in the interval, and for each candidate T value, merge the values ​​that satisfy Adjacent objects, calculate the proportion of merged objects , based on the normal distribution assumption, calculate the confidence interval of the natural merger rate, and select the candidate that meets the screening conditions from all candidate T values Value as target threshold , the screening conditions include: The proportion of merged objects at the target threshold T is within the statistical confidence interval of the natural merging rate of the sample set D; The average JM distance between objects after merging by the target threshold T is the largest; If there are multiple candidate T values ​​that meet the screening conditions, the one closest to the peak point is selected first. The candidate threshold T is used as the target threshold ; When the candidate threshold interval No valid target threshold within When the Box-Cox transformation is triggered, the sample set Perform Box-Cox transformation to obtain a new sample set , re-execute steps S1.3.2-S1.3.5 to generate a new candidate threshold interval. If the recalculated candidate threshold interval still does not have a valid target threshold When , it is forced to select As the target threshold T, and output a warning log.

[0032] See also Figure 3 As a preferred embodiment of the present invention, the method for extracting the urban underlying surface type specifically includes: Construct a feature training set based on object-oriented data, and set the number L of LightGBM base classifiers and the number K of feature subsets; randomly dividing the feature training set into K mutually exclusive training subsets according to the number K of the feature subsets, resampling each training subset, and applying a non-negative matrix factorization (NMF) method to the object-oriented data in the resampled training subset to generate rotation coefficients; The rotation coefficients are reorganized in the order corresponding to the corresponding training subsets and then multiplied with the feature training set to construct an enhanced feature set; The enhanced feature set is trained in parallel by L LightGBM base classifiers to obtain L trained LightGBM base classifiers; For the image objects to be classified, the L trained LightGBM base classifier models are used to perform predictions respectively to obtain L preliminary classification results; A weighted vote is performed on the L preliminary classification results based on the classification confidence, and the urban underlying surface type with the highest vote is output as the final extracted urban underlying surface type.

[0033] As a preferred embodiment of the present invention, the division of the training subsets satisfies: (6) In formula (6), is the total number of features in the feature training set, is the number of preset feature subsets, To obtain a random perturbation term that obeys uniform distribution, K mutually exclusive training subsets are generated by Fisher-Yates shuffling algorithm, each training subset contains Features.

[0034] Specifically, the Fisher-Yates algorithm is used to randomly rearrange the feature training set to generate K mutually exclusive training subsets, including: Generate random indices starting from the end of the feature training set , swap the current feature with the random position feature; repeat this process until all features are processed, ensuring that each feature subset contains features, among which is a uniformly distributed random perturbation term used to balance the feature overlap between subsets; when the total number of features in the feature training set is Cannot be When divisible, by Term compensation remainder, e.g. hour, , Value is an integer and the number of features in the training subset differs .

[0035] See also Figure 6 As a preferred embodiment of the present invention, the object-oriented data in the training subset is subjected to a non-negative matrix factorization (NMF) method to generate rotation coefficients, including: Construct a non-negative data matrix for the object-oriented data in each training subset , the non-negative data matrix is ​​transformed into Decompose into basis matrix Coefficient matrix ; Define the objective function of non-negative matrix factorization (NMF): (7) In formula (7), is a non-negative data matrix, is the basis matrix, is the coefficient matrix, is the sparse regularization coefficient, is the number of features in the training subset, is the feature dimension of the training subset, is the preset implicit feature dimension; The objective function is minimized by alternating least squares method, and the basis matrix is ​​iteratively updated. and coefficient matrix ; The error of non-negative matrix decomposition after each iteration is measured by calculating KL divergence. When the error change is less than the preset threshold or the maximum number of iterations is reached, the decomposition of the non-negative data matrix is ​​terminated and the final converged basis matrix is ​​extracted. as the rotation factor.

[0036] Example 2 This embodiment provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, the method for extracting urban underlying surface types based on high-resolution remote sensing images described in Example 1 is implemented.

[0037] This embodiment also provides a computer-readable storage medium having a computer program stored thereon. When the program is executed by a processor, the method for extracting urban underlying surface types based on high-resolution remote sensing images described in the first embodiment is implemented.

[0038] Example 3 For a typical urban area in Henan Province, the 2024 Sentinel-2 high-resolution remote sensing imagery was used, covering the main urban area of ​​Zhengzhou and surrounding suburbs.

[0039] Image data preprocessing includes: using the Sen2Cor tool for atmospheric correction to eliminate aerosol interference and obtain surface reflectance data; achieving multi-phase image registration through ground control points (GCPs) with a root mean square error (RMSE) of <0.5 pixels; and calculating the normalized difference vegetation index (NDVI) to assist in distinguishing between vegetation-covered areas and non-vegetated areas.

[0040] Based on the urban characteristics of Henan Province, including densely built-up areas, farmland patches, and water areas, the optimal segmentation scale is determined: Select the segmentation scale: Settings pixels, covering typical urban features (e.g., average building size 20-50m, farmland patches 50-100m); Compute texture and spectral features: Calculate spatial texture information entropy: Calculate 8 types of texture features (contrast, correlation, etc.) based on GLCM, and calculate the entropy value after normalization. For example, the texture complexity of the building area is high (entropy value ), suburban farmland texture is simple (entropy value ); Calculate spectral heterogeneity indicators: Select the B2 (blue), B3 (green), B4 (red), and B8 (near-infrared) bands and calculate the absolute value of the mean difference between adjacent objects. For example, the spectra of roads and built-up areas are significantly different, while the spectra of farmland and green land are similar. A comprehensive optimization objective function is constructed based on the homogeneity evaluation index and the heterogeneity evaluation index, and the weights of the homogeneity evaluation index and the heterogeneity evaluation index are adjusted through a dynamic weighting strategy. In dense urban areas, the weight of the homogeneity evaluation index is increased; in suburban farmland, the weight of the heterogeneity evaluation index is enhanced.

[0041] The non-dominated sorting genetic algorithm is used to perform multi-objective optimization on each selected segmentation scale, and the segmentation scale with the smallest comprehensive optimization objective function is selected as the optimal segmentation scale. The high-resolution remote sensing image is segmented according to the optimal segmentation scale to generate an initial segmentation result.

[0042] Calculate the JM distance between each image object and the adjacent objects in the initial segmentation result. When the JM distance is less than the target threshold T, merge the objects. Otherwise, retain the segmentation result. The target threshold T is determined as: Generate an initial object set based on the initial segmentation results, and calculate the JM distance between each initial image object and its adjacent objects to form a sample set. , the Gaussian kernel density estimation method is used to construct the sample set A probability distribution model is obtained by determining a peak point of the probability distribution model by maximum likelihood estimation; Peak point The candidate threshold interval is generated with T=0.95 as the center. Through grid search, it is found that the proportion of merged objects R(T)=68% when T=0.95 is in the confidence interval of natural merging rate (65%-70%), and the average JM distance is the largest (0.89), which meets the screening conditions. T=0.95 is selected as the target threshold T; When the JM distance is less than the target threshold T, the initial image object and the adjacent object are merged, otherwise the segmentation result is retained to obtain the final segmentation result.

[0043] A dynamic filling algorithm is used for the shadow areas of buildings, and images at different times of the year are compared to restore the underlying surface type of the blocked areas.

[0044] Based on the object-oriented data, multiple training subsets are generated by random sampling without replacement. Non-negative matrix factorization (NMF) is used to perform feature optimization and data dimensionality reduction on each training subset. The processed data is then input into the LightGBM base classifier for integrated classification to complete the extraction of urban underlying surface types: The original feature training set contains 12 dimensions (spectral, texture, NDVI, etc.), which are divided into training subsets, and the number of features in each training subset dimension, non-negative matrix factorization (NMF) extracts features and sets sparse coefficients ; The Brunet algorithm is used to perform the non-negative matrix decomposition process. The error of the non-negative matrix decomposition is measured by calculating the KL divergence. When the error change is less than the preset threshold or the maximum number of iterations is reached, the non-negative matrix decomposition is terminated and the basis matrix is ​​extracted. as the rotation factor; The rotation coefficients are reorganized in the order of the corresponding original training subsets and then multiplied with the feature training set to construct an enhanced feature set.

[0045] Number of LightGBM base classifiers and the number of feature subsets The LightGBM base classifier has a learning rate of 0.1 and a maximum depth of 8. By using feature importance feedback to improve the classification threshold, the classification accuracy is improved by 12%. The enhanced feature set is trained in parallel by L=50 LightGBM base classifiers, and the prediction results of L=50 LightGBM base classifiers are weighted voted to output the final underlying surface type of the main urban area (types include buildings, roads, vegetation, water areas, and bare soil). The vector data format is GeoJSON, and the overall accuracy is calculated using the confusion matrix. The Kappa coefficient is 0.89, the user accuracy of the building class reaches 94.1%, and the missed classification error of the road class is only 5.7%, which is better than the traditional random forest method.

[0046] The present invention can adaptively adjust the segmentation scale and feature weight according to the urban functional areas, fuse multi-temporal images to eliminate shadow interference, and improve the classification accuracy of shadow areas. The NMF feature dimensionality reduction shortens the LightGBM training time by 40%, making it suitable for large-scale image processing.

[0047] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may take the form of a fully hardware embodiment, a fully software embodiment, or an embodiment combining software and hardware. Furthermore, the present application may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code. The solutions in the embodiments of the present application may be implemented in various computer languages, such as the object-oriented programming language Java and the interpreted scripting language JavaScript.

[0048] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0049] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0050] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0051] Although the preferred embodiments of the present application have been described, those skilled in the art may make additional changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present application.

[0052] Obviously, those skilled in the art may make various changes and modifications to this application without departing from the spirit and scope of this application. Thus, if these modifications and variations of this application fall within the scope of the claims of this application and their equivalents, this application is intended to include these modifications and variations.

Claims

1. A method for extracting urban underlying surface types based on high-resolution remote sensing images, characterized in that: include: Determine the optimal segmentation scale of high-resolution remote sensing images, segment the high-resolution remote sensing images according to the optimal segmentation scale through multi-scale segmentation methods, and optimize the segmentation results with probabilistic distance to obtain object-oriented data; Based on the object-oriented data, multiple training subsets are generated by random sampling without replacement. Non-negative matrix factorization (NMF) is used to perform feature optimization and data dimensionality reduction on each training subset to obtain the processed training subset. The processed training subset is input into the LightGBM base classifier for image classification to complete the extraction of urban underlying surface types.

2. The method for extracting urban underlying surface types based on high-resolution remote sensing images according to claim 1, characterized in that: The method for obtaining object-oriented data specifically includes: A comprehensive optimization objective function is constructed by using a homogeneity evaluation index and a heterogeneity evaluation index, and the optimal segmentation scale of the high-resolution remote sensing image is determined according to the comprehensive optimization objective function by using the RMNE method; Segmenting the high-resolution remote sensing image according to the optimal segmentation scale using a multi-scale segmentation method to generate an initial segmentation result, which includes multiple initial image objects; Based on the initial segmentation result, the JM distance between two adjacent initial image objects is calculated. When the JM distance is less than the target threshold T, the two adjacent initial image objects are merged. Otherwise, the segmentation result is retained to obtain the final segmentation result. According to the final segmentation results, object-oriented data is constructed.

3. The method for extracting urban underlying surface types based on high-resolution remote sensing images according to claim 2, wherein: The RMNE method determines the optimal segmentation scale of the high-resolution remote sensing image according to the comprehensive optimization objective function, including: Generate a set of candidate segmentation scales for high-resolution remote sensing images , each candidate segmentation scale Corresponding to different spatial texture information entropy and spectral heterogeneity index , the spatial texture information entropy As an evaluation index of internal homogeneity of objects, the calculation formula is: （1） In formula (1), The first Probability distribution of texture-like features; The spectral heterogeneity index As an evaluation index of heterogeneity between adjacent objects, the calculation formula is: （2） In formula (2), For the current object Band mean, is the mean value of the corresponding band of adjacent objects, is the spectral standard deviation, is the weight coefficient and satisfies ; According to formula (1), the internal homogeneity evaluation index of each candidate segmentation scale is obtained, and according to formula (2), the heterogeneity evaluation index between adjacent objects corresponding to each candidate segmentation scale is obtained; Constructing a comprehensive optimization objective function based on the homogeneity evaluation index and the heterogeneity evaluation index , the comprehensive optimization objective function The calculation formula is: （3） In formula (3), is the maximum value of spatial texture information entropy in the candidate segmentation scale set, is the minimum value of spectral heterogeneity in the set of candidate segmentation scales, is the weight coefficient and satisfies ; The non-dominated sorting genetic algorithm is used to perform multi-objective optimization on each candidate segmentation scale to minimize the comprehensive optimization objective function. As the goal, generate the Pareto front solution set; According to the comprehensive optimization objective function The value of is selected from the Pareto front solution set to make the comprehensive optimization objective function Minimum candidate segmentation scale as the optimal segmentation scale.

4. The method for extracting urban underlying surface types based on high-resolution remote sensing images according to claim 2, wherein: The JM distance is calculated using the Gaussian kernel function for probability density estimation, and the calculation formula is: （4） In formula (4), is the Gaussian kernel bandwidth parameter, For two adjacent initial image objects in the The difference in feature dimensions; is the probability density function based on Gaussian distribution.

5. The method for extracting urban underlying surface types based on high-resolution remote sensing images according to claim 2, wherein: The method for determining the target threshold T includes: Based on the initial segmentation results, the initial object set is constructed and the JM distance between two adjacent initial image objects is calculated to form a sample set. ; The sample set is constructed using Gaussian kernel density estimation method The probability distribution model is used to determine the peak point of the probability distribution model by maximum likelihood estimation. ; Peak point As the center, generate candidate threshold intervals : （5） In formula (5), For the sample set The standard deviation of is the adaptive coefficient, The value of is based on the sample set Dynamic adjustment of skewness; Through grid search in the candidate threshold interval Step length Traversing candidates Calculate the value for each candidate The proportion of merged objects corresponding to the value; Select the candidate that meets the screening conditions from all candidate T values Value as target threshold , the screening conditions include: the proportion of merged objects of the target threshold T is within the statistical confidence interval of the natural merging rate of the sample set D; and the average JM distance between objects merged through the target threshold T is the largest; When there is no target threshold that meets the screening conditions in the candidate threshold range When , the sample set is reconstructed by Box-Cox transformation distribution after recalculation; when the recalculated candidate threshold interval still has no target threshold When as the target threshold T.

6. The method for extracting urban underlying surface types based on high-resolution remote sensing images according to claim 1, wherein: The method for extracting the urban underlying surface type specifically includes: Construct a feature training set based on object-oriented data, and set the number L of LightGBM base classifiers and the number K of feature subsets; randomly dividing the feature training set into K mutually exclusive training subsets according to the number K of the feature subsets, resampling each training subset, and applying a non-negative matrix factorization (NMF) method to the object-oriented data in the resampled training subset to generate rotation coefficients; The rotation coefficients are reorganized in the order corresponding to the corresponding training subsets and then multiplied with the feature training set to construct an enhanced feature set; The enhanced feature set is trained in parallel by L LightGBM base classifiers to obtain L trained LightGBM base classifiers; For the image objects to be classified, the L trained LightGBM base classifier models are used to perform predictions respectively to obtain L preliminary classification results; A weighted vote is performed on the L preliminary classification results based on the classification confidence, and the urban underlying surface type with the highest vote is output as the final extracted urban underlying surface type.

7. The method for extracting urban underlying surface types based on high-resolution remote sensing images according to claim 6, characterized in that: The division of the training subsets satisfies: （6） In formula (6), is the total number of features in the feature training set, is the number of preset feature subsets, To obtain a random perturbation term that obeys uniform distribution, K mutually exclusive training subsets are generated by Fisher-Yates shuffling algorithm, each training subset contains Features.

8. The method for extracting urban underlying surface types based on high-resolution remote sensing images according to claim 6, wherein: The object-oriented data in the training subset is subjected to a non-negative matrix factorization (NMF) method to generate rotation coefficients, including: Construct a non-negative data matrix for the object-oriented data in each training subset , the non-negative data matrix is ​​transformed into Decompose into basis matrix Coefficient matrix ; Define the objective function of non-negative matrix factorization (NMF): （7） In formula (7), is a non-negative data matrix, is the basis matrix, is the coefficient matrix, is the sparse regularization coefficient, is the number of features in the training subset, is the feature dimension of the training subset, is the preset implicit feature dimension; The objective function is minimized by alternating least squares method, and the basis matrix is ​​iteratively updated. and coefficient matrix ; The error of non-negative matrix decomposition after each iteration is measured by calculating KL divergence. When the error change is less than the preset threshold or the maximum number of iterations is reached, the decomposition of the non-negative data matrix is ​​terminated and the final converged basis matrix is ​​extracted. as the rotation factor.

9. An electronic device, characterized in that: include: A memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the method for extracting urban underlying surface types based on high-resolution remote sensing images according to any one of claims 1 to 8 is implemented.

10. A computer-readable storage medium storing a computer program, wherein the computer program enables a computer to execute the method for extracting urban underlying surface types based on high-resolution remote sensing images according to any one of claims 1 to 8.

Images