A method and system for optimizing the scheduling of integrated water and fertilizer irrigation

Abstract

This invention belongs to the field of optimization technology, specifically relating to a method and system for optimizing integrated water and fertilizer irrigation scheduling. The method includes: acquiring data from the water and fertilizer system; extracting the cumulative effect and rate of change values ​​of multi-scale time windows to generate composite time-series features; constructing an irrigation demand prediction model; quantifying the coupling degree between candidate features and path features using a feature adjacency knowledge graph and reducing the sampling probability accordingly; generating an initial subset of candidate features; when the maximum splitting evaluation index is lower than a preset lower limit; incorporating the strongly correlated adjacent features of the optimal feature according to the feature adjacency knowledge graph; updating the connection weights of the feature adjacency knowledge graph; and outputting the expected value vector and covariance matrix; inputting the data into the irrigation demand prediction model to obtain the expected value vector and covariance matrix at the current time; and solving for the optimal irrigation duration, fertilization timing, and ratio command. This invention takes into account prediction uncertainty, improving prediction accuracy and the stability of control command output.

Description

Technical Field

[0001] This invention relates to the field of optimization technology. More specifically, this invention relates to a method and system for optimizing the scheduling of integrated water and fertilizer irrigation. Background Technology

[0002] Integrated water and fertilizer management technology combines irrigation and fertilization, applying it to the high-yield and high-quality cultivation process of modern agriculture. It plays a core role in precise water and fertilizer management based on environmental and crop conditions. In actual agricultural production, there are nonlinear interactions between environmental factors such as temperature and humidity, light, soil conditions, and crop growth processes. Furthermore, the conduction and absorption of water and fertilizer in the soil and crops often exhibit time-cumulative effects and dynamic lags. Early water and fertilizer scheduling relied primarily on simple threshold judgments based on single-dimensional historical data or conventional shallow rule models. These conventional data-driven models lack the comprehensive ability to represent the cumulative trends and transient changes of environmental factors across different periods when facing complex, multi-scale dynamic environments. This makes it difficult for the system to accurately identify the true water and fertilizer needs of crops under complex environmental alternations, thus easily leading to lagging or biased irrigation and fertilization strategies.

[0003] To overcome the shortcomings of conventional models in dynamically representing complex environmental factors, improved schemes using high-dimensional data-driven machine learning models for water and fertilizer demand prediction and scheduling have gradually evolved in existing technologies. This improved scheme typically inputs collected multidimensional environmental data and massive historical features directly into prediction models such as decision trees for training, leveraging the model's data fitting capabilities to uncover the impact of various parameters on water and fertilizer demand. In practical applications, the control system, based on the single deterministic target value of water and fertilizer demand output by this prediction model, directly maps and generates corresponding low-level control commands such as irrigation duration, fertilization sequence, and fertilizer ratio, thereby achieving automated scheduling and control of the integrated water and fertilizer system.

[0004] While the aforementioned improvements enhance multidimensional data fitting, they still have shortcomings in the depth of feature interaction mining and the security of command scheduling. High-dimensional time-series data collected by actual agricultural IoT contains valuable high-order correlation features, and agricultural physical systems are inherently subject to random environmental fluctuations and predictive uncertainties. However, existing high-dimensional prediction models mostly employ isolated feature evaluation mechanisms during training and optimization, with the control end relying solely on a single, definite prediction target value to generate scheduling commands. This mechanism, on the one hand, ignores high-order interactions and coupling relationships along complex feature paths, causing the model to easily fall into local optima when facing sudden environmental changes or complex growth cycles due to redundant feature interference or omission of key features, resulting in a significant decrease in the accuracy of front-end demand prediction. On the other hand, due to the lack of effective quantification methods for the uncertainty of prediction results, the control commands cannot establish a risk avoidance mechanism for prediction deviations during the generation process. When actual environmental conditions deviate from the ideal distribution assumptions of the prediction model, scheduling commands generated based on absolutely deterministic targets will experience severe supply-demand mismatches, leading to frequent excessive or insufficient water and fertilizer supply actions by the system, ultimately causing unnecessary waste of water and fertilizer resources and even substantial obstruction of crop growth. Summary of the Invention

[0005] To address the technical challenges of overcoming the local optima of models caused by the omission of high-dimensional correlation features in the water and fertilizer demand forecasting process, and to quantify the forecast uncertainty at the scheduling decision-making stage to avoid the risk of supply and demand mismatch caused by environmental fluctuations, this invention provides a water and fertilizer integrated irrigation scheduling optimization method and system.

[0006] In a first aspect, the present invention provides a method for optimizing the scheduling of integrated water and fertilizer irrigation, comprising: S1, acquiring historical environmental data, crop growth data, and irrigation data of the integrated water and fertilizer system, extracting cumulative effect values ​​and rate of change values ​​under multi-scale time windows, and generating composite time-series features; S2, constructing an irrigation demand prediction model, performing initialization of a feature adjacency knowledge graph representing high-order interaction relationships during training, tracking existing splitting paths at decision tree splitting nodes, quantifying the coupling degree between candidate features and path features and decaying the sampling probability accordingly, generating an initial candidate feature subset, and when the maximum splitting evaluation index based on the initial candidate feature subset is lower than a preset lower limit, according to the feature adjacency knowledge... The graph incorporates the strongly correlated adjacent features of the best features in the initial candidate feature subset to obtain an expanded candidate feature subset and re-finds the optimal split point. The connection weights of the feature adjacency knowledge graph are updated according to the changes in the evaluation index until the model is generated. The leaf nodes of the irrigation demand prediction model output the expected value vector and covariance matrix. S3: Input the real-time collected data into the irrigation demand prediction model to obtain the expected value vector and covariance matrix at the current moment. Establish a system response mapping model. With the expected value vector as the optimization center, use the covariance matrix to quantify uncertainty. Combine the system response mapping model to construct a risk avoidance optimization problem and solve the optimal irrigation duration, fertilization timing, and fertilizer ratio combination instructions.

[0007] This invention extracts composite time-series features under multi-scale time windows, and uses knowledge graphs to quantify feature coupling and decay sampling probabilities when constructing a prediction model. It dynamically expands candidate subsets and updates graph weights by combining strongly correlated features, and the model's leaf nodes output expectations and covariance. A system response mapping model is established, and the covariance is used to quantify uncertainty and construct a risk avoidance optimization problem to solve for the optimal command. This invention can identify high-order correlations between multi-dimensional time-series features, avoid the model training from getting stuck in local optima, and fully quantify the uncertainty risk in the prediction process. This improves the stability and accuracy of water and fertilizer control commands and avoids resource waste and crop growth hindrance.

[0008] Preferably, the method for obtaining the composite time series features includes: setting time windows of lengths of 1 day, 7 days, and 30 days, and truncating the continuous historical time series observation sequence into three benchmark observation stage windows with different periods; for each benchmark observation stage window, for flux-type variables, summing the historical time series observation data within the window, and for state-type variables, calculating the average value of the historical time series observation data within the window, and unifying the summation result with the average value as the cumulative effect value of the corresponding time scale; subtracting the same type of time series observation data from the time series observation data of the 1st day before the current time from the time series observation data of the 2nd day before the current time to obtain the difference, dividing the difference by the corrected benchmark value to obtain the corresponding data change rate value, where the corrected benchmark value is the sum of the same type of time series observation data of the 2nd day before the current time and a preset small positive number of the same dimension; and concatenating all extracted cumulative effect values ​​and change rate values ​​to generate composite time series features.

