Crop planting structure extraction method based on hierarchical extraction and multi-feature integration

Abstract

The invention relates to the field of remote sensing application, and particularly discloses a crop planting structure extraction method based on hierarchical extraction and multi-feature integration. The method comprises the steps of data acquisition and preprocessing, vegetation region extraction, classification model training and crop structure classification, and performs vegetation and non-vegetation division on a target region by adopting an OTSU algorithm and an NDVI hierarchical extraction method, so that the sample demand quantity is reduced, and the crop structure classification efficiency is improved. And the random forest classification model is trained by using the sample data, the spectral features, the vegetation features and the texture features, so that the precision and reliability of model classification are improved.

Description

Technical field: The present invention relates to the field of remote sensing applications, and in particular to a method for extracting crop planting structures based on hierarchical extraction and multi-feature integration. Background technology: Traditional crop statistics rely on ground surveys, which are costly and subjective, making them inadequate for modern agriculture. Satellite remote sensing technology can efficiently monitor changes in farmland and provide more accurate crop information. Crop classification research falls into two main categories: one uses single-phase satellite spectral information for classification, but this can be limited in accuracy due to similar spectral characteristics; the other uses multi-phase satellite data, leveraging temporal variation to improve classification accuracy and is suitable for areas with similar crop growth periods and seasonality.

[0003] In the prior art, such as Chinese Patent Publication No. CN119540780B, a monitoring and early warning method for the distribution of Spartina alterniflora based on image recognition, discloses a method for extracting vegetation features using Sentinel-1 radar data and Sentinel-2 optical remote sensing data. This technical solution can only identify a single type of vegetation. The spectrum and phenological characteristics of Spartina alterniflora are relatively stable, that is, the spectral characteristic difference of the vegetation of Spartina alterniflora in one year is relatively small, and the spectral characteristic difference with other vegetation around it is relatively large. For wheat and corn, the spectral characteristic similarity is higher during the growth period, such as the green band reflectivity of the heading stage and the jointing stage is close, but the phenological characteristics of wheat and corn are more complicated to change, relying only on a single phase or simple multi-phase classification to easily cause misjudgment. In addition, this method relies on Sentinel-1 radar data, and although it is sensitive to vegetation structure, it is limited in distinguishing the subtle phenological changes of crops, such as the expansion of corn leaves and the tillering of wheat, and classification robustness can decline when optical data is missing in cloudy areas.

[0004] In addition, Chinese Patent Publication No. CN119992184A, a method and system for extracting peanut planting distribution based on phenological characteristics, discloses a method for extracting peanut planting distribution using Sentinel-2 image data. This method relies on simple vegetation indices such as NDVI for classification. However, the NDVI curves of wheat and corn highly overlap during the mid-growth period, making it difficult to effectively distinguish them using spectral threshold segmentation alone. Furthermore, if this method does not fully utilize the temporal characteristics of the crop's entire growth period (such as the rapid canopy closure of corn and the slow tillering dynamics of wheat), the classification accuracy will be limited by the similarity of the growth period, making it difficult to meet the needs of large-scale, precise monitoring. Summary of the invention: The purpose of the present invention is to provide a crop planting structure extraction method based on hierarchical extraction and multi-feature integration to solve the problems existing in the prior art.

[0006] The present invention is implemented by the following technical solution: a method for extracting crop planting structure based on hierarchical extraction and multi-feature integration, comprising the following steps: S1 Data acquisition and preprocessing: Based on the GEE platform, Sentinel-2A satellite image data of the two target areas were downloaded and de-clouded to obtain the de-clouded image data of the two target areas. S2 extracts vegetation areas: Calculate the NDVI frequency distribution histogram for the two cloud-removed image data obtained in step S1, determine the optimal thresholds for the two images, fuse the two images using the OTSU algorithm with a maximum value strategy, mask them after fusion to obtain fused image data, determine the vegetation area in the fused image data, and obtain image data with the vegetation area marked; S3 classification model training: Sample data is obtained by combining GNSS RTK positioning with drone aerial photography. The sample data, spectral features, vegetation features, and texture features are used as model input information and fed into a random forest (RF) model for training. The trained random forest (RF) model is obtained. S4 Crop structure classification: The image data with the marked vegetation area obtained in step S2 is input into the trained Random Forest (RF) model to classify the vegetation areas marked in the image and calculate the area, thereby determining the crop structure type and area of ​​each crop in the vegetation area in the image data with the marked vegetation area.

