Rice yield prediction method coupling crop model and remote sensing data

Abstract

The invention discloses a crop model and remote sensing data coupled rice yield prediction method, which comprises the steps of fusing a time increment and a space increment based on a vegetation index spatio-temporal data fusion algorithm and a constrained least square method to obtain a final combined increment, inputting the final combined increment into a trained LAI remote sensing estimation model, outputting a remote sensing inversion LAI index, and predicting the yield of rice. The method comprises the following steps: constructing a cost function according to a WOFOST model simulation LAI index and a remote sensing inversion LAI index, utilizing an SCE-UA algorithm to iteratively optimize WOFOST model to-be-optimized parameters by taking cost function minimization as a target to obtain optimal input parameters of the WOFOST model, inputting the optimal input parameters and other to-be-optimized parameters into the WOFOST model, and outputting a predicted rice yield result. The problem that the field size and the remote sensing size cannot be matched in space is solved.

Description

Technical Field

[0001] The invention relates to the technical field of rice yield prediction, and in particular to a rice yield prediction method by coupling 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 of crop growth, 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 product prediction. With the continuous development of remote sensing technology, high-frequency, multi-band, and multi-spatial resolution remote sensing products, with their characteristics of strong timeliness and wide monitoring range, form a good complementary relationship with the advantages of crop models such as strong mechanism and time continuity. At present, almost all crop models are established at the field scale, resulting in the inability to match the field scale with the remote sensing scale in space.

[0003] Therefore, a rice yield prediction method that couples crop models and remote sensing data was developed to solve the above problems. Summary of the invention

[0004] The present invention proposes a rice yield prediction method by coupling crop models and remote sensing data, so as to solve the problem that almost all current crop models in the prior art are established at the field scale, resulting in the inability to match the field scale with the remote sensing scale in space.

[0005] The present invention achieves the above-mentioned purpose through the following technical solutions:

[0006] The present invention provides a rice yield prediction method by coupling crop models and remote sensing data, comprising:

[0007] 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;

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

[0009] According to the normalized vegetation index data set, a time increment is obtained based on a linear unmixing analysis;

[0010] According to the normalized vegetation index data set, a spatial increment is obtained based on a thin plate spline interpolation method;

[0011] 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;

[0012] 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;

[0013] 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;

[0014] 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;

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

[0016] A cost function is constructed according to the LAI index simulated by the WOFOST model and the LAI index inverted by remote sensing. 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.

[0017] Furthermore, the calculation formulas for calculating the normalized difference vegetation index with high temporal resolution and low spatial resolution and high spatial resolution and low temporal resolution according to the multi-temporal remote sensing data are as follows: ,

[0018] In the formula, 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.

[0019] Further, according to the normalized vegetation index data set, a time increment is obtained based on a linear unmixing analysis, including: ,

[0020] ,

[0021] Time increment ,include:

[0022]

[0023] Where n is the number of low-resolution pixels, 1 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. The NDVI increment is 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 .

[0024] Further, according to the normalized vegetation index data set, a spatial increment is obtained based on a thin plate spline interpolation method, including:

[0025] ,

[0026] In the formula, represents the space increment, and They are respectively the thin plate spline interpolation method and and The interpolated value at time instant, Indicates low spatial resolution ( ) is the j-th high-resolution pixel in .

[0027] Furthermore, 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:

[0028] ,

[0029] ,

[0030] ,

[0031] 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 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.

[0032] Further, the remote sensing estimated vegetation index value at the LAI measured point is calculated according to the latitude and longitude of the LAI measured point and the multi-temporal remote sensing data, including:

[0033] 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;

[0034] Construct a 3*3 window with the pixel where the measured point is located as the center pixel;

[0035] 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.

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

[0037] ,

[0038] ,

[0039] 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.

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

[0041] Furthermore, the cost function is:

[0042] ,

[0043] 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.

[0044] The beneficial effects of the present invention are:

[0045] The present invention proposes a rice yield prediction method by coupling crop models with remote sensing data, which solves the problem that almost all current crop models are established at the field scale, resulting in the incompatibility between the field scale and the remote sensing scale in space. By coupling the crop model with remote sensing data, the prediction of rice yield at the field scale is realized. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0047] Figure 2 Generate a combined incremental technology roadmap based on the CLS algorithm in the embodiment of this application;

