Heterogeneous landscape land use mixed ground object classification method based on homogeneous patch decomposition

Abstract

The application discloses a kind of based on the heterogeneous landscape land use mixed ground object classification method of homogeneous patch decomposition, comprising: reading original remote sensing image data and segmentation, generate multiscale patch set and construct patch hierarchical relationship diagram;Calculate patch heterogeneity index, extract gradient structure feature and realize bidirectional feature transmission;Build feature dictionary and based on sparse representation Quantification mixed degree of ground object;Determine the optimal scale by structure consistency index, and according to the heterogeneity degree executes adaptive classification strategy.The application solves the technical problems of optimal scale determination and mixed ground object representation of heterogeneous landscape, improves the classification precision of complex landscape land use, realizes fine land use classification.

Description

Technical Field

[0001] This invention belongs to the field of land use classification, and in particular, it is a method for classifying mixed land use features in heterogeneous landscapes based on homogeneous image decomposition. Background Technology

[0002] Land use classification of heterogeneous landscapes is a crucial step in remote sensing applications, playing a significant role in urban and rural planning, resource management, ecological monitoring, and environmental change research. With the improvement of high-resolution remote sensing data acquisition capabilities, traditional classification methods struggle to effectively handle landscape heterogeneity, particularly the complex structural characteristics of mixed land cover areas such as urban fringe areas and agroforestry transition zones. This leads to decreased classification accuracy and limits the potential of remote sensing technology in refined land use monitoring.

[0003] Currently, pixel-based statistical classification methods and object-based image analysis (OBIA) are two main paradigms for remote sensing image classification. The former labels individual pixels using clustering or supervised classification algorithms, ignoring spatial context information; the latter generates homogeneous patches through image segmentation and classifies them based on shape, texture, and contextual relationships, but its segmentation scale is singular, limiting its ability to handle complex landscapes with internal heterogeneity. In recent years, multi-scale analysis methods have been increasingly applied to remote sensing image processing, such as wavelet decomposition, scale-space filtering, and hierarchical segmentation. However, these methods often process information at different scales independently, lacking an effective mechanism for information transfer between scales. Furthermore, existing methods for extracting internal structural features of patches are mostly based on simple statistics, making it difficult to characterize gradient changes and directional features. In addition, the solution to the mixed pixel problem mainly relies on spectral mixing analysis, failing to fully utilize spatial structural information.

[0004] The core problems with existing technologies are twofold: firstly, the lack of an effective mechanism to determine the optimal analysis scale for heterogeneous landscapes, particularly the failure to quantitatively assess the consistency between the internal structure of a patch and the structure of its sub-patterns, leading to blind scale selection; secondly, the failure to establish a bidirectional feature transfer mechanism between different scales, preventing the effective fusion of global contextual information at higher scales and local detail information at lower scales, resulting in insufficient classification accuracy for highly heterogeneous regions. Furthermore, existing technologies fail to incorporate sparse representation theory to accurately quantify the composition of land cover in mixed land cover regions, making it difficult to implement adaptive classification strategies based on the degree of heterogeneity, thus limiting the application value of classification results in complex landscapes. Summary of the Invention

[0005] The purpose of this invention is to provide a method for classifying mixed land use features in heterogeneous landscapes based on homogeneous image decomposition, in order to solve at least one technical problem existing in the prior art.

[0006] The technical solution, a mixed land use classification method for heterogeneous landscapes based on homogeneous image patch decomposition, includes the following steps:

[0007] Read and segment the raw remote sensing image data to generate a multi-scale image patch set, construct the hierarchical relationship between the images, and obtain the image patch relationship map data;

[0008] Based on multi-scale speckle sets and speckle relationship diagram data, we calculate the heterogeneity index, gradient structure features, and sub-spot distribution features within speckles to achieve feature transfer and fusion between different scales and obtain a coupled feature set.

[0009] Read the multi-scale patch set and the original remote sensing image data, construct a feature dictionary and generate a sparse representation of the patches to obtain the sparse coefficient matrix;

[0010] Based on the image patch relationship map data, coupled feature set, sparse coefficient matrix and pre-stored training sample data, the optimal feature scale is determined and adaptive classification is performed. Multi-scale results are fused and classification results are optimized to generate a fine land use classification map.

[0011] Beneficial effects: This invention solves the technical problems of determining the optimal scale of heterogeneous landscapes and representing mixed land features, improves the accuracy of land use classification in complex landscapes, and realizes refined land use classification. Attached Figure Description

[0012] Figure 1 A flowchart illustrating the steps of a heterogeneous landscape land use mixed feature classification method based on homogeneous image patch decomposition, provided in this application embodiment.

[0013] Figure 2 A flowchart illustrating the steps for feature transfer and fusion between different scales provided in this application embodiment.

[0014] Figure 3 A flowchart illustrating the steps for constructing a feature dictionary and generating sparse representations of image spots, as provided in this application embodiment.

[0015] Figure 4 A flowchart illustrating the steps for fusing multi-scale results and optimizing classification results in an embodiment of this application. Detailed Implementation

[0016] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0017] It should be noted that, for the purpose of clearly demonstrating the steps of this application, each step has been numbered in the specification. These numbers are for ease of explanation only and do not limit the execution order of the steps. In actual operation, depending on the technical requirements of the specific implementation scenario, the steps may be executed in a different order than that shown in the specification, and in some cases, parallel processing between steps can also be achieved.

[0018] like Figure 1 As shown, the heterogeneous landscape land use mixed feature classification method based on homogeneous image patch decomposition includes the following steps:

[0019] S1. Read and segment the original remote sensing image data, generate a multi-scale image patch set, construct the hierarchical relationship between the image patches, and obtain the image patch relationship map data;

[0020] Specifically, raw remote sensing imagery data can capture detailed information about the Earth's surface, potentially including rich content such as topography, vegetation, and water bodies. This data is divided into many smaller patches, each corresponding to a specific area (such as a forest or a field), called "pixels," with the aim of analyzing local characteristics in greater detail. The different scales within a multi-scale patch set mean dividing patches by size, for example, from large "villages" to small "houses." Through such multi-scale analysis, we can understand the characteristics of ground features from macroscopic to microscopic levels.

[0021] S2. Based on multi-scale image patch set and image patch relationship map data, calculate the heterogeneity index, gradient structure features and sub-image patch distribution features within the image patch, realize feature transfer and fusion between different scales, and obtain a coupled feature set.

[0022] Specifically, heterogeneity indices refer to the "diversity" within a patch, such as the species diversity of plants in a forest area or the differences in crop distribution in farmland, aiming to understand the complexity within the patch. Gradient structure refers to certain trends of change, such as a gradual transition from flat terrain to hills, or from arid regions to wetlands; this feature helps identify transitional areas in the natural environment. Sub-patterns are smaller regions distributed within larger patches; analyzing these distribution characteristics can reveal details within large regions, such as the distribution of residential, commercial, and industrial areas in a town.

[0023] S3. Read the multi-scale patch set and the original remote sensing image data, construct the feature dictionary and generate the sparse representation of the patches to obtain the sparse coefficient matrix;

[0024] Specifically, the feature dictionary is like a "large dictionary" that records the characteristics of each image patch, such as color, texture, and shape. This allows for the classification and organization of various image patch properties for subsequent analysis. Sparse representation is a data compression method that emphasizes expressing the key features of an image patch with as little information as possible. For example, a complex image may contain thousands of details, but through sparse representation, a few key data points can be condensed, saving storage space and improving processing efficiency.

[0025] S4. Based on the image patch relationship map data, coupled feature set, sparse coefficient matrix and pre-stored training sample data, determine the optimal feature scale and perform adaptive classification, fuse multi-scale results and optimize classification results to generate a fine land use classification map.

[0026] Specifically, by analyzing the input data, a scale (size) best suited for the classification task is selected. This is because some features are more pronounced over a larger scale, while others may be more significant over a smaller scale. Using the selected scale and training sample data, a classification model is trained, assigning a label (e.g., city, forest, grassland) to each image patch. To improve classification accuracy, the classification results from different scales are integrated; and the results are further adjusted and improved to address potential classification errors, ensuring higher accuracy in the classification map.

[0027] This embodiment ensures that the classification results can reflect the macroscopic and microscopic characteristics of the region, avoiding the limitations of single-scale analysis; it improves the efficiency and computational performance of data processing, while retaining key features and reducing the impact of redundant information; it can dynamically adjust the classification strategy according to different patch characteristics, improving classification accuracy; and it fully explores the information potential of different data sources, thereby achieving a more comprehensive expression of ground feature characteristics.

[0028] According to one aspect of this application, the steps of constructing hierarchical relationships between image spots to obtain image spot relationship map data include:

