Salinization farmland crop yield prediction method

Abstract

The invention discloses a salinization farmland crop yield prediction method, which is applied to the technical field of crop yield monitoring and prediction, and comprises the following steps: based on a vegetation index and a salinity index, inverting salinity by using a partial least square regression method; inverting the soil moisture based on an empirical relation model of the conditional vegetation temperature index and the actually measured soil moisture; inverting the leaf area index based on an empirical relation model of the normalized vegetation index and the actually measured leaf area index; the method comprises the following steps: constructing an HYDRUS-1D-EPIC coupling model considering salinity, carrying out sensitivity analysis by using an extended Fourier amplitude sensitivity test method, carrying out parameter calibration and posterior distribution inference on the coupling model by using an MCMC method, and sampling posterior distribution by using a DREAM method; and assimilating the inverted salinity, soil moisture and leaf area index to the coupling model by using an ensemble Kalman filtering method, and predicting the crop yield of the salinized farmland. According to the method, the yield prediction precision of the salinized crops is effectively improved.

Description

Technical Field

[0001] The present invention relates to the technical field of crop yield monitoring and prediction, and in particular to a method for predicting crop yield in salinized farmland. Background Art

[0002] Monitoring and predicting crop yields is a key step in addressing food security issues on a global scale. Crop yields are affected by a combination of climate change, population growth, soil erosion, and natural variability in weather. Therefore, providing a method that can accurately and timely assess crop growth and yield is of great significance for improving the sustainability of agricultural food production.

[0003] In recent years, crop models based on physical processes have used multi-source data such as meteorology, soil, and management measures to quantitatively simulate the crop growth and development process and yield formation. Accurate crop yield prediction can provide agricultural managers with timely and reliable decision-making basis. Generally, crop models have good application effects at a single point scale, but due to the limitations of observation data at the regional scale, the simulation effect at the regional scale needs to be improved. Remote sensing technology can be used to obtain crop growth and land use information over a large area, and remote sensing data products are rich in types. Data assimilation algorithms can effectively integrate observation results into the simulation process to improve the accuracy of model simulation results. By assimilating crop models with observation data, the applicability and generalization ability of the model under different environmental conditions can be improved, and more reliable crop yield predictions can be obtained. Therefore, data assimilation technology is increasingly being used to improve the simulation ability of crop models at large spatial scales.

[0004] However, in the existing studies on the combination of data assimilation and crop models, most of the crop models used are WOFOST, DSSAT, AquaCrop and SAFY. In addition, most data assimilation studies focus only on a single state variable, and rarely explore the joint assimilation of multiple variables. In the joint assimilation studies, there are even fewer studies that use soil salinity as a data assimilation variable. In fact, soil salinization is a key factor restricting crop growth in arid areas. Timely and accurate monitoring and simulation of soil salt content is the basis and key to the management and improvement of salinized soils and the formulation of agricultural irrigation systems.

[0005] Therefore, how to provide a method for predicting crop yield in salinized farmland that can jointly assimilate multiple variables and use soil salinity as a data assimilation variable in the joint assimilation research, effectively improve the accuracy of coupling models in simulating crop growth and development in salinized areas, and thus effectively improve the accuracy of crop yield prediction in salinized areas is an urgent problem that technicians in this field need to solve. Summary of the invention

[0006] In view of this, the present invention proposes a method for predicting crop yield in salinized farmland.

[0007] In order to achieve the above object, the present invention adopts the following technical solution: A method for predicting crop yield in salinized farmland, comprising: Step 1: Based on the vegetation index and salt index, the partial least squares regression method is used to invert the salt content; the soil moisture content is inverted based on the empirical relationship model between the conditional vegetation temperature index and the measured soil moisture content; the leaf area index is inverted based on the empirical relationship model between the normalized difference vegetation index and the measured leaf area index; Step 2: Construct the HYDRUS-1D-EPIC coupling model considering salt content, and use the extended Fourier amplitude sensitivity test method to perform sensitivity analysis, use the MCMC method to calibrate the parameters of the coupling model and infer the posterior distribution, and use the DREAM method to sample the posterior distribution; Step 3: The inverted salt, soil moisture, and leaf area index are assimilated into the HYDRUS-1D-EPIC coupled model using the ensemble Kalman filter method to predict crop yields in salinized farmland.

[0008] Optionally, in step 1, the vegetation index includes: normalized vegetation index, enhanced vegetation index, soil adjusted vegetation index, ratio vegetation index, and difference vegetation index.

[0009] Optionally, in step 1, the salt index includes: ; in, is the blue band reflectivity; is the red band reflectivity; is the reflectivity in the near-infrared band; is the green band reflectivity.