[0007] Furthermore, the declouding processing of the image data in step S1 includes calculating the values ​​of Bit10 and Bit11 of the band QA60, setting the image when the values ​​of Bit10 and Bit11 are both 0 as a cloud mask, using the cloud mask to remove cloud information in the image, obtaining a cloud-free image, and cropping the image according to the target area; if there is no image that meets the conditions, screening images with a cloud amount less than 10 and performing declouding processing, cropping according to the target area and then splicing to obtain declouded image data of the target area for two periods.

[0008] Furthermore, in step S3, the B2, B3, B4, B5, B6, B7, B8, B8A, B11, and B12 bands of Sentinel-2A are selected as spectral features.

[0009] Furthermore, in step S3, some band data from Sentinel-2A are selected to construct vegetation index as vegetation features, wherein the vegetation features include: Normalized Difference Vegetation Index (NDVI), which is calculated as follows: , Normalized Difference Vegetation Index 1 (NDVI1), which is calculated as , Normalized Difference Vegetation Index 2 (NDVI2), the calculation formula is , Normalized Difference Vegetation Index 3 (NDVI3), the calculation formula is , the red edge normalized vegetation index 1 (NDre1), which is calculated as , the red edge normalized vegetation index 2 (NDre2), which is calculated as , Green Normalized Difference Vegetation Index (GNDVI), which is calculated as , Land Surface Water Index (LSWI), which is calculated as , bare soil index (BSI), which is calculated as , the modified chlorophyll absorption reflectance index (MCARI), which is calculated as follows: , Enhanced Vegetation Index (EVI), which is calculated as , Difference Vegetation Index (DVI), which is calculated as , Soil Adjusted Vegetation Index (SAVI), which is calculated as , Ratio Vegetation Index (RVI), which is calculated as , Ratio Vegetation Index 1 (RVI1), which is calculated as , Green Chlorophyll Vegetation Index (GCVI), which is calculated as .

[0010] Furthermore, in step S3, texture features are extracted based on the B8 band in the Sentinel-2A data through the gray level co-occurrence matrix (GLCM) method, which includes calculating the second-order moment, mean sum, variance, contrast, difference, correlation, entropy and inverse moment of the B8 band.

[0011] Furthermore, the number of decision trees (ntree parameter) of the random forest model in step S3 is set to 120.

[0012] Advantages of the present invention: 1. The method of the present invention is based on the OTSU algorithm and the NDVI layered extraction method to reduce sample dependence and improve data processing efficiency. By using the OTSU algorithm and the NDVI layered extraction method for vegetation and non-vegetation segmentation, it can reduce sample requirements and improve classification efficiency. The strategy of first masking vegetation and then fine-classifying is adopted to reduce non-vegetation noise interference and improve classification efficiency. Selecting the characteristic data of two periods of vegetation for masking can eliminate misclassification and classification efficiency caused by similar phenological characteristics between the classified vegetation types.

[0013] 2. Use the Random Forest feature importance selection algorithm to reduce data dimensionality. Random Forest plays an important role in remote sensing image classification and can perform both feature selection and dimensionality reduction. The importance of different feature variables is calculated using the out-of-bag error rate, which is then used for feature selection. Compared to traditional classification methods (such as maximum likelihood), it can better capture nonlinear feature relationships and improve classification accuracy.

[0014] 3. The sample data, 10 spectral features, 16 vegetation features, and 8 texture features were fed into a random forest as training data. The importance score of each feature variable was calculated to assess its contribution to the classification results. Because the random nature of the random forest (RF) model can cause fluctuations in feature importance scores, five iterations were performed to minimize the impact of randomness on the results. By integrating multi-dimensional features (spectral, vegetation, and texture) and a hierarchical classification strategy, the accuracy of crop recognition was significantly improved. Description of the drawings: Figure 1 It is a flow chart of the technical solution of the present invention.

[0016] Figure 2 Schematic diagram of the target crop growth stages.

[0017] Figure 3 This is the frequency distribution histogram of the cloud-removed image data for the two periods.

[0018] Figure 4 This is the inter-class variance result diagram of the two-phase cloud-removed image data.

[0019] Figure 5 This is a diagram showing the classification results of vegetation and non-vegetation according to an embodiment of the present invention.

[0020] Figure 6 This is a diagram of the training process of the random forest model.