[0029] S11. Read the raw remote sensing image data and preprocess it using Gaussian filtering to obtain the preprocessed image;

[0030] S12. Read the preprocessed image and use a multi-scale segmentation algorithm to generate image patches under different parameters, including scale, shape and compactness, to obtain a multi-scale image patch set, including {S1, S2, ..., Sm}, where S1 is the finest scale and Sm is the coarsest scale.

[0031] S13. Read the multi-scale image patch set, analyze the inclusion relationship between images at each scale, and construct the image patch hierarchy relationship map G = (V, E) by calculating the spatial overlap, where V is the image patch set at all scales and E is the inclusion relationship between images, thus obtaining the image patch relationship map data.

[0032] S14. Read the multi-scale image spot set and image spot relationship map data, calculate the spatial distribution characteristics of the image spots at each scale level, including area distribution, shape index and spatial clustering, and obtain the spatial characteristic data of the image spots.

[0033] like Figure 2 As shown, according to one aspect of this application, the steps for achieving feature transfer and fusion between different scales to obtain a coupled feature set include:

[0034] S21. Based on the multi-scale patch set and the original remote sensing image data, calculate the internal heterogeneity index of each patch to obtain the patch heterogeneity index matrix;

[0035] S22. Construct the gradient structure features for each image spot to obtain the gradient structure feature set;

[0036] S23. Based on the multi-scale image patch set and image patch relationship map data, analyze the distribution characteristics of child images in the parent image patch and the inclusion relationship between images of different scales defined in the image patch relationship map data, and obtain the child image patch distribution characteristic matrix.

[0037] S24. Based on the gradient structure feature set, the sub-pixel distribution feature matrix, and the speckle heterogeneity index matrix, feature transfer and fusion between different scales are realized to obtain the coupled feature set.

[0038] Specifically, the multi-scale patch set and the original remote sensing image data are read, and the internal heterogeneity index of each patch is calculated. The formula is: IHI(S) = w1·SD(S) + w2·Entropy(S) + w3·EdgeDensity(S) + w4·SpectralVariance(S); where SD is the standard deviation, Entropy is the entropy value, EdgeDensity is the edge density, SpectralVariance is the spectral variance, and w1 to w4 are weighting coefficients, thus obtaining the patch heterogeneity index matrix.

[0039] Read the multi-scale patch set and the original remote sensing image data, and construct the gradient structure features of each patch. The calculation steps are as follows: Apply the multi-directional gradient operator {▽1, ▽2, ..., ▽} to each patch S. k};Calculate the gradient structure features of the image patch: PGST(S) = ∑[▽ i S·▽ i S T ], T This indicates transpose; extract the principal gradient direction and intensity to obtain the gradient structure feature set.

[0040] Read the multi-scale patch set and patch relationship map data, and analyze the distribution characteristics of child patches in the parent patch: For each patch S, obtain the set of all its child patches C(S) = {c1, c2, ..., c n}; Extract sub-image patch distribution features: SPDF(S) = {f(dist(c1)), f(dist(c2)), ..., f(dist(c... n ))};where dist(c) describes the feature distribution of sub-image spot c, and f is the feature transformation function, which yields the feature matrix of sub-image spot distribution.

[0041] Read the gradient structure feature set, sub-pixel distribution feature matrix, and speckle heterogeneity index matrix to achieve feature transfer between different scales: calculate the scale transfer weight matrix W, where W(i,j) represents the information transfer weight from scale i to scale j; adopt a bidirectional feature transfer mechanism, including top-down and bottom-up transfer; update speckle features: F'(S) = F(S) + ∑[W(i,j)·F(S')], where S' is the related speckle of S, and F(S) is the original speckle feature; obtain the coupled feature set.

[0042] According to one aspect of this application, the steps of constructing gradient structure features for each image spot to obtain a gradient structure feature set include:

[0043] Extract the image sub-region corresponding to each image spot from the multi-scale image spot set;

[0044] A multi-directional gradient operator set is constructed and applied to each band to generate a multi-band gradient map. Based on this, the gradient intensity and direction of each image spot in each direction are calculated.

[0045] Gradient covariance matrix is ​​constructed based on gradient intensity and direction to extract image spot gradient structure features;

[0046] Based on the gradient structure characteristics of the image spot, the spatial distribution characteristics and structural change characteristics of the internal structure of the image spot are extracted to form a gradient structure feature vector.

[0047] By combining the image heterogeneity index matrix with the feature enhancement coefficients, a gradient structure feature set is generated.

[0048] Specifically, the multi-scale patch set and the original remote sensing image data are read, and a sub-region is extracted for each patch to ensure the integrity of the analysis boundary, thus obtaining the patch sub-region set.

[0049] Read the set of image subregions and construct a multi-directional gradient operator set {▽1, ▽2, ..., ▽ k Generate Sobel operators for the four principal directions: 0°, 45°, 90°, and 135°; for each direction j, compute the gradient operator ▽. j= [Gx_j, Gy_j], where Gx_j and Gy_j are the gradient components in the x and y directions, respectively; apply the gradient operator to each band to obtain a multi-band gradient set; obtain a multi-directional gradient operator set and a multi-band gradient map.

[0050] Read the multi-band gradient atlas and calculate the gradient intensity of each image spot S in each direction j: For each band b, calculate the gradient intensity G_bj(S) = sqrt(Gx_bj) 2 (S) + Gy_bj 2 Gx_bj(S) is the gradient component in the x-direction (horizontal direction); Gy_bj(S) is the gradient component in the y-direction (vertical direction); calculate the gradient direction θ_bj(S) = arctan(Gy_bj(S) / Gx_bj(S)); perform weighted fusion on all bands: G_j(S) = ∑wb·G_bj(S), where wb is the band weight; obtain the image spot gradient intensity matrix and the image spot gradient direction matrix.

[0051] Read the image spot gradient intensity matrix and image spot gradient direction matrix to construct the image spot gradient structure features: calculate the gradient covariance matrix M(S) = [∑G_j(S)·cos 2 θ_j(S), ∑G_j(S)·cosθ_j(S)·sinθ_j(S); ∑G_j(S)·cosθ_j(S)·sinθ_j(S), ∑G_j(S)·sin 2 θ_j(S)], where θ_j(S) is the overall gradient direction; perform eigenvalue decomposition on the covariance matrix M(S) to obtain eigenvalues ​​λ1≥λ2 and corresponding eigenvectors v1, v2; calculate gradient structure features: directionality D(S) = (λ1-λ2) / (λ1+λ2), intensity I(S) = λ1+λ2, coherence C(S) = λ1 / (λ1+λ2); obtain the basic set of image spot gradient structure features.

[0052] Read the basic set of image spot gradient structure features and extract the spatial distribution features of the internal structure of the image spot: divide the image spot into an n×n local grid {g1, g2, ..., g...} n 2}; Calculate the local gradient structure features {D(gi), I(gi), C(gi)} for each mesh gi; Extract the inter-mesh structure variation features: GV(S) = sqrt(∑(D(gi)-D*) 2 / n 2), GE(S) = -∑p(D(gi))·log(p(D(gi))), where GV(S) is the gradient change, GE(S) is the gradient entropy, D* is the average gradient directionality, and p(D(gi)) is the probability distribution of the gradient directionality D(gi) within grid gi; construct the complete gradient structure feature vector F(S) = [D(S), I(S), C(S), GV(S), GE(S), ...] to obtain the gradient structure feature set.

[0053] Read the gradient structure feature set and the speckle heterogeneity index matrix, and perform adaptive feature enhancement: Based on the speckle heterogeneity index IHI(S), calculate the adaptive enhancement coefficient β(S) = f(IHI(S)); enhance the structural features of high heterogeneity specks (IHI(S)>θ): F'(S) = F(S)·(1+β(S)), where θ is a preset threshold; keep the original features unchanged for low heterogeneity specks, and obtain the enhanced gradient structure feature set as the final gradient structure feature set output.

[0054] According to one aspect of this application, the steps for achieving feature transfer and fusion between different scales to obtain a coupled feature set include:

[0055] Based on multi-scale speckle sets and speckle relationship graph data, the spatial overlap and feature similarity between speckles of different scales are calculated, and a scale relationship weight matrix is ​​constructed.

[0056] Based on the scale relationship weight matrix, gradient structure feature set, sub-pixel distribution feature matrix, and speckle heterogeneity index matrix, top-down feature transfer is implemented, enabling parent speckle features to be transferred to child speckles, and the transfer intensity is adjusted according to the differences between parent and child speckle features to obtain the down-transferred feature set; bottom-up feature transfer is implemented, enabling sub-pixel information with significant features to be transferred to parent speckles to obtain the up-transferred feature set.

