Irrigation optimization method based on crop model

Abstract

The application discloses an irrigation optimization method based on a crop model, and comprises the following steps: acquiring input data such as meteorological data, soil data and crop data; running a crop growth model to simulate crop growth and water and salt dynamics under an irrigation scheme, wherein the model dynamically calculates root water absorption stress according to calculated crop growth stage indexes; and based on simulation results of the model and according to multiple optimization targets such as water saving, yield maintenance and salt control, an optimal irrigation scheme is obtained by solving the multiple optimization targets by using a multi-objective optimization algorithm. By introducing a stress calculation mechanism that changes with the growth stage, the application improves the simulation accuracy, and can efficiently solve the problem in combination with a proxy model, so that the obtained irrigation scheme is more reliable and practical in balancing the targets.

Description

Technical Field

[0001] This invention relates to the field of agricultural water and soil resource management and crop modeling technology, and in particular to an irrigation optimization method based on crop models. Background Technology

[0002] Against the backdrop of increasingly prominent water scarcity and soil salinization problems, there is a need to develop precision irrigation technologies to achieve synergistic optimization of water conservation, yield preservation, and salt control. Process-oriented crop growth models, as powerful numerical tools capable of quantitatively simulating the dynamic interaction between crop growth and soil water and salt, have become a core technical means for researching and formulating irrigation plans. By combining crop models with optimization algorithms, rapid simulation and optimization can be performed in various irrigation scenarios, providing scientific decision support for field irrigation management. Therefore, it is necessary to improve the simulation accuracy of crop models themselves and the coupling efficiency between them and optimization algorithms.

[0003] Currently, various crop growth models, such as the soil-water-atmosphere-plant system models SWAP, DSSAT, and AquaCrop, are widely used in irrigation optimization research. These models can comprehensively simulate complex processes such as soil moisture movement, solute transport, crop growth and development, and yield formation based on physical laws such as the Richards equation. In practical applications, these models are often used as the objective function of an unknown analytical expression and coupled with optimization algorithms such as genetic algorithms and particle swarm optimization. Within the model, the crop response to drought and salinity stress is usually calculated based on a fixed set of stress threshold parameters using empirical stress response functions (such as the Feddes model or the Maas-Hoffman model).

[0004] However, existing technologies have limitations in simulating the dynamic characteristics of crop stress responses, which in turn affects the accuracy and practicality of optimization schemes. Specifically, existing models generally use fixed stress response parameters to describe crop tolerance throughout the entire growth period. This static approach ignores the objective fact that crops exhibit differences in physiological functions and sensitivity to water and salt stress at different growth stages (such as seedling, flowering, and maturity), resulting in insufficient accuracy in simulating crop-water-salt interactions. Consequently, when optimization algorithms solve based on biased simulation results, the resulting optimal irrigation scheme may differ significantly from the actual optimal solution, making it difficult to effectively guide production practices. Furthermore, since crop models themselves are computationally time-consuming, directly coupling them with multi-objective optimization algorithms that require extensive function evaluations incurs enormous computational costs, limiting the full exploration of complex irrigation schemes in practice. Summary of the Invention

[0005] The purpose of this invention is to provide an irrigation optimization method based on crop models, in order to solve the aforementioned problems existing in the prior art.

[0006] Technical solution: An irrigation optimization method based on crop models, comprising:

[0007] Acquire input data for irrigation optimization, including meteorological, soil, and crop data for the target area;

[0008] Based on the input data, the crop growth model is run to simulate crop growth and soil water and salt dynamics under the irrigation scheme to be evaluated, and the corresponding simulation results are obtained. During the operation of the crop growth model, the root water uptake stress is dynamically calculated according to the crop growth stage.

[0009] Based on the simulation results and according to multiple preset optimization objectives, the optimal irrigation scheme is obtained by iteratively solving the irrigation scheme through a multi-objective optimization algorithm.

[0010] Output the optimal irrigation plan.

[0011] Beneficial Effects: Addressing the problem that existing models, employing fixed parameters, cannot accurately simulate the dynamic characteristics of crop stress responses, this invention proposes a complete dynamic stress calculation framework. This framework involves constructing continuously changing crop growth stage indices. These indices, as independent variables, drive key parameters in the stress response function (such as salt tolerance thresholds or osmotic regulation sensitivity coefficients) to change dynamically in real time throughout the crop's physiological development. By embedding such biologically more realistic dynamic response mechanisms into the model's physical calculations, the accuracy of the model's simulation of crop growth processes under complex water and salt stress conditions is improved, providing a more reliable simulation foundation for optimization solutions.

[0012] To address the issue of excessively high computational costs associated with directly coupling crop models with multi-objective optimization algorithms, this invention introduces an optimization strategy based on a surrogate model. This strategy approximates the high-cost original crop growth model using a surrogate model with extremely low computational cost (such as a Kriging model), enabling the optimization algorithm to complete massive iterative calculations within an acceptable timeframe. Furthermore, this invention proposes a sequential self-addition strategy to intelligently construct and iteratively optimize this surrogate model. This strategy can efficiently improve the prediction accuracy of the surrogate model in key regions with minimal calls to the original model, ensuring both efficiency and accuracy of the final optimization results.

[0013] This invention achieves efficient solutions to complex irrigation optimization problems by using dynamic stress and surrogate model optimization, while improving simulation accuracy. This results in an irrigation scheme that is both scientific and practically feasible. Attached Figure Description

[0014] Figure 1 This is a general framework diagram of an irrigation optimization method based on a crop model in an embodiment of this application.

[0015] Figure 2 This is a flowchart illustrating the steps involved in implementing the split-type stress frame in an embodiment of this application.

[0016] Figure 3 This is a flowchart illustrating the steps for calculating the salt stress reduction coefficient in an embodiment of this application.

[0017] Figure 4 This is a flowchart illustrating the steps involved in implementing a unified and effective water potential framework in the embodiments of this application. Detailed Implementation

[0018] Example 1: This example provides an overall framework for an irrigation optimization method based on a crop model, such as... Figure 1 As shown.

[0019] To clearly illustrate the present invention, the following embodiments will use the optimization of irrigation in saline-alkali land for winter wheat in a specific region as an example. However, those skilled in the art will understand that the method of the present invention is not limited to a specific crop type or region, but is also applicable to other crops such as corn and cotton, as well as different water and salt environments. In this embodiment, the method can be implemented based on an improved crop growth and water and salt transport model, for example, by improving the kernel level of the well-known SWAP model.

[0020] The overall technical architecture of the method of this invention can be summarized as a three-stage progressive process, specifically including: the first stage, improving the kernel of the crop growth model to enable it to respond to dynamic stress calculation mechanisms in response to crop growth stages; the second stage, constructing a computationally efficient surrogate model based on the improved model; and the third stage, using the surrogate model and a multi-objective optimization algorithm to solve irrigation schemes. Subsequent embodiments will elaborate on the specific technical details of each stage.

[0021] Step 101: Obtain input data for irrigation optimization, including meteorological data, soil data, and crop data for the target area.

[0022] Specifically, input data can further include the following types:

[0023] Meteorological data can specifically include time series data from standard meteorological stations, including daily maximum and minimum temperatures, rainfall, average wind speed, sunshine duration, and relative humidity.

[0024] Soil data can specifically be a stratified set of parameters describing the physical properties of a soil profile. For example, it may include the total depth of the soil profile, the thickness of each soil layer, bulk density, saturated water content, field capacity, wilting point water content, as well as van Genuchten model parameters and saturated hydraulic conductivity used to describe soil moisture dynamics.

[0025] Crop data can specifically be a set of parameters describing the physiological and developmental characteristics of a particular crop variety. For example, it may include the crop's sowing date, harvest date, potential root growth depth, crop coefficient, and developmental parameters calculated for subsequent dynamic stresses, such as the crop's basal temperature for development and the total effective accumulated temperature required to complete the entire growth cycle.

[0026] In addition, to achieve more accurate simulation and optimization, the input data may optionally include engineering parameters of the irrigation and drainage system, as well as historical observation data for parameter calibration.

