Dynamic crop moisture condition estimation method for assimilating unmanned aerial vehicle observation and crop model

Abstract

The invention discloses a dynamic crop moisture condition estimation method for assimilating unmanned aerial vehicle (UAV) observation and a crop model, which relates to the technical field of agriculture and comprises the following steps of: inverting key crop morphological parameters through a machine learning algorithm by using UAV remote sensing data with high temporal-spatial resolution; and assimilating the inverted key crop morphological parameters into a crop growth model to optimize a model simulation process, further constructing a CWS diagnosis model and generating a customized irrigation scheme, assimilating LAI and AGB data obtained by UAV by adopting an EnKF method, remarkably reducing simulation errors of LAI and AGB by data assimilation, verifying feasibility of an assimilation framework, realizing complementary advantages of the method, and improving the accuracy of the method. On one hand, a true value is provided for correcting model deviation, initial conditions of modeling and parameter uncertainty are reduced, on the other hand, the capability of the model for analyzing a crop growth and development mechanism is enhanced, and the DA framework has regional popularization potential in combination with the simplicity of the SAFYE model.

Description

Technical Field

[0001] The present invention relates to the field of agricultural technology, and in particular to a dynamic crop moisture estimation method that assimilates unmanned aerial vehicle observations and crop models. Background Art

[0002] According to the United Nations, the global population is projected to grow to 10.3 billion by 2080, significantly increasing food demand and exacerbating water shortages in irrigated agriculture. Efficient water use is a key factor influencing crop growth, development, and yield formation. Therefore, accurately identifying crop water stress and implementing precise irrigation management are crucial for ensuring normal crop growth and maintaining sustainable agricultural development. Drone remote sensing, with its high temporal and spatial resolution, offers a rapid, non-destructive, cost-effective, and highly promising solution for monitoring CWS. CWS, the primary abiotic stress factor, stems from insufficient water supply to crops due to soil moisture deficit, directly disrupting crop physiological and biochemical processes and morphological development. In other words, remote sensing technology primarily identifies stress by monitoring crop responses to water stress, rather than directly detecting the stress itself.

[0003] Currently, UAV-based CWS monitoring methods can be divided into two main categories: assessment methods based on water stress indices and direct monitoring methods based on phenotypic and physiological characteristics. The former primarily assesses the degree of stress by calculating a crop water stress index. These include a theoretical crop water stress index based on the principle of energy balance, an empirical crop water stress index based on the empirical relationship between canopy-air temperature difference and saturated water vapor pressure difference, and a statistical crop water stress index calculated based on temperature histogram quantiles. The latter directly explores the relationship between underlying data and CWS using traditional linear regression or machine learning algorithms. Key crop water indicators such as stomatal conductance, chlorophyll fluorescence, stem / leaf water potential, leaf water content, and plant water content are used to quantitatively assess the degree of water stress by comparing them with crops under well-watered conditions. It should be noted that CWS can be divided into two types: instantaneous and long-term. Instantaneous stress is mainly manifested as short-term fluctuations in leaf water content, while long-term stress can lead to significant changes in crop physical parameters. However, the discrete observation characteristics of UAVs make it difficult to achieve continuous monitoring. Existing studies usually use single or a few flight data to represent the water conditions throughout the day or even the entire growing period, which may lead to insufficient representation on the time scale and misdiagnosis of CWS. Therefore, the advantages of UAVs should be combined to develop a more stable CWS dynamic monitoring model.

[0004] Crop models are important tools for optimizing field management. They can simulate the physiological processes of crop growth under complex environmental constraints, thereby providing more complete information on crop water dynamics. However, in regional-scale modeling, uncertainties in initial conditions, crop parameters, and soil properties limit the spatial applicability of the models. Data assimilation technology integrates remote sensing observations with crop models and uses remote sensing inversion parameters to calibrate the models, significantly improving not only the estimation accuracy of crop state variables but also data utilization. Its core is to adjust model parameters or state variables through optimization algorithms to minimize the difference between observed and simulated data and achieve optimal fusion of the two. In recent years, researchers have made full use of the technical advantages of UAV remote sensing data and DA to make significant progress in improving the estimation accuracy of indicators such as leaf area index and yield. However, relatively insufficient attention has been paid to the dynamic monitoring of CWS.

[0005] It's important to note that while existing, established crop models can effectively assess yields and predict irrigation needs, and through assimilation with UAV data, dynamically simulate crop growth and development from various perspectives, offering significant mechanistic advantages and simulation accuracy, these models require complex input data, including weather, soil, crop, and field management data, making this DA approach less applicable. In contrast, simpler crop models are more suitable for DA than more complex ones. They require fewer parameters to describe crop growth, require relatively low input data, and offer fast computational speeds. Combined with the portability of UAVs, this significantly enhances the practicality and potential for widespread adoption of this approach.

[0006] The Simple Algorithm for Yield (SAFYE) model, a light-driven model that focuses on the relationship between light energy utilization and crop dry matter accumulation, has been widely used in biomass and yield estimation studies for crops such as wheat and corn. Based on the FAO-56 methodology, the SAFYE model, which includes a water balance module, has been further developed and applied in areas such as crop evapotranspiration simulation and irrigation scheduling. However, the SAFYE model's potential for simulating CWS has not been fully explored, especially in studies integrating UAV observations for dynamic monitoring and diagnosis of CWS throughout the growing season.

[0007] In view of the above technical defects, a solution is now proposed. Summary of the Invention

