Soybean planting area extraction method based on multi-temporal Sentinel-2 data

Abstract

The present invention relates to a soybean planting area extraction method based on multi-phase Sentinel-2 data, comprising: obtaining Sentinel-2 image data and auxiliary data, and performing preprocessing; removing non-crop pixels in the image of the study area to obtain the overall vegetation distribution of the study area; generating a set of all features, fusing the data together, and performing masking; performing feature optimization, screening out the best feature subsets corresponding to each classifier, and selecting the best classifier; forming a soybean optimal extraction model by the obtained best classifier and the best feature subset corresponding to the classifier, and evaluating the soybean extraction effect of the soybean optimal extraction model, and obtaining the best soybean mapping effect in the study area. The present invention improves accuracy, reduces the probability of misclassification and omission; enriches spectral features, and also extracts some ground objects that are difficult to distinguish spectrally for use as auxiliary data; greatly reduces workload, reduces feature redundancy and noise, and improves work efficiency.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing image classification, and in particular to a soybean planting area extraction method based on multi-temporal Sentinel-2 data. Background Art

[0002] Soybeans are an important oilseed crop, accounting for a large proportion of the food structure and playing a vital role in the national security system. China is the country of origin of soybeans and was once the largest producer. It is now the largest importer and consumer of soybeans. According to the Ministry of Agriculture and Rural Affairs, my country's soybean production in 2020 was 19.6 million tons, an increase of 1.5 million tons year-on-year, reaching an increase of 8.3%. The sown area was 98,667 km 2 , an increase of 5,500 km 2 , reaching a growth rate of 5.9%. At the same time, starting in 2019, the Ministry of Agriculture and Rural Affairs implemented the Soybean Revitalization Plan, which is a rare opportunity for my country's soybean industry. Therefore, accurately grasping the spatial distribution information of soybeans is of great significance for regulating the soybean market, strengthening soybean production management, and ensuring sustainable development.

[0003] Currently, traditional agricultural survey methods have significant limitations for estimating soybean acreage. They are costly, time-consuming, and susceptible to human interference. Unmanned aerial vehicle (UAV) remote sensing technology is suitable for studying rapidly changing and seasonal agricultural production and offers flexibility. However, it is difficult to collect data in remote areas, and the data collected is limited, making it unsuitable for large-scale remote sensing monitoring. Satellite remote sensing technology, on the other hand, is suitable for large-scale monitoring and offers advantages such as efficiency, convenience, and low cost. Furthermore, most soybean research has focused on the United States, Argentina, Brazil, and Northeast China, characterized by simple planting structures, a single crop variety, and large-scale cultivation. However, research on soybeans in Anhui Province, with its mountainous terrain, unpredictable weather, fragmented farmland, and low level of mechanization, is rare. Summary of the Invention

[0004] The purpose of the present invention is to provide a soybean planting area extraction method based on multi-temporal Sentinel-2 data, which can effectively and accurately classify the ground objects in the study area based on remote sensing data, effectively improve the classification accuracy by integrating the features of multiple temporal phases, and greatly reduce the workload.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a soybean planting area extraction method based on multi-temporal Sentinel-2 data, the method comprising the following steps in sequence:

[0006] (1) Acquire Sentinel-2 image data and auxiliary data, preprocess the Sentinel-2 image data and auxiliary data to obtain images of the study area;

[0007] (2) Using the decision tree method, non-crop pixels in the image of the study area were removed to obtain the overall distribution of vegetation in the study area;

[0008] (3) Based on the Sentinel-2 image data preprocessed in step (1), a set of all features is generated, and the vegetation index, texture features, microwave features, Sentinel-1IW GRD level data and the Sentinel-2 image data preprocessed in step (1) are fused together, and the fused data is masked using the overall distribution of vegetation in the study area to obtain a set of all masked features;

[0009] (4) Perform feature optimization on the set of all masked features to screen out the best feature subsets corresponding to each classifier, and select the best classifier based on the classification results;