[0048] Figure 3 This is a flow chart of sensitivity analysis of WOFOST model parameters in the embodiment of this application;

[0049] Figure 4 This is a flow chart of rice yield prediction using the WOFOST model and remote sensing data assimilation in an embodiment of the present application. DETAILED DESCRIPTION

[0050] In order to make the purpose, 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 in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings here can be arranged and designed in various different configurations.

[0051] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0052] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, further definition and explanation thereof is not required in subsequent drawings.

[0053] The specific implementation modes of the present invention are described in detail below in conjunction with the accompanying drawings.

[0054] like Figure 4 As shown, a rice yield prediction method coupling crop model and remote sensing data includes:

[0055] 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;

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

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

[0058] According to the normalized vegetation index data set, a time increment is obtained based on a linear unmixing analysis;

[0059] According to the normalized vegetation index data set, a spatial increment is obtained based on a thin plate spline interpolation method;

[0060] 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;

[0061] 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;

[0062] 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;

[0063] 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;

[0064] Simulating the LAI index based on the WOFOST model according to the parameters to be optimized;

[0065] A cost function is constructed according to the LAI index simulated by the WOFOST model and the LAI index inverted by remote sensing. 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.

[0066] A rice yield prediction method coupling crop model and remote sensing data specifically includes:

[0067] Step 1: Remote sensing data preprocessing

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

[0069] 1. Radiation calibration

[0070] Through radiation calibration, the systematic error caused by the sensor itself and other factors in the remote sensing image data can be eliminated. The radiation calibration formula can be used to calculate DN to obtain the calibrated emissivity or reflectivity.

[0071] The formula for converting the DN value of remote sensing images into radiance using the calibration coefficient is:

[0072] ,

[0073] L is the equivalent radiance after conversion, in units of DN is the digital quantization output value of the sensor; gain and offset are gain and offset respectively, which can be obtained by querying the calibration information of the original remote sensing image.

[0074] 2. Atmospheric correction

[0075] Atmospheric correction is used to eliminate the influence of factors such as atmosphere and light on the emission of ground objects, and obtain real physical model parameters such as ground object reflectivity, emissivity, and surface temperature. These are used to eliminate the influence of water vapor, oxygen, etc. in the atmosphere on the reflection of ground objects, and eliminate the influence of atmospheric molecules and aerosol scattering.

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

[0077] 3. Geometric Correction

[0078] Image deformation caused by systematic and non-systematic factors is solved by geometric correction, so as to achieve geometric alignment with a standard image or map. The present invention uses a processed high-resolution image with accurate geometric position information in the same area as a reference image to correct the high-resolution image to be calibrated. If the geometric accuracy of the image is required to be high, the geometric correction can be completed by using ground control points and a geometric correction model to construct the geometric relationship between the image and the ground coordinates / image.

[0079] (IV) Inlay cutting

[0080] Remote sensing image mosaicking means that when the study area exceeds the scope covered by a single remote sensing image, two images are usually spliced ​​together to form one or a series of larger images covering the entire area; the purpose of cropping is to remove the area outside the project, and the image is usually cropped according to the administrative boundary or natural boundary. The present invention uses remote sensing processing software such as ARCGIS and ENVI to achieve mosaicking and cropping of remote sensing images.

[0081] Step 2: Generate high temporal and spatial resolution NDVI based on vegetation index fusion algorithm

[0082] 1. Calculation of medium and high resolution NDVI datasets

[0083] The medium-resolution NDVI data used in the present invention is calculated by 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 by the HLS surface reflectance product according to formula (1), with a temporal resolution of about 16 days and a spatial resolution of 30 meters.

[0084] (1)

[0085] 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.

[0086] (II) Obtaining time increment based on linear unmixing

[0087] According to the linear spectral mixing theory, the temporal variation of the NDVI of a low-resolution pixel can be considered as a linear combination of the NDVI increments of all high-resolution pixels within the pixel over a short period of time. Therefore, assuming that high-resolution pixels of the same land cover type have similar increments in a local area, a linear mixing model is used to convert the low-resolution pixels from the baseline date to the baseline date. Decompose into forecast dates , and its calculation formula is shown in (2) and (3):