[0057] Read the downlink feature set and the uplink feature set, fuse the bidirectional transmitted features and normalize and reduce their dimensionality to obtain the coupled feature set; among them, downlink features are preferred for highly heterogeneous regions and uplink features are preferred for low heterogeneous regions.

[0058] According to one aspect of this application, the steps of constructing the scale relation weight matrix and implementing top-down feature transfer include:

[0059] When constructing the scale relationship weight matrix, the spatial overlap and feature similarity between images are taken into account. The weight calculation formula is: W(p,c) = OR(p,c)·FS(p,c);

[0060] Where OR(p,c) represents spatial overlap and FS(p,c) represents feature similarity;

[0061] When implementing top-down feature transfer, a differentiated transfer strategy is adopted. The transfer intensity is proportional to the difference in features between parent and child image patches. The transfer formula is: FD(ci) = F(ci) + γD·WD(p,ci)·[F(p) – F*(C(p))];

[0062] Where FD(ci) is the feature vector of sub-spot ci after top-down feature transfer; F(ci) is the original feature vector of sub-spot ci; γD is the downlink adjustment coefficient, which controls the intensity of top-down transfer; WD(p,ci) is the transfer weight from parent spot p to sub-spot ci; F(p) is the feature vector of parent spot p; F*(C(p)) is the average feature of all sub-spots of parent spot p, and C(p) represents the set of all sub-spots of p.

[0063] Specifically, the multi-scale patch set and patch relationship map data are read, and the scale relationship weights are calculated: a scale transition matrix T is constructed, where T(i,j) represents the transition probability from scale i to scale j; for each pair of parent and child patches (p,c), the spatial overlap OR(p,c) = Area(p∩c) / Area(c) is calculated; the feature similarity FS(p,c) = exp(-||F(p)-F(c)|| 2 / σ 2 The scale relationship weights W(p,c) = OR(p,c)·FS(p,c) are calculated to obtain the scale relationship weight matrix. Where σ 2 is the smoothing factor, and c is the sub-image spot.

[0064] By reading the scale relation weight matrix, gradient structure feature set, sub-pixel distribution feature matrix, and piece heterogeneity index matrix, a top-down feature transfer mechanism is constructed: for each parent piece p, collect the set of all its sub-pixels C(p) = {c1, c2, ..., c...} n}; Calculate the top-down weight WD(p, ci) = W(p, ci) / ∑W(p, cj); Implement feature downlink: FD(ci) = F(ci) + γD·WD(p, ci)·[F(p) –F*(C(p))], and obtain the downlink feature set.

[0065] Read the scale relationship weight matrix, gradient structure feature set, sub-pixel distribution feature matrix, and speckle heterogeneity index matrix to construct a bottom-up feature transfer mechanism: for each sub-pixel c, find its parent speckle p; calculate the bottom-up transfer weight WU(c, p) = W(c, p)·IHI(c) / ∑IHI(cj); realize feature upload: FU(p) = F(p) + γU·∑WU(ci, p)·[F(ci) - F(p)]; where γU is the upload adjustment coefficient, and W(c, p) is the scale relationship weight from sub-pixel c to parent speckle p; adaptively adjust according to the heterogeneity of the parent speckle to obtain the uploaded feature set.

[0066] Read the downlink feature set and the uplink feature set, and fuse the bidirectional transfer features: calculate the heterogeneity adaptive fusion coefficient α(S) = sigmoid(IHI(S)-θ); for each patch S, the fused feature F'(S) = α(S)·FD(S) + (1-α(S))·FU(S); for patches with high heterogeneity, more emphasis is placed on downlink features; for patches with low heterogeneity, more emphasis is placed on uplink features to obtain the fused feature set.

[0067] Read the fused feature set and perform feature normalization and optimization: standardize each feature dimension using Z-score; use principal component analysis (PCA) to reduce dimensionality and retain principal components that explain more than 95% of the variance; reconstruct the feature space to generate the final coupled features, and obtain the coupled feature set as the output.

[0068] like Figure 3 As shown, according to one aspect of this application, the steps of constructing a feature dictionary and generating a sparse representation of image patches to obtain a sparse coefficient matrix include:

[0069] S31. Based on the multi-scale patch set and the original remote sensing image data, extract the spectral, texture and shape features of each patch to form a patch feature vector set;

[0070] S32. Based on the image patch feature vector set, create a land cover feature dictionary that can represent the feature patterns of different land cover types;

[0071] S33. Use the feature dictionary of land features to perform sparse decomposition on the feature vector set of the image patch. By optimizing the solution of the sparse coefficients, each image patch can be represented as a linear combination of dictionary elements, generating a sparse coefficient matrix that reflects the composition ratio of land features.

[0072] S34. Calculate the degree of land cover mixing in each image patch based on the sparse coefficient matrix, quantify the mixing status of land cover types in heterogeneous landscapes, and form a land cover mixing index.

[0073] Specifically, the multi-scale patch set and the original remote sensing image data are read, and the standard feature vector of each patch is extracted, including spectral features, texture features and shape features, to obtain the patch feature vector set.

[0074] Read the feature vector set of image patches and construct a dictionary of land cover types D = {d1, d2, ..., d...} n}, where each dictionary element represents a feature pattern of a land cover type, thus obtaining a land cover feature dictionary.

[0075] Read the feature vector set of the image patches and the dictionary of land cover features, and perform a sparse representation of each image patch S: represent the image patch as a sparse linear combination of the dictionary: S≈D·α; use the Least Absolute Value Shrinking and Selection (LASSO) algorithm to solve for the sparsity coefficient α, with the optimization objective being: min ||S - D·α||2 2 + λ||α||1; λ is the regularization parameter in the sparse representation optimization process; the sparse coefficient α reflects the proportion of different land cover types within the image patch, resulting in the sparse coefficient matrix.

[0076] Read the sparse coefficient matrix, analyze the distribution characteristics of the sparse coefficients, and calculate the ground cover mixing index in the image patch. The formula is: HMI(S) = 1 - max(α) / sum(α); where HMI is the mixing index, and the larger the value, the higher the mixing degree, thus obtaining the ground cover mixing index.

[0077] like Figure 4 As shown, according to one aspect of this application, the steps of fusing multi-scale results and optimizing classification results to generate a refined land use classification map include:

[0078] S41. Read the image patch relationship map data and image patch heterogeneity index matrix, calculate the consistency between the internal structure of the image patch and the structure of the sub-image patch, generate the structure consistency index matrix, and determine the optimal analysis scale for each region to form the optimal scale mapping map.

[0079] S42. Combining the optimal scale mapping map, coupled feature set, sparse coefficient matrix and training sample data, implement a differentiated classification strategy based on the degree of regional heterogeneity to obtain preliminary classification results;

[0080] S43. Optimize the spatial context relationship based on the preliminary classification results to generate an optimized classification map;

[0081] S44. By comprehensively optimizing the classification map and the sparse coefficient matrix, fine-grained classification adjustments are made to obtain a refined land use classification map.

[0082] Specifically, the image patch relationship map data, image patch heterogeneity index matrix, and gradient structure feature set are read, and the structural consistency index of each image patch is calculated: the consistency between the internal structure of the image patch and the structure of its sub-image patches is evaluated: SCI(S) = sim(PGST(S), {PGST(c)| c ∈ C(S)}); based on the structural consistency index matrix SCI and the heterogeneity index, the optimal feature scale λ(S) of each region is determined, and the optimal scale mapping map is obtained.

[0083] Read the coupled feature set, sparse coefficient matrix, optimal scale map, and training sample data, and perform multi-scale classification: apply different classification strategies for patches with different degrees of heterogeneity; use random forest classification based on coupled features for patches with high homogeneity; perform hybrid classification by combining sparse coefficients for patches with high heterogeneity; select the best classification result for each region based on the optimal scale map to obtain the preliminary classification result.

[0084] Read the preliminary classification results, the image heterogeneity index matrix, and the gradient structure feature set, and optimize the classification results using a conditional random field model: Construct the conditional random field energy function: E(X) = ∑ψ u (x i ) + ∑ψ p (x i x j ); where x i and x j These represent the category labels for pixels or image spots, respectively. u For the potential energy of a single pixel, ψ p The energy is the potential energy of paired pixels; the classification result is optimized by minimizing the energy function to obtain the optimized classification map.

[0085] Read the optimized classification map and sparse coefficient matrix, and calculate the classification confidence of each pixel: assess the classification reliability based on the consistency of classification probability and sparse coefficient; for mixed patches, generate a sub-pixel level classification ratio map; integrate the classification results and confidence assessment to obtain the final refined land use classification map and classification confidence map.