[0010] (5) The best classifier obtained and the best feature subset corresponding to the classifier are used to form the best soybean extraction model, and the soybean extraction effect of the best soybean extraction model is evaluated, and the best soybean mapping effect in the study area is obtained.

[0011] In step (1), the Sentinel-2 image data is obtained by downloading 6 L2A-level Sentinel-2 image data from the European Space Agency Copernicus data, which are respectively on August 18, 2019, August 28, 2019 and September 7, 2019; the auxiliary data include Sentinel-1IW GRD-level data, Planet verification data, vector data, surface cover map and statistical yearbook; the preprocessing includes: for the Sentinel-2 image data, all its bands are resampled to 10m using bilinear interpolation method, and the obtained 10 bands with a resolution of 10m are exported to ENVI format, and band synthesis, mosaicking, cropping and reflectance calculation are performed using ENVI software to obtain the image of the study area; for the Sentinel-1IW GRD-level data in the auxiliary data, the Graph Graph in SNAP software is used to generate the image of the study area. The Builder process tool performs thermal noise removal, orbit correction, radiometric calibration, speckle filtering, terrain correction, and decibel conversion, and then mosaics and crops the results to obtain images of the study area.

[0012] The step (2) specifically refers to: constructing a decision tree using the normalized building index (NDBI), the improved normalized water index (MNDWI), the B8 band of the Sentinel-2 image data, and the FROM-GLC10 file, and sequentially removing buildings, roads, water bodies, and bare soil from the image of the study area through the decision tree; the expression of the normalized building index (NDBI) is:

[0013] NDBI=(B11-B8) / (B11+B8)

[0014] Where, B8 is the reflectivity value of the near-infrared band, and B11 is the reflectivity value of the short-wave infrared band;

[0015] The improved normalized water index MNDWI is an improvement on the wavelength combination mode based on the normalized difference index NDWI, and its expression is:

[0016] MNDWI=(B3-B11) / (B3+B11)

[0017] Where, B3 is the reflectance value of the green band, and B11 is the reflectance value of the short-wave infrared band;

[0018] The B8 band is the 10m near-infrared band that comes with Sentinel-2 image data;

[0019] The elimination of non-crop pixels in the image of the study area refers to masking the non-crop pixels in the image of the study area through a decision tree to obtain the overall distribution of vegetation in the study area.

[0020] In step (3), the set of all features includes vegetation index, texture feature, microwave feature, Sentinel-1 IWGRD level data and Sentinel-2 image data after preprocessing in step (1); the vegetation index is a linear or nonlinear combination between bands, reflecting the active radiation value of vegetation and the richness of vegetation, and the vegetation indices of EVI, MTCI, NDVIre2, REP, and SAVI of the three phases of Sentinel-2 image data, the vegetation index obtained by fusion of August 18 and August 28, and the vegetation index obtained by fusion of August 28 and August 28 are selected respectively. The vegetation indices obtained after mutual fusion on September 7th were 25 in total. The texture features were extracted by using the gray-level co-occurrence matrix in ENVI software and a 3×3 window to extract 8 texture features of the three phases of Sentinel-2 image data, namely mean, variance, contrast, second-order moment, homogeneity, information entropy, correlation, heterogeneity, subtraction of the texture features corresponding to August 18 and August 28, and subtraction of the texture features corresponding to August 28 and September 7. A total of 40 texture features were obtained, and 8 texture features were randomly screened and retained. Through Sentinel-1IW GRD-level data obtains data in two polarization modes, VV and VH. Based on the two data, algebraic operations are performed to obtain VH+VV, VH-VV, VH×VV, and VH / VV of the three phases of Sentinel-2 image data, respectively, for a total of 30 microwave features. Then, the VH, VV, and VH+VV features of the three phases of Sentinel-2 image data are randomly screened, resulting in a total of 9 microwave features.

[0021] The step (4) specifically includes: using ENVI software to mark the ground feature samples in the study area, selecting the features of the three phases of Sentinel-2 image data, and the features generated by the fusion of the phases, randomly screening 85 features, and using the ReliefF algorithm to perform feature optimization on these 85 features, thereby obtaining a ranking diagram of all feature weights. The formula of the ReliefF algorithm is as follows:

[0022]

[0023] Where H(x) is the nearest sample point of the same type as sample x; M(x) is the nearest sample point of a different type from sample x; θ is the hypothesis interval, that is, the maximum distance the decision surface can move while keeping the sample classification unchanged;

[0024] To screen the best feature subset, the sequential forward selection method is used. In a coupled method with the classifier, features are first added to the classifier in the order of feature weights. Each time a feature is added, an overall classification accuracy OA is given until all 85 features are input and the overall classification accuracy curve is drawn. Judging by the overall classification accuracy curve, if the overall classification accuracy increases with the addition of features, the feature is retained. When the overall classification accuracy decreases with the addition of features, the feature is removed. When the accuracy reaches the maximum value, all subsequent features are removed, and the best feature subset is finally screened. The soybean classification result map is finally generated by using three classifiers: random forest, BP neural network, and support vector machine (SVM) and their corresponding best feature subsets. Planet images are used to classify soybeans in each sample plot. As the real soybean distribution map in the sample plot, the classification results of the three classifiers are verified. The best classifier is judged based on the kappa coefficient. The higher the kappa coefficient, the higher the classification result.

[0025] The calculation formula of the kappa coefficient is as follows:

[0026]

[0027] Where N is the total number of pixels, m is the number of categories, and x kk is the number of pixels on the diagonal of the confusion matrix, x k+ and x +k are the total number of pixels in the k-th row and k-th column respectively.

[0028] In step (5), the evaluation of the soybean extraction effect of the optimal soybean extraction model specifically refers to the evaluation by the following two methods:

[0029] (6a) Input the Sentinel-2 image data after preprocessing in step (1) into the optimal classifier to obtain the classification results. Use the 3m resolution planet image to verify the classification results and generate the confusion matrix and kappa coefficient;

[0030] (6b) Input the best feature subset into the best classifier to obtain the classification results. Use 3m resolution planet images to verify the classification results and generate confusion matrix and kappa coefficient;

[0031] The confusion matrix includes overall accuracy OA, mapping accuracy PA and user accuracy UA, wherein the overall accuracy OA is the ratio of the number of correctly classified pixels to the total number of categories; the mapping accuracy PA is the probability that the classifier classifies all pixels into N categories assuming that the pixels are actually N categories; the user accuracy UA is the probability that the classifier classifies the pixels into N categories when the pixel is judged to be N categories.

[0032] The kappa coefficient obtained by steps (6a) and (6b) is used to test the classification effect of the optimal soybean extraction model. The higher the kappa coefficient, the better the classification effect of the optimal soybean extraction model. The optimal soybean mapping effect in the study area is obtained by the optimal soybean extraction model.

[0033] It can be seen from the above technical solution that the beneficial effects of the present invention are: first, a decision tree is constructed by normalizing the water body index, the normalized building index and the FROM-GLC10 data, and the non-vegetation pixels in the study area are masked, thereby improving the accuracy and reducing the probability of misclassification and omission; second, the present invention abandons the original method of calculating the optimal phase by JM distance, and adopts a method of combining multiple phases to extract the favorable features of other phases, which not only enriches the spectral features, but also extracts some ground objects that are difficult to distinguish spectrally for use as auxiliary data; third, by establishing the optimal soybean extraction model and screening the optimal feature subset, the workload is greatly reduced, the feature redundancy and noise are reduced, and the work efficiency is improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 is a flow chart of the method of the present invention;

[0035] Figure 2 Schematic diagram of using decision tree to mask non-vegetation pixels;

[0036] Figure 3 、 4 They are all ranking diagrams of all feature weights;

[0037] Figure 5 、 6 , 7 are line graphs between feature dimension and classification accuracy. DETAILED DESCRIPTION

[0038] like Figure 1 As shown, a soybean planting area extraction method based on multi-temporal Sentinel-2 data includes the following steps in sequence:

[0039] (1) Acquire Sentinel-2 image data and auxiliary data, preprocess the Sentinel-2 image data and auxiliary data to obtain images of the study area;

[0040] (2) Using the decision tree method, non-crop pixels in the image of the study area were removed to obtain the overall distribution of vegetation in the study area;

[0041] (3) Based on the Sentinel-2 image data preprocessed in step (1), a set of all features is generated, and the vegetation index, texture features, microwave features, Sentinel-1IW GRD level data and the Sentinel-2 image data preprocessed in step (1) are fused together, and the fused data is masked using the overall distribution of vegetation in the study area to obtain a set of all masked features;

[0042] (4) Perform feature optimization on the set of all masked features to screen out the best feature subsets corresponding to each classifier, and select the best classifier based on the classification results;

[0043] (5) The best classifier obtained and the best feature subset corresponding to the classifier are used to form the best soybean extraction model, and the soybean extraction effect of the best soybean extraction model is evaluated, and the best soybean mapping effect in the study area is obtained.

[0044] In step (1), the Sentinel-2 image data is obtained by downloading 6 L2A-level Sentinel-2 image data from the European Space Agency Copernicus data, which are respectively on August 18, 2019, August 28, 2019 and September 7, 2019; the auxiliary data include Sentinel-1IW GRD-level data, Planet verification data, vector data, surface cover map and statistical yearbook; the preprocessing includes: for the Sentinel-2 image data, all its bands are resampled to 10m using bilinear interpolation method, and the obtained 10 bands with a resolution of 10m are exported to ENVI format, and band synthesis, mosaicking, cropping and reflectance calculation are performed using ENVI software to obtain the image of the study area; for the Sentinel-1IW GRD-level data in the auxiliary data, the Graph Graph in SNAP software is used to generate the image of the study area. The Builder process tool performs thermal noise removal, orbit correction, radiometric calibration, speckle filtering, terrain correction, and decibel conversion, and then mosaics and crops the results to obtain images of the study area. Specific Sentinel-2 band information is shown in Table 1 below:

[0045] Table 1 Sentinel-2 band information

[0046]

[0047] like Figure 2As shown, step (2) specifically refers to: constructing a decision tree through the normalized building index NDBI, the improved normalized water index MNDWI, the B8 band of Sentinel-2 image data and the FROM-GLC10 file, and sequentially eliminating buildings, roads, water bodies and bare soil in the image of the study area through the decision tree; the expression of the normalized building index NDBI is:

[0048] NDBI=(B11-B8) / (B11+B8)

[0049] Where, B8 is the reflectivity value of the near-infrared band, and B11 is the reflectivity value of the short-wave infrared band;

[0050] The improved normalized water index MNDWI is an improvement on the wavelength combination mode based on the normalized difference index NDWI, and its expression is:

[0051] MNDWI=(B3-B11) / (B3+B11)

[0052] Where, B3 is the reflectance value of the green band, and B11 is the reflectance value of the short-wave infrared band;

[0053] The B8 band is the 10m near-infrared band that comes with Sentinel-2 image data;

[0054] The elimination of non-crop pixels in the image of the study area refers to masking the non-crop pixels in the image of the study area through a decision tree to obtain the overall distribution of vegetation in the study area.

[0055] In step (3), the set of all features includes vegetation index, texture feature, microwave feature, Sentinel-1IW GRD level data and Sentinel-2 image data after preprocessing in step (1); the vegetation index is a linear or nonlinear combination between bands, reflecting the active radiation value of vegetation and the richness of vegetation, and the vegetation indices of EVI, MTCI, NDVIre2, REP, and SAVI of the three phases of Sentinel-2 image data, the vegetation index obtained by fusion of August 18 and August 28, and the vegetation index obtained by fusion of August 28 and September 7 are selected respectively, and a total of 25 vegetation indices are obtained; the specific vegetation index data are shown in Table 2 below:

[0056] Table 2 Vegetation index used in the present invention

[0057]