[0009] This invention, by dividing windows into different time scales and calculating the accumulation and rate of change of flux-type and state-type variables respectively for variable types, can capture the dynamic cumulative effect and short-term fluctuation characteristics of the water and fertilizer system under different cycles in a more granular manner, thereby enhancing the model's ability to express the evolution of historical environment and crop state.

[0010] Preferably, the initialization of the feature adjacency knowledge graph representing high-order interaction relationships during training includes: using the generated composite temporal features as node entities to construct an initial feature adjacency knowledge graph; calculating the Pearson correlation coefficient between any two feature nodes using a calculation formula, which is the product of the covariance of the two feature sequences and the standard deviation of the two feature sequences; determining whether the absolute value of the Pearson correlation coefficient between any two feature nodes is greater than a preset correlation threshold, and if so, establishing an undirected edge between the two feature nodes and using the undirected edge as the initial connection relationship of the feature adjacency knowledge graph.

[0011] This invention utilizes the Pearson correlation coefficient to construct an initial feature undirected graph, which can intuitively and cost-effectively map the initial linear correlation strength between time-series features, providing a reliable prior topological structure for subsequent exploration of complex nonlinear interaction effects.

[0012] Preferably, the step of quantifying the coupling degree between the candidate features and the path features and thereby attenuating the sampling probability to generate an initial candidate feature subset includes: backtracking upwards from the current split node to the root node to obtain all used split features on the split path; querying the non-negative connection weights between each candidate feature and the used split features in the feature adjacency knowledge graph, summing the non-negative connection weights between a single candidate feature and all used split features to obtain the coupling degree of the candidate feature; normalizing the coupling degree of each candidate feature to a closed interval of 0 to 1, calculating 1 minus the normalized coupling degree to obtain a probability attenuation coefficient, multiplying the initial sampling probability of each candidate feature by the probability attenuation coefficient to obtain the adjusted sampling probability, performing sum-normalization on the adjusted sampling probability, and performing random sampling based on the normalized adjusted sampling probability to generate an initial candidate feature subset.

[0013] This invention calculates the coupling degree between candidate features and existing splitting path features by graph backtracking and converts it into a probability decay coefficient to adjust the sampling probability. This reduces the probability of highly redundant features being selected and forces the decision tree to explore a feature space with novel regression contributions during splitting.

[0014] Preferably, the step of incorporating the strongly correlated adjacent features of the optimal feature into the feature adjacency knowledge graph to obtain an expanded candidate feature subset includes: when the maximum splitting evaluation index of the initial candidate feature subset is lower than a preset splitting quality lower limit, locking the feature with the highest evaluation index in the initial candidate feature subset as the core feature; querying all adjacent features that have a connection relationship with the core feature in the feature adjacency knowledge graph; sorting all adjacent features in descending order of connection weight, and extracting the top 3 adjacent features as strongly correlated adjacent features; adding the strongly correlated adjacent features to the initial candidate feature subset and performing a deduplication operation to obtain the expanded candidate feature subset.

[0015] When the quality of the current candidate feature set is low, this invention directly utilizes the high-weight adjacency features of the core features in the graph for directional expansion, avoiding invalid blind search caused by random sampling and accelerating the convergence of high-quality split nodes.

[0016] Preferably, updating the connection weights of the feature adjacency knowledge graph based on changes in evaluation metrics includes: recording a first maximum split evaluation metric obtained based on an initial subset of candidate features, and a second maximum split evaluation metric obtained by re-finding the optimal split point based on an expanded subset of candidate features; calculating the difference between the second maximum split evaluation metric and the first maximum split evaluation metric; if the difference is greater than a preset enhancement threshold, then in the feature adjacency knowledge graph, increasing the connection weight between the core feature and the included strongly associated adjacency features by 1; if the difference is less than or equal to the preset enhancement threshold, then in the feature adjacency knowledge graph, decreasing the connection weight between the core feature and the included strongly associated adjacency features by 0.1, and taking the maximum value between the reduced connection weight and 0, so that the connection weight remains non-negative.

[0017] This invention dynamically adjusts the graph connection weights based on the difference in the gain of the split evaluation index before and after the introduction of new features, which can achieve adaptive evolution of the graph topology, gradually filter out false associations and accurately retain the real high-order feature interaction rules.

[0018] Preferably, the leaf node output expected value vector and covariance matrix of the irrigation demand prediction model include: calculating the mean of the target water and nutrient concentrations as the expected value vector based on the training sample set falling into the leaf node, and calculating the corresponding sample covariance matrix; when the actual number of samples falling into the leaf node is less than the dimension of the target state variable plus 1, backtracking to the parent node sample set to calculate the sample covariance matrix, or using the diagonal elements of the global training sample covariance matrix to form a diagonal covariance matrix as a degenerate replacement; superimposing a regularization term on the main diagonal of the obtained sample covariance matrix or diagonal covariance matrix to generate the final covariance matrix as the output result of the leaf node.

[0019] When leaf node samples are sparse, this invention prevents the covariance matrix from losing its positive definiteness by backing up the parent node or using a global diagonal matrix as a degenerate alternative, and by superimposing a regularization term on the diagonal. This ensures the numerical stability of the subsequent uncertainty quantification and optimization solution process based on covariance.

[0020] Preferably, the step of constructing the risk avoidance optimization problem by combining the system response mapping model includes: extracting the expected values ​​of target soil moisture content, target nitrogen concentration, target phosphorus concentration, and target potassium concentration from the expected value vector output by the irrigation demand prediction model, and using them as optimization centers; extracting the variance elements on the main diagonal from the covariance matrix, and performing square root operations on the variance elements to obtain the standard deviation of the corresponding target state variables; calculating the risk safety margin based on the standard deviation of each target state variable, wherein the risk safety margin is the product of the corresponding standard deviation and a preset confidence coefficient; shrinking the preset safety constraint boundary of each target state variable based on the risk safety margin to obtain a conservative constraint boundary, and truncating the lower boundary of the shrunken boundary by taking the maximum value of the result with 0; and using the conservative constraint boundary as the boundary of the constraint conditions in the risk avoidance optimization problem.

[0021] This invention analyzes the covariance matrix into standard deviation and transforms it into a risk safety margin to shrink the constraint boundary inward. It directly maps the probability distribution predicted by the model into a conservative feasible region. When the confidence of the model prediction is low, it automatically limits the amplitude of the control action to prevent over-irrigation or fertilizer damage.

[0022] Preferably, the step of solving for the optimal combination of irrigation duration, fertilization timing, and fertilizer ratio includes: constructing an optimization objective function based on risk avoidance principles; the optimization objective function includes a first part and a second part; the first part is a weighted evaluation term of the difference between the water state and nutrient concentration state predicted by the system response mapping model and the expected value vector; the second part is a risk penalty evaluation term based on the deviation of the conservative constraint boundary and the standard deviation of the corresponding target state variable; iteratively searching for the minimum value of the optimization objective function; extracting the array of independent variables returned by the optimization solver after convergence; and reversing the array of independent variables to obtain the optimal combination of irrigation duration, fertilization timing, and fertilizer ratio.

