CN110728446A - County scale crop yield estimation method based on CNN-LSTM - Google Patents

Abstract

The invention relates to the technical field of remote sensing image information extraction, in particular to a county-level crop yield estimation method based on CNN-LSTM, which comprises the following steps: acquiring and processing S1 data, superimposing and filtering S2 data, and S3 acquiring county Level-scale feature tensor data, S4 builds and trains a CNN-LSTM model, and S5 applies the CNN-LSTM model trained in S4 to estimate the yield of target crops. The method of the invention is based on the remote sensing data reflecting the growth state of the crops and the environmental data affecting the growth of the crops, and extracts more features of the crops in the county area through the method of histogram statistics and converts them into tensors, so as to train the extracted CNN-LSTM model, Effectively improve the small-scale crop yield estimation accuracy.

Description

A county-level crop yield estimation method based on CNN-LSTM

technical field

The invention relates to the technical field of remote sensing image information extraction, in particular to a county-level crop yield estimation method based on CNN-LSTM.

Background technique

Crop yield is the most important indicator of agriculture and has many connections with human society. Yield forecasting is one of the most challenging tasks in precision agriculture and has important implications for yield mapping, crop market planning, crop insurance and harvest management.

Remote sensing technology has been widely used in crop yield estimation. Various relevant information can be extracted from remote sensing data to assist in yield estimation. In particular, various vegetation indices (VI), such as the normalized vegetation index (NDVI), have been widely used. Other indices such as Green Leaf Area Index (GLAI), Crop Water Stress Index (CWSI), Normalized Water Index (NDWI), Green Vegetation Index (GVI), Soil Adjusted Vegetation Index (SAVI) etc. have also been used for crop yield prediction . In addition, meteorological variables such as precipitation, air temperature and some soil condition data, including soil moisture, temperature and quality are often used in yield prediction as indicators of crop growth environment. Such feature extraction methods usually rely on computing the regional mean, but ignore detailed features.

Based on remote sensing data, there are two main methods of crop yield prediction: crop simulation and empirical statistical models. Although crop simulation models accurately simulate the physical process of crop growth, these models can hardly be applied on a large spatiotemporal scale due to insufficient data. In contrast, empirical statistical models are simple and require less input data, and thus have been widely used as a general-purpose alternative to process-based models. Machine learning algorithms including Support Vector Machines (SVMs), Decision Trees (DTs), Multilayer Perceptions (MLPs), and Restricted Boltzmann Machines (RBMs) can provide alternatives to traditional regression methods and overcome many limit. Furthermore, Artificial Neural Networks (ANNs) are also considered as alternative models. The traditional artificial neural network, the multi-layer perceptron model, has been successfully applied to crop yield estimation of various crops.

In recent years, deep learning (DL) has been regarded as a breakthrough technology in machine learning and data mining in agricultural remote sensing. Most DL algorithms including Stacked Sparse Autoencoder (SSAE), Convolutional Neural Network (CNN) and Recurrent Neural Network (RNN) have been used for yield prediction. Researchers generally believe that CNN can explore more spatial features, and LSTM, as a special RNN, has the ability to reveal phenological features, however, little attention has been paid to combining the advantages of CNN and LSTM for yield prediction.

SUMMARY OF THE INVENTION

In order to overcome the above-mentioned shortcomings of the prior art, the present invention provides a county-level crop yield estimation method based on CNN-LSTM. This method proposes a CNN-LSTM model for end-of-season crop yield prediction at the county level.

The present invention provides a county-level crop yield estimation method based on CNN-LSTM, which is characterized by comprising the following steps:

S1. Data acquisition and processing: respectively obtain the crop remote sensing data RS, environmental remote sensing data ENV, crop classification data D1, county-level soybean yield data D3, and county boundary vector data D2 of the target year of the production area to be estimated, and analyze them respectively. A plurality of described crop remote sensing data RS and a plurality of described environmental remote sensing data ENV of the same year are preprocessed to obtain the annual multi-temporal crop remote sensing data set I _RS and environmental remote sensing data set I _ENV respectively;

S2. Superposition and filtering of data: The multi-temporal crop remote sensing data set I _RS in S1 and the I _ENV of the corresponding year are superimposed and synthesized to obtain the annual data set I _RS & _ENV , combined with S1 In the described D1 of the corresponding year, the non-target crop pixels are filtered out to the data collection I _RS & _ENV , to obtain the annual data collection DI _RS & _ENV that only comprises the target crop pixel information;