[0008] The purpose of this paper is to simulate the dynamic process of maize growth and ETc act under different water conditions by assimilating crop phenotypic parameters inverted from UAV multispectral data into the SAFYE model, and to evaluate the application potential of this technical framework for diagnosing CWS throughout the entire growth period.

[0009] To achieve the above objectives, the present invention adopts the following technical solution: a dynamic crop moisture estimation method that assimilates UAV observations and crop models, comprising the following steps:

[0010] Step 1: Set up several corn planting experimental plots and manage them in groups. Set up ground monitoring points at the four corners and the center of each corn planting experimental plot. Use the LAI-2200C canopy analyzer to measure once at each ground monitoring point. Take the average value as the leaf area index value to be processed. Deploy an automated continuous monitoring system at each ground monitoring point to continuously monitor the soil volumetric moisture content.

[0011] The management process of the corn planting experimental area in step 1 specifically includes:

[0012] Corn is planted in a wide-narrow row pattern, with a wide row spacing of 0.8m, a narrow row spacing of 0.4m, and a plant spacing of 0.2m. The planting density is approximately 67,000 plants / hm2. 2 ,The irrigation method is drip irrigation, and the spacing between drip irrigation tapes is 1.2m;

[0013] Based on the irrigation system of the corn planting experimental area, four irrigation treatments were set up: severe water shortage W1, moderate water shortage W2, mild water shortage W3, and local experience irrigation quota W4. W1, W2, and W3 were set at 40%, 60%, and 80% of the irrigation quota of W4, respectively. Two replicate corn planting experimental areas were set up for each treatment. Each corn planting experimental area was equipped with a turbine flowmeter for independent water control. A 1.5-meter separation zone was set between the corn planting experimental areas. In early November of the year before planting, both corn planting experimental areas were irrigated with high-quota Yellow River water. The irrigation amount was based on the local farmers' experience. After irrigation, the fields were frozen and sown after melting in April of the following year.

[0014] Step 2: Use a Phantom 4 Multispectral drone to collect image data. DJI Terra 3.6.8 software reconstructs and stitches the image data and adds black and white ground plate photos for radiometric correction. This creates a multispectral orthophoto map of the corn planting experimental area, which is then georeferenced by adding GCPs.

[0015] Step 3: Based on the core theoretical framework of the SAFY model, the SAFYE model is expanded and coupled with the FAO-56 water balance module to achieve joint simulation of crop growth and LAI development, biomass formation, evapotranspiration and soil moisture dynamics. The biomass accumulation module of the SAFYE model is based on Monteith's light energy utilization theory, and its daily biomass increment ΔAGB is determined by the global radiation R g , light energy utilization efficiency LUE, climate efficiency ε c and light interception efficiency ε IThe temperature stress coefficient FT is jointly determined and introduced to explain the effect of temperature on crop growth and development;

[0016] Step 4: Obtain observable biophysical parameters to calibrate the SAFYE model. These parameters include leaf area index, yield, and biomass distribution among various organs. Based on field observations of these biophysical parameters, the SAFYE model parameters are optimized at the regional scale. Four key crop morphological parameters that control crop growth and development are optimized, ultimately outputting a global optimal solution.

[0017] Step 5: Obtain the leaf area index and aboveground biomass of corn based on the multispectral orthophotos of the corn planting experimental area, and invert key crop morphological parameters using a machine learning algorithm. The specific process is as follows:

[0018] Screening preset multispectral indices as key estimation indicators, wherein the multispectral indices are calculated based on reflectance data of five characteristic bands: 450, 560, 650, 730, and 840 nm;

[0019] The machine learning regression function in the Scikit-Learn package was used to build and invert the model in the Python 3.10 environment;

[0020] Step 6: Obtain key crop morphological parameters as external assimilation data for the SAFYE model DA. Use the Monte Carlo variant of the Kalman filter (EnKF) to calculate the prediction error covariance matrix using the predicted value set. Assimilate the key crop morphological parameters into the crop growth SAFYE model to optimize the model simulation process and obtain the optimized SAFYE model.

[0021] Step 7: Use the optimized SAFYE model to simulate ET throughout the corn growth period c act Dynamic changes are used to build a CWS diagnostic model, and then the crop water stress index is calculated to obtain the dynamic changes of crop water stress conditions, and a customized irrigation plan is generated based on the dynamic simulation of the crop water stress index difference.

[0022] Furthermore, the Phantom 4 Multispectral UAV in step 2 is equipped with one RGB sensor and five monochrome multispectral sensors, and the flight parameters are set as follows: an altitude of 20 m, a ground sampling distance of 1.06 cm / pixel, and an 80% overlap rate between heading and sideways.

[0023] Furthermore, the SAFYE model in step 3 is expressed as follows:

[0024] ΔAGB=APAR×LUE×F T (T a )×Ks ;

[0025] The light and effective radiation APAR are calculated by the Beer-Lanmbert law, and the light and effective radiation APAR are related to the leaf area index LAI and the global radiation R g It is an exponential relationship, as shown below:

[0026] APAR=1-e -k×LAI ×R g ×ε c ;

[0027]

[0028] The base of the natural logarithm in the e exponential function, T a is the temperature, T min is the minimum temperature suitable for corn growth, T max is the highest temperature suitable for corn growth, T opt is the optimal temperature for corn growth, β0 is the preset polynomial degree;

[0029] The aboveground biomass AGB allocation process is based on the following rules: the daily biomass increment ΔAGB is allocated to the leaf biomass through the leaf allocation function Pl, and then converted into the leaf area increment ΔLAI+ through the specific leaf area SLA. When the accumulated temperature AT exceeds the senescence threshold STT, the leaf senescence ΔLAI- is triggered, which is expressed by the following formula:

[0030] ΔLAI + =ΔAGB×Pl×SLA,Pl>0;

[0031]

[0032] Where Pl represents the leaf allocation function, Pl a represents the leaf allocation function with a value range of (0.05-0.5), Pl b represents the leaf allocation function with a value range of (0.001-0.01);

[0033] S TT Indicates the aging temperature threshold, which is used to determine the time when LAI begins to decrease. The STT value range is 500℃-1200℃;

[0034] R s Indicates the aging rate, with a value range of 2000℃-15000℃;

[0035] H I represents the harvest index, which is the ratio of final grain yield to final DAM, and is set to 0.5;

[0036] The process of calculating the water balance part based on the water balance equation is as follows:

[0037] Evapotranspiration is determined by the double crop coefficient method, and the water stress coefficient Ks is added to reflect the water stress of the crop:

[0038] ET=(K e +K cb ×K s )×ET0;

[0039] The soil evaporation coefficient Ke is determined based on the crop coverage CC, and is further determined based on the evaporation reduction coefficient β and the relative soil water depletion RH of the soil layer;

[0040] The crop coverage CC is calculated according to the following formula:

[0041] K e =(1-CC)×[1-(1-RH1) β ], where β is the evaporation reduction coefficient, which is set to 0.94;

[0042] CC=76.78×(1-e -0.8105×LAI ) 0.9567 , where e is the base of the natural logarithm preset in the exponential function;

[0043] RH x is the relative humidity of the soil layer, CAW x is the current available water content in the soil layer;

[0044] SC x =(θ fc -θ wp )×SLT x , SC x is the current soil water storage capacity, SLT x is the soil layer thickness, θ fc and is the field capacity, θ wp is the wilting moisture content;

[0045] Through the basic crop coefficient K cb Indicates the potential transpiration rate of crops, basic crop coefficient K cb Depending on the development stage of the crop, the maximum transpiration and leaf area index of the crop are used to Kcb Simplify, where K cb,max is the maximum transpiration coefficient, K trp is the transpiration index, expressed as follows:

[0046] Where e is the base of the natural logarithm preset in the exponential function;

[0047] The crop water stress coefficient Ks is expressed in the form of a simple linear model as follows:

[0048]

[0049] θ t =(1-p)(θ fc -θ wp )+θ wp ;

[0050] Where θ is the soil water content in the root zone, θ t is the threshold water content, p is the average fraction of available soil water that can be removed from the root zone before water stress;

[0051] The daily change of soil water in the area is calculated based on the soil water balance principle and is expressed as follows:

[0052] ΔCAW = P + I + GF-ET-SR-DP;

[0053] where P, I, GF, ET, SR, and DP represent precipitation, irrigation, groundwater recharge, evapotranspiration, runoff, and deep seepage, respectively, and ΔCAW is the daily change in soil water in the root zone;

[0054] The corn planting experimental area in the present invention is an arid plain, and groundwater, runoff and leakage are not considered. The simplified formula is:

[0055] ΔCAW=P+I-ET

[0056] For example:

[0057] The soil water content on the first day is: CAW1 = θ0 × 0.2 × 1000 + ΔCAW1;

[0058] Note: CAW(mm)=θ×d×1000;

[0059] θ is the volumetric water content of soil, θ0 is the initial soil water content, and d is the thickness of the soil layer, which is taken as 0.2 m;

[0060] The soil water content on the second day is: CAW2 = CAW1 + ΔCAW2;

[0061] ······;

[0062] Soil water content on day x is: CAW x =CAW x-1 +ΔCAW x ;

[0063] The initial soil moisture content θ0 is set to 0.2;

[0064] According to the above formula, corn evapotranspiration is calculated through continuous monitoring of soil moisture.

[0065] Furthermore, the specific process of optimizing the regional-scale SAFYE model parameters is as follows:

[0066] refer to Figure 3 , obtain the minimum measured leaf area index and aboveground biomass time series data, and obtain the root mean square error between the leaf area index LAIsim and aboveground biomass AGB simulated according to the SAFYE model to optimize specific parameters. The parameter calibration cost function taking the leaf area index as an example is as follows:

[0067]

[0068] Where i = 1, 2, 3, n, n is LAI mea The amount of data, LAI mea is the leaf area index observation value, LAI sim It is the leaf area index simulation value obtained according to the SAFYE model;

[0069] The leaf area index was determined between 08:00 and 10:00 on the day of image data collection. The LAI-2200C canopy analyzer was used to measure once at the four corners and the center of each corn planting experiment. A total of five measurement points were obtained, and the average value was taken as the LAI value of the treatment.

[0070] The aboveground biomass (AGB) was obtained by destructive sampling. The aboveground parts of corn plants were removed, crushed using a grinder, dried at 105°C for 30 minutes to inhibit metabolic activity, and further dried at 80°C to constant weight. The aboveground biomass was then estimated based on the planting density.

[0071] The PSO algorithm is used to achieve parameter optimization by simulating group intelligent behavior. A population of m particles is initialized in the D-dimensional search space. Each particle dynamically adjusts its motion state according to its individual historical optimal position and the group's optimal position to form a collaborative evolution mechanism. During the iteration process, the particle continuously updates its position until the fitness function converges or reaches the preset maximum number of iterations. The four key crop morphological parameters that control crop growth and development, namely Pl a Pl b 、S TT and R s Optimize and finally output the global optimal solution.

