Cultivated land extraction method based on time sequence multispectral image features

Abstract

The invention discloses a cultivated land extraction method based on time sequence multispectral image features, and the method comprises the steps: obtaining multi-stage Sentinel-2 image data, carrying out the preprocessing, and making an image data cultivated land sample label; carrying out index feature, spatial gray feature and texture feature calculation and dimension reduction on the image data, and calculating time sequence features; based on the index features, the spatial gray features, the texture features and the time sequence features, screening the features through a feature classifier, establishing a cultivated land classification model, and training the cultivated land classification model; and carrying out farmland pattern spot extraction on the Sentinel-2 image data. According to the method, the optimal remote sensing image data space, spectrum, texture and time sequence features are extracted from cultivated land, the relation among multi-dimensional feature spaces such as the spectrum, the vegetation index, the texture and statistical features in the time dimension related to cultivated land extraction is mined, the effective cultivated land space-time spectrum feature extraction method is provided, the feature generalization and migration ability is achieved, and robustness is good.

Description

Technical Field

[0001] The present invention relates to the technical field of cultivated land extraction based on multi-spectral images, and in particular to a method for extracting cultivated land based on the characteristics of time-series multi-spectral images. Background Art

[0002] The extraction of cultivated land scope and spatial distribution has important influences and functions on crop yield assessment, food safety, etc. Remote sensing technology is considered an effective tool that can repeatedly provide accurate and objective cultivated land information in terms of time and space. Currently, there are various satellite remote sensing data sources that can be used for the extraction of cultivated land information, such as MODIS, Landsat, and Sentinel-2. Remote sensing cultivated land detection mainly focuses on the investigation of cultivated land quantity and spatial distribution, and its technical methods have developed from visual interpretation to semi-automation based on statistical models, and gradually turned to data-driven classification methods mainly based on machine learning and deep learning. Different from general surface types, such as buildings, roads, and impervious surfaces, the cultivated land concept includes a relatively wide range of ground objects. Different crops may be planted within the cultivated land scope at the same time, and the imaging differences of their images are relatively obvious, resulting in relatively large inter-class differences. In addition, the phenomenon that cultivated land changes significantly between years leads to relatively low classification accuracy of existing methods in local areas, and the degree of manual participation is still relatively high during the implementation process. The selection of classification features and labels as the input of the model directly affects the accuracy of cultivated land extraction. Currently, there is no effective fixed feature combination method, making it difficult to perform feature generalization and migration, and lacking robustness.

[0003] For example: The invention application with the application number 202410084966.1 discloses a method for extracting cultivated land information from high-resolution images based on object orientation. This application comprehensively considers the image spatial resolution, the complexity of ground object distribution, the fragmentation degree of the region, etc., and based on the shape features of roads and field ridges, extracts the road and field ridge areas, which are not only used to distinguish the boundaries of each cultivated land, but also obtain a more refined cultivated land area. However, its solution also has the problem that different crops may be planted within the cultivated land scope at the same time, the imaging differences of their images are relatively obvious, and the classification of cultivated land is poor.

[0004] Therefore, aiming at the problem of misclassification and missed classification of cultivated land categories caused by the inability of single-phase time-series data to reflect the differences of cultivated land such as surface vegetation, objectively determining the optimal remote sensing image data spatial, spectral, texture, and time-series features for cultivated land extraction, the present invention explores the relationships between multi-dimensional feature spaces such as spectra, vegetation indices, textures, and statistical features in the time dimension related to cultivated land extraction, and discloses an effective method for extracting the spatio-temporal spectral features of cultivated land based on Sentinel-2, which can achieve the ability of cultivated land feature generalization and migration and has good robustness. Summary of the Invention

[0005] In view of the above problems, the purpose of the present invention is to provide a cultivated land extraction method based on temporal multi-spectral image features, which can realize the extraction of cultivated land from satellite remote sensing images with large-scale and heterogeneous data, improve the classification accuracy of cultivated land, determine the optimal spatio-temporal spectral feature combination, and have regional migration ability.

