A method for identifying land cover in a highly heterogeneous region

Abstract

This invention discloses a method for identifying land cover in highly heterogeneous areas. The method includes: collecting heterogeneous historical land cover products of the target area and obtaining stable regions of different cover types through overlay analysis; dividing the target area into blocks based on its latitude and longitude information with a spatial constraint of 1°×1°, then setting three thresholds of 50, 70, and 100, and using a stratified sampling strategy to obtain uniformly distributed sample data of the target area; acquiring temporal optical, radar, and digital elevation data of the target area, and addressing the influence of clouds and rain through data synthesis; and finally constructing a tile-based classification model for classification. The tile-based classification model constructed in this invention addresses the heterogeneity of plateau regions, fully considering the similarity of adjacent land features and the large climate span of the plateau. Simultaneously, a strategy for stable, rapid, and uniform sample data generation is developed, saving material and financial resources. The mapping results provide data support for ecological civilization construction and biodiversity conservation.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image recognition, and specifically to a method for identifying land cover in highly heterogeneous areas. Background Technology

[0002] Land cover mapping plays a crucial role in natural resource management. It provides essential information for scientific research in areas such as climate change, biodiversity conservation, ecosystem protection, and urban planning. The emergence and development of remote sensing and computer science technologies have supported large-scale land cover mapping. With the advent of the Google Earth Engine (GEE) cloud platform, users can extract regional or national land cover information on their personal computers, and image compositing methods can overcome the problem of low availability of optical data in plateau regions.

[0003] Generally, the efficiency and quality of remote sensing image classification depend on image quality, sample representativeness, and method availability. Supervised classification methods heavily rely on real ground sample data, which is a major drawback of large-scale land cover mapping. Collecting real ground sample data requires significant human and material resources. While unsupervised classification does not require sample data, its poor classification performance limits its application in large-scale land cover mapping. Furthermore, influenced by a combination of factors such as genetic characteristics, environmental factors, and cultivation conditions, the spectral characteristics of ground features exhibit complex characteristics such as time lag and abrupt changes in large-scale, heterogeneous spaces, leading to biases in the representation of key phenological information of crops in remote sensing images. A key issue is how to effectively consider large-scale regional differences and accurately identify land cover types in plateau regions using discrete remote sensing spatiotemporal image data and representative sample data. This, in turn, ensures that the identification model used has good universality and generalization in plateau regions. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method for identifying land cover in highly heterogeneous areas. To solve the aforementioned technical problems, this invention adopts the following technical solution:

[0005] A method for identifying land cover in highly heterogeneous areas is characterized by the following steps:

[0006] Step 1: Collect heterogeneous historical land cover data for the target area, and obtain stable areas for each land cover category in the target area through overlay analysis;

[0007] Step 2: Based on the latitude and longitude of the target area, generate a 1°×1° rectangular grid for spatial constraint to obtain the tile information of the target area;

[0008] Step 3: Use a stratified block sampling strategy to automatically generate sample data for each grid with different coverage categories. Three basic thresholds for sample generation were designed, namely 50, 70 and 100.

[0009] Step 4: Select optical remote sensing data with a cloud pixel ratio of less than 10% for the target year, perform cloud removal processing, and then extract spectral features and index features;

[0010] Step 5: Filter radar data for the target year and preprocess it. Overlay the time-series synthesized optical data, radar data and DEM elevation data to obtain a composite image. Finally, extract polarization features and terrain features.

[0011] Step 6: Combining the spectral features, exponential features, polarization features, and topographic features from Steps 4 and 5, a recursive feature elimination algorithm is used to optimize features and obtain the best feature combination.

[0012] Step 7: Select the best features from Step 6 and the sample data from Step 3, and use a tile-based classification model to classify the data to obtain a spatial distribution map of land cover information in the target area.

[0013] Preferably, in step 1, a stable region refers to a region where the type of cover that grows in the area has remained stable for many years without any replacement.

[0014] As a preferred option, in step 3, the stratified sampling strategy generates sample data using the stratifiedSample() function within the stable region, based on different cover types, and iterates through all 1°×1° grid tiles.