[0072] Furthermore, the specific model parameters in step five are set as follows: the number of decision trees ntree is 500, the minimum leaf tree mtry is 1, the data in the experimental process are obtained as training samples, and the training samples are divided into a training set and a test set in a ratio of 7:3.

[0073] Furthermore, the expression of the Monte Carlo variant EnKF of the Kalman filter in step 6 is as follows:

[0074] A represents the observation data analysis matrix of set member i, i represents the observation data simulation matrix of set member i, P t represents the simulated variance of the crop parameter, R e represents the variance of the observed data, D i Used to update collection members.

[0075] Furthermore, the specific process of calculating the crop water stress index in step seven is as follows:

[0076] The crop water stress index was calculated according to the following formula:

[0077] The crop water stress index is dimensionless and ranges from 0 to 1, where 0 and 1 represent no water stress and the most severe water stress, respectively. c act is the actual evapotranspiration, ET p is the potential evapotranspiration.

[0078] Furthermore, the optimized SAFYE model obtained in step 7 is used to evaluate the performance of the SAFYE model under the model assimilation framework using ground data from field observations. The evaluation process is as follows:

[0079]

[0080] Among them, O i is the observed value, S i is the predicted value, is the mean of the observed values, n is the number of samples, R is the coefficient of determination, RMSE is the root mean square error, NRMSE is the normalized root mean square error, and d is the consistency index;

[0081] The evaluation process uses the consistency index d, root mean square error, and normalized root mean square error as evaluation indices. The smaller the values ​​of the root mean square error and the normalized root mean square error, the closer the value of d is to 1, indicating that the SAFYE model has higher performance accuracy and better assimilation effect in simulating crop growth and development.

[0082] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0083] This dynamic crop moisture estimation method assimilates UAV observations and crop models. It uses high-temporal and spatial resolution UAV remote sensing data to invert key crop morphological parameters through a machine learning algorithm. The inverted key crop morphological parameters are assimilated into the crop growth model to optimize the model simulation process, and then a CWS diagnostic model is constructed to generate customized irrigation plans. The EnKF method is used to assimilate the LAI and AGB data acquired by the UAV. The effectiveness of the assimilation framework in simulating the entire crop growth period under different water stress conditions is evaluated. Data assimilation significantly reduces the simulation errors of LAI and AGB, verifying the feasibility of the assimilation framework. It can achieve complementary advantages of the methods. On the one hand, it provides real values ​​to correct model bias and reduce the uncertainty of initial conditions and parameters of modeling. On the other hand, it enhances the model's ability to analyze the crop growth and development mechanisms. Combined with the simplicity of the SAFYE model, the DA framework has the potential for regional promotion. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 Shows the overall process schematic diagram of the present invention;

[0085] Figure 2 The input and output parameters and model setting diagram of the SAFY model of the present invention are shown;

[0086] Figure 3 A schematic diagram showing the results of dynamically estimating the moisture status of crops during the growing season using CWSI in the present invention is shown. DETAILED DESCRIPTION

[0087] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. 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 making creative efforts are within the scope of protection of the present invention.

[0088] Example:

[0089] like Figure 1-3 As shown in Figure 2, the dynamic crop water status estimation method that assimilates UAV observations and crop models includes the following steps:

[0090] Step 1: Set up several corn planting experimental plots and manage them in groups. Set up ground monitoring points at the four corners and the center of each corn planting experimental plot. Use the LAI-2200C canopy analyzer to measure once at each ground monitoring point. Take the average value as the leaf area index value to be processed. Deploy an automated continuous monitoring system at each ground monitoring point to continuously monitor the soil volumetric moisture content.

[0091] The management process of the corn planting experimental area in step 1 specifically includes:

[0092] Corn is planted in a wide-narrow row pattern, with a wide row spacing of 0.8m, a narrow row spacing of 0.4m, and a plant spacing of 0.2m. The planting density is approximately 67,000 plants / hm2. 2 ,The irrigation method is drip irrigation, and the spacing between drip irrigation tapes is 1.2m;

[0093] Based on the irrigation system of the corn planting experimental area, four irrigation treatments were set up: severe water shortage W1, moderate water shortage W2, mild water shortage W3, and local experience irrigation quota W4. W1, W2, and W3 were set at 40%, 60%, and 80% of the irrigation quota of W4, respectively. Two replicate corn planting experimental areas were set up for each treatment. Each corn planting experimental area was equipped with a turbine flowmeter for independent water control. A 1.5-meter separation zone was set between the corn planting experimental areas. In early November of the year before planting, both corn planting experimental areas were irrigated with high-quota Yellow River water. The irrigation amount was based on the local farmers' experience. After irrigation, the fields were frozen and sown after melting in April of the following year.

[0094] Step 2: Use a Phantom 4 Multispectral drone to collect image data. DJI Terra 3.6.8 software reconstructs and stitches the image data and adds black and white ground plate photos for radiometric correction. This creates a multispectral orthophoto map of the corn planting experimental area, which is then georeferenced by adding GCPs.

[0095] The Phantom 4 Multispectral UAV in step 2 is equipped with one RGB sensor and five monochrome multispectral sensors. The flight parameters are set as follows: an altitude of 20 m, a ground sampling distance of 1.06 cm / pixel, and an 80% overlap rate between heading and sideways.

[0096] Step 3: Based on the core theoretical framework of the SAFY model, the SAFYE model is expanded and coupled with the FAO-56 water balance module to achieve joint simulation of crop growth and LAI development, biomass formation, evapotranspiration and soil moisture dynamics. The biomass accumulation module of the SAFYE model is based on Monteith's light energy utilization theory, and its daily biomass increment ΔAGB is determined by the global radiation R g , light energy utilization efficiency LUE, climate efficiency ε c and light interception efficiency ε I Jointly determine and introduce the temperature stress coefficient F T To explain the impact of temperature on crop growth and development, the input and output parameters of the SAFYE model and the model settings are referenced Figure 2 ;