S3. Obtain county-level feature tensor data: Use the county boundary vector data D2 in S1 to cut each of the data sets DI _RS & _ENV to obtain the annual data set for each county DI _RS & _ENV , count the distribution of the band pixel histogram of each band in the data set DI _RS & _ENV for each county each year, and calculate the band pixels corresponding to the data set DI _RS & _ENV for each county each year Convert the histogram into a deep learning tensor to obtain the annual feature tensor data of each county;

S4. Build and train a CNN-LSTM model: build a CNN-LSTM model, and train the CNN-LSTM model by using a plurality of the feature tensor data in S3 and the county-level soybean yield data D3 in the corresponding year of the corresponding county in S1 , and evaluate the model accuracy;

S5: Use the CNN-LSTM model trained in S4 to estimate the yield of the target crop.

Further, the crop remote sensing data RS includes t time-phase data, each time-phase data has m bands, and the environmental remote sensing data ENV includes u time-phase data, each time-phase data has n bands.

Further, the S1 also includes:

S11, carry out cloud removal processing to each described crop remote sensing data RS, classify the crop remote sensing data RS through the cloud removal processing by year, and superimpose and synthesize multiple described crop remote sensing data RS of the same year processing, the annual multi-temporal crop remote sensing dataset _IRS can be obtained;

The multiple pieces of the environmental remote sensing data ENV in the same year are superimposed and synthesized to obtain the annual environmental remote sensing data set I _ENV .

Further, the environmental remote sensing data ENV includes surface temperature data and weather parameter data, and obtaining the environmental remote sensing data set I _ENV also includes the following steps:

S11a, aligning the surface temperature data and the weather parameter data with the crop remote sensing data RS, respectively, to obtain the aligned surface temperature data and the weather parameter data;

S11b, in units of years, respectively process the aligned surface temperature data and the weather parameter data in S12 to obtain an annual multi-temporal surface temperature data set and an annual multi-temporal weather parameter data set;

S11c: Process the annual multi-temporal surface temperature data set in S13 and the multi-temporal weather parameter data set of the corresponding year to obtain the environmental remote sensing data set I _ENV of the corresponding year.

Further, the data set I _RS & _ENV in S2 has t time phases, and each time phase data has m+n frequency bands.

Further, the S3 also includes:

S31, confirm the true distribution limit of the target crop pixel data in each band: the target crop pixel data of each band in each of the crop remote sensing data RS and the environmental remote sensing data ENV in the production area to be estimated respectively form a corresponding band pixel histogram Figure, obtain the maximum value and the minimum value of the target crop pixel data of each waveband in each described crop remote sensing data RS and described environmental remote sensing data ENV, take the maximum value and minimum value of each waveband target crop pixel data as the limit, the The target crop pixel data range of each band in each of the crop remote sensing data RS and the environmental remote sensing data ENV is divided into w intervals;

S32, using the county boundary vector data D2 in S1 to cut the data collection DI _RS & _ENV in S2 respectively, to obtain the annual data collection DI _RS & _ENV of each county, count each The band pixel histogram of each band in the annual data set DI _RS & _ENV of the county area, the interval of each band is obtained in S31, and the bands of each band in the annual data set DI _RS & _ENV of each county area are obtained. The pixels are converted into the interval of the corresponding band, and the feature tensor data corresponding to each county can be obtained each year, and the shape of the feature tensor data is t*(m+n)*1*w.

Further, the S4 also includes the following steps:

S41, set training parameters: set the verification data split ratio of the CNN-LSTM model to 0.2, the training period to 100 times, and the data block size to 16;

S42. Input multiple feature tensor data into the set CNN-LSTM model, and use the fit function in keras to train the model, using the county-level soybean yield data D3 in the corresponding year of the county in S1 as training data to train the CNN-LSTM model.

The beneficial effects brought by the technical solution provided by the present invention are as follows: the method of the present invention extracts more features of county crops by means of histogram statistics and converts them into tensors based on the remote sensing data reflecting the growth state of crops and the environmental data affecting the growth of crops , to train the extracted CNN-LSTM model, which effectively improves the small-scale crop yield estimation accuracy, which is of great significance to crop yield estimation, crop market planning, crop insurance and harvest management.