[0021] Figure 7 This is a diagram of crop classification results according to an embodiment of the present invention. Specific implementation method: To specifically illustrate the crop planting structure extraction method based on hierarchical extraction and multi-feature integration of the present invention, the present invention is further described below in conjunction with an embodiment. The specific flow chart is as follows: Figure 1 As shown: A method for extracting crop planting structure based on hierarchical extraction and multi-feature integration is based on the GEE platform and directly calls Sentinel-2A data. The data has undergone orthorectification, geometric correction, and atmospheric correction. The Hetao Irrigation District in Bayannur City, western Inner Mongolia Autonomous Region, China, is the target observation area of ​​this embodiment. According to preliminary investigations and screening, the types of crops planted in this area are mostly corn, sunflower, wheat, and tomato. Since wheat is in the jointing stage in May and other crops are still in the seedling stage, wheat can be extracted using images from May. Corn and sunflower are in the milky stage in August, and the two crops have obvious spectral differences in the images. Tomatoes have reached the harvest period in August, so these three crops can be accurately distinguished in August (such as Figure 2 Therefore, at least two periods of imagery are required for classification. Based on crop information and image quality, the final time ranges are May 26-29, 2022, and August 8-9, 2022, respectively.

[0023] S1 Data Acquisition and Preprocessing: After determining the target area and the time interval between the two imagery phases, Sentinel-2A satellite imagery of the Hetao Irrigation District from May 26 to May 29, 2022, and from August 8 to September 9, 2022, was downloaded from the GEE platform. Cloud removal was performed on each of the two imagery phases. Sentinel-2A, launched in 2015 by the European Space Agency (ESA), is a high-resolution optical satellite. It is the only satellite to include three red-edge bands, making it advantageous for monitoring vegetation characteristics. Its Multispectral Imager (MSI) includes 13 bands and offers spatial resolutions of 10, 20, and 60 meters (Table 1). It has a swath width of up to 290 kilometers and a revisit time of five days. This provides continuous time series data for surface monitoring and analysis.

[0024] Table 1 Spectral bands of Sentinel-2 sensors

[0025] Image cloud removal: Bits 10 and 11 of band QA60 were calculated on the GEE platform and set to 0 to create a cloud mask. This mask was used to remove cloud information from the image, resulting in a cloud-free image. The image was then cropped based on the study area. If an image did not meet the criteria, images with a cloud cover less than 10 were screened and clouded. The images were then cropped based on the study area and stitched together. A total of 37 images were selected, with specific dates shown in Table 2. Table 2 Time and number of images of irrigation areas in 2022

[0026] S2 Extract vegetation area: Calculate the NDVI frequency distribution histogram for the two periods of cloud-free image data obtained in step S1 (e.g. Figure 3 ), Figure 3 (a) is the first phase image and (b) is the second phase image. Based on this, the inter-class variance result diagram is obtained by exhaustive calculation of the intra-class variance (such as Figure 4 ), Figure 4 (c) is the first-period image, and (d) is the second-period image. The optimal thresholds for the influence of the first-period image and the second-period image are 0.22 and 0.33, respectively. The two-period image is fused using the OTSU algorithm with the maximum strategy. After fusion, it is masked to obtain the fused image data. The fusion of the two-period image data can eliminate the misjudgment problem caused by the small differences in the spectral characteristics of the various vegetation observed in the same time period. For example, the jointing period of corn and the heading period of wheat are close in time, and the reflectance of the green band is similar. If only a single-period image data of this time period is selected, the classification effect of corn and wheat will be poor. If the vegetation to be identified and classified is more complex and diverse, more periods of image data can be selected for fusion to reduce the influence of the close spectral characteristics between vegetation. The vegetation area in the fused image data is determined. Among them, the first-period image with a value greater than 0.22 is a vegetation area, and less than 0.22 is a non-vegetation area. The second-period image with a value greater than 0.33 is a vegetation area, and less than 0.33 is a non-vegetation area. Figure 5 This figure shows the spatial distribution of vegetation and non-vegetation within the Hetao Irrigation District. To illustrate the classification results, four specific areas were selected: areas a, b, c, and d represent bare land, residential areas, roads, and water bodies, respectively. Vegetation is represented in green, and non-vegetation in blue, clearly highlighting the differences between the two types of land features.

[0027] S3 classification model training: Sample data is obtained by combining GNSS RTK positioning with drone aerial photography. The sample data, spectral features, vegetation features, and texture features are used as model input information and fed into a random forest (RF) model for training. The trained random forest (RF) model is obtained. The planting structure classification sample data were obtained from June to August 2022 by combining GNSS RTK positioning and drone aerial photography. The main collection types and the number of samples of each type are shown in Table 3.

[0028] Table 3: Sample size type information

[0029] Others include woodlands, grasslands, and other crop types. The survey primarily covered land use types, crop species and growth cycles, on-site photographs, and coordinate information. Sample points for each land use type were selected randomly, evenly, and as widely distributed as possible across the entire area. All samples were polygonal in shape, resulting in a total of 1,070 samples. These samples were divided into training and precision verification samples based on the polygonal shape of the plots, with a ratio of 7:3.

[0030] The B2, B3, B4, B5, B6, B7, B8, B8A, B11, and B12 bands in Sentinel-2A were selected as spectral features. The spectral features are the original bands B2, B3, B4, B5, B6, B7, B8, B8A, B11, and B12 data in Sentinel-2A, which can be downloaded and used directly. The B1, B9, and B10 bands that are not suitable for crop classification are excluded.

[0031] In step S3, some band data from Sentinel-2A are selected to construct vegetation indices as vegetation features. The selected vegetation features and their calculation formulas are shown in Table 4.

[0032] Table 4

[0033] Vegetation characteristics can be used to quantitatively and qualitatively assess vegetation cover and growth. NDVI and EVI are the most commonly used indices, reflecting comprehensive changes in crop growth. NDVI1, NDVI2, NDVI3, NDre1, and NDre2 are derived from Sentinel-2's unique red-edge bands (B5, B6, and B7), enabling more accurate extraction of vegetation information. Different vegetation indices are suitable for monitoring specific crop characteristics. For example, GNDVI is commonly used to detect wilting in dense vegetation or mature crops. LSWI accurately reflects vegetation moisture content and is suitable for drought monitoring. BSI distinguishes crop from non-crop areas and assesses soil exposure. MCARI analyzes chlorophyll content. GCVI addresses NDVI saturation in high-cover areas. RVI is suitable for estimating green biomass. SAVI reduces the impact of soil background on NDVI. DVI is sensitive to changes in soil characteristics and is more suitable for low to medium vegetation cover. The introduction of these indices overcomes the limitations of single spectral features and improves classification accuracy in areas with complex crop structures.

[0034] Texture features are extracted based on the B8 band in the Sentinel-2A data using the gray-level co-occurrence matrix (GLCM) method, which includes calculating the second-order moment, mean sum, variance, contrast, dissimilarity, correlation, entropy, and inverse moment of the B8 band. The eight most common texture features are obtained through the above calculations.

[0035] The aforementioned sample data, along with 10 spectral features, 16 vegetation features, and 8 texture features, were fed into a random forest. The importance score of each feature variable was calculated to assess its contribution to the classification results. Because the random nature of the random forest (RF) model can cause fluctuations in feature importance scores, five iterations were performed to minimize the impact of randomness on the results. Finally, feature variables were ranked and selected based on their average importance scores, ensuring robust and reliable feature selection.

[0036] Classification Method: Random Forest (RF) is an ensemble classifier composed of multiple decision trees. Its main idea is to improve the model's generalization ability and classification accuracy by constructing multiple decision trees and aggregating their predictions. RF is an evolution of the traditional decision tree, used to classify new data by taking a majority vote on the classification results of all constructed decision trees. In RF, each node is divided by the optimal combination of a randomly selected subset of feature variables. This randomness effectively reduces the risk of model overfitting. When distinguishing crops such as corn, wheat, sunflower, and tomato, each decision tree may focus on different feature combinations (such as vegetation index or texture features) to classify the crops. This ensemble approach not only improves classification accuracy but also enhances the model's robustness, enabling it to effectively cope with data noise and complex feature distributions. In the RF model, ntree and mtry are two key hyperparameters. The ntree parameter determines the number of decision trees in the ensemble learning, while the mtry parameter controls the number of randomly selected features in the distribution for each decision tree node. According to relevant research, the recommended value of the mtry parameter is the square root of the total number of features. This value is also the default setting of the RF classifier in the GEE platform. Regarding the optimization of the ntree parameter, when the number of decision trees is 120, the model accuracy reaches its peak. Further increasing the number of decision trees does not significantly improve the performance of the model. Therefore, this embodiment will use 120 as the optimal value of ntree. Figure 6 shown.

[0037] S4 Crop structure classification: The fused image of the marked vegetation and non-vegetation areas obtained in step S2 is input into the trained Random Forest (RF) model to classify the areas marked as vegetation in the image and calculate the area. The crop structure type and area of ​​each crop in the vegetation area in the fused image of the marked vegetation and non-vegetation areas are determined. Figure 7 This is the final crop structure distribution map of the Hetao Irrigation District.

[0038] In order to verify the accuracy of the crop classification results of the random forest (RF) algorithm, a comparative analysis was conducted with the sown area proportions of each crop in the "2022 Bayannur Statistical Yearbook", and the results are shown in Table 5.

[0039] Table 5

[0040] Compared to the statistical yearbook data, the classification results show errors in area share of 3.45%, 2.66%, 0.11%, and 0.92% for corn, sunflower, wheat, and tomato, respectively. These errors may be due to the fact that the Hetao Irrigation District's boundaries do not completely coincide with the administrative boundaries of its affiliated banners and counties. This boundary deviation may contribute to the discrepancy in spatial extent between the statistical yearbook data and the classification results, thus affecting the accuracy of the comparison. However, since the cultivated land in Bayannur City is primarily concentrated in the Hetao Irrigation District, errors in boundary extent have a relatively small impact on the area calculation, further demonstrating the high reliability of the classification results.

Claims

1. A crop planting structure extraction method based on hierarchical extraction and multi-feature integration, characterized in that: The following steps are involved: S1 Data Acquisition and Preprocessing: Based on the GEE platform, two or more Sentinel-2A satellite imagery data of the target area are downloaded and de-clouded to obtain the de-clouded imagery data of the target area. S2 extracts vegetation areas: the NDVI frequency distribution histograms are calculated for the de-clouded image data obtained in step S1, the optimal thresholds of the de-clouded images are determined, the de-clouded images are fused using the OTSU algorithm with a maximum value strategy, and after fusion, the images are masked to obtain fused image data, the vegetation areas in the fused image data are determined, and image data with the vegetation areas marked is obtained; S3 classification model training: Sample data is obtained by combining GNSS RTK positioning with drone aerial photography. The sample data, spectral features, vegetation features, and texture features are used as model input information and fed into a random forest (RF) model for training. The trained random forest (RF) model is obtained. S4 Crop structure classification: The image data with the marked vegetation area obtained in step S2 is input into the trained Random Forest (RF) model to classify the vegetation areas marked in the image and calculate the area, thereby determining the crop structure type and area of ​​each crop in the vegetation area in the image data with the marked vegetation area.

2. The crop planting structure extraction method based on hierarchical extraction and multi-feature integration according to claim 1 is characterized in that: The de-clouding processing of the image data in step S1 includes calculating the values ​​of Bit10 and Bit11 of the band QA60, setting the image when the values ​​of Bit10 and Bit11 are both 0 as a cloud mask, using the cloud mask to remove cloud information in the image, obtaining a cloud-free image, and cropping the image according to the target area; if there is no image that meets the conditions, screening images with a cloud amount less than 10 and performing de-clouding processing, cropping according to the target area and then splicing, to obtain de-clouded image data of the target area for two periods.

3. The crop planting structure extraction method based on hierarchical extraction and multi-feature integration according to claim 1 is characterized in that: In step S3, bands B2, B3, B4, B5, B6, B7, B8, B8A, B11, and B12 of Sentinel-2A are selected as spectral features.

4. The crop planting structure extraction method based on hierarchical extraction and multi-feature integration according to claim 1 is characterized in that: In step S3, some band data from Sentinel-2A are selected to construct vegetation index as vegetation features, wherein the vegetation features include: Normalized Difference Vegetation Index (NDVI), which is calculated by the formula: , Normalized Difference Vegetation Index 1 (NDVI1), which is calculated as , Normalized Difference Vegetation Index 2 (NDVI2), the calculation formula is , Normalized Difference Vegetation Index 3 (NDVI3), the calculation formula is , the red edge normalized vegetation index 1 (NDre1), which is calculated as , the red edge normalized vegetation index 2 (NDre2), which is calculated as , Green Normalized Difference Vegetation Index (GNDVI), which is calculated as , Land Surface Water Index (LSWI), which is calculated as , bare soil index (BSI), which is calculated as , the modified chlorophyll absorption reflectance index (MCARI), which is calculated as follows: , Enhanced Vegetation Index (EVI), which is calculated as , Difference Vegetation Index (DVI), which is calculated as , Soil Adjusted Vegetation Index (SAVI), which is calculated as , Ratio Vegetation Index (RVI), which is calculated as , Ratio Vegetation Index 1 (RVI1), which is calculated as , Green Chlorophyll Vegetation Index (GCVI), which is calculated as .

5. The crop planting structure extraction method based on hierarchical extraction and multi-feature integration according to claim 1 is characterized in that: In step S3, texture features are extracted based on the B8 band in the Sentinel-2A data through the gray level co-occurrence matrix (GLCM) method, which includes calculating the second-order moment, mean sum, variance, contrast, difference, correlation, entropy and inverse moment of the B8 band.

6. The crop planting structure extraction method based on hierarchical extraction and multi-feature integration according to claim 1 is characterized in that: The number of decision trees of the random forest model in step S3 is set to 120.

Images