A rapid crop mapping method for large areas based on the integration of multi-source remote sensing data

Abstract

The present invention discloses a large-area crop rapid mapping method integrating multi-source remote sensing data, obtaining multi-source remote sensing satellite images and preprocessing them, amplifying samples based on the spectral similarity between field-collected sample points and amplified sample points, and screening amplified sample points based on dynamic time warping distance; drawing a reference map of crop distribution in a study sub-area using collaborative multi-source data; constructing a sub-pixel mapping sample set for the study area based on the reference map of crop distribution in the study sub-area and predicting sub-pixel distribution maps of various types of crops. The present invention makes full use of existing field-collected sample points, explores the mapping potential of multi-source remote sensing data and different classifiers, and realizes rapid sub-pixel mapping of crops in a large-scale study area.

Description

Technical Field

[0001] The invention belongs to the technical field of agricultural remote sensing, and in particular relates to a large-area crop rapid mapping method integrating multi-source remote sensing data. Background Art

[0002] Remote sensing data has the advantages of large area, long time series, and synchronous earth observation, and is an important basis for accurately extracting spatial distribution information of crops. However, due to technical limitations and interference from factors such as cloud and rain weather, a single remote sensing data source is difficult to meet the needs of time series images in crop mapping. The existing optical remote sensing images can be divided into medium and high spatial resolution images and low spatial resolution images. Medium and high spatial resolution images (such as Landsat series satellite images, Sentinel-2 satellite images, and domestic high-resolution series satellite images) are the main data sources for remote sensing crop mapping, but they are easily affected by clouds and rain and lack historical data, which limits the accuracy of crop mapping and the time series range of crop mapping in related studies. Low spatial resolution images (such as MODIS satellite images) have the advantages of short re-entry period, large imaging width, and rich accumulation of historical image data. They have been widely used in the field of global land surface environment monitoring, but their low spatial resolution causes them to be greatly affected by mixed pixels. Therefore, how to effectively coordinate multi-source remote sensing data sources to achieve large-area crop mapping is a key scientific problem to be solved.

[0003] In addition, the traditional method of obtaining field information is the field sampling method, which is to collect sample points on the spot through field surveys and GPS equipment. The samples obtained by this method are of high quality, but the sampling cost is too high, and it is impossible to collect field samples in the past, so it is difficult to meet the needs of long-term crop mapping in large areas. Therefore, how to generate high-quality samples of the target year and target area based on existing samples is a key problem that needs to be solved in the field of crop mapping.

[0004] In summary, large-area crop mapping currently faces two main problems. First, in terms of sample acquisition, it is difficult to quickly and cheaply obtain mapping samples in a large-scale study area. Second, in terms of remote sensing image data utilization, medium and high spatial resolution images face the problem of long revisit cycles, small coverage area of ​​single images, and susceptibility to cloudy and rainy weather images. At the same time, it is difficult to efficiently obtain sufficient and evenly distributed samples. Low spatial resolution images face the problem of various types of crops being highly mixed within the pixels, and the use of hard classification methods is prone to misclassification and omission. Summary of the invention

[0005] The purpose of the present invention is to provide a method for rapid mapping of large-area crops by integrating multi-source remote sensing data, in order to address the problem that high mixing of pixels in medium and high resolution images may easily lead to misclassification or omission of crops, and that traditional sample acquisition methods are difficult to efficiently obtain samples for the target year.

[0006] The above-mentioned purpose of the present invention is achieved by the following technical means:

[0007] The method for rapid mapping of large-area crops by integrating multi-source remote sensing data includes the following steps:

[0008] Step 1: Obtain the original low spatial resolution satellite images of the study area corresponding to the time span of the study area and the original medium and high spatial resolution satellite images of the study sub-area corresponding to the time span of the study sub-area, and preprocess them respectively to obtain the preprocessed study area time series image dataset and the preprocessed study sub-area time series image dataset;

[0009] Step 2: Obtain the field sample point data of each research sub-area in historical years and the field sample point data of the target year; obtain the crop distribution map of each research sub-area in historical years; obtain the land use data set of each research sub-area in the target year, and extract the cultivated land layer from the land use data set; obtain the pre-processed high spatial resolution satellite images corresponding to the phenological key period of each research sub-area in historical years and the pre-processed high spatial resolution satellite images corresponding to the phenological key period of the target year;

[0010] Step 3: For the research sub-areas where the number of field sample points corresponding to the crop type in the target year is less than the field sample point threshold, sample points are amplified based on pixel spectral similarity, and the reliability of the amplified sample points is verified and screened to obtain reliable amplified sample points;

[0011] Step 4: Based on the field sample points of the target year and the amplified sample points obtained in step 3, combined with the pre-processed time series image dataset of the study sub-area, obtain a reference map of crop distribution in the study sub-area in the target year;

[0012] Step 5: randomly generate multiple sub-pixel mapping sample points, and construct a sub-pixel mapping sample set of the study area based on the sub-pixel mapping sample points, the pre-processed time series image dataset of the study area, the land use dataset, and the crop distribution reference map;

[0013] Step 6: Construct a mapping model, and use the sub-pixel mapping sample set of the study area to train and verify the screening mapping model to obtain the crop mapping optimization model corresponding to each crop type. Input the preprocessed study area time series image dataset into the crop mapping optimization model to obtain the sub-pixel distribution map corresponding to each crop type.

[0014] The pretreatment is performed in step 1 as described above, specifically comprising the following steps:

[0015] The original low spatial resolution satellite images were declouded and median synthesized to obtain a pre-processed time series image dataset of the study area.

[0016] For original medium and high spatial resolution satellite images,

[0017] When one kind of original medium-high spatial resolution satellite image is obtained, the original medium-high spatial resolution satellite image is de-clouded and median-synthesized to obtain a time series image dataset of the research sub-area after preprocessing;

[0018] When a variety of original medium and high spatial resolution satellite images are obtained, the various original medium and high spatial resolution satellite images are declouded to obtain the corresponding medium and high spatial resolution satellite images after declouding, the medium and high spatial resolution satellite images after declouding are processed by reflectivity standardization to obtain medium and high spatial resolution satellite images after reflectivity standardization, and the medium and high spatial resolution satellite images after reflectivity standardization are synthesized by median value to obtain the time series image dataset of the research sub-area after preprocessing;

[0019] Among them, the reflectivity standardization processing is performed on the cloud-removed medium- and high-spatial-resolution satellite images to obtain the medium- and high-spatial-resolution satellite images after reflectivity standardization processing, which specifically includes the following steps:

[0020] Taking a declouded medium-to-high spatial resolution satellite image as a reference, the spectral bands of other declouded medium-to-high spatial resolution satellite images are transformed through the conversion coefficient calculated by least squares regression to obtain the medium-to-high spatial resolution satellite image after spectral band conversion; the collection of the declouded medium-to-high spatial resolution satellite image as a reference and the medium-to-high spatial resolution satellite image after spectral band conversion is used as the medium-to-high spatial resolution satellite image after reflectance standardization.

[0021] As mentioned above, step 3 includes the following steps:

[0022] Step 3.1 For the study sub-areas where the number of field sample points corresponding to the crop type in the target year is less than the field sample point threshold, randomly generate preliminary amplification sample points, and the minimum distance between the preliminary amplification sample points is not less than the spatial resolution of the original medium- and high-spatial-resolution satellite image;

[0023] Overlay the crop distribution map of each research sub-area in historical years obtained in step 2 with the prepared amplified sample points to obtain the crop types of the prepared amplified sample points in historical years;

[0024] Step 3.2: Perform false color synthesis on the pre-processed high spatial resolution satellite images corresponding to the phenological critical period of each study sub-area in the historical years and the pre-processed high spatial resolution satellite images corresponding to the phenological critical period of the target year, respectively, to obtain the false color median composite images corresponding to the phenological critical period of each study sub-area in the historical years and the false color median composite images corresponding to the phenological critical period of the target year; compare the changes of each prepared amplification sample point from the pixels in the false color median composite images corresponding to the phenological critical period of the historical years to the pixels in the false color median composite images corresponding to the phenological critical period of the target year; if the changes are not obvious, the corresponding prepared amplification sample point is retained as the amplification sample point, and the crop type corresponding to the amplification sample point is retained; if the changes are obvious, the corresponding prepared amplification sample point is deleted;