[0027] Step 102: Based on the input data, run the crop growth model to simulate crop growth and soil water and salt dynamics under the irrigation scheme to be evaluated, and obtain the corresponding simulation results. During the operation of the crop growth model, the root water uptake stress is dynamically calculated according to the crop growth stage.

[0028] Specifically, the mechanism used in the crop growth model to dynamically calculate root water uptake stress includes a set of pre-configured model parameters. These parameters are obtained through an automated inversion calibration method, which includes: providing a pre-stored observation dataset containing field observation results; defining an objective function to calculate the error between the simulation results of the crop growth model and the observation dataset; and using an automated optimization algorithm to solve for the set of pre-configured model parameters with the objective function as the goal.

[0029] The crop growth model, acting as a simulation engine, receives input data and the irrigation scheme to be evaluated. The irrigation scheme can be represented as a time-series vector, with elements representing the amount of irrigation water applied at preset irrigation decision points. After the model runs, the output simulation results are one or more sets of performance indicators used to evaluate the merits of the irrigation scheme, such as the total crop yield at the end of the season, the total water consumption throughout the growing season, and the average salinity of the root zone soil at the end of the season.

[0030] Specifically, the crop growth model dynamically calculates root water uptake stress. Most existing crop models use a fixed set of stress response parameters when calculating the stress of water and salt on root water uptake, failing to accurately reflect the varying sensitivities of crops to drought and salt stress at different growth stages. For example, crops are extremely sensitive to water stress during flowering but have higher tolerance during maturity. This invention, by introducing a computational mechanism that dynamically changes with crop growth stages, can more realistically simulate stress effects, providing a more accurate physical basis for subsequent optimization solutions. Subsequent embodiments will detail two preferred technical paths for implementing this dynamic calculation.

[0031] Step 103: Based on the simulation results and according to multiple preset optimization objectives, the optimal irrigation scheme is obtained by iteratively solving the irrigation scheme through a multi-objective optimization algorithm.

[0032] In this embodiment, the simulation results are used as the evaluation basis for the objective function. In this embodiment, multiple optimization objectives aim to find the best balance between water conservation, yield preservation, and salt control. For example, they can be specifically defined as minimizing total irrigation water consumption, maximizing relative crop yield, and minimizing soil salinity in the root zone at the end of the season.

[0033] Since multiple optimization objectives are often conflicting, multi-objective optimization algorithms are used to solve them. The solution, the optimal irrigation scheme, is technically not a single solution, but rather a set of schemes known as the Pareto Optimal Set. Each scheme in this set represents a unique trade-off strategy, meaning that under this scheme, no single optimization objective can be further improved without sacrificing at least one other objective. Subsequent examples will elaborate on the specific strategies for using surrogate models to assist multi-objective optimization algorithms.

[0034] Step 104: Output the optimal irrigation plan.

[0035] The output can be a complete set of Pareto optimal solutions for decision-makers to choose from based on their actual preferences. In a preferred embodiment, this set of solutions can be further analyzed to output a single optimal compromise. The specific form of the optimal irrigation plan is an executable irrigation schedule, which explicitly specifies how much irrigation water should be applied on which dates.

[0036] Example 2: This example details the specific calculation methods and physical meanings of the key independent variables used to achieve dynamic calculations in the crop growth model, namely the crop growth stage indicators.

[0037] Step 201: During the operation of the crop growth model, crop growth stage indicators corresponding to the crop growth stages are calculated based on the simulated crop development process. These crop growth stage indicators are then substituted into one or more pre-defined continuous functions used to characterize stress response to determine the root water uptake stress under the current growth stage.

[0038] In this embodiment, the biologically discrete or qualitative growth stages of crops, such as seedling stage, jointing stage, flowering stage, and grain-filling stage, are transformed into continuously changing, dimensionless numerical indicators. These indicators can accurately reflect the developmental process of crops at any given time and serve as the core independent variable for all subsequent dynamic stress response functions, enabling the model to adjust its response to environmental stresses in a continuous and smooth manner.

[0039] Step 202, the crop growth stage index, is calculated through the following steps: based on the daily average temperature in the meteorological data and the preset basic temperature for crop development in the crop data, the effective accumulated temperature is accumulated daily; and the daily accumulated effective accumulated temperature is normalized according to the preset cumulative effective accumulated temperature for the entire growth period in the crop data to obtain the crop growth stage index.

[0040] This embodiment provides a specific and preferred implementation method. Specifically, the calculation process includes two core steps: effective accumulated temperature accumulation and normalization processing.

[0041] Specifically, this involves daily calculation of the effective accumulated temperature. Effective accumulated temperature is a core environmental driver guiding crop growth and development; its physical meaning is the accumulation of heat exceeding the minimum temperature threshold required for crop development. The crop development baseline temperature is a preset parameter in the crop data, representing the minimum temperature at which a specific crop variety begins effective physiological development. For example, for winter wheat, this baseline temperature can be set to 0 degrees Celsius. The formula for calculating the daily effective accumulated temperature can be expressed as:

[0042] GDD daily (i)=max(T avg (i)-T base ,0);

[0043] Among them, GDD daily (i) represents the daily effective accumulated temperature on day i within the simulation period, in degrees Celsius per day; max() is the function to find the maximum value; T avg (i) represents the daily average temperature of day i, which can be obtained by averaging the daily maximum and minimum temperatures from the meteorological data. base This is the preset base temperature for crop development.

[0044] During the model simulation, GDD is calculated daily starting from the crop sowing date. daily The accumulated temperature is then summed to obtain the cumulative effective accumulated temperature up to the current simulation date. The calculation formula can be expressed as:

[0045] ;

[0046] Among them, GDD d Let be the cumulative effective accumulated temperature up to day d, where i is the date index from the sowing date, d is the date index of the current simulation day, and ∑ is the summation function.

[0047] After obtaining the accumulated effective temperature, it is normalized. This normalization process eliminates differences in total accumulated temperature requirements among different crop varieties or varieties with different maturity periods, resulting in a universally applicable developmental progress index ranging from 0 to 1. The accumulated effective temperature over the entire growth period is also a pre-defined parameter in the crop data, representing the total heat required for a specific crop variety from sowing to physiological maturity. The normalization formula can be expressed as:

[0048] GS d =GDD d / GDD total ;

[0049] Among them, GS d GDD represents the crop growth stage indicator up to day d. total This is the accumulated effective temperature for the entire reproductive period as preset.

[0050] As a concrete example, suppose the basal temperature T for the crop development of a certain winter wheat variety is... base The cumulative effective temperature (GDD) over the entire growth period is 0 degrees Celsius. total The temperature is 1600 degrees Celsius per day. On a certain day in the simulation, designated day 60, the average daily temperature is 15 degrees Celsius, and the cumulative effective accumulated temperature from sowing to day 59 is 385 degrees Celsius per day. On day 60, calculate the effective accumulated temperature GDD for that day. daily (60) = max(15-0,0) = 15 degrees Celsius per day. Update the cumulative effective accumulated temperature (GDD). 60 =385 + 15 = 400 degrees Celsius per day. Calculate the crop growth stage index GS for that day. 60 =400 / 1600=0.25. This result is from GS. d =0.25 precisely calculates that at this moment, the crop has completed 25 percent of the heat accumulation required for its entire life cycle.

[0051] In some alternative implementations, the crop development basal temperature T base and cumulative effective accumulated temperature (GDD) throughout the reproductive period total Instead of fixed values, segmented parameter values ​​can be used based on different developmental stages of the crop, such as the vegetative growth stage and the reproductive growth stage, to further improve the accuracy of the development process simulation.

[0052] Furthermore, it should be noted that the method of normalization using effective accumulated temperature is the preferred method of this invention, because this method can more accurately reflect the actual biological development process of crops in different years and under different meteorological conditions than simply using the number of days in the growing season, thus providing a more reliable basis for the accurate calculation of subsequent dynamic stress.

[0053] Example 3: This example details a dynamic stress calculation framework that calculates water stress and salt stress separately, introduces salt tolerance and interaction terms that change dynamically with the growth stage, and finally merges them.