[0088] (2)

[0089] (3)

[0090] Time increment ,include:

[0091]

[0092] 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 The fraction of land classes within low-resolution pixels (x, y), is the set of all low-resolution NDVI increments within the moving window; , and Respectively represent The minimum, maximum, and standard deviation of .

[0093] (III) Obtaining spatial increments based on thin plate spline interpolation

[0094] The thin plate spline interpolation method (TPS) is used to and The low-resolution NDVI at the moment is interpolated to high resolution. Then use and The difference between the interpolation results is used to obtain the spatial increment. Since this increment only uses the spatial dependency between low-resolution pixels, it is called spatial increment. , as shown in Announcement (4).

[0095] (4)

[0096] In the formula, and TPS in and The interpolated value at time instant, Represents the jth high-resolution pixel in the low spatial resolution (x, y).

[0097] 4. Combining time increment and space increment to predict high temporal and spatial resolution NDVI dataset

[0098] like Figure 2 As shown, the time increment ( ) and spatial increment ( ) can improve the performance and robustness of the fusion algorithm. and are the time increments ( ) and spatial increment ( ) The spatial increment and temporal increment after downscaling, is the increment of the low spatial resolution pixel, that is, the low resolution increment. The simplest and most effective combination method of the two increments is to add them together through a reasonable weight coefficient based on the least squares method. To further improve the accuracy of the combined increment, the residual needs to be assigned to each high resolution pixel within the low resolution pixel (x, y). As shown in formula (5):

[0099] (5)

[0100] In the formula, for High-resolution pixels at each moment ( ) of the NDVI value, Predicted by the algorithm NDVI value of high-resolution pixel at the moment; High resolution pixels The corresponding combined increment is calculated by formula (6); Low resolution pixels The residual within is calculated by formula (7), and the final combined increment is for:

[0101] (6)

[0102] (7)

[0103] In formula (6), and are the weight coefficients of space increment and time increment, and are the spatial increment and the time increment respectively. In formula (7), For low spatial resolution pixels arrive The real increment of NDVI. According to formula (5), High-resolution NDVI forecast at any time High-resolution NDVI at all times ( ).

[0104] Step 3: 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, and train and test the LAI remote sensing estimation model according to the training sample set and the test sample set, specifically including:

[0105] ARCGIS software was used to match the position of the point in the remote sensing index image based on the latitude and longitude information of the LAI measured point; then, a 3*3 window was constructed with the pixel where the measured point was located as the central pixel ( Figure 1 As shown in the figure, "O" is the pixel where the measured point is located, and the gray shadow is the 3*3 window); finally, the average value of the 9 pixels in the 3*3 window is calculated as the vegetation index value corresponding to the measured point. Then, the LAI measured value and the remote sensing vegetation index value are divided into a training sample set and a test sample set in a ratio of 7:3; the training sample set is mainly used to build the 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.

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

[0107] (8)

[0108] (9)

[0109] 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.

[0110] The final combined increment is input into the trained LAI remote sensing estimation model, and the remote sensing inversion LAI index is output. The WOFOST model is constructed according to the meteorological, soil and crop parameters and the remote sensing inversion LAI index.

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

[0112] like Figure 3As shown, the WOFOST model has more input parameters, and in the actual simulation process, we found that by changing the values ​​of different input parameters, there are certain differences in the degree of 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 parameter changes that have less 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 to improving the operation rate of the model. The present invention selects 26 input parameters related to rice growth from the numerous input parameters of the WOFOST model, uses python language to drive the model operation, and uses the Sobol index in the SALib module to perform a global sensitivity analysis on the WOFOST model input parameters to obtain first-order sensitivity, second-order sensitivity and total-order sensitivity. According to the first-order sensitivity, second-order sensitivity and total-order sensitivity, the parameters to be optimized are obtained. The parameters to be optimized and the main value ranges are shown in Table 1:

[0113] Table 1

[0114] 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]

