A method for predicting rice yield by coupling a crop model with remote sensing data

Abstract

The present invention discloses a rice yield prediction method for coupling a crop model and remote sensing data, which includes fusing the time increment and the space increment based on the vegetation index spatio-temporal data fusion algorithm and the constrained least squares method to obtain the final combined increment. The final combined increment is input into the trained LAI remote sensing estimation model to output the remotely sensed LAI index. A cost function is constructed according to the LAI index simulated by the WOFOST model and the remotely sensed LAI index. The SCE-UA algorithm is used to iteratively optimize the parameters to be optimized of the WOFOST model with the goal of minimizing the cost function, so as to obtain the optimal input parameters of the WOFOST model. The optimal input parameters and other parameters to be optimized are input into the WOFOST model to output the predicted rice yield result. The present invention solves the problem that the spatial mismatch between the field scale and the remote sensing scale cannot be matched.

Description

Technical Field

[0001] The present invention relates to the technical field of rice yield prediction, and in particular to a rice yield prediction method that couples a crop model with remote sensing data. Background Art

[0002] Traditional rice yield prediction mainly relies on historical statistical data and empirical judgment. Although this method is simple and easy to implement, it is difficult to fully consider the complexity and uncertainty in the crop growth process, and the prediction accuracy is limited. With the rapid development of information technology, the emergence of crop models has opened up a new path for rice yield prediction. With the continuous development of remote sensing technology, remote sensing products with high frequency, multi-band, and multi-spatial resolution, with their characteristics of strong timeliness and wide monitoring range, and the advantages of strong mechanism and continuous time of crop models, form a good complementary relationship. Currently, almost all crop models are established at the field scale, resulting in a mismatch in space between the field scale and the remote sensing scale.

[0003] Therefore, a rice yield prediction method that couples a crop model with remote sensing data is developed to solve the above problems. Summary of the Invention

[0004] The present invention proposes a rice yield prediction method that couples a crop model with remote sensing data to solve the problem that currently almost all crop models are established at the field scale, resulting in a mismatch in space between the field scale and the remote sensing scale in the prior art.

[0005] The present invention achieves the above object through the following technical solutions:

[0006] A rice yield prediction method that couples a crop model with remote sensing data according to the present invention includes:

[0007] Obtaining information, where the information includes multi-temporal remote sensing data in the rice field, ground-measured LAI data, the longitude and latitude of the LAI measurement points, and the meteorological, soil, and crop parameters required in the WOFOST model. The multi-temporal remote sensing data includes medium- and high-resolution remote sensing data sets;

[0008] Calculating a normalized difference vegetation index data set according to the multi-temporal remote sensing data;

[0009] Based on the normalized difference vegetation index data set, obtaining a time increment through linear unmixing analysis;

[0010] Based on the normalized difference vegetation index data set, obtaining a spatial increment through thin plate spline interpolation;

[0011] Fusing the time increment and the spatial increment based on a vegetation index spatio-temporal data fusion algorithm and a constrained least squares method to obtain a final combined increment;

[0012] Calculate the vegetation index value estimated by remote sensing at the LAI measurement point based on the longitude and latitude of the LAI measurement point and the multi-temporal remote sensing data;

[0013] Divide the ground-measured LAI data and the vegetation index value estimated by remote sensing at the LAI measurement point into a training sample set and a test sample set according to a proportion. Construct an LAI remote sensing estimation model based on the random forest algorithm. Train and test the LAI remote sensing estimation model according to the training sample set and the test sample set to obtain a trained LAI remote sensing estimation model. Input the final combined increment into the trained LAI remote sensing estimation model to output the remotely sensed LAI index;

[0014] Construct a WOFOST model based on the meteorological, soil and crop parameters and the remotely sensed LAI index. Conduct a sensitivity analysis on the input parameters of the WOFOST model based on the Sobol index in the SALib module to obtain the parameters to be optimized;

[0015] Simulate the LAI index based on the WOFOST model according to the parameters to be optimized;

[0016] Construct a cost function based on the LAI index simulated by the WOFOST model and the remotely sensed LAI index. Use the SCE-UA algorithm to iteratively optimize the parameters to be optimized with the goal of minimizing the cost function to obtain the optimal input parameters of the WOFOST model. Input the optimal input parameters and other parameters to be optimized into the WOFOST model to output the predicted rice yield result.

