A method for extracting cultivated land based on temporal multispectral image features

Abstract

The present invention discloses a method for extracting cultivated land based on time-series multispectral image features. The method comprises the following steps: obtaining multiple periods of Sentinel-2 image data for preprocessing, and producing cultivated land sample labels for the image data; calculating index features, spatial grayscale features, and texture features for the image data, performing dimensionality reduction, and calculating time-series features; screening features using a feature classifier based on the index features, spatial grayscale features, texture features, and time-series features, establishing a cultivated land classification model, and training the model; and extracting cultivated land patches from the Sentinel-2 image data. The present invention extracts optimal spatial, spectral, textural, and temporal features of remote sensing image data for cultivated land, explores the relationships between multidimensional feature spaces such as spectra, vegetation indices, textures, and statistical features in the time dimension related to cultivated land extraction, and provides an effective method for extracting spatiotemporal spectral features of cultivated land. The method has feature generalization and migration capabilities and good robustness.

Description

Technical Field

[0001] The present invention relates to the technical field of farmland extraction based on multispectral images, and in particular to a farmland extraction method based on time-series multispectral image features. Background Art

[0002] The extraction and spatial distribution of cultivated land have a significant impact on crop yield assessment and food safety. Remote sensing technology is considered an effective tool for providing accurate and objective cultivated land information reproducibly across time and space. Currently, a variety of satellite remote sensing data sources are available for extracting cultivated land information, such as MODIS, Landsat, and Sentinel-2. Remote sensing cultivated land detection primarily focuses on surveying the quantity and spatial distribution of cultivated land. Technical methods have evolved from visual interpretation to semi-automated methods based on statistical models, and are gradually shifting to data-driven classification methods primarily based on machine learning and deep learning. Unlike general land surface types such as buildings, roads, and impervious surfaces, the concept of cultivated land encompasses a wider range of land features. Different crops may be grown on cultivated land at the same time, resulting in significant image differences and large inter-class differences. In addition, the significant interannual changes in cultivated land lead to low classification accuracy in local areas of existing methods. The degree of manual participation in the implementation process is still high. The selection of classification features and labels as model inputs directly affects the accuracy of cultivated land extraction. There is currently no effective fixed feature combination method, which makes feature generalization and migration difficult and lacks robustness.

[0003] For example, invention application number 202410084966.1 discloses a method for extracting cultivated land information from high-resolution object-oriented images. This method comprehensively considers the spatial resolution of the image, the complexity of the distribution of land features, and the degree of regional fragmentation. Based on the shape characteristics of roads and ridges, it extracts road and ridge areas, which not only distinguish the boundaries of individual cultivated land but also obtains more detailed cultivated land areas. However, this solution also has the problem that different crops may be planted on the cultivated land at the same time, resulting in significant image differences and poor classification of cultivated land.

[0004] Therefore, in order to solve the problem of misclassification and omission of cultivated land categories caused by the inability of single-phase data to reflect the differences in cultivated land such as surface vegetation, the optimal remote sensing image data space, spectrum, texture, and time series features for cultivated land extraction are objectively determined. The present invention mines the relationship between multidimensional feature spaces such as spectrum, vegetation index, texture, and statistical features in the time dimension related to cultivated land extraction, and publishes an effective extraction method for spatiotemporal spectral features of cultivated land based on Sentinel-2, which can achieve the generalization and migration capabilities of cultivated land features and has good robustness. Summary of the Invention

[0005] In response to the above-mentioned problems, the purpose of the present invention is to provide a method for extracting cultivated land based on time-series multispectral image features, which can realize large-scale, heterogeneous data cultivated land extraction from satellite remote sensing images, improve the accuracy of cultivated land classification, determine the optimal combination of spatiotemporal spectral features, and have regional migration capabilities.

[0006] The embodiment of the present invention provides a method for extracting cultivated land based on time-series multispectral image features, comprising the steps of:

[0007] S1. Acquire multiple phases of Sentinel-2 image data, pre-process the image data, mark the cultivated land areas, and create cultivated land sample labels for the image data;

[0008] S2. Calculate index features, spatial grayscale features, and texture features for the preprocessed image data, perform dimensionality reduction, and calculate the time series features after dimensionality reduction;