[0058] The texture features are extracted by using the gray level co-occurrence matrix in the ENVI software and a 3×3 window, respectively, for the three phases of the Sentinel-2 image data, namely, mean, variance, contrast, second-order moment, homogeneity, information entropy, correlation, heterogeneity, subtraction of the texture features corresponding to August 18 and August 28, and subtraction of the texture features corresponding to August 28 and September 7, obtaining a total of 40 texture features, and randomly screening and retaining 8 texture features; using the Sentinel-1IW GRD-level data to obtain data of two polarization modes, VV and VH, and performing algebraic operations on the basis of the two data to obtain VH+VV, VH-VV, VH×VV, and VH / VV of the three phases of the Sentinel-2 image data, respectively, obtaining a total of 30 microwave features, and then randomly screening the VH, VV, and VH+VV features of the three phases of the Sentinel-2 image data to obtain a total of 9 microwave features. In the present invention, the random screening can be performed by visual interpretation or by optical instrument screening to screen out features with clear features and less noise and eliminate blurred images to ensure better screening effect.

[0059] The step (4) specifically refers to: using ENVI software, marking the ground feature samples in the study area, selecting the features of the three phases of Sentinel-2 image data, and the features generated after the fusion of each phase, randomly screening 85 features, and using the ReliefF algorithm to perform feature optimization on these 85 features, thereby obtaining a ranking diagram of all feature weights, as shown in the figure below: Figure 3 、 Figure 4 As shown; the formula of the ReliefF algorithm is as follows:

[0060]

[0061] Where H)x) is the nearest sample point of the same type as sample x; M)x) is the nearest sample point of a different type from sample x; θ is the hypothesis interval, that is, the maximum distance the decision surface can move while keeping the sample classification unchanged; Tables 3 and 4 contain the numbers of all features:

[0062] Table 3 Characteristic variable numbers

[0063]

[0064]

[0065] Table 4 Characteristic variable numbers

[0066]

[0067] Bands B8A, B8, B7, and B6 rank highly, while vegetation indices like EVI and SAVI also have high weights, highlighting the importance of the red edge and shortwave infrared bands for soybean identification. The texture feature MEAN ranks highly in the overall feature ranking, indicating that compared to other texture features, it is most beneficial for soybean extraction. However, while microwave data is unaffected by weather and sunlight, it ranks low overall in the ranking, contributing little to soybean extraction. Within the overall feature ranking, features from August 18th rank highly, followed closely by most features from August 28th. Features from other time phases generally have low weights, indicating that, among the five time phases, August 18th, early in the soybean pod-setting period, is the most favorable for soybean classification. While the shortwave infrared and red edge bands rank highly, the importance of features from other time phases cannot be underestimated.

[0068] To screen the best feature subset, the sequential forward selection method is used. In a coupled method with the classifier, features are first added to the classifier in the order of feature weights. An overall classification accuracy OA is given each time a feature is added until all 85 features are input and an overall classification accuracy curve is drawn. Judging by the overall classification accuracy curve, if the overall classification accuracy increases with the addition of features, the feature is retained. If the overall classification accuracy decreases with the addition of features, the feature is removed. When the accuracy reaches the maximum value, all subsequent features are removed, and finally the best feature subset is screened.

[0069] As shown in Table 5 below, the ReliefF-RF model discards the feature subsets ranked 8th, 10th, 12th, 14th-15th, and 19th-15th, and ultimately retains the top 26 13 feature subsets, such as Figure 5 As shown in the figure, the ReliefF-SVM model discards the features numbered 8, 10, 12, -16, 18-25, and finally retains the first 26 feature subsets, as shown in the figure. Figure 6 As shown in Figure 2, the ReliefF-BPNN model discards the features ranked 8, 10-15, 17-19, and 21-23, and ultimately retains the first 24 feature subsets, as shown in Figure 2. Figure 7 As shown in Table 5, the best feature subsets of the three models are listed. Combining the best feature subsets with each classifier, we can get the ReliefF-RF model, ReliefF-SVM model, and ReliefF-BPNN model. Figure 5 、 6 ,The triangle in 7 indicates that the classification accuracy reaches the highest,at this time.