[0017] Further, the calculation formulas for the normalized vegetation index with high temporal resolution and low spatial resolution and high spatial resolution and low temporal resolution based on the multi-temporal remote sensing data are as follows:

[0018] ,

[0019] In the formula, 、 are the reflectances of the red band and the near-infrared band of the HLS data with a resolution of 30 meters and the MODIS data with a resolution of 1000 meters, respectively.

[0020] Further, based on the normalized vegetation index dataset, a time increment is obtained through linear unmixing analysis, including:

[0021] ,

[0022] ,

[0023] Time increment , including:

[0024]

[0025] where n is the number of low - resolution pixels, and l is the number of land - cover types within the moving window. is the low - resolution pixel that can be directly obtained from the low - resolution NDVI time - series image. of the NDVI increment. is the high - resolution NDVI increment corresponding to the c land type within the moving window. is the fraction of the land type within the low - resolution pixel inside. is the set of all low - resolution NDVI increments within the moving window. , and respectively represent the minimum value, maximum value, and standard deviation.

[0026] Furthermore, based on the normalized vegetation index dataset, a spatial increment is obtained based on the thin - plate spline interpolation method, including:

[0027] ,

[0028] where represents the spatial increment, and are the interpolation values of the thin - plate spline interpolation method at and times respectively, represents the j - th high - resolution pixel in the low - spatial resolution ( ).

[0029] Furthermore, the time increment and the spatial increment are fused based on the vegetation - index spatio - temporal data fusion algorithm and the constrained least - squares method to obtain the final combined increment, including:

[0030] ,

[0031] ,

[0032] ,

[0033] where is the NDVI value of the high - resolution pixel at time is the NDVI value of the high - resolution pixel predicted by the algorithm at time, is the high - resolution pixel The corresponding combined increment is the residual within the low-resolution pixel ; and are the weight coefficients of the spatial increment and the temporal increment, and are the spatial increment and the temporal increment respectively, and m is the number of high-resolution pixels within the low-resolution pixel.

[0034] Furthermore, calculating the remotely sensed estimated vegetation index value at the LAI measured point according to the longitude and latitude of the LAI measured point and the multi-temporal remote sensing data includes:

[0035] Using ARCGIS software to match the position of the measured point in the remote sensing index image according to the longitude and latitude information of the LAI measured point;

[0036] Constructing a 3×3 window with the pixel where the measured point is located as the central pixel;

[0037] Calculating the average value of the 9 pixels within the 3×3 window as the remotely sensed estimated vegetation index value at the measured point.

[0038] Furthermore, the test sample set is used to test the accuracy of the random forest model, and the formulas for the test indexes bias (Bias) and root mean square error (RMSE) are as follows:

[0039] ,

[0040] ,

[0041] In the formula, is the measured value of LAI, is the remotely sensed estimated LAI, n1 is the number of samples, and i' is the i'th sample.

[0042] Furthermore, the input parameters include: the effective accumulated temperature from emergence to flowering, the effective accumulated temperature from flowering to maturity, the initial total dry weight, the leaf area index at emergence, the specific leaf area at the crop development stage, the leaf senescence index, the initial light energy utilization efficiency of a single leaf, the maximum photosynthesis rate of the leaf, the maximum synthesis rate reduction factor (TMPF) of leaf CO 2 2, the conversion rate of leaf growth assimilates, the conversion rate of storage organ growth assimilates, the conversion rate of root growth assimilates, and the conversion rate of stem growth assimilates.

[0043] Furthermore, the cost function is:

[0044] ,

[0045] In the formula, is the cost function; x1 is the parameter to be optimized, is the expected value of the parameter to be optimized, B is the error covariance matrix of x1, i is the number of remote sensing observations, and n2 is the total number of remote sensing observations. is the vegetation index value of the i-th remote sensing inversion. is the error covariance matrix of is the LAI of the i-th remote sensing observation simulated by the WOFOST model, and the superscript T represents the transpose of the matrix.

[0046] The beneficial effects of the present invention are as follows:

[0047] A rice yield prediction method coupling a crop model and remote sensing data proposed by the present invention solves the problem that almost all current crop models are established at the field scale, resulting in the inability to match the field scale and the remote sensing scale spatially. Through the coupling of the crop model and remote sensing data, the prediction of rice yield at the field scale is realized. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 is a 3*3 window centered on the pixel where the LAI measurement site is located in the embodiment of the present application;