[0009] S3. Based on index features, spatial grayscale features, texture features, and temporal features, a feature classifier is used to select the features that have the greatest impact on the classification results, and a cultivated land classification model is established and trained.

[0010] S4. Extract the cultivated land patches from the Sentinel-2 image data based on the trained cultivated land classification model and save the results.

[0011] Furthermore, the image data is pre-processed in S1, including the following steps:

[0012] S11. Perform multi-period image spatial registration, overlap region cropping, and farmland labeling on Sentinel-2 remote sensing images.

[0013] S12, resampling the cropped remote sensing image data of each band, reading metadata of the central wavelength of each band, and performing multi-band rearrangement;

[0014] S13. Calculate multi-band index features of remote sensing images based on true color band and multi-band combination data.

[0015] Furthermore, the multi-band index features include:

[0016] NDVI = (NIR-R) / (NIR+R)

[0017] NDWI=(G-NIR) / (G+NIR)

[0018] EVI=2.5*(NIR-R) / (NIR+6R-7.5B+1)

[0019] MSAVI=(2NIR+1-sqrt((2NIR-1) 2 -8(NIR-R)) / 2

[0020] NDBI = (MIR - NIR) / (MIR + NIR)

[0021] LSWI=(NIR-SWIR) / (NIR-SWIR)

[0022] NDTI=(SWIR1-SWIR2) / (SWIR1+SWIR2)

[0023] REP=(705+35*(0.5*(RE3+R)-RE1) / (RE2-RE1)

[0024] Among them, B, G, and R are true color bands, which are the B2, B3, and B4 bands of Sentinel-2 remote sensing images respectively. NIR, SWIR1, SWIR2, RE1, RE2, and RE3 are the B2, B3, B4, B8A, B11, B12, B5, B6, and B7 bands of Sentinel-2 remote sensing images respectively. MIR is the B11 band of Sentinel-2 remote sensing images.

[0025] Furthermore, the spatial grayscale features and texture features include:

[0026] Homogeneity:

[0027] Contrast ratio:

[0028] Mean:

[0029] variance:

[0030] Dependencies:

[0031] Dissimilarity:

[0032] Angular second moment:

[0033] The distance in each feature calculation is set to 1, the grayscale level is 16, and each result is the mean value in the four directions of 0°, 45°, 90°, and 135°, which is expressed in a unified formula:

[0034] p(i,j,d,θ)={[(i,j),(x+dx,y+dy)]|f(x,y)=i,f(x+dx,y+dx=j)

[0035] Where d is the specified distance between two gray levels and θ is the direction.

[0036] Furthermore, the time series characteristic formula is expressed as:

[0037]

[0038] Wherein, T is the total number of image time series, t is the current phase of the image, and F is the calculated feature value.

[0039] Furthermore, the exponential features, spatial grayscale features, and texture features are calculated and dimensionally reduced in S2, and the time series features after dimension reduction are calculated, including the following steps:

[0040] S21. Calculate the index features of the image data. Based on the Pearson correlation coefficient, consider the index features with a correlation higher than 97% as the same feature, and reduce the dimensionality of the index features.

[0041] S22. Calculate texture features based on the exponential features after dimensionality reduction. Based on the Pearson correlation coefficient, consider texture features with a correlation greater than 97% as the same feature and reduce the dimensionality of the texture features. Simultaneously, grayscale the true color image and calculate the spatial grayscale features.

[0042] S23, superimposing the single features after correlation screening to calculate the time series features, and using the coefficient of variation less than 15% as the threshold to reduce the dimension of the time series features;

[0043] S24. Obtain features after dimensionality reduction, including spatial grayscale features, texture features, exponential features, and time series features.

[0044] Furthermore, the feature dimension reduction index includes the Pearson correlation coefficient and the coefficient of variation, wherein,

[0045] Pearson correlation coefficient, the formula is expressed as:

[0046]

[0047] Among them, r is the correlation coefficient, X and Y are pixel values.

[0048] The coefficient of variation is expressed as:

[0049]

[0050] CV=σ / μ*100%

[0051] μ is the characteristic mean, σ is the characteristic standard deviation, and CV is the coefficient of variation.

[0052] Furthermore, the cultivated land classification model is established and trained in S3, including:

[0053] LightGBM is selected as the classifier, and IOU and F1-score are selected as the accuracy evaluation indicators. The accuracy score is obtained based on the accuracy evaluation indicators. The model is iteratively trained according to the accuracy score until the average absolute change of the accuracy score in every three rounds is less than 3%;

[0054] Based on the trained model, the contribution of each feature to the classification accuracy is counted, and feature screening is performed using gain as a measurement indicator. High-gain features are retained, and the model is re-trained iteratively until the accuracy score of each five rounds of model training iterations increases by less than 1%. The iteration is stopped, and the feature calculation function is recorded.

[0055] Furthermore, obtaining the accuracy score Score based on the accuracy evaluation index includes:

[0056] The accuracy evaluation index IOU is expressed as follows:

[0057]

[0058] The accuracy evaluation index F1-score is expressed as follows:

[0059]

[0060] The precision score is expressed as:

[0061] Score = (0.5*IOU+0.5*F1)*100

[0062] Among them, P (Precision) and R (Recall), the formula is expressed as:

[0063]

[0064] Among them, TP is a positive sample predicted to be positive, TN is a negative sample predicted to be negative, FP is a negative sample predicted to be positive, and FN is a positive sample predicted to be negative.

[0065] Furthermore, the S4 includes: extracting cultivated land patches through the trained model, saving the patches as binary GeoTiff data, and designing the grid GeoTiff data to be saved as vector shp format data.

[0066] Beneficial effects of the present invention:

[0067] 1. The cultivated land feature extraction method of the present invention is used for cultivated land extraction from satellite time-series image data. It solves the problem that cultivated land extraction in the past was limited by satellite sources and regions, making it impossible to use heterogeneous satellite data and realize large-scale cultivated land extraction. At the same time, it has been verified that the method of the present invention can not only improve the accuracy of cultivated land extraction, but also has the ability to generalize land objects. It can also realize large-scale, heterogeneous satellite remote sensing image cultivated land extraction, improve the accuracy of cultivated land classification, determine the optimal combination of spatiotemporal spectral features, and have regional migration capabilities.

[0068] 2. The present invention addresses the problem of misclassification and omission of cultivated land categories caused by the inability of single-phase data to reflect the differences in cultivated land such as surface vegetation. It objectively determines the optimal spatial, spectral, texture, and temporal characteristics of remote sensing image data for cultivated land extraction, explores the relationship between multidimensional feature spaces such as spectra, vegetation indices, textures, and statistical characteristics in the time dimension related to cultivated land extraction, and provides an effective method for extracting spatiotemporal spectral features of cultivated land. The method has feature generalization and migration capabilities and good robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 Schematic diagram of the process of the cultivated land extraction method based on time series multispectral image features of the present invention;

[0070] Figure 2 This is a structural flow chart of the cultivated land extraction method based on time series multispectral image features of the present invention;

[0071] Figure 3 A flow chart showing the structure of preprocessing image data according to the present invention;

[0072] Figure 4 This is a structural flow chart of the present invention for performing various feature calculations and dimensionality reduction;

[0073] Figure 5 This is a structural flow chart of the model training performed by the present invention;

[0074] Figure 6 It is a structural flow chart of the application of the model of the present invention;

[0075] Figure 7 This is a sampling diagram of the test area of ​​the present invention;

[0076] Figure 8 This is a comparative display of local details of the test area of ​​the present invention;

[0077] Figure 9 Schematic diagram of the structure of the electronic device of the present invention. DETAILED DESCRIPTION

[0078] The following describes embodiments of the present invention in detail. Examples of the embodiments are shown in the accompanying drawings, wherein the same or similar symbols throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended only to explain the present invention and are not to be construed as limiting the present invention.

[0079] The existing cultivated land extraction model has defects in extracting cultivated land from remote sensing images, resulting in low cultivated land classification accuracy.

[0080] In order to solve the above problems, the present invention provides a method for extracting cultivated land based on time series multispectral image features. Figure 1 and 2 A flowchart of a method is provided for an embodiment of the present invention.

[0081] The extraction method disclosed in the present invention aims to improve the automation and accuracy of cultivated land interpretation from multispectral remote sensing images. It mainly includes core functions such as remote sensing image and sample preprocessing, feature extraction and dimensionality reduction, cultivated land feature classification, and cultivated land results output and storage. The method includes:

[0082] S1. Acquire multiple phases of Sentinel-2 image data, preprocess the image data, mark the cultivated land areas, and create cultivated land sample labels for the image data.

[0083] like Figure 3 As shown in the figure, multiple phases of Sentinel-2 complete remote sensing image data are obtained and preprocessed, including multi-phase image spatial registration, overlapping area clipping, and cultivated land labeling of Sentinel-2 remote sensing images; resampling the cropped remote sensing image data of each band, reading the central wavelength metadata of each band, and performing multi-band rearrangement; and calculating the multi-band index characteristics of the remote sensing image based on the true color band and multi-band combination data.

[0084] The cropped band data can be resampled to 10m, the central wavelength metadata of each band can be read, the true color bands can be combined, and saved as Geotiff format data, with the coordinate system selected as WGS-84.

[0085] S2. Calculate and reduce the dimension of the preprocessed image data by performing index features, spatial grayscale features, and texture features, and calculate the time series features after dimensionality reduction.

[0086] Perform multi-band feature calculations on remote sensing image data, including vegetation index features, spatial grayscale features, and texture features.

[0087] Among them, the multi-band index characteristics are shown in Table 1:

[0088] Table 1. Characteristics of multi-band vegetation index of remote sensing data

[0089]

[0090]

[0091] Among them, B, G, and R are true color bands, which are the B2, B3, and B4 bands of Sentinel-2 remote sensing images respectively. NIR, SWIR1, SWIR2, RE1, RE2, and RE3 are the B2, B3, B4, B8A, B11, B12, B5, B6, and B7 bands of Sentinel-2 remote sensing images respectively. MIR is the B11 band of Sentinel-2 remote sensing images.

[0092] The spatial grayscale feature is the texture feature obtained after the true color band is grayed. It adopts the same method as other band texture features. According to the vegetation index feature, seven different texture feature indices are extracted: homogeneity, contrast, mean, equation, correlation, dissimilarity and angular second moment (ASM). The texture features are shown in Table 2:

[0093] Table 2 Remote sensing data texture feature index

[0094]

[0095] The distance in each feature calculation is set to 1, the grayscale level is 16, and each result is the mean in the four directions of 0°, 45°, 90°, and 135°.

[0096] The above texture feature calculation is expressed in a unified formula:

[0097] p(i,j,d,θ)={[(i,j),(x+dx,y+dy)]|f(x,y)=i,f(x+dx,y+dx=j)

[0098] Where d is the specified distance between two gray levels and θ is the direction.

[0099] Time series features are statistical features of the maximum, minimum, and median values ​​of related features over multiple time periods. The calculation formula for the time series statistical features of a specific feature is as follows:

[0100]

[0101] Wherein, T is the total number of image time series, t is the current phase of the image, and F is the calculated feature value.

[0102] like Figure 4As shown, the present invention uses the Pearson correlation coefficient and the coefficient of variation as feature dimensionality reduction indicators to perform dimensionality reduction processing on various features, mainly including:

[0103] First, the index features of the image data are calculated. Based on the Pearson correlation coefficient, the index features with correlations higher than 97% are considered to be the same feature, and the dimensionality of the index features is reduced.

[0104] The Pearson correlation coefficient, as an indicator to measure the correlation between different variables, is:

[0105]

[0106] Among them, r is the correlation coefficient, X and Y are pixel values.

[0107] Then, based on the exponential features after dimensionality reduction, the texture features are calculated. Based on the Pearson correlation coefficient, texture features with correlation higher than 97% are regarded as the same feature, and the dimensionality of the texture features is reduced. At the same time, the RGB synthesized true color image is grayscaled, and the spatial grayscale features are calculated to further reduce the dimensionality of the spatial texture features.

[0108] Then, the single features after correlation screening are superimposed to calculate the time series features, and the coefficient of variation less than 15% is used as the threshold to reduce the dimension of the time series features;

[0109] The coefficient of variation, as an indicator to measure the difference in the time dimension, is:

[0110]

[0111] CV=σ / μ*100%

[0112] Among them, μ is the characteristic mean, σ is the characteristic standard deviation, and CV is the coefficient of variation.

[0113] Finally, the features obtained after dimensionality reduction include: spatial grayscale features, texture features, exponential features and time series features.

[0114] S3. Based on index features, spatial grayscale features, texture features and time series features, the feature classifier is used to screen the features that have the greatest impact on the classification results, and a cultivated land classification model is established and trained.

[0115] like Figure 5 As shown in the figure, based on the acquired spatial grayscale features, texture features, index features and time series features, combined with the prepared farmland sample label true value, as the input sample of the classifier, the above high-dimensional features are classified. The classifier is LightGBM, and the accuracy evaluation indicators are IOU and F1-score.

[0116] Obtain the accuracy score based on the accuracy evaluation indicators, including:

[0117] The accuracy evaluation index IOU is expressed as follows:

[0118]

[0119] The accuracy evaluation index F1-score is expressed as follows:

[0120]

[0121] Among them, P (Precision) and R (Recall), the formula is expressed as:

[0122]

[0123]

[0124] Among them, TP is a positive sample predicted to be positive, TN is a negative sample predicted to be negative, FP is a negative sample predicted to be positive, and FN is a positive sample predicted to be negative.

[0125] The final accuracy score is Score = (0.5*IOU+0.5*F1)*100

[0126] First, set a low number of iterations to pre-train the sample data, and stop the iteration when the average absolute change in the model accuracy score after three rounds is less than 3%.

[0127] Based on the trained model, the contribution of each input feature to classification accuracy is statistically analyzed. Gain is used as a measure of feature importance to screen features. Features whose classification accuracy decreases by no more than 5% after removal are removed as low-gain features. The final high-gain features are retained as input and the model is retrained until the maximum accuracy score improves by less than 1% after every five model training iterations. Iterations are then stopped and the input features are restored to the feature quantities calculated in step 2. A calculation function record is then made. A new farmland feature extraction model is reconstructed using the trained model, the input feature calculation function model, the classification category system, and the evaluation metrics. The newly generated model features low computational effort and high accuracy, and feature dimensionality reduction accelerates model convergence.

[0128] S4. Extract the cultivated land patches from the Sentinel-2 image data based on the trained cultivated land classification model and save the results.

[0129] like Figure 6As shown in the figure, it is mainly for the further application of the cultivated land extraction model. Based on the trained cultivated land extraction model, the time series Sentinel-2 images are input in batches. According to the recorded characteristic function operator, the image data is firstly feature extracted to extract the cultivated land patches, and then exported as binary GeoTiff to save the results. At the same time, the function of converting the raster GeoTiff data into vector shp format data for saving is added.

[0130] Application examples:

[0131] Take the local area of ​​XX City, XX Province as the experimental area, such as Figure 7 As shown in the figure, rice, winter wheat and other crops are mainly grown in the area. Multiple Sentinel-2 L2A data were acquired from the imaging platform in three areas and used for process example verification.

[0132] Data acquisition for the three regions was primarily concentrated between April and November. Seven periods of data were used for modeling in Hefei, seven periods of data for Region 1, and six periods of data for Region 2. See Table 3 for a detailed data set.

[0133] Table 3 Verification area data table

[0134]

[0135]

[0136] The original samples used for classification in the example come from experts with manual verification, and a multi-period 10-meter resolution sample database is produced; this sample database covers the entire country and includes 9 major land use and land cover types, including water bodies, grasslands, impervious surfaces, and cultivated land. It has high accuracy and can meet the needs of general surface classification research.

[0137] First, the image data was preprocessed and feature selected to obtain six vegetation index features: NDVI, NDWI, EVI, LSWI, NDTI, and REP; grayscale spatial features of 2-, 3-, and 4-band combinations; and NDVI texture and temporal features. After feature selection, the remaining features were the NDTI, NDVI, and LSWI vegetation index features; spatial features; NDVI contrast and homogeneous texture features; and MDTI temporal statistical features. The cultivated land extraction model was then trained.

[0138] like Figure 8 As shown, the data of the verification area was input and compared with the results published by ESRI. It was found that there were a large number of omissions in the ESRI labels, while the method of the present invention had a small number of misclassifications in areas with concentrated buildings. However, in terms of the precision of the classification results, the results of the present invention were better than the ESRI classification results.

[0139] The present invention also provides an electronic device, Figure 9 A schematic diagram of the structure of an electronic device provided by an embodiment of the present invention, such as Figure 9 As shown, the electronic device may include: a processor, a communications interface, a memory, and a communication bus, wherein the processor, the communications interface, and the memory communicate with each other via the communication bus. The processor may call logic instructions in the memory, for example, to execute the following method:

[0140] S1. Acquire multiple phases of Sentinel-2 image data, pre-process the image data, mark the cultivated land areas, and create cultivated land sample labels for the image data;

[0141] S2. Calculate index features, spatial grayscale features, and texture features for the preprocessed image data, perform dimensionality reduction, and calculate the time series features after dimensionality reduction;

[0142] S3. Based on index features, spatial grayscale features, texture features, and temporal features, a feature classifier is used to select the features that have the greatest impact on the classification results, and a cultivated land classification model is established and trained.

[0143] S4. Extract the cultivated land patches from the Sentinel-2 image data based on the trained cultivated land classification model and save the results.

[0144] In addition, the logical instructions in the above-mentioned memory can be implemented in the form of a software functional unit and can be stored in a computer-readable storage medium when sold or used as an independent product. Based on this understanding, the technical solution of the present invention is essentially or the part that contributes to the prior art or the part of the technical solution can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, server, or network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, and other media that can store program code.

[0145] An embodiment of the present invention further provides a non-transitory computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a processor, the method provided in each of the above embodiments is implemented, for example, including:

[0146] S1. Acquire multiple phases of Sentinel-2 image data, pre-process the image data, mark the cultivated land areas, and create cultivated land sample labels for the image data;

[0147] S2. Calculate index features, spatial grayscale features, and texture features for the preprocessed image data, perform dimensionality reduction, and calculate the time series features after dimensionality reduction;

[0148] S3. Based on index features, spatial grayscale features, texture features, and temporal features, a feature classifier is used to select the features that have the greatest impact on the classification results, and a cultivated land classification model is established and trained.

[0149] S4. Extract the cultivated land patches from the Sentinel-2 image data based on the trained cultivated land classification model and save the results.

[0150] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, i.e., they may be located in one location or distributed across multiple network units. Some or all of the modules may be selected based on actual needs to achieve the objectives of the present embodiment. Persons of ordinary skill in the art will be able to understand and implement the present invention without inventive effort.

[0151] Through the description of the above embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus a necessary general hardware platform, or of course, by hardware. Based on this understanding, the essence of the above technical solution or the part that contributes to the existing technology can be embodied in the form of a software product. The computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a magnetic disk, an optical disk, etc., and includes a number of instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment or certain parts of the embodiments.

[0152] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the various embodiments of the present invention.

Claims

1. A method for extracting cultivated land based on time series multispectral image features, characterized in that: include: S1. Acquire multiple phases of Sentinel-2 image data, pre-process the image data, mark the cultivated land areas, and create cultivated land sample labels for the image data; S2. Calculate index features, spatial grayscale features, and texture features for the preprocessed image data, perform dimensionality reduction, and calculate the time series features after dimensionality reduction; S3. Based on index features, spatial grayscale features, texture features, and temporal features, a feature classifier is used to select the features that have the greatest impact on the classification results, and a cultivated land classification model is established and trained. S4. Extracting cultivated land patches from Sentinel-2 image data based on the trained cultivated land classification model and saving the results; The step S2 calculates and reduces the dimension of the exponential feature, the spatial grayscale feature, and the texture feature, and calculates the time series feature after the dimension reduction, including the following steps: S21. Calculate the index features of the image data. Based on the Pearson correlation coefficient, consider the index features with a correlation higher than 97% as the same feature, and reduce the dimensionality of the index features. S22. Calculate texture features based on the exponential features after dimensionality reduction. Based on the Pearson correlation coefficient, consider texture features with a correlation greater than 97% as the same feature and reduce the dimensionality of the texture features. Simultaneously, grayscale the true color image and calculate the spatial grayscale features. S23, superimposing the single features after correlation screening to calculate the time series features, and using the coefficient of variation less than 15% as the threshold to reduce the dimension of the time series features; S24, obtaining features after dimensionality reduction, including spatial grayscale features, texture features, index features, and time series features; The cultivated land classification model is established and trained in S3, including: LightGBM is selected as the classifier, and IOU and F1-score are selected as the accuracy evaluation indicators. The accuracy score is obtained based on the accuracy evaluation indicators. The model is iteratively trained according to the accuracy score until the average absolute change of the accuracy score in every three rounds is less than 3%; Based on the trained model, the contribution of each feature to the classification accuracy is counted, and feature screening is performed using gain as a measurement indicator. High-gain features are retained, and the model is re-trained until the accuracy score of each five rounds of model training iterations increases by less than 1%. The iteration is stopped, and the feature calculation function is recorded.

2. The method for extracting cultivated land based on temporal multispectral image features according to claim 1, characterized in that: The image data is pre-processed in S1, including the following steps: S11. Perform multi-period image spatial registration, overlap region cropping, and farmland labeling on Sentinel-2 remote sensing images. S12, resampling the cropped remote sensing image data of each band, reading metadata of the central wavelength of each band, and performing multi-band rearrangement; S13. Calculate multi-band index features of remote sensing images based on true color band and multi-band combination data.

3. The method for extracting cultivated land based on temporal multispectral image features according to claim 2, characterized in that: Multi-band index features include: NDVI = (NIR-R) / (NIR+R) NDWI=(G-NIR) / (G+NIR) EVI=2.5*(NIR-R) / (NIR+6R-7.5B+1) <h2 style=";text-align:left;direction:ltr">MSAVI = (2NIR + 1 - sqrt ((2NIR - 1))<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> -8(NIR-R)) / 2 NDBI = (MIR - NIR) / (MIR + NIR) LSWI=(NIR-SWIR) / (NIR-SWIR) NDTI=(SWIR1-SWIR2) / (SWIR1+SWIR2) REP=(705+35*(0.5*(RE3+R)-RE1) / (RE2-RE1) Among them, B, G, and R are true color bands, which are the B2, B3, and B4 bands of Sentinel-2 remote sensing images respectively. NIR, SWIR1, SWIR2, RE1, RE2, and RE3 are the B2, B3, B4, B8A, B11, B12, B5, B6, and B7 bands of Sentinel-2 remote sensing images respectively. MIR is the B11 band of Sentinel-2 remote sensing images.

4. The method for extracting cultivated land based on temporal multispectral image features according to claim 1, characterized in that: The spatial grayscale features and texture features include: Homogeneity: Contrast ratio: Mean: variance: Dependencies: Dissimilarity: Angular second moment: The distance in each feature calculation is set to 1, the grayscale level is 16, and each result is the mean value in the four directions of 0°, 45°, 90°, and 135°, which is expressed in a unified formula: p(i,j,d,θ)={[(i,j),(x+dx,y+dy)]|f(x,y)=i,f(x+dx,y+dx=j) Where d is the specified distance between two gray levels and θ is the direction.

5. The method for extracting cultivated land based on time series multispectral image features according to claim 1, characterized in that: The time series characteristic formula is expressed as: Wherein, T is the total number of image time series, t is the current phase of the image, and F is the calculated feature value.

6. The method for extracting cultivated land based on time series multispectral image features according to claim 5, characterized in that: The feature dimension reduction index includes the Pearson correlation coefficient and the coefficient of variation, wherein, Pearson correlation coefficient, the formula is expressed as: Among them, r is the correlation coefficient, X and Y are pixel values; The coefficient of variation is expressed as: CV=σ / μ*100% μ is the characteristic mean, σ is the characteristic standard deviation, and CV is the coefficient of variation.

7. The method for extracting cultivated land based on time series multispectral image features according to claim 6, characterized in that: The step of obtaining the accuracy score based on the accuracy evaluation index includes: The accuracy evaluation index IOU is expressed as follows: The accuracy evaluation index F1-score is expressed as follows: The precision score is expressed as: Score = (0.5*IOU+0.5*F1)*100 Among them, P (Precision) and R (Recall), the formula is expressed as: Among them, TP is a positive sample predicted to be positive, TN is a negative sample predicted to be negative, FP is a negative sample predicted to be positive, and FN is a positive sample predicted to be negative.

8. The method for extracting cultivated land based on time series multispectral image features according to claim 1, characterized in that: The S4 includes: extracting cultivated land patches through the trained model, saving them as binary GeoTiff data, and designing the grid GeoTiff data to save them as vector shp format data.

Images