[0097] The SAFYE model in step 3 is expressed as follows:

[0098] ΔAGB=APAR×LUE×F T (T a )×K s ;

[0099] The light and effective radiation APAR are calculated by the Beer-Lanmbert law, and the light and effective radiation APAR are related to the leaf area index LAI and the global radiation R g It is an exponential relationship, as shown below:

[0100] APAR=1-e -k×LAI ×R g ×ε c ;

[0101]

[0102] The base of the natural logarithm in the e exponential function, T a is the temperature, T min is the minimum temperature suitable for corn growth, T max is the highest temperature suitable for corn growth, T opt is the optimal temperature for corn growth, β0 is the preset polynomial degree;

[0103] The aboveground biomass AGB allocation process is based on the following rules: the daily biomass increment ΔAGB is allocated to the leaf biomass through the leaf allocation function Pl, and then converted into the leaf area increment ΔLAI+ through the specific leaf area SLA. When the accumulated temperature AT exceeds the senescence threshold S TT Leaf senescence is triggered when ΔLAI- is expressed by the following formula:

[0104] ΔLAI + =ΔAGB×Pl×SLA,Pl>0;

[0105]

[0106] Where Pl represents the leaf allocation function, Pl a represents the leaf allocation function with a value range of (0.05-0.5), Pl b represents the leaf allocation function with a value range of (0.001-0.01);

[0107] S TT Represents the aging temperature threshold, which is used to determine the moment when LAI begins to decrease. TT The value range is 500℃-1200℃;

[0108] R s Indicates the aging rate, with a value range of 2000℃-15000℃;

[0109] HI represents the harvest index, which is the ratio of final grain yield to final DAM, and is set to 0.5;

[0110] The process of calculating the water balance part based on the water balance equation is as follows:

[0111] Evapotranspiration is determined by the double crop coefficient method, and the water stress coefficient Ks is added to reflect the water stress of the crop:

[0112] ET=(K e +K cb ×K s )×ET0;

[0113] The soil evaporation coefficient Ke is determined based on the crop coverage CC, and is further determined based on the evaporation reduction coefficient β and the relative soil water depletion RH of the soil layer;

[0114] The crop coverage CC is calculated according to the following formula:

[0115] K e =(1-CC)×[1-(1-RH1) β ], where β is the evaporation reduction coefficient, which is set to 0.94;

[0116] CC=76.78×(1-e -0.8105×LAI ) 0.9567 , where e is the base of the natural logarithm preset in the exponential function;

[0117] RH x is the relative humidity of the soil layer, CAW x is the current available water content in the soil layer;

[0118] SC x =(θ fc -θ wp )×SLT x , SC x is the current soil water storage capacity, SLT x is the soil layer thickness, θ fc and is the field capacity, θ wp is the wilting moisture content;

[0119] Through the basic crop coefficient K cb Indicates the potential transpiration rate of crops, basic crop coefficient K cb Depending on the development stage of the crop, the maximum transpiration and leaf area index of the crop are used to determine K cb Simplify, where K cb,max is the maximum transpiration coefficient, K trp is the transpiration index, expressed as follows:

[0120] Where e is the base of the natural logarithm preset in the exponential function;

[0121] The crop water stress coefficient Ks is expressed in the form of a simple linear model as follows:

[0122]

[0123] θ t =(1-p)(θ fc -θ wp )+θ wp ;

[0124] Where θ is the soil water content in the root zone, θ t is the threshold water content, p is the average fraction of available soil water that can be removed from the root zone before water stress;

[0125] The daily change of soil water in the area is calculated based on the soil water balance principle and is expressed as follows:

[0126] ΔCAW = P + I + GF-ET-SR-DP;

[0127] where P, I, GF, ET, SR, and DP represent precipitation, irrigation, groundwater recharge, evapotranspiration, runoff, and deep seepage, respectively, and ΔCAW is the daily change in soil water in the root zone;

[0128] The corn planting experimental area in the present invention is an arid plain, and groundwater, runoff and leakage are not considered. The simplified formula is:

[0129] ΔCAW=P+I-ET

[0130] For example:

[0131] The soil water content on the first day is: CAW1 = θ0 × 0.2 × 1000 + ΔCAW1;

[0132] Note: CAW(mm)=θ×d×1000;

[0133] θ is the volumetric water content of soil, θ0 is the initial soil water content, and d is the thickness of the soil layer, which is taken as 0.2 m;

[0134] The soil water content on the second day is: CAW2 = CAW1 + ΔCAW2;

[0135] ······;

[0136] Soil water content on day x is: CAW x =CAW x-1 +ΔCAW x ;

[0137] The initial soil moisture content θ0 is set to 0.2;

[0138] According to the above formula, corn evapotranspiration is calculated through continuous monitoring of soil moisture.

[0139] Step 4: Obtain observable biophysical parameters to calibrate the SAFYE model. These parameters include leaf area index, yield, and biomass distribution among various organs. Based on field observations of these biophysical parameters, the SAFYE model parameters are optimized at the regional scale. Four key crop morphological parameters that control crop growth and development are optimized, ultimately outputting a global optimal solution.

[0140] The specific process of optimizing the regional-scale SAFYE model parameters is as follows:

[0141] Obtain the time series data of leaf area index and aboveground biomass obtained by minimizing the actual measurement, and obtain the leaf area index LAI simulated according to the SAFYE model sim The root mean square error between the AGB and the aboveground biomass is used to optimize specific parameters. The parameter calibration cost function taking the leaf area index as an example is as follows:

[0142]

[0143] Where i = 1, 2, 3, n, n is LAI mea The amount of data, LAI mea is the leaf area index observation value, LAI sim It is the leaf area index simulation value obtained according to the SAFYE model;

[0144] The leaf area index was determined between 08:00 and 10:00 on the day of image data collection. The LAI-2200C canopy analyzer was used to measure once at the four corners and the center of each corn planting experiment. A total of five measurement points were obtained, and the average value was taken as the LAI value of the treatment.

[0145] The aboveground biomass (AGB) was obtained by destructive sampling. The aboveground parts of corn plants were removed, crushed using a grinder, dried at 105°C for 30 minutes to inhibit metabolic activity, and further dried at 80°C to constant weight. The aboveground biomass was then estimated based on the planting density.

[0146] The PSO algorithm is used to achieve parameter optimization by simulating group intelligent behavior. A population of m particles is initialized in the D-dimensional search space. Each particle dynamically adjusts its motion state according to its individual historical optimal position and the group's optimal position to form a collaborative evolution mechanism. During the iteration process, the particle continuously updates its position until the fitness function converges or reaches the preset maximum number of iterations. The four key crop morphological parameters that control crop growth and development, namely Pla Pl b 、S TT and R s Optimize and finally output the global optimal solution;

[0147] Step 5: Obtain the leaf area index and aboveground biomass of corn based on the multispectral orthophoto map of the corn planting experimental area, and invert key crop morphological parameters through machine learning algorithms. The specific process is as follows, refer to Figure 3 :

[0148] Screening preset multispectral indices as key estimation indicators, wherein the multispectral indices are calculated based on reflectance data of five characteristic bands: 450, 560, 650, 730, and 840 nm;

[0149] The machine learning regression function in the Scikit-Learn package was used to build and invert the model in the Python 3.10 environment;

[0150] The specific model parameter settings are as follows: the number of decision trees ntree is 500, the minimum leaf tree mtry is 1, the data obtained during the experiment are used as training samples, and the training samples are divided into training set and test set in a ratio of 7:3;

[0151] Step 6: Obtain key crop morphological parameters as external assimilation data for the SAFYE model DA. Use the Monte Carlo variant of the Kalman filter (EnKF) to calculate the prediction error covariance matrix using the predicted value set. Assimilate the key crop morphological parameters into the crop growth SAFYE model to optimize the model simulation process and obtain the optimized SAFYE model.

[0152] The expression of the Monte Carlo variant EnKF of the Kalman filter in step 6 is as follows:

[0153] A represents the observation data analysis matrix of set member i, i represents the observation data simulation matrix of set member i, P t represents the simulated variance of the crop parameter, R e represents the variance of the observed data, D i Used to update collection members;

[0154] Step 7: Use the optimized SAFYE model to simulate ET throughout the corn growth period c act Dynamic changes are used to build a CWS diagnostic model, and then the crop water stress index is calculated to obtain the dynamic changes of crop water stress conditions. Based on the dynamic simulation of the difference in crop water stress index, a customized irrigation plan is generated. Figure 3 .

[0155] The specific process of calculating the crop water stress index in step seven is as follows:

[0156] The crop water stress index was calculated according to the following formula:

[0157] The crop water stress index is dimensionless and ranges from 0 to 1, where 0 and 1 represent no water stress and the most severe water stress, respectively. c act is the actual evapotranspiration, ET p is the potential evapotranspiration;

[0158] The optimized SAFYE model obtained in step 7 is used to evaluate the performance of the SAFYE model under the model assimilation framework using ground data observed in the field. The evaluation process is as follows:

[0159]

[0160] Among them, O i is the observed value, S i is the predicted value, is the mean of the observed values, n is the number of samples, R is the coefficient of determination, RMSE is the root mean square error, NRMSE is the normalized root mean square error, and d is the consistency index;

[0161] The evaluation process uses the consistency index d, root mean square error, and normalized root mean square error as evaluation indices. The smaller the values ​​of the root mean square error and the normalized root mean square error, the closer the value of d is to 1, indicating that the SAFYE model has higher performance accuracy and better assimilation effect in simulating crop growth and development.

[0162] The specific process of generating a customized irrigation plan is as follows:

[0163] ET simulated by integrating the optimized SAFYE model c act The dynamic calculation of ET0 with the Penman-Monteith method achieved the dynamic quantification of the crop water stress index over multiple corn growing seasons, indicating that the dynamic changes in the crop water stress index exhibited typical irrigation response characteristics.

[0164] It drops significantly after each irrigation or rainfall event, then gradually recovers until the next water replenishment. From the perspective of the growth period, the crop water stress index maintains a high level in the early growth period, fluctuates and decreases in the middle growth period due to irrigation and the arrival of the rainy season, and gradually recovers and stabilizes in the late growth period.

[0165] The differences in crop water stress index values ​​between different treatments revealed that as the amount of irrigation decreased, the decrease in the crop water stress index also decreased. For example, compared with the severely water-deficient W1 treatment, the well-irrigated W4 treatment had the largest decrease and more days with a crop water stress index of 0.

[0166] Comparative analysis among different irrigation treatments showed that the response of crop water stress index to irrigation amount had a dose effect: the greater the irrigation amount, the more significant the decrease in crop water stress index. Taking the fully irrigated W4 treatment (set according to local experience) as the benchmark, the crop water stress index had more days with a 0 value and the largest decrease in the crop water stress index.