Description of drawings

1 is a flowchart of a county-level crop yield estimation method based on CNN-LSTM of the present invention;

Fig. 2 is the framework diagram of CNN-LSTM proposed by the present invention;

Figure 3 is a schematic diagram of a tensor of a temporal data feature in a county;

Fig. 4 is the result graph that adopts the method of the present invention to measure soybean yield and the comparison graph of the official announcement yield graph;

Fig. 5 is a graph of the precision result of using the method of the present invention to predict soybean yield in 2011-2015;

Figure 6 shows the relationship between the five-year predicted value and the observed value.

Detailed ways

In order to make the objectives, technical solutions and advantages of the present invention clearer, the embodiments of the present invention will be further described below with reference to the accompanying drawings.

In the invention, taking 15 states in the central United States as the case area, based on the remote sensing data and historical yield data from 2003 to 2015 as the data source, the soybean yield in the region from 2011 to 2015 was estimated, and the estimated results were compared with the real ones. Compare with official data.

Selection of production areas to be estimated: According to the soybean planting distribution published by the United States Department of Agriculture (USDA), soybeans are grown in 31 states, in this case, 15 states were selected as examples, including North Dakota, South Dakota, Nebraska California, Minnesota, Iowa, Kansas, Missouri, Arkansas, Mississippi, Tennessee, Illinois, Indiana, Ohio, Michigan and Wisconsin. The soybean planting area in the above 15 states accounts for 88.75% of the national soybean planting area. According to the soybean growing period, the time period of remote sensing data is selected from April to December.

Here, the target crops appearing below are all soybeans.

Get data on the region to be estimated, including:

1. County-level soybean yield data: County-level soybean yield data from 2003 to 2015 were collected from the USDA, and the yield unit is kg/ha. Source: https://www.nass.usda.gov/Quick_Stats/ Lite/index.php;

2. Crop classification remote sensing data: Collect the crop classification remote sensing data from 2003 to 2015 from the crop data layer (CDL) issued by the USDA, one per year, which is a crop-specific land cover data layer with a high resolution 30m, obtained from: https://nassgeodata.gmu.edu/CropScape/;

3. The county boundary vector remote sensing data of the production area to be estimated: its format is shp, obtained from: https://catalog.data.gov/dataset/tiger-line-shapefile-2016-nation-u-s-current-county-and- equivalent-national-shapefile;