[0010] Optionally, in step 1, the determination coefficient, root mean square error and mean absolute error are used to determine the inversion accuracy of the partial least squares regression method; The coefficient of determination is as follows: The root mean square error is as follows: The mean absolute error is as follows: in, is the coefficient of determination; is the root mean square error; is the mean absolute error; is the sample size; is the i-th observation value; is the i-th observation value predicted by the model; for The sample mean of .

[0011] Optionally, in step 1, the calculation of the conditional vegetation temperature index is as follows: in, is the conditional vegetation temperature index; is the normalized vegetation index of the i-th day; for The maximum LST value of all pixels when it is equal to a certain value; The NDVI of a pixel is LST at ; for The minimum LST value of all pixels when it is equal to a certain value; and are the coefficients obtained by fitting the dry edge data points; and are the coefficients obtained by fitting the wet edge data points.

[0012] Optionally, in step 1, the normalized vegetation index is calculated as follows: in, is the normalized difference vegetation index; is the red band reflectivity; is the reflectivity in the near-infrared band.

[0013] Optionally, in step 2, a HYDRUS-1D-EPIC coupling model considering salinity is constructed, specifically: Water movement, as follows: The one-dimensional Richards equation is used to describe the water movement: in, is the volumetric moisture content; For time; is the pressure head; is the spatial coordinate, with the upward orientation being positive; is the soil hydraulic conductivity; is the root source and sink item; The crop roots absorb water as follows: Soil water absorbed by crop roots from unit volume of soil per unit time: in, is the soil moisture stress coefficient; is the soil salt stress coefficient; is the potential root water absorption rate at a certain depth; Soil moisture stress coefficient: in, It is the anaerobic point for crops; The optimum moisture point for crops; It is a water stress point; It is the point where crops wither; Soil salt stress coefficient: in, is the salt concentration of soil solution; It is the salt concentration of the soil solution corresponding to the reduction of root water absorption rate by 50%; is an empirical constant; Potential root water uptake rate at a certain depth: in, is the normal root water absorption distribution term; is the root system depth; is the soil layer depth; Soil hydraulic characteristic parameters are as follows: in, is the saturated moisture content; is the residual moisture content; is the saturated hydraulic conductivity; is the fitting parameter; is the effective saturation; is the derivative of the air intake value; is the pore size distribution index; The potential transpiration rate of crops is as follows: in, is the solar radiation extinction coefficient; is the crop leaf area index; is the reference crop evapotranspiration; The convection-diffusion equation is used to describe the salt migration and transformation process in the soil, as follows: in, is the salt concentration of soil solution; is the dry bulk density of soil; is the soil water flux; is the amount of solute adsorbed by soil particles; is the hydrodynamic diffusion coefficient of soil solute; in, is the adsorption distribution coefficient; is the radial dispersion; is the flow rate; is the molecular diffusion coefficient of free water; is the distortion coefficient in the liquid phase; The EPIC model calculates and simulates production on a daily basis; The accumulated temperature is calculated as follows: in, is the accumulated temperature on the i-th day; , are the maximum and minimum temperatures on the i-th day, respectively; It is the base temperature for crop growth; The heat unit coefficient is used to indicate the stage of crop growth and development: in, is the heat unit coefficient of the i-th day; The maximum heat unit required for crops to mature; Leaf area changes are as follows: in, is the leaf area index on the i-th day; is the minimum crop growth stress factor value; is the maximum leaf area index; , is the crop parameter; From the beginning of the decline in leaf area to the end of the growing season : in, is the parameter of LAI attenuation rate; is the actual maximum leaf area index; is the heat unit coefficient corresponding to the decrease of leaf area index; The potential biomass growth is as follows: Beer's law is used to calculate the solar radiation received by crops: in, is the photosynthetically active radiation received on day i; is the total solar radiation on the i-th day; in, is the potential daily biomass growth; Conversion factors for crops converting energy into biomass; Plant height, as follows: in, is the crop height on the i-th day; is the maximum plant height; Root growth, as follows: Maximum daily root weight change: Root Depth: in, is the change in root depth on the i-th day; is the root depth on day i; is the maximum root depth; Crop yields are as follows: Calculate crop yield using harvest index: in, is the grain yield; To accumulate aboveground biomass; is the harvest index; in, is the harvest index on the i-th day; is the heat unit factor affecting the harvest index; in, It is the actual harvest index of crops; It is the sensitivity index of crops to drought; Actual crop yield: ; It is the actual crop yield after considering environmental stress; Environmental stress, including: Water and salt stress: in, is the water-salt stress coefficient; is the crop potential transpiration rate; is the actual transpiration rate of the crop; Temperature stress: in, is the plant temperature stress factor; is the daily average temperature; The optimal temperature for crop growth; It is the basic temperature for crop growth; Effects of environmental stress on crop growth: in, is the actual biomass; Environmental impact factors; Optionally, in step 2, the DREAM method is used to sample the posterior distribution, specifically: Before sampling using the DREAM method, the parameter prior distribution and likelihood function are clarified to calculate the posterior distribution probability density , and set: the total number of chains required for sampling ; Maximum total number of sampling times before sampling convergence , , is the maximum number of sampling times of a single chain before convergence; the total number of continued sampling times after sampling convergence , , is the number of single chain sampling after convergence; Randomly select initial sampling points for each chain from the prior distribution to form an initial sampling set , and randomly select one from all Markov chains of the current state, and update the state through mutation and crossover of the DREAM method; Calculate the likelihood function of each candidate state, and calculate the posterior probability distribution of each candidate state based on the likelihood function; the posterior probability distribution is used to compare the relative probabilities of different states to determine whether to accept the state; Determine whether to accept the proposed point; Calculating the acceptance rate ,as follows: in, is the posterior probability of the t+1th iteration of the i-th chain parameter; To express a proposal; is the posterior probability of the tth iteration of the parameters of the i-th chain; is the tth sampling result of the i-th chain; From uniform distribution Random sampling, assuming the value is , and the acceptance rate and For comparison, if , then accept As the t+1th sampling result of the i-th chain, otherwise, The value of is taken as the t+1th sampling result of the i-th chain; Repeat the above steps until all N chains have completed the t+1th sampling; When all chains have completed the t+1th sampling, it is determined whether the sampling has converged; calculate The statistics are as follows: in, , , The parameters are The intra-chain and inter-chain variances of the j-th dimension, as well as the posterior estimated standard deviation; is the rounding operator; is the value of the jth dimension of the nth parameter sample taken by the current i-th chain; is the mean value of the jth dimension of all parameter samples taken by the current i-th chain; is the mean value of the jth dimension of all parameter samples taken by all current chains; for The jth term of the statistic; Determine whether the sampling is convergent. If If each component of is less than or equal to 1.2, the sampling reaches convergence until the maximum total number of sampling times before the sampling convergence is reached.