[0049] Figure 2 is a technical roadmap for generating combined increments based on the CLS algorithm in the embodiment of the present application;

[0050] Figure 3 is a flowchart for sensitivity analysis of WOFOST model parameters in the embodiment of the present application;

[0051] Figure 4 is a flowchart for predicting rice yield by assimilating the WOFOST model and remote sensing data in the embodiment of the present application. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0052] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. The components of the embodiments of the present invention described and illustrated herein can be arranged and designed in various different configurations.

[0053] Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed invention, but merely represents selected embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention.

[0054] It should be noted that similar reference numerals and letters refer to similar items in the following figures. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0055] The following will, with reference to the accompanying drawings, elaborate on the specific embodiments of the present invention in detail.

[0056] As Figure 4 shown, a rice yield prediction method for coupling a crop model with remote sensing data includes:

[0057] Obtaining information, where the information includes multi-temporal remote sensing data in the rice field, ground-measured LAI data, the longitude and latitude of the LAI measurement points, and meteorological, soil, and crop parameters required in the WOFOST model. The multi-temporal remote sensing data includes medium- and high-resolution remote sensing data sets;

[0058] Among them, the required meteorological parameters include radiation value, air temperature, water vapor pressure, wind speed, and precipitation. The required soil parameters include saturated water content, soil water content, field capacity, and wilting coefficient. The crop parameters include leaf area index, biomass, accumulated temperature, effective temperature, and light energy utilization rate of a single leaf;

[0059] Calculating a normalized difference vegetation index data set based on the multi-temporal remote sensing data;

[0060] Based on the normalized difference vegetation index data set, obtaining a time increment through linear unmixing analysis;

[0061] Based on the normalized difference vegetation index data set, obtaining a spatial increment through thin plate spline interpolation;

[0062] Fusing the time increment and the spatial increment based on a vegetation index spatio-temporal data fusion algorithm and a constrained least squares method to obtain a final combined increment;

[0063] Calculating the remotely sensed vegetation index value at the LAI measurement points based on the longitude and latitude of the LAI measurement points and the multi-temporal remote sensing data;

[0064] Dividing the ground-measured LAI data and the remotely sensed vegetation index value at the LAI measurement points into a training sample set and a test sample set according to a proportion. Based on the random forest algorithm, constructing an LAI remote sensing estimation model, training and testing the LAI remote sensing estimation model according to the training sample set and the test sample set to obtain a trained LAI remote sensing estimation model, and inputting the final combined increment into the trained LAI remote sensing estimation model to output a remotely sensed retrieved LAI index;

[0065] Construct a WOFOST model based on the meteorological, soil and crop parameters and the remotely sensed LAI index retrieved. Conduct a sensitivity analysis on the input parameters of the WOFOST model based on the Sobol index in the SALib module to obtain the parameters to be optimized.

[0066] Simulate the LAI index based on the WOFOST model according to the parameters to be optimized.

[0067] Construct a cost function based on the LAI index simulated by the WOFOST model and the remotely sensed LAI index retrieved. Use the SCE-UA algorithm to iteratively optimize the parameters to be optimized with the goal of minimizing the cost function, obtain the optimal input parameters of the WOFOST model, input the optimal input parameters and other parameters to be optimized into the WOFOST model, and output the predicted rice yield results.

[0068] A method for predicting rice yield by coupling a crop model and remote sensing data specifically includes:

[0069] Step 1. Preprocessing of remote sensing data

[0070] The preprocessing of remote sensing data mainly includes radiometric calibration, atmospheric correction, geometric correction, and image mosaicking and cropping.

[0071] (I) Radiometric calibration

[0072] Eliminate the systematic errors in the remote sensing image data caused by factors such as the sensor itself through radiometric calibration. The calibrated radiance or reflectance can be calculated for DN using the radiometric calibration formula.

[0073] The formula for converting the DN value of the remote sensing image to radiance using the calibration coefficient:

[0074] ,

[0075] L is the converted equivalent radiance, with the unit of ; DN is the digital quantization output value of the sensor; gain and offset are the gain and offset respectively, which can be obtained by querying the calibration information of the original remote sensing image.

[0076] (II) Atmospheric correction