4. Crop remote sensing data: Download the crop remote sensing data (RS) in soybean growing period from 2003 to 2015 from the designated website (access address: https://doi.org/10.5067/MODIS/MOD09A1.006), among which , the crop remote sensing data is MOD09A1 V6 product, and the estimated value of the surface spectral reflectance of Terra MODIS bands 1-7 at 500m resolution is provided in the RS each year, and the gas, aerosol and Rayleigh scattering etc. Corrected for atmospheric conditions, this guarantees the absence of clouds or cloud shadows and aerosol loadings in each of the RSs, for each pixel, based on high observational coverage, low viewing angles.

5. Environmental remote sensing data: Environmental data includes surface temperature data and weather parameter data. Download the surface temperature data and weather parameter data during the soybean growing period from 2003 to 2015 from the designated website.

Among them, the surface temperature data adopts MOD11A2 V6 product, and the address is https://doi.org/10.5067/MODIS/MOD11A2.006. The product provides average 8-day surface temperature (LST) in a 1km*1km grid. Each pixel value in MOD11A2 is a simple average of all corresponding MOD11A1 LST pixels collected over this 8-day period. The daytime and nighttime surface temperature bands of the two data bands were used as soil long-term factors.

The weather parameter data adopts the Daymet product, available at https://daymet.ornl.gov, which is a collection of grid estimates of daily weather parameters generated by interpolation and extrapolation of daily meteorological observations. . In the present invention, two important weather parameters in Daymet are selected: two characteristic bands of precipitation and air pressure, and the resolution is 1km.

Please refer to Figure 1, a county-level crop yield estimation method based on CNN-LSTM, including the following steps:

S1. Data acquisition and processing: According to the above method, multiple pieces of crop remote sensing data RS, multiple pieces of environmental remote sensing data ENV, multiple pieces of crop classification data and multiple pieces of county-level soybean yield data D3 were obtained from 2003 to 2015, respectively. And one piece of county boundary vector data D2, and preprocesses multiple pieces of the RS and multiple pieces of the ENV in the same year, respectively, to obtain the annual multi-temporal crop remote sensing data set I _RS and environmental remote sensing data set I _ENV .

The S1 further includes: S11, in the google earth engine, perform cloud removal processing on each piece of the crop remote sensing data RS, classify the crop remote sensing data RS after the cloud removal processing by year, and classify the crop remote sensing data RS by year, and the same year A plurality of described crop remote sensing data _RSs are superimposed and synthesized, and the annual described multi-temporal crop remote sensing data sets IRS can be obtained;

A plurality of described environmental remote sensing data ENVs in the same year are carried out superimposed synthesis processing, so that the annual described I _ENVs can be obtained;

Wherein, during the growth period of soybean, the multi-phase crop remote sensing data set _{IRS includes t3 time-phase data, each time-phase data has m bands, and the I ENV includes t4} time _- phase data, each time-phase data The time-phase data has n bands, t ₃ =t ₄ =34, m=7, n=4.

Here, it should be noted that each year's multi-temporal crop remote sensing dataset I _RS and environmental remote sensing dataset I _ENV should cover all crop remote sensing data RS, surface temperature data, and weather parameter data throughout the growing season of the target crop. Since the periods of different remote sensing data are different, in the present invention, the remote sensing data with the longest return visit period shall prevail, and other remote sensing data are aligned with the remote sensing data with the longest return visit period in time by taking the mean value of multi-temporal data. The two types of data of the last crop remote sensing data RS and the environmental remote sensing data ENV are both t ₂ , and t ₂ =8.

The specific operation is: the product cycle of crop remote sensing data RS is 8 days. During the growth period, 34 crop remote sensing data need to be acquired every year, and the crop remote sensing data acquired in the same year are classified and superimposed. engine) to obtain a multi-temporal crop remote sensing dataset _IRS every year; the environmental data includes surface temperature data and weather parameter data, wherein the product of the surface temperature data provides an average 8-day surface temperature, correspondingly, in the growth During the period, 34 pieces of surface temperature data need to be obtained every year, and the surface temperature data obtained in the same year are classified and superimposed (in Google Earth Engine) to obtain a multi-temporal surface temperature data set for each year. , and the product of weather parameter data, each day corresponds to one phase data, because the product cycle of crop remote sensing data RS and surface temperature data is 8 days, and the product cycle of weather parameter data is 1 day, so it is necessary to first convert the weather parameter data The product is averaged and resampled to unify its cycle to 8 days to align the product cycle of weather parameter data with the product cycle of crop remote sensing data RS and surface temperature data. Therefore, during the growth period, the aligned weather parameters Each year also contains 34 weather parameter data. By classifying the aligned weather parameter data in the same year, the multi-temporal weather parameter data set for each year can be obtained. The annual multi-temporal surface temperature dataset and the same-year multi-temporal weather parameter dataset are classified as the same-year environmental remote sensing dataset I _ENV .

S2. Superposition and filtering of data: The multi-temporal crop remote sensing data set I _RS in S1 and the I _ENV of the corresponding year are superimposed and synthesized to obtain the annual data set I _RS & _ENV , combined with S1 In the D1 of the corresponding year in the data set I _RS & _ENV , the non-target crop pixels are filtered out, so as to obtain the data set DI _RS & _ENV which only contains the target crop pixel information every year.

Wherein, in the growth period, the data set I _RS & _ENV and the data set DI _RS & _ENV both include t ₁ time-phase data, each time-phase data has m+n bands, t ₁ =34, m+n=11.

S3, using the county boundary vector data D2 in S1 to cut each piece of the data set DI _RS & _ENV to obtain the annual data set DI _RS & _ENV of each county, and count each county every year The band pixel histogram distribution of each band in the data collection DI _RS & _ENV , and converting the band pixel histogram corresponding to the DI _RS & _ENV in each county each year into a deep learning tensor, you can obtain each Yearly feature tensor data for each county;