[0167] Based on this, it is assumed that corn in other treatments with a crop water stress index greater than that of the W4 treatment will suffer from CWS. That is, the difference in crop water stress index between each treatment and the W4 treatment is used as a quantitative indicator of water stress. The higher the difference in crop water stress index, the more severe the CWS suffered by the corn. Analysis shows that the more severe the water deficit (W1>W2>W3), the earlier the difference in crop water stress index appears and the higher its peak value, which is consistent with the actual CWS situation of the crop. Therefore, the dynamic simulation based on the difference in crop water stress index can provide a basis for crop irrigation adjustment.

[0168] This paper uses UAV remote sensing data with high temporal and spatial resolution to invert key crop morphological parameters through a machine learning algorithm. The inverted key crop morphological parameters are assimilated into the crop growth model to optimize the model simulation process, and then a CWS diagnostic model is constructed to generate customized irrigation plans. The LAI and AGB data acquired by UAV are assimilated using the EnKF method. The effectiveness of the assimilation framework in simulating the complete growth period of crops under different water stress conditions is evaluated. Data assimilation significantly reduces the simulation errors of LAI and AGB, verifying the feasibility of the assimilation framework. It can achieve complementary advantages of the methods. On the one hand, it provides real values ​​to correct model deviations and reduce the uncertainty of initial conditions and parameters of modeling. On the other hand, it enhances the model's ability to analyze the crop growth and development mechanisms. Combined with the simplicity of the SAFYE model, the DA framework has the potential for regional promotion.

[0169] The size of the interval and threshold is set to facilitate comparison. The size of the threshold depends on the amount of sample data and the number of bases set by technical personnel in this field for each set of sample data; as long as it does not affect the proportional relationship between the parameter and the quantized value.

[0170] The above formulas are all dimensionless and numerical calculations. The formulas are obtained by collecting a large amount of data and performing software simulation to obtain the most recent real situation. The preset parameters in the formulas are set by those skilled in the art according to actual conditions.

[0171] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with the technical field, within the technical scope disclosed by the present invention, who makes equivalent replacements or changes based on the technical solution and inventive concept of the present invention, should be covered by the scope of protection of the present invention.

Claims

1. A dynamic crop moisture estimation method that assimilates UAV observations and crop models, characterized by: The following steps are involved: Step 1: Set up several corn planting experimental plots and manage them in groups. Set up ground monitoring points at the four corners and the center of each corn planting experimental plot. Use the LAI-2200C canopy analyzer to measure once at each ground monitoring point. Take the average value as the leaf area index value to be processed. And deploy an automated continuous monitoring system at the ground monitoring points. Step 2: Use a Phantom 4 Multispectral drone to collect image data. DJI Terra 3.6.8 software reconstructs and stitches the image data and adds ground black and white plate photos for radiometric correction. This produces a multispectral orthophoto map of the corn planting experimental area, which is then georeferenced by adding GCPs. Step 3: Based on the core theoretical framework of the SAFY model, the SAFYE model is expanded and coupled with the FAO-56 water balance module. The biomass accumulation module of the SAFYE model is based on the Monteith light energy utilization theory, and its daily biomass increment ΔAGB is determined by the global radiation R g , light energy utilization efficiency LUE, climate efficiency ε c and light interception efficiency ε I Determine together and introduce the temperature stress coefficient F T To explain the effects of temperature on crop growth and development; Step 4: Obtain observable biophysical parameters to calibrate the SAFYE model. Based on field observations of biophysical parameters, optimize the SAFYE model parameters at the regional scale. Optimize key crop morphological parameters that control crop growth and development, and ultimately output a global optimal solution. Step 5: Obtain the leaf area index and aboveground biomass of corn based on the multispectral orthophoto map of the corn planting experimental area, and invert key crop morphological parameters using a machine learning algorithm; Step 6: Obtain key crop morphological parameters as external assimilation data for the SAFYE model DA. Use the Monte Carlo variant of the Kalman filter (EnKF) to calculate the prediction error covariance matrix using the predicted value set. Assimilate the key crop morphological parameters into the crop growth SAFYE model to optimize the model simulation process and obtain the optimized SAFYE model. Step 7: Use the optimized SAFYE model to simulate ET throughout the corn growth period cact Dynamic changes are used to build a CWS diagnostic model, and then the crop water stress index is calculated to obtain the dynamic changes of crop water stress conditions, and a customized irrigation plan is generated based on the dynamic simulation of the crop water stress index difference.

2. The method for dynamic crop moisture estimation based on UAV observations and crop models according to claim 1, characterized in that: The Phantom 4 Multispectral drone in step 2 is equipped with one RGB sensor and five monochrome multispectral sensors. The flight parameters are set as follows: an altitude of 20 meters, a ground sampling distance of 1.06 cm / pixel, and an 80% overlap rate between heading and sideways.