[0086] According to one aspect of this application, the steps of determining the optimal analysis scale for each region and forming an optimal scale mapping map include:

[0087] Analyze the image patch relationship map data, identify the sub-image patch sets contained in each image patch and their spatial distribution characteristics, and generate a simplified image patch relationship map;

[0088] Based on the simplified image patch relationship diagram, the structural feature similarity between the image patch and its sub-image patches is calculated to form a structural similarity matrix; the quantity diversity, spatial distribution and type differences of the sub-image patches are analyzed to assess the compositional complexity within the image patch and establish a compositional complexity matrix; the structural similarity matrix and the compositional complexity matrix are combined to construct a structural consistency index matrix.

[0089] Analyze the structural consistency index matrix and the patch heterogeneity index matrix to determine the scale corresponding to the maximum structural consistency for each region and generate the optimal scale mapping map.

[0090] According to one aspect of this application, the steps for obtaining the composition of the complexity matrix and the construction of the structural consistency index matrix include:

[0091] By combining the sub-image patch quantity diversity index, spatial distribution entropy index, and type diversity index, a comprehensive evaluation is formed through weighted fusion, resulting in a composition complexity matrix.

[0092] The formula for calculating the structural consistency index matrix is: SCI(S) = ω(S)·SSavg(S) - (1-ω(S))·CC(S);

[0093] Where ω(S) is the adaptive weight function, SSavg(S) is the average structural similarity matrix, and CC(S) is the compositional complexity.

[0094] When determining the optimal scale, the optimal combination of structural consistency index and heterogeneity threshold among all patches covering that location is found to ensure that the selected scale maintains structural consistency without overcomplicating the structure.

[0095] Specifically, the image patch relationship map data, gradient structure feature set, and image patch heterogeneity index matrix are read, and the image patch relationship data is preprocessed: the set of all sub-image patches C(S) = {c1, c2, ..., c...} for each image patch S is identified. n}; Calculate the spatial distribution statistical characteristics of sub-image spots: center distance, directional distribution, size ratio, etc.; screen important sub-image spots, remove noise and small sub-image spots, and obtain a simplified image spot relationship diagram.

[0096] Read the simplified image patch relationship diagram and gradient structure feature set, and calculate the similarity between the internal structure of the image patch and the structure of the sub-image patches: extract the gradient structure feature vector F_g(S) = [D(S), I(S), C(S), ...] of image patch S; extract the gradient structure feature vector F_g(ci) of each sub-image patch ci; calculate the structure similarity SS(S, ci) = cos(F_g(S), F_g(ci)); calculate the weighted average structure similarity: SSavg(S) = ∑[Area(ci) / Area(S)]·SS(S, ci), and obtain the structure similarity matrix.

[0097] Read the simplified patch relationship diagram and sub-patch distribution feature matrix, and calculate the compositional complexity of the patch: Calculate the quantity diversity of sub-patches within patch S: NDC(S) = n / log(Area(S)); Calculate the entropy of the spatial distribution of sub-patches: SPDE(S) = -∑p(ci)·log(p(ci)); Calculate the type diversity of sub-patches: SPDT(S) = 1-∑[Area(ci) / Area(S)] 2 The overall compositional complexity is calculated as: CC(S) = w1·NDC(S) + w2·SPDE(S) + w3·SPDT(S), yielding the compositional complexity matrix. Here, n represents the number of sub-pixels within plaque S, Area(S) is the area of ​​plaque S, and p(ci) is the proportion of the area occupied by sub-pixels ci within the entire plaque S.

[0098] Read the structural similarity matrix, composition complexity matrix, and speckle heterogeneity index matrix, and calculate the structural consistency index: Design an adaptive weight function: ω(S) = exp(-IHI(S) / σ), where σ is a smoothing parameter; Calculate the structural consistency index: SCI(S) = ω(S)·SSavg(S) - (1-ω(S))·CC(S); The higher the SCI value, the more consistent the internal structure of the speckle; the lower the SCI value, the less consistent the internal structure, thus obtaining the structural consistency index matrix.

[0099] Read the structural consistency index matrix, the speckle heterogeneity index matrix, and the multi-scale speckle set to determine the optimal feature scale for each region: for each location (x, y), find the maximum structural consistency index across all scales: SCI*(x, y) = max{SCI(S)|S contains (x, y)}; simultaneously consider the heterogeneity threshold constraint: S is considered a valid scale only when IHI(S) < θ; the optimal scale λ*(x, y) = arg max{SCI(S)·(θ-IHI(S)) / θ|S contains (x, y)}; perform spatial smoothing on the optimal scale map to eliminate scale breaks and obtain the optimal scale mapping map.

[0100] According to one aspect of this application, the steps for implementing a differentiated classification strategy based on the degree of regional heterogeneity to obtain preliminary classification results include:

[0101] By reading the coupled feature set, sparse coefficient matrix, optimal scale mapping map and training sample data, the image patch regions are divided into three categories: low heterogeneity, medium heterogeneity and high heterogeneity.

[0102] Constructing a hierarchical classifier: For low heterogeneity regions, a random forest classifier based on coupling features is used; for medium heterogeneity regions, a support vector machine combining coupling features and sparse coefficients is used; and for high heterogeneity regions, a Gaussian mixture model is built using sparse coefficients and gradient structure features, resulting in a hierarchical classification model set.

[0103] Based on the optimal scale mapping map and hierarchical classification model set, the optimal scale patch to which each pixel belongs is determined, and the corresponding classifier is selected according to its heterogeneity. For highly heterogeneous regions, pixel classification label map and classification probability map are generated.

[0104] By combining the pixel classification label map, classification probability map, and sparse coefficient matrix, the proportion of mixed land features is calculated, the classification probability of highly heterogeneous areas is adjusted, classification fragmentation is eliminated, and preliminary classification results are formed.

[0105] Specifically, the coupled feature set, sparse coefficient matrix, optimal scale mapping map, and training sample data are read and partitioned according to the degree of heterogeneity: the images are divided into three categories: low heterogeneity (IHI < θ1), medium heterogeneity (θ1 ≤ IHI < θ2), and high heterogeneity (IHI ≥ θ2); a feature subset is constructed for each category, and the most relevant features are extracted; the feature weights of each category are optimized: w(f) = MI(f, Y) / H(f), where MI(f, Y) is the mutual information between feature f and category label Y, and H is the entropy, thus obtaining the partitioned feature set and feature weight matrix.

[0106] Read the partition feature set, feature weight matrix, and training sample data to construct a hierarchical classifier: for low heterogeneity regions, use the Random Forest classifier RF1 with coupled features as the feature set; for medium heterogeneity regions, use the Support Vector Machine classifier SVM with coupled features plus sparse coefficients as the feature set; for high heterogeneity regions, use the Gaussian Mixture Model GMM with sparse coefficients plus gradient structure features as the feature set; train each classifier to obtain the classification model parameters and classification accuracy, thus obtaining a hierarchical classification model set.

[0107] Read the hierarchical classification model set, partition feature set, and optimal scale mapping map, and perform multi-scale adaptive classification: For each pixel p, determine its optimal scale patch S*(p); select the corresponding classifier based on the heterogeneity index of the optimal scale patch S*(p); classify pixel p using the selected classifier to obtain the class label y(p) and classification probability P(y|p); for highly heterogeneous regions, generate a multi-class probability distribution {P(y1|p), P(y2|p), ..., P(y...}}. k |p)}, thus obtaining the pixel classification label map and classification probability map.

[0108] Read the pixel classification label map, classification probability map, sparse coefficient matrix, and gradient structure feature set to optimize the classification results: For each highly heterogeneous patch S, calculate the mixed ground cover ratio: MP(S, yi) = ∑α(S)i / ∑α(S)j; Adjust the classification probability using the sparse coefficient α: P'(yi|p) = P(yi|p)·MP(S, yi) / ∑P(yj|p)·MP(S, yj); Apply Markov random field smoothing, with the energy function: E(y) = ∑[1-P'(yi|p)] + β·∑[y(p)≠y(q)]; where β is the smoothing adjustment parameter, y(p) is the class label of pixel p, and y(q) is the class label of pixel q; Use the graph cut algorithm to optimize the label configuration and obtain the preliminary classification results.