The S3 further includes: S31. Confirm the real distribution limits of the target crop pixel data in each band: use the ui.Chart.image.histogram function in the google earth engine to separately treat each piece of the crop remote sensing data RS and the crop in the estimated production area. The target crop pixel data statistics of each band in the environmental remote sensing data ENV form a corresponding band pixel histogram, and obtain the maximum value of the target crop pixel data of each band in each of the described crop remote sensing data RS and the environmental remote sensing data ENV and the minimum value, and then obtain the real distribution limit of the target crop pixel data of each band in each described crop remote sensing data RS and described environmental remote sensing data ENV, with each described crop remote sensing data RS and described environmental remote sensing data ENV The maximum value and the minimum value of the target crop pixel data in each band are the limits, and the range of each band target crop pixel data in each of the crop remote sensing data RS and the environmental remote sensing data ENV is divided into w intervals;

S32, using the county boundary vector data D2 in S1 to cut the data set DI _RS & _ENV in S2 respectively to obtain the data set DI _RS & _ENV for each county each year, using google earth The ui.Chart.image.histogram function in the engine counts the band pixel histogram of each band in the data set DI _RS & E _NV in each county every year, using the ee.Reducer.fixed Histogram function in the google earth engine, Using the interval of each band obtained in S31, convert the band pixels of each band in the data collection DI _RS & _ENV of each county each year into the interval of the corresponding band, and then the feature tensor corresponding to each county year can be obtained. data, the shape of the feature tensor data is t*(m+n)*1*w.

Here, it should be noted that this method assumes that when crop pixel information is used to estimate yield, it has nothing to do with the geographic location of the pixel, so the image can be reconstructed by counting the band pixel histogram of each band, that is, the histogram can be used to represent the original image. image. In the case of a given number of intervals, the larger the range of target crop pixel data in each band, the larger the obtained histogram interval, but usually the theoretical limit value of each band target crop pixel data in the obtained remote sensing data will be too large. Therefore, in order to generate regularized feature tensor data, it is necessary to find out the real distribution limit of the target crop pixel data in each band before transforming the pixel histogram of each band. Therefore, the purpose of S31 is to obtain the real distribution of the target crop pixel range in each band, and then obtain the minimum and maximum values of each band target crop pixel data in each of the crop remote sensing data RS and the environmental remote sensing data ENV, to Generate regularized tensors. In the present invention, the pixel histogram of each band is counted by adopting the method of histogram conversion, and the specific operation is as follows: to convert the maximum value of the target crop pixel data of each band in each of the crop remote sensing data RS and the environmental remote sensing data ENV The value and the minimum value are the limits, and the target crop pixel data range of each band in each of the crop remote sensing data RS and the environmental remote sensing data ENV is divided into w intervals, and then the target crop pixel value in each band is assigned to Corresponding w intervals, and count the number of target crop pixels in each interval, so each band corresponds to a 1*w feature tensor data. After converting the band pixel histogram of each band in the annual data set DI _RS & _ENV of each county into the interval of the corresponding band, the corresponding annual data DI _RS & _ENV of each county will correspond to a tensor T _RS & _ENV , the shape of the feature tensor data is t*(m+n)*1*w, in the present invention, w is 32, t=t ₃ =t ₄ =34, that is, the shape of the feature tensor data in each county per year is 34*11*1*32. Here, it should be noted that w=32 is only a specific embodiment of the present invention, and in practical applications, w may also be 64, 128, or the like.

Figure 3 is a schematic diagram of each band being converted into a 32-box histogram at a certain moment in Marion County in a Kansas state.

S4. Construct a CNN-LSTM model, and use the plurality of the feature tensor data in S3 and the county-level soybean yield data D3 in the corresponding county in S1 to train the CNN-LSTM model, and evaluate the model accuracy. When the accuracy of the model after training meets the requirements, the final CNN-LSTM model is obtained. If not, the number of training times is increased until the accuracy of the model after training meets the requirements; wherein, the CNN-LSTM model includes two parts: CNN and LSTM , CNN is used for feature extraction on feature tensor data, and LSTM is used to learn the features extracted by CNN. When training the model, the early stopping method is used to monitor the val_loss parameter. If the parameter is not optimized for 15 consecutive epochs, the training will be stopped.