[0025] Step 3.3: Based on the preprocessed time series image dataset of the study sub-area obtained in step 1, for each crop type in each study sub-area, the time series data of the normalized vegetation index NDVI corresponding to the field collection sample points and the time series data of the normalized vegetation index NDVI corresponding to the amplified sample points are obtained respectively, and the time series data of the normalized vegetation index NDVI includes the NDVI value at each time series point;

[0026] Set the upper and lower quantiles;

[0027] Step 3.4: For each crop type in each study sub-area, count the values ​​of the NDVI values ​​at the upper quantile of all field-collected sample points, and record them as the NDVI upper limit of the crop type; and count the values ​​of the NDVI values ​​at the lower quantile of all field-collected sample points, and record them as the NDVI lower limit of the crop type;

[0028] Step 3.5 For each crop type in each study sub-area, for each current amplified sample point, compare the NDVI value of the amplified sample point at each time point with the crop type NDVI upper limit value and the crop type NDVI lower limit value. If the NDVI value of the amplified sample point at any time point is lower than the crop type NDVI lower limit value or higher than the crop type NDVI upper limit value, then remove the amplified sample point;

[0029] Step 3.6 In each study sub-area: for each crop type, the corresponding dynamic time warping distance is calculated based on the time series data of the normalized vegetation index NDVI at two different field-collected sample points, which is recorded as the dynamic time warping distance DTW_cr_c&c; corresponding to the pairwise combination of all field-collected sample points of the same crop type, all dynamic time warping distances DTW_cr_c&c corresponding to the crop type are obtained; at the same time, for each crop type, the corresponding dynamic time warping distance is calculated based on the time series data of the normalized vegetation index NDVI at a field-collected sample point and the time series data of the normalized vegetation index NDVI at an amplified sample point, which is recorded as the dynamic time warping distance DTW_cr_c&e, corresponding to the pairwise combination of all field-collected sample points and amplified sample points of the same crop type, all dynamic time warping distances DTW_cr_c&e corresponding to the crop type are obtained;

[0030] For each crop type in each study sub-area: determine whether the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e. If so, it is considered that the amplified sample points in the corresponding crop type in the corresponding study sub-area are reliable, and then proceed to the subsequent step 4; if not, reduce the upper quantile and increase the lower quantile, and return to step 3.4.

[0031] Step 4 as described above includes the following steps:

[0032] Step 4.1 For each crop type in each study sub-area, combine the field collection sample points and the amplification sample points in the target year to obtain the crop sample points of the corresponding crop type in each study sub-area;

[0033] Step 4.2: Calculate the time series characteristics of vegetation indices at the pixel scale based on the preprocessed time series image dataset of the study sub-area. Vegetation indices include: Normalized Difference Vegetation Index NDVI, Enhanced Vegetation Index EVI, Land Surface Water Index LSWI, Normalized Thermal Difference Index NDTI, and Normalized Difference Decay Index NDSVI.

[0034] Step 4.3: For each crop type in each study sub-area, the time series characteristics of the vegetation index at the crop sample point and pixel scale are combined, and a pixel-by-pixel prediction is performed on the cultivated land layer obtained in step 2 based on the random forest model to obtain a reference map of crop distribution in the study sub-area in the target year.

[0035] Step 5 as described above includes the following steps:

[0036] Step 5.1 For each study sub-area, randomly generate n sub-pixel mapping sample points, and the spacing between sub-pixel mapping sample points should not be less than the spatial resolution of the original low spatial resolution satellite image;

[0037] Step 5.2 uses the location information of the sub-pixel mapping sample points to extract the spectral band characteristics of the corresponding sub-pixel mapping sample points from the pre-processed time series image dataset of the study area;

[0038] Step 5.3 Overlay the land use dataset of the target year in step 2 and the preprocessed time series image dataset of the study area. Use the sub-pixels of the preprocessed time series image dataset of the study area to count the percentage of natural vegetation sub-pixels and non-vegetation sub-pixels in the sub-pixels where the sub-pixel mapping sample points are located. The percentage of natural vegetation sub-pixels and the percentage of non-vegetation sub-pixels are used as non-spectral features.

[0039] Step 5.4 Overlay the crop distribution reference map of the study sub-area in the target year obtained in step 4 with the pre-processed study area time series image dataset, and for each crop type, calculate the percentage of crop sub-pixel at each sub-pixel mapping sample point in turn;

[0040] Step 5.5: Combine the spectral band characteristics of the sub-pixel mapping sample points obtained in step 5.2, the non-spectral characteristics of the sub-pixel mapping sample points obtained in step 5.3, and the percentage of crop sub-pixels at the sub-pixel mapping sample points corresponding to each crop type obtained in step 5.4 to obtain the sub-pixel mapping sample set of the study area.

[0041] Step 6 as described above includes the following steps:

[0042] Step 6.1 construct a mapping model so that the spectral band characteristics and non-spectral characteristics of the sub-pixel mapping sample points are input into the mapping model to obtain the predicted crop sub-pixel percentage at the sub-pixel mapping sample points;

[0043] Step 6.2: Use ten-fold cross validation to divide the sub-pixel mapping sample set of the study area into a training set and a validation set. Use the training set to train the mapping model, and use the validation set to cross-validate the trained mapping model. For each crop type, select the trained mapping model with the best precision evaluation index in the cross validation as the corresponding crop mapping optimization model.

[0044] Step 6.3: Input the preprocessed study area time series image dataset of the target year into the crop mapping optimization model corresponding to each crop type to obtain the sub-pixel distribution map corresponding to the crop type. The sub-pixel distribution map shows the final predicted crop sub-pixel percentage at each sub-pixel mapping sample point.

[0045] As described above in steps 3.3 and 3.6, the field sample points are preferentially selected from the field sample points of the target year in step 2.1. If the number of field sample points of the target year is less than the field sample point threshold, field sample points of historical years are also included.

[0046] As described above in step 5.4, for each crop type, the percentage of crop sub-pixel at each sub-pixel mapping sample point is calculated in turn, specifically including the following steps:

[0047] For each sub-pixel mapping sample point corresponding to the pre-processed time series image dataset of the study area, the total number of reference pixels N and the number of reference pixels a containing the corresponding crop type are counted. The percentage P of crop sub-pixels at the sub-pixel mapping sample point is: P =a / N; the reference pixel is the pixel of the crop distribution reference map of the study sub-area in the target year.

[0048] In step 3.2 as described above, the pixel changes of each prepared amplified sample point from the pixel in the false color median composite image corresponding to the phenological critical period of the historical year to the pixel in the false color median composite image corresponding to the phenological critical period of the target year are compared one by one, specifically including the following steps:

[0049] The pixel values ​​of the near infrared band, red band, and green band of the pixels in the false color median composite image corresponding to the phenological key period of the historical year at the prepared amplification sample point are recorded as NIR1 pixel value, R1 pixel value, and G1 pixel value, respectively;

[0050] The pixel values ​​of the near infrared band, red band, and green band of the pixel in the false color median composite image corresponding to the phenological key period of the target year at the prepared amplification sample point are NIR2 pixel value, R2 pixel value, and G2 pixel value respectively;

[0051] Comparison thresholds are set for the near infrared band, red band, and green band, respectively, and are recorded as NIR_th comparison threshold, R_th comparison threshold, and G_th comparison threshold;

[0052] If |NIR1- NIR2|≤NIR_th, and |R1- R2|≤R_th, and |G1- G2|≤G_th, it is considered that the change of the corresponding prepared amplification sample point from the pixel in the false color median composite image corresponding to the phenological critical period of the historical year to the pixel in the false color median composite image corresponding to the phenological critical period of the target year is not obvious;

[0053] If |NIR1- NIR2|>NIR_th, or |R1- R2|>R_th, or |G1- G2|>G_th, it is considered that the corresponding preliminary amplification sample point corresponds to a significant change from the pixel in the false color median composite image corresponding to the phenological critical period of the historical year to the pixel in the false color median composite image corresponding to the phenological critical period of the target year.

[0054] In step 3.5 described above, judging whether the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e specifically includes the following steps:

[0055] For each crop type in each research sub-area, the median of the dynamic time warping distance DTW_cr_c&c and the median of the dynamic time warping distance DTW_cr_c&e are counted, a threshold interval is set, and the median of the dynamic time warping distance DTW_cr_c&c is included in the threshold interval. If the median of the dynamic time warping distance DTW_cr_c&e is within the threshold interval, it is considered that the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e, otherwise it is considered that the numerical distribution of the dynamic time warping distance DTW_cr_c&c is not close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e;