[0015] As a preferred option, the specific process of processing and synthesizing optical data in step 4 is as follows: screening optical images with cloud cover of less than 10%; performing cloud masking on the optical data; and synthesizing the cloud-masked images into a clear, cloudless image of the target area through median synthesis.

[0016] Preferably, in step 4, the calculated spectral features include shortwave infrared band, red-edge bands 1-3, red band, green band, blue band, and near-infrared band; the calculated index features include normalized difference vegetation index, ground chlorophyll vegetation index, enhanced vegetation index, normalized difference water index, improved normalized difference water index, and surface water index (LSWI), and the calculation formulas are as follows:

[0017]

[0018]

[0019]

[0020]

[0021]

[0022]

[0023] Wherein, mir represents the mid-infrared band; swir2 represents the second short-wave infrared band.

[0024] Preferably, in step 5, the polarization characteristics are the radar backscattering coefficients in the VV and VH bands; the terrain characteristics include the ratio, second-order angular moment, and entropy, calculated using the following formulas:

[0025]

[0026]

[0027]

[0028] Where n represents the order of the gray-level co-occurrence matrix, i,j represent the pixel determinant coordinate information in the image, and P(i,j) represents the gray-level joint probability matrix.

[0029] Preferably, in step 6, the recursive feature elimination algorithm evaluates the importance of classification features based on the Gini coefficient, uses all features for RFE, and employs five-fold cross-validation for accuracy evaluation, as shown in the following formula:

[0030]

[0031]

[0032] Where P(j / h) represents the probability that the test sample h belongs to jth when sampling from the sample data; n(h) represents the total number of training samples under the test sample h; nj(h) represents the total number of test samples h belonging to jth; and j represents the number of categories.

[0033] As a preferred embodiment, step 7, based on the tile classification model, involves the following steps: extracting sample information from the central region and its surrounding 8 neighborhoods using a 1°×1° grid of tiles; classifying the central tile using all sample information within a 3×3 tile and a random forest (RF) algorithm; performing this operation on each tile, and then calling the mosaic() function to stitch the tiles together; finally, obtaining the land cover mapping results for the target area.

[0034] The advantages of this invention over the prior art are:

[0035] (1) This method is simple and practical. It can be implemented by coding on the GEE and Python platforms to extract large-scale land cover information, saving the complicated and cumbersome operation steps of single-machine operation.

[0036] (2) This method saves manpower, material resources and financial resources for obtaining sample data. By superimposing analysis of heterogeneous historical land cover products and introducing 1°×1° tiles, sample data can be generated automatically, stably and quickly, avoiding local clustering of data.

[0037] (3) This method solves the problem of heterogeneity in plateau areas. The tile-based classification model fully considers the spectral similarity of geographically adjacent features, thus improving the limitations of conventional classification. Attached Figure Description

[0038] Figure 1 This is a schematic diagram of the overall process of the present invention;

[0039] Figure 2 This is a portion of the sample data generated by an example of the present invention;

[0040] Figure 3 This is the result of remote sensing mapping of land cover information at a resolution of 10 meters for the target area of ​​this invention. Detailed Implementation

[0041] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0042] Example 1

[0043] A method for identifying land cover in highly heterogeneous areas includes the following steps:

[0044] Step 1: Collect heterogeneous historical land cover data for the target area, and obtain stable areas for each land cover category in the target area through overlay analysis;

[0045] A stable region refers to an area where the land cover type has remained stable over many years without any change. Specifically, if the historical land cover data from different sources over many years all belong to the same land type, then the intersection of these data points is defined as a stable region.

[0046] This example uses five data points from four historical land cover products covering 2015-2020. By overlaying the same land type from different data points geographically, the intersecting part is the source of the sample database.

[0047] Step 2: Based on the latitude and longitude of the target area, generate a 1°×1° rectangular grid for spatial constraint to obtain the tile information of the target area, totaling 56 tiles;

[0048] Step 3: Use a stratified block sampling strategy to automatically generate sample data for each grid with different coverage categories. Three basic thresholds for sample generation were designed, namely 50, 70 and 100.

[0049] The stratified sampling strategy generates sample data separately within a stable region using the `stratifiedSample()` function based on different cover types, iterating through all 1°×1° grid tiles. Training data generated with thresholds of 50, 70, and 100 are 19600, 27440, and 39200 respectively. Some sample data are shown below. Figure 2 As shown.