Among them, CNN contains two Conv2Ds, the first Conv2D has 32 filters, the second has 64 filters, and the convolution kernel size of each filter is 1*2. Each convolutional layer is immediately followed by a pooling layer with a kernel size of 1*2 and is processed using batch normalization. Via a temporal distribution wrapper, two stacked Conv2D layers are applied to each temporal slice of the input for feature extraction. The Conv2D-processed timing tensors are then flattened and batch-normalized before being fed into the LSTM layers. There is only one LSTM layer in the LSTM part. The number of neurons of the LSTM is set to 256, followed by a dense layer with 64 neurons. Afterwards, all time outputs are flattened into a long vector, which is sent with a dropout probability of 0.5 to a Dropout layer that can randomly turn off a certain percentage of neurons during training, which can help prevent groups of neurons from becoming overwhelmed fit. Finally, a dense layer of neurons is used to output the predicted output. The activation function of the model adopts the modified linear unit ReLU. The model is built with keras. The structure of the CNN-LSTM model constructed by the present invention is shown in FIG. 2 .

Training the CNN-LSTM model includes the following steps:

S41. Set training parameters: set the validation data split ratio (validation_split) of the CNN-LSTM model to 0.2, the training period to 100 times, and the data block size to 16;

S42. Input multiple feature tensor data into the set CNN-LSTM model, and use the fit function in keras to train the model, using the county-level soybean yield data D3 in the corresponding year of the county in S1 as training data to train the CNN-LSTM model.

Here, it should be noted that the input of the CNN-LSTM model is the annual feature tensor data of each county, the label data is the historical output data of the corresponding year in the corresponding county in S1, and the output is the predicted output of the corresponding year in the corresponding county. Input data to the CNN-LSTM model, use the fit function in keras to train the model, use the historical output data of the corresponding year in the corresponding county in S1 as the training data, and set the validation data split ratio (validation_split) of the CNN-LSTM model to 0.2, The purpose is to randomly separate 20% of the training data as validation data. The training process is set for a total of 100 epoch training cycles, the patch size (data block size) is set to 16, and the mean square error crop training metric is used. To improve practical performance, an early stopping method is adopted to reduce the generalization error of deep learning systems. When the monitored "val_loss" metric stops improving after 10 consecutive epochs, the training will end and the model will be saved. The optimizer uses Adaptive Momentum (ADAM).

S5: Use the CNN-LSTM model trained in S4 to estimate the yield of the target crop, so as to obtain the target crop yield at the county level in the production area to be estimated.

In addition, in order to verify the effectiveness of this method, the method of the present invention was used to estimate the soybean yields of the above-mentioned 15 states from 2011 to 2015, respectively, and the estimated values were compared with the output values of the corresponding counties published by the USDA in the corresponding years. For comparison, RMSE and R ² are used as evaluation indicators. The accuracy evaluation standard is the determination coefficient R ² and the root mean square error RMSE. The value range of R ² is [0, 1]. The larger the R ² is, the better the model is. The better the fitting effect, the worse the model fitting effect; the smaller the RMSE, the better the model fitting effect.

Figure 4 shows the comparison results between the yield estimation results obtained by the method of the present invention and the results published by the USDA.

The RMSE is shown in Figure 5 by evaluating the production accuracy of each year.

Comparing the 5-year forecast data with the data published by USDA, we get R ² =0.78 (when R ² >0.6, it can represent the high accuracy of the forecast result), as shown in Figure 6, which proves the effectiveness of this method.

In this document, the related terms such as front, rear, upper and lower are defined by the positions of the components in the drawings and the positions between the components, which are only for the clarity and convenience of expressing the technical solution. It should be understood that the use of the locative words should not limit the scope of protection claimed in this application.

The above-described embodiments and features of the embodiments herein may be combined with each other without conflict.

The above are only preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included in the protection of the present invention. within the range.

Claims (6)