[0014] Optionally, in step 3, the inverted salinity, soil moisture and leaf area index are assimilated into the HYDRUS-1D-EPIC coupled model using the ensemble Kalman filter method, specifically: Initialize the state variables of the model simulation and generate the initial simulation set ; Generate observation set and forecast set: When the model runs to time t-1, there are observations of state variables at time t, then add a mean of 0 and a covariance matrix of Gaussian perturbation of , generating an observation set consistent with the dimension of the model state variable simulation set ; Continue to run the model simulation set for one time step to generate the forecast set at time t ,as follows: in, It is a prediction of the state of the coupled model at the next moment based on the simulation state at the current moment; is the forecast set at time t; is the simulation result of the jth element of the simulation set at time t-1; is the forecast set corresponding to the j-th element at time t; Update the forecast set based on the observation error and forecast error as follows: Further simplified to:

[0015] in, is the updated result of the jth element in the forecast set at time t; is the Kalman gain; For Corresponding observations; is the identity matrix; is the forecast error; is the observation error; To analyze the error, it is used to calculate the error of the updated state variables and quantify the uncertainty of the model; Replace the model simulation set with the updated forecast set as follows: in, is the model simulation set; is the forecast set after the model is updated; is the updated result of the jth element in the forecast set at time t; Continue running the model until a new observation appears at the next moment, then repeat the above generation, update, replacement, and running operations until the model simulation ends.

[0016] It can be seen from the above technical solutions that, compared with the prior art, the present invention proposes a method for predicting crop yield in salinized farmland. By constructing a HYDRUS-1D-EPIC coupling model that takes salt into consideration, and using the extended Fourier amplitude sensitivity test method to perform sensitivity analysis, the MCMC method is used to calibrate the parameters of the coupling model and infer the posterior distribution, the DREAM method is used to sample the posterior distribution, and the ensemble Kalman filter method is used to assimilate the inverted salt, soil moisture, and leaf area index into the HYDRUS-1D-EPIC coupling model to predict the crop yield in salinized farmland, realizing the joint assimilation of multiple variables, and in the study of joint assimilation, soil salinity is used as a data assimilation variable, which effectively improves the accuracy of the coupling model in simulating crop growth and development in salinized areas, thereby effectively improving the prediction accuracy of the salinized crop yield of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0018] Figure 1 It is a schematic diagram of the method flow of the present invention. DETAILED DESCRIPTION

[0019] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. 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.

[0020] Embodiment 1: Embodiment 1 of the present invention discloses a method for predicting crop yield in salinized farmland, such as Figure 1 As shown, including: Step 1: Based on the vegetation index and salt index, the partial least squares regression method is used to invert the salt content; the soil moisture is inverted based on the empirical relationship model between the conditional vegetation temperature index and the measured soil moisture; the leaf area index is inverted based on the empirical relationship model between the normalized difference vegetation index and the measured leaf area index.