[0050] Step 4: Select optical remote sensing data with a cloud pixel ratio of less than 10% for the target year, perform cloud removal processing, and then extract spectral features and index features;

[0051] The specific process of processing and synthesizing optical data is as follows: screening optical images with cloud cover of less than 10%; performing cloud masking on the optical data; and using the median() function to synthesize a clear, cloudless image of the target area from the cloud-masked image through median synthesis.

[0052] The calculated spectral characteristics include shortwave infrared band, red-edge bands 1-3, red band, green band, blue band, and near-infrared band; the calculated index characteristics include normalized difference vegetation index (NDVI), ground chlorophyll vegetation index (GCVI), enhanced vegetation index (EVI), normalized difference water index (NDWI), modified normalized difference water index (MNDVI), and surface water index (LSWI), and the calculation formulas are as follows:

[0053]

[0054]

[0055]

[0056]

[0057]

[0058]

[0059] Wherein, mir represents the mid-infrared band; swir2 represents the second short-wave infrared band.

[0060] Step 5: Filter radar data for the target year and preprocess it. Overlay the time-series synthesized optical data, radar data and DEM elevation data to obtain a composite image. Finally, extract polarization features and terrain features.

[0061] Polarization characteristics are represented by radar backscattering coefficients in the VV and VH bands; terrain features include specificity (CONTRAST), second-order angular moment (ASM), and entropy (ENT), calculated using the following formulas:

[0062]

[0063]

[0064]

[0065] Where n represents the order of the gray-level co-occurrence matrix, i,j represent the determinant coordinates of pixels in the image, and P(i,j) represents the gray-level joint probability matrix. Data acquisition, processing, and cropping are all implemented in the GEE remote sensing big data platform, which is fast and efficient, avoiding the complex operations of stand-alone processing.

[0066] Step 6: Combining the spectral features, exponential features, polarization features, and topographic features from Steps 4 and 5, use the Recursive Feature Elimination (RFE) algorithm to perform feature optimization and obtain the optimal feature combination;

[0067] The recursive feature elimination algorithm evaluates the importance of categorical features based on the Gini coefficient, applies all features to the RFE, and uses five-fold cross-validation to evaluate accuracy, as shown in the following formula:

[0068]

[0069]

[0070] Where P(j / h) represents the probability that test sample h belongs to jth when sampling from the sample data; n(h) represents the total number of training samples under test sample h; nj(h) represents the total number of test samples h belonging to jth; and j represents the number of classes. The recursive feature elimination algorithm was implemented in Python 3.7. Sample data generated with three thresholds of 50, 70, and 100 were used as input to the algorithm, obtaining three optimal feature combination spaces. Simultaneously, we used 5x cross-validation to evaluate the best prediction accuracy of the classifier. Figure 3 As shown, when performing RFE using sample pools with thresholds of 50 and 70, the test accuracy increases with the number of features, eventually stabilizing. When the number of features reaches 8, the test accuracy does not significantly improve; the test accuracy is highest when the number of features reaches 18. When performing RFE using a sample pool with a threshold of 100, the test accuracy does not significantly improve when the number of features reaches 11; the test accuracy is highest when the number of features reaches 22, and then gradually decreases with further increases in the number of features.

[0071] Step 7: Select the best features from Step 6 and the sample data from Step 3, and use a tile-based classification model to classify the data to obtain a spatial distribution map of land cover information in the target area.

[0072] The specific steps based on the tile classification model are as follows: Based on a 1°×1° grid of tiles, extract sample information from the central area and its surrounding 8 neighborhoods; use all sample information within a 3×3 tile and the Random Forest (RF) algorithm to classify the central tile; perform this operation on each tile, and then call the mosaic() function to stitch them together; finally, obtain the land cover mapping results for the target area.

[0073] To assess the accuracy of the land cover map, we performed pixel-by-pixel validation and quantitative analysis using validation samples, calculating the overall accuracy (OA), production accuracy (PA), user accuracy (UA), and F1 score. As shown in Table 1, overall, the OA reached 91.92%, and the more stringent F1 score reached 88.39%, which is relatively high accuracy for a 10m provincial land cover map.