[0109] Read the preliminary classification results and classification probability map, and refine the results: detect uncertain regions (P'(yi|p)<τ) in the classification results, where τ is the threshold for uncertain regions; apply a region growing algorithm to uncertain regions and fuse information from adjacent high-confidence regions; filter out small isolated regions to maintain the spatial continuity of the classification results; generate multi-level classification results, including pixel-level and object-level results, and output the corrected preliminary classification results.

[0110] In a specific embodiment of this application, Xingyang City, Henan Province, is taken as the study area. Based on the 2017 Gaofen-1 (GF-1) remote sensing image data, and considering the complex and diverse land use types in this region, the heterogeneous landscape land use mixed feature classification method based on homogeneous patch decomposition proposed in this application is applied for research practice. The geographical coordinates of the study area are 34°36′—34°59′ N and 113°7′—113°30′ E, with a total area of ​​943 square kilometers. The terrain is mainly plains and hills, and the land use types include nine types: cultivated land, forest land, grassland, orchard, urban construction land, other construction land, transportation land, water bodies, and other forest land, exhibiting typical complex and heterogeneous landscape characteristics.

[0111] Step 1: Data Acquisition and Preprocessing.

[0112] 1.1 Research data.

[0113] The main data used in this embodiment include: Gaofen-1 (GF-1) satellite remote sensing imagery: acquired in August 2017, containing panchromatic bands (spatial resolution 2 meters) and multispectral bands (red, green, blue, and near-infrared bands, spatial resolution 8 meters); Xingyang City 2016 land use vector data: used to provide prior knowledge for classification and accuracy verification benchmarks.

[0114] 1.2 Data preprocessing.

[0115] The GF-1 remote sensing image was preprocessed as follows: Radiometric correction: The original (DN) values ​​in the remote sensing image were converted into reflectance to eliminate the effects of atmospheric scattering and absorption; Geometric correction: The image was geometrically fine-corrected based on control points to ensure spatial positioning accuracy better than 1 pixel; Panchromatic and multispectral fusion: A Gram-Schmidt fusion method based on high-pass filtering was used to generate a 2-meter resolution multispectral image; Study area cropping: The fused image was cropped based on administrative boundary vector data to obtain the remote sensing image of the study area.

[0116] Step 2: Generate multi-scale layered image spots.

[0117] 2.1 Perform multi-scale segmentation.

[0118] A multi-scale segmentation algorithm was employed, with different scale parameters to generate multi-level image patches: a scale parameter sequence {10, 20, 40, 80, 160} was set, generating image patch sets {S1, S2, S3, S4, S5} at five scale levels, where S1 is the finest scale (scale parameter 10) and S5 is the coarsest scale (scale parameter 160). A shape factor of 0.3 and a compactness factor of 0.5 were set during segmentation to ensure that shape features are considered while taking spectral information into account, making it suitable for heterogeneous landscape analysis. Weighting coefficients were set for different bands: 0.4 for near-infrared band, 0.3 for red band, 0.2 for green band, and 0.1 for blue band, enhancing the separation between vegetated and non-vegetated areas. The results of the multi-scale segmentation showed that scale S1 generated 28,645 image patches with an average area of ​​approximately 700 square meters; scale S5 generated 1,253 image patches with an average area of ​​approximately 15,000 square meters.

[0119] 2.2 Constructing the hierarchical relationship of image spots.

[0120] Based on a multi-scale patch set, a patch hierarchy graph G = (V, E) is constructed: the spatial overlap OR(p, c) between any two patches p and c of different scales is calculated: OR(p, c) = Area(p∩c) / Area(c), where Area(p∩c) represents the spatial overlap area of ​​patches p and c, and Area(c) represents the area of ​​patch c. When OR(p, c) > 0.9 and the scale parameter of p is greater than the scale parameter of c, a hierarchical inclusion relationship from p to c is established, i.e., a directed edge (p, c) is added to the hierarchy graph G. The patch hierarchy graph G is constructed by traversing all patch pairs. The constructed graph contains 52,487 nodes (representing patches of all scales) and 123,654 edges (representing inclusion relationships between patches).

[0121] Compared with traditional single-scale segmentation methods, the multi-scale image patch hierarchy map constructed in this embodiment can simultaneously preserve the macroscopic structure and microscopic details of the landscape, providing multi-level information support for subsequent heterogeneous landscape analysis and effectively solving the problem that a single scale is difficult to adapt to complex landscape features.

[0122] Step 3: Conduct a comprehensive analysis of image spot characteristics.

[0123] 3.1 Calculate the image patch heterogeneity index.

[0124] For each patch S, the internal heterogeneity index is calculated using the formula: IHI(S) = w1·SD(S) + w2·Entropy(S) + w3·EdgeDensity(S) + w4·SpectralVariance(S); where SD(S) is the standard deviation of pixels within patch S, normalized to the [0, 1] interval; Entropy(S) is the entropy value of pixels within patch S, measuring the degree of disorder in the spectral distribution; EdgeDensity(S) is the edge density within patch S, calculated after extracting the edges using the Canny operator; SpectralVariance(S) is the spectral variance within patch S, measuring spectral heterogeneity; w1 to w4 are weighting coefficients, set to 0.3, 0.25, 0.2, and 0.25, respectively.

[0125] In this embodiment, the heterogeneity index was normalized using the maximum entropy method to obtain a patch heterogeneity index matrix ranging from [0, 1]. The results show that the patch heterogeneity index of the urban-rural transition zone and the agroforestry ecotone in the study area is generally higher than 0.7, while the heterogeneity index of the pure farmland area and the urban built-up area is mostly lower than 0.3.

[0126] 3.2 Extracting gradient structure features.

[0127] Compared to traditional methods that rely solely on statistical features, multi-directional gradient analysis can capture the structural and directional characteristics within image patches, making it particularly suitable for heterogeneous landscapes with distinct textures and structures. The specific steps are as follows:

[0128] Constructing multi-directional gradient operators: Generate Sobel operators for four principal directions: 0°, 45°, 90°, and 135°; define the gradient operator ▽ for each direction j. j = [Gx_j, Gy_j]; Apply these gradient operators individually to each band b of the multispectral image to obtain a multiband gradient atlas.

[0129] Calculate gradient strength and direction: For each band b, calculate the gradient strength G_bj(S) = sqrt(Gx_bj). 2 (S) +Gy_bj 2(S)); Calculate the gradient direction θ_bj(S) = arctan(Gy_bj(S) / Gx_bj(S)); Calculate the gradient direction based on the importance of the bands: G_j(S) = ∑wb·G_bj(S), where wb is the band weight, the weight of the near-infrared band is 0.4, the weight of the red band is 0.3, the weight of the green band is 0.2, and the weight of the blue band is 0.1.

[0130] Constructing gradient structure features: Calculate the gradient covariance matrix M(S) of the image patch S: M(S) = [∑G_j(S)·cos 2 θ_j(S), ∑G_j(S)·cosθ_j(S)·sinθ_j(S); ∑G_j(S)·cosθ_j(S)·sinθ_j(S), ∑G_j(S)·sin 2 [θ_j(S)]; Perform eigenvalue decomposition on the covariance matrix M(S) to obtain eigenvalues ​​λ1≥λ2 and corresponding eigenvectors v1 and v2; Calculate three basic gradient structure indices: Directionality D(S) = (λ1-λ2) / (λ1+λ2), with a value range of [0, 1], where a larger value indicates stronger structural directionality; Intensity I(S) = λ1+λ2, representing the overall intensity of gradient changes; and Continuity C(S) = λ1 / (λ1+λ2), representing the continuity of the gradient structure.

[0131] Extracting spatial distribution features: Divide the image spot S into a 4×4 local grid {g1, g2, ..., g...} 16}; Calculate the local gradient structure features {D(gi), I(gi), C(gi)} for each grid; Extract the structural variation features between grids: structural variability GV(S) = sqrt(∑(D(gi)-D*) 2 / 16); Structural entropy GE(S) = -∑p(D(gi))·log(p(D(gi))); Construct the complete gradient structure feature vector F(S) = [D(S), I(S), C(S), GV(S), GE(S), ...].

[0132] Adaptive feature enhancement: Based on the speckle heterogeneity index IHI(S), the adaptive enhancement coefficient β(S) = 2 / (1+e (-3·(IHI(S)-0.5)) For high heterogeneity patches (IHI(S) > 0.6), enhance structural features: F'(S) = F(S)·(1+β(S)); for low heterogeneity patches, maintain the original features.

[0133] Analysis of the study area revealed that the directional index D(S) in agroforestry ecotones typically ranged from 0.6 to 0.8, with a high structural variability GV(S). Urban areas, on the other hand, exhibited a larger intensity index I(S) but weaker directionality. Cultivated land areas, however, demonstrated strong directional consistency and high continuity. This embodiment effectively distinguished different types of mixed land cover areas, improving classification accuracy.

[0134] 3.3 Analyze the distribution characteristics of sub-image spots.

[0135] Based on the image patch relationship diagram, the distribution characteristics of child images in the parent image patch are analyzed: For each image patch S, the set of all its child images C(S) = {c1, c2, ..., c3} is obtained. n}; Calculate the spatial distribution characteristics of sub-spots: Sub-spot quantity diversity: NDC(S) = n / log(Area(S)); Sub-spot area distribution entropy: AAE(S) = -∑[Area(ci) / Area(S)]·log[Area(ci) / Area(S)]; Sub-spot distance variance: DV(S) = Var(dist(ci,cj)), where dist is the Euclidean distance; Sub-spot shape heterogeneity: SH(S) = Var(Shape(ci)), where Shape is the shape index; Construct the sub-spot distribution feature matrix: SPDF(S) = [NDC(S), AAE(S), DV(S), SH(S)].