[0021] Using remote sensing to invert salinity, Sentinel-1 and Sentinel-2 were selected as the remote sensing data sources. The spectrum consists of four bands: blue, green, red, and near infrared.

[0022] Vegetation index, including: Normalized Difference Vegetation Index, Enhanced Vegetation Index, Soil Adjusted Vegetation Index, Ratio Vegetation Index, and Difference Vegetation Index.

[0023] Salt index, including: ; in, is the blue band reflectivity; is the red band reflectivity; is the reflectivity in the near-infrared band; is the green band reflectivity.

[0024] Partial least squares regression (PLSR) is a statistical method for building regression models, which is particularly suitable for situations where the independent variables are highly correlated (multicollinearity) or the number of independent variables is greater than the number of observations. Therefore, the present invention chooses to use the PLSR model for inversion. The PLSR model is established using MATLAB software. PLSR introduces data dimensionality reduction, information synthesis and screening technology in the modeling process, and improves the explanatory power of the model by extracting new comprehensive components that can best explain the data system.

[0025] Using the coefficient of determination , Root Mean Square Error And the mean absolute error To judge the inversion accuracy of partial least squares regression method; The coefficient of determination is as follows: The root mean square error is as follows: The mean absolute error is as follows: in, is the coefficient of determination; is the root mean square error; is the mean absolute error; is the sample size; is the i-th observation value; is the i-th observation value predicted by the model; for The sample mean of .

[0026] Optionally, in step 1, the calculation of the conditional vegetation temperature index is as follows: in, is the conditional vegetation temperature index; is the normalized vegetation index of the i-th day; for The maximum LST value of all pixels when it is equal to a certain value; The NDVI of a pixel is LST at ; for The minimum LST value of all pixels when it is equal to a certain value; and are the coefficients obtained by fitting the dry edge data points; and are the coefficients obtained by fitting the wet edge data points.

[0027] The VTCI inverted from Landsat data is a relative dryness and wetness index in theoretical terms. The smaller the VTCI, the more severe the drought, the less soil moisture, and the more severe the water stress on crops; the larger the VTCI, the more soil moisture, the less severe the drought or even no drought. Regression analysis was performed on the VTCI inverted from Landsat data and the measured soil volume moisture content during the same period to estimate the soil moisture content.

[0028] The calculation of the normalized vegetation index is as follows: in, is the normalized difference vegetation index; is the red band reflectivity; is the reflectivity in the near-infrared band.

[0029] The vegetation index is calculated using the multispectral reflectance data of the Sentinel-2 satellite and the leaf area index data of the corresponding spatial resolution is calculated in combination with the empirical model. When establishing the empirical relationship between NDVI and LAI based on the measured data, three functional relationships, linear, power exponential and exponential, are tried, and the relationship with the best accuracy is selected for the calculation of LAI.

[0030] Step 2: Construct the HYDRUS-1D-EPIC coupling model considering salt content, and use the extended Fourier amplitude sensitivity test method to perform sensitivity analysis, use the MCMC method to calibrate the parameters of the coupling model and infer the posterior distribution, and use the DREAM method to sample the posterior distribution.