[0074] Table 1. Accuracy Evaluation of 10m Land Cover Map of Yunnan Plateau

[0075]

[0076] Although the invention has been described herein with reference to several illustrative embodiments, it should be understood that many other modifications and implementations can be devised by those skilled in the art, which will fall within the scope and spirit of the principles disclosed herein. More specifically, various variations and modifications can be made to the components and / or layout of the subject matter arrangement within the scope of the disclosure, drawings, and claims. Besides variations and modifications to the components and / or layout, other uses will be apparent to those skilled in the art.

Claims

1. A method for identifying land cover in highly heterogeneous regions, the method comprising: The method includes the following steps: ​ Step 1: Collect heterogeneous historical land cover data for the target area, and obtain stable areas for each land cover category in the target area through overlay analysis; Step 2: Based on the latitude and longitude of the target area, generate a 1°×1° rectangular grid for spatial constraint to obtain the tile information of the target area; Step 3: Use a stratified block sampling strategy to automatically generate sample data for each grid with different coverage categories. Three basic thresholds for sample generation were designed, namely 50, 70 and 100. Step 4: Select optical remote sensing data with a cloud pixel ratio of less than 10% for the target year, perform cloud removal processing, and then extract spectral features and index features; Step 5: Filter radar data for the target year and preprocess it. Overlay the time-series synthesized optical data, radar data and DEM elevation data to obtain a composite image. Finally, extract polarization features and terrain features. Step 6: Combining the spectral features, exponential features, polarization features, and topographic features from Steps 4 and 5, a recursive feature elimination algorithm is used to optimize features and obtain the best feature combination. Step 7: Select the best features from Step 6 and the sample data from Step 3, and use a tile-based classification model to classify the data to obtain a spatial distribution map of land cover information in the target area.

2. The method of claim 1, wherein: In step 1, a stable region refers to a region where the type of cover that grows in the area has remained stable for many years without any replacement.

3. The method of claim 1, wherein: In step 3, the stratified sampling strategy generates sample data using the stratifiedSample() function within the stable region, based on different cover types, and iterates through all 1°×1° grid tiles.

4. The method of claim 1, wherein: In step 4, the specific process of processing and synthesizing optical data is as follows: screening optical images with cloud cover of less than 10%; performing cloud masking on the optical data; and synthesizing the cloud-masked images into a clear, cloudless image of the target area through median synthesis.

5. The method of claim 1, wherein: In step 4, the calculated spectral features include shortwave infrared band, red-edge bands 1-3, red band, green band, blue band, and near-infrared band; the calculated index features include normalized difference vegetation index, ground chlorophyll vegetation index, enhanced vegetation index, normalized difference water index, improved normalized difference water index, and surface water index (LSWI), and the calculation formulas are as follows: Wherein, mir represents the mid-infrared band; swir2 represents the second short-wave infrared band.

6. The method of claim 1, wherein: In step 5, the polarization characteristics are the radar backscattering coefficients VV and VH bands; the terrain characteristics include the ratio, second-order angular moment, and entropy, calculated using the following formulas: Where n represents the order of the gray-level co-occurrence matrix, i,j represent the pixel determinant coordinate information in the image, and P(i,j) represents the gray-level joint probability matrix.

7. The method of claim 1, wherein: In step 6, the recursive feature elimination algorithm evaluates the importance of classification features based on the Gini coefficient, uses all features for RFE, and employs five-fold cross-validation to evaluate accuracy, as shown in the following formula: Where P(j / h) represents the probability that the test sample h belongs to jth when sampling from the sample data; n(h) represents the total number of training samples under the test sample h; nj(h) represents the total number of test samples h belonging to jth; and j represents the number of categories.

8. The method of claim 1, wherein: In step 7, the specific steps based on the tile classification model are as follows: Based on 1°×1° grid tiles, extract sample information of the central area and its surrounding 8 neighborhoods; use all sample information within 3×3 tiles and the random forest algorithm to classify the central tile; perform this operation on each tile, and then call the mosaic() function to stitch them together; finally, obtain the land cover mapping result of the target area.

Images