[0056] The step of setting the threshold range specifically includes the following steps:

[0057] Set the threshold percentage equal to A, A<50%, then the upper limit of the threshold interval is , the lower limit of the threshold interval is .

[0058] Compared with the prior art, the present invention has the following beneficial effects:

[0059] The present invention firstly expands the sample points based on the performance of different crops in the remote sensing data source and the similarity of the vegetation index distribution of the sample points. Secondly, by combining multi-source remote sensing data, the high spatial resolution advantage of medium- and high-spatial-resolution satellite images and the medium- and high-temporal resolution advantage of low-spatial-resolution satellite images are brought into play. In a small area, a high-precision crop distribution reference map is obtained by combining machine learning, and the sub-pixel crop recognition capability of medium- and high-resolution images is improved. At the same time, a sufficient number of uniformly distributed sub-pixel mapping sample sets of the study area are generated through the crop distribution reference map of a small-scale study sub-area, thereby improving the efficiency of obtaining the sub-pixel distribution map of crops in a large-scale study area.

[0060] The present invention adopts a soft classification method. By considering the spectral changes within the pixel and the spatial similarity, a regression model is established to obtain the percentages of different end members in each pixel, determine the proportions of different crops at the sub-pixel scale, fully mine data information to achieve accurate extraction of crops in large areas, effectively improve the problem of mixed pixels, and has great application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 It is a schematic flow chart of the method of the present invention.

[0062] Figure 2a It is a consistency assessment diagram of the amplified sample points and field-collected sample points corresponding to corn in the study sub-area.

[0063] Figure 2b It is a consistency assessment diagram of the amplified sample points and field-collected sample points corresponding to rice in the study sub-area.

[0064] Figure 2c It is a consistency assessment diagram of the amplified sample points and field-collected sample points corresponding to soybean in the study sub-area.

[0065] Among them, the red area represents the numerical distribution of the dynamic time warping distance calculated based on the time series data of the normalized vegetation index NDVI at two different field-collected sample points; the blue area represents the numerical distribution of the dynamic time warping distance calculated based on the time series data of the normalized vegetation index NDVI at a field-collected sample point and the time series data of the normalized vegetation index NDVI at an amplified sample point. DETAILED DESCRIPTION

[0066] In order to facilitate those skilled in the art to understand and implement the present invention, the present invention is further described in detail below with reference to examples. The implementation examples described herein are only used to illustrate and explain the present invention and are not intended to limit the present invention.

[0067] Example

[0068] A rapid crop mapping method for large areas integrating multi-source remote sensing data, the flow chart is as follows Figure 1 As shown, the following steps are included:

[0069] Step 1, obtain multi-source remote sensing satellite images and preprocess them, the multi-source remote sensing satellite images include original medium and high spatial resolution satellite images corresponding to the time span of the study sub-area and original low spatial resolution satellite images corresponding to the time span of the study area; the original medium and high spatial resolution satellite images are preprocessed to obtain a preprocessed study sub-area time series image dataset; the original low spatial resolution satellite images are preprocessed to obtain a preprocessed study area time series image dataset;

[0070] Among them, the study sub-area is included in the study area, the time span of the study sub-area and the time span of the study area are both included in the range of the target year, and the time span of the study sub-area is included in the range of the time span of the study area. The original medium-high spatial resolution satellite images include: one or more of the original Landsat 8 series satellite images, the original Sentinel-2 satellite images and the original domestic high-resolution series satellite images; the original low spatial resolution satellite images include: the original MODIS satellite images; the embodiment provides a method, which is applicable to the remaining satellite images, and when used, it is only limited that the spatial resolution of the original medium-high spatial resolution satellite images corresponding to the study sub-area is higher than the spatial resolution of the original low spatial resolution satellite images corresponding to the study area.

[0071] Acquire multi-source remote sensing satellite images and perform preprocessing, which includes the following steps:

[0072] Step 1.1 Select the study area and the study sub-area in the study area, and specify the time span of the study area and the time span of the study sub-area;

[0073] Obtain original medium- and high-spatial-resolution satellite images of each study sub-area corresponding to the time span of the study sub-area;

[0074] Obtain the original low spatial resolution satellite images of the study area corresponding to the time span of the study area.