[0031] Construct the HYDRUS-1D-EPIC coupling model considering salt, specifically: Water movement, as follows: The one-dimensional Richards equation is used to describe the water movement: in, is the volumetric moisture content; For time; is the pressure head; is the spatial coordinate, with the upward orientation being positive; is the soil hydraulic conductivity; is the root source and sink item; The crop roots absorb water as follows: Soil water absorbed by crop roots from unit volume of soil per unit time: in, is the soil moisture stress coefficient; is the soil salt stress coefficient; is the potential root water absorption rate at a certain depth; Soil moisture stress coefficient: in, It is the anaerobic point for crops; The optimum moisture point for crops; It is a water stress point; It is the point where crops wither; Soil salt stress coefficient: in, is the salt concentration of soil solution; It is the salt concentration of the soil solution corresponding to the reduction of root water absorption rate by 50%; is an empirical constant; Potential root water uptake rate at a certain depth: in, is the normal root water absorption distribution term; is the root system depth; is the soil layer depth; Soil hydraulic characteristic parameters are as follows: in, is the saturated moisture content; is the residual moisture content; is the saturated hydraulic conductivity; is the fitting parameter; is the effective saturation; is the derivative of the air intake value; is the pore size distribution index; The potential transpiration rate of crops is as follows: in, is the solar radiation extinction coefficient; is the crop leaf area index; is the reference crop evapotranspiration; The convection-diffusion equation is used to describe the salt migration and transformation process in the soil, as follows: in, is the salt concentration of soil solution; is the dry bulk density of soil; is the soil water flux; is the amount of solute adsorbed by soil particles; is the hydrodynamic diffusion coefficient of soil solute; in, is the adsorption distribution coefficient; is the radial dispersion; is the flow rate; is the molecular diffusion coefficient of free water; is the distortion coefficient in the liquid phase; The EPIC model calculates and simulates production on a daily basis; The accumulated temperature is calculated as follows: in, is the accumulated temperature on the i-th day; , are the maximum and minimum temperatures on the i-th day, respectively; It is the base temperature for crop growth; The heat unit coefficient is used to indicate the stage of crop growth and development: in, is the heat unit coefficient of the i-th day; The maximum heat unit required for crops to mature; Leaf area changes are as follows: in, is the leaf area index on the i-th day; is the minimum crop growth stress factor value; is the maximum leaf area index; , is the crop parameter; From the beginning of the decline in leaf area to the end of the growing season : in, is the parameter of LAI attenuation rate; is the actual maximum leaf area index; is the heat unit coefficient corresponding to the decrease of leaf area index; The potential biomass growth is as follows: Beer's law is used to calculate the solar radiation received by crops: in, is the photosynthetically active radiation received on day i; is the total solar radiation on the i-th day; in, is the potential daily biomass growth; Conversion factors for crops converting energy into biomass; Plant height, as follows: in, is the crop height on the i-th day; is the maximum plant height; Root growth, as follows: Maximum daily root weight change: Root Depth: in, is the change in root depth on the i-th day; is the root depth on day i; is the maximum root depth; Crop yields are as follows: Calculate crop yield using harvest index: in, is the grain yield; To accumulate aboveground biomass; is the harvest index; in, is the harvest index on the i-th day; is the heat unit factor affecting the harvest index; in, It is the actual harvest index of crops; It is the sensitivity index of crops to drought; Actual crop yield: ; It is the actual crop yield after considering environmental stress; Environmental stress, including: Water and salt stress: in, is the water-salt stress coefficient; is the crop potential transpiration rate; is the actual transpiration rate of the crop; Temperature stress: in, is the plant temperature stress factor; is the daily average temperature; The optimal temperature for crop growth; It is the basic temperature for crop growth; Effects of environmental stress on crop growth: in, is the actual biomass; Environmental impact factors; The model simulation step is on a daily time scale, with meteorological data (temperature, wind speed, radiation, sunshine time, relative humidity, rainfall, etc.) as input and soil moisture, soil salinity, leaf area, yield, etc. as output.

[0032] Extended Fourier amplitude sensitivity test (Extended FAST) decomposes the variance of model parameters to model output variables and divides the sensitivity of parameters into two categories: the sensitivity of a single parameter to the model output variable, measured by the first-order sensitivity index, and the sensitivity of the interaction between parameters to the model output variable, measured by the total sensitivity index. This process is performed using Simlab software. The sensitivity analysis is performed using the extended Fourier amplitude sensitivity test method, specifically: Select the main parameters of the model and give the initial values ​​of the parameters; Set the selected parameters to fluctuate by ±30% and set them to be evenly distributed; Simlab software was used to sample each parameter based on the Extended FAST method and set the number of samples; Input the parameter group generated by Simlab software into the coupling model to drive the model to perform Monte Carlo simulation; Use the model to output the variables each time; The generated variables are input into the software, and the software uses the Extended FAST method to solve the parameter sensitivity index corresponding to different output variables.

[0033] The purpose of this step is to determine the main parameters that affect the assimilation variables, because only a few parameters have a decisive effect on the output of the model.

[0034] The MCMC method is used to calibrate the parameters of the coupled model and infer the posterior distribution. The MCMC method is a numerical computing technique for sampling from complex probability distributions. The MCMC method constructs a Markov chain to make it wander in the state space of the probability distribution and eventually converge to the target distribution. This method is very useful in statistics, machine learning, Bayesian inference and other fields. The main idea of ​​the MCMC method is to construct a Markov chain by defining a transition core that describes how to transition from the current state to the next state. This core needs to satisfy the Markov property, that is, the probability of the next state depends only on the current state, not on the past state. Then, by repeatedly applying the transition core in the state space, a state sequence can be generated, which will eventually converge to the target probability distribution. The converged state sequence can be used to estimate the statistical properties of the probability distribution, such as the mean, variance, etc.