[0023] This invention introduces a weighted evaluation of the expected difference in state and a risk penalty term based on the deviation of the constraint boundary into the objective function. This enables the optimization solver to actively avoid high-risk parameter regions that are prone to causing harm while pursuing the approximation of the target state.

[0024] Secondly, the present invention provides an integrated water and fertilizer irrigation scheduling optimization system, including a processor and a memory, wherein the memory stores computer program instructions, and when the computer program instructions are executed by the processor, the aforementioned integrated water and fertilizer irrigation scheduling optimization method is implemented.

[0025] By adopting the above technical solution, a computer program is generated from the above-mentioned water and fertilizer integrated irrigation scheduling optimization method and stored in the memory so that it can be loaded and executed by the processor. In this way, a terminal device can be made based on the memory and the processor for convenient use.

[0026] The technical solution of the present invention has the following beneficial technical effects:

[0027] This invention acquires historical environmental, crop growth, and irrigation data of an integrated water and fertilizer system, extracts cumulative effect values ​​and rate of change values ​​under multi-scale time windows, and generates composite time-series features. The above implementation method extracts cumulative effect values ​​and rate of change values ​​of historical environmental, crop growth, and irrigation data through multi-scale time windows, thereby generating composite time-series features. This improves upon the problem of not being able to generate composite time-series features due to the lack of extraction of cumulative effect values ​​and rate of change values ​​under multi-scale time windows, providing composite time-series feature data as input for subsequent model training.

[0028] This invention constructs an irrigation demand prediction model. During training, it initializes a feature adjacency knowledge graph representing high-order interaction relationships among features. At decision tree split nodes, it tracks existing splitting paths, calculates the coupling degree between candidate features and path features, and decays the sampling probability accordingly to generate an initial subset of candidate features. When the maximum splitting evaluation index based on the initial subset of candidate features falls below a preset lower limit, it incorporates strongly correlated adjacent features of the optimal feature from the initial subset of candidate features into the feature adjacency knowledge graph, resulting in an expanded subset of candidate features. It then searches for the optimal splitting point again, updating the connection weights of the feature adjacency knowledge graph based on changes in the evaluation index until the model is generated. The leaf node outputs of the irrigation demand prediction model represent the expected value vector and covariance matrix of the target water and nutrient concentrations. This implementation utilizes a feature adjacency knowledge graph to represent high-order interaction relationships among features and calculates coupling degree. Based on the maximum splitting evaluation index, strongly correlated adjacent features are incorporated into the initial subset of candidate features. The optimal splitting point is found and connection weights are updated based on changes in the evaluation index, thereby generating an irrigation demand prediction model with high-order interaction relationships among features, and obtaining the expected value vector and covariance matrix of the target water and nutrient concentrations.

[0029] This invention inputs real-time collected data into an irrigation demand prediction model to obtain the expected value vector and covariance matrix at the current moment. It then establishes a system response mapping model that maps irrigation duration, fertilization timing, and fertilizer ratio combination commands to water and nutrient concentration states. This invention uses the expected value vector as the optimization center, utilizes the covariance matrix to represent uncertainty, and combines it with the system response mapping model to construct a risk-avoidance optimization problem. Under preset constraints, it solves for the optimal combination of irrigation duration, fertilization timing, and fertilizer ratio commands. The above implementation combines the expected value vector and covariance matrix with the system response mapping model, using the covariance matrix to represent uncertainty and construct a risk-avoidance optimization problem. This allows for the solution of the optimal combination of irrigation duration, fertilization timing, and fertilizer ratio commands under preset constraints, obtaining the optimal combination commands to complete the optimization of integrated water and fertilizer irrigation scheduling. Attached Figure Description

[0030] Figure 1This is a flowchart of a water and fertilizer integrated irrigation scheduling optimization method; Figure 2 This is a time-series broken line diagram of soil moisture content. Detailed Implementation

[0031] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0032] This invention discloses a method for optimizing the scheduling of integrated water and fertilizer irrigation, referring to... Figure 1 This includes the following steps: S1, extract data to generate composite time-series features. Specifically, acquire historical environmental, crop growth, and irrigation data of the fertigation system, extract cumulative effect values ​​and rate of change values ​​under multi-scale time windows, and generate composite time-series features.

[0033] Data on air temperature and humidity, light intensity, and soil electrical conductivity were collected by soil temperature and humidity sensors deployed in farmland and by weather stations. A connection was established with a relational database, and historical data was retrieved using a structured query language. The data was loaded into an in-memory data frame. For missing observations caused by sensor communication delays, linear interpolation was used to complete the data. Time windows of 24h, 168h, and 720h were set, corresponding to 1 day, 7 days, and 30 days respectively. Cumulative rainfall and cumulative photosynthetic radiation were calculated as cumulative effect values ​​within the multi-scale time windows. For state variables such as soil moisture and soil electrical conductivity, the mean of the corresponding window was calculated as the cumulative effect value. The relative change rates of soil moisture and electrical conductivity were obtained by dividing the difference between statistical values ​​from adjacent days by the sum of preset small positive numbers consistent with the units of the previous day's statistical values ​​and corresponding variables. The original environmental data, crop growth data, and the newly extracted cumulative effect and change rate values ​​were concatenated along the column dimension, and feature standardization processing with mean removal and variance normalization was performed to generate a multi-dimensional composite time-series feature matrix. In this feature matrix, each row corresponds to a sampling time or a training sample, and each column corresponds to an original variable, cumulative effect value, or rate of change value. Subsequent model training uses the feature matrix concatenated by this column dimension as input.

[0034] In an optional embodiment, the cumulative effect value and rate of change value under multi-scale time windows are extracted to generate composite time series features, including: By setting time windows of 1 day, 7 days, and 30 days, the continuous historical time series observation sequence is truncated into three benchmark observation phase windows with different periods. For each baseline observation window, the data is processed according to the variable type of the observed data. For flux-type variables, the historical time series observation data within the window are summed. For state-type variables, the average value of the historical time series observation data within the window is calculated. The summation result and the mean value are combined as the cumulative effect value of the corresponding time scale. Extract the time series observation data from the day before the current time, subtract the same time series observation data from the day before the current time to obtain the difference, divide the difference by the correction benchmark value, which is the sum of the same time series observation data from the day before the current time and a preset small positive number of the same dimension, to obtain the corresponding data change rate value. All extracted cumulative effect values ​​and rate of change values ​​are concatenated to generate composite time-series features.