[0136] 3.4 Implement the feature transfer mechanism.

[0137] By designing a bidirectional feature transfer mechanism, information exchange between different scales is achieved, overcoming the limitation of independent information between scales in traditional methods. The specific implementation steps are as follows:

[0138] Calculate the scale relationship weights: For each pair of parent-child image patches (p, c), calculate the spatial overlap OR(p, c) = Area(p∩c) / Area(c); calculate the feature similarity FS(p, c) = exp(-||F(p)-F(c)|| 2 / σ 2 ), where σ 2 =2.0, F is the gradient structure feature vector; the comprehensive calculation of the scale relationship weight W(p,c) = OR(p,c)·FS(p,c).

[0139] Perform top-down feature transfer: For each parent image patch p, calculate the top-down transfer weight WD(p, ci) = W(p, ci) / ∑W(p, cj); implement the feature transfer formula: FD(ci) = F(ci) + γD·WD(p, ci)·[F(p) – F*(C(p))]; where F*(C(p)) is the average feature of all sub-image patches, and γD is the transfer adjustment coefficient, set to 0.6.

[0140] Perform bottom-up feature transfer: Calculate the bottom-up transfer weight WU(c,p) = W(c,p)·IHI(c) / ∑IHI(cj); Implement the feature upload formula: FU(p) = F(p) + γU·∑WU(ci,p)·[F(ci) - F(p)], where γU is the upload adjustment coefficient, set to 0.4.

[0141] Fusion of bidirectional transport characteristics: Calculation of heterogeneous adaptive fusion coefficient α(S) = 1 / (1+e (-(IHI(S)-0.5) / 0.1) The fusion feature formula is: F'(S) = α(S)·FD(S) + (1-α(S))·FU(S). For pixels with high heterogeneity (IHI > 0.6), the α value approaches 1, and more emphasis is placed on downlink features; for pixels with low heterogeneity (IHI < 0.4), the α value approaches 0, and more emphasis is placed on uplink features.

[0142] Feature normalization and optimization are performed: Z-score standardization is applied to each feature dimension; Principal component analysis (PCA) is used for dimensionality reduction, retaining principal components that explain more than 95% of the variance (in this embodiment, the first 8 principal components are retained); the feature space is reconstructed to generate the final coupled feature set.

[0143] By implementing a bidirectional feature transfer mechanism, information at different scales was effectively fused. Results show that for complex terrain areas in Xingyang City, the classification accuracy improved by an average of 8.7% after feature transfer. Compared to traditional unidirectional or independent scale analysis methods, the bidirectional feature transfer mechanism in this embodiment can simultaneously consider global structure and local details, making it particularly suitable for processing mixed landform areas in heterogeneous landscapes and solving the technical challenge of effectively fusing information at different scales.

[0144] Step 4: Sparse representation based on dictionary learning.

[0145] 4.1 Extracting the feature vector of the image patch.

[0146] Standard feature vectors are extracted from the multi-scale image patch set: spectral features: calculate the mean, standard deviation, maximum, minimum and median of each band, a total of 20 features; texture features: calculate homogeneity, contrast, entropy, correlation and other 8 features based on the gray-level co-occurrence matrix; shape features: calculate area, perimeter, compactness, major axis and minor axis ratio and other 6 features; context features: calculate the boundary length ratio with adjacent images and the average spectral difference and other 4 features; finally, a 38-dimensional feature vector is constructed for each image patch to form an image patch feature vector set.

[0147] 4.2 Construct a dictionary of geographic features.

[0148] Based on the feature vector set of image spots, a land cover type dictionary is constructed: representative image spots are selected for each land cover type from the training samples, with a total of 270 sample image spots selected; the dictionary elements are iteratively optimized using the K-SVD algorithm, with the initial dictionary size set to 45 (approximately 1 / 6 of the number of training samples); the stopping condition is when the reconstruction error is lower than the threshold (set to 0.01) or the maximum number of iterations is reached (set to 100); the final land cover feature dictionary D contains 45 dictionary elements, each of which is a 38-dimensional feature vector.

[0149] 4.3 Solve for the sparsity coefficient.

[0150] For each patch S, a sparse representation is performed based on the land cover feature dictionary: the patch is represented as a sparse linear combination of the dictionary: S ≈ D·α; the sparse coefficient α is solved using the LASSO algorithm, and the optimization objective is: min ||S - D·α||2 2 + λ||α||1, where λ is the regularization parameter, set to 0.1. The solution is obtained iteratively using the coordinate descent method, with a maximum number of iterations of 500 and a convergence threshold of 1e-6; this yields a sparse coefficient matrix reflecting the proportion of different land cover types within each image patch.

[0151] 4.4 Analyze the degree of land cover mixing.

[0152] Based on the sparse coefficient matrix, the degree of vegetation mixing within image patches is analyzed: the vegetation mixing degree index is calculated as: HMI(S) = 1 - max(α) / sum(α); where HMI is the mixing index, and a larger value indicates a higher degree of mixing. The analysis results show that: the HMI value of images in the urban-rural transition zone is generally higher than 0.7, indicating high mixing; the HMI value of images in pure farmland areas is mostly lower than 0.3, indicating relatively pure farmland; and the HMI value of images in agroforestry ecotones is mostly between 0.4 and 0.6, indicating a moderate degree of mixing.

[0153] Compared to traditional spectral mixing analysis, this embodiment employs a dictionary-based sparse representation method that can more accurately quantify the degree of land cover mixing within image patches, providing a new technical perspective for heterogeneous landscape classification. By analyzing the sparsity coefficients, the proportion of different land cover components within each image patch can be quantitatively assessed, effectively solving the problem of difficulty in quantifying and representing mixed land cover areas in traditional methods.

[0154] Step 5: Determine the optimal scale and perform multi-scale fusion classification.

[0155] 5.1 Calculate the structural consistency index.

[0156] The structural consistency assessment addresses the problem of blind scale selection in traditional methods by performing speckle relationship preprocessing: identifying all sub-species sets C(S) = {c1, c2, ..., c...} for each speckle S. n}; Filter out sub-image spots with spatial overlap greater than 0.1, remove noise and small sub-image spots, and obtain a simplified image spot relationship map.

[0157] Calculate structural similarity. Extract the gradient structural feature vector F_g(S) = [D(S), I(S), C(S), GV(S), GE(S)] of the image spot S; extract the corresponding feature vector F_g(ci) for each sub-image spot ci; calculate the structural similarity SS(S, ci) = cos(F_g(S), F_g(ci)); calculate the weighted average structural similarity: SSavg(S) = ∑[Area(ci) / Area(S)]·SS(S, ci).

[0158] Calculate the compositional complexity. Calculate the numerical diversity of sub-patterns within patch S: NDC(S) = n / log(Area(S)); calculate the entropy of the spatial distribution of sub-patterns: SPDE(S) = -∑p(ci)·log(p(ci)); calculate the type diversity of sub-patterns: SPDT(S) = 1 - ∑[Area(ci) / Area(S)] 2 The overall complexity is calculated as follows: CC(S) = 0.3·NDC(S) +0.4·SPDE(S) + 0.3·SPDT(S).

[0159] Calculate the structural consistency index. Design an adaptive weighting function: ω(S) = exp(-IHI(S) / 0.3); Calculate the structural consistency index: SCI(S) = ω(S)·SSavg(S) - (1-ω(S))·CC(S); The higher the SCI value, the more consistent the internal structure of the image patch; the lower the SCI value, the less consistent the internal structure.

[0160] Determine the optimal scale: For each location (x, y), find the maximum structural consistency index among all patches containing that location: SCI*(x, y) = max{SCI(S)|S contains (x, y)}; Simultaneously consider the heterogeneity threshold constraint: S is considered an effective scale only when IHI(S) < 0.7; Determine the optimal scale: λ*(x, y) = arg max{SCI(S)·(0.7-IHI(S)) / 0.7|S contains (x, y)}; Perform a 3×3 kernel median filter on the optimal scale map to eliminate scale breaks.