[0035] The DREAM method is used to sample the posterior distribution, specifically: Before sampling using the DREAM method, the parameter prior distribution and likelihood function are clarified to calculate the posterior distribution probability density , and set: the total number of chains required for sampling ; Maximum total number of sampling times before sampling convergence , , is the maximum number of sampling times of a single chain before convergence; the total number of continued sampling times after sampling convergence , , is the number of single chain sampling after convergence; Randomly select initial sampling points for each chain from the prior distribution to form an initial sampling set , and randomly select one from all Markov chains of the current state, and update the state through mutation and crossover of the DREAM method; Calculate the likelihood function of each candidate state, and calculate the posterior probability distribution of each candidate state based on the likelihood function; the posterior probability distribution is used to compare the relative probabilities of different states to determine whether to accept the state; Determine whether to accept the proposed point; Calculating the acceptance rate ,as follows: in, is the posterior probability of the t+1th iteration of the i-th chain parameter; To express a proposal; is the posterior probability of the tth iteration of the parameters of the i-th chain; is the tth sampling result of the i-th chain; From uniform distribution Random sampling, assuming the value is , and the acceptance rate and For comparison, if , then accept As the t+1th sampling result of the i-th chain, otherwise, The value of is taken as the t+1th sampling result of the i-th chain; Repeat the above steps until all N chains have completed the t+1th sampling; When all chains have completed the t+1th sampling, it is determined whether the sampling has converged; calculate The statistics are as follows: in, , , The parameters are The intra-chain and inter-chain variances of the j-th dimension, as well as the posterior estimated standard deviation; is the rounding operator; is the value of the jth dimension of the nth parameter sample taken by the current i-th chain; is the mean value of the jth dimension of all parameter samples taken by the current i-th chain; is the mean value of the jth dimension of all parameter samples taken by all current chains; for The jth term of the statistic; Determine whether the sampling is convergent. If If each component of is less than or equal to 1.2, the sampling reaches convergence until the maximum total number of sampling times before the sampling convergence is reached.

[0036] Step 3: The inverted salt, soil moisture, and leaf area index are assimilated into the HYDRUS-1D-EPIC coupled model using the ensemble Kalman filter method to predict crop yields in salinized farmland.

[0037] Ensemble Kalman Filter (EnKF) is developed on the basis of Kalman Filter. It introduces finite sets to approximate the probability distribution of state variables, thus adapting to the dynamic model and the Bayesian filtering data assimilation method when the observation operator is a nonlinear model.

[0038] The inverted salt, soil moisture and leaf area index are assimilated into the HYDRUS-1D-EPIC coupled model using the ensemble Kalman filter method, specifically: Initialize the state variables of the model simulation and generate the initial simulation set ; Generate observation set and forecast set: When the model runs to time t-1, there are observations of state variables at time t, then add a mean of 0 and a covariance matrix of Gaussian perturbation of , generating an observation set consistent with the dimension of the model state variable simulation set ; Continue to run the model simulation set for one time step to generate the forecast set at time t ,as follows: in, It is a prediction of the state of the coupled model at the next moment based on the simulation state at the current moment; is the forecast set at time t; is the simulation result of the jth element of the simulation set at time t-1; is the forecast set corresponding to the j-th element at time t; Update the forecast set based on the observation error and forecast error as follows: Further simplified to:

[0039] in, is the updated result of the jth element in the forecast set at time t; is the Kalman gain; For Corresponding observations; is the identity matrix; is the forecast error; is the observation error; To analyze the error, it is used to calculate the error of the updated state variables and quantify the uncertainty of the model; Replace the model simulation set with the updated forecast set as follows: in, is the model simulation set; is the forecast set after the model is updated; is the updated result of the jth element in the forecast set at time t; Continue running the model until a new observation appears at the next moment, then repeat the above generation, update, replacement, and running operations until the model simulation ends.

[0040] The embodiment of the present invention discloses a method for predicting crop yield in salinized farmland. By constructing a HYDRUS-1D-EPIC coupling model that takes salt into consideration, and using the extended Fourier amplitude sensitivity test method to perform sensitivity analysis, using the MCMC method to calibrate the parameters of the coupling model and infer the posterior distribution, using the DREAM method to sample the posterior distribution, and using the ensemble Kalman filter method to assimilate the inverted salt, soil moisture, and leaf area index into the HYDRUS-1D-EPIC coupling model, the yield of crops in salinized farmland is predicted, and the joint assimilation of multiple variables is achieved. In the study of joint assimilation, soil salinity is used as a data assimilation variable, which effectively improves the accuracy of the coupling model in simulating crop growth and development in salinized areas, thereby effectively improving the prediction accuracy of the salinized crop yield of the present invention.

[0041] In this specification, each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0042] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for predicting crop yield in salinized farmland, characterized in that: include: Step 1: Based on the vegetation index and salt index, the partial least squares regression method is used to invert the salt content; Soil moisture is inverted based on the empirical relationship model between the conditional vegetation temperature index and the measured soil moisture; leaf area index is inverted based on the empirical relationship model between the normalized difference vegetation index and the measured leaf area index; Step 2: Construct the HYDRUS-1D-EPIC coupling model considering salt content, and use the extended Fourier amplitude sensitivity test method to perform sensitivity analysis, use the MCMC method to calibrate the parameters of the coupling model and infer the posterior distribution, and use the DREAM method to sample the posterior distribution; Step 3: Use the ensemble Kalman filter method to assimilate the inverted salt, soil moisture and leaf area index into the HYDRUS-1D-EPIC coupling model to predict the crop yield in salinized farmland.