[0035] Historical observation sequences were truncated using hourly sampling (e.g., once per hour, 24 data points per day). Time windows of 1 day, 7 days, and 30 days corresponded to observation data segments from the past 24, 168, and 720 time steps, respectively. During data processing, flux-type variables, such as cumulative irrigation, natural precipitation, and total fertilizer application, were distinguished from state-type variables, such as soil temperature, relative soil humidity, and air temperature and humidity. For flux-type variables, a summation operation was used to obtain the total input; for example, assuming the hourly rainfall sequence for the past 24 hours is [0, 2, 0, 0, 3, 0, 0, 0, 0, 0, 5, 6, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5] mm, then the cumulative effect value on a 1-day scale is the sequence summation of 25 mm. For state-type variables, mean calculation is used to smooth fluctuations; for example, if the soil moisture sequence over the past 24 hours is [20%, 21%, 20%, 19%, 20.5%, 20.5%, 20.2%, 19.8%, 19.5%, 19%, 21.5%, 23%, 24%, 24%, 23.5%, 22.8%, 22%, 21%, 20%, 19%, 18.2%, 17.5%, 17%, 19%], then the mean of the sequence, 20.5%, is calculated as the cumulative effect value.

[0036] When extracting the rate of change of data, the relative difference is calculated using the daily average values ​​of two adjacent days. To avoid the denominator being zero, a preset small positive number ε is added to the denominator when the similar time-series observation data of the second day before the current time is zero or close to zero. For example, if the daily average soil moisture of the first day before the current time is 21%, and the daily average soil moisture of the second day before the current time is 20%, the calculated relative rate of change is 5%, which is used to represent the short-term trend. The cumulative effect values ​​of each variable at the three time scales, such as the mean and sum, as well as the extracted rate of change values, are converted into numerical vectors, and the column dimensions are concatenated sequentially according to the variable order to form, for example, a composite time-series feature vector with a length of 50 dimensions. In this embodiment, the measurement accuracy of the soil moisture sensor is ±1%, and the preset small positive number ε is 0.5%.

[0037] Reference Figure 2 The time-series line graph of soil moisture content reflects the trend of environmental data fluctuations, demonstrating the dynamic accumulation and consumption of moisture over continuous time steps. Based on the fluctuation patterns of moisture data, the technical effectiveness of extracting multi-dimensional time-series features in accurately characterizing historical environmental states is proven.

[0038] S2, Construct an irrigation demand prediction model and output the expected value and covariance. Specifically, during training, the irrigation demand prediction model is constructed. Initialization is performed on the feature adjacency knowledge graph representing high-order interactions between features. Existing splitting paths are tracked at decision tree splitting nodes. The coupling degree between candidate features and path features is quantified, and the sampling probability is decayed accordingly to generate an initial candidate feature subset. When the maximum splitting evaluation index based on the initial candidate feature subset is lower than a preset lower limit, the strongly correlated adjacent features of the optimal feature in the initial candidate feature subset are incorporated according to the feature adjacency knowledge graph to obtain an expanded candidate feature subset. The optimal splitting point is then re-searched. The connection weights of the feature adjacency knowledge graph are updated according to the changes in the evaluation index until the model is generated. The leaf node outputs of the irrigation demand prediction model represent the expected value vector and covariance matrix of the target water and nutrient concentrations.

[0039] A multivariate Gaussian regression tree model is used to construct an irrigation demand prediction model. During the model training initialization phase, an undirected graph is constructed as a feature adjacency knowledge graph. Pearson correlation coefficients are calculated for each pair of composite time-series features. When the absolute value of the Pearson correlation coefficient is greater than a preset correlation threshold of 0.6, an undirected edge is established between the corresponding two feature nodes, and the absolute value of the Pearson correlation coefficient is used as the initial non-negative connection weight of this undirected edge. At the split nodes of the regression tree growth, all feature indices traversed from the root node to the current node are recorded to form an existing split path. The non-negative connection weights between each candidate feature and each feature in the path feature set are calculated, and the coupling degree is obtained by summing the non-negative connection weights. After normalizing the coupling degree to a closed interval of 0 to 1, the probability decay coefficient is calculated by subtracting the normalized coupling degree from 1, and the sampling probability is updated using the probability decay coefficient. The correlation threshold is determined through a grid search on the validation set, with a search range of [0.4, 0.8] and a step size of 0.05, using the minimization of prediction error on the validation set as the selection criterion. For each candidate relevance threshold, an initial feature adjacency knowledge graph is constructed, a decision tree model is trained, and the prediction error on the validation set is recorded. The relevance threshold that minimizes the validation error is selected as the final value.

[0040] An initial subset of candidate features is generated by extracting feature indices without replacement based on the decayed probability distribution. The reduction in normalized weighted variance of the features within this subset is calculated as the maximum split evaluation index. Before calculating the variance reduction, soil moisture and nitrogen, phosphorus, and potassium concentrations are dimensionless according to the training set standard deviation of each target state variable, a preset safety margin width, or a preset scale factor. Then, the variance reduction of each target state variable is weighted and summed. If this index is lower than the set split quality lower limit threshold of 0.05, the top three features with the highest connection weights to the optimal feature with the largest normalized weighted variance reduction in the current subset are queried from the feature adjacency knowledge graph and selected as strongly correlated adjacency features. The split quality lower limit threshold is determined based on the training set size and the dimension of the target variables. When the sample size is sufficient, a lower threshold can be set to increase the split depth; when the sample size is small, the threshold should be increased to prevent overfitting.

[0041] Strongly correlated adjacent features are merged to generate an expanded candidate feature subset. This expanded subset is then re-traversed to calculate the reduction in normalized weighted variance and find the optimal split point. After splitting, the gain magnitude of the reduction in normalized weighted variance of the expanded candidate feature subset compared to the initial candidate feature subset is calculated. The connection weights of the corresponding edges in the feature adjacency knowledge graph are updated based on this gain magnitude. This recursive splitting continues until the maximum tree depth limit or the leaf node sample number limit is met. High-order feature interaction relationships refer to the combined interaction relationships formed by the cumulative path of multiple feature nodes through pairwise connection weights under the condition of multi-level splitting paths in the decision tree. The initial graph uses pairwise correlations to build edges. During training, the interaction influence under multi-feature conditions is expressed by the cumulative coupling degree of multiple connections between the path feature set and the candidate features.

[0042] Upon reaching a leaf node, based on the training sample set falling into that leaf node, the mean values ​​of the target moisture and nutrient concentrations are calculated as the expected value vector. The corresponding sample covariance matrix is ​​then calculated and persistently stored as the output of that leaf node. To avoid covariance matrix degradation due to insufficient leaf node samples, a minimum sample size for each leaf node is preset to be no less than the dimension of the target state variable plus one. When the actual number of samples falling into a leaf node is insufficient, further splitting is stopped, and the calculation of the covariance matrix is ​​back to the parent node's sample set. Alternatively, a diagonal covariance matrix is ​​constructed from the diagonal elements of the global training sample covariance matrix as a degradation substitute. After calculating the sample covariance matrix, a regularization term λI is superimposed on its main diagonal, where λ is a preset small positive number to ensure that the output covariance matrix is ​​a positive definite or approximately positive definite matrix, facilitating subsequent risk avoidance optimization calculations. In this embodiment, the regularization term coefficient λ is preset to 1×10⁻¹⁰. -5 Specifically, the mean variance of each target variable in the training set is calculated, and λ is taken as 0.1% of this mean. In this embodiment, the mean is 1×10⁻⁶. -2 Therefore, λ = 1 × 10 -5 When the condition number exceeds 1×10 4 At that time, λ dynamically increases to 10 times its original value, until the condition number is less than 1 × 10. 4 .