[0070] Table 5 Best feature subsets for different models

[0071]

[0072] The soybean classification result map was finally generated by using three classifiers: random forest, BP neural network, and support vector machine (SVM) and their corresponding optimal feature subsets. The planet image was used to classify the soybeans in each sample plot. The actual distribution map of soybeans in the sample plot was used to verify the classification results of the three classifiers. The best classifier was determined based on the kappa coefficient. The higher the kappa coefficient, the better the classification result.

[0073] The calculation formula of the kappa coefficient is as follows:

[0074]

[0075] Where N is the total number of pixels, m is the number of categories, and x kk is the number of pixels on the diagonal of the confusion matrix, x k+ and x +k are the total number of pixels in the k-th row and k-th column respectively.

[0076] In step (5), the evaluation of the soybean extraction effect of the optimal soybean extraction model specifically refers to the evaluation by the following two methods:

[0077] (6a) Input the Sentinel-2 image data after preprocessing in step (1) into the optimal classifier to obtain the classification results. Use the 3m resolution planet image to verify the classification results and generate the confusion matrix and kappa coefficient;

[0078] (6b) Input the best feature subset into the best classifier to obtain the classification results. Use 3m resolution planet images to verify the classification results and generate confusion matrix and kappa coefficient;

[0079] The confusion matrix includes overall accuracy OA, mapping accuracy PA and user accuracy UA, wherein the overall accuracy OA is the ratio of the number of correctly classified pixels to the total number of categories; the mapping accuracy PA is the probability that the classifier classifies all pixels into N categories assuming that the pixels are actually N categories; the user accuracy UA is the probability that the classifier classifies the pixels into N categories when the pixel is judged to be N categories.

[0080] The kappa coefficient obtained by steps (6a) and (6b) is used to test the classification effect of the optimal soybean extraction model. The higher the kappa coefficient, the better the classification effect of the optimal soybean extraction model. The optimal soybean mapping effect in the study area is obtained by the optimal soybean extraction model.

[0081] The optimal soybean extraction models were evaluated. The optimal feature subsets for each model were obtained from Table 5. Each optimized optimal feature subset was classified and validated using planet imagery. The cartographic accuracy, user precision, overall precision, and Kappa coefficients of the classification results for each model were calculated, as shown in Table 6. The results show that the ReliefF-RF model's Kappa coefficient was higher than the other two models in all seven plots, and its overall precision was slightly higher than the other models, with this performance being most pronounced in plot 7. The ReliefF-SVM model's cartographic accuracy was generally higher than the other two plots, but its user precision was lower, indicating that it was more likely to classify other landforms as soybeans than the other two models. The ReliefF-BPNN model had the lowest kappa coefficient of the three models, with average cartographic accuracy and user precision. In the classification results using only the reflectance of the ten raw bands of Sentinel-2 without optimization, its kappa coefficient was the lowest among the three models. At the same time, it was found that sample 5 had the lowest kappa coefficient among all the sample plots, and the overall mapping accuracy was also low. Through observation, it was found that the spatial distribution of soybeans in sample 5 was fragmented, and its area accounted for the smallest proportion, making it easy to classify soybeans into other land features.

[0082] Table 6

[0083]

[0084]

[0085] In summary, the present invention abandons the original method of calculating the optimal phase through JM distance, and adopts a method combining multiple phases to extract the favorable features of other phases, explore the powerful features of other phases, and greatly enrich the spectral features; using other auxiliary data, it also extracts some ground objects that are difficult to distinguish spectrally, greatly improving the classification accuracy; through the normalized building index NDBI, the improved normalized water index MNDWI, and the FROM-GLC10 file to construct a decision tree, the non-vegetation pixels in the study area are masked, while reducing the sample requirement, thereby improving the accuracy and reducing the probability of misclassification and omission; through the establishment of the optimal soybean extraction model and the screening of the optimal feature subset, the workload is greatly reduced, the feature redundancy and noise are reduced, and the work efficiency is improved.