[0054] Understandably, in some existing crop models, root water uptake stress is obtained by multiplying the water stress reduction coefficient by the salt stress reduction coefficient. The limitation of this method is that the parameters describing salt stress, such as the salt tolerance threshold, are typically set as fixed constants throughout the entire growth period, and it fails to adequately consider the complex synergistic stress effects between water and salt stress that vary with crop development stages. The technical solution proposed in this embodiment aims to address the above problems.

[0055] Step 301: Determine the root water absorption stress under the current growth stage, specifically through a separated stress framework, such as... Figure 2 As shown, the separated stress framework includes: independently calculating the water stress reduction coefficient based on the soil moisture state simulated by the crop growth model; and independently calculating the salt stress reduction coefficient based on the soil salinity state and crop growth stage indicators simulated by the crop growth model, wherein the calculation process of the salt stress reduction coefficient reflects the influence of the growth stage; and merging the water stress reduction coefficient and the salt stress reduction coefficient to obtain the final root water uptake stress.

[0056] In this embodiment, the decomposition method refers to breaking down the total root water uptake stress into multiple components caused by different stress sources, mainly soil moisture deficit and excessive soil salinity. Each component is modeled and calculated independently, and then the components are merged in a specific way.

[0057] Specifically, the calculation of the water stress reduction coefficient can employ the Feddes model, which is well-known in the art. In this embodiment, the basic method for combining the water stress reduction coefficient and the salinity stress reduction coefficient is multiplicative merging, i.e., α total =α w ×α s , where α w Let α be the water stress reduction coefficient. s This represents the salt stress reduction coefficient. This merging method is a conventional technique in the field for handling multi-source stress. In subsequent steps, a salt-drought interaction correction term will be introduced to further enhance this basic merging method. The model calculates the reduction coefficient, ranging from 0 to 1, based on the soil water potential simulated in real time by the crop growth model and compared with a set of preset water potential threshold parameters.

[0058] Step 302: Based on the soil salinity status and crop growth stage indicators simulated by the crop growth model, independently calculate the salinity stress reduction coefficient, such as... Figure 3As shown, the specific steps include: using crop growth stage indicators as independent variables and inputting them into a first preset continuous function to calculate the salt tolerance threshold that dynamically changes with the growth stage; and using the dynamically changing salt tolerance threshold and soil salinity status as inputs into the salt stress response function to calculate the salt stress reduction coefficient.

[0059] In this embodiment, the salt stress response function is preferably the van Genuchten-Hoffman function, which can be expressed in the following form:

[0060] α s =1 / [1+(ECe / EC 50 (GS)) p ];

[0061] Where, α s The salt stress reduction coefficient ranges from 0 to 1. ECe represents the current root zone soil saturated extract conductivity simulated by the crop growth model, in decibels per meter. 50 (GS) is the salt tolerance threshold dynamically calculated based on the first preset continuous function. Its physical meaning is the soil salt concentration that causes a 50% decrease in crop yield. p is an empirical shape parameter used to control the rate at which yield decreases with increasing salt content. It is usually between 2 and 5 and is a pre-configured model parameter.

[0062] This embodiment details the core mechanism for dynamic calculation of salt stress in this invention. Specifically, a crop's salt tolerance changes dynamically throughout its life cycle; for example, it is typically more sensitive during the seedling stage, while tolerance increases at certain growth stages.

[0063] The specific implementation of this mechanism involves no longer using a fixed salt tolerance threshold, but instead mapping the crop growth stage index GS to a dynamically changing salt tolerance threshold EC through a preset continuous function. 50 (GS). The dynamic threshold EC calculated in real time will be used. 50 (GS), together with the current root zone soil salinity state simulated by the crop model, such as the soil saturated extract conductivity ECe, are substituted into the salt stress response function, such as the van Genuchten-Hoffman function, to calculate the final salt stress reduction coefficient.

[0064] Step 303: The first preset continuous function is a polynomial function or a radial basis function. Through a set of pre-configured coefficients or parameters, the crop growth stage index is mapped to a dynamically changing salt tolerance threshold.

[0065] This embodiment provides two specific, optional implementation methods for calculating the first preset continuous function of the dynamic salt tolerance threshold.

[0066] The first implementation method is to use a polynomial function. For example, a quadratic polynomial can be used to describe the relationship between the salt tolerance threshold and the gestational age index (GS), and its formula can be expressed as:

[0067] EC 50 (GS)=p1*GS 2 +p2*GS+p3;

[0068] Among them, EC 50 (GS) is the dynamically calculated salt tolerance threshold, in decimens per meter. GS is a dimensionless crop growth stage index. p1, p2 and p3 are model parameters predetermined by parameter calibration methods.

[0069] The second implementation, and also a preferred method, uses radial basis functions (RBF). Compared to polynomial functions, RBF can more flexibly fit complex nonlinear relationships, such as the phenomenon of crops exhibiting specific high tolerance windows during the mid-growth stage. Furthermore, parameter constraints can ensure that their output values ​​remain positive, better reflecting the physical meaning of salt tolerance thresholds. The form of using a superposition of multiple Gaussian RBF functions can be expressed as:

[0070] ;

[0071] Among them, EC 50 (GS) is the dynamically calculated salt tolerance threshold, EC base Where A is the basic salt tolerance threshold, N is the total number of radial basis functions used, and A is the base salt tolerance threshold. i Let μ be the amplitude coefficient of the i-th radial basis function, exp() be the natural exponential function, GS be the crop growth stage index, and μ be the amplitude coefficient of the i-th radial basis function. i Let σ be the center position of the i-th radial basis function on the reproductive stage axis. i Let A be the width parameter of the i-th radial basis function. i μ i and σ i These are all pre-configured model parameters. For example, when N=1, this function can simulate μ. i Near the reproductive stage it represents, an amplitude of A appears. i Improved salt tolerance.

[0072] As a concrete numerical example, suppose a single Gaussian radial basis function (N=1) is used to describe the dynamic change of the salt tolerance threshold of a certain winter wheat variety. Based on the parameter calibration method, the following set of typical parameters can be obtained: EC base =4.0 dS / m, A1=3.0 dS / m, μ1=0.50, σ1=0.15. The physical meaning of this parameter set is: throughout the entire growth period, the crop's basic salt tolerance threshold is 4.0 dS / m, while in the mid-growth stage represented by μ1=0.50, the tolerance threshold reaches a peak, reaching a maximum of EC... base +A1=7.0dS / m, reflecting the physiological characteristics of enhanced crop osmotic regulation capacity during this stage.

[0073] Based on the above parameters, the calculation process for three typical reproductive stages is shown:

[0074] During the seedling stage, when GS=0.15: EC 50 (0.15) = 4.0 + 3.0 × exp(-(0.15 - 0.50) 2 / (2×0.15 2 =4.0+3.0×exp(-2.72)=4.0+3.0×0.066=4.20dS / m, indicating that the salt tolerance of seedlings is low, close to the baseline level.

[0075] During mid-fertility, when GS=0.50: EC 50 (0.50) = 4.0 + 3.0 × exp(-(0.50 - 0.50) 2 / (2×0.15 2 The value of 4.0 + 3.0 × 1.0 = 7.0 dS / m indicates that the salt tolerance reaches its peak at this stage.

[0076] During the grouting period, when GS=0.80: EC 50 (0.80) = 4.0 + 3.0 × exp(-(0.80 - 0.50) 2 / (2×0.15 2 =4.0+3.0×exp(-2.0)=4.0+3.0×0.135=4.41dS / m, indicating that the salt tolerance during the grouting period has dropped back to near the base level.

[0077] It can be seen that, compared to using a fixed salt tolerance threshold, such as an average of approximately 5.2 dS / m over the entire growth period, the dynamic salt tolerance threshold of this invention can accurately capture the differences in salt tolerance at different growth stages. This difference will directly affect the calculation of the salt stress reduction coefficient. For example, under the same soil salinity conditions at the seedling stage, the dynamic threshold will calculate a stronger salt stress, prompting the optimization algorithm to allocate more irrigation and leaching water at this stage. This aligns with the experience in agronomic practice of emphasizing salt prevention during the seedling stage. It should be noted that the above numerical examples are used to illustrate the working principle of the dynamic stress calculation mechanism of this invention and its differences from traditional fixed-parameter methods. In practical applications, those skilled in the art can calibrate and verify the improved model of this invention and the traditional model using fixed parameters on the same field observation dataset, respectively, by comparing the simulation accuracy indicators (such as root mean square error RMSE, coefficient of determination R) of the two models on key variables such as soil moisture, salinity, and crop yield. 2 The accuracy improvement resulting from the dynamic stress calculation mechanism of this invention is calculated and evaluated using this method. Those skilled in the art will understand that the specific improvement may vary depending on crop variety, soil type, and climatic conditions.

[0078] Step 304, the separated stress framework further includes: taking crop growth stage indicators as independent variables and inputting them into a second preset continuous function to calculate the interaction intensity coefficient that dynamically changes with the growth stage; calculating the salt-drought interaction correction term based on the interaction intensity coefficient, water stress reduction coefficient, and salt stress reduction coefficient; wherein, the water stress reduction coefficient and the salt stress reduction coefficient are combined, specifically, the water stress reduction coefficient, the salt stress reduction coefficient, and the salt-drought interaction correction term are combined to obtain the final root water uptake stress.

[0079] This embodiment introduces a dynamic description of the salt-drought interaction in the separate stress framework of the present invention, which solves the problem that the traditional simple multiplication merging method may underestimate the synergistic stress effect.

[0080] Specifically, the interaction strength coefficient that dynamically changes with the fertility stage index GS is calculated using a second preset continuous function, such as a Gaussian function. The formula can be expressed as:

[0081] γ int (GS)=γ max *exp(-(GS-GS peak ) 2 / (2*w 2 ));

[0082] Where, γ int (GS) is the dynamically calculated interaction strength coefficient, whose value range is typically between 0 and 1, γ max For maximum interaction strength, GSpeak γ represents the reproductive stage with the strongest interaction, w is the width parameter describing the duration of the reproductive stage sensitive to the interaction, and γ is the reproductive stage index. max GS peak Both and w are pre-configured model parameters.

[0083] Based on the interaction strength coefficient and the independently calculated reduction coefficients for water stress and salinity stress, a salt-drought interaction correction term is calculated. This correction term aims to apply an additional reduction effect when water and salinity stress coexist. Its calculation formula can be expressed as:

[0084] α int =1-γ int (GS)*(1-α w )*(1-α s );

[0085] Where, α int For the salt-drought interaction correction term, α w α is the independently calculated water stress reduction coefficient. s This represents the salt stress reduction coefficient.

[0086] The three components are combined to obtain the final root water uptake stress. The combination method is to multiply the three components, and the calculation formula is as follows: α total =α w *α s *α int ;

[0087] Where, α total This is the comprehensive stress reduction coefficient ultimately used for root water uptake. In this way, the present invention not only achieves the dynamics of salt stress itself, but also the dynamics of the salt-drought interaction, enabling more accurate simulation of crop growth responses under complex stress environments.

[0088] As a concrete numerical example, let's illustrate the effect of the interaction correction term. Assume that during mid-fertility (GS=0.50), the interaction strength coefficient γ... int (0.50) Based on the Gaussian function above, it is calculated to be 0.6, assuming γ max =0.6, GS peak =0.50, and the model simulation yielded a water stress reduction coefficient α. w =0.7, salt stress reduction coefficient α s =0.8.

[0089] Calculate the salt-drought interaction correction term: α int =1-0.6×(1-0.7)×(1-0.8)=1-0.6×0.3×0.2=1-0.036=0.964.

[0090] Calculate the final overall stress reduction factor: α total =α w ×α s ×α int =0.7×0.8×0.964=0.540.

[0091] In contrast, if the traditional simple multiplication merging method is used, without interactive correction terms, then: α total_simple =α w ×α s =0.7×0.8=0.560.

[0092] In this example, the interaction correction term further reduces the final stress from 0.560 to 0.540, meaning that at this reproductive stage, when water and salt stress are present simultaneously, their synergistic effect leads to an additional decrease in water uptake capacity of approximately 3.6%. This synergistic effect is particularly significant in the mid-reproductive stage, when the interaction is strongest, while it is weaker or negligible in other stages.

[0093] Example 4: This example serves as an alternative and preferred implementation of the dynamic stress calculation mechanism. It should be noted that the separate stress framework relies somewhat on empirical methods when merging the various stress components. This example provides a more physically clear framework, unifying soil matrix potential (representing drought stress) and osmotic potential (representing salinity stress) under a single effective water potential variable.

[0094] Step 401: Determine the root water absorption stress under the current growth stage, specifically through a unified effective water potential framework, such as... Figure 4 As shown, the unified effective water potential framework includes: calculating the osmotic potential based on the soil salinity state simulated by the crop growth model; determining the osmotic sensitivity coefficient that dynamically changes with the growth stage based on the crop growth stage index and through a preset continuous function; weighting the osmotic potential using the osmotic sensitivity coefficient and merging the weighted osmotic potential with the soil matrix potential simulated by the crop growth model to construct the unified effective water potential; and calculating the final root water uptake stress based on the unified effective water potential.

[0095] In this embodiment, it is necessary to convert soil salinity status indicators, such as soil saturated extract conductivity ECe, which are typically output by crop growth models, into water potential units with the same dimensions, i.e., osmotic potential. This conversion can be achieved based on empirical relationships known in the art, for example:

[0096] h o =-36×ECe;

[0097] Among them, h oThe calculated osmotic potential is expressed in centimeters of water column (cmHg). ECe is the electrical conductivity of the saturated soil extract, expressed in centi-Siemens per meter (dSm). -36 is an empirical conversion factor known in the art, expressed in centimeters of water column (dSm) per meter (dSm), which incorporates the dimensional conversion from electrical conductivity to water potential.

[0098] After obtaining the osmotic potential h o and the soil matrix potential h obtained by model simulation m Then, the two are merged to construct a unified effective water potential. Simple linear addition cannot reflect complex plant physiological responses, which will be explained in detail in subsequent embodiments.

[0099] Step 402: The osmotic potential is merged with the soil matrix potential simulated by the crop growth model to construct a unified effective water potential. Specifically, this includes: taking the crop growth stage index as an independent variable and inputting it into a third preset continuous function to calculate the osmotic sensitivity coefficient that dynamically changes with the growth stage; using the osmotic sensitivity coefficient to weight the osmotic potential to obtain a weighted osmotic potential; and merging the weighted osmotic potential with the soil matrix potential to obtain a unified effective water potential.

[0100] This embodiment focuses on the plant's inherent osmotic regulation capacity, which involves accumulating solutes within root cells to lower its own water potential, counteracting the high osmotic potential of the external soil, and maintaining root water absorption. The strength of this osmotic regulation capacity varies at different growth stages of the crop. To represent this physiological mechanism in the model, this embodiment introduces the concept of an osmotic sensitivity coefficient.

[0101] Specifically, the infiltration sensitivity coefficient is a dimensionless coefficient that dynamically changes with the crop growth stage index GS, used to calculate the actual response of the crop to soil infiltration potential at a specific growth stage. This coefficient is calculated using a third preset continuous function, designed to reflect the physiological characteristics of default high sensitivity and the emergence of a tolerance window at a specific stage. For example, it can be constructed using a baseline of 1 superimposed with a negative Gaussian function.

[0102] ω(GS)=1-A ω *exp(-(GS-GS peak_ω ) 2 / (2*w ω 2 ));

[0103] Where ω(GS) is the permeability sensitivity coefficient under GS during the reproductive stage, and its value range is usually between 0 and 1. ω The amplitude of the maximum osmotic regulation capacity, GS peak_ω This is the reproductive stage where osmotic regulation is strongest. ω A is the duration width parameter for the capability phase of this emphasis section. ω GSpeak_ω and w ω All parameters are pre-configured model parameters. When ω(GS) is 1, it indicates that the crop is sensitive and cannot perform osmotic regulation; when it is less than 1, it indicates that the crop has performed partial osmotic regulation, effectively offsetting part of the soil osmotic potential.

[0104] Once the dynamic permeability sensitivity coefficient is obtained, a unified effective water potential can be constructed, and its calculation formula is as follows:

[0105] h eff =h m +ω(GS)*h o ;

[0106] Among them, h eff To obtain the final unified effective water potential, h m h represents the soil matrix potential simulated by the model. o This represents the soil permeability potential.

[0107] As a concrete example, suppose that during the mid-growth stage, GS is 0.6, at which point the crop's osmotic regulation capacity is strongest, and ω(0.6) is calculated to be 0.5 according to the above function. If the model simulates the soil matrix potential h at this moment... m At a depth of -600 cm, the soil infiltration potential h o The value is -500 cm. Therefore, the uniform effective water potential h actually experienced by the crop is... eff =-600+0.5×(-500)=-850 cm. It can be seen that this result is not a simple arithmetic sum of matrix potential and osmotic potential, but rather because the crop's strong osmotic regulation ability effectively alleviates the stress it feels.

[0108] Step 403: Based on the unified effective water potential, the final root water uptake stress is calculated. Specifically, the unified effective water potential is used as the only independent variable and substituted into a single, preset stress response function to calculate the final root water uptake stress.

[0109] This embodiment unifies the effects of two stressors, drought and salinity, into a single variable h. eff Subsequently, the calculation of the final stress reduction coefficient no longer requires complex empirical merging, but instead adopts a standard, single stress response function.

[0110] In this embodiment, the single, pre-defined stress response function can preferably be a FedDes model function. The independent variable input to this function is no longer the traditional soil matrix potential h. m Instead, it unifies the effective water potential h. eff The specific form of this function can be a piecewise function as follows:

[0111] If heff If >= h1, then α total = 0;

[0112] If h2 < h eff < h1, then α total = (h eff - h1) / (h2 - h1);

[0113] If h3 <= h eff <= h2, then α total = 1;

[0114] If h4 < h eff < h3, then α total = (h eff - h4) / (h3 - h4);

[0115] If h eff <= h4, then α total = 0;

[0116] Where, α total is the final root water absorption stress reduction coefficient, h eff is the unified effective water potential, h1 is the upper water potential threshold at which crop root water absorption stops, h2 and h3 are the optimal water potential intervals for optimal water absorption of the crop, h4 is the lower water potential threshold at which root water absorption completely stops, usually corresponding to the permanent wilting point. This set of threshold parameters are pre-configured model parameters.

[0117] In this embodiment, a unified effective water potential is constructed through a dynamic osmotic sensitivity coefficient, and combined with a single stress response function, a dynamic stress calculation method with clearer physical mechanism and more intrinsic parameter coupling is provided compared to the separated framework.

[0118] Embodiment 5. As a dynamic stress calculation mechanism, this embodiment provides a specific and preferred implementation method at the numerical calculation level. It clarifies that the dynamic stress calculation logic is deeply and synchronously coupled with the numerical solver of the physical equation at the bottom layer of the crop growth model, ensuring the physical self-consistency and numerical stability of the simulation results.

[0119] It can be understood that most crop growth models based on physical processes solve one or more partial differential equations. For example, the Richards equation used to describe soil water movement. Since this equation is non-linear, in the process of advancing from the current time step to the next time step, an iterative algorithm is usually required within the time step, such as the Picard iteration method or the Newton-Raphson iteration method, and the calculation is repeated multiple times until the solution converges. Root water absorption is expressed as a sink term in the Richards equation, and its magnitude is affected by the root water absorption stress reduction coefficient.

[0120] A simplified coupling approach involves calculating root water uptake stress at the beginning of each time step, based on the soil water and salt state at the end of the previous time step, and then applying this fixed stress value to all numerical iterations within the current time step. This method ignores the fact that the soil water and salt state itself undergoes drastic changes during the iterations of the current time step, while stress should respond to these changes in real time. This temporal lag introduces physical inconsistencies, especially during events that drastically alter soil moisture conditions, such as irrigation or rainfall, potentially leading to distorted simulation results or numerical oscillations.

[0121] Step 501, determining the root water uptake stress under the current growth stage, is performed within the numerical iterative loop of the crop growth model used to solve the hydrodynamic equation. Specifically, it includes: in each iteration of the numerical iterative loop, obtaining the soil water potential and soil salinity status under the current iteration; calculating the root water uptake stress under the current iteration based on the soil water potential, soil salinity status, and crop growth stage indicators; and using the root water uptake stress under the current iteration as the pool term of the hydrodynamic equation to solve the updated soil water potential in the same iteration, until the numerical iterative loop converges.

[0122] This embodiment details the iterative internal coupling scheme adopted by the present invention to solve the aforementioned technical problems. This scheme seamlessly embeds the entire dynamic stress calculation process into the numerical iterative loop of the flow motion equation solver.

[0123] As a concrete example, taking the time step of solving the Richards equation using the Picard iteration method as an example, the detailed execution flow of the method of this invention is as follows:

[0124] At the start of the simulation at time step t, the soil water potential profile at time step t+Δt is solved. The solver begins its first iteration, denoted as iteration number k=1.

[0125] At the start of the k-th iteration, obtain the soil water potential h under the current iteration. m The estimated values ​​of (k) and soil salinity state ECe(k). For the first iteration, this value can be taken from the result at the end of the previous time step.

[0126] Based on the current iteration h m Using (k) and ECe(k), along with the crop growth stage index GS for the day (which remains constant within a time step), a complete dynamic stress calculation is performed. Specifically, if a separate framework is used, water stress, dynamic salinity stress, and dynamic interaction terms are recalculated and merged; if a unified effective water potential framework is used, the unified effective water potential is recalculated and substituted into the stress function. Regardless of the framework used, this calculation yields a root uptake stress reduction coefficient synchronized with the current iteration state, denoted as α. total(k), the latest calculated stress reduction coefficient α total (k) is used to update the root water absorption sink term in the Richards equation.

[0127] Solving the updated system of equations yields a new estimate of soil water potential, denoted as h. m (k+1).

[0128] Check convergence. For example, you can calculate h. m (k+1) and h m (k) The difference norm across the entire soil profile. If this difference is less than a preset convergence threshold, the calculation at the current time step is considered converged, the iteration loop ends, and h is set. m (k+1) is taken as the final result. If convergence fails, let k = k+1, and set h... m (k+1) is used as the input for the next iteration, and all the above steps are repeated.

[0129] By placing stress calculations within an iterative loop, this invention ensures that in each iteration, the root water uptake sink term used to calculate water flow accurately reflects the soil water and salt conditions under that iteration. In other words, soil water potential (which determines the cause of stress) and root water uptake (the result of stress influence) are numerically synchronized, coupled, and self-consistent. This approach enhances the physical realism of the model and the robustness of numerical calculations, enabling more accurate capture of the rapid dynamic feedback process between root water uptake and the soil environment.

[0130] Example 6: This example describes an optimization strategy based on a surrogate model, providing a specific and efficient implementation strategy for obtaining the optimal irrigation scheme through a multi-objective optimization algorithm.

[0131] The multi-objective optimization problem in this embodiment is explicitly defined mathematically: it involves finding an optimal set of decision variables while simultaneously optimizing multiple conflicting objective functions.

[0132] The decision variable is the irrigation scheme, which can be represented as a multidimensional vector I=[i1,i2,...,i...]. D ], where D is the preset number of irrigations during the entire crop growth period, i k Let be the irrigation water volume for the k-th irrigation event.

[0133] In this embodiment, the multiple optimization objectives may specifically include the following three:

[0134] Objective 1: Minimize total irrigation water consumption. This objective aims to conserve water resources, and its mathematical expression is:

[0135] f1(I)=min(∑ k=1 D i k);

[0136] Objective two: Minimize crop yield loss. This objective is equivalent to maximizing relative crop yield, aiming to ensure food security. Its mathematical expression is:

[0137] f2(I)=min(1-Y a (I) / Y p );

[0138] Among them, Y a (I) represents the actual yield, Y, obtained from crop growth model simulation under irrigation scheme I. p This represents the potential yield obtained from model simulation under ideal water and fertilizer conditions.

[0139] Objective 3: Minimize soil salinity in the root zone at the end of the season. This objective aims to control and improve soil salinization, ensuring the sustainable use of land resources. Its mathematical expression is:

[0140] f3(I)=min(EC rz_end (I));

[0141] Among them, EC rz_end (I) represents the average soil salinity in the root zone at harvest time, as simulated by the crop growth model under irrigation scheme I.

[0142] It is understandable that coupling multi-objective optimization algorithms, such as the second generation of non-dominated sorting genetic algorithms NSGA-II or NSGA-III, with crop growth models presents technical challenges. This is because a single run of a crop growth model typically takes several seconds to several minutes, while the optimization algorithm requires tens of thousands or even more model calls to obtain a convergent solution set. The high computational cost makes direct coupling virtually impossible in practice. To address these technical problems, this embodiment employs a surrogate model strategy.

[0143] Step 601 involves obtaining the optimal irrigation scheme through a multi-objective optimization algorithm. Specifically, this includes: constructing a surrogate model based on the crop growth model to approximate the input-output relationship of the crop growth model. The input of the surrogate model is the irrigation scheme, and the output is the corresponding simulation result. The multi-objective optimization algorithm is then applied to the surrogate model to obtain the optimal irrigation scheme.

[0144] In this embodiment, the surrogate model, also known as the alternative model or meta-model, is a low-computational-cost mathematical or statistical model used to learn and simulate the behavior of the high-cost original crop growth model. The input of the surrogate model is the same as that of the original model, namely the irrigation scheme vector I; its output is a rapid prediction of the simulation results of the original model, namely, the estimated values ​​of the three optimization objectives [f1(I), f2(I), f3(I)].

[0145] As a preferred implementation, the surrogate model can be a Kriging model, also known as a Gaussian process regression model. This model not only provides predicted values ​​for unknown points, but also gives the uncertainty or confidence level of those predicted values.

[0146] The core mathematical idea of ​​the Kriging model can be briefly described as follows: For an unknown irrigation scheme I*, the predicted value and prediction uncertainty of the simulation results are expressed as follows:

[0147] =μ+r T (I*)·R (-1) ·(y-μ·1);

[0148] σ²(I*)=σ²·[1-r T (I*)·R (-1) ·r(I*)];

[0149] in, Let σ²(I) be the predicted value of the simulation result for irrigation scheme I, used to quantify the uncertainty of the prediction. μ is the global mean parameter, σ² is the process variance parameter, y is the vector of the true simulation results in the existing training dataset, l is a one-dimensional vector, R is the correlation matrix between the sample points in the training dataset, whose elements are calculated by a parameterized correlation function (e.g., a Gaussian correlation function), and r(I*) is the correlation vector between the unknown point I and each training sample point. The model's hyperparameters, including μ, σ², and the correlation function parameters, can be automatically determined from the training dataset using the maximum likelihood estimation method. The key advantage of this model is that the value of σ²(I) is larger in regions far from the training samples and smaller in regions with dense training samples, providing a direct measure of informational value for subsequent adaptive point addition strategies.

[0150] Specifically, the optimization process of this invention is transformed into: generating a set of input-output data pairs for training by running a high-cost original crop growth model a finite number of times; constructing a Kriging surrogate model based on this training dataset; and combining a multi-objective optimization algorithm with the surrogate model to solve the problem. Since the time cost of a single call to the surrogate model is extremely low, the optimization algorithm can complete a sufficient number of iterations within an acceptable time to find a high-quality optimal solution set.

[0151] Step 602, constructing the surrogate model, specifically through a sequential adaptive addition strategy, includes: generating an initial training dataset and constructing an initial surrogate model; and iteratively executing the following steps until a preset termination condition is met: based on the predictions of the current surrogate model, using preset addition criteria, determining the next addition scheme from the candidate irrigation schemes defined by the decision variable space of the irrigation schemes that can maximize the improvement of the multi-objective optimization solution set; running the crop growth model to obtain the real simulation results corresponding to the next addition scheme; and adding the next addition scheme and its real simulation results to the training dataset to update the surrogate model.

[0152] This embodiment provides a more efficient and intelligent implementation of the proxy model construction process, aiming to build the most accurate proxy model possible with the fewest calls to the original model. This strategy includes two phases: initialization and iterative enhancement.

[0153] During the initialization phase, a small initial training dataset is generated using a space-filling design method. Preferably, Latin Hypercube Sampling (LHS) can be employed because it ensures that the sampling points are relatively evenly distributed across each dimension of the decision variables. For example, for an optimization problem involving 5 irrigation decisions, an initial design with 50 sample points is generated. The original crop growth model is run 50 times to obtain the actual simulation results corresponding to 50 irrigation schemes, and an initial Kriging surrogate model is constructed based on this data.

[0154] During the iterative enhancement phase, the following operations are performed cyclically:

[0155] Based on the existing surrogate model, a pre-defined addition criterion is used to guide the selection of the next sample point. Preferably, this criterion is the Expected Hypervolume Improvement (EHVI). The physical meaning of this criterion is to comprehensively consider the surrogate model's predictions for unknown points and their uncertainties, to find the candidate point that has the greatest potential to expand the hypervolume dominated by the current Pareto optimal solution set. This candidate point is usually located at the forefront of the current optimal solution set or in a region with very high model uncertainty; it is both a deep optimization of the optimal region and an exploration of the unknown space.

[0156] In the entire space of candidate irrigation schemes, the point that maximizes the EHVI value is found through an optimization algorithm, and this point is determined as the next addition scheme.

[0157] For the selected addition scheme, a high-cost original crop growth model is invoked once to obtain its accurate and realistic simulation results.

[0158] Add new input-output data pairs to the existing training dataset and retrain or update the Kriging surrogate model using the updated dataset to further improve its accuracy.

[0159] The above iterative process will continue until a preset termination condition is met, such as the total number of calls to the original model reaching a preset computational budget limit, or the improvement in EHVI being less than a very small threshold. After the iteration ends, the final proxy model obtained will be used for subsequent optimization solutions.

[0160] Example 7: This example details how, after obtaining a Pareto optimal solution set containing multiple non-dominated solutions through multi-objective optimization, a systematic multi-criteria decision analysis method is used to select a unique, balanced compromise solution as the final executable irrigation scheme recommended to the user.

[0161] The output of a multi-objective optimization algorithm is a Pareto optimal solution set. It's important to note that every solution in this set, i.e., every irrigation scheme, is mathematically optimal; no other solution is superior to it across all optimization objectives. However, these solutions represent different trade-off strategies. For example, scheme A in the solution set might achieve acceptable yields with very low irrigation water usage, but lead to slight salt accumulation at the end of the season; while scheme B leaches salt through slightly more irrigation, ensuring soil health, but at the cost of increased water consumption. Faced with a set of choices containing multiple trade-offs, decision-makers need clear and computable methods to assist them in making their final choice.

[0162] Step 701: The optimal irrigation scheme obtained by solving the multi-objective optimization algorithm is the Pareto optimal solution set. The method also includes: using the ideal point approximation sorting method to select a compromise scheme from the Pareto optimal solution set as the final output optimal irrigation scheme.

[0163] In this embodiment, the ideal point approximation ranking method, or TOPSIS (Technique for Order of Preference by Similarity to Ideal Solution), is used as a multi-criteria decision analysis tool for selection. The optimal compromise should be geometrically closest to the positive ideal solution and farthest from the negative ideal solution. The following describes in detail the specific process of selecting the final solution from the Pareto optimal solution set using the TOPSIS method:

[0164] Construct a decision matrix. Assume the Pareto optimal solution set contains m irrigation schemes, each with n optimization objectives. In this example, n=3, representing total irrigation water volume, yield loss, and end-of-season salinity, respectively. An m x n decision matrix X can be constructed, where elements x... ijThis represents the performance value of the i-th irrigation scheme on the j-th optimization objective.

[0165] The decision matrix needs to be normalized. Since the three optimization objectives have different dimensions and numerical ranges, direct comparison is inappropriate. Therefore, the matrix needs to be normalized to eliminate the influence of dimensions. A commonly used normalization method is vector normalization, and its calculation formula is as follows:

[0166] r ij =x ij / sqrt(∑ i=1 m (x ij 2 ));

[0167] Where, r ij For the elements of the normalized decision matrix, x ij represents the elements of the original decision matrix, and m represents the total number of alternatives.

[0168] Construct a weighted normalized decision matrix. To reflect the decision-maker's preference for different optimization objectives, a weight w can be assigned to each optimization objective. j , where ∑w j =1. For example, in extremely water-scarce regions, policymakers may assign a higher weight to total irrigation water targets. The weighted normalization calculation is as follows:

[0169] v ij =w j *r ij ;

[0170] Among them, v ij These are the elements of the weighted normalized decision matrix.

[0171] Determine the positive and negative ideal solutions. A positive ideal solution A+ is a virtual optimal solution composed of the optimal values ​​of all optimization objectives in the decision matrix. A negative ideal solution A- is a virtual worst solution composed of the worst values ​​of all optimization objectives. Specifically:

[0172] A+={min i (v i1 ),min i (v i2 ),min i (v i3 )};

[0173] A-={max i (v i1 ),max i (v i2 ),max i (v i3 )};

[0174] Since all three objectives in this invention are cost-related objectives, the smaller the value, the better. Therefore, the positive ideal solution takes the minimum value of each column, and the negative ideal solution takes the maximum value of each column.

[0175] Calculate the distances from each solution to the positive and negative ideal solutions. Euclidean distance can be used as a metric. For the i-th solution, its distance S to the positive ideal solution... i+ and the distance S to the negative ideal solution i- The calculation formula is:

[0176] S i+ =sqrt(∑ j=1 n (v ij -v j+ ) 2 );

[0177] S i- =sqrt(∑ j=1 n (v ij -v j- ) 2 );

[0178] Among them, v j+ Let v be the j-th component of the positive ideal solution A+. j- Let j be the j-th component of the negative ideal solution A-.

[0179] Calculate and rank the relative closeness of each option. Relative closeness C i This indicates the degree to which the i-th solution is close to the ideal solution. Its value ranges from 0 to 1, and the larger the value, the better the solution.

[0180] C i =S i- / (S i+ +S i- );

[0181] After calculating the relative similarity C of all the schemes i Then, all plans were processed according to C. i The values ​​are sorted in descending order. The top-ranked solution is C. i The solution with the highest value is selected as the optimal compromise solution and becomes the unique and recommended irrigation solution in the final output of this invention.

[0182] In some alternative implementations, an interactive interface can be provided to decision-makers, allowing them to dynamically adjust the weights of each optimization objective. jThe system will recalculate and recommend the optimal compromise solution in real time based on the new weights, realizing a dynamic and personalized decision support process. Alternatively, other multi-criteria decision analysis methods, such as the Analytic Hierarchy Process (AHP) or the VIKOR method (a multi-criteria compromise ranking method), can be used to replace TOPSIS to achieve similar functionality.

[0183] Example 8: In this example, before irrigation optimization is performed, the key parameters inside the crop growth model, especially those related to the dynamic stress calculation mechanism, can be accurately determined, so that the model can accurately reflect the growth and water-salt response patterns of a specific plot and a specific crop variety.

[0184] It is understandable that the simulation accuracy of any crop growth model is highly dependent on the accuracy of its internal parameters. For the purposes of this invention, in addition to the original parameters of the crop model, parameters describing the dynamic stress response function, such as the dynamic salt tolerance threshold EC, are also included. 50 The function parameters of (GS), or the function parameters describing the dynamic permeability sensitivity coefficient ω(GS), have a decisive impact on the reliability of the final optimization results. Therefore, when applying the method of this invention to new regions or new crop varieties, these key parameters must be calibrated.

[0185] Step 801, the method further includes: before running the crop growth model, determining at least one preset parameter in the crop growth model through a model calibration method.

[0186] In this embodiment, the model calibration method is a systematic process aimed at finding a set of optimal model parameter values ​​by comparing them with real-world observation data. At least one preset parameter, preferably, is the new parameter introduced in this invention for calculating the dynamic stress response characteristics of crops.

[0187] Specifically, when using a split-type stress frame, the parameters that need to be calibrated may include: for EC 50 The polynomial function coefficients p1, p2, p3 of (GS); or used for EC 50 The radial basis function parameter A of (GS) i ,μ i ,σ i ; and γ for the salt-drought interaction correction term max GS peak,w .

[0188] When using a unified effective water potential framework, the main parameters that need to be calibrated are: the function parameter A used for the permeability sensitivity coefficient ω(GS). ω GS peak_ω ,w ω .

[0189] The following provides two optional implementation methods for determining this type of parameter.

[0190] The first approach is the literature data fitting method. This method is suitable for scenarios where detailed local experimental data is lacking, but relevant research literature is abundant. The specific process is as follows: Scientific literature is retrieved to collect and organize experimental data on salt tolerance of specific crops (e.g., winter wheat) at different growth stages. This data is typically presented as pairs of crop growth stages and salt tolerance thresholds. The qualitative growth stages in the literature are converted into crop growth stage indices (GS). Statistical fitting techniques, such as nonlinear least squares, are used to fit the function to be calibrated (e.g., radial basis function) to the data points extracted from the literature. The output of the fitting process is a set of optimal parameter values ​​after calibration.

[0191] The second implementation method, and also a preferred approach, is the model inversion calibration method based on local experimental data. This method can provide the most accurate and targeted parameters for specific plots and varieties. Its specific process is as follows:

[0192] A historical dataset for calibration needs to be prepared. This dataset typically comes from field control experiments over one or more growing seasons and must contain two types of information: first, complete model input data for that growing season, including daily meteorological data, detailed soil parameters, crop parameters, and precisely recorded irrigation and fertilization management practices; second, field observations of key state variables for that growing season, such as time-series observations of soil moisture content at different depths, final observations of crop biomass or yield, and measured values ​​of soil salinity profiles at the end of the season.

[0193] Define an objective function to calculate the model's simulation accuracy and guide parameter optimization. Preferably, this objective function can be set to minimize the root mean square error between the model's simulated values ​​and the field observations. This objective function comprehensively evaluates the performance of the selected parameter set in reproducing historical reality.

[0194] A global optimization algorithm is employed to automatically search for the optimal combination of parameters that minimizes the objective function. Since the relationship between model parameters and the objective function is often highly nonlinear and complex, robust global optimization algorithms are preferred. For example, the SCE-UA algorithm (Shuffled Complex Evolution - University of Arizona), Bayesian optimization, or Particle Swarm Optimization can be used.

[0195] The inversion calibration process is executed. In each iteration, the optimization algorithm automatically generates a set of candidate parameter values, inputs these values ​​into the crop growth model, and runs a complete simulation based on historical input data. Then, the root mean square error (RMSE) between the simulation results and historical observations is calculated. Based on the calculated RMSE, the optimization algorithm adjusts its search strategy and generates the next set of candidate parameters, repeating this process until a preset number of iterations is reached or the objective function converges to an acceptable low level. When the optimization process terminates, the parameter combination that minimizes the RMSE is determined as the final calibration result.

[0196] By performing the above-described offline calibration process, the crop growth model of the present invention can provide a highly reliable and accurate simulation basis for subsequent irrigation optimization.

[0197] Example 9: The process described in this example is an enhancement measure to improve the reliability of the final output scheme. The final irrigation scheme, which is optimized based on a computationally efficient surrogate model and selected through decision-making, is finally verified by a high-precision original crop growth model. This ensures that the prediction accuracy of the surrogate model is reliable in the region of the optimal solution and provides a systematic callback and correction mechanism when the prediction is inaccurate.

[0198] In practical applications, surrogate models are essentially mathematical approximations of complex original crop growth models. Although sequential adaptive addition strategies can improve accuracy, non-negligible prediction errors may still exist in certain complex response regions. The verification and callback process proposed in this embodiment is used to manage uncertainty and serves as a quality verification mechanism before the final solution is output.

[0199] The specific execution process of this embodiment may include the following steps:

[0200] Obtain the scheme to be verified and its predictive performance. The scheme to be verified is the unique optimal compromise irrigation scheme selected from the Pareto optimal solution set through the Ideal Point Approximation Ranking Method (TOPSIS) or other multi-criteria decision-making methods. At the same time, obtain the predictive performance value of this scheme on the surrogate model, that is, the values ​​of the three optimization objectives predicted by the surrogate model for it, denoted as the prediction objective vector.

[0201] Perform a high-precision baseline simulation. Using the proposed scheme as input for irrigation management, run a complete, high-precision, parameter-calibrated original crop growth model once. Since only one run is required, the time cost is acceptable. The output of this simulation, including total irrigation water volume, actual yield, and end-of-season root zone salinity, is considered the baseline target vector that best approximates the actual physical processes.

[0202] Perform performance comparison and error assessment. Compare the predicted target vector with the baseline target vector item by item, and calculate the prediction error of the surrogate model at the optimal solution point. For example, calculate the relative error of each target:

[0203] Error yield =|Y a_true -Y a_proxy | / Y a_true ;

[0204] Error salt =|EC true -EC proxy | / EC true ;

[0205] Among them, Error yield Y represents the relative error in production forecasting. a_true As the baseline simulated output, Y a_proxy Predicting output for the surrogate model. Error salt EC represents the relative error in the quarter-end salinity forecast. true As a baseline for simulated salinity, EC proxy To predict salinity using a surrogate model.

[0206] Decisions are made based on the error assessment results. The calculated relative errors are compared with a preset acceptable error threshold, such as five percent.

[0207] If the relative errors of all targets are less than the acceptable error threshold, it indicates that the surrogate model's predictions in that region are accurate and reliable. Therefore, the scheme to be validated passes validation and is confirmed as the final optimal irrigation scheme, which can be output for practical application.

[0208] If the relative error of at least one target exceeds the acceptable error threshold, it indicates that the surrogate model's accuracy in that critical region is insufficient, and its recommended solution may not be the true optimal solution. In this case, a callback mechanism is triggered. Specifically, the callback mechanism involves adding the solution to be validated and its corresponding baseline target vector as new, high-value training data points to the training dataset. The surrogate model is then retrained or updated using the augmented dataset. The updated surrogate model, having learned accurate information in its previously poorly performing regions, will see improvements in both its overall accuracy and local accuracy near the optimal region.

[0209] In an optional, more refined implementation, after executing the callback mechanism and updating the surrogate model, the multi-objective optimization algorithm can be re-executed. Since the surrogate model has been updated, this optimization may generate a revised, higher-quality Pareto optimal solution set, thereby yielding a revised, more reliable optimal compromise. This solution can then re-enter the verification process of this embodiment, forming a closed-loop, iterative optimization process until a final, verified solution is found.

[0210] This embodiment enhances the robustness and reliability of the invention in practical applications by introducing a closed-loop process of prediction-verification-callback, ensuring that the final output irrigation scheme is not only theoretically optimal based on the surrogate model, but also a reliable and feasible practical scheme confirmed by a high-precision physical model.

[0211] It should be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately.

Claims

1. An irrigation optimization method based on crop models, characterized in that, include: Acquire input data for irrigation optimization, including meteorological, soil, and crop data for the target area; Based on the input data, the crop growth model is run to simulate crop growth and soil water and salt dynamics under the irrigation scheme to be evaluated, and the corresponding simulation results are obtained. During the operation of the crop growth model, the root water uptake stress is dynamically calculated according to the crop growth stage. Based on the simulation results and according to multiple preset optimization objectives, the optimal irrigation scheme is obtained by iteratively solving the irrigation scheme through a multi-objective optimization algorithm. Output the optimal irrigation plan.

2. The method according to claim 1, characterized in that, Dynamic calculations of root water uptake stress based on crop growth stages are performed, specifically including: During the operation of the crop growth model, crop growth stage indicators corresponding to the crop growth stages are calculated based on the simulated crop development process. By substituting crop growth stage indicators as independent variables into one or more pre-defined continuous functions used to characterize stress response, the root water absorption stress under the current growth stage can be determined.

3. The method according to claim 2, characterized in that, Determining root water uptake stress at the current growth stage is achieved through a split-stress framework, which includes: Based on the soil moisture state simulated by the crop growth model, the water stress reduction coefficient is calculated independently. Based on soil salinity status and crop growth stage indicators simulated by crop growth models, the salt stress reduction coefficient is calculated independently. The calculation process of the salt stress reduction coefficient reflects the influence of growth stage. The reduction coefficients of water stress and salt stress are combined to obtain the final root water uptake stress.

4. The method according to claim 3, characterized in that, Based on soil salinity status and crop growth stage indicators simulated by crop growth models, the salinity stress reduction coefficient is calculated independently, specifically including: The crop growth stage index is used as an independent variable and input into the first preset continuous function to calculate the salt tolerance threshold that changes dynamically with the growth stage. The dynamically changing salt tolerance threshold and soil salinity status are used as inputs and substituted into the salt stress response function to calculate the salt stress reduction coefficient.

5. The method according to claim 4, characterized in that, The first preset continuous function is a polynomial function or a radial basis function, which maps crop growth stage indicators to dynamically changing salt tolerance thresholds through a set of pre-configured coefficients or parameters.

6. The method according to claim 3, characterized in that, The split-stress framework also includes: The crop growth stage index is used as an independent variable and input into the second preset continuous function to calculate the interaction intensity coefficient that changes dynamically with the growth stage. Based on the interaction intensity coefficient, water stress reduction coefficient, and salt stress reduction coefficient, the salt-drought interaction correction term was calculated. Specifically, the water stress reduction coefficient and the salt stress reduction coefficient are combined. The water stress reduction coefficient, the salt stress reduction coefficient, and the salt-drought interaction correction term are combined to obtain the final root water uptake stress.

7. The method according to claim 2, characterized in that, Determining root water uptake stress at the current growth stage is achieved through a unified effective water potential framework, which includes: The osmotic potential was calculated based on the soil salinity state simulated by the crop growth model. Based on crop growth stage indicators, the penetration sensitivity coefficient that dynamically changes with the growth stage is determined through a preset continuous function. After weighting the osmotic potential using the permeability sensitivity coefficient, the weighted osmotic potential is combined with the soil matrix potential simulated by the crop growth model to construct a unified effective water potential. Based on the unified effective water potential, the final root water absorption stress was calculated.

8. The method according to claim 7, characterized in that, The penetration sensitivity coefficient, which dynamically changes with the reproductive stage, is determined by a pre-defined continuous function, specifically including: By using crop growth stage indicators as independent variables and inputting them into the third preset continuous function, the penetration sensitivity coefficient that dynamically changes with the growth stage can be calculated. The permeability potential is weighted using the permeability sensitivity coefficient to obtain the weighted permeability potential; By merging the weighted osmotic potential with the soil matrix potential, a unified effective water potential is obtained.

9. The method according to claim 8, characterized in that, Based on the unified effective water potential, the final root water uptake stress was calculated as follows: By substituting the uniform effective water potential as the sole independent variable into a single, pre-defined stress response function, the final root water absorption stress is calculated.

10. The method according to claim 2, characterized in that, Determining root water uptake stress at the current growth stage is performed within the numerical iterative loop of the crop growth model used to solve the hydrodynamic equations, specifically including: In each iteration of the numerical iteration loop, the soil water potential and soil salinity status under the current iteration are obtained; Based on the soil water potential, soil salinity, and crop growth stage indicators under the current iteration, the root water uptake stress under the current iteration is calculated. The root water absorption stress of the current iteration is used as the sink term of the water flow equation to solve the updated soil water potential in the same iteration until the numerical iteration cycle converges.

Images