[0043] In an optional embodiment, initializing a feature adjacency knowledge graph representing high-order interaction relationships of features includes: The generated composite temporal features are used as node entities to construct an initial feature adjacency knowledge graph. The Pearson correlation coefficient between any two feature nodes is calculated using the formula: the covariance of the two feature sequences divided by the product of the standard deviations of the two feature sequences. Determine whether the absolute value of the Pearson correlation coefficient between any two feature nodes is greater than a preset correlation threshold. If so, establish an undirected edge between the two feature nodes and use the undirected edge as the initial connection relationship of the feature adjacency knowledge graph.

[0044] The graph construction process maps the extracted composite temporal feature set to a node entity set of the knowledge graph. Assuming the preceding steps generated 50 composite features, an undirected graph with 50 vertices is constructed. For any two feature nodes, the Pearson correlation coefficient is calculated using historical training data. This coefficient is the product of the covariances of the two variables divided by their respective standard deviations, representing the linear relationship between the two feature sequences.

[0045] When building the topology, all 1225 feature combinations are traversed. For example, assuming features Cumulative daily rainfall and its characteristics The Pearson correlation coefficient calculated from the 1-day soil moisture change rate was 0.78. Since the absolute value of 0.78 is greater than the set correlation threshold of 0.6, a correlation was determined between the two, and this was established at the node of the graph. and Instantiate an undirected edge between them and initialize the absolute value of the Pearson correlation coefficient, 0.78, as the non-negative connection weight of the undirected edge. Conversely, if the feature 30-day average temperature and If the absolute value of the correlation coefficient of the 1-day phosphorus change rate is 0.25 and not greater than 0.6, then the two nodes remain disconnected. The generated map can be stored using a 50×50 sparse adjacency matrix, where non-zero elements indicate initial feature associations, improving the guidance of feature interaction modeling. For two feature nodes without an established edge, their connection weight is recorded as 0 for subsequent coupling degree summation calculation.

[0046] In an optional embodiment, existing splitting paths are traced at the decision tree splitting nodes, the coupling degree between candidate features and path features is quantified, and the sampling probability is decayed accordingly to generate an initial subset of candidate features, including: At the current split node, backtrack upwards to the root node to obtain all used split features along the split path; In the feature adjacency knowledge graph, query the non-negative connection weights between each candidate feature and the used split features in the candidate feature set, and sum the non-negative connection weights between a single candidate feature and all used split features to obtain the coupling degree of the candidate feature. The coupling degree of each candidate feature is normalized to a closed interval of 0 to 1. The probability decay coefficient is obtained by subtracting the normalized coupling degree from 1. The initial sampling probability of each candidate feature is multiplied by the probability decay coefficient to obtain a non-negative adjusted sampling probability. The adjusted sampling probabilities are then summed and normalized. Random sampling is performed based on the normalized adjusted sampling probabilities to generate an initial subset of candidate features. When there are no used split features on the split path, or when the coupling degree of each candidate feature is the same, making normalization impossible, the adjusted sampling probability of each candidate feature is set to a uniform distribution.

[0047] When constructing an irrigation demand prediction model based on a multivariate Gaussian regression tree or an ensemble regression tree model composed of multiple multivariate Gaussian regression trees, it is assumed that the current optimization process is at a split node with a depth of 3. Tracing back up the tree branches to the root node, the feature sets used as splitting conditions along the way are extracted. For example, the used feature set includes features... and For all unused candidate features at the current node (e.g., 20 candidate features), index and query them one by one in the pre-constructed feature adjacency knowledge graph. After several rounds of splitting feedback updates, some connection weights may be greater than the absolute value of the initial Pearson correlation coefficient. (The remaining text appears to be incomplete and requires further context.) For example, suppose the query finds the map... and The connection weight is 1.2. and If the connection weight is 0.8, then The total coupling degree of the existing path features is calculated to be 2.

[0048] To map the coupling degree to a probability decay coefficient, the coupling degree set calculated from all candidate features is extracted, and Min-Max normalization is performed to map the coupling degree set to the closed interval [0,1]. If the features... The normalized coupling degree is 0.6, so the calculated difference, i.e., the attenuation multiplier, is 0.4. If we assume that the initial sampling probabilities of all candidate features are equal before adjustment, i.e., the initial probability is 0.05, then... The adjusted sampling probability becomes 0.02. After normalizing the adjusted probabilities of all candidate features to a sum of 1, a roulette wheel random sampling mechanism with no replacement is executed to select a fixed number of features to form an initial subset of candidate features. This step reduces the probability of features redundant with already used splitting features being selected, allowing the model to explore new features that contribute to regression splitting, thus improving the diversity and overfitting resistance of the tree model.

[0049] In an optional embodiment, the strongly correlated adjacent features of the best feature in the initial candidate feature subset are incorporated into the feature adjacency knowledge graph to obtain an expanded candidate feature subset, including: When the maximum splitting evaluation index of the initial candidate feature subset is lower than the preset splitting quality lower limit, the feature with the highest evaluation index in the initial candidate feature subset is locked as the core feature. In the feature adjacency knowledge graph, query all adjacent features that are connected to the core feature; All adjacency features are sorted in descending order of connection weight, and the top 3 adjacency features are extracted as strongly associated adjacency features. Strongly correlated adjacency features are added to the initial candidate feature subset, and deduplication is performed to obtain an expanded candidate feature subset.

[0050] In the local feature optimization process of decision tree nodes, the reduction in normalized weighted variance is used as the split evaluation index to avoid bias in the split evaluation due to differences in numerical scale of target state variables with different dimensions. Assume the system sets a preset lower limit for split quality, i.e., the threshold for the reduction in normalized weighted variance is 0.05. When the algorithm calculates the optimal split point in the initial candidate feature subset containing 7 features, if it finds that the maximum reduction in normalized weighted variance is 0.032, which is lower than the 0.05 threshold, it indicates that the features in the current candidate feature subset do not have regression discriminative ability. In this case, the feature with the highest reduction in normalized weighted variance (0.032) among these 7 features is locked as the core feature for this expansion, and expansion is carried out using the topological relationships of the features in the knowledge graph.

[0051] Using the core feature as an index, all first-degree neighbor nodes connected to the core feature are retrieved in the feature adjacency knowledge graph. Assume five adjacency features are retrieved, with corresponding edge connection weights of 1.5, 1.2, 0.9, 0.7, and 0.4. These adjacency features are sorted in descending order of their weights, and the top three features (those with weights of 1.5, 1.2, and 0.9) are selected and marked as strongly correlated adjacency features. These three features are merged into an initial candidate feature subset. If a feature already exists in the initial candidate feature subset, it is removed using hash deduplication logic. This operation generates an expanded candidate feature subset, which expands the feature space locally using graph knowledge, incorporating incremental features.