[0161] The optimal scale analysis results show that in the study area of ​​Xingyang City, the optimal analysis scale for pure farmland areas is mostly S4 (scale parameter 80), for urban built-up areas it is mostly S3 (scale parameter 40), while for urban-rural transition zones and agroforestry ecotones it is mostly S1 or S2 (scale parameter 10 or 20), which is highly consistent with the spatial variation characteristics of actual land features.

[0162] Compared with traditional fixed-scale or experience-based selection methods, the structural consistency index proposed in this embodiment can automatically determine the optimal analysis scale for each region, effectively solving the problem of blind scale selection and improving the accuracy of heterogeneous landscape classification.

[0163] 5.2 Multi-scale classification strategy.

[0164] Implement an adaptive classification strategy based on the optimal scale mapping:

[0165] Heterogeneity partitioning: Low heterogeneity region: IHI < 0.4; Medium heterogeneity region: 0.4 ≤ IHI < 0.7; High heterogeneity region: IHI ≥ 0.7.

[0166] Constructing hierarchical classifiers: Low heterogeneity regions: Use a random forest classifier, selecting coupled features as input, with the number of trees set to 200; Medium heterogeneity regions: Use a support vector machine (SVM) classifier, selecting coupled features + sparse coefficients as input, with the kernel function being RBF, regularization parameter C=100, and kernel function parameter γ=0.01; High heterogeneity regions: Use a Gaussian mixture model (GMM), selecting sparse coefficients + gradient structure features as input, with the number of Gaussian components being 3 times that of the land use category.

[0167] Perform multi-scale adaptive classification: For each pixel p, determine its optimal scale patch S*(p); select the corresponding classifier based on the heterogeneity index IHI of S*(p); classify pixel p using the selected classifier to obtain the class label y(p) and classification probability P(y|p); for highly heterogeneous regions, generate a multi-class probability distribution {P(y1|p), P(y2|p), ..., P(y... n |p)}.

[0168] Optimize classification results: For highly heterogeneous patches S, adjust the classification probability using a sparse coefficient: P'(yi|p) = P(yi|p)·MP(S,yi) / ∑P(yj|p)·MP(S,yj); where MP(S,yi) = ∑α(S)i / ∑α(S)j, representing the mixing ratio of class yi in patch S. Apply conditional random field smoothing, with the energy function: E(y) = ∑[1-P'(yi|p)] +β·∑[y(p)≠y(q)]; where β is the smoothing coefficient, set to 0.8, and y(p) and y(q) are the labels of adjacent pixels. Use a graph cut algorithm to optimize the label configuration.

[0169] By employing a multi-scale adaptive classification strategy, this embodiment can select the most suitable classifier and feature set based on the regional heterogeneity characteristics, solving the technical challenge that different heterogeneous regions require different classification strategies and improving the classification accuracy of mixed land cover regions.

[0170] 5.3. Evaluate the classification accuracy.

[0171] Based on experimental results from Xingyang City, the accuracy of this embodiment was compared with that of traditional methods. A stratified random sampling method was used to select validation samples, totaling 1500 sample points covering 9 land use types. Four indicators were calculated: overall accuracy (OA), Kappa coefficient, user accuracy (UA), and producer accuracy (PA). Comparisons were made with the following methods: traditional pixel-based maximum likelihood classification (MLC); object-based image analysis (OBIA); and multi-scale segmented classification (MSS). The results showed that the overall accuracy of this embodiment reached 87.6%, with a Kappa coefficient of 0.842; the overall accuracy of the traditional MLC method was 72.3%, with a Kappa coefficient of 0.667; the overall accuracy of the OBIA method was 78.9%, with a Kappa coefficient of 0.745; and the overall accuracy of the MSS method was 82.4%, with a Kappa coefficient of 0.793.

[0172] In the accuracy analysis of mixed land cover areas, the accuracy of this embodiment is 83.5% in the urban-rural transition zone, which is 15.7% higher than the OBIA method; the accuracy of this embodiment is 85.3% in the agroforestry ecotone, which is 12.9% higher than the OBIA method. This shows that this embodiment has significant advantages in handling heterogeneous landscape areas.

[0173] Experimental results demonstrate that the heterogeneous landscape land use mixed feature classification method based on homogeneous patch decomposition proposed in this embodiment can effectively improve the accuracy of land use classification in complex landscapes, especially showing significant advantages in highly heterogeneous urban-rural transition zones and agroforestry ecotones. By employing gradient structure feature extraction, a two-way feature transfer mechanism, and structural consistency index calculation, key technical challenges in heterogeneous landscape classification using traditional methods are addressed. The parameters in this embodiment achieved good results in the Xingyang City study area, but for other study areas, appropriate adjustments may be necessary based on specific circumstances. In particular, the scale parameters for multi-scale segmentation and the heterogeneity zoning threshold should be adjusted according to the complexity and resolution of the land features in the study area.

[0174] This embodiment uses Xingyang City, Henan Province as the study area and successfully implemented a heterogeneous landscape land use mixed feature classification method based on homogeneous patch decomposition. The main conclusions are as follows: Multi-scale hierarchical patch generation and patch hierarchy construction provide multi-level information support for heterogeneous landscape analysis, overcoming the problem that single-scale methods are insufficient to adapt to complex landscape features; gradient structure feature extraction successfully captures the structural and directional features within patches, especially demonstrating significantly better representation capabilities for complex feature boundaries and texture structures than traditional statistical features, showing a significant advantage in classification; the bidirectional feature transfer mechanism achieves effective information fusion between different scales, considering both global structure and local details, making it particularly suitable for... The classification accuracy was improved by an average of 8.7% when dealing with mixed land cover areas in heterogeneous landscapes. The structural consistency index can automatically determine the optimal analysis scale for each region, solving the problem of blind scale selection in traditional methods and significantly improving the accuracy of heterogeneous landscape classification. The sparse representation method based on dictionary learning successfully quantified the degree of land cover mixing in images, providing a new technical perspective for heterogeneous landscape classification, especially improving the classification accuracy of complex areas such as urban-rural transition zones by 15.7%. The multi-scale adaptive classification strategy selects the most suitable classifier and feature set according to the regional heterogeneity characteristics, solving the technical problem that different heterogeneous regions require different classification strategies, and the overall classification accuracy reached 87.6%.

[0175] This invention is applicable to the following scenarios: Dynamic monitoring of land use in urban fringe areas: accurately identifying subtle changes in urban-rural transition zones, providing a basis for decision-making in urban expansion analysis and planning. Ecological environment assessment in agroforestry ecotones: precisely distinguishing land use types in mixed-feature areas, providing technical support for ecological protection red line delineation and ecological compensation. Assessment of ecosystem services in complex landscapes: supporting the assessment of ecosystem service function value and ecological protection planning through refined land use classification. High-precision remote sensing mapping: improving the classification accuracy of existing land use classification systems, especially for more accurate representation of mixed-feature areas.

[0176] The promotional value of this invention lies in: Technical versatility: Although this embodiment is based on Gaofen-1 satellite data, the method can be extended to other medium- and high-resolution remote sensing images, such as GF-2, Sentinel-2, and WorldView. Computational efficiency: By using speckle decomposition and multi-scale fusion, the computational complexity of complex landscape classification is reduced, and the processing efficiency is improved by approximately 40% compared to pixel-level fusion analysis. Application flexibility: This invention can be integrated with existing classification systems and workflows to improve the accuracy and reliability of land use classification.

[0177] This invention addresses the blindness of scale selection through the Structural Consistency Index (SCI). This index quantitatively assesses the consistency between the internal structure of a patch and the structure of its sub-patterns, combining compositional complexity and heterogeneity indicators to determine the most suitable analysis scale for each region. In particular, the SCI calculation formula achieves a balanced assessment of patch structural integrity and complexity, thus avoiding over-segmentation or under-segmentation. By constructing a bidirectional feature transfer mechanism, effective information exchange between different scales is achieved. Top-down transfer transmits global contextual information from parent patches to sub-patterns, with the transfer formula ensuring effective downward transmission of global information while preserving local differences. Bottom-up transfer transmits local details from sub-patterns to parent patches, especially in highly heterogeneous regions, ensuring effective fusion of multi-scale information. A dictionary-based sparse representation method is introduced to accurately quantify mixed feature regions. By constructing a feature dictionary, each patch can be represented as a linear combination of dictionary elements. The sparsity coefficient directly reflects the proportion of feature composition, and the mixing index accurately quantifies the degree of mixing, providing a scientific basis for subsequent classification. An adaptive classification strategy based on heterogeneity was constructed, dividing regions into low, medium, and high heterogeneity levels, and applying different classifiers and feature sets to each level. For high-heterogeneity regions, a Gaussian mixture model was established using sparse coefficients and gradient structure features, and the classification accuracy for mixed land cover regions was improved by adjusting the classification probabilities. A gradient structure feature extraction method was proposed, analyzing the internal structural changes of image patches through multi-directional gradient operators, capturing directional and textural variation information that traditional statistical features cannot express. In particular, the directional, intensity, and coherence features obtained through gradient covariance matrix decomposition enable the system to accurately characterize the boundaries and internal structures of land covers.