[0006] An embodiment of the present invention provides a cultivated land extraction method based on temporal multi-spectral image features, including the steps of:

[0007] S1. Obtain multi-temporal Sentinel-2 image data, preprocess the image data, label the cultivated land area, and make the cultivated land sample labels of the image data;

[0008] S2. Calculate and reduce the dimensionality of the index features, spatial gray-scale features, and texture features of the preprocessed image data, and calculate the temporal features after dimensionality reduction;

[0009] S3. Based on the index features, spatial gray-scale features, texture features, and temporal features, screen the features that have the greatest impact on the classification results through a feature classifier, establish a cultivated land classification model and train it;

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

[0011] Further, the preprocessing of the image data in S1 includes the steps of:

[0012] S11. Perform multi-temporal image spatial registration, overlapping area cropping, and cultivated land annotation on the Sentinel-2 remote sensing image;

[0013] S12. Resample the cropped remote sensing image data of each band, read the metadata of the central wavelength of each band, and perform multi-band rearrangement;

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

[0015] Further, 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] Wherein, B, G, and R are true color bands, which are the B2, B3, and B4 bands of the Sentinel-2 remote sensing image respectively. NIR, SWIR1, SWIR2, RE1, RE2, and RE3 are the B2, B3, B4, B8A, B11, B12, B5, B6, and B7 bands of the Sentinel-2 remote sensing image respectively. MIR is the B11 band of the Sentinel-2 remote sensing image.

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

[0026] Homogeneity:

[0027] Contrast:

[0028] Mean:

[0029] Variance:

[0030] Correlation:

[0031] Dissimilarity:

[0032] Angular second moment:

[0033] Wherein, the distance is set to 1 in the calculation of each feature, the gray level is 16, and each result is the mean value in four directions of 0°, 45°, 90°, and 135°. It is expressed by the unified formula as:

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

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

[0036] Further, the formula for the temporal feature is expressed as:

[0037]

[0038] Among them, T is the total number of image timings, t is the current phase of the image, and F is the calculated feature quantity.

[0039] Further, in S2, the exponential feature, spatial gray-level feature, and texture feature are calculated and dimensionally reduced, and the temporal feature after dimensional reduction is calculated, including the steps:

[0040] S21. Calculate the exponential feature of the image data. Based on the Pearson correlation coefficient, the exponential features with a correlation higher than 97% are regarded as the same feature, and the exponential feature is dimensionally reduced;

[0041] S22. Based on the dimensionally reduced exponential feature, calculate the texture feature. Based on the Pearson correlation coefficient, the texture features with a correlation higher than 97% are regarded as the same feature, and the texture feature is dimensionally reduced; at the same time, the true-color image is grayed out, and the spatial gray-level feature is calculated;

[0042] S23. Superimpose and calculate the temporal feature for the single feature after correlation screening, and use a coefficient of variation lower than 15% as the threshold to dimensionally reduce the temporal feature;

[0043] S24. Obtain the features after dimensional reduction, including the spatial gray-level feature, texture feature, exponential feature, and temporal feature.

[0044] Further, the feature dimensional reduction indicators include the Pearson correlation coefficient and the coefficient of variation. Among them,

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

[0046]

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

[0048] The coefficient of variation, the formula is expressed as:

[0049]

[0050] CV = σ / μ * 100%

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

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

[0053] Select LightGBM as the classifier, choose IOU and F1-score as the accuracy evaluation metrics, obtain the accuracy score Score based on the accuracy evaluation metrics, and iteratively train the model according to the accuracy score Score until the average absolute change of the accuracy score Score is less than 3% every three rounds;

[0054] According to the trained model, statistically analyze the contribution degree of each feature to the classification accuracy, use gain as the measurement index for feature screening, retain the high-gain features, and re-iteratively train the model until the improvement of the accuracy score Score of each five-round model training iteration is less than 1%, then stop the iteration and record the feature calculation function.

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

[0056] The accuracy evaluation metric IOU is expressed by the formula:

[0057]

[0058] The accuracy evaluation metric F1-score is expressed by the formula:

[0059]

[0060] The accuracy score is expressed by the formula:

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

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

[0063]

[0064] Among them, TP is the positive sample predicted as positive, TN is the negative sample predicted as negative, FP is the negative sample predicted as positive, and FN is the positive sample predicted as negative.

[0065] Furthermore, in S4, it includes: extracting cultivated land patches through the trained model, saving them as binary GeoTiff data, and designing to save the raster GeoTiff data as vector shp format data.

[0066] The beneficial effects of the present invention:

[0067] 1. The cultivated land feature extraction method of the present invention is used for extracting cultivated land from satellite time-series image data, solving the problems that the previous cultivated land extraction was restricted by satellite sources and regions, unable to use heterogeneous satellite data and achieve large-scale cultivated land extraction. At the same time, through verification, the method of the present invention can not only improve the accuracy of cultivated land extraction, but also has the ability of feature generalization of ground objects, and can also achieve large-scale and heterogeneous data satellite remote sensing image cultivated land extraction, improve the classification accuracy of cultivated land, determine the best spatio-temporal spectral feature combination, and has the ability of regional migration.

[0068] 2. The present invention aims at the problem of misclassification and omission of cultivated land categories caused by the inability of single-phase time-series data to reflect the differences of cultivated land such as surface vegetation. It objectively determines the spatial, spectral, texture, and time-series features of remote sensing image data that are optimal for cultivated land extraction, explores the relationships between multi-dimensional feature spaces such as spectral, vegetation index, texture, and statistical features in the time dimension related to cultivated land extraction, provides an effective method for extracting spatio-temporal spectral features of cultivated land, has the ability of feature generalization and migration, and has good robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 It is a flow chart of the cultivated land extraction method based on time-series multi-spectral image features of the present invention;

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

[0071] Figure 3 It is a structural flow chart for preprocessing the image data of the present invention;

[0072] Figure 4 It is a structural flow chart for calculating and dimension-reducing various features of the present invention;

[0073] Figure 5 It is a structural flow chart for model training of the present invention;

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

[0075] Figure 7 It is a sampling map of the example test area of the present invention;

[0076] Figure 8 It is a local detail comparison display map of the example test area of the present invention;

[0077] Figure 9 It is a structural schematic diagram of the electronic device of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0078] Embodiments of the present invention will be described in detail below. Examples of the embodiments are shown in the accompanying drawings, where the same or similar symbols represent the same or similar elements or elements with the same or similar functions throughout. The embodiments described below by referring to the accompanying drawings are exemplary and are only used to explain the present invention and should not be construed as a limitation to the present invention.

[0079] Existing arable land extraction models have defects in the extraction of arable land from remote sensing images, resulting in low classification accuracy of arable land.

[0080] In view of the above problems, the present invention provides a method for extracting arable land based on the characteristics of multi-temporal multispectral images. Figure 1 and 2 is a schematic flowchart of the method provided by the embodiment of the present invention.

[0081] The extraction method disclosed by the present invention aims to improve the automation degree and accuracy of the interpretation of arable land in multi-spectral remote sensing images, and mainly includes core functions such as remote sensing image and sample preprocessing functions, feature extraction and dimensionality reduction functions, arable land feature classification, arable land result output and storage, etc. The method includes:

[0082] S1. Obtain multi-temporal Sentinel-2 image data, preprocess the image data, label the arable land area, and make arable land sample labels for the image data.

[0083] As Figure 3 shown, obtain multi-temporal complete Sentinel-2 remote sensing image data, and preprocess the remote sensing image data, including: performing multi-temporal image spatial registration, overlapping area cropping and arable land annotation on the Sentinel-2 remote sensing image; resampling the cropped remote sensing image data of each band, reading the central wavelength metadata of each band, and performing multi-band rearrangement; calculating multi-band index features for the remote sensing image based on the true color band and multi-band combined data.

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