2. The method for predicting crop yield in salinized farmland according to claim 1, characterized in that: In step 1, the vegetation index includes: normalized vegetation index, enhanced vegetation index, soil adjusted vegetation index, ratio vegetation index, and difference vegetation index.

3. The method for predicting crop yield in salinized farmland according to claim 1, characterized in that: In step 1, the salt index includes: ; in, is the blue band reflectivity; is the red band reflectivity; is the reflectivity in the near-infrared band; is the green band reflectivity.

4. The method for predicting crop yield in salinized farmland according to claim 1, characterized in that: In step 1, the determination coefficient, root mean square error and mean absolute error are used to determine the inversion accuracy of the partial least squares regression method; The coefficient of determination is as follows: The root mean square error is as follows: The mean absolute error is as follows: in, is the coefficient of determination; is the root mean square error; is the mean absolute error; is the sample size; is the i-th observation value; is the i-th observation value predicted by the model; for The sample mean of .

5. The method for predicting crop yield in salinized farmland according to claim 1, characterized in that: In step 1, the calculation of the conditional vegetation temperature index is as follows: in, is the conditional vegetation temperature index; is the normalized vegetation index of the i-th day; for The maximum LST value of all pixels when it is equal to a certain value; The NDVI of a pixel is LST at ; for The minimum LST value of all pixels when it is equal to a certain value; and are the coefficients obtained by fitting the dry edge data points; and are the coefficients obtained by fitting the wet edge data points.

6. The method for predicting crop yield in salinized farmland according to claim 1, characterized in that: In step 1, the calculation of the normalized vegetation index is as follows: in, is the normalized difference vegetation index; is the red band reflectivity; is the reflectivity in the near-infrared band.