[0178] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details of the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for classifying mixed land use features in heterogeneous landscapes based on homogeneous image patch decomposition, characterized in that, Includes the following steps: Read and segment the raw remote sensing image data to generate a multi-scale image patch set, construct the hierarchical relationship between the images, and obtain the image patch relationship map data; Based on multi-scale speckle sets and speckle relationship diagram data, we calculate the heterogeneity index, gradient structure features, and sub-spot distribution features within speckles to achieve feature transfer and fusion between different scales and obtain a coupled feature set. Read the multi-scale patch set and the original remote sensing image data, construct a feature dictionary and generate a sparse representation of the patches to obtain the sparse coefficient matrix; Based on the image patch relationship map data, coupled feature set, sparse coefficient matrix and pre-stored training sample data, the optimal feature scale is determined and adaptive classification is performed. Multi-scale results are fused and classification results are optimized to generate a fine land use classification map. The steps to achieve feature transfer and fusion across different scales to obtain a coupled feature set include: Based on the multi-scale patch set and the original remote sensing image data, the internal heterogeneity index of each patch is calculated to obtain the patch heterogeneity index matrix; Construct the gradient structure features for each image spot to obtain the gradient structure feature set; Based on multi-scale image patch set and image patch relationship map data, the distribution characteristics of child images in parent images are analyzed to obtain the child image patch distribution characteristic matrix; Based on the gradient structure feature set, the sub-pixel distribution feature matrix, and the speckle heterogeneity index matrix, feature transfer and fusion between different scales are achieved to obtain a coupled feature set.

2. The method according to claim 1, characterized in that, The steps for constructing a feature dictionary and generating a sparse representation of image spots to obtain a sparse coefficient matrix include: Based on the multi-scale patch set and the original remote sensing image data, the spectral, texture and shape features of each patch are extracted to form a patch feature vector set; Based on the image patch feature vector set, a land cover feature dictionary that can represent the feature patterns of different land cover types is created; The feature vector set of the image patch is sparsely decomposed using the feature dictionary of land cover. By optimizing the solution of the sparse coefficients, each image patch can be represented as a linear combination of dictionary elements, generating a sparse coefficient matrix that reflects the composition ratio of land cover. The degree of land cover mixing in each image patch is calculated based on the sparse coefficient matrix, which quantifies the mixing status of land cover types in heterogeneous landscapes and forms a land cover mixing index.

3. The method according to claim 1, characterized in that, The steps for generating a detailed land use classification map by integrating multi-scale results and optimizing classification results include: Read the image patch relationship map data and image patch heterogeneity index matrix, calculate the consistency between the internal structure of the image patch and the sub-image patch structure, determine the optimal analysis scale for each region, and form an optimal scale mapping map; By combining the optimal scale mapping map, coupled feature set, sparse coefficient matrix and training sample data, a differentiated classification strategy is implemented according to the degree of regional heterogeneity to obtain preliminary classification results; The spatial context is optimized based on the preliminary classification results to generate an optimized classification map; By comprehensively optimizing the classification map and the sparse coefficient matrix, and making fine-grained classification adjustments, a refined land use classification map is obtained.

4. The method according to claim 1, characterized in that, The steps to construct the gradient structure features for each image spot and obtain the gradient structure feature set include: Extract the image sub-regions corresponding to each image spot from the multi-scale image spot set to obtain the image spot sub-region set; Read the image spot sub-region set, construct a multi-directional gradient operator set, apply it to each band respectively, generate a multi-band gradient map set, and calculate the gradient intensity and direction of each image spot in each direction. Gradient covariance matrix is ​​constructed based on gradient intensity and direction to extract image spot gradient structure features; Based on the gradient structure characteristics of the image spot, the spatial distribution characteristics and structural change characteristics of the internal structure of the image spot are extracted to form a gradient structure feature vector. A gradient structure feature set is generated based on the gradient structure feature vector, the image heterogeneity index, and the adaptive enhancement coefficient.

5. The method according to claim 1, characterized in that, The steps to achieve feature transfer and fusion across different scales to obtain a coupled feature set also include: Based on multi-scale speckle sets and speckle relationship graph data, the spatial overlap and feature similarity between speckles of different scales are calculated, and a scale relationship weight matrix is ​​constructed. Based on the scale relationship weight matrix, gradient structure feature set, sub-pixel distribution feature matrix, and speckle heterogeneity index matrix, top-down feature transfer is achieved, enabling parent speckle features to be transferred to child speckles, and the transfer intensity is adjusted according to the differences between parent and child speckle features to obtain the down-transferred feature set; bottom-up feature transfer is implemented, enabling child speckle information to be transferred to parent speckles to obtain the up-transferred feature set. Read the downlink feature set and the uplink feature set, fuse the bidirectional transmitted features, and normalize and reduce their dimensionality to obtain the coupled feature set.

6. The method according to claim 3, characterized in that, The steps to determine the optimal analysis scale for each region and generate the optimal scale mapping include: Analyze the patch relationship map data, identify the sub-patterns contained in each patch and their spatial distribution characteristics, and generate a simplified patch relationship map; Based on the simplified image patch relationship diagram, the structural feature similarity between the image patch and its sub-image patches is calculated to form a structural similarity matrix; the quantity diversity, spatial distribution and type differences of the sub-image patches are analyzed to assess the compositional complexity within the image patch and establish a compositional complexity matrix; the structural similarity matrix and the compositional complexity matrix are combined to construct a structural consistency index matrix. Analyze the structural consistency index matrix and the patch heterogeneity index matrix to determine the scale corresponding to the maximum structural consistency for each region and generate the optimal scale mapping map.

7. The method according to claim 3, characterized in that, The steps to obtain preliminary classification results by implementing a differentiated classification strategy based on the degree of regional heterogeneity include: By reading the coupled feature set, sparse coefficient matrix, optimal scale mapping map and training sample data, the image patch regions are divided into three categories: low heterogeneity, medium heterogeneity and high heterogeneity. Constructing a hierarchical classifier: For low heterogeneity regions, a random forest classifier based on coupling features is used; for medium heterogeneity regions, a support vector machine combining coupling features and sparse coefficients is used; and for high heterogeneity regions, a Gaussian mixture model is built using sparse coefficients and gradient structure features, resulting in a hierarchical classification model set. Based on the optimal scale mapping map and hierarchical classification model set, the optimal scale patch to which each pixel belongs and its corresponding classifier are determined, and pixel classification label map and classification probability map are generated for highly heterogeneous regions. By combining the pixel classification label map, classification probability map, and sparse coefficient matrix, the proportion of mixed land features is calculated, the classification probability of highly heterogeneous areas is adjusted, and preliminary classification results are formed.

8. The method according to claim 5, characterized in that, The steps for constructing the scale-related weight matrix and implementing top-down feature transfer include: The formula for constructing the scale relationship weight matrix is: W(p,c) = OR(p,c)·FS(p,c); Where OR(p,c) represents spatial overlap and FS(p,c) represents feature similarity; When implementing top-down feature transfer, a differentiated transfer strategy is adopted. The transfer intensity is proportional to the difference in features between parent and child image patches. The transfer formula is: FD(ci) = F(ci) + γD·WD(p,ci)·[F(p) – F*(C(p))]; Where FD(ci) is the feature vector of sub-image spot ci after top-down feature transfer; F(ci) is the original feature vector of sub-image spot ci; γD is the downlink adjustment coefficient, which controls the intensity of top-down transfer; WD(p,ci) is the transfer weight from parent image spot p to sub-image spot ci; F(p) is the feature vector of parent image spot p; F*(C(p)) is the average feature of all sub-image spots of parent image spot p, and C(p) represents the set of all sub-image spots of p.

9. The method according to claim 6, characterized in that, The steps to obtain the composition complexity matrix and construct the structure consistency index matrix include: By combining the sub-image patch quantity diversity index, spatial distribution entropy index, and type diversity index, a comprehensive evaluation is formed through weighted fusion, resulting in a composition complexity matrix. The formula for calculating the structural consistency index matrix is: SCI(S) = ω(S)·SSavg(S) - (1-ω(S))·CC(S); Where ω(S) is the adaptive weight function, SSavg(S) is the average structural similarity matrix, and CC(S) is the compositional complexity.

Images