[0077] Eliminate the influence of the atmosphere and light on the emission of ground objects through atmospheric correction, obtain the true physical model parameters such as the reflectance, radiance, and surface temperature of ground objects, and use them to eliminate the influence of water vapor, oxygen, etc. in the atmosphere on the reflection of ground objects, and eliminate the influence of atmospheric molecular and aerosol scattering.

[0078] The present invention uses the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) module in the ENVI remote sensing processing software for atmospheric correction. The FLAASH module has the advantages of high spectral reduction accuracy and good fidelity of the spectral information of ground objects after atmospheric correction.

[0079] (III) Geometric correction

[0080] Geometric correction is used to solve the image deformation caused by systematic and non-systematic factors, so as to achieve geometric registration with the standard image or map. The present invention uses a high-resolution image with accurate geometric position information processed in the same area as the reference image to correct the to-be-calibrated high-resolution image. If high requirements are placed on the geometric accuracy of the image, ground control points can be used, and a geometric correction model can be adopted to construct the geometric relationship between the image and the ground coordinates / between images to complete the geometric correction.

[0081] (IV) Mosaic and cropping

[0082] Remote sensing image mosaic means that when the study area exceeds the range covered by a single remote sensing image, usually two images need to be stitched together to form one or a series of larger images covering the whole area; the purpose of cropping is to remove the area outside the project, and usually the image is cropped by administrative division boundaries or natural division boundaries. The present invention uses remote sensing processing software such as ARCGIS and ENVI to realize the mosaic and cropping of remote sensing images.

[0083] Step 2: Generate high spatio-temporal resolution NDVI based on the vegetation index fusion algorithm

[0084] (I) Calculate medium- and high-resolution NDVI datasets

[0085] The medium-resolution NDVI data used in the present invention is calculated from MODIS MOD09GA using the normalized difference vegetation index NDVI formula (1), with a temporal resolution of 1 day and a spatial resolution of 1000 meters; the high-resolution NDVI data is calculated from the HLS surface reflectance product according to formula (1), with a temporal resolution of about 16 days and a spatial resolution of 30 meters.

[0086] (1)

[0087] In the formula, 、 are the reflectances of the red band and the near-infrared band of the 30-meter resolution HLS data and the 1000-meter resolution MODIS data, respectively.

[0088] (II) Obtain the time increment based on linear unmixing

[0089] According to the linear spectral mixture theory, the change of low-resolution pixel NDVI over time can be considered as a linear combination of the NDVI increments of all high-resolution pixels within this pixel in a short period. Therefore, assuming that the high-resolution pixels of the same land cover type have similar increments in the local area, the low-resolution pixel is decomposed from the reference date to the prediction date using the linear mixture model, and its calculation formula is shown in (2) and (3):

[0090] (2)

[0091] (3)

[0092] Time increment , including:

[0093]

[0094] In the formula, n is the number of low-resolution pixels, is the number of land cover types within the moving window, is the NDVI increment of the low-resolution pixel ( ) that can be directly obtained from the low-resolution NDVI time series image, is the high-resolution NDVI increment corresponding to land type c within the moving window, is the fraction of land type in the low-resolution pixel (x, y), is the set of all low-resolution NDVI increments within the moving window; , and respectively represent the minimum value, maximum value, and standard deviation of

[0095] (III) Obtaining the spatial increment based on thin plate spline interpolation

[0096] Using the thin plate spline interpolation method (TPS), the low-resolution NDVI at and moments is interpolated to high resolution. Then, the difference between the interpolation results of and is used to obtain the spatial increment. Since this increment only uses the spatial dependence between low-resolution pixels, it is called the spatial increment , as shown in formula (4).

[0097] (4)

[0098] In the formula, and are the interpolated values of TPS at and moments respectively, represents the j-th high-resolution pixel in the low spatial resolution (x, y).

[0099] (4) Predicting high spatio-temporal resolution NDVI dataset by combining time increment and spatial increment

[0100] As Figure 2 shown, a reasonable combination of the time increment ( ) and the spatial increment ( ) can improve the performance and robustness of the fusion algorithm. and are the spatial increment and time increment after downscaling of the time increment ( ) and the spatial increment ( ) respectively. is the increment of the low spatial resolution pixel, i.e., the low-resolution increment. The simplest and most effective combination method of the two increments is to add them with reasonable weight coefficients based on the least squares method. To further improve the accuracy of the combined increment, the residuals need to be distributed to each high-resolution pixel within the low-resolution pixel (x, y) . As shown in formula (5):