[0075] In this embodiment, the study area includes Heilongjiang Province, Jilin Province and Liaoning Province; the study sub-area includes 8 counties and county-level cities in the study area, including 4 in Heilongjiang Province (Longjiang County of Qiqihar City, Nehe City of Qiqihar City, Yilan County of Harbin City, and Fujin City of Jiamusi City), 2 in Jilin Province (Fuyu City of Songyuan City, Qian'an County of Songyuan City), and 2 in Liaoning Province (Changtu County of Tieling City and Xinmin City of Shenyang City);

[0076] The time span of the study area and the study sub-area is from March 15 to November 15, 2020;

[0077] The original medium and high spatial resolution satellite images obtained include: original Landsat 8 satellite images and original Sentinel-2 satellite images. This embodiment uses both the original Landsat 8 satellite images with medium spatial resolution and the original Sentinel-2 satellite images with high spatial resolution to ensure the mapping effect; the original low spatial resolution satellite images obtained are original MODIS satellite images.

[0078] In this embodiment, the original Landsat 8 satellite image, the original Sentinel-2 satellite image, and the original MODIS satellite image are all obtained from the Google Earth Engine platform.

[0079] Step 1.2 Preprocess the original medium- and high-spatial-resolution satellite images and the original low-spatial-resolution satellite images, respectively, to obtain the preprocessed study sub-area time series image dataset and the preprocessed study area time series image dataset, respectively, which specifically includes the following steps:

[0080] The original low spatial resolution satellite images were preprocessed by cloud removal and median synthesis to obtain a time series image dataset of the study area after preprocessing.

[0081] For original medium and high spatial resolution satellite images,

[0082] When an original medium-high spatial resolution satellite image is obtained, the preprocessing process includes: cloud removal and median synthesis, and a time series image dataset of the research sub-area is obtained after preprocessing;

[0083] When a variety of original medium and high spatial resolution satellite images are obtained, the preprocessing process includes: declouding various original medium and high spatial resolution satellite images respectively to obtain the corresponding medium and high spatial resolution satellite images after declouding, performing reflectivity standardization on the medium and high spatial resolution satellite images after declouding to obtain medium and high spatial resolution satellite images after reflectivity standardization, and median synthesizing the medium and high spatial resolution satellite images after reflectivity standardization to obtain the time series image dataset of the research sub-area after preprocessing.

[0084] The cloud removal processing for different remote sensing satellite images is as follows:

[0085] For the original Landsat 8 satellite image, the Quality Assessment (QA) Band was used to remove clouds, shadows, and SLC-off strips (SLC-off strips are image strips after the abnormality of the airborne scan line corrector, where SLC is the abbreviation of the airborne scan line corrector Scan Lines Corrector) in the original Landsat 8 satellite image using the Google Earth Engine platform. The confidence level when removing clouds and shadows on the Google Earth Engine platform was 67-100%, and the Landsat 8 satellite image after cloud removal was obtained;

[0086] The original Sentinel-2 satellite image presents surface reflectance data. For the original Sentinel-2 satellite image, the Sentinel-2 cloud probability product provided by the Google Earth Engine platform is used to remove pixels with cloud probability higher than 65% in the original Sentinel-2 satellite image to obtain the de-clouded Sentinel-2 satellite image.

[0087] For the original MODIS satellite images, the Google Earth Engine platform is used in combination with the quality control band (StateQA) of the MODIS images. Bit operations are used to remove clouds, snow cover and other interfering pixels, and then a 30-day median synthesis is performed to obtain the preprocessed MODIS satellite images.

[0088] When multiple original medium and high spatial resolution satellite images are selected, reflectance standardization processing is performed on all cloud-removed medium and high spatial resolution satellite images to obtain medium and high spatial resolution satellite images after reflectance standardization processing, which specifically includes the following steps:

[0089] Taking a declouded medium-to-high spatial resolution satellite image as a reference, the spectral bands of other declouded medium-to-high spatial resolution satellite images are transformed through the conversion coefficient calculated by least squares regression to obtain the medium-to-high spatial resolution satellite image after spectral band conversion; the collection of the declouded medium-to-high spatial resolution satellite image as a reference and the medium-to-high spatial resolution satellite image after spectral band conversion is used as the medium-to-high spatial resolution satellite image after reflectance standardization.

[0090] In this embodiment, there are differences in band range and spectral reflectance between the Landsat 8 OLI remote sensing satellite sensor used to obtain the original Landsat 8 satellite image and the Sentinel-2 MSI remote sensing satellite sensor used to obtain the original Sentinel-2 satellite image. In order to use the Landsat 8 OLI remote sensing satellite sensor and the Sentinel-2 MSI remote sensing satellite sensor in a coordinated manner and obtain comparable time series images, the conversion coefficient calculated based on ordinary least squares regression is used to convert the spectral band of the de-clouded Sentinel-2 satellite image with reference to the de-clouded Landsat 8 satellite image to obtain the band-converted Sentinel-2 satellite image. Then, the collection of the de-clouded Landsat 8 satellite image and the band-converted Sentinel-2 satellite image is used as a medium-to-high spatial resolution satellite image after reflectance standardization.

[0091] In this embodiment, a collection of Landsat 8 satellite images after cloud removal and Sentinel-2 satellite images after band conversion in the study sub-area is subjected to median synthesis in a fixed time window to obtain a time series image dataset of the study sub-area after preprocessing; in this embodiment, the time window is selected as 20 days, that is, the time interval is 20 days.

[0092] Step 2: Obtain the field sample point data of each study sub-area in historical years and the field sample point data of the target year; obtain the crop distribution map of each study sub-area in historical years; obtain the land use data set of each study sub-area in the target year, and extract the cultivated land layer from the land use data set; obtain the pre-processed high spatial resolution satellite images corresponding to the phenological key period of each study sub-area in historical years and the pre-processed high spatial resolution satellite images corresponding to the phenological key period of the target year. The detailed operation steps are as follows:

[0093] Step 2.1 Obtain the field sample point data for each study sub-area in historical years and target years.

[0094] In this embodiment, the historical year is 2019, the target year is 2020, and the field-collected sample point data in 2019 and 2020 in the study sub-area are obtained; the field-collected sample point data is a point vector file, and the data attributes corresponding to each field-collected sample point include: the location information of the field-collected sample point (i.e., the latitude and longitude of the sample point) and the crop type at the field-collected sample point. The field-collected sample point data of the target year is used to determine whether to perform sample point amplification on the corresponding study sub-area in step 3.1, and then obtain the pre-amplified sample points; in step 3.3, the time series data of the normalized vegetation index NDVI corresponding to the field-collected sample points is obtained based on the field-collected sample points in the field-collected sample point data and the pre-processed study sub-area time series image data set obtained in step 1.2; in step 4, the field-collected sample point data of the target year participates in the construction of crop sample points.

[0095] Step 2.2 Obtain the crop distribution map of each study sub-area in historical years, which is used to obtain the crop types of the amplified sample points in historical years in step 3.1

[0096] In this embodiment, the crop distribution map of the research sub-area in 2019 is extracted from the crop distribution map of the Northeast region in 2019. The crop distribution map includes the crop types corresponding to each region. In this embodiment, the classification system of crop types in the crop distribution map includes: rice, corn, and soybeans, etc., and the distribution information of rice, corn, and soybeans is in the same crop distribution map;

[0097] Step 2.3 Obtain the land use dataset of each study sub-area in the target year, and extract the cultivated land layer of the study sub-area from the land use dataset within the study sub-area;

[0098] In this embodiment, the GlobeLand30 data of the study sub-area in 2020 is selected as the land use dataset of the study sub-area in the target year, and the cultivated land layer of the study sub-area is extracted from the GlobeLand30 data in the study sub-area; the GlobeLand30 data is also used to be combined with the original MODIS satellite image in step 5.3 to extract the non-spectral characteristics of the sub-pixel mapping sample points; the cultivated land layer is used to predict the crop distribution reference map in combination with the crop sample points in step 4.3.

[0099] Step 2.4: obtain and preprocess the original high spatial resolution satellite images corresponding to the phenological critical period of each study sub-area in the historical year and the target year, and obtain the preprocessed high spatial resolution satellite images corresponding to the phenological critical period of the study sub-area in the historical year and the preprocessed high spatial resolution satellite images corresponding to the phenological critical period of the target year; which are used to determine whether to retain the corresponding sample points for preliminary amplification in step 3.2;

[0100] In this embodiment, the critical phenological period of the historical year and the critical phenological period of the target year are selected as: July 15 to August 15, 2019, and July 15 to August 15, 2020, respectively.

[0101] The process of obtaining the original Sentinel-2 satellite images corresponding to the phenological key period of the historical year and the original Sentinel-2 satellite images corresponding to the phenological key period of the target year and preprocessing them in step 2.4 is the same as that in step 1 (i.e., downloading the original Sentinel-2 satellite images of the phenological key period of the historical year and the target year from the Google Earth Engine platform, and then performing cloud removal and median synthesis on the original Sentinel-2 satellite images corresponding to the phenological key period of the historical year and the target year), only the selected time range is different. In step 2.4, different original high spatial resolution satellite images from step 1 can also be used, such as domestic high-resolution series satellite images.

[0102] In this embodiment, in order to ensure the mapping accuracy, the Sentinel-2 satellite image is selected in step 2.4.

[0103] Step 3: For the research sub-areas where the number of field sample points corresponding to the crop type in the target year is less than the field sample point threshold, the sample points are expanded based on the pixel spectral similarity, and the dynamic time warping distance is calculated based on the time series data of the normalized vegetation index NDVI of the field sample points and the time series data of the normalized vegetation index NDVI of the expanded sample points, so as to verify and screen the expanded sample points and obtain reliable expanded sample points:

[0104] Step 3.1 For the study sub-areas where the number of field sample points corresponding to the crop type in the target year is less than the field sample point number threshold, randomly generate preliminary amplification sample points, and then superimpose the crop distribution map of each study sub-area in the historical year obtained in step 2.2 with the preliminary amplification sample points, and use the location information of the preliminary amplification sample points to extract the crop type of the preliminary amplification sample points in the historical year from the crop distribution map of the study sub-area in the historical year;

[0105] In this embodiment, the threshold of the number of field sample points is set to 150. In 2020, there are sufficient field sample points collected in Fujin City, Jiamusi City (the number of sample points for rice, corn, soybeans and other crops is greater than 150), and there is no need to expand the sample points. In the remaining 7 research sub-areas, 600 random points are randomly generated (i.e., rice, corn, soybeans, and other crops, a total of four categories, 150 randomly generated in each category) as pre-amplified sample points, where the minimum distance between two pre-amplified sample points is not less than the spatial resolution of the medium-high spatial resolution satellite sensor (i.e., the spatial resolution of the original medium-high spatial resolution satellite image). In this embodiment, the minimum distance between two pre-amplified sample points is set to the spatial resolution of the Landsat 8 OLI remote sensing satellite sensor, i.e., 30 meters. Then, using the location information of the pre-amplified sample points, the crop type of each pre-amplified sample point in the historical year is extracted from the crop distribution map of each research sub-area in the historical year obtained in step 2.2.

[0106] Step 3.2 Perform false color synthesis on the preprocessed high spatial resolution satellite images corresponding to the phenological critical period of each study sub-area in historical years and the preprocessed high spatial resolution satellite images corresponding to the phenological critical period of the target year obtained in step 2.4, respectively, to obtain the false color median composite images corresponding to the phenological critical period of each study sub-area in historical years and the false color median composite images corresponding to the phenological critical period of the target year respectively; compare the changes of each prepared amplification sample point from the pixels in the false color median composite image corresponding to the phenological critical period of the historical years to the pixels in the false color median composite image corresponding to the phenological critical period of the target year one by one; if the change is not obvious, the corresponding prepared amplification sample point is retained as the amplification sample point, and the crop type corresponding to the amplification sample point is retained (the crop type is obtained in step 3.1); if the change is obvious, the corresponding prepared amplification sample point is deleted.

[0107] Determine whether the change of each prepared amplified sample point from the pixel in the false color median composite image corresponding to the phenological key period of the historical year to the pixel in the false color median composite image corresponding to the phenological key period of the target year is obvious, specifically including the following steps:

[0108] The pixel values ​​of the near infrared band, red band, and green band of the pixels in the false color median composite image corresponding to the phenological key period of the historical year at the prepared amplification sample point are recorded as NIR1 pixel value, R1 pixel value, and G1 pixel value, respectively;

[0109] The pixel values ​​of the near infrared band, red band, and green band of the pixel in the false color median composite image corresponding to the phenological key period of the target year at the prepared amplification sample point are NIR2 pixel value, R2 pixel value, and G2 pixel value respectively;

[0110] Comparison thresholds are set for the near infrared band, red band, and green band, respectively, and are recorded as NIR_th comparison threshold, R_th comparison threshold, and G_th comparison threshold;

[0111] If |NIR1- NIR2|≤NIR_th, and |R1- R2|≤R_th, and |G1- G2|≤G_th, it is considered that the change of the corresponding prepared amplification sample point from the pixel in the false color median composite image corresponding to the phenological critical period of the historical year to the pixel in the false color median composite image corresponding to the phenological critical period of the target year is not obvious;

[0112] If |NIR1- NIR2|>NIR_th, or |R1- R2|>R_th, or |G1- G2|>G_th, it is considered that the corresponding preliminary amplification sample point corresponds to a significant change from the pixel in the false color median composite image corresponding to the phenological critical period of the historical year to the pixel in the false color median composite image corresponding to the phenological critical period of the target year.

[0113] In this embodiment, the pre-processed Sentinel-2 satellite image corresponding to the phenological critical period in 2019 obtained in step 2.4 and the pre-processed Sentinel-2 satellite image corresponding to the phenological critical period in 2020 are respectively subjected to false color synthesis to obtain the false color median synthetic image corresponding to the phenological critical period in 2019 and the false color median synthetic image corresponding to the phenological critical period in 2020, and all the prepared amplified sample points are superimposed; the change of each prepared amplified sample point from the pixel in the false color median synthetic image corresponding to the phenological critical period in 2019 to the pixel in the false color median synthetic image corresponding to the phenological critical period in 2020 is compared one by one, and the change of the false color median synthetic image corresponding to the prepared amplified sample point represents the change of the reflectivity of the pixel where the prepared amplified sample point is located. If there is no obvious change, the corresponding prepared amplified sample point is retained as the amplified sample point, and the crop type corresponding to the amplified sample point is retained; if there is an obvious change, the corresponding prepared amplified sample point is deleted; finally, the amplified sample point in each research sub-area is obtained. The comparison threshold of near infrared is NIR_th=0.15, the comparison threshold of red light band is R_th=0.015, and the comparison threshold of green light band is G_th=0.015.

[0114] Step 3.3 Based on the preprocessed time series image dataset of the study sub-area obtained in step 1.2, for each crop type in each study sub-area, obtain the time series data of the normalized vegetation index NDVI corresponding to the field collection sample points and the time series data of the normalized vegetation index NDVI corresponding to the amplified sample points respectively; in step 3.3, the field collection sample points are preferentially selected from the field collection sample points of the target year. If the number of field collection sample points in the target year is less than the field sample point number threshold, the field collection sample points of the historical years are also included; the time series data of the normalized vegetation index NDVI includes the NDVI value at each time series point;

[0115] Set the upper quantile and lower quantile limits; used for further screening of amplification sample points in subsequent steps;

[0116] Step 3.4: For each crop type in each study sub-area, the values ​​of the upper quantile of the NDVI values ​​at all field-collected sample points are counted, and recorded as the upper quantile of the NDVI of the crop type; and the values ​​of the lower quantile of the NDVI values ​​at all field-collected sample points are counted, and recorded as the lower quantile of the NDVI of the crop type (in this embodiment, the upper quantile is the 90th quantile, and the lower quantile is the 10th quantile);

[0117] Step 3.5: For each crop type in each study sub-area, for each current amplified sample point, compare the NDVI value of the amplified sample point at each time point with the crop type NDVI upper limit value and the crop type NDVI lower limit value. If the NDVI value of the amplified sample point at any time point is lower than the crop type NDVI lower limit value or higher than the crop type NDVI upper limit value, the amplified sample point will be eliminated.

[0118] In this embodiment, after steps 3.2- 3.5 After two-step screening: the number of expanded sample points for rice, corn, soybean and other crops in Longjiang County were 124, 133, 96 and 124 respectively; the number of expanded sample points for rice, corn, soybean and other crops in Nehe County were 140, 126, 129 and 136 respectively; the number of expanded sample points for rice, corn, soybean and other crops in Yilan County were 132, 136, 103 and 126 respectively; the number of expanded sample points for rice, corn, soybean and other crops in Qian'an County were 91, 142, 94 and 114 respectively; the number of expanded sample points for rice, corn, soybean and other crops in Fuyu County were 118, 144, 135 and 106 respectively; the number of expanded sample points for rice, corn, soybean and other crops in Changtu County were 96, 141, 118 and 115 respectively; the number of expanded sample points for rice, corn, soybean and other crops in Xinmin County were 129, 134, 107 and 129 respectively.

[0119] Step 3.6: Verify the reliability of the amplified sample points in the study sub-area:

[0120] In each research sub-area: for each crop type, the corresponding dynamic time warping distance is calculated based on the time series data of the normalized vegetation index NDVI at two different field-collected sample points, which is recorded as the dynamic time warping distance DTW_cr_c&c; corresponding to the pairwise combination of all field-collected sample points of the same crop type, all dynamic time warping distances DTW_cr_c&c corresponding to the crop type are obtained; at the same time, for each crop type, the corresponding dynamic time warping distance is calculated based on the time series data of the normalized vegetation index NDVI at a field-collected sample point and the time series data of the normalized vegetation index NDVI at an amplified sample point, which is recorded as the dynamic time warping distance DTW_cr_c&e, corresponding to the pairwise combination of all field-collected sample points and amplified sample points of the same crop type, all dynamic time warping distances DTW_cr_c&e corresponding to the crop type are obtained;

[0121] For each crop type in each study sub-area: determine whether the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e. If so, it is considered that the amplified sample points in the corresponding crop type in the corresponding study sub-area are reliable, and then proceed to the subsequent step 4; if not, reduce the upper quantile and increase the lower quantile, and return to step 3.4.

[0122] Determining whether the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e specifically includes the following steps:

[0123] For each crop type in each study sub-area, the median (i.e., the median) of the dynamic time warping distance DTW_cr_c&c and the median of the dynamic time warping distance DTW_cr_c&e are counted, a threshold interval is set, and the median of the dynamic time warping distance DTW_cr_c&c is included in the threshold interval. If the median of the dynamic time warping distance DTW_cr_c&e is within the threshold interval, it is considered that the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e, otherwise it is considered that the numerical distribution of the dynamic time warping distance DTW_cr_c&c is not close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e.

[0124] Setting the threshold range includes the following steps:

[0125] If the threshold percentage is set equal to A (A<50%), the upper limit of the threshold interval is , the lower limit of the threshold interval is .

[0126] In this embodiment, in order to facilitate intuitive comparison of the numerical distribution of the dynamic time warping distance DTW_cr_c&c and the numerical distribution of the dynamic time warping distance DTW_cr_c&e, in each research sub-area, the corresponding field sample vs field sample DTW distribution curve is drawn based on all the dynamic time warping distances DTW_cr_c&c corresponding to the crop type; the field sample vs field sample DTW distribution curve is used as a reference standard for the reliability of the amplified sample points; the corresponding field sample vs amplified sample DTW distribution curve is drawn based on all the dynamic time warping distances DTW_cr_c&e corresponding to the crop type; and the reliability of the amplified sample points is judged by comparing the field sample vs field sample DTW distribution curve with the field sample vs amplified sample DTW distribution curve.

[0127] In the dynamic time warping distance DTW_cr_c&c and the dynamic time warping distance DTW_cr_c&e, cr represents the crop type, c&c represents "field sample vs. field sample", and c&e represents "field sample vs. amplified sample". In this embodiment, the crop type is rice, corn or soybean. Taking rice as an example, cr is "rice":

[0128] The dynamic time warping distance DTW_rice_c&c is calculated by using the time series data of the normalized vegetation index NDVI at a field-collected sample point of rice and the time series data of the normalized vegetation index NDVI at another field-collected sample point of rice;

[0129] The dynamic time warping distance DTW_rice_c&c is to calculate the dynamic time warping distance using the time series data of the normalized vegetation index NDVI at a field sample point of rice and the time series data of the normalized vegetation index NDVI at an amplified sample point of rice. Calculating the dynamic time warping distance is an existing technology.

[0130] Similar to step 3.3, in step 3.6, the field-collected sample points are preferentially selected from the field-collected sample points of the target year. If the number of field-collected sample points of the target year is less than the field-collected sample point threshold, the field-collected sample points of the historical years are also included.

[0131] Figure 2a~Figure 2c In order to obtain the consistency evaluation result diagram of the amplified sample points and the field collected sample points in the study sub-area (i.e., the reliability evaluation diagram of the amplified sample points in the study sub-area), Figure 2a~Figure 2c The long dashed line in the middle represents the median, and the two short dashed lines above and below represent 1.5 times the interquartile range. Figure 2a For example, if there are m field sample points for rice, the corresponding The dynamic time warping distance DTW_Rice_c&c, the median and interquartile range refer to the median and interquartile range of the entire dynamic time warping distance DTW_Rice_c&c value set. The vertical axis in the figure represents the dynamic time warping distance value; for the horizontal axis, for each research sub-area, the field sample vs field sample DTW distribution curve and the field sample vs amplified sample DTW distribution curve are located on both sides of the corresponding vertical line, and the distance from the vertical line represents the number of sample points with the corresponding dynamic time warping distance value (the red area corresponds to the field sample vs field sample DTW distribution curve and the area surrounded by the vertical line, and the blue area corresponds to the field sample vs amplified sample DTW distribution curve and the area surrounded by the vertical line). Figure 2a~Figure 2c Each crop type in the study sub-area corresponds to Figure 2a~Figure 2cThe lengths of the same type of dashed lines on the left and right sides are the same. The closer the values ​​of the three dashed lines are, the more reliable the amplification sample point set of the corresponding research sub-area is. In this embodiment, 20% is used as the threshold percentage, that is, if the median of the field sample vs. amplification sample DTW distribution curve is within 20% (i.e., 80%-120%) of the median of the field sample vs. field sample DTW distribution curve, then the amplification sample point set of the corresponding research sub-area is considered to be reliable. Figure 2a~Figure 2c It can be seen that the median of the dynamic time warping distance DTW_cr_c&c value distribution is similar to the median of the dynamic time warping distance DTW _cr_c&e value distribution, and the amplified sample points are highly consistent with the field collected sample points, indicating that the amplified sample points are reliable.

[0132] Step 4: Based on the field sample points of the target year and the amplified sample points obtained in step 3, combined with the pre-processed time series image dataset of the study sub-area, obtain the reference map of crop distribution in the study sub-area in the target year. The detailed operation steps are as follows:

[0133] Step 4.1 For each crop type in each study sub-area, combine the field sample points of the target year obtained in step 2.1 and the amplified sample points obtained in step 3 as the crop sample points of the corresponding crop type in each study sub-area.

[0134] In this embodiment, the field sample points collected in 2020 in step 2 and the amplified sample points obtained in step 3 are combined as the crop sample points of the eight study sub-areas in 2020;

[0135] Step 4.2 calculates the time series characteristics of the vegetation index at the pixel scale based on the preprocessed study sub-area time series image dataset obtained in step 1.2, that is, based on each pixel in the preprocessed study sub-area time series image dataset at each time point in the time span of the study sub-area, the following vegetation indices are calculated respectively. This embodiment selects: Normalized Difference Vegetation Index NDVI, Enhanced Vegetation Index EVI, Surface Water Index LSWI, Normalized Thermal Difference Index NDTI and Normalized Difference Decay Index NDSVI.

[0136] Step 4.3 For each crop type in each study sub-area, the crop sample points obtained in step 4.1 and the pixel-scale vegetation index time series characteristics obtained in step 4.2 are combined, and a pixel-by-pixel prediction is performed on the cultivated land layer obtained in step 2.3 based on the random forest model to obtain a reference map of crop distribution in the study sub-area in the target year. The pixels of the reference map of crop distribution in the study sub-area in the target year are used as reference pixels.

[0137] Step 5: randomly generate multiple sub-pixel mapping sample points, and construct a sub-pixel mapping sample set for the study area based on the sub-pixel mapping sample points, the pre-processed time series image dataset of the study area, the land use dataset, and the crop distribution reference map. Specifically, the following steps are included:

[0138] Step 5.1, randomly generate 2500 sub-pixel mapping sample points in each study sub-area (a total of 20,000 sub-pixel mapping sample points are randomly generated in 8 study sub-areas), and the spacing between sub-pixel mapping sample points is not less than the spatial resolution of the low spatial resolution sensor (i.e., the original low spatial resolution satellite image), so in this embodiment, the spacing between each randomly generated sub-pixel mapping sample point is greater than or equal to 500m;

[0139] Step 5.2: Extracting spectral features of sub-pixel mapping sample points: Using the location information of the sub-pixel mapping sample points, extracting spectral band features of the corresponding sub-pixel mapping sample points from the pre-processed study area time series image dataset (in this embodiment, the pre-processed study area time series image dataset is the pre-processed MODIS satellite image), a total of 7 bands;

[0140] Step 5.3 Extract non-spectral features of sub-pixel mapping sample points:

[0141] The land use dataset of the target year in step 2.3 and the preprocessed time series image dataset of the study area are superimposed. The sub-pixels of the preprocessed time series image dataset of the study area are used to count the percentage of natural vegetation sub-pixels and non-vegetation sub-pixels in the sub-pixels where the sub-pixel mapping sample points are located. The percentage of natural vegetation sub-pixels and the percentage of non-vegetation sub-pixels are used as non-spectral features.

[0142] Considering that the percentage of crop sub-pixels is also affected by the proportion of other non-cultivated vegetation and non-vegetation in the sub-pixels, the non-spectral features select the percentage of natural vegetation (grassland, wetland, woodland, shrub) sub-pixels and the percentage of non-vegetation (water body, artificial surface, bare land, tundra, permanent snow) sub-pixels. In this embodiment, GlobeLand30 and the pre-processed MODIS satellite image are superimposed, and the percentage of natural vegetation sub-pixels and non-vegetation sub-pixels in the sub-pixels (MODIS pixels in this embodiment) where the sub-pixel mapping sample point is located are counted through the sub-pixels of the pre-processed MODIS satellite image, and used as non-spectral features for subsequent model training;

[0143] Step 5.4 Extract the percentage of sub-pixel crops at the sub-pixel mapping sample point:

[0144] The crop distribution reference map of the study sub-area in the target year obtained in step 4 is superimposed with the preprocessed study area time series image dataset. For each crop type (such as corn, rice and soybean in this embodiment), the percentage of crop sub-pixel at each sub-pixel mapping sample point is calculated in turn, specifically including the following steps:

[0145] Overlay the crop distribution reference map of the study sub-area in the target year obtained in step 4.3 with the pre-processed time series image dataset of the study area;

[0146] For each sub-pixel mapping sample point corresponding to the pre-processed time series image dataset of the study area, the total number of reference pixels N and the number of reference pixels a containing the corresponding crop type (e.g., soybean) are counted. The percentage P of crop sub-pixels at the sub-pixel mapping sample point is: P =a / N. The reference pixel is the pixel of the reference map of crop distribution in the study sub-area in the target year.

[0147] Based on this method, the sub-pixel percentages of crops for rice, corn, and soybean can be calculated respectively.

[0148] Step 5.5: Combine the spectral band characteristics of the sub-pixel mapping sample points obtained in step 5.2, the non-spectral characteristics of the sub-pixel mapping sample points obtained in step 5.3, and the percentage of crop sub-pixels at the sub-pixel mapping sample points corresponding to each crop type obtained in step 5.4 to obtain the sub-pixel mapping sample set of the study area.

[0149] In this embodiment, the number of sub-pixel mapping sample points is 20,000, each sub-pixel mapping sample point includes 7 spectral features and 2 non-spectral features, and the number of time points is 8, that is, time points with equal time intervals selected from the time span (March 15, 2020 to November 15, 2020) selected in step 1. Since the MODIS satellite image is synthesized for 30 days after cloud removal, the time points here change from 30 (the time interval from March 15, 2020 to November 15 is 8) to 8 (the time interval from March 15, 2020 to November 15 is 30).

[0150] Step 6: Construct a mapping model, and use the sub-pixel mapping sample set of the study area to train and verify the screening mapping model to obtain the crop mapping optimization model corresponding to each crop type. Input the pre-processed study area time series image dataset into the crop mapping optimization model to obtain the sub-pixel distribution map corresponding to each crop type. The specific steps include the following:

[0151] Step 6.1 constructs a mapping model so that the spectral band characteristics and non-spectral characteristics of the sub-pixel mapping sample points are input into the mapping model to obtain the predicted crop sub-pixel percentage at the sub-pixel mapping sample points: In this embodiment, the mapping model uses a Bi-LSTM regression model, and the input data of the Bi-LSTM regression model are the spectral band characteristics and non-spectral characteristics of the sub-pixel mapping sample points in the sub-pixel mapping sample set of the study area obtained in step 5.5. The Bi-LSTM regression model outputs the predicted crop sub-pixel percentage at the sub-pixel mapping sample points. The training label is the crop sub-pixel percentage at the sub-pixel mapping sample points corresponding to each crop type in step 5.4. The sub-pixel mapping sample set of the study area has a total of 9 features and 8 time points.

[0152] Step 6.2 Use ten-fold cross validation to divide the sub-pixel mapping sample set of the study area into training set and validation set. Use the training set to train the mapping model, and use the validation set to cross-validate the trained mapping model. For each crop type, select the trained mapping model with the best precision evaluation index in the cross-validation as the corresponding crop mapping preferred model. Since the sub-pixel mapping sample set of the study area contains the percentage of crop sub-pixels at the sub-pixel mapping sample points of multiple crop types, cross-validation is performed on the validation set each time, and the corresponding precision evaluation index is calculated for each crop type.

[0153] In this embodiment, the number of sub-pixel mapping sample points in the training set is 18,000 (20,000 0.9), the number of sub-pixel mapping sample points in the validation set is 2000 (20000 0.1), train the mapping model through the training set, and save the mapping model and accuracy evaluation results after each training; the batch size of each batch of training sample points is selected as 128, the learning rate is 0.001, and the training round epoch is 50 rounds; the determination coefficient R² and the root mean square error RMSE are used as the accuracy evaluation index of the validation set cross-validation.

[0154] In this embodiment, a corn mapping optimization model, a rice mapping optimization model and a soybean mapping optimization model are obtained respectively.

[0155] Step 6.3 Input the preprocessed study area time series image dataset of the target year into the crop mapping optimization model corresponding to each crop type to obtain the sub-pixel distribution map corresponding to the crop type. The sub-pixel distribution map shows the final predicted crop sub-pixel percentage at each sub-pixel mapping sample point.

[0156] In this embodiment, the rice mapping optimization model is loaded, and the MODIS satellite image preprocessed in the study area is input to obtain the sub-pixel distribution map of rice; the corn mapping optimization model is loaded, and the MODIS satellite image preprocessed in the study area is input to obtain the sub-pixel distribution map of corn; the soybean mapping optimization model is loaded, and the MODIS satellite image preprocessed in the study area is input to obtain the sub-pixel distribution map of soybean.

[0157] The present invention is not limited to the above-mentioned implementation modes. The above-mentioned embodiments are only descriptions of the preferred embodiments of the present invention, and do not limit the concept of the present invention. The implementation schemes in the above-mentioned embodiments can be further combined or replaced. Various changes and improvements made to the technical solutions of the present invention by those skilled in the art all belong to the protection scope of the present invention.

Claims

1. A method for rapid mapping of large-area crops by integrating multi-source remote sensing data, characterized in that: The steps include: Step 1: Obtain the original low spatial resolution satellite images of the study area corresponding to the time span of the study area and the original medium and high spatial resolution satellite images of the study sub-area corresponding to the time span of the study sub-area, and preprocess them respectively to obtain the preprocessed study area time series image dataset and the preprocessed study sub-area time series image dataset; Step 2: Obtain the field sample point data collected in historical years and the field sample point data collected in the target year for each study sub-area; Obtain crop distribution maps of each research sub-area in historical years; obtain land use datasets of each research sub-area in the target year, and extract cultivated land layers from the land use datasets; obtain pre-processed high spatial resolution satellite images corresponding to the phenological key periods of each research sub-area in historical years and pre-processed high spatial resolution satellite images corresponding to the phenological key periods of the target year; Step 3: For the research sub-areas where the number of field sample points corresponding to the crop type in the target year is less than the field sample point threshold, sample points are amplified based on pixel spectral similarity, and the reliability of the amplified sample points is verified and screened to obtain reliable amplified sample points; Step 4: Based on the field sample points of the target year and the amplified sample points obtained in step 3, combined with the pre-processed time series image dataset of the study sub-area, obtain a reference map of crop distribution in the study sub-area in the target year; Step 5: randomly generate multiple sub-pixel mapping sample points, and construct a sub-pixel mapping sample set of the study area based on the sub-pixel mapping sample points, the pre-processed time series image dataset of the study area, the land use dataset, and the crop distribution reference map; Step 6: Construct a mapping model, and use the sub-pixel mapping sample set of the study area to train and verify the screening mapping model to obtain the crop mapping optimization model corresponding to each crop type. Input the preprocessed study area time series image dataset into the crop mapping optimization model to obtain the sub-pixel distribution map corresponding to each crop type.

2. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 1 is characterized in that: In step 1, pretreatment is performed, which specifically includes the following steps: The original low spatial resolution satellite images were declouded and median synthesized to obtain a pre-processed time series image dataset of the study area. For original medium and high spatial resolution satellite images, When one kind of original medium-high spatial resolution satellite image is obtained, the original medium-high spatial resolution satellite image is de-clouded and median-synthesized to obtain a time series image dataset of the research sub-area after preprocessing; When a variety of original medium and high spatial resolution satellite images are obtained, the various original medium and high spatial resolution satellite images are declouded to obtain the corresponding medium and high spatial resolution satellite images after declouding, the medium and high spatial resolution satellite images after declouding are processed by reflectivity standardization to obtain medium and high spatial resolution satellite images after reflectivity standardization, and the medium and high spatial resolution satellite images after reflectivity standardization are synthesized by median value to obtain the time series image dataset of the research sub-area after preprocessing; Among them, the reflectivity standardization processing is performed on the cloud-removed medium- and high-spatial-resolution satellite images to obtain the medium- and high-spatial-resolution satellite images after reflectivity standardization processing, which specifically includes the following steps: Taking a declouded medium-to-high spatial resolution satellite image as a reference, the spectral bands of other declouded medium-to-high spatial resolution satellite images are transformed through the conversion coefficient calculated by least squares regression to obtain the medium-to-high spatial resolution satellite image after spectral band conversion; the collection of the declouded medium-to-high spatial resolution satellite image as a reference and the medium-to-high spatial resolution satellite image after spectral band conversion is used as the medium-to-high spatial resolution satellite image after reflectance standardization.

3. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 1 is characterized in that: The step 3 comprises the following steps: Step 3.1 For the study sub-areas where the number of field sample points corresponding to the crop type in the target year is less than the field sample point threshold, randomly generate preliminary amplification sample points, and the minimum distance between the preliminary amplification sample points is not less than the spatial resolution of the original medium- and high-spatial-resolution satellite image; Overlay the crop distribution map of each research sub-area in historical years obtained in step 2 with the prepared amplified sample points to obtain the crop types of the prepared amplified sample points in historical years; Step 3.2: Perform false color synthesis on the pre-processed high spatial resolution satellite images corresponding to the phenological critical period of each study sub-area in the historical years and the pre-processed high spatial resolution satellite images corresponding to the phenological critical period of the target year, respectively, to obtain the false color median composite images corresponding to the phenological critical period of each study sub-area in the historical years and the false color median composite images corresponding to the phenological critical period of the target year; compare the changes of each prepared amplification sample point from the pixels in the false color median composite images corresponding to the phenological critical period of the historical years to the pixels in the false color median composite images corresponding to the phenological critical period of the target year; if the changes are not obvious, the corresponding prepared amplification sample point is retained as the amplification sample point, and the crop type corresponding to the amplification sample point is retained; if the changes are obvious, the corresponding prepared amplification sample point is deleted; Step 3.3: Based on the preprocessed time series image dataset of the study sub-area obtained in step 1, for each crop type in each study sub-area, the time series data of the normalized vegetation index NDVI corresponding to the field collection sample points and the time series data of the normalized vegetation index NDVI corresponding to the amplified sample points are obtained respectively, and the time series data of the normalized vegetation index NDVI includes the NDVI value at each time series point; Set the upper and lower quantiles; Step 3.4: For each crop type in each study sub-area, count the values ​​of the NDVI values ​​at the upper quantile of all field-collected sample points, and record them as the NDVI upper limit of the crop type; and count the values ​​of the NDVI values ​​at the lower quantile of all field-collected sample points, and record them as the NDVI lower limit of the crop type; Step 3.5 For each crop type in each study sub-area, for each current amplified sample point, compare the NDVI value of the amplified sample point at each time point with the crop type NDVI upper limit value and the crop type NDVI lower limit value. If the NDVI value of the amplified sample point at any time point is lower than the crop type NDVI lower limit value or higher than the crop type NDVI upper limit value, then remove the amplified sample point; Step 3.6 In each study sub-area: for each crop type, the corresponding dynamic time warping distance is calculated based on the time series data of the normalized vegetation index NDVI at two different field-collected sample points, which is recorded as the dynamic time warping distance DTW_cr_c&c; corresponding to the pairwise combination of all field-collected sample points of the same crop type, all dynamic time warping distances DTW_cr_c&c corresponding to the crop type are obtained; at the same time, for each crop type, the corresponding dynamic time warping distance is calculated based on the time series data of the normalized vegetation index NDVI at a field-collected sample point and the time series data of the normalized vegetation index NDVI at an amplified sample point, which is recorded as the dynamic time warping distance DTW_cr_c&e, corresponding to the pairwise combination of all field-collected sample points and amplified sample points of the same crop type, all dynamic time warping distances DTW_cr_c&e corresponding to the crop type are obtained; For each crop type in each study sub-area: determine whether the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e. If so, it is considered that the amplified sample points in the corresponding crop type in the corresponding study sub-area are reliable, and then proceed to the subsequent step 4; if not, reduce the upper quantile and increase the lower quantile, and return to step 3.

4.

4. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 3 is characterized in that: The step 4 comprises the following steps: Step 4.1 For each crop type in each study sub-area, combine the field collection sample points and the amplification sample points in the target year to obtain the crop sample points of the corresponding crop type in each study sub-area; Step 4.2: Calculate the time series characteristics of vegetation indices at the pixel scale based on the preprocessed time series image dataset of the study sub-area. Vegetation indices include: Normalized Difference Vegetation Index NDVI, Enhanced Vegetation Index EVI, Land Surface Water Index LSWI, Normalized Thermal Difference Index NDTI, and Normalized Difference Decay Index NDSVI. Step 4.3: For each crop type in each study sub-area, the time series characteristics of the vegetation index at the crop sample point and pixel scale are combined, and a pixel-by-pixel prediction is performed on the cultivated land layer obtained in step 2 based on the random forest model to obtain a reference map of crop distribution in the study sub-area in the target year.

5. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 4 is characterized in that: The step 5 comprises the following steps: Step 5.1 For each study sub-area, randomly generate n sub-pixel mapping sample points, and the spacing between sub-pixel mapping sample points should not be less than the spatial resolution of the original low spatial resolution satellite image; Step 5.2 uses the location information of the sub-pixel mapping sample points to extract the spectral band characteristics of the corresponding sub-pixel mapping sample points from the pre-processed time series image dataset of the study area; Step 5.3 Overlay the land use dataset of the target year in step 2 and the preprocessed time series image dataset of the study area. Use the sub-pixels of the preprocessed time series image dataset of the study area to count the percentage of natural vegetation sub-pixels and non-vegetation sub-pixels in the sub-pixels where the sub-pixel mapping sample points are located. The percentage of natural vegetation sub-pixels and the percentage of non-vegetation sub-pixels are used as non-spectral features. Step 5.4 Overlay the crop distribution reference map of the study sub-area in the target year obtained in step 4 with the pre-processed study area time series image dataset, and for each crop type, calculate the percentage of crop sub-pixel at each sub-pixel mapping sample point in turn; Step 5.5: Combine the spectral band characteristics of the sub-pixel mapping sample points obtained in step 5.2, the non-spectral characteristics of the sub-pixel mapping sample points obtained in step 5.3, and the percentage of crop sub-pixels at the sub-pixel mapping sample points corresponding to each crop type obtained in step 5.4 to obtain the sub-pixel mapping sample set of the study area.

6. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 5 is characterized in that: The step 6 comprises the following steps: Step 6.1 construct a mapping model so that the spectral band characteristics and non-spectral characteristics of the sub-pixel mapping sample points are input into the mapping model to obtain the predicted crop sub-pixel percentage at the sub-pixel mapping sample points; Step 6.2: Use ten-fold cross validation to divide the sub-pixel mapping sample set of the study area into a training set and a validation set. Use the training set to train the mapping model, and use the validation set to cross-validate the trained mapping model. For each crop type, select the trained mapping model with the best precision evaluation index in the cross validation as the corresponding crop mapping optimization model. Step 6.3: Input the preprocessed study area time series image dataset of the target year into the crop mapping optimization model corresponding to each crop type to obtain the sub-pixel distribution map corresponding to the crop type. The sub-pixel distribution map shows the final predicted crop sub-pixel percentage at each sub-pixel mapping sample point.

7. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 3 is characterized in that: In step 3.3 and step 3.6, the field sample points are preferably selected from the field sample points of the target year in step 2.

1. If the number of field sample points of the target year is less than the field sample point threshold, field sample points of historical years are also included.

8. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 5 is characterized in that: In step 5.4, for each crop type, the percentage of crop sub-pixel at each sub-pixel mapping sample point is calculated in turn, specifically including the following steps: For each sub-pixel mapping sample point corresponding to the preprocessed time-series image dataset of the study area, the total number of reference pixels N and the number of reference pixels a containing the corresponding crop type are counted. The percentage P of crop sub-pixels at the sub-pixel mapping sample point is: P=a / N; where the reference pixel is the pixel of the crop distribution reference map of the study sub-area in the target year.

9. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 3 is characterized in that: In step 3.2, the pixel changes of each prepared amplified sample point from the pixel in the false color median composite image corresponding to the phenological key period of the historical year to the pixel in the false color median composite image corresponding to the phenological key period of the target year are compared one by one, specifically including the following steps: The pixel values ​​of the near infrared band, red band, and green band of the pixels in the false color median composite image corresponding to the phenological key period of the historical year at the prepared amplification sample point are recorded as NIR1 pixel value, R1 pixel value, and G1 pixel value, respectively; The pixel values ​​of the near infrared band, red band, and green band of the pixel in the false color median composite image corresponding to the phenological key period of the target year at the prepared amplification sample point are NIR2 pixel value, R2 pixel value, and G2 pixel value respectively; For the near infrared band, red light band and green light band, comparison thresholds are set respectively, which are recorded as NIR_th comparison threshold, R_th comparison threshold and G_th comparison threshold; If |NIR1- NIR2|≤NIR_th, and |R1- R2|≤R_th, and |G1- G2|≤G_th, it is considered that the change of the corresponding prepared amplification sample point from the pixel in the false color median composite image corresponding to the phenological critical period of the historical year to the pixel in the false color median composite image corresponding to the phenological critical period of the target year is not obvious; If |NIR1- NIR2|>NIR_th, or |R1- R2|>R_th, or |G1- G2|>G_th, it is considered that the corresponding preliminary amplification sample point corresponds to a significant change from the pixel in the false color median composite image corresponding to the phenological critical period of the historical year to the pixel in the false color median composite image corresponding to the phenological critical period of the target year.

10. The large-area crop rapid mapping method integrating multi-source remote sensing data according to claim 3 is characterized in that: In the step 3.5, judging whether the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e specifically includes the following steps: For each crop type in each research sub-area, the median of the dynamic time warping distance DTW_cr_c&c and the median of the dynamic time warping distance DTW_cr_c&e are counted, a threshold interval is set, and the median of the dynamic time warping distance DTW_cr_c&c is included in the threshold interval. If the median of the dynamic time warping distance DTW_cr_c&e is within the threshold interval, it is considered that the numerical distribution of the dynamic time warping distance DTW_cr_c&c is close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e, otherwise it is considered that the numerical distribution of the dynamic time warping distance DTW_cr_c&c is not close to the numerical distribution of the dynamic time warping distance DTW_cr_c&e; The step of setting the threshold range specifically includes the following steps: Set the threshold percentage equal to A, A<50%, then the upper limit of the threshold interval is , the lower limit of the threshold interval is .

Images