Claims

1. A soybean planting area extraction method based on multi-temporal Sentinel-2 data, characterized by: The method comprises the following steps in sequence: (1) Acquire Sentinel-2 image data and auxiliary data, preprocess the Sentinel-2 image data and auxiliary data to obtain images of the study area; (2) Using the decision tree method, non-crop pixels in the image of the study area were removed to obtain the overall distribution of vegetation in the study area; (3) Based on the Sentinel-2 image data preprocessed in step (1), a set of all features is generated, vegetation index, texture features, microwave features, Sentinel-1 IW GRD level data and Sentinel-2 image data preprocessed in step (1) are fused together, and the fused data is masked using the overall distribution of vegetation in the study area to obtain a set of all masked features; (4) Perform feature optimization on the set of all masked features to screen out the best feature subsets corresponding to each classifier, and select the best classifier based on the classification results; (5) The best soybean extraction model is formed by the obtained best classifier and the best feature subset corresponding to the classifier, and the soybean extraction effect of the best soybean extraction model is evaluated, and the best soybean mapping effect in the study area is obtained; The step (4) specifically includes: using ENVI software to mark the ground feature samples in the study area, selecting the features of the three phases of Sentinel-2 image data, and the features generated by the fusion of the phases, randomly screening 85 features, and using the ReliefF algorithm to perform feature optimization on these 85 features, thereby obtaining a ranking diagram of all feature weights. The formula of the ReliefF algorithm is as follows: Where H(x) is the nearest sample point of the same type as sample x; M(x) is the nearest sample point of a different type from sample x; θ is the hypothesis interval, that is, the maximum distance the decision surface can move while keeping the sample classification unchanged; To screen the best feature subset, the sequential forward selection method is used. In a coupled method with the classifier, features are first added to the classifier in the order of feature weights. Each time a feature is added, an overall classification accuracy OA is given until all 85 features are input and the overall classification accuracy curve is drawn. Judging by the overall classification accuracy curve, if the overall classification accuracy increases with the addition of features, the feature is retained. When the overall classification accuracy decreases with the addition of features, the feature is removed. When the accuracy reaches the maximum value, all subsequent features are removed, and the best feature subset is finally screened. The soybean classification result map is finally generated by using three classifiers: random forest, BP neural network, and support vector machine (SVM) and their corresponding best feature subsets. Planet images are used to classify soybeans in each sample plot. As the real soybean distribution map in the sample plot, the classification results of the three classifiers are verified. The best classifier is judged based on the kappa coefficient. The higher the kappa coefficient, the higher the classification result. The calculation formula of the kappa coefficient is as follows: Where N is the total number of pixels, m is the number of categories, and x kk is the number of pixels on the diagonal of the confusion matrix, x k+ and x +k are the total number of pixels in the k-th row and k-th column respectively.

2. The soybean planting area extraction method based on multi-temporal Sentinel-2 data according to claim 1, characterized in that: In step (1), the Sentinel-2 image data is obtained by downloading 6 L2A-level Sentinel-2 image data from the European Space Agency Copernicus data, which are respectively on August 18, 2019, August 28, 2019 and September 7, 2019; the auxiliary data include Sentinel-1 IW GRD-level data, Planet verification data, vector data, surface cover map and statistical yearbook; the preprocessing includes: for the Sentinel-2 image data, all its bands are resampled to 10m using bilinear interpolation method, and the obtained 10 bands with a resolution of 10m are exported to ENVI format, and the bands are synthesized, mosaicked, cropped and the reflectivity is calculated using ENVI software to obtain the image of the study area; for the Sentinel-1 IWGRD-level data in the auxiliary data, the Graph Graph in SNAP software is used to generate the image of the study area. The Builder process tool performs thermal noise removal, orbit correction, radiometric calibration, speckle filtering, terrain correction, and decibel conversion, and then mosaics and crops the results to obtain images of the study area.