[0101] (5)

[0102] In the formula, is the NDVI value of the high-resolution pixel ( ) at moment, is the NDVI value of the high-resolution pixel predicted by the algorithm at moment; is the combined increment corresponding to the high-resolution pixel , calculated by formula (6); is the residual within the low-resolution pixel , calculated by formula (7), and the final combined increment

[0103] (6)

[0104] (7)

[0105] In formula (6), and are the weight coefficients of the spatial increment and the time increment, and are the spatial increment and the time increment respectively. In formula (7), is the low spatial resolution pixel from to the true increment of NDVI. According to formula (5), it is possible to predict the high-resolution NDVI at time based on the high-resolution NDVI at time ( ).

[0106] Step 3: Calculate the remotely sensed vegetation index value estimated by remote sensing at the LAI measurement points according to the longitude and latitude of the LAI measurement points and the multi-temporal remote sensing data. Train and test the LAI remote sensing estimation model according to the training sample set and the test sample set, specifically including:

[0107] Use ARCGIS software to match the position of the point in the remote sensing index image according to the longitude and latitude information of the LAI measurement points; then, construct a 3*3 window with the pixel where the measurement point is located as the central pixel ( Figure 1 shown, "O" is the pixel where the measurement point is located, and the gray shadow is the 3*3 window); finally, calculate the average value of the 9 pixels within the 3*3 window as the vegetation index value corresponding to the measurement point. Then, divide the measured LAI value and the remotely sensed vegetation index value into a training sample set and a test sample set according to a ratio of 7:3; the training sample set is mainly used for constructing a random forest model, and the test sample set is used to test the accuracy of the random forest model, thereby obtaining a trained LAI remote sensing estimation model.

[0108] The formulas for the bias (Bias) and root mean square error (RMSE) test indicators are shown in (8)-(9):

[0109] (8)

[0110] (9)

[0111] In the formula, is the measured LAI value, is the LAI estimated by remote sensing, n1 is the number of samples, and i' is the i'th sample.

[0112] Input the final combined increment into the trained LAI remote sensing estimation model, output the remotely sensed retrieved LAI index, and construct a WOFOST model according to the meteorological, soil and crop parameters and the remotely sensed retrieved LAI index.

[0113] Step 4: Determination of parameters to be optimized

[0114] As Figure 3As shown, the WOFOST model has many input parameters. However, in the actual simulation process, we found that by changing the values of different input parameters, there are certain differences in the degree of their influence on the output results. The change of some input parameters has a greater impact on the output results of WOFOST, but there are also some input parameters whose change has a smaller impact on the model output results. Therefore, on the basis of ensuring the biological significance and simulation effect of the input parameters, appropriately reducing the input parameters is of great significance for improving the running speed of the model. In the present invention, 26 input parameters related to rice growth are selected from the numerous input parameters of the WOFOST model. The model is driven to run using the Python language, and the Sobol index in the SALib module is used to conduct a global sensitivity analysis of the input parameters of the WOFOST model to obtain the first-order sensitivity, second-order sensitivity, and total-order sensitivity. The parameters to be optimized are obtained based on the first-order sensitivity, second-order sensitivity, and total-order sensitivity. The parameters to be optimized and their main value ranges are shown in Table 1:

[0115] Table 1

[0116] Parameter Value range Parameter Value range Parameter Value range TSUM1 [1300,1500] EFFTB2 [0,1] TMPF2 [0,1] TSUM2 [500,800] EFFTB3 [0,1] TMPF3 [0,1] TDWI [100,300] EFFTB4 [0,1] TMPF4 [0,1] LAIEM [0.0007,0.30] EFFTB5 [0,1] TMPF5 [0,1] SLATB1 [0.001,0.004] AMAX1 [25,50] CVL [0.5,1] SLATB2 [0.001,0.004] AMAX2 [30,60] CVO [0.5,1] SLATB3 [0.001,0.004] AMAX3 [40,80] CVR [0.5,1] SPAN [20,30] AMAX4 [30,60] CVS [0.5,1] EFFTB1 [0,1] TMPF1 [0,1]