1. a county-level crop yield estimation method based on CNN-LSTM, is characterized in that: comprise the steps: S1. Data acquisition and processing: respectively obtain the crop remote sensing data RS, environmental remote sensing data ENV, crop classification data D1, county-level soybean yield data D3, and county boundary vector data D2 of the target year of the production area to be estimated, and analyze them respectively. A plurality of described crop remote sensing data RS and a plurality of described environmental remote sensing data ENV of the same year are preprocessed to obtain the annual multi-temporal crop remote sensing data set I _RS and environmental remote sensing data set I _ENV respectively; S2. Superposition and filtering of data: The multi-temporal crop remote sensing data set I _RS in S1 and the environmental remote sensing data set I _ENV of the corresponding year are superimposed and synthesized to obtain an annual data set I _RS & _ENV , in conjunction with the crop classification data D1 described in S1, filter out the non-target crop pixels on the data aggregate I _RS & _ENV of the corresponding year, to obtain an annual data aggregate DI _RS &_ENV; S3. Obtain county-level feature tensor data: Use the county boundary vector data D2 in S1 to cut each of the data sets DI _RS & _ENV to obtain the annual data set for each county DI _RS & _ENV , calculate the distribution of the band pixel histogram of each band in the data set DI _RS & _ENV for each county each year, and calculate the band pixel histogram corresponding to the data set DI _RS & _ENV for each county each year Convert to deep learning tensor to get the annual feature tensor data of each county; S4. Build and train a CNN-LSTM model: build a CNN-LSTM model, and train the CNN-LSTM model by using a plurality of the feature tensor data in S3 and the county-level soybean yield data D3 in the corresponding year of the corresponding county in S1 ; S5: Use the CNN-LSTM model trained in S4 to estimate the yield of the target crop, so as to obtain the target crop yield at the county level in the production area to be estimated. 2. a kind of county-level crop yield estimation method based on CNN-LSTM according to claim 1, is characterized in that: described S1 also comprises: S11, carry out cloud removal processing to each described crop remote sensing data RS, classify the crop remote sensing data RS through the cloud removal processing by year, and superimpose and synthesize multiple described crop remote sensing data RS of the same year processing, the annual multi-temporal crop remote sensing dataset _IRS can be obtained; The multiple pieces of the environmental remote sensing data ENV in the same year are superimposed and synthesized to obtain the annual environmental remote sensing data set I _ENV . 3. a kind of county-level crop yield estimation method based on CNN-LSTM according to claim 2, is characterized in that: described ENV comprises surface temperature data and weather parameter data, obtains described environmental remote sensing data set I _ENV also comprises The following steps: S11a, aligning the surface temperature data and the weather parameter data with the crop remote sensing data RS, respectively, to obtain the aligned surface temperature data and the weather parameter data; S11b, in units of years, respectively process the aligned surface temperature data and the weather parameter data in S11a to obtain an annual multi-temporal surface temperature data set and an annual multi-temporal weather parameter data set; S11c: Process the annual multi-temporal surface temperature data set in S13 and the multi-temporal weather parameter data set of the corresponding year to obtain the environmental remote sensing data set I _ENV of the corresponding year. 4. a kind of county-level crop yield estimation method based on CNN-LSTM according to claim 1, is characterized in that: described S3 also comprises: S31, confirm the true distribution limit of the target crop pixel data in each band: the target crop pixel data of each band in each of the crop remote sensing data RS and the environmental remote sensing data ENV in the production area to be estimated respectively form a corresponding band pixel histogram Figure, obtain the maximum value and the minimum value of the target crop pixel data of each waveband in each described crop remote sensing data RS and described environmental remote sensing data ENV, take the maximum value and minimum value of each waveband target crop pixel data as the limit, the The target crop pixel data range of each band in each of the crop remote sensing data RS and the environmental remote sensing data ENV is divided into multiple intervals; S32, using the county boundary vector data D2 in S1 to cut the data collection DI _RS & _ENV in S2 respectively, to obtain the annual data collection DI _RS & _ENV of each county, count each The band pixel histogram of each band in the annual data set DI _RS & _ENV of the county area, the interval of each band is obtained in S31, and the bands of each band in the annual data set DI _RS & _ENV of each county area are obtained. Convert the pixels to the corresponding band range, and you can get the feature tensor data corresponding to each county each year. 5. a county-level crop yield estimation method based on CNN-LSTM according to claim 1, is characterized in that: described CNN-LSTM model comprises CNN and LSTM, wherein, CNN is used for characteristic tensor data to be characterized Extraction, LSTM is used to learn the features extracted by CNN. 6. a kind of county-level crop yield estimation method based on CNN-LSTM according to claim 1, is characterized in that: described S4 also comprises the following steps: S41, set training parameters: set the verification data split ratio of the CNN-LSTM model to 0.2, the training period to 100 times, and the data block size to 16; S42. Input multiple feature tensor data into the CNN-LSTM model set in S41, and use the fit function in keras to train the model. The county-level soybean yield data D3 in the corresponding year of the corresponding county in S1 is Training data to train the CNN-LSTM model.

Images