[0115] The input parameters in Table 1 represent the following contents: TSUM1 and TSUM2 are the effective accumulated temperatures from seedling to flowering and from flowering to maturity in the WOFOST model, respectively; TDWI is the initial total dry matter weight; LAIEM is the leaf area index at seedling emergence; SLATB1, SLATB2, and SLATB3 are the specific leaf areas at the crop development stage (DVS) = 0.0, DVS = 0.5, and DVS = 2.0; SPAN is the leaf senescence index; EFFTB is the initial utilization efficiency of light energy per leaf, and EFFTB1, EFFTB2, and EFFTB3 are respectively EFFTB1, EFFTB2, and EFFTB3 at 0℃, 10℃, 20℃, 30℃, and 40℃. B3, EFFTB4, EFFTB5; AMAX is the maximum rate of leaf photosynthesis, and AMAX when DVS=0.0, DVS=1.0, DVS=1.3 and DVS=2.0 are named AMAX1, AMAX2, AMAX3 and AMAX4; TMPF under 0℃, 10℃, 20℃, 30℃ and 40℃ are recorded as TMPF1, TMPF2, TMPF3, TMPF4 and TMPF5; CVL is the assimilate conversion rate of leaf growth; CVO is the assimilate conversion rate of storage organ growth; CVR is the assimilate conversion rate of root growth; CVS is the assimilate conversion rate of stem growth.

[0116] The specific steps for sensitivity analysis of input parameters are as follows:

[0117] The present invention uses the Sobol index in the python SALib module to perform 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 contribution of input variables and their interactions to the uncertainty of model output. This method is suitable for complex nonlinear and high-dimensional systems.

[0118] The core idea of ​​the Sobol method is to decompose the total variance of the model output into the contribution of each input variable and the contribution of the interaction between the variables. The mathematical idea is that for the input variable X=( , , …, ), the variance of the model output Y=f(X) can be decomposed into:

[0119] (10)

[0120] In the formula, For variables Independent contribution of For variables and The interactive contribution of Contribute to higher-order interactions of all variables.

[0121] The first-order sensitivity index Represents input variables The contribution to the output result. The larger the value, the greater the input variable. The greater the impact on the output results.

[0122] (11)

[0123] Second-order interaction sensitivity index Represented as a variable and The contribution ratio of interaction to the output result.

[0124] (12) Total sensitivity index Represented as input variable The total contribution of the main effect of and the interaction effects of other variables, The higher the value, the more The higher the importance to the output result.

[0125] (13)

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

[0127] Based on the results of the global sensitivity analysis of Sobol input parameters, the input parameters with greater sensitivity were adjusted, and the remaining parameters with less sensitivity, parameters with less impact on the output results, and parameters that were not subjected to sensitivity analysis were obtained through field measurements, reference to the research results of others, or the default values ​​of the WOFOST model to obtain the parameters to be optimized, so that the WOFOST model can more accurately simulate the growth and development process of rice and predict yield.

[0128] Step 5: Rice yield prediction at the surface level by coupling crop models with remote sensing data

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

[0130] (14)

[0131] 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.

[0132] The present invention selects LAI as the optimization comparison object, uses the LAI estimated by remote sensing and the LAI data simulated by the WOFOST crop model to construct a cost function (Formula 14), and uses the SCE-UA algorithm to iteratively optimize the WOFOST model input parameters to ultimately minimize the cost function, thereby obtaining the optimal input parameters of the model. The optimal input parameters are then input into the WOFOST model together with other parameters to be optimized to simulate the yield results after data assimilation. The other input parameters represent other parameters to be optimized except the optimal input parameters.

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

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

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

[0136] 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, and the spatial increment is obtained by using the thin plate spline interpolation method. The two increments are combined using the vegetation index spatiotemporal data fusion algorithm, and the final combined increment is obtained through the constrained least squares method (CLS). A high spatiotemporal resolution (30 meters per day) NDVI product data set is generated, which effectively solves the problem of spatial matching between the field scale and the remote sensing scale.

[0137] 3. Use Python language to drive the model, 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 the WOFOST model parameters are "localized".

[0138] 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 assimilates remote sensing LAI data; and the applicability of the algorithm proposed in this invention is tested by comparing the rice yield results simulated by a single WOFOST model with the actual yield data.

[0139] 5. Coupling crop models and remote sensing data can make it possible to quantitatively predict crop growth and yield at the field scale.

[0140] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection 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