[0117] The characterization content of each input parameter in Table 1 is as follows: TSUM1 and TSUM2 are the effective accumulated temperatures from emergence to flowering and from flowering to maturity in the WOFOST model, respectively; TDWI is the initial total dry weight; LAIEM is the leaf area index at emergence; SLATB1, SLATB2, and SLATB3 are the specific leaf areas at crop development stage (DVS) = 0.0, DVS = 0.5, and DVS = 2.0; SPAN is the leaf senescence index; EFFTB is the initial light energy utilization efficiency of a single leaf. The EFFTB values at 0°C, 10°C, 20°C, 30°C, and 40°C are EFFTB1, EFFTB2, EFFTB3, EFFTB4, and EFFTB5, respectively; AMAX is the maximum rate of leaf photosynthesis. The AMAX values at DVS = 0.0, DVS = 1.0, DVS = 1.3, and DVS = 2.0 are named AMAX1, AMAX2, AMAX3, and AMAX4, respectively; TMPF is denoted as TMPF1, TMPF2, TMPF3, TMPF4, and TMPF5 at 0°C, 10°C, 20°C, 30°C, and 40°C, respectively; CVL is the conversion rate of leaf growth assimilates; CVO is the conversion rate of assimilates for storage organ growth; CVR is the conversion rate of assimilates for root growth; CVS is the conversion rate of assimilates for stem growth.

[0118] Among them, the specific steps for conducting a sensitivity analysis of the input parameters are as follows:

[0119] In this invention, the Sobol indices in the python SALib module are used to conduct global sensitivity analysis on the important input parameters of the WOFOST model. The Sobol method is a global sensitivity analysis method based on variance decomposition, which is used to quantify the contributions of input variables and their interactions to the uncertainty of model outputs. This method is applicable to complex non-linear and high-dimensional systems.

[0120] The core idea of the Sobol method is to decompose the total variance of the model output into the contributions of each input variable and the contributions of the interactions between variables. Its mathematical concept is that for the input variable X = ( , , …, ), the variance of the model output Y = f(X) can be decomposed as:

[0121] (10)

[0122] In the formula, is the independent contribution of variable ; is the interaction contribution of variables and ; is the high-order interaction contribution of all variables.

[0123] Among them, the first-order sensitivity index represents the proportion of the contribution of the input variable to the output result. The larger the value, the greater the impact of the input variable on the output result.

[0124] (11)

[0125] The second-order interaction sensitivity index represents the proportion of the contribution of the interaction of variables and to the output result.

[0126] (12)

[0127] The total sensitivity index represents the total contribution of the main effect of the input variable and the interaction effects of other variables. The higher the value, the more important the input variable is to the output result.

[0128] (13)

[0129] The implementation process can be summarized as follows: based on the determined model and input variables, Monte Carlo sampling or Latin hypercube sampling is used to generate input variable samples, the WOFOST model is called by Python to run, the yield results are calculated, and then the total variance of the model output is calculated based on the Sobol index. Finally, the sensitivity index of each input variable is calculated by mathematical methods.

[0130] Based on the results of the global sensitivity analysis of the Sobol input parameters, the input parameters with greater sensitivity are adjusted. For the remaining parameters with smaller sensitivity and less impact on the output results, and the parameters without sensitivity analysis, they are obtained through field measurements, referring to the research results of others or the default values of the WOFOST model, so as to obtain the parameters to be optimized, making the WOFOST model more accurately simulate the growth and development process of rice and yield prediction.

[0131] Step 5: Yield prediction at the rice field scale by coupling the crop model and remote sensing data

[0132] The SCE-UA algorithm is a machine learning algorithm developed on the basis of the controlled random search algorithm and the genetic algorithm. This algorithm incorporates the complex shape segmentation and hybrid ideas, that is, it inherits the characteristics of global search and also has the characteristics of complex evolution, with the advantages of flexibility, extensiveness, and high accuracy for nonlinear optimization problems. Its cost function formula is shown in (14):

[0133] (14)

[0134] In the formula, is the cost function; x1 is the parameter to be optimized, is the expected value of the parameter to be optimized, B is the error covariance matrix of x1, i is the number of remote sensing observations, n2 is the total number of remote sensing observations, is the vegetation index value of the i-th remote sensing inversion, is the error covariance matrix of, is the LAI of the i-th remote sensing observation simulated by the WOFOST model, and the superscript T represents the transpose of the matrix.