[0052] In an optional embodiment, updating the connection weights of the feature adjacency knowledge graph based on changes in evaluation metrics includes: Record the first maximum split evaluation index obtained based on the initial candidate feature subset, and the second maximum split evaluation index obtained by re-finding the optimal split point based on the expanded candidate feature subset; Calculate the difference between the second maximum splitting assessment index and the first maximum splitting assessment index; If the difference is greater than the preset increase threshold, the connection weight between the core feature that triggered this expansion and the included strongly associated neighboring features will be increased by 1 in the feature adjacency knowledge graph. If the difference is less than or equal to the preset improvement threshold, then in the feature adjacency knowledge graph, the connection weight between the core feature that triggered this expansion and the included strongly associated adjacency feature will be reduced by 0.1, and the reduced connection weight will be taken as the maximum value of 0, so that the connection weight remains non-negative.

[0053] Update the edge weights of the feature adjacency knowledge graph. Before feature expansion, record the first maximum split evaluation metric, i.e., the reduction in normalized weighted variance, as 0.032. After inputting strongly correlated adjacency features and re-finding the optimal split point based on the expanded candidate feature subset, assuming a newly input feature provides a reduction in normalized weighted variance of 0.094 at the new split point, record the second maximum split evaluation metric as 0.094. Perform a subtraction operation; the calculated difference is 0.062.

[0054] Since the calculated difference of 0.062 is greater than the system's preset boost threshold of 0.05, it indicates that the included strongly correlated adjacency features improved the regression split contribution. In this embodiment, the boost threshold is set to 0.05, consistent with the lower limit of split quality. When the training sample size N>5000, the boost threshold is lowered to 0.03 to retain weak gain feature associations; when N<1000, the boost threshold is raised to 0.06 to filter noisy associations. The boost threshold is determined through cross-validation, with a search range of [0.01, 0.10] and a step size of 0.01. Therefore, in the adjacency matrix of the graph, the weights of the undirected edges between the core feature and the three included strongly correlated adjacency features are increased by a step size of 1; for example, the original weight of 1.5 will be updated to 2.5. Conversely, if the second maximum split evaluation index is 0.038, and the calculated difference of 0.006≤0.05, it indicates that the included strongly correlated adjacency features did not make a substantial contribution. The connection weights between the core feature and these three features are reduced by 0.1, for example, from 1.5 to 1.4. Through iteration, the graph can filter out low-contribution feature associations and retain feature interaction rules that contribute to the splitting process.

[0055] S3, Input Data to Solve for Combined Instructions. Specifically, real-time collected data is input into the irrigation demand prediction model to obtain the expected value vector and covariance matrix at the current moment. A system response mapping model is established, mapping the combination instructions of irrigation duration, fertilization timing, and fertilizer ratio to the water state and nutrient concentration state. Using the expected value vector as the optimization center, the covariance matrix is ​​used to quantify uncertainty. Combined with the system response mapping model, a risk avoidance optimization problem is constructed, and the optimal combination instructions of irrigation duration, fertilization timing, and fertilizer ratio are solved under preset constraints.

[0056] The system receives the real-time data stream reported by the sensor at the current moment, calls the same feature engineering pipeline as the previous steps to convert it into a real-time composite time-series feature vector, performs prediction inference, parses the leaf node output to obtain the expected value vector and covariance matrix of the target water state and nutrient concentration state, and constructs a multilayer perceptron as the system response mapping model based on historical operation commands and corresponding state change data. This model contains three fully connected layers and uses the ReLU activation function for forward propagation to fit the nonlinear response relationship from control commands to water state and nutrient concentration state.

[0057] Based on the theory of integrated water and fertilizer management and the principles of fluid dynamics, the control command variable is defined as follows: ,in, Indicates irrigation duration, , , These represent the normalized time-series parameters representing the start time of nitrogen, phosphorus, and potassium fertilizer injection relative to the irrigation duration T, respectively. , , These represent the ratio adjustment coefficients or equivalent opening proportions of the injection channels for nitrogen, phosphorus, and potassium fertilizers, respectively. The actual injection start times are respectively... , and The constraints include 0 ≤ ≤ , 0≤ , , ≤1, 0≤ , , ≤1, and the ratio adjustment coefficient of each fertilizer injection channel satisfies ,in The water pump's rated flow rate, soil infiltration capacity, and single irrigation limit are preset. In this embodiment, the mother liquor concentration and basic flow rate of the fertilizer injection pump in each fertilizer channel are preset fixed parameters. The total fertilizer injection volume is determined by the irrigation duration, fertilizer injection sequence parameters, and the proportioning adjustment coefficient of each fertilizer injection channel. , , These are used to adjust the pumping duty cycle or equivalent opening of the corresponding fertilizer injection channel. The actual injection volume of the i-th fertilizer is calculated based on the corresponding mother liquor concentration, the base flow rate of the fertilizer injection pump, the proportioning adjustment coefficient of the corresponding fertilizer injection channel, and the injection duration, and is expressed as follows: ,in Let be the concentration of the mother liquor of the i-th fertilizer. The base flow rate corresponding to the fertilizer injection channel is i, which can be N, P, or K.

[0058] When the fertilizer injection pump supports independent flow rate adjustment, the injection flow rate of each fertilizer channel can be incorporated into the combined command as an extended control variable. Without setting a separate injection duration variable, each fertilizer channel continuously injects fertilizer from the corresponding injection start time until the end time of the irrigation. When independent control of the injection duration is required, the injection duration of each fertilizer channel can be incorporated into the combined command as an extended control variable.

[0059] An optimization objective function is constructed based on the concept of risk aversion. This objective function consists of two parts: the first part is a weighted evaluation term representing the difference between the water and nutrient concentration states predicted by the system response mapping model and the expected value vector; the second part is a risk penalty evaluation term based on the standard deviation of the main diagonal of the covariance matrix and the deviation from the conservative constraint boundary. Specifically, the state vector output by the system response mapping model is denoted as g(u), the expected value vector output by the irrigation demand prediction model is denoted as μ, and the regularized covariance matrix is ​​denoted as... The g(u)-μ is dimensionless according to the preset scaling factor or corresponding safety boundary width of each state variable, resulting in the normalized state deviation vector e(u). Based on optimization theory and risk aversion principles, the optimization objective can be adopted... The sum of the risk penalty term, where W is the preset state weight matrix; and the covariance matrix. This is used to calculate the risk safety margin of each state variable and shrink the safety constraint boundary, instead of directly forming the Jacobian matrix of the system response mapping model. The product of the form. The Jacobian matrix of the system response mapping model is only used as an intermediate quantity when the optimizer performs gradient iteration. Its dimension is the number of state variables multiplied by the number of control variables, avoiding matrix multiplication with the state covariance matrix that is dimension mismatched.

[0060] Set irrigation duration to be greater than or equal to 0 and not exceed the maximum irrigation duration. The sum of the ratio adjustment coefficients of each fertilizer injection channel is equal to 100%, the fertilizer injection timing parameters are in a closed interval of 0 to 1, and the water pump operating flow rate, fertilizer pump opening degree, fertilizer pump duty cycle, or equipment operating power calculated from each control parameter are set to not exceed the rated hardware threshold of the corresponding equipment as preset constraints. The minimum value of the objective function is iteratively searched, and the array of independent variables returned after the optimization solver converges is extracted. This array of independent variables is then reverse-analyzed into the optimal irrigation duration and the combined instruction of fertilizer injection timing and fertilizer injection channel ratio adjustment, and sent to the programmable logic controller for equipment control.