3. The soybean planting area extraction method based on multi-temporal Sentinel-2 data according to claim 1, characterized in that: The step (2) specifically refers to: constructing a decision tree using the normalized building index (NDBI), the improved normalized water index (MNDWI), the B8 band of the Sentinel-2 image data, and the FROM-GLC10 file, and sequentially removing buildings, roads, water bodies, and bare soil from the image of the study area through the decision tree; the expression of the normalized building index (NDBI) is: NDBI=(B11-B8) / (B11+B8) Where, B8 is the reflectivity value of the near-infrared band, and B11 is the reflectivity value of the short-wave infrared band; The improved normalized water index MNDWI is an improvement on the wavelength combination mode based on the normalized difference index NDWI, and its expression is: MNDWI=(B3-B11) / (B3+B11) Where, B3 is the reflectance value of the green band, and B11 is the reflectance value of the short-wave infrared band; The B8 band is the 10m near-infrared band that comes with Sentinel-2 image data; The elimination of non-crop pixels in the image of the study area refers to masking the non-crop pixels in the image of the study area through a decision tree to obtain the overall distribution of vegetation in the study area.

4. The soybean planting area extraction method based on multi-temporal Sentinel-2 data according to claim 1, characterized in that: In step (3), the set of all features includes vegetation index, texture feature, microwave feature, Sentinel-1 IWGRD-level data and Sentinel-2 image data after preprocessing in step (1); the vegetation index is a linear or nonlinear combination between bands, reflecting the active radiation value of vegetation and the richness of vegetation, and the vegetation indices of EVI, MTCI, NDVIre2, REP, and SAVI of the three phases of Sentinel-2 image data, the vegetation index obtained by mutual fusion of August 18 and August 28, and the vegetation index obtained by mutual fusion of August 28 and September 7 are selected respectively, for a total of 25 vegetation indices were obtained; the texture features were extracted by using the gray level co-occurrence matrix in the ENVI software and a 3×3 window, respectively, for the three phases of the Sentinel-2 image data, namely, mean, variance, contrast, second-order moment, homogeneity, information entropy, correlation, heterogeneity, subtraction of the texture features corresponding to August 18 and August 28, and subtraction of the texture features corresponding to August 28 and September 7, a total of 40 texture features were obtained, and 8 texture features were randomly screened and retained; VV and VH polarization data were obtained by Sentinel-1 IW GRD-level data, and algebraic operations were performed on the basis of the two data to obtain VH+VV, VH-VV, VH×VV, and VH / VV of the three phases of the Sentinel-2 image data, respectively, for a total of 30 microwave features, and then the VH, VV, and VH+VV features of the three phases of the Sentinel-2 image data were randomly screened to obtain a total of 9 microwave features.

5. The soybean planting area extraction method based on multi-temporal Sentinel-2 data according to claim 1, characterized in that: In step (5), the evaluation of the soybean extraction effect of the optimal soybean extraction model specifically refers to the evaluation by the following two methods: (6a) Input the Sentinel-2 image data after preprocessing in step (1) into the optimal classifier to obtain the classification results. Use the 3m resolution planet image to verify the classification results and generate the confusion matrix and kappa coefficient; (6b) Input the best feature subset into the best classifier to obtain the classification results. Use 3m resolution planet images to verify the classification results and generate confusion matrix and kappa coefficient; The confusion matrix includes overall accuracy OA, mapping accuracy PA and user accuracy UA, wherein the overall accuracy OA is the ratio of the number of correctly classified pixels to the total number of categories; the mapping accuracy PA is the probability that the classifier classifies all pixels into N categories assuming that the pixels are actually N categories; the user accuracy UA is the probability that the classifier classifies the pixels into N categories when the pixel is judged to be N categories. The kappa coefficient obtained by steps (6a) and (6b) is used to test the classification effect of the optimal soybean extraction model. The higher the kappa coefficient, the better the classification effect of the optimal soybean extraction model. The optimal soybean mapping effect in the study area is obtained by the optimal soybean extraction model.

Images