[0135] In the present invention, LAI is selected as the optimization comparison object, and the cost function (formula 14) is constructed by using the LAI estimated by remote sensing and the LAI data simulated by the WOFOST crop model. The SCE-UA algorithm is used to iteratively optimize the input parameters of the WOFOST model, finally minimizing the cost function, and then obtaining the optimal input parameters of the model. Then, the optimal input parameters, together with other parameters to be optimized, are input into the WOFOST model to simulate the yield results after data assimilation. Other input parameters refer to other parameters to be optimized except for the optimal input parameters.

[0136] The WOFOST model was run to obtain the optimal fitting values of LAI in the study area and at the measured points, as well as LAI and TWSO. The accuracy and applicability of the rice yield prediction method based on the assimilation of crop model and remote sensing data proposed by the present invention were evaluated by comparing the R2 and RMSE between the LAI output by WOFOST and the LAI inverted from remote sensing data and the LAI measured in the field, and the R2 and RMSE values between the TWSO output by WOFOST and the yield data measured in the field.

[0137] The beneficial effects of the present invention compared with the prior art are as follows:

[0138] 1. The method of the present invention collects basic data including high-resolution remote sensing data, meteorological data, soil data, crop data, yield and LAI data measured in the field, and completes the preprocessing of high-resolution remote sensing data.

[0139] 2. Give full play to the advantages of medium- and high-resolution remote sensing data in terms of time or space scale. Based on the unmixing analysis, the time increment of medium- and high-resolution NDVI data is obtained, the thin plate spline interpolation method is used to obtain the spatial increment, the two increments are combined by the vegetation index spatio-temporal data fusion algorithm, and the final combined increment is obtained through the constrained least squares method (CLS) to generate a high spatio-temporal resolution (30 meters per day) NDVI product dataset, which well solves the problem of spatial matching between the field scale and the remote sensing scale.

[0140] 3. Use the python language to drive the model to run, and use the Sobol index in the SALib module to perform sensitivity analysis on the input parameters of the WOFOST model to obtain the first-order sensitivity, second-order sensitivity and total-order sensitivity, so as to determine the factors that need to be adjusted and considered first when "localizing" the WOFOST model parameters.

[0141] 4. Taking the LAI index as the combination point, the SCE-UA machine learning algorithm is used to construct a rice surface scale yield prediction model that couples the WOFOST model and remote sensing LAI data assimilation; and the applicability of the algorithm proposed by the present invention is verified by comparing with the rice yield results simulated by a single WOFOST model and the measured yield data.

[0142] 5. Coupling the crop model and remote sensing data can provide the possibility for realizing the quantitative estimation of crop growth and yield at the field scale.

[0143] The above are only the preferred embodiments of the present invention. It should be pointed out that for those of ordinary skill in the art, without departing from the technical principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. A rice yield prediction method by coupling crop models and remote sensing data, characterized in that: include: Acquiring information, the information including multi-temporal remote sensing data of rice fields, ground-measured LAI data, latitude and longitude of LAI measurement points, and meteorological, soil and crop parameters required by the WOFOST model, the multi-temporal remote sensing data including medium- and high-resolution remote sensing data sets; Calculating a normalized vegetation index data set based on the multi-temporal remote sensing data; According to the normalized vegetation index data set, a time increment is obtained based on a linear unmixing analysis; According to the normalized vegetation index data set, a spatial increment is obtained based on a thin plate spline interpolation method; Based on the vegetation index spatiotemporal data fusion algorithm and the constrained least squares method, the time increment and the space increment are fused to obtain a final combined increment; Calculate the remote sensing estimated vegetation index value at the LAI measured point according to the latitude and longitude of the LAI measured point and the multi-temporal remote sensing data; Divide the LAI data measured on the ground and the vegetation index value estimated by remote sensing at the LAI measuring point into a training sample set and a test sample set in proportion, build an LAI remote sensing estimation model based on a random forest algorithm, train and test the LAI remote sensing estimation model according to the training sample set and the test sample set to obtain a trained LAI remote sensing estimation model, input the final combined increment into the trained LAI remote sensing estimation model, and output the remote sensing inversion LAI index; A WOFOST model is constructed according to the meteorological, soil and crop parameters and the remote sensing inversion LAI index, and a sensitivity analysis is performed on the input parameters of the WOFOST model based on the Sobol index in the SALib module to obtain parameters to be optimized; Simulating the LAI index based on the WOFOST model according to the parameters to be optimized; A cost function is constructed according to the LAI index simulated by the WOFOST model and the LAI index inverted by remote sensing, and the parameters to be optimized are iteratively optimized using the SCE-UA algorithm with the goal of minimizing the cost function to obtain the optimal input parameters of the WOFOST model. The optimal input parameters and other parameters to be optimized are input into the WOFOST model to output the predicted rice yield result.