[0085] S2. Calculate and reduce the dimensionality of index features, spatial gray features and texture features for the preprocessed image data, and calculate the temporal features after dimensionality reduction.

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

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

[0088] Table 1 Multi-band vegetation index feature table 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 the Sentinel-2 remote sensing image respectively. NIR, SWIR1, SWIR2, RE1, RE2, and RE3 are the B2, B3, B4, B8A, B11, B12, B5, B6, and B7 bands of the Sentinel-2 remote sensing image respectively. MIR is the B11 band of the Sentinel-2 remote sensing image.

[0092] The spatial gray-level feature is the texture feature obtained after graying the true-color bands. It uses the same method as the texture features of other bands. Based on the vegetation index feature, seven different texture feature indices, namely homogeneity, contrast, mean, variance, correlation, dissimilarity, and angular second moment (ASM), are extracted. The texture features are shown in Table 2:

[0093] Table 2 Texture Feature Index Table of Remote Sensing Data

[0094]

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

[0096] The above texture feature calculation is represented by a unified formula as:

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

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

[0099] The temporal feature is the statistical feature of the maximum value, minimum value, and median of the relevant features in the multi-temporal dimension. The calculation formula for the temporal statistical feature of a specific above-mentioned feature is:

[0100]

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

[0102] Such as Figure 4As shown in the figure, 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, calculate the index features of the image data. Based on the Pearson correlation coefficient, the index features with a correlation higher than 97% are regarded as the same feature, and the index features are dimensionally reduced;

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

[0105]

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

[0107] Then, based on the dimensionally reduced index features, calculate the texture features. Based on the Pearson correlation coefficient, the texture features with a correlation higher than 97% are regarded as the same feature, and the texture features are dimensionally reduced; at the same time, gray-scale the true-color image synthesized by RGB, calculate the spatial gray-scale features, and further reduce the dimensionality of the spatial texture features.

[0108] Then, superimpose and calculate the time-series features of the single features screened by correlation, and use a coefficient of variation lower than 15% as the threshold to reduce the dimensionality of the time-series features;

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

[0110]

[0111] CV = σ / μ * 100%

[0112] where μ is the feature mean, σ is the feature standard deviation, and CV is the coefficient of variation.

[0113] Finally, the features after dimensionality reduction obtained include: spatial gray-scale features, texture features, index features, and time-series features.

[0114] S3. Based on the index features, spatial gray-scale features, texture features, and time-series features, screen the features that have the greatest impact on the classification result through a feature classifier, establish a cultivated land classification model, and perform training.

[0115] As Figure 5 shown, based on the obtained spatial gray-scale features, texture features, index features, and time-series features, combined with the true values of the cultivated land sample labels made, as the input samples of the classifier, classify the above high-dimensional features. The classifier selects LightGBM, and the accuracy evaluation indicators select IOU and F1-score.

[0116] Based on the accuracy evaluation indicators, obtain the accuracy score Score, including:

[0117] The precision evaluation index IOU, which is expressed by the formula as follows:

[0118]

[0119] The precision evaluation index F1-score, which is expressed by the formula as follows:

[0120]

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

[0122]

[0123]

[0124] Among them, TP is the positive sample predicted as positive, TN is the negative sample predicted as negative, FP is the negative sample predicted as positive, and FN is the positive sample predicted as negative.

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

[0126] First, set a relatively low number of iterations to pre-train the sample data. When the average absolute change of the model precision Score per three rounds is lower than 3%, stop the iteration.

[0127] According to the trained model, statistically analyze the contribution degree of each input feature to the classification precision. Use gain as the index to measure the feature importance, and screen the features. Remove the features whose classification precision index drops by no more than 5% after removing this feature as low-gain features, and retain the final high-gain features as the input. Retrain the model until the maximum precision Score of the model training iteration per five rounds increases by less than 1%. Then stop the iteration, restore the input features to the feature quantity calculated in step two, and record the calculation function. The trained model, the input feature calculation function model, the classification category system, and the evaluation index are reconstituted into a new cultivated land feature extraction model. The newly generated model has the characteristics of low computational complexity and high precision. After feature dimensionality reduction, the model convergence speed will be accelerated.

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