7. The method for predicting crop yield in salinized farmland according to claim 1, characterized in that: In step 2, a HYDRUS-1D-EPIC coupling model considering salinity is constructed, specifically: Water movement, as follows: The one-dimensional Richards equation is used to describe the water movement: in, is the volumetric moisture content; For time; is the pressure head; is the spatial coordinate, with the upward orientation being positive; is the soil hydraulic conductivity; is the root source and sink item; The crop roots absorb water as follows: Soil water absorbed by crop roots from unit volume of soil per unit time: in, is the soil moisture stress coefficient; is the soil salt stress coefficient; is the potential root water absorption rate at a certain depth; Soil moisture stress coefficient: in, It is the anaerobic point for crops; The optimum moisture point for crops; It is a water stress point; It is the point where crops wither; Soil salt stress coefficient: in, is the salt concentration of soil solution; It is the salt concentration of the soil solution corresponding to the reduction of root water absorption rate by 50%; is an empirical constant; Potential root water uptake rate at a certain depth: in, is the normal root water absorption distribution term; is the root system depth; is the soil depth; Soil hydraulic characteristic parameters are as follows: in, is the saturated moisture content; is the residual moisture content; is the saturated hydraulic conductivity; is the fitting parameter; is the effective saturation; is the derivative of the air intake value; is the pore size distribution index; The potential transpiration rate of crops is as follows: in, is the solar radiation extinction coefficient; is the crop leaf area index; is the reference crop evapotranspiration; The convection-diffusion equation is used to describe the salt migration and transformation process in the soil, as follows: in, is the salt concentration of soil solution; is the soil dry bulk density; is the soil water flux; is the amount of solute adsorbed by soil particles; is the hydrodynamic diffusion coefficient of soil solute; in, is the adsorption distribution coefficient; is the radial dispersion; is the flow rate; is the molecular diffusion coefficient of free water; is the distortion coefficient in the liquid phase; The EPIC model calculates and simulates production on a daily basis; The accumulated temperature is calculated as follows: in, is the accumulated temperature on the i-th day; , are the maximum and minimum temperatures on the i-th day, respectively; It is the base temperature for crop growth; The heat unit coefficient is used to indicate the stage of crop growth and development: in, is the heat unit coefficient of the i-th day; The maximum heat unit required for crops to mature; Leaf area changes are as follows: in, is the leaf area index on the i-th day; is the minimum crop growth stress factor value; is the maximum leaf area index; , is the crop parameter; From the beginning of the decline in leaf area to the end of the growing season : in, is the parameter of LAI attenuation rate; is the actual maximum leaf area index; is the heat unit coefficient corresponding to the decrease of leaf area index; The potential biomass growth is as follows: Beer's law is used to calculate the solar radiation received by crops: in, is the photosynthetically active radiation received on day i; is the total solar radiation on the i-th day; in, is the potential daily biomass growth; Conversion factors for crops converting energy into biomass; Plant height, as follows: in, is the crop height on the i-th day; is the maximum plant height; Root growth, as follows: Maximum daily root weight change: Root Depth: in, is the change in root depth on the i-th day; is the root depth on day i; is the maximum root depth; Crop yields are as follows: Calculate crop yield using harvest index: in, is the grain yield; To accumulate aboveground biomass; is the harvest index; in, is the harvest index on the i-th day; is the heat unit factor affecting the harvest index; in, It is the actual harvest index of crops; It is the sensitivity index of crops to drought; Actual crop yield: ; It is the actual crop yield after considering environmental stress; Environmental stresses, including: Water and salt stress: in, is the water-salt stress coefficient; is the crop potential transpiration rate; is the actual transpiration rate of the crop; Temperature stress: in, is the plant temperature stress factor; is the daily average temperature; The optimal temperature for crop growth; It is the basic temperature for crop growth; Effects of environmental stress on crop growth: in, is the actual biomass; Environmental impact factors; 8. The method for predicting crop yield in salinized farmland according to claim 1, characterized in that: In step 2, the DREAM method is used to sample the posterior distribution, specifically: Before sampling using the DREAM method, the parameter prior distribution and likelihood function are clarified to calculate the posterior distribution probability density. , and set: the total number of chains required for sampling ; Maximum total number of sampling times before sampling convergence , , is the maximum number of sampling times of a single chain before convergence; the total number of continued sampling times after sampling convergence , , is the number of single chain sampling after convergence; Randomly select initial sampling points for each chain from the prior distribution to form an initial sampling set , and randomly select one from all Markov chains of the current state, and update the state through mutation and crossover of the DREAM method; Calculate the likelihood function of each candidate state, and calculate the posterior probability distribution of each candidate state based on the likelihood function; the posterior probability distribution is used to compare the relative probabilities of different states to determine whether to accept the state; Determine whether to accept the proposed point; Calculating the acceptance rate ,as follows: in, is the posterior probability of the t+1th iteration of the i-th chain parameter; To express a proposal; is the posterior probability of the tth iteration of the parameters of the i-th chain; is the tth sampling result of the i-th chain; From uniform distribution Random sampling, assuming the value is , and the acceptance rate and For comparison, if , then accept As the t+1th sampling result of the i-th chain, otherwise, The value of is taken as the t+1th sampling result of the i-th chain; Repeat the above steps until all N chains have completed the t+1th sampling; When all chains have completed the t+1th sampling, it is determined whether the sampling has converged; calculate The statistics are as follows: in, , , The parameters are The intra-chain and inter-chain variances of the j-th dimension, as well as the posterior estimated standard deviation; is the rounding operator; is the value of the jth dimension of the nth parameter sample taken by the current i-th chain; is the mean value of the jth dimension of all parameter samples taken by the current i-th chain; is the mean value of the jth dimension of all parameter samples taken by all current chains; for The jth term of the statistic; Determine whether the sampling is convergent. If If each component of is less than or equal to 1.2, the sampling reaches convergence until the maximum total number of sampling times before the sampling convergence is reached.

9. The method for predicting crop yield in salinized farmland according to claim 1, characterized in that: In step 3, the inverted salinity, soil moisture and leaf area index are assimilated into the HYDRUS-1D-EPIC coupling model using the ensemble Kalman filter method, specifically: Initialize the state variables of the model simulation and generate the initial simulation set ; Generate observation set and forecast set: When the model runs to time t-1, there are observations of state variables at time t, then add a mean of 0 and a covariance matrix of Gaussian perturbation of , generating an observation set consistent with the dimension of the model state variable simulation set ; Continue to run the model simulation set for one time step to generate the forecast set at time t ,as follows: in, It is a prediction of the state of the coupled model at the next moment based on the simulation state at the current moment; is the forecast set at time t; is the simulation result of the jth element of the simulation set at time t-1; is the forecast set corresponding to the j-th element at time t; Update the forecast set based on the observation error and forecast error as follows: Further simplified to: in, is the updated result of the jth element in the forecast set at time t; is the Kalman gain; For Corresponding observations; is the identity matrix; is the forecast error; is the observation error; To analyze the error, it is used to calculate the error of the updated state variables and quantify the uncertainty of the model; Replace the model simulation set with the updated forecast set as follows: in, is the model simulation set; is the forecast set after the model is updated; is the updated result of the jth element in the forecast set at time t; Continue running the model until a new observation appears at the next moment, then repeat the above generation, update, replacement, and running operations until the model simulation ends.

Images