2. The rice yield prediction method of coupling crop model and remote sensing data according to claim 1, characterized in that: The calculation formula for calculating the normalized vegetation index based on the multi-temporal remote sensing data is as follows: ， In the formula, , These are the reflectances of the red band and near-infrared band of the HLS data with a resolution of 30 meters and the MODIS data with a resolution of 1000 meters, respectively.

3. The rice yield prediction method of coupling crop model and remote sensing data according to claim 1, characterized in that: According to the normalized vegetation index data set, a time increment is obtained based on a linear unmixing analysis, including: ， ， Time increment ,include: ， Where n is the number of low-resolution pixels, is the number of land cover types in the moving window, It is a low-resolution pixel that can be directly obtained from the low-resolution NDVI time series image ( ), is the high-resolution NDVI increment corresponding to land type c within the moving window, for Land class in low resolution pixels The score within is the set of all low-resolution NDVI increments within the moving window, , and Respectively represent The minimum, maximum and standard deviation of Indicates low spatial resolution The j-th high-resolution pixel in .

4. The rice yield prediction method of coupling crop model and remote sensing data according to claim 3 is characterized in that: According to the normalized vegetation index dataset, a spatial increment is obtained based on a thin plate spline interpolation method, including: ， In the formula, represents the space increment, and They are respectively the thin plate spline interpolation method and and The interpolated value at that moment.

5. The rice yield prediction method of coupling crop model and remote sensing data according to claim 4, characterized in that: The time increment and the space increment are fused based on the vegetation index spatiotemporal data fusion algorithm and the constrained least squares method to obtain the final combined increment, including: ， ， ， In the formula, for High resolution pixels at all times The NDVI value, Predicted by the algorithm The NDVI value of the high-resolution pixel at the moment, High resolution pixels The corresponding combined increment is the final combined increment, Low resolution pixels The residual inside and is the weight coefficient of space increment and time increment, and are the spatial increment and temporal increment respectively, and m is the number of high-resolution pixels in the low-resolution pixel.

6. The rice yield prediction method of coupling crop model and remote sensing data according to claim 1, characterized in that: Calculating the remote sensing estimated vegetation index value at the LAI measured point according to the latitude and longitude of the LAI measured point and the multi-temporal remote sensing data, including: Using ARCGIS software to match the position of the measured point in the remote sensing index image according to the latitude and longitude information of the LAI measured point; Construct a 3*3 window with the pixel where the measured point is located as the center pixel; The average value of the 9 pixels in the 3*3 window is calculated as the vegetation index value estimated by remote sensing at the measured point.

7. The rice yield prediction method of coupling crop model and remote sensing data according to claim 1, characterized in that: The test sample set is used to test the accuracy of the random forest model. The formulas for testing the index bias Bias and root mean square error RMSE are as follows: ， ， In the formula, is the measured value of LAI, is the LAI estimated by remote sensing, n1 is the number of samples, and i' is the i'th sample.

8. The rice yield prediction method of coupling crop model and remote sensing data according to claim 1, characterized in that: The input parameters include: effective accumulated temperature from seedling emergence to flowering, effective accumulated temperature from flowering to maturity, initial total dry matter weight, leaf area index at seedling emergence, specific leaf area at the crop development stage, leaf senescence index, initial utilization efficiency of single leaf light energy, maximum leaf photosynthesis rate, leaf CO2 maximum synthesis rate reduction factor, leaf growth assimilate conversion rate, storage organ growth assimilate conversion rate, root growth assimilate conversion rate, and stem growth assimilate conversion rate.

9. The rice yield prediction method of coupling crop model and remote sensing data according to claim 1, characterized in that: The cost function is: ， In the formula, is the cost function; x1 is the parameter to be optimized, is the expected value of the parameter to be optimized, B is the error covariance matrix of x1, i is the number of remote sensing observations, n2 is the total number of remote sensing observations, is the vegetation index value obtained by remote sensing inversion for the i-th time, for The error covariance matrix of is the LAI of the i-th remote sensing observation simulated by the WOFOST model. The superscript T represents the transpose of the matrix.

Images