[0129] Such as Figure 6As shown, it is mainly for the further application of the cultivated land extraction model. Based on the trained cultivated land extraction model, time-series Sentinel-2 images are batch-input. According to the recorded feature function operators, the feature extraction of the image data is first carried out to extract cultivated land patches, which are exported as binary GeoTiff to save the results. At the same time, the function of converting raster GeoTiff data into vector shp format data for saving is added.

[0130] Example of application:

[0131] Taking a local area of XX City in XX Province as the experimental area, as Figure 7 shown, the main crops planted in the area are rice, winter wheat and other food crops. Multi-period Sentinel-2, L2A-level data are obtained in three areas through the image platform and are respectively used for the verification of the process instance.

[0132] The time for obtaining data in the three areas mainly focuses on April to November. Among them, there are 7 periods of data for Hefei area used for modeling, 7 periods of data for area 1 for verification, and 6 periods of data for area 2 for verification. The specific data list is shown in Table 3.

[0133] Table 3 Data table of verification areas

[0134]

[0135]

[0136] In the example, the original samples used for classification are sourced from experts supplemented by manual verification to produce a multi-period sample database with a resolution of 10 meters; its sample database covers the whole country and includes 9 main land use and land cover types such as water bodies, grasslands, impervious surfaces, and cultivated land, with high accuracy and can meet the needs of general surface classification research.

[0137] First, the image data is preprocessed, and feature selection is carried out on the image data to obtain vegetation index features such as NDVI, NDWI, EVI, LSWI, NDTI, and REP 6, grayscale spatial features of the 2, 3, and 4 band combinations, NDVI texture features, and time-series features. After feature screening, finally, three vegetation index features of NDTI, NDVI, and LSWI, spatial features, contrast and homogeneity texture features of NDVI, and time-series statistical features of MDTI are retained, and the cultivated land extraction model is trained.

[0138] As Figure 8 shown, the data of the verification area is input and compared and analyzed with the results published by ESRI. It is found that there are a large number of missed classification phenomena in the ESRI labels, while in the areas where buildings are concentrated in the method of the present invention, there are a small number of misclassification cases, but in terms of the fineness of the classification results, the results of the present invention are better than the ESRI classification results.

[0139] The present invention also provides an electronic device. Figure 9 It is a schematic structural diagram of the electronic device provided by an embodiment of the present invention. As Figure 9 shown, the electronic device may include: a processor, a communications interface, a memory, and a communication bus. Among them, the processor, the communications interface, and the memory complete mutual communication through the communication bus. The processor can call logical instructions in the memory to execute the following methods, for example:

[0140] S1. Obtain multi-temporal Sentinel-2 image data, preprocess the image data, label the cultivated land area, and make cultivated land sample labels for the image data;

[0141] S2. Calculate and reduce the dimensionality of the index features, spatial gray features, and texture features of the preprocessed image data, and calculate the temporal features after dimensionality reduction;

[0142] S3. Based on the index features, spatial gray features, texture features, and temporal features, screen the features that have the greatest impact on the classification result through a feature classifier, establish a cultivated land classification model, and train it;

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

[0144] In addition, when the logical instructions in the above-mentioned memory are implemented in the form of software functional units and sold or used as an independent product, they can be stored in a computer-readable storage medium. Based on such an understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or a part of this 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 causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in various embodiments of the present invention. The aforementioned storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memories (ROMs), random access memories (RAMs), magnetic disks, or optical discs that can store program codes.

[0145] An embodiment of the present invention also provides a non-transitory computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, it is implemented to execute the methods provided in the above-mentioned embodiments, for example, including:

[0146] S1. Obtain multi - period Sentinel - 2 image data, pre - process the image data, label the cultivated land areas, and make the cultivated land sample labels for the image data;