3. The method for dynamic crop moisture estimation based on UAV observations and crop models according to claim 1, characterized in that: The SAFYE model in step 3 is expressed as follows: ΔAGB=APAR×LUE×F T (T a )×K s ; The light and effective radiation APAR are calculated by the Beer-Lanmbert law, and the light and effective radiation APAR are related to the leaf area index LAI and the global radiation R g It is an exponential relationship, as shown below: APAR=1-e -k×LAI ×R g ×ε c ; The base of the natural logarithm in the e exponential function, T a is the temperature, T min is the minimum temperature suitable for corn growth, T max is the highest temperature suitable for corn growth, T opt is the optimal temperature for corn growth, β0 is the preset polynomial degree; The aboveground biomass AGB allocation process is based on the following rules: the daily biomass increment ΔAGB is allocated to the leaf biomass through the leaf allocation function Pl, and then converted into the leaf area increment ΔLAI+ through the specific leaf area SLA. When the accumulated temperature AT exceeds the senescence threshold STT, the leaf senescence ΔLAI- is triggered, which is expressed by the following formula: ΔLAI + =ΔAGB×Pl×SLA,Pl>0; Where Pl represents the leaf allocation function, Pl a represents the leaf allocation function with a value range of (0.05-0.5), Pl b represents the leaf allocation function with a value range of (0.001-0.01); S TT Indicates the aging temperature threshold, which is used to determine the time when LAI begins to decrease. The STT value range is 500℃-1200℃; R s Indicates the aging rate, with a value range of 2000℃-15000℃; H I It represents the harvest index, which is the ratio of final grain yield to final DAM, and is set to 0.5; The process of calculating the water balance part based on the water balance equation is as follows: Evapotranspiration is determined by the double crop coefficient method, and the water stress coefficient Ks is added to reflect the water stress of the crop: AND=(K e +K cb ×K s )×ET0; The soil evaporation coefficient Ke is determined based on the crop coverage CC, and is further determined based on the evaporation reduction coefficient β and the relative soil water depletion RH of the soil layer; The crop coverage CC is calculated according to the following formula: K e =(1-CC)×[1-(1-RH1) β ], where β is the evaporation reduction coefficient, which is set to 0.94; CC=76.78×(1-e -0.8105×LAI ) 0.9567 , where e is the base of the natural logarithm preset in the exponential function; RH x is the relative humidity of the soil layer, CAW x is the current available water content in the soil layer; SC x =(θ fc -θ wp )×SLT x , SC x is the current soil water storage capacity, SLT x is the soil layer thickness, θ fc and is the field capacity, θ wp is the wilting moisture content; Through the basic crop coefficient K cb Indicates the potential transpiration rate of crops, basic crop coefficient K cb Depending on the development stage of the crop, the maximum transpiration and leaf area index of the crop are used to Kcb Simplify, where K cb,max is the maximum transpiration coefficient, K trp is the transpiration index, expressed as follows: Where e is the base of the natural logarithm preset in the exponential function; The crop water stress coefficient Ks is expressed in the form of a simple linear model as follows: i t =(1-p)(θ) fc -θ wp )+θ wp ; Where θ is the soil water content in the root zone, θ t is the threshold water content, p is the average fraction of available soil water that can be removed from the root zone before water stress; The daily change of soil water in the area is calculated based on the soil water balance principle and is expressed as follows: ΔCAW = P + I + GF-ET-SR-DP; where P, I, GF, ET, SR, and DP represent precipitation, irrigation, groundwater recharge, evapotranspiration, runoff, and deep seepage, respectively, and ΔCAW is the daily change in soil water in the root zone.

4. The method for dynamic crop moisture estimation based on UAV observations and crop models according to claim 1, characterized in that: The specific process of optimizing the SAFYE model parameters at the regional scale is as follows: Referring to Figure 3, the cost function for parameter calibration using the leaf area index as an example is as follows: Where i = 1, 2, 3, n, n is LAI mea The amount of data, LAI mea is the leaf area index observation value, LAI sim It is the leaf area index simulation value obtained according to the SAFYE model; The PSO algorithm is used to achieve parameter optimization by simulating group intelligent behavior. A population of m particles is initialized in the D-dimensional search space. Each particle dynamically adjusts its motion state according to its individual historical optimal position and the group's optimal position to form a collaborative evolution mechanism. During the iteration process, the particle continuously updates its position until the fitness function converges or reaches the preset maximum number of iterations. The four key crop morphological parameters that control crop growth and development, namely Pl a Pl b 、S TT and R s Optimize and finally output the global optimal solution.

5. The method for dynamic crop moisture estimation based on UAV observations and crop models according to claim 1, characterized in that: The specific model parameters in step five are set as follows: the number of decision trees ntree is 500, the minimum leaf tree mtry is 1, the data in the experimental process are obtained as training samples, and the training samples are divided into a training set and a test set in a ratio of 7:

3.

6. The method for dynamic crop moisture estimation based on UAV observations and crop models according to claim 1, characterized in that: The expression of the Monte Carlo variant EnKF of the Kalman filter in step 6 is as follows: A represents the observation data analysis matrix of set member i, i represents the observation data simulation matrix of set member i, P t represents the simulated variance of the crop parameter, R e represents the variance of the observed data, D i Used to update collection members.

7. The method for dynamic crop moisture estimation based on UAV observations and crop models according to claim 1, wherein: The specific process of calculating the crop water stress index in step seven is as follows: The crop water stress index was calculated according to the following formula: The crop water stress index is dimensionless and ranges from 0 to 1, where 0 and 1 represent no water stress and the most severe water stress, respectively. cact is the actual evapotranspiration, ET p is the potential evapotranspiration.

8. The method for dynamic crop moisture estimation based on UAV observations and crop models according to claim 1, wherein: The optimized SAFYE model obtained in step 7 is used to evaluate the performance of the SAFYE model under the model assimilation framework using ground data observed in the field. The evaluation process is as follows: Among them, O i is the observed value, S i is the predicted value, is the mean of the observed values, n is the number of samples, R is the coefficient of determination, RMSE is the root mean square error, NRMSE is the normalized root mean square error, and d is the consistency index.

Images