[0061] In an optional embodiment, using the expected value vector as the optimization center, the uncertainty is quantified using the covariance matrix, and a risk aversion optimization problem is constructed by combining the system response mapping model, including: From the expected value vector output by the irrigation demand prediction model, the expected values ​​of the target soil moisture content, target nitrogen concentration, target phosphorus concentration and target potassium concentration are extracted as the optimization center. From the covariance matrix output by the irrigation demand prediction model, extract the four variance elements on the main diagonal, and perform square root operation on each variance element to obtain the standard deviation of the corresponding target state variable. The risk safety margin is calculated based on the standard deviation of each target state variable. The risk safety margin is the product of the corresponding standard deviation and the preset confidence coefficient. Based on the risk safety margin, the preset safety constraint boundary of each target state variable is shrunk to obtain the conservative constraint boundary. The lower boundary of the shrunk boundary is non-negatively truncated by taking the maximum value of the result with 0. Using conservative constraint boundaries as the boundaries of constraints in risk-avoidance optimization problems tightens the feasible constraint range of the target state variables as the prediction uncertainty increases.

[0062] Irrigation demand prediction models, such as those using multivariate Gaussian regression tree models, output a probability distribution during forward inference, containing a 4-dimensional expectation vector and a 4×4 covariance matrix. From the expectation vector, four key parameter scalars for the future target time are extracted sequentially, for example: expected soil moisture content of 25%, expected nitrogen concentration of 120 mg / L, expected phosphorus concentration of 30 mg / L, and expected potassium concentration of 150 mg / L. These four scalar values ​​are then assembled into the target points, or optimization centers, of subsequent optimization control algorithms, such as the sequential quadratic programming (SQP) optimizer.

[0063] To map the model's uncertainty to constraint boundaries, a 4×4 covariance matrix is ​​analyzed, and four variance elements located on the main diagonal are extracted, for example, 4, 25, 4, and 16 respectively. The square roots of these four variances are then performed sequentially to obtain the predicted standard deviations of each variable: 2%, 5 mg / L, 2 mg / L, and 4 mg / L respectively. Based on the statistical properties of the Gaussian distribution, the risk margin for each target state variable is calculated. The risk margin is the product of the corresponding standard deviation and a pre-set confidence coefficient of 2.58. The confidence coefficient is determined based on the expected risk tolerance, with 2.58 corresponding to a 99% confidence level. This risk margin is not used to expand outwards from the expected value to form a feasible interval, but rather to contract the pre-set safety constraint boundaries inwards. Specifically, according to the principle of probability boundary contraction, if the pre-set safety constraint boundary of a certain state variable is... The corresponding risk safety margin is Then the conservative constraint boundary is .when As the value increases, the lower boundary rises and the upper boundary falls, tightening the feasible constraint range as uncertainty increases. If... Greater than In such cases, a conservative protection strategy is triggered; when the system has a preset minimum allowable width... When this happens, the constraint boundary of the variable is set to the safety center value. Centered on, with a width of The minimum allowable range; when no minimum allowable width is preset or the range still exceeds the equipment safety boundary, output a manual review prompt to the controller and suspend automatic optimization distribution.

[0064] Taking nitrogen concentration as an example, if the preset safety constraint boundary range is 100 mg / L to 140 mg / L, the prediction standard deviation is 5 mg / L, and the confidence coefficient is 2.58, then the risk safety margin is 12.9 mg / L, and the conservative constraint boundary range is 112.9 mg / L to 127.1 mg / L. The calculated lower bound value is then compared to 0, taking the maximum value. If the calculated lower bound for phosphorus at a certain low concentration is -1.5 mg / L, it is corrected to 0 mg / L after non-negative truncation. These four sets of conservative constraint boundaries are set as the nonlinear inequality constraint boundaries of the optimizer, ensuring crop water and fertilizer supply while avoiding the risk of over-irrigation or fertilizer damage caused by model prediction errors.

[0065] The ablation experiment used historical observation data of greenhouse tomatoes from a large agricultural demonstration park over the past two years, divided into training, validation, and test sets in a 7:2:1 ratio. Input features included environmental variables and historical water and fertilizer records. The prediction targets were soil moisture and nitrogen, phosphorus, and potassium concentrations. Evaluation indicators included normalized mean absolute error and risk incidence rate. The baseline group used a traditional random forest model and single-timescale original features. Experimental group one introduced a multi-scale composite time-series feature extraction module based on the baseline group. Experimental group two added a probability decay sampling and candidate feature expansion mechanism guided by a feature adjacency knowledge graph, building upon experimental group one. The complete group further introduced risk avoidance constraint control based on the covariance matrix safety margin shrinkage boundary, building upon experimental group two.

[0066] Test results show that the normalized mean absolute error of the baseline group was 5.2%, and the risk incidence rate was 12.5%; the normalized mean absolute error of experimental group 1 decreased to 4.3%, and the risk incidence rate decreased to 10.6%; the normalized mean absolute error of experimental group 2 further decreased to 3.1%, and the risk incidence rate decreased to 7.8%; the normalized mean absolute error of the complete group was 3%, and the risk incidence rate decreased to 2.1%.

[0067] Experimental results show that multi-scale composite time-series features can characterize the cumulative effects and trends under different periods, thereby improving the basic prediction accuracy; the sampling and expansion mechanism guided by the feature adjacency knowledge graph can reduce feature redundancy in the split path, improve the diversity of the tree model, reduce the risk of the model getting trapped in local optima, and further reduce the prediction error; the safety margin contraction constraint based on the covariance matrix can transform the model prediction uncertainty into a conservative optimization boundary, reducing the risk of over-irrigation and over-fertilization while maintaining prediction accuracy.

[0068] This invention also discloses an integrated water and fertilizer irrigation scheduling optimization system, including a processor and a memory. The memory stores computer program instructions, which, when executed by the processor, implement an integrated water and fertilizer irrigation scheduling optimization method according to the present invention.

[0069] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for optimizing the scheduling of integrated water and fertilizer irrigation, characterized in that, include: S1. Obtain historical environmental data, crop growth data, and irrigation data of the integrated water and fertilizer system, extract the cumulative effect value and rate of change value under multi-scale time windows, and generate composite time series features; S2. Construct an irrigation demand prediction model. During training, initialize the feature adjacency knowledge graph representing the high-order interaction relationships of features. Track existing splitting paths at the decision tree splitting nodes, quantify the coupling degree between candidate features and path features, and reduce the sampling probability accordingly to generate an initial candidate feature subset. When the maximum splitting evaluation index based on the initial candidate feature subset is lower than the preset lower limit, incorporate the strongly correlated adjacency features of the best feature in the initial candidate feature subset according to the feature adjacency knowledge graph to obtain an expanded candidate feature subset and find the optimal splitting point again. Update the connection weights of the feature adjacency knowledge graph according to the change of evaluation index until the model is generated. The leaf nodes of the irrigation demand prediction model output the expected value vector and covariance matrix. S3. Input the real-time collected data into the irrigation demand prediction model to obtain the expected value vector and covariance matrix at the current moment, establish a system response mapping model, take the expected value vector as the optimization center, use the covariance matrix to quantify uncertainty, combine the system response mapping model to construct a risk avoidance optimization problem, and solve the optimal irrigation duration, fertilization sequence and fertilizer ratio combination instructions.