[0147] S2. Calculate and reduce the dimensionality of the exponential features, spatial gray - level features, and texture features for the pre - processed image data, and calculate the time - series features after dimensionality reduction;

[0148] S3. Based on the exponential features, spatial gray - level features, texture features, and time - series features, screen the features that have the greatest impact on the classification results through a feature classifier, establish a cultivated land classification model, and train it;

[0149] S4. Extract the cultivated land patches from the Sentinel - 2 image data according to 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 separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed to multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of this embodiment. Those of ordinary skill in the art can understand and implement it without creative work.

[0151] Through the description of the above - mentioned implementation manners, those skilled in the art can clearly understand that each implementation manner can be realized by means of software plus a necessary general - purpose hardware platform, and of course, it can also be realized by hardware. Based on such an understanding, the essence of the above - mentioned technical solution, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer - readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to enable 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 some 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 them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the 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 area, and create cultivated land sample labels for the image data; S2, calculating and reducing the dimension of the index feature, spatial grayscale feature and texture feature of the preprocessed image data, and calculating the time series feature after the dimension reduction; S3, based on index features, spatial grayscale features, texture features and time series features, the features that have the greatest impact on the classification results are screened through feature classifiers, and a cultivated land classification model is established and trained; S4. Extract the cultivated land patches from the Sentinel-2 image data based on the trained cultivated land classification model and save the results.

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

3. The method for extracting cultivated land based on time series 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; and MIR is the B11 band of Sentinel-2 remote sensing images.

4. The method for extracting cultivated land based on time series multispectral image features according to claim 1, characterized in that: The spatial grayscale features and texture features include: Homogeneity: Contrast ratio: Mean: variance: Relevance: Dissimilarity: Angular second moment: Among them, the distance in each feature calculation is set to 1, the gray level is 16, and each result is the mean value in the four directions of 0°, 45°, 90°, and 135°, which is expressed by a unified formula: p(i,j,d,θ)={[(i,j),(x+dx,y+dy)]|f(x,y)=i,f(x+dx,y+dx=j) Among them, 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 timing characteristic formula is expressed as: Among them, 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 1, characterized in that: In the step S2, the exponential feature, the spatial grayscale feature and the texture feature are calculated and dimensionally reduced, and the time series feature after dimension reduction is calculated, including the steps of: S21, calculating the index features of the image data, based on the Pearson correlation coefficient, treating the index features with correlation higher than 97% as the same feature, and reducing the dimension of the index features; S22, based on the index features after dimensionality reduction, the texture features are calculated, and based on the Pearson correlation coefficient, the texture features with correlation higher than 97% are regarded as the same feature, and the texture features are reduced in dimensionality; at the same time, the true color image is grayed and the spatial grayscale features are calculated; S23, superimposing the single features after correlation screening to calculate the time series features, and using the coefficient of variation less than 15% as a threshold to reduce the dimension of the time series features; S24. Obtain features after dimensionality reduction, including spatial grayscale features, texture features, exponential features, and time series features.

7. The method for extracting cultivated land based on time series multispectral image features according to claim 6, characterized in that: The feature dimension reduction index includes the Pearson correlation coefficient and the coefficient of variation, where: Pearson correlation coefficient, the formula is: 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.

8. The method for extracting cultivated land based on time series multispectral image features according to claim 1, characterized in that: The cultivated land classification model is established and trained in S3, including: Select LightGBM as the classifier, select IOU and F1-score as the accuracy evaluation indicators, obtain the accuracy score based on the accuracy evaluation indicators, and iteratively train the model according to the accuracy score until the average absolute change of the accuracy score in every three rounds is less than 3%; According to the trained model, the contribution of each feature to the classification accuracy is counted, and gain is used as a measurement indicator for feature screening. High-gain features are retained, and the model is re-trained until the accuracy score Score of every five rounds of model training iterations increases by less than 1%. Then, the iteration is stopped and the feature calculation function is recorded.

9. The method for extracting cultivated land based on time series multispectral image features according to claim 8, characterized in that: The step of obtaining the accuracy score 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.

10. 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