2. The water and fertilizer integrated irrigation scheduling optimization method according to claim 1, characterized in that, The method for obtaining the composite time series features includes: setting time windows of lengths of 1 day, 7 days, and 30 days, and truncating the continuous historical time series observation sequence into three benchmark observation stage windows with different periods; for each benchmark observation stage window, for flux-type variables, summing the historical time series observation data within the window, and for state-type variables, calculating the average value of the historical time series observation data within the window, and unifying the summation result with the average value as the cumulative effect value of the corresponding time scale; subtracting the same type of time series observation data from the time series observation data of the 1st day before the current time from the same type of time series observation data of the 2nd day before the current time to obtain the difference, dividing the difference by the corrected benchmark value to obtain the corresponding data change rate value, where the corrected benchmark value is the sum of the same type of time series observation data of the 2nd day before the current time and a preset small positive number of the same dimension; and concatenating all extracted cumulative effect values ​​and change rate values ​​to generate composite time series features.

3. The water and fertilizer integrated irrigation scheduling optimization method according to claim 2, characterized in that, The initialization of the feature adjacency knowledge graph representing high-order interaction relationships during training includes: The generated composite temporal features are used as node entities to construct an initial feature adjacency knowledge graph. The Pearson correlation coefficient between any two feature nodes is calculated using the formula: the covariance of the two feature sequences divided by the product of the standard deviations of the two feature sequences. Determine whether the absolute value of the Pearson correlation coefficient between any two feature nodes is greater than a preset correlation threshold. If so, establish an undirected edge between the two feature nodes and use the undirected edge as the initial connection relationship of the feature adjacency knowledge graph.

4. The method for optimizing the scheduling of integrated water and fertilizer irrigation according to claim 1, characterized in that, The process of quantifying the coupling degree between candidate features and path features and thereby attenuating the sampling probability to generate an initial candidate feature subset includes: backtracking upwards from the current split node to the root node to obtain all used split features on the split path; querying the non-negative connection weights between each candidate feature and used split features in the candidate feature set in the feature adjacency knowledge graph; and summing the non-negative connection weights between a single candidate feature and all used split features to obtain the coupling degree of the candidate feature. The coupling degree of each candidate feature is normalized to a closed interval of 0 to 1. The probability decay coefficient is obtained by subtracting the normalized coupling degree from 1. The initial sampling probability of each candidate feature is multiplied by the probability decay coefficient to obtain the adjusted sampling probability. The adjusted sampling probabilities are then summed and normalized. Random sampling is performed based on the normalized adjusted sampling probabilities to generate an initial subset of candidate features.

5. The method for optimizing the scheduling of integrated water and fertilizer irrigation according to claim 1, characterized in that, The step of incorporating the strongly correlated adjacent features of the optimal feature into the feature adjacency knowledge graph to obtain an expanded candidate feature subset includes: when the maximum splitting evaluation index of the initial candidate feature subset is lower than the preset splitting quality lower limit, locking the feature with the highest evaluation index in the initial candidate feature subset as the core feature. In the feature adjacency knowledge graph, query all adjacent features that are connected to the core feature; sort all adjacent features in descending order of connection weight, and extract the top 3 adjacent features as strongly associated adjacent features; add the strongly associated adjacent features to the initial candidate feature subset and perform deduplication to obtain the expanded candidate feature subset.

6. The water and fertilizer integrated irrigation scheduling optimization method according to claim 5, characterized in that, The step of updating the connection weights of the feature adjacency knowledge graph based on changes in evaluation metrics includes: Record the first maximum split evaluation index obtained based on the initial candidate feature subset, and the second maximum split evaluation index obtained by re-finding the optimal split point based on the expanded candidate feature subset; Calculate the difference between the second maximum split evaluation index and the first maximum split evaluation index; if the difference is greater than the preset improvement threshold, then in the feature adjacency knowledge graph, increase the connection weight between the core feature and the included strongly associated adjacency features by 1; if the difference is less than or equal to the preset improvement threshold, then in the feature adjacency knowledge graph, decrease the connection weight between the core feature and the included strongly associated adjacency features by 0.1, and take the maximum value between the reduced connection weight and 0 to keep the connection weight non-negative.

7. The method for optimizing the scheduling of integrated water and fertilizer irrigation according to claim 1, characterized in that, The leaf node output expected value vector and covariance matrix of the irrigation demand prediction model include: Based on the training sample set of the falling leaf nodes, the mean values ​​of the target water and nutrient concentrations are calculated as the expected value vector, and the corresponding sample covariance matrix is ​​calculated. When the number of samples actually falling into a leaf node is less than the dimension of the target state variable plus 1, the sample covariance matrix is ​​calculated by backtracking to the parent node sample set, or a diagonal covariance matrix is ​​constructed from the diagonal elements of the global training sample covariance matrix as a degenerate alternative; a regularization term is superimposed on the main diagonal of the obtained sample covariance matrix or diagonal covariance matrix to generate the final covariance matrix as the output result of the leaf node.

8. The method for optimizing the scheduling of integrated water and fertilizer irrigation according to claim 1, characterized in that, The risk avoidance optimization problem constructed by combining the system response mapping model includes: extracting the expected values ​​of target soil moisture content, target nitrogen concentration, target phosphorus concentration, and target potassium concentration from the expected value vector output by the irrigation demand prediction model, and using them as optimization centers; extracting the variance elements on the main diagonal from the covariance matrix, and performing square root operations on the variance elements to obtain the standard deviation of the corresponding target state variables; calculating the risk safety margin based on the standard deviation of each target state variable, where the risk safety margin is the product of the corresponding standard deviation and a pre-set confidence coefficient; Based on the risk safety margin, the preset safety constraint boundary of each objective state variable is shrunk to obtain the conservative constraint boundary. The lower boundary after shrinkage is truncated by taking the maximum value of 0. The conservative constraint boundary is used as the boundary of the constraint conditions in the risk avoidance optimization problem.

9. The method for optimizing the scheduling of integrated water and fertilizer irrigation according to claim 8, characterized in that, The instructions for finding the optimal combination of irrigation duration, fertilization timing, and fertilizer ratio include: An optimization objective function is constructed based on the idea of ​​risk aversion. The optimization objective function includes a first part and a second part. The first part is a weighted evaluation term of the difference between the water state and nutrient concentration state predicted by the system response mapping model and the expected value vector. The second part is a risk penalty evaluation term based on the deviation of the conservative constraint boundary and the standard deviation of the corresponding target state variable. The algorithm iteratively searches for the minimum value of the objective function and extracts the array of independent variables returned by the optimizer after convergence. The array of independent variables is then reverse-analyzed into the optimal combination of irrigation duration, fertilization timing, and fertilizer ratio instructions.

10. A water and fertilizer integrated irrigation scheduling and optimization system, characterized in that, include: A processor and a memory, wherein the memory stores computer program instructions that, when executed by the processor, implement the water and fertilizer integrated irrigation scheduling optimization method according to any